ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'develop' into pair-lj-lepton-sphere

Browse Source
This commit is contained in:
Axel Kohlmeyer
2023-04-03 17:46:34 -04:00
parent 3f40d1ef5d c9b6fad6a0
commit dfd16c95ff
85 changed files with 594 additions and 3352 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
cmake/CMakeLists.txt
View File

@ -257,7 +257,6 @@ set(STANDARD_PACKAGES
KIM
KSPACE
LATBOLTZ
LATTE
LEPTON
MACHDYN
MANIFOLD
@ -441,7 +440,7 @@ if(BUILD_OMP)
target_link_libraries(lmp PRIVATE OpenMP::OpenMP_CXX)
endif()
if(PKG_MSCG OR PKG_ATC OR PKG_AWPMD OR PKG_ML-QUIP OR PKG_ML-POD OR PKG_LATTE OR PKG_ELECTRODE)
if(PKG_MSCG OR PKG_ATC OR PKG_AWPMD OR PKG_ML-QUIP OR PKG_ML-POD OR PKG_ELECTRODE)
enable_language(C)
if (NOT USE_INTERNAL_LINALG)
find_package(LAPACK)
@ -521,7 +520,7 @@ else()
endif()
foreach(PKG_WITH_INCL KSPACE PYTHON ML-IAP VORONOI COLVARS ML-HDNNP MDI MOLFILE NETCDF
PLUMED QMMM ML-QUIP SCAFACOS MACHDYN VTK KIM LATTE MSCG COMPRESS ML-PACE LEPTON)
PLUMED QMMM ML-QUIP SCAFACOS MACHDYN VTK KIM MSCG COMPRESS ML-PACE LEPTON)
if(PKG_${PKG_WITH_INCL})
include(Packages/${PKG_WITH_INCL})
endif()

54
cmake/Modules/Packages/LATTE.cmake
View File

@ -1,54 +0,0 @@
enable_language(Fortran)
# using lammps in a super-build setting
if(TARGET LATTE::latte)
target_link_libraries(lammps PRIVATE LATTE::latte)
return()
endif()
find_package(LATTE 1.2.2 CONFIG)
if(LATTE_FOUND)
set(DOWNLOAD_LATTE_DEFAULT OFF)
else()
set(DOWNLOAD_LATTE_DEFAULT ON)
endif()
option(DOWNLOAD_LATTE "Download the LATTE library instead of using an already installed one" ${DOWNLOAD_LATTE_DEFAULT})
if(DOWNLOAD_LATTE)
message(STATUS "LATTE download requested - we will build our own")
set(LATTE_URL "https://github.com/lanl/LATTE/archive/v1.2.2.tar.gz" CACHE STRING "URL for LATTE tarball")
set(LATTE_MD5 "820e73a457ced178c08c71389a385de7" CACHE STRING "MD5 checksum of LATTE tarball")
mark_as_advanced(LATTE_URL)
mark_as_advanced(LATTE_MD5)
GetFallbackURL(LATTE_URL LATTE_FALLBACK)
# CMake cannot pass BLAS or LAPACK library variable to external project if they are a list
list(LENGTH BLAS_LIBRARIES} NUM_BLAS)
list(LENGTH LAPACK_LIBRARIES NUM_LAPACK)
if((NUM_BLAS GREATER 1) OR (NUM_LAPACK GREATER 1) AND NOT USE_INTERNAL_LINALG)
message(FATAL_ERROR "Cannot compile downloaded LATTE library due to a technical limitation. "
"Try to configure LAMMPS with '-D USE_INTERNAL_LINALG=on' added as a workaround.")
endif()
include(ExternalProject)
ExternalProject_Add(latte_build
URL ${LATTE_URL} ${LATTE_FALLBACK}
URL_MD5 ${LATTE_MD5}
SOURCE_SUBDIR cmake
CMAKE_ARGS -DCMAKE_INSTALL_PREFIX=<INSTALL_DIR> ${CMAKE_REQUEST_PIC} -DCMAKE_INSTALL_LIBDIR=lib
-DBLAS_LIBRARIES=${BLAS_LIBRARIES} -DLAPACK_LIBRARIES=${LAPACK_LIBRARIES}
-DCMAKE_Fortran_COMPILER=${CMAKE_Fortran_COMPILER} -DCMAKE_Fortran_FLAGS=${CMAKE_Fortran_FLAGS}
-DCMAKE_Fortran_FLAGS_${BTYPE}=${CMAKE_Fortran_FLAGS_${BTYPE}} -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}
-DCMAKE_MAKE_PROGRAM=${CMAKE_MAKE_PROGRAM} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TOOLCHAIN_FILE}
BUILD_BYPRODUCTS <INSTALL_DIR>/lib/liblatte.a
)
ExternalProject_get_property(latte_build INSTALL_DIR)
add_library(LAMMPS::LATTE UNKNOWN IMPORTED)
set_target_properties(LAMMPS::LATTE PROPERTIES
IMPORTED_LOCATION "${INSTALL_DIR}/lib/liblatte.a"
INTERFACE_LINK_LIBRARIES "${LAPACK_LIBRARIES}")
target_link_libraries(lammps PRIVATE LAMMPS::LATTE)
add_dependencies(LAMMPS::LATTE latte_build)
else()
find_package(LATTE 1.2.2 REQUIRED CONFIG)
target_link_libraries(lammps PRIVATE LATTE::latte)
endif()

1
cmake/presets/all_off.cmake
View File

@ -43,7 +43,6 @@ set(ALL_PACKAGES
KOKKOS
KSPACE
LATBOLTZ
LATTE
LEPTON
MACHDYN
MANIFOLD

1
cmake/presets/all_on.cmake
View File

@ -45,7 +45,6 @@ set(ALL_PACKAGES
KOKKOS
KSPACE
LATBOLTZ
LATTE
LEPTON
MACHDYN
MANIFOLD

3
cmake/presets/download.cmake
View File

@ -1,14 +1,13 @@
# Preset that turns on packages with automatic downloads of sources or potentials.
# Compilation of libraries like Plumed or ScaFaCoS can take a considerable amount of time.
set(ALL_PACKAGES KIM LATTE MSCG VORONOI PLUMED SCAFACOS MACHDYN MESONT MDI ML-PACE)
set(ALL_PACKAGES KIM MSCG VORONOI PLUMED SCAFACOS MACHDYN MESONT MDI ML-PACE)
foreach(PKG ${ALL_PACKAGES})
set(PKG_${PKG} ON CACHE BOOL "" FORCE)
endforeach()
set(DOWNLOAD_KIM ON CACHE BOOL "" FORCE)
set(DOWNLOAD_LATTE ON CACHE BOOL "" FORCE)
set(DOWNLOAD_MDI ON CACHE BOOL "" FORCE)
set(DOWNLOAD_MSCG ON CACHE BOOL "" FORCE)
set(DOWNLOAD_VORO ON CACHE BOOL "" FORCE)

1
cmake/presets/mingw-cross.cmake
View File

@ -35,7 +35,6 @@ set(WIN_PACKAGES
INTEL
INTERLAYER
KSPACE
LATTE
LEPTON
MACHDYN
MANIFOLD

1
cmake/presets/nolib.cmake
View File

@ -12,7 +12,6 @@ set(PACKAGES_WITH_LIB
KIM
KOKKOS
LATBOLTZ
LATTE
LEPTON
MACHDYN
MDI

16
doc/src/Build_basics.rst
View File

@ -128,14 +128,14 @@ and adds vectorization support when compiled with compatible compilers,
in particular the Intel compilers on top of OpenMP. Also, the ``KOKKOS``
package can be compiled to include OpenMP threading.
In addition, there are a few commands in LAMMPS that have native OpenMP
support included as well. These are commands in the ``MPIIO``,
``ML-SNAP``, ``DIFFRACTION``, and ``DPD-REACT`` packages. Furthermore,
some packages support OpenMP threading indirectly through the libraries
they interface to: e.g. ``LATTE``, ``KSPACE``, and ``COLVARS``. See the
:doc:`Packages details <Packages_details>` page for more info on these
packages, and the pages for their respective commands for OpenMP
threading info.
In addition, there are a few commands in LAMMPS that have native
OpenMP support included as well. These are commands in the ``MPIIO``,
``ML-SNAP``, ``DIFFRACTION``, and ``DPD-REACT`` packages.
Furthermore, some packages support OpenMP threading indirectly through
the libraries they interface to: e.g. ``KSPACE``, and ``COLVARS``.
See the :doc:`Packages details <Packages_details>` page for more info
on these packages, and the pages for their respective commands for
OpenMP threading info.
For CMake, if you use ``BUILD_OMP=yes``, you can use these packages
and turn on their native OpenMP support and turn on their native OpenMP

64
doc/src/Build_extras.rst
View File

@ -43,7 +43,6 @@ This is the list of packages that may require additional steps.
* :ref:`INTEL <intel>`
* :ref:`KIM <kim>`
* :ref:`KOKKOS <kokkos>`
* :ref:`LATTE <latte>`
* :ref:`LEPTON <lepton>`
* :ref:`MACHDYN <machdyn>`
* :ref:`MDI <mdi>`
@ -832,63 +831,6 @@ will thus always enable it.
----------
.. _latte:
LATTE package
-------------------------
To build with this package, you must download and build the LATTE
library.
.. tabs::
.. tab:: CMake build
.. code-block:: bash
-D DOWNLOAD_LATTE=value # download LATTE for build, value = no (default) or yes
-D LATTE_LIBRARY=path # LATTE library file (only needed if a custom location)
-D USE_INTERNAL_LINALG=value # Use the internal linear algebra library instead of LAPACK
# value = no (default) or yes
If ``DOWNLOAD_LATTE`` is set, the LATTE library will be downloaded
and built inside the CMake build directory. If the LATTE library
is already on your system (in a location CMake cannot find it),
``LATTE_LIBRARY`` is the filename (plus path) of the LATTE library
file, not the directory the library file is in.
The LATTE library requires LAPACK (and BLAS) and CMake can identify
their locations and pass that info to the LATTE build script. But
on some systems this triggers a (current) limitation of CMake and
the configuration will fail. Try enabling ``USE_INTERNAL_LINALG`` in
those cases to use the bundled linear algebra library and work around
the limitation.
.. tab:: Traditional make
You can download and build the LATTE library manually if you
prefer; follow the instructions in ``lib/latte/README``\ . You
can also do it in one step from the ``lammps/src`` dir, using a
command like these, which simply invokes the
``lib/latte/Install.py`` script with the specified args:
.. code-block:: bash
make lib-latte # print help message
make lib-latte args="-b" # download and build in lib/latte/LATTE-master
make lib-latte args="-p $HOME/latte" # use existing LATTE installation in $HOME/latte
make lib-latte args="-b -m gfortran" # download and build in lib/latte and
# copy Makefile.lammps.gfortran to Makefile.lammps
Note that 3 symbolic (soft) links, ``includelink`` and ``liblink``
and ``filelink.o``, are created in ``lib/latte`` to point to
required folders and files in the LATTE home directory. When
LAMMPS itself is built it will use these links. You should also
check that the ``Makefile.lammps`` file you create is appropriate
for the compiler you use on your system to build LATTE.
----------
.. _lepton:
LEPTON package
@ -1413,9 +1355,9 @@ This package depends on the KSPACE package.
KSPACE package so the latter one *must* be enabled.
The ELECTRODE package also requires LAPACK (and BLAS) and CMake
can identify their locations and pass that info to the LATTE build
script. But on some systems this may cause problems when linking
or the dependency is not desired. Try enabling
can identify their locations and pass that info to the ELECTRODE
build script. But on some systems this may cause problems when
linking or the dependency is not desired. Try enabling
``USE_INTERNAL_LINALG`` in those cases to use the bundled linear
algebra library and work around the limitation.

1
doc/src/Build_package.rst
View File

@ -46,7 +46,6 @@ packages:
* :ref:`INTEL <intel>`
* :ref:`KIM <kim>`
* :ref:`KOKKOS <kokkos>`
* :ref:`LATTE <latte>`
* :ref:`LEPTON <lepton>`
* :ref:`MACHDYN <machdyn>`
* :ref:`MDI <mdi>`

1
doc/src/Commands_fix.rst
View File

@ -104,7 +104,6 @@ OPT.
* :doc:`langevin/drude <fix_langevin_drude>`
* :doc:`langevin/eff <fix_langevin_eff>`
* :doc:`langevin/spin <fix_langevin_spin>`
* :doc:`latte <fix_latte>`
* :doc:`lb/fluid <fix_lb_fluid>`
* :doc:`lb/momentum <fix_lb_momentum>`
* :doc:`lb/viscous <fix_lb_viscous>`

14
doc/src/Commands_removed.rst
View File

@ -38,6 +38,20 @@ been folded into the :doc:`reset_atoms <reset_atoms>` command. If
present, LAMMPS will replace the commands accordingly and print a
warning.
LATTE package
-------------
.. deprecated:: TBD
The LATTE package with the fix latte command was removed from LAMMPS.
This functionality has been superseded by :doc:`fix mdi/qm <fix_mdi_qm>`
and :doc:`fix mdi/qmmm <fix_mdi_qmmm>` from the :ref:`MDI package
<PKG-MDI>`. These fixes are compatible with several quantum software
packages, including LATTE. See the ``examples/QUANTUM`` dir and the
:doc:`MDI coupling HOWTO <Howto_mdi>` page. MDI supports running LAMMPS
with LATTE as a plugin library (similar to the way fix latte worked), as
well as on a different set of MPI processors.
MEAM package
------------

2
doc/src/Examples.rst
View File

@ -94,8 +94,6 @@ Lowercase directories
+-------------+------------------------------------------------------------------+
| kim | use of potentials from the `OpenKIM Repository <openkim_>`_ |
+-------------+------------------------------------------------------------------+
| latte | examples for using fix latte for DFTB via the LATTE library |
+-------------+------------------------------------------------------------------+
| mdi | use of the MDI package and MolSSI MDI code coupling library |
+-------------+------------------------------------------------------------------+
| meam | MEAM test for SiC and shear (same as shear examples) |

18
doc/src/Howto_couple.rst
View File

@ -12,16 +12,16 @@ LAMMPS can be coupled to other codes in at least 4 different ways. Each
has advantages and disadvantages, which you will have to think about in
the context of your application.
1. Define a new :doc:`fix <fix>` command that calls the other code. In
this scenario, LAMMPS is the driver code. During timestepping, the
fix is invoked, and can make library calls to the other code, which
has been linked to LAMMPS as a library. This is the way the
:ref:`LATTE <PKG-LATTE>` package, which performs density-functional
tight-binding calculations using the `LATTE software
<https://github.com/lanl/LATTE>`_ to compute forces, is interfaced to
LAMMPS. See the :doc:`fix latte <fix_latte>` command for more
1. Define a new :doc:`fix <fix>` or :doc:`compute <compute>` command
that calls the other code. In this scenario, LAMMPS is the driver
code. During timestepping, the fix or compute is invoked, and can
make library calls to the other code, which has been linked to LAMMPS
as a library. This is the way the :ref:`VORONOI <PKG-VORONOI>`
package, which computes Voronoi tesselations using the `Voro++
library <http://math.lbl.gov/voro++>`_, is interfaced to LAMMPS. See
the :doc:`compute voronoi <compute_voronoi_atom>` command for more
details. Also see the :doc:`Modify <Modify>` pages for information
on how to add a new fix to LAMMPS.
on how to add a new fix or compute to LAMMPS.
.. spacer

2
doc/src/Howto_tip4p.rst
View File

@ -101,7 +101,7 @@ not as part of the pair coefficients.
- 0.52422
* - LJ :math:`\epsilon` of OO (kcal/mole)
- 0.1550
- 0.1577
- 0.21084
- 0.1852
- 0.16275
* - LJ :math:`\sigma` of OO (:math:`\AA`)

2
doc/src/Install_mac.rst
View File

@ -6,7 +6,7 @@ LAMMPS can be downloaded, built, and configured for OS X on a Mac with
instructions for :doc:`downloading an executable via Conda
<Install_conda>`.) The following LAMMPS packages are unavailable at
this time because of additional requirements not yet met: GPU, KOKKOS,
LATTE, MSCG, MPIIO, POEMS, VORONOI.
MSCG, MPIIO, POEMS, VORONOI.
After installing Homebrew, you can install LAMMPS on your system with
the following commands:

2
doc/src/Intro_features.rst
View File

@ -88,7 +88,7 @@ commands)
* charge equilibration (QEq via dynamic, point, shielded, Slater methods)
* coarse-grained potentials: DPD, GayBerne, REsquared, colloidal, DLVO, oxDNA / oxRNA, SPICA
* mesoscopic potentials: granular, Peridynamics, SPH, mesoscopic tubular potential (MESONT)
* semi-empirical potentials: multi-ion generalized pseudopotential theory (MGPT), second moment tight binding + QEq (SMTB-Q), density functional tight-binding (LATTE)
* semi-empirical potentials: multi-ion generalized pseudopotential theory (MGPT), second moment tight binding + QEq (SMTB-Q)
* electron force field (eFF, AWPMD)
* bond potentials: harmonic, FENE, Morse, nonlinear, Class II (COMPASS), quartic (breakable), tabulated, scripted
* angle potentials: harmonic, CHARMM, cosine, cosine/squared, cosine/periodic, Class II (COMPASS), tabulated, scripted

38
doc/src/Packages_details.rst
View File

@ -67,7 +67,6 @@ page gives those details.
* :ref:`KOKKOS <PKG-KOKKOS>`
* :ref:`KSPACE <PKG-KSPACE>`
* :ref:`LATBOLTZ <PKG-LATBOLTZ>`
* :ref:`LATTE <PKG-LATTE>`
* :ref:`LEPTON <PKG-LEPTON>`
* :ref:`MACHDYN <PKG-MACHDYN>`
* :ref:`MANIFOLD <PKG-MANIFOLD>`
@ -1357,43 +1356,6 @@ The LATBOLTZ package requires that LAMMPS is build in :ref:`MPI parallel mode <s
----------
.. _PKG-LATTE:
LATTE package
-------------
**Contents:**
A fix command which wraps the LATTE DFTB code, so that molecular
dynamics can be run with LAMMPS using density-functional tight-binding
quantum forces calculated by LATTE.
More information on LATTE can be found at this website:
`https://github.com/lanl/LATTE <latte-home_>`_. A brief technical
description is given with the :doc:`fix latte <fix_latte>` command.
.. _latte-home: https://github.com/lanl/LATTE
**Authors:** Christian Negre (LANL) and Steve Plimpton (Sandia). LATTE
itself is developed at Los Alamos National Laboratory by Marc
Cawkwell, Anders Niklasson, and Christian Negre.
**Install:**
This package has :ref:`specific installation instructions <latte>` on
the :doc:`Build extras <Build_extras>` page.
**Supporting info:**
* src/LATTE: filenames -> commands
* src/LATTE/README
* lib/latte/README
* :doc:`fix latte <fix_latte>`
* examples/latte
* `LAMMPS-LATTE tutorial <https://github.com/lanl/LATTE/wiki/Using-LATTE-through-LAMMPS>`_
----------
.. _PKG-LEPTON:
LEPTON package

5
doc/src/Packages_list.rst
View File

@ -233,11 +233,6 @@ whether an extra library is needed to build and use the package:
- :doc:`fix lb/fluid <fix_lb_fluid>`
- PACKAGES/latboltz
- no
* - :ref:`LATTE <PKG-LATTE>`
- quantum DFTB forces via LATTE
- :doc:`fix latte <fix_latte>`
- latte
- ext
* - :ref:`LEPTON <PKG-LEPTON>`
- evaluate strings as potential function
- :doc:`pair_style lepton <pair_lepton>`

1
doc/src/fix.rst
View File

@ -256,7 +256,6 @@ accelerated styles exist.
* :doc:`langevin/drude <fix_langevin_drude>` - Langevin temperature control of Drude oscillators
* :doc:`langevin/eff <fix_langevin_eff>` - Langevin temperature control for the electron force field model
* :doc:`langevin/spin <fix_langevin_spin>` - Langevin temperature control for a spin or spin-lattice system
* :doc:`latte <fix_latte>` - wrapper on LATTE density-functional tight-binding code
* :doc:`lb/fluid <fix_lb_fluid>` - lattice-Boltzmann fluid on a uniform mesh
* :doc:`lb/momentum <fix_lb_momentum>` - :doc:`fix momentum <fix_momentum>` replacement for use with a lattice-Boltzmann fluid
* :doc:`lb/viscous <fix_lb_viscous>` - :doc:`fix viscous <fix_viscous>` replacement for use with a lattice-Boltzmann fluid

2
doc/src/fix_indent.rst
View File

@ -19,7 +19,7 @@ Syntax
.. parsed-literal::
*sphere* args = x y z R
x,y,z = initial position of center of indenter (distance units)
x,y,z = position of center of indenter (distance units)
R = sphere radius of indenter (distance units)
any of x,y,z,R can be a variable (see below)
*cylinder* args = dim c1 c2 R

269
doc/src/fix_latte.rst
View File

@ -1,269 +0,0 @@
.. index:: fix latte
fix latte command
=================
Syntax
""""""
.. parsed-literal::
fix ID group-ID latte keyword value ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* latte = style name of this fix command
* zero or more keyword/value pairs may be appended
.. parsed-literal::
keyword = *coulomb* or *exclude*
*coulomb* value = peID
peID = ID of compute used to calculate per-atom energy
*exclude* value = fixID
fixID = ID of fix which potentially excludes atoms before calling LATTE
Examples
""""""""
.. code-block:: LAMMPS
fix dftb all latte
fix dftb all exclude GCMC
Description
"""""""""""
This fix style is a wrapper on the self-consistent charge transfer
density functional based tight binding (DFTB) code LATTE. If you
download and build LATTE, it can be called as a library by LAMMPS via
this fix to run dynamics or perform energy minimization using DFTB
forces and energies computed by LATTE.
LATTE is principally developed and supported by Marc Cawkwell and
co-workers at Los Alamos National Laboratory (LANL). See the full
list of contributors in the src/LATTE/README file.
To use this fix, the LATTE program needs to be compiled as a library
and linked with LAMMPS. LATTE can be downloaded (or cloned) from
`https://github.com/lanl/LATTE <https://github.com/lanl/LATTE>`_.
Instructions on how to download and build LATTE on your system can be
found in the lib/latte/README. Note that you can also use the "make
lib-latte" command from the LAMMPS src directory to automate this
process.
Once LAMMPS is built with the LATTE package, you can run the example
input scripts for molecular dynamics or energy minimization that are
found in examples/latte.
A step-by-step tutorial can be followed at: `LAMMPS-LATTE tutorial <https://github.com/lanl/LATTE/wiki/Using-LATTE-through-LAMMPS>`_
Currently, LAMMPS must be run in serial or as a single MPI task, to
use this fix. This is because the version of the LATTE library LAMMPS
uses does not support MPI. On the LAMMPS size, this is typically not
a bottleneck, since LATTE will be doing 99% or more of the work to
compute quantum-accurate forces. On the LATTE side, the LATTE library
does support threaded parallelism via OpenMP. You must build the
LATTE library with OpenMP support, then set the OMP_NUM_THREADS
environment variable before performing a LAMMPS + LATTE simulation to
tell LATTE how many threads to invoke.
.. note::
NEB calculations can be done using this fix using multiple
replicas and running LAMMPS in parallel. However, each replica must
be run on a single MPI task. For details, see the :doc:`neb <neb>`
command page and the :doc:`-partition command-line switch <Run_options>`
----------
The *coulomb* argument is not yet supported by fix latte (as of Sept
2022). Eventually it will be used to enable LAMMPS to calculate a
Coulomb potential as an alternative to LATTE performing the
calculation.
.. versionadded:: 15Sep2022
The *exclude* argument allows this fix to work in tandem with another
fix which may decide to delete one or more atoms of molecules. The
specified fixID is the ID of the other fix.
The one current example of such a fix is the :doc:`fix gcmc
<fix_gcmc>` command which performs Monte Carlo insertions and
deletions. If a trial deletion is performed, then LAMMPS needs to
only pass LATTE the atoms which remain. Fix gcmc does not actually
remove any atoms until after the new energy is computed (in this case
by LATTE), and a Monte Carlo accept/reject decision is made for the
trial deletion.
----------
LATTE is a code for performing self-consistent charge transfer
tight-binding (SC-TB) calculations of total energies and the forces
acting on atoms in molecules and solids. This tight-binding method is
becoming more and more popular and widely used in chemistry,
biochemistry, material science, etc.
The SC-TB formalism is derived from an expansion of the Kohn-Sham
density functional to second order in charge fluctuations about a
reference charge of overlapping atom-centered densities and bond
integrals are parameterized using a Slater-Koster tight-binding
approach. This procedure, which usually is referred to as the DFTB
method has been described in detail by (:ref:`Elstner <Elstner>`) and
(:ref:`Finnis <Finnis2>`) and coworkers.
The work of the LATTE developers follows that of Elstner closely with
respect to the physical model. However, the development of LATTE is
geared principally toward large-scale, long duration, microcanonical
quantum-based Born-Oppenheimer molecular dynamics (QMD) simulations.
One of the main bottlenecks of an electronic structure calculation is
the solution of the generalized eigenvalue problem which scales with
the cube of the system size O(N\^3).
The Theoretical and Computer sciences divisions at Los Alamos National
Laboratory have accumulated large experience addressing this issue by
calculating the density matrix directly instead of using
diagonalization. We typically use a recursive sparse Fermi-operator
expansion using second-order spectral projection functions
(SP2-algorithm), which was introduced by Niklasson in 2002
(:ref:`Niklasson2002 <Niklasson2002>`), (:ref:`Rubensson <Rubensson>`),
(:ref:`Mniszewski <Mniszewski>`). When the matrices involved in the
recursive expansion are sufficiently sparse, the calculation of the
density matrix scales linearly as a function of the system size O(N).
Another important feature is the extended Lagrangian framework for
Born-Oppenheimer molecular dynamics (XL-BOMD)
(:ref:`Niklasson2008 <Niklasson2008>`) (:ref:`Niklasson2014 <Niklasson2014>`),
(:ref:`Niklasson2017 <Niklasson2017>`) that allows for a drastic reduction
or even a complete removal of the iterative self-consistent field
optimization. Often only a single density matrix calculation per
molecular dynamics time step is required, yet total energy stability
is well maintained. The SP2 and XL-BOMD techniques enables stable
linear scaling MD simulations with a very small computational
overhead. This opens a number of opportunities in many different
areas of chemistry and materials science, as we now can simulate
larger system sizes and longer time scales
(:ref:`Cawkwell2012 <Cawkwell2012>`), (:ref:`Negre2016 <Negre2016>`).
----------
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files
<restart>`.
The :doc:`fix_modify <fix_modify>` *energy* option is supported by
this fix to add the potential energy computed by LATTE to the global
potential energy of the system as part of :doc:`thermodynamic output
<thermo_style>`. The default setting for this fix is :doc:`fix_modify
energy yes <fix_modify>`.
The :doc:`fix_modify <fix_modify>` *virial* option is supported by
this fix to add the contribution computed by LATTE to the global
pressure of the system as part of :doc:`thermodynamic output
<thermo_style>`. The default setting for this fix is :doc:`fix_modify
virial yes <fix_modify>`.
This fix computes a global scalar which can be accessed by various
:doc:`output commands <Howto_output>`. The scalar is the potential
energy discussed above. The scalar value calculated by this fix is
"extensive".
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command.
The DFTB forces computed by LATTE via this fix are used during an
energy minimization, invoked by the :doc:`minimize <minimize>`
command.
.. note::
If you want the potential energy associated with the DFTB
forces to be included in the total potential energy of the system (the
quantity being minimized), you MUST not disable the
:doc:`fix_modify <fix_modify>` *energy* option for this fix.
Restrictions
""""""""""""
This fix is part of the LATTE package. It is only enabled if LAMMPS
was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
You must use metal units, as set by the :doc:`units <units>` command to
use this fix.
LATTE does not currently compute per-atom energy or per-atom virial
contributions. So they will not show up as part of the calculations
performed by the :doc:`compute pe/atom <compute_pe_atom>` or
:doc:`compute stress/atom <compute_stress_atom>` commands.
Related commands
""""""""""""""""
none
Default
"""""""
none
----------
.. _Elstner:
**(Elstner)** M. Elstner, D. Poresag, G. Jungnickel, J. Elsner,
M. Haugk, T. Frauenheim, S. Suhai, and G. Seifert, Phys. Rev. B, 58,
7260 (1998).
.. _Elstner1:
**(Elstner)** M. Elstner, D. Poresag, G. Jungnickel, J. Elsner,
M. Haugk, T. Frauenheim, S. Suhai, and G. Seifert, Phys. Rev. B, 58,
7260 (1998).
.. _Finnis2:
**(Finnis)** M. W. Finnis, A. T. Paxton, M. Methfessel, and M. van
Schilfgarde, Phys. Rev. Lett., 81, 5149 (1998).
.. _Mniszewski:
**(Mniszewski)** S. M. Mniszewski, M. J. Cawkwell, M. E. Wall,
J. Mohd-Yusof, N. Bock, T. C. Germann, and A. M. N. Niklasson,
J. Chem. Theory Comput., 11, 4644 (2015).
.. _Niklasson2002:
**(Niklasson2002)** A. M. N. Niklasson, Phys. Rev. B, 66, 155115 (2002).
.. _Rubensson:
**(Rubensson)** E. H. Rubensson, A. M. N. Niklasson, SIAM
J. Sci. Comput. 36 (2), 147-170, (2014).
.. _Niklasson2008:
**(Niklasson2008)** A. M. N. Niklasson, Phys. Rev. Lett., 100, 123004
(2008).
.. _Niklasson2014:
**(Niklasson2014)** A. M. N. Niklasson and M. Cawkwell, J. Chem. Phys.,
141, 164123, (2014).
.. _Niklasson2017:
**(Niklasson2017)** A. M. N. Niklasson, J. Chem. Phys., 147, 054103 (2017).
.. _Cawkwell2012:
**(Cawkwell2012)** A. M. N. Niklasson, M. J. Cawkwell, Phys. Rev. B, 86
(17), 174308 (2012).
.. _Negre2016:
**(Negre2016)** C. F. A. Negre, S. M. Mniszewski, M. J. Cawkwell,
N. Bock, M. E. Wall, and A. M. N. Niklasson, J. Chem. Theory Comp.,
12, 3063 (2016).

1
doc/utils/sphinx-config/false_positives.txt
View File

@ -3545,6 +3545,7 @@ ters
tersoff
Tersoff
tesselation
tesselations
Tetot
tex
tfac

1
examples/README
View File

@ -83,7 +83,6 @@ hugoniostat: Hugoniostat shock dynamics
hyper: global and local hyperdynamics of diffusion on Pt surface
indent: spherical indenter into a 2d solid
kim: use of potentials in Knowledge Base for Interatomic Models (KIM)
latte: use of LATTE density-functional tight-binding quantum code
mc: MC package models: GCMC, Widom, fix mol/swap
mdi: use of the MDI package and MolSSI MDI code coupling library
meam: MEAM test for SiC and shear (same as shear examples)

36
examples/latte/TBparam/bondints.nonortho
View File

@ -1,36 +0,0 @@
Noints= 34
Element1 Element2 Kind H0 B1 B2 B3 B4 B5 R1 Rcut H0 B1 B2 B3 B4 B5 R1 Rcut
N O sss -11.430028 -2.257346 -1.152844 0.000000 0.000000 1.200000 3.500000 4.000000 0.340064 -1.703613 -0.622348 0.036738 -0.040158 1.200000 3.500000 4.000000
N O sps 11.597479 -1.382001 -0.765170 0.000000 0.000000 1.200000 3.500000 4.000000 -0.370946 -1.040947 -0.931097 0.252441 -0.115450 1.200000 3.500000 4.000000
O N sps 12.143744 -0.822913 -0.676127 0.000000 0.000000 1.200000 3.500000 4.000000 -0.420014 -1.107918 -0.905594 0.188424 -0.088365 1.200000 3.500000 4.000000
N O pps 9.465191 -1.082032 -0.769214 0.000000 0.000000 1.200000 3.500000 4.000000 -0.314073 0.499050 -2.914288 2.067657 -0.738439 1.200000 3.500000 4.000000
N O ppp -4.676789 -2.171480 -0.288002 0.000000 0.000000 1.200000 3.500000 4.000000 0.223937 -1.991867 -0.537630 -0.081270 -0.004130 1.200000 3.500000 4.000000
C O sss -14.369472 -2.077439 -0.875471 0.000000 0.000000 1.200000 3.500000 4.000000 0.375339 -1.547372 -0.642492 0.020614 -0.026699 1.200000 3.500000 4.000000
C O sps 9.576296 -1.156217 -0.494803 0.000000 0.000000 1.200000 3.500000 4.000000 -0.373027 -0.776043 -1.019920 0.257539 -0.102838 1.200000 3.500000 4.000000
O C sps 14.037374 -1.192632 -0.654572 0.000000 0.000000 1.200000 3.500000 4.000000 -0.458068 -1.035067 -0.937868 0.190562 -0.077841 1.200000 3.500000 4.000000
C O pps 9.331152 -0.718120 -0.822100 0.000000 0.000000 1.200000 3.500000 4.000000 -0.322293 0.795473 -3.476601 2.589965 -0.897800 1.200000 3.500000 4.000000
C O ppp -5.334367 -2.263939 -0.204910 0.000000 0.000000 1.200000 3.500000 4.000000 0.244570 -1.922717 -0.573671 -0.057280 -0.004108 1.200000 3.500000 4.000000
C N sss -7.010061 -1.730597 -0.575559 0.000000 0.000000 1.500000 3.500000 4.000000 0.263438 -1.754525 -0.584215 -0.007801 -0.021729 1.500000 3.500000 4.000000
C N sps 7.543283 -1.293768 -0.624363 0.000000 0.000000 1.500000 3.500000 4.000000 -0.326609 -1.197485 -0.807786 0.134891 -0.084373 1.500000 3.500000 4.000000
N C sps 9.090970 -1.494255 -0.616711 0.000000 0.000000 1.500000 3.500000 4.000000 -0.337943 -1.335442 -0.769693 0.119373 -0.079493 1.500000 3.500000 4.000000
C N pps 6.892240 -0.931920 -0.769164 0.000000 0.000000 1.500000 3.500000 4.000000 -0.350240 -0.467439 -1.849316 1.854403 -0.988471 1.500000 3.500000 4.000000
C N ppp -2.903346 -2.149349 -0.253006 0.000000 0.000000 1.500000 3.500000 4.000000 0.158424 -2.114409 -0.582346 -0.051076 -0.006183 1.500000 3.500000 4.000000
C C sss -9.404207 -1.363297 -0.507128 0.000000 0.000000 1.400000 3.500000 4.000000 0.346977 -1.519820 -0.570812 -0.013518 -0.015829 1.400000 3.500000 4.000000
C C sps 8.662429 -1.047410 -0.661999 0.000000 0.000000 1.400000 3.500000 4.000000 -0.400467 -0.984048 -0.853949 0.157178 -0.073381 1.400000 3.500000 4.000000
C C pps 6.811512 -0.552299 -0.776890 0.000000 0.000000 1.400000 3.500000 4.000000 -0.382417 0.102889 -2.786680 2.646356 -1.134320 1.400000 3.500000 4.000000
C C ppp -3.550127 -1.925572 -0.132715 0.000000 0.000000 1.400000 3.500000 4.000000 0.214357 -1.948923 -0.578323 -0.034356 -0.007257 1.400000 3.500000 4.000000
H C sss -9.072577 -1.393093 -0.430611 0.000000 0.000000 1.100000 3.500000 4.000000 0.416003 -1.459596 -0.654874 0.009140 -0.012658 1.100000 3.500000 4.000000
H C sps 8.176008 -0.985177 -0.427403 0.000000 0.000000 1.100000 3.500000 4.000000 -0.495695 -0.901626 -1.007214 0.189808 -0.057087 1.100000 3.500000 4.000000
H H sss -9.340000 -1.145903 -0.391777 0.000000 0.000000 0.750000 3.500000 4.000000 0.575007 -1.391261 -0.778831 0.080209 -0.017759 0.750000 3.500000 4.000000
O O sss -12.737687 -1.851608 -0.666621 0.000000 0.000000 1.200000 3.500000 4.000000 0.296445 -1.911896 -0.663451 0.038054 -0.046608 1.200000 3.500000 4.000000
O O sps 13.683050 -1.684554 -0.468349 0.000000 0.000000 1.200000 3.500000 4.000000 -0.362143 -1.285274 -0.939591 0.204641 -0.106438 1.200000 3.500000 4.000000
O O pps 9.460772 -1.211748 -0.581016 0.000000 0.000000 1.200000 3.500000 4.000000 -0.312044 0.121814 -2.519352 1.681266 -0.644566 1.200000 3.500000 4.000000
O O ppp -4.494595 -2.709223 -0.284124 0.000000 0.000000 1.200000 3.500000 4.000000 0.193010 -2.168462 -0.580629 -0.105104 0.004891 1.200000 3.500000 4.000000
H O sss -12.230931 -1.808632 -0.421164 0.000000 0.000000 1.000000 3.500000 4.000000 0.404725 -1.702546 -0.707938 0.074904 -0.039922 1.000000 3.500000 4.000000
H O sps 9.466088 -1.321262 -0.386336 0.000000 0.000000 1.000000 3.500000 4.000000 -0.447660 -0.952979 -1.163537 0.400616 -0.156965 1.000000 3.500000 4.000000
N N sss -7.710330 -2.365312 -0.525527 0.000000 0.000000 1.500000 3.500000 4.000000 0.231654 -1.879002 -0.572765 -0.004579 -0.031106 1.500000 3.500000 4.000000
N N sps 8.222314 -1.612118 -0.690081 0.000000 0.000000 1.500000 3.500000 4.000000 -0.305271 -1.385158 -0.751032 0.114531 -0.090839 1.500000 3.500000 4.000000
N N pps 7.178570 -1.176467 -0.571049 0.000000 0.000000 1.500000 3.500000 4.000000 -0.324668 -0.547805 -1.638658 1.495168 -0.827868 1.500000 3.500000 4.000000
N N ppp -2.829344 -2.408049 -0.387709 0.000000 0.000000 1.500000 3.500000 4.000000 0.142909 -2.162036 -0.571942 -0.071640 -0.004682 1.500000 3.500000 4.000000
H N sss -12.095890 -1.519057 -0.277247 0.000000 0.000000 1.000000 3.500000 4.000000 0.446693 -1.500463 -0.657448 0.065741 -0.037004 1.000000 3.500000 4.000000
H N sps 9.851338 -1.231616 -0.370836 0.000000 0.000000 1.000000 3.500000 4.000000 -0.501530 -0.785734 -1.123232 0.394878 -0.148501 1.000000 3.500000 4.000000

7
examples/latte/TBparam/electrons.dat
View File

@ -1,7 +0,0 @@
Noelem= 5
Element basis Numel Es Ep Ed Ef Mass HubbardU Wss Wpp Wdd Wff
N sp 5.0 -18.58 -7.09 0.0 0.0 14.0067 15.93 0.0 -0.6950 0.0 0.0
O sp 6.0 -23.96 -9.02 0.0 0.0 15.9994 12.15 0.0 -0.7577 0.0 0.0
H s 1.0 -6.35 0.0 0.0 0.0 1.0079 12.85 -1.7937 0.0 0.0 0.0
C sp 4.0 -13.75 -5.28 0.0 0.0 12.01 10.0 0.0 -0.621 0.0 0.0
Ti sd 4.0 -5.5 0.0 -3.0 0.0 47.867 10.0 0.0 0.0 0.0 0.0

12
examples/latte/TBparam/ppots.nonortho
View File

@ -1,12 +0,0 @@
Nopps= 10
Ele1 Ele2 A0 A1 A2 A3 A4 A5 A6 C R1 Rcut
N O 13.182426 20.050322 -46.806321 38.206953 -12.319656 0.000000 0.000000 0.000000 1.600000 1.700000
C N 88.953762 10.294988 -27.706877 22.101434 -6.836438 0.000000 0.000000 0.000000 1.600000 1.700000
C O 0.944093 30.116337 -59.608215 45.107654 -13.178839 0.000000 0.000000 0.000000 1.600000 1.700000
C H 104.889589 3.971095 -23.823043 26.408093 -11.317522 0.000000 0.000000 0.000000 1.200000 1.300000
C C 3.962931 24.467772 -51.156024 39.031644 -11.342979 0.000000 0.000000 0.000000 1.600000 1.700000
H H 38.512100 3.887860 -37.769100 57.083500 -34.512200 0.000000 0.000000 0.000000 0.900000 1.000000
N N 43.228899 15.004605 -36.621777 29.234888 -8.912743 0.000000 0.000000 0.000000 1.600000 1.700000
N H 0.625470 28.081241 -63.414297 53.286361 -17.352234 0.000000 0.000000 0.000000 1.300000 1.400000
O O 10.999870 19.303033 -45.747853 37.946431 -11.935755 0.000000 0.000000 0.000000 1.500000 1.600000
O H 0.481176 33.175383 -81.158683 74.935408 -26.792315 0.000000 0.000000 0.000000 1.200000 1.300000

23
examples/latte/data.ch4
View File

@ -1,23 +0,0 @@
LAMMPS Description
5 atoms
2 atom types
0.0000000000000000 19.523000000000000 xlo xhi
0.0000000000000000 12.757999999999999 ylo yhi
0.0000000000000000 11.692000000000000 zlo zhi
0.0000000000000000 0.0000000000000000 0.0000000000000000 xy xz yz
Masses
1 12.010000000000000
2 1.0078250169754028
Atoms
1 1 1 0.0 -9.16600 2.05200 0.00000
2 1 2 0.0 -8.09600 2.05200 0.00000
3 1 2 0.0 -9.52300 2.75800 -0.72000
4 1 2 0.0 -9.52300 2.32200 0.97200
5 1 2 0.0 -9.52300 1.07500 -0.25200

49
examples/latte/data.graphene
View File

@ -1,49 +0,0 @@
LAMMPS Description
32 atoms
1 atom types
0.0000000000000000 10.000000000000000 xlo xhi
0.0000000000000000 8.0000000000000000 ylo yhi
0.0000000000000000 20.000000000000000 zlo zhi
4.8985871965894128E-016 1.2246467991473533E-015 1.2246467991473533E-015 xy xz yz
Masses
1 12.010000000000000
Atoms
1 1 1 0.0 4.93100 4.25000 0.00500
2 1 1 0.0 8.62100 2.12100 0.14000
3 1 1 0.0 3.70700 2.12600 0.14700
4 1 1 0.0 7.38200 4.25400 0.07800
5 1 1 0.0 2.47900 4.25400 0.08000
6 1 1 0.0 6.15800 6.37400 -0.01000
7 1 1 0.0 1.23700 6.38300 0.06600
8 1 1 0.0 1.24000 2.12100 0.14600
9 1 1 0.0 6.15500 2.12600 0.12900
10 1 1 0.0 0.00700 4.25200 0.12200
11 1 1 0.0 8.62100 6.38500 0.04100
12 1 1 0.0 3.70000 6.37400 -0.01000
13 1 1 0.0 0.00600 1.41600 0.13000
14 1 1 0.0 4.93000 1.40800 0.14700
15 1 1 0.0 8.61800 3.54600 0.11500
16 1 1 0.0 3.70800 3.55300 0.08400
17 1 1 0.0 7.39400 5.68000 0.03500
18 1 1 0.0 2.46500 5.68000 0.03500
19 1 1 0.0 6.16000 7.80500 0.02700
20 1 1 0.0 1.23800 7.81100 0.06000
21 1 1 0.0 2.47300 1.41800 0.16100
22 1 1 0.0 7.38900 1.41700 0.14800
23 1 1 0.0 1.24200 3.54700 0.12600
24 1 1 0.0 6.15300 3.55300 0.07400
25 1 1 0.0 0.00700 5.67800 0.09700
26 1 1 0.0 4.93100 5.66800 -0.03100
27 1 1 0.0 8.62000 7.81300 0.03900
28 1 1 0.0 3.70100 7.80200 0.03700
29 1 1 0.0 0.00700 -0.01000 0.08900
30 1 1 0.0 4.93100 -0.01500 0.16100
31 1 1 0.0 2.47300 -0.01200 0.14400
32 1 1 0.0 7.38900 -0.01300 0.14800

63
examples/latte/data.sucrose
View File

@ -1,63 +0,0 @@
LAMMPS Description
45 atoms
3 atom types
0.0000000000000000 17.202999999999999 xlo xhi
0.0000000000000000 18.009000000000000 ylo yhi
0.0000000000000000 21.643000000000001 zlo zhi
Masses
1 15.9994
2 12.01
3 1.0079
Atoms
1 1 1 0.0 8.62700 8.66700 12.48600
2 1 1 0.0 9.11200 9.11800 10.27300
3 1 1 0.0 8.45700 11.33100 10.49000
4 1 1 0.0 11.72600 8.36500 10.66700
5 1 1 0.0 8.06500 8.99400 7.93600
6 1 2 0.0 9.62800 9.07200 11.59100
7 1 3 0.0 9.90900 10.08300 11.89200
8 1 2 0.0 9.07000 10.40400 9.64000
9 1 1 0.0 6.14600 11.61000 8.00500
10 1 1 0.0 11.07200 10.13000 8.37600
11 1 1 0.0 6.10200 10.00900 11.62100
12 1 2 0.0 8.14000 10.29100 8.45100
13 1 3 0.0 8.49000 10.91200 7.62300
14 1 1 0.0 7.41500 7.08400 14.43400
15 1 2 0.0 10.75100 8.07000 11.65100
16 1 3 0.0 11.24000 8.11800 12.63400
17 1 2 0.0 7.09000 11.63400 10.17000
18 1 3 0.0 7.06900 12.69800 9.91100
19 1 2 0.0 7.97200 7.44200 12.14000
20 1 3 0.0 7.54700 7.58300 11.13800
21 1 1 0.0 11.24900 5.73000 11.78600
22 1 2 0.0 10.26800 6.65900 11.37300
23 1 3 0.0 10.12300 6.53400 10.29200
24 1 2 0.0 6.78400 10.79500 8.95500
25 1 3 0.0 6.12100 9.95500 9.19600
26 1 2 0.0 10.47500 10.88300 9.39800
27 1 3 0.0 10.49500 11.92100 9.06900
28 1 3 0.0 11.09100 10.82000 10.30900
29 1 2 0.0 8.99100 6.32000 12.11700
30 1 3 0.0 9.23100 6.01100 13.14400
31 1 2 0.0 6.86600 7.25300 13.14500
32 1 3 0.0 6.17700 8.10100 13.15700
33 1 3 0.0 6.28900 6.35300 12.94300
34 1 2 0.0 6.24000 11.39400 11.39300
35 1 3 0.0 6.66500 11.86500 12.28300
36 1 3 0.0 5.23100 11.78100 11.26000
37 1 1 0.0 8.34300 5.24100 11.48000
38 1 3 0.0 12.00100 9.28600 10.78200
39 1 3 0.0 12.06300 5.97500 11.33000
40 1 3 0.0 6.99600 9.67600 11.79700
41 1 3 0.0 7.93700 7.87600 14.60900
42 1 3 0.0 10.95500 9.19800 8.60700
43 1 3 0.0 5.94400 11.05900 7.24100
44 1 3 0.0 7.94900 8.39500 8.68400
45 1 3 0.0 8.96400 4.50300 11.48800

41
examples/latte/data.water
View File

@ -1,41 +0,0 @@
LAMMPS Description
24 atoms
2 atom types
0.0000000000000000 6.2670000000000003 xlo xhi
0.0000000000000000 6.2670000000000003 ylo yhi
0.0000000000000000 6.2670000000000003 zlo zhi
Masses
1 15.994915008544922
2 1.0078250169754028
Atoms
1 1 1 0.0 3.08800 3.70000 3.12400
2 1 2 0.0 4.05800 3.70000 3.12400
3 1 2 0.0 2.76400 3.13200 3.84100
4 1 1 0.0 2.47000 0.39000 1.36000
5 1 2 0.0 1.54000 0.37000 1.73000
6 1 2 0.0 2.48000 0.00000 0.44000
7 1 1 0.0 1.99300 0.41700 5.25000
8 1 2 0.0 2.39300 1.32700 5.16000
9 1 2 0.0 0.99300 0.49700 5.31000
10 1 1 0.0 2.05300 6.09700 3.48000
11 1 2 0.0 2.12300 5.20700 3.02000
12 1 2 0.0 1.11300 0.17000 3.40000
13 1 1 0.0 4.90000 5.37700 2.14000
14 1 2 0.0 5.51000 6.17700 2.18000
15 1 2 0.0 3.95000 5.68700 2.21000
16 1 1 0.0 0.92000 3.82700 0.56000
17 1 2 0.0 0.00000 3.54700 0.27000
18 1 2 0.0 1.23000 4.59700 0.00000
19 1 1 0.0 0.89000 2.03700 3.41000
20 1 2 0.0 0.72000 2.86700 2.87000
21 1 2 0.0 1.79000 1.66700 3.19000
22 1 1 0.0 4.45000 4.61700 5.43000
23 1 2 0.0 4.75000 3.89700 4.81000
24 1 2 0.0 4.06000 4.21700 6.26000

44
examples/latte/in.latte.graphene.boxrelax
View File

@ -1,44 +0,0 @@
# Simple water model with LATTE
units metal
atom_style full
atom_modify sort 0 0.0 # turn off sorting of the coordinates
read_data data.graphene
# replicate system if requested
variable x index 1
variable y index 1
variable z index 1
variable nrep equal v_x*v_y*v_z
if "${nrep} > 1" then "replicate $x $y $z"
# initialize system
velocity all create 0.0 87287 loop geom
pair_style zero 1.0
pair_coeff * *
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.00025
fix 1 all box/relax iso 0.0 vmax 0.001
fix 2 all latte
fix_modify 2 energy yes
thermo_style custom etotal
# minimization
thermo 1
fix 3 all print 1 "Total Energy ="
min_style cg
min_modify dmax 0.1
min_modify line quadratic
minimize 1.0e-4 1.0e-4 10000 10000

65
examples/latte/in.latte.multiple
View File

@ -1,65 +0,0 @@
units metal
atom_style full
atom_modify sort 0 0.0 # turn off sorting of the coordinates
read_data data.water
# initialize system
velocity all create 0.0 87287 loop geom
pair_style zero 1.0
pair_coeff * *
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.00025
fix 1 all nve
fix 2 all latte
fix_modify 2 energy yes
thermo_style custom step temp pe etotal press
# dynamics
thermo 10
run 10
# Clear up previus calculation
clear
# simple CH4 molecule with LATTE
units metal
atom_style full
atom_modify sort 0 0.0 # turn off sorting of the coordinates
read_data data.ch4
# initialize system
velocity all create 0.0 87287 loop geom
pair_style zero 1.0
pair_coeff * *
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.00025
fix 1 all nve
fix 2 all latte
fix_modify 2 energy yes
thermo_style custom step temp pe etotal press
# dynamics
thermo 10
run 10

40
examples/latte/in.latte.sucrose.md
View File

@ -1,40 +0,0 @@
# simple sucrose model with LATTE
units metal
atom_style full
atom_modify sort 0 0.0 # turn off sorting of the coordinates
read_data data.sucrose
# replicate system if requested
variable x index 1
variable y index 1
variable z index 1
variable nrep equal v_x*v_y*v_z
if "${nrep} > 1" then "replicate $x $y $z"
# initialize system
velocity all create 0.0 87287 loop geom
pair_style zero 1.0
pair_coeff * *
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.00025
fix 1 all nve
fix 2 all latte
fix_modify 2 energy yes
thermo_style custom step temp pe etotal press
# dynamics
thermo 10
run 100

40
examples/latte/in.latte.water.md
View File

@ -1,40 +0,0 @@
# simple water model with LATTE
units metal
atom_style full
atom_modify sort 0 0.0 # turn off sorting of the coordinates
read_data data.water
# replicate system if requested
variable x index 1
variable y index 1
variable z index 1
variable nrep equal v_x*v_y*v_z
if "${nrep} > 1" then "replicate $x $y $z"
# initialize system
velocity all create 0.0 87287 loop geom
pair_style zero 1.0
pair_coeff * *
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.00025
fix 1 all nve
fix 2 all latte
fix_modify 2 energy yes
thermo_style custom step temp pe etotal press
# dynamics
thermo 10
run 100

42
examples/latte/in.latte.water.min
View File

@ -1,42 +0,0 @@
# simple water model with LATTE
units metal
atom_style full
atom_modify sort 0 0.0 # turn off sorting of the coordinates
read_data data.water
# replicate system if requested
variable x index 1
variable y index 1
variable z index 1
variable nrep equal v_x*v_y*v_z
if "${nrep} > 1" then "replicate $x $y $z"
# initialize system
velocity all create 0.0 87287 loop geom
pair_style zero 1.0
pair_coeff * *
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.00025
fix 1 all nve
fix 2 all latte
fix_modify 2 energy yes
thermo_style custom step temp pe etotal press
# minimization
thermo 10
min_style fire
minimize 1.0e-4 1.0e-4 500 500

32
examples/latte/latte.in
View File

@ -1,32 +0,0 @@
LATTE INPUT FILE
================
#This input file resumes the content of MDcontroller and TBparam/control.in
#The parser will only read it is present inside the running folder.
#In case this file is not present Latte will read the two files as original.
#The order of the keywords is not important in this file.
#To get a full description of these keywords please see:
## https://github.com/lanl/LATTE/blob/master/Manual/LATTE_manual.pdf
#General controls
CONTROL{
XCONTROL= 1
BASISTYPE= NONORTHO
PARAMPATH= "./TBparam"
KBT= 0.0
ENTROPYKIND= 1
PPOTON= 1
SPINON= 0 SPINTOL= 1.0e-4
ELECTRO= 1 ELECMETH= 0 ELEC_QTOL= 1.0e-8
MAXSCF= 450
BREAKTOL= 1.0E-6 MINSP2ITER= 22 SP2CONV= REL
FULLQCONV= 1 QITER= 3
QMIX= 0.25 SPINMIX= 0.25 MDMIX= 0.25
SPARSEON= 0 THRESHOLDON= 1 NUMTHRESH= 1.0e-6 FILLINSTOP= 100 BLKSZ= 4
MSPARSE= 1500
SKIN= 1.0
CHARGE= 0
XBO= 1
XBODISON= 1
XBODISORDER= 5
KON= 0
}

179
examples/latte/log.13Sep22.latte.graphene.boxrelax.g++.1
View File

@ -1,179 +0,0 @@
LAMMPS (3 Aug 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
# Simple water model with LATTE
units metal
atom_style full
atom_modify sort 0 0.0 # turn off sorting of the coordinates
read_data data.graphene
Reading data file ...
triclinic box = (0 0 0) to (10 8 20) with tilt (4.8985872e-16 1.2246468e-15 1.2246468e-15)
1 by 1 by 1 MPI processor grid
reading atoms ...
32 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000 seconds
read_data CPU = 0.001 seconds
# replicate system if requested
variable x index 1
variable y index 1
variable z index 1
variable nrep equal v_x*v_y*v_z
if "${nrep} > 1" then "replicate $x $y $z"
# initialize system
velocity all create 0.0 87287 loop geom
pair_style zero 1.0
pair_coeff * *
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.00025
fix 1 all box/relax iso 0.0 vmax 0.001
fix 2 all latte
fix_modify 2 energy yes
thermo_style custom etotal
# minimization
thermo 1
fix 3 all print 1 "Total Energy ="
min_style cg
min_modify dmax 0.1
min_modify line quadratic
minimize 1.0e-4 1.0e-4 10000 10000
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2
ghost atom cutoff = 2
binsize = 1, bins = 11 9 20
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair zero, perpetual
attributes: half, newton on
pair build: half/bin/newton/tri
stencil: half/bin/3d/tri
bin: standard
Per MPI rank memory allocation (min/avg/max) = 7.03 | 7.03 | 7.03 Mbytes
TotEng
-247.46002
-247.67224
-247.87937
-248.08148
-248.27865
-248.47096
-248.65851
-248.84137
-249.01964
-249.19342
-249.36281
-249.52791
-249.68883
-249.8457
-249.99865
-250.1478
-250.29332
-250.43535
-250.57409
-250.70972
-250.84247
-250.97258
-251.10035
-251.2261
-251.35021
-251.47314
-251.59543
-251.71776
-251.84096
-251.9661
-252.09459
-252.22833
-252.37003
-252.52371
-252.69578
-252.89752
-253.15197
-253.52044
-254.31418
-255.6175
-256.8162
-258.1227
-259.38401
-260.74831
-262.03991
-263.5463
-264.70486
-267.69143
-267.88682
-269.0352
-270.602
-270.65395
-270.7429
-271.55831
-271.81159
-271.87447
-273.03096
-273.23109
-273.27869
-273.34621
-273.4082
-273.45599
-273.53849
-273.57478
-273.71381
-273.74092
Loop time of 20.5496 on 1 procs for 65 steps with 32 atoms
99.7% CPU use with 1 MPI tasks x 1 OpenMP threads
Minimization stats:
Stopping criterion = energy tolerance
Energy initial, next-to-last, final =
-247.460020562055 -273.713813242259 -273.740918498854
Force two-norm initial, final = 201.60784 9.4927634
Force max component initial, final = 188.92406 2.4327308
Final line search alpha, max atom move = 0.00022885545 0.0005567437
Iterations, force evaluations = 65 65
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 7.7869e-05 | 7.7869e-05 | 7.7869e-05 | 0.0 | 0.00
Bond | 3.531e-06 | 3.531e-06 | 3.531e-06 | 0.0 | 0.00
Neigh | 1.3988e-05 | 1.3988e-05 | 1.3988e-05 | 0.0 | 0.00
Comm | 0.00014355 | 0.00014355 | 0.00014355 | 0.0 | 0.00
Output | 0.00071475 | 0.00071475 | 0.00071475 | 0.0 | 0.00
Modify | 20.547 | 20.547 | 20.547 | 0.0 | 99.99
Other | | 0.001683 | | | 0.01
Nlocal: 32 ave 32 max 32 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 100 ave 100 max 100 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 48 ave 48 max 48 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 48
Ave neighs/atom = 1.5
Ave special neighs/atom = 0
Neighbor list builds = 1
Dangerous builds = 0
Total wall time: 0:00:20

189
examples/latte/log.13Sep22.latte.multiple.g++.1
View File

@ -1,189 +0,0 @@
LAMMPS (3 Aug 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
units metal
atom_style full
atom_modify sort 0 0.0 # turn off sorting of the coordinates
read_data data.water
Reading data file ...
orthogonal box = (0 0 0) to (6.267 6.267 6.267)
1 by 1 by 1 MPI processor grid
reading atoms ...
24 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000 seconds
read_data CPU = 0.001 seconds
# initialize system
velocity all create 0.0 87287 loop geom
pair_style zero 1.0
pair_coeff * *
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.00025
fix 1 all nve
fix 2 all latte
fix_modify 2 energy yes
thermo_style custom step temp pe etotal press
# dynamics
thermo 10
run 10
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2
ghost atom cutoff = 2
binsize = 1, bins = 7 7 7
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair zero, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 5.88 | 5.88 | 5.88 Mbytes
Step Temp PotEng TotEng Press
0 0 -104.95596 -104.95596 48235.442
10 336.53107 -105.96027 -104.95977 97996.851
Loop time of 0.334108 on 1 procs for 10 steps with 24 atoms
Performance: 0.646 ns/day, 37.123 hours/ns, 29.930 timesteps/s
99.2% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 3.714e-06 | 3.714e-06 | 3.714e-06 | 0.0 | 0.00
Bond | 5.02e-07 | 5.02e-07 | 5.02e-07 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 1.1209e-05 | 1.1209e-05 | 1.1209e-05 | 0.0 | 0.00
Output | 2.3638e-05 | 2.3638e-05 | 2.3638e-05 | 0.0 | 0.01
Modify | 0.33404 | 0.33404 | 0.33404 | 0.0 | 99.98
Other | | 2.795e-05 | | | 0.01
Nlocal: 24 ave 24 max 24 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 71 ave 71 max 71 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 37 ave 37 max 37 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 37
Ave neighs/atom = 1.5416667
Ave special neighs/atom = 0
Neighbor list builds = 0
Dangerous builds = 0
# Clear up previus calculation
clear
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
# simple CH4 molecule with LATTE
units metal
atom_style full
atom_modify sort 0 0.0 # turn off sorting of the coordinates
read_data data.ch4
Reading data file ...
triclinic box = (0 0 0) to (19.523 12.758 11.692) with tilt (0 0 0)
1 by 1 by 1 MPI processor grid
reading atoms ...
5 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000 seconds
read_data CPU = 0.007 seconds
# initialize system
velocity all create 0.0 87287 loop geom
pair_style zero 1.0
pair_coeff * *
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.00025
fix 1 all nve
fix 2 all latte
fix_modify 2 energy yes
thermo_style custom step temp pe etotal press
# dynamics
thermo 10
run 10
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2
ghost atom cutoff = 2
binsize = 1, bins = 20 13 12
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair zero, perpetual
attributes: half, newton on
pair build: half/bin/newton/tri
stencil: half/bin/3d/tri
bin: standard
Per MPI rank memory allocation (min/avg/max) = 5.902 | 5.902 | 5.902 Mbytes
Step Temp PotEng TotEng Press
0 0 -23.980353 -23.980353 348.02716
10 19.123149 -23.990297 -23.98041 18.774332
Loop time of 0.0121573 on 1 procs for 10 steps with 5 atoms
Performance: 17.767 ns/day, 1.351 hours/ns, 822.549 timesteps/s
99.7% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.224e-06 | 1.224e-06 | 1.224e-06 | 0.0 | 0.01
Bond | 2.93e-07 | 2.93e-07 | 2.93e-07 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 3.845e-06 | 3.845e-06 | 3.845e-06 | 0.0 | 0.03
Output | 8.633e-06 | 8.633e-06 | 8.633e-06 | 0.0 | 0.07
Modify | 0.012132 | 0.012132 | 0.012132 | 0.0 | 99.80
Other | | 1.089e-05 | | | 0.09
Nlocal: 5 ave 5 max 5 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 7 ave 7 max 7 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 10 ave 10 max 10 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 10
Ave neighs/atom = 2
Ave special neighs/atom = 0
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:00

112
examples/latte/log.13Sep22.latte.sucrose.md.g++.1
View File

@ -1,112 +0,0 @@
LAMMPS (3 Aug 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
# simple sucrose model with LATTE
units metal
atom_style full
atom_modify sort 0 0.0 # turn off sorting of the coordinates
read_data data.sucrose
Reading data file ...
orthogonal box = (0 0 0) to (17.203 18.009 21.643)
1 by 1 by 1 MPI processor grid
reading atoms ...
45 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000 seconds
read_data CPU = 0.001 seconds
# replicate system if requested
variable x index 1
variable y index 1
variable z index 1
variable nrep equal v_x*v_y*v_z
if "${nrep} > 1" then "replicate $x $y $z"
# initialize system
velocity all create 0.0 87287 loop geom
pair_style zero 1.0
pair_coeff * *
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.00025
fix 1 all nve
fix 2 all latte
fix_modify 2 energy yes
thermo_style custom step temp pe etotal press
# dynamics
thermo 10
run 100
Generated 0 of 3 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2
ghost atom cutoff = 2
binsize = 1, bins = 18 19 22
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair zero, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 5.93 | 5.93 | 5.93 Mbytes
Step Temp PotEng TotEng Press
0 0 -251.26617 -251.26617 16.617233
10 0.025263738 -251.26631 -251.26617 8.0576369
20 0.034232485 -251.26636 -251.26617 1.6672772
30 0.059079585 -251.2665 -251.26617 11.058355
40 0.055499785 -251.26648 -251.26617 14.837599
50 0.058499496 -251.2665 -251.26617 6.7180488
60 0.071094531 -251.26657 -251.26617 6.6131215
70 0.084309398 -251.26665 -251.26617 12.372502
80 0.1089929 -251.26679 -251.26617 8.8352747
90 0.11378255 -251.26681 -251.26617 5.1175071
100 0.13003967 -251.26691 -251.26617 8.2429118
Loop time of 14.4456 on 1 procs for 100 steps with 45 atoms
Performance: 0.150 ns/day, 160.507 hours/ns, 6.923 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 7.5758e-05 | 7.5758e-05 | 7.5758e-05 | 0.0 | 0.00
Bond | 6.748e-06 | 6.748e-06 | 6.748e-06 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 9.0137e-05 | 9.0137e-05 | 9.0137e-05 | 0.0 | 0.00
Output | 0.00025976 | 0.00025976 | 0.00025976 | 0.0 | 0.00
Modify | 14.445 | 14.445 | 14.445 | 0.0 | 99.99
Other | | 0.0005283 | | | 0.00
Nlocal: 45 ave 45 max 45 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 59 ave 59 max 59 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 59
Ave neighs/atom = 1.3111111
Ave special neighs/atom = 0
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:14

112
examples/latte/log.13Sep22.latte.water.md.g++.1
View File

@ -1,112 +0,0 @@
LAMMPS (3 Aug 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
# simple water model with LATTE
units metal
atom_style full
atom_modify sort 0 0.0 # turn off sorting of the coordinates
read_data data.water
Reading data file ...
orthogonal box = (0 0 0) to (6.267 6.267 6.267)
1 by 1 by 1 MPI processor grid
reading atoms ...
24 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000 seconds
read_data CPU = 0.001 seconds
# replicate system if requested
variable x index 1
variable y index 1
variable z index 1
variable nrep equal v_x*v_y*v_z
if "${nrep} > 1" then "replicate $x $y $z"
# initialize system
velocity all create 0.0 87287 loop geom
pair_style zero 1.0
pair_coeff * *
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.00025
fix 1 all nve
fix 2 all latte
fix_modify 2 energy yes
thermo_style custom step temp pe etotal press
# dynamics
thermo 10
run 100
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2
ghost atom cutoff = 2
binsize = 1, bins = 7 7 7
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair zero, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 5.88 | 5.88 | 5.88 Mbytes
Step Temp PotEng TotEng Press
0 0 -104.95596 -104.95596 48235.442
10 336.53107 -105.96027 -104.95977 97996.851
20 529.06408 -106.53023 -104.95733 131519.85
30 753.62603 -107.19952 -104.959 49296.66
40 716.65648 -107.08803 -104.95742 28307.121
50 824.04392 -107.40823 -104.95836 102167.59
60 933.56146 -107.73479 -104.95933 92508.517
70 851.18489 -107.48767 -104.95711 13993.262
80 999.8028 -107.93147 -104.95907 36700.736
90 998.77488 -107.9257 -104.95636 107233.54
100 1281.4438 -108.76963 -104.95992 49702.386
Loop time of 3.14578 on 1 procs for 100 steps with 24 atoms
Performance: 0.687 ns/day, 34.953 hours/ns, 31.789 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 3.0818e-05 | 3.0818e-05 | 3.0818e-05 | 0.0 | 0.00
Bond | 4.704e-06 | 4.704e-06 | 4.704e-06 | 0.0 | 0.00
Neigh | 1.8668e-05 | 1.8668e-05 | 1.8668e-05 | 0.0 | 0.00
Comm | 0.00010831 | 0.00010831 | 0.00010831 | 0.0 | 0.00
Output | 0.00021087 | 0.00021087 | 0.00021087 | 0.0 | 0.01
Modify | 3.1452 | 3.1452 | 3.1452 | 0.0 | 99.98
Other | | 0.0002339 | | | 0.01
Nlocal: 24 ave 24 max 24 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 77 ave 77 max 77 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 31 ave 31 max 31 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 31
Ave neighs/atom = 1.2916667
Ave special neighs/atom = 0
Neighbor list builds = 2
Dangerous builds = 0
Total wall time: 0:00:03

122
examples/latte/log.13Sep22.latte.water.min.g++.1
View File

@ -1,122 +0,0 @@
LAMMPS (3 Aug 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
# simple water model with LATTE
units metal
atom_style full
atom_modify sort 0 0.0 # turn off sorting of the coordinates
read_data data.water
Reading data file ...
orthogonal box = (0 0 0) to (6.267 6.267 6.267)
1 by 1 by 1 MPI processor grid
reading atoms ...
24 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000 seconds
read_data CPU = 0.001 seconds
# replicate system if requested
variable x index 1
variable y index 1
variable z index 1
variable nrep equal v_x*v_y*v_z
if "${nrep} > 1" then "replicate $x $y $z"
# initialize system
velocity all create 0.0 87287 loop geom
pair_style zero 1.0
pair_coeff * *
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.00025
fix 1 all nve
fix 2 all latte
fix_modify 2 energy yes
thermo_style custom step temp pe etotal press
# minimization
thermo 10
min_style fire
minimize 1.0e-4 1.0e-4 500 500
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2
ghost atom cutoff = 2
binsize = 1, bins = 7 7 7
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair zero, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d
bin: standard
Parameters for fire:
dmax delaystep dtgrow dtshrink alpha0 alphashrink tmax tmin integrator halfstepback
0.1 20 1.1 0.5 0.25 0.99 10 0.02 eulerimplicit yes
Per MPI rank memory allocation (min/avg/max) = 5.88 | 5.88 | 5.88 Mbytes
Step Temp PotEng TotEng Press
0 0 -104.95596 -104.95596 48235.442
10 853.69689 -106.31143 -103.7734 79191.444
20 1112.0893 -107.2723 -103.96607 82675.468
30 1897.6249 -108.36769 -102.72608 71447.508
40 3068.3491 -110.06452 -100.94237 47627.967
50 3.730935 -110.16042 -110.14932 5913.0643
60 28.603141 -110.18885 -110.10381 5778.8586
66 54.717686 -110.21503 -110.05236 5739.5831
Loop time of 2.48723 on 1 procs for 66 steps with 24 atoms
99.5% CPU use with 1 MPI tasks x 1 OpenMP threads
Minimization stats:
Stopping criterion = energy tolerance
Energy initial, next-to-last, final =
-104.955957263186 -110.209885831179 -110.215033825672
Force two-norm initial, final = 19.119006 0.51695213
Force max component initial, final = 11.775801 0.1663917
Final line search alpha, max atom move = 0 0
Iterations, force evaluations = 66 69
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 2.7159e-05 | 2.7159e-05 | 2.7159e-05 | 0.0 | 0.00
Bond | 3.124e-06 | 3.124e-06 | 3.124e-06 | 0.0 | 0.00
Neigh | 1.0201e-05 | 1.0201e-05 | 1.0201e-05 | 0.0 | 0.00
Comm | 0.000109 | 0.000109 | 0.000109 | 0.0 | 0.00
Output | 0.00010568 | 0.00010568 | 0.00010568 | 0.0 | 0.00
Modify | 2.4866 | 2.4866 | 2.4866 | 0.0 | 99.98
Other | | 0.0003552 | | | 0.01
Nlocal: 24 ave 24 max 24 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 75 ave 75 max 75 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 33 ave 33 max 33 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 33
Ave neighs/atom = 1.375
Ave special neighs/atom = 0
Neighbor list builds = 1
Dangerous builds = 0
Total wall time: 0:00:02

171
examples/latte/log.21Jun18.latte.multiple.g++.1
View File

@ -1,171 +0,0 @@
LAMMPS (22 Jun 2018)
units metal
atom_style full
atom_modify sort 0 0.0 # turn off sorting of the coordinates
read_data data.water
orthogonal box = (0 0 0) to (6.267 6.267 6.267)
1 by 1 by 1 MPI processor grid
reading atoms ...
24 atoms
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
# initialize system
velocity all create 0.0 87287 loop geom
pair_style zero 1.0
pair_coeff * *
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.00025
fix 1 all nve
fix 2 all latte NULL
fix_modify 2 energy yes
thermo_style custom step temp pe etotal press
# dynamics
thermo 10
run 10
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2
ghost atom cutoff = 2
binsize = 1, bins = 7 7 7
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair zero, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 5.629 | 5.629 | 5.629 Mbytes
Step Temp PotEng TotEng Press
0 0 -104.95596 -104.95596 48235.442
10 336.53107 -105.96027 -104.95977 97996.851
Loop time of 0.912921 on 1 procs for 10 steps with 24 atoms
Performance: 0.237 ns/day, 101.436 hours/ns, 10.954 timesteps/s
6614.5% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 4.3392e-05 | 4.3392e-05 | 4.3392e-05 | 0.0 | 0.00
Bond | 3.8862e-05 | 3.8862e-05 | 3.8862e-05 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.00036955 | 0.00036955 | 0.00036955 | 0.0 | 0.04
Output | 9.799e-05 | 9.799e-05 | 9.799e-05 | 0.0 | 0.01
Modify | 0.91199 | 0.91199 | 0.91199 | 0.0 | 99.90
Other | | 0.0003812 | | | 0.04
Nlocal: 24 ave 24 max 24 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 71 ave 71 max 71 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 37 ave 37 max 37 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 37
Ave neighs/atom = 1.54167
Ave special neighs/atom = 0
Neighbor list builds = 0
Dangerous builds = 0
# Clear up previus calculation
clear
# simple CH4 molecule with LATTE
units metal
atom_style full
atom_modify sort 0 0.0 # turn off sorting of the coordinates
read_data data.ch4
triclinic box = (0 0 0) to (19.523 12.758 11.692) with tilt (0 0 0)
1 by 1 by 1 MPI processor grid
reading atoms ...
5 atoms
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
# initialize system
velocity all create 0.0 87287 loop geom
pair_style zero 1.0
pair_coeff * *
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
timestep 0.00025
fix 1 all nve
fix 2 all latte NULL
fix_modify 2 energy yes
thermo_style custom step temp pe etotal press
# dynamics
thermo 10
run 10
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2
ghost atom cutoff = 2
binsize = 1, bins = 20 13 12
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair zero, perpetual
attributes: half, newton on
pair build: half/bin/newton/tri
stencil: half/bin/3d/newton/tri
bin: standard
Per MPI rank memory allocation (min/avg/max) = 5.651 | 5.651 | 5.651 Mbytes
Step Temp PotEng TotEng Press
0 0 -23.980353 -23.980353 348.02716
10 19.123149 -23.990297 -23.98041 18.774332
Loop time of 0.0415399 on 1 procs for 10 steps with 5 atoms
Performance: 5.200 ns/day, 4.616 hours/ns, 240.732 timesteps/s
6674.9% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 5.722e-06 | 5.722e-06 | 5.722e-06 | 0.0 | 0.01
Bond | 7.6294e-06 | 7.6294e-06 | 7.6294e-06 | 0.0 | 0.02
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 1.8358e-05 | 1.8358e-05 | 1.8358e-05 | 0.0 | 0.04
Output | 2.0981e-05 | 2.0981e-05 | 2.0981e-05 | 0.0 | 0.05
Modify | 0.041394 | 0.041394 | 0.041394 | 0.0 | 99.65
Other | | 9.322e-05 | | | 0.22
Nlocal: 5 ave 5 max 5 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 7 ave 7 max 7 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 10 ave 10 max 10 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 10
Ave neighs/atom = 2
Ave special neighs/atom = 0
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:01

1
lib/.gitignore vendored
View File

@ -1,3 +1,4 @@
Makefile.lammps
.depend
Makefile.auto
/latte

5
lib/latte/.gitignore vendored
View File

@ -1,5 +0,0 @@
# files and folders to ignore
/filelink
/liblink
/includelink
/LATTE-*

146
lib/latte/Install.py
View File

@ -1,146 +0,0 @@
#!/usr/bin/env python
"""
Install.py tool to download, unpack, build, and link to the LATTE library
used to automate the steps described in the README file in this dir
"""
from __future__ import print_function
import sys, os, subprocess, shutil, tarfile
from argparse import ArgumentParser
sys.path.append('..')
from install_helpers import fullpath, geturl, checkmd5sum, getfallback
parser = ArgumentParser(prog='Install.py',
description="LAMMPS library build wrapper script")
# settings
version = '1.2.2'
suffix = 'gfortran'
# known checksums for different LATTE versions. used to validate the download.
checksums = { \
'1.1.0' : '533635721ee222d0ed2925a18fb5b294', \
'1.2.0' : '68bf0db879da5e068a71281020239ae7', \
'1.2.1' : '85ac414fdada2d04619c8f936344df14', \
'1.2.2' : '820e73a457ced178c08c71389a385de7', \
}
# help message
HELP = """
Syntax from src dir: make lib-latte args="-b"
or: make lib-latte args="-p /usr/local/latte"
or: make lib-latte args="-m gfortran"
or: make lib-latte args="-b -v 1.2.2"
Syntax from lib dir: python Install.py -b
or: python Install.py -p /usr/local/latte
or: python Install.py -m gfortran
or: python Install.py -v 1.2.2 -b
Example:
make lib-latte args="-b -m gfortran" # download/build in lib/latte
make lib-latte args="-p $HOME/latte" # use existing LATTE installation
"""
pgroup = parser.add_mutually_exclusive_group()
pgroup.add_argument("-b", "--build", action="store_true",
help="download and build the LATTE library")
pgroup.add_argument("-p", "--path",
help="specify folder of existing LATTE installation")
parser.add_argument("-m", "--machine", choices=['gfortran', 'ifort', 'linalg', 'serial', 'mpi'],
help="suffix of a Makefile.lammps.* file used for linking LAMMPS with this library")
parser.add_argument("-v", "--version", default=version,
help="set version of LATTE to download and build (default: %s)" % version)
args = parser.parse_args()
# print help message and exit, if neither build nor path options are given
if not args.build and not args.path:
parser.print_help()
sys.exit(HELP)
homepath = fullpath(".")
buildflag = args.build
pathflag = args.path is not None
version = args.version
suffixflag = args.machine is not None
if suffixflag:
suffix = args.machine
if pathflag:
lattedir = args.path
if not os.path.isdir(lattedir):
sys.exit("LATTE path %s does not exist" % lattedir)
lattedir = fullpath(lattedir)
homedir = "LATTE-%s" % version
if buildflag:
url = "https://github.com/lanl/LATTE/archive/v%s.tar.gz" % version
lattepath = fullpath(homepath)
lattedir = os.path.join(lattepath, homedir)
fallback = getfallback('latte', url)
filename = 'LATTE.tar.gz'
# download and unpack LATTE tarball
if buildflag:
print("Downloading LATTE ...")
try:
geturl(url, filename)
except:
geturl(fallback, filename)
# verify downloaded archive integrity via md5 checksum, if known.
if version in checksums:
if not checkmd5sum(checksums[version], filename):
print("Checksum did not match. Trying fallback URL", fallback)
geturl(fallback, filename)
if not checkmd5sum(checksums[version], filename):
sys.exit("Checksum for LATTE library does not match for fallback, too.")
print("Unpacking LATTE ...")
if os.path.exists(lattedir):
shutil.rmtree(lattedir)
if tarfile.is_tarfile('LATTE.tar.gz'):
tgz = tarfile.open('LATTE.tar.gz')
tgz.extractall()
os.remove('LATTE.tar.gz')
else:
sys.exit("File LATTE.tar.gz is not a supported archive")
# build LATTE
print("Building LATTE ...")
cmd = 'cd "%s"; make' % lattedir
try:
txt = subprocess.check_output(cmd, stderr=subprocess.STDOUT, shell=True)
print(txt.decode('UTF-8'))
except subprocess.CalledProcessError as e:
sys.exit("Make failed with:\n %s" % e.output.decode('UTF-8'))
# create 3 links in lib/latte to LATTE dirs
print("Creating links to LATTE files")
if os.path.isfile("includelink") or os.path.islink("includelink"):
os.remove("includelink")
if os.path.isfile("liblink") or os.path.islink("liblink"):
os.remove("liblink")
if os.path.isfile("filelink.o") or os.path.islink("filelink.o"):
os.remove("filelink.o")
os.symlink(os.path.join(lattedir, 'src'), 'includelink')
os.symlink(lattedir, 'liblink')
os.symlink(os.path.join(lattedir, 'src', 'latte_c_bind.o'), 'filelink.o')
# copy Makefile.lammps.suffix to Makefile.lammps
if suffixflag or not os.path.exists("Makefile.lammps"):
print("Creating Makefile.lammps")
if os.path.exists("Makefile.lammps.%s" % suffix):
shutil.copyfile("Makefile.lammps.%s" % suffix, 'Makefile.lammps')

7
lib/latte/Makefile.lammps.gfortran
View File

@ -1,7 +0,0 @@
# Settings that the LAMMPS build will import when this package library is used
# GNU Fortran settings
latte_SYSINC =
latte_SYSLIB = ../../lib/latte/filelink.o -llatte -lgfortran -llapack -lblas
latte_SYSPATH = -fopenmp

12
lib/latte/Makefile.lammps.ifort
View File

@ -1,12 +0,0 @@
# Settings that the LAMMPS build will import when this package library is used
# Intel ifort settings
latte_SYSINC =
latte_SYSLIB = ../../lib/latte/filelink.o \
-llatte -lifport -lifcore -lsvml -lompstub -limf \
-lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core \
-lmkl_intel_thread -lpthread -openmp
latte_SYSPATH = -openmp -L${MKLROOT}/lib/intel64 -lmkl_lapack95_lp64 \
-L/opt/intel/composer_xe_2013_sp1.2.144/compiler/lib/intel64

7
lib/latte/Makefile.lammps.linalg
View File

@ -1,7 +0,0 @@
# Settings that the LAMMPS build will import when this package library is used
# GNU Fortran settings for use with bundled linalg lib
latte_SYSINC =
latte_SYSLIB = ../../lib/latte/filelink.o -llatte -llinalg -lgfortran
latte_SYSPATH = -L../../lib/linalg -fopenmp

1
lib/latte/Makefile.lammps.mpi
View File

@ -1 +0,0 @@
Makefile.lammps.linalg

1
lib/latte/Makefile.lammps.serial
View File

@ -1 +0,0 @@
Makefile.lammps.linalg

56
lib/latte/README
View File

@ -1,56 +0,0 @@
This directory contains links to the LATTE library which is required
to use the LATTE package and its fix latte command in a LAMMPS input
script.
Information about the LATTE DFTB code can be found at:
https://github.com/lanl/LATTE
The LATTE development effort is led by Marc Cawkwell and
Anders Niklasson at Los Alamos National Laboratory.
You can type "make lib-latte" from the src directory to see help on
how to download and build this library via make commands, or you can
do the same thing by typing "python Install.py" from within this
directory, or you can do it manually by following the instructions
below.
-----------------
Instructions:
1. Download or clone the LATTE source code from
https://github.com/lanl/LATTE. If you download a zipfile
or tarball, unpack the tarball either in this /lib/latte
directory or somewhere else on your system.
2. Modify the makefile.CHOICES according to your system architecture
and compilers. Check that the MAKELIB flag is ON in makefile.CHOICES
and finally, build the code via the make command
% make
3. Create three symbolic links in this dir (lib/latte)
E.g if you built LATTE in this dir:
% ln -s ./LATTE-master/src includelink
% ln -s ./LATTE-master liblink
% ln -s ./LATTE-master/src/latte_c_bind.o filelink.o
4. Choose a Makefile.lammps.* file appropriate for your compiler
(GNU gfortran with external BLAS, GNU gfortran with local liblinalg,
or Intel ifort with MKL) and copy it to Makefile.lammps.
Note that you may need to edit Makefile.lammps for paths
and compiler options appropriate to your system.
-----------------
When these steps are complete you can build LAMMPS
with the LATTE package installed:
% cd lammps/src
% make yes-latte
% make g++ (or whatever target you wish)
Note that if you download and unpack a new LAMMPS tarball, the
"includelink" and "liblink" and "filelink.o" symbolic links will be
lost and you will need to re-create them (step 3). If you built LATTE
in this directory (as opposed to somewhere else on your system), you
will also need to repeat steps 1,2,4.

2
src/.gitignore vendored
View File

@ -814,8 +814,6 @@
/fix_lambdah_calc.h
/fix_langevin_eff.cpp
/fix_langevin_eff.h
/fix_latte.cpp
/fix_latte.h
/latboltz_const.h
/fix_lb_fluid.cpp
/fix_lb_fluid.h

2
src/ASPHERE/pair_ylz.h
View File

@ -31,7 +31,7 @@ class PairYLZ : public Pair {
~PairYLZ() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **);
void coeff(int, char **) override;
void init_style() override;
double init_one(int, int) override;
void write_restart(FILE *) override;

15
src/COLVARS/fix_colvars.cpp
View File

@ -324,7 +324,7 @@ FixColvars::FixColvars(LAMMPS *lmp, int narg, char **arg) :
if (!out_name) out_name = utils::strdup("out");
/* initialize various state variables. */
tstat_id = -1;
tstat_fix = nullptr;
energy = 0.0;
nlevels_respa = 0;
init_flag = 0;
@ -432,12 +432,13 @@ void FixColvars::one_time_init()
double t_target = 0.0;
if (tmp_name) {
if (strcmp(tmp_name,"NULL") == 0) {
tstat_id = -1;
tstat_fix = nullptr;
} else {
tstat_id = modify->find_fix(tmp_name);
if (tstat_id < 0) error->one(FLERR,"Could not find tstat fix ID");
double *tt = (double*)modify->fix[tstat_id]->extract("t_target",tmp);
tstat_fix = modify->get_fix_by_id(tmp_name);
if (!tstat_fix) error->one(FLERR, "Could not find thermostat fix ID {}", tmp_name);
double *tt = (double*) tstat_fix->extract("t_target", tmp);
if (tt) t_target = *tt;
else error->one(FLERR, "Fix ID {} is not a thermostat fix", tmp_name);
}
}
@ -675,13 +676,13 @@ void FixColvars::post_force(int /*vflag*/)
if (proxy->want_exit())
error->one(FLERR,"Run aborted on request from colvars module.\n");
if (tstat_id < 0) {
if (!tstat_fix) {
proxy->set_temperature(0.0);
} else {
int tmp;
// get thermostat target temperature from corresponding fix,
// if the fix supports extraction.
double *tt = (double *) modify->fix[tstat_id]->extract("t_target",tmp);
double *tt = (double *) tstat_fix->extract("t_target", tmp);
if (tt)
proxy->set_temperature(*tt);
else

2
src/COLVARS/fix_colvars.h
View File

@ -68,7 +68,7 @@ class FixColvars : public Fix {
char *out_name; // prefix string for all output files
char *tmp_name; // name of thermostat fix.
int rng_seed; // seed to initialize random number generator
int tstat_id; // id of the thermostat fix
Fix *tstat_fix; // pointer to thermostat fix
double energy; // biasing energy of the fix
int me; // my MPI rank in this "world".

10
src/EXTRA-COMPUTE/compute_hma.cpp
View File

@ -188,11 +188,11 @@ void ComputeHMA::init_list(int /* id */, NeighList *ptr)
void ComputeHMA::setup()
{
int dummy=0;
int ifix = modify->find_fix(id_temp);
if (ifix < 0) error->all(FLERR,"Could not find compute hma temperature ID");
auto temperat = (double *) modify->fix[ifix]->extract("t_target",dummy);
if (temperat==nullptr) error->all(FLERR,"Could not find compute hma temperature ID");
finaltemp = * temperat;
Fix *ifix = modify->get_fix_by_id(id_temp);
if (!ifix) error->all(FLERR,"Could not find compute hma temperature fix ID {}", id_temp);
auto temperat = (double *) ifix->extract("t_target",dummy);
if (temperat == nullptr) error->all(FLERR,"Fix ID {} is not a thermostat {}", id_temp);
finaltemp = *temperat;
// set fix which stores original atom coords

39
src/EXTRA-FIX/fix_controller.cpp
View File

@ -71,28 +71,22 @@ FixController::FixController(LAMMPS *lmp, int narg, char **arg) :
// error check
if (pvwhich == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(pvID);
if (icompute < 0)
error->all(FLERR,"Compute ID for fix controller does not exist");
Compute *c = modify->compute[icompute];
Compute *c = modify->get_compute_by_id(pvID);
if (!c) error->all(FLERR,"Compute ID {} for fix controller does not exist", pvID);
int flag = 0;
if (c->scalar_flag && pvindex == 0) flag = 1;
else if (c->vector_flag && pvindex > 0) flag = 1;
if (!flag) error->all(FLERR,"Fix controller compute does not "
"calculate a global scalar or vector");
if (!flag)
error->all(FLERR,"Fix controller compute does not calculate a global scalar or vector");
if (pvindex && pvindex > c->size_vector)
error->all(FLERR,"Fix controller compute vector is "
"accessed out-of-range");
error->all(FLERR,"Fix controller compute vector is accessed out-of-range");
} else if (pvwhich == ArgInfo::FIX) {
int ifix = modify->find_fix(pvID);
if (ifix < 0)
error->all(FLERR,"Fix ID for fix controller does not exist");
Fix *f = modify->fix[ifix];
Fix *f = modify->get_fix_by_id(pvID);
if (!f) error->all(FLERR,"Fix ID {} for fix controller does not exist", pvID);
int flag = 0;
if (f->scalar_flag && pvindex == 0) flag = 1;
else if (f->vector_flag && pvindex > 0) flag = 1;
if (!flag) error->all(FLERR,"Fix controller fix does not "
"calculate a global scalar or vector");
if (!flag) error->all(FLERR,"Fix controller fix does not calculate a global scalar or vector");
if (pvindex && pvindex > f->size_vector)
error->all(FLERR,"Fix controller fix vector is accessed out-of-range");
} else if (pvwhich == ArgInfo::VARIABLE) {
@ -135,25 +129,20 @@ int FixController::setmask()
void FixController::init()
{
if (pvwhich == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(pvID);
if (icompute < 0)
error->all(FLERR,"Compute ID for fix controller does not exist");
pcompute = modify->compute[icompute];
pcompute = modify->get_compute_by_id(pvID);
if (!pcompute) error->all(FLERR,"Compute ID {} for fix controller does not exist", pvID);
} else if (pvwhich == ArgInfo::FIX) {
int ifix = modify->find_fix(pvID);
if (ifix < 0) error->all(FLERR,"Fix ID for fix controller does not exist");
pfix = modify->fix[ifix];
pfix = modify->get_fix_by_id(pvID);
if (!pfix) error->all(FLERR,"Fix ID {} for fix controller does not exist", pvID);
} else if (pvwhich == ArgInfo::VARIABLE) {
pvar = input->variable->find(pvID);
if (pvar < 0)
error->all(FLERR,"Variable name for fix controller does not exist");
if (pvar < 0) error->all(FLERR,"Variable name for fix controller does not exist");
}
cvar = input->variable->find(cvID);
if (cvar < 0)
error->all(FLERR,"Variable name for fix controller does not exist");
if (cvar < 0) error->all(FLERR,"Variable name for fix controller does not exist");
// set sampling time

4
src/EXTRA-PAIR/pair_born_gauss.cpp
View File

@ -31,7 +31,9 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairBornGauss::PairBornGauss(LAMMPS *lmp) : Pair(lmp)
PairBornGauss::PairBornGauss(LAMMPS *lmp) :
Pair(lmp), cut(nullptr), biga0(nullptr), alpha(nullptr), biga1(nullptr), beta(nullptr),
r0(nullptr), offset(nullptr)
{
single_enable = 1;
respa_enable = 0;

526
src/KOKKOS/fix_shake_kokkos.cpp
View File

@ -144,6 +144,7 @@ FixShakeKokkos<DeviceType>::~FixShakeKokkos()
memoryKK->destroy_kokkos(k_shake_type,shake_type);
memoryKK->destroy_kokkos(k_xshake,xshake);
memoryKK->destroy_kokkos(k_list,list);
memoryKK->destroy_kokkos(k_closest_list,closest_list);
memoryKK->destroy_kokkos(k_vatom,vatom);
}
@ -204,13 +205,10 @@ void FixShakeKokkos<DeviceType>::pre_neighbor()
// local copies of atom quantities
// used by SHAKE until next re-neighboring
atomKK->sync(execution_space,X_MASK);
ebond = 0.0;
x = atom->x;
v = atom->v;
f = atom->f;
mass = atom->mass;
rmass = atom->rmass;
type = atom->type;
d_x = atomKK->k_x.view<DeviceType>();
nlocal = atom->nlocal;
map_style = atom->map_style;
@ -231,6 +229,10 @@ void FixShakeKokkos<DeviceType>::pre_neighbor()
memoryKK->destroy_kokkos(k_list,list);
memoryKK->create_kokkos(k_list,list,maxlist,"shake:list");
d_list = k_list.view<DeviceType>();
memoryKK->destroy_kokkos(k_closest_list,closest_list);
memoryKK->create_kokkos(k_closest_list,closest_list,maxlist,4,"shake:closest_list");
d_closest_list = k_closest_list.view<DeviceType>();
}
// Atom Map
@ -244,60 +246,22 @@ void FixShakeKokkos<DeviceType>::pre_neighbor()
k_map_hash = atomKK->k_map_hash;
}
k_sametag = atomKK->k_sametag;
k_sametag.template sync<DeviceType>();
d_sametag = k_sametag.view<DeviceType>();
// build list of SHAKE clusters I compute
Kokkos::deep_copy(d_scalars,0);
{
// local variables for lambda capture
copymode = 1;
auto d_shake_flag = this->d_shake_flag;
auto d_shake_atom = this->d_shake_atom;
auto d_list = this->d_list;
auto d_error_flag = this->d_error_flag;
auto d_nlist = this->d_nlist;
auto map_style = atom->map_style;
auto k_map_array = this->k_map_array;
auto k_map_hash = this->k_map_hash;
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagFixShakePreNeighbor>(0,nlocal),*this);
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType>(0,nlocal),
LAMMPS_LAMBDA(const int& i) {
if (d_shake_flag[i]) {
if (d_shake_flag[i] == 2) {
const int atom1 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,0),map_style,k_map_array,k_map_hash);
const int atom2 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,1),map_style,k_map_array,k_map_hash);
if (atom1 == -1 || atom2 == -1) {
d_error_flag() = 1;
}
if (i <= atom1 && i <= atom2) {
const int nlist = Kokkos::atomic_fetch_add(&d_nlist(),1);
d_list[nlist] = i;
}
} else if (d_shake_flag[i] % 2 == 1) {
const int atom1 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,0),map_style,k_map_array,k_map_hash);
const int atom2 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,1),map_style,k_map_array,k_map_hash);
const int atom3 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,2),map_style,k_map_array,k_map_hash);
if (atom1 == -1 || atom2 == -1 || atom3 == -1)
d_error_flag() = 1;
if (i <= atom1 && i <= atom2 && i <= atom3) {
const int nlist = Kokkos::atomic_fetch_add(&d_nlist(),1);
d_list[nlist] = i;
}
} else {
const int atom1 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,0),map_style,k_map_array,k_map_hash);
const int atom2 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,1),map_style,k_map_array,k_map_hash);
const int atom3 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,2),map_style,k_map_array,k_map_hash);
const int atom4 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,3),map_style,k_map_array,k_map_hash);
if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1)
d_error_flag() = 1;
if (i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4) {
const int nlist = Kokkos::atomic_fetch_add(&d_nlist(),1);
d_list[nlist] = i;
}
}
}
});
}
copymode = 0;
k_list.modify<DeviceType>();
k_closest_list.modify<DeviceType>();
Kokkos::deep_copy(h_scalars,d_scalars);
nlist = h_nlist();
@ -308,6 +272,65 @@ void FixShakeKokkos<DeviceType>::pre_neighbor()
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixShakeKokkos<DeviceType>::operator()(TagFixShakePreNeighbor, const int &i) const {
if (d_shake_flag[i]) {
if (d_shake_flag[i] == 2) {
int atom1 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,0),map_style,k_map_array,k_map_hash);
int atom2 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,1),map_style,k_map_array,k_map_hash);
if (atom1 == -1 || atom2 == -1) {
d_error_flag() = 1;
}
atom1 = closest_image(i, atom1);
atom2 = closest_image(i, atom2);
if (i <= atom1 && i <= atom2) {
const int nlist = Kokkos::atomic_fetch_add(&d_nlist(),1);
d_list[nlist] = i;
d_closest_list(nlist,0) = atom1;
d_closest_list(nlist,1) = atom2;
}
} else if (d_shake_flag[i] % 2 == 1) {
int atom1 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,0),map_style,k_map_array,k_map_hash);
int atom2 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,1),map_style,k_map_array,k_map_hash);
int atom3 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,2),map_style,k_map_array,k_map_hash);
if (atom1 == -1 || atom2 == -1 || atom3 == -1)
d_error_flag() = 1;
atom1 = closest_image(i, atom1);
atom2 = closest_image(i, atom2);
atom3 = closest_image(i, atom3);
if (i <= atom1 && i <= atom2 && i <= atom3) {
const int nlist = Kokkos::atomic_fetch_add(&d_nlist(),1);
d_list[nlist] = i;
d_closest_list(nlist,0) = atom1;
d_closest_list(nlist,1) = atom2;
d_closest_list(nlist,2) = atom3;
}
} else {
int atom1 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,0),map_style,k_map_array,k_map_hash);
int atom2 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,1),map_style,k_map_array,k_map_hash);
int atom3 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,2),map_style,k_map_array,k_map_hash);
int atom4 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(i,3),map_style,k_map_array,k_map_hash);
if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1)
d_error_flag() = 1;
atom1 = closest_image(i, atom1);
atom2 = closest_image(i, atom2);
atom3 = closest_image(i, atom3);
atom4 = closest_image(i, atom4);
if (i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4) {
const int nlist = Kokkos::atomic_fetch_add(&d_nlist(),1);
d_list[nlist] = i;
d_closest_list(nlist,0) = atom1;
d_closest_list(nlist,1) = atom2;
d_closest_list(nlist,2) = atom3;
d_closest_list(nlist,3) = atom4;
}
}
}
}
/* ----------------------------------------------------------------------
compute the force adjustment for SHAKE constraint
------------------------------------------------------------------------- */
@ -339,22 +362,18 @@ void FixShakeKokkos<DeviceType>::post_force(int vflag)
atomKK->sync(execution_space,X_MASK|F_MASK|TYPE_MASK);
k_shake_flag.sync<DeviceType>();
k_shake_atom.sync<DeviceType>();
k_shake_type.sync<DeviceType>();
k_list.sync<DeviceType>();
k_closest_list.sync<DeviceType>();
if (update->ntimestep == next_output) {
atomKK->sync(Host,X_MASK);
k_shake_flag.sync_host();
k_shake_atom.sync_host();
k_shake_type.sync_host();
if (update->ntimestep == next_output)
stats();
}
// xshake = unconstrained move with current v,f
// communicate results if necessary
unconstrained_update();
if (comm->nprocs > 1) comm->forward_comm(this);
comm->forward_comm(this);
k_xshake.sync<DeviceType>();
// virial setup
@ -369,7 +388,6 @@ void FixShakeKokkos<DeviceType>::post_force(int vflag)
d_vatom = k_vatom.template view<KKDeviceType>();
}
neighflag = lmp->kokkos->neighflag;
// FULL neighlist still needs atomics in fix shake
@ -397,8 +415,6 @@ void FixShakeKokkos<DeviceType>::post_force(int vflag)
Kokkos::deep_copy(d_error_flag,0);
update_domain_variables();
EV_FLOAT ev;
// loop over clusters to add constraint forces
@ -457,6 +473,26 @@ void FixShakeKokkos<DeviceType>::post_force(int vflag)
}
}
/* ----------------------------------------------------------------------
substitute shake constraints with very strong bonds
------------------------------------------------------------------------- */
template<class DeviceType>
void FixShakeKokkos<DeviceType>::min_post_force(int vflag)
{
// not yet ported to Kokkos
atomKK->sync(Host,X_MASK | F_MASK);
k_shake_flag.sync_host();
k_shake_type.sync_host();
k_list.sync_host();
k_closest_list.sync_host();
FixShake::min_post_force(vflag);
atomKK->modified(Host,F_MASK);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
@ -464,10 +500,10 @@ template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void FixShakeKokkos<DeviceType>::operator()(TagFixShakePostForce<NEIGHFLAG,EVFLAG>, const int &i, EV_FLOAT& ev) const {
const int m = d_list[i];
if (d_shake_flag[m] == 2) shake<NEIGHFLAG,EVFLAG>(m,ev);
else if (d_shake_flag[m] == 3) shake3<NEIGHFLAG,EVFLAG>(m,ev);
else if (d_shake_flag[m] == 4) shake4<NEIGHFLAG,EVFLAG>(m,ev);
else shake3angle<NEIGHFLAG,EVFLAG>(m,ev);
if (d_shake_flag[m] == 2) shake<NEIGHFLAG,EVFLAG>(i,ev);
else if (d_shake_flag[m] == 3) shake3<NEIGHFLAG,EVFLAG>(i,ev);
else if (d_shake_flag[m] == 4) shake4<NEIGHFLAG,EVFLAG>(i,ev);
else shake3angle<NEIGHFLAG,EVFLAG>(i,ev);
}
template<class DeviceType>
@ -485,7 +521,6 @@ void FixShakeKokkos<DeviceType>::operator()(TagFixShakePostForce<NEIGHFLAG,EVFLA
template<class DeviceType>
int FixShakeKokkos<DeviceType>::dof(int igroup)
{
d_mask = atomKK->k_mask.view<DeviceType>();
d_tag = atomKK->k_tag.view<DeviceType>();
nlocal = atom->nlocal;
@ -545,7 +580,6 @@ void FixShakeKokkos<DeviceType>::unconstrained_update()
else
atomKK->sync(execution_space,X_MASK|V_MASK|F_MASK|TYPE_MASK);
k_shake_flag.sync<DeviceType>();
k_xshake.sync<DeviceType>();
@ -593,12 +627,14 @@ void FixShakeKokkos<DeviceType>::unconstrained_update()
k_xshake.modify<DeviceType>();
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
calculate SHAKE constraint forces for size 2 cluster = single bond
------------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void FixShakeKokkos<DeviceType>::shake(int m, EV_FLOAT& ev) const
void FixShakeKokkos<DeviceType>::shake(int ilist, EV_FLOAT& ev) const
{
// The f array is duplicated for OpenMP, atomic for CUDA, and neither for Serial
@ -606,33 +642,30 @@ void FixShakeKokkos<DeviceType>::shake(int m, EV_FLOAT& ev) const
auto v_f = ScatterViewHelper<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_f),decltype(ndup_f)>::get(dup_f,ndup_f);
auto a_f = v_f.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
int nlist,list[2];
int atomlist[2];
double v[6];
double invmass0,invmass1;
// local atom IDs and constraint distances
int i0 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(m,0),map_style,k_map_array,k_map_hash);
int i1 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(m,1),map_style,k_map_array,k_map_hash);
int m = d_list[ilist];
int i0 = d_closest_list(ilist,0);
int i1 = d_closest_list(ilist,1);
double bond1 = d_bond_distance[d_shake_type(m,0)];
// r01 = distance vec between atoms, with PBC
// r01 = distance vec between atoms
double r01[3];
r01[0] = d_x(i0,0) - d_x(i1,0);
r01[1] = d_x(i0,1) - d_x(i1,1);
r01[2] = d_x(i0,2) - d_x(i1,2);
minimum_image(r01);
// s01 = distance vec after unconstrained update, with PBC
// use Domain::minimum_image_once(), not minimum_image()
// b/c xshake values might be huge, due to e.g. fix gcmc
// s01 = distance vec after unconstrained update
double s01[3];
s01[0] = d_xshake(i0,0) - d_xshake(i1,0);
s01[1] = d_xshake(i0,1) - d_xshake(i1,1);
s01[2] = d_xshake(i0,2) - d_xshake(i1,2);
minimum_image_once(s01);
// scalar distances between atoms
@ -689,9 +722,9 @@ void FixShakeKokkos<DeviceType>::shake(int m, EV_FLOAT& ev) const
}
if (EVFLAG) {
nlist = 0;
if (i0 < nlocal) list[nlist++] = i0;
if (i1 < nlocal) list[nlist++] = i1;
int count = 0;
if (i0 < nlocal) atomlist[count++] = i0;
if (i1 < nlocal) atomlist[count++] = i1;
v[0] = lamda*r01[0]*r01[0];
v[1] = lamda*r01[1]*r01[1];
@ -700,16 +733,18 @@ void FixShakeKokkos<DeviceType>::shake(int m, EV_FLOAT& ev) const
v[4] = lamda*r01[0]*r01[2];
v[5] = lamda*r01[1]*r01[2];
v_tally<NEIGHFLAG>(ev,nlist,list,2.0,v);
v_tally<NEIGHFLAG>(ev,count,atomlist,2.0,v);
}
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
calculate SHAKE constraint forces for size 3 cluster = two bonds
------------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void FixShakeKokkos<DeviceType>::shake3(int m, EV_FLOAT& ev) const
void FixShakeKokkos<DeviceType>::shake3(int ilist, EV_FLOAT& ev) const
{
// The f array is duplicated for OpenMP, atomic for CUDA, and neither for Serial
@ -717,47 +752,42 @@ void FixShakeKokkos<DeviceType>::shake3(int m, EV_FLOAT& ev) const
auto v_f = ScatterViewHelper<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_f),decltype(ndup_f)>::get(dup_f,ndup_f);
auto a_f = v_f.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
int nlist,list[3];
int atomlist[3];
double v[6];
double invmass0,invmass1,invmass2;
// local atom IDs and constraint distances
int i0 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(m,0),map_style,k_map_array,k_map_hash);
int i1 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(m,1),map_style,k_map_array,k_map_hash);
int i2 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(m,2),map_style,k_map_array,k_map_hash);
int m = d_list[ilist];
int i0 = d_closest_list(ilist,0);
int i1 = d_closest_list(ilist,1);
int i2 = d_closest_list(ilist,2);
double bond1 = d_bond_distance[d_shake_type(m,0)];
double bond2 = d_bond_distance[d_shake_type(m,1)];
// r01,r02 = distance vec between atoms, with PBC
// r01,r02 = distance vec between atoms
double r01[3];
r01[0] = d_x(i0,0) - d_x(i1,0);
r01[1] = d_x(i0,1) - d_x(i1,1);
r01[2] = d_x(i0,2) - d_x(i1,2);
minimum_image(r01);
double r02[3];
r02[0] = d_x(i0,0) - d_x(i2,0);
r02[1] = d_x(i0,1) - d_x(i2,1);
r02[2] = d_x(i0,2) - d_x(i2,2);
minimum_image(r02);
// s01,s02 = distance vec after unconstrained update, with PBC
// use Domain::minimum_image_once(), not minimum_image()
// b/c xshake values might be huge, due to e.g. fix gcmc
// s01,s02 = distance vec after unconstrained update
double s01[3];
s01[0] = d_xshake(i0,0) - d_xshake(i1,0);
s01[1] = d_xshake(i0,1) - d_xshake(i1,1);
s01[2] = d_xshake(i0,2) - d_xshake(i1,2);
minimum_image_once(s01);
double s02[3];
s02[0] = d_xshake(i0,0) - d_xshake(i2,0);
s02[1] = d_xshake(i0,1) - d_xshake(i2,1);
s02[2] = d_xshake(i0,2) - d_xshake(i2,2);
minimum_image_once(s02);
// scalar distances between atoms
@ -871,10 +901,10 @@ void FixShakeKokkos<DeviceType>::shake3(int m, EV_FLOAT& ev) const
}
if (EVFLAG) {
nlist = 0;
if (i0 < nlocal) list[nlist++] = i0;
if (i1 < nlocal) list[nlist++] = i1;
if (i2 < nlocal) list[nlist++] = i2;
int count = 0;
if (i0 < nlocal) atomlist[count++] = i0;
if (i1 < nlocal) atomlist[count++] = i1;
if (i2 < nlocal) atomlist[count++] = i2;
v[0] = lamda01*r01[0]*r01[0] + lamda02*r02[0]*r02[0];
v[1] = lamda01*r01[1]*r01[1] + lamda02*r02[1]*r02[1];
@ -883,16 +913,18 @@ void FixShakeKokkos<DeviceType>::shake3(int m, EV_FLOAT& ev) const
v[4] = lamda01*r01[0]*r01[2] + lamda02*r02[0]*r02[2];
v[5] = lamda01*r01[1]*r01[2] + lamda02*r02[1]*r02[2];
v_tally<NEIGHFLAG>(ev,nlist,list,3.0,v);
v_tally<NEIGHFLAG>(ev,count,atomlist,3.0,v);
}
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
calculate SHAKE constraint forces for size 4 cluster = three bonds
------------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void FixShakeKokkos<DeviceType>::shake4(int m, EV_FLOAT& ev) const
void FixShakeKokkos<DeviceType>::shake4(int ilist, EV_FLOAT& ev) const
{
// The f array is duplicated for OpenMP, atomic for CUDA, and neither for Serial
@ -900,61 +932,54 @@ void FixShakeKokkos<DeviceType>::shake4(int m, EV_FLOAT& ev) const
auto v_f = ScatterViewHelper<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_f),decltype(ndup_f)>::get(dup_f,ndup_f);
auto a_f = v_f.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
int nlist,list[4];
int atomlist[4];
double v[6];
double invmass0,invmass1,invmass2,invmass3;
// local atom IDs and constraint distances
int i0 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(m,0),map_style,k_map_array,k_map_hash);
int i1 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(m,1),map_style,k_map_array,k_map_hash);
int i2 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(m,2),map_style,k_map_array,k_map_hash);
int i3 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(m,3),map_style,k_map_array,k_map_hash);
int m = d_list[ilist];
int i0 = d_closest_list(ilist,0);
int i1 = d_closest_list(ilist,1);
int i2 = d_closest_list(ilist,2);
int i3 = d_closest_list(ilist,3);
double bond1 = d_bond_distance[d_shake_type(m,0)];
double bond2 = d_bond_distance[d_shake_type(m,1)];
double bond3 = d_bond_distance[d_shake_type(m,2)];
// r01,r02,r03 = distance vec between atoms, with PBC
// r01,r02,r03 = distance vec between atoms
double r01[3];
r01[0] = d_x(i0,0) - d_x(i1,0);
r01[1] = d_x(i0,1) - d_x(i1,1);
r01[2] = d_x(i0,2) - d_x(i1,2);
minimum_image(r01);
double r02[3];
r02[0] = d_x(i0,0) - d_x(i2,0);
r02[1] = d_x(i0,1) - d_x(i2,1);
r02[2] = d_x(i0,2) - d_x(i2,2);
minimum_image(r02);
double r03[3];
r03[0] = d_x(i0,0) - d_x(i3,0);
r03[1] = d_x(i0,1) - d_x(i3,1);
r03[2] = d_x(i0,2) - d_x(i3,2);
minimum_image(r03);
// s01,s02,s03 = distance vec after unconstrained update, with PBC
// use Domain::minimum_image_once(), not minimum_image()
// b/c xshake values might be huge, due to e.g. fix gcmc
// s01,s02,s03 = distance vec after unconstrained update
double s01[3];
s01[0] = d_xshake(i0,0) - d_xshake(i1,0);
s01[1] = d_xshake(i0,1) - d_xshake(i1,1);
s01[2] = d_xshake(i0,2) - d_xshake(i1,2);
minimum_image_once(s01);
double s02[3];
s02[0] = d_xshake(i0,0) - d_xshake(i2,0);
s02[1] = d_xshake(i0,1) - d_xshake(i2,1);
s02[2] = d_xshake(i0,2) - d_xshake(i2,2);
minimum_image_once(s02);
double s03[3];
s03[0] = d_xshake(i0,0) - d_xshake(i3,0);
s03[1] = d_xshake(i0,1) - d_xshake(i3,1);
s03[2] = d_xshake(i0,2) - d_xshake(i3,2);
minimum_image_once(s03);
// scalar distances between atoms
@ -1132,11 +1157,11 @@ void FixShakeKokkos<DeviceType>::shake4(int m, EV_FLOAT& ev) const
}
if (EVFLAG) {
nlist = 0;
if (i0 < nlocal) list[nlist++] = i0;
if (i1 < nlocal) list[nlist++] = i1;
if (i2 < nlocal) list[nlist++] = i2;
if (i3 < nlocal) list[nlist++] = i3;
int count = 0;
if (i0 < nlocal) atomlist[count++] = i0;
if (i1 < nlocal) atomlist[count++] = i1;
if (i2 < nlocal) atomlist[count++] = i2;
if (i3 < nlocal) atomlist[count++] = i3;
v[0] = lamda01*r01[0]*r01[0]+lamda02*r02[0]*r02[0]+lamda03*r03[0]*r03[0];
v[1] = lamda01*r01[1]*r01[1]+lamda02*r02[1]*r02[1]+lamda03*r03[1]*r03[1];
@ -1145,16 +1170,18 @@ void FixShakeKokkos<DeviceType>::shake4(int m, EV_FLOAT& ev) const
v[4] = lamda01*r01[0]*r01[2]+lamda02*r02[0]*r02[2]+lamda03*r03[0]*r03[2];
v[5] = lamda01*r01[1]*r01[2]+lamda02*r02[1]*r02[2]+lamda03*r03[1]*r03[2];
v_tally<NEIGHFLAG>(ev,nlist,list,4.0,v);
v_tally<NEIGHFLAG>(ev,count,atomlist,4.0,v);
}
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
calculate SHAKE constraint forces for size 3 cluster = two bonds + angle
------------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void FixShakeKokkos<DeviceType>::shake3angle(int m, EV_FLOAT& ev) const
void FixShakeKokkos<DeviceType>::shake3angle(int ilist, EV_FLOAT& ev) const
{
// The f array is duplicated for OpenMP, atomic for CUDA, and neither for Serial
@ -1162,60 +1189,53 @@ void FixShakeKokkos<DeviceType>::shake3angle(int m, EV_FLOAT& ev) const
auto v_f = ScatterViewHelper<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_f),decltype(ndup_f)>::get(dup_f,ndup_f);
auto a_f = v_f.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
int nlist,list[3];
int atomlist[3];
double v[6];
double invmass0,invmass1,invmass2;
// local atom IDs and constraint distances
int i0 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(m,0),map_style,k_map_array,k_map_hash);
int i1 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(m,1),map_style,k_map_array,k_map_hash);
int i2 = AtomKokkos::map_kokkos<DeviceType>(d_shake_atom(m,2),map_style,k_map_array,k_map_hash);
int m = d_list[ilist];
int i0 = d_closest_list(ilist,0);
int i1 = d_closest_list(ilist,1);
int i2 = d_closest_list(ilist,2);
double bond1 = d_bond_distance[d_shake_type(m,0)];
double bond2 = d_bond_distance[d_shake_type(m,1)];
double bond12 = d_angle_distance[d_shake_type(m,2)];
// r01,r02,r12 = distance vec between atoms, with PBC
// r01,r02,r12 = distance vec between atoms
double r01[3];
r01[0] = d_x(i0,0) - d_x(i1,0);
r01[1] = d_x(i0,1) - d_x(i1,1);
r01[2] = d_x(i0,2) - d_x(i1,2);
minimum_image(r01);
double r02[3];
r02[0] = d_x(i0,0) - d_x(i2,0);
r02[1] = d_x(i0,1) - d_x(i2,1);
r02[2] = d_x(i0,2) - d_x(i2,2);
minimum_image(r02);
double r12[3];
r12[0] = d_x(i1,0) - d_x(i2,0);
r12[1] = d_x(i1,1) - d_x(i2,1);
r12[2] = d_x(i1,2) - d_x(i2,2);
minimum_image(r12);
// s01,s02,s12 = distance vec after unconstrained update, with PBC
// use Domain::minimum_image_once(), not minimum_image()
// b/c xshake values might be huge, due to e.g. fix gcmc
// s01,s02,s12 = distance vec after unconstrained update
double s01[3];
s01[0] = d_xshake(i0,0) - d_xshake(i1,0);
s01[1] = d_xshake(i0,1) - d_xshake(i1,1);
s01[2] = d_xshake(i0,2) - d_xshake(i1,2);
minimum_image_once(s01);
double s02[3];
s02[0] = d_xshake(i0,0) - d_xshake(i2,0);
s02[1] = d_xshake(i0,1) - d_xshake(i2,1);
s02[2] = d_xshake(i0,2) - d_xshake(i2,2);
minimum_image_once(s02);
double s12[3];
s12[0] = d_xshake(i1,0) - d_xshake(i2,0);
s12[1] = d_xshake(i1,1) - d_xshake(i2,1);
s12[2] = d_xshake(i1,2) - d_xshake(i2,2);
minimum_image_once(s12);
// scalar distances between atoms
@ -1386,10 +1406,10 @@ void FixShakeKokkos<DeviceType>::shake3angle(int m, EV_FLOAT& ev) const
}
if (EVFLAG) {
nlist = 0;
if (i0 < nlocal) list[nlist++] = i0;
if (i1 < nlocal) list[nlist++] = i1;
if (i2 < nlocal) list[nlist++] = i2;
int count = 0;
if (i0 < nlocal) atomlist[count++] = i0;
if (i1 < nlocal) atomlist[count++] = i1;
if (i2 < nlocal) atomlist[count++] = i2;
v[0] = lamda01*r01[0]*r01[0]+lamda02*r02[0]*r02[0]+lamda12*r12[0]*r12[0];
v[1] = lamda01*r01[1]*r01[1]+lamda02*r02[1]*r02[1]+lamda12*r12[1]*r12[1];
@ -1398,10 +1418,27 @@ void FixShakeKokkos<DeviceType>::shake3angle(int m, EV_FLOAT& ev) const
v[4] = lamda01*r01[0]*r01[2]+lamda02*r02[0]*r02[2]+lamda12*r12[0]*r12[2];
v[5] = lamda01*r01[1]*r01[2]+lamda02*r02[1]*r02[2]+lamda12*r12[1]*r12[2];
v_tally<NEIGHFLAG>(ev,nlist,list,3.0,v);
v_tally<NEIGHFLAG>(ev,count,atomlist,3.0,v);
}
}
/* ----------------------------------------------------------------------
print-out bond & angle statistics
------------------------------------------------------------------------- */
template<class DeviceType>
void FixShakeKokkos<DeviceType>::stats()
{
atomKK->sync(Host,X_MASK);
k_shake_flag.sync_host();
k_shake_atom.sync_host();
k_shake_type.sync_host();
k_list.sync_host();
k_closest_list.sync_host();
FixShake::stats();
}
/* ----------------------------------------------------------------------
allocate local atom-based arrays
------------------------------------------------------------------------- */
@ -1648,6 +1685,13 @@ void FixShakeKokkos<DeviceType>::shake_end_of_step(int vflag) {
template<class DeviceType>
void FixShakeKokkos<DeviceType>::correct_coordinates(int vflag) {
x = atom->x;
v = atom->v;
f = atom->f;
mass = atom->mass;
rmass = atom->rmass;
type = atom->type;
atomKK->sync(Host,X_MASK|V_MASK|F_MASK);
// save current forces and velocities so that you
@ -1715,11 +1759,9 @@ void FixShakeKokkos<DeviceType>::correct_coordinates(int vflag) {
double **xtmp = xshake;
xshake = x;
if (comm->nprocs > 1) {
forward_comm_device = 0;
comm->forward_comm(this);
forward_comm_device = 1;
}
forward_comm_device = 0;
comm->forward_comm(this);
forward_comm_device = 1;
xshake = xtmp;
atomKK->modified(Host,X_MASK|V_MASK|F_MASK);
@ -1739,7 +1781,7 @@ void FixShakeKokkos<DeviceType>::correct_coordinates(int vflag) {
template<class DeviceType>
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
void FixShakeKokkos<DeviceType>::v_tally(EV_FLOAT &ev, int n, int *list, double total,
void FixShakeKokkos<DeviceType>::v_tally(EV_FLOAT &ev, int n, int *atomlist, double total,
double *v) const
{
int m;
@ -1759,7 +1801,7 @@ void FixShakeKokkos<DeviceType>::v_tally(EV_FLOAT &ev, int n, int *list, double
for (int i = 0; i < n; i++) {
auto v_vatom = ScatterViewHelper<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_vatom),decltype(ndup_vatom)>::get(dup_vatom,ndup_vatom);
auto a_vatom = v_vatom.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
m = list[i];
m = atomlist[i];
a_vatom(m,0) += fraction*v[0];
a_vatom(m,1) += fraction*v[1];
a_vatom(m,2) += fraction*v[2];
@ -1770,163 +1812,41 @@ void FixShakeKokkos<DeviceType>::v_tally(EV_FLOAT &ev, int n, int *list, double
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void FixShakeKokkos<DeviceType>::update_domain_variables()
{
triclinic = domain->triclinic;
xperiodic = domain->xperiodic;
xprd_half = domain->xprd_half;
xprd = domain->xprd;
yperiodic = domain->yperiodic;
yprd_half = domain->yprd_half;
yprd = domain->yprd;
zperiodic = domain->zperiodic;
zprd_half = domain->zprd_half;
zprd = domain->zprd;
xy = domain->xy;
xz = domain->xz;
yz = domain->yz;
}
/* ----------------------------------------------------------------------
minimum image convention in periodic dimensions
use 1/2 of box size as test
for triclinic, also add/subtract tilt factors in other dims as needed
changed "if" to "while" to enable distance to
far-away ghost atom returned by atom->map() to be wrapped back into box
could be problem for looking up atom IDs when cutoff > boxsize
this should not be used if atom has moved infinitely far outside box
b/c while could iterate forever
e.g. fix shake prediction of new position with highly overlapped atoms
use minimum_image_once() instead
return local index of atom J or any of its images that is closest to atom I
if J is not a valid index like -1, just return it
copied from domain.cpp
------------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixShakeKokkos<DeviceType>::minimum_image(double *delta) const
int FixShakeKokkos<DeviceType>::closest_image(const int i, int j) const
{
if (triclinic == 0) {
if (xperiodic) {
while (fabs(delta[0]) > xprd_half) {
if (delta[0] < 0.0) delta[0] += xprd;
else delta[0] -= xprd;
}
}
if (yperiodic) {
while (fabs(delta[1]) > yprd_half) {
if (delta[1] < 0.0) delta[1] += yprd;
else delta[1] -= yprd;
}
}
if (zperiodic) {
while (fabs(delta[2]) > zprd_half) {
if (delta[2] < 0.0) delta[2] += zprd;
else delta[2] -= zprd;
}
}
if (j < 0) return j;
} else {
if (zperiodic) {
while (fabs(delta[2]) > zprd_half) {
if (delta[2] < 0.0) {
delta[2] += zprd;
delta[1] += yz;
delta[0] += xz;
} else {
delta[2] -= zprd;
delta[1] -= yz;
delta[0] -= xz;
}
}
}
if (yperiodic) {
while (fabs(delta[1]) > yprd_half) {
if (delta[1] < 0.0) {
delta[1] += yprd;
delta[0] += xy;
} else {
delta[1] -= yprd;
delta[0] -= xy;
}
}
}
if (xperiodic) {
while (fabs(delta[0]) > xprd_half) {
if (delta[0] < 0.0) delta[0] += xprd;
else delta[0] -= xprd;
}
}
}
}
/* ----------------------------------------------------------------------
minimum image convention in periodic dimensions
use 1/2 of box size as test
for triclinic, also add/subtract tilt factors in other dims as needed
only shift by one box length in each direction
this should not be used if multiple box shifts are required
copied from domain.cpp
------------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixShakeKokkos<DeviceType>::minimum_image_once(double *delta) const
{
if (triclinic == 0) {
if (xperiodic) {
if (fabs(delta[0]) > xprd_half) {
if (delta[0] < 0.0) delta[0] += xprd;
else delta[0] -= xprd;
}
}
if (yperiodic) {
if (fabs(delta[1]) > yprd_half) {
if (delta[1] < 0.0) delta[1] += yprd;
else delta[1] -= yprd;
}
}
if (zperiodic) {
if (fabs(delta[2]) > zprd_half) {
if (delta[2] < 0.0) delta[2] += zprd;
else delta[2] -= zprd;
}
}
} else {
if (zperiodic) {
if (fabs(delta[2]) > zprd_half) {
if (delta[2] < 0.0) {
delta[2] += zprd;
delta[1] += yz;
delta[0] += xz;
} else {
delta[2] -= zprd;
delta[1] -= yz;
delta[0] -= xz;
}
}
}
if (yperiodic) {
if (fabs(delta[1]) > yprd_half) {
if (delta[1] < 0.0) {
delta[1] += yprd;
delta[0] += xy;
} else {
delta[1] -= yprd;
delta[0] -= xy;
}
}
}
if (xperiodic) {
if (fabs(delta[0]) > xprd_half) {
if (delta[0] < 0.0) delta[0] += xprd;
else delta[0] -= xprd;
}
const X_FLOAT xi0 = d_x(i,0);
const X_FLOAT xi1 = d_x(i,1);
const X_FLOAT xi2 = d_x(i,2);
int closest = j;
X_FLOAT delx = xi0 - d_x(j,0);
X_FLOAT dely = xi1 - d_x(j,1);
X_FLOAT delz = xi2 - d_x(j,2);
X_FLOAT rsqmin = delx*delx + dely*dely + delz*delz;
X_FLOAT rsq;
while (d_sametag[j] >= 0) {
j = d_sametag[j];
delx = xi0 - d_x(j,0);
dely = xi1 - d_x(j,1);
delz = xi2 - d_x(j,2);
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < rsqmin) {
rsqmin = rsq;
closest = j;
}
}
return closest;
}
/* ---------------------------------------------------------------------- */

19
src/KOKKOS/fix_shake_kokkos.h
View File

@ -30,6 +30,8 @@ FixStyle(shake/kk/host,FixShakeKokkos<LMPHostType>);
namespace LAMMPS_NS {
struct TagFixShakePreNeighbor{};
template<int NEIGHFLAG, int EVFLAG>
struct TagFixShakePostForce{};
@ -54,6 +56,7 @@ class FixShakeKokkos : public FixShake, public KokkosBase {
void min_setup(int) override;
void pre_neighbor() override;
void post_force(int) override;
void min_post_force(int) override;
void grow_arrays(int) override;
void copy_arrays(int, int, int) override;
@ -76,6 +79,9 @@ class FixShakeKokkos : public FixShake, public KokkosBase {
void unconstrained_update() override;
KOKKOS_INLINE_FUNCTION
void operator()(TagFixShakePreNeighbor, const int&) const;
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagFixShakePostForce<NEIGHFLAG,EVFLAG>, const int&, EV_FLOAT&) const;
@ -132,9 +138,13 @@ class FixShakeKokkos : public FixShake, public KokkosBase {
DAT::tdual_int_1d k_list;
typename AT::t_int_1d d_list; // list of clusters to SHAKE
DAT::tdual_int_2d k_closest_list;
typename AT::t_int_2d d_closest_list; // list of closest atom indices in SHAKE clusters
DAT::tdual_int_scalar k_error_flag;
DAT::tdual_int_scalar k_nlist;
void stats() override;
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
@ -191,6 +201,8 @@ class FixShakeKokkos : public FixShake, public KokkosBase {
tagint **shake_atom_tmp;
int **shake_type_tmp;
DAT::tdual_int_1d k_sametag;
typename AT::t_int_1d d_sametag;
int map_style;
DAT::tdual_int_1d k_map_array;
dual_hash_type k_map_hash;
@ -198,12 +210,7 @@ class FixShakeKokkos : public FixShake, public KokkosBase {
// copied from Domain
KOKKOS_INLINE_FUNCTION
void minimum_image(double *) const;
KOKKOS_INLINE_FUNCTION
void minimum_image_once(double *) const;
void update_domain_variables();
int closest_image(const int, int) const;
int triclinic;
int xperiodic,yperiodic,zperiodic;

67
src/LATTE/Install.sh
View File

@ -1,67 +0,0 @@
# Install/unInstall package files in LAMMPS
# mode = 0/1/2 for uninstall/install/update
mode=$1
# enforce using portable C locale
LC_ALL=C
export LC_ALL
# arg1 = file, arg2 = file it depends on
action () {
if (test $mode = 0) then
rm -f ../$1
elif (! cmp -s $1 ../$1) then
if (test -z "$2" || test -e ../$2) then
cp $1 ..
if (test $mode = 2) then
echo " updating src/$1"
fi
fi
elif (test -n "$2") then
if (test ! -e ../$2) then
rm -f ../$1
fi
fi
}
# all package files with no dependencies
for file in *.cpp *.h; do
test -f ${file} && action $file
done
# edit 2 Makefile.package files to include/exclude package info
if (test $1 = 1) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*latte[^ \t]* //' ../Makefile.package
sed -i -e 's|^PKG_INC =[ \t]*|&-I../../lib/latte/includelink |' ../Makefile.package
sed -i -e 's|^PKG_PATH =[ \t]*|&-L../../lib/latte/liblink |' ../Makefile.package
#sed -i -e 's|^PKG_LIB =[ \t]*|&-llatte |' ../Makefile.package
sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(latte_SYSINC) |' ../Makefile.package
sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(latte_SYSLIB) |' ../Makefile.package
sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(latte_SYSPATH) |' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then
sed -i -e '/^[ \t]*include.*latte.*$/d' ../Makefile.package.settings
# multiline form needed for BSD sed on Macs
sed -i -e '4 i \
include ..\/..\/lib\/latte\/Makefile.lammps
' ../Makefile.package.settings
fi
elif (test $1 = 0) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*latte[^ \t]* //' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then
sed -i -e '/^[ \t]*include.*latte.*$/d' ../Makefile.package.settings
fi
fi

23
src/LATTE/README
View File

@ -1,23 +0,0 @@
This package provides a fix latte command which is a wrapper on the
LATTE DFTB code, so that molecular dynamics can be run with LAMMPS
using density-functional tight-binding quantum forces calculated by
LATTE. More information on LATTE can be found at this web site:
https://github.com/lanl/LATTE. Its authors are Marc Cawkwell and
Anders Niklasson and Christian Negre from Los Alamos National
Laboratory (LANL). A brief technical description of LATTE is given
on the fix_latte doc page.
Using this package requires the LATTE code to first be downloaded and
built as a library on your system. This can be done in lib/latte or
elsewhere on your system. Details of the download and build process
for LATTE are given in the lib/latte/README file and it can also be
done via the make lib-latte command from the LAMMPS src directory.
Once you have successfully built LAMMPS with this package and the
LATTE library you can test it using an input file from the examples
latte dir, e.g.
lmp_serial -in lammps/examples/latte/in.latte.water
This pair style was written in collaboration with the LATTE
developers.

469
src/LATTE/fix_latte.cpp
View File

@ -1,469 +0,0 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Christian Negre (LANL)
------------------------------------------------------------------------- */
#include "fix_latte.h"
#include "atom.h"
#include "comm.h"
#include "compute.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "memory.h"
#include "modify.h"
#include "update.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
extern "C" {
void latte(int *, int *, double *, int *, int *,
double *, double *, double *, double *,
double *, double *, double *, int *,
double *, double *, double *, double *, int * , bool *);
int latte_abiversion();
}
// the ABIVERSION number here must be kept consistent
// with its counterpart in the LATTE library and the
// prototype above. We want to catch mismatches with
// a meaningful error messages, as they can cause
// difficult to debug crashes or memory corruption.
#define LATTE_ABIVERSION 20180622
/* ---------------------------------------------------------------------- */
FixLatte::FixLatte(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 3) utils::missing_cmd_args(FLERR, "fix latte", error);
if (strcmp(update->unit_style,"metal") != 0)
error->all(FLERR,"Must use units metal with fix latte command");
if (comm->nprocs != 1)
error->all(FLERR,"Fix latte currently runs only in serial");
if (LATTE_ABIVERSION != latte_abiversion())
error->all(FLERR,"LAMMPS is linked against incompatible LATTE library");
scalar_flag = 1;
global_freq = 1;
extscalar = 1;
energy_global_flag = 1;
virial_global_flag = 1;
thermo_energy = thermo_virial = 1;
// process optional args
coulomb = 0;
id_pe = nullptr;
exclude = 0;
id_exclude = nullptr;
int iarg = 3;
while (iarg < narg) {
if (strcmp(arg[iarg],"coulomb") == 0) {
if (iarg+2 > narg)
utils::missing_cmd_args(FLERR, "fix latte coulomb", error);
coulomb = 1;
error->all(FLERR,"Fix latte does not yet support LAMMPS calculation of a Coulomb potential");
delete[] id_pe;
id_pe = utils::strdup(arg[iarg+1]);
c_pe = modify->get_compute_by_id(id_pe);
if (!c_pe) error->all(FLERR,"Could not find fix latte compute ID {}", id_pe);
if (c_pe->peatomflag == 0) error->all(FLERR,"Fix latte compute ID does not compute pe/atom");
iarg += 2;
} else if (strcmp(arg[iarg],"exclude") == 0) {
if (iarg+2 > narg)
utils::missing_cmd_args(FLERR, "fix latte exclude", error);
exclude = 1;
delete[] id_exclude;
id_exclude = utils::strdup(arg[iarg+1]);
iarg += 2;
} else
error->all(FLERR, "Unknown fix latte keyword: {}", arg[iarg]);
}
// initializations
nmax = 0;
qpotential = nullptr;
flatte = nullptr;
latte_energy = 0.0;
}
/* ---------------------------------------------------------------------- */
FixLatte::~FixLatte()
{
delete[] id_pe;
delete[] id_exclude;
memory->destroy(qpotential);
memory->destroy(flatte);
}
/* ---------------------------------------------------------------------- */
int FixLatte::setmask()
{
int mask = 0;
mask |= PRE_REVERSE;
mask |= POST_FORCE;
mask |= MIN_POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixLatte::init()
{
// error checks
if (domain->dimension == 2)
error->all(FLERR,"Fix latte requires 3d problem");
if (coulomb) {
if (atom->q_flag == 0 || force->pair == nullptr || force->kspace == nullptr)
error->all(FLERR,"Fix latte cannot compute Coulomb potential");
c_pe = modify->get_compute_by_id(id_pe);
if (!c_pe) error->all(FLERR,"Fix latte could not find Coulomb compute ID {}",id_pe);
}
// must be fully periodic or fully non-periodic
if (domain->nonperiodic == 0) pbcflag = 1;
else if (!domain->xperiodic && !domain->yperiodic && !domain->zperiodic)
pbcflag = 0;
else error->all(FLERR,"Fix latte requires 3d simulation");
// create qpotential & flatte if needed
// for now, assume nlocal will never change
if (coulomb && qpotential == nullptr) {
memory->create(qpotential,atom->nlocal,"latte:qpotential");
memory->create(flatte,atom->nlocal,3,"latte:flatte");
}
// extract pointer to exclusion_group variable from id_exclude
// exclusion_group is index of a group the Fix defines
if (exclude) {
Fix *f_exclude = modify->get_fix_by_id(id_exclude);
if (!f_exclude) error->all(FLERR,"Fix latte could not find exclude fix ID {}", id_exclude);
int dim;
exclusion_group_ptr = (int *) f_exclude->extract("exclusion_group", dim);
if (!exclusion_group_ptr || dim != 0)
error->all(FLERR,"Fix latte could not query exclude_group of fix ID {}", id_exclude);
}
}
/* ---------------------------------------------------------------------- */
void FixLatte::init_list(int /*id*/, NeighList * /*ptr*/)
{
// list = ptr;
}
/* ---------------------------------------------------------------------- */
void FixLatte::setup(int vflag)
{
natoms_last = -1;
setupflag = 1;
post_force(vflag);
setupflag = 0;
}
/* ---------------------------------------------------------------------- */
void FixLatte::min_setup(int vflag)
{
natoms_last = -1;
setupflag = 1;
post_force(vflag);
setupflag = 0;
}
/* ---------------------------------------------------------------------- */
void FixLatte::setup_pre_reverse(int eflag, int vflag)
{
pre_reverse(eflag,vflag);
}
/* ----------------------------------------------------------------------
integrate electronic degrees of freedom
------------------------------------------------------------------------- */
void FixLatte::initial_integrate(int /*vflag*/) {}
/* ----------------------------------------------------------------------
store eflag, so can use it in post_force to tally per-atom energies
------------------------------------------------------------------------- */
void FixLatte::pre_reverse(int eflag, int /*vflag*/)
{
eflag_caller = eflag;
}
/* ---------------------------------------------------------------------- */
void FixLatte::post_force(int vflag)
{
int eflag = eflag_caller;
ev_init(eflag,vflag);
// compute Coulombic potential = pe[i]/q[i]
// invoke compute pe/atom
// wrap with clear/add and trigger pe/atom calculation every step
if (coulomb) {
modify->clearstep_compute();
if (!(c_pe->invoked_flag & Compute::INVOKED_PERATOM)) {
c_pe->compute_peratom();
c_pe->invoked_flag |= Compute::INVOKED_PERATOM;
}
modify->addstep_compute(update->ntimestep+1);
double *pe = c_pe->vector_atom;
double *q = atom->q;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (q[i]) qpotential[i] = pe[i]/q[i];
else qpotential[i] = 0.0;
}
// hardwire these unsupported flags for now
int coulombflag = 0;
neighflag = 0;
// set flags used by LATTE
// note that LATTE does not compute per-atom energies or virials
flags_latte[0] = pbcflag; // 1 for fully periodic, 0 for fully non-periodic
flags_latte[1] = coulombflag; // 1 for LAMMPS computes Coulombics, 0 for LATTE
flags_latte[2] = eflag_atom; // 1 to return per-atom energies, 0 for no
flags_latte[3] = vflag_global && thermo_virial; // 1 to return global/per-atom
flags_latte[4] = vflag_atom && thermo_virial; // virial, 0 for no
flags_latte[5] = neighflag; // 1 to pass neighbor list to LATTE, 0 for no
// newsystem flag determines whether LATTE treats snapshot
// as new system (more work) or increment to last system
// if setup or atom count changed then newsystem = 1
// else newsystem = 0
if (setupflag || atom->natoms != natoms_last) newsystem = 1;
else newsystem = 0;
// setup arguments for latte() function and invoke it
// either for all atoms or excluding some atoms
// in latter case, need to construct reduced-size per-atom vectors/arrays
if (!exclude) latte_wrapper_all();
else {
int anyexclude = 0;
int exclusion_group = *exclusion_group_ptr;
if (exclusion_group) {
int excludebit = group->bitmask[exclusion_group];
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & excludebit) anyexclude = 1;
}
if (!anyexclude) latte_wrapper_all();
else latte_wrapper_exclude();
}
newsystem = 0;
natoms_last = atom->natoms;
// sum LATTE forces to LAMMPS forces
// e.g. LAMMPS may compute Coulombics at some point
if (coulomb) {
double **f = atom->f;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
f[i][0] += flatte[i][0];
f[i][1] += flatte[i][1];
f[i][2] += flatte[i][2];
}
}
}
/* ----------------------------------------------------------------------
invoke LATTE on all LAMMPS atoms
------------------------------------------------------------------------- */
void FixLatte::latte_wrapper_all()
{
int natoms = atom->nlocal;
double *coords = &atom->x[0][0];
int *types = atom->type;
int ntypes = atom->ntypes;
double *mass = &atom->mass[1];
double *boxlo = domain->boxlo;
double *boxhi = domain->boxhi;
double *forces;
bool latteerror = false;
if (coulomb) forces = &flatte[0][0];
else forces = &atom->f[0][0];
int maxiter = -1;
latte(flags_latte,&natoms,coords,types,&ntypes,mass,boxlo,boxhi,
&domain->xy,&domain->xz,&domain->yz,forces,&maxiter,&latte_energy,
&atom->v[0][0],&update->dt,virial,&newsystem,&latteerror);
if (latteerror) error->all(FLERR,"Internal LATTE problem");
}
/* ----------------------------------------------------------------------
invoke LATTE on only LAMMPS atoms not in exclude group
------------------------------------------------------------------------- */
void FixLatte::latte_wrapper_exclude()
{
int m;
int exclusion_group = *exclusion_group_ptr;
int excludebit = group->bitmask[exclusion_group];
int *mask = atom->mask;
int nlocal = atom->nlocal;
// nlatte = number of non-excluded atoms to pass to LATTE
int nlatte = 0;
for (int i = 0; i < nlocal; i++)
if (!(mask[i] & excludebit)) nlatte++;
// created compressed type vector and coords array
int *typeinclude;
double **xinclude,**finclude;
memory->create(typeinclude,nlatte,"latte:typeinclude");
memory->create(xinclude,nlatte,3,"latte:xinclude");
memory->create(finclude,nlatte,3,"latte:finclude");
double *coords = &xinclude[0][0];
int *types = typeinclude;
double *forces = &finclude[0][0];
double **x = atom->x;
int *type = atom->type;
nlatte = 0;
m = 0;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & excludebit) continue;
types[nlatte] = type[i];
nlatte++;
coords[m+0] = x[i][0];
coords[m+1] = x[i][1];
coords[m+2] = x[i][2];
m += 3;
}
int ntypes = atom->ntypes;
double *mass = &atom->mass[1];
double *boxlo = domain->boxlo;
double *boxhi = domain->boxhi;
bool latteerror = false;
int maxiter = -1;
latte(flags_latte,&nlatte,coords,types,&ntypes,mass,boxlo,boxhi,
&domain->xy,&domain->xz,&domain->yz,forces,&maxiter,&latte_energy,
&atom->v[0][0],&update->dt,virial,&newsystem,&latteerror);
if (latteerror) error->all(FLERR,"Internal LATTE problem");
// expand compressed forces array
double **f = atom->f;
m = 0;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & excludebit) continue;
f[i][0] = forces[m+0];
f[i][1] = forces[m+1];
f[i][2] = forces[m+2];
m += 3;
}
memory->destroy(typeinclude);
memory->destroy(xinclude);
memory->destroy(finclude);
}
/* ---------------------------------------------------------------------- */
void FixLatte::min_post_force(int vflag)
{
post_force(vflag);
}
/* ----------------------------------------------------------------------
integrate electronic degrees of freedom
------------------------------------------------------------------------- */
void FixLatte::final_integrate() {}
/* ---------------------------------------------------------------------- */
void FixLatte::reset_dt()
{
//dtv = update->dt;
//dtf = 0.5 * update->dt * force->ftm2v;
}
/* ----------------------------------------------------------------------
DFTB energy from LATTE
------------------------------------------------------------------------- */
double FixLatte::compute_scalar()
{
return latte_energy;
}
/* ----------------------------------------------------------------------
memory usage of local arrays
------------------------------------------------------------------------- */
double FixLatte::memory_usage()
{
double bytes = 0.0;
if (coulomb) bytes += (double)nmax * sizeof(double);
if (coulomb) bytes += (double)nmax*3 * sizeof(double);
return bytes;
}

72
src/LATTE/fix_latte.h
View File

@ -1,72 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
// clang-format off
FixStyle(latte,FixLatte);
// clang-format on
#else
#ifndef LMP_FIX_LATTE_H
#define LMP_FIX_LATTE_H
#include "fix.h"
namespace LAMMPS_NS {
class FixLatte : public Fix {
public:
FixLatte(class LAMMPS *, int, char **);
~FixLatte() override;
int setmask() override;
void init() override;
void init_list(int, class NeighList *) override;
void setup(int) override;
void min_setup(int) override;
void setup_pre_reverse(int, int) override;
void initial_integrate(int) override;
void pre_reverse(int, int) override;
void post_force(int) override;
void min_post_force(int) override;
void final_integrate() override;
void reset_dt() override;
double compute_scalar() override;
double memory_usage() override;
protected:
int coulomb, pbcflag, pe_peratom, virial_global, virial_peratom, neighflag;
int exclude, excludebit;
int eflag_caller;
char *id_pe, *id_exclude;
int *exclusion_group_ptr;
int setupflag, newsystem;
bigint natoms_last;
int flags_latte[6];
int nmax;
double *qpotential;
double **flatte;
double latte_energy;
class NeighList *list;
class Compute *c_pe;
void latte_wrapper_all();
void latte_wrapper_exclude();
};
} // namespace LAMMPS_NS
#endif
#endif

35
src/MC/fix_charge_regulation.cpp
View File

@ -79,6 +79,8 @@ FixChargeRegulation::FixChargeRegulation(LAMMPS *lmp, int narg, char **arg) :
c_pe(nullptr), random_equal(nullptr), random_unequal(nullptr),
idftemp(nullptr)
{
if (narg < 5) utils::missing_cmd_args(FLERR, "fix charge/regulation", error);
if (lmp->citeme) lmp->citeme->add(cite_fix_charge_regulation);
// Region restrictions not yet implemented ..
@ -151,6 +153,16 @@ FixChargeRegulation::~FixChargeRegulation() {
delete[] pHstr;
delete[] idftemp;
// delete exclusion group created in init()
// unset neighbor exclusion settings made in init()
// not necessary if group and neighbor classes already destroyed
// when LAMMPS exits
if (exclusion_group_bit && group) {
auto group_id = std::string("FixChargeRegulation:gcmc_exclusion_group:") + id;
group->assign(group_id + " delete");
}
if (group) {
int igroupall = group->find("all");
neighbor->exclusion_group_group_delete(exclusion_group, igroupall);
@ -211,12 +223,11 @@ void FixChargeRegulation::init() {
// create unique group name for atoms to be excluded
auto group_id = fmt::format("FixChargeRegulation:exclusion_group:{}",id);
auto group_id = std::string("FixChargeRegulation:exclusion_group:") + id;
group->assign(group_id + " subtract all all");
exclusion_group = group->find(group_id);
if (exclusion_group == -1)
error->all(FLERR,"Could not find fix charge/regulation exclusion "
"group ID");
error->all(FLERR,"Could not find fix charge/regulation exclusion group ID");
exclusion_group_bit = group->bitmask[exclusion_group];
// neighbor list exclusion setup
@ -240,8 +251,7 @@ void FixChargeRegulation::init() {
MPI_Allreduce(&flag, &flagall, 1, MPI_INT, MPI_SUM, world);
if (flagall)
error->all(FLERR, "Cannot use fix charge/regulation on atoms "
"in atom_modify first group");
error->all(FLERR, "Cannot use fix charge/regulation on atoms in atom_modify first group");
}
// construct group bitmask for all new atoms
@ -1284,16 +1294,13 @@ void FixChargeRegulation::restart(char *buf)
/* ---------------------------------------------------------------------- */
void FixChargeRegulation::setThermoTemperaturePointer() {
int ifix = -1;
ifix = modify->find_fix(idftemp);
if (ifix == -1) {
error->all(FLERR, "fix charge/regulation regulation could not find "
"a temperature fix id provided by tempfixid\n");
}
Fix *temperature_fix = modify->fix[ifix];
int dim;
target_temperature_tcp = (double *) temperature_fix->extract("t_target", dim);
Fix *ifix = modify->get_fix_by_id(idftemp);
if (!ifix)
error->all(FLERR, "fix charge/regulation could not find thermostat fix id {}", idftemp);
int dim;
target_temperature_tcp = (double *) ifix->extract("t_target", dim);
if (!target_temperature_tcp) error->all(FLERR, "Fix id {} does not control temperature", idftemp);
}
/* ---------------------------------------------------------------------- */

58
src/MC/fix_gcmc.cpp
View File

@ -61,19 +61,19 @@ using namespace MathConst;
#define MAXENERGYTEST 1.0e50
enum{EXCHATOM,EXCHMOL}; // exchmode
enum{NONE,MOVEATOM,MOVEMOL}; // movemode
enum { EXCHATOM, EXCHMOL }; // exchmode
enum { NONE, MOVEATOM, MOVEMOL }; // movemode
/* ---------------------------------------------------------------------- */
FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
region(nullptr), idregion(nullptr), full_flag(false), groupstrings(nullptr),
grouptypestrings(nullptr), grouptypebits(nullptr), grouptypes(nullptr), local_gas_list(nullptr),
molcoords(nullptr), molq(nullptr), molimage(nullptr), random_equal(nullptr), random_unequal(nullptr),
fixrigid(nullptr), fixshake(nullptr), idrigid(nullptr), idshake(nullptr)
Fix(lmp, narg, arg), region(nullptr), idregion(nullptr), full_flag(false),
groupstrings(nullptr), grouptypestrings(nullptr), grouptypebits(nullptr), grouptypes(nullptr),
local_gas_list(nullptr), molcoords(nullptr), molq(nullptr), molimage(nullptr),
random_equal(nullptr), random_unequal(nullptr), fixrigid(nullptr), fixshake(nullptr),
idrigid(nullptr), idshake(nullptr)
{
if (narg < 11) error->all(FLERR,"Illegal fix gcmc command");
if (narg < 11) utils::missing_cmd_args(FLERR, "fix gcmc", error);
if (atom->molecular == Atom::TEMPLATE)
error->all(FLERR,"Fix gcmc does not (yet) work with atom_style template");
@ -418,11 +418,27 @@ FixGCMC::~FixGCMC()
delete[] grouptypestrings[igroup];
memory->sfree(grouptypestrings);
}
if (full_flag && group) {
// delete exclusion group created in init()
// delete molecule group created in init()
// unset neighbor exclusion settings made in init()
// not necessary if group and neighbor classes already destroyed
// when LAMMPS exits
if (exclusion_group_bit && group) {
auto group_id = std::string("FixGCMC:gcmc_exclusion_group:") + id;
group->assign(group_id + " delete");
}
if (molecule_group_bit && group) {
auto group_id = std::string("FixGCMC:rotation_gas_atoms:") + id;
group->assign(group_id + " delete");
}
if (full_flag && group && neighbor) {
int igroupall = group->find("all");
neighbor->exclusion_group_group_delete(exclusion_group,igroupall);
}
}
/* ---------------------------------------------------------------------- */
@ -438,7 +454,6 @@ int FixGCMC::setmask()
void FixGCMC::init()
{
// set index and check validity of region
if (idregion) {
@ -450,14 +465,14 @@ void FixGCMC::init()
if (triclinic) {
if ((region_xlo < domain->boxlo_bound[0]) || (region_xhi > domain->boxhi_bound[0]) ||
(region_ylo < domain->boxlo_bound[1]) || (region_yhi > domain->boxhi_bound[1]) ||
(region_zlo < domain->boxlo_bound[2]) || (region_zhi > domain->boxhi_bound[2])) {
(region_ylo < domain->boxlo_bound[1]) || (region_yhi > domain->boxhi_bound[1]) ||
(region_zlo < domain->boxlo_bound[2]) || (region_zhi > domain->boxhi_bound[2])) {
error->all(FLERR,"Fix gcmc region extends outside simulation box");
}
} else {
if ((region_xlo < domain->boxlo[0]) || (region_xhi > domain->boxhi[0]) ||
(region_ylo < domain->boxlo[1]) || (region_yhi > domain->boxhi[1]) ||
(region_zlo < domain->boxlo[2]) || (region_zhi > domain->boxhi[2]))
(region_ylo < domain->boxlo[1]) || (region_yhi > domain->boxhi[1]) ||
(region_zlo < domain->boxlo[2]) || (region_zhi > domain->boxhi[2]))
error->all(FLERR,"Fix gcmc region extends outside simulation box");
}
@ -550,9 +565,8 @@ void FixGCMC::init()
fixrigid = nullptr;
if (rigidflag) {
int ifix = modify->find_fix(idrigid);
if (ifix < 0) error->all(FLERR,"Fix gcmc rigid fix does not exist");
fixrigid = modify->fix[ifix];
fixrigid = modify->get_fix_by_id(idrigid);
if (!fixrigid) error->all(FLERR,"Fix gcmc rigid fix ID {} does not exist", idrigid);
int tmp;
if (&onemols[imol] != (Molecule **) fixrigid->extract("onemol",tmp))
error->all(FLERR, "Fix gcmc and fix rigid/small not using same molecule template ID");
@ -563,13 +577,11 @@ void FixGCMC::init()
fixshake = nullptr;
if (shakeflag) {
int ifix = modify->find_fix(idshake);
if (ifix < 0) error->all(FLERR,"Fix gcmc shake fix does not exist");
fixshake = modify->fix[ifix];
fixshake = modify->get_fix_by_id(idshake);
if (!fixshake) error->all(FLERR,"Fix gcmc shake fix ID {} does not exist", idshake);
int tmp;
if (&onemols[imol] != (Molecule **) fixshake->extract("onemol",tmp))
error->all(FLERR,"Fix gcmc and fix shake not using "
"same molecule template ID");
error->all(FLERR,"Fix gcmc and fix shake not using same molecule template ID");
}
if (domain->dimension == 2)

25
src/MC/fix_widom.cpp
View File

@ -62,7 +62,7 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) :
if (narg < 8) utils::missing_cmd_args(FLERR, "fix widom", error);
if (atom->molecular == Atom::TEMPLATE)
error->all(FLERR,"Fix widom does not (yet) work with atom_style template");
error->all(FLERR, "Fix widom does not (yet) work with atom_style template");
dynamic_group_allow = 1;
@ -247,6 +247,26 @@ FixWidom::~FixWidom()
memory->destroy(molq);
memory->destroy(molimage);
// delete exclusion group created in init()
// delete molecule group created in init()
// unset neighbor exclusion settings made in init()
// not necessary if group and neighbor classes already destroyed
// when LAMMPS exits
if (exclusion_group_bit && group) {
auto group_id = std::string("FixWidom:widom_exclusion_group:") + id;
group->assign(group_id + " delete");
}
if (molecule_group_bit && group) {
auto group_id = std::string("FixWidom:rotation_gas_atoms:") + id;
group->assign(group_id + " delete");
}
if (full_flag && group && neighbor) {
int igroupall = group->find("all");
neighbor->exclusion_group_group_delete(exclusion_group,igroupall);
}
}
/* ---------------------------------------------------------------------- */
@ -403,7 +423,8 @@ void FixWidom::init()
// full_flag on molecules on more than one processor.
// Print error if this is the current mode
if (full_flag && (exchmode == EXCHMOL) && comm->nprocs > 1)
error->all(FLERR,"fix widom does currently not support full_energy option with molecules on more than 1 MPI process.");
error->all(FLERR,"fix widom does currently not support full_energy option with "
"molecules on more than 1 MPI process.");
}

11
src/MDI/fix_mdi_qmmm.cpp
View File

@ -585,7 +585,6 @@ void FixMDIQMMM::pre_force(int vflag)
{
int ilocal;
double rsq;
double delta[3];
// invoke pair hybrid sub-style pair coulomb and Kspace directly
// set eflag = 2 so they calculate per-atom energy
@ -655,11 +654,11 @@ void FixMDIQMMM::pre_force(int vflag)
if (ilocal >= 0) {
for (int j = 0; j < nqm; j++) {
if (j == i) continue;
delta[0] = xqm[i][0] - xqm[j][0];
delta[1] = xqm[i][1] - xqm[j][1];
delta[2] = xqm[i][2] - xqm[j][2];
domain->minimum_image_once(delta);
rsq = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2];
double delx = xqm[i][0] - xqm[j][0];
double dely = xqm[i][1] - xqm[j][1];
double delz = xqm[i][2] - xqm[j][2];
domain->minimum_image(delx, dely, delz);
rsq = delx*delx + dely*dely + delz*delz;
qpotential_mine[i] -= qqrd2e * qqm[j] / sqrt(rsq);
}
}

12
src/MISC/fix_imd.cpp
View File

@ -1198,11 +1198,13 @@ void * imdsock_create() {
int imdsock_bind(void * v, int port) {
auto s = (imdsocket *) v;
memset(&(s->addr), 0, sizeof(s->addr));
s->addr.sin_family = PF_INET;
s->addr.sin_port = htons(port);
auto *addr = &(s->addr);
s->addrlen = sizeof(s->addr);
memset(addr, 0, s->addrlen);
addr->sin_family = PF_INET;
addr->sin_port = htons(port);
return bind(s->sd, (struct sockaddr *) &s->addr, sizeof(s->addr));
return bind(s->sd, (struct sockaddr *) addr, s->addrlen);
}
int imdsock_listen(void * v) {
@ -1227,7 +1229,7 @@ void *imdsock_accept(void * v) {
int len;
#endif
len = sizeof(s->addr);
len = s->addrlen;
rc = accept(s->sd, (struct sockaddr *) &s->addr, ( _SOCKLEN_TYPE * ) &len);
if (rc >= 0) {
new_s = (imdsocket *) malloc(sizeof(imdsocket));

3
src/Makefile
View File

@ -87,7 +87,6 @@ PACKAGE = \
kokkos \
kspace \
latboltz \
latte \
lepton \
machdyn \
manifold \
@ -213,7 +212,6 @@ PACKLIB = \
gpu \
kim \
kokkos \
latte \
lepton \
mpiio \
mscg \
@ -249,7 +247,6 @@ PACKEXT = \
adios \
h5md \
kim \
latte \
machdyn \
ml-hdnnp \
ml-pace \

42
src/OPENMP/thr_omp.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* -------------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -17,42 +16,37 @@
OpenMP based threading support for LAMMPS
------------------------------------------------------------------------- */
#include <cstring>
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "modify.h"
#include "neighbor.h"
#include "thr_omp.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "atom.h"
#include "bond.h"
#include "comm.h"
#include "compute.h"
#include "dihedral.h"
#include "error.h"
#include "force.h"
#include "improper.h"
#include "math_const.h"
#include "modify.h"
#include "neighbor.h"
#include "pair.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
using MathConst::THIRD;
/* ---------------------------------------------------------------------- */
ThrOMP::ThrOMP(LAMMPS *ptr, int style)
: lmp(ptr), fix(nullptr), thr_style(style), thr_error(0)
ThrOMP::ThrOMP(LAMMPS *ptr, int style) : lmp(ptr), fix(nullptr), thr_style(style), thr_error(0)
{
// register fix omp with this class
int ifix = lmp->modify->find_fix("package_omp");
if (ifix < 0)
lmp->error->all(FLERR,"The 'package omp' command is required for /omp styles");
fix = static_cast<FixOMP *>(lmp->modify->fix[ifix]);
fix = static_cast<FixOMP *>(lmp->modify->get_fix_by_id("package_omp"));
if (!fix) lmp->error->all(FLERR, "The 'package omp' command is required for /omp styles");
}
// clang-format off
/* ----------------------------------------------------------------------
Hook up per thread per atom arrays into the tally infrastructure
---------------------------------------------------------------------- */
@ -154,7 +148,6 @@ void ThrOMP::ev_setup_thr(int eflag, int vflag, int nall, double *eatom,
memset(&(thr->cvatom_imprp[0][0]),0,nall*9*sizeof(double));
}
}
// nothing to do for THR_KSPACE
}
@ -1312,7 +1305,6 @@ void ThrOMP::ev_tally_thr(Dihedral * const dihed, const int i1, const int i2,
if (i4 < nlocal) v_tally9(thr->cvatom_dihed[i4],v4);
}
}
}
/* ----------------------------------------------------------------------

3
src/Purge.list
View File

@ -51,6 +51,9 @@ lmpinstalledpkgs.h
lmpgitversion.h
mliap_model_python_couple.cpp
mliap_model_python_couple.h
# removed on 29 March 2023
fix_latte.cpp
fix_latte.h
# renamed on 23 February 2023
fix_pimd.cpp
fix_pimd.h

256
src/RIGID/fix_shake.cpp
View File

@ -55,7 +55,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
b_count(nullptr), b_count_all(nullptr), b_ave(nullptr), b_max(nullptr), b_min(nullptr),
b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), a_count(nullptr),
a_count_all(nullptr), a_ave(nullptr), a_max(nullptr), a_min(nullptr), a_ave_all(nullptr),
a_max_all(nullptr), a_min_all(nullptr), atommols(nullptr), onemols(nullptr)
a_max_all(nullptr), a_min_all(nullptr), atommols(nullptr), onemols(nullptr), closest_list(nullptr)
{
energy_global_flag = energy_peratom_flag = 1;
virial_global_flag = virial_peratom_flag = 1;
@ -244,6 +244,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
maxlist = 0;
list = nullptr;
closest_list = nullptr;
}
/* ---------------------------------------------------------------------- */
@ -319,6 +320,7 @@ FixShake::~FixShake()
}
memory->destroy(list);
memory->destroy(closest_list);
}
/* ---------------------------------------------------------------------- */
@ -546,6 +548,8 @@ void FixShake::pre_neighbor()
maxlist = nlocal;
memory->destroy(list);
memory->create(list,maxlist,"shake:list");
memory->destroy(closest_list);
memory->create(closest_list,maxlist,4,"shake:closest_list");
}
// build list of SHAKE clusters I compute
@ -560,7 +564,14 @@ void FixShake::pre_neighbor()
if (atom1 == -1 || atom2 == -1)
error->one(FLERR,"Shake atoms {} {} missing on proc {} at step {}",shake_atom[i][0],
shake_atom[i][1],comm->me,update->ntimestep);
if (i <= atom1 && i <= atom2) list[nlist++] = i;
atom1 = domain->closest_image(i, atom1);
atom2 = domain->closest_image(i, atom2);
if (i <= atom1 && i <= atom2) {
list[nlist] = i;
closest_list[nlist][0] = atom1;
closest_list[nlist][1] = atom2;
nlist++;
}
} else if (shake_flag[i] % 2 == 1) {
atom1 = atom->map(shake_atom[i][0]);
atom2 = atom->map(shake_atom[i][1]);
@ -569,7 +580,16 @@ void FixShake::pre_neighbor()
error->one(FLERR,"Shake atoms {} {} {} missing on proc {} at step {}",shake_atom[i][0],
shake_atom[i][1],shake_atom[i][2],
comm->me,update->ntimestep);
if (i <= atom1 && i <= atom2 && i <= atom3) list[nlist++] = i;
atom1 = domain->closest_image(i, atom1);
atom2 = domain->closest_image(i, atom2);
atom3 = domain->closest_image(i, atom3);
if (i <= atom1 && i <= atom2 && i <= atom3) {
list[nlist] = i;
closest_list[nlist][0] = atom1;
closest_list[nlist][1] = atom2;
closest_list[nlist][2] = atom3;
nlist++;
}
} else {
atom1 = atom->map(shake_atom[i][0]);
atom2 = atom->map(shake_atom[i][1]);
@ -579,8 +599,18 @@ void FixShake::pre_neighbor()
error->one(FLERR,"Shake atoms {} {} {} {} missing on proc {} at step {}",shake_atom[i][0],
shake_atom[i][1],shake_atom[i][2],
shake_atom[i][3],comm->me,update->ntimestep);
if (i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4)
list[nlist++] = i;
atom1 = domain->closest_image(i, atom1);
atom2 = domain->closest_image(i, atom2);
atom3 = domain->closest_image(i, atom3);
atom4 = domain->closest_image(i, atom4);
if (i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4) {
list[nlist] = i;
closest_list[nlist][0] = atom1;
closest_list[nlist][1] = atom2;
closest_list[nlist][2] = atom3;
closest_list[nlist][3] = atom4;
nlist++;
}
}
}
}
@ -597,7 +627,7 @@ void FixShake::post_force(int vflag)
// communicate results if necessary
unconstrained_update();
if (comm->nprocs > 1) comm->forward_comm(this);
comm->forward_comm(this);
// virial setup
@ -610,10 +640,10 @@ void FixShake::post_force(int vflag)
int m;
for (int i = 0; i < nlist; i++) {
m = list[i];
if (shake_flag[m] == 2) shake(m);
else if (shake_flag[m] == 3) shake3(m);
else if (shake_flag[m] == 4) shake4(m);
else shake3angle(m);
if (shake_flag[m] == 2) shake(i);
else if (shake_flag[m] == 3) shake3(i);
else if (shake_flag[m] == 4) shake4(i);
else shake3angle(i);
}
// store vflag for coordinate_constraints_end_of_step()
@ -643,7 +673,7 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop)
// communicate results if necessary
unconstrained_update_respa(ilevel);
if (comm->nprocs > 1) comm->forward_comm(this);
comm->forward_comm(this);
// virial setup only needed on last iteration of innermost level
// and if pressure is requested
@ -658,10 +688,10 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop)
int m;
for (int i = 0; i < nlist; i++) {
m = list[i];
if (shake_flag[m] == 2) shake(m);
else if (shake_flag[m] == 3) shake3(m);
else if (shake_flag[m] == 4) shake4(m);
else shake3angle(m);
if (shake_flag[m] == 2) shake(i);
else if (shake_flag[m] == 3) shake3(i);
else if (shake_flag[m] == 4) shake4(i);
else shake3angle(i);
}
// store vflag for coordinate_constraints_end_of_step()
@ -705,18 +735,18 @@ void FixShake::min_post_force(int vflag)
for (int i = 0; i < nlist; i++) {
int m = list[i];
if (shake_flag[m] == 2) {
bond_force(shake_atom[m][0], shake_atom[m][1], bond_distance[shake_type[m][0]]);
bond_force(closest_list[i][0], closest_list[i][1], bond_distance[shake_type[m][0]]);
} else if (shake_flag[m] == 3) {
bond_force(shake_atom[m][0], shake_atom[m][1], bond_distance[shake_type[m][0]]);
bond_force(shake_atom[m][0], shake_atom[m][2], bond_distance[shake_type[m][1]]);
bond_force(closest_list[i][0], closest_list[i][1], bond_distance[shake_type[m][0]]);
bond_force(closest_list[i][0], closest_list[i][2], bond_distance[shake_type[m][1]]);
} else if (shake_flag[m] == 4) {
bond_force(shake_atom[m][0], shake_atom[m][1], bond_distance[shake_type[m][0]]);
bond_force(shake_atom[m][0], shake_atom[m][2], bond_distance[shake_type[m][1]]);
bond_force(shake_atom[m][0], shake_atom[m][3], bond_distance[shake_type[m][2]]);
bond_force(closest_list[i][0], closest_list[i][1], bond_distance[shake_type[m][0]]);
bond_force(closest_list[i][0], closest_list[i][2], bond_distance[shake_type[m][1]]);
bond_force(closest_list[i][0], closest_list[i][3], bond_distance[shake_type[m][2]]);
} else {
bond_force(shake_atom[m][0], shake_atom[m][1], bond_distance[shake_type[m][0]]);
bond_force(shake_atom[m][0], shake_atom[m][2], bond_distance[shake_type[m][1]]);
bond_force(shake_atom[m][1], shake_atom[m][2], angle_distance[shake_type[m][2]]);
bond_force(closest_list[i][0], closest_list[i][1], bond_distance[shake_type[m][0]]);
bond_force(closest_list[i][0], closest_list[i][2], bond_distance[shake_type[m][1]]);
bond_force(closest_list[i][1], closest_list[i][2], angle_distance[shake_type[m][2]]);
}
}
}
@ -1737,37 +1767,36 @@ void FixShake::unconstrained_update_respa(int ilevel)
}
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
calculate SHAKE constraint forces for size 2 cluster = single bond
------------------------------------------------------------------------- */
void FixShake::shake(int m)
void FixShake::shake(int ilist)
{
int nlist,list[2];
int atomlist[2];
double v[6];
double invmass0,invmass1;
// local atom IDs and constraint distances
int i0 = atom->map(shake_atom[m][0]);
int i1 = atom->map(shake_atom[m][1]);
int m = list[ilist];
int i0 = closest_list[ilist][0];
int i1 = closest_list[ilist][1];
double bond1 = bond_distance[shake_type[m][0]];
// r01 = distance vec between atoms, with PBC
// r01 = distance vec between atoms
double r01[3];
r01[0] = x[i0][0] - x[i1][0];
r01[1] = x[i0][1] - x[i1][1];
r01[2] = x[i0][2] - x[i1][2];
domain->minimum_image(r01);
// s01 = distance vec after unconstrained update, with PBC
// use Domain::minimum_image_once(), not minimum_image()
// b/c xshake values might be huge, due to e.g. fix gcmc
// s01 = distance vec after unconstrained update
double s01[3];
s01[0] = xshake[i0][0] - xshake[i1][0];
s01[1] = xshake[i0][1] - xshake[i1][1];
s01[2] = xshake[i0][2] - xshake[i1][2];
domain->minimum_image_once(s01);
// scalar distances between atoms
@ -1823,9 +1852,9 @@ void FixShake::shake(int m)
}
if (evflag) {
nlist = 0;
if (i0 < nlocal) list[nlist++] = i0;
if (i1 < nlocal) list[nlist++] = i1;
int count = 0;
if (i0 < nlocal) atomlist[count++] = i0;
if (i1 < nlocal) atomlist[count++] = i1;
v[0] = lamda*r01[0]*r01[0];
v[1] = lamda*r01[1]*r01[1];
@ -1837,55 +1866,52 @@ void FixShake::shake(int m)
double fpairlist[] = {lamda};
double dellist[][3] = {{r01[0], r01[1], r01[2]}};
int pairlist[][2] = {{i0,i1}};
v_tally(nlist,list,2.0,v,nlocal,1,pairlist,fpairlist,dellist);
v_tally(count,atomlist,2.0,v,nlocal,1,pairlist,fpairlist,dellist);
}
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
calculate SHAKE constraint forces for size 3 cluster = two bonds
------------------------------------------------------------------------- */
void FixShake::shake3(int m)
void FixShake::shake3(int ilist)
{
int nlist,list[3];
int atomlist[3];
double v[6];
double invmass0,invmass1,invmass2;
// local atom IDs and constraint distances
int i0 = atom->map(shake_atom[m][0]);
int i1 = atom->map(shake_atom[m][1]);
int i2 = atom->map(shake_atom[m][2]);
int m = list[ilist];
int i0 = closest_list[ilist][0];
int i1 = closest_list[ilist][1];
int i2 = closest_list[ilist][2];
double bond1 = bond_distance[shake_type[m][0]];
double bond2 = bond_distance[shake_type[m][1]];
// r01,r02 = distance vec between atoms, with PBC
// r01,r02 = distance vec between atoms
double r01[3];
r01[0] = x[i0][0] - x[i1][0];
r01[1] = x[i0][1] - x[i1][1];
r01[2] = x[i0][2] - x[i1][2];
domain->minimum_image(r01);
double r02[3];
r02[0] = x[i0][0] - x[i2][0];
r02[1] = x[i0][1] - x[i2][1];
r02[2] = x[i0][2] - x[i2][2];
domain->minimum_image(r02);
// s01,s02 = distance vec after unconstrained update, with PBC
// use Domain::minimum_image_once(), not minimum_image()
// b/c xshake values might be huge, due to e.g. fix gcmc
// s01,s02 = distance vec after unconstrained update
double s01[3];
s01[0] = xshake[i0][0] - xshake[i1][0];
s01[1] = xshake[i0][1] - xshake[i1][1];
s01[2] = xshake[i0][2] - xshake[i1][2];
domain->minimum_image_once(s01);
double s02[3];
s02[0] = xshake[i0][0] - xshake[i2][0];
s02[1] = xshake[i0][1] - xshake[i2][1];
s02[2] = xshake[i0][2] - xshake[i2][2];
domain->minimum_image_once(s02);
// scalar distances between atoms
@ -1998,10 +2024,10 @@ void FixShake::shake3(int m)
}
if (evflag) {
nlist = 0;
if (i0 < nlocal) list[nlist++] = i0;
if (i1 < nlocal) list[nlist++] = i1;
if (i2 < nlocal) list[nlist++] = i2;
int count = 0;
if (i0 < nlocal) atomlist[count++] = i0;
if (i1 < nlocal) atomlist[count++] = i1;
if (i2 < nlocal) atomlist[count++] = i2;
v[0] = lamda01*r01[0]*r01[0] + lamda02*r02[0]*r02[0];
v[1] = lamda01*r01[1]*r01[1] + lamda02*r02[1]*r02[1];
@ -2014,69 +2040,64 @@ void FixShake::shake3(int m)
double dellist[][3] = {{r01[0], r01[1], r01[2]},
{r02[0], r02[1], r02[2]}};
int pairlist[][2] = {{i0,i1}, {i0,i2}};
v_tally(nlist,list,3.0,v,nlocal,2,pairlist,fpairlist,dellist);
v_tally(count,atomlist,3.0,v,nlocal,2,pairlist,fpairlist,dellist);
}
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
calculate SHAKE constraint forces for size 4 cluster = three bonds
------------------------------------------------------------------------- */
void FixShake::shake4(int m)
void FixShake::shake4(int ilist)
{
int nlist,list[4];
int atomlist[4];
double v[6];
double invmass0,invmass1,invmass2,invmass3;
// local atom IDs and constraint distances
int i0 = atom->map(shake_atom[m][0]);
int i1 = atom->map(shake_atom[m][1]);
int i2 = atom->map(shake_atom[m][2]);
int i3 = atom->map(shake_atom[m][3]);
int m = list[ilist];
int i0 = closest_list[ilist][0];
int i1 = closest_list[ilist][1];
int i2 = closest_list[ilist][2];
int i3 = closest_list[ilist][3];
double bond1 = bond_distance[shake_type[m][0]];
double bond2 = bond_distance[shake_type[m][1]];
double bond3 = bond_distance[shake_type[m][2]];
// r01,r02,r03 = distance vec between atoms, with PBC
// r01,r02,r03 = distance vec between atoms
double r01[3];
r01[0] = x[i0][0] - x[i1][0];
r01[1] = x[i0][1] - x[i1][1];
r01[2] = x[i0][2] - x[i1][2];
domain->minimum_image(r01);
double r02[3];
r02[0] = x[i0][0] - x[i2][0];
r02[1] = x[i0][1] - x[i2][1];
r02[2] = x[i0][2] - x[i2][2];
domain->minimum_image(r02);
double r03[3];
r03[0] = x[i0][0] - x[i3][0];
r03[1] = x[i0][1] - x[i3][1];
r03[2] = x[i0][2] - x[i3][2];
domain->minimum_image(r03);
// s01,s02,s03 = distance vec after unconstrained update, with PBC
// use Domain::minimum_image_once(), not minimum_image()
// b/c xshake values might be huge, due to e.g. fix gcmc
// s01,s02,s03 = distance vec after unconstrained update
double s01[3];
s01[0] = xshake[i0][0] - xshake[i1][0];
s01[1] = xshake[i0][1] - xshake[i1][1];
s01[2] = xshake[i0][2] - xshake[i1][2];
domain->minimum_image_once(s01);
double s02[3];
s02[0] = xshake[i0][0] - xshake[i2][0];
s02[1] = xshake[i0][1] - xshake[i2][1];
s02[2] = xshake[i0][2] - xshake[i2][2];
domain->minimum_image_once(s02);
double s03[3];
s03[0] = xshake[i0][0] - xshake[i3][0];
s03[1] = xshake[i0][1] - xshake[i3][1];
s03[2] = xshake[i0][2] - xshake[i3][2];
domain->minimum_image_once(s03);
// scalar distances between atoms
@ -2253,11 +2274,11 @@ void FixShake::shake4(int m)
}
if (evflag) {
nlist = 0;
if (i0 < nlocal) list[nlist++] = i0;
if (i1 < nlocal) list[nlist++] = i1;
if (i2 < nlocal) list[nlist++] = i2;
if (i3 < nlocal) list[nlist++] = i3;
int count = 0;
if (i0 < nlocal) atomlist[count++] = i0;
if (i1 < nlocal) atomlist[count++] = i1;
if (i2 < nlocal) atomlist[count++] = i2;
if (i3 < nlocal) atomlist[count++] = i3;
v[0] = lamda01*r01[0]*r01[0]+lamda02*r02[0]*r02[0]+lamda03*r03[0]*r03[0];
v[1] = lamda01*r01[1]*r01[1]+lamda02*r02[1]*r02[1]+lamda03*r03[1]*r03[1];
@ -2271,68 +2292,63 @@ void FixShake::shake4(int m)
{r02[0], r02[1], r02[2]},
{r03[0], r03[1], r03[2]}};
int pairlist[][2] = {{i0,i1}, {i0,i2}, {i0,i3}};
v_tally(nlist,list,4.0,v,nlocal,3,pairlist,fpairlist,dellist);
v_tally(count,atomlist,4.0,v,nlocal,3,pairlist,fpairlist,dellist);
}
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
calculate SHAKE constraint forces for size 3 cluster = two bonds + angle
------------------------------------------------------------------------- */
void FixShake::shake3angle(int m)
void FixShake::shake3angle(int ilist)
{
int nlist,list[3];
int atomlist[3];
double v[6];
double invmass0,invmass1,invmass2;
// local atom IDs and constraint distances
int i0 = atom->map(shake_atom[m][0]);
int i1 = atom->map(shake_atom[m][1]);
int i2 = atom->map(shake_atom[m][2]);
int m = list[ilist];
int i0 = closest_list[ilist][0];
int i1 = closest_list[ilist][1];
int i2 = closest_list[ilist][2];
double bond1 = bond_distance[shake_type[m][0]];
double bond2 = bond_distance[shake_type[m][1]];
double bond12 = angle_distance[shake_type[m][2]];
// r01,r02,r12 = distance vec between atoms, with PBC
// r01,r02,r12 = distance vec between atoms
double r01[3];
r01[0] = x[i0][0] - x[i1][0];
r01[1] = x[i0][1] - x[i1][1];
r01[2] = x[i0][2] - x[i1][2];
domain->minimum_image(r01);
double r02[3];
r02[0] = x[i0][0] - x[i2][0];
r02[1] = x[i0][1] - x[i2][1];
r02[2] = x[i0][2] - x[i2][2];
domain->minimum_image(r02);
double r12[3];
r12[0] = x[i1][0] - x[i2][0];
r12[1] = x[i1][1] - x[i2][1];
r12[2] = x[i1][2] - x[i2][2];
domain->minimum_image(r12);
// s01,s02,s12 = distance vec after unconstrained update, with PBC
// use Domain::minimum_image_once(), not minimum_image()
// b/c xshake values might be huge, due to e.g. fix gcmc
// s01,s02,s12 = distance vec after unconstrained update
double s01[3];
s01[0] = xshake[i0][0] - xshake[i1][0];
s01[1] = xshake[i0][1] - xshake[i1][1];
s01[2] = xshake[i0][2] - xshake[i1][2];
domain->minimum_image_once(s01);
double s02[3];
s02[0] = xshake[i0][0] - xshake[i2][0];
s02[1] = xshake[i0][1] - xshake[i2][1];
s02[2] = xshake[i0][2] - xshake[i2][2];
domain->minimum_image_once(s02);
double s12[3];
s12[0] = xshake[i1][0] - xshake[i2][0];
s12[1] = xshake[i1][1] - xshake[i2][1];
s12[2] = xshake[i1][2] - xshake[i2][2];
domain->minimum_image_once(s12);
// scalar distances between atoms
@ -2502,10 +2518,10 @@ void FixShake::shake3angle(int m)
}
if (evflag) {
nlist = 0;
if (i0 < nlocal) list[nlist++] = i0;
if (i1 < nlocal) list[nlist++] = i1;
if (i2 < nlocal) list[nlist++] = i2;
int count = 0;
if (i0 < nlocal) atomlist[count++] = i0;
if (i1 < nlocal) atomlist[count++] = i1;
if (i2 < nlocal) atomlist[count++] = i2;
v[0] = lamda01*r01[0]*r01[0]+lamda02*r02[0]*r02[0]+lamda12*r12[0]*r12[0];
v[1] = lamda01*r01[1]*r01[1]+lamda02*r02[1]*r02[1]+lamda12*r12[1]*r12[1];
@ -2519,19 +2535,16 @@ void FixShake::shake3angle(int m)
{r02[0], r02[1], r02[2]},
{r12[0], r12[1], r12[2]}};
int pairlist[][2] = {{i0,i1}, {i0,i2}, {i1,i2}};
v_tally(nlist,list,3.0,v,nlocal,3,pairlist,fpairlist,dellist);
v_tally(count,atomlist,3.0,v,nlocal,3,pairlist,fpairlist,dellist);
}
}
/* ----------------------------------------------------------------------
apply bond force for minimization
apply bond force for minimization between atom indices i1 and i2
------------------------------------------------------------------------- */
void FixShake::bond_force(tagint id1, tagint id2, double length)
void FixShake::bond_force(int i1, int i2, double length)
{
int i1 = atom->map(id1);
int i2 = atom->map(id2);
if ((i1 < 0) || (i2 < 0)) return;
// distance vec between atoms, with PBC
@ -2539,7 +2552,6 @@ void FixShake::bond_force(tagint id1, tagint id2, double length)
double delx = x[i1][0] - x[i2][0];
double dely = x[i1][1] - x[i2][1];
double delz = x[i1][2] - x[i2][2];
domain->minimum_image(delx, dely, delz);
// compute and apply force
@ -2548,21 +2560,21 @@ void FixShake::bond_force(tagint id1, tagint id2, double length)
const double rk = kbond * dr;
const double fbond = (r > 0.0) ? -2.0 * rk / r : 0.0;
const double eb = rk*dr;
int list[2];
int nlist = 0;
int atomlist[2];
int count = 0;
if (i1 < nlocal) {
f[i1][0] += delx * fbond;
f[i1][1] += dely * fbond;
f[i1][2] += delz * fbond;
list[nlist++] = i1;
atomlist[count++] = i1;
ebond += 0.5*eb;
}
if (i2 < nlocal) {
f[i2][0] -= delx * fbond;
f[i2][1] -= dely * fbond;
f[i2][2] -= delz * fbond;
list[nlist++] = i2;
atomlist[count++] = i2;
ebond += 0.5*eb;
}
if (evflag) {
@ -2573,7 +2585,7 @@ void FixShake::bond_force(tagint id1, tagint id2, double length)
v[3] = 0.5 * delx * dely * fbond;
v[4] = 0.5 * delx * delz * fbond;
v[5] = 0.5 * dely * delz * fbond;
ev_tally(nlist, list, 2.0, eb, v);
ev_tally(count, atomlist, 2.0, eb, v);
}
}
@ -2616,14 +2628,13 @@ void FixShake::stats()
// bond stats
if (n == 1) n = 3;
int iatom = atom->map(shake_atom[i][0]);
int iatom = closest_list[ii][0];
for (int j = 1; j < n; j++) {
int jatom = atom->map(shake_atom[i][j]);
int jatom = closest_list[ii][j];
if (jatom >= nlocal) continue;
delx = x[iatom][0] - x[jatom][0];
dely = x[iatom][1] - x[jatom][1];
delz = x[iatom][2] - x[jatom][2];
domain->minimum_image(delx,dely,delz);
r = sqrt(delx*delx + dely*dely + delz*delz);
int m = shake_type[i][j-1];
@ -2636,9 +2647,9 @@ void FixShake::stats()
// angle stats
if (shake_flag[i] == 1) {
int iatom = atom->map(shake_atom[i][0]);
int jatom = atom->map(shake_atom[i][1]);
int katom = atom->map(shake_atom[i][2]);
int iatom = closest_list[ii][0];
int jatom = closest_list[ii][1];
int katom = closest_list[ii][2];
int n = 0;
if (iatom < nlocal) ++n;
if (jatom < nlocal) ++n;
@ -2647,19 +2658,16 @@ void FixShake::stats()
delx = x[iatom][0] - x[jatom][0];
dely = x[iatom][1] - x[jatom][1];
delz = x[iatom][2] - x[jatom][2];
domain->minimum_image(delx,dely,delz);
r1 = sqrt(delx*delx + dely*dely + delz*delz);
delx = x[iatom][0] - x[katom][0];
dely = x[iatom][1] - x[katom][1];
delz = x[iatom][2] - x[katom][2];
domain->minimum_image(delx,dely,delz);
r2 = sqrt(delx*delx + dely*dely + delz*delz);
delx = x[jatom][0] - x[katom][0];
dely = x[jatom][1] - x[katom][1];
delz = x[jatom][2] - x[katom][2];
domain->minimum_image(delx,dely,delz);
r3 = sqrt(delx*delx + dely*dely + delz*delz);
angle = acos((r1*r1 + r2*r2 - r3*r3) / (2.0*r1*r2));
@ -3258,8 +3266,6 @@ void FixShake::correct_coordinates(int vflag) {
double **xtmp = xshake;
xshake = x;
if (comm->nprocs > 1) {
comm->forward_comm(this);
}
comm->forward_comm(this);
xshake = xtmp;
}

5
src/RIGID/fix_shake.h
View File

@ -113,6 +113,7 @@ class FixShake : public Fix {
double dtf_inner, dtf_innerhalf; // timesteps for rRESPA trial move
int *list; // list of clusters to SHAKE
int **closest_list; // list of closest atom indices in SHAKE clusters
int nlist, maxlist; // size and max-size of list
// stat quantities
@ -140,8 +141,8 @@ class FixShake : public Fix {
void shake3(int);
void shake4(int);
void shake3angle(int);
void bond_force(tagint, tagint, double);
void stats();
void bond_force(int, int, double);
virtual void stats();
int bondtype_findset(int, tagint, tagint, int);
int angletype_findset(int, tagint, tagint, int);

68
src/domain.cpp
View File

@ -974,10 +974,6 @@ void Domain::subbox_too_small_check(double thresh)
changed "if" to "while" to enable distance to
far-away ghost atom returned by atom->map() to be wrapped back into box
could be problem for looking up atom IDs when cutoff > boxsize
this should not be used if atom has moved infinitely far outside box
b/c while could iterate forever
e.g. fix shake prediction of new position with highly overlapped atoms
uses minimum_image_once() instead
------------------------------------------------------------------------- */
static constexpr double MAXIMGCOUNT = 16;
@ -1050,70 +1046,6 @@ void Domain::minimum_image(double &dx, double &dy, double &dz) const
}
}
/* ----------------------------------------------------------------------
minimum image convention in periodic dimensions
use 1/2 of box size as test
for triclinic, also add/subtract tilt factors in other dims as needed
only shift by one box length in each direction
this should not be used if multiple box shifts are required
------------------------------------------------------------------------- */
void Domain::minimum_image_once(double *delta) const
{
if (triclinic == 0) {
if (xperiodic) {
if (fabs(delta[0]) > xprd_half) {
if (delta[0] < 0.0) delta[0] += xprd;
else delta[0] -= xprd;
}
}
if (yperiodic) {
if (fabs(delta[1]) > yprd_half) {
if (delta[1] < 0.0) delta[1] += yprd;
else delta[1] -= yprd;
}
}
if (zperiodic) {
if (fabs(delta[2]) > zprd_half) {
if (delta[2] < 0.0) delta[2] += zprd;
else delta[2] -= zprd;
}
}
} else {
if (zperiodic) {
if (fabs(delta[2]) > zprd_half) {
if (delta[2] < 0.0) {
delta[2] += zprd;
delta[1] += yz;
delta[0] += xz;
} else {
delta[2] -= zprd;
delta[1] -= yz;
delta[0] -= xz;
}
}
}
if (yperiodic) {
if (fabs(delta[1]) > yprd_half) {
if (delta[1] < 0.0) {
delta[1] += yprd;
delta[0] += xy;
} else {
delta[1] -= yprd;
delta[0] -= xy;
}
}
}
if (xperiodic) {
if (fabs(delta[0]) > xprd_half) {
if (delta[0] < 0.0) delta[0] += xprd;
else delta[0] -= xprd;
}
}
}
}
/* ----------------------------------------------------------------------
return local index of atom J or any of its images that is closest to atom I
if J is not a valid index like -1, just return it

1
src/domain.h
View File

@ -119,7 +119,6 @@ class Domain : protected Pointers {
void subbox_too_small_check(double);
void minimum_image(double &, double &, double &) const;
void minimum_image(double *delta) const { minimum_image(delta[0], delta[1], delta[2]); }
void minimum_image_once(double *) const;
int closest_image(int, int);
int closest_image(const double *const, int);
void closest_image(const double *const, const double *const, double *const);

21
src/modify.cpp
View File

@ -1085,7 +1085,7 @@ int Modify::find_fix(const std::string &id)
{
if (id.empty()) return -1;
for (int ifix = 0; ifix < nfix; ifix++)
if (id == fix[ifix]->id) return ifix;
if (fix[ifix] && (id == fix[ifix]->id)) return ifix;
return -1;
}
@ -1098,7 +1098,7 @@ Fix *Modify::get_fix_by_id(const std::string &id) const
{
if (id.empty()) return nullptr;
for (int ifix = 0; ifix < nfix; ifix++)
if (id == fix[ifix]->id) return fix[ifix];
if (fix[ifix] && (id == fix[ifix]->id)) return fix[ifix];
return nullptr;
}
@ -1112,9 +1112,9 @@ const std::vector<Fix *> Modify::get_fix_by_style(const std::string &style) cons
std::vector<Fix *> matches;
if (style.empty()) return matches;
for (int ifix = 0; ifix < nfix; ifix++)
if (utils::strmatch(fix[ifix]->style, style)) matches.push_back(fix[ifix]);
for (int ifix = 0; ifix < nfix; ifix++) {
if (fix[ifix] && utils::strmatch(fix[ifix]->style, style)) matches.push_back(fix[ifix]);
}
return matches;
}
@ -1349,7 +1349,7 @@ int Modify::find_compute(const std::string &id)
{
if (id.empty()) return -1;
for (int icompute = 0; icompute < ncompute; icompute++)
if (id == compute[icompute]->id) return icompute;
if (compute[icompute] && (id == compute[icompute]->id)) return icompute;
return -1;
}
@ -1362,7 +1362,7 @@ Compute *Modify::get_compute_by_id(const std::string &id) const
{
if (id.empty()) return nullptr;
for (int icompute = 0; icompute < ncompute; icompute++)
if (id == compute[icompute]->id) return compute[icompute];
if (compute[icompute] && (id == compute[icompute]->id)) return compute[icompute];
return nullptr;
}
@ -1376,9 +1376,10 @@ const std::vector<Compute *> Modify::get_compute_by_style(const std::string &sty
std::vector<Compute *> matches;
if (style.empty()) return matches;
for (int icompute = 0; icompute < ncompute; icompute++)
if (utils::strmatch(compute[icompute]->style, style)) matches.push_back(compute[icompute]);
for (int icompute = 0; icompute < ncompute; icompute++) {
if (compute[icompute] && utils::strmatch(compute[icompute]->style, style))
matches.push_back(compute[icompute]);
}
return matches;
}

8
src/read_data.cpp
View File

@ -1154,6 +1154,14 @@ void ReadData::header(int firstpass)
if (me == 0) {
char *eof = utils::fgets_trunc(line, MAXLINE, fp);
if (eof == nullptr) error->one(FLERR, "Unexpected end of data file");
// check for units keyword in first line and print warning on mismatch
auto units = Tokenizer(utils::strfind(line, "units = \\w+")).as_vector();
if (units.size() > 2) {
if (units[2] != update->unit_style)
error->warning(FLERR, "Inconsistent units in data file: current = {}, data file = {}",
update->unit_style, units[2]);
}
}
while (true) {

2
src/version.h
View File

@ -1,2 +1,2 @@
#define LAMMPS_VERSION "28 Mar 2023"
#define LAMMPS_UPDATE "Update 1"
#define LAMMPS_UPDATE "Development"

4
src/write_data.cpp
View File

@ -245,8 +245,8 @@ void WriteData::write(const std::string &file)
void WriteData::header()
{
fmt::print(fp,"LAMMPS data file via write_data, version {}, "
"timestep = {}\n\n",lmp->version,update->ntimestep);
fmt::print(fp,"LAMMPS data file via write_data, version {}, timestep = {}, units = {}\n\n",
lmp->version, update->ntimestep, update->unit_style);
fmt::print(fp,"{} atoms\n{} atom types\n",atom->natoms,atom->ntypes);

7
tools/singularity/fedora35_mingw.def
View File

@ -35,12 +35,17 @@ From: fedora:35
texlive-latex-fonts texlive-pslatex texlive-collection-latexrecommended \
texlive-latex texlive-latexconfig doxygen-latex texlive-collection-latex \
texlive-latex-bin texlive-lualatex-math texlive-fncychap texlive-tabulary \
texlive-framed texlive-wrapfig texlive-upquote texlive-capt-of \
texlive-framed texlive-wrapfig texlive-upquote texlive-capt-of texlive-pict2e \
texlive-needspace texlive-titlesec texlive-anysize texlive-dvipng texlive-xindy \
blas-devel lapack-devel libyaml-devel openkim-models kim-api-devel \
zstd libzstd-devel yaml-cpp-devel
dnf clean all
# install NSIS EnVar plugin
curl -L -o EnVar_plugin.zip https://nsis.sourceforge.io/mediawiki/images/7/7f/EnVar_plugin.zip
unzip -d /usr/share/nsis EnVar_plugin.zip
rm EnVar_plugin.zip
# enable Lmod and load MPI
source /usr/share/lmod/lmod/init/profile
module purge

7
tools/singularity/fedora36_mingw.def
View File

@ -35,12 +35,17 @@ From: fedora:36
texlive-latex-fonts texlive-pslatex texlive-collection-latexrecommended \
texlive-latex texlive-latexconfig doxygen-latex texlive-collection-latex \
texlive-latex-bin texlive-lualatex-math texlive-fncychap texlive-tabulary \
texlive-framed texlive-wrapfig texlive-upquote texlive-capt-of \
texlive-framed texlive-wrapfig texlive-upquote texlive-capt-of texlive-pict2e \
texlive-needspace texlive-titlesec texlive-anysize texlive-dvipng texlive-xindy \
blas-devel lapack-devel libyaml-devel openkim-models kim-api-devel \
zstd libzstd-devel yaml-cpp-devel
dnf clean all
# install NSIS EnVar plugin
curl -L -o EnVar_plugin.zip https://nsis.sourceforge.io/mediawiki/images/7/7f/EnVar_plugin.zip
unzip -d /usr/share/nsis EnVar_plugin.zip
rm EnVar_plugin.zip
# enable Lmod and load MPI
source /usr/share/lmod/lmod/init/profile
module purge

2
unittest/CMakeLists.txt
View File

@ -54,7 +54,7 @@ add_test(NAME RunLammps
COMMAND $<TARGET_FILE:lmp> -log none -echo none -in in.empty)
set_tests_properties(RunLammps PROPERTIES
ENVIRONMENT "TSAN_OPTIONS=ignore_noninstrumented_modules=1;HWLOC_HIDE_ERRORS=2"
PASS_REGULAR_EXPRESSION "LAMMPS \\([0-9]+ [A-Za-z]+ 2[0-9][0-9][0-9]( - Update [0-9]+)?\\)")
PASS_REGULAR_EXPRESSION "LAMMPS \\([0-9]+ [A-Za-z]+ 2[0-9][0-9][0-9]( - Update [0-9]+)?( - Development)?( - Maintenance)?\\)")
# check if the compiled executable will print the help message
add_test(NAME HelpMessage

6
unittest/c-library/CMakeLists.txt
View File

@ -80,7 +80,7 @@ if(BUILD_MPI)
COMMAND $<TARGET_FILE:simple-plugin> 1 ${LAMMPS_DIR}/examples/COUPLE/plugin/in.lj $<TARGET_FILE:lammps>)
set_tests_properties(RunCoupleSimplePlugin PROPERTIES
ENVIRONMENT "TSAN_OPTIONS=ignore_noninstrumented_modules=1;HWLOC_HIDE_ERRORS=2"
PASS_REGULAR_EXPRESSION "LAMMPS \\([0-9]+ [A-Za-z]+ 2[0-9][0-9][0-9]( - Update [0-9]+)?\\)")
PASS_REGULAR_EXPRESSION "LAMMPS \\([0-9]+ [A-Za-z]+ 2[0-9][0-9][0-9]( - Update [0-9]+)?( - Development)?( - Maintenance)?\\)")
endif()
add_subdirectory(${LAMMPS_DIR}/examples/COUPLE/simple ${CMAKE_BINARY_DIR}/build-simple)
add_test(NAME RunCoupleSimpleC
@ -93,10 +93,10 @@ if(BUILD_MPI)
COMMAND $<TARGET_FILE:simpleF90> 1 ${LAMMPS_DIR}/examples/COUPLE/simple/in.lj)
set_tests_properties(RunCoupleSimpleF90 PROPERTIES
ENVIRONMENT "TSAN_OPTIONS=ignore_noninstrumented_modules=1;HWLOC_HIDE_ERRORS=2"
PASS_REGULAR_EXPRESSION "LAMMPS \\([0-9]+ [A-Za-z]+ 2[0-9][0-9][0-9]( - Update [0-9]+)?\\)")
PASS_REGULAR_EXPRESSION "LAMMPS \\([0-9]+ [A-Za-z]+ 2[0-9][0-9][0-9]( - Update [0-9]+)?( - Development)?( - Maintenance)?\\)")
endif()
set_tests_properties(RunCoupleSimpleC RunCoupleSimpleCC PROPERTIES
ENVIRONMENT "TSAN_OPTIONS=ignore_noninstrumented_modules=1;HWLOC_HIDE_ERRORS=2"
PASS_REGULAR_EXPRESSION "LAMMPS \\([0-9]+ [A-Za-z]+ 2[0-9][0-9][0-9]( - Update [0-9]+)?\\)")
PASS_REGULAR_EXPRESSION "LAMMPS \\([0-9]+ [A-Za-z]+ 2[0-9][0-9][0-9]( - Update [0-9]+)?( - Development)?( - Maintenance)?\\)")
endif()