Time step rescaling recovers continuous-time dynamical properties for discrete-time Langevin integration of nonequilibrium systems.

When simulating molecular systems using deterministic equations of motion (e.g., Newtonian dynamics), such equations are generally numerically integrated according to a well-developed set of algorithms that share commonly agreed-upon desirable properties. However, for stochastic equations of motion (e.g., Langevin dynamics), there is still broad disagreement over which integration algorithms are most appropriate. While multiple desiderata have been proposed throughout the literature, consensus on which criteria are important is absent, and no published integration scheme satisfies all desiderata simultaneously. Additional nontrivial complications stem from simulating systems driven out of equilibrium using existing stochastic integration schemes in conjunction with recently developed nonequilibrium fluctuation theorems. Here, we examine a family of discrete time integration schemes for Langevin dynamics, assessing how each member satisfies a variety of desiderata that have been enumerated in prior efforts to construct suitable Langevin integrators. We show that the incorporation of a novel time step rescaling in the deterministic updates of position and velocity can correct a number of dynamical defects in these integrators. Finally, we identify a particular splitting (related to the velocity Verlet discretization) that has essentially universally appropriate properties for the simulation of Langevin dynamics for molecular systems in equilibrium, nonequilibrium, and path sampling contexts.

Introduction

Simulating the dynamics of molecular systems on a digital computer requires that the equations of motion be discretized. The resulting discrete-time integration algorithm that governs the updates of particle positions and velocities will necessarily have properties that differ from the continuous equations of motion on which it was based. To construct such an algorithm, one must decide which properties of the original dynamics should be preserved. Even then, a multitude of integration schemes may satisfy these properties and still recover the continuous stochastic equations of motion in the limit of an infinitesimally small time step.

For integrating the deterministic classical equations of motion prescribed by Newtonian dynamics, explicit symplectic integration schemes such as velocity Verlet are now widely regarded as being optimal for condensed matter systems for a number of reasons: they are reversible, simple to implement, preserve phase space volume, require minimal force evaluations, and are generally stable over long integration times.[1−5]

For stochastic equations of motion—in which the influence of some system components (often the solvent) are not represented explicitly, but instead by random collisions with fictitious particles—no such generally adopted integrator yet exists. In particular, the dynamics produced by existing algorithms differ (in a time step dependent manner) in important respects from the dynamics of the continuous equations of motion. Nevertheless, significant effort has been devoted to developing such stochastic integrators due to their utility in simulating many systems of interest to the chemical, biophysical, and physical sciences.

Driven nonequilibrium systems present their own set of special challenges. For example, a powerful set of nonequilibrium work fluctuation theorems[6] permit the computation of equilibrium properties of systems from their nonequilibrium statistics, but require as input the distribution of work values associated with the ensemble of trajectories. Calculations that use naive analogies of the work for continuous dynamics—failing to take into account details of the discrete integration scheme—possess systematic biases.[7]

Here, we consider various choices that could be made in constructing discrete-time integration schemes for Langevin dynamics. By examining the various possible Strang splittings of the Langevin Liouville operator, incorporating a novel time step rescaling (eq 15), and comparing their resulting properties to several desiderata that have been enumerated in the literature over recent decades, we show how it is possible to simultaneously satisfy nearly all these criteria with a single integration scheme (eq 7) that is generally applicable, simple to implement, computationally efficient, and produces thermodynamically consistent nonequilibrium statistics.

Theory

The Langevin Equation

A standard framework for the stochastic simulation of molecular systems assumes that the variables of interest evolve according to Langevin dynamics with uncorrelated Gaussian noise,[8] which represents interactions with the surrounding environment through frictional drag and stochastic collisions with fictitious bath particles:Here r and v are time-dependent position and velocity, m is mass, β = 1/kBT, kB is Boltzmann’s constant, T is the temperature of the environment, γ is a friction coefficient (with dimensions of inverse time), and W(t) is a standard Wiener process. The force f(t) is due to the (in general time-dependent) Hamiltonian (t) on the system with position r(t), as determined by the derivative of the potential energy, –∂(t)/∂r, evaluated at r(t). For multidimensional, multiparticle systems, r, v, f, and dW are vectors, and m is a diagonal matrix (see Supporting Information, “Multiple dimensions”).

While these equations of motion can be solved exactly for some simple systems, nearly all complex systems of interest require computational techniques in order to generate dynamical trajectories. On digital computers, this requires discretization of time. The selection of an appropriate discrete time integration scheme is made difficult by the fact that many discretizations may exist that recover the same continuous stochastic differential equations of motion in the limit of an infinitesimally small time step, but these schemes may possess very different properties for finite time steps.

Desiderata

Finite time step integrators for molecular systems cannot hope to exactly reproduce snapshots from the dynamical trajectories of the continuous equations of motion of a real physical system. Imprecision of experimental measurements ensures that simulated initial conditions necessarily deviate from ‘true’ ones, and Lyapunov instability of even deterministic dynamics ensures the rapid chaotic growth of such deviations, making the exact reproduction of a particular trajectory impossible even were it desirable.[3] Moreover, artifacts are inevitably introduced when one discretizes continuous equations of motion in a straightforward manner: dynamical motion increasingly diverges from that of continuous equations of motion with increasing friction and/or time step. Instead, the goal is often to reproduce certain statistical (often observable) properties of the system, especially in terms of correlation functions and (possibly time-dependent) ensemble expectations given a set of initial conditions. Thus, a desirable approximation scheme should share certain statistical and dynamical properties of the ensemble of trajectories associated with the exact equations of motion, in lieu of being able to exactly integrate trajectories.

Pastor, et al.[9] proposed that a useful discrete-time integrator should reproduce seven quantities associated with the continuous-time equations of motion: for a free particle (zero-force), the mean-squared displacement (MSD) as a function of time, the mean-squared velocity (MSV), and the velocity autocorrelation function (VAC); for a uniform external force, the terminal velocity; and for a harmonic potential (linear force), the MSD, MSV, and the virial. In Table 1, we define these desired dynamical properties and list their analytically computed values for the continuous equations of motion. In the Supporting Information, under “Determination of rescaling parameters,” we describe in detail the calculation of these quantities for our family of integrators.

Table 1

Definition of Dynamical Propertiesa

external force	quantity		expression		continuous-limit value
zero	mean-squared displacement	⟨r2(n)⟩	or	⟨r2(n + 1/2)⟩	2/(βmγ) nΔt
	mean-squared velocity	⟨v2(n)⟩	or	⟨v2(n + 1/2)⟩	1/(βm)
	velocity autocorrelation	⟨v(n)v(n + Δn)⟩	or	⟨v(n + 1/2) v(n + 1/2 + Δn)⟩	1/(βm) e–γΔnΔt
uniform, f	terminal drift	⟨r(n + 1) – r(n)⟩ / Δt	=	⟨r(n + 1/2) – r(n – 1/2)⟩ / Δt	f/(mγ)
linear, −kr	mean-squared displacement	⟨r2(n)⟩	or	⟨r2(n + 1/2)⟩	1/(βk)
	mean-squared velocity	⟨v2(n)⟩	or	⟨v2(n + 1/2)⟩	1/(βm)
	virial	m⟨v2(n)⟩ – k⟨r2(n)⟩	or	m⟨v2(n + 1/2)⟩ – k⟨r2(n + 1/2)⟩	0

Angled brackets denote an average over the ensemble of phase-space trajectories produced by a given integration scheme.

One should not expect a discrete algorithm to give meaningful results on time scales less than a single time step. Those users who prefer to treat a discrete algorithm as a black box will be most interested in integrators that satisfy given dynamical properties at integer time steps. Those willing to look further under the hood would need to know at which specific point within a given time step to measure a given dynamical quantity to recover (or most closely approximate) the continuous-limit value. Thus, we examine dynamical properties at both integer and fractional time steps.

There are several other criteria that one may want an integrator to satisfy, not the least of which being ease of implementation and analysis. Additionally, an integrator that is computationally efficient should have an accuracy that scales reasonably with time step length (here, quadratically, the same accuracy order as popular symplectic integration schemes for deterministic dynamics such as velocity Verlet), permitting relatively large time steps; minimize the number of force evaluations (one per time step) so as to minimize computational effort; and easily incorporate constraints (typically reflecting covalent chemical bonds to light elements such as hydrogen) that push the integrator stability limit to larger time step. Path sampling[10] or path reweighting[11,12] strategies often require an integrator that induces an irreducible Markov chain (i.e., it is possible to transition from any phase space point to any other in a single time step through specific choice of the random variables[13]), and that for a given trajectory has a readily evaluated path action that governs the probability of that trajectory within the dynamical ensemble.

Indeed, the task is further complicated when the integrator must produce thermodynamically consistent nonequilibrium simulations. To facilitate the use of nonequilibrium fluctuation relations[6] and estimators derived from them, the integrator must distinguish between the heat, work, and shadow work, by properly splitting dynamics into stochastic, explicit Hamiltonian-update, and deterministic substeps.[14] For the calculated works to be thermodynamically meaningful, the integrator must also have a form symmetric under time-reversal.

Pastor, et al. demonstrate that no member of their family of overdamped integrators can simultaneously satisfy more than four of their seven desired dynamical properties. Many other underdamped discretizations of the Langevin equation have been proposed[15−25] that achieve some subset of these enumerated desiderata, yet there still is no widespread agreement on an integrator that performs satisfactorily for a broad set of purposes.

Resolution: Time Step Rescaling

Drawing inspiration from several popular integrators, in this paper we derive a simple family of integrators that split the different update types to permit the definition of thermodynamically meaningful quantities for work and heat. With a novel rescaling of the time step, the resulting dynamics preserves five of Pastor et al.’s seven dynamical properties for any values of friction and time step (six if fractional time step values are also considered), and furthermore satisfies all of the other desiderata enumerated above, including its utility for nonequilibrium simulations and schemes involving path sampling or reweighting. The resulting integrator represents a stochastic generalization of velocity Verlet, is simple to implement, and could be a general all-purpose replacement for the various discrete-time Langevin integrators now in use.

In sampling contexts, this time step rescaling does slow the exploration of conformational space if the raw integration time step Δt is not concomitantly increased. We believe that such a time step rescaling permits a larger raw time step, but exactly how much larger Δt can be increased (and hence to what extent sampling efficiency can be recovered or even improved) before reaching the stability limit (or some other integration pathology) remains an open question. Addressing this potentially system-specific issue requires further theoretical and numerical investigation, which is beyond the scope of this paper.

Integrator Splitting

We derive a family of integrators by splitting the time evolution operator into stochastic and deterministic components[20] and choosing the adjustable parameters to match dynamical quantities from the continuous equations of motion. We write the one-dimensional version here, but the generalization to multiple dimensions is straightforward (see the Supporting Information, “Multiple dimensions”). The Langevin Liouville operator (sometimes termed the Liouvillian)[26] can be naturally written as a sum of four parts = + + + . The first operator represents stochastic thermalization[27] via an Ornstein–Uehlenbeck operator,the next two operators represent deterministic Newtonian evolution of velocity and position,and the last operator represents the time evolution of the system Hamiltonian according to the predetermined schedule (or protocol) Λ,where n denotes the time step index and t = nΔt the simulation clock time for time step Δt. (Note that in the case of a time-independent Hamiltonian, e is the identity operator.) Similar to several other integrators,[17,20,21,28,29] we approximate the dynamics over a time Δt by applying a series of Strang (symmetric Trotter) operator splittings:[30]where (A, B, C, D) represents a permutation of 0, , , .

There are six Strang splittings of the Liouville operators 0, , and . The Hamiltonian update operator commutes with 0 and with , so for each of these six splittings, there are only two unique placements of . For each of the six splittings, one placement of interleaves the position, Hamiltonian, and deterministic velocity updates in such a way as to require multiple force evaluations per step, making the scheme computationally inefficient. Thus there are six distinct splittings that each give rise to different finite time step dynamics and require only one force evaluation per step. Notably, because the error in each Strang splitting is (Δt3), all are identical to the true Liouville operator in the limit Δt → 0+.

One such splitting is a stochastic generalization of Velocity Verlet that, in analogy to the nomenclature of Leimkuhler and Matthews,[25] we call OVRVO (to denote the respective ordering of Ornstein–Uehlenbeck (O), deterministic velocity (V), and deterministic position (R) updates),The Hamiltonian-update step (d) is placed to minimize the number of force evaluations per time step.

For this operator splitting, a single update step that advances the simulation clock by Δt is given explicitly byHere, a = exp(−γΔt), and + and – are independent, normally distributed random variables with zero mean and unit variance (hence, when scaled by (βm)−1/2, distributed according to the equilibrium Maxwell–Boltzmann velocity distribution).

The substeps in eq 7 are the finite difference expressions of the corresponding suboperators in eq 6. The initial (a) and final (g) operators randomize the velocity and leave the position unchanged, mixing the old velocity with a Maxwell–Boltzmann random variate (with old velocity weighted according to a = exp[−γΔt]). These operators can be analytically integrated to give the first (7a) and last (7g) substeps that are stochastic, Markovian, and detailed-balanced (with respect to the true canonical measure[31]).[20,32] The operators (b) and (f) correspond to deterministic velocity updates, while (c) and (e) correspond to deterministic position updates. Together they are approximated by the finite difference expressions of eqs 7b, 7c, 7e, and 7f, which together constitute the deterministic and symplectic velocity Verlet integrator,[1,28] but slightly altered by an effective time step rescaling 0 ≤ b ≤ 1, chosen to maintain the continuous-limit zero-force diffusion coefficient and terminal drift under a uniform external force, regardless of friction coefficient or time step (derived in Supporting Information, “Determination of rescaling parameters”). In the limit Δt → 0 this rescaling factor b converges to unity and the splitting converges to the continuous-time equation of motion (eq 1). Operator (d) and its finite difference expression as eq 7d makes explicit the midpoint Hamiltonian update. Note that several other popular integrators do not explicitly include the update of the system Hamiltonian, presumably because they are not concerned with calculating distributions of the work generated by explicit Hamiltonian changes.

The alternative splittings with one force evaluation per time step include ORVRO (a stochastic generalization of position Verlet[28]),RVOVR (an explicit Hamiltonian-update generalization of Leimkuhler and Matthew’s ‘ABOBA’),VRORV (an explicit Hamiltonian-update generalization of Leimkuhler and Matthew’s ‘BAOAB’),ROVOR,and VOROV,

Since for a single time step Δt the error is O(Δt3) for any of these Strang splittings, when applied over N = t/Δt time steps the global error is O(Δt2). Figure 1 confirms that OVRVO errors in the energy are second order in the time step Δt.

Figure 1

Numerical demonstration that the errors in energy of the OVRVO integrator (eq 7) are second order in Δt. Here, we use a previously described model system[19] of a harmonic potential, with unit spring constant, friction coefficient, temperature, and mass, with initial conditions r(0) = v(0) = 0. The error is the absolute deviation of the estimate of ⟨r2(1) + v2(1)⟩ = 0.9796111900... (twice the energy) computed by ensemble averaging over 108 independent realizations. Standard errors of the mean are substantially smaller than the symbol size. The line is the graph of the function Δt2.

Time Step Rescaling Recovers Dynamical Properties

Standard integrators implicitly set the parameter b to unity in 7b, 7c, 7e, and 7f. However, we show later that a nonunit b recovers, for arbitrarily large time step, the continuous-limit values of important dynamical quantities. The time step rescaling can be most simply derived by noting that in the absence of a potential, for any of the six splittings, a trajectory is isomorphic to a semiflexible Gaussian polymer chain:[33] the set of positions correspond to the beads on the polymer chain, and the displacement vectors during single time steps correspond to the intermonomer bond vectors in the chain. In one dimension, the displacement is a normally distributed random variable with zero mean and variance σ02 = (bΔt)2/(βm) (in accordance with Maxwell–Boltzmann velocity statistics and the time interval bΔt), and the autocorrelation between velocities separated by N steps decays exponentially asFor this system, the mean-squared displacement in the large time (γt ≫ 1) limit is[33]

The time step rescaling b results from equating this expression to the mean-squared displacement of a freely diffusing particle in one dimension, 2Dt = 2t/(βmγ), for a total simulation time over N steps, t = N Δt. In particular, the time step used in the position update step is rescaled by the factorensuring that the effective free-particle diffusion constant is independent of time step length (see Figure 2). In the low friction limit γΔt ≪ 1, b = 1 – O([γΔt]2), and in the high friction limit γΔt ≫ 1, b = [2/(γΔt)]1/2.

Figure 2

Time step rescaling recovers correct force-free diffusion as a function of time step. Root mean-squared displacement versus relative time step length Δt at time t = 64 for a freely diffusing particle in one dimension, with unit mass, temperature, and friction coefficient, subject to the OVRVO integrator (eq 7) without time step rescaling, b = 1 (○), or with time step rescaling, eq 15 (×).

Note that even though the position update utilizes an effective time step of bΔt, the simulation clock is still advanced by the full time step Δt. The zero-force MSV and VAC and the uniform-force terminal drift are unaffected by the choice of b. We derive the time step rescaling from a different perspective and in more detail in Supporting Information, “Determination of rescaling parameters.”

Integrator Properties

OVRVO Generalizes Other Popular Integrators

For an explicitly time-independent system Hamiltonian, this family of integrators reduces to various other schemes in certain limits or approximations. At zero friction, γ = 0 and a = b = 1, thus stochastic substeps have no effect, so OVRVO, VRORV, and VOROV are identical to the deterministic velocity Verlet integrator, whereas ORVRO, RVOVR, and ROVOR are identical to position Verlet. In the high-friction or long-time limit, a = 0 and b = [2/(γΔt)]1/2, and OVRVO reduces to the Euler integrator for overdamped Langevin dynamics[34] (also known as the Euler-Maruyama method[35]):ROVOR and VOROV interpose a velocity randomization substep between the deterministic velocity and position updates. In this γΔt ≫ 1 limit, the velocity is completely randomized before each position-update substep, and thus the position updates are completely independent of the Hamiltonian. The other four splittings preserve the influence of the Hamiltonian on the dynamics even in this limit of large friction (or large time step).

OVRVO also reduces to several other popular integrators in other limits or approximations. If the effective time step rescaling for the deterministic substeps is omitted, such that b = 1, then OVRVO is equivalent to an integrator described by Adjanor, Athènes, and Calvo;[18] and by Bussi and Parrinello.[20] If we also combine all stochastic and deterministic velocity updates (eqs 7f, 7g, 7a, and 7b), we recover the integrator of Athènes;[16] and recasting the Athènes integrator as a Verlet-style integrator (only monitoring position) we converge with the Brünger-Brooks-Karplus (BBK) integrator[15] in the low friction limit. If instead we only combine adjacent pairs of stochastic and deterministic velocity updates (eqs 7a with 7b, 7f with 7g) (still with no time step rescaling) we produce the low friction limit of the Langevin Leapfrog integrator of Izaguirre, Sweet, and Pande.[11,22]

Nonequilibrium Work

There is significant interest in probing the probability distribution of work required during a nonequilibrium driving process, which via the work fluctuation theorems[6] can report on various system properties. Such usage requires careful splitting of thermodynamically distinct energy changes.[7]

The total energy change ΔE during the nth full time step of OVRVO can be cleanly separated into heat Q, protocol work Wprot, and shadow work Wshad:[7]Here, (r,n) is the potential energy for configuration r under Hamiltonian (n), and (v) = (1/2)mv2 is the kinetic energy for velocity v. The five other splittings permit similar decompositions of energy changes into heat, protocol work, and shadow work.

Constraints

Constraints, such as rigid bond lengths, can be readily incorporated into the dynamics using standard techniques.[24] The symplectic part of the integrator can be constrained with the standard RATTLE [7b, 7f] algorithm.[3,36] Since RATTLE is symplectic if iterated to convergence,[37] adding constraints does not interfere with the underlying reversibility of the dynamics. Similarly, the velocity randomization substeps, eqs 7a and 7g, can be constrained with RATTLE, which modifies the heat flow, but preserves detailed balance.[24] Consequently, constrained versions of this family of integrators still obey the precepts of nonequilibrium thermodynamics,[7] with the same definitions of heat, protocol work, and shadow work (eq 17a), provided that the definition of free energy is altered to account for the constrained degrees of freedom.[24]

Computational Efficiency

All six splittings require one force evaluation per time step. For OVRVO, for example, the force in substep 7f is identical to the force in substep 7b of the next time interval. Measuring heat requires two evaluations of kinetic energy per time step for all six splittings, for OVRVO just after substep 7a and just before substep 7g. Separately measuring protocol work and shadow work requires two potential energy evaluations per time step, once each just before and after the Hamiltonian-update substep. Shadow work measurement also requires the kinetic energy evaluations already needed to measure heat, as shadow work and heat are the only processes that change the kinetic energy.

However, if only the total thermodynamic work W = ΔE – Q is of interest, this can be calculated given the heat Q( during each step n → n + 1 (easily accumulated during integration using eq 17c) and the total energy at the beginning and end of the simulation,

The OVRVO integrator requires two normal random numbers per velocity per time step, one each for the initial 7a and final 7g velocity randomizations. Splitting the velocity randomization across time steps ensures that the dynamics is microscopically reversible and Markovian, and that the induced Markov chain is irreducible. (The two separate randomization steps permit the independent adjustment of the velocity and the position to arbitrary values.[13]) ORVRO and VOROV also induce irreducible Markov chains. RVOVR, VRORV, and ROVOR, by contrast, do not generate irreducible Markov chains: they effectively agglomerate the two velocity randomizations into a single randomization involving one random number, so for a particular new velocity only one position is possible. Irreducibility has utility for path sampling[10,38−40] and path reweighting[11,12] schemes, since any proposed discrete time trajectory through phase space is a valid trajectory of an irreducible Markov chain. However, many practical applications do not require strict irreducibility, so one can halve the number of required random variables for OVRVO and ORVRO by combining the last stochastic substep of one full step with the first stochastic substep of the next full step, and for RVOVR and VRORV by combining the two stochastic substeps of a given full step. ROVOR and VOROV separate their stochastic substeps such that they cannot be easily combined. Leimkuhler and Matthews show that in the high friction limit and at medium time step VRORV with this single velocity randomization (and time-independent Hamiltonian) is second-order accurate when other integrators become first-order.[25]

When only the total thermodynamic work is of interest, we can combine the last two velocity updates of eq 7 with the first two updates from the next step, and combine the two position updates, to give a three substep stochastic Leapfrog integrator:Under these circumstances, RVOVR, ROVOR, and VOROV reduce to similar three substep integrators, but, due to their sequencing of substeps, ORVRO and VRORV each only reduce to a five substep integrator.

Path Action

The path action S[X] is a necessary quantity for many path sampling[10,38−40] and path reweighting[11,12] techniques. The conditional path probability functional is a product of single time step probabilities,Here, X is a trajectory through phase space between x(0) ≡ {r(0),v(0)} at time 0 and x(NΔt) at NΔt. Each time step probability is determined by the probability of the requisite random variables, which for OVRVO isThe first factor βm/[bΔt(1 – a)] is the Jacobian for the change of variables from {r(n + 1), v(n + 1)} to {+(n), –(n + 1)}, and the probabilities are normal with zero mean and unit variance,whereThe intermediate velocities can be determined by the initial and final positions and forces:Combining eqs 20–24 gives the action as a function of position and velocity at the unit time steps,The path probability obeys the expected symmetry under time-reversal,[7] where the work is defined as in eqs 17d and 17e.

VOROV has a similarly simple expression for the path action. The path action for ORVRO additionally requires the evaluation of the force and its derivative at the half-step position, r(n + (1/2)), hence requires paths of twice as many points, and therefore is of lesser utility. RVOVR, VRORV, and ROVOR induce reducible Markov chains and thus these splittings have infinite path actions for the vast preponderance of paths.

Results

All of our Strang splittings lead to seven-substep integration schemes that are time-symmetric; are second-order accurate in Δt; make Hamiltonian changes explicit; distinguish between heat, protocol work, and shadow work; and easily incorporate constraints. Six of the 12 unique splittings require a single force evaluation per time step and thus are computationally efficient. Setting a ≡ e–γΔ and b ≡ [(2/(γΔt)) tanh (γΔt/2)]1/2 gives for all six one-force-evaluation splittings the continuous-limit MSD, MSV, and VAC in the zero-force case, and the linear-force virial with asymptotic error O(Δt2). The six splittings differ in the remaining desiderata. We summarize their properties in Table 2.

Table 2

Comparison of Properties for Different Splittingsa

desideratum	OVRVO	ORVRO	RVOVR	VRORV	VOROV	ROVOR
All Six Splittings Perform Identically
form is time-reversal symmetric	yes	yes	yes	yes	yes	yes
splits heat, work, and shadow work	yes	yes	yes	yes	yes	yes
easily incorporates constraints	yes	yes	yes	yes	yes	yes
force evaluations per time step	one	one	one	one	one	one
zero-force MSV	exact	exact	exact	exact	exact	exact
zero-force VAC	exact	exact	exact	exact	exact	exact
zero-force MSD	exact	exact	exact	exact	exact	exact
linear-force virial	O(Δt2)	O(Δt2)	O(Δt2)	O(Δt2)	O(Δt2)	O(Δt2)
Splittings Differ in Performance
uniform-force terminal drift	exact	exact	exact	exact	O(Δt2)	O(Δt2)
linear-force MSD	O(Δt2) at nexact at n + 1/2	O(Δt2) at nexact at n + 1/2	exact at n	exact at n	O(Δt2)	O(Δt2)
linear-force MSV	exact at n	exact at n	O(Δt2) at nexact at n + 1/2	O(Δt2) at nexact at n + 1/2	O(Δt4) at n	O(Δt2) at nO(Δt4) at n + 1/2
irreducible Markov chain	yes	yes	no	no	yes	no
path action	simple	requires values atn + 1/2	may be ∞	may be ∞	simple	may be ∞
Hamiltonian dependence for large γΔt	yes	yes	yes	yes	no	no
can halve number of random variates	yes	yes	yes	yes	no	no
generalizes several popular integrators	yes	no	no	no	no	no

Desiderata are grouped into those satisfied by all six splittings and those where the splittings differ in their performance.

Conclusions

As a stochastic generalization of a standard deterministic technique, OVRVO implements what can be considered a form of velocity Verlet with velocity randomization (VVVR). It is an integrator of general utility, satisfying five of Pastor et al.’s seven dynamical properties (six if fractional time step quantities are considered), as well as the remaining enumerated desiderata. It is well-suited for the study of nonequilibrium thermodynamics: work is easily measured because the Hamiltonian changes are made explicit, and these measured works are thermodynamically meaningful because OVRVO distinguishes between heat, protocol work, and shadow work. OVRVO’s simple form for the path action facilitates its use in trajectory reweighting or path sampling methods. Our novel time step rescaling maintains (for an arbitrary time step) various continuous-limit dynamical quantities, in particular the uniform-force terminal drift and linear-force fluctuations in position and velocity. Further research is required to determine to what extent this rescaling permits a larger raw integration time step, and hence how it affects sampling efficiency. Finally, OVRVO generalizes several popular integration schemes and thus relates naturally to the existing literature. Its extension to a multiple time step integrator[41] is straightforward, requiring only replacement of the inner symplectic step (7b–7f) with a corresponding symplectic multiple time step integrator for deterministic dynamics.

By contrast, alternative splittings suffer from shortcomings of varying severity and number. ROVOR and VOROV produce uniform-force terminal drift and linear-force fluctuations of position and velocity that differ from continuous-limit values, lose Hamiltonian dependence in the limit of large γΔt, and require two random numbers per time step even for reducible Markov chains. RVOVR and VRORV induce Markov chains that are not irreducible and thus these splittings have limited utility for path-sampling schemes. ORVRO seems similarly useful to OVRVO, but in path sampling applications requires paths with twice the number of points as OVRVO.