2 % This file is part of the GROMACS molecular simulation package.
4 % Copyright (c) 2013,2014, 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}
118 \label{eqn:lorentzberthelot}
121 finally an geometric average for both parameters can be used (type 3):
124 \sigma_{ij} &=& \left({\sigma_{ii} \, \sigma_{jj}}\right)^{1/2} \\
125 \epsilon_{ij} &=& \left({\epsilon_{ii} \, \epsilon_{jj}}\right)^{1/2}
128 This last rule is used by the OPLS force field.
131 %\ifthenelse{\equal{\gmxlite}{1}}{}{
132 \subsection{\normindex{Buckingham potential}}
134 potential has a more flexible and realistic repulsion term
135 than the Lennard-Jones interaction, but is also more expensive to
136 compute. The potential form is:
138 V_{bh}(\rij) = A_{ij} \exp(-B_{ij} \rij) -
139 \frac{C_{ij}}{\rij^6}
142 \centerline{\includegraphics[width=8cm]{plots/f-bham}}
143 \caption {The Buckingham interaction.}
147 See also \figref{bham}. The force derived from this is:
149 \ve{F}_i(\rij) = \left[ A_{ij}B_{ij}\exp(-B_{ij} \rij) -
150 6\frac{C_{ij}}{\rij^7} \right] \rnorm
153 %} % Brace matches ifthenelse test for gmxlite
155 \subsection{Coulomb interaction}
157 \newcommand{\epsr}{\varepsilon_r}
158 \newcommand{\epsrf}{\varepsilon_{rf}}
159 The \normindex{Coulomb} interaction between two charge particles is given by:
161 V_c(\rij) = f \frac{q_i q_j}{\epsr \rij}
164 See also \figref{coul}, where $f = \frac{1}{4\pi \varepsilon_0} =
165 138.935\,485$ (see \chref{defunits})
168 \centerline{\includegraphics[width=8cm]{plots/vcrf}}
169 \caption[The Coulomb interaction with and without reaction field.]{The
170 Coulomb interaction (for particles with equal signed charge) with and
171 without reaction field. In the latter case $\epsr$ was 1, $\epsrf$ was 78,
172 and $r_c$ was 0.9 nm.
173 The dot-dashed line is the same as the dashed line, except for a constant.}
177 The force derived from this potential is:
179 \ve{F}_i(\rvij) = f \frac{q_i q_j}{\epsr\rij^2}\rnorm
182 A plain Coulomb interaction should only be used without cut-off or when all pairs fall within the cut-off, since there is an abrupt, large change in the force at the cut-off. In case you do want to use a cut-off, the potential can be shifted by a constant to make the potential the integral of the force. With the group cut-off scheme, this shift is only applied to non-excluded pairs. With the Verlet cut-off scheme, the shift is also applied to excluded pairs and self interactions, which makes the potential equivalent to a reaction field with $\epsrf=1$ (see below).
184 In {\gromacs} the relative \swapindex{dielectric}{constant}
186 may be set in the in the input for {\tt grompp}.
188 %\ifthenelse{\equal{\gmxlite}{1}}{}{
189 \subsection{Coulomb interaction with \normindex{reaction field}}
191 The Coulomb interaction can be modified for homogeneous systems by
192 assuming a constant dielectric environment beyond the cut-off $r_c$
193 with a dielectric constant of {$\epsrf$}. The interaction then reads:
196 f \frac{q_i q_j}{\epsr\rij}\left[1+\frac{\epsrf-\epsr}{2\epsrf+\epsr}
197 \,\frac{\rij^3}{r_c^3}\right]
198 - f\frac{q_i q_j}{\epsr r_c}\,\frac{3\epsrf}{2\epsrf+\epsr}
201 in which the constant expression on the right makes the potential
202 zero at the cut-off $r_c$. For charged cut-off spheres this corresponds
203 to neutralization with a homogeneous background charge.
204 We can rewrite \eqnref{vcrf} for simplicity as
206 V_{crf} ~=~ f \frac{q_i q_j}{\epsr}\left[\frac{1}{\rij} + k_{rf}~ \rij^2 -c_{rf}\right]
210 k_{rf} &=& \frac{1}{r_c^3}\,\frac{\epsrf-\epsr}{(2\epsrf+\epsr)} \label{eqn:krf}\\
211 c_{rf} &=& \frac{1}{r_c}+k_{rf}\,r_c^2 ~=~ \frac{1}{r_c}\,\frac{3\epsrf}{(2\epsrf+\epsr)}
214 For large $\epsrf$ the $k_{rf}$ goes to $r_c^{-3}/2$,
215 while for $\epsrf$ = $\epsr$ the correction vanishes.
217 the modified interaction is plotted, and it is clear that the derivative
218 with respect to $\rij$ (= -force) goes to zero at the cut-off distance.
219 The force derived from this potential reads:
221 \ve{F}_i(\rvij) = f \frac{q_i q_j}{\epsr}\left[\frac{1}{\rij^2} - 2 k_{rf}\rij\right]\rnorm
224 The reaction-field correction should also be applied to all excluded
225 atoms pairs, including self pairs, in which case the normal Coulomb
226 term in \eqnsref{vcrf}{fcrf} is absent.
228 Tironi {\etal} have introduced a generalized reaction field in which
229 the dielectric continuum beyond the cut-off $r_c$ also has an ionic strength
230 $I$~\cite{Tironi95}. In this case we can rewrite the constants $k_{rf}$ and
231 $c_{rf}$ using the inverse Debye screening length $\kappa$:
234 \frac{2 I \,F^2}{\varepsilon_0 \epsrf RT}
235 = \frac{F^2}{\varepsilon_0 \epsrf RT}\sum_{i=1}^{K} c_i z_i^2 \\
236 k_{rf} &=& \frac{1}{r_c^3}\,
237 \frac{(\epsrf-\epsr)(1 + \kappa r_c) + \half\epsrf(\kappa r_c)^2}
238 {(2\epsrf + \epsr)(1 + \kappa r_c) + \epsrf(\kappa r_c)^2}
240 c_{rf} &=& \frac{1}{r_c}\,
241 \frac{3\epsrf(1 + \kappa r_c + \half(\kappa r_c)^2)}
242 {(2\epsrf+\epsr)(1 + \kappa r_c) + \epsrf(\kappa r_c)^2}
245 where $F$ is Faraday's constant, $R$ is the ideal gas constant, $T$
246 the absolute temperature, $c_i$ the molar concentration for species
247 $i$ and $z_i$ the charge number of species $i$ where we have $K$
248 different species. In the limit of zero ionic strength ($\kappa=0$)
249 \eqnsref{kgrf}{cgrf} reduce to the simple forms of \eqnsref{krf}{crf}
252 \subsection{Modified non-bonded interactions}
253 \label{sec:mod_nb_int}
254 In {\gromacs}, the non-bonded potentials can be
255 modified by a shift function. The purpose of this is to replace the
256 truncated forces by forces that are continuous and have continuous
257 derivatives at the \normindex{cut-off} radius. With such forces the
258 timestep integration produces much smaller errors and there are no
259 such complications as creating charges from dipoles by the truncation
260 procedure. In fact, by using shifted forces there is no need for
261 charge groups in the construction of neighbor lists. However, the
262 shift function produces a considerable modification of the Coulomb
263 potential. Unless the ``missing'' long-range potential is properly
264 calculated and added (through the use of PPPM, Ewald, or PME), the
265 effect of such modifications must be carefully evaluated. The
266 modification of the Lennard-Jones dispersion and repulsion is only
267 minor, but it does remove the noise caused by cut-off effects.
269 There is {\em no} fundamental difference between a switch function
270 (which multiplies the potential with a function) and a shift function
271 (which adds a function to the force or potential)~\cite{Spoel2006a}. The switch
272 function is a special case of the shift function, which we apply to
273 the {\em force function} $F(r)$, related to the electrostatic or
274 van der Waals force acting on particle $i$ by particle $j$ as:
276 \ve{F}_i = c F(r_{ij}) \frac{\rvij}{r_{ij}}
278 For pure Coulomb or Lennard-Jones interactions
279 $F(r)=F_\alpha(r)=r^{-(\alpha+1)}$.
280 The shifted force $F_s(r)$ can generally be written as:
284 F_s(r)~=&~F_\alpha(r) & r < r_1 \\
286 F_s(r)~=&~F_\alpha(r)+S(r) & r_1 \le r < r_c \\
287 F_s(r)~=&~0 & r_c \le r
290 When $r_1=0$ this is a traditional shift function, otherwise it acts as a
291 switch function. The corresponding shifted coulomb potential then reads:
293 V_s(r_{ij}) = f \Phi_s (r_{ij}) q_i q_j
295 where $\Phi(r)$ is the potential function
297 \Phi_s(r) = \int^{\infty}_r~F_s(x)\, dx
300 The {\gromacs} shift function should be smooth at the boundaries, therefore
301 the following boundary conditions are imposed on the shift function:
306 S(r_c) &=&-F_\alpha(r_c) \\
307 S'(r_c) &=&-F_\alpha'(r_c)
310 A 3$^{rd}$ degree polynomial of the form
312 S(r) = A(r-r_1)^2 + B(r-r_1)^3
314 fulfills these requirements. The constants A and B are given by the
315 boundary condition at $r_c$:
319 A &~=~& -\displaystyle
320 \frac{(\alpha+4)r_c~-~(\alpha+1)r_1} {r_c^{\alpha+2}~(r_c-r_1)^2} \\
321 B &~=~& \displaystyle
322 \frac{(\alpha+3)r_c~-~(\alpha+1)r_1}{r_c^{\alpha+2}~(r_c-r_1)^3}
325 Thus the total force function is:
327 F_s(r) = \frac{\alpha}{r^{\alpha+1}} + A(r-r_1)^2 + B(r-r_1)^3
329 and the potential function reads:
331 \Phi(r) = \frac{1}{r^\alpha} - \frac{A}{3} (r-r_1)^3 - \frac{B}{4} (r-r_1)^4 - C
335 C = \frac{1}{r_c^\alpha} - \frac{A}{3} (r_c-r_1)^3 - \frac{B}{4} (r_c-r_1)^4
338 When $r_1$ = 0, the modified Coulomb force function is
340 F_s(r) = \frac{1}{r^2} - \frac{5 r^2}{r_c^4} + \frac{4 r^3}{r_c^5}
342 which is identical to the {\em \index{parabolic force}}
343 function recommended to be used as a short-range function in
344 conjunction with a \swapindex{Poisson}{solver}
345 for the long-range part~\cite{Berendsen93a}.
346 The modified Coulomb potential function is:
348 \Phi(r) = \frac{1}{r} - \frac{5}{3r_c} + \frac{5r^3}{3r_c^4} - \frac{r^4}{r_c^5}
350 See also \figref{shift}.
353 \centerline{\includegraphics[width=10cm]{plots/shiftf}}
354 \caption[The Coulomb Force, Shifted Force and Shift Function
355 $S(r)$,.]{The Coulomb Force, Shifted Force and Shift Function $S(r)$,
356 using r$_1$ = 2 and r$_c$ = 4.}
360 \subsection{Modified short-range interactions with Ewald summation}
361 When Ewald summation\index{Ewald sum} or \seeindex{particle-mesh
362 Ewald}{PME}\index{Ewald, particle-mesh} is used to calculate the
363 long-range interactions, the
364 short-range Coulomb potential must also be modified, similar to the
365 switch function above. In this case the short range potential is given
368 V(r) = f \frac{\mbox{erfc}(\beta r_{ij})}{r_{ij}} q_i q_j,
370 where $\beta$ is a parameter that determines the relative weight
371 between the direct space sum and the reciprocal space sum and erfc$(x)$ is
372 the complementary error function. For further
373 details on long-range electrostatics, see \secref{lr_elstat}.
374 %} % Brace matches ifthenelse test for gmxlite
377 \section{Bonded interactions}
378 Bonded interactions are based on a fixed list of atoms. They are not
379 exclusively pair interactions, but include 3- and 4-body interactions
380 as well. There are {\em bond stretching} (2-body), {\em bond angle}
381 (3-body), and {\em dihedral angle} (4-body) interactions. A special
382 type of dihedral interaction (called {\em improper dihedral}) is used
383 to force atoms to remain in a plane or to prevent transition to a
384 configuration of opposite chirality (a mirror image).
386 \subsection{Bond stretching}
388 \subsubsection{Harmonic potential}
389 \label{subsec:harmonicbond}
390 The \swapindex{bond}{stretching} between two covalently bonded atoms
391 $i$ and $j$ is represented by a harmonic potential:
394 \centerline{\raisebox{2cm}{\includegraphics[width=5cm]{plots/bstretch}}\includegraphics[width=7cm]{plots/f-bond}}
395 \caption[Bond stretching.]{Principle of bond stretching (left), and the bond
396 stretching potential (right).}
397 \label{fig:bstretch1}
401 V_b~(\rij) = \half k^b_{ij}(\rij-b_{ij})^2
403 See also \figref{bstretch1}, with the force given by:
405 \ve{F}_i(\rvij) = k^b_{ij}(\rij-b_{ij}) \rnorm
408 %\ifthenelse{\equal{\gmxlite}{1}}{}{
409 \subsubsection{Fourth power potential}
410 \label{subsec:G96bond}
411 In the \gromosv{96} force field~\cite{gromos96}, the covalent bond potential
412 is, for reasons of computational efficiency, written as:
414 V_b~(\rij) = \frac{1}{4}k^b_{ij}\left(\rij^2-b_{ij}^2\right)^2
416 The corresponding force is:
418 \ve{F}_i(\rvij) = k^b_{ij}(\rij^2-b_{ij}^2)~\rvij
420 The force constants for this form of the potential are related to the usual
421 harmonic force constant $k^{b,\mathrm{harm}}$ (\secref{bondpot}) as
423 2 k^b b_{ij}^2 = k^{b,\mathrm{harm}}
425 The force constants are mostly derived from the harmonic ones used in
426 \gromosv{87}~\cite{biomos}. Although this form is computationally more
428 (because no square root has to be evaluated), it is conceptually more
429 complex. One particular disadvantage is that since the form is not harmonic,
430 the average energy of a single bond is not equal to $\half kT$ as it is for
431 the normal harmonic potential.
433 \subsection{\normindex{Morse potential} bond stretching}
434 \label{subsec:Morsebond}
435 %\author{F.P.X. Everdij}
437 For some systems that require an anharmonic bond stretching potential,
438 the Morse potential~\cite{Morse29}
439 between two atoms {\it i} and {\it j} is available
440 in {\gromacs}. This potential differs from the harmonic potential in
441 that it has an asymmetric potential well and a zero force at infinite
442 distance. The functional form is:
444 \displaystyle V_{morse} (r_{ij}) = D_{ij} [1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))]^2,
446 See also \figref{morse}, and the corresponding force is:
449 \displaystyle {\bf F}_{morse} ({\bf r}_{ij})&=&2 D_{ij} \beta_{ij} r_{ij} \exp(-\beta_{ij}(r_{ij}-b_{ij})) * \\
450 \displaystyle \: & \: &[1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))] \frac{\displaystyle {\bf r}_{ij}}{\displaystyle r_{ij}},
453 where \( \displaystyle D_{ij} \) is the depth of the well in kJ/mol,
454 \( \displaystyle \beta_{ij} \) defines the steepness of the well (in
455 nm\(^{-1} \)), and \( \displaystyle b_{ij} \) is the equilibrium
456 distance in nm. The steepness parameter \( \displaystyle \beta_{ij}
457 \) can be expressed in terms of the reduced mass of the atoms {\it i}
458 and {\it j}, the fundamental vibration frequency \( \displaystyle
459 \omega_{ij} \) and the well depth \( \displaystyle D_{ij} \):
461 \displaystyle \beta_{ij}= \omega_{ij} \sqrt{\frac{\mu_{ij}}{2 D_{ij}}}
463 and because \( \displaystyle \omega = \sqrt{k/\mu} \), one can rewrite \( \displaystyle \beta_{ij} \) in terms of the harmonic force constant \( \displaystyle k_{ij} \):
465 \displaystyle \beta_{ij}= \sqrt{\frac{k_{ij}}{2 D_{ij}}}
468 For small deviations \( \displaystyle (r_{ij}-b_{ij}) \), one can
469 approximate the \( \displaystyle \exp \)-term to first-order using a
472 \displaystyle \exp(-x) \approx 1-x
475 and substituting \eqnref{betaij} and \eqnref{expminx} in the functional form:
478 \displaystyle V_{morse} (r_{ij})&=&D_{ij} [1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))]^2\\
479 \displaystyle \:&=&D_{ij} [1 - (1 -\sqrt{\frac{k_{ij}}{2 D_{ij}}}(r_{ij}-b_{ij}))]^2\\
480 \displaystyle \:&=&\frac{1}{2} k_{ij} (r_{ij}-b_{ij}))^2
483 we recover the harmonic bond stretching potential.
486 \centerline{\includegraphics[width=7cm]{plots/f-morse}}
487 \caption{The Morse potential well, with bond length 0.15 nm.}
491 \subsection{Cubic bond stretching potential}
492 \label{subsec:cubicbond}
493 Another anharmonic bond stretching potential that is slightly simpler
494 than the Morse potential adds a cubic term in the distance to the
495 simple harmonic form:
497 V_b~(\rij) = k^b_{ij}(\rij-b_{ij})^2 + k^b_{ij}k^{cub}_{ij}(\rij-b_{ij})^3
499 A flexible \normindex{water} model (based on
500 the SPC water model~\cite{Berendsen81}) including
501 a cubic bond stretching potential for the O-H bond
502 was developed by Ferguson~\cite{Ferguson95}. This model was found
503 to yield a reasonable infrared spectrum. The Ferguson water model is
504 available in the {\gromacs} library ({\tt flexwat-ferguson.itp}).
505 It should be noted that the potential is asymmetric: overstretching leads to
506 infinitely low energies. The \swapindex{integration}{timestep} is therefore
509 The force corresponding to this potential is:
511 \ve{F}_i(\rvij) = 2k^b_{ij}(\rij-b_{ij})~\rnorm + 3k^b_{ij}k^{cub}_{ij}(\rij-b_{ij})^2~\rnorm
514 \subsection{FENE bond stretching potential\index{FENE potential}}
515 \label{subsec:FENEbond}
516 In coarse-grained polymer simulations the beads are often connected
517 by a FENE (finitely extensible nonlinear elastic) potential~\cite{Warner72}:
519 V_{\mbox{\small FENE}}(\rij) =
520 -\half k^b_{ij} b^2_{ij} \log\left(1 - \frac{\rij^2}{b^2_{ij}}\right)
522 The potential looks complicated, but the expression for the force is simpler:
524 F_{\mbox{\small FENE}}(\rvij) =
525 -k^b_{ij} \left(1 - \frac{\rij^2}{b^2_{ij}}\right)^{-1} \rvij
527 At short distances the potential asymptotically goes to a harmonic
528 potential with force constant $k^b$, while it diverges at distance $b$.
529 %} % Brace matches ifthenelse test for gmxlite
531 \subsection{Harmonic angle potential}
532 \label{subsec:harmonicangle}
533 \newcommand{\tijk}{\theta_{ijk}}
534 The bond-\swapindex{angle}{vibration} between a triplet of atoms $i$ - $j$ - $k$
535 is also represented by a harmonic potential on the angle $\tijk$
538 \centerline{\raisebox{1cm}{\includegraphics[width=5cm]{plots/angle}}\includegraphics[width=7cm]{plots/f-angle}}
539 \caption[Angle vibration.]{Principle of angle vibration (left) and the
540 bond angle potential (right).}
545 V_a(\tijk) = \half k^{\theta}_{ijk}(\tijk-\tijk^0)^2
547 As the bond-angle vibration is represented by a harmonic potential, the
548 form is the same as the bond stretching (\figref{bstretch1}).
550 The force equations are given by the chain rule:
553 \Fvi ~=~ -\displaystyle\frac{d V_a(\tijk)}{d \rvi} \\
554 \Fvk ~=~ -\displaystyle\frac{d V_a(\tijk)}{d \rvk} \\
557 ~ \mbox{ ~ where ~ } ~
558 \tijk = \arccos \frac{(\rvij \cdot \ve{r}_{kj})}{r_{ij}r_{kj}}
560 The numbering $i,j,k$ is in sequence of covalently bonded atoms. Atom
561 $j$ is in the middle; atoms $i$ and $k$ are at the ends (see \figref{angle}).
562 {\bf Note} that in the input in topology files, angles are given in degrees and
563 force constants in kJ/mol/rad$^2$.
565 %\ifthenelse{\equal{\gmxlite}{1}}{}{
566 \subsection{Cosine based angle potential}
567 \label{subsec:G96angle}
568 In the \gromosv{96} force field a simplified function is used to represent angle
571 V_a(\tijk) = \half k^{\theta}_{ijk}\left(\cos(\tijk) - \cos(\tijk^0)\right)^2
576 \cos(\tijk) = \frac{\rvij\cdot\ve{r}_{kj}}{\rij r_{kj}}
578 The corresponding force can be derived by partial differentiation with respect
579 to the atomic positions. The force constants in this function are related
580 to the force constants in the harmonic form $k^{\theta,\mathrm{harm}}$
581 (\ssecref{harmonicangle}) by:
583 k^{\theta} \sin^2(\tijk^0) = k^{\theta,\mathrm{harm}}
585 In the \gromosv{96} manual there is a much more complicated conversion formula
586 which is temperature dependent. The formulas are equivalent at 0 K
587 and the differences at 300 K are on the order of 0.1 to 0.2\%.
588 {\bf Note} that in the input in topology files, angles are given in degrees and
589 force constants in kJ/mol.
591 \subsection{Restricted bending potential}
593 The restricted bending (ReB) potential~\cite{MonicaGoga2013} prevents the bending angle $\theta$
594 from reaching the $180^{\circ}$ value. In this way, the numerical instabilities
595 due to the calculation of the torsion angle and potential are eliminated when
596 performing coarse-grained molecular dynamics simulations.
598 To systematically hinder the bending angles from reaching the $180^{\circ}$ value,
599 the bending potential \ref{eq:G96angle} is divided by a $\sin^2\theta$ factor:
602 V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta} \frac{(\cos\theta_i - \cos\theta_0)^2}{\sin^2\theta_i}.
606 Figure ~\figref{ReB} shows the comparison between the ReB potential, \ref{eq:ReB},
607 and the standard one \ref{eq:G96angle}.
610 \centerline{\includegraphics[width=10cm]{plots/fig-02}}
612 \caption{Bending angle potentials: cosine harmonic (solid black line), angle harmonic
613 (dashed black line) and restricted bending (red) with the same bending constant
614 $k_{\theta}=85$ kJ mol$^{-1}$ and equilibrium angle $\theta_0=130^{\circ}$.
615 The orange line represents the sum of a cosine harmonic ($k =50$ kJ mol$^{-1}$)
616 with a restricted bending ($k =25$ kJ mol$^{-1}$) potential, both with $\theta_0=130^{\circ}$.}
620 The wall of the ReB potential is very repulsive in the region close to $180^{\circ}$ and,
621 as a result, the bending angles are kept within a safe interval, far from instabilities.
622 The power $2$ of $\sin\theta_i$ in the denominator has been chosen to guarantee this behavior
623 and allows an elegant differentiation:
626 F_{\rm ReB}(\theta_i) = \frac{2k_{\theta}}{\sin^4\theta_i}(\cos\theta_i - \cos\theta_0) (1 - \cos\theta_i\cos\theta_0) \frac{\partial \cos\theta_i}{\partial \vec r_{k}}.
630 Due to its construction, the restricted bending potential cannot be used for equilibrium
631 $\theta_0$ values too close to $0^{\circ}$ or $180^{\circ}$ (from experience, at least $10^{\circ}$
632 difference is recommended). It is very important that, in the starting configuration,
633 all the bending angles have to be in the safe interval to avoid initial instabilities.
634 This bending potential can be used in combination with any form of torsion potential.
635 It will always prevent three consecutive particles from becoming collinear and,
636 as a result, any torsion potential will remain free of singularities.
637 It can be also added to a standard bending potential to affect the angle around $180^{\circ}$,
638 but to keep its original form around the minimum (see the orange curve in \figref{ReB}).
641 \subsection{Urey-Bradley potential}
642 \label{subsec:Urey-Bradley}
643 The \swapindex{Urey-Bradley bond-angle}{vibration} between a triplet
644 of atoms $i$ - $j$ - $k$ is represented by a harmonic potential on the
645 angle $\tijk$ and a harmonic correction term on the distance between
646 the atoms $i$ and $k$. Although this can be easily written as a simple
647 sum of two terms, it is convenient to have it as a single entry in the
648 topology file and in the output as a separate energy term. It is used mainly
649 in the CHARMm force field~\cite{BBrooks83}. The energy is given by:
652 V_a(\tijk) = \half k^{\theta}_{ijk}(\tijk-\tijk^0)^2 + \half k^{UB}_{ijk}(r_{ik}-r_{ik}^0)^2
655 The force equations can be deduced from sections~\ssecref{harmonicbond}
656 and~\ssecref{harmonicangle}.
658 \subsection{Bond-Bond cross term}
659 \label{subsec:bondbondcross}
660 The bond-bond cross term for three particles $i, j, k$ forming bonds
661 $i-j$ and $k-j$ is given by~\cite{Lawrence2003b}:
663 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)
666 where $k_{rr'}$ is the force constant, and $r_{1e}$ and $r_{2e}$ are the
667 equilibrium bond lengths of the $i-j$ and $k-j$ bonds respectively. The force
668 associated with this potential on particle $i$ is:
670 \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|}
672 The force on atom $k$ can be obtained by swapping $i$ and $k$ in the above
673 equation. Finally, the force on atom $j$ follows from the fact that the sum
674 of internal forces should be zero: $\ve{F}_j = -\ve{F}_i-\ve{F}_k$.
676 \subsection{Bond-Angle cross term}
677 \label{subsec:bondanglecross}
678 The bond-angle cross term for three particles $i, j, k$ forming bonds
679 $i-j$ and $k-j$ is given by~\cite{Lawrence2003b}:
681 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)
683 where $k_{r\theta}$ is the force constant, $r_{3e}$ is the $i-k$ distance,
684 and the other constants are the same as in Equation~\ref{crossbb}. The force
685 associated with the potential on atom $i$ is:
687 \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|} \\
688 + \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]
691 \subsection{Quartic angle potential}
692 \label{subsec:quarticangle}
693 For special purposes there is an angle potential
694 that uses a fourth order polynomial:
696 V_q(\tijk) ~=~ \sum_{n=0}^5 C_n (\tijk-\tijk^0)^n
698 %} % Brace matches ifthenelse test for gmxlite
700 %% new commands %%%%%%%%%%%%%%%%%%%%%%
701 \newcommand{\rvkj}{{\bf r}_{kj}}
702 \newcommand{\rkj}{r_{kj}}
703 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
705 \subsection{Improper dihedrals\swapindexquiet{improper}{dihedral}}
707 Improper dihedrals are meant to keep \swapindex{planar}{group}s ({\eg}
708 aromatic rings) planar, or to prevent molecules from flipping over to their
709 \normindex{mirror image}s, see \figref{imp}.
712 \centerline{\includegraphics[width=4cm]{plots/ring-imp}\hspace{1cm}
713 \includegraphics[width=3cm]{plots/subst-im}\hspace{1cm}\includegraphics[width=3cm]{plots/tetra-im}}
714 \caption[Improper dihedral angles.]{Principle of improper
715 dihedral angles. Out of plane bending for rings (left), substituents
716 of rings (middle), out of tetrahedral (right). The improper dihedral
717 angle $\xi$ is defined as the angle between planes (i,j,k) and (j,k,l)
722 \subsubsection{Improper dihedrals: harmonic type}
723 \label{subsec:harmonicimproperdihedral}
724 The simplest improper dihedral potential is a harmonic potential; it is plotted in
727 V_{id}(\xi_{ijkl}) = \half k_{\xi}(\xi_{ijkl}-\xi_0)^2
729 Since the potential is harmonic it is discontinuous,
730 but since the discontinuity is chosen at 180$^\circ$ distance from $\xi_0$
731 this will never cause problems.
732 {\bf Note} that in the input in topology files, angles are given in degrees and
733 force constants in kJ/mol/rad$^2$.
736 \centerline{\includegraphics[width=8cm]{plots/f-imps}}
737 \caption{Improper dihedral potential.}
741 \subsubsection{Improper dihedrals: periodic type}
742 \label{subsec:periodicimproperdihedral}
743 This potential is identical to the periodic proper dihedral (see below).
744 There is a separate dihedral type for this (type 4) only to be able
745 to distinguish improper from proper dihedrals in the parameter section
748 \subsection{Proper dihedrals\swapindexquiet{proper}{dihedral}}
749 For the normal \normindex{dihedral} interaction there is a choice of
750 either the {\gromos} periodic function or a function based on
751 expansion in powers of $\cos \phi$ (the so-called Ryckaert-Bellemans
752 potential). This choice has consequences for the inclusion of special
753 interactions between the first and the fourth atom of the dihedral
754 quadruple. With the periodic {\gromos} potential a special 1-4
755 LJ-interaction must be included; with the Ryckaert-Bellemans potential
756 {\em for alkanes} the \normindex{1-4 interaction}s must be excluded
757 from the non-bonded list. {\bf Note:} Ryckaert-Bellemans potentials
758 are also used in {\eg} the OPLS force field in combination with 1-4
759 interactions. You should therefore not modify topologies generated by
760 {\tt \normindex{pdb2gmx}} in this case.
762 \subsubsection{Proper dihedrals: periodic type}
763 \label{subsec:properdihedral}
764 Proper dihedral angles are defined according to the IUPAC/IUB
765 convention, where $\phi$ is the angle between the $ijk$ and the $jkl$
766 planes, with {\bf zero} corresponding to the {\em cis} configuration
767 ($i$ and $l$ on the same side). There are two dihedral function types
768 in {\gromacs} topology files. There is the standard type 1 which behaves
769 like any other bonded interactions. For certain force fields, type 9
770 is useful. Type 9 allows multiple potential functions to be applied
771 automatically to a single dihedral in the {\tt [ dihedral ]} section
772 when multiple parameters are defined for the same atomtypes
773 in the {\tt [ dihedraltypes ]} section.
776 \centerline{\raisebox{1cm}{\includegraphics[width=5cm]{plots/dih}}\includegraphics[width=7cm]{plots/f-dih}}
777 \caption[Proper dihedral angle.]{Principle of proper dihedral angle
778 (left, in {\em trans} form) and the dihedral angle potential (right).}
782 V_d(\phi_{ijkl}) = k_{\phi}(1 + \cos(n \phi - \phi_s))
785 %\ifthenelse{\equal{\gmxlite}{1}}{}{
786 \subsubsection{Proper dihedrals: Ryckaert-Bellemans function}
787 \label{subsec:RBdihedral}
788 For alkanes, the following proper dihedral potential is often used
789 (see \figref{rbdih}):
791 V_{rb}(\phi_{ijkl}) = \sum_{n=0}^5 C_n( \cos(\psi ))^n,
793 where $\psi = \phi - 180^\circ$. \\
794 {\bf Note:} A conversion from one convention to another can be achieved by
795 multiplying every coefficient \( \displaystyle C_n \)
796 by \( \displaystyle (-1)^n \).
798 An example of constants for $C$ is given in \tabref{crb}.
802 \begin{tabular}{|lr|lr|lr|}
804 $C_0$ & 9.28 & $C_2$ & -13.12 & $C_4$ & 26.24 \\
805 $C_1$ & 12.16 & $C_3$ & -3.06 & $C_5$ & -31.5 \\
809 \caption{Constants for Ryckaert-Bellemans potential (kJ mol$^{-1}$).}
814 \centerline{\includegraphics[width=8cm]{plots/f-rbs}}
815 \caption{Ryckaert-Bellemans dihedral potential.}
819 ({\bf Note:} The use of this potential implies exclusion of LJ interactions
820 between the first and the last atom of the dihedral, and $\psi$ is defined
821 according to the ``polymer convention'' ($\psi_{trans}=0$).)
823 The RB dihedral function can also be used to include Fourier dihedrals
826 V_{rb} (\phi_{ijkl}) ~=~ \frac{1}{2} \left[F_1(1+\cos(\phi)) + F_2(
827 1-\cos(2\phi)) + F_3(1+\cos(3\phi)) + F_4(1-\cos(4\phi))\right]
829 Because of the equalities \( \cos(2\phi) = 2\cos^2(\phi) - 1 \),
830 \( \cos(3\phi) = 4\cos^3(\phi) - 3\cos(\phi) \) and
831 \( \cos(4\phi) = 8\cos^4(\phi) - 8\cos^2(\phi) + 1 \)
832 one can translate the OPLS parameters to
833 Ryckaert-Bellemans parameters as follows:
837 \displaystyle C_0&=&F_2 + \frac{1}{2} (F_1 + F_3)\\
838 \displaystyle C_1&=&\frac{1}{2} (- F_1 + 3 \, F_3)\\
839 \displaystyle C_2&=& -F_2 + 4 \, F_4\\
840 \displaystyle C_3&=&-2 \, F_3\\
841 \displaystyle C_4&=&-4 \, F_4\\
842 \displaystyle C_5&=&0
845 with OPLS parameters in protein convention and RB parameters in
846 polymer convention (this yields a minus sign for the odd powers of
848 \noindent{\bf Note:} Mind the conversion from {\bf kcal mol$^{-1}$} for
849 literature OPLS and RB parameters to {\bf kJ mol$^{-1}$} in {\gromacs}.\\
850 %} % Brace matches ifthenelse test for gmxlite
852 \subsubsection{Proper dihedrals: Fourier function}
853 \label{subsec:Fourierdihedral}
854 The OPLS potential function is given as the first three
855 or four~\cite{Jorgensen2005a} cosine terms of a Fourier series.
856 In {\gromacs} the four term function is implemented:
858 V_{F} (\phi_{ijkl}) ~=~ \frac{1}{2} \left[C_1(1+\cos(\phi)) + C_2(
859 1-\cos(2\phi)) + C_3(1+\cos(3\phi)) + C_4(1+\cos(4\phi))\right],
861 %\ifthenelse{\equal{\gmxlite}{1}}{}{
862 Internally, {\gromacs}
863 uses the Ryckaert-Bellemans code
864 to compute Fourier dihedrals (see above), because this is more efficient.\\
865 \noindent{\bf Note:} Mind the conversion from {\emph kcal mol$^{-1}$} for
866 literature OPLS parameters to {\bf kJ mol$^{-1}$} in {\gromacs}.\\
868 \subsubsection{Proper dihedrals: Restricted torsion potential}
870 In a manner very similar to the restricted bending potential (see \ref{subsec:ReB}),
871 a restricted torsion/dihedral potential is introduced:
874 V_{\rm ReT}(\phi_i) = \frac{1}{2} k_{\phi} \frac{(\cos\phi_i - \cos\phi_0)^2}{\sin^2\phi_i}
878 with the advantages of being a function of $\cos\phi$ (no problems taking the derivative of $\sin\phi$)
879 and of keeping the torsion angle at only one minimum value. In this case, the factor $\sin^2\phi$ does
880 not allow the dihedral angle to move from the [$-180^{\circ}$:0] to [0:$180^{\circ}$] interval, i.e. it cannot have maxima both at $-\phi_0$ and $+\phi_0$ maxima, but only one of them.
881 For this reason, all the dihedral angles of the starting configuration should have their values in the
882 desired angles interval and the the equilibrium $\phi_0$ value should not be too close to the interval limits
883 (as for the restricted bending potential, described in \ref{subsec:ReB}, at least $10^{\circ}$ difference is recommended).
885 \subsubsection{Proper dihedrals: Combined bending-torsion potential}
887 When the four particles forming the dihedral angle become collinear (this situation will never happen in
888 atomistic simulations, but it can occur in coarse-grained simulations) the calculation of the
889 torsion angle and potential leads to numerical instabilities.
890 One way to avoid this is to use the restricted bending potential (see \ref{subsec:ReB})
891 that prevents the dihedral
892 from reaching the $180^{\circ}$ value.
894 Another way is to disregard any effects of the dihedral becoming ill-defined,
895 keeping the dihedral force and potential calculation continuous in entire angle range
896 by coupling the torsion potential (in a cosine form) with the bending potentials of the
897 adjacent bending angles in a unique expression:
900 V_{\rm CBT}(\theta_{i-1}, \theta_i, \phi_i) = k_{\phi} \sin^3\theta_{i-1} \sin^3\theta_{i} \sum_{n=0}^4 { a_n \cos^n\phi_i}.
904 This combined bending-torsion (CBT) potential has been proposed by~\cite{BulacuGiessen2005}
905 for polymer melt simulations and is extensively described in~\cite{MonicaGoga2013}.
907 This potential has two main advantages:
910 it does not only depend on the dihedral angle $\phi_i$ (between the $i-2$, $i-1$, $i$ and $i+1$ beads)
911 but also on the bending angles $\theta_{i-1}$ and $\theta_i$ defined from three adjacent beads
912 ($i-2$, $i-1$ and $i$, and $i-1$, $i$ and $i+1$, respectively).
913 The two $\sin^3\theta$ pre-factors, tentatively suggested by~\cite{ScottScheragator1966} and theoretically
914 discussed by~\cite{PaulingBond}, cancel the torsion potential and force when either of the two bending angles
915 approaches the value of $180^\circ$.
917 its dependence on $\phi_i$ is expressed through a polynomial in $\cos\phi_i$ that avoids the singularities in
918 $\phi=0^\circ$ or $180^\circ$ in calculating the torsional force.
921 These two properties make the CBT potential well-behaved for MD simulations with weak constraints
922 on the bending angles or even for steered / non-equilibrium MD in which the bending and torsion angles suffer major
924 When using the CBT potential, the bending potentials for the adjacent $\theta_{i-1}$ and $\theta_i$ may have any form.
925 It is also possible to leave out the two angle bending terms ($\theta_{i-1}$ and $\theta_{i}$) completely.
926 \figref{CBT} illustrates the difference between a torsion potential with and without the $\sin^{3}\theta$ factors
927 (blue and gray curves, respectively).
930 \centerline{\includegraphics[width=10cm]{plots/fig-04}}
931 \caption{Blue: surface plot of the combined bending-torsion potential
932 (\ref{eq:CBT} with $k = 10$ kJ mol$^{-1}$, $a_0=2.41$, $a_1=-2.95$, $a_2=0.36$, $a_3=1.33$)
933 when, for simplicity, the bending angles behave the same ($\theta_1=\theta_2=\theta$).
934 Gray: the same torsion potential without the $\sin^{3}\theta$ terms (Ryckaert-Bellemans type).
935 $\phi$ is the dihedral angle.}
939 Additionally, the derivative of $V_{CBT}$ with respect to the Cartesian variables is straightforward:
942 \frac{\partial V_{\rm CBT}(\theta_{i-1},\theta_i,\phi_i)} {\partial \vec r_{l}} = \frac{\partial V_{\rm CBT}}{\partial \theta_{i-1}} \frac{\partial \theta_{i-1}}{\partial \vec r_{l}} +
943 \frac{\partial V_{\rm CBT}}{\partial \theta_{i }} \frac{\partial \theta_{i }}{\partial \vec r_{l}} +
944 \frac{\partial V_{\rm CBT}}{\partial \phi_{i }} \frac{\partial \phi_{i }}{\partial \vec r_{l}}
948 The CBT is based on a cosine form without multiplicity, so it can only be symmetrical around $0^{\circ}$.
949 To obtain an asymmetrical dihedral angle distribution (e.g. only one maximum in [$-180^{\circ}$:$180^{\circ}$] interval),
950 a standard torsion potential such as harmonic angle or periodic cosine potentials should be used instead of a CBT potential.
951 However, these two forms have the inconveniences of the force derivation ($1/\sin\phi$) and of the alignment of beads
952 ($\theta_i$ or $\theta_{i-1} = 0^{\circ}, 180^{\circ}$).
953 Coupling such non-$\cos\phi$ potentials with $\sin^3\theta$ factors does not improve simulation stability since there are
954 cases in which $\theta$ and $\phi$ are simultaneously $180^{\circ}$. The integration at this step would be possible
955 (due to the cancelling of the torsion potential) but the next step would be singular
956 ($\theta$ is not $180^{\circ}$ and $\phi$ is very close to $180^{\circ}$).
958 %\ifthenelse{\equal{\gmxlite}{1}}{}{
959 \subsection{Tabulated bonded interaction functions\index{tabulated bonded interaction function}}
960 \label{subsec:tabulatedinteraction}
961 For full flexibility, any functional shape can be used for
962 bonds, angles and dihedrals through user-supplied tabulated functions.
963 The functional shapes are:
965 V_b(r_{ij}) &=& k \, f^b_n(r_{ij}) \\
966 V_a(\tijk) &=& k \, f^a_n(\tijk) \\
967 V_d(\phi_{ijkl}) &=& k \, f^d_n(\phi_{ijkl})
969 where $k$ is a force constant in units of energy
970 and $f$ is a cubic spline function; for details see \ssecref{cubicspline}.
971 For each interaction, the force constant $k$ and the table number $n$
972 are specified in the topology.
973 There are two different types of bonds, one that generates exclusions (type 8)
974 and one that does not (type 9).
975 For details see \tabref{topfile2}.
976 The table files are supplied to the {\tt mdrun} program.
977 After the table file name an underscore, the letter ``b'' for bonds,
978 ``a'' for angles or ``d'' for dihedrals and the table number are appended.
979 For example, for a bond with $n=0$ (and using the default table file name)
980 the table is read from the file {\tt table_b0.xvg}. Multiple tables can be
981 supplied simply by using different values of $n$, and are applied to the appropriate
982 bonds, as specified in the topology (\tabref{topfile2}).
983 The format for the table files is three columns with $x$, $f(x)$, $-f'(x)$,
984 where $x$ should be uniformly-spaced. Requirements for entries in the topology
985 are given in~\tabref{topfile2}.
986 The setup of the tables is as follows:
988 $x$ is the distance in nm. For distances beyond the table length,
989 {\tt mdrun} will quit with an error message.
991 $x$ is the angle in degrees. The table should go from
992 0 up to and including 180 degrees; the derivative is taken in degrees.
994 $x$ is the dihedral angle in degrees. The table should go from
995 -180 up to and including 180 degrees;
996 the IUPAC/IUB convention is used, {\ie} zero is cis,
997 the derivative is taken in degrees.
998 %} % Brace matches ifthenelse test for gmxlite
1000 \section{Restraints}
1001 Special potentials are used for imposing restraints on the motion of
1002 the system, either to avoid disastrous deviations, or to include
1003 knowledge from experimental data. In either case they are not really
1004 part of the force field and the reliability of the parameters is not
1005 important. The potential forms, as implemented in {\gromacs}, are
1006 mentioned just for the sake of completeness. Restraints and constraints
1007 refer to quite different algorithms in {\gromacs}.
1009 \subsection{Position restraints\swapindexquiet{position}{restraint}}
1010 \label{subsec:positionrestraint}
1011 These are used to restrain particles to fixed reference positions
1012 $\ve{R}_i$. They can be used during equilibration in order to avoid
1013 drastic rearrangements of critical parts ({\eg} to restrain motion
1014 in a protein that is subjected to large solvent forces when the
1015 solvent is not yet equilibrated). Another application is the
1016 restraining of particles in a shell around a region that is simulated
1017 in detail, while the shell is only approximated because it lacks
1018 proper interaction from missing particles outside the
1019 shell. Restraining will then maintain the integrity of the inner
1020 part. For spherical shells, it is a wise procedure to make the force
1021 constant depend on the radius, increasing from zero at the inner
1022 boundary to a large value at the outer boundary. This feature has
1023 not, however, been implemented in {\gromacs}.
1024 \newcommand{\unitv}[1]{\hat{\bf #1}}
1025 \newcommand{\halfje}[1]{\frac{#1}{2}}
1027 The following form is used:
1029 V_{pr}(\ve{r}_i) = \halfje{1}k_{pr}|\rvi-\ve{R}_i|^2
1031 The potential is plotted in \figref{positionrestraint}.
1034 \centerline{\includegraphics[width=8cm]{plots/f-pr}}
1035 \caption{Position restraint potential.}
1036 \label{fig:positionrestraint}
1039 The potential form can be rewritten without loss of generality as:
1041 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]
1047 F_i^x &=& -k_{pr}^x~(x_i - X_i) \\
1048 F_i^y &=& -k_{pr}^y~(y_i - Y_i) \\
1049 F_i^z &=& -k_{pr}^z~(z_i - Z_i)
1052 Using three different force constants the position
1053 restraints can be turned on or off
1054 in each spatial dimension; this means that atoms can be harmonically
1055 restrained to a plane or a line.
1056 Position restraints are applied to a special fixed list of atoms. Such a
1057 list is usually generated by the {\tt \normindex{pdb2gmx}} program.
1059 \subsection{\swapindex{Flat-bottomed}{position restraint}s}
1060 \label{subsec:fbpositionrestraint}
1061 Flat-bottomed position restraints can be used to restrain particles to
1062 part of the simulation volume. No force acts on the restrained
1063 particle within the flat-bottomed region of the potential, however a
1064 harmonic force acts to move the particle to the flat-bottomed region
1065 if it is outside it. It is possible to apply normal and
1066 flat-bottomed position restraints on the same particle (however, only
1067 with the same reference position $\ve{R}_i$). The following general potential
1068 is used (Figure~\ref{fig:fbposres}A):
1070 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}],
1072 where $\ve{R}_i$ is the reference position, $r_\mathrm{fb}$ is the distance
1073 from the center with a flat potential, $k_\mathrm{fb}$ the force constant, and $H$ is the Heaviside step
1074 function. The distance $d_g(\ve{r}_i;\ve{R}_i)$ from the reference
1075 position depends on the geometry $g$ of the flat-bottomed potential.
1078 \centerline{\includegraphics[width=10cm]{plots/fbposres}}
1079 \caption{Flat-bottomed position restraint potential. (A) Not
1080 inverted, (B) inverted.}
1081 \label{fig:fbposres}
1084 The following geometries for the flat-bottomed potential are supported:\newline
1085 {\bfseries Sphere} ($g =1$): The particle is kept in a sphere of given
1086 radius. The force acts towards the center of the sphere. The following distance calculation is used:
1088 d_g(\ve{r}_i;\ve{R}_i) = |\ve{r}_i-\ve{R}_i|
1090 {\bfseries Cylinder} ($g=6,7,8$): The particle is kept in a cylinder of given radius
1091 parallel to the $x$ ($g=6$), $y$ ($g=7$), or $z$-axis ($g=8$). For backwards compatibility, setting
1092 $g=2$ is mapped to $g=8$ in the code so that old {\tt .tpr} files and topologies work.
1093 The force from the flat-bottomed potential acts towards the axis of the cylinder.
1094 The component of the force parallel to the cylinder axis is zero.
1095 For a cylinder aligned along the $z$-axis:
1097 d_g(\ve{r}_i;\ve{R}_i) = \sqrt{ (x_i-X_i)^2 + (y_i - Y_i)^2 }
1099 {\bfseries Layer} ($g=3,4,5$): The particle is kept in a layer defined by the
1100 thickness and the normal of the layer. The layer normal can be parallel to the $x$, $y$, or
1101 $z$-axis. The force acts parallel to the layer normal.\\
1103 d_g(\ve{r}_i;\ve{R}_i) = |x_i-X_i|, \;\;\;\mbox{or}\;\;\;
1104 d_g(\ve{r}_i;\ve{R}_i) = |y_i-Y_i|, \;\;\;\mbox{or}\;\;\;
1105 d_g(\ve{r}_i;\ve{R}_i) = |z_i-Z_i|.
1108 It is possible to apply multiple independent flat-bottomed position
1109 restraints of different geometry on one particle. For example, applying
1110 a cylinder and a layer in $z$ keeps a particle within a
1111 disk. Applying three layers in $x$, $y$, and $z$ keeps the particle within a cuboid.
1113 In addition, it is possible to invert the restrained region with the
1114 unrestrained region, leading to a potential that acts to keep the particle {\it outside} of the volume
1115 defined by $\ve{R}_i$, $g$, and $r_\mathrm{fb}$. That feature is
1116 switched on by defining a negative $r_\mathrm{fb}$ in the
1117 topology. The following potential is used (Figure~\ref{fig:fbposres}B):
1119 V_\mathrm{fb}^{\mathrm{inv}}(\ve{r}_i) = \frac{1}{2}k_\mathrm{fb}
1120 [d_g(\ve{r}_i;\ve{R}_i) - |r_\mathrm{fb}|]^2\,
1121 H[ -(d_g(\ve{r}_i;\ve{R}_i) - |r_\mathrm{fb}|)].
1126 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1127 \subsection{Angle restraints\swapindexquiet{angle}{restraint}}
1128 \label{subsec:anglerestraint}
1129 These are used to restrain the angle between two pairs of particles
1130 or between one pair of particles and the $z$-axis.
1131 The functional form is similar to that of a proper dihedral.
1132 For two pairs of atoms:
1134 V_{ar}(\ve{r}_i,\ve{r}_j,\ve{r}_k,\ve{r}_l)
1135 = k_{ar}(1 - \cos(n (\theta - \theta_0))
1138 \theta = \arccos\left(\frac{\ve{r}_j -\ve{r}_i}{\|\ve{r}_j -\ve{r}_i\|}
1139 \cdot \frac{\ve{r}_l -\ve{r}_k}{\|\ve{r}_l -\ve{r}_k\|} \right)
1141 For one pair of atoms and the $z$-axis:
1143 V_{ar}(\ve{r}_i,\ve{r}_j) = k_{ar}(1 - \cos(n (\theta - \theta_0))
1146 \theta = \arccos\left(\frac{\ve{r}_j -\ve{r}_i}{\|\ve{r}_j -\ve{r}_i\|}
1147 \cdot \left( \begin{array}{c} 0 \\ 0 \\ 1 \\ \end{array} \right) \right)
1149 A multiplicity ($n$) of 2 is useful when you do not want to distinguish
1150 between parallel and anti-parallel vectors.
1151 The equilibrium angle $\theta$ should be between 0 and 180 degrees
1152 for multiplicity 1 and between 0 and 90 degrees for multiplicity 2.
1155 \subsection{Dihedral restraints\swapindexquiet{dihedral}{restraint}}
1156 \label{subsec:dihedralrestraint}
1157 These are used to restrain the dihedral angle $\phi$ defined by four particles
1158 as in an improper dihedral (sec.~\ref{sec:imp}) but with a slightly
1159 modified potential. Using:
1161 \phi' = \left(\phi-\phi_0\right) ~{\rm MOD}~ 2\pi
1164 where $\phi_0$ is the reference angle, the potential is defined as:
1166 V_{dihr}(\phi') ~=~ \left\{
1167 \begin{array}{lcllll}
1168 \half k_{dihr}(\phi'-\phi_0-\Delta\phi)^2
1169 &\mbox{for}& \phi' & > & \Delta\phi \\[1.5ex]
1170 0 &\mbox{for}& \phi' & \le & \Delta\phi \\[1.5ex]
1174 where $\Delta\phi$ is a user defined angle and $k_{dihr}$ is the force
1176 {\bf Note} that in the input in topology files, angles are given in degrees and
1177 force constants in kJ/mol/rad$^2$.
1178 %} % Brace matches ifthenelse test for gmxlite
1180 \subsection{Distance restraints\swapindexquiet{distance}{restraint}}
1181 \label{subsec:distancerestraint}
1183 add a penalty to the potential when the distance between specified
1184 pairs of atoms exceeds a threshold value. They are normally used to
1185 impose experimental restraints from, for instance, experiments in nuclear
1186 magnetic resonance (NMR), on the motion of the system. Thus, MD can be
1187 used for structure refinement using NMR data\index{nmr
1188 refinement}\index{refinement,nmr}.
1189 In {\gromacs} there are three ways to impose restraints on pairs of atoms:
1191 \item Simple harmonic restraints: use {\tt [ bonds ]} type 6
1192 %\ifthenelse{\equal{\gmxlite}{1}}
1194 {(see \secref{excl}).}
1195 \item\label{subsec:harmonicrestraint}Piecewise linear/harmonic restraints: {\tt [ bonds ]} type 10.
1196 \item Complex NMR distance restraints, optionally with pair, time and/or
1199 The last two options will be detailed now.
1201 The potential form for distance restraints is quadratic below a specified
1202 lower bound and between two specified upper bounds, and linear beyond the
1203 largest bound (see \figref{dist}).
1205 V_{dr}(r_{ij}) ~=~ \left\{
1206 \begin{array}{lcllllll}
1207 \half k_{dr}(r_{ij}-r_0)^2
1208 &\mbox{for}& & & r_{ij} & < & r_0 \\[1.5ex]
1209 0 &\mbox{for}& r_0 & \le & r_{ij} & < & r_1 \\[1.5ex]
1210 \half k_{dr}(r_{ij}-r_1)^2
1211 &\mbox{for}& r_1 & \le & r_{ij} & < & r_2 \\[1.5ex]
1212 \half k_{dr}(r_2-r_1)(2r_{ij}-r_2-r_1)
1213 &\mbox{for}& r_2 & \le & r_{ij} & &
1219 \centerline{\includegraphics[width=8cm]{plots/f-dr}}
1220 \caption{Distance Restraint potential.}
1227 \begin{array}{lcllllll}
1228 -k_{dr}(r_{ij}-r_0)\frac{\rvij}{r_{ij}}
1229 &\mbox{for}& & & r_{ij} & < & r_0 \\[1.5ex]
1230 0 &\mbox{for}& r_0 & \le & r_{ij} & < & r_1 \\[1.5ex]
1231 -k_{dr}(r_{ij}-r_1)\frac{\rvij}{r_{ij}}
1232 &\mbox{for}& r_1 & \le & r_{ij} & < & r_2 \\[1.5ex]
1233 -k_{dr}(r_2-r_1)\frac{\rvij}{r_{ij}}
1234 &\mbox{for}& r_2 & \le & r_{ij} & &
1238 For restraints not derived from NMR data, this functionality
1239 will usually suffice and a section of {\tt [ bonds ]} type 10
1240 can be used to apply individual restraints between pairs of
1241 %\ifthenelse{\equal{\gmxlite}{1}}{atoms.}{
1242 atoms, see \ssecref{topfile}.
1243 %} % Brace matches ifthenelse test for gmxlite
1244 For applying restraints derived from NMR measurements, more complex
1245 functionality might be required, which is provided through
1246 the {\tt [~distance_restraints~]} section and is described below.
1248 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1249 \subsubsection{Time averaging\swapindexquiet{time-averaged}{distance restraint}}
1250 Distance restraints based on instantaneous distances can potentially reduce
1251 the fluctuations in a molecule significantly. This problem can be overcome by restraining
1252 to a {\em time averaged} distance~\cite{Torda89}.
1253 The forces with time averaging are:
1256 \begin{array}{lcllllll}
1257 -k^a_{dr}(\bar{r}_{ij}-r_0)\frac{\rvij}{r_{ij}}
1258 &\mbox{for}& & & \bar{r}_{ij} & < & r_0 \\[1.5ex]
1259 0 &\mbox{for}& r_0 & \le & \bar{r}_{ij} & < & r_1 \\[1.5ex]
1260 -k^a_{dr}(\bar{r}_{ij}-r_1)\frac{\rvij}{r_{ij}}
1261 &\mbox{for}& r_1 & \le & \bar{r}_{ij} & < & r_2 \\[1.5ex]
1262 -k^a_{dr}(r_2-r_1)\frac{\rvij}{r_{ij}}
1263 &\mbox{for}& r_2 & \le & \bar{r}_{ij} & &
1266 where $\bar{r}_{ij}$ is given by an exponential running average with decay time $\tau$:
1268 \bar{r}_{ij} ~=~ < r_{ij}^{-3} >^{-1/3}
1271 The force constant $k^a_{dr}$ is switched on slowly to compensate for
1272 the lack of history at the beginning of the simulation:
1274 k^a_{dr} = k_{dr} \left(1-\exp\left(-\frac{t}{\tau}\right)\right)
1276 Because of the time averaging, we can no longer speak of a distance restraint
1279 This way an atom can satisfy two incompatible distance restraints
1280 {\em on average} by moving between two positions.
1281 An example would be an amino acid side-chain that is rotating around
1282 its $\chi$ dihedral angle, thereby coming close to various other groups.
1283 Such a mobile side chain can give rise to multiple NOEs that can not be
1284 fulfilled by a single structure.
1286 The computation of the time
1287 averaged distance in the {\tt mdrun} program is done in the following fashion:
1290 \overline{r^{-3}}_{ij}(0) &=& r_{ij}(0)^{-3} \\
1291 \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]
1292 \label{eqn:ravdisre}
1296 When a pair is within the bounds, it can still feel a force
1297 because the time averaged distance can still be beyond a bound.
1298 To prevent the protons from being pulled too close together, a mixed
1299 approach can be used. In this approach, the penalty is zero when the
1300 instantaneous distance is within the bounds, otherwise the violation is
1301 the square root of the product of the instantaneous violation and the
1302 time averaged violation:
1305 \begin{array}{lclll}
1306 k^a_{dr}\sqrt{(r_{ij}-r_0)(\bar{r}_{ij}-r_0)}\frac{\rvij}{r_{ij}}
1307 & \mbox{for} & r_{ij} < r_0 & \mbox{and} & \bar{r}_{ij} < r_0 \\[1.5ex]
1309 \mbox{min}\left(\sqrt{(r_{ij}-r_1)(\bar{r}_{ij}-r_1)},r_2-r_1\right)
1310 \frac{\rvij}{r_{ij}}
1311 & \mbox{for} & r_{ij} > r_1 & \mbox{and} & \bar{r}_{ij} > r_1 \\[1.5ex]
1316 \subsubsection{Averaging over multiple pairs\swapindexquiet{ensemble-averaged}{distance restraint}}
1318 Sometimes it is unclear from experimental data which atom pair
1319 gives rise to a single NOE, in other occasions it can be obvious that
1320 more than one pair contributes due to the symmetry of the system, {\eg} a
1321 methyl group with three protons. For such a group, it is not possible
1322 to distinguish between the protons, therefore they should all be taken into
1323 account when calculating the distance between this methyl group and another
1324 proton (or group of protons).
1325 Due to the physical nature of magnetic resonance, the intensity of the
1326 NOE signal is inversely proportional to the sixth power of the inter-atomic
1328 Thus, when combining atom pairs,
1329 a fixed list of $N$ restraints may be taken together,
1330 where the apparent ``distance'' is given by:
1332 r_N(t) = \left [\sum_{n=1}^{N} \bar{r}_{n}(t)^{-6} \right]^{-1/6}
1335 where we use $r_{ij}$ or \eqnref{rav} for the $\bar{r}_{n}$.
1336 The $r_N$ of the instantaneous and time-averaged distances
1337 can be combined to do a mixed restraining, as indicated above.
1338 As more pairs of protons contribute to the same NOE signal, the intensity
1339 will increase, and the summed ``distance'' will be shorter than any of
1340 its components due to the reciprocal summation.
1342 There are two options for distributing the forces over the atom pairs.
1343 In the conservative option, the force is defined as the derivative of the
1344 restraint potential with respect to the coordinates. This results in
1345 a conservative potential when time averaging is not used.
1346 The force distribution over the pairs is proportional to $r^{-6}$.
1347 This means that a close pair feels a much larger force than a distant pair,
1348 which might lead to a molecule that is ``too rigid.''
1349 The other option is an equal force distribution. In this case each pair
1350 feels $1/N$ of the derivative of the restraint potential with respect to
1351 $r_N$. The advantage of this method is that more conformations might be
1352 sampled, but the non-conservative nature of the forces can lead to
1353 local heating of the protons.
1355 It is also possible to use {\em ensemble averaging} using multiple
1356 (protein) molecules. In this case the bounds should be lowered as in:
1359 r_1 &~=~& r_1 * M^{-1/6} \\
1360 r_2 &~=~& r_2 * M^{-1/6}
1363 where $M$ is the number of molecules. The {\gromacs} preprocessor {\tt grompp}
1364 can do this automatically when the appropriate option is given.
1365 The resulting ``distance'' is
1366 then used to calculate the scalar force according to:
1370 ~& 0 \hspace{4cm} & r_{N} < r_1 \\
1371 & k_{dr}(r_{N}-r_1)\frac{\rvij}{r_{ij}} & r_1 \le r_{N} < r_2 \\
1372 & k_{dr}(r_2-r_1)\frac{\rvij}{r_{ij}} & r_{N} \ge r_2
1375 where $i$ and $j$ denote the atoms of all the
1376 pairs that contribute to the NOE signal.
1378 \subsubsection{Using distance restraints}
1380 A list of distance restrains based on NOE data can be added to a molecule
1381 definition in your topology file, like in the following example:
1384 [ distance_restraints ]
1385 ; ai aj type index type' low up1 up2 fac
1386 10 16 1 0 1 0.0 0.3 0.4 1.0
1387 10 28 1 1 1 0.0 0.3 0.4 1.0
1388 10 46 1 1 1 0.0 0.3 0.4 1.0
1389 16 22 1 2 1 0.0 0.3 0.4 2.5
1390 16 34 1 3 1 0.0 0.5 0.6 1.0
1393 In this example a number of features can be found. In columns {\tt
1394 ai} and {\tt aj} you find the atom numbers of the particles to be
1395 restrained. The {\tt type} column should always be 1. As explained in
1396 ~\ssecref{distancerestraint}, multiple distances can contribute to a single NOE
1397 signal. In the topology this can be set using the {\tt index}
1398 column. In our example, the restraints 10-28 and 10-46 both have index
1399 1, therefore they are treated simultaneously. An extra requirement
1400 for treating restraints together is that the restraints must be on
1401 successive lines, without any other intervening restraint. The {\tt
1402 type'} column will usually be 1, but can be set to 2 to obtain a
1403 distance restraint that will never be time- and ensemble-averaged;
1404 this can be useful for restraining hydrogen bonds. The columns {\tt
1405 low}, {\tt up1}, and {\tt up2} hold the values of $r_0$, $r_1$, and
1406 $r_2$ from ~\eqnref{disre}. In some cases it can be useful to have
1407 different force constants for some restraints; this is controlled by
1408 the column {\tt fac}. The force constant in the parameter file is
1409 multiplied by the value in the column {\tt fac} for each restraint.
1410 %} % Brace matches ifthenelse test for gmxlite
1412 \newcommand{\SSS}{{\mathbf S}}
1413 \newcommand{\DD}{{\mathbf D}}
1414 \newcommand{\RR}{{\mathbf R}}
1416 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1417 \subsection{Orientation restraints\swapindexquiet{orientation}{restraint}}
1418 \label{subsec:orientationrestraint}
1419 This section describes how orientations between vectors,
1420 as measured in certain NMR experiments, can be calculated
1421 and restrained in MD simulations.
1422 The presented refinement methodology and a comparison of results
1423 with and without time and ensemble averaging have been
1424 published~\cite{Hess2003}.
1425 \subsubsection{Theory}
1426 In an NMR experiment, orientations of vectors can be measured when a
1427 molecule does not tumble completely isotropically in the solvent.
1428 Two examples of such orientation measurements are
1429 residual \normindex{dipolar couplings}
1430 (between two nuclei) or chemical shift anisotropies.
1431 An observable for a vector $\ve{r}_i$ can be written as follows:
1433 \delta_i = \frac{2}{3} \mbox{tr}(\SSS\DD_i)
1435 where $\SSS$ is the dimensionless order tensor of the molecule.
1436 The tensor $\DD_i$ is given by:
1439 \DD_i = \frac{c_i}{\|\ve{r}_i\|^\alpha} \left(
1441 %3 r_x r_x - \ve{r}\cdot\ve{r} & 3 r_x r_y & 3 r_x r_z \\
1442 %3 r_x r_y & 3 r_y r_y - \ve{r}\cdot\ve{r} & 3yz \\
1443 %3 r_x r_z & 3 r_y r_z & 3 r_z r_z - \ve{r}\cdot\ve{r}
1444 %\end{array} \right)
1446 3 x x - 1 & 3 x y & 3 x z \\
1447 3 x y & 3 y y - 1 & 3 y z \\
1448 3 x z & 3 y z & 3 z z - 1 \\
1453 x=\frac{r_{i,x}}{\|\ve{r}_i\|}, \quad
1454 y=\frac{r_{i,y}}{\|\ve{r}_i\|}, \quad
1455 z=\frac{r_{i,z}}{\|\ve{r}_i\|}
1457 For a dipolar coupling $\ve{r}_i$ is the vector connecting the two
1458 nuclei, $\alpha=3$ and the constant $c_i$ is given by:
1460 c_i = \frac{\mu_0}{4\pi} \gamma_1^i \gamma_2^i \frac{\hbar}{4\pi}
1462 where $\gamma_1^i$ and $\gamma_2^i$ are the gyromagnetic ratios of the
1465 The order tensor is symmetric and has trace zero. Using a rotation matrix
1466 ${\mathbf T}$ it can be transformed into the following form:
1468 {\mathbf T}^T \SSS {\mathbf T} = s \left( \begin{array}{ccc}
1469 -\frac{1}{2}(1-\eta) & 0 & 0 \\
1470 0 & -\frac{1}{2}(1+\eta) & 0 \\
1474 where $-1 \leq s \leq 1$ and $0 \leq \eta \leq 1$.
1475 $s$ is called the order parameter and $\eta$ the asymmetry of the
1476 order tensor $\SSS$. When the molecule tumbles isotropically in the
1477 solvent, $s$ is zero, and no orientational effects can be observed
1478 because all $\delta_i$ are zero.
1482 \subsubsection{Calculating orientations in a simulation}
1483 For reasons which are explained below, the $\DD$ matrices are calculated
1484 which respect to a reference orientation of the molecule. The orientation
1485 is defined by a rotation matrix $\RR$, which is needed to least-squares fit
1486 the current coordinates of a selected set of atoms onto
1487 a reference conformation. The reference conformation is the starting
1488 conformation of the simulation. In case of ensemble averaging, which will
1489 be treated later, the structure is taken from the first subsystem.
1490 The calculated $\DD_i^c$ matrix is given by:
1493 \DD_i^c(t) = \RR(t) \DD_i(t) \RR^T(t)
1495 The calculated orientation for vector $i$ is given by:
1497 \delta^c_i(t) = \frac{2}{3} \mbox{tr}(\SSS(t)\DD_i^c(t))
1499 The order tensor $\SSS(t)$ is usually unknown.
1500 A reasonable choice for the order tensor is the tensor
1501 which minimizes the (weighted) mean square difference between the calculated
1502 and the observed orientations:
1505 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
1507 To properly combine different types of measurements, the unit of $w_i$ should
1508 be such that all terms are dimensionless. This means the unit of $w_i$
1509 is the unit of $\delta_i$ to the power $-2$.
1510 {\bf Note} that scaling all $w_i$ with a constant factor does not influence
1513 \subsubsection{Time averaging}
1514 Since the tensors $\DD_i$ fluctuate rapidly in time, much faster than can
1515 be observed in an experiment, they should be averaged over time in the simulation.
1516 However, in a simulation the time and the number of copies of
1517 a molecule are limited. Usually one can not obtain a converged average
1518 of the $\DD_i$ tensors over all orientations of the molecule.
1519 If one assumes that the average orientations of the $\ve{r}_i$ vectors
1520 within the molecule converge much faster than the tumbling time of
1521 the molecule, the tensor can be averaged in an axis system that
1522 rotates with the molecule, as expressed by equation~(\ref{D_rot}).
1523 The time-averaged tensors are calculated
1524 using an exponentially decaying memory function:
1526 \DD^a_i(t) = \frac{\displaystyle
1527 \int_{u=t_0}^t \DD^c_i(u) \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
1529 \int_{u=t_0}^t \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
1532 Assuming that the order tensor $\SSS$ fluctuates slower than the
1533 $\DD_i$, the time-averaged orientation can be calculated as:
1535 \delta_i^a(t) = \frac{2}{3} \mbox{tr}(\SSS(t) \DD_i^a(t))
1537 where the order tensor $\SSS(t)$ is calculated using expression~(\ref{S_msd})
1538 with $\delta_i^c(t)$ replaced by $\delta_i^a(t)$.
1540 \subsubsection{Restraining}
1541 The simulated structure can be restrained by applying a force proportional
1542 to the difference between the calculated and the experimental orientations.
1543 When no time averaging is applied, a proper potential can be defined as:
1545 V = \frac{1}{2} k \sum_{i=1}^N w_i (\delta_i^c (t) -\delta_i^{exp})^2
1547 where the unit of $k$ is the unit of energy.
1548 Thus the effective force constant for restraint $i$ is $k w_i$.
1549 The forces are given by minus the gradient of $V$.
1550 The force $\ve{F}\!_i$ working on vector $\ve{r}_i$ is:
1553 & = & - \frac{\mbox{d} V}{\mbox{d}\ve{r}_i} \\
1554 & = & -k w_i (\delta_i^c (t) -\delta_i^{exp}) \frac{\mbox{d} \delta_i (t)}{\mbox{d}\ve{r}_i} \\
1555 & = & -k w_i (\delta_i^c (t) -\delta_i^{exp})
1556 \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)
1559 \subsubsection{Ensemble averaging}
1560 Ensemble averaging can be applied by simulating a system of $M$ subsystems
1561 that each contain an identical set of orientation restraints. The systems only
1562 interact via the orientation restraint potential which is defined as:
1564 V = M \frac{1}{2} k \sum_{i=1}^N w_i
1565 \langle \delta_i^c (t) -\delta_i^{exp} \rangle^2
1567 The force on vector $\ve{r}_{i,m}$ in subsystem $m$ is given by:
1569 \ve{F}\!_{i,m}(t) = - \frac{\mbox{d} V}{\mbox{d}\ve{r}_{i,m}} =
1570 -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}} \\
1573 \subsubsection{Time averaging}
1574 When using time averaging it is not possible to define a potential.
1575 We can still define a quantity that gives a rough idea of the energy
1576 stored in the restraints:
1578 V = M \frac{1}{2} k^a \sum_{i=1}^N w_i
1579 \langle \delta_i^a (t) -\delta_i^{exp} \rangle^2
1581 The force constant $k_a$ is switched on slowly to compensate for the
1582 lack of history at times close to $t_0$. It is exactly proportional
1583 to the amount of average that has been accumulated:
1586 k \, \frac{1}{\tau}\int_{u=t_0}^t \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
1588 What really matters is the definition of the force. It is chosen to
1589 be proportional to the square root of the product of the time-averaged
1590 and the instantaneous deviation.
1591 Using only the time-averaged deviation induces large oscillations.
1592 The force is given by:
1595 %\left\{ \begin{array}{ll}
1596 %0 & \mbox{for} \quad \langle \delta_i^a (t) -\delta_i^{exp} \rangle \langle \delta_i (t) -\delta_i^{exp} \rangle \leq 0 \\
1597 %... & \mbox{for} \quad \langle \delta_i^a (t) -\delta_i^{exp} \rangle \langle \delta_i (t) -\delta_i^{exp} \rangle > 0
1600 \left\{ \begin{array}{ll}
1601 0 & \quad \mbox{for} \quad a\, b \leq 0 \\
1603 k^a w_i \frac{a}{|a|} \sqrt{a\, b} \, \frac{\mbox{d} \delta_{i,m}^c (t)}{\mbox{d}\ve{r}_{i,m}}
1604 & \quad \mbox{for} \quad a\, b > 0
1609 a &=& \langle \delta_i^a (t) -\delta_i^{exp} \rangle \\
1610 b &=& \langle \delta_i^c (t) -\delta_i^{exp} \rangle
1613 \subsubsection{Using orientation restraints}
1614 Orientation restraints can be added to a molecule definition in
1615 the topology file in the section {\tt [~orientation_restraints~]}.
1616 Here we give an example section containing five N-H residual dipolar
1617 coupling restraints:
1620 [ orientation_restraints ]
1621 ; ai aj type exp. label alpha const. obs. weight
1623 31 32 1 1 3 3 6.083 -6.73 1.0
1624 43 44 1 1 4 3 6.083 -7.87 1.0
1625 55 56 1 1 5 3 6.083 -7.13 1.0
1626 65 66 1 1 6 3 6.083 -2.57 1.0
1627 73 74 1 1 7 3 6.083 -2.10 1.0
1630 The unit of the observable is Hz, but one can choose any other unit.
1632 ai} and {\tt aj} you find the atom numbers of the particles to be
1633 restrained. The {\tt type} column should always be 1.
1634 The {\tt exp.} column denotes the experiment number, starting
1635 at 1. For each experiment a separate order tensor $\SSS$
1636 is optimized. The label should be a unique number larger than zero
1637 for each restraint. The {\tt alpha} column contains the power $\alpha$
1638 that is used in equation~(\ref{orient_def}) to calculate the orientation.
1639 The {\tt const.} column contains the constant $c_i$ used in the same
1640 equation. The constant should have the unit of the observable times
1641 nm$^\alpha$. The column {\tt obs.} contains the observable, in any
1642 unit you like. The last column contains the weights $w_i$; the unit
1643 should be the inverse of the square of the unit of the observable.
1645 Some parameters for orientation restraints can be specified in the
1646 {\tt grompp.mdp} file, for a study of the effect of different
1647 force constants and averaging times and ensemble averaging see~\cite{Hess2003}.
1648 %} % Brace matches ifthenelse test for gmxlite
1650 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1651 \section{Polarization}
1652 Polarization can be treated by {\gromacs} by attaching
1653 \normindex{shell} (\normindex{Drude}) particles to atoms and/or
1654 virtual sites. The energy of the shell particle is then minimized at
1655 each time step in order to remain on the Born-Oppenheimer surface.
1657 \subsection{Simple polarization}
1658 This is merely a harmonic potential with equilibrium distance 0.
1660 \subsection{Water polarization}
1661 A special potential for water that allows anisotropic polarization of
1662 a single shell particle~\cite{Maaren2001a}.
1664 \subsection{Thole polarization}
1665 Based on early work by \normindex{Thole}~\cite{Thole81}, Roux and
1666 coworkers have implemented potentials for molecules like
1667 ethanol~\cite{Lamoureux2003a,Lamoureux2003b,Noskov2005a}. Within such
1668 molecules, there are intra-molecular interactions between shell
1669 particles, however these must be screened because full Coulomb would
1670 be too strong. The potential between two shell particles $i$ and $j$ is:
1671 \newcommand{\rbij}{\bar{r}_{ij}}
1673 V_{thole} ~=~ \frac{q_i q_j}{r_{ij}}\left[1-\left(1+\frac{\rbij}{2}\right){\rm exp}^{-\rbij}\right]
1675 {\bf Note} that there is a sign error in Equation~1 of Noskov {\em et al.}~\cite{Noskov2005a}:
1677 \rbij ~=~ a\frac{r_{ij}}{(\alpha_i \alpha_j)^{1/6}}
1679 where $a$ is a magic (dimensionless) constant, usually chosen to be
1680 2.6~\cite{Noskov2005a}; $\alpha_i$ and $\alpha_j$ are the polarizabilities
1681 of the respective shell particles.
1683 %} % Brace matches ifthenelse test for gmxlite
1685 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1686 \section{Free energy interactions}
1688 \index{free energy interactions}
1689 \newcommand{\LAM}{\lambda}
1690 \newcommand{\LL}{(1-\LAM)}
1691 \newcommand{\dvdl}[1]{\frac{\partial #1}{\partial \LAM}}
1692 This section describes the $\lambda$-dependence of the potentials
1693 used for free energy calculations (see \secref{fecalc}).
1694 All common types of potentials and constraints can be
1695 interpolated smoothly from state A ($\lambda=0$) to state B
1696 ($\lambda=1$) and vice versa.
1697 All bonded interactions are interpolated by linear interpolation
1698 of the interaction parameters. Non-bonded interactions can be
1699 interpolated linearly or via soft-core interactions.
1701 Starting in {\gromacs} 4.6, $\lambda$ is a vector, allowing different
1702 components of the free energy transformation to be carried out at
1703 different rates. Coulomb, Lennard-Jones, bonded, and restraint terms
1704 can all be controlled independently, as described in the {\tt .mdp}
1707 \subsubsection{Harmonic potentials}
1708 The example given here is for the bond potential, which is harmonic
1709 in {\gromacs}. However, these equations apply to the angle potential
1710 and the improper dihedral potential as well.
1712 V_b &=&\half\left[\LL k_b^A +
1713 \LAM k_b^B\right] \left[b - \LL b_0^A - \LAM b_0^B\right]^2 \\
1714 \dvdl{V_b}&=&\half(k_b^B-k_b^A)
1715 \left[b - \LL b_0^A + \LAM b_0^B\right]^2 +
1717 & & \phantom{\half}(b_0^A-b_0^B) \left[b - \LL b_0^A -\LAM b_0^B\right]
1718 \left[\LL k_b^A + \LAM k_b^B \right]
1721 \subsubsection{\gromosv{96} bonds and angles}
1722 Fourth-power bond stretching and cosine-based angle potentials
1723 are interpolated by linear interpolation of the force constant
1724 and the equilibrium position. Formulas are not given here.
1726 \subsubsection{Proper dihedrals}
1727 For the proper dihedrals, the equations are somewhat more complicated:
1729 V_d &=&\left[\LL k_d^A + \LAM k_d^B \right]
1730 \left( 1+ \cos\left[n_{\phi} \phi -
1731 \LL \phi_s^A - \LAM \phi_s^B
1733 \dvdl{V_d}&=&(k_d^B-k_d^A)
1736 n_{\phi} \phi- \LL \phi_s^A - \LAM \phi_s^B
1740 &&(\phi_s^B - \phi_s^A) \left[\LL k_d^A - \LAM k_d^B\right]
1741 \sin\left[ n_{\phi}\phi - \LL \phi_s^A - \LAM \phi_s^B \right]
1743 {\bf Note:} that the multiplicity $n_{\phi}$ can not be parameterized
1744 because the function should remain periodic on the interval $[0,2\pi]$.
1746 \subsubsection{Tabulated bonded interactions}
1747 For tabulated bonded interactions only the force constant can interpolated:
1749 V &=& (\LL k^A + \LAM k^B) \, f \\
1750 \dvdl{V} &=& (k^B - k^A) \, f
1753 \subsubsection{Coulomb interaction}
1754 The \normindex{Coulomb} interaction between two particles
1755 of which the charge varies with $\LAM$ is:
1757 V_c &=& \frac{f}{\epsrf \rij}\left[\LL q_i^A q_j^A + \LAM\, q_i^B q_j^B\right] \\
1758 \dvdl{V_c}&=& \frac{f}{\epsrf \rij}\left[- q_i^A q_j^A + q_i^B q_j^B\right]
1760 where $f = \frac{1}{4\pi \varepsilon_0} = 138.935\,485$ (see \chref{defunits}).
1762 \subsubsection{Coulomb interaction with \normindex{reaction field}}
1763 The Coulomb interaction including a reaction field, between two particles
1764 of which the charge varies with $\LAM$ is:
1766 V_c &=& f\left[\frac{1}{\rij} + k_{rf}~ \rij^2 -c_{rf}\right]
1767 \left[\LL q_i^A q_j^A + \LAM\, q_i^B q_j^B\right] \\
1768 \dvdl{V_c}&=& f\left[\frac{1}{\rij} + k_{rf}~ \rij^2 -c_{rf}\right]
1769 \left[- q_i^A q_j^A + q_i^B q_j^B\right]
1770 \label{eq:dVcoulombdlambda}
1772 {\bf Note} that the constants $k_{rf}$ and $c_{rf}$ are
1773 defined using the dielectric
1774 constant $\epsrf$ of the medium (see \secref{coulrf}).
1776 \subsubsection{Lennard-Jones interaction}
1777 For the \normindex{Lennard-Jones} interaction between two particles
1778 of which the {\em atom type} varies with $\LAM$ we can write:
1780 V_{LJ} &=& \frac{\LL C_{12}^A + \LAM\, C_{12}^B}{\rij^{12}} -
1781 \frac{\LL C_6^A + \LAM\, C_6^B}{\rij^6} \\
1782 \dvdl{V_{LJ}}&=&\frac{C_{12}^B - C_{12}^A}{\rij^{12}} -
1783 \frac{C_6^B - C_6^A}{\rij^6}
1784 \label{eq:dVljdlambda}
1786 It should be noted that it is also possible to express a pathway from
1787 state A to state B using $\sigma$ and $\epsilon$ (see \eqnref{sigeps}).
1788 It may seem to make sense physically to vary the force field parameters
1789 $\sigma$ and $\epsilon$ rather
1790 than the derived parameters $C_{12}$ and $C_{6}$.
1791 However, the difference between the pathways in parameter space
1792 is not large, and the free energy itself
1793 does not depend on the pathway, so we use the simple formulation
1796 \subsubsection{Kinetic Energy}
1797 When the mass of a particle changes, there is also a contribution of
1798 the kinetic energy to the free energy (note that we can not write
1799 the momentum \ve{p} as m\ve{v}, since that would result
1800 in the sign of $\dvdl{E_k}$ being incorrect~\cite{Gunsteren98a}):
1803 E_k &=& \half\frac{\ve{p}^2}{\LL m^A + \LAM m^B} \\
1804 \dvdl{E_k}&=& -\half\frac{\ve{p}^2(m^B-m^A)}{(\LL m^A + \LAM m^B)^2}
1806 after taking the derivative, we {\em can} insert \ve{p} = m\ve{v}, such that:
1808 \dvdl{E_k}~=~ -\half\ve{v}^2(m^B-m^A)
1811 \subsubsection{Constraints}
1812 \label{subsubsec:constraints}
1813 \newcommand{\clam}{C_{\lambda}}
1814 The constraints are formally part of the Hamiltonian, and therefore
1815 they give a contribution to the free energy. In {\gromacs} this can be
1816 calculated using the \normindex{LINCS} or the \normindex{SHAKE} algorithm.
1817 If we have a number of constraint equations $g_k$:
1819 g_k = \ve{r}_{k} - d_{k}
1821 where $\ve{r}_k$ is the distance vector between two particles and
1822 $d_k$ is the constraint distance between the two particles, we can write
1823 this using a $\LAM$-dependent distance as
1825 g_k = \ve{r}_{k} - \left(\LL d_{k}^A + \LAM d_k^B\right)
1827 the contribution $\clam$
1828 to the Hamiltonian using Lagrange multipliers $\lambda$:
1830 \clam &=& \sum_k \lambda_k g_k \\
1831 \dvdl{\clam} &=& \sum_k \lambda_k \left(d_k^B-d_k^A\right)
1835 \subsection{Soft-core interactions\index{soft-core interactions}}
1837 \centerline{\includegraphics[height=6cm]{plots/softcore}}
1838 \caption{Soft-core interactions at $\LAM=0.5$, with $p=2$ and
1839 $C_6^A=C_{12}^A=C_6^B=C_{12}^B=1$.}
1840 \label{fig:softcore}
1842 In a free-energy calculation where particles grow out of nothing, or
1843 particles disappear, using the the simple linear interpolation of the
1844 Lennard-Jones and Coulomb potentials as described in Equations~\ref{eq:dVljdlambda}
1845 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
1846 close to each other, leading to large fluctuations in the measured values of
1847 $\partial V/\partial \LAM$ (which, because of the simple linear
1848 interpolation, depends on the potentials at both the endpoints of $\LAM$).
1850 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
1851 involved at $\LAM$ values between 0 and 1, but not \emph{at} $\LAM=0$
1854 In {\gromacs} the soft-core potentials $V_{sc}$ are shifted versions of the
1855 regular potentials, so that the singularity in the potential and its
1856 derivatives at $r=0$ is never reached:
1858 V_{sc}(r) &=& \LL V^A(r_A) + \LAM V^B(r_B)
1860 r_A &=& \left(\alpha \sigma_A^6 \LAM^p + r^6 \right)^\frac{1}{6}
1862 r_B &=& \left(\alpha \sigma_B^6 \LL^p + r^6 \right)^\frac{1}{6}
1864 where $V^A$ and $V^B$ are the normal ``hard core'' Van der Waals or
1865 electrostatic potentials in state A ($\LAM=0$) and state B ($\LAM=1$)
1866 respectively, $\alpha$ is the soft-core parameter (set with {\tt sc_alpha}
1867 in the {\tt .mdp} file), $p$ is the soft-core $\LAM$ power (set with
1868 {\tt sc_power}), $\sigma$ is the radius of the interaction, which is
1869 $(C_{12}/C_6)^{1/6}$ or an input parameter ({\tt sc_sigma}) when $C_6$
1870 or $C_{12}$ is zero.
1872 For intermediate $\LAM$, $r_A$ and $r_B$ alter the interactions very little
1873 for $r > \alpha^{1/6} \sigma$ and quickly switch the soft-core
1874 interaction to an almost constant value for smaller $r$ (\figref{softcore}).
1877 F_{sc}(r) = -\frac{\partial V_{sc}(r)}{\partial r} =
1878 \LL F^A(r_A) \left(\frac{r}{r_A}\right)^5 +
1879 \LAM F^B(r_B) \left(\frac{r}{r_B}\right)^5
1881 where $F^A$ and $F^B$ are the ``hard core'' forces.
1882 The contribution to the derivative of the free energy is:
1884 \dvdl{V_{sc}(r)} & = &
1885 V^B(r_B) -V^A(r_A) +
1886 \LL \frac{\partial V^A(r_A)}{\partial r_A}
1887 \frac{\partial r_A}{\partial \LAM} +
1888 \LAM\frac{\partial V^B(r_B)}{\partial r_B}
1889 \frac{\partial r_B}{\partial \LAM}
1892 V^B(r_B) -V^A(r_A) + \nonumber \\
1895 \left[ \LAM F^B(r_B) r^{-5}_B \sigma_B^6 \LL^{p-1} -
1896 \LL F^A(r_A) r^{-5}_A \sigma_A^6 \LAM^{p-1} \right]
1899 The original GROMOS Lennard-Jones soft-core function~\cite{Beutler94}
1900 uses $p=2$, but $p=1$ gives a smoother $\partial H/\partial\LAM$ curve.
1901 %When the changes between the two states involve both the disappearing
1902 %and appearing of atoms, it is important that the overlapping of atoms
1903 %happens around $\LAM=0.5$. This can usually be achieved with
1904 %$\alpha$$\approx0.7$ for $p=1$ and $\alpha$$\approx1.5$ for $p=2$.
1905 %MRS: this is now eliminated as of 4.6, since changes between atoms are done linearly.
1907 Another issue that should be considered is the soft-core effect of hydrogens
1908 without Lennard-Jones interaction. Their soft-core $\sigma$ is
1909 set with {\tt sc-sigma} in the {\tt .mdp} file. These hydrogens
1910 produce peaks in $\partial H/\partial\LAM$ at $\LAM$ is 0 and/or 1 for $p=1$
1911 and close to 0 and/or 1 with $p=2$. Lowering {\tt\mbox{sc-sigma}} will decrease
1912 this effect, but it will also increase the interactions with hydrogens
1913 relative to the other interactions in the soft-core state.
1915 When soft core potentials are selected (by setting {\tt sc-alpha} \textgreater
1916 0), and the Coulomb and Lennard-Jones potentials are turned on or off
1917 sequentially, then the Coulombic interaction is turned off linearly,
1918 rather than using soft core interactions, which should be less
1919 statistically noisy in most cases. This behavior can be overwritten
1920 by using the mdp option {\tt sc-coul} to {\tt yes}. Additionally, the
1921 soft-core interaction potential is only applied when either the A or B
1922 state has zero interaction potential. If both A and B states have
1923 nonzero interaction potential, default linear scaling described above
1924 is used. When both Coulombic and Lennard-Jones interactions are turned
1925 off simultaneously, a soft-core potential is used, and a hydrogen is
1926 being introduced or deleted, the sigma is set to {\tt sc-sigma-min},
1927 which itself defaults to {\tt sc-sigma-default}.
1929 Recently, a new formulation of the soft-core approach has been derived
1930 that in most cases gives lower and more even statistical variance than
1931 the standard soft-core path described above.~\cite{Pham2011,Pham2012}
1932 Specifically, we have:
1934 V_{sc}(r) &=& \LL V^A(r_A) + \LAM V^B(r_B)
1936 r_A &=& \left(\alpha \sigma_A^{48} \LAM^p + r^{48} \right)^\frac{1}{48}
1938 r_B &=& \left(\alpha \sigma_B^{48} \LL^p + r^{48} \right)^\frac{1}{48}
1940 This ``1-1-48'' path is also implemented in {\gromacs}. Note that for this path the soft core $\alpha$
1941 should satisfy $0.001 < \alpha < 0.003$, rather than $\alpha \approx
1944 %} % Brace matches ifthenelse test for gmxlite
1946 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1948 \subsection{Exclusions and 1-4 Interactions.}
1949 Atoms within a molecule that are close by in the chain,
1950 {\ie} atoms that are covalently bonded, or linked by one or two
1951 atoms are called {\em first neighbors, second neighbors} and
1952 {\em \swapindex{third}{neighbor}s}, respectively (see \figref{chain}).
1953 Since the interactions of atom {\bf i} with atoms {\bf i+1} and {\bf i+2}
1954 are mainly quantum mechanical, they can not be modeled by a Lennard-Jones potential.
1955 Instead it is assumed that these interactions are adequately modeled
1956 by a harmonic bond term or constraint ({\bf i, i+1}) and a harmonic angle term
1957 ({\bf i, i+2}). The first and second neighbors (atoms {\bf i+1} and {\bf i+2})
1959 {\em excluded} from the Lennard-Jones \swapindex{interaction}{list}
1961 atoms {\bf i+1} and {\bf i+2} are called {\em \normindex{exclusions}} of atom {\bf i}.
1964 \centerline{\includegraphics[width=8cm]{plots/chain}}
1965 \caption{Atoms along an alkane chain.}
1969 For third neighbors, the normal Lennard-Jones repulsion is sometimes
1970 still too strong, which means that when applied to a molecule, the
1971 molecule would deform or break due to the internal strain. This is
1972 especially the case for carbon-carbon interactions in a {\em
1973 cis}-conformation ({\eg} {\em cis}-butane). Therefore, for some of these
1974 interactions, the Lennard-Jones repulsion has been reduced in the
1975 {\gromos} force field, which is implemented by keeping a separate list of
1976 1-4 and normal Lennard-Jones parameters. In other force fields, such
1977 as OPLS~\cite{Jorgensen88}, the standard Lennard-Jones parameters are reduced
1978 by a factor of two, but in that case also the dispersion (r$^{-6}$)
1979 and the Coulomb interaction are scaled.
1980 {\gromacs} can use either of these methods.
1982 \subsection{Charge Groups\index{charge group}}
1984 In principle, the force calculation in MD is an $O(N^2)$ problem.
1985 Therefore, we apply a \normindex{cut-off} for non-bonded force (NBF)
1986 calculations; only the particles within a certain distance of each
1987 other are interacting. This reduces the cost to $O(N)$ (typically
1988 $100N$ to $200N$) of the NBF. It also introduces an error, which is,
1989 in most cases, acceptable, except when applying the cut-off implies
1990 the creation of charges, in which case you should consider using the
1991 lattice sum methods provided by {\gromacs}.
1993 Consider a water molecule interacting with another atom. If we would apply
1994 a plain cut-off on an atom-atom basis we might include the atom-oxygen
1995 interaction (with a charge of $-0.82$) without the compensating charge
1996 of the protons, and as a result, induce a large dipole moment over the system.
1997 Therefore, we have to keep groups of atoms with total charge
1998 0 together. These groups are called {\em charge groups}. Note that with
1999 a proper treatment of long-range electrostatics (e.g. particle-mesh Ewald
2000 (\secref{pme}), keeping charge groups together is not required.
2002 \subsection{Treatment of Cut-offs in the group scheme\index{cut-off}}
2003 \newcommand{\rs}{$r_{short}$}
2004 \newcommand{\rl}{$r_{long}$}
2005 {\gromacs} is quite flexible in treating cut-offs, which implies
2006 there can be quite a number of parameters to set. These parameters are
2007 set in the input file for {\tt grompp}. There are two sort of parameters
2008 that affect the cut-off interactions; you can select which type
2009 of interaction to use in each case, and which cut-offs should be
2010 used in the neighbor searching.
2012 For both Coulomb and van der Waals interactions there are interaction
2013 type selectors (termed {\tt vdwtype} and {\tt coulombtype}) and two
2014 parameters, for a total of six non-bonded interaction parameters. See
2015 \secref{mdpopt} for a complete description of these parameters.
2017 The neighbor searching (NS) can be performed using a single-range, or a twin-range
2018 approach. Since the former is merely a special case of the latter, we will
2019 discuss the more general twin-range. In this case, NS is described by two
2020 radii: {\tt rlist} and max({\tt rcoulomb},{\tt rvdw}).
2021 Usually one builds the neighbor list every 10 time steps
2022 or every 20 fs (parameter {\tt nstlist}). In the neighbor list, all interaction
2023 pairs that fall within {\tt rlist} are stored. Furthermore, the
2024 interactions between pairs that do not
2025 fall within {\tt rlist} but do fall within max({\tt rcoulomb},{\tt rvdw})
2026 are computed during NS. The
2027 forces and energy are stored separately and added to short-range forces
2028 at every time step between successive NS. If {\tt rlist} =
2029 max({\tt rcoulomb},{\tt rvdw}), no forces
2030 are evaluated during neighbor list generation.
2031 The \normindex{virial} is calculated from the sum of the short- and
2033 This means that the virial can be slightly asymmetrical at non-NS steps.
2034 When mdrun is compiled to use mixed precision, the virial is almost always asymmetrical because the
2035 off-diagonal elements are about as large as each element in the sum.
2036 In most cases this is not really a problem, since the fluctuations in the
2037 virial can be 2 orders of magnitude larger than the average.
2039 Except for the plain cut-off,
2040 all of the interaction functions in \tabref{funcparm}
2041 require that neighbor searching be done with a larger radius than the $r_c$
2042 specified for the functional form, because of the use of charge groups.
2043 The extra radius is typically of the order of 0.25 nm (roughly the
2044 largest distance between two atoms in a charge group plus the distance a
2045 charge group can diffuse within neighbor list updates).
2047 %If your charge groups are very large it may be interesting to turn off charge
2048 %groups, by setting the option
2049 %{\tt bAtomList = yes} in your {\tt grompp.mdp} file.
2050 %In this case only a small extra radius to account for diffusion needs to be
2051 %added (0.1 nm). Do not however use this together with the plain cut-off
2052 %method, as it will generate large artifacts (\secref{cg}).
2053 %In summary, there are four parameters that describe NS behavior:
2054 %{\tt nstlist} (update frequency in number of time steps),
2055 %{\tt bAtomList} (whether or not charge groups are used to generate neighbor list, the default is to use charge groups, so {\tt bAtomList = no}),
2056 %{\tt rshort} and {\tt rlong} which are the two radii {\rs} and {\rl}
2061 \begin{tabular}{|ll|l|}
2063 \multicolumn{2}{|c|}{Type} & Parameters \\
2065 Coulomb&Plain cut-off & $r_c$, $\epsr$ \\
2066 &Reaction field & $r_c$, $\epsrf$ \\
2067 &Shift function & $r_1$, $r_c$, $\epsr$ \\
2068 &Switch function & $r_1$, $r_c$, $\epsr$ \\
2070 VdW&Plain cut-off & $r_c$ \\
2071 &Shift function & $r_1$, $r_c$ \\
2072 &Switch function & $r_1$, $r_c$ \\
2075 \caption[Parameters for the different functional forms of the
2076 non-bonded interactions.]{Parameters for the different functional
2077 forms of the non-bonded interactions.}
2078 \label{tab:funcparm}
2080 %} % Brace matches ifthenelse test for gmxlite
2083 \newcommand{\vvis}{\ve{r}_s}
2084 \newcommand{\Fi}{\ve{F}_i'}
2085 \newcommand{\Fj}{\ve{F}_j'}
2086 \newcommand{\Fk}{\ve{F}_k'}
2087 \newcommand{\Fl}{\ve{F}_l'}
2088 \newcommand{\Fvis}{\ve{F}_{s}}
2089 \newcommand{\rvik}{\ve{r}_{ik}}
2090 \newcommand{\rvis}{\ve{r}_{is}}
2091 \newcommand{\rvjk}{\ve{r}_{jk}}
2092 \newcommand{\rvjl}{\ve{r}_{jl}}
2094 %\ifthenelse{\equal{\gmxlite}{1}}{}{
2095 \section{Virtual interaction sites\index{virtual interaction sites}}
2096 \label{sec:virtual_sites}
2097 Virtual interaction sites (called \seeindex{dummy atoms}{virtual interaction sites} in {\gromacs} versions before 3.3)
2098 can be used in {\gromacs} in a number of ways.
2099 We write the position of the virtual site $\ve{r}_s$ as a function of
2100 the positions of other particles \ve{r}$_i$: $\ve{r}_s =
2101 f(\ve{r}_1..\ve{r}_n)$. The virtual site, which may carry charge or be
2102 involved in other interactions, can now be used in the force
2103 calculation. The force acting on the virtual site must be
2104 redistributed over the particles with mass in a consistent way.
2105 A good way to do this can be found in ref.~\cite{Berendsen84b}.
2106 We can write the potential energy as:
2108 V = V(\vvis,\ve{r}_1,\ldots,\ve{r}_n) = V^*(\ve{r}_1,\ldots,\ve{r}_n)
2110 The force on the particle $i$ is then:
2112 \ve{F}_i = -\frac{\partial V^*}{\partial \ve{r}_i}
2113 = -\frac{\partial V}{\partial \ve{r}_i} -
2114 \frac{\partial V}{\partial \vvis}
2115 \frac{\partial \vvis}{\partial \ve{r}_i}
2116 = \ve{F}_i^{direct} + \Fi
2118 The first term is the normal force.
2119 The second term is the force on particle $i$ due to the virtual site, which
2120 can be written in tensor notation:
2121 \newcommand{\partd}[2]{\displaystyle\frac{\partial #1}{\partial #2_i}}
2123 \Fi = \left[\begin{array}{ccc}
2124 \partd{x_s}{x} & \partd{y_s}{x} & \partd{z_s}{x} \\[1ex]
2125 \partd{x_s}{y} & \partd{y_s}{y} & \partd{z_s}{y} \\[1ex]
2126 \partd{x_s}{z} & \partd{y_s}{z} & \partd{z_s}{z}
2127 \end{array}\right]\Fvis
2130 where $\Fvis$ is the force on the virtual site and $x_s$, $y_s$ and
2131 $z_s$ are the coordinates of the virtual site. In this way, the total
2132 force and the total torque are conserved~\cite{Berendsen84b}.
2134 The computation of the \normindex{virial}
2135 (\eqnref{Xi}) is non-trivial when virtual sites are used. Since the
2136 virial involves a summation over all the atoms (rather than virtual
2137 sites), the forces must be redistributed from the virtual sites to the
2138 atoms (using ~\eqnref{fvsite}) {\em before} computation of the
2139 virial. In some special cases where the forces on the atoms can be
2140 written as a linear combination of the forces on the virtual sites (types 2
2141 and 3 below) there is no difference between computing the virial
2142 before and after the redistribution of forces. However, in the
2143 general case redistribution should be done first.
2146 \centerline{\includegraphics[width=15cm]{plots/dummies}}
2147 \caption[Virtual site construction.]{The six different types of virtual
2148 site construction in \protect{\gromacs}. The constructing atoms are
2149 shown as black circles, the virtual sites in gray.}
2153 There are six ways to construct virtual sites from surrounding atoms in
2154 {\gromacs}, which we classify by the number of constructing
2155 atoms. {\bf Note} that all site types mentioned can be constructed from
2156 types 3fd (normalized, in-plane) and 3out (non-normalized, out of
2157 plane). However, the amount of computation involved increases sharply
2158 along this list, so we strongly recommended using the first adequate
2159 virtual site type that will be sufficient for a certain purpose.
2160 \figref{vsites} depicts 6 of the available virtual site constructions.
2161 The conceptually simplest construction types are linear combinations:
2163 \vvis = \sum_{i=1}^N w_i \, \ve{r}_i
2165 The force is then redistributed using the same weights:
2170 The types of virtual sites supported in {\gromacs} are given in the list below.
2171 Constructing atoms in virtual sites can be virtual sites themselves, but
2172 only if they are higher in the list, i.e. virtual sites can be
2173 constructed from ``particles'' that are simpler virtual sites.
2175 \item[{\bf\sf 2.}]\label{subsec:vsite2}As a linear combination of two atoms
2176 (\figref{vsites} 2):
2178 w_i = 1 - a ~,~~ w_j = a
2180 In this case the virtual site is on the line through atoms $i$ and
2183 \item[{\bf\sf 3.}]\label{subsec:vsite3}As a linear combination of three atoms
2184 (\figref{vsites} 3):
2186 w_i = 1 - a - b ~,~~ w_j = a ~,~~ w_k = b
2188 In this case the virtual site is in the plane of the other three
2191 \item[{\bf\sf 3fd.}]\label{subsec:vsite3fd}In the plane of three atoms, with a fixed distance
2192 (\figref{vsites} 3fd):
2194 \vvis ~=~ \ve{r}_i + b \frac{ \rvij + a \rvjk }
2195 {| \rvij + a \rvjk |}
2197 In this case the virtual site is in the plane of the other three
2198 particles at a distance of $|b|$ from $i$.
2199 The force on particles $i$, $j$ and $k$ due to the force on the virtual
2200 site can be computed as:
2203 \Fi &=& \displaystyle \Fvis - \gamma ( \Fvis - \ve{p} ) \\[1ex]
2204 \Fj &=& \displaystyle (1-a)\gamma (\Fvis - \ve{p}) \\[1ex]
2205 \Fk &=& \displaystyle a \gamma (\Fvis - \ve{p}) \\
2209 \displaystyle \gamma = \frac{b}{| \rvij + a \rvjk |} \\[2ex]
2210 \displaystyle \ve{p} = \frac{ \rvis \cdot \Fvis }
2211 { \rvis \cdot \rvis } \rvis
2215 \item[{\bf\sf 3fad.}]\label{subsec:vsite3fad}In the plane of three atoms, with a fixed angle and
2216 distance (\figref{vsites} 3fad):
2218 \label{eqn:vsite2fad-F}
2219 \vvis ~=~ \ve{r}_i +
2220 d \cos \theta \frac{\rvij}{|\rvij|} +
2221 d \sin \theta \frac{\ve{r}_\perp}{|\ve{r}_\perp|}
2223 \ve{r}_\perp ~=~ \rvjk -
2224 \frac{ \rvij \cdot \rvjk }
2225 { \rvij \cdot \rvij }
2228 In this case the virtual site is in the plane of the other three
2229 particles at a distance of $|d|$ from $i$ at an angle of
2230 $\alpha$ with $\rvij$. Atom $k$ defines the plane and the
2231 direction of the angle. {\bf Note} that in this case $b$ and
2232 $\alpha$ must be specified, instead of $a$ and $b$ (see also
2233 \secref{vsitetop}). The force on particles $i$, $j$ and $k$
2234 due to the force on the virtual site can be computed as (with
2235 $\ve{r}_\perp$ as defined in \eqnref{vsite2fad-F}):
2236 \newcommand{\dfrac}{\displaystyle\frac}
2239 \begin{array}{lclllll}
2241 \dfrac{d \cos \theta}{|\rvij|} \ve{F}_1 &+&
2242 \dfrac{d \sin \theta}{|\ve{r}_\perp|} \left(
2243 \dfrac{ \rvij \cdot \rvjk }
2244 { \rvij \cdot \rvij } \ve{F}_2 +
2245 \ve{F}_3 \right) \\[3ex]
2247 \dfrac{d \cos \theta}{|\rvij|} \ve{F}_1 &-&
2248 \dfrac{d \sin \theta}{|\ve{r}_\perp|} \left(
2250 \dfrac{ \rvij \cdot \rvjk }
2251 { \rvij \cdot \rvij } \ve{F}_2 +
2252 \ve{F}_3 \right) \\[3ex]
2254 \dfrac{d \sin \theta}{|\ve{r}_\perp|} \ve{F}_2 \\[3ex]
2258 \dfrac{ \rvij \cdot \Fvis }
2259 { \rvij \cdot \rvij } \rvij
2261 \ve{F}_2 = \ve{F}_1 -
2262 \dfrac{ \ve{r}_\perp \cdot \Fvis }
2263 { \ve{r}_\perp \cdot \ve{r}_\perp } \ve{r}_\perp
2265 \ve{F}_3 = \dfrac{ \rvij \cdot \Fvis }
2266 { \rvij \cdot \rvij } \ve{r}_\perp
2270 \item[{\bf\sf 3out.}]\label{subsec:vsite3out}As a non-linear combination of three atoms, out of plane
2271 (\figref{vsites} 3out):
2273 \vvis ~=~ \ve{r}_i + a \rvij + b \rvik +
2274 c (\rvij \times \rvik)
2276 This enables the construction of virtual sites out of the plane of the
2278 The force on particles $i,j$ and $k$ due to the force on the virtual
2279 site can be computed as:
2283 \Fj &=& \left[\begin{array}{ccc}
2284 a & -c\,z_{ik} & c\,y_{ik} \\[0.5ex]
2285 c\,z_{ik} & a & -c\,x_{ik} \\[0.5ex]
2286 -c\,y_{ik} & c\,x_{ik} & a
2287 \end{array}\right]\Fvis \\
2289 \Fk &=& \left[\begin{array}{ccc}
2290 b & c\,z_{ij} & -c\,y_{ij} \\[0.5ex]
2291 -c\,z_{ij} & b & c\,x_{ij} \\[0.5ex]
2292 c\,y_{ij} & -c\,x_{ij} & b
2293 \end{array}\right]\Fvis \\
2294 \Fi &=& \Fvis - \Fj - \Fk
2298 \item[{\bf\sf 4fdn.}]\label{subsec:vsite4fdn}From four atoms, with a fixed distance, see separate Fig.\ \ref{fig:vsite-4fdn}.
2299 This construction is a bit
2300 complex, in particular since the previous type (4fd) could be unstable which forced us
2301 to introduce a more elaborate construction:
2304 \centerline{\includegraphics[width=5cm]{plots/vsite-4fdn}}
2305 \caption {The new 4fdn virtual site construction, which is stable even when all constructing
2306 atoms are in the same plane.}
2307 \label{fig:vsite-4fdn}
2311 \mathbf{r}_{ja} &=& a\, \mathbf{r}_{ik} - \mathbf{r}_{ij} = a\, (\mathbf{x}_k - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\
2312 \mathbf{r}_{jb} &=& b\, \mathbf{r}_{il} - \mathbf{r}_{ij} = b\, (\mathbf{x}_l - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\
2313 \mathbf{r}_m &=& \mathbf{r}_{ja} \times \mathbf{r}_{jb} \nonumber \\
2314 \mathbf{x}_s &=& \mathbf{x}_i + c \frac{\mathbf{r}_m}{|\mathbf{r}_m|}
2318 In this case the virtual site is at a distance of $|c|$ from $i$, while $a$ and $b$ are
2319 parameters. {\bf Note} that the vectors $\mathbf{r}_{ik}$ and $\mathbf{r}_{ij}$ are not normalized
2320 to save floating-point operations.
2321 The force on particles $i$, $j$, $k$ and $l$ due to the force
2322 on the virtual site are computed through chain rule derivatives
2323 of the construction expression. This is exact and conserves energy,
2324 but it does lead to relatively lengthy expressions that we do not
2325 include here (over 200 floating-point operations). The interested reader can
2326 look at the source code in \verb+vsite.c+. Fortunately, this vsite type is normally
2327 only used for chiral centers such as $C_{\alpha}$ atoms in proteins.
2329 The new 4fdn construct is identified with a `type' value of 2 in the topology. The earlier 4fd
2330 type is still supported internally (`type' value 1), but it should not be used for
2331 new simulations. All current {\gromacs} tools will automatically generate type 4fdn instead.
2334 \item[{\bf\sf N.}]\label{subsec:vsiteN} A linear combination of $N$ atoms with relative
2335 weights $a_i$. The weight for atom $i$ is:
2337 w_i = a_i \left(\sum_{j=1}^N a_j \right)^{-1}
2339 There are three options for setting the weights:
2341 \item[COG] center of geometry: equal weights
2342 \item[COM] center of mass: $a_i$ is the mass of atom $i$;
2343 when in free-energy simulations the mass of the atom is changed,
2344 only the mass of the A-state is used for the weight
2345 \item[COW] center of weights: $a_i$ is defined by the user
2349 %} % Brace matches ifthenelse test for gmxlite
2351 \newcommand{\dr}{{\rm d}r}
2352 \newcommand{\avcsix}{\left< C_6 \right>}
2354 %\ifthenelse{\equal{\gmxlite}{1}}{}{
2355 \section{Long Range Electrostatics}
2356 \label{sec:lr_elstat}
2357 \subsection{Ewald summation\index{Ewald sum}}
2359 The total electrostatic energy of $N$ particles and their periodic
2360 images\index{periodic boundary conditions} is given by
2362 V=\frac{f}{2}\sum_{n_x}\sum_{n_y}
2363 \sum_{n_{z}*} \sum_{i}^{N} \sum_{j}^{N}
2364 \frac{q_i q_j}{{\bf r}_{ij,{\bf n}}}.
2365 \label{eqn:totalcoulomb}
2367 $(n_x,n_y,n_z)={\bf n}$ is the box index vector, and the star indicates that
2368 terms with $i=j$ should be omitted when $(n_x,n_y,n_z)=(0,0,0)$. The
2369 distance ${\bf r}_{ij,{\bf n}}$ is the real distance between the charges and
2370 not the minimum-image. This sum is conditionally convergent, but
2373 Ewald summation was first introduced as a method to calculate
2374 long-range interactions of the periodic images in
2375 crystals~\cite{Ewald21}. The idea is to convert the single
2376 slowly-converging sum \eqnref{totalcoulomb} into two
2377 quickly-converging terms and a constant term:
2379 V &=& V_{\mathrm{dir}} + V_{\mathrm{rec}} + V_{0} \\[0.5ex]
2380 V_{\mathrm{dir}} &=& \frac{f}{2} \sum_{i,j}^{N}
2381 \sum_{n_x}\sum_{n_y}
2382 \sum_{n_{z}*} q_i q_j \frac{\mbox{erfc}(\beta {r}_{ij,{\bf n}} )}{{r}_{ij,{\bf n}}} \\[0.5ex]
2383 V_{\mathrm{rec}} &=& \frac{f}{2 \pi V} \sum_{i,j}^{N} q_i q_j
2384 \sum_{m_x}\sum_{m_y}
2385 \sum_{m_{z}*} \frac{\exp{\left( -(\pi {\bf m}/\beta)^2 + 2 \pi i
2386 {\bf m} \cdot ({\bf r}_i - {\bf r}_j)\right)}}{{\bf m}^2} \\[0.5ex]
2387 V_{0} &=& -\frac{f \beta}{\sqrt{\pi}}\sum_{i}^{N} q_i^2,
2389 where $\beta$ is a parameter that determines the relative weight of the
2390 direct and reciprocal sums and ${\bf m}=(m_x,m_y,m_z)$.
2391 In this way we can use a short cut-off (of the order of $1$~nm) in the direct space sum and a
2392 short cut-off in the reciprocal space sum ({\eg} 10 wave vectors in each
2393 direction). Unfortunately, the computational cost of the reciprocal
2394 part of the sum increases as $N^2$
2395 (or $N^{3/2}$ with a slightly better algorithm) and it is therefore not
2396 realistic for use in large systems.
2398 \subsubsection{Using Ewald}
2399 Don't use Ewald unless you are absolutely sure this is what you want -
2400 for almost all cases the PME method below will perform much better.
2401 If you still want to employ classical Ewald summation enter this in
2402 your {\tt .mdp} file, if the side of your box is about $3$~nm:
2409 fourierspacing = 0.6
2413 The ratio of the box dimensions and the {\tt fourierspacing} parameter determines
2414 the highest magnitude of wave vectors $m_x,m_y,m_z$ to use in each
2415 direction. With a 3-nm cubic box this example would use $11$ wave vectors
2416 (from $-5$ to $5$) in each direction. The {\tt ewald-rtol} parameter
2417 is the relative strength of the electrostatic interaction at the
2418 cut-off. Decreasing this gives you a more accurate direct sum, but a
2419 less accurate reciprocal sum.
2421 \subsection{\normindex{PME}}
2423 Particle-mesh Ewald is a method proposed by Tom
2424 Darden~\cite{Darden93} to improve the performance of the
2425 reciprocal sum. Instead of directly summing wave vectors, the charges
2426 are assigned to a grid using interpolation. The implementation in
2427 {\gromacs} uses cardinal B-spline interpolation~\cite{Essmann95},
2428 which is referred to as smooth PME (SPME).
2429 The grid is then Fourier transformed with a 3D FFT algorithm and the
2430 reciprocal energy term obtained by a single sum over the grid in
2433 The potential at the grid points is calculated by inverse
2434 transformation, and by using the interpolation factors we get the
2435 forces on each atom.
2437 The PME algorithm scales as $N \log(N)$, and is substantially faster
2438 than ordinary Ewald summation on medium to large systems. On very
2439 small systems it might still be better to use Ewald to avoid the
2440 overhead in setting up grids and transforms.
2441 For the parallelization of PME see the section on MPMD PME (\ssecref{mpmd_pme}).
2443 With the Verlet cut-off scheme, the PME direct space potential is
2444 shifted by a constant such that the potential is zero at the
2445 cut-off. This shift is small and since the net system charge is close
2446 to zero, the total shift is very small, unlike in the case of the
2447 Lennard-Jones potential where all shifts add up. We apply the shift
2448 anyhow, such that the potential is the exact integral of the force.
2450 \subsubsection{Using PME}
2451 As an example for using Particle-mesh Ewald summation in {\gromacs}, specify the
2452 following lines in your {\tt .mdp} file:
2459 fourierspacing = 0.12
2464 In this case the {\tt fourierspacing} parameter determines the maximum
2465 spacing for the FFT grid (i.e. minimum number of grid points),
2466 and {\tt pme-order} controls the
2467 interpolation order. Using fourth-order (cubic) interpolation and this
2468 spacing should give electrostatic energies accurate to about
2469 $5\cdot10^{-3}$. Since the Lennard-Jones energies are not this
2470 accurate it might even be possible to increase this spacing slightly.
2472 Pressure scaling works with PME, but be aware of the fact that
2473 anisotropic scaling can introduce artificial ordering in some systems.
2475 \subsection{\normindex{P3M-AD}}
2477 The \seeindex{Particle-Particle Particle-Mesh}{P3M} methods of
2478 Hockney \& Eastwood can also be applied in {\gromacs} for the
2479 treatment of long range electrostatic interactions~\cite{Hockney81}.
2480 Although the P3M method was the first efficient long-range electrostatics
2481 method for molecular simulation, the smooth PME (SPME) method has largely
2482 replaced P3M as the method of choice in atomistic simulations. One performance
2483 disadvantage of the original P3M method was that it required 3 3D-FFT
2484 back transforms to obtain the forces on the particles. But this is not
2485 required for P3M and the forces can be derived through analytical differentiation
2486 of the potential, as done in PME. The resulting method is termed P3M-AD.
2487 The only remaining difference between P3M-AD and PME is the optimization
2488 of the lattice Green influence function for error minimization that P3M uses.
2489 However, in 2012 it has been shown that the SPME influence function can be
2490 modified to obtain P3M~\cite{Ballenegger2012}.
2491 This means that the advantage of error minimization in P3M-AD can be used
2492 at the same computational cost and with the same code as PME,
2493 just by adding a few lines to modify the influence function.
2494 However, at optimal parameter setting the effect of error minimization
2495 in P3M-AD is less than 10\%. P3M-AD does show large accuracy gains with
2496 interlaced (also known as staggered) grids, but that is not supported
2497 in {\gromacs} (yet).
2499 P3M is used in {\gromacs} with exactly the same options as used with PME
2500 by selecting the electrostatics type:
2502 coulombtype = P3M-AD
2505 \subsection{Optimizing Fourier transforms and PME calculations}
2506 It is recommended to optimize the parameters for calculation of
2507 electrostatic interaction such as PME grid dimensions and cut-off radii.
2508 This is particularly relevant to do before launching long production runs.
2510 {\gromacs} includes a special tool, {\tt g_tune_pme}, which automates the
2511 process of selecting the optimal size of the grid and number of PME-only
2515 % Temporarily removed since I am not sure about the state of the testlr
2518 %It is possible to test the accuracy of your settings using the program
2519 %{\tt\normindex{testlr}} in the {\tt src/gmxlib} dir. This program computes
2520 %forces and potentials using PPPM and an Ewald implementation and gives the
2521 %absolute and RMS errors in both:
2526 %Potential 0.113 0.035
2528 %{\bf Note:} these numbers were generated using a grid spacing of
2529 %0.058 nm and $r_c$ = 1.0 nm.
2531 %You can see what the accuracy is without optimizing the
2532 %$\hat{G}(k)$ function, if you pass the {\tt -ghat} option to {\tt
2533 %testlr}. Try it if you think the {\tt mk_ghat} procedure is a waste
2535 %} % Brace matches ifthenelse test for gmxlite
2538 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2539 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2540 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2542 %\ifthenelse{\equal{\gmxlite}{1}}{}{
2543 \section{Long Range Van der Waals interactions}
2544 \subsection{Dispersion correction\index{dispersion correction}}
2545 In this section, we derive long-range corrections due to the use of a
2546 cut-off for Lennard-Jones or Buckingham interactions.
2547 We assume that the cut-off is
2548 so long that the repulsion term can safely be neglected, and therefore
2549 only the dispersion term is taken into account. Due to the nature of
2550 the dispersion interaction (we are truncating a potential proportional
2551 to $-r^{-6}$), energy and pressure corrections are both negative. While
2552 the energy correction is usually small, it may be important for free
2553 energy calculations where differences between two different Hamiltonians
2554 are considered. In contrast, the pressure correction is very large and
2555 can not be neglected under any circumstances where a correct pressure is
2556 required, especially for any NPT simulations. Although it is, in
2557 principle, possible to parameterize a force field such that the pressure
2558 is close to the desired experimental value without correction, such a
2559 method makes the parameterization dependent on the cut-off and is therefore
2562 \subsubsection{Energy}
2564 The long-range contribution of the dispersion interaction to the
2565 virial can be derived analytically, if we assume a homogeneous
2566 system beyond the cut-off distance $r_c$. The dispersion energy
2567 between two particles is written as:
2569 V(\rij) ~=~- C_6\,\rij^{-6}
2571 and the corresponding force is:
2573 \Fvij ~=~- 6\,C_6\,\rij^{-8}\rvij
2575 In a periodic system it is not easy to calculate the full potentials,
2576 so usually a cut-off is applied, which can be abrupt or smooth.
2577 We will call the potential and force with cut-off $V_c$ and $\ve{F}_c$.
2578 The long-range contribution to the dispersion energy
2579 in a system with $N$ particles and particle density $\rho$ = $N/V$ is:
2581 \label{eqn:enercorr}
2582 V_{lr} ~=~ \half N \rho\int_0^{\infty} 4\pi r^2 g(r) \left( V(r) -V_c(r) \right) {\dr}
2584 We will integrate this for the shift function, which is the most general
2585 form of van der Waals interaction available in {\gromacs}.
2586 The shift function has a constant difference $S$ from 0 to $r_1$
2587 and is 0 beyond the cut-off distance $r_c$.
2588 We can integrate \eqnref{enercorr}, assuming that the density in the sphere
2589 within $r_1$ is equal to the global density and
2590 the radial distribution function $g(r)$ is 1 beyond $r_1$:
2593 V_{lr} &=& \half N \left(
2594 \rho\int_0^{r_1} 4\pi r^2 g(r) \, C_6 \,S\,{\dr}
2595 + \rho\int_{r_1}^{r_c} 4\pi r^2 \left( V(r) -V_c(r) \right) {\dr}
2596 + \rho\int_{r_c}^{\infty} 4\pi r^2 V(r) \, {\dr}
2598 & = & \half N \left(\left(\frac{4}{3}\pi \rho r_1^{3} - 1\right) C_6 \,S
2599 + \rho\int_{r_1}^{r_c} 4\pi r^2 \left( V(r) -V_c(r) \right) {\dr}
2600 -\frac{4}{3} \pi N \rho\, C_6\,r_c^{-3}
2603 where the term $-1$ corrects for the self-interaction.
2604 For a plain cut-off we only need to assume that $g(r)$ is 1 beyond $r_c$
2605 and the correction reduces to~\cite{Allen87}:
2607 V_{lr} & = & -\frac{2}{3} \pi N \rho\, C_6\,r_c^{-3}
2609 If we consider, for example, a box of pure water, simulated with a cut-off
2610 of 0.9 nm and a density of 1 g cm$^{-3}$ this correction is
2611 $-0.75$ kJ mol$^{-1}$ per molecule.
2613 For a homogeneous mixture we need to define
2614 an {\em average dispersion constant}:
2617 \avcsix = \frac{2}{N(N-1)}\sum_i^N\sum_{j>i}^N C_6(i,j)\\
2619 In {\gromacs}, excluded pairs of atoms do not contribute to the average.
2621 In the case of inhomogeneous simulation systems, {\eg} a system with a
2622 lipid interface, the energy correction can be applied if
2623 $\avcsix$ for both components is comparable.
2625 \subsubsection{Virial and pressure}
2626 The scalar virial of the system due to the dispersion interaction between
2627 two particles $i$ and $j$ is given by:
2629 \Xi~=~-\half \rvij \cdot \Fvij ~=~ 3\,C_6\,\rij^{-6}
2631 The pressure is given by:
2633 P~=~\frac{2}{3\,V}\left(E_{kin} - \Xi\right)
2635 The long-range correction to the virial is given by:
2637 \Xi_{lr} ~=~ \half N \rho \int_0^{\infty} 4\pi r^2 g(r) (\Xi -\Xi_c) \,\dr
2639 We can again integrate the long-range contribution to the
2640 virial assuming $g(r)$ is 1 beyond $r_1$:
2642 \Xi_{lr}&=& \half N \rho \left(
2643 \int_{r_1}^{r_c} 4 \pi r^2 (\Xi -\Xi_c) \,\dr
2644 + \int_{r_c}^{\infty} 4 \pi r^2 3\,C_6\,\rij^{-6}\, \dr
2646 &=& \half N \rho \left(
2647 \int_{r_1}^{r_c} 4 \pi r^2 (\Xi -\Xi_c) \, \dr
2648 + 4 \pi C_6 \, r_c^{-3} \right)
2650 For a plain cut-off the correction to the pressure is~\cite{Allen87}:
2652 P_{lr}~=~-\frac{4}{3} \pi C_6\, \rho^2 r_c^{-3}
2654 Using the same example of a water box, the correction to the virial is
2655 0.75 kJ mol$^{-1}$ per molecule,
2656 the corresponding correction to the pressure for
2657 SPC water is approximately $-280$ bar.
2659 For homogeneous mixtures, we can again use the average dispersion constant
2660 $\avcsix$ (\eqnref{avcsix}):
2662 P_{lr}~=~-\frac{4}{3} \pi \avcsix \rho^2 r_c^{-3}
2665 For inhomogeneous systems, \eqnref{pcorr} can be applied under the same
2666 restriction as holds for the energy (see \secref{ecorr}).
2668 \subsection{Lennard-Jones PME\index{LJ-PME}}
2670 In order to treat systems, using Lennard-Jones potentials, that are
2671 non-homogeneous outside of the cut-off distance, we can instead use
2672 the Particle-mesh Ewald method as discussed for electrostatics above.
2673 In this case the modified Ewald equations become
2675 V &=& V_{\mathrm{dir}} + V_{\mathrm{rec}} + V_{0} \\[0.5ex]
2676 V_{\mathrm{dir}} &=& -\frac{1}{2} \sum_{i,j}^{N}
2677 \sum_{n_x}\sum_{n_y}
2678 \sum_{n_{z}*} \frac{C^{ij}_6 g(\beta {r}_{ij,{\bf n}})}{{r_{ij,{\bf n}}}^6}
2679 \label{eqn:ljpmerealspace}\\[0.5ex]
2680 V_{\mathrm{rec}} &=& \frac{{\pi}^{\frac{3}{2}} \beta^{3}}{2V} \sum_{m_x}\sum_{m_y}\sum_{m_{z}*}
2681 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]
2682 V_{0} &=& -\frac{\beta^{6}}{12}\sum_{i}^{N} C^{ii}_6
2685 where ${\bf m}=(m_x,m_y,m_z)$, $\beta$ is the parameter determining the weight between
2686 direct and reciprocal space, and ${C^{ij}_6}$ is the combined dispersion
2687 parameter for particle $i$ and $j$. The star indicates that terms
2688 with $i = j$ should be omitted when $((n_x,n_y,n_z)=(0,0,0))$, and
2689 ${\bf r}_{ij,{\bf n}}$ is the real distance between the particles.
2690 Following the derivation by Essmann~\cite{Essmann95}, the functions $f$ and $g$ introduced above are defined as
2692 f(x)&=&1/3\left[(1-2x^2){\mathrm{exp}}(-x^2) + 2{x^3}\sqrt{\pi}\,{\mathrm{erfc}}(x) \right] \\
2693 g(x)&=&{\mathrm{exp}}(-x^2)(1+x^2+\frac{x^4}{2}).
2696 The above methodology works fine as long as the dispersion parameters can be combined geometrically (\eqnref{comb}) in the same
2697 way as the charges for electrostatics
2699 C^{ij}_{6,\mathrm{geom}} = \left(C^{ii}_6 \, C^{jj}_6\right)^{1/2}
2701 For Lorentz-Berthelot combination rules (\eqnref{lorentzberthelot}), the reciprocal part of this sum has to be calculated
2702 seven times due to the splitting of the dispersion parameter according to
2704 C^{ij}_{6,\mathrm{L-B}} = (\sigma_i+\sigma_j)^6=\sum_{n=0}^{6} P_{n}\sigma_{i}^{n}\sigma_{j}^{(6-n)},
2706 for $P_{n}$ the Pascal triangle coefficients. This introduces a
2707 non-negligible cost to the reciprocal part, requiring seven separate
2708 FFTs, and therefore this has been the limiting factor in previous
2709 attempts to implement LJ-PME. A solution to this problem is to use
2710 geometrical combination rules in order to calculate an approximate
2711 interaction parameter for the reciprocal part of the potential,
2712 yielding a total interaction of
2714 V(r<r_c) & = & \underbrace{C^{\mathrm{dir}}_6 g(\beta r) r^{-6}}_{\mathrm{Direct \ space}} + \underbrace{C^\mathrm{recip}_{6,\mathrm{geom}} [1 - g(\beta r)] r^{-6}}_{\mathrm{Reciprocal \ space}} \nonumber \\
2715 &=& C^\mathrm{recip}_{6,\mathrm{geom}}r^{-6} + \left(C^{\mathrm{dir}}_6-C^\mathrm{recip}_{6,\mathrm{geom}}\right)g(\beta r)r^{-6} \\
2716 V(r>r_c) & = & \underbrace{C^\mathrm{recip}_{6,\mathrm{geom}} [1 - g(\beta r)] r^{-6}}_{\mathrm{Reciprocal \ space}}.
2718 This will preserve a well-defined Hamiltonian and significantly increase
2719 the performance of the simulations. The approximation does introduce
2720 some errors, but since the difference is located in the interactions
2721 calculated in reciprocal space, the effect will be very small compared
2722 to the total interaction energy. In a simulation of a lipid bilayer,
2723 using a cut-off of 1.0 nm, the relative error in total dispersion
2724 energy was below 0.5\%. A more thorough discussion of this can be
2725 found in \cite{Wennberg13}.
2727 In {\gromacs} we now perform the proper calculation of this interaction
2728 by subtracting, from the direct-space interactions, the contribution
2729 made by the approximate potential that is used in the reciprocal part
2731 V_\mathrm{dir} = C^{\mathrm{dir}}_6 r^{-6} - C^\mathrm{recip}_6 [1 - g(\beta r)] r^{-6}.
2732 \label{eqn:ljpmedirectspace}
2734 This potential will reduce to the expression in \eqnref{ljpmerealspace} when $C^{\mathrm{dir}}_6 = C^\mathrm{recip}_6$,
2735 and the total interaction is given by
2737 \nonumber V(r<r_c) &=& \underbrace{C^{\mathrm{dir}}_6 r^{-6} - C^\mathrm{recip}_6 [1 - g(\beta r)] r^{-6}}_{\mathrm{Direct \ space}} + \underbrace{C^\mathrm{recip}_6 [1 - g(\beta r)] r^{-6}}_{\mathrm{Reciprocal \ space}} \\
2738 &=&C^{\mathrm{dir}}_6 r^{-6}
2739 \label {eqn:ljpmecorr2} \\
2740 V(r>r_c) &=& C^\mathrm{recip}_6 [1 - g(\beta r)] r^{-6}.
2742 For the case when $C^{\mathrm{dir}}_6 \neq C^\mathrm{recip}_6$ this
2743 will retain an unmodified LJ force up to the cut-off, and the error
2744 is an order of magnitude smaller than in simulations where the
2745 direct-space interactions do not account for the approximation used in
2746 reciprocal space. When using a VdW interaction modifier of
2747 potential-shift, the constant
2749 \left(-C^{\mathrm{dir}}_6 + C^\mathrm{recip}_6 [1 - g(\beta r_c)]\right) r_c^{-6}
2751 is added to \eqnref{ljpmecorr2} in order to ensure that the potential
2752 is continuous at the cutoff. Note that, in the same way as \eqnref{ljpmedirectspace}, this degenerates into the
2753 expected $-C_6g(\beta r_c)r^{-6}_c$ when $C^{\mathrm{dir}}_6 =
2754 C^\mathrm{recip}_6$. In addition to this, a long-range dispersion
2755 correction can be applied to correct for the approximation using a
2756 combination rule in reciprocal space. This correction assumes, as for
2757 the cut-off LJ potential, a uniform particle distribution. But since
2758 the error of the combination rule approximation is very small this
2759 long-range correction is not necessary in most cases. Also note that
2760 this homogenous correction does not correct the surface tension, which
2761 is an inhomogeneous property.
2763 \subsubsection{Using LJ-PME}
2764 As an example for using Particle-mesh Ewald summation for Lennard-Jones interactions in {\gromacs}, specify the
2765 following lines in your {\tt .mdp} file:
2769 vdw-modifier = Potential-Shift
2772 fourierspacing = 0.12
2774 ewald-rtol-lj = 0.001
2775 lj-pme-comb-rule = geometric
2778 The same Fourier grid and interpolation order are used if both
2779 LJ-PME and electrostatic PME are active, so the settings for
2780 {\tt fourierspacing} and {\tt pme-order} are common to both.
2781 {\tt ewald-rtol-lj} controls the
2782 splitting between direct and reciprocal space in the same way as
2783 {\tt ewald-rtol}. In addition to this, the combination rule to be used
2784 in reciprocal space is determined by {\tt lj-pme-comb-rule}. If the
2785 current force field uses Lorentz-Berthelot combination rules, it is
2786 possible to set {\tt lj-pme-comb-rule = geometric} in order to gain a
2787 significant increase in performance for a small loss in accuracy. The
2788 details of this approximation can be found in the section above.
2790 Note that the use of a complete long-range dispersion correction means
2791 that as with Coulomb PME, {\tt rvdw} is now a free parameter in the
2792 method, rather than being necessarily restricted by the force-field
2793 parameterization scheme. Thus it is now possible to optimize the
2794 cutoff, spacing, order and tolerance terms for accuracy and best
2797 Naturally, the use of LJ-PME rather than LJ cut-off adds computation
2798 and communication done for the reciprocal-space part, so for best
2799 performance in balancing the load of parallel simulations using
2800 PME-only ranks, more such ranks should be used. It may be possible to
2801 improve upon the automatic load-balancing used by {\tt mdrun}.
2803 %} % Brace matches ifthenelse test for gmxlite
2805 \section{Force field\index{force field}}
2807 A force field is built up from two distinct components:
2809 \item The set of equations (called the {\em
2810 \index{potential function}s}) used to generate the potential
2811 energies and their derivatives, the forces. These are described in
2812 detail in the previous chapter.
2813 \item The parameters used in this set of equations. These are not
2814 given in this manual, but in the data files corresponding to your
2815 {\gromacs} distribution.
2817 Within one set of equations various sets of parameters can be
2818 used. Care must be taken that the combination of equations and
2819 parameters form a consistent set. It is in general dangerous to make
2820 {\em ad hoc} changes in a subset of parameters, because the various
2821 contributions to the total force are usually interdependent. This
2822 means in principle that every change should be documented, verified by
2823 comparison to experimental data and published in a peer-reviewed
2824 journal before it can be used.
2826 {\gromacs} {\gmxver} includes several force fields, and additional
2827 ones are available on the website. If you do not know which one to
2828 select we recommend \gromosv{96} for united-atom setups and OPLS-AA/L for
2829 all-atom parameters. That said, we describe the available options in
2832 \subsubsection{All-hydrogen force field}
2833 The \gromosv{87}-based all-hydrogen force field is almost identical to the
2834 normal \gromosv{87} force field, since the extra hydrogens have no
2835 Lennard-Jones interaction and zero charge. The only differences are in
2836 the bond angle and improper dihedral angle terms. This force field is
2837 only useful when you need the exact hydrogen positions, for instance
2838 for distance restraints derived from NMR measurements. When citing
2839 this force field please read the previous paragraph.
2841 \subsection{\gromosv{96}\index{GROMOS96 force field}}
2842 {\gromacs} supports the \gromosv{96} force fields~\cite{gromos96}.
2843 All parameters for the 43A1, 43A2 (development, improved alkane
2844 dihedrals), 45A3, 53A5, and 53A6 parameter sets are included. All standard
2845 building blocks are included and topologies can be built automatically
2848 The \gromosv{96} force field is a further development of the \gromosv{87} force field.
2849 It has improvements over the \gromosv{87} force field for proteins and small molecules.
2850 {\bf Note} that the sugar parameters present in 53A6 do correspond to those published in
2851 2004\cite{Oostenbrink2004}, which are different from those present in 45A4, which
2852 is not included in {\gromacs} at this time. The 45A4 parameter set corresponds to a later
2853 revision of these parameters.
2854 The \gromosv{96} force field is not, however, recommended for use with long alkanes and
2855 lipids. The \gromosv{96} force field differs from the \gromosv{87}
2856 force field in a few respects:
2858 \item the force field parameters
2859 \item the parameters for the bonded interactions are not linked to atom types
2860 \item a fourth power bond stretching potential (\ssecref{G96bond})
2861 \item an angle potential based on the cosine of the angle (\ssecref{G96angle})
2863 There are two differences in implementation between {\gromacs} and \gromosv{96}
2864 which can lead to slightly different results when simulating the same system
2867 \item in \gromosv{96} neighbor searching for solvents is performed on the
2868 first atom of the solvent molecule. This is not implemented in {\gromacs},
2869 but the difference with searching by centers of charge groups is very small
2870 \item the virial in \gromosv{96} is molecule-based. This is not implemented in
2871 {\gromacs}, which uses atomic virials
2873 The \gromosv{96} force field was parameterized with a Lennard-Jones cut-off
2874 of 1.4 nm, so be sure to use a Lennard-Jones cut-off ({\tt rvdw}) of at least 1.4.
2875 A larger cut-off is possible because the Lennard-Jones potential and forces
2876 are almost zero beyond 1.4 nm.
2878 \subsubsection{\gromosv{96} files\swapindexquiet{GROMOS96}{files}}
2879 {\gromacs} can read and write \gromosv{96} coordinate and trajectory files.
2880 These files should have the extension {\tt .g96}.
2881 Such a file can be a \gromosv{96} initial/final
2882 configuration file, a coordinate trajectory file, or a combination of both.
2883 The file is fixed format; all floats are written as 15.9, and as such, files can get huge.
2884 {\gromacs} supports the following data blocks in the given order:
2894 POSITION/POSITIONRED (mandatory)
2895 VELOCITY/VELOCITYRED (optional)
2900 See the \gromosv{96} manual~\cite{gromos96} for a complete description
2901 of the blocks. {\bf Note} that all {\gromacs} programs can read compressed
2902 (.Z) or gzipped (.gz) files.
2904 \subsection{OPLS/AA\index{OPLS/AA force field}}
2906 \subsection{AMBER\index{AMBER force field}}
2908 {\gromacs} provides native support for the following AMBER force fields:
2911 \item AMBER94~\cite{Cornell1995}
2912 \item AMBER96~\cite{Kollman1996}
2913 \item AMBER99~\cite{Wang2000}
2914 \item AMBER99SB~\cite{Hornak2006}
2915 \item AMBER99SB-ILDN~\cite{Lindorff2010}
2916 \item AMBER03~\cite{Duan2003}
2917 \item AMBERGS~\cite{Garcia2002}
2920 \subsection{CHARMM\index{CHARMM force field}}
2921 \label{subsec:charmmff}
2923 {\gromacs} supports the CHARMM force field for proteins~\cite{mackerell04, mackerell98}, lipids~\cite{feller00} and nucleic acids~\cite{foloppe00,Mac2000}. 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.
2925 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.
2927 A port of the CHARMM36 force field for use with GROMACS is also available at \url{http://mackerell.umaryland.edu/charmm_ff.shtml#gromacs}.
2929 \subsection{Coarse-grained force fields}
2930 \index{force-field, coarse-grained}
2931 \label{sec:cg-forcefields}
2932 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.
2934 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:
2936 \item Conserving free energies
2938 \item Simplex method
2939 \item MARTINI force field (see next section)
2941 \item Conserving distributions (like the radial distribution function), so-called structure-based coarse-graining
2943 \item (iterative) Boltzmann inversion
2944 \item Inverse Monte Carlo
2946 \item Conversing forces
2948 \item Force matching
2952 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}.
2954 \subsection{MARTINI\index{Martini force field}}
2956 The MARTINI force field is a coarse-grain parameter set that allows for the construction
2957 of many systems, including proteins and membranes.
2959 \subsection{PLUM\index{PLUM force field}}
2961 The \normindex{PLUM force field}~\cite{bereau12} is an example of a solvent-free
2962 protein-membrane model for which the membrane was derived from structure-based
2963 coarse-graining~\cite{wang_jpcb10}. A {\gromacs} implementation can be found at
2964 \href{http://code.google.com/p/plumx/}{code.google.com/p/plumx}.
2966 % LocalWords: dihedrals centro ij dV dr LJ lj rcl jj Bertelot OPLS bh bham rf
2967 % LocalWords: coul defunits grompp crf vcrf fcrf Tironi Debye kgrf cgrf krf dx
2968 % LocalWords: PPPM der Waals erfc lr elstat chirality bstretch bondpot kT kJ
2969 % LocalWords: anharmonic morse mol betaij expminx SPC timestep fs FENE ijk kj
2970 % LocalWords: anglepot CHARMm UB ik rr substituents ijkl Ryckaert Bellemans rb
2971 % LocalWords: alkanes pdb gmx IUPAC IUB jkl cis rbdih crb kcal cubicspline xvg
2972 % LocalWords: topfile mdrun posres ar dihr lcllll NMR nmr lcllllll NOEs lclll
2973 % LocalWords: rav preprocessor ccccccccc ai aj fac disre mdp multi topol tpr
2974 % LocalWords: fc ravdisre nstdisreout dipolar lll ccc orientational MSD const
2975 % LocalWords: orire fitgrp nstorireout Drude intra Noskov et al fecalc coulrf
2976 % LocalWords: polarizabilities parameterized sigeps Ek sc softcore GROMOS NBF
2977 % LocalWords: hydrogens alkane vdwtype coulombtype mdpopt rlist rcoulomb rvdw
2978 % LocalWords: nstlist virial funcparm VdW jk jl fvsite fd vsites lcr vsitetop
2979 % LocalWords: vsite lclllll lcl parameterize parameterization enercorr avcsix
2980 % LocalWords: pcorr ecorr totalcoulomb dir fourierspacing ewald rtol Darden gz
2981 % LocalWords: FFT parallelization MPMD mpmd pme fft hoc Gromos gromos oxygens
2982 % LocalWords: virials POSITIONRED VELOCITYRED gzipped Charmm Larsson Bjelkmar
2983 % LocalWords: Cuendet CMAP nocmap dihedral Lennard covalent Verlet
2984 % LocalWords: Berthelot nm flexwat ferguson itp harmonicangle versa
2985 % LocalWords: harmonicbond atomtypes dihedraltypes equilibrated fdn
2986 % LocalWords: distancerestraint LINCS Coulombic ja jb il SPME ILDN
2987 % LocalWords: Hamiltonians atomtype AMBERGS rtp cmap graining VOTCA