Reversibility of free energy simulations: Slow growth may have a unique advantage. (With a note on use of Ewald summation)

Abstract

We review the slow-growth method for computing free energy changes for processes in conformation space or in "chemical" space, in which a system parameter, x is changed at each integration time step, and the free energy, ΔA is approximated by accumulating the work performed at each step. The method is simple to implement and use, convergence can be monitored by performing longer simulations and by performing the simulations changing x in both directions, and statistical error can be evaluated by performing multiple independent simulations. Because slow growth simulates a continuous process, it closely approximates the ideal isothermal quasi-static process used in defining the free energy in thermodynamics, and thus a small hysteresis in slow-growth results practically guarantees that the process is reversible, which is of course a prerequisite for the results to represent a free energy change. Whenever hysteresis is not negligible (which happens when the required long simulation times are unattainable), Boltzmann exponential averaging of slow growth results should be used to produce an upper bound on the free energy change (Jarzynski, C., Phys. Rev. Lett., 78, 2690-2693, 1997), with exponential averaging of results for change in the opposite direction giving a lower bound; it is then reasonable to choose the mean of the bounds as the best estimate. The work, W sg for transfer of benzamidine from water to vacuum has been computed by insertion and extraction simulations, at different switching times. (As implemented, the molecular transformation calculation requires two evaluations of the Ewald sum; the increase in computer time required for this has been reduced by use of a multiple time step scheme in which the Ewald summations are executed at intervals of several integration time steps.) For the longest switching time, the distribution of values of W sg is narrow, hysteresis is small and all methods produce a similar result for ΔA, As the switching time is reduced, (i) the distribution becomes non-Gaussian, (ii) the frictional portions and the distributions for insertion and extraction differ, (iii) the mean of the linear averages and the mean of the exponential averages for insertion and extraction both fail to give an accurate estimate of ΔA. © 2002 Taylor & Francis Ltd.