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:

$${\mathbf{v}}_i(t+\mathrm{\Delta }t/2)={\mathbf{v}}_i(t-\mathrm{\Delta }t/2)+{\mathbf{f}}_i(t)\mathrm{\Delta }t/m_i$$

$${\mathbf{r}}_i(t+\mathrm{\Delta }t)={\mathbf{r}}_i(t)+{\mathbf{v}}_i(t+\mathrm{\Delta }t/2)\mathrm{\Delta }t$$

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 $\mathrm{\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

$${\mathbf{v}}_i(t)=\frac{{\mathbf{v}}_i(t-\mathrm{\Delta }t/2)+{\mathbf{v}}_i(t+\mathrm{\Delta }t/2)}{2}$$

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:

$$m_i\frac{d{\mathbf{v}}_i}{dt}={\mathbf{f}}_i-\gamma m_i{\mathbf{v}}_i+{\mathbf{R}}_i$$

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 unit variance. T is the temperature of the heat bath.

The integration is done using the Langevin leap-frog method. [57] In each step, the positions and velocities are updated as follows:

$${\mathbf{v}}_i(t+\mathrm{\Delta }t/2)={\mathbf{v}}_i(t-\mathrm{\Delta }t/2)\alpha +{\mathbf{f}}_i(t)(1-\alpha )/(\gamma m_i)+\sqrt{k_{B}T(1-{\alpha }^2)/m_i}R$$

$${\mathbf{r}}_i(t+\mathrm{\Delta }t)={\mathbf{r}}_i(t)+{\mathbf{v}}_i(t+\mathrm{\Delta }t/2)\mathrm{\Delta }t$$

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 \mathrm{\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 similar to LangevinIntegerator, but it instead uses the LFMiddle discretization. [58] In each step, the positions and velocities are updated as follows:

$${\mathbf{v}}_i(t+\mathrm{\Delta }t/2)={\mathbf{v}}_i(t-\mathrm{\Delta }t/2)+{\mathbf{f}}_i(t)\mathrm{\Delta }t/m_i$$

$${\mathbf{r}}_i(t+\mathrm{\Delta }t/2)={\mathbf{r}}_i(t)+{\mathbf{v}}_i(t+\mathrm{\Delta }t/2)\mathrm{\Delta }t/2$$

$${{\mathbf{v}}^{\prime }}_i(t+\mathrm{\Delta }t/2)={\mathbf{v}}_i(t+\mathrm{\Delta }t/2)\alpha +\sqrt{kT(1-{\alpha }^2)/m}R$$

$${\mathbf{r}}_i(t+\mathrm{\Delta }t)={\mathbf{r}}_i(t+\mathrm{\Delta }t/2)+{{\mathbf{v}}^{\prime }}_i(t+\mathrm{\Delta }t/2)\mathrm{\Delta }t/2$$

This tends to produce more accurate sampling of configurational properties (such as free energies), but less accurate sampling of kinetic properties. Because configurational properties are much more important than kinetic ones in most simulations, this integrator is generally preferred over LangevinIntegrator. It often allows one to use a larger time step while still maintaining similar or better accuracy.

One disadvantage of this integrator is that it requires applying constraints twice per time step, compared to only once for LangevinIntegrator. This can make it slightly slower for systems that involve constraints. However, this usually is more than compensated by allowing you to use a larger time step.

21.4. NoseHooverIntegrator

Like LangevinMiddleIntegerator, this uses the LFMiddle discretization. [58] In each step, the positions and velocities are updated as follows:

$${\mathbf{v}}_i(t+\mathrm{\Delta }t/2)={\mathbf{v}}_i(t-\mathrm{\Delta }t/2)+{\mathbf{f}}_i(t)\mathrm{\Delta }t/m_i$$

$${\mathbf{r}}_i(t+\mathrm{\Delta }t/2)={\mathbf{r}}_i(t)+{\mathbf{v}}_i(t+\mathrm{\Delta }t/2)\mathrm{\Delta }t/2$$

$${{\mathbf{v}}^{\prime }}_i(t+\mathrm{\Delta }t/2)=\mathrm{scale}\times {\mathbf{v}}_i(t+\mathrm{\Delta }t/2)$$

$${\mathbf{r}}_i(t+\mathrm{\Delta }t)={\mathbf{r}}_i(t+\mathrm{\Delta }t/2)+{{\mathbf{v}}^{\prime }}_i(t+\mathrm{\Delta }t/2)\mathrm{\Delta }t/2$$

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. [59] 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}(\mathrm{\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:

$$\frac{d{\mathbf{r}}_i}{dt}=\frac{1}{\gamma m_i}{\mathbf{f}}_i+{\mathbf{R}}_i$$

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 $2k_{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:

$$error={\left(\mathrm{\Delta }t\right)}^2\sum _i\frac{|{\mathbf{f}}_i|}{m_i}$$

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 $\mathrm{\Delta }t$ that makes the error exactly equal the specified error tolerance:

$$\mathrm{\Delta }t=\sqrt{\frac{\delta }{\sum _i\frac{|{\mathbf{f}}_i|}{m_i}}}$$

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 $\mathrm{\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 $\mathrm{\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.13. 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.