Merge pull request lammps#4050 from akohlmey/bioff-doc-updates

Bio force field documentation updates

.gitignore

@@ -60,3 +60,4 @@ src/Makefile.package.settings-e
 /cmake/build/x64-Debug-Clang
 /install/x64-GUI-MSVC
 /install
+.Rhistory

doc/src/Howto_bioFF.rst

@@ -11,12 +11,30 @@ commands like :doc:`pair_coeff <pair_coeff>` or :doc:`bond_coeff
 <bond_coeff>` and so on.  See the :doc:`Tools <Tools>` doc page for
 additional tools that can use CHARMM, AMBER, or Materials Studio
 generated files to assign force field coefficients and convert their
-output into LAMMPS input.
+output into LAMMPS input. LAMMPS input scripts can also be generated by `charmm-gui.org <https://charmm-gui.org/>`_.
 
-See :ref:`(MacKerell) <howto-MacKerell>` for a description of the CHARMM
-force field.  See :ref:`(Cornell) <howto-Cornell>` for a description of
-the AMBER force field.  See :ref:`(Sun) <howto-Sun>` for a description
-of the COMPASS force field.
+CHARMM and AMBER
+----------------
+
+The `CHARMM force field <https://mackerell.umaryland.edu/charmm_ff.shtml>`_ :ref:`(MacKerell) <howto-MacKerell>` and `AMBER force field <https://ambermd.org/AmberModels.php>`_ :ref:`(Cornell) <howto-Cornell>` have potential energy function of the form
+
+.. math::
+
+  V & = \sum_{bonds} E_b + \sum_{angles} \!E_a + \!\overbrace{\sum_{dihedral} \!\!E_d}^{\substack{
+         \text{charmm} \\
+        \text{charmmfsw}
+      }} +\!\!\! \sum_{impropers} \!\!\!E_i \\[.6em]
+      & \quad + \!\!\!\!\!\!\!\!\!\!\underbrace{~\sum_{pairs} \left(E_{LJ}+E_{coul}\right)}_{\substack{
+         \text{lj/charmm/coul/charmm} \\
+        \text{lj/charmm/coul/charmm/implicit} \\
+        \text{lj/charmm/coul/long} \\
+        \text{lj/charmm/coul/msm} \\
+         \text{lj/charmmfsw/coul/charmmfsh} \\
+        \text{lj/charmmfsw/coul/long}
+      }} \!\!\!\!\!\!\!\!+ \!\!\sum_{special}\! E_s + \!\!\!\!\sum_{residues} \!\!\!{\scriptstyle\mathrm{CMAP}(\phi,\psi)}
+
+
+The terms are computed by bond styles (relationship between 2 atoms), angle styles (between 3 atoms) , dihedral/improper styles (between 4 atoms), pair styles (non-covalently bonded pair interactions) and special bonds. The CMAP term (see :doc:`fix cmap <fix_cmap>` command for details) corrects for pairs of dihedral angles ("Correction MAP") to significantly improve the structural and dynamic properties of proteins in crystalline and solution environments :ref:`(Brooks) <howto-Brooks>`. The AMBER force field does not include the CMAP term.
 
 The interaction styles listed below compute force field formulas that
 are consistent with common options in CHARMM or AMBER.  See each
@@ -31,10 +49,61 @@ command's documentation for the formula it computes.
 * :doc:`pair_style <pair_charmm>` lj/charmm/coul/charmm
 * :doc:`pair_style <pair_charmm>` lj/charmm/coul/charmm/implicit
 * :doc:`pair_style <pair_charmm>` lj/charmm/coul/long
-
 * :doc:`special_bonds <special_bonds>` charmm
 * :doc:`special_bonds <special_bonds>` amber
 
+The pair styles compute Lennard Jones (LJ) and Coulombic interactions with additional switching or shifting functions that ramp the energy and/or force smoothly to zero between an inner :math:`(a)` and outer :math:`(b)` cutoff. The older styles with *charmm* (not *charmmfsw* or *charmmfsh*\ ) in their name compute the LJ and Coulombic interactions with an energy switching function (esw) S(r) which ramps the energy smoothly to zero between the inner and outer cutoff. This can cause irregularities in pairwise forces (due to the discontinuous second derivative of energy at the boundaries of the switching region), which in some cases can result in complications in energy minimization and detectable artifacts in MD simulations.
+
+.. math::
+
+   LJ(r) &= 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} -
+           \left(\frac{\sigma}{r}\right)^6 \right]\\[.6em]
+   C(r) &= \frac{C q_i q_j}{ \epsilon r}\\[.6em]
+   S(r) &=  \frac{ \left(b^2 - r^2\right)^2
+     \left(b^2 + 2r^2 - 3{a^2}\right)}
+   { \left(b^2 - a^2\right)^3 }\\[.6em]
+   E_{LJ}(r) &=  \begin{cases}
+        LJ(r), & r \leq a \\
+        LJ(r) S(r), & a < r \leq b \\
+        0, &r > b
+      \end{cases} \\[.6em]
+   E_{coul}(r) &=  \begin{cases}
+        C(r), & r \leq a \\
+        C(r) S(r), & a < r \leq b \\
+        0, & r > b
+      \end{cases}
+
+.. image:: img/howto_charmm_ELJ.png
+  :align: center
+
+|
+
+The newer styles with *charmmfsw* or *charmmfsh* in their name replace energy switching with force switching (fsw) for LJ interactions and force shifting (fsh) functions for Coulombic interactions :ref:`(Steinbach) <howto-Steinbach>`
+
+.. math::
+
+     E_{LJ}(r) = & \begin{cases}
+ 4  \epsilon \sigma^6  \left(\frac{\displaystyle\sigma
+   ^6-r^6}{\displaystyle r^{12}}-\frac{\displaystyle\sigma ^6}{\displaystyle a^6
+   b^6}+\frac{\displaystyle 1}{\displaystyle a^3 b^3}\right) & r\leq a \\
+ \frac{\displaystyle 4 \epsilon \sigma^6   \left(\sigma ^6
+   \left(b^6-r^6\right)^2-b^3 r^6 \left(a^3+b^3\right)
+   \left(b^3-r^3\right)^2\right)}{\displaystyle b^6 r^{12}
+   \left(b^6-a^6\right)} & a<r \leq b\\
+   0, & r>b
+  \end{cases}\\[.6em]
+   E_{coul}(r) & =  \begin{cases}
+        C(r) \frac{\displaystyle (b-r)^2}{\displaystyle r b^2}, &  r \leq b \\
+        0, & r > b
+      \end{cases}
+
+.. image:: img/howto_charmmfsw_ELJ.png
+  :align: center
+
+|
+
+These styles are used by LAMMPS input scripts generated by `charmm-gui.org <https://charmm-gui.org/>`_ :ref:`(Brooks) <howto-Brooks>`. A `minimal PDB example 1HVN <https://www.rcsb.org/structure/1HVN>`_ with at least one protein segment, at least one DNA segment, and no modified engineered residues is available in the ``lammps/examples/charmm/1hvn`` directory. A better example is `PDB 2CV5 <https://www.rcsb.org/structure/2CV5>`_ with size too big to include in lammps examples, which is left as an exercise to the reader (go to charmm-gui.org and type in 2CV5 in PDB field of Solution Builder to generate LAMMPS scripts to simulate a solvated human nucleosome with histone octamer and dsDNA wrapped around it).
+
 .. note::
 
    For CHARMM, newer *charmmfsw* or *charmmfsh* styles were released in
@@ -43,9 +112,16 @@ command's documentation for the formula it computes.
    <pair_charmm>` and :doc:`dihedral charmm <dihedral_charmm>` doc
    pages.
 
+.. note::
+
+  TIP3P water model MUST be used with CHARMM force field not TIP4P, TIP5P or SPC. In fact, `"using the SPC model with CHARMM parameters is a bad idea" <https://matsci.org/t/using-spc-water-with-charmm-ff/24715>`_ and `"to enable TIP4P style water in CHARMM, you would have to write a new pair style" <https://matsci.org/t/hybrid-pair-styles-for-charmm-and-tip4p-ew/32609>`_ . LAMMPS input scripts generated by Solution Builder on charmm-gui.org use TIP3P molecules for solvation. Any other water model can and probably will lead to false conclusions.
+
+COMPASS
+-------
+
 COMPASS is a general force field for atomistic simulation of common
 organic molecules, inorganic small molecules, and polymers which was
-developed using ab initio and empirical parameterization techniques.
+developed using ab initio and empirical parameterization techniques :ref:`(Sun) <howto-Sun>`.
 See the :doc:`Tools <Tools>` page for the msi2lmp tool for creating
 LAMMPS template input and data files from BIOVIA's Materials Studio
 files.  Please note that the msi2lmp tool is very old and largely
@@ -70,6 +146,9 @@ documentation for the formula it computes.
 
 * :doc:`special_bonds <special_bonds>` lj/coul 0 0 1
 
+DREIDING
+--------
+
 DREIDING is a generic force field developed by the `Goddard group <http://www.wag.caltech.edu>`_ at Caltech and is useful for
 predicting structures and dynamics of organic, biological and main-group
 inorganic molecules. The philosophy in DREIDING is to use general force
@@ -113,18 +192,25 @@ documentation for the formula it computes.
 .. _howto-MacKerell:
 
 **(MacKerell)** MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field,
-Fischer, Gao, Guo, Ha, et al, J Phys Chem, 102, 3586 (1998).
+Fischer, Gao, Guo, Ha, et al (1998). All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J Phys Chem, 102, 3586 . https://doi.org/10.1021/jp973084f
 
 .. _howto-Cornell:
 
 **(Cornell)** Cornell, Cieplak, Bayly, Gould, Merz, Ferguson,
-Spellmeyer, Fox, Caldwell, Kollman, JACS 117, 5179-5197 (1995).
+Spellmeyer, Fox, Caldwell, Kollman (1995). A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. JACS 117, 5179-5197. https://doi.org/10.1021/ja00124a002
+
+.. _howto-Steinbach:
+
+**(Steinbach)** Steinbach, Brooks (1994). New spherical-cutoff methods for long-range forces in macromolecular simulation. J Comput Chem, 15, 667. https://doi.org/10.1002/jcc.540150702
+
+.. _howto-Brooks:
+
+**(Brooks)** Brooks, et al (2009). CHARMM: The biomolecular simulation program. J Comput Chem, 30, 1545. https://onlinelibrary.wiley.com/doi/10.1002/jcc.21287
 
 .. _howto-Sun:
 
-**(Sun)** Sun, J. Phys. Chem. B, 102, 7338-7364 (1998).
+**(Sun)** Sun (1998). COMPASS: An ab Initio Force-Field Optimized for Condensed-Phase ApplicationsOverview with Details on Alkane and Benzene Compounds. J. Phys. Chem. B, 102, 7338-7364. https://doi.org/10.1021/jp980939v
 
 .. _howto-Mayo:
 
-**(Mayo)** Mayo, Olfason, Goddard III, J Phys Chem, 94, 8897-8909
-(1990).
+**(Mayo)** Mayo, Olfason, Goddard III (1990). DREIDING: a generic force field for molecular simulations. J Phys Chem, 94, 8897-8909. https://doi.org/10.1021/j100389a010

doc/src/pair_charmm.rst

@@ -112,26 +112,22 @@ Description
 These pair styles compute Lennard Jones (LJ) and Coulombic
 interactions with additional switching or shifting functions that ramp
 the energy and/or force smoothly to zero between an inner and outer
-cutoff.  They are implementations of the widely used CHARMM force
-field used in the `CHARMM <https://www.charmm.org>`_ MD code (and
-others).  See :ref:`(MacKerell) <pair-MacKerell>` for a description of the
-CHARMM force field.
+cutoff.  They implement the widely used CHARMM force field, see
+:doc:`Howto discussion on biomolecular force fields <Howto_bioFF>` for
+details.
 
 The styles with *charmm* (not *charmmfsw* or *charmmfsh*\ ) in their
 name are the older, original LAMMPS implementations.  They compute the
-LJ and Coulombic interactions with an energy switching function (esw,
-shown in the formula below as S(r)), which ramps the energy smoothly
-to zero between the inner and outer cutoff.  This can cause
-irregularities in pairwise forces (due to the discontinuous second
-derivative of energy at the boundaries of the switching region), which
-in some cases can result in detectable artifacts in an MD simulation.
+LJ and Coulombic interactions with an energy switching function which
+ramps the energy smoothly to zero between the inner and outer cutoff.
+This can cause irregularities in pairwise forces (due to the discontinuous
+second derivative of energy at the boundaries of the switching region),
+which in some cases can result in detectable artifacts in an MD simulation.
 
 The newer styles with *charmmfsw* or *charmmfsh* in their name replace
 the energy switching with force switching (fsw) and force shifting
 (fsh) functions, for LJ and Coulombic interactions respectively.
-These follow the formulas and description given in
-:ref:`(Steinbach) <Steinbach>` and :ref:`(Brooks) <Brooks1>` to minimize these
-artifacts.
+
 
 .. note::
 
@@ -152,26 +148,6 @@ artifacts.
    the CHARMM force field energies and forces, when using one of these
    two CHARMM pair styles.
 
-.. math::
-
-   E = & LJ(r) \qquad \qquad \qquad r < r_{\rm in} \\
-     = & S(r) * LJ(r) \qquad \qquad r_{\rm in} < r < r_{\rm out} \\
-     = & 0 \qquad \qquad \qquad \qquad r > r_{\rm out} \\
-   E = & C(r) \qquad \qquad \qquad r < r_{\rm in} \\
-     = & S(r) * C(r) \qquad \qquad r_{\rm in} < r < r_{\rm out} \\
-     = & 0 \qquad \qquad \qquad \qquad r > r_{\rm out} \\
-   LJ(r) = & 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} -
-           \left(\frac{\sigma}{r}\right)^6 \right] \\
-   C(r) = & \frac{C q_i q_j}{ \epsilon r} \\
-   S(r) = & \frac{ \left[r_{\rm out}^2 - r^2\right]^2
-     \left[r_{\rm out}^2 + 2r^2 - 3{r_{\rm in}^2}\right]}
-   { \left[r_{\rm out}^2 - {r_{\rm in}}^2\right]^3 }
-
-where S(r) is the energy switching function mentioned above for the
-*charmm* styles.  See the :ref:`(Steinbach) <Steinbach>` paper for the
-functional forms of the force switching and force shifting functions
-used in the *charmmfsw* and *charmmfsh* styles.
-
 When using the *lj/charmm/coul/charmm styles*, both the LJ and
 Coulombic terms require an inner and outer cutoff. They can be the
 same for both formulas or different depending on whether 2 or 4