Simulation Acceleration

From AlchemistryWiki
Jump to: navigation, search

In many cases of interest, carrying out robust free energy calculations may require significant investment of computational resources, beyond that which can be obtained by most researchers. This page examines additional tools for accelerating the sampling. As these begin to branch out from the fundamentals of free energy, we will not go deeply into all these methods. They are not needed to carry out free energy calculations, but they may be extremely useful to converge calculations in complex systems with slow dynamics and we encourage readers to examine these topics to see what will work for them.

Umbrella Sampling

One standard method for improving sampling in atomistic simulations is umbrella sampling,[1] where bias terms are added to constrain the simulation in some way, then the restraints are removed. This method can lower barriers on the potential energy surface, or restrain simulations to slow-interconverting configurations that are relevant to the binding affinities (e.g. different torsional states). Doing this allows for the free energy components to be properly computed and then combined.[2][3][4] If a system has slow degrees of freedom, like some hydration free energies,[5] this method allows one to sample more frequently in the slow state.

One can also envision using umbrella sampling to compute the free energy of constraining the free ligand into the bound conformation directly before computing the free energy of binding, and then computing the free energy of releasing these restraints in solution. Doing so usually decreases correlation times for sampling intermediate states, thereby increasing the simulation efficiency. [2][6]

Hamiltonian Replica Exchange

Although it is perfectly valid to run each intermediate state as its own simulation, there are tools that let you run parallel, coupled simulations and swap information between them to improve sampling. One such method is called Hamiltonian Exchange (HEX) or Hamiltonian Replica Exchange (HREX). In HEX, the parallel simulations at each alchemical variable are allowed to swap atoms/molecules (under certain conditions) from state A which has different energy barriers than state B. This enables sampling configurations that may take significant amounts of time to access normally when simulating only one state. [7][8][9][10][11][12]

We highly recommend doing HEX methods as most codes do this on top of well-tested temperature replica exchange methods, and so there is very little overhead. the outputs of Hamiltonian exchange simulations can be analyzed in the same way as the outputs of [math]K[/math] uncoupled simulations. These simulations are guaranteed to decorrelate as fast or faster than standard simulations, though the exact amount of improvement depends on the system; at least they will not be any worse!

Correlation times with HEX is a bit more complicated since atoms/molecules are swapping along very different trajectories. In this case, you will need to compute along trajectories that sample different states as opposed to a single state.

Some Suggestions for Replica Exchange Simulations

  • The level of acceleration one gets is a function both of phase space mixing between the replicas and the acceleration provided by of replicas kinetics that have faster sampling than the original distribution in the degree of freedom of interest.
  • Choice of Hamiltonian/temperature states depends quite a bit on what one actually wants to do: is it just accelerating kinetics, or does one actually care about properties as a function of temperature or calculating free energies? That can change the choice. Identifying a choice of replica that accelerates kinetics in desired way is nontrivial problem.
  • Space replicas to be even in phase space overlap. The easiest way to do this is to guess a spacing, then look at the results pf MBAR /TI after a short run. You can then respace the replicas so the error in the free energies between consecutive states is close to being equal. It is most important to get rid of very bad spacing; there is not much efficiency gain in going from decent spacing to optimal spacing [13].
  • Acceptance rates should be around 38% to maximize mixing in state space, but anywhere in the 25-50% range is nearly as good, so it is not that useful to tune it too much [14].
    • More precise quantitative measurements of state space mixing are determined by taking the transition matrix between states, and maximizing [math]1-\lambda_2[/math], where [math]\lambda_2[/math] is the second largest eigenvalue of the transition matrix. [15] This quantity is known as the spectral gap, and is related to the rate of mixing in state space. Remember that mixing in state space is insufficient to accelerate sampling; there must be some state where configurations are sampled more quickly.
  • To improve mixing between states, perform multiple replica switches each time one pauses to exchange. [15] Anecdotal evidence is that effectiveness of multiple switches decreases past 10-100 exchanges, but this is system dependent.
  • Exchange should be as frequent as is permissible given the networking and code [15], as long as the velocities are not reassigned on switches. Reassigning velocities too frequently slows kinetics and negates the advantages gained by performing rapid replica switches [16].

Other Methods

The other methods that will be mentioned briefly are Expanded Ensemble and [math]\lambda[/math]-dynamics. Expanded ensemble works by running a single simulation and sampling both the intermediate states and separate coordinates simultaneously. [math]\lambda[/math]-dynamics works by treating the alchemical variable as dynamical, introducing a fictitious mass corresponding to the [math]\lambda[/math] degree of freedom; although this seems rather novel, it is essentially identical to Monte Carlo techniques.[10][17][18][19] [math]\lambda[/math]-dynamics is showing promise, but is still in preliminary stages of development[20][21][22][19][23]

In general, we do NOT recommend Expanded Ensemble or [math]\lambda[/math]-dynamics to beginners. The methods and implementations are not up to the same robustness as HEX yet and there are tweaks and extra parameters that have to be coded to obtain proper convergence. Given more time and development, these methods may become more accessible to beginners and we will be able to recommend them in the future.

Slow Growth

In slow-growth methods, the coupling parameter, [math]\lambda[/math], is slowly varied over the course of one or several simulations to take a system gradually from [math]\lambda=0[/math] to [math]\lambda=1[/math]; from this, the free energy difference between endpoints is estimated. In equilibrium methods, on the other hand, separate simulations are run at multiple different lambda values and information from the individual simulations is used to estimate free energy differences between adjoining lambda values.

Slow-growth methods have a number of well-documented problems, such as Hamiltonian lag and hysteresis (add references). Additionally, if a slow-growth calculation does not lead to a sufficiently precise free energy estimate, the only way to improve the free energy estimate is to repeat the calculation with a slower rate of change in lambda – the simulation cannot be extended to gain additional precision, meaning that significant data can be wasted. Fast-growth methods, while conceptually appealing, do not appear to offer substantial advantages over equilibrium methods.[24][25]

Nonequilibrium (Fast Growth) methds

Fast-growth methods are based on a theorem by Jarzynski in 1997[26] that the free energy difference associated with a particular equilibrium process can be computed by taking an the exponential average of the irreversible work done in performing many (potentially rapid) non-equilibrium trials of the same process.

[math] \Delta G = -k_B T \ln \left \langle e^{-W/k_b T} \right \rangle [/math]

Where [math] W[/math] is the work to execute the transition between the initial and final thermodynamic state. When applied to alchemical free energy calculations, this simply amounts to estimating free energy differences by performing many different rapid versions of a slow-growth trial – that is, rapidly changing lambda between two values (i.e. 0 and 1) in many separate trials, and monitoring the work done each time. Equilibrium free energy calculations can be thought of as "instantaneous growth" as they rely on estimating the free energy difference between neighboring [math]\lambda[/math] values based only on instantaneous evaluations of potential energy differences between [math]\lambda[/math] values (or by evaluation of and extrapolation of [math]dV/d\lambda[/math] at a particular [math]\lambda[/math] value).


  1. Torrie, G. M., and Valleau, J. P. (1977) Non-physical Sampling Distributions in Monte-Carlo Free-Energy Estimation : Umbrella Sampling. J. Comput. Phys. 23, 187–199. - Find at Cite-U-Like
  2. 2.0 2.1 Wang, J., Deng, Y., and Roux, B. (2006) Absolute Binding Free Energy Calculations Using Molecular Dynamics Simulations with Restraining Potentials. Biophys. J. 91, 2798–2814. - Find at Cite-U-Like
  3. Mobley, D. L., Chodera, J. D., and Dill, K. A. (2007) Confine-and-release method: Obtaining correct binding free energies in the presence of protein conformational change. J. Chem. Theory Comput. 3, 1231–1235. - Find at Cite-U-Like
  4. Mobley, D. L., Graves, A. P., Chodera, J. D., McReynolds, A. C., Shoichet, B. K., and Dill, K. A. (2007) Predicting absolute ligand binding free energies to a simple model site. J. Mol. Biol. 371, 1118–1134.\
  5. Klimovich, P. V., and Mobley, D. L. (2010) Predicting hydration free energies using all-atom molecular dynamics simulations and multiple starting conformations. J. Comp. Aided Mol. Design 24, 307–316. - Find at Cite-U-Like
  6. Woo, H.-J., and Roux, B. (2005) Calculation of absolute protein-ligand binding free energy from computer simulation. Proc. Natl. Acad. Sci. 102, 6825–6830. - Find at Cite-U-Like
  7. Okamoto, Y. (2004) Generalized-ensemble algorithms: Enhanced sampling techniques for Monte Carlo and molecular dynamics simulations. J. Mol. Graph. Model. 22, 425–439. - Find at Cite-U-Like
  8. Roux, B., and Faraldo-Gómez, J. D. (2007) Characterization of conformational equilibria through Hamiltonian and temperature replica- exchange simulations: Assessing entropic and environmental effects. J. Comput. Chem. 28, 1634–1647. - Find at Cite-U-Like
  9. Woods, C. J., Essex, J. W., and King, M. A. (2003) Enhanced Configurational Sampling in Binding Free Energy Calculations. J. Phys. Chem. B 107, 13711–13718. - Find at Cite-U-Like
  10. 10.0 10.1 Banba, S., Guo, Z., and Brooks III, C. L. (2000) Efficient sampling of ligand orientations and conformations in free energy calculations using the lambda-dynamics method. J. Phys. Chem. B 104, 6903–6910. - Find at Cite-U-Like
  11. Bitetti-Putzer, R., Yang, W., and Karplus, M. (2003) Generalized ensembles serve to improve the convergence of free energy simulations. Chem. Phys. Lett. 377, 633–641. - Find at Cite-U-Like
  12. Hritz, J., and Oostenbrink, C. (2008) Hamiltonian replica exchange molecular dynamics using soft-core interactions. J. Chem. Phys. 128, 144121. - Find at Cite-U-Like
  13. Naden, L. N., Pham T. T., and Shirts, M. R. (2014) Linear Basis Function Approach to Efficient Alchemical Free Energy Calculations. 1. Removal of Uncharged Atomic Sites. J. Chem. Theory Comput." 10, 1128-1149 - Find at Cite-U-Like
  14. Predescu, C., Predescu, M., Ciobanu, C. V. (2005) On the Efficiency of Exchange in Parallel Tempering Monte Carlo Simulations J. Phys. Chem. B, 109, 4189-4196 - Find at Cite-U-Like
  15. 15.0 15.1 15.2 Chodera, J.D. and Shirts, M. R. (2011) Replica exchange and expanded ensemble simulations as Gibbs sampling: Simple improvements for enhanced mixing. J. Chem. Phys. 135, 194110 - Find at Cite-U-Like
  16. Basconi, J. E. and Shirts, M. R. (2013) Effects of Temperature Control Algorithms on Transport Properties and Kinetics in Molecular Dynamics Simulations. J. Chem. Theory Comput. 9, 2887-2899 - Find at Cite-U-Like
  17. Guo, Z., Brooks III, C. L., and Kong, X. (1998) Efficient and flexible algorithm for free energy calculations using the [math]\lambda[/math]-dynamics approach. J. Phys. Chem. B 102, 2032–2036. - Find at Cite-U-Like
  18. Kong, X., and Brooks III, C. L. (1996) [math]\lambda[/math]-dynamics: A new approach to free energy calculations. J. Chem. Phys. 105, 2414–2423. - Find at Cite-U-Like
  19. 19.0 19.1 Li, H., and Yang, W. (2007) Forging the missing link in free energy estimations: lambda-WHAM in thermodynamic integration, overlap histogramming, and free energy perturbation. Chem. Phys. Lett. 440, 155–159. - Find at Cite-U-Like
  20. Zheng, L., Carbone, I. O., Lugovskoy, A., Berg, B. A., and Yang, W. (2008) A hybrid recursion method to robustly ensure convergence efficiencies in the simulated scaling based free energy simulations. J. Chem. Phys. 129, 034105. - Find at Cite-U-Like
  21. Zheng, L., and Yang, W. (2008) Essential energy space random walks to accelerate molecular dynamics simulations: Convergence improvements via an adaptive-length selfhealing strategy. J. Chem. Phys. 129, 014105.
  22. Min, D., and Yang, W. (2008) Energy difference space random walk to achieve fast free energy calculations. J. Chem. Phys. 128, 191102.
  23. Min, D., Li, H., Li, G., Bitetti-Putzer, R., and Yang, W. (2007) Synergistic approach to improve “alchemical” free energy calculation in rugged energy surface. J. Chem. Phys. 126, 144109. - Find at Cite-U-Like
  24. Jarzynski, C. Phys. Rev. E 2006, 73, 046105. - Find at Cite-U-Like
  25. Oostenbrink, C.; van Gunsteren, W. F. Chemical Physics 2006, 323. - Find at Cite-U-Like
  26. Jarzynski, C. Physical Review Letters 1997, 78. - Find at Cite-U-Like