2 % This file is part of the GROMACS molecular simulation package.
4 % Copyright (c) 2013, 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.
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.
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.
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.
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.
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.
35 \chapter{Interaction function and force fields\index{force field}}
37 To accommodate the potential functions used
38 in some popular force fields (see \ref{sec:ff}), {\gromacs} offers a choice of functions,
39 both for non-bonded interaction and for dihedral interactions. They
40 are described in the appropriate subsections.
42 The potential functions can be subdivided into three parts
44 \item {\em Non-bonded}: Lennard-Jones or Buckingham, and Coulomb or
45 modified Coulomb. The non-bonded interactions are computed on the
46 basis of a neighbor list (a list of non-bonded atoms within a certain
47 radius), in which exclusions are already removed.
48 \item {\em Bonded}: covalent bond-stretching, angle-bending,
49 improper dihedrals, and proper dihedrals. These are computed on the
51 \item {\em Restraints}: position restraints, angle restraints,
52 distance restraints, orientation restraints and dihedral restraints, all
56 \section{Non-bonded interactions}
57 Non-bonded interactions in {\gromacs} are pair-additive and centro-symmetric:
59 V(\ve{r}_1,\ldots \ve{r}_N) = \sum_{i<j}V_{ij}(\rvij);
62 \ve{F}_i = -\sum_j \frac{dV_{ij}(r_{ij})}{dr_{ij}} \frac{\rvij}{r_{ij}} = -\ve{F}_j
64 The non-bonded interactions contain a \normindex{repulsion} term,
65 a \normindex{dispersion}
66 term, and a Coulomb term. The repulsion and dispersion term are
67 combined in either the Lennard-Jones (or 6-12 interaction), or the
68 Buckingham (or exp-6 potential). In addition, (partially) charged atoms
69 act through the Coulomb term.
71 \subsection{The Lennard-Jones interaction}
73 The \normindex{Lennard-Jones} potential $V_{LJ}$ between two atoms equals:
75 V_{LJ}(\rij) = \frac{C_{ij}^{(12)}}{\rij^{12}} -
76 \frac{C_{ij}^{(6)}}{\rij^6}
79 The parameters $C^{(12)}_{ij}$ and $C^{(6)}_{ij}$ depend on pairs of
80 {\em atom types}; consequently they are taken from a matrix of
81 LJ-parameters. In the Verlet cut-off scheme, the potential is shifted
82 by a constant such that it is zero at the cut-off distance.
85 \centerline{\includegraphics[width=8cm]{plots/f-lj}}
86 \caption {The Lennard-Jones interaction.}
90 The force derived from this potential is:
92 \ve{F}_i(\rvij) = \left( 12~\frac{C_{ij}^{(12)}}{\rij^{13}} -
93 6~\frac{C_{ij}^{(6)}}{\rij^7} \right) \rnorm
96 The LJ potential may also be written in the following form:
98 V_{LJ}(\rvij) = 4\epsilon_{ij}\left(\left(\frac{\sigma_{ij}} {\rij}\right)^{12}
99 - \left(\frac{\sigma_{ij}}{\rij}\right)^{6} \right)
103 In constructing the parameter matrix for the non-bonded LJ-parameters,
104 two types of \normindex{combination rule}s can be used within {\gromacs},
105 only geometric averages (type 1 in the input section of the force field file):
108 C_{ij}^{(6)} &=& \left({C_{ii}^{(6)} \, C_{jj}^{(6)}}\right)^{1/2} \\
109 C_{ij}^{(12)} &=& \left({C_{ii}^{(12)} \, C_{jj}^{(12)}}\right)^{1/2}
113 or, alternatively the Lorentz-Berthelot rules can be used. An arithmetic average is used to calculate $\sigma_{ij}$, while a geometric average is used to calculate $\epsilon_{ij}$ (type 2):
116 \sigma_{ij} &=& \frac{1}{ 2}(\sigma_{ii} + \sigma_{jj}) \\
117 \epsilon_{ij} &=& \left({\epsilon_{ii} \, \epsilon_{jj}}\right)^{1/2}
120 finally an geometric average for both parameters can be used (type 3):
123 \sigma_{ij} &=& \left({\sigma_{ii} \, \sigma_{jj}}\right)^{1/2} \\
124 \epsilon_{ij} &=& \left({\epsilon_{ii} \, \epsilon_{jj}}\right)^{1/2}
127 This last rule is used by the OPLS force field.
130 %\ifthenelse{\equal{\gmxlite}{1}}{}{
131 \subsection{\normindex{Buckingham potential}}
133 potential has a more flexible and realistic repulsion term
134 than the Lennard-Jones interaction, but is also more expensive to
135 compute. The potential form is:
137 V_{bh}(\rij) = A_{ij} \exp(-B_{ij} \rij) -
138 \frac{C_{ij}}{\rij^6}
141 \centerline{\includegraphics[width=8cm]{plots/f-bham}}
142 \caption {The Buckingham interaction.}
146 See also \figref{bham}. The force derived from this is:
148 \ve{F}_i(\rij) = \left[ A_{ij}B_{ij}\exp(-B_{ij} \rij) -
149 6\frac{C_{ij}}{\rij^7} \right] \rnorm
152 %} % Brace matches ifthenelse test for gmxlite
154 \subsection{Coulomb interaction}
156 \newcommand{\epsr}{\varepsilon_r}
157 \newcommand{\epsrf}{\varepsilon_{rf}}
158 The \normindex{Coulomb} interaction between two charge particles is given by:
160 V_c(\rij) = f \frac{q_i q_j}{\epsr \rij}
163 See also \figref{coul}, where $f = \frac{1}{4\pi \varepsilon_0} =
164 138.935\,485$ (see \chref{defunits})
167 \centerline{\includegraphics[width=8cm]{plots/vcrf}}
168 \caption[The Coulomb interaction with and without reaction field.]{The
169 Coulomb interaction (for particles with equal signed charge) with and
170 without reaction field. In the latter case $\epsr$ was 1, $\epsrf$ was 78,
171 and $r_c$ was 0.9 nm.
172 The dot-dashed line is the same as the dashed line, except for a constant.}
176 The force derived from this potential is:
178 \ve{F}_i(\rvij) = f \frac{q_i q_j}{\epsr\rij^2}\rnorm
181 In {\gromacs} the relative \swapindex{dielectric}{constant}
183 may be set in the in the input for {\tt grompp}.
185 %\ifthenelse{\equal{\gmxlite}{1}}{}{
186 \subsection{Coulomb interaction with \normindex{reaction field}}
188 The Coulomb interaction can be modified for homogeneous systems by
189 assuming a constant dielectric environment beyond the cut-off $r_c$
190 with a dielectric constant of {$\epsrf$}. The interaction then reads:
193 f \frac{q_i q_j}{\epsr\rij}\left[1+\frac{\epsrf-\epsr}{2\epsrf+\epsr}
194 \,\frac{\rij^3}{r_c^3}\right]
195 - f\frac{q_i q_j}{\epsr r_c}\,\frac{3\epsrf}{2\epsrf+\epsr}
198 in which the constant expression on the right makes the potential
199 zero at the cut-off $r_c$. For charged cut-off spheres this corresponds
200 to neutralization with a homogeneous background charge.
201 We can rewrite \eqnref{vcrf} for simplicity as
203 V_{crf} ~=~ f \frac{q_i q_j}{\epsr}\left[\frac{1}{\rij} + k_{rf}~ \rij^2 -c_{rf}\right]
207 k_{rf} &=& \frac{1}{r_c^3}\,\frac{\epsrf-\epsr}{(2\epsrf+\epsr)} \label{eqn:krf}\\
208 c_{rf} &=& \frac{1}{r_c}+k_{rf}\,r_c^2 ~=~ \frac{1}{r_c}\,\frac{3\epsrf}{(2\epsrf+\epsr)}
211 For large $\epsrf$ the $k_{rf}$ goes to $r_c^{-3}/2$,
212 while for $\epsrf$ = $\epsr$ the correction vanishes.
214 the modified interaction is plotted, and it is clear that the derivative
215 with respect to $\rij$ (= -force) goes to zero at the cut-off distance.
216 The force derived from this potential reads:
218 \ve{F}_i(\rvij) = f \frac{q_i q_j}{\epsr}\left[\frac{1}{\rij^2} - 2 k_{rf}\rij\right]\rnorm
221 The reaction-field correction should also be applied to all excluded
222 atoms pairs, including self pairs, in which case the normal Coulomb
223 term in \eqnsref{vcrf}{fcrf} is absent.
225 Tironi {\etal} have introduced a generalized reaction field in which
226 the dielectric continuum beyond the cut-off $r_c$ also has an ionic strength
227 $I$~\cite{Tironi95}. In this case we can rewrite the constants $k_{rf}$ and
228 $c_{rf}$ using the inverse Debye screening length $\kappa$:
231 \frac{2 I \,F^2}{\varepsilon_0 \epsrf RT}
232 = \frac{F^2}{\varepsilon_0 \epsrf RT}\sum_{i=1}^{K} c_i z_i^2 \\
233 k_{rf} &=& \frac{1}{r_c^3}\,
234 \frac{(\epsrf-\epsr)(1 + \kappa r_c) + \half\epsrf(\kappa r_c)^2}
235 {(2\epsrf + \epsr)(1 + \kappa r_c) + \epsrf(\kappa r_c)^2}
237 c_{rf} &=& \frac{1}{r_c}\,
238 \frac{3\epsrf(1 + \kappa r_c + \half(\kappa r_c)^2)}
239 {(2\epsrf+\epsr)(1 + \kappa r_c) + \epsrf(\kappa r_c)^2}
242 where $F$ is Faraday's constant, $R$ is the ideal gas constant, $T$
243 the absolute temperature, $c_i$ the molar concentration for species
244 $i$ and $z_i$ the charge number of species $i$ where we have $K$
245 different species. In the limit of zero ionic strength ($\kappa=0$)
246 \eqnsref{kgrf}{cgrf} reduce to the simple forms of \eqnsref{krf}{crf}
249 \subsection{Modified non-bonded interactions}
250 \label{sec:mod_nb_int}
251 In {\gromacs}, the non-bonded potentials can be
252 modified by a shift function. The purpose of this is to replace the
253 truncated forces by forces that are continuous and have continuous
254 derivatives at the \normindex{cut-off} radius. With such forces the
255 timestep integration produces much smaller errors and there are no
256 such complications as creating charges from dipoles by the truncation
257 procedure. In fact, by using shifted forces there is no need for
258 charge groups in the construction of neighbor lists. However, the
259 shift function produces a considerable modification of the Coulomb
260 potential. Unless the ``missing'' long-range potential is properly
261 calculated and added (through the use of PPPM, Ewald, or PME), the
262 effect of such modifications must be carefully evaluated. The
263 modification of the Lennard-Jones dispersion and repulsion is only
264 minor, but it does remove the noise caused by cut-off effects.
266 There is {\em no} fundamental difference between a switch function
267 (which multiplies the potential with a function) and a shift function
268 (which adds a function to the force or potential)~\cite{Spoel2006a}. The switch
269 function is a special case of the shift function, which we apply to
270 the {\em force function} $F(r)$, related to the electrostatic or
271 van der Waals force acting on particle $i$ by particle $j$ as:
273 \ve{F}_i = c F(r_{ij}) \frac{\rvij}{r_{ij}}
275 For pure Coulomb or Lennard-Jones interactions
276 $F(r)=F_\alpha(r)=r^{-(\alpha+1)}$.
277 The shifted force $F_s(r)$ can generally be written as:
281 F_s(r)~=&~F_\alpha(r) & r < r_1 \\
283 F_s(r)~=&~F_\alpha(r)+S(r) & r_1 \le r < r_c \\
284 F_s(r)~=&~0 & r_c \le r
287 When $r_1=0$ this is a traditional shift function, otherwise it acts as a
288 switch function. The corresponding shifted coulomb potential then reads:
290 V_s(r_{ij}) = f \Phi_s (r_{ij}) q_i q_j
292 where $\Phi(r)$ is the potential function
294 \Phi_s(r) = \int^{\infty}_r~F_s(x)\, dx
297 The {\gromacs} shift function should be smooth at the boundaries, therefore
298 the following boundary conditions are imposed on the shift function:
303 S(r_c) &=&-F_\alpha(r_c) \\
304 S'(r_c) &=&-F_\alpha'(r_c)
307 A 3$^{rd}$ degree polynomial of the form
309 S(r) = A(r-r_1)^2 + B(r-r_1)^3
311 fulfills these requirements. The constants A and B are given by the
312 boundary condition at $r_c$:
316 A &~=~& -\displaystyle
317 \frac{(\alpha+4)r_c~-~(\alpha+1)r_1} {r_c^{\alpha+2}~(r_c-r_1)^2} \\
318 B &~=~& \displaystyle
319 \frac{(\alpha+3)r_c~-~(\alpha+1)r_1}{r_c^{\alpha+2}~(r_c-r_1)^3}
322 Thus the total force function is:
324 F_s(r) = \frac{\alpha}{r^{\alpha+1}} + A(r-r_1)^2 + B(r-r_1)^3
326 and the potential function reads:
328 \Phi(r) = \frac{1}{r^\alpha} - \frac{A}{3} (r-r_1)^3 - \frac{B}{4} (r-r_1)^4 - C
332 C = \frac{1}{r_c^\alpha} - \frac{A}{3} (r_c-r_1)^3 - \frac{B}{4} (r_c-r_1)^4
335 When $r_1$ = 0, the modified Coulomb force function is
337 F_s(r) = \frac{1}{r^2} - \frac{5 r^2}{r_c^4} + \frac{4 r^3}{r_c^5}
339 which is identical to the {\em \index{parabolic force}}
340 function recommended to be used as a short-range function in
341 conjunction with a \swapindex{Poisson}{solver}
342 for the long-range part~\cite{Berendsen93a}.
343 The modified Coulomb potential function is:
345 \Phi(r) = \frac{1}{r} - \frac{5}{3r_c} + \frac{5r^3}{3r_c^4} - \frac{r^4}{r_c^5}
347 See also \figref{shift}.
350 \centerline{\includegraphics[angle=270,width=10cm]{plots/shiftf}}
351 \caption[The Coulomb Force, Shifted Force and Shift Function
352 $S(r)$,.]{The Coulomb Force, Shifted Force and Shift Function $S(r)$,
353 using r$_1$ = 2 and r$_c$ = 4.}
357 \subsection{Modified short-range interactions with Ewald summation}
358 When Ewald summation\index{Ewald sum} or \seeindex{particle-mesh
359 Ewald}{PME}\index{Ewald, particle-mesh} is used to calculate the
360 long-range interactions, the
361 short-range Coulomb potential must also be modified, similar to the
362 switch function above. In this case the short range potential is given
365 V(r) = f \frac{\mbox{erfc}(\beta r_{ij})}{r_{ij}} q_i q_j,
367 where $\beta$ is a parameter that determines the relative weight
368 between the direct space sum and the reciprocal space sum and erfc$(x)$ is
369 the complementary error function. For further
370 details on long-range electrostatics, see \secref{lr_elstat}.
371 %} % Brace matches ifthenelse test for gmxlite
374 \section{Bonded interactions}
375 Bonded interactions are based on a fixed list of atoms. They are not
376 exclusively pair interactions, but include 3- and 4-body interactions
377 as well. There are {\em bond stretching} (2-body), {\em bond angle}
378 (3-body), and {\em dihedral angle} (4-body) interactions. A special
379 type of dihedral interaction (called {\em improper dihedral}) is used
380 to force atoms to remain in a plane or to prevent transition to a
381 configuration of opposite chirality (a mirror image).
383 \subsection{Bond stretching}
385 \subsubsection{Harmonic potential}
386 \label{subsec:harmonicbond}
387 The \swapindex{bond}{stretching} between two covalently bonded atoms
388 $i$ and $j$ is represented by a harmonic potential:
391 \centerline{\raisebox{4cm}{\includegraphics[angle=270,width=5cm]{plots/bstretch}}\includegraphics[width=7cm]{plots/f-bond}}
392 \caption[Bond stretching.]{Principle of bond stretching (left), and the bond
393 stretching potential (right).}
394 \label{fig:bstretch1}
398 V_b~(\rij) = \half k^b_{ij}(\rij-b_{ij})^2
400 See also \figref{bstretch1}, with the force given by:
402 \ve{F}_i(\rvij) = k^b_{ij}(\rij-b_{ij}) \rnorm
405 %\ifthenelse{\equal{\gmxlite}{1}}{}{
406 \subsubsection{Fourth power potential}
407 \label{subsec:G96bond}
408 In the \gromosv{96} force field~\cite{gromos96}, the covalent bond potential
409 is, for reasons of computational efficiency, written as:
411 V_b~(\rij) = \frac{1}{4}k^b_{ij}\left(\rij^2-b_{ij}^2\right)^2
413 The corresponding force is:
415 \ve{F}_i(\rvij) = k^b_{ij}(\rij^2-b_{ij}^2)~\rvij
417 The force constants for this form of the potential are related to the usual
418 harmonic force constant $k^{b,harm}$ (\secref{bondpot}) as
420 2 k^b b_{ij}^2 = k^{b,harm}
422 The force constants are mostly derived from the harmonic ones used in
423 \gromosv{87}~\cite{biomos}. Although this form is computationally more
425 (because no square root has to be evaluated), it is conceptually more
426 complex. One particular disadvantage is that since the form is not harmonic,
427 the average energy of a single bond is not equal to $\half kT$ as it is for
428 the normal harmonic potential.
430 \subsection{\normindex{Morse potential} bond stretching}
431 \label{subsec:Morsebond}
432 %\author{F.P.X. Everdij}
434 For some systems that require an anharmonic bond stretching potential,
435 the Morse potential~\cite{Morse29}
436 between two atoms {\it i} and {\it j} is available
437 in {\gromacs}. This potential differs from the harmonic potential in
438 that it has an asymmetric potential well and a zero force at infinite
439 distance. The functional form is:
441 \displaystyle V_{morse} (r_{ij}) = D_{ij} [1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))]^2,
443 See also \figref{morse}, and the corresponding force is:
446 \displaystyle {\bf F}_{morse} ({\bf r}_{ij})&=&2 D_{ij} \beta_{ij} r_{ij} \exp(-\beta_{ij}(r_{ij}-b_{ij})) * \\
447 \displaystyle \: & \: &[1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))] \frac{\displaystyle {\bf r}_{ij}}{\displaystyle r_{ij}},
450 where \( \displaystyle D_{ij} \) is the depth of the well in kJ/mol,
451 \( \displaystyle \beta_{ij} \) defines the steepness of the well (in
452 nm\(^{-1} \)), and \( \displaystyle b_{ij} \) is the equilibrium
453 distance in nm. The steepness parameter \( \displaystyle \beta_{ij}
454 \) can be expressed in terms of the reduced mass of the atoms {\it i}
455 and {\it j}, the fundamental vibration frequency \( \displaystyle
456 \omega_{ij} \) and the well depth \( \displaystyle D_{ij} \):
458 \displaystyle \beta_{ij}= \omega_{ij} \sqrt{\frac{\mu_{ij}}{2 D_{ij}}}
460 and because \( \displaystyle \omega = \sqrt{k/\mu} \), one can rewrite \( \displaystyle \beta_{ij} \) in terms of the harmonic force constant \( \displaystyle k_{ij} \):
462 \displaystyle \beta_{ij}= \sqrt{\frac{k_{ij}}{2 D_{ij}}}
465 For small deviations \( \displaystyle (r_{ij}-b_{ij}) \), one can
466 approximate the \( \displaystyle \exp \)-term to first-order using a
469 \displaystyle \exp(-x) \approx 1-x
472 and substituting \eqnref{betaij} and \eqnref{expminx} in the functional form:
475 \displaystyle V_{morse} (r_{ij})&=&D_{ij} [1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))]^2\\
476 \displaystyle \:&=&D_{ij} [1 - (1 -\sqrt{\frac{k_{ij}}{2 D_{ij}}}(r_{ij}-b_{ij}))]^2\\
477 \displaystyle \:&=&\frac{1}{2} k_{ij} (r_{ij}-b_{ij}))^2
480 we recover the harmonic bond stretching potential.
483 \centerline{\includegraphics[width=7cm]{plots/f-morse}}
484 \caption{The Morse potential well, with bond length 0.15 nm.}
488 \subsection{Cubic bond stretching potential}
489 \label{subsec:cubicbond}
490 Another anharmonic bond stretching potential that is slightly simpler
491 than the Morse potential adds a cubic term in the distance to the
492 simple harmonic form:
494 V_b~(\rij) = k^b_{ij}(\rij-b_{ij})^2 + k^b_{ij}k^{cub}_{ij}(\rij-b_{ij})^3
496 A flexible \normindex{water} model (based on
497 the SPC water model~\cite{Berendsen81}) including
498 a cubic bond stretching potential for the O-H bond
499 was developed by Ferguson~\cite{Ferguson95}. This model was found
500 to yield a reasonable infrared spectrum. The Ferguson water model is
501 available in the {\gromacs} library ({\tt flexwat-ferguson.itp}).
502 It should be noted that the potential is asymmetric: overstretching leads to
503 infinitely low energies. The \swapindex{integration}{timestep} is therefore
506 The force corresponding to this potential is:
508 \ve{F}_i(\rvij) = 2k^b_{ij}(\rij-b_{ij})~\rnorm + 3k^b_{ij}k^{cub}_{ij}(\rij-b_{ij})^2~\rnorm
511 \subsection{FENE bond stretching potential\index{FENE potential}}
512 \label{subsec:FENEbond}
513 In coarse-grained polymer simulations the beads are often connected
514 by a FENE (finitely extensible nonlinear elastic) potential~\cite{Warner72}:
516 V_{\mbox{\small FENE}}(\rij) =
517 -\half k^b_{ij} b^2_{ij} \log\left(1 - \frac{\rij^2}{b^2_{ij}}\right)
519 The potential looks complicated, but the expression for the force is simpler:
521 F_{\mbox{\small FENE}}(\rvij) =
522 -k^b_{ij} \left(1 - \frac{\rij^2}{b^2_{ij}}\right)^{-1} \rvij
524 At short distances the potential asymptotically goes to a harmonic
525 potential with force constant $k^b$, while it diverges at distance $b$.
526 %} % Brace matches ifthenelse test for gmxlite
528 \subsection{Harmonic angle potential}
529 \label{subsec:harmonicangle}
530 \newcommand{\tijk}{\theta_{ijk}}
531 The bond-\swapindex{angle}{vibration} between a triplet of atoms $i$ - $j$ - $k$
532 is also represented by a harmonic potential on the angle $\tijk$
535 \centerline{\raisebox{4cm}{\includegraphics[angle=270,width=5cm]{plots/angle}}\includegraphics[width=7cm]{plots/f-angle}}
536 \caption[Angle vibration.]{Principle of angle vibration (left) and the
537 bond angle potential (right).}
542 V_a(\tijk) = \half k^{\theta}_{ijk}(\tijk-\tijk^0)^2
544 As the bond-angle vibration is represented by a harmonic potential, the
545 form is the same as the bond stretching (\figref{bstretch1}).
547 The force equations are given by the chain rule:
550 \Fvi ~=~ -\displaystyle\frac{d V_a(\tijk)}{d \rvi} \\
551 \Fvk ~=~ -\displaystyle\frac{d V_a(\tijk)}{d \rvk} \\
554 ~ \mbox{ ~ where ~ } ~
555 \tijk = \arccos \frac{(\rvij \cdot \ve{r}_{kj})}{r_{ij}r_{kj}}
557 The numbering $i,j,k$ is in sequence of covalently bonded atoms. Atom
558 $j$ is in the middle; atoms $i$ and $k$ are at the ends (see \figref{angle}).
559 {\bf Note} that in the input in topology files, angles are given in degrees and
560 force constants in kJ/mol/rad$^2$.
562 %\ifthenelse{\equal{\gmxlite}{1}}{}{
563 \subsection{Cosine based angle potential}
564 \label{subsec:G96angle}
565 In the \gromosv{96} force field a simplified function is used to represent angle
568 V_a(\tijk) = \half k^{\theta}_{ijk}\left(\cos(\tijk) - \cos(\tijk^0)\right)^2
572 \cos(\tijk) = \frac{\rvij\cdot\ve{r}_{kj}}{\rij r_{kj}}
574 The corresponding force can be derived by partial differentiation with respect
575 to the atomic positions. The force constants in this function are related
576 to the force constants in the harmonic form $k^{\theta,harm}$
577 (\ssecref{harmonicangle}) by:
579 k^{\theta} \sin^2(\tijk^0) = k^{\theta,harm}
581 In the \gromosv{96} manual there is a much more complicated conversion formula
582 which is temperature dependent. The formulas are equivalent at 0 K
583 and the differences at 300 K are on the order of 0.1 to 0.2\%.
584 {\bf Note} that in the input in topology files, angles are given in degrees and
585 force constants in kJ/mol.
587 \subsection{Urey-Bradley potential}
588 \label{subsec:Urey-Bradley}
589 The \swapindex{Urey-Bradley bond-angle}{vibration} between a triplet
590 of atoms $i$ - $j$ - $k$ is represented by a harmonic potential on the
591 angle $\tijk$ and a harmonic correction term on the distance between
592 the atoms $i$ and $k$. Although this can be easily written as a simple
593 sum of two terms, it is convenient to have it as a single entry in the
594 topology file and in the output as a separate energy term. It is used mainly
595 in the CHARMm force field~\cite{BBrooks83}. The energy is given by:
598 V_a(\tijk) = \half k^{\theta}_{ijk}(\tijk-\tijk^0)^2 + \half k^{UB}_{ijk}(r_{ik}-r_{ik}^0)^2
601 The force equations can be deduced from sections~\ssecref{harmonicbond}
602 and~\ssecref{harmonicangle}.
604 \subsection{Bond-Bond cross term}
605 \label{subsec:bondbondcross}
606 The bond-bond cross term for three particles $i, j, k$ forming bonds
607 $i-j$ and $k-j$ is given by~\cite{Lawrence2003b}:
609 V_{rr'} ~=~ k_{rr'} \left(\left|\ve{r}_{i}-\ve{r}_j\right|-r_{1e}\right) \left(\left|\ve{r}_{k}-\ve{r}_j\right|-r_{2e}\right)
612 where $k_{rr'}$ is the force constant, and $r_{1e}$ and $r_{2e}$ are the
613 equilibrium bond lengths of the $i-j$ and $k-j$ bonds respectively. The force
614 associated with this potential on particle $i$ is:
616 \ve{F}_{i} = -k_{rr'}\left(\left|\ve{r}_{k}-\ve{r}_j\right|-r_{2e}\right)\frac{\ve{r}_i-\ve{r}_j}{\left|\ve{r}_{i}-\ve{r}_j\right|}
618 The force on atom $k$ can be obtained by swapping $i$ and $k$ in the above
619 equation. Finally, the force on atom $j$ follows from the fact that the sum
620 of internal forces should be zero: $\ve{F}_j = -\ve{F}_i-\ve{F}_k$.
622 \subsection{Bond-Angle cross term}
623 \label{subsec:bondanglecross}
624 The bond-angle cross term for three particles $i, j, k$ forming bonds
625 $i-j$ and $k-j$ is given by~\cite{Lawrence2003b}:
627 V_{r\theta} ~=~ k_{r\theta} \left(\left|\ve{r}_{i}-\ve{r}_k\right|-r_{3e} \right) \left(\left|\ve{r}_{i}-\ve{r}_j\right|-r_{1e} + \left|\ve{r}_{k}-\ve{r}_j\right|-r_{2e}\right)
629 where $k_{r\theta}$ is the force constant, $r_{3e}$ is the $i-k$ distance,
630 and the other constants are the same as in Equation~\ref{crossbb}. The force
631 associated with the potential on atom $i$ is:
633 \ve{F}_{i} ~=~ -k_{r\theta}\left[\left(\left|\ve{r}_{i}-\ve{r}_{k}\right|-r_{3e}\right)\frac{\ve{r}_i-\ve{r}_j}{\left|\ve{r}_{i}-\ve{r}_j\right|} \\
634 + \left(\left|\ve{r}_{i}-\ve{r}_j\right|-r_{1e} + \left|\ve{r}_{k}-\ve{r}_j\right|-r_{2e}\right)\frac{\ve{r}_i-\ve{r}_k}{\left|\ve{r}_{i}-\ve{r}_k\right|}\right]
637 \subsection{Quartic angle potential}
638 \label{subsec:quarticangle}
639 For special purposes there is an angle potential
640 that uses a fourth order polynomial:
642 V_q(\tijk) ~=~ \sum_{n=0}^5 C_n (\tijk-\tijk^0)^n
644 %} % Brace matches ifthenelse test for gmxlite
646 %% new commands %%%%%%%%%%%%%%%%%%%%%%
647 \newcommand{\rvkj}{{\bf r}_{kj}}
648 \newcommand{\rkj}{r_{kj}}
649 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
651 \subsection{Improper dihedrals\swapindexquiet{improper}{dihedral}}
653 Improper dihedrals are meant to keep \swapindex{planar}{group}s ({\eg}
654 aromatic rings) planar, or to prevent molecules from flipping over to their
655 \normindex{mirror image}s, see \figref{imp}.
658 \centerline{\includegraphics[angle=270,width=4cm]{plots/ring-imp}\hspace{1cm}
659 \includegraphics[angle=270,width=3cm]{plots/subst-im}\hspace{1cm}\includegraphics[angle=270,width=3cm]{plots/tetra-im}}
660 \caption[Improper dihedral angles.]{Principle of improper
661 dihedral angles. Out of plane bending for rings (left), substituents
662 of rings (middle), out of tetrahedral (right). The improper dihedral
663 angle $\xi$ is defined as the angle between planes (i,j,k) and (j,k,l)
668 \subsubsection{Improper dihedrals: harmonic type}
669 \label{subsec:harmonicimproperdihedral}
670 The simplest improper dihedral potential is a harmonic potential; it is plotted in
673 V_{id}(\xi_{ijkl}) = \half k_{\xi}(\xi_{ijkl}-\xi_0)^2
675 Since the potential is harmonic it is discontinuous,
676 but since the discontinuity is chosen at 180$^\circ$ distance from $\xi_0$
677 this will never cause problems.
678 {\bf Note} that in the input in topology files, angles are given in degrees and
679 force constants in kJ/mol/rad$^2$.
682 \centerline{\includegraphics[width=8cm]{plots/f-imps}}
683 \caption{Improper dihedral potential.}
687 \subsubsection{Improper dihedrals: periodic type}
688 \label{subsec:periodicimproperdihedral}
689 This potential is identical to the periodic proper dihedral (see below).
690 There is a separate dihedral type for this (type 4) only to be able
691 to distinguish improper from proper dihedrals in the parameter section
694 \subsection{Proper dihedrals\swapindexquiet{proper}{dihedral}}
695 For the normal \normindex{dihedral} interaction there is a choice of
696 either the {\gromos} periodic function or a function based on
697 expansion in powers of $\cos \phi$ (the so-called Ryckaert-Bellemans
698 potential). This choice has consequences for the inclusion of special
699 interactions between the first and the fourth atom of the dihedral
700 quadruple. With the periodic {\gromos} potential a special 1-4
701 LJ-interaction must be included; with the Ryckaert-Bellemans potential
702 {\em for alkanes} the \normindex{1-4 interaction}s must be excluded
703 from the non-bonded list. {\bf Note:} Ryckaert-Bellemans potentials
704 are also used in {\eg} the OPLS force field in combination with 1-4
705 interactions. You should therefore not modify topologies generated by
706 {\tt \normindex{pdb2gmx}} in this case.
708 \subsubsection{Proper dihedrals: periodic type}
709 \label{subsec:properdihedral}
710 Proper dihedral angles are defined according to the IUPAC/IUB
711 convention, where $\phi$ is the angle between the $ijk$ and the $jkl$
712 planes, with {\bf zero} corresponding to the {\em cis} configuration
713 ($i$ and $l$ on the same side). There are two dihedral function types
714 in {\gromacs} topology files. There is the standard type 1 which behaves
715 like any other bonded interactions. For certain force fields, type 9
716 is useful. Type 9 allows multiple potential functions to be applied
717 automatically to a single dihedral in the {\tt [ dihedral ]} section
718 when multiple parameters are defined for the same atomtypes
719 in the {\tt [ dihedraltypes ]} section.
722 \centerline{\raisebox{4.5cm}{\includegraphics[angle=270,width=5cm]{plots/dih}}\includegraphics[width=7cm]{plots/f-dih}}
723 \caption[Proper dihedral angle.]{Principle of proper dihedral angle
724 (left, in {\em trans} form) and the dihedral angle potential (right).}
728 V_d(\phi_{ijkl}) = k_{\phi}(1 + \cos(n \phi - \phi_s))
731 %\ifthenelse{\equal{\gmxlite}{1}}{}{
732 \subsubsection{Proper dihedrals: Ryckaert-Bellemans function}
733 \label{subsec:RBdihedral}
734 For alkanes, the following proper dihedral potential is often used
735 (see \figref{rbdih}):
737 V_{rb}(\phi_{ijkl}) = \sum_{n=0}^5 C_n( \cos(\psi ))^n,
739 where $\psi = \phi - 180^\circ$. \\
740 {\bf Note:} A conversion from one convention to another can be achieved by
741 multiplying every coefficient \( \displaystyle C_n \)
742 by \( \displaystyle (-1)^n \).
744 An example of constants for $C$ is given in \tabref{crb}.
748 \begin{tabular}{|lr|lr|lr|}
750 $C_0$ & 9.28 & $C_2$ & -13.12 & $C_4$ & 26.24 \\
751 $C_1$ & 12.16 & $C_3$ & -3.06 & $C_5$ & -31.5 \\
755 \caption{Constants for Ryckaert-Bellemans potential (kJ mol$^{-1}$).}
760 \centerline{\includegraphics[width=8cm]{plots/f-rbs}}
761 \caption{Ryckaert-Bellemans dihedral potential.}
765 ({\bf Note:} The use of this potential implies exclusion of LJ interactions
766 between the first and the last atom of the dihedral, and $\psi$ is defined
767 according to the ``polymer convention'' ($\psi_{trans}=0$).)
769 The RB dihedral function can also be used to include Fourier dihedrals
772 V_{rb} (\phi_{ijkl}) ~=~ \frac{1}{2} \left[F_1(1+\cos(\phi)) + F_2(
773 1-\cos(2\phi)) + F_3(1+\cos(3\phi)) + F_4(1-\cos(4\phi))\right]
775 Because of the equalities \( \cos(2\phi) = 2\cos^2(\phi) - 1 \),
776 \( \cos(3\phi) = 4\cos^3(\phi) - 3\cos(\phi) \) and
777 \( \cos(4\phi) = 8\cos^4(\phi) - 8\cos^2(\phi) + 1 \)
778 one can translate the OPLS parameters to
779 Ryckaert-Bellemans parameters as follows:
783 \displaystyle C_0&=&F_2 + \frac{1}{2} (F_1 + F_3)\\
784 \displaystyle C_1&=&\frac{1}{2} (- F_1 + 3 \, F_3)\\
785 \displaystyle C_2&=& -F_2 + 4 \, F_4\\
786 \displaystyle C_3&=&-2 \, F_3\\
787 \displaystyle C_4&=&-4 \, F_4\\
788 \displaystyle C_5&=&0
791 with OPLS parameters in protein convention and RB parameters in
792 polymer convention (this yields a minus sign for the odd powers of
794 \noindent{\bf Note:} Mind the conversion from {\bf kcal mol$^{-1}$} for
795 literature OPLS and RB parameters to {\bf kJ mol$^{-1}$} in {\gromacs}.\\
796 %} % Brace matches ifthenelse test for gmxlite
798 \subsubsection{Proper dihedrals: Fourier function}
799 \label{subsec:Fourierdihedral}
800 The OPLS potential function is given as the first three
801 or four~\cite{Jorgensen2005a} cosine terms of a Fourier series.
802 In {\gromacs} the four term function is implemented:
804 V_{F} (\phi_{ijkl}) ~=~ \frac{1}{2} \left[C_1(1+\cos(\phi)) + C_2(
805 1-\cos(2\phi)) + C_3(1+\cos(3\phi)) + C_4(1+\cos(4\phi))\right],
807 %\ifthenelse{\equal{\gmxlite}{1}}{}{
808 Internally, {\gromacs}
809 uses the Ryckaert-Bellemans code
810 to compute Fourier dihedrals (see above), because this is more efficient.\\
811 \noindent{\bf Note:} Mind the conversion from {\emph kcal mol$^{-1}$} for
812 literature OPLS parameters to {\bf kJ mol$^{-1}$} in {\gromacs}.\\
814 %\ifthenelse{\equal{\gmxlite}{1}}{}{
815 \subsection{Tabulated bonded interaction functions\index{tabulated bonded interaction function}}
816 \label{subsec:tabulatedinteraction}
817 For full flexibility, any functional shape can be used for
818 bonds, angles and dihedrals through user-supplied tabulated functions.
819 The functional shapes are:
821 V_b(r_{ij}) &=& k \, f^b_n(r_{ij}) \\
822 V_a(\tijk) &=& k \, f^a_n(\tijk) \\
823 V_d(\phi_{ijkl}) &=& k \, f^d_n(\phi_{ijkl})
825 where $k$ is a force constant in units of energy
826 and $f$ is a cubic spline function; for details see \ssecref{cubicspline}.
827 For each interaction, the force constant $k$ and the table number $n$
828 are specified in the topology.
829 There are two different types of bonds, one that generates exclusions (type 8)
830 and one that does not (type 9).
831 For details see \tabref{topfile2}.
832 The table files are supplied to the {\tt mdrun} program.
833 After the table file name an underscore, the letter ``b'' for bonds,
834 ``a'' for angles or ``d'' for dihedrals and the table number are appended.
835 For example, for a bond with $n=0$ (and using the default table file name)
836 the table is read from the file {\tt table_b0.xvg}. Multiple tables can be
837 supplied simply by using different values of $n$, and are applied to the appropriate
838 bonds, as specified in the topology (\tabref{topfile2}).
839 The format for the table files is three columns with $x$, $f(x)$, $-f'(x)$,
840 where $x$ should be uniformly-spaced. Requirements for entries in the topology
841 are given in~\tabref{topfile2}.
842 The setup of the tables is as follows:
844 $x$ is the distance in nm. For distances beyond the table length,
845 {\tt mdrun} will quit with an error message.
847 $x$ is the angle in degrees. The table should go from
848 0 up to and including 180 degrees; the derivative is taken in degrees.
850 $x$ is the dihedral angle in degrees. The table should go from
851 -180 up to and including 180 degrees;
852 the IUPAC/IUB convention is used, {\ie} zero is cis,
853 the derivative is taken in degrees.
854 %} % Brace matches ifthenelse test for gmxlite
857 Special potentials are used for imposing restraints on the motion of
858 the system, either to avoid disastrous deviations, or to include
859 knowledge from experimental data. In either case they are not really
860 part of the force field and the reliability of the parameters is not
861 important. The potential forms, as implemented in {\gromacs}, are
862 mentioned just for the sake of completeness. Restraints and constraints
863 refer to quite different algorithms in {\gromacs}.
865 \subsection{Position restraints\swapindexquiet{position}{restraint}}
866 \label{subsec:positionrestraint}
867 These are used to restrain particles to fixed reference positions
868 $\ve{R}_i$. They can be used during equilibration in order to avoid
869 drastic rearrangements of critical parts ({\eg} to restrain motion
870 in a protein that is subjected to large solvent forces when the
871 solvent is not yet equilibrated). Another application is the
872 restraining of particles in a shell around a region that is simulated
873 in detail, while the shell is only approximated because it lacks
874 proper interaction from missing particles outside the
875 shell. Restraining will then maintain the integrity of the inner
876 part. For spherical shells, it is a wise procedure to make the force
877 constant depend on the radius, increasing from zero at the inner
878 boundary to a large value at the outer boundary. This feature has
879 not, however, been implemented in {\gromacs}.
880 \newcommand{\unitv}[1]{\hat{\bf #1}}
881 \newcommand{\halfje}[1]{\frac{#1}{2}}
883 The following form is used:
885 V_{pr}(\ve{r}_i) = \halfje{1}k_{pr}|\rvi-\ve{R}_i|^2
887 The potential is plotted in \figref{positionrestraint}.
890 \centerline{\includegraphics[width=8cm]{plots/f-pr}}
891 \caption{Position restraint potential.}
892 \label{fig:positionrestraint}
895 The potential form can be rewritten without loss of generality as:
897 V_{pr}(\ve{r}_i) = \halfje{1} \left[ k_{pr}^x (x_i-X_i)^2 ~\unitv{x} + k_{pr}^y (y_i-Y_i)^2 ~\unitv{y} + k_{pr}^z (z_i-Z_i)^2 ~\unitv{z}\right]
903 F_i^x &=& -k_{pr}^x~(x_i - X_i) \\
904 F_i^y &=& -k_{pr}^y~(y_i - Y_i) \\
905 F_i^z &=& -k_{pr}^z~(z_i - Z_i)
908 Using three different force constants the position
909 restraints can be turned on or off
910 in each spatial dimension; this means that atoms can be harmonically
911 restrained to a plane or a line.
912 Position restraints are applied to a special fixed list of atoms. Such a
913 list is usually generated by the {\tt \normindex{pdb2gmx}} program.
915 \subsection{\swapindex{Flat-bottomed}{position restraint}s}
916 \label{subsec:fbpositionrestraint}
917 Flat-bottomed position restraints can be used to restrain particles to
918 part of the simulation volume. No force acts on the restrained
919 particle within the flat-bottomed region of the potential, however a
920 harmonic force acts to move the particle to the flat-bottomed region
921 if it is outside it. It is possible to apply normal and
922 flat-bottomed position restraints on the same particle (however, only
923 with the same reference position $\ve{R}_i$). The following general potential
924 is used (Figure~\ref{fig:fbposres}A):
926 V_\mathrm{fb}(\ve{r}_i) = \frac{1}{2}k_\mathrm{fb} [d_g(\ve{r}_i;\ve{R}_i) - r_\mathrm{fb}]^2\,H[d_g(\ve{r}_i;\ve{R}_i) - r_\mathrm{fb}],
928 where $\ve{R}_i$ is the reference position, $r_\mathrm{fb}$ is the distance
929 from the center with a flat potential, $k_\mathrm{fb}$ the force constant, and $H$ is the Heaviside step
930 function. The distance $d_g(\ve{r}_i;\ve{R}_i)$ from the reference
931 position depends on the geometry $g$ of the flat-bottomed potential.
934 \centerline{\includegraphics[width=10cm]{plots/fbposres}}
935 \caption{Flat-bottomed position restraint potential. (A) Not
936 inverted, (B) inverted.}
940 The following geometries for the flat-bottomed potential are supported:\newline
941 {\bfseries Sphere} ($g =1$): The particle is kept in a sphere of given
942 radius. The force acts towards the center of the sphere. The following distance calculation is used:
944 d_g(\ve{r}_i;\ve{R}_i) = |\ve{r}_i-\ve{R}_i|
946 {\bfseries Cylinder} ($g=2$): The particle is kept in a cylinder of given radius
947 parallel to the $z$-axis. The force from the flat-bottomed potential acts
948 towards the axis of the cylinder. The $z$-component of the force is zero.
950 d_g(\ve{r}_i;\ve{R}_i) = \sqrt{ (x_i-X_i)^2 + (y_i - Y_i)^2 }
952 {\bfseries Layer} ($g=3,4,5$): The particle is kept in a layer defined by the
953 thickness and the normal of the layer. The layer normal can be parallel to the $x$, $y$, or
954 $z$-axis. The force acts parallel to the layer normal.\\
956 d_g(\ve{r}_i;\ve{R}_i) = |x_i-X_i|, \;\;\;\mbox{or}\;\;\;
957 d_g(\ve{r}_i;\ve{R}_i) = |y_i-Y_i|, \;\;\;\mbox{or}\;\;\;
958 d_g(\ve{r}_i;\ve{R}_i) = |z_i-Z_i|.
961 It is possible to apply multiple independent flat-bottomed position
962 restraints of different geometry on one particle. For example, applying
963 a cylinder and a layer in $z$ keeps a particle within a
964 disk. Applying three layers in $x$, $y$, and $z$ keeps the particle within a cuboid.
966 In addition, it is possible to invert the restrained region with the
967 unrestrained region, leading to a potential that acts to keep the particle {\it outside} of the volume
968 defined by $\ve{R}_i$, $g$, and $r_\mathrm{fb}$. That feature is
969 switched on by defining a negative $r_\mathrm{fb}$ in the
970 topology. The following potential is used (Figure~\ref{fig:fbposres}B):
972 V_\mathrm{fb}^{\mathrm{inv}}(\ve{r}_i) = \frac{1}{2}k_\mathrm{fb}
973 [d_g(\ve{r}_i;\ve{R}_i) - |r_\mathrm{fb}|]^2\,
974 H[ -(d_g(\ve{r}_i;\ve{R}_i) - |r_\mathrm{fb}|)].
979 %\ifthenelse{\equal{\gmxlite}{1}}{}{
980 \subsection{Angle restraints\swapindexquiet{angle}{restraint}}
981 \label{subsec:anglerestraint}
982 These are used to restrain the angle between two pairs of particles
983 or between one pair of particles and the $z$-axis.
984 The functional form is similar to that of a proper dihedral.
985 For two pairs of atoms:
987 V_{ar}(\ve{r}_i,\ve{r}_j,\ve{r}_k,\ve{r}_l)
988 = k_{ar}(1 - \cos(n (\theta - \theta_0))
991 \theta = \arccos\left(\frac{\ve{r}_j -\ve{r}_i}{\|\ve{r}_j -\ve{r}_i\|}
992 \cdot \frac{\ve{r}_l -\ve{r}_k}{\|\ve{r}_l -\ve{r}_k\|} \right)
994 For one pair of atoms and the $z$-axis:
996 V_{ar}(\ve{r}_i,\ve{r}_j) = k_{ar}(1 - \cos(n (\theta - \theta_0))
999 \theta = \arccos\left(\frac{\ve{r}_j -\ve{r}_i}{\|\ve{r}_j -\ve{r}_i\|}
1000 \cdot \left( \begin{array}{c} 0 \\ 0 \\ 1 \\ \end{array} \right) \right)
1002 A multiplicity ($n$) of 2 is useful when you do not want to distinguish
1003 between parallel and anti-parallel vectors.
1004 The equilibrium angle $\theta$ should be between 0 and 180 degrees
1005 for multiplicity 1 and between 0 and 90 degrees for multiplicity 2.
1008 \subsection{Dihedral restraints\swapindexquiet{dihedral}{restraint}}
1009 \label{subsec:dihedralrestraint}
1010 These are used to restrain the dihedral angle $\phi$ defined by four particles
1011 as in an improper dihedral (sec.~\ref{sec:imp}) but with a slightly
1012 modified potential. Using:
1014 \phi' = \left(\phi-\phi_0\right) ~{\rm MOD}~ 2\pi
1016 where $\phi_0$ is the reference angle, the potential is defined as:
1018 V_{dihr}(\phi') ~=~ \left\{
1019 \begin{array}{lcllll}
1020 \half k_{dihr}(\phi'-\phi_0-\Delta\phi)^2
1021 &\mbox{for}& \phi' & > & \Delta\phi \\[1.5ex]
1022 0 &\mbox{for}& \phi' & \le & \Delta\phi \\[1.5ex]
1026 where $\Delta\phi$ is a user defined angle and $k_{dihr}$ is the force
1028 {\bf Note} that in the input in topology files, angles are given in degrees and
1029 force constants in kJ/mol/rad$^2$.
1030 %} % Brace matches ifthenelse test for gmxlite
1032 \subsection{Distance restraints\swapindexquiet{distance}{restraint}}
1033 \label{subsec:distancerestraint}
1035 add a penalty to the potential when the distance between specified
1036 pairs of atoms exceeds a threshold value. They are normally used to
1037 impose experimental restraints from, for instance, experiments in nuclear
1038 magnetic resonance (NMR), on the motion of the system. Thus, MD can be
1039 used for structure refinement using NMR data\index{nmr
1040 refinement}\index{refinement,nmr}.
1041 In {\gromacs} there are three ways to impose restraints on pairs of atoms:
1043 \item Simple harmonic restraints: use {\tt [ bonds ]} type 6
1044 %\ifthenelse{\equal{\gmxlite}{1}}
1046 {(see \secref{excl}).}
1047 \item\label{subsec:harmonicrestraint}Piecewise linear/harmonic restraints: {\tt [ bonds ]} type 10.
1048 \item Complex NMR distance restraints, optionally with pair, time and/or
1051 The last two options will be detailed now.
1053 The potential form for distance restraints is quadratic below a specified
1054 lower bound and between two specified upper bounds, and linear beyond the
1055 largest bound (see \figref{dist}).
1057 V_{dr}(r_{ij}) ~=~ \left\{
1058 \begin{array}{lcllllll}
1059 \half k_{dr}(r_{ij}-r_0)^2
1060 &\mbox{for}& & & r_{ij} & < & r_0 \\[1.5ex]
1061 0 &\mbox{for}& r_0 & \le & r_{ij} & < & r_1 \\[1.5ex]
1062 \half k_{dr}(r_{ij}-r_1)^2
1063 &\mbox{for}& r_1 & \le & r_{ij} & < & r_2 \\[1.5ex]
1064 \half k_{dr}(r_2-r_1)(2r_{ij}-r_2-r_1)
1065 &\mbox{for}& r_2 & \le & r_{ij} & &
1071 \centerline{\includegraphics[width=8cm]{plots/f-dr}}
1072 \caption{Distance Restraint potential.}
1079 \begin{array}{lcllllll}
1080 -k_{dr}(r_{ij}-r_0)\frac{\rvij}{r_{ij}}
1081 &\mbox{for}& & & r_{ij} & < & r_0 \\[1.5ex]
1082 0 &\mbox{for}& r_0 & \le & r_{ij} & < & r_1 \\[1.5ex]
1083 -k_{dr}(r_{ij}-r_1)\frac{\rvij}{r_{ij}}
1084 &\mbox{for}& r_1 & \le & r_{ij} & < & r_2 \\[1.5ex]
1085 -k_{dr}(r_2-r_1)\frac{\rvij}{r_{ij}}
1086 &\mbox{for}& r_2 & \le & r_{ij} & &
1090 For restraints not derived from NMR data, this functionality
1091 will usually suffice and a section of {\tt [ bonds ]} type 10
1092 can be used to apply individual restraints between pairs of
1093 %\ifthenelse{\equal{\gmxlite}{1}}{atoms.}{
1094 atoms, see \ssecref{topfile}.
1095 %} % Brace matches ifthenelse test for gmxlite
1096 For applying restraints derived from NMR measurements, more complex
1097 functionality might be required, which is provided through
1098 the {\tt [~distance_restraints~]} section and is described below.
1100 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1101 \subsubsection{Time averaging\swapindexquiet{time-averaged}{distance restraint}}
1102 Distance restraints based on instantaneous distances can potentially reduce
1103 the fluctuations in a molecule significantly. This problem can be overcome by restraining
1104 to a {\em time averaged} distance~\cite{Torda89}.
1105 The forces with time averaging are:
1108 \begin{array}{lcllllll}
1109 -k^a_{dr}(\bar{r}_{ij}-r_0)\frac{\rvij}{r_{ij}}
1110 &\mbox{for}& & & \bar{r}_{ij} & < & r_0 \\[1.5ex]
1111 0 &\mbox{for}& r_0 & \le & \bar{r}_{ij} & < & r_1 \\[1.5ex]
1112 -k^a_{dr}(\bar{r}_{ij}-r_1)\frac{\rvij}{r_{ij}}
1113 &\mbox{for}& r_1 & \le & \bar{r}_{ij} & < & r_2 \\[1.5ex]
1114 -k^a_{dr}(r_2-r_1)\frac{\rvij}{r_{ij}}
1115 &\mbox{for}& r_2 & \le & \bar{r}_{ij} & &
1118 where $\bar{r}_{ij}$ is given by an exponential running average with decay time $\tau$:
1120 \bar{r}_{ij} ~=~ < r_{ij}^{-3} >^{-1/3}
1123 The force constant $k^a_{dr}$ is switched on slowly to compensate for
1124 the lack of history at the beginning of the simulation:
1126 k^a_{dr} = k_{dr} \left(1-\exp\left(-\frac{t}{\tau}\right)\right)
1128 Because of the time averaging, we can no longer speak of a distance restraint
1131 This way an atom can satisfy two incompatible distance restraints
1132 {\em on average} by moving between two positions.
1133 An example would be an amino acid side-chain that is rotating around
1134 its $\chi$ dihedral angle, thereby coming close to various other groups.
1135 Such a mobile side chain can give rise to multiple NOEs that can not be
1136 fulfilled by a single structure.
1138 The computation of the time
1139 averaged distance in the {\tt mdrun} program is done in the following fashion:
1142 \overline{r^{-3}}_{ij}(0) &=& r_{ij}(0)^{-3} \\
1143 \overline{r^{-3}}_{ij}(t) &=& \overline{r^{-3}}_{ij}(t-\Delta t)~\exp{\left(-\frac{\Delta t}{\tau}\right)} + r_{ij}(t)^{-3}\left[1-\exp{\left(-\frac{\Delta t}{\tau}\right)}\right]
1144 \label{eqn:ravdisre}
1148 When a pair is within the bounds, it can still feel a force
1149 because the time averaged distance can still be beyond a bound.
1150 To prevent the protons from being pulled too close together, a mixed
1151 approach can be used. In this approach, the penalty is zero when the
1152 instantaneous distance is within the bounds, otherwise the violation is
1153 the square root of the product of the instantaneous violation and the
1154 time averaged violation:
1157 \begin{array}{lclll}
1158 k^a_{dr}\sqrt{(r_{ij}-r_0)(\bar{r}_{ij}-r_0)}\frac{\rvij}{r_{ij}}
1159 & \mbox{for} & r_{ij} < r_0 & \mbox{and} & \bar{r}_{ij} < r_0 \\[1.5ex]
1161 \mbox{min}\left(\sqrt{(r_{ij}-r_1)(\bar{r}_{ij}-r_1)},r_2-r_1\right)
1162 \frac{\rvij}{r_{ij}}
1163 & \mbox{for} & r_{ij} > r_1 & \mbox{and} & \bar{r}_{ij} > r_1 \\[1.5ex]
1168 \subsubsection{Averaging over multiple pairs\swapindexquiet{ensemble-averaged}{distance restraint}}
1170 Sometimes it is unclear from experimental data which atom pair
1171 gives rise to a single NOE, in other occasions it can be obvious that
1172 more than one pair contributes due to the symmetry of the system, {\eg} a
1173 methyl group with three protons. For such a group, it is not possible
1174 to distinguish between the protons, therefore they should all be taken into
1175 account when calculating the distance between this methyl group and another
1176 proton (or group of protons).
1177 Due to the physical nature of magnetic resonance, the intensity of the
1178 NOE signal is inversely proportional to the sixth power of the inter-atomic
1180 Thus, when combining atom pairs,
1181 a fixed list of $N$ restraints may be taken together,
1182 where the apparent ``distance'' is given by:
1184 r_N(t) = \left [\sum_{n=1}^{N} \bar{r}_{n}(t)^{-6} \right]^{-1/6}
1187 where we use $r_{ij}$ or \eqnref{rav} for the $\bar{r}_{n}$.
1188 The $r_N$ of the instantaneous and time-averaged distances
1189 can be combined to do a mixed restraining, as indicated above.
1190 As more pairs of protons contribute to the same NOE signal, the intensity
1191 will increase, and the summed ``distance'' will be shorter than any of
1192 its components due to the reciprocal summation.
1194 There are two options for distributing the forces over the atom pairs.
1195 In the conservative option, the force is defined as the derivative of the
1196 restraint potential with respect to the coordinates. This results in
1197 a conservative potential when time averaging is not used.
1198 The force distribution over the pairs is proportional to $r^{-6}$.
1199 This means that a close pair feels a much larger force than a distant pair,
1200 which might lead to a molecule that is ``too rigid.''
1201 The other option is an equal force distribution. In this case each pair
1202 feels $1/N$ of the derivative of the restraint potential with respect to
1203 $r_N$. The advantage of this method is that more conformations might be
1204 sampled, but the non-conservative nature of the forces can lead to
1205 local heating of the protons.
1207 It is also possible to use {\em ensemble averaging} using multiple
1208 (protein) molecules. In this case the bounds should be lowered as in:
1211 r_1 &~=~& r_1 * M^{-1/6} \\
1212 r_2 &~=~& r_2 * M^{-1/6}
1215 where $M$ is the number of molecules. The {\gromacs} preprocessor {\tt grompp}
1216 can do this automatically when the appropriate option is given.
1217 The resulting ``distance'' is
1218 then used to calculate the scalar force according to:
1222 ~& 0 \hspace{4cm} & r_{N} < r_1 \\
1223 & k_{dr}(r_{N}-r_1)\frac{\rvij}{r_{ij}} & r_1 \le r_{N} < r_2 \\
1224 & k_{dr}(r_2-r_1)\frac{\rvij}{r_{ij}} & r_{N} \ge r_2
1227 where $i$ and $j$ denote the atoms of all the
1228 pairs that contribute to the NOE signal.
1230 \subsubsection{Using distance restraints}
1232 A list of distance restrains based on NOE data can be added to a molecule
1233 definition in your topology file, like in the following example:
1236 [ distance_restraints ]
1237 ; ai aj type index type' low up1 up2 fac
1238 10 16 1 0 1 0.0 0.3 0.4 1.0
1239 10 28 1 1 1 0.0 0.3 0.4 1.0
1240 10 46 1 1 1 0.0 0.3 0.4 1.0
1241 16 22 1 2 1 0.0 0.3 0.4 2.5
1242 16 34 1 3 1 0.0 0.5 0.6 1.0
1245 In this example a number of features can be found. In columns {\tt
1246 ai} and {\tt aj} you find the atom numbers of the particles to be
1247 restrained. The {\tt type} column should always be 1. As explained in
1248 ~\ssecref{distancerestraint}, multiple distances can contribute to a single NOE
1249 signal. In the topology this can be set using the {\tt index}
1250 column. In our example, the restraints 10-28 and 10-46 both have index
1251 1, therefore they are treated simultaneously. An extra requirement
1252 for treating restraints together is that the restraints must be on
1253 successive lines, without any other intervening restraint. The {\tt
1254 type'} column will usually be 1, but can be set to 2 to obtain a
1255 distance restraint that will never be time- and ensemble-averaged;
1256 this can be useful for restraining hydrogen bonds. The columns {\tt
1257 low}, {\tt up1}, and {\tt up2} hold the values of $r_0$, $r_1$, and
1258 $r_2$ from ~\eqnref{disre}. In some cases it can be useful to have
1259 different force constants for some restraints; this is controlled by
1260 the column {\tt fac}. The force constant in the parameter file is
1261 multiplied by the value in the column {\tt fac} for each restraint.
1262 %} % Brace matches ifthenelse test for gmxlite
1264 \newcommand{\SSS}{{\mathbf S}}
1265 \newcommand{\DD}{{\mathbf D}}
1266 \newcommand{\RR}{{\mathbf R}}
1268 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1269 \subsection{Orientation restraints\swapindexquiet{orientation}{restraint}}
1270 \label{subsec:orientationrestraint}
1271 This section describes how orientations between vectors,
1272 as measured in certain NMR experiments, can be calculated
1273 and restrained in MD simulations.
1274 The presented refinement methodology and a comparison of results
1275 with and without time and ensemble averaging have been
1276 published~\cite{Hess2003}.
1277 \subsubsection{Theory}
1278 In an NMR experiment, orientations of vectors can be measured when a
1279 molecule does not tumble completely isotropically in the solvent.
1280 Two examples of such orientation measurements are
1281 residual \normindex{dipolar couplings}
1282 (between two nuclei) or chemical shift anisotropies.
1283 An observable for a vector $\ve{r}_i$ can be written as follows:
1285 \delta_i = \frac{2}{3} \mbox{tr}(\SSS\DD_i)
1287 where $\SSS$ is the dimensionless order tensor of the molecule.
1288 The tensor $\DD_i$ is given by:
1291 \DD_i = \frac{c_i}{\|\ve{r}_i\|^\alpha} \left(
1293 %3 r_x r_x - \ve{r}\cdot\ve{r} & 3 r_x r_y & 3 r_x r_z \\
1294 %3 r_x r_y & 3 r_y r_y - \ve{r}\cdot\ve{r} & 3yz \\
1295 %3 r_x r_z & 3 r_y r_z & 3 r_z r_z - \ve{r}\cdot\ve{r}
1296 %\end{array} \right)
1298 3 x x - 1 & 3 x y & 3 x z \\
1299 3 x y & 3 y y - 1 & 3 y z \\
1300 3 x z & 3 y z & 3 z z - 1 \\
1305 x=\frac{r_{i,x}}{\|\ve{r}_i\|}, \quad
1306 y=\frac{r_{i,y}}{\|\ve{r}_i\|}, \quad
1307 z=\frac{r_{i,z}}{\|\ve{r}_i\|}
1309 For a dipolar coupling $\ve{r}_i$ is the vector connecting the two
1310 nuclei, $\alpha=3$ and the constant $c_i$ is given by:
1312 c_i = \frac{\mu_0}{4\pi} \gamma_1^i \gamma_2^i \frac{\hbar}{4\pi}
1314 where $\gamma_1^i$ and $\gamma_2^i$ are the gyromagnetic ratios of the
1317 The order tensor is symmetric and has trace zero. Using a rotation matrix
1318 ${\mathbf T}$ it can be transformed into the following form:
1320 {\mathbf T}^T \SSS {\mathbf T} = s \left( \begin{array}{ccc}
1321 -\frac{1}{2}(1-\eta) & 0 & 0 \\
1322 0 & -\frac{1}{2}(1+\eta) & 0 \\
1326 where $-1 \leq s \leq 1$ and $0 \leq \eta \leq 1$.
1327 $s$ is called the order parameter and $\eta$ the asymmetry of the
1328 order tensor $\SSS$. When the molecule tumbles isotropically in the
1329 solvent, $s$ is zero, and no orientational effects can be observed
1330 because all $\delta_i$ are zero.
1334 \subsubsection{Calculating orientations in a simulation}
1335 For reasons which are explained below, the $\DD$ matrices are calculated
1336 which respect to a reference orientation of the molecule. The orientation
1337 is defined by a rotation matrix $\RR$, which is needed to least-squares fit
1338 the current coordinates of a selected set of atoms onto
1339 a reference conformation. The reference conformation is the starting
1340 conformation of the simulation. In case of ensemble averaging, which will
1341 be treated later, the structure is taken from the first subsystem.
1342 The calculated $\DD_i^c$ matrix is given by:
1345 \DD_i^c(t) = \RR(t) \DD_i(t) \RR^T(t)
1347 The calculated orientation for vector $i$ is given by:
1349 \delta^c_i(t) = \frac{2}{3} \mbox{tr}(\SSS(t)\DD_i^c(t))
1351 The order tensor $\SSS(t)$ is usually unknown.
1352 A reasonable choice for the order tensor is the tensor
1353 which minimizes the (weighted) mean square difference between the calculated
1354 and the observed orientations:
1357 MSD(t) = \left(\sum_{i=1}^N w_i\right)^{-1} \sum_{i=1}^N w_i (\delta_i^c (t) -\delta_i^{exp})^2
1359 To properly combine different types of measurements, the unit of $w_i$ should
1360 be such that all terms are dimensionless. This means the unit of $w_i$
1361 is the unit of $\delta_i$ to the power $-2$.
1362 {\bf Note} that scaling all $w_i$ with a constant factor does not influence
1365 \subsubsection{Time averaging}
1366 Since the tensors $\DD_i$ fluctuate rapidly in time, much faster than can
1367 be observed in an experiment, they should be averaged over time in the simulation.
1368 However, in a simulation the time and the number of copies of
1369 a molecule are limited. Usually one can not obtain a converged average
1370 of the $\DD_i$ tensors over all orientations of the molecule.
1371 If one assumes that the average orientations of the $\ve{r}_i$ vectors
1372 within the molecule converge much faster than the tumbling time of
1373 the molecule, the tensor can be averaged in an axis system that
1374 rotates with the molecule, as expressed by equation~(\ref{D_rot}).
1375 The time-averaged tensors are calculated
1376 using an exponentially decaying memory function:
1378 \DD^a_i(t) = \frac{\displaystyle
1379 \int_{u=t_0}^t \DD^c_i(u) \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
1381 \int_{u=t_0}^t \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
1384 Assuming that the order tensor $\SSS$ fluctuates slower than the
1385 $\DD_i$, the time-averaged orientation can be calculated as:
1387 \delta_i^a(t) = \frac{2}{3} \mbox{tr}(\SSS(t) \DD_i^a(t))
1389 where the order tensor $\SSS(t)$ is calculated using expression~(\ref{S_msd})
1390 with $\delta_i^c(t)$ replaced by $\delta_i^a(t)$.
1392 \subsubsection{Restraining}
1393 The simulated structure can be restrained by applying a force proportional
1394 to the difference between the calculated and the experimental orientations.
1395 When no time averaging is applied, a proper potential can be defined as:
1397 V = \frac{1}{2} k \sum_{i=1}^N w_i (\delta_i^c (t) -\delta_i^{exp})^2
1399 where the unit of $k$ is the unit of energy.
1400 Thus the effective force constant for restraint $i$ is $k w_i$.
1401 The forces are given by minus the gradient of $V$.
1402 The force $\ve{F}\!_i$ working on vector $\ve{r}_i$ is:
1405 & = & - \frac{\mbox{d} V}{\mbox{d}\ve{r}_i} \\
1406 & = & -k w_i (\delta_i^c (t) -\delta_i^{exp}) \frac{\mbox{d} \delta_i (t)}{\mbox{d}\ve{r}_i} \\
1407 & = & -k w_i (\delta_i^c (t) -\delta_i^{exp})
1408 \frac{2 c_i}{\|\ve{r}\|^{2+\alpha}} \left(2 \RR^T \SSS \RR \ve{r}_i - \frac{2+\alpha}{\|\ve{r}\|^2} \mbox{tr}(\RR^T \SSS \RR \ve{r}_i \ve{r}_i^T) \ve{r}_i \right)
1411 \subsubsection{Ensemble averaging}
1412 Ensemble averaging can be applied by simulating a system of $M$ subsystems
1413 that each contain an identical set of orientation restraints. The systems only
1414 interact via the orientation restraint potential which is defined as:
1416 V = M \frac{1}{2} k \sum_{i=1}^N w_i
1417 \langle \delta_i^c (t) -\delta_i^{exp} \rangle^2
1419 The force on vector $\ve{r}_{i,m}$ in subsystem $m$ is given by:
1421 \ve{F}\!_{i,m}(t) = - \frac{\mbox{d} V}{\mbox{d}\ve{r}_{i,m}} =
1422 -k w_i \langle \delta_i^c (t) -\delta_i^{exp} \rangle \frac{\mbox{d} \delta_{i,m}^c (t)}{\mbox{d}\ve{r}_{i,m}} \\
1425 \subsubsection{Time averaging}
1426 When using time averaging it is not possible to define a potential.
1427 We can still define a quantity that gives a rough idea of the energy
1428 stored in the restraints:
1430 V = M \frac{1}{2} k^a \sum_{i=1}^N w_i
1431 \langle \delta_i^a (t) -\delta_i^{exp} \rangle^2
1433 The force constant $k_a$ is switched on slowly to compensate for the
1434 lack of history at times close to $t_0$. It is exactly proportional
1435 to the amount of average that has been accumulated:
1438 k \, \frac{1}{\tau}\int_{u=t_0}^t \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
1440 What really matters is the definition of the force. It is chosen to
1441 be proportional to the square root of the product of the time-averaged
1442 and the instantaneous deviation.
1443 Using only the time-averaged deviation induces large oscillations.
1444 The force is given by:
1447 %\left\{ \begin{array}{ll}
1448 %0 & \mbox{for} \quad \langle \delta_i^a (t) -\delta_i^{exp} \rangle \langle \delta_i (t) -\delta_i^{exp} \rangle \leq 0 \\
1449 %... & \mbox{for} \quad \langle \delta_i^a (t) -\delta_i^{exp} \rangle \langle \delta_i (t) -\delta_i^{exp} \rangle > 0
1452 \left\{ \begin{array}{ll}
1453 0 & \quad \mbox{for} \quad a\, b \leq 0 \\
1455 k^a w_i \frac{a}{|a|} \sqrt{a\, b} \, \frac{\mbox{d} \delta_{i,m}^c (t)}{\mbox{d}\ve{r}_{i,m}}
1456 & \quad \mbox{for} \quad a\, b > 0
1461 a &=& \langle \delta_i^a (t) -\delta_i^{exp} \rangle \\
1462 b &=& \langle \delta_i^c (t) -\delta_i^{exp} \rangle
1465 \subsubsection{Using orientation restraints}
1466 Orientation restraints can be added to a molecule definition in
1467 the topology file in the section {\tt [~orientation_restraints~]}.
1468 Here we give an example section containing five N-H residual dipolar
1469 coupling restraints:
1472 [ orientation_restraints ]
1473 ; ai aj type exp. label alpha const. obs. weight
1475 31 32 1 1 3 3 6.083 -6.73 1.0
1476 43 44 1 1 4 3 6.083 -7.87 1.0
1477 55 56 1 1 5 3 6.083 -7.13 1.0
1478 65 66 1 1 6 3 6.083 -2.57 1.0
1479 73 74 1 1 7 3 6.083 -2.10 1.0
1482 The unit of the observable is Hz, but one can choose any other unit.
1484 ai} and {\tt aj} you find the atom numbers of the particles to be
1485 restrained. The {\tt type} column should always be 1.
1486 The {\tt exp.} column denotes the experiment number, starting
1487 at 1. For each experiment a separate order tensor $\SSS$
1488 is optimized. The label should be a unique number larger than zero
1489 for each restraint. The {\tt alpha} column contains the power $\alpha$
1490 that is used in equation~(\ref{orient_def}) to calculate the orientation.
1491 The {\tt const.} column contains the constant $c_i$ used in the same
1492 equation. The constant should have the unit of the observable times
1493 nm$^\alpha$. The column {\tt obs.} contains the observable, in any
1494 unit you like. The last column contains the weights $w_i$; the unit
1495 should be the inverse of the square of the unit of the observable.
1497 Some parameters for orientation restraints can be specified in the
1498 {\tt grompp.mdp} file, for a study of the effect of different
1499 force constants and averaging times and ensemble averaging see~\cite{Hess2003}.
1500 %} % Brace matches ifthenelse test for gmxlite
1502 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1503 \section{Polarization}
1504 Polarization can be treated by {\gromacs} by attaching
1505 \normindex{shell} (\normindex{Drude}) particles to atoms and/or
1506 virtual sites. The energy of the shell particle is then minimized at
1507 each time step in order to remain on the Born-Oppenheimer surface.
1509 \subsection{Simple polarization}
1510 This is merely a harmonic potential with equilibrium distance 0.
1512 \subsection{Water polarization}
1513 A special potential for water that allows anisotropic polarization of
1514 a single shell particle~\cite{Maaren2001a}.
1516 \subsection{Thole polarization}
1517 Based on early work by \normindex{Thole}~\cite{Thole81}, Roux and
1518 coworkers have implemented potentials for molecules like
1519 ethanol~\cite{Lamoureux2003a,Lamoureux2003b,Noskov2005a}. Within such
1520 molecules, there are intra-molecular interactions between shell
1521 particles, however these must be screened because full Coulomb would
1522 be too strong. The potential between two shell particles $i$ and $j$ is:
1523 \newcommand{\rbij}{\bar{r}_{ij}}
1525 V_{thole} ~=~ \frac{q_i q_j}{r_{ij}}\left[1-\left(1+\frac{\rbij}{2}\right){\rm exp}^{-\rbij}\right]
1527 {\bf Note} that there is a sign error in Equation~1 of Noskov {\em et al.}~\cite{Noskov2005a}:
1529 \rbij ~=~ a\frac{r_{ij}}{(\alpha_i \alpha_j)^{1/6}}
1531 where $a$ is a magic (dimensionless) constant, usually chosen to be
1532 2.6~\cite{Noskov2005a}; $\alpha_i$ and $\alpha_j$ are the polarizabilities
1533 of the respective shell particles.
1535 %} % Brace matches ifthenelse test for gmxlite
1537 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1538 \section{Free energy interactions}
1540 \index{free energy interactions}
1541 \newcommand{\LAM}{\lambda}
1542 \newcommand{\LL}{(1-\LAM)}
1543 \newcommand{\dvdl}[1]{\frac{\partial #1}{\partial \LAM}}
1544 This section describes the $\lambda$-dependence of the potentials
1545 used for free energy calculations (see \secref{fecalc}).
1546 All common types of potentials and constraints can be
1547 interpolated smoothly from state A ($\lambda=0$) to state B
1548 ($\lambda=1$) and vice versa.
1549 All bonded interactions are interpolated by linear interpolation
1550 of the interaction parameters. Non-bonded interactions can be
1551 interpolated linearly or via soft-core interactions.
1553 Starting in {\gromacs} 4.6, $\lambda$ is a vector, allowing different
1554 components of the free energy transformation to be carried out at
1555 different rates. Coulomb, Lennard-Jones, bonded, and restraint terms
1556 can all be controlled independently, as described in the {\tt .mdp}
1559 \subsubsection{Harmonic potentials}
1560 The example given here is for the bond potential, which is harmonic
1561 in {\gromacs}. However, these equations apply to the angle potential
1562 and the improper dihedral potential as well.
1564 V_b &=&\half\left[\LL k_b^A +
1565 \LAM k_b^B\right] \left[b - \LL b_0^A - \LAM b_0^B\right]^2 \\
1566 \dvdl{V_b}&=&\half(k_b^B-k_b^A)
1567 \left[b - \LL b_0^A + \LAM b_0^B\right]^2 +
1569 & & \phantom{\half}(b_0^A-b_0^B) \left[b - \LL b_0^A -\LAM b_0^B\right]
1570 \left[\LL k_b^A + \LAM k_b^B \right]
1573 \subsubsection{\gromosv{96} bonds and angles}
1574 Fourth-power bond stretching and cosine-based angle potentials
1575 are interpolated by linear interpolation of the force constant
1576 and the equilibrium position. Formulas are not given here.
1578 \subsubsection{Proper dihedrals}
1579 For the proper dihedrals, the equations are somewhat more complicated:
1581 V_d &=&\left[\LL k_d^A + \LAM k_d^B \right]
1582 \left( 1+ \cos\left[n_{\phi} \phi -
1583 \LL \phi_s^A - \LAM \phi_s^B
1585 \dvdl{V_d}&=&(k_d^B-k_d^A)
1588 n_{\phi} \phi- \LL \phi_s^A - \LAM \phi_s^B
1592 &&(\phi_s^B - \phi_s^A) \left[\LL k_d^A - \LAM k_d^B\right]
1593 \sin\left[ n_{\phi}\phi - \LL \phi_s^A - \LAM \phi_s^B \right]
1595 {\bf Note:} that the multiplicity $n_{\phi}$ can not be parameterized
1596 because the function should remain periodic on the interval $[0,2\pi]$.
1598 \subsubsection{Tabulated bonded interactions}
1599 For tabulated bonded interactions only the force constant can interpolated:
1601 V &=& (\LL k^A + \LAM k^B) \, f \\
1602 \dvdl{V} &=& (k^B - k^A) \, f
1605 \subsubsection{Coulomb interaction}
1606 The \normindex{Coulomb} interaction between two particles
1607 of which the charge varies with $\LAM$ is:
1609 V_c &=& \frac{f}{\epsrf \rij}\left[\LL q_i^A q_j^A + \LAM\, q_i^B q_j^B\right] \\
1610 \dvdl{V_c}&=& \frac{f}{\epsrf \rij}\left[- q_i^A q_j^A + q_i^B q_j^B\right]
1612 where $f = \frac{1}{4\pi \varepsilon_0} = 138.935\,485$ (see \chref{defunits}).
1614 \subsubsection{Coulomb interaction with \normindex{reaction field}}
1615 The Coulomb interaction including a reaction field, between two particles
1616 of which the charge varies with $\LAM$ is:
1618 V_c &=& f\left[\frac{1}{\rij} + k_{rf}~ \rij^2 -c_{rf}\right]
1619 \left[\LL q_i^A q_j^A + \LAM\, q_i^B q_j^B\right] \\
1620 \dvdl{V_c}&=& f\left[\frac{1}{\rij} + k_{rf}~ \rij^2 -c_{rf}\right]
1621 \left[- q_i^A q_j^A + q_i^B q_j^B\right]
1622 \label{eq:dVcoulombdlambda}
1624 {\bf Note} that the constants $k_{rf}$ and $c_{rf}$ are
1625 defined using the dielectric
1626 constant $\epsrf$ of the medium (see \secref{coulrf}).
1628 \subsubsection{Lennard-Jones interaction}
1629 For the \normindex{Lennard-Jones} interaction between two particles
1630 of which the {\em atom type} varies with $\LAM$ we can write:
1632 V_{LJ} &=& \frac{\LL C_{12}^A + \LAM\, C_{12}^B}{\rij^{12}} -
1633 \frac{\LL C_6^A + \LAM\, C_6^B}{\rij^6} \\
1634 \dvdl{V_{LJ}}&=&\frac{C_{12}^B - C_{12}^A}{\rij^{12}} -
1635 \frac{C_6^B - C_6^A}{\rij^6}
1636 \label{eq:dVljdlambda}
1638 It should be noted that it is also possible to express a pathway from
1639 state A to state B using $\sigma$ and $\epsilon$ (see \eqnref{sigeps}).
1640 It may seem to make sense physically to vary the force field parameters
1641 $\sigma$ and $\epsilon$ rather
1642 than the derived parameters $C_{12}$ and $C_{6}$.
1643 However, the difference between the pathways in parameter space
1644 is not large, and the free energy itself
1645 does not depend on the pathway, so we use the simple formulation
1648 \subsubsection{Kinetic Energy}
1649 When the mass of a particle changes, there is also a contribution of
1650 the kinetic energy to the free energy (note that we can not write
1651 the momentum \ve{p} as m\ve{v}, since that would result
1652 in the sign of $\dvdl{E_k}$ being incorrect~\cite{Gunsteren98a}):
1655 E_k &=& \half\frac{\ve{p}^2}{\LL m^A + \LAM m^B} \\
1656 \dvdl{E_k}&=& -\half\frac{\ve{p}^2(m^B-m^A)}{(\LL m^A + \LAM m^B)^2}
1658 after taking the derivative, we {\em can} insert \ve{p} = m\ve{v}, such that:
1660 \dvdl{E_k}~=~ -\half\ve{v}^2(m^B-m^A)
1663 \subsubsection{Constraints}
1664 \label{subsubsec:constraints}
1665 \newcommand{\clam}{C_{\lambda}}
1666 The constraints are formally part of the Hamiltonian, and therefore
1667 they give a contribution to the free energy. In {\gromacs} this can be
1668 calculated using the \normindex{LINCS} or the \normindex{SHAKE} algorithm.
1669 If we have a number of constraint equations $g_k$:
1671 g_k = \ve{r}_{k} - d_{k}
1673 where $\ve{r}_k$ is the distance vector between two particles and
1674 $d_k$ is the constraint distance between the two particles, we can write
1675 this using a $\LAM$-dependent distance as
1677 g_k = \ve{r}_{k} - \left(\LL d_{k}^A + \LAM d_k^B\right)
1679 the contribution $\clam$
1680 to the Hamiltonian using Lagrange multipliers $\lambda$:
1682 \clam &=& \sum_k \lambda_k g_k \\
1683 \dvdl{\clam} &=& \sum_k \lambda_k \left(d_k^B-d_k^A\right)
1687 \subsection{Soft-core interactions\index{soft-core interactions}}
1689 \centerline{\includegraphics[height=6cm]{plots/softcore}}
1690 \caption{Soft-core interactions at $\LAM=0.5$, with $p=2$ and
1691 $C_6^A=C_{12}^A=C_6^B=C_{12}^B=1$.}
1692 \label{fig:softcore}
1694 In a free-energy calculation where particles grow out of nothing, or
1695 particles disappear, using the the simple linear interpolation of the
1696 Lennard-Jones and Coulomb potentials as described in Equations~\ref{eq:dVljdlambda}
1697 and \ref{eq:dVcoulombdlambda} may lead to poor convergence. When the particles have nearly disappeared, or are close to appearing (at $\LAM$ close to 0 or 1), the interaction energy will be weak enough for particles to get very
1698 close to each other, leading to large fluctuations in the measured values of
1699 $\partial V/\partial \LAM$ (which, because of the simple linear
1700 interpolation, depends on the potentials at both the endpoints of $\LAM$).
1702 To circumvent these problems, the singularities in the potentials need to be removed. This can be done by modifying the regular Lennard-Jones and Coulomb potentials with ``soft-core'' potentials that limit the energies and forces
1703 involved at $\LAM$ values between 0 and 1, but not \emph{at} $\LAM=0$
1706 In {\gromacs} the soft-core potentials $V_{sc}$ are shifted versions of the
1707 regular potentials, so that the singularity in the potential and its
1708 derivatives at $r=0$ is never reached:
1710 V_{sc}(r) &=& \LL V^A(r_A) + \LAM V^B(r_B)
1712 r_A &=& \left(\alpha \sigma_A^6 \LAM^p + r^6 \right)^\frac{1}{6}
1714 r_B &=& \left(\alpha \sigma_B^6 \LL^p + r^6 \right)^\frac{1}{6}
1716 where $V^A$ and $V^B$ are the normal ``hard core'' Van der Waals or
1717 electrostatic potentials in state A ($\LAM=0$) and state B ($\LAM=1$)
1718 respectively, $\alpha$ is the soft-core parameter (set with {\tt sc_alpha}
1719 in the {\tt .mdp} file), $p$ is the soft-core $\LAM$ power (set with
1720 {\tt sc_power}), $\sigma$ is the radius of the interaction, which is
1721 $(C_{12}/C_6)^{1/6}$ or an input parameter ({\tt sc_sigma}) when $C_6$
1722 or $C_{12}$ is zero.
1724 For intermediate $\LAM$, $r_A$ and $r_B$ alter the interactions very little
1725 for $r > \alpha^{1/6} \sigma$ and quickly switch the soft-core
1726 interaction to an almost constant value for smaller $r$ (\figref{softcore}).
1729 F_{sc}(r) = -\frac{\partial V_{sc}(r)}{\partial r} =
1730 \LL F^A(r_A) \left(\frac{r}{r_A}\right)^5 +
1731 \LAM F^B(r_B) \left(\frac{r}{r_B}\right)^5
1733 where $F^A$ and $F^B$ are the ``hard core'' forces.
1734 The contribution to the derivative of the free energy is:
1736 \dvdl{V_{sc}(r)} & = &
1737 V^B(r_B) -V^A(r_A) +
1738 \LL \frac{\partial V^A(r_A)}{\partial r_A}
1739 \frac{\partial r_A}{\partial \LAM} +
1740 \LAM\frac{\partial V^B(r_B)}{\partial r_B}
1741 \frac{\partial r_B}{\partial \LAM}
1744 V^B(r_B) -V^A(r_A) + \nonumber \\
1747 \left[ \LAM F^B(r_B) r^{-5}_B \sigma_B^6 \LL^{p-1} -
1748 \LL F^A(r_A) r^{-5}_A \sigma_A^6 \LAM^{p-1} \right]
1751 The original GROMOS Lennard-Jones soft-core function~\cite{Beutler94}
1752 uses $p=2$, but $p=1$ gives a smoother $\partial H/\partial\LAM$ curve.
1753 %When the changes between the two states involve both the disappearing
1754 %and appearing of atoms, it is important that the overlapping of atoms
1755 %happens around $\LAM=0.5$. This can usually be achieved with
1756 %$\alpha$$\approx0.7$ for $p=1$ and $\alpha$$\approx1.5$ for $p=2$.
1757 %MRS: this is now eliminated as of 4.6, since changes between atoms are done linearly.
1759 Another issue that should be considered is the soft-core effect of hydrogens
1760 without Lennard-Jones interaction. Their soft-core $\sigma$ is
1761 set with {\tt sc-sigma} in the {\tt .mdp} file. These hydrogens
1762 produce peaks in $\partial H/\partial\LAM$ at $\LAM$ is 0 and/or 1 for $p=1$
1763 and close to 0 and/or 1 with $p=2$. Lowering {\tt\mbox{sc-sigma}} will decrease
1764 this effect, but it will also increase the interactions with hydrogens
1765 relative to the other interactions in the soft-core state.
1767 When soft core potentials are selected (by setting {\tt sc-alpha} \textgreater
1768 0), and the Coulomb and Lennard-Jones potentials are turned on or off
1769 sequentially, then the Coulombic interaction is turned off linearly,
1770 rather than using soft core interactions, which should be less
1771 statistically noisy in most cases. This behavior can be overwritten
1772 by using the mdp option {\tt sc-coul} to 'yes'. Additionally, the
1773 soft-core interaction potential is only applied when either the A or B
1774 state has zero interaction potential. If both A and B states have
1775 nonzero interaction potential, default linear scaling described above
1776 is used. When both Coulombic and Lennard-Jones interactions are turned
1777 off simultaneously, a soft-core potential is used, and a hydrogen is
1778 being introduced or deleted, the sigma is set to {\tt sc-sigma-min},
1779 which itself defaults to {\tt sc-sigma-default}.
1781 Recently, a new formulation of the soft-core approach has been derived
1782 that in most cases gives lower and more even statistical variance than
1783 the standard soft-core path described above.~\cite{Pham2011,Pham2012}
1784 Specifically, we have:
1786 V_{sc}(r) &=& \LL V^A(r_A) + \LAM V^B(r_B)
1788 r_A &=& \left(\alpha \sigma_A^{48} \LAM^p + r^{48} \right)^\frac{1}{48}
1790 r_B &=& \left(\alpha \sigma_B^{48} \LL^p + r^{48} \right)^\frac{1}{48}
1792 This ``1-1-48'' path is also implemented in {\gromacs}. Note that for this path the soft core $\alpha$
1793 should satisfy $0.001 < \alpha < 0.003$,rather than $\alpha \approx
1796 %} % Brace matches ifthenelse test for gmxlite
1798 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1800 \subsection{Exclusions and 1-4 Interactions.}
1801 Atoms within a molecule that are close by in the chain,
1802 {\ie} atoms that are covalently bonded, or linked by one or two
1803 atoms are called {\em first neighbors, second neighbors} and
1804 {\em \swapindex{third}{neighbor}s}, respectively (see \figref{chain}).
1805 Since the interactions of atom {\bf i} with atoms {\bf i+1} and {\bf i+2}
1806 are mainly quantum mechanical, they can not be modeled by a Lennard-Jones potential.
1807 Instead it is assumed that these interactions are adequately modeled
1808 by a harmonic bond term or constraint ({\bf i, i+1}) and a harmonic angle term
1809 ({\bf i, i+2}). The first and second neighbors (atoms {\bf i+1} and {\bf i+2})
1811 {\em excluded} from the Lennard-Jones \swapindex{interaction}{list}
1813 atoms {\bf i+1} and {\bf i+2} are called {\em \normindex{exclusions}} of atom {\bf i}.
1816 \centerline{\includegraphics[angle=270,width=8cm]{plots/chain}}
1817 \caption{Atoms along an alkane chain.}
1821 For third neighbors, the normal Lennard-Jones repulsion is sometimes
1822 still too strong, which means that when applied to a molecule, the
1823 molecule would deform or break due to the internal strain. This is
1824 especially the case for carbon-carbon interactions in a {\em
1825 cis}-conformation ({\eg} {\em cis}-butane). Therefore, for some of these
1826 interactions, the Lennard-Jones repulsion has been reduced in the
1827 {\gromos} force field, which is implemented by keeping a separate list of
1828 1-4 and normal Lennard-Jones parameters. In other force fields, such
1829 as OPLS~\cite{Jorgensen88}, the standard Lennard-Jones parameters are reduced
1830 by a factor of two, but in that case also the dispersion (r$^{-6}$)
1831 and the Coulomb interaction are scaled.
1832 {\gromacs} can use either of these methods.
1834 \subsection{Charge Groups\index{charge group}}
1836 In principle, the force calculation in MD is an $O(N^2)$ problem.
1837 Therefore, we apply a \normindex{cut-off} for non-bonded force (NBF)
1838 calculations; only the particles within a certain distance of each
1839 other are interacting. This reduces the cost to $O(N)$ (typically
1840 $100N$ to $200N$) of the NBF. It also introduces an error, which is,
1841 in most cases, acceptable, except when applying the cut-off implies
1842 the creation of charges, in which case you should consider using the
1843 lattice sum methods provided by {\gromacs}.
1845 Consider a water molecule interacting with another atom. When we would apply
1846 the cut-off on an atom-atom basis we might include the atom-oxygen
1847 interaction (with a charge of $-0.82$) without the compensating charge
1848 of the protons, and as a result, induce a large dipole moment over the system.
1849 Therefore, we have to keep groups of atoms with total charge
1850 0 together. These groups are called {\em charge groups}.
1852 \subsection{Treatment of Cut-offs\index{cut-off}}
1853 \newcommand{\rs}{$r_{short}$}
1854 \newcommand{\rl}{$r_{long}$}
1855 {\gromacs} is quite flexible in treating cut-offs, which implies
1856 there can be quite a number of parameters to set. These parameters are
1857 set in the input file for {\tt grompp}. There are two sort of parameters
1858 that affect the cut-off interactions; you can select which type
1859 of interaction to use in each case, and which cut-offs should be
1860 used in the neighbor searching.
1862 For both Coulomb and van der Waals interactions there are interaction
1863 type selectors (termed {\tt vdwtype} and {\tt coulombtype}) and two
1864 parameters, for a total of six non-bonded interaction parameters. See
1865 \secref{mdpopt} for a complete description of these parameters.
1867 The neighbor searching (NS) can be performed using a single-range, or a twin-range
1868 approach. Since the former is merely a special case of the latter, we will
1869 discuss the more general twin-range. In this case, NS is described by two
1870 radii: {\tt rlist} and max({\tt rcoulomb},{\tt rvdw}).
1871 Usually one builds the neighbor list every 10 time steps
1872 or every 20 fs (parameter {\tt nstlist}). In the neighbor list, all interaction
1873 pairs that fall within {\tt rlist} are stored. Furthermore, the
1874 interactions between pairs that do not
1875 fall within {\tt rlist} but do fall within max({\tt rcoulomb},{\tt rvdw})
1876 are computed during NS. The
1877 forces and energy are stored separately and added to short-range forces
1878 at every time step between successive NS. If {\tt rlist} =
1879 max({\tt rcoulomb},{\tt rvdw}), no forces
1880 are evaluated during neighbor list generation.
1881 The \normindex{virial} is calculated from the sum of the short- and
1883 This means that the virial can be slightly asymmetrical at non-NS steps.
1884 In single precision, the virial is almost always asymmetrical because the
1885 off-diagonal elements are about as large as each element in the sum.
1886 In most cases this is not really a problem, since the fluctuations in the
1887 virial can be 2 orders of magnitude larger than the average.
1889 Except for the plain cut-off,
1890 all of the interaction functions in \tabref{funcparm}
1891 require that neighbor searching be done with a larger radius than the $r_c$
1892 specified for the functional form, because of the use of charge groups.
1893 The extra radius is typically of the order of 0.25 nm (roughly the
1894 largest distance between two atoms in a charge group plus the distance a
1895 charge group can diffuse within neighbor list updates).
1897 %If your charge groups are very large it may be interesting to turn off charge
1898 %groups, by setting the option
1899 %{\tt bAtomList = yes} in your {\tt grompp.mdp} file.
1900 %In this case only a small extra radius to account for diffusion needs to be
1901 %added (0.1 nm). Do not however use this together with the plain cut-off
1902 %method, as it will generate large artifacts (\secref{cg}).
1903 %In summary, there are four parameters that describe NS behavior:
1904 %{\tt nstlist} (update frequency in number of time steps),
1905 %{\tt bAtomList} (whether or not charge groups are used to generate neighbor list, the default is to use charge groups, so {\tt bAtomList = no}),
1906 %{\tt rshort} and {\tt rlong} which are the two radii {\rs} and {\rl}
1911 \begin{tabular}{|ll|l|}
1913 \multicolumn{2}{|c|}{Type} & Parameters \\
1915 Coulomb&Plain cut-off & $r_c$, $\epsr$ \\
1916 &Reaction field & $r_c$, $\epsrf$ \\
1917 &Shift function & $r_1$, $r_c$, $\epsr$ \\
1918 &Switch function & $r_1$, $r_c$, $\epsr$ \\
1920 VdW&Plain cut-off & $r_c$ \\
1921 &Shift function & $r_1$, $r_c$ \\
1922 &Switch function & $r_1$, $r_c$ \\
1925 \caption[Parameters for the different functional forms of the
1926 non-bonded interactions.]{Parameters for the different functional
1927 forms of the non-bonded interactions.}
1928 \label{tab:funcparm}
1930 %} % Brace matches ifthenelse test for gmxlite
1933 \newcommand{\vvis}{\ve{r}_s}
1934 \newcommand{\Fi}{\ve{F}_i'}
1935 \newcommand{\Fj}{\ve{F}_j'}
1936 \newcommand{\Fk}{\ve{F}_k'}
1937 \newcommand{\Fl}{\ve{F}_l'}
1938 \newcommand{\Fvis}{\ve{F}_{s}}
1939 \newcommand{\rvik}{\ve{r}_{ik}}
1940 \newcommand{\rvis}{\ve{r}_{is}}
1941 \newcommand{\rvjk}{\ve{r}_{jk}}
1942 \newcommand{\rvjl}{\ve{r}_{jl}}
1944 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1945 \section{Virtual interaction sites\index{virtual interaction sites}}
1946 \label{sec:virtual_sites}
1947 Virtual interaction sites (called \seeindex{dummy atoms}{virtual interaction sites} in {\gromacs} versions before 3.3)
1948 can be used in {\gromacs} in a number of ways.
1949 We write the position of the virtual site $\ve{r}_s$ as a function of
1950 the positions of other particles \ve{r}$_i$: $\ve{r}_s =
1951 f(\ve{r}_1..\ve{r}_n)$. The virtual site, which may carry charge or be
1952 involved in other interactions, can now be used in the force
1953 calculation. The force acting on the virtual site must be
1954 redistributed over the particles with mass in a consistent way.
1955 A good way to do this can be found in ref.~\cite{Berendsen84b}.
1956 We can write the potential energy as:
1958 V = V(\vvis,\ve{r}_1,\ldots,\ve{r}_n) = V^*(\ve{r}_1,\ldots,\ve{r}_n)
1960 The force on the particle $i$ is then:
1962 \ve{F}_i = -\frac{\partial V^*}{\partial \ve{r}_i}
1963 = -\frac{\partial V}{\partial \ve{r}_i} -
1964 \frac{\partial V}{\partial \vvis}
1965 \frac{\partial \vvis}{\partial \ve{r}_i}
1966 = \ve{F}_i^{direct} + \Fi
1968 The first term is the normal force.
1969 The second term is the force on particle $i$ due to the virtual site, which
1970 can be written in tensor notation:
1971 \newcommand{\partd}[2]{\displaystyle\frac{\partial #1}{\partial #2_i}}
1973 \Fi = \left[\begin{array}{ccc}
1974 \partd{x_s}{x} & \partd{y_s}{x} & \partd{z_s}{x} \\[1ex]
1975 \partd{x_s}{y} & \partd{y_s}{y} & \partd{z_s}{y} \\[1ex]
1976 \partd{x_s}{z} & \partd{y_s}{z} & \partd{z_s}{z}
1977 \end{array}\right]\Fvis
1980 where $\Fvis$ is the force on the virtual site and $x_s$, $y_s$ and
1981 $z_s$ are the coordinates of the virtual site. In this way, the total
1982 force and the total torque are conserved~\cite{Berendsen84b}.
1984 The computation of the \normindex{virial}
1985 (\eqnref{Xi}) is non-trivial when virtual sites are used. Since the
1986 virial involves a summation over all the atoms (rather than virtual
1987 sites), the forces must be redistributed from the virtual sites to the
1988 atoms (using ~\eqnref{fvsite}) {\em before} computation of the
1989 virial. In some special cases where the forces on the atoms can be
1990 written as a linear combination of the forces on the virtual sites (types 2
1991 and 3 below) there is no difference between computing the virial
1992 before and after the redistribution of forces. However, in the
1993 general case redistribution should be done first.
1996 \centerline{\includegraphics[width=15cm]{plots/dummies}}
1997 \caption[Virtual site construction.]{The six different types of virtual
1998 site construction in \protect{\gromacs}. The constructing atoms are
1999 shown as black circles, the virtual sites in gray.}
2003 There are six ways to construct virtual sites from surrounding atoms in
2004 {\gromacs}, which we classify by the number of constructing
2005 atoms. {\bf Note} that all site types mentioned can be constructed from
2006 types 3fd (normalized, in-plane) and 3out (non-normalized, out of
2007 plane). However, the amount of computation involved increases sharply
2008 along this list, so we strongly recommended using the first adequate
2009 virtual site type that will be sufficient for a certain purpose.
2010 \figref{vsites} depicts 6 of the available virtual site constructions.
2011 The conceptually simplest construction types are linear combinations:
2013 \vvis = \sum_{i=1}^N w_i \, \ve{r}_i
2015 The force is then redistributed using the same weights:
2020 The types of virtual sites supported in {\gromacs} are given in the list below.
2021 Constructing atoms in virtual sites can be virtual sites themselves, but
2022 only if they are higher in the list, i.e. virtual sites can be
2023 constructed from ``particles'' that are simpler virtual sites.
2025 \item[{\bf\sf 2.}]\label{subsec:vsite2}As a linear combination of two atoms
2026 (\figref{vsites} 2):
2028 w_i = 1 - a ~,~~ w_j = a
2030 In this case the virtual site is on the line through atoms $i$ and
2033 \item[{\bf\sf 3.}]\label{subsec:vsite3}As a linear combination of three atoms
2034 (\figref{vsites} 3):
2036 w_i = 1 - a - b ~,~~ w_j = a ~,~~ w_k = b
2038 In this case the virtual site is in the plane of the other three
2041 \item[{\bf\sf 3fd.}]\label{subsec:vsite3fd}In the plane of three atoms, with a fixed distance
2042 (\figref{vsites} 3fd):
2044 \vvis ~=~ \ve{r}_i + b \frac{ \rvij + a \rvjk }
2045 {| \rvij + a \rvjk |}
2047 In this case the virtual site is in the plane of the other three
2048 particles at a distance of $|b|$ from $i$.
2049 The force on particles $i$, $j$ and $k$ due to the force on the virtual
2050 site can be computed as:
2053 \Fi &=& \displaystyle \Fvis - \gamma ( \Fvis - \ve{p} ) \\[1ex]
2054 \Fj &=& \displaystyle (1-a)\gamma (\Fvis - \ve{p}) \\[1ex]
2055 \Fk &=& \displaystyle a \gamma (\Fvis - \ve{p}) \\
2059 \displaystyle \gamma = \frac{b}{| \rvij + a \rvjk |} \\[2ex]
2060 \displaystyle \ve{p} = \frac{ \rvis \cdot \Fvis }
2061 { \rvis \cdot \rvis } \rvis
2065 \item[{\bf\sf 3fad.}]\label{subsec:vsite3fad}In the plane of three atoms, with a fixed angle and
2066 distance (\figref{vsites} 3fad):
2068 \label{eqn:vsite2fad-F}
2069 \vvis ~=~ \ve{r}_i +
2070 d \cos \theta \frac{\rvij}{|\rvij|} +
2071 d \sin \theta \frac{\ve{r}_\perp}{|\ve{r}_\perp|}
2073 \ve{r}_\perp ~=~ \rvjk -
2074 \frac{ \rvij \cdot \rvjk }
2075 { \rvij \cdot \rvij }
2078 In this case the virtual site is in the plane of the other three
2079 particles at a distance of $|d|$ from $i$ at an angle of
2080 $\alpha$ with $\rvij$. Atom $k$ defines the plane and the
2081 direction of the angle. {\bf Note} that in this case $b$ and
2082 $\alpha$ must be specified, instead of $a$ and $b$ (see also
2083 \secref{vsitetop}). The force on particles $i$, $j$ and $k$
2084 due to the force on the virtual site can be computed as (with
2085 $\ve{r}_\perp$ as defined in \eqnref{vsite2fad-F}):
2086 \newcommand{\dfrac}{\displaystyle\frac}
2089 \begin{array}{lclllll}
2091 \dfrac{d \cos \theta}{|\rvij|} \ve{F}_1 &+&
2092 \dfrac{d \sin \theta}{|\ve{r}_\perp|} \left(
2093 \dfrac{ \rvij \cdot \rvjk }
2094 { \rvij \cdot \rvij } \ve{F}_2 +
2095 \ve{F}_3 \right) \\[3ex]
2097 \dfrac{d \cos \theta}{|\rvij|} \ve{F}_1 &-&
2098 \dfrac{d \sin \theta}{|\ve{r}_\perp|} \left(
2100 \dfrac{ \rvij \cdot \rvjk }
2101 { \rvij \cdot \rvij } \ve{F}_2 +
2102 \ve{F}_3 \right) \\[3ex]
2104 \dfrac{d \sin \theta}{|\ve{r}_\perp|} \ve{F}_2 \\[3ex]
2108 \dfrac{ \rvij \cdot \Fvis }
2109 { \rvij \cdot \rvij } \rvij
2111 \ve{F}_2 = \ve{F}_1 -
2112 \dfrac{ \ve{r}_\perp \cdot \Fvis }
2113 { \ve{r}_\perp \cdot \ve{r}_\perp } \ve{r}_\perp
2115 \ve{F}_3 = \dfrac{ \rvij \cdot \Fvis }
2116 { \rvij \cdot \rvij } \ve{r}_\perp
2120 \item[{\bf\sf 3out.}]\label{subsec:vsite3out}As a non-linear combination of three atoms, out of plane
2121 (\figref{vsites} 3out):
2123 \vvis ~=~ \ve{r}_i + a \rvij + b \rvik +
2124 c (\rvij \times \rvik)
2126 This enables the construction of virtual sites out of the plane of the
2128 The force on particles $i,j$ and $k$ due to the force on the virtual
2129 site can be computed as:
2133 \Fj &=& \left[\begin{array}{ccc}
2134 a & -c\,z_{ik} & c\,y_{ik} \\[0.5ex]
2135 c\,z_{ik} & a & -c\,x_{ik} \\[0.5ex]
2136 -c\,y_{ik} & c\,x_{ik} & a
2137 \end{array}\right]\Fvis \\
2139 \Fk &=& \left[\begin{array}{ccc}
2140 b & c\,z_{ij} & -c\,y_{ij} \\[0.5ex]
2141 -c\,z_{ij} & b & c\,x_{ij} \\[0.5ex]
2142 c\,y_{ij} & -c\,x_{ij} & b
2143 \end{array}\right]\Fvis \\
2144 \Fi &=& \Fvis - \Fj - \Fk
2148 \item[{\bf\sf 4fdn.}]\label{subsec:vsite4fdn}From four atoms, with a fixed distance, see separate Fig.\ \ref{fig:vsite-4fdn}.
2149 This construction is a bit
2150 complex, in particular since the previous type (4fd) could be unstable which forced us
2151 to introduce a more elaborate construction:
2154 \centerline{\includegraphics[width=5cm]{plots/vsite-4fdn}}
2155 \caption {The new 4fdn virtual site construction, which is stable even when all constructing
2156 atoms are in the same plane.}
2157 \label{fig:vsite-4fdn}
2161 \mathbf{r}_{ja} &=& a\, \mathbf{r}_{ik} - \mathbf{r}_{ij} = a\, (\mathbf{x}_k - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\
2162 \mathbf{r}_{jb} &=& b\, \mathbf{r}_{il} - \mathbf{r}_{ij} = b\, (\mathbf{x}_l - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\
2163 \mathbf{r}_m &=& \mathbf{r}_{ja} \times \mathbf{r}_{jb} \nonumber \\
2164 \mathbf{x}_s &=& \mathbf{x}_i + c \frac{\mathbf{r}_m}{|\mathbf{r}_m|}
2168 In this case the virtual site is at a distance of $|c|$ from $i$, while $a$ and $b$ are
2169 parameters. {\bf Note} that the vectors $\mathbf{r}_{ik}$ and $\mathbf{r}_{ij}$ are not normalized
2170 to save floating-point operations.
2171 The force on particles $i$, $j$, $k$ and $l$ due to the force
2172 on the virtual site are computed through chain rule derivatives
2173 of the construction expression. This is exact and conserves energy,
2174 but it does lead to relatively lengthy expressions that we do not
2175 include here (over 200 floating-point operations). The interested reader can
2176 look at the source code in \verb+vsite.c+. Fortunately, this vsite type is normally
2177 only used for chiral centers such as $C_{\alpha}$ atoms in proteins.
2179 The new 4fdn construct is identified with a `type' value of 2 in the topology. The earlier 4fd
2180 type is still supported internally (`type' value 1), but it should not be used for
2181 new simulations. All current {\gromacs} tools will automatically generate type 4fdn instead.
2184 \item[{\bf\sf N.}]\label{subsec:vsiteN} A linear combination of $N$ atoms with relative
2185 weights $a_i$. The weight for atom $i$ is:
2187 w_i = a_i \left(\sum_{j=1}^N a_j \right)^{-1}
2189 There are three options for setting the weights:
2191 \item[COG] center of geometry: equal weights
2192 \item[COM] center of mass: $a_i$ is the mass of atom $i$;
2193 when in free-energy simulations the mass of the atom is changed,
2194 only the mass of the A-state is used for the weight
2195 \item[COW] center of weights: $a_i$ is defined by the user
2199 %} % Brace matches ifthenelse test for gmxlite
2201 \newcommand{\dr}{{\rm d}r}
2202 \newcommand{\avcsix}{\left< C_6 \right>}
2204 %\ifthenelse{\equal{\gmxlite}{1}}{}{
2205 \section{Long Range Electrostatics}
2206 \label{sec:lr_elstat}
2207 \subsection{Ewald summation\index{Ewald sum}}
2209 The total electrostatic energy of $N$ particles and their periodic
2210 images\index{periodic boundary conditions} is given by
2212 V=\frac{f}{2}\sum_{n_x}\sum_{n_y}
2213 \sum_{n_{z}*} \sum_{i}^{N} \sum_{j}^{N}
2214 \frac{q_i q_j}{{\bf r}_{ij,{\bf n}}}.
2215 \label{eqn:totalcoulomb}
2217 $(n_x,n_y,n_z)={\bf n}$ is the box index vector, and the star indicates that
2218 terms with $i=j$ should be omitted when $(n_x,n_y,n_z)=(0,0,0)$. The
2219 distance ${\bf r}_{ij,{\bf n}}$ is the real distance between the charges and
2220 not the minimum-image. This sum is conditionally convergent, but
2223 Ewald summation was first introduced as a method to calculate
2224 long-range interactions of the periodic images in
2225 crystals~\cite{Ewald21}. The idea is to convert the single
2226 slowly-converging sum \eqnref{totalcoulomb} into two
2227 quickly-converging terms and a constant term:
2229 V &=& V_{dir} + V_{rec} + V_{0} \\[0.5ex]
2230 V_{dir} &=& \frac{f}{2} \sum_{i,j}^{N}
2231 \sum_{n_x}\sum_{n_y}
2232 \sum_{n_{z}*} q_i q_j \frac{\mbox{erfc}(\beta {r}_{ij,{\bf n}} )}{{r}_{ij,{\bf n}}} \\[0.5ex]
2233 V_{rec} &=& \frac{f}{2 \pi V} \sum_{i,j}^{N} q_i q_j
2234 \sum_{m_x}\sum_{m_y}
2235 \sum_{m_{z}*} \frac{\exp{\left( -(\pi {\bf m}/\beta)^2 + 2 \pi i
2236 {\bf m} \cdot ({\bf r}_i - {\bf r}_j)\right)}}{{\bf m}^2} \\[0.5ex]
2237 V_{0} &=& -\frac{f \beta}{\sqrt{\pi}}\sum_{i}^{N} q_i^2,
2239 where $\beta$ is a parameter that determines the relative weight of the
2240 direct and reciprocal sums and ${\bf m}=(m_x,m_y,m_z)$.
2241 In this way we can use a short cut-off (of the order of $1$~nm) in the direct space sum and a
2242 short cut-off in the reciprocal space sum ({\eg} 10 wave vectors in each
2243 direction). Unfortunately, the computational cost of the reciprocal
2244 part of the sum increases as $N^2$
2245 (or $N^{3/2}$ with a slightly better algorithm) and it is therefore not
2246 realistic for use in large systems.
2248 \subsubsection{Using Ewald}
2249 Don't use Ewald unless you are absolutely sure this is what you want -
2250 for almost all cases the PME method below will perform much better.
2251 If you still want to employ classical Ewald summation enter this in
2252 your {\tt .mdp} file, if the side of your box is about $3$~nm:
2259 fourierspacing = 0.6
2263 The ratio of the box dimensions and the {\tt fourierspacing} parameter determines
2264 the highest magnitude of wave vectors $m_x,m_y,m_z$ to use in each
2265 direction. With a 3-nm cubic box this example would use $11$ wave vectors
2266 (from $-5$ to $5$) in each direction. The {\tt ewald-rtol} parameter
2267 is the relative strength of the electrostatic interaction at the
2268 cut-off. Decreasing this gives you a more accurate direct sum, but a
2269 less accurate reciprocal sum.
2271 \subsection{\normindex{PME}}
2273 Particle-mesh Ewald is a method proposed by Tom
2274 Darden~\cite{Darden93} to improve the performance of the
2275 reciprocal sum. Instead of directly summing wave vectors, the charges
2276 are assigned to a grid using interpolation. The implementation in
2277 {\gromacs} uses cardinal B-spline interpolation~\cite{Essmann95},
2278 which is referred to as smooth PME (SPME).
2279 The grid is then Fourier transformed with a 3D FFT algorithm and the
2280 reciprocal energy term obtained by a single sum over the grid in
2283 The potential at the grid points is calculated by inverse
2284 transformation, and by using the interpolation factors we get the
2285 forces on each atom.
2287 The PME algorithm scales as $N \log(N)$, and is substantially faster
2288 than ordinary Ewald summation on medium to large systems. On very
2289 small systems it might still be better to use Ewald to avoid the
2290 overhead in setting up grids and transforms.
2291 For the parallelization of PME see the section on MPMD PME (\ssecref{mpmd_pme}).
2293 With the Verlet cut-off scheme, the PME direct space potential is
2294 shifted by a constant such that the potential is zero at the
2295 cut-off. This shift is small and since the net system charge is close
2296 to zero, the total shift is very small, unlike in the case of the
2297 Lennard-Jones potential where all shifts add up. We apply the shift
2298 anyhow, such that the potential is the exact integral of the force.
2300 \subsubsection{Using PME}
2301 To use Particle-mesh Ewald summation in {\gromacs}, specify the
2302 following lines in your {\tt .mdp} file:
2309 fourierspacing = 0.12
2314 In this case the {\tt fourierspacing} parameter determines the maximum
2315 spacing for the FFT grid (i.e. minimum number of grid points),
2316 and {\tt pme-order} controls the
2317 interpolation order. Using fourth-order (cubic) interpolation and this
2318 spacing should give electrostatic energies accurate to about
2319 $5\cdot10^{-3}$. Since the Lennard-Jones energies are not this
2320 accurate it might even be possible to increase this spacing slightly.
2322 Pressure scaling works with PME, but be aware of the fact that
2323 anisotropic scaling can introduce artificial ordering in some systems.
2325 \subsection{\normindex{P3M-AD}}
2327 The \seeindex{Particle-Particle Particle-Mesh}{P3M} methods of
2328 Hockney \& Eastwood can also be applied in {\gromacs} for the
2329 treatment of long range electrostatic interactions~\cite{Hockney81}.
2330 Although the P3M method was the first efficient long-range electrostatics
2331 method for molecular simulation, the smooth PME (SPME) method has largely
2332 replaced P3M as the method of choice in atomistic simulations. One performance
2333 disadvantage of the original P3M method was that it required 3 3D-FFT
2334 back transforms to obtain the forces on the particles. But this is not
2335 required for P3M and the forces can be derived through analytical differentiation
2336 of the potential, as done in PME. The resulting method is termed P3M-AD.
2337 The only remaining difference between P3M-AD and PME is the optimization
2338 of the lattice Green influence function for error minimization that P3M uses.
2339 However, in 2012 it has been shown that the SPME influence function can be
2340 modified to obtain P3M~\cite{Ballenegger2012}.
2341 This means that the advantage of error minimization in P3M-AD can be used
2342 at the same computational cost and with the same code as PME,
2343 just by adding a few lines to modify the influence function.
2344 However, at optimal parameter setting the effect of error minimization
2345 in P3M-AD is less than 10\%. P3M-AD does show large accuracy gains with
2346 interlaced (also known as staggered) grids, but that is not supported
2347 in {\gromacs} (yet).
2349 P3M is used in {\gromacs} with exactly the same options as used with PME
2350 by selecting the electrostatics type:
2352 coulombtype = P3M-AD
2355 \subsection{Optimizing Fourier transforms}
2356 To get the best possible performance you should try to avoid large
2357 prime numbers for grid dimensions.
2358 The FFT code used in {\gromacs} is
2359 optimized for grid sizes of the form $2^a 3^b 5^c 7^d 11^e 13^f$,
2360 where $e+f$ is $0$ or $1$ and the other exponents arbitrary. (See
2361 further the documentation of the FFT algorithms at
2362 \href{http://www.fftw.org}{www.fftw.org}.
2364 It is also possible to optimize the transforms for the current problem
2365 by performing some calculations at the start of the run. This is not
2366 done by default since it takes a couple of minutes, but for large
2367 runs it will save time. Turn it on by specifying
2372 in your {\tt .mdp} file.
2374 When running in parallel, the grid must be communicated several times,
2375 thus hurting scaling performance. With PME you can improve this
2376 by increasing grid spacing while simultaneously increasing the
2377 interpolation to {\eg} sixth order.
2378 Since the interpolation is entirely local, doing so will
2379 improve the scaling in most cases.
2382 % Temporarily removed since I am not sure about the state of the testlr
2385 %It is possible to test the accuracy of your settings using the program
2386 %{\tt\normindex{testlr}} in the {\tt src/gmxlib} dir. This program computes
2387 %forces and potentials using PPPM and an Ewald implementation and gives the
2388 %absolute and RMS errors in both:
2393 %Potential 0.113 0.035
2395 %{\bf Note:} these numbers were generated using a grid spacing of
2396 %0.058 nm and $r_c$ = 1.0 nm.
2398 %You can see what the accuracy is without optimizing the
2399 %$\hat{G}(k)$ function, if you pass the {\tt -ghat} option to {\tt
2400 %testlr}. Try it if you think the {\tt mk_ghat} procedure is a waste
2402 %} % Brace matches ifthenelse test for gmxlite
2405 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2406 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2407 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2409 %\ifthenelse{\equal{\gmxlite}{1}}{}{
2410 \section{Long Range Van der Waals interactions}
2411 \subsection{Dispersion correction\index{dispersion correction}}
2412 In this section, we derive long-range corrections due to the use of a
2413 cut-off for Lennard-Jones or Buckingham interactions.
2414 We assume that the cut-off is
2415 so long that the repulsion term can safely be neglected, and therefore
2416 only the dispersion term is taken into account. Due to the nature of
2417 the dispersion interaction (we are truncating a potential proportional
2418 to $-r^{-6}$), energy and pressure corrections are both negative. While
2419 the energy correction is usually small, it may be important for free
2420 energy calculations where differences between two different Hamiltonians
2421 are considered. In contrast, the pressure correction is very large and
2422 can not be neglected under any circumstances where a correct pressure is
2423 required, especially for any NPT simulations. Although it is, in
2424 principle, possible to parameterize a force field such that the pressure
2425 is close to the desired experimental value without correction, such a
2426 method makes the parameterization dependent on the cut-off and is therefore
2429 \subsubsection{Energy}
2431 The long-range contribution of the dispersion interaction to the
2432 virial can be derived analytically, if we assume a homogeneous
2433 system beyond the cut-off distance $r_c$. The dispersion energy
2434 between two particles is written as:
2436 V(\rij) ~=~- C_6\,\rij^{-6}
2438 and the corresponding force is:
2440 \Fvij ~=~- 6\,C_6\,\rij^{-8}\rvij
2442 In a periodic system it is not easy to calculate the full potentials,
2443 so usually a cut-off is applied, which can be abrupt or smooth.
2444 We will call the potential and force with cut-off $V_c$ and $\ve{F}_c$.
2445 The long-range contribution to the dispersion energy
2446 in a system with $N$ particles and particle density $\rho$ = $N/V$ is:
2448 \label{eqn:enercorr}
2449 V_{lr} ~=~ \half N \rho\int_0^{\infty} 4\pi r^2 g(r) \left( V(r) -V_c(r) \right) {\dr}
2451 We will integrate this for the shift function, which is the most general
2452 form of van der Waals interaction available in {\gromacs}.
2453 The shift function has a constant difference $S$ from 0 to $r_1$
2454 and is 0 beyond the cut-off distance $r_c$.
2455 We can integrate \eqnref{enercorr}, assuming that the density in the sphere
2456 within $r_1$ is equal to the global density and
2457 the radial distribution function $g(r)$ is 1 beyond $r_1$:
2460 V_{lr} &=& \half N \left(
2461 \rho\int_0^{r_1} 4\pi r^2 g(r) \, C_6 \,S\,{\dr}
2462 + \rho\int_{r_1}^{r_c} 4\pi r^2 \left( V(r) -V_c(r) \right) {\dr}
2463 + \rho\int_{r_c}^{\infty} 4\pi r^2 V(r) \, {\dr}
2465 & = & \half N \left(\left(\frac{4}{3}\pi \rho r_1^{3} - 1\right) C_6 \,S
2466 + \rho\int_{r_1}^{r_c} 4\pi r^2 \left( V(r) -V_c(r) \right) {\dr}
2467 -\frac{4}{3} \pi N \rho\, C_6\,r_c^{-3}
2470 where the term $-1$ corrects for the self-interaction.
2471 For a plain cut-off we only need to assume that $g(r)$ is 1 beyond $r_c$
2472 and the correction reduces to~\cite{Allen87}:
2474 V_{lr} & = & -\frac{2}{3} \pi N \rho\, C_6\,r_c^{-3}
2476 If we consider, for example, a box of pure water, simulated with a cut-off
2477 of 0.9 nm and a density of 1 g cm$^{-3}$ this correction is
2478 $-0.75$ kJ mol$^{-1}$ per molecule.
2480 For a homogeneous mixture we need to define
2481 an {\em average dispersion constant}:
2484 \avcsix = \frac{2}{N(N-1)}\sum_i^N\sum_{j>i}^N C_6(i,j)\\
2486 In {\gromacs}, excluded pairs of atoms do not contribute to the average.
2488 In the case of inhomogeneous simulation systems, {\eg} a system with a
2489 lipid interface, the energy correction can be applied if
2490 $\avcsix$ for both components is comparable.
2492 \subsubsection{Virial and pressure}
2493 The scalar virial of the system due to the dispersion interaction between
2494 two particles $i$ and $j$ is given by:
2496 \Xi~=~-\half \rvij \cdot \Fvij ~=~ 3\,C_6\,\rij^{-6}
2498 The pressure is given by:
2500 P~=~\frac{2}{3\,V}\left(E_{kin} - \Xi\right)
2502 The long-range correction to the virial is given by:
2504 \Xi_{lr} ~=~ \half N \rho \int_0^{\infty} 4\pi r^2 g(r) (\Xi -\Xi_c) \,\dr
2506 We can again integrate the long-range contribution to the
2507 virial assuming $g(r)$ is 1 beyond $r_1$:
2509 \Xi_{lr}&=& \half N \rho \left(
2510 \int_{r_1}^{r_c} 4 \pi r^2 (\Xi -\Xi_c) \,\dr
2511 + \int_{r_c}^{\infty} 4 \pi r^2 3\,C_6\,\rij^{-6}\, \dr
2513 &=& \half N \rho \left(
2514 \int_{r_1}^{r_c} 4 \pi r^2 (\Xi -\Xi_c) \, \dr
2515 + 4 \pi C_6 \, r_c^{-3} \right)
2517 For a plain cut-off the correction to the pressure is~\cite{Allen87}:
2519 P_{lr}~=~-\frac{4}{3} \pi C_6\, \rho^2 r_c^{-3}
2521 Using the same example of a water box, the correction to the virial is
2522 0.75 kJ mol$^{-1}$ per molecule,
2523 the corresponding correction to the pressure for
2524 SPC water is approximately $-280$ bar.
2526 For homogeneous mixtures, we can again use the average dispersion constant
2527 $\avcsix$ (\eqnref{avcsix}):
2529 P_{lr}~=~-\frac{4}{3} \pi \avcsix \rho^2 r_c^{-3}
2532 For inhomogeneous systems, \eqnref{pcorr} can be applied under the same
2533 restriction as holds for the energy (see \secref{ecorr}).
2535 \subsection{Lennard-Jones PME\index{LJ-PME}}
2537 In order to treat systems, using Lennard-Jones potentials, that are
2538 non-homogeneous outside of the cut-off distance, we can instead use
2539 the Particle-mesh Ewald method as discussed for electrostatics above.
2540 In this case the modified Ewald equations become
2542 V &=& V_{dir} + V_{rec} + V_{0} \\[0.5ex]
2543 V_{dir} &=& -\frac{1}{2} \sum_{i,j}^{N}
2544 \sum_{n_x}\sum_{n_y}
2545 \sum_{n_{z}*} \frac{C_{ij}^{(6)}g(\beta {r}_{ij,{\bf n}})}{{r_{ij,{\bf n}}}^6} \\[0.5ex]
2546 V_{rec} &=& \frac{{\pi}^{\frac{3}{2}} \beta^{3}}{2V} \sum_{m_x}\sum_{m_y}\sum_{m_{z}*}
2547 f(\pi |{\mathbf m}|/\beta) \times \sum_{i,j}^{N} C_{ij}^{(6)} {\mathrm{exp}}\left[-2\pi i {\bf m}\cdot({\bf r_i}-{\bf r_j})\right] \\[0.5ex]
2548 V_{0} &=& -\frac{\beta^{6}}{12}\sum_{i}^{N} C_{ii}^{(6)},
2551 where ${\bf m}=(m_x,m_y,m_z)$, $\beta$ is the parameter determining the weight between
2552 direct and reciprocal space, and ${C_{ij}}^{(6)}$ is the combined dispersion
2553 parameter for particle $i$ and $j$. The star indicates that terms
2554 with $i = j$ should be omitted when $((n_x,n_y,n_z)=(0,0,0))$, and
2555 ${\bf r}_{ij,{\bf n}}$ is the real distance between the particles.
2556 Following the derivation by Essmann~\cite{Essmann95}, the functions $f$ and $g$ introduced above are defined as
2558 f(x)&=&1/3\left[(1-2x^2){\mathrm{exp}}(-x^2) + 2{x^3}\sqrt{\pi}\,{\mathrm{erfc}}(x) \right] \\
2559 g(x)&=&{\mathrm{exp}}(-x^2)(1+x^2+\frac{x^4}{2}).
2562 The above methodology works fine as long as the dispersion parameters can be factorized in the same
2563 way as the charges for electrostatics
2565 C_{ij}^{(6)} = \left({C_{ii}^{(6)} \, C_{jj}^{(6)}}\right)^{1/2}
2567 For Lorentz-Berthelot combination rules, the reciprocal part of this sum has to be calculated
2568 seven times due to the splitting of the dispersion parameter according to
2570 C_{ij}^{(6)}=(\sigma_i+\sigma_j)^6=\sum_{n=0}^{6} P_{n}\sigma_{i}^{n}\sigma_{j}^{(6-n)},
2572 for $P_{n}$ the Pascal triangle coefficients. This introduces a
2573 non-negligible cost to the reciprocal part, requiring seven separate
2574 FFTs, and therefore this have been the limiting factor in previous
2575 attempts to implement LJ-PME. A solution to this problem is to
2576 approximate the interaction parameters in reciprocal space using
2577 geometrical combination rules. This will preserve a well-defined
2578 Hamiltonian and significantly increase the performance of the
2582 \centerline{\includegraphics[width=15cm]{plots/ljpmedifference}}
2583 \caption {Dispersion potential between phosphorous and oxygen, The total, real and
2584 reciprocal parts of the total dispersion are shown. The reciprocal parts are calculated
2585 using either LB or geometric rules. The difference introduced by the use of the
2586 geometric approximation in the reciprocal part is small compared to the total
2587 interaction energy.}
2588 \label{fig:ljpmedifference}
2591 This approximation does introduce some errors, but since the
2592 difference is located in the interactions calculated in reciprocal
2593 space, the effect will be very small compared to the total interaction
2594 energy (see \figref{ljpmedifference}). The relative error in
2595 the total dispersion energy will stay below 0.5\% in a lipid bilayer
2596 simulation, when using a real space cut-off of 1.0 nm. A more
2597 thorough discussion of this can be found in \cite{Wennberg13}.
2599 \subsubsection{Using LJ-PME}
2600 To use Particle-mesh Ewald summation for Lennard-Jones interactions in {\gromacs}, specify the
2601 following lines in your {\tt .mdp} file:
2607 fourierspacing = 0.12
2609 ewald-rtol-lj = 0.001
2610 lj-pme-comb-rule = geometric
2613 The {\tt fourierspacing} and {\tt pme-order} are the same parameters
2614 as is used for electrostatic PME, and {\tt ewald-rtol-lj} controls
2615 splitting between real and reciprocal space in the same way as
2616 {\tt ewald-rtol}. In addition to this, the combination rule to be used
2617 in reciprocal space is determined by {\tt lj-pme-comb-rule}. If the
2618 current force field uses Lorentz-Berthelot combination rules, it is
2619 possible to set {\tt lj-pme-comb-rule = geometric} in order to gain a
2620 significant increase in performance for a small loss in accuracy. The
2621 details of this approximation can be found in the section above. The
2622 implementation of LJ-PME is currently unsupported together with the
2623 Verlet cut-off scheme and/or free energy calculations. These features
2624 will be added in upcoming releases
2625 %} % Brace matches ifthenelse test for gmxlite
2627 \section{Force field\index{force field}}
2629 A force field is built up from two distinct components:
2631 \item The set of equations (called the {\em
2632 \index{potential function}s}) used to generate the potential
2633 energies and their derivatives, the forces. These are described in
2634 detail in the previous chapter.
2635 \item The parameters used in this set of equations. These are not
2636 given in this manual, but in the data files corresponding to your
2637 {\gromacs} distribution.
2639 Within one set of equations various sets of parameters can be
2640 used. Care must be taken that the combination of equations and
2641 parameters form a consistent set. It is in general dangerous to make
2642 {\em ad hoc} changes in a subset of parameters, because the various
2643 contributions to the total force are usually interdependent. This
2644 means in principle that every change should be documented, verified by
2645 comparison to experimental data and published in a peer-reviewed
2646 journal before it can be used.
2648 {\gromacs} {\gmxver} includes several force fields, and additional
2649 ones are available on the website. If you do not know which one to
2650 select we recommend \gromosv{96} for united-atom setups and OPLS-AA/L for
2651 all-atom parameters. That said, we describe the available options in
2654 \subsection{GROMOS87\index{GROMOS87 force field}}
2655 The \gromosv{87} suite of programs and corresponding force
2656 field~\cite{biomos} formed the basis for the development of {\gromacs}
2657 in the early 1990s. The original GROMOS87 force field is not
2658 available in {\gromacs}. In previous versions ($<$ 3.3.2) there used
2659 to be the so-called ``{\gromacs} force field,'' which was based on
2660 \gromosv{87}~\cite{biomos}\index{GROMOS87}, with a small modification
2661 concerning the interaction between water oxygens and carbon
2662 atoms~\cite{Buuren93b,Mark94}, as well as 10 extra atom
2663 types~\cite{Jorgensen83,Buuren93a,Buuren93b,Mark94,Liu95}. Whenever
2664 using this force field, please cite the above references, and {\bf do not}
2665 call it the ``{\gromacs} force field,'' instead name it
2666 \gromosv{87}~\cite{biomos} with corrections as detailed
2667 in~\cite{Buuren93b,Mark94}. As noted by {\tt pdb2gmx}, this force field is
2668 ``deprecated,'' indicating that newer, perhaps more reliable, versions of this
2669 parameter set are available. For backwards compatibility, it is maintained
2670 in the current release. Should you have a justifiable reason to use this
2671 force field, all necessary files are provided in the {\tt gmx.ff} sub-directory
2672 of the {\gromacs} library. See also the note in~\ssecref{atomtype}.
2674 \subsubsection{All-hydrogen force field}
2675 The \gromosv{87}-based all-hydrogen force field is almost identical to the
2676 normal \gromosv{87} force field, since the extra hydrogens have no
2677 Lennard-Jones interaction and zero charge. The only differences are in
2678 the bond angle and improper dihedral angle terms. This force field is
2679 only useful when you need the exact hydrogen positions, for instance
2680 for distance restraints derived from NMR measurements. When citing
2681 this force field please read the previous paragraph.
2683 \subsection{\gromosv{96}\index{GROMOS96 force field}}
2684 {\gromacs} supports the \gromosv{96} force fields~\cite{gromos96}.
2685 All parameters for the 43A1, 43A2 (development, improved alkane
2686 dihedrals), 45A3, 53A5, and 53A6 parameter sets are included. All standard
2687 building blocks are included and topologies can be built automatically
2690 The \gromosv{96} force field is a further development of the \gromosv{87} force field.
2691 It has improvements over the \gromosv{87} force field for proteins and small molecules.
2692 {\bf Note} that the sugar parameters present in 53A6 do correspond to those published in
2693 2004\cite{Oostenbrink2004}, which are different from those present in 45A4, which
2694 is not included in {\gromacs} at this time. The 45A4 parameter set corresponds to a later
2695 revision of these parameters.
2696 The \gromosv{96} force field is not, however, recommended for use with long alkanes and
2697 lipids. The \gromosv{96} force field differs from the \gromosv{87}
2698 force field in a few respects:
2700 \item the force field parameters
2701 \item the parameters for the bonded interactions are not linked to atom types
2702 \item a fourth power bond stretching potential (\ssecref{G96bond})
2703 \item an angle potential based on the cosine of the angle (\ssecref{G96angle})
2705 There are two differences in implementation between {\gromacs} and \gromosv{96}
2706 which can lead to slightly different results when simulating the same system
2709 \item in \gromosv{96} neighbor searching for solvents is performed on the
2710 first atom of the solvent molecule. This is not implemented in {\gromacs},
2711 but the difference with searching by centers of charge groups is very small
2712 \item the virial in \gromosv{96} is molecule-based. This is not implemented in
2713 {\gromacs}, which uses atomic virials
2715 The \gromosv{96} force field was parameterized with a Lennard-Jones cut-off
2716 of 1.4 nm, so be sure to use a Lennard-Jones cut-off ({\tt rvdw}) of at least 1.4.
2717 A larger cut-off is possible because the Lennard-Jones potential and forces
2718 are almost zero beyond 1.4 nm.
2720 \subsubsection{\gromosv{96} files\swapindexquiet{GROMOS96}{files}}
2721 {\gromacs} can read and write \gromosv{96} coordinate and trajectory files.
2722 These files should have the extension {\tt .g96}.
2723 Such a file can be a \gromosv{96} initial/final
2724 configuration file, a coordinate trajectory file, or a combination of both.
2725 The file is fixed format; all floats are written as 15.9, and as such, files can get huge.
2726 {\gromacs} supports the following data blocks in the given order:
2736 POSITION/POSITIONRED (mandatory)
2737 VELOCITY/VELOCITYRED (optional)
2742 See the \gromosv{96} manual~\cite{gromos96} for a complete description
2743 of the blocks. {\bf Note} that all {\gromacs} programs can read compressed
2744 (.Z) or gzipped (.gz) files.
2746 \subsection{OPLS/AA\index{OPLS/AA force field}}
2748 \subsection{AMBER\index{AMBER force field}}
2750 As of version 4.5, {\gromacs} provides native support for the following AMBER force fields:
2753 \item AMBER94~\cite{Cornell1995}
2754 \item AMBER96~\cite{Kollman1996}
2755 \item AMBER99~\cite{Wang2000}
2756 \item AMBER99SB~\cite{Hornak2006}
2757 \item AMBER99SB-ILDN~\cite{Lindorff2010}
2758 \item AMBER03~\cite{Duan2003}
2759 \item AMBERGS~\cite{Garcia2002}
2762 \subsection{CHARMM\index{CHARMM force field}}
2764 As of version 4.5, {\gromacs} supports the CHARMM27 force field for proteins~\cite{mackerell04, mackerell98}, lipids~\cite{feller00} and nucleic acids~\cite{foloppe00}. The protein parameters (and to some extent the lipid and nucleic acid parameters) were thoroughly tested -- both by comparing potential energies between the port and the standard parameter set in the CHARMM molecular simulation package, as well by how the protein force field behaves together with {\gromacs}-specific techniques such as virtual sites (enabling long time steps) and a fast implicit solvent recently implemented~\cite{Larsson10} -- and the details and results are presented in the paper by Bjelkmar et al.~\cite{Bjelkmar10}. The nucleic acid parameters, as well as the ones for HEME, were converted and tested by Michel Cuendet.
2766 When selecting the CHARMM force field in {\tt \normindex{pdb2gmx}} the default option is to use \normindex{CMAP} (for torsional correction map). To exclude CMAP, use {\tt -nocmap}. The basic form of the CMAP term implemented in {\gromacs} is a function of the $\phi$ and $\psi$ backbone torsion angles. This term is defined in the {\tt .rtp} file by a {\tt [ cmap ]} statement at the end of each residue supporting CMAP. The following five atom names define the two torsional angles. Atoms 1-4 define $\phi$, and atoms 2-5 define $\psi$. The corresponding atom types are then matched to the correct CMAP type in the {\tt cmap.itp} file that contains the correction maps.
2768 \subsection{Coarse-grained force-fields}
2769 \index{force-field, coarse-grained}
2770 \label{sec:cg-forcefields}
2771 Coarse-graining is a systematic way of reducing the number of degrees of freedom representing a system of interest. To achieve this, typically whole groups of atoms are represented by single beads and the coarse-grained force fields describes their effective interactions. Depending on the choice of parameterization, the functional form of such an interaction can be complicated and often tabulated potentials are used.
2773 Coarse-grained models are designed to reproduce certain properties of a reference system. This can be either a full atomistic model or even experimental data. Depending on the properties to reproduce there are different methods to derive such force fields. An incomplete list of methods is given below:
2775 \item Conserving free energies
2777 \item Simplex method
2778 \item MARTINI force-field (see next section)
2780 \item Conserving distributions (like the radial distribution function), so-called structure-based coarse-graining
2782 \item (iterative) Boltzmann inversion
2783 \item Inverse Monte Carlo
2785 \item Conversing forces
2787 \item Force matching
2791 Note that coarse-grained potentials are state dependent (e.g. temperature, density,...) and should be re-parametrized depending on the system of interest and the simulation conditions. This can for example be done using the \normindex{Versatile Object-oriented Toolkit for Coarse-Graining Applications (VOTCA)}~\cite{ruehle2009}. The package was designed to assists in systematic coarse-graining, provides implementations for most of the algorithms mentioned above and has a well tested interface to {\gromacs}. It is available as open source and further information can be found at \href{http://www.votca.org}{www.votca.org}.
2793 \subsection{MARTINI\index{Martini force field}}
2795 The MARTINI force field is a coarse-grain parameter set that allows for the construction
2796 of many systems, including proteins and membranes.
2798 \subsection{PLUM\index{PLUM force field}}
2800 The \normindex{PLUM force field}~\cite{bereau12} is an example of a solvent-free
2801 protein-membrane model for which the membrane was derived from structure-based
2802 coarse-graining~\cite{wang_jpcb10}. A {\gromacs} implementation can be found at
2803 \href{http://code.google.com/p/plumx/}{code.google.com/p/plumx}.
2805 % LocalWords: dihedrals centro ij dV dr LJ lj rcl jj Bertelot OPLS bh bham rf
2806 % LocalWords: coul defunits grompp crf vcrf fcrf Tironi Debye kgrf cgrf krf dx
2807 % LocalWords: PPPM der Waals erfc lr elstat chirality bstretch bondpot kT kJ
2808 % LocalWords: anharmonic morse mol betaij expminx SPC timestep fs FENE ijk kj
2809 % LocalWords: anglepot CHARMm UB ik rr substituents ijkl Ryckaert Bellemans rb
2810 % LocalWords: alkanes pdb gmx IUPAC IUB jkl cis rbdih crb kcal cubicspline xvg
2811 % LocalWords: topfile mdrun posres ar dihr lcllll NMR nmr lcllllll NOEs lclll
2812 % LocalWords: rav preprocessor ccccccccc ai aj fac disre mdp multi topol tpr
2813 % LocalWords: fc ravdisre nstdisreout dipolar lll ccc orientational MSD const
2814 % LocalWords: orire fitgrp nstorireout Drude intra Noskov et al fecalc coulrf
2815 % LocalWords: polarizabilities parameterized sigeps Ek sc softcore GROMOS NBF
2816 % LocalWords: hydrogens alkane vdwtype coulombtype mdpopt rlist rcoulomb rvdw
2817 % LocalWords: nstlist virial funcparm VdW jk jl fvsite fd vsites lcr vsitetop
2818 % LocalWords: vsite lclllll lcl parameterize parameterization enercorr avcsix
2819 % LocalWords: pcorr ecorr totalcoulomb dir fourierspacing ewald rtol Darden gz
2820 % LocalWords: FFT parallelization MPMD mpmd pme fft hoc Gromos gromos oxygens
2821 % LocalWords: virials POSITIONRED VELOCITYRED gzipped Charmm Larsson Bjelkmar
2822 % LocalWords: Cuendet CMAP nocmap dihedral Lennard covalent Verlet
2823 % LocalWords: Berthelot nm flexwat ferguson itp harmonicangle versa
2824 % LocalWords: harmonicbond atomtypes dihedraltypes equilibrated fdn
2825 % LocalWords: distancerestraint LINCS Coulombic ja jb il SPME ILDN
2826 % LocalWords: Hamiltonians atomtype AMBERGS rtp cmap graining VOTCA