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 The constraints are formally part of the Hamiltonian, and therefore
1814 they give a contribution to the free energy. In {\gromacs} this can be
1815 calculated using the \normindex{LINCS} or the \normindex{SHAKE} algorithm.
1816 If we have $k = 1 \ldots K$ constraint equations $g_k$ for LINCS, then
1818 g_k = |\ve{r}_{k}| - d_{k}
1820 where $\ve{r}_k$ is the displacement vector between two particles and
1821 $d_k$ is the constraint distance between the two particles. We can express
1822 the fact that the constraint distance has a $\LAM$ dependency by
1824 d_k = \LL d_{k}^A + \LAM d_k^B
1827 Thus the $\LAM$-dependent constraint equation is
1829 g_k = |\ve{r}_{k}| - \left(\LL d_{k}^A + \LAM d_k^B\right).
1832 The (zero) contribution $G$ to the Hamiltonian from the constraints
1833 (using Lagrange multipliers $\lambda_k$, which are logically distinct
1834 from the free-energy $\LAM$) is
1836 G &=& \sum^K_k \lambda_k g_k \\
1837 \dvdl{G} &=& \frac{\partial G}{\partial d_k} \dvdl{d_k} \\
1838 &=& - \sum^K_k \lambda_k \left(d_k^B-d_k^A\right)
1841 For SHAKE, the constraint equations are
1843 g_k = \ve{r}_{k}^2 - d_{k}^2
1845 with $d_k$ as before, so
1847 \dvdl{G} &=& -2 \sum^K_k \lambda_k \left(d_k^B-d_k^A\right)
1850 \subsection{Soft-core interactions\index{soft-core interactions}}
1852 \centerline{\includegraphics[height=6cm]{plots/softcore}}
1853 \caption{Soft-core interactions at $\LAM=0.5$, with $p=2$ and
1854 $C_6^A=C_{12}^A=C_6^B=C_{12}^B=1$.}
1855 \label{fig:softcore}
1857 In a free-energy calculation where particles grow out of nothing, or
1858 particles disappear, using the the simple linear interpolation of the
1859 Lennard-Jones and Coulomb potentials as described in Equations~\ref{eq:dVljdlambda}
1860 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
1861 close to each other, leading to large fluctuations in the measured values of
1862 $\partial V/\partial \LAM$ (which, because of the simple linear
1863 interpolation, depends on the potentials at both the endpoints of $\LAM$).
1865 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
1866 involved at $\LAM$ values between 0 and 1, but not \emph{at} $\LAM=0$
1869 In {\gromacs} the soft-core potentials $V_{sc}$ are shifted versions of the
1870 regular potentials, so that the singularity in the potential and its
1871 derivatives at $r=0$ is never reached:
1873 V_{sc}(r) &=& \LL V^A(r_A) + \LAM V^B(r_B)
1875 r_A &=& \left(\alpha \sigma_A^6 \LAM^p + r^6 \right)^\frac{1}{6}
1877 r_B &=& \left(\alpha \sigma_B^6 \LL^p + r^6 \right)^\frac{1}{6}
1879 where $V^A$ and $V^B$ are the normal ``hard core'' Van der Waals or
1880 electrostatic potentials in state A ($\LAM=0$) and state B ($\LAM=1$)
1881 respectively, $\alpha$ is the soft-core parameter (set with {\tt sc_alpha}
1882 in the {\tt .mdp} file), $p$ is the soft-core $\LAM$ power (set with
1883 {\tt sc_power}), $\sigma$ is the radius of the interaction, which is
1884 $(C_{12}/C_6)^{1/6}$ or an input parameter ({\tt sc_sigma}) when $C_6$
1885 or $C_{12}$ is zero.
1887 For intermediate $\LAM$, $r_A$ and $r_B$ alter the interactions very little
1888 for $r > \alpha^{1/6} \sigma$ and quickly switch the soft-core
1889 interaction to an almost constant value for smaller $r$ (\figref{softcore}).
1892 F_{sc}(r) = -\frac{\partial V_{sc}(r)}{\partial r} =
1893 \LL F^A(r_A) \left(\frac{r}{r_A}\right)^5 +
1894 \LAM F^B(r_B) \left(\frac{r}{r_B}\right)^5
1896 where $F^A$ and $F^B$ are the ``hard core'' forces.
1897 The contribution to the derivative of the free energy is:
1899 \dvdl{V_{sc}(r)} & = &
1900 V^B(r_B) -V^A(r_A) +
1901 \LL \frac{\partial V^A(r_A)}{\partial r_A}
1902 \frac{\partial r_A}{\partial \LAM} +
1903 \LAM\frac{\partial V^B(r_B)}{\partial r_B}
1904 \frac{\partial r_B}{\partial \LAM}
1907 V^B(r_B) -V^A(r_A) + \nonumber \\
1910 \left[ \LAM F^B(r_B) r^{-5}_B \sigma_B^6 \LL^{p-1} -
1911 \LL F^A(r_A) r^{-5}_A \sigma_A^6 \LAM^{p-1} \right]
1914 The original GROMOS Lennard-Jones soft-core function~\cite{Beutler94}
1915 uses $p=2$, but $p=1$ gives a smoother $\partial H/\partial\LAM$ curve.
1916 %When the changes between the two states involve both the disappearing
1917 %and appearing of atoms, it is important that the overlapping of atoms
1918 %happens around $\LAM=0.5$. This can usually be achieved with
1919 %$\alpha$$\approx0.7$ for $p=1$ and $\alpha$$\approx1.5$ for $p=2$.
1920 %MRS: this is now eliminated as of 4.6, since changes between atoms are done linearly.
1922 Another issue that should be considered is the soft-core effect of hydrogens
1923 without Lennard-Jones interaction. Their soft-core $\sigma$ is
1924 set with {\tt sc-sigma} in the {\tt .mdp} file. These hydrogens
1925 produce peaks in $\partial H/\partial\LAM$ at $\LAM$ is 0 and/or 1 for $p=1$
1926 and close to 0 and/or 1 with $p=2$. Lowering {\tt\mbox{sc-sigma}} will decrease
1927 this effect, but it will also increase the interactions with hydrogens
1928 relative to the other interactions in the soft-core state.
1930 When soft core potentials are selected (by setting {\tt sc-alpha} \textgreater
1931 0), and the Coulomb and Lennard-Jones potentials are turned on or off
1932 sequentially, then the Coulombic interaction is turned off linearly,
1933 rather than using soft core interactions, which should be less
1934 statistically noisy in most cases. This behavior can be overwritten
1935 by using the mdp option {\tt sc-coul} to {\tt yes}. Additionally, the
1936 soft-core interaction potential is only applied when either the A or B
1937 state has zero interaction potential. If both A and B states have
1938 nonzero interaction potential, default linear scaling described above
1939 is used. When both Coulombic and Lennard-Jones interactions are turned
1940 off simultaneously, a soft-core potential is used, and a hydrogen is
1941 being introduced or deleted, the sigma is set to {\tt sc-sigma-min},
1942 which itself defaults to {\tt sc-sigma-default}.
1944 Recently, a new formulation of the soft-core approach has been derived
1945 that in most cases gives lower and more even statistical variance than
1946 the standard soft-core path described above.~\cite{Pham2011,Pham2012}
1947 Specifically, we have:
1949 V_{sc}(r) &=& \LL V^A(r_A) + \LAM V^B(r_B)
1951 r_A &=& \left(\alpha \sigma_A^{48} \LAM^p + r^{48} \right)^\frac{1}{48}
1953 r_B &=& \left(\alpha \sigma_B^{48} \LL^p + r^{48} \right)^\frac{1}{48}
1955 This ``1-1-48'' path is also implemented in {\gromacs}. Note that for this path the soft core $\alpha$
1956 should satisfy $0.001 < \alpha < 0.003$, rather than $\alpha \approx
1959 %} % Brace matches ifthenelse test for gmxlite
1961 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1963 \subsection{Exclusions and 1-4 Interactions.}
1964 Atoms within a molecule that are close by in the chain,
1965 {\ie} atoms that are covalently bonded, or linked by one or two
1966 atoms are called {\em first neighbors, second neighbors} and
1967 {\em \swapindex{third}{neighbor}s}, respectively (see \figref{chain}).
1968 Since the interactions of atom {\bf i} with atoms {\bf i+1} and {\bf i+2}
1969 are mainly quantum mechanical, they can not be modeled by a Lennard-Jones potential.
1970 Instead it is assumed that these interactions are adequately modeled
1971 by a harmonic bond term or constraint ({\bf i, i+1}) and a harmonic angle term
1972 ({\bf i, i+2}). The first and second neighbors (atoms {\bf i+1} and {\bf i+2})
1974 {\em excluded} from the Lennard-Jones \swapindex{interaction}{list}
1976 atoms {\bf i+1} and {\bf i+2} are called {\em \normindex{exclusions}} of atom {\bf i}.
1979 \centerline{\includegraphics[width=8cm]{plots/chain}}
1980 \caption{Atoms along an alkane chain.}
1984 For third neighbors, the normal Lennard-Jones repulsion is sometimes
1985 still too strong, which means that when applied to a molecule, the
1986 molecule would deform or break due to the internal strain. This is
1987 especially the case for carbon-carbon interactions in a {\em
1988 cis}-conformation ({\eg} {\em cis}-butane). Therefore, for some of these
1989 interactions, the Lennard-Jones repulsion has been reduced in the
1990 {\gromos} force field, which is implemented by keeping a separate list of
1991 1-4 and normal Lennard-Jones parameters. In other force fields, such
1992 as OPLS~\cite{Jorgensen88}, the standard Lennard-Jones parameters are reduced
1993 by a factor of two, but in that case also the dispersion (r$^{-6}$)
1994 and the Coulomb interaction are scaled.
1995 {\gromacs} can use either of these methods.
1997 \subsection{Charge Groups\index{charge group}}
1999 In principle, the force calculation in MD is an $O(N^2)$ problem.
2000 Therefore, we apply a \normindex{cut-off} for non-bonded force (NBF)
2001 calculations; only the particles within a certain distance of each
2002 other are interacting. This reduces the cost to $O(N)$ (typically
2003 $100N$ to $200N$) of the NBF. It also introduces an error, which is,
2004 in most cases, acceptable, except when applying the cut-off implies
2005 the creation of charges, in which case you should consider using the
2006 lattice sum methods provided by {\gromacs}.
2008 Consider a water molecule interacting with another atom. If we would apply
2009 a plain cut-off on an atom-atom basis we might include the atom-oxygen
2010 interaction (with a charge of $-0.82$) without the compensating charge
2011 of the protons, and as a result, induce a large dipole moment over the system.
2012 Therefore, we have to keep groups of atoms with total charge
2013 0 together. These groups are called {\em charge groups}. Note that with
2014 a proper treatment of long-range electrostatics (e.g. particle-mesh Ewald
2015 (\secref{pme}), keeping charge groups together is not required.
2017 \subsection{Treatment of Cut-offs in the group scheme\index{cut-off}}
2018 \newcommand{\rs}{$r_{short}$}
2019 \newcommand{\rl}{$r_{long}$}
2020 {\gromacs} is quite flexible in treating cut-offs, which implies
2021 there can be quite a number of parameters to set. These parameters are
2022 set in the input file for {\tt grompp}. There are two sort of parameters
2023 that affect the cut-off interactions; you can select which type
2024 of interaction to use in each case, and which cut-offs should be
2025 used in the neighbor searching.
2027 For both Coulomb and van der Waals interactions there are interaction
2028 type selectors (termed {\tt vdwtype} and {\tt coulombtype}) and two
2029 parameters, for a total of six non-bonded interaction parameters. See
2030 \secref{mdpopt} for a complete description of these parameters.
2032 The neighbor searching (NS) can be performed using a single-range, or a twin-range
2033 approach. Since the former is merely a special case of the latter, we will
2034 discuss the more general twin-range. In this case, NS is described by two
2035 radii: {\tt rlist} and max({\tt rcoulomb},{\tt rvdw}).
2036 Usually one builds the neighbor list every 10 time steps
2037 or every 20 fs (parameter {\tt nstlist}). In the neighbor list, all interaction
2038 pairs that fall within {\tt rlist} are stored. Furthermore, the
2039 interactions between pairs that do not
2040 fall within {\tt rlist} but do fall within max({\tt rcoulomb},{\tt rvdw})
2041 are computed during NS. The
2042 forces and energy are stored separately and added to short-range forces
2043 at every time step between successive NS. If {\tt rlist} =
2044 max({\tt rcoulomb},{\tt rvdw}), no forces
2045 are evaluated during neighbor list generation.
2046 The \normindex{virial} is calculated from the sum of the short- and
2048 This means that the virial can be slightly asymmetrical at non-NS steps.
2049 When mdrun is compiled to use mixed precision, the virial is almost always asymmetrical because the
2050 off-diagonal elements are about as large as each element in the sum.
2051 In most cases this is not really a problem, since the fluctuations in the
2052 virial can be 2 orders of magnitude larger than the average.
2054 Except for the plain cut-off,
2055 all of the interaction functions in \tabref{funcparm}
2056 require that neighbor searching be done with a larger radius than the $r_c$
2057 specified for the functional form, because of the use of charge groups.
2058 The extra radius is typically of the order of 0.25 nm (roughly the
2059 largest distance between two atoms in a charge group plus the distance a
2060 charge group can diffuse within neighbor list updates).
2062 %If your charge groups are very large it may be interesting to turn off charge
2063 %groups, by setting the option
2064 %{\tt bAtomList = yes} in your {\tt grompp.mdp} file.
2065 %In this case only a small extra radius to account for diffusion needs to be
2066 %added (0.1 nm). Do not however use this together with the plain cut-off
2067 %method, as it will generate large artifacts (\secref{cg}).
2068 %In summary, there are four parameters that describe NS behavior:
2069 %{\tt nstlist} (update frequency in number of time steps),
2070 %{\tt bAtomList} (whether or not charge groups are used to generate neighbor list, the default is to use charge groups, so {\tt bAtomList = no}),
2071 %{\tt rshort} and {\tt rlong} which are the two radii {\rs} and {\rl}
2076 \begin{tabular}{|ll|l|}
2078 \multicolumn{2}{|c|}{Type} & Parameters \\
2080 Coulomb&Plain cut-off & $r_c$, $\epsr$ \\
2081 &Reaction field & $r_c$, $\epsrf$ \\
2082 &Shift function & $r_1$, $r_c$, $\epsr$ \\
2083 &Switch function & $r_1$, $r_c$, $\epsr$ \\
2085 VdW&Plain cut-off & $r_c$ \\
2086 &Shift function & $r_1$, $r_c$ \\
2087 &Switch function & $r_1$, $r_c$ \\
2090 \caption[Parameters for the different functional forms of the
2091 non-bonded interactions.]{Parameters for the different functional
2092 forms of the non-bonded interactions.}
2093 \label{tab:funcparm}
2095 %} % Brace matches ifthenelse test for gmxlite
2098 \newcommand{\vvis}{\ve{r}_s}
2099 \newcommand{\Fi}{\ve{F}_i'}
2100 \newcommand{\Fj}{\ve{F}_j'}
2101 \newcommand{\Fk}{\ve{F}_k'}
2102 \newcommand{\Fl}{\ve{F}_l'}
2103 \newcommand{\Fvis}{\ve{F}_{s}}
2104 \newcommand{\rvik}{\ve{r}_{ik}}
2105 \newcommand{\rvis}{\ve{r}_{is}}
2106 \newcommand{\rvjk}{\ve{r}_{jk}}
2107 \newcommand{\rvjl}{\ve{r}_{jl}}
2109 %\ifthenelse{\equal{\gmxlite}{1}}{}{
2110 \section{Virtual interaction sites\index{virtual interaction sites}}
2111 \label{sec:virtual_sites}
2112 Virtual interaction sites (called \seeindex{dummy atoms}{virtual interaction sites} in {\gromacs} versions before 3.3)
2113 can be used in {\gromacs} in a number of ways.
2114 We write the position of the virtual site $\ve{r}_s$ as a function of
2115 the positions of other particles \ve{r}$_i$: $\ve{r}_s =
2116 f(\ve{r}_1..\ve{r}_n)$. The virtual site, which may carry charge or be
2117 involved in other interactions, can now be used in the force
2118 calculation. The force acting on the virtual site must be
2119 redistributed over the particles with mass in a consistent way.
2120 A good way to do this can be found in ref.~\cite{Berendsen84b}.
2121 We can write the potential energy as:
2123 V = V(\vvis,\ve{r}_1,\ldots,\ve{r}_n) = V^*(\ve{r}_1,\ldots,\ve{r}_n)
2125 The force on the particle $i$ is then:
2127 \ve{F}_i = -\frac{\partial V^*}{\partial \ve{r}_i}
2128 = -\frac{\partial V}{\partial \ve{r}_i} -
2129 \frac{\partial V}{\partial \vvis}
2130 \frac{\partial \vvis}{\partial \ve{r}_i}
2131 = \ve{F}_i^{direct} + \Fi
2133 The first term is the normal force.
2134 The second term is the force on particle $i$ due to the virtual site, which
2135 can be written in tensor notation:
2136 \newcommand{\partd}[2]{\displaystyle\frac{\partial #1}{\partial #2_i}}
2138 \Fi = \left[\begin{array}{ccc}
2139 \partd{x_s}{x} & \partd{y_s}{x} & \partd{z_s}{x} \\[1ex]
2140 \partd{x_s}{y} & \partd{y_s}{y} & \partd{z_s}{y} \\[1ex]
2141 \partd{x_s}{z} & \partd{y_s}{z} & \partd{z_s}{z}
2142 \end{array}\right]\Fvis
2145 where $\Fvis$ is the force on the virtual site and $x_s$, $y_s$ and
2146 $z_s$ are the coordinates of the virtual site. In this way, the total
2147 force and the total torque are conserved~\cite{Berendsen84b}.
2149 The computation of the \normindex{virial}
2150 (\eqnref{Xi}) is non-trivial when virtual sites are used. Since the
2151 virial involves a summation over all the atoms (rather than virtual
2152 sites), the forces must be redistributed from the virtual sites to the
2153 atoms (using ~\eqnref{fvsite}) {\em before} computation of the
2154 virial. In some special cases where the forces on the atoms can be
2155 written as a linear combination of the forces on the virtual sites (types 2
2156 and 3 below) there is no difference between computing the virial
2157 before and after the redistribution of forces. However, in the
2158 general case redistribution should be done first.
2161 \centerline{\includegraphics[width=15cm]{plots/dummies}}
2162 \caption[Virtual site construction.]{The six different types of virtual
2163 site construction in \protect{\gromacs}. The constructing atoms are
2164 shown as black circles, the virtual sites in gray.}
2168 There are six ways to construct virtual sites from surrounding atoms in
2169 {\gromacs}, which we classify by the number of constructing
2170 atoms. {\bf Note} that all site types mentioned can be constructed from
2171 types 3fd (normalized, in-plane) and 3out (non-normalized, out of
2172 plane). However, the amount of computation involved increases sharply
2173 along this list, so we strongly recommended using the first adequate
2174 virtual site type that will be sufficient for a certain purpose.
2175 \figref{vsites} depicts 6 of the available virtual site constructions.
2176 The conceptually simplest construction types are linear combinations:
2178 \vvis = \sum_{i=1}^N w_i \, \ve{r}_i
2180 The force is then redistributed using the same weights:
2185 The types of virtual sites supported in {\gromacs} are given in the list below.
2186 Constructing atoms in virtual sites can be virtual sites themselves, but
2187 only if they are higher in the list, i.e. virtual sites can be
2188 constructed from ``particles'' that are simpler virtual sites.
2190 \item[{\bf\sf 2.}]\label{subsec:vsite2}As a linear combination of two atoms
2191 (\figref{vsites} 2):
2193 w_i = 1 - a ~,~~ w_j = a
2195 In this case the virtual site is on the line through atoms $i$ and
2198 \item[{\bf\sf 3.}]\label{subsec:vsite3}As a linear combination of three atoms
2199 (\figref{vsites} 3):
2201 w_i = 1 - a - b ~,~~ w_j = a ~,~~ w_k = b
2203 In this case the virtual site is in the plane of the other three
2206 \item[{\bf\sf 3fd.}]\label{subsec:vsite3fd}In the plane of three atoms, with a fixed distance
2207 (\figref{vsites} 3fd):
2209 \vvis ~=~ \ve{r}_i + b \frac{ \rvij + a \rvjk }
2210 {| \rvij + a \rvjk |}
2212 In this case the virtual site is in the plane of the other three
2213 particles at a distance of $|b|$ from $i$.
2214 The force on particles $i$, $j$ and $k$ due to the force on the virtual
2215 site can be computed as:
2218 \Fi &=& \displaystyle \Fvis - \gamma ( \Fvis - \ve{p} ) \\[1ex]
2219 \Fj &=& \displaystyle (1-a)\gamma (\Fvis - \ve{p}) \\[1ex]
2220 \Fk &=& \displaystyle a \gamma (\Fvis - \ve{p}) \\
2224 \displaystyle \gamma = \frac{b}{| \rvij + a \rvjk |} \\[2ex]
2225 \displaystyle \ve{p} = \frac{ \rvis \cdot \Fvis }
2226 { \rvis \cdot \rvis } \rvis
2230 \item[{\bf\sf 3fad.}]\label{subsec:vsite3fad}In the plane of three atoms, with a fixed angle and
2231 distance (\figref{vsites} 3fad):
2233 \label{eqn:vsite2fad-F}
2234 \vvis ~=~ \ve{r}_i +
2235 d \cos \theta \frac{\rvij}{|\rvij|} +
2236 d \sin \theta \frac{\ve{r}_\perp}{|\ve{r}_\perp|}
2238 \ve{r}_\perp ~=~ \rvjk -
2239 \frac{ \rvij \cdot \rvjk }
2240 { \rvij \cdot \rvij }
2243 In this case the virtual site is in the plane of the other three
2244 particles at a distance of $|d|$ from $i$ at an angle of
2245 $\alpha$ with $\rvij$. Atom $k$ defines the plane and the
2246 direction of the angle. {\bf Note} that in this case $b$ and
2247 $\alpha$ must be specified, instead of $a$ and $b$ (see also
2248 \secref{vsitetop}). The force on particles $i$, $j$ and $k$
2249 due to the force on the virtual site can be computed as (with
2250 $\ve{r}_\perp$ as defined in \eqnref{vsite2fad-F}):
2251 \newcommand{\dfrac}{\displaystyle\frac}
2254 \begin{array}{lclllll}
2256 \dfrac{d \cos \theta}{|\rvij|} \ve{F}_1 &+&
2257 \dfrac{d \sin \theta}{|\ve{r}_\perp|} \left(
2258 \dfrac{ \rvij \cdot \rvjk }
2259 { \rvij \cdot \rvij } \ve{F}_2 +
2260 \ve{F}_3 \right) \\[3ex]
2262 \dfrac{d \cos \theta}{|\rvij|} \ve{F}_1 &-&
2263 \dfrac{d \sin \theta}{|\ve{r}_\perp|} \left(
2265 \dfrac{ \rvij \cdot \rvjk }
2266 { \rvij \cdot \rvij } \ve{F}_2 +
2267 \ve{F}_3 \right) \\[3ex]
2269 \dfrac{d \sin \theta}{|\ve{r}_\perp|} \ve{F}_2 \\[3ex]
2273 \dfrac{ \rvij \cdot \Fvis }
2274 { \rvij \cdot \rvij } \rvij
2276 \ve{F}_2 = \ve{F}_1 -
2277 \dfrac{ \ve{r}_\perp \cdot \Fvis }
2278 { \ve{r}_\perp \cdot \ve{r}_\perp } \ve{r}_\perp
2280 \ve{F}_3 = \dfrac{ \rvij \cdot \Fvis }
2281 { \rvij \cdot \rvij } \ve{r}_\perp
2285 \item[{\bf\sf 3out.}]\label{subsec:vsite3out}As a non-linear combination of three atoms, out of plane
2286 (\figref{vsites} 3out):
2288 \vvis ~=~ \ve{r}_i + a \rvij + b \rvik +
2289 c (\rvij \times \rvik)
2291 This enables the construction of virtual sites out of the plane of the
2293 The force on particles $i,j$ and $k$ due to the force on the virtual
2294 site can be computed as:
2298 \Fj &=& \left[\begin{array}{ccc}
2299 a & -c\,z_{ik} & c\,y_{ik} \\[0.5ex]
2300 c\,z_{ik} & a & -c\,x_{ik} \\[0.5ex]
2301 -c\,y_{ik} & c\,x_{ik} & a
2302 \end{array}\right]\Fvis \\
2304 \Fk &=& \left[\begin{array}{ccc}
2305 b & c\,z_{ij} & -c\,y_{ij} \\[0.5ex]
2306 -c\,z_{ij} & b & c\,x_{ij} \\[0.5ex]
2307 c\,y_{ij} & -c\,x_{ij} & b
2308 \end{array}\right]\Fvis \\
2309 \Fi &=& \Fvis - \Fj - \Fk
2313 \item[{\bf\sf 4fdn.}]\label{subsec:vsite4fdn}From four atoms, with a fixed distance, see separate Fig.\ \ref{fig:vsite-4fdn}.
2314 This construction is a bit
2315 complex, in particular since the previous type (4fd) could be unstable which forced us
2316 to introduce a more elaborate construction:
2319 \centerline{\includegraphics[width=5cm]{plots/vsite-4fdn}}
2320 \caption {The new 4fdn virtual site construction, which is stable even when all constructing
2321 atoms are in the same plane.}
2322 \label{fig:vsite-4fdn}
2326 \mathbf{r}_{ja} &=& a\, \mathbf{r}_{ik} - \mathbf{r}_{ij} = a\, (\mathbf{x}_k - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\
2327 \mathbf{r}_{jb} &=& b\, \mathbf{r}_{il} - \mathbf{r}_{ij} = b\, (\mathbf{x}_l - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\
2328 \mathbf{r}_m &=& \mathbf{r}_{ja} \times \mathbf{r}_{jb} \nonumber \\
2329 \mathbf{x}_s &=& \mathbf{x}_i + c \frac{\mathbf{r}_m}{|\mathbf{r}_m|}
2333 In this case the virtual site is at a distance of $|c|$ from $i$, while $a$ and $b$ are
2334 parameters. {\bf Note} that the vectors $\mathbf{r}_{ik}$ and $\mathbf{r}_{ij}$ are not normalized
2335 to save floating-point operations.
2336 The force on particles $i$, $j$, $k$ and $l$ due to the force
2337 on the virtual site are computed through chain rule derivatives
2338 of the construction expression. This is exact and conserves energy,
2339 but it does lead to relatively lengthy expressions that we do not
2340 include here (over 200 floating-point operations). The interested reader can
2341 look at the source code in \verb+vsite.c+. Fortunately, this vsite type is normally
2342 only used for chiral centers such as $C_{\alpha}$ atoms in proteins.
2344 The new 4fdn construct is identified with a `type' value of 2 in the topology. The earlier 4fd
2345 type is still supported internally (`type' value 1), but it should not be used for
2346 new simulations. All current {\gromacs} tools will automatically generate type 4fdn instead.
2349 \item[{\bf\sf N.}]\label{subsec:vsiteN} A linear combination of $N$ atoms with relative
2350 weights $a_i$. The weight for atom $i$ is:
2352 w_i = a_i \left(\sum_{j=1}^N a_j \right)^{-1}
2354 There are three options for setting the weights:
2356 \item[COG] center of geometry: equal weights
2357 \item[COM] center of mass: $a_i$ is the mass of atom $i$;
2358 when in free-energy simulations the mass of the atom is changed,
2359 only the mass of the A-state is used for the weight
2360 \item[COW] center of weights: $a_i$ is defined by the user
2364 %} % Brace matches ifthenelse test for gmxlite
2366 \newcommand{\dr}{{\rm d}r}
2367 \newcommand{\avcsix}{\left< C_6 \right>}
2369 %\ifthenelse{\equal{\gmxlite}{1}}{}{
2370 \section{Long Range Electrostatics}
2371 \label{sec:lr_elstat}
2372 \subsection{Ewald summation\index{Ewald sum}}
2374 The total electrostatic energy of $N$ particles and their periodic
2375 images\index{periodic boundary conditions} is given by
2377 V=\frac{f}{2}\sum_{n_x}\sum_{n_y}
2378 \sum_{n_{z}*} \sum_{i}^{N} \sum_{j}^{N}
2379 \frac{q_i q_j}{{\bf r}_{ij,{\bf n}}}.
2380 \label{eqn:totalcoulomb}
2382 $(n_x,n_y,n_z)={\bf n}$ is the box index vector, and the star indicates that
2383 terms with $i=j$ should be omitted when $(n_x,n_y,n_z)=(0,0,0)$. The
2384 distance ${\bf r}_{ij,{\bf n}}$ is the real distance between the charges and
2385 not the minimum-image. This sum is conditionally convergent, but
2388 Ewald summation was first introduced as a method to calculate
2389 long-range interactions of the periodic images in
2390 crystals~\cite{Ewald21}. The idea is to convert the single
2391 slowly-converging sum \eqnref{totalcoulomb} into two
2392 quickly-converging terms and a constant term:
2394 V &=& V_{\mathrm{dir}} + V_{\mathrm{rec}} + V_{0} \\[0.5ex]
2395 V_{\mathrm{dir}} &=& \frac{f}{2} \sum_{i,j}^{N}
2396 \sum_{n_x}\sum_{n_y}
2397 \sum_{n_{z}*} q_i q_j \frac{\mbox{erfc}(\beta {r}_{ij,{\bf n}} )}{{r}_{ij,{\bf n}}} \\[0.5ex]
2398 V_{\mathrm{rec}} &=& \frac{f}{2 \pi V} \sum_{i,j}^{N} q_i q_j
2399 \sum_{m_x}\sum_{m_y}
2400 \sum_{m_{z}*} \frac{\exp{\left( -(\pi {\bf m}/\beta)^2 + 2 \pi i
2401 {\bf m} \cdot ({\bf r}_i - {\bf r}_j)\right)}}{{\bf m}^2} \\[0.5ex]
2402 V_{0} &=& -\frac{f \beta}{\sqrt{\pi}}\sum_{i}^{N} q_i^2,
2404 where $\beta$ is a parameter that determines the relative weight of the
2405 direct and reciprocal sums and ${\bf m}=(m_x,m_y,m_z)$.
2406 In this way we can use a short cut-off (of the order of $1$~nm) in the direct space sum and a
2407 short cut-off in the reciprocal space sum ({\eg} 10 wave vectors in each
2408 direction). Unfortunately, the computational cost of the reciprocal
2409 part of the sum increases as $N^2$
2410 (or $N^{3/2}$ with a slightly better algorithm) and it is therefore not
2411 realistic for use in large systems.
2413 \subsubsection{Using Ewald}
2414 Don't use Ewald unless you are absolutely sure this is what you want -
2415 for almost all cases the PME method below will perform much better.
2416 If you still want to employ classical Ewald summation enter this in
2417 your {\tt .mdp} file, if the side of your box is about $3$~nm:
2424 fourierspacing = 0.6
2428 The ratio of the box dimensions and the {\tt fourierspacing} parameter determines
2429 the highest magnitude of wave vectors $m_x,m_y,m_z$ to use in each
2430 direction. With a 3-nm cubic box this example would use $11$ wave vectors
2431 (from $-5$ to $5$) in each direction. The {\tt ewald-rtol} parameter
2432 is the relative strength of the electrostatic interaction at the
2433 cut-off. Decreasing this gives you a more accurate direct sum, but a
2434 less accurate reciprocal sum.
2436 \subsection{\normindex{PME}}
2438 Particle-mesh Ewald is a method proposed by Tom
2439 Darden~\cite{Darden93} to improve the performance of the
2440 reciprocal sum. Instead of directly summing wave vectors, the charges
2441 are assigned to a grid using interpolation. The implementation in
2442 {\gromacs} uses cardinal B-spline interpolation~\cite{Essmann95},
2443 which is referred to as smooth PME (SPME).
2444 The grid is then Fourier transformed with a 3D FFT algorithm and the
2445 reciprocal energy term obtained by a single sum over the grid in
2448 The potential at the grid points is calculated by inverse
2449 transformation, and by using the interpolation factors we get the
2450 forces on each atom.
2452 The PME algorithm scales as $N \log(N)$, and is substantially faster
2453 than ordinary Ewald summation on medium to large systems. On very
2454 small systems it might still be better to use Ewald to avoid the
2455 overhead in setting up grids and transforms.
2456 For the parallelization of PME see the section on MPMD PME (\ssecref{mpmd_pme}).
2458 With the Verlet cut-off scheme, the PME direct space potential is
2459 shifted by a constant such that the potential is zero at the
2460 cut-off. This shift is small and since the net system charge is close
2461 to zero, the total shift is very small, unlike in the case of the
2462 Lennard-Jones potential where all shifts add up. We apply the shift
2463 anyhow, such that the potential is the exact integral of the force.
2465 \subsubsection{Using PME}
2466 As an example for using Particle-mesh Ewald summation in {\gromacs}, specify the
2467 following lines in your {\tt .mdp} file:
2474 fourierspacing = 0.12
2479 In this case the {\tt fourierspacing} parameter determines the maximum
2480 spacing for the FFT grid (i.e. minimum number of grid points),
2481 and {\tt pme-order} controls the
2482 interpolation order. Using fourth-order (cubic) interpolation and this
2483 spacing should give electrostatic energies accurate to about
2484 $5\cdot10^{-3}$. Since the Lennard-Jones energies are not this
2485 accurate it might even be possible to increase this spacing slightly.
2487 Pressure scaling works with PME, but be aware of the fact that
2488 anisotropic scaling can introduce artificial ordering in some systems.
2490 \subsection{\normindex{P3M-AD}}
2492 The \seeindex{Particle-Particle Particle-Mesh}{P3M} methods of
2493 Hockney \& Eastwood can also be applied in {\gromacs} for the
2494 treatment of long range electrostatic interactions~\cite{Hockney81}.
2495 Although the P3M method was the first efficient long-range electrostatics
2496 method for molecular simulation, the smooth PME (SPME) method has largely
2497 replaced P3M as the method of choice in atomistic simulations. One performance
2498 disadvantage of the original P3M method was that it required 3 3D-FFT
2499 back transforms to obtain the forces on the particles. But this is not
2500 required for P3M and the forces can be derived through analytical differentiation
2501 of the potential, as done in PME. The resulting method is termed P3M-AD.
2502 The only remaining difference between P3M-AD and PME is the optimization
2503 of the lattice Green influence function for error minimization that P3M uses.
2504 However, in 2012 it has been shown that the SPME influence function can be
2505 modified to obtain P3M~\cite{Ballenegger2012}.
2506 This means that the advantage of error minimization in P3M-AD can be used
2507 at the same computational cost and with the same code as PME,
2508 just by adding a few lines to modify the influence function.
2509 However, at optimal parameter setting the effect of error minimization
2510 in P3M-AD is less than 10\%. P3M-AD does show large accuracy gains with
2511 interlaced (also known as staggered) grids, but that is not supported
2512 in {\gromacs} (yet).
2514 P3M is used in {\gromacs} with exactly the same options as used with PME
2515 by selecting the electrostatics type:
2517 coulombtype = P3M-AD
2520 \subsection{Optimizing Fourier transforms and PME calculations}
2521 It is recommended to optimize the parameters for calculation of
2522 electrostatic interaction such as PME grid dimensions and cut-off radii.
2523 This is particularly relevant to do before launching long production runs.
2525 {\gromacs} includes a special tool, {\tt g_tune_pme}, which automates the
2526 process of selecting the optimal size of the grid and number of PME-only
2530 % Temporarily removed since I am not sure about the state of the testlr
2533 %It is possible to test the accuracy of your settings using the program
2534 %{\tt\normindex{testlr}} in the {\tt src/gmxlib} dir. This program computes
2535 %forces and potentials using PPPM and an Ewald implementation and gives the
2536 %absolute and RMS errors in both:
2541 %Potential 0.113 0.035
2543 %{\bf Note:} these numbers were generated using a grid spacing of
2544 %0.058 nm and $r_c$ = 1.0 nm.
2546 %You can see what the accuracy is without optimizing the
2547 %$\hat{G}(k)$ function, if you pass the {\tt -ghat} option to {\tt
2548 %testlr}. Try it if you think the {\tt mk_ghat} procedure is a waste
2550 %} % Brace matches ifthenelse test for gmxlite
2553 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2554 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2555 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2557 %\ifthenelse{\equal{\gmxlite}{1}}{}{
2558 \section{Long Range Van der Waals interactions}
2559 \subsection{Dispersion correction\index{dispersion correction}}
2560 In this section, we derive long-range corrections due to the use of a
2561 cut-off for Lennard-Jones or Buckingham interactions.
2562 We assume that the cut-off is
2563 so long that the repulsion term can safely be neglected, and therefore
2564 only the dispersion term is taken into account. Due to the nature of
2565 the dispersion interaction (we are truncating a potential proportional
2566 to $-r^{-6}$), energy and pressure corrections are both negative. While
2567 the energy correction is usually small, it may be important for free
2568 energy calculations where differences between two different Hamiltonians
2569 are considered. In contrast, the pressure correction is very large and
2570 can not be neglected under any circumstances where a correct pressure is
2571 required, especially for any NPT simulations. Although it is, in
2572 principle, possible to parameterize a force field such that the pressure
2573 is close to the desired experimental value without correction, such a
2574 method makes the parameterization dependent on the cut-off and is therefore
2577 \subsubsection{Energy}
2579 The long-range contribution of the dispersion interaction to the
2580 virial can be derived analytically, if we assume a homogeneous
2581 system beyond the cut-off distance $r_c$. The dispersion energy
2582 between two particles is written as:
2584 V(\rij) ~=~- C_6\,\rij^{-6}
2586 and the corresponding force is:
2588 \Fvij ~=~- 6\,C_6\,\rij^{-8}\rvij
2590 In a periodic system it is not easy to calculate the full potentials,
2591 so usually a cut-off is applied, which can be abrupt or smooth.
2592 We will call the potential and force with cut-off $V_c$ and $\ve{F}_c$.
2593 The long-range contribution to the dispersion energy
2594 in a system with $N$ particles and particle density $\rho$ = $N/V$ is:
2596 \label{eqn:enercorr}
2597 V_{lr} ~=~ \half N \rho\int_0^{\infty} 4\pi r^2 g(r) \left( V(r) -V_c(r) \right) {\dr}
2599 We will integrate this for the shift function, which is the most general
2600 form of van der Waals interaction available in {\gromacs}.
2601 The shift function has a constant difference $S$ from 0 to $r_1$
2602 and is 0 beyond the cut-off distance $r_c$.
2603 We can integrate \eqnref{enercorr}, assuming that the density in the sphere
2604 within $r_1$ is equal to the global density and
2605 the radial distribution function $g(r)$ is 1 beyond $r_1$:
2608 V_{lr} &=& \half N \left(
2609 \rho\int_0^{r_1} 4\pi r^2 g(r) \, C_6 \,S\,{\dr}
2610 + \rho\int_{r_1}^{r_c} 4\pi r^2 \left( V(r) -V_c(r) \right) {\dr}
2611 + \rho\int_{r_c}^{\infty} 4\pi r^2 V(r) \, {\dr}
2613 & = & \half N \left(\left(\frac{4}{3}\pi \rho r_1^{3} - 1\right) C_6 \,S
2614 + \rho\int_{r_1}^{r_c} 4\pi r^2 \left( V(r) -V_c(r) \right) {\dr}
2615 -\frac{4}{3} \pi N \rho\, C_6\,r_c^{-3}
2618 where the term $-1$ corrects for the self-interaction.
2619 For a plain cut-off we only need to assume that $g(r)$ is 1 beyond $r_c$
2620 and the correction reduces to~\cite{Allen87}:
2622 V_{lr} & = & -\frac{2}{3} \pi N \rho\, C_6\,r_c^{-3}
2624 If we consider, for example, a box of pure water, simulated with a cut-off
2625 of 0.9 nm and a density of 1 g cm$^{-3}$ this correction is
2626 $-0.75$ kJ mol$^{-1}$ per molecule.
2628 For a homogeneous mixture we need to define
2629 an {\em average dispersion constant}:
2632 \avcsix = \frac{2}{N(N-1)}\sum_i^N\sum_{j>i}^N C_6(i,j)\\
2634 In {\gromacs}, excluded pairs of atoms do not contribute to the average.
2636 In the case of inhomogeneous simulation systems, {\eg} a system with a
2637 lipid interface, the energy correction can be applied if
2638 $\avcsix$ for both components is comparable.
2640 \subsubsection{Virial and pressure}
2641 The scalar virial of the system due to the dispersion interaction between
2642 two particles $i$ and $j$ is given by:
2644 \Xi~=~-\half \rvij \cdot \Fvij ~=~ 3\,C_6\,\rij^{-6}
2646 The pressure is given by:
2648 P~=~\frac{2}{3\,V}\left(E_{kin} - \Xi\right)
2650 The long-range correction to the virial is given by:
2652 \Xi_{lr} ~=~ \half N \rho \int_0^{\infty} 4\pi r^2 g(r) (\Xi -\Xi_c) \,\dr
2654 We can again integrate the long-range contribution to the
2655 virial assuming $g(r)$ is 1 beyond $r_1$:
2657 \Xi_{lr}&=& \half N \rho \left(
2658 \int_{r_1}^{r_c} 4 \pi r^2 (\Xi -\Xi_c) \,\dr
2659 + \int_{r_c}^{\infty} 4 \pi r^2 3\,C_6\,\rij^{-6}\, \dr
2661 &=& \half N \rho \left(
2662 \int_{r_1}^{r_c} 4 \pi r^2 (\Xi -\Xi_c) \, \dr
2663 + 4 \pi C_6 \, r_c^{-3} \right)
2665 For a plain cut-off the correction to the pressure is~\cite{Allen87}:
2667 P_{lr}~=~-\frac{4}{3} \pi C_6\, \rho^2 r_c^{-3}
2669 Using the same example of a water box, the correction to the virial is
2670 0.75 kJ mol$^{-1}$ per molecule,
2671 the corresponding correction to the pressure for
2672 SPC water is approximately $-280$ bar.
2674 For homogeneous mixtures, we can again use the average dispersion constant
2675 $\avcsix$ (\eqnref{avcsix}):
2677 P_{lr}~=~-\frac{4}{3} \pi \avcsix \rho^2 r_c^{-3}
2680 For inhomogeneous systems, \eqnref{pcorr} can be applied under the same
2681 restriction as holds for the energy (see \secref{ecorr}).
2683 \subsection{Lennard-Jones PME\index{LJ-PME}}
2685 In order to treat systems, using Lennard-Jones potentials, that are
2686 non-homogeneous outside of the cut-off distance, we can instead use
2687 the Particle-mesh Ewald method as discussed for electrostatics above.
2688 In this case the modified Ewald equations become
2690 V &=& V_{\mathrm{dir}} + V_{\mathrm{rec}} + V_{0} \\[0.5ex]
2691 V_{\mathrm{dir}} &=& -\frac{1}{2} \sum_{i,j}^{N}
2692 \sum_{n_x}\sum_{n_y}
2693 \sum_{n_{z}*} \frac{C^{ij}_6 g(\beta {r}_{ij,{\bf n}})}{{r_{ij,{\bf n}}}^6}
2694 \label{eqn:ljpmerealspace}\\[0.5ex]
2695 V_{\mathrm{rec}} &=& \frac{{\pi}^{\frac{3}{2}} \beta^{3}}{2V} \sum_{m_x}\sum_{m_y}\sum_{m_{z}*}
2696 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]
2697 V_{0} &=& -\frac{\beta^{6}}{12}\sum_{i}^{N} C^{ii}_6
2700 where ${\bf m}=(m_x,m_y,m_z)$, $\beta$ is the parameter determining the weight between
2701 direct and reciprocal space, and ${C^{ij}_6}$ is the combined dispersion
2702 parameter for particle $i$ and $j$. The star indicates that terms
2703 with $i = j$ should be omitted when $((n_x,n_y,n_z)=(0,0,0))$, and
2704 ${\bf r}_{ij,{\bf n}}$ is the real distance between the particles.
2705 Following the derivation by Essmann~\cite{Essmann95}, the functions $f$ and $g$ introduced above are defined as
2707 f(x)&=&1/3\left[(1-2x^2){\mathrm{exp}}(-x^2) + 2{x^3}\sqrt{\pi}\,{\mathrm{erfc}}(x) \right] \\
2708 g(x)&=&{\mathrm{exp}}(-x^2)(1+x^2+\frac{x^4}{2}).
2711 The above methodology works fine as long as the dispersion parameters can be combined geometrically (\eqnref{comb}) in the same
2712 way as the charges for electrostatics
2714 C^{ij}_{6,\mathrm{geom}} = \left(C^{ii}_6 \, C^{jj}_6\right)^{1/2}
2716 For Lorentz-Berthelot combination rules (\eqnref{lorentzberthelot}), the reciprocal part of this sum has to be calculated
2717 seven times due to the splitting of the dispersion parameter according to
2719 C^{ij}_{6,\mathrm{L-B}} = (\sigma_i+\sigma_j)^6=\sum_{n=0}^{6} P_{n}\sigma_{i}^{n}\sigma_{j}^{(6-n)},
2721 for $P_{n}$ the Pascal triangle coefficients. This introduces a
2722 non-negligible cost to the reciprocal part, requiring seven separate
2723 FFTs, and therefore this has been the limiting factor in previous
2724 attempts to implement LJ-PME. A solution to this problem is to use
2725 geometrical combination rules in order to calculate an approximate
2726 interaction parameter for the reciprocal part of the potential,
2727 yielding a total interaction of
2729 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 \\
2730 &=& 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} \\
2731 V(r>r_c) & = & \underbrace{C^\mathrm{recip}_{6,\mathrm{geom}} [1 - g(\beta r)] r^{-6}}_{\mathrm{Reciprocal \ space}}.
2733 This will preserve a well-defined Hamiltonian and significantly increase
2734 the performance of the simulations. The approximation does introduce
2735 some errors, but since the difference is located in the interactions
2736 calculated in reciprocal space, the effect will be very small compared
2737 to the total interaction energy. In a simulation of a lipid bilayer,
2738 using a cut-off of 1.0 nm, the relative error in total dispersion
2739 energy was below 0.5\%. A more thorough discussion of this can be
2740 found in \cite{Wennberg13}.
2742 In {\gromacs} we now perform the proper calculation of this interaction
2743 by subtracting, from the direct-space interactions, the contribution
2744 made by the approximate potential that is used in the reciprocal part
2746 V_\mathrm{dir} = C^{\mathrm{dir}}_6 r^{-6} - C^\mathrm{recip}_6 [1 - g(\beta r)] r^{-6}.
2747 \label{eqn:ljpmedirectspace}
2749 This potential will reduce to the expression in \eqnref{ljpmerealspace} when $C^{\mathrm{dir}}_6 = C^\mathrm{recip}_6$,
2750 and the total interaction is given by
2752 \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}} \\
2753 &=&C^{\mathrm{dir}}_6 r^{-6}
2754 \label {eqn:ljpmecorr2} \\
2755 V(r>r_c) &=& C^\mathrm{recip}_6 [1 - g(\beta r)] r^{-6}.
2757 For the case when $C^{\mathrm{dir}}_6 \neq C^\mathrm{recip}_6$ this
2758 will retain an unmodified LJ force up to the cut-off, and the error
2759 is an order of magnitude smaller than in simulations where the
2760 direct-space interactions do not account for the approximation used in
2761 reciprocal space. When using a VdW interaction modifier of
2762 potential-shift, the constant
2764 \left(-C^{\mathrm{dir}}_6 + C^\mathrm{recip}_6 [1 - g(\beta r_c)]\right) r_c^{-6}
2766 is added to \eqnref{ljpmecorr2} in order to ensure that the potential
2767 is continuous at the cutoff. Note that, in the same way as \eqnref{ljpmedirectspace}, this degenerates into the
2768 expected $-C_6g(\beta r_c)r^{-6}_c$ when $C^{\mathrm{dir}}_6 =
2769 C^\mathrm{recip}_6$. In addition to this, a long-range dispersion
2770 correction can be applied to correct for the approximation using a
2771 combination rule in reciprocal space. This correction assumes, as for
2772 the cut-off LJ potential, a uniform particle distribution. But since
2773 the error of the combination rule approximation is very small this
2774 long-range correction is not necessary in most cases. Also note that
2775 this homogenous correction does not correct the surface tension, which
2776 is an inhomogeneous property.
2778 \subsubsection{Using LJ-PME}
2779 As an example for using Particle-mesh Ewald summation for Lennard-Jones interactions in {\gromacs}, specify the
2780 following lines in your {\tt .mdp} file:
2784 vdw-modifier = Potential-Shift
2787 fourierspacing = 0.12
2789 ewald-rtol-lj = 0.001
2790 lj-pme-comb-rule = geometric
2793 The same Fourier grid and interpolation order are used if both
2794 LJ-PME and electrostatic PME are active, so the settings for
2795 {\tt fourierspacing} and {\tt pme-order} are common to both.
2796 {\tt ewald-rtol-lj} controls the
2797 splitting between direct and reciprocal space in the same way as
2798 {\tt ewald-rtol}. In addition to this, the combination rule to be used
2799 in reciprocal space is determined by {\tt lj-pme-comb-rule}. If the
2800 current force field uses Lorentz-Berthelot combination rules, it is
2801 possible to set {\tt lj-pme-comb-rule = geometric} in order to gain a
2802 significant increase in performance for a small loss in accuracy. The
2803 details of this approximation can be found in the section above.
2805 Note that the use of a complete long-range dispersion correction means
2806 that as with Coulomb PME, {\tt rvdw} is now a free parameter in the
2807 method, rather than being necessarily restricted by the force-field
2808 parameterization scheme. Thus it is now possible to optimize the
2809 cutoff, spacing, order and tolerance terms for accuracy and best
2812 Naturally, the use of LJ-PME rather than LJ cut-off adds computation
2813 and communication done for the reciprocal-space part, so for best
2814 performance in balancing the load of parallel simulations using
2815 PME-only ranks, more such ranks should be used. It may be possible to
2816 improve upon the automatic load-balancing used by {\tt mdrun}.
2818 %} % Brace matches ifthenelse test for gmxlite
2820 \section{Force field\index{force field}}
2822 A force field is built up from two distinct components:
2824 \item The set of equations (called the {\em
2825 \index{potential function}s}) used to generate the potential
2826 energies and their derivatives, the forces. These are described in
2827 detail in the previous chapter.
2828 \item The parameters used in this set of equations. These are not
2829 given in this manual, but in the data files corresponding to your
2830 {\gromacs} distribution.
2832 Within one set of equations various sets of parameters can be
2833 used. Care must be taken that the combination of equations and
2834 parameters form a consistent set. It is in general dangerous to make
2835 {\em ad hoc} changes in a subset of parameters, because the various
2836 contributions to the total force are usually interdependent. This
2837 means in principle that every change should be documented, verified by
2838 comparison to experimental data and published in a peer-reviewed
2839 journal before it can be used.
2841 {\gromacs} {\gmxver} includes several force fields, and additional
2842 ones are available on the website. If you do not know which one to
2843 select we recommend \gromosv{96} for united-atom setups and OPLS-AA/L for
2844 all-atom parameters. That said, we describe the available options in
2847 \subsubsection{All-hydrogen force field}
2848 The \gromosv{87}-based all-hydrogen force field is almost identical to the
2849 normal \gromosv{87} force field, since the extra hydrogens have no
2850 Lennard-Jones interaction and zero charge. The only differences are in
2851 the bond angle and improper dihedral angle terms. This force field is
2852 only useful when you need the exact hydrogen positions, for instance
2853 for distance restraints derived from NMR measurements. When citing
2854 this force field please read the previous paragraph.
2856 \subsection{\gromosv{96}\index{GROMOS96 force field}}
2857 {\gromacs} supports the \gromosv{96} force fields~\cite{gromos96}.
2858 All parameters for the 43A1, 43A2 (development, improved alkane
2859 dihedrals), 45A3, 53A5, and 53A6 parameter sets are included. All standard
2860 building blocks are included and topologies can be built automatically
2863 The \gromosv{96} force field is a further development of the \gromosv{87} force field.
2864 It has improvements over the \gromosv{87} force field for proteins and small molecules.
2865 {\bf Note} that the sugar parameters present in 53A6 do correspond to those published in
2866 2004\cite{Oostenbrink2004}, which are different from those present in 45A4, which
2867 is not included in {\gromacs} at this time. The 45A4 parameter set corresponds to a later
2868 revision of these parameters.
2869 The \gromosv{96} force field is not, however, recommended for use with long alkanes and
2870 lipids. The \gromosv{96} force field differs from the \gromosv{87}
2871 force field in a few respects:
2873 \item the force field parameters
2874 \item the parameters for the bonded interactions are not linked to atom types
2875 \item a fourth power bond stretching potential (\ssecref{G96bond})
2876 \item an angle potential based on the cosine of the angle (\ssecref{G96angle})
2878 There are two differences in implementation between {\gromacs} and \gromosv{96}
2879 which can lead to slightly different results when simulating the same system
2882 \item in \gromosv{96} neighbor searching for solvents is performed on the
2883 first atom of the solvent molecule. This is not implemented in {\gromacs},
2884 but the difference with searching by centers of charge groups is very small
2885 \item the virial in \gromosv{96} is molecule-based. This is not implemented in
2886 {\gromacs}, which uses atomic virials
2888 The \gromosv{96} force field was parameterized with a Lennard-Jones cut-off
2889 of 1.4 nm, so be sure to use a Lennard-Jones cut-off ({\tt rvdw}) of at least 1.4.
2890 A larger cut-off is possible because the Lennard-Jones potential and forces
2891 are almost zero beyond 1.4 nm.
2893 \subsubsection{\gromosv{96} files\swapindexquiet{GROMOS96}{files}}
2894 {\gromacs} can read and write \gromosv{96} coordinate and trajectory files.
2895 These files should have the extension {\tt .g96}.
2896 Such a file can be a \gromosv{96} initial/final
2897 configuration file, a coordinate trajectory file, or a combination of both.
2898 The file is fixed format; all floats are written as 15.9, and as such, files can get huge.
2899 {\gromacs} supports the following data blocks in the given order:
2909 POSITION/POSITIONRED (mandatory)
2910 VELOCITY/VELOCITYRED (optional)
2915 See the \gromosv{96} manual~\cite{gromos96} for a complete description
2916 of the blocks. {\bf Note} that all {\gromacs} programs can read compressed
2917 (.Z) or gzipped (.gz) files.
2919 \subsection{OPLS/AA\index{OPLS/AA force field}}
2921 \subsection{AMBER\index{AMBER force field}}
2923 {\gromacs} provides native support for the following AMBER force fields:
2926 \item AMBER94~\cite{Cornell1995}
2927 \item AMBER96~\cite{Kollman1996}
2928 \item AMBER99~\cite{Wang2000}
2929 \item AMBER99SB~\cite{Hornak2006}
2930 \item AMBER99SB-ILDN~\cite{Lindorff2010}
2931 \item AMBER03~\cite{Duan2003}
2932 \item AMBERGS~\cite{Garcia2002}
2935 \subsection{CHARMM\index{CHARMM force field}}
2936 \label{subsec:charmmff}
2938 {\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.
2940 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.
2942 A port of the CHARMM36 force field for use with GROMACS is also available at \url{http://mackerell.umaryland.edu/charmm_ff.shtml#gromacs}.
2944 \subsection{Coarse-grained force fields}
2945 \index{force-field, coarse-grained}
2946 \label{sec:cg-forcefields}
2947 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.
2949 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:
2951 \item Conserving free energies
2953 \item Simplex method
2954 \item MARTINI force field (see next section)
2956 \item Conserving distributions (like the radial distribution function), so-called structure-based coarse-graining
2958 \item (iterative) Boltzmann inversion
2959 \item Inverse Monte Carlo
2961 \item Conversing forces
2963 \item Force matching
2967 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}.
2969 \subsection{MARTINI\index{Martini force field}}
2971 The MARTINI force field is a coarse-grain parameter set that allows for the construction
2972 of many systems, including proteins and membranes.
2974 \subsection{PLUM\index{PLUM force field}}
2976 The \normindex{PLUM force field}~\cite{bereau12} is an example of a solvent-free
2977 protein-membrane model for which the membrane was derived from structure-based
2978 coarse-graining~\cite{wang_jpcb10}. A {\gromacs} implementation can be found at
2979 \href{http://code.google.com/p/plumx/}{code.google.com/p/plumx}.
2981 % LocalWords: dihedrals centro ij dV dr LJ lj rcl jj Bertelot OPLS bh bham rf
2982 % LocalWords: coul defunits grompp crf vcrf fcrf Tironi Debye kgrf cgrf krf dx
2983 % LocalWords: PPPM der Waals erfc lr elstat chirality bstretch bondpot kT kJ
2984 % LocalWords: anharmonic morse mol betaij expminx SPC timestep fs FENE ijk kj
2985 % LocalWords: anglepot CHARMm UB ik rr substituents ijkl Ryckaert Bellemans rb
2986 % LocalWords: alkanes pdb gmx IUPAC IUB jkl cis rbdih crb kcal cubicspline xvg
2987 % LocalWords: topfile mdrun posres ar dihr lcllll NMR nmr lcllllll NOEs lclll
2988 % LocalWords: rav preprocessor ccccccccc ai aj fac disre mdp multi topol tpr
2989 % LocalWords: fc ravdisre nstdisreout dipolar lll ccc orientational MSD const
2990 % LocalWords: orire fitgrp nstorireout Drude intra Noskov et al fecalc coulrf
2991 % LocalWords: polarizabilities parameterized sigeps Ek sc softcore GROMOS NBF
2992 % LocalWords: hydrogens alkane vdwtype coulombtype mdpopt rlist rcoulomb rvdw
2993 % LocalWords: nstlist virial funcparm VdW jk jl fvsite fd vsites lcr vsitetop
2994 % LocalWords: vsite lclllll lcl parameterize parameterization enercorr avcsix
2995 % LocalWords: pcorr ecorr totalcoulomb dir fourierspacing ewald rtol Darden gz
2996 % LocalWords: FFT parallelization MPMD mpmd pme fft hoc Gromos gromos oxygens
2997 % LocalWords: virials POSITIONRED VELOCITYRED gzipped Charmm Larsson Bjelkmar
2998 % LocalWords: Cuendet CMAP nocmap dihedral Lennard covalent Verlet
2999 % LocalWords: Berthelot nm flexwat ferguson itp harmonicangle versa
3000 % LocalWords: harmonicbond atomtypes dihedraltypes equilibrated fdn
3001 % LocalWords: distancerestraint LINCS Coulombic ja jb il SPME ILDN
3002 % LocalWords: Hamiltonians atomtype AMBERGS rtp cmap graining VOTCA