Using nonequilibrium fluctuation theorems to understand and correct errors in equilibrium and nonequilibrium discrete Langevin dynamics simulations

Common algorithms for computationally simulating Langevin dynamics must discretize the stochastic differential equations of motion. These resulting finite time step integrators necessarily have several practical issues in common: Microscopic reversibility is violated, the sampled stationary distribution differs from the desired equilibrium distribution, and the work accumulated in nonequilibrium simulations is not directly usable in estimators based on nonequilibrium work theorems. Here, we show that even with a time-independent Hamiltonian, finite time step Langevin integrators can be thought of as a driven, nonequilibrium physical process. Once an appropriate work-like quantity is defined -- here called the shadow work -- recently developed nonequilibrium fluctuation theorems can be used to measure or correct for the errors introduced by the use of finite time steps. In particular, we demonstrate that amending estimators based on nonequilibrium work theorems to include this shadow work removes the time step dependent error from estimates of free energies. We also quantify, for the first time, the magnitude of deviations between the sampled stationary distribution and the desired equilibrium distribution for equilibrium Langevin simulations of solvated systems of varying size. While these deviations can be large, they can be eliminated altogether by Metropolization or greatly diminished by small reductions in the time step. Through this connection with driven processes, further developments in nonequilibrium fluctuation theorems can provide additional analytical tools for dealing with errors in finite time step integrators.


I. INTRODUCTION
In the computational natural sciences, dynamic properties of a stochastic system are often calculated using simple numerical integrators for Langevin dynamics [1], where the system is driven from equilibrium by a timedependent Hamiltonian H(t). In the simplest case of a single stochastic particle, r and v are time-dependent position and velocity, m is mass, f is force, β = 1/k B T , k B 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 is determined by the derivative of the potential energy, f ≡ − ∂H/∂r. For multidimensional, multiparticle systems, r, v, f , and d W are vectors, and m is a diagonal matrix. In order to simulate Langevin dynamics on a digital computer, it is necessary to adopt some approximate algorithm that divides time into discrete steps [2]. However, most such schemes have an inherent problem: Even * david.sivak@ucsf.edu; Current address: Center for Systems and Synthetic Biology, University of California, San Francisco, California 94158, USA with a time-independent Hamiltonian, they do not preserve the canonical equilibrium distribution determined by H nor do they satisfy microscopic reversibility. (By reversibility we mean that the probability of sampling a particular trajectory starting from equilibrium is equal to the probability of sampling the trajectory's time reversal, reversing velocities if necessary.) We show that these pathologies arise because discrete time step integrators of Langevin dynamics can be viewed as simulations of artificial driven nonequilibrium dynamics. This perspective has the advantage that the complications generated by this unwanted but inevitable breaking of timereversal symmetry can be understood and remedied in a controlled and systematic fashion with insights from nonequilibrium statistical thermodynamics [3][4][5][6].
We can appreciate some of the problems inherent in finite time step Langevin dynamics by first considering the zero friction limit, γ = 0, with a time-independent Hamiltonian, where Langevin dynamics reduces to deterministic Newtonian dynamics. A simple, popular integrator for Newtonian dynamics is the velocity Verlet algorithm [7,8], Because of the finite time step, the trajectories generated by this algorithm are inaccurate: They do not faithfully follow the precepts of Newtonian mechanics. Also, the actual energy of the system is not conserved, but rather it fluctuates from one time step to the next. However, the velocity Verlet integration scheme is symplectic (in that the Jacobian of the transformation from old to new positions and velocities is unity, and therefore the phasespace volume is conserved [9]), which ameliorates some problems due to the finite time step. For example, although a finite time step symplectic integrator does not conserve the energy of the system Hamiltonian, it does conserve the energy of a shadow Hamiltonian, which is close to the desired Hamiltonian if the time step is not too large [2,10]. For sufficiently small timesteps, this conservation of the shadow Hamiltonian prevents longterm drift in the system Hamiltonian over the duration of the simulation. Essentially, a finite time step dynamics performs work on the system, over-and-above any work due to intentional perturbations from a time-dependent Hamiltonian [6]. We can imagine this finite time step integration scheme in the following way. At the beginning of each time step, we first perturb the system Hamiltonian such that it becomes the shadow Hamiltonian, changing the energy of the system. The symplectic integrator then updates the position and velocity [Eq. (2)], perfectly preserving the shadow energy of the shadow Hamiltonian. We then switch the Hamiltonian back to the original one, again perturbing the energy. The net change in the energy of the system during this time step is due to work performed on the system by perturbing back and forth between the system and shadow Hamiltonian. We can determine this shadow work (also known as error work [6] or an effective energy change [11]) during each time step by measuring the difference in energy using the system Hamiltonian, so we do not need to know the form of the shadow Hamiltonian. This shadow work is distinct from any protocol work applied to the system due to explicit, time-dependent perturbations of the system Hamiltonian. Note that Markov-chain Monte Carlo (MCMC) simulations do not generate shadow work [12] because the dynamics satisfies detailed balance explicitly, which ensures that the trajectories are microscopically reversible [13] and that the appropriate equilibrium ensemble is preserved for a time-independent Hamiltonian [2].
Discretizations of continuous-time Langevin dynamics are essentially a combination of deterministic and stochastic dynamics, and, as a result, they suffer from a combination of problems. With a finite time step, the deterministic parts of the dynamics tend to pump energy into the system in the form of shadow work, driving the system away from equilibrium, whereas the stochastic parts of the dynamics relax the velocities back toward the equilibrium Maxwell-Boltzmann distribution, removing energy from the system in the form of heat. It follows that, even for a system with a Hamiltonian that is explicitly time-independent, a finite-time-step Langevin dynamics has an effective Hamiltonian alternating between the system Hamiltonian and the shadow Hamilto-nian, and thus actually simulates a driven, nonequilibrium system, with a net energy flow. Microscopic timereversal symmetry is broken, and in general we can not determine the steady-state, nonequilibrium distribution. These difficulties may be circumvented by reducing the time step, but at the cost of increasing the computational effort required to simulate the same interval of time; this hardly constitutes a satisfactory resolution.
The main point of this paper is this interpretation of the errors induced by discrete simulation of Langevin dynamics, in terms of a driven thermodynamic process. This perspective forms a bridge between the study of numerical integrators and the rapidly expanding field of nonequilibrium statistical mechanics, permitting the invocation of a wide array of nonequilibrium work fluctuation relations to characterize and correct for biases in estimates of equilibrium and nonequilibrium thermodynamic quantities.

II. CONCRETE INTEGRATOR
We demonstrate the utility of this perspective for an integration scheme that is explicitly time-symmetric, that cleanly separates the stochastic and deterministic parts of the dynamics, and for which the deterministic parts are symplectic and the stochastic parts are detailed balanced. This construction allows a clean separation of the system's energy change into work, shadow work, and heat, simplifying our analysis in terms of a driven nonequilibrium process. Fortunately, integrators with these properties have received recent attention [4,5,[14][15][16]. As a concrete example, we consider the integrator used by Bussi and Parrinello [11], where we make the Hamiltonian update explicit: Here, ∆t is the time step by which the simulation clock is advanced, f (n) is the force at position r(n) due to the Hamiltonian H(n), a = exp(−γ ∆t), and N + and N − 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 . The stochastic substep (3a) randomizes the velocity, transferring heat between the system and environment, while the Hamiltonian is fixed and the position unchanged. We then switch from the system to the shadow Hamiltonian, performing shadow work on the system. Substeps (3b) and (3c) update the velocity and then the position according to the symplectic dynamics of the shadow Hamiltonian, exactly conserving the energy. We next switch back to the system Hamiltonian (performing shadow work), and in (3d) update the system Hamiltonian from H(n) to H(n + 1), according to the prescribed protocol Λ. This action performs protocol work on the system. We switch back to the shadow Hamiltonian (doing shadow work), symplectically update position and then velocity (3e,3f), and then restore the system Hamiltonian (again performing shadow work). Finally, we conclude with another velocity-randomization substep (3g). first and last substeps (3a,3g) are stochastic, Markovian, and detailed-balanced (with respect to the canonical measure) velocity randomizations, which leave the position unchanged. The five middle substeps (3b-3f) constitute the deterministic velocity Verlet integrator (2), with the midpoint Hamiltonian update made explicit. The order of substeps and the effective Hamiltonian switches are illustrated in Fig. 1. Note that the deterministic substeps (3b,3c,3e,3f) are each individually symplectic.

III. NONEQUILIBRIUM THERMODYNAMICS
A central relation of driven, nonequilibrium thermodynamics [17][18][19][20] relates the microscopic irreversibility of trajectories to the work W [X, Λ] performed on the system during the forward protocol [21][22][23]: Here, X is a trajectory through phase space between time 0 and N ∆t, Λ represents a protocol for perturbing the system (typically through the time dependence of the system Hamiltonian), ∆F eq [Λ] is the free energy difference between the equilibrium distributions for the initial and final values of the system Hamiltonian, and P X Λ is the probability of the trajectory, given the protocol and an initial equilibrium ensemble. The time-reversed pro-tocolΛ (time-reversed trajectoryX) retraces the same series of perturbations (phase-space transitions) as the forward protocol Λ (forward trajectory X), but under time inversion and hence in reverse. Subject to a protocol, a driven system is microscopically reversible if the probability of a trajectory and its time reversal are identical, and therefore the work imposed by the protocol equals the free energy change [24]. It is straightforward to extend this fluctuation theorem to mixed stochastic-deterministic dynamics, such as the Langevin integrator, Eq. (3), provided that the individual substeps satisfy this symmetry. It is for this reason that we insist on a clean separation of the deterministic and stochastic substeps.
The total work W = n W (n) is the sum of the contributions W (n) from individual steps. The total change in energy ∆E during the step n → n + 1 can be cleanly separated into heat Q, protocol work W prot , and shadow work W shad : Here, ∆E a-g are the energy changes during the corresponding substeps of Eq. (3). Heat is the energy exchanged with the thermal environment, protocol work is the energy change due to deliberate manipulation of the Hamiltonian (i.e., the explicit time-dependence of the system Hamiltonian), and shadow work is the energy change due to alternation between the system and shadow Hamiltonians, resulting from the finite time step of the symplectic part of the integrator. The essential distinction between heat and work is that heat flow is change of the system energy due to change in the current distribution over microstates, whereas work is change of energy due to change in the equilibrium distribution over microstates.
The stochastic velocity randomization substeps obey Eq. (4) since they are balanced, in that they preserve the canonical equilibrium distribution [21]. The set of deterministic velocity Verlet substeps also obeys Eq. (4), so long as the total work includes the shadow work [3,6], since the dynamics is symplectic and microscopically reversible with respect to the shadow Hamiltonian [2,10]. Since both the deterministic and stochastic substeps are Markovian, it follows that we can safely intermix the two dynamics, and (4) still holds.
It therefore follows that the Langevin integrator obeys various derived relations of nonequilibrium statistical dynamics, such as the Jarzynski equality [25], fluctuation relations [21,26], interrelations between path ensemble averages [22,27] and various interrelations between dissipation and time asymmetry [28][29][30][31]. Furthermore, by its separation of protocol work and shadow work, the Langevin integrator permits the separation of the respective contributions to microscopic irreversibility of deliberate perturbation (physically meaningful) and the finite time step (a discretization artifact). Notably, the statistics of the protocol work alone systematically deviate from those of the total work, and hence lead to biased inference when using the machinery of nonequilibrium thermodynamics. In Secs. VI and VII we explicitly demonstrate this underappreciated point.

IV. "EQUILIBRIUM" SIMULATIONS SAMPLE PERTURBED DISTRIBUTIONS
It is common practice in the study of the equilibrium properties of molecular systems to use a single finitetime-step mixed stochastic and deterministic dynamical simulation to sample from an equilibrium distribution. However, this distribution departs from the true equilibrium distribution for the system Hamiltonian, a distribution that we can now understand as the steady state due to driving by the finite time step. Thus a question of significant practical interest presents itself: How far from equilibrium is the effective nonequilibrium steady state induced by this time discretization for a system with a time-independent Hamiltonian? Since the explicit system Hamiltonian is unchanging, no protocol work is performed, and thus our analysis in this section focuses on the shadow work alone. Practitioners commonly estimate artifactual errors by monitoring some essentially arbitrary, yet easily measured, observable of the system, such as the total energy. However, we can exploit recent advances in nonequilibrium statistical dynamics to provide a principled characterization of how far the system is driven from equilibrium [32].
The natural measure of this instantaneous distance that the system has been driven away from equilibrium is the difference between a nonequilibrium free energy [33,34] F neq ≡ E −T S and the corresponding equilibrium free energy F eq for the given Hamiltonian H. If the Hamiltonian were held constant and the (previously driven) system were allowed to relax to equilibrium, this deviation from the equilibrium free energy would represent the heat that would be lost to the environment, or equivalently the maximum work that could be imparted to a mechanically coupled system. For the perturbations imposed by the discrete dynamics, this nonequilibrium free-energy deviation is approximated near equilibrium by [32] ∆F neq ≡ F neq − F eq ≈ 1 2 W shad − (t f − t i )P ss , (6) where W shad is the shadow work over the whole simulation, P ss is the power (work per unit of time) once transients have died off and the system has settled into a nonequilibrium steady state, and t f − t i is the total simulation time. Normalizing this nonequilibrium free-energy deviation by the size of the system (number of degrees of freedom) provides a natural measure of how far from equilibrium each degree of freedom is on average.
To estimate the nonequilibrium steady-state freeenergy deviation for a molecular system, we simulate cubic boxes of TIP3P waters of various sizes, both with and without constraints on the water O-H and H-H interatomic distances. (See the Appendix for simulation details.) Initial coordinates and momenta are sampled from equilibrium in an isothermal-isobaric (NPT) ensemble (that is, an ensemble that maintains constant number of particles, constant pressure, and constant temperature) at 1 atm and 298 K using the generalized hybrid Monte Carlo (GHMC) integrator [16,35]. These initial conditions are simulated for M steps with the Langevin integrator [Eq. (3)] at constant volume (using a collision rate γ = 9.1/ps) to measure the nonequilibrium work to reach steady state, followed by an additional M steps to measure the steady-state power. We have determined that, for all systems and time steps simulated, M = 1028 steps is sufficient to reach steady state (see Fig. 4). We have also calculated the statistical uncertainty according to Eq. (A.2).
Because the system (a periodic water box) is homogeneous, it is possible to collapse all system sizes onto universal curves describing the nonequilibrium freeenergy deviation per molecule as a function of time step for unconstrained and constrained systems, respectively (Fig. 2). For the unconstrained system, whose numerical integration becomes unstable beyond ∆t = 1.5 fs, the nonequilibrium free-energy deviation ∆F neq rapidly rises as the time step surpasses the typical time step employed for flexible systems, ∆t ≈ 1 fs. For a system of 220 waters, for example, ∆F neq = 11.4 ± 0.2 k B T at ∆t = 1 fs. For constrained water boxes, however, ∆F neq reaches this magnitude only at large time steps-here, ∆t ≈ 5 fs, not far from the stability limit at 6 fs and well beyond 2 fs, the standard time step for biomolecular simulations.
Empirically, the nonequilibrium free-energy deviation ∆F neq for both unconstrained and constrained systems appears to show a quartic dependence on the time step ∆t (Fig. 2, gray curves), such that where the prefactor a depends strongly on whether constraints are employed; see the caption of Fig. 2. This trend is consistent with earlier work observing the strong dependence of Metropolization acceptance probabilities on time step [36] and highlights how small reductions in time step can rapidly reduce the deviations of the sampled steady-state distribution from the desired equilibrium distribution defined by the system Hamiltonian p eq (x) ∝ exp[−βH(x)], without unduly burdensome computational cost. We detail in Sec. VI some methods that correct for these nonequilibrium perturbations. Even in the absence of correction procedures, the above calculation represents a thermodynamically meaningful determination of the deviation from the desired equilibrium sampling associated with the continuous Langevin equation of motion, as a function of simulation parameters.

V. MULTIVARIATE FLUCTUATION THEOREM
We seek an analytical framework that describes the correlation between the shadow work (performed by integration) and the protocol work (due to explicit Hamiltonian changes). We want this framework to provide a generic method to characterize the effect that shadow work has on the distribution of protocol work, and specifically on the time-reversal symmetry [Eq. (4)] that protocol work would satisfy in its absence. Furthermore, we want this framework to suggest systematic techniques to correct for these distorting effects. We propose such a framework through the generalization of work fluctuation theorems to the context of two sources of work. These results, although formulated specifically for our situation of explicit and artifactual work, are entirely general to situations involving any two sources of work.
Rearrangement of Eq. (4) and splitting the work into two distinct work contributions W 1 , W 2 gives Multiplication by delta functions of the two works, δ(W 1 [X, Λ] − W prot )δ(W 2 [X, Λ] − W shad ), and integration over all trajectories produces what we refer to as the multivariate fluctuation theorem, This is a special case of the generalized detailed fluctuation theorem for joint probabilities of Garcìa-Garcìa, et al. [37,38]. Equation (9) gives an expression in terms of the excess work W prot + W shad − ∆F eq for the ratio of the joint probability distributions over protocol and shadow works realized during the forward and reverse protocols, respectively. Equation (9) can be trivially extended to arbitrary decompositions of the total work, where each component corresponds to a group of individual work steps. It thus represents a generalization of the work fluctuation theorem [22] to contexts with multiple sources of work. From Eq. (9), several other modified fluctuation theorems can be derived that modify a standard fluctuation theorem for one of the works with an exponential average over the other work. For example, in Sec. VI, we derive a Jarzynski equation modified by the presence of shadow work [Eq. (12)], and, in Sec. VII, we derive a similarly modified integrated transient fluctuation theorem [Eq. (15)].

VI. RECOVERING EQUILIBRIUM STATISTICS FROM NONEQUILIBRIUM SIMULATIONS
Now that we are equipped with our new interpretation of finite-time-step Langevin dynamics as a driven nonequilibrium process even in the absence of an explicit driving force, nonequilibrium thermodynamics affords various approaches for recovering true equilibrium properties of the system.
One approach is to maintain the simulation at equilibrium by incorporating Monte Carlo moves that conditionally accept or reject candidate trajectory segments or single time steps, for example by using the Metropolis criterion P accept = min(1, exp{−βW shad }) [39]. In order to maintain detailed balance, the velocity must be inverted if the proposed state is rejected [12], which may lead to increased correlation times. Applied to single time steps, this is essentially the idea behind the GHMC integrator [12,35], and when applied to trajectory segments, this is the idea behind work-bias Monte Carlo [40] and nonequilibrium candidate Monte Carlo [41] simulations. In either case, Metropolization results in an MCMC process that samples the true equilibrium distribution.
Another approach to recovering equilibrium statistics is to perform a Monte Carlo sampling of trajectories [42,43], generating an ensemble of trajectories weighted by the Boltzmann-weighted work over the entire trajectory, exp{−βW shad }. This approach allows both accurate equilibrium statistics and realistic dynamics, albeit at a potentially high computational cost.
Instead of sampling equilibrium trajectories, we can alternatively apply nonequilibrium relations, such as the Jarzynski equality [25] and path ensemble averages [22,27,44,45], to directly recover equilibrium properties from the statistics of a driven system, essentially by reweighting trajectories by exp{−βW tot }, where it is important that the work includes both the protocol work and the shadow work. Note that the initial configurations must be sampled from the correct equilibrium ensemble, which can be accomplished with a standard MCMC process, or with one of the approaches discussed above, such as GHMC simulation.
We now demonstrate the importance of including the shadow work by using the Jarzynski equality to estimate free energy changes in a simple model system. The Jarzynski equality [25] relates the equilibrium free energy change, resulting from some perturbation of the system, to the exponential average of the work incurred during many realizations of the system response to that perturbation, (10b) In the second line, we have explicitly split the effective thermodynamic work into protocol and shadow work.
Here, angled brackets with subscript Λ indicate expectations over trajectories starting in the equilibrium distribution for the initial value of the Hamiltonian H(0) and integrated according to Eq. (3), with the Hamiltonian evolving according to Λ. Although standard Langevin integrators are used in myriad multidimensional contexts, we examine in Fig. 3 the shadow work contribution in a simple one-dimensional system to suggest the ubiquity of the issues raised here. In particular, we consider a particle in thermal contact with the environment, subject to a quartic potential that is initially stationary and then translated at a constant velocity. The exact free energy change is zero. When one uses only the protocol work (neglecting the shadow work), the Jarzynski free energy estimate empirically shows a systematic error that scales roughly as ∆t 2 [Figs. 3a,b, circles]. Using the total thermodynamic work (including the shadow work) eliminates this error, and the Jarzynski estimator gives the correct free energy change [ Figs. 3a,b, ×s]. In Fig. 3, standard errors are calculated from 10 8 independent simulations and are smaller than the symbol size. The y axis is the same in the left and right sub-figures. We can understand the origin of this error by analyzing our estimator in terms of the multivariate fluctuation theorem [Eq. (9)] derived above in Sec. V. Rearranging Eq. (9), decomposing the joint probability into the marginal and conditional probabilities, and integrating over the shadow work, we find that when ignoring the contributions of shadow work, the Jarzynski estimator of the free energy β ∆F eq ≡ − ln e −βWprot Λ has a systematic bias from the true free energy change β∆F eq that is a function of the distribution of shadow works: β ∆F eq = β∆F eq − ln e −βW shad Λ .
Empirically, the correction term − ln e −βW shad Λ [Figs. 3a,b, + signs] reproduces the error in the Jarzynski estimator without shadow work, β ∆F eq . The correction factor γ ≡ e −βW shad Λ is analogous to the correction factor that appears in the Jarzynski equality with feedback [46]. Curiously, the correction to the Jarzynski estimator is solely a function of the shadow work distribution, and, in particular, does not explicitly depend on correlations between the shadow work and the protocol work.

VII. CORRECTING NONEQUILIBRIUM FLUCTUATION THEOREMS
In addition to these errors for equilibrium estimators during simulations with an explicitly time-independent Hamiltonian, ignoring the contribution of shadow work leads to systematic errors in estimates of nonequilibrium quantities when the Hamiltonian is explicitly time dependent: The simulated system is actually subject to a different Hamiltonian than the system one, and thus the probability distribution of protocol works does not obey the relevant time-reversal symmetry (4). We quantitate this time-reversal asymmetry by examining violations of the integrated transient fluctuation theorem (ITFT) [47], which for time-symmetric protocols relates the ratio of the probabilities of realizing a negative and a positive total work, respectively, to the exponentially-weighted total work, conditional on the total work being positive: This relation follows directly from Eq. (4). Manipulating Eq. (9) to a similar form produces 4 (x − xmin) 4 uniformly translating with velocity 1/2, starting from equilibrium, with unit temperature, mass, spring constant, and friction coefficient. Standard errors are smaller than symbol size. (a) Error in free energy calculated from the Jarzynski equality (10a) as a function of position of the quartic potential, neglecting shadow work (circles) and including shadow work (×s), for ∆t of 1/4 (red) and 1/8 (blue). The exact free energy change is zero. The error in the naive Jarzynski estimator (circles) is entirely captured by the correction term − ln e −βW shad Λ (+ signs) from Eq. (12), as can be seen by the agreement of these symbols to within statistical error. (b) Semilog plot of error in the Jarzynski free energy estimate after the quartic potential has moved to r = 2.5, as a function of time step length, neglecting shadow work (circles) and including shadow work (×s). Also shown is the correction term from Eq.  13)], for ∆t of 1/4 (red) and 1/8 (blue). The error in the naive ITFT ratio is entirely captured by the correction factor e −βW tot W prot >0 / e −βW prot W prot >0 (+ signs), as can be seen by the agreement of these symbols to within statistical error. (d) Semilog plot of the ITFT ratio after the quartic potential has moved to r = 2.5, as a function of time step length ∆t, neglecting shadow work (circles) and including shadow work (×s). Also shown is the correction factor e −βW tot W prot >0 / e −βW prot W prot >0 (+ signs).
For this relation to hold, the work in the exponential must be the total work, not the protocol work that appears elsewhere in the equation. When one ignores the shadow work and measures only the protocol work, the ratio of the left-hand side and right-hand side, departs from unity to the extent that the protocol work fluctuations do not obey the relevant time-reversal symmetry that the total work fluctuations do. Departure from unity in Eq. (15) quantifies the violation of the nonequilibrium time-reversal symmetry obeyed by a proper thermodynamic work encompassing all energy changes not related to heat.
Figures 3c,d show that, for the simple system described in Sec. VI, the protocol work alone (circles) does not obey the nonequilibrium fluctuation relation required of a thermodynamic work (with an error that empirically scales with the square of the time step), but the sum of the protocol and shadow works (×s) does obey it. The correction factor e −βWtot Wprot>0 / e −βWprot Wprot>0 (+ signs) reproduces the error in the ITFT ratio neglecting shadow work. Thus, ignoring the shadow work and using the protocol work rather than the total work produces systematic biases in estimators of nonequilibrium quantities (such as the nonequilibrium free energy [32] or the nonequilibrium energetic efficiency [48]).

VIII. EPILOGUE
For Hamiltonian dynamics, a finite-time-step symplectic integrator conserves a shadow Hamiltonian and is microscopically reversible. But, as we have seen, for Langevin dynamics, discretization of the dynamics leads (even for a time-independent Hamiltonian) to a mixed deterministic-stochastic nonequilibrium dynamics, which preserves the equilibrium distribution of neither the system nor the shadow Hamiltonian and which is not timereversal symmetric. However, we can measure the work, heat, and shadow work, and thereby separate the respective contributions to time-reversal symmetry breaking of the finite time step and deliberate perturbation. This procedure allows us to apply results from nonequilibrium thermodynamics to characterize in a thermodynamically meaningful way the error produced by finite-time-step integration and to correct for such errors to recover equilibrium and nonequilibrium properties of the system. While we focus in this paper on work distributions, we note that discrete integrators can also introduce artifacts into other aspects of a system's dynamical evolution, for example, producing erroneous free-particle diffusion coefficients and uniform force-field terminal drifts. These artifacts can be mitigated through time-step rescaling, as discussed in Ref. [49]. Where measurements of work and heat are not required, correct statistics of nonequilibrium trajectories through phase space can be recovered using the Metropolis-adjusted geometric Langevin algorithm of Bou-Rabee and Vanden-Eijnden, which under reasonable conditions on the potential energy is pathwise convergent to the distribution of trajectories for the continuous equations of motion [50].
to Langevin simulation [Eq. (3)] at fixed volume using a collision rate of 9.1/ps. We integrated these initial conditions for a total of 4096 steps using a variety of different time steps from 0.25 fs to 7 fs, with the accumulated shadow work after 2 n steps stored (n = 0, 1, . . . , 12). The limit of stability was determined by the largest time step that did not generate infinite cumulative work values in 4096 time steps in any sample; we determined the limit to be 2 fs for unconstrained simulations and 6 fs for constrained simulations.
To estimate, using Eq. (6), the nonequilibrium free energy of the steady-state ensemble sampled by discrete Langevin integration, we used the average accumulated shadow work after M steps as the work to switch into steady state, while we used the average dissipated power in the next M steps as an average steady-state power: Here, the · GHMC notation denotes averages computed over Langevin simulations initiated from GHMC-sampled initial configurations and momenta. Through analysis of M = 2 n for n = 0, 1, . . . , 11, we found that the steadystate power, and hence the estimated nonequilibrium free energy, converged after M = 1024 steps (see Fig. 4), so we used this value for all subsequent analysis.
We estimated the squared uncertainty in the nonequilibrium free energy as