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}
120 finally an geometric average for both parameters can be used (type 3):
123 \sigma_{ij} &=& \left({\sigma_{ii} \, \sigma_{jj}}\right)^{1/2} \\
124 \epsilon_{ij} &=& \left({\epsilon_{ii} \, \epsilon_{jj}}\right)^{1/2}
127 This last rule is used by the OPLS force field.
130 %\ifthenelse{\equal{\gmxlite}{1}}{}{
131 \subsection{\normindex{Buckingham potential}}
133 potential has a more flexible and realistic repulsion term
134 than the Lennard-Jones interaction, but is also more expensive to
135 compute. The potential form is:
137 V_{bh}(\rij) = A_{ij} \exp(-B_{ij} \rij) -
138 \frac{C_{ij}}{\rij^6}
141 \centerline{\includegraphics[width=8cm]{plots/f-bham}}
142 \caption {The Buckingham interaction.}
146 See also \figref{bham}. The force derived from this is:
148 \ve{F}_i(\rij) = \left[ A_{ij}B_{ij}\exp(-B_{ij} \rij) -
149 6\frac{C_{ij}}{\rij^7} \right] \rnorm
152 %} % Brace matches ifthenelse test for gmxlite
154 \subsection{Coulomb interaction}
156 \newcommand{\epsr}{\varepsilon_r}
157 \newcommand{\epsrf}{\varepsilon_{rf}}
158 The \normindex{Coulomb} interaction between two charge particles is given by:
160 V_c(\rij) = f \frac{q_i q_j}{\epsr \rij}
163 See also \figref{coul}, where $f = \frac{1}{4\pi \varepsilon_0} =
164 138.935\,485$ (see \chref{defunits})
167 \centerline{\includegraphics[width=8cm]{plots/vcrf}}
168 \caption[The Coulomb interaction with and without reaction field.]{The
169 Coulomb interaction (for particles with equal signed charge) with and
170 without reaction field. In the latter case $\epsr$ was 1, $\epsrf$ was 78,
171 and $r_c$ was 0.9 nm.
172 The dot-dashed line is the same as the dashed line, except for a constant.}
176 The force derived from this potential is:
178 \ve{F}_i(\rvij) = f \frac{q_i q_j}{\epsr\rij^2}\rnorm
181 In {\gromacs} the relative \swapindex{dielectric}{constant}
183 may be set in the in the input for {\tt grompp}.
185 %\ifthenelse{\equal{\gmxlite}{1}}{}{
186 \subsection{Coulomb interaction with \normindex{reaction field}}
188 The Coulomb interaction can be modified for homogeneous systems by
189 assuming a constant dielectric environment beyond the cut-off $r_c$
190 with a dielectric constant of {$\epsrf$}. The interaction then reads:
193 f \frac{q_i q_j}{\epsr\rij}\left[1+\frac{\epsrf-\epsr}{2\epsrf+\epsr}
194 \,\frac{\rij^3}{r_c^3}\right]
195 - f\frac{q_i q_j}{\epsr r_c}\,\frac{3\epsrf}{2\epsrf+\epsr}
198 in which the constant expression on the right makes the potential
199 zero at the cut-off $r_c$. For charged cut-off spheres this corresponds
200 to neutralization with a homogeneous background charge.
201 We can rewrite \eqnref{vcrf} for simplicity as
203 V_{crf} ~=~ f \frac{q_i q_j}{\epsr}\left[\frac{1}{\rij} + k_{rf}~ \rij^2 -c_{rf}\right]
207 k_{rf} &=& \frac{1}{r_c^3}\,\frac{\epsrf-\epsr}{(2\epsrf+\epsr)} \label{eqn:krf}\\
208 c_{rf} &=& \frac{1}{r_c}+k_{rf}\,r_c^2 ~=~ \frac{1}{r_c}\,\frac{3\epsrf}{(2\epsrf+\epsr)}
211 For large $\epsrf$ the $k_{rf}$ goes to $r_c^{-3}/2$,
212 while for $\epsrf$ = $\epsr$ the correction vanishes.
214 the modified interaction is plotted, and it is clear that the derivative
215 with respect to $\rij$ (= -force) goes to zero at the cut-off distance.
216 The force derived from this potential reads:
218 \ve{F}_i(\rvij) = f \frac{q_i q_j}{\epsr}\left[\frac{1}{\rij^2} - 2 k_{rf}\rij\right]\rnorm
221 The reaction-field correction should also be applied to all excluded
222 atoms pairs, including self pairs, in which case the normal Coulomb
223 term in \eqnsref{vcrf}{fcrf} is absent.
225 Tironi {\etal} have introduced a generalized reaction field in which
226 the dielectric continuum beyond the cut-off $r_c$ also has an ionic strength
227 $I$~\cite{Tironi95}. In this case we can rewrite the constants $k_{rf}$ and
228 $c_{rf}$ using the inverse Debye screening length $\kappa$:
231 \frac{2 I \,F^2}{\varepsilon_0 \epsrf RT}
232 = \frac{F^2}{\varepsilon_0 \epsrf RT}\sum_{i=1}^{K} c_i z_i^2 \\
233 k_{rf} &=& \frac{1}{r_c^3}\,
234 \frac{(\epsrf-\epsr)(1 + \kappa r_c) + \half\epsrf(\kappa r_c)^2}
235 {(2\epsrf + \epsr)(1 + \kappa r_c) + \epsrf(\kappa r_c)^2}
237 c_{rf} &=& \frac{1}{r_c}\,
238 \frac{3\epsrf(1 + \kappa r_c + \half(\kappa r_c)^2)}
239 {(2\epsrf+\epsr)(1 + \kappa r_c) + \epsrf(\kappa r_c)^2}
242 where $F$ is Faraday's constant, $R$ is the ideal gas constant, $T$
243 the absolute temperature, $c_i$ the molar concentration for species
244 $i$ and $z_i$ the charge number of species $i$ where we have $K$
245 different species. In the limit of zero ionic strength ($\kappa=0$)
246 \eqnsref{kgrf}{cgrf} reduce to the simple forms of \eqnsref{krf}{crf}
249 \subsection{Modified non-bonded interactions}
250 \label{sec:mod_nb_int}
251 In {\gromacs}, the non-bonded potentials can be
252 modified by a shift function. The purpose of this is to replace the
253 truncated forces by forces that are continuous and have continuous
254 derivatives at the \normindex{cut-off} radius. With such forces the
255 timestep integration produces much smaller errors and there are no
256 such complications as creating charges from dipoles by the truncation
257 procedure. In fact, by using shifted forces there is no need for
258 charge groups in the construction of neighbor lists. However, the
259 shift function produces a considerable modification of the Coulomb
260 potential. Unless the ``missing'' long-range potential is properly
261 calculated and added (through the use of PPPM, Ewald, or PME), the
262 effect of such modifications must be carefully evaluated. The
263 modification of the Lennard-Jones dispersion and repulsion is only
264 minor, but it does remove the noise caused by cut-off effects.
266 There is {\em no} fundamental difference between a switch function
267 (which multiplies the potential with a function) and a shift function
268 (which adds a function to the force or potential)~\cite{Spoel2006a}. The switch
269 function is a special case of the shift function, which we apply to
270 the {\em force function} $F(r)$, related to the electrostatic or
271 van der Waals force acting on particle $i$ by particle $j$ as:
273 \ve{F}_i = c F(r_{ij}) \frac{\rvij}{r_{ij}}
275 For pure Coulomb or Lennard-Jones interactions
276 $F(r)=F_\alpha(r)=r^{-(\alpha+1)}$.
277 The shifted force $F_s(r)$ can generally be written as:
281 F_s(r)~=&~F_\alpha(r) & r < r_1 \\
283 F_s(r)~=&~F_\alpha(r)+S(r) & r_1 \le r < r_c \\
284 F_s(r)~=&~0 & r_c \le r
287 When $r_1=0$ this is a traditional shift function, otherwise it acts as a
288 switch function. The corresponding shifted coulomb potential then reads:
290 V_s(r_{ij}) = f \Phi_s (r_{ij}) q_i q_j
292 where $\Phi(r)$ is the potential function
294 \Phi_s(r) = \int^{\infty}_r~F_s(x)\, dx
297 The {\gromacs} shift function should be smooth at the boundaries, therefore
298 the following boundary conditions are imposed on the shift function:
303 S(r_c) &=&-F_\alpha(r_c) \\
304 S'(r_c) &=&-F_\alpha'(r_c)
307 A 3$^{rd}$ degree polynomial of the form
309 S(r) = A(r-r_1)^2 + B(r-r_1)^3
311 fulfills these requirements. The constants A and B are given by the
312 boundary condition at $r_c$:
316 A &~=~& -\displaystyle
317 \frac{(\alpha+4)r_c~-~(\alpha+1)r_1} {r_c^{\alpha+2}~(r_c-r_1)^2} \\
318 B &~=~& \displaystyle
319 \frac{(\alpha+3)r_c~-~(\alpha+1)r_1}{r_c^{\alpha+2}~(r_c-r_1)^3}
322 Thus the total force function is:
324 F_s(r) = \frac{\alpha}{r^{\alpha+1}} + A(r-r_1)^2 + B(r-r_1)^3
326 and the potential function reads:
328 \Phi(r) = \frac{1}{r^\alpha} - \frac{A}{3} (r-r_1)^3 - \frac{B}{4} (r-r_1)^4 - C
332 C = \frac{1}{r_c^\alpha} - \frac{A}{3} (r_c-r_1)^3 - \frac{B}{4} (r_c-r_1)^4
335 When $r_1$ = 0, the modified Coulomb force function is
337 F_s(r) = \frac{1}{r^2} - \frac{5 r^2}{r_c^4} + \frac{4 r^3}{r_c^5}
339 which is identical to the {\em \index{parabolic force}}
340 function recommended to be used as a short-range function in
341 conjunction with a \swapindex{Poisson}{solver}
342 for the long-range part~\cite{Berendsen93a}.
343 The modified Coulomb potential function is:
345 \Phi(r) = \frac{1}{r} - \frac{5}{3r_c} + \frac{5r^3}{3r_c^4} - \frac{r^4}{r_c^5}
347 See also \figref{shift}.
350 \centerline{\includegraphics[angle=270,width=10cm]{plots/shiftf}}
351 \caption[The Coulomb Force, Shifted Force and Shift Function
352 $S(r)$,.]{The Coulomb Force, Shifted Force and Shift Function $S(r)$,
353 using r$_1$ = 2 and r$_c$ = 4.}
357 \subsection{Modified short-range interactions with Ewald summation}
358 When Ewald summation\index{Ewald sum} or \seeindex{particle-mesh
359 Ewald}{PME}\index{Ewald, particle-mesh} is used to calculate the
360 long-range interactions, the
361 short-range Coulomb potential must also be modified, similar to the
362 switch function above. In this case the short range potential is given
365 V(r) = f \frac{\mbox{erfc}(\beta r_{ij})}{r_{ij}} q_i q_j,
367 where $\beta$ is a parameter that determines the relative weight
368 between the direct space sum and the reciprocal space sum and erfc$(x)$ is
369 the complementary error function. For further
370 details on long-range electrostatics, see \secref{lr_elstat}.
371 %} % Brace matches ifthenelse test for gmxlite
374 \section{Bonded interactions}
375 Bonded interactions are based on a fixed list of atoms. They are not
376 exclusively pair interactions, but include 3- and 4-body interactions
377 as well. There are {\em bond stretching} (2-body), {\em bond angle}
378 (3-body), and {\em dihedral angle} (4-body) interactions. A special
379 type of dihedral interaction (called {\em improper dihedral}) is used
380 to force atoms to remain in a plane or to prevent transition to a
381 configuration of opposite chirality (a mirror image).
383 \subsection{Bond stretching}
385 \subsubsection{Harmonic potential}
386 \label{subsec:harmonicbond}
387 The \swapindex{bond}{stretching} between two covalently bonded atoms
388 $i$ and $j$ is represented by a harmonic potential:
391 \centerline{\raisebox{4cm}{\includegraphics[angle=270,width=5cm]{plots/bstretch}}\includegraphics[width=7cm]{plots/f-bond}}
392 \caption[Bond stretching.]{Principle of bond stretching (left), and the bond
393 stretching potential (right).}
394 \label{fig:bstretch1}
398 V_b~(\rij) = \half k^b_{ij}(\rij-b_{ij})^2
400 See also \figref{bstretch1}, with the force given by:
402 \ve{F}_i(\rvij) = k^b_{ij}(\rij-b_{ij}) \rnorm
405 %\ifthenelse{\equal{\gmxlite}{1}}{}{
406 \subsubsection{Fourth power potential}
407 \label{subsec:G96bond}
408 In the \gromosv{96} force field~\cite{gromos96}, the covalent bond potential
409 is, for reasons of computational efficiency, written as:
411 V_b~(\rij) = \frac{1}{4}k^b_{ij}\left(\rij^2-b_{ij}^2\right)^2
413 The corresponding force is:
415 \ve{F}_i(\rvij) = k^b_{ij}(\rij^2-b_{ij}^2)~\rvij
417 The force constants for this form of the potential are related to the usual
418 harmonic force constant $k^{b,harm}$ (\secref{bondpot}) as
420 2 k^b b_{ij}^2 = k^{b,harm}
422 The force constants are mostly derived from the harmonic ones used in
423 \gromosv{87}~\cite{biomos}. Although this form is computationally more
425 (because no square root has to be evaluated), it is conceptually more
426 complex. One particular disadvantage is that since the form is not harmonic,
427 the average energy of a single bond is not equal to $\half kT$ as it is for
428 the normal harmonic potential.
430 \subsection{\normindex{Morse potential} bond stretching}
431 \label{subsec:Morsebond}
432 %\author{F.P.X. Everdij}
434 For some systems that require an anharmonic bond stretching potential,
435 the Morse potential~\cite{Morse29}
436 between two atoms {\it i} and {\it j} is available
437 in {\gromacs}. This potential differs from the harmonic potential in
438 that it has an asymmetric potential well and a zero force at infinite
439 distance. The functional form is:
441 \displaystyle V_{morse} (r_{ij}) = D_{ij} [1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))]^2,
443 See also \figref{morse}, and the corresponding force is:
446 \displaystyle {\bf F}_{morse} ({\bf r}_{ij})&=&2 D_{ij} \beta_{ij} r_{ij} \exp(-\beta_{ij}(r_{ij}-b_{ij})) * \\
447 \displaystyle \: & \: &[1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))] \frac{\displaystyle {\bf r}_{ij}}{\displaystyle r_{ij}},
450 where \( \displaystyle D_{ij} \) is the depth of the well in kJ/mol,
451 \( \displaystyle \beta_{ij} \) defines the steepness of the well (in
452 nm\(^{-1} \)), and \( \displaystyle b_{ij} \) is the equilibrium
453 distance in nm. The steepness parameter \( \displaystyle \beta_{ij}
454 \) can be expressed in terms of the reduced mass of the atoms {\it i}
455 and {\it j}, the fundamental vibration frequency \( \displaystyle
456 \omega_{ij} \) and the well depth \( \displaystyle D_{ij} \):
458 \displaystyle \beta_{ij}= \omega_{ij} \sqrt{\frac{\mu_{ij}}{2 D_{ij}}}
460 and because \( \displaystyle \omega = \sqrt{k/\mu} \), one can rewrite \( \displaystyle \beta_{ij} \) in terms of the harmonic force constant \( \displaystyle k_{ij} \):
462 \displaystyle \beta_{ij}= \sqrt{\frac{k_{ij}}{2 D_{ij}}}
465 For small deviations \( \displaystyle (r_{ij}-b_{ij}) \), one can
466 approximate the \( \displaystyle \exp \)-term to first-order using a
469 \displaystyle \exp(-x) \approx 1-x
472 and substituting \eqnref{betaij} and \eqnref{expminx} in the functional form:
475 \displaystyle V_{morse} (r_{ij})&=&D_{ij} [1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))]^2\\
476 \displaystyle \:&=&D_{ij} [1 - (1 -\sqrt{\frac{k_{ij}}{2 D_{ij}}}(r_{ij}-b_{ij}))]^2\\
477 \displaystyle \:&=&\frac{1}{2} k_{ij} (r_{ij}-b_{ij}))^2
480 we recover the harmonic bond stretching potential.
483 \centerline{\includegraphics[width=7cm]{plots/f-morse}}
484 \caption{The Morse potential well, with bond length 0.15 nm.}
488 \subsection{Cubic bond stretching potential}
489 \label{subsec:cubicbond}
490 Another anharmonic bond stretching potential that is slightly simpler
491 than the Morse potential adds a cubic term in the distance to the
492 simple harmonic form:
494 V_b~(\rij) = k^b_{ij}(\rij-b_{ij})^2 + k^b_{ij}k^{cub}_{ij}(\rij-b_{ij})^3
496 A flexible \normindex{water} model (based on
497 the SPC water model~\cite{Berendsen81}) including
498 a cubic bond stretching potential for the O-H bond
499 was developed by Ferguson~\cite{Ferguson95}. This model was found
500 to yield a reasonable infrared spectrum. The Ferguson water model is
501 available in the {\gromacs} library ({\tt flexwat-ferguson.itp}).
502 It should be noted that the potential is asymmetric: overstretching leads to
503 infinitely low energies. The \swapindex{integration}{timestep} is therefore
506 The force corresponding to this potential is:
508 \ve{F}_i(\rvij) = 2k^b_{ij}(\rij-b_{ij})~\rnorm + 3k^b_{ij}k^{cub}_{ij}(\rij-b_{ij})^2~\rnorm
511 \subsection{FENE bond stretching potential\index{FENE potential}}
512 \label{subsec:FENEbond}
513 In coarse-grained polymer simulations the beads are often connected
514 by a FENE (finitely extensible nonlinear elastic) potential~\cite{Warner72}:
516 V_{\mbox{\small FENE}}(\rij) =
517 -\half k^b_{ij} b^2_{ij} \log\left(1 - \frac{\rij^2}{b^2_{ij}}\right)
519 The potential looks complicated, but the expression for the force is simpler:
521 F_{\mbox{\small FENE}}(\rvij) =
522 -k^b_{ij} \left(1 - \frac{\rij^2}{b^2_{ij}}\right)^{-1} \rvij
524 At short distances the potential asymptotically goes to a harmonic
525 potential with force constant $k^b$, while it diverges at distance $b$.
526 %} % Brace matches ifthenelse test for gmxlite
528 \subsection{Harmonic angle potential}
529 \label{subsec:harmonicangle}
530 \newcommand{\tijk}{\theta_{ijk}}
531 The bond-\swapindex{angle}{vibration} between a triplet of atoms $i$ - $j$ - $k$
532 is also represented by a harmonic potential on the angle $\tijk$
535 \centerline{\raisebox{4cm}{\includegraphics[angle=270,width=5cm]{plots/angle}}\includegraphics[width=7cm]{plots/f-angle}}
536 \caption[Angle vibration.]{Principle of angle vibration (left) and the
537 bond angle potential (right).}
542 V_a(\tijk) = \half k^{\theta}_{ijk}(\tijk-\tijk^0)^2
544 As the bond-angle vibration is represented by a harmonic potential, the
545 form is the same as the bond stretching (\figref{bstretch1}).
547 The force equations are given by the chain rule:
550 \Fvi ~=~ -\displaystyle\frac{d V_a(\tijk)}{d \rvi} \\
551 \Fvk ~=~ -\displaystyle\frac{d V_a(\tijk)}{d \rvk} \\
554 ~ \mbox{ ~ where ~ } ~
555 \tijk = \arccos \frac{(\rvij \cdot \ve{r}_{kj})}{r_{ij}r_{kj}}
557 The numbering $i,j,k$ is in sequence of covalently bonded atoms. Atom
558 $j$ is in the middle; atoms $i$ and $k$ are at the ends (see \figref{angle}).
559 {\bf Note} that in the input in topology files, angles are given in degrees and
560 force constants in kJ/mol/rad$^2$.
562 %\ifthenelse{\equal{\gmxlite}{1}}{}{
563 \subsection{Cosine based angle potential}
564 \label{subsec:G96angle}
565 In the \gromosv{96} force field a simplified function is used to represent angle
568 V_a(\tijk) = \half k^{\theta}_{ijk}\left(\cos(\tijk) - \cos(\tijk^0)\right)^2
573 \cos(\tijk) = \frac{\rvij\cdot\ve{r}_{kj}}{\rij r_{kj}}
575 The corresponding force can be derived by partial differentiation with respect
576 to the atomic positions. The force constants in this function are related
577 to the force constants in the harmonic form $k^{\theta,harm}$
578 (\ssecref{harmonicangle}) by:
580 k^{\theta} \sin^2(\tijk^0) = k^{\theta,harm}
582 In the \gromosv{96} manual there is a much more complicated conversion formula
583 which is temperature dependent. The formulas are equivalent at 0 K
584 and the differences at 300 K are on the order of 0.1 to 0.2\%.
585 {\bf Note} that in the input in topology files, angles are given in degrees and
586 force constants in kJ/mol.
588 \subsection{Restricted bending potential}
590 The restricted bending (ReB) potential~\cite{MonicaGoga2013} prevents the bending angle $\theta$
591 from reaching the $180^{\circ}$ value. In this way, the numerical instabilities
592 due to the calculation of the torsion angle and potential are eliminated when
593 performing coarse-grained molecular dynamics simulations.
595 To systematically hinder the bending angles from reaching the $180^{\circ}$ value,
596 the bending potential \ref{eq:G96angle} is divided by a $\sin^2\theta$ factor:
599 V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta} \frac{(\cos\theta_i - \cos\theta_0)^2}{\sin^2\theta_i}.
603 Figure ~\figref{ReB} shows the comparison between the ReB potential, \ref{eq:ReB},
604 and the standard one \ref{eq:G96angle}.
607 \centerline{\includegraphics[width=10cm]{plots/fig-02}}
609 \caption{Bending angle potentials: cosine harmonic (solid black line), angle harmonic
610 (dashed black line) and restricted bending (red) with the same bending constant
611 $k_{\theta}=85$ kJ mol$^{-1}$ and equilibrium angle $\theta_0=130^{\circ}$.
612 The orange line represents the sum of a cosine harmonic ($k =50$ kJ mol$^{-1}$)
613 with a restricted bending ($k =25$ kJ mol$^{-1}$) potential, both with $\theta_0=130^{\circ}$.}
617 The wall of the ReB potential is very repulsive in the region close to $180^{\circ}$ and,
618 as a result, the bending angles are kept within a safe interval, far from instabilities.
619 The power $2$ of $\sin\theta_i$ in the denominator has been chosen to guarantee this behavior
620 and allows an elegant differentiation:
623 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}}.
627 Due to its construction, the restricted bending potential cannot be used for equilibrium
628 $\theta_0$ values too close to $0^{\circ}$ or $180^{\circ}$ (from experience, at least $10^{\circ}$
629 difference is recommended). It is very important that, in the starting configuration,
630 all the bending angles have to be in the safe interval to avoid initial instabilities.
631 This bending potential can be used in combination with any form of torsion potential.
632 It will always prevent three consecutive particles from becoming collinear and,
633 as a result, any torsion potential will remain free of singularities.
634 It can be also added to a standard bending potential to affect the angle around $180^{\circ}$,
635 but to keep its original form around the minimum (see the orange curve in \figref{ReB}).
638 \subsection{Urey-Bradley potential}
639 \label{subsec:Urey-Bradley}
640 The \swapindex{Urey-Bradley bond-angle}{vibration} between a triplet
641 of atoms $i$ - $j$ - $k$ is represented by a harmonic potential on the
642 angle $\tijk$ and a harmonic correction term on the distance between
643 the atoms $i$ and $k$. Although this can be easily written as a simple
644 sum of two terms, it is convenient to have it as a single entry in the
645 topology file and in the output as a separate energy term. It is used mainly
646 in the CHARMm force field~\cite{BBrooks83}. The energy is given by:
649 V_a(\tijk) = \half k^{\theta}_{ijk}(\tijk-\tijk^0)^2 + \half k^{UB}_{ijk}(r_{ik}-r_{ik}^0)^2
652 The force equations can be deduced from sections~\ssecref{harmonicbond}
653 and~\ssecref{harmonicangle}.
655 \subsection{Bond-Bond cross term}
656 \label{subsec:bondbondcross}
657 The bond-bond cross term for three particles $i, j, k$ forming bonds
658 $i-j$ and $k-j$ is given by~\cite{Lawrence2003b}:
660 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)
663 where $k_{rr'}$ is the force constant, and $r_{1e}$ and $r_{2e}$ are the
664 equilibrium bond lengths of the $i-j$ and $k-j$ bonds respectively. The force
665 associated with this potential on particle $i$ is:
667 \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|}
669 The force on atom $k$ can be obtained by swapping $i$ and $k$ in the above
670 equation. Finally, the force on atom $j$ follows from the fact that the sum
671 of internal forces should be zero: $\ve{F}_j = -\ve{F}_i-\ve{F}_k$.
673 \subsection{Bond-Angle cross term}
674 \label{subsec:bondanglecross}
675 The bond-angle cross term for three particles $i, j, k$ forming bonds
676 $i-j$ and $k-j$ is given by~\cite{Lawrence2003b}:
678 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)
680 where $k_{r\theta}$ is the force constant, $r_{3e}$ is the $i-k$ distance,
681 and the other constants are the same as in Equation~\ref{crossbb}. The force
682 associated with the potential on atom $i$ is:
684 \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|} \\
685 + \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]
688 \subsection{Quartic angle potential}
689 \label{subsec:quarticangle}
690 For special purposes there is an angle potential
691 that uses a fourth order polynomial:
693 V_q(\tijk) ~=~ \sum_{n=0}^5 C_n (\tijk-\tijk^0)^n
695 %} % Brace matches ifthenelse test for gmxlite
697 %% new commands %%%%%%%%%%%%%%%%%%%%%%
698 \newcommand{\rvkj}{{\bf r}_{kj}}
699 \newcommand{\rkj}{r_{kj}}
700 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
702 \subsection{Improper dihedrals\swapindexquiet{improper}{dihedral}}
704 Improper dihedrals are meant to keep \swapindex{planar}{group}s ({\eg}
705 aromatic rings) planar, or to prevent molecules from flipping over to their
706 \normindex{mirror image}s, see \figref{imp}.
709 \centerline{\includegraphics[angle=270,width=4cm]{plots/ring-imp}\hspace{1cm}
710 \includegraphics[angle=270,width=3cm]{plots/subst-im}\hspace{1cm}\includegraphics[angle=270,width=3cm]{plots/tetra-im}}
711 \caption[Improper dihedral angles.]{Principle of improper
712 dihedral angles. Out of plane bending for rings (left), substituents
713 of rings (middle), out of tetrahedral (right). The improper dihedral
714 angle $\xi$ is defined as the angle between planes (i,j,k) and (j,k,l)
719 \subsubsection{Improper dihedrals: harmonic type}
720 \label{subsec:harmonicimproperdihedral}
721 The simplest improper dihedral potential is a harmonic potential; it is plotted in
724 V_{id}(\xi_{ijkl}) = \half k_{\xi}(\xi_{ijkl}-\xi_0)^2
726 Since the potential is harmonic it is discontinuous,
727 but since the discontinuity is chosen at 180$^\circ$ distance from $\xi_0$
728 this will never cause problems.
729 {\bf Note} that in the input in topology files, angles are given in degrees and
730 force constants in kJ/mol/rad$^2$.
733 \centerline{\includegraphics[width=8cm]{plots/f-imps}}
734 \caption{Improper dihedral potential.}
738 \subsubsection{Improper dihedrals: periodic type}
739 \label{subsec:periodicimproperdihedral}
740 This potential is identical to the periodic proper dihedral (see below).
741 There is a separate dihedral type for this (type 4) only to be able
742 to distinguish improper from proper dihedrals in the parameter section
745 \subsection{Proper dihedrals\swapindexquiet{proper}{dihedral}}
746 For the normal \normindex{dihedral} interaction there is a choice of
747 either the {\gromos} periodic function or a function based on
748 expansion in powers of $\cos \phi$ (the so-called Ryckaert-Bellemans
749 potential). This choice has consequences for the inclusion of special
750 interactions between the first and the fourth atom of the dihedral
751 quadruple. With the periodic {\gromos} potential a special 1-4
752 LJ-interaction must be included; with the Ryckaert-Bellemans potential
753 {\em for alkanes} the \normindex{1-4 interaction}s must be excluded
754 from the non-bonded list. {\bf Note:} Ryckaert-Bellemans potentials
755 are also used in {\eg} the OPLS force field in combination with 1-4
756 interactions. You should therefore not modify topologies generated by
757 {\tt \normindex{pdb2gmx}} in this case.
759 \subsubsection{Proper dihedrals: periodic type}
760 \label{subsec:properdihedral}
761 Proper dihedral angles are defined according to the IUPAC/IUB
762 convention, where $\phi$ is the angle between the $ijk$ and the $jkl$
763 planes, with {\bf zero} corresponding to the {\em cis} configuration
764 ($i$ and $l$ on the same side). There are two dihedral function types
765 in {\gromacs} topology files. There is the standard type 1 which behaves
766 like any other bonded interactions. For certain force fields, type 9
767 is useful. Type 9 allows multiple potential functions to be applied
768 automatically to a single dihedral in the {\tt [ dihedral ]} section
769 when multiple parameters are defined for the same atomtypes
770 in the {\tt [ dihedraltypes ]} section.
773 \centerline{\raisebox{4.5cm}{\includegraphics[angle=270,width=5cm]{plots/dih}}\includegraphics[width=7cm]{plots/f-dih}}
774 \caption[Proper dihedral angle.]{Principle of proper dihedral angle
775 (left, in {\em trans} form) and the dihedral angle potential (right).}
779 V_d(\phi_{ijkl}) = k_{\phi}(1 + \cos(n \phi - \phi_s))
782 %\ifthenelse{\equal{\gmxlite}{1}}{}{
783 \subsubsection{Proper dihedrals: Ryckaert-Bellemans function}
784 \label{subsec:RBdihedral}
785 For alkanes, the following proper dihedral potential is often used
786 (see \figref{rbdih}):
788 V_{rb}(\phi_{ijkl}) = \sum_{n=0}^5 C_n( \cos(\psi ))^n,
790 where $\psi = \phi - 180^\circ$. \\
791 {\bf Note:} A conversion from one convention to another can be achieved by
792 multiplying every coefficient \( \displaystyle C_n \)
793 by \( \displaystyle (-1)^n \).
795 An example of constants for $C$ is given in \tabref{crb}.
799 \begin{tabular}{|lr|lr|lr|}
801 $C_0$ & 9.28 & $C_2$ & -13.12 & $C_4$ & 26.24 \\
802 $C_1$ & 12.16 & $C_3$ & -3.06 & $C_5$ & -31.5 \\
806 \caption{Constants for Ryckaert-Bellemans potential (kJ mol$^{-1}$).}
811 \centerline{\includegraphics[width=8cm]{plots/f-rbs}}
812 \caption{Ryckaert-Bellemans dihedral potential.}
816 ({\bf Note:} The use of this potential implies exclusion of LJ interactions
817 between the first and the last atom of the dihedral, and $\psi$ is defined
818 according to the ``polymer convention'' ($\psi_{trans}=0$).)
820 The RB dihedral function can also be used to include Fourier dihedrals
823 V_{rb} (\phi_{ijkl}) ~=~ \frac{1}{2} \left[F_1(1+\cos(\phi)) + F_2(
824 1-\cos(2\phi)) + F_3(1+\cos(3\phi)) + F_4(1-\cos(4\phi))\right]
826 Because of the equalities \( \cos(2\phi) = 2\cos^2(\phi) - 1 \),
827 \( \cos(3\phi) = 4\cos^3(\phi) - 3\cos(\phi) \) and
828 \( \cos(4\phi) = 8\cos^4(\phi) - 8\cos^2(\phi) + 1 \)
829 one can translate the OPLS parameters to
830 Ryckaert-Bellemans parameters as follows:
834 \displaystyle C_0&=&F_2 + \frac{1}{2} (F_1 + F_3)\\
835 \displaystyle C_1&=&\frac{1}{2} (- F_1 + 3 \, F_3)\\
836 \displaystyle C_2&=& -F_2 + 4 \, F_4\\
837 \displaystyle C_3&=&-2 \, F_3\\
838 \displaystyle C_4&=&-4 \, F_4\\
839 \displaystyle C_5&=&0
842 with OPLS parameters in protein convention and RB parameters in
843 polymer convention (this yields a minus sign for the odd powers of
845 \noindent{\bf Note:} Mind the conversion from {\bf kcal mol$^{-1}$} for
846 literature OPLS and RB parameters to {\bf kJ mol$^{-1}$} in {\gromacs}.\\
847 %} % Brace matches ifthenelse test for gmxlite
849 \subsubsection{Proper dihedrals: Fourier function}
850 \label{subsec:Fourierdihedral}
851 The OPLS potential function is given as the first three
852 or four~\cite{Jorgensen2005a} cosine terms of a Fourier series.
853 In {\gromacs} the four term function is implemented:
855 V_{F} (\phi_{ijkl}) ~=~ \frac{1}{2} \left[C_1(1+\cos(\phi)) + C_2(
856 1-\cos(2\phi)) + C_3(1+\cos(3\phi)) + C_4(1+\cos(4\phi))\right],
858 %\ifthenelse{\equal{\gmxlite}{1}}{}{
859 Internally, {\gromacs}
860 uses the Ryckaert-Bellemans code
861 to compute Fourier dihedrals (see above), because this is more efficient.\\
862 \noindent{\bf Note:} Mind the conversion from {\emph kcal mol$^{-1}$} for
863 literature OPLS parameters to {\bf kJ mol$^{-1}$} in {\gromacs}.\\
865 \subsubsection{Proper dihedrals: Restricted torsion potential}
867 In a manner very similar to the restricted bending potential (see \ref{subsec:ReB}),
868 a restricted torsion/dihedral potential is introduced:
871 V_{\rm ReT}(\phi_i) = \frac{1}{2} k_{\phi} \frac{(\cos\phi_i - \cos\phi_0)^2}{\sin^2\phi_i}
875 with the advantages of being a function of $\cos\phi$ (no problems taking the derivative of $\sin\phi$)
876 and of keeping the torsion angle at only one minimum value. In this case, the factor $\sin^2\phi$ does
877 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.
878 For this reason, all the dihedral angles of the starting configuration should have their values in the
879 desired angles interval and the the equilibrium $\phi_0$ value should not be too close to the interval limits
880 (as for the restricted bending potential, described in \ref{subsec:ReB}, at least $10^{\circ}$ difference is recommended).
882 \subsubsection{Proper dihedrals: Combined bending-torsion potential}
884 When the four particles forming the dihedral angle become collinear (this situation will never happen in
885 atomistic simulations, but it can occur in coarse-grained simulations) the calculation of the
886 torsion angle and potential leads to numerical instabilities.
887 One way to avoid this is to use the restricted bending potential (see \ref{subsec:ReB})
888 that prevents the dihedral
889 from reaching the $180^{\circ}$ value.
891 Another way is to disregard any effects of the dihedral becoming ill-defined,
892 keeping the dihedral force and potential calculation continuous in entire angle range
893 by coupling the torsion potential (in a cosine form) with the bending potentials of the
894 adjacent bending angles in a unique expression:
897 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}.
901 This combined bending-torsion (CBT) potential has been proposed by~\cite{BulacuGiessen2005}
902 for polymer melt simulations and is extensively described in~\cite{MonicaGoga2013}.
904 This potential has two main advantages:
907 it does not only depend on the dihedral angle $\phi_i$ (between the $i-2$, $i-1$, $i$ and $i+1$ beads)
908 but also on the bending angles $\theta_{i-1}$ and $\theta_i$ defined from three adjacent beads
909 ($i-2$, $i-1$ and $i$, and $i-1$, $i$ and $i+1$, respectively).
910 The two $\sin^3\theta$ pre-factors, tentatively suggested by~\cite{ScottScheragator1966} and theoretically
911 discussed by~\cite{PaulingBond}, cancel the torsion potential and force when either of the two bending angles
912 approaches the value of $180^\circ$.
914 its dependence on $\phi_i$ is expressed through a polynomial in $\cos\phi_i$ that avoids the singularities in
915 $\phi=0^\circ$ or $180^\circ$ in calculating the torsional force.
918 These two properties make the CBT potential well-behaved for MD simulations with weak constraints
919 on the bending angles or even for steered / non-equilibrium MD in which the bending and torsion angles suffer major
921 When using the CBT potential, the bending potentials for the adjacent $\theta_{i-1}$ and $\theta_i$ may have any form.
922 It is also possible to leave out the two angle bending terms ($\theta_{i-1}$ and $\theta_{i}$) completely.
923 \figref{CBT} illustrates the difference between a torsion potential with and without the $\sin^{3}\theta$ factors
924 (blue and gray curves, respectively).
927 \centerline{\includegraphics[width=10cm]{plots/fig-04}}
928 \caption{Blue: surface plot of the combined bending-torsion potential
929 (\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$)
930 when, for simplicity, the bending angles behave the same ($\theta_1=\theta_2=\theta$).
931 Gray: the same torsion potential without the $\sin^{3}\theta$ terms (Ryckaert-Bellemans type).
932 $\phi$ is the dihedral angle.}
936 Additionally, the derivative of $V_{CBT}$ with respect to the Cartesian variables is straightforward:
939 \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}} +
940 \frac{\partial V_{\rm CBT}}{\partial \theta_{i }} \frac{\partial \theta_{i }}{\partial \vec r_{l}} +
941 \frac{\partial V_{\rm CBT}}{\partial \phi_{i }} \frac{\partial \phi_{i }}{\partial \vec r_{l}}
945 The CBT is based on a cosine form without multiplicity, so it can only be symmetrical around $0^{\circ}$.
946 To obtain an asymmetrical dihedral angle distribution (e.g. only one maximum in [$-180^{\circ}$:$180^{\circ}$] interval),
947 a standard torsion potential such as harmonic angle or periodic cosine potentials should be used instead of a CBT potential.
948 However, these two forms have the inconveniences of the force derivation ($1/\sin\phi$) and of the alignment of beads
949 ($\theta_i$ or $\theta_{i-1} = 0^{\circ}, 180^{\circ}$).
950 Coupling such non-$\cos\phi$ potentials with $\sin^3\theta$ factors does not improve simulation stability since there are
951 cases in which $\theta$ and $\phi$ are simultaneously $180^{\circ}$. The integration at this step would be possible
952 (due to the cancelling of the torsion potential) but the next step would be singular
953 ($\theta$ is not $180^{\circ}$ and $\phi$ is very close to $180^{\circ}$).
955 %\ifthenelse{\equal{\gmxlite}{1}}{}{
956 \subsection{Tabulated bonded interaction functions\index{tabulated bonded interaction function}}
957 \label{subsec:tabulatedinteraction}
958 For full flexibility, any functional shape can be used for
959 bonds, angles and dihedrals through user-supplied tabulated functions.
960 The functional shapes are:
962 V_b(r_{ij}) &=& k \, f^b_n(r_{ij}) \\
963 V_a(\tijk) &=& k \, f^a_n(\tijk) \\
964 V_d(\phi_{ijkl}) &=& k \, f^d_n(\phi_{ijkl})
966 where $k$ is a force constant in units of energy
967 and $f$ is a cubic spline function; for details see \ssecref{cubicspline}.
968 For each interaction, the force constant $k$ and the table number $n$
969 are specified in the topology.
970 There are two different types of bonds, one that generates exclusions (type 8)
971 and one that does not (type 9).
972 For details see \tabref{topfile2}.
973 The table files are supplied to the {\tt mdrun} program.
974 After the table file name an underscore, the letter ``b'' for bonds,
975 ``a'' for angles or ``d'' for dihedrals and the table number are appended.
976 For example, for a bond with $n=0$ (and using the default table file name)
977 the table is read from the file {\tt table_b0.xvg}. Multiple tables can be
978 supplied simply by using different values of $n$, and are applied to the appropriate
979 bonds, as specified in the topology (\tabref{topfile2}).
980 The format for the table files is three columns with $x$, $f(x)$, $-f'(x)$,
981 where $x$ should be uniformly-spaced. Requirements for entries in the topology
982 are given in~\tabref{topfile2}.
983 The setup of the tables is as follows:
985 $x$ is the distance in nm. For distances beyond the table length,
986 {\tt mdrun} will quit with an error message.
988 $x$ is the angle in degrees. The table should go from
989 0 up to and including 180 degrees; the derivative is taken in degrees.
991 $x$ is the dihedral angle in degrees. The table should go from
992 -180 up to and including 180 degrees;
993 the IUPAC/IUB convention is used, {\ie} zero is cis,
994 the derivative is taken in degrees.
995 %} % Brace matches ifthenelse test for gmxlite
998 Special potentials are used for imposing restraints on the motion of
999 the system, either to avoid disastrous deviations, or to include
1000 knowledge from experimental data. In either case they are not really
1001 part of the force field and the reliability of the parameters is not
1002 important. The potential forms, as implemented in {\gromacs}, are
1003 mentioned just for the sake of completeness. Restraints and constraints
1004 refer to quite different algorithms in {\gromacs}.
1006 \subsection{Position restraints\swapindexquiet{position}{restraint}}
1007 \label{subsec:positionrestraint}
1008 These are used to restrain particles to fixed reference positions
1009 $\ve{R}_i$. They can be used during equilibration in order to avoid
1010 drastic rearrangements of critical parts ({\eg} to restrain motion
1011 in a protein that is subjected to large solvent forces when the
1012 solvent is not yet equilibrated). Another application is the
1013 restraining of particles in a shell around a region that is simulated
1014 in detail, while the shell is only approximated because it lacks
1015 proper interaction from missing particles outside the
1016 shell. Restraining will then maintain the integrity of the inner
1017 part. For spherical shells, it is a wise procedure to make the force
1018 constant depend on the radius, increasing from zero at the inner
1019 boundary to a large value at the outer boundary. This feature has
1020 not, however, been implemented in {\gromacs}.
1021 \newcommand{\unitv}[1]{\hat{\bf #1}}
1022 \newcommand{\halfje}[1]{\frac{#1}{2}}
1024 The following form is used:
1026 V_{pr}(\ve{r}_i) = \halfje{1}k_{pr}|\rvi-\ve{R}_i|^2
1028 The potential is plotted in \figref{positionrestraint}.
1031 \centerline{\includegraphics[width=8cm]{plots/f-pr}}
1032 \caption{Position restraint potential.}
1033 \label{fig:positionrestraint}
1036 The potential form can be rewritten without loss of generality as:
1038 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]
1044 F_i^x &=& -k_{pr}^x~(x_i - X_i) \\
1045 F_i^y &=& -k_{pr}^y~(y_i - Y_i) \\
1046 F_i^z &=& -k_{pr}^z~(z_i - Z_i)
1049 Using three different force constants the position
1050 restraints can be turned on or off
1051 in each spatial dimension; this means that atoms can be harmonically
1052 restrained to a plane or a line.
1053 Position restraints are applied to a special fixed list of atoms. Such a
1054 list is usually generated by the {\tt \normindex{pdb2gmx}} program.
1056 \subsection{\swapindex{Flat-bottomed}{position restraint}s}
1057 \label{subsec:fbpositionrestraint}
1058 Flat-bottomed position restraints can be used to restrain particles to
1059 part of the simulation volume. No force acts on the restrained
1060 particle within the flat-bottomed region of the potential, however a
1061 harmonic force acts to move the particle to the flat-bottomed region
1062 if it is outside it. It is possible to apply normal and
1063 flat-bottomed position restraints on the same particle (however, only
1064 with the same reference position $\ve{R}_i$). The following general potential
1065 is used (Figure~\ref{fig:fbposres}A):
1067 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}],
1069 where $\ve{R}_i$ is the reference position, $r_\mathrm{fb}$ is the distance
1070 from the center with a flat potential, $k_\mathrm{fb}$ the force constant, and $H$ is the Heaviside step
1071 function. The distance $d_g(\ve{r}_i;\ve{R}_i)$ from the reference
1072 position depends on the geometry $g$ of the flat-bottomed potential.
1075 \centerline{\includegraphics[width=10cm]{plots/fbposres}}
1076 \caption{Flat-bottomed position restraint potential. (A) Not
1077 inverted, (B) inverted.}
1078 \label{fig:fbposres}
1081 The following geometries for the flat-bottomed potential are supported:\newline
1082 {\bfseries Sphere} ($g =1$): The particle is kept in a sphere of given
1083 radius. The force acts towards the center of the sphere. The following distance calculation is used:
1085 d_g(\ve{r}_i;\ve{R}_i) = |\ve{r}_i-\ve{R}_i|
1087 {\bfseries Cylinder} ($g=2$): The particle is kept in a cylinder of given radius
1088 parallel to the $z$-axis. The force from the flat-bottomed potential acts
1089 towards the axis of the cylinder. The $z$-component of the force is zero.
1091 d_g(\ve{r}_i;\ve{R}_i) = \sqrt{ (x_i-X_i)^2 + (y_i - Y_i)^2 }
1093 {\bfseries Layer} ($g=3,4,5$): The particle is kept in a layer defined by the
1094 thickness and the normal of the layer. The layer normal can be parallel to the $x$, $y$, or
1095 $z$-axis. The force acts parallel to the layer normal.\\
1097 d_g(\ve{r}_i;\ve{R}_i) = |x_i-X_i|, \;\;\;\mbox{or}\;\;\;
1098 d_g(\ve{r}_i;\ve{R}_i) = |y_i-Y_i|, \;\;\;\mbox{or}\;\;\;
1099 d_g(\ve{r}_i;\ve{R}_i) = |z_i-Z_i|.
1102 It is possible to apply multiple independent flat-bottomed position
1103 restraints of different geometry on one particle. For example, applying
1104 a cylinder and a layer in $z$ keeps a particle within a
1105 disk. Applying three layers in $x$, $y$, and $z$ keeps the particle within a cuboid.
1107 In addition, it is possible to invert the restrained region with the
1108 unrestrained region, leading to a potential that acts to keep the particle {\it outside} of the volume
1109 defined by $\ve{R}_i$, $g$, and $r_\mathrm{fb}$. That feature is
1110 switched on by defining a negative $r_\mathrm{fb}$ in the
1111 topology. The following potential is used (Figure~\ref{fig:fbposres}B):
1113 V_\mathrm{fb}^{\mathrm{inv}}(\ve{r}_i) = \frac{1}{2}k_\mathrm{fb}
1114 [d_g(\ve{r}_i;\ve{R}_i) - |r_\mathrm{fb}|]^2\,
1115 H[ -(d_g(\ve{r}_i;\ve{R}_i) - |r_\mathrm{fb}|)].
1120 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1121 \subsection{Angle restraints\swapindexquiet{angle}{restraint}}
1122 \label{subsec:anglerestraint}
1123 These are used to restrain the angle between two pairs of particles
1124 or between one pair of particles and the $z$-axis.
1125 The functional form is similar to that of a proper dihedral.
1126 For two pairs of atoms:
1128 V_{ar}(\ve{r}_i,\ve{r}_j,\ve{r}_k,\ve{r}_l)
1129 = k_{ar}(1 - \cos(n (\theta - \theta_0))
1132 \theta = \arccos\left(\frac{\ve{r}_j -\ve{r}_i}{\|\ve{r}_j -\ve{r}_i\|}
1133 \cdot \frac{\ve{r}_l -\ve{r}_k}{\|\ve{r}_l -\ve{r}_k\|} \right)
1135 For one pair of atoms and the $z$-axis:
1137 V_{ar}(\ve{r}_i,\ve{r}_j) = k_{ar}(1 - \cos(n (\theta - \theta_0))
1140 \theta = \arccos\left(\frac{\ve{r}_j -\ve{r}_i}{\|\ve{r}_j -\ve{r}_i\|}
1141 \cdot \left( \begin{array}{c} 0 \\ 0 \\ 1 \\ \end{array} \right) \right)
1143 A multiplicity ($n$) of 2 is useful when you do not want to distinguish
1144 between parallel and anti-parallel vectors.
1145 The equilibrium angle $\theta$ should be between 0 and 180 degrees
1146 for multiplicity 1 and between 0 and 90 degrees for multiplicity 2.
1149 \subsection{Dihedral restraints\swapindexquiet{dihedral}{restraint}}
1150 \label{subsec:dihedralrestraint}
1151 These are used to restrain the dihedral angle $\phi$ defined by four particles
1152 as in an improper dihedral (sec.~\ref{sec:imp}) but with a slightly
1153 modified potential. Using:
1155 \phi' = \left(\phi-\phi_0\right) ~{\rm MOD}~ 2\pi
1157 where $\phi_0$ is the reference angle, the potential is defined as:
1159 V_{dihr}(\phi') ~=~ \left\{
1160 \begin{array}{lcllll}
1161 \half k_{dihr}(\phi'-\phi_0-\Delta\phi)^2
1162 &\mbox{for}& \phi' & > & \Delta\phi \\[1.5ex]
1163 0 &\mbox{for}& \phi' & \le & \Delta\phi \\[1.5ex]
1167 where $\Delta\phi$ is a user defined angle and $k_{dihr}$ is the force
1169 {\bf Note} that in the input in topology files, angles are given in degrees and
1170 force constants in kJ/mol/rad$^2$.
1171 %} % Brace matches ifthenelse test for gmxlite
1173 \subsection{Distance restraints\swapindexquiet{distance}{restraint}}
1174 \label{subsec:distancerestraint}
1176 add a penalty to the potential when the distance between specified
1177 pairs of atoms exceeds a threshold value. They are normally used to
1178 impose experimental restraints from, for instance, experiments in nuclear
1179 magnetic resonance (NMR), on the motion of the system. Thus, MD can be
1180 used for structure refinement using NMR data\index{nmr
1181 refinement}\index{refinement,nmr}.
1182 In {\gromacs} there are three ways to impose restraints on pairs of atoms:
1184 \item Simple harmonic restraints: use {\tt [ bonds ]} type 6
1185 %\ifthenelse{\equal{\gmxlite}{1}}
1187 {(see \secref{excl}).}
1188 \item\label{subsec:harmonicrestraint}Piecewise linear/harmonic restraints: {\tt [ bonds ]} type 10.
1189 \item Complex NMR distance restraints, optionally with pair, time and/or
1192 The last two options will be detailed now.
1194 The potential form for distance restraints is quadratic below a specified
1195 lower bound and between two specified upper bounds, and linear beyond the
1196 largest bound (see \figref{dist}).
1198 V_{dr}(r_{ij}) ~=~ \left\{
1199 \begin{array}{lcllllll}
1200 \half k_{dr}(r_{ij}-r_0)^2
1201 &\mbox{for}& & & r_{ij} & < & r_0 \\[1.5ex]
1202 0 &\mbox{for}& r_0 & \le & r_{ij} & < & r_1 \\[1.5ex]
1203 \half k_{dr}(r_{ij}-r_1)^2
1204 &\mbox{for}& r_1 & \le & r_{ij} & < & r_2 \\[1.5ex]
1205 \half k_{dr}(r_2-r_1)(2r_{ij}-r_2-r_1)
1206 &\mbox{for}& r_2 & \le & r_{ij} & &
1212 \centerline{\includegraphics[width=8cm]{plots/f-dr}}
1213 \caption{Distance Restraint potential.}
1220 \begin{array}{lcllllll}
1221 -k_{dr}(r_{ij}-r_0)\frac{\rvij}{r_{ij}}
1222 &\mbox{for}& & & r_{ij} & < & r_0 \\[1.5ex]
1223 0 &\mbox{for}& r_0 & \le & r_{ij} & < & r_1 \\[1.5ex]
1224 -k_{dr}(r_{ij}-r_1)\frac{\rvij}{r_{ij}}
1225 &\mbox{for}& r_1 & \le & r_{ij} & < & r_2 \\[1.5ex]
1226 -k_{dr}(r_2-r_1)\frac{\rvij}{r_{ij}}
1227 &\mbox{for}& r_2 & \le & r_{ij} & &
1231 For restraints not derived from NMR data, this functionality
1232 will usually suffice and a section of {\tt [ bonds ]} type 10
1233 can be used to apply individual restraints between pairs of
1234 %\ifthenelse{\equal{\gmxlite}{1}}{atoms.}{
1235 atoms, see \ssecref{topfile}.
1236 %} % Brace matches ifthenelse test for gmxlite
1237 For applying restraints derived from NMR measurements, more complex
1238 functionality might be required, which is provided through
1239 the {\tt [~distance_restraints~]} section and is described below.
1241 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1242 \subsubsection{Time averaging\swapindexquiet{time-averaged}{distance restraint}}
1243 Distance restraints based on instantaneous distances can potentially reduce
1244 the fluctuations in a molecule significantly. This problem can be overcome by restraining
1245 to a {\em time averaged} distance~\cite{Torda89}.
1246 The forces with time averaging are:
1249 \begin{array}{lcllllll}
1250 -k^a_{dr}(\bar{r}_{ij}-r_0)\frac{\rvij}{r_{ij}}
1251 &\mbox{for}& & & \bar{r}_{ij} & < & r_0 \\[1.5ex]
1252 0 &\mbox{for}& r_0 & \le & \bar{r}_{ij} & < & r_1 \\[1.5ex]
1253 -k^a_{dr}(\bar{r}_{ij}-r_1)\frac{\rvij}{r_{ij}}
1254 &\mbox{for}& r_1 & \le & \bar{r}_{ij} & < & r_2 \\[1.5ex]
1255 -k^a_{dr}(r_2-r_1)\frac{\rvij}{r_{ij}}
1256 &\mbox{for}& r_2 & \le & \bar{r}_{ij} & &
1259 where $\bar{r}_{ij}$ is given by an exponential running average with decay time $\tau$:
1261 \bar{r}_{ij} ~=~ < r_{ij}^{-3} >^{-1/3}
1264 The force constant $k^a_{dr}$ is switched on slowly to compensate for
1265 the lack of history at the beginning of the simulation:
1267 k^a_{dr} = k_{dr} \left(1-\exp\left(-\frac{t}{\tau}\right)\right)
1269 Because of the time averaging, we can no longer speak of a distance restraint
1272 This way an atom can satisfy two incompatible distance restraints
1273 {\em on average} by moving between two positions.
1274 An example would be an amino acid side-chain that is rotating around
1275 its $\chi$ dihedral angle, thereby coming close to various other groups.
1276 Such a mobile side chain can give rise to multiple NOEs that can not be
1277 fulfilled by a single structure.
1279 The computation of the time
1280 averaged distance in the {\tt mdrun} program is done in the following fashion:
1283 \overline{r^{-3}}_{ij}(0) &=& r_{ij}(0)^{-3} \\
1284 \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]
1285 \label{eqn:ravdisre}
1289 When a pair is within the bounds, it can still feel a force
1290 because the time averaged distance can still be beyond a bound.
1291 To prevent the protons from being pulled too close together, a mixed
1292 approach can be used. In this approach, the penalty is zero when the
1293 instantaneous distance is within the bounds, otherwise the violation is
1294 the square root of the product of the instantaneous violation and the
1295 time averaged violation:
1298 \begin{array}{lclll}
1299 k^a_{dr}\sqrt{(r_{ij}-r_0)(\bar{r}_{ij}-r_0)}\frac{\rvij}{r_{ij}}
1300 & \mbox{for} & r_{ij} < r_0 & \mbox{and} & \bar{r}_{ij} < r_0 \\[1.5ex]
1302 \mbox{min}\left(\sqrt{(r_{ij}-r_1)(\bar{r}_{ij}-r_1)},r_2-r_1\right)
1303 \frac{\rvij}{r_{ij}}
1304 & \mbox{for} & r_{ij} > r_1 & \mbox{and} & \bar{r}_{ij} > r_1 \\[1.5ex]
1309 \subsubsection{Averaging over multiple pairs\swapindexquiet{ensemble-averaged}{distance restraint}}
1311 Sometimes it is unclear from experimental data which atom pair
1312 gives rise to a single NOE, in other occasions it can be obvious that
1313 more than one pair contributes due to the symmetry of the system, {\eg} a
1314 methyl group with three protons. For such a group, it is not possible
1315 to distinguish between the protons, therefore they should all be taken into
1316 account when calculating the distance between this methyl group and another
1317 proton (or group of protons).
1318 Due to the physical nature of magnetic resonance, the intensity of the
1319 NOE signal is inversely proportional to the sixth power of the inter-atomic
1321 Thus, when combining atom pairs,
1322 a fixed list of $N$ restraints may be taken together,
1323 where the apparent ``distance'' is given by:
1325 r_N(t) = \left [\sum_{n=1}^{N} \bar{r}_{n}(t)^{-6} \right]^{-1/6}
1328 where we use $r_{ij}$ or \eqnref{rav} for the $\bar{r}_{n}$.
1329 The $r_N$ of the instantaneous and time-averaged distances
1330 can be combined to do a mixed restraining, as indicated above.
1331 As more pairs of protons contribute to the same NOE signal, the intensity
1332 will increase, and the summed ``distance'' will be shorter than any of
1333 its components due to the reciprocal summation.
1335 There are two options for distributing the forces over the atom pairs.
1336 In the conservative option, the force is defined as the derivative of the
1337 restraint potential with respect to the coordinates. This results in
1338 a conservative potential when time averaging is not used.
1339 The force distribution over the pairs is proportional to $r^{-6}$.
1340 This means that a close pair feels a much larger force than a distant pair,
1341 which might lead to a molecule that is ``too rigid.''
1342 The other option is an equal force distribution. In this case each pair
1343 feels $1/N$ of the derivative of the restraint potential with respect to
1344 $r_N$. The advantage of this method is that more conformations might be
1345 sampled, but the non-conservative nature of the forces can lead to
1346 local heating of the protons.
1348 It is also possible to use {\em ensemble averaging} using multiple
1349 (protein) molecules. In this case the bounds should be lowered as in:
1352 r_1 &~=~& r_1 * M^{-1/6} \\
1353 r_2 &~=~& r_2 * M^{-1/6}
1356 where $M$ is the number of molecules. The {\gromacs} preprocessor {\tt grompp}
1357 can do this automatically when the appropriate option is given.
1358 The resulting ``distance'' is
1359 then used to calculate the scalar force according to:
1363 ~& 0 \hspace{4cm} & r_{N} < r_1 \\
1364 & k_{dr}(r_{N}-r_1)\frac{\rvij}{r_{ij}} & r_1 \le r_{N} < r_2 \\
1365 & k_{dr}(r_2-r_1)\frac{\rvij}{r_{ij}} & r_{N} \ge r_2
1368 where $i$ and $j$ denote the atoms of all the
1369 pairs that contribute to the NOE signal.
1371 \subsubsection{Using distance restraints}
1373 A list of distance restrains based on NOE data can be added to a molecule
1374 definition in your topology file, like in the following example:
1377 [ distance_restraints ]
1378 ; ai aj type index type' low up1 up2 fac
1379 10 16 1 0 1 0.0 0.3 0.4 1.0
1380 10 28 1 1 1 0.0 0.3 0.4 1.0
1381 10 46 1 1 1 0.0 0.3 0.4 1.0
1382 16 22 1 2 1 0.0 0.3 0.4 2.5
1383 16 34 1 3 1 0.0 0.5 0.6 1.0
1386 In this example a number of features can be found. In columns {\tt
1387 ai} and {\tt aj} you find the atom numbers of the particles to be
1388 restrained. The {\tt type} column should always be 1. As explained in
1389 ~\ssecref{distancerestraint}, multiple distances can contribute to a single NOE
1390 signal. In the topology this can be set using the {\tt index}
1391 column. In our example, the restraints 10-28 and 10-46 both have index
1392 1, therefore they are treated simultaneously. An extra requirement
1393 for treating restraints together is that the restraints must be on
1394 successive lines, without any other intervening restraint. The {\tt
1395 type'} column will usually be 1, but can be set to 2 to obtain a
1396 distance restraint that will never be time- and ensemble-averaged;
1397 this can be useful for restraining hydrogen bonds. The columns {\tt
1398 low}, {\tt up1}, and {\tt up2} hold the values of $r_0$, $r_1$, and
1399 $r_2$ from ~\eqnref{disre}. In some cases it can be useful to have
1400 different force constants for some restraints; this is controlled by
1401 the column {\tt fac}. The force constant in the parameter file is
1402 multiplied by the value in the column {\tt fac} for each restraint.
1403 %} % Brace matches ifthenelse test for gmxlite
1405 \newcommand{\SSS}{{\mathbf S}}
1406 \newcommand{\DD}{{\mathbf D}}
1407 \newcommand{\RR}{{\mathbf R}}
1409 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1410 \subsection{Orientation restraints\swapindexquiet{orientation}{restraint}}
1411 \label{subsec:orientationrestraint}
1412 This section describes how orientations between vectors,
1413 as measured in certain NMR experiments, can be calculated
1414 and restrained in MD simulations.
1415 The presented refinement methodology and a comparison of results
1416 with and without time and ensemble averaging have been
1417 published~\cite{Hess2003}.
1418 \subsubsection{Theory}
1419 In an NMR experiment, orientations of vectors can be measured when a
1420 molecule does not tumble completely isotropically in the solvent.
1421 Two examples of such orientation measurements are
1422 residual \normindex{dipolar couplings}
1423 (between two nuclei) or chemical shift anisotropies.
1424 An observable for a vector $\ve{r}_i$ can be written as follows:
1426 \delta_i = \frac{2}{3} \mbox{tr}(\SSS\DD_i)
1428 where $\SSS$ is the dimensionless order tensor of the molecule.
1429 The tensor $\DD_i$ is given by:
1432 \DD_i = \frac{c_i}{\|\ve{r}_i\|^\alpha} \left(
1434 %3 r_x r_x - \ve{r}\cdot\ve{r} & 3 r_x r_y & 3 r_x r_z \\
1435 %3 r_x r_y & 3 r_y r_y - \ve{r}\cdot\ve{r} & 3yz \\
1436 %3 r_x r_z & 3 r_y r_z & 3 r_z r_z - \ve{r}\cdot\ve{r}
1437 %\end{array} \right)
1439 3 x x - 1 & 3 x y & 3 x z \\
1440 3 x y & 3 y y - 1 & 3 y z \\
1441 3 x z & 3 y z & 3 z z - 1 \\
1446 x=\frac{r_{i,x}}{\|\ve{r}_i\|}, \quad
1447 y=\frac{r_{i,y}}{\|\ve{r}_i\|}, \quad
1448 z=\frac{r_{i,z}}{\|\ve{r}_i\|}
1450 For a dipolar coupling $\ve{r}_i$ is the vector connecting the two
1451 nuclei, $\alpha=3$ and the constant $c_i$ is given by:
1453 c_i = \frac{\mu_0}{4\pi} \gamma_1^i \gamma_2^i \frac{\hbar}{4\pi}
1455 where $\gamma_1^i$ and $\gamma_2^i$ are the gyromagnetic ratios of the
1458 The order tensor is symmetric and has trace zero. Using a rotation matrix
1459 ${\mathbf T}$ it can be transformed into the following form:
1461 {\mathbf T}^T \SSS {\mathbf T} = s \left( \begin{array}{ccc}
1462 -\frac{1}{2}(1-\eta) & 0 & 0 \\
1463 0 & -\frac{1}{2}(1+\eta) & 0 \\
1467 where $-1 \leq s \leq 1$ and $0 \leq \eta \leq 1$.
1468 $s$ is called the order parameter and $\eta$ the asymmetry of the
1469 order tensor $\SSS$. When the molecule tumbles isotropically in the
1470 solvent, $s$ is zero, and no orientational effects can be observed
1471 because all $\delta_i$ are zero.
1475 \subsubsection{Calculating orientations in a simulation}
1476 For reasons which are explained below, the $\DD$ matrices are calculated
1477 which respect to a reference orientation of the molecule. The orientation
1478 is defined by a rotation matrix $\RR$, which is needed to least-squares fit
1479 the current coordinates of a selected set of atoms onto
1480 a reference conformation. The reference conformation is the starting
1481 conformation of the simulation. In case of ensemble averaging, which will
1482 be treated later, the structure is taken from the first subsystem.
1483 The calculated $\DD_i^c$ matrix is given by:
1486 \DD_i^c(t) = \RR(t) \DD_i(t) \RR^T(t)
1488 The calculated orientation for vector $i$ is given by:
1490 \delta^c_i(t) = \frac{2}{3} \mbox{tr}(\SSS(t)\DD_i^c(t))
1492 The order tensor $\SSS(t)$ is usually unknown.
1493 A reasonable choice for the order tensor is the tensor
1494 which minimizes the (weighted) mean square difference between the calculated
1495 and the observed orientations:
1498 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
1500 To properly combine different types of measurements, the unit of $w_i$ should
1501 be such that all terms are dimensionless. This means the unit of $w_i$
1502 is the unit of $\delta_i$ to the power $-2$.
1503 {\bf Note} that scaling all $w_i$ with a constant factor does not influence
1506 \subsubsection{Time averaging}
1507 Since the tensors $\DD_i$ fluctuate rapidly in time, much faster than can
1508 be observed in an experiment, they should be averaged over time in the simulation.
1509 However, in a simulation the time and the number of copies of
1510 a molecule are limited. Usually one can not obtain a converged average
1511 of the $\DD_i$ tensors over all orientations of the molecule.
1512 If one assumes that the average orientations of the $\ve{r}_i$ vectors
1513 within the molecule converge much faster than the tumbling time of
1514 the molecule, the tensor can be averaged in an axis system that
1515 rotates with the molecule, as expressed by equation~(\ref{D_rot}).
1516 The time-averaged tensors are calculated
1517 using an exponentially decaying memory function:
1519 \DD^a_i(t) = \frac{\displaystyle
1520 \int_{u=t_0}^t \DD^c_i(u) \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
1522 \int_{u=t_0}^t \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
1525 Assuming that the order tensor $\SSS$ fluctuates slower than the
1526 $\DD_i$, the time-averaged orientation can be calculated as:
1528 \delta_i^a(t) = \frac{2}{3} \mbox{tr}(\SSS(t) \DD_i^a(t))
1530 where the order tensor $\SSS(t)$ is calculated using expression~(\ref{S_msd})
1531 with $\delta_i^c(t)$ replaced by $\delta_i^a(t)$.
1533 \subsubsection{Restraining}
1534 The simulated structure can be restrained by applying a force proportional
1535 to the difference between the calculated and the experimental orientations.
1536 When no time averaging is applied, a proper potential can be defined as:
1538 V = \frac{1}{2} k \sum_{i=1}^N w_i (\delta_i^c (t) -\delta_i^{exp})^2
1540 where the unit of $k$ is the unit of energy.
1541 Thus the effective force constant for restraint $i$ is $k w_i$.
1542 The forces are given by minus the gradient of $V$.
1543 The force $\ve{F}\!_i$ working on vector $\ve{r}_i$ is:
1546 & = & - \frac{\mbox{d} V}{\mbox{d}\ve{r}_i} \\
1547 & = & -k w_i (\delta_i^c (t) -\delta_i^{exp}) \frac{\mbox{d} \delta_i (t)}{\mbox{d}\ve{r}_i} \\
1548 & = & -k w_i (\delta_i^c (t) -\delta_i^{exp})
1549 \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)
1552 \subsubsection{Ensemble averaging}
1553 Ensemble averaging can be applied by simulating a system of $M$ subsystems
1554 that each contain an identical set of orientation restraints. The systems only
1555 interact via the orientation restraint potential which is defined as:
1557 V = M \frac{1}{2} k \sum_{i=1}^N w_i
1558 \langle \delta_i^c (t) -\delta_i^{exp} \rangle^2
1560 The force on vector $\ve{r}_{i,m}$ in subsystem $m$ is given by:
1562 \ve{F}\!_{i,m}(t) = - \frac{\mbox{d} V}{\mbox{d}\ve{r}_{i,m}} =
1563 -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}} \\
1566 \subsubsection{Time averaging}
1567 When using time averaging it is not possible to define a potential.
1568 We can still define a quantity that gives a rough idea of the energy
1569 stored in the restraints:
1571 V = M \frac{1}{2} k^a \sum_{i=1}^N w_i
1572 \langle \delta_i^a (t) -\delta_i^{exp} \rangle^2
1574 The force constant $k_a$ is switched on slowly to compensate for the
1575 lack of history at times close to $t_0$. It is exactly proportional
1576 to the amount of average that has been accumulated:
1579 k \, \frac{1}{\tau}\int_{u=t_0}^t \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
1581 What really matters is the definition of the force. It is chosen to
1582 be proportional to the square root of the product of the time-averaged
1583 and the instantaneous deviation.
1584 Using only the time-averaged deviation induces large oscillations.
1585 The force is given by:
1588 %\left\{ \begin{array}{ll}
1589 %0 & \mbox{for} \quad \langle \delta_i^a (t) -\delta_i^{exp} \rangle \langle \delta_i (t) -\delta_i^{exp} \rangle \leq 0 \\
1590 %... & \mbox{for} \quad \langle \delta_i^a (t) -\delta_i^{exp} \rangle \langle \delta_i (t) -\delta_i^{exp} \rangle > 0
1593 \left\{ \begin{array}{ll}
1594 0 & \quad \mbox{for} \quad a\, b \leq 0 \\
1596 k^a w_i \frac{a}{|a|} \sqrt{a\, b} \, \frac{\mbox{d} \delta_{i,m}^c (t)}{\mbox{d}\ve{r}_{i,m}}
1597 & \quad \mbox{for} \quad a\, b > 0
1602 a &=& \langle \delta_i^a (t) -\delta_i^{exp} \rangle \\
1603 b &=& \langle \delta_i^c (t) -\delta_i^{exp} \rangle
1606 \subsubsection{Using orientation restraints}
1607 Orientation restraints can be added to a molecule definition in
1608 the topology file in the section {\tt [~orientation_restraints~]}.
1609 Here we give an example section containing five N-H residual dipolar
1610 coupling restraints:
1613 [ orientation_restraints ]
1614 ; ai aj type exp. label alpha const. obs. weight
1616 31 32 1 1 3 3 6.083 -6.73 1.0
1617 43 44 1 1 4 3 6.083 -7.87 1.0
1618 55 56 1 1 5 3 6.083 -7.13 1.0
1619 65 66 1 1 6 3 6.083 -2.57 1.0
1620 73 74 1 1 7 3 6.083 -2.10 1.0
1623 The unit of the observable is Hz, but one can choose any other unit.
1625 ai} and {\tt aj} you find the atom numbers of the particles to be
1626 restrained. The {\tt type} column should always be 1.
1627 The {\tt exp.} column denotes the experiment number, starting
1628 at 1. For each experiment a separate order tensor $\SSS$
1629 is optimized. The label should be a unique number larger than zero
1630 for each restraint. The {\tt alpha} column contains the power $\alpha$
1631 that is used in equation~(\ref{orient_def}) to calculate the orientation.
1632 The {\tt const.} column contains the constant $c_i$ used in the same
1633 equation. The constant should have the unit of the observable times
1634 nm$^\alpha$. The column {\tt obs.} contains the observable, in any
1635 unit you like. The last column contains the weights $w_i$; the unit
1636 should be the inverse of the square of the unit of the observable.
1638 Some parameters for orientation restraints can be specified in the
1639 {\tt grompp.mdp} file, for a study of the effect of different
1640 force constants and averaging times and ensemble averaging see~\cite{Hess2003}.
1641 %} % Brace matches ifthenelse test for gmxlite
1643 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1644 \section{Polarization}
1645 Polarization can be treated by {\gromacs} by attaching
1646 \normindex{shell} (\normindex{Drude}) particles to atoms and/or
1647 virtual sites. The energy of the shell particle is then minimized at
1648 each time step in order to remain on the Born-Oppenheimer surface.
1650 \subsection{Simple polarization}
1651 This is merely a harmonic potential with equilibrium distance 0.
1653 \subsection{Water polarization}
1654 A special potential for water that allows anisotropic polarization of
1655 a single shell particle~\cite{Maaren2001a}.
1657 \subsection{Thole polarization}
1658 Based on early work by \normindex{Thole}~\cite{Thole81}, Roux and
1659 coworkers have implemented potentials for molecules like
1660 ethanol~\cite{Lamoureux2003a,Lamoureux2003b,Noskov2005a}. Within such
1661 molecules, there are intra-molecular interactions between shell
1662 particles, however these must be screened because full Coulomb would
1663 be too strong. The potential between two shell particles $i$ and $j$ is:
1664 \newcommand{\rbij}{\bar{r}_{ij}}
1666 V_{thole} ~=~ \frac{q_i q_j}{r_{ij}}\left[1-\left(1+\frac{\rbij}{2}\right){\rm exp}^{-\rbij}\right]
1668 {\bf Note} that there is a sign error in Equation~1 of Noskov {\em et al.}~\cite{Noskov2005a}:
1670 \rbij ~=~ a\frac{r_{ij}}{(\alpha_i \alpha_j)^{1/6}}
1672 where $a$ is a magic (dimensionless) constant, usually chosen to be
1673 2.6~\cite{Noskov2005a}; $\alpha_i$ and $\alpha_j$ are the polarizabilities
1674 of the respective shell particles.
1676 %} % Brace matches ifthenelse test for gmxlite
1678 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1679 \section{Free energy interactions}
1681 \index{free energy interactions}
1682 \newcommand{\LAM}{\lambda}
1683 \newcommand{\LL}{(1-\LAM)}
1684 \newcommand{\dvdl}[1]{\frac{\partial #1}{\partial \LAM}}
1685 This section describes the $\lambda$-dependence of the potentials
1686 used for free energy calculations (see \secref{fecalc}).
1687 All common types of potentials and constraints can be
1688 interpolated smoothly from state A ($\lambda=0$) to state B
1689 ($\lambda=1$) and vice versa.
1690 All bonded interactions are interpolated by linear interpolation
1691 of the interaction parameters. Non-bonded interactions can be
1692 interpolated linearly or via soft-core interactions.
1694 Starting in {\gromacs} 4.6, $\lambda$ is a vector, allowing different
1695 components of the free energy transformation to be carried out at
1696 different rates. Coulomb, Lennard-Jones, bonded, and restraint terms
1697 can all be controlled independently, as described in the {\tt .mdp}
1700 \subsubsection{Harmonic potentials}
1701 The example given here is for the bond potential, which is harmonic
1702 in {\gromacs}. However, these equations apply to the angle potential
1703 and the improper dihedral potential as well.
1705 V_b &=&\half\left[\LL k_b^A +
1706 \LAM k_b^B\right] \left[b - \LL b_0^A - \LAM b_0^B\right]^2 \\
1707 \dvdl{V_b}&=&\half(k_b^B-k_b^A)
1708 \left[b - \LL b_0^A + \LAM b_0^B\right]^2 +
1710 & & \phantom{\half}(b_0^A-b_0^B) \left[b - \LL b_0^A -\LAM b_0^B\right]
1711 \left[\LL k_b^A + \LAM k_b^B \right]
1714 \subsubsection{\gromosv{96} bonds and angles}
1715 Fourth-power bond stretching and cosine-based angle potentials
1716 are interpolated by linear interpolation of the force constant
1717 and the equilibrium position. Formulas are not given here.
1719 \subsubsection{Proper dihedrals}
1720 For the proper dihedrals, the equations are somewhat more complicated:
1722 V_d &=&\left[\LL k_d^A + \LAM k_d^B \right]
1723 \left( 1+ \cos\left[n_{\phi} \phi -
1724 \LL \phi_s^A - \LAM \phi_s^B
1726 \dvdl{V_d}&=&(k_d^B-k_d^A)
1729 n_{\phi} \phi- \LL \phi_s^A - \LAM \phi_s^B
1733 &&(\phi_s^B - \phi_s^A) \left[\LL k_d^A - \LAM k_d^B\right]
1734 \sin\left[ n_{\phi}\phi - \LL \phi_s^A - \LAM \phi_s^B \right]
1736 {\bf Note:} that the multiplicity $n_{\phi}$ can not be parameterized
1737 because the function should remain periodic on the interval $[0,2\pi]$.
1739 \subsubsection{Tabulated bonded interactions}
1740 For tabulated bonded interactions only the force constant can interpolated:
1742 V &=& (\LL k^A + \LAM k^B) \, f \\
1743 \dvdl{V} &=& (k^B - k^A) \, f
1746 \subsubsection{Coulomb interaction}
1747 The \normindex{Coulomb} interaction between two particles
1748 of which the charge varies with $\LAM$ is:
1750 V_c &=& \frac{f}{\epsrf \rij}\left[\LL q_i^A q_j^A + \LAM\, q_i^B q_j^B\right] \\
1751 \dvdl{V_c}&=& \frac{f}{\epsrf \rij}\left[- q_i^A q_j^A + q_i^B q_j^B\right]
1753 where $f = \frac{1}{4\pi \varepsilon_0} = 138.935\,485$ (see \chref{defunits}).
1755 \subsubsection{Coulomb interaction with \normindex{reaction field}}
1756 The Coulomb interaction including a reaction field, between two particles
1757 of which the charge varies with $\LAM$ is:
1759 V_c &=& f\left[\frac{1}{\rij} + k_{rf}~ \rij^2 -c_{rf}\right]
1760 \left[\LL q_i^A q_j^A + \LAM\, q_i^B q_j^B\right] \\
1761 \dvdl{V_c}&=& f\left[\frac{1}{\rij} + k_{rf}~ \rij^2 -c_{rf}\right]
1762 \left[- q_i^A q_j^A + q_i^B q_j^B\right]
1763 \label{eq:dVcoulombdlambda}
1765 {\bf Note} that the constants $k_{rf}$ and $c_{rf}$ are
1766 defined using the dielectric
1767 constant $\epsrf$ of the medium (see \secref{coulrf}).
1769 \subsubsection{Lennard-Jones interaction}
1770 For the \normindex{Lennard-Jones} interaction between two particles
1771 of which the {\em atom type} varies with $\LAM$ we can write:
1773 V_{LJ} &=& \frac{\LL C_{12}^A + \LAM\, C_{12}^B}{\rij^{12}} -
1774 \frac{\LL C_6^A + \LAM\, C_6^B}{\rij^6} \\
1775 \dvdl{V_{LJ}}&=&\frac{C_{12}^B - C_{12}^A}{\rij^{12}} -
1776 \frac{C_6^B - C_6^A}{\rij^6}
1777 \label{eq:dVljdlambda}
1779 It should be noted that it is also possible to express a pathway from
1780 state A to state B using $\sigma$ and $\epsilon$ (see \eqnref{sigeps}).
1781 It may seem to make sense physically to vary the force field parameters
1782 $\sigma$ and $\epsilon$ rather
1783 than the derived parameters $C_{12}$ and $C_{6}$.
1784 However, the difference between the pathways in parameter space
1785 is not large, and the free energy itself
1786 does not depend on the pathway, so we use the simple formulation
1789 \subsubsection{Kinetic Energy}
1790 When the mass of a particle changes, there is also a contribution of
1791 the kinetic energy to the free energy (note that we can not write
1792 the momentum \ve{p} as m\ve{v}, since that would result
1793 in the sign of $\dvdl{E_k}$ being incorrect~\cite{Gunsteren98a}):
1796 E_k &=& \half\frac{\ve{p}^2}{\LL m^A + \LAM m^B} \\
1797 \dvdl{E_k}&=& -\half\frac{\ve{p}^2(m^B-m^A)}{(\LL m^A + \LAM m^B)^2}
1799 after taking the derivative, we {\em can} insert \ve{p} = m\ve{v}, such that:
1801 \dvdl{E_k}~=~ -\half\ve{v}^2(m^B-m^A)
1804 \subsubsection{Constraints}
1805 \label{subsubsec:constraints}
1806 \newcommand{\clam}{C_{\lambda}}
1807 The constraints are formally part of the Hamiltonian, and therefore
1808 they give a contribution to the free energy. In {\gromacs} this can be
1809 calculated using the \normindex{LINCS} or the \normindex{SHAKE} algorithm.
1810 If we have a number of constraint equations $g_k$:
1812 g_k = \ve{r}_{k} - d_{k}
1814 where $\ve{r}_k$ is the distance vector between two particles and
1815 $d_k$ is the constraint distance between the two particles, we can write
1816 this using a $\LAM$-dependent distance as
1818 g_k = \ve{r}_{k} - \left(\LL d_{k}^A + \LAM d_k^B\right)
1820 the contribution $\clam$
1821 to the Hamiltonian using Lagrange multipliers $\lambda$:
1823 \clam &=& \sum_k \lambda_k g_k \\
1824 \dvdl{\clam} &=& \sum_k \lambda_k \left(d_k^B-d_k^A\right)
1828 \subsection{Soft-core interactions\index{soft-core interactions}}
1830 \centerline{\includegraphics[height=6cm]{plots/softcore}}
1831 \caption{Soft-core interactions at $\LAM=0.5$, with $p=2$ and
1832 $C_6^A=C_{12}^A=C_6^B=C_{12}^B=1$.}
1833 \label{fig:softcore}
1835 In a free-energy calculation where particles grow out of nothing, or
1836 particles disappear, using the the simple linear interpolation of the
1837 Lennard-Jones and Coulomb potentials as described in Equations~\ref{eq:dVljdlambda}
1838 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
1839 close to each other, leading to large fluctuations in the measured values of
1840 $\partial V/\partial \LAM$ (which, because of the simple linear
1841 interpolation, depends on the potentials at both the endpoints of $\LAM$).
1843 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
1844 involved at $\LAM$ values between 0 and 1, but not \emph{at} $\LAM=0$
1847 In {\gromacs} the soft-core potentials $V_{sc}$ are shifted versions of the
1848 regular potentials, so that the singularity in the potential and its
1849 derivatives at $r=0$ is never reached:
1851 V_{sc}(r) &=& \LL V^A(r_A) + \LAM V^B(r_B)
1853 r_A &=& \left(\alpha \sigma_A^6 \LAM^p + r^6 \right)^\frac{1}{6}
1855 r_B &=& \left(\alpha \sigma_B^6 \LL^p + r^6 \right)^\frac{1}{6}
1857 where $V^A$ and $V^B$ are the normal ``hard core'' Van der Waals or
1858 electrostatic potentials in state A ($\LAM=0$) and state B ($\LAM=1$)
1859 respectively, $\alpha$ is the soft-core parameter (set with {\tt sc_alpha}
1860 in the {\tt .mdp} file), $p$ is the soft-core $\LAM$ power (set with
1861 {\tt sc_power}), $\sigma$ is the radius of the interaction, which is
1862 $(C_{12}/C_6)^{1/6}$ or an input parameter ({\tt sc_sigma}) when $C_6$
1863 or $C_{12}$ is zero.
1865 For intermediate $\LAM$, $r_A$ and $r_B$ alter the interactions very little
1866 for $r > \alpha^{1/6} \sigma$ and quickly switch the soft-core
1867 interaction to an almost constant value for smaller $r$ (\figref{softcore}).
1870 F_{sc}(r) = -\frac{\partial V_{sc}(r)}{\partial r} =
1871 \LL F^A(r_A) \left(\frac{r}{r_A}\right)^5 +
1872 \LAM F^B(r_B) \left(\frac{r}{r_B}\right)^5
1874 where $F^A$ and $F^B$ are the ``hard core'' forces.
1875 The contribution to the derivative of the free energy is:
1877 \dvdl{V_{sc}(r)} & = &
1878 V^B(r_B) -V^A(r_A) +
1879 \LL \frac{\partial V^A(r_A)}{\partial r_A}
1880 \frac{\partial r_A}{\partial \LAM} +
1881 \LAM\frac{\partial V^B(r_B)}{\partial r_B}
1882 \frac{\partial r_B}{\partial \LAM}
1885 V^B(r_B) -V^A(r_A) + \nonumber \\
1888 \left[ \LAM F^B(r_B) r^{-5}_B \sigma_B^6 \LL^{p-1} -
1889 \LL F^A(r_A) r^{-5}_A \sigma_A^6 \LAM^{p-1} \right]
1892 The original GROMOS Lennard-Jones soft-core function~\cite{Beutler94}
1893 uses $p=2$, but $p=1$ gives a smoother $\partial H/\partial\LAM$ curve.
1894 %When the changes between the two states involve both the disappearing
1895 %and appearing of atoms, it is important that the overlapping of atoms
1896 %happens around $\LAM=0.5$. This can usually be achieved with
1897 %$\alpha$$\approx0.7$ for $p=1$ and $\alpha$$\approx1.5$ for $p=2$.
1898 %MRS: this is now eliminated as of 4.6, since changes between atoms are done linearly.
1900 Another issue that should be considered is the soft-core effect of hydrogens
1901 without Lennard-Jones interaction. Their soft-core $\sigma$ is
1902 set with {\tt sc-sigma} in the {\tt .mdp} file. These hydrogens
1903 produce peaks in $\partial H/\partial\LAM$ at $\LAM$ is 0 and/or 1 for $p=1$
1904 and close to 0 and/or 1 with $p=2$. Lowering {\tt\mbox{sc-sigma}} will decrease
1905 this effect, but it will also increase the interactions with hydrogens
1906 relative to the other interactions in the soft-core state.
1908 When soft core potentials are selected (by setting {\tt sc-alpha} \textgreater
1909 0), and the Coulomb and Lennard-Jones potentials are turned on or off
1910 sequentially, then the Coulombic interaction is turned off linearly,
1911 rather than using soft core interactions, which should be less
1912 statistically noisy in most cases. This behavior can be overwritten
1913 by using the mdp option {\tt sc-coul} to 'yes'. Additionally, the
1914 soft-core interaction potential is only applied when either the A or B
1915 state has zero interaction potential. If both A and B states have
1916 nonzero interaction potential, default linear scaling described above
1917 is used. When both Coulombic and Lennard-Jones interactions are turned
1918 off simultaneously, a soft-core potential is used, and a hydrogen is
1919 being introduced or deleted, the sigma is set to {\tt sc-sigma-min},
1920 which itself defaults to {\tt sc-sigma-default}.
1922 Recently, a new formulation of the soft-core approach has been derived
1923 that in most cases gives lower and more even statistical variance than
1924 the standard soft-core path described above.~\cite{Pham2011,Pham2012}
1925 Specifically, we have:
1927 V_{sc}(r) &=& \LL V^A(r_A) + \LAM V^B(r_B)
1929 r_A &=& \left(\alpha \sigma_A^{48} \LAM^p + r^{48} \right)^\frac{1}{48}
1931 r_B &=& \left(\alpha \sigma_B^{48} \LL^p + r^{48} \right)^\frac{1}{48}
1933 This ``1-1-48'' path is also implemented in {\gromacs}. Note that for this path the soft core $\alpha$
1934 should satisfy $0.001 < \alpha < 0.003$,rather than $\alpha \approx
1937 %} % Brace matches ifthenelse test for gmxlite
1939 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1941 \subsection{Exclusions and 1-4 Interactions.}
1942 Atoms within a molecule that are close by in the chain,
1943 {\ie} atoms that are covalently bonded, or linked by one or two
1944 atoms are called {\em first neighbors, second neighbors} and
1945 {\em \swapindex{third}{neighbor}s}, respectively (see \figref{chain}).
1946 Since the interactions of atom {\bf i} with atoms {\bf i+1} and {\bf i+2}
1947 are mainly quantum mechanical, they can not be modeled by a Lennard-Jones potential.
1948 Instead it is assumed that these interactions are adequately modeled
1949 by a harmonic bond term or constraint ({\bf i, i+1}) and a harmonic angle term
1950 ({\bf i, i+2}). The first and second neighbors (atoms {\bf i+1} and {\bf i+2})
1952 {\em excluded} from the Lennard-Jones \swapindex{interaction}{list}
1954 atoms {\bf i+1} and {\bf i+2} are called {\em \normindex{exclusions}} of atom {\bf i}.
1957 \centerline{\includegraphics[angle=270,width=8cm]{plots/chain}}
1958 \caption{Atoms along an alkane chain.}
1962 For third neighbors, the normal Lennard-Jones repulsion is sometimes
1963 still too strong, which means that when applied to a molecule, the
1964 molecule would deform or break due to the internal strain. This is
1965 especially the case for carbon-carbon interactions in a {\em
1966 cis}-conformation ({\eg} {\em cis}-butane). Therefore, for some of these
1967 interactions, the Lennard-Jones repulsion has been reduced in the
1968 {\gromos} force field, which is implemented by keeping a separate list of
1969 1-4 and normal Lennard-Jones parameters. In other force fields, such
1970 as OPLS~\cite{Jorgensen88}, the standard Lennard-Jones parameters are reduced
1971 by a factor of two, but in that case also the dispersion (r$^{-6}$)
1972 and the Coulomb interaction are scaled.
1973 {\gromacs} can use either of these methods.
1975 \subsection{Charge Groups\index{charge group}}
1977 In principle, the force calculation in MD is an $O(N^2)$ problem.
1978 Therefore, we apply a \normindex{cut-off} for non-bonded force (NBF)
1979 calculations; only the particles within a certain distance of each
1980 other are interacting. This reduces the cost to $O(N)$ (typically
1981 $100N$ to $200N$) of the NBF. It also introduces an error, which is,
1982 in most cases, acceptable, except when applying the cut-off implies
1983 the creation of charges, in which case you should consider using the
1984 lattice sum methods provided by {\gromacs}.
1986 Consider a water molecule interacting with another atom. When we would apply
1987 the cut-off on an atom-atom basis we might include the atom-oxygen
1988 interaction (with a charge of $-0.82$) without the compensating charge
1989 of the protons, and as a result, induce a large dipole moment over the system.
1990 Therefore, we have to keep groups of atoms with total charge
1991 0 together. These groups are called {\em charge groups}.
1993 \subsection{Treatment of Cut-offs\index{cut-off}}
1994 \newcommand{\rs}{$r_{short}$}
1995 \newcommand{\rl}{$r_{long}$}
1996 {\gromacs} is quite flexible in treating cut-offs, which implies
1997 there can be quite a number of parameters to set. These parameters are
1998 set in the input file for {\tt grompp}. There are two sort of parameters
1999 that affect the cut-off interactions; you can select which type
2000 of interaction to use in each case, and which cut-offs should be
2001 used in the neighbor searching.
2003 For both Coulomb and van der Waals interactions there are interaction
2004 type selectors (termed {\tt vdwtype} and {\tt coulombtype}) and two
2005 parameters, for a total of six non-bonded interaction parameters. See
2006 \secref{mdpopt} for a complete description of these parameters.
2008 The neighbor searching (NS) can be performed using a single-range, or a twin-range
2009 approach. Since the former is merely a special case of the latter, we will
2010 discuss the more general twin-range. In this case, NS is described by two
2011 radii: {\tt rlist} and max({\tt rcoulomb},{\tt rvdw}).
2012 Usually one builds the neighbor list every 10 time steps
2013 or every 20 fs (parameter {\tt nstlist}). In the neighbor list, all interaction
2014 pairs that fall within {\tt rlist} are stored. Furthermore, the
2015 interactions between pairs that do not
2016 fall within {\tt rlist} but do fall within max({\tt rcoulomb},{\tt rvdw})
2017 are computed during NS. The
2018 forces and energy are stored separately and added to short-range forces
2019 at every time step between successive NS. If {\tt rlist} =
2020 max({\tt rcoulomb},{\tt rvdw}), no forces
2021 are evaluated during neighbor list generation.
2022 The \normindex{virial} is calculated from the sum of the short- and
2024 This means that the virial can be slightly asymmetrical at non-NS steps.
2025 In single precision, the virial is almost always asymmetrical because the
2026 off-diagonal elements are about as large as each element in the sum.
2027 In most cases this is not really a problem, since the fluctuations in the
2028 virial can be 2 orders of magnitude larger than the average.
2030 Except for the plain cut-off,
2031 all of the interaction functions in \tabref{funcparm}
2032 require that neighbor searching be done with a larger radius than the $r_c$
2033 specified for the functional form, because of the use of charge groups.
2034 The extra radius is typically of the order of 0.25 nm (roughly the
2035 largest distance between two atoms in a charge group plus the distance a
2036 charge group can diffuse within neighbor list updates).
2038 %If your charge groups are very large it may be interesting to turn off charge
2039 %groups, by setting the option
2040 %{\tt bAtomList = yes} in your {\tt grompp.mdp} file.
2041 %In this case only a small extra radius to account for diffusion needs to be
2042 %added (0.1 nm). Do not however use this together with the plain cut-off
2043 %method, as it will generate large artifacts (\secref{cg}).
2044 %In summary, there are four parameters that describe NS behavior:
2045 %{\tt nstlist} (update frequency in number of time steps),
2046 %{\tt bAtomList} (whether or not charge groups are used to generate neighbor list, the default is to use charge groups, so {\tt bAtomList = no}),
2047 %{\tt rshort} and {\tt rlong} which are the two radii {\rs} and {\rl}
2052 \begin{tabular}{|ll|l|}
2054 \multicolumn{2}{|c|}{Type} & Parameters \\
2056 Coulomb&Plain cut-off & $r_c$, $\epsr$ \\
2057 &Reaction field & $r_c$, $\epsrf$ \\
2058 &Shift function & $r_1$, $r_c$, $\epsr$ \\
2059 &Switch function & $r_1$, $r_c$, $\epsr$ \\
2061 VdW&Plain cut-off & $r_c$ \\
2062 &Shift function & $r_1$, $r_c$ \\
2063 &Switch function & $r_1$, $r_c$ \\
2066 \caption[Parameters for the different functional forms of the
2067 non-bonded interactions.]{Parameters for the different functional
2068 forms of the non-bonded interactions.}
2069 \label{tab:funcparm}
2071 %} % Brace matches ifthenelse test for gmxlite
2074 \newcommand{\vvis}{\ve{r}_s}
2075 \newcommand{\Fi}{\ve{F}_i'}
2076 \newcommand{\Fj}{\ve{F}_j'}
2077 \newcommand{\Fk}{\ve{F}_k'}
2078 \newcommand{\Fl}{\ve{F}_l'}
2079 \newcommand{\Fvis}{\ve{F}_{s}}
2080 \newcommand{\rvik}{\ve{r}_{ik}}
2081 \newcommand{\rvis}{\ve{r}_{is}}
2082 \newcommand{\rvjk}{\ve{r}_{jk}}
2083 \newcommand{\rvjl}{\ve{r}_{jl}}
2085 %\ifthenelse{\equal{\gmxlite}{1}}{}{
2086 \section{Virtual interaction sites\index{virtual interaction sites}}
2087 \label{sec:virtual_sites}
2088 Virtual interaction sites (called \seeindex{dummy atoms}{virtual interaction sites} in {\gromacs} versions before 3.3)
2089 can be used in {\gromacs} in a number of ways.
2090 We write the position of the virtual site $\ve{r}_s$ as a function of
2091 the positions of other particles \ve{r}$_i$: $\ve{r}_s =
2092 f(\ve{r}_1..\ve{r}_n)$. The virtual site, which may carry charge or be
2093 involved in other interactions, can now be used in the force
2094 calculation. The force acting on the virtual site must be
2095 redistributed over the particles with mass in a consistent way.
2096 A good way to do this can be found in ref.~\cite{Berendsen84b}.
2097 We can write the potential energy as:
2099 V = V(\vvis,\ve{r}_1,\ldots,\ve{r}_n) = V^*(\ve{r}_1,\ldots,\ve{r}_n)
2101 The force on the particle $i$ is then:
2103 \ve{F}_i = -\frac{\partial V^*}{\partial \ve{r}_i}
2104 = -\frac{\partial V}{\partial \ve{r}_i} -
2105 \frac{\partial V}{\partial \vvis}
2106 \frac{\partial \vvis}{\partial \ve{r}_i}
2107 = \ve{F}_i^{direct} + \Fi
2109 The first term is the normal force.
2110 The second term is the force on particle $i$ due to the virtual site, which
2111 can be written in tensor notation:
2112 \newcommand{\partd}[2]{\displaystyle\frac{\partial #1}{\partial #2_i}}
2114 \Fi = \left[\begin{array}{ccc}
2115 \partd{x_s}{x} & \partd{y_s}{x} & \partd{z_s}{x} \\[1ex]
2116 \partd{x_s}{y} & \partd{y_s}{y} & \partd{z_s}{y} \\[1ex]
2117 \partd{x_s}{z} & \partd{y_s}{z} & \partd{z_s}{z}
2118 \end{array}\right]\Fvis
2121 where $\Fvis$ is the force on the virtual site and $x_s$, $y_s$ and
2122 $z_s$ are the coordinates of the virtual site. In this way, the total
2123 force and the total torque are conserved~\cite{Berendsen84b}.
2125 The computation of the \normindex{virial}
2126 (\eqnref{Xi}) is non-trivial when virtual sites are used. Since the
2127 virial involves a summation over all the atoms (rather than virtual
2128 sites), the forces must be redistributed from the virtual sites to the
2129 atoms (using ~\eqnref{fvsite}) {\em before} computation of the
2130 virial. In some special cases where the forces on the atoms can be
2131 written as a linear combination of the forces on the virtual sites (types 2
2132 and 3 below) there is no difference between computing the virial
2133 before and after the redistribution of forces. However, in the
2134 general case redistribution should be done first.
2137 \centerline{\includegraphics[width=15cm]{plots/dummies}}
2138 \caption[Virtual site construction.]{The six different types of virtual
2139 site construction in \protect{\gromacs}. The constructing atoms are
2140 shown as black circles, the virtual sites in gray.}
2144 There are six ways to construct virtual sites from surrounding atoms in
2145 {\gromacs}, which we classify by the number of constructing
2146 atoms. {\bf Note} that all site types mentioned can be constructed from
2147 types 3fd (normalized, in-plane) and 3out (non-normalized, out of
2148 plane). However, the amount of computation involved increases sharply
2149 along this list, so we strongly recommended using the first adequate
2150 virtual site type that will be sufficient for a certain purpose.
2151 \figref{vsites} depicts 6 of the available virtual site constructions.
2152 The conceptually simplest construction types are linear combinations:
2154 \vvis = \sum_{i=1}^N w_i \, \ve{r}_i
2156 The force is then redistributed using the same weights:
2161 The types of virtual sites supported in {\gromacs} are given in the list below.
2162 Constructing atoms in virtual sites can be virtual sites themselves, but
2163 only if they are higher in the list, i.e. virtual sites can be
2164 constructed from ``particles'' that are simpler virtual sites.
2166 \item[{\bf\sf 2.}]\label{subsec:vsite2}As a linear combination of two atoms
2167 (\figref{vsites} 2):
2169 w_i = 1 - a ~,~~ w_j = a
2171 In this case the virtual site is on the line through atoms $i$ and
2174 \item[{\bf\sf 3.}]\label{subsec:vsite3}As a linear combination of three atoms
2175 (\figref{vsites} 3):
2177 w_i = 1 - a - b ~,~~ w_j = a ~,~~ w_k = b
2179 In this case the virtual site is in the plane of the other three
2182 \item[{\bf\sf 3fd.}]\label{subsec:vsite3fd}In the plane of three atoms, with a fixed distance
2183 (\figref{vsites} 3fd):
2185 \vvis ~=~ \ve{r}_i + b \frac{ \rvij + a \rvjk }
2186 {| \rvij + a \rvjk |}
2188 In this case the virtual site is in the plane of the other three
2189 particles at a distance of $|b|$ from $i$.
2190 The force on particles $i$, $j$ and $k$ due to the force on the virtual
2191 site can be computed as:
2194 \Fi &=& \displaystyle \Fvis - \gamma ( \Fvis - \ve{p} ) \\[1ex]
2195 \Fj &=& \displaystyle (1-a)\gamma (\Fvis - \ve{p}) \\[1ex]
2196 \Fk &=& \displaystyle a \gamma (\Fvis - \ve{p}) \\
2200 \displaystyle \gamma = \frac{b}{| \rvij + a \rvjk |} \\[2ex]
2201 \displaystyle \ve{p} = \frac{ \rvis \cdot \Fvis }
2202 { \rvis \cdot \rvis } \rvis
2206 \item[{\bf\sf 3fad.}]\label{subsec:vsite3fad}In the plane of three atoms, with a fixed angle and
2207 distance (\figref{vsites} 3fad):
2209 \label{eqn:vsite2fad-F}
2210 \vvis ~=~ \ve{r}_i +
2211 d \cos \theta \frac{\rvij}{|\rvij|} +
2212 d \sin \theta \frac{\ve{r}_\perp}{|\ve{r}_\perp|}
2214 \ve{r}_\perp ~=~ \rvjk -
2215 \frac{ \rvij \cdot \rvjk }
2216 { \rvij \cdot \rvij }
2219 In this case the virtual site is in the plane of the other three
2220 particles at a distance of $|d|$ from $i$ at an angle of
2221 $\alpha$ with $\rvij$. Atom $k$ defines the plane and the
2222 direction of the angle. {\bf Note} that in this case $b$ and
2223 $\alpha$ must be specified, instead of $a$ and $b$ (see also
2224 \secref{vsitetop}). The force on particles $i$, $j$ and $k$
2225 due to the force on the virtual site can be computed as (with
2226 $\ve{r}_\perp$ as defined in \eqnref{vsite2fad-F}):
2227 \newcommand{\dfrac}{\displaystyle\frac}
2230 \begin{array}{lclllll}
2232 \dfrac{d \cos \theta}{|\rvij|} \ve{F}_1 &+&
2233 \dfrac{d \sin \theta}{|\ve{r}_\perp|} \left(
2234 \dfrac{ \rvij \cdot \rvjk }
2235 { \rvij \cdot \rvij } \ve{F}_2 +
2236 \ve{F}_3 \right) \\[3ex]
2238 \dfrac{d \cos \theta}{|\rvij|} \ve{F}_1 &-&
2239 \dfrac{d \sin \theta}{|\ve{r}_\perp|} \left(
2241 \dfrac{ \rvij \cdot \rvjk }
2242 { \rvij \cdot \rvij } \ve{F}_2 +
2243 \ve{F}_3 \right) \\[3ex]
2245 \dfrac{d \sin \theta}{|\ve{r}_\perp|} \ve{F}_2 \\[3ex]
2249 \dfrac{ \rvij \cdot \Fvis }
2250 { \rvij \cdot \rvij } \rvij
2252 \ve{F}_2 = \ve{F}_1 -
2253 \dfrac{ \ve{r}_\perp \cdot \Fvis }
2254 { \ve{r}_\perp \cdot \ve{r}_\perp } \ve{r}_\perp
2256 \ve{F}_3 = \dfrac{ \rvij \cdot \Fvis }
2257 { \rvij \cdot \rvij } \ve{r}_\perp
2261 \item[{\bf\sf 3out.}]\label{subsec:vsite3out}As a non-linear combination of three atoms, out of plane
2262 (\figref{vsites} 3out):
2264 \vvis ~=~ \ve{r}_i + a \rvij + b \rvik +
2265 c (\rvij \times \rvik)
2267 This enables the construction of virtual sites out of the plane of the
2269 The force on particles $i,j$ and $k$ due to the force on the virtual
2270 site can be computed as:
2274 \Fj &=& \left[\begin{array}{ccc}
2275 a & -c\,z_{ik} & c\,y_{ik} \\[0.5ex]
2276 c\,z_{ik} & a & -c\,x_{ik} \\[0.5ex]
2277 -c\,y_{ik} & c\,x_{ik} & a
2278 \end{array}\right]\Fvis \\
2280 \Fk &=& \left[\begin{array}{ccc}
2281 b & c\,z_{ij} & -c\,y_{ij} \\[0.5ex]
2282 -c\,z_{ij} & b & c\,x_{ij} \\[0.5ex]
2283 c\,y_{ij} & -c\,x_{ij} & b
2284 \end{array}\right]\Fvis \\
2285 \Fi &=& \Fvis - \Fj - \Fk
2289 \item[{\bf\sf 4fdn.}]\label{subsec:vsite4fdn}From four atoms, with a fixed distance, see separate Fig.\ \ref{fig:vsite-4fdn}.
2290 This construction is a bit
2291 complex, in particular since the previous type (4fd) could be unstable which forced us
2292 to introduce a more elaborate construction:
2295 \centerline{\includegraphics[width=5cm]{plots/vsite-4fdn}}
2296 \caption {The new 4fdn virtual site construction, which is stable even when all constructing
2297 atoms are in the same plane.}
2298 \label{fig:vsite-4fdn}
2302 \mathbf{r}_{ja} &=& a\, \mathbf{r}_{ik} - \mathbf{r}_{ij} = a\, (\mathbf{x}_k - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\
2303 \mathbf{r}_{jb} &=& b\, \mathbf{r}_{il} - \mathbf{r}_{ij} = b\, (\mathbf{x}_l - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\
2304 \mathbf{r}_m &=& \mathbf{r}_{ja} \times \mathbf{r}_{jb} \nonumber \\
2305 \mathbf{x}_s &=& \mathbf{x}_i + c \frac{\mathbf{r}_m}{|\mathbf{r}_m|}
2309 In this case the virtual site is at a distance of $|c|$ from $i$, while $a$ and $b$ are
2310 parameters. {\bf Note} that the vectors $\mathbf{r}_{ik}$ and $\mathbf{r}_{ij}$ are not normalized
2311 to save floating-point operations.
2312 The force on particles $i$, $j$, $k$ and $l$ due to the force
2313 on the virtual site are computed through chain rule derivatives
2314 of the construction expression. This is exact and conserves energy,
2315 but it does lead to relatively lengthy expressions that we do not
2316 include here (over 200 floating-point operations). The interested reader can
2317 look at the source code in \verb+vsite.c+. Fortunately, this vsite type is normally
2318 only used for chiral centers such as $C_{\alpha}$ atoms in proteins.
2320 The new 4fdn construct is identified with a `type' value of 2 in the topology. The earlier 4fd
2321 type is still supported internally (`type' value 1), but it should not be used for
2322 new simulations. All current {\gromacs} tools will automatically generate type 4fdn instead.
2325 \item[{\bf\sf N.}]\label{subsec:vsiteN} A linear combination of $N$ atoms with relative
2326 weights $a_i$. The weight for atom $i$ is:
2328 w_i = a_i \left(\sum_{j=1}^N a_j \right)^{-1}
2330 There are three options for setting the weights:
2332 \item[COG] center of geometry: equal weights
2333 \item[COM] center of mass: $a_i$ is the mass of atom $i$;
2334 when in free-energy simulations the mass of the atom is changed,
2335 only the mass of the A-state is used for the weight
2336 \item[COW] center of weights: $a_i$ is defined by the user
2340 %} % Brace matches ifthenelse test for gmxlite
2342 \newcommand{\dr}{{\rm d}r}
2343 \newcommand{\avcsix}{\left< C_6 \right>}
2345 %\ifthenelse{\equal{\gmxlite}{1}}{}{
2346 \section{Long Range Electrostatics}
2347 \label{sec:lr_elstat}
2348 \subsection{Ewald summation\index{Ewald sum}}
2350 The total electrostatic energy of $N$ particles and their periodic
2351 images\index{periodic boundary conditions} is given by
2353 V=\frac{f}{2}\sum_{n_x}\sum_{n_y}
2354 \sum_{n_{z}*} \sum_{i}^{N} \sum_{j}^{N}
2355 \frac{q_i q_j}{{\bf r}_{ij,{\bf n}}}.
2356 \label{eqn:totalcoulomb}
2358 $(n_x,n_y,n_z)={\bf n}$ is the box index vector, and the star indicates that
2359 terms with $i=j$ should be omitted when $(n_x,n_y,n_z)=(0,0,0)$. The
2360 distance ${\bf r}_{ij,{\bf n}}$ is the real distance between the charges and
2361 not the minimum-image. This sum is conditionally convergent, but
2364 Ewald summation was first introduced as a method to calculate
2365 long-range interactions of the periodic images in
2366 crystals~\cite{Ewald21}. The idea is to convert the single
2367 slowly-converging sum \eqnref{totalcoulomb} into two
2368 quickly-converging terms and a constant term:
2370 V &=& V_{dir} + V_{rec} + V_{0} \\[0.5ex]
2371 V_{dir} &=& \frac{f}{2} \sum_{i,j}^{N}
2372 \sum_{n_x}\sum_{n_y}
2373 \sum_{n_{z}*} q_i q_j \frac{\mbox{erfc}(\beta {r}_{ij,{\bf n}} )}{{r}_{ij,{\bf n}}} \\[0.5ex]
2374 V_{rec} &=& \frac{f}{2 \pi V} \sum_{i,j}^{N} q_i q_j
2375 \sum_{m_x}\sum_{m_y}
2376 \sum_{m_{z}*} \frac{\exp{\left( -(\pi {\bf m}/\beta)^2 + 2 \pi i
2377 {\bf m} \cdot ({\bf r}_i - {\bf r}_j)\right)}}{{\bf m}^2} \\[0.5ex]
2378 V_{0} &=& -\frac{f \beta}{\sqrt{\pi}}\sum_{i}^{N} q_i^2,
2380 where $\beta$ is a parameter that determines the relative weight of the
2381 direct and reciprocal sums and ${\bf m}=(m_x,m_y,m_z)$.
2382 In this way we can use a short cut-off (of the order of $1$~nm) in the direct space sum and a
2383 short cut-off in the reciprocal space sum ({\eg} 10 wave vectors in each
2384 direction). Unfortunately, the computational cost of the reciprocal
2385 part of the sum increases as $N^2$
2386 (or $N^{3/2}$ with a slightly better algorithm) and it is therefore not
2387 realistic for use in large systems.
2389 \subsubsection{Using Ewald}
2390 Don't use Ewald unless you are absolutely sure this is what you want -
2391 for almost all cases the PME method below will perform much better.
2392 If you still want to employ classical Ewald summation enter this in
2393 your {\tt .mdp} file, if the side of your box is about $3$~nm:
2400 fourierspacing = 0.6
2404 The ratio of the box dimensions and the {\tt fourierspacing} parameter determines
2405 the highest magnitude of wave vectors $m_x,m_y,m_z$ to use in each
2406 direction. With a 3-nm cubic box this example would use $11$ wave vectors
2407 (from $-5$ to $5$) in each direction. The {\tt ewald-rtol} parameter
2408 is the relative strength of the electrostatic interaction at the
2409 cut-off. Decreasing this gives you a more accurate direct sum, but a
2410 less accurate reciprocal sum.
2412 \subsection{\normindex{PME}}
2414 Particle-mesh Ewald is a method proposed by Tom
2415 Darden~\cite{Darden93} to improve the performance of the
2416 reciprocal sum. Instead of directly summing wave vectors, the charges
2417 are assigned to a grid using interpolation. The implementation in
2418 {\gromacs} uses cardinal B-spline interpolation~\cite{Essmann95},
2419 which is referred to as smooth PME (SPME).
2420 The grid is then Fourier transformed with a 3D FFT algorithm and the
2421 reciprocal energy term obtained by a single sum over the grid in
2424 The potential at the grid points is calculated by inverse
2425 transformation, and by using the interpolation factors we get the
2426 forces on each atom.
2428 The PME algorithm scales as $N \log(N)$, and is substantially faster
2429 than ordinary Ewald summation on medium to large systems. On very
2430 small systems it might still be better to use Ewald to avoid the
2431 overhead in setting up grids and transforms.
2432 For the parallelization of PME see the section on MPMD PME (\ssecref{mpmd_pme}).
2434 With the Verlet cut-off scheme, the PME direct space potential is
2435 shifted by a constant such that the potential is zero at the
2436 cut-off. This shift is small and since the net system charge is close
2437 to zero, the total shift is very small, unlike in the case of the
2438 Lennard-Jones potential where all shifts add up. We apply the shift
2439 anyhow, such that the potential is the exact integral of the force.
2441 \subsubsection{Using PME}
2442 To use Particle-mesh Ewald summation in {\gromacs}, specify the
2443 following lines in your {\tt .mdp} file:
2450 fourierspacing = 0.12
2455 In this case the {\tt fourierspacing} parameter determines the maximum
2456 spacing for the FFT grid (i.e. minimum number of grid points),
2457 and {\tt pme-order} controls the
2458 interpolation order. Using fourth-order (cubic) interpolation and this
2459 spacing should give electrostatic energies accurate to about
2460 $5\cdot10^{-3}$. Since the Lennard-Jones energies are not this
2461 accurate it might even be possible to increase this spacing slightly.
2463 Pressure scaling works with PME, but be aware of the fact that
2464 anisotropic scaling can introduce artificial ordering in some systems.
2466 \subsection{\normindex{P3M-AD}}
2468 The \seeindex{Particle-Particle Particle-Mesh}{P3M} methods of
2469 Hockney \& Eastwood can also be applied in {\gromacs} for the
2470 treatment of long range electrostatic interactions~\cite{Hockney81}.
2471 Although the P3M method was the first efficient long-range electrostatics
2472 method for molecular simulation, the smooth PME (SPME) method has largely
2473 replaced P3M as the method of choice in atomistic simulations. One performance
2474 disadvantage of the original P3M method was that it required 3 3D-FFT
2475 back transforms to obtain the forces on the particles. But this is not
2476 required for P3M and the forces can be derived through analytical differentiation
2477 of the potential, as done in PME. The resulting method is termed P3M-AD.
2478 The only remaining difference between P3M-AD and PME is the optimization
2479 of the lattice Green influence function for error minimization that P3M uses.
2480 However, in 2012 it has been shown that the SPME influence function can be
2481 modified to obtain P3M~\cite{Ballenegger2012}.
2482 This means that the advantage of error minimization in P3M-AD can be used
2483 at the same computational cost and with the same code as PME,
2484 just by adding a few lines to modify the influence function.
2485 However, at optimal parameter setting the effect of error minimization
2486 in P3M-AD is less than 10\%. P3M-AD does show large accuracy gains with
2487 interlaced (also known as staggered) grids, but that is not supported
2488 in {\gromacs} (yet).
2490 P3M is used in {\gromacs} with exactly the same options as used with PME
2491 by selecting the electrostatics type:
2493 coulombtype = P3M-AD
2496 \subsection{Optimizing Fourier transforms}
2497 To get the best possible performance you should try to avoid large
2498 prime numbers for grid dimensions.
2499 The FFT code used in {\gromacs} is
2500 optimized for grid sizes of the form $2^a 3^b 5^c 7^d 11^e 13^f$,
2501 where $e+f$ is $0$ or $1$ and the other exponents arbitrary. (See
2502 further the documentation of the FFT algorithms at
2503 \href{http://www.fftw.org}{www.fftw.org}.
2505 It is also possible to optimize the transforms for the current problem
2506 by performing some calculations at the start of the run. This is not
2507 done by default since it takes a couple of minutes, but for large
2508 runs it will save time. Turn it on by specifying
2513 in your {\tt .mdp} file.
2515 When running in parallel, the grid must be communicated several times,
2516 thus hurting scaling performance. With PME you can improve this
2517 by increasing grid spacing while simultaneously increasing the
2518 interpolation to {\eg} sixth order.
2519 Since the interpolation is entirely local, doing so will
2520 improve the scaling in most cases.
2523 % Temporarily removed since I am not sure about the state of the testlr
2526 %It is possible to test the accuracy of your settings using the program
2527 %{\tt\normindex{testlr}} in the {\tt src/gmxlib} dir. This program computes
2528 %forces and potentials using PPPM and an Ewald implementation and gives the
2529 %absolute and RMS errors in both:
2534 %Potential 0.113 0.035
2536 %{\bf Note:} these numbers were generated using a grid spacing of
2537 %0.058 nm and $r_c$ = 1.0 nm.
2539 %You can see what the accuracy is without optimizing the
2540 %$\hat{G}(k)$ function, if you pass the {\tt -ghat} option to {\tt
2541 %testlr}. Try it if you think the {\tt mk_ghat} procedure is a waste
2543 %} % Brace matches ifthenelse test for gmxlite
2546 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2547 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2548 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2550 %\ifthenelse{\equal{\gmxlite}{1}}{}{
2551 \section{Long Range Van der Waals interactions}
2552 \subsection{Dispersion correction\index{dispersion correction}}
2553 In this section, we derive long-range corrections due to the use of a
2554 cut-off for Lennard-Jones or Buckingham interactions.
2555 We assume that the cut-off is
2556 so long that the repulsion term can safely be neglected, and therefore
2557 only the dispersion term is taken into account. Due to the nature of
2558 the dispersion interaction (we are truncating a potential proportional
2559 to $-r^{-6}$), energy and pressure corrections are both negative. While
2560 the energy correction is usually small, it may be important for free
2561 energy calculations where differences between two different Hamiltonians
2562 are considered. In contrast, the pressure correction is very large and
2563 can not be neglected under any circumstances where a correct pressure is
2564 required, especially for any NPT simulations. Although it is, in
2565 principle, possible to parameterize a force field such that the pressure
2566 is close to the desired experimental value without correction, such a
2567 method makes the parameterization dependent on the cut-off and is therefore
2570 \subsubsection{Energy}
2572 The long-range contribution of the dispersion interaction to the
2573 virial can be derived analytically, if we assume a homogeneous
2574 system beyond the cut-off distance $r_c$. The dispersion energy
2575 between two particles is written as:
2577 V(\rij) ~=~- C_6\,\rij^{-6}
2579 and the corresponding force is:
2581 \Fvij ~=~- 6\,C_6\,\rij^{-8}\rvij
2583 In a periodic system it is not easy to calculate the full potentials,
2584 so usually a cut-off is applied, which can be abrupt or smooth.
2585 We will call the potential and force with cut-off $V_c$ and $\ve{F}_c$.
2586 The long-range contribution to the dispersion energy
2587 in a system with $N$ particles and particle density $\rho$ = $N/V$ is:
2589 \label{eqn:enercorr}
2590 V_{lr} ~=~ \half N \rho\int_0^{\infty} 4\pi r^2 g(r) \left( V(r) -V_c(r) \right) {\dr}
2592 We will integrate this for the shift function, which is the most general
2593 form of van der Waals interaction available in {\gromacs}.
2594 The shift function has a constant difference $S$ from 0 to $r_1$
2595 and is 0 beyond the cut-off distance $r_c$.
2596 We can integrate \eqnref{enercorr}, assuming that the density in the sphere
2597 within $r_1$ is equal to the global density and
2598 the radial distribution function $g(r)$ is 1 beyond $r_1$:
2601 V_{lr} &=& \half N \left(
2602 \rho\int_0^{r_1} 4\pi r^2 g(r) \, C_6 \,S\,{\dr}
2603 + \rho\int_{r_1}^{r_c} 4\pi r^2 \left( V(r) -V_c(r) \right) {\dr}
2604 + \rho\int_{r_c}^{\infty} 4\pi r^2 V(r) \, {\dr}
2606 & = & \half N \left(\left(\frac{4}{3}\pi \rho r_1^{3} - 1\right) C_6 \,S
2607 + \rho\int_{r_1}^{r_c} 4\pi r^2 \left( V(r) -V_c(r) \right) {\dr}
2608 -\frac{4}{3} \pi N \rho\, C_6\,r_c^{-3}
2611 where the term $-1$ corrects for the self-interaction.
2612 For a plain cut-off we only need to assume that $g(r)$ is 1 beyond $r_c$
2613 and the correction reduces to~\cite{Allen87}:
2615 V_{lr} & = & -\frac{2}{3} \pi N \rho\, C_6\,r_c^{-3}
2617 If we consider, for example, a box of pure water, simulated with a cut-off
2618 of 0.9 nm and a density of 1 g cm$^{-3}$ this correction is
2619 $-0.75$ kJ mol$^{-1}$ per molecule.
2621 For a homogeneous mixture we need to define
2622 an {\em average dispersion constant}:
2625 \avcsix = \frac{2}{N(N-1)}\sum_i^N\sum_{j>i}^N C_6(i,j)\\
2627 In {\gromacs}, excluded pairs of atoms do not contribute to the average.
2629 In the case of inhomogeneous simulation systems, {\eg} a system with a
2630 lipid interface, the energy correction can be applied if
2631 $\avcsix$ for both components is comparable.
2633 \subsubsection{Virial and pressure}
2634 The scalar virial of the system due to the dispersion interaction between
2635 two particles $i$ and $j$ is given by:
2637 \Xi~=~-\half \rvij \cdot \Fvij ~=~ 3\,C_6\,\rij^{-6}
2639 The pressure is given by:
2641 P~=~\frac{2}{3\,V}\left(E_{kin} - \Xi\right)
2643 The long-range correction to the virial is given by:
2645 \Xi_{lr} ~=~ \half N \rho \int_0^{\infty} 4\pi r^2 g(r) (\Xi -\Xi_c) \,\dr
2647 We can again integrate the long-range contribution to the
2648 virial assuming $g(r)$ is 1 beyond $r_1$:
2650 \Xi_{lr}&=& \half N \rho \left(
2651 \int_{r_1}^{r_c} 4 \pi r^2 (\Xi -\Xi_c) \,\dr
2652 + \int_{r_c}^{\infty} 4 \pi r^2 3\,C_6\,\rij^{-6}\, \dr
2654 &=& \half N \rho \left(
2655 \int_{r_1}^{r_c} 4 \pi r^2 (\Xi -\Xi_c) \, \dr
2656 + 4 \pi C_6 \, r_c^{-3} \right)
2658 For a plain cut-off the correction to the pressure is~\cite{Allen87}:
2660 P_{lr}~=~-\frac{4}{3} \pi C_6\, \rho^2 r_c^{-3}
2662 Using the same example of a water box, the correction to the virial is
2663 0.75 kJ mol$^{-1}$ per molecule,
2664 the corresponding correction to the pressure for
2665 SPC water is approximately $-280$ bar.
2667 For homogeneous mixtures, we can again use the average dispersion constant
2668 $\avcsix$ (\eqnref{avcsix}):
2670 P_{lr}~=~-\frac{4}{3} \pi \avcsix \rho^2 r_c^{-3}
2673 For inhomogeneous systems, \eqnref{pcorr} can be applied under the same
2674 restriction as holds for the energy (see \secref{ecorr}).
2676 \subsection{Lennard-Jones PME\index{LJ-PME}}
2678 In order to treat systems, using Lennard-Jones potentials, that are
2679 non-homogeneous outside of the cut-off distance, we can instead use
2680 the Particle-mesh Ewald method as discussed for electrostatics above.
2681 In this case the modified Ewald equations become
2683 V &=& V_{dir} + V_{rec} + V_{0} \\[0.5ex]
2684 V_{dir} &=& -\frac{1}{2} \sum_{i,j}^{N}
2685 \sum_{n_x}\sum_{n_y}
2686 \sum_{n_{z}*} \frac{C_{ij}^{(6)}g(\beta {r}_{ij,{\bf n}})}{{r_{ij,{\bf n}}}^6} \\[0.5ex]
2687 V_{rec} &=& \frac{{\pi}^{\frac{3}{2}} \beta^{3}}{2V} \sum_{m_x}\sum_{m_y}\sum_{m_{z}*}
2688 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]
2689 V_{0} &=& -\frac{\beta^{6}}{12}\sum_{i}^{N} C_{ii}^{(6)},
2692 where ${\bf m}=(m_x,m_y,m_z)$, $\beta$ is the parameter determining the weight between
2693 direct and reciprocal space, and ${C_{ij}}^{(6)}$ is the combined dispersion
2694 parameter for particle $i$ and $j$. The star indicates that terms
2695 with $i = j$ should be omitted when $((n_x,n_y,n_z)=(0,0,0))$, and
2696 ${\bf r}_{ij,{\bf n}}$ is the real distance between the particles.
2697 Following the derivation by Essmann~\cite{Essmann95}, the functions $f$ and $g$ introduced above are defined as
2699 f(x)&=&1/3\left[(1-2x^2){\mathrm{exp}}(-x^2) + 2{x^3}\sqrt{\pi}\,{\mathrm{erfc}}(x) \right] \\
2700 g(x)&=&{\mathrm{exp}}(-x^2)(1+x^2+\frac{x^4}{2}).
2703 The above methodology works fine as long as the dispersion parameters can be factorized in the same
2704 way as the charges for electrostatics
2706 C_{ij}^{(6)} = \left({C_{ii}^{(6)} \, C_{jj}^{(6)}}\right)^{1/2}
2708 For Lorentz-Berthelot combination rules, the reciprocal part of this sum has to be calculated
2709 seven times due to the splitting of the dispersion parameter according to
2711 C_{ij}^{(6)}=(\sigma_i+\sigma_j)^6=\sum_{n=0}^{6} P_{n}\sigma_{i}^{n}\sigma_{j}^{(6-n)},
2713 for $P_{n}$ the Pascal triangle coefficients. This introduces a
2714 non-negligible cost to the reciprocal part, requiring seven separate
2715 FFTs, and therefore this have been the limiting factor in previous
2716 attempts to implement LJ-PME. A solution to this problem is to
2717 approximate the interaction parameters in reciprocal space using
2718 geometrical combination rules. This will preserve a well-defined
2719 Hamiltonian and significantly increase the performance of the
2723 \centerline{\includegraphics[width=15cm]{plots/ljpmedifference}}
2724 \caption {Dispersion potential between phosphorous and oxygen, The total, real and
2725 reciprocal parts of the total dispersion are shown. The reciprocal parts are calculated
2726 using either LB or geometric rules. The difference introduced by the use of the
2727 geometric approximation in the reciprocal part is small compared to the total
2728 interaction energy.}
2729 \label{fig:ljpmedifference}
2732 This approximation does introduce some errors, but since the
2733 difference is located in the interactions calculated in reciprocal
2734 space, the effect will be very small compared to the total interaction
2735 energy (see \figref{ljpmedifference}). The relative error in
2736 the total dispersion energy will stay below 0.5\% in a lipid bilayer
2737 simulation, when using a real space cut-off of 1.0 nm. A more
2738 thorough discussion of this can be found in \cite{Wennberg13}.
2740 \subsubsection{Using LJ-PME}
2741 To use Particle-mesh Ewald summation for Lennard-Jones interactions in {\gromacs}, specify the
2742 following lines in your {\tt .mdp} file:
2748 fourierspacing = 0.12
2750 ewald-rtol-lj = 0.001
2751 lj-pme-comb-rule = geometric
2754 The {\tt fourierspacing} and {\tt pme-order} are the same parameters
2755 as is used for electrostatic PME, and {\tt ewald-rtol-lj} controls
2756 splitting between real and reciprocal space in the same way as
2757 {\tt ewald-rtol}. In addition to this, the combination rule to be used
2758 in reciprocal space is determined by {\tt lj-pme-comb-rule}. If the
2759 current force field uses Lorentz-Berthelot combination rules, it is
2760 possible to set {\tt lj-pme-comb-rule = geometric} in order to gain a
2761 significant increase in performance for a small loss in accuracy. The
2762 details of this approximation can be found in the section above. The
2763 implementation of LJ-PME is currently unsupported together with the
2764 Verlet cut-off scheme and/or free energy calculations. These features
2765 will be added in upcoming releases
2766 %} % Brace matches ifthenelse test for gmxlite
2768 \section{Force field\index{force field}}
2770 A force field is built up from two distinct components:
2772 \item The set of equations (called the {\em
2773 \index{potential function}s}) used to generate the potential
2774 energies and their derivatives, the forces. These are described in
2775 detail in the previous chapter.
2776 \item The parameters used in this set of equations. These are not
2777 given in this manual, but in the data files corresponding to your
2778 {\gromacs} distribution.
2780 Within one set of equations various sets of parameters can be
2781 used. Care must be taken that the combination of equations and
2782 parameters form a consistent set. It is in general dangerous to make
2783 {\em ad hoc} changes in a subset of parameters, because the various
2784 contributions to the total force are usually interdependent. This
2785 means in principle that every change should be documented, verified by
2786 comparison to experimental data and published in a peer-reviewed
2787 journal before it can be used.
2789 {\gromacs} {\gmxver} includes several force fields, and additional
2790 ones are available on the website. If you do not know which one to
2791 select we recommend \gromosv{96} for united-atom setups and OPLS-AA/L for
2792 all-atom parameters. That said, we describe the available options in
2795 \subsection{GROMOS87\index{GROMOS87 force field}}
2796 The \gromosv{87} suite of programs and corresponding force
2797 field~\cite{biomos} formed the basis for the development of {\gromacs}
2798 in the early 1990s. The original GROMOS87 force field is not
2799 available in {\gromacs}. In previous versions ($<$ 3.3.2) there used
2800 to be the so-called ``{\gromacs} force field,'' which was based on
2801 \gromosv{87}~\cite{biomos}\index{GROMOS87}, with a small modification
2802 concerning the interaction between water oxygens and carbon
2803 atoms~\cite{Buuren93b,Mark94}, as well as 10 extra atom
2804 types~\cite{Jorgensen83,Buuren93a,Buuren93b,Mark94,Liu95}.
2806 Since version 5.0 this force field has been ``deprecated''. Should
2807 you have a justifiable reason to use this force field please
2808 use eariler versions of {\gromacs}.
2810 \subsubsection{All-hydrogen force field}
2811 The \gromosv{87}-based all-hydrogen force field is almost identical to the
2812 normal \gromosv{87} force field, since the extra hydrogens have no
2813 Lennard-Jones interaction and zero charge. The only differences are in
2814 the bond angle and improper dihedral angle terms. This force field is
2815 only useful when you need the exact hydrogen positions, for instance
2816 for distance restraints derived from NMR measurements. When citing
2817 this force field please read the previous paragraph.
2819 \subsection{\gromosv{96}\index{GROMOS96 force field}}
2820 {\gromacs} supports the \gromosv{96} force fields~\cite{gromos96}.
2821 All parameters for the 43A1, 43A2 (development, improved alkane
2822 dihedrals), 45A3, 53A5, and 53A6 parameter sets are included. All standard
2823 building blocks are included and topologies can be built automatically
2826 The \gromosv{96} force field is a further development of the \gromosv{87} force field.
2827 It has improvements over the \gromosv{87} force field for proteins and small molecules.
2828 {\bf Note} that the sugar parameters present in 53A6 do correspond to those published in
2829 2004\cite{Oostenbrink2004}, which are different from those present in 45A4, which
2830 is not included in {\gromacs} at this time. The 45A4 parameter set corresponds to a later
2831 revision of these parameters.
2832 The \gromosv{96} force field is not, however, recommended for use with long alkanes and
2833 lipids. The \gromosv{96} force field differs from the \gromosv{87}
2834 force field in a few respects:
2836 \item the force field parameters
2837 \item the parameters for the bonded interactions are not linked to atom types
2838 \item a fourth power bond stretching potential (\ssecref{G96bond})
2839 \item an angle potential based on the cosine of the angle (\ssecref{G96angle})
2841 There are two differences in implementation between {\gromacs} and \gromosv{96}
2842 which can lead to slightly different results when simulating the same system
2845 \item in \gromosv{96} neighbor searching for solvents is performed on the
2846 first atom of the solvent molecule. This is not implemented in {\gromacs},
2847 but the difference with searching by centers of charge groups is very small
2848 \item the virial in \gromosv{96} is molecule-based. This is not implemented in
2849 {\gromacs}, which uses atomic virials
2851 The \gromosv{96} force field was parameterized with a Lennard-Jones cut-off
2852 of 1.4 nm, so be sure to use a Lennard-Jones cut-off ({\tt rvdw}) of at least 1.4.
2853 A larger cut-off is possible because the Lennard-Jones potential and forces
2854 are almost zero beyond 1.4 nm.
2856 \subsubsection{\gromosv{96} files\swapindexquiet{GROMOS96}{files}}
2857 {\gromacs} can read and write \gromosv{96} coordinate and trajectory files.
2858 These files should have the extension {\tt .g96}.
2859 Such a file can be a \gromosv{96} initial/final
2860 configuration file, a coordinate trajectory file, or a combination of both.
2861 The file is fixed format; all floats are written as 15.9, and as such, files can get huge.
2862 {\gromacs} supports the following data blocks in the given order:
2872 POSITION/POSITIONRED (mandatory)
2873 VELOCITY/VELOCITYRED (optional)
2878 See the \gromosv{96} manual~\cite{gromos96} for a complete description
2879 of the blocks. {\bf Note} that all {\gromacs} programs can read compressed
2880 (.Z) or gzipped (.gz) files.
2882 \subsection{OPLS/AA\index{OPLS/AA force field}}
2884 \subsection{AMBER\index{AMBER force field}}
2886 As of version 4.5, {\gromacs} provides native support for the following AMBER force fields:
2889 \item AMBER94~\cite{Cornell1995}
2890 \item AMBER96~\cite{Kollman1996}
2891 \item AMBER99~\cite{Wang2000}
2892 \item AMBER99SB~\cite{Hornak2006}
2893 \item AMBER99SB-ILDN~\cite{Lindorff2010}
2894 \item AMBER03~\cite{Duan2003}
2895 \item AMBERGS~\cite{Garcia2002}
2898 \subsection{CHARMM\index{CHARMM force field}}
2900 As of version 4.5, {\gromacs} supports the CHARMM27 force field for proteins~\cite{mackerell04, mackerell98}, lipids~\cite{feller00} and nucleic acids~\cite{foloppe00}. The protein parameters (and to some extent the lipid and nucleic acid parameters) were thoroughly tested -- both by comparing potential energies between the port and the standard parameter set in the CHARMM molecular simulation package, as well by how the protein force field behaves together with {\gromacs}-specific techniques such as virtual sites (enabling long time steps) and a fast implicit solvent recently implemented~\cite{Larsson10} -- and the details and results are presented in the paper by Bjelkmar et al.~\cite{Bjelkmar10}. The nucleic acid parameters, as well as the ones for HEME, were converted and tested by Michel Cuendet.
2902 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.
2904 \subsection{Coarse-grained force-fields}
2905 \index{force-field, coarse-grained}
2906 \label{sec:cg-forcefields}
2907 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.
2909 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:
2911 \item Conserving free energies
2913 \item Simplex method
2914 \item MARTINI force-field (see next section)
2916 \item Conserving distributions (like the radial distribution function), so-called structure-based coarse-graining
2918 \item (iterative) Boltzmann inversion
2919 \item Inverse Monte Carlo
2921 \item Conversing forces
2923 \item Force matching
2927 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}.
2929 \subsection{MARTINI\index{Martini force field}}
2931 The MARTINI force field is a coarse-grain parameter set that allows for the construction
2932 of many systems, including proteins and membranes.
2934 \subsection{PLUM\index{PLUM force field}}
2936 The \normindex{PLUM force field}~\cite{bereau12} is an example of a solvent-free
2937 protein-membrane model for which the membrane was derived from structure-based
2938 coarse-graining~\cite{wang_jpcb10}. A {\gromacs} implementation can be found at
2939 \href{http://code.google.com/p/plumx/}{code.google.com/p/plumx}.
2941 % LocalWords: dihedrals centro ij dV dr LJ lj rcl jj Bertelot OPLS bh bham rf
2942 % LocalWords: coul defunits grompp crf vcrf fcrf Tironi Debye kgrf cgrf krf dx
2943 % LocalWords: PPPM der Waals erfc lr elstat chirality bstretch bondpot kT kJ
2944 % LocalWords: anharmonic morse mol betaij expminx SPC timestep fs FENE ijk kj
2945 % LocalWords: anglepot CHARMm UB ik rr substituents ijkl Ryckaert Bellemans rb
2946 % LocalWords: alkanes pdb gmx IUPAC IUB jkl cis rbdih crb kcal cubicspline xvg
2947 % LocalWords: topfile mdrun posres ar dihr lcllll NMR nmr lcllllll NOEs lclll
2948 % LocalWords: rav preprocessor ccccccccc ai aj fac disre mdp multi topol tpr
2949 % LocalWords: fc ravdisre nstdisreout dipolar lll ccc orientational MSD const
2950 % LocalWords: orire fitgrp nstorireout Drude intra Noskov et al fecalc coulrf
2951 % LocalWords: polarizabilities parameterized sigeps Ek sc softcore GROMOS NBF
2952 % LocalWords: hydrogens alkane vdwtype coulombtype mdpopt rlist rcoulomb rvdw
2953 % LocalWords: nstlist virial funcparm VdW jk jl fvsite fd vsites lcr vsitetop
2954 % LocalWords: vsite lclllll lcl parameterize parameterization enercorr avcsix
2955 % LocalWords: pcorr ecorr totalcoulomb dir fourierspacing ewald rtol Darden gz
2956 % LocalWords: FFT parallelization MPMD mpmd pme fft hoc Gromos gromos oxygens
2957 % LocalWords: virials POSITIONRED VELOCITYRED gzipped Charmm Larsson Bjelkmar
2958 % LocalWords: Cuendet CMAP nocmap dihedral Lennard covalent Verlet
2959 % LocalWords: Berthelot nm flexwat ferguson itp harmonicangle versa
2960 % LocalWords: harmonicbond atomtypes dihedraltypes equilibrated fdn
2961 % LocalWords: distancerestraint LINCS Coulombic ja jb il SPME ILDN
2962 % LocalWords: Hamiltonians atomtype AMBERGS rtp cmap graining VOTCA