docs/manual/implement.tex

   1 %
   2 % This file is part of the GROMACS molecular simulation package.
   3 %
   4 % Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
   5 % Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6 % and including many others, as listed in the AUTHORS file in the
   7 % top-level source directory and at http://www.gromacs.org.
   8 %
   9 % GROMACS is free software; you can redistribute it and/or
  10 % modify it under the terms of the GNU Lesser General Public License
  11 % as published by the Free Software Foundation; either version 2.1
  12 % of the License, or (at your option) any later version.
  13 %
  14 % GROMACS is distributed in the hope that it will be useful,
  15 % but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17 % Lesser General Public License for more details.
  18 %
  19 % You should have received a copy of the GNU Lesser General Public
  20 % License along with GROMACS; if not, see
  21 % http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22 % Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23 %
  24 % If you want to redistribute modifications to GROMACS, please
  25 % consider that scientific software is very special. Version
  26 % control is crucial - bugs must be traceable. We will be happy to
  27 % consider code for inclusion in the official distribution, but
  28 % derived work must not be called official GROMACS. Details are found
  29 % in the README & COPYING files - if they are missing, get the
  30 % official version at http://www.gromacs.org.
  31 %
  32 % To help us fund GROMACS development, we humbly ask that you cite
  33 % the research papers on the package. Check out http://www.gromacs.org.
  34
  35 \chapter{Some implementation details}
  36 In this chapter we will present some implementation details. This is
  37 far from complete, but we deemed it necessary to clarify some things
  38 that would otherwise be hard to understand.
  39
  40 \section{Single Sum Virial in {\gromacs}}
  41 \label{sec:virial}
  42 The \normindex{virial} $\Xi$ can be written in full tensor form as:
  43 \beq
  44 \Xi~=~-\half~\sum_{i < j}^N~\rvij\otimes\Fvij
  45 \eeq
  46 where $\otimes$ denotes the {\em direct product} of two vectors.\footnote
  47 {$({\bf u}\otimes{\bf v})^{\ab}~=~{\bf u}_{\al}{\bf v}_{\be}$} When this is
  48 computed in the inner loop of an MD program 9 multiplications and 9
  49 additions are needed.\footnote{The calculation of
  50 Lennard-Jones and Coulomb forces is about 50 floating point operations.}
  51
  52 Here it is shown how it is possible to extract the virial calculation
  53 from the inner loop~\cite{Bekker93b}.
  54
  55 \subsection{Virial}
  56 In a system with periodic boundary conditions\index{periodic boundary
  57 conditions}, the
  58 periodicity must be taken into account for the virial:
  59 \beq
  60 \Xi~=~-\half~\sum_{i < j}^{N}~\rnij\otimes\Fvij
  61 \eeq
  62 where $\rnij$ denotes the distance vector of the
  63 {\em nearest image} of atom $i$ from atom $j$. In this definition we add
  64 a {\em shift vector} $\delta_i$ to the position vector $\rvi$
  65 of atom $i$. The difference vector $\rnij$ is thus equal to:
  66 \beq
  67 \rnij~=~\rvi+\delta_i-\rvj
  68 \eeq
  69 or in shorthand:
  70 \beq
  71 \rnij~=~\rni-\rvj
  72 \eeq
  73 In a triclinic system, there are 27 possible images of $i$; when a truncated
  74 octahedron is used, there are 15 possible images.
  75
  76 \subsection{Virial from non-bonded forces}
  77 Here the derivation for the single sum virial in the {\em non-bonded force}
  78 routine is given. $i \neq j$ in all formulae below.
  79 \newcommand{\di}{\delta_{i}}
  80 \newcommand{\qrt}{\frac{1}{4}}
  81 \bea
  82 \Xi
  83 &~=~&-\half~\sum_{i < j}^{N}~\rnij\otimes\Fvij                          \\
  84 &~=~&-\qrt\sum_{i=1}^N~\sum_{j=1}^N ~(\rvi+\di-\rvj)\otimes\Fvij        \\
  85 &~=~&-\qrt\sum_{i=1}^N~\sum_{j=1}^N ~(\rvi+\di)\otimes\Fvij-\rvj\otimes\Fvij    \\
  86 &~=~&-\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)  \\
  87 &~=~&-\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)   \\
  88 &~=~&-\qrt\left(\sum_{i=1}^N~(\rvi+\di)\otimes\Fvi~+~\sum_{j=1}^N ~\rvj\otimes\Fvj\right)       \\
  89 &~=~&-\qrt\left(2~\sum_{i=1}^N~\rvi\otimes\Fvi+\sum_{i=1}^N~\di\otimes\Fvi\right)
  90 \eea
  91 In these formulae we introduced:
  92 \bea
  93 \Fvi&~=~&\sum_{j=1}^N~\Fvij                                     \\
  94 \Fvj&~=~&\sum_{i=1}^N~\Fvji
  95 \eea
  96 which is the total force on $i$ with respect to $j$. Because we use Newton's Third Law:
  97 \beq
  98 \Fvij~=~-\Fvji
  99 \eeq
 100 we must, in the implementation, double the term containing the shift $\delta_i$.
 101
 102 \subsection{The intra-molecular shift (mol-shift)}
 103 For the bonded forces and SHAKE it is possible to make a {\em mol-shift}
 104 list, in which the periodicity is stored. We simple have an array {\tt mshift}
 105 in which for each atom an index in the {\tt shiftvec} array is stored.
 106
 107 The algorithm to generate such a list can be derived from graph theory,
 108 considering each particle in a molecule as a bead in a graph, the bonds
 109 as edges.
 110 \begin{enumerate}
 111 \item[1]        Represent the bonds and atoms as bidirectional graph
 112 \item[2]        Make all atoms white
 113 \item[3]        Make one of the white atoms black (atom $i$) and put it in the
 114                 central box
 115 \item[4]        Make all of the neighbors of $i$ that are currently
 116                 white, gray
 117 \item[5]        Pick one of the gray atoms (atom $j$), give it the
 118                 correct periodicity with respect to any of
 119                 its black neighbors
 120                 and make it black
 121 \item[6]        Make all of the neighbors of $j$ that are currently
 122                 white, gray
 123 \item[7]        If any gray atom remains, go to [5]
 124 \item[8]        If any white atom remains, go to [3]
 125 \end{enumerate}
 126 Using this algorithm we can
 127 \begin{itemize}
 128 \item   optimize the bonded force calculation as well as SHAKE
 129 \item   calculate the virial from the bonded forces
 130         in the single sum method again
 131 \end{itemize}
 132
 133 Find a representation of the bonds as a bidirectional graph.
 134
 135 \subsection{Virial from Covalent Bonds}
 136 Since the covalent bond force gives a contribution to the virial, we have:
 137 \bea
 138 b       &~=~&   \|\rnij\|                                       \\
 139 V_b     &~=~&   \half k_b(b-b_0)^2                              \\
 140 \Fvi    &~=~&   -\nabla V_b                                     \\
 141         &~=~&   k_b(b-b_0)\frac{\rnij}{b}                       \\
 142 \Fvj    &~=~&   -\Fvi
 143 \eea
 144 The virial contribution from the bonds then is:
 145 \bea
 146 \Xi_b   &~=~&   -\half(\rni\otimes\Fvi~+~\rvj\otimes\Fvj)       \\
 147         &~=~&   -\half\rnij\otimes\Fvi
 148 \eea
 149
 150 \subsection{Virial from SHAKE}
 151 An important contribution to the virial comes from shake. Satisfying
 152 the constraints a force {\bf G} that is exerted on the particles ``shaken.'' If this
 153 force does not come out of the algorithm (as in standard SHAKE) it can be
 154 calculated afterward (when using {\em leap-frog}) by:
 155 \bea
 156 \Delta\rvi&~=~&\rvi(t+\Dt)-
 157 [\rvi(t)+{\bf v}_i(t-\frac{\Dt}{2})\Dt+\frac{\Fvi}{m_i}\Dt^2]   \\
 158 {\bf G}_i&~=~&\frac{m_i\Delta\rvi}{\Dt^2}
 159 \eea
 160 This does not help us in the general case. Only when no periodicity
 161 is needed (like in rigid water) this can be used, otherwise
 162 we must add the virial calculation in the inner loop of SHAKE.
 163
 164 When it {\em is} applicable the virial can be calculated in the single sum way:
 165 \beq
 166 \Xi~=~-\half\sum_i^{N_c}~\rvi\otimes\Fvi
 167 \eeq
 168 where $N_c$ is the number of constrained atoms.
 169
 170 %Another method is the Non-Iterative shake as proposed (and implemented)
 171 %by Yoneya. In this algorithm the Lagrangian multipliers are solved in a
 172 %matrix equation, and given these multipliers it is easy to get the periodicity
 173 %in the virial afterwards.
 174
 175 %More...
 176
 177
 178 \section{Optimizations}
 179 Here we describe some of the algorithmic optimizations used
 180 in {\gromacs}, apart from parallelism.
 181
 182 \subsection{Inner Loops for Water}
 183 \label{sec:waterloops}
 184 {\gromacs} uses special inner loops to calculate non-bonded
 185 interactions for water molecules with other atoms, and yet
 186 another set of loops for interactions between pairs of
 187 water molecules. There highly optimized loops for two types of water models.
 188 For three site models similar to
 189 SPC~\cite{Berendsen81}, {\ie}:
 190 \begin{enumerate}
 191 \item   There are three atoms in the molecule.
 192 \item   The whole molecule is a single charge group.
 193 \item   The first atom has Lennard-Jones (\secref{lj}) and
 194         Coulomb (\secref{coul}) interactions.
 195 \item   Atoms two and three have only Coulomb interactions,
 196         and equal charges.
 197 \end{enumerate}
 198 These loops also works for the SPC/E~\cite{Berendsen87} and
 199 TIP3P~\cite{Jorgensen83} water models.
 200 And for four site water models similar to TIP4P~\cite{Jorgensen83}:
 201 \begin{enumerate}
 202 \item   There are four atoms in the molecule.
 203 \item   The whole molecule is a single charge group.
 204 \item   The first atom has only Lennard-Jones (\secref{lj}) interactions.
 205 \item   Atoms two and three have only Coulomb (\secref{coul}) interactions,
 206         and equal charges.
 207 \item   Atom four has only Coulomb interactions.
 208 \end{enumerate}
 209
 210 The benefit of these implementations is that there are more floating-point
 211 operations in a single loop, which implies that some compilers
 212 can schedule the code better. However, it turns out that even
 213 some of the most advanced compilers have problems with scheduling,
 214 implying that manual tweaking is necessary to get optimum
 215 \normindex{performance}.
 216 This may include common-sub-expression elimination, or moving
 217 code around.
 218
 219 % LocalWords:  Virial virial triclinic intra mol mshift shiftvec sqrt SPC lj yf
 220 % LocalWords:  coul