tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into nm_split_styles

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-10-08 15:57:17 -04:00
parent be3468ae07 ab51c1bd3d
commit fc42992cdf
230 changed files with 2545 additions and 2032 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
cmake/CMakeLists.txt
View File

@ -90,8 +90,11 @@ if((CMAKE_CXX_COMPILER_ID STREQUAL "Intel") OR (CMAKE_CXX_COMPILER_ID STREQUAL "
endif()
endif()
# we require C++11 without extensions
# we require C++11 without extensions. Kokkos requires at least C++14 (currently)
set(CMAKE_CXX_STANDARD 11)
if(PKG_KOKKOS AND (CMAKE_CXX_STANDARD LESS 14))
set(CMAKE_CXX_STANDARD 14)
endif()
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF CACHE BOOL "Use compiler extensions")
# ugly hack for MSVC which by default always reports an old C++ standard in the __cplusplus macro

4
cmake/Modules/Packages/KOKKOS.cmake
View File

@ -1,6 +1,8 @@
########################################################################
# As of version 3.3.0 Kokkos requires C++14
set(CMAKE_CXX_STANDARD 14)
if(CMAKE_CXX_STANDARD LESS 14)
message(FATAL_ERROR "The KOKKOS package requires the C++ standard to be set to at least C++14")
endif()
########################################################################
# consistency checks and Kokkos options/settings required by LAMMPS
if(Kokkos_ENABLE_CUDA)

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

@ -19,6 +19,14 @@ if(DOWNLOAD_LATTE)
set(LATTE_MD5 "820e73a457ced178c08c71389a385de7" CACHE STRING "MD5 checksum of LATTE tarball")
mark_as_advanced(LATTE_URL)
mark_as_advanced(LATTE_MD5)
# 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))
message(FATAL_ERROR "Cannot compile downloaded LATTE library due to a technical limitation")
endif()
include(ExternalProject)
ExternalProject_Add(latte_build
URL ${LATTE_URL}

5
cmake/Modules/Packages/MACHDYN.cmake
View File

@ -7,8 +7,9 @@ endif()
option(DOWNLOAD_EIGEN3 "Download Eigen3 instead of using an already installed one)" ${DOWNLOAD_EIGEN3_DEFAULT})
if(DOWNLOAD_EIGEN3)
message(STATUS "Eigen3 download requested - we will build our own")
set(EIGEN3_URL "https://gitlab.com/libeigen/eigen/-/archive/3.3.9/eigen-3.3.9.tar.gz" CACHE STRING "URL for Eigen3 tarball")
set(EIGEN3_MD5 "609286804b0f79be622ccf7f9ff2b660" CACHE STRING "MD5 checksum of Eigen3 tarball")
set(EIGEN3_URL "https://download.lammps.org/thirdparty/eigen-3.4.0.tar.gz" CACHE STRING "URL for Eigen3 tarball")
set(EIGEN3_MD5 "4c527a9171d71a72a9d4186e65bea559" CACHE STRING "MD5 checksum of Eigen3 tarball")
mark_as_advanced(EIGEN3_URL)
mark_as_advanced(EIGEN3_MD5)
include(ExternalProject)

10
cmake/Modules/Packages/ML-HDNNP.cmake
View File

@ -45,12 +45,12 @@ if(DOWNLOAD_N2P2)
# get path to MPI include directory when cross-compiling to windows
if((CMAKE_SYSTEM_NAME STREQUAL Windows) AND CMAKE_CROSSCOMPILING)
get_target_property(N2P2_MPI_INCLUDE MPI::MPI_CXX INTERFACE_INCLUDE_DIRECTORIES)
set(N2P2_PROJECT_OPTIONS "-I ${N2P2_MPI_INCLUDE} -DMPICH_SKIP_MPICXX=1")
set(N2P2_PROJECT_OPTIONS "-I${N2P2_MPI_INCLUDE}")
set(MPI_CXX_COMPILER ${CMAKE_CXX_COMPILER})
endif()
if(CMAKE_CXX_COMPILER_ID STREQUAL "Intel")
get_target_property(N2P2_MPI_INCLUDE MPI::MPI_CXX INTERFACE_INCLUDE_DIRECTORIES)
set(N2P2_PROJECT_OPTIONS "-I ${N2P2_MPI_INCLUDE} -DMPICH_SKIP_MPICXX=1")
set(N2P2_PROJECT_OPTIONS "-I${N2P2_MPI_INCLUDE}")
set(MPI_CXX_COMPILER ${CMAKE_CXX_COMPILER})
endif()
endif()
@ -69,6 +69,12 @@ if(DOWNLOAD_N2P2)
# echo final flag for debugging
message(STATUS "N2P2 BUILD OPTIONS: ${N2P2_BUILD_OPTIONS}")
# must have "sed" command to compile n2p2 library (for now)
find_program(HAVE_SED sed)
if(NOT HAVE_SED)
message(FATAL_ERROR "Must have 'sed' program installed to compile 'n2p2' library for ML-HDNNP package")
endif()
# download compile n2p2 library. much patch MPI calls in LAMMPS interface to accommodate MPI-2 (e.g. for cross-compiling)
include(ExternalProject)
ExternalProject_Add(n2p2_build

2
cmake/Modules/Packages/ML-QUIP.cmake
View File

@ -50,7 +50,7 @@ if(DOWNLOAD_QUIP)
GIT_TAG origin/public
GIT_SHALLOW YES
GIT_PROGRESS YES
PATCH_COMMAND cp ${CMAKE_BINARY_DIR}/quip.config <SOURCE_DIR>/arch/Makefile.lammps
PATCH_COMMAND ${CMAKE_COMMAND} -E copy_if_different ${CMAKE_BINARY_DIR}/quip.config <SOURCE_DIR>/arch/Makefile.lammps
CONFIGURE_COMMAND env QUIP_ARCH=lammps make config
BUILD_COMMAND env QUIP_ARCH=lammps make libquip
INSTALL_COMMAND ""

7
cmake/Modules/Packages/MSCG.cmake
View File

@ -12,6 +12,13 @@ if(DOWNLOAD_MSCG)
mark_as_advanced(MSCG_URL)
mark_as_advanced(MSCG_MD5)
# 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))
message(FATAL_ERROR "Cannot compile downloaded MSCG library due to a technical limitation")
endif()
include(ExternalProject)
ExternalProject_Add(mscg_build
URL ${MSCG_URL}

5
cmake/Modules/Packages/SCAFACOS.cmake
View File

@ -23,6 +23,11 @@ if(DOWNLOAD_SCAFACOS)
file(DOWNLOAD ${LAMMPS_THIRDPARTY_URL}/scafacos-1.0.1-fix.diff ${CMAKE_CURRENT_BINARY_DIR}/scafacos-1.0.1.fix.diff
EXPECTED_HASH MD5=4baa1333bb28fcce102d505e1992d032)
find_program(HAVE_PATCH patch)
if(NOT HAVE_PATCH)
message(FATAL_ERROR "The 'patch' program is required to build the ScaFaCoS library")
endif()
include(ExternalProject)
ExternalProject_Add(scafacos_build
URL ${SCAFACOS_URL}

5
cmake/Modules/Packages/VORONOI.cmake
View File

@ -26,6 +26,11 @@ if(DOWNLOAD_VORO)
set(VORO_BUILD_OPTIONS CXX=${CMAKE_CXX_COMPILER} CFLAGS=${VORO_BUILD_CFLAGS})
endif()
find_program(HAVE_PATCH patch)
if(NOT HAVE_PATCH)
message(FATAL_ERROR "The 'patch' program is required to build the voro++ library")
endif()
ExternalProject_Add(voro_build
URL ${VORO_URL}
URL_MD5 ${VORO_MD5}

2
doc/lammps.1
View File

@ -1,4 +1,4 @@
.TH LAMMPS "20 September 2021" "2021-09-20"
.TH LAMMPS "29 September 2021" "2021-09-29"
.SH NAME
.B LAMMPS
\- Molecular Dynamics Simulator.

3
doc/src/Build_settings.rst
View File

@ -71,7 +71,8 @@ LAMMPS can use them if they are available on your system.
-D FFTW3_INCLUDE_DIR=path # path to FFTW3 include files
-D FFTW3_LIBRARY=path # path to FFTW3 libraries
-D FFT_FFTW_THREADS=on # enable using threaded FFTW3 libraries
-D FFTW3_OMP_LIBRARY=path # path to FFTW3 OpenMP wrapper libraries
-D FFT_FFTW_THREADS=on # enable using OpenMP threaded FFTW3 libraries
-D MKL_INCLUDE_DIR=path # ditto for Intel MKL library
-D FFT_MKL_THREADS=on # enable using threaded FFTs with MKL libraries
-D MKL_LIBRARY=path # path to MKL libraries

1
doc/src/Developer.rst
View File

@ -11,6 +11,7 @@ of time and requests from the LAMMPS user community.
:maxdepth: 1
Developer_org
Developer_parallel
Developer_flow
Developer_write
Developer_notes

120
doc/src/Developer_par_comm.rst Normal file
View File

@ -0,0 +1,120 @@
Communication
^^^^^^^^^^^^^
Following the partitioning scheme in use all per-atom data is
distributed across the MPI processes, which allows LAMMPS to handle very
large systems provided it uses a correspondingly large number of MPI
processes. Since The per-atom data (atom IDs, positions, velocities,
types, etc.) To be able to compute the short-range interactions MPI
processes need not only access to data of atoms they "own" but also
information about atoms from neighboring sub-domains, in LAMMPS referred
to as "ghost" atoms. These are copies of atoms storing required
per-atom data for up to the communication cutoff distance. The green
dashed-line boxes in the :ref:`domain-decomposition` figure illustrate
the extended ghost-atom sub-domain for one processor.
This approach is also used to implement periodic boundary
conditions: atoms that lie within the cutoff distance across a periodic
boundary are also stored as ghost atoms and taken from the periodic
replication of the sub-domain, which may be the same sub-domain, e.g. if
running in serial. As a consequence of this, force computation in
LAMMPS is not subject to minimum image conventions and thus cutoffs may
be larger than half the simulation domain.
.. _ghost-atom-comm:
.. figure:: img/ghost-comm.png
:align: center
ghost atom communication
This figure shows the ghost atom communication patterns between
sub-domains for "brick" (left) and "tiled" communication styles for
2d simulations. The numbers indicate MPI process ranks. Here the
sub-domains are drawn spatially separated for clarity. The
dashed-line box is the extended sub-domain of processor 0 which
includes its ghost atoms. The red- and blue-shaded boxes are the
regions of communicated ghost atoms.
Efficient communication patterns are needed to update the "ghost" atom
data, since that needs to be done at every MD time step or minimization
step. The diagrams of the `ghost-atom-comm` figure illustrate how ghost
atom communication is performed in two stages for a 2d simulation (three
in 3d) for both a regular and irregular partitioning of the simulation
box. For the regular case (left) atoms are exchanged first in the
*x*-direction, then in *y*, with four neighbors in the grid of processor
sub-domains.
In the *x* stage, processor ranks 1 and 2 send owned atoms in their
red-shaded regions to rank 0 (and vice versa). Then in the *y* stage,
ranks 3 and 4 send atoms in their blue-shaded regions to rank 0, which
includes ghost atoms they received in the *x* stage. Rank 0 thus
acquires all its ghost atoms; atoms in the solid blue corner regions
are communicated twice before rank 0 receives them.
For the irregular case (right) the two stages are similar, but a
processor can have more than one neighbor in each direction. In the
*x* stage, MPI ranks 1,2,3 send owned atoms in their red-shaded regions to
rank 0 (and vice versa). These include only atoms between the lower
and upper *y*-boundary of rank 0's sub-domain. In the *y* stage, ranks
4,5,6 send atoms in their blue-shaded regions to rank 0. This may
include ghost atoms they received in the *x* stage, but only if they
are needed by rank 0 to fill its extended ghost atom regions in the
+/-*y* directions (blue rectangles). Thus in this case, ranks 5 and
6 do not include ghost atoms they received from each other (in the *x*
stage) in the atoms they send to rank 0. The key point is that while
the pattern of communication is more complex in the irregular
partitioning case, it can still proceed in two stages (three in 3d)
via atom exchanges with only neighboring processors.
When attributes of owned atoms are sent to neighboring processors to
become attributes of their ghost atoms, LAMMPS calls this a "forward"
communication. On timesteps when atoms migrate to new owning processors
and neighbor lists are rebuilt, each processor creates a list of its
owned atoms which are ghost atoms in each of its neighbor processors.
These lists are used to pack per-atom coordinates (for example) into
message buffers in subsequent steps until the next reneighboring.
A "reverse" communication is when computed ghost atom attributes are
sent back to the processor who owns the atom. This is used (for
example) to sum partial forces on ghost atoms to the complete force on
owned atoms. The order of the two stages described in the
:ref:`ghost-atom-comm` figure is inverted and the same lists of atoms
are used to pack and unpack message buffers with per-atom forces. When
a received buffer is unpacked, the ghost forces are summed to owned atom
forces. As in forward communication, forces on atoms in the four blue
corners of the diagrams are sent, received, and summed twice (once at
each stage) before owning processors have the full force.
These two operations are used many places within LAMMPS aside from
exchange of coordinates and forces, for example by manybody potentials
to share intermediate per-atom values, or by rigid-body integrators to
enable each atom in a body to access body properties. Here are
additional details about how these communication operations are
performed in LAMMPS:
- When exchanging data with different processors, forward and reverse
communication is done using ``MPI_Send()`` and ``MPI_IRecv()`` calls.
If a processor is "exchanging" atoms with itself, only the pack and
unpack operations are performed, e.g. to create ghost atoms across
periodic boundaries when running on a single processor.
- For forward communication of owned atom coordinates, periodic box
lengths are added and subtracted when the receiving processor is
across a periodic boundary from the sender. There is then no need to
apply a minimum image convention when calculating distances between
atom pairs when building neighbor lists or computing forces.
- The cutoff distance for exchanging ghost atoms is typically equal to
the neighbor cutoff. But it can also chosen to be longer if needed,
e.g. half the diameter of a rigid body composed of multiple atoms or
over 3x the length of a stretched bond for dihedral interactions. It
can also exceed the periodic box size. For the regular communication
pattern (left), if the cutoff distance extends beyond a neighbor
processor's sub-domain, then multiple exchanges are performed in the
same direction. Each exchange is with the same neighbor processor,
but buffers are packed/unpacked using a different list of atoms. For
forward communication, in the first exchange a processor sends only
owned atoms. In subsequent exchanges, it sends ghost atoms received
in previous exchanges. For the irregular pattern (right) overlaps of
a processor's extended ghost-atom sub-domain with all other processors
in each dimension are detected.

188
doc/src/Developer_par_long.rst Normal file
View File

@ -0,0 +1,188 @@
Long-range interactions
^^^^^^^^^^^^^^^^^^^^^^^
For charged systems, LAMMPS can compute long-range Coulombic
interactions via the FFT-based particle-particle/particle-mesh (PPPM)
method implemented in :doc:`kspace style pppm and its variants
<kspace_style>`. For that Coulombic interactions are partitioned into
short- and long-range components. The short-ranged portion is computed
in real space as a loop over pairs of charges within a cutoff distance,
using neighbor lists. The long-range portion is computed in reciprocal
space using a kspace style. For the PPPM implementation the simulation
cell is overlaid with a regular FFT grid in 3d. It proceeds in several stages:
a) each atom's point charge is interpolated to nearby FFT grid points,
b) a forward 3d FFT is performed,
c) a convolution operation is performed in reciprocal space,
d) one or more inverse 3d FFTs are performed, and
e) electric field values from grid points near each atom are interpolated to compute
its forces.
For any of the spatial-decomposition partitioning schemes each processor
owns the brick-shaped portion of FFT grid points contained within its
sub-domain. The two interpolation operations use a stencil of grid
points surrounding each atom. To accommodate the stencil size, each
processor also stores a few layers of ghost grid points surrounding its
brick. Forward and reverse communication of grid point values is
performed similar to the corresponding :doc:`atom data communication
<Developer_par_comm>`. In this case, electric field values on owned
grid points are sent to neighboring processors to become ghost point
values. Likewise charge values on ghost points are sent and summed to
values on owned points.
For triclinic simulation boxes, the FFT grid planes are parallel to
the box faces, but the mapping of charge and electric field values
to/from grid points is done in reduced coordinates where the tilted
box is conceptually a unit cube, so that the stencil and FFT
operations are unchanged. However the FFT grid size required for a
given accuracy is larger for triclinic domains than it is for
orthogonal boxes.
.. _fft-parallel:
.. figure:: img/fft-decomp-parallel.png
:align: center
parallel FFT in PPPM
Stages of a parallel FFT for a simulation domain overlaid
with an 8x8x8 3d FFT grid, partitioned across 64 processors.
Within each of the 4 diagrams, grid cells of the same color are
owned by a single processor; for simplicity only cells owned by 4
or 8 of the 64 processors are colored. The two images on the left
illustrate brick-to-pencil communication. The two images on the
right illustrate pencil-to-pencil communication, which in this
case transposes the *y* and *z* dimensions of the grid.
Parallel 3d FFTs require substantial communication relative to their
computational cost. A 3d FFT is implemented by a series of 1d FFTs
along the *x-*, *y-*, and *z-*\ direction of the FFT grid. Thus the FFT
grid cannot be decomposed like atoms into 3 dimensions for parallel
processing of the FFTs but only in 1 (as planes) or 2 (as pencils)
dimensions and in between the steps the grid needs to be transposed to
have the FFT grid portion "owned" by each MPI process complete in the
direction of the 1d FFTs it has to perform. LAMMPS uses the
pencil-decomposition algorithm as shown in the :ref:`fft-parallel` figure.
Initially (far left), each processor owns a brick of same-color grid
cells (actually grid points) contained within in its sub-domain. A
brick-to-pencil communication operation converts this layout to 1d
pencils in the *x*-dimension (center left). Again, cells of the same
color are owned by the same processor. Each processor can then compute
a 1d FFT on each pencil of data it wholly owns using a call to the
configured FFT library. A pencil-to-pencil communication then converts
this layout to pencils in the *y* dimension (center right) which
effectively transposes the *x* and *y* dimensions of the grid, followed
by 1d FFTs in *y*. A final transpose of pencils from *y* to *z* (far
right) followed by 1d FFTs in *z* completes the forward FFT. The data
is left in a *z*-pencil layout for the convolution operation. One or
more inverse FFTs then perform the sequence of 1d FFTs and communication
steps in reverse order; the final layout of resulting grid values is the
same as the initial brick layout.
Each communication operation within the FFT (brick-to-pencil or
pencil-to-pencil or pencil-to-brick) converts one tiling of the 3d grid
to another, where a tiling in this context means an assignment of a
small brick-shaped subset of grid points to each processor, the union of
which comprise the entire grid. The parallel `fftMPI library
<https://lammps.github.io/fftmpi/>`_ written for LAMMPS allows arbitrary
definitions of the tiling so that an irregular partitioning of the
simulation domain can use it directly. Transforming data from one
tiling to another is implemented in `fftMPI` using point-to-point
communication, where each processor sends data to a few other
processors, since each tile in the initial tiling overlaps with a
handful of tiles in the final tiling.
The transformations could also be done using collective communication
across all $P$ processors with a single call to ``MPI_Alltoall()``, but
this is typically much slower. However, for the specialized brick and
pencil tiling illustrated in :ref:`fft-parallel` figure, collective
communication across the entire MPI communicator is not required. In
the example an :math:`8^3` grid with 512 grid cells is partitioned
across 64 processors; each processor owns a 2x2x2 3d brick of grid
cells. The initial brick-to-pencil communication (upper left to upper
right) only requires collective communication within subgroups of 4
processors, as illustrated by the 4 colors. More generally, a
brick-to-pencil communication can be performed by partitioning *P*
processors into :math:`P^{\frac{2}{3}}` subgroups of
:math:`P^{\frac{1}{3}}` processors each. Each subgroup performs
collective communication only within its subgroup. Similarly,
pencil-to-pencil communication can be performed by partitioning *P*
processors into :math:`P^{\frac{1}{2}}` subgroups of
:math:`P^{\frac{1}{2}}` processors each. This is illustrated in the
figure for the :math:`y \Rightarrow z` communication (center). An
eight-processor subgroup owns the front *yz* plane of data and performs
collective communication within the subgroup to transpose from a
*y*-pencil to *z*-pencil layout.
LAMMPS invokes point-to-point communication by default, but also
provides the option of partitioned collective communication when using a
:doc:`kspace_modify collective yes <kspace_modify>` command to switch to
that mode. In the latter case, the code detects the size of the
disjoint subgroups and partitions the single *P*-size communicator into
multiple smaller communicators, each of which invokes collective
communication. Testing on a large IBM Blue Gene/Q machine at Argonne
National Labs showed a significant improvement in FFT performance for
large processor counts; partitioned collective communication was faster
than point-to-point communication or global collective communication
involving all *P* processors.
Here are some additional details about FFTs for long-range and related
grid/particle operations that LAMMPS supports:
- The fftMPI library allows each grid dimension to be a multiple of
small prime factors (2,3,5), and allows any number of processors to
perform the FFT. The resulting brick and pencil decompositions are
thus not always as well-aligned but the size of subgroups of
processors for the two modes of communication (brick/pencil and
pencil/pencil) still scale as :math:`O(P^{\frac{1}{3}})` and
:math:`O(P^{\frac{1}{2}})`.
- For efficiency in performing 1d FFTs, the grid transpose
operations illustrated in Figure \ref{fig:fft} also involve
reordering the 3d data so that a different dimension is contiguous
in memory. This reordering can be done during the packing or
unpacking of buffers for MPI communication.
- For large systems and particularly a large number of MPI processes,
the dominant cost for parallel FFTs is often the communication, not
the computation of 1d FFTs, even though the latter scales as :math:`N
\log(N)` in the number of grid points *N* per grid direction. This is
due to the fact that only a 2d decomposition into pencils is possible
while atom data (and their corresponding short-range force and energy
computations) can be decomposed efficiently in 3d.
This can be addressed by reducing the number of MPI processes involved
in the MPI communication by using :doc:`hybrid MPI + OpenMP
parallelization <Speed_omp>`. This will use OpenMP parallelization
inside the MPI domains and while that may have a lower parallel
efficiency, it reduces the communication overhead.
As an alternative it is also possible to start a :ref:`multi-partition
<partition>` calculation and then use the :doc:`verlet/split
integrator <run_style>` to perform the PPPM computation on a
dedicated, separate partition of MPI processes. This uses an integer
"1:*p*" mapping of *p* sub-domains of the atom decomposition to one
sub-domain of the FFT grid decomposition and where pairwise non-bonded
and bonded forces and energies are computed on the larger partition
and the PPPM kspace computation concurrently on the smaller partition.
- LAMMPS also implements PPPM-based solvers for other long-range
interactions, dipole and dispersion (Lennard-Jones), which can be used
in conjunction with long-range Coulombics for point charges.
- LAMMPS implements a ``GridComm`` class which overlays the simulation
domain with a regular grid, partitions it across processors in a
manner consistent with processor sub-domains, and provides methods for
forward and reverse communication of owned and ghost grid point
values. It is used for PPPM as an FFT grid (as outlined above) and
also for the MSM algorithm which uses a cascade of grid sizes from
fine to coarse to compute long-range Coulombic forces. The GridComm
class is also useful for models where continuum fields interact with
particles. For example, the two-temperature model (TTM) defines heat
transfer between atoms (particles) and electrons (continuum gas) where
spatial variations in the electron temperature are computed by finite
differences of a discretized heat equation on a regular grid. The
:doc:`fix ttm/grid <fix_ttm>` command uses the ``GridComm`` class
internally to perform its grid operations on a distributed grid
instead of the original :doc:`fix ttm <fix_ttm>` which uses a
replicated grid.

159
doc/src/Developer_par_neigh.rst Normal file
View File

@ -0,0 +1,159 @@
Neighbor lists
^^^^^^^^^^^^^^
To compute forces efficiently, each processor creates a Verlet-style
neighbor list which enumerates all pairs of atoms *i,j* (*i* = owned,
*j* = owned or ghost) with separation less than the applicable
neighbor list cutoff distance. In LAMMPS the neighbor lists are stored
in a multiple-page data structure; each page is a contiguous chunk of
memory which stores vectors of neighbor atoms *j* for many *i* atoms.
This allows pages to be incrementally allocated or deallocated in blocks
as needed. Neighbor lists typically consume the most memory of any data
structure in LAMMPS. The neighbor list is rebuilt (from scratch) once
every few timesteps, then used repeatedly each step for force or other
computations. The neighbor cutoff distance is :math:`R_n = R_f +
\Delta_s`, where :math:`R_f` is the (largest) force cutoff defined by
the interatomic potential for computing short-range pairwise or manybody
forces and :math:`\Delta_s` is a "skin" distance that allows the list to
be used for multiple steps assuming that atoms do not move very far
between consecutive time steps. Typically the code triggers
reneighboring when any atom has moved half the skin distance since the
last reneighboring; this and other options of the neighbor list rebuild
can be adjusted with the :doc:`neigh_modify <neigh_modify>` command.
On steps when reneighboring is performed, atoms which have moved outside
their owning processor's sub-domain are first migrated to new processors
via communication. Periodic boundary conditions are also (only)
enforced on these steps to ensure each atom is re-assigned to the
correct processor. After migration, the atoms owned by each processor
are stored in a contiguous vector. Periodically each processor
spatially sorts owned atoms within its vector to reorder it for improved
cache efficiency in force computations and neighbor list building. For
that atoms are spatially binned and then reordered so that atoms in the
same bin are adjacent in the vector. Atom sorting can be disabled or
its settings modified with the :doc:`atom_modify <atom_modify>` command.
.. _neighbor-stencil:
.. figure:: img/neigh-stencil.png
:align: center
neighbor list stencils
A 2d simulation sub-domain (thick black line) and the corresponding
ghost atom cutoff region (dashed blue line) for both orthogonal
(left) and triclinic (right) domains. A regular grid of neighbor
bins (thin lines) overlays the entire simulation domain and need not
align with sub-domain boundaries; only the portion overlapping the
augmented sub-domain is shown. In the triclinic case it overlaps the
bounding box of the tilted rectangle. The blue- and red-shaded bins
represent a stencil of bins searched to find neighbors of a particular
atom (black dot).
To build a local neighbor list in linear time, the simulation domain is
overlaid (conceptually) with a regular 3d (or 2d) grid of neighbor bins,
as shown in the :ref:`neighbor-stencil` figure for 2d models and a
single MPI processor's sub-domain. Each processor stores a set of
neighbor bins which overlap its sub-domain extended by the neighbor
cutoff distance :math:`R_n`. As illustrated, the bins need not align
with processor boundaries; an integer number in each dimension is fit to
the size of the entire simulation box.
Most often LAMMPS builds what it calls a "half" neighbor list where
each *i,j* neighbor pair is stored only once, with either atom *i* or
*j* as the central atom. The build can be done efficiently by using a
pre-computed "stencil" of bins around a central origin bin which
contains the atom whose neighbors are being searched for. A stencil
is simply a list of integer offsets in *x,y,z* of nearby bins
surrounding the origin bin which are close enough to contain any
neighbor atom *j* within a distance :math:`R_n` from any atom *i* in the
origin bin. Note that for a half neighbor list, the stencil can be
asymmetric since each atom only need store half its nearby neighbors.
These stencils are illustrated in the figure for a half list and a bin
size of :math:`\frac{1}{2} R_n`. There are 13 red+blue stencil bins in
2d (for the orthogonal case, 15 for triclinic). In 3d there would be
63, 13 in the plane of bins that contain the origin bin and 25 in each
of the two planes above it in the *z* direction (75 for triclinic). The
reason the triclinic stencil has extra bins is because the bins tile the
bounding box of the entire triclinic domain and thus are not periodic
with respect to the simulation box itself. The stencil and logic for
determining which *i,j* pairs to include in the neighbor list are
altered slightly to account for this.
To build a neighbor list, a processor first loops over its "owned" plus
"ghost" atoms and assigns each to a neighbor bin. This uses an integer
vector to create a linked list of atom indices within each bin. It then
performs a triply-nested loop over its owned atoms *i*, the stencil of
bins surrounding atom *i*'s bin, and the *j* atoms in each stencil bin
(including ghost atoms). If the distance :math:`r_{ij} < R_n`, then
atom *j* is added to the vector of atom *i*'s neighbors.
Here are additional details about neighbor list build options LAMMPS
supports:
- The choice of bin size is an option; a size half of :math:`R_n` has
been found to be optimal for many typical cases. Smaller bins incur
additional overhead to loop over; larger bins require more distance
calculations. Note that for smaller bin sizes, the 2d stencil in the
figure would be more semi-circular in shape (hemispherical in 3d),
with bins near the corners of the square eliminated due to their
distance from the origin bin.
- Depending on the interatomic potential(s) and other commands used in
an input script, multiple neighbor lists and stencils with different
attributes may be needed. This includes lists with different cutoff
distances, e.g. for force computation versus occasional diagnostic
computations such as a radial distribution function, or for the
r-RESPA time integrator which can partition pairwise forces by
distance into subsets computed at different time intervals. It
includes "full" lists (as opposed to half lists) where each *i,j* pair
appears twice, stored once with *i* and *j*, and which use a larger
symmetric stencil. It also includes lists with partial enumeration of
ghost atom neighbors. The full and ghost-atom lists are used by
various manybody interatomic potentials. Lists may also use different
criteria for inclusion of a pair interaction. Typically this simply
depends only on the distance between two atoms and the cutoff
distance. But for finite-size coarse-grained particles with
individual diameters (e.g. polydisperse granular particles), it can
also depend on the diameters of the two particles.
- When using :doc:`pair style hybrid <pair_hybrid>` multiple sub-lists
of the master neighbor list for the full system need to be generated,
one for each sub-style, which contains only the *i,j* pairs needed to
compute interactions between subsets of atoms for the corresponding
potential. This means not all *i* or *j* atoms owned by a processor
are included in a particular sub-list.
- Some models use different cutoff lengths for pairwise interactions
between different kinds of particles which are stored in a single
neighbor list. One example is a solvated colloidal system with large
colloidal particles where colloid/colloid, colloid/solvent, and
solvent/solvent interaction cutoffs can be dramatically different.
Another is a model of polydisperse finite-size granular particles;
pairs of particles interact only when they are in contact with each
other. Mixtures with particle size ratios as high as 10-100x may be
used to model realistic systems. Efficient neighbor list building
algorithms for these kinds of systems are available in LAMMPS. They
include a method which uses different stencils for different cutoff
lengths and trims the stencil to only include bins that straddle the
cutoff sphere surface. More recently a method which uses both
multiple stencils and multiple bin sizes was developed; it builds
neighbor lists efficiently for systems with particles of any size
ratio, though other considerations (timestep size, force computations)
may limit the ability to model systems with huge polydispersity.
- For small and sparse systems and as a fallback method, LAMMPS also
supports neighbor list construction without binning by using a full
:math:`O(N^2)` loop over all *i,j* atom pairs in a sub-domain when
using the :doc:`neighbor nsq <neighbor>` command.
- Dependent on the "pair" setting of the :doc:`newton <newton>` command,
the "half" neighbor lists may contain **all** pairs of atoms where
atom *j* is a ghost atom (i.e. when the newton pair setting is *off*)
For the newton pair *on* setting the atom *j* is only added to the
list if its *z* coordinate is larger, or if equal the *y* coordinate
is larger, and that is equal, too, the *x* coordinate is larger. For
homogeneously dense systems that will result in picking neighbors from
a same size sector in always the same direction relative to the
"owned" atom and thus it should lead to similar length neighbor lists
and thus reduce the chance of a load imbalance.

114
doc/src/Developer_par_openmp.rst Normal file
View File

@ -0,0 +1,114 @@
OpenMP Parallelism
^^^^^^^^^^^^^^^^^^
The styles in the INTEL, KOKKOS, and OPENMP package offer to use OpenMP
thread parallelism to predominantly distribute loops over local data
and thus follow an orthogonal parallelization strategy to the
decomposition into spatial domains used by the :doc:`MPI partitioning
<Developer_par_part>`. For clarity, this section discusses only the
implementation in the OPENMP package as it is the simplest. The INTEL
and KOKKOS package offer additional options and are more complex since
they support more features and different hardware like co-processors
or GPUs.
One of the key decisions when implementing the OPENMP package was to
keep the changes to the source code small, so that it would be easier to
maintain the code and keep it in sync with the non-threaded standard
implementation. this is achieved by a) making the OPENMP version a
derived class from the regular version (e.g. ``PairLJCutOMP`` from
``PairLJCut``) and overriding only methods that are multi-threaded or
need to be modified to support multi-threading (similar to what was done
in the OPT package), b) keeping the structure in the modified code very
similar so that side-by-side comparisons are still useful, and c)
offloading additional functionality and multi-thread support functions
into three separate classes ``ThrOMP``, ``ThrData``, and ``FixOMP``.
``ThrOMP`` provides additional, multi-thread aware functionality not
available in the corresponding base class (e.g. ``Pair`` for
``PairLJCutOMP``) like multi-thread aware variants of the "tally"
functions. Those functions are made available through multiple
inheritance so those new functions have to have unique names to avoid
ambiguities; typically ``_thr`` is appended to the name of the function.
``ThrData`` is a classes that manages per-thread data structures.
It is used instead of extending the corresponding storage to per-thread
arrays to avoid slowdowns due to "false sharing" when multiple threads
update adjacent elements in an array and thus force the CPU cache lines
to be reset and re-fetched. ``FixOMP`` finally manages the "multi-thread
state" like settings and access to per-thread storage, it is activated
by the :doc:`package omp <package>` command.
Avoiding data races
"""""""""""""""""""
A key problem when implementing thread parallelism in an MD code is
to avoid data races when updating accumulated properties like forces,
energies, and stresses. When interactions are computed, they always
involve multiple atoms and thus there are race conditions when multiple
threads want to update per-atom data of the same atoms. Five possible
strategies have been considered to avoid this:
1) restructure the code so that there is no overlapping access possible
when computing in parallel, e.g. by breaking lists into multiple
parts and synchronizing threads in between.
2) have each thread be "responsible" for a specific group of atoms and
compute these interactions multiple times, once on each thread that
is responsible for a given atom and then have each thread only update
the properties of this atom.
3) use mutexes around functions and regions of code where the data race
could happen
4) use atomic operations when updating per-atom properties
5) use replicated per-thread data structures to accumulate data without
conflicts and then use a reduction to combine those results into the
data structures used by the regular style.
Option 5 was chosen for the OPENMP package because it would retain the
performance for the case of 1 thread and the code would be more
maintainable. Option 1 would require extensive code changes,
particularly to the neighbor list code; options 2 would have incurred a
2x or more performance penalty for the serial case; option 3 causes
significant overhead and would enforce serialization of operations in
inner loops and thus defeat the purpose of multi-threading; option 4
slows down the serial case although not quite as bad as option 2. The
downside of option 5 is that the overhead of the reduction operations
grows with the number of threads used, so there would be a crossover
point where options 2 or 4 would result in faster executing. That is
why option 2 for example is used in the GPU package because a GPU is a
processor with a massive number of threads. However, since the MPI
parallelization is generally more effective for typical MD systems, the
expectation is that thread parallelism is only used for a smaller number
of threads (2-8). At the time of its implementation, that number was
equivalent to the number of CPU cores per CPU socket on high-end
supercomputers.
Thus arrays like the force array are dimensioned to the number of atoms
times the number of threads when enabling OpenMP support and inside the
compute functions a pointer to a different chunk is obtained by each thread.
Similarly, accumulators like potential energy or virial are kept in
per-thread instances of the ``ThrData`` class and then only reduced and
stored in their global counterparts at the end of the force computation.
Loop scheduling
"""""""""""""""
Multi-thread parallelization is applied by distributing (outer) loops
statically across threads. Typically this would be the loop over local
atoms *i* when processing *i,j* pairs of atoms from a neighbor list.
The design of the neighbor list code results in atoms having a similar
number of neighbors for homogeneous systems and thus load imbalances
across threads are not common and typically happen for systems where
also the MPI parallelization would be unbalanced, which would typically
have a more pronounced impact on the performance. This same loop
scheduling scheme can also be applied to the reduction operations on
per-atom data to try and reduce the overhead of the reduction operation.
Neighbor list parallelization
"""""""""""""""""""""""""""""
In addition to the parallelization of force computations, also the
generation of the neighbor lists is parallelized. As explained
previously, neighbor lists are built by looping over "owned" atoms and
storing the neighbors in "pages". In the OPENMP variants of the
neighbor list code, each thread operates on a different chunk of "owned"
atoms and allocates and fills its own set of pages with neighbor list
data. This is achieved by each thread keeping its own instance of the
:cpp:class:`MyPage <LAMMPS_NS::MyPage>` page allocator class.

89
doc/src/Developer_par_part.rst Normal file
View File

@ -0,0 +1,89 @@
Partitioning
^^^^^^^^^^^^
The underlying spatial decomposition strategy used by LAMMPS for
distributed-memory parallelism is set with the :doc:`comm_style command
<comm_style>` and can be either "brick" (a regular grid) or "tiled".
.. _domain-decomposition:
.. figure:: img/domain-decomp.png
:align: center
domain decomposition
This figure shows the different kinds of domain decomposition used
for MPI parallelization: "brick" on the left with an orthogonal
(left) and a triclinic (middle) simulation domain, and a "tiled"
decomposition (right). The black lines show the division into
sub-domains and the contained atoms are "owned" by the corresponding
MPI process. The green dashed lines indicate how sub-domains are
extended with "ghost" atoms up to the communication cutoff distance.
The LAMMPS simulation box is a 3d or 2d volume, which can be orthogonal
or triclinic in shape, as illustrated in the :ref:`domain-decomposition`
figure for the 2d case. Orthogonal means the box edges are aligned with
the *x*, *y*, *z* Cartesian axes, and the box faces are thus all
rectangular. Triclinic allows for a more general parallelepiped shape
in which edges are aligned with three arbitrary vectors and the box
faces are parallelograms. In each dimension box faces can be periodic,
or non-periodic with fixed or shrink-wrapped boundaries. In the fixed
case, atoms which move outside the face are deleted; shrink-wrapped
means the position of the box face adjusts continuously to enclose all
the atoms.
For distributed-memory MPI parallelism, the simulation box is spatially
decomposed (partitioned) into non-overlapping sub-domains which fill the
box. The default partitioning, "brick", is most suitable when atom
density is roughly uniform, as shown in the left-side images of the
:ref:`domain-decomposition` figure. The sub-domains comprise a regular
grid and all sub-domains are identical in size and shape. Both the
orthogonal and triclinic boxes can deform continuously during a
simulation, e.g. to compress a solid or shear a liquid, in which case
the processor sub-domains likewise deform.
For models with non-uniform density, the number of particles per
processor can be load-imbalanced with the default partitioning. This
reduces parallel efficiency, as the overall simulation rate is limited
by the slowest processor, i.e. the one with the largest computational
load. For such models, LAMMPS supports multiple strategies to reduce
the load imbalance:
- The processor grid decomposition is by default based on the simulation
cell volume and tries to optimize the volume to surface ratio for the sub-domains.
This can be changed with the :doc:`processors command <processors>`.
- The parallel planes defining the size of the sub-domains can be shifted
with the :doc:`balance command <balance>`. Which can be done in addition
to choosing a more optimal processor grid.
- The recursive bisectioning algorithm in combination with the "tiled"
communication style can produce a partitioning with equal numbers of
particles in each sub-domain.
.. |decomp1| image:: img/decomp-regular.png
:width: 24%
.. |decomp2| image:: img/decomp-processors.png
:width: 24%
.. |decomp3| image:: img/decomp-balance.png
:width: 24%
.. |decomp4| image:: img/decomp-rcb.png
:width: 24%
|decomp1| |decomp2| |decomp3| |decomp4|
The pictures above demonstrate different decompositions for a 2d system
with 12 MPI ranks. The atom colors indicate the load imbalance of each
sub-domain with green being optimal and red the least optimal.
Due to the vacuum in the system, the default decomposition is unbalanced
with several MPI ranks without atoms (left). By forcing a 1x12x1
processor grid, every MPI rank does computations now, but number of
atoms per sub-domain is still uneven and the thin slice shape increases
the amount of communication between sub-domains (center left). With a
2x6x1 processor grid and shifting the sub-domain divisions, the load
imbalance is further reduced and the amount of communication required
between sub-domains is less (center right). And using the recursive
bisectioning leads to further improved decomposition (right).

28
doc/src/Developer_parallel.rst Normal file
View File

@ -0,0 +1,28 @@
Parallel algorithms
-------------------
LAMMPS is designed to enable running simulations in parallel using the
MPI parallel communication standard with distributed data via domain
decomposition. The parallelization aims to be efficient result in good
strong scaling (= good speedup for the same system) and good weak
scaling (= the computational cost of enlarging the system is
proportional to the system size). Additional parallelization using GPUs
or OpenMP can also be applied within the sub-domain assigned to an MPI
process. For clarity, most of the following illustrations show the 2d
simulation case. The underlying algorithms in those cases, however,
apply to both 2d and 3d cases equally well.
.. note::
The text and most of the figures in this chapter were adapted
for the manual from the section on parallel algorithms in the
:ref:`new LAMMPS paper <lammps_paper>`.
.. toctree::
:maxdepth: 1
Developer_par_part
Developer_par_comm
Developer_par_neigh
Developer_par_long
Developer_par_openmp

6
doc/src/Developer_utils.rst
View File

@ -60,6 +60,9 @@ silently returning the result of a partial conversion or zero in cases
where the string is not a valid number. This behavior allows to more
easily detect typos or issues when processing input files.
Similarly the :cpp:func:`logical() <LAMMPS_NS::utils::logical>` function
will convert a string into a boolean and will only accept certain words.
The *do_abort* flag should be set to ``true`` in case this function
is called only on a single MPI rank, as that will then trigger the
a call to ``Error::one()`` for errors instead of ``Error::all()``
@ -83,6 +86,9 @@ strings for compliance without conversion.
.. doxygenfunction:: tnumeric
:project: progguide
.. doxygenfunction:: logical
:project: progguide
String processing
^^^^^^^^^^^^^^^^^

80
doc/src/Howto_thermostat.rst
View File

@ -2,8 +2,8 @@ Thermostats
===========
Thermostatting means controlling the temperature of particles in an MD
simulation. :doc:`Barostatting <Howto_barostat>` means controlling the
pressure. Since the pressure includes a kinetic component due to
simulation. :doc:`Barostatting <Howto_barostat>` means controlling
the pressure. Since the pressure includes a kinetic component due to
particle velocities, both these operations require calculation of the
temperature. Typically a target temperature (T) and/or pressure (P)
is specified by the user, and the thermostat or barostat attempts to
@ -26,11 +26,13 @@ can be invoked via the *dpd/tstat* pair style:
* :doc:`pair_style dpd/tstat <pair_dpd>`
:doc:`Fix nvt <fix_nh>` only thermostats the translational velocity of
particles. :doc:`Fix nvt/sllod <fix_nvt_sllod>` also does this, except
that it subtracts out a velocity bias due to a deforming box and
integrates the SLLOD equations of motion. See the :doc:`Howto nemd <Howto_nemd>` page for further details. :doc:`Fix nvt/sphere <fix_nvt_sphere>` and :doc:`fix nvt/asphere <fix_nvt_asphere>` thermostat not only translation
velocities but also rotational velocities for spherical and aspherical
particles.
particles. :doc:`Fix nvt/sllod <fix_nvt_sllod>` also does this,
except that it subtracts out a velocity bias due to a deforming box
and integrates the SLLOD equations of motion. See the :doc:`Howto
nemd <Howto_nemd>` page for further details. :doc:`Fix nvt/sphere
<fix_nvt_sphere>` and :doc:`fix nvt/asphere <fix_nvt_asphere>`
thermostat not only translation velocities but also rotational
velocities for spherical and aspherical particles.
.. note::
@ -40,25 +42,31 @@ particles.
e.g. molecular systems. The latter can be tricky to do correctly.
DPD thermostatting alters pairwise interactions in a manner analogous
to the per-particle thermostatting of :doc:`fix langevin <fix_langevin>`.
to the per-particle thermostatting of :doc:`fix langevin
<fix_langevin>`.
Any of the thermostatting fixes can be instructed to use custom temperature
computes that remove bias which has two effects: first, the current
calculated temperature, which is compared to the requested target temperature,
is calculated with the velocity bias removed; second, the thermostat adjusts
only the thermal temperature component of the particle's velocities, which are
the velocities with the bias removed. The removed bias is then added back
to the adjusted velocities. See the doc pages for the individual
fixes and for the :doc:`fix_modify <fix_modify>` command for
instructions on how to assign a temperature compute to a
thermostatting fix. For example, you can apply a thermostat to only
the x and z components of velocity by using it in conjunction with
:doc:`compute temp/partial <compute_temp_partial>`. Of you could
thermostat only the thermal temperature of a streaming flow of
particles without affecting the streaming velocity, by using
:doc:`compute temp/profile <compute_temp_profile>`.
Any of the thermostatting fixes can be instructed to use custom
temperature computes that remove bias which has two effects: first,
the current calculated temperature, which is compared to the requested
target temperature, is calculated with the velocity bias removed;
second, the thermostat adjusts only the thermal temperature component
of the particle's velocities, which are the velocities with the bias
removed. The removed bias is then added back to the adjusted
velocities. See the doc pages for the individual fixes and for the
:doc:`fix_modify <fix_modify>` command for instructions on how to
assign a temperature compute to a thermostatting fix.
Below is a list of some custom temperature computes that can be used like that:
For example, you can apply a thermostat only to atoms in a spatial
region by using it in conjunction with :doc:`compute temp/region
<compute_temp_region>`. Or you can apply a thermostat to only the x
and z components of velocity by using it with :doc:`compute
temp/partial <compute_temp_partial>`. Of you could thermostat only
the thermal temperature of a streaming flow of particles without
affecting the streaming velocity, by using :doc:`compute temp/profile
<compute_temp_profile>`.
Below is a list of custom temperature computes that can be used like
that:
* :doc:`compute_temp_asphere`
* :doc:`compute_temp_body`
@ -72,8 +80,6 @@ Below is a list of some custom temperature computes that can be used like that:
* :doc:`compute_temp_rotate`
* :doc:`compute_temp_sphere`
.. note::
Only the nvt fixes perform time integration, meaning they update
@ -86,17 +92,17 @@ Below is a list of some custom temperature computes that can be used like that:
* :doc:`fix nve/sphere <fix_nve_sphere>`
* :doc:`fix nve/asphere <fix_nve_asphere>`
Thermodynamic output, which can be setup via the
:doc:`thermo_style <thermo_style>` command, often includes temperature
values. As explained on the page for the
:doc:`thermo_style <thermo_style>` command, the default temperature is
setup by the thermo command itself. It is NOT the temperature
associated with any thermostatting fix you have defined or with any
compute you have defined that calculates a temperature. The doc pages
for the thermostatting fixes explain the ID of the temperature compute
they create. Thus if you want to view these temperatures, you need to
specify them explicitly via the :doc:`thermo_style custom <thermo_style>` command. Or you can use the
:doc:`thermo_modify <thermo_modify>` command to re-define what
Thermodynamic output, which can be setup via the :doc:`thermo_style
<thermo_style>` command, often includes temperature values. As
explained on the page for the :doc:`thermo_style <thermo_style>`
command, the default temperature is setup by the thermo command
itself. It is NOT the temperature associated with any thermostatting
fix you have defined or with any compute you have defined that
calculates a temperature. The doc pages for the thermostatting fixes
explain the ID of the temperature compute they create. Thus if you
want to view these temperatures, you need to specify them explicitly
via the :doc:`thermo_style custom <thermo_style>` command. Or you can
use the :doc:`thermo_modify <thermo_modify>` command to re-define what
temperature compute is used for default thermodynamic output.
----------

59
doc/src/Intro_citing.rst
View File

@ -4,28 +4,41 @@ Citing LAMMPS
Core Algorithms
^^^^^^^^^^^^^^^
Since LAMMPS is a community project, there is not a single one
publication or reference that describes **all** of LAMMPS.
The canonical publication that describes the foundation, that is
the basic spatial decomposition approach, the neighbor finding,
and basic communications algorithms used in LAMMPS is:
The paper mentioned below is the best overview of LAMMPS, but there are
also publications describing particular models or algorithms implemented
in LAMMPS or complementary software that is has interfaces to. Please
see below for how to cite contributions to LAMMPS.
`S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, J Comp Phys, 117, 1-19 (1995). <http://www.sandia.gov/~sjplimp/papers/jcompphys95.pdf>`_
.. _lammps_paper:
So any project using LAMMPS (or a derivative application using LAMMPS as
a simulation engine) should cite this paper. A new publication
describing the developments and improvements of LAMMPS in the 25 years
since then is currently in preparation.
The latest canonical publication that describes the basic features, the
source code design, the program structure, the spatial decomposition
approach, the neighbor finding, basic communications algorithms, and how
users and developers have contributed to LAMMPS is:
`LAMMPS - A flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comp. Phys. Comm. (accepted 09/2021), DOI:10.1016/j.cpc.2021.108171 <https://doi.org/10.1016/j.cpc.2021.108171>`_
So a project using LAMMPS or a derivative application that uses LAMMPS
as a simulation engine should cite this paper. The paper is expected to
be published in its final form under the same DOI in the first half
of 2022. Please also give the URL of the LAMMPS website in your paper,
namely https://www.lammps.org.
The original publication describing the parallel algorithms used in the
initial versions of LAMMPS is:
`S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, J Comp Phys, 117, 1-19 (1995). <http://www.sandia.gov/~sjplimp/papers/jcompphys95.pdf>`_
DOI for the LAMMPS code
^^^^^^^^^^^^^^^^^^^^^^^
LAMMPS developers use the `Zenodo service at CERN
<https://zenodo.org/>`_ to create digital object identifies (DOI) for
stable releases of the LAMMPS code. There are two types of DOIs for the
LAMMPS source code: the canonical DOI for **all** versions of LAMMPS,
which will always point to the **latest** stable release version is:
LAMMPS developers use the `Zenodo service at CERN <https://zenodo.org/>`_
to create digital object identifies (DOI) for stable releases of the
LAMMPS source code. There are two types of DOIs for the LAMMPS source code.
The canonical DOI for **all** versions of LAMMPS, which will always
point to the **latest** stable release version is:
- DOI: `10.5281/zenodo.3726416 <https://dx.doi.org/10.5281/zenodo.3726416>`_
@ -45,11 +58,13 @@ about LAMMPS and its features.
Citing contributions
^^^^^^^^^^^^^^^^^^^^
LAMMPS has many features and that use either previously published
methods and algorithms or novel features. It also includes potential
parameter filed for specific models. Where available, a reminder about
references for optional features used in a specific run is printed to
the screen and log file. Style and output location can be selected with
the :ref:`-cite command-line switch <cite>`. Additional references are
LAMMPS has many features that use either previously published methods
and algorithms or novel features. It also includes potential parameter
files for specific models. Where available, a reminder about references
for optional features used in a specific run is printed to the screen
and log file. Style and output location can be selected with the
:ref:`-cite command-line switch <cite>`. Additional references are
given in the documentation of the :doc:`corresponding commands
<Commands_all>` or in the :doc:`Howto tutorials <Howto>`.
<Commands_all>` or in the :doc:`Howto tutorials <Howto>`. So please
make certain, that you provide the proper acknowledgments and citations
in any published works using LAMMPS.

2
doc/src/Library_create.rst
View File

@ -34,7 +34,7 @@ simple example demonstrating its use:
int lmpargc = sizeof(lmpargv)/sizeof(const char *);
/* create LAMMPS instance */
handle = lammps_open_no_mpi(lmpargc, lmpargv, NULL);
handle = lammps_open_no_mpi(lmpargc, (char **)lmpargv, NULL);
if (handle == NULL) {
printf("LAMMPS initialization failed");
lammps_mpi_finalize();

BIN
doc/src/PDF/colvars-refman-lammps.pdf
View File

Binary file not shown.

16
doc/src/Run_basics.rst
View File

@ -2,17 +2,25 @@ Basics of running LAMMPS
========================
LAMMPS is run from the command line, reading commands from a file via
the -in command line flag, or from standard input.
Using the "-in in.file" variant is recommended:
the -in command line flag, or from standard input. Using the "-in
in.file" variant is recommended (see note below). The name of the
LAMMPS executable is either ``lmp`` or ``lmp_<machine>`` with
`<machine>` being the machine string used when compiling LAMMPS. This
is required when compiling LAMMPS with the traditional build system
(e.g. with ``make mpi``), but optional when using CMake to configure and
build LAMMPS:
.. code-block:: bash
$ lmp_serial -in in.file
$ lmp_serial < in.file
$ lmp -in in.file
$ lmp < in.file
$ /path/to/lammps/src/lmp_serial -i in.file
$ mpirun -np 4 lmp_mpi -in in.file
$ mpiexec -np 4 lmp -in in.file
$ mpirun -np 8 /path/to/lammps/src/lmp_mpi -in in.file
$ mpirun -np 6 /usr/local/bin/lmp -in in.file
$ mpiexec -n 6 /usr/local/bin/lmp -in in.file
You normally run the LAMMPS command in the directory where your input
script is located. That is also where output files are produced by
@ -23,7 +31,7 @@ executable itself can be placed elsewhere.
.. note::
The redirection operator "<" will not always work when running
in parallel with mpirun; for those systems the -in form is required.
in parallel with mpirun or mpiexec; for those systems the -in form is required.
As LAMMPS runs it prints info to the screen and a logfile named
*log.lammps*\ . More info about output is given on the

25
doc/src/fix_langevin.rst
View File

@ -138,16 +138,18 @@ temperature with optional time-dependence as well.
Like other fixes that perform thermostatting, this fix can be used
with :doc:`compute commands <compute>` that remove a "bias" from the
atom velocities. E.g. removing the center-of-mass velocity from a
group of atoms or removing the x-component of velocity from the
calculation. This is not done by default, but only if the
:doc:`fix_modify <fix_modify>` command is used to assign a temperature
compute to this fix that includes such a bias term. See the doc pages
for individual :doc:`compute commands <compute>` to determine which ones
include a bias. In this case, the thermostat works in the following
manner: bias is removed from each atom, thermostatting is performed on
the remaining thermal degrees of freedom, and the bias is added back
in.
atom velocities. E.g. to apply the thermostat only to atoms within a
spatial :doc:`region <region>`, or to remove the center-of-mass
velocity from a group of atoms, or to remove the x-component of
velocity from the calculation.
This is not done by default, but only if the :doc:`fix_modify
<fix_modify>` command is used to assign a temperature compute to this
fix that includes such a bias term. See the doc pages for individual
:doc:`compute temp commands <compute>` to determine which ones include
a bias. In this case, the thermostat works in the following manner:
bias is removed from each atom, thermostatting is performed on the
remaining thermal degrees of freedom, and the bias is added back in.
The *damp* parameter is specified in time units and determines how
rapidly the temperature is relaxed. For example, a value of 100.0 means
@ -183,7 +185,8 @@ omega (which is derived from the angular momentum in the case of
aspherical particles).
The rotational temperature of the particles can be monitored by the
:doc:`compute temp/sphere <compute_temp_sphere>` and :doc:`compute temp/asphere <compute_temp_asphere>` commands with their rotate
:doc:`compute temp/sphere <compute_temp_sphere>` and :doc:`compute
temp/asphere <compute_temp_asphere>` commands with their rotate
options.
For the *omega* keyword there is also a scale factor of

25
doc/src/fix_langevin_drude.rst
View File

@ -167,17 +167,20 @@ functions, and include :doc:`thermo_style <thermo_style>` command
keywords for the simulation box parameters and timestep and elapsed
time. Thus it is easy to specify a time-dependent temperature.
Like other fixes that perform thermostatting, this fix can be used with
:doc:`compute commands <compute>` that remove a "bias" from the atom
velocities. E.g. removing the center-of-mass velocity from a group of
atoms. This is not done by default, but only if the
:doc:`fix_modify <fix_modify>` command is used to assign a temperature
compute to this fix that includes such a bias term. See the doc pages
for individual :doc:`compute commands <compute>` to determine which ones
include a bias. In this case, the thermostat works in the following
manner: bias is removed from each atom, thermostatting is performed on
the remaining thermal degrees of freedom, and the bias is added back
in. NOTE: this feature has not been tested.
Like other fixes that perform thermostatting, this fix can be used
with :doc:`compute commands <compute>` that remove a "bias" from the
atom velocities. E.g. to apply the thermostat only to atoms within a
spatial :doc:`region <region>`, or to remove the center-of-mass
velocity from a group of atoms, or to remove the x-component of
velocity from the calculation.
This is not done by default, but only if the :doc:`fix_modify
<fix_modify>` command is used to assign a temperature compute to this
fix that includes such a bias term. See the doc pages for individual
:doc:`compute temp commands <compute>` to determine which ones include
a bias. In this case, the thermostat works in the following manner:
bias is removed from each atom, thermostatting is performed on the
remaining thermal degrees of freedom, and the bias is added back in.
Note: The temperature thermostatting the core-Drude particle pairs
should be chosen low enough, so as to mimic as closely as possible the

27
doc/src/fix_nh.rst
View File

@ -486,19 +486,20 @@ temperature or pressure during thermodynamic output via the
compute-ID. It also means that changing attributes of *thermo_temp*
or *thermo_press* will have no effect on this fix.
Like other fixes that perform thermostatting, fix nvt and fix npt can
be used with :doc:`compute commands <compute>` that calculate a
temperature after removing a "bias" from the atom velocities.
E.g. removing the center-of-mass velocity from a group of atoms or
only calculating temperature on the x-component of velocity or only
calculating temperature for atoms in a geometric region. This is not
done by default, but only if the :doc:`fix_modify <fix_modify>` command
is used to assign a temperature compute to this fix that includes such
a bias term. See the doc pages for individual :doc:`compute commands <compute>` to determine which ones include a bias. In
this case, the thermostat works in the following manner: the current
temperature is calculated taking the bias into account, bias is
removed from each atom, thermostatting is performed on the remaining
thermal degrees of freedom, and the bias is added back in.
Like other fixes that perform thermostatting, this fix can be used
with :doc:`compute commands <compute>` that remove a "bias" from the
atom velocities. E.g. to apply the thermostat only to atoms within a
spatial :doc:`region <region>`, or to remove the center-of-mass
velocity from a group of atoms, or to remove the x-component of
velocity from the calculation.
This is not done by default, but only if the :doc:`fix_modify
<fix_modify>` command is used to assign a temperature compute to this
fix that includes such a bias term. See the doc pages for individual
:doc:`compute temp commands <compute>` to determine which ones include
a bias. In this case, the thermostat works in the following manner:
bias is removed from each atom, thermostatting is performed on the
remaining thermal degrees of freedom, and the bias is added back in.
----------

30
doc/src/fix_npt_asphere.rst
View File

@ -48,8 +48,9 @@ can also have a bias velocity removed from them before thermostatting
takes place; see the description below.
Additional parameters affecting the thermostat and barostat are
specified by keywords and values documented with the :doc:`fix npt <fix_nh>` command. See, for example, discussion of the *temp*,
*iso*, *aniso*, and *dilate* keywords.
specified by keywords and values documented with the :doc:`fix npt
<fix_nh>` command. See, for example, discussion of the *temp*, *iso*,
*aniso*, and *dilate* keywords.
The particles in the fix group are the only ones whose velocities and
positions are updated by the velocity/position update portion of the
@ -89,18 +90,19 @@ It also means that changing attributes of *thermo_temp* or
*thermo_press* will have no effect on this fix.
Like other fixes that perform thermostatting, this fix can be used
with :doc:`compute commands <compute>` that calculate a temperature
after removing a "bias" from the atom velocities. E.g. removing the
center-of-mass velocity from a group of atoms or only calculating
temperature on the x-component of velocity or only calculating
temperature for atoms in a geometric region. This is not done by
default, but only if the :doc:`fix_modify <fix_modify>` command is used
to assign a temperature compute to this fix that includes such a bias
term. See the doc pages for individual :doc:`compute commands <compute>` to determine which ones include a bias. In
this case, the thermostat works in the following manner: the current
temperature is calculated taking the bias into account, bias is
removed from each atom, thermostatting is performed on the remaining
thermal degrees of freedom, and the bias is added back in.
with :doc:`compute commands <compute>` that remove a "bias" from the
atom velocities. E.g. to apply the thermostat only to atoms within a
spatial :doc:`region <region>`, or to remove the center-of-mass
velocity from a group of atoms, or to remove the x-component of
velocity from the calculation.
This is not done by default, but only if the :doc:`fix_modify
<fix_modify>` command is used to assign a temperature compute to this
fix that includes such a bias term. See the doc pages for individual
:doc:`compute temp commands <compute>` to determine which ones include
a bias. In this case, the thermostat works in the following manner:
bias is removed from each atom, thermostatting is performed on the
remaining thermal degrees of freedom, and the bias is added back in.
----------

25
doc/src/fix_npt_body.rst
View File

@ -87,18 +87,19 @@ It also means that changing attributes of *thermo_temp* or
*thermo_press* will have no effect on this fix.
Like other fixes that perform thermostatting, this fix can be used
with :doc:`compute commands <compute>` that calculate a temperature
after removing a "bias" from the atom velocities. E.g. removing the
center-of-mass velocity from a group of atoms or only calculating
temperature on the x-component of velocity or only calculating
temperature for atoms in a geometric region. This is not done by
default, but only if the :doc:`fix_modify <fix_modify>` command is used
to assign a temperature compute to this fix that includes such a bias
term. See the doc pages for individual :doc:`compute commands <compute>` to determine which ones include a bias. In
this case, the thermostat works in the following manner: the current
temperature is calculated taking the bias into account, bias is
removed from each atom, thermostatting is performed on the remaining
thermal degrees of freedom, and the bias is added back in.
with :doc:`compute commands <compute>` that remove a "bias" from the
atom velocities. E.g. to apply the thermostat only to atoms within a
spatial :doc:`region <region>`, or to remove the center-of-mass
velocity from a group of atoms, or to remove the x-component of
velocity from the calculation.
This is not done by default, but only if the :doc:`fix_modify
<fix_modify>` command is used to assign a temperature compute to this
fix that includes such a bias term. See the doc pages for individual
:doc:`compute temp commands <compute>` to determine which ones include
a bias. In this case, the thermostat works in the following manner:
bias is removed from each atom, thermostatting is performed on the
remaining thermal degrees of freedom, and the bias is added back in.
----------

27
doc/src/fix_npt_cauchy.rst
View File

@ -400,19 +400,20 @@ temperature or pressure during thermodynamic output via the
compute-ID. It also means that changing attributes of *thermo_temp*
or *thermo_press* will have no effect on this fix.
Like other fixes that perform thermostatting, fix npt/cauchy can
be used with :doc:`compute commands <compute>` that calculate a
temperature after removing a "bias" from the atom velocities.
E.g. removing the center-of-mass velocity from a group of atoms or
only calculating temperature on the x-component of velocity or only
calculating temperature for atoms in a geometric region. This is not
done by default, but only if the :doc:`fix_modify <fix_modify>` command
is used to assign a temperature compute to this fix that includes such
a bias term. See the doc pages for individual :doc:`compute commands <compute>` to determine which ones include a bias. In
this case, the thermostat works in the following manner: the current
temperature is calculated taking the bias into account, bias is
removed from each atom, thermostatting is performed on the remaining
thermal degrees of freedom, and the bias is added back in.
Like other fixes that perform thermostatting, this fix can be used
with :doc:`compute commands <compute>` that remove a "bias" from the
atom velocities. E.g. to apply the thermostat only to atoms within a
spatial :doc:`region <region>`, or to remove the center-of-mass
velocity from a group of atoms, or to remove the x-component of
velocity from the calculation.
This is not done by default, but only if the :doc:`fix_modify
<fix_modify>` command is used to assign a temperature compute to this
fix that includes such a bias term. See the doc pages for individual
:doc:`compute temp commands <compute>` to determine which ones include
a bias. In this case, the thermostat works in the following manner:
bias is removed from each atom, thermostatting is performed on the
remaining thermal degrees of freedom, and the bias is added back in.
----------

25
doc/src/fix_npt_sphere.rst
View File

@ -103,18 +103,19 @@ appropriate compute-ID. It also means that changing attributes of
*thermo_temp* or *thermo_press* will have no effect on this fix.
Like other fixes that perform thermostatting, this fix can be used
with :doc:`compute commands <compute>` that calculate a temperature
after removing a "bias" from the atom velocities. E.g. removing the
center-of-mass velocity from a group of atoms or only calculating
temperature on the x-component of velocity or only calculating
temperature for atoms in a geometric region. This is not done by
default, but only if the :doc:`fix_modify <fix_modify>` command is used
to assign a temperature compute to this fix that includes such a bias
term. See the doc pages for individual :doc:`compute commands <compute>` to determine which ones include a bias. In
this case, the thermostat works in the following manner: the current
temperature is calculated taking the bias into account, bias is
removed from each atom, thermostatting is performed on the remaining
thermal degrees of freedom, and the bias is added back in.
with :doc:`compute commands <compute>` that remove a "bias" from the
atom velocities. E.g. to apply the thermostat only to atoms within a
spatial :doc:`region <region>`, or to remove the center-of-mass
velocity from a group of atoms, or to remove the x-component of
velocity from the calculation.
This is not done by default, but only if the :doc:`fix_modify
<fix_modify>` command is used to assign a temperature compute to this
fix that includes such a bias term. See the doc pages for individual
:doc:`compute temp commands <compute>` to determine which ones include
a bias. In this case, the thermostat works in the following manner:
bias is removed from each atom, thermostatting is performed on the
remaining thermal degrees of freedom, and the bias is added back in.
----------

25
doc/src/fix_nvt_asphere.rst
View File

@ -72,18 +72,19 @@ It also means that changing attributes of *thermo_temp* will have no
effect on this fix.
Like other fixes that perform thermostatting, this fix can be used
with :doc:`compute commands <compute>` that calculate a temperature
after removing a "bias" from the atom velocities. E.g. removing the
center-of-mass velocity from a group of atoms or only calculating
temperature on the x-component of velocity or only calculating
temperature for atoms in a geometric region. This is not done by
default, but only if the :doc:`fix_modify <fix_modify>` command is used
to assign a temperature compute to this fix that includes such a bias
term. See the doc pages for individual :doc:`compute commands <compute>` to determine which ones include a bias. In
this case, the thermostat works in the following manner: the current
temperature is calculated taking the bias into account, bias is
removed from each atom, thermostatting is performed on the remaining
thermal degrees of freedom, and the bias is added back in.
with :doc:`compute commands <compute>` that remove a "bias" from the
atom velocities. E.g. to apply the thermostat only to atoms within a
spatial :doc:`region <region>`, or to remove the center-of-mass
velocity from a group of atoms, or to remove the x-component of
velocity from the calculation.
This is not done by default, but only if the :doc:`fix_modify
<fix_modify>` command is used to assign a temperature compute to this
fix that includes such a bias term. See the doc pages for individual
:doc:`compute temp commands <compute>` to determine which ones include
a bias. In this case, the thermostat works in the following manner:
bias is removed from each atom, thermostatting is performed on the
remaining thermal degrees of freedom, and the bias is added back in.
----------

25
doc/src/fix_nvt_body.rst
View File

@ -69,18 +69,19 @@ It also means that changing attributes of *thermo_temp* will have no
effect on this fix.
Like other fixes that perform thermostatting, this fix can be used
with :doc:`compute commands <compute>` that calculate a temperature
after removing a "bias" from the atom velocities. E.g. removing the
center-of-mass velocity from a group of atoms or only calculating
temperature on the x-component of velocity or only calculating
temperature for atoms in a geometric region. This is not done by
default, but only if the :doc:`fix_modify <fix_modify>` command is used
to assign a temperature compute to this fix that includes such a bias
term. See the doc pages for individual :doc:`compute commands <compute>` to determine which ones include a bias. In
this case, the thermostat works in the following manner: the current
temperature is calculated taking the bias into account, bias is
removed from each atom, thermostatting is performed on the remaining
thermal degrees of freedom, and the bias is added back in.
with :doc:`compute commands <compute>` that remove a "bias" from the
atom velocities. E.g. to apply the thermostat only to atoms within a
spatial :doc:`region <region>`, or to remove the center-of-mass
velocity from a group of atoms, or to remove the x-component of
velocity from the calculation.
This is not done by default, but only if the :doc:`fix_modify
<fix_modify>` command is used to assign a temperature compute to this
fix that includes such a bias term. See the doc pages for individual
:doc:`compute temp commands <compute>` to determine which ones include
a bias. In this case, the thermostat works in the following manner:
bias is removed from each atom, thermostatting is performed on the
remaining thermal degrees of freedom, and the bias is added back in.
----------

74
doc/src/fix_nvt_sllod.rst
View File

@ -37,15 +37,16 @@ trajectory consistent with the canonical ensemble.
This thermostat is used for a simulation box that is changing size
and/or shape, for example in a non-equilibrium MD (NEMD) simulation.
The size/shape change is induced by use of the :doc:`fix deform <fix_deform>` command, so each point in the simulation box
can be thought of as having a "streaming" velocity. This
position-dependent streaming velocity is subtracted from each atom's
actual velocity to yield a thermal velocity which is used for
temperature computation and thermostatting. For example, if the box
is being sheared in x, relative to y, then points at the bottom of the
box (low y) have a small x velocity, while points at the top of the
box (hi y) have a large x velocity. These velocities do not
contribute to the thermal "temperature" of the atom.
The size/shape change is induced by use of the :doc:`fix deform
<fix_deform>` command, so each point in the simulation box can be
thought of as having a "streaming" velocity. This position-dependent
streaming velocity is subtracted from each atom's actual velocity to
yield a thermal velocity which is used for temperature computation and
thermostatting. For example, if the box is being sheared in x,
relative to y, then points at the bottom of the box (low y) have a
small x velocity, while points at the top of the box (hi y) have a
large x velocity. These velocities do not contribute to the thermal
"temperature" of the atom.
.. note::
@ -60,13 +61,15 @@ contribute to the thermal "temperature" of the atom.
consistent.
The SLLOD equations of motion, originally proposed by Hoover and Ladd
(see :ref:`(Evans and Morriss) <Evans3>`), were proven to be equivalent to
Newton's equations of motion for shear flow by :ref:`(Evans and Morriss) <Evans3>`. They were later shown to generate the desired
velocity gradient and the correct production of work by stresses for
all forms of homogeneous flow by :ref:`(Daivis and Todd) <Daivis>`. As
implemented in LAMMPS, they are coupled to a Nose/Hoover chain
thermostat in a velocity Verlet formulation, closely following the
implementation used for the :doc:`fix nvt <fix_nh>` command.
(see :ref:`(Evans and Morriss) <Evans3>`), were proven to be
equivalent to Newton's equations of motion for shear flow by
:ref:`(Evans and Morriss) <Evans3>`. They were later shown to generate
the desired velocity gradient and the correct production of work by
stresses for all forms of homogeneous flow by :ref:`(Daivis and Todd)
<Daivis>`. As implemented in LAMMPS, they are coupled to a
Nose/Hoover chain thermostat in a velocity Verlet formulation, closely
following the implementation used for the :doc:`fix nvt <fix_nh>`
command.
.. note::
@ -94,27 +97,28 @@ underscore + "temp", and the group for the new compute is the same as
the fix group.
Note that this is NOT the compute used by thermodynamic output (see
the :doc:`thermo_style <thermo_style>` command) with ID = *thermo_temp*.
This means you can change the attributes of this fix's temperature
(e.g. its degrees-of-freedom) via the
:doc:`compute_modify <compute_modify>` command or print this temperature
during thermodynamic output via the :doc:`thermo_style custom <thermo_style>` command using the appropriate compute-ID.
It also means that changing attributes of *thermo_temp* will have no
effect on this fix.
the :doc:`thermo_style <thermo_style>` command) with ID =
*thermo_temp*. This means you can change the attributes of this fix's
temperature (e.g. its degrees-of-freedom) via the :doc:`compute_modify
<compute_modify>` command or print this temperature during
thermodynamic output via the :doc:`thermo_style custom <thermo_style>`
command using the appropriate compute-ID. It also means that changing
attributes of *thermo_temp* will have no effect on this fix.
Like other fixes that perform thermostatting, this fix can be used
with :doc:`compute commands <compute>` that calculate a temperature
after removing a "bias" from the atom velocities. E.g. removing the
center-of-mass velocity from a group of atoms or only calculating
temperature on the x-component of velocity or only calculating
temperature for atoms in a geometric region. This is not done by
default, but only if the :doc:`fix_modify <fix_modify>` command is used
to assign a temperature compute to this fix that includes such a bias
term. See the doc pages for individual :doc:`compute commands <compute>` to determine which ones include a bias. In
this case, the thermostat works in the following manner: the current
temperature is calculated taking the bias into account, bias is
removed from each atom, thermostatting is performed on the remaining
thermal degrees of freedom, and the bias is added back in.
with :doc:`compute commands <compute>` that remove a "bias" from the
atom velocities. E.g. to apply the thermostat only to atoms within a
spatial :doc:`region <region>`, or to remove the center-of-mass
velocity from a group of atoms, or to remove the x-component of
velocity from the calculation.
This is not done by default, but only if the :doc:`fix_modify
<fix_modify>` command is used to assign a temperature compute to this
fix that includes such a bias term. See the doc pages for individual
:doc:`compute temp commands <compute>` to determine which ones include
a bias. In this case, the thermostat works in the following manner:
bias is removed from each atom, thermostatting is performed on the
remaining thermal degrees of freedom, and the bias is added back in.
----------

25
doc/src/fix_nvt_sphere.rst
View File

@ -86,18 +86,19 @@ It also means that changing attributes of *thermo_temp* will have no
effect on this fix.
Like other fixes that perform thermostatting, this fix can be used
with :doc:`compute commands <compute>` that calculate a temperature
after removing a "bias" from the atom velocities. E.g. removing the
center-of-mass velocity from a group of atoms or only calculating
temperature on the x-component of velocity or only calculating
temperature for atoms in a geometric region. This is not done by
default, but only if the :doc:`fix_modify <fix_modify>` command is used
to assign a temperature compute to this fix that includes such a bias
term. See the doc pages for individual :doc:`compute commands <compute>` to determine which ones include a bias. In
this case, the thermostat works in the following manner: the current
temperature is calculated taking the bias into account, bias is
removed from each atom, thermostatting is performed on the remaining
thermal degrees of freedom, and the bias is added back in.
with :doc:`compute commands <compute>` that remove a "bias" from the
atom velocities. E.g. to apply the thermostat only to atoms within a
spatial :doc:`region <region>`, or to remove the center-of-mass
velocity from a group of atoms, or to remove the x-component of
velocity from the calculation.
This is not done by default, but only if the :doc:`fix_modify
<fix_modify>` command is used to assign a temperature compute to this
fix that includes such a bias term. See the doc pages for individual
:doc:`compute temp commands <compute>` to determine which ones include
a bias. In this case, the thermostat works in the following manner:
bias is removed from each atom, thermostatting is performed on the
remaining thermal degrees of freedom, and the bias is added back in.
----------

25
doc/src/fix_temp_berendsen.rst
View File

@ -102,18 +102,19 @@ It also means that changing attributes of *thermo_temp* will have no
effect on this fix.
Like other fixes that perform thermostatting, this fix can be used
with :doc:`compute commands <compute>` that calculate a temperature
after removing a "bias" from the atom velocities. E.g. removing the
center-of-mass velocity from a group of atoms or only calculating
temperature on the x-component of velocity or only calculating
temperature for atoms in a geometric region. This is not done by
default, but only if the :doc:`fix_modify <fix_modify>` command is used
to assign a temperature compute to this fix that includes such a bias
term. See the doc pages for individual :doc:`compute commands <compute>` to determine which ones include a bias. In
this case, the thermostat works in the following manner: the current
temperature is calculated taking the bias into account, bias is
removed from each atom, thermostatting is performed on the remaining
thermal degrees of freedom, and the bias is added back in.
with :doc:`compute commands <compute>` that remove a "bias" from the
atom velocities. E.g. to apply the thermostat only to atoms within a
spatial :doc:`region <region>`, or to remove the center-of-mass
velocity from a group of atoms, or to remove the x-component of
velocity from the calculation.
This is not done by default, but only if the :doc:`fix_modify
<fix_modify>` command is used to assign a temperature compute to this
fix that includes such a bias term. See the doc pages for individual
:doc:`compute temp commands <compute>` to determine which ones include
a bias. In this case, the thermostat works in the following manner:
bias is removed from each atom, thermostatting is performed on the
remaining thermal degrees of freedom, and the bias is added back in.
----------

41
doc/src/fix_temp_csvr.rst
View File

@ -110,28 +110,29 @@ during thermodynamic output via the :doc:`thermo_style custom <thermo_style>` co
It also means that changing attributes of *thermo_temp* will have no
effect on this fix.
Like other fixes that perform thermostatting, these fixes can be used
with :doc:`compute commands <compute>` that calculate a temperature
after removing a "bias" from the atom velocities. E.g. removing the
center-of-mass velocity from a group of atoms or only calculating
temperature on the x-component of velocity or only calculating
temperature for atoms in a geometric region. This is not done by
default, but only if the :doc:`fix_modify <fix_modify>` command is used
to assign a temperature compute to this fix that includes such a bias
term. See the doc pages for individual :doc:`compute commands <compute>` to determine which ones include a bias. In
this case, the thermostat works in the following manner: the current
temperature is calculated taking the bias into account, bias is
removed from each atom, thermostatting is performed on the remaining
thermal degrees of freedom, and the bias is added back in.
Like other fixes that perform thermostatting, this fix can be used
with :doc:`compute commands <compute>` that remove a "bias" from the
atom velocities. E.g. to apply the thermostat only to atoms within a
spatial :doc:`region <region>`, or to remove the center-of-mass
velocity from a group of atoms, or to remove the x-component of
velocity from the calculation.
This is not done by default, but only if the :doc:`fix_modify
<fix_modify>` command is used to assign a temperature compute to this
fix that includes such a bias term. See the doc pages for individual
:doc:`compute temp commands <compute>` to determine which ones include
a bias. In this case, the thermostat works in the following manner:
bias is removed from each atom, thermostatting is performed on the
remaining thermal degrees of freedom, and the bias is added back in.
An important feature of these thermostats is that they have an
associated effective energy that is a constant of motion.
The effective energy is the total energy (kinetic + potential) plus
the accumulated kinetic energy changes due to the thermostat. The
latter quantity is the global scalar computed by these fixes. This
feature is useful to check the integration of the equations of motion
against discretization errors. In other words, the conservation of
the effective energy can be used to choose an appropriate integration
associated effective energy that is a constant of motion. The
effective energy is the total energy (kinetic + potential) plus the
accumulated kinetic energy changes due to the thermostat. The latter
quantity is the global scalar computed by these fixes. This feature is
useful to check the integration of the equations of motion against
discretization errors. In other words, the conservation of the
effective energy can be used to choose an appropriate integration
:doc:`timestep <timestep>`. This is similar to the usual paradigm of
checking the conservation of the total energy in the microcanonical
ensemble.

26
doc/src/fix_temp_rescale.rst
View File

@ -109,19 +109,19 @@ command using the appropriate compute-ID. It also means that changing
attributes of *thermo_temp* will have no effect on this fix.
Like other fixes that perform thermostatting, this fix can be used
with :doc:`compute commands <compute>` that calculate a temperature
after removing a "bias" from the atom velocities. E.g. removing the
center-of-mass velocity from a group of atoms or only calculating
temperature on the x-component of velocity or only calculating
temperature for atoms in a geometric region. This is not done by
default, but only if the :doc:`fix_modify <fix_modify>` command is
used to assign a temperature compute to this fix that includes such a
bias term. See the doc pages for individual :doc:`compute commands
<compute>` to determine which ones include a bias. In this case, the
thermostat works in the following manner: the current temperature is
calculated taking the bias into account, bias is removed from each
atom, thermostatting is performed on the remaining thermal degrees of
freedom, and the bias is added back in.
with :doc:`compute commands <compute>` that remove a "bias" from the
atom velocities. E.g. to apply the thermostat only to atoms within a
spatial :doc:`region <region>`, or to remove the center-of-mass
velocity from a group of atoms, or to remove the x-component of
velocity from the calculation.
This is not done by default, but only if the :doc:`fix_modify
<fix_modify>` command is used to assign a temperature compute to this
fix that includes such a bias term. See the doc pages for individual
:doc:`compute temp commands <compute>` to determine which ones include
a bias. In this case, the thermostat works in the following manner:
bias is removed from each atom, thermostatting is performed on the
remaining thermal degrees of freedom, and the bias is added back in.
----------

40
doc/src/fix_tgnh_drude.rst
View File

@ -187,26 +187,32 @@ barostatting.
----------
Like other fixes that perform thermostatting, these fixes can
be used with :doc:`compute commands <compute>` that calculate a
temperature after removing a "bias" from the atom velocities.
This is not done by default, but only if the :doc:`fix_modify <fix_modify>` command
is used to assign a temperature compute to this fix that includes such
a bias term. See the doc pages for individual :doc:`compute commands <compute>` to determine which ones include a bias. In
this case, the thermostat works in the following manner: the current
temperature is calculated taking the bias into account, bias is
removed from each atom, thermostatting is performed on the remaining
thermal DOF, and the bias is added back in.
Like other fixes that perform thermostatting, this fix can be used
with :doc:`compute commands <compute>` that remove a "bias" from the
atom velocities. E.g. to apply the thermostat only to atoms within a
spatial :doc:`region <region>`, or to remove the center-of-mass
velocity from a group of atoms, or to remove the x-component of
velocity from the calculation.
This is not done by default, but only if the :doc:`fix_modify
<fix_modify>` command is used to assign a temperature compute to this
fix that includes such a bias term. See the doc pages for individual
:doc:`compute temp commands <compute>` to determine which ones include
a bias. In this case, the thermostat works in the following manner:
bias is removed from each atom, thermostatting is performed on the
remaining thermal degrees of freedom, and the bias is added back in.
.. note::
However, not all temperature compute commands are valid to be used with these fixes.
Precisely, only temperature compute that does not modify the DOF of the group can be used.
E.g. :doc:`compute temp/ramp <compute_temp_ramp>` and :doc:`compute viscosity/cos <compute_viscosity_cos>`
compute the kinetic energy after remove a velocity gradient without affecting the DOF of the group,
then they can be invoked in this way.
In contrast, :doc:`compute temp/partial <compute_temp_partial>` may remove the DOF at one or more dimensions,
therefore it cannot be used with these fixes.
However, not all temperature compute commands are valid to be used
with these fixes. Precisely, only temperature compute that does
not modify the DOF of the group can be used. E.g. :doc:`compute
temp/ramp <compute_temp_ramp>` and :doc:`compute viscosity/cos
<compute_viscosity_cos>` compute the kinetic energy after remove a
velocity gradient without affecting the DOF of the group, then they
can be invoked in this way. In contrast, :doc:`compute
temp/partial <compute_temp_partial>` may remove the DOF at one or
more dimensions, therefore it cannot be used with these fixes.
----------

2
doc/src/group.rst
View File

@ -38,7 +38,7 @@ Syntax
*intersect* args = two or more group IDs
*dynamic* args = parent-ID keyword value ...
one or more keyword/value pairs may be appended
keyword = *region* or *var* or *every*
keyword = *region* or *var* or *property* or *every*
*region* value = region-ID
*var* value = name of variable
*property* value = name of custom integer or floating point vector

BIN
doc/src/img/decomp-balance.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 129 KiB

BIN
doc/src/img/decomp-processors.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 121 KiB

BIN
doc/src/img/decomp-rcb.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 121 KiB

BIN
doc/src/img/decomp-regular.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 121 KiB

BIN
doc/src/img/domain-decomp.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 547 KiB

BIN
doc/src/img/fft-decomp-parallel.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 35 KiB

BIN
doc/src/img/ghost-comm.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 33 KiB

BIN
doc/src/img/neigh-stencil.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 53 KiB

2
doc/utils/requirements.txt
View File

@ -1,4 +1,4 @@
Sphinx==4.0.3
Sphinx
sphinxcontrib-spelling
git+git://github.com/akohlmey/sphinx-fortran@parallel-read
sphinx_tabs

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

@ -2262,6 +2262,7 @@ Nmols
nn
nnodes
Nocedal
nO
nocite
nocoeff
nodeless
@ -3657,6 +3658,7 @@ Yc
ycm
Yeh
yellowgreen
yEs
Yethiraj
yflag
yhi

2
examples/PACKAGES/charge_regulation/in.chreg-polymer
View File

@ -8,7 +8,7 @@ bond_style harmonic
bond_coeff 1 100 1.122462 # K R0
velocity all create 1.0 8008 loop geom
pair_style lj/cut/coul/long 1.122462 20
pair_style lj/cut/coul/long/soft 2 0.5 10.0 1.122462 20
pair_coeff * * 1.0 1.0 1.122462 # charges
kspace_style pppm 1.0e-3
pair_modify shift yes

4
lib/colvars/colvarmodule.cpp
View File

@ -1476,7 +1476,9 @@ int colvarmodule::write_output_files()
bi != biases.end();
bi++) {
// Only write output files if they have not already been written this time step
if ((*bi)->output_freq == 0 || (cvm::step_absolute() % (*bi)->output_freq) != 0) {
if ((*bi)->output_freq == 0 ||
cvm::step_relative() == 0 ||
(cvm::step_absolute() % (*bi)->output_freq) != 0) {
error_code |= (*bi)->write_output_files();
}
error_code |= (*bi)->write_state_to_replicas();

2
lib/colvars/colvars_version.h
View File

@ -1,3 +1,3 @@
#ifndef COLVARS_VERSION
#define COLVARS_VERSION "2021-08-06"
#define COLVARS_VERSION "2021-09-21"
#endif

14
lib/gpu/geryon/ocl_device.h
View File

@ -556,16 +556,22 @@ void UCL_Device::add_properties(cl_device_id device_list) {
sizeof(float_width),&float_width,nullptr));
op.preferred_vector_width32=float_width;
// Determine if double precision is supported
cl_uint double_width;
CL_SAFE_CALL(clGetDeviceInfo(device_list,
CL_DEVICE_PREFERRED_VECTOR_WIDTH_DOUBLE,
sizeof(double_width),&double_width,nullptr));
op.preferred_vector_width64=double_width;
if (double_width==0)
op.double_precision=false;
else
// Determine if double precision is supported: All bits in the mask must be set.
cl_device_fp_config double_mask = (CL_FP_FMA|CL_FP_ROUND_TO_NEAREST|CL_FP_ROUND_TO_ZERO|
CL_FP_ROUND_TO_INF|CL_FP_INF_NAN|CL_FP_DENORM);
cl_device_fp_config double_avail;
CL_SAFE_CALL(clGetDeviceInfo(device_list,CL_DEVICE_DOUBLE_FP_CONFIG,
sizeof(double_avail),&double_avail,nullptr));
if ((double_avail & double_mask) == double_mask)
op.double_precision=true;
else
op.double_precision=false;
CL_SAFE_CALL(clGetDeviceInfo(device_list,
CL_DEVICE_PROFILING_TIMER_RESOLUTION,

2
lib/gpu/lal_born_coul_long.cpp
View File

@ -34,7 +34,7 @@ BornCoulLongT::BornCoulLong() : BaseCharge<numtyp,acctyp>(),
}
template <class numtyp, class acctyp>
BornCoulLongT::~BornCoulLongT() {
BornCoulLongT::~BornCoulLong() {
clear();
}

2
lib/gpu/lal_born_coul_wolf.cpp
View File

@ -34,7 +34,7 @@ BornCoulWolfT::BornCoulWolf() : BaseCharge<numtyp,acctyp>(),
}
template <class numtyp, class acctyp>
BornCoulWolfT::~BornCoulWolfT() {
BornCoulWolfT::~BornCoulWolf() {
clear();
}

2
lib/gpu/lal_buck_coul_long.cpp
View File

@ -34,7 +34,7 @@ BuckCoulLongT::BuckCoulLong() : BaseCharge<numtyp,acctyp>(),
}
template <class numtyp, class acctyp>
BuckCoulLongT::~BuckCoulLongT() {
BuckCoulLongT::~BuckCoulLong() {
clear();
}

13
lib/gpu/lal_device.cpp
View File

@ -333,6 +333,12 @@ int DeviceT::init_device(MPI_Comm world, MPI_Comm replica, const int ngpu,
gpu_barrier();
}
// check if double precision support is available
#if defined(_SINGLE_DOUBLE) || defined(_DOUBLE_DOUBLE)
if (!gpu->double_precision())
return -16;
#endif
// Setup auto bin size calculation for calls from atom::sort
// - This is repeated in neighbor init with additional info
if (_user_cell_size<0.0) {
@ -546,14 +552,9 @@ int DeviceT::init_nbor(Neighbor *nbor, const int nlocal,
return -3;
if (_user_cell_size<0.0) {
#ifndef LAL_USE_OLD_NEIGHBOR
_neighbor_shared.setup_auto_cell_size(true,cutoff,nbor->simd_size());
#else
_neighbor_shared.setup_auto_cell_size(false,cutoff,nbor->simd_size());
#endif
} else
_neighbor_shared.setup_auto_cell_size(false,_user_cell_size,
nbor->simd_size());
_neighbor_shared.setup_auto_cell_size(false,_user_cell_size,nbor->simd_size());
nbor->set_cutoff(cutoff);
return 0;

7
lib/machdyn/Install.py
View File

@ -17,11 +17,12 @@ parser = ArgumentParser(prog='Install.py',
# settings
version = '3.3.9'
version = '3.4.0'
tarball = "eigen.tar.gz"
# known checksums for different Eigen versions. used to validate the download.
checksums = { \
'3.4.0' : '4c527a9171d71a72a9d4186e65bea559', \
'3.3.9' : '609286804b0f79be622ccf7f9ff2b660', \
'3.3.7' : '9e30f67e8531477de4117506fe44669b' \
}
@ -35,7 +36,7 @@ Syntax from src dir: make lib-smd args="-b"
Syntax from lib dir: python Install.py -b
or: python Install.py -p /usr/include/eigen3"
or: python Install.py -v 3.3.7 -b
or: python Install.py -v 3.4.0 -b
Example:
@ -77,7 +78,7 @@ if pathflag:
if buildflag:
print("Downloading Eigen ...")
eigentar = os.path.join(homepath, tarball)
url = "https://gitlab.com/libeigen/eigen/-/archive/%s/eigen-%s.tar.gz" % (version,version)
url = "https://download.lammps.org/thirdparty/eigen-%s.tar.gz" % version
geturl(url, eigentar)
# verify downloaded archive integrity via md5 checksum, if known.

6
lib/pace/Makefile
View File

@ -2,8 +2,8 @@ SHELL = /bin/sh
# ------ FILES ------
SRC_FILES = $(wildcard src/ML-PACE/*.cpp)
SRC = $(filter-out src/ML-PACE/pair_pace.cpp, $(SRC_FILES))
SRC_FILES = $(wildcard src/USER-PACE/*.cpp)
SRC = $(filter-out src/USER-PACE/pair_pace.cpp, $(SRC_FILES))
# ------ DEFINITIONS ------
@ -12,7 +12,7 @@ OBJ = $(SRC:.cpp=.o)
# ------ SETTINGS ------
CXXFLAGS = -O3 -fPIC -Isrc/ML-PACE
CXXFLAGS = -O3 -fPIC -Isrc/USER-PACE
ARCHIVE = ar
ARCHFLAG = -rc

2
lib/pace/Makefile.lammps
View File

@ -1,3 +1,3 @@
pace_SYSINC =-I../../lib/pace/src/ML-PACE
pace_SYSINC =-I../../lib/pace/src/USER-PACE
pace_SYSLIB = -L../../lib/pace/ -lpace
pace_SYSPATH =

26
src/.gitignore vendored
View File

@ -27,6 +27,9 @@
/*_ssa.h
/*_ssa.cpp
!accelerator_kokkos.h
!accelerator_omp.h
/fix_mdi_engine.cpp
/fix_mdi_engine.h
/library_mdi.cpp
@ -202,7 +205,6 @@
/plugin.cpp
/plugin.h
/lammpsplugin.h
/atom_vec_spin.cpp
/atom_vec_spin.h
@ -265,8 +267,6 @@
/fix_drag.h
/fix_numdiff.cpp
/fix_numdiff.h
/fix_nve_noforce.cpp
/fix_nve_noforce.h
/fix_spring_rg.cpp
/fix_spring_rg.h
/fix_temp_csld.cpp
@ -367,8 +367,6 @@
/atom_vec_dpd.h
/atom_vec_electron.cpp
/atom_vec_electron.h
/atom_vec_ellipsoid.cpp
/atom_vec_ellipsoid.h
/atom_vec_full.cpp
/atom_vec_full.h
/atom_vec_full_hars.cpp
@ -535,8 +533,6 @@
/dihedral_harmonic.h
/dihedral_helix.cpp
/dihedral_helix.h
/dihedral_hybrid.cpp
/dihedral_hybrid.h
/dihedral_multi_harmonic.cpp
/dihedral_multi_harmonic.h
/dihedral_nharmonic.cpp
@ -858,8 +854,6 @@
/fix_ti_rs.h
/fix_ti_spring.cpp
/fix_ti_spring.h
/fix_ttm.cpp
/fix_ttm.h
/fix_tune_kspace.cpp
/fix_tune_kspace.h
/fix_wall_body_polygon.cpp
@ -885,8 +879,6 @@
/fix_widom.cpp
/fix_widom.h
/gpu_extra.h
/gridcomm.cpp
/gridcomm.h
/group_ndx.cpp
/group_ndx.h
/gz_file_writer.cpp
@ -911,14 +903,13 @@
/improper_fourier.h
/improper_harmonic.cpp
/improper_harmonic.h
/improper_hybrid.cpp
/improper_hybrid.h
/improper_inversion_harmonic.cpp
/improper_inversion_harmonic.h
/improper_ring.cpp
/improper_ring.h
/improper_umbrella.cpp
/improper_umbrella.h
/interlayer_taper.h
/kissfft.h
/lj_sdk_common.h
/math_complex.h
@ -933,7 +924,6 @@
/msm_cg.h
/neb.cpp
/neb.h
/pair_adp.cpp
/pair_adp.h
/pair_agni.cpp
@ -994,6 +984,8 @@
/pair_cosine_squared.h
/pair_coul_diel.cpp
/pair_coul_diel.h
/pair_coul_exclude.cpp
/pair_coul_exclude.h
/pair_coul_long.cpp
/pair_coul_long.h
/pair_coul_msm.cpp
@ -1332,8 +1324,6 @@
/thr_data.h
/verlet_split.cpp
/verlet_split.h
/write_dump.cpp
/write_dump.h
/xdr_compat.cpp
/xdr_compat.h
/zstd_file_writer.cpp
@ -1431,6 +1421,10 @@
/fix_srp.h
/fix_tfmc.cpp
/fix_tfmc.h
/fix_ttm.cpp
/fix_ttm.h
/fix_ttm_grid.cpp
/fix_ttm_grid.h
/fix_ttm_mod.cpp
/fix_ttm_mod.h
/pair_born_coul_long_cs.cpp

4
src/BOCS/fix_bocs.cpp
View File

@ -233,9 +233,7 @@ FixBocs::FixBocs(LAMMPS *lmp, int narg, char **arg) :
iarg += 2;
} else if (strcmp(arg[iarg],"mtk") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bocs command");
if (strcmp(arg[iarg+1],"yes") == 0) mtk_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) mtk_flag = 0;
else error->all(FLERR,"Illegal fix bocs command");
mtk_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"tloop") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bocs command");

34
src/COLVARS/fix_colvars.cpp
View File

@ -312,32 +312,26 @@ FixColvars::FixColvars(LAMMPS *lmp, int narg, char **arg) :
tmp_name = nullptr;
/* parse optional arguments */
int argsdone = 4;
while (argsdone < narg) {
int iarg = 4;
while (iarg < narg) {
// we have keyword/value pairs. check if value is missing
if (argsdone+1 == narg)
if (iarg+1 == narg)
error->all(FLERR,"Missing argument to keyword");
if (0 == strcmp(arg[argsdone], "input")) {
inp_name = strdup(arg[argsdone+1]);
} else if (0 == strcmp(arg[argsdone], "output")) {
out_name = strdup(arg[argsdone+1]);
} else if (0 == strcmp(arg[argsdone], "seed")) {
rng_seed = utils::inumeric(FLERR,arg[argsdone+1],false,lmp);
} else if (0 == strcmp(arg[argsdone], "unwrap")) {
if (0 == strcmp(arg[argsdone+1], "yes")) {
unwrap_flag = 1;
} else if (0 == strcmp(arg[argsdone+1], "no")) {
unwrap_flag = 0;
} else {
error->all(FLERR,"Incorrect fix colvars unwrap flag");
}
} else if (0 == strcmp(arg[argsdone], "tstat")) {
tmp_name = strdup(arg[argsdone+1]);
if (0 == strcmp(arg[iarg], "input")) {
inp_name = strdup(arg[iarg+1]);
} else if (0 == strcmp(arg[iarg], "output")) {
out_name = strdup(arg[iarg+1]);
} else if (0 == strcmp(arg[iarg], "seed")) {
rng_seed = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
} else if (0 == strcmp(arg[iarg], "unwrap")) {
unwrap_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
} else if (0 == strcmp(arg[iarg], "tstat")) {
tmp_name = strdup(arg[iarg+1]);
} else {
error->all(FLERR,"Unknown fix colvars parameter");
}
++argsdone; ++argsdone;
++iarg; ++iarg;
}
if (!out_name) out_name = strdup("out");

10
src/COMPRESS/dump_atom_zstd.cpp
View File

@ -188,17 +188,11 @@ int DumpAtomZstd::modify_param(int narg, char **arg)
try {
if (strcmp(arg[0], "checksum") == 0) {
if (narg < 2) error->all(FLERR, "Illegal dump_modify command");
if (strcmp(arg[1], "yes") == 0)
writer.setChecksum(true);
else if (strcmp(arg[1], "no") == 0)
writer.setChecksum(false);
else
error->all(FLERR, "Illegal dump_modify command");
writer.setChecksum(utils::logical(FLERR, arg[1], false, lmp) == 1);
return 2;
} else if (strcmp(arg[0], "compression_level") == 0) {
if (narg < 2) error->all(FLERR, "Illegal dump_modify command");
int compression_level = utils::inumeric(FLERR, arg[1], false, lmp);
writer.setCompressionLevel(compression_level);
writer.setCompressionLevel(utils::inumeric(FLERR, arg[1], false, lmp));
return 2;
}
} catch (FileWriterException &e) {

10
src/COMPRESS/dump_cfg_zstd.cpp
View File

@ -233,17 +233,11 @@ int DumpCFGZstd::modify_param(int narg, char **arg)
try {
if (strcmp(arg[0], "checksum") == 0) {
if (narg < 2) error->all(FLERR, "Illegal dump_modify command");
if (strcmp(arg[1], "yes") == 0)
writer.setChecksum(true);
else if (strcmp(arg[1], "no") == 0)
writer.setChecksum(false);
else
error->all(FLERR, "Illegal dump_modify command");
writer.setChecksum(utils::logical(FLERR, arg[1], false, lmp) == 1);
return 2;
} else if (strcmp(arg[0], "compression_level") == 0) {
if (narg < 2) error->all(FLERR, "Illegal dump_modify command");
int compression_level = utils::inumeric(FLERR, arg[1], false, lmp);
writer.setCompressionLevel(compression_level);
writer.setCompressionLevel(utils::inumeric(FLERR, arg[1], false, lmp));
return 2;
}
} catch (FileWriterException &e) {

15
src/COMPRESS/dump_custom_zstd.cpp
View File

@ -203,16 +203,13 @@ int DumpCustomZstd::modify_param(int narg, char **arg)
int consumed = DumpCustom::modify_param(narg, arg);
if (consumed == 0) {
try {
if (strcmp(arg[0],"checksum") == 0) {
if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
if (strcmp(arg[1],"yes") == 0) writer.setChecksum(true);
else if (strcmp(arg[1],"no") == 0) writer.setChecksum(false);
else error->all(FLERR,"Illegal dump_modify command");
if (strcmp(arg[0], "checksum") == 0) {
if (narg < 2) error->all(FLERR, "Illegal dump_modify command");
writer.setChecksum(utils::logical(FLERR, arg[1], false, lmp) == 1);
return 2;
} else if (strcmp(arg[0],"compression_level") == 0) {
if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
int compression_level = utils::inumeric(FLERR, arg[1], false, lmp);
writer.setCompressionLevel(compression_level);
} else if (strcmp(arg[0], "compression_level") == 0) {
if (narg < 2) error->all(FLERR, "Illegal dump_modify command");
writer.setCompressionLevel(utils::inumeric(FLERR, arg[1], false, lmp));
return 2;
}
} catch (FileWriterException &e) {

10
src/COMPRESS/dump_local_zstd.cpp
View File

@ -190,17 +190,11 @@ int DumpLocalZstd::modify_param(int narg, char **arg)
try {
if (strcmp(arg[0], "checksum") == 0) {
if (narg < 2) error->all(FLERR, "Illegal dump_modify command");
if (strcmp(arg[1], "yes") == 0)
writer.setChecksum(true);
else if (strcmp(arg[1], "no") == 0)
writer.setChecksum(false);
else
error->all(FLERR, "Illegal dump_modify command");
writer.setChecksum(utils::logical(FLERR, arg[1], false, lmp) == 1);
return 2;
} else if (strcmp(arg[0], "compression_level") == 0) {
if (narg < 2) error->all(FLERR, "Illegal dump_modify command");
int compression_level = utils::inumeric(FLERR, arg[1], false, lmp);
writer.setCompressionLevel(compression_level);
writer.setCompressionLevel(utils::inumeric(FLERR, arg[1], false, lmp));
return 2;
}
} catch (FileWriterException &e) {

10
src/COMPRESS/dump_xyz_zstd.cpp
View File

@ -156,17 +156,11 @@ int DumpXYZZstd::modify_param(int narg, char **arg)
try {
if (strcmp(arg[0], "checksum") == 0) {
if (narg < 2) error->all(FLERR, "Illegal dump_modify command");
if (strcmp(arg[1], "yes") == 0)
writer.setChecksum(true);
else if (strcmp(arg[1], "no") == 0)
writer.setChecksum(false);
else
error->all(FLERR, "Illegal dump_modify command");
writer.setChecksum(utils::logical(FLERR, arg[1], false, lmp) == 1);
return 2;
} else if (strcmp(arg[0], "compression_level") == 0) {
if (narg < 2) error->all(FLERR, "Illegal dump_modify command");
int compression_level = utils::inumeric(FLERR, arg[1], false, lmp);
writer.setCompressionLevel(compression_level);
writer.setCompressionLevel(utils::inumeric(FLERR, arg[1], false, lmp));
return 2;
}
} catch (FileWriterException &e) {

7
src/DIELECTRIC/fix_polarize_bem_gmres.cpp
View File

@ -796,12 +796,7 @@ int FixPolarizeBEMGMRES::modify_param(int narg, char **arg)
iarg += 2;
} else if (strcmp(arg[iarg], "kspace") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix_modify command");
if (strcmp(arg[iarg + 1], "yes") == 0)
kspaceflag = 1;
else if (strcmp(arg[iarg + 1], "no") == 0)
kspaceflag = 0;
else
error->all(FLERR, "Illegal fix_modify command for fix polarize");
kspaceflag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "dielectrics") == 0) {
if (iarg + 6 > narg) error->all(FLERR, "Illegal fix_modify command");

7
src/DIELECTRIC/fix_polarize_bem_icc.cpp
View File

@ -355,12 +355,7 @@ int FixPolarizeBEMICC::modify_param(int narg, char **arg)
iarg += 2;
} else if (strcmp(arg[iarg], "kspace") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix_modify command");
if (strcmp(arg[iarg + 1], "yes") == 0)
kspaceflag = 1;
else if (strcmp(arg[iarg + 1], "no") == 0)
kspaceflag = 0;
else
error->all(FLERR, "Illegal fix_modify command for fix polarize");
kspaceflag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "dielectrics") == 0) {
if (iarg + 6 > narg) error->all(FLERR, "Illegal fix_modify command");

7
src/DIELECTRIC/fix_polarize_functional.cpp
View File

@ -478,12 +478,7 @@ int FixPolarizeFunctional::modify_param(int narg, char **arg)
while (iarg < narg) {
if (strcmp(arg[iarg], "kspace") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix_modify command");
if (strcmp(arg[iarg + 1], "yes") == 0)
kspaceflag = 1;
else if (strcmp(arg[iarg + 1], "no") == 0)
kspaceflag = 0;
else
error->all(FLERR, "Illegal fix_modify command for fix polarize/functional");
kspaceflag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "dielectrics") == 0) {
if (iarg + 6 > narg) error->all(FLERR, "Illegal fix_modify command");

4
src/DPD-REACT/random_external_state.h
View File

@ -78,8 +78,8 @@
namespace random_external_state {
typedef uint64_t es_RNG_t;
enum { MAX_URAND = 0xffffffffU };
enum { MAX_URAND64 = 0xffffffffffffffffULL - 1 };
constexpr uint32_t MAX_URAND = 0xffffffffU;
constexpr uint64_t MAX_URAND64 = 0xffffffffffffffffULL - 1;
LAMMPS_INLINE
uint32_t es_urand(es_RNG_t &state_)

14
src/DRUDE/fix_drude_transform.cpp
View File

@ -13,16 +13,18 @@
------------------------------------------------------------------------- */
/** Fix Drude Transform ******************************************************/
#include "fix_drude_transform.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "fix_drude.h"
#include "modify.h"
#include <cmath>
#include <cstring>
#include "fix_drude.h"
#include "atom.h"
#include "domain.h"
#include "comm.h"
#include "error.h"
#include "modify.h"
using namespace LAMMPS_NS;
using namespace FixConst;

6
src/DRUDE/fix_drude_transform.h
View File

@ -25,10 +25,10 @@ FixStyle(drude/transform/inverse,FixDrudeTransform<true>);
namespace LAMMPS_NS {
template <bool inverse> class FixDrudeTransform : public Fix {
template <bool inverse> class FixDrudeTransform: public Fix {
public:
FixDrudeTransform<inverse>(class LAMMPS *, int, char **);
~FixDrudeTransform<inverse>();
FixDrudeTransform(class LAMMPS *, int, char **);
~FixDrudeTransform();
int setmask();
void init();
void setup(int vflag);

4
src/DRUDE/fix_langevin_drude.cpp
View File

@ -91,9 +91,7 @@ FixLangevinDrude::FixLangevinDrude(LAMMPS *lmp, int narg, char **arg) :
while (iarg < narg) {
if (strcmp(arg[iarg],"zero") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin/drude command");
if (strcmp(arg[iarg+1],"no") == 0) zero = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) zero = 1;
else error->all(FLERR,"Illegal fix langevin/drude command");
zero = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else error->all(FLERR,"Illegal fix langevin/drude command");
}

20
src/DRUDE/fix_tgnh_drude.cpp
View File

@ -256,9 +256,7 @@ FixTGNHDrude::FixTGNHDrude(LAMMPS *lmp, int narg, char **arg) :
iarg += 2;
} else if (strcmp(arg[iarg],"mtk") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
if (strcmp(arg[iarg+1],"yes") == 0) mtk_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) mtk_flag = 0;
else error->all(FLERR,"Illegal fix nvt/npt/nph command");
mtk_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"tloop") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
@ -277,27 +275,19 @@ FixTGNHDrude::FixTGNHDrude(LAMMPS *lmp, int narg, char **arg) :
iarg += 2;
} else if (strcmp(arg[iarg],"scalexy") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
if (strcmp(arg[iarg+1],"yes") == 0) scalexy = 1;
else if (strcmp(arg[iarg+1],"no") == 0) scalexy = 0;
else error->all(FLERR,"Illegal fix nvt/npt/nph command");
scalexy = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"scalexz") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
if (strcmp(arg[iarg+1],"yes") == 0) scalexz = 1;
else if (strcmp(arg[iarg+1],"no") == 0) scalexz = 0;
else error->all(FLERR,"Illegal fix nvt/npt/nph command");
scalexz = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"scaleyz") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
if (strcmp(arg[iarg+1],"yes") == 0) scaleyz = 1;
else if (strcmp(arg[iarg+1],"no") == 0) scaleyz = 0;
else error->all(FLERR,"Illegal fix nvt/npt/nph command");
scaleyz = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"flip") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
if (strcmp(arg[iarg+1],"yes") == 0) flipflag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) flipflag = 0;
else error->all(FLERR,"Illegal fix nvt/npt/nph command");
flipflag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"fixedpoint") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");

30
src/EXTRA-COMPUTE/compute_ackland_atom.cpp
View File

@ -18,20 +18,22 @@
Updated algorithm by: Brian Barnes, brian.c.barnes11.civ@mail.mil
------------------------------------------------------------------------- */
#include <cmath>
#include <cstring>
#include "compute_ackland_atom.h"
#include "atom.h"
#include "update.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "force.h"
#include "neighbor.h"
#include "pair.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
@ -60,16 +62,10 @@ ComputeAcklandAtom::ComputeAcklandAtom(LAMMPS *lmp, int narg, char **arg) :
int iarg = 3;
while (narg > iarg) {
if (strcmp("legacy",arg[iarg]) == 0) {
++iarg;
if (iarg >= narg)
error->all(FLERR,"Invalid compute ackland/atom command");
if (strcmp("yes",arg[iarg]) == 0)
legacy = 1;
else if (strcmp("no",arg[iarg]) == 0)
legacy = 0;
else error->all(FLERR,"Invalid compute ackland/atom command");
if (iarg+2 > narg) error->all(FLERR,"Invalid compute ackland/atom command");
legacy = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
}
++iarg;
}
}

54
src/EXTRA-COMPUTE/compute_entropy_atom.cpp
View File

@ -17,21 +17,23 @@
------------------------------------------------------------------------- */
#include "compute_entropy_atom.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "force.h"
#include "pair.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "domain.h"
#include "modify.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include "pair.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
@ -44,8 +46,7 @@ ComputeEntropyAtom(LAMMPS *lmp, int narg, char **arg) :
pair_entropy(nullptr), pair_entropy_avg(nullptr)
{
if (narg < 5 || narg > 10)
error->all(FLERR,"Illegal compute entropy/atom command; wrong number"
" of arguments");
error->all(FLERR,"Illegal compute entropy/atom command; wrong number of arguments");
// Arguments are: sigma cutoff avg yes/no cutoff2 local yes/no
// sigma is the gaussian width
@ -57,11 +58,9 @@ ComputeEntropyAtom(LAMMPS *lmp, int narg, char **arg) :
// the g(r)
sigma = utils::numeric(FLERR,arg[3],false,lmp);
if (sigma <= 0.0) error->all(FLERR,"Illegal compute entropy/atom"
" command; sigma must be positive");
if (sigma <= 0.0) error->all(FLERR,"Illegal compute entropy/atom command; sigma must be positive");
cutoff = utils::numeric(FLERR,arg[4],false,lmp);
if (cutoff <= 0.0) error->all(FLERR,"Illegal compute entropy/atom"
" command; cutoff must be positive");
if (cutoff <= 0.0) error->all(FLERR,"Illegal compute entropy/atom command; cutoff must be positive");
cutoff2 = 0.;
avg_flag = 0;
@ -71,26 +70,17 @@ ComputeEntropyAtom(LAMMPS *lmp, int narg, char **arg) :
int iarg = 5;
while (iarg < narg) {
if (strcmp(arg[iarg],"avg") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute entropy/atom;"
" missing arguments after avg");
if (strcmp(arg[iarg+1],"yes") == 0) avg_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) avg_flag = 0;
else error->all(FLERR,"Illegal compute entropy/atom;"
" argument after avg should be yes or no");
if (iarg+3 > narg) error->all(FLERR,"Illegal compute entropy/atom command");
avg_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
cutoff2 = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (cutoff2 < 0.0) error->all(FLERR,"Illegal compute entropy/atom"
" command; negative cutoff2");
if (cutoff2 < 0.0) error->all(FLERR,"Illegal compute entropy/atom command; negative cutoff2");
cutsq2 = cutoff2*cutoff2;
iarg += 3;
} else if (strcmp(arg[iarg],"local") == 0) {
if (strcmp(arg[iarg+1],"yes") == 0) local_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) local_flag = 0;
else error->all(FLERR,"Illegal compute entropy/atom;"
" argument after local should be yes or no");
if (iarg+2 > narg) error->all(FLERR,"Illegal compute entropy/atom command");
local_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else error->all(FLERR,"Illegal compute entropy/atom; argument after"
" sigma and cutoff should be avg or local");
} else error->all(FLERR,"Illegal compute entropy/atom command");
}

4
src/EXTRA-DUMP/dump_dcd.cpp
View File

@ -258,9 +258,7 @@ int DumpDCD::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"unwrap") == 0) {
if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
if (strcmp(arg[1],"yes") == 0) unwrap_flag = 1;
else if (strcmp(arg[1],"no") == 0) unwrap_flag = 0;
else error->all(FLERR,"Illegal dump_modify command");
unwrap_flag = utils::logical(FLERR,arg[1],false,lmp);
return 2;
}
return 0;

23
src/EXTRA-DUMP/dump_xtc.cpp
View File

@ -24,19 +24,20 @@
------------------------------------------------------------------------- */
#include "dump_xtc.h"
#include <cmath>
#include <cstring>
#include <climits>
#include "domain.h"
#include "atom.h"
#include "update.h"
#include "group.h"
#include "output.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "memory.h"
#include "output.h"
#include "update.h"
#include <climits>
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
@ -278,9 +279,7 @@ int DumpXTC::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"unwrap") == 0) {
if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
if (strcmp(arg[1],"yes") == 0) unwrap_flag = 1;
else if (strcmp(arg[1],"no") == 0) unwrap_flag = 0;
else error->all(FLERR,"Illegal dump_modify command");
unwrap_flag = utils::logical(FLERR,arg[1],false,lmp);
return 2;
} else if (strcmp(arg[0],"precision") == 0) {
if (narg < 2) error->all(FLERR,"Illegal dump_modify command");

4
src/EXTRA-FIX/fix_flow_gauss.cpp
View File

@ -91,9 +91,7 @@ FixFlowGauss::FixFlowGauss(LAMMPS *lmp, int narg, char **arg) :
while (iarg < narg) {
if (strcmp(arg[iarg],"energy") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal energy keyword");
if (strcmp(arg[iarg+1],"yes") == 0) workflag = true;
else if (strcmp(arg[iarg+1],"no") != 0)
error->all(FLERR,"Illegal energy keyword");
workflag = utils::logical(FLERR,arg[iarg+1],false,lmp) == 1;
iarg += 2;
} else error->all(FLERR,"Illegal fix flow/gauss command");
}

43
src/EXTRA-FIX/fix_gld.cpp
View File

@ -19,17 +19,18 @@
#include "fix_gld.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "memory.h"
#include "random_mars.h"
#include "respa.h"
#include "update.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "force.h"
#include "update.h"
#include "respa.h"
#include "comm.h"
#include "random_mars.h"
#include "memory.h"
#include "error.h"
#include "group.h"
#define GLD_UNIFORM_DISTRO
@ -128,26 +129,14 @@ FixGLD::FixGLD(LAMMPS *lmp, int narg, char **arg) :
while (iarg < narg) {
if (strcmp(arg[iarg],"zero") == 0) {
if (iarg+2 > narg) {
error->all(FLERR, "Illegal fix gld command");
}
if (strcmp(arg[iarg+1],"no") == 0) {
zeroflag = 0;
} else if (strcmp(arg[iarg+1],"yes") == 0) {
zeroflag = 1;
} else {
error->all(FLERR,"Illegal fix gld command");
}
if (iarg+2 > narg) error->all(FLERR, "Illegal fix gld command");
zeroflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
}
else if (strcmp(arg[iarg],"frozen") == 0) {
if (iarg+2 > narg) {
error->all(FLERR, "Illegal fix gld command");
}
if (strcmp(arg[iarg+1],"no") == 0) {
freezeflag = 0;
} else if (strcmp(arg[iarg+1],"yes") == 0) {
freezeflag = 1;
if (iarg+2 > narg) error->all(FLERR, "Illegal fix gld command");
freezeflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
if (freezeflag) {
for (int i = 0; i < atom->nlocal; i++) {
if (atom->mask[i] & groupbit) {
for (int k = 0; k < 3*prony_terms; k=k+3)
@ -158,8 +147,6 @@ FixGLD::FixGLD(LAMMPS *lmp, int narg, char **arg) :
}
}
}
} else {
error->all(FLERR, "Illegal fix gld command");
}
iarg += 2;
}

30
src/EXTRA-FIX/fix_npt_cauchy.cpp
View File

@ -297,9 +297,7 @@ FixNPTCauchy::FixNPTCauchy(LAMMPS *lmp, int narg, char **arg) :
iarg += 2;
} else if (strcmp(arg[iarg],"mtk") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
if (strcmp(arg[iarg+1],"yes") == 0) mtk_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) mtk_flag = 0;
else error->all(FLERR,"Illegal fix npt/cauchy command");
mtk_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"tloop") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
@ -318,27 +316,19 @@ FixNPTCauchy::FixNPTCauchy(LAMMPS *lmp, int narg, char **arg) :
iarg += 2;
} else if (strcmp(arg[iarg],"scalexy") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
if (strcmp(arg[iarg+1],"yes") == 0) scalexy = 1;
else if (strcmp(arg[iarg+1],"no") == 0) scalexy = 0;
else error->all(FLERR,"Illegal fix npt/cauchy command");
scalexy = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"scalexz") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
if (strcmp(arg[iarg+1],"yes") == 0) scalexz = 1;
else if (strcmp(arg[iarg+1],"no") == 0) scalexz = 0;
else error->all(FLERR,"Illegal fix npt/cauchy command");
scalexz = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"scaleyz") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
if (strcmp(arg[iarg+1],"yes") == 0) scaleyz = 1;
else if (strcmp(arg[iarg+1],"no") == 0) scaleyz = 0;
else error->all(FLERR,"Illegal fix npt/cauchy command");
scaleyz = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"flip") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
if (strcmp(arg[iarg+1],"yes") == 0) flipflag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) flipflag = 0;
else error->all(FLERR,"Illegal fix npt/cauchy command");
flipflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"update") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
@ -352,10 +342,8 @@ FixNPTCauchy::FixNPTCauchy(LAMMPS *lmp, int narg, char **arg) :
alpha = utils::numeric(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"continue") == 0) {
if (strcmp(arg[iarg+1],"yes") != 0 && strcmp(arg[iarg+1],"no") != 0)
error->all(FLERR,"Illegal cauchystat continue value. "
"Must be 'yes' or 'no'");
restartPK = !strcmp(arg[iarg+1],"yes");
if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
restartPK = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"fixedpoint") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
@ -2454,7 +2442,7 @@ double FixNPTCauchy::memory_usage()
void FixNPTCauchy::CauchyStat_init()
{
if (comm->me == 0) {
std::string mesg = fmt::format("Using fix npt/cauchy with alpha={:f.8}\n",alpha);
std::string mesg = fmt::format("Using fix npt/cauchy with alpha={:.8f}\n",alpha);
if (restartPK==1) {
mesg += " (this is a continuation run)\n";
} else {
@ -2475,7 +2463,7 @@ void FixNPTCauchy::CauchyStat_init()
error->all(FLERR,"Illegal fix npt/cauchy command: Alpha cannot be zero or negative.");
if (restart_stored < 0) {
modify->add_fix(std::string(id_store) + "all STORE global 1 6");
modify->add_fix(std::string(id_store) + " all STORE global 1 6");
restart_stored = modify->find_fix(id_store);
}
init_store = (FixStore *)modify->fix[restart_stored];

12
src/EXTRA-FIX/fix_pafi.cpp
View File

@ -103,16 +103,12 @@ FixPAFI::FixPAFI(LAMMPS *lmp, int narg, char **arg) :
int iarg = 7;
while (iarg < narg) {
if (strcmp(arg[iarg],"overdamped") == 0) {
if (strcmp(arg[iarg+1],"no") == 0) od_flag = 0;
else if (strcmp(arg[iarg+1],"0") == 0) od_flag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) od_flag = 1;
else if (strcmp(arg[iarg+1],"1") == 0) od_flag = 1;
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pafi command");
od_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"com") == 0) {
if (strcmp(arg[iarg+1],"no") == 0) com_flag = 0;
else if (strcmp(arg[iarg+1],"0") == 0) com_flag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) com_flag = 1;
else if (strcmp(arg[iarg+1],"1") == 0) com_flag = 1;
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pafi command");
com_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else error->all(FLERR,"Illegal fix pafi command");
}

189
src/EXTRA-MOLECULE/angle_gaussian.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
@ -34,9 +33,9 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */
AngleGaussian::AngleGaussian(LAMMPS *lmp)
: Angle(lmp), nterms(nullptr), angle_temperature(nullptr),
alpha(nullptr), width(nullptr), theta0(nullptr)
AngleGaussian::AngleGaussian(LAMMPS *lmp) :
Angle(lmp), nterms(nullptr), angle_temperature(nullptr), alpha(nullptr), width(nullptr),
theta0(nullptr)
{
}
@ -49,13 +48,13 @@ AngleGaussian::~AngleGaussian()
memory->destroy(nterms);
memory->destroy(angle_temperature);
for (int i = 1; i <= atom->nangletypes; i++) {
delete [] alpha[i];
delete [] width[i];
delete [] theta0[i];
delete[] alpha[i];
delete[] width[i];
delete[] theta0[i];
}
delete [] alpha;
delete [] width;
delete [] theta0;
delete[] alpha;
delete[] width;
delete[] theta0;
}
}
@ -63,15 +62,15 @@ AngleGaussian::~AngleGaussian()
void AngleGaussian::compute(int eflag, int vflag)
{
int i1,i2,i3,n,type;
double delx1,dely1,delz1,delx2,dely2,delz2;
double eangle,f1[3],f3[3];
int i1, i2, i3, n, type;
double delx1, dely1, delz1, delx2, dely2, delz2;
double eangle, f1[3], f3[3];
double dtheta;
double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22;
double rsq1, rsq2, r1, r2, c, s, a, a11, a12, a22;
double prefactor, exponent, g_i, sum_g_i, sum_numerator;
eangle = 0.0;
ev_init(eflag,vflag);
ev_init(eflag, vflag);
double **x = atom->x;
double **f = atom->f;
@ -92,7 +91,7 @@ void AngleGaussian::compute(int eflag, int vflag)
dely1 = x[i1][1] - x[i2][1];
delz1 = x[i1][2] - x[i2][2];
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
rsq1 = delx1 * delx1 + dely1 * dely1 + delz1 * delz1;
r1 = sqrt(rsq1);
// 2nd bond
@ -101,20 +100,20 @@ void AngleGaussian::compute(int eflag, int vflag)
dely2 = x[i3][1] - x[i2][1];
delz2 = x[i3][2] - x[i2][2];
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
rsq2 = delx2 * delx2 + dely2 * dely2 + delz2 * delz2;
r2 = sqrt(rsq2);
// angle (cos and sin)
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
c /= r1*r2;
c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
c /= r1 * r2;
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
s = sqrt(1.0 - c*c);
s = sqrt(1.0 - c * c);
if (s < SMAL) s = SMAL;
s = 1.0/s;
s = 1.0 / s;
// force & energy
double theta = acos(c);
@ -123,28 +122,28 @@ void AngleGaussian::compute(int eflag, int vflag)
sum_numerator = 0.0;
for (int i = 0; i < nterms[type]; i++) {
dtheta = theta - theta0[type][i];
prefactor = (alpha[type][i]/(width[type][i]*sqrt(MY_PI2)));
exponent = -2*dtheta*dtheta/(width[type][i]*width[type][i]);
g_i = prefactor*exp(exponent);
prefactor = (alpha[type][i] / (width[type][i] * sqrt(MY_PI2)));
exponent = -2 * dtheta * dtheta / (width[type][i] * width[type][i]);
g_i = prefactor * exp(exponent);
sum_g_i += g_i;
sum_numerator += g_i*dtheta/(width[type][i]*width[type][i]);
sum_numerator += g_i * dtheta / (width[type][i] * width[type][i]);
}
if (sum_g_i < SMALL) sum_g_i = SMALL;
if (eflag) eangle = -(force->boltz*angle_temperature[type])*log(sum_g_i);
if (eflag) eangle = -(force->boltz * angle_temperature[type]) * log(sum_g_i);
// I should check about the sign of this expression
a = -4.0*(force->boltz*angle_temperature[type])*(sum_numerator/sum_g_i)*s;
a11 = a*c / rsq1;
a12 = -a / (r1*r2);
a22 = a*c / rsq2;
a = -4.0 * (force->boltz * angle_temperature[type]) * (sum_numerator / sum_g_i) * s;
a11 = a * c / rsq1;
a12 = -a / (r1 * r2);
a22 = a * c / rsq2;
f1[0] = a11*delx1 + a12*delx2;
f1[1] = a11*dely1 + a12*dely2;
f1[2] = a11*delz1 + a12*delz2;
f3[0] = a22*delx2 + a12*delx1;
f3[1] = a22*dely2 + a12*dely1;
f3[2] = a22*delz2 + a12*delz1;
f1[0] = a11 * delx1 + a12 * delx2;
f1[1] = a11 * dely1 + a12 * dely2;
f1[2] = a11 * delz1 + a12 * delz2;
f3[0] = a22 * delx2 + a12 * delx1;
f3[1] = a22 * dely2 + a12 * dely1;
f3[2] = a22 * delz2 + a12 * delz1;
// apply force to each of 3 atoms
@ -166,8 +165,9 @@ void AngleGaussian::compute(int eflag, int vflag)
f[i3][2] += f3[2];
}
if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
delx1,dely1,delz1,delx2,dely2,delz2);
if (evflag)
ev_tally(i1, i2, i3, nlocal, newton_bond, eangle, f1, f3, delx1, dely1, delz1, delx2, dely2,
delz2);
}
}
@ -176,20 +176,20 @@ void AngleGaussian::compute(int eflag, int vflag)
void AngleGaussian::allocate()
{
allocated = 1;
int n = atom->nangletypes;
int n = atom->nangletypes + 1;
memory->create(nterms,n+1,"angle:nterms");
memory->create(angle_temperature,n+1,"angle:angle_temperature");
memory->create(nterms, n, "angle:nterms");
memory->create(angle_temperature, n, "angle:angle_temperature");
alpha = new double *[n+1];
width = new double *[n+1];
theta0 = new double *[n+1];
memset(alpha,0,sizeof(double)*(n+1));
memset(width,0,sizeof(double)*(n+1));
memset(theta0,0,sizeof(double)*(n+1));
alpha = new double *[n];
width = new double *[n];
theta0 = new double *[n];
memset(alpha, 0, sizeof(double *) * n);
memset(width, 0, sizeof(double *) * n);
memset(theta0, 0, sizeof(double *) * n);
memory->create(setflag,n+1,"angle:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0;
memory->create(setflag, n, "angle:setflag");
memset(setflag, 0, sizeof(int) * n);
}
/* ----------------------------------------------------------------------
@ -198,15 +198,14 @@ void AngleGaussian::allocate()
void AngleGaussian::coeff(int narg, char **arg)
{
if (narg < 6) error->all(FLERR,"Incorrect args for angle coefficients");
if (narg < 6) error->all(FLERR, "Incorrect args for angle coefficients");
int ilo,ihi;
utils::bounds(FLERR,arg[0],1,atom->nangletypes,ilo,ihi,error);
int ilo, ihi;
utils::bounds(FLERR, arg[0], 1, atom->nangletypes, ilo, ihi, error);
double angle_temperature_one = utils::numeric(FLERR,arg[1],false,lmp);
int n = utils::inumeric(FLERR,arg[2],false,lmp);
if (narg != 3*n + 3)
error->all(FLERR,"Incorrect args for angle coefficients");
double angle_temperature_one = utils::numeric(FLERR, arg[1], false, lmp);
int n = utils::inumeric(FLERR, arg[2], false, lmp);
if (narg != 3 * n + 3) error->all(FLERR, "Incorrect args for angle coefficients");
if (!allocated) allocate();
@ -217,21 +216,21 @@ void AngleGaussian::coeff(int narg, char **arg)
angle_temperature[i] = angle_temperature_one;
nterms[i] = n;
delete[] alpha[i];
alpha[i] = new double [n];
alpha[i] = new double[n];
delete[] width[i];
width[i] = new double [n];
width[i] = new double[n];
delete[] theta0[i];
theta0[i] = new double [n];
theta0[i] = new double[n];
for (int j = 0; j < n; j++) {
alpha[i][j] = utils::numeric(FLERR,arg[3+3*j],false,lmp);
width[i][j] = utils::numeric(FLERR,arg[4+3*j],false,lmp);
theta0[i][j] = utils::numeric(FLERR,arg[5+3*j],false,lmp)* MY_PI / 180.0;
alpha[i][j] = utils::numeric(FLERR, arg[3 + 3 * j], false, lmp);
width[i][j] = utils::numeric(FLERR, arg[4 + 3 * j], false, lmp);
theta0[i][j] = utils::numeric(FLERR, arg[5 + 3 * j], false, lmp) * MY_PI / 180.0;
setflag[i] = 1;
}
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
if (count == 0) error->all(FLERR, "Incorrect args for angle coefficients");
}
/* ---------------------------------------------------------------------- */
@ -247,12 +246,12 @@ double AngleGaussian::equilibrium_angle(int i)
void AngleGaussian::write_restart(FILE *fp)
{
fwrite(&angle_temperature[1],sizeof(double),atom->nangletypes,fp);
fwrite(&nterms[1],sizeof(int),atom->nangletypes,fp);
fwrite(&angle_temperature[1], sizeof(double), atom->nangletypes, fp);
fwrite(&nterms[1], sizeof(int), atom->nangletypes, fp);
for (int i = 1; i <= atom->nangletypes; i++) {
fwrite(alpha[i],sizeof(double),nterms[i],fp);
fwrite(width[i],sizeof(double),nterms[i],fp);
fwrite(theta0[i],sizeof(double),nterms[i],fp);
fwrite(alpha[i], sizeof(double), nterms[i], fp);
fwrite(width[i], sizeof(double), nterms[i], fp);
fwrite(theta0[i], sizeof(double), nterms[i], fp);
}
}
@ -265,31 +264,32 @@ void AngleGaussian::read_restart(FILE *fp)
allocate();
if (comm->me == 0) {
utils::sfread(FLERR,&angle_temperature[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
utils::sfread(FLERR,&nterms[1],sizeof(int),atom->nangletypes,fp,nullptr,error);
utils::sfread(FLERR, &angle_temperature[1], sizeof(double), atom->nangletypes, fp, nullptr,
error);
utils::sfread(FLERR, &nterms[1], sizeof(int), atom->nangletypes, fp, nullptr, error);
}
MPI_Bcast(&angle_temperature[1],atom->nangletypes,MPI_DOUBLE,0,world);
MPI_Bcast(&nterms[1],atom->nangletypes,MPI_INT,0,world);
MPI_Bcast(&angle_temperature[1], atom->nangletypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&nterms[1], atom->nangletypes, MPI_INT, 0, world);
// allocate
for (int i = 1; i <= atom->nangletypes; i++) {
alpha[i] = new double [nterms[i]];
width[i] = new double [nterms[i]];
theta0[i] = new double [nterms[i]];
alpha[i] = new double[nterms[i]];
width[i] = new double[nterms[i]];
theta0[i] = new double[nterms[i]];
}
if (comm->me == 0) {
for (int i = 1; i <= atom->nangletypes; i++) {
utils::sfread(FLERR,alpha[i],sizeof(double),nterms[i],fp,nullptr,error);
utils::sfread(FLERR,width[i],sizeof(double),nterms[i],fp,nullptr,error);
utils::sfread(FLERR,theta0[i],sizeof(double),nterms[i],fp,nullptr,error);
utils::sfread(FLERR, alpha[i], sizeof(double), nterms[i], fp, nullptr, error);
utils::sfread(FLERR, width[i], sizeof(double), nterms[i], fp, nullptr, error);
utils::sfread(FLERR, theta0[i], sizeof(double), nterms[i], fp, nullptr, error);
}
}
for (int i = 1; i <= atom->nangletypes; i++) {
MPI_Bcast(alpha[i],nterms[i],MPI_DOUBLE,0,world);
MPI_Bcast(width[i],nterms[i],MPI_DOUBLE,0,world);
MPI_Bcast(theta0[i],nterms[i],MPI_DOUBLE,0,world);
MPI_Bcast(alpha[i], nterms[i], MPI_DOUBLE, 0, world);
MPI_Bcast(width[i], nterms[i], MPI_DOUBLE, 0, world);
MPI_Bcast(theta0[i], nterms[i], MPI_DOUBLE, 0, world);
}
for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
@ -302,13 +302,12 @@ void AngleGaussian::read_restart(FILE *fp)
void AngleGaussian::write_data(FILE *fp)
{
for (int i = 1; i <= atom->nangletypes; i++) {
fprintf(fp,"%d %g %d",i,angle_temperature[i],nterms[i]);
fprintf(fp, "%d %g %d", i, angle_temperature[i], nterms[i]);
for (int j = 0; j < nterms[i]; j++) {
fprintf(fp," %g %g %g",alpha[i][j],width[i][j],(theta0[i][j]/MY_PI)*180.0);
fprintf(fp, " %g %g %g", alpha[i][j], width[i][j], (theta0[i][j] / MY_PI) * 180.0);
}
fprintf(fp, "\n");
}
}
/* ---------------------------------------------------------------------- */
@ -320,30 +319,30 @@ double AngleGaussian::single(int type, int i1, int i2, int i3)
double delx1 = x[i1][0] - x[i2][0];
double dely1 = x[i1][1] - x[i2][1];
double delz1 = x[i1][2] - x[i2][2];
domain->minimum_image(delx1,dely1,delz1);
double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
domain->minimum_image(delx1, dely1, delz1);
double r1 = sqrt(delx1 * delx1 + dely1 * dely1 + delz1 * delz1);
double delx2 = x[i3][0] - x[i2][0];
double dely2 = x[i3][1] - x[i2][1];
double delz2 = x[i3][2] - x[i2][2];
domain->minimum_image(delx2,dely2,delz2);
double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
domain->minimum_image(delx2, dely2, delz2);
double r2 = sqrt(delx2 * delx2 + dely2 * dely2 + delz2 * delz2);
double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
c /= r1*r2;
double c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
c /= r1 * r2;
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
double theta = acos(c) ;
double theta = acos(c);
double sum_g_i = 0.0;
for (int i = 0; i < nterms[type]; i++) {
double dtheta = theta - theta0[type][i];
double prefactor = (alpha[type][i]/(width[type][i]*sqrt(MY_PI2)));
double exponent = -2*dtheta*dtheta/(width[type][i]*width[type][i]);
double g_i = prefactor*exp(exponent);
double prefactor = (alpha[type][i] / (width[type][i] * sqrt(MY_PI2)));
double exponent = -2 * dtheta * dtheta / (width[type][i] * width[type][i]);
double g_i = prefactor * exp(exponent);
sum_g_i += g_i;
}
if (sum_g_i < SMALL) sum_g_i = SMALL;
return -(force->boltz*angle_temperature[type])*log(sum_g_i);
return -(force->boltz * angle_temperature[type]) * log(sum_g_i);
}

182
src/EXTRA-MOLECULE/bond_gaussian.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
@ -14,16 +13,16 @@
#include "bond_gaussian.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "neighbor.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
@ -32,9 +31,9 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */
BondGaussian::BondGaussian(LAMMPS *lmp)
: Bond(lmp), nterms(nullptr), bond_temperature(nullptr),
alpha(nullptr), width(nullptr), r0(nullptr)
BondGaussian::BondGaussian(LAMMPS *lmp) :
Bond(lmp), nterms(nullptr), bond_temperature(nullptr), alpha(nullptr), width(nullptr),
r0(nullptr)
{
reinitflag = 1;
}
@ -48,13 +47,13 @@ BondGaussian::~BondGaussian()
memory->destroy(nterms);
memory->destroy(bond_temperature);
for (int i = 1; i <= atom->nbondtypes; i++) {
delete [] alpha[i];
delete [] width[i];
delete [] r0[i];
delete[] alpha[i];
delete[] width[i];
delete[] r0[i];
}
delete [] alpha;
delete [] width;
delete [] r0;
delete[] alpha;
delete[] width;
delete[] r0;
}
}
@ -62,13 +61,13 @@ BondGaussian::~BondGaussian()
void BondGaussian::compute(int eflag, int vflag)
{
int i1,i2,n,type;
double delx,dely,delz,ebond,fbond;
double rsq,r,dr;
int i1, i2, n, type;
double delx, dely, delz, ebond, fbond;
double rsq, r, dr;
double prefactor, exponent, g_i, sum_g_i, sum_numerator;
ebond = 0.0;
ev_init(eflag,vflag);
ev_init(eflag, vflag);
double **x = atom->x;
double **f = atom->f;
@ -86,43 +85,45 @@ void BondGaussian::compute(int eflag, int vflag)
dely = x[i1][1] - x[i2][1];
delz = x[i1][2] - x[i2][2];
rsq = delx*delx + dely*dely + delz*delz;
rsq = delx * delx + dely * dely + delz * delz;
r = sqrt(rsq);
sum_g_i = 0.0;
sum_numerator = 0.0;
for (int i = 0; i < nterms[type]; i++) {
dr = r - r0[type][i];
prefactor = (alpha[type][i]/(width[type][i]*sqrt(MY_PI2)));
exponent = -2*dr*dr/(width[type][i]*width[type][i]);
g_i = prefactor*exp(exponent);
prefactor = (alpha[type][i] / (width[type][i] * sqrt(MY_PI2)));
exponent = -2 * dr * dr / (width[type][i] * width[type][i]);
g_i = prefactor * exp(exponent);
sum_g_i += g_i;
sum_numerator += g_i*dr/(width[type][i]*width[type][i]);
sum_numerator += g_i * dr / (width[type][i] * width[type][i]);
}
// force & energy
if (sum_g_i < SMALL) sum_g_i = SMALL;
if (r > 0.0) fbond = -4.0*(force->boltz*bond_temperature[type])*(sum_numerator/sum_g_i)/r;
else fbond = 0.0;
if (r > 0.0)
fbond = -4.0 * (force->boltz * bond_temperature[type]) * (sum_numerator / sum_g_i) / r;
else
fbond = 0.0;
if (eflag) ebond = -(force->boltz*bond_temperature[type])*log(sum_g_i);
if (eflag) ebond = -(force->boltz * bond_temperature[type]) * log(sum_g_i);
// apply force to each of 2 atoms
if (newton_bond || i1 < nlocal) {
f[i1][0] += delx*fbond;
f[i1][1] += dely*fbond;
f[i1][2] += delz*fbond;
f[i1][0] += delx * fbond;
f[i1][1] += dely * fbond;
f[i1][2] += delz * fbond;
}
if (newton_bond || i2 < nlocal) {
f[i2][0] -= delx*fbond;
f[i2][1] -= dely*fbond;
f[i2][2] -= delz*fbond;
f[i2][0] -= delx * fbond;
f[i2][1] -= dely * fbond;
f[i2][2] -= delz * fbond;
}
if (evflag) ev_tally(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz);
if (evflag) ev_tally(i1, i2, nlocal, newton_bond, ebond, fbond, delx, dely, delz);
}
}
@ -131,20 +132,20 @@ void BondGaussian::compute(int eflag, int vflag)
void BondGaussian::allocate()
{
allocated = 1;
int n = atom->nbondtypes;
int n = atom->nbondtypes + 1;
memory->create(nterms,n+1,"bond:nterms");
memory->create(bond_temperature,n+1,"bond:bond_temperature");
memory->create(nterms, n, "bond:nterms");
memory->create(bond_temperature, n, "bond:bond_temperature");
alpha = new double *[n+1];
width = new double *[n+1];
r0 = new double *[n+1];
memset(alpha,0,sizeof(double)*(n+1));
memset(width,0,sizeof(double)*(n+1));
memset(r0,0,sizeof(double)*(n+1));
alpha = new double *[n];
width = new double *[n];
r0 = new double *[n];
memset(alpha, 0, sizeof(double *) * n);
memset(width, 0, sizeof(double *) * n);
memset(r0, 0, sizeof(double *) * n);
memory->create(setflag,n+1,"bond:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0;
memory->create(setflag, n, "bond:setflag");
memset(setflag, 0, sizeof(int) * n);
}
/* ----------------------------------------------------------------------
@ -153,15 +154,14 @@ void BondGaussian::allocate()
void BondGaussian::coeff(int narg, char **arg)
{
if (narg < 6) error->all(FLERR,"Incorrect args for bond coefficients");
if (narg < 6) error->all(FLERR, "Incorrect args for bond coefficients");
int ilo,ihi;
utils::bounds(FLERR,arg[0],1,atom->nbondtypes,ilo,ihi,error);
int ilo, ihi;
utils::bounds(FLERR, arg[0], 1, atom->nbondtypes, ilo, ihi, error);
double bond_temp_one = utils::numeric(FLERR,arg[1],false,lmp);
int n = utils::inumeric(FLERR,arg[2],false,lmp);
if (narg != 3*n + 3)
error->all(FLERR,"Incorrect args for bond coefficients");
double bond_temp_one = utils::numeric(FLERR, arg[1], false, lmp);
int n = utils::inumeric(FLERR, arg[2], false, lmp);
if (narg != 3 * n + 3) error->all(FLERR, "Incorrect args for bond coefficients");
if (!allocated) allocate();
@ -170,21 +170,21 @@ void BondGaussian::coeff(int narg, char **arg)
bond_temperature[i] = bond_temp_one;
nterms[i] = n;
delete[] alpha[i];
alpha[i] = new double [n];
alpha[i] = new double[n];
delete[] width[i];
width[i] = new double [n];
width[i] = new double[n];
delete[] r0[i];
r0[i] = new double [n];
r0[i] = new double[n];
for (int j = 0; j < n; j++) {
alpha[i][j] = utils::numeric(FLERR,arg[3+3*j],false,lmp);
width[i][j] = utils::numeric(FLERR,arg[4+3*j],false,lmp);
r0[i][j] = utils::numeric(FLERR,arg[5+3*j],false,lmp);
alpha[i][j] = utils::numeric(FLERR, arg[3 + 3 * j], false, lmp);
width[i][j] = utils::numeric(FLERR, arg[4 + 3 * j], false, lmp);
r0[i][j] = utils::numeric(FLERR, arg[5 + 3 * j], false, lmp);
setflag[i] = 1;
}
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients");
if (count == 0) error->all(FLERR, "Incorrect args for bond coefficients");
}
/* ----------------------------------------------------------------------
@ -202,12 +202,12 @@ double BondGaussian::equilibrium_distance(int i)
void BondGaussian::write_restart(FILE *fp)
{
fwrite(&bond_temperature[1],sizeof(double),atom->nbondtypes,fp);
fwrite(&nterms[1],sizeof(int),atom->nbondtypes,fp);
fwrite(&bond_temperature[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&nterms[1], sizeof(int), atom->nbondtypes, fp);
for (int i = 1; i <= atom->nbondtypes; i++) {
fwrite(alpha[i],sizeof(double),nterms[i],fp);
fwrite(width[i],sizeof(double),nterms[i],fp);
fwrite(r0[i],sizeof(double),nterms[i],fp);
fwrite(alpha[i], sizeof(double), nterms[i], fp);
fwrite(width[i], sizeof(double), nterms[i], fp);
fwrite(r0[i], sizeof(double), nterms[i], fp);
}
}
@ -220,31 +220,32 @@ void BondGaussian::read_restart(FILE *fp)
allocate();
if (comm->me == 0) {
utils::sfread(FLERR,&bond_temperature[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
utils::sfread(FLERR,&nterms[1],sizeof(int),atom->nbondtypes,fp,nullptr,error);
utils::sfread(FLERR, &bond_temperature[1], sizeof(double), atom->nbondtypes, fp, nullptr,
error);
utils::sfread(FLERR, &nterms[1], sizeof(int), atom->nbondtypes, fp, nullptr, error);
}
MPI_Bcast(&bond_temperature[1],atom->nbondtypes,MPI_DOUBLE,0,world);
MPI_Bcast(&nterms[1],atom->nbondtypes,MPI_INT,0,world);
MPI_Bcast(&bond_temperature[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&nterms[1], atom->nbondtypes, MPI_INT, 0, world);
// allocate
for (int i = 1; i <= atom->nbondtypes; i++) {
alpha[i] = new double [nterms[i]];
width[i] = new double [nterms[i]];
r0[i] = new double [nterms[i]];
alpha[i] = new double[nterms[i]];
width[i] = new double[nterms[i]];
r0[i] = new double[nterms[i]];
}
if (comm->me == 0) {
for (int i = 1; i <= atom->nbondtypes; i++) {
utils::sfread(FLERR,alpha[i],sizeof(double),nterms[i],fp,nullptr,error);
utils::sfread(FLERR,width[i],sizeof(double),nterms[i],fp,nullptr,error);
utils::sfread(FLERR,r0[i],sizeof(double),nterms[i],fp,nullptr,error);
utils::sfread(FLERR, alpha[i], sizeof(double), nterms[i], fp, nullptr, error);
utils::sfread(FLERR, width[i], sizeof(double), nterms[i], fp, nullptr, error);
utils::sfread(FLERR, r0[i], sizeof(double), nterms[i], fp, nullptr, error);
}
}
for (int i = 1; i <= atom->nbondtypes; i++) {
MPI_Bcast(alpha[i],nterms[i],MPI_DOUBLE,0,world);
MPI_Bcast(width[i],nterms[i],MPI_DOUBLE,0,world);
MPI_Bcast(r0[i],nterms[i],MPI_DOUBLE,0,world);
MPI_Bcast(alpha[i], nterms[i], MPI_DOUBLE, 0, world);
MPI_Bcast(width[i], nterms[i], MPI_DOUBLE, 0, world);
MPI_Bcast(r0[i], nterms[i], MPI_DOUBLE, 0, world);
}
for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1;
@ -257,9 +258,9 @@ void BondGaussian::read_restart(FILE *fp)
void BondGaussian::write_data(FILE *fp)
{
for (int i = 1; i <= atom->nbondtypes; i++) {
fprintf(fp,"%d %g %d",i,bond_temperature[i],nterms[i]);
fprintf(fp, "%d %g %d", i, bond_temperature[i], nterms[i]);
for (int j = 0; j < nterms[i]; j++) {
fprintf(fp," %g %g %g",alpha[i][j],width[i][j],r0[i][j]);
fprintf(fp, " %g %g %g", alpha[i][j], width[i][j], r0[i][j]);
}
fprintf(fp, "\n");
}
@ -267,26 +268,25 @@ void BondGaussian::write_data(FILE *fp)
/* ---------------------------------------------------------------------- */
double BondGaussian::single(int type, double rsq, int /*i*/, int /*j*/,
double &fforce)
double BondGaussian::single(int type, double rsq, int /*i*/, int /*j*/, double &fforce)
{
double r = sqrt(rsq);
fforce = 0;
double sum_g_i = 0.0;
double sum_numerator = 0.0;
for (int i = 0; i < nterms[type]; i++) {
double dr = r - r0[type][i];
double prefactor = (alpha[type][i]/(width[type][i]*sqrt(MY_PI2)));
double exponent = -2*dr*dr/(width[type][i]*width[type][i]);
double g_i = prefactor*exp(exponent);
sum_g_i += g_i;
sum_numerator += g_i*dr/(width[type][i]*width[type][i]);
}
for (int i = 0; i < nterms[type]; i++) {
double dr = r - r0[type][i];
double prefactor = (alpha[type][i] / (width[type][i] * sqrt(MY_PI2)));
double exponent = -2 * dr * dr / (width[type][i] * width[type][i]);
double g_i = prefactor * exp(exponent);
sum_g_i += g_i;
sum_numerator += g_i * dr / (width[type][i] * width[type][i]);
}
if (sum_g_i < SMALL) sum_g_i = SMALL;
if (r > 0.0) fforce = -4.0*(force->boltz*bond_temperature[type])*(sum_numerator/sum_g_i)/r;
if (r > 0.0)
fforce = -4.0 * (force->boltz * bond_temperature[type]) * (sum_numerator / sum_g_i) / r;
return -(force->boltz*bond_temperature[type])*log(sum_g_i);
return -(force->boltz * bond_temperature[type]) * log(sum_g_i);
}

2
src/EXTRA-PAIR/pair_coul_exclude.cpp
View File

@ -189,7 +189,7 @@ void PairCoulExclude::init_style()
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairCoulExclude::init_one(int i, int j)
double PairCoulExclude::init_one(int /*i*/, int /*j*/)
{
return cut_global;
}

14
src/FEP/compute_fep.cpp
View File

@ -124,18 +124,12 @@ ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) :
while (iarg < narg) {
if (strcmp(arg[iarg],"tail") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal optional keyword "
"in compute fep");
if (strcmp(arg[iarg+1],"no") == 0) tailflag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) tailflag = 1;
else error->all(FLERR,"Illegal optional keyword in compute fep");
if (iarg+2 > narg) error->all(FLERR,"Illegal optional keyword in compute fep");
tailflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"volume") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal optional keyword "
"in compute fep");
if (strcmp(arg[iarg+1],"no") == 0) volumeflag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) volumeflag = 1;
else error->all(FLERR,"Illegal optional keyword in compute fep");
if (iarg+2 > narg) error->all(FLERR,"Illegal optional keyword in compute fep");
volumeflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else
error->all(FLERR,"Illegal optional keyword in compute fep");

12
src/FEP/fix_adapt_fep.cpp
View File

@ -138,21 +138,15 @@ FixAdaptFEP::FixAdaptFEP(LAMMPS *lmp, int narg, char **arg) :
while (iarg < narg) {
if (strcmp(arg[iarg],"reset") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt/fep command");
if (strcmp(arg[iarg+1],"no") == 0) resetflag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) resetflag = 1;
else error->all(FLERR,"Illegal fix adapt/fep command");
resetflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"scale") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt/fep command");
if (strcmp(arg[iarg+1],"no") == 0) scaleflag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) scaleflag = 1;
else error->all(FLERR,"Illegal fix adapt/fep command");
scaleflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"after") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt/fep command");
if (strcmp(arg[iarg+1],"no") == 0) afterflag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) afterflag = 1;
else error->all(FLERR,"Illegal fix adapt/fep command");
afterflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else error->all(FLERR,"Illegal fix adapt/fep command");
}

17
src/GPU/fix_gpu.cpp
View File

@ -136,16 +136,17 @@ FixGPU::FixGPU(LAMMPS *lmp, int narg, char **arg) :
while (iarg < narg) {
if (strcmp(arg[iarg],"neigh") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal package gpu command");
if (strcmp(arg[iarg+1],"yes") == 0) _gpu_mode = GPU_NEIGH;
else if (strcmp(arg[iarg+1],"no") == 0) _gpu_mode = GPU_FORCE;
else if (strcmp(arg[iarg+1],"hybrid") == 0) _gpu_mode = GPU_HYB_NEIGH;
const std::string modearg = arg[iarg+1];
if ((modearg == "yes") || (modearg == "on") || (modearg == "true"))
_gpu_mode = GPU_NEIGH;
else if ((modearg == "no") || (modearg == "off") || (modearg == "false"))
_gpu_mode = GPU_FORCE;
else if (modearg == "hybrid") _gpu_mode = GPU_HYB_NEIGH;
else error->all(FLERR,"Illegal package gpu command");
iarg += 2;
} else if (strcmp(arg[iarg],"newton") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal package gpu command");
if (strcmp(arg[iarg+1],"off") == 0) newtonflag = 0;
else if (strcmp(arg[iarg+1],"on") == 0) newtonflag = 1;
else error->all(FLERR,"Illegal package gpu command");
newtonflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"binsize") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal package gpu command");
@ -185,9 +186,7 @@ FixGPU::FixGPU(LAMMPS *lmp, int narg, char **arg) :
iarg += 2;
} else if (strcmp(arg[iarg],"pair/only") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal package gpu command");
if (strcmp(arg[iarg+1],"off") == 0) pair_only_flag = 0;
else if (strcmp(arg[iarg+1],"on") == 0) pair_only_flag = 1;
else error->all(FLERR,"Illegal package gpu command");
pair_only_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"ocl_args") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal package gpu command");

6
src/GPU/gpu_extra.h
View File

@ -78,7 +78,11 @@ inline void check_flag(int error_flag, LAMMPS_NS::Error *error, MPI_Comm &world)
else if (all_success == -13)
error->all(FLERR, "Invalid device configuration.");
else if (all_success == -15)
error->all(FLERR, "P3M built for FP64 and GPU device is FP32 only.");
error->all(FLERR, "PPPM was compiled for double precision floating point "
"but GPU device supports single precision only.");
else if (all_success == -16)
error->all(FLERR, "GPU library was compiled for double or mixed precision "
"floating point but GPU device supports single precision only.");
else
error->all(FLERR, "Unknown error in GPU library");
}

40
src/H5MD/dump_h5md.cpp
View File

@ -16,23 +16,24 @@
Contributing author: Pierre de Buyl (KU Leuven)
------------------------------------------------------------------------- */
#include <cmath>
#include <cstdio>
#include <cstring>
#include <climits>
#include "ch5md.h"
#include "dump_h5md.h"
#include "domain.h"
#include "atom.h"
#include "update.h"
#include "group.h"
#include "output.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "memory.h"
#include "output.h"
#include "update.h"
#include "version.h"
#include "ch5md.h"
#include <climits>
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
#define MYMIN(a,b) ((a) < (b) ? (a) : (b))
@ -152,25 +153,14 @@ DumpH5MD::DumpH5MD(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg)
error->all(FLERR, "Invalid number of arguments in dump h5md");
}
box_is_set = true;
if (strcmp(arg[iarg+1], "yes")==0)
do_box=true;
else if (strcmp(arg[iarg+1], "no")==0)
do_box=false;
else
error->all(FLERR, "Illegal dump h5md command");
do_box = utils::logical(FLERR,arg[iarg+1],false,lmp) == 1;
iarg+=2;
} else if (strcmp(arg[iarg], "create_group")==0) {
if (iarg+1>=narg) {
error->all(FLERR, "Invalid number of arguments in dump h5md");
}
create_group_is_set = true;
if (strcmp(arg[iarg+1], "yes")==0)
create_group=true;
else if (strcmp(arg[iarg+1], "no")==0) {
create_group=false;
}
else
error->all(FLERR, "Illegal dump h5md command");
create_group = utils::logical(FLERR,arg[iarg+1],false,lmp) == 1;
iarg+=2;
} else if (strcmp(arg[iarg], "author")==0) {
if (iarg+1>=narg) {
@ -470,9 +460,7 @@ int DumpH5MD::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"unwrap") == 0) {
if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
if (strcmp(arg[1],"yes") == 0) unwrap_flag = 1;
else if (strcmp(arg[1],"no") == 0) unwrap_flag = 0;
else error->all(FLERR,"Illegal dump_modify command");
unwrap_flag = utils::logical(FLERR, arg[1], false, lmp);
return 2;
}
return 0;

8
src/INTEL/fix_intel.cpp
View File

@ -118,9 +118,7 @@ FixIntel::FixIntel(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
iarg += 2;
} else if (strcmp(arg[iarg], "ghost") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal package intel command");
if (strcmp(arg[iarg+1],"yes") == 0) _offload_ghost = 1;
else if (strcmp(arg[iarg+1],"no") == 0) _offload_ghost = 0;
else error->all(FLERR,"Illegal package intel command");
_offload_ghost = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "tpc") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal package intel command");
@ -135,9 +133,7 @@ FixIntel::FixIntel(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
iarg++;
} else if (strcmp(arg[iarg], "lrt") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal package intel command");
if (strcmp(arg[iarg+1],"yes") == 0) _lrt = 1;
else if (strcmp(arg[iarg+1],"no") == 0) _lrt = 0;
else error->all(FLERR,"Illegal package intel command");
_lrt = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
}

19
src/INTERLAYER/pair_drip.cpp
View File

@ -31,7 +31,6 @@
#include "neigh_request.h"
#include "neighbor.h"
#include "potential_file_reader.h"
#include "tokenizer.h"
#include <cmath>
#include <cstring>
@ -241,17 +240,17 @@ void PairDRIP::read_file(char *filename)
nparams++;
}
MPI_Bcast(&nparams, 1, MPI_INT, 0, world);
MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);
if (comm->me != 0) {
params = (Param *) memory->srealloc(params, maxparam * sizeof(Param), "pair:params");
}
MPI_Bcast(params, maxparam * sizeof(Param), MPI_BYTE, 0, world);
}
MPI_Bcast(&nparams, 1, MPI_INT, 0, world);
MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);
if (comm->me != 0) {
params = (Param *) memory->srealloc(params, maxparam * sizeof(Param), "pair:params");
}
MPI_Bcast(params, maxparam * sizeof(Param), MPI_BYTE, 0, world);
memory->destroy(elem2param);
memory->create(elem2param, nelements, nelements, "pair:elem2param");
for (int i = 0; i < nelements; i++) {

279
src/KIM/kim_init.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
@ -85,21 +84,23 @@ void KimInit::command(int narg, char **arg)
if ((narg < 2) || (narg > 3)) error->all(FLERR, "Illegal 'kim init' command");
if (domain->box_exist)
error->all(FLERR, "Must use 'kim init' command before "
"simulation box is defined");
error->all(FLERR, "Must use 'kim init' command before simulation box is defined");
char *model_name = utils::strdup(arg[0]);
char *user_units = utils::strdup(arg[1]);
if (narg == 3) {
auto arg_str = std::string(arg[2]);
if (arg_str == "unit_conversion_mode") unit_conversion_mode = true;
if (arg_str == "unit_conversion_mode")
unit_conversion_mode = true;
else {
error->all(FLERR, "Illegal 'kim init' command.\nThe argument "
"followed by unit_style {} is an optional "
"argument and when is used must "
"be unit_conversion_mode", user_units);
error->all(FLERR,
"Illegal 'kim init' command.\n"
"The argument followed by unit_style {} is an optional argument and when "
"is used must be unit_conversion_mode",
user_units);
}
} else unit_conversion_mode = false;
} else
unit_conversion_mode = false;
char *model_units;
KIM_Model *pkim = nullptr;
@ -117,14 +118,9 @@ void KimInit::command(int narg, char **arg)
/* ---------------------------------------------------------------------- */
namespace {
void get_kim_unit_names(
char const * const system,
KIM_LengthUnit & lengthUnit,
KIM_EnergyUnit & energyUnit,
KIM_ChargeUnit & chargeUnit,
KIM_TemperatureUnit & temperatureUnit,
KIM_TimeUnit & timeUnit,
Error * error)
void get_kim_unit_names(char const *const system, KIM_LengthUnit &lengthUnit,
KIM_EnergyUnit &energyUnit, KIM_ChargeUnit &chargeUnit,
KIM_TemperatureUnit &temperatureUnit, KIM_TimeUnit &timeUnit, Error *error)
{
const std::string system_str(system);
if (system_str == "real") {
@ -157,20 +153,64 @@ void get_kim_unit_names(
chargeUnit = KIM_CHARGE_UNIT_e;
temperatureUnit = KIM_TEMPERATURE_UNIT_K;
timeUnit = KIM_TIME_UNIT_fs;
} else if ((system_str == "lj") ||
(system_str == "micro") ||
(system_str == "nano")) {
error->all(FLERR, "LAMMPS unit_style {} not supported "
"by KIM models", system_str);
} else if ((system_str == "lj") || (system_str == "micro") || (system_str == "nano")) {
error->all(FLERR, "LAMMPS unit_style {} not supported by KIM models", system_str);
} else {
error->all(FLERR, "Unknown unit_style");
}
}
} // namespace
} // namespace
void KimInit::determine_model_type_and_units(char * model_name,
char * user_units,
char ** model_units,
void KimInit::print_dirs(struct KIM_Collections *const collections) const
{
int kim_error = 0;
int dirListExtent = 0;
int dirCounter = 0;
std::string mesg = "#=== KIM is looking for 'Portable Models' in these directories ===\n";
std::vector<struct KIM_Collection> collection_list;
collection_list.push_back(KIM_COLLECTION_currentWorkingDirectory);
collection_list.push_back(KIM_COLLECTION_environmentVariable);
collection_list.push_back(KIM_COLLECTION_user);
collection_list.push_back(KIM_COLLECTION_system);
for (auto col : collection_list) {
kim_error = KIM_Collections_CacheListOfDirectoryNames(
collections, col, KIM_COLLECTION_ITEM_TYPE_portableModel, &dirListExtent);
if (!kim_error) {
for (int i = 0; i < dirListExtent; ++i) {
char const *name;
kim_error = KIM_Collections_GetDirectoryName(collections, i, &name);
// Don't check for error due to bug in kim-api-2.2.1 and below.
#if ((KIM_VERSION_MAJOR * 1000 + KIM_VERSION_MINOR) * 1000 + KIM_VERSION_PATCH) <= 2002001
kim_error = 0;
#endif
if (!kim_error) mesg += fmt::format("# {:2}: {}\n", ++dirCounter, name);
}
}
}
dirCounter = 0;
mesg += "#=== KIM is looking for 'Simulator Models' in these directories ===\n";
for (auto col : collection_list) {
kim_error = KIM_Collections_CacheListOfDirectoryNames(
collections, col, KIM_COLLECTION_ITEM_TYPE_simulatorModel, &dirListExtent);
if (!kim_error) {
for (int i = 0; i < dirListExtent; ++i) {
char const *name;
kim_error = KIM_Collections_GetDirectoryName(collections, i, &name);
// Don't check for error due to bug in kim-api-2.2.1 and below.
#if ((KIM_VERSION_MAJOR * 1000 + KIM_VERSION_MINOR) * 1000 + KIM_VERSION_PATCH) <= 2002001
kim_error = 0;
#endif
if (!kim_error) mesg += fmt::format("# {:2}: {}\n", ++dirCounter, name);
}
}
}
input->write_echo(mesg);
}
void KimInit::determine_model_type_and_units(char *model_name, char *user_units, char **model_units,
KIM_Model *&pkim)
{
KIM_LengthUnit lengthUnit;
@ -179,33 +219,26 @@ void KimInit::determine_model_type_and_units(char * model_name,
KIM_TemperatureUnit temperatureUnit;
KIM_TimeUnit timeUnit;
int units_accepted;
KIM_Collections * collections;
KIM_Collections *collections;
KIM_CollectionItemType itemType;
int kim_error = KIM_Collections_Create(&collections);
if (kim_error)
error->all(FLERR, "Unable to access KIM Collections to find Model");
if (kim_error) error->all(FLERR, "Unable to access KIM Collections to find Model");
auto logID = fmt::format("{}_Collections", comm->me);
KIM_Collections_SetLogID(collections, logID.c_str());
print_dirs(collections);
kim_error = KIM_Collections_GetItemType(collections, model_name, &itemType);
if (kim_error) error->all(FLERR, "KIM Model name not found");
KIM_Collections_Destroy(&collections);
if (KIM_CollectionItemType_Equal(itemType,
KIM_COLLECTION_ITEM_TYPE_portableModel)) {
get_kim_unit_names(user_units, lengthUnit, energyUnit,
chargeUnit, temperatureUnit, timeUnit, error);
int kim_error = KIM_Model_Create(KIM_NUMBERING_zeroBased,
lengthUnit,
energyUnit,
chargeUnit,
temperatureUnit,
timeUnit,
model_name,
&units_accepted,
&pkim);
if (KIM_CollectionItemType_Equal(itemType, KIM_COLLECTION_ITEM_TYPE_portableModel)) {
get_kim_unit_names(user_units, lengthUnit, energyUnit, chargeUnit, temperatureUnit, timeUnit,
error);
int kim_error = KIM_Model_Create(KIM_NUMBERING_zeroBased, lengthUnit, energyUnit, chargeUnit,
temperatureUnit, timeUnit, model_name, &units_accepted, &pkim);
if (kim_error) error->all(FLERR, "Unable to load KIM Simulator Model");
@ -219,20 +252,12 @@ void KimInit::determine_model_type_and_units(char * model_name,
} else if (unit_conversion_mode) {
KIM_Model_Destroy(&pkim);
int const num_systems = 5;
char const * const systems[num_systems]
= {"metal", "real", "si", "cgs", "electron"};
for (int i=0; i < num_systems; ++i) {
get_kim_unit_names(systems[i], lengthUnit, energyUnit,
chargeUnit, temperatureUnit, timeUnit, error);
kim_error = KIM_Model_Create(KIM_NUMBERING_zeroBased,
lengthUnit,
energyUnit,
chargeUnit,
temperatureUnit,
timeUnit,
model_name,
&units_accepted,
&pkim);
char const *const systems[num_systems] = {"metal", "real", "si", "cgs", "electron"};
for (int i = 0; i < num_systems; ++i) {
get_kim_unit_names(systems[i], lengthUnit, energyUnit, chargeUnit, temperatureUnit,
timeUnit, error);
kim_error = KIM_Model_Create(KIM_NUMBERING_zeroBased, lengthUnit, energyUnit, chargeUnit,
temperatureUnit, timeUnit, model_name, &units_accepted, &pkim);
if (units_accepted) {
logID = fmt::format("{}_Model", comm->me);
KIM_Model_SetLogID(pkim, logID.c_str());
@ -246,12 +271,10 @@ void KimInit::determine_model_type_and_units(char * model_name,
KIM_Model_Destroy(&pkim);
error->all(FLERR, "KIM Model does not support the requested unit system");
}
} else if (KIM_CollectionItemType_Equal(
itemType, KIM_COLLECTION_ITEM_TYPE_simulatorModel)) {
KIM_SimulatorModel * simulatorModel;
} else if (KIM_CollectionItemType_Equal(itemType, KIM_COLLECTION_ITEM_TYPE_simulatorModel)) {
KIM_SimulatorModel *simulatorModel;
kim_error = KIM_SimulatorModel_Create(model_name, &simulatorModel);
if (kim_error)
error->all(FLERR, "Unable to load KIM Simulator Model");
if (kim_error) error->all(FLERR, "Unable to load KIM Simulator Model");
model_type = SM;
logID = fmt::format("{}_SimulatorModel", comm->me);
@ -264,13 +287,11 @@ void KimInit::determine_model_type_and_units(char * model_name,
KIM_SimulatorModel_GetNumberOfSimulatorFields(simulatorModel, &sim_fields);
KIM_SimulatorModel_CloseTemplateMap(simulatorModel);
for (int i = 0; i < sim_fields; ++i) {
KIM_SimulatorModel_GetSimulatorFieldMetadata(
simulatorModel, i, &sim_lines, &sim_field);
KIM_SimulatorModel_GetSimulatorFieldMetadata(simulatorModel, i, &sim_lines, &sim_field);
const std::string sim_field_str(sim_field);
if (sim_field_str == "units") {
KIM_SimulatorModel_GetSimulatorFieldLine(
simulatorModel, i, 0, &sim_value);
KIM_SimulatorModel_GetSimulatorFieldLine(simulatorModel, i, 0, &sim_value);
*model_units = utils::strdup(sim_value);
break;
}
@ -280,16 +301,15 @@ void KimInit::determine_model_type_and_units(char * model_name,
const std::string model_units_str(*model_units);
const std::string user_units_str(user_units);
if ((!unit_conversion_mode) && (model_units_str != user_units_str)) {
error->all(FLERR, "Incompatible units for KIM Simulator Model"
", required units = {}", model_units_str);
error->all(FLERR, "Incompatible units for KIM Simulator Model, required units = {}",
model_units_str);
}
}
}
/* ---------------------------------------------------------------------- */
void KimInit::do_init(char *model_name, char *user_units, char *model_units,
KIM_Model *&pkim)
void KimInit::do_init(char *model_name, char *user_units, char *model_units, KIM_Model *&pkim)
{
// create storage proxy fix. delete existing fix, if needed.
@ -304,8 +324,7 @@ void KimInit::do_init(char *model_name, char *user_units, char *model_units,
fix_store->setptr("model_units", (void *) model_units);
// Begin output to log file
input->write_echo("#=== BEGIN kim init ==================================="
"=======\n");
input->write_echo("#=== BEGIN kim init ==========================================\n");
KIM_SimulatorModel *simulatorModel;
if (model_type == SM) {
@ -316,18 +335,16 @@ void KimInit::do_init(char *model_name, char *user_units, char *model_units,
KIM_SimulatorModel_SetLogID(simulatorModel, logID.c_str());
char const *sim_name, *sim_version;
KIM_SimulatorModel_GetSimulatorNameAndVersion(
simulatorModel, &sim_name, &sim_version);
KIM_SimulatorModel_GetSimulatorNameAndVersion(simulatorModel, &sim_name, &sim_version);
const std::string sim_name_str(sim_name);
if (sim_name_str != "LAMMPS")
error->all(FLERR, "Incompatible KIM Simulator Model");
if (sim_name_str != "LAMMPS") error->all(FLERR, "Incompatible KIM Simulator Model");
if (comm->me == 0) {
auto mesg = fmt::format("# Using KIM Simulator Model : {}\n"
"# For Simulator : {} {}\n"
"# Running on : LAMMPS {}\n#\n", model_name,
sim_name_str, sim_version, lmp->version);
"# For Simulator : {} {}\n"
"# Running on : LAMMPS {}\n#\n",
model_name, sim_name_str, sim_version, lmp->version);
utils::logmesg(lmp, mesg);
}
@ -350,18 +367,16 @@ void KimInit::do_init(char *model_name, char *user_units, char *model_units,
// Set the skin and timestep default values as
// 2.0 Angstroms and 1.0 femtosecond
const std::string skin_cmd =
(model_units_str == "real") ? "neighbor 2.0 bin # Angstroms":
(model_units_str == "metal") ? "neighbor 2.0 bin # Angstroms":
(model_units_str == "si") ? "neighbor 2e-10 bin # meters":
(model_units_str == "cgs") ? "neighbor 2e-8 bin # centimeters":
"neighbor 3.77945224 bin # Bohr";
const std::string step_cmd =
(model_units_str == "real") ? "timestep 1.0 # femtoseconds":
(model_units_str == "metal") ? "timestep 1.0e-3 # picoseconds":
(model_units_str == "si") ? "timestep 1e-15 # seconds":
(model_units_str == "cgs") ? "timestep 1e-15 # seconds":
"timestep 1.0 # femtoseconds";
const std::string skin_cmd = (model_units_str == "real") ? "neighbor 2.0 bin # Angstroms"
: (model_units_str == "metal") ? "neighbor 2.0 bin # Angstroms"
: (model_units_str == "si") ? "neighbor 2e-10 bin # meters"
: (model_units_str == "cgs") ? "neighbor 2e-8 bin # centimeters"
: "neighbor 3.77945224 bin # Bohr";
const std::string step_cmd = (model_units_str == "real") ? "timestep 1.0 # femtoseconds"
: (model_units_str == "metal") ? "timestep 1.0e-3 # picoseconds"
: (model_units_str == "si") ? "timestep 1e-15 # seconds"
: (model_units_str == "cgs") ? "timestep 1e-15 # seconds"
: "timestep 1.0 # femtoseconds";
input->one(skin_cmd);
input->one(step_cmd);
@ -373,14 +388,12 @@ void KimInit::do_init(char *model_name, char *user_units, char *model_units,
// init model
for (int i = 0; i < sim_fields; ++i) {
KIM_SimulatorModel_GetSimulatorFieldMetadata(
simulatorModel, i, &sim_lines, &sim_field);
KIM_SimulatorModel_GetSimulatorFieldMetadata(simulatorModel, i, &sim_lines, &sim_field);
const std::string sim_field_str(sim_field);
if (sim_field_str == "model-init") {
for (int j = 0; j < sim_lines; ++j) {
KIM_SimulatorModel_GetSimulatorFieldLine(
simulatorModel, i, j, &sim_value);
KIM_SimulatorModel_GetSimulatorFieldLine(simulatorModel, i, j, &sim_value);
input->one(sim_value);
}
break;
@ -404,31 +417,28 @@ void KimInit::do_init(char *model_name, char *user_units, char *model_units,
int max_len(0);
for (int i = 0; i < numberOfParameters; ++i) {
KIM_Model_GetParameterMetadata(pkim, i, &kim_DataType,
&extent, &str_name, &str_desc);
max_len = MAX(max_len, (int)strlen(str_name));
KIM_Model_GetParameterMetadata(pkim, i, &kim_DataType, &extent, &str_name, &str_desc);
max_len = MAX(max_len, (int) strlen(str_name));
}
max_len = MAX(18, max_len + 1);
mesg += fmt::format(" No. | {:<{}} | data type | extent\n",
"Parameter name", max_len);
mesg += fmt::format(" No. | {:<{}} | data type | extent\n", "Parameter name", max_len);
mesg += fmt::format("{:-<{}}\n", "-", max_len + 35);
for (int i = 0; i < numberOfParameters; ++i) {
KIM_Model_GetParameterMetadata(pkim, i, &kim_DataType,
&extent, &str_name, &str_desc);
KIM_Model_GetParameterMetadata(pkim, i, &kim_DataType, &extent, &str_name, &str_desc);
auto data_type = std::string("\"");
data_type += KIM_DataType_ToString(kim_DataType) + std::string("\"");
mesg += fmt::format(" {:<8} | {:<{}} | {:<10} | {}\n", i + 1, str_name,
max_len, data_type, extent);
mesg += fmt::format(" {:<8} | {:<{}} | {:<10} | {}\n", i + 1, str_name, max_len, data_type,
extent);
}
} else mesg += "No mutable parameters.\n";
} else
mesg += "No mutable parameters.\n";
KIM_Model_Destroy(&pkim);
input->write_echo(mesg);
}
// End output to log file
input->write_echo("#=== END kim init ====================================="
"=======\n\n");
input->write_echo("#=== END kim init ============================================\n\n");
}
/* ---------------------------------------------------------------------- */
@ -446,24 +456,11 @@ void KimInit::do_variables(const std::string &from, const std::string &to)
int ier;
std::string var_str;
int v_unit;
const char *units[] = {"mass",
"distance",
"time",
"energy",
"velocity",
"force",
"torque",
"temperature",
"pressure",
"viscosity",
"charge",
"dipole",
"efield",
"density",
nullptr};
const char *units[] = {"mass", "distance", "time", "energy", "velocity",
"force", "torque", "temperature", "pressure", "viscosity",
"charge", "dipole", "efield", "density", nullptr};
input->write_echo(fmt::format("# Conversion factors from {} to {}:\n",
from, to));
input->write_echo(fmt::format("# Conversion factors from {} to {}:\n", from, to));
auto variable = input->variable;
for (int i = 0; units[i] != nullptr; ++i) {
@ -473,24 +470,23 @@ void KimInit::do_variables(const std::string &from, const std::string &to)
variable->set(var_str + " internal 1.0");
v_unit = variable->find(var_str.c_str());
}
ier = lammps_unit_conversion(units[i], from, to,
conversion_factor);
ier = lammps_unit_conversion(units[i], from, to, conversion_factor);
if (ier != 0)
error->all(FLERR, "Unable to obtain conversion factor: "
"unit = {}; from = {}; to = {}",
units[i], from, to);
error->all(FLERR,
"Unable to obtain conversion factor: "
"unit = {}; from = {}; to = {}",
units[i], from, to);
variable->internal_set(v_unit, conversion_factor);
input->write_echo(fmt::format("variable {:<15s} internal {:<15.12e}\n",
var_str, conversion_factor));
input->write_echo(
fmt::format("variable {:<15s} internal {:<15.12e}\n", var_str, conversion_factor));
}
input->write_echo("#\n");
}
/* ---------------------------------------------------------------------- */
void KimInit::write_log_cite(class LAMMPS *lmp,
KimInit::model_type_enum model_type,
void KimInit::write_log_cite(class LAMMPS *lmp, KimInit::model_type_enum model_type,
char *model_name)
{
if (!lmp->citeme) return;
@ -501,7 +497,7 @@ void KimInit::write_log_cite(class LAMMPS *lmp,
std::string cite_id;
if (kim_id.empty()) {
cite_id = fmt::format("KIM potential: unpublished, \"{}\"\n",model_name_str);
cite_id = fmt::format("KIM potential: unpublished, \"{}\"\n", model_name_str);
} else {
KIM_Collections *collections;
int err = KIM_Collections_Create(&collections);
@ -513,12 +509,10 @@ void KimInit::write_log_cite(class LAMMPS *lmp,
int extent;
if (model_type == MO) {
err = KIM_Collections_CacheListOfItemMetadataFiles(
collections, KIM_COLLECTION_ITEM_TYPE_portableModel,
model_name, &extent);
collections, KIM_COLLECTION_ITEM_TYPE_portableModel, model_name, &extent);
} else if (model_type == SM) {
err = KIM_Collections_CacheListOfItemMetadataFiles(
collections, KIM_COLLECTION_ITEM_TYPE_simulatorModel,
model_name, &extent);
collections, KIM_COLLECTION_ITEM_TYPE_simulatorModel, model_name, &extent);
} else {
lmp->error->all(FLERR, "Unknown model type");
}
@ -529,19 +523,18 @@ void KimInit::write_log_cite(class LAMMPS *lmp,
}
cite_id = fmt::format("OpenKIM potential: https://openkim.org/cite/"
"{}#item-citation\n\n",kim_id);
"{}#item-citation\n\n",
kim_id);
for (int i = 0; i < extent; ++i) {
char const *fileName;
int availableAsString;
char const *fileString;
err = KIM_Collections_GetItemMetadataFile(
collections, i, &fileName, nullptr, nullptr,
&availableAsString, &fileString);
err = KIM_Collections_GetItemMetadataFile(collections, i, &fileName, nullptr, nullptr,
&availableAsString, &fileString);
if (err) continue;
if (utils::strmatch(fileName, "^kimcite") && availableAsString)
cite_id += fileString;
if (utils::strmatch(fileName, "^kimcite") && availableAsString) cite_id += fileString;
}
KIM_Collections_Destroy(&collections);
}

5
src/KIM/kim_init.h
View File

@ -62,7 +62,8 @@
#include "pointers.h"
// Forward declaration.
typedef struct KIM_Model KIM_Model;
struct KIM_Model;
struct KIM_Collections;
namespace LAMMPS_NS {
@ -80,6 +81,8 @@ class KimInit : protected Pointers {
void determine_model_type_and_units(char *, char *, char **, KIM_Model *&);
void do_init(char *, char *, char *, KIM_Model *&);
void do_variables(const std::string &, const std::string &);
void print_dirs(struct KIM_Collections * const collections) const;
};
} // namespace LAMMPS_NS

2
src/KOKKOS/compute_coord_atom_kokkos.cpp
View File

@ -59,7 +59,7 @@ ComputeCoordAtomKokkos<DeviceType>::ComputeCoordAtomKokkos(LAMMPS *lmp, int narg
/* ---------------------------------------------------------------------- */
template<class DeviceType>
ComputeCoordAtomKokkos<DeviceType>::~ComputeCoordAtomKokkos<DeviceType>()
ComputeCoordAtomKokkos<DeviceType>::~ComputeCoordAtomKokkos()
{
if (copymode) return;

4
src/KOKKOS/fix_rx_kokkos.h
View File

@ -196,7 +196,7 @@ class FixRxKokkos : public FixRX {
double& h0, VectorType& y, VectorType& rwk, UserDataType& userData) const;
//!< ODE Solver diagnostics.
void odeDiagnostics(void);
void odeDiagnostics();
//!< Special counters per-ode.
int *diagnosticCounterPerODEnSteps;
@ -231,7 +231,7 @@ class FixRxKokkos : public FixRX {
bool update_kinetics_data;
void create_kinetics_data(void);
void create_kinetics_data();
// Need a dual-view and device-view for dpdThetaLocal and sumWeights since they're used in several callbacks.
DAT::tdual_efloat_1d k_dpdThetaLocal, k_sumWeights;

16
src/KOKKOS/kokkos.cpp
View File

@ -381,9 +381,7 @@ void KokkosLMP::accelerator(int narg, char **arg)
iarg += 2;
} else if (strcmp(arg[iarg],"newton") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal package kokkos command");
if (strcmp(arg[iarg+1],"off") == 0) newtonflag = 0;
else if (strcmp(arg[iarg+1],"on") == 0) newtonflag = 1;
else error->all(FLERR,"Illegal package kokkos command");
newtonflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"comm") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal package kokkos command");
@ -459,21 +457,15 @@ void KokkosLMP::accelerator(int narg, char **arg)
} else if ((strcmp(arg[iarg],"gpu/aware") == 0)
|| (strcmp(arg[iarg],"cuda/aware") == 0)) {
if (iarg+2 > narg) error->all(FLERR,"Illegal package kokkos command");
if (strcmp(arg[iarg+1],"off") == 0) gpu_aware_flag = 0;
else if (strcmp(arg[iarg+1],"on") == 0) gpu_aware_flag = 1;
else error->all(FLERR,"Illegal package kokkos command");
gpu_aware_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"pair/only") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal package kokkos command");
if (strcmp(arg[iarg+1],"off") == 0) pair_only_flag = 0;
else if (strcmp(arg[iarg+1],"on") == 0) pair_only_flag = 1;
else error->all(FLERR,"Illegal package kokkos command");
pair_only_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"neigh/thread") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal package kokkos command");
if (strcmp(arg[iarg+1],"off") == 0) neigh_thread = 0;
else if (strcmp(arg[iarg+1],"on") == 0) neigh_thread = 1;
else error->all(FLERR,"Illegal package kokkos command");
neigh_thread = utils::logical(FLERR,arg[iarg+1],false,lmp);
neigh_thread_set = 1;
iarg += 2;
} else error->all(FLERR,"Illegal package kokkos command");

Some files were not shown because too many files have changed in this diff Show More