docs/reference-manual/details.rst

   1 Some implementation details
   2 ===========================
   3
   4 In this chapter we will present some implementation details. This is far
   5 from complete, but we deemed it necessary to clarify some things that
   6 would otherwise be hard to understand.
   7
   8 Single Sum Virial in |Gromacs|
   9 ------------------------------
  10
  11 The virial :math:`\Xi` can be written in full tensor form as:
  12
  13 .. math:: \Xi~=~-\frac{1}{2}~\sum_{i < j}^N~\mathbf{r}_{ij}\otimes\mathbf{F}_{ij}
  14           :label: eqnvirialfulltensorform
  15
  16 where :math:`\otimes` denotes the *direct product* of two vectors. [1]_
  17 When this is computed in the inner loop of an MD program 9
  18 multiplications and 9 additions are needed. [2]_
  19
  20 Here it is shown how it is possible to extract the virial calculation
  21 from the inner loop \ :ref:`177 <refBekker93b>`.
  22
  23 Virial
  24 ~~~~~~
  25
  26 In a system with periodic boundary conditions, the periodicity must be
  27 taken into account for the virial:
  28
  29 .. math:: \Xi~=~-\frac{1}{2}~\sum_{i < j}^{N}~\mathbf{r}_{ij}^n\otimes\mathbf{F}_{ij}
  30           :label: eqnvirialperiodicity
  31
  32 where :math:`\mathbf{r}_{ij}^n` denotes the distance
  33 vector of the *nearest image* of atom :math:`i` from atom :math:`j`. In
  34 this definition we add a *shift vector* :math:`\delta_i` to the position
  35 vector :math:`\mathbf{r}_i` of atom :math:`i`. The
  36 difference vector :math:`\mathbf{r}_{ij}^n` is thus equal
  37 to:
  38
  39 .. math:: \mathbf{r}_{ij}^n~=~\mathbf{r}_i+\delta_i-\mathbf{r}_j
  40           :label: eqnvirialdiffvector
  41
  42 or in shorthand:
  43
  44 .. math:: \mathbf{r}_{ij}^n~=~\mathbf{r}_i^n-\mathbf{r}_j
  45           :label: eqnvirialdiffvecshort
  46
  47 In a triclinic system, there are 27 possible images of :math:`i`; when
  48 a truncated octahedron is used, there are 15 possible images.
  49
  50 Virial from non-bonded forces
  51 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  52
  53 Here the derivation for the single sum virial in the *non-bonded force*
  54 routine is given. There are a couple of considerations that are special
  55 to |Gromacs| that we take into account:
  56
  57 -  When calculating short-range interactions, we apply the *minimum
  58    image convention* and only consider the closest image of each
  59    neighbor - and in particular we never allow interactions between a
  60    particle and any of its periodic images. For all the equations below,
  61    this means :math:`i \neq j`.
  62
  63 -  In general, either the :math:`i` or :math:`j` particle might be
  64    shifted to a neighbor cell to get the closest interaction (shift
  65    :math:`\delta_{ij}`). However, with minimum image convention there
  66    can be at most 27 different shifts for particles in the central cell,
  67    and for typical (very short-ranged) biomolecular interactions there
  68    are typically only a few different shifts involved for each particle,
  69    not to mention that each interaction can only be present for one
  70    shift.
  71
  72 -  For the |Gromacs| nonbonded interactions we use this to split the
  73    neighborlist of each :math:`i` particle into multiple separate lists,
  74    where each list has a constant shift :math:`\delta_i` for the
  75    :math:`i` partlcle. We can represent this as a sum over shifts (for
  76    which we use index :math:`s`), with the constraint that each particle
  77    interaction can only contribute to one of the terms in this sum, and
  78    the shift is no longer dependent on the :math:`j` particles. For any
  79    sum that does not contain complex dependence on :math:`s`, this means
  80    the sum trivially reduces to just the sum over :math:`i` and/or
  81    :math:`j`.
  82
  83 -  To simplify some of the sums, we replace sums over :math:`j<i` with
  84    double sums over all particles (remember, :math:`i \neq j`) and
  85    divide by 2.
  86
  87 Starting from the above definition of the virial, we then get
  88
  89 .. math:: \begin{aligned}
  90           \Xi
  91           &~=~&-{\frac{1}{2}}~\sum_{i < j}^{N}~{\mathbf r}^n_{ij} \otimes {\mathbf F}_{ij} \nonumber \\
  92           &~=~&-{\frac{1}{2}}~\sum_{i < j}^{N}~\left( {\mathbf r}_i + \delta_{ij} - {\mathbf r}_j \right) \otimes {\mathbf F}_{ij} \nonumber \\
  93           &~=~&-{\frac{1}{4}}~\sum_{i=1}^{N}~\sum_{j=1}^{N}~\left( {\mathbf r}_i + \delta_{ij} - {\mathbf r}_j \right) \otimes {\mathbf F}_{ij} \nonumber \\
  94           &~=~&-{\frac{1}{4}}~\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 \\
  95           &~=~&-{\frac{1}{4}}~\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 \\
  96           &~=~&-{\frac{1}{4}}~\sum_{i=1}^{N}~\sum_{s}~\sum_{j=1}^N ~\left( {\mathbf r}_i + \delta_{i,s} \right) \otimes {\mathbf F}_{ij,s} + {\frac{1}{4}}\sum_{i=1}^{N}~\sum_{s}~\sum_{j=1}^{N} {\mathbf r}_j \otimes {\mathbf F}_{ij,s} \nonumber \\
  97           &~=~&-{\frac{1}{4}}~\sum_{i=1}^{N}~\sum_{s}~\sum_{j=1}^N ~\left( {\mathbf r}_i + \delta_{i,s} \right) \otimes {\mathbf F}_{ij,s} + {\frac{1}{4}}\sum_{i=1}^{N}~\sum_{j=1}^{N} {\mathbf r}_j \otimes {\mathbf F}_{ij} \nonumber \\
  98           &~=~&-{\frac{1}{4}}~\sum_{s}~\sum_{i=1}^{N}~\left( {\mathbf r}_i + \delta_{i,s} \right) \otimes ~\sum_{j=1}^N {\mathbf F}_{ij,s} + {\frac{1}{4}}\sum_{j=1}^N {\mathbf r}_j \otimes \sum_{i=1}^{N} {\mathbf F}_{ij} \nonumber \\
  99           &~=~&-{\frac{1}{4}}~\sum_{s}~\sum_{i=1}^{N}~\left( {\mathbf r}_i + \delta_{i,s} \right) \otimes ~\sum_{j=1}^N {\mathbf F}_{ij,s} - {\frac{1}{4}}\sum_{j=1}^N {\mathbf r}_j \otimes \sum_{i=1}^{N} {\mathbf F}_{ji} \nonumber \\
 100           &~=~&-{\frac{1}{4}}~\sum_{s}~\sum_{i=1}^{N}~\left( {\mathbf r}_i + \delta_{i,s} \right) \otimes {\mathbf F}_{i,s} - {\frac{1}{4}}\sum_{j=1}^N~{\mathbf r}_j \otimes {\mathbf F}_{j}  \nonumber \\
 101           &~=~&-{\frac{1}{4}}~\left(\sum_{i=1}^{N}~{\mathbf r}_i  \otimes {\mathbf F}_{i} + \sum_{j=1}^N~{\mathbf r}_j \otimes {\mathbf F}_{j} \right) - {\frac{1}{4}}\sum_{s}~\sum_{i=1}^{N} \delta_{i,s} \otimes {\mathbf F}_{i,s}  \nonumber \\
 102           &~=~&-{\frac{1}{2}}\sum_{i=1}^{N}~{\mathbf r}_i \otimes {\mathbf F}_{i} -{\frac{1}{4}}\sum_{s}~\sum_{i=1}^{N}~\delta_{i,s} \otimes {\mathbf F}_{i,s} \nonumber \\
 103           &~=~&-{\frac{1}{2}}\sum_{i=1}^{N}~{\mathbf r}_i \otimes {\mathbf F}_{i} -{\frac{1}{4}}\sum_{s}~\delta_{s} \otimes {\mathbf F}_{s} \nonumber \\
 104           &~=~&\Xi_0 + \Xi_1\end{aligned}
 105           :label: eqnviriallong
 106
 107 In the second-last stage, we have used the property that each shift
 108 vector itself does not depend on the coordinates of particle :math:`i`,
 109 so it is possible to sum up all forces corresponding to each shift
 110 vector (in the nonbonded kernels), and then just use a sum over the
 111 different shift vectors outside the kernels. We have also used
 112
 113 .. math:: \begin{aligned}
 114           \mathbf{F}_i&~=~&\sum_{j=1}^N~\mathbf{F}_{ij}                                 \\
 115           \mathbf{F}_j&~=~&\sum_{i=1}^N~\mathbf{F}_{ji}\end{aligned}
 116           :label: eqnvirialtotalforce
 117
 118 which is the total force on :math:`i` with respect to :math:`j`.
 119 Because we use Newton’s Third Law:
 120
 121 .. math:: \mathbf{F}_{ij}~=~-\mathbf{F}_{ji}
 122           :label: eqnnewtonsthird
 123
 124 we must, in the implementation, double the term containing the shift
 125 :math:`\delta_i`. Similarly, in a few places we have summed the
 126 shift-dependent force over all shifts to come up with the total force
 127 per interaction or particle.
 128
 129 This separates the total virial :math:`\Xi` into a component
 130 :math:`\Xi_0` that is a single sum over particles, and a second
 131 component :math:`\Xi_1` that describes the influence of the particle
 132 shifts, and that is only a sum over the different shift vectors.
 133
 134 The intra-molecular shift (mol-shift)
 135 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 136
 137 For the bonded forces and SHAKE it is possible to make a *mol-shift*
 138 list, in which the periodicity is stored. We simple have an array mshift
 139 in which for each atom an index in the shiftvec array is stored.
 140
 141 The algorithm to generate such a list can be derived from graph theory,
 142 considering each particle in a molecule as a bead in a graph, the bonds
 143 as edges.
 144
 145 #. Represent the bonds and atoms as bidirectional graph
 146
 147 #. Make all atoms white
 148
 149 #. Make one of the white atoms black (atom :math:`i`) and put it in the
 150    central box
 151
 152 #. Make all of the neighbors of :math:`i` that are currently white, gray
 153
 154 #. Pick one of the gray atoms (atom :math:`j`), give it the correct
 155    periodicity with respect to any of its black neighbors and make it
 156    black
 157
 158 #. Make all of the neighbors of :math:`j` that are currently white, gray
 159
 160 #. If any gray atom remains, go to [5]
 161
 162 #. If any white atom remains, go to [3]
 163
 164 Using this algorithm we can
 165
 166 -  optimize the bonded force calculation as well as SHAKE
 167
 168 -  calculate the virial from the bonded forces in the single sum method
 169    again
 170
 171 Find a representation of the bonds as a bidirectional graph.
 172
 173 Virial from Covalent Bonds
 174 ~~~~~~~~~~~~~~~~~~~~~~~~~~
 175
 176 Since the covalent bond force gives a contribution to the virial, we
 177 have:
 178
 179 .. math:: \begin{aligned}
 180           b     &~=~&   \|\mathbf{r}_{ij}^n\|                                   \\
 181           V_b   &~=~&   \frac{1}{2} k_b(b-b_0)^2                                \\
 182           \mathbf{F}_i  &~=~&   -\nabla V_b                                     \\
 183                 &~=~&   k_b(b-b_0)\frac{\mathbf{r}_{ij}^n}{b}                   \\
 184           \mathbf{F}_j  &~=~&   -\mathbf{F}_i\end{aligned}
 185           :label: eqncovbondvirial
 186
 187 The virial contribution from the bonds then is:
 188
 189 .. math:: \begin{aligned}
 190           \Xi_b &~=~&   -\frac{1}{2}(\mathbf{r}_i^n\otimes\mathbf{F}_i~+~\mathbf{r}_j\otimes\mathbf{F}_j)       \\
 191           &~=~& -\frac{1}{2}\mathbf{r}_{ij}^n\otimes\mathbf{F}_i\end{aligned}
 192           :label: eqncovbondvirialcontri
 193
 194 Virial from SHAKE
 195 ~~~~~~~~~~~~~~~~~
 196
 197 An important contribution to the virial comes from shake. Satisfying the
 198 constraints a force **G** that is exerted on the particles “shaken.” If
 199 this force does not come out of the algorithm (as in standard SHAKE) it
 200 can be calculated afterward (when using *leap-frog*) by:
 201
 202 .. math:: \begin{aligned}
 203           \Delta\mathbf{r}_i&~=~&{\mathbf{r}_i}(t+{\Delta t})-
 204           [\mathbf{r}_i(t)+{\bf v}_i(t-\frac{{\Delta t}}{2}){\Delta t}+\frac{\mathbf{F}_i}{m_i}{\Delta t}^2]    \\
 205           {\bf G}_i&~=~&\frac{m_i{\Delta}{\mathbf{r}_i}}{{\Delta t}^2i}\end{aligned}
 206           :label: eqnshakevirial
 207
 208 This does not help us in the general case. Only when no periodicity is
 209 needed (like in rigid water) this can be used, otherwise we must add the
 210 virial calculation in the inner loop of SHAKE.
 211
 212 When it *is* applicable the virial can be calculated in the single sum
 213 way:
 214
 215 .. math:: \Xi~=~-\frac{1}{2}\sum_i^{N_c}~\mathbf{r}_i\otimes\mathbf{F}_i
 216           :label: eqnshakevirialsinglesum
 217
 218 where :math:`N_c` is the number of constrained atoms.
 219
 220 Optimizations
 221 -------------
 222
 223 Here we describe some of the algorithmic optimizations used in |Gromacs|,
 224 apart from parallelism.
 225
 226 .. _waterloops:
 227
 228 Inner Loops for Water
 229 ~~~~~~~~~~~~~~~~~~~~~
 230
 231 |Gromacs| uses special inner loops to calculate non-bonded interactions
 232 for water molecules with other atoms, and yet another set of loops for
 233 interactions between pairs of water molecules. There highly optimized
 234 loops for two types of water models. For three site models similar to
 235 SPC \ :ref:`80 <refBerendsen81>`, *i.e.*:
 236
 237 #. There are three atoms in the molecule.
 238
 239 #. The whole molecule is a single charge group.
 240
 241 #. The first atom has Lennard-Jones (sec. :ref:`lj`) and Coulomb
 242    (sec. :ref:`coul`) interactions.
 243
 244 #. Atoms two and three have only Coulomb interactions, and equal
 245    charges.
 246
 247 These loops also works for the SPC/E \ :ref:`178 <refBerendsen87>` and
 248 TIP3P \ :ref:`128 <refJorgensen83>` water models. And for four site water
 249 models similar to TIP4P \ :ref:`128 <refJorgensen83>`:
 250
 251 #. There are four atoms in the molecule.
 252
 253 #. The whole molecule is a single charge group.
 254
 255 #. The first atom has only Lennard-Jones (sec. :ref:`lj`) interactions.
 256
 257 #. Atoms two and three have only Coulomb (sec. :ref:`coul`) interactions,
 258    and equal charges.
 259
 260 #. Atom four has only Coulomb interactions.
 261
 262 The benefit of these implementations is that there are more
 263 floating-point operations in a single loop, which implies that some
 264 compilers can schedule the code better. However, it turns out that even
 265 some of the most advanced compilers have problems with scheduling,
 266 implying that manual tweaking is necessary to get optimum performance.
 267 This may include common-sub-expression elimination, or moving code
 268 around.
 269
 270 .. raw:: latex
 271
 272     \clearpage
 273
 274
 275 .. [1]
 276    Note that some derivations, an alternative notation
 277    :math:`\xi_{\mathrm{alt}} = v_{\xi} = p_{\xi}/Q` is used.
 278
 279 .. [2]
 280    The calculation of Lennard-Jones and Coulomb forces is about 50
 281    floating point operations.