Density-Dependent Analysis of Nonequilibrium Paths Improves Free Energy Estimates II. A Feynman-Kac Formalism
The nonequilibrium fluctuation theorems have paved the way for estimating equilibrium thermodynamic properties, such as free energy differences, using trajectories from driven nonequilibrium processes. While many statistical estimators may be derived from these identities, some are more efficient than others. It has recently been suggested that trajectories sampled using a particular time-dependent protocol for perturbing the Hamiltonian may be analyzed with another one. Choosing an analysis protocol based on the nonequilibrium density was empirically demonstrated to reduce the variance and bias of free energy estimates. Here, we present an alternate mathematical formalism for protocol postprocessing based on the Feynmac-Kac theorem. The estimator that results from this formalism is demonstrated on a few low-dimensional model systems. It is found to have reduced bias compared to both the standard form of Jarzynski’s equality and the previous protocol postprocessing formalism.
A key goal in computational thermodynamics is the estimation of free energy differences between equilibrium states. Challenges in efficiently obtaining accurate values, however, continue to motivate the development of novel methods.Chipot2007 () The discovery of new theorems in nonequilibrium statistical mechanics Jarzynski1997a (); Jarzynski1997b (); Crooks1998 (); Crooks1999 (); Crooks2000 () have opened up a promising direction: free energy calculations based on simulations of driven nonequilibrium processes. Frenkel2002 (); Chipot2007 () The most straightforward implementation of this approach involves performing multiple repetitions of a process in which a system is driven out of equilibrium by switching an external parameter according to a protocol , where . If the free energy difference of interest is between thermodynamic states defined by setting to and , then the protocol is defined so that and are the end states, and . The free energy difference between the initial thermodynamic state and the equilibrium state at any time , , may then computed by applying, Jarzynski1997a (); Jarzynski1997b ()
where is an average over all possible trajectories (realizations of the process), and denotes the work done on the system up to time during a particular trajectory, . For a finite sample of trajectories for n = 1, 2, …, N, the sample mean provides an estimator for this expectation.
Unfortunately, this estimator often suffers from poor convergence. The expected number of realizations needed to obtain a reliable estimate of grows rapidly with the dissipation, , that invariably accompanies driven nonequilibrium processes. Gore2003 (); Jarzynski2006 (); Kofke2006 () In turn, dissipation reflects the lag that develops between the state of the system and the equilibrium state corresponding to the instantaneous value of the external parameter.Vaikuntanathan2008 (); Vaikuntanathan2009 () (See Fig 1). This lag is ultimately responsible for the poor convergence of Eq. 1.
The connection between lag and the convergence of the free energy estimator can be better understood by considering two limiting cases. First, consider the case of an infinitely slow reversible process. As the system remains in equilibrium throughout the process, there is no lag. In this case, convergence only requires a single sample because the work performed along any isothermal quasi-static trajectory, , is equal to the free energy difference, . LandauLifshitz () The opposite limit is that of an infinitely fast process, in which Eq. 1 reduces to the more familiar free energy perturbation identity. Free energy estimates based on this identity converge quickly only if there is significant overlap in the important phase space regions of the end states,Chipot2007 (); Kofke2006 () which in turn reflects the lag. Likewise, in the intermediate situation of a finite-time process, the convergence of free energy estimates depends on overlap between the sampled phase space and the important phase space of an equilibrium state of interest.
In order to reduce lag and improve convergence, strategies such as importance sampling of trajectories Athenes2002 (); Sun2003 (); Atilgan2004 (); Ytreberg2004 (); Athenes2004 (); Oberhofer2005 (); Oberhofer2008 () and “escorted” free energy simulations Vaikuntanathan2008 () have been introduced. (For a brief overview, see Ref. Minh2009 ().) In this paper, we consider an alternate strategy, protocol postprocessing. This strategy involves introducing a function with , which we will refer to as the analysis protocol. The central result of this paper (Eq. 13) is an expression for the free energy difference using trajectories generated in the original process (in which the work parameter is switched according to the protocol ). While this result is valid for any choice of and reduces to Eq. 1 for , we will argue that Eq. 13 provides efficient estimates of the free energy difference whenever the equilibrium densities corresponding to the analysis protocol have a high degree of overlap with density of the system (See Fig 1).
Protocol postprocessing was previous introduced Minh2009 () in the context of importance sampling in path-space. Athenes2002 (); Sun2003 (); Atilgan2004 (); Ytreberg2004 (); Athenes2004 (); Oberhofer2005 (); Oberhofer2008 () In the present work, we utilize an alternate mathematical formalism, the Feynman-Kac theorem. Kac1949 () This formalism is similar to that used in the escorted free energy simulation method, Vaikuntanathan2008 () and indeed, the two methodologies may be used in conjunction with one another. The new formalism has at least two advantages over the previous method: first, in certain special cases, it is analytically a zero-variance estimator. Secondly, for a few simple model systems, we find that the bias and variance of free energy estimates are substantially reduced.
Ii Feynman-Kac Formalism
The derivation of Jarzynski’s equality Jarzynski1997a (); Jarzynski1997b () from the Feynman-Kac theorem has been well-documented. Hummer2001a (); Hummer2005 (); Ge2008 () In this section, we recapitulate Hummer and Szabo’s derivation Hummer2001a () and extend it to protocol postprocessing.
ii.1 Jarzynski’s equality
Consider a classical system with a Hamiltonian, , that depends on its position in -dimensional phase (or configuration) space, , and a parameter vector, . The system evolves according to dynamics which, if the temperature and are held constant, preserve the canonical equilibrium distribution , where is a partition function. These conditions are satisfied by several dynamics such as Hamilton’s equations, Langevin dynamics, and the Andersen and Nosé-Hoover thermostats.
We are interested in driven nonequilibrium processes in which the system is first prepared in equilibrium with and temperature , after which the external parameters are switched according to the protocol . Each realization of this process can be described by the trajectory, . The phase space density of an ensemble of such trajectories evolves according to a Liouville-type equation,
Hummer and Szabo recognized that the Feynman-Kac theorem provides a solution to the “sink equation”,
We remind the reader that the angled brackets denote a path-ensemble average, or expectation, over all possible realizations of the described driven nonequilibrium process.
Another solution to Eq. 3 is given by an improperly normalized Boltzmann distribution, . Equating this solution to that from the Feynman-Kac theorem immediately gives,
in which the work is defined as,
By integrating both sides over , we obtain
ii.2 Protocol Postprocessing
In the protocol postprocessing strategy, trajectories are first generated according to the sampling protocol . Next, a potentially distinct analysis protocol , with , is introduced. This analysis protocol is not used to generate any new trajectories. Rather, the previously generated trajectories are used as samples for estimating the free energy difference . The standard form of Jarzynski’s equality can be seen as a special case where the sampling and analysis protocols are identical. While the formalism described below is valid for any , it will not always be advantageous. In Section IV, however, we will describe how to choose a that leads to an efficient free energy estimate.
We begin the derivation by formally separating the evolution operator into two terms,
where the auxiliary operator represents the difference between the evolution operators given the sampling and analysis protocols.
Now consider a sink equation analogous to Eq. 3,
where the function not only includes a time-derivative of the Hamiltonian, but also a term with the operator ,
Here, the operator only acts on the term in the numerator. One solution to Eq. 9 is .
By equating this solution to the path integral solution obtained from the Feynman-Kac theorem, we obtain an equation analogous to Eq. 5:
where the work has the modified form,
Integrating over , we obtain a protocol postprocessing form of Jarzynski’s equality,
Again, the angled brackets denote a path-ensemble average, or expectation, over all possible realizations of the driven nonequilibrium process with the protocol ; the protocol has nothing to do with sampling.
To be more concrete, let us consider a system moving with overdamped Langevin (Brownian) dynamics in a one-dimensional potential . The density evolves according to the Smoluchowski equation,
where is the diffusion coefficient and the prime symbol represents a derivative with respect to .
Given an analysis protocol , the auxiliary operator for this example system is defined as,
Substituting this expression into Eq. 12, we obtain a modified form of the work,
Using this expression for in Eq. 13, we can now estimate the free energy difference from trajectories generated in the process in which external parameter is switched according to the protocol .
For dimensions indexed by , the Smoluchowski equation is,
where is the diffusion coefficient in dimension . Following steps analogous to those above, we obtain the modified work,
where all are implicitly functions of the position at time , and is also a function of .
Iii Importance Sampling Formalism
Section II.2 is not the first description of protocol postprocessing; it was preceded by a formalism based on importance sampling. Minh2009 () In this section, we describe the previous formalism in the current notation and compare it with the present results.
Explicitly in terms of path integrals, we may rewrite Eq. 1 as,
where denotes the work performed on the system as it evolves along a particular trajectory in which the external parameter is changed according to the protocol , is the probability density associated with the trajectory , and is a metric over paths.
Now suppose that the external parameter is changed according to the protocol for which the associated probability density of a trajectory is . The same free energy difference may be computed by estimating different path integrals, Ytreberg2004 (); Minh2009 ()
where is the ratio of densities. If the two protocols sampling are equivalent, then .
This expression differs from Eq. 13 in that it includes two expectations, the definitions of work are different, and it requires a ratio of probabilities, . The ratio is different from a “modification” to the work term. For example, in overdamped Langevin dynamics, this ratio is, MinhAdib2009 (); Minh2009 ()
Now suppose that we break down in Eq. 17, into one term with and a “modification” term. If we multiply this modification term by and take the exponent, we obtain a term which is used similarly to ,
but is quite distinct.
For multiple dimensions of overdamped Langevin dynamics, the ratio is,
where again, are implicitly functions of and time . As in the 1D case, this expression is not equivalent to the Feynman-Kac estimator.
In later sections, we will describe several advantages of the new formalism.
Iv Dissipation and Lag
As protocol postprocessing is merely another mathematical formalism for computing free energies, there is no a priori reason to expect that it will perform any better or worse than the usual nonequilibrium work estimator, Eq. 1. For clever choices of the analysis protocol, however, we can show that Eq. 13 leads to a highly efficient estimator for .
iv.1 Exactly solved models
Suppose that we construct a “perfect” analysis protocol whose instantaneous equilibrium density is equivalent to the nonequilibrium density, so that , where denotes the equilibrium distribution corresponding to and . When a perfect analysis protocol is used, then
for every trajectory! This may be seen by first substituting in the evolution equation,
where we have used Eq. 8. Since for this , we obtain,
By substituting this into the modified work, Eq. 12 and integrating, we obtain Eq. 25. As this equation is valid for every trajectory, Eq. 13 is a zero variance estimator of . As a demonstration of this principle, consider two exactly solved Minh2009 () models: a Brownian particle in a harmonic oscillator that either (i) has its center moving at a constant velocity, or (ii) has a time-dependent natural frequency. In both cases, the potential has the general time-dependent form where the vector denotes the set of external parameters. The Smoluchowski equation describing the evolution of the phase space density can be solved to give MinhAdib2009 (); Minh2009 ()
where , , and denotes an average over the distribution . In case (i), the spring coefficient is held fixed at while is switched according to (). In this case, the free energy difference is always zero and is a constant, . The most typical path is,
In case (ii), is held fixed at and the spring coefficient is switched according to (). In this case, , and
In either case, we may choose the analysis protocol such that . With this choice, the Boltzmann distribution corresponding to the analysis protocol is equal to . Hence, the modified work calculated from Eq. 17 is always equal to the free energy difference . In contrast, the importance sampling form of protocol postprocessing yields different work values for each trajectory.
iv.2 Dissipation Bounds Lag
In general, it is not feasible to find a perfect analysis protocol. Indeed, in most cases, the nonequilibrium densities will not belong to the family of equilibrium distributions indexed by , . However, Eq. 25 suggests that efficient estimators of free energy energies can be obtained if we can find an analysis protocol such that closely resembles the nonequilibrium density . In the following paragraphs, we will make this argument more rigorous.
The convergence of the protocol postprocessing form of Jarzynski’s equality will depend on a criterion analogous to that in the original form: obtaining trajectories in which the modified work, , is less than the free energy difference, . Jarzynski2006 (); Kofke2006 () Chances of obtaining such trajectories are improved when the average dissipation, , is small. Jarzynski2006 (); Kofke2006 () This dissipation can be related to an information theoretic measure of overlap between the distributions describing the state of the system and the equilibrium state corresponding to the , . To obtain this relation, we note that the properties of the delta function enable the path-ensemble average in Eq. 11 to be written as,
where the double subscript indicates a path-ensemble average for trajectories driven with the protocol and which pass through at time . Since the nonequilibrium density at time is , we may rearrange Eq. 11 to obtain,
As in Ref. Vaikuntanathan2009 (), we then take the logarithm of both sides of the equation, invoke Jensen’s inequality, multiply both sides by , and integrate over . Our final result is,
where is the Kullback-Leibler divergence, or relative entropy, between the nonequilibrium density and the equilibrium density corresponding to the analysis protocol. The relative entropy is zero when two distributions are identical and grows larger when they diverge Cover1991 (). Eq. 33 suggests, but does not prove (the inequality goes the wrong way), that a reasonable strategy for reducing dissipation and improving the convergence of the free energy estimator is to choose an analysis protocol in which the “analysis” density closely resembles the evolving state of the system.
V General Case
Based on the results in Section IV, we speculate that a reasonable strategy for minimizing dissipation and improving the efficiency of the free energy estimator is to choose an analysis protocol so that the Kullback-Leibler divergence is small for all . Obtaining such a protocol will usually entail a search over the space of to find an equilibrium distribution which is similar to . While the nonequilibrium distribution is not analytically tractable for most systems, it is possible to use sampled trajectories to compare the relative entropy between and for different values of . Specifically, given a set of trajectories and several candidate values of , the relative entropy is minimized by the parameter vector that minimizes , which may be estimated by the sample average, Minh2009 ()
where denotes the state of system in phase space at time as it evolves along the trajectory . (Note that in , the integral does not depend on .) A reasonable choice for the search space of is the range of the sampling protocol . This choice has the advantage that may be estimated via Jarzynski’s equality; for distributions that are not accessed during the sampling protocol, it may be more difficult to estimate corresponding free energies.
As noted in Section II, the flexibility in choosing means that the free energy may be different from . Indeed, unless there is no lag, an analysis protocol which minimizes the lag will always have different states than the sampling protocol. Since we are typically interested in free energies between the end states of the sampling protocol ( and ), this discrepancy was addressed by introducing an adaptive algorithm, nonequilibrium density-dependent sampling (NEDDS). Minh2009 () NEDDS is equally applicable to the current formalism.
In brief, NEDDS entails running all desired simulations of the nonequilibrium process simultaneously. The sampling protocol initially involves an interpolation between the desired end states and . After reaching state , the protocol extrapolates past it until an adaptively determined stopping time. (While such an extrapolation may not always be physically meaningful, it is nearly always computationally feasible.) Without loss of generality, let us assume that . The stopping time is decided by performing the following calculations while the simulations are in progress:
The free energy difference, , between the initial and instantaneous state at the current time step, , is estimated using Eq. 1.
is evaluated with values from the current state and all preceding states using Eq. 34.
If the choice of that minimizes , , is between and , , then it is appended to the analysis protocol, . Otherwise, if it is at or beyond , , then the final value of the analysis protocol is set to , .
Lastly, is incremented and is evaluated by protocol postprocessing.
This procedure ensures that protocol postprocessing estimators can compute the free energy difference between the states and .
Vi Model Systems
We now demonstrate NEDDS with protocol postprocessing (both importance sampling and Feynman-Kac) formalisms and compare its efficiency to standard sample mean estimates from Jarzynski’s equality, Eq. 1, on a few toy model systems. First, consider an overdamped Brownian particle evolving on an one-dimensional surface, , as studied by Sun. Sun2003 () In this system, the free energy difference between the states and at was analytically found to be . Oberhofer2005 ()
As described, Minh2009 () simulations of nonequilibrium driven processes were performed in which was switched between 0 and 1 according to the equation of motion,
where is the position at time , is the diffusion coefficient, is the time step, and is a standard normal random variable. was incremented at each time step by . NEDDS was used to obtain the analysis protocol concluding at , and the free energy difference was computed using either Eq. 13 or Eq. 21. For comparison, the standard Jarzynski estimate was applied to two types of simulations taking the same amount of simulation time as the analysis protocol obtained from NEDDS: either (i) was switched between 0 and 1 at a slower velocity, or (ii) the NEDDS analysis protocol was used as a new sampling protocol.
While the importance sampling formalism was found to be an improvement over the standard form of Jarzynski’s equality, Minh2009 () we find that the estimator based on Eq. 13 is even better (Fig. 2). Even for the fastest switching rates, where dissipation is expected to be high, the systematic bias is largely eliminated. No benefit was found from using the analysis protocol from NEDDS as a new sampling protocol; in fact, the bias was worse than with the constant velocity protocol.
We also performed similar tests on another one-dimensional surface, , first described by Hummer. Hummer2007 () Hummer’s surface, a double well potential that includes a harmonic bias, mimics the setup of a single-molecule pulling experiment, and hence has been used to demonstrate estimators of free energies MinhAdib2008 (); MinhChodera2009 () and other quantities MinhChodera2010 () in the context of these experiments. The simulations were performed using the same equation of motion, diffusion coefficient, and time step as described above for Sun’s system. was switched between -1.5 and 1.5.
The performance trends with Hummer’s system are similar to those with Sun’s (Fig. 3). Results from the standard form of Jarzynski’s equality are more biased than with NEDDS and the importance sampling formalism, which in turn is more biased than the Feynman-Kac formalism. In contrast to Sun’s system, however, the estimates from Eq. (13) are noticeably biased at the fastest switching rates. Another distinction between the trends from the two systems is that results from using a constant velocity protocol and using the analysis protocol as a new sampling protocol are rather similar.
As a final demonstration, we consider a two-dimensional surface,
in which dictates the progress of a harmonic bias along a curve (Fig. 4). Simulations were performed as with the 1D potentials, using Eq. 35 along each dimension, as well as the same diffusion coefficient and time step. was switched between 0 and 1.
The performance trends in this system are the same as in Hummer’s system (Fig. 5).
Vii Discussion and Conclusion
We have presented a method for analyzing nonequilibrium trajectories which borrows from a similar philosophy as previous work Minh2009 () but is based on a distinct mathematical formalism. The new formalism has the advantages that it analytically is a zero-variance estimator if a “perfect” analysis protocol is obtained, and it improves the convergence of free energy estimates in all our tested model systems. Further tests on more complex multidimensional systems are a potential future research direction.
We expect that protocol postprocessing will be most useful when (i) there is little phase space overlap between the end states of interest (otherwise free energy differences can be computed without nonequilibrium work identities), (ii) estimates of from the nonequilibrium work relation suffer from poor convergence for a given nonequilbrium process in which the system is driven between the end states of interest and (iii) it is reasonable to speculate that the nonequilibrium driven process has a nonequilibrium density that always resembles an equilibrium density parameterized by a vector along the protocol. Exact convergence properties, of course, will depend on the system.
The observed convergence benefits hint that many improved sampling and analysis algorithms based on nonequilibrium driven processes still remain to be discovered.
We thank Andy Ballard, John Chodera, and Christopher Jarzynski for helpful comments on the manuscript. D. Minh is funded by a Director’s Postdoctoral Fellowship at Argonne and S. Vaikuntanathan acknowledges support from the National Science Foundation (USA) under CHE-0841557 and the University of Maryland, College Park.
- C. Chipot and A. Pohorille, Free Energy Calculations (Springer, Berlin, 2007).
- C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997).
- C. Jarzynski, Phys. Rev. E 56, 5018 (1997).
- G. E. Crooks, J. Stat. Phys. 90, 1481 (1998).
- G. E. Crooks, Phys. Rev. E 60, 2721 (Sep 1999).
- G. E. Crooks, Phys. Rev. E 61, 2361 (2000).
- D. Frenkel and B. Smit, Understanding Molecular Simulation (Elsevier, 2002) pp. 183–189.
- J. Gore, F. Ritort, and C. Bustamante, Proc. Natl. Acad. Sci. U.S.A. 100, 12564 (2003).
- C. Jarzynski, Phys. Rev. E 73, 046105 (2006).
- D. A. Kofke, Mol. Phys. 104, 3701 (2006).
- S. Vaikuntanathan and C. Jarzynski, Phys. Rev. Lett. 100, 190601 (2008).
- S. Vaikuntanathan and C. Jarzynski, Europhys. Lett. 87, 60005 (2009).
- L. Landau and E. Lifshitz, Statistical Physics, 3rd ed. (Pergamon Press, Oxford, 1990).
- M. Athenes, Phys. Rev. E 66, 046705 (2002).
- S. Sun, J. Chem. Phys. 118, 5769 (2003).
- E. Atilgan and S. X. Sun, J. Chem. Phys. 121, 10392 (2004).
- F. M. Ytreberg and D. M. Zuckerman, J. Chem. Phys. 120, 10876 (Jan 2004).
- M. Athenes, Eur. Phys. J. B 38, 651 (2004).
- H. Oberhofer, C. Dellago, and P. Geissler, J. Phys. Chem. B 109, 6902 (Jan 2005).
- H. Oberhofer and C. Dellago, Comput. Phys. Commun. 179, 41 (2008).
- D. D. L. Minh, J. Chem. Phys. 130, 204102 (2009).
- M. Kac, Trans. Am. Math Soc. 65, 1Ð13 (1949).
- G. Hummer and A. Szabo, Proc. Natl. Acad. Sci. U.S.A. 98, 3658 (2001).
- G. Hummer and A. Szabo, Acc. Chem. Res. 38, 504 (2005).
- H. Ge and D. Q. Jiang, J. Stat. Phys. 131, 675 (2008).
- D. D. L. Minh and A. B. Adib, Phys. Rev. E 79, 021122 (2009).
- T. M. Cover and J. A. Thomas, Elements of Information Theory (John Wiley and Sons, Inc., 1991).
- G. Hummer, in Free Energy Calculations, Vol. 86, edited by C. Chipot and A. Pohorille (Springer, Berlin, 2007).
- D. D. L. Minh and A. B. Adib, Phys. Rev. Lett. 100, 180602 (2008).
- D. D. L. Minh and J. D. Chodera, J. Chem. Phys. 131, 134110 (2009).
- D. D. L. Minh and J. D. Chodera, J. Chem. Phys. in press (2010).