redesign the Howto pages on water models and add inputs and molecule files
This commit is contained in:
Axel Kohlmeyer
@ -70,6 +70,7 @@ Force fields howto
|
||||
Howto_amoeba
|
||||
Howto_tip3p
|
||||
Howto_tip4p
|
||||
Howto_tip5p
|
||||
Howto_spc
|
||||
|
||||
Packages howto
|
||||
|
||||
@ -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).
|
||||
|
||||
@ -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>`_.
|
||||
|
||||
|
||||
@ -2,19 +2,33 @@ 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:
|
||||
|
||||
#. 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.
|
||||
|
||||
This can be done with the following pair styles for Coulomb with a cutoff:
|
||||
|
||||
* :doc:`pair_style tip4p/cut <pair_lj_cut_tip4p>`
|
||||
* :doc:`pair_style lj/cut/tip4p/cut <pair_lj_cut_tip4p>`
|
||||
|
||||
or these commands for a long-range model:
|
||||
or these commands for a long-range Coulomb treatment:
|
||||
|
||||
* :doc:`pair_style tip4p/long <pair_coul>`
|
||||
* :doc:`pair_style lj/cut/tip4p/long <pair_lj_cut_tip4p>`
|
||||
@ -31,71 +45,95 @@ 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.
|
||||
|
||||
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.
|
||||
#. 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.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 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 TIP4/Ice model (J Chem Phys, 122, 234511 (2005);
|
||||
https://doi.org/10.1063/1.1931662) these values can be used:
|
||||
.. list-table::
|
||||
:header-rows: 1
|
||||
:widths: auto
|
||||
|
||||
| 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
|
||||
|
|
||||
|
||||
For the TIP4P/2005 model (J Chem Phys, 123, 234505 (2005);
|
||||
https://doi.org/10.1063/1.2121687), these values can be used:
|
||||
|
||||
| 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
|
||||
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 <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
161
doc/src/Howto_tip5p.rst
Normal 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)
|
||||
@ -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
|
||||
|
||||
Reference in New Issue
Block a user