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 \newcommand{\nproc}{\mbox{$M$
}}
36 \newcommand{\natom}{\mbox{$N$
}}
37 \newcommand{\nx}{\mbox{$n_x$
}}
38 \newcommand{\ny}{\mbox{$n_y$
}}
39 \newcommand{\nz}{\mbox{$n_z$
}}
40 \newcommand{\nsgrid}{NS grid
}
41 \newcommand{\fftgrid}{FFT grid
}
42 \newcommand{\dgrid}{\mbox{$
\delta_{grid
}$
}}
43 \newcommand{\bfv}[1]{{\mbox{\boldmath{$
#1$
}}}}
44 % non-italicized boldface for math (e.g. matrices)
45 \newcommand{\bfm}[1]{{\bf #1}}
46 \newcommand{\dt}{\Delta t
}
47 \newcommand{\rv}{\bfv{r
}}
48 \newcommand{\vv}{\bfv{v
}}
49 \newcommand{\F}{\bfv{F
}}
50 \newcommand{\pb}{\bfv{p
}}
51 \newcommand{\veps}{v_
{\epsilon}}
52 \newcommand{\peps}{p_
{\epsilon}}
53 \newcommand{\sinhx}[1]{\frac{\sinh{\left(
#1\right)
}}{#1}}
56 \section{Introduction
}
57 In this chapter we first give describe some general concepts used in
58 {\gromacs}:
{\em periodic boundary conditions
} (
\secref{pbc
})
59 and the
{\em group concept
} (
\secref{groupconcept
}). The MD algorithm is
60 described in
\secref{MD
}: first a global form of the algorithm is
61 given, which is refined in subsequent subsections. The (simple) EM
62 (Energy Minimization) algorithm is described in
\secref{EM
}. Some
63 other algorithms for special purpose dynamics are described after
66 %\ifthenelse{\equal{\gmxlite}{1}}{}{
67 %In the final \secref{par} of this chapter a few principles are
68 %given on which parallelization of {\gromacs} is based. The
69 %parallelization is hardly visible for the user and is therefore not
71 %} % Brace matches ifthenelse test for gmxlite
73 A few issues are of general interest. In all cases the
{\em system
}
74 must be defined, consisting of molecules. Molecules again consist of
75 particles with defined interaction functions. The detailed
76 description of the
{\em topology
} of the molecules and of the
{\em force
77 field
} and the calculation of forces is given in
78 \chref{ff
}. In the present chapter we describe
79 other aspects of the algorithm, such as pair list generation, update of
80 velocities and positions, coupling to external temperature and
81 pressure, conservation of constraints.
82 \ifthenelse{\equal{\gmxlite}{1}}{}{
83 The
{\em analysis
} of the data generated by an MD simulation is treated in
\chref{analysis
}.
84 } % Brace matches ifthenelse test for gmxlite
86 \section{Periodic boundary conditions
\index{periodic boundary conditions
}}
89 \centerline{\includegraphics[width=
9cm
]{plots/pbctric
}}
90 \caption {Periodic boundary conditions in two dimensions.
}
93 The classical way to minimize edge effects in a finite system is to
94 apply
{\em periodic boundary conditions
}. The atoms of the system to
95 be simulated are put into a space-filling box, which is surrounded by
96 translated copies of itself (
\figref{pbc
}). Thus there are no
97 boundaries of the system; the artifact caused by unwanted boundaries
98 in an isolated cluster is now replaced by the artifact of periodic
99 conditions. If the system is crystalline, such boundary conditions are
100 desired (although motions are naturally restricted to periodic motions
101 with wavelengths fitting into the box). If one wishes to simulate
102 non-periodic systems, such as liquids or solutions, the periodicity by
103 itself causes errors. The errors can be evaluated by comparing various
104 system sizes; they are expected to be less severe than the errors
105 resulting from an unnatural boundary with vacuum.
107 There are several possible shapes for space-filling unit cells. Some,
108 like the
{\em \normindex{rhombic dodecahedron
}} and the
109 {\em \normindex{truncated octahedron
}}~
\cite{Adams79
} are closer to being a sphere
110 than a cube is, and are therefore better suited to the
111 study of an approximately spherical macromolecule in solution, since
112 fewer solvent molecules are required to fill the box given a minimum
113 distance between macromolecular images. At the same time, rhombic
114 dodecahedra and truncated octahedra are special cases of
{\em triclinic
}
115 unit cells
\index{triclinic unit cell
}; the most general space-filling unit cells
116 that comprise all possible space-filling shapes~
\cite{Bekker95
}.
117 For this reason,
{\gromacs} is based on the triclinic unit cell.
119 {\gromacs} uses periodic boundary conditions, combined with the
{\em
120 \normindex{minimum image convention
}}: only one -- the nearest -- image of each
121 particle is considered for short-range non-bonded interaction terms.
122 For long-range electrostatic interactions this is not always accurate
123 enough, and
{\gromacs} therefore also incorporates lattice sum methods
124 such as Ewald Sum, PME and PPPM.
126 {\gromacs} supports triclinic boxes of any shape.
127 The simulation box (unit cell) is defined by the
3 box vectors
128 $
{\bf a
}$,$
{\bf b
}$ and $
{\bf c
}$.
129 The box vectors must satisfy the following conditions:
135 \label{eqn:box_shift1
}
136 a_x>
0,~~~~b_y>
0,~~~~c_z>
0
139 \label{eqn:box_shift2
}
140 |b_x|
\leq \frac{1}{2} \, a_x,~~~~
141 |c_x|
\leq \frac{1}{2} \, a_x,~~~~
142 |c_y|
\leq \frac{1}{2} \, b_y
144 Equations
\ref{eqn:box_rot
} can always be satisfied by rotating the box.
145 Inequalities (
\ref{eqn:box_shift1
}) and (
\ref{eqn:box_shift2
}) can always be
146 satisfied by adding and subtracting box vectors.
148 Even when simulating using a triclinic box,
{\gromacs} always keeps the
149 particles in a brick-shaped volume for efficiency,
150 as illustrated in
\figref{pbc
} for a
2-dimensional system.
151 Therefore, from the output trajectory it might seem that the simulation was
152 done in a rectangular box. The program
{\tt trjconv
} can be used to convert
153 the trajectory to a different unit-cell representation.
155 It is also possible to simulate without periodic boundary conditions,
156 but it is usually more efficient to simulate an isolated cluster of molecules
157 in a large periodic box, since fast grid searching can only be used
158 in a periodic system.
162 \includegraphics[width=
5cm
]{plots/rhododec
}
163 ~~~~
\includegraphics[width=
5cm
]{plots/truncoct
}
165 \caption {A rhombic dodecahedron and truncated octahedron
166 (arbitrary orientations).
}
167 \label{fig:boxshapes
}
170 \subsection{Some useful box types
}
173 \begin{tabular
}{|c|c|c|ccc|ccc|
}
175 box type & image & box &
\multicolumn{3}{c|
}{box vectors
} &
\multicolumn{3}{c|
}{box vector angles
} \\
176 & distance & volume & ~
{\bf a
}~ &
{\bf b
} &
{\bf c
} &
177 $
\angle{\bf bc
}$ & $
\angle{\bf ac
}$ & $
\angle{\bf ab
}$ \\
179 & & & $d$ &
0 &
0 & & & \\
180 cubic & $d$ & $d^
3$ &
0 & $d$ &
0 & $
90^
\circ$ & $
90^
\circ$ & $
90^
\circ$ \\
181 & & &
0 &
0 & $d$ & & & \\
183 rhombic & & & $d$ &
0 & $
\frac{1}{2}\,d$ & & & \\
184 dodecahedron & $d$ & $
\frac{1}{2}\sqrt{2}\,d^
3$ &
0 & $d$ & $
\frac{1}{2}\,d$ & $
60^
\circ$ & $
60^
\circ$ & $
90^
\circ$ \\
185 (xy-square) & & $
0.707\,d^
3$ &
0 &
0 & $
\frac{1}{2}\sqrt{2}\,d$ & & & \\
187 rhombic & & & $d$ & $
\frac{1}{2}\,d$ & $
\frac{1}{2}\,d$ & & & \\
188 dodecahedron & $d$ & $
\frac{1}{2}\sqrt{2}\,d^
3$ &
0 & $
\frac{1}{2}\sqrt{3}\,d$ & $
\frac{1}{6}\sqrt{3}\,d$ & $
60^
\circ$ & $
60^
\circ$ & $
60^
\circ$ \\
189 (xy-hexagon) & & $
0.707\,d^
3$ &
0 &
0 & $
\frac{1}{3}\sqrt{6}\,d$ & & & \\
191 truncated & & & $d$ & $
\frac{1}{3}\,d$ & $-
\frac{1}{3}\,d$ & & &\\
192 octahedron & $d$ & $
\frac{4}{9}\sqrt{3}\,d^
3$ &
0 & $
\frac{2}{3}\sqrt{2}\,d$ & $
\frac{1}{3}\sqrt{2}\,d$ & $
71.53^
\circ$ & $
109.47^
\circ$ & $
71.53^
\circ$ \\
193 & & $
0.770\,d^
3$ &
0 &
0 & $
\frac{1}{3}\sqrt{6}\,d$ & & & \\
197 \caption{The cubic box, the rhombic
\normindex{dodecahedron
} and the truncated
198 \normindex{octahedron
}.
}
201 The three most useful box types for simulations of solvated systems
202 are described in
\tabref{boxtypes
}. The rhombic dodecahedron
203 (
\figref{boxshapes
}) is the smallest and most regular space-filling
204 unit cell. Each of the
12 image cells is at the same distance. The
205 volume is
71\% of the volume of a cube having the same image
206 distance. This saves about
29\% of CPU-time when simulating a
207 spherical or flexible molecule in solvent. There are two different
208 orientations of a rhombic dodecahedron that satisfy equations
209 \ref{eqn:box_rot
},
\ref{eqn:box_shift1
} and
\ref{eqn:box_shift2
}.
210 The program
{\tt editconf
} produces the orientation
211 which has a square intersection with the xy-plane. This orientation
212 was chosen because the first two box vectors coincide with the x and
213 y-axis, which is easier to comprehend. The other orientation can be
214 useful for simulations of membrane proteins. In this case the
215 cross-section with the xy-plane is a hexagon, which has an area which
216 is
14\% smaller than the area of a square with the same image
217 distance. The height of the box ($c_z$) should be changed to obtain
218 an optimal spacing. This box shape not only saves CPU time, it
219 also results in a more uniform arrangement of the proteins.
221 \subsection{Cut-off restrictions
}
222 The
\normindex{minimum image convention
} implies that the cut-off radius used to
223 truncate non-bonded interactions may not exceed half the shortest box
226 \label{eqn:physicalrc
}
227 R_c <
\half \min(\|
{\bf a
}\|,\|
{\bf b
}\|,\|
{\bf c
}\|),
229 because otherwise more than one image would be within the cut-off distance
230 of the force. When a macromolecule, such as a protein, is studied in
231 solution, this restriction alone is not sufficient: in principle, a single
232 solvent molecule should not be able
233 to `see' both sides of the macromolecule. This means that the length of
234 each box vector must exceed the length of the macromolecule in the
235 direction of that edge
{\em plus
} two times the cut-off radius $R_c$.
236 It is, however, common to compromise in this respect, and make the solvent
237 layer somewhat smaller in order to reduce the computational cost.
238 For efficiency reasons the cut-off with triclinic boxes is more restricted.
239 For grid search the extra restriction is weak:
242 R_c <
\min(a_x,b_y,c_z)
244 For simple search the extra restriction is stronger:
247 R_c <
\half \min(a_x,b_y,c_z)
250 Each unit cell (cubic, rectangular or triclinic)
251 is surrounded by
26 translated images. A
252 particular image can therefore always be identified by an index pointing to one
253 of
27 {\em translation vectors
} and constructed by applying a
254 translation with the indexed vector (see
\ssecref{forces
}).
255 Restriction (
\ref{eqn:gridrc
}) ensures that only
26 images need to be
258 %\ifthenelse{\equal{\gmxlite}{1}}{}{
259 \section{The group concept
}
260 \label{sec:groupconcept
}\index{group
}
261 The
{\gromacs} MD and analysis programs use user-defined
{\em groups
} of
262 atoms to perform certain actions on. The maximum number of groups is
263 256, but each atom can only belong to six different groups, one
264 each of the following:
266 \item[temperature-coupling group
\swapindex{temperature-coupling
}{group
}]
267 The
\normindex{temperature coupling
} parameters (reference
268 temperature, time constant, number of degrees of freedom, see
269 \ssecref{update
}) can be defined for each T-coupling group
270 separately. For example, in a solvated macromolecule the solvent (that
271 tends to generate more heating by force and integration errors) can be
272 coupled with a shorter time constant to a bath than is a macromolecule,
273 or a surface can be kept cooler than an adsorbing molecule. Many
274 different T-coupling groups may be defined. See also center of mass
277 \item[\swapindex{freeze
}{group
}\index{frozen atoms
}]
278 Atoms that belong to a freeze group are kept stationary in the
279 dynamics. This is useful during equilibration,
{\eg} to avoid badly
280 placed solvent molecules giving unreasonable kicks to protein atoms,
281 although the same effect can also be obtained by putting a restraining
282 potential on the atoms that must be protected. The freeze option can
283 be used, if desired, on just one or two coordinates of an atom,
284 thereby freezing the atoms in a plane or on a line. When an atom is
285 partially frozen, constraints will still be able to move it, even in a
286 frozen direction. A fully frozen atom can not be moved by constraints.
287 Many freeze groups can be defined. Frozen coordinates are unaffected
288 by pressure scaling; in some cases this can produce unwanted results,
289 particularly when constraints are also used (in this case you will
290 get very large pressures). Accordingly, it is recommended to avoid
291 combining freeze groups with constraints and pressure coupling. For the
292 sake of equilibration it could suffice to start with freezing in a
293 constant volume simulation, and afterward use position restraints in
294 conjunction with constant pressure.
296 \item[\swapindex{accelerate
}{group
}]
297 On each atom in an ``accelerate group'' an acceleration
298 $
\ve{a
}^g$ is imposed. This is equivalent to an external
299 force. This feature makes it possible to drive the system into a
300 non-equilibrium state and enables the performance of
301 \swapindex{non-equilibrium
}{MD
} and hence to obtain transport properties.
303 \item[\swapindex{energy-monitor
}{group
}]
304 Mutual interactions between all energy-monitor groups are compiled
305 during the simulation. This is done separately for Lennard-Jones and
306 Coulomb terms. In principle up to
256 groups could be defined, but
307 that would lead to
256$
\times$
256 items! Better use this concept
310 All non-bonded interactions between pairs of energy-monitor groups can
311 be excluded
\index{exclusions
}
312 \ifthenelse{\equal{\gmxlite}{1}}
314 {(see
\secref{mdpopt
}).
}
315 Pairs of particles from excluded pairs of energy-monitor groups
316 are not put into the pair list.
317 This can result in a significant speedup
318 for simulations where interactions within or between parts of the system
321 \item[\swapindex{center of mass
}{group
}\index{removing COM motion
}]
322 In
\gromacs\ the center of mass (COM) motion can be removed, for
323 either the complete system or for groups of atoms. The latter is
324 useful,
{\eg} for systems where there is limited friction (
{\eg} gas
325 systems) to prevent center of mass motion to occur. It makes sense to
326 use the same groups for temperature coupling and center of mass motion
329 \item[\swapindex{Compressed position output
}{group
}]
331 In order to further reduce the size of the compressed trajectory file
332 (
{\tt .xtc
{\index{XTC
}}} or
{\tt .tng
{\index{TNG
}}}), it is possible
333 to store only a subset of all particles. All x-compression groups that
334 are specified are saved, the rest are not. If no such groups are
335 specified, than all atoms are saved to the compressed trajectory file.
338 The use of groups in
{\gromacs} tools is described in
339 \secref{usinggroups
}.
340 %} % Brace matches ifthenelse test for gmxlite
342 \section{Molecular Dynamics
}
346 \addtolength{\fboxsep}{0.5cm
}
347 \begin{shadowenv
}[12cm
]
348 {\large \bf THE GLOBAL MD ALGORITHM
}
349 \rule{\textwidth}{2pt
} \\
350 {\bf 1. Input initial conditions
}\\
[2ex
]
351 Potential interaction $V$ as a function of atom positions\\
352 Positions $
\ve{r
}$ of all atoms in the system\\
353 Velocities $
\ve{v
}$ of all atoms in the system \\
355 \rule{\textwidth}{1pt
}\\
356 {\bf repeat
2,
3,
4} for the required number of steps:\\
357 \rule{\textwidth}{1pt
}\\
358 {\bf 2. Compute forces
} \\
[1ex
]
359 The force on any atom \\
[1ex
]
360 $
\ve{F
}_i = -
\displaystyle\frac{\partial V
}{\partial \ve{r
}_i
}$ \\
[1ex
]
361 is computed by calculating the force between non-bonded atom pairs: \\
362 $
\ve{F
}_i =
\sum_j \ve{F
}_
{ij
}$ \\
363 plus the forces due to bonded interactions (which may depend on
1,
2,
364 3, or
4 atoms), plus restraining and/or external forces. \\
365 The potential and kinetic energies and the pressure tensor are computed. \\
367 {\bf 3. Update configuration
} \\
[1ex
]
368 The movement of the atoms is simulated by numerically solving Newton's
369 equations of motion \\
[1ex
]
371 \frac {\de^
2\ve{r
}_i
}{\de t^
2} =
\frac{\ve{F
}_i
}{m_i
} $ \\
374 \frac{\de\ve{r
}_i
}{\de t
} =
\ve{v
}_i ; \;\;
375 \frac{\de\ve{v
}_i
}{\de t
} =
\frac{\ve{F
}_i
}{m_i
} $ \\
[1ex
]
377 {\bf 4.
} if required:
{\bf Output step
} \\
378 write positions, velocities, energies, temperature, pressure, etc. \\
380 \caption{The global MD algorithm
}
384 A global flow scheme for MD is given in
\figref{global
}. Each
385 MD or EM run requires as input a set of initial coordinates and --
386 optionally -- initial velocities of all particles involved. This
387 chapter does not describe how these are obtained; for the setup of an
388 actual MD run check the online manual at
{\wwwpage}.
390 \subsection{Initial conditions
}
391 \subsubsection{Topology and force field
}
392 The system topology, including a description of the force field, must
394 \ifthenelse{\equal{\gmxlite}{1}}
396 {Force fields and topologies are described in
\chref{ff
}
397 and
\ref{ch:top
}, respectively.
}
398 All this information is static; it is never modified during the run.
400 \subsubsection{Coordinates and velocities
}
402 \centerline{\includegraphics[width=
8cm
]{plots/maxwell
}}
403 \caption{A Maxwell-Boltzmann velocity distribution, generated from
408 Then, before a run starts, the box size and the coordinates and
409 velocities of all particles are required. The box size and shape is
410 determined by three vectors (nine numbers) $
\ve{b
}_1,
\ve{b
}_2,
\ve{b
}_3$,
411 which represent the three basis vectors of the periodic box.
413 If the run starts at $t=t_0$, the coordinates at $t=t_0$ must be
414 known. The
{\em leap-frog algorithm
}, the default algorithm used to
415 update the time step with $
\Dt$ (see
\ssecref{update
}), also requires
416 that the velocities at $t=t_0 -
\hDt$ are known. If velocities are not
417 available, the program can generate initial atomic velocities
418 $v_i, i=
1\ldots 3N$ with a
\index{Maxwell-Boltzmann distribution
}
419 (
\figref{maxwell
}) at a given absolute temperature $T$:
421 p(v_i) =
\sqrt{\frac{m_i
}{2 \pi kT
}}\exp\left(-
\frac{m_i v_i^
2}{2kT
}\right)
423 where $k$ is Boltzmann's constant (see
\chref{defunits
}).
424 To accomplish this, normally distributed random numbers are generated
425 by adding twelve random numbers $R_k$ in the range $
0 \le R_k <
1$ and
426 subtracting
6.0 from their sum. The result is then multiplied by the
427 standard deviation of the velocity distribution $
\sqrt{kT/m_i
}$. Since
428 the resulting total energy will not correspond exactly to the required
429 temperature $T$, a correction is made: first the center-of-mass motion
430 is removed and then all velocities are scaled so that the total
431 energy corresponds exactly to $T$ (see
\eqnref{E-T
}).
432 % Why so complicated? What's wrong with Box-Mueller transforms?
434 \subsubsection{Center-of-mass motion
\index{removing COM motion
}}
435 The
\swapindex{center-of-mass
}{velocity
} is normally set to zero at
436 every step; there is (usually) no net external force acting on the
437 system and the center-of-mass velocity should remain constant. In
438 practice, however, the update algorithm introduces a very slow change in
439 the center-of-mass velocity, and therefore in the total kinetic energy of
440 the system -- especially when temperature coupling is used. If such
441 changes are not quenched, an appreciable center-of-mass motion
442 can develop in long runs, and the temperature will be
443 significantly misinterpreted. Something similar may happen due to overall
444 rotational motion, but only when an isolated cluster is simulated. In
445 periodic systems with filled boxes, the overall rotational motion is
446 coupled to other degrees of freedom and does not cause such problems.
449 \subsection{Neighbor searching
\swapindexquiet{neighbor
}{searching
}}
451 As mentioned in
\chref{ff
}, internal forces are
452 either generated from fixed (static) lists, or from dynamic lists.
453 The latter consist of non-bonded interactions between any pair of particles.
454 When calculating the non-bonded forces, it is convenient to have all
455 particles in a rectangular box.
456 As shown in
\figref{pbc
}, it is possible to transform a
457 triclinic box into a rectangular box.
458 The output coordinates are always in a rectangular box, even when a
459 dodecahedron or triclinic box was used for the simulation.
460 Equation
\ref{eqn:box_rot
} ensures that we can reset particles
461 in a rectangular box by first shifting them with
462 box vector $
{\bf c
}$, then with $
{\bf b
}$ and finally with $
{\bf a
}$.
463 Equations
\ref{eqn:box_shift2
},
\ref{eqn:physicalrc
} and
\ref{eqn:gridrc
}
464 ensure that we can find the
14 nearest triclinic images within
465 a linear combination that does not involve multiples of box vectors.
467 \subsubsection{Pair lists generation
}
468 The non-bonded pair forces need to be calculated only for those pairs
469 $i,j$ for which the distance $r_
{ij
}$ between $i$ and the
470 \swapindex{nearest
}{image
}
471 of $j$ is less than a given cut-off radius $R_c$. Some of the particle
472 pairs that fulfill this criterion are excluded, when their interaction
473 is already fully accounted for by bonded interactions.
{\gromacs}
474 employs a
{\em pair list
} that contains those particle pairs for which
475 non-bonded forces must be calculated. The pair list contains atoms
476 $i$, a displacement vector for atom $i$, and all particles $j$ that
477 are within
\verb'rlist' of this particular image of atom $i$. The
478 list is updated every
\verb'nstlist' steps, where
\verb'nstlist' is
479 typically
10. There is an option to calculate the total non-bonded
480 force on each particle due to all particle in a shell around the
481 list cut-off,
{\ie} at a distance between
\verb'rlist' and
482 \verb'rlistlong'. This force is calculated during the pair list update
483 and retained during
\verb'nstlist' steps.
485 To make the
\normindex{neighbor list
}, all particles that are close
486 (
{\ie} within the neighbor list cut-off) to a given particle must be found.
487 This searching, usually called neighbor search (NS) or pair search,
488 involves periodic boundary conditions and determining the
{\em image
}
489 (see
\secref{pbc
}). The search algorithm is $O(N)$, although a simpler
490 $O(N^
2)$ algorithm is still available under some conditions.
492 \subsubsection{\normindex{Cut-off schemes
}: group versus Verlet
}
493 From version
4.6,
{\gromacs} supports two different cut-off scheme
494 setups: the original one based on atom groups and one using a Verlet
495 buffer. There are some important differences that affect results,
496 performance and feature support. The group scheme can be made to work
497 (almost) like the Verlet scheme, but this will lead to a decrease in
498 performance. The group scheme is especially fast for water molecules,
499 which are abundant in many simulations.
501 In the group scheme, a neighbor list is generated consisting of pairs
502 of groups of at least one atom. These groups were originally
503 \swapindex{charge
}{group
}s
\ifthenelse{\equal{\gmxlite}{1}}{}{(see
504 \secref{chargegroup
})
}, but with a proper treatment of long-range
505 electrostatics, performance is their only advantage. A pair of groups
506 is put into the neighbor list when their center of geometry is within
507 the cut-off distance. Interactions between all atom pairs (one from
508 each charge group) are calculated for a certain number of MD steps,
509 until the neighbor list is updated. This setup is efficient, as the
510 neighbor search only checks distance between charge group pair, not
511 atom pairs (saves a factor of $
3 \times 3 =
9$ with a three-atom water
512 model) and the non-bonded force kernels can be optimized for, say, a
513 water molecule ``group''. Without explicit buffering, this setup leads
514 to energy drift as some atom pairs which are within the cut-off don't
515 interact and some outside the cut-off do interact. This can be caused
518 \item atoms moving across the cut-off between neighbor search steps, and/or
519 \item for charge groups consisting of more than one atom, atom pairs
520 moving in/out of the cut-off when their charge group center of
521 geometry distance is outside/inside of the cut-off.
523 Explicitly adding a buffer to the neighbor list will remove such
524 artifacts, but this comes at a high computational cost. How severe the
525 artifacts are depends on the system, the properties in which you are
526 interested, and the cut-off setup.
528 The Verlet cut-off scheme uses a buffered pair list by default. It
529 also uses clusters of atoms, but these are not static as in the group
530 scheme. Rather, the clusters are defined spatially and consist of
4 or
531 8 atoms, which is convenient for stream computing, using e.g. SSE, AVX
532 or CUDA on GPUs. At neighbor search steps, an atom pair list (or
533 cluster pair list, but that's an implementation detail) is created
534 with a Verlet buffer. Thus the pair-list cut-off is larger than the
535 interaction cut-off. In the non-bonded force kernels, forces are only
536 added when an atom pair is within the cut-off distance at that
537 particular time step. This ensures that as atoms move between pair
538 search steps, forces between nearly all atoms within the cut-off
539 distance are calculated. We say
{\em nearly
} all atoms, because
540 {\gromacs} uses a fixed pair list update frequency for
541 efficiency. There is a small chance that an atom pair distance is
542 decreased to within the cut-off in this fixed number of steps. This
543 small chance results in a small energy drift. When temperature
544 coupling is used, the buffer size can be determined automatically,
545 given a certain limit on the energy drift.
547 The Verlet scheme specific settings in the
{\tt mdp
} file are:
549 cutoff-scheme = Verlet
550 verlet-buffer-tolerance =
0.005
552 The Verlet buffer size is determined from the latter option, which is
553 by default set to
0.005 kJ/mol/ps pair energy error per atom. Note that
554 errors in pair energies cancel and the effect on the total energy drift
555 is usually at least an order of magnitude smaller than the tolerance.
556 Furthermore, the drift of the total energy is affected by many other
557 factors, the constraint contribution is often the dominating one.
558 For constant energy (NVE) simulations, this drift should be set to -
1
559 and a buffer has to be set manually by specifying
{\tt rlist
} $>$
{\tt
560 rcoulomb
}. The simplest way to get a reasonable buffer size is to
561 use an NVT
{\tt mdp
} file with the target temperature set to what you
562 expect in your NVE simulation, and transfer the buffer size printed by
563 {\tt grompp
} to your NVE
{\tt mdp
} file.
565 The Verlet cut-off scheme is implemented in a very efficient fashion
566 based on clusters of particles. The simplest example is a cluster size
567 of
4 particles. The pair list is then constructed based on cluster
568 pairs. The cluster-pair search is much faster searching based on
569 particle pairs, because $
4 \times 4 =
16$ particle pairs are put in
570 the list at once. The non-bonded force calculation kernel can then
571 calculate all
16 particle-pair interactions at once, which maps nicely
572 to SIMD units which can perform multiple floating operations at once
573 (e.g. SSE, AVX, CUDA on GPUs, BlueGene FPUs). These non-bonded kernels
574 are much faster than the kernels used in the group scheme for most
575 types of systems, except for water molecules when not using a buffered
576 pair list. This latter case is quite common for (bio-)molecular
577 simulations, so for greatest speed, it is worth comparing the
578 performance of both schemes.
580 As the Verlet cut-off scheme was introduced in version
4.6, not
581 all features of the group scheme are supported yet. The Verlet scheme
582 supports a few new features which the group scheme does not support.
583 A list of features not (fully) supported in both cut-off schemes is
584 given in
\tabref{cutoffschemesupport
}.
588 \begin{tabular
}{|l|c|c|
}
590 Non-bonded interaction feature & group & Verlet \\
592 unbuffered cut-off scheme & $
\surd$ & \\
593 exact cut-off & shift/switch & $
\surd$ \\
594 shifted interactions & force+energy & energy \\
595 switched forces & $
\surd$ & \\
596 non-periodic systems & $
\surd$ & Z + walls \\
597 implicit solvent & $
\surd$ & \\
598 free energy perturbed non-bondeds & $
\surd$ & \\
599 group energy contributions & $
\surd$ & CPU (not on GPU) \\
600 energy group exclusions & $
\surd$ & \\
601 AdResS multi-scale & $
\surd$ & \\
602 OpenMP multi-threading & only PME & $
\surd$ \\
603 native GPU support & & $
\surd$ \\
607 \caption{Differences (only) in the support of non-bonded features
608 between the group and Verlet cut-off schemes.
}
609 \label{tab:cutoffschemesupport
}
612 \ifthenelse{\equal{\gmxlite}{1}}{}{
613 \subsubsection{Energy drift and pair-list buffering
}
614 For a canonical ensemble, the average energy error caused by the
615 finite Verlet buffer size can be determined from the atomic
616 displacements and the shape of the potential at the cut-off.
617 %Since we are interested in the small drift regime, we will assume
618 %#that atoms will only move within the cut-off distance in the last step,
619 %$n_\mathrm{ps}-1$, of the pair list update interval $n_\mathrm{ps}$.
620 %Over this number of steps the displacment of an atom with mass $m$
621 The displacement distribution along one dimension for a freely moving
622 particle with mass $m$ over time $t$ at temperature $T$ is Gaussian
623 with zero mean and variance $
\sigma^
2 = t\,k_B T/m$. For the distance
624 between two atoms, the variance changes to $
\sigma^
2 =
\sigma_{12}^
2 =
625 t\,k_B T(
1/m_1+
1/m_2)$. Note that in practice particles usually
626 interact with other particles over time $t$ and therefore the real
627 displacement distribution is much narrower. Given a non-bonded
628 interaction cut-off distance of $r_c$ and a pair-list cut-off
629 $r_
\ell=r_c+r_b$, we can then write the average energy error after
630 time $t$ for pair interactions between one particle of type
1
631 surrounded by particles of type
2 with number density $
\rho_2$, when
632 the inter particle distance changes from $r_0$ to $r_t$, as:
635 \langle \Delta V
\rangle \! &=&
636 \int_{0}^
{r_c
} \int_{r_
\ell}^
\infty 4 \pi r_0^
2 \rho_2 V(r_t) G\!
\left(
\frac{r_t-r_0
}{\sigma}\right) d r_0\, d r_t \\
638 \int_{-
\infty}^
{r_c
} \int_{r_
\ell}^
\infty 4 \pi r_0^
2 \rho_2 \Big[ V'(r_c) (r_t - r_c) +
641 \phantom{\int_{-
\infty}^
{r_c
} \int_{r_
\ell}^
\infty 4 \pi r_0^
2 \rho_2 \Big[}
642 V''(r_c)
\frac{1}{2}(r_t - r_c)^
2 \Big] G\!
\left(
\frac{r_t-r_0
}{\sigma}\right) d r_0 \, d r_t\\
644 4 \pi (r_
\ell+
\sigma)^
2 \rho_2
645 \int_{-
\infty}^
{r_c
} \int_{r_
\ell}^
\infty \Big[ V'(r_c) (r_t - r_c) +
648 \phantom{4 \pi (r_
\ell+
\sigma)^
2 \rho_2 \int_{-
\infty}^
{r_c
} \int_{r_
\ell}^
\infty \Big[}
649 V''(r_c)
\frac{1}{2}(r_t - r_c)^
2 \Big] G\!
\left(
\frac{r_t-r_0
}{\sigma}\right) d r_0 \, d r_t\\
651 4 \pi (r_
\ell+
\sigma)^
2 \rho_2 \bigg\
{
652 \frac{1}{2}V'(r_c)
\left[r_b
\sigma G\!
\left(
\frac{r_b
}{\sigma}\right) - (r_b^
2+
\sigma^
2)E\!
\left(
\frac{r_b
}{\sigma}\right)
\right] +
655 \phantom{4 \pi (r_
\ell+
\sigma)^
2 \rho_2 \bigg\
{ }
656 \frac{1}{6}V''(r_c)
\left[ \sigma(r_b^
2+
\sigma^
2)G\!
\left(
\frac{r_b
}{\sigma}\right) - r_b(r_b^
2+
3\sigma^
2 ) E\!
\left(
\frac{r_b
}{\sigma}\right)
\right]
660 where $G$ is a Gaussian distribution with
0 mean and unit variance and
661 $E(x)=
\frac{1}{2}\mathrm{erfc
}(x/
\sqrt{2})$. We always want to achieve
662 small energy error, so $
\sigma$ will be small compared to both $r_c$
663 and $r_
\ell$, thus the approximations in the equations above are good,
664 since the Gaussian distribution decays rapidly. The energy error needs
665 to be averaged over all particle pair types and weighted with the
666 particle counts. In
{\gromacs} we don't allow cancellation of error
667 between pair types, so we average the absolute values. To obtain the
668 average energy error per unit time, it needs to be divided by the
669 neighbor-list life time $t = (
{\tt nstlist
} -
1)
\times{\tt dt
}$. This
670 function can not be inverted analytically, so we use bisection to
671 obtain the buffer size $r_b$ for a target drift. Again we note that
672 in practice the error we usually be much smaller than this estimate,
673 as in the condensed phase particle displacements will be much smaller
674 than for freely moving particles, which is the assumption used here.
676 When (bond) constraints are present, some particles will have fewer
677 degrees of freedom. This will reduce the energy errors. The
678 displacement in an arbitrary direction of a particle with
2 degrees of
679 freedom is not Gaussian, but rather follows the complementary error
681 \frac{\sqrt{\pi}}{2\sqrt{2}\sigma}\,
\mathrm{erfc
}\left(
\frac{|r|
}{\sqrt{2}\,
\sigma}\right)
682 \eeq where $
\sigma^
2$ is again $k_B T/m$. This distribution can no
683 longer be integrated analytically to obtain the energy error. But we
684 can generate a tight upper bound using a scaled and shifted Gaussian
685 distribution (not shown). This Gaussian distribution can then be used
686 to calculate the energy error as described above. We consider
687 particles constrained, i.e. having
2 degrees of freedom or fewer, when
688 they are connected by constraints to particles with a total mass of at
689 least
1.5 times the mass of the particles itself. For a particle with
690 a single constraint this would give a total mass along the constraint
691 direction of at least
2.5, which leads to a reduction in the variance
692 of the displacement along that direction by at least a factor of
6.25.
693 As the Gaussian distribution decays very rapidly, this effectively
694 removes one degree of freedom from the displacement. Multiple
695 constraints would reduce the displacement even more, but as this gets
696 very complex, we consider those as particles with
2 degrees of
699 There is one important implementation detail that reduces the energy
700 errors caused by the finite Verlet buffer list size. The derivation
701 above assumes a particle pair-list. However, the
{\gromacs}
702 implementation uses a cluster pair-list for efficiency. The pair list
703 consists of pairs of clusters of
4 particles in most cases, also
704 called a $
4 \times 4$ list, but the list can also be $
4 \times 8$ (GPU
705 CUDA kernels and AVX
256-bit single precision kernels) or $
4 \times 2$
706 (SSE double-precision kernels). This means that the pair-list is
707 effectively much larger than the corresponding $
1 \times 1$ list. Thus
708 slightly beyond the pair-list cut-off there will still be a large
709 fraction of particle pairs present in the list. This fraction can be
710 determined in a simulation and accurately estimated under some
711 reasonable assumptions. The fraction decreases with increasing
712 pair-list range, meaning that a smaller buffer can be used. For
713 typical all-atom simulations with a cut-off of
0.9 nm this fraction is
714 around
0.9, which gives a reduction in the energy errors of a factor of
715 10. This reduction is taken into account during the automatic Verlet
716 buffer calculation and results in a smaller buffer size.
719 \centerline{\includegraphics[width=
9cm
]{plots/verlet-drift
}}
720 \caption {Energy drift per atom for an SPC/E water system at
300K with
721 a time step of
2 fs and a pair-list update period of
10 steps
722 (pair-list life time:
18 fs). PME was used with
{\tt ewald-rtol
} set
723 to
10$^
{-
5}$; this parameter affects the shape of the potential at
724 the cut-off. Error estimates due to finite Verlet buffer size are
725 shown for a $
1 \times 1$ atom pair list and $
4 \times 4$ atom pair
726 list without and with (dashed line) cancellation of positive and
727 negative errors. Real energy drift is shown for double- and
728 single-precision simulations. Single-precision rounding errors in
729 the SETTLE constraint algorithm cause the drift to become negative
730 at large buffer size. Note that at zero buffer size, the real drift
731 is small because positive (H-H) and negative (O-H) energy errors
733 \label{fig:verletdrift
}
736 In
\figref{verletdrift
} one can see that for small buffer sizes the drift
737 of the total energy is much smaller than the pair energy error tolerance,
738 due to cancellation of errors. For larger buffer size, the error estimate
739 is a factor of
6 higher than drift of the total energy, or alternatively
740 the buffer estimate is
0.024 nm too large. This is because the protons
741 don't move freely over
18 fs, but rather vibrate.
742 %At a buffer size of zero there is cancellation of
743 %drift due to repulsive (H-H) and attractive (O-H) interactions.
745 \subsubsection{Cut-off artifacts and switched interactions
}
746 With the Verlet scheme, the pair potentials are shifted to be zero at
747 the cut-off, such that the potential is the integral of the force.
748 Note that in the group scheme this is not possible, because no exact
749 cut-off distance is used. There can still be energy drift from
750 non-zero forces at the cut-off. This effect is extremely small and
751 often not noticeable, as other integration errors may dominate. To
752 completely avoid cut-off artifacts, the non-bonded forces can be
753 switched exactly to zero at some distance smaller than the neighbor
754 list cut-off (there are several ways to do this in
{\gromacs}, see
755 \secref{mod_nb_int
}). One then has a buffer with the size equal to the
756 neighbor list cut-off less the longest interaction cut-off. With the
757 group cut-off scheme, one can then also choose to let
{\tt mdrun
} only
758 update the neighbor list when required. That is when one or more
759 particles have moved more than half the buffer size from the center of
760 geometry of the
\swapindex{charge
}{group
} to which they belong (see
761 \secref{chargegroup
}), as determined at the previous neighbor search.
762 This option guarantees that there are no cut-off artifacts.
{\bf
763 Note
} that for larger systems this comes at a high computational
764 cost, since the neighbor list update frequency will be determined by
765 just one or two particles moving slightly beyond the half buffer
766 length (which not even necessarily implies that the neighbor list is
767 invalid), while
99.99\% of the particles are fine.
} % Brace matches
768 ifthenelse test for gmxlite
770 \subsubsection{Simple search
\swapindexquiet{simple
}{search
}}
771 Due to
\eqnsref{box_rot
}{simplerc
}, the vector $
\rvij$
772 connecting images within the cut-off $R_c$ can be found by constructing:
774 \ve{r
}''' & = &
\ve{r
}_j-
\ve{r
}_i \\
775 \ve{r
}'' & = &
\ve{r
}''' -
{\bf c
}*
\verb'round'(r'''_z/c_z) \\
776 \ve{r
}' & = &
\ve{r
}'' -
{\bf b
}*
\verb'round'(r''_y/b_y) \\
777 \ve{r
}_
{ij
} & = &
\ve{r
}' -
{\bf a
}*
\verb'round'(r'_x/a_x)
779 When distances between two particles in a triclinic box are needed
780 that do not obey
\eqnref{box_rot
},
781 many shifts of combinations of box vectors need to be considered to find
784 \ifthenelse{\equal{\gmxlite}{1}}{}{
787 \centerline{\includegraphics[width=
8cm
]{plots/nstric
}}
788 \caption {Grid search in two dimensions. The arrows are the box vectors.
}
792 \subsubsection{Grid search
\swapindexquiet{grid
}{search
}}
794 The grid search is schematically depicted in
\figref{grid
}. All
795 particles are put on the
{\nsgrid}, with the smallest spacing $
\ge$
796 $R_c/
2$ in each of the directions. In the direction of each box
797 vector, a particle $i$ has three images. For each direction the image
798 may be -
1,
0 or
1, corresponding to a translation over -
1,
0 or +
1 box
799 vector. We do not search the surrounding
{\nsgrid} cells for neighbors
800 of $i$ and then calculate the image, but rather construct the images
801 first and then search neighbors corresponding to that image of $i$.
802 As
\figref{grid
} shows, some grid cells may be searched more than once
803 for different images of $i$. This is not a problem, since, due to the
804 minimum image convention, at most one image will ``see'' the
805 $j$-particle. For every particle, fewer than
125 (
5$^
3$) neighboring
806 cells are searched. Therefore, the algorithm scales linearly with the
807 number of particles. Although the prefactor is large, the scaling
808 behavior makes the algorithm far superior over the standard $O(N^
2)$
809 algorithm when there are more than a few hundred particles. The
810 grid search is equally fast for rectangular and triclinic boxes. Thus
811 for most protein and peptide simulations the rhombic dodecahedron will
812 be the preferred box shape.
813 } % Brace matches ifthenelse test for gmxlite
815 \ifthenelse{\equal{\gmxlite}{1}}{}{
816 \subsubsection{Charge groups
}
817 \label{sec:chargegroup
}\swapindexquiet{charge
}{group
}%
818 Charge groups were originally introduced to reduce cut-off artifacts
819 of Coulomb interactions. When a plain cut-off is used, significant
820 jumps in the potential and forces arise when atoms with (partial) charges
821 move in and out of the cut-off radius. When all chemical moieties have
822 a net charge of zero, these jumps can be reduced by moving groups
823 of atoms with net charge zero, called charge groups, in and
824 out of the neighbor list. This reduces the cut-off effects from
825 the charge-charge level to the dipole-dipole level, which decay
826 much faster. With the advent of full range electrostatics methods,
827 such as particle mesh Ewald (
\secref{pme
}), the use of charge groups is
828 no longer required for accuracy. It might even have a slight negative effect
829 on the accuracy or efficiency, depending on how the neighbor list is made
830 and the interactions are calculated.
832 But there is still an important reason for using ``charge groups'': efficiency.
833 Where applicable, neighbor searching is carried out on the basis of
834 charge groups which are defined in the molecular topology.
835 If the nearest image distance between the
{\em
836 geometrical centers
} of the atoms of two charge groups is less than
837 the cut-off radius, all atom pairs between the charge groups are
838 included in the pair list.
839 The neighbor searching for a water system, for instance,
840 is $
3^
2=
9$ times faster when each molecule is treated as a charge group.
841 Also the highly optimized water force loops (see
\secref{waterloops
})
842 only work when all atoms in a water molecule form a single charge group.
843 Currently the name
{\em neighbor-search group
} would be more appropriate,
844 but the name charge group is retained for historical reasons.
845 When developing a new force field, the advice is to use charge groups
846 of
3 to
4 atoms for optimal performance. For all-atom force fields
847 this is relatively easy, as one can simply put hydrogen atoms, and in some
848 case oxygen atoms, in the same charge group as the heavy atom they
849 are connected to; for example: CH$_3$, CH$_2$, CH, NH$_2$, NH, OH, CO$_2$, CO.
850 } % Brace matches ifthenelse test for gmxlite
852 \subsection{Compute forces
}
853 \label{subsec:forces
}
855 \subsubsection{Potential energy
}
856 When forces are computed, the
\swapindex{potential
}{energy
} of each
857 interaction term is computed as well. The total potential energy is
858 summed for various contributions, such as Lennard-Jones, Coulomb, and
859 bonded terms. It is also possible to compute these contributions for
860 {\em energy-monitor groups
} of atoms that are separately defined (see
861 \secref{groupconcept
}).
863 \subsubsection{Kinetic energy and temperature
}
864 The
\normindex{temperature
} is given by the total
865 \swapindex{kinetic
}{energy
} of the $N$-particle system:
867 E_
{kin
} =
\half \sum_{i=
1}^N m_i v_i^
2
869 From this the absolute temperature $T$ can be computed using:
871 \half N_
{df
} kT = E_
{kin
}
874 where $k$ is Boltzmann's constant and $N_
{df
}$ is the number of
875 degrees of freedom which can be computed from:
877 N_
{df
} ~=~
3 N - N_c - N_
{com
}
879 Here $N_c$ is the number of
{\em \normindex{constraints
}} imposed on the system.
880 When performing molecular dynamics $N_
{com
}=
3$ additional degrees of
881 freedom must be removed, because the three
882 center-of-mass velocities are constants of the motion, which are usually
883 set to zero. When simulating in vacuo, the rotation around the center of mass
884 can also be removed, in this case $N_
{com
}=
6$.
885 When more than one temperature-coupling group
\index{temperature-coupling group
} is used, the number of degrees
886 of freedom for group $i$ is:
888 N^i_
{df
} ~=~ (
3 N^i - N^i_c)
\frac{3 N - N_c - N_
{com
}}{3 N - N_c
}
891 The kinetic energy can also be written as a tensor, which is necessary
892 for pressure calculation in a triclinic system, or systems where shear
895 {\bf E
}_
{kin
} =
\half \sum_i^N m_i
\vvi \otimes \vvi
898 \subsubsection{Pressure and virial
}
899 The
\normindex{pressure
}
900 tensor
{\bf P
} is calculated from the difference between
901 kinetic energy $E_
{kin
}$ and the
\normindex{virial
} $
{\bf \Xi}$:
903 {\bf P
} =
\frac{2}{V
} (
{\bf E
}_
{kin
}-
{\bf \Xi})
906 where $V$ is the volume of the computational box.
907 The scalar pressure $P$, which can be used for pressure coupling in the case
908 of isotropic systems, is computed as:
910 P =
{\rm trace
}(
{\bf P
})/
3
913 The virial $
{\bf \Xi}$ tensor is defined as:
915 {\bf \Xi} = -
\half \sum_{i<j
} \rvij \otimes \Fvij
919 \ifthenelse{\equal{\gmxlite}{1}}{}{
920 The
{\gromacs} implementation of the virial computation is described
922 } % Brace matches ifthenelse test for gmxlite
925 \subsection{The
\swapindex{leap-frog
}{integrator
}}
926 \label{subsec:update
}
928 \centerline{\includegraphics[width=
8cm
]{plots/leapfrog
}}
929 \caption[The Leap-Frog integration method.
]{The Leap-Frog integration method. The algorithm is called Leap-Frog because $
\ve{r
}$ and $
\ve{v
}$ are leaping
930 like frogs over each other's backs.
}
934 The default MD integrator in
{\gromacs} is the so-called
{\em leap-frog
}
935 algorithm~
\cite{Hockney74
} for the integration of the equations of
936 motion. When extremely accurate integration with temperature
937 and/or pressure coupling is required, the velocity Verlet integrators
938 are also present and may be preferable (see
\ssecref{vverlet
}). The leap-frog
939 algorithm uses positions $
\ve{r
}$ at time $t$ and
940 velocities $
\ve{v
}$ at time $t-
\hDt$; it updates positions and
941 velocities using the forces
942 $
\ve{F
}(t)$ determined by the positions at time $t$ using these relations:
944 \label{eqn:leapfrogv
}
945 \ve{v
}(t+
\hDt) &~=~&
\ve{v
}(t-
\hDt)+
\frac{\Dt}{m
}\ve{F
}(t) \\
946 \ve{r
}(t+
\Dt) &~=~&
\ve{r
}(t)+
\Dt\ve{v
}(t+
\hDt)
948 The algorithm is visualized in
\figref{leapfrog
}.
949 It produces trajectories that are identical to the Verlet~
\cite{Verlet67
} algorithm,
950 whose position-update relation is
952 \ve{r
}(t+
\Dt)~=~
2\ve{r
}(t) -
\ve{r
}(t-
\Dt) +
\frac{1}{m
}\ve{F
}(t)
\Dt^
2+O(
\Dt^
4)
954 The algorithm is of third order in $
\ve{r
}$ and is time-reversible.
955 See ref.~
\cite{Berendsen86b
} for the merits of this algorithm and comparison
956 with other time integration algorithms.
958 The
\swapindex{equations of
}{motion
} are modified for temperature
959 coupling and pressure coupling, and extended to include the
960 conservation of constraints, all of which are described below.
962 \subsection{The
\swapindex{velocity Verlet
}{integrator
}}
963 \label{subsec:vverlet
}
964 The velocity Verlet algorithm~
\cite{Swope82
} is also implemented in
965 {\gromacs}, though it is not yet fully integrated with all sets of
966 options. In velocity Verlet, positions $
\ve{r
}$ and velocities
967 $
\ve{v
}$ at time $t$ are used to integrate the equations of motion;
968 velocities at the previous half step are not required.
\bea
969 \label{eqn:velocityverlet1
}
970 \ve{v
}(t+
\hDt) &~=~&
\ve{v
}(t)+
\frac{\Dt}{2m
}\ve{F
}(t) \\
971 \ve{r
}(t+
\Dt) &~=~&
\ve{r
}(t)+
\Dt\,
\ve{v
}(t+
\hDt) \\
972 \ve{v
}(t+
\Dt) &~=~&
\ve{v
}(t+
\hDt)+
\frac{\Dt}{2m
}\ve{F
}(t+
\Dt)
976 \label{eqn:velocityverlet2
}
977 \ve{r
}(t+
\Dt) &~=~&
\ve{r
}(t)+
\Dt\,
\ve{v
} +
\frac{\Dt^
2}{2m
}\ve{F
}(t) \\
978 \ve{v
}(t+
\Dt) &~=~&
\ve{v
}(t)+
\frac{\Dt}{2m
}\left[\ve{F
}(t) +
\ve{F
}(t+
\Dt)
\right]
980 With no temperature or pressure coupling, and with
{\em corresponding
}
981 starting points, leap-frog and velocity Verlet will generate identical
982 trajectories, as can easily be verified by hand from the equations
983 above. Given a single starting file with the
{\em same
} starting
984 point $
\ve{x
}(
0)$ and $
\ve{v
}(
0)$, leap-frog and velocity Verlet will
985 {\em not
} give identical trajectories, as leap-frog will interpret the
986 velocities as corresponding to $t=-
\hDt$, while velocity Verlet will
987 interpret them as corresponding to the timepoint $t=
0$.
989 \subsection{Understanding reversible integrators: The Trotter decomposition
}
990 To further understand the relationship between velocity Verlet and
991 leap-frog integration, we introduce the reversible Trotter formulation
992 of dynamics, which is also useful to understanding implementations of
993 thermostats and barostats in
{\gromacs}.
995 A system of coupled, first-order differential equations can be evolved
996 from time $t =
0$ to time $t$ by applying the evolution operator
998 \Gamma(t) &=&
\exp(iLt)
\Gamma(
0)
\nonumber \\
999 iL &=&
\dot{\Gamma}\cdot \nabla_{\Gamma},
1001 where $L$ is the Liouville operator, and $
\Gamma$ is the
1002 multidimensional vector of independent variables (positions and
1004 A short-time approximation to the true operator, accurate at time $
\Dt
1005 = t/P$, is applied $P$ times in succession to evolve the system as
1007 \Gamma(t) =
\prod_{i=
1}^P
\exp(iL
\Dt)
\Gamma(
0)
1009 For NVE dynamics, the Liouville operator is
1011 iL =
\sum_{i=
1}^
{N
} \vv_i \cdot \nabla_{\rv_i} +
\sum_{i=
1}^N
\frac{1}{m_i
}\F(r_i)
\cdot \nabla_{\vv_i}.
1013 This can be split into two additive operators
1015 iL_1 &=&
\sum_{i=
1}^N
\frac{1}{m_i
}\F(r_i)
\cdot \nabla_{\vv_i} \nonumber \\
1016 iL_2 &=&
\sum_{i=
1}^
{N
} \vv_i \cdot \nabla_{\rv_i}
1018 Then a short-time, symmetric, and thus reversible approximation of the true dynamics will be
1020 \exp(iL
\Dt) =
\exp(iL_2
\hDt)
\exp(iL_1
\Dt)
\exp(iL_2
\hDt) +
\mathcal{O
}(
\Dt^
3).
1021 \label{eq:NVE_Trotter
}
1023 This corresponds to velocity Verlet integration. The first
1024 exponential term over $
\hDt$ corresponds to a velocity half-step, the
1025 second exponential term over $
\Dt$ corresponds to a full velocity
1026 step, and the last exponential term over $
\hDt$ is the final velocity
1027 half step. For future times $t = n
\Dt$, this becomes
1029 \exp(iLn
\Dt) &
\approx&
\left(
\exp(iL_2
\hDt)
\exp(iL_1
\Dt)
\exp(iL_2
\hDt)
\right)^n
\nonumber \\
1030 &
\approx&
\exp(iL_2
\hDt)
\bigg(
\exp(iL_1
\Dt)
\exp(iL_2
\Dt)
\bigg)^
{n-
1} \nonumber \\
1031 & & \;\;\;\;
\exp(iL_1
\Dt)
\exp(iL_2
\hDt)
1033 This formalism allows us to easily see the difference between the
1034 different flavors of Verlet integrators. The leap-frog integrator can
1035 be seen as starting with Eq.~
\ref{eq:NVE_Trotter
} with the
1036 $
\exp\left(iL_1
\dt\right)$ term, instead of the half-step velocity
1039 \exp(iLn
\dt) &=&
\exp\left(iL_1
\dt\right)
\exp\left(iL_2
\Dt \right) +
\mathcal{O
}(
\Dt^
3).
1041 Here, the full step in velocity is between $t-
\hDt$ and $t+
\hDt$,
1042 since it is a combination of the velocity half steps in velocity
1043 Verlet. For future times $t = n
\Dt$, this becomes
1045 \exp(iLn
\dt) &
\approx&
\bigg(
\exp\left(iL_1
\dt\right)
\exp\left(iL_2
\Dt \right)
\bigg)^
{n
}.
1047 Although at first this does not appear symmetric, as long as the full velocity
1048 step is between $t-
\hDt$ and $t+
\hDt$, then this is simply a way of
1049 starting velocity Verlet at a different place in the cycle.
1051 Even though the trajectory and thus potential energies are identical
1052 between leap-frog and velocity Verlet, the kinetic energy and
1053 temperature will not necessarily be the same. Standard velocity
1054 Verlet uses the velocities at the $t$ to calculate the kinetic energy
1055 and thus the temperature only at time $t$; the kinetic energy is then a sum over all particles
1057 KE_
{\mathrm{full
}}(t) &=&
\sum_i \left(
\frac{1}{2m_i
}\ve{v
}_i(t)
\right)^
2 \nonumber\\
1058 &=&
\sum_i \frac{1}{2m_i
}\left(
\frac{1}{2}\ve{v
}_i(t-
\hDt)+
\frac{1}{2}\ve{v
}_i(t+
\hDt)
\right)^
2,
1060 with the square on the
{\em outside
} of the average. Standard
1061 leap-frog calculates the kinetic energy at time $t$ based on the
1062 average kinetic energies at the timesteps $t+
\hDt$ and $t-
\hDt$, or
1063 the sum over all particles
1065 KE_
{\mathrm{average
}}(t) &=&
\sum_i \frac{1}{2m_i
}\left(
\frac{1}{2}\ve{v
}_i(t-
\hDt)^
2+
\frac{1}{2}\ve{v
}_i(t+
\hDt)^
2\right),
1067 where the square is
{\em inside
} the average.
1069 A non-standard variant of velocity Verlet which averages the kinetic
1070 energies $KE(t+
\hDt)$ and $KE(t-
\hDt)$, exactly like leap-frog, is also
1071 now implemented in
{\gromacs} (as
{\tt .mdp
} file option
{\tt md-vv-avek
}). Without
1072 temperature and pressure coupling, velocity Verlet with
1073 half-step-averaged kinetic energies and leap-frog will be identical up
1074 to numerical precision. For temperature- and pressure-control schemes,
1075 however, velocity Verlet with half-step-averaged kinetic energies and
1076 leap-frog will be different, as will be discussed in the section in
1077 thermostats and barostats.
1079 The half-step-averaged kinetic energy and temperature are slightly more
1080 accurate for a given step size; the difference in average kinetic
1081 energies using the half-step-averaged kinetic energies (
{\em md
} and
1082 {\em md-vv-avek
}) will be closer to the kinetic energy obtained in the
1083 limit of small step size than will the full-step kinetic energy (using
1084 {\em md-vv
}). For NVE simulations, this difference is usually not
1085 significant, since the positions and velocities of the particles are
1086 still identical; it makes a difference in the way the the temperature
1087 of the simulations are
{\em interpreted
}, but
{\em not
} in the
1088 trajectories that are produced. Although the kinetic energy is more
1089 accurate with the half-step-averaged method, meaning that it changes
1090 less as the timestep gets large, it is also more noisy. The RMS deviation
1091 of the total energy of the system (sum of kinetic plus
1092 potential) in the half-step-averaged kinetic energy case will be
1093 higher (about twice as high in most cases) than the full-step kinetic
1094 energy. The drift will still be the same, however, as again, the
1095 trajectories are identical.
1097 For NVT simulations, however, there
{\em will
} be a difference, as
1098 discussed in the section on temperature control, since the velocities
1099 of the particles are adjusted such that kinetic energies of the
1100 simulations, which can be calculated either way, reach the
1101 distribution corresponding to the set temperature. In this case, the
1102 three methods will not give identical results.
1104 Because the velocity and position are both defined at the same time
1105 $t$ the velocity Verlet integrator can be used for some methods,
1106 especially rigorously correct pressure control methods, that are not
1107 actually possible with leap-frog. The integration itself takes
1108 negligibly more time than leap-frog, but twice as many communication
1109 calls are currently required. In most cases, and especially for large
1110 systems where communication speed is important for parallelization and
1111 differences between thermodynamic ensembles vanish in the $
1/N$ limit,
1112 and when only NVT ensembles are required, leap-frog will likely be the
1113 preferred integrator. For pressure control simulations where the fine
1114 details of the thermodynamics are important, only velocity Verlet
1115 allows the true ensemble to be calculated. In either case, simulation
1116 with double precision may be required to get fine details of
1117 thermodynamics correct.
1119 \subsection{Twin-range cut-offs
\index{twin-range!cut-off
}}
1120 To save computation time, slowly varying forces can be calculated
1121 less often than rapidly varying forces. In
{\gromacs}
1122 such a
\normindex{multiple time step
} splitting is possible between
1123 short and long range non-bonded interactions.
1124 In
{\gromacs} versions up to
4.0, an irreversible integration scheme
1125 was used which is also used by the
{\gromos} simulation package:
1126 every $n$ steps the long range forces are determined and these are
1127 then also used (without modification) for the next $n-
1$ integration steps
1128 in
\eqnref{leapfrogv
}. Such an irreversible scheme can result in bad energy
1129 conservation and, possibly, bad sampling.
1130 Since version
4.5, a leap-frog version of the reversible Trotter decomposition scheme~
\cite{Tuckerman1992a
} is used.
1131 In this integrator the long-range forces are determined every $n$ steps
1132 and are then integrated into the velocity in
\eqnref{leapfrogv
} using
1133 a time step of $
\Dt_\mathrm{LR
} = n
\Dt$:
1136 \left\
{ \begin{array
}{lll
} \displaystyle
1137 \ve{v
}(t-
\hDt) +
\frac{1}{m
}\left[\ve{F
}_
\mathrm{SR
}(t) + n
\ve{F
}_
\mathrm{LR
}(t)
\right] \Dt &,&
\mathrm{step
} ~\%~ n =
0 \\
\noalign{\medskip} \displaystyle
1138 \ve{v
}(t-
\hDt) +
\frac{1}{m
}\ve{F
}_
\mathrm{SR
}(t)
\Dt &,&
\mathrm{step
} ~\%~ n
\neq 0 \\
1142 The parameter $n$ is equal to the neighbor list update frequency. In
1143 4.5, the velocity Verlet version of multiple time-stepping is not yet
1146 Several other simulation packages uses multiple time stepping for
1147 bonds and/or the PME mesh forces. In
{\gromacs} we have not implemented
1148 this (yet), since we use a different philosophy. Bonds can be constrained
1149 (which is also a more sound approximation of a physical quantum
1150 oscillator), which allows the smallest time step to be increased
1151 to the larger one. This not only halves the number of force calculations,
1152 but also the update calculations. For even larger time steps, angle vibrations
1153 involving hydrogen atoms can be removed using virtual interaction
1154 \ifthenelse{\equal{\gmxlite}{1}}
1156 {sites (see
\secref{rmfast
}),
}
1157 which brings the shortest time step up to
1158 PME mesh update frequency of a multiple time stepping scheme.
1160 As an example we show the energy conservation for integrating
1161 the equations of motion for SPC/E water at
300 K. To avoid cut-off
1162 effects, reaction-field electrostatics with $
\epsilon_{RF
}=
\infty$ and
1163 shifted Lennard-Jones interactions are used, both with a buffer region.
1164 The long-range interactions were evaluated between
1.0 and
1.4 nm.
1165 In
\figref{leapfrog
} one can see that for electrostatics the Trotter scheme
1166 does an order of magnitude better up to $
\Dt_{LR
}$ =
16 fs.
1167 The electrostatics depends strongly on the orientation of the water molecules,
1168 which changes rapidly.
1169 For Lennard-Jones interactions, the energy drift is linear in $
\Dt_{LR
}$
1170 and roughly two orders of magnitude smaller than for the electrostatics.
1171 Lennard-Jones forces are smaller than Coulomb forces and
1172 they are mainly affected by translation of water molecules, not rotation.
1175 \centerline{\includegraphics[width=
12cm
]{plots/drift-all
}}
1176 \caption{Energy drift per degree of freedom in SPC/E water
1177 with twin-range cut-offs
1178 for reaction field (left) and Lennard-Jones interaction (right)
1179 as a function of the long-range time step length for the irreversible
1180 ``
\gromos'' scheme and a reversible Trotter scheme.
}
1181 \label{fig:twinrangeener
}
1184 \subsection{Temperature coupling
\index{temperature coupling
}}
1185 While direct use of molecular dynamics gives rise to the NVE (constant
1186 number, constant volume, constant energy ensemble), most quantities
1187 that we wish to calculate are actually from a constant temperature
1188 (NVT) ensemble, also called the canonical ensemble.
{\gromacs} can use
1189 the
{\em weak-coupling
} scheme of Berendsen~
\cite{Berendsen84
},
1190 stochastic randomization through the Andersen
1191 thermostat~
\cite{Andersen80
}, the extended ensemble Nos
{\'e
}-Hoover
1192 scheme~
\cite{Nose84,Hoover85
}, or a velocity-rescaling
1193 scheme~
\cite{Bussi2007a
} to simulate constant temperature, with
1194 advantages of each of the schemes laid out below.
1196 There are several other reasons why it might be necessary to control
1197 the temperature of the system (drift during equilibration, drift as a
1198 result of force truncation and integration errors, heating due to
1199 external or frictional forces), but this is not entirely correct to do
1200 from a thermodynamic standpoint, and in some cases only masks the
1201 symptoms (increase in temperature of the system) rather than the
1202 underlying problem (deviations from correct physics in the dynamics).
1203 For larger systems, errors in ensemble averages and structural
1204 properties incurred by using temperature control to remove slow drifts
1205 in temperature appear to be negligible, but no completely
1206 comprehensive comparisons have been carried out, and some caution must
1207 be taking in interpreting the results.
1209 \subsubsection{Berendsen temperature coupling
\pawsindexquiet{Berendsen
}{temperature coupling
}\index{weak coupling
}}
1210 The Berendsen algorithm mimics weak coupling with first-order
1211 kinetics to an external heat bath with given temperature $T_0$.
1212 See ref.~
\cite{Berendsen91
} for a comparison with the
1213 Nos
{\'e
}-Hoover scheme. The effect of this algorithm is
1214 that a deviation of the system temperature from $T_0$ is slowly
1215 corrected according to:
1217 \frac{\de T
}{\de t
} =
\frac{T_0-T
}{\tau}
1218 \label{eqn:Tcoupling
}
1220 which means that a temperature deviation decays exponentially with a
1221 time constant $
\tau$.
1222 This method of coupling has the advantage that the strength of the
1223 coupling can be varied and adapted to the user requirement: for
1224 equilibration purposes the coupling time can be taken quite short
1225 (
{\eg} 0.01 ps), but for reliable equilibrium runs it can be taken much
1226 longer (
{\eg} 0.5 ps) in which case it hardly influences the
1227 conservative dynamics.
1229 The Berendsen thermostat suppresses the fluctuations of the kinetic
1230 energy. This means that one does not generate a proper canonical
1231 ensemble, so rigorously, the sampling will be incorrect. This
1232 error scales with $
1/N$, so for very large systems most ensemble
1233 averages will not be affected significantly, except for the
1234 distribution of the kinetic energy itself. However, fluctuation
1235 properties, such as the heat capacity, will be affected. A similar
1236 thermostat which does produce a correct ensemble is the velocity
1237 rescaling thermostat~
\cite{Bussi2007a
} described below.
1239 The heat flow into or out of the system is affected by scaling the
1240 velocities of each particle every step, or every $n_
\mathrm{TC
}$ steps,
1241 with a time-dependent factor $
\lambda$, given by:
1243 \lambda =
\left[ 1 +
\frac{n_
\mathrm{TC
} \Delta t
}{\tau_T}
1244 \left\
{\frac{T_0
}{T(t -
\hDt)
} -
1 \right\
} \right]^
{1/
2}
1247 The parameter $
\tau_T$ is close, but not exactly equal, to the time constant
1248 $
\tau$ of the temperature coupling (
\eqnref{Tcoupling
}):
1250 \tau =
2 C_V
\tau_T / N_
{df
} k
1252 where $C_V$ is the total heat capacity of the system, $k$ is Boltzmann's
1253 constant, and $N_
{df
}$ is the total number of degrees of freedom. The
1254 reason that $
\tau \neq \tau_T$ is that the kinetic energy change
1255 caused by scaling the velocities is partly redistributed between
1256 kinetic and potential energy and hence the change in temperature is
1257 less than the scaling energy. In practice, the ratio $
\tau /
\tau_T$
1258 ranges from
1 (gas) to
2 (harmonic solid) to
3 (water). When we use
1259 the term ``temperature coupling time constant,'' we mean the parameter
1260 \normindex{$
\tau_T$
}.
1261 {\bf Note
} that in practice the scaling factor $
\lambda$ is limited to
1262 the range of
0.8 $<=
\lambda <=$
1.25, to avoid scaling by very large
1263 numbers which may crash the simulation. In normal use,
1264 $
\lambda$ will always be much closer to
1.0.
1266 \subsubsection{Velocity-rescaling temperature coupling
\pawsindexquiet{velocity-rescaling
}{temperature coupling
}}
1267 The velocity-rescaling thermostat~
\cite{Bussi2007a
} is essentially a Berendsen
1268 thermostat (see above) with an additional stochastic term that ensures
1269 a correct kinetic energy distribution by modifying it according to
1271 \de K = (K_0 - K)
\frac{\de t
}{\tau_T} +
2 \sqrt{\frac{K K_0
}{N_f
}} \frac{\de W
}{\sqrt{\tau_T}},
1272 \label{eqn:vrescale
}
1274 where $K$ is the kinetic energy, $N_f$ the number of degrees of freedom and $
\de W$ a Wiener process.
1275 There are no additional parameters, except for a random seed.
1276 This thermostat produces a correct canonical ensemble and still has
1277 the advantage of the Berendsen thermostat: first order decay of
1278 temperature deviations and no oscillations.
1279 When an $NVT$ ensemble is used, the conserved energy quantity
1280 is written to the energy and log file.
1282 \subsubsection{\normindex{Andersen thermostat
}}
1283 One simple way to maintain a thermostatted ensemble is to take an
1284 $NVE$ integrator and periodically re-select the velocities of the
1285 particles from a Maxwell-Boltzmann distribution.~
\cite{Andersen80
}.
1286 This can either be done by randomizing all the velocities
1287 simultaneously (massive collision) every $
\tau_T/
\Dt$ steps, or by
1288 randomizing every particle with some small probability every timestep,
1289 equal to $
\Dt/
\tau$, where in both cases $
\Dt$ is the timestep and
1290 $
\tau_T$ is a characteristic coupling time scale.
1292 Because of the way constraints operate, all particles in the same
1293 constraint group must be re-randomized simultaneously. This
1294 thermostat is also only possible with velocity Verlet algorithms,
1295 because it operates directly on the velocities at each timestep.
1297 This algorithm avoids some of the ergodicity issues of other
1298 algorithms, as energy cannot flow back and forth between energetically
1299 decoupled components of the system as in velocity scaling motions.
1300 However, it can slow down the kinetics of system by randomizing
1301 correlated motions of the system, including slowing sampling when
1302 $
\tau_T$ is at moderate levels (less than
10 ps). This algorithm
1303 should therefore generally not be used when examining kinetics of the
1304 system, but can avoid ergodicity problems of scaling problems when
1305 examining thermodynamic properties.
1307 % \ifthenelse{\equal{\gmxlite}{1}}{}{
1308 \subsubsection{Nos
{\'e
}-Hoover temperature coupling
\index{Nose-Hoover temperature coupling@Nos
{\'e
}-Hoover temperature coupling|see
{temperature coupling, Nos
{\'e
}-Hoover
}}{\index{temperature coupling Nose-Hoover@temperature coupling Nos
{\'e
}-Hoover
}}\index{extended ensemble
}}
1310 The Berendsen weak-coupling algorithm is
1311 extremely efficient for relaxing a system to the target temperature,
1312 but once the system has reached equilibrium it might be more
1313 important to probe a correct canonical ensemble. This is unfortunately
1314 not the case for the weak-coupling scheme.
1316 To enable canonical ensemble simulations,
{\gromacs} also supports the
1317 extended-ensemble approach first proposed by Nos
{\'e
}~
\cite{Nose84
}
1318 and later modified by Hoover~
\cite{Hoover85
}. The system Hamiltonian is
1319 extended by introducing a thermal reservoir and a friction term in the
1320 equations of motion. The friction force is proportional to the
1321 product of each particle's velocity and a friction parameter, $
\xi$.
1322 This friction parameter (or ``heat bath'' variable) is a fully
1323 dynamic quantity with its own momentum ($p_
{\xi}$) and equation of
1324 motion; the time derivative is calculated from the difference between
1325 the current kinetic energy and the reference temperature.
1327 In this formulation, the particles' equations of motion in
1328 \figref{global
} are replaced by:
1330 \frac {\de^
2\ve{r
}_i
}{\de t^
2} =
\frac{\ve{F
}_i
}{m_i
} -
1331 \frac{p_
{\xi}}{Q
}\frac{\de \ve{r
}_i
}{\de t
} ,
1332 \label{eqn:NH-eqn-of-motion
}
1333 \eeq where the equation of motion for the heat bath parameter $
\xi$ is:
1334 \beq \frac {\de p_
{\xi}}{\de t
} =
\left( T - T_0
\right).
\eeq The
1335 reference temperature is denoted $T_0$, while $T$ is the current
1336 instantaneous temperature of the system. The strength of the coupling
1337 is determined by the constant $Q$ (usually called the ``mass parameter''
1338 of the reservoir) in combination with the reference
1339 temperature.~
\footnote{Note that some derivations, an alternative
1340 notation $
\xi_{\mathrm{alt
}} = v_
{\xi} = p_
{\xi}/Q$ is used.
}
1342 The conserved quantity for the Nos
{\'e
}-Hoover equations of motion is not
1343 the total energy, but rather
1345 H =
\sum_{i=
1}^
{N
} \frac{\pb_i}{2m_i
} + U
\left(
\rv_1,
\rv_2,
\ldots,
\rv_N\right) +
\frac{p_
{\xi}^
2}{2Q
} + N_fkT
\xi,
1347 where $N_f$ is the total number of degrees of freedom.
1349 In our opinion, the mass parameter is a somewhat awkward way of
1350 describing coupling strength, especially due to its dependence on
1351 reference temperature (and some implementations even include the
1352 number of degrees of freedom in your system when defining $Q$). To
1353 maintain the coupling strength, one would have to change $Q$ in
1354 proportion to the change in reference temperature. For this reason, we
1355 prefer to let the
{\gromacs} user work instead with the period
1356 $
\tau_T$ of the oscillations of kinetic energy between the system and
1357 the reservoir instead. It is directly related to $Q$ and $T_0$ via:
1359 Q =
\frac {\tau_T^
2 T_0
}{4 \pi^
2}.
1361 This provides a much more intuitive way of selecting the
1362 Nos
{\'e
}-Hoover coupling strength (similar to the weak-coupling
1363 relaxation), and in addition $
\tau_T$ is independent of system size
1364 and reference temperature.
1366 It is however important to keep the difference between the
1367 weak-coupling scheme and the Nos
{\'e
}-Hoover algorithm in mind:
1368 Using weak coupling you get a
1369 strongly damped
{\em exponential relaxation
},
1370 while the Nos
{\'e
}-Hoover approach
1371 produces an
{\em oscillatory relaxation
}.
1372 The actual time it takes to relax with Nos
{\'e
}-Hoover coupling is
1373 several times larger than the period of the
1374 oscillations that you select. These oscillations (in contrast
1375 to exponential relaxation) also means that
1376 the time constant normally should be
4--
5 times larger
1377 than the relaxation time used with weak coupling, but your
1380 Nos
{\'e
}-Hoover dynamics in simple systems such as collections of
1381 harmonic oscillators, can be
{\em nonergodic
}, meaning that only a
1382 subsection of phase space is ever sampled, even if the simulations
1383 were to run for infinitely long. For this reason, the Nos
{\'e
}-Hoover
1384 chain approach was developed, where each of the Nos
{\'e
}-Hoover
1385 thermostats has its own Nos
{\'e
}-Hoover thermostat controlling its
1386 temperature. In the limit of an infinite chain of thermostats, the
1387 dynamics are guaranteed to be ergodic. Using just a few chains can
1388 greatly improve the ergodicity, but recent research has shown that the
1389 system will still be nonergodic, and it is still not entirely clear
1390 what the practical effect of this~
\cite{Cooke2008
}. Currently, the
1391 default number of chains is
10, but this can be controlled by the
1392 user. In the case of chains, the equations are modified in the
1393 following way to include a chain of thermostatting
1394 particles~
\cite{Martyna1992
}:
1397 \frac {\de^
2\ve{r
}_i
}{\de t^
2} &~=~&
\frac{\ve{F
}_i
}{m_i
} -
\frac{p_
{{\xi}_1
}}{Q_1
} \frac{\de \ve{r
}_i
}{\de t
} \nonumber \\
1398 \frac {\de p_
{{\xi}_1
}}{\de t
} &~=~&
\left( T - T_0
\right) - p_
{{\xi}_1
} \frac{p_
{{\xi}_2
}}{Q_2
} \nonumber \\
1399 \frac {\de p_
{{\xi}_
{i=
2\ldots N
}}}{\de t
} &~=~&
\left(
\frac{p_
{\xi_{i-
1}}^
2}{Q_
{i-
1}} -kT
\right) - p_
{\xi_i} \frac{p_
{\xi_{i+
1}}}{Q_
{i+
1}} \nonumber \\
1400 \frac {\de p_
{\xi_N}}{\de t
} &~=~&
\left(
\frac{p_
{\xi_{N-
1}}^
2}{Q_
{N-
1}}-kT
\right)
1401 \label{eqn:NH-chain-eqn-of-motion
}
1403 The conserved quantity for Nos
{\'e
}-Hoover chains is
1405 H =
\sum_{i=
1}^
{N
} \frac{\pb_i}{2m_i
} + U
\left(
\rv_1,
\rv_2,
\ldots,
\rv_N\right) +
\sum_{k=
1}^M
\frac{p^
2_
{\xi_k}}{2Q^
{\prime}_k
} + N_fkT
\xi_1 + kT
\sum_{k=
2}^M
\xi_k
1407 The values and velocities of the Nos
{\'e
}-Hoover thermostat variables
1408 are generally not included in the output, as they take up a fair
1409 amount of space and are generally not important for analysis of
1410 simulations, but this can be overridden by defining the environment
1411 variable
{\tt GMX_NOSEHOOVER_CHAINS
}, which will print the values of all
1412 the positions and velocities of all Nos
{\'e
}-Hoover particles in the
1413 chain to the
{\tt .edr
} file. Leap-frog simulations currently can only have
1414 Nos
{\'e
}-Hoover chain lengths of
1, but this will likely be updated in
1417 As described in the integrator section, for temperature coupling, the
1418 temperature that the algorithm attempts to match to the reference
1419 temperature is calculated differently in velocity Verlet and leap-frog
1420 dynamics. Velocity Verlet (
{\em md-vv
}) uses the full-step kinetic
1421 energy, while leap-frog and
{\em md-vv-avek
} use the half-step-averaged
1424 We can examine the Trotter decomposition again to better understand
1425 the differences between these constant-temperature integrators. In
1426 the case of Nos
{\'e
}-Hoover dynamics (for simplicity, using a chain
1427 with $N=
1$, with more details in Ref.~
\cite{Martyna1996
}), we split
1428 the Liouville operator as
1430 iL = iL_1 + iL_2 + iL_
{\mathrm{NHC
}},
1434 iL_1 &=&
\sum_{i=
1}^N
\left[\frac{\pb_i}{m_i
}\right]\cdot \frac{\partial}{\partial \rv_i} \nonumber \\
1435 iL_2 &=&
\sum_{i=
1}^N
\F_i\cdot \frac{\partial}{\partial \pb_i} \nonumber \\
1436 iL_
{\mathrm{NHC
}} &=&
\sum_{i=
1}^N-
\frac{p_
{\xi}}{Q
}\vv_i\cdot \nabla_{\vv_i} +
\frac{p_
{\xi}}{Q
}\frac{\partial }{\partial \xi} +
\left( T - T_0
\right)
\frac{\partial }{\partial p_
{\xi}}
1438 For standard velocity Verlet with Nos
{\'e
}-Hoover temperature control, this becomes
1440 \exp(iL
\dt) &=&
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right)
\exp\left(iL_2
\dt/
2\right)
\nonumber \\
1441 &&
\exp\left(iL_1
\dt\right)
\exp\left(iL_2
\dt/
2\right)
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right) +
\mathcal{O
}(
\Dt^
3).
1443 For half-step-averaged temperature control using
{\em md-vv-avek
},
1444 this decomposition will not work, since we do not have the full step
1445 temperature until after the second velocity step. However, we can
1446 construct an alternate decomposition that is still reversible, by
1447 switching the place of the NHC and velocity portions of the
1450 \exp(iL
\dt) &=&
\exp\left(iL_2
\dt/
2\right)
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right)
\exp\left(iL_1
\dt\right)
\nonumber \\
1451 &&
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right)
\exp\left(iL_2
\dt/
2\right)+
\mathcal{O
}(
\Dt^
3)
1452 \label{eq:half_step_NHC_integrator
}
1454 This formalism allows us to easily see the difference between the
1455 different flavors of velocity Verlet integrator. The leap-frog
1456 integrator can be seen as starting with
1457 Eq.~
\ref{eq:half_step_NHC_integrator
} just before the $
\exp\left(iL_1
1458 \dt\right)$ term, yielding:
1460 \exp(iL
\dt) &=&
\exp\left(iL_1
\dt\right)
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right)
\nonumber \\
1461 &&
\exp\left(iL_2
\dt\right)
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right) +
\mathcal{O
}(
\Dt^
3)
1463 and then using some algebra tricks to solve for some quantities are
1464 required before they are actually calculated~
\cite{Holian95
}.
1468 \subsubsection{Group temperature coupling
}\index{temperature-coupling group
}%
1469 In
{\gromacs} temperature coupling can be performed on groups of
1470 atoms, typically a protein and solvent. The reason such algorithms
1471 were introduced is that energy exchange between different components
1472 is not perfect, due to different effects including cut-offs etc. If
1473 now the whole system is coupled to one heat bath, water (which
1474 experiences the largest cut-off noise) will tend to heat up and the
1475 protein will cool down. Typically
100 K differences can be obtained.
1476 With the use of proper electrostatic methods (PME) these difference
1477 are much smaller but still not negligible. The parameters for
1478 temperature coupling in groups are given in the
{\tt mdp
} file.
1479 Recent investigation has shown that small temperature differences
1480 between protein and water may actually be an artifact of the way
1481 temperature is calculated when there are finite timesteps, and very
1482 large differences in temperature are likely a sign of something else
1483 seriously going wrong with the system, and should be investigated
1484 carefully~
\cite{Eastwood2010
}.
1486 One special case should be mentioned: it is possible to temperature-couple only
1487 part of the system, leaving other parts without temperature
1488 coupling. This is done by specifying $
{-
1}$ for the time constant
1489 $
\tau_T$ for the group that should not be thermostatted. If only
1490 part of the system is thermostatted, the system will still eventually
1491 converge to an NVT system. In fact, one suggestion for minimizing
1492 errors in the temperature caused by discretized timesteps is that if
1493 constraints on the water are used, then only the water degrees of
1494 freedom should be thermostatted, not protein degrees of freedom, as
1495 the higher frequency modes in the protein can cause larger deviations
1496 from the ``true'' temperature, the temperature obtained with small
1497 timesteps~
\cite{Eastwood2010
}.
1499 \subsection{Pressure coupling
\index{pressure coupling
}}
1500 In the same spirit as the temperature coupling, the system can also be
1501 coupled to a ``pressure bath.''
{\gromacs} supports both the Berendsen
1502 algorithm~
\cite{Berendsen84
} that scales coordinates and box vectors
1503 every step, the extended-ensemble Parrinello-Rahman approach~
\cite{Parrinello81,Nose83
}, and for
1504 the velocity Verlet variants, the Martyna-Tuckerman-Tobias-Klein
1505 (MTTK) implementation of pressure
1506 control~
\cite{Martyna1996
}. Parrinello-Rahman and Berendsen can be
1507 combined with any of the temperature coupling methods above; MTTK can
1508 only be used with Nos
{\'e
}-Hoover temperature control.
1510 \subsubsection{Berendsen pressure coupling
\pawsindexquiet{Berendsen
}{pressure coupling
}\index{weak coupling
}}
1511 \label{sec:berendsen_pressure_coupling
}
1512 The Berendsen algorithm rescales the
1513 coordinates and box vectors every step, or every $n_
\mathrm{PC
}$ steps,
1514 with a matrix
{\boldmath $
\mu$
},
1515 which has the effect of a first-order kinetic relaxation of the pressure
1516 towards a given reference pressure $
{\bf P
}_0$ according to
1518 \frac{\de {\bf P
}}{\de t
} =
\frac{{\bf P
}_0-
{\bf P
}}{\tau_p}.
1520 The scaling matrix
{\boldmath $
\mu$
} is given by
1523 =
\delta_{ij
} -
\frac{n_
\mathrm{PC
}\Delta t
}{3\,
\tau_p} \beta_{ij
} \
{P_
{0ij
} - P_
{ij
}(t) \
}.
1526 \index{isothermal compressibility
}
1527 \index{compressibility
}
1528 Here,
{\boldmath $
\beta$
} is the isothermal compressibility of the system.
1529 In most cases this will be a diagonal matrix, with equal elements on the
1530 diagonal, the value of which is generally not known.
1531 It suffices to take a rough estimate because the value of
{\boldmath $
\beta$
}
1532 only influences the non-critical time constant of the
1533 pressure relaxation without affecting the average pressure itself.
1534 For water at
1 atm and
300 K
1535 $
\beta =
4.6 \times 10^
{-
10}$ Pa$^
{-
1} =
4.6 \times 10^
{-
5}$ bar$^
{-
1}$,
1536 which is $
7.6 \times 10^
{-
4}$ MD units (see
\chref{defunits
}).
1537 Most other liquids have similar values.
1538 When scaling completely anisotropically, the system has to be rotated in
1539 order to obey
\eqnref{box_rot
}.
1540 This rotation is approximated in first order in the scaling, which is usually
1541 less than $
10^
{-
4}$. The actual scaling matrix
{\boldmath $
\mu'$
} is
1543 \mbox{\boldmath $
\mu'$
} =
1544 \left(
\begin{array
}{ccc
}
1545 \mu_{xx
} &
\mu_{xy
} +
\mu_{yx
} &
\mu_{xz
} +
\mu_{zx
} \\
1546 0 &
\mu_{yy
} &
\mu_{yz
} +
\mu_{zy
} \\
1550 The velocities are neither scaled nor rotated.
1552 In
{\gromacs}, the Berendsen scaling can also be done isotropically,
1553 which means that instead of $
\ve{P
}$ a diagonal matrix with elements of size
1554 trace$(
\ve{P
})/
3$ is used. For systems with interfaces, semi-isotropic
1555 scaling can be useful.
1556 In this case, the $x/y$-directions are scaled isotropically and the $z$
1557 direction is scaled independently. The compressibility in the $x/y$ or
1558 $z$-direction can be set to zero, to scale only in the other direction(s).
1560 If you allow full anisotropic deformations and use constraints you
1561 might have to scale more slowly or decrease your timestep to avoid
1562 errors from the constraint algorithms. It is important to note that
1563 although the Berendsen pressure control algorithm yields a simulation
1564 with the correct average pressure, it does not yield the exact NPT
1565 ensemble, and it is not yet clear exactly what errors this approximation
1568 % \ifthenelse{\equal{\gmxlite}{1}}{}{
1569 \subsubsection{Parrinello-Rahman pressure coupling
\pawsindexquiet{Parrinello-Rahman
}{pressure coupling
}}
1571 In cases where the fluctuations in pressure or volume are important
1572 {\em per se
} (
{\eg} to calculate thermodynamic properties), especially
1573 for small systems, it may be a problem that the exact ensemble is not
1574 well defined for the weak-coupling scheme, and that it does not
1575 simulate the true NPT ensemble.
1577 {\gromacs} also supports constant-pressure simulations using the
1578 Parrinello-Rahman approach~
\cite{Parrinello81,Nose83
}, which is similar
1579 to the Nos
{\'e
}-Hoover temperature coupling, and in theory gives the
1580 true NPT ensemble. With the Parrinello-Rahman barostat, the box
1581 vectors as represented by the matrix
\ve{b
} obey the matrix equation
1582 of motion
\footnote{The box matrix representation
\ve{b
} in
{\gromacs}
1583 corresponds to the transpose of the box matrix representation
\ve{h
}
1584 in the paper by Nos
{\'e
} and Klein. Because of this, some of our
1585 equations will look slightly different.
}
1587 \frac{\de \ve{b
}^
2}{\de t^
2}= V
\ve{W
}^
{-
1} \ve{b
}'^
{-
1} \left(
\ve{P
} -
\ve{P
}_
{ref
}\right).
1590 The volume of the box is denoted $V$, and $
\ve{W
}$ is a matrix parameter that determines
1591 the strength of the coupling. The matrices
\ve{P
} and
\ve{P
}$_
{ref
}$ are the
1592 current and reference pressures, respectively.
1594 The equations of motion for the particles are also changed, just as
1595 for the Nos
{\'e
}-Hoover coupling. In most cases you would combine the
1596 Parrinello-Rahman barostat with the Nos
{\'e
}-Hoover
1597 thermostat, but to keep it simple we only show the Parrinello-Rahman
1600 \bea \frac {\de^
2\ve{r
}_i
}{\de t^
2} & = &
\frac{\ve{F
}_i
}{m_i
} -
1601 \ve{M
} \frac{\de \ve{r
}_i
}{\de t
} , \\
\ve{M
} & = &
\ve{b
}^
{-
1} \left[
1602 \ve{b
} \frac{\de \ve{b
}'
}{\de t
} +
\frac{\de \ve{b
}}{\de t
} \ve{b
}'
1603 \right] \ve{b
}'^
{-
1}.
\eea The (inverse) mass parameter matrix
1604 $
\ve{W
}^
{-
1}$ determines the strength of the coupling, and how the box
1605 can be deformed. The box restriction (
\ref{eqn:box_rot
}) will be
1606 fulfilled automatically if the corresponding elements of $
\ve{W
}^
{-
1}$
1607 are zero. Since the coupling strength also depends on the size of your
1608 box, we prefer to calculate it automatically in
{\gromacs}. You only
1609 have to provide the approximate isothermal compressibilities
1610 {\boldmath $
\beta$
} and the pressure time constant $
\tau_p$ in the
1611 input file ($L$ is the largest box matrix element):
\beq \left(
1612 \ve{W
}^
{-
1} \right)_
{ij
} =
\frac{4 \pi^
2 \beta_{ij
}}{3 \tau_p^
2 L
}.
1613 \eeq Just as for the Nos
{\'e
}-Hoover thermostat, you should realize
1614 that the Parrinello-Rahman time constant is
{\em not
} equivalent to
1615 the relaxation time used in the Berendsen pressure coupling algorithm.
1616 In most cases you will need to use a
4--
5 times larger time constant
1617 with Parrinello-Rahman coupling. If your pressure is very far from
1618 equilibrium, the Parrinello-Rahman coupling may result in very large
1619 box oscillations that could even crash your run. In that case you
1620 would have to increase the time constant, or (better) use the weak-coupling
1621 scheme to reach the target pressure, and then switch to
1622 Parrinello-Rahman coupling once the system is in equilibrium.
1623 Additionally, using the leap-frog algorithm, the pressure at time $t$
1624 is not available until after the time step has completed, and so the
1625 pressure from the previous step must be used, which makes the algorithm
1626 not directly reversible, and may not be appropriate for high precision
1627 thermodynamic calculations.
1629 \subsubsection{Surface-tension coupling
\pawsindexquiet{surface-tension
}{pressure coupling
}}
1630 When a periodic system consists of more than one phase, separated by
1631 surfaces which are parallel to the $xy$-plane,
1632 the surface tension and the $z$-component of the pressure can be coupled
1633 to a pressure bath. Presently, this only works with the Berendsen
1634 pressure coupling algorithm in
{\gromacs}.
1635 The average surface tension $
\gamma(t)$ can be calculated from
1636 the difference between the normal and the lateral pressure
1639 \frac{1}{n
} \int_0^
{L_z
}
1640 \left\
{ P_
{zz
}(z,t) -
\frac{P_
{xx
}(z,t) + P_
{yy
}(z,t)
}{2} \right\
} \mbox{d
}z \\
1642 \frac{L_z
}{n
} \left\
{ P_
{zz
}(t) -
\frac{P_
{xx
}(t) + P_
{yy
}(t)
}{2} \right\
},
1644 where $L_z$ is the height of the box and $n$ is the number of surfaces.
1645 The pressure in the z-direction is corrected by scaling the height of
1646 the box with $
\mu_z$
1648 \Delta P_
{zz
} =
\frac{\Delta t
}{\tau_p} \
{ P_
{0zz
} - P_
{zz
}(t) \
}
1651 \mu_{zz
} =
1 +
\beta_{zz
} \Delta P_
{zz
}
1653 This is similar to normal pressure coupling, except that the power
1654 of $
1/
3$ is missing.
1655 The pressure correction in the $z$-direction is then used to get the
1656 correct convergence for the surface tension to the reference value $
\gamma_0$.
1657 The correction factor for the box length in the $x$/$y$-direction is
1659 \mu_{x/y
} =
1 +
\frac{\Delta t
}{2\,
\tau_p} \beta_{x/y
}
1660 \left(
\frac{n
\gamma_0}{\mu_{zz
} L_z
}
1661 -
\left\
{ P_
{zz
}(t)+
\Delta P_
{zz
} -
\frac{P_
{xx
}(t) + P_
{yy
}(t)
}{2} \right\
}
1664 The value of $
\beta_{zz
}$ is more critical than with normal pressure
1665 coupling. Normally an incorrect compressibility will just scale $
\tau_p$,
1666 but with surface tension coupling it affects the convergence of the surface
1668 When $
\beta_{zz
}$ is set to zero (constant box height), $
\Delta P_z$ is also set
1669 to zero, which is necessary for obtaining the correct surface tension.
1671 \subsubsection{MTTK pressure control algorithms
}
1673 As mentioned in the previous section, one weakness of leap-frog
1674 integration is in constant pressure simulations, since the pressure
1675 requires a calculation of both the virial and the kinetic energy at
1676 the full time step; for leap-frog, this information is not available
1677 until
{\em after
} the full timestep. Velocity Verlet does allow the
1678 calculation, at the cost of an extra round of global communication,
1679 and can compute, mod any integration errors, the true NPT ensemble.
1681 The full equations, combining both pressure coupling and temperature
1682 coupling, are taken from Martyna
{\em et al.
}~
\cite{Martyna1996
} and
1683 Tuckerman~
\cite{Tuckerman2006
} and are referred to here as MTTK
1684 equations (Martyna-Tuckerman-Tobias-Klein). We introduce for
1685 convenience $
\epsilon = (
1/
3)
\ln (V/V_0)$, where $V_0$ is a reference
1686 volume. The momentum of $
\epsilon$ is $
\veps = p_
{\epsilon}/W =
1687 \dot{\epsilon} =
\dot{V
}/
3V$, and define $
\alpha =
1 +
3/N_
{dof
}$ (see
1688 Ref~
\cite{Tuckerman2006
})
1690 The isobaric equations are
1692 \dot{\rv}_i &=&
\frac{\pb_i}{m_i
} +
\frac{\peps}{W
} \rv_i \nonumber \\
1693 \frac{\dot{\pb}_i
}{m_i
} &=&
\frac{1}{m_i
}\F_i -
\alpha\frac{\peps}{W
} \frac{\pb_i}{m_i
} \nonumber \\
1694 \dot{\epsilon} &=&
\frac{\peps}{W
} \nonumber \\
1695 \frac{\dot{\peps}}{W
} &=&
\frac{3V
}{W
}(P_
{\mathrm{int
}} - P) + (
\alpha-
1)
\left(
\sum_{n=
1}^N
\frac{\pb_i^
2}{m_i
}\right),\\
1699 P_
{\mathrm{int
}} &=& P_
{\mathrm{kin
}} -P_
{\mathrm{vir
}} =
\frac{1}{3V
}\left[\sum_{i=
1}^N
\left(
\frac{\pb_i^
2}{2m_i
} -
\rv_i \cdot \F_i\
1702 The terms including $
\alpha$ are required to make phase space
1703 incompressible~
\cite{Tuckerman2006
}. The $
\epsilon$ acceleration term
1706 \frac{\dot{\peps}}{W
} &=&
\frac{3V
}{W
}\left(
\alpha P_
{\mathrm{kin
}} - P_
{\mathrm{vir
}} - P
\right)
1708 In terms of velocities, these equations become
1710 \dot{\rv}_i &=&
\vv_i +
\veps \rv_i \nonumber \\
1711 \dot{\vv}_i &=&
\frac{1}{m_i
}\F_i -
\alpha\veps \vv_i \nonumber \\
1712 \dot{\epsilon} &=&
\veps \nonumber \\
1713 \dot{\veps} &=&
\frac{3V
}{W
}(P_
{\mathrm{int
}} - P) + (
\alpha-
1)
\left(
\sum_{n=
1}^N
\frac{1}{2} m_i
\vv_i^
2\right)
\nonumber \\
1714 P_
{\mathrm{int
}} &=& P_
{\mathrm{kin
}} - P_
{\mathrm{vir
}} =
\frac{1}{3V
}\left[\sum_{i=
1}^N
\left(
\frac{1}{2} m_i
\vv_i^
2 -
\rv_i \cdot \F_i\right)
\right]
1716 For these equations, the conserved quantity is
1718 H =
\sum_{i=
1}^
{N
} \frac{\pb_i^
2}{2m_i
} + U
\left(
\rv_1,
\rv_2,
\ldots,
\rv_N\right) +
\frac{p_
\epsilon}{2W
} + PV
1720 The next step is to add temperature control. Adding Nos
{\'e
}-Hoover
1721 chains, including to the barostat degree of freedom, where we use
1722 $
\eta$ for the barostat Nos
{\'e
}-Hoover variables, and $Q^
{\prime}$
1723 for the coupling constants of the thermostats of the barostats, we get
1725 \dot{\rv}_i &=&
\frac{\pb_i}{m_i
} +
\frac{\peps}{W
} \rv_i \nonumber \\
1726 \frac{\dot{\pb}_i
}{m_i
} &=&
\frac{1}{m_i
}\F_i -
\alpha\frac{\peps}{W
} \frac{\pb_i}{m_i
} -
\frac{p_
{\xi_1}}{Q_1
}\frac{\pb_i}{m_i
}\nonumber \\
1727 \dot{\epsilon} &=&
\frac{\peps}{W
} \nonumber \\
1728 \frac{\dot{\peps}}{W
} &=&
\frac{3V
}{W
}(
\alpha P_
{\mathrm{kin
}} - P_
{\mathrm{vir
}} - P) -
\frac{p_
{\eta_1}}{Q^
{\prime}_1
}\peps \nonumber \\
1729 \dot{\xi}_k &=&
\frac{p_
{\xi_k}}{Q_k
} \nonumber \\
1730 \dot{\eta}_k &=&
\frac{p_
{\eta_k}}{Q^
{\prime}_k
} \nonumber \\
1731 \dot{p
}_
{\xi_k} &=& G_k -
\frac{p_
{\xi_{k+
1}}}{Q_
{k+
1}} \;\;\;\; k=
1,
\ldots, M-
1 \nonumber \\
1732 \dot{p
}_
{\eta_k} &=& G^
\prime_k -
\frac{p_
{\eta_{k+
1}}}{Q^
\prime_{k+
1}} \;\;\;\; k=
1,
\ldots, M-
1 \nonumber \\
1733 \dot{p
}_
{\xi_M} &=& G_M
\nonumber \\
1734 \dot{p
}_
{\eta_M} &=& G^
\prime_M,
\nonumber \\
1738 P_
{\mathrm{int
}} &=& P_
{\mathrm{kin
}} - P_
{\mathrm{vir
}} =
\frac{1}{3V
}\left[\sum_{i=
1}^N
\left(
\frac{\pb_i^
2}{2m_i
} -
\rv_i \cdot \F_i\right)
\right] \nonumber \\
1739 G_1 &=&
\sum_{i=
1}^N
\frac{\pb^
2_i
}{m_i
} - N_f kT
\nonumber \\
1740 G_k &=&
\frac{p^
2_
{\xi_{k-
1}}}{2Q_
{k-
1}} - kT \;\; k =
2,
\ldots,M
\nonumber \\
1741 G^
\prime_1 &=&
\frac{\peps^
2}{2W
} - kT
\nonumber \\
1742 G^
\prime_k &=&
\frac{p^
2_
{\eta_{k-
1}}}{2Q^
\prime_{k-
1}} - kT \;\; k =
2,
\ldots,M
1744 The conserved quantity is now
1746 H =
\sum_{i=
1}^
{N
} \frac{\pb_i}{2m_i
} + U
\left(
\rv_1,
\rv_2,
\ldots,
\rv_N\right) +
\frac{p^
2_
\epsilon}{2W
} + PV +
\nonumber \\
1747 \sum_{k=
1}^M
\frac{p^
2_
{\xi_k}}{2Q_k
} +
\sum_{k=
1}^M
\frac{p^
2_
{\eta_k}}{2Q^
{\prime}_k
} + N_fkT
\xi_1 + kT
\sum_{i=
2}^M
\xi_k + kT
\sum_{k=
1}^M
\eta_k
1749 Returning to the Trotter decomposition formalism, for pressure control and temperature control~
\cite{Martyna1996
} we get:
1751 iL = iL_1 + iL_2 + iL_
{\epsilon,
1} + iL_
{\epsilon,
2} + iL_
{\mathrm{NHC-baro
}} + iL_
{\mathrm{NHC
}}
1753 where ``NHC-baro'' corresponds to the Nos
{\`e
}-Hoover chain of the barostat,
1754 and NHC corresponds to the NHC of the particles,
1756 iL_1 &=&
\sum_{i=
1}^N
\left[\frac{\pb_i}{m_i
} +
\frac{\peps}{W
}\rv_i\right]\cdot \frac{\partial}{\partial \rv_i} \\
1757 iL_2 &=&
\sum_{i=
1}^N
\F_i -
\alpha \frac{\peps}{W
}\pb_i \cdot \frac{\partial}{\partial \pb_i} \\
1758 iL_
{\epsilon,
1} &=&
\frac{p_
\epsilon}{W
} \frac{\partial}{\partial \epsilon}\\
1759 iL_
{\epsilon,
2} &=& G_
{\epsilon} \frac{\partial}{\partial p_
\epsilon}
1763 G_
{\epsilon} =
3V
\left(
\alpha P_
{\mathrm{kin
}} - P_
{\mathrm{vir
}} - P
\right)
1765 Using the Trotter decomposition, we get
1767 \exp(iL
\dt) &=&
\exp\left(iL_
{\mathrm{NHC-baro
}}\dt/
2\right)
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right)
\nonumber \nonumber \\
1768 &&
\exp\left(iL_
{\epsilon,
2}\dt/
2\right)
\exp\left(iL_2
\dt/
2\right)
\nonumber \nonumber \\
1769 &&
\exp\left(iL_
{\epsilon,
1}\dt\right)
\exp\left(iL_1
\dt\right)
\nonumber \nonumber \\
1770 &&
\exp\left(iL_2
\dt/
2\right)
\exp\left(iL_
{\epsilon,
2}\dt/
2\right)
\nonumber \nonumber \\
1771 &&
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right)
\exp\left(iL_
{\mathrm{NHC-baro
}}\dt/
2\right) +
\mathcal{O
}(
\dt^
3)
1773 The action of $
\exp\left(iL_1
\dt\right)$ comes from the solution of
1774 the the differential equation
1775 $
\dot{\rv}_i =
\vv_i +
\veps \rv_i$
1776 with $
\vv_i =
\pb_i/m_i$ and $
\veps$ constant with initial condition
1777 $
\rv_i(
0)$, evaluate at $t=
\Delta t$. This yields the evolution
1779 \rv_i(
\dt) =
\rv_i(
0)e^
{\veps \dt} +
\Delta t
\vv_i(
0) e^
{\veps \dt/
2} \sinhx{\veps \dt/
2}.
1781 The action of $
\exp\left(iL_2
\dt/
2\right)$ comes from the solution
1782 of the differential equation $
\dot{\vv}_i =
\frac{\F_i}{m_i
} -
1783 \alpha\veps\vv_i$, yielding
1785 \vv_i(
\dt/
2) =
\vv_i(
0)e^
{-
\alpha\veps \dt/
2} +
\frac{\Delta t
}{2m_i
}\F_i(
0) e^
{-
\alpha\veps \dt/
4}\sinhx{\alpha\veps \dt/
4}.
1787 {\em md-vv-avek
} uses the full step kinetic energies for determining the pressure with the pressure control,
1788 but the half-step-averaged kinetic energy for the temperatures, which can be written as a Trotter decomposition as
1790 \exp(iL
\dt) &=&
\exp\left(iL_
{\mathrm{NHC-baro
}}\dt/
2\right)
\nonumber \exp\left(iL_
{\epsilon,
2}\dt/
2\right)
\exp\left(iL_2
\dt/
2\right)
\nonumber \\
1791 &&
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right)
\exp\left(iL_
{\epsilon,
1}\dt\right)
\exp\left(iL_1
\dt\right)
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right)
\nonumber \\
1792 &&
\exp\left(iL_2
\dt/
2\right)
\exp\left(iL_
{\epsilon,
2}\dt/
2\right)
\exp\left(iL_
{\mathrm{NHC-baro
}}\dt/
2\right) +
\mathcal{O
}(
\dt^
3)
1794 With constraints, the equations become significantly more
1795 complicated, in that each of these equations need to be solved
1796 iteratively for the constraint forces. The discussion of the details of the iteration
1797 is beyond the scope of this manual; readers are encouraged to see the
1798 implementation described in~
\cite{Yu2010
}.
1801 \subsubsection{Infrequent evaluation of temperature and pressure coupling
}
1803 Temperature and pressure control require global communication to
1804 compute the kinetic energy and virial, which can become costly if
1805 performed every step for large systems. We can rearrange the Trotter
1806 decomposition to give alternate symplectic, reversible integrator with
1807 the coupling steps every $n$ steps instead of every steps. These new
1808 integrators will diverge if the coupling time step is too large, as
1809 the auxiliary variable integrations will not converge. However, in
1810 most cases, long coupling times are more appropriate, as they disturb
1811 the dynamics less~
\cite{Martyna1996
}.
1813 Standard velocity Verlet with Nos
{\'e
}-Hoover temperature control has a Trotter expansion
1815 \exp(iL
\dt) &
\approx&
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right)
\exp\left(iL_2
\dt/
2\right)
\nonumber \\
1816 &&
\exp\left(iL_1
\dt\right)
\exp\left(iL_2
\dt/
2\right)
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right).
1818 If the Nos
{\'e
}-Hoover chain is sufficiently slow with respect to the motions of the system, we can
1819 write an alternate integrator over $n$ steps for velocity Verlet as
1821 \exp(iL
\dt) &
\approx& (
\exp\left(iL_
{\mathrm{NHC
}}(n
\dt/
2)
\right)
\left[\exp\left(iL_2
\dt/
2\right)
\right.
\nonumber \\
1822 &&
\left.
\exp\left(iL_1
\dt\right)
\exp\left(iL_2
\dt/
2\right)
\right]^n
\exp\left(iL_
{\mathrm{NHC
}}(n
\dt/
2)
\right).
1824 For pressure control, this becomes
1826 \exp(iL
\dt) &
\approx&
\exp\left(iL_
{\mathrm{NHC-baro
}}(n
\dt/
2)
\right)
\exp\left(iL_
{\mathrm{NHC
}}(n
\dt/
2)
\right)
\nonumber \nonumber \\
1827 &&
\exp\left(iL_
{\epsilon,
2}(n
\dt/
2)
\right)
\left[\exp\left(iL_2
\dt/
2\right)
\right.
\nonumber \nonumber \\
1828 &&
\exp\left(iL_
{\epsilon,
1}\dt\right)
\exp\left(iL_1
\dt\right)
\nonumber \nonumber \\
1829 &&
\left.
\exp\left(iL_2
\dt/
2\right)
\right]^n
\exp\left(iL_
{\epsilon,
2}(n
\dt/
2)
\right)
\nonumber \nonumber \\
1830 &&
\exp\left(iL_
{\mathrm{NHC
}}(n
\dt/
2)
\right)
\exp\left(iL_
{\mathrm{NHC-baro
}}(n
\dt/
2)
\right),
1832 where the box volume integration occurs every step, but the auxiliary variable
1833 integrations happen every $n$ steps.
1835 % } % Brace matches ifthenelse test for gmxlite
1838 \subsection{The complete update algorithm
}
1841 \addtolength{\fboxsep}{0.5cm
}
1842 \begin{shadowenv
}[12cm
]
1843 {\large \bf THE UPDATE ALGORITHM
}
1844 \rule{\textwidth}{2pt
} \\
1846 Positions $
\ve{r
}$ of all atoms at time $t$ \\
1847 Velocities $
\ve{v
}$ of all atoms at time $t-
\hDt$ \\
1848 Accelerations $
\ve{F
}/m$ on all atoms at time $t$.\\
1849 (Forces are computed disregarding any constraints)\\
1850 Total kinetic energy and virial at $t-
\Dt$\\
1852 {\bf 1.
} Compute the scaling factors $
\lambda$ and $
\mu$\\
1853 according to
\eqnsref{lambda
}{mu
}\\
1855 {\bf 2.
} Update and scale velocities: $
\ve{v
}' =
\lambda (
\ve{v
} +
1856 \ve{a
} \Delta t)$ \\
1858 {\bf 3.
} Compute new unconstrained coordinates: $
\ve{r
}' =
\ve{r
} +
\ve{v
}'
1861 {\bf 4.
} Apply constraint algorithm to coordinates: constrain($
\ve{r
}^
{'
} \rightarrow \ve{r
}'';
1864 {\bf 5.
} Correct velocities for constraints: $
\ve{v
} = (
\ve{r
}'' -
1865 \ve{r
}) /
\Delta t$ \\
1867 {\bf 6.
} Scale coordinates and box: $
\ve{r
} =
\mu \ve{r
}'';
\ve{b
} =
1870 \caption{The MD update algorithm with the leap-frog integrator
}
1871 \label{fig:complete-update
}
1874 The complete algorithm for the update of velocities and coordinates is
1875 given using leap-frog in
\figref{complete-update
}. The SHAKE algorithm of step
1876 4 is explained below.
1878 {\gromacs} has a provision to ``freeze'' (prevent motion of) selected
1879 particles
\index{frozen atoms
}, which must be defined as a ``
\swapindex{freeze
}{group
}.'' This is implemented
1880 using a
{\em freeze factor $
\ve{f
}_g$
}, which is a vector, and differs for each
1881 freeze group (see
\secref{groupconcept
}). This vector contains only
1882 zero (freeze) or one (don't freeze).
1883 When we take this freeze factor and the external acceleration $
\ve{a
}_h$ into
1884 account the update algorithm for the velocities becomes
1886 \ve{v
}(t+
\hdt)~=~
\ve{f
}_g *
\lambda *
\left[ \ve{v
}(t-
\hdt) +
\frac{\ve{F
}(t)
}{m
}\Delta t +
\ve{a
}_h
\Delta t
\right],
1888 where $g$ and $h$ are group indices which differ per atom.
1890 \subsection{Output step
}
1891 The most important output of the MD run is the
{\em
1892 \swapindex{trajectory
}{file
}}, which contains particle coordinates
1893 and (optionally) velocities at regular intervals.
1894 The trajectory file contains frames that could include positions,
1895 velocities and/or forces, as well as information about the dimensions
1896 of the simulation volume, integration step, integration time, etc. The
1897 interpretation of the time varies with the integrator chosen, as
1898 described above. For velocity-Verlet integrators, velocities labeled
1899 at time $t$ are for that time. For other integrators (e.g. leap-frog,
1900 stochastic dynamics), the velocities labeled at time $t$ are for time
1903 Since the trajectory
1904 files are lengthy, one should not save every step! To retain all
1905 information it suffices to write a frame every
15 steps, since at
1906 least
30 steps are made per period of the highest frequency in the
1907 system, and Shannon's
\normindex{sampling
} theorem states that two samples per
1908 period of the highest frequency in a band-limited signal contain all
1909 available information. But that still gives very long files! So, if
1910 the highest frequencies are not of interest,
10 or
20 samples per ps
1911 may suffice. Be aware of the distortion of high-frequency motions by
1912 the
{\em stroboscopic effect
}, called
{\em aliasing
}: higher frequencies
1913 are mirrored with respect to the sampling frequency and appear as
1916 {\gromacs} can also write reduced-precision coordinates for a subset of
1917 the simulation system to a special compressed trajectory file
1918 format. All the other tools can read and write this format. See
1919 \secref{mdpopt
} for details on how to set up your
{\tt .mdp
} file
1920 to have
{\tt mdrun
} use this feature.
1922 % \ifthenelse{\equal{\gmxlite}{1}}{}{
1923 \section{Shell molecular dynamics
}
1924 {\gromacs} can simulate
\normindex{polarizability
} using the
1925 \normindex{shell model
} of Dick and Overhauser~
\cite{Dick58
}. In such models
1926 a shell particle representing the electronic degrees of freedom is
1927 attached to a nucleus by a spring. The potential energy is minimized with
1928 respect to the shell position at every step of the simulation (see below).
1929 Successful applications of shell models in
{\gromacs} have been published
1930 for $N_2$~
\cite{Jordan95
} and water~
\cite{Maaren2001a
}.
1932 \subsection{Optimization of the shell positions
}
1933 The force
\ve{F
}$_S$ on a shell particle $S$ can be decomposed into two
1936 \ve{F
}_S ~=~
\ve{F
}_
{bond
} +
\ve{F
}_
{nb
}
1938 where
\ve{F
}$_
{bond
}$ denotes the component representing the
1939 polarization energy, usually represented by a harmonic potential and
1940 \ve{F
}$_
{nb
}$ is the sum of Coulomb and van der Waals interactions. If we
1941 assume that
\ve{F
}$_
{nb
}$ is almost constant we can analytically derive the
1942 optimal position of the shell, i.e. where
\ve{F
}$_S$ =
0. If we have the
1943 shell S connected to atom A we have
1945 \ve{F
}_
{bond
} ~=~ k_b
\left(
\ve{x
}_S -
\ve{x
}_A
\right).
1947 In an iterative solver, we have positions
\ve{x
}$_S(n)$ where $n$ is
1948 the iteration count. We now have at iteration $n$
1950 \ve{F
}_
{nb
} ~=~
\ve{F
}_S - k_b
\left(
\ve{x
}_S(n) -
\ve{x
}_A
\right)
1952 and the optimal position for the shells $x_S(n+
1)$ thus follows from
1954 \ve{F
}_S - k_b
\left(
\ve{x
}_S(n) -
\ve{x
}_A
\right) + k_b
\left(
\ve{x
}_S(n+
1) -
\ve{x
}_A
\right) =
0
1958 \Delta \ve{x
}_S =
\ve{x
}_S(n+
1) -
\ve{x
}_S(n)
1962 \Delta \ve{x
}_S =
\ve{F
}_S/k_b
1964 which then yields the algorithm to compute the next trial in the optimization
1967 \ve{x
}_S(n+
1) ~=~
\ve{x
}_S(n) +
\ve{F
}_S/k_b.
1969 % } % Brace matches ifthenelse test for gmxlite
1971 \section{Constraint algorithms
\index{constraint algorithms
}}
1972 Constraints can be imposed in
{\gromacs} using LINCS (default) or
1973 the traditional SHAKE method.
1975 \subsection{\normindex{SHAKE
}}
1976 \label{subsec:SHAKE
}
1977 The SHAKE~
\cite{Ryckaert77
} algorithm changes a set of unconstrained
1978 coordinates $
\ve{r
}^
{'
}$ to a set of coordinates $
\ve{r
}''$ that
1979 fulfill a list of distance constraints, using a set $
\ve{r
}$
1982 {\rm SHAKE
}(
\ve{r
}^
{'
} \rightarrow \ve{r
}'';\,
\ve{r
})
1984 This action is consistent with solving a set of Lagrange multipliers
1985 in the constrained equations of motion. SHAKE needs a
{\em relative tolerance
};
1986 it will continue until all constraints are satisfied within
1987 that relative tolerance. An error message is
1988 given if SHAKE cannot reset the coordinates because the deviation is
1989 too large, or if a given number of iterations is surpassed.
1991 Assume the equations of motion must fulfill $K$ holonomic constraints,
1994 \sigma_k(
\ve{r
}_1
\ldots \ve{r
}_N) =
0; \;\; k=
1 \ldots K.
1996 For example, $(
\ve{r
}_1 -
\ve{r
}_2)^
2 - b^
2 =
0$.
1997 Then the forces are defined as
1999 -
\frac{\partial}{\partial \ve{r
}_i
} \left( V +
\sum_{k=
1}^K
\lambda_k
2002 where $
\lambda_k$ are Lagrange multipliers which must be solved to
2003 fulfill the constraint equations. The second part of this sum
2004 determines the
{\em constraint forces
} $
\ve{G
}_i$, defined by
2006 \ve{G
}_i = -
\sum_{k=
1}^K
\lambda_k \frac{\partial \sigma_k}{\partial
2009 The displacement due to the constraint forces in the leap-frog or
2010 Verlet algorithm is equal to $(
\ve{G
}_i/m_i)(
\Dt)^
2$. Solving the
2011 Lagrange multipliers (and hence the displacements) requires the
2012 solution of a set of coupled equations of the second degree. These are
2013 solved iteratively by SHAKE.
2014 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2015 \label{subsec:SETTLE
}
2016 For the special case of rigid water molecules, that often make up more
2017 than
80\% of the simulation system we have implemented the
2019 algorithm~
\cite{Miyamoto92
} (
\secref{constraints
}).
2021 For velocity Verlet, an additional round of constraining must be
2022 done, to constrain the velocities of the second velocity half step,
2023 removing any component of the velocity parallel to the bond vector.
2024 This step is called RATTLE, and is covered in more detail in the
2025 original Andersen paper~
\cite{Andersen1983a
}.
2027 % } % Brace matches ifthenelse test for gmxlite
2032 \newcommand{\fs}[1]{\begin{equation
} \label{eqn:
#1}}
2033 \newcommand{\fe}{\end{equation
}}
2034 \newcommand{\p}{\partial}
2035 \newcommand{\Bm}{\ve{B
}}
2036 \newcommand{\M}{\ve{M
}}
2037 \newcommand{\iM}{\M^
{-
1}}
2038 \newcommand{\Tm}{\ve{T
}}
2039 \newcommand{\Sm}{\ve{S
}}
2040 \newcommand{\fo}{\ve{f
}}
2041 \newcommand{\con}{\ve{g
}}
2042 \newcommand{\lenc}{\ve{d
}}
2044 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2045 \subsection{\normindex{LINCS
}}
2046 \label{subsec:lincs
}
2048 \subsubsection{The LINCS algorithm
}
2049 LINCS is an algorithm that resets bonds to their correct lengths
2050 after an unconstrained update~
\cite{Hess97
}.
2051 The method is non-iterative, as it always uses two steps.
2052 Although LINCS is based on matrices, no matrix-matrix multiplications are
2053 needed. The method is more stable and faster than SHAKE,
2054 but it can only be used with bond constraints and
2055 isolated angle constraints, such as the proton angle in OH.
2056 Because of its stability, LINCS is especially useful for Brownian dynamics.
2057 LINCS has two parameters, which are explained in the subsection parameters.
2058 The parallel version of LINCS, P-LINCS, is described
2059 in subsection
\ssecref{plincs
}.
2061 \subsubsection{The LINCS formulas
}
2062 We consider a system of $N$ particles, with positions given by a
2063 $
3N$ vector $
\ve{r
}(t)$.
2064 For molecular dynamics the equations of motion are given by Newton's Law
2066 {\de^
2 \ve{r
} \over \de t^
2} =
\iM \ve{F
},
2068 where $
\ve{F
}$ is the $
3N$ force vector
2069 and $
\M$ is a $
3N
\times 3N$ diagonal matrix,
2070 containing the masses of the particles.
2071 The system is constrained by $K$ time-independent constraint equations
2073 g_i(
\ve{r
}) = |
\ve{r
}_
{i_1
}-
\ve{r
}_
{i_2
} | - d_i =
0 ~~~~~~i=
1,
\ldots,K.
2076 In a numerical integration scheme, LINCS is applied after an
2077 unconstrained update, just like SHAKE. The algorithm works in two
2078 steps (see figure
\figref{lincs
}). In the first step, the projections
2079 of the new bonds on the old bonds are set to zero. In the second step,
2080 a correction is applied for the lengthening of the bonds due to
2081 rotation. The numerics for the first step and the second step are very
2082 similar. A complete derivation of the algorithm can be found in
2083 \cite{Hess97
}. Only a short description of the first step is given
2087 \centerline{\includegraphics[height=
50mm
]{plots/lincs
}}
2088 \caption[The three position updates needed for one time step.
]{The
2089 three position updates needed for one time step. The dashed line is
2090 the old bond of length $d$, the solid lines are the new bonds. $l=d
2091 \cos \theta$ and $p=(
2 d^
2 - l^
2)^
{1 \over 2}$.
}
2095 A new notation is introduced for the gradient matrix of the constraint
2096 equations which appears on the right hand side of this equation:
2098 B_
{hi
} =
{\p g_h
\over \p r_i
}
2100 Notice that $
\Bm$ is a $K
\times 3N$ matrix, it contains the directions
2102 The following equation shows how the new constrained coordinates
2103 $
\ve{r
}_
{n+
1}$ are related to the unconstrained coordinates
2104 $
\ve{r
}_
{n+
1}^
{unc
}$ by
2107 \ve{r
}_
{n+
1}=(
\ve{I
}-
\Tm_n \ve{B
}_n)
\ve{r
}_
{n+
1}^
{unc
} +
\Tm_n \lenc=
2109 \ve{r
}_
{n+
1}^
{unc
} -
2110 \iM \Bm_n (
\Bm_n \iM \Bm_n^T)^
{-
1} (
\Bm_n \ve{r
}_
{n+
1}^
{unc
} -
\lenc)
2113 where $
\Tm =
\iM \Bm^T (
\Bm \iM \Bm^T)^
{-
1}$.
2114 The derivation of this equation from
\eqnsref{c1
}{c2
} can be found
2117 This first step does not set the real bond lengths to the prescribed lengths,
2118 but the projection of the new bonds onto the old directions of the bonds.
2119 To correct for the rotation of bond $i$, the projection of the
2120 bond, $p_i$, on the old direction is set to
2122 p_i=
\sqrt{2 d_i^
2 - l_i^
2},
2124 where $l_i$ is the bond length after the first projection.
2125 The corrected positions are
2127 \ve{r
}_
{n+
1}^*=(
\ve{I
}-
\Tm_n \Bm_n)
\ve{r
}_
{n+
1} +
\Tm_n \ve{p
}.
2129 This correction for rotational effects is actually an iterative process,
2130 but during MD only one iteration is applied.
2131 The relative constraint deviation after this procedure will be less than
2132 0.0001 for every constraint.
2133 In energy minimization, this might not be accurate enough, so the number
2134 of iterations is equal to the order of the expansion (see below).
2136 Half of the CPU time goes to inverting the constraint coupling
2137 matrix $
\Bm_n \iM \Bm_n^T$, which has to be done every time step.
2138 This $K
\times K$ matrix
2139 has $
1/m_
{i_1
} +
1/m_
{i_2
}$ on the diagonal.
2140 The off-diagonal elements are only non-zero when two bonds are connected,
2142 $
\cos \phi /m_c$, where $m_c$ is
2143 the mass of the atom connecting the
2144 two bonds and $
\phi$ is the angle between the bonds.
2146 The matrix $
\Tm$ is inverted through a power expansion.
2147 A $K
\times K$ matrix $
\ve{S
}$ is
2148 introduced which is the inverse square root of
2149 the diagonal of $
\Bm_n \iM \Bm_n^T$.
2150 This matrix is used to convert the diagonal elements
2151 of the coupling matrix to one:
2154 (
\Bm_n \iM \Bm_n^T)^
{-
1}
2155 =
\Sm \Sm^
{-
1} (
\Bm_n \iM \Bm_n^T)^
{-
1} \Sm^
{-
1} \Sm \\
[2mm
]
2156 =
\Sm (
\Sm \Bm_n \iM \Bm_n^T
\Sm)^
{-
1} \Sm =
2157 \Sm (
\ve{I
} -
\ve{A
}_n)^
{-
1} \Sm
2160 The matrix $
\ve{A
}_n$ is symmetric and sparse and has zeros on the diagonal.
2161 Thus a simple trick can be used to calculate the inverse:
2163 (
\ve{I
}-
\ve{A
}_n)^
{-
1}=
2164 \ve{I
} +
\ve{A
}_n +
\ve{A
}_n^
2 +
\ve{A
}_n^
3 +
\ldots
2167 This inversion method is only valid if the absolute values of all the
2168 eigenvalues of $
\ve{A
}_n$ are smaller than one.
2169 In molecules with only bond constraints, the connectivity is so low
2170 that this will always be true, even if ring structures are present.
2171 Problems can arise in angle-constrained molecules.
2172 By constraining angles with additional distance constraints,
2173 multiple small ring structures are introduced.
2174 This gives a high connectivity, leading to large eigenvalues.
2175 Therefore LINCS should NOT be used with coupled angle-constraints.
2177 For molecules with all bonds constrained the eigenvalues of $A$
2178 are around
0.4. This means that with each additional order
2179 in the expansion
\eqnref{m3
} the deviations decrease by a factor
0.4.
2180 But for relatively isolated triangles of constraints the largest
2181 eigenvalue is around
0.7.
2182 Such triangles can occur when removing hydrogen angle vibrations
2183 with an additional angle constraint in alcohol groups
2184 or when constraining water molecules with LINCS, for instance
2185 with flexible constraints.
2186 The constraints in such triangles converge twice as slow as
2187 the other constraints. Therefore, starting with
{\gromacs} 4,
2188 additional terms are added to the expansion for such triangles
2190 (
\ve{I
}-
\ve{A
}_n)^
{-
1} \approx
2191 \ve{I
} +
\ve{A
}_n +
\ldots +
\ve{A
}_n^
{N_i
} +
2192 \left(
\ve{A
}^*_n +
\ldots +
{\ve{A
}_n^*
}^
{N_i
} \right)
\ve{A
}_n^
{N_i
}
2194 where $N_i$ is the normal order of the expansion and
2195 $
\ve{A
}^*$ only contains the elements of $
\ve{A
}$ that couple
2196 constraints within rigid triangles, all other elements are zero.
2197 In this manner, the accuracy of angle constraints comes close
2198 to that of the other constraints, while the series of matrix vector
2199 multiplications required for determining the expansion
2200 only needs to be extended for a few constraint couplings.
2201 This procedure is described in the P-LINCS paper
\cite{Hess2008a
}.
2203 \subsubsection{The LINCS Parameters
}
2204 The accuracy of LINCS depends on the number of matrices used
2205 in the expansion
\eqnref{m3
}. For MD calculations a fourth order
2206 expansion is enough. For Brownian dynamics with
2207 large time steps an eighth order expansion may be necessary.
2208 The order is a parameter in the
{\tt *.mdp
} file.
2209 The implementation of LINCS is done in such a way that the
2210 algorithm will never crash. Even when it is impossible to
2211 to reset the constraints LINCS will generate a conformation
2212 which fulfills the constraints as well as possible.
2213 However, LINCS will generate a warning when in one step a bond
2214 rotates over more than a predefined angle.
2215 This angle is set by the user in the
{\tt *.mdp
} file.
2217 % } % Brace matches ifthenelse test for gmxlite
2220 \section{Simulated Annealing
}
2222 The well known
\swapindex{simulated
}{annealing
}
2223 (SA) protocol is supported in
{\gromacs}, and you can even couple multiple
2224 groups of atoms separately with an arbitrary number of reference temperatures
2225 that change during the simulation. The annealing is implemented by simply
2226 changing the current reference temperature for each group in the temperature
2227 coupling, so the actual relaxation and coupling properties depends on the
2228 type of thermostat you use and how hard you are coupling it. Since we are
2229 changing the reference temperature it is important to remember that the system
2230 will NOT instantaneously reach this value - you need to allow for the inherent
2231 relaxation time in the coupling algorithm too. If you are changing the
2232 annealing reference temperature faster than the temperature relaxation you
2233 will probably end up with a crash when the difference becomes too large.
2235 The annealing protocol is specified as a series of corresponding times and
2236 reference temperatures for each group, and you can also choose whether you only
2237 want a single sequence (after which the temperature will be coupled to the
2238 last reference value), or if the annealing should be periodic and restart at
2239 the first reference point once the sequence is completed. You can mix and
2240 match both types of annealing and non-annealed groups in your simulation.
2242 \newcommand{\vrond}{\stackrel{\circ}{\ve{r
}}}
2243 \newcommand{\rond}{\stackrel{\circ}{r
}}
2244 \newcommand{\ruis}{\ve{r
}^G
}
2246 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2247 \section{Stochastic Dynamics
\swapindexquiet{stochastic
}{dynamics
}}
2249 Stochastic or velocity
\swapindex{Langevin
}{dynamics
} adds a friction
2250 and a noise term to Newton's equations of motion, as
2253 m_i
{\de^
2 \ve{r
}_i
\over \de t^
2} =
2254 - m_i
\gamma_i {\de \ve{r
}_i
\over \de t
} +
\ve{F
}_i(
\ve{r
}) +
\vrond_i,
2256 where $
\gamma_i$ is the friction constant $
[1/
\mbox{ps
}]$ and
2257 $
\vrond_i\!\!(t)$ is a noise process with
2258 $
\langle \rond_i\!\!(t)
\rond_j\!\!(t+s)
\rangle =
2259 2 m_i
\gamma_i k_B T
\delta(s)
\delta_{ij
}$.
2260 When $
1/
\gamma_i$ is large compared to the time scales present in the system,
2261 one could see stochastic dynamics as molecular dynamics with stochastic
2262 temperature-coupling. The advantage compared to MD with Berendsen
2263 temperature-coupling is that in case of SD the generated ensemble is known.
2264 For simulating a system in vacuum there is the additional advantage that there is no
2265 accumulation of errors for the overall translational and rotational
2267 When $
1/
\gamma_i$ is small compared to the time scales present in the system,
2268 the dynamics will be completely different from MD, but the sampling is
2271 In
{\gromacs} there are two algorithms to integrate equation (
\ref{SDeq
}):
2272 a simple and efficient one
2273 and a more complex leap-frog algorithm~
\cite{Gunsteren88
}.
2274 The accuracy of both integrators is equivalent to the normal MD leap-frog and
2275 velocity-Verlet integrator, except with constraints where the simple
2276 SD integrator is significantly less accurate. There is a proper way
2277 of applying constraints with the simple integrator, but that requires
2278 a second constraining step~
\cite{Goga2012
}, which diminishes the gain.
2279 The simple integrator is:
2282 \ve{v
}(t+
\hDt) &~=~&
\alpha \,
\ve{v
}(t-
\hDt) +
\frac{1 -
\alpha}{m
\gamma}\ve{F
}(t) +
\sqrt{\frac{k_B T
}{m
}(
1 -
\alpha^
2)
} \,
\ruis_i \\
2283 \ve{r
}(t+
\Dt) &~=~&
\ve{r
}(t)+
\Dt \,
\ve{v
}(t+
\hDt) \\
2284 \alpha &~=~&
\left(
1 -
\frac{\gamma \Dt}{m
} \right)
2286 where $
\ruis_i$ is Gaussian distributed noise with $
\mu =
0$, $
\sigma =
1$.
2287 With constraints you should only consider using the simple integrator when $
\gamma \Dt/m
\ll 0.01$.
2289 In the complex algorithm four Gaussian random numbers are required
2290 per integration step per degree of freedom, and with constraints the
2291 coordinates need to be constrained twice per integration step.
2292 Depending on the computational cost of the force calculation,
2293 this can take a significant part of the simulation time.
2294 Exact continuation of a stochastic dynamics simulation is not possible,
2295 because the state of the random number generator is not stored.
2296 When using SD as a thermostat, an appropriate value for $
\gamma$ is
0.5 ps$^
{-
1}$,
2297 since this results in a friction that is lower than the internal friction
2298 of water, while it is high enough to remove excess heat
2299 (unless plain cut-off or reaction-field electrostatics is used).
2300 With this value of $
\gamma$ the efficient algorithm will usually be accurate
2303 \section{Brownian Dynamics
\swapindexquiet{Brownian
}{dynamics
}}
2305 In the limit of high friction, stochastic dynamics reduces to
2306 Brownian dynamics, also called position Langevin dynamics.
2307 This applies to over-damped systems,
2308 {\ie} systems in which the inertia effects are negligible.
2311 {\de \ve{r
}_i
\over \de t
} =
\frac{1}{\gamma_i} \ve{F
}_i(
\ve{r
}) +
\vrond_i
2313 where $
\gamma_i$ is the friction coefficient $
[\mbox{amu/ps
}]$ and
2314 $
\vrond_i\!\!(t)$ is a noise process with
2315 $
\langle \rond_i\!\!(t)
\rond_j\!\!(t+s)
\rangle =
2316 2 \delta(s)
\delta_{ij
} k_B T /
\gamma_i$.
2317 In
{\gromacs} the equations are integrated with a simple, explicit scheme
2319 \ve{r
}_i(t+
\Delta t) =
\ve{r
}_i(t) +
2320 {\Delta t
\over \gamma_i} \ve{F
}_i(
\ve{r
}(t))
2321 +
\sqrt{2 k_B T
{\Delta t
\over \gamma_i}}\,
\ruis_i,
2323 where $
\ruis_i$ is Gaussian distributed noise with $
\mu =
0$, $
\sigma =
1$.
2324 The friction coefficients $
\gamma_i$ can be chosen the same for all
2325 particles or as $
\gamma_i = m_i\,
\gamma_i$, where the friction constants
2326 $
\gamma_i$ can be different for different groups of atoms.
2327 Because the system is assumed to be over-damped, large timesteps
2328 can be used. LINCS should be used for the constraints since SHAKE
2329 will not converge for large atomic displacements.
2330 BD is an option of the
{\tt mdrun
} program.
2331 % } % Brace matches ifthenelse test for gmxlite
2333 \section{Energy Minimization
}
2334 \label{sec:EM
}\index{energy minimization
}%
2335 Energy minimization in
{\gromacs} can be done using steepest descent,
2336 conjugate gradients, or l-bfgs (limited-memory
2337 Broyden-Fletcher-Goldfarb-Shanno quasi-Newtonian minimizer...we
2338 prefer the abbreviation). EM is just an option of the
{\tt mdrun
}
2341 \subsection{Steepest Descent
\index{steepest descent
}}
2342 Although steepest descent is certainly not the most efficient
2343 algorithm for searching, it is robust and easy to implement.
2345 We define the vector $
\ve{r
}$ as the vector of all $
3N$ coordinates.
2346 Initially a maximum displacement $h_0$ (
{\eg} 0.01 nm) must be given.
2348 First the forces $
\ve{F
}$ and potential energy are calculated.
2349 New positions are calculated by
2351 \ve{r
}_
{n+
1} =
\ve{r
}_n +
\frac{\ve{F
}_n
}{\max (|
\ve{F
}_n|)
} h_n,
2353 where $h_n$ is the maximum displacement and $
\ve{F
}_n$ is the force,
2354 or the negative gradient of the potential $V$. The notation $
\max
2355 (|
\ve{F
}_n|)$ means the largest of the absolute values of the force
2356 components. The forces and energy are again computed for the new positions \\
2357 If ($V_
{n+
1} < V_n$) the new positions are accepted and $h_
{n+
1} =
1.2
2359 If ($V_
{n+
1} \geq V_n$) the new positions are rejected and $h_n =
0.2 h_n$.
2361 The algorithm stops when either a user-specified number of force
2362 evaluations has been performed (
{\eg} 100), or when the maximum of the absolute
2363 values of the force (gradient) components is smaller than a specified
2365 Since force truncation produces some noise in the
2366 energy evaluation, the stopping criterion should not be made too tight
2367 to avoid endless iterations. A reasonable value for $
\epsilon$ can be
2368 estimated from the root mean square force $f$ a harmonic oscillator would exhibit at a
2369 temperature $T$. This value is
2371 f =
2 \pi \nu \sqrt{ 2mkT
},
2373 where $
\nu$ is the oscillator frequency, $m$ the (reduced) mass, and
2374 $k$ Boltzmann's constant. For a weak oscillator with a wave number of
2375 100 cm$^
{-
1}$ and a mass of
10 atomic units, at a temperature of
1 K,
2376 $f=
7.7$ kJ~mol$^
{-
1}$~nm$^
{-
1}$. A value for $
\epsilon$ between
1 and
2379 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2380 \subsection{Conjugate Gradient
\index{conjugate gradient
}}
2381 Conjugate gradient is slower than steepest descent in the early stages
2382 of the minimization, but becomes more efficient closer to the energy
2383 minimum. The parameters and stop criterion are the same as for
2384 steepest descent. In
{\gromacs} conjugate gradient can not be used
2385 with constraints, including the SETTLE algorithm for
2386 water~
\cite{Miyamoto92
}, as this has not been implemented. If water is
2387 present it must be of a flexible model, which can be specified in the
2388 {\tt *.mdp
} file by
{\tt define = -DFLEXIBLE
}.
2390 This is not really a restriction, since the accuracy of conjugate
2391 gradient is only required for minimization prior to a normal-mode
2392 analysis, which cannot be performed with constraints. For most other
2393 purposes steepest descent is efficient enough.
2394 % } % Brace matches ifthenelse test for gmxlite
2396 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2397 \subsection{\normindex{L-BFGS
}}
2398 The original BFGS algorithm works by successively creating better
2399 approximations of the inverse Hessian matrix, and moving the system to
2400 the currently estimated minimum. The memory requirements for this are
2401 proportional to the square of the number of particles, so it is not
2402 practical for large systems like biomolecules. Instead, we use the
2403 L-BFGS algorithm of Nocedal~
\cite{Byrd95a,Zhu97a
}, which approximates
2404 the inverse Hessian by a fixed number of corrections from previous
2405 steps. This sliding-window technique is almost as efficient as the
2406 original method, but the memory requirements are much lower -
2407 proportional to the number of particles multiplied with the correction
2408 steps. In practice we have found it to converge faster than conjugate
2409 gradients, but due to the correction steps it is not yet parallelized.
2410 It is also noteworthy that switched or shifted interactions usually
2411 improve the convergence, since sharp cut-offs mean the potential
2412 function at the current coordinates is slightly different from the
2413 previous steps used to build the inverse Hessian approximation.
2414 % } % Brace matches ifthenelse test for gmxlite
2416 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2417 \section{Normal-Mode Analysis
\index{normal-mode analysis
}\index{NMA
}}
2418 Normal-mode analysis~
\cite{Levitt83,Go83,BBrooks83b
}
2419 can be performed using
{\gromacs}, by diagonalization of the mass-weighted
2420 \normindex{Hessian
} $H$:
2422 R^T M^
{-
1/
2} H M^
{-
1/
2} R &=&
\mbox{diag
}(
\lambda_1,
\ldots,
\lambda_{3N
})
2424 \lambda_i &=& (
2 \pi \omega_i)^
2
2426 where $M$ contains the atomic masses, $R$ is a matrix that contains
2427 the eigenvectors as columns, $
\lambda_i$ are the eigenvalues
2428 and $
\omega_i$ are the corresponding frequencies.
2430 First the Hessian matrix, which is a $
3N
\times 3N$ matrix where $N$
2431 is the number of atoms, needs to be calculated:
2433 H_
{ij
} &=&
\frac{\partial^
2 V
}{\partial x_i
\partial x_j
}
2435 where $x_i$ and $x_j$ denote the atomic x, y or z coordinates.
2436 In practice, this equation is not used, but the Hessian is
2437 calculated numerically from the force as:
2440 \frac{f_i(
{\bf x
}+h
{\bf e
}_j) - f_i(
{\bf x
}-h
{\bf e
}_j)
}{2h
}
2442 f_i &=& -
\frac{\partial V
}{\partial x_i
}
2444 where $
{\bf e
}_j$ is the unit vector in direction $j$.
2445 It should be noted that
2446 for a usual normal-mode calculation, it is necessary to completely minimize
2447 the energy prior to computation of the Hessian.
2448 The tolerance required depends on the type of system,
2449 but a rough indication is
0.001 kJ mol$^
{-
1}$.
2450 Minimization should be done with conjugate gradients or L-BFGS in double precision.
2452 A number of
{\gromacs} programs are involved in these
2453 calculations. First, the energy should be minimized using
{\tt mdrun
}.
2454 Then,
{\tt mdrun
} computes the Hessian.
{\bf Note
} that for generating
2455 the run input file, one should use the minimized conformation from
2456 the full precision trajectory file, as the structure file is not
2458 {\tt \normindex{g_nmeig
}} does the diagonalization and
2459 the sorting of the normal modes according to their frequencies.
2460 Both
{\tt mdrun
} and
{\tt g_nmeig
} should be run in double precision.
2461 The normal modes can be analyzed with the program
{\tt g_anaeig
}.
2462 Ensembles of structures at any temperature and for any subset of
2463 normal modes can be generated with
{\tt \normindex{g_nmens
}}.
2464 An overview of normal-mode analysis and the related principal component
2465 analysis (see
\secref{covanal
}) can be found in~
\cite{Hayward95b
}.
2466 % } % Brace matches ifthenelse test for gmxlite
2468 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2470 \section{Free energy calculations
\index{free energy calculations
}}
2472 \subsection{Slow-growth methods
\index{slow-growth methods
}}
2473 Free energy calculations can be performed
2474 in
{\gromacs} using a number of methods, including ``slow-growth.'' An example problem
2475 might be calculating the difference in free energy of binding of an inhibitor
{\bf I
}
2476 to an enzyme
{\bf E
} and to a mutated enzyme
{\bf E$^
{\prime}$
}. It
2477 is not feasible with computer simulations to perform a docking
2478 calculation for such a large complex, or even releasing the inhibitor from
2479 the enzyme in a reasonable amount of computer time with reasonable accuracy.
2480 However, if we consider the free energy cycle in~
\figref{free
}A
2483 \Delta G_1 -
\Delta G_2 =
\Delta G_3 -
\Delta G_4
2486 If we are interested in the left-hand term we can equally well compute
2487 the right-hand term.
2489 \centerline{\includegraphics[width=
6cm,angle=
270]{plots/free1
}\hspace{2cm
}\includegraphics[width=
6cm,angle=
270]{plots/free2
}}
2490 \caption[Free energy cycles.
]{Free energy cycles.
{\bf A:
} to
2491 calculate $
\Delta G_
{12}$, the free energy difference between the
2492 binding of inhibitor
{\bf I
} to enzymes
{\bf E
} respectively
{\bf
2493 E$^
{\prime}$
}.
{\bf B:
} to calculate $
\Delta G_
{12}$, the free energy
2494 difference for binding of inhibitors
{\bf I
} respectively
{\bf I$^
{\prime}$
} to
2499 If we want to compute the difference in free energy of binding of two
2500 inhibitors
{\bf I
} and
{\bf I$^
{\prime}$
} to an enzyme
{\bf E
} (
\figref{free
}B)
2501 we can again use
\eqnref{ddg
} to compute the desired property.
2503 \newcommand{\sA}{^
{\mathrm{A
}}}
2504 \newcommand{\sB}{^
{\mathrm{B
}}}
2505 Free energy differences between two molecular species can
2506 be calculated in
{\gromacs} using the ``slow-growth'' method.
2507 Such free energy differences between different molecular species are
2508 physically meaningless, but they can be used to obtain meaningful
2509 quantities employing a thermodynamic cycle.
2510 The method requires a simulation during which the Hamiltonian of the
2511 system changes slowly from that describing one system (A) to that
2512 describing the other system (B). The change must be so slow that the
2513 system remains in equilibrium during the process; if that requirement
2514 is fulfilled, the change is reversible and a slow-growth simulation from B to A
2515 will yield the same results (but with a different sign) as a slow-growth
2516 simulation from A to B. This is a useful check, but the user should be
2517 aware of the danger that equality of forward and backward growth results does
2518 not guarantee correctness of the results.
2520 The required modification of the Hamiltonian $H$ is realized by making
2521 $H$ a function of a
\textit{coupling parameter
} $
\lambda:
2522 H=H(p,q;
\lambda)$ in such a way that $
\lambda=
0$ describes system A
2523 and $
\lambda=
1$ describes system B:
2525 H(p,q;
0)=H
\sA (p,q);~~~~ H(p,q;
1)=H
\sB (p,q).
2527 In
{\gromacs}, the functional form of the $
\lambda$-dependence is
2528 different for the various force-field contributions and is described
2529 in section
\secref{feia
}.
2531 The Helmholtz free energy $A$ is related to the
2532 partition function $Q$ of an $N,V,T$ ensemble, which is assumed to be
2533 the equilibrium ensemble generated by a MD simulation at constant
2534 volume and temperature. The generally more useful Gibbs free energy
2535 $G$ is related to the partition function $
\Delta$ of an $N,p,T$
2536 ensemble, which is assumed to be the equilibrium ensemble generated by
2537 a MD simulation at constant pressure and temperature:
2539 A(
\lambda) &=& -k_BT
\ln Q \\
2540 Q &=& c
\int\!\!
\int \exp[-
\beta H(p,q;
\lambda)
]\,dp\,dq \\
2541 G(
\lambda) &=& -k_BT
\ln \Delta \\
2542 \Delta &=& c
\int\!\!
\int\!\!
\int \exp[-
\beta H(p,q;
\lambda) -
\beta
2546 where $
\beta =
1/(k_BT)$ and $c = (N! h^
{3N
})^
{-
1}$.
2547 These integrals over phase space cannot be evaluated from a
2548 simulation, but it is possible to evaluate the derivative with
2549 respect to $
\lambda$ as an ensemble average:
2551 \frac{dA
}{d
\lambda} =
\frac{\int\!\!
\int (
\partial H/
\partial
2552 \lambda)
\exp[-
\beta H(p,q;
\lambda)
]\,dp\,dq
}{\int\!\!
\int \exp[-
\beta
2553 H(p,q;
\lambda)
]\,dp\,dq
} =
2554 \left\langle \frac{\partial H
}{\partial \lambda} \right\rangle_{NVT;
\lambda},
2556 with a similar relation for $dG/d
\lambda$ in the $N,p,T$
2557 ensemble. The difference in free energy between A and B can be found
2558 by integrating the derivative over $
\lambda$:
2560 A
\sB(V,T)-A
\sA(V,T) &=&
\int_0^
1 \left\langle \frac{\partial
2561 H
}{\partial \lambda} \right\rangle_{NVT;
\lambda} \,d
\lambda
2563 G
\sB(p,T)-G
\sA(p,T) &=&
\int_0^
1 \left\langle \frac{\partial
2564 H
}{\partial \lambda} \right\rangle_{NpT;
\lambda} \,d
\lambda.
2567 If one wishes to evaluate $G
\sB(p,T)-G
\sA(p,T)$,
2568 the natural choice is a constant-pressure simulation. However, this
2569 quantity can also be obtained from a slow-growth simulation at
2570 constant volume, starting with system A at pressure $p$ and volume $V$
2571 and ending with system B at pressure $p_B$, by applying the following
2572 small (but, in principle, exact) correction:
2575 A
\sB(V)-A
\sA(V) -
\int_p^
{p
\sB}[V
\sB(p')-V
]\,dp'
2577 Here we omitted the constant $T$ from the notation. This correction is
2578 roughly equal to $-
\frac{1}{2} (p
\sB-p)
\Delta V=(
\Delta V)^
2/(
2
2579 \kappa V)$, where $
\Delta V$ is the volume change at $p$ and $
\kappa$
2580 is the isothermal compressibility. This is usually
2581 small; for example, the growth of a water molecule from nothing
2582 in a bath of
1000 water molecules at constant volume would produce an
2583 additional pressure of as much as
22 bar, but a correction to the
2584 Helmholtz free energy of just -
1 kJ mol$^
{-
1}$.
%-20 J/mol.
2586 In Cartesian coordinates, the kinetic energy term in the Hamiltonian
2587 depends only on the momenta, and can be separately integrated and, in
2588 fact, removed from the equations. When masses do not change, there is
2589 no contribution from the kinetic energy at all; otherwise the
2590 integrated contribution to the free energy is $-
\frac{3}{2} k_BT
\ln
2591 (m
\sB/m
\sA)$.
{\bf Note
} that this is only true in the absence of constraints.
2593 \subsection{Thermodynamic integration
\index{thermodynamic integration
}\index{BAR
}\index{Bennett's acceptance ratio
}}
2594 {\gromacs} offers the possibility to integrate eq.~
\ref{eq:delA
} or
2595 eq.
\ref{eq:delG
} in one simulation over the full range from A to
2596 B. However, if the change is large and insufficient sampling can be
2597 expected, the user may prefer to determine the value of $
\langle
2598 dG/d
\lambda \rangle$ accurately at a number of well-chosen
2599 intermediate values of $
\lambda$. This can easily be done by setting
2600 the stepsize
{\tt delta_lambda
} to zero. Each simulation can be
2601 equilibrated first, and a proper error estimate can be made for each
2602 value of $dG/d
\lambda$ from the fluctuation of $
\partial H/
\partial
2603 \lambda$. The total free energy change is then determined afterward
2604 by an appropriate numerical integration procedure.
2606 {\gromacs} now also supports the use of Bennett's Acceptance Ratio~
\cite{Bennett1976
}
2607 for calculating values of $
\Delta$G for transformations from state A to state B using
2608 the program
{\tt \normindex{g_bar
}}. The same data can also be used to calculate free
2609 energies using MBAR~
\cite{Shirts2008
}, though the analysis currently requires external tools from
2610 the external
{\tt pymbar
} package, at https://SimTK.org/home/pymbar.
2612 The $
\lambda$-dependence for the force-field contributions is
2613 described in detail in section
\secref{feia
}.
2614 % } % Brace matches ifthenelse test for gmxlite
2616 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2617 \section{Replica exchange
\index{replica exchange
}}
2618 Replica exchange molecular dynamics (
\normindex{REMD
})
2619 is a method that can be used to speed up
2620 the sampling of any type of simulation, especially if
2621 conformations are separated by relatively high energy barriers.
2622 It involves simulating multiple replicas of the same system
2623 at different temperatures and randomly exchanging the complete state
2624 of two replicas at regular intervals with the probability:
2626 P(
1 \leftrightarrow 2)=
\min\left(
1,
\exp\left[
2627 \left(
\frac{1}{k_B T_1
} -
\frac{1}{k_B T_2
}\right)(U_1 - U_2)
2630 where $T_1$ and $T_2$ are the reference temperatures and $U_1$ and $U_2$
2631 are the instantaneous potential energies of replicas
1 and
2 respectively.
2632 After exchange the velocities are scaled by $(T_1/T_2)^
{\pm0.5
}$
2633 and a neighbor search is performed the next step.
2634 This combines the fast sampling and frequent barrier-crossing
2635 of the highest temperature with correct Boltzmann sampling at
2636 all the different temperatures~
\cite{Hukushima96a,Sugita99
}.
2637 We only attempt exchanges for neighboring temperatures as the probability
2638 decreases very rapidly with the temperature difference.
2639 One should not attempt exchanges for all possible pairs in one step.
2640 If, for instance, replicas
1 and
2 would exchange, the chance of
2641 exchange for replicas
2 and
3 not only depends on the energies of
2642 replicas
2 and
3, but also on the energy of replica
1.
2643 In
{\gromacs} this is solved by attempting exchange for all ``odd''
2644 pairs on ``odd'' attempts and for all ``even'' pairs on ``even'' attempts.
2645 If we have four replicas:
0,
1,
2 and
3, ordered in temperature
2646 and we attempt exchange every
1000 steps, pairs
0-
1 and
2-
3
2647 will be tried at steps
1000,
3000 etc. and pair
1-
2 at steps
2000,
4000 etc.
2649 How should one choose the temperatures?
2650 The energy difference can be written as:
2652 U_1 - U_2 = N_
{df
} \frac{c
}{2} k_B (T_1 - T_2)
2654 where $N_
{df
}$ is the total number of degrees of freedom of one replica
2655 and $c$ is
1 for harmonic potentials and around
2 for protein/water systems.
2656 If $T_2 = (
1+
\epsilon) T_1$ the probability becomes:
2658 P(
1 \leftrightarrow 2)
2659 =
\exp\left( -
\frac{\epsilon^
2 c\,N_
{df
}}{2 (
1+
\epsilon)
} \right)
2660 \approx \exp\left(-
\epsilon^
2 \frac{c
}{2} N_
{df
} \right)
2662 Thus for a probability of $e^
{-
2}\approx 0.135$
2663 one obtains $
\epsilon \approx 2/
\sqrt{c\,N_
{df
}}$.
2664 With all bonds constrained one has $N_
{df
} \approx 2\, N_
{atoms
}$
2665 and thus for $c$ =
2 one should choose $
\epsilon$ as $
1/
\sqrt{N_
{atoms
}}$.
2666 However there is one problem when using pressure coupling. The density at
2667 higher temperatures will decrease, leading to higher energy~
\cite{Seibert2005a
},
2668 which should be taken into account. The
{\gromacs} website features a
2669 so-called ``REMD calculator,'' that lets you type in the temperature range and
2670 the number of atoms, and based on that proposes a set of temperatures.
2672 An extension to the REMD for the isobaric-isothermal ensemble was
2673 proposed by Okabe
{\em et al.
}~
\cite{Okabe2001a
}. In this work the
2674 exchange probability is modified to:
2676 P(
1 \leftrightarrow 2)=
\min\left(
1,
\exp\left[
2677 \left(
\frac{1}{k_B T_1
} -
\frac{1}{k_B T_2
}\right)(U_1 - U_2) +
2678 \left(
\frac{P_1
}{k_B T_1
} -
\frac{P_2
}{k_B T_2
}\right)
\left(V_1-V_2
\right)
2681 where $P_1$ and $P_2$ are the respective reference pressures and $V_1$ and
2682 $V_2$ are the respective instantaneous volumes in the simulations.
2683 In most cases the differences in volume are so small that the second
2684 term is negligible. It only plays a role when the difference between
2685 $P_1$ and $P_2$ is large or in phase transitions.
2687 Hamiltonian replica exchange is also supported in
{\gromacs}. In
2688 Hamiltonian replica exchange, each replica has a different
2689 Hamiltonian, defined by the free energy pathway specified for the simulation. The
2690 exchange probability to maintain the correct ensemble probabilities is:
2691 \beq P(
1 \leftrightarrow 2)=
\min\left(
1,
\exp\left[
2692 \left(
\frac{1}{k_B T
} -
\frac{1}{k_B T
}\right)((U_1(x_2) - U_1(x_1)) + (U_2(x_1) - U_2(x_2)))
2696 The separate Hamiltonians are defined by the free energy functionality
2697 of
{\gromacs}, with swaps made between the different values of
2698 $
\lambda$ defined in the mdp file.
2700 Hamiltonian and temperature replica exchange can also be performed
2701 simultaneously, using the acceptance criteria:
2703 P(
1 \leftrightarrow 2)=
\min\left(
1,
\exp\left[
2704 \left(
\frac{1}{k_B T
} -
\right)(
\frac{U_1(x_2) - U_1(x_1)
}{k_B T_1
} +
\frac{U_2(x_1) - U_2(x_2)
}{k_B T_2
})
2708 Gibbs sampling replica exchange has also been implemented in
2709 {\gromacs}~
\cite{Chodera2011
}. In Gibbs sampling replica exchange, all
2710 possible pairs are tested for exchange, allowing swaps between
2711 replicas that are not neighbors.
2713 Gibbs sampling replica exchange requires no additional potential
2714 energy calculations. However there is an additional communication
2715 cost in Gibbs sampling replica exchange, as for some permutations,
2716 more than one round of swaps must take place. In some cases, this
2717 extra communication cost might affect the efficiency.
2719 All replica exchange variants are options of the
{\tt mdrun
}
2720 program. It will only work when MPI is installed, due to the inherent
2721 parallelism in the algorithm. For efficiency each replica can run on a
2722 separate node. See the manual page of
{\tt mdrun
} on how to use these
2725 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2727 \section{Essential Dynamics sampling
\index{essential dynamics
}\index{principal component analysis
}\seeindexquiet{PCA
}{covariance analysis
}}
2728 The results from Essential Dynamics (see
\secref{covanal
})
2729 of a protein can be used to guide MD simulations. The idea is that
2730 from an initial MD simulation (or from other sources) a definition of
2731 the collective fluctuations with largest amplitude is obtained. The
2732 position along one or more of these collective modes can be
2733 constrained in a (second) MD simulation in a number of ways for
2734 several purposes. For example, the position along a certain mode may
2735 be kept fixed to monitor the average force (free-energy gradient) on
2736 that coordinate in that position. Another application is to enhance
2737 sampling efficiency with respect to usual MD
2738 \cite{Degroot96a,Degroot96b
}. In this case, the system is encouraged
2739 to sample its available configuration space more systematically than
2740 in a diffusion-like path that proteins usually take.
2742 Another possibility to enhance sampling is
\normindex{flooding
}.
2743 Here a flooding potential is added to certain
2744 (collective) degrees of freedom to expel the system out
2745 of a region of phase space
\cite{Lange2006a
}.
2747 The procedure for essential dynamics sampling or flooding is as follows.
2748 First, the eigenvectors and eigenvalues need to be determined
2749 using covariance analysis (
{\tt g_covar
})
2750 or normal-mode analysis (
{\tt g_nmeig
}).
2751 Then, this information is fed into
{\tt make_edi
},
2752 which has many options for selecting vectors and setting parameters,
2753 see
{\tt gmx make_edi -h
}.
2754 The generated
{\tt edi
} input file is then passed to
{\tt mdrun
}.
2756 % } % Brace matches ifthenelse test for gmxlite
2758 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2759 \section{\normindex{Expanded Ensemble
}}
2761 In an expanded ensemble simulation~
\cite{Lyubartsev1992
}, both the coordinates and the
2762 thermodynamic ensemble are treated as configuration variables that can
2763 be sampled over. The probability of any given state can be written as:
2765 P(
\vec{x
},k)
\propto \exp\left(-
\beta_k U_k + g_k
\right),
2767 where $
\beta_k =
\frac{1}{k_B T_k
}$ is the $
\beta$ corresponding to the $k$th
2768 thermodynamic state, and $g_k$ is a user-specified weight factor corresponding
2769 to the $k$th state. This space is therefore a
{\em mixed
},
{\em generalized
}, or
{\em
2770 expanded
} ensemble which samples from multiple thermodynamic
2771 ensembles simultaneously. $g_k$ is chosen to give a specific weighting
2772 of each subensemble in the expanded ensemble, and can either be fixed,
2773 or determined by an iterative procedure. The set of $g_k$ is
2774 frequently chosen to give each thermodynamic ensemble equal
2775 probability, in which case $g_k$ is equal to the free energy in
2776 non-dimensional units, but they can be set to arbitrary values as
2777 desired. Several different algorithms can be used to equilibrate
2778 these weights, described in the mdp option listings.
2779 % } % Brace matches ifthenelse test for gmxlite
2781 In
{\gromacs}, this space is sampled by alternating sampling in the $k$
2782 and $
\vec{x
}$ directions. Sampling in the $
\vec{x
}$ direction is done
2783 by standard molecular dynamics sampling; sampling between the
2784 different thermodynamics states is done by Monte Carlo, with several
2785 different Monte Carlo moves supported. The $k$ states can be defined
2786 by different temperatures, or choices of the free energy $
\lambda$
2787 variable, or both. Expanded ensemble simulations thus represent a
2788 serialization of the replica exchange formalism, allowing a single
2789 simulation to explore many thermodynamic states.
2793 % This stuff is quite outdated now.
2794 %\ifthenelse{\equal{\gmxlite}{1}}{}{
2796 %} % Brace matches ifthenelse test for gmxlite
2798 \section{Parallelization
\index{parallelization
}}
2799 The CPU time required for a simulation can be reduced by running the simulation
2800 in parallel over more than one processor or processor core.
2801 Ideally one would want to have linear scaling: running on $N$ processors/cores
2802 makes the simulation $N$ times faster. In practice this can only be
2803 achieved for a small number of processors. The scaling will depend
2804 a lot on the algorithms used. Also, different algorithms can have different
2805 restrictions on the interaction ranges between atoms.
2806 In
{\gromacs} we have two types of parallelization: particle decomposition
2807 and domain decomposition. Particle decomposition is only useful for
2808 a few special cases. Domain decomposition, which is the default algorithm,
2809 will always be faster and scale better.
2811 \section{Particle decomposition
\index{particle decomposition
}}
2812 Particle decomposition, also called
\index{force decomposition
},
2813 is the simplest type of decomposition. At the start of the simulation,
2814 particles are assigned to processors. Then forces between particles
2815 need to be assigned to processors such that the force load is evenly balanced.
2816 This decomposition requires that each processor know the coordinates
2817 of at least half of the particles in the system.
2818 Thus for a high number of processors $N$, about $N
\times N/
2$ coordinates
2819 need to be communicated. Because of this quadratic relation
2820 particle decomposition does not scale well.
2822 Particle decomposition was the only method available before version
4
2823 of
{\gromacs}. Now it is only useful in cases where domain decomposition
2824 does not work, such as systems with long-range bonded interactions,
2825 especially NMR distance or orientation restraints.
2826 With particle decomposition only whole molecules can be assigned to a processor.
2828 \section{Domain decomposition
\index{domain decomposition
}}
2829 Since most interactions in molecular simulations are local,
2830 domain decomposition is a natural way to decompose the system.
2831 In domain decomposition, a spatial domain is assigned to each processor,
2832 which will then integrate the equations of motion for the particles
2833 that currently reside in its local domain. With domain decomposition,
2834 there are two choices that have to be made: the division of the unit cell
2835 into domains and the assignment of the forces to processors.
2836 Most molecular simulation packages use the half-shell method for assigning
2837 the forces. But there are two methods that always require less communication:
2838 the eighth shell~
\cite{Liem1991
} and the midpoint~
\cite{Shaw2006
} method.
2839 {\gromacs} currently uses the eighth shell method, but for certain systems
2840 or hardware architectures it might be advantageous to use the midpoint
2841 method. Therefore, we might implement the midpoint method in the future.
2842 Most of the details of the domain decomposition can be found
2843 in the
{\gromacs} 4 paper~
\cite{Hess2008b
}.
2845 \subsection{Coordinate and force communication
}
2846 In the most general case of a triclinic unit cell,
2847 the space in divided with a
1-,
2-, or
3-D grid in parallelepipeds
2848 that we call domain decomposition cells.
2849 Each cell is assigned to a processor.
2850 The system is partitioned over the processors at the beginning
2851 of each MD step in which neighbor searching is performed.
2852 Since the neighbor searching is based on charge groups, charge groups
2853 are also the units for the domain decomposition.
2854 Charge groups are assigned to the cell where their center of geometry resides.
2855 Before the forces can be calculated, the coordinates from some
2856 neighboring cells need to be communicated,
2857 and after the forces are calculated, the forces need to be communicated
2858 in the other direction.
2859 The communication and force assignment is based on zones that
2860 can cover one or multiple cells.
2861 An example of a zone setup is shown in
\figref{ddcells
}.
2864 \centerline{\includegraphics[width=
6cm
]{plots/dd-cells
}}
2866 A non-staggered domain decomposition grid of
3$
\times$
2$
\times$
2 cells.
2867 Coordinates in zones
1 to
7 are communicated to the corner cell
2868 that has its home particles in zone
0.
2869 $r_c$ is the cut-off radius.
2874 The coordinates are communicated by moving data along the ``negative''
2875 direction in $x$, $y$ or $z$ to the next neighbor. This can be done in one
2876 or multiple pulses. In
\figref{ddcells
} two pulses in $x$ are required,
2877 then one in $y$ and then one in $z$. The forces are communicated by
2878 reversing this procedure. See the
{\gromacs} 4 paper~
\cite{Hess2008b
}
2879 for details on determining which non-bonded and bonded forces
2880 should be calculated on which node.
2882 \subsection{Dynamic load balancing
\swapindexquiet{dynamic
}{load balancing
}}
2883 When different processors have a different computational load
2884 (load imbalance), all processors will have to wait for the one
2885 that takes the most time. One would like to avoid such a situation.
2886 Load imbalance can occur due to three reasons:
2888 \item inhomogeneous particle distribution
2889 \item inhomogeneous interaction cost distribution (charged/uncharged,
2890 water/non-water due to
{\gromacs} water innerloops)
2891 \item statistical fluctuation (only with small particle numbers)
2893 So we need a dynamic load balancing algorithm
2894 where the volume of each domain decomposition cell
2895 can be adjusted
{\em independently
}.
2896 To achieve this, the
2- or
3-D domain decomposition grids need to be
2897 staggered.
\figref{ddtric
} shows the most general case in
2-D.
2898 Due to the staggering, one might require two distance checks
2899 for deciding if a charge group needs to be communicated:
2900 a non-bonded distance and a bonded distance check.
2903 \centerline{\includegraphics[width=
7cm
]{plots/dd-tric
}}
2905 The zones to communicate to the processor of zone
0,
2906 see the text for details. $r_c$ and $r_b$ are the non-bonded
2907 and bonded cut-off radii respectively, $d$ is an example
2908 of a distance between following, staggered boundaries of cells.
2913 By default,
{\tt mdrun
} automatically turns on the dynamic load
2914 balancing during a simulation when the total performance loss
2915 due to the force calculation imbalance is
5\% or more.
2916 {\bf Note
} that the reported force load imbalance numbers might be higher,
2917 since the force calculation is only part of work that needs to be done
2918 during an integration step.
2919 The load imbalance is reported in the log file at log output steps
2920 and when the
{\tt -v
} option is used also on screen.
2921 The average load imbalance and the total performance loss
2922 due to load imbalance are reported at the end of the log file.
2924 There is one important parameter for the dynamic load balancing,
2925 which is the minimum allowed scaling. By default, each dimension
2926 of the domain decomposition cell can scale down by at least
2927 a factor of
0.8. For
3-D domain decomposition this allows cells
2928 to change their volume by about a factor of
0.5, which should allow
2929 for compensation of a load imbalance of
100\%.
2930 The required scaling can be changed with the
{\tt -dds
} option of
{\tt mdrun
}.
2932 \subsection{Constraints in parallel
\index{constraints
}}
2933 \label{subsec:plincs
}
2934 Since with domain decomposition parts of molecules can reside
2935 on different processors, bond constraints can cross cell boundaries.
2936 Therefore a parallel constraint algorithm is required.
2937 {\gromacs} uses the
\normindex{P-LINCS
} algorithm~
\cite{Hess2008a
},
2938 which is the parallel version of the
\normindex{LINCS
} algorithm~
\cite{Hess97
}
2939 % \ifthenelse{\equal{\gmxlite}{1}}
2941 {(see
\ssecref{lincs
}).
}
2942 The P-LINCS procedure is illustrated in
\figref{plincs
}.
2943 When molecules cross the cell boundaries, atoms in such molecules
2944 up to (
{\tt lincs_order +
1}) bonds away are communicated over the cell boundaries.
2945 Then, the normal LINCS algorithm can be applied to the local bonds
2946 plus the communicated ones. After this procedure, the local bonds
2947 are correctly constrained, even though the extra communicated ones are not.
2948 One coordinate communication step is required for the initial LINCS step
2949 and one for each iteration. Forces do not need to be communicated.
2952 \centerline{\includegraphics[width=
6cm
]{plots/par-lincs2
}}
2954 Example of the parallel setup of P-LINCS with one molecule
2955 split over three domain decomposition cells, using a matrix
2956 expansion order of
3.
2957 The top part shows which atom coordinates need to be communicated
2958 to which cells. The bottom parts show the local constraints (solid)
2959 and the non-local constraints (dashed) for each of the three cells.
2964 \subsection{Interaction ranges
}
2965 Domain decomposition takes advantage of the locality of interactions.
2966 This means that there will be limitations on the range of interactions.
2967 By default,
{\tt mdrun
} tries to find the optimal balance between
2968 interaction range and efficiency. But it can happen that a simulation
2969 stops with an error message about missing interactions,
2970 or that a simulation might run slightly faster with shorter
2971 interaction ranges. A list of interaction ranges
2972 and their default values is given in
\tabref{dd_ranges
}.
2976 \begin{tabular
}{|c|c|ll|
}
2978 interaction & range & option & default \\
2980 non-bonded & $r_c$ = max($r_
{list
}$,$r_
{VdW
}$,$r_
{Coul
}$) &
{\tt mdp
} file & \\
2981 two-body bonded & max($r_
{mb
}$,$r_c$) &
{\tt mdrun -rdd
} & starting conf. +
10\% \\
2982 multi-body bonded & $r_
{mb
}$ &
{\tt mdrun -rdd
} & starting conf. +
10\% \\
2983 constraints & $r_
{con
}$ &
{\tt mdrun -rcon
} & est. from bond lengths \\
2984 virtual sites & $r_
{con
}$ &
{\tt mdrun -rcon
} &
0 \\
2988 \caption{The interaction ranges with domain decomposition.
}
2989 \label{tab:dd_ranges
}
2992 In most cases the defaults of
{\tt mdrun
} should not cause the simulation
2993 to stop with an error message of missing interactions.
2994 The range for the bonded interactions is determined from the distance
2995 between bonded charge-groups in the starting configuration, with
10\% added
2996 for headroom. For the constraints, the value of $r_
{con
}$ is determined by
2997 taking the maximum distance that (
{\tt lincs_order +
1}) bonds can cover
2998 when they all connect at angles of
120 degrees.
2999 The actual constraint communication is not limited by $r_
{con
}$,
3000 but by the minimum cell size $L_C$, which has the following lower limit:
3002 L_C
\geq \max(r_
{mb
},r_
{con
})
3004 Without dynamic load balancing the system is actually allowed to scale
3005 beyond this limit when pressure scaling is used.
3006 {\bf Note
} that for triclinic boxes, $L_C$ is not simply the box diagonal
3007 component divided by the number of cells in that direction,
3008 rather it is the shortest distance between the triclinic cells borders.
3009 For rhombic dodecahedra this is a factor of $
\sqrt{3/
2}$ shorter
3012 When $r_
{mb
} > r_c$,
{\tt mdrun
} employs a smart algorithm to reduce
3013 the communication. Simply communicating all charge groups within
3014 $r_
{mb
}$ would increase the amount of communication enormously.
3015 Therefore only charge-groups that are connected by bonded interactions
3016 to charge groups which are not locally present are communicated.
3017 This leads to little extra communication, but also to a slightly
3018 increased cost for the domain decomposition setup.
3019 In some cases,
{\eg} coarse-grained simulations with a very short cut-off,
3020 one might want to set $r_
{mb
}$ by hand to reduce this cost.
3022 \subsection{Multiple-Program, Multiple-Data PME parallelization
\index{PME
}}
3023 \label{subsec:mpmd_pme
}
3024 Electrostatics interactions are long-range, therefore special
3025 algorithms are used to avoid summation over many atom pairs.
3026 In
{\gromacs} this is usually
3027 % \ifthenelse{\equal{\gmxlite}{1}}
3029 {PME (
\secref{pme
}).
}
3030 Since with PME all particles interact with each other, global communication
3031 is required. This will usually be the limiting factor for
3032 scaling with domain decomposition.
3033 To reduce the effect of this problem, we have come up with
3034 a Multiple-Program, Multiple-Data approach~
\cite{Hess2008b
}.
3035 Here, some processors are selected to do only the PME mesh calculation,
3036 while the other processors, called particle-particle (PP) nodes,
3037 do all the rest of the work.
3038 For rectangular boxes the optimal PP to PME node ratio is usually
3:
1,
3039 for rhombic dodecahedra usually
2:
1.
3040 When the number of PME nodes is reduced by a factor of
4, the number
3041 of communication calls is reduced by about a factor of
16.
3042 Or put differently, we can now scale to
4 times more nodes.
3043 In addition, for modern
4 or
8 core machines in a network,
3044 the effective network bandwidth for PME is quadrupled,
3045 since only a quarter of the cores will be using the network connection
3046 on each machine during the PME calculations.
3049 \centerline{\includegraphics[width=
12cm
]{plots/mpmd-pme
}}
3051 Example of
8 nodes without (left) and with (right) MPMD.
3052 The PME communication (red arrows) is much higher on the left
3053 than on the right. For MPMD additional PP - PME coordinate
3054 and force communication (blue arrows) is required,
3055 but the total communication complexity is lower.
3056 \label{fig:mpmd_pme
}
3060 {\tt mdrun
} will by default interleave the PP and PME nodes.
3061 If the processors are not number consecutively inside the machines,
3062 one might want to use
{\tt mdrun -ddorder pp_pme
}.
3063 For machines with a real
3-D torus and proper communication software
3064 that assigns the processors accordingly one should use
3065 {\tt mdrun -ddorder cartesian
}.
3067 To optimize the performance one should usually set up the cut-offs
3068 and the PME grid such that the PME load is
25 to
33\% of the total
3069 calculation load.
{\tt grompp
} will print an estimate for this load
3070 at the end and also
{\tt mdrun
} calculates the same estimate
3071 to determine the optimal number of PME nodes to use.
3072 For high parallelization it might be worthwhile to optimize
3073 the PME load with the
{\tt mdp
} settings and/or the number
3074 of PME nodes with the
{\tt -npme
} option of
{\tt mdrun
}.
3075 For changing the electrostatics settings it is useful to know
3076 the accuracy of the electrostatics remains nearly constant
3077 when the Coulomb cut-off and the PME grid spacing are scaled
3079 {\bf Note
} that it is usually better to overestimate than to underestimate
3080 the number of PME nodes, since the number of PME nodes is smaller
3081 than the number of PP nodes, which leads to less total waiting time.
3083 The PME domain decomposition can be
1-D or
2-D along the $x$ and/or
3084 $y$ axis.
2-D decomposition is also known as
\normindex{pencil decomposition
} because of
3085 the shape of the domains at high parallelization.
3086 1-D decomposition along the $y$ axis can only be used when
3087 the PP decomposition has only
1 domain along $x$.
2-D PME decomposition
3088 has to have the number of domains along $x$ equal to the number of
3089 the PP decomposition.
{\tt mdrun
} automatically chooses
1-D or
2-D
3090 PME decomposition (when possible with the total given number of nodes),
3091 based on the minimum amount of communication for the coordinate redistribution
3092 in PME plus the communication for the grid overlap and transposes.
3093 To avoid superfluous communication of coordinates and forces
3094 between the PP and PME nodes, the number of DD cells in the $x$
3095 direction should ideally be the same or a multiple of the number
3096 of PME nodes. By default,
{\tt mdrun
} takes care of this issue.
3098 \subsection{Domain decomposition flow chart
}
3099 In
\figref{dd_flow
} a flow chart is shown for domain decomposition
3100 with all possible communication for different algorithms.
3101 For simpler simulations, the same flow chart applies,
3102 without the algorithms and communication for
3103 the algorithms that are not used.
3106 \centerline{\includegraphics[width=
12cm
]{plots/flowchart
}}
3108 Flow chart showing the algorithms and communication (arrows)
3109 for a standard MD simulation with virtual sites, constraints
3110 and separate PME-mesh nodes.
3116 \section{Implicit solvation
\index{implicit solvation
}\index{Generalized Born methods
}}
3118 Implicit solvent models provide an efficient way of representing
3119 the electrostatic effects of solvent molecules, while saving a
3120 large piece of the computations involved in an accurate, aqueous
3121 description of the surrounding water in molecular dynamics simulations.
3122 Implicit solvation models offer several advantages compared with
3123 explicit solvation, including eliminating the need for the equilibration of water
3124 around the solute, and the absence of viscosity, which allows the protein
3125 to more quickly explore conformational space.
3127 Implicit solvent calculations in
{\gromacs} can be done using the
3128 generalized Born-formalism, and the Still~
\cite{Still97
}, HCT~
\cite{Truhlar96
},
3129 and OBC~
\cite{Case04
} models are available for calculating the Born radii.
3131 Here, the free energy $G_
{solv
}$ of solvation is the sum of three terms,
3132 a solvent-solvent cavity term ($G_
{cav
}$), a solute-solvent van der
3133 Waals term ($G_
{vdw
}$), and finally a solvent-solute electrostatics
3134 polarization term ($G_
{pol
}$).
3136 The sum of $G_
{cav
}$ and $G_
{vdw
}$ corresponds to the (non-polar)
3137 free energy of solvation for a molecule from which all charges
3138 have been removed, and is commonly called $G_
{np
}$,
3139 calculated from the total solvent accessible surface area
3140 multiplied with a surface tension.
3141 The total expression for the solvation free energy then becomes:
3144 G_
{solv
} = G_
{np
} + G_
{pol
}
3148 Under the generalized Born model, $G_
{pol
}$ is calculated from the generalized Born equation~
\cite{Still97
}:
3151 G_
{pol
} =
\left(
1-
\frac{1}{\epsilon}\right)
\sum_{i=
1}^n
\sum_{j>i
}^n
\frac {q_i q_j
}{\sqrt{r^
2_
{ij
} + b_i b_j
\exp\left(
\frac{-r^
2_
{ij
}}{4 b_i b_j
}\right)
}}
3152 \label{eqn:gb_still
}
3155 In
{\gromacs}, we have introduced the substitution~
\cite{Larsson10
}:
3158 c_i=
\frac{1}{\sqrt{b_i
}}
3159 \label{eqn:gb_subst
}
3162 which makes it possible to introduce a cheap transformation to a new
3163 variable $x$ when evaluating each interaction, such that:
3166 x=
\frac{r_
{ij
}}{\sqrt{b_i b_j
}} = r_
{ij
} c_i c_j
3167 \label{eqn:gb_subst2
}
3170 In the end, the full re-formulation of~
\ref{eqn:gb_still
} becomes:
3173 G_
{pol
} =
\left(
1-
\frac{1}{\epsilon}\right)
\sum_{i=
1}^n
\sum_{j>i
}^n
\frac{q_i q_j
}{\sqrt{b_i b_j
}} ~
\xi (x) =
\left(
1-
\frac{1}{\epsilon}\right)
\sum_{i=
1}^n q_i c_i
\sum_{j>i
}^n q_j c_j~
\xi (x)
3174 \label{eqn:gb_final
}
3177 The non-polar part ($G_
{np
}$) of Equation~
\ref{eqn:gb_solv
} is calculated
3178 directly from the Born radius of each atom using a simple ACE type
3179 approximation by Schaefer
{\em et al.
}~
\cite{Karplus98
}, including a
3180 simple loop over all atoms.
3181 This requires only one extra solvation parameter, independent of atom type,
3182 but differing slightly between the three Born radii models.
3184 % LocalWords: GROningen MAchine BIOSON Groningen GROMACS Berendsen der Spoel
3185 % LocalWords: Drunen Comp Phys Comm ROck NS FFT pbc EM ifthenelse gmxlite ff
3186 % LocalWords: octahedra triclinic Ewald PME PPPM trjconv xy solvated
3187 % LocalWords: boxtypes boxshapes editconf Lennard mdpopt COM XTC TNG kT defunits
3188 % LocalWords: Boltzmann's Mueller nb int mdrun chargegroup simplerc prefactor
3189 % LocalWords: pme waterloops CH NH CO df com virial integrator Verlet vverlet
3190 % LocalWords: integrators ref timepoint timestep timesteps mdp md vv avek NVE
3191 % LocalWords: NVT off's leapfrogv lll LR rmfast SPC fs Nos physicality ps GMX
3192 % LocalWords: Tcoupling nonergodic thermostatting NOSEHOOVER algorithmes ij yx
3193 % LocalWords: Parrinello Rahman rescales atm anisotropically ccc xz zx yy yz
3194 % LocalWords: zy zz se barostat compressibilities MTTK NPT Martyna al isobaric
3195 % LocalWords: Tuckerman vir PV fkT iLt iL Liouville NHC Eq baro mu trj mol bc
3196 % LocalWords: freezegroup Shannon's polarizability Overhauser barostats iLn KE
3197 % LocalWords: negligibly thermostatted Tobias rhombic maxwell et xtc tng TC rlist
3198 % LocalWords: waals LINCS holonomic plincs lincs unc ang SA Langevin SD amu BD
3199 % LocalWords: bfgs Broyden Goldfarb Shanno mkT kJ DFLEXIBLE Nocedal diag nmeig
3200 % LocalWords: diagonalization anaeig nmens covanal ddg feia BT dp dq pV dV dA
3201 % LocalWords: NpT eq stepsize REMD constrainted website Okabe MPI covar edi dd
3202 % LocalWords: progman NMR ddcells innerloops ddtric tric dds rdd conf rcon est
3203 % LocalWords: mb PP MPMD ddorder pp cartesian grompp npme parallelizable edr
3204 % LocalWords: macromolecule nstlist vacuo parallelization dof indices MBAR AVX
3205 % LocalWords: TOL numerics parallelized eigenvectors dG parallelepipeds VdW np
3206 % LocalWords: Coul multi solvation HCT OBC solv cav vdw Schaefer symplectic dt
3207 % LocalWords: pymbar multinode subensemble Monte solute subst groupconcept GPU
3208 % LocalWords: dodecahedron octahedron dodecahedra equilibration usinggroups nm
3209 % LocalWords: topologies rlistlong CUDA GPUs rcoulomb SIMD BlueGene FPUs erfc
3210 % LocalWords: cutoffschemesupport unbuffered bondeds AdResS OpenMP ewald rtol
3211 % LocalWords: verletdrift peptide RMS rescaling ergodicity ergodic discretized
3212 % LocalWords: isothermal compressibility isotropically anisotropic iteratively
3213 % LocalWords: incompressible integrations translational biomolecules NMA PCA
3214 % LocalWords: Bennett's equilibrated Hamiltonians covariance equilibrate
3215 % LocalWords: inhomogeneous conformational online other's th