%
% This file is part of the GROMACS molecular simulation package.
%
- % Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
+ % Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by
% Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
% and including many others, as listed in the AUTHORS file in the
% top-level source directory and at http://www.gromacs.org.
non-bonded forces must be calculated. The pair list contains particles
$i$, a displacement vector for particle $i$, and all particles $j$ that
are within \verb'rlist' of this particular image of particle $i$. The
-list is updated every \verb'nstlist' steps, where \verb'nstlist' is
-typically 10. There is an option to calculate the total non-bonded
-force on each particle due to all particle in a shell around the
-list cut-off, {\ie} at a distance between \verb'rlist' and
-\verb'rlistlong'. This force is calculated during the pair list update
-and retained during \verb'nstlist' steps.
+list is updated every \verb'nstlist' steps.
To make the \normindex{neighbor list}, all particles that are close
({\ie} within the neighbor list cut-off) to a given particle must be found.
\ifthenelse{\equal{\gmxlite}{1}}{}{
\subsubsection{Energy drift and pair-list buffering}
- For a canonical (NVT) ensemble, the average energy error caused by the
- finite Verlet buffer size can be determined from the atomic
+ For a canonical (NVT) ensemble, the average energy error caused by
+ diffusion of $j$ particles from outside the pair-list cut-off
+ $r_\ell$ to inside the interaction cut-off $r_c$ over the lifetime
+ of the list can be determined from the atomic
displacements and the shape of the potential at the cut-off.
%Since we are interested in the small drift regime, we will assume
%#that atoms will only move within the cut-off distance in the last step,
%$n_\mathrm{ps}-1$, of the pair list update interval $n_\mathrm{ps}$.
%Over this number of steps the displacment of an atom with mass $m$
The displacement distribution along one dimension for a freely moving
- particle with mass $m$ over time $t$ at temperature $T$ is Gaussian
- with zero mean and variance $\sigma^2 = t\,k_B T/m$. For the distance
+ particle with mass $m$ over time $t$ at temperature $T$ is
+ a Gaussian $G(x)$
+ of zero mean and variance $\sigma^2 = t^2 k_B T/m$. For the distance
between two particles, the variance changes to $\sigma^2 = \sigma_{12}^2 =
- t\,k_B T(1/m_1+1/m_2)$. Note that in practice particles usually
+ t^2 k_B T(1/m_1+1/m_2)$. Note that in practice particles usually
interact with other particles over time $t$ and therefore the real
displacement distribution is much narrower. Given a non-bonded
interaction cut-off distance of $r_c$ and a pair-list cut-off
- $r_\ell=r_c+r_b$, we can then write the average energy error after
- time $t$ for pair interactions between one particle of type 1
- surrounded by particles of type 2 with number density $\rho_2$, when
- the inter particle distance changes from $r_0$ to $r_t$, as:
-
+ $r_\ell=r_c+r_b$ for $r_b$ the Verlet buffer size, we can then
+ write the average energy error after time $t$ for all missing pair
+ interactions between a single $i$ particle of type 1 surrounded
+ by all $j$ particles that are of type 2 with number density $\rho_2$,
+ when the inter-particle distance changes from $r_0$ to $r_t$, as:
+ \beq
+ \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_t-r_0}{\sigma}\right) d r_0\, d r_t
+ \eeq
+ 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:
\begin{eqnarray}
- \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_t-r_0}{\sigma}\right) d r_0\, d r_t \\
- &\approx&
+ \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 \Big] G\!\left(\frac{r_t-r_0}{\sigma}\right) d r_0 \, d r_t\\
+ 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_t-r_0}{\sigma}\right) d r_0 \, d r_t
+ \end{eqnarray}
+ Replacing the factor $r_0^2$ by $(r_\ell + \sigma)^2$, which results in a slight overestimate, allows us to calculate the integrals analytically:
+ \begin{eqnarray}
+ \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 \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] +
+ \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)\left[ r_b\sigma(r_b^2+5\sigma^2)G\!\left(\frac{r_b}{\sigma}\right) - (r_b^4+6r_b^2\sigma^2+3\sigma^4 ) E\!\left(\frac{r_b}{\sigma}\right) \right]
+ \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{eqnarray}
- where $G$ is a Gaussian distribution with 0 mean and unit variance and
+ where $G(x)$ is a Gaussian distribution with 0 mean and unit variance and
$E(x)=\frac{1}{2}\mathrm{erfc}(x/\sqrt{2})$. We always want to achieve
small energy error, so $\sigma$ will be small compared to both $r_c$
and $r_\ell$, thus the approximations in the equations above are good,
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
- neighbor-list life time $t = ({\tt nstlist} - 1)\times{\tt dt}$. This
+ neighbor-list life time $t = ({\tt nstlist} - 1)\times{\tt dt}$. The
function can not be inverted analytically, so we use bisection to
obtain the buffer size $r_b$ for a target drift. Again we note that
in practice the error we usually be much smaller than this estimate,
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. The
- displacement in an arbitrary direction of a particle with 2 degrees of
- freedom is not Gaussian, but rather follows the complementary error
- function: \beq
+ 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:
+ \beq
\frac{\sqrt{\pi}}{2\sqrt{2}\sigma}\,\mathrm{erfc}\left(\frac{|r|}{\sqrt{2}\,\sigma}\right)
- \eeq where $\sigma^2$ is again $k_B T/m$. This distribution can no
+ \label{eqn:2D_disp}
+ \eeq
+ where $\sigma^2$ is again $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. We consider
- particles constrained, i.e. having 2 degrees of freedom or fewer, when
- they are connected by constraints to particles with a total mass of at
- least 1.5 times the mass of the particles itself. For a particle with
- a single constraint this would give a total mass along the constraint
- direction of at least 2.5, which leads to a reduction in the variance
- of the displacement along that direction by at least a factor of 6.25.
- As the Gaussian distribution decays very rapidly, this effectively
- removes one degree of freedom from the displacement. Multiple
- constraints would reduce the displacement even more, but as this gets
- very complex, we consider those as particles with 2 degrees of
- freedom.
+ 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 $\sigma$ in \eqnref{2D_disp} (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
with double precision may be required to get fine details of
thermodynamics correct.
-\subsection{Twin-range cut-offs\index{twin-range!cut-off}}
-To save computation time, slowly varying forces can be calculated
-less often than rapidly varying forces. In {\gromacs}
-such a \normindex{multiple time step} splitting is possible between
-short and long range non-bonded interactions.
-In {\gromacs} versions up to 4.0, an irreversible integration scheme
-was used which is also used by the {\gromos} simulation package:
-every $n$ steps the long range forces are determined and these are
-then also used (without modification) for the next $n-1$ integration steps
-in \eqnref{leapfrogv}. Such an irreversible scheme can result in bad energy
-conservation and, possibly, bad sampling.
-Since version 4.5, a leap-frog version of the reversible Trotter decomposition scheme~\cite{Tuckerman1992a} is used.
-In this integrator the long-range forces are determined every $n$ steps
-and are then integrated into the velocity in \eqnref{leapfrogv} using
-a time step of $\Dt_\mathrm{LR} = n \Dt$:
-\beq
-\ve{v}(t+\hDt) =
-\left\{ \begin{array}{lll} \displaystyle
- \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
- \ve{v}(t-\hDt) + \frac{1}{m}\ve{F}_\mathrm{SR}(t)\Dt &,& \mathrm{step} ~\%~ n \neq 0 \\
-\end{array} \right.
-\eeq
-
-The parameter $n$ is equal to the neighbor list update frequency. In
-4.5, the velocity Verlet version of multiple time-stepping is not yet
-fully implemented.
-
+\subsection{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 brings the shortest time step up to
PME mesh update frequency of a multiple time stepping scheme.
-As an example we show the energy conservation for integrating
-the equations of motion for SPC/E water at 300 K. To avoid cut-off
-effects, reaction-field electrostatics with $\epsilon_{RF}=\infty$ and
-shifted Lennard-Jones interactions are used, both with a buffer region.
-The long-range interactions were evaluated between 1.0 and 1.4 nm.
-In \figref{leapfrog} one can see that for electrostatics the Trotter scheme
-does an order of magnitude better up to $\Dt_{LR}$ = 16 fs.
-The electrostatics depends strongly on the orientation of the water molecules,
-which changes rapidly.
-For Lennard-Jones interactions, the energy drift is linear in $\Dt_{LR}$
-and roughly two orders of magnitude smaller than for the electrostatics.
-Lennard-Jones forces are smaller than Coulomb forces and
-they are mainly affected by translation of water molecules, not rotation.
-
-\begin{figure}
-\centerline{\includegraphics[width=12cm]{plots/drift-all}}
-\caption{Energy drift per degree of freedom in SPC/E water
-with twin-range cut-offs
-for reaction field (left) and Lennard-Jones interaction (right)
-as a function of the long-range time step length for the irreversible
-``\gromos'' scheme and a reversible Trotter scheme.}
-\label{fig:twinrangeener}
-\end{figure}
-
\subsection{Temperature coupling\index{temperature coupling}}
While direct use of molecular dynamics gives rise to the NVE (constant
number, constant volume, constant energy ensemble), most quantities
the dynamics will be completely different from MD, but the sampling is
still correct.
-In {\gromacs} there are two algorithms to integrate equation (\ref{SDeq}):
-a simple and efficient one
-and a more complex leap-frog algorithm~\cite{Gunsteren88}, which is now deprecated.
-The accuracy of both integrators is equivalent to the normal MD leap-frog and
-Velocity Verlet integrator, except with constraints where the complex
-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}.
-The simple integrator is:
+In {\gromacs} there is one simple and efficient implementation. Its
+accuracy is equivalent to the normal MD leap-frog and
+Velocity Verlet integrator. It 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}.
+It can be described as:
\bea
\label{eqn:sd_int1}
\ve{v}' &~=~& \ve{v}(t-\hDt) + \frac{1}{m}\ve{F}(t)\Dt \\
where $\ruis_i$ is Gaussian distributed noise with $\mu = 0$, $\sigma = 1$.
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.
-In the deprecated complex algorithm, four Gaussian random numbers are required
-per integration step per degree of freedom, and with constraints the
-coordinates need to be constrained twice per integration step.
-Depending on the computational cost of the force calculation,
-this can take a significant part of the simulation time.
-Exact continuation of a stochastic dynamics simulation is not possible,
-because the state of the random number generator is not stored.
-
When using SD as a thermostat, an appropriate value for $\gamma$ is e.g. 0.5 ps$^{-1}$,
since this results in a friction that is lower than the internal friction
of water, while it still provides efficient thermostatting.
% LocalWords: pymbar multinode subensemble Monte solute subst groupconcept GPU
% LocalWords: dodecahedron octahedron dodecahedra equilibration usinggroups nm
% LocalWords: topologies rlistlong CUDA GPUs rcoulomb SIMD BlueGene FPUs erfc
-% LocalWords: cutoffschemesupport unbuffered bondeds AdResS OpenMP ewald rtol
+% LocalWords: cutoffschemesupport unbuffered bondeds OpenMP ewald rtol
% LocalWords: verletdrift peptide RMS rescaling ergodicity ergodic discretized
% LocalWords: isothermal compressibility isotropically anisotropic iteratively
% LocalWords: incompressible integrations translational biomolecules NMA PCA
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff