{'fraction_non_alphanumeric': 0.03261216669327798, 'fraction_numerical': 0.00826281334823567, 'mean_word_length': 4.666817462112644, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Density-Dependent Analysis of Nonequilibrium Paths Improves Free Energy Estimates

When a system is driven out of equilibrium by a time-dependent protocol that modifies the Hamiltonian, it follows a nonequilibrium path. Samples of these paths can be used in nonequilibrium work theorems to estimate equilibrium quantities, such as free energy differences. Here, we consider analyzing paths generated with one protocol using another one. It is posited that analysis protocols which minimize the lag, the difference between the nonequilibrium and the instantaneous equilibrium densities, will reduce the dissipation of reprocessed trajectories and lead to better free energy estimates. Indeed, when minimal lag analysis protocols based on exactly soluble propagators or relative entropies are applied to several test cases, substantial gains in the accuracy and precision of estimated free energy differences are observed.


I. INTRODUCTION
The accurate and efficient estimation of free energy differences is an important goal in chemical physics and remains an active area of research. One promising approach to free energy estimation entails measuring the work done on a system over repetitions of an irreversible process. According to the second law of thermodynamics, the mean work is greater than the free energy difference between the end states of the process, F Λ . Nonequilibrium work theorems 1,2,3,4,5 supplement this upper bound by rigorously equating F Λ with other averaged functions of the work. These theorems have been empirically validated in single-molecule pulling experiments 6,7 and computer simulations (e.g. Ref. 8 ).
Jarzynski's equality, 1,2 a unidirectional nonequilibrium work theorem, relates the free energy difference to an exponential average of the work. Unfortunately, because it uses a nonlinear (specifically, a logarithmic) function of the average, the free energy estimator based on this equality suffers from a systematic finite-sampling bias. 9,10,11 While accurate in the limit of infinite sampling, this estimator is usually dominated by rare events where the work is less than the free energy difference, and thereby converges slowly. 12 If the average amount of work dissipated as heat is reduced, these low-work events will be more frequent and accurate free energy estimation will usually require fewer work samples. The most straightforward way to reduce heat dissipation is to slow the rate of the process; in the limit of infinitely slow switching, the process is reversible and the work is equal to the free energy difference. Unfortunately, reducing the switching rate requires additional time and lowers the signalto-noise ratio in single-molecule pulling experiments. 13 Under the constraint of constant experiment length, it is possible to reduce heat dissipation by optimizing the switching protocol that controls how the thermodynamic state changes with time. Protocol variation predates Jarzynski's equality, having been applied to tightening free energy bounds from the second law of thermodynamics. 14,15,16,17,18 More recently, variational calculus has been applied to find optimal protocols that minimize the mean work. 19,20,21 While protocol variation is, in principle, feasible in laboratory experiments, many more approaches to improving nonequilibrium-based free energy estimation are possible in computer simulations. For example, Wu and Kofke were inspired by the Rosenbluth chain sampling scheme to develop methods for generating low-work nonequilibrium paths. 22 Vaikuntanathan and Jarzynski took another approach, altering the system dynamics, to reduce heat dissipation and improve free energy estimates. 23 The approach most mathematically similar to this work, however, is importance sampling in nonequilibrium path space. 24,25,26,27,28 In importance sampling, samples from one distribution are used to estimate expectations in another. The technique is often applied to Markov chain Monte Carlo and molecular dynamics simulations (where it is usually called umbrella sampling): after applying a configurational bias to overcome energy barriers and promote ergodicity, expectations are calculated for the unbiased ensemble. Importance sampling has been extended to transition path sampling 29,30 with nonequilibrium trajectories. In this algorithm, a biasing function modifies the Monte Carlo acceptance criteria of proposed paths in a way that improves the convergence of free energy estimates. 24,25,26,27,28 Here, we apply the importance sampling formalism in a completely different way. Instead of sampling nonequilibrium trajectories in a biased manner, we focus on the analysis of previously generated paths. Instead of asking which path-ensemble we would like to sample from, we ask which path-ensemble average we would like to evaluate. This is accomplished by processing paths generated using one protocol -the sampling protocol, Λ s -using another -the analysis protocol, Λ.
While we have infinite freedom in selecting an analysis protocol, not all choices will improve the convergence of free energy estimates. One reasonable strategy for choosing Λ is to minimize the lag, the difference between the nonequilibrium and instantaneous equilibrium densities; Vaikuntanathan and Jarzynski found that under certain dynamics, dissipation is eliminated if there is no lag, leading to a zero-variance estimator of F Λ . 23 To reduce the lag, they modified their equation of motion with an additional flow-field term that "escorts" the system along a near-equilibrium path. Essentially, this strategy modifies the nonequilibrium density. In this paper, we take the opposite approach: using the analysis protocol to choose an instantaneous equilibrium density that closely matches the sampled nonequilibrium density.
As an illustrative case, consider a Brownian particle in a harmonic oscillator, or spring, which moves at a constant velocity (Fig. 1). If the system starts in thermal equilibrium, its density is a Gaussian about the initial spring position. When the spring starts moving, the density remains a Gaussian with the same variance, but its mean position, x T (t), lags behind the spring position. 31,32 For this particular system, an analysis protocol based on x T (t) will have no lag. We shall further explore this system in Section III.
One complication with using an analysis protocol that minimizes the lag is that its end state is usually not the same as in the sampling protocol. Thus, the free energy difference being estimated differs. To estimate the same F Λ with a minimal lag analysis protocol, it may be necessary to extend or otherwise modify the sampling. To distinguish the two situations, we shall refer to the former as protocol postprocessing and the latter as nonequilibrium density-dependent sampling (NEDDS). Both fall under the aegis of density-dependent analysis.
The structure of this paper is as follows: in Section II, the importance sampling form of Jarzynski's equality is detailed; in Section III, density-dependent analysis is demonstrated on two cases in which the propagator is analytically known; in Section IV, a general method for finding minimal lag analysis protocols is described, ap-plied to an adaptive algorithm for NEDDS, and tested on the model system; and lastly, implications of this method and possible future directions are discussed.

II. FREE ENERGY FORMALISM
Consider a system whose Hamiltonian, H = H(x; λ), depends on x, its position in phase space (or configuration space), and a control parameter, λ. Initially, the system is prepared in thermal equilibrium at λ(0). The parameter λ is perturbed according to a protocol Λ = λ(t) until it reaches a final state at λ(τ ). Jarznyski's equality, 1,2 relates the free energy difference between the initial and final states of the protocol, F Λ , to an average over all possible paths, X = x(t), resulting from this nonequilibrium procedure. Specifically, this expectation (denoted by the angled brackets ... Λ ), is a path integral over infinitesimal elements dX with the protocol-dependent density where N s is the sample size. In a standard Jarzynski estimate, r =  27 By variation of the asymptotic variance with respect to the sampling bias, Oberhofer and Dellago found that optimal work-weighted sampling is given by 28 Unfortunately, this optimal choice is impractical because it includes the sought quantity F Λ .
In the present method, which applies Eq. (3) in a novel manner, ρ s [X] = ρ Λs [X] depends on the sampling protocol and r differs from unity when the analysis protocol Λ varies from Λ s . Notably, the relevant work is W [X|Λ], not W [X|Λ s ], meaning that different choices of Λ will result in various work distributions and convergence properties. This new way of applying importance sampling leads to different, albeit analogous, asymptotic variance expressions. 33 The present approach is more general than previous applications of Eq. (3), which require transition path sampling, because it does not require biased sampling and paths can be generated by ordinary dynamical equations. Indeed, under certain assumptions, such as those suggested by Nummela and Andricioaei, 34 it should be possible to apply the present method to laboratory experiments.
We note, as a caveat, that the importance sampling form of Jarzynski's equality will only be useful for stochastic dynamics where r can be computed. Under deterministic dynamics, r is a delta function and having different sampling and analysis protocols will not improve free energy estimates.

III. CASES WITH AN ANALYTICAL PROPAGATOR
As mentioned earlier, we would like to choose an analysis protocol that minimizes the lag, such that the instantaneous equilibrium density corresponds with the sampled nonequilibrium density. This is particularly tractable when the propagator is exactly known. Here, we demonstrate Eq. (3) on two such cases: a Brownian particle in a harmonic oscillator (i) moving at a constant velocity or (ii) with a time-dependent natural frequency. With both, the potential energy has the general form U (x) = k(x −x) 2 /2 and the nonequilibrium density is where x T (t) and k T (t) are the most typical paths and spring coefficients, respectively. As these propagators can be obtained by close analogy to the path integral derivation of work-weighted propagators, 32 their derivations are not detailed here. In case (i), k is constant and λ moves the spring position according tox = Λ s = vt, such that ∆F = 0, k T (t) = k, and In case (ii),x is zero and λ controls the spring coefficient, k = Λ s , such that ∆F = 1 2 ln[k(0)/k(τ )]. In the corresponding nonequilibrium density, x T (t) = 0, and Based on these expressions, it is evident that, for case (i), the minimal lag analysis protocol is Λ ml = x T (t) from Eq. (5), and for case (ii), it is Λ ml = k T (t) from Eq. (6). In these special cases, the nonequilibrium density is exactly the equilibrium density corresponding to Λ ml and there is no lag.
To test whether density-dependent analysis leads to improved free energy estimates, one-dimensional Brownian dynamics simulations were run with the equation of motion, where x j is the position at time j∆t, ∆t is the time step, and R j is a standard normal random variable. The primes denote spatial derivatives such that The work was evaluted with the discrete formula, . This action is valid in the continuum limit, as J → ∞ and ∆t → 0. To approach this limit, we chose D = 1 and a time step of ∆t = 0.001.
The simulations were performed over 10 m steps (truncated to be an integer), where m refers to 7 evenly spaced numbers between 1.5 and 3. In case (i), k was set to 25 and Λ s was chosen to start from λ 0 = 0 and linearly progress to the target state at λ f = 1. With case (ii), Λ s is a linear interpolation between 1 and 100. Afterwards, the trajectories both analyzed with the standard Jarzynski estimate and subjected to protocol postprocessing with Λ ml .
For comparison, NEDDS was implemented by switching λ at a faster rate such that the final state went beyond λ f and the final nonequilibrium density, according to the propagators, corresponded to the target state. This is illustrated in Fig. 1d, where moving the harmonic oscillator at a faster rate than in Fig. 1c allows for the final density to correspond to the equilibrium state with λ = 1. These trajectories, which took the same amount of simulation time for the same number of steps, were then reanalyzed with the appropriate Λ ml .
In case (i), we find that protocol postprocessing with Λ ml leads to a desirable result: most work values are reduced such that a larger fraction of them are less than the free energy difference (Fig. 2). Of these negative dissipation trajectories, most have a probability ratio less than one. Conversely, several positive dissipation trajectories have a probability ratio greater than one. For this set of trajectories, the modified work distribution leads to a more accurate free energy estimate.
Over a large number of repetitions and range of switching speeds, we find that free energy estimates based on Λ = Λ ml are vastly improved over the standard procedure, Λ = Λ s , having significantly less variance and systematic bias (Fig. 3). The standard estimator only approaches the accuracy and precision of protocol processing for slow switches. Clearly, these trajectories are much better at estimating the end state free energy differences for Λ ml than for Λ s . The estimates of F Λ from NEDDS also require considerably less sampling than the standard procedure, although the effect is somewhat less dramatic.
Similarly, in case (ii), density-dependent methods also show improvement over the standard Jarzynski estimate. For the time-dependent natural frequency, the systematic bias of the standard estimate is relatively small but nonetheless evident at all sampled switching rates (Fig. 4). Estimates from both density-dependent methods have reduced bias and variance, and are found to be of similar quality to each other.

IV. GENERAL CASE
In most practical situations, unfortunately, the propagator is not known ahead of time. Thus, prior to simulations, it is unclear how long paths need to be gener- ated before the nonequilibrium density matches a density characteristic of the target state. While paths are being generated, however, it is possible to estimate the difference between the sampled density and arbitrary equilibrium states. States which minimize this difference can be be collected in an analysis protocol with minimal lag. One measure of the distance between two probability distributions is the Kullback-Leibler divergence, or the relative entropy. The relative entropy between the nonequilibrium density and an arbitrary equilibrium state T is, When the integral is separated into two at the logarithm, one part is a constant with respect to T . The divergence is minimized by finding a state where the other, − dx p neq (x, t) ln p T (x), is the least. Using sampled discrete paths, this integral can be estimated by Ns n=1 ln p T (x jn ), where x jn is the position at step j of path n. For a state T , the equilibrium density is is the test state Hamiltonian and F T is its free energy. Thus, the relative entropy is minimized by the smallest value of, among different states T . Generally, the free energy, F T , is unknown, but for states which occur along the switching protocol, F T − F 0 (where F 0 is the free energy at λ 0 ) can be estimated using the standard form of Jarzynski's equality. These states constitute our search space for minimizing the lag. Suppose we are interested in the free energy difference between the states defined by λ 0 and λ f . We can use D T to estimate Λ ml on the fly and determine when to stop sampling via the following adaptive algorithm: 1. Start with j = 0 and the work W 0 = 0. For each of N s paths, obtain x 0 by drawing samples from the equilibrium ensemble at λ 0 .
2. Propagate each path, calculating x j+1 using a dynamical equation such as Eq. 7. To obtain W j+1 , calculate the work done on the system during the time step and add it to W j . The next step in the sampling protocol, λ j+1 , is found by adding a predetermined value, µ, to λ j . The sign of µ must be the same as λ f − λ 0 . Increment j by one.
3. Using W j values in the standard form of Jarzynski's equality, estimate F j −F 0 , the free energy difference between the states with λ j and λ 0 .
4. For each state T corresponding to {λ 0 , λ 1 , ... λ j }, use H T (x) and the free energy difference estimated in the previous algorithm step to calculate D T −F 0 . The λ which minimizes D T − F 0 is λ ml . Add λ ml to the minimal lag protocol Λ ml .
5. If λ ml hasn't crossed λ f , repeat from algorithm step 2. Otherwise, set the final value in Λ ml to λ f . 6. Estimate the free energy difference using Eq. (3) with Λ = Λ ml .
This algorithm was tested on Sun's system, 24 where the potential energy is U (x) = x 4 − 16λx 2 . Using Eq.
(3), the free energy difference was estimated between the initial state with λ 0 = 0, where the potential is a single well, and the target state λ f = 1, where it is a double well, such that ∆F = −62.9407. 27 Brownian dynamics simulations were performed with the same diffusion coefficient, time step, and equation of motion as in Section III. The increment of λ at each time step was µ = v∆t, where v = 10 m and m refers to 9 evenly spaced values between 0 and 2. For comparison, the standard Jarzynski  estimate was applied to simulations where λ is switched between 0 and 1 at a slower velocity, taking the same total time as in the corresponding NEDDS simulations.
In a representative set of simulations, the density most noticeably lags behind the sampling state at the beginning (Fig. 5). Around the state defined by λ = 0.9, the lag quickly diminishes. However, the minima of D T does not reach the target state until the sampling λ is beyond 1.
Based on many repetitions of this procedure at different pulling speeds, we find that our NEDDS algorithm converges much more quickly than the standard Jarzynski estimate (Fig. 6). The systematic bias is largely eliminated with simulations that are switched nearly an order of magnitude faster. At the fastest switching rates, NEDDS remains biased but still outperforms the standard Jarzynski estimate. With these fast switchings, it is possible that the nonequilibrium density does not correspond well to any traversed equilibrium state.

V. DISCUSSION AND CONCLUSION
With the goal of minimizing the lag via the choice of analysis protocol, we have developed density-dependent methods to analyze nonequilbrium paths, to estimate which states may constitute a protocol that minimizes the lag, and to adaptively sample paths until the desired density is achieved. Our promising results validate the strategy and provide further evidence for the link between lag and heat dissipation. They also hint that the accurate estimation of free energy differences may require adequate sampling in the important regions of both end states.
Analysis protocols provide another degree of freedom for lag reduction, and can be used in conjunction with other methods, such as sampling protocol optimization or biased path sampling. Furthermore, their use should extend beyond Jarzynski's equality; they can poten-tially be applied in bidirectional nonequilibrium work expressions 33 or any relationship between a nonequilibrium process and a state function, such as Hummer and Szabo's expression for the potential of mean force. 35 Quite possibly, our results are just the tip of an iceberg and this paper will open up new research directions for sampling and analyzing nonequilibrium trajectories.

VI. ACKNOWLEDGMENTS
The author thanks Artur Adib, Christopher Jarzynski, Attila Szabo, and Suriyanarayanan Vaikuntanathan for pertinent discussions, and Gerhard Hummer for suggesting that he considers the lag. He also thanks Artur Adib for supporting a postdoctoral fellowship. This research was supported by the Intramural Research Program of the NIH, NIDDK.
In these appendices, we derive asymptotic variance and bias expressions for free energies estimated using protocol postprocessing. Our derivations are similar to that of Oberhofer et. al. for biased sampling of nonequilibrium trajectories. 27,28 Although these expressions are not directly relevant to the main text, they support the line of inquiry explored in the paper and may be useful to those who wish to expand upon the results. Both unidirectional and bidirectional expressions are considered.

APPENDIX A: UNIDIRECTIONAL
Expressed in importance sampling form, Jarzynski's equality is, This is Eq. (3) in the main text, reproduced here for convenience. Towards deriving the asymptotic expressions, we first define, The sample mean estimators for the expectations of A and B are, These form an estimator for the free energy,F However, the sample mean estimators deviate from their true expectations, Here, subscripts on the angle brackets are omitted for notational simplicity. Deviations in these expectations lead to variance and bias inF Λ . The magnitude of the error can be estimated by a Taylor series expansion about F Λ , which becomes an increasingly reasonable approximation in the large sampling, or asymptotic, limit. To the first order, this expansion is, The partial derivatives are evaluated at their mean values. The variance is defined as Using the first-order Taylor series expansion, this is, The variance of a sample mean is the variance of the variable over the number of samples. Since A and B are functions of the same data points, the covariance of their sample means is, similarly, σ 2 [Ā,B] = 1 N σ 2 [A, B]. Thus, the asymptotic variance is, The bias is defined as B N ≡ F Λ − F Λ . If we approximate this error with a first-order Taylor series expansion, it is always zero. In order to obtain a nonzero bias expression, we use a second-order Taylor series expansion about F Λ , Using this approximation, the bias is found to be, Notably, when the dissipated work is zero, W [X|Λ] − F Λ = 0, expressions for both the variance and bias are likewise zero.

APPENDIX B: BIDIRECTIONAL
Here, we consider the possibility of reanalyzing bidirectional data, collected using both a protocol Λ and its time reversalΛ. The results derived in this section suggest that for bidirectional data, the optimal analysis protocol is actually the sampling protocol. Thus, protocol postprocessing is less promising when applied to bidirectional data than to unidirectional data.
For notational consistency, we start our discussion with the Crooks Fluctuation Theorem, 3,4 As in the main text, ρ Λ [X] is the probability of observing trajectory X, given the protocol Λ. Analogously, ρΛ[X] is the probability of observing the time reversal, or conjugate twin, of X, using the reverse protocolΛ.
This theorem can be used to derive a relationship between forward and reverse path-ensemble averages, 5 In the above, F Λ [X] is an arbitrary functional. As the last path-ensemble average is over trajectoriesX, trajectories sampled fromΛ must be reversed prior to being evaluated with the functional. Next, we rearrange the path-ensemble average theorem into an expression for the free energy. 5 As in the unidirectional case, these path-ensemble averages can be written as reweighed samples from other sampling densities.
The probability ratior is defined similarly to r,r The ratio of r s / r s can be shown to be unity by converting r s into a forward path-ensemble average using Eqs. (B4) and (B9), and applying the importance sampling form of Jarzynski's equality, Eq. (A1). The asymptotic variance of Eq. (B7) can be calculated by a similar procedure to the unidirectional case. We start by defining, Replacing A and B with C and D, we follow the same logic as in the unidirectional case from Eqs. (A6) to (A14). Next, we note that C and D are independent samples drawn from different ensembles and their correlation is zero. Thus, the variance estimate is, We would like to combine the two terms including C and D in a single path-ensemble average. In order to do that we need to convert the ensemble averages containing D to the forward direction. For D , this is,   We would like to express this equation in a form in which it is clear how to include data sampled from forward and reverse path-ensembles. The path-ensemble average in the numerator is, As with unidirectional data, a natural question to ask is how to find the optimal analysis protocol when processing bidirectional data. It no longer makes sense to reduce the lag; a protocol which reduces the lag for forward trajectories may increase it for trajectories from the reverse protocol. Furthermore, when Eq. (B34) is variationally optimized with respect to r, the optimal r is found to be constant. Since this only occurs when Λ = Λ s , this result suggests that, compared with unidirectional data, protocol processing less likely to be useful for treating bidirectional data.