diff --git a/doc/src/Howto.rst b/doc/src/Howto.rst index f229deafb3..1366ecb839 100644 --- a/doc/src/Howto.rst +++ b/doc/src/Howto.rst @@ -70,6 +70,7 @@ Force fields howto Howto_amoeba Howto_tip3p Howto_tip4p + Howto_tip5p Howto_spc Packages howto diff --git a/doc/src/Howto_spc.rst b/doc/src/Howto_spc.rst index 520942547e..6414d3b846 100644 --- a/doc/src/Howto_spc.rst +++ b/doc/src/Howto_spc.rst @@ -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) ` 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 +` 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 `_. ---------- .. _howto-Berendsen: -**(Berendsen)** Berendsen, Grigera, Straatsma, J Phys Chem, 91, -6269-6271 (1987). +**(Berendsen)** Berendsen, Grigera, Straatsma, J Phys Chem, 91, 6269-6271 (1987). diff --git a/doc/src/Howto_tip3p.rst b/doc/src/Howto_tip3p.rst index 5571b3291a..682c7f2640 100644 --- a/doc/src/Howto_tip3p.rst +++ b/doc/src/Howto_tip3p.rst @@ -1,53 +1,211 @@ TIP3P water model ================= -The TIP3P water model as implemented in CHARMM -:ref:`(MacKerell) ` 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 ` 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) +` 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) `. +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 ` -These are the parameters to use for TIP3P with a long-range Coulomb -solver (e.g. Ewald or PPPM in LAMMPS), see :ref:`(Price) ` for -details: +or these commands for a long-range Coulomb model: + +* :doc:`pair_style lj/cut/coul/long ` +* :doc:`pair_style lj/cut/coul/long/soft ` +* :doc:`kspace_style pppm ` +* :doc:`kspace_style pppm/disp ` + +In LAMMPS the :doc:`fix shake or fix rattle ` command can be +used to hold the two O-H bonds and the H-O-H angle rigid. A bond style +of :doc:`harmonic ` and an angle style of :doc:`harmonic +` or :doc:`charmm ` should also be used. +In case of rigid bonds also bond style :doc:`zero ` and angle +style :doc:`zero ` can be used. + +The table below lists the force field parameters (in real :doc:`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) `, or a TIP3P model with parameters +optimized for a long-range Coulomb solver (e.g. Ewald or PPPM in LAMMPS) +:ref:`(Price) `. 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 +` 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 `_. diff --git a/doc/src/Howto_tip4p.rst b/doc/src/Howto_tip4p.rst index 19b811a282..7912557e61 100644 --- a/doc/src/Howto_tip4p.rst +++ b/doc/src/Howto_tip4p.rst @@ -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 ` 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 ` and an angle style of :doc:`harmonic +` or :doc:`charmm ` should also be used. +In case of rigid bonds also bond style :doc:`zero ` and angle +style :doc:`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 ` -* :doc:`pair_style lj/cut/tip4p/cut ` +#. Use a specially written pair style that uses the :ref:`TIP3P geometry + ` 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 ` -* :doc:`pair_style lj/cut/tip4p/long ` -* :doc:`pair_style lj/long/tip4p/long ` -* :doc:`pair_style tip4p/long/soft ` -* :doc:`pair_style lj/cut/tip4p/long/soft ` -* :doc:`kspace_style pppm/tip4p ` -* :doc:`kspace_style pppm/disp/tip4p ` + * :doc:`pair_style tip4p/cut ` + * :doc:`pair_style lj/cut/tip4p/cut ` -The bond lengths and bond angles should be held fixed using the -:doc:`fix shake ` or :doc:`fix rattle ` 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) `. Note that the OM distance is specified in -the :doc:`pair_style ` command, not as part of the pair -coefficients. + * :doc:`pair_style tip4p/long ` + * :doc:`pair_style lj/cut/tip4p/long ` + * :doc:`pair_style lj/long/tip4p/long ` + * :doc:`pair_style tip4p/long/soft ` + * :doc:`pair_style lj/cut/tip4p/long/soft ` + * :doc:`kspace_style pppm/tip4p ` + * :doc:`kspace_style pppm/disp/tip4p ` -| 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 ` or :doc:`fix rattle ` 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 ` where + the oxygen atom carries no charge and the M point no Lennard-Jones + interactions. Since :doc:`fix shake ` or :doc:`fix rattle + ` may not be applied to this kind of geometry, :doc:`fix + rigid or fix rigid/small ` 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 +`) to for a selection of popular variants of the TIP4P model. +There is the rigid TIP4P model with a cutoff :ref:`(Jorgensen) +`, the TIP4/Ice model :ref:`(Abascal1) `, the +TIP4P/2005 model :ref:`(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 ` 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 ` command. +Below is the code for a LAMMPS input file using the implicit method and +the :ref:`TIP3P molecule file `. 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 +` 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 ` 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 + 2 -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 `_. ---------- @@ -116,3 +265,13 @@ Wikipedia also has a nice article on `water models ` 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 +`. + +A suitable pair style with cutoff Coulomb would be: + +* :doc:`pair_style lj/cut/coul/cut ` + +or these commands for a long-range model: + +* :doc:`pair_style lj/cut/coul/long ` +* :doc:`pair_style lj/cut/coul/long/soft ` +* :doc:`kspace_style pppm ` +* :doc:`kspace_style pppm/disp ` + +A TIP5P model *must* be run using a :doc:`rigid fix ` 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 +`) to for a the TIP5P model with a cutoff :ref:`(Mahoney) +` and the TIP5P-E model :ref:`(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 ` 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 ` +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 `_. + +---------- + +.. _Mahoney: + +**(Mahoney)** Mahoney, Jorgensen, J Chem Phys 112, 8910 (2000) + +.. _Rick: + +**(Rick)** Rick, J Chem Phys 120, 6085 (2004) diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index fdacbf9c4c..0349c894a2 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -1,5 +1,6 @@ aa aat +Abascal abc abf ABI @@ -1985,6 +1986,7 @@ magelec Maginn magneton magnetons +Mahoney mainboard mainboards makefile @@ -3186,6 +3188,7 @@ Sandia sandybrown sanitizer Sanyal +Sanz Sarath sc scafacos