21. Integrators¶
21.1. VerletIntegrator¶
VerletIntegrator implements the leap-frog Verlet integration method. The positions and velocities stored in the context are offset from each other by half a time step. In each step, they are updated as follows:
where \(\mathbf{v}_i\) is the velocity of particle i, \(\mathbf{r}_i\) is its position, \(\mathbf{f}_i\) is the force acting on it, \(m_i\) is its mass, and \(\Delta t\) is the time step.
Because the positions are always half a time step later than the velocities, care must be used when calculating the energy of the system. In particular, the potential energy and kinetic energy in a State correspond to different times, and you cannot simply add them to get the total energy of the system. Instead, it is better to retrieve States after two successive time steps, calculate the on-step velocities as
then use those velocities to calculate the kinetic energy at time t.
21.2. LangevinIntegator¶
LangevinIntegator simulates a system in contact with a heat bath by integrating the Langevin equation of motion:
where \(\mathbf{v}_i\) is the velocity of particle i, \(\mathbf{f}_i\) is the force acting on it, \(m_i\) is its mass, \(\gamma\) is the friction coefficient, and \(\mathbf{R}_i\) is an uncorrelated random force whose components are chosen from a normal distribution with mean zero and variance \(2 \gamma k_B T\). T is the temperature of the heat bath.
The integration is done using the LFMiddle discretization. [59] In each step, the positions and velocities are updated as follows:
where \(k\) is Boltzmann’s constant, \(T\) is the temperature, \(\gamma\) is the friction coefficient, \(R\) is a normally distributed random number, and \(\alpha=\exp(-\gamma\Delta t)\).
The same comments about the offset between positions and velocities apply to this integrator as to VerletIntegrator.
21.3. LangevinMiddleIntegrator¶
This integrator is identical to LangevinIntegerator. Separate classes exist only for historical reasons.
21.4. NoseHooverIntegrator¶
Like LangevinMiddleIntegerator, this uses the LFMiddle discretization. [59] In each step, the positions and velocities are updated as follows:
The universal scale factor used in the third step is determined by propagating auxilliary degrees of freedom alongside the regular particles. The original Nosé-Hoover formulation used a single harmonic oscillator for the heat bath, but this is problematic in small or stiff systems, which are non-ergodic, so the chain formulation extends this by replacing the single oscillator thermostat with a chain of connected oscillators. [60] For large systems a single oscillator (i.e. a chain length of one) will suffice, but longer chains are necessary to properly thermostat non-ergodic systems. The OpenMM default is to use a chain length of three to cover the latter case, but this can be safely reduced to increase efficiency in large systems.
The heat bath propagation is performed using a multi-timestep algorithm. Each propagation step is discretized into substeps using a factorization from Yoshida and Suzuki; the default discretization uses a \(\mathcal{O}(\Delta t^6)\) approach that uses 7 points, but 1, 3 or 5 points may also be used to increase performace, at the expense of accuracy. Each step is further subdivided into multi-timesteps with a default of 3 multi time steps per propagation; as with the number of Yoshida-Suziki points this value may be increase to increase accuracy but with additional computational expense.
21.5. BrownianIntegrator¶
BrownianIntegrator simulates a system in contact with a heat bath by integrating the Brownian equation of motion:
where \(\mathbf{r}_i\) is the position of particle i, \(\mathbf{f}_i\) is the force acting on it, \(\gamma\) is the friction coefficient, and \(\mathbf{R}_i\) is an uncorrelated random force whose components are chosen from a normal distribution with mean zero and variance \(2 k_B T/m_i \gamma\), where T is the temperature of the heat bath.
The Brownian equation of motion is derived from the Langevin equation of motion in the limit of large \(\gamma\). In that case, the velocity of a particle is determined entirely by the instantaneous force acting on it, and kinetic energy ceases to have much meaning, since it disappears as soon as the applied force is removed.
21.6. VariableVerletIntegrator¶
This is very similar to VerletIntegrator, but instead of using the same step size for every time step, it continuously adjusts the step size to keep the integration error below a user-specified tolerance. It compares the positions generated by Verlet integration with those that would be generated by an explicit Euler integrator, and takes the difference between them as an estimate of the integration error:
where \(\mathbf{f}_i\) is the force acting on particle i and \(m_i\) is its mass. (In practice, the error made by the Euler integrator is usually larger than that made by the Verlet integrator, so this tends to overestimate the true error. Even so, it can provide a useful mechanism for step size control.)
It then selects the value of \(\Delta t\) that makes the error exactly equal the specified error tolerance:
where \(\delta\) is the error tolerance. This is the largest step that may be taken consistent with the user-specified accuracy requirement.
(Note that the integrator may sometimes choose to use a smaller value for \(\Delta t\) than given above. For example, it might restrict how much the step size can grow from one step to the next, or keep the step size constant rather than increasing it by a very small amount. This behavior is not specified and may vary between Platforms. It is required, however, that \(\Delta t\) never be larger than the value given above.)
A variable time step integrator is generally superior to a fixed time step one in both stability and efficiency. It can take larger steps on average, but will automatically reduce the step size to preserve accuracy and avoid instability when unusually large forces occur. Conversely, when each uses the same step size on average, the variable time step one will usually be more accurate since the time steps are concentrated in the most difficult areas of the trajectory.
Unlike a fixed step size Verlet integrator, variable step size Verlet is not symplectic. This means that for a given average step size, it will not conserve energy as precisely over long time periods, even though each local region of the trajectory is more accurate. For this reason, it is most appropriate when precise energy conservation is not important, such as when simulating a system at constant temperature. For constant energy simulations that must maintain the energy accurately over long time periods, the fixed step size Verlet may be more appropriate.
21.7. VariableLangevinIntegrator¶
This is similar to LangevinIntegrator, but it continuously adjusts the step size using the same method as VariableVerletIntegrator. It is usually preferred over the fixed step size Langevin integrator for the reasons given above. Furthermore, because Langevin dynamics involves a random force, it can never be symplectic and therefore the fixed step size Verlet integrator’s advantages do not apply to the Langevin integrator.
21.8. CustomIntegrator¶
CustomIntegrator is a very flexible class that can be used to implement a wide range of integration methods. This includes both deterministic and stochastic integrators; Metropolized integrators; multiple time step integrators; and algorithms that must integrate additional quantities along with the particle positions and momenta.
The algorithm is specified as a series of computations that are executed in order to perform a single time step. Each computation computes the value (or values) of a variable. There are two types of variables: global variables have a single value, while per-DOF variables have a separate value for every degree of freedom (that is, every x, y, or z component of a particle). CustomIntegrator defines lots of variables you can compute and/or use in computing other variables. Some examples include the step size (global), the particle positions (per-DOF), and the force acting on each particle (per-DOF). In addition, you can define as many variables as you want for your own use.
The actual computations are defined by mathematical expressions as described in section 20.14. Several types of computations are supported:
Global: the expression is evaluated once, and the result is stored into a global variable.
Per-DOF: the expression is evaluated once for every degree of freedom, and the results are stored into a per-DOF variable.
Sum: the expression is evaluated once for every degree of freedom. The results for all degrees of freedom are added together, and the sum is stored into a global variable.
There also are other, more specialized types of computations that do not involve mathematical expressions. For example, there are computations that apply distance constraints, modifying the particle positions or velocities accordingly.
CustomIntegrator is a very powerful tool, and this description only gives a vague idea of the scope of its capabilities. For full details and examples, consult the API documentation.
21.9. DPDIntegrator¶
This implements dissipative particle dynamics (DPD). [61] It is similar to a Langevin integrator, but instead of applying friction and noise to the Cartesian coordinates of particles, it applies them to inter-particle displacements. This allows it to conserve momentum and efficiently model the hydrodynamics of coarse grained models.
The system evolves according to the equation of motion
where \(\mathbf{v}_i\) is the velocity of particle i, \(\mathbf{f}_i\) is the force acting on it, \(m_i\) is its mass, \(r_{ij}\) is the distance between particles i and j, \(\mathbf{e}_{ij}\) is a unit vector pointing from particle i to particle j, \(\gamma\) is the friction coefficient, and \(\mathbf{R}_{ij} = \mathbf{R}_{ji}\) is an uncorrelated random force whose components are chosen from a normal distribution with mean zero and unit variance. T is the temperature of the heat bath. The friction and noise are limited to pairs closer than a cutoff distance \(r_c\) using the function \(\omega(r) = 1-r/r_c\). The friction coefficient \(\gamma\) and cutoff distance \(r_c\) may be constants, or alternatively they can vary depending on the types of the particles.
The integration is done using the same LFMiddle discretization used for LangevinIntegrator. [59]
21.10. QTBIntegrator¶
This integrator implements the Adaptive Quantum Thermal Bath (adQTB) algorithm. [62] This is a fast method for approximating nuclear quantum effects by applying a Langevin thermostat whose random force varies with frequency to match the expected energy of a quantum harmonic oscillator.
It integrates the Langevin equation
Unlike LangevinIntegrator, for which the random noise force \(\mathbf{R}_i\) is uncorrelated white noise, in this case its cross correlation function is given by
where
In the limit of high temperature, \(\theta(\omega, T) \approx k_B T\), reproducing the classical behavior. In the limit of low temperature, \(\theta(\omega, T) \approx \hbar \omega / 2\), which is the ground state energy of a quantum harmonic oscillator with frequency \(\omega\).
A problem that can arise with this method is zero-point energy leakage. The thermostat drives the system toward the desired quantum energy distribution, but the classical dynamics of the system causes it to continuously relax toward the classical distribution. The result is too much energy in the low frequency modes and too little energy in the high frequency modes. A solution is to allow the two friction coefficients to differ. The coefficient \(\gamma_f\) appearing in the friction term of the Langevin equation is fixed as a constant. The coefficient \(\gamma_r\) that determines the magnitude of the random force becomes frequency dependent, \(\gamma_r(\omega)\). Its spectrum is dynamically adjusted to produce the correct distribution of energy between modes.
In practice, the simulation is divided into short segments, typically on the order of 1000 time steps. At the end of each segment, we calculate the velocity autocorrelation function \(C_{vv}(\omega)\) of each degree of freedom, as well as the correlation \(C_{vR}(\omega)\) between the velocity and random force \(\mathbf{R}\). From these we can calculate the amount by which the fluctuation-dissipation theorem is violated at each frequency:
We then select new friction coefficients according to
where \(k\) is the index of the segment and the adaptation rate \(A\) may vary between degrees of freedom. Random noise for the next segment is generated based on the new spectrum, and the simulation continues. The value of \(A\) generally needs to be determined by trial and error to find a value that produces fast adaptation and good convergence.
This process is much more robust if data is pooled over all degrees of freedom that are expected to be equivalent. One specifies a type index for each particle. \(\Delta_{FDT}(\omega)\) is computed as the average over all three coordinates of all particles with the same type. The more particles that are averaged over, the larger \(A\) can be, leading to faster adaptation.
When running simulations with adQTB, it is critical to equilibrate the system long enough for \(\gamma_r(\omega)\) to converge. Until it does, the simulation does not explore the correct distribution of states.
Simulations with adQTB tend to require a larger friction coefficient \(\gamma_f\) than is usual in fully classical simulations, typically 10 ps-1 or more. If the friction is too low, it is impossible to fully compensate for zero-point energy leakage. Even reducing the random force to zero may still leave too much energy in low frequency modes. It is important to check the converged spectrum of \(\gamma_r(\omega)\). If it is close to zero at some frequencies, that suggests \(\gamma_f\) needs to be increased.
A possible consequence of the large friction coefficient is unwanted broadening of spectral peaks, leading to a less accurate energy spectrum. This is compensated through a deconvolution procedure. [63] The target spectrum \(\theta(\omega, T)\) is replaced with a deconvolved version \(\tilde{\theta}(\omega, T)\) related to it by
An iterative procedure is used to solve for \(\tilde{\theta}(\omega, T)\), which is then used in place of \(\theta(\omega, T)\) when computing the random force.
One must be very careful when trying to compute velocity dependent thermodynamic quantities, such as the instantaneous temperature or instantaneous pressure. The standard calculations for these quantities assume the velocities follow a classical distribution. They do not produce correct results for an adQTB simulation, in which the velocities follow a quantum distribution.
