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 +
652 \phantom{4 \pi (r_
\ell+
\sigma)^
2 \rho_2 \int_{-
\infty}^
{r_c
} \int_{r_
\ell}^
\infty \Big[}
653 V'''(r_c)
\frac{1}{6}(r_t - r_c)^
3 \Big] G\!
\left(
\frac{r_t-r_0
}{\sigma}\right)
656 4 \pi (r_
\ell+
\sigma)^
2 \rho_2 \bigg\
{
657 \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] +
660 \phantom{4 \pi (r_
\ell+
\sigma)^
2 \rho_2 \bigg\
{ }
661 \frac{1}{6}V''(r_c)
\left[ \sigma(r_b^
2+
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] +
664 \phantom{4 \pi (r_
\ell+
\sigma)^
2 \rho_2 \bigg\
{ }
665 \frac{1}{24}V'''(r_c)
\left[ r_b
\sigma(r_b^
2+
5\sigma^
2)G\!
\left(
\frac{r_b
}{\sigma}\right) - (r_b^
4+
6r_b^
2\sigma^
2+
3\sigma^
4 ) E\!
\left(
\frac{r_b
}{\sigma}\right)
\right]
669 where $G$ is a Gaussian distribution with
0 mean and unit variance and
670 $E(x)=
\frac{1}{2}\mathrm{erfc
}(x/
\sqrt{2})$. We always want to achieve
671 small energy error, so $
\sigma$ will be small compared to both $r_c$
672 and $r_
\ell$, thus the approximations in the equations above are good,
673 since the Gaussian distribution decays rapidly. The energy error needs
674 to be averaged over all particle pair types and weighted with the
675 particle counts. In
{\gromacs} we don't allow cancellation of error
676 between pair types, so we average the absolute values. To obtain the
677 average energy error per unit time, it needs to be divided by the
678 neighbor-list life time $t = (
{\tt nstlist
} -
1)
\times{\tt dt
}$. This
679 function can not be inverted analytically, so we use bisection to
680 obtain the buffer size $r_b$ for a target drift. Again we note that
681 in practice the error we usually be much smaller than this estimate,
682 as in the condensed phase particle displacements will be much smaller
683 than for freely moving particles, which is the assumption used here.
685 When (bond) constraints are present, some particles will have fewer
686 degrees of freedom. This will reduce the energy errors. The
687 displacement in an arbitrary direction of a particle with
2 degrees of
688 freedom is not Gaussian, but rather follows the complementary error
690 \frac{\sqrt{\pi}}{2\sqrt{2}\sigma}\,
\mathrm{erfc
}\left(
\frac{|r|
}{\sqrt{2}\,
\sigma}\right)
691 \eeq where $
\sigma^
2$ is again $k_B T/m$. This distribution can no
692 longer be integrated analytically to obtain the energy error. But we
693 can generate a tight upper bound using a scaled and shifted Gaussian
694 distribution (not shown). This Gaussian distribution can then be used
695 to calculate the energy error as described above. We consider
696 particles constrained, i.e. having
2 degrees of freedom or fewer, when
697 they are connected by constraints to particles with a total mass of at
698 least
1.5 times the mass of the particles itself. For a particle with
699 a single constraint this would give a total mass along the constraint
700 direction of at least
2.5, which leads to a reduction in the variance
701 of the displacement along that direction by at least a factor of
6.25.
702 As the Gaussian distribution decays very rapidly, this effectively
703 removes one degree of freedom from the displacement. Multiple
704 constraints would reduce the displacement even more, but as this gets
705 very complex, we consider those as particles with
2 degrees of
708 There is one important implementation detail that reduces the energy
709 errors caused by the finite Verlet buffer list size. The derivation
710 above assumes a particle pair-list. However, the
{\gromacs}
711 implementation uses a cluster pair-list for efficiency. The pair list
712 consists of pairs of clusters of
4 particles in most cases, also
713 called a $
4 \times 4$ list, but the list can also be $
4 \times 8$ (GPU
714 CUDA kernels and AVX
256-bit single precision kernels) or $
4 \times 2$
715 (SSE double-precision kernels). This means that the pair-list is
716 effectively much larger than the corresponding $
1 \times 1$ list. Thus
717 slightly beyond the pair-list cut-off there will still be a large
718 fraction of particle pairs present in the list. This fraction can be
719 determined in a simulation and accurately estimated under some
720 reasonable assumptions. The fraction decreases with increasing
721 pair-list range, meaning that a smaller buffer can be used. For
722 typical all-atom simulations with a cut-off of
0.9 nm this fraction is
723 around
0.9, which gives a reduction in the energy errors of a factor of
724 10. This reduction is taken into account during the automatic Verlet
725 buffer calculation and results in a smaller buffer size.
728 \centerline{\includegraphics[width=
9cm
]{plots/verlet-drift
}}
729 \caption {Energy drift per atom for an SPC/E water system at
300K with
730 a time step of
2 fs and a pair-list update period of
10 steps
731 (pair-list life time:
18 fs). PME was used with
{\tt ewald-rtol
} set
732 to
10$^
{-
5}$; this parameter affects the shape of the potential at
733 the cut-off. Error estimates due to finite Verlet buffer size are
734 shown for a $
1 \times 1$ atom pair list and $
4 \times 4$ atom pair
735 list without and with (dashed line) cancellation of positive and
736 negative errors. Real energy drift is shown for double- and
737 single-precision simulations. Single-precision rounding errors in
738 the SETTLE constraint algorithm cause the drift to become negative
739 at large buffer size. Note that at zero buffer size, the real drift
740 is small because positive (H-H) and negative (O-H) energy errors
742 \label{fig:verletdrift
}
745 In
\figref{verletdrift
} one can see that for small buffer sizes the drift
746 of the total energy is much smaller than the pair energy error tolerance,
747 due to cancellation of errors. For larger buffer size, the error estimate
748 is a factor of
6 higher than drift of the total energy, or alternatively
749 the buffer estimate is
0.024 nm too large. This is because the protons
750 don't move freely over
18 fs, but rather vibrate.
751 %At a buffer size of zero there is cancellation of
752 %drift due to repulsive (H-H) and attractive (O-H) interactions.
754 \subsubsection{Cut-off artifacts and switched interactions
}
755 With the Verlet scheme, the pair potentials are shifted to be zero at
756 the cut-off, such that the potential is the integral of the force.
757 Note that in the group scheme this is not possible, because no exact
758 cut-off distance is used. There can still be energy drift from
759 non-zero forces at the cut-off. This effect is extremely small and
760 often not noticeable, as other integration errors may dominate. To
761 completely avoid cut-off artifacts, the non-bonded forces can be
762 switched exactly to zero at some distance smaller than the neighbor
763 list cut-off (there are several ways to do this in
{\gromacs}, see
764 \secref{mod_nb_int
}). One then has a buffer with the size equal to the
765 neighbor list cut-off less the longest interaction cut-off. With the
766 group cut-off scheme, one can then also choose to let
{\tt mdrun
} only
767 update the neighbor list when required. That is when one or more
768 particles have moved more than half the buffer size from the center of
769 geometry of the
\swapindex{charge
}{group
} to which they belong (see
770 \secref{chargegroup
}), as determined at the previous neighbor search.
771 This option guarantees that there are no cut-off artifacts.
{\bf
772 Note
} that for larger systems this comes at a high computational
773 cost, since the neighbor list update frequency will be determined by
774 just one or two particles moving slightly beyond the half buffer
775 length (which not even necessarily implies that the neighbor list is
776 invalid), while
99.99\% of the particles are fine.
} % Brace matches
777 ifthenelse test for gmxlite
779 \subsubsection{Simple search
\swapindexquiet{simple
}{search
}}
780 Due to
\eqnsref{box_rot
}{simplerc
}, the vector $
\rvij$
781 connecting images within the cut-off $R_c$ can be found by constructing:
783 \ve{r
}''' & = &
\ve{r
}_j-
\ve{r
}_i \\
784 \ve{r
}'' & = &
\ve{r
}''' -
{\bf c
}*
\verb'round'(r'''_z/c_z) \\
785 \ve{r
}' & = &
\ve{r
}'' -
{\bf b
}*
\verb'round'(r''_y/b_y) \\
786 \ve{r
}_
{ij
} & = &
\ve{r
}' -
{\bf a
}*
\verb'round'(r'_x/a_x)
788 When distances between two particles in a triclinic box are needed
789 that do not obey
\eqnref{box_rot
},
790 many shifts of combinations of box vectors need to be considered to find
793 \ifthenelse{\equal{\gmxlite}{1}}{}{
796 \centerline{\includegraphics[width=
8cm
]{plots/nstric
}}
797 \caption {Grid search in two dimensions. The arrows are the box vectors.
}
801 \subsubsection{Grid search
\swapindexquiet{grid
}{search
}}
803 The grid search is schematically depicted in
\figref{grid
}. All
804 particles are put on the
{\nsgrid}, with the smallest spacing $
\ge$
805 $R_c/
2$ in each of the directions. In the direction of each box
806 vector, a particle $i$ has three images. For each direction the image
807 may be -
1,
0 or
1, corresponding to a translation over -
1,
0 or +
1 box
808 vector. We do not search the surrounding
{\nsgrid} cells for neighbors
809 of $i$ and then calculate the image, but rather construct the images
810 first and then search neighbors corresponding to that image of $i$.
811 As
\figref{grid
} shows, some grid cells may be searched more than once
812 for different images of $i$. This is not a problem, since, due to the
813 minimum image convention, at most one image will ``see'' the
814 $j$-particle. For every particle, fewer than
125 (
5$^
3$) neighboring
815 cells are searched. Therefore, the algorithm scales linearly with the
816 number of particles. Although the prefactor is large, the scaling
817 behavior makes the algorithm far superior over the standard $O(N^
2)$
818 algorithm when there are more than a few hundred particles. The
819 grid search is equally fast for rectangular and triclinic boxes. Thus
820 for most protein and peptide simulations the rhombic dodecahedron will
821 be the preferred box shape.
822 } % Brace matches ifthenelse test for gmxlite
824 \ifthenelse{\equal{\gmxlite}{1}}{}{
825 \subsubsection{Charge groups
}
826 \label{sec:chargegroup
}\swapindexquiet{charge
}{group
}%
827 Charge groups were originally introduced to reduce cut-off artifacts
828 of Coulomb interactions. When a plain cut-off is used, significant
829 jumps in the potential and forces arise when atoms with (partial) charges
830 move in and out of the cut-off radius. When all chemical moieties have
831 a net charge of zero, these jumps can be reduced by moving groups
832 of atoms with net charge zero, called charge groups, in and
833 out of the neighbor list. This reduces the cut-off effects from
834 the charge-charge level to the dipole-dipole level, which decay
835 much faster. With the advent of full range electrostatics methods,
836 such as particle mesh Ewald (
\secref{pme
}), the use of charge groups is
837 no longer required for accuracy. It might even have a slight negative effect
838 on the accuracy or efficiency, depending on how the neighbor list is made
839 and the interactions are calculated.
841 But there is still an important reason for using ``charge groups'': efficiency.
842 Where applicable, neighbor searching is carried out on the basis of
843 charge groups which are defined in the molecular topology.
844 If the nearest image distance between the
{\em
845 geometrical centers
} of the atoms of two charge groups is less than
846 the cut-off radius, all atom pairs between the charge groups are
847 included in the pair list.
848 The neighbor searching for a water system, for instance,
849 is $
3^
2=
9$ times faster when each molecule is treated as a charge group.
850 Also the highly optimized water force loops (see
\secref{waterloops
})
851 only work when all atoms in a water molecule form a single charge group.
852 Currently the name
{\em neighbor-search group
} would be more appropriate,
853 but the name charge group is retained for historical reasons.
854 When developing a new force field, the advice is to use charge groups
855 of
3 to
4 atoms for optimal performance. For all-atom force fields
856 this is relatively easy, as one can simply put hydrogen atoms, and in some
857 case oxygen atoms, in the same charge group as the heavy atom they
858 are connected to; for example: CH$_3$, CH$_2$, CH, NH$_2$, NH, OH, CO$_2$, CO.
859 } % Brace matches ifthenelse test for gmxlite
861 \subsection{Compute forces
}
862 \label{subsec:forces
}
864 \subsubsection{Potential energy
}
865 When forces are computed, the
\swapindex{potential
}{energy
} of each
866 interaction term is computed as well. The total potential energy is
867 summed for various contributions, such as Lennard-Jones, Coulomb, and
868 bonded terms. It is also possible to compute these contributions for
869 {\em energy-monitor groups
} of atoms that are separately defined (see
870 \secref{groupconcept
}).
872 \subsubsection{Kinetic energy and temperature
}
873 The
\normindex{temperature
} is given by the total
874 \swapindex{kinetic
}{energy
} of the $N$-particle system:
876 E_
{kin
} =
\half \sum_{i=
1}^N m_i v_i^
2
878 From this the absolute temperature $T$ can be computed using:
880 \half N_
{df
} kT = E_
{kin
}
883 where $k$ is Boltzmann's constant and $N_
{df
}$ is the number of
884 degrees of freedom which can be computed from:
886 N_
{df
} ~=~
3 N - N_c - N_
{com
}
888 Here $N_c$ is the number of
{\em \normindex{constraints
}} imposed on the system.
889 When performing molecular dynamics $N_
{com
}=
3$ additional degrees of
890 freedom must be removed, because the three
891 center-of-mass velocities are constants of the motion, which are usually
892 set to zero. When simulating in vacuo, the rotation around the center of mass
893 can also be removed, in this case $N_
{com
}=
6$.
894 When more than one temperature-coupling group
\index{temperature-coupling group
} is used, the number of degrees
895 of freedom for group $i$ is:
897 N^i_
{df
} ~=~ (
3 N^i - N^i_c)
\frac{3 N - N_c - N_
{com
}}{3 N - N_c
}
900 The kinetic energy can also be written as a tensor, which is necessary
901 for pressure calculation in a triclinic system, or systems where shear
904 {\bf E
}_
{kin
} =
\half \sum_i^N m_i
\vvi \otimes \vvi
907 \subsubsection{Pressure and virial
}
908 The
\normindex{pressure
}
909 tensor
{\bf P
} is calculated from the difference between
910 kinetic energy $E_
{kin
}$ and the
\normindex{virial
} $
{\bf \Xi}$:
912 {\bf P
} =
\frac{2}{V
} (
{\bf E
}_
{kin
}-
{\bf \Xi})
915 where $V$ is the volume of the computational box.
916 The scalar pressure $P$, which can be used for pressure coupling in the case
917 of isotropic systems, is computed as:
919 P =
{\rm trace
}(
{\bf P
})/
3
922 The virial $
{\bf \Xi}$ tensor is defined as:
924 {\bf \Xi} = -
\half \sum_{i<j
} \rvij \otimes \Fvij
928 \ifthenelse{\equal{\gmxlite}{1}}{}{
929 The
{\gromacs} implementation of the virial computation is described
931 } % Brace matches ifthenelse test for gmxlite
934 \subsection{The
\swapindex{leap-frog
}{integrator
}}
935 \label{subsec:update
}
937 \centerline{\includegraphics[width=
8cm
]{plots/leapfrog
}}
938 \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
939 like frogs over each other's backs.
}
943 The default MD integrator in
{\gromacs} is the so-called
{\em leap-frog
}
944 algorithm~
\cite{Hockney74
} for the integration of the equations of
945 motion. When extremely accurate integration with temperature
946 and/or pressure coupling is required, the velocity Verlet integrators
947 are also present and may be preferable (see
\ssecref{vverlet
}). The leap-frog
948 algorithm uses positions $
\ve{r
}$ at time $t$ and
949 velocities $
\ve{v
}$ at time $t-
\hDt$; it updates positions and
950 velocities using the forces
951 $
\ve{F
}(t)$ determined by the positions at time $t$ using these relations:
953 \label{eqn:leapfrogv
}
954 \ve{v
}(t+
\hDt) &~=~&
\ve{v
}(t-
\hDt)+
\frac{\Dt}{m
}\ve{F
}(t) \\
955 \ve{r
}(t+
\Dt) &~=~&
\ve{r
}(t)+
\Dt\ve{v
}(t+
\hDt)
957 The algorithm is visualized in
\figref{leapfrog
}.
958 It produces trajectories that are identical to the Verlet~
\cite{Verlet67
} algorithm,
959 whose position-update relation is
961 \ve{r
}(t+
\Dt)~=~
2\ve{r
}(t) -
\ve{r
}(t-
\Dt) +
\frac{1}{m
}\ve{F
}(t)
\Dt^
2+O(
\Dt^
4)
963 The algorithm is of third order in $
\ve{r
}$ and is time-reversible.
964 See ref.~
\cite{Berendsen86b
} for the merits of this algorithm and comparison
965 with other time integration algorithms.
967 The
\swapindex{equations of
}{motion
} are modified for temperature
968 coupling and pressure coupling, and extended to include the
969 conservation of constraints, all of which are described below.
971 \subsection{The
\swapindex{velocity Verlet
}{integrator
}}
972 \label{subsec:vverlet
}
973 The velocity Verlet algorithm~
\cite{Swope82
} is also implemented in
974 {\gromacs}, though it is not yet fully integrated with all sets of
975 options. In velocity Verlet, positions $
\ve{r
}$ and velocities
976 $
\ve{v
}$ at time $t$ are used to integrate the equations of motion;
977 velocities at the previous half step are not required.
\bea
978 \label{eqn:velocityverlet1
}
979 \ve{v
}(t+
\hDt) &~=~&
\ve{v
}(t)+
\frac{\Dt}{2m
}\ve{F
}(t) \\
980 \ve{r
}(t+
\Dt) &~=~&
\ve{r
}(t)+
\Dt\,
\ve{v
}(t+
\hDt) \\
981 \ve{v
}(t+
\Dt) &~=~&
\ve{v
}(t+
\hDt)+
\frac{\Dt}{2m
}\ve{F
}(t+
\Dt)
985 \label{eqn:velocityverlet2
}
986 \ve{r
}(t+
\Dt) &~=~&
\ve{r
}(t)+
\Dt\,
\ve{v
} +
\frac{\Dt^
2}{2m
}\ve{F
}(t) \\
987 \ve{v
}(t+
\Dt) &~=~&
\ve{v
}(t)+
\frac{\Dt}{2m
}\left[\ve{F
}(t) +
\ve{F
}(t+
\Dt)
\right]
989 With no temperature or pressure coupling, and with
{\em corresponding
}
990 starting points, leap-frog and velocity Verlet will generate identical
991 trajectories, as can easily be verified by hand from the equations
992 above. Given a single starting file with the
{\em same
} starting
993 point $
\ve{x
}(
0)$ and $
\ve{v
}(
0)$, leap-frog and velocity Verlet will
994 {\em not
} give identical trajectories, as leap-frog will interpret the
995 velocities as corresponding to $t=-
\hDt$, while velocity Verlet will
996 interpret them as corresponding to the timepoint $t=
0$.
998 \subsection{Understanding reversible integrators: The Trotter decomposition
}
999 To further understand the relationship between velocity Verlet and
1000 leap-frog integration, we introduce the reversible Trotter formulation
1001 of dynamics, which is also useful to understanding implementations of
1002 thermostats and barostats in
{\gromacs}.
1004 A system of coupled, first-order differential equations can be evolved
1005 from time $t =
0$ to time $t$ by applying the evolution operator
1007 \Gamma(t) &=&
\exp(iLt)
\Gamma(
0)
\nonumber \\
1008 iL &=&
\dot{\Gamma}\cdot \nabla_{\Gamma},
1010 where $L$ is the Liouville operator, and $
\Gamma$ is the
1011 multidimensional vector of independent variables (positions and
1013 A short-time approximation to the true operator, accurate at time $
\Dt
1014 = t/P$, is applied $P$ times in succession to evolve the system as
1016 \Gamma(t) =
\prod_{i=
1}^P
\exp(iL
\Dt)
\Gamma(
0)
1018 For NVE dynamics, the Liouville operator is
1020 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}.
1022 This can be split into two additive operators
1024 iL_1 &=&
\sum_{i=
1}^N
\frac{1}{m_i
}\F(r_i)
\cdot \nabla_{\vv_i} \nonumber \\
1025 iL_2 &=&
\sum_{i=
1}^
{N
} \vv_i \cdot \nabla_{\rv_i}
1027 Then a short-time, symmetric, and thus reversible approximation of the true dynamics will be
1029 \exp(iL
\Dt) =
\exp(iL_2
\hDt)
\exp(iL_1
\Dt)
\exp(iL_2
\hDt) +
\mathcal{O
}(
\Dt^
3).
1030 \label{eq:NVE_Trotter
}
1032 This corresponds to velocity Verlet integration. The first
1033 exponential term over $
\hDt$ corresponds to a velocity half-step, the
1034 second exponential term over $
\Dt$ corresponds to a full velocity
1035 step, and the last exponential term over $
\hDt$ is the final velocity
1036 half step. For future times $t = n
\Dt$, this becomes
1038 \exp(iLn
\Dt) &
\approx&
\left(
\exp(iL_2
\hDt)
\exp(iL_1
\Dt)
\exp(iL_2
\hDt)
\right)^n
\nonumber \\
1039 &
\approx&
\exp(iL_2
\hDt)
\bigg(
\exp(iL_1
\Dt)
\exp(iL_2
\Dt)
\bigg)^
{n-
1} \nonumber \\
1040 & & \;\;\;\;
\exp(iL_1
\Dt)
\exp(iL_2
\hDt)
1042 This formalism allows us to easily see the difference between the
1043 different flavors of Verlet integrators. The leap-frog integrator can
1044 be seen as starting with Eq.~
\ref{eq:NVE_Trotter
} with the
1045 $
\exp\left(iL_1
\dt\right)$ term, instead of the half-step velocity
1048 \exp(iLn
\dt) &=&
\exp\left(iL_1
\dt\right)
\exp\left(iL_2
\Dt \right) +
\mathcal{O
}(
\Dt^
3).
1050 Here, the full step in velocity is between $t-
\hDt$ and $t+
\hDt$,
1051 since it is a combination of the velocity half steps in velocity
1052 Verlet. For future times $t = n
\Dt$, this becomes
1054 \exp(iLn
\dt) &
\approx&
\bigg(
\exp\left(iL_1
\dt\right)
\exp\left(iL_2
\Dt \right)
\bigg)^
{n
}.
1056 Although at first this does not appear symmetric, as long as the full velocity
1057 step is between $t-
\hDt$ and $t+
\hDt$, then this is simply a way of
1058 starting velocity Verlet at a different place in the cycle.
1060 Even though the trajectory and thus potential energies are identical
1061 between leap-frog and velocity Verlet, the kinetic energy and
1062 temperature will not necessarily be the same. Standard velocity
1063 Verlet uses the velocities at the $t$ to calculate the kinetic energy
1064 and thus the temperature only at time $t$; the kinetic energy is then a sum over all particles
1066 KE_
{\mathrm{full
}}(t) &=&
\sum_i \left(
\frac{1}{2m_i
}\ve{v
}_i(t)
\right)^
2 \nonumber\\
1067 &=&
\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,
1069 with the square on the
{\em outside
} of the average. Standard
1070 leap-frog calculates the kinetic energy at time $t$ based on the
1071 average kinetic energies at the timesteps $t+
\hDt$ and $t-
\hDt$, or
1072 the sum over all particles
1074 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),
1076 where the square is
{\em inside
} the average.
1078 A non-standard variant of velocity Verlet which averages the kinetic
1079 energies $KE(t+
\hDt)$ and $KE(t-
\hDt)$, exactly like leap-frog, is also
1080 now implemented in
{\gromacs} (as
{\tt .mdp
} file option
{\tt md-vv-avek
}). Without
1081 temperature and pressure coupling, velocity Verlet with
1082 half-step-averaged kinetic energies and leap-frog will be identical up
1083 to numerical precision. For temperature- and pressure-control schemes,
1084 however, velocity Verlet with half-step-averaged kinetic energies and
1085 leap-frog will be different, as will be discussed in the section in
1086 thermostats and barostats.
1088 The half-step-averaged kinetic energy and temperature are slightly more
1089 accurate for a given step size; the difference in average kinetic
1090 energies using the half-step-averaged kinetic energies (
{\em md
} and
1091 {\em md-vv-avek
}) will be closer to the kinetic energy obtained in the
1092 limit of small step size than will the full-step kinetic energy (using
1093 {\em md-vv
}). For NVE simulations, this difference is usually not
1094 significant, since the positions and velocities of the particles are
1095 still identical; it makes a difference in the way the the temperature
1096 of the simulations are
{\em interpreted
}, but
{\em not
} in the
1097 trajectories that are produced. Although the kinetic energy is more
1098 accurate with the half-step-averaged method, meaning that it changes
1099 less as the timestep gets large, it is also more noisy. The RMS deviation
1100 of the total energy of the system (sum of kinetic plus
1101 potential) in the half-step-averaged kinetic energy case will be
1102 higher (about twice as high in most cases) than the full-step kinetic
1103 energy. The drift will still be the same, however, as again, the
1104 trajectories are identical.
1106 For NVT simulations, however, there
{\em will
} be a difference, as
1107 discussed in the section on temperature control, since the velocities
1108 of the particles are adjusted such that kinetic energies of the
1109 simulations, which can be calculated either way, reach the
1110 distribution corresponding to the set temperature. In this case, the
1111 three methods will not give identical results.
1113 Because the velocity and position are both defined at the same time
1114 $t$ the velocity Verlet integrator can be used for some methods,
1115 especially rigorously correct pressure control methods, that are not
1116 actually possible with leap-frog. The integration itself takes
1117 negligibly more time than leap-frog, but twice as many communication
1118 calls are currently required. In most cases, and especially for large
1119 systems where communication speed is important for parallelization and
1120 differences between thermodynamic ensembles vanish in the $
1/N$ limit,
1121 and when only NVT ensembles are required, leap-frog will likely be the
1122 preferred integrator. For pressure control simulations where the fine
1123 details of the thermodynamics are important, only velocity Verlet
1124 allows the true ensemble to be calculated. In either case, simulation
1125 with double precision may be required to get fine details of
1126 thermodynamics correct.
1128 \subsection{Twin-range cut-offs
\index{twin-range!cut-off
}}
1129 To save computation time, slowly varying forces can be calculated
1130 less often than rapidly varying forces. In
{\gromacs}
1131 such a
\normindex{multiple time step
} splitting is possible between
1132 short and long range non-bonded interactions.
1133 In
{\gromacs} versions up to
4.0, an irreversible integration scheme
1134 was used which is also used by the
{\gromos} simulation package:
1135 every $n$ steps the long range forces are determined and these are
1136 then also used (without modification) for the next $n-
1$ integration steps
1137 in
\eqnref{leapfrogv
}. Such an irreversible scheme can result in bad energy
1138 conservation and, possibly, bad sampling.
1139 Since version
4.5, a leap-frog version of the reversible Trotter decomposition scheme~
\cite{Tuckerman1992a
} is used.
1140 In this integrator the long-range forces are determined every $n$ steps
1141 and are then integrated into the velocity in
\eqnref{leapfrogv
} using
1142 a time step of $
\Dt_\mathrm{LR
} = n
\Dt$:
1145 \left\
{ \begin{array
}{lll
} \displaystyle
1146 \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
1147 \ve{v
}(t-
\hDt) +
\frac{1}{m
}\ve{F
}_
\mathrm{SR
}(t)
\Dt &,&
\mathrm{step
} ~\%~ n
\neq 0 \\
1151 The parameter $n$ is equal to the neighbor list update frequency. In
1152 4.5, the velocity Verlet version of multiple time-stepping is not yet
1155 Several other simulation packages uses multiple time stepping for
1156 bonds and/or the PME mesh forces. In
{\gromacs} we have not implemented
1157 this (yet), since we use a different philosophy. Bonds can be constrained
1158 (which is also a more sound approximation of a physical quantum
1159 oscillator), which allows the smallest time step to be increased
1160 to the larger one. This not only halves the number of force calculations,
1161 but also the update calculations. For even larger time steps, angle vibrations
1162 involving hydrogen atoms can be removed using virtual interaction
1163 \ifthenelse{\equal{\gmxlite}{1}}
1165 {sites (see
\secref{rmfast
}),
}
1166 which brings the shortest time step up to
1167 PME mesh update frequency of a multiple time stepping scheme.
1169 As an example we show the energy conservation for integrating
1170 the equations of motion for SPC/E water at
300 K. To avoid cut-off
1171 effects, reaction-field electrostatics with $
\epsilon_{RF
}=
\infty$ and
1172 shifted Lennard-Jones interactions are used, both with a buffer region.
1173 The long-range interactions were evaluated between
1.0 and
1.4 nm.
1174 In
\figref{leapfrog
} one can see that for electrostatics the Trotter scheme
1175 does an order of magnitude better up to $
\Dt_{LR
}$ =
16 fs.
1176 The electrostatics depends strongly on the orientation of the water molecules,
1177 which changes rapidly.
1178 For Lennard-Jones interactions, the energy drift is linear in $
\Dt_{LR
}$
1179 and roughly two orders of magnitude smaller than for the electrostatics.
1180 Lennard-Jones forces are smaller than Coulomb forces and
1181 they are mainly affected by translation of water molecules, not rotation.
1184 \centerline{\includegraphics[width=
12cm
]{plots/drift-all
}}
1185 \caption{Energy drift per degree of freedom in SPC/E water
1186 with twin-range cut-offs
1187 for reaction field (left) and Lennard-Jones interaction (right)
1188 as a function of the long-range time step length for the irreversible
1189 ``
\gromos'' scheme and a reversible Trotter scheme.
}
1190 \label{fig:twinrangeener
}
1193 \subsection{Temperature coupling
\index{temperature coupling
}}
1194 While direct use of molecular dynamics gives rise to the NVE (constant
1195 number, constant volume, constant energy ensemble), most quantities
1196 that we wish to calculate are actually from a constant temperature
1197 (NVT) ensemble, also called the canonical ensemble.
{\gromacs} can use
1198 the
{\em weak-coupling
} scheme of Berendsen~
\cite{Berendsen84
},
1199 stochastic randomization through the Andersen
1200 thermostat~
\cite{Andersen80
}, the extended ensemble Nos
{\'e
}-Hoover
1201 scheme~
\cite{Nose84,Hoover85
}, or a velocity-rescaling
1202 scheme~
\cite{Bussi2007a
} to simulate constant temperature, with
1203 advantages of each of the schemes laid out below.
1205 There are several other reasons why it might be necessary to control
1206 the temperature of the system (drift during equilibration, drift as a
1207 result of force truncation and integration errors, heating due to
1208 external or frictional forces), but this is not entirely correct to do
1209 from a thermodynamic standpoint, and in some cases only masks the
1210 symptoms (increase in temperature of the system) rather than the
1211 underlying problem (deviations from correct physics in the dynamics).
1212 For larger systems, errors in ensemble averages and structural
1213 properties incurred by using temperature control to remove slow drifts
1214 in temperature appear to be negligible, but no completely
1215 comprehensive comparisons have been carried out, and some caution must
1216 be taking in interpreting the results.
1218 \subsubsection{Berendsen temperature coupling
\pawsindexquiet{Berendsen
}{temperature coupling
}\index{weak coupling
}}
1219 The Berendsen algorithm mimics weak coupling with first-order
1220 kinetics to an external heat bath with given temperature $T_0$.
1221 See ref.~
\cite{Berendsen91
} for a comparison with the
1222 Nos
{\'e
}-Hoover scheme. The effect of this algorithm is
1223 that a deviation of the system temperature from $T_0$ is slowly
1224 corrected according to:
1226 \frac{\de T
}{\de t
} =
\frac{T_0-T
}{\tau}
1227 \label{eqn:Tcoupling
}
1229 which means that a temperature deviation decays exponentially with a
1230 time constant $
\tau$.
1231 This method of coupling has the advantage that the strength of the
1232 coupling can be varied and adapted to the user requirement: for
1233 equilibration purposes the coupling time can be taken quite short
1234 (
{\eg} 0.01 ps), but for reliable equilibrium runs it can be taken much
1235 longer (
{\eg} 0.5 ps) in which case it hardly influences the
1236 conservative dynamics.
1238 The Berendsen thermostat suppresses the fluctuations of the kinetic
1239 energy. This means that one does not generate a proper canonical
1240 ensemble, so rigorously, the sampling will be incorrect. This
1241 error scales with $
1/N$, so for very large systems most ensemble
1242 averages will not be affected significantly, except for the
1243 distribution of the kinetic energy itself. However, fluctuation
1244 properties, such as the heat capacity, will be affected. A similar
1245 thermostat which does produce a correct ensemble is the velocity
1246 rescaling thermostat~
\cite{Bussi2007a
} described below.
1248 The heat flow into or out of the system is affected by scaling the
1249 velocities of each particle every step, or every $n_
\mathrm{TC
}$ steps,
1250 with a time-dependent factor $
\lambda$, given by:
1252 \lambda =
\left[ 1 +
\frac{n_
\mathrm{TC
} \Delta t
}{\tau_T}
1253 \left\
{\frac{T_0
}{T(t -
\hDt)
} -
1 \right\
} \right]^
{1/
2}
1256 The parameter $
\tau_T$ is close, but not exactly equal, to the time constant
1257 $
\tau$ of the temperature coupling (
\eqnref{Tcoupling
}):
1259 \tau =
2 C_V
\tau_T / N_
{df
} k
1261 where $C_V$ is the total heat capacity of the system, $k$ is Boltzmann's
1262 constant, and $N_
{df
}$ is the total number of degrees of freedom. The
1263 reason that $
\tau \neq \tau_T$ is that the kinetic energy change
1264 caused by scaling the velocities is partly redistributed between
1265 kinetic and potential energy and hence the change in temperature is
1266 less than the scaling energy. In practice, the ratio $
\tau /
\tau_T$
1267 ranges from
1 (gas) to
2 (harmonic solid) to
3 (water). When we use
1268 the term ``temperature coupling time constant,'' we mean the parameter
1269 \normindex{$
\tau_T$
}.
1270 {\bf Note
} that in practice the scaling factor $
\lambda$ is limited to
1271 the range of
0.8 $<=
\lambda <=$
1.25, to avoid scaling by very large
1272 numbers which may crash the simulation. In normal use,
1273 $
\lambda$ will always be much closer to
1.0.
1275 \subsubsection{Velocity-rescaling temperature coupling
\pawsindexquiet{velocity-rescaling
}{temperature coupling
}}
1276 The velocity-rescaling thermostat~
\cite{Bussi2007a
} is essentially a Berendsen
1277 thermostat (see above) with an additional stochastic term that ensures
1278 a correct kinetic energy distribution by modifying it according to
1280 \de K = (K_0 - K)
\frac{\de t
}{\tau_T} +
2 \sqrt{\frac{K K_0
}{N_f
}} \frac{\de W
}{\sqrt{\tau_T}},
1281 \label{eqn:vrescale
}
1283 where $K$ is the kinetic energy, $N_f$ the number of degrees of freedom and $
\de W$ a Wiener process.
1284 There are no additional parameters, except for a random seed.
1285 This thermostat produces a correct canonical ensemble and still has
1286 the advantage of the Berendsen thermostat: first order decay of
1287 temperature deviations and no oscillations.
1288 When an $NVT$ ensemble is used, the conserved energy quantity
1289 is written to the energy and log file.
1291 \subsubsection{\normindex{Andersen thermostat
}}
1292 One simple way to maintain a thermostatted ensemble is to take an
1293 $NVE$ integrator and periodically re-select the velocities of the
1294 particles from a Maxwell-Boltzmann distribution.~
\cite{Andersen80
}
1295 This can either be done by randomizing all the velocities
1296 simultaneously (massive collision) every $
\tau_T/
\Dt$ steps (
{\tt andersen-massive
}), or by
1297 randomizing every particle with some small probability every timestep (
{\tt andersen
}),
1298 equal to $
\Dt/
\tau$, where in both cases $
\Dt$ is the timestep and
1299 $
\tau_T$ is a characteristic coupling time scale.
1300 Because of the way constraints operate, all particles in the same
1301 constraint group must be randomized simultaneously. Because of
1302 parallelization issues, the
{\tt andersen
} version cannot currently (
5.0) be
1303 used in systems with constraints.
{\tt andersen-massive
} can be used regardless of constraints.
1304 This thermostat is also currently only possible with velocity Verlet algorithms,
1305 because it operates directly on the velocities at each timestep.
1307 This algorithm completely avoids some of the ergodicity issues of other thermostatting
1308 algorithms, as energy cannot flow back and forth between energetically
1309 decoupled components of the system as in velocity scaling motions.
1310 However, it can slow down the kinetics of system by randomizing
1311 correlated motions of the system, including slowing sampling when
1312 $
\tau_T$ is at moderate levels (less than
10 ps). This algorithm
1313 should therefore generally not be used when examining kinetics or
1314 transport properties of the system.~
\cite{Basconi2013
}
1316 % \ifthenelse{\equal{\gmxlite}{1}}{}{
1317 \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
}}
1319 The Berendsen weak-coupling algorithm is
1320 extremely efficient for relaxing a system to the target temperature,
1321 but once the system has reached equilibrium it might be more
1322 important to probe a correct canonical ensemble. This is unfortunately
1323 not the case for the weak-coupling scheme.
1325 To enable canonical ensemble simulations,
{\gromacs} also supports the
1326 extended-ensemble approach first proposed by Nos
{\'e
}~
\cite{Nose84
}
1327 and later modified by Hoover~
\cite{Hoover85
}. The system Hamiltonian is
1328 extended by introducing a thermal reservoir and a friction term in the
1329 equations of motion. The friction force is proportional to the
1330 product of each particle's velocity and a friction parameter, $
\xi$.
1331 This friction parameter (or ``heat bath'' variable) is a fully
1332 dynamic quantity with its own momentum ($p_
{\xi}$) and equation of
1333 motion; the time derivative is calculated from the difference between
1334 the current kinetic energy and the reference temperature.
1336 In this formulation, the particles' equations of motion in
1337 \figref{global
} are replaced by:
1339 \frac {\de^
2\ve{r
}_i
}{\de t^
2} =
\frac{\ve{F
}_i
}{m_i
} -
1340 \frac{p_
{\xi}}{Q
}\frac{\de \ve{r
}_i
}{\de t
} ,
1341 \label{eqn:NH-eqn-of-motion
}
1342 \eeq where the equation of motion for the heat bath parameter $
\xi$ is:
1343 \beq \frac {\de p_
{\xi}}{\de t
} =
\left( T - T_0
\right).
\eeq The
1344 reference temperature is denoted $T_0$, while $T$ is the current
1345 instantaneous temperature of the system. The strength of the coupling
1346 is determined by the constant $Q$ (usually called the ``mass parameter''
1347 of the reservoir) in combination with the reference
1348 temperature.~
\footnote{Note that some derivations, an alternative
1349 notation $
\xi_{\mathrm{alt
}} = v_
{\xi} = p_
{\xi}/Q$ is used.
}
1351 The conserved quantity for the Nos
{\'e
}-Hoover equations of motion is not
1352 the total energy, but rather
1354 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,
1356 where $N_f$ is the total number of degrees of freedom.
1358 In our opinion, the mass parameter is a somewhat awkward way of
1359 describing coupling strength, especially due to its dependence on
1360 reference temperature (and some implementations even include the
1361 number of degrees of freedom in your system when defining $Q$). To
1362 maintain the coupling strength, one would have to change $Q$ in
1363 proportion to the change in reference temperature. For this reason, we
1364 prefer to let the
{\gromacs} user work instead with the period
1365 $
\tau_T$ of the oscillations of kinetic energy between the system and
1366 the reservoir instead. It is directly related to $Q$ and $T_0$ via:
1368 Q =
\frac {\tau_T^
2 T_0
}{4 \pi^
2}.
1370 This provides a much more intuitive way of selecting the
1371 Nos
{\'e
}-Hoover coupling strength (similar to the weak-coupling
1372 relaxation), and in addition $
\tau_T$ is independent of system size
1373 and reference temperature.
1375 It is however important to keep the difference between the
1376 weak-coupling scheme and the Nos
{\'e
}-Hoover algorithm in mind:
1377 Using weak coupling you get a
1378 strongly damped
{\em exponential relaxation
},
1379 while the Nos
{\'e
}-Hoover approach
1380 produces an
{\em oscillatory relaxation
}.
1381 The actual time it takes to relax with Nos
{\'e
}-Hoover coupling is
1382 several times larger than the period of the
1383 oscillations that you select. These oscillations (in contrast
1384 to exponential relaxation) also means that
1385 the time constant normally should be
4--
5 times larger
1386 than the relaxation time used with weak coupling, but your
1389 Nos
{\'e
}-Hoover dynamics in simple systems such as collections of
1390 harmonic oscillators, can be
{\em nonergodic
}, meaning that only a
1391 subsection of phase space is ever sampled, even if the simulations
1392 were to run for infinitely long. For this reason, the Nos
{\'e
}-Hoover
1393 chain approach was developed, where each of the Nos
{\'e
}-Hoover
1394 thermostats has its own Nos
{\'e
}-Hoover thermostat controlling its
1395 temperature. In the limit of an infinite chain of thermostats, the
1396 dynamics are guaranteed to be ergodic. Using just a few chains can
1397 greatly improve the ergodicity, but recent research has shown that the
1398 system will still be nonergodic, and it is still not entirely clear
1399 what the practical effect of this~
\cite{Cooke2008
}. Currently, the
1400 default number of chains is
10, but this can be controlled by the
1401 user. In the case of chains, the equations are modified in the
1402 following way to include a chain of thermostatting
1403 particles~
\cite{Martyna1992
}:
1406 \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 \\
1407 \frac {\de p_
{{\xi}_1
}}{\de t
} &~=~&
\left( T - T_0
\right) - p_
{{\xi}_1
} \frac{p_
{{\xi}_2
}}{Q_2
} \nonumber \\
1408 \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 \\
1409 \frac {\de p_
{\xi_N}}{\de t
} &~=~&
\left(
\frac{p_
{\xi_{N-
1}}^
2}{Q_
{N-
1}}-kT
\right)
1410 \label{eqn:NH-chain-eqn-of-motion
}
1412 The conserved quantity for Nos
{\'e
}-Hoover chains is
1414 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
1416 The values and velocities of the Nos
{\'e
}-Hoover thermostat variables
1417 are generally not included in the output, as they take up a fair
1418 amount of space and are generally not important for analysis of
1419 simulations, but this can be overridden by defining the environment
1420 variable
{\tt GMX_NOSEHOOVER_CHAINS
}, which will print the values of all
1421 the positions and velocities of all Nos
{\'e
}-Hoover particles in the
1422 chain to the
{\tt .edr
} file. Leap-frog simulations currently can only have
1423 Nos
{\'e
}-Hoover chain lengths of
1, but this will likely be updated in
1426 As described in the integrator section, for temperature coupling, the
1427 temperature that the algorithm attempts to match to the reference
1428 temperature is calculated differently in velocity Verlet and leap-frog
1429 dynamics. Velocity Verlet (
{\em md-vv
}) uses the full-step kinetic
1430 energy, while leap-frog and
{\em md-vv-avek
} use the half-step-averaged
1433 We can examine the Trotter decomposition again to better understand
1434 the differences between these constant-temperature integrators. In
1435 the case of Nos
{\'e
}-Hoover dynamics (for simplicity, using a chain
1436 with $N=
1$, with more details in Ref.~
\cite{Martyna1996
}), we split
1437 the Liouville operator as
1439 iL = iL_1 + iL_2 + iL_
{\mathrm{NHC
}},
1443 iL_1 &=&
\sum_{i=
1}^N
\left[\frac{\pb_i}{m_i
}\right]\cdot \frac{\partial}{\partial \rv_i} \nonumber \\
1444 iL_2 &=&
\sum_{i=
1}^N
\F_i\cdot \frac{\partial}{\partial \pb_i} \nonumber \\
1445 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}}
1447 For standard velocity Verlet with Nos
{\'e
}-Hoover temperature control, this becomes
1449 \exp(iL
\dt) &=&
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right)
\exp\left(iL_2
\dt/
2\right)
\nonumber \\
1450 &&
\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).
1452 For half-step-averaged temperature control using
{\em md-vv-avek
},
1453 this decomposition will not work, since we do not have the full step
1454 temperature until after the second velocity step. However, we can
1455 construct an alternate decomposition that is still reversible, by
1456 switching the place of the NHC and velocity portions of the
1459 \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 \\
1460 &&
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right)
\exp\left(iL_2
\dt/
2\right)+
\mathcal{O
}(
\Dt^
3)
1461 \label{eq:half_step_NHC_integrator
}
1463 This formalism allows us to easily see the difference between the
1464 different flavors of velocity Verlet integrator. The leap-frog
1465 integrator can be seen as starting with
1466 Eq.~
\ref{eq:half_step_NHC_integrator
} just before the $
\exp\left(iL_1
1467 \dt\right)$ term, yielding:
1469 \exp(iL
\dt) &=&
\exp\left(iL_1
\dt\right)
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right)
\nonumber \\
1470 &&
\exp\left(iL_2
\dt\right)
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right) +
\mathcal{O
}(
\Dt^
3)
1472 and then using some algebra tricks to solve for some quantities are
1473 required before they are actually calculated~
\cite{Holian95
}.
1477 \subsubsection{Group temperature coupling
}\index{temperature-coupling group
}%
1478 In
{\gromacs} temperature coupling can be performed on groups of
1479 atoms, typically a protein and solvent. The reason such algorithms
1480 were introduced is that energy exchange between different components
1481 is not perfect, due to different effects including cut-offs etc. If
1482 now the whole system is coupled to one heat bath, water (which
1483 experiences the largest cut-off noise) will tend to heat up and the
1484 protein will cool down. Typically
100 K differences can be obtained.
1485 With the use of proper electrostatic methods (PME) these difference
1486 are much smaller but still not negligible. The parameters for
1487 temperature coupling in groups are given in the
{\tt mdp
} file.
1488 Recent investigation has shown that small temperature differences
1489 between protein and water may actually be an artifact of the way
1490 temperature is calculated when there are finite timesteps, and very
1491 large differences in temperature are likely a sign of something else
1492 seriously going wrong with the system, and should be investigated
1493 carefully~
\cite{Eastwood2010
}.
1495 One special case should be mentioned: it is possible to temperature-couple only
1496 part of the system, leaving other parts without temperature
1497 coupling. This is done by specifying $
{-
1}$ for the time constant
1498 $
\tau_T$ for the group that should not be thermostatted. If only
1499 part of the system is thermostatted, the system will still eventually
1500 converge to an NVT system. In fact, one suggestion for minimizing
1501 errors in the temperature caused by discretized timesteps is that if
1502 constraints on the water are used, then only the water degrees of
1503 freedom should be thermostatted, not protein degrees of freedom, as
1504 the higher frequency modes in the protein can cause larger deviations
1505 from the ``true'' temperature, the temperature obtained with small
1506 timesteps~
\cite{Eastwood2010
}.
1508 \subsection{Pressure coupling
\index{pressure coupling
}}
1509 In the same spirit as the temperature coupling, the system can also be
1510 coupled to a ``pressure bath.''
{\gromacs} supports both the Berendsen
1511 algorithm~
\cite{Berendsen84
} that scales coordinates and box vectors
1512 every step, the extended-ensemble Parrinello-Rahman approach~
\cite{Parrinello81,Nose83
}, and for
1513 the velocity Verlet variants, the Martyna-Tuckerman-Tobias-Klein
1514 (MTTK) implementation of pressure
1515 control~
\cite{Martyna1996
}. Parrinello-Rahman and Berendsen can be
1516 combined with any of the temperature coupling methods above; MTTK can
1517 only be used with Nos
{\'e
}-Hoover temperature control.
1519 \subsubsection{Berendsen pressure coupling
\pawsindexquiet{Berendsen
}{pressure coupling
}\index{weak coupling
}}
1520 \label{sec:berendsen_pressure_coupling
}
1521 The Berendsen algorithm rescales the
1522 coordinates and box vectors every step, or every $n_
\mathrm{PC
}$ steps,
1523 with a matrix
{\boldmath $
\mu$
},
1524 which has the effect of a first-order kinetic relaxation of the pressure
1525 towards a given reference pressure $
{\bf P
}_0$ according to
1527 \frac{\de {\bf P
}}{\de t
} =
\frac{{\bf P
}_0-
{\bf P
}}{\tau_p}.
1529 The scaling matrix
{\boldmath $
\mu$
} is given by
1532 =
\delta_{ij
} -
\frac{n_
\mathrm{PC
}\Delta t
}{3\,
\tau_p} \beta_{ij
} \
{P_
{0ij
} - P_
{ij
}(t) \
}.
1535 \index{isothermal compressibility
}
1536 \index{compressibility
}
1537 Here,
{\boldmath $
\beta$
} is the isothermal compressibility of the system.
1538 In most cases this will be a diagonal matrix, with equal elements on the
1539 diagonal, the value of which is generally not known.
1540 It suffices to take a rough estimate because the value of
{\boldmath $
\beta$
}
1541 only influences the non-critical time constant of the
1542 pressure relaxation without affecting the average pressure itself.
1543 For water at
1 atm and
300 K
1544 $
\beta =
4.6 \times 10^
{-
10}$ Pa$^
{-
1} =
4.6 \times 10^
{-
5}$ bar$^
{-
1}$,
1545 which is $
7.6 \times 10^
{-
4}$ MD units (see
\chref{defunits
}).
1546 Most other liquids have similar values.
1547 When scaling completely anisotropically, the system has to be rotated in
1548 order to obey
\eqnref{box_rot
}.
1549 This rotation is approximated in first order in the scaling, which is usually
1550 less than $
10^
{-
4}$. The actual scaling matrix
{\boldmath $
\mu'$
} is
1552 \mbox{\boldmath $
\mu'$
} =
1553 \left(
\begin{array
}{ccc
}
1554 \mu_{xx
} &
\mu_{xy
} +
\mu_{yx
} &
\mu_{xz
} +
\mu_{zx
} \\
1555 0 &
\mu_{yy
} &
\mu_{yz
} +
\mu_{zy
} \\
1559 The velocities are neither scaled nor rotated.
1561 In
{\gromacs}, the Berendsen scaling can also be done isotropically,
1562 which means that instead of $
\ve{P
}$ a diagonal matrix with elements of size
1563 trace$(
\ve{P
})/
3$ is used. For systems with interfaces, semi-isotropic
1564 scaling can be useful.
1565 In this case, the $x/y$-directions are scaled isotropically and the $z$
1566 direction is scaled independently. The compressibility in the $x/y$ or
1567 $z$-direction can be set to zero, to scale only in the other direction(s).
1569 If you allow full anisotropic deformations and use constraints you
1570 might have to scale more slowly or decrease your timestep to avoid
1571 errors from the constraint algorithms. It is important to note that
1572 although the Berendsen pressure control algorithm yields a simulation
1573 with the correct average pressure, it does not yield the exact NPT
1574 ensemble, and it is not yet clear exactly what errors this approximation
1577 % \ifthenelse{\equal{\gmxlite}{1}}{}{
1578 \subsubsection{Parrinello-Rahman pressure coupling
\pawsindexquiet{Parrinello-Rahman
}{pressure coupling
}}
1580 In cases where the fluctuations in pressure or volume are important
1581 {\em per se
} (
{\eg} to calculate thermodynamic properties), especially
1582 for small systems, it may be a problem that the exact ensemble is not
1583 well defined for the weak-coupling scheme, and that it does not
1584 simulate the true NPT ensemble.
1586 {\gromacs} also supports constant-pressure simulations using the
1587 Parrinello-Rahman approach~
\cite{Parrinello81,Nose83
}, which is similar
1588 to the Nos
{\'e
}-Hoover temperature coupling, and in theory gives the
1589 true NPT ensemble. With the Parrinello-Rahman barostat, the box
1590 vectors as represented by the matrix
\ve{b
} obey the matrix equation
1591 of motion
\footnote{The box matrix representation
\ve{b
} in
{\gromacs}
1592 corresponds to the transpose of the box matrix representation
\ve{h
}
1593 in the paper by Nos
{\'e
} and Klein. Because of this, some of our
1594 equations will look slightly different.
}
1596 \frac{\de \ve{b
}^
2}{\de t^
2}= V
\ve{W
}^
{-
1} \ve{b
}'^
{-
1} \left(
\ve{P
} -
\ve{P
}_
{ref
}\right).
1599 The volume of the box is denoted $V$, and $
\ve{W
}$ is a matrix parameter that determines
1600 the strength of the coupling. The matrices
\ve{P
} and
\ve{P
}$_
{ref
}$ are the
1601 current and reference pressures, respectively.
1603 The equations of motion for the particles are also changed, just as
1604 for the Nos
{\'e
}-Hoover coupling. In most cases you would combine the
1605 Parrinello-Rahman barostat with the Nos
{\'e
}-Hoover
1606 thermostat, but to keep it simple we only show the Parrinello-Rahman
1609 \bea \frac {\de^
2\ve{r
}_i
}{\de t^
2} & = &
\frac{\ve{F
}_i
}{m_i
} -
1610 \ve{M
} \frac{\de \ve{r
}_i
}{\de t
} , \\
\ve{M
} & = &
\ve{b
}^
{-
1} \left[
1611 \ve{b
} \frac{\de \ve{b
}'
}{\de t
} +
\frac{\de \ve{b
}}{\de t
} \ve{b
}'
1612 \right] \ve{b
}'^
{-
1}.
\eea The (inverse) mass parameter matrix
1613 $
\ve{W
}^
{-
1}$ determines the strength of the coupling, and how the box
1614 can be deformed. The box restriction (
\ref{eqn:box_rot
}) will be
1615 fulfilled automatically if the corresponding elements of $
\ve{W
}^
{-
1}$
1616 are zero. Since the coupling strength also depends on the size of your
1617 box, we prefer to calculate it automatically in
{\gromacs}. You only
1618 have to provide the approximate isothermal compressibilities
1619 {\boldmath $
\beta$
} and the pressure time constant $
\tau_p$ in the
1620 input file ($L$ is the largest box matrix element):
\beq \left(
1621 \ve{W
}^
{-
1} \right)_
{ij
} =
\frac{4 \pi^
2 \beta_{ij
}}{3 \tau_p^
2 L
}.
1622 \eeq Just as for the Nos
{\'e
}-Hoover thermostat, you should realize
1623 that the Parrinello-Rahman time constant is
{\em not
} equivalent to
1624 the relaxation time used in the Berendsen pressure coupling algorithm.
1625 In most cases you will need to use a
4--
5 times larger time constant
1626 with Parrinello-Rahman coupling. If your pressure is very far from
1627 equilibrium, the Parrinello-Rahman coupling may result in very large
1628 box oscillations that could even crash your run. In that case you
1629 would have to increase the time constant, or (better) use the weak-coupling
1630 scheme to reach the target pressure, and then switch to
1631 Parrinello-Rahman coupling once the system is in equilibrium.
1632 Additionally, using the leap-frog algorithm, the pressure at time $t$
1633 is not available until after the time step has completed, and so the
1634 pressure from the previous step must be used, which makes the algorithm
1635 not directly reversible, and may not be appropriate for high precision
1636 thermodynamic calculations.
1638 \subsubsection{Surface-tension coupling
\pawsindexquiet{surface-tension
}{pressure coupling
}}
1639 When a periodic system consists of more than one phase, separated by
1640 surfaces which are parallel to the $xy$-plane,
1641 the surface tension and the $z$-component of the pressure can be coupled
1642 to a pressure bath. Presently, this only works with the Berendsen
1643 pressure coupling algorithm in
{\gromacs}.
1644 The average surface tension $
\gamma(t)$ can be calculated from
1645 the difference between the normal and the lateral pressure
1648 \frac{1}{n
} \int_0^
{L_z
}
1649 \left\
{ P_
{zz
}(z,t) -
\frac{P_
{xx
}(z,t) + P_
{yy
}(z,t)
}{2} \right\
} \mbox{d
}z \\
1651 \frac{L_z
}{n
} \left\
{ P_
{zz
}(t) -
\frac{P_
{xx
}(t) + P_
{yy
}(t)
}{2} \right\
},
1653 where $L_z$ is the height of the box and $n$ is the number of surfaces.
1654 The pressure in the z-direction is corrected by scaling the height of
1655 the box with $
\mu_{zz
}$
1657 \Delta P_
{zz
} =
\frac{\Delta t
}{\tau_p} \
{ P_
{0zz
} - P_
{zz
}(t) \
}
1660 \mu_{zz
} =
1 +
\beta_{zz
} \Delta P_
{zz
}
1662 This is similar to normal pressure coupling, except that the factor
1663 of $
1/
3$ is missing.
1664 The pressure correction in the $z$-direction is then used to get the
1665 correct convergence for the surface tension to the reference value $
\gamma_0$.
1666 The correction factor for the box length in the $x$/$y$-direction is
1668 \mu_{x/y
} =
1 +
\frac{\Delta t
}{2\,
\tau_p} \beta_{x/y
}
1669 \left(
\frac{n
\gamma_0}{\mu_{zz
} L_z
}
1670 -
\left\
{ P_
{zz
}(t)+
\Delta P_
{zz
} -
\frac{P_
{xx
}(t) + P_
{yy
}(t)
}{2} \right\
}
1673 The value of $
\beta_{zz
}$ is more critical than with normal pressure
1674 coupling. Normally an incorrect compressibility will just scale $
\tau_p$,
1675 but with surface tension coupling it affects the convergence of the surface
1677 When $
\beta_{zz
}$ is set to zero (constant box height), $
\Delta P_
{zz
}$ is also set
1678 to zero, which is necessary for obtaining the correct surface tension.
1680 \subsubsection{MTTK pressure control algorithms
}
1682 As mentioned in the previous section, one weakness of leap-frog
1683 integration is in constant pressure simulations, since the pressure
1684 requires a calculation of both the virial and the kinetic energy at
1685 the full time step; for leap-frog, this information is not available
1686 until
{\em after
} the full timestep. Velocity Verlet does allow the
1687 calculation, at the cost of an extra round of global communication,
1688 and can compute, mod any integration errors, the true NPT ensemble.
1690 The full equations, combining both pressure coupling and temperature
1691 coupling, are taken from Martyna
{\em et al.
}~
\cite{Martyna1996
} and
1692 Tuckerman~
\cite{Tuckerman2006
} and are referred to here as MTTK
1693 equations (Martyna-Tuckerman-Tobias-Klein). We introduce for
1694 convenience $
\epsilon = (
1/
3)
\ln (V/V_0)$, where $V_0$ is a reference
1695 volume. The momentum of $
\epsilon$ is $
\veps = p_
{\epsilon}/W =
1696 \dot{\epsilon} =
\dot{V
}/
3V$, and define $
\alpha =
1 +
3/N_
{dof
}$ (see
1697 Ref~
\cite{Tuckerman2006
})
1699 The isobaric equations are
1701 \dot{\rv}_i &=&
\frac{\pb_i}{m_i
} +
\frac{\peps}{W
} \rv_i \nonumber \\
1702 \frac{\dot{\pb}_i
}{m_i
} &=&
\frac{1}{m_i
}\F_i -
\alpha\frac{\peps}{W
} \frac{\pb_i}{m_i
} \nonumber \\
1703 \dot{\epsilon} &=&
\frac{\peps}{W
} \nonumber \\
1704 \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),\\
1708 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\
1711 The terms including $
\alpha$ are required to make phase space
1712 incompressible~
\cite{Tuckerman2006
}. The $
\epsilon$ acceleration term
1715 \frac{\dot{\peps}}{W
} &=&
\frac{3V
}{W
}\left(
\alpha P_
{\mathrm{kin
}} - P_
{\mathrm{vir
}} - P
\right)
1717 In terms of velocities, these equations become
1719 \dot{\rv}_i &=&
\vv_i +
\veps \rv_i \nonumber \\
1720 \dot{\vv}_i &=&
\frac{1}{m_i
}\F_i -
\alpha\veps \vv_i \nonumber \\
1721 \dot{\epsilon} &=&
\veps \nonumber \\
1722 \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 \\
1723 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]
1725 For these equations, the conserved quantity is
1727 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
1729 The next step is to add temperature control. Adding Nos
{\'e
}-Hoover
1730 chains, including to the barostat degree of freedom, where we use
1731 $
\eta$ for the barostat Nos
{\'e
}-Hoover variables, and $Q^
{\prime}$
1732 for the coupling constants of the thermostats of the barostats, we get
1734 \dot{\rv}_i &=&
\frac{\pb_i}{m_i
} +
\frac{\peps}{W
} \rv_i \nonumber \\
1735 \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 \\
1736 \dot{\epsilon} &=&
\frac{\peps}{W
} \nonumber \\
1737 \frac{\dot{\peps}}{W
} &=&
\frac{3V
}{W
}(
\alpha P_
{\mathrm{kin
}} - P_
{\mathrm{vir
}} - P) -
\frac{p_
{\eta_1}}{Q^
{\prime}_1
}\peps \nonumber \\
1738 \dot{\xi}_k &=&
\frac{p_
{\xi_k}}{Q_k
} \nonumber \\
1739 \dot{\eta}_k &=&
\frac{p_
{\eta_k}}{Q^
{\prime}_k
} \nonumber \\
1740 \dot{p
}_
{\xi_k} &=& G_k -
\frac{p_
{\xi_{k+
1}}}{Q_
{k+
1}} \;\;\;\; k=
1,
\ldots, M-
1 \nonumber \\
1741 \dot{p
}_
{\eta_k} &=& G^
\prime_k -
\frac{p_
{\eta_{k+
1}}}{Q^
\prime_{k+
1}} \;\;\;\; k=
1,
\ldots, M-
1 \nonumber \\
1742 \dot{p
}_
{\xi_M} &=& G_M
\nonumber \\
1743 \dot{p
}_
{\eta_M} &=& G^
\prime_M,
\nonumber \\
1747 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 \\
1748 G_1 &=&
\sum_{i=
1}^N
\frac{\pb^
2_i
}{m_i
} - N_f kT
\nonumber \\
1749 G_k &=&
\frac{p^
2_
{\xi_{k-
1}}}{2Q_
{k-
1}} - kT \;\; k =
2,
\ldots,M
\nonumber \\
1750 G^
\prime_1 &=&
\frac{\peps^
2}{2W
} - kT
\nonumber \\
1751 G^
\prime_k &=&
\frac{p^
2_
{\eta_{k-
1}}}{2Q^
\prime_{k-
1}} - kT \;\; k =
2,
\ldots,M
1753 The conserved quantity is now
1755 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 \\
1756 \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
1758 Returning to the Trotter decomposition formalism, for pressure control and temperature control~
\cite{Martyna1996
} we get:
1760 iL = iL_1 + iL_2 + iL_
{\epsilon,
1} + iL_
{\epsilon,
2} + iL_
{\mathrm{NHC-baro
}} + iL_
{\mathrm{NHC
}}
1762 where ``NHC-baro'' corresponds to the Nos
{\`e
}-Hoover chain of the barostat,
1763 and NHC corresponds to the NHC of the particles,
1765 iL_1 &=&
\sum_{i=
1}^N
\left[\frac{\pb_i}{m_i
} +
\frac{\peps}{W
}\rv_i\right]\cdot \frac{\partial}{\partial \rv_i} \\
1766 iL_2 &=&
\sum_{i=
1}^N
\F_i -
\alpha \frac{\peps}{W
}\pb_i \cdot \frac{\partial}{\partial \pb_i} \\
1767 iL_
{\epsilon,
1} &=&
\frac{p_
\epsilon}{W
} \frac{\partial}{\partial \epsilon}\\
1768 iL_
{\epsilon,
2} &=& G_
{\epsilon} \frac{\partial}{\partial p_
\epsilon}
1772 G_
{\epsilon} =
3V
\left(
\alpha P_
{\mathrm{kin
}} - P_
{\mathrm{vir
}} - P
\right)
1774 Using the Trotter decomposition, we get
1776 \exp(iL
\dt) &=&
\exp\left(iL_
{\mathrm{NHC-baro
}}\dt/
2\right)
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right)
\nonumber \nonumber \\
1777 &&
\exp\left(iL_
{\epsilon,
2}\dt/
2\right)
\exp\left(iL_2
\dt/
2\right)
\nonumber \nonumber \\
1778 &&
\exp\left(iL_
{\epsilon,
1}\dt\right)
\exp\left(iL_1
\dt\right)
\nonumber \nonumber \\
1779 &&
\exp\left(iL_2
\dt/
2\right)
\exp\left(iL_
{\epsilon,
2}\dt/
2\right)
\nonumber \nonumber \\
1780 &&
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right)
\exp\left(iL_
{\mathrm{NHC-baro
}}\dt/
2\right) +
\mathcal{O
}(
\dt^
3)
1782 The action of $
\exp\left(iL_1
\dt\right)$ comes from the solution of
1783 the the differential equation
1784 $
\dot{\rv}_i =
\vv_i +
\veps \rv_i$
1785 with $
\vv_i =
\pb_i/m_i$ and $
\veps$ constant with initial condition
1786 $
\rv_i(
0)$, evaluate at $t=
\Delta t$. This yields the evolution
1788 \rv_i(
\dt) =
\rv_i(
0)e^
{\veps \dt} +
\Delta t
\vv_i(
0) e^
{\veps \dt/
2} \sinhx{\veps \dt/
2}.
1790 The action of $
\exp\left(iL_2
\dt/
2\right)$ comes from the solution
1791 of the differential equation $
\dot{\vv}_i =
\frac{\F_i}{m_i
} -
1792 \alpha\veps\vv_i$, yielding
1794 \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}.
1796 {\em md-vv-avek
} uses the full step kinetic energies for determining the pressure with the pressure control,
1797 but the half-step-averaged kinetic energy for the temperatures, which can be written as a Trotter decomposition as
1799 \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 \\
1800 &&
\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 \\
1801 &&
\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)
1803 With constraints, the equations become significantly more
1804 complicated, in that each of these equations need to be solved
1805 iteratively for the constraint forces. The discussion of the details of the iteration
1806 is beyond the scope of this manual; readers are encouraged to see the
1807 implementation described in~
\cite{Yu2010
}.
1810 \subsubsection{Infrequent evaluation of temperature and pressure coupling
}
1812 Temperature and pressure control require global communication to
1813 compute the kinetic energy and virial, which can become costly if
1814 performed every step for large systems. We can rearrange the Trotter
1815 decomposition to give alternate symplectic, reversible integrator with
1816 the coupling steps every $n$ steps instead of every steps. These new
1817 integrators will diverge if the coupling time step is too large, as
1818 the auxiliary variable integrations will not converge. However, in
1819 most cases, long coupling times are more appropriate, as they disturb
1820 the dynamics less~
\cite{Martyna1996
}.
1822 Standard velocity Verlet with Nos
{\'e
}-Hoover temperature control has a Trotter expansion
1824 \exp(iL
\dt) &
\approx&
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right)
\exp\left(iL_2
\dt/
2\right)
\nonumber \\
1825 &&
\exp\left(iL_1
\dt\right)
\exp\left(iL_2
\dt/
2\right)
\exp\left(iL_
{\mathrm{NHC
}}\dt/
2\right).
1827 If the Nos
{\'e
}-Hoover chain is sufficiently slow with respect to the motions of the system, we can
1828 write an alternate integrator over $n$ steps for velocity Verlet as
1830 \exp(iL
\dt) &
\approx& (
\exp\left(iL_
{\mathrm{NHC
}}(n
\dt/
2)
\right)
\left[\exp\left(iL_2
\dt/
2\right)
\right.
\nonumber \\
1831 &&
\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).
1833 For pressure control, this becomes
1835 \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 \\
1836 &&
\exp\left(iL_
{\epsilon,
2}(n
\dt/
2)
\right)
\left[\exp\left(iL_2
\dt/
2\right)
\right.
\nonumber \nonumber \\
1837 &&
\exp\left(iL_
{\epsilon,
1}\dt\right)
\exp\left(iL_1
\dt\right)
\nonumber \nonumber \\
1838 &&
\left.
\exp\left(iL_2
\dt/
2\right)
\right]^n
\exp\left(iL_
{\epsilon,
2}(n
\dt/
2)
\right)
\nonumber \nonumber \\
1839 &&
\exp\left(iL_
{\mathrm{NHC
}}(n
\dt/
2)
\right)
\exp\left(iL_
{\mathrm{NHC-baro
}}(n
\dt/
2)
\right),
1841 where the box volume integration occurs every step, but the auxiliary variable
1842 integrations happen every $n$ steps.
1844 % } % Brace matches ifthenelse test for gmxlite
1847 \subsection{The complete update algorithm
}
1850 \addtolength{\fboxsep}{0.5cm
}
1851 \begin{shadowenv
}[12cm
]
1852 {\large \bf THE UPDATE ALGORITHM
}
1853 \rule{\textwidth}{2pt
} \\
1855 Positions $
\ve{r
}$ of all atoms at time $t$ \\
1856 Velocities $
\ve{v
}$ of all atoms at time $t-
\hDt$ \\
1857 Accelerations $
\ve{F
}/m$ on all atoms at time $t$.\\
1858 (Forces are computed disregarding any constraints)\\
1859 Total kinetic energy and virial at $t-
\Dt$\\
1861 {\bf 1.
} Compute the scaling factors $
\lambda$ and $
\mu$\\
1862 according to
\eqnsref{lambda
}{mu
}\\
1864 {\bf 2.
} Update and scale velocities: $
\ve{v
}' =
\lambda (
\ve{v
} +
1865 \ve{a
} \Delta t)$ \\
1867 {\bf 3.
} Compute new unconstrained coordinates: $
\ve{r
}' =
\ve{r
} +
\ve{v
}'
1870 {\bf 4.
} Apply constraint algorithm to coordinates: constrain($
\ve{r
}^
{'
} \rightarrow \ve{r
}'';
1873 {\bf 5.
} Correct velocities for constraints: $
\ve{v
} = (
\ve{r
}'' -
1874 \ve{r
}) /
\Delta t$ \\
1876 {\bf 6.
} Scale coordinates and box: $
\ve{r
} =
\mu \ve{r
}'';
\ve{b
} =
1879 \caption{The MD update algorithm with the leap-frog integrator
}
1880 \label{fig:complete-update
}
1883 The complete algorithm for the update of velocities and coordinates is
1884 given using leap-frog in
\figref{complete-update
}. The SHAKE algorithm of step
1885 4 is explained below.
1887 {\gromacs} has a provision to ``freeze'' (prevent motion of) selected
1888 particles
\index{frozen atoms
}, which must be defined as a ``
\swapindex{freeze
}{group
}.'' This is implemented
1889 using a
{\em freeze factor $
\ve{f
}_g$
}, which is a vector, and differs for each
1890 freeze group (see
\secref{groupconcept
}). This vector contains only
1891 zero (freeze) or one (don't freeze).
1892 When we take this freeze factor and the external acceleration $
\ve{a
}_h$ into
1893 account the update algorithm for the velocities becomes
1895 \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],
1897 where $g$ and $h$ are group indices which differ per atom.
1899 \subsection{Output step
}
1900 The most important output of the MD run is the
{\em
1901 \swapindex{trajectory
}{file
}}, which contains particle coordinates
1902 and (optionally) velocities at regular intervals.
1903 The trajectory file contains frames that could include positions,
1904 velocities and/or forces, as well as information about the dimensions
1905 of the simulation volume, integration step, integration time, etc. The
1906 interpretation of the time varies with the integrator chosen, as
1907 described above. For velocity-Verlet integrators, velocities labeled
1908 at time $t$ are for that time. For other integrators (e.g. leap-frog,
1909 stochastic dynamics), the velocities labeled at time $t$ are for time
1912 Since the trajectory
1913 files are lengthy, one should not save every step! To retain all
1914 information it suffices to write a frame every
15 steps, since at
1915 least
30 steps are made per period of the highest frequency in the
1916 system, and Shannon's
\normindex{sampling
} theorem states that two samples per
1917 period of the highest frequency in a band-limited signal contain all
1918 available information. But that still gives very long files! So, if
1919 the highest frequencies are not of interest,
10 or
20 samples per ps
1920 may suffice. Be aware of the distortion of high-frequency motions by
1921 the
{\em stroboscopic effect
}, called
{\em aliasing
}: higher frequencies
1922 are mirrored with respect to the sampling frequency and appear as
1925 {\gromacs} can also write reduced-precision coordinates for a subset of
1926 the simulation system to a special compressed trajectory file
1927 format. All the other tools can read and write this format. See
1928 \secref{mdpopt
} for details on how to set up your
{\tt .mdp
} file
1929 to have
{\tt mdrun
} use this feature.
1931 % \ifthenelse{\equal{\gmxlite}{1}}{}{
1932 \section{Shell molecular dynamics
}
1933 {\gromacs} can simulate
\normindex{polarizability
} using the
1934 \normindex{shell model
} of Dick and Overhauser~
\cite{Dick58
}. In such models
1935 a shell particle representing the electronic degrees of freedom is
1936 attached to a nucleus by a spring. The potential energy is minimized with
1937 respect to the shell position at every step of the simulation (see below).
1938 Successful applications of shell models in
{\gromacs} have been published
1939 for $N_2$~
\cite{Jordan95
} and water~
\cite{Maaren2001a
}.
1941 \subsection{Optimization of the shell positions
}
1942 The force
\ve{F
}$_S$ on a shell particle $S$ can be decomposed into two
1945 \ve{F
}_S ~=~
\ve{F
}_
{bond
} +
\ve{F
}_
{nb
}
1947 where
\ve{F
}$_
{bond
}$ denotes the component representing the
1948 polarization energy, usually represented by a harmonic potential and
1949 \ve{F
}$_
{nb
}$ is the sum of Coulomb and van der Waals interactions. If we
1950 assume that
\ve{F
}$_
{nb
}$ is almost constant we can analytically derive the
1951 optimal position of the shell, i.e. where
\ve{F
}$_S$ =
0. If we have the
1952 shell S connected to atom A we have
1954 \ve{F
}_
{bond
} ~=~ k_b
\left(
\ve{x
}_S -
\ve{x
}_A
\right).
1956 In an iterative solver, we have positions
\ve{x
}$_S(n)$ where $n$ is
1957 the iteration count. We now have at iteration $n$
1959 \ve{F
}_
{nb
} ~=~
\ve{F
}_S - k_b
\left(
\ve{x
}_S(n) -
\ve{x
}_A
\right)
1961 and the optimal position for the shells $x_S(n+
1)$ thus follows from
1963 \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
1967 \Delta \ve{x
}_S =
\ve{x
}_S(n+
1) -
\ve{x
}_S(n)
1971 \Delta \ve{x
}_S =
\ve{F
}_S/k_b
1973 which then yields the algorithm to compute the next trial in the optimization
1976 \ve{x
}_S(n+
1) ~=~
\ve{x
}_S(n) +
\ve{F
}_S/k_b.
1978 % } % Brace matches ifthenelse test for gmxlite
1980 \section{Constraint algorithms
\index{constraint algorithms
}}
1981 Constraints can be imposed in
{\gromacs} using LINCS (default) or
1982 the traditional SHAKE method.
1984 \subsection{\normindex{SHAKE
}}
1985 \label{subsec:SHAKE
}
1986 The SHAKE~
\cite{Ryckaert77
} algorithm changes a set of unconstrained
1987 coordinates $
\ve{r
}^
{'
}$ to a set of coordinates $
\ve{r
}''$ that
1988 fulfill a list of distance constraints, using a set $
\ve{r
}$
1991 {\rm SHAKE
}(
\ve{r
}^
{'
} \rightarrow \ve{r
}'';\,
\ve{r
})
1993 This action is consistent with solving a set of Lagrange multipliers
1994 in the constrained equations of motion. SHAKE needs a
{\em relative tolerance
};
1995 it will continue until all constraints are satisfied within
1996 that relative tolerance. An error message is
1997 given if SHAKE cannot reset the coordinates because the deviation is
1998 too large, or if a given number of iterations is surpassed.
2000 Assume the equations of motion must fulfill $K$ holonomic constraints,
2003 \sigma_k(
\ve{r
}_1
\ldots \ve{r
}_N) =
0; \;\; k=
1 \ldots K.
2005 For example, $(
\ve{r
}_1 -
\ve{r
}_2)^
2 - b^
2 =
0$.
2006 Then the forces are defined as
2008 -
\frac{\partial}{\partial \ve{r
}_i
} \left( V +
\sum_{k=
1}^K
\lambda_k
2011 where $
\lambda_k$ are Lagrange multipliers which must be solved to
2012 fulfill the constraint equations. The second part of this sum
2013 determines the
{\em constraint forces
} $
\ve{G
}_i$, defined by
2015 \ve{G
}_i = -
\sum_{k=
1}^K
\lambda_k \frac{\partial \sigma_k}{\partial
2018 The displacement due to the constraint forces in the leap-frog or
2019 Verlet algorithm is equal to $(
\ve{G
}_i/m_i)(
\Dt)^
2$. Solving the
2020 Lagrange multipliers (and hence the displacements) requires the
2021 solution of a set of coupled equations of the second degree. These are
2022 solved iteratively by SHAKE.
2023 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2024 \label{subsec:SETTLE
}
2025 For the special case of rigid water molecules, that often make up more
2026 than
80\% of the simulation system we have implemented the
2028 algorithm~
\cite{Miyamoto92
} (
\secref{constraints
}).
2030 For velocity Verlet, an additional round of constraining must be
2031 done, to constrain the velocities of the second velocity half step,
2032 removing any component of the velocity parallel to the bond vector.
2033 This step is called RATTLE, and is covered in more detail in the
2034 original Andersen paper~
\cite{Andersen1983a
}.
2036 % } % Brace matches ifthenelse test for gmxlite
2041 \newcommand{\fs}[1]{\begin{equation
} \label{eqn:
#1}}
2042 \newcommand{\fe}{\end{equation
}}
2043 \newcommand{\p}{\partial}
2044 \newcommand{\Bm}{\ve{B
}}
2045 \newcommand{\M}{\ve{M
}}
2046 \newcommand{\iM}{\M^
{-
1}}
2047 \newcommand{\Tm}{\ve{T
}}
2048 \newcommand{\Sm}{\ve{S
}}
2049 \newcommand{\fo}{\ve{f
}}
2050 \newcommand{\con}{\ve{g
}}
2051 \newcommand{\lenc}{\ve{d
}}
2053 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2054 \subsection{\normindex{LINCS
}}
2055 \label{subsec:lincs
}
2057 \subsubsection{The LINCS algorithm
}
2058 LINCS is an algorithm that resets bonds to their correct lengths
2059 after an unconstrained update~
\cite{Hess97
}.
2060 The method is non-iterative, as it always uses two steps.
2061 Although LINCS is based on matrices, no matrix-matrix multiplications are
2062 needed. The method is more stable and faster than SHAKE,
2063 but it can only be used with bond constraints and
2064 isolated angle constraints, such as the proton angle in OH.
2065 Because of its stability, LINCS is especially useful for Brownian dynamics.
2066 LINCS has two parameters, which are explained in the subsection parameters.
2067 The parallel version of LINCS, P-LINCS, is described
2068 in subsection
\ssecref{plincs
}.
2070 \subsubsection{The LINCS formulas
}
2071 We consider a system of $N$ particles, with positions given by a
2072 $
3N$ vector $
\ve{r
}(t)$.
2073 For molecular dynamics the equations of motion are given by Newton's Law
2075 {\de^
2 \ve{r
} \over \de t^
2} =
\iM \ve{F
},
2077 where $
\ve{F
}$ is the $
3N$ force vector
2078 and $
\M$ is a $
3N
\times 3N$ diagonal matrix,
2079 containing the masses of the particles.
2080 The system is constrained by $K$ time-independent constraint equations
2082 g_i(
\ve{r
}) = |
\ve{r
}_
{i_1
}-
\ve{r
}_
{i_2
} | - d_i =
0 ~~~~~~i=
1,
\ldots,K.
2085 In a numerical integration scheme, LINCS is applied after an
2086 unconstrained update, just like SHAKE. The algorithm works in two
2087 steps (see figure
\figref{lincs
}). In the first step, the projections
2088 of the new bonds on the old bonds are set to zero. In the second step,
2089 a correction is applied for the lengthening of the bonds due to
2090 rotation. The numerics for the first step and the second step are very
2091 similar. A complete derivation of the algorithm can be found in
2092 \cite{Hess97
}. Only a short description of the first step is given
2096 \centerline{\includegraphics[height=
50mm
]{plots/lincs
}}
2097 \caption[The three position updates needed for one time step.
]{The
2098 three position updates needed for one time step. The dashed line is
2099 the old bond of length $d$, the solid lines are the new bonds. $l=d
2100 \cos \theta$ and $p=(
2 d^
2 - l^
2)^
{1 \over 2}$.
}
2104 A new notation is introduced for the gradient matrix of the constraint
2105 equations which appears on the right hand side of this equation:
2107 B_
{hi
} =
{\p g_h
\over \p r_i
}
2109 Notice that $
\Bm$ is a $K
\times 3N$ matrix, it contains the directions
2111 The following equation shows how the new constrained coordinates
2112 $
\ve{r
}_
{n+
1}$ are related to the unconstrained coordinates
2113 $
\ve{r
}_
{n+
1}^
{unc
}$ by
2116 \ve{r
}_
{n+
1}=(
\ve{I
}-
\Tm_n \ve{B
}_n)
\ve{r
}_
{n+
1}^
{unc
} +
\Tm_n \lenc=
2118 \ve{r
}_
{n+
1}^
{unc
} -
2119 \iM \Bm_n (
\Bm_n \iM \Bm_n^T)^
{-
1} (
\Bm_n \ve{r
}_
{n+
1}^
{unc
} -
\lenc)
2122 where $
\Tm =
\iM \Bm^T (
\Bm \iM \Bm^T)^
{-
1}$.
2123 The derivation of this equation from
\eqnsref{c1
}{c2
} can be found
2126 This first step does not set the real bond lengths to the prescribed lengths,
2127 but the projection of the new bonds onto the old directions of the bonds.
2128 To correct for the rotation of bond $i$, the projection of the
2129 bond, $p_i$, on the old direction is set to
2131 p_i=
\sqrt{2 d_i^
2 - l_i^
2},
2133 where $l_i$ is the bond length after the first projection.
2134 The corrected positions are
2136 \ve{r
}_
{n+
1}^*=(
\ve{I
}-
\Tm_n \Bm_n)
\ve{r
}_
{n+
1} +
\Tm_n \ve{p
}.
2138 This correction for rotational effects is actually an iterative process,
2139 but during MD only one iteration is applied.
2140 The relative constraint deviation after this procedure will be less than
2141 0.0001 for every constraint.
2142 In energy minimization, this might not be accurate enough, so the number
2143 of iterations is equal to the order of the expansion (see below).
2145 Half of the CPU time goes to inverting the constraint coupling
2146 matrix $
\Bm_n \iM \Bm_n^T$, which has to be done every time step.
2147 This $K
\times K$ matrix
2148 has $
1/m_
{i_1
} +
1/m_
{i_2
}$ on the diagonal.
2149 The off-diagonal elements are only non-zero when two bonds are connected,
2151 $
\cos \phi /m_c$, where $m_c$ is
2152 the mass of the atom connecting the
2153 two bonds and $
\phi$ is the angle between the bonds.
2155 The matrix $
\Tm$ is inverted through a power expansion.
2156 A $K
\times K$ matrix $
\ve{S
}$ is
2157 introduced which is the inverse square root of
2158 the diagonal of $
\Bm_n \iM \Bm_n^T$.
2159 This matrix is used to convert the diagonal elements
2160 of the coupling matrix to one:
2163 (
\Bm_n \iM \Bm_n^T)^
{-
1}
2164 =
\Sm \Sm^
{-
1} (
\Bm_n \iM \Bm_n^T)^
{-
1} \Sm^
{-
1} \Sm \\
[2mm
]
2165 =
\Sm (
\Sm \Bm_n \iM \Bm_n^T
\Sm)^
{-
1} \Sm =
2166 \Sm (
\ve{I
} -
\ve{A
}_n)^
{-
1} \Sm
2169 The matrix $
\ve{A
}_n$ is symmetric and sparse and has zeros on the diagonal.
2170 Thus a simple trick can be used to calculate the inverse:
2172 (
\ve{I
}-
\ve{A
}_n)^
{-
1}=
2173 \ve{I
} +
\ve{A
}_n +
\ve{A
}_n^
2 +
\ve{A
}_n^
3 +
\ldots
2176 This inversion method is only valid if the absolute values of all the
2177 eigenvalues of $
\ve{A
}_n$ are smaller than one.
2178 In molecules with only bond constraints, the connectivity is so low
2179 that this will always be true, even if ring structures are present.
2180 Problems can arise in angle-constrained molecules.
2181 By constraining angles with additional distance constraints,
2182 multiple small ring structures are introduced.
2183 This gives a high connectivity, leading to large eigenvalues.
2184 Therefore LINCS should NOT be used with coupled angle-constraints.
2186 For molecules with all bonds constrained the eigenvalues of $A$
2187 are around
0.4. This means that with each additional order
2188 in the expansion
\eqnref{m3
} the deviations decrease by a factor
0.4.
2189 But for relatively isolated triangles of constraints the largest
2190 eigenvalue is around
0.7.
2191 Such triangles can occur when removing hydrogen angle vibrations
2192 with an additional angle constraint in alcohol groups
2193 or when constraining water molecules with LINCS, for instance
2194 with flexible constraints.
2195 The constraints in such triangles converge twice as slow as
2196 the other constraints. Therefore, starting with
{\gromacs} 4,
2197 additional terms are added to the expansion for such triangles
2199 (
\ve{I
}-
\ve{A
}_n)^
{-
1} \approx
2200 \ve{I
} +
\ve{A
}_n +
\ldots +
\ve{A
}_n^
{N_i
} +
2201 \left(
\ve{A
}^*_n +
\ldots +
{\ve{A
}_n^*
}^
{N_i
} \right)
\ve{A
}_n^
{N_i
}
2203 where $N_i$ is the normal order of the expansion and
2204 $
\ve{A
}^*$ only contains the elements of $
\ve{A
}$ that couple
2205 constraints within rigid triangles, all other elements are zero.
2206 In this manner, the accuracy of angle constraints comes close
2207 to that of the other constraints, while the series of matrix vector
2208 multiplications required for determining the expansion
2209 only needs to be extended for a few constraint couplings.
2210 This procedure is described in the P-LINCS paper
\cite{Hess2008a
}.
2212 \subsubsection{The LINCS Parameters
}
2213 The accuracy of LINCS depends on the number of matrices used
2214 in the expansion
\eqnref{m3
}. For MD calculations a fourth order
2215 expansion is enough. For Brownian dynamics with
2216 large time steps an eighth order expansion may be necessary.
2217 The order is a parameter in the
{\tt *.mdp
} file.
2218 The implementation of LINCS is done in such a way that the
2219 algorithm will never crash. Even when it is impossible to
2220 to reset the constraints LINCS will generate a conformation
2221 which fulfills the constraints as well as possible.
2222 However, LINCS will generate a warning when in one step a bond
2223 rotates over more than a predefined angle.
2224 This angle is set by the user in the
{\tt *.mdp
} file.
2226 % } % Brace matches ifthenelse test for gmxlite
2229 \section{Simulated Annealing
}
2231 The well known
\swapindex{simulated
}{annealing
}
2232 (SA) protocol is supported in
{\gromacs}, and you can even couple multiple
2233 groups of atoms separately with an arbitrary number of reference temperatures
2234 that change during the simulation. The annealing is implemented by simply
2235 changing the current reference temperature for each group in the temperature
2236 coupling, so the actual relaxation and coupling properties depends on the
2237 type of thermostat you use and how hard you are coupling it. Since we are
2238 changing the reference temperature it is important to remember that the system
2239 will NOT instantaneously reach this value - you need to allow for the inherent
2240 relaxation time in the coupling algorithm too. If you are changing the
2241 annealing reference temperature faster than the temperature relaxation you
2242 will probably end up with a crash when the difference becomes too large.
2244 The annealing protocol is specified as a series of corresponding times and
2245 reference temperatures for each group, and you can also choose whether you only
2246 want a single sequence (after which the temperature will be coupled to the
2247 last reference value), or if the annealing should be periodic and restart at
2248 the first reference point once the sequence is completed. You can mix and
2249 match both types of annealing and non-annealed groups in your simulation.
2251 \newcommand{\vrond}{\stackrel{\circ}{\ve{r
}}}
2252 \newcommand{\rond}{\stackrel{\circ}{r
}}
2253 \newcommand{\ruis}{\ve{r
}^G
}
2255 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2256 \section{Stochastic Dynamics
\swapindexquiet{stochastic
}{dynamics
}}
2258 Stochastic or velocity
\swapindex{Langevin
}{dynamics
} adds a friction
2259 and a noise term to Newton's equations of motion, as
2262 m_i
{\de^
2 \ve{r
}_i
\over \de t^
2} =
2263 - m_i
\gamma_i {\de \ve{r
}_i
\over \de t
} +
\ve{F
}_i(
\ve{r
}) +
\vrond_i,
2265 where $
\gamma_i$ is the friction constant $
[1/
\mbox{ps
}]$ and
2266 $
\vrond_i\!\!(t)$ is a noise process with
2267 $
\langle \rond_i\!\!(t)
\rond_j\!\!(t+s)
\rangle =
2268 2 m_i
\gamma_i k_B T
\delta(s)
\delta_{ij
}$.
2269 When $
1/
\gamma_i$ is large compared to the time scales present in the system,
2270 one could see stochastic dynamics as molecular dynamics with stochastic
2271 temperature-coupling. The advantage compared to MD with Berendsen
2272 temperature-coupling is that in case of SD the generated ensemble is known.
2273 For simulating a system in vacuum there is the additional advantage that there is no
2274 accumulation of errors for the overall translational and rotational
2276 When $
1/
\gamma_i$ is small compared to the time scales present in the system,
2277 the dynamics will be completely different from MD, but the sampling is
2280 In
{\gromacs} there are two algorithms to integrate equation (
\ref{SDeq
}):
2281 a simple and efficient one
2282 and a more complex leap-frog algorithm~
\cite{Gunsteren88
}.
2283 The accuracy of both integrators is equivalent to the normal MD leap-frog and
2284 velocity-Verlet integrator, except with constraints where the simple
2285 SD integrator is significantly less accurate. There is a proper way
2286 of applying constraints with the simple integrator, but that requires
2287 a second constraining step~
\cite{Goga2012
}, which diminishes the gain.
2288 The simple integrator is:
2291 \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 \\
2292 \ve{r
}(t+
\Dt) &~=~&
\ve{r
}(t)+
\Dt \,
\ve{v
}(t+
\hDt) \\
2293 \alpha &~=~&
\left(
1 -
\frac{\gamma \Dt}{m
} \right)
2295 where $
\ruis_i$ is Gaussian distributed noise with $
\mu =
0$, $
\sigma =
1$.
2296 With constraints you should only consider using the simple integrator when $
\gamma \Dt/m
\ll 0.01$.
2298 In the complex algorithm four Gaussian random numbers are required
2299 per integration step per degree of freedom, and with constraints the
2300 coordinates need to be constrained twice per integration step.
2301 Depending on the computational cost of the force calculation,
2302 this can take a significant part of the simulation time.
2303 Exact continuation of a stochastic dynamics simulation is not possible,
2304 because the state of the random number generator is not stored.
2305 When using SD as a thermostat, an appropriate value for $
\gamma$ is
0.5 ps$^
{-
1}$,
2306 since this results in a friction that is lower than the internal friction
2307 of water, while it is high enough to remove excess heat
2308 (unless plain cut-off or reaction-field electrostatics is used).
2309 With this value of $
\gamma$ the efficient algorithm will usually be accurate
2312 \section{Brownian Dynamics
\swapindexquiet{Brownian
}{dynamics
}}
2314 In the limit of high friction, stochastic dynamics reduces to
2315 Brownian dynamics, also called position Langevin dynamics.
2316 This applies to over-damped systems,
2317 {\ie} systems in which the inertia effects are negligible.
2320 {\de \ve{r
}_i
\over \de t
} =
\frac{1}{\gamma_i} \ve{F
}_i(
\ve{r
}) +
\vrond_i
2322 where $
\gamma_i$ is the friction coefficient $
[\mbox{amu/ps
}]$ and
2323 $
\vrond_i\!\!(t)$ is a noise process with
2324 $
\langle \rond_i\!\!(t)
\rond_j\!\!(t+s)
\rangle =
2325 2 \delta(s)
\delta_{ij
} k_B T /
\gamma_i$.
2326 In
{\gromacs} the equations are integrated with a simple, explicit scheme
2328 \ve{r
}_i(t+
\Delta t) =
\ve{r
}_i(t) +
2329 {\Delta t
\over \gamma_i} \ve{F
}_i(
\ve{r
}(t))
2330 +
\sqrt{2 k_B T
{\Delta t
\over \gamma_i}}\,
\ruis_i,
2332 where $
\ruis_i$ is Gaussian distributed noise with $
\mu =
0$, $
\sigma =
1$.
2333 The friction coefficients $
\gamma_i$ can be chosen the same for all
2334 particles or as $
\gamma_i = m_i\,
\gamma_i$, where the friction constants
2335 $
\gamma_i$ can be different for different groups of atoms.
2336 Because the system is assumed to be over-damped, large timesteps
2337 can be used. LINCS should be used for the constraints since SHAKE
2338 will not converge for large atomic displacements.
2339 BD is an option of the
{\tt mdrun
} program.
2340 % } % Brace matches ifthenelse test for gmxlite
2342 \section{Energy Minimization
}
2343 \label{sec:EM
}\index{energy minimization
}%
2344 Energy minimization in
{\gromacs} can be done using steepest descent,
2345 conjugate gradients, or l-bfgs (limited-memory
2346 Broyden-Fletcher-Goldfarb-Shanno quasi-Newtonian minimizer...we
2347 prefer the abbreviation). EM is just an option of the
{\tt mdrun
}
2350 \subsection{Steepest Descent
\index{steepest descent
}}
2351 Although steepest descent is certainly not the most efficient
2352 algorithm for searching, it is robust and easy to implement.
2354 We define the vector $
\ve{r
}$ as the vector of all $
3N$ coordinates.
2355 Initially a maximum displacement $h_0$ (
{\eg} 0.01 nm) must be given.
2357 First the forces $
\ve{F
}$ and potential energy are calculated.
2358 New positions are calculated by
2360 \ve{r
}_
{n+
1} =
\ve{r
}_n +
\frac{\ve{F
}_n
}{\max (|
\ve{F
}_n|)
} h_n,
2362 where $h_n$ is the maximum displacement and $
\ve{F
}_n$ is the force,
2363 or the negative gradient of the potential $V$. The notation $
\max
2364 (|
\ve{F
}_n|)$ means the largest of the absolute values of the force
2365 components. The forces and energy are again computed for the new positions \\
2366 If ($V_
{n+
1} < V_n$) the new positions are accepted and $h_
{n+
1} =
1.2
2368 If ($V_
{n+
1} \geq V_n$) the new positions are rejected and $h_n =
0.2 h_n$.
2370 The algorithm stops when either a user-specified number of force
2371 evaluations has been performed (
{\eg} 100), or when the maximum of the absolute
2372 values of the force (gradient) components is smaller than a specified
2374 Since force truncation produces some noise in the
2375 energy evaluation, the stopping criterion should not be made too tight
2376 to avoid endless iterations. A reasonable value for $
\epsilon$ can be
2377 estimated from the root mean square force $f$ a harmonic oscillator would exhibit at a
2378 temperature $T$. This value is
2380 f =
2 \pi \nu \sqrt{ 2mkT
},
2382 where $
\nu$ is the oscillator frequency, $m$ the (reduced) mass, and
2383 $k$ Boltzmann's constant. For a weak oscillator with a wave number of
2384 100 cm$^
{-
1}$ and a mass of
10 atomic units, at a temperature of
1 K,
2385 $f=
7.7$ kJ~mol$^
{-
1}$~nm$^
{-
1}$. A value for $
\epsilon$ between
1 and
2388 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2389 \subsection{Conjugate Gradient
\index{conjugate gradient
}}
2390 Conjugate gradient is slower than steepest descent in the early stages
2391 of the minimization, but becomes more efficient closer to the energy
2392 minimum. The parameters and stop criterion are the same as for
2393 steepest descent. In
{\gromacs} conjugate gradient can not be used
2394 with constraints, including the SETTLE algorithm for
2395 water~
\cite{Miyamoto92
}, as this has not been implemented. If water is
2396 present it must be of a flexible model, which can be specified in the
2397 {\tt *.mdp
} file by
{\tt define = -DFLEXIBLE
}.
2399 This is not really a restriction, since the accuracy of conjugate
2400 gradient is only required for minimization prior to a normal-mode
2401 analysis, which cannot be performed with constraints. For most other
2402 purposes steepest descent is efficient enough.
2403 % } % Brace matches ifthenelse test for gmxlite
2405 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2406 \subsection{\normindex{L-BFGS
}}
2407 The original BFGS algorithm works by successively creating better
2408 approximations of the inverse Hessian matrix, and moving the system to
2409 the currently estimated minimum. The memory requirements for this are
2410 proportional to the square of the number of particles, so it is not
2411 practical for large systems like biomolecules. Instead, we use the
2412 L-BFGS algorithm of Nocedal~
\cite{Byrd95a,Zhu97a
}, which approximates
2413 the inverse Hessian by a fixed number of corrections from previous
2414 steps. This sliding-window technique is almost as efficient as the
2415 original method, but the memory requirements are much lower -
2416 proportional to the number of particles multiplied with the correction
2417 steps. In practice we have found it to converge faster than conjugate
2418 gradients, but due to the correction steps it is not yet parallelized.
2419 It is also noteworthy that switched or shifted interactions usually
2420 improve the convergence, since sharp cut-offs mean the potential
2421 function at the current coordinates is slightly different from the
2422 previous steps used to build the inverse Hessian approximation.
2423 % } % Brace matches ifthenelse test for gmxlite
2425 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2426 \section{Normal-Mode Analysis
\index{normal-mode analysis
}\index{NMA
}}
2427 Normal-mode analysis~
\cite{Levitt83,Go83,BBrooks83b
}
2428 can be performed using
{\gromacs}, by diagonalization of the mass-weighted
2429 \normindex{Hessian
} $H$:
2431 R^T M^
{-
1/
2} H M^
{-
1/
2} R &=&
\mbox{diag
}(
\lambda_1,
\ldots,
\lambda_{3N
})
2433 \lambda_i &=& (
2 \pi \omega_i)^
2
2435 where $M$ contains the atomic masses, $R$ is a matrix that contains
2436 the eigenvectors as columns, $
\lambda_i$ are the eigenvalues
2437 and $
\omega_i$ are the corresponding frequencies.
2439 First the Hessian matrix, which is a $
3N
\times 3N$ matrix where $N$
2440 is the number of atoms, needs to be calculated:
2442 H_
{ij
} &=&
\frac{\partial^
2 V
}{\partial x_i
\partial x_j
}
2444 where $x_i$ and $x_j$ denote the atomic x, y or z coordinates.
2445 In practice, this equation is not used, but the Hessian is
2446 calculated numerically from the force as:
2449 \frac{f_i(
{\bf x
}+h
{\bf e
}_j) - f_i(
{\bf x
}-h
{\bf e
}_j)
}{2h
}
2451 f_i &=& -
\frac{\partial V
}{\partial x_i
}
2453 where $
{\bf e
}_j$ is the unit vector in direction $j$.
2454 It should be noted that
2455 for a usual normal-mode calculation, it is necessary to completely minimize
2456 the energy prior to computation of the Hessian.
2457 The tolerance required depends on the type of system,
2458 but a rough indication is
0.001 kJ mol$^
{-
1}$.
2459 Minimization should be done with conjugate gradients or L-BFGS in double precision.
2461 A number of
{\gromacs} programs are involved in these
2462 calculations. First, the energy should be minimized using
{\tt mdrun
}.
2463 Then,
{\tt mdrun
} computes the Hessian.
{\bf Note
} that for generating
2464 the run input file, one should use the minimized conformation from
2465 the full precision trajectory file, as the structure file is not
2467 {\tt \normindex{g_nmeig
}} does the diagonalization and
2468 the sorting of the normal modes according to their frequencies.
2469 Both
{\tt mdrun
} and
{\tt g_nmeig
} should be run in double precision.
2470 The normal modes can be analyzed with the program
{\tt g_anaeig
}.
2471 Ensembles of structures at any temperature and for any subset of
2472 normal modes can be generated with
{\tt \normindex{g_nmens
}}.
2473 An overview of normal-mode analysis and the related principal component
2474 analysis (see
\secref{covanal
}) can be found in~
\cite{Hayward95b
}.
2475 % } % Brace matches ifthenelse test for gmxlite
2477 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2479 \section{Free energy calculations
\index{free energy calculations
}}
2481 \subsection{Slow-growth methods
\index{slow-growth methods
}}
2482 Free energy calculations can be performed
2483 in
{\gromacs} using a number of methods, including ``slow-growth.'' An example problem
2484 might be calculating the difference in free energy of binding of an inhibitor
{\bf I
}
2485 to an enzyme
{\bf E
} and to a mutated enzyme
{\bf E$^
{\prime}$
}. It
2486 is not feasible with computer simulations to perform a docking
2487 calculation for such a large complex, or even releasing the inhibitor from
2488 the enzyme in a reasonable amount of computer time with reasonable accuracy.
2489 However, if we consider the free energy cycle in~
\figref{free
}A
2492 \Delta G_1 -
\Delta G_2 =
\Delta G_3 -
\Delta G_4
2495 If we are interested in the left-hand term we can equally well compute
2496 the right-hand term.
2498 \centerline{\includegraphics[width=
6cm,angle=
270]{plots/free1
}\hspace{2cm
}\includegraphics[width=
6cm,angle=
270]{plots/free2
}}
2499 \caption[Free energy cycles.
]{Free energy cycles.
{\bf A:
} to
2500 calculate $
\Delta G_
{12}$, the free energy difference between the
2501 binding of inhibitor
{\bf I
} to enzymes
{\bf E
} respectively
{\bf
2502 E$^
{\prime}$
}.
{\bf B:
} to calculate $
\Delta G_
{12}$, the free energy
2503 difference for binding of inhibitors
{\bf I
} respectively
{\bf I$^
{\prime}$
} to
2508 If we want to compute the difference in free energy of binding of two
2509 inhibitors
{\bf I
} and
{\bf I$^
{\prime}$
} to an enzyme
{\bf E
} (
\figref{free
}B)
2510 we can again use
\eqnref{ddg
} to compute the desired property.
2512 \newcommand{\sA}{^
{\mathrm{A
}}}
2513 \newcommand{\sB}{^
{\mathrm{B
}}}
2514 Free energy differences between two molecular species can
2515 be calculated in
{\gromacs} using the ``slow-growth'' method.
2516 Such free energy differences between different molecular species are
2517 physically meaningless, but they can be used to obtain meaningful
2518 quantities employing a thermodynamic cycle.
2519 The method requires a simulation during which the Hamiltonian of the
2520 system changes slowly from that describing one system (A) to that
2521 describing the other system (B). The change must be so slow that the
2522 system remains in equilibrium during the process; if that requirement
2523 is fulfilled, the change is reversible and a slow-growth simulation from B to A
2524 will yield the same results (but with a different sign) as a slow-growth
2525 simulation from A to B. This is a useful check, but the user should be
2526 aware of the danger that equality of forward and backward growth results does
2527 not guarantee correctness of the results.
2529 The required modification of the Hamiltonian $H$ is realized by making
2530 $H$ a function of a
\textit{coupling parameter
} $
\lambda:
2531 H=H(p,q;
\lambda)$ in such a way that $
\lambda=
0$ describes system A
2532 and $
\lambda=
1$ describes system B:
2534 H(p,q;
0)=H
\sA (p,q);~~~~ H(p,q;
1)=H
\sB (p,q).
2536 In
{\gromacs}, the functional form of the $
\lambda$-dependence is
2537 different for the various force-field contributions and is described
2538 in section
\secref{feia
}.
2540 The Helmholtz free energy $A$ is related to the
2541 partition function $Q$ of an $N,V,T$ ensemble, which is assumed to be
2542 the equilibrium ensemble generated by a MD simulation at constant
2543 volume and temperature. The generally more useful Gibbs free energy
2544 $G$ is related to the partition function $
\Delta$ of an $N,p,T$
2545 ensemble, which is assumed to be the equilibrium ensemble generated by
2546 a MD simulation at constant pressure and temperature:
2548 A(
\lambda) &=& -k_BT
\ln Q \\
2549 Q &=& c
\int\!\!
\int \exp[-
\beta H(p,q;
\lambda)
]\,dp\,dq \\
2550 G(
\lambda) &=& -k_BT
\ln \Delta \\
2551 \Delta &=& c
\int\!\!
\int\!\!
\int \exp[-
\beta H(p,q;
\lambda) -
\beta
2555 where $
\beta =
1/(k_BT)$ and $c = (N! h^
{3N
})^
{-
1}$.
2556 These integrals over phase space cannot be evaluated from a
2557 simulation, but it is possible to evaluate the derivative with
2558 respect to $
\lambda$ as an ensemble average:
2560 \frac{dA
}{d
\lambda} =
\frac{\int\!\!
\int (
\partial H/
\partial
2561 \lambda)
\exp[-
\beta H(p,q;
\lambda)
]\,dp\,dq
}{\int\!\!
\int \exp[-
\beta
2562 H(p,q;
\lambda)
]\,dp\,dq
} =
2563 \left\langle \frac{\partial H
}{\partial \lambda} \right\rangle_{NVT;
\lambda},
2565 with a similar relation for $dG/d
\lambda$ in the $N,p,T$
2566 ensemble. The difference in free energy between A and B can be found
2567 by integrating the derivative over $
\lambda$:
2569 A
\sB(V,T)-A
\sA(V,T) &=&
\int_0^
1 \left\langle \frac{\partial
2570 H
}{\partial \lambda} \right\rangle_{NVT;
\lambda} \,d
\lambda
2572 G
\sB(p,T)-G
\sA(p,T) &=&
\int_0^
1 \left\langle \frac{\partial
2573 H
}{\partial \lambda} \right\rangle_{NpT;
\lambda} \,d
\lambda.
2576 If one wishes to evaluate $G
\sB(p,T)-G
\sA(p,T)$,
2577 the natural choice is a constant-pressure simulation. However, this
2578 quantity can also be obtained from a slow-growth simulation at
2579 constant volume, starting with system A at pressure $p$ and volume $V$
2580 and ending with system B at pressure $p_B$, by applying the following
2581 small (but, in principle, exact) correction:
2584 A
\sB(V)-A
\sA(V) -
\int_p^
{p
\sB}[V
\sB(p')-V
]\,dp'
2586 Here we omitted the constant $T$ from the notation. This correction is
2587 roughly equal to $-
\frac{1}{2} (p
\sB-p)
\Delta V=(
\Delta V)^
2/(
2
2588 \kappa V)$, where $
\Delta V$ is the volume change at $p$ and $
\kappa$
2589 is the isothermal compressibility. This is usually
2590 small; for example, the growth of a water molecule from nothing
2591 in a bath of
1000 water molecules at constant volume would produce an
2592 additional pressure of as much as
22 bar, but a correction to the
2593 Helmholtz free energy of just -
1 kJ mol$^
{-
1}$.
%-20 J/mol.
2595 In Cartesian coordinates, the kinetic energy term in the Hamiltonian
2596 depends only on the momenta, and can be separately integrated and, in
2597 fact, removed from the equations. When masses do not change, there is
2598 no contribution from the kinetic energy at all; otherwise the
2599 integrated contribution to the free energy is $-
\frac{3}{2} k_BT
\ln
2600 (m
\sB/m
\sA)$.
{\bf Note
} that this is only true in the absence of constraints.
2602 \subsection{Thermodynamic integration
\index{thermodynamic integration
}\index{BAR
}\index{Bennett's acceptance ratio
}}
2603 {\gromacs} offers the possibility to integrate eq.~
\ref{eq:delA
} or
2604 eq.
\ref{eq:delG
} in one simulation over the full range from A to
2605 B. However, if the change is large and insufficient sampling can be
2606 expected, the user may prefer to determine the value of $
\langle
2607 dG/d
\lambda \rangle$ accurately at a number of well-chosen
2608 intermediate values of $
\lambda$. This can easily be done by setting
2609 the stepsize
{\tt delta_lambda
} to zero. Each simulation can be
2610 equilibrated first, and a proper error estimate can be made for each
2611 value of $dG/d
\lambda$ from the fluctuation of $
\partial H/
\partial
2612 \lambda$. The total free energy change is then determined afterward
2613 by an appropriate numerical integration procedure.
2615 {\gromacs} now also supports the use of Bennett's Acceptance Ratio~
\cite{Bennett1976
}
2616 for calculating values of $
\Delta$G for transformations from state A to state B using
2617 the program
{\tt \normindex{g_bar
}}. The same data can also be used to calculate free
2618 energies using MBAR~
\cite{Shirts2008
}, though the analysis currently requires external tools from
2619 the external
{\tt pymbar
} package, at https://SimTK.org/home/pymbar.
2621 The $
\lambda$-dependence for the force-field contributions is
2622 described in detail in section
\secref{feia
}.
2623 % } % Brace matches ifthenelse test for gmxlite
2625 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2626 \section{Replica exchange
\index{replica exchange
}}
2627 Replica exchange molecular dynamics (
\normindex{REMD
})
2628 is a method that can be used to speed up
2629 the sampling of any type of simulation, especially if
2630 conformations are separated by relatively high energy barriers.
2631 It involves simulating multiple replicas of the same system
2632 at different temperatures and randomly exchanging the complete state
2633 of two replicas at regular intervals with the probability:
2635 P(
1 \leftrightarrow 2)=
\min\left(
1,
\exp\left[
2636 \left(
\frac{1}{k_B T_1
} -
\frac{1}{k_B T_2
}\right)(U_1 - U_2)
2639 where $T_1$ and $T_2$ are the reference temperatures and $U_1$ and $U_2$
2640 are the instantaneous potential energies of replicas
1 and
2 respectively.
2641 After exchange the velocities are scaled by $(T_1/T_2)^
{\pm0.5
}$
2642 and a neighbor search is performed the next step.
2643 This combines the fast sampling and frequent barrier-crossing
2644 of the highest temperature with correct Boltzmann sampling at
2645 all the different temperatures~
\cite{Hukushima96a,Sugita99
}.
2646 We only attempt exchanges for neighboring temperatures as the probability
2647 decreases very rapidly with the temperature difference.
2648 One should not attempt exchanges for all possible pairs in one step.
2649 If, for instance, replicas
1 and
2 would exchange, the chance of
2650 exchange for replicas
2 and
3 not only depends on the energies of
2651 replicas
2 and
3, but also on the energy of replica
1.
2652 In
{\gromacs} this is solved by attempting exchange for all ``odd''
2653 pairs on ``odd'' attempts and for all ``even'' pairs on ``even'' attempts.
2654 If we have four replicas:
0,
1,
2 and
3, ordered in temperature
2655 and we attempt exchange every
1000 steps, pairs
0-
1 and
2-
3
2656 will be tried at steps
1000,
3000 etc. and pair
1-
2 at steps
2000,
4000 etc.
2658 How should one choose the temperatures?
2659 The energy difference can be written as:
2661 U_1 - U_2 = N_
{df
} \frac{c
}{2} k_B (T_1 - T_2)
2663 where $N_
{df
}$ is the total number of degrees of freedom of one replica
2664 and $c$ is
1 for harmonic potentials and around
2 for protein/water systems.
2665 If $T_2 = (
1+
\epsilon) T_1$ the probability becomes:
2667 P(
1 \leftrightarrow 2)
2668 =
\exp\left( -
\frac{\epsilon^
2 c\,N_
{df
}}{2 (
1+
\epsilon)
} \right)
2669 \approx \exp\left(-
\epsilon^
2 \frac{c
}{2} N_
{df
} \right)
2671 Thus for a probability of $e^
{-
2}\approx 0.135$
2672 one obtains $
\epsilon \approx 2/
\sqrt{c\,N_
{df
}}$.
2673 With all bonds constrained one has $N_
{df
} \approx 2\, N_
{atoms
}$
2674 and thus for $c$ =
2 one should choose $
\epsilon$ as $
1/
\sqrt{N_
{atoms
}}$.
2675 However there is one problem when using pressure coupling. The density at
2676 higher temperatures will decrease, leading to higher energy~
\cite{Seibert2005a
},
2677 which should be taken into account. The
{\gromacs} website features a
2678 so-called ``REMD calculator,'' that lets you type in the temperature range and
2679 the number of atoms, and based on that proposes a set of temperatures.
2681 An extension to the REMD for the isobaric-isothermal ensemble was
2682 proposed by Okabe
{\em et al.
}~
\cite{Okabe2001a
}. In this work the
2683 exchange probability is modified to:
2685 P(
1 \leftrightarrow 2)=
\min\left(
1,
\exp\left[
2686 \left(
\frac{1}{k_B T_1
} -
\frac{1}{k_B T_2
}\right)(U_1 - U_2) +
2687 \left(
\frac{P_1
}{k_B T_1
} -
\frac{P_2
}{k_B T_2
}\right)
\left(V_1-V_2
\right)
2690 where $P_1$ and $P_2$ are the respective reference pressures and $V_1$ and
2691 $V_2$ are the respective instantaneous volumes in the simulations.
2692 In most cases the differences in volume are so small that the second
2693 term is negligible. It only plays a role when the difference between
2694 $P_1$ and $P_2$ is large or in phase transitions.
2696 Hamiltonian replica exchange is also supported in
{\gromacs}. In
2697 Hamiltonian replica exchange, each replica has a different
2698 Hamiltonian, defined by the free energy pathway specified for the simulation. The
2699 exchange probability to maintain the correct ensemble probabilities is:
2700 \beq P(
1 \leftrightarrow 2)=
\min\left(
1,
\exp\left[
2701 \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)))
2705 The separate Hamiltonians are defined by the free energy functionality
2706 of
{\gromacs}, with swaps made between the different values of
2707 $
\lambda$ defined in the mdp file.
2709 Hamiltonian and temperature replica exchange can also be performed
2710 simultaneously, using the acceptance criteria:
2712 P(
1 \leftrightarrow 2)=
\min\left(
1,
\exp\left[
2713 \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
})
2717 Gibbs sampling replica exchange has also been implemented in
2718 {\gromacs}~
\cite{Chodera2011
}. In Gibbs sampling replica exchange, all
2719 possible pairs are tested for exchange, allowing swaps between
2720 replicas that are not neighbors.
2722 Gibbs sampling replica exchange requires no additional potential
2723 energy calculations. However there is an additional communication
2724 cost in Gibbs sampling replica exchange, as for some permutations,
2725 more than one round of swaps must take place. In some cases, this
2726 extra communication cost might affect the efficiency.
2728 All replica exchange variants are options of the
{\tt mdrun
}
2729 program. It will only work when MPI is installed, due to the inherent
2730 parallelism in the algorithm. For efficiency each replica can run on a
2731 separate node. See the manual page of
{\tt mdrun
} on how to use these
2734 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2736 \section{Essential Dynamics sampling
\index{essential dynamics
}\index{principal component analysis
}\seeindexquiet{PCA
}{covariance analysis
}}
2737 The results from Essential Dynamics (see
\secref{covanal
})
2738 of a protein can be used to guide MD simulations. The idea is that
2739 from an initial MD simulation (or from other sources) a definition of
2740 the collective fluctuations with largest amplitude is obtained. The
2741 position along one or more of these collective modes can be
2742 constrained in a (second) MD simulation in a number of ways for
2743 several purposes. For example, the position along a certain mode may
2744 be kept fixed to monitor the average force (free-energy gradient) on
2745 that coordinate in that position. Another application is to enhance
2746 sampling efficiency with respect to usual MD
2747 \cite{Degroot96a,Degroot96b
}. In this case, the system is encouraged
2748 to sample its available configuration space more systematically than
2749 in a diffusion-like path that proteins usually take.
2751 Another possibility to enhance sampling is
\normindex{flooding
}.
2752 Here a flooding potential is added to certain
2753 (collective) degrees of freedom to expel the system out
2754 of a region of phase space
\cite{Lange2006a
}.
2756 The procedure for essential dynamics sampling or flooding is as follows.
2757 First, the eigenvectors and eigenvalues need to be determined
2758 using covariance analysis (
{\tt g_covar
})
2759 or normal-mode analysis (
{\tt g_nmeig
}).
2760 Then, this information is fed into
{\tt make_edi
},
2761 which has many options for selecting vectors and setting parameters,
2762 see
{\tt gmx make_edi -h
}.
2763 The generated
{\tt edi
} input file is then passed to
{\tt mdrun
}.
2765 % } % Brace matches ifthenelse test for gmxlite
2767 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2768 \section{\normindex{Expanded Ensemble
}}
2770 In an expanded ensemble simulation~
\cite{Lyubartsev1992
}, both the coordinates and the
2771 thermodynamic ensemble are treated as configuration variables that can
2772 be sampled over. The probability of any given state can be written as:
2774 P(
\vec{x
},k)
\propto \exp\left(-
\beta_k U_k + g_k
\right),
2776 where $
\beta_k =
\frac{1}{k_B T_k
}$ is the $
\beta$ corresponding to the $k$th
2777 thermodynamic state, and $g_k$ is a user-specified weight factor corresponding
2778 to the $k$th state. This space is therefore a
{\em mixed
},
{\em generalized
}, or
{\em
2779 expanded
} ensemble which samples from multiple thermodynamic
2780 ensembles simultaneously. $g_k$ is chosen to give a specific weighting
2781 of each subensemble in the expanded ensemble, and can either be fixed,
2782 or determined by an iterative procedure. The set of $g_k$ is
2783 frequently chosen to give each thermodynamic ensemble equal
2784 probability, in which case $g_k$ is equal to the free energy in
2785 non-dimensional units, but they can be set to arbitrary values as
2786 desired. Several different algorithms can be used to equilibrate
2787 these weights, described in the mdp option listings.
2788 % } % Brace matches ifthenelse test for gmxlite
2790 In
{\gromacs}, this space is sampled by alternating sampling in the $k$
2791 and $
\vec{x
}$ directions. Sampling in the $
\vec{x
}$ direction is done
2792 by standard molecular dynamics sampling; sampling between the
2793 different thermodynamics states is done by Monte Carlo, with several
2794 different Monte Carlo moves supported. The $k$ states can be defined
2795 by different temperatures, or choices of the free energy $
\lambda$
2796 variable, or both. Expanded ensemble simulations thus represent a
2797 serialization of the replica exchange formalism, allowing a single
2798 simulation to explore many thermodynamic states.
2802 \section{Parallelization
\index{parallelization
}}
2803 The CPU time required for a simulation can be reduced by running the simulation
2804 in parallel over more than one processor or processor core.
2805 Ideally one would want to have linear scaling: running on $N$ processors/cores
2806 makes the simulation $N$ times faster. In practice this can only be
2807 achieved for a small number of processors. The scaling will depend
2808 a lot on the algorithms used. Also, different algorithms can have different
2809 restrictions on the interaction ranges between atoms.
2811 \section{Domain decomposition
\index{domain decomposition
}}
2812 Since most interactions in molecular simulations are local,
2813 domain decomposition is a natural way to decompose the system.
2814 In domain decomposition, a spatial domain is assigned to each processor,
2815 which will then integrate the equations of motion for the particles
2816 that currently reside in its local domain. With domain decomposition,
2817 there are two choices that have to be made: the division of the unit cell
2818 into domains and the assignment of the forces to processors.
2819 Most molecular simulation packages use the half-shell method for assigning
2820 the forces. But there are two methods that always require less communication:
2821 the eighth shell~
\cite{Liem1991
} and the midpoint~
\cite{Shaw2006
} method.
2822 {\gromacs} currently uses the eighth shell method, but for certain systems
2823 or hardware architectures it might be advantageous to use the midpoint
2824 method. Therefore, we might implement the midpoint method in the future.
2825 Most of the details of the domain decomposition can be found
2826 in the
{\gromacs} 4 paper~
\cite{Hess2008b
}.
2828 \subsection{Coordinate and force communication
}
2829 In the most general case of a triclinic unit cell,
2830 the space in divided with a
1-,
2-, or
3-D grid in parallelepipeds
2831 that we call domain decomposition cells.
2832 Each cell is assigned to a processor.
2833 The system is partitioned over the processors at the beginning
2834 of each MD step in which neighbor searching is performed.
2835 Since the neighbor searching is based on charge groups, charge groups
2836 are also the units for the domain decomposition.
2837 Charge groups are assigned to the cell where their center of geometry resides.
2838 Before the forces can be calculated, the coordinates from some
2839 neighboring cells need to be communicated,
2840 and after the forces are calculated, the forces need to be communicated
2841 in the other direction.
2842 The communication and force assignment is based on zones that
2843 can cover one or multiple cells.
2844 An example of a zone setup is shown in
\figref{ddcells
}.
2847 \centerline{\includegraphics[width=
6cm
]{plots/dd-cells
}}
2849 A non-staggered domain decomposition grid of
3$
\times$
2$
\times$
2 cells.
2850 Coordinates in zones
1 to
7 are communicated to the corner cell
2851 that has its home particles in zone
0.
2852 $r_c$ is the cut-off radius.
2857 The coordinates are communicated by moving data along the ``negative''
2858 direction in $x$, $y$ or $z$ to the next neighbor. This can be done in one
2859 or multiple pulses. In
\figref{ddcells
} two pulses in $x$ are required,
2860 then one in $y$ and then one in $z$. The forces are communicated by
2861 reversing this procedure. See the
{\gromacs} 4 paper~
\cite{Hess2008b
}
2862 for details on determining which non-bonded and bonded forces
2863 should be calculated on which node.
2865 \subsection{Dynamic load balancing
\swapindexquiet{dynamic
}{load balancing
}}
2866 When different processors have a different computational load
2867 (load imbalance), all processors will have to wait for the one
2868 that takes the most time. One would like to avoid such a situation.
2869 Load imbalance can occur due to three reasons:
2871 \item inhomogeneous particle distribution
2872 \item inhomogeneous interaction cost distribution (charged/uncharged,
2873 water/non-water due to
{\gromacs} water innerloops)
2874 \item statistical fluctuation (only with small particle numbers)
2876 So we need a dynamic load balancing algorithm
2877 where the volume of each domain decomposition cell
2878 can be adjusted
{\em independently
}.
2879 To achieve this, the
2- or
3-D domain decomposition grids need to be
2880 staggered.
\figref{ddtric
} shows the most general case in
2-D.
2881 Due to the staggering, one might require two distance checks
2882 for deciding if a charge group needs to be communicated:
2883 a non-bonded distance and a bonded distance check.
2886 \centerline{\includegraphics[width=
7cm
]{plots/dd-tric
}}
2888 The zones to communicate to the processor of zone
0,
2889 see the text for details. $r_c$ and $r_b$ are the non-bonded
2890 and bonded cut-off radii respectively, $d$ is an example
2891 of a distance between following, staggered boundaries of cells.
2896 By default,
{\tt mdrun
} automatically turns on the dynamic load
2897 balancing during a simulation when the total performance loss
2898 due to the force calculation imbalance is
5\% or more.
2899 {\bf Note
} that the reported force load imbalance numbers might be higher,
2900 since the force calculation is only part of work that needs to be done
2901 during an integration step.
2902 The load imbalance is reported in the log file at log output steps
2903 and when the
{\tt -v
} option is used also on screen.
2904 The average load imbalance and the total performance loss
2905 due to load imbalance are reported at the end of the log file.
2907 There is one important parameter for the dynamic load balancing,
2908 which is the minimum allowed scaling. By default, each dimension
2909 of the domain decomposition cell can scale down by at least
2910 a factor of
0.8. For
3-D domain decomposition this allows cells
2911 to change their volume by about a factor of
0.5, which should allow
2912 for compensation of a load imbalance of
100\%.
2913 The required scaling can be changed with the
{\tt -dds
} option of
{\tt mdrun
}.
2915 \subsection{Constraints in parallel
\index{constraints
}}
2916 \label{subsec:plincs
}
2917 Since with domain decomposition parts of molecules can reside
2918 on different processors, bond constraints can cross cell boundaries.
2919 Therefore a parallel constraint algorithm is required.
2920 {\gromacs} uses the
\normindex{P-LINCS
} algorithm~
\cite{Hess2008a
},
2921 which is the parallel version of the
\normindex{LINCS
} algorithm~
\cite{Hess97
}
2922 % \ifthenelse{\equal{\gmxlite}{1}}
2924 {(see
\ssecref{lincs
}).
}
2925 The P-LINCS procedure is illustrated in
\figref{plincs
}.
2926 When molecules cross the cell boundaries, atoms in such molecules
2927 up to (
{\tt lincs_order +
1}) bonds away are communicated over the cell boundaries.
2928 Then, the normal LINCS algorithm can be applied to the local bonds
2929 plus the communicated ones. After this procedure, the local bonds
2930 are correctly constrained, even though the extra communicated ones are not.
2931 One coordinate communication step is required for the initial LINCS step
2932 and one for each iteration. Forces do not need to be communicated.
2935 \centerline{\includegraphics[width=
6cm
]{plots/par-lincs2
}}
2937 Example of the parallel setup of P-LINCS with one molecule
2938 split over three domain decomposition cells, using a matrix
2939 expansion order of
3.
2940 The top part shows which atom coordinates need to be communicated
2941 to which cells. The bottom parts show the local constraints (solid)
2942 and the non-local constraints (dashed) for each of the three cells.
2947 \subsection{Interaction ranges
}
2948 Domain decomposition takes advantage of the locality of interactions.
2949 This means that there will be limitations on the range of interactions.
2950 By default,
{\tt mdrun
} tries to find the optimal balance between
2951 interaction range and efficiency. But it can happen that a simulation
2952 stops with an error message about missing interactions,
2953 or that a simulation might run slightly faster with shorter
2954 interaction ranges. A list of interaction ranges
2955 and their default values is given in
\tabref{dd_ranges
}.
2959 \begin{tabular
}{|c|c|ll|
}
2961 interaction & range & option & default \\
2963 non-bonded & $r_c$ = max($r_
{list
}$,$r_
{VdW
}$,$r_
{Coul
}$) &
{\tt mdp
} file & \\
2964 two-body bonded & max($r_
{mb
}$,$r_c$) &
{\tt mdrun -rdd
} & starting conf. +
10\% \\
2965 multi-body bonded & $r_
{mb
}$ &
{\tt mdrun -rdd
} & starting conf. +
10\% \\
2966 constraints & $r_
{con
}$ &
{\tt mdrun -rcon
} & est. from bond lengths \\
2967 virtual sites & $r_
{con
}$ &
{\tt mdrun -rcon
} &
0 \\
2971 \caption{The interaction ranges with domain decomposition.
}
2972 \label{tab:dd_ranges
}
2975 In most cases the defaults of
{\tt mdrun
} should not cause the simulation
2976 to stop with an error message of missing interactions.
2977 The range for the bonded interactions is determined from the distance
2978 between bonded charge-groups in the starting configuration, with
10\% added
2979 for headroom. For the constraints, the value of $r_
{con
}$ is determined by
2980 taking the maximum distance that (
{\tt lincs_order +
1}) bonds can cover
2981 when they all connect at angles of
120 degrees.
2982 The actual constraint communication is not limited by $r_
{con
}$,
2983 but by the minimum cell size $L_C$, which has the following lower limit:
2985 L_C
\geq \max(r_
{mb
},r_
{con
})
2987 Without dynamic load balancing the system is actually allowed to scale
2988 beyond this limit when pressure scaling is used.
2989 {\bf Note
} that for triclinic boxes, $L_C$ is not simply the box diagonal
2990 component divided by the number of cells in that direction,
2991 rather it is the shortest distance between the triclinic cells borders.
2992 For rhombic dodecahedra this is a factor of $
\sqrt{3/
2}$ shorter
2995 When $r_
{mb
} > r_c$,
{\tt mdrun
} employs a smart algorithm to reduce
2996 the communication. Simply communicating all charge groups within
2997 $r_
{mb
}$ would increase the amount of communication enormously.
2998 Therefore only charge-groups that are connected by bonded interactions
2999 to charge groups which are not locally present are communicated.
3000 This leads to little extra communication, but also to a slightly
3001 increased cost for the domain decomposition setup.
3002 In some cases,
{\eg} coarse-grained simulations with a very short cut-off,
3003 one might want to set $r_
{mb
}$ by hand to reduce this cost.
3005 \subsection{Multiple-Program, Multiple-Data PME parallelization
\index{PME
}}
3006 \label{subsec:mpmd_pme
}
3007 Electrostatics interactions are long-range, therefore special
3008 algorithms are used to avoid summation over many atom pairs.
3009 In
{\gromacs} this is usually
3010 % \ifthenelse{\equal{\gmxlite}{1}}
3012 {PME (
\secref{pme
}).
}
3013 Since with PME all particles interact with each other, global communication
3014 is required. This will usually be the limiting factor for
3015 scaling with domain decomposition.
3016 To reduce the effect of this problem, we have come up with
3017 a Multiple-Program, Multiple-Data approach~
\cite{Hess2008b
}.
3018 Here, some processors are selected to do only the PME mesh calculation,
3019 while the other processors, called particle-particle (PP) nodes,
3020 do all the rest of the work.
3021 For rectangular boxes the optimal PP to PME node ratio is usually
3:
1,
3022 for rhombic dodecahedra usually
2:
1.
3023 When the number of PME nodes is reduced by a factor of
4, the number
3024 of communication calls is reduced by about a factor of
16.
3025 Or put differently, we can now scale to
4 times more nodes.
3026 In addition, for modern
4 or
8 core machines in a network,
3027 the effective network bandwidth for PME is quadrupled,
3028 since only a quarter of the cores will be using the network connection
3029 on each machine during the PME calculations.
3032 \centerline{\includegraphics[width=
12cm
]{plots/mpmd-pme
}}
3034 Example of
8 nodes without (left) and with (right) MPMD.
3035 The PME communication (red arrows) is much higher on the left
3036 than on the right. For MPMD additional PP - PME coordinate
3037 and force communication (blue arrows) is required,
3038 but the total communication complexity is lower.
3039 \label{fig:mpmd_pme
}
3043 {\tt mdrun
} will by default interleave the PP and PME nodes.
3044 If the processors are not number consecutively inside the machines,
3045 one might want to use
{\tt mdrun -ddorder pp_pme
}.
3046 For machines with a real
3-D torus and proper communication software
3047 that assigns the processors accordingly one should use
3048 {\tt mdrun -ddorder cartesian
}.
3050 To optimize the performance one should usually set up the cut-offs
3051 and the PME grid such that the PME load is
25 to
33\% of the total
3052 calculation load.
{\tt grompp
} will print an estimate for this load
3053 at the end and also
{\tt mdrun
} calculates the same estimate
3054 to determine the optimal number of PME nodes to use.
3055 For high parallelization it might be worthwhile to optimize
3056 the PME load with the
{\tt mdp
} settings and/or the number
3057 of PME nodes with the
{\tt -npme
} option of
{\tt mdrun
}.
3058 For changing the electrostatics settings it is useful to know
3059 the accuracy of the electrostatics remains nearly constant
3060 when the Coulomb cut-off and the PME grid spacing are scaled
3062 {\bf Note
} that it is usually better to overestimate than to underestimate
3063 the number of PME nodes, since the number of PME nodes is smaller
3064 than the number of PP nodes, which leads to less total waiting time.
3066 The PME domain decomposition can be
1-D or
2-D along the $x$ and/or
3067 $y$ axis.
2-D decomposition is also known as
\normindex{pencil decomposition
} because of
3068 the shape of the domains at high parallelization.
3069 1-D decomposition along the $y$ axis can only be used when
3070 the PP decomposition has only
1 domain along $x$.
2-D PME decomposition
3071 has to have the number of domains along $x$ equal to the number of
3072 the PP decomposition.
{\tt mdrun
} automatically chooses
1-D or
2-D
3073 PME decomposition (when possible with the total given number of nodes),
3074 based on the minimum amount of communication for the coordinate redistribution
3075 in PME plus the communication for the grid overlap and transposes.
3076 To avoid superfluous communication of coordinates and forces
3077 between the PP and PME nodes, the number of DD cells in the $x$
3078 direction should ideally be the same or a multiple of the number
3079 of PME nodes. By default,
{\tt mdrun
} takes care of this issue.
3081 \subsection{Domain decomposition flow chart
}
3082 In
\figref{dd_flow
} a flow chart is shown for domain decomposition
3083 with all possible communication for different algorithms.
3084 For simpler simulations, the same flow chart applies,
3085 without the algorithms and communication for
3086 the algorithms that are not used.
3089 \centerline{\includegraphics[width=
12cm
]{plots/flowchart
}}
3091 Flow chart showing the algorithms and communication (arrows)
3092 for a standard MD simulation with virtual sites, constraints
3093 and separate PME-mesh nodes.
3099 \section{Implicit solvation
\index{implicit solvation
}\index{Generalized Born methods
}}
3101 Implicit solvent models provide an efficient way of representing
3102 the electrostatic effects of solvent molecules, while saving a
3103 large piece of the computations involved in an accurate, aqueous
3104 description of the surrounding water in molecular dynamics simulations.
3105 Implicit solvation models offer several advantages compared with
3106 explicit solvation, including eliminating the need for the equilibration of water
3107 around the solute, and the absence of viscosity, which allows the protein
3108 to more quickly explore conformational space.
3110 Implicit solvent calculations in
{\gromacs} can be done using the
3111 generalized Born-formalism, and the Still~
\cite{Still97
}, HCT~
\cite{Truhlar96
},
3112 and OBC~
\cite{Case04
} models are available for calculating the Born radii.
3114 Here, the free energy $G_
{solv
}$ of solvation is the sum of three terms,
3115 a solvent-solvent cavity term ($G_
{cav
}$), a solute-solvent van der
3116 Waals term ($G_
{vdw
}$), and finally a solvent-solute electrostatics
3117 polarization term ($G_
{pol
}$).
3119 The sum of $G_
{cav
}$ and $G_
{vdw
}$ corresponds to the (non-polar)
3120 free energy of solvation for a molecule from which all charges
3121 have been removed, and is commonly called $G_
{np
}$,
3122 calculated from the total solvent accessible surface area
3123 multiplied with a surface tension.
3124 The total expression for the solvation free energy then becomes:
3127 G_
{solv
} = G_
{np
} + G_
{pol
}
3131 Under the generalized Born model, $G_
{pol
}$ is calculated from the generalized Born equation~
\cite{Still97
}:
3134 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)
}}
3135 \label{eqn:gb_still
}
3138 In
{\gromacs}, we have introduced the substitution~
\cite{Larsson10
}:
3141 c_i=
\frac{1}{\sqrt{b_i
}}
3142 \label{eqn:gb_subst
}
3145 which makes it possible to introduce a cheap transformation to a new
3146 variable $x$ when evaluating each interaction, such that:
3149 x=
\frac{r_
{ij
}}{\sqrt{b_i b_j
}} = r_
{ij
} c_i c_j
3150 \label{eqn:gb_subst2
}
3153 In the end, the full re-formulation of~
\ref{eqn:gb_still
} becomes:
3156 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)
3157 \label{eqn:gb_final
}
3160 The non-polar part ($G_
{np
}$) of Equation~
\ref{eqn:gb_solv
} is calculated
3161 directly from the Born radius of each atom using a simple ACE type
3162 approximation by Schaefer
{\em et al.
}~
\cite{Karplus98
}, including a
3163 simple loop over all atoms.
3164 This requires only one extra solvation parameter, independent of atom type,
3165 but differing slightly between the three Born radii models.
3167 % LocalWords: GROningen MAchine BIOSON Groningen GROMACS Berendsen der Spoel
3168 % LocalWords: Drunen Comp Phys Comm ROck NS FFT pbc EM ifthenelse gmxlite ff
3169 % LocalWords: octahedra triclinic Ewald PME PPPM trjconv xy solvated
3170 % LocalWords: boxtypes boxshapes editconf Lennard mdpopt COM XTC TNG kT defunits
3171 % LocalWords: Boltzmann's Mueller nb int mdrun chargegroup simplerc prefactor
3172 % LocalWords: pme waterloops CH NH CO df com virial integrator Verlet vverlet
3173 % LocalWords: integrators ref timepoint timestep timesteps mdp md vv avek NVE
3174 % LocalWords: NVT off's leapfrogv lll LR rmfast SPC fs Nos physicality ps GMX
3175 % LocalWords: Tcoupling nonergodic thermostatting NOSEHOOVER algorithmes ij yx
3176 % LocalWords: Parrinello Rahman rescales atm anisotropically ccc xz zx yy yz
3177 % LocalWords: zy zz se barostat compressibilities MTTK NPT Martyna al isobaric
3178 % LocalWords: Tuckerman vir PV fkT iLt iL Liouville NHC Eq baro mu trj mol bc
3179 % LocalWords: freezegroup Shannon's polarizability Overhauser barostats iLn KE
3180 % LocalWords: negligibly thermostatted Tobias rhombic maxwell et xtc tng TC rlist
3181 % LocalWords: waals LINCS holonomic plincs lincs unc ang SA Langevin SD amu BD
3182 % LocalWords: bfgs Broyden Goldfarb Shanno mkT kJ DFLEXIBLE Nocedal diag nmeig
3183 % LocalWords: diagonalization anaeig nmens covanal ddg feia BT dp dq pV dV dA
3184 % LocalWords: NpT eq stepsize REMD constrainted website Okabe MPI covar edi dd
3185 % LocalWords: progman NMR ddcells innerloops ddtric tric dds rdd conf rcon est
3186 % LocalWords: mb PP MPMD ddorder pp cartesian grompp npme parallelizable edr
3187 % LocalWords: macromolecule nstlist vacuo parallelization dof indices MBAR AVX
3188 % LocalWords: TOL numerics parallelized eigenvectors dG parallelepipeds VdW np
3189 % LocalWords: Coul multi solvation HCT OBC solv cav vdw Schaefer symplectic dt
3190 % LocalWords: pymbar multinode subensemble Monte solute subst groupconcept GPU
3191 % LocalWords: dodecahedron octahedron dodecahedra equilibration usinggroups nm
3192 % LocalWords: topologies rlistlong CUDA GPUs rcoulomb SIMD BlueGene FPUs erfc
3193 % LocalWords: cutoffschemesupport unbuffered bondeds AdResS OpenMP ewald rtol
3194 % LocalWords: verletdrift peptide RMS rescaling ergodicity ergodic discretized
3195 % LocalWords: isothermal compressibility isotropically anisotropic iteratively
3196 % LocalWords: incompressible integrations translational biomolecules NMA PCA
3197 % LocalWords: Bennett's equilibrated Hamiltonians covariance equilibrate
3198 % LocalWords: inhomogeneous conformational online other's th