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 A plain Coulomb interaction should only be used without cut-off or when all pairs fall within the cut-off, since there is an abrupt, large change in the force at the cut-off. In case you do want to use a cut-off, the potential can be shifted by a constant to make the potential the integral of the force. With the group cut-off scheme, this shift is only applied to non-excluded pairs. With the Verlet cut-off scheme, the shift is also applied to excluded pairs and self interactions, which makes the potential equivalent to a reaction field with $\epsrf=1$ (see below).
183 In {\gromacs} the relative \swapindex{dielectric}{constant}
185 may be set in the in the input for {\tt grompp}.
187 %\ifthenelse{\equal{\gmxlite}{1}}{}{
188 \subsection{Coulomb interaction with \normindex{reaction field}}
190 The Coulomb interaction can be modified for homogeneous systems by
191 assuming a constant dielectric environment beyond the cut-off $r_c$
192 with a dielectric constant of {$\epsrf$}. The interaction then reads:
195 f \frac{q_i q_j}{\epsr\rij}\left[1+\frac{\epsrf-\epsr}{2\epsrf+\epsr}
196 \,\frac{\rij^3}{r_c^3}\right]
197 - f\frac{q_i q_j}{\epsr r_c}\,\frac{3\epsrf}{2\epsrf+\epsr}
200 in which the constant expression on the right makes the potential
201 zero at the cut-off $r_c$. For charged cut-off spheres this corresponds
202 to neutralization with a homogeneous background charge.
203 We can rewrite \eqnref{vcrf} for simplicity as
205 V_{crf} ~=~ f \frac{q_i q_j}{\epsr}\left[\frac{1}{\rij} + k_{rf}~ \rij^2 -c_{rf}\right]
209 k_{rf} &=& \frac{1}{r_c^3}\,\frac{\epsrf-\epsr}{(2\epsrf+\epsr)} \label{eqn:krf}\\
210 c_{rf} &=& \frac{1}{r_c}+k_{rf}\,r_c^2 ~=~ \frac{1}{r_c}\,\frac{3\epsrf}{(2\epsrf+\epsr)}
213 For large $\epsrf$ the $k_{rf}$ goes to $r_c^{-3}/2$,
214 while for $\epsrf$ = $\epsr$ the correction vanishes.
216 the modified interaction is plotted, and it is clear that the derivative
217 with respect to $\rij$ (= -force) goes to zero at the cut-off distance.
218 The force derived from this potential reads:
220 \ve{F}_i(\rvij) = f \frac{q_i q_j}{\epsr}\left[\frac{1}{\rij^2} - 2 k_{rf}\rij\right]\rnorm
223 The reaction-field correction should also be applied to all excluded
224 atoms pairs, including self pairs, in which case the normal Coulomb
225 term in \eqnsref{vcrf}{fcrf} is absent.
227 Tironi {\etal} have introduced a generalized reaction field in which
228 the dielectric continuum beyond the cut-off $r_c$ also has an ionic strength
229 $I$~\cite{Tironi95}. In this case we can rewrite the constants $k_{rf}$ and
230 $c_{rf}$ using the inverse Debye screening length $\kappa$:
233 \frac{2 I \,F^2}{\varepsilon_0 \epsrf RT}
234 = \frac{F^2}{\varepsilon_0 \epsrf RT}\sum_{i=1}^{K} c_i z_i^2 \\
235 k_{rf} &=& \frac{1}{r_c^3}\,
236 \frac{(\epsrf-\epsr)(1 + \kappa r_c) + \half\epsrf(\kappa r_c)^2}
237 {(2\epsrf + \epsr)(1 + \kappa r_c) + \epsrf(\kappa r_c)^2}
239 c_{rf} &=& \frac{1}{r_c}\,
240 \frac{3\epsrf(1 + \kappa r_c + \half(\kappa r_c)^2)}
241 {(2\epsrf+\epsr)(1 + \kappa r_c) + \epsrf(\kappa r_c)^2}
244 where $F$ is Faraday's constant, $R$ is the ideal gas constant, $T$
245 the absolute temperature, $c_i$ the molar concentration for species
246 $i$ and $z_i$ the charge number of species $i$ where we have $K$
247 different species. In the limit of zero ionic strength ($\kappa=0$)
248 \eqnsref{kgrf}{cgrf} reduce to the simple forms of \eqnsref{krf}{crf}
251 \subsection{Modified non-bonded interactions}
252 \label{sec:mod_nb_int}
253 In {\gromacs}, the non-bonded potentials can be
254 modified by a shift function. The purpose of this is to replace the
255 truncated forces by forces that are continuous and have continuous
256 derivatives at the \normindex{cut-off} radius. With such forces the
257 timestep integration produces much smaller errors and there are no
258 such complications as creating charges from dipoles by the truncation
259 procedure. In fact, by using shifted forces there is no need for
260 charge groups in the construction of neighbor lists. However, the
261 shift function produces a considerable modification of the Coulomb
262 potential. Unless the ``missing'' long-range potential is properly
263 calculated and added (through the use of PPPM, Ewald, or PME), the
264 effect of such modifications must be carefully evaluated. The
265 modification of the Lennard-Jones dispersion and repulsion is only
266 minor, but it does remove the noise caused by cut-off effects.
268 There is {\em no} fundamental difference between a switch function
269 (which multiplies the potential with a function) and a shift function
270 (which adds a function to the force or potential)~\cite{Spoel2006a}. The switch
271 function is a special case of the shift function, which we apply to
272 the {\em force function} $F(r)$, related to the electrostatic or
273 van der Waals force acting on particle $i$ by particle $j$ as:
275 \ve{F}_i = c F(r_{ij}) \frac{\rvij}{r_{ij}}
277 For pure Coulomb or Lennard-Jones interactions
278 $F(r)=F_\alpha(r)=r^{-(\alpha+1)}$.
279 The shifted force $F_s(r)$ can generally be written as:
283 F_s(r)~=&~F_\alpha(r) & r < r_1 \\
285 F_s(r)~=&~F_\alpha(r)+S(r) & r_1 \le r < r_c \\
286 F_s(r)~=&~0 & r_c \le r
289 When $r_1=0$ this is a traditional shift function, otherwise it acts as a
290 switch function. The corresponding shifted coulomb potential then reads:
292 V_s(r_{ij}) = f \Phi_s (r_{ij}) q_i q_j
294 where $\Phi(r)$ is the potential function
296 \Phi_s(r) = \int^{\infty}_r~F_s(x)\, dx
299 The {\gromacs} shift function should be smooth at the boundaries, therefore
300 the following boundary conditions are imposed on the shift function:
305 S(r_c) &=&-F_\alpha(r_c) \\
306 S'(r_c) &=&-F_\alpha'(r_c)
309 A 3$^{rd}$ degree polynomial of the form
311 S(r) = A(r-r_1)^2 + B(r-r_1)^3
313 fulfills these requirements. The constants A and B are given by the
314 boundary condition at $r_c$:
318 A &~=~& -\displaystyle
319 \frac{(\alpha+4)r_c~-~(\alpha+1)r_1} {r_c^{\alpha+2}~(r_c-r_1)^2} \\
320 B &~=~& \displaystyle
321 \frac{(\alpha+3)r_c~-~(\alpha+1)r_1}{r_c^{\alpha+2}~(r_c-r_1)^3}
324 Thus the total force function is:
326 F_s(r) = \frac{\alpha}{r^{\alpha+1}} + A(r-r_1)^2 + B(r-r_1)^3
328 and the potential function reads:
330 \Phi(r) = \frac{1}{r^\alpha} - \frac{A}{3} (r-r_1)^3 - \frac{B}{4} (r-r_1)^4 - C
334 C = \frac{1}{r_c^\alpha} - \frac{A}{3} (r_c-r_1)^3 - \frac{B}{4} (r_c-r_1)^4
337 When $r_1$ = 0, the modified Coulomb force function is
339 F_s(r) = \frac{1}{r^2} - \frac{5 r^2}{r_c^4} + \frac{4 r^3}{r_c^5}
341 which is identical to the {\em \index{parabolic force}}
342 function recommended to be used as a short-range function in
343 conjunction with a \swapindex{Poisson}{solver}
344 for the long-range part~\cite{Berendsen93a}.
345 The modified Coulomb potential function is:
347 \Phi(r) = \frac{1}{r} - \frac{5}{3r_c} + \frac{5r^3}{3r_c^4} - \frac{r^4}{r_c^5}
349 See also \figref{shift}.
352 \centerline{\includegraphics[width=10cm]{plots/shiftf}}
353 \caption[The Coulomb Force, Shifted Force and Shift Function
354 $S(r)$,.]{The Coulomb Force, Shifted Force and Shift Function $S(r)$,
355 using r$_1$ = 2 and r$_c$ = 4.}
359 \subsection{Modified short-range interactions with Ewald summation}
360 When Ewald summation\index{Ewald sum} or \seeindex{particle-mesh
361 Ewald}{PME}\index{Ewald, particle-mesh} is used to calculate the
362 long-range interactions, the
363 short-range Coulomb potential must also be modified, similar to the
364 switch function above. In this case the short range potential is given
367 V(r) = f \frac{\mbox{erfc}(\beta r_{ij})}{r_{ij}} q_i q_j,
369 where $\beta$ is a parameter that determines the relative weight
370 between the direct space sum and the reciprocal space sum and erfc$(x)$ is
371 the complementary error function. For further
372 details on long-range electrostatics, see \secref{lr_elstat}.
373 %} % Brace matches ifthenelse test for gmxlite
376 \section{Bonded interactions}
377 Bonded interactions are based on a fixed list of atoms. They are not
378 exclusively pair interactions, but include 3- and 4-body interactions
379 as well. There are {\em bond stretching} (2-body), {\em bond angle}
380 (3-body), and {\em dihedral angle} (4-body) interactions. A special
381 type of dihedral interaction (called {\em improper dihedral}) is used
382 to force atoms to remain in a plane or to prevent transition to a
383 configuration of opposite chirality (a mirror image).
385 \subsection{Bond stretching}
387 \subsubsection{Harmonic potential}
388 \label{subsec:harmonicbond}
389 The \swapindex{bond}{stretching} between two covalently bonded atoms
390 $i$ and $j$ is represented by a harmonic potential:
393 \centerline{\raisebox{2cm}{\includegraphics[width=5cm]{plots/bstretch}}\includegraphics[width=7cm]{plots/f-bond}}
394 \caption[Bond stretching.]{Principle of bond stretching (left), and the bond
395 stretching potential (right).}
396 \label{fig:bstretch1}
400 V_b~(\rij) = \half k^b_{ij}(\rij-b_{ij})^2
402 See also \figref{bstretch1}, with the force given by:
404 \ve{F}_i(\rvij) = k^b_{ij}(\rij-b_{ij}) \rnorm
407 %\ifthenelse{\equal{\gmxlite}{1}}{}{
408 \subsubsection{Fourth power potential}
409 \label{subsec:G96bond}
410 In the \gromosv{96} force field~\cite{gromos96}, the covalent bond potential
411 is, for reasons of computational efficiency, written as:
413 V_b~(\rij) = \frac{1}{4}k^b_{ij}\left(\rij^2-b_{ij}^2\right)^2
415 The corresponding force is:
417 \ve{F}_i(\rvij) = k^b_{ij}(\rij^2-b_{ij}^2)~\rvij
419 The force constants for this form of the potential are related to the usual
420 harmonic force constant $k^{b,\mathrm{harm}}$ (\secref{bondpot}) as
422 2 k^b b_{ij}^2 = k^{b,\mathrm{harm}}
424 The force constants are mostly derived from the harmonic ones used in
425 \gromosv{87}~\cite{biomos}. Although this form is computationally more
427 (because no square root has to be evaluated), it is conceptually more
428 complex. One particular disadvantage is that since the form is not harmonic,
429 the average energy of a single bond is not equal to $\half kT$ as it is for
430 the normal harmonic potential.
432 \subsection{\normindex{Morse potential} bond stretching}
433 \label{subsec:Morsebond}
434 %\author{F.P.X. Everdij}
436 For some systems that require an anharmonic bond stretching potential,
437 the Morse potential~\cite{Morse29}
438 between two atoms {\it i} and {\it j} is available
439 in {\gromacs}. This potential differs from the harmonic potential in
440 that it has an asymmetric potential well and a zero force at infinite
441 distance. The functional form is:
443 \displaystyle V_{morse} (r_{ij}) = D_{ij} [1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))]^2,
445 See also \figref{morse}, and the corresponding force is:
448 \displaystyle {\bf F}_{morse} ({\bf r}_{ij})&=&2 D_{ij} \beta_{ij} r_{ij} \exp(-\beta_{ij}(r_{ij}-b_{ij})) * \\
449 \displaystyle \: & \: &[1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))] \frac{\displaystyle {\bf r}_{ij}}{\displaystyle r_{ij}},
452 where \( \displaystyle D_{ij} \) is the depth of the well in kJ/mol,
453 \( \displaystyle \beta_{ij} \) defines the steepness of the well (in
454 nm\(^{-1} \)), and \( \displaystyle b_{ij} \) is the equilibrium
455 distance in nm. The steepness parameter \( \displaystyle \beta_{ij}
456 \) can be expressed in terms of the reduced mass of the atoms {\it i}
457 and {\it j}, the fundamental vibration frequency \( \displaystyle
458 \omega_{ij} \) and the well depth \( \displaystyle D_{ij} \):
460 \displaystyle \beta_{ij}= \omega_{ij} \sqrt{\frac{\mu_{ij}}{2 D_{ij}}}
462 and because \( \displaystyle \omega = \sqrt{k/\mu} \), one can rewrite \( \displaystyle \beta_{ij} \) in terms of the harmonic force constant \( \displaystyle k_{ij} \):
464 \displaystyle \beta_{ij}= \sqrt{\frac{k_{ij}}{2 D_{ij}}}
467 For small deviations \( \displaystyle (r_{ij}-b_{ij}) \), one can
468 approximate the \( \displaystyle \exp \)-term to first-order using a
471 \displaystyle \exp(-x) \approx 1-x
474 and substituting \eqnref{betaij} and \eqnref{expminx} in the functional form:
477 \displaystyle V_{morse} (r_{ij})&=&D_{ij} [1 - \exp(-\beta_{ij}(r_{ij}-b_{ij}))]^2\\
478 \displaystyle \:&=&D_{ij} [1 - (1 -\sqrt{\frac{k_{ij}}{2 D_{ij}}}(r_{ij}-b_{ij}))]^2\\
479 \displaystyle \:&=&\frac{1}{2} k_{ij} (r_{ij}-b_{ij}))^2
482 we recover the harmonic bond stretching potential.
485 \centerline{\includegraphics[width=7cm]{plots/f-morse}}
486 \caption{The Morse potential well, with bond length 0.15 nm.}
490 \subsection{Cubic bond stretching potential}
491 \label{subsec:cubicbond}
492 Another anharmonic bond stretching potential that is slightly simpler
493 than the Morse potential adds a cubic term in the distance to the
494 simple harmonic form:
496 V_b~(\rij) = k^b_{ij}(\rij-b_{ij})^2 + k^b_{ij}k^{cub}_{ij}(\rij-b_{ij})^3
498 A flexible \normindex{water} model (based on
499 the SPC water model~\cite{Berendsen81}) including
500 a cubic bond stretching potential for the O-H bond
501 was developed by Ferguson~\cite{Ferguson95}. This model was found
502 to yield a reasonable infrared spectrum. The Ferguson water model is
503 available in the {\gromacs} library ({\tt flexwat-ferguson.itp}).
504 It should be noted that the potential is asymmetric: overstretching leads to
505 infinitely low energies. The \swapindex{integration}{timestep} is therefore
508 The force corresponding to this potential is:
510 \ve{F}_i(\rvij) = 2k^b_{ij}(\rij-b_{ij})~\rnorm + 3k^b_{ij}k^{cub}_{ij}(\rij-b_{ij})^2~\rnorm
513 \subsection{FENE bond stretching potential\index{FENE potential}}
514 \label{subsec:FENEbond}
515 In coarse-grained polymer simulations the beads are often connected
516 by a FENE (finitely extensible nonlinear elastic) potential~\cite{Warner72}:
518 V_{\mbox{\small FENE}}(\rij) =
519 -\half k^b_{ij} b^2_{ij} \log\left(1 - \frac{\rij^2}{b^2_{ij}}\right)
521 The potential looks complicated, but the expression for the force is simpler:
523 F_{\mbox{\small FENE}}(\rvij) =
524 -k^b_{ij} \left(1 - \frac{\rij^2}{b^2_{ij}}\right)^{-1} \rvij
526 At short distances the potential asymptotically goes to a harmonic
527 potential with force constant $k^b$, while it diverges at distance $b$.
528 %} % Brace matches ifthenelse test for gmxlite
530 \subsection{Harmonic angle potential}
531 \label{subsec:harmonicangle}
532 \newcommand{\tijk}{\theta_{ijk}}
533 The bond-\swapindex{angle}{vibration} between a triplet of atoms $i$ - $j$ - $k$
534 is also represented by a harmonic potential on the angle $\tijk$
537 \centerline{\raisebox{1cm}{\includegraphics[width=5cm]{plots/angle}}\includegraphics[width=7cm]{plots/f-angle}}
538 \caption[Angle vibration.]{Principle of angle vibration (left) and the
539 bond angle potential (right).}
544 V_a(\tijk) = \half k^{\theta}_{ijk}(\tijk-\tijk^0)^2
546 As the bond-angle vibration is represented by a harmonic potential, the
547 form is the same as the bond stretching (\figref{bstretch1}).
549 The force equations are given by the chain rule:
552 \Fvi ~=~ -\displaystyle\frac{d V_a(\tijk)}{d \rvi} \\
553 \Fvk ~=~ -\displaystyle\frac{d V_a(\tijk)}{d \rvk} \\
556 ~ \mbox{ ~ where ~ } ~
557 \tijk = \arccos \frac{(\rvij \cdot \ve{r}_{kj})}{r_{ij}r_{kj}}
559 The numbering $i,j,k$ is in sequence of covalently bonded atoms. Atom
560 $j$ is in the middle; atoms $i$ and $k$ are at the ends (see \figref{angle}).
561 {\bf Note} that in the input in topology files, angles are given in degrees and
562 force constants in kJ/mol/rad$^2$.
564 %\ifthenelse{\equal{\gmxlite}{1}}{}{
565 \subsection{Cosine based angle potential}
566 \label{subsec:G96angle}
567 In the \gromosv{96} force field a simplified function is used to represent angle
570 V_a(\tijk) = \half k^{\theta}_{ijk}\left(\cos(\tijk) - \cos(\tijk^0)\right)^2
575 \cos(\tijk) = \frac{\rvij\cdot\ve{r}_{kj}}{\rij r_{kj}}
577 The corresponding force can be derived by partial differentiation with respect
578 to the atomic positions. The force constants in this function are related
579 to the force constants in the harmonic form $k^{\theta,\mathrm{harm}}$
580 (\ssecref{harmonicangle}) by:
582 k^{\theta} \sin^2(\tijk^0) = k^{\theta,\mathrm{harm}}
584 In the \gromosv{96} manual there is a much more complicated conversion formula
585 which is temperature dependent. The formulas are equivalent at 0 K
586 and the differences at 300 K are on the order of 0.1 to 0.2\%.
587 {\bf Note} that in the input in topology files, angles are given in degrees and
588 force constants in kJ/mol.
590 \subsection{Restricted bending potential}
592 The restricted bending (ReB) potential~\cite{MonicaGoga2013} prevents the bending angle $\theta$
593 from reaching the $180^{\circ}$ value. In this way, the numerical instabilities
594 due to the calculation of the torsion angle and potential are eliminated when
595 performing coarse-grained molecular dynamics simulations.
597 To systematically hinder the bending angles from reaching the $180^{\circ}$ value,
598 the bending potential \ref{eq:G96angle} is divided by a $\sin^2\theta$ factor:
601 V_{\rm ReB}(\theta_i) = \frac{1}{2} k_{\theta} \frac{(\cos\theta_i - \cos\theta_0)^2}{\sin^2\theta_i}.
605 Figure ~\figref{ReB} shows the comparison between the ReB potential, \ref{eq:ReB},
606 and the standard one \ref{eq:G96angle}.
609 \centerline{\includegraphics[width=10cm]{plots/fig-02}}
611 \caption{Bending angle potentials: cosine harmonic (solid black line), angle harmonic
612 (dashed black line) and restricted bending (red) with the same bending constant
613 $k_{\theta}=85$ kJ mol$^{-1}$ and equilibrium angle $\theta_0=130^{\circ}$.
614 The orange line represents the sum of a cosine harmonic ($k =50$ kJ mol$^{-1}$)
615 with a restricted bending ($k =25$ kJ mol$^{-1}$) potential, both with $\theta_0=130^{\circ}$.}
619 The wall of the ReB potential is very repulsive in the region close to $180^{\circ}$ and,
620 as a result, the bending angles are kept within a safe interval, far from instabilities.
621 The power $2$ of $\sin\theta_i$ in the denominator has been chosen to guarantee this behavior
622 and allows an elegant differentiation:
625 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}}.
629 Due to its construction, the restricted bending potential cannot be used for equilibrium
630 $\theta_0$ values too close to $0^{\circ}$ or $180^{\circ}$ (from experience, at least $10^{\circ}$
631 difference is recommended). It is very important that, in the starting configuration,
632 all the bending angles have to be in the safe interval to avoid initial instabilities.
633 This bending potential can be used in combination with any form of torsion potential.
634 It will always prevent three consecutive particles from becoming collinear and,
635 as a result, any torsion potential will remain free of singularities.
636 It can be also added to a standard bending potential to affect the angle around $180^{\circ}$,
637 but to keep its original form around the minimum (see the orange curve in \figref{ReB}).
640 \subsection{Urey-Bradley potential}
641 \label{subsec:Urey-Bradley}
642 The \swapindex{Urey-Bradley bond-angle}{vibration} between a triplet
643 of atoms $i$ - $j$ - $k$ is represented by a harmonic potential on the
644 angle $\tijk$ and a harmonic correction term on the distance between
645 the atoms $i$ and $k$. Although this can be easily written as a simple
646 sum of two terms, it is convenient to have it as a single entry in the
647 topology file and in the output as a separate energy term. It is used mainly
648 in the CHARMm force field~\cite{BBrooks83}. The energy is given by:
651 V_a(\tijk) = \half k^{\theta}_{ijk}(\tijk-\tijk^0)^2 + \half k^{UB}_{ijk}(r_{ik}-r_{ik}^0)^2
654 The force equations can be deduced from sections~\ssecref{harmonicbond}
655 and~\ssecref{harmonicangle}.
657 \subsection{Bond-Bond cross term}
658 \label{subsec:bondbondcross}
659 The bond-bond cross term for three particles $i, j, k$ forming bonds
660 $i-j$ and $k-j$ is given by~\cite{Lawrence2003b}:
662 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)
665 where $k_{rr'}$ is the force constant, and $r_{1e}$ and $r_{2e}$ are the
666 equilibrium bond lengths of the $i-j$ and $k-j$ bonds respectively. The force
667 associated with this potential on particle $i$ is:
669 \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|}
671 The force on atom $k$ can be obtained by swapping $i$ and $k$ in the above
672 equation. Finally, the force on atom $j$ follows from the fact that the sum
673 of internal forces should be zero: $\ve{F}_j = -\ve{F}_i-\ve{F}_k$.
675 \subsection{Bond-Angle cross term}
676 \label{subsec:bondanglecross}
677 The bond-angle cross term for three particles $i, j, k$ forming bonds
678 $i-j$ and $k-j$ is given by~\cite{Lawrence2003b}:
680 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)
682 where $k_{r\theta}$ is the force constant, $r_{3e}$ is the $i-k$ distance,
683 and the other constants are the same as in Equation~\ref{crossbb}. The force
684 associated with the potential on atom $i$ is:
686 \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|} \\
687 + \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]
690 \subsection{Quartic angle potential}
691 \label{subsec:quarticangle}
692 For special purposes there is an angle potential
693 that uses a fourth order polynomial:
695 V_q(\tijk) ~=~ \sum_{n=0}^5 C_n (\tijk-\tijk^0)^n
697 %} % Brace matches ifthenelse test for gmxlite
699 %% new commands %%%%%%%%%%%%%%%%%%%%%%
700 \newcommand{\rvkj}{{\bf r}_{kj}}
701 \newcommand{\rkj}{r_{kj}}
702 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
704 \subsection{Improper dihedrals\swapindexquiet{improper}{dihedral}}
706 Improper dihedrals are meant to keep \swapindex{planar}{group}s ({\eg}
707 aromatic rings) planar, or to prevent molecules from flipping over to their
708 \normindex{mirror image}s, see \figref{imp}.
711 \centerline{\includegraphics[width=4cm]{plots/ring-imp}\hspace{1cm}
712 \includegraphics[width=3cm]{plots/subst-im}\hspace{1cm}\includegraphics[width=3cm]{plots/tetra-im}}
713 \caption[Improper dihedral angles.]{Principle of improper
714 dihedral angles. Out of plane bending for rings (left), substituents
715 of rings (middle), out of tetrahedral (right). The improper dihedral
716 angle $\xi$ is defined as the angle between planes (i,j,k) and (j,k,l)
721 \subsubsection{Improper dihedrals: harmonic type}
722 \label{subsec:harmonicimproperdihedral}
723 The simplest improper dihedral potential is a harmonic potential; it is plotted in
726 V_{id}(\xi_{ijkl}) = \half k_{\xi}(\xi_{ijkl}-\xi_0)^2
728 Since the potential is harmonic it is discontinuous,
729 but since the discontinuity is chosen at 180$^\circ$ distance from $\xi_0$
730 this will never cause problems.
731 {\bf Note} that in the input in topology files, angles are given in degrees and
732 force constants in kJ/mol/rad$^2$.
735 \centerline{\includegraphics[width=8cm]{plots/f-imps}}
736 \caption{Improper dihedral potential.}
740 \subsubsection{Improper dihedrals: periodic type}
741 \label{subsec:periodicimproperdihedral}
742 This potential is identical to the periodic proper dihedral (see below).
743 There is a separate dihedral type for this (type 4) only to be able
744 to distinguish improper from proper dihedrals in the parameter section
747 \subsection{Proper dihedrals\swapindexquiet{proper}{dihedral}}
748 For the normal \normindex{dihedral} interaction there is a choice of
749 either the {\gromos} periodic function or a function based on
750 expansion in powers of $\cos \phi$ (the so-called Ryckaert-Bellemans
751 potential). This choice has consequences for the inclusion of special
752 interactions between the first and the fourth atom of the dihedral
753 quadruple. With the periodic {\gromos} potential a special 1-4
754 LJ-interaction must be included; with the Ryckaert-Bellemans potential
755 {\em for alkanes} the \normindex{1-4 interaction}s must be excluded
756 from the non-bonded list. {\bf Note:} Ryckaert-Bellemans potentials
757 are also used in {\eg} the OPLS force field in combination with 1-4
758 interactions. You should therefore not modify topologies generated by
759 {\tt \normindex{pdb2gmx}} in this case.
761 \subsubsection{Proper dihedrals: periodic type}
762 \label{subsec:properdihedral}
763 Proper dihedral angles are defined according to the IUPAC/IUB
764 convention, where $\phi$ is the angle between the $ijk$ and the $jkl$
765 planes, with {\bf zero} corresponding to the {\em cis} configuration
766 ($i$ and $l$ on the same side). There are two dihedral function types
767 in {\gromacs} topology files. There is the standard type 1 which behaves
768 like any other bonded interactions. For certain force fields, type 9
769 is useful. Type 9 allows multiple potential functions to be applied
770 automatically to a single dihedral in the {\tt [ dihedral ]} section
771 when multiple parameters are defined for the same atomtypes
772 in the {\tt [ dihedraltypes ]} section.
775 \centerline{\raisebox{1cm}{\includegraphics[width=5cm]{plots/dih}}\includegraphics[width=7cm]{plots/f-dih}}
776 \caption[Proper dihedral angle.]{Principle of proper dihedral angle
777 (left, in {\em trans} form) and the dihedral angle potential (right).}
781 V_d(\phi_{ijkl}) = k_{\phi}(1 + \cos(n \phi - \phi_s))
784 %\ifthenelse{\equal{\gmxlite}{1}}{}{
785 \subsubsection{Proper dihedrals: Ryckaert-Bellemans function}
786 \label{subsec:RBdihedral}
787 For alkanes, the following proper dihedral potential is often used
788 (see \figref{rbdih}):
790 V_{rb}(\phi_{ijkl}) = \sum_{n=0}^5 C_n( \cos(\psi ))^n,
792 where $\psi = \phi - 180^\circ$. \\
793 {\bf Note:} A conversion from one convention to another can be achieved by
794 multiplying every coefficient \( \displaystyle C_n \)
795 by \( \displaystyle (-1)^n \).
797 An example of constants for $C$ is given in \tabref{crb}.
801 \begin{tabular}{|lr|lr|lr|}
803 $C_0$ & 9.28 & $C_2$ & -13.12 & $C_4$ & 26.24 \\
804 $C_1$ & 12.16 & $C_3$ & -3.06 & $C_5$ & -31.5 \\
808 \caption{Constants for Ryckaert-Bellemans potential (kJ mol$^{-1}$).}
813 \centerline{\includegraphics[width=8cm]{plots/f-rbs}}
814 \caption{Ryckaert-Bellemans dihedral potential.}
818 ({\bf Note:} The use of this potential implies exclusion of LJ interactions
819 between the first and the last atom of the dihedral, and $\psi$ is defined
820 according to the ``polymer convention'' ($\psi_{trans}=0$).)
822 The RB dihedral function can also be used to include Fourier dihedrals
825 V_{rb} (\phi_{ijkl}) ~=~ \frac{1}{2} \left[F_1(1+\cos(\phi)) + F_2(
826 1-\cos(2\phi)) + F_3(1+\cos(3\phi)) + F_4(1-\cos(4\phi))\right]
828 Because of the equalities \( \cos(2\phi) = 2\cos^2(\phi) - 1 \),
829 \( \cos(3\phi) = 4\cos^3(\phi) - 3\cos(\phi) \) and
830 \( \cos(4\phi) = 8\cos^4(\phi) - 8\cos^2(\phi) + 1 \)
831 one can translate the OPLS parameters to
832 Ryckaert-Bellemans parameters as follows:
836 \displaystyle C_0&=&F_2 + \frac{1}{2} (F_1 + F_3)\\
837 \displaystyle C_1&=&\frac{1}{2} (- F_1 + 3 \, F_3)\\
838 \displaystyle C_2&=& -F_2 + 4 \, F_4\\
839 \displaystyle C_3&=&-2 \, F_3\\
840 \displaystyle C_4&=&-4 \, F_4\\
841 \displaystyle C_5&=&0
844 with OPLS parameters in protein convention and RB parameters in
845 polymer convention (this yields a minus sign for the odd powers of
847 \noindent{\bf Note:} Mind the conversion from {\bf kcal mol$^{-1}$} for
848 literature OPLS and RB parameters to {\bf kJ mol$^{-1}$} in {\gromacs}.\\
849 %} % Brace matches ifthenelse test for gmxlite
851 \subsubsection{Proper dihedrals: Fourier function}
852 \label{subsec:Fourierdihedral}
853 The OPLS potential function is given as the first three
854 or four~\cite{Jorgensen2005a} cosine terms of a Fourier series.
855 In {\gromacs} the four term function is implemented:
857 V_{F} (\phi_{ijkl}) ~=~ \frac{1}{2} \left[C_1(1+\cos(\phi)) + C_2(
858 1-\cos(2\phi)) + C_3(1+\cos(3\phi)) + C_4(1+\cos(4\phi))\right],
860 %\ifthenelse{\equal{\gmxlite}{1}}{}{
861 Internally, {\gromacs}
862 uses the Ryckaert-Bellemans code
863 to compute Fourier dihedrals (see above), because this is more efficient.\\
864 \noindent{\bf Note:} Mind the conversion from {\emph kcal mol$^{-1}$} for
865 literature OPLS parameters to {\bf kJ mol$^{-1}$} in {\gromacs}.\\
867 \subsubsection{Proper dihedrals: Restricted torsion potential}
869 In a manner very similar to the restricted bending potential (see \ref{subsec:ReB}),
870 a restricted torsion/dihedral potential is introduced:
873 V_{\rm ReT}(\phi_i) = \frac{1}{2} k_{\phi} \frac{(\cos\phi_i - \cos\phi_0)^2}{\sin^2\phi_i}
877 with the advantages of being a function of $\cos\phi$ (no problems taking the derivative of $\sin\phi$)
878 and of keeping the torsion angle at only one minimum value. In this case, the factor $\sin^2\phi$ does
879 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.
880 For this reason, all the dihedral angles of the starting configuration should have their values in the
881 desired angles interval and the the equilibrium $\phi_0$ value should not be too close to the interval limits
882 (as for the restricted bending potential, described in \ref{subsec:ReB}, at least $10^{\circ}$ difference is recommended).
884 \subsubsection{Proper dihedrals: Combined bending-torsion potential}
886 When the four particles forming the dihedral angle become collinear (this situation will never happen in
887 atomistic simulations, but it can occur in coarse-grained simulations) the calculation of the
888 torsion angle and potential leads to numerical instabilities.
889 One way to avoid this is to use the restricted bending potential (see \ref{subsec:ReB})
890 that prevents the dihedral
891 from reaching the $180^{\circ}$ value.
893 Another way is to disregard any effects of the dihedral becoming ill-defined,
894 keeping the dihedral force and potential calculation continuous in entire angle range
895 by coupling the torsion potential (in a cosine form) with the bending potentials of the
896 adjacent bending angles in a unique expression:
899 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}.
903 This combined bending-torsion (CBT) potential has been proposed by~\cite{BulacuGiessen2005}
904 for polymer melt simulations and is extensively described in~\cite{MonicaGoga2013}.
906 This potential has two main advantages:
909 it does not only depend on the dihedral angle $\phi_i$ (between the $i-2$, $i-1$, $i$ and $i+1$ beads)
910 but also on the bending angles $\theta_{i-1}$ and $\theta_i$ defined from three adjacent beads
911 ($i-2$, $i-1$ and $i$, and $i-1$, $i$ and $i+1$, respectively).
912 The two $\sin^3\theta$ pre-factors, tentatively suggested by~\cite{ScottScheragator1966} and theoretically
913 discussed by~\cite{PaulingBond}, cancel the torsion potential and force when either of the two bending angles
914 approaches the value of $180^\circ$.
916 its dependence on $\phi_i$ is expressed through a polynomial in $\cos\phi_i$ that avoids the singularities in
917 $\phi=0^\circ$ or $180^\circ$ in calculating the torsional force.
920 These two properties make the CBT potential well-behaved for MD simulations with weak constraints
921 on the bending angles or even for steered / non-equilibrium MD in which the bending and torsion angles suffer major
923 When using the CBT potential, the bending potentials for the adjacent $\theta_{i-1}$ and $\theta_i$ may have any form.
924 It is also possible to leave out the two angle bending terms ($\theta_{i-1}$ and $\theta_{i}$) completely.
925 \figref{CBT} illustrates the difference between a torsion potential with and without the $\sin^{3}\theta$ factors
926 (blue and gray curves, respectively).
929 \centerline{\includegraphics[width=10cm]{plots/fig-04}}
930 \caption{Blue: surface plot of the combined bending-torsion potential
931 (\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$)
932 when, for simplicity, the bending angles behave the same ($\theta_1=\theta_2=\theta$).
933 Gray: the same torsion potential without the $\sin^{3}\theta$ terms (Ryckaert-Bellemans type).
934 $\phi$ is the dihedral angle.}
938 Additionally, the derivative of $V_{CBT}$ with respect to the Cartesian variables is straightforward:
941 \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}} +
942 \frac{\partial V_{\rm CBT}}{\partial \theta_{i }} \frac{\partial \theta_{i }}{\partial \vec r_{l}} +
943 \frac{\partial V_{\rm CBT}}{\partial \phi_{i }} \frac{\partial \phi_{i }}{\partial \vec r_{l}}
947 The CBT is based on a cosine form without multiplicity, so it can only be symmetrical around $0^{\circ}$.
948 To obtain an asymmetrical dihedral angle distribution (e.g. only one maximum in [$-180^{\circ}$:$180^{\circ}$] interval),
949 a standard torsion potential such as harmonic angle or periodic cosine potentials should be used instead of a CBT potential.
950 However, these two forms have the inconveniences of the force derivation ($1/\sin\phi$) and of the alignment of beads
951 ($\theta_i$ or $\theta_{i-1} = 0^{\circ}, 180^{\circ}$).
952 Coupling such non-$\cos\phi$ potentials with $\sin^3\theta$ factors does not improve simulation stability since there are
953 cases in which $\theta$ and $\phi$ are simultaneously $180^{\circ}$. The integration at this step would be possible
954 (due to the cancelling of the torsion potential) but the next step would be singular
955 ($\theta$ is not $180^{\circ}$ and $\phi$ is very close to $180^{\circ}$).
957 %\ifthenelse{\equal{\gmxlite}{1}}{}{
958 \subsection{Tabulated bonded interaction functions\index{tabulated bonded interaction function}}
959 \label{subsec:tabulatedinteraction}
960 For full flexibility, any functional shape can be used for
961 bonds, angles and dihedrals through user-supplied tabulated functions.
962 The functional shapes are:
964 V_b(r_{ij}) &=& k \, f^b_n(r_{ij}) \\
965 V_a(\tijk) &=& k \, f^a_n(\tijk) \\
966 V_d(\phi_{ijkl}) &=& k \, f^d_n(\phi_{ijkl})
968 where $k$ is a force constant in units of energy
969 and $f$ is a cubic spline function; for details see \ssecref{cubicspline}.
970 For each interaction, the force constant $k$ and the table number $n$
971 are specified in the topology.
972 There are two different types of bonds, one that generates exclusions (type 8)
973 and one that does not (type 9).
974 For details see \tabref{topfile2}.
975 The table files are supplied to the {\tt mdrun} program.
976 After the table file name an underscore, the letter ``b'' for bonds,
977 ``a'' for angles or ``d'' for dihedrals and the table number are appended.
978 For example, for a bond with $n=0$ (and using the default table file name)
979 the table is read from the file {\tt table_b0.xvg}. Multiple tables can be
980 supplied simply by using different values of $n$, and are applied to the appropriate
981 bonds, as specified in the topology (\tabref{topfile2}).
982 The format for the table files is three columns with $x$, $f(x)$, $-f'(x)$,
983 where $x$ should be uniformly-spaced. Requirements for entries in the topology
984 are given in~\tabref{topfile2}.
985 The setup of the tables is as follows:
987 $x$ is the distance in nm. For distances beyond the table length,
988 {\tt mdrun} will quit with an error message.
990 $x$ is the angle in degrees. The table should go from
991 0 up to and including 180 degrees; the derivative is taken in degrees.
993 $x$ is the dihedral angle in degrees. The table should go from
994 -180 up to and including 180 degrees;
995 the IUPAC/IUB convention is used, {\ie} zero is cis,
996 the derivative is taken in degrees.
997 %} % Brace matches ifthenelse test for gmxlite
1000 Special potentials are used for imposing restraints on the motion of
1001 the system, either to avoid disastrous deviations, or to include
1002 knowledge from experimental data. In either case they are not really
1003 part of the force field and the reliability of the parameters is not
1004 important. The potential forms, as implemented in {\gromacs}, are
1005 mentioned just for the sake of completeness. Restraints and constraints
1006 refer to quite different algorithms in {\gromacs}.
1008 \subsection{Position restraints\swapindexquiet{position}{restraint}}
1009 \label{subsec:positionrestraint}
1010 These are used to restrain particles to fixed reference positions
1011 $\ve{R}_i$. They can be used during equilibration in order to avoid
1012 drastic rearrangements of critical parts ({\eg} to restrain motion
1013 in a protein that is subjected to large solvent forces when the
1014 solvent is not yet equilibrated). Another application is the
1015 restraining of particles in a shell around a region that is simulated
1016 in detail, while the shell is only approximated because it lacks
1017 proper interaction from missing particles outside the
1018 shell. Restraining will then maintain the integrity of the inner
1019 part. For spherical shells, it is a wise procedure to make the force
1020 constant depend on the radius, increasing from zero at the inner
1021 boundary to a large value at the outer boundary. This feature has
1022 not, however, been implemented in {\gromacs}.
1023 \newcommand{\unitv}[1]{\hat{\bf #1}}
1024 \newcommand{\halfje}[1]{\frac{#1}{2}}
1026 The following form is used:
1028 V_{pr}(\ve{r}_i) = \halfje{1}k_{pr}|\rvi-\ve{R}_i|^2
1030 The potential is plotted in \figref{positionrestraint}.
1033 \centerline{\includegraphics[width=8cm]{plots/f-pr}}
1034 \caption{Position restraint potential.}
1035 \label{fig:positionrestraint}
1038 The potential form can be rewritten without loss of generality as:
1040 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]
1046 F_i^x &=& -k_{pr}^x~(x_i - X_i) \\
1047 F_i^y &=& -k_{pr}^y~(y_i - Y_i) \\
1048 F_i^z &=& -k_{pr}^z~(z_i - Z_i)
1051 Using three different force constants the position
1052 restraints can be turned on or off
1053 in each spatial dimension; this means that atoms can be harmonically
1054 restrained to a plane or a line.
1055 Position restraints are applied to a special fixed list of atoms. Such a
1056 list is usually generated by the {\tt \normindex{pdb2gmx}} program.
1058 \subsection{\swapindex{Flat-bottomed}{position restraint}s}
1059 \label{subsec:fbpositionrestraint}
1060 Flat-bottomed position restraints can be used to restrain particles to
1061 part of the simulation volume. No force acts on the restrained
1062 particle within the flat-bottomed region of the potential, however a
1063 harmonic force acts to move the particle to the flat-bottomed region
1064 if it is outside it. It is possible to apply normal and
1065 flat-bottomed position restraints on the same particle (however, only
1066 with the same reference position $\ve{R}_i$). The following general potential
1067 is used (Figure~\ref{fig:fbposres}A):
1069 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}],
1071 where $\ve{R}_i$ is the reference position, $r_\mathrm{fb}$ is the distance
1072 from the center with a flat potential, $k_\mathrm{fb}$ the force constant, and $H$ is the Heaviside step
1073 function. The distance $d_g(\ve{r}_i;\ve{R}_i)$ from the reference
1074 position depends on the geometry $g$ of the flat-bottomed potential.
1077 \centerline{\includegraphics[width=10cm]{plots/fbposres}}
1078 \caption{Flat-bottomed position restraint potential. (A) Not
1079 inverted, (B) inverted.}
1080 \label{fig:fbposres}
1083 The following geometries for the flat-bottomed potential are supported:\newline
1084 {\bfseries Sphere} ($g =1$): The particle is kept in a sphere of given
1085 radius. The force acts towards the center of the sphere. The following distance calculation is used:
1087 d_g(\ve{r}_i;\ve{R}_i) = |\ve{r}_i-\ve{R}_i|
1089 {\bfseries Cylinder} ($g=2$): The particle is kept in a cylinder of given radius
1090 parallel to the $z$-axis. The force from the flat-bottomed potential acts
1091 towards the axis of the cylinder. The $z$-component of the force is zero.
1093 d_g(\ve{r}_i;\ve{R}_i) = \sqrt{ (x_i-X_i)^2 + (y_i - Y_i)^2 }
1095 {\bfseries Layer} ($g=3,4,5$): The particle is kept in a layer defined by the
1096 thickness and the normal of the layer. The layer normal can be parallel to the $x$, $y$, or
1097 $z$-axis. The force acts parallel to the layer normal.\\
1099 d_g(\ve{r}_i;\ve{R}_i) = |x_i-X_i|, \;\;\;\mbox{or}\;\;\;
1100 d_g(\ve{r}_i;\ve{R}_i) = |y_i-Y_i|, \;\;\;\mbox{or}\;\;\;
1101 d_g(\ve{r}_i;\ve{R}_i) = |z_i-Z_i|.
1104 It is possible to apply multiple independent flat-bottomed position
1105 restraints of different geometry on one particle. For example, applying
1106 a cylinder and a layer in $z$ keeps a particle within a
1107 disk. Applying three layers in $x$, $y$, and $z$ keeps the particle within a cuboid.
1109 In addition, it is possible to invert the restrained region with the
1110 unrestrained region, leading to a potential that acts to keep the particle {\it outside} of the volume
1111 defined by $\ve{R}_i$, $g$, and $r_\mathrm{fb}$. That feature is
1112 switched on by defining a negative $r_\mathrm{fb}$ in the
1113 topology. The following potential is used (Figure~\ref{fig:fbposres}B):
1115 V_\mathrm{fb}^{\mathrm{inv}}(\ve{r}_i) = \frac{1}{2}k_\mathrm{fb}
1116 [d_g(\ve{r}_i;\ve{R}_i) - |r_\mathrm{fb}|]^2\,
1117 H[ -(d_g(\ve{r}_i;\ve{R}_i) - |r_\mathrm{fb}|)].
1122 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1123 \subsection{Angle restraints\swapindexquiet{angle}{restraint}}
1124 \label{subsec:anglerestraint}
1125 These are used to restrain the angle between two pairs of particles
1126 or between one pair of particles and the $z$-axis.
1127 The functional form is similar to that of a proper dihedral.
1128 For two pairs of atoms:
1130 V_{ar}(\ve{r}_i,\ve{r}_j,\ve{r}_k,\ve{r}_l)
1131 = k_{ar}(1 - \cos(n (\theta - \theta_0))
1134 \theta = \arccos\left(\frac{\ve{r}_j -\ve{r}_i}{\|\ve{r}_j -\ve{r}_i\|}
1135 \cdot \frac{\ve{r}_l -\ve{r}_k}{\|\ve{r}_l -\ve{r}_k\|} \right)
1137 For one pair of atoms and the $z$-axis:
1139 V_{ar}(\ve{r}_i,\ve{r}_j) = k_{ar}(1 - \cos(n (\theta - \theta_0))
1142 \theta = \arccos\left(\frac{\ve{r}_j -\ve{r}_i}{\|\ve{r}_j -\ve{r}_i\|}
1143 \cdot \left( \begin{array}{c} 0 \\ 0 \\ 1 \\ \end{array} \right) \right)
1145 A multiplicity ($n$) of 2 is useful when you do not want to distinguish
1146 between parallel and anti-parallel vectors.
1147 The equilibrium angle $\theta$ should be between 0 and 180 degrees
1148 for multiplicity 1 and between 0 and 90 degrees for multiplicity 2.
1151 \subsection{Dihedral restraints\swapindexquiet{dihedral}{restraint}}
1152 \label{subsec:dihedralrestraint}
1153 These are used to restrain the dihedral angle $\phi$ defined by four particles
1154 as in an improper dihedral (sec.~\ref{sec:imp}) but with a slightly
1155 modified potential. Using:
1157 \phi' = \left(\phi-\phi_0\right) ~{\rm MOD}~ 2\pi
1160 where $\phi_0$ is the reference angle, the potential is defined as:
1162 V_{dihr}(\phi') ~=~ \left\{
1163 \begin{array}{lcllll}
1164 \half k_{dihr}(\phi'-\phi_0-\Delta\phi)^2
1165 &\mbox{for}& \phi' & > & \Delta\phi \\[1.5ex]
1166 0 &\mbox{for}& \phi' & \le & \Delta\phi \\[1.5ex]
1170 where $\Delta\phi$ is a user defined angle and $k_{dihr}$ is the force
1172 {\bf Note} that in the input in topology files, angles are given in degrees and
1173 force constants in kJ/mol/rad$^2$.
1174 %} % Brace matches ifthenelse test for gmxlite
1176 \subsection{Distance restraints\swapindexquiet{distance}{restraint}}
1177 \label{subsec:distancerestraint}
1179 add a penalty to the potential when the distance between specified
1180 pairs of atoms exceeds a threshold value. They are normally used to
1181 impose experimental restraints from, for instance, experiments in nuclear
1182 magnetic resonance (NMR), on the motion of the system. Thus, MD can be
1183 used for structure refinement using NMR data\index{nmr
1184 refinement}\index{refinement,nmr}.
1185 In {\gromacs} there are three ways to impose restraints on pairs of atoms:
1187 \item Simple harmonic restraints: use {\tt [ bonds ]} type 6
1188 %\ifthenelse{\equal{\gmxlite}{1}}
1190 {(see \secref{excl}).}
1191 \item\label{subsec:harmonicrestraint}Piecewise linear/harmonic restraints: {\tt [ bonds ]} type 10.
1192 \item Complex NMR distance restraints, optionally with pair, time and/or
1195 The last two options will be detailed now.
1197 The potential form for distance restraints is quadratic below a specified
1198 lower bound and between two specified upper bounds, and linear beyond the
1199 largest bound (see \figref{dist}).
1201 V_{dr}(r_{ij}) ~=~ \left\{
1202 \begin{array}{lcllllll}
1203 \half k_{dr}(r_{ij}-r_0)^2
1204 &\mbox{for}& & & r_{ij} & < & r_0 \\[1.5ex]
1205 0 &\mbox{for}& r_0 & \le & r_{ij} & < & r_1 \\[1.5ex]
1206 \half k_{dr}(r_{ij}-r_1)^2
1207 &\mbox{for}& r_1 & \le & r_{ij} & < & r_2 \\[1.5ex]
1208 \half k_{dr}(r_2-r_1)(2r_{ij}-r_2-r_1)
1209 &\mbox{for}& r_2 & \le & r_{ij} & &
1215 \centerline{\includegraphics[width=8cm]{plots/f-dr}}
1216 \caption{Distance Restraint potential.}
1223 \begin{array}{lcllllll}
1224 -k_{dr}(r_{ij}-r_0)\frac{\rvij}{r_{ij}}
1225 &\mbox{for}& & & r_{ij} & < & r_0 \\[1.5ex]
1226 0 &\mbox{for}& r_0 & \le & r_{ij} & < & r_1 \\[1.5ex]
1227 -k_{dr}(r_{ij}-r_1)\frac{\rvij}{r_{ij}}
1228 &\mbox{for}& r_1 & \le & r_{ij} & < & r_2 \\[1.5ex]
1229 -k_{dr}(r_2-r_1)\frac{\rvij}{r_{ij}}
1230 &\mbox{for}& r_2 & \le & r_{ij} & &
1234 For restraints not derived from NMR data, this functionality
1235 will usually suffice and a section of {\tt [ bonds ]} type 10
1236 can be used to apply individual restraints between pairs of
1237 %\ifthenelse{\equal{\gmxlite}{1}}{atoms.}{
1238 atoms, see \ssecref{topfile}.
1239 %} % Brace matches ifthenelse test for gmxlite
1240 For applying restraints derived from NMR measurements, more complex
1241 functionality might be required, which is provided through
1242 the {\tt [~distance_restraints~]} section and is described below.
1244 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1245 \subsubsection{Time averaging\swapindexquiet{time-averaged}{distance restraint}}
1246 Distance restraints based on instantaneous distances can potentially reduce
1247 the fluctuations in a molecule significantly. This problem can be overcome by restraining
1248 to a {\em time averaged} distance~\cite{Torda89}.
1249 The forces with time averaging are:
1252 \begin{array}{lcllllll}
1253 -k^a_{dr}(\bar{r}_{ij}-r_0)\frac{\rvij}{r_{ij}}
1254 &\mbox{for}& & & \bar{r}_{ij} & < & r_0 \\[1.5ex]
1255 0 &\mbox{for}& r_0 & \le & \bar{r}_{ij} & < & r_1 \\[1.5ex]
1256 -k^a_{dr}(\bar{r}_{ij}-r_1)\frac{\rvij}{r_{ij}}
1257 &\mbox{for}& r_1 & \le & \bar{r}_{ij} & < & r_2 \\[1.5ex]
1258 -k^a_{dr}(r_2-r_1)\frac{\rvij}{r_{ij}}
1259 &\mbox{for}& r_2 & \le & \bar{r}_{ij} & &
1262 where $\bar{r}_{ij}$ is given by an exponential running average with decay time $\tau$:
1264 \bar{r}_{ij} ~=~ < r_{ij}^{-3} >^{-1/3}
1267 The force constant $k^a_{dr}$ is switched on slowly to compensate for
1268 the lack of history at the beginning of the simulation:
1270 k^a_{dr} = k_{dr} \left(1-\exp\left(-\frac{t}{\tau}\right)\right)
1272 Because of the time averaging, we can no longer speak of a distance restraint
1275 This way an atom can satisfy two incompatible distance restraints
1276 {\em on average} by moving between two positions.
1277 An example would be an amino acid side-chain that is rotating around
1278 its $\chi$ dihedral angle, thereby coming close to various other groups.
1279 Such a mobile side chain can give rise to multiple NOEs that can not be
1280 fulfilled by a single structure.
1282 The computation of the time
1283 averaged distance in the {\tt mdrun} program is done in the following fashion:
1286 \overline{r^{-3}}_{ij}(0) &=& r_{ij}(0)^{-3} \\
1287 \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]
1288 \label{eqn:ravdisre}
1292 When a pair is within the bounds, it can still feel a force
1293 because the time averaged distance can still be beyond a bound.
1294 To prevent the protons from being pulled too close together, a mixed
1295 approach can be used. In this approach, the penalty is zero when the
1296 instantaneous distance is within the bounds, otherwise the violation is
1297 the square root of the product of the instantaneous violation and the
1298 time averaged violation:
1301 \begin{array}{lclll}
1302 k^a_{dr}\sqrt{(r_{ij}-r_0)(\bar{r}_{ij}-r_0)}\frac{\rvij}{r_{ij}}
1303 & \mbox{for} & r_{ij} < r_0 & \mbox{and} & \bar{r}_{ij} < r_0 \\[1.5ex]
1305 \mbox{min}\left(\sqrt{(r_{ij}-r_1)(\bar{r}_{ij}-r_1)},r_2-r_1\right)
1306 \frac{\rvij}{r_{ij}}
1307 & \mbox{for} & r_{ij} > r_1 & \mbox{and} & \bar{r}_{ij} > r_1 \\[1.5ex]
1312 \subsubsection{Averaging over multiple pairs\swapindexquiet{ensemble-averaged}{distance restraint}}
1314 Sometimes it is unclear from experimental data which atom pair
1315 gives rise to a single NOE, in other occasions it can be obvious that
1316 more than one pair contributes due to the symmetry of the system, {\eg} a
1317 methyl group with three protons. For such a group, it is not possible
1318 to distinguish between the protons, therefore they should all be taken into
1319 account when calculating the distance between this methyl group and another
1320 proton (or group of protons).
1321 Due to the physical nature of magnetic resonance, the intensity of the
1322 NOE signal is inversely proportional to the sixth power of the inter-atomic
1324 Thus, when combining atom pairs,
1325 a fixed list of $N$ restraints may be taken together,
1326 where the apparent ``distance'' is given by:
1328 r_N(t) = \left [\sum_{n=1}^{N} \bar{r}_{n}(t)^{-6} \right]^{-1/6}
1331 where we use $r_{ij}$ or \eqnref{rav} for the $\bar{r}_{n}$.
1332 The $r_N$ of the instantaneous and time-averaged distances
1333 can be combined to do a mixed restraining, as indicated above.
1334 As more pairs of protons contribute to the same NOE signal, the intensity
1335 will increase, and the summed ``distance'' will be shorter than any of
1336 its components due to the reciprocal summation.
1338 There are two options for distributing the forces over the atom pairs.
1339 In the conservative option, the force is defined as the derivative of the
1340 restraint potential with respect to the coordinates. This results in
1341 a conservative potential when time averaging is not used.
1342 The force distribution over the pairs is proportional to $r^{-6}$.
1343 This means that a close pair feels a much larger force than a distant pair,
1344 which might lead to a molecule that is ``too rigid.''
1345 The other option is an equal force distribution. In this case each pair
1346 feels $1/N$ of the derivative of the restraint potential with respect to
1347 $r_N$. The advantage of this method is that more conformations might be
1348 sampled, but the non-conservative nature of the forces can lead to
1349 local heating of the protons.
1351 It is also possible to use {\em ensemble averaging} using multiple
1352 (protein) molecules. In this case the bounds should be lowered as in:
1355 r_1 &~=~& r_1 * M^{-1/6} \\
1356 r_2 &~=~& r_2 * M^{-1/6}
1359 where $M$ is the number of molecules. The {\gromacs} preprocessor {\tt grompp}
1360 can do this automatically when the appropriate option is given.
1361 The resulting ``distance'' is
1362 then used to calculate the scalar force according to:
1366 ~& 0 \hspace{4cm} & r_{N} < r_1 \\
1367 & k_{dr}(r_{N}-r_1)\frac{\rvij}{r_{ij}} & r_1 \le r_{N} < r_2 \\
1368 & k_{dr}(r_2-r_1)\frac{\rvij}{r_{ij}} & r_{N} \ge r_2
1371 where $i$ and $j$ denote the atoms of all the
1372 pairs that contribute to the NOE signal.
1374 \subsubsection{Using distance restraints}
1376 A list of distance restrains based on NOE data can be added to a molecule
1377 definition in your topology file, like in the following example:
1380 [ distance_restraints ]
1381 ; ai aj type index type' low up1 up2 fac
1382 10 16 1 0 1 0.0 0.3 0.4 1.0
1383 10 28 1 1 1 0.0 0.3 0.4 1.0
1384 10 46 1 1 1 0.0 0.3 0.4 1.0
1385 16 22 1 2 1 0.0 0.3 0.4 2.5
1386 16 34 1 3 1 0.0 0.5 0.6 1.0
1389 In this example a number of features can be found. In columns {\tt
1390 ai} and {\tt aj} you find the atom numbers of the particles to be
1391 restrained. The {\tt type} column should always be 1. As explained in
1392 ~\ssecref{distancerestraint}, multiple distances can contribute to a single NOE
1393 signal. In the topology this can be set using the {\tt index}
1394 column. In our example, the restraints 10-28 and 10-46 both have index
1395 1, therefore they are treated simultaneously. An extra requirement
1396 for treating restraints together is that the restraints must be on
1397 successive lines, without any other intervening restraint. The {\tt
1398 type'} column will usually be 1, but can be set to 2 to obtain a
1399 distance restraint that will never be time- and ensemble-averaged;
1400 this can be useful for restraining hydrogen bonds. The columns {\tt
1401 low}, {\tt up1}, and {\tt up2} hold the values of $r_0$, $r_1$, and
1402 $r_2$ from ~\eqnref{disre}. In some cases it can be useful to have
1403 different force constants for some restraints; this is controlled by
1404 the column {\tt fac}. The force constant in the parameter file is
1405 multiplied by the value in the column {\tt fac} for each restraint.
1406 %} % Brace matches ifthenelse test for gmxlite
1408 \newcommand{\SSS}{{\mathbf S}}
1409 \newcommand{\DD}{{\mathbf D}}
1410 \newcommand{\RR}{{\mathbf R}}
1412 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1413 \subsection{Orientation restraints\swapindexquiet{orientation}{restraint}}
1414 \label{subsec:orientationrestraint}
1415 This section describes how orientations between vectors,
1416 as measured in certain NMR experiments, can be calculated
1417 and restrained in MD simulations.
1418 The presented refinement methodology and a comparison of results
1419 with and without time and ensemble averaging have been
1420 published~\cite{Hess2003}.
1421 \subsubsection{Theory}
1422 In an NMR experiment, orientations of vectors can be measured when a
1423 molecule does not tumble completely isotropically in the solvent.
1424 Two examples of such orientation measurements are
1425 residual \normindex{dipolar couplings}
1426 (between two nuclei) or chemical shift anisotropies.
1427 An observable for a vector $\ve{r}_i$ can be written as follows:
1429 \delta_i = \frac{2}{3} \mbox{tr}(\SSS\DD_i)
1431 where $\SSS$ is the dimensionless order tensor of the molecule.
1432 The tensor $\DD_i$ is given by:
1435 \DD_i = \frac{c_i}{\|\ve{r}_i\|^\alpha} \left(
1437 %3 r_x r_x - \ve{r}\cdot\ve{r} & 3 r_x r_y & 3 r_x r_z \\
1438 %3 r_x r_y & 3 r_y r_y - \ve{r}\cdot\ve{r} & 3yz \\
1439 %3 r_x r_z & 3 r_y r_z & 3 r_z r_z - \ve{r}\cdot\ve{r}
1440 %\end{array} \right)
1442 3 x x - 1 & 3 x y & 3 x z \\
1443 3 x y & 3 y y - 1 & 3 y z \\
1444 3 x z & 3 y z & 3 z z - 1 \\
1449 x=\frac{r_{i,x}}{\|\ve{r}_i\|}, \quad
1450 y=\frac{r_{i,y}}{\|\ve{r}_i\|}, \quad
1451 z=\frac{r_{i,z}}{\|\ve{r}_i\|}
1453 For a dipolar coupling $\ve{r}_i$ is the vector connecting the two
1454 nuclei, $\alpha=3$ and the constant $c_i$ is given by:
1456 c_i = \frac{\mu_0}{4\pi} \gamma_1^i \gamma_2^i \frac{\hbar}{4\pi}
1458 where $\gamma_1^i$ and $\gamma_2^i$ are the gyromagnetic ratios of the
1461 The order tensor is symmetric and has trace zero. Using a rotation matrix
1462 ${\mathbf T}$ it can be transformed into the following form:
1464 {\mathbf T}^T \SSS {\mathbf T} = s \left( \begin{array}{ccc}
1465 -\frac{1}{2}(1-\eta) & 0 & 0 \\
1466 0 & -\frac{1}{2}(1+\eta) & 0 \\
1470 where $-1 \leq s \leq 1$ and $0 \leq \eta \leq 1$.
1471 $s$ is called the order parameter and $\eta$ the asymmetry of the
1472 order tensor $\SSS$. When the molecule tumbles isotropically in the
1473 solvent, $s$ is zero, and no orientational effects can be observed
1474 because all $\delta_i$ are zero.
1478 \subsubsection{Calculating orientations in a simulation}
1479 For reasons which are explained below, the $\DD$ matrices are calculated
1480 which respect to a reference orientation of the molecule. The orientation
1481 is defined by a rotation matrix $\RR$, which is needed to least-squares fit
1482 the current coordinates of a selected set of atoms onto
1483 a reference conformation. The reference conformation is the starting
1484 conformation of the simulation. In case of ensemble averaging, which will
1485 be treated later, the structure is taken from the first subsystem.
1486 The calculated $\DD_i^c$ matrix is given by:
1489 \DD_i^c(t) = \RR(t) \DD_i(t) \RR^T(t)
1491 The calculated orientation for vector $i$ is given by:
1493 \delta^c_i(t) = \frac{2}{3} \mbox{tr}(\SSS(t)\DD_i^c(t))
1495 The order tensor $\SSS(t)$ is usually unknown.
1496 A reasonable choice for the order tensor is the tensor
1497 which minimizes the (weighted) mean square difference between the calculated
1498 and the observed orientations:
1501 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
1503 To properly combine different types of measurements, the unit of $w_i$ should
1504 be such that all terms are dimensionless. This means the unit of $w_i$
1505 is the unit of $\delta_i$ to the power $-2$.
1506 {\bf Note} that scaling all $w_i$ with a constant factor does not influence
1509 \subsubsection{Time averaging}
1510 Since the tensors $\DD_i$ fluctuate rapidly in time, much faster than can
1511 be observed in an experiment, they should be averaged over time in the simulation.
1512 However, in a simulation the time and the number of copies of
1513 a molecule are limited. Usually one can not obtain a converged average
1514 of the $\DD_i$ tensors over all orientations of the molecule.
1515 If one assumes that the average orientations of the $\ve{r}_i$ vectors
1516 within the molecule converge much faster than the tumbling time of
1517 the molecule, the tensor can be averaged in an axis system that
1518 rotates with the molecule, as expressed by equation~(\ref{D_rot}).
1519 The time-averaged tensors are calculated
1520 using an exponentially decaying memory function:
1522 \DD^a_i(t) = \frac{\displaystyle
1523 \int_{u=t_0}^t \DD^c_i(u) \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
1525 \int_{u=t_0}^t \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
1528 Assuming that the order tensor $\SSS$ fluctuates slower than the
1529 $\DD_i$, the time-averaged orientation can be calculated as:
1531 \delta_i^a(t) = \frac{2}{3} \mbox{tr}(\SSS(t) \DD_i^a(t))
1533 where the order tensor $\SSS(t)$ is calculated using expression~(\ref{S_msd})
1534 with $\delta_i^c(t)$ replaced by $\delta_i^a(t)$.
1536 \subsubsection{Restraining}
1537 The simulated structure can be restrained by applying a force proportional
1538 to the difference between the calculated and the experimental orientations.
1539 When no time averaging is applied, a proper potential can be defined as:
1541 V = \frac{1}{2} k \sum_{i=1}^N w_i (\delta_i^c (t) -\delta_i^{exp})^2
1543 where the unit of $k$ is the unit of energy.
1544 Thus the effective force constant for restraint $i$ is $k w_i$.
1545 The forces are given by minus the gradient of $V$.
1546 The force $\ve{F}\!_i$ working on vector $\ve{r}_i$ is:
1549 & = & - \frac{\mbox{d} V}{\mbox{d}\ve{r}_i} \\
1550 & = & -k w_i (\delta_i^c (t) -\delta_i^{exp}) \frac{\mbox{d} \delta_i (t)}{\mbox{d}\ve{r}_i} \\
1551 & = & -k w_i (\delta_i^c (t) -\delta_i^{exp})
1552 \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)
1555 \subsubsection{Ensemble averaging}
1556 Ensemble averaging can be applied by simulating a system of $M$ subsystems
1557 that each contain an identical set of orientation restraints. The systems only
1558 interact via the orientation restraint potential which is defined as:
1560 V = M \frac{1}{2} k \sum_{i=1}^N w_i
1561 \langle \delta_i^c (t) -\delta_i^{exp} \rangle^2
1563 The force on vector $\ve{r}_{i,m}$ in subsystem $m$ is given by:
1565 \ve{F}\!_{i,m}(t) = - \frac{\mbox{d} V}{\mbox{d}\ve{r}_{i,m}} =
1566 -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}} \\
1569 \subsubsection{Time averaging}
1570 When using time averaging it is not possible to define a potential.
1571 We can still define a quantity that gives a rough idea of the energy
1572 stored in the restraints:
1574 V = M \frac{1}{2} k^a \sum_{i=1}^N w_i
1575 \langle \delta_i^a (t) -\delta_i^{exp} \rangle^2
1577 The force constant $k_a$ is switched on slowly to compensate for the
1578 lack of history at times close to $t_0$. It is exactly proportional
1579 to the amount of average that has been accumulated:
1582 k \, \frac{1}{\tau}\int_{u=t_0}^t \exp\left(-\frac{t-u}{\tau}\right)\mbox{d} u
1584 What really matters is the definition of the force. It is chosen to
1585 be proportional to the square root of the product of the time-averaged
1586 and the instantaneous deviation.
1587 Using only the time-averaged deviation induces large oscillations.
1588 The force is given by:
1591 %\left\{ \begin{array}{ll}
1592 %0 & \mbox{for} \quad \langle \delta_i^a (t) -\delta_i^{exp} \rangle \langle \delta_i (t) -\delta_i^{exp} \rangle \leq 0 \\
1593 %... & \mbox{for} \quad \langle \delta_i^a (t) -\delta_i^{exp} \rangle \langle \delta_i (t) -\delta_i^{exp} \rangle > 0
1596 \left\{ \begin{array}{ll}
1597 0 & \quad \mbox{for} \quad a\, b \leq 0 \\
1599 k^a w_i \frac{a}{|a|} \sqrt{a\, b} \, \frac{\mbox{d} \delta_{i,m}^c (t)}{\mbox{d}\ve{r}_{i,m}}
1600 & \quad \mbox{for} \quad a\, b > 0
1605 a &=& \langle \delta_i^a (t) -\delta_i^{exp} \rangle \\
1606 b &=& \langle \delta_i^c (t) -\delta_i^{exp} \rangle
1609 \subsubsection{Using orientation restraints}
1610 Orientation restraints can be added to a molecule definition in
1611 the topology file in the section {\tt [~orientation_restraints~]}.
1612 Here we give an example section containing five N-H residual dipolar
1613 coupling restraints:
1616 [ orientation_restraints ]
1617 ; ai aj type exp. label alpha const. obs. weight
1619 31 32 1 1 3 3 6.083 -6.73 1.0
1620 43 44 1 1 4 3 6.083 -7.87 1.0
1621 55 56 1 1 5 3 6.083 -7.13 1.0
1622 65 66 1 1 6 3 6.083 -2.57 1.0
1623 73 74 1 1 7 3 6.083 -2.10 1.0
1626 The unit of the observable is Hz, but one can choose any other unit.
1628 ai} and {\tt aj} you find the atom numbers of the particles to be
1629 restrained. The {\tt type} column should always be 1.
1630 The {\tt exp.} column denotes the experiment number, starting
1631 at 1. For each experiment a separate order tensor $\SSS$
1632 is optimized. The label should be a unique number larger than zero
1633 for each restraint. The {\tt alpha} column contains the power $\alpha$
1634 that is used in equation~(\ref{orient_def}) to calculate the orientation.
1635 The {\tt const.} column contains the constant $c_i$ used in the same
1636 equation. The constant should have the unit of the observable times
1637 nm$^\alpha$. The column {\tt obs.} contains the observable, in any
1638 unit you like. The last column contains the weights $w_i$; the unit
1639 should be the inverse of the square of the unit of the observable.
1641 Some parameters for orientation restraints can be specified in the
1642 {\tt grompp.mdp} file, for a study of the effect of different
1643 force constants and averaging times and ensemble averaging see~\cite{Hess2003}.
1644 %} % Brace matches ifthenelse test for gmxlite
1646 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1647 \section{Polarization}
1648 Polarization can be treated by {\gromacs} by attaching
1649 \normindex{shell} (\normindex{Drude}) particles to atoms and/or
1650 virtual sites. The energy of the shell particle is then minimized at
1651 each time step in order to remain on the Born-Oppenheimer surface.
1653 \subsection{Simple polarization}
1654 This is merely a harmonic potential with equilibrium distance 0.
1656 \subsection{Water polarization}
1657 A special potential for water that allows anisotropic polarization of
1658 a single shell particle~\cite{Maaren2001a}.
1660 \subsection{Thole polarization}
1661 Based on early work by \normindex{Thole}~\cite{Thole81}, Roux and
1662 coworkers have implemented potentials for molecules like
1663 ethanol~\cite{Lamoureux2003a,Lamoureux2003b,Noskov2005a}. Within such
1664 molecules, there are intra-molecular interactions between shell
1665 particles, however these must be screened because full Coulomb would
1666 be too strong. The potential between two shell particles $i$ and $j$ is:
1667 \newcommand{\rbij}{\bar{r}_{ij}}
1669 V_{thole} ~=~ \frac{q_i q_j}{r_{ij}}\left[1-\left(1+\frac{\rbij}{2}\right){\rm exp}^{-\rbij}\right]
1671 {\bf Note} that there is a sign error in Equation~1 of Noskov {\em et al.}~\cite{Noskov2005a}:
1673 \rbij ~=~ a\frac{r_{ij}}{(\alpha_i \alpha_j)^{1/6}}
1675 where $a$ is a magic (dimensionless) constant, usually chosen to be
1676 2.6~\cite{Noskov2005a}; $\alpha_i$ and $\alpha_j$ are the polarizabilities
1677 of the respective shell particles.
1679 %} % Brace matches ifthenelse test for gmxlite
1681 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1682 \section{Free energy interactions}
1684 \index{free energy interactions}
1685 \newcommand{\LAM}{\lambda}
1686 \newcommand{\LL}{(1-\LAM)}
1687 \newcommand{\dvdl}[1]{\frac{\partial #1}{\partial \LAM}}
1688 This section describes the $\lambda$-dependence of the potentials
1689 used for free energy calculations (see \secref{fecalc}).
1690 All common types of potentials and constraints can be
1691 interpolated smoothly from state A ($\lambda=0$) to state B
1692 ($\lambda=1$) and vice versa.
1693 All bonded interactions are interpolated by linear interpolation
1694 of the interaction parameters. Non-bonded interactions can be
1695 interpolated linearly or via soft-core interactions.
1697 Starting in {\gromacs} 4.6, $\lambda$ is a vector, allowing different
1698 components of the free energy transformation to be carried out at
1699 different rates. Coulomb, Lennard-Jones, bonded, and restraint terms
1700 can all be controlled independently, as described in the {\tt .mdp}
1703 \subsubsection{Harmonic potentials}
1704 The example given here is for the bond potential, which is harmonic
1705 in {\gromacs}. However, these equations apply to the angle potential
1706 and the improper dihedral potential as well.
1708 V_b &=&\half\left[\LL k_b^A +
1709 \LAM k_b^B\right] \left[b - \LL b_0^A - \LAM b_0^B\right]^2 \\
1710 \dvdl{V_b}&=&\half(k_b^B-k_b^A)
1711 \left[b - \LL b_0^A + \LAM b_0^B\right]^2 +
1713 & & \phantom{\half}(b_0^A-b_0^B) \left[b - \LL b_0^A -\LAM b_0^B\right]
1714 \left[\LL k_b^A + \LAM k_b^B \right]
1717 \subsubsection{\gromosv{96} bonds and angles}
1718 Fourth-power bond stretching and cosine-based angle potentials
1719 are interpolated by linear interpolation of the force constant
1720 and the equilibrium position. Formulas are not given here.
1722 \subsubsection{Proper dihedrals}
1723 For the proper dihedrals, the equations are somewhat more complicated:
1725 V_d &=&\left[\LL k_d^A + \LAM k_d^B \right]
1726 \left( 1+ \cos\left[n_{\phi} \phi -
1727 \LL \phi_s^A - \LAM \phi_s^B
1729 \dvdl{V_d}&=&(k_d^B-k_d^A)
1732 n_{\phi} \phi- \LL \phi_s^A - \LAM \phi_s^B
1736 &&(\phi_s^B - \phi_s^A) \left[\LL k_d^A - \LAM k_d^B\right]
1737 \sin\left[ n_{\phi}\phi - \LL \phi_s^A - \LAM \phi_s^B \right]
1739 {\bf Note:} that the multiplicity $n_{\phi}$ can not be parameterized
1740 because the function should remain periodic on the interval $[0,2\pi]$.
1742 \subsubsection{Tabulated bonded interactions}
1743 For tabulated bonded interactions only the force constant can interpolated:
1745 V &=& (\LL k^A + \LAM k^B) \, f \\
1746 \dvdl{V} &=& (k^B - k^A) \, f
1749 \subsubsection{Coulomb interaction}
1750 The \normindex{Coulomb} interaction between two particles
1751 of which the charge varies with $\LAM$ is:
1753 V_c &=& \frac{f}{\epsrf \rij}\left[\LL q_i^A q_j^A + \LAM\, q_i^B q_j^B\right] \\
1754 \dvdl{V_c}&=& \frac{f}{\epsrf \rij}\left[- q_i^A q_j^A + q_i^B q_j^B\right]
1756 where $f = \frac{1}{4\pi \varepsilon_0} = 138.935\,485$ (see \chref{defunits}).
1758 \subsubsection{Coulomb interaction with \normindex{reaction field}}
1759 The Coulomb interaction including a reaction field, between two particles
1760 of which the charge varies with $\LAM$ is:
1762 V_c &=& f\left[\frac{1}{\rij} + k_{rf}~ \rij^2 -c_{rf}\right]
1763 \left[\LL q_i^A q_j^A + \LAM\, q_i^B q_j^B\right] \\
1764 \dvdl{V_c}&=& f\left[\frac{1}{\rij} + k_{rf}~ \rij^2 -c_{rf}\right]
1765 \left[- q_i^A q_j^A + q_i^B q_j^B\right]
1766 \label{eq:dVcoulombdlambda}
1768 {\bf Note} that the constants $k_{rf}$ and $c_{rf}$ are
1769 defined using the dielectric
1770 constant $\epsrf$ of the medium (see \secref{coulrf}).
1772 \subsubsection{Lennard-Jones interaction}
1773 For the \normindex{Lennard-Jones} interaction between two particles
1774 of which the {\em atom type} varies with $\LAM$ we can write:
1776 V_{LJ} &=& \frac{\LL C_{12}^A + \LAM\, C_{12}^B}{\rij^{12}} -
1777 \frac{\LL C_6^A + \LAM\, C_6^B}{\rij^6} \\
1778 \dvdl{V_{LJ}}&=&\frac{C_{12}^B - C_{12}^A}{\rij^{12}} -
1779 \frac{C_6^B - C_6^A}{\rij^6}
1780 \label{eq:dVljdlambda}
1782 It should be noted that it is also possible to express a pathway from
1783 state A to state B using $\sigma$ and $\epsilon$ (see \eqnref{sigeps}).
1784 It may seem to make sense physically to vary the force field parameters
1785 $\sigma$ and $\epsilon$ rather
1786 than the derived parameters $C_{12}$ and $C_{6}$.
1787 However, the difference between the pathways in parameter space
1788 is not large, and the free energy itself
1789 does not depend on the pathway, so we use the simple formulation
1792 \subsubsection{Kinetic Energy}
1793 When the mass of a particle changes, there is also a contribution of
1794 the kinetic energy to the free energy (note that we can not write
1795 the momentum \ve{p} as m\ve{v}, since that would result
1796 in the sign of $\dvdl{E_k}$ being incorrect~\cite{Gunsteren98a}):
1799 E_k &=& \half\frac{\ve{p}^2}{\LL m^A + \LAM m^B} \\
1800 \dvdl{E_k}&=& -\half\frac{\ve{p}^2(m^B-m^A)}{(\LL m^A + \LAM m^B)^2}
1802 after taking the derivative, we {\em can} insert \ve{p} = m\ve{v}, such that:
1804 \dvdl{E_k}~=~ -\half\ve{v}^2(m^B-m^A)
1807 \subsubsection{Constraints}
1808 \label{subsubsec:constraints}
1809 \newcommand{\clam}{C_{\lambda}}
1810 The constraints are formally part of the Hamiltonian, and therefore
1811 they give a contribution to the free energy. In {\gromacs} this can be
1812 calculated using the \normindex{LINCS} or the \normindex{SHAKE} algorithm.
1813 If we have a number of constraint equations $g_k$:
1815 g_k = \ve{r}_{k} - d_{k}
1817 where $\ve{r}_k$ is the distance vector between two particles and
1818 $d_k$ is the constraint distance between the two particles, we can write
1819 this using a $\LAM$-dependent distance as
1821 g_k = \ve{r}_{k} - \left(\LL d_{k}^A + \LAM d_k^B\right)
1823 the contribution $\clam$
1824 to the Hamiltonian using Lagrange multipliers $\lambda$:
1826 \clam &=& \sum_k \lambda_k g_k \\
1827 \dvdl{\clam} &=& \sum_k \lambda_k \left(d_k^B-d_k^A\right)
1831 \subsection{Soft-core interactions\index{soft-core interactions}}
1833 \centerline{\includegraphics[height=6cm]{plots/softcore}}
1834 \caption{Soft-core interactions at $\LAM=0.5$, with $p=2$ and
1835 $C_6^A=C_{12}^A=C_6^B=C_{12}^B=1$.}
1836 \label{fig:softcore}
1838 In a free-energy calculation where particles grow out of nothing, or
1839 particles disappear, using the the simple linear interpolation of the
1840 Lennard-Jones and Coulomb potentials as described in Equations~\ref{eq:dVljdlambda}
1841 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
1842 close to each other, leading to large fluctuations in the measured values of
1843 $\partial V/\partial \LAM$ (which, because of the simple linear
1844 interpolation, depends on the potentials at both the endpoints of $\LAM$).
1846 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
1847 involved at $\LAM$ values between 0 and 1, but not \emph{at} $\LAM=0$
1850 In {\gromacs} the soft-core potentials $V_{sc}$ are shifted versions of the
1851 regular potentials, so that the singularity in the potential and its
1852 derivatives at $r=0$ is never reached:
1854 V_{sc}(r) &=& \LL V^A(r_A) + \LAM V^B(r_B)
1856 r_A &=& \left(\alpha \sigma_A^6 \LAM^p + r^6 \right)^\frac{1}{6}
1858 r_B &=& \left(\alpha \sigma_B^6 \LL^p + r^6 \right)^\frac{1}{6}
1860 where $V^A$ and $V^B$ are the normal ``hard core'' Van der Waals or
1861 electrostatic potentials in state A ($\LAM=0$) and state B ($\LAM=1$)
1862 respectively, $\alpha$ is the soft-core parameter (set with {\tt sc_alpha}
1863 in the {\tt .mdp} file), $p$ is the soft-core $\LAM$ power (set with
1864 {\tt sc_power}), $\sigma$ is the radius of the interaction, which is
1865 $(C_{12}/C_6)^{1/6}$ or an input parameter ({\tt sc_sigma}) when $C_6$
1866 or $C_{12}$ is zero.
1868 For intermediate $\LAM$, $r_A$ and $r_B$ alter the interactions very little
1869 for $r > \alpha^{1/6} \sigma$ and quickly switch the soft-core
1870 interaction to an almost constant value for smaller $r$ (\figref{softcore}).
1873 F_{sc}(r) = -\frac{\partial V_{sc}(r)}{\partial r} =
1874 \LL F^A(r_A) \left(\frac{r}{r_A}\right)^5 +
1875 \LAM F^B(r_B) \left(\frac{r}{r_B}\right)^5
1877 where $F^A$ and $F^B$ are the ``hard core'' forces.
1878 The contribution to the derivative of the free energy is:
1880 \dvdl{V_{sc}(r)} & = &
1881 V^B(r_B) -V^A(r_A) +
1882 \LL \frac{\partial V^A(r_A)}{\partial r_A}
1883 \frac{\partial r_A}{\partial \LAM} +
1884 \LAM\frac{\partial V^B(r_B)}{\partial r_B}
1885 \frac{\partial r_B}{\partial \LAM}
1888 V^B(r_B) -V^A(r_A) + \nonumber \\
1891 \left[ \LAM F^B(r_B) r^{-5}_B \sigma_B^6 \LL^{p-1} -
1892 \LL F^A(r_A) r^{-5}_A \sigma_A^6 \LAM^{p-1} \right]
1895 The original GROMOS Lennard-Jones soft-core function~\cite{Beutler94}
1896 uses $p=2$, but $p=1$ gives a smoother $\partial H/\partial\LAM$ curve.
1897 %When the changes between the two states involve both the disappearing
1898 %and appearing of atoms, it is important that the overlapping of atoms
1899 %happens around $\LAM=0.5$. This can usually be achieved with
1900 %$\alpha$$\approx0.7$ for $p=1$ and $\alpha$$\approx1.5$ for $p=2$.
1901 %MRS: this is now eliminated as of 4.6, since changes between atoms are done linearly.
1903 Another issue that should be considered is the soft-core effect of hydrogens
1904 without Lennard-Jones interaction. Their soft-core $\sigma$ is
1905 set with {\tt sc-sigma} in the {\tt .mdp} file. These hydrogens
1906 produce peaks in $\partial H/\partial\LAM$ at $\LAM$ is 0 and/or 1 for $p=1$
1907 and close to 0 and/or 1 with $p=2$. Lowering {\tt\mbox{sc-sigma}} will decrease
1908 this effect, but it will also increase the interactions with hydrogens
1909 relative to the other interactions in the soft-core state.
1911 When soft core potentials are selected (by setting {\tt sc-alpha} \textgreater
1912 0), and the Coulomb and Lennard-Jones potentials are turned on or off
1913 sequentially, then the Coulombic interaction is turned off linearly,
1914 rather than using soft core interactions, which should be less
1915 statistically noisy in most cases. This behavior can be overwritten
1916 by using the mdp option {\tt sc-coul} to {\tt yes}. Additionally, the
1917 soft-core interaction potential is only applied when either the A or B
1918 state has zero interaction potential. If both A and B states have
1919 nonzero interaction potential, default linear scaling described above
1920 is used. When both Coulombic and Lennard-Jones interactions are turned
1921 off simultaneously, a soft-core potential is used, and a hydrogen is
1922 being introduced or deleted, the sigma is set to {\tt sc-sigma-min},
1923 which itself defaults to {\tt sc-sigma-default}.
1925 Recently, a new formulation of the soft-core approach has been derived
1926 that in most cases gives lower and more even statistical variance than
1927 the standard soft-core path described above.~\cite{Pham2011,Pham2012}
1928 Specifically, we have:
1930 V_{sc}(r) &=& \LL V^A(r_A) + \LAM V^B(r_B)
1932 r_A &=& \left(\alpha \sigma_A^{48} \LAM^p + r^{48} \right)^\frac{1}{48}
1934 r_B &=& \left(\alpha \sigma_B^{48} \LL^p + r^{48} \right)^\frac{1}{48}
1936 This ``1-1-48'' path is also implemented in {\gromacs}. Note that for this path the soft core $\alpha$
1937 should satisfy $0.001 < \alpha < 0.003$, rather than $\alpha \approx
1940 %} % Brace matches ifthenelse test for gmxlite
1942 %\ifthenelse{\equal{\gmxlite}{1}}{}{
1944 \subsection{Exclusions and 1-4 Interactions.}
1945 Atoms within a molecule that are close by in the chain,
1946 {\ie} atoms that are covalently bonded, or linked by one or two
1947 atoms are called {\em first neighbors, second neighbors} and
1948 {\em \swapindex{third}{neighbor}s}, respectively (see \figref{chain}).
1949 Since the interactions of atom {\bf i} with atoms {\bf i+1} and {\bf i+2}
1950 are mainly quantum mechanical, they can not be modeled by a Lennard-Jones potential.
1951 Instead it is assumed that these interactions are adequately modeled
1952 by a harmonic bond term or constraint ({\bf i, i+1}) and a harmonic angle term
1953 ({\bf i, i+2}). The first and second neighbors (atoms {\bf i+1} and {\bf i+2})
1955 {\em excluded} from the Lennard-Jones \swapindex{interaction}{list}
1957 atoms {\bf i+1} and {\bf i+2} are called {\em \normindex{exclusions}} of atom {\bf i}.
1960 \centerline{\includegraphics[width=8cm]{plots/chain}}
1961 \caption{Atoms along an alkane chain.}
1965 For third neighbors, the normal Lennard-Jones repulsion is sometimes
1966 still too strong, which means that when applied to a molecule, the
1967 molecule would deform or break due to the internal strain. This is
1968 especially the case for carbon-carbon interactions in a {\em
1969 cis}-conformation ({\eg} {\em cis}-butane). Therefore, for some of these
1970 interactions, the Lennard-Jones repulsion has been reduced in the
1971 {\gromos} force field, which is implemented by keeping a separate list of
1972 1-4 and normal Lennard-Jones parameters. In other force fields, such
1973 as OPLS~\cite{Jorgensen88}, the standard Lennard-Jones parameters are reduced
1974 by a factor of two, but in that case also the dispersion (r$^{-6}$)
1975 and the Coulomb interaction are scaled.
1976 {\gromacs} can use either of these methods.
1978 \subsection{Charge Groups\index{charge group}}
1980 In principle, the force calculation in MD is an $O(N^2)$ problem.
1981 Therefore, we apply a \normindex{cut-off} for non-bonded force (NBF)
1982 calculations; only the particles within a certain distance of each
1983 other are interacting. This reduces the cost to $O(N)$ (typically
1984 $100N$ to $200N$) of the NBF. It also introduces an error, which is,
1985 in most cases, acceptable, except when applying the cut-off implies
1986 the creation of charges, in which case you should consider using the
1987 lattice sum methods provided by {\gromacs}.
1989 Consider a water molecule interacting with another atom. If we would apply
1990 a plain cut-off on an atom-atom basis we might include the atom-oxygen
1991 interaction (with a charge of $-0.82$) without the compensating charge
1992 of the protons, and as a result, induce a large dipole moment over the system.
1993 Therefore, we have to keep groups of atoms with total charge
1994 0 together. These groups are called {\em charge groups}. Note that with
1995 a proper treatment of long-range electrostatics (e.g. particle-mesh Ewald
1996 (\secref{pme}), keeping charge groups together is not required.
1998 \subsection{Treatment of Cut-offs in the group scheme\index{cut-off}}
1999 \newcommand{\rs}{$r_{short}$}
2000 \newcommand{\rl}{$r_{long}$}
2001 {\gromacs} is quite flexible in treating cut-offs, which implies
2002 there can be quite a number of parameters to set. These parameters are
2003 set in the input file for {\tt grompp}. There are two sort of parameters
2004 that affect the cut-off interactions; you can select which type
2005 of interaction to use in each case, and which cut-offs should be
2006 used in the neighbor searching.
2008 For both Coulomb and van der Waals interactions there are interaction
2009 type selectors (termed {\tt vdwtype} and {\tt coulombtype}) and two
2010 parameters, for a total of six non-bonded interaction parameters. See
2011 \secref{mdpopt} for a complete description of these parameters.
2013 The neighbor searching (NS) can be performed using a single-range, or a twin-range
2014 approach. Since the former is merely a special case of the latter, we will
2015 discuss the more general twin-range. In this case, NS is described by two
2016 radii: {\tt rlist} and max({\tt rcoulomb},{\tt rvdw}).
2017 Usually one builds the neighbor list every 10 time steps
2018 or every 20 fs (parameter {\tt nstlist}). In the neighbor list, all interaction
2019 pairs that fall within {\tt rlist} are stored. Furthermore, the
2020 interactions between pairs that do not
2021 fall within {\tt rlist} but do fall within max({\tt rcoulomb},{\tt rvdw})
2022 are computed during NS. The
2023 forces and energy are stored separately and added to short-range forces
2024 at every time step between successive NS. If {\tt rlist} =
2025 max({\tt rcoulomb},{\tt rvdw}), no forces
2026 are evaluated during neighbor list generation.
2027 The \normindex{virial} is calculated from the sum of the short- and
2029 This means that the virial can be slightly asymmetrical at non-NS steps.
2030 In single precision, the virial is almost always asymmetrical because the
2031 off-diagonal elements are about as large as each element in the sum.
2032 In most cases this is not really a problem, since the fluctuations in the
2033 virial can be 2 orders of magnitude larger than the average.
2035 Except for the plain cut-off,
2036 all of the interaction functions in \tabref{funcparm}
2037 require that neighbor searching be done with a larger radius than the $r_c$
2038 specified for the functional form, because of the use of charge groups.
2039 The extra radius is typically of the order of 0.25 nm (roughly the
2040 largest distance between two atoms in a charge group plus the distance a
2041 charge group can diffuse within neighbor list updates).
2043 %If your charge groups are very large it may be interesting to turn off charge
2044 %groups, by setting the option
2045 %{\tt bAtomList = yes} in your {\tt grompp.mdp} file.
2046 %In this case only a small extra radius to account for diffusion needs to be
2047 %added (0.1 nm). Do not however use this together with the plain cut-off
2048 %method, as it will generate large artifacts (\secref{cg}).
2049 %In summary, there are four parameters that describe NS behavior:
2050 %{\tt nstlist} (update frequency in number of time steps),
2051 %{\tt bAtomList} (whether or not charge groups are used to generate neighbor list, the default is to use charge groups, so {\tt bAtomList = no}),
2052 %{\tt rshort} and {\tt rlong} which are the two radii {\rs} and {\rl}
2057 \begin{tabular}{|ll|l|}
2059 \multicolumn{2}{|c|}{Type} & Parameters \\
2061 Coulomb&Plain cut-off & $r_c$, $\epsr$ \\
2062 &Reaction field & $r_c$, $\epsrf$ \\
2063 &Shift function & $r_1$, $r_c$, $\epsr$ \\
2064 &Switch function & $r_1$, $r_c$, $\epsr$ \\
2066 VdW&Plain cut-off & $r_c$ \\
2067 &Shift function & $r_1$, $r_c$ \\
2068 &Switch function & $r_1$, $r_c$ \\
2071 \caption[Parameters for the different functional forms of the
2072 non-bonded interactions.]{Parameters for the different functional
2073 forms of the non-bonded interactions.}
2074 \label{tab:funcparm}
2076 %} % Brace matches ifthenelse test for gmxlite
2079 \newcommand{\vvis}{\ve{r}_s}
2080 \newcommand{\Fi}{\ve{F}_i'}
2081 \newcommand{\Fj}{\ve{F}_j'}
2082 \newcommand{\Fk}{\ve{F}_k'}
2083 \newcommand{\Fl}{\ve{F}_l'}
2084 \newcommand{\Fvis}{\ve{F}_{s}}
2085 \newcommand{\rvik}{\ve{r}_{ik}}
2086 \newcommand{\rvis}{\ve{r}_{is}}
2087 \newcommand{\rvjk}{\ve{r}_{jk}}
2088 \newcommand{\rvjl}{\ve{r}_{jl}}
2090 %\ifthenelse{\equal{\gmxlite}{1}}{}{
2091 \section{Virtual interaction sites\index{virtual interaction sites}}
2092 \label{sec:virtual_sites}
2093 Virtual interaction sites (called \seeindex{dummy atoms}{virtual interaction sites} in {\gromacs} versions before 3.3)
2094 can be used in {\gromacs} in a number of ways.
2095 We write the position of the virtual site $\ve{r}_s$ as a function of
2096 the positions of other particles \ve{r}$_i$: $\ve{r}_s =
2097 f(\ve{r}_1..\ve{r}_n)$. The virtual site, which may carry charge or be
2098 involved in other interactions, can now be used in the force
2099 calculation. The force acting on the virtual site must be
2100 redistributed over the particles with mass in a consistent way.
2101 A good way to do this can be found in ref.~\cite{Berendsen84b}.
2102 We can write the potential energy as:
2104 V = V(\vvis,\ve{r}_1,\ldots,\ve{r}_n) = V^*(\ve{r}_1,\ldots,\ve{r}_n)
2106 The force on the particle $i$ is then:
2108 \ve{F}_i = -\frac{\partial V^*}{\partial \ve{r}_i}
2109 = -\frac{\partial V}{\partial \ve{r}_i} -
2110 \frac{\partial V}{\partial \vvis}
2111 \frac{\partial \vvis}{\partial \ve{r}_i}
2112 = \ve{F}_i^{direct} + \Fi
2114 The first term is the normal force.
2115 The second term is the force on particle $i$ due to the virtual site, which
2116 can be written in tensor notation:
2117 \newcommand{\partd}[2]{\displaystyle\frac{\partial #1}{\partial #2_i}}
2119 \Fi = \left[\begin{array}{ccc}
2120 \partd{x_s}{x} & \partd{y_s}{x} & \partd{z_s}{x} \\[1ex]
2121 \partd{x_s}{y} & \partd{y_s}{y} & \partd{z_s}{y} \\[1ex]
2122 \partd{x_s}{z} & \partd{y_s}{z} & \partd{z_s}{z}
2123 \end{array}\right]\Fvis
2126 where $\Fvis$ is the force on the virtual site and $x_s$, $y_s$ and
2127 $z_s$ are the coordinates of the virtual site. In this way, the total
2128 force and the total torque are conserved~\cite{Berendsen84b}.
2130 The computation of the \normindex{virial}
2131 (\eqnref{Xi}) is non-trivial when virtual sites are used. Since the
2132 virial involves a summation over all the atoms (rather than virtual
2133 sites), the forces must be redistributed from the virtual sites to the
2134 atoms (using ~\eqnref{fvsite}) {\em before} computation of the
2135 virial. In some special cases where the forces on the atoms can be
2136 written as a linear combination of the forces on the virtual sites (types 2
2137 and 3 below) there is no difference between computing the virial
2138 before and after the redistribution of forces. However, in the
2139 general case redistribution should be done first.
2142 \centerline{\includegraphics[width=15cm]{plots/dummies}}
2143 \caption[Virtual site construction.]{The six different types of virtual
2144 site construction in \protect{\gromacs}. The constructing atoms are
2145 shown as black circles, the virtual sites in gray.}
2149 There are six ways to construct virtual sites from surrounding atoms in
2150 {\gromacs}, which we classify by the number of constructing
2151 atoms. {\bf Note} that all site types mentioned can be constructed from
2152 types 3fd (normalized, in-plane) and 3out (non-normalized, out of
2153 plane). However, the amount of computation involved increases sharply
2154 along this list, so we strongly recommended using the first adequate
2155 virtual site type that will be sufficient for a certain purpose.
2156 \figref{vsites} depicts 6 of the available virtual site constructions.
2157 The conceptually simplest construction types are linear combinations:
2159 \vvis = \sum_{i=1}^N w_i \, \ve{r}_i
2161 The force is then redistributed using the same weights:
2166 The types of virtual sites supported in {\gromacs} are given in the list below.
2167 Constructing atoms in virtual sites can be virtual sites themselves, but
2168 only if they are higher in the list, i.e. virtual sites can be
2169 constructed from ``particles'' that are simpler virtual sites.
2171 \item[{\bf\sf 2.}]\label{subsec:vsite2}As a linear combination of two atoms
2172 (\figref{vsites} 2):
2174 w_i = 1 - a ~,~~ w_j = a
2176 In this case the virtual site is on the line through atoms $i$ and
2179 \item[{\bf\sf 3.}]\label{subsec:vsite3}As a linear combination of three atoms
2180 (\figref{vsites} 3):
2182 w_i = 1 - a - b ~,~~ w_j = a ~,~~ w_k = b
2184 In this case the virtual site is in the plane of the other three
2187 \item[{\bf\sf 3fd.}]\label{subsec:vsite3fd}In the plane of three atoms, with a fixed distance
2188 (\figref{vsites} 3fd):
2190 \vvis ~=~ \ve{r}_i + b \frac{ \rvij + a \rvjk }
2191 {| \rvij + a \rvjk |}
2193 In this case the virtual site is in the plane of the other three
2194 particles at a distance of $|b|$ from $i$.
2195 The force on particles $i$, $j$ and $k$ due to the force on the virtual
2196 site can be computed as:
2199 \Fi &=& \displaystyle \Fvis - \gamma ( \Fvis - \ve{p} ) \\[1ex]
2200 \Fj &=& \displaystyle (1-a)\gamma (\Fvis - \ve{p}) \\[1ex]
2201 \Fk &=& \displaystyle a \gamma (\Fvis - \ve{p}) \\
2205 \displaystyle \gamma = \frac{b}{| \rvij + a \rvjk |} \\[2ex]
2206 \displaystyle \ve{p} = \frac{ \rvis \cdot \Fvis }
2207 { \rvis \cdot \rvis } \rvis
2211 \item[{\bf\sf 3fad.}]\label{subsec:vsite3fad}In the plane of three atoms, with a fixed angle and
2212 distance (\figref{vsites} 3fad):
2214 \label{eqn:vsite2fad-F}
2215 \vvis ~=~ \ve{r}_i +
2216 d \cos \theta \frac{\rvij}{|\rvij|} +
2217 d \sin \theta \frac{\ve{r}_\perp}{|\ve{r}_\perp|}
2219 \ve{r}_\perp ~=~ \rvjk -
2220 \frac{ \rvij \cdot \rvjk }
2221 { \rvij \cdot \rvij }
2224 In this case the virtual site is in the plane of the other three
2225 particles at a distance of $|d|$ from $i$ at an angle of
2226 $\alpha$ with $\rvij$. Atom $k$ defines the plane and the
2227 direction of the angle. {\bf Note} that in this case $b$ and
2228 $\alpha$ must be specified, instead of $a$ and $b$ (see also
2229 \secref{vsitetop}). The force on particles $i$, $j$ and $k$
2230 due to the force on the virtual site can be computed as (with
2231 $\ve{r}_\perp$ as defined in \eqnref{vsite2fad-F}):
2232 \newcommand{\dfrac}{\displaystyle\frac}
2235 \begin{array}{lclllll}
2237 \dfrac{d \cos \theta}{|\rvij|} \ve{F}_1 &+&
2238 \dfrac{d \sin \theta}{|\ve{r}_\perp|} \left(
2239 \dfrac{ \rvij \cdot \rvjk }
2240 { \rvij \cdot \rvij } \ve{F}_2 +
2241 \ve{F}_3 \right) \\[3ex]
2243 \dfrac{d \cos \theta}{|\rvij|} \ve{F}_1 &-&
2244 \dfrac{d \sin \theta}{|\ve{r}_\perp|} \left(
2246 \dfrac{ \rvij \cdot \rvjk }
2247 { \rvij \cdot \rvij } \ve{F}_2 +
2248 \ve{F}_3 \right) \\[3ex]
2250 \dfrac{d \sin \theta}{|\ve{r}_\perp|} \ve{F}_2 \\[3ex]
2254 \dfrac{ \rvij \cdot \Fvis }
2255 { \rvij \cdot \rvij } \rvij
2257 \ve{F}_2 = \ve{F}_1 -
2258 \dfrac{ \ve{r}_\perp \cdot \Fvis }
2259 { \ve{r}_\perp \cdot \ve{r}_\perp } \ve{r}_\perp
2261 \ve{F}_3 = \dfrac{ \rvij \cdot \Fvis }
2262 { \rvij \cdot \rvij } \ve{r}_\perp
2266 \item[{\bf\sf 3out.}]\label{subsec:vsite3out}As a non-linear combination of three atoms, out of plane
2267 (\figref{vsites} 3out):
2269 \vvis ~=~ \ve{r}_i + a \rvij + b \rvik +
2270 c (\rvij \times \rvik)
2272 This enables the construction of virtual sites out of the plane of the
2274 The force on particles $i,j$ and $k$ due to the force on the virtual
2275 site can be computed as:
2279 \Fj &=& \left[\begin{array}{ccc}
2280 a & -c\,z_{ik} & c\,y_{ik} \\[0.5ex]
2281 c\,z_{ik} & a & -c\,x_{ik} \\[0.5ex]
2282 -c\,y_{ik} & c\,x_{ik} & a
2283 \end{array}\right]\Fvis \\
2285 \Fk &=& \left[\begin{array}{ccc}
2286 b & c\,z_{ij} & -c\,y_{ij} \\[0.5ex]
2287 -c\,z_{ij} & b & c\,x_{ij} \\[0.5ex]
2288 c\,y_{ij} & -c\,x_{ij} & b
2289 \end{array}\right]\Fvis \\
2290 \Fi &=& \Fvis - \Fj - \Fk
2294 \item[{\bf\sf 4fdn.}]\label{subsec:vsite4fdn}From four atoms, with a fixed distance, see separate Fig.\ \ref{fig:vsite-4fdn}.
2295 This construction is a bit
2296 complex, in particular since the previous type (4fd) could be unstable which forced us
2297 to introduce a more elaborate construction:
2300 \centerline{\includegraphics[width=5cm]{plots/vsite-4fdn}}
2301 \caption {The new 4fdn virtual site construction, which is stable even when all constructing
2302 atoms are in the same plane.}
2303 \label{fig:vsite-4fdn}
2307 \mathbf{r}_{ja} &=& a\, \mathbf{r}_{ik} - \mathbf{r}_{ij} = a\, (\mathbf{x}_k - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\
2308 \mathbf{r}_{jb} &=& b\, \mathbf{r}_{il} - \mathbf{r}_{ij} = b\, (\mathbf{x}_l - \mathbf{x}_i) - (\mathbf{x}_j - \mathbf{x}_i) \nonumber \\
2309 \mathbf{r}_m &=& \mathbf{r}_{ja} \times \mathbf{r}_{jb} \nonumber \\
2310 \mathbf{x}_s &=& \mathbf{x}_i + c \frac{\mathbf{r}_m}{|\mathbf{r}_m|}
2314 In this case the virtual site is at a distance of $|c|$ from $i$, while $a$ and $b$ are
2315 parameters. {\bf Note} that the vectors $\mathbf{r}_{ik}$ and $\mathbf{r}_{ij}$ are not normalized
2316 to save floating-point operations.
2317 The force on particles $i$, $j$, $k$ and $l$ due to the force
2318 on the virtual site are computed through chain rule derivatives
2319 of the construction expression. This is exact and conserves energy,
2320 but it does lead to relatively lengthy expressions that we do not
2321 include here (over 200 floating-point operations). The interested reader can
2322 look at the source code in \verb+vsite.c+. Fortunately, this vsite type is normally
2323 only used for chiral centers such as $C_{\alpha}$ atoms in proteins.
2325 The new 4fdn construct is identified with a `type' value of 2 in the topology. The earlier 4fd
2326 type is still supported internally (`type' value 1), but it should not be used for
2327 new simulations. All current {\gromacs} tools will automatically generate type 4fdn instead.
2330 \item[{\bf\sf N.}]\label{subsec:vsiteN} A linear combination of $N$ atoms with relative
2331 weights $a_i$. The weight for atom $i$ is:
2333 w_i = a_i \left(\sum_{j=1}^N a_j \right)^{-1}
2335 There are three options for setting the weights:
2337 \item[COG] center of geometry: equal weights
2338 \item[COM] center of mass: $a_i$ is the mass of atom $i$;
2339 when in free-energy simulations the mass of the atom is changed,
2340 only the mass of the A-state is used for the weight
2341 \item[COW] center of weights: $a_i$ is defined by the user
2345 %} % Brace matches ifthenelse test for gmxlite
2347 \newcommand{\dr}{{\rm d}r}
2348 \newcommand{\avcsix}{\left< C_6 \right>}
2350 %\ifthenelse{\equal{\gmxlite}{1}}{}{
2351 \section{Long Range Electrostatics}
2352 \label{sec:lr_elstat}
2353 \subsection{Ewald summation\index{Ewald sum}}
2355 The total electrostatic energy of $N$ particles and their periodic
2356 images\index{periodic boundary conditions} is given by
2358 V=\frac{f}{2}\sum_{n_x}\sum_{n_y}
2359 \sum_{n_{z}*} \sum_{i}^{N} \sum_{j}^{N}
2360 \frac{q_i q_j}{{\bf r}_{ij,{\bf n}}}.
2361 \label{eqn:totalcoulomb}
2363 $(n_x,n_y,n_z)={\bf n}$ is the box index vector, and the star indicates that
2364 terms with $i=j$ should be omitted when $(n_x,n_y,n_z)=(0,0,0)$. The
2365 distance ${\bf r}_{ij,{\bf n}}$ is the real distance between the charges and
2366 not the minimum-image. This sum is conditionally convergent, but
2369 Ewald summation was first introduced as a method to calculate
2370 long-range interactions of the periodic images in
2371 crystals~\cite{Ewald21}. The idea is to convert the single
2372 slowly-converging sum \eqnref{totalcoulomb} into two
2373 quickly-converging terms and a constant term:
2375 V &=& V_{\mathrm{dir}} + V_{\mathrm{rec}} + V_{0} \\[0.5ex]
2376 V_{\mathrm{dir}} &=& \frac{f}{2} \sum_{i,j}^{N}
2377 \sum_{n_x}\sum_{n_y}
2378 \sum_{n_{z}*} q_i q_j \frac{\mbox{erfc}(\beta {r}_{ij,{\bf n}} )}{{r}_{ij,{\bf n}}} \\[0.5ex]
2379 V_{\mathrm{rec}} &=& \frac{f}{2 \pi V} \sum_{i,j}^{N} q_i q_j
2380 \sum_{m_x}\sum_{m_y}
2381 \sum_{m_{z}*} \frac{\exp{\left( -(\pi {\bf m}/\beta)^2 + 2 \pi i
2382 {\bf m} \cdot ({\bf r}_i - {\bf r}_j)\right)}}{{\bf m}^2} \\[0.5ex]
2383 V_{0} &=& -\frac{f \beta}{\sqrt{\pi}}\sum_{i}^{N} q_i^2,
2385 where $\beta$ is a parameter that determines the relative weight of the
2386 direct and reciprocal sums and ${\bf m}=(m_x,m_y,m_z)$.
2387 In this way we can use a short cut-off (of the order of $1$~nm) in the direct space sum and a
2388 short cut-off in the reciprocal space sum ({\eg} 10 wave vectors in each
2389 direction). Unfortunately, the computational cost of the reciprocal
2390 part of the sum increases as $N^2$
2391 (or $N^{3/2}$ with a slightly better algorithm) and it is therefore not
2392 realistic for use in large systems.
2394 \subsubsection{Using Ewald}
2395 Don't use Ewald unless you are absolutely sure this is what you want -
2396 for almost all cases the PME method below will perform much better.
2397 If you still want to employ classical Ewald summation enter this in
2398 your {\tt .mdp} file, if the side of your box is about $3$~nm:
2405 fourierspacing = 0.6
2409 The ratio of the box dimensions and the {\tt fourierspacing} parameter determines
2410 the highest magnitude of wave vectors $m_x,m_y,m_z$ to use in each
2411 direction. With a 3-nm cubic box this example would use $11$ wave vectors
2412 (from $-5$ to $5$) in each direction. The {\tt ewald-rtol} parameter
2413 is the relative strength of the electrostatic interaction at the
2414 cut-off. Decreasing this gives you a more accurate direct sum, but a
2415 less accurate reciprocal sum.
2417 \subsection{\normindex{PME}}
2419 Particle-mesh Ewald is a method proposed by Tom
2420 Darden~\cite{Darden93} to improve the performance of the
2421 reciprocal sum. Instead of directly summing wave vectors, the charges
2422 are assigned to a grid using interpolation. The implementation in
2423 {\gromacs} uses cardinal B-spline interpolation~\cite{Essmann95},
2424 which is referred to as smooth PME (SPME).
2425 The grid is then Fourier transformed with a 3D FFT algorithm and the
2426 reciprocal energy term obtained by a single sum over the grid in
2429 The potential at the grid points is calculated by inverse
2430 transformation, and by using the interpolation factors we get the
2431 forces on each atom.
2433 The PME algorithm scales as $N \log(N)$, and is substantially faster
2434 than ordinary Ewald summation on medium to large systems. On very
2435 small systems it might still be better to use Ewald to avoid the
2436 overhead in setting up grids and transforms.
2437 For the parallelization of PME see the section on MPMD PME (\ssecref{mpmd_pme}).
2439 With the Verlet cut-off scheme, the PME direct space potential is
2440 shifted by a constant such that the potential is zero at the
2441 cut-off. This shift is small and since the net system charge is close
2442 to zero, the total shift is very small, unlike in the case of the
2443 Lennard-Jones potential where all shifts add up. We apply the shift
2444 anyhow, such that the potential is the exact integral of the force.
2446 \subsubsection{Using PME}
2447 To use Particle-mesh Ewald summation in {\gromacs}, specify the
2448 following lines in your {\tt .mdp} file:
2455 fourierspacing = 0.12
2460 In this case the {\tt fourierspacing} parameter determines the maximum
2461 spacing for the FFT grid (i.e. minimum number of grid points),
2462 and {\tt pme-order} controls the
2463 interpolation order. Using fourth-order (cubic) interpolation and this
2464 spacing should give electrostatic energies accurate to about
2465 $5\cdot10^{-3}$. Since the Lennard-Jones energies are not this
2466 accurate it might even be possible to increase this spacing slightly.
2468 Pressure scaling works with PME, but be aware of the fact that
2469 anisotropic scaling can introduce artificial ordering in some systems.
2471 \subsection{\normindex{P3M-AD}}
2473 The \seeindex{Particle-Particle Particle-Mesh}{P3M} methods of
2474 Hockney \& Eastwood can also be applied in {\gromacs} for the
2475 treatment of long range electrostatic interactions~\cite{Hockney81}.
2476 Although the P3M method was the first efficient long-range electrostatics
2477 method for molecular simulation, the smooth PME (SPME) method has largely
2478 replaced P3M as the method of choice in atomistic simulations. One performance
2479 disadvantage of the original P3M method was that it required 3 3D-FFT
2480 back transforms to obtain the forces on the particles. But this is not
2481 required for P3M and the forces can be derived through analytical differentiation
2482 of the potential, as done in PME. The resulting method is termed P3M-AD.
2483 The only remaining difference between P3M-AD and PME is the optimization
2484 of the lattice Green influence function for error minimization that P3M uses.
2485 However, in 2012 it has been shown that the SPME influence function can be
2486 modified to obtain P3M~\cite{Ballenegger2012}.
2487 This means that the advantage of error minimization in P3M-AD can be used
2488 at the same computational cost and with the same code as PME,
2489 just by adding a few lines to modify the influence function.
2490 However, at optimal parameter setting the effect of error minimization
2491 in P3M-AD is less than 10\%. P3M-AD does show large accuracy gains with
2492 interlaced (also known as staggered) grids, but that is not supported
2493 in {\gromacs} (yet).
2495 P3M is used in {\gromacs} with exactly the same options as used with PME
2496 by selecting the electrostatics type:
2498 coulombtype = P3M-AD
2501 \subsection{Optimizing Fourier transforms and PME calculations}
2502 It is recommended to optimize the parameters for calculation of
2503 electrostatic interaction such as PME grid dimensions and cut-off radii.
2504 This is particularly relevant to do before launching long production runs.
2506 {\gromacs} includes a special tool, {\tt g_tune_pme}, which automates the
2507 process of selecting the optimal size of the grid and number of PME-only
2511 % Temporarily removed since I am not sure about the state of the testlr
2514 %It is possible to test the accuracy of your settings using the program
2515 %{\tt\normindex{testlr}} in the {\tt src/gmxlib} dir. This program computes
2516 %forces and potentials using PPPM and an Ewald implementation and gives the
2517 %absolute and RMS errors in both:
2522 %Potential 0.113 0.035
2524 %{\bf Note:} these numbers were generated using a grid spacing of
2525 %0.058 nm and $r_c$ = 1.0 nm.
2527 %You can see what the accuracy is without optimizing the
2528 %$\hat{G}(k)$ function, if you pass the {\tt -ghat} option to {\tt
2529 %testlr}. Try it if you think the {\tt mk_ghat} procedure is a waste
2531 %} % Brace matches ifthenelse test for gmxlite
2534 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2535 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2536 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2538 %\ifthenelse{\equal{\gmxlite}{1}}{}{
2539 \section{Long Range Van der Waals interactions}
2540 \subsection{Dispersion correction\index{dispersion correction}}
2541 In this section, we derive long-range corrections due to the use of a
2542 cut-off for Lennard-Jones or Buckingham interactions.
2543 We assume that the cut-off is
2544 so long that the repulsion term can safely be neglected, and therefore
2545 only the dispersion term is taken into account. Due to the nature of
2546 the dispersion interaction (we are truncating a potential proportional
2547 to $-r^{-6}$), energy and pressure corrections are both negative. While
2548 the energy correction is usually small, it may be important for free
2549 energy calculations where differences between two different Hamiltonians
2550 are considered. In contrast, the pressure correction is very large and
2551 can not be neglected under any circumstances where a correct pressure is
2552 required, especially for any NPT simulations. Although it is, in
2553 principle, possible to parameterize a force field such that the pressure
2554 is close to the desired experimental value without correction, such a
2555 method makes the parameterization dependent on the cut-off and is therefore
2558 \subsubsection{Energy}
2560 The long-range contribution of the dispersion interaction to the
2561 virial can be derived analytically, if we assume a homogeneous
2562 system beyond the cut-off distance $r_c$. The dispersion energy
2563 between two particles is written as:
2565 V(\rij) ~=~- C_6\,\rij^{-6}
2567 and the corresponding force is:
2569 \Fvij ~=~- 6\,C_6\,\rij^{-8}\rvij
2571 In a periodic system it is not easy to calculate the full potentials,
2572 so usually a cut-off is applied, which can be abrupt or smooth.
2573 We will call the potential and force with cut-off $V_c$ and $\ve{F}_c$.
2574 The long-range contribution to the dispersion energy
2575 in a system with $N$ particles and particle density $\rho$ = $N/V$ is:
2577 \label{eqn:enercorr}
2578 V_{lr} ~=~ \half N \rho\int_0^{\infty} 4\pi r^2 g(r) \left( V(r) -V_c(r) \right) {\dr}
2580 We will integrate this for the shift function, which is the most general
2581 form of van der Waals interaction available in {\gromacs}.
2582 The shift function has a constant difference $S$ from 0 to $r_1$
2583 and is 0 beyond the cut-off distance $r_c$.
2584 We can integrate \eqnref{enercorr}, assuming that the density in the sphere
2585 within $r_1$ is equal to the global density and
2586 the radial distribution function $g(r)$ is 1 beyond $r_1$:
2589 V_{lr} &=& \half N \left(
2590 \rho\int_0^{r_1} 4\pi r^2 g(r) \, C_6 \,S\,{\dr}
2591 + \rho\int_{r_1}^{r_c} 4\pi r^2 \left( V(r) -V_c(r) \right) {\dr}
2592 + \rho\int_{r_c}^{\infty} 4\pi r^2 V(r) \, {\dr}
2594 & = & \half N \left(\left(\frac{4}{3}\pi \rho r_1^{3} - 1\right) C_6 \,S
2595 + \rho\int_{r_1}^{r_c} 4\pi r^2 \left( V(r) -V_c(r) \right) {\dr}
2596 -\frac{4}{3} \pi N \rho\, C_6\,r_c^{-3}
2599 where the term $-1$ corrects for the self-interaction.
2600 For a plain cut-off we only need to assume that $g(r)$ is 1 beyond $r_c$
2601 and the correction reduces to~\cite{Allen87}:
2603 V_{lr} & = & -\frac{2}{3} \pi N \rho\, C_6\,r_c^{-3}
2605 If we consider, for example, a box of pure water, simulated with a cut-off
2606 of 0.9 nm and a density of 1 g cm$^{-3}$ this correction is
2607 $-0.75$ kJ mol$^{-1}$ per molecule.
2609 For a homogeneous mixture we need to define
2610 an {\em average dispersion constant}:
2613 \avcsix = \frac{2}{N(N-1)}\sum_i^N\sum_{j>i}^N C_6(i,j)\\
2615 In {\gromacs}, excluded pairs of atoms do not contribute to the average.
2617 In the case of inhomogeneous simulation systems, {\eg} a system with a
2618 lipid interface, the energy correction can be applied if
2619 $\avcsix$ for both components is comparable.
2621 \subsubsection{Virial and pressure}
2622 The scalar virial of the system due to the dispersion interaction between
2623 two particles $i$ and $j$ is given by:
2625 \Xi~=~-\half \rvij \cdot \Fvij ~=~ 3\,C_6\,\rij^{-6}
2627 The pressure is given by:
2629 P~=~\frac{2}{3\,V}\left(E_{kin} - \Xi\right)
2631 The long-range correction to the virial is given by:
2633 \Xi_{lr} ~=~ \half N \rho \int_0^{\infty} 4\pi r^2 g(r) (\Xi -\Xi_c) \,\dr
2635 We can again integrate the long-range contribution to the
2636 virial assuming $g(r)$ is 1 beyond $r_1$:
2638 \Xi_{lr}&=& \half N \rho \left(
2639 \int_{r_1}^{r_c} 4 \pi r^2 (\Xi -\Xi_c) \,\dr
2640 + \int_{r_c}^{\infty} 4 \pi r^2 3\,C_6\,\rij^{-6}\, \dr
2642 &=& \half N \rho \left(
2643 \int_{r_1}^{r_c} 4 \pi r^2 (\Xi -\Xi_c) \, \dr
2644 + 4 \pi C_6 \, r_c^{-3} \right)
2646 For a plain cut-off the correction to the pressure is~\cite{Allen87}:
2648 P_{lr}~=~-\frac{4}{3} \pi C_6\, \rho^2 r_c^{-3}
2650 Using the same example of a water box, the correction to the virial is
2651 0.75 kJ mol$^{-1}$ per molecule,
2652 the corresponding correction to the pressure for
2653 SPC water is approximately $-280$ bar.
2655 For homogeneous mixtures, we can again use the average dispersion constant
2656 $\avcsix$ (\eqnref{avcsix}):
2658 P_{lr}~=~-\frac{4}{3} \pi \avcsix \rho^2 r_c^{-3}
2661 For inhomogeneous systems, \eqnref{pcorr} can be applied under the same
2662 restriction as holds for the energy (see \secref{ecorr}).
2664 \subsection{Lennard-Jones PME\index{LJ-PME}}
2666 In order to treat systems, using Lennard-Jones potentials, that are
2667 non-homogeneous outside of the cut-off distance, we can instead use
2668 the Particle-mesh Ewald method as discussed for electrostatics above.
2669 In this case the modified Ewald equations become
2671 V &=& V_{\mathrm{dir}} + V_{\mathrm{rec}} + V_{0} \\[0.5ex]
2672 V_{\mathrm{dir}} &=& -\frac{1}{2} \sum_{i,j}^{N}
2673 \sum_{n_x}\sum_{n_y}
2674 \sum_{n_{z}*} \frac{C_{ij}^{(6)}g(\beta {r}_{ij,{\bf n}})}{{r_{ij,{\bf n}}}^6} \\[0.5ex]
2675 V_{\mathrm{rec}} &=& \frac{{\pi}^{\frac{3}{2}} \beta^{3}}{2V} \sum_{m_x}\sum_{m_y}\sum_{m_{z}*}
2676 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]
2677 V_{0} &=& -\frac{\beta^{6}}{12}\sum_{i}^{N} C_{ii}^{(6)},
2680 where ${\bf m}=(m_x,m_y,m_z)$, $\beta$ is the parameter determining the weight between
2681 direct and reciprocal space, and ${C_{ij}}^{(6)}$ is the combined dispersion
2682 parameter for particle $i$ and $j$. The star indicates that terms
2683 with $i = j$ should be omitted when $((n_x,n_y,n_z)=(0,0,0))$, and
2684 ${\bf r}_{ij,{\bf n}}$ is the real distance between the particles.
2685 Following the derivation by Essmann~\cite{Essmann95}, the functions $f$ and $g$ introduced above are defined as
2687 f(x)&=&1/3\left[(1-2x^2){\mathrm{exp}}(-x^2) + 2{x^3}\sqrt{\pi}\,{\mathrm{erfc}}(x) \right] \\
2688 g(x)&=&{\mathrm{exp}}(-x^2)(1+x^2+\frac{x^4}{2}).
2691 The above methodology works fine as long as the dispersion parameters can be factorized in the same
2692 way as the charges for electrostatics
2694 C_{ij}^{(6)} = \left({C_{ii}^{(6)} \, C_{jj}^{(6)}}\right)^{1/2}
2696 For Lorentz-Berthelot combination rules, the reciprocal part of this sum has to be calculated
2697 seven times due to the splitting of the dispersion parameter according to
2699 C_{ij}^{(6)}=(\sigma_i+\sigma_j)^6=\sum_{n=0}^{6} P_{n}\sigma_{i}^{n}\sigma_{j}^{(6-n)},
2701 for $P_{n}$ the Pascal triangle coefficients. This introduces a
2702 non-negligible cost to the reciprocal part, requiring seven separate
2703 FFTs, and therefore this have been the limiting factor in previous
2704 attempts to implement LJ-PME. A solution to this problem is to
2705 approximate the interaction parameters in reciprocal space using
2706 geometrical combination rules. This will preserve a well-defined
2707 Hamiltonian and significantly increase the performance of the
2711 \centerline{\includegraphics[width=15cm]{plots/ljpmedifference}}
2712 \caption {Dispersion potential between phosphorous and oxygen, The total, real and
2713 reciprocal parts of the total dispersion are shown. The reciprocal parts are calculated
2714 using either LB or geometric rules. The difference introduced by the use of the
2715 geometric approximation in the reciprocal part is small compared to the total
2716 interaction energy.}
2717 \label{fig:ljpmedifference}
2720 This approximation does introduce some errors, but since the
2721 difference is located in the interactions calculated in reciprocal
2722 space, the effect will be very small compared to the total interaction
2723 energy (see \figref{ljpmedifference}). The relative error in
2724 the total dispersion energy will stay below 0.5\% in a lipid bilayer
2725 simulation, when using a real space cut-off of 1.0 nm. A more
2726 thorough discussion of this can be found in \cite{Wennberg13}.
2728 \subsubsection{Using LJ-PME}
2729 To use Particle-mesh Ewald summation for Lennard-Jones interactions in {\gromacs}, specify the
2730 following lines in your {\tt .mdp} file:
2736 fourierspacing = 0.12
2738 ewald-rtol-lj = 0.001
2739 lj-pme-comb-rule = geometric
2742 The {\tt fourierspacing} and {\tt pme-order} are the same parameters
2743 as is used for electrostatic PME, and {\tt ewald-rtol-lj} controls
2744 splitting between real and reciprocal space in the same way as
2745 {\tt ewald-rtol}. In addition to this, the combination rule to be used
2746 in reciprocal space is determined by {\tt lj-pme-comb-rule}. If the
2747 current force field uses Lorentz-Berthelot combination rules, it is
2748 possible to set {\tt lj-pme-comb-rule = geometric} in order to gain a
2749 significant increase in performance for a small loss in accuracy. The
2750 details of this approximation can be found in the section above. The
2751 implementation of LJ-PME is currently unsupported together with the
2752 Verlet cut-off scheme and/or free energy calculations. These features
2753 will be added in upcoming releases
2754 %} % Brace matches ifthenelse test for gmxlite
2756 \section{Force field\index{force field}}
2758 A force field is built up from two distinct components:
2760 \item The set of equations (called the {\em
2761 \index{potential function}s}) used to generate the potential
2762 energies and their derivatives, the forces. These are described in
2763 detail in the previous chapter.
2764 \item The parameters used in this set of equations. These are not
2765 given in this manual, but in the data files corresponding to your
2766 {\gromacs} distribution.
2768 Within one set of equations various sets of parameters can be
2769 used. Care must be taken that the combination of equations and
2770 parameters form a consistent set. It is in general dangerous to make
2771 {\em ad hoc} changes in a subset of parameters, because the various
2772 contributions to the total force are usually interdependent. This
2773 means in principle that every change should be documented, verified by
2774 comparison to experimental data and published in a peer-reviewed
2775 journal before it can be used.
2777 {\gromacs} {\gmxver} includes several force fields, and additional
2778 ones are available on the website. If you do not know which one to
2779 select we recommend \gromosv{96} for united-atom setups and OPLS-AA/L for
2780 all-atom parameters. That said, we describe the available options in
2783 \subsubsection{All-hydrogen force field}
2784 The \gromosv{87}-based all-hydrogen force field is almost identical to the
2785 normal \gromosv{87} force field, since the extra hydrogens have no
2786 Lennard-Jones interaction and zero charge. The only differences are in
2787 the bond angle and improper dihedral angle terms. This force field is
2788 only useful when you need the exact hydrogen positions, for instance
2789 for distance restraints derived from NMR measurements. When citing
2790 this force field please read the previous paragraph.
2792 \subsection{\gromosv{96}\index{GROMOS96 force field}}
2793 {\gromacs} supports the \gromosv{96} force fields~\cite{gromos96}.
2794 All parameters for the 43A1, 43A2 (development, improved alkane
2795 dihedrals), 45A3, 53A5, and 53A6 parameter sets are included. All standard
2796 building blocks are included and topologies can be built automatically
2799 The \gromosv{96} force field is a further development of the \gromosv{87} force field.
2800 It has improvements over the \gromosv{87} force field for proteins and small molecules.
2801 {\bf Note} that the sugar parameters present in 53A6 do correspond to those published in
2802 2004\cite{Oostenbrink2004}, which are different from those present in 45A4, which
2803 is not included in {\gromacs} at this time. The 45A4 parameter set corresponds to a later
2804 revision of these parameters.
2805 The \gromosv{96} force field is not, however, recommended for use with long alkanes and
2806 lipids. The \gromosv{96} force field differs from the \gromosv{87}
2807 force field in a few respects:
2809 \item the force field parameters
2810 \item the parameters for the bonded interactions are not linked to atom types
2811 \item a fourth power bond stretching potential (\ssecref{G96bond})
2812 \item an angle potential based on the cosine of the angle (\ssecref{G96angle})
2814 There are two differences in implementation between {\gromacs} and \gromosv{96}
2815 which can lead to slightly different results when simulating the same system
2818 \item in \gromosv{96} neighbor searching for solvents is performed on the
2819 first atom of the solvent molecule. This is not implemented in {\gromacs},
2820 but the difference with searching by centers of charge groups is very small
2821 \item the virial in \gromosv{96} is molecule-based. This is not implemented in
2822 {\gromacs}, which uses atomic virials
2824 The \gromosv{96} force field was parameterized with a Lennard-Jones cut-off
2825 of 1.4 nm, so be sure to use a Lennard-Jones cut-off ({\tt rvdw}) of at least 1.4.
2826 A larger cut-off is possible because the Lennard-Jones potential and forces
2827 are almost zero beyond 1.4 nm.
2829 \subsubsection{\gromosv{96} files\swapindexquiet{GROMOS96}{files}}
2830 {\gromacs} can read and write \gromosv{96} coordinate and trajectory files.
2831 These files should have the extension {\tt .g96}.
2832 Such a file can be a \gromosv{96} initial/final
2833 configuration file, a coordinate trajectory file, or a combination of both.
2834 The file is fixed format; all floats are written as 15.9, and as such, files can get huge.
2835 {\gromacs} supports the following data blocks in the given order:
2845 POSITION/POSITIONRED (mandatory)
2846 VELOCITY/VELOCITYRED (optional)
2851 See the \gromosv{96} manual~\cite{gromos96} for a complete description
2852 of the blocks. {\bf Note} that all {\gromacs} programs can read compressed
2853 (.Z) or gzipped (.gz) files.
2855 \subsection{OPLS/AA\index{OPLS/AA force field}}
2857 \subsection{AMBER\index{AMBER force field}}
2859 {\gromacs} provides native support for the following AMBER force fields:
2862 \item AMBER94~\cite{Cornell1995}
2863 \item AMBER96~\cite{Kollman1996}
2864 \item AMBER99~\cite{Wang2000}
2865 \item AMBER99SB~\cite{Hornak2006}
2866 \item AMBER99SB-ILDN~\cite{Lindorff2010}
2867 \item AMBER03~\cite{Duan2003}
2868 \item AMBERGS~\cite{Garcia2002}
2871 \subsection{CHARMM\index{CHARMM force field}}
2872 \label{subsec:charmmff}
2874 {\gromacs} supports the CHARMM force field for proteins~\cite{mackerell04, mackerell98}, lipids~\cite{feller00} and nucleic acids~\cite{foloppe00,Mac2000}. The protein parameters (and to some extent the lipid and nucleic acid parameters) were thoroughly tested -- both by comparing potential energies between the port and the standard parameter set in the CHARMM molecular simulation package, as well by how the protein force field behaves together with {\gromacs}-specific techniques such as virtual sites (enabling long time steps) and a fast implicit solvent recently implemented~\cite{Larsson10} -- and the details and results are presented in the paper by Bjelkmar et al.~\cite{Bjelkmar10}. The nucleic acid parameters, as well as the ones for HEME, were converted and tested by Michel Cuendet.
2876 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.
2878 A port of the CHARMM36 force field for use with GROMACS is also available at \url{http://mackerell.umaryland.edu/charmm_ff.shtml#gromacs}.
2880 \subsection{Coarse-grained force fields}
2881 \index{force-field, coarse-grained}
2882 \label{sec:cg-forcefields}
2883 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.
2885 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:
2887 \item Conserving free energies
2889 \item Simplex method
2890 \item MARTINI force field (see next section)
2892 \item Conserving distributions (like the radial distribution function), so-called structure-based coarse-graining
2894 \item (iterative) Boltzmann inversion
2895 \item Inverse Monte Carlo
2897 \item Conversing forces
2899 \item Force matching
2903 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}.
2905 \subsection{MARTINI\index{Martini force field}}
2907 The MARTINI force field is a coarse-grain parameter set that allows for the construction
2908 of many systems, including proteins and membranes.
2910 \subsection{PLUM\index{PLUM force field}}
2912 The \normindex{PLUM force field}~\cite{bereau12} is an example of a solvent-free
2913 protein-membrane model for which the membrane was derived from structure-based
2914 coarse-graining~\cite{wang_jpcb10}. A {\gromacs} implementation can be found at
2915 \href{http://code.google.com/p/plumx/}{code.google.com/p/plumx}.
2917 % LocalWords: dihedrals centro ij dV dr LJ lj rcl jj Bertelot OPLS bh bham rf
2918 % LocalWords: coul defunits grompp crf vcrf fcrf Tironi Debye kgrf cgrf krf dx
2919 % LocalWords: PPPM der Waals erfc lr elstat chirality bstretch bondpot kT kJ
2920 % LocalWords: anharmonic morse mol betaij expminx SPC timestep fs FENE ijk kj
2921 % LocalWords: anglepot CHARMm UB ik rr substituents ijkl Ryckaert Bellemans rb
2922 % LocalWords: alkanes pdb gmx IUPAC IUB jkl cis rbdih crb kcal cubicspline xvg
2923 % LocalWords: topfile mdrun posres ar dihr lcllll NMR nmr lcllllll NOEs lclll
2924 % LocalWords: rav preprocessor ccccccccc ai aj fac disre mdp multi topol tpr
2925 % LocalWords: fc ravdisre nstdisreout dipolar lll ccc orientational MSD const
2926 % LocalWords: orire fitgrp nstorireout Drude intra Noskov et al fecalc coulrf
2927 % LocalWords: polarizabilities parameterized sigeps Ek sc softcore GROMOS NBF
2928 % LocalWords: hydrogens alkane vdwtype coulombtype mdpopt rlist rcoulomb rvdw
2929 % LocalWords: nstlist virial funcparm VdW jk jl fvsite fd vsites lcr vsitetop
2930 % LocalWords: vsite lclllll lcl parameterize parameterization enercorr avcsix
2931 % LocalWords: pcorr ecorr totalcoulomb dir fourierspacing ewald rtol Darden gz
2932 % LocalWords: FFT parallelization MPMD mpmd pme fft hoc Gromos gromos oxygens
2933 % LocalWords: virials POSITIONRED VELOCITYRED gzipped Charmm Larsson Bjelkmar
2934 % LocalWords: Cuendet CMAP nocmap dihedral Lennard covalent Verlet
2935 % LocalWords: Berthelot nm flexwat ferguson itp harmonicangle versa
2936 % LocalWords: harmonicbond atomtypes dihedraltypes equilibrated fdn
2937 % LocalWords: distancerestraint LINCS Coulombic ja jb il SPME ILDN
2938 % LocalWords: Hamiltonians atomtype AMBERGS rtp cmap graining VOTCA