In MD simulations the potential energy function of the atomic positions is solved by numerical integration of Newton’s equations of motion [135,136]. The equations of motion are determin-istic, in the sense that from a set of initial coordinates and an initial distribution of velocities a trajectory can be calculated that predicts the state of the system (positions, velocities and ac-celerations of the particles) at all other times. The required initial coordinates may be obtained from structures of biomolecules, which were solved experimentally by X-ray crystallography or NMR spectroscopy.
There has been a development of numerous specific MD integration algorithms most im-portantly aimed at computational efficiency and high accuracy (conservation of energy and mo-mentum). The mathematically equivalent Verlet-type algorithms (the Verlet, velocity-Verlet, and the leap-frog algorithm) are the most commonly used in classical MD today [174]. The integration in GROMACS [175–177] is performed using the leap-frog algorithm [178] and can be summarized in the following scheme.
Here ∆t denotes the time step, r(t) the particle’s coordinate vector and v(t) the respective velocities.
2.2.1 Time step, constraints and virtual sites
Integration time step. The integration time step size depends on the type of propagation algorithm, as well as the steepness of the potential. The maximal time step is in general in-versely proportional to the frequency of the fastest degree of freedom present in the system.
For an accurate and stable propagation, the value of the time step must be an order of mag-nitude smaller than the fastest motion. Exceeding this length of the time step will lead to a rapid accumulation of errors in the numerical integration and will ultimately break down the propagation, as indicated by a drift in the total energy.
The dynamics of complex biomolecular systems take place on different time scales at physiological temperature [9,179]. The fastest motions proceed on a femtosecond time scale, e.g. intra-molecular fluctuations of bonds, as well as angular and torsional oscillations.
Side chain rotations or loop and domain motions occur on a picosecond and nano- to millisecond time scale, respectively. The biologically relevant conformational changes of biomolecules happen mostly on longer time scales, which range from at least several nanoseconds to often beyond even seconds [9].
Constrained bond dynamics. The fast bond vibrations require a small integra-tion step size and thereby severely limit the accessible simulaintegra-tion time. A common way to alleviate this problem is to apply constraint algorithms, like SHAKE [180] or LINCS (LINear Constraint Solver) [181, 182] that implement geometric (holonomic) constraints, while advancing the particle coordinates. Replacing bonds with constraints is a usual procedure in simulations with the GROMACS software package, i.e. fixing the distances between bonded atoms to their equilibrium values. Thereby, the bond vibrations are removed and an increased time step of up to 2 fs can be applied. Furthermore, it has been argued that a constraint is a more faithful representation of the quantum-mechanical ground state of the bond-stretching vibration than the description in terms of a harmonic oscillator [183].
Virtual interaction sites. Most of the fastest motions in a biomolecular simula-tion necessarily involve hydrogen atoms because of their small mass. After constraining all bond-lengths, the next shortest oscillation period in a simulation is therefore the hydrogen bond-angle vibration with 13 fs [179].
Yet another way to achieve an effective increase of the integration time step is thus to define all hydrogen atoms as virtual interaction sites. The virtual site’s position is reconstructed from three predefined and nearby heavy atoms that have a fixed orientation with respect to each other, thereby removing all internal high-frequency degrees of freedom. Since only the heavy atom positions are integrated, all forces acting on the hydrogen atom will be redistributed over these particular atoms. However, a slightly different virtual site construction approach is required for hydroxyl or amine groups, since there the rotational freedom has to be persevered. By disregarding these very fast oscillations, the next shortest periods are around 20 fs, which in turn allow for a maximum time step of 5 fs, while still integrating with reasonable accuracy [175,179].
Note that for the popular explicit solvent models (SPC [184], TIP3P [185], TIP4P [186]) a completely rigid valence geometry is used as a good approximation, therefore the libration of the water molecule with a frequency of 28 fs corresponds to the fastest degree of freedom [179].
2.2.2 Ensembles and temperature
The MD algorithm presented so far describes the dynamics of an isolated system solving New-ton’s equation of motion and therefore should in principle generate a constant NVE (parti-cle number, volume, energy) or microcanonical ensemble of conformations [136]. However, biomolecular processes occur in systems which are in thermal (T) and mechanical (P for pres-sure) equilibrium with their environment. Consequently, the sampling of canonical (NVT) or isothermal-isobaric (NPT) ensembles is more suitable. In order to achieve these thermodynamic conditions in the simulation system, the solution of the equations of motion has to be modified.
Algorithms for constant temperature MD are called thermostats. Several approaches of temperature control have been presented, while each of them has certain advantages, some of them are more sophisticated and rigorous than others [187–190]. The Berendsen thermostat [187] affects the heat flow by rescaling all particle velocities to adjust the instantaneous kinetic energy of the system to the desired temperature. However, the weak coupling (first-order kinetics) to an external heat bath with given reference temperature suppresses the fluctuations of the kinetic energy and therefore does not generate proper canonical ensembles [191]. The velocity rescaling thermostat [188] is similar to the Berendsen thermostat in the sense that it imposes an exponential relaxation of temperature deviations on the system and will not produce oscillations. At the same time this algorithm generates a correct ensemble, while an additional stochastic term ensures that the correct kinetic energy distribution is obtained.
For constant pressure MD similar algorithms have been developed, such as the weak coupling scheme of Berendsen [187] but they will not be discussed here.