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

Merge branch 'develop' into bond-harmonic-restrain

Browse Source
This commit is contained in:
Axel Kohlmeyer
2023-03-13 19:43:16 -04:00
parent c369c0252f 2d8936068e
commit 3ba37fc954
131 changed files with 7424 additions and 3152 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
doc/src/Build_extras.rst
View File

@ -274,7 +274,7 @@ To enable GPU binning via CUDA performance primitives set the Makefile variable
most modern GPUs.
To support the CUDA multiprocessor server you can set the define
``-DCUDA_PROXY``. Please note that in this case you must **not** use
``-DCUDA_MPS_SUPPORT``. Please note that in this case you must **not** use
the CUDA performance primitives and thus set the variable ``CUDPP_OPT``
to empty.
@ -1964,10 +1964,10 @@ OPENMP package
Apple offers the `Xcode package and IDE
<https://developer.apple.com/xcode/>`_ for compiling software on
macOS, so you have likely installed it to compile LAMMPS. Their
compiler is based on `Clang <https://clang.llvm.org/>`, but while it
compiler is based on `Clang <https://clang.llvm.org/>`_, but while it
is capable of processing OpenMP directives, the necessary header
files and OpenMP runtime library are missing. The `R developers
<https://www.r-project.org/>` have figured out a way to build those
<https://www.r-project.org/>`_ have figured out a way to build those
in a compatible fashion. One can download them from
`https://mac.r-project.org/openmp/
<https://mac.r-project.org/openmp/>`_. Simply adding those files as

1
doc/src/Howto.rst
View File

@ -70,6 +70,7 @@ Force fields howto
Howto_amoeba
Howto_tip3p
Howto_tip4p
Howto_tip5p
Howto_spc
Packages howto

22
doc/src/Howto_mdi.rst
View File

@ -63,18 +63,18 @@ The package also provides a :doc:`mdi plugin <mdi>` command, which
enables LAMMPS to operate as an MDI driver and load an MDI engine as a
plugin library.
The package furthermore includes a `fix mdi/qm <fix_mdi_qm>` command, in
which LAMMPS operates as an MDI driver in conjunction with a quantum
mechanics code as an MDI engine. The post_force() method of the
``fix_mdi_qm.cpp`` file shows how a driver issues MDI commands to another
code. This command can be used to couple to an MDI engine, which is
either a stand-alone code or a plugin library.
The package furthermore includes a :doc:`fix mdi/qm <fix_mdi_qm>`
command, in which LAMMPS operates as an MDI driver in conjunction with a
quantum mechanics code as an MDI engine. The post_force() method of the
``fix_mdi_qm.cpp`` file shows how a driver issues MDI commands to
another code. This command can be used to couple to an MDI engine,
which is either a stand-alone code or a plugin library.
As explained in the `fix mdi/qm <fix_mdi_qm>` command documentation, it
can be used to perform *ab initio* MD simulations or energy
minimizations, or to evaluate the quantum energy and forces for a series
of independent systems. The ``examples/mdi`` directory has example
input scripts for all of these use cases.
As explained in the :doc:`fix mdi/qm <fix_mdi_qm>` command
documentation, it can be used to perform *ab initio* MD simulations or
energy minimizations, or to evaluate the quantum energy and forces for a
series of independent systems. The ``examples/mdi`` directory has
example input scripts for all of these use cases.
----------

114
doc/src/Howto_spc.rst
View File

@ -20,7 +20,6 @@ atoms and the water molecule to run a rigid SPC model.
| LJ :math:`\epsilon`, :math:`\sigma` of OH, HH = 0.0
| :math:`r_0` of OH bond = 1.0
| :math:`\theta_0` of HOH angle = 109.47\ :math:`^{\circ}`
|
Note that as originally proposed, the SPC model was run with a 9
Angstrom cutoff for both LJ and Coulomb terms. It can also be used
@ -33,16 +32,123 @@ the partial charge assignments change:
| O charge = -0.8476
| H charge = 0.4238
|
See the :ref:`(Berendsen) <howto-Berendsen>` reference for more details on both
the SPC and SPC/E models.
Below is the code for a LAMMPS input file and a molecule file
(``spce.mol``) of SPC/E water for use with the :doc:`molecule command
<molecule>` demonstrating how to set up a small bulk water system for
SPC/E with rigid bonds.
.. 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/coul/cut 10.0
pair_coeff 1 1 0.1553 3.166
pair_coeff 1 2 0.0 1.0
pair_coeff 2 2 0.0 1.0
bond_style zero
bond_coeff 1 1.0
angle_style zero
angle_coeff 1 109.47
molecule water spce.mol
create_atoms 0 random 33 34564 NULL mol water 25367 overlap 1.33
timestep 1.0
fix rigid all shake 0.0001 10 10000 b 1 a 1
minimize 0.0 0.0 1000 10000
run 0 post no
reset_timestep 0
velocity all create 300.0 5463576
fix integrate all nvt temp 300.0 300.0 1.0
thermo_style custom step temp press etotal density pe ke
thermo 1000
run 20000 upto
write_data tip4p.data nocoeff
.. _spce_molecule:
.. code-block::
# Water molecule. SPC/E geometry
3 atoms
2 bonds
1 angles
Coords
1 0.00000 -0.06461 0.00000
2 0.81649 0.51275 0.00000
3 -0.81649 0.51275 0.00000
Types
1 1 # O
2 2 # H
3 2 # H
Charges
1 -0.8476
2 0.4238
3 0.4238
Bonds
1 1 1 2
2 1 1 3
Angles
1 1 2 1 3
Shake Flags
1 1
2 1
3 1
Shake Atoms
1 1 2 3
2 1 2 3
3 1 2 3
Shake Bond Types
1 1 1 1
2 1 1 1
3 1 1 1
Special Bond Counts
1 2 0 0
2 1 1 0
3 1 1 0
Special Bonds
1 2 3
2 1 3
3 1 2
Wikipedia also has a nice article on `water models <https://en.wikipedia.org/wiki/Water_model>`_.
----------
.. _howto-Berendsen:
**(Berendsen)** Berendsen, Grigera, Straatsma, J Phys Chem, 91,
6269-6271 (1987).
**(Berendsen)** Berendsen, Grigera, Straatsma, J Phys Chem, 91, 6269-6271 (1987).

244
doc/src/Howto_tip3p.rst
View File

@ -1,53 +1,211 @@
TIP3P water model
=================
The TIP3P water model as implemented in CHARMM
:ref:`(MacKerell) <howto-tip3p>` specifies a 3-site rigid water molecule with
charges and Lennard-Jones parameters assigned to each of the 3 atoms.
In LAMMPS the :doc:`fix shake <fix_shake>` command can be used to hold
the two O-H bonds and the H-O-H angle rigid. A bond style of
*harmonic* and an angle style of *harmonic* or *charmm* should also be
used.
The TIP3P water model as implemented in CHARMM :ref:`(MacKerell)
<howto-tip3p>` specifies a 3-site rigid water molecule with charges and
Lennard-Jones parameters assigned to each of the 3 atoms.
These are the additional parameters (in real units) to set for O and H
atoms and the water molecule to run a rigid TIP3P-CHARMM model with a
cutoff. The K values can be used if a flexible TIP3P model (without
fix shake) is desired. If the LJ epsilon and sigma for HH and OH are
set to 0.0, it corresponds to the original 1983 TIP3P model
:ref:`(Jorgensen) <Jorgensen1>`.
A suitable pair style with cutoff Coulomb would be:
| O mass = 15.9994
| H mass = 1.008
| O charge = -0.834
| H charge = 0.417
| LJ :math:`\epsilon` of OO = 0.1521
| LJ :math:`\sigma` of OO = 3.1507
| LJ :math:`\epsilon` of HH = 0.0460
| LJ :math:`\sigma` of HH = 0.4000
| LJ :math:`\epsilon` of OH = 0.0836
| LJ :math:`\sigma` of OH = 1.7753
| K of OH bond = 450
| :math:`r_0` of OH bond = 0.9572
| K of HOH angle = 55
| :math:`\theta` of HOH angle = 104.52\ :math:`^{\circ}`
|
* :doc:`pair_style lj/cut/coul/cut <pair_lj_cut_coul>`
These are the parameters to use for TIP3P with a long-range Coulomb
solver (e.g. Ewald or PPPM in LAMMPS), see :ref:`(Price) <Price1>` for
details:
or these commands for a long-range Coulomb model:
* :doc:`pair_style lj/cut/coul/long <pair_lj_cut_coul>`
* :doc:`pair_style lj/cut/coul/long/soft <pair_fep_soft>`
* :doc:`kspace_style pppm <kspace_style>`
* :doc:`kspace_style pppm/disp <kspace_style>`
In LAMMPS the :doc:`fix shake or fix rattle <fix_shake>` command can be
used to hold the two O-H bonds and the H-O-H angle rigid. 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.
The table below lists the force field parameters (in real :doc:`units
<units>`) to for the water molecule atoms to run a rigid or flexible
TIP3P-CHARMM model with a cutoff, the original 1983 TIP3P model
:ref:`(Jorgensen) <Jorgensen1>`, or a TIP3P model with parameters
optimized for a long-range Coulomb solver (e.g. Ewald or PPPM in LAMMPS)
:ref:`(Price) <Price1>`. The K values can be used if a flexible TIP3P
model (without fix shake) is desired, for rigid bonds/angles they are
ignored.
.. list-table::
:header-rows: 1
:widths: auto
* - Parameter
- TIP3P-CHARMM
- TIP3P (original)
- TIP3P (Ewald)
* - O mass (amu)
- 15.9994
- 15.9994
- 15.9994
* - H mass (amu)
- 1.008
- 1.008
- 1.008
* - O charge (:math:`e`)
- -0.834
- -0.834
- -0.834
* - H charge (:math:`e`)
- 0.417
- 0.417
- 0.417
* - LJ :math:`\epsilon` of OO (kcal/mole)
- 0.1521
- 0.1521
- 0.1020
* - LJ :math:`\sigma` of OO (:math:`\AA`)
- 3.1507
- 3.1507
- 3.188
* - LJ :math:`\epsilon` of HH (kcal/mole)
- 0.0460
- 0.0
- 0.0
* - LJ :math:`\sigma` of HH (:math:`\AA`)
- 0.4
- 1.0
- 1.0
* - LJ :math:`\epsilon` of OH (kcal/mole)
- 0.0836
- 0.0
- 0.0
* - LJ :math:`\sigma` of OH (:math:`\AA`)
- 1.7753
- 1.0
- 1.0
* - K of OH bond (kcal/mole/:math:`\AA^2`)
- 450
- 450
- 450
* - :math:`r_0` of OH bond (:math:`\AA`)
- 0.9572
- 0.9572
- 0.9572
* - K of HOH angle (kcal/mole)
- 55.0
- 55.0
- 55.0
* - :math:`\theta_0` of HOH angle
- 104.52\ :math:`^{\circ}`
- 104.52\ :math:`^{\circ}`
- 104.52\ :math:`^{\circ}`
Below is the code for a LAMMPS input file and a molecule file
(``tip3p.mol``) of TIP3P water for use with the :doc:`molecule command
<molecule>` demonstrating how to set up a small bulk water system for
TIP3P with rigid bonds.
.. 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/coul/cut 8.0
pair_coeff 1 1 0.1521 3.1507
pair_coeff 2 2 0.0 1.0
bond_style zero
bond_coeff 1 0.9574
angle_style zero
angle_coeff 1 104.52
molecule water tip3p.mol
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
run 0 post no
reset_timestep 0
velocity all create 300.0 5463576
fix integrate all nvt temp 300 300 1.0
thermo_style custom step temp press etotal pe
thermo 1000
run 20000
write_data tip3p.data nocoeff
.. _tip3p_molecule:
.. code-block::
# Water molecule. TIP3P geometry
3 atoms
2 bonds
1 angles
Coords
1 0.00000 -0.06556 0.00000
2 0.75695 0.52032 0.00000
3 -0.75695 0.52032 0.00000
Types
1 1 # O
2 2 # H
3 2 # H
Charges
1 -0.834
2 0.417
3 0.417
Bonds
1 1 1 2
2 1 1 3
Angles
1 1 2 1 3
Shake Flags
1 1
2 1
3 1
Shake Atoms
1 1 2 3
2 1 2 3
3 1 2 3
Shake Bond Types
1 1 1 1
2 1 1 1
3 1 1 1
Special Bond Counts
1 2 0 0
2 1 1 0
3 1 1 0
Special Bonds
1 2 3
2 1 3
3 1 2
| O mass = 15.9994
| H mass = 1.008
| O charge = -0.830
| H charge = 0.415
| LJ :math:`\epsilon` of OO = 0.102
| LJ :math:`\sigma` of OO = 3.188
| LJ :math:`\epsilon`, :math:`\sigma` of OH, HH = 0.0
| K of OH bond = 450
| :math:`r_0` of OH bond = 0.9572
| K of HOH angle = 55
| :math:`\theta` of HOH angle = 104.52\ :math:`^{\circ}`
|
Wikipedia also has a nice article on `water models <https://en.wikipedia.org/wiki/Water_model>`_.

325
doc/src/Howto_tip4p.rst
View File

@ -2,100 +2,138 @@ TIP4P water model
=================
The four-point TIP4P rigid water model extends the traditional
three-point TIP3P model by adding an additional site, usually
massless, where the charge associated with the oxygen atom is placed.
This site M is located at a fixed distance away from the oxygen along
the bisector of the HOH bond angle. A bond style of *harmonic* and an
angle style of *harmonic* or *charmm* should also be used.
:doc:`three-point TIP3P <Howto_tip3p>` model by adding an additional
site M, usually massless, where the charge associated with the oxygen
atom is placed. This site M is located at a fixed distance away from
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.
A TIP4P model is run with LAMMPS using either these commands
for a cutoff model:
There are two ways to implement TIP4P water in LAMMPS:
* :doc:`pair_style tip4p/cut <pair_lj_cut_tip4p>`
* :doc:`pair_style lj/cut/tip4p/cut <pair_lj_cut_tip4p>`
#. Use a specially written pair style that uses the :ref:`TIP3P geometry
<tip3p_molecule>` without the point M. The point M location is then
implicitly derived from the other atoms or each water molecule and
used during the force computation. The forces on M are then
projected on the oxygen and the two hydrogen atoms. This is
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.
or these commands for a long-range model:
This can be done with the following pair styles for Coulomb with a cutoff:
* :doc:`pair_style tip4p/long <pair_coul>`
* :doc:`pair_style lj/cut/tip4p/long <pair_lj_cut_tip4p>`
* :doc:`pair_style lj/long/tip4p/long <pair_lj_long>`
* :doc:`pair_style tip4p/long/soft <pair_fep_soft>`
* :doc:`pair_style lj/cut/tip4p/long/soft <pair_fep_soft>`
* :doc:`kspace_style pppm/tip4p <kspace_style>`
* :doc:`kspace_style pppm/disp/tip4p <kspace_style>`
* :doc:`pair_style tip4p/cut <pair_lj_cut_tip4p>`
* :doc:`pair_style lj/cut/tip4p/cut <pair_lj_cut_tip4p>`
The bond lengths and bond angles should be held fixed using the
:doc:`fix shake <fix_shake>` or :doc:`fix rattle <fix_shake>` command,
unless a parameterization for a flexible TIP4P model is used. The
parameter sets listed below are all for rigid TIP4P model variants and
thus the bond and angle force constants are not used and can be set to
any legal value; only equilibrium length and angle are used.
or these commands for a long-range Coulomb treatment:
These are the additional parameters (in real units) to set for O and H
atoms and the water molecule to run a rigid TIP4P model with a cutoff
:ref:`(Jorgensen) <Jorgensen5>`. Note that the OM distance is specified in
the :doc:`pair_style <pair_style>` command, not as part of the pair
coefficients.
* :doc:`pair_style tip4p/long <pair_coul>`
* :doc:`pair_style lj/cut/tip4p/long <pair_lj_cut_tip4p>`
* :doc:`pair_style lj/long/tip4p/long <pair_lj_long>`
* :doc:`pair_style tip4p/long/soft <pair_fep_soft>`
* :doc:`pair_style lj/cut/tip4p/long/soft <pair_fep_soft>`
* :doc:`kspace_style pppm/tip4p <kspace_style>`
* :doc:`kspace_style pppm/disp/tip4p <kspace_style>`
| O mass = 15.9994
| H mass = 1.008
| O charge = -1.040
| H charge = 0.520
| :math:`r_0` of OH bond = 0.9572
| :math:`\theta` of HOH angle = 104.52\ :math:`^{\circ}`
| OM distance = 0.15
| LJ :math:`\epsilon` of O-O = 0.1550
| LJ :math:`\sigma` of O-O = 3.1536
| LJ :math:`\epsilon`, :math:`\sigma` of OH, HH = 0.0
| Coulomb cutoff = 8.5
|
The bond lengths and bond angles should be held fixed using the
:doc:`fix shake <fix_shake>` or :doc:`fix rattle <fix_shake>` command,
unless a parameterization for a flexible TIP4P model is used. The
parameter sets listed below are all for rigid TIP4P model variants and
thus the bond and angle force constants are not used and can be set to
any legal value; only equilibrium length and angle are used.
For the TIP4/Ice model (J Chem Phys, 122, 234511 (2005);
https://doi.org/10.1063/1.1931662) these values can be used:
#. Use an :ref:`explicit 4 point TIP4P geometry <tip4p_molecule>` where
the oxygen atom carries no charge and the M point no Lennard-Jones
interactions. Since :doc:`fix shake <fix_shake>` or :doc:`fix rattle
<fix_shake>` may not be applied to this kind of geometry, :doc:`fix
rigid or fix rigid/small <fix_rigid>` or its thermostatted variants
are required to maintain a rigid geometry. This avoids some of the
issues with respect to analysis and non-tip4p styles, but it is a
more costly force computation (more atoms in the same volume and thus
more neighbors in the neighbor lists) and requires a much shorter
timestep for stable integration of the rigid body motion. Since no
bonds or angles are required, they do not need to be defined and atom
style charge would be sufficient for a bulk TIP4P water system. In
order to avoid that LAMMPS produces an error due to the massless M
site a tiny non-zero mass needs to be assigned.
| O mass = 15.9994
| H mass = 1.008
| O charge = -1.1794
| H charge = 0.5897
| :math:`r_0` of OH bond = 0.9572
| :math:`\theta` of HOH angle = 104.52\ :math:`^{\circ}`
| OM distance = 0.1577
| LJ :math:`\epsilon` of O-O = 0.21084
| LJ :math:`\sigma` of O-O = 3.1668
| LJ :math:`\epsilon`, :math:`\sigma` of OH, HH = 0.0
| Coulomb cutoff = 8.5
|
The table below lists the force field parameters (in real :doc:`units
<units>`) to for a selection of popular variants of the TIP4P model.
There is the rigid TIP4P model with a cutoff :ref:`(Jorgensen)
<Jorgensen5>`, the TIP4/Ice model :ref:`(Abascal1) <Abascal1>`, the
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.
For the TIP4P/2005 model (J Chem Phys, 123, 234505 (2005);
https://doi.org/10.1063/1.2121687), these values can be used:
.. list-table::
:header-rows: 1
:widths: auto
| O mass = 15.9994
| H mass = 1.008
| O charge = -1.1128
| H charge = 0.5564
| :math:`r_0` of OH bond = 0.9572
| :math:`\theta` of HOH angle = 104.52\ :math:`^{\circ}`
| OM distance = 0.1546
| LJ :math:`\epsilon` of O-O = 0.1852
| LJ :math:`\sigma` of O-O = 3.1589
| LJ :math:`\epsilon`, :math:`\sigma` of OH, HH = 0.0
| Coulomb cutoff = 8.5
|
These are the parameters to use for TIP4P with a long-range Coulombic
solver (e.g. Ewald or PPPM in LAMMPS):
| O mass = 15.9994
| H mass = 1.008
| O charge = -1.0484
| H charge = 0.5242
| :math:`r_0` of OH bond = 0.9572
| :math:`\theta` of HOH angle = 104.52\ :math:`^{\circ}`
| OM distance = 0.1250
| LJ :math:`\epsilon` of O-O = 0.16275
| LJ :math:`\sigma` of O-O = 3.16435
| LJ :math:`\epsilon`, :math:`\sigma` of OH, HH = 0.0
|
* - Parameter
- TIP4P (original)
- TIP4P/Ice
- TIP4P/2005
- TIP4P (Ewald)
* - O mass (amu)
- 15.9994
- 15.9994
- 15.9994
- 15.9994
* - H mass (amu)
- 1.008
- 1.008
- 1.008
- 1.008
* - O or M charge (:math:`e`)
- -1.040
- -1.1794
- -1.1128
- -1.04844
* - H charge (:math:`e`)
- 0.520
- 0.5897
- 0.5564
- 0.52422
* - LJ :math:`\epsilon` of OO (kcal/mole)
- 0.1550
- 0.1577
- 0.1852
- 0.16275
* - LJ :math:`\sigma` of OO (:math:`\AA`)
- 3.1536
- 3.1668
- 3.1589
- 3.16435
* - LJ :math:`\epsilon` of HH, MM, OH, OM, HM (kcal/mole)
- 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
* - :math:`r_0` of OH bond (:math:`\AA`)
- 0.9572
- 0.9572
- 0.9572
- 0.9572
* - :math:`\theta_0` of HOH angle
- 104.52\ :math:`^{\circ}`
- 104.52\ :math:`^{\circ}`
- 104.52\ :math:`^{\circ}`
- 104.52\ :math:`^{\circ}`
* - OM distance (:math:`\AA`)
- 0.15
- 0.1577
- 0.1546
- 0.1250
Note that the when using the TIP4P pair style, the neighbor list cutoff
for Coulomb interactions is effectively extended by a distance 2 \* (OM
@ -108,6 +146,117 @@ trade-off for your model. The OM distance and the LJ and Coulombic
cutoffs are set in the :doc:`pair_style lj/cut/tip4p/long
<pair_lj_cut_tip4p>` command.
Below is the code for a LAMMPS input file using the implicit method and
the :ref:`TIP3P molecule file <tip3p_molecule>`. Because the TIP4P
charges are different from TIP3P they need to be reset (or the molecule
file changed):
.. 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/cut 1 2 1 1 0.15 8.0
pair_coeff 1 1 0.1550 3.1536
pair_coeff 2 2 0.0 1.0
bond_style zero
bond_coeff 1 0.9574
angle_style zero
angle_coeff 1 104.52
molecule water tip3p.mol # this uses the TIP3P geometry
create_atoms 0 random 33 34564 NULL mol water 25367 overlap 1.33
# must change charges for TIP4P
set type 1 charge -1.040
set type 2 charge 0.520
fix rigid all shake 0.001 10 10000 b 1 a 1
minimize 0.0 0.0 1000 10000
run 0 post no
reset_timestep 0
velocity all create 300.0 5463576
fix integrate all nvt temp 300 300 1.0
thermo_style custom step temp press etotal pe
thermo 1000
run 20000
write_data tip3p.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/nvt/small
<fix_rigid>` no bonds need to be defined and thus no extra storage needs
to be reserved for them, but we need to switch to atom style full or use
:doc:`fix property/atom mol <fix_property_atom>` so that fix
rigid/nvt/small can identify rigid bodies by their molecule ID:
.. code-block:: LAMMPS
units real
atom_style charge
region box block -5 5 -5 5 -5 5
create_box 3 box
mass 1 15.9994
mass 2 1.008
mass 3 1.0e-100
pair_style lj/cut/coul/cut 8.0
pair_coeff 1 1 0.1550 3.1536
pair_coeff 2 2 0.0 1.0
pair_coeff 3 3 0.0 1.0
fix mol all property/atom mol
molecule water tip4p.mol
create_atoms 0 random 33 34564 NULL mol water 25367 overlap 1.33
timestep 0.1
fix integrate all rigid/nvt/small molecule temp 300.0 300.0 1.0
velocity all create 300.0 5463576
thermo_style custom step temp press etotal density pe ke
thermo 1000
run 20000
write_data tip4p.data nocoeff
.. _tip4p_molecule:
.. code-block::
# Water molecule. Explicit TIP4P geometry for use with fix rigid
4 atoms
Coords
1 0.00000 -0.06556 0.00000
2 0.75695 0.52032 0.00000
3 -0.75695 0.52032 0.00000
4 0.00000 0.08444 0.00000
Types
1 1 # O
2 2 # H
3 2 # H
4 3 # M
Charges
1 0.000
2 0.520
3 0.520
4 -1.040
Wikipedia also has a nice article on `water models <https://en.wikipedia.org/wiki/Water_model>`_.
----------
@ -116,3 +265,13 @@ Wikipedia also has a nice article on `water models <https://en.wikipedia.org/wik
**(Jorgensen)** Jorgensen, Chandrasekhar, Madura, Impey, Klein, J Chem
Phys, 79, 926 (1983).
.. _Abascal1:
**(Abascal1)** Abascal, Sanz, Fernandez, Vega, J Chem Phys, 122, 234511 (2005)
https://doi.org/10.1063/1.1931662
.. _Abascal2:
**(Abascal2)** Abascal, J Chem Phys, 123, 234505 (2005)
https://doi.org/10.1063/1.2121687

161
doc/src/Howto_tip5p.rst Normal file
View File

@ -0,0 +1,161 @@
TIP5P water model
=================
The five-point TIP5P rigid water model extends the :doc:`three-point
TIP3P model <Howto_tip3p>` by adding two additional sites L, usually
massless, where the charge associated with the oxygen atom is placed.
These sites L are located at a fixed distance away from the oxygen atom,
forming a tetrahedral angle that is rotated by 90 degrees from the HOH
plane. Those sites thus somewhat approximate lone pairs of the oxygen
and consequently improve the water structure to become even more
"tetrahedral" in comparison to the :doc:`four-point TIP4P model
<Howto_tip4p>`.
A suitable pair style with cutoff Coulomb would be:
* :doc:`pair_style lj/cut/coul/cut <pair_lj_cut_coul>`
or these commands for a long-range model:
* :doc:`pair_style lj/cut/coul/long <pair_lj_cut_coul>`
* :doc:`pair_style lj/cut/coul/long/soft <pair_fep_soft>`
* :doc:`kspace_style pppm <kspace_style>`
* :doc:`kspace_style pppm/disp <kspace_style>`
A TIP5P model *must* be run using a :doc:`rigid fix <fix_rigid>` since
there is no other option to keep this kind of structure rigid in LAMMPS.
In order to avoid that LAMMPS produces an error due to the massless L
sites, those need to be assigned a tiny non-zero mass.
The table below lists the force field parameters (in real :doc:`units
<units>`) to for a the TIP5P model with a cutoff :ref:`(Mahoney)
<Mahoney>` and the TIP5P-E model :ref:`(Rick) <Rick>` for use with a
long-range Coulombic solver (e.g. Ewald or PPPM in LAMMPS).
.. list-table::
:header-rows: 1
:widths: auto
* - Parameter
- TIP5P
- TIP5P-E
* - O mass (amu)
- 15.9994
- 15.9994
* - H mass (amu)
- 1.008
- 1.008
* - O charge (:math:`e`)
- 0.0
- 0.0
* - L charge (:math:`e`)
- -0.241
- -0.241
* - H charge (:math:`e`)
- 0.241
- 0.241
* - LJ :math:`\epsilon` of OO (kcal/mole)
- 0.1600
- 0.1780
* - LJ :math:`\sigma` of OO (:math:`\AA`)
- 3.1200
- 3.0970
* - LJ :math:`\epsilon` of HH, LL, OH, OL, HL (kcal/mole)
- 0.0
- 0.0
* - LJ :math:`\sigma` of HH, LL, OH, OL, HL (:math:`\AA`)
- 1.0
- 1.0
* - :math:`r_0` of OH bond (:math:`\AA`)
- 0.9572
- 0.9572
* - :math:`\theta_0` of HOH angle
- 104.52\ :math:`^{\circ}`
- 104.52\ :math:`^{\circ}`
* - OL distance (:math:`\AA`)
- 0.70
- 0.70
* - :math:`\theta_0` of LOL angle
- 109.47\ :math:`^{\circ}`
- 109.47\ :math:`^{\circ}`
Below is the code for a LAMMPS input file for setting up a simulation of
TIP5P water with a molecule file. Because of using :doc:`fix
rigid/nvt/small <fix_rigid>` no bonds need to be defined and thus no
extra storage needs to be reserved for them, but we need to switch to
atom style full or use :doc:`fix property/atom mol <fix_property_atom>`
so that fix rigid/nvt/small can identify rigid bodies by their molecule
ID:
.. code-block:: LAMMPS
units real
atom_style charge
region box block -5 5 -5 5 -5 5
create_box 3 box
mass 1 15.9994
mass 2 1.008
mass 3 1.0e-100
pair_style lj/cut/coul/cut 8.0
pair_coeff 1 1 0.160 3.12
pair_coeff 2 2 0.0 1.0
pair_coeff 3 3 0.0 1.0
fix mol all property/atom mol
molecule water tip5p.mol
create_atoms 0 random 33 34564 NULL mol water 25367 overlap 1.33
timestep 0.20
fix integrate all rigid/nvt/small molecule temp 300.0 300.0 1.0
reset_timestep 0
velocity all create 300.0 5463576
thermo_style custom step temp press etotal density pe ke
thermo 1000
run 20000
write_data tip5p.data nocoeff
.. _tip5p_molecule:
.. code-block::
# Water molecule. Explicit TIP5P geometry for use with fix rigid
5 atoms
Coords
1 0.00000 -0.06556 0.00000
2 0.75695 0.52032 0.00000
3 -0.75695 0.52032 0.00000
4 0.00000 -0.46971 0.57154
5 0.00000 -0.46971 -0.57154
Types
1 1 # O
2 2 # H
3 2 # H
4 3 # L
5 3 # L
Charges
1 0.000
2 0.241
3 0.241
4 -0.241
5 -0.241
Wikipedia also has a nice article on `water models <https://en.wikipedia.org/wiki/Water_model>`_.
----------
.. _Mahoney:
**(Mahoney)** Mahoney, Jorgensen, J Chem Phys 112, 8910 (2000)
.. _Rick:
**(Rick)** Rick, J Chem Phys 120, 6085 (2004)

2
doc/src/Install_linux.rst
View File

@ -42,7 +42,7 @@ static linkage, there is no ``liblammps.so`` library file and thus also the
LAMMPS python module, which depends on it, is not included.
The compressed tar archives available for download have names following
the pattern `lammps-linux-x86_64-<version>.tar.gz` and will all unpack
the pattern ``lammps-linux-x86_64-<version>.tar.gz`` and will all unpack
into a ``lammps-static`` folder. The executables are then in the
``lammps-static/bin/`` folder. Since they do not depend on any other
software, they may be freely moved or copied around.

11
doc/src/Install_windows.rst
View File

@ -17,11 +17,12 @@ install the Windows MPI package (MPICH2 from Argonne National Labs),
needed to run in parallel with MPI.
The LAMMPS binaries contain *all* :doc:`optional packages <Packages>`
included in the source distribution except: KIM, KOKKOS, MSCG, PYTHON,
ADIOS, H5MD, NETCDF, QMMM, ML-QUIP, and VTK.
The serial version also does not include the MPIIO and
LATBOLTZ packages. The GPU package is compiled for OpenCL with
mixed precision kernels.
included in the source distribution except: ADIOS, H5MD, KIM, ML-PACE,
ML-QUIP, MSCG, NETCDF, PLUMED, QMMM, SCAFACOS, and VTK. The serial
version also does not include the MPIIO and LATBOLTZ packages. The
PYTHON package is only available in the Python installers that bundle a
Python runtime. The GPU package is compiled for OpenCL with mixed
precision kernels.
The LAMMPS library is compiled as a shared library and the
:doc:`LAMMPS Python module <Python_module>` is installed, so that

2
doc/src/Intro_features.rst
View File

@ -195,7 +195,7 @@ Multi-replica models
* :doc:`parallel replica dynamics <prd>`
* :doc:`temperature accelerated dynamics <tad>`
* :doc:`parallel tempering <temper>`
* path-integral MD: `first variant <fix_pimd>`, `second variant <fix_ipi>`
* path-integral MD: :doc:`first variant <fix_pimd>`, :doc:`second variant <fix_ipi>`
* multi-walker collective variables with :doc:`Colvars <fix_colvars>` and :doc:`Plumed <fix_plumed>`
.. _prepost:

2
doc/src/balance.rst
View File

@ -498,7 +498,7 @@ data to other processors during load-balancing will be random or
deterministic. Random is generally faster; deterministic will ensure
the new ordering of atoms on each processor is the same each time the
same simulation is run. This can be useful for debugging purposes.
Since the balance commmand is a one-time operation, the default is
Since the balance command is a one-time operation, the default is
*yes* to perform sorting.
The *out* keyword writes a text file to the specified *filename* with

2
doc/src/bond_fene.rst
View File

@ -51,7 +51,7 @@ in the same form as in pair style :doc:`nm/cut <pair_nm>`. The bond energy is th
.. math::
E = -0.5 K r_0^2 \ln \left[ 1 - \left(\frac{r}{R_0}\right)^2\right] + \frac{E_0}{(n-m)} \left[ m \left(\frac{r_0}{r}\right)^n - n \left(\frac{r_0}{r}\right)^m \right]
E = -0.5 K R_0^2 \ln \left[ 1 - \left(\frac{r}{R_0}\right)^2\right] + \frac{E_0}{(n-m)} \left[ m \left(\frac{r_0}{r}\right)^n - n \left(\frac{r_0}{r}\right)^m \right]
Similar to the *fene* style, the generalized Lennard-Jones is cut off at
the potential minimum, :math:`r_0`, to be repulsive only. The following

22
doc/src/bond_table.rst
View File

@ -112,8 +112,9 @@ are estimated (less accurately) by the first two and last two force
values in the table.
The "EQ" parameter is also optional. If used, it is followed by a the
equilibrium bond length, which is used, for example, by the :doc:`fix shake <fix_shake>` command. If not used, the equilibrium bond
length is to the distance in the table with the lowest potential energy.
equilibrium bond length, which is used, for example, by the :doc:`fix
shake <fix_shake>` command. If not used, the equilibrium bond length is
to the distance in the table with the lowest potential energy.
Following a blank line, the next N lines list the tabulated values.
On each line, the first value is the index from 1 to N, the second value is
@ -135,16 +136,15 @@ one that matches the specified keyword.
----------
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
Restart info
""""""""""""
This bond style writes the settings for the "bond_style table"
command to :doc:`binary restart files <restart>`, so a bond_style
command does not need to specified in an input script that reads a
restart file. However, the coefficient information is not stored in
the restart file, since it is tabulated in the potential files. Thus,
bond_coeff commands do need to be specified in the restart input
script.
This bond style writes the settings for the "bond_style table" command
to :doc:`binary restart files <restart>`, so a bond_style command does
not need to specified in an input script that reads a restart file.
However, the coefficient information is not stored in the restart file,
since it is tabulated in the potential files. Thus, bond_coeff commands
do need to be specified in the restart input script.
Restrictions
""""""""""""

2
doc/src/bond_write.rst
View File

@ -70,7 +70,7 @@ be specified even if the potential has a finite value at r = 0.0.
Related commands
""""""""""""""""
:doc:`bond_style table <bond_table>`, `angle_write <angle_write>`,
:doc:`bond_style table <bond_table>`, :doc:`angle_write <angle_write>`,
:doc:`bond_style <bond_style>`, :doc:`bond_coeff <bond_coeff>`
Default

2
doc/src/dump.rst
View File

@ -375,7 +375,7 @@ output with each snapshot:
nx ny nz
The value dim will be 2 or 3 for 2d or 3d simulations. It is included
so that post-processing tools like `OVITO <https://www.ovito.org>`,
so that post-processing tools like `OVITO <https://www.ovito.org>`_,
which can visualize grid-based quantities know how to draw each grid
cell. The grid size will match the input script parameters for
grid(s) created by the computes or fixes which are referenced by the

16
doc/src/fix_adapt.rst
View File

@ -139,7 +139,7 @@ formulas for the meaning of these parameters:
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`buck/mdf <pair_mdf>` | a,c | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`coul/cut <pair_coul>` | scale | type pairs |
| :doc:`coul/cut, coul/cut/global <pair_coul>` | scale | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`coul/cut/soft <pair_fep_soft>` | lambda | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
@ -151,10 +151,16 @@ formulas for the meaning of these parameters:
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`coul/long/soft <pair_fep_soft>` | scale, lambda, coulombic_cutoff | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`coul/slater/long <pair_coul_slater>` | scale | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`coul/streitz <pair_coul>` | scale | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`eam, eam/alloy, eam/fs <pair_eam>` | scale | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`gauss <pair_gauss>` | a | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`harmonic/cut <pair_harmonic_cut>` | k, cutoff | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`lennard/mdf <pair_mdf>` | A,B | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`lj/class2 <pair_class2>` | epsilon,sigma | type pairs |
@ -181,6 +187,8 @@ formulas for the meaning of these parameters:
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`lubricate <pair_lubricate>` | mu | global |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`meam <pair_meam>` | scale | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`mie/cut <pair_mie>` | epsilon,sigma,gamma_repulsive,gamma_attractive | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`morse, morse/smooth/linear <pair_morse>` | D0,R0,alpha | type pairs |
@ -191,7 +199,7 @@ formulas for the meaning of these parameters:
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`nm/cut/coul/cut, nm/cut/coul/long <pair_nm>` | E0,R0,m,n,coulombic_cutoff | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`reaxff <pair_reaxff>` | chi, eta, gamma | type global |
| :doc:`pace, pace/extrapolation <pair_pace>` | scale | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`snap <pair_snap>` | scale | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
@ -203,11 +211,13 @@ formulas for the meaning of these parameters:
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`spin/neel <pair_spin_neel>` | coulombic_cutoff | type global |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`soft <pair_soft>` | a | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`table <pair_table>` | table_cutoff | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`ufm <pair_ufm>` | epsilon,sigma | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`soft <pair_soft>` | a | type pairs |
| :doc:`wf/cut <pair_wf_cut>` | epsilon,sigma,nu,mu | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
.. note::

144
doc/src/fix_alchemy.rst
View File

@ -12,7 +12,7 @@ Syntax
* ID, group-ID are documented in :doc:`fix <fix>` command
* alchemy = style name of this fix command
* v_name = variable with name that determines the :math:`\lambda_p` value
* v_name = variable with name that determines the :math:`\lambda_R` value
Examples
""""""""
@ -26,87 +26,126 @@ Description
.. versionadded:: TBD
This fix command enables running an "alchemical transformation" MD
simulation between two topologies (i.e. the same number and positions of
atoms, but differences in atom parameters like type, charge, bonds,
angles and so on). For this a :ref:`multi-partition run <partition>` is
required with exactly two partitions. During the MD run, the fix will
will determine a factor, :math:`\lambda_p`, for each partition *p* that
will be taken from an equal style or equivalent :doc:`variable
<variable>`. Typically, this variable would be chose to linearly ramp
*down* from 1.0 to 0.0 for the *first* partition (*p=0*) and linearly
ramp *up* from 0.0 to 1.0 for the *second* partition (*p=1*). The
forces used for the propagation of the atoms will be the sum of the
forces of the two systems combined and scaled with their respective
:math:`\lambda_p` factor. This allows to perform transformations that
are not easily possible with :doc:`pair style hybrid/scaled
<pair_hybrid>`, :doc:`fix adapt <fix_adapt>` or :doc:`fix adapt/fep
<fix_adapt_fep>`.
This fix command enables an "alchemical transformation" to be performed
between two systems, whereby one system slowly transforms into the other
over the course of a molecular dynamics run. This is useful for
measuring thermodynamic differences between two different systems. It
also allows transformations that are not easily possible with the
:doc:`pair style hybrid/scaled <pair_hybrid>`, :doc:`fix adapt
<fix_adapt>` or :doc:`fix adapt/fep <fix_adapt_fep>` commands.
.. note::
Example inputs are included in the ``examples/PACKAGES/alchemy``
directory for (a) transforming a pure copper system into a
copper/aluminum bronze alloy and (b) transforming two water molecules
in a box of water into a hydronium and a hydroxyl ion.
Since the definition of the variable to provide the :math:`\lambda_p` is
independent in the two partitions, no check is made that the two values
remain between 0.0 and 1.0 and that they add up to 1.0. So care needs to
be taken when defining those variables that this is the case.
The two systems must be defined as :doc:`separate replica
<Howto_replica>` and run in separate partitions of processors using the
:doc:`-partition <Run_options>` command-line switch. Exactly two
partitions must be specified, and each partition must use the same number
of processors and the same domain decomposition.
Due to the specifics of the implementation, the initial geometry and
dimensions of the system must be exactly the same and the fix will
synchronize them during the run. It is thus not possible to initialize
the two partitions by reading different data files or creating different
systems from scratch, but rather they have to be started from the same
system and then the desired modifications need to be applied to the
system of the second partition. The commands :doc:`pair style
hybrid/scaled <pair_hybrid>`, :doc:`fix adapt <fix_adapt>` or :doc:`fix
adapt/fep <fix_adapt_fep>` could be used for simulations where the
requirements for fix alchemy are not given.
Because the forces applied to the atoms are the same mix of the forces
from each partition and the simulation starts with the same atom
positions across both partitions, they will generate the same trajectory
of coordinates for each atom, and the same simulation box size and
shape. The latter two conditions are *enforced* by this fix; it
exchanges coordinates and box information between the replicas. This is
not strictly required, but since MD simulations are an example of a
chaotic system, even the tiniest random difference will eventually grow
exponentially into an unwanted divergence.
The commands below demonstrate how the setup for the second partition
can be done for the example of transforming a pure copper system into a
copper/aluminum bronze.
Otherwise, the properties of each atom (type, charge, bond and angle
partners, etc.), as well as energy and forces between interacting atoms
(pair, bond, angle styles, etc.) can be different in the two systems.
This can be initialized in the same input script by using commands which
only apply to one or the other replica. The example scripts use a
world-style :doc:`variable <variable>` command along with
:doc:`if/then/else <if>` commands for this purpose. The
:doc:`partition <partition>` command can also be used.
.. code-block:: LAMMPS
variable name world pure alloy
create_box 2 box
create_atoms 1 box
pair_style eam/alloy
pair_coeff * * AlCu.eam.alloy Cu Al
# replace 5% of copper with aluminum on the second partition only
variable name world pure alloy
if "${name} == alloy" then &
"set type 1 type/fraction 2 0.05 6745234"
# define ramp variable to combine the two different partitions
if "${name} == pure" then &
"variable ramp equal ramp(1.0,0.0)" &
else &
"variable ramp equal ramp(0.0,1.0)"
Both replicas must define an instance of this fix, but with a different
*v_name* variable. The named variable must be an equal-style or
equivalent :doc:`variable <variable>`. The two variables should be
defined so that one ramps *down* from 1.0 to 0.0 for the *first* replica
(*R=0*) and the other ramps *up* from 0.0 to 1.0 for the *second*
replica (*R=1*). A simple way is to do this is linearly, which can be
done using the ramp() function of the :doc:`variable <variable>`
command. You could also define a variable which returns a value between
0.0 and 1.0 as a non-linear function of the timestep. Here is a linear
example:
.. code-block:: LAMMPS
partition yes 1 variable ramp equal ramp(1.0,0.0)
partition yes 2 variable ramp equal ramp(0.0,1.0)
fix 2 all alchemy v_ramp
.. note::
The ``examples/PACKAGES/alchemy`` folder contains complete example
inputs for this command.
For an alchemical transformation, the two variables should sum to
exactly 1.0 at any timestep. LAMMPS does *NOT* check that this is
the case.
If you use the ``ramp()`` function to define the two variables, this fix
can easily be used across successive runs in the same input script by
ensuring each instance of the :doc:`run <run>` command specifies the
appropriate *start* or *stop* options.
At each timestep of an MD run, the two instances of this fix evaluate
their respective variables as a :math:`\lambda_R` factor, where *R* = 0
or 1 for each replica. The forces used by each system for the
propagation of their atoms is set to the sum of the forces for the two
systems, each scaled by their respective :math:`\lambda_R` factor. Thus,
during the MD run, the system will transform incrementally from the
first system to the second system.
.. note::
As mentioned above, the coordinates of the atoms and box size/shape
must be exactly the same in the two replicas. Therefore, it is
generally not a good idea to initialize the two replicas by reading
different data files or creating them individually from scratch.
Rather, a single system should be initialized and then desired
modifications applied to the system to either replica. If your
input script somehow induces the two systems to become different
(e.g. by performing :doc:`atom_modify sort <atom_modify>`
differently, or by adding or depositing a different number of atoms),
then LAMMPS will detect the mismatch and generate an error. This is
done by ensuring that each step the number and ordering of atoms is
identical within each pair of processors in the two replicas.
----------
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files <restart>`.
None of the :doc:`fix_modify <fix_modify>` options are relevant to this fix.
No information about this fix is written to :doc:`binary restart files
<restart>`. None of the :doc:`fix_modify <fix_modify>` options are
relevant to this fix.
This fix stores a global scalar (the current value of :math:`\lambda_p`)
This fix stores a global scalar (the current value of :math:`\lambda_R`)
and a global vector of length 3 which contains the potential energy of
the first partition, the second partition and the combined value,
respectively. The global scalar is unitless and "intensive", the vector
is in :doc:`energy units <units>` and "extensive". These values can be
used by any command that uses a global value from a fix as input. See
the :doc:`Howto output <Howto_output>` doc page for an overview of
LAMMPS output options.
the :doc:`output howto <Howto_output>` page for an overview of LAMMPS
output options.
This fix is not invoked during :doc:`energy minimization <minimize>`.
@ -117,12 +156,9 @@ This fix is part of the REPLICA package. It is only enabled if LAMMPS
was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
There may be only one instance of this fix in use at any time.
There may be only one instance of this fix in use at a time within
each replica.
This fix requires to perform a :ref:`multi-partition run <partition>`
with *exactly* two partitions.
This fix is *not* compatible with :doc:`load balancing <fix_balance>`.
Related commands
""""""""""""""""

2
doc/src/fix_balance.rst
View File

@ -314,7 +314,7 @@ data to other processors during load-balancing will be random or
deterministic. Random is generally faster; deterministic will ensure
the new ordering of atoms on each processor is the same each time the
same simulation is run. This can be useful for debugging purposes.
Since the fix balance commmand is performed during timestepping, the
Since the fix balance command is performed during timestepping, the
default is *no* so that sorting is not performed.
The *out* keyword writes text to the specified *filename* with the

6
doc/src/fix_deposit.rst
View File

@ -221,9 +221,9 @@ options choose a z-coordinate for insertion independently.
The vx, vy, and vz components of velocity for the inserted particle
are set by sampling a uniform distribution between the bounds set by
the values specified for the *vx*, *vy*, and *vz* keywords. Note that
normally, new particles should be a assigned a negative vertical
velocity so that they move towards the surface. For molecules, the
the values specified for the *vx*, *vy*, and *vz* keywords. Note that
normally, new particles should be a assigned a negative vertical
velocity so that they move towards the surface. For molecules, the
same velocity is given to every particle (no rotation or bond vibration).
If the *target* option is used, the velocity vector of the inserted

66
doc/src/fix_mscg.rst
View File

@ -40,15 +40,15 @@ Description
This fix applies the Multi-Scale Coarse-Graining (MSCG) method to
snapshots from a dump file to generate potentials for coarse-grained
simulations from all-atom simulations, using a force-matching
technique (:ref:`Izvekov <Izvekov>`, :ref:`Noid <Noid>`).
simulations from all-atom simulations, using a force-matching technique
(:ref:`Izvekov <Izvekov>`, :ref:`Noid <Noid>`).
It makes use of the MS-CG library, written and maintained by Greg
Voth's group at the University of Chicago, which is freely available
on their `MS-CG GitHub site <https://github.com/uchicago-voth/MSCG-release>`_. See instructions
on obtaining and installing the MS-CG library in the src/MSCG/README
file, which must be done before you build LAMMPS with this fix command
and use the command in a LAMMPS input script.
It makes use of the MS-CG library, written and maintained by Greg Voth's
group at the University of Chicago, which is freely available on their
`MS-CG GitHub site <https://github.com/uchicago-voth/MSCG-release>`_.
See instructions on obtaining and installing the MS-CG library in the
src/MSCG/README file, which must be done before you build LAMMPS with
this fix command and use the command in a LAMMPS input script.
An example script using this fix is provided the examples/mscg
directory.
@ -65,15 +65,18 @@ simulations is as follows:
6. Check the results of the force matching.
7. Run coarse-grained simulations using the new coarse-grained potentials.
This fix can perform the range finding and force matching steps 4 and
5 of the above workflow when used in conjunction with the
:doc:`rerun <rerun>` command. It does not perform steps 1-3 and 6-7.
This fix can perform the range finding and force matching steps 4 and 5
of the above workflow when used in conjunction with the :doc:`rerun
<rerun>` command. It does not perform steps 1-3 and 6-7.
Step 2 can be performed using a Python script (what is the name?)
provided with the MS-CG library which defines the coarse-grained model
and converts a standard LAMMPS dump file for an all-atom simulation
(step 1) into a LAMMPS dump file which has the positions of and forces
on the coarse-grained beads.
Step 2 can be performed using a Python script (cgmap), which defines the
coarse-grained model and converts a standard LAMMPS dump file for an
all-atom simulation (step 1) into a LAMMPS dump file which has the
positions of and forces on the coarse-grained beads. To use cgmap the
following repositories need to be downloaded and installed.
#. The custom lammpsdata branch of mdtraj from https://github.com/hockyg/mdtraj/tree/lammpsdata
#. The master branch of cgmap from https://github.com/uchicago-voth/cgmap
In step 3, an input file named "control.in" is needed by the MS-CG
library which sets parameters for the range finding and force matching
@ -83,12 +86,12 @@ info on this file.
When this fix is used to perform steps 4 and 5, the MS-CG library also
produces additional output files. The range finder functionality
(step 4) outputs files defining pair and bonded interaction ranges.
The force matching functionality (step 5) outputs tabulated force
files for every interaction in the system. Other diagnostic files can
also be output depending on the parameters in the MS-CG library input
script. Again, see the documentation provided with the MS-CG library
for more info.
(step 4) outputs files defining pair and bonded interaction ranges. The
force matching functionality (step 5) outputs tabulated force files for
every interaction in the system. Other diagnostic files can also be
output depending on the parameters in the MS-CG library input script.
Again, see the documentation provided with the MS-CG library for more
info.
----------
@ -97,8 +100,8 @@ be invoked. If *on*, the step 4 range finder functionality is invoked.
*off*, the step 5 force matching functionality is invoked.
If the *name* keyword is used, string names are defined to associate
with the integer atom types in LAMMPS. *Ntype* names must be
provided, one for each atom type (1-Ntype).
with the integer atom types in LAMMPS. *Ntype* names must be provided,
one for each atom type (1-Ntype).
The *max* keyword specifies the maximum number of bonds, angles, and
dihedrals a bead can have in the coarse-grained model.
@ -107,16 +110,13 @@ Restrictions
""""""""""""
This fix is part of the MSCG package. It is only enabled if LAMMPS was
built with that package. See the :doc:`Build package <Build_package>`
doc page for more info.
built with that package. Building the MSCG package also requires
external libraries. See the :doc:`Build_package` and :doc:`Build_extras`
pages for more info.
The MS-CG library uses C++11, which may not be supported by older
compilers. The MS-CG library also has some additional numeric library
dependencies, which are described in its documentation.
Currently, the MS-CG library is not setup to run in parallel with MPI,
so this fix can only be used in a serial LAMMPS build and run
on a single processor.
Currently, the MS-CG library is not set up to run in parallel with MPI,
so this fix can only be used in a serial LAMMPS build and run on a
single processor.
Related commands
""""""""""""""""

30
doc/src/fix_mvv_dpd.rst
View File

@ -65,33 +65,37 @@ a default value of 0.5 is used, which effectively reproduces the
standard velocity-Verlet (VV) scheme. For more details, see
:ref:`Groot <Groot2>`.
Fix *mvv/dpd* updates the position and velocity of each atom. It can
be used with the :doc:`pair_style mdpd <pair_mesodpd>` command or other
Fix *mvv/dpd* updates the position and velocity of each atom. It can be
used with the :doc:`pair_style mdpd <pair_mesodpd>` command or other
pair styles such as :doc:`pair dpd <pair_dpd>`.
Fix *mvv/edpd* updates the per-atom temperature, in addition to
position and velocity, and must be used with the :doc:`pair_style edpd <pair_mesodpd>` command.
Fix *mvv/edpd* updates the per-atom temperature, in addition to position
and velocity, and must be used with the :doc:`pair_style edpd
<pair_mesodpd>` command.
Fix *mvv/tdpd* updates the per-atom chemical concentration, in
addition to position and velocity, and must be used with the
:doc:`pair_style tdpd <pair_mesodpd>` command.
Fix *mvv/tdpd* updates the per-atom chemical concentration, in addition
to position and velocity, and must be used with the :doc:`pair_style
tdpd <pair_mesodpd>` command.
----------
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files <restart>`. None of the :doc:`fix_modify <fix_modify>` options
are relevant to this fix. No global or per-atom quantities are stored
by this fix for access by various :doc:`output commands <Howto_output>`.
No information about this fix is written to :doc:`binary restart files
<restart>`. None of the :doc:`fix_modify <fix_modify>` options are
relevant to this fix. No global or per-atom quantities are stored by
this fix for access by various :doc:`output commands <Howto_output>`.
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>`.
the :doc:`run <run>` command. This fix is not invoked during
:doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
This fix is part of the DPD-MESO package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package <Build_package>` page for more info.
These fixes are part of the DPD-MESO package. They are only enabled if
LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
Related commands
""""""""""""""""

9
doc/src/fix_nh_uef.rst
View File

@ -218,10 +218,11 @@ use :doc:`change_box <change_box>` before invoking the fix.
Related commands
""""""""""""""""
:doc:`fix nvt <fix_nh>`, :doc:`fix npt <fix_nh>`, `fix nvt/sllod
:doc:<fix_nvt_sllod>`, `compute temp/uef <compute_temp_uef>`,
:doc::doc:`compute pressure/uef <compute_pressure_uef>`, `dump cfg/uef
:doc:<dump_cfg_uef>`
:doc:`fix nvt <fix_nh>`, :doc:`fix npt <fix_nh>`,
:doc:`fix nvt/sllod <fix_nvt_sllod>`,
:doc:`compute temp/uef <compute_temp_uef>`,
:doc:`compute pressure/uef <compute_pressure_uef>`,
:doc:`dump cfg/uef <dump_cfg_uef>`
Default
"""""""

2
doc/src/fix_pimd.rst
View File

@ -1,7 +1,7 @@
.. index:: fix pimd/nvt
fix pimd/nvt command
================
====================
Syntax
""""""

11
doc/src/fix_polarize.rst
View File

@ -31,7 +31,7 @@ Examples
fix 2 interface polarize/bem/gmres 5 0.0001
fix 1 interface polarize/bem/icc 1 0.0001
fix 3 interface polarize/functional 1 0.001
fix 3 interface polarize/functional 1 0.0001
Used in input scripts:
@ -69,8 +69,9 @@ along the normal vector is then 78 - 4 = 74, the mean dielectric value
is (78 + 4) / 2 = 41. Each boundary element also has its area and the
local mean curvature, which is used by these fixes for computing a
correction term in the local electric field. To model charged
interfaces, the interface particle will have a non-zero charge value,
coming from its area and surface charge density.
interfaces, an interface particle will have a non-zero charge value,
coming from its area and surface charge density, and its local dielectric
constant set to the mean dielectric value.
For non-interface particles such as atoms and charged particles, the
interface normal vectors, element area, and dielectric mismatch are
@ -211,6 +212,8 @@ Note that the *polarize/bem/gmres* and *polarize/bem/icc* fixes only
support :doc:`units <units>` *lj*, *real*, *metal*, *si* and *nano* at
the moment.
Note that *polarize/functional* does not yet support charged interfaces.
Related commands
""""""""""""""""
@ -223,7 +226,7 @@ Related commands
Default
"""""""
*iter_max* = 20
*iter_max* = 50
*kspace* = yes

22
doc/src/fix_shake.rst
View File

@ -63,7 +63,7 @@ however, can *only* be applied during molecular dynamics runs.
.. versionchanged:: 15Sep2022
These fixes may still be used during minimization. In that case the
These fixes may now also be used during minimization. In that case the
constraints are *approximated* by strong harmonic restraints.
**SHAKE vs RATTLE:**
@ -133,9 +133,9 @@ constraint lists atom types. All bonds connected to an atom of the
specified type will be constrained. The *m* constraint lists atom
masses. All bonds connected to atoms of the specified masses will be
constrained (within a fudge factor of MASSDELTA specified in
fix_shake.cpp). The *a* constraint lists angle types. If both bonds
in the angle are constrained then the angle will also be constrained
if its type is in the list.
``src/RIGID/fix_shake.cpp``). The *a* constraint lists angle types. If
both bonds in the angle are constrained then the angle will also be
constrained if its type is in the list.
For all constraints, a particular bond is only constrained if both
atoms in the bond are in the group specified with the SHAKE fix.
@ -205,11 +205,11 @@ LAMMPS closely follows (:ref:`Andersen (1983) <Andersen3>`).
The *fix rattle* command modifies forces and velocities and thus
should be defined after all other integration fixes in your input
script. If you define other fixes that modify velocities or forces
after *fix rattle* operates, then *fix rattle* will not take them into
account and the overall time integration will typically not satisfy
the RATTLE constraints. You can check whether the constraints work
correctly by setting the value of RATTLE_DEBUG in src/fix_rattle.cpp
to 1 and recompiling LAMMPS.
after *fix rattle* operates, then *fix rattle* will not take them
into account and the overall time integration will typically not
satisfy the RATTLE constraints. You can check whether the
constraints work correctly by setting the value of RATTLE_DEBUG in
``src/RIGID/fix_rattle.cpp`` to 1 and recompiling LAMMPS.
----------
@ -275,8 +275,8 @@ reducing the :doc:`timestep <timestep>`.
Related commands
""""""""""""""""
`fix rigid <fix_rigid>`, `fix ehex <fix_ehex>`,
`fix nve/manifold/rattle <fix_nve_manifold_rattle>`
:doc:`fix rigid <fix_rigid>`, :doc:`fix ehex <fix_ehex>`,
:doc:`fix nve/manifold/rattle <fix_nve_manifold_rattle>`
Default

6
doc/src/fix_wall.rst
View File

@ -194,6 +194,8 @@ For style *wall/morse*, the energy E is given by a Morse potential:
E = D_0 \left[ e^{- 2 \alpha (r - r_0)} - 2 e^{- \alpha (r - r_0)} \right]
\qquad r < r_c
.. versionadded:: TBD
For style *wall/lepton*, the energy E is provided as an Lepton
expression string using "r" as the distance variable. The `Lepton
library <https://simtk.org/projects/lepton>`_, that the *wall/lepton*
@ -213,6 +215,8 @@ spring as in fix *wall/harmonic* with a force constant *K* (same as
:math:`\epsilon` above) of 100 energy units. More details on the Lepton
expression strings are given below.
.. versionadded:: TBD
For style *wall/table*, the energy E and forces are determined from
interpolation tables listed in one or more files as a function of
distance. The interpolation tables are used to evaluate energy and
@ -371,7 +375,7 @@ is *no*, which means the system must be non-periodic when using a wall.
But you may wish to use a periodic box. E.g. to allow some particles to
interact with the wall via the fix group-ID, and others to pass through
it and wrap around a periodic box. In this case you should ensure that
the wall if sufficiently far enough away from the box boundary. If you
the wall is sufficiently far enough away from the box boundary. If you
do not, then particles may interact with both the wall and with periodic
images on the other side of the box, which is probably not what you
want.

14
doc/src/minimize.rst
View File

@ -49,19 +49,19 @@ and forces) by pushing the atoms off of each other.
The distance that atoms can move during individual minimization steps
can be quite large, especially at the beginning of a minimization.
Thus `neighbor list settings <neigh_modify>` of *every = 1* and
Thus :doc:`neighbor list settings <neigh_modify>` of *every = 1* and
*delay = 0* are **required**. This may be combined with either
*check = no* (always update the neighbor list) or *check = yes* (only
update the neighbor list if at least one atom has moved more than
half the `neighbor list skin <neighbor>` distance since the last
half the :doc:`neighbor list skin <neighbor>` distance since the last
reneighboring). Using *check = yes* is recommended since it avoids
unneeded reneighboring steps when the system is closer to the minimum
and thus atoms move only small distances. Using *check = no* may
be required for debugging or when coupling LAMMPS with external
codes that require a predictable sequence of neighbor list updates.
and thus atoms move only small distances. Using *check = no* may be
required for debugging or when coupling LAMMPS with external codes
that require a predictable sequence of neighbor list updates.
If the settings are **not** *every = 1* and *delay = 0*, LAMMPS
will temporarily apply a `neigh_modify every 1 delay 0 check yes
If the settings are **not** *every = 1* and *delay = 0*, LAMMPS will
temporarily apply a :doc:`neigh_modify every 1 delay 0 check yes
<neigh_modify>` setting during the minimization and restore the
original setting at the end of the minimization. A corresponding
message will be printed to the screen and log file, if this happens.

5
doc/src/molecule.rst
View File

@ -121,6 +121,11 @@ molecule (header keyword = inertia).
ensure space is allocated for storing topology info for molecules that
are added later.
----------
Format of a molecule file
"""""""""""""""""""""""""
The format of an individual molecule file is similar but
(not identical) to the data file read by the :doc:`read_data <read_data>`
commands, and is as follows.

11
doc/src/pair_eam.rst
View File

@ -116,11 +116,12 @@ are parameterized in terms of LAMMPS :doc:`metal units <units>`.
.. note::
Note that unlike for other potentials, cutoffs for EAM
potentials are not set in the pair_style or pair_coeff command; they
are specified in the EAM potential files themselves. Likewise, the
EAM potential files list atomic masses; thus you do not need to use
the :doc:`mass <mass>` command to specify them.
Note that unlike for other potentials, cutoffs for EAM potentials are not
set in the pair_style or pair_coeff command; they are specified in the EAM
potential files themselves. Likewise, valid EAM potential files usually
contain atomic masses; thus you may not need to use the :doc:`mass <mass>`
command to specify them, unless the potential file uses a dummy value (e.g.
0.0). LAMMPS will print a warning, if this is the case.
There are web sites that distribute and document EAM potentials stored
in DYNAMO or other formats:

13
doc/src/pair_gauss.rst
View File

@ -103,12 +103,13 @@ Mixing, shift, table, tail correction, restart, rRESPA info
For atom type pairs I,J and I != J, the A, B, H, sigma_h, r_mh
parameters, and the cutoff distance for these pair styles can be mixed:
A (energy units)
sqrt(1/B) (distance units, see below)
H (energy units)
sigma_h (distance units)
r_mh (distance units)
cutoff (distance units):ul
* A (energy units)
* :math:`\sqrt{\frac{1}{B}}` (distance units, see below)
* H (energy units)
* :math:`r_{mh}` (distance units)
* :math:`\sigma_h` (distance units)
* cutoff (distance units)
The default mix value is *geometric*\ .
Only *arithmetic* and *geometric* mix values are supported.

3
doc/src/pair_hybrid.rst
View File

@ -91,7 +91,8 @@ contributions from sub-styles are weighted by their scale factors, which
may be fractional or even negative. Furthermore the scale factors may
be variables that may change during a simulation. This enables
switching smoothly between two different pair styles or two different
parameter sets during a run.
parameter sets during a run in a similar fashion as could be done
with :doc:`fix adapt <fix_adapt>` or :doc:`fix alchemy <fix_alchemy>`.
All pair styles that will be used are listed as "sub-styles" following
the *hybrid* or *hybrid/overlay* keyword, in any order. In case of the

2
doc/src/pair_mesodpd.rst
View File

@ -299,7 +299,7 @@ Restrictions
""""""""""""
The pair styles *edpd*, *mdpd*, *mdpd/rhosum* and *tdpd* are part of
the DPD-MESO package. It is only enabled if LAMMPS was built with
the DPD-MESO package. They are only enabled if LAMMPS was built with
that package. See the :doc:`Build package <Build_package>` page for
more info.

2
doc/src/pair_tersoff_mod.rst
View File

@ -145,7 +145,7 @@ coefficients in the formulae above:
* c3
* c4
* c5
* c0 (energy units, tersoff/mod/c only):ul
* c0 (energy units, tersoff/mod/c only)
The n, :math:`\eta`, :math:`\lambda_2`, B, :math:`\lambda_1`, and A parameters are only used for
two-body interactions. The :math:`\beta`, :math:`\alpha`, c1, c2, c3, c4, c5, h

9
doc/src/pair_ylz.rst
View File

@ -142,10 +142,11 @@ the :doc:`atom_style ellipsoid <atom_style>` command.
Related commands
""""""""""""""""
:doc:`pair_coeff <pair_coeff>`, :doc:`fix nve/asphere
:doc:<fix_nve_asphere>`, `compute temp/asphere <compute_temp_asphere>`,
:doc::doc:`pair_style resquared <pair_resquared>`, :doc:`pair_style
:doc:gayberne <pair_gayberne>`
:doc:`pair_coeff <pair_coeff>`,
:doc:`fix nve/asphere <fix_nve_asphere>`,
:doc:`compute temp/asphere <compute_temp_asphere>`,
:doc:`pair_style resquared <pair_resquared>`,
:doc:`pair_style gayberne <pair_gayberne>`
Default
"""""""

10
doc/src/read_data.rst
View File

@ -1477,10 +1477,12 @@ The *Triangles* section must appear after the *Atoms* section.
where the keywords have these meanings:
vx,vy,vz = translational velocity of atom
lx,ly,lz = angular momentum of aspherical atom
wx,wy,wz = angular velocity of spherical atom
ervel = electron radial velocity (0 for fixed-core):ul
.. parsed-literal::
vx,vy,vz = translational velocity of atom
lx,ly,lz = angular momentum of aspherical atom
wx,wy,wz = angular velocity of spherical atom
ervel = electron radial velocity (0 for fixed-core)
The velocity lines can appear in any order. This section can only be
used after an *Atoms* section. This is because the *Atoms* section

154
doc/src/run_style.rst
View File

@ -78,27 +78,41 @@ processors. See the :doc:`-partition command-line switch <Run_options>`
for info on how to run LAMMPS with multiple partitions.
Specifically, this style performs all computation except the
:doc:`kspace_style <kspace_style>` portion of the force field on the first
partition. This include the :doc:`pair style <pair_style>`, :doc:`bond style <bond_style>`, :doc:`neighbor list building <neighbor>`,
:doc:`fixes <fix>` including time integration, and output. The
:doc:`kspace_style <kspace_style>` portion of the calculation is
:doc:`kspace_style <kspace_style>` portion of the force field on the
first partition. This include the :doc:`pair style <pair_style>`,
:doc:`bond style <bond_style>`, :doc:`neighbor list building
<neighbor>`, :doc:`fixes <fix>` including time integration, and output.
The :doc:`kspace_style <kspace_style>` portion of the calculation is
performed on the second partition.
This is most useful for the PPPM kspace_style when its performance on
a large number of processors degrades due to the cost of communication
in its 3d FFTs. In this scenario, splitting your P total processors
into 2 subsets of processors, P1 in the first partition and P2 in the
second partition, can enable your simulation to run faster. This is
because the long-range forces in PPPM can be calculated at the same
time as pairwise and bonded forces are being calculated, and the FFTs
can actually speed up when running on fewer processors.
This can lead to a significant speedup, if the number of processors can
be easily increased and the fraction of time is spent in computing
Kspace interactions is significant, too. The two partitions may have a
different number of processors. This is most useful for the PPPM
kspace_style when its performance on a large number of processors
degrades due to the cost of communication in its 3d FFTs. In this
scenario, splitting your P total processors into 2 subsets of
processors, P1 in the first partition and P2 in the second partition,
can enable your simulation to run faster. This is because the
long-range forces in PPPM can be calculated at the same time as pairwise
and bonded forces are being calculated *and* the parallel 3d FFTs can be
faster to compute when running on fewer processors. Please note that
the scenario of using fewer MPI processes to reduce communication
overhead can also be implemented through using MPI with OpenMP threads
via the INTEL, KOKKOS, or OPENMP package. This alternative option is
typically more effective in case of a fixed number of available
processors and less complex to execute.
To use this style, you must define 2 partitions where P1 is a multiple
of P2. Typically having P1 be 3x larger than P2 is a good choice.
The 3d processor layouts in each partition must overlay in the
following sense. If P1 is a Px1 by Py1 by Pz1 grid, and P2 = Px2 by
Py2 by Pz2, then Px1 must be an integer multiple of Px2, and similarly
for Py1 a multiple of Py2, and Pz1 a multiple of Pz2.
To use the *verlet/split* style, you must define 2 partitions with the
:doc:`-partition command-line switch <Run_options>`, where partition P1
is either the same size or an integer multiple of the size of the
partition P2. Typically having P1 be 3x larger than P2 is a good
choice, since the (serial) performance of LAMMPS is often best if the
time spent in the ``Pair`` computation versus ``Kspace`` is a 3:1 split.
The 3d processor layouts in each partition must overlay in the following
sense. If P1 is a Px1 by Py1 by Pz1 grid, and P2 = Px2 by Py2 by Pz2,
then Px1 must be an integer multiple of Px2, and similarly for Py1 a
multiple of Py2, and Pz1 a multiple of Pz2.
Typically the best way to do this is to let the first partition choose
its own optimal layout, then require the second partition's layout to
@ -122,9 +136,10 @@ of 60 and 15 processors each:
When you run in 2-partition mode with the *verlet/split* style, the
thermodynamic data for the entire simulation will be output to the log
and screen file of the first partition, which are log.lammps.0 and
screen.0 by default; see the :doc:`-plog and -pscreen command-line switches <Run_options>` to change this. The log and screen file
for the second partition will not contain thermodynamic output beyond the
first timestep of the run.
screen.0 by default; see the :doc:`-plog and -pscreen command-line
switches <Run_options>` to change this. The log and screen file for the
second partition will not contain thermodynamic output beyond the first
timestep of the run.
See the :doc:`Accelerator packages <Speed_packages>` page for
performance details of the speed-up offered by the *verlet/split*
@ -137,70 +152,73 @@ options to support this, and strategies are discussed in :doc:`Section
----------
The *respa* style implements the rRESPA multi-timescale integrator
:ref:`(Tuckerman) <Tuckerman3>` with N hierarchical levels, where level 1 is
the innermost loop (shortest timestep) and level N is the outermost
:ref:`(Tuckerman) <Tuckerman3>` with N hierarchical levels, where level
1 is the innermost loop (shortest timestep) and level N is the outermost
loop (largest timestep). The loop factor arguments specify what the
looping factor is between levels. N1 specifies the number of
iterations of level 1 for a single iteration of level 2, N2 is the
iterations of level 2 per iteration of level 3, etc. N-1 looping
parameters must be specified.
looping factor is between levels. N1 specifies the number of iterations
of level 1 for a single iteration of level 2, N2 is the iterations of
level 2 per iteration of level 3, etc. N-1 looping parameters must be
specified.
Thus with a 4-level respa setting of "2 2 2" for the 3 loop factors,
you could choose to have bond interactions computed 8x per large
timestep, angle interactions computed 4x, pair interactions computed
2x, and long-range interactions once per large timestep.
Thus with a 4-level respa setting of "2 2 2" for the 3 loop factors, you
could choose to have bond interactions computed 8x per large timestep,
angle interactions computed 4x, pair interactions computed 2x, and
long-range interactions once per large timestep.
The :doc:`timestep <timestep>` command sets the large timestep for the
outermost rRESPA level. Thus if the 3 loop factors are "2 2 2" for
4-level rRESPA, and the outer timestep is set to 4.0 fs, then the
inner timestep would be 8x smaller or 0.5 fs. All other LAMMPS
commands that specify number of timesteps (e.g. :doc:`thermo <thermo>`
for thermo output every N steps, :doc:`neigh_modify delay/every <neigh_modify>` parameters, :doc:`dump <dump>` every N
steps, etc) refer to the outermost timesteps.
4-level rRESPA, and the outer timestep is set to 4.0 fs, then the inner
timestep would be 8x smaller or 0.5 fs. All other LAMMPS commands that
specify number of timesteps (e.g. :doc:`thermo <thermo>` for thermo
output every N steps, :doc:`neigh_modify delay/every <neigh_modify>`
parameters, :doc:`dump <dump>` every N steps, etc) refer to the
outermost timesteps.
The rRESPA keywords enable you to specify at what level of the
hierarchy various forces will be computed. If not specified, the
defaults are that bond forces are computed at level 1 (innermost
loop), angle forces are computed where bond forces are, dihedral
forces are computed where angle forces are, improper forces are
computed where dihedral forces are, pair forces are computed at the
outermost level, and kspace forces are computed where pair forces are.
The inner, middle, outer forces have no defaults.
The rRESPA keywords enable you to specify at what level of the hierarchy
various forces will be computed. If not specified, the defaults are
that bond forces are computed at level 1 (innermost loop), angle forces
are computed where bond forces are, dihedral forces are computed where
angle forces are, improper forces are computed where dihedral forces
are, pair forces are computed at the outermost level, and kspace forces
are computed where pair forces are. The inner, middle, outer forces
have no defaults.
For fixes that support it, the rRESPA level at which a given fix is
active, can be selected through the :doc:`fix_modify <fix_modify>` command.
active, can be selected through the :doc:`fix_modify <fix_modify>`
command.
The *inner* and *middle* keywords take additional arguments for
cutoffs that are used by the pairwise force computations. If the 2
cutoffs for *inner* are 5.0 and 6.0, this means that all pairs up to
6.0 apart are computed by the inner force. Those between 5.0 and 6.0
have their force go ramped to 0.0 so the overlap with the next regime
(middle or outer) is smooth. The next regime (middle or outer) will
compute forces for all pairs from 5.0 outward, with those from 5.0 to
6.0 having their value ramped in an inverse manner.
The *inner* and *middle* keywords take additional arguments for cutoffs
that are used by the pairwise force computations. If the 2 cutoffs for
*inner* are 5.0 and 6.0, this means that all pairs up to 6.0 apart are
computed by the inner force. Those between 5.0 and 6.0 have their force
go ramped to 0.0 so the overlap with the next regime (middle or outer)
is smooth. The next regime (middle or outer) will compute forces for
all pairs from 5.0 outward, with those from 5.0 to 6.0 having their
value ramped in an inverse manner.
Note that you can use *inner* and *outer* without using *middle* to
split the pairwise computations into two portions instead of three.
Unless you are using a very long pairwise cutoff, a 2-way split is
often faster than a 3-way split, since it avoids too much duplicate
Unless you are using a very long pairwise cutoff, a 2-way split is often
faster than a 3-way split, since it avoids too much duplicate
computation of pairwise interactions near the intermediate cutoffs.
Also note that only a few pair potentials support the use of the
*inner* and *middle* and *outer* keywords. If not, only the *pair*
keyword can be used with that pair style, meaning all pairwise forces
are computed at the same rRESPA level. See the doc pages for
individual pair styles for details.
Also note that only a few pair potentials support the use of the *inner*
and *middle* and *outer* keywords. If not, only the *pair* keyword can
be used with that pair style, meaning all pairwise forces are computed
at the same rRESPA level. See the doc pages for individual pair styles
for details.
Another option for using pair potentials with rRESPA is with the
*hybrid* keyword, which requires the use of the :doc:`pair_style hybrid or hybrid/overlay <pair_hybrid>` command. In this scenario, different
*hybrid* keyword, which requires the use of the :doc:`pair_style hybrid
or hybrid/overlay <pair_hybrid>` command. In this scenario, different
sub-styles of the hybrid pair style are evaluated at different rRESPA
levels. This can be useful, for example, to set different timesteps
for hybrid coarse-grained/all-atom models. The *hybrid* keyword
requires as many level assignments as there are hybrid sub-styles,
which assigns each sub-style to a rRESPA level, following their order
of definition in the pair_style command. Since the *hybrid* keyword
operates on pair style computations, it is mutually exclusive with
either the *pair* or the *inner*\ /\ *middle*\ /\ *outer* keywords.
levels. This can be useful, for example, to set different timesteps for
hybrid coarse-grained/all-atom models. The *hybrid* keyword requires as
many level assignments as there are hybrid sub-styles, which assigns
each sub-style to a rRESPA level, following their order of definition in
the pair_style command. Since the *hybrid* keyword operates on pair
style computations, it is mutually exclusive with either the *pair* or
the *inner*\ /\ *middle*\ /\ *outer* keywords.
When using rRESPA (or for any MD simulation) care must be taken to
choose a timestep size(s) that ensures the Hamiltonian for the chosen