.. _md:
Molecular Dynamics

.. _gmxmdscheme:
**THE GLOBAL MD ALGORITHM**


 **1. Input initial conditions**
 Potential interaction :math:`V` as a function of atom positions
 Positions :math:`\mathbf{r}` of all atoms in the system
 Velocities :math:`\mathbf{v}` of all atoms in the system
 :math:`\Downarrow`


 **repeat 2,3,4** for the required number of steps:


 **2. Compute forces**
 The force on any atom
 :math:`\mathbf{F}_i =  \displaystyle\frac{\partial V}{\partial \mathbf{r}_i}`
 is computed by calculating the force between nonbonded atom pairs:
 :math:`\mathbf{F}_i = \sum_j \mathbf{F}_{ij}`
 plus the forces due to bonded interactions (which may depend on 1, 2,
3, or 4 atoms), plus restraining and/or external forces.
 The potential and kinetic energies and the pressure tensor may be
computed.
 :math:`\Downarrow`
 **3. Update configuration**
 The movement of the atoms is simulated by numerically solving Newton’s
equations of motion
 :math:`\displaystyle \frac {{\mbox{d}}^2\mathbf{r}_i}{{\mbox{d}}t^2} = \frac{\mathbf{F}_i}{m_i}`
 or
 :math:`\displaystyle \frac{{\mbox{d}}\mathbf{r}_i}{{\mbox{d}}t} = \mathbf{v}_i ; \;\; \frac{{\mbox{d}}\mathbf{v}_i}{{\mbox{d}}t} = \frac{\mathbf{F}_i}{m_i}`
 :math:`\Downarrow`
 **4.** if required: **Output step**
 write positions, velocities, energies, temperature, pressure, etc.
A global flow scheme for MD is given above.
Each MD or EM run requires as input
a set of initial coordinates and – optionally – initial velocities of
all particles involved. This chapter does not describe how these are
obtained; for the setup of an actual MD run check the :ref:`user guide`
in Sections :ref:`gmxsysprep` and :ref:`gmxgettingstarted`.
Initial conditions
~~~~~~~~~~~~~~~~~~
Topology and force field
^^^^^^^^^^^^^^^^^^^^^^^^
The system topology, including a description of the force field, must be
read in. Force fields and topologies are described in chapter :ref:`ff`
and :ref:`top`, respectively. All this information is static; it is never
modified during the run.
Coordinates and velocities
^^^^^^^^^^^^^^^^^^^^^^^^^^
.. _figmaxwell:
.. figure:: plots/maxwell.*
:width: 8.00000cm
A MaxwellBoltzmann velocity distribution, generated from
random numbers.
Then, before a run starts, the box size and the coordinates and
velocities of all particles are required. The box size and shape is
determined by three vectors (nine numbers)
:math:`\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3`,
which represent the three basis vectors of the periodic box.
If the run starts at :math:`t=t_0`, the coordinates at :math:`t=t_0`
must be known. The *leapfrog algorithm*, the default algorithm used to
update the time step with :math:`{{\Delta t}}` (see :ref:`update`),
also requires that the velocities at
:math:`t=t_0  {{\frac{1}{2}}{{\Delta t}}}` are known. If velocities are
not available, the program can generate initial atomic velocities
:math:`v_i, i=1\ldots 3N` with a MaxwellBoltzmann distribution
(:numref:`Fig. %s `) at a given absolute temperature
:math:`T`:
.. math:: p(v_i) = \sqrt{\frac{m_i}{2 \pi kT}}\exp\left(\frac{m_i v_i^2}{2kT}\right)
:label: eqnmaxwellboltzman
where :math:`k` is Boltzmann’s constant (see chapter :ref:`defunits`). To
accomplish this, normally distributed random numbers are generated by
adding twelve random numbers :math:`R_k` in the range
:math:`0 \le R_k < 1` and subtracting 6.0 from their sum. The result is
then multiplied by the standard deviation of the velocity distribution
:math:`\sqrt{kT/m_i}`. Since the resulting total energy will not
correspond exactly to the required temperature :math:`T`, a correction
is made: first the centerofmass motion is removed and then all
velocities are scaled so that the total energy corresponds exactly to
:math:`T` (see :eq:`eqn. %s `).
Centerofmass motion
^^^^^^^^^^^^^^^^^^^^^
The centerofmass velocity is normally set to zero at every step; there
is (usually) no net external force acting on the system and the
centerofmass velocity should remain constant. In practice, however,
the update algorithm introduces a very slow change in the centerofmass
velocity, and therefore in the total kinetic energy of the system –
especially when temperature coupling is used. If such changes are not
quenched, an appreciable centerofmass motion can develop in long runs,
and the temperature will be significantly misinterpreted. Something
similar may happen due to overall rotational motion, but only when an
isolated cluster is simulated. In periodic systems with filled boxes,
the overall rotational motion is coupled to other degrees of freedom and
does not cause such problems.
Neighbor searching
~~~~~~~~~~~~~~~~~~
As mentioned in chapter :ref:`ff`, internal forces are either generated
from fixed (static) lists, or from dynamic lists. The latter consist of
nonbonded interactions between any pair of particles. When calculating
the nonbonded forces, it is convenient to have all particles in a
rectangular box. As shown in :numref:`Fig. %s `, it is possible to transform
a triclinic box into a rectangular box. The output coordinates are
always in a rectangular box, even when a dodecahedron or triclinic box
was used for the simulation. :eq:`Equation %s ` ensures that we can
reset particles in a rectangular box by first shifting them with box
vector :math:`{\bf c}`, then with :math:`{\bf b}` and finally with
:math:`{\bf a}`. Equations :eq:`%s `,
:eq:`%s ` and :eq:`%s `
ensure that we can find the 14 nearest triclinic images within a linear
combination that does not involve multiples of box vectors.
Pair lists generation
^^^^^^^^^^^^^^^^^^^^^
The nonbonded pair forces need to be calculated only for those pairs
:math:`i,j` for which the distance :math:`r_{ij}` between :math:`i` and
the nearest image of :math:`j` is less than a given cutoff radius
:math:`R_c`. Some of the particle pairs that fulfill this criterion are
excluded, when their interaction is already fully accounted for by
bonded interactions. But for most electrostatic treatments, correction
forces also need to be computed for such excluded atom pairs.
Gromacs employs a *pair list* that contains those
particle pairs for which nonbonded forces must be calculated. The pair
list contains particles :math:`i`, a displacement vector for particle
:math:`i`, and all particles :math:`j` that are within ``rlist`` of this
particular image of particle :math:`i`. The list is updated every
``nstlist`` steps.
To make the pair list, all atom pairs that are within the pair list
cutoff distance need to be found and stored in a list. Note that
such a list generally does not store all neighbors for each atom,
since each atom pair should appear only once in the list. This
searching, usually called neighbor search (NS) or pair search, involves
periodic boundary conditions and determining the *image* (see
sec. :ref:`pbc`). The search algorithm employed in Gromacs is :math:`O(N)`.
As pair searching is an expensive operation, a generated pair list
is retained for a certain number of integration steps.
A buffer is needed to account for relative displacements
of atoms over the steps where a fixed pair list is retained.
Gromacs uses a buffered pair list by default. It also
uses clusters of particles, but these are not static as in the old charge group
scheme. Rather, the clusters are defined spatially and consist of 4 or 8
particles, which is convenient for stream computing, using e.g. SSE, AVX
or CUDA on GPUs. At neighbor search steps, a pair list is created with a
Verlet buffer, i.e. the pairlist cutoff is larger than the interaction
cutoff. In the nonbonded kernels, interactions are only computed when
a particle pair is within the cutoff distance at that particular time
step. This ensures that as particles move between pair search steps,
forces between nearly all particles within the cutoff distance are
calculated. We say *nearly* all particles, because Gromacs uses a fixed
pair list update frequency for efficiency. A particlepair, whose
distance was outside the cutoff, could possibly move enough during this
fixed number of steps that its distance is now within the cutoff. This
small chance results in a small energy drift, and the size of the chance
depends on the temperature. When temperature coupling is used, the
buffer size can be determined automatically, given a certain tolerance
on the energy drift. The default tolerance is 0.005 kJ/mol/ns per
particle, but in practice the energy drift is usually an order of
magnitude smaller. Note that in single precision for normal atomistic
simulations constraints cause a drift somewhere around 0.0001 kJ/mol/ns
per particle, so it doesn't make sense to go much lower than that.
The pair list is implemented in a very efficient fashion
based on clusters of particles. The simplest example is a cluster size
of 4 particles. The pair list is then constructed based on cluster
pairs. The clusterpair search is much faster searching based on
particle pairs, because :math:`4 \times 4 = 16` particle pairs are put
in the list at once. The nonbonded force calculation kernel can then
calculate many particlepair interactions at once, which maps nicely to
SIMD or SIMT units on modern hardware, which can perform multiple
floating operations at once. These nonbonded kernels are much faster
than the kernels used in the group scheme for most types of systems,
particularly on newer hardware. For further information on algorithmic
and implementation details of the Verlet cutoff scheme and the NxM
kernels, as well as detailed performance analysis, please consult the
following article: :ref:`182 `.
Additionally, when the list buffer is determined automatically as
described below, we also apply dynamic pair list pruning. The pair list
can be constructed infrequently, but that can lead to a lot of pairs in
the list that are outside the cutoff range for all or most of the life
time of this pair list. Such pairs can be pruned out by applying a
clusterpair kernel that only determines which clusters are in range.
Because of the way the nonbonded data is regularized in Gromacs, this
kernel is an order of magnitude faster than the search and the
interaction kernel. On the GPU this pruning is overlapped with the
integration on the CPU, so it is free in most cases. Therefore we can
prune every 410 integration steps with little overhead and
significantly reduce the number of cluster pairs in the interaction
kernel. This procedure is applied automatically, unless the user set the
pairlist buffer size manually.
Energy drift and pairlist buffering
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
For a canonical (NVT) ensemble, the average energy error caused by
diffusion of :math:`j` particles from outside the pairlist cutoff
:math:`r_\ell` to inside the interaction cutoff :math:`r_c` over the
lifetime of the list can be determined from the atomic displacements and
the shape of the potential at the cutoff. The displacement distribution
along one dimension for a freely moving particle with mass :math:`m`
over time :math:`t` at temperature :math:`T` is a Gaussian :math:`G(x)`
of zero mean and variance :math:`\sigma^2 = t^2 k_B T/m`. For the
distance between two particles, the variance changes to
:math:`\sigma^2 = \sigma_{12}^2 =
t^2 k_B T(1/m_1+1/m_2)`. Note that in practice particles usually
interact with (bump into) other particles over time :math:`t` and
therefore the real displacement distribution is much narrower. Given a
nonbonded interaction cutoff distance of :math:`r_c` and a pairlist
cutoff :math:`r_\ell=r_c+r_b` for :math:`r_b` the Verlet buffer size,
we can then write the average energy error after time :math:`t` for all
missing pair interactions between a single :math:`i` particle of type 1
surrounded by all :math:`j` particles that are of type 2 with number
density :math:`\rho_2`, when the interparticle distance changes from
:math:`r_0` to :math:`r_t`, as:
.. math:: \langle \Delta V \rangle =
\int_{0}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 V(r_t) G\!\left(\frac{r_tr_0}{\sigma}\right) d r_0\, d r_t
:label: eqnverletbufenergy
To evaluate this analytically, we need to make some approximations.
First we replace :math:`V(r_t)` by a Taylor expansion around
:math:`r_c`, then we can move the lower bound of the integral over
:math:`r_0` to :math:`\infty` which will simplify the result:
.. math:: \begin{aligned}
\langle \Delta V \rangle &\approx&
\int_{\infty}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 \Big[ V'(r_c) (r_t  r_c) +
\nonumber\\
& &
\phantom{\int_{\infty}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 \Big[}
V''(r_c)\frac{1}{2}(r_t  r_c)^2 +
\nonumber\\
& &
\phantom{\int_{\infty}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 \Big[}
V'''(r_c)\frac{1}{6}(r_t  r_c)^3 +
\nonumber\\
& &
\phantom{\int_{\infty}^{r_c} \int_{r_\ell}^\infty 4 \pi r_0^2 \rho_2 \Big[}
O \! \left((r_t  r_c)^4 \right)\Big] G\!\left(\frac{r_tr_0}{\sigma}\right) d r_0 \, d r_t\end{aligned}
:label: eqnverletaylor
Replacing the factor :math:`r_0^2` by :math:`(r_\ell + \sigma)^2`,
which results in a slight overestimate, allows us to calculate the
integrals analytically:
.. math:: \begin{aligned}
\langle \Delta V \rangle \!
&\approx&
4 \pi (r_\ell+\sigma)^2 \rho_2
\int_{\infty}^{r_c} \int_{r_\ell}^\infty \Big[ V'(r_c) (r_t  r_c) +
\nonumber\\
& &
\phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \int_{\infty}^{r_c} \int_{r_\ell}^\infty \Big[}
V''(r_c)\frac{1}{2}(r_t  r_c)^2 +
\nonumber\\
& &
\phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \int_{\infty}^{r_c} \int_{r_\ell}^\infty \Big[}
V'''(r_c)\frac{1}{6}(r_t  r_c)^3 \Big] G\!\left(\frac{r_tr_0}{\sigma}\right)
d r_0 \, d r_t\\
&=&
4 \pi (r_\ell+\sigma)^2 \rho_2 \bigg\{
\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] +
\nonumber\\
& &
\phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \bigg\{ }
\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] +
\nonumber\\
& &
\phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \bigg\{ }
\frac{1}{24}V'''(r_c)\bigg[ r_b\sigma(r_b^2+5\sigma^2) G\!\left(\frac{r_b}{\sigma}\right)
\nonumber\\
& &
\phantom{4 \pi (r_\ell+\sigma)^2 \rho_2 \bigg\{ \frac{1}{24}V'''(r_c)\bigg[ }
 (r_b^4+6r_b^2\sigma^2+3\sigma^4 ) E\!\left(\frac{r_b}{\sigma}\right) \bigg]
\bigg\}\end{aligned}
:label: eqnverletanalytical
where :math:`G(x)` is a Gaussian distribution with 0 mean and unit
variance and :math:`E(x)=\frac{1}{2}\mathrm{erfc}(x/\sqrt{2})`. We
always want to achieve small energy error, so :math:`\sigma` will be
small compared to both :math:`r_c` and :math:`r_\ell`, thus the
approximations in the equations above are good, since the Gaussian
distribution decays rapidly. The energy error needs to be averaged over
all particle pair types and weighted with the particle counts. In
Gromacs we don’t allow cancellation of error between pair types, so we
average the absolute values. To obtain the average energy error per unit
time, it needs to be divided by the neighborlist life time
:math:`t = ({\tt nstlist}  1)\times{\tt dt}`. The function can not be
inverted analytically, so we use bisection to obtain the buffer size
:math:`r_b` for a target drift. Again we note that in practice the error
we usually be much smaller than this estimate, as in the condensed phase
particle displacements will be much smaller than for freely moving
particles, which is the assumption used here.
When (bond) constraints are present, some particles will have fewer
degrees of freedom. This will reduce the energy errors. For simplicity,
we only consider one constraint per particle, the heaviest particle in
case a particle is involved in multiple constraints. This simplification
overestimates the displacement. The motion of a constrained particle is
a superposition of the 3D motion of the center of mass of both particles
and a 2D rotation around the center of mass. The displacement in an
arbitrary direction of a particle with 2 degrees of freedom is not
Gaussian, but rather follows the complementary error function:
.. math:: \frac{\sqrt{\pi}}{2\sqrt{2}\sigma}\,\mathrm{erfc}\left(\frac{r}{\sqrt{2}\,\sigma}\right)
:label: eqn2Ddisp
where :math:`\sigma^2` is again :math:`t^2 k_B T/m`. This distribution
can no longer be integrated analytically to obtain the energy error. But
we can generate a tight upper bound using a scaled and shifted Gaussian
distribution (not shown). This Gaussian distribution can then be used to
calculate the energy error as described above. The rotation displacement
around the center of mass can not be more than the length of the arm. To
take this into account, we scale :math:`\sigma` in
:eq:`eqn. %s ` (details not presented here) to
obtain an overestimate of the real displacement. This latter effect
significantly reduces the buffer size for longer neighborlist lifetimes
in e.g. water, as constrained hydrogens are by far the fastest
particles, but they can not move further than 0.1 nm from the heavy atom
they are connected to.
There is one important implementation detail that reduces the energy
errors caused by the finite Verlet buffer list size. The derivation
above assumes a particle pairlist. However, the Gromacs implementation
uses a cluster pairlist for efficiency. The pair list consists of pairs
of clusters of 4 particles in most cases, also called a
:math:`4 \times 4` list, but the list can also be :math:`4 \times 8`
(GPU CUDA kernels and AVX 256bit single precision kernels) or
:math:`4 \times 2` (SSE doubleprecision kernels). This means that the
pairlist is effectively much larger than the corresponding
:math:`1 \times 1` list. Thus slightly beyond the pairlist cutoff
there will still be a large fraction of particle pairs present in the
list. This fraction can be determined in a simulation and accurately
estimated under some reasonable assumptions. The fraction decreases with
increasing pairlist range, meaning that a smaller buffer can be used.
For typical allatom simulations with a cutoff of 0.9 nm this fraction
is around 0.9, which gives a reduction in the energy errors of a factor
of 10. This reduction is taken into account during the automatic Verlet
buffer calculation and results in a smaller buffer size.
.. _figverletdrift:
.. figure:: plots/verletdrift.*
:width: 9.00000cm
Energy drift per atom for an SPC/E water system at 300K with a
time step of 2 fs and a pairlist update period of 10 steps
(pairlist life time: 18 fs). PME was used with
``ewaldrtol`` set to 10\ :math:`^{5}`; this parameter
affects the shape of the potential at the cutoff. Error estimates
due to finite Verlet buffer size are shown for a :math:`1 \times 1`
atom pair list and :math:`4 \times 4` atom pair list without and with
(dashed line) cancellation of positive and negative errors. Real
energy drift is shown for simulations using double and
mixedprecision settings. Rounding errors in the SETTLE constraint
algorithm from the use of single precision causes the drift to become
negative at large buffer size. Note that at zero buffer size, the
real drift is small because positive (HH) and negative (OH) energy
errors cancel.
In :numref:`Fig. %s ` one can see that for small
buffer sizes the drift of the total energy is much smaller than the pair
energy error tolerance, due to cancellation of errors. For larger buffer
size, the error estimate is a factor of 6 higher than drift of the total
energy, or alternatively the buffer estimate is 0.024 nm too large. This
is because the protons don’t move freely over 18 fs, but rather vibrate.
Cutoff artifacts and switched interactions
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
By default, the pair potentials are shifted to be zero at the cutoff,
which makes the potential the integral of the force. However, there can
still be energy drift when the forces are nonzero at the cutoff. This
effect is extremely small and often not noticeable, as other integration
errors (e.g. from constraints) may dominate. To completely avoid cutoff
artifacts, the nonbonded forces can be switched exactly to zero at some
distance smaller than the neighbor list cutoff (there are several ways
to do this in Gromacs, see sec. :ref:`modnbint`). One then has a
buffer with the size equal to the neighbor list cutoff less the longest
interaction cutoff.
Simple search
^^^^^^^^^^^^^
Due to :eq:`eqns. %s ` and
:eq:`%s `, the vector
:math:`{\mathbf{r}_{ij}}` connecting images within the
cutoff :math:`R_c` can be found by constructing:
.. math:: \begin{aligned}
\mathbf{r}''' & = & \mathbf{r}_j\mathbf{r}_i \\
\mathbf{r}'' & = & \mathbf{r}'''  \mathbf{c}*\mathrm{round}(r'''_z/c_z) \\
\mathbf{r}' & = & \mathbf{r}''  \mathbf{b}*\mathrm{round}(r''_y/b_y) \\
\mathbf{r}_{ij} & = & \mathbf{r}'  \mathbf{a}*\mathrm{round}(r'_x/a_x)
\end{aligned}
:label: eqnsearchvec
When distances between two particles in a triclinic box are needed that
do not obey :eq:`eqn. %s `, many shifts of
combinations of box vectors need to be considered to find the nearest
image.
.. _figgrid:
.. figure:: plots/nstric.*
:width: 8.00000cm
Grid search in two dimensions. The arrows are the box vectors.
Grid search
^^^^^^^^^^^
The grid search is schematically depicted in
:numref:`Fig. %s `. All particles are put on the NS grid,
with the smallest spacing :math:`\ge` :math:`R_c/2` in each of the
directions. In the direction of each box vector, a particle :math:`i`
has three images. For each direction the image may be 1,0 or 1,
corresponding to a translation over 1, 0 or +1 box vector. We do not
search the surrounding NS grid cells for neighbors of :math:`i` and then
calculate the image, but rather construct the images first and then
search neighbors corresponding to that image of :math:`i`. As
:numref:`Fig. %s ` shows, some grid cells may be searched
more than once for different images of :math:`i`. This is not a problem,
since, due to the minimum image convention, at most one image will “see”
the :math:`j`particle. For every particle, fewer than 125 (5\ :math:`^3`)
neighboring cells are searched. Therefore, the algorithm scales linearly
with the number of particles. Although the prefactor is large, the
scaling behavior makes the algorithm far superior over the standard
:math:`O(N^2)` algorithm when there are more than a few hundred
particles. The grid search is equally fast for rectangular and triclinic
boxes. Thus for most protein and peptide simulations the rhombic
dodecahedron will be the preferred box shape.
.. _chargegroup:
Charge groups
^^^^^^^^^^^^^
Charge groups were originally introduced to reduce cutoff artifacts of
Coulomb interactions. This concept has been superseded by exact atomistic
cutoff treatments. For historical reasons charge groups are still
defined in the atoms section for each moleculetype in the topology,
but they are no longer used.
.. _forces:
Compute forces
~~~~~~~~~~~~~~
Potential energy
^^^^^^^^^^^^^^^^
When forces are computed, the potential energy of each interaction term
is computed as well. The total potential energy is summed for various
contributions, such as LennardJones, Coulomb, and bonded terms. It is
also possible to compute these contributions for *energymonitor groups*
of atoms that are separately defined (see sec. :ref:`groupconcept`).
Kinetic energy and temperature
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The temperature is given by the total kinetic energy of the
:math:`N`particle system:
.. math:: E_{kin} = {\frac{1}{2}}\sum_{i=1}^N m_i v_i^2
:label: eqntempEkin
From this the absolute temperature :math:`T` can be computed using:
.. math:: {\frac{1}{2}}N_{\mathrm{df}} kT = E_{\mathrm{kin}}
:label: eqnET
where :math:`k` is Boltzmann’s constant and :math:`N_{df}` is the
number of degrees of freedom which can be computed from:
.. math:: N_{\mathrm{df}} ~=~ 3 N  N_c  N_{\mathrm{com}}
:label: eqndofcoupling
Here :math:`N_c` is the number of *constraints* imposed on the system.
When performing molecular dynamics :math:`N_{\mathrm{com}}=3` additional
degrees of freedom must be removed, because the three centerofmass
velocities are constants of the motion, which are usually set to zero.
When simulating in vacuo, the rotation around the center of mass can
also be removed, in this case :math:`N_{\mathrm{com}}=6`. When more than
one temperaturecoupling group is used, the number of degrees of freedom
for group :math:`i` is:
.. math:: N^i_{\mathrm{df}} ~=~ (3 N^i  N^i_c) \frac{3 N  N_c  N_{\mathrm{com}}}{3 N  N_c}
:label: eqndofonecouplgroup
The kinetic energy can also be written as a tensor, which is necessary
for pressure calculation in a triclinic system, or systems where shear
forces are imposed:
.. math:: {\bf E}_{\mathrm{kin}} = {\frac{1}{2}}\sum_i^N m_i {\mathbf{v}_i}\otimes {\mathbf{v}_i}
:label: eqnEkintensor
Pressure and virial
^^^^^^^^^^^^^^^^^^^
The pressure tensor **P** is calculated from the difference between
kinetic energy :math:`E_{\mathrm{kin}}` and the virial
:math:`{\bf \Xi}`:
.. math:: {\bf P} = \frac{2}{V} ({\bf E}_{\mathrm{kin}}{\bf \Xi})
:label: eqnP
where :math:`V` is the volume of the computational box. The scalar
pressure :math:`P`, which can be used for pressure coupling in the case
of isotropic systems, is computed as:
.. math:: P = {\rm trace}({\bf P})/3
The virial :math:`{\bf \Xi}` tensor is defined as:
.. math:: {\bf \Xi} = {\frac{1}{2}}\sum_{i` for the integration of the
equations of motion. When extremely accurate integration with
temperature and/or pressure coupling is required, the velocity Verlet
integrators are also present and may be preferable (see
:ref:`vverlet`). The leapfrog algorithm uses positions
:math:`\mathbf{r}` at time :math:`t` and velocities
:math:`\mathbf{v}` at time
:math:`t{{\frac{1}{2}}{{\Delta t}}}`; it updates positions and
velocities using the forces :math:`\mathbf{F}(t)`
determined by the positions at time :math:`t` using these relations:
.. math:: \begin{aligned}
\mathbf{v}(t+{{\frac{1}{2}}{{\Delta t}}}) &~=~& \mathbf{v}(t{{\frac{1}{2}}{{\Delta t}}})+\frac{{{\Delta t}}}{m}\mathbf{F}(t) \\
\mathbf{r}(t+{{\Delta t}}) &~=~& \mathbf{r}(t)+{{\Delta t}}\mathbf{v}(t+{{\frac{1}{2}}{{\Delta t}}})\end{aligned}
:label: eqnleapfrogv
The algorithm is visualized in :numref:`Fig. %s `. It
produces trajectories that are identical to the Verlet \ :ref:`23 `
algorithm, whose positionupdate relation is
.. math:: \mathbf{r}(t+{{\Delta t}})~=~2\mathbf{r}(t)  \mathbf{r}(t{{\Delta t}}) + \frac{1}{m}\mathbf{F}(t){{\Delta t}}^2+O({{\Delta t}}^4)
:label: eqnleapfrogp
The algorithm is of third order in :math:`\mathbf{r}` and
is timereversible. See ref. \ :ref:`24 ` for the
merits of this algorithm and comparison with other time integration
algorithms.
The equations of motion are modified for temperature coupling and
pressure coupling, and extended to include the conservation of
constraints, all of which are described below.
.. _vverlet:
The velocity Verlet integrator
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The velocity Verlet algorithm\ :ref:`25 ` is also implemented in
Gromacs, though it is not yet fully integrated with all sets of options.
In velocity Verlet, positions :math:`\mathbf{r}` and
velocities :math:`\mathbf{v}` at time :math:`t` are used
to integrate the equations of motion; velocities at the previous half
step are not required.
.. math:: \begin{aligned}
\mathbf{v}(t+{{\frac{1}{2}}{{\Delta t}}}) &~=~& \mathbf{v}(t)+\frac{{{\Delta t}}}{2m}\mathbf{F}(t) \\
\mathbf{r}(t+{{\Delta t}}) &~=~& \mathbf{r}(t)+{{\Delta t}}\,\mathbf{v}(t+{{\frac{1}{2}}{{\Delta t}}}) \\
\mathbf{v}(t+{{\Delta t}}) &~=~& \mathbf{v}(t+{{\frac{1}{2}}{{\Delta t}}})+\frac{{{\Delta t}}}{2m}\mathbf{F}(t+{{\Delta t}})\end{aligned}
:label: eqnvelocityverlet1
or, equivalently,
.. math:: \begin{aligned}
\mathbf{r}(t+{{\Delta t}}) &~=~& \mathbf{r}(t)+ {{\Delta t}}\,\mathbf{v} + \frac{{{\Delta t}}^2}{2m}\mathbf{F}(t) \\
\mathbf{v}(t+{{\Delta t}}) &~=~& \mathbf{v}(t)+ \frac{{{\Delta t}}}{2m}\left[\mathbf{F}(t) + \mathbf{F}(t+{{\Delta t}})\right]\end{aligned}
:label: eqnvelocityverlet2
With no temperature or pressure coupling, and with *corresponding*
starting points, leapfrog and velocity Verlet will generate identical
trajectories, as can easily be verified by hand from the equations
above. Given a single starting file with the *same* starting point
:math:`\mathbf{x}(0)` and
:math:`\mathbf{v}(0)`, leapfrog and velocity Verlet will
*not* give identical trajectories, as leapfrog will interpret the
velocities as corresponding to :math:`t={{\frac{1}{2}}{{\Delta t}}}`,
while velocity Verlet will interpret them as corresponding to the
timepoint :math:`t=0`.
Understanding reversible integrators: The Trotter decomposition
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
To further understand the relationship between velocity Verlet and
leapfrog integration, we introduce the reversible Trotter formulation
of dynamics, which is also useful to understanding implementations of
thermostats and barostats in Gromacs.
A system of coupled, firstorder differential equations can be evolved
from time :math:`t = 0` to time :math:`t` by applying the evolution
operator
.. math:: \begin{aligned}
\Gamma(t) &=& \exp(iLt) \Gamma(0) \nonumber \\
iL &=& \dot{\Gamma}\cdot \nabla_{\Gamma},\end{aligned}
:label: eqnevoluoperator
where :math:`L` is the Liouville operator, and :math:`\Gamma` is the
multidimensional vector of independent variables (positions and
velocities). A shorttime approximation to the true operator, accurate
at time :math:`{{\Delta t}}= t/P`, is applied :math:`P` times in
succession to evolve the system as
.. math:: \Gamma(t) = \prod_{i=1}^P \exp(iL{{\Delta t}}) \Gamma(0)
:label: eqnevolvesystem
For NVE dynamics, the Liouville operator is
.. math:: \begin{aligned}
iL = \sum_{i=1}^{N} {{\mathbf{v}}}_i \cdot \nabla_{{{\mathbf{r}}}_i} + \sum_{i=1}^N \frac{1}{m_i}{{\mathbf{F}}}(r_i) \cdot \nabla_{{{\mathbf{v}}}_i}.\end{aligned}
:label: eqnliouvilleoperator
This can be split into two additive operators
.. math:: \begin{aligned}
iL_1 &=& \sum_{i=1}^N \frac{1}{m_i}{{\mathbf{F}}}(r_i) \cdot \nabla_{{{\mathbf{v}}}_i} \nonumber \\
iL_2 &=& \sum_{i=1}^{N} {{\mathbf{v}}}_i \cdot \nabla_{{{\mathbf{r}}}_i} \end{aligned}
:label: eqnlotwoadditive
Then a shorttime, symmetric, and thus reversible approximation of the
true dynamics will be
.. math:: \begin{aligned}
\exp(iL{{\Delta t}}) = \exp(iL_2{{\frac{1}{2}}{{\Delta t}}}) \exp(iL_1{{\Delta t}}) \exp(iL_2{{\frac{1}{2}}{{\Delta t}}}) + \mathcal{O}({{\Delta t}}^3).
\end{aligned}
:label: eqNVETrotter
This corresponds to velocity Verlet integration. The first exponential
term over :math:`{{\frac{1}{2}}{{\Delta t}}}` corresponds to a velocity
halfstep, the second exponential term over :math:`{{\Delta t}}`
corresponds to a full velocity step, and the last exponential term over
:math:`{{\frac{1}{2}}{{\Delta t}}}` is the final velocity half step. For
future times :math:`t = n{{\Delta t}}`, this becomes
.. math:: \begin{aligned}
\exp(iLn{{\Delta t}}) &\approx& \left(\exp(iL_2{{\frac{1}{2}}{{\Delta t}}}) \exp(iL_1{{\Delta t}}) \exp(iL_2{{\frac{1}{2}}{{\Delta t}}})\right)^n \nonumber \\
&\approx& \exp(iL_2{{\frac{1}{2}}{{\Delta t}}}) \bigg(\exp(iL_1{{\Delta t}}) \exp(iL_2{{\Delta t}})\bigg)^{n1} \nonumber \\
& & \;\;\;\; \exp(iL_1{{\Delta t}}) \exp(iL_2{{\frac{1}{2}}{{\Delta t}}}) \end{aligned}
:label: eqntrottertimestep
This formalism allows us to easily see the difference between the
different flavors of Verlet integrators. The leapfrog integrator can be
seen as starting with :eq:`Eq. %s ` with the
:math:`\exp\left(iL_1 {\Delta t}\right)` term, instead of the halfstep
velocity term, yielding
.. math:: \begin{aligned}
\exp(iLn{\Delta t}) &=& \exp\left(iL_1 {\Delta t}\right) \exp\left(iL_2 {{\Delta t}}\right) + \mathcal{O}({{\Delta t}}^3).\end{aligned}
:label: eqnleapfroghalfvel
Here, the full step in velocity is between
:math:`t{{\frac{1}{2}}{{\Delta t}}}` and
:math:`t+{{\frac{1}{2}}{{\Delta t}}}`, since it is a combination of the
velocity half steps in velocity Verlet. For future times
:math:`t = n{{\Delta t}}`, this becomes
.. math:: \begin{aligned}
\exp(iLn{\Delta t}) &\approx& \bigg(\exp\left(iL_1 {\Delta t}\right) \exp\left(iL_2 {{\Delta t}}\right) \bigg)^{n}.\end{aligned}
:label: eqnvelverlethalfvel
Although at first this does not appear symmetric, as long as the full
velocity step is between :math:`t{{\frac{1}{2}}{{\Delta t}}}` and
:math:`t+{{\frac{1}{2}}{{\Delta t}}}`, then this is simply a way of
starting velocity Verlet at a different place in the cycle.
Even though the trajectory and thus potential energies are identical
between leapfrog and velocity Verlet, the kinetic energy and
temperature will not necessarily be the same. Standard velocity Verlet
uses the velocities at the :math:`t` to calculate the kinetic energy and
thus the temperature only at time :math:`t`; the kinetic energy is then
a sum over all particles
.. math:: \begin{aligned}
KE_{\mathrm{full}}(t) &=& \sum_i \left(\frac{1}{2m_i}\mathbf{v}_i(t)\right)^2 \nonumber\\
&=& \sum_i \frac{1}{2m_i}\left(\frac{1}{2}\mathbf{v}_i(t{{\frac{1}{2}}{{\Delta t}}})+\frac{1}{2}\mathbf{v}_i(t+{{\frac{1}{2}}{{\Delta t}}})\right)^2,\end{aligned}
:label: eqnTrotterEkin
with the square on the *outside* of the average. Standard leapfrog
calculates the kinetic energy at time :math:`t` based on the average
kinetic energies at the timesteps :math:`t+{{\frac{1}{2}}{{\Delta t}}}`
and :math:`t{{\frac{1}{2}}{{\Delta t}}}`, or the sum over all particles
.. math:: \begin{aligned}
KE_{\mathrm{average}}(t) &=& \sum_i \frac{1}{2m_i}\left(\frac{1}{2}\mathbf{v}_i(t{{\frac{1}{2}}{{\Delta t}}})^2+\frac{1}{2}\mathbf{v}_i(t+{{\frac{1}{2}}{{\Delta t}}})^2\right),\end{aligned}
:label: eqnTrottersumparticles
where the square is *inside* the average.
A nonstandard variant of velocity Verlet which averages the kinetic
energies :math:`KE(t+{{\frac{1}{2}}{{\Delta t}}})` and
:math:`KE(t{{\frac{1}{2}}{{\Delta t}}})`, exactly like leapfrog, is
also now implemented in Gromacs (as :ref:`mdp` file option
:mdpvalue:`integrator=mdvvavek`). Without temperature and pressure coupling,
velocity Verlet with halfstepaveraged kinetic energies and leapfrog
will be identical up to numerical precision. For temperature and
pressurecontrol schemes, however, velocity Verlet with
halfstepaveraged kinetic energies and leapfrog will be different, as
will be discussed in the section in thermostats and barostats.
The halfstepaveraged kinetic energy and temperature are slightly more
accurate for a given step size; the difference in average kinetic
energies using the halfstepaveraged kinetic energies (
:mdpvalue:`integrator=md` and :mdpvalue:`integrator=mdvvavek`
) will be closer to the kinetic energy obtained in the limit
of small step size than will the fullstep kinetic energy (using
:mdpvalue:`integrator=mdvv`). For NVE simulations, this difference is usually not
significant, since the positions and velocities of the particles are
still identical; it makes a difference in the way the temperature of
the simulations are **interpreted**, but **not** in the trajectories that
are produced. Although the kinetic energy is more accurate with the
halfstepaveraged method, meaning that it changes less as the timestep
gets large, it is also more noisy. The RMS deviation of the total energy
of the system (sum of kinetic plus potential) in the halfstepaveraged
kinetic energy case will be higher (about twice as high in most cases)
than the fullstep kinetic energy. The drift will still be the same,
however, as again, the trajectories are identical.
For NVT simulations, however, there **will** be a difference, as discussed
in the section on temperature control, since the velocities of the
particles are adjusted such that kinetic energies of the simulations,
which can be calculated either way, reach the distribution corresponding
to the set temperature. In this case, the three methods will not give
identical results.
Because the velocity and position are both defined at the same time
:math:`t` the velocity Verlet integrator can be used for some methods,
especially rigorously correct pressure control methods, that are not
actually possible with leapfrog. The integration itself takes
negligibly more time than leapfrog, but twice as many communication
calls are currently required. In most cases, and especially for large
systems where communication speed is important for parallelization and
differences between thermodynamic ensembles vanish in the :math:`1/N`
limit, and when only NVT ensembles are required, leapfrog will likely
be the preferred integrator. For pressure control simulations where the
fine details of the thermodynamics are important, only velocity Verlet
allows the true ensemble to be calculated. In either case, simulation
with double precision may be required to get fine details of
thermodynamics correct.
Multiple time stepping
~~~~~~~~~~~~~~~~~~~~~~
Several other simulation packages uses multiple time stepping for bonds
and/or the PME mesh forces. In Gromacs we have not implemented this
(yet), since we use a different philosophy. Bonds can be constrained
(which is also a more sound approximation of a physical quantum
oscillator), which allows the smallest time step to be increased to the
larger one. This not only halves the number of force calculations, but
also the update calculations. For even larger time steps, angle
vibrations involving hydrogen atoms can be removed using virtual
interaction sites (see sec. :ref:`rmfast`), which brings the shortest
time step up to PME mesh update frequency of a multiple time stepping
scheme.
.. _tempcoupling:
Temperature coupling
~~~~~~~~~~~~~~~~~~~~
While direct use of molecular dynamics gives rise to the NVE (constant
number, constant volume, constant energy ensemble), most quantities that
we wish to calculate are actually from a constant temperature (NVT)
ensemble, also called the canonical ensemble. Gromacs can use the
*weakcoupling* scheme of Berendsen \ :ref:`26 `, stochastic
randomization through the Andersen thermostat \ :ref:`27 `, the
extended ensemble NoséHoover scheme \ :ref:`28 `, :ref:`29 `, or a
velocityrescaling scheme \ :ref:`30 ` to
simulate constant temperature, with advantages of each of the schemes
laid out below.
There are several other reasons why it might be necessary to control the
temperature of the system (drift during equilibration, drift as a result
of force truncation and integration errors, heating due to external or
frictional forces), but this is not entirely correct to do from a
thermodynamic standpoint, and in some cases only masks the symptoms
(increase in temperature of the system) rather than the underlying
problem (deviations from correct physics in the dynamics). For larger
systems, errors in ensemble averages and structural properties incurred
by using temperature control to remove slow drifts in temperature appear
to be negligible, but no completely comprehensive comparisons have been
carried out, and some caution must be taking in interpreting the
results.
When using temperature and/or pressure coupling the total energy is no
longer conserved. Instead there is a conserved energy quantity the
formula of which will depend on the combination or temperature and
pressure coupling algorithm used. For all coupling algorithms, except
for Andersen temperature coupling and ParrinelloRahman pressure
coupling combined with shear stress, the conserved energy quantity is
computed and stored in the energy and log file. Note that this quantity
will not be conserved when external forces are applied to the system,
such as pulling on group with a changing distance or an electric field.
Furthermore, how well the energy is conserved depends on the accuracy of
all algorithms involved in the simulation. Usually the algorithms that
cause most drift are constraints and the pairlist buffer, depending on
the parameters used.
Berendsen temperature coupling
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The Berendsen algorithm mimics weak coupling with firstorder kinetics
to an external heat bath with given temperature :math:`T_0`. See
ref. \ :ref:`31 ` for a comparison with the NoséHoover scheme. The
effect of this algorithm is that a deviation of the system temperature
from :math:`T_0` is slowly corrected according to:
.. math:: \frac{{\mbox{d}}T}{{\mbox{d}}t} = \frac{T_0T}{\tau}
:label: eqnTcoupling
which means that a temperature deviation decays exponentially with a
time constant :math:`\tau`. This method of coupling has the advantage
that the strength of the coupling can be varied and adapted to the user
requirement: for equilibration purposes the coupling time can be taken
quite short (*e.g.* 0.01 ps), but for reliable equilibrium runs it can
be taken much longer (*e.g.* 0.5 ps) in which case it hardly influences
the conservative dynamics.
The Berendsen thermostat suppresses the fluctuations of the kinetic
energy. This means that one does not generate a proper canonical
ensemble, so rigorously, the sampling will be incorrect. This error
scales with :math:`1/N`, so for very large systems most ensemble
averages will not be affected significantly, except for the distribution
of the kinetic energy itself. However, fluctuation properties, such as
the heat capacity, will be affected. A similar thermostat which does
produce a correct ensemble is the velocity rescaling
thermostat \ :ref:`30 ` described below.
The heat flow into or out of the system is affected by scaling the
velocities of each particle every step, or every :math:`n_\mathrm{TC}`
steps, with a timedependent factor :math:`\lambda`, given by:
.. math:: \lambda = \left[ 1 + \frac{n_\mathrm{TC} \Delta t}{\tau_T}
\left\{\frac{T_0}{T(t  {{\frac{1}{2}}{{\Delta t}}})}  1 \right\} \right]^{1/2}
:label: eqnlambda
The parameter :math:`\tau_T` is close, but not exactly equal, to the
time constant :math:`\tau` of the temperature coupling
(:eq:`eqn. %s `):
.. math:: \tau = 2 C_V \tau_T / N_{df} k
:label: eqnTcoupltau
where :math:`C_V` is the total heat capacity of the system, :math:`k`
is Boltzmann’s constant, and :math:`N_{df}` is the total number of
degrees of freedom. The reason that :math:`\tau \neq \tau_T` is that the
kinetic energy change caused by scaling the velocities is partly
redistributed between kinetic and potential energy and hence the change
in temperature is less than the scaling energy. In practice, the ratio
:math:`\tau / \tau_T` ranges from 1 (gas) to 2 (harmonic solid) to 3
(water). When we use the term *temperature coupling time constant*, we
mean the parameter :math:`\tau_T`. **Note** that in practice the scaling
factor :math:`\lambda` is limited to the range of 0.8
:math:`<= \lambda <=` 1.25, to avoid scaling by very large numbers which
may crash the simulation. In normal use, :math:`\lambda` will always be
much closer to 1.0.
The thermostat modifies the kinetic energy at each scaling step by:
.. math:: \Delta E_k = (\lambda  1)^2 E_k
:label: eqnThermostat
The sum of these changes over the run needs to subtracted from the
total energy to obtain the conserved energy quantity.
Velocityrescaling temperature coupling
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The velocityrescaling thermostat \ :ref:`30 `
is essentially a Berendsen thermostat (see above) with an additional
stochastic term that ensures a correct kinetic energy distribution by
modifying it according to
.. math:: {\mbox{d}}K = (K_0  K) \frac{{\mbox{d}}t}{\tau_T} + 2 \sqrt{\frac{K K_0}{N_f}} \frac{{\mbox{d}}W}{\sqrt{\tau_T}},
:label: eqnvrescale
where :math:`K` is the kinetic energy, :math:`N_f` the number of
degrees of freedom and :math:`{\mbox{d}}W` a Wiener process. There are
no additional parameters, except for a random seed. This thermostat
produces a correct canonical ensemble and still has the advantage of the
Berendsen thermostat: first order decay of temperature deviations and no
oscillations.
Andersen thermostat
^^^^^^^^^^^^^^^^^^^
One simple way to maintain a thermostatted ensemble is to take an
:math:`NVE` integrator and periodically reselect the velocities of the
particles from a MaxwellBoltzmann distribution \ :ref:`27 `. This
can either be done by randomizing all the velocities simultaneously
(massive collision) every :math:`\tau_T/{{\Delta t}}` steps
(``andersenmassive``), or by randomizing every particle
with some small probability every timestep (``andersen``),
equal to :math:`{{\Delta t}}/\tau`, where in both cases
:math:`{{\Delta t}}` is the timestep and :math:`\tau_T` is a
characteristic coupling time scale. Because of the way constraints
operate, all particles in the same constraint group must be randomized
simultaneously. Because of parallelization issues, the
``andersen`` version cannot currently (5.0) be used in
systems with constraints. ``andersenmassive`` can be used
regardless of constraints. This thermostat is also currently only
possible with velocity Verlet algorithms, because it operates directly
on the velocities at each timestep.
This algorithm completely avoids some of the ergodicity issues of other
thermostatting algorithms, as energy cannot flow back and forth between
energetically decoupled components of the system as in velocity scaling
motions. However, it can slow down the kinetics of system by randomizing
correlated motions of the system, including slowing sampling when
:math:`\tau_T` is at moderate levels (less than 10 ps). This algorithm
should therefore generally not be used when examining kinetics or
transport properties of the system \ :ref:`32 `.
NoséHoover temperature coupling
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The Berendsen weakcoupling algorithm is extremely efficient for
relaxing a system to the target temperature, but once the system has
reached equilibrium it might be more important to probe a correct
canonical ensemble. This is unfortunately not the case for the
weakcoupling scheme.
To enable canonical ensemble simulations, Gromacs also supports the
extendedensemble approach first proposed by Nosé :ref:`28 ` and later
modified by Hoover \ :ref:`29 `. The system Hamiltonian is extended by
introducing a thermal reservoir and a friction term in the equations of
motion. The friction force is proportional to the product of each
particle’s velocity and a friction parameter, :math:`\xi`. This friction
parameter (or *heat bath* variable) is a fully dynamic quantity with its
own momentum (:math:`p_{\xi}`) and equation of motion; the time
derivative is calculated from the difference between the current kinetic
energy and the reference temperature.
In this formulation, the particles´ equations of motion in
the global :ref:`MD scheme ` are replaced by:
.. math:: \frac {{\mbox{d}}^2\mathbf{r}_i}{{\mbox{d}}t^2} = \frac{\mathbf{F}_i}{m_i} 
\frac{p_{\xi}}{Q}\frac{{\mbox{d}}\mathbf{r}_i}{{\mbox{d}}t} ,
:label: eqnNHeqnofmotion
where the equation of motion for the heat bath parameter :math:`\xi` is:
.. math:: \frac {{\mbox{d}}p_{\xi}}{{\mbox{d}}t} = \left( T  T_0 \right).
:label: eqnNHheatbath
The reference temperature is denoted :math:`T_0`, while :math:`T` is
the current instantaneous temperature of the system. The strength of the
coupling is determined by the constant :math:`Q` (usually called the
*mass parameter* of the reservoir) in combination with the reference
temperature. [1]_
The conserved quantity for the NoséHoover equations of motion is not
the total energy, but rather
.. math:: \begin{aligned}
H = \sum_{i=1}^{N} \frac{{{\mathbf{p}}}_i}{2m_i} + U\left({{\mathbf{r}}}_1,{{\mathbf{r}}}_2,\ldots,{{\mathbf{r}}}_N\right) +\frac{p_{\xi}^2}{2Q} + N_{f}kT\xi,\end{aligned}
:label: eqnNHconservedbasic
where :math:`N_f` is the total number of degrees of freedom.
In our opinion, the mass parameter is a somewhat awkward way of
describing coupling strength, especially due to its dependence on
reference temperature (and some implementations even include the number
of degrees of freedom in your system when defining :math:`Q`). To
maintain the coupling strength, one would have to change :math:`Q` in
proportion to the change in reference temperature. For this reason, we
prefer to let the Gromacs user work instead with the period
:math:`\tau_T` of the oscillations of kinetic energy between the system
and the reservoir instead. It is directly related to :math:`Q` and
:math:`T_0` via:
.. math:: Q = \frac {\tau_T^2 T_0}{4 \pi^2}.
:label: eqnNHQ
This provides a much more intuitive way of selecting the NoséHoover
coupling strength (similar to the weakcoupling relaxation), and in
addition :math:`\tau_T` is independent of system size and reference
temperature.
It is however important to keep the difference between the weakcoupling
scheme and the NoséHoover algorithm in mind: Using weak coupling you
get a strongly damped *exponential relaxation*, while the NoséHoover
approach produces an *oscillatory relaxation*. The actual time it takes
to relax with NoséHoover coupling is several times larger than the
period of the oscillations that you select. These oscillations (in
contrast to exponential relaxation) also means that the time constant
normally should be 4–5 times larger than the relaxation time used with
weak coupling, but your mileage may vary.
NoséHoover dynamics in simple systems such as collections of harmonic
oscillators, can be *nonergodic*, meaning that only a subsection of
phase space is ever sampled, even if the simulations were to run for
infinitely long. For this reason, the NoséHoover chain approach was
developed, where each of the NoséHoover thermostats has its own
NoséHoover thermostat controlling its temperature. In the limit of an
infinite chain of thermostats, the dynamics are guaranteed to be
ergodic. Using just a few chains can greatly improve the ergodicity, but
recent research has shown that the system will still be nonergodic, and
it is still not entirely clear what the practical effect of
this \ :ref:`33 `. Currently, the default number of chains is 10, but
this can be controlled by the user. In the case of chains, the equations
are modified in the following way to include a chain of thermostatting
particles \ :ref:`34 `:
.. math:: \begin{aligned}
\frac {{\mbox{d}}^2\mathbf{r}_i}{{\mbox{d}}t^2} &~=~& \frac{\mathbf{F}_i}{m_i}  \frac{p_{{\xi}_1}}{Q_1} \frac{{\mbox{d}}\mathbf{r}_i}{{\mbox{d}}t} \nonumber \\
\frac {{\mbox{d}}p_{{\xi}_1}}{{\mbox{d}}t} &~=~& \left( T  T_0 \right)  p_{{\xi}_1} \frac{p_{{\xi}_2}}{Q_2} \nonumber \\
\frac {{\mbox{d}}p_{{\xi}_{i=2\ldots N}}}{{\mbox{d}}t} &~=~& \left(\frac{p_{\xi_{i1}}^2}{Q_{i1}} kT\right)  p_{\xi_i} \frac{p_{\xi_{i+1}}}{Q_{i+1}} \nonumber \\
\frac {{\mbox{d}}p_{\xi_N}}{{\mbox{d}}t} &~=~& \left(\frac{p_{\xi_{N1}}^2}{Q_{N1}}kT\right)
\end{aligned}
:label: eqnNHchaineqnofmotion
The conserved quantity for NoséHoover chains is
.. math:: \begin{aligned}
H = \sum_{i=1}^{N} \frac{{{\mathbf{p}}}_i}{2m_i} + U\left({{\mathbf{r}}}_1,{{\mathbf{r}}}_2,\ldots,{{\mathbf{r}}}_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 \end{aligned}
:label: eqnNHconservedquantity
The values and velocities of the NoséHoover thermostat variables are
generally not included in the output, as they take up a fair amount of
space and are generally not important for analysis of simulations, but
by setting an :ref:`mdp` option the values of all the positions and velocities
of all NoséHoover particles in the chain are written to the :ref:`edr` file.
Leapfrog simulations currently can only have NoséHoover chain lengths
of 1, but this will likely be updated in later version.
As described in the integrator section, for temperature coupling, the
temperature that the algorithm attempts to match to the reference
temperature is calculated differently in velocity Verlet and leapfrog
dynamics. Velocity Verlet (*mdvv*) uses the fullstep kinetic energy,
while leapfrog and *mdvvavek* use the halfstepaveraged kinetic
energy.
We can examine the Trotter decomposition again to better understand the
differences between these constanttemperature integrators. In the case
of NoséHoover dynamics (for simplicity, using a chain with :math:`N=1`,
with more details in Ref. \ :ref:`35 `), we split the Liouville
operator as
.. math:: iL = iL_1 + iL_2 + iL_{\mathrm{NHC}},
:label: eqnNHTrotter
where
.. math:: \begin{aligned}
iL_1 &=& \sum_{i=1}^N \left[\frac{{{\mathbf{p}}}_i}{m_i}\right]\cdot \frac{\partial}{\partial {{\mathbf{r}}}_i} \nonumber \\
iL_2 &=& \sum_{i=1}^N {{\mathbf{F}}}_i\cdot \frac{\partial}{\partial {{\mathbf{p}}}_i} \nonumber \\
iL_{\mathrm{NHC}} &=& \sum_{i=1}^N\frac{p_{\xi}}{Q}{{\mathbf{v}}}_i\cdot \nabla_{{{\mathbf{v}}}_i} +\frac{p_{\xi}}{Q}\frac{\partial }{\partial \xi} + \left( T  T_0 \right)\frac{\partial }{\partial p_{\xi}}\end{aligned}
:label: eqnNHTrotter2
For standard velocity Verlet with NoséHoover temperature control, this
becomes
.. math:: \begin{aligned}
\exp(iL{\Delta t}) &=& \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \exp\left(iL_2 {\Delta t}/2\right) \nonumber \\
&&\exp\left(iL_1 {\Delta t}\right) \exp\left(iL_2 {\Delta t}/2\right) \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) + \mathcal{O}({{\Delta t}}^3).\end{aligned}
:label: eqnNHTrotter3
For halfstepaveraged temperature control using *mdvvavek*, this
decomposition will not work, since we do not have the full step
temperature until after the second velocity step. However, we can
construct an alternate decomposition that is still reversible, by
switching the place of the NHC and velocity portions of the
decomposition:
.. math:: \begin{aligned}
\exp(iL{\Delta t}) &=& \exp\left(iL_2 {\Delta t}/2\right) \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right)\exp\left(iL_1 {\Delta t}\right)\nonumber \\
&&\exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \exp\left(iL_2 {\Delta t}/2\right)+ \mathcal{O}({{\Delta t}}^3)
\end{aligned}
:label: eqhalfstepNHCintegrator
This formalism allows us to easily see the difference between the
different flavors of velocity Verlet integrator. The leapfrog
integrator can be seen as starting with
:eq:`Eq. %s ` just before the
:math:`\exp\left(iL_1 {\Delta t}\right)` term, yielding:
.. math:: \begin{aligned}
\exp(iL{\Delta t}) &=& \exp\left(iL_1 {\Delta t}\right) \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \nonumber \\
&&\exp\left(iL_2 {\Delta t}\right) \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) + \mathcal{O}({{\Delta t}}^3)\end{aligned}
:label: eqnNHleapfrog
and then using some algebra tricks to solve for some quantities are
required before they are actually calculated \ :ref:`36 `.
Group temperature coupling
^^^^^^^^^^^^^^^^^^^^^^^^^^
In Gromacs temperature coupling can be performed on groups of atoms,
typically a protein and solvent. The reason such algorithms were
introduced is that energy exchange between different components is not
perfect, due to different effects including cutoffs etc. If now the
whole system is coupled to one heat bath, water (which experiences the
largest cutoff noise) will tend to heat up and the protein will cool
down. Typically 100 K differences can be obtained. With the use of
proper electrostatic methods (PME) these difference are much smaller but
still not negligible. The parameters for temperature coupling in groups
are given in the :ref:`mdp` file. Recent investigation has shown that small
temperature differences between protein and water may actually be an
artifact of the way temperature is calculated when there are finite
timesteps, and very large differences in temperature are likely a sign
of something else seriously going wrong with the system, and should be
investigated carefully \ :ref:`37 `.
One special case should be mentioned: it is possible to
temperaturecouple only part of the system, leaving other parts without
temperature coupling. This is done by specifying :math:`{1}` for the
time constant :math:`\tau_T` for the group that should not be
thermostatted. If only part of the system is thermostatted, the system
will still eventually converge to an NVT system. In fact, one suggestion
for minimizing errors in the temperature caused by discretized timesteps
is that if constraints on the water are used, then only the water
degrees of freedom should be thermostatted, not protein degrees of
freedom, as the higher frequency modes in the protein can cause larger
deviations from the *true* temperature, the temperature obtained with
small timesteps \ :ref:`37 `.
Pressure coupling
~~~~~~~~~~~~~~~~~
In the same spirit as the temperature coupling, the system can also be
coupled to a *pressure bath.* Gromacs supports both the Berendsen
algorithm \ :ref:`26 ` that scales coordinates and box
vectors every step, the extendedensemble ParrinelloRahman
approach \ :ref:`38 `, :ref:`39 `, and for the
velocity Verlet variants, the MartynaTuckermanTobiasKlein (MTTK)
implementation of pressure control \ :ref:`35 `.
ParrinelloRahman and Berendsen can be combined with any of the
temperature coupling methods above. MTTK can only be used with
NoséHoover temperature control. From 5.1 afterwards, it can only used
when the system does not have constraints.
Berendsen pressure coupling
^^^^^^^^^^^^^^^^^^^^^^^^^^^
The Berendsen algorithm rescales the coordinates and box vectors every
step, or every :math:`n_\mathrm{PC}` steps, with a matrix :math:`\mu`,
which has the effect of a firstorder kinetic relaxation of the pressure
towards a given reference pressure :math:`{\bf P}_0` according to
.. math:: \frac{{\mbox{d}}{\bf P}}{{\mbox{d}}t} = \frac{{\bf P}_0{\bf P}}{\tau_p}.
:label: eqnberendsenpressure
The scaling matrix :math:`\mu` is given by
.. math:: \mu_{ij}
= \delta_{ij}  \frac{n_\mathrm{PC}\Delta t}{3\, \tau_p} \beta_{ij} \{P_{0ij}  P_{ij}(t) \}.
:label: eqnmu
Here, :math:`\beta` is the isothermal compressibility of the system. In
most cases this will be a diagonal matrix, with equal elements on the
diagonal, the value of which is generally not known. It suffices to take
a rough estimate because the value of :math:`\beta` only influences the
noncritical time constant of the pressure relaxation without affecting
the average pressure itself. For water at 1 atm and 300 K
:math:`\beta = 4.6 \times 10^{10}`
Pa\ :math:`^{1} = 4.6 \times 10^{5}` bar\ :math:`^{1}`, which is
:math:`7.6 \times 10^{4}` MD units (see chapter :ref:`defunits`). Most
other liquids have similar values. When scaling completely
anisotropically, the system has to be rotated in order to obey
:eq:`eqn. %s `. This rotation is approximated in first order in the
scaling, which is usually less than :math:`10^{4}`. The actual scaling
matrix :math:`\mu'` is
.. math:: \mathbf{\mu'} =
\left(\begin{array}{ccc}
\mu_{xx} & \mu_{xy} + \mu_{yx} & \mu_{xz} + \mu_{zx} \\
0 & \mu_{yy} & \mu_{yz} + \mu_{zy} \\
0 & 0 & \mu_{zz}
\end{array}\right).
:label: eqnberendsenpressurescaling
The velocities are neither scaled nor rotated. Since the equations of
motion are modified by pressure coupling, the conserved energy quantity
also needs to be modified. For first order pressure coupling, the work
the barostat applies to the system every step needs to be subtracted
from the total energy to obtain the conserved energy quantity:
.. math::  \sum_{i,j} (\mu_{ij} \delta_{ij}) P_{ij} V =
\sum_{i,j} 2(\mu_{ij} \delta_{ij}) \Xi_{ij}
:label: eqnberendsenpressureconserved
where :math:`\delta_{ij}` is the Kronecker delta and :math:`{\bf \Xi}`
is the virial. Note that the factor 2 originates from the factor
:math:`\frac{1}{2}` in the virial definition
(:eq:`eqn. %s `).
In Gromacs, the Berendsen scaling can also be done isotropically, which
means that instead of :math:`\mathbf{P}` a diagonal matrix
with elements of size trace\ :math:`(\mathbf{P})/3` is
used. For systems with interfaces, semiisotropic scaling can be useful.
In this case, the :math:`x/y`directions are scaled isotropically and
the :math:`z` direction is scaled independently. The compressibility in
the :math:`x/y` or :math:`z`direction can be set to zero, to scale only
in the other direction(s).
If you allow full anisotropic deformations and use constraints you might
have to scale more slowly or decrease your timestep to avoid errors from
the constraint algorithms. It is important to note that although the
Berendsen pressure control algorithm yields a simulation with the
correct average pressure, it does not yield the exact NPT ensemble, and
it is not yet clear exactly what errors this approximation may yield.
Stochastic cell rescaling
^^^^^^^^^^^^^^^^^^^^^^^^^
The stochastic cell rescaling algorithm is a variant of the Berendsen
algorithm that allows correct fluctuations to be sampled. Similarly
to the Berendsen algorithm, it rescales the coordinates and box vectors
every step, or every :math:`n_\mathrm{PC}` steps
with the effect of a firstorder kinetic relaxation of the
pressure towards a given reference pressure :math:`P_0`.
At variance with the Berendsen algorithm, the rescaling matrix is calculated
including a stochastic term that makes volume fluctuations correct.
The isotropic version can be easily written in term of the strain
:math:`\epsilon=\log(V/V_0)` that is evolved according to the following equation
of motion
.. math:: \mbox{d}\epsilon=\frac{\beta}{\tau_p}(P_0P)\mbox{d}t + \sqrt{\frac{2k_BT\beta}{V\tau_p}}\mbox{d}W
:label: eqnstochasticcellrescaling
Here, :math:`\beta` is the isothermal compressibility of the system.
It suffices to take a rough estimate because the value of :math:`\beta` only influences
the noncritical time constant of the pressure relaxation without affecting
the volume distribution itself. For water at 1 atm and 300 K
:math:`\beta = 4.6 \times 10^{10}`
Pa\ :math:`^{1} = 4.6 \times 10^{5}` bar\ :math:`^{1}`, which is
:math:`7.6 \times 10^{4}` MD units (see chapter :ref:`defunits`). Most
other liquids have similar values.
Another difference with respect to the Berendsen algorithm is that
velocities are scaled with a factor that is the reciprocal of the
scaling factor for positions.
A semiisotropic implementation is also provided. By defining the variables
:math:`\epsilon_{xy}=\log(A/A_0)` and :math:`\epsilon_z=\log(L/L_0)`,
where :math:`A` and :math:`L` are the area of the simulation box in the
:math:`xy` plane and its height, respectively, the following equations
can be obtained:
.. math:: \mbox{d}\epsilon_{xy}=\frac{2\beta}{3\tau_p}(P_0\frac{\gamma}{L}\frac{P_{xx}+P_{yy}}{2})\mbox{d}t+\sqrt{\frac{4k_BT\beta}{3V\tau_p}}\mbox{d}W_{xy}
:label: eqnstochasticcellrescalingxy
.. math:: \mbox{d}\epsilon_z=\frac{\beta}{3\tau_p}(P_0P_{zz})\mbox{d}t + \sqrt{\frac{2k_BT\beta}{3V\tau_p}}\mbox{d}W_z
:label: eqnstochasticcellrescalingz
Here :math:`\gamma` is the external surface tension and :math:`P_{xx}`,
:math:`P_{yy}`, and :math:`P_{zz}` the components of the internal pressure.
More detailed explanations can be found in the original reference \ :ref:`184 `.
ParrinelloRahman pressure coupling
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
In cases where the fluctuations in pressure or volume are important *per
se* (*e.g.* to calculate thermodynamic properties), especially for small
systems, it may be a problem that the exact ensemble is not well defined
for the weakcoupling scheme, and that it does not simulate the true NPT
ensemble.
Gromacs also supports constantpressure simulations using the
ParrinelloRahman approach \ :ref:`38 `,
:ref:`39 `, which is similar to the NoséHoover temperature
coupling, and in theory gives the true NPT ensemble. With the
ParrinelloRahman barostat, the box vectors as represented by the matrix
obey the matrix equation of motion [2]_
.. math:: \frac{{\mbox{d}}\mathbf{b}^2}{{\mbox{d}}t^2}= V \mathbf{W}^{1} \mathbf{b}'^{1} \left( \mathbf{P}  \mathbf{P}_{ref}\right).
:label: eqnPRpressure
The volume of the box is denoted :math:`V`, and
:math:`\mathbf{W}` is a matrix parameter that determines
the strength of the coupling. The matrices and :math:`_{ref}` are the
current and reference pressures, respectively.
The equations of motion for the particles are also changed, just as for
the NoséHoover coupling. In most cases you would combine the
ParrinelloRahman barostat with the NoséHoover thermostat, but to keep
it simple we only show the ParrinelloRahman modification here. The
modified Hamiltonian, which will be conserved, is:
.. math:: E_\mathrm{pot} + E_\mathrm{kin} + \sum_i P_{ii} V +
\sum_{i,j} \frac{1}{2} W_{ij} \left( \frac{{\mbox{d}}b_{ij}}{{\mbox{d}}t} \right)^2
:label: eqnPRpressureconserved
The equations of motion for the atoms, obtained from the Hamiltonian
are:
.. math:: \begin{aligned}
\frac {{\mbox{d}}^2\mathbf{r}_i}{{\mbox{d}}t^2} & = & \frac{\mathbf{F}_i}{m_i} 
\mathbf{M} \frac{{\mbox{d}}\mathbf{r}_i}{{\mbox{d}}t} , \\ \mathbf{M} & = & \mathbf{b}^{1} \left[
\mathbf{b} \frac{{\mbox{d}}\mathbf{b}'}{{\mbox{d}}t} + \frac{{\mbox{d}}\mathbf{b}}{{\mbox{d}}t} \mathbf{b}'
\right] \mathbf{b}'^{1}.
\end{aligned}
:label: eqnPRpressuremotion
This extra term has the appearance of a friction, but it should be
noted that it is ficticious, and rather an effect of the
ParrinelloRahman equations of motion being defined with all particle
coordinates represented relative to the box vectors, while Gromacs uses
normal Cartesian coordinates for positions, velocities and forces. It is
worth noting that the kinetic energy too should formally be calculated
based on velocities relative to the box vectors. This can have an effect
e.g. for external constant stress, but for now we only support coupling
to constant external pressures, and for any normal simulation the
velocities of box vectors should be extremely small compared to particle
velocities. Gang Liu has done some work on deriving this for Cartesian
coordinates\ :ref:`40 ` that we will try to implement at some
point in the future together with support for external stress.
The (inverse) mass parameter matrix
:math:`\mathbf{W}^{1}` determines the strength of the
coupling, and how the box can be deformed. The box restriction
(:eq:`%s `) will be fulfilled automatically if the corresponding
elements of :math:`\mathbf{W}^{1}` are zero. Since the
coupling strength also depends on the size of your box, we prefer to
calculate it automatically in Gromacs. You only have to provide the
approximate isothermal compressibilities :math:`\beta` and the pressure
time constant :math:`\tau_p` in the input file (:math:`L` is the largest
box matrix element):
.. math:: \left(
\mathbf{W}^{1} \right)_{ij} = \frac{4 \pi^2 \beta_{ij}}{3 \tau_p^2 L}.
:label: eqnPRpressuretimeconst
Just as for the NoséHoover thermostat, you should realize that the
ParrinelloRahman time constant is *not* equivalent to the relaxation
time used in the Berendsen pressure coupling algorithm. In most cases
you will need to use a 4–5 times larger time constant with
ParrinelloRahman coupling. If your pressure is very far from
equilibrium, the ParrinelloRahman coupling may result in very large box
oscillations that could even crash your run. In that case you would have
to increase the time constant, or (better) use the weakcoupling scheme
to reach the target pressure, and then switch to ParrinelloRahman
coupling once the system is in equilibrium. Additionally, using the
leapfrog algorithm, the pressure at time :math:`t` is not available
until after the time step has completed, and so the pressure from the
previous step must be used, which makes the algorithm not directly
reversible, and may not be appropriate for high precision thermodynamic
calculations.
Surfacetension coupling
^^^^^^^^^^^^^^^^^^^^^^^^
When a periodic system consists of more than one phase, separated by
surfaces which are parallel to the :math:`xy`plane, the surface tension
and the :math:`z`component of the pressure can be coupled to a pressure
bath. Presently, this only works with the Berendsen pressure coupling
algorithm in Gromacs. The average surface tension :math:`\gamma(t)` can
be calculated from the difference between the normal and the lateral
pressure
.. math:: \begin{aligned}
\gamma(t) & = &
\frac{1}{n} \int_0^{L_z}
\left\{ P_{zz}(z,t)  \frac{P_{xx}(z,t) + P_{yy}(z,t)}{2} \right\} \mbox{d}z \\
& = &
\frac{L_z}{n} \left\{ P_{zz}(t)  \frac{P_{xx}(t) + P_{yy}(t)}{2} \right\},\end{aligned}
:label: eqnsurftenscoupl
where :math:`L_z` is the height of the box and :math:`n` is the number
of surfaces. The pressure in the zdirection is corrected by scaling the
height of the box with :math:`\mu_{zz}`
.. math:: \Delta P_{zz} = \frac{\Delta t}{\tau_p} \{ P_{0zz}  P_{zz}(t) \}
:label: eqnzpressure
.. math:: \mu_{zz} = 1 + \beta_{zz} \Delta P_{zz}
:label: eqnzpressure2
This is similar to normal pressure coupling, except that the factor of
:math:`1/3` is missing. The pressure correction in the
:math:`z`direction is then used to get the correct convergence for the
surface tension to the reference value :math:`\gamma_0`. The correction
factor for the box length in the :math:`x`/:math:`y`direction is
.. math:: \mu_{x/y} = 1 + \frac{\Delta t}{2\,\tau_p} \beta_{x/y}
\left( \frac{n \gamma_0}{\mu_{zz} L_z}
 \left\{ P_{zz}(t)+\Delta P_{zz}  \frac{P_{xx}(t) + P_{yy}(t)}{2} \right\}
\right)
:label: eqnboxlengthcorr
The value of :math:`\beta_{zz}` is more critical than with normal
pressure coupling. Normally an incorrect compressibility will just scale
:math:`\tau_p`, but with surface tension coupling it affects the
convergence of the surface tension. When :math:`\beta_{zz}` is set to
zero (constant box height), :math:`\Delta P_{zz}` is also set to zero,
which is necessary for obtaining the correct surface tension.
MTTK pressure control algorithms
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
As mentioned in the previous section, one weakness of leapfrog
integration is in constant pressure simulations, since the pressure
requires a calculation of both the virial and the kinetic energy at the
full time step; for leapfrog, this information is not available until
*after* the full timestep. Velocity Verlet does allow the calculation,
at the cost of an extra round of global communication, and can compute,
mod any integration errors, the true NPT ensemble.
The full equations, combining both pressure coupling and temperature
coupling, are taken from Martyna *et al.* :ref:`35 ` and
Tuckerman \ :ref:`41 ` and are referred to here as MTTK
equations (MartynaTuckermanTobiasKlein). We introduce for convenience
:math:`\epsilon = (1/3)\ln (V/V_0)`, where :math:`V_0` is a reference
volume. The momentum of :math:`\epsilon` is
:math:`{v_{\epsilon}}= p_{\epsilon}/W =
\dot{\epsilon} = \dot{V}/3V`, and define :math:`\alpha = 1 + 3/N_{dof}`
(see Ref \ :ref:`41 `)
The isobaric equations are
.. math:: \begin{aligned}
\dot{{{\mathbf{r}}}}_i &=& \frac{{{\mathbf{p}}}_i}{m_i} + \frac{{p_{\epsilon}}}{W} {{\mathbf{r}}}_i \nonumber \\
\frac{\dot{{{\mathbf{p}}}}_i}{m_i} &=& \frac{1}{m_i}{{\mathbf{F}}}_i  \alpha\frac{{p_{\epsilon}}}{W} \frac{{{\mathbf{p}}}_i}{m_i} \nonumber \\
\dot{\epsilon} &=& \frac{{p_{\epsilon}}}{W} \nonumber \\
\frac{\dot{{p_{\epsilon}}}}{W} &=& \frac{3V}{W}(P_{\mathrm{int}}  P) + (\alpha1)\left(\sum_{n=1}^N\frac{{{\mathbf{p}}}_i^2}{m_i}\right),\\\end{aligned}
:label: eqnMTTKisobaric
where
.. math:: \begin{aligned}
P_{\mathrm{int}} &=& P_{\mathrm{kin}} P_{\mathrm{vir}} = \frac{1}{3V}\left[\sum_{i=1}^N \left(\frac{{{\mathbf{p}}}_i^2}{2m_i}  {{\mathbf{r}}}_i \cdot {{\mathbf{F}}}_i\
\right)\right].\end{aligned}
:label: eqnMTTKisobaric2
The terms including :math:`\alpha` are required to make phase space
incompressible \ :ref:`41 `. The :math:`\epsilon`
acceleration term can be rewritten as
.. math:: \begin{aligned}
\frac{\dot{{p_{\epsilon}}}}{W} &=& \frac{3V}{W}\left(\alpha P_{\mathrm{kin}}  P_{\mathrm{vir}}  P\right)\end{aligned}
:label: eqnMTTKaccel
In terms of velocities, these equations become
.. math:: \begin{aligned}
\dot{{{\mathbf{r}}}}_i &=& {{\mathbf{v}}}_i + {v_{\epsilon}}{{\mathbf{r}}}_i \nonumber \\
\dot{{{\mathbf{v}}}}_i &=& \frac{1}{m_i}{{\mathbf{F}}}_i  \alpha{v_{\epsilon}}{{\mathbf{v}}}_i \nonumber \\
\dot{\epsilon} &=& {v_{\epsilon}}\nonumber \\
\dot{{v_{\epsilon}}} &=& \frac{3V}{W}(P_{\mathrm{int}}  P) + (\alpha1)\left( \sum_{n=1}^N \frac{1}{2} m_i {{\mathbf{v}}}_i^2\right)\nonumber \\
P_{\mathrm{int}} &=& P_{\mathrm{kin}}  P_{\mathrm{vir}} = \frac{1}{3V}\left[\sum_{i=1}^N \left(\frac{1}{2} m_i{{\mathbf{v}}}_i^2  {{\mathbf{r}}}_i \cdot {{\mathbf{F}}}_i\right)\right]\end{aligned}
:label: eqnMTTKvel
For these equations, the conserved quantity is
.. math:: \begin{aligned}
H = \sum_{i=1}^{N} \frac{{{\mathbf{p}}}_i^2}{2m_i} + U\left({{\mathbf{r}}}_1,{{\mathbf{r}}}_2,\ldots,{{\mathbf{r}}}_N\right) + \frac{p_\epsilon}{2W} + PV\end{aligned}
:label: eqnMTTKconserved
The next step is to add temperature control. Adding NoséHoover chains,
including to the barostat degree of freedom, where we use :math:`\eta`
for the barostat NoséHoover variables, and :math:`Q^{\prime}` for the
coupling constants of the thermostats of the barostats, we get
.. math:: \begin{aligned}
\dot{{{\mathbf{r}}}}_i &=& \frac{{{\mathbf{p}}}_i}{m_i} + \frac{{p_{\epsilon}}}{W} {{\mathbf{r}}}_i \nonumber \\
\frac{\dot{{{\mathbf{p}}}}_i}{m_i} &=& \frac{1}{m_i}{{\mathbf{F}}}_i  \alpha\frac{{p_{\epsilon}}}{W} \frac{{{\mathbf{p}}}_i}{m_i}  \frac{p_{\xi_1}}{Q_1}\frac{{{\mathbf{p}}}_i}{m_i}\nonumber \\
\dot{\epsilon} &=& \frac{{p_{\epsilon}}}{W} \nonumber \\
\frac{\dot{{p_{\epsilon}}}}{W} &=& \frac{3V}{W}(\alpha P_{\mathrm{kin}}  P_{\mathrm{vir}}  P) \frac{p_{\eta_1}}{Q^{\prime}_1}{p_{\epsilon}}\nonumber \\
\dot{\xi}_k &=& \frac{p_{\xi_k}}{Q_k} \nonumber \\
\dot{\eta}_k &=& \frac{p_{\eta_k}}{Q^{\prime}_k} \nonumber \\
\dot{p}_{\xi_k} &=& G_k  \frac{p_{\xi_{k+1}}}{Q_{k+1}} \;\;\;\; k=1,\ldots, M1 \nonumber \\
\dot{p}_{\eta_k} &=& G^\prime_k  \frac{p_{\eta_{k+1}}}{Q^\prime_{k+1}} \;\;\;\; k=1,\ldots, M1 \nonumber \\
\dot{p}_{\xi_M} &=& G_M \nonumber \\
\dot{p}_{\eta_M} &=& G^\prime_M, \nonumber \\\end{aligned}
:label: eqnMTTKthermandbar
where
.. math:: \begin{aligned}
P_{\mathrm{int}} &=& P_{\mathrm{kin}}  P_{\mathrm{vir}} = \frac{1}{3V}\left[\sum_{i=1}^N \left(\frac{{{\mathbf{p}}}_i^2}{2m_i}  {{\mathbf{r}}}_i \cdot {{\mathbf{F}}}_i\right)\right] \nonumber \\
G_1 &=& \sum_{i=1}^N \frac{{{\mathbf{p}}}^2_i}{m_i}  N_f kT \nonumber \\
G_k &=& \frac{p^2_{\xi_{k1}}}{2Q_{k1}}  kT \;\; k = 2,\ldots,M \nonumber \\
G^\prime_1 &=& \frac{{p_{\epsilon}}^2}{2W}  kT \nonumber \\
G^\prime_k &=& \frac{p^2_{\eta_{k1}}}{2Q^\prime_{k1}}  kT \;\; k = 2,\ldots,M\end{aligned}
:label: eqnMTTKthermandbar2
The conserved quantity is now
.. math:: \begin{aligned}
H = \sum_{i=1}^{N} \frac{{{\mathbf{p}}}_i}{2m_i} + U\left({{\mathbf{r}}}_1,{{\mathbf{r}}}_2,\ldots,{{\mathbf{r}}}_N\right) + \frac{p^2_\epsilon}{2W} + PV + \nonumber \\
\sum_{k=1}^M\frac{p^2_{\xi_k}}{2Q_k} +\sum_{k=1}^M\frac{p^2_{\eta_k}}{2Q^{\prime}_k} + N_{f}kT\xi_1 + kT\sum_{i=2}^M \xi_k + kT\sum_{k=1}^M \eta_k\end{aligned}
:label: eqnMTTKthermandbarconserved
Returning to the Trotter decomposition formalism, for pressure control
and temperature control \ :ref:`35 ` we get:
.. math:: \begin{aligned}
iL = iL_1 + iL_2 + iL_{\epsilon,1} + iL_{\epsilon,2} + iL_{\mathrm{NHCbaro}} + iL_{\mathrm{NHC}}\end{aligned}
:label: eqnMTTKthermandbarTrotter
where “NHCbaro” corresponds to the NosèHoover chain of the barostat,
and NHC corresponds to the NHC of the particles,
.. math:: \begin{aligned}
iL_1 &=& \sum_{i=1}^N \left[\frac{{{\mathbf{p}}}_i}{m_i} + \frac{{p_{\epsilon}}}{W}{{\mathbf{r}}}_i\right]\cdot \frac{\partial}{\partial {{\mathbf{r}}}_i} \\
iL_2 &=& \sum_{i=1}^N {{\mathbf{F}}}_i  \alpha \frac{{p_{\epsilon}}}{W}{{\mathbf{p}}}_i \cdot \frac{\partial}{\partial {{\mathbf{p}}}_i} \\
iL_{\epsilon,1} &=& \frac{p_\epsilon}{W} \frac{\partial}{\partial \epsilon}\\
iL_{\epsilon,2} &=& G_{\epsilon} \frac{\partial}{\partial p_\epsilon}\end{aligned}
:label: eqnMTTKthermandbarTrotter2
and where
.. math:: \begin{aligned}
G_{\epsilon} = 3V\left(\alpha P_{\mathrm{kin}}  P_{\mathrm{vir}}  P\right)\end{aligned}
:label: eqnMTTKthermandbarTrotter3
Using the Trotter decomposition, we get
.. math:: \begin{aligned}
\exp(iL{\Delta t}) &=& \exp\left(iL_{\mathrm{NHCbaro}}{\Delta t}/2\right)\exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \nonumber \nonumber \\
&&\exp\left(iL_{\epsilon,2}{\Delta t}/2\right) \exp\left(iL_2 {\Delta t}/2\right) \nonumber \nonumber \\
&&\exp\left(iL_{\epsilon,1}{\Delta t}\right) \exp\left(iL_1 {\Delta t}\right) \nonumber \nonumber \\
&&\exp\left(iL_2 {\Delta t}/2\right) \exp\left(iL_{\epsilon,2}{\Delta t}/2\right) \nonumber \nonumber \\
&&\exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right)\exp\left(iL_{\mathrm{NHCbaro}}{\Delta t}/2\right) + \mathcal{O}({\Delta t}^3)\end{aligned}
:label: eqnMTTKthermandbarTrotterdecomp
The action of :math:`\exp\left(iL_1 {\Delta t}\right)` comes from the
solution of the differential equation
:math:`\dot{{{\mathbf{r}}}}_i = {{\mathbf{v}}}_i + {v_{\epsilon}}{{\mathbf{r}}}_i`
with
:math:`{{\mathbf{v}}}_i = {{\mathbf{p}}}_i/m_i`
and :math:`{v_{\epsilon}}` constant with initial condition
:math:`{{\mathbf{r}}}_i(0)`, evaluate at
:math:`t=\Delta t`. This yields the evolution
.. math:: {{\mathbf{r}}}_i({\Delta t}) = {{\mathbf{r}}}_i(0)e^{{v_{\epsilon}}{\Delta t}} + \Delta t {{\mathbf{v}}}_i(0) e^{{v_{\epsilon}}{\Delta t}/2} {\frac{\sinh{\left( {v_{\epsilon}}{\Delta t}/2\right)}}{{v_{\epsilon}}{\Delta t}/2}}.
:label: eqnMTTKthermandbarTrotterevol
The action of :math:`\exp\left(iL_2 {\Delta t}/2\right)` comes from the
solution of the differential equation
:math:`\dot{{{\mathbf{v}}}}_i = \frac{{{\mathbf{F}}}_i}{m_i} 
\alpha{v_{\epsilon}}{{\mathbf{v}}}_i`, yielding
.. math:: {{\mathbf{v}}}_i({\Delta t}/2) = {{\mathbf{v}}}_i(0)e^{\alpha{v_{\epsilon}}{\Delta t}/2} + \frac{\Delta t}{2m_i}{{\mathbf{F}}}_i(0) e^{\alpha{v_{\epsilon}}{\Delta t}/4}{\frac{\sinh{\left( \alpha{v_{\epsilon}}{\Delta t}/4\right)}}{\alpha{v_{\epsilon}}{\Delta t}/4}}.
:label: eqnMTTKthermandbarTrotterevol2
*mdvvavek* uses the full step kinetic energies for determining the
pressure with the pressure control, but the halfstepaveraged kinetic
energy for the temperatures, which can be written as a Trotter
decomposition as
.. math:: \begin{aligned}
\exp(iL{\Delta t}) &=& \exp\left(iL_{\mathrm{NHCbaro}}{\Delta t}/2\right)\nonumber \exp\left(iL_{\epsilon,2}{\Delta t}/2\right) \exp\left(iL_2 {\Delta t}/2\right) \nonumber \\
&&\exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \exp\left(iL_{\epsilon,1}{\Delta t}\right) \exp\left(iL_1 {\Delta t}\right) \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \nonumber \\
&&\exp\left(iL_2 {\Delta t}/2\right) \exp\left(iL_{\epsilon,2}{\Delta t}/2\right) \exp\left(iL_{\mathrm{NHCbaro}}{\Delta t}/2\right) + \mathcal{O}({\Delta t}^3)\end{aligned}
:label: eqnvvavekTrotterdecomp
With constraints, the equations become significantly more complicated,
in that each of these equations need to be solved iteratively for the
constraint forces. Before Gromacs 5.1, these iterative constraints were
solved as described in \ :ref:`42 `. From Gromacs 5.1 onward,
MTTK with constraints has been removed because of numerical stability
issues with the iterations.
Infrequent evaluation of temperature and pressure coupling
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Temperature and pressure control require global communication to compute
the kinetic energy and virial, which can become costly if performed
every step for large systems. We can rearrange the Trotter decomposition
to give alternate symplectic, reversible integrator with the coupling
steps every :math:`n` steps instead of every steps. These new
integrators will diverge if the coupling time step is too large, as the
auxiliary variable integrations will not converge. However, in most
cases, long coupling times are more appropriate, as they disturb the
dynamics less \ :ref:`35 `.
Standard velocity Verlet with NoséHoover temperature control has a
Trotter expansion
.. math:: \begin{aligned}
\exp(iL{\Delta t}) &\approx& \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right) \exp\left(iL_2 {\Delta t}/2\right) \nonumber \\
&&\exp\left(iL_1 {\Delta t}\right) \exp\left(iL_2 {\Delta t}/2\right) \exp\left(iL_{\mathrm{NHC}}{\Delta t}/2\right).\end{aligned}
:label: eqnVVNHTrotter
If the NoséHoover chain is sufficiently slow with respect to the
motions of the system, we can write an alternate integrator over
:math:`n` steps for velocity Verlet as
.. math:: \begin{aligned}
\exp(iL{\Delta t}) &\approx& (\exp\left(iL_{\mathrm{NHC}}(n{\Delta t}/2)\right)\left[\exp\left(iL_2 {\Delta t}/2\right)\right. \nonumber \\
&&\left.\exp\left(iL_1 {\Delta t}\right) \exp\left(iL_2 {\Delta t}/2\right)\right]^n \exp\left(iL_{\mathrm{NHC}}(n{\Delta t}/2)\right).\end{aligned}
:label: eqnVVNHTrotter2
For pressure control, this becomes
.. math:: \begin{aligned}
\exp(iL{\Delta t}) &\approx& \exp\left(iL_{\mathrm{NHCbaro}}(n{\Delta t}/2)\right)\exp\left(iL_{\mathrm{NHC}}(n{\Delta t}/2)\right) \nonumber \nonumber \\
&&\exp\left(iL_{\epsilon,2}(n{\Delta t}/2)\right) \left[\exp\left(iL_2 {\Delta t}/2\right)\right. \nonumber \nonumber \\
&&\exp\left(iL_{\epsilon,1}{\Delta t}\right) \exp\left(iL_1 {\Delta t}\right) \nonumber \nonumber \\
&&\left.\exp\left(iL_2 {\Delta t}/2\right)\right]^n \exp\left(iL_{\epsilon,2}(n{\Delta t}/2)\right) \nonumber \nonumber \\
&&\exp\left(iL_{\mathrm{NHC}}(n{\Delta t}/2)\right)\exp\left(iL_{\mathrm{NHCbaro}}(n{\Delta t}/2)\right),\end{aligned}
:label: eqnVVNpressure
where the box volume integration occurs every step, but the auxiliary
variable integrations happen every :math:`n` steps.
The complete update algorithm
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
.. _gmx_md_update:
**THE UPDATE ALGORITHM**

Given:
Positions :math:`\mathbf{r}` of all atoms at time
:math:`t`
Velocities :math:`\mathbf{v}` of all atoms at time
:math:`t{{\frac{1}{2}}{{\Delta t}}}`
Accelerations :math:`\mathbf{F}/m` on all atoms at time
:math:`t`.
(Forces are computed disregarding any constraints)
Total kinetic energy and virial at :math:`t{{\Delta t}}`
:math:`\Downarrow`
1. Compute the scaling factors :math:`\lambda` and :math:`\mu`
according to :eq:`eqns. %s ` and :eq:`%s `
:math:`\Downarrow`
2. Update and scale velocities:
:math:`\mathbf{v}' = \lambda (\mathbf{v} +
\mathbf{a} \Delta t)`
:math:`\Downarrow`
3. Compute new unconstrained coordinates:
:math:`\mathbf{r}' = \mathbf{r} + \mathbf{v}'
\Delta t`
:math:`\Downarrow`
4. Apply constraint algorithm to coordinates:
constrain(\ :math:`\mathbf{r}^{'} \rightarrow \mathbf{r}'';
\, \mathbf{r}`)
:math:`\Downarrow`
5. Correct velocities for constraints:
:math:`\mathbf{v} = (\mathbf{r}'' 
\mathbf{r}) / \Delta t`
:math:`\Downarrow`
6. Scale coordinates and box:
:math:`\mathbf{r} = \mu \mathbf{r}''; \mathbf{b} =
\mu \mathbf{b}`
The complete algorithm for the update of velocities and coordinates is
given using leapfrog in :ref:`the outline above `
The SHAKE algorithm of step 4 is explained below.
Gromacs has a provision to *freeze* (prevent motion of) selected
particles, which must be defined as a *freeze group*. This is
implemented using a *freeze factor* :math:`\mathbf{f}_g`,
which is a vector, and differs for each freeze group (see
sec. :ref:`groupconcept`). This vector contains only zero (freeze) or one
(don’t freeze). When we take this freeze factor and the external
acceleration :math:`\mathbf{a}_h` into account the update
algorithm for the velocities becomes
.. math:: \mathbf{v}(t+{\frac{\Delta t}{2}})~=~\mathbf{f}_g * \lambda * \left[ \mathbf{v}(t{\frac{\Delta t}{2}}) +\frac{\mathbf{F}(t)}{m}\Delta t + \mathbf{a}_h \Delta t \right],
:label: eqntotalupdate
where :math:`g` and :math:`h` are group indices which differ per atom.
Output step
~~~~~~~~~~~
The most important output of the MD run is the *trajectory file*, which
contains particle coordinates and (optionally) velocities at regular
intervals. The trajectory file contains frames that could include
positions, velocities and/or forces, as well as information about the
dimensions of the simulation volume, integration step, integration time,
etc. The interpretation of the time varies with the integrator chosen,
as described above. For Velocity Verlet integrators, velocities labeled
at time :math:`t` are for that time. For other integrators (e.g.
leapfrog, stochastic dynamics), the velocities labeled at time
:math:`t` are for time :math:`t  {{\frac{1}{2}}{{\Delta t}}}`.
Since the trajectory files are lengthy, one should not save every step!
To retain all information it suffices to write a frame every 15 steps,
since at least 30 steps are made per period of the highest frequency in
the system, and Shannon’s sampling theorem states that two samples per
period of the highest frequency in a bandlimited signal contain all
available information. But that still gives very long files! So, if the
highest frequencies are not of interest, 10 or 20 samples per ps may
suffice. Be aware of the distortion of highfrequency motions by the
*stroboscopic effect*, called *aliasing*: higher frequencies are
mirrored with respect to the sampling frequency and appear as lower
frequencies.
Gromacs can also write reducedprecision coordinates for a subset of the
simulation system to a special compressed trajectory file format. All
the other tools can read and write this format. See the User Guide for
details on how to set up your :ref:`mdp` file to have :ref:`mdrun ` use this feature.
.. [1]
Note that some derivations, an alternative notation
:math:`\xi_{\mathrm{alt}} = v_{\xi} = p_{\xi}/Q` is used.
.. [2]
The box matrix representation in corresponds to the transpose of the
box matrix representation in the paper by Nosé and Klein. Because of
this, some of our equations will look slightly different.