other algorithms for special purpose dynamics are described after
this.
-%\ifthenelse{\equal{\gmxlite}{1}}{}{
-%In the final \secref{par} of this chapter a few principles are
-%given on which parallelization of {\gromacs} is based. The
-%parallelization is hardly visible for the user and is therefore not
-%treated in detail.
-%} % Brace matches ifthenelse test for gmxlite
-
A few issues are of general interest. In all cases the {\em system}
must be defined, consisting of molecules. Molecules again consist of
particles with defined interaction functions. The detailed
other aspects of the algorithm, such as pair list generation, update of
velocities and positions, coupling to external temperature and
pressure, conservation of constraints.
-\ifthenelse{\equal{\gmxlite}{1}}{}{
The {\em analysis} of the data generated by an MD simulation is treated in \chref{analysis}.
-} % Brace matches ifthenelse test for gmxlite
\section{Periodic boundary conditions\index{periodic boundary conditions}}
\label{sec:pbc}
Restriction (\ref{eqn:gridrc}) ensures that only 26 images need to be
considered.
-%\ifthenelse{\equal{\gmxlite}{1}}{}{
\section{The group concept}
\label{sec:groupconcept}\index{group}
The {\gromacs} MD and analysis programs use user-defined {\em groups} of
sparingly.
All non-bonded interactions between pairs of energy-monitor groups can
-be excluded\index{exclusions}
-\ifthenelse{\equal{\gmxlite}{1}}
-{.}
-{(see details in the User Guide).}
+be excluded\index{exclusions} (see details in the User Guide).
Pairs of particles from excluded pairs of energy-monitor groups
are not put into the pair list.
This can result in a significant speedup
\end{description}
The use of groups in {\gromacs} tools is described in
\secref{usinggroups}.
-%} % Brace matches ifthenelse test for gmxlite
\section{Molecular Dynamics}
\label{sec:MD}
\subsubsection{Topology and force field}
The system topology, including a description of the force field, must
be read in.
-\ifthenelse{\equal{\gmxlite}{1}}
-{.}
-{Force fields and topologies are described in \chref{ff}
-and \ref{ch:top}, respectively.}
+Force fields and topologies are described in \chref{ff}
+and \ref{ch:top}, respectively.
All this information is static; it is never modified during the run.
\subsubsection{Coordinates and velocities}
In the group scheme, a neighbor list is generated consisting of pairs
of groups of at least one particle. These groups were originally
-\swapindex{charge}{group}s \ifthenelse{\equal{\gmxlite}{1}}{}{(see
- \secref{chargegroup})}, but with a proper treatment of long-range
+\swapindex{charge}{group}s (see
+ \secref{chargegroup}), but with a proper treatment of long-range
electrostatics, performance in unbuffered simulations is their only advantage. A pair of groups
is put into the neighbor list when their center of geometry is within
the cut-off distance. Interactions between all particle pairs (one from
are much faster than the kernels used in the group scheme for most
types of systems, particularly on newer hardware.
-\ifthenelse{\equal{\gmxlite}{1}}{}{
+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 cut-off range for all or most of
+the life time of this pair list. Such pairs can be pruned out by
+applying a cluster-pair kernel that only determines which clusters
+are in range. Because of the way the non-bonded 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 4-10 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 pair-list buffer size manually.
+
\subsubsection{Energy drift and pair-list buffering}
For a canonical (NVT) ensemble, the average energy error caused by
diffusion of $j$ particles from outside the pair-list cut-off
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^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
+interact with (bump into) 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$ for $r_b$ the Verlet buffer size, we can then
\secref{mod_nb_int}). One then has a buffer with the size equal to the
neighbor list cut-off less the longest interaction cut-off.
-} % Brace matches ifthenelse test for gmxlite
\subsubsection{Simple search\swapindexquiet{simple}{search}}
Due to \eqnsref{box_rot}{simplerc}, the vector $\rvij$
many shifts of combinations of box vectors need to be considered to find
the nearest image.
-\ifthenelse{\equal{\gmxlite}{1}}{}{
\begin{figure}
\centerline{\includegraphics[width=8cm]{plots/nstric}}
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.
-} % Brace matches ifthenelse test for gmxlite
-\ifthenelse{\equal{\gmxlite}{1}}{}{
\subsubsection{Charge groups}
\label{sec:chargegroup}\swapindexquiet{charge}{group}%
Charge groups were originally introduced to reduce cut-off artifacts
With the Verlet cut-off scheme, charge groups are ignored.
-} % Brace matches ifthenelse test for gmxlite
\subsection{Compute forces}
\label{subsec:forces}
\label{eqn:Xi}
\eeq
-\ifthenelse{\equal{\gmxlite}{1}}{}{
The {\gromacs} implementation of the virial computation is described
in \secref{virial}.
-} % Brace matches ifthenelse test for gmxlite
-
\subsection{The \swapindex{leap-frog}{integrator}}
\label{subsec:update}
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
-\ifthenelse{\equal{\gmxlite}{1}}
-{sites,}
-{sites (see \secref{rmfast}),}
+sites (see \secref{rmfast}),
which brings the shortest time step up to
PME mesh update frequency of a multiple time stepping scheme.
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 \normindex{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 Parrinello-Rahman 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 pair-list buffer, depending
+on the parameters used.
+
\subsubsection{Berendsen temperature coupling\pawsindexquiet{Berendsen}{temperature coupling}\index{weak coupling}}
The Berendsen algorithm mimics weak coupling with first-order
kinetics to an external heat bath with given temperature $T_0$.
numbers which may crash the simulation. In normal use,
$\lambda$ will always be much closer to 1.0.
+The thermostat modifies the kinetic energy at each scaling step by:
+\beq
+\Delta E_k = (\lambda - 1)^2 E_k
+\eeq
+The sum of these changes over the run needs to subtracted from the total energy
+to obtain the conserved energy quantity.
+
\subsubsection{Velocity-rescaling temperature coupling\pawsindexquiet{velocity-rescaling}{temperature coupling}}
The velocity-rescaling thermostat~\cite{Bussi2007a} is essentially a Berendsen
thermostat (see above) with an additional stochastic term that ensures
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.
-When an $NVT$ ensemble is used, the conserved energy quantity
-is written to the energy and log file.
\subsubsection{\normindex{Andersen thermostat}}
One simple way to maintain a thermostatted ensemble is to take an
should therefore generally not be used when examining kinetics or
transport properties of the system.~\cite{Basconi2013}
-% \ifthenelse{\equal{\gmxlite}{1}}{}{
\subsubsection{Nos{\'e}-Hoover temperature coupling\index{Nose-Hoover temperature coupling@Nos{\'e}-Hoover temperature coupling|see{temperature coupling, Nos{\'e}-Hoover}}{\index{temperature coupling Nose-Hoover@temperature coupling Nos{\'e}-Hoover}}\index{extended ensemble}}
The Berendsen weak-coupling algorithm is
and then using some algebra tricks to solve for some quantities are
required before they are actually calculated~\cite{Holian95}.
-% }
\subsubsection{Group temperature coupling}\index{temperature-coupling group}%
In {\gromacs} temperature coupling can be performed on groups of
\end{array}\right).
\eeq
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:
+\beq
+- \sum_{i,j} (\mu_{ij} -\delta_{ij}) P_{ij} V =
+\sum_{i,j} 2(\mu_{ij} -\delta_{ij}) \Xi_{ij}
+\eeq
+where $\delta_{ij}$ is the Kronecker delta and ${\bf \Xi}$ is the virial.
+Note that the factor 2 originates from the factor $\frac{1}{2}$
+in the virial definition (\eqnref{Xi}).
+
In {\gromacs}, the Berendsen scaling can also be done isotropically,
which means that instead of $\ve{P}$ a diagonal matrix with elements of size
ensemble, and it is not yet clear exactly what errors this approximation
may yield.
-% \ifthenelse{\equal{\gmxlite}{1}}{}{
\subsubsection{Parrinello-Rahman pressure coupling\pawsindexquiet{Parrinello-Rahman}{pressure coupling}}
In cases where the fluctuations in pressure or volume are important
for the Nos{\'e}-Hoover coupling. In most cases you would combine the
Parrinello-Rahman barostat with the Nos{\'e}-Hoover
thermostat, but to keep it simple we only show the Parrinello-Rahman
-modification here:
-
+modification here. The modified Hamiltonian, which will be conserved, is:
+\beq
+E_\mathrm{pot} + E_\mathrm{kin} + \sum_i P_{ii} V +
+\sum_{i,j} \frac{1}{2} W_{ij} \left( \frac{\de b_{ij}}{\de t} \right)^2
+\eeq
+The equations of motion for the atoms, obtained from the Hamiltonian are:
\bea \frac {\de^2\ve{r}_i}{\de t^2} & = & \frac{\ve{F}_i}{m_i} -
\ve{M} \frac{\de \ve{r}_i}{\de t} , \\ \ve{M} & = & \ve{b}^{-1} \left[
\ve{b} \frac{\de \ve{b}'}{\de t} + \frac{\de \ve{b}}{\de t} \ve{b}'
where the box volume integration occurs every step, but the auxiliary variable
integrations happen every $n$ steps.
-% } % Brace matches ifthenelse test for gmxlite
\subsection{The complete update algorithm}
the User Guide for details on how to set up your {\tt .mdp} file
to have {\tt mdrun} use this feature.
-% \ifthenelse{\equal{\gmxlite}{1}}{}{
\section{Shell molecular dynamics}
{\gromacs} can simulate \normindex{polarizability} using the
\normindex{shell model} of Dick and Overhauser~\cite{Dick58}. In such models
\begin{equation}
\ve{x}_S(n+1) ~=~ \ve{x}_S(n) + \ve{F}_S/k_b.
\end{equation}
-% } % Brace matches ifthenelse test for gmxlite
\section{Constraint algorithms\index{constraint algorithms}}
Constraints can be imposed in {\gromacs} using LINCS (default) or
Lagrange multipliers (and hence the displacements) requires the
solution of a set of coupled equations of the second degree. These are
solved iteratively by SHAKE.
-% \ifthenelse{\equal{\gmxlite}{1}}{}{
\label{subsec:SETTLE}
For the special case of rigid water molecules, that often make up more
than 80\% of the simulation system we have implemented the
This step is called RATTLE, and is covered in more detail in the
original Andersen paper~\cite{Andersen1983a}.
-% } % Brace matches ifthenelse test for gmxlite
\newcommand{\con}{\ve{g}}
\newcommand{\lenc}{\ve{d}}
-% \ifthenelse{\equal{\gmxlite}{1}}{}{
\subsection{\normindex{LINCS}}
\label{subsec:lincs}
rotates over more than a predefined angle.
This angle is set by the user in the {\tt *.mdp} file.
-% } % Brace matches ifthenelse test for gmxlite
\section{Simulated Annealing}
\newcommand{\rond}{\stackrel{\circ}{r}}
\newcommand{\ruis}{\ve{r}^G}
-% \ifthenelse{\equal{\gmxlite}{1}}{}{
\section{Stochastic Dynamics\swapindexquiet{stochastic}{dynamics}}
\label{sec:SD}
Stochastic or velocity \swapindex{Langevin}{dynamics} adds a friction
can be used. LINCS should be used for the constraints since SHAKE
will not converge for large atomic displacements.
BD is an option of the {\tt mdrun} program.
-% } % Brace matches ifthenelse test for gmxlite
\section{Energy Minimization}
\label{sec:EM}\index{energy minimization}%
$f=7.7$ kJ~mol$^{-1}$~nm$^{-1}$. A value for $\epsilon$ between 1 and
10 is acceptable.
-% \ifthenelse{\equal{\gmxlite}{1}}{}{
\subsection{Conjugate Gradient\index{conjugate gradient}}
Conjugate gradient is slower than steepest descent in the early stages
of the minimization, but becomes more efficient closer to the energy
gradient is only required for minimization prior to a normal-mode
analysis, which cannot be performed with constraints. For most other
purposes steepest descent is efficient enough.
-% } % Brace matches ifthenelse test for gmxlite
-% \ifthenelse{\equal{\gmxlite}{1}}{}{
\subsection{\normindex{L-BFGS}}
The original BFGS algorithm works by successively creating better
approximations of the inverse Hessian matrix, and moving the system to
improve the convergence, since sharp cut-offs mean the potential
function at the current coordinates is slightly different from the
previous steps used to build the inverse Hessian approximation.
-% } % Brace matches ifthenelse test for gmxlite
-% \ifthenelse{\equal{\gmxlite}{1}}{}{
\section{Normal-Mode Analysis\index{normal-mode analysis}\index{NMA}}
Normal-mode analysis~\cite{Levitt83,Go83,BBrooks83b}
can be performed using {\gromacs}, by diagonalization of the mass-weighted
normal modes can be generated with {\tt \normindex{gmx nmens}}.
An overview of normal-mode analysis and the related principal component
analysis (see \secref{covanal}) can be found in~\cite{Hayward95b}.
-% } % Brace matches ifthenelse test for gmxlite
-% \ifthenelse{\equal{\gmxlite}{1}}{}{
\section{Free energy calculations\index{free energy calculations}}
\label{sec:fecalc}
The $\lambda$-dependence for the force-field contributions is
described in detail in section \secref{feia}.
-% } % Brace matches ifthenelse test for gmxlite
-% \ifthenelse{\equal{\gmxlite}{1}}{}{
\section{Replica exchange\index{replica exchange}}
Replica exchange molecular dynamics (\normindex{REMD})
is a method that can be used to speed up
separate rank. See the manual page of {\tt mdrun} on how to use these
multinode features.
-% \ifthenelse{\equal{\gmxlite}{1}}{}{
\section{Essential Dynamics sampling\index{essential dynamics}\index{principal component analysis}\seeindexquiet{PCA}{covariance analysis}}
The results from Essential Dynamics (see \secref{covanal})
see {\tt gmx make_edi -h}.
The generated {\tt edi} input file is then passed to {\tt mdrun}.
-% } % Brace matches ifthenelse test for gmxlite
-% \ifthenelse{\equal{\gmxlite}{1}}{}{
\section{\normindex{Expanded Ensemble}}
In an expanded ensemble simulation~\cite{Lyubartsev1992}, both the coordinates and the
non-dimensional units, but they can be set to arbitrary values as
desired. Several different algorithms can be used to equilibrate
these weights, described in the mdp option listings.
-% } % Brace matches ifthenelse test for gmxlite
In {\gromacs}, this space is sampled by alternating sampling in the $k$
and $\vec{x}$ directions. Sampling in the $\vec{x}$ direction is done
with the change in domain volume and the load outside
the measured region does not depend strongly on the domain volume.
In CPU-only simulations, the load is measured between the coordinate
-and the force communication. In hybrid CPU-GPU simulations we overlap
-communication on the CPU with calculation on the GPU. Therefore we
+and the force communication. In simulations with non-bonded
+work on GPUs, we overlap communication and
+work on the CPU with calculation on the GPU. Therefore we
measure from the last communication before the force calculation to
when the CPU or GPU is finished, whichever is last.
When not using PME ranks, we subtract the time in PME from the CPU time,
Therefore a parallel constraint algorithm is required.
{\gromacs} uses the \normindex{P-LINCS} algorithm~\cite{Hess2008a},
which is the parallel version of the \normindex{LINCS} algorithm~\cite{Hess97}
-% \ifthenelse{\equal{\gmxlite}{1}}
-{.}
-{(see \ssecref{lincs}).}
+(see \ssecref{lincs}).
The P-LINCS procedure is illustrated in \figref{plincs}.
When molecules cross the cell boundaries, atoms in such molecules
up to ({\tt lincs_order + 1}) bonds away are communicated over the cell boundaries.
Electrostatics interactions are long-range, therefore special
algorithms are used to avoid summation over many atom pairs.
In {\gromacs} this is usually
-% \ifthenelse{\equal{\gmxlite}{1}}
-{.}
-{PME (\secref{pme}).}
+PME (\secref{pme}).
Since with PME all particles interact with each other, global communication
is required. This will usually be the limiting factor for
scaling with domain decomposition.
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff