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 . 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 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 
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 .
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 , which represents interactions 
 with the surrounding environment through frictional drag and stochastic 
 collisions with fictitious bath particles:1a1bhere 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 ,
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 .
finite time step integrators 
 for molecular 
 systems can not 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 .
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 . 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 
we describe in detail the calculation of these quantities for our 
 family of integrators .
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 .
angled brackets denote an average 
 over the ensemble of phase - space trajectories produced by a given 
 integration scheme .
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 or path reweighting 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 ) , 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 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 .
for the calculated 
 works to be thermodynamically meaningful , the integrator must also 
 have a form symmetric under time - reversal .
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 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 .
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 .
we derive a family of integrators by splitting the time evolution 
 operator into stochastic and deterministic components 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 )
can 
 be naturally written as a sum of four parts   = o + v + r + h. the first 
 operator represents stochastic thermalization via an ornstein  uehlenbeck operator,2the next two operators 
 represent deterministic 
 newtonian evolution of velocity and position,3and the last operator represents the time 
 evolution of the system hamiltonian according to the predetermined 
 schedule ( or protocol ) ,4where n denotes the time 
 step index and t = nt 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 , we approximate the dynamics over a time 
 t by applying a series of strang ( symmetric 
 trotter ) operator splittings:5a5b5c5dwhere ( a , b , c , d ) represents a permutation of 0 , v , r , h. there are six strang splittings of the liouville 
 operators 0 , r , and v. the hamiltonian 
 update operator h commutes with 0 and with r , so for each 
 of these six splittings , there are only two unique placements of h. for each 
 of the six splittings , one placement of h 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 ( t ) , 
 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 , 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 
 by7a7b7c7d7e7f7ghere , a = exp(t ) , and   and   are independent , normally 
 distributed random variables with zero mean and unit variance ( hence , 
 when scaled by ( m ) , 
 distributed according to the equilibrium maxwell 
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 ( 7 g ) substeps that are stochastic , markovian , and detailed - balanced 
 ( with respect to the true canonical measure ) .
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 , 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 , 
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),8rvovr ( an 
 explicit hamiltonian - update 
 generalization of leimkuhler and matthew s  aboba),9vrorv ( an explicit hamiltonian - update 
 generalization of leimkuhler and matthew s 
baoab),10rovor,11and vorov,12 since for a single time step t the 
 error 
 is o(t ) for any 
 of these strang splittings , when applied over n = t/t time steps the global error 
 is o(t ) .
figure 1 confirms that ovrvo errors in the energy are second 
 order in the time step t .
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 : 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 = ( bt)/(m ) ( in accordance with maxwell 
boltzmann 
 velocity statistics and the time interval bt ) , and the autocorrelation between velocities separated 
 by n steps decays exponentially as13for this system , the mean - squared displacement 
 in the large time ( t  1 ) limit is 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 factor15ensuring 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 ] ) , and in the high 
 friction limit t  1 , b = [ 2/(t ) ] .
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 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 r(1 ) + v(1) = 0.9796111900 ... 
 ( twice the energy ) computed by ensemble averaging over 10 independent realizations .
note that even though the position update utilizes an effective 
 time step of bt , 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 , 
 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 (  ) .
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 ) ] , and ovrvo reduces to the euler integrator for overdamped langevin 
 dynamics ( also known as the euler - maruyama 
 method):16rovor 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 , athnes , 
 and calvo ; and by bussi and parrinello .
if we also combine all stochastic and deterministic 
 velocity updates ( eqs 7f , 7 g , 7a , and 7b ) , we recover 
 the integrator of athnes ; and 
 recasting the athnes integrator as a verlet - style integrator 
 ( only monitoring position ) we converge with the brnger - brooks - karplus 
 ( bbk ) integrator in the low friction 
 limit .
if instead we only combine adjacent pairs of stochastic and 
 deterministic velocity updates ( eqs 7a with 7b , 7f with 7 g ) ( still with no time step rescaling ) we produce the low friction 
 limit of the langevin leapfrog integrator of izaguirre , sweet , and 
 pande .
there is significant 
 interest in 
 probing the probability distribution of work required during a nonequilibrium 
 driving process , which via the work fluctuation theorems can report on various system properties .
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:17c17d17ehere , ( r , n ) 
 is the potential energy for configuration r under 
 hamiltonian ( n ) , and ( v ) = ( 1/2)mv 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 , such as rigid bond lengths , 
 can be readily incorporated into the dynamics using standard techniques .
the symplectic part of the integrator can be 
 constrained with the standard rattle [ 7b , 7f ] algorithm .
since rattle is symplectic 
 if iterated to convergence , adding constraints 
 does not interfere with the underlying reversibility of the dynamics . 

similarly , the velocity randomization substeps , eqs 7a and 7 g , can be constrained with rattle , 
 which modifies the heat flow , but preserves detailed balance .
consequently , constrained versions of this family 
 of integrators still obey the precepts of nonequilibrium thermodynamics , 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 .
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 7 g . 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,18 the ovrvo integrator requires two normal random numbers per 
 velocity 
 per time step , one each for the initial 7a and 
 final 7 g 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 . )
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 and path reweighting 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 can not 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 .
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 .
the path action s[x ] is a necessary quantity for many path 
 sampling and path reweighting 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(nt ) at nt .
each time step probability is determined by the probability 
 of the requisite random variables , which for ovrvo is21the first factor m/[bt(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,22wherethe intermediate velocities can 
 be determined 
 by the initial and final positions and forces:24a24bcombining eqs 2024 
 gives the action as a function of position and velocity at the unit 
 time steps,25the path probability obeys the expected symmetry 
 under time - reversal , where the work is 
 defined as in eqs 17d and 17e .
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 .
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 ) ]
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(t ) .
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 .
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 is straightforward , requiring 
 only replacement of the inner symplectic step ( 7b7f ) with a corresponding symplectic 
 multiple time step integrator for deterministic dynamics .
desiderata are 
 grouped into those 
 satisfied by all six splittings and those where the splittings differ 
 in their performance . by 
 contrast
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 .
 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
Theory
Integrator Splitting
Time Step Rescaling
Recovers Dynamical Properties
Integrator Properties
Results
Conclusions