Fix SIMD detection on new AMD AVX CPUs w/o fma
[gromacs/AngularHB.git] / docs / manual / algorithms.tex
blobcaaafabcf242086498f0f3b87448bb0937697554
2 % This file is part of the GROMACS molecular simulation package.
4 % Copyright (c) 2013,2014,2015,2016, 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}}
54 \chapter{Algorithms}
55 \label{ch:algorithms}
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
64 this.
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
70 %treated in detail.
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}}
87 \label{sec:pbc}
88 \begin{figure}
89 \centerline{\includegraphics[width=9cm]{plots/pbctric}}
90 \caption {Periodic boundary conditions in two dimensions.}
91 \label{fig:pbc}
92 \end{figure}
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:
130 \beq
131 \label{eqn:box_rot}
132 a_y = a_z = b_z = 0
133 \eeq
134 \beq
135 \label{eqn:box_shift1}
136 a_x>0,~~~~b_y>0,~~~~c_z>0
137 \eeq
138 \beq
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
143 \eeq
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.
160 \begin{figure}
161 \centerline{
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}
168 \end{figure}
170 \subsection{Some useful box types}
171 \begin{table}
172 \centerline{
173 \begin{tabular}{|c|c|c|ccc|ccc|}
174 \dline
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}$ \\
178 \dline
179 & & & $d$ & 0 & 0 & & & \\
180 cubic & $d$ & $d^3$ & 0 & $d$ & 0 & $90^\circ$ & $90^\circ$ & $90^\circ$ \\
181 & & & 0 & 0 & $d$ & & & \\
182 \hline
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$ & & & \\
186 \hline
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$ & & & \\
190 \hline
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$ & & & \\
194 \dline
195 \end{tabular}
197 \caption{The cubic box, the rhombic \normindex{dodecahedron} and the truncated
198 \normindex{octahedron}.}
199 \label{tab:boxtypes}
200 \end{table}
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
224 vector:
225 \beq
226 \label{eqn:physicalrc}
227 R_c < \half \min(\|{\bf a}\|,\|{\bf b}\|,\|{\bf c}\|),
228 \eeq
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:
240 \beq
241 \label{eqn:gridrc}
242 R_c < \min(a_x,b_y,c_z)
243 \eeq
244 For simple search the extra restriction is stronger:
245 \beq
246 \label{eqn:simplerc}
247 R_c < \half \min(a_x,b_y,c_z)
248 \eeq
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
256 considered.
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:
265 \begin{description}
266 \item[\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
275 groups below.
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
308 sparingly.
310 All non-bonded interactions between pairs of energy-monitor groups can
311 be excluded\index{exclusions}
312 \ifthenelse{\equal{\gmxlite}{1}}
314 {(see details in the User Guide).}
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
319 are not required.
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
327 removal.
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.
337 \end{description}
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}
343 \label{sec:MD}
344 \begin{figure}
345 \begin{center}
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 \\
354 $\Downarrow$\\
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 may be computed. \\
366 $\Downarrow$\\
367 {\bf 3. Update configuration} \\[1ex]
368 The movement of the atoms is simulated by numerically solving Newton's
369 equations of motion \\[1ex]
370 $\displaystyle
371 \frac {\de^2\ve{r}_i}{\de t^2} = \frac{\ve{F}_i}{m_i} $ \\
372 or \\
373 $\displaystyle
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]
376 $\Downarrow$ \\
377 {\bf 4.} if required: {\bf Output step} \\
378 write positions, velocities, energies, temperature, pressure, etc. \\
379 \end{shadowenv}
380 \caption{The global MD algorithm}
381 \label{fig:global}
382 \end{center}
383 \end{figure}
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
393 be read in.
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}
401 \begin{figure}
402 \centerline{\includegraphics[width=8cm]{plots/maxwell}}
403 \caption{A Maxwell-Boltzmann velocity distribution, generated from
404 random numbers.}
405 \label{fig:maxwell}
406 \end{figure}
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$:
420 \beq
421 p(v_i) = \sqrt{\frac{m_i}{2 \pi kT}}\exp\left(-\frac{m_i v_i^2}{2kT}\right)
422 \eeq
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}}
450 \label{subsec:ns}
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 particles
476 $i$, a displacement vector for particle $i$, and all particles $j$ that
477 are within \verb'rlist' of this particular image of particle $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 particle 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, but on the most recent x86
500 processors, this advantage is negated by the better instruction-level
501 parallelism available in the Verlet-scheme implementation. The group
502 scheme is deprecated in version 5.0, and will be removed in a future
503 version. For practical details of choosing and setting up
504 cut-off schemes, please see the User Guide.
506 In the group scheme, a neighbor list is generated consisting of pairs
507 of groups of at least one particle. These groups were originally
508 \swapindex{charge}{group}s \ifthenelse{\equal{\gmxlite}{1}}{}{(see
509 \secref{chargegroup})}, but with a proper treatment of long-range
510 electrostatics, performance in unbuffered simulations is their only advantage. A pair of groups
511 is put into the neighbor list when their center of geometry is within
512 the cut-off distance. Interactions between all particle pairs (one from
513 each charge group) are calculated for a certain number of MD steps,
514 until the neighbor list is updated. This setup is efficient, as the
515 neighbor search only checks distance between charge-group pair, not
516 particle pairs (saves a factor of $3 \times 3 = 9$ with a three-particle water
517 model) and the non-bonded force kernels can be optimized for, say, a
518 water molecule ``group''. Without explicit buffering, this setup leads
519 to energy drift as some particle pairs which are within the cut-off don't
520 interact and some outside the cut-off do interact. This can be caused
522 \begin{itemize}
523 \item particles moving across the cut-off between neighbor search steps, and/or
524 \item for charge groups consisting of more than one particle, particle pairs
525 moving in/out of the cut-off when their charge group center of
526 geometry distance is outside/inside of the cut-off.
527 \end{itemize}
528 Explicitly adding a buffer to the neighbor list will remove such
529 artifacts, but this comes at a high computational cost. How severe the
530 artifacts are depends on the system, the properties in which you are
531 interested, and the cut-off setup.
533 The Verlet cut-off scheme uses a buffered pair list by default. It
534 also uses clusters of particles, but these are not static as in the group
535 scheme. Rather, the clusters are defined spatially and consist of 4 or
536 8 particles, which is convenient for stream computing, using e.g. SSE, AVX
537 or CUDA on GPUs. At neighbor search steps, a pair list is created
538 with a Verlet buffer, ie. the pair-list cut-off is larger than the
539 interaction cut-off. In the non-bonded kernels, interactions are only
540 computed when a particle pair is within the cut-off distance at that
541 particular time step. This ensures that as particles move between pair
542 search steps, forces between nearly all particles within the cut-off
543 distance are calculated. We say {\em nearly} all particles, because
544 {\gromacs} uses a fixed pair list update frequency for
545 efficiency. A particle-pair, whose distance was outside the cut-off,
546 could possibly move enough during this fixed number of
547 steps that its distance is now within the cut-off. This
548 small chance results in a small energy drift, and the size of the
549 chance depends on the temperature. When temperature
550 coupling is used, the buffer size can be determined automatically,
551 given a certain tolerance on the energy drift.
553 The Verlet cut-off scheme is implemented in a very efficient fashion
554 based on clusters of particles. The simplest example is a cluster size
555 of 4 particles. The pair list is then constructed based on cluster
556 pairs. The cluster-pair search is much faster searching based on
557 particle pairs, because $4 \times 4 = 16$ particle pairs are put in
558 the list at once. The non-bonded force calculation kernel can then
559 calculate many particle-pair interactions at once, which maps nicely
560 to SIMD or SIMT units on modern hardware, which can perform multiple
561 floating operations at once. These non-bonded kernels
562 are much faster than the kernels used in the group scheme for most
563 types of systems, particularly on newer hardware.
565 \ifthenelse{\equal{\gmxlite}{1}}{}{
566 \subsubsection{Energy drift and pair-list buffering}
567 For a canonical (NVT) ensemble, the average energy error caused by
568 diffusion of $j$ particles from outside the pair-list cut-off
569 $r_\ell$ to inside the interaction cut-off $r_c$ over the lifetime
570 of the list can be determined from the atomic
571 displacements and the shape of the potential at the cut-off.
572 %Since we are interested in the small drift regime, we will assume
573 %#that atoms will only move within the cut-off distance in the last step,
574 %$n_\mathrm{ps}-1$, of the pair list update interval $n_\mathrm{ps}$.
575 %Over this number of steps the displacment of an atom with mass $m$
576 The displacement distribution along one dimension for a freely moving
577 particle with mass $m$ over time $t$ at temperature $T$ is
578 a Gaussian $G(x)$
579 of zero mean and variance $\sigma^2 = t^2 k_B T/m$. For the distance
580 between two particles, the variance changes to $\sigma^2 = \sigma_{12}^2 =
581 t^2 k_B T(1/m_1+1/m_2)$. Note that in practice particles usually
582 interact with other particles over time $t$ and therefore the real
583 displacement distribution is much narrower. Given a non-bonded
584 interaction cut-off distance of $r_c$ and a pair-list cut-off
585 $r_\ell=r_c+r_b$ for $r_b$ the Verlet buffer size, we can then
586 write the average energy error after time $t$ for all missing pair
587 interactions between a single $i$ particle of type 1 surrounded
588 by all $j$ particles that are of type 2 with number density $\rho_2$,
589 when the inter-particle distance changes from $r_0$ to $r_t$, as:
590 \beq
591 \langle \Delta V \rangle =
592 \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
593 \eeq
594 To evaluate this analytically, we need to make some approximations. First we replace $V(r_t)$ by a Taylor expansion around $r_c$, then we can move the lower bound of the integral over $r_0$ to $-\infty$ which will simplify the result:
595 \begin{eqnarray}
596 \langle \Delta V \rangle &\approx&
597 \int_{-\infty}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 \Big[ V'(r_c) (r_t - r_c) +
598 \nonumber\\
600 \phantom{\int_{-\infty}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 \Big[}
601 V''(r_c)\frac{1}{2}(r_t - r_c)^2 +
602 \nonumber\\
604 \phantom{\int_{-\infty}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 \Big[}
605 V'''(r_c)\frac{1}{6}(r_t - r_c)^3 +
606 \nonumber\\
608 \phantom{\int_{-\infty}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 \Big[}
609 O \! \left((r_t - r_c)^4 \right)\Big] G\!\left(\frac{r_t-r_0}{\sigma}\right) d r_0 \, d r_t
610 \end{eqnarray}
611 Replacing the factor $r_0^2$ by $(r_\ell + \sigma)^2$, which results in a slight overestimate, allows us to calculate the integrals analytically:
612 \begin{eqnarray}
613 \langle \Delta V \rangle \!
614 &\approx&
615 4 \pi (r_\ell+\sigma)^2 \rho_2
616 \int_{-\infty}^{r_c} \int_{r_\ell}^\infty \Big[ V'(r_c) (r_t - r_c) +
617 \nonumber\\
619 \phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \int_{-\infty}^{r_c} \int_{r_\ell}^\infty \Big[}
620 V''(r_c)\frac{1}{2}(r_t - r_c)^2 +
621 \nonumber\\
623 \phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \int_{-\infty}^{r_c} \int_{r_\ell}^\infty \Big[}
624 V'''(r_c)\frac{1}{6}(r_t - r_c)^3 \Big] G\!\left(\frac{r_t-r_0}{\sigma}\right)
625 d r_0 \, d r_t\\
627 4 \pi (r_\ell+\sigma)^2 \rho_2 \bigg\{
628 \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] +
629 \nonumber\\
631 \phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \bigg\{ }
632 \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] +
633 \nonumber\\
635 \phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \bigg\{ }
636 \frac{1}{24}V'''(r_c)\bigg[ r_b\sigma(r_b^2+5\sigma^2) G\!\left(\frac{r_b}{\sigma}\right)
637 \nonumber\\
639 \phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \bigg\{ \frac{1}{24}V'''(r_c)\bigg[ }
640 - (r_b^4+6r_b^2\sigma^2+3\sigma^4 ) E\!\left(\frac{r_b}{\sigma}\right) \bigg]
641 \bigg\}
642 \end{eqnarray}
644 where $G(x)$ is a Gaussian distribution with 0 mean and unit variance and
645 $E(x)=\frac{1}{2}\mathrm{erfc}(x/\sqrt{2})$. We always want to achieve
646 small energy error, so $\sigma$ will be small compared to both $r_c$
647 and $r_\ell$, thus the approximations in the equations above are good,
648 since the Gaussian distribution decays rapidly. The energy error needs
649 to be averaged over all particle pair types and weighted with the
650 particle counts. In {\gromacs} we don't allow cancellation of error
651 between pair types, so we average the absolute values. To obtain the
652 average energy error per unit time, it needs to be divided by the
653 neighbor-list life time $t = ({\tt nstlist} - 1)\times{\tt dt}$. The
654 function can not be inverted analytically, so we use bisection to
655 obtain the buffer size $r_b$ for a target drift. Again we note that
656 in practice the error we usually be much smaller than this estimate,
657 as in the condensed phase particle displacements will be much smaller
658 than for freely moving particles, which is the assumption used here.
660 When (bond) constraints are present, some particles will have fewer
661 degrees of freedom. This will reduce the energy errors. For simplicity,
662 we only consider one constraint per particle, the heaviest particle
663 in case a particle is involved in multiple constraints.
664 This simplification overestimates the displacement. The motion of
665 a constrained particle is a superposition of the 3D motion of the
666 center of mass of both particles and a 2D rotation around the center of mass.
667 The displacement in an arbitrary direction of a particle with 2 degrees
668 of freedom is not Gaussian, but rather follows the complementary error
669 function:
670 \beq
671 \frac{\sqrt{\pi}}{2\sqrt{2}\sigma}\,\mathrm{erfc}\left(\frac{|r|}{\sqrt{2}\,\sigma}\right)
672 \label{eqn:2D_disp}
673 \eeq
674 where $\sigma^2$ is again $t^2 k_B T/m$. This distribution can no
675 longer be integrated analytically to obtain the energy error. But we
676 can generate a tight upper bound using a scaled and shifted Gaussian
677 distribution (not shown). This Gaussian distribution can then be used
678 to calculate the energy error as described above. The rotation displacement
679 around the center of mass can not be more than the length of the arm.
680 To take this into account, we scale $\sigma$ in \eqnref{2D_disp} (details
681 not presented here) to obtain an overestimate of the real displacement.
682 This latter effect significantly reduces the buffer size for longer
683 neighborlist lifetimes in e.g. water, as constrained hydrogens are by far
684 the fastest particles, but they can not move further than 0.1 nm
685 from the heavy atom they are connected to.
688 There is one important implementation detail that reduces the energy
689 errors caused by the finite Verlet buffer list size. The derivation
690 above assumes a particle pair-list. However, the {\gromacs}
691 implementation uses a cluster pair-list for efficiency. The pair list
692 consists of pairs of clusters of 4 particles in most cases, also
693 called a $4 \times 4$ list, but the list can also be $4 \times 8$ (GPU
694 CUDA kernels and AVX 256-bit single precision kernels) or $4 \times 2$
695 (SSE double-precision kernels). This means that the pair-list is
696 effectively much larger than the corresponding $1 \times 1$ list. Thus
697 slightly beyond the pair-list cut-off there will still be a large
698 fraction of particle pairs present in the list. This fraction can be
699 determined in a simulation and accurately estimated under some
700 reasonable assumptions. The fraction decreases with increasing
701 pair-list range, meaning that a smaller buffer can be used. For
702 typical all-atom simulations with a cut-off of 0.9 nm this fraction is
703 around 0.9, which gives a reduction in the energy errors of a factor of
704 10. This reduction is taken into account during the automatic Verlet
705 buffer calculation and results in a smaller buffer size.
707 \begin{figure}
708 \centerline{\includegraphics[width=9cm]{plots/verlet-drift}}
709 \caption {Energy drift per atom for an SPC/E water system at 300K with
710 a time step of 2 fs and a pair-list update period of 10 steps
711 (pair-list life time: 18 fs). PME was used with {\tt ewald-rtol} set
712 to 10$^{-5}$; this parameter affects the shape of the potential at
713 the cut-off. Error estimates due to finite Verlet buffer size are
714 shown for a $1 \times 1$ atom pair list and $4 \times 4$ atom pair
715 list without and with (dashed line) cancellation of positive and
716 negative errors. Real energy drift is shown for simulations using
717 double- and mixed-precision settings. Rounding errors in the SETTLE
718 constraint algorithm from the use of single precision causes
719 the drift to become negative
720 at large buffer size. Note that at zero buffer size, the real drift
721 is small because positive (H-H) and negative (O-H) energy errors
722 cancel.}
723 \label{fig:verletdrift}
724 \end{figure}
726 In \figref{verletdrift} one can see that for small buffer sizes the drift
727 of the total energy is much smaller than the pair energy error tolerance,
728 due to cancellation of errors. For larger buffer size, the error estimate
729 is a factor of 6 higher than drift of the total energy, or alternatively
730 the buffer estimate is 0.024 nm too large. This is because the protons
731 don't move freely over 18 fs, but rather vibrate.
732 %At a buffer size of zero there is cancellation of
733 %drift due to repulsive (H-H) and attractive (O-H) interactions.
735 \subsubsection{Cut-off artifacts and switched interactions}
736 With the Verlet scheme, the pair potentials are shifted to be zero at
737 the cut-off, which makes the potential the integral of the force.
738 This is only possible in the group scheme if the shape of the potential
739 is such that its value is zero at the cut-off distance.
740 However, there can still be energy drift when the
741 forces are non-zero at the cut-off. This effect is extremely small and
742 often not noticeable, as other integration errors (e.g. from constraints)
743 may dominate. To
744 completely avoid cut-off artifacts, the non-bonded forces can be
745 switched exactly to zero at some distance smaller than the neighbor
746 list cut-off (there are several ways to do this in {\gromacs}, see
747 \secref{mod_nb_int}). One then has a buffer with the size equal to the
748 neighbor list cut-off less the longest interaction cut-off.
750 } % Brace matches ifthenelse test for gmxlite
752 \subsubsection{Simple search\swapindexquiet{simple}{search}}
753 Due to \eqnsref{box_rot}{simplerc}, the vector $\rvij$
754 connecting images within the cut-off $R_c$ can be found by constructing:
755 \bea
756 \ve{r}''' & = & \ve{r}_j-\ve{r}_i \\
757 \ve{r}'' & = & \ve{r}''' - {\bf c}*\verb'round'(r'''_z/c_z) \\
758 \ve{r}' & = & \ve{r}'' - {\bf b}*\verb'round'(r''_y/b_y) \\
759 \ve{r}_{ij} & = & \ve{r}' - {\bf a}*\verb'round'(r'_x/a_x)
760 \eea
761 When distances between two particles in a triclinic box are needed
762 that do not obey \eqnref{box_rot},
763 many shifts of combinations of box vectors need to be considered to find
764 the nearest image.
766 \ifthenelse{\equal{\gmxlite}{1}}{}{
768 \begin{figure}
769 \centerline{\includegraphics[width=8cm]{plots/nstric}}
770 \caption {Grid search in two dimensions. The arrows are the box vectors.}
771 \label{fig:grid}
772 \end{figure}
774 \subsubsection{Grid search\swapindexquiet{grid}{search}}
775 \label{sec:nsgrid}
776 The grid search is schematically depicted in \figref{grid}. All
777 particles are put on the {\nsgrid}, with the smallest spacing $\ge$
778 $R_c/2$ in each of the directions. In the direction of each box
779 vector, a particle $i$ has three images. For each direction the image
780 may be -1,0 or 1, corresponding to a translation over -1, 0 or +1 box
781 vector. We do not search the surrounding {\nsgrid} cells for neighbors
782 of $i$ and then calculate the image, but rather construct the images
783 first and then search neighbors corresponding to that image of $i$.
784 As \figref{grid} shows, some grid cells may be searched more than once
785 for different images of $i$. This is not a problem, since, due to the
786 minimum image convention, at most one image will ``see'' the
787 $j$-particle. For every particle, fewer than 125 (5$^3$) neighboring
788 cells are searched. Therefore, the algorithm scales linearly with the
789 number of particles. Although the prefactor is large, the scaling
790 behavior makes the algorithm far superior over the standard $O(N^2)$
791 algorithm when there are more than a few hundred particles. The
792 grid search is equally fast for rectangular and triclinic boxes. Thus
793 for most protein and peptide simulations the rhombic dodecahedron will
794 be the preferred box shape.
795 } % Brace matches ifthenelse test for gmxlite
797 \ifthenelse{\equal{\gmxlite}{1}}{}{
798 \subsubsection{Charge groups}
799 \label{sec:chargegroup}\swapindexquiet{charge}{group}%
800 Charge groups were originally introduced to reduce cut-off artifacts
801 of Coulomb interactions. When a plain cut-off is used, significant
802 jumps in the potential and forces arise when atoms with (partial) charges
803 move in and out of the cut-off radius. When all chemical moieties have
804 a net charge of zero, these jumps can be reduced by moving groups
805 of atoms with net charge zero, called charge groups, in and
806 out of the neighbor list. This reduces the cut-off effects from
807 the charge-charge level to the dipole-dipole level, which decay
808 much faster. With the advent of full range electrostatics methods,
809 such as particle-mesh Ewald (\secref{pme}), the use of charge groups is
810 no longer required for accuracy. It might even have a slight negative effect
811 on the accuracy or efficiency, depending on how the neighbor list is made
812 and the interactions are calculated.
814 But there is still an important reason for using ``charge groups'': efficiency with the group cut-off scheme.
815 Where applicable, neighbor searching is carried out on the basis of
816 charge groups which are defined in the molecular topology.
817 If the nearest image distance between the {\em
818 geometrical centers} of the atoms of two charge groups is less than
819 the cut-off radius, all atom pairs between the charge groups are
820 included in the pair list.
821 The neighbor searching for a water system, for instance,
822 is $3^2=9$ times faster when each molecule is treated as a charge group.
823 Also the highly optimized water force loops (see \secref{waterloops})
824 only work when all atoms in a water molecule form a single charge group.
825 Currently the name {\em neighbor-search group} would be more appropriate,
826 but the name charge group is retained for historical reasons.
827 When developing a new force field, the advice is to use charge groups
828 of 3 to 4 atoms for optimal performance. For all-atom force fields
829 this is relatively easy, as one can simply put hydrogen atoms, and in some
830 case oxygen atoms, in the same charge group as the heavy atom they
831 are connected to; for example: CH$_3$, CH$_2$, CH, NH$_2$, NH, OH, CO$_2$, CO.
833 With the Verlet cut-off scheme, charge groups are ignored.
835 } % Brace matches ifthenelse test for gmxlite
837 \subsection{Compute forces}
838 \label{subsec:forces}
840 \subsubsection{Potential energy}
841 When forces are computed, the \swapindex{potential}{energy} of each
842 interaction term is computed as well. The total potential energy is
843 summed for various contributions, such as Lennard-Jones, Coulomb, and
844 bonded terms. It is also possible to compute these contributions for
845 {\em energy-monitor groups} of atoms that are separately defined (see
846 \secref{groupconcept}).
848 \subsubsection{Kinetic energy and temperature}
849 The \normindex{temperature} is given by the total
850 \swapindex{kinetic}{energy} of the $N$-particle system:
851 \beq
852 E_{kin} = \half \sum_{i=1}^N m_i v_i^2
853 \eeq
854 From this the absolute temperature $T$ can be computed using:
855 \beq
856 \half N_{\mathrm{df}} kT = E_{\mathrm{kin}}
857 \label{eqn:E-T}
858 \eeq
859 where $k$ is Boltzmann's constant and $N_{df}$ is the number of
860 degrees of freedom which can be computed from:
861 \beq
862 N_{\mathrm{df}} ~=~ 3 N - N_c - N_{\mathrm{com}}
863 \eeq
864 Here $N_c$ is the number of {\em \normindex{constraints}} imposed on the system.
865 When performing molecular dynamics $N_{\mathrm{com}}=3$ additional degrees of
866 freedom must be removed, because the three
867 center-of-mass velocities are constants of the motion, which are usually
868 set to zero. When simulating in vacuo, the rotation around the center of mass
869 can also be removed, in this case $N_{\mathrm{com}}=6$.
870 When more than one temperature-coupling group\index{temperature-coupling group} is used, the number of degrees
871 of freedom for group $i$ is:
872 \beq
873 N^i_{\mathrm{df}} ~=~ (3 N^i - N^i_c) \frac{3 N - N_c - N_{\mathrm{com}}}{3 N - N_c}
874 \eeq
876 The kinetic energy can also be written as a tensor, which is necessary
877 for pressure calculation in a triclinic system, or systems where shear
878 forces are imposed:
879 \beq
880 {\bf E}_{\mathrm{kin}} = \half \sum_i^N m_i \vvi \otimes \vvi
881 \eeq
883 \subsubsection{Pressure and virial}
884 The \normindex{pressure}
885 tensor {\bf P} is calculated from the difference between
886 kinetic energy $E_{\mathrm{kin}}$ and the \normindex{virial} ${\bf \Xi}$:
887 \beq
888 {\bf P} = \frac{2}{V} ({\bf E}_{\mathrm{kin}}-{\bf \Xi})
889 \label{eqn:P}
890 \eeq
891 where $V$ is the volume of the computational box.
892 The scalar pressure $P$, which can be used for pressure coupling in the case
893 of isotropic systems, is computed as:
894 \beq
895 P = {\rm trace}({\bf P})/3
896 \eeq
898 The virial ${\bf \Xi}$ tensor is defined as:
899 \beq
900 {\bf \Xi} = -\half \sum_{i<j} \rvij \otimes \Fvij
901 \label{eqn:Xi}
902 \eeq
904 \ifthenelse{\equal{\gmxlite}{1}}{}{
905 The {\gromacs} implementation of the virial computation is described
906 in \secref{virial}.
907 } % Brace matches ifthenelse test for gmxlite
910 \subsection{The \swapindex{leap-frog}{integrator}}
911 \label{subsec:update}
912 \begin{figure}
913 \centerline{\includegraphics[width=8cm]{plots/leapfrog}}
914 \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
915 like frogs over each other's backs.}
916 \label{fig:leapfrog}
917 \end{figure}
919 The default MD integrator in {\gromacs} is the so-called {\em leap-frog}
920 algorithm~\cite{Hockney74} for the integration of the equations of
921 motion. When extremely accurate integration with temperature
922 and/or pressure coupling is required, the velocity Verlet integrators
923 are also present and may be preferable (see \ssecref{vverlet}). The leap-frog
924 algorithm uses positions $\ve{r}$ at time $t$ and
925 velocities $\ve{v}$ at time $t-\hDt$; it updates positions and
926 velocities using the forces
927 $\ve{F}(t)$ determined by the positions at time $t$ using these relations:
928 \bea
929 \label{eqn:leapfrogv}
930 \ve{v}(t+\hDt) &~=~& \ve{v}(t-\hDt)+\frac{\Dt}{m}\ve{F}(t) \\
931 \ve{r}(t+\Dt) &~=~& \ve{r}(t)+\Dt\ve{v}(t+\hDt)
932 \eea
933 The algorithm is visualized in \figref{leapfrog}.
934 It produces trajectories that are identical to the Verlet~\cite{Verlet67} algorithm,
935 whose position-update relation is
936 \beq
937 \ve{r}(t+\Dt)~=~2\ve{r}(t) - \ve{r}(t-\Dt) + \frac{1}{m}\ve{F}(t)\Dt^2+O(\Dt^4)
938 \eeq
939 The algorithm is of third order in $\ve{r}$ and is time-reversible.
940 See ref.~\cite{Berendsen86b} for the merits of this algorithm and comparison
941 with other time integration algorithms.
943 The \swapindex{equations of}{motion} are modified for temperature
944 coupling and pressure coupling, and extended to include the
945 conservation of constraints, all of which are described below.
947 \subsection{The \swapindex{velocity Verlet}{integrator}}
948 \label{subsec:vverlet}
949 The velocity Verlet algorithm~\cite{Swope82} is also implemented in
950 {\gromacs}, though it is not yet fully integrated with all sets of
951 options. In velocity Verlet, positions $\ve{r}$ and velocities
952 $\ve{v}$ at time $t$ are used to integrate the equations of motion;
953 velocities at the previous half step are not required. \bea
954 \label{eqn:velocityverlet1}
955 \ve{v}(t+\hDt) &~=~& \ve{v}(t)+\frac{\Dt}{2m}\ve{F}(t) \\
956 \ve{r}(t+\Dt) &~=~& \ve{r}(t)+\Dt\,\ve{v}(t+\hDt) \\
957 \ve{v}(t+\Dt) &~=~& \ve{v}(t+\hDt)+\frac{\Dt}{2m}\ve{F}(t+\Dt)
958 \eea
959 or, equivalently,
960 \bea
961 \label{eqn:velocityverlet2}
962 \ve{r}(t+\Dt) &~=~& \ve{r}(t)+ \Dt\,\ve{v} + \frac{\Dt^2}{2m}\ve{F}(t) \\
963 \ve{v}(t+\Dt) &~=~& \ve{v}(t)+ \frac{\Dt}{2m}\left[\ve{F}(t) + \ve{F}(t+\Dt)\right]
964 \eea
965 With no temperature or pressure coupling, and with {\em corresponding}
966 starting points, leap-frog and velocity Verlet will generate identical
967 trajectories, as can easily be verified by hand from the equations
968 above. Given a single starting file with the {\em same} starting
969 point $\ve{x}(0)$ and $\ve{v}(0)$, leap-frog and velocity Verlet will
970 {\em not} give identical trajectories, as leap-frog will interpret the
971 velocities as corresponding to $t=-\hDt$, while velocity Verlet will
972 interpret them as corresponding to the timepoint $t=0$.
974 \subsection{Understanding reversible integrators: The Trotter decomposition}
975 To further understand the relationship between velocity Verlet and
976 leap-frog integration, we introduce the reversible Trotter formulation
977 of dynamics, which is also useful to understanding implementations of
978 thermostats and barostats in {\gromacs}.
980 A system of coupled, first-order differential equations can be evolved
981 from time $t = 0$ to time $t$ by applying the evolution operator
982 \bea
983 \Gamma(t) &=& \exp(iLt) \Gamma(0) \nonumber \\
984 iL &=& \dot{\Gamma}\cdot \nabla_{\Gamma},
985 \eea
986 where $L$ is the Liouville operator, and $\Gamma$ is the
987 multidimensional vector of independent variables (positions and
988 velocities).
989 A short-time approximation to the true operator, accurate at time $\Dt
990 = t/P$, is applied $P$ times in succession to evolve the system as
991 \beq
992 \Gamma(t) = \prod_{i=1}^P \exp(iL\Dt) \Gamma(0)
993 \eeq
994 For NVE dynamics, the Liouville operator is
995 \bea
996 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}.
997 \eea
998 This can be split into two additive operators
999 \bea
1000 iL_1 &=& \sum_{i=1}^N \frac{1}{m_i}\F(r_i) \cdot \nabla_{\vv_i} \nonumber \\
1001 iL_2 &=& \sum_{i=1}^{N} \vv_i \cdot \nabla_{\rv_i}
1002 \eea
1003 Then a short-time, symmetric, and thus reversible approximation of the true dynamics will be
1004 \bea
1005 \exp(iL\Dt) = \exp(iL_2\hDt) \exp(iL_1\Dt) \exp(iL_2\hDt) + \mathcal{O}(\Dt^3).
1006 \label{eq:NVE_Trotter}
1007 \eea
1008 This corresponds to velocity Verlet integration. The first
1009 exponential term over $\hDt$ corresponds to a velocity half-step, the
1010 second exponential term over $\Dt$ corresponds to a full velocity
1011 step, and the last exponential term over $\hDt$ is the final velocity
1012 half step. For future times $t = n\Dt$, this becomes
1013 \bea
1014 \exp(iLn\Dt) &\approx& \left(\exp(iL_2\hDt) \exp(iL_1\Dt) \exp(iL_2\hDt)\right)^n \nonumber \\
1015 &\approx& \exp(iL_2\hDt) \bigg(\exp(iL_1\Dt) \exp(iL_2\Dt)\bigg)^{n-1} \nonumber \\
1016 & & \;\;\;\; \exp(iL_1\Dt) \exp(iL_2\hDt)
1017 \eea
1018 This formalism allows us to easily see the difference between the
1019 different flavors of Verlet integrators. The leap-frog integrator can
1020 be seen as starting with Eq.~\ref{eq:NVE_Trotter} with the
1021 $\exp\left(iL_1 \dt\right)$ term, instead of the half-step velocity
1022 term, yielding
1023 \bea
1024 \exp(iLn\dt) &=& \exp\left(iL_1 \dt\right) \exp\left(iL_2 \Dt \right) + \mathcal{O}(\Dt^3).
1025 \eea
1026 Here, the full step in velocity is between $t-\hDt$ and $t+\hDt$,
1027 since it is a combination of the velocity half steps in velocity
1028 Verlet. For future times $t = n\Dt$, this becomes
1029 \bea
1030 \exp(iLn\dt) &\approx& \bigg(\exp\left(iL_1 \dt\right) \exp\left(iL_2 \Dt \right) \bigg)^{n}.
1031 \eea
1032 Although at first this does not appear symmetric, as long as the full velocity
1033 step is between $t-\hDt$ and $t+\hDt$, then this is simply a way of
1034 starting velocity Verlet at a different place in the cycle.
1036 Even though the trajectory and thus potential energies are identical
1037 between leap-frog and velocity Verlet, the kinetic energy and
1038 temperature will not necessarily be the same. Standard velocity
1039 Verlet uses the velocities at the $t$ to calculate the kinetic energy
1040 and thus the temperature only at time $t$; the kinetic energy is then a sum over all particles
1041 \bea
1042 KE_{\mathrm{full}}(t) &=& \sum_i \left(\frac{1}{2m_i}\ve{v}_i(t)\right)^2 \nonumber\\
1043 &=& \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,
1044 \eea
1045 with the square on the {\em outside} of the average. Standard
1046 leap-frog calculates the kinetic energy at time $t$ based on the
1047 average kinetic energies at the timesteps $t+\hDt$ and $t-\hDt$, or
1048 the sum over all particles
1049 \bea
1050 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),
1051 \eea
1052 where the square is {\em inside} the average.
1054 A non-standard variant of velocity Verlet which averages the kinetic
1055 energies $KE(t+\hDt)$ and $KE(t-\hDt)$, exactly like leap-frog, is also
1056 now implemented in {\gromacs} (as {\tt .mdp} file option {\tt md-vv-avek}). Without
1057 temperature and pressure coupling, velocity Verlet with
1058 half-step-averaged kinetic energies and leap-frog will be identical up
1059 to numerical precision. For temperature- and pressure-control schemes,
1060 however, velocity Verlet with half-step-averaged kinetic energies and
1061 leap-frog will be different, as will be discussed in the section in
1062 thermostats and barostats.
1064 The half-step-averaged kinetic energy and temperature are slightly more
1065 accurate for a given step size; the difference in average kinetic
1066 energies using the half-step-averaged kinetic energies ({\em md} and
1067 {\em md-vv-avek}) will be closer to the kinetic energy obtained in the
1068 limit of small step size than will the full-step kinetic energy (using
1069 {\em md-vv}). For NVE simulations, this difference is usually not
1070 significant, since the positions and velocities of the particles are
1071 still identical; it makes a difference in the way the the temperature
1072 of the simulations are {\em interpreted}, but {\em not} in the
1073 trajectories that are produced. Although the kinetic energy is more
1074 accurate with the half-step-averaged method, meaning that it changes
1075 less as the timestep gets large, it is also more noisy. The RMS deviation
1076 of the total energy of the system (sum of kinetic plus
1077 potential) in the half-step-averaged kinetic energy case will be
1078 higher (about twice as high in most cases) than the full-step kinetic
1079 energy. The drift will still be the same, however, as again, the
1080 trajectories are identical.
1082 For NVT simulations, however, there {\em will} be a difference, as
1083 discussed in the section on temperature control, since the velocities
1084 of the particles are adjusted such that kinetic energies of the
1085 simulations, which can be calculated either way, reach the
1086 distribution corresponding to the set temperature. In this case, the
1087 three methods will not give identical results.
1089 Because the velocity and position are both defined at the same time
1090 $t$ the velocity Verlet integrator can be used for some methods,
1091 especially rigorously correct pressure control methods, that are not
1092 actually possible with leap-frog. The integration itself takes
1093 negligibly more time than leap-frog, but twice as many communication
1094 calls are currently required. In most cases, and especially for large
1095 systems where communication speed is important for parallelization and
1096 differences between thermodynamic ensembles vanish in the $1/N$ limit,
1097 and when only NVT ensembles are required, leap-frog will likely be the
1098 preferred integrator. For pressure control simulations where the fine
1099 details of the thermodynamics are important, only velocity Verlet
1100 allows the true ensemble to be calculated. In either case, simulation
1101 with double precision may be required to get fine details of
1102 thermodynamics correct.
1104 \subsection{Twin-range cut-offs\index{twin-range!cut-off}}
1105 To save computation time, slowly varying forces can be calculated
1106 less often than rapidly varying forces. In {\gromacs}
1107 such a \normindex{multiple time step} splitting is possible between
1108 short and long range non-bonded interactions.
1109 In {\gromacs} versions up to 4.0, an irreversible integration scheme
1110 was used which is also used by the {\gromos} simulation package:
1111 every $n$ steps the long range forces are determined and these are
1112 then also used (without modification) for the next $n-1$ integration steps
1113 in \eqnref{leapfrogv}. Such an irreversible scheme can result in bad energy
1114 conservation and, possibly, bad sampling.
1115 Since version 4.5, a leap-frog version of the reversible Trotter decomposition scheme~\cite{Tuckerman1992a} is used.
1116 In this integrator the long-range forces are determined every $n$ steps
1117 and are then integrated into the velocity in \eqnref{leapfrogv} using
1118 a time step of $\Dt_\mathrm{LR} = n \Dt$:
1119 \beq
1120 \ve{v}(t+\hDt) =
1121 \left\{ \begin{array}{lll} \displaystyle
1122 \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
1123 \ve{v}(t-\hDt) + \frac{1}{m}\ve{F}_\mathrm{SR}(t)\Dt &,& \mathrm{step} ~\%~ n \neq 0 \\
1124 \end{array} \right.
1125 \eeq
1127 The parameter $n$ is equal to the neighbor list update frequency. In
1128 4.5, the velocity Verlet version of multiple time-stepping is not yet
1129 fully implemented.
1131 Several other simulation packages uses multiple time stepping for
1132 bonds and/or the PME mesh forces. In {\gromacs} we have not implemented
1133 this (yet), since we use a different philosophy. Bonds can be constrained
1134 (which is also a more sound approximation of a physical quantum
1135 oscillator), which allows the smallest time step to be increased
1136 to the larger one. This not only halves the number of force calculations,
1137 but also the update calculations. For even larger time steps, angle vibrations
1138 involving hydrogen atoms can be removed using virtual interaction
1139 \ifthenelse{\equal{\gmxlite}{1}}
1140 {sites,}
1141 {sites (see \secref{rmfast}),}
1142 which brings the shortest time step up to
1143 PME mesh update frequency of a multiple time stepping scheme.
1145 As an example we show the energy conservation for integrating
1146 the equations of motion for SPC/E water at 300 K. To avoid cut-off
1147 effects, reaction-field electrostatics with $\epsilon_{RF}=\infty$ and
1148 shifted Lennard-Jones interactions are used, both with a buffer region.
1149 The long-range interactions were evaluated between 1.0 and 1.4 nm.
1150 In \figref{leapfrog} one can see that for electrostatics the Trotter scheme
1151 does an order of magnitude better up to $\Dt_{LR}$ = 16 fs.
1152 The electrostatics depends strongly on the orientation of the water molecules,
1153 which changes rapidly.
1154 For Lennard-Jones interactions, the energy drift is linear in $\Dt_{LR}$
1155 and roughly two orders of magnitude smaller than for the electrostatics.
1156 Lennard-Jones forces are smaller than Coulomb forces and
1157 they are mainly affected by translation of water molecules, not rotation.
1159 \begin{figure}
1160 \centerline{\includegraphics[width=12cm]{plots/drift-all}}
1161 \caption{Energy drift per degree of freedom in SPC/E water
1162 with twin-range cut-offs
1163 for reaction field (left) and Lennard-Jones interaction (right)
1164 as a function of the long-range time step length for the irreversible
1165 ``\gromos'' scheme and a reversible Trotter scheme.}
1166 \label{fig:twinrangeener}
1167 \end{figure}
1169 \subsection{Temperature coupling\index{temperature coupling}}
1170 While direct use of molecular dynamics gives rise to the NVE (constant
1171 number, constant volume, constant energy ensemble), most quantities
1172 that we wish to calculate are actually from a constant temperature
1173 (NVT) ensemble, also called the canonical ensemble. {\gromacs} can use
1174 the {\em weak-coupling} scheme of Berendsen~\cite{Berendsen84},
1175 stochastic randomization through the Andersen
1176 thermostat~\cite{Andersen80}, the extended ensemble Nos{\'e}-Hoover
1177 scheme~\cite{Nose84,Hoover85}, or a velocity-rescaling
1178 scheme~\cite{Bussi2007a} to simulate constant temperature, with
1179 advantages of each of the schemes laid out below.
1181 There are several other reasons why it might be necessary to control
1182 the temperature of the system (drift during equilibration, drift as a
1183 result of force truncation and integration errors, heating due to
1184 external or frictional forces), but this is not entirely correct to do
1185 from a thermodynamic standpoint, and in some cases only masks the
1186 symptoms (increase in temperature of the system) rather than the
1187 underlying problem (deviations from correct physics in the dynamics).
1188 For larger systems, errors in ensemble averages and structural
1189 properties incurred by using temperature control to remove slow drifts
1190 in temperature appear to be negligible, but no completely
1191 comprehensive comparisons have been carried out, and some caution must
1192 be taking in interpreting the results.
1194 \subsubsection{Berendsen temperature coupling\pawsindexquiet{Berendsen}{temperature coupling}\index{weak coupling}}
1195 The Berendsen algorithm mimics weak coupling with first-order
1196 kinetics to an external heat bath with given temperature $T_0$.
1197 See ref.~\cite{Berendsen91} for a comparison with the
1198 Nos{\'e}-Hoover scheme. The effect of this algorithm is
1199 that a deviation of the system temperature from $T_0$ is slowly
1200 corrected according to:
1201 \beq
1202 \frac{\de T}{\de t} = \frac{T_0-T}{\tau}
1203 \label{eqn:Tcoupling}
1204 \eeq
1205 which means that a temperature deviation decays exponentially with a
1206 time constant $\tau$.
1207 This method of coupling has the advantage that the strength of the
1208 coupling can be varied and adapted to the user requirement: for
1209 equilibration purposes the coupling time can be taken quite short
1210 ({\eg} 0.01 ps), but for reliable equilibrium runs it can be taken much
1211 longer ({\eg} 0.5 ps) in which case it hardly influences the
1212 conservative dynamics.
1214 The Berendsen thermostat suppresses the fluctuations of the kinetic
1215 energy. This means that one does not generate a proper canonical
1216 ensemble, so rigorously, the sampling will be incorrect. This
1217 error scales with $1/N$, so for very large systems most ensemble
1218 averages will not be affected significantly, except for the
1219 distribution of the kinetic energy itself. However, fluctuation
1220 properties, such as the heat capacity, will be affected. A similar
1221 thermostat which does produce a correct ensemble is the velocity
1222 rescaling thermostat~\cite{Bussi2007a} described below.
1224 The heat flow into or out of the system is affected by scaling the
1225 velocities of each particle every step, or every $n_\mathrm{TC}$ steps,
1226 with a time-dependent factor $\lambda$, given by:
1227 \beq
1228 \lambda = \left[ 1 + \frac{n_\mathrm{TC} \Delta t}{\tau_T}
1229 \left\{\frac{T_0}{T(t - \hDt)} - 1 \right\} \right]^{1/2}
1230 \label{eqn:lambda}
1231 \eeq
1232 The parameter $\tau_T$ is close, but not exactly equal, to the time constant
1233 $\tau$ of the temperature coupling (\eqnref{Tcoupling}):
1234 \beq
1235 \tau = 2 C_V \tau_T / N_{df} k
1236 \eeq
1237 where $C_V$ is the total heat capacity of the system, $k$ is Boltzmann's
1238 constant, and $N_{df}$ is the total number of degrees of freedom. The
1239 reason that $\tau \neq \tau_T$ is that the kinetic energy change
1240 caused by scaling the velocities is partly redistributed between
1241 kinetic and potential energy and hence the change in temperature is
1242 less than the scaling energy. In practice, the ratio $\tau / \tau_T$
1243 ranges from 1 (gas) to 2 (harmonic solid) to 3 (water). When we use
1244 the term ``temperature coupling time constant,'' we mean the parameter
1245 \normindex{$\tau_T$}.
1246 {\bf Note} that in practice the scaling factor $\lambda$ is limited to
1247 the range of 0.8 $<= \lambda <=$ 1.25, to avoid scaling by very large
1248 numbers which may crash the simulation. In normal use,
1249 $\lambda$ will always be much closer to 1.0.
1251 \subsubsection{Velocity-rescaling temperature coupling\pawsindexquiet{velocity-rescaling}{temperature coupling}}
1252 The velocity-rescaling thermostat~\cite{Bussi2007a} is essentially a Berendsen
1253 thermostat (see above) with an additional stochastic term that ensures
1254 a correct kinetic energy distribution by modifying it according to
1255 \beq
1256 \de K = (K_0 - K) \frac{\de t}{\tau_T} + 2 \sqrt{\frac{K K_0}{N_f}} \frac{\de W}{\sqrt{\tau_T}},
1257 \label{eqn:vrescale}
1258 \eeq
1259 where $K$ is the kinetic energy, $N_f$ the number of degrees of freedom and $\de W$ a Wiener process.
1260 There are no additional parameters, except for a random seed.
1261 This thermostat produces a correct canonical ensemble and still has
1262 the advantage of the Berendsen thermostat: first order decay of
1263 temperature deviations and no oscillations.
1264 When an $NVT$ ensemble is used, the conserved energy quantity
1265 is written to the energy and log file.
1267 \subsubsection{\normindex{Andersen thermostat}}
1268 One simple way to maintain a thermostatted ensemble is to take an
1269 $NVE$ integrator and periodically re-select the velocities of the
1270 particles from a Maxwell-Boltzmann distribution.~\cite{Andersen80}
1271 This can either be done by randomizing all the velocities
1272 simultaneously (massive collision) every $\tau_T/\Dt$ steps ({\tt andersen-massive}), or by
1273 randomizing every particle with some small probability every timestep ({\tt andersen}),
1274 equal to $\Dt/\tau$, where in both cases $\Dt$ is the timestep and
1275 $\tau_T$ is a characteristic coupling time scale.
1276 Because of the way constraints operate, all particles in the same
1277 constraint group must be randomized simultaneously. Because of
1278 parallelization issues, the {\tt andersen} version cannot currently (5.0) be
1279 used in systems with constraints. {\tt andersen-massive} can be used regardless of constraints.
1280 This thermostat is also currently only possible with velocity Verlet algorithms,
1281 because it operates directly on the velocities at each timestep.
1283 This algorithm completely avoids some of the ergodicity issues of other thermostatting
1284 algorithms, as energy cannot flow back and forth between energetically
1285 decoupled components of the system as in velocity scaling motions.
1286 However, it can slow down the kinetics of system by randomizing
1287 correlated motions of the system, including slowing sampling when
1288 $\tau_T$ is at moderate levels (less than 10 ps). This algorithm
1289 should therefore generally not be used when examining kinetics or
1290 transport properties of the system.~\cite{Basconi2013}
1292 % \ifthenelse{\equal{\gmxlite}{1}}{}{
1293 \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}}
1295 The Berendsen weak-coupling algorithm is
1296 extremely efficient for relaxing a system to the target temperature,
1297 but once the system has reached equilibrium it might be more
1298 important to probe a correct canonical ensemble. This is unfortunately
1299 not the case for the weak-coupling scheme.
1301 To enable canonical ensemble simulations, {\gromacs} also supports the
1302 extended-ensemble approach first proposed by Nos{\'e}~\cite{Nose84}
1303 and later modified by Hoover~\cite{Hoover85}. The system Hamiltonian is
1304 extended by introducing a thermal reservoir and a friction term in the
1305 equations of motion. The friction force is proportional to the
1306 product of each particle's velocity and a friction parameter, $\xi$.
1307 This friction parameter (or ``heat bath'' variable) is a fully
1308 dynamic quantity with its own momentum ($p_{\xi}$) and equation of
1309 motion; the time derivative is calculated from the difference between
1310 the current kinetic energy and the reference temperature.
1312 In this formulation, the particles' equations of motion in
1313 \figref{global} are replaced by:
1314 \beq
1315 \frac {\de^2\ve{r}_i}{\de t^2} = \frac{\ve{F}_i}{m_i} -
1316 \frac{p_{\xi}}{Q}\frac{\de \ve{r}_i}{\de t} ,
1317 \label{eqn:NH-eqn-of-motion}
1318 \eeq where the equation of motion for the heat bath parameter $\xi$ is:
1319 \beq \frac {\de p_{\xi}}{\de t} = \left( T - T_0 \right). \eeq The
1320 reference temperature is denoted $T_0$, while $T$ is the current
1321 instantaneous temperature of the system. The strength of the coupling
1322 is determined by the constant $Q$ (usually called the ``mass parameter''
1323 of the reservoir) in combination with the reference
1324 temperature.~\footnote{Note that some derivations, an alternative
1325 notation $\xi_{\mathrm{alt}} = v_{\xi} = p_{\xi}/Q$ is used.}
1327 The conserved quantity for the Nos{\'e}-Hoover equations of motion is not
1328 the total energy, but rather
1329 \bea
1330 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,
1331 \eea
1332 where $N_f$ is the total number of degrees of freedom.
1334 In our opinion, the mass parameter is a somewhat awkward way of
1335 describing coupling strength, especially due to its dependence on
1336 reference temperature (and some implementations even include the
1337 number of degrees of freedom in your system when defining $Q$). To
1338 maintain the coupling strength, one would have to change $Q$ in
1339 proportion to the change in reference temperature. For this reason, we
1340 prefer to let the {\gromacs} user work instead with the period
1341 $\tau_T$ of the oscillations of kinetic energy between the system and
1342 the reservoir instead. It is directly related to $Q$ and $T_0$ via:
1343 \beq
1344 Q = \frac {\tau_T^2 T_0}{4 \pi^2}.
1345 \eeq
1346 This provides a much more intuitive way of selecting the
1347 Nos{\'e}-Hoover coupling strength (similar to the weak-coupling
1348 relaxation), and in addition $\tau_T$ is independent of system size
1349 and reference temperature.
1351 It is however important to keep the difference between the
1352 weak-coupling scheme and the Nos{\'e}-Hoover algorithm in mind:
1353 Using weak coupling you get a
1354 strongly damped {\em exponential relaxation},
1355 while the Nos{\'e}-Hoover approach
1356 produces an {\em oscillatory relaxation}.
1357 The actual time it takes to relax with Nos{\'e}-Hoover coupling is
1358 several times larger than the period of the
1359 oscillations that you select. These oscillations (in contrast
1360 to exponential relaxation) also means that
1361 the time constant normally should be 4--5 times larger
1362 than the relaxation time used with weak coupling, but your
1363 mileage may vary.
1365 Nos{\'e}-Hoover dynamics in simple systems such as collections of
1366 harmonic oscillators, can be {\em nonergodic}, meaning that only a
1367 subsection of phase space is ever sampled, even if the simulations
1368 were to run for infinitely long. For this reason, the Nos{\'e}-Hoover
1369 chain approach was developed, where each of the Nos{\'e}-Hoover
1370 thermostats has its own Nos{\'e}-Hoover thermostat controlling its
1371 temperature. In the limit of an infinite chain of thermostats, the
1372 dynamics are guaranteed to be ergodic. Using just a few chains can
1373 greatly improve the ergodicity, but recent research has shown that the
1374 system will still be nonergodic, and it is still not entirely clear
1375 what the practical effect of this~\cite{Cooke2008}. Currently, the
1376 default number of chains is 10, but this can be controlled by the
1377 user. In the case of chains, the equations are modified in the
1378 following way to include a chain of thermostatting
1379 particles~\cite{Martyna1992}:
1381 \bea
1382 \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 \\
1383 \frac {\de p_{{\xi}_1}}{\de t} &~=~& \left( T - T_0 \right) - p_{{\xi}_1} \frac{p_{{\xi}_2}}{Q_2} \nonumber \\
1384 \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 \\
1385 \frac {\de p_{\xi_N}}{\de t} &~=~& \left(\frac{p_{\xi_{N-1}}^2}{Q_{N-1}}-kT\right)
1386 \label{eqn:NH-chain-eqn-of-motion}
1387 \eea
1388 The conserved quantity for Nos{\'e}-Hoover chains is
1389 \bea
1390 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
1391 \eea
1392 The values and velocities of the Nos{\'e}-Hoover thermostat variables
1393 are generally not included in the output, as they take up a fair
1394 amount of space and are generally not important for analysis of
1395 simulations, but this can be overridden by defining the environment
1396 variable {\tt GMX_NOSEHOOVER_CHAINS}, which will print the values of all
1397 the positions and velocities of all Nos{\'e}-Hoover particles in the
1398 chain to the {\tt .edr} file. Leap-frog simulations currently can only have
1399 Nos{\'e}-Hoover chain lengths of 1, but this will likely be updated in
1400 later version.
1402 As described in the integrator section, for temperature coupling, the
1403 temperature that the algorithm attempts to match to the reference
1404 temperature is calculated differently in velocity Verlet and leap-frog
1405 dynamics. Velocity Verlet ({\em md-vv}) uses the full-step kinetic
1406 energy, while leap-frog and {\em md-vv-avek} use the half-step-averaged
1407 kinetic energy.
1409 We can examine the Trotter decomposition again to better understand
1410 the differences between these constant-temperature integrators. In
1411 the case of Nos{\'e}-Hoover dynamics (for simplicity, using a chain
1412 with $N=1$, with more details in Ref.~\cite{Martyna1996}), we split
1413 the Liouville operator as
1414 \beq
1415 iL = iL_1 + iL_2 + iL_{\mathrm{NHC}},
1416 \eeq
1417 where
1418 \bea
1419 iL_1 &=& \sum_{i=1}^N \left[\frac{\pb_i}{m_i}\right]\cdot \frac{\partial}{\partial \rv_i} \nonumber \\
1420 iL_2 &=& \sum_{i=1}^N \F_i\cdot \frac{\partial}{\partial \pb_i} \nonumber \\
1421 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}}
1422 \eea
1423 For standard velocity Verlet with Nos{\'e}-Hoover temperature control, this becomes
1424 \bea
1425 \exp(iL\dt) &=& \exp\left(iL_{\mathrm{NHC}}\dt/2\right) \exp\left(iL_2 \dt/2\right) \nonumber \\
1426 &&\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).
1427 \eea
1428 For half-step-averaged temperature control using {\em md-vv-avek},
1429 this decomposition will not work, since we do not have the full step
1430 temperature until after the second velocity step. However, we can
1431 construct an alternate decomposition that is still reversible, by
1432 switching the place of the NHC and velocity portions of the
1433 decomposition:
1434 \bea
1435 \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 \\
1436 &&\exp\left(iL_{\mathrm{NHC}}\dt/2\right) \exp\left(iL_2 \dt/2\right)+ \mathcal{O}(\Dt^3)
1437 \label{eq:half_step_NHC_integrator}
1438 \eea
1439 This formalism allows us to easily see the difference between the
1440 different flavors of velocity Verlet integrator. The leap-frog
1441 integrator can be seen as starting with
1442 Eq.~\ref{eq:half_step_NHC_integrator} just before the $\exp\left(iL_1
1443 \dt\right)$ term, yielding:
1444 \bea
1445 \exp(iL\dt) &=& \exp\left(iL_1 \dt\right) \exp\left(iL_{\mathrm{NHC}}\dt/2\right) \nonumber \\
1446 &&\exp\left(iL_2 \dt\right) \exp\left(iL_{\mathrm{NHC}}\dt/2\right) + \mathcal{O}(\Dt^3)
1447 \eea
1448 and then using some algebra tricks to solve for some quantities are
1449 required before they are actually calculated~\cite{Holian95}.
1453 \subsubsection{Group temperature coupling}\index{temperature-coupling group}%
1454 In {\gromacs} temperature coupling can be performed on groups of
1455 atoms, typically a protein and solvent. The reason such algorithms
1456 were introduced is that energy exchange between different components
1457 is not perfect, due to different effects including cut-offs etc. If
1458 now the whole system is coupled to one heat bath, water (which
1459 experiences the largest cut-off noise) will tend to heat up and the
1460 protein will cool down. Typically 100 K differences can be obtained.
1461 With the use of proper electrostatic methods (PME) these difference
1462 are much smaller but still not negligible. The parameters for
1463 temperature coupling in groups are given in the {\tt mdp} file.
1464 Recent investigation has shown that small temperature differences
1465 between protein and water may actually be an artifact of the way
1466 temperature is calculated when there are finite timesteps, and very
1467 large differences in temperature are likely a sign of something else
1468 seriously going wrong with the system, and should be investigated
1469 carefully~\cite{Eastwood2010}.
1471 One special case should be mentioned: it is possible to temperature-couple only
1472 part of the system, leaving other parts without temperature
1473 coupling. This is done by specifying ${-1}$ for the time constant
1474 $\tau_T$ for the group that should not be thermostatted. If only
1475 part of the system is thermostatted, the system will still eventually
1476 converge to an NVT system. In fact, one suggestion for minimizing
1477 errors in the temperature caused by discretized timesteps is that if
1478 constraints on the water are used, then only the water degrees of
1479 freedom should be thermostatted, not protein degrees of freedom, as
1480 the higher frequency modes in the protein can cause larger deviations
1481 from the ``true'' temperature, the temperature obtained with small
1482 timesteps~\cite{Eastwood2010}.
1484 \subsection{Pressure coupling\index{pressure coupling}}
1485 In the same spirit as the temperature coupling, the system can also be
1486 coupled to a ``pressure bath.'' {\gromacs} supports both the Berendsen
1487 algorithm~\cite{Berendsen84} that scales coordinates and box vectors
1488 every step, the extended-ensemble Parrinello-Rahman approach~\cite{Parrinello81,Nose83}, and for
1489 the velocity Verlet variants, the Martyna-Tuckerman-Tobias-Klein
1490 (MTTK) implementation of pressure
1491 control~\cite{Martyna1996}. Parrinello-Rahman and Berendsen can be
1492 combined with any of the temperature coupling methods above. MTTK can
1493 only be used with Nos{\'e}-Hoover temperature control. From 5.1 afterwards,
1494 it can only used when the system does not have constraints.
1496 \subsubsection{Berendsen pressure coupling\pawsindexquiet{Berendsen}{pressure coupling}\index{weak coupling}}
1497 \label{sec:berendsen_pressure_coupling}
1498 The Berendsen algorithm rescales the
1499 coordinates and box vectors every step, or every $n_\mathrm{PC}$ steps,
1500 with a matrix {\boldmath $\mu$},
1501 which has the effect of a first-order kinetic relaxation of the pressure
1502 towards a given reference pressure ${\bf P}_0$ according to
1503 \beq
1504 \frac{\de {\bf P}}{\de t} = \frac{{\bf P}_0-{\bf P}}{\tau_p}.
1505 \eeq
1506 The scaling matrix {\boldmath $\mu$} is given by
1507 \beq
1508 \mu_{ij}
1509 = \delta_{ij} - \frac{n_\mathrm{PC}\Delta t}{3\, \tau_p} \beta_{ij} \{P_{0ij} - P_{ij}(t) \}.
1510 \label{eqn:mu}
1511 \eeq
1512 \index{isothermal compressibility}
1513 \index{compressibility}
1514 Here, {\boldmath $\beta$} is the isothermal compressibility of the system.
1515 In most cases this will be a diagonal matrix, with equal elements on the
1516 diagonal, the value of which is generally not known.
1517 It suffices to take a rough estimate because the value of {\boldmath $\beta$}
1518 only influences the non-critical time constant of the
1519 pressure relaxation without affecting the average pressure itself.
1520 For water at 1 atm and 300 K
1521 $\beta = 4.6 \times 10^{-10}$ Pa$^{-1} = 4.6 \times 10^{-5}$ bar$^{-1}$,
1522 which is $7.6 \times 10^{-4}$ MD units (see \chref{defunits}).
1523 Most other liquids have similar values.
1524 When scaling completely anisotropically, the system has to be rotated in
1525 order to obey \eqnref{box_rot}.
1526 This rotation is approximated in first order in the scaling, which is usually
1527 less than $10^{-4}$. The actual scaling matrix {\boldmath $\mu'$} is
1528 \beq
1529 \mbox{\boldmath $\mu'$} =
1530 \left(\begin{array}{ccc}
1531 \mu_{xx} & \mu_{xy} + \mu_{yx} & \mu_{xz} + \mu_{zx} \\
1532 0 & \mu_{yy} & \mu_{yz} + \mu_{zy} \\
1533 0 & 0 & \mu_{zz}
1534 \end{array}\right).
1535 \eeq
1536 The velocities are neither scaled nor rotated.
1538 In {\gromacs}, the Berendsen scaling can also be done isotropically,
1539 which means that instead of $\ve{P}$ a diagonal matrix with elements of size
1540 trace$(\ve{P})/3$ is used. For systems with interfaces, semi-isotropic
1541 scaling can be useful.
1542 In this case, the $x/y$-directions are scaled isotropically and the $z$
1543 direction is scaled independently. The compressibility in the $x/y$ or
1544 $z$-direction can be set to zero, to scale only in the other direction(s).
1546 If you allow full anisotropic deformations and use constraints you
1547 might have to scale more slowly or decrease your timestep to avoid
1548 errors from the constraint algorithms. It is important to note that
1549 although the Berendsen pressure control algorithm yields a simulation
1550 with the correct average pressure, it does not yield the exact NPT
1551 ensemble, and it is not yet clear exactly what errors this approximation
1552 may yield.
1554 % \ifthenelse{\equal{\gmxlite}{1}}{}{
1555 \subsubsection{Parrinello-Rahman pressure coupling\pawsindexquiet{Parrinello-Rahman}{pressure coupling}}
1557 In cases where the fluctuations in pressure or volume are important
1558 {\em per se} ({\eg} to calculate thermodynamic properties), especially
1559 for small systems, it may be a problem that the exact ensemble is not
1560 well defined for the weak-coupling scheme, and that it does not
1561 simulate the true NPT ensemble.
1563 {\gromacs} also supports constant-pressure simulations using the
1564 Parrinello-Rahman approach~\cite{Parrinello81,Nose83}, which is similar
1565 to the Nos{\'e}-Hoover temperature coupling, and in theory gives the
1566 true NPT ensemble. With the Parrinello-Rahman barostat, the box
1567 vectors as represented by the matrix \ve{b} obey the matrix equation
1568 of motion\footnote{The box matrix representation \ve{b} in {\gromacs}
1569 corresponds to the transpose of the box matrix representation \ve{h}
1570 in the paper by Nos{\'e} and Klein. Because of this, some of our
1571 equations will look slightly different.}
1572 \beq
1573 \frac{\de \ve{b}^2}{\de t^2}= V \ve{W}^{-1} \ve{b}'^{-1} \left( \ve{P} - \ve{P}_{ref}\right).
1574 \eeq
1576 The volume of the box is denoted $V$, and $\ve{W}$ is a matrix parameter that determines
1577 the strength of the coupling. The matrices \ve{P} and \ve{P}$_{ref}$ are the
1578 current and reference pressures, respectively.
1580 The equations of motion for the particles are also changed, just as
1581 for the Nos{\'e}-Hoover coupling. In most cases you would combine the
1582 Parrinello-Rahman barostat with the Nos{\'e}-Hoover
1583 thermostat, but to keep it simple we only show the Parrinello-Rahman
1584 modification here:
1586 \bea \frac {\de^2\ve{r}_i}{\de t^2} & = & \frac{\ve{F}_i}{m_i} -
1587 \ve{M} \frac{\de \ve{r}_i}{\de t} , \\ \ve{M} & = & \ve{b}^{-1} \left[
1588 \ve{b} \frac{\de \ve{b}'}{\de t} + \frac{\de \ve{b}}{\de t} \ve{b}'
1589 \right] \ve{b}'^{-1}. \eea The (inverse) mass parameter matrix
1590 $\ve{W}^{-1}$ determines the strength of the coupling, and how the box
1591 can be deformed. The box restriction (\ref{eqn:box_rot}) will be
1592 fulfilled automatically if the corresponding elements of $\ve{W}^{-1}$
1593 are zero. Since the coupling strength also depends on the size of your
1594 box, we prefer to calculate it automatically in {\gromacs}. You only
1595 have to provide the approximate isothermal compressibilities
1596 {\boldmath $\beta$} and the pressure time constant $\tau_p$ in the
1597 input file ($L$ is the largest box matrix element): \beq \left(
1598 \ve{W}^{-1} \right)_{ij} = \frac{4 \pi^2 \beta_{ij}}{3 \tau_p^2 L}.
1599 \eeq Just as for the Nos{\'e}-Hoover thermostat, you should realize
1600 that the Parrinello-Rahman time constant is {\em not} equivalent to
1601 the relaxation time used in the Berendsen pressure coupling algorithm.
1602 In most cases you will need to use a 4--5 times larger time constant
1603 with Parrinello-Rahman coupling. If your pressure is very far from
1604 equilibrium, the Parrinello-Rahman coupling may result in very large
1605 box oscillations that could even crash your run. In that case you
1606 would have to increase the time constant, or (better) use the weak-coupling
1607 scheme to reach the target pressure, and then switch to
1608 Parrinello-Rahman coupling once the system is in equilibrium.
1609 Additionally, using the leap-frog algorithm, the pressure at time $t$
1610 is not available until after the time step has completed, and so the
1611 pressure from the previous step must be used, which makes the algorithm
1612 not directly reversible, and may not be appropriate for high precision
1613 thermodynamic calculations.
1615 \subsubsection{Surface-tension coupling\pawsindexquiet{surface-tension}{pressure coupling}}
1616 When a periodic system consists of more than one phase, separated by
1617 surfaces which are parallel to the $xy$-plane,
1618 the surface tension and the $z$-component of the pressure can be coupled
1619 to a pressure bath. Presently, this only works with the Berendsen
1620 pressure coupling algorithm in {\gromacs}.
1621 The average surface tension $\gamma(t)$ can be calculated from
1622 the difference between the normal and the lateral pressure
1623 \bea
1624 \gamma(t) & = &
1625 \frac{1}{n} \int_0^{L_z}
1626 \left\{ P_{zz}(z,t) - \frac{P_{xx}(z,t) + P_{yy}(z,t)}{2} \right\} \mbox{d}z \\
1627 & = &
1628 \frac{L_z}{n} \left\{ P_{zz}(t) - \frac{P_{xx}(t) + P_{yy}(t)}{2} \right\},
1629 \eea
1630 where $L_z$ is the height of the box and $n$ is the number of surfaces.
1631 The pressure in the z-direction is corrected by scaling the height of
1632 the box with $\mu_{zz}$
1633 \beq
1634 \Delta P_{zz} = \frac{\Delta t}{\tau_p} \{ P_{0zz} - P_{zz}(t) \}
1635 \eeq
1636 \beq
1637 \mu_{zz} = 1 + \beta_{zz} \Delta P_{zz}
1638 \eeq
1639 This is similar to normal pressure coupling, except that the factor
1640 of $1/3$ is missing.
1641 The pressure correction in the $z$-direction is then used to get the
1642 correct convergence for the surface tension to the reference value $\gamma_0$.
1643 The correction factor for the box length in the $x$/$y$-direction is
1644 \beq
1645 \mu_{x/y} = 1 + \frac{\Delta t}{2\,\tau_p} \beta_{x/y}
1646 \left( \frac{n \gamma_0}{\mu_{zz} L_z}
1647 - \left\{ P_{zz}(t)+\Delta P_{zz} - \frac{P_{xx}(t) + P_{yy}(t)}{2} \right\}
1648 \right)
1649 \eeq
1650 The value of $\beta_{zz}$ is more critical than with normal pressure
1651 coupling. Normally an incorrect compressibility will just scale $\tau_p$,
1652 but with surface tension coupling it affects the convergence of the surface
1653 tension.
1654 When $\beta_{zz}$ is set to zero (constant box height), $\Delta P_{zz}$ is also set
1655 to zero, which is necessary for obtaining the correct surface tension.
1657 \subsubsection{MTTK pressure control algorithms}
1659 As mentioned in the previous section, one weakness of leap-frog
1660 integration is in constant pressure simulations, since the pressure
1661 requires a calculation of both the virial and the kinetic energy at
1662 the full time step; for leap-frog, this information is not available
1663 until {\em after} the full timestep. Velocity Verlet does allow the
1664 calculation, at the cost of an extra round of global communication,
1665 and can compute, mod any integration errors, the true NPT ensemble.
1667 The full equations, combining both pressure coupling and temperature
1668 coupling, are taken from Martyna {\em et al.}~\cite{Martyna1996} and
1669 Tuckerman~\cite{Tuckerman2006} and are referred to here as MTTK
1670 equations (Martyna-Tuckerman-Tobias-Klein). We introduce for
1671 convenience $\epsilon = (1/3)\ln (V/V_0)$, where $V_0$ is a reference
1672 volume. The momentum of $\epsilon$ is $\veps = p_{\epsilon}/W =
1673 \dot{\epsilon} = \dot{V}/3V$, and define $\alpha = 1 + 3/N_{dof}$ (see
1674 Ref~\cite{Tuckerman2006})
1676 The isobaric equations are
1677 \bea
1678 \dot{\rv}_i &=& \frac{\pb_i}{m_i} + \frac{\peps}{W} \rv_i \nonumber \\
1679 \frac{\dot{\pb}_i}{m_i} &=& \frac{1}{m_i}\F_i - \alpha\frac{\peps}{W} \frac{\pb_i}{m_i} \nonumber \\
1680 \dot{\epsilon} &=& \frac{\peps}{W} \nonumber \\
1681 \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),\\
1682 \eea
1683 where
1684 \bea
1685 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\
1686 \right)\right].
1687 \eea
1688 The terms including $\alpha$ are required to make phase space
1689 incompressible~\cite{Tuckerman2006}. The $\epsilon$ acceleration term
1690 can be rewritten as
1691 \bea
1692 \frac{\dot{\peps}}{W} &=& \frac{3V}{W}\left(\alpha P_{\mathrm{kin}} - P_{\mathrm{vir}} - P\right)
1693 \eea
1694 In terms of velocities, these equations become
1695 \bea
1696 \dot{\rv}_i &=& \vv_i + \veps \rv_i \nonumber \\
1697 \dot{\vv}_i &=& \frac{1}{m_i}\F_i - \alpha\veps \vv_i \nonumber \\
1698 \dot{\epsilon} &=& \veps \nonumber \\
1699 \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 \\
1700 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]
1701 \eea
1702 For these equations, the conserved quantity is
1703 \bea
1704 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
1705 \eea
1706 The next step is to add temperature control. Adding Nos{\'e}-Hoover
1707 chains, including to the barostat degree of freedom, where we use
1708 $\eta$ for the barostat Nos{\'e}-Hoover variables, and $Q^{\prime}$
1709 for the coupling constants of the thermostats of the barostats, we get
1710 \bea
1711 \dot{\rv}_i &=& \frac{\pb_i}{m_i} + \frac{\peps}{W} \rv_i \nonumber \\
1712 \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 \\
1713 \dot{\epsilon} &=& \frac{\peps}{W} \nonumber \\
1714 \frac{\dot{\peps}}{W} &=& \frac{3V}{W}(\alpha P_{\mathrm{kin}} - P_{\mathrm{vir}} - P) -\frac{p_{\eta_1}}{Q^{\prime}_1}\peps \nonumber \\
1715 \dot{\xi}_k &=& \frac{p_{\xi_k}}{Q_k} \nonumber \\
1716 \dot{\eta}_k &=& \frac{p_{\eta_k}}{Q^{\prime}_k} \nonumber \\
1717 \dot{p}_{\xi_k} &=& G_k - \frac{p_{\xi_{k+1}}}{Q_{k+1}} \;\;\;\; k=1,\ldots, M-1 \nonumber \\
1718 \dot{p}_{\eta_k} &=& G^\prime_k - \frac{p_{\eta_{k+1}}}{Q^\prime_{k+1}} \;\;\;\; k=1,\ldots, M-1 \nonumber \\
1719 \dot{p}_{\xi_M} &=& G_M \nonumber \\
1720 \dot{p}_{\eta_M} &=& G^\prime_M, \nonumber \\
1721 \eea
1722 where
1723 \bea
1724 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 \\
1725 G_1 &=& \sum_{i=1}^N \frac{\pb^2_i}{m_i} - N_f kT \nonumber \\
1726 G_k &=& \frac{p^2_{\xi_{k-1}}}{2Q_{k-1}} - kT \;\; k = 2,\ldots,M \nonumber \\
1727 G^\prime_1 &=& \frac{\peps^2}{2W} - kT \nonumber \\
1728 G^\prime_k &=& \frac{p^2_{\eta_{k-1}}}{2Q^\prime_{k-1}} - kT \;\; k = 2,\ldots,M
1729 \eea
1730 The conserved quantity is now
1731 \bea
1732 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 \\
1733 \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
1734 \eea
1735 Returning to the Trotter decomposition formalism, for pressure control and temperature control~\cite{Martyna1996} we get:
1736 \bea
1737 iL = iL_1 + iL_2 + iL_{\epsilon,1} + iL_{\epsilon,2} + iL_{\mathrm{NHC-baro}} + iL_{\mathrm{NHC}}
1738 \eea
1739 where ``NHC-baro'' corresponds to the Nos{\`e}-Hoover chain of the barostat,
1740 and NHC corresponds to the NHC of the particles,
1741 \bea
1742 iL_1 &=& \sum_{i=1}^N \left[\frac{\pb_i}{m_i} + \frac{\peps}{W}\rv_i\right]\cdot \frac{\partial}{\partial \rv_i} \\
1743 iL_2 &=& \sum_{i=1}^N \F_i - \alpha \frac{\peps}{W}\pb_i \cdot \frac{\partial}{\partial \pb_i} \\
1744 iL_{\epsilon,1} &=& \frac{p_\epsilon}{W} \frac{\partial}{\partial \epsilon}\\
1745 iL_{\epsilon,2} &=& G_{\epsilon} \frac{\partial}{\partial p_\epsilon}
1746 \eea
1747 and where
1748 \bea
1749 G_{\epsilon} = 3V\left(\alpha P_{\mathrm{kin}} - P_{\mathrm{vir}} - P\right)
1750 \eea
1751 Using the Trotter decomposition, we get
1752 \bea
1753 \exp(iL\dt) &=& \exp\left(iL_{\mathrm{NHC-baro}}\dt/2\right)\exp\left(iL_{\mathrm{NHC}}\dt/2\right) \nonumber \nonumber \\
1754 &&\exp\left(iL_{\epsilon,2}\dt/2\right) \exp\left(iL_2 \dt/2\right) \nonumber \nonumber \\
1755 &&\exp\left(iL_{\epsilon,1}\dt\right) \exp\left(iL_1 \dt\right) \nonumber \nonumber \\
1756 &&\exp\left(iL_2 \dt/2\right) \exp\left(iL_{\epsilon,2}\dt/2\right) \nonumber \nonumber \\
1757 &&\exp\left(iL_{\mathrm{NHC}}\dt/2\right)\exp\left(iL_{\mathrm{NHC-baro}}\dt/2\right) + \mathcal{O}(\dt^3)
1758 \eea
1759 The action of $\exp\left(iL_1 \dt\right)$ comes from the solution of
1760 the the differential equation
1761 $\dot{\rv}_i = \vv_i + \veps \rv_i$
1762 with $\vv_i = \pb_i/m_i$ and $\veps$ constant with initial condition
1763 $\rv_i(0)$, evaluate at $t=\Delta t$. This yields the evolution
1764 \beq
1765 \rv_i(\dt) = \rv_i(0)e^{\veps \dt} + \Delta t \vv_i(0) e^{\veps \dt/2} \sinhx{\veps \dt/2}.
1766 \eeq
1767 The action of $\exp\left(iL_2 \dt/2\right)$ comes from the solution
1768 of the differential equation $\dot{\vv}_i = \frac{\F_i}{m_i} -
1769 \alpha\veps\vv_i$, yielding
1770 \beq
1771 \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}.
1772 \eeq
1773 {\em md-vv-avek} uses the full step kinetic energies for determining the pressure with the pressure control,
1774 but the half-step-averaged kinetic energy for the temperatures, which can be written as a Trotter decomposition as
1775 \bea
1776 \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 \\
1777 &&\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 \\
1778 &&\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)
1779 \eea
1781 With constraints, the equations become significantly more complicated,
1782 in that each of these equations need to be solved iteratively for the
1783 constraint forces. Before {\gromacs} 5.1, these iterative
1784 constraints were solved as described in~\cite{Yu2010}. From {\gromacs}
1785 5.1 onward, MTTK with constraints has been removed because of
1786 numerical stability issues with the iterations.
1788 \subsubsection{Infrequent evaluation of temperature and pressure coupling}
1790 Temperature and pressure control require global communication to
1791 compute the kinetic energy and virial, which can become costly if
1792 performed every step for large systems. We can rearrange the Trotter
1793 decomposition to give alternate symplectic, reversible integrator with
1794 the coupling steps every $n$ steps instead of every steps. These new
1795 integrators will diverge if the coupling time step is too large, as
1796 the auxiliary variable integrations will not converge. However, in
1797 most cases, long coupling times are more appropriate, as they disturb
1798 the dynamics less~\cite{Martyna1996}.
1800 Standard velocity Verlet with Nos{\'e}-Hoover temperature control has a Trotter expansion
1801 \bea
1802 \exp(iL\dt) &\approx& \exp\left(iL_{\mathrm{NHC}}\dt/2\right) \exp\left(iL_2 \dt/2\right) \nonumber \\
1803 &&\exp\left(iL_1 \dt\right) \exp\left(iL_2 \dt/2\right) \exp\left(iL_{\mathrm{NHC}}\dt/2\right).
1804 \eea
1805 If the Nos{\'e}-Hoover chain is sufficiently slow with respect to the motions of the system, we can
1806 write an alternate integrator over $n$ steps for velocity Verlet as
1807 \bea
1808 \exp(iL\dt) &\approx& (\exp\left(iL_{\mathrm{NHC}}(n\dt/2)\right)\left[\exp\left(iL_2 \dt/2\right)\right. \nonumber \\
1809 &&\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).
1810 \eea
1811 For pressure control, this becomes
1812 \bea
1813 \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 \\
1814 &&\exp\left(iL_{\epsilon,2}(n\dt/2)\right) \left[\exp\left(iL_2 \dt/2\right)\right. \nonumber \nonumber \\
1815 &&\exp\left(iL_{\epsilon,1}\dt\right) \exp\left(iL_1 \dt\right) \nonumber \nonumber \\
1816 &&\left.\exp\left(iL_2 \dt/2\right)\right]^n \exp\left(iL_{\epsilon,2}(n\dt/2)\right) \nonumber \nonumber \\
1817 &&\exp\left(iL_{\mathrm{NHC}}(n\dt/2)\right)\exp\left(iL_{\mathrm{NHC-baro}}(n\dt/2)\right),
1818 \eea
1819 where the box volume integration occurs every step, but the auxiliary variable
1820 integrations happen every $n$ steps.
1822 % } % Brace matches ifthenelse test for gmxlite
1825 \subsection{The complete update algorithm}
1826 \begin{figure}
1827 \begin{center}
1828 \addtolength{\fboxsep}{0.5cm}
1829 \begin{shadowenv}[12cm]
1830 {\large \bf THE UPDATE ALGORITHM}
1831 \rule{\textwidth}{2pt} \\
1832 Given:\\
1833 Positions $\ve{r}$ of all atoms at time $t$ \\
1834 Velocities $\ve{v}$ of all atoms at time $t-\hDt$ \\
1835 Accelerations $\ve{F}/m$ on all atoms at time $t$.\\
1836 (Forces are computed disregarding any constraints)\\
1837 Total kinetic energy and virial at $t-\Dt$\\
1838 $\Downarrow$ \\
1839 {\bf 1.} Compute the scaling factors $\lambda$ and $\mu$\\
1840 according to \eqnsref{lambda}{mu}\\
1841 $\Downarrow$ \\
1842 {\bf 2.} Update and scale velocities: $\ve{v}' = \lambda (\ve{v} +
1843 \ve{a} \Delta t)$ \\
1844 $\Downarrow$ \\
1845 {\bf 3.} Compute new unconstrained coordinates: $\ve{r}' = \ve{r} + \ve{v}'
1846 \Delta t$ \\
1847 $\Downarrow$ \\
1848 {\bf 4.} Apply constraint algorithm to coordinates: constrain($\ve{r}^{'} \rightarrow \ve{r}'';
1849 \, \ve{r}$) \\
1850 $\Downarrow$ \\
1851 {\bf 5.} Correct velocities for constraints: $\ve{v} = (\ve{r}'' -
1852 \ve{r}) / \Delta t$ \\
1853 $\Downarrow$ \\
1854 {\bf 6.} Scale coordinates and box: $\ve{r} = \mu \ve{r}''; \ve{b} =
1855 \mu \ve{b}$ \\
1856 \end{shadowenv}
1857 \caption{The MD update algorithm with the leap-frog integrator}
1858 \label{fig:complete-update}
1859 \end{center}
1860 \end{figure}
1861 The complete algorithm for the update of velocities and coordinates is
1862 given using leap-frog in \figref{complete-update}. The SHAKE algorithm of step
1863 4 is explained below.
1865 {\gromacs} has a provision to ``freeze'' (prevent motion of) selected
1866 particles\index{frozen atoms}, which must be defined as a ``\swapindex{freeze}{group}.'' This is implemented
1867 using a {\em freeze factor $\ve{f}_g$}, which is a vector, and differs for each
1868 freeze group (see \secref{groupconcept}). This vector contains only
1869 zero (freeze) or one (don't freeze).
1870 When we take this freeze factor and the external acceleration $\ve{a}_h$ into
1871 account the update algorithm for the velocities becomes
1872 \beq
1873 \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],
1874 \eeq
1875 where $g$ and $h$ are group indices which differ per atom.
1877 \subsection{Output step}
1878 The most important output of the MD run is the {\em
1879 \swapindex{trajectory}{file}}, which contains particle coordinates
1880 and (optionally) velocities at regular intervals.
1881 The trajectory file contains frames that could include positions,
1882 velocities and/or forces, as well as information about the dimensions
1883 of the simulation volume, integration step, integration time, etc. The
1884 interpretation of the time varies with the integrator chosen, as
1885 described above. For Velocity Verlet integrators, velocities labeled
1886 at time $t$ are for that time. For other integrators (e.g. leap-frog,
1887 stochastic dynamics), the velocities labeled at time $t$ are for time
1888 $t - \hDt$.
1890 Since the trajectory
1891 files are lengthy, one should not save every step! To retain all
1892 information it suffices to write a frame every 15 steps, since at
1893 least 30 steps are made per period of the highest frequency in the
1894 system, and Shannon's \normindex{sampling} theorem states that two samples per
1895 period of the highest frequency in a band-limited signal contain all
1896 available information. But that still gives very long files! So, if
1897 the highest frequencies are not of interest, 10 or 20 samples per ps
1898 may suffice. Be aware of the distortion of high-frequency motions by
1899 the {\em stroboscopic effect}, called {\em aliasing}: higher frequencies
1900 are mirrored with respect to the sampling frequency and appear as
1901 lower frequencies.
1903 {\gromacs} can also write reduced-precision coordinates for a subset of
1904 the simulation system to a special compressed trajectory file
1905 format. All the other tools can read and write this format. See
1906 the User Guide for details on how to set up your {\tt .mdp} file
1907 to have {\tt mdrun} use this feature.
1909 % \ifthenelse{\equal{\gmxlite}{1}}{}{
1910 \section{Shell molecular dynamics}
1911 {\gromacs} can simulate \normindex{polarizability} using the
1912 \normindex{shell model} of Dick and Overhauser~\cite{Dick58}. In such models
1913 a shell particle representing the electronic degrees of freedom is
1914 attached to a nucleus by a spring. The potential energy is minimized with
1915 respect to the shell position at every step of the simulation (see below).
1916 Successful applications of shell models in {\gromacs} have been published
1917 for $N_2$~\cite{Jordan95} and water~\cite{Maaren2001a}.
1919 \subsection{Optimization of the shell positions}
1920 The force \ve{F}$_S$ on a shell particle $S$ can be decomposed into two
1921 components
1922 \begin{equation}
1923 \ve{F}_S ~=~ \ve{F}_{bond} + \ve{F}_{nb}
1924 \end{equation}
1925 where \ve{F}$_{bond}$ denotes the component representing the
1926 polarization energy, usually represented by a harmonic potential and
1927 \ve{F}$_{nb}$ is the sum of Coulomb and van der Waals interactions. If we
1928 assume that \ve{F}$_{nb}$ is almost constant we can analytically derive the
1929 optimal position of the shell, i.e. where \ve{F}$_S$ = 0. If we have the
1930 shell S connected to atom A we have
1931 \begin{equation}
1932 \ve{F}_{bond} ~=~ k_b \left( \ve{x}_S - \ve{x}_A\right).
1933 \end{equation}
1934 In an iterative solver, we have positions \ve{x}$_S(n)$ where $n$ is
1935 the iteration count. We now have at iteration $n$
1936 \begin{equation}
1937 \ve{F}_{nb} ~=~ \ve{F}_S - k_b \left( \ve{x}_S(n) - \ve{x}_A\right)
1938 \end{equation}
1939 and the optimal position for the shells $x_S(n+1)$ thus follows from
1940 \begin{equation}
1941 \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
1942 \end{equation}
1943 if we write
1944 \begin{equation}
1945 \Delta \ve{x}_S = \ve{x}_S(n+1) - \ve{x}_S(n)
1946 \end{equation}
1947 we finally obtain
1948 \begin{equation}
1949 \Delta \ve{x}_S = \ve{F}_S/k_b
1950 \end{equation}
1951 which then yields the algorithm to compute the next trial in the optimization
1952 of shell positions
1953 \begin{equation}
1954 \ve{x}_S(n+1) ~=~ \ve{x}_S(n) + \ve{F}_S/k_b.
1955 \end{equation}
1956 % } % Brace matches ifthenelse test for gmxlite
1958 \section{Constraint algorithms\index{constraint algorithms}}
1959 Constraints can be imposed in {\gromacs} using LINCS (default) or
1960 the traditional SHAKE method.
1962 \subsection{\normindex{SHAKE}}
1963 \label{subsec:SHAKE}
1964 The SHAKE~\cite{Ryckaert77} algorithm changes a set of unconstrained
1965 coordinates $\ve{r}^{'}$ to a set of coordinates $\ve{r}''$ that
1966 fulfill a list of distance constraints, using a set $\ve{r}$
1967 reference, as
1968 \beq
1969 {\rm SHAKE}(\ve{r}^{'} \rightarrow \ve{r}'';\, \ve{r})
1970 \eeq
1971 This action is consistent with solving a set of Lagrange multipliers
1972 in the constrained equations of motion. SHAKE needs a {\em relative tolerance};
1973 it will continue until all constraints are satisfied within
1974 that relative tolerance. An error message is
1975 given if SHAKE cannot reset the coordinates because the deviation is
1976 too large, or if a given number of iterations is surpassed.
1978 Assume the equations of motion must fulfill $K$ holonomic constraints,
1979 expressed as
1980 \beq
1981 \sigma_k(\ve{r}_1 \ldots \ve{r}_N) = 0; \;\; k=1 \ldots K.
1982 \eeq
1983 For example, $(\ve{r}_1 - \ve{r}_2)^2 - b^2 = 0$.
1984 Then the forces are defined as
1985 \beq
1986 - \frac{\partial}{\partial \ve{r}_i} \left( V + \sum_{k=1}^K \lambda_k
1987 \sigma_k \right),
1988 \eeq
1989 where $\lambda_k$ are Lagrange multipliers which must be solved to
1990 fulfill the constraint equations. The second part of this sum
1991 determines the {\em constraint forces} $\ve{G}_i$, defined by
1992 \beq
1993 \ve{G}_i = -\sum_{k=1}^K \lambda_k \frac{\partial \sigma_k}{\partial
1994 \ve{r}_i}
1995 \eeq
1996 The displacement due to the constraint forces in the leap-frog or
1997 Verlet algorithm is equal to $(\ve{G}_i/m_i)(\Dt)^2$. Solving the
1998 Lagrange multipliers (and hence the displacements) requires the
1999 solution of a set of coupled equations of the second degree. These are
2000 solved iteratively by SHAKE.
2001 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2002 \label{subsec:SETTLE}
2003 For the special case of rigid water molecules, that often make up more
2004 than 80\% of the simulation system we have implemented the
2005 \normindex{SETTLE}
2006 algorithm~\cite{Miyamoto92} (\secref{constraints}).
2008 For velocity Verlet, an additional round of constraining must be
2009 done, to constrain the velocities of the second velocity half step,
2010 removing any component of the velocity parallel to the bond vector.
2011 This step is called RATTLE, and is covered in more detail in the
2012 original Andersen paper~\cite{Andersen1983a}.
2014 % } % Brace matches ifthenelse test for gmxlite
2019 \newcommand{\fs}[1]{\begin{equation} \label{eqn:#1}}
2020 \newcommand{\fe}{\end{equation}}
2021 \newcommand{\p}{\partial}
2022 \newcommand{\Bm}{\ve{B}}
2023 \newcommand{\M}{\ve{M}}
2024 \newcommand{\iM}{\M^{-1}}
2025 \newcommand{\Tm}{\ve{T}}
2026 \newcommand{\Sm}{\ve{S}}
2027 \newcommand{\fo}{\ve{f}}
2028 \newcommand{\con}{\ve{g}}
2029 \newcommand{\lenc}{\ve{d}}
2031 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2032 \subsection{\normindex{LINCS}}
2033 \label{subsec:lincs}
2035 \subsubsection{The LINCS algorithm}
2036 LINCS is an algorithm that resets bonds to their correct lengths
2037 after an unconstrained update~\cite{Hess97}.
2038 The method is non-iterative, as it always uses two steps.
2039 Although LINCS is based on matrices, no matrix-matrix multiplications are
2040 needed. The method is more stable and faster than SHAKE,
2041 but it can only be used with bond constraints and
2042 isolated angle constraints, such as the proton angle in OH.
2043 Because of its stability, LINCS is especially useful for Brownian dynamics.
2044 LINCS has two parameters, which are explained in the subsection parameters.
2045 The parallel version of LINCS, P-LINCS, is described
2046 in subsection \ssecref{plincs}.
2048 \subsubsection{The LINCS formulas}
2049 We consider a system of $N$ particles, with positions given by a
2050 $3N$ vector $\ve{r}(t)$.
2051 For molecular dynamics the equations of motion are given by Newton's Law
2052 \fs{c1}
2053 {\de^2 \ve{r} \over \de t^2} = \iM \ve{F},
2055 where $\ve{F}$ is the $3N$ force vector
2056 and $\M$ is a $3N \times 3N$ diagonal matrix,
2057 containing the masses of the particles.
2058 The system is constrained by $K$ time-independent constraint equations
2059 \fs{c2}
2060 g_i(\ve{r}) = | \ve{r}_{i_1}-\ve{r}_{i_2} | - d_i = 0 ~~~~~~i=1,\ldots,K.
2063 In a numerical integration scheme, LINCS is applied after an
2064 unconstrained update, just like SHAKE. The algorithm works in two
2065 steps (see figure \figref{lincs}). In the first step, the projections
2066 of the new bonds on the old bonds are set to zero. In the second step,
2067 a correction is applied for the lengthening of the bonds due to
2068 rotation. The numerics for the first step and the second step are very
2069 similar. A complete derivation of the algorithm can be found in
2070 \cite{Hess97}. Only a short description of the first step is given
2071 here.
2073 \begin{figure}
2074 \centerline{\includegraphics[height=50mm]{plots/lincs}}
2075 \caption[The three position updates needed for one time step.]{The
2076 three position updates needed for one time step. The dashed line is
2077 the old bond of length $d$, the solid lines are the new bonds. $l=d
2078 \cos \theta$ and $p=(2 d^2 - l^2)^{1 \over 2}$.}
2079 \label{fig:lincs}
2080 \end{figure}
2082 A new notation is introduced for the gradient matrix of the constraint
2083 equations which appears on the right hand side of this equation:
2084 \fs{c3}
2085 B_{hi} = {\p g_h \over \p r_i}
2087 Notice that $\Bm$ is a $K \times 3N$ matrix, it contains the directions
2088 of the constraints.
2089 The following equation shows how the new constrained coordinates
2090 $\ve{r}_{n+1}$ are related to the unconstrained coordinates
2091 $\ve{r}_{n+1}^{unc}$ by
2092 \fs{m0}
2093 \begin{array}{c}
2094 \ve{r}_{n+1}=(\ve{I}-\Tm_n \ve{B}_n) \ve{r}_{n+1}^{unc} + \Tm_n \lenc=
2095 \\[2mm]
2096 \ve{r}_{n+1}^{unc} -
2097 \iM \Bm_n (\Bm_n \iM \Bm_n^T)^{-1} (\Bm_n \ve{r}_{n+1}^{unc} - \lenc)
2098 \end{array}
2100 where $\Tm = \iM \Bm^T (\Bm \iM \Bm^T)^{-1}$.
2101 The derivation of this equation from \eqnsref{c1}{c2} can be found
2102 in \cite{Hess97}.
2104 This first step does not set the real bond lengths to the prescribed lengths,
2105 but the projection of the new bonds onto the old directions of the bonds.
2106 To correct for the rotation of bond $i$, the projection of the
2107 bond, $p_i$, on the old direction is set to
2108 \fs{m1a}
2109 p_i=\sqrt{2 d_i^2 - l_i^2},
2111 where $l_i$ is the bond length after the first projection.
2112 The corrected positions are
2113 \fs{m1b}
2114 \ve{r}_{n+1}^*=(\ve{I}-\Tm_n \Bm_n)\ve{r}_{n+1} + \Tm_n \ve{p}.
2116 This correction for rotational effects is actually an iterative process,
2117 but during MD only one iteration is applied.
2118 The relative constraint deviation after this procedure will be less than
2119 0.0001 for every constraint.
2120 In energy minimization, this might not be accurate enough, so the number
2121 of iterations is equal to the order of the expansion (see below).
2123 Half of the CPU time goes to inverting the constraint coupling
2124 matrix $\Bm_n \iM \Bm_n^T$, which has to be done every time step.
2125 This $K \times K$ matrix
2126 has $1/m_{i_1} + 1/m_{i_2}$ on the diagonal.
2127 The off-diagonal elements are only non-zero when two bonds are connected,
2128 then the element is
2129 $\cos \phi /m_c$, where $m_c$ is
2130 the mass of the atom connecting the
2131 two bonds and $\phi$ is the angle between the bonds.
2133 The matrix $\Tm$ is inverted through a power expansion.
2134 A $K \times K$ matrix $\ve{S}$ is
2135 introduced which is the inverse square root of
2136 the diagonal of $\Bm_n \iM \Bm_n^T$.
2137 This matrix is used to convert the diagonal elements
2138 of the coupling matrix to one:
2139 \fs{m2}
2140 \begin{array}{c}
2141 (\Bm_n \iM \Bm_n^T)^{-1}
2142 = \Sm \Sm^{-1} (\Bm_n \iM \Bm_n^T)^{-1} \Sm^{-1} \Sm \\[2mm]
2143 = \Sm (\Sm \Bm_n \iM \Bm_n^T \Sm)^{-1} \Sm =
2144 \Sm (\ve{I} - \ve{A}_n)^{-1} \Sm
2145 \end{array}
2147 The matrix $\ve{A}_n$ is symmetric and sparse and has zeros on the diagonal.
2148 Thus a simple trick can be used to calculate the inverse:
2149 \fs{m3}
2150 (\ve{I}-\ve{A}_n)^{-1}=
2151 \ve{I} + \ve{A}_n + \ve{A}_n^2 + \ve{A}_n^3 + \ldots
2154 This inversion method is only valid if the absolute values of all the
2155 eigenvalues of $\ve{A}_n$ are smaller than one.
2156 In molecules with only bond constraints, the connectivity is so low
2157 that this will always be true, even if ring structures are present.
2158 Problems can arise in angle-constrained molecules.
2159 By constraining angles with additional distance constraints,
2160 multiple small ring structures are introduced.
2161 This gives a high connectivity, leading to large eigenvalues.
2162 Therefore LINCS should NOT be used with coupled angle-constraints.
2164 For molecules with all bonds constrained the eigenvalues of $A$
2165 are around 0.4. This means that with each additional order
2166 in the expansion \eqnref{m3} the deviations decrease by a factor 0.4.
2167 But for relatively isolated triangles of constraints the largest
2168 eigenvalue is around 0.7.
2169 Such triangles can occur when removing hydrogen angle vibrations
2170 with an additional angle constraint in alcohol groups
2171 or when constraining water molecules with LINCS, for instance
2172 with flexible constraints.
2173 The constraints in such triangles converge twice as slow as
2174 the other constraints. Therefore, starting with {\gromacs} 4,
2175 additional terms are added to the expansion for such triangles
2176 \fs{m3_ang}
2177 (\ve{I}-\ve{A}_n)^{-1} \approx
2178 \ve{I} + \ve{A}_n + \ldots + \ve{A}_n^{N_i} +
2179 \left(\ve{A}^*_n + \ldots + {\ve{A}_n^*}^{N_i} \right) \ve{A}_n^{N_i}
2181 where $N_i$ is the normal order of the expansion and
2182 $\ve{A}^*$ only contains the elements of $\ve{A}$ that couple
2183 constraints within rigid triangles, all other elements are zero.
2184 In this manner, the accuracy of angle constraints comes close
2185 to that of the other constraints, while the series of matrix vector
2186 multiplications required for determining the expansion
2187 only needs to be extended for a few constraint couplings.
2188 This procedure is described in the P-LINCS paper\cite{Hess2008a}.
2190 \subsubsection{The LINCS Parameters}
2191 The accuracy of LINCS depends on the number of matrices used
2192 in the expansion \eqnref{m3}. For MD calculations a fourth order
2193 expansion is enough. For Brownian dynamics with
2194 large time steps an eighth order expansion may be necessary.
2195 The order is a parameter in the {\tt *.mdp} file.
2196 The implementation of LINCS is done in such a way that the
2197 algorithm will never crash. Even when it is impossible to
2198 to reset the constraints LINCS will generate a conformation
2199 which fulfills the constraints as well as possible.
2200 However, LINCS will generate a warning when in one step a bond
2201 rotates over more than a predefined angle.
2202 This angle is set by the user in the {\tt *.mdp} file.
2204 % } % Brace matches ifthenelse test for gmxlite
2207 \section{Simulated Annealing}
2208 \label{sec:SA}
2209 The well known \swapindex{simulated}{annealing}
2210 (SA) protocol is supported in {\gromacs}, and you can even couple multiple
2211 groups of atoms separately with an arbitrary number of reference temperatures
2212 that change during the simulation. The annealing is implemented by simply
2213 changing the current reference temperature for each group in the temperature
2214 coupling, so the actual relaxation and coupling properties depends on the
2215 type of thermostat you use and how hard you are coupling it. Since we are
2216 changing the reference temperature it is important to remember that the system
2217 will NOT instantaneously reach this value - you need to allow for the inherent
2218 relaxation time in the coupling algorithm too. If you are changing the
2219 annealing reference temperature faster than the temperature relaxation you
2220 will probably end up with a crash when the difference becomes too large.
2222 The annealing protocol is specified as a series of corresponding times and
2223 reference temperatures for each group, and you can also choose whether you only
2224 want a single sequence (after which the temperature will be coupled to the
2225 last reference value), or if the annealing should be periodic and restart at
2226 the first reference point once the sequence is completed. You can mix and
2227 match both types of annealing and non-annealed groups in your simulation.
2229 \newcommand{\vrond}{\stackrel{\circ}{\ve{r}}}
2230 \newcommand{\rond}{\stackrel{\circ}{r}}
2231 \newcommand{\ruis}{\ve{r}^G}
2233 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2234 \section{Stochastic Dynamics\swapindexquiet{stochastic}{dynamics}}
2235 \label{sec:SD}
2236 Stochastic or velocity \swapindex{Langevin}{dynamics} adds a friction
2237 and a noise term to Newton's equations of motion, as
2238 \beq
2239 \label{SDeq}
2240 m_i {\de^2 \ve{r}_i \over \de t^2} =
2241 - m_i \gamma_i {\de \ve{r}_i \over \de t} + \ve{F}_i(\ve{r}) + \vrond_i,
2242 \eeq
2243 where $\gamma_i$ is the friction constant $[1/\mbox{ps}]$ and
2244 $\vrond_i\!\!(t)$ is a noise process with
2245 $\langle \rond_i\!\!(t) \rond_j\!\!(t+s) \rangle =
2246 2 m_i \gamma_i k_B T \delta(s) \delta_{ij}$.
2247 When $1/\gamma_i$ is large compared to the time scales present in the system,
2248 one could see stochastic dynamics as molecular dynamics with stochastic
2249 temperature-coupling. The advantage compared to MD with Berendsen
2250 temperature-coupling is that in case of SD the generated ensemble is known.
2251 For simulating a system in vacuum there is the additional advantage that there is no
2252 accumulation of errors for the overall translational and rotational
2253 degrees of freedom.
2254 When $1/\gamma_i$ is small compared to the time scales present in the system,
2255 the dynamics will be completely different from MD, but the sampling is
2256 still correct.
2258 In {\gromacs} there are two algorithms to integrate equation (\ref{SDeq}):
2259 a simple and efficient one
2260 and a more complex leap-frog algorithm~\cite{Gunsteren88}, which is now deprecated.
2261 The accuracy of both integrators is equivalent to the normal MD leap-frog and
2262 Velocity Verlet integrator, except with constraints where the complex
2263 SD integrator samples at a temperature that is slightly too high (although that error is smaller than the one from the Velocity Verlet integrator that uses the kinetic energy from the full-step velocity). The simple integrator is nearly identical to the common way of discretizing the Langevin equation, but the friction and velocity term are applied in an impulse fashion~\cite{Goga2012}.
2264 The simple integrator is:
2265 \bea
2266 \label{eqn:sd_int1}
2267 \ve{v}' &~=~& \ve{v}(t-\hDt) + \frac{1}{m}\ve{F}(t)\Dt \\
2268 \Delta\ve{v} &~=~& -\alpha \, \ve{v}'(t+\hDt) + \sqrt{\frac{k_B T}{m}(1 - \alpha^2)} \, \ruis_i \\
2269 \ve{r}(t+\Dt) &~=~& \ve{r}(t)+\left(\ve{v}' +\frac{1}{2}\Delta \ve{v}\right)\Dt \label{eqn:sd1_x_upd}\\
2270 \ve{v}(t+\hDt) &~=~& \ve{v}' + \Delta \ve{v} \\
2271 \alpha &~=~& 1 - e^{-\gamma \Dt}
2272 \eea
2273 where $\ruis_i$ is Gaussian distributed noise with $\mu = 0$, $\sigma = 1$.
2274 The velocity is first updated a full time step without friction and noise to get $\ve{v}'$, identical to the normal update in leap-frog. The friction and noise are then applied as an impulse at step $t+\Dt$. The advantage of this scheme is that the velocity-dependent terms act at the full time step, which makes the correct integration of forces that depend on both coordinates and velocities, such as constraints and dissipative particle dynamics (DPD, not implented yet), straightforward. With constraints, the coordinate update \eqnref{sd1_x_upd} is split into a normal leap-frog update and a $\Delta \ve{v}$. After both of these updates the constraints are applied to coordinates and velocities.
2276 In the deprecated complex algorithm, four Gaussian random numbers are required
2277 per integration step per degree of freedom, and with constraints the
2278 coordinates need to be constrained twice per integration step.
2279 Depending on the computational cost of the force calculation,
2280 this can take a significant part of the simulation time.
2281 Exact continuation of a stochastic dynamics simulation is not possible,
2282 because the state of the random number generator is not stored.
2284 When using SD as a thermostat, an appropriate value for $\gamma$ is e.g. 0.5 ps$^{-1}$,
2285 since this results in a friction that is lower than the internal friction
2286 of water, while it still provides efficient thermostatting.
2289 \section{Brownian Dynamics\swapindexquiet{Brownian}{dynamics}}
2290 \label{sec:BD}
2291 In the limit of high friction, stochastic dynamics reduces to
2292 Brownian dynamics, also called position Langevin dynamics.
2293 This applies to over-damped systems,
2294 {\ie} systems in which the inertia effects are negligible.
2295 The equation is
2296 \beq
2297 {\de \ve{r}_i \over \de t} = \frac{1}{\gamma_i} \ve{F}_i(\ve{r}) + \vrond_i
2298 \eeq
2299 where $\gamma_i$ is the friction coefficient $[\mbox{amu/ps}]$ and
2300 $\vrond_i\!\!(t)$ is a noise process with
2301 $\langle \rond_i\!\!(t) \rond_j\!\!(t+s) \rangle =
2302 2 \delta(s) \delta_{ij} k_B T / \gamma_i$.
2303 In {\gromacs} the equations are integrated with a simple, explicit scheme
2304 \beq
2305 \ve{r}_i(t+\Delta t) = \ve{r}_i(t) +
2306 {\Delta t \over \gamma_i} \ve{F}_i(\ve{r}(t))
2307 + \sqrt{2 k_B T {\Delta t \over \gamma_i}}\, \ruis_i,
2308 \eeq
2309 where $\ruis_i$ is Gaussian distributed noise with $\mu = 0$, $\sigma = 1$.
2310 The friction coefficients $\gamma_i$ can be chosen the same for all
2311 particles or as $\gamma_i = m_i\,\gamma_i$, where the friction constants
2312 $\gamma_i$ can be different for different groups of atoms.
2313 Because the system is assumed to be over-damped, large timesteps
2314 can be used. LINCS should be used for the constraints since SHAKE
2315 will not converge for large atomic displacements.
2316 BD is an option of the {\tt mdrun} program.
2317 % } % Brace matches ifthenelse test for gmxlite
2319 \section{Energy Minimization}
2320 \label{sec:EM}\index{energy minimization}%
2321 Energy minimization in {\gromacs} can be done using steepest descent,
2322 conjugate gradients, or l-bfgs (limited-memory
2323 Broyden-Fletcher-Goldfarb-Shanno quasi-Newtonian minimizer...we
2324 prefer the abbreviation). EM is just an option of the {\tt mdrun}
2325 program.
2327 \subsection{Steepest Descent\index{steepest descent}}
2328 Although steepest descent is certainly not the most efficient
2329 algorithm for searching, it is robust and easy to implement.
2331 We define the vector $\ve{r}$ as the vector of all $3N$ coordinates.
2332 Initially a maximum displacement $h_0$ ({\eg} 0.01 nm) must be given.
2334 First the forces $\ve{F}$ and potential energy are calculated.
2335 New positions are calculated by
2336 \beq
2337 \ve{r}_{n+1} = \ve{r}_n + \frac{\ve{F}_n}{\max (|\ve{F}_n|)} h_n,
2338 \eeq
2339 where $h_n$ is the maximum displacement and $\ve{F}_n$ is the force,
2340 or the negative gradient of the potential $V$. The notation $\max
2341 (|\ve{F}_n|)$ means the largest of the absolute values of the force
2342 components. The forces and energy are again computed for the new positions \\
2343 If ($V_{n+1} < V_n$) the new positions are accepted and $h_{n+1} = 1.2
2344 h_n$. \\
2345 If ($V_{n+1} \geq V_n$) the new positions are rejected and $h_n = 0.2 h_n$.
2347 The algorithm stops when either a user-specified number of force
2348 evaluations has been performed ({\eg} 100), or when the maximum of the absolute
2349 values of the force (gradient) components is smaller than a specified
2350 value $\epsilon$.
2351 Since force truncation produces some noise in the
2352 energy evaluation, the stopping criterion should not be made too tight
2353 to avoid endless iterations. A reasonable value for $\epsilon$ can be
2354 estimated from the root mean square force $f$ a harmonic oscillator would exhibit at a
2355 temperature $T$. This value is
2356 \beq
2357 f = 2 \pi \nu \sqrt{ 2mkT},
2358 \eeq
2359 where $\nu$ is the oscillator frequency, $m$ the (reduced) mass, and
2360 $k$ Boltzmann's constant. For a weak oscillator with a wave number of
2361 100 cm$^{-1}$ and a mass of 10 atomic units, at a temperature of 1 K,
2362 $f=7.7$ kJ~mol$^{-1}$~nm$^{-1}$. A value for $\epsilon$ between 1 and
2363 10 is acceptable.
2365 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2366 \subsection{Conjugate Gradient\index{conjugate gradient}}
2367 Conjugate gradient is slower than steepest descent in the early stages
2368 of the minimization, but becomes more efficient closer to the energy
2369 minimum. The parameters and stop criterion are the same as for
2370 steepest descent. In {\gromacs} conjugate gradient can not be used
2371 with constraints, including the SETTLE algorithm for
2372 water~\cite{Miyamoto92}, as this has not been implemented. If water is
2373 present it must be of a flexible model, which can be specified in the
2374 {\tt *.mdp} file by {\tt define = -DFLEXIBLE}.
2376 This is not really a restriction, since the accuracy of conjugate
2377 gradient is only required for minimization prior to a normal-mode
2378 analysis, which cannot be performed with constraints. For most other
2379 purposes steepest descent is efficient enough.
2380 % } % Brace matches ifthenelse test for gmxlite
2382 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2383 \subsection{\normindex{L-BFGS}}
2384 The original BFGS algorithm works by successively creating better
2385 approximations of the inverse Hessian matrix, and moving the system to
2386 the currently estimated minimum. The memory requirements for this are
2387 proportional to the square of the number of particles, so it is not
2388 practical for large systems like biomolecules. Instead, we use the
2389 L-BFGS algorithm of Nocedal~\cite{Byrd95a,Zhu97a}, which approximates
2390 the inverse Hessian by a fixed number of corrections from previous
2391 steps. This sliding-window technique is almost as efficient as the
2392 original method, but the memory requirements are much lower -
2393 proportional to the number of particles multiplied with the correction
2394 steps. In practice we have found it to converge faster than conjugate
2395 gradients, but due to the correction steps it is not yet parallelized.
2396 It is also noteworthy that switched or shifted interactions usually
2397 improve the convergence, since sharp cut-offs mean the potential
2398 function at the current coordinates is slightly different from the
2399 previous steps used to build the inverse Hessian approximation.
2400 % } % Brace matches ifthenelse test for gmxlite
2402 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2403 \section{Normal-Mode Analysis\index{normal-mode analysis}\index{NMA}}
2404 Normal-mode analysis~\cite{Levitt83,Go83,BBrooks83b}
2405 can be performed using {\gromacs}, by diagonalization of the mass-weighted
2406 \normindex{Hessian} $H$:
2407 \bea
2408 R^T M^{-1/2} H M^{-1/2} R &=& \mbox{diag}(\lambda_1,\ldots,\lambda_{3N})
2410 \lambda_i &=& (2 \pi \omega_i)^2
2411 \eea
2412 where $M$ contains the atomic masses, $R$ is a matrix that contains
2413 the eigenvectors as columns, $\lambda_i$ are the eigenvalues
2414 and $\omega_i$ are the corresponding frequencies.
2416 First the Hessian matrix, which is a $3N \times 3N$ matrix where $N$
2417 is the number of atoms, needs to be calculated:
2418 \bea
2419 H_{ij} &=& \frac{\partial^2 V}{\partial x_i \partial x_j}
2420 \eea
2421 where $x_i$ and $x_j$ denote the atomic x, y or z coordinates.
2422 In practice, this equation is not used, but the Hessian is
2423 calculated numerically from the force as:
2424 \bea
2425 H_{ij} &=& -
2426 \frac{f_i({\bf x}+h{\bf e}_j) - f_i({\bf x}-h{\bf e}_j)}{2h}
2428 f_i &=& - \frac{\partial V}{\partial x_i}
2429 \eea
2430 where ${\bf e}_j$ is the unit vector in direction $j$.
2431 It should be noted that
2432 for a usual normal-mode calculation, it is necessary to completely minimize
2433 the energy prior to computation of the Hessian.
2434 The tolerance required depends on the type of system,
2435 but a rough indication is 0.001 kJ mol$^{-1}$.
2436 Minimization should be done with conjugate gradients or L-BFGS in double precision.
2438 A number of {\gromacs} programs are involved in these
2439 calculations. First, the energy should be minimized using {\tt mdrun}.
2440 Then, {\tt mdrun} computes the Hessian. {\bf Note} that for generating
2441 the run input file, one should use the minimized conformation from
2442 the full precision trajectory file, as the structure file is not
2443 accurate enough.
2444 {\tt \normindex{g_nmeig}} does the diagonalization and
2445 the sorting of the normal modes according to their frequencies.
2446 Both {\tt mdrun} and {\tt g_nmeig} should be run in double precision.
2447 The normal modes can be analyzed with the program {\tt g_anaeig}.
2448 Ensembles of structures at any temperature and for any subset of
2449 normal modes can be generated with {\tt \normindex{g_nmens}}.
2450 An overview of normal-mode analysis and the related principal component
2451 analysis (see \secref{covanal}) can be found in~\cite{Hayward95b}.
2452 % } % Brace matches ifthenelse test for gmxlite
2454 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2456 \section{Free energy calculations\index{free energy calculations}}
2457 \label{sec:fecalc}
2458 \subsection{Slow-growth methods\index{slow-growth methods}}
2459 Free energy calculations can be performed
2460 in {\gromacs} using a number of methods, including ``slow-growth.'' An example problem
2461 might be calculating the difference in free energy of binding of an inhibitor {\bf I}
2462 to an enzyme {\bf E} and to a mutated enzyme {\bf E$^{\prime}$}. It
2463 is not feasible with computer simulations to perform a docking
2464 calculation for such a large complex, or even releasing the inhibitor from
2465 the enzyme in a reasonable amount of computer time with reasonable accuracy.
2466 However, if we consider the free energy cycle in~\figref{free}A
2467 we can write:
2468 \beq
2469 \Delta G_1 - \Delta G_2 = \Delta G_3 - \Delta G_4
2470 \label{eqn:ddg}
2471 \eeq
2472 If we are interested in the left-hand term we can equally well compute
2473 the right-hand term.
2474 \begin{figure}
2475 \centerline{\includegraphics[width=6cm,angle=270]{plots/free1}\hspace{2cm}\includegraphics[width=6cm,angle=270]{plots/free2}}
2476 \caption[Free energy cycles.]{Free energy cycles. {\bf A:} to
2477 calculate $\Delta G_{12}$, the free energy difference between the
2478 binding of inhibitor {\bf I} to enzymes {\bf E} respectively {\bf
2479 E$^{\prime}$}. {\bf B:} to calculate $\Delta G_{12}$, the free energy
2480 difference for binding of inhibitors {\bf I} respectively {\bf I$^{\prime}$} to
2481 enzyme {\bf E}.}
2482 \label{fig:free}
2483 \end{figure}
2485 If we want to compute the difference in free energy of binding of two
2486 inhibitors {\bf I} and {\bf I$^{\prime}$} to an enzyme {\bf E} (\figref{free}B)
2487 we can again use \eqnref{ddg} to compute the desired property.
2489 \newcommand{\sA}{^{\mathrm{A}}}
2490 \newcommand{\sB}{^{\mathrm{B}}}
2491 Free energy differences between two molecular species can
2492 be calculated in {\gromacs} using the ``slow-growth'' method.
2493 Such free energy differences between different molecular species are
2494 physically meaningless, but they can be used to obtain meaningful
2495 quantities employing a thermodynamic cycle.
2496 The method requires a simulation during which the Hamiltonian of the
2497 system changes slowly from that describing one system (A) to that
2498 describing the other system (B). The change must be so slow that the
2499 system remains in equilibrium during the process; if that requirement
2500 is fulfilled, the change is reversible and a slow-growth simulation from B to A
2501 will yield the same results (but with a different sign) as a slow-growth
2502 simulation from A to B. This is a useful check, but the user should be
2503 aware of the danger that equality of forward and backward growth results does
2504 not guarantee correctness of the results.
2506 The required modification of the Hamiltonian $H$ is realized by making
2507 $H$ a function of a \textit{coupling parameter} $\lambda:
2508 H=H(p,q;\lambda)$ in such a way that $\lambda=0$ describes system A
2509 and $\lambda=1$ describes system B:
2510 \beq
2511 H(p,q;0)=H\sA (p,q);~~~~ H(p,q;1)=H\sB (p,q).
2512 \eeq
2513 In {\gromacs}, the functional form of the $\lambda$-dependence is
2514 different for the various force-field contributions and is described
2515 in section \secref{feia}.
2517 The Helmholtz free energy $A$ is related to the
2518 partition function $Q$ of an $N,V,T$ ensemble, which is assumed to be
2519 the equilibrium ensemble generated by a MD simulation at constant
2520 volume and temperature. The generally more useful Gibbs free energy
2521 $G$ is related to the partition function $\Delta$ of an $N,p,T$
2522 ensemble, which is assumed to be the equilibrium ensemble generated by
2523 a MD simulation at constant pressure and temperature:
2524 \bea
2525 A(\lambda) &=& -k_BT \ln Q \\
2526 Q &=& c \int\!\!\int \exp[-\beta H(p,q;\lambda)]\,dp\,dq \\
2527 G(\lambda) &=& -k_BT \ln \Delta \\
2528 \Delta &=& c \int\!\!\int\!\!\int \exp[-\beta H(p,q;\lambda) -\beta
2529 pV]\,dp\,dq\,dV \\
2530 G &=& A + pV,
2531 \eea
2532 where $\beta = 1/(k_BT)$ and $c = (N! h^{3N})^{-1}$.
2533 These integrals over phase space cannot be evaluated from a
2534 simulation, but it is possible to evaluate the derivative with
2535 respect to $\lambda$ as an ensemble average:
2536 \beq
2537 \frac{dA}{d\lambda} = \frac{\int\!\!\int (\partial H/ \partial
2538 \lambda) \exp[-\beta H(p,q;\lambda)]\,dp\,dq}{\int\!\!\int \exp[-\beta
2539 H(p,q;\lambda)]\,dp\,dq} =
2540 \left\langle \frac{\partial H}{\partial \lambda} \right\rangle_{NVT;\lambda},
2541 \eeq
2542 with a similar relation for $dG/d\lambda$ in the $N,p,T$
2543 ensemble. The difference in free energy between A and B can be found
2544 by integrating the derivative over $\lambda$:
2545 \bea
2546 A\sB(V,T)-A\sA(V,T) &=& \int_0^1 \left\langle \frac{\partial
2547 H}{\partial \lambda} \right\rangle_{NVT;\lambda} \,d\lambda
2548 \label{eq:delA} \\
2549 G\sB(p,T)-G\sA(p,T) &=& \int_0^1 \left\langle \frac{\partial
2550 H}{\partial \lambda} \right\rangle_{NpT;\lambda} \,d\lambda.
2551 \label{eq:delG}
2552 \eea
2553 If one wishes to evaluate $G\sB(p,T)-G\sA(p,T)$,
2554 the natural choice is a constant-pressure simulation. However, this
2555 quantity can also be obtained from a slow-growth simulation at
2556 constant volume, starting with system A at pressure $p$ and volume $V$
2557 and ending with system B at pressure $p_B$, by applying the following
2558 small (but, in principle, exact) correction:
2559 \beq
2560 G\sB(p)-G\sA(p) =
2561 A\sB(V)-A\sA(V) - \int_p^{p\sB}[V\sB(p')-V]\,dp'
2562 \eeq
2563 Here we omitted the constant $T$ from the notation. This correction is
2564 roughly equal to $-\frac{1}{2} (p\sB-p)\Delta V=(\Delta V)^2/(2
2565 \kappa V)$, where $\Delta V$ is the volume change at $p$ and $\kappa$
2566 is the isothermal compressibility. This is usually
2567 small; for example, the growth of a water molecule from nothing
2568 in a bath of 1000 water molecules at constant volume would produce an
2569 additional pressure of as much as 22 bar, but a correction to the
2570 Helmholtz free energy of just -1 kJ mol$^{-1}$. %-20 J/mol.
2572 In Cartesian coordinates, the kinetic energy term in the Hamiltonian
2573 depends only on the momenta, and can be separately integrated and, in
2574 fact, removed from the equations. When masses do not change, there is
2575 no contribution from the kinetic energy at all; otherwise the
2576 integrated contribution to the free energy is $-\frac{3}{2} k_BT \ln
2577 (m\sB/m\sA)$. {\bf Note} that this is only true in the absence of constraints.
2579 \subsection{Thermodynamic integration\index{thermodynamic integration}\index{BAR}\index{Bennett's acceptance ratio}}
2580 {\gromacs} offers the possibility to integrate eq.~\ref{eq:delA} or
2581 eq. \ref{eq:delG} in one simulation over the full range from A to
2582 B. However, if the change is large and insufficient sampling can be
2583 expected, the user may prefer to determine the value of $\langle
2584 dG/d\lambda \rangle$ accurately at a number of well-chosen
2585 intermediate values of $\lambda$. This can easily be done by setting
2586 the stepsize {\tt delta_lambda} to zero. Each simulation can be
2587 equilibrated first, and a proper error estimate can be made for each
2588 value of $dG/d\lambda$ from the fluctuation of $\partial H/\partial
2589 \lambda$. The total free energy change is then determined afterward
2590 by an appropriate numerical integration procedure.
2592 {\gromacs} now also supports the use of Bennett's Acceptance Ratio~\cite{Bennett1976}
2593 for calculating values of $\Delta$G for transformations from state A to state B using
2594 the program {\tt \normindex{g_bar}}. The same data can also be used to calculate free
2595 energies using MBAR~\cite{Shirts2008}, though the analysis currently requires external tools from
2596 the external {\tt pymbar} package, at https://SimTK.org/home/pymbar.
2598 The $\lambda$-dependence for the force-field contributions is
2599 described in detail in section \secref{feia}.
2600 % } % Brace matches ifthenelse test for gmxlite
2602 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2603 \section{Replica exchange\index{replica exchange}}
2604 Replica exchange molecular dynamics (\normindex{REMD})
2605 is a method that can be used to speed up
2606 the sampling of any type of simulation, especially if
2607 conformations are separated by relatively high energy barriers.
2608 It involves simulating multiple replicas of the same system
2609 at different temperatures and randomly exchanging the complete state
2610 of two replicas at regular intervals with the probability:
2611 \beq
2612 P(1 \leftrightarrow 2)=\min\left(1,\exp\left[
2613 \left(\frac{1}{k_B T_1} - \frac{1}{k_B T_2}\right)(U_1 - U_2)
2614 \right] \right)
2615 \eeq
2616 where $T_1$ and $T_2$ are the reference temperatures and $U_1$ and $U_2$
2617 are the instantaneous potential energies of replicas 1 and 2 respectively.
2618 After exchange the velocities are scaled by $(T_1/T_2)^{\pm0.5}$
2619 and a neighbor search is performed the next step.
2620 This combines the fast sampling and frequent barrier-crossing
2621 of the highest temperature with correct Boltzmann sampling at
2622 all the different temperatures~\cite{Hukushima96a,Sugita99}.
2623 We only attempt exchanges for neighboring temperatures as the probability
2624 decreases very rapidly with the temperature difference.
2625 One should not attempt exchanges for all possible pairs in one step.
2626 If, for instance, replicas 1 and 2 would exchange, the chance of
2627 exchange for replicas 2 and 3 not only depends on the energies of
2628 replicas 2 and 3, but also on the energy of replica 1.
2629 In {\gromacs} this is solved by attempting exchange for all ``odd''
2630 pairs on ``odd'' attempts and for all ``even'' pairs on ``even'' attempts.
2631 If we have four replicas: 0, 1, 2 and 3, ordered in temperature
2632 and we attempt exchange every 1000 steps, pairs 0-1 and 2-3
2633 will be tried at steps 1000, 3000 etc. and pair 1-2 at steps 2000, 4000 etc.
2635 How should one choose the temperatures?
2636 The energy difference can be written as:
2637 \beq
2638 U_1 - U_2 = N_{df} \frac{c}{2} k_B (T_1 - T_2)
2639 \eeq
2640 where $N_{df}$ is the total number of degrees of freedom of one replica
2641 and $c$ is 1 for harmonic potentials and around 2 for protein/water systems.
2642 If $T_2 = (1+\epsilon) T_1$ the probability becomes:
2643 \beq
2644 P(1 \leftrightarrow 2)
2645 = \exp\left( -\frac{\epsilon^2 c\,N_{df}}{2 (1+\epsilon)} \right)
2646 \approx \exp\left(-\epsilon^2 \frac{c}{2} N_{df} \right)
2647 \eeq
2648 Thus for a probability of $e^{-2}\approx 0.135$
2649 one obtains $\epsilon \approx 2/\sqrt{c\,N_{df}}$.
2650 With all bonds constrained one has $N_{df} \approx 2\, N_{atoms}$
2651 and thus for $c$ = 2 one should choose $\epsilon$ as $1/\sqrt{N_{atoms}}$.
2652 However there is one problem when using pressure coupling. The density at
2653 higher temperatures will decrease, leading to higher energy~\cite{Seibert2005a},
2654 which should be taken into account. The {\gromacs} website features a
2655 so-called ``REMD calculator,'' that lets you type in the temperature range and
2656 the number of atoms, and based on that proposes a set of temperatures.
2658 An extension to the REMD for the isobaric-isothermal ensemble was
2659 proposed by Okabe {\em et al.}~\cite{Okabe2001a}. In this work the
2660 exchange probability is modified to:
2661 \beq
2662 P(1 \leftrightarrow 2)=\min\left(1,\exp\left[
2663 \left(\frac{1}{k_B T_1} - \frac{1}{k_B T_2}\right)(U_1 - U_2) +
2664 \left(\frac{P_1}{k_B T_1} - \frac{P_2}{k_B T_2}\right)\left(V_1-V_2\right)
2665 \right] \right)
2666 \eeq
2667 where $P_1$ and $P_2$ are the respective reference pressures and $V_1$ and
2668 $V_2$ are the respective instantaneous volumes in the simulations.
2669 In most cases the differences in volume are so small that the second
2670 term is negligible. It only plays a role when the difference between
2671 $P_1$ and $P_2$ is large or in phase transitions.
2673 Hamiltonian replica exchange is also supported in {\gromacs}. In
2674 Hamiltonian replica exchange, each replica has a different
2675 Hamiltonian, defined by the free energy pathway specified for the simulation. The
2676 exchange probability to maintain the correct ensemble probabilities is:
2677 \beq P(1 \leftrightarrow 2)=\min\left(1,\exp\left[
2678 \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)))
2679 \right]
2680 \right)
2681 \eeq
2682 The separate Hamiltonians are defined by the free energy functionality
2683 of {\gromacs}, with swaps made between the different values of
2684 $\lambda$ defined in the mdp file.
2686 Hamiltonian and temperature replica exchange can also be performed
2687 simultaneously, using the acceptance criteria:
2688 \beq
2689 P(1 \leftrightarrow 2)=\min\left(1,\exp\left[
2690 \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})
2691 \right] \right)
2692 \eeq
2694 Gibbs sampling replica exchange has also been implemented in
2695 {\gromacs}~\cite{Chodera2011}. In Gibbs sampling replica exchange, all
2696 possible pairs are tested for exchange, allowing swaps between
2697 replicas that are not neighbors.
2699 Gibbs sampling replica exchange requires no additional potential
2700 energy calculations. However there is an additional communication
2701 cost in Gibbs sampling replica exchange, as for some permutations,
2702 more than one round of swaps must take place. In some cases, this
2703 extra communication cost might affect the efficiency.
2705 All replica exchange variants are options of the {\tt mdrun}
2706 program. It will only work when MPI is installed, due to the inherent
2707 parallelism in the algorithm. For efficiency each replica can run on a
2708 separate rank. See the manual page of {\tt mdrun} on how to use these
2709 multinode features.
2711 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2713 \section{Essential Dynamics sampling\index{essential dynamics}\index{principal component analysis}\seeindexquiet{PCA}{covariance analysis}}
2714 The results from Essential Dynamics (see \secref{covanal})
2715 of a protein can be used to guide MD simulations. The idea is that
2716 from an initial MD simulation (or from other sources) a definition of
2717 the collective fluctuations with largest amplitude is obtained. The
2718 position along one or more of these collective modes can be
2719 constrained in a (second) MD simulation in a number of ways for
2720 several purposes. For example, the position along a certain mode may
2721 be kept fixed to monitor the average force (free-energy gradient) on
2722 that coordinate in that position. Another application is to enhance
2723 sampling efficiency with respect to usual MD
2724 \cite{Degroot96a,Degroot96b}. In this case, the system is encouraged
2725 to sample its available configuration space more systematically than
2726 in a diffusion-like path that proteins usually take.
2728 Another possibility to enhance sampling is \normindex{flooding}.
2729 Here a flooding potential is added to certain
2730 (collective) degrees of freedom to expel the system out
2731 of a region of phase space \cite{Lange2006a}.
2733 The procedure for essential dynamics sampling or flooding is as follows.
2734 First, the eigenvectors and eigenvalues need to be determined
2735 using covariance analysis ({\tt g_covar})
2736 or normal-mode analysis ({\tt g_nmeig}).
2737 Then, this information is fed into {\tt make_edi},
2738 which has many options for selecting vectors and setting parameters,
2739 see {\tt gmx make_edi -h}.
2740 The generated {\tt edi} input file is then passed to {\tt mdrun}.
2742 % } % Brace matches ifthenelse test for gmxlite
2744 % \ifthenelse{\equal{\gmxlite}{1}}{}{
2745 \section{\normindex{Expanded Ensemble}}
2747 In an expanded ensemble simulation~\cite{Lyubartsev1992}, both the coordinates and the
2748 thermodynamic ensemble are treated as configuration variables that can
2749 be sampled over. The probability of any given state can be written as:
2750 \beq
2751 P(\vec{x},k) \propto \exp\left(-\beta_k U_k + g_k\right),
2752 \eeq
2753 where $\beta_k = \frac{1}{k_B T_k}$ is the $\beta$ corresponding to the $k$th
2754 thermodynamic state, and $g_k$ is a user-specified weight factor corresponding
2755 to the $k$th state. This space is therefore a {\em mixed}, {\em generalized}, or {\em
2756 expanded} ensemble which samples from multiple thermodynamic
2757 ensembles simultaneously. $g_k$ is chosen to give a specific weighting
2758 of each subensemble in the expanded ensemble, and can either be fixed,
2759 or determined by an iterative procedure. The set of $g_k$ is
2760 frequently chosen to give each thermodynamic ensemble equal
2761 probability, in which case $g_k$ is equal to the free energy in
2762 non-dimensional units, but they can be set to arbitrary values as
2763 desired. Several different algorithms can be used to equilibrate
2764 these weights, described in the mdp option listings.
2765 % } % Brace matches ifthenelse test for gmxlite
2767 In {\gromacs}, this space is sampled by alternating sampling in the $k$
2768 and $\vec{x}$ directions. Sampling in the $\vec{x}$ direction is done
2769 by standard molecular dynamics sampling; sampling between the
2770 different thermodynamics states is done by Monte Carlo, with several
2771 different Monte Carlo moves supported. The $k$ states can be defined
2772 by different temperatures, or choices of the free energy $\lambda$
2773 variable, or both. Expanded ensemble simulations thus represent a
2774 serialization of the replica exchange formalism, allowing a single
2775 simulation to explore many thermodynamic states.
2779 \section{Parallelization\index{parallelization}}
2780 The CPU time required for a simulation can be reduced by running the simulation
2781 in parallel over more than one core.
2782 Ideally, one would want to have linear scaling: running on $N$ cores
2783 makes the simulation $N$ times faster. In practice this can only be
2784 achieved for a small number of cores. The scaling will depend
2785 a lot on the algorithms used. Also, different algorithms can have different
2786 restrictions on the interaction ranges between atoms.
2788 \section{Domain decomposition\index{domain decomposition}}
2789 Since most interactions in molecular simulations are local,
2790 domain decomposition is a natural way to decompose the system.
2791 In domain decomposition, a spatial domain is assigned to each rank,
2792 which will then integrate the equations of motion for the particles
2793 that currently reside in its local domain. With domain decomposition,
2794 there are two choices that have to be made: the division of the unit cell
2795 into domains and the assignment of the forces to domains.
2796 Most molecular simulation packages use the half-shell method for assigning
2797 the forces. But there are two methods that always require less communication:
2798 the eighth shell~\cite{Liem1991} and the midpoint~\cite{Shaw2006} method.
2799 {\gromacs} currently uses the eighth shell method, but for certain systems
2800 or hardware architectures it might be advantageous to use the midpoint
2801 method. Therefore, we might implement the midpoint method in the future.
2802 Most of the details of the domain decomposition can be found
2803 in the {\gromacs} 4 paper~\cite{Hess2008b}.
2805 \subsection{Coordinate and force communication}
2806 In the most general case of a triclinic unit cell,
2807 the space in divided with a 1-, 2-, or 3-D grid in parallelepipeds
2808 that we call domain decomposition cells.
2809 Each cell is assigned to a particle-particle rank.
2810 The system is partitioned over the ranks at the beginning
2811 of each MD step in which neighbor searching is performed.
2812 Since the neighbor searching is based on charge groups, charge groups
2813 are also the units for the domain decomposition.
2814 Charge groups are assigned to the cell where their center of geometry resides.
2815 Before the forces can be calculated, the coordinates from some
2816 neighboring cells need to be communicated,
2817 and after the forces are calculated, the forces need to be communicated
2818 in the other direction.
2819 The communication and force assignment is based on zones that
2820 can cover one or multiple cells.
2821 An example of a zone setup is shown in \figref{ddcells}.
2823 \begin{figure}
2824 \centerline{\includegraphics[width=6cm]{plots/dd-cells}}
2825 \caption{
2826 A non-staggered domain decomposition grid of 3$\times$2$\times$2 cells.
2827 Coordinates in zones 1 to 7 are communicated to the corner cell
2828 that has its home particles in zone 0.
2829 $r_c$ is the cut-off radius.
2830 \label{fig:ddcells}
2832 \end{figure}
2834 The coordinates are communicated by moving data along the ``negative''
2835 direction in $x$, $y$ or $z$ to the next neighbor. This can be done in one
2836 or multiple pulses. In \figref{ddcells} two pulses in $x$ are required,
2837 then one in $y$ and then one in $z$. The forces are communicated by
2838 reversing this procedure. See the {\gromacs} 4 paper~\cite{Hess2008b}
2839 for details on determining which non-bonded and bonded forces
2840 should be calculated on which rank.
2842 \subsection{Dynamic load balancing\swapindexquiet{dynamic}{load balancing}}
2843 When different ranks have a different computational load
2844 (load imbalance), all ranks will have to wait for the one
2845 that takes the most time. One would like to avoid such a situation.
2846 Load imbalance can occur due to three reasons:
2847 \begin{itemize}
2848 \item inhomogeneous particle distribution
2849 \item inhomogeneous interaction cost distribution (charged/uncharged,
2850 water/non-water due to {\gromacs} water innerloops)
2851 \item statistical fluctuation (only with small particle numbers)
2852 \end{itemize}
2853 So we need a dynamic load balancing algorithm
2854 where the volume of each domain decomposition cell
2855 can be adjusted {\em independently}.
2856 To achieve this, the 2- or 3-D domain decomposition grids need to be
2857 staggered. \figref{ddtric} shows the most general case in 2-D.
2858 Due to the staggering, one might require two distance checks
2859 for deciding if a charge group needs to be communicated:
2860 a non-bonded distance and a bonded distance check.
2862 \begin{figure}
2863 \centerline{\includegraphics[width=7cm]{plots/dd-tric}}
2864 \caption{
2865 The zones to communicate to the rank of zone 0,
2866 see the text for details. $r_c$ and $r_b$ are the non-bonded
2867 and bonded cut-off radii respectively, $d$ is an example
2868 of a distance between following, staggered boundaries of cells.
2869 \label{fig:ddtric}
2871 \end{figure}
2873 By default, {\tt mdrun} automatically turns on the dynamic load
2874 balancing during a simulation when the total performance loss
2875 due to the force calculation imbalance is 5\% or more.
2876 {\bf Note} that the reported force load imbalance numbers might be higher,
2877 since the force calculation is only part of work that needs to be done
2878 during an integration step.
2879 The load imbalance is reported in the log file at log output steps
2880 and when the {\tt -v} option is used also on screen.
2881 The average load imbalance and the total performance loss
2882 due to load imbalance are reported at the end of the log file.
2884 There is one important parameter for the dynamic load balancing,
2885 which is the minimum allowed scaling. By default, each dimension
2886 of the domain decomposition cell can scale down by at least
2887 a factor of 0.8. For 3-D domain decomposition this allows cells
2888 to change their volume by about a factor of 0.5, which should allow
2889 for compensation of a load imbalance of 100\%.
2890 The minimum allowed scaling can be changed with the {\tt -dds}
2891 option of {\tt mdrun}.
2893 \subsection{Constraints in parallel\index{constraints}}
2894 \label{subsec:plincs}
2895 Since with domain decomposition parts of molecules can reside
2896 on different ranks, bond constraints can cross cell boundaries.
2897 Therefore a parallel constraint algorithm is required.
2898 {\gromacs} uses the \normindex{P-LINCS} algorithm~\cite{Hess2008a},
2899 which is the parallel version of the \normindex{LINCS} algorithm~\cite{Hess97}
2900 % \ifthenelse{\equal{\gmxlite}{1}}
2902 {(see \ssecref{lincs}).}
2903 The P-LINCS procedure is illustrated in \figref{plincs}.
2904 When molecules cross the cell boundaries, atoms in such molecules
2905 up to ({\tt lincs_order + 1}) bonds away are communicated over the cell boundaries.
2906 Then, the normal LINCS algorithm can be applied to the local bonds
2907 plus the communicated ones. After this procedure, the local bonds
2908 are correctly constrained, even though the extra communicated ones are not.
2909 One coordinate communication step is required for the initial LINCS step
2910 and one for each iteration. Forces do not need to be communicated.
2912 \begin{figure}
2913 \centerline{\includegraphics[width=6cm]{plots/par-lincs2}}
2914 \caption{
2915 Example of the parallel setup of P-LINCS with one molecule
2916 split over three domain decomposition cells, using a matrix
2917 expansion order of 3.
2918 The top part shows which atom coordinates need to be communicated
2919 to which cells. The bottom parts show the local constraints (solid)
2920 and the non-local constraints (dashed) for each of the three cells.
2921 \label{fig:plincs}
2923 \end{figure}
2925 \subsection{Interaction ranges}
2926 Domain decomposition takes advantage of the locality of interactions.
2927 This means that there will be limitations on the range of interactions.
2928 By default, {\tt mdrun} tries to find the optimal balance between
2929 interaction range and efficiency. But it can happen that a simulation
2930 stops with an error message about missing interactions,
2931 or that a simulation might run slightly faster with shorter
2932 interaction ranges. A list of interaction ranges
2933 and their default values is given in \tabref{dd_ranges}.
2935 \begin{table}
2936 \centerline{
2937 \begin{tabular}{|c|c|ll|}
2938 \dline
2939 interaction & range & option & default \\
2940 \dline
2941 non-bonded & $r_c$ = max($r_{\mathrm{list}}$,$r_{\mathrm{VdW}}$,$r_{\mathrm{Coul}}$) & {\tt mdp} file & \\
2942 two-body bonded & max($r_{\mathrm{mb}}$,$r_c$) & {\tt mdrun -rdd} & starting conf. + 10\% \\
2943 multi-body bonded & $r_{\mathrm{mb}}$ & {\tt mdrun -rdd} & starting conf. + 10\% \\
2944 constraints & $r_{\mathrm{con}}$ & {\tt mdrun -rcon} & est. from bond lengths \\
2945 virtual sites & $r_{\mathrm{con}}$ & {\tt mdrun -rcon} & 0 \\
2946 \dline
2947 \end{tabular}
2949 \caption{The interaction ranges with domain decomposition.}
2950 \label{tab:dd_ranges}
2951 \end{table}
2953 In most cases the defaults of {\tt mdrun} should not cause the simulation
2954 to stop with an error message of missing interactions.
2955 The range for the bonded interactions is determined from the distance
2956 between bonded charge-groups in the starting configuration, with 10\% added
2957 for headroom. For the constraints, the value of $r_{\mathrm{con}}$ is determined by
2958 taking the maximum distance that ({\tt lincs_order + 1}) bonds can cover
2959 when they all connect at angles of 120 degrees.
2960 The actual constraint communication is not limited by $r_{\mathrm{con}}$,
2961 but by the minimum cell size $L_C$, which has the following lower limit:
2962 \beq
2963 L_C \geq \max(r_{\mathrm{mb}},r_{\mathrm{con}})
2964 \eeq
2965 Without dynamic load balancing the system is actually allowed to scale
2966 beyond this limit when pressure scaling is used.
2967 {\bf Note} that for triclinic boxes, $L_C$ is not simply the box diagonal
2968 component divided by the number of cells in that direction,
2969 rather it is the shortest distance between the triclinic cells borders.
2970 For rhombic dodecahedra this is a factor of $\sqrt{3/2}$ shorter
2971 along $x$ and $y$.
2973 When $r_{\mathrm{mb}} > r_c$, {\tt mdrun} employs a smart algorithm to reduce
2974 the communication. Simply communicating all charge groups within
2975 $r_{\mathrm{mb}}$ would increase the amount of communication enormously.
2976 Therefore only charge-groups that are connected by bonded interactions
2977 to charge groups which are not locally present are communicated.
2978 This leads to little extra communication, but also to a slightly
2979 increased cost for the domain decomposition setup.
2980 In some cases, {\eg} coarse-grained simulations with a very short cut-off,
2981 one might want to set $r_{\mathrm{mb}}$ by hand to reduce this cost.
2983 \subsection{Multiple-Program, Multiple-Data PME parallelization\index{PME}}
2984 \label{subsec:mpmd_pme}
2985 Electrostatics interactions are long-range, therefore special
2986 algorithms are used to avoid summation over many atom pairs.
2987 In {\gromacs} this is usually
2988 % \ifthenelse{\equal{\gmxlite}{1}}
2990 {PME (\secref{pme}).}
2991 Since with PME all particles interact with each other, global communication
2992 is required. This will usually be the limiting factor for
2993 scaling with domain decomposition.
2994 To reduce the effect of this problem, we have come up with
2995 a Multiple-Program, Multiple-Data approach~\cite{Hess2008b}.
2996 Here, some ranks are selected to do only the PME mesh calculation,
2997 while the other ranks, called particle-particle (PP) ranks,
2998 do all the rest of the work.
2999 For rectangular boxes the optimal PP to PME rank ratio is usually 3:1,
3000 for rhombic dodecahedra usually 2:1.
3001 When the number of PME ranks is reduced by a factor of 4, the number
3002 of communication calls is reduced by about a factor of 16.
3003 Or put differently, we can now scale to 4 times more ranks.
3004 In addition, for modern 4 or 8 core machines in a network,
3005 the effective network bandwidth for PME is quadrupled,
3006 since only a quarter of the cores will be using the network connection
3007 on each machine during the PME calculations.
3009 \begin{figure}
3010 \centerline{\includegraphics[width=12cm]{plots/mpmd-pme}}
3011 \caption{
3012 Example of 8 ranks without (left) and with (right) MPMD.
3013 The PME communication (red arrows) is much higher on the left
3014 than on the right. For MPMD additional PP - PME coordinate
3015 and force communication (blue arrows) is required,
3016 but the total communication complexity is lower.
3017 \label{fig:mpmd_pme}
3019 \end{figure}
3021 {\tt mdrun} will by default interleave the PP and PME ranks.
3022 If the ranks are not number consecutively inside the machines,
3023 one might want to use {\tt mdrun -ddorder pp_pme}.
3024 For machines with a real 3-D torus and proper communication software
3025 that assigns the ranks accordingly one should use
3026 {\tt mdrun -ddorder cartesian}.
3028 To optimize the performance one should usually set up the cut-offs
3029 and the PME grid such that the PME load is 25 to 33\% of the total
3030 calculation load. {\tt grompp} will print an estimate for this load
3031 at the end and also {\tt mdrun} calculates the same estimate
3032 to determine the optimal number of PME ranks to use.
3033 For high parallelization it might be worthwhile to optimize
3034 the PME load with the {\tt mdp} settings and/or the number
3035 of PME ranks with the {\tt -npme} option of {\tt mdrun}.
3036 For changing the electrostatics settings it is useful to know
3037 the accuracy of the electrostatics remains nearly constant
3038 when the Coulomb cut-off and the PME grid spacing are scaled
3039 by the same factor.
3040 {\bf Note} that it is usually better to overestimate than to underestimate
3041 the number of PME ranks, since the number of PME ranks is smaller
3042 than the number of PP ranks, which leads to less total waiting time.
3044 The PME domain decomposition can be 1-D or 2-D along the $x$ and/or
3045 $y$ axis. 2-D decomposition is also known as \normindex{pencil decomposition} because of
3046 the shape of the domains at high parallelization.
3047 1-D decomposition along the $y$ axis can only be used when
3048 the PP decomposition has only 1 domain along $x$. 2-D PME decomposition
3049 has to have the number of domains along $x$ equal to the number of
3050 the PP decomposition. {\tt mdrun} automatically chooses 1-D or 2-D
3051 PME decomposition (when possible with the total given number of ranks),
3052 based on the minimum amount of communication for the coordinate redistribution
3053 in PME plus the communication for the grid overlap and transposes.
3054 To avoid superfluous communication of coordinates and forces
3055 between the PP and PME ranks, the number of DD cells in the $x$
3056 direction should ideally be the same or a multiple of the number
3057 of PME ranks. By default, {\tt mdrun} takes care of this issue.
3059 \subsection{Domain decomposition flow chart}
3060 In \figref{dd_flow} a flow chart is shown for domain decomposition
3061 with all possible communication for different algorithms.
3062 For simpler simulations, the same flow chart applies,
3063 without the algorithms and communication for
3064 the algorithms that are not used.
3066 \begin{figure}
3067 \centerline{\includegraphics[width=12cm]{plots/flowchart}}
3068 \caption{
3069 Flow chart showing the algorithms and communication (arrows)
3070 for a standard MD simulation with virtual sites, constraints
3071 and separate PME-mesh ranks.
3072 \label{fig:dd_flow}
3074 \end{figure}
3077 \section{Implicit solvation\index{implicit solvation}\index{Generalized Born methods}}
3078 \label{sec:gbsa}
3079 Implicit solvent models provide an efficient way of representing
3080 the electrostatic effects of solvent molecules, while saving a
3081 large piece of the computations involved in an accurate, aqueous
3082 description of the surrounding water in molecular dynamics simulations.
3083 Implicit solvation models offer several advantages compared with
3084 explicit solvation, including eliminating the need for the equilibration of water
3085 around the solute, and the absence of viscosity, which allows the protein
3086 to more quickly explore conformational space.
3088 Implicit solvent calculations in {\gromacs} can be done using the
3089 generalized Born-formalism, and the Still~\cite{Still97}, HCT~\cite{Truhlar96},
3090 and OBC~\cite{Case04} models are available for calculating the Born radii.
3092 Here, the free energy $G_{\mathrm{solv}}$ of solvation is the sum of three terms,
3093 a solvent-solvent cavity term ($G_{\mathrm{cav}}$), a solute-solvent van der
3094 Waals term ($G_{\mathrm{vdw}}$), and finally a solvent-solute electrostatics
3095 polarization term ($G_{\mathrm{pol}}$).
3097 The sum of $G_{\mathrm{cav}}$ and $G_{\mathrm{vdw}}$ corresponds to the (non-polar)
3098 free energy of solvation for a molecule from which all charges
3099 have been removed, and is commonly called $G_{\mathrm{np}}$,
3100 calculated from the total solvent accessible surface area
3101 multiplied with a surface tension.
3102 The total expression for the solvation free energy then becomes:
3104 \beq
3105 G_{\mathrm{solv}} = G_{\mathrm{np}} + G_{\mathrm{pol}}
3106 \label{eqn:gb_solv}
3107 \eeq
3109 Under the generalized Born model, $G_{\mathrm{pol}}$ is calculated from the generalized Born equation~\cite{Still97}:
3111 \beq
3112 G_{\mathrm{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)}}
3113 \label{eqn:gb_still}
3114 \eeq
3116 In {\gromacs}, we have introduced the substitution~\cite{Larsson10}:
3118 \beq
3119 c_i=\frac{1}{\sqrt{b_i}}
3120 \label{eqn:gb_subst}
3121 \eeq
3123 which makes it possible to introduce a cheap transformation to a new
3124 variable $x$ when evaluating each interaction, such that:
3126 \beq
3127 x=\frac{r_{ij}}{\sqrt{b_i b_j }} = r_{ij} c_i c_j
3128 \label{eqn:gb_subst2}
3129 \eeq
3131 In the end, the full re-formulation of~\ref{eqn:gb_still} becomes:
3133 \beq
3134 G_{\mathrm{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)
3135 \label{eqn:gb_final}
3136 \eeq
3138 The non-polar part ($G_{\mathrm{np}}$) of Equation~\ref{eqn:gb_solv} is calculated
3139 directly from the Born radius of each atom using a simple ACE type
3140 approximation by Schaefer {\em et al.}~\cite{Karplus98}, including a
3141 simple loop over all atoms.
3142 This requires only one extra solvation parameter, independent of atom type,
3143 but differing slightly between the three Born radii models.
3145 % LocalWords: GROningen MAchine BIOSON Groningen GROMACS Berendsen der Spoel
3146 % LocalWords: Drunen Comp Phys Comm ROck NS FFT pbc EM ifthenelse gmxlite ff
3147 % LocalWords: octahedra triclinic Ewald PME PPPM trjconv xy solvated
3148 % LocalWords: boxtypes boxshapes editconf Lennard COM XTC TNG kT defunits
3149 % LocalWords: Boltzmann's Mueller nb int mdrun chargegroup simplerc prefactor
3150 % LocalWords: pme waterloops CH NH CO df com virial integrator Verlet vverlet
3151 % LocalWords: integrators ref timepoint timestep timesteps mdp md vv avek NVE
3152 % LocalWords: NVT off's leapfrogv lll LR rmfast SPC fs Nos physicality ps GMX
3153 % LocalWords: Tcoupling nonergodic thermostatting NOSEHOOVER algorithmes ij yx
3154 % LocalWords: Parrinello Rahman rescales atm anisotropically ccc xz zx yy yz
3155 % LocalWords: zy zz se barostat compressibilities MTTK NPT Martyna al isobaric
3156 % LocalWords: Tuckerman vir PV fkT iLt iL Liouville NHC Eq baro mu trj mol bc
3157 % LocalWords: freezegroup Shannon's polarizability Overhauser barostats iLn KE
3158 % LocalWords: negligibly thermostatted Tobias rhombic maxwell et xtc tng TC rlist
3159 % LocalWords: waals LINCS holonomic plincs lincs unc ang SA Langevin SD amu BD
3160 % LocalWords: bfgs Broyden Goldfarb Shanno mkT kJ DFLEXIBLE Nocedal diag nmeig
3161 % LocalWords: diagonalization anaeig nmens covanal ddg feia BT dp dq pV dV dA
3162 % LocalWords: NpT eq stepsize REMD constrainted website Okabe MPI covar edi dd
3163 % LocalWords: progman NMR ddcells innerloops ddtric tric dds rdd conf rcon est
3164 % LocalWords: mb PP MPMD ddorder pp cartesian grompp npme parallelizable edr
3165 % LocalWords: macromolecule nstlist vacuo parallelization dof indices MBAR AVX
3166 % LocalWords: TOL numerics parallelized eigenvectors dG parallelepipeds VdW np
3167 % LocalWords: Coul multi solvation HCT OBC solv cav vdw Schaefer symplectic dt
3168 % LocalWords: pymbar multinode subensemble Monte solute subst groupconcept GPU
3169 % LocalWords: dodecahedron octahedron dodecahedra equilibration usinggroups nm
3170 % LocalWords: topologies rlistlong CUDA GPUs rcoulomb SIMD BlueGene FPUs erfc
3171 % LocalWords: cutoffschemesupport unbuffered bondeds AdResS OpenMP ewald rtol
3172 % LocalWords: verletdrift peptide RMS rescaling ergodicity ergodic discretized
3173 % LocalWords: isothermal compressibility isotropically anisotropic iteratively
3174 % LocalWords: incompressible integrations translational biomolecules NMA PCA
3175 % LocalWords: Bennett's equilibrated Hamiltonians covariance equilibrate
3176 % LocalWords: inhomogeneous conformational online other's th