Improve manual from comments by Gang Liu
authorErik Lindahl <erik@kth.se>
Sat, 30 Dec 2017 21:28:47 +0000 (22:28 +0100)
committerMark Abraham <mark.j.abraham@gmail.com>
Tue, 2 Jan 2018 16:10:03 +0000 (17:10 +0100)
Fixes #2078.

Change-Id: If1267822be4963a59bac4dc07f9731c9489cc3e3

docs/manual/algorithms.tex
docs/manual/forcefield.tex
docs/manual/implement.tex
docs/manual/monster.bib

index 7e748bc534a7c2358c6c20b89e680f5e7ad2e19e..514ae907c047f531acd39d8ba9a42dba9d3f16a3 100644 (file)
@@ -1,7 +1,7 @@
 %
 % This file is part of the GROMACS molecular simulation package.
 %
-% Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
+% Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
 % Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
 % and including many others, as listed in the AUTHORS file in the
 % top-level source directory and at http://www.gromacs.org.
@@ -1550,7 +1550,24 @@ The equations of motion for the atoms, obtained from the Hamiltonian are:
 \bea \frac {\de^2\ve{r}_i}{\de t^2} & = & \frac{\ve{F}_i}{m_i} -
 \ve{M} \frac{\de \ve{r}_i}{\de t} , \\ \ve{M} & = & \ve{b}^{-1} \left[
   \ve{b} \frac{\de \ve{b}'}{\de t} + \frac{\de \ve{b}}{\de t} \ve{b}'
-  \right] \ve{b}'^{-1}.  \eea The (inverse) mass parameter matrix
+  \right] \ve{b}'^{-1}.
+  \eea 
+This extra term has the appearance of a friction, but it should be 
+noted that it is ficticious, and rather an effect of the
+Parrinello-Rahman equations of motion being defined with all
+particle coordinates represented relative to the box vectors, while
+{\gromacs] uses normal Cartesian coordinates for positions, 
+velocities and forces. It is worth noting that the kinetic energy too
+should formally be calculated based on velocities relative to the
+box vectors. This can have an effect e.g. for external constant stress,
+but for now we only support coupling to constant external 
+pressures, and for any normal simulation the velocities of box
+vectors should be extremely small compared to particle velocities.
+Gang Liu has done some work on deriving this for Cartesian
+coordinates\cite{Liu2015} that we will try to implement at
+some point in the future together with support for external stress.
+
+The (inverse) mass parameter matrix
 $\ve{W}^{-1}$ determines the strength of the coupling, and how the box
 can be deformed.  The box restriction (\ref{eqn:box_rot}) will be
 fulfilled automatically if the corresponding elements of $\ve{W}^{-1}$
index 8574bc409459fe2b117bfe8df3b7c4d1ee615ba4..f50da64e64a21cd5a34c936b76221430f4905fce 100644 (file)
@@ -1,7 +1,7 @@
 %
 % This file is part of the GROMACS molecular simulation package.
 %
-% Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
+% Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
 % Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
 % and including many others, as listed in the AUTHORS file in the
 % top-level source directory and at http://www.gromacs.org.
@@ -56,14 +56,19 @@ externally applied forces, see \chref{special}.
 \end{enumerate}
 
 \section{Non-bonded interactions}
-Non-bonded interactions in {\gromacs} are pair-additive and centro-symmetric:
+Non-bonded interactions in {\gromacs} are pair-additive:
 \beq
 V(\ve{r}_1,\ldots \ve{r}_N) = \sum_{i<j}V_{ij}(\rvij);
 \eeq
 \beq
-\ve{F}_i = -\sum_j \frac{dV_{ij}(r_{ij})}{dr_{ij}} \frac{\rvij}{r_{ij}} = -\ve{F}_j
+\ve{F}_i = -\sum_j \frac{dV_{ij}(r_{ij})}{dr_{ij}} \frac{\rvij}{r_{ij}}
 \eeq
-The non-bonded interactions contain a \normindex{repulsion} term, 
+Since the potential only depends on the scalar distance, interactions
+will be centro-symmetric, i.e.\ the vectorial partial force on particle $i$ from
+the pairwise interaction $V_{ij}(r_{ij})$ has the opposite direction of the partial force on
+particle $j$. For efficiency reasons, interactions are calculated by loops over interactions and
+updating both partial forces rather than summing one complete nonbonded force at 
+a time. The non-bonded interactions contain a \normindex{repulsion} term, 
 a \normindex{dispersion}
 term, and a Coulomb term. The repulsion and dispersion term are
 combined in either the Lennard-Jones (or 6-12 interaction), or the
index d72d122ca01d320efafc950d314520398bb56db6..25c3d9798f822faf95e7c3de266d900a1e44c18a 100644 (file)
@@ -1,7 +1,7 @@
 %
 % This file is part of the GROMACS molecular simulation package.
 %
-% Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
+% Copyright (c) 2013,2014,2015,2017,2018, by the GROMACS development team, led by
 % Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
 % and including many others, as listed in the AUTHORS file in the
 % top-level source directory and at http://www.gromacs.org.
@@ -75,20 +75,55 @@ octahedron is used, there are 15 possible images.
 
 \subsection{Virial from non-bonded forces}
 Here the derivation for the single sum virial in the {\em non-bonded force} 
-routine is given. $i \neq j$ in all formulae below.
+routine is given. There are a couple of considerations that are special
+to {\gromacs} that we take into account:
+\begin{itemize}
+\item When calculating short-range interactions, we apply the 
+{\em minimum image convention} and only consider the closest 
+image of each neighbor - and in particular we never allow interactions
+between a particle and any of its periodic images. For all the 
+equations below, this means $i \neq j$.
+\item In general, either the $i$ or $j$ particle might be shifted to a neighbor
+cell to get the closest interaction (shift $\delta_{ij}$). However, with minimum image
+convention there can be at most 27 different shifts for particles in the central cell,
+and for typical (very short-ranged) biomolecular interactions there are typically only a few
+different shifts involved for each particle, not to mention that each interaction can
+only be present for one shift. 
+\item For the {\gromacs} nonbonded interactions
+we use this to split the neighborlist of each $i$ particle into multiple 
+separate lists, where each list has a constant shift $\delta_i$ for the $i$ partlcle. We
+can represent this as a sum over shifts (for which we use index $s$), with the constraint that
+each particle interaction can only contribute to one of the terms in this sum, and the
+shift is no longer dependent on the $j$ particles. For any sum that does not contain
+complex dependence on $s$, this means the sum trivially reduces to just the sum
+over $i$ and/or $j$.
+\item To simplify some of the sums, we replace sums over $j<i$ with double sums over
+all particles (remember, $i \neq j$) and divide by 2.
+\end{itemize}
+
+Starting from the above definition of the virial, we then get
 \newcommand{\di}{\delta_{i}}
 \newcommand{\qrt}{\frac{1}{4}}
 \bea
-\Xi    
-&~=~&-\half~\sum_{i < j}^{N}~\rnij\otimes\Fvij                         \\
-&~=~&-\qrt\sum_{i=1}^N~\sum_{j=1}^N ~(\rvi+\di-\rvj)\otimes\Fvij       \\
-&~=~&-\qrt\sum_{i=1}^N~\sum_{j=1}^N ~(\rvi+\di)\otimes\Fvij-\rvj\otimes\Fvij   \\
-&~=~&-\qrt\left(\sum_{i=1}^N~\sum_{j=1}^N ~(\rvi+\di)\otimes\Fvij~-~\sum_{i=1}^N~\sum_{j=1}^N ~\rvj\otimes\Fvij\right) \\
-&~=~&-\qrt\left(\sum_{i=1}^N~(\rvi+\di)\otimes\sum_{j=1}^N~\Fvij~-~\sum_{j=1}^N ~\rvj\otimes\sum_{i=1}^N~\Fvij\right)  \\
-&~=~&-\qrt\left(\sum_{i=1}^N~(\rvi+\di)\otimes\Fvi~+~\sum_{j=1}^N ~\rvj\otimes\Fvj\right)      \\
-&~=~&-\qrt\left(2~\sum_{i=1}^N~\rvi\otimes\Fvi+\sum_{i=1}^N~\di\otimes\Fvi\right)
+\Xi
+&~=~&-\half~\sum_{i < j}^{N}~{\mathbf r}^n_{ij} \otimes {\mathbf F}_{ij} \nonumber \\
+&~=~&-\half~\sum_{i < j}^{N}~\left( {\mathbf r}_i + \delta_{ij} - {\mathbf r}_j \right) \otimes {\mathbf F}_{ij} \nonumber \\
+&~=~&-\qrt~\sum_{i=1}^{N}~\sum_{j=1}^{N}~\left( {\mathbf r}_i + \delta_{ij} - {\mathbf r}_j \right) \otimes {\mathbf F}_{ij} \nonumber \\
+&~=~&-\qrt~\sum_{i=1}^{N}~\sum_{s}~\sum_{j=1}^{N}~\left( {\mathbf r}_i + \delta_{i,s} - {\mathbf r}_j \right) \otimes {\mathbf F}_{ij,s} \nonumber \\
+&~=~&-\qrt~\sum_{i=}^{N}~\sum_{s}~\sum_{j=1}^{N}~\left( \left( {\mathbf r}_i + \delta_{i,s} \right) \otimes {\mathbf F}_{ij,s} -{\mathbf r}_j \otimes {\mathbf F}_{ij,s} \right) \nonumber \\
+&~=~&-\qrt~\sum_{i=1}^{N}~\sum_{s}~\sum_{j=1}^N ~\left( {\mathbf r}_i + \delta_{i,s} \right) \otimes {\mathbf F}_{ij,s} + \qrt \sum_{i=1}^{N}~\sum_{s}~\sum_{j=1}^{N} {\mathbf r}_j \otimes {\mathbf F}_{ij,s} \nonumber \\
+&~=~&-\qrt~\sum_{i=1}^{N}~\sum_{s}~\sum_{j=1}^N ~\left( {\mathbf r}_i + \delta_{i,s} \right) \otimes {\mathbf F}_{ij,s} + \qrt \sum_{i=1}^{N}~\sum_{j=1}^{N} {\mathbf r}_j \otimes {\mathbf F}_{ij} \nonumber \\
+&~=~&-\qrt~\sum_{s}~\sum_{i=1}^{N}~\left( {\mathbf r}_i + \delta_{i,s} \right) \otimes ~\sum_{j=1}^N {\mathbf F}_{ij,s} + \qrt \sum_{j=1}^N {\mathbf r}_j \otimes \sum_{i=1}^{N} {\mathbf F}_{ij} \nonumber \\
+&~=~&-\qrt~\sum_{s}~\sum_{i=1}^{N}~\left( {\mathbf r}_i + \delta_{i,s} \right) \otimes ~\sum_{j=1}^N {\mathbf F}_{ij,s} - \qrt \sum_{j=1}^N {\mathbf r}_j \otimes \sum_{i=1}^{N} {\mathbf F}_{ji} \nonumber \\
+&~=~&-\qrt~\sum_{s}~\sum_{i=1}^{N}~\left( {\mathbf r}_i + \delta_{i,s} \right) \otimes {\mathbf F}_{i,s} - \qrt \sum_{j=1}^N~{\mathbf r}_j \otimes {\mathbf F}_{j}  \nonumber \\
+&~=~&-\qrt~\left(\sum_{i=1}^{N}~{\mathbf r}_i  \otimes {\mathbf F}_{i} + \sum_{j=1}^N~{\mathbf r}_j \otimes {\mathbf F}_{j} \right) - \qrt \sum_{s}~\sum_{i=1}^{N} \delta_{i,s} \otimes {\mathbf F}_{i,s}  \nonumber \\
+&~=~&-\half \sum_{i=1}^{N}~{\mathbf r}_i \otimes {\mathbf F}_{i} -\qrt \sum_{s}~\sum_{i=1}^{N}~\delta_{i,s} \otimes {\mathbf F}_{i,s} \nonumber \\
+&~=~&-\half \sum_{i=1}^{N}~{\mathbf r}_i \otimes {\mathbf F}_{i} -\qrt \sum_{s}~\delta_{s} \otimes {\mathbf F}_{s} \nonumber \\
+&~=~&\Xi_0 + \Xi_1
 \eea
-In these formulae we introduced:
+In the second-last stage, we have used the property that each shift vector itself does not depend on the coordinates of particle $i$, so it is possible to sum
+up all forces corresponding to each shift vector (in the nonbonded kernels), and then just use a sum over the different shift vectors outside the kernels.
+We have also used
 \bea
 \Fvi&~=~&\sum_{j=1}^N~\Fvij                                    \\
 \Fvj&~=~&\sum_{i=1}^N~\Fvji
@@ -97,7 +132,11 @@ which is the total force on $i$ with respect to $j$. Because we use Newton's Thi
 \beq
 \Fvij~=~-\Fvji
 \eeq
-we must, in the implementation, double the term containing the shift $\delta_i$.
+we must, in the implementation, double the term containing the shift $\delta_i$. Similarly, in a few places we have summed the shift-dependent force
+over all shifts to come up with the total force per interaction or particle.
+
+This separates the total virial $\Xi$ into a component $\Xi_0$ that is a single sum over particles, and a second component $\Xi_1$ that describes the influence of
+the particle shifts, and that is only a sum over the different shift vectors.
 
 \subsection{The intra-molecular shift (mol-shift)}
 For the bonded forces and SHAKE it is possible to make a {\em mol-shift}
index 661963a182a6bd85065c83033f5ac7d691f7e052..90f844faa7b8471fd1437f6812581c03aac6807e 100644 (file)
   pages =       {1727--1739}
 }
 
+@Article{Liu2015,
+  author =       {Gang Liu},
+  title =        {Dynamical Equations For The Period Vectors In A Periodic System Under Constant External Stress},
+  journal =      {Can. J. Phys.},
+  year =         2015,
+  volume =       93,
+  pages =        {974--978}
+}
+
 @Article{Loof92,
   author =       "Hans de Loof and Lennart Nilsson and Rudolf Rigler",
   title =        "Molecular Dynamics Simulations of Galanin in Aqueous