1Methods
2-------
3
4Exclusions and 1-4 Interactions.
5~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
6
7Atoms within a molecule that are close by in the chain, *i.e.* atoms
8that are covalently bonded, or linked by one or two atoms are called
9*first neighbors, second neighbors* and *third neighbors*, respectively
10(see :numref:`Fig. %s <fig-chain>`). Since the interactions of atom **i** with atoms
11**i+1** and **i+2** are mainly quantum mechanical, they can not be
12modeled by a Lennard-Jones potential. Instead it is assumed that these
13interactions are adequately modeled by a harmonic bond term or
14constraint (**i, i+1**) and a harmonic angle term (**i, i+2**). The
15first and second neighbors (atoms **i+1** and **i+2**) are therefore
16*excluded* from the Lennard-Jones interaction list of atom **i**; atoms
17**i+1** and **i+2** are called *exclusions* of atom **i**.
18
19.. _fig-chain:
20
21.. figure:: plots/chain.*
22   :width: 8.00000cm
23
24   Atoms along an alkane chain.
25
26For third neighbors, the normal Lennard-Jones repulsion is sometimes
27still too strong, which means that when applied to a molecule, the
28molecule would deform or break due to the internal strain. This is
29especially the case for carbon-carbon interactions in a
30*cis*-conformation (*e.g.* *cis*-butane). Therefore, for some of these
31interactions, the Lennard-Jones repulsion has been reduced in the GROMOS
32force field, which is implemented by keeping a separate list of 1-4 and
33normal Lennard-Jones parameters. In other force fields, such as
34OPLS \ :ref:`103 <refJorgensen88>`, the standard Lennard-Jones
35parameters are reduced by a factor of two, but in that case also the
36dispersion (r\ :math:`^{-6}`) and the Coulomb interaction are scaled.
37|Gromacs| can use either of these methods.
38
39Charge Groups
40~~~~~~~~~~~~~
41
42In principle, the force calculation in MD is an :math:`O(N^2)` problem.
43Therefore, we apply a cut-off for non-bonded force (NBF) calculations;
44only the particles within a certain distance of each other are
45interacting. This reduces the cost to :math:`O(N)` (typically
46:math:`100N` to :math:`200N`) of the NBF. It also introduces an error,
47which is, in most cases, acceptable, except when applying the cut-off
48implies the creation of charges, in which case you should consider using
49the lattice sum methods provided by |Gromacs|.
50
51Consider a water molecule interacting with another atom. If we would
52apply a plain cut-off on an atom-atom basis we might include the
53atom-oxygen interaction (with a charge of :math:`-0.82`) without the
54compensating charge of the protons, and as a result, induce a large
55dipole moment over the system. Therefore, we have to keep groups of
56atoms with total charge 0 together. These groups are called *charge
57groups*. Note that with a proper treatment of long-range electrostatics
58(e.g. particle-mesh Ewald (sec. :ref:`pme`), keeping charge groups
59together is not required.
60
61Treatment of Cut-offs in the group scheme
62~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
63
64|Gromacs| is quite flexible in treating cut-offs, which implies there can
65be quite a number of parameters to set. These parameters are set in the
66input file for grompp. There are two sort of parameters that affect the
67cut-off interactions; you can select which type of interaction to use in
68each case, and which cut-offs should be used in the neighbor searching.
69
70For both Coulomb and van der Waals interactions there are interaction
71type selectors (termed vdwtype and coulombtype) and two parameters, for
72a total of six non-bonded interaction parameters. See the User Guide for
73a complete description of these parameters.
74
75In the group cut-off scheme, all of the interaction functions in
76:numref:`Table %s <tab-funcparm>` require that neighbor searching be done with a
77radius at least as large as the :math:`r_c` specified for the functional
78form, because of the use of charge groups. The extra radius is typically
79of the order of 0.25 nm (roughly the largest distance between two atoms
80in a charge group plus the distance a charge group can diffuse within
81neighbor list updates).
82
83.. |CPCOP| replace:: :math:`r_c`, :math:`{\varepsilon}_{r}`
84.. |CRFP|  replace:: :math:`r_c`, :math:`{\varepsilon}_{rf}`
85.. |CSHFP| replace:: :math:`r_1`, :math:`r_c`, :math:`{\varepsilon}_{r}`
86.. |CSWFP| replace:: :math:`r_1`, :math:`r_c`, :math:`{\varepsilon}_{r}`
87.. |VPCOP| replace:: :math:`r_c`
88.. |VSHFP| replace:: :math:`r_1`, :math:`r_c`
89.. |VSWFP| replace:: :math:`r_1`, :math:`r_c`
90
91.. _tab-funcparm:
92
93.. table:: Parameters for the different functional forms of the
94           non-bonded interactions.
95
96           +----------------------------+------------+
97           | Type                       | Parameters |
98           +=========+==================+============+
99           | Coulomb | Plain cut-off    | |CPCOP|    |
100           |         +------------------+------------+
101           |         | Reaction field   | |CRFP|     |
102           |         +------------------+------------+
103           |         | Shift function   | |CSHFP|    |
104           |         +------------------+------------+
105           |         | Switch function  | |CSWFP|    |
106           +---------+------------------+------------+
107           | VdW     | Plain cut-off    | |VPCOP|    |
108           |         +------------------+------------+
109           |         | Shift function   | |VSHFP|    |
110           |         +------------------+------------+
111           |         | Switch function  | |VSWFP|    |
112           +---------+------------------+------------+
113
114.. _virtualsites:
115
116Virtual interaction sites
117-------------------------
118
119Virtual interaction sites (called dummy atoms in
120|Gromacs| versions before 3.3) can be used in |Gromacs| in a number of ways.
121We write the position of the virtual site :math:`\mathbf{r}_s` as a function
122of the positions of other particles
123:math:`\mathbf{r}`\ :math:`_i`: :math:`\mathbf{r}_s =
124f(\mathbf{r}_1..\mathbf{r}_n)`. The virtual site, which may carry charge or be
125involved in other interactions, can now be used in the force
126calculation. The force acting on the virtual site must be redistributed
127over the particles with mass in a consistent way. A good way to do this
128can be found in ref. \ :ref:`104 <refBerendsen84b>`. We can write the
129potential energy as:
130
131.. math:: V = V(\mathbf{r}_s,\mathbf{r}_1,\ldots,\mathbf{r}_n) = V^*(\mathbf{r}_1,\ldots,\mathbf{r}_n)
132          :label: eqnvsiteepot
133
134The force on the particle :math:`i` is then:
135
136.. math:: \mathbf{F}_i = -\frac{\partial V^*}{\partial \mathbf{r}_i}
137          = -\frac{\partial V}{\partial \mathbf{r}_i} -
138             \frac{\partial V}{\partial \mathbf{r}_s}
139             \frac{\partial \mathbf{r}_s}{\partial \mathbf{r}_i}
140          = \mathbf{F}_i^{direct} + \mathbf{F}_i
141          :label: eqnvsiteforce
142
143The first term is the normal force. The second term is the force on
144particle :math:`i` due to the virtual site, which can be written in
145tensor notation:
146
147.. math::  \mathbf{F}_i = \left[\begin{array}{ccc}
148           {\displaystyle\frac{\partial x_s}{\partial x_i}} & {\displaystyle\frac{\partial y_s}{\partial x_i}} & {\displaystyle\frac{\partial z_s}{\partial x_i}} \\[1ex]
149           {\displaystyle\frac{\partial x_s}{\partial y_i}} & {\displaystyle\frac{\partial y_s}{\partial y_i}} & {\displaystyle\frac{\partial z_s}{\partial y_i}} \\[1ex]
150           {\displaystyle\frac{\partial x_s}{\partial z_i}} & {\displaystyle\frac{\partial y_s}{\partial z_i}} & {\displaystyle\frac{\partial z_s}{\partial z_i}} \end{array}\right]\mathbf{F}_{s}
151           :label: eqnfvsite
152
153where :math:`\mathbf{F}_{s}` is the force on the virtual site and
154:math:`x_s`, :math:`y_s` and :math:`z_s` are the coordinates of the
155virtual site. In this way, the total force and the total torque are
156conserved \ :ref:`104 <refBerendsen84b>`.
157
158The computation of the virial (:eq:`eqn. %s <eqnXi>`) is non-trivial when
159virtual sites are used. Since the virial involves a summation over all
160the atoms (rather than virtual sites), the forces must be redistributed
161from the virtual sites to the atoms (using  :eq:`eqn. %s <eqnfvsite>`) *before*
162computation of the virial. In some special cases where the forces on the
163atoms can be written as a linear combination of the forces on the
164virtual sites (types 2 and 3 below) there is no difference between
165computing the virial before and after the redistribution of forces.
166However, in the general case redistribution should be done first.
167
168.. _fig-vsites:
169
170.. figure:: plots/dummies.*
171   :width: 15.00000cm
172
173   The seven different types of virtual site construction. The
174   constructing atoms are shown as black circles, the virtual sites in
175   gray.
176
177There are six ways to construct virtual sites from surrounding atoms in
178|Gromacs|, which we classify by the number of constructing atoms. **Note**
179that all site types mentioned can be constructed from types 3fd
180(normalized, in-plane) and 3out (non-normalized, out of plane). However,
181the amount of computation involved increases sharply along this list, so
182we strongly recommended using the first adequate virtual site type that
183will be sufficient for a certain purpose. :numref:`Fig. %s <fig-vsites>` depicts 6 of
184the available virtual site constructions. The conceptually simplest
185construction types are linear combinations:
186
187.. math:: \mathbf{r}_s = \sum_{i=1}^N w_i \, \mathbf{r}_i
188          :label: eqnvsitelincomb
189
190The force is then redistributed using the same weights:
191
192.. math:: \mathbf{F}_i = w_i \, \mathbf{F}_{s}
193          :label: eqnvsitelincombforce
194
195The types of virtual sites supported in |Gromacs| are given in the list
196below. Constructing atoms in virtual sites can be virtual sites
197themselves, but only if they are higher in the list, i.e. virtual sites
198can be constructed from “particles” that are simpler virtual sites.
199
200-  On top of an atom. This allows giving an atom multiple atom types and
201   with that also assigned multiple, different bonded interactions. This
202   can espically be of use in free-energy calculations.
203
204-  The coordinates of the virtual site equal that of the constructing atom:
205
206   .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i
207             :label: eqnvsite1
208
209-  The force is moved to the constructing atom:
210
211   .. math:: \mathbf{F}_i ~=~ \mathbf{F}_{s}
212             :label: eqnvsite1force
213
214-  As a linear combination of two atoms
215   (:numref:`Fig. %s <fig-vsites>` 2):
216
217   .. math:: w_i = 1 - a ~,~~ w_j = a
218             :label: eqnvsitelin2atom
219
220-  In this case the virtual site is on the line through atoms :math:`i`
221   and :math:`j`.
222
223-  On the line through two atoms, with a fixed distance
224   (:numref:`Fig. %s <fig-vsites>` 2fd):
225
226   .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + a \frac{ \mathbf{r}_{ij} }
227                                                  { | \mathbf{r}_{ij} | }
228             :label: eqnvsite2fdatom
229
230-  In this case the virtual site is on the line through the other two
231   particles at a distance of :math:`|a|` from :math:`i`. The force on
232   particles :math:`i` and :math:`j` due to the force on the virtual site
233   can be computed as:
234
235   .. math:: \begin{array}{lcr}
236                     \mathbf{F}_i &=& \displaystyle \mathbf{F}_{s} - \gamma ( \mathbf{F}_{is} - \mathbf{p} ) \\[1ex]
237                     \mathbf{F}_j &=& \displaystyle \gamma (\mathbf{F}_{s} - \mathbf{p})      \\[1ex]
238                     \end{array}
239                     ~\mbox{ where }~
240                     \begin{array}{c}
241             \displaystyle \gamma = \frac{a}{ | \mathbf{r}_{ij} | } \\[2ex]
242             \displaystyle \mathbf{p} = \frac{ \mathbf{r}_{is} \cdot \mathbf{F}_{s} }
243                                   { \mathbf{r}_{is} \cdot \mathbf{r}_{is} } \mathbf{r}_{is}
244             \end{array}
245             :label: eqnvsite2fdforce
246
247-  As a linear combination of three atoms
248   (:numref:`Fig. %s <fig-vsites>` 3):
249
250   .. math:: w_i = 1 - a - b ~,~~ w_j = a ~,~~ w_k = b
251             :label: eqnvsitelin3atom
252
253-  In this case the virtual site is in the plane of the other three
254   particles.
255
256-  In the plane of three atoms, with a fixed distance
257   (:numref:`Fig. %s <fig-vsites>` 3fd):
258
259   .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + b \frac{  (1 - a) \mathbf{r}_{ij} + a \mathbf{r}_{jk}  }
260                                                  { | (1 - a) \mathbf{r}_{ij} + a \mathbf{r}_{jk} | }
261             :label: eqnvsiteplane3atom
262
263-  In this case the virtual site is in the plane of the other three
264   particles at a distance of :math:`|b|` from :math:`i`. The force on
265   particles :math:`i`, :math:`j` and :math:`k` due to the force on the
266   virtual site can be computed as:
267
268   .. math:: \begin{array}{lcr}
269                     \mathbf{F}_i &=& \displaystyle \mathbf{F}_{s} - \gamma ( \mathbf{F}_{is} - \mathbf{p} ) \\[1ex]
270                     \mathbf{F}_j &=& \displaystyle (1-a)\gamma (\mathbf{F}_{s} - \mathbf{p})      \\[1ex]
271                     \mathbf{F}_k &=& \displaystyle a \gamma (\mathbf{F}_{s} - \mathbf{p})         \\
272                     \end{array}
273                     ~\mbox{ where }~
274                     \begin{array}{c}
275             \displaystyle \gamma = \frac{b}{ | \mathbf{r}_{ij} + a \mathbf{r}_{jk} | } \\[2ex]
276             \displaystyle \mathbf{p} = \frac{ \mathbf{r}_{is} \cdot \mathbf{F}_{s} }
277                                   { \mathbf{r}_{is} \cdot \mathbf{r}_{is} } \mathbf{r}_{is}
278             \end{array}
279             :label: eqnvsiteplane3atomforce
280
281-  In the plane of three atoms, with a fixed angle and
282   distance (:numref:`Fig. %s <fig-vsites>` 3fad):
283
284   .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i +
285             d \cos \theta \frac{\mathbf{r}_{ij}}{ | \mathbf{r}_{ij} | } +
286             d \sin \theta \frac{\mathbf{r}_\perp}{ | \mathbf{r}_\perp | }
287             ~\mbox{ where }~
288             \mathbf{r}_\perp ~=~ \mathbf{r}_{jk} -
289             \frac{ \mathbf{r}_{ij} \cdot \mathbf{r}_{jk} }
290             { \mathbf{r}_{ij} \cdot \mathbf{r}_{ij} }
291             \mathbf{r}_{ij}
292             :label: eqnvsite2fadF
293
294-  In this case the virtual site is in the plane of the other three
295   particles at a distance of :math:`|d|` from :math:`i` at an angle of
296   :math:`\alpha` with :math:`\mathbf{r}_{ij}`. Atom
297   :math:`k` defines the plane and the direction of the angle. **Note**
298   that in this case :math:`b` and :math:`\alpha` must be specified,
299   instead of :math:`a` and :math:`b` (see also sec. :ref:`vsitetop`).
300   The force on particles :math:`i`, :math:`j` and :math:`k` due to the
301   force on the virtual site can be computed as (with
302   :math:`\mathbf{r}_\perp` as defined in
303   :eq:`eqn. %s <eqnvsite2fadF>`):
304
305   .. math:: \begin{array}{c}
306                     \begin{array}{lclllll}
307                     \mathbf{F}_i &=& \mathbf{F}_{s} &-&
308                             \dfrac{d \cos \theta}{ | \mathbf{r}_{ij} | } \mathbf{F}_1 &+&
309                             \dfrac{d \sin \theta}{ | \mathbf{r}_\perp | } \left(
310                             \dfrac{ \mathbf{r}_{ij} \cdot \mathbf{r}_{jk} }
311                                  { \mathbf{r}_{ij} \cdot \mathbf{r}_{ij} } \mathbf{F}_2     +
312                             \mathbf{F}_3 \right)                                \\[3ex]
313                     \mathbf{F}_j &=& &&
314                             \dfrac{d \cos \theta}{ | \mathbf{r}_{ij} | } \mathbf{F}_1 &-&
315                             \dfrac{d \sin \theta}{ | \mathbf{r}_\perp | } \left(
316                              \mathbf{F}_2 +
317                              \dfrac{ \mathbf{r}_{ij} \cdot \mathbf{r}_{jk} }
318                                     { \mathbf{r}_{ij} \cdot \mathbf{r}_{ij} } \mathbf{F}_2 +
319                             \mathbf{F}_3 \right)                                \\[3ex]
320                     \mathbf{F}_k &=& && &&
321                             \dfrac{d \sin \theta}{ | \mathbf{r}_\perp | } \mathbf{F}_2  \\[3ex]
322                     \end{array}                                             \\[5ex]
323                     ~\mbox{where }~
324                     \mathbf{F}_1 = \mathbf{F}_{s} -
325                               \dfrac{ \mathbf{r}_{ij} \cdot \mathbf{F}_{s} }
326                                     { \mathbf{r}_{ij} \cdot \mathbf{r}_{ij} } \mathbf{r}_{ij}
327                     ~\mbox{, }~
328                     \mathbf{F}_2 = \mathbf{F}_1 -
329                               \dfrac{ \mathbf{r}_\perp \cdot \mathbf{F}_{s} }
330                                     { \mathbf{r}_\perp \cdot \mathbf{r}_\perp } \mathbf{r}_\perp
331                     ~\mbox{and }~
332                     \mathbf{F}_3 = \dfrac{ \mathbf{r}_{ij} \cdot \mathbf{F}_{s} }
333                                      { \mathbf{r}_{ij} \cdot \mathbf{r}_{ij} } \mathbf{r}_\perp
334             \end{array}
335             :label: eqnvsite2fadFforce
336
337-  As a non-linear combination of three atoms, out of
338   plane (:numref:`Fig. %s <fig-vsites>` 3out):
339
340   .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + a \mathbf{r}_{ij} + b \mathbf{r}_{ik} +
341                              c (\mathbf{r}_{ij} \times \mathbf{r}_{ik})
342             :label: eqnvsitenonlin3atom
343
344-  This enables the construction of virtual sites out of the plane of
345   the other atoms. The force on particles :math:`i,j` and :math:`k` due
346   to the force on the virtual site can be computed as:
347
348   .. math:: \begin{array}{lcl}
349             \mathbf{F}_j &=& \left[\begin{array}{ccc}
350              a              &  -c\,z_{ik}   & c\,y_{ik}     \\[0.5ex]
351              c\,z_{ik}      &   a           & -c\,x_{ik}    \\[0.5ex]
352             -c\,y_{ik}      &   c\,x_{ik}   & a
353             \end{array}\right]\mathbf{F}_{s}                                 \\
354             \mathbf{F}_k &=& \left[\begin{array}{ccc}
355              b              &   c\,z_{ij}   & -c\,y_{ij}    \\[0.5ex]
356             -c\,z_{ij}      &   b           & c\,x_{ij}     \\[0.5ex]
357              c\,y_{ij}      &  -c\,x_{ij}   & b
358             \end{array}\right]\mathbf{F}_{s}                                 \\
359             \mathbf{F}_i &=& \mathbf{F}_{s} - \mathbf{F}_j - \mathbf{F}_k
360             \end{array}
361             :label: eqnvsitenonlin3atomforce
362
363-  From four atoms, with a fixed distance, see
364   separate :numref:`Fig. %s <fig-vsite4fdn>`. This construction is a bit complex,
365   in particular since the previous type (4fd) could be unstable which
366   forced us to introduce a more elaborate construction:
367
368.. _fig-vsite4fdn:
369
370.. figure:: plots/vsite-4fdn.*
371      :width: 5.00000cm
372
373      The new 4fdn virtual site construction, which is stable even when
374      all constructing atoms are in the same plane.
375
376-
377      .. math::   \begin{aligned}
378                  \mathbf{r}_{ja} &=& a\, \mathbf{r}_{ik} - \mathbf{r}_{ij} = a\, (\mathbf{x}_k - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\
379                  \mathbf{r}_{jb} &=& b\, \mathbf{r}_{il} - \mathbf{r}_{ij} = b\, (\mathbf{x}_l - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\
380                  \mathbf{r}_m &=& \mathbf{r}_{ja} \times \mathbf{r}_{jb} \nonumber \\
381                  \mathbf{x}_s &=& \mathbf{x}_i + c \frac{\mathbf{r}_m}{ | \mathbf{r}_m | }
382                  \end{aligned}
383                  :label: eqnvsite
384
385-  In this case the virtual site is at a distance of :math:`|c|` from
386   :math:`i`, while :math:`a` and :math:`b` are parameters. **Note**
387   that the vectors :math:`\mathbf{r}_{ik}` and :math:`\mathbf{r}_{ij}`
388   are not normalized to save floating-point operations. The force on
389   particles :math:`i`, :math:`j`, :math:`k` and :math:`l` due to the
390   force on the virtual site are computed through chain rule derivatives
391   of the construction expression. This is exact and conserves energy,
392   but it does lead to relatively lengthy expressions that we do not
393   include here (over 200 floating-point operations). The interested
394   reader can look at the source code in ``vsite.c``. Fortunately, this
395   vsite type is normally only used for chiral centers such as
396   :math:`C_{\alpha}` atoms in proteins.
397
398   The new 4fdn construct is identified with a ‘type’ value of 2 in the
399   topology. The earlier 4fd type is still supported internally (‘type’
400   value 1), but it should not be used for new simulations. All current
401   |Gromacs| tools will automatically generate type 4fdn instead.
402
403-  A linear combination of :math:`N` atoms with relative
404   weights :math:`a_i`. The weight for atom :math:`i` is:
405
406   .. math:: w_i = a_i \left(\sum_{j=1}^N a_j \right)^{-1}
407             :label: eqnvsiterelweight
408
409-   There are three options for setting the weights:
410
411   -  center of geometry: equal weights
412
413   -  center of mass: :math:`a_i` is the mass of atom :math:`i`; when in
414      free-energy simulations the mass of the atom is changed, only the
415      mass of the A-state is used for the weight
416
417   -  center of weights: :math:`a_i` is defined by the user
418
419