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

Merge branch 'develop' into fix-set-command

Browse Source
This commit is contained in:
Axel Kohlmeyer
2025-05-29 05:21:46 -04:00
parent 1668bcffcf 2744647c75
commit 129a3a83e5
207 changed files with 12636 additions and 4187 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
doc/src/Build.rst
View File

@ -14,32 +14,10 @@ As an alternative, you can download a package with pre-built executables
or automated build trees, as described in the :doc:`Install <Install>`
section of the manual.
Prerequisites
-------------
Which software you need to compile and use LAMMPS strongly depends on
which :doc:`features and settings <Build_settings>` and which
:doc:`optional packages <Packages_list>` you are trying to include.
Common to all is that you need a C++ and C compiler, where the C++
compiler has to support at least the C++11 standard (note that some
compilers require command-line flag to activate C++11 support).
Furthermore, if you are building with CMake, you need at least CMake
version 3.20 and a compatible build tool (make or ninja-build); if you
are building the the legacy GNU make based build system you need GNU
make (other make variants are not going to work since the build system
uses features unique to GNU make) and a Unix-like build environment with
a Bourne shell, and shell tools like "sed", "grep", "touch", "test",
"tr", "cp", "mv", "rm", "ln", "diff" and so on. Parts of LAMMPS
interface with or use Python version 3.6 or later.
The LAMMPS developers aim to keep LAMMPS very portable and usable -
at least in parts - on most operating systems commonly used for
running MD simulations. Please see the :doc:`section on portablility
<Intro_portability>` for more details.
.. toctree::
:maxdepth: 1
Build_prerequisites
Build_cmake
Build_make
Build_link

22
doc/src/Build_prerequisites.rst Normal file
View File

@ -0,0 +1,22 @@
Prerequisites
-------------
Which software you need to compile and use LAMMPS strongly depends on
which :doc:`features and settings <Build_settings>` and which
:doc:`optional packages <Packages_list>` you are trying to include.
Common to all is that you need a C++ and C compiler, where the C++
compiler has to support at least the C++11 standard (note that some
compilers require command-line flag to activate C++11 support).
Furthermore, if you are building with CMake, you need at least CMake
version 3.20 and a compatible build tool (make or ninja-build); if you
are building the the legacy GNU make based build system you need GNU
make (other make variants are not going to work since the build system
uses features unique to GNU make) and a Unix-like build environment with
a Bourne shell, and shell tools like "sed", "grep", "touch", "test",
"tr", "cp", "mv", "rm", "ln", "diff" and so on. Parts of LAMMPS
interface with or use Python version 3.6 or later.
The LAMMPS developers aim to keep LAMMPS very portable and usable -
at least in parts - on most operating systems commonly used for
running MD simulations. Please see the :doc:`section on portablility
<Intro_portability>` for more details.

2
doc/src/Commands_fix.rst
View File

@ -77,6 +77,7 @@ OPT.
* :doc:`flow/gauss <fix_flow_gauss>`
* :doc:`freeze (k) <fix_freeze>`
* :doc:`gcmc <fix_gcmc>`
* :doc:`gjf <fix_gjf>`
* :doc:`gld <fix_gld>`
* :doc:`gle <fix_gle>`
* :doc:`gravity (ko) <fix_gravity>`
@ -111,6 +112,7 @@ OPT.
* :doc:`mvv/tdpd <fix_mvv_dpd>`
* :doc:`neb <fix_neb>`
* :doc:`neb/spin <fix_neb_spin>`
* :doc:`neighbor/swap <fix_neighbor_swap>`
* :doc:`nonaffine/displacement <fix_nonaffine_displacement>`
* :doc:`nph (ko) <fix_nh>`
* :doc:`nph/asphere (o) <fix_nph_asphere>`

17
doc/src/Commands_removed.rst
View File

@ -1,7 +1,7 @@
Removed commands and packages
=============================
.. contents:: \
.. contents::
------
@ -12,10 +12,21 @@ stop LAMMPS and print a suitable error message in most cases, when a
style/command is used that has been removed or will replace the command
with the direct alternative (if available) and print a warning.
GJF formulation in fix langevin
-------------------------------
.. deprecated:: TBD
The *gjf* keyword in fix langevin is deprecated and will be removed
soon. The GJF functionality has been moved to its own fix style
:doc:`fix gjf <fix_gjf>` and it is strongly recommended to use that
fix instead.
LAMMPS shell
------------
.. versionchanged:: 29Aug2024
.. deprecated:: 29Aug2024
The LAMMPS shell has been removed from the LAMMPS distribution. Users
are encouraged to use the :ref:`LAMMPS-GUI <lammps_gui>` tool instead.
@ -23,7 +34,7 @@ are encouraged to use the :ref:`LAMMPS-GUI <lammps_gui>` tool instead.
i-PI tool
---------
.. versionchanged:: 27Jun2024
.. deprecated:: 27Jun2024
The i-PI tool has been removed from the LAMMPS distribution. Instead,
instructions to install i-PI from PyPI via pip are provided.

18
doc/src/Errors_details.rst
View File

@ -159,13 +159,17 @@ angle, dihedral, or improper with just one atom in the actual
sub-domain. Typically, this cutoff is set to the largest cutoff from
the :doc:`pair style(s) <pair_style>` plus the :doc:`neighbor list skin
distance <neighbor>` and will typically be sufficient for all bonded
interactions. But if the pair style cutoff is small, this may not be
enough. LAMMPS will print a warning in this case using some heuristic
based on the equilibrium bond length, but that still may not be
sufficient for cases where the force constants are small and thus bonds
may be stretched very far. The communication cutoff can be adjusted
with :doc:`comm_modify cutoff \<value\> <comm_modify>`, but setting this
too large will waste CPU time and memory.
interactions. But if the pair style cutoff is small (e.g. with a
repulsive-only Lennard-Jones potential) this may not be enough. It is
even worse if there is no pair style defined (or the pair style is set
to "none"), since then there will be no ghost atoms created at all.
The communication cutoff can be set or adjusted with :doc:`comm_modify
cutoff \<value\> <comm_modify>`, but setting this too large will waste
CPU time and memory. LAMMPS will print warnings in these cases. For
bonds it uses some heuristic based on the equilibrium bond length, but
that still may not be sufficient for cases where the force constants are
small and thus bonds may be stretched very far.
.. _hint09:

1
doc/src/Howto.rst
View File

@ -66,6 +66,7 @@ Force fields howto
:name: force_howto
:maxdepth: 1
Howto_FFgeneral
Howto_bioFF
Howto_amoeba
Howto_tip3p

55
doc/src/Howto_FFgeneral.rst Normal file
View File

@ -0,0 +1,55 @@
Some general force field considerations
=======================================
A compact summary of the concepts, definitions, and properties of force
fields with explicit bonded interactions (like the ones discussed in
this HowTo) is given in :ref:`(Gissinger) <Typelabel2>`.
A force field has 2 parts: the formulas that define its potential
functions and the coefficients used for a particular system. To assign
parameters it is first required to assign atom types. Those are not
only based on the elements, but also on the chemical environment due to
the atoms bound to them. This often follows the chemical concept of
*functional groups*. Example: a carbon atom bound with a single bond to
a single OH-group (alcohol) would be a different atom type than a carbon
atom bound to a methyl CH3 group (aliphatic carbon). The atom types
usually then determine the non-bonded Lennard-Jones parameters and the
parameters for bonds, angles, dihedrals, and impropers. On top of that,
partial charges have to be applied. Those are usually independent of
the atom types and are determined either for groups of atoms called
residues with some fitting procedure based on quantum mechanical
calculations, or based on some increment system that add or subtract
increments from the partial charge of an atom based on the types of
the neighboring atoms.
Force fields differ in the strategies they employ to determine the
parameters and charge distribution in how generic or specific they are
which in turn has an impact on the accuracy (compare for example
CGenFF to CHARMM and GAFF to Amber). Because of the different
strategies, it is not a good idea to use a mix of parameters from
different force field *families* (like CHARMM, Amber, or GROMOS)
and that extends to the parameters for the solvent, especially
water. The publication describing the parameterization of a force
field will describe which water model to use. Changing the water
model usually leads to overall worse results (even if it may improve
on the water itself).
In addition, one has to consider that *families* of force fields like
CHARMM, Amber, OPLS, or GROMOS have evolved over time and thus provide
different *revisions* of the force field parameters. These often
corresponds to changes in the functional form or the parameterization
strategies. This may also result in changes required for simulation
settings like the preferred cutoff or how Coulomb interactions are
computed (cutoff, smoothed/shifted cutoff, or long-range with Ewald
summation or equivalent). Unless explicitly stated in the publication
describing the force field, the Coulomb interaction cannot be chosen at
will but must match the revision of the force field. That said,
liberties may be taken during the initial equilibration of a system to
speed up the process, but not for production simulations.
----------
.. _Typelabel2:
**(Gissinger)** J. R. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M. Choi, D. S. Karls, A. Stukowski, W. Im, H. Heinz, A. Kohlmeyer, and E. B. Tadmor, J Phys Chem B, 128, 3282-3297 (2024).

38
doc/src/Howto_bioFF.rst
View File

@ -1,22 +1,16 @@
CHARMM, AMBER, COMPASS, DREIDING, and OPLS force fields
=======================================================
A compact summary of the concepts, definitions, and properties of
force fields with explicit bonded interactions (like the ones discussed
in this HowTo) is given in :ref:`(Gissinger) <Typelabel2>`.
A force field has 2 parts: the formulas that define it and the
coefficients used for a particular system. Here we only discuss
formulas implemented in LAMMPS that correspond to formulas commonly used
in the CHARMM, AMBER, COMPASS, and DREIDING force fields. Setting
coefficients is done either from special sections in an input data file
via the :doc:`read_data <read_data>` command or in the input script with
commands like :doc:`pair_coeff <pair_coeff>` or :doc:`bond_coeff
<bond_coeff>` and so on. See the :doc:`Tools <Tools>` doc page for
additional tools that can use CHARMM, AMBER, or Materials Studio
generated files to assign force field coefficients and convert their
output into LAMMPS input. LAMMPS input scripts can also be generated by
`charmm-gui.org <https://charmm-gui.org/>`_.
Here we only discuss formulas implemented in LAMMPS that correspond to
formulas commonly used in the CHARMM, AMBER, COMPASS, and DREIDING force
fields. Setting coefficients is done either from special sections in an
input data file via the :doc:`read_data <read_data>` command or in the
input script with commands like :doc:`pair_coeff <pair_coeff>` or
:doc:`bond_coeff <bond_coeff>` and so on. See the :doc:`Tools <Tools>`
doc page for additional tools that can use CHARMM, AMBER, or Materials
Studio generated files to assign force field coefficients and convert
their output into LAMMPS input. LAMMPS input scripts can also be
generated by `charmm-gui.org <https://charmm-gui.org/>`_.
CHARMM and AMBER
----------------
@ -203,9 +197,11 @@ rather than individual force constants and geometric parameters that
depend on the particular combinations of atoms involved in the bond,
angle, or torsion terms. DREIDING has an :doc:`explicit hydrogen bond
term <pair_hbond_dreiding>` to describe interactions involving a
hydrogen atom on very electronegative atoms (N, O, F). Unlike CHARMM
or AMBER, the DREIDING force field has not been parameterized for
considering solvents (like water).
hydrogen atom on very electronegative atoms (N, O, F). Unlike CHARMM or
AMBER, the DREIDING force field has not been parameterized for
considering solvents (like water) and has no rules for assigning
(partial) charges. That will seriously limit its accuracy when used for
simulating systems where those matter.
See :ref:`(Mayo) <howto-Mayo>` for a description of the DREIDING force field
@ -272,10 +268,6 @@ compatible with a subset of OPLS interactions.
----------
.. _Typelabel2:
**(Gissinger)** J. R. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M. Choi, D. S. Karls, A. Stukowski, W. Im, H. Heinz, A. Kohlmeyer, and E. B. Tadmor, J Phys Chem B, 128, 3282-3297 (2024).
.. _howto-MacKerell:
**(MacKerell)** MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field, Fischer, Gao, Guo, Ha, et al (1998). J Phys Chem, 102, 3586 . https://doi.org/10.1021/jp973084f

5
doc/src/Howto_lammps_gui.rst
View File

@ -308,7 +308,10 @@ of the *Output* window showing how many warnings and errors were
detected and how many lines the entire output has. By clicking on the
button on the right with the warning symbol or by using the keyboard
shortcut `Ctrl-N` (`Command-N` on macOS), you can jump to the next
line with a warning or error.
line with a warning or error. If there is a URL pointing to additional
explanations in the online manual, that URL will be highlighted and
double-clicking on it shall open the corresponding manual page in
the web browser. The option is also available from the context menu.
By default, the *Output* window is replaced each time a run is started.
The runs are counted and the run number for the current run is displayed

4
doc/src/Howto_spc.rst
View File

@ -1,5 +1,5 @@
SPC water model
===============
SPC and SPC/E water model
=========================
The SPC water model specifies a 3-site rigid water molecule with
charges and Lennard-Jones parameters assigned to each of the three atoms.

9
doc/src/Howto_thermostat.rst
View File

@ -21,9 +21,14 @@ can be invoked via the *dpd/tstat* pair style:
* :doc:`fix nvt/sllod <fix_nvt_sllod>`
* :doc:`fix temp/berendsen <fix_temp_berendsen>`
* :doc:`fix temp/csvr <fix_temp_csvr>`
* :doc:`fix ffl <fix_ffl>`
* :doc:`fix gjf <fix_gjf>`
* :doc:`fix gld <fix_gld>`
* :doc:`fix gle <fix_gle>`
* :doc:`fix langevin <fix_langevin>`
* :doc:`fix temp/rescale <fix_temp_rescale>`
* :doc:`pair_style dpd/tstat <pair_dpd>`
* :doc:`pair_style dpd/ext/tstat <pair_dpd_ext>`
:doc:`Fix nvt <fix_nh>` only thermostats the translational velocity of
particles. :doc:`Fix nvt/sllod <fix_nvt_sllod>` also does this,
@ -82,10 +87,10 @@ that:
.. note::
Only the nvt fixes perform time integration, meaning they update
Not all thermostat fixes perform time integration, meaning they update
the velocities and positions of particles due to forces and velocities
respectively. The other thermostat fixes only adjust velocities; they
do NOT perform time integration updates. Thus they should be used in
do NOT perform time integration updates. Thus, they should be used in
conjunction with a constant NVE integration fix such as these:
* :doc:`fix nve <fix_nve>`

129
doc/src/Howto_tip4p.rst
View File

@ -1,5 +1,5 @@
TIP4P water model
=================
TIP4P and OPC water models
==========================
The four-point TIP4P rigid water model extends the traditional
:doc:`three-point TIP3P <Howto_tip3p>` model by adding an additional
@ -9,9 +9,11 @@ the oxygen along the bisector of the HOH bond angle. A bond style of
:doc:`harmonic <bond_harmonic>` and an angle style of :doc:`harmonic
<angle_harmonic>` or :doc:`charmm <angle_charmm>` should also be used.
In case of rigid bonds also bond style :doc:`zero <bond_zero>` and angle
style :doc:`zero <angle_zero>` can be used.
style :doc:`zero <angle_zero>` can be used. Very similar to the TIP4P
model is the OPC water model. It can be realized the same way as TIP4P
but has different geometry and force field parameters.
There are two ways to implement TIP4P water in LAMMPS:
There are two ways to implement TIP4P-like water in LAMMPS:
#. Use a specially written pair style that uses the :ref:`TIP3P geometry
<tip3p_molecule>` without the point M. The point M location is then
@ -21,7 +23,10 @@ There are two ways to implement TIP4P water in LAMMPS:
computationally very efficient, but the charge distribution in space
is only correct within the tip4p labeled styles. So all other
computations using charges will "see" the negative charge incorrectly
on the oxygen atom.
located on the oxygen atom unless they are specially written for using
the TIP4P geometry internally as well, e.g. :doc:`compute dipole/tip4p
<compute_dipole>`, :doc:`fix efield/tip4p <fix_efield>`, or
:doc:`kspace_style pppm/tip4p <kspace_style>`.
This can be done with the following pair styles for Coulomb with a cutoff:
@ -68,77 +73,90 @@ TIP4P/2005 model :ref:`(Abascal2) <Abascal2>` and a version of TIP4P
parameters adjusted for use with a long-range Coulombic solver
(e.g. Ewald or PPPM in LAMMPS). Note that for implicit TIP4P models the
OM distance is specified in the :doc:`pair_style <pair_style>` command,
not as part of the pair coefficients.
not as part of the pair coefficients. Also parameters for the OPC
model (:ref:`Izadi <Izadi>`) are provided.
.. list-table::
:header-rows: 1
:widths: 36 19 13 15 17
:widths: 40 12 12 14 11 11
* - Parameter
- TIP4P (original)
- TIP4P/Ice
- TIP4P/2005
- TIP4P (Ewald)
- OPC
* - O mass (amu)
- 15.9994
- 15.9994
- 15.9994
- 15.9994
- 15.9994
* - H mass (amu)
- 1.008
- 1.008
- 1.008
- 1.008
- 1.008
* - O or M charge (:math:`e`)
- -1.040
- -1.1794
- -1.1128
- -1.04844
- -1.3582
* - H charge (:math:`e`)
- 0.520
- 0.5897
- 0.5564
- 0.52422
- 0.6791
* - LJ :math:`\epsilon` of OO (kcal/mole)
- 0.1550
- 0.21084
- 0.1852
- 0.16275
- 0.21280
* - LJ :math:`\sigma` of OO (:math:`\AA`)
- 3.1536
- 3.1668
- 3.1589
- 3.16435
- 3.1660
* - LJ :math:`\epsilon` of HH, MM, OH, OM, HM (kcal/mole)
- 0.0
- 0.0
- 0.0
- 0.0
- 0.0
* - LJ :math:`\sigma` of HH, MM, OH, OM, HM (:math:`\AA`)
- 1.0
- 1.0
- 1.0
- 1.0
- 1.0
* - :math:`r_0` of OH bond (:math:`\AA`)
- 0.9572
- 0.9572
- 0.9572
- 0.9572
- 0.8724
* - :math:`\theta_0` of HOH angle
- 104.52\ :math:`^{\circ}`
- 104.52\ :math:`^{\circ}`
- 104.52\ :math:`^{\circ}`
- 104.52\ :math:`^{\circ}`
- 103.60\ :math:`^{\circ}`
* - OM distance (:math:`\AA`)
- 0.15
- 0.1577
- 0.1546
- 0.1250
- 0.1594
Note that the when using the TIP4P pair style, the neighbor list cutoff
Note that the when using a TIP4P pair style, the neighbor list cutoff
for Coulomb interactions is effectively extended by a distance 2 \* (OM
distance), to account for the offset distance of the fictitious charges
on O atoms in water molecules. Thus it is typically best in an
on O atoms in water molecules. Thus, it is typically best in an
efficiency sense to use a LJ cutoff >= Coulomb cutoff + 2\*(OM
distance), to shrink the size of the neighbor list. This leads to
slightly larger cost for the long-range calculation, so you can test the
@ -192,6 +210,94 @@ file changed):
run 20000
write_data tip4p-implicit.data nocoeff
When constructing an OPC model, we cannot use the ``tip3p.mol`` file due
to the different geometry. Below is a molecule file providing the 3
sites of an implicit OPC geometry for use with TIP4P styles. Note, that
the "Shake" and "Special" sections are missing here. Those will be
auto-generated by LAMMPS when the molecule file is loaded *after* the
simulation box has been created. These sections are required only when
the molecule file is loaded *before*.
.. _opc3p_molecule:
.. code-block::
# Water molecule. 3 point geometry for OPC model
3 atoms
2 bonds
1 angles
Coords
1 0.00000 -0.06037 0.00000
2 0.68558 0.50250 0.00000
3 -0.68558 0.50250 0.00000
Types
1 1 # O
2 2 # H
3 2 # H
Charges
1 -1.3582
2 0.6791
3 0.6791
Bonds
1 1 1 2
2 1 1 3
Angles
1 1 2 1 3
Below is a LAMMPS input file using the implicit method to implement
the OPC model using the molecule file from above and including the
PPPM long-range Coulomb solver.
.. code-block:: LAMMPS
units real
atom_style full
region box block -5 5 -5 5 -5 5
create_box 2 box bond/types 1 angle/types 1 &
extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2
mass 1 15.9994
mass 2 1.008
pair_style lj/cut/tip4p/long 1 2 1 1 0.1594 12.0
pair_coeff 1 1 0.2128 3.166
pair_coeff 2 2 0.0 1.0
bond_style zero
bond_coeff 1 0.8724
angle_style zero
angle_coeff 1 103.6
kspace_style pppm/tip4p 1.0e-5
molecule water opc3p.mol # this file has the OPC geometry but is without M
create_atoms 0 random 33 34564 NULL mol water 25367 overlap 1.33
fix rigid all shake 0.001 10 10000 b 1 a 1
minimize 0.0 0.0 1000 10000
reset_timestep 0
timestep 1.0
velocity all create 300.0 5463576
fix integrate all nvt temp 300 300 100.0
thermo_style custom step temp press etotal pe
thermo 1000
run 20000
write_data opc-implicit.data nocoeff
Below is the code for a LAMMPS input file using the explicit method and
a TIP4P molecule file. Because of using :doc:`fix rigid/small
<fix_rigid>` no bonds need to be defined and thus no extra storage needs
@ -279,3 +385,8 @@ Phys, 79, 926 (1983).
**(Abascal2)** Abascal, J Chem Phys, 123, 234505 (2005)
https://doi.org/10.1063/1.2121687
.. _Izadi:
**(Izadi)** Izadi, Anandakrishnan, Onufriev, J. Phys. Chem. Lett., 5, 21, 3863 (2014)
https://doi.org/10.1021/jz501780a

5
doc/src/Intro_authors.rst
View File

@ -84,8 +84,9 @@ lammps.org". General questions about LAMMPS should be posted in the
\normalsize
Past developers include Paul Crozier and Mark Stevens, both at SNL,
and Ray Shan, now at Materials Design.
Past core developers include Paul Crozier and Mark Stevens, both at SNL,
and Ray Shan while at SNL and later at Materials Design, now at Thermo
Fisher Scientific.
----------

25
doc/src/Intro_portability.rst
View File

@ -28,8 +28,9 @@ Build systems
LAMMPS can be compiled from source code using a (traditional) build
system based on shell scripts, a few shell utilities (grep, sed, cat,
tr) and the GNU make program. This requires running within a Bourne
shell (``/bin/sh``). Alternatively, a build system with different back
ends can be created using CMake. CMake must be at least version 3.16.
shell (``/bin/sh`` or ``/bin/bash``). Alternatively, a build system
with different back ends can be created using CMake. CMake must be
at least version 3.16.
Operating systems
^^^^^^^^^^^^^^^^^
@ -40,11 +41,18 @@ Also, compilation and correct execution on macOS and Windows (using
Microsoft Visual C++) is checked automatically for the largest part of
the source code. Some (optional) features are not compatible with all
operating systems, either through limitations of the corresponding
LAMMPS source code or through incompatibilities of source code or
build system of required external libraries or packages.
LAMMPS source code or through incompatibilities or build system
limitations of required external libraries or packages.
Executables for Windows may be created natively using either Cygwin or
Visual Studio or with a Linux to Windows MinGW cross-compiler.
Executables for Windows may be created either natively using Cygwin,
MinGW, Intel, Clang, or Microsoft Visual C++ compilers, or with a Linux
to Windows MinGW cross-compiler. Native compilation is supported using
Microsoft Visual Studio or a terminal window (using the CMake build
system).
Executables for macOS may be created either using Xcode or GNU compilers
installed with Homebrew. In the latter case, building of LAMMPS through
Homebrew instead of a manual compile is also possible.
Additionally, FreeBSD and Solaris have been tested successfully to
run LAMMPS and produce results consistent with those on Linux.
@ -61,8 +69,9 @@ CPU architectures
^^^^^^^^^^^^^^^^^
The primary CPU architecture for running LAMMPS is 64-bit x86, but also
32-bit x86, and 64-bit ARM and PowerPC (64-bit, Little Endian) are
regularly tested.
64-bit ARM and PowerPC (64-bit, Little Endian) are currently regularly
tested. Further architectures are tested by Linux distributions that
bundle LAMMPS.
Portability compliance
^^^^^^^^^^^^^^^^^^^^^^

BIN
doc/src/JPG/lammps-gui-screen.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 361 KiB

BIN
doc/src/JPG/lammps-releases.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 68 KiB

After

Width:  |  Height:  |  Size: 70 KiB

4
doc/src/Modify_pair.rst
View File

@ -46,6 +46,8 @@ Here is a brief list of some the class methods in the Pair class that
+---------------------------------+------------------------------------------------------------------------+
| compute_inner/middle/outer | versions of compute used by rRESPA |
+---------------------------------+------------------------------------------------------------------------+
| compute_atomic_energy | energy of one atom, equivalent to per-atom energy |
+---------------------------------+------------------------------------------------------------------------+
| memory_usage | return estimated amount of memory used by the pair style |
+---------------------------------+------------------------------------------------------------------------+
| modify_params | process arguments to pair_modify command |
@ -122,3 +124,5 @@ setting.
+---------------------------------+-------------------------------------------------------------+---------+
| spinflag | 1 if compatible with spin kspace_style | 0 |
+---------------------------------+-------------------------------------------------------------+---------+
| atomic_energy_enable | 1 if compute_atomic_energy() routine exists | 0 |
+---------------------------------+-------------------------------------------------------------+---------+

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

Binary file not shown.

16
doc/src/Tools.rst
View File

@ -475,9 +475,13 @@ beginners to start with LAMMPS, it is also the expectation that
LAMMPS-GUI users will eventually transition to workflows that most
experienced LAMMPS users employ.
All features have been extensively exposed to keyboard shortcuts, so
that there is also appeal for experienced LAMMPS users for prototyping
and testing simulation setups.
.. image:: JPG/lammps-gui-screen.png
:align: center
:scale: 50%
Features have been extensively exposed to keyboard shortcuts, so that
there is also appeal for experienced LAMMPS users for prototyping and
testing simulation setups.
Features
^^^^^^^^
@ -502,7 +506,7 @@ Here are a few highlights of LAMMPS-GUI
- Visualization of current state in Image Viewer (via calling :doc:`write_dump image <dump_image>`)
- Capture of images created via :doc:`dump image <dump_image>` in Slide show window
- Dialog to set variables, similar to the LAMMPS command-line flag '-v' / '-var'
- Support for GPU, INTEL, KOKKOS/OpenMP, OPENMAP, and OPT and accelerator packages
- Support for GPU, INTEL, KOKKOS/OpenMP, OPENMP, and OPT accelerator packages
Parallelization
^^^^^^^^^^^^^^^
@ -523,8 +527,8 @@ with CMake is required.
The LAMMPS-GUI has been successfully compiled and tested on:
- Ubuntu Linux 20.04LTS x86_64 using GCC 9, Qt version 5.12
- Fedora Linux 40 x86\_64 using GCC 14 and Clang 17, Qt version 5.15LTS
- Fedora Linux 40 x86\_64 using GCC 14, Qt version 6.7
- Fedora Linux 41 x86\_64 using GCC 14 and Clang 17, Qt version 5.15LTS
- Fedora Linux 41 x86\_64 using GCC 14, Qt version 6.8
- Apple macOS 12 (Monterey) and macOS 13 (Ventura) with Xcode on arm64 and x86\_64, Qt version 5.15LTS
- Windows 10 and 11 x86_64 with Visual Studio 2022 and Visual C++ 14.36, Qt version 5.15LTS
- Windows 10 and 11 x86_64 with Visual Studio 2022 and Visual C++ 14.40, Qt version 6.7

30
doc/src/compute_bond_local.rst
View File

@ -64,20 +64,32 @@ All these properties are computed for the pair of atoms in a bond,
whether the two atoms represent a simple diatomic molecule, or are part
of some larger molecule.
The value *dist* is the current length of the bond.
The values *dx*, *dy*, and *dz* are the xyz components of the
*distance* between the pair of atoms. This value is always the
distance from the atom of lower to the one with the higher id.
.. versionchanged:: TBD
The sign of *dx*, *dy*, *dz* is no longer determined by the atom IDs
of the bonded atoms but by their order in the bond list to be
consistent with *fx*, *fy*, and *fz*.
The value *dist* is the current length of the bond. The values *dx*,
*dy*, and *dz* are the :math:`(x,y,z)` components of the distance vector
:math:`\vec{x_i} - \vec{x_j}` between the atoms in the bond. The order
of the atoms is determined by the bond list and the respective atom-IDs
can be output with :doc:`compute property/local
<compute_property_local>`.
The value *engpot* is the potential energy for the bond,
based on the current separation of the pair of atoms in the bond.
The value *force* is the magnitude of the force acting between the
pair of atoms in the bond.
The value *force* is the magnitude of the force acting between the pair
of atoms in the bond, which is positive for a repulsive force and
negative for an attractive force.
The values *fx*, *fy*, and *fz* are the xyz components of
*force* between the pair of atoms in the bond. For bond styles that apply
non-central forces, such as :doc:`bond_style bpm/rotational
The values *fx*, *fy*, and *fz* are the :math:`(x,y,z)` components of
the force on the first atom *i* in the bond due to the second atom *j*.
Mathematically, they are obtained by multiplying the value of *force*
from above with a unit vector created from the *dx*, *dy*, and *dz*
components of the distance vector also described above. For bond styles
that apply non-central forces, such as :doc:`bond_style bpm/rotational
<bond_bpm_rotational>`, these values only include the :math:`(x,y,z)`
components of the normal force component.

30
doc/src/compute_pair_local.rst
View File

@ -56,19 +56,33 @@ force cutoff distance for that interaction, as defined by the
:doc:`pair_style <pair_style>` and :doc:`pair_coeff <pair_coeff>`
commands.
The value *dist* is the distance between the pair of atoms.
The values *dx*, *dy*, and *dz* are the :math:`(x,y,z)` components of the
*distance* between the pair of atoms. This value is always the
distance from the atom of higher to the one with the lower atom ID.
.. versionchanged:: TBD
The sign of *dx*, *dy*, *dz* is no longer determined by the value of
their atom-IDs but by their order in the neighbor list to be
consistent with *fx*, *fy*, and *fz*.
The value *dist* is the distance between the pair of atoms. The values
*dx*, *dy*, and *dz* are the :math:`(x,y,z)` components of the distance
vector :math:`\vec{x_i} - \vec{x_j}` between the pair of atoms. The
order of the atoms is determined by the neighbor list and the respective
atom-IDs can be output with :doc:`compute property/local
<compute_property_local>`.
The value *eng* is the interaction energy for the pair of atoms.
The value *force* is the force acting between the pair of atoms, which
is positive for a repulsive force and negative for an attractive
force. The values *fx*, *fy*, and *fz* are the :math:`(x,y,z)` components of
*force* on atom I. For pair styles that apply non-central forces,
such as :doc:`granular pair styles <pair_gran>`, these values only include
the :math:`(x,y,z)` components of the normal force component.
force.
The values *fx*, *fy*, and *fz* are the :math:`(x,y,z)` components of
the force vector on the first atom *i* of a pair in the neighbor list
due to the second atom *j*. Mathematically, they are obtained by
multiplying the value of *force* from above with a unit vector created
from the *dx*, *dy*, and *dz* components of the distance vector also
described above. For pair styles that apply non-central forces, such as
:doc:`granular pair styles <pair_gran>`, these values only include the
:math:`(x,y,z)` components of the normal force component.
A pair style may define additional pairwise quantities which can be
accessed as *p1* to *pN*, where :math:`N` is defined by the pair style.

2
doc/src/fix.rst
View File

@ -256,6 +256,7 @@ accelerated styles exist.
* :doc:`flow/gauss <fix_flow_gauss>` - Gaussian dynamics for constant mass flux
* :doc:`freeze <fix_freeze>` - freeze atoms in a granular simulation
* :doc:`gcmc <fix_gcmc>` - grand canonical insertions/deletions
* :doc:`gjf <fix_gjf>` - statistically correct Langevin temperature control using the GJ methods
* :doc:`gld <fix_gld>` - generalized Langevin dynamics integrator
* :doc:`gle <fix_gle>` - generalized Langevin equation thermostat
* :doc:`gravity <fix_gravity>` - add gravity to atoms in a granular simulation
@ -290,6 +291,7 @@ accelerated styles exist.
* :doc:`mvv/tdpd <fix_mvv_dpd>` - constant temperature DPD using the modified velocity-Verlet algorithm
* :doc:`neb <fix_neb>` - nudged elastic band (NEB) spring forces
* :doc:`neb/spin <fix_neb_spin>` - nudged elastic band (NEB) spring forces for spins
* :doc:`neighbor/swap <fix_neighbor_swap>` - kinetic Monte Carlo (kMC) atom swapping
* :doc:`nonaffine/displacement <fix_nonaffine_displacement>` - calculate nonaffine displacement of atoms
* :doc:`nph <fix_nh>` - constant NPH time integration via Nose/Hoover
* :doc:`nph/asphere <fix_nph_asphere>` - NPH for aspherical particles

208
doc/src/fix_gjf.rst Normal file
View File

@ -0,0 +1,208 @@
.. index:: fix gjf
fix gjf command
========================
Syntax
""""""
.. code-block:: LAMMPS
fix ID group-ID gjf Tstart Tstop damp seed keyword values ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* gjf = style name of this fix command
* Tstart,Tstop = desired temperature at start/end of run (temperature units)
* Tstart can be a variable (see below)
* damp = damping parameter (time units)
* seed = random number seed to use for white noise (positive integer)
* zero or more keyword/value pairs may be appended
* keyword = *vel* or *method*
.. parsed-literal::
*vel* value = *vfull* or *vhalf*
*vfull* = use on-site velocity
*vhalf* = use half-step velocity
*method* value = *1-8*
*1-8* = choose one of the many GJ formulations
*7* = requires input of additional scalar between 0 and 1
Examples
""""""""
.. code-block:: LAMMPS
fix 3 boundary gjf 10.0 10.0 1.0 699483
fix 1 all gjf 10.0 100.0 100.0 48279 vel vfull method 4
fix 2 all gjf 10.0 10.0 1.0 26488 method 7 0.95
Description
"""""""""""
.. versionadded:: TBD
Apply a Langevin thermostat as described in :ref:`(Gronbech-Jensen-2020) <Gronbech-Jensen-2020>`
to a group of atoms which models an interaction with a background
implicit solvent. As described in the papers cited below, the GJ methods
provide exact diffusion, drift, and Boltzmann sampling for linear systems for
any time step within the stability limit. The purpose of this set of methods
is therefore to significantly improve statistical accuracy at longer time steps
compared to other thermostats.
The current implementation provides the user with the option to output
the velocity in one of two forms: *vfull* or *vhalf*. The option *vhalf*
outputs the 2GJ half-step velocity given in :ref:`Gronbech Jensen/Gronbech-Jensen
<Gronbech-Jensen-2019>`; for linear systems, this velocity is shown to not
have any statistical errors for any stable time step. The option *vfull*
outputs the on-site velocity given in :ref:`Gronbech-Jensen/Farago
<Gronbech-Jensen-Farago>`; this velocity is shown to be systematically lower
than the target temperature by a small amount, which grows
quadratically with the timestep. An overview of statistically correct Boltzmann
and Maxwell-Boltzmann sampling of true on-site and true half-step velocities is
given in :ref:`Gronbech-Jensen-2020 <Gronbech-Jensen-2020>`.
This fix allows the use of several GJ methods as listed in :ref:`Gronbech-Jensen-2020 <Gronbech-Jensen-2020>`.
The GJ-VII method is described in :ref:`Finkelstein <Finkelstein>` and GJ-VIII
is described in :ref:`Gronbech-Jensen-2024 <Gronbech-Jensen-2024>`.
The implementation follows the splitting form provided in Eqs. (24) and (25)
in :ref:`Gronbech-Jensen-2024 <Gronbech-Jensen-2024>`, including the application
of Gaussian noise values, per the description in
:ref:`Gronbech-Jensen-2023 <Gronbech-Jensen-2023>`.
.. note::
Unlike the :doc:`fix langevin <fix_langevin>` command which performs force
modifications only, this fix performs thermostatting and time integration.
Thus you no longer need a separate time integration fix, like :doc:`fix nve <fix_nve>`.
See the :doc:`Howto thermostat <Howto_thermostat>` page for
a discussion of different ways to compute temperature and perform
thermostatting.
The desired temperature at each timestep is a ramped value during the
run from *Tstart* to *Tstop*\ .
*Tstart* can be specified as an equal-style or atom-style
:doc:`variable <variable>`. In this case, the *Tstop* setting is
ignored. If the value is a variable, it should be specified as
v_name, where name is the variable name. In this case, the variable
will be evaluated each timestep, and its value used to determine the
target temperature.
Equal-style variables can specify formulas with various mathematical
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.
Atom-style variables can specify the same formulas as equal-style
variables but can also include per-atom values, such as atom
coordinates. Thus it is easy to specify a spatially-dependent
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. 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.
The *damp* parameter is specified in time units and determines how
rapidly the temperature is relaxed. For example, a value of 100.0 means
to relax the temperature in a timespan of (roughly) 100 time units
(:math:`\tau` or fs or ps - see the :doc:`units <units>` command). The
damp factor can be thought of as inversely related to the viscosity of
the solvent. I.e. a small relaxation time implies a high-viscosity
solvent and vice versa. See the discussion about :math:`\gamma` and
viscosity in the documentation for the :doc:`fix viscous <fix_viscous>`
command for more details.
The random # *seed* must be a positive integer. A Marsaglia random
number generator is used. Each processor uses the input seed to
generate its own unique seed and its own stream of random numbers.
Thus the dynamics of the system will not be identical on two runs on
different numbers of processors.
----------
The keyword/value option pairs are used in the following ways.
The keyword *vel* determines which velocity is used to determine
quantities of interest in the simulation.
The keyword *method* selects one of the eight GJ-methods implemented in LAMMPS.
----------
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files <restart>`.
Because the state of the random number generator is not saved in restart files,
this means you cannot do "exact" restarts with this fix, where the simulation
continues on the same as if no restart had taken place. However, in a
statistical sense, a restarted simulation should produce the same behavior.
Additionally, the GJ methods implement noise exclusively within each time step
(unlike the BBK thermostat of the fix-langevin). The restart is done with
either vfull or vhalf velocity output for as long as the choice of vfull/vhalf
is the same for the simulation as it is in the restart file.
The :doc:`fix_modify <fix_modify>` *temp* option is supported by this
fix. You can use it to assign a temperature :doc:`compute <compute>`
you have defined to this fix which will be used in its thermostatting
procedure, as described above. For consistency, the group used by
this fix and by the compute should be the same.
This fix can ramp its target temperature over multiple runs, using the
*start* and *stop* keywords of the :doc:`run <run>` command. See the
:doc:`run <run>` command for details of how to do this.
This fix is not invoked during :doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
This fix is not compatible with run_style respa. It is not compatible with
accelerated packages such as KOKKOS.
Related commands
""""""""""""""""
:doc:`fix langevin <fix_langevin>`, :doc:`fix nvt <fix_nh>`
Default
"""""""
The option defaults are vel = vhalf, method = 1.
----------
.. _Gronbech-Jensen-2020:
**(Gronbech-Jensen-2020)** Gronbech-Jensen, Mol Phys 118, e1662506 (2020).
.. _Gronbech-Jensen-2019:
**(Gronbech Jensen/Gronbech-Jensen)** Gronbech Jensen and Gronbech-Jensen, Mol Phys, 117, 2511 (2019)
.. _Gronbech-Jensen-Farago:
**(Gronbech-Jensen/Farago)** Gronbech-Jensen and Farago, Mol Phys, 111, 983 (2013).
.. _Finkelstein:
**(Finkelstein)** Finkelstein, Cheng, Florin, Seibold, Gronbech-Jensen, J. Chem. Phys., 155, 18 (2021)
.. _Gronbech-Jensen-2024:
**(Gronbech-Jensen-2024)** Gronbech-Jensen, J. Stat. Phys. 191, 137 (2024).
.. _Gronbech-Jensen-2023:
**(Gronbech-Jensen-2023)** Gronbech-Jensen, J. Stat. Phys. 190, 96 (2023).

19
doc/src/fix_langevin.rst
View File

@ -56,7 +56,7 @@ Examples
Description
"""""""""""
Apply a Langevin thermostat as described in :ref:`(Schneider) <Schneider1>`
Apply a Langevin thermostat as described in :ref:`(Bruenger) <Bruenger1>`
to a group of atoms which models an interaction with a background
implicit solvent. Used with :doc:`fix nve <fix_nve>`, this command
performs Brownian dynamics (BD), since the total force on each atom
@ -241,6 +241,13 @@ to zero by subtracting off an equal part of it from each atom in the
group. As a result, the center-of-mass of a system with zero initial
momentum will not drift over time.
.. deprecated:: TDB
The *gjf* keyword in fix langevin is deprecated and will be removed
soon. The GJF functionality has been moved to its own fix style
:doc:`fix gjf <fix_gjf>` and it is strongly recommended to use that
fix instead.
The keyword *gjf* can be used to run the :ref:`Gronbech-Jensen/Farago
<Gronbech-Jensen>` time-discretization of the Langevin model. As
described in the papers cited below, the purpose of this method is to
@ -324,14 +331,16 @@ types, tally = no, zero = no, gjf = no.
----------
.. _Bruenger1:
**(Bruenger)** Bruenger, Brooks, and Karplus, Chem. Phys. Lett. 105, 495 (1982).
[Previously attributed to Schneider and Stoll, Phys. Rev. B 17, 1302 (1978).
Implementation remains unchanged.]
.. _Dunweg1:
**(Dunweg)** Dunweg and Paul, Int J of Modern Physics C, 2, 817-27 (1991).
.. _Schneider1:
**(Schneider)** Schneider and Stoll, Phys Rev B, 17, 1302 (1978).
.. _Gronbech-Jensen:
**(Gronbech-Jensen)** Gronbech-Jensen and Farago, Mol Phys, 111, 983

264
doc/src/fix_neighbor_swap.rst Normal file
View File

@ -0,0 +1,264 @@
.. index:: fix neighbor/swap
fix neighbor/swap command
=========================
Syntax
""""""
.. code-block:: LAMMPS
fix ID group-ID neighbor/swap N X seed T R0 voro-ID keyword values ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* neighbor/swap = style name of this fix command
* N = invoke this fix every N steps
* X = number of swaps to attempt every N steps
* seed = random # seed (positive integer)
* T = scaling temperature of the MC swaps (temperature units)
* R0 = scaling swap probability of the MC swaps (distance units)
* voro-ID = valid voronoi compute id (compute voronoi/atom)
* one or more keyword/value pairs may be appended to args
* keywords *types* and *diff* are mutually exclusive, but one must be specified
* keyword = *types* or *diff* or *ke* or *region* or *rates*
.. parsed-literal::
*types* values = two or more atom types (Integers in range [1,Ntypes] or type labels)
*diff* values = one atom type
*ke* value = *yes* or *no*
*yes* = kinetic energy is conserved after atom swaps
*no* = no conservation of kinetic energy after atom swaps
*region* value = region-ID
region-ID = ID of region to use as an exchange/move volume
*rates* values = V1 V2 . . . Vntypes values to conduct variable diffusion for different atom types (unitless)
Examples
""""""""
.. code-block:: LAMMPS
compute voroN all voronoi/atom neighbors yes
fix mc all neighbor/swap 10 160 15238 1000.0 3.0 voroN diff 2
fix myFix all neighbor/swap 100 1 12345 298.0 3.0 voroN region my_swap_region types 5 6
fix kmc all neighbor/swap 1 100 345 1.0 3.0 voroN diff 3 rates 3 1 6
Description
"""""""""""
.. versionadded:: TBD
This fix performs Monte-Carlo (MC) evaluations to enable kinetic
Monte Carlo (kMC)-type behavior during MD simulation by allowing
neighboring atoms to swap their positions. In contrast to the :doc:`fix
atom/swap <fix_atom_swap>` command which swaps pairs of atoms anywhere
in the simulation domain, the restriction of the MC swapping to
neighbors enables a hybrid MD/kMC-like simulation.
Neighboring atoms are defined by using a Voronoi tesselation performed
by the :doc:`compute voronoi/atom <compute_voronoi_atom>` command.
Two atoms are neighbors if their Voronoi cells share a common face
(3d) or edge (2d).
The selection of a swap neighbor is made using a distance-based
criterion for weighting the selection probability of each swap, in the
same manner as kMC selects a next event using relative probabilities.
The acceptance or rejection of each swap is determined via the
Metropolis criterion after evaluating the change in system energy due
to the swap.
A detailed explanation of the original implementation of this
algorithm can be found in :ref:`(Tavenner 2023) <TavennerMDkMC>`
where it was used to simulated accelerated diffusion in an MD context.
Simulating inherently kinetically-limited behaviors which rely on rare
events (such as atomic diffusion in a solid) is challenging for
traditional MD since its relatively short timescale will not naturally
sample many events. This fix addresses this challenge by allowing rare
neighbor hopping events to be sampled in a kMC-like fashion at a much
faster rate (set by the specified *N* and *X* parameters). This enables
the processes of atomic diffusion to be approximated during an MD
simulation, effectively decoupling the MD atomic vibrational timescale
and the atomic hopping (kMC event) timescale.
The algorithm implemented by this fix is as follows:
- The MD simulation is paused every *N* steps
- A Voronoi tesselation is performed for the current atom configuration.
- Then *X* atom swaps are attempted, one after the other.
- For each swap, an atom *I* is selected randomly from the list of
atom types specified by either the *types* or *diff* keywords.
- One of *I*'s Voronoi neighbors *J* is selected using the
distance-weighted probability for each neighbor detailed below.
- The *I,J* atom IDs are communicated to all processors so that a
global energy evaluation can be performed for the post-swap state
of the system.
- The swap is accepted or rejected based on the Metropolis criterion
using the energy change of the system and the specified temperature
*T*.
Here are a few comments on the computational cost of the swapping
algorithm.
1. The cost of a global energy evaluation is similar to that of an MD
timestep.
2. Similar to other MC algorithms in LAMMPS, improved parallel
efficiency is achieved with a smaller number of atoms per
processor than would typically be used in an standard MD
simulation. This is because the per-energy evaluation cost
increases relative to the balance of MD/MC steps as indicated by
1., but the communication cost remains relatively constant for a
given number of MD steps.
3. The MC portion of the simulation will run dramatically slower if
the pair style uses different cutoffs for different atom types (or
type pairs). This is because each atom swap then requires a
rebuild of the neighbor list to ensure the post-swap global energy
can be computed correctly.
Limitations are imposed on selection of *I,J* atom pairs to avoid
swapping of atoms which are outside of a reasonable cutoff (e.g. due to
a Voronoi tesselation near free surfaces) though the use of a
distance-weighted probability scaling.
----------
This section gives more details on other arguments and keywords.
The random number generator (RNG) used by all the processors for MC
operations is initialized with the specified *seed*.
The distance-based probability is weighted by the specified *R0* which
sets the radius :math:`r_0` in this formula
.. math::
p_{ij} = e^{(\frac{r_{ij}}{r_0})^2}
where :math:`p_{ij}` is the probability of selecting atom :math:`j` to
swap with atom :math:`i`. Typically, a value for *R0* around the
average nearest-neighbor spacing is appropriate. Since this is simply a
probability weighting, the swapping behavior is not very sensitive to
the exact value of *R0*.
The required *voro-ID* value is the compute-ID of a
:doc:`compute voronoi/atom <compute_voronoi_atom>` command like
this:
.. code-block:: LAMMPS
compute compute-ID group-ID voronoi/atom neighbors yes
It must return per-atom list of valid neighbor IDs as in the
:doc:`compute voronoi/atom <compute_voronoi_atom>` command.
The keyword *types* takes two or more atom types as its values. Only
atoms *I* of the first atom type will be selected. Only atoms *J* of the
remaining atom types will be considered as potential swap partners.
The keyword *diff* take a single atom type as its value. Only atoms
*I* of the that atom type will be selected. Atoms *J* of all
remaining atom types will be considered as potential swap partners.
This includes the atom type specified with the *diff* keyword to
account for self-diffusive hops between two atoms of the same type.
Note that the *neighbors yes* option must be enabled for use with this
fix. The group-ID should include all the atoms which this fix will
potentially select. I.e. the group-ID used in the voronoi compute should
include the same atoms as that indicated by the *types* keyword. If the
*diff* keyword is used, the group-ID should include atoms of all types
in the simulation.
The keyword *ke* takes *yes* (default) or *no* as its value. It two
atoms are swapped with different masses, then a value of *yes* will
rescale their respective velocities to conserve the kinetic energy of
the system. A value of *no* will perform no rescaling, so that
kinetic energy is not conserved. See the restriction on this keyword
below.
The *region* keyword takes a *region-ID* as its value. If specified,
then only atoms *I* and *J* within the geometric region will be
considered as swap partners. See the :doc:`region <region>` command
for details. This means the group-ID for the :doc:`compute
voronoi/atom <compute_voronoi_atom>` command also need only contain
atoms within the region.
The keyword *rates* can modify the swap rate based on the type of atom
*J*. Ntype values must be specified, where Ntype = the number of atom
types in the system. Each value is used to scale the probability
weighting given by the equation above. In the third example command
above, a simulation has 3 atoms types. Atom *I*s of type 1 are
eligible for swapping. Swaps may occur with atom *J*s of all 3 types.
Assuming all *J* atoms are equidistant from an atom *I*, *J* atoms of
type 1 will be 3x more likely to be selected as a swap partner than
atoms of type 2. And *J* atoms of type 3 will be 6.5x more likely to
be selected than atoms of type 2. If the *rates* keyword is not used,
all atom types will be treated with the same probability during selection
of swap attempts.
Restart, fix_modify, output, run start/stop, minimize info
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This fix writes the state of the fix to :doc:`binary restart files
<restart>`. This includes information about the random number generator
seed, the next timestep for MC exchanges, and the number of exchange
attempts and successes. See the :doc:`read_restart <read_restart>`
command for info on how to re-specify a fix in an input script that
reads a restart file, so that the operation of the fix continues in an
uninterrupted fashion.
None of the :doc:`fix_modify <fix_modify>` options are relevant to this
fix.
This fix computes a global vector of length 2, which can be accessed
by various :doc:`output commands <Howto_output>`. The vector values are
the following global cumulative quantities:
#. swap attempts
#. swap accepts
The vector values calculated by this fix are "intensive".
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is not invoked during
:doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
This fix is part of the MC package. It is only enabled if LAMMPS was
built with that package. See the :doc:`Build package <Build_package>`
doc page for more info. Also this fix requires that the :ref:`VORONOI
package <PKG-VORONOI>` is installed, otherwise the fix will not be
compiled.
The :doc:`compute voronoi/atom <compute_voronoi_atom>` command
referenced by the required voro-ID must return neighboring atoms as
illustrated in the examples above.
If this fix is used with systems that do not have per-type masses
(e.g. atom style sphere), the *ke* keyword must be set to *off* since
the implemented algorithm will not be able to re-scale velocities
properly.
Related commands
""""""""""""""""
:doc:`fix nvt <fix_nh>`, :doc:`compute voronoi/atom <compute_voronoi_atom>`
:doc:`delete_atoms <delete_atoms>`, :doc:`fix gcmc <fix_gcmc>`,
:doc:`fix atom/swap <fix_atom_swap>`, :doc:`fix mol/swap <fix_mol_swap>`,
:doc:`fix sgcmc <fix_sgcmc>`
Default
"""""""
The option defaults are *ke* = yes and *rates* = 1 for all atom types.
----------
.. _TavennerMDkMC:
**(Tavenner 2023)** J Tavenner, M Mendelev, J Lawson, Computational
Materials Science, 218, 111929 (2023).

31
doc/src/fix_qeq_rel_reaxff.rst
View File

@ -37,18 +37,18 @@ Examples
Description
"""""""""""
.. versionadded:: 19Nov2024
.. versionadded:: 2Apr2025
This fix implements the QEqR method for charge equilibration, which
differs from the QEq charge equilibration method :ref:`(Rappe and
Goddard) <Rappe4>` only in how external electric fields are accounted
for. This fix therefore raises a warning when used without :doc:`fix
efield <fix_efield>` since :doc:`fix qeq/reaxff <fix_qeq_reaxff>` should
be used without an external electric field. Charges are computed with
the QEqR method by minimizing the electrostatic energy of the system in
the same way as the QEq method but where the absolute electronegativity,
:math:`\chi_i`, of each atom in the QEq method is replaced with an
effective electronegativity given by
This fix implements the QEqR method :ref:`(Lalli) <lalli2>` for charge
equilibration, which differs from the QEq charge equilibration method
:ref:`(Rappe and Goddard) <Rappe4>` only in how external electric fields
are accounted for. This fix therefore raises a warning when used without
:doc:`fix efield <fix_efield>` since :doc:`fix qeq/reaxff <fix_qeq_reaxff>`
should be used when no external electric field is present. Charges are
computed with the QEqR method by minimizing the electrostatic energy of
the system in the same way as the QEq method but where the absolute
electronegativity, :math:`\chi_i`, of each atom in the QEq method is
replaced with an effective electronegativity given by
.. math::
\chi_{\mathrm{r}i} = \chi_i + \frac{\sum_{j=1}^{N} \beta(\phi_i - \phi_j) S_{ij}}
@ -61,8 +61,9 @@ external electric field and :math:`S_{ij}` is the overlap integral
between atoms :math:`i` and :math:`j`. This formulation is advantageous
over the method used by :doc:`fix qeq/reaxff <fix_qeq_reaxff>` to
account for an external electric field in that it permits periodic
boundaries in the direction of an external electric field and in that it
does not worsen long-range charge transfer seen with QEq.
boundaries in the direction of an external electric field and in
that it does not worsen long-range charge transfer seen with
QEq. See :ref:`Lalli <lalli2>` for further details.
This fix is typically used in conjunction with the ReaxFF force field
model as implemented in the :doc:`pair_style reaxff <pair_reaxff>`
@ -184,6 +185,10 @@ scale = 1.0 and maxiter = 200
----------
.. _lalli2:
**(Lalli)** Lalli and Giusti, Journal of Chemical Physics, 162, 174311 (2025).
.. _Rappe4:
**(Rappe)** Rappe and Goddard III, Journal of Physical Chemistry, 95,

11
doc/src/fix_qtpie_reaxff.rst
View File

@ -59,8 +59,7 @@ and atom :math:`j`.
The effect of an external electric field can be incorporated into the QTPIE
method by modifying the absolute or effective electronegativities of each
atom :ref:`(Chen) <qtpie-Chen>`. This fix models the effect of an external
electric field by using the effective electronegativity given in
:ref:`(Gergs) <Gergs>`:
electric field by using the effective electronegativity :ref:`(Lalli) <lalli>`
.. math::
\tilde{\chi}_{\mathrm{r}i} = \frac{\sum_{j=1}^{N} (\chi_i - \chi_j + \beta(\phi_i - \phi_j)) S_{ij}}
@ -68,7 +67,8 @@ electric field by using the effective electronegativity given in
where :math:`\beta` is a scaling factor and :math:`\phi_i` and :math:`\phi_j`
are the electric potentials at the positions of atoms :math:`i` and :math:`j`
due to the external electric field.
due to the external electric field. Additional details regarding the
implementation and performance of this fix are provided in :ref:`Lalli <lalli>`.
This fix is typically used in conjunction with the ReaxFF force
field model as implemented in the :doc:`pair_style reaxff <pair_reaxff>`
@ -206,10 +206,9 @@ scale = 1.0 and maxiter = 200
**(Chen)** Chen, Jiahao. Theory and applications of fluctuating-charge models.
University of Illinois at Urbana-Champaign, 2009.
.. _Gergs:
.. _lalli:
**(Gergs)** Gergs, Dirkmann and Mussenbrock.
Journal of Applied Physics 123.24 (2018).
**(Lalli)** Lalli and Giusti, Journal of Chemical Physics, 162, 174311 (2025).
.. _qeq-Aktulga2:

43
doc/src/fix_sgcmc.rst
View File

@ -30,7 +30,9 @@ Syntax
N = number of times sampling window is moved during one MC cycle
*window_size* frac
frac = size of sampling window (must be between 0.5 and 1.0)
*atomic/energy* yes/no
yes = use the atomic energy method to calculate energy changes
no = use the default method to calculate energy changes
Examples
""""""""
@ -127,6 +129,14 @@ The number of times the window is moved during a MC cycle is set using
the parameter *window_moves* (see Sect. III.B in :ref:`Sadigh1
<Sadigh1>` for details).
The *atomic/energy* keyword controls which method is used for calculating
the energy change when atom types are swapped. A value of *no*
uses the default method, see discussion below in Restrictions section.
A value of *yes* uses the atomic energy method,
if the method has been implemented for the LAMMPS energy model,
otherwise LAMMPS will exit with an error message.
So far this has only been implemented for EAM type potentials.
------------
Restart, fix_modify, output, run start/stop, minimize info
@ -159,16 +169,26 @@ page for more info.
This fix style requires an :doc:`atom style <atom_style>` with per atom
type masses.
At present the fix provides optimized subroutines for EAM type
potentials (see above) that calculate potential energy changes due to
*local* atom type swaps very efficiently. Other potentials are
supported by using the generic potential functions. This, however, will
lead to exceedingly slow simulations since it implies that the
energy of the *entire* system is recomputed at each MC trial step. If
other potentials are to be used it is strongly recommended to modify and
optimize the existing generic potential functions for this purpose.
Also, the generic energy calculation can not be used for parallel
execution i.e. it only works with a single MPI process.
The fix provides three methods for calculating the potential energy
change due to atom type swaps. For EAM type potentials, the default
method is a carefully optimized local energy change calculation that
is part of the source code for this fix. It takes advantage of the
specific computational and communication requirements of EAM. Customizing
the local method to handle other energy models such as Tersoff has been done,
but these cases are not supported in the public LAMMPS code.
For all other LAMMPS energy models, the default method calculates
the *total* potential energy of the system before and after each
atom type swap. This method does not depend on the details of the
energy model and so is guaranteed to be correct. It is also
orders of magnitude slower than the custom EAM calculation.
In addition, it can not be used with parallel execution i.e. only
a single MPI process is allowed.
The third method uses the *atomic/energy* keyword described above.
This allows parallel execution and it is also a local calculation,
making it only a bit slower than a fully-optimized local calculation.
So far, this has been implemented for EAM type potentials.
It is straightforward to extend this to other potentials,
requiring adding an atomic energy method to the pair style.
------------
@ -180,6 +200,7 @@ The optional parameters default to the following values:
* *randseed* = 324234
* *window_moves* = 8
* *window_size* = automatic
* *atomic/energy* = no
------------

64
doc/src/pair_granular.rst
View File

@ -44,7 +44,7 @@ Examples
pair_coeff * * hertz 1000.0 50.0 tangential mindlin 1000.0 1.0 0.4 heat area 0.1
pair_style granular
pair_coeff * * mdr 5e6 0.4 1.9e5 2.0 0.5 0.5 tangential linear_history 940.0 0.0 0.7 rolling sds 2.7e5 0.0 0.6 damping none
pair_coeff * * mdr 5e6 0.4 1.9e5 2.0 0.5 0.5 tangential linear_history 940.0 1.0 0.7 rolling sds 2.7e5 0.0 0.6 damping mdr 1
Description
"""""""""""
@ -88,7 +88,8 @@ and their required arguments are:
3. *hertz/material* : E, :math:`\eta_{n0}` (or :math:`e`), :math:`\nu`
4. *dmt* : E, :math:`\eta_{n0}` (or :math:`e`), :math:`\nu`, :math:`\gamma`
5. *jkr* : E, :math:`\eta_{n0}` (or :math:`e`), :math:`\nu`, :math:`\gamma`
6. *mdr* : :math:`E`, :math:`\nu`, :math:`Y`, :math:`\Delta\gamma`, :math:`\psi_b`, :math:`e`
6. *mdr* : :math:`E`, :math:`\nu`, :math:`Y`, :math:`\Delta\gamma`,
:math:`\psi_b`, :math:`\eta_{n0}`
Here, :math:`k_n` is spring stiffness (with units that depend on model
choice, see below); :math:`\eta_{n0}` is a damping prefactor (or, in its
@ -177,6 +178,8 @@ two-part series :ref:`Zunker and Kamrin Part I <Zunker2024I>` and
:ref:`Zunker and Kamrin Part II <Zunker2024II>`. Further development
and demonstrations of its application to industrially relevant powder
compaction processes are presented in :ref:`Zunker et al. <Zunker2025>`.
If you use the *mdr* normal model the only supported damping option is
the *mdr* damping class described below.
The model requires the following inputs:
@ -200,8 +203,8 @@ The model requires the following inputs:
triggered. Lower values of :math:`\psi_b` delay the onset of the bulk elastic
response.
6. *Coefficient of restitution* :math:`0 \le e \le 1` : The coefficient of
restitution is a tunable parameter that controls damping in the normal direction.
6. *Damping coefficent* :math:`\eta_{n0} \ge 0` : The damping coefficient
is a tunable parameter that controls damping in the normal direction.
.. note::
@ -213,18 +216,12 @@ The *mdr* model produces a nonlinear force-displacement response, therefore the
critical timestep :math:`\Delta t` depends on the inputs and level of
deformation. As a conservative starting point the timestep can be assumed to be
dictated by the bulk elastic response such that
:math:`\Delta t = 0.35\sqrt{m/k_\textrm{bulk}}`, where :math:`m` is the mass of
:math:`\Delta t = 0.08\sqrt{m/k_\textrm{bulk}}`, where :math:`m` is the mass of
the smallest particle and :math:`k_\textrm{bulk} = \kappa R_\textrm{min}` is an
effective stiffness related to the bulk elastic response.
Here, :math:`\kappa = E/(3(1-2\nu))` is the bulk modulus and
:math:`R_\textrm{min}` is the radius of the smallest particle.
.. note::
The *mdr* model requires some specific settings to function properly,
please read the following text carefully to ensure all requirements are
followed.
The *atom_style* must be set to *sphere 1* to enable dynamic particle
radii. The *mdr* model is designed to respect the incompressibility of
plastic deformation and inherently tracks free surface displacements
@ -253,13 +250,6 @@ algorithm see :ref:`Zunker et al. <Zunker2025>`.
newton off
The damping model must be set to *none*. The *mdr* model already has a built
in damping model.
.. code-block:: LAMMPS
pair_coeff * * mdr 5e6 0.4 1.9e5 2 0.5 0.5 damping none
The definition of multiple *mdr* models in the *pair_style* is currently not
supported. Similarly, the *mdr* model cannot be combined with a different normal
model in the *pair_style*. Physically this means that only one homogeneous
@ -270,7 +260,7 @@ The *mdr* model currently only supports *fix wall/gran/region*, not
any *fix wall/gran/region* commands must also use the *mdr* model.
Additionally, the following *mdr* inputs must match between the
*pair_style* and *fix wall/gran/region* definitions: :math:`E`,
:math:`\nu`, :math:`Y`, :math:`\psi_b`, and :math:`e`. The exception
:math:`\nu`, :math:`Y`, :math:`\psi_b`, and :math:`\eta_{n0}`. The exception
is :math:`\Delta\gamma`, which may vary, permitting different
adhesive behaviors between particle-particle and particle-wall interactions.
@ -336,6 +326,7 @@ for the damping model currently supported are:
3. *viscoelastic*
4. *tsuji*
5. *coeff_restitution*
6. *mdr* (class) : :math:`d_{type}`
If the *damping* keyword is not specified, the *viscoelastic* model is
used by default.
@ -425,6 +416,37 @@ the damping coefficient, it accurately reproduces the specified coefficient of
restitution for both monodisperse and polydisperse particle pairs. This damping
model is not compatible with cohesive normal models such as *JKR* or *DMT*.
The *mdr* damping class contains multiple damping models that can be toggled between
by specifying different integer values for the :math:`d_{type}` input parameter. This
damping option is only compatible with the normal *mdr* contact model.
Setting :math:`d_{type} = 1` is the suggested damping option. This specifies a damping
model that takes into account the contact stiffness :math:`k_{mdr}` calculated
by the normal *mdr* contact model to determine the damping coefficient:
.. math::
\eta_n = \eta_{n0} (m_{eff}k_{mdr})^{1/2},
where :math:`k_{mdr}` is proportional to contact radius :math:`a_{mdr}` tracked by the
normal *mdr* contact model:
.. math::
k_{mdr} = 2 E_{eff} a_{mdr}.
In this case, :math:`\eta_{n0}` is simply a dimensionless coefficient that scales the
the overall damping coefficient.
The other supported option is :math:`d_{type} = 2`, which defines a simple damping model
similar to the *velocity* option
.. math::
\eta_n = \eta_{n0},
but has additional checks to avoid non-physical damping after plastic deformation.
The total normal force is computed as the sum of the elastic and
damping components:
@ -1068,8 +1090,8 @@ a bulk elastic response. Journal of the Mechanics and Physics of Solids,
**(Zunker et al, 2025)** Zunker, W., Dunatunga, S., Thakur, S.,
Tang, P., & Kamrin, K. (2025). Experimentally validated DEM for large
deformation powder compaction: mechanically-derived contact model and
screening of non-physical contacts.
deformation powder compaction: Mechanically-derived contact model and
screening of non-physical contacts. Powder Technology, 120972.
.. _Luding2008:

23
doc/src/pair_lj_smooth.rst
View File

@ -48,13 +48,19 @@ At the inner cutoff the force and its first derivative
will match the non-smoothed LJ formula. At the outer cutoff the force
and its first derivative will be 0.0. The inner cutoff cannot be 0.0.
Explicit expressions for the coefficients C1, C2, C3, C4, as well as the
energy discontinuity at the cutoff can be found here :ref:`(Leoni_1) <Leoni_1>`
and here :ref:`(Leoni_2) <Leoni_2>`
.. note::
this force smoothing causes the energy to be discontinuous both
in its values and first derivative. This can lead to poor energy
conservation and may require the use of a thermostat. Plot the energy
and force resulting from this formula via the
:doc:`pair_write <pair_write>` command to see the effect.
conservation and may require the use of a thermostat. The energy
value discontinuity can be eliminated by shifting the potential
energy to be zero at the outer cutoff using the pair_modify shift
option. With or without shifting, you can plot the resulting energy
and force via the :doc:`pair_write <pair_write>` command to see the effect.
The following coefficients must be defined for each pair of atoms
types via the :doc:`pair_coeff <pair_coeff>` command as in the examples
@ -122,3 +128,14 @@ Default
"""""""
none
----------
.. _Leoni_1:
**(Leoni_1)** F. Leoni et al., Phys Rev Lett, 134, 128201 (2025).
.. _Leoni_2:
**(Leoni_2)** F. Leoni et al., Phys Rev Lett, 134, Supplementary Material (2025).