Difference between revisions of "Best Practices"

From AlchemistryWiki
Jump to navigation Jump to search
Line 1: Line 1:
= Best Practices =
+
{{Fundamentals | cTopic=Theory}}
Here, we provide an overview of some of what we believe are the best practices for free energy calculations, with references, wherever possible. This document assumes you already have a basic idea of what these calculations are and what they do (if not, start with our FreeEnergyIntroduction), and that you already have a working knowledge of molecular simulations and basic terms like convergence and equilibration (if not, start with a textbook like Leach's [http://www.google.com/search?hl=en&lr=&client=safari&rls=en&q=Leach+Molecular+Modelling:+Principles+and+Applications&btnG=Search Molecular Modelling: Principles and Applications]). Simulation-package specific issues will not be addressed here.
 
  
Feel free to edit this document or the attached discussion pages. We hope this to be a place storing the community consensus on these issues. Do back up any edits with appropriate references.
+
To understand the advantages and disadvantages of different free energy methods, it is important to begin with a review of the underlying principles. This page is dedicated to the most fundamental concepts of free energy calculations and is designed to give an in-depth view of the approaches, starting from the basics. This page also contains some of the common nomenclature and symbols that are seen throughout the rest of the [[:Category:Free Energy Fundamentals|free energy fundamentals]] pages and in the literature as a whole.
  
== Introduction ==
+
This page is not meant to be an end-all repository of the background mathematics and principals required for free energy calculations, but it will serve as a good start point and hopefully a quick reference.
Free energy calculations are appealing, in that, in principle, they allow calculation of rigorously correct free energies, '''given a particular set of parameters and physical assumptions'''. We firmly believe that the goal of these calculations, then, should be to obtain these correct free energies given the particular assumptions -- and ''not'' necessarily to match experiment. Only then can the underlying parameters and physical assumptions really be tested. In other words, there exists a single ''right'' answer for a free energy calculation, given this particular set of parameters and physical assumptions, and our goal is to obtain it. Here, we will call these free energies "correct" free energies. "Correct", then, implies at the very least that the computed free energies have converged, and that there are no underlying methodological problems.
 
  
It is also important to remember in what follows that free energy is a function of state, so there are many possible choices of pathway for a thermodynamic cycle connecting the same two endpoints. Some pathways may be more efficient than others (sometimes by many orders of magnitude) so methods which are in principle correct may not always be practical.
+
=Why the name "Alchemical"?=
 +
One of the first questions, and often most confusing points, about a number of free energy calculations is why we refer to them as "alchemical" changes in a large number of computational methods. When people first hear the word "alchemy," they usually think of the medieval alchemists who were trying to turn lead into gold, or other such transformations. How computational free energy adopted the name takes a bit of understanding of simulation limitations and the properties of what we are calculating.
  
Here, we discuss best practices in the context of several practical examples: Solvation free energies, and binding free energies. Solvation free energies provide a basic starting point for considering a number of the key issues, and, for binding free energies, there are many more choices of pathway, which raises additional complications. In both cases, there are also different methods for computing the relevant free energy differences, which may differ in efficiency.
+
Considering that the natural evolution of some of the processes we try to model are well beyond reasonable simulation time scales, we must come up with very efficient ways to compute the free energy differences. Because free energy is a state variable (path independent), we can design simulations that provide a convenient pathway to computing free energy. Furthermore, since we are doing simulations, we are not limited by experimentally observable conditions so long as we are meticulous with our calculations.
  
== Guidelines for free energy calculation pathways ==
+
With these observations in mind, it was found that we can simulate modifications to atoms which can change their properties to reflect non-physical or entirely different species. This is roughly the same definition of alchemy from old in that we are changing the properties of the atoms to be something else, although we do it in a mathematically sound and rigorous manner; hence, the term "alchemical" was coined as many free energy calculations are "like alchemy" in their pathways and methods.
Alchemical free energies almost involve some insertion or deletion of atoms -- both in hydration free energy calculations, and in absolute and relative binding free energy calculations. By insertion and deletion, we mean decoupling or annihilation (see [[Decoupling and annihilation]] for definitions) of the interactions of the atoms in question. We believe that a basic list of rules and guidelines should be followed in any calculation that involves insertion or deletion:
 
  
* '''Rule 1''': Always use soft-core potentials while decoupling or annihilating Lennard-Jones interactions
+
There are of course many factors and checks that must be done to ensure accuracy and robustness of the calculations, many of which will be shown in the rest of the [[:Category:Free Energy Fundamentals|free energy fundamentals pages]].
* '''Rule 2: '''Never leave a partial atomic charge on an atom while its Lennard-Jones interactions are being removed
 
* '''Guideline 3: '''It is usually more efficient to perform electrostatic and Lennard-Jones transformations separately
 
* '''Guideline 4: '''Inserting or deleting atoms is usually less efficient than mutating them, so transformations should involve as few insertions and deletions as possible.
 
* '''Guideline 5: '''Keep configuration space in mind and think about convergence.
 
  
It is worth looking at each of these in more detail.
+
=Assumptions for the Fundamentals=
 +
The assumptions listed here are carried out through the rest of the fundamentals sections. Free energy calculations can usually be set up without these assumptions, but we'll make these assumptions to simplify the explanations.  
  
 +
* '''A standard molecular mechanics model will be assumed'''; this includes:
 +
** Harmonic bond angle terms
 +
** Periodic dihedral terms
 +
** Non-bonded terms made up of point charges and Lennard-Jones repulsion/dispersion terms
 +
* constant temperature is maintained:  (discussed [[#Facts from the Free Energy Difference Definition | below]])
 +
* '''Masses do not alchemically change.''' If one wishes to do this, substitute all potential energies, <math>U</math>, with the more general Hamiltonian, <math>\mathcal{H}</math>.
 +
* '''QM/MM will ''not'' be considered''' for fundamentals because of the field has yet to be sufficiently developed.{{Cite|Woods2008|>Woods, C. J., Manby, F. R., and Mulholland, A. J. (2008) An efficient method for the calculation of quantum mechanics/molecular mechanics free energies. J. Chem. Phys. 128, 014109.|http://www.citeulike.org/group/14929/article/9052667}}
  
=== Rule 1: Soft core potentials ===
+
=Free Energy Difference Equation=
Rule 1 -- that soft core potentials should always be used when turning on or off Lennard-Jones interactions -- really should be a rule observed by all free energy calculations, we believe. Lennard-Jones interactions between particles have a really steep (<math>1/r^{12}</math>) rise in potential energy. This prevents particles from overlapping. However, to delete particles (atoms or molecules) these interactions need to be turned off somehow, and this is not as straightforward as it might seem.
+
Most free energy calculations and methods started from a single core equation derived from statistical mechanics. The free energy difference between two states is directly related to the ratio of probabilities of those states through the partition functions. This free energy difference for an NVT ensemble is
  
One relatively-commonly used choice ([http://dx.doi.org/10.1098/rsta.2005.1625 Fowler et al., 2005], [http://dx.doi.org/10.1007/s10822-005-9021-3 Chipot, Rozanska and Dixit, 2005], [http://amber.scripps.edu/doc9/index.html AMBER], and others) for turning these off is simple linear-scaling of that term in the potential energy or Hamiltonian: That is,  <math>V(\lambda) = (1 - \lambda) V_0 + \lambda V_1</math>, where <math>V_0</math> is the potential energy with full Lennard-Jones interactions, and <math>V_1</math> is the potential energy where the Lennard-Jones interactions have been turned off for the atoms which are being deleted. This means that, for the atoms being deleted, Lennard-Jones interactions scale at small <math>r</math> as <math>(1-\lambda)/r^{12}</math>. This has two unfortunate and interconnected consequences. First, there is a discontinuous change in the form of the interaction potential when going from <math>\lambda=1-\epsilon</math> (where <math>\epsilon</math> is a very small number) to <math>\lambda=1</math>, as the <math>1/r^{12}</math> term still is fairly important even at <math>\lambda</math> very near 1, but is entirely turned off at <math>\lambda=1</math>. Secondly, it leads to large forces, numerical instabilities, and other problems in simulations near <math>\lambda=1</math>. Formally, it has been shown that this leads to a integrable singularity in <math>dV/(d\lambda)</math>, which means that computing correct free energies with this scheme using thermodynamic integration is impossible using numerical techniques ([http://dx.doi.org/10.1063/1.432264 Mruzik et al., 1976],[http://fulcrum.physbio.mssm.edu/~mezei/ms/cv54mmc.pdf Mezei and Beveridge, 1986], [http://fulcrum.physbio.mssm.edu/~mezei/ms/cv85.pdf Resat and Mezei, 1993] and especially [http://dx.doi.org/10.1016/0009-2614(94)00397-1 Beutler et al., 1994], [http://dx.doi.org/10.1080/08927020211973 Pitera and van Gunsteren, 2002] and [http://www.dillgroup.ucsf.edu/~dmobley/papers/steinbrecher.pdf Steinbrecher et al., 2007] and references therein) and similar problems plague free energy perturbation schemes.
+
:<math>\displaystyle \Delta A_{ij} = -k_B T \ln \frac{Q_j}{Q_i} = -k_B T \ln \frac{\int_{\Gamma_j} e^{-\frac{U_j(\vec{q})}{k_B T}} d\vec{q}}{\int_{\Gamma_i} e^{-\frac{U_i(\vec{q})}{k_B T}}d\vec{q}}</math>
  
In an attempt to get around this, some have suggested scaling the potential energies with <math>(1-\lambda)^k</math>, where <math>k</math> is an integer greater than 1. It can be shown that, for <math>k >= 4</math>, this leads to an integrable singularity in <math>dV/d\lambda</math>, so thermodynamic integration can in principle be done ([http://fulcrum.physbio.mssm.edu/~mezei/ms/cv54mmc.pdf Mezei and Beveridge, 1986], [http://www.sciencedirect.com/science?_ob=ArticleURL&_udi=B6TFN-44WCYB3-T8&_user=4430&_coverDate=06/03/1994&_rdoc=1&_fmt=&_orig=search&_sort=d&view=c&_acct=C000059594&_version=1&_urlVersion=0&_userid=4430&md5=43ac38d23e3e595dc97c273e012406e6 Beutler et al., 1994]). But integrable singularities still pose very substantial problems for molecular simulation, and this approach can still lead to large forces, numerical instabilities and energy conservation problems ([http://dx.doi.org/10.1016/0009-2614(94)00397-1 Beutler et al., 1994] and [http://www.dillgroup.ucsf.edu/~dmobley/papers/steinbrecher.pdf Steinbrecher et al., 2007]) and make free energy differences extremely difficult to converge ([D. Mobley, unpub. data]).
+
where <math>\Delta A_{ij}</math> is the Helmholtz free energy difference between state <math>j</math> and state <math>i</math>, <math>k_B</math> the Boltzmann constant, <math>Q</math> the canonical partition function, <math>T</math> is the temperature, <math>U_i</math> and <math>U_j</math> are the potential energies as a function of the coordinates and momenta <math>\vec{q}</math> for two states, and <math>\Gamma_i</math> and <math>\Gamma_j</math> are the ''phase space volumes'' of <math>\vec{q}</math> over which we sample, or the total set of all allowed positions and momenta of the system.
  
Since free energies are path-independent, an elegant solution to this problem was developed ([http://dx.doi.org/10.1016/0009-2614(94)00397-1 Beutler et al., 1994]) – to modify the Lennard-Jones functional form to gradually smooth out the <math>1/r^{12}</math> term as a function of <math>\lambda</math>, rather than simply multiplying it by a prefactor. This removes problems with numerical instabilities and singularities, and improves convergence properties ([http://dx.doi.org/10.1016/0009-2614(94)00397-1 Beutler et al., 1994], [http://dx.doi.org/10.1063/1.466707 Zacharias et al., 1994], [http://dx.doi.org/10.1080/08927020211973 Pitera and van Gunsteren, 2002]).  The basic idea is that it allows particles to gradually begin to overlap as <math>\lambda</math> is changed, rather than saving a drastic change in interactions for the point going from <math>\lambda=1-epsilon<math> to <math>\lambda=1</math>. This approach is known as soft core potentials (and, alternately, "separation-shifted scaling"), and has subsequently been shown to be a nearly optimal path for modifying Lennard-Jones interactions ([http://dx.doi.org/10.1002/jcc.20025 Blondel, 2004]. In some work, several groups have further tested this approach and found a slightly modified functional form and set of parameters from that originally proposed ([http://dx.doi.org/10.1016/0009-2614(94)00397-1 Beutler et al., 1994]) which leads to improved efficiency for free energy calculations ([http://dx.doi.org/10.1063/1.1877132 Shirts and Pande, 2005], [D. Mobley, unpublished data]); we recommend that the soft core potentials and parameters from that work be employed in all free energy calculations involving insertion or deletion of particles.
+
From this equation, all free energy calculations are derived.
  
Some testing has suggested that the <math>(1-\lambda)^k</math> scaling approach may be essentially adequate for hydration free energy calculations ([D. Mobley, unpublished data],[http://www.dillgroup.ucsf.edu/~dmobley/papers/steinbrecher.pdf Steinbrecher et al., 2007]) but it still less efficient there than soft-core potentials, so this does not affect our recommendation.
+
=Nomenclature and Variables=
 +
In an effort to keep things uniform, this section contains all the common constants and variables that are seen throughout the [[:Category:Free Energy Fundamentals|fundamentals sections.]] Although the Table itself is not in any particular order, the sorting buttons should help you find what you are looking for. The alternate listings of some variables are things one may encounter in literature, although this site will strive to be uniform in its naming process.  
  
In summary: Linearly scaling Lennard-Jones interactions back as a function of <math>\lambda</math> for insertion/deletion of particles is formally incorrect for numerical integration and leads to wrong estimates of free energy differences. While more complicated schemes involving <math>\lambda^k</math> scaling can be formally correct, there are serious concerns regarding their accuracy. Soft-core potentials provide a rigorously correct, efficient alternative to these and should be used whenever particles are inserted or deleted, preferably with the functional form and parameters of (([http://dx.doi.org/10.1063/1.1877132 Shirts and Pande, 2005]), unless future work finds a still more efficient set of parameters.
+
{|class="wikitable sortable" style="text-align: center"
 +
|-
 +
! scope="col" | Variable, Acronym, Term
 +
! scope="col" | Definition
 +
! scope="col" class="unsortable" |Notes
 +
|-
 +
| <math>k_B</math>
 +
| [http://en.wikipedia.org/wiki/Boltzmann_constant Boltzmann's Constant]
 +
|
 +
|-
 +
| <math>U</math>
 +
| Potential Energy
 +
| This is the total internal energy  calculated a classical simulation, the sum of the potential and kinetic energies.
 +
|
 +
|-
 +
| <math>T</math>
 +
| Temperature
 +
| This is often shown with <math>k_B</math> as well, and so common in fact that there is a common variable affiliated with it: <math>\beta</math>.
 +
|-
 +
| <math>\beta</math>
 +
| Inverse Boltzmann Temperature <math>\beta=(k_BT)^{-1}</math>
 +
| Because it is so common in most the equations, one will see <math>\beta^{-1}</math> more often than <math>k_BT</math>
 +
|-
 +
| <math>P</math>
 +
| Pressure
 +
|
 +
|-
 +
| <math>u</math>
 +
| Reduced potential, expressed as <math>u=\beta(U+PV-\sum \mu N_i)</math>. For many equations, it is simpler to work with the reduced potential.
 +
| In the case of NVT ensembles, the PV and <math>\mu N_i</math> terms are zero. For NPT ensembles, the <math>\sum \mu N_i</math> is zero, and for <math>NV\mu</math> ensembles, the PV term is zero.
 +
|-
  
===Rule 2: Turn off partial charges===
+
| <math>V</math>
 +
| Volume (traditional)
 +
| '''CAUTION:''' Volume in this sense is the kind you think of with box size. Not to be confused with the ''Phase Space Volume.'' This particular volume is not often considered by itself and should only show up on this site within <math>PV</math> terms, or ''NVT''.
 +
|-
 +
| <math>\Gamma</math>
 +
| Phase Space Volume: Total set of coordinates and momenta where the system has a nonzero probability of being found
 +
| '''CAUTION:''' This is not volume in standard physical sense, meaning the volume of the box, but the volume space times the momentum space.
 +
|-
 +
| <math>q</math> or <math>\vec{q}</math>
 +
| Collective variable for coordinates and momentum.
 +
| Can also be seen as <math>\mathbf{x}</math>, <math>\mathbf{\vec{x}}</math>, and a number of other terms in literature.
 +
|-
 +
| <math>\epsilon</math> and <math>\sigma</math>
 +
| Lennard-Jones parameters for the well depth and the radius of the point particle respectively.
 +
|
 +
|-
 +
| <math>Q</math>
 +
| Partition function for the canonical or NVT ensemble, but sometimes also signifies an arbitrary partition function of any ensemble. <br/><math>Q=\int_{\Gamma}e^{-\beta U (\vec{q})}\,\mathrm{d}\vec{q}</math>
 +
| In the examples here, this will represent the canonical partition function although most of the theory shown here works for an arbitrary partition function as well, though sometimes with slight modifications, where <math>Q=\int_{\Gamma}e^{-u(\vec{q})}\,\mathrm{d}\vec{q}</math>
 +
|-
 +
| <math>\tau</math>
 +
| [[Simulation_Information_Gathering#Correlation|Autocorrelation time]]
 +
|
 +
|-
 +
| ''NVT''
 +
| System held at constant '''N'''umber of particles, '''V'''olume, and '''T'''emperature. Free energy calculations performed in this ensemble yield Helmholtz free energies.
 +
|
 +
|-
 +
| ''NPT''
 +
| System held at constant '''N'''umber of particles, '''P'''ressue, and '''T'''emperature. Free energy calculations performed in this ensemble yield Gibbs free energies.
 +
|
 +
|-
 +
| ''TI''
 +
| [[Thermodynamic Integration]]
 +
|
 +
|-
 +
| ''EXP''
 +
| [[Exponential Averaging]]
 +
| Also called "free energy perturbation" and the "Zwanzig relationship."
 +
|-
 +
| ''BAR''
 +
| [[Bennett Acceptance Ratio]]
 +
|
 +
|-
 +
| ''MBAR''
 +
| [[Multistate Bennett Acceptance Ratio]]
 +
|
 +
|-
 +
| Boltzmann Weight
 +
| The probability of a microstate <math>\vec{q}</math>, defined <math>e^{-u(\vec{q})}/Q</math>, where <math>u(\vec{q})</math> is the reduced potential.  It can also refer to simply <math>e^{-u(\vec{q})}</math>, although this is not a full weight because it is not normalized by the partition function. 
 +
|-
 +
| State (variable <math>K</math>)
 +
| A unique set of conditions and parameters that completely describe a thermodynamic system.
 +
| For all examples here, the upper case variable <math>K</math> will denote all states (or the count thereof), and the lower case variable <math>k</math> will refer to a specific state.
 +
|-
 +
| <math>\lambda</math>
 +
| Alchemical variable, thermodynamic path coordinate
 +
| The most common symbol for denoting/tracking the progress along some alchemical pathway. Since most free energy calculations involve an [[Intermediate States | alchemical thermodynamic path]], [[Main Page| this site]] was named accordingly.
 +
|-
 +
| Annihilating/Annihilated
 +
| The act of removing ALL of an atoms interactions with the system, this includes the molecule's interactions with itself.
 +
| Intermolecular force = OFF, Intramolecular forces = OFF. Also known as disappeared or destroyed. Not to be confused with "decoupling."
 +
|-
 +
| Decoupling
 +
| The act of removing an atom or molecules interactions with its surroundings ONLY; the interactions between the molecule itself remain active
 +
| Intermolecular force = OFF, Intramolecular forces = ON. Not to be confused with "annihilation."
 +
|}
  
Rule 2 states that a partial atomic charge should never be allowed to remain on an atom while its Lennard-Jones interactions are being removed. To understand the reason for this, consider two atoms of opposite charge, A and B. Lennard-Jones interactions of atom A are being scaled back. Regardless of the scaling scheme used, at some lambda value, atoms A and B will begin to overlap occasionally, since the final state allows A and B to overlap totally. If A has a remaining partial atomic charge when these overlaps become possible, the two point charges assigned to A and B can actually overlap as well. Since the potential energy of Coulomb interactions between point charges scales as <math>q_{A} q_{B}/r_{AB}</math>, where <math>r_{AB}</math> is the distance between A and B, this presents significant problems when <math>q_A</math> and <math>q_B</math> have opposite signs. In particular, there is an infinite energy minimum at <math>r_{AB}=0</math>, so the two particles would in principle get trapped on top of one another.  
+
==Facts from the Free Energy Difference Definition==
 +
There are a number of key notes we can learn from the definition of free energy differences. Each of these can be important  in interpreting simulation results.
  
In practice, what usually happens in molecular dynamics simulations in these circumstances is that the forces get extremely large as A and B begin to overlap, and the simulation will crash. Constraint algorithms are often the first to fail, so this may lead to a warning about constraints (i.e. LINCS or SHAKE) and then a crash. This issue is discussed briefly by [http://dx.doi.org/10.1080/08927020211973 Pitera and van Gunsteren] and in more detail by [http://dx.doi.org/10.1063/1.1924449 Anwar and Heyes].
+
# '''Only free energy ''differences'' are ever calculated.''' There is never a calculation where absolute free energies are needed (and rarely can they be calculated at all) as all of the biological or thermodynamical quantities of interest are based on a free energy difference. As such, there must always be a minimum of two defined thermodynamic states.
 +
#* Even ''absolute'' free energies of binding are still free energy differences between two states: the ligand restricted to the binding site, and the ligand free to explore all other configurations.
 +
# '''Free energy differences between states at different temperatures are usually not what you want to be calculating''' for problems of interest.  If it did, you would get <math>\Delta A_{ij} = -k_B T_j \ln Q_j + k_B T_i \ln Q_i</math>, which is no longer a ratio calculation and not needed for biological systems of interest. Temperature dependence on free energy is more likely to be "what is <math> \Delta A_{ij} </math> change at two different temperatures?"
 +
# '''There are two different phase space volumes.''' <math>\Gamma_i</math> and <math>\Gamma_j</math> are often the same, but they are not required to be. The methods presented here almost always assumes that the two phase spaces overlap. However, when the spaces do not overlap, these methods break down and it is difficult to identify this problem without in depth knowledge of your system. Consider the example of a hard sphere solute with radius <math>\sigma</math> at state <math>i</math> and a Lennard Jones repulsion/dispersion potential, with the same <math>\sigma</math> at state <math>j</math>. Since <math>\Gamma _i</math> will not have molecules at a distance less than <math>\sigma</math>, but <math>\Gamma_j</math> will, the two phase spaces are not the same and these methods will either break down or return very error-prone results.
 +
#* The degree to which the phase spaces are shared is called the "phase space overlap". Efficient free energy calculations require significant phase space overlap. There are  [[Intermediate_States | a number of strategies to address]] lack of overlap between target spaces, but determining the best way for any given situation is still a research question.
 +
#* It should also be noted that "near zero probability" and "always zero probability" are two distinct things when considering phase space. So long as there is a chance for an observation to be made, no matter how small, it is considered part of the phase space.
  
In view of this problem, we recommend always turning off partial charges for any atoms for which Lennard-Jones interactions are being removed before doing the Lennard-Jones transformation. Additionally, when Lennard-Jones parameters for an atom are being substantially modified during a free energy calculation (i.e. for relative free energy calculations involving mutation of an atom) and soft-core potentials are employed, similar problems may arise, so it may be useful to remove partial charges on atoms which are being mutated, as well.
+
=References=
 +
<References />
  
Several groups have developed modified electrostatics scaling methods in an attempt to bypass this problem and allow electrostatics interactions and Lennard-Jones interactions to be turned off in only one set of calculations (for example, [http://dx.doi.org/10.1063/1.1924449 Anwar and Heyes]), but since electrostatics transformations are usually so smooth a function of <math>\lambda</math> and need only few <math>\lambda</math> values for good overlap ([http://dx.doi.org/10.1063/1.1877132 Shirts et al., 2005]; [http://dx.doi.org/10.1021/jp0667442 Mobley et al., 2007], and others) it is unclear that this results in any significant efficiency gain over performing the transformations separately.
+
[[Category:Free Energy Fundamentals]]
 
 
In view of this, our recommendation is that either (a) partial charges on any particles being inserted or deleted be turned off prior to the use of soft core potentials for those particles, or (b) a soft core scheme for electrostatics be implemented to allow simultaneous changes.
 
 
 
===Guideline 3: Perform electrostatics transformations separately from Lennard-Jones ===
 
 
 
As noted in Rule 2, above, electrostatics transformations are typically smooth functions of lambda with good phase-space overlap between even coarsely-spaced lambda values([http://dx.doi.org/10.1063/1.1877132 Shirts et al., 2005]; [http://dx.doi.org/10.1021/jp0667442 Mobley et al., 2007], and others). As a consequence, these are quite efficient compared to Lennard-Jones calculations. As established above, when particles are being inserted or deleted, the electrostatic interactions of these particles should be set to zero before turning off their Lennard-Jones interactions. But what about electrostatic interactions on atoms which are merely being mutated (i.e. a change of partial charge and Lennard-Jones radius), as in relative free energy calculations?
 
 
 
We are not aware of any study which has looked at this in detail, but given the efficiency of free energy calculations modifying electrostatics interactions relative to those significantly modifying Lennard-Jones interactions, we believe it makes sense to perform the two sets of calculations separately. Given that the two transformations have different lambda-dependences, it might actually be less efficient to perform them together than separately. Performing them separately has an additional advantage, as well: Uncertainties in the two components can be assessed separately, and computational effort focused on reducing the largest uncertainty (i.e. by extending some simulations to get additional sampling).
 
 
 
Further testing should be focused in this area, to determine whether alternative scaling approaches which can modify Lennard-Jones and electrostatic interactions simultaneously ([http://dx.doi.org/10.1063/1.1924449 Anwar and Heyes]) are actually more efficient than the approach of separate modification that we propose.
 
 
 
 
 
===Guideline 4: Use few insertions/deletions===
 
 
 
 
 
Electrostatics transformations are usually smooth functions of lambda (without soft core potentials), and require few lambda values, while Lennard-Jones transformations – especially those involving insertion and deletions – are difficult transformations which require substantially more lambda values (even when using soft core potentials) to obtain good phase-space overlap and accurate free energy differences ([http://dx.doi.org/10.1063/1.1877132 Shirts et al., 2005]; [http://dx.doi.org/10.1021/jp0667442 Mobley et al., 2007], [http://dx.doi.org/10.1016/j.jmb.2007.06.002 Mobley et al., 2007b], and others). Thus, insertions and deletions of particles can be thought of as “difficult” transformations (i.e. [http://dx.doi.org/10.1103/PhysRevE.73.046105 Jarzynski, 2006]). Consequently, it is far more computationally efficient to modify existing particles (atoms) than to insert or delete new atoms; this should be kept in mind when constructing mutation pathways for relative free energy calculations, since multiple choices of mutation pathways between a set of molecules are typically possible.
 
 
 
This guideline is not at all helpful for absolute free energy calculations, since these by design involve inserting or deleting entire molecules.
 
 
 
 
 
===Guideline 5: Think about configuration space and convergence===
 
 
 
Given that many choices of pathway are possible, it can often be helpful to think about whether a particular choice of pathway makes convergence easier or more difficult.
 
 
 
For example, in absolute binding free energies, one can incorporate the standard state using either simple distance restraints between the ligand and the protein, or by restraining the ligand orientation as well. At the fully noninteracting state, the amount of configuration space the ligand will need to sample is dictated by this choice. Hence, a ligand with only a single reference distance restrained relative to the protein will need to sample a spherical shell in configuration space, while a ligand with all six relative degrees of freedom restrained would need to sample only a very small region of configuration space. These two can take drastically different amounts of time, so in fact it can be much more efficient, at least in some cases, to use the additional restraints ([http://dx.doi.org/10.1063/1.2221683 Mobley et al., 2006]).
 
 
 
==Performing free energy calculations==
 
 
 
Here, we attempt to separate methods for analyzing the results of free energy calculations, from methods for performing free energy calculations. For example, thermodynamic integration and exponential averaging (among others) are methods for analyzing free energy calculations, while equilibrium, slow-growth, and fast-growth methods can be used for performing free energy calculations. In this section, we focus on methods for performing free energy calculations, and address analysis methods below.
 
 
 
There are several basic groups of methods for performing free energy calculations: Slow-growth, fast-growth, and equilibrium (or instantaneous growth) free energy methods. Slow-growth and equilibrium methods are more traditional, while fast-growth (non-equilibrium) methods have gained considerable recent interest. 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=</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. Fast-growth methods are based on the demonstration by [http://dx.doi.org/10.1103/PhysRevLett.78.2690 Jarzynski] in 1997 that the free energy difference associated with a particular equilibrium process can be computed by taking an appropriate average of the irreversible work done in performing many (potentially rapid) non-equilibrium trials of the same process. 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).
 
 
 
Which method should be used for alchemical free energy calculations? At this point, we believe the evidence is in favor of equilibrium methods. 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 ([http://dx.doi.org/10.1103/PhysRevE.73.046105 Jarzynski, 2006], [http://dx.doi.org/10.1016/j.chemphys.2005.08.054 Oostenbrink and van Gunsteren, 2006]).
 
 
 
In view of these facts, our recommendation is to use equilibrium simulations at a set of separate <math>\lambda</math> values to estimate free energy differences. But how should these simulations be performed? Should a <math>lambda=0</math> simulation be performed first, for example, and then the ending conformation used to seed a <math>lambda=0.1</math> simulation? Or should they be independent? There is no reason in principle that simulations cannot be performed serially – but this leads to several potential pitfalls and disadvantages. First, if simulations are performed serially, and it is later concluded that some of the intermediate simulations were not long enough, the entire set of calculations performed after those intermediate simulations must be repeated. Second, results at one <math>\lambda</math> value become coupled to those at previous lambda values, potentially leading to hysteresis (''add references'') and difficulties in troubleshooting. Independent simulations, on the other hand, provide the great practical advantage that if, at the analysis stage, it appears that data is insufficient at some range of <math>\lambda</math> values, these particular simulations can simply be extended without requiring repetition of other component simulations. Additionally, with large computer clusters now becoming relatively common, independent simulations can often be performed in parallel, leading to significant real-time savings over running the simulations serially.
 
 
 
One requirement for running equilibrium simulations at independent lambda values is that simulations must separately be equilibrated at each lambda value. Significant literature exists on assessing equilibration in molecular dynamics simulations (''add references''), so we refer the reader there. Suffice it to say that meaningful results depend crucially on having adequate equilibration at each lambda value.
 
 
 
Some workers have attempted to address the need for equilibration by performing a long (many nanosecond) pre-equilibration at <math>\lambda=0</math> before beginning free energy calculations at a variety of lambda values ([http://dx.doi.org/10.1063/1.1999637 Fujitani et al.]). It is possible that this approach is a mistake. For example, if this pre-equilibration is substantially longer than simulation times, it may allow conformational changes in the system at <math>\lambda=0</math> which are then not sampled on simulation timescales, trapping the system in a configuration that (while appropriate at <math>\lambda=0</math>) may not be appropriate at other lambda values. For an example of what could happen, see ([http://dx.doi.org/10.1021/ct700032n Mobley et al., 2007b]), where the protein remains kinetically trapped in its starting conformational state throughout all component free energy calculations. Our recommendation is that, if long equilibration is deemed necessary, equilibration periods should be equally long at each lambda value.
 
 
 
=== Phase space overlap ===
 
 
 
Phase space overlap is an additional consideration in performing free energy calculations. Essentially, all of the equilibrium approaches for estimating free energy differences rely on having simulations run at separate lambda values that are sufficiently closely-spaced that adjoining lambda values sample similar configurations. There are three basic ways to attempt to assess whether this requirement is met, listed in order of ascending sophistication:
 
# Trial and error: For a particular system, begin with very closely spaced lambda values, and then gradually increase separation until results begin to deteriorate.
 
# Error analysis:
 
## Using autocorrelation analysis (i.e. [http://dx.doi.org/10.1021/ct0502864 Chodera et al., 2006])
 
## If necessary, using block bootstrap approaches (i.e. [http://dx.doi.org/10.1063/1.2221683 Mobley et al., 2006])
 
# Phase space overlap measures (i.e. [http://dx.doi.org/10.1063/1.1992483 Wu and Kofke, 2005]).
 
 
 
Different types of transformations may have different overlap properties. For example, above, we noted that electrostatics transformations are usually smooth, which, more rigorously, means that relatively good phase space overlap is maintained even with fairly coarsely-spaced lambda values. Insertion and deletion of particles, on the other hand, requires more closely-spaced lambda values to ensure good overlap.
 

Revision as of 21:20, 14 February 2013



To understand the advantages and disadvantages of different free energy methods, it is important to begin with a review of the underlying principles. This page is dedicated to the most fundamental concepts of free energy calculations and is designed to give an in-depth view of the approaches, starting from the basics. This page also contains some of the common nomenclature and symbols that are seen throughout the rest of the free energy fundamentals pages and in the literature as a whole.

This page is not meant to be an end-all repository of the background mathematics and principals required for free energy calculations, but it will serve as a good start point and hopefully a quick reference.

Why the name "Alchemical"?

One of the first questions, and often most confusing points, about a number of free energy calculations is why we refer to them as "alchemical" changes in a large number of computational methods. When people first hear the word "alchemy," they usually think of the medieval alchemists who were trying to turn lead into gold, or other such transformations. How computational free energy adopted the name takes a bit of understanding of simulation limitations and the properties of what we are calculating.

Considering that the natural evolution of some of the processes we try to model are well beyond reasonable simulation time scales, we must come up with very efficient ways to compute the free energy differences. Because free energy is a state variable (path independent), we can design simulations that provide a convenient pathway to computing free energy. Furthermore, since we are doing simulations, we are not limited by experimentally observable conditions so long as we are meticulous with our calculations.

With these observations in mind, it was found that we can simulate modifications to atoms which can change their properties to reflect non-physical or entirely different species. This is roughly the same definition of alchemy from old in that we are changing the properties of the atoms to be something else, although we do it in a mathematically sound and rigorous manner; hence, the term "alchemical" was coined as many free energy calculations are "like alchemy" in their pathways and methods.

There are of course many factors and checks that must be done to ensure accuracy and robustness of the calculations, many of which will be shown in the rest of the free energy fundamentals pages.

Assumptions for the Fundamentals

The assumptions listed here are carried out through the rest of the fundamentals sections. Free energy calculations can usually be set up without these assumptions, but we'll make these assumptions to simplify the explanations.

  • A standard molecular mechanics model will be assumed; this includes:
    • Harmonic bond angle terms
    • Periodic dihedral terms
    • Non-bonded terms made up of point charges and Lennard-Jones repulsion/dispersion terms
  • constant temperature is maintained: (discussed below)
  • Masses do not alchemically change. If one wishes to do this, substitute all potential energies, [math]\displaystyle{ U }[/math], with the more general Hamiltonian, [math]\displaystyle{ \mathcal{H} }[/math].
  • QM/MM will not be considered for fundamentals because of the field has yet to be sufficiently developed.[1]

Free Energy Difference Equation

Most free energy calculations and methods started from a single core equation derived from statistical mechanics. The free energy difference between two states is directly related to the ratio of probabilities of those states through the partition functions. This free energy difference for an NVT ensemble is

[math]\displaystyle{ \displaystyle \Delta A_{ij} = -k_B T \ln \frac{Q_j}{Q_i} = -k_B T \ln \frac{\int_{\Gamma_j} e^{-\frac{U_j(\vec{q})}{k_B T}} d\vec{q}}{\int_{\Gamma_i} e^{-\frac{U_i(\vec{q})}{k_B T}}d\vec{q}} }[/math]

where [math]\displaystyle{ \Delta A_{ij} }[/math] is the Helmholtz free energy difference between state [math]\displaystyle{ j }[/math] and state [math]\displaystyle{ i }[/math], [math]\displaystyle{ k_B }[/math] the Boltzmann constant, [math]\displaystyle{ Q }[/math] the canonical partition function, [math]\displaystyle{ T }[/math] is the temperature, [math]\displaystyle{ U_i }[/math] and [math]\displaystyle{ U_j }[/math] are the potential energies as a function of the coordinates and momenta [math]\displaystyle{ \vec{q} }[/math] for two states, and [math]\displaystyle{ \Gamma_i }[/math] and [math]\displaystyle{ \Gamma_j }[/math] are the phase space volumes of [math]\displaystyle{ \vec{q} }[/math] over which we sample, or the total set of all allowed positions and momenta of the system.

From this equation, all free energy calculations are derived.

Nomenclature and Variables

In an effort to keep things uniform, this section contains all the common constants and variables that are seen throughout the fundamentals sections. Although the Table itself is not in any particular order, the sorting buttons should help you find what you are looking for. The alternate listings of some variables are things one may encounter in literature, although this site will strive to be uniform in its naming process.

Variable, Acronym, Term Definition Notes
[math]\displaystyle{ k_B }[/math] Boltzmann's Constant
[math]\displaystyle{ U }[/math] Potential Energy This is the total internal energy calculated a classical simulation, the sum of the potential and kinetic energies.
[math]\displaystyle{ T }[/math] Temperature This is often shown with [math]\displaystyle{ k_B }[/math] as well, and so common in fact that there is a common variable affiliated with it: [math]\displaystyle{ \beta }[/math].
[math]\displaystyle{ \beta }[/math] Inverse Boltzmann Temperature [math]\displaystyle{ \beta=(k_BT)^{-1} }[/math] Because it is so common in most the equations, one will see [math]\displaystyle{ \beta^{-1} }[/math] more often than [math]\displaystyle{ k_BT }[/math]
[math]\displaystyle{ P }[/math] Pressure
[math]\displaystyle{ u }[/math] Reduced potential, expressed as [math]\displaystyle{ u=\beta(U+PV-\sum \mu N_i) }[/math]. For many equations, it is simpler to work with the reduced potential. In the case of NVT ensembles, the PV and [math]\displaystyle{ \mu N_i }[/math] terms are zero. For NPT ensembles, the [math]\displaystyle{ \sum \mu N_i }[/math] is zero, and for [math]\displaystyle{ NV\mu }[/math] ensembles, the PV term is zero.
[math]\displaystyle{ V }[/math] Volume (traditional) CAUTION: Volume in this sense is the kind you think of with box size. Not to be confused with the Phase Space Volume. This particular volume is not often considered by itself and should only show up on this site within [math]\displaystyle{ PV }[/math] terms, or NVT.
[math]\displaystyle{ \Gamma }[/math] Phase Space Volume: Total set of coordinates and momenta where the system has a nonzero probability of being found CAUTION: This is not volume in standard physical sense, meaning the volume of the box, but the volume space times the momentum space.
[math]\displaystyle{ q }[/math] or [math]\displaystyle{ \vec{q} }[/math] Collective variable for coordinates and momentum. Can also be seen as [math]\displaystyle{ \mathbf{x} }[/math], [math]\displaystyle{ \mathbf{\vec{x}} }[/math], and a number of other terms in literature.
[math]\displaystyle{ \epsilon }[/math] and [math]\displaystyle{ \sigma }[/math] Lennard-Jones parameters for the well depth and the radius of the point particle respectively.
[math]\displaystyle{ Q }[/math] Partition function for the canonical or NVT ensemble, but sometimes also signifies an arbitrary partition function of any ensemble.
[math]\displaystyle{ Q=\int_{\Gamma}e^{-\beta U (\vec{q})}\,\mathrm{d}\vec{q} }[/math]
In the examples here, this will represent the canonical partition function although most of the theory shown here works for an arbitrary partition function as well, though sometimes with slight modifications, where [math]\displaystyle{ Q=\int_{\Gamma}e^{-u(\vec{q})}\,\mathrm{d}\vec{q} }[/math]
[math]\displaystyle{ \tau }[/math] Autocorrelation time
NVT System held at constant Number of particles, Volume, and Temperature. Free energy calculations performed in this ensemble yield Helmholtz free energies.
NPT System held at constant Number of particles, Pressue, and Temperature. Free energy calculations performed in this ensemble yield Gibbs free energies.
TI Thermodynamic Integration
EXP Exponential Averaging Also called "free energy perturbation" and the "Zwanzig relationship."
BAR Bennett Acceptance Ratio
MBAR Multistate Bennett Acceptance Ratio
Boltzmann Weight The probability of a microstate [math]\displaystyle{ \vec{q} }[/math], defined [math]\displaystyle{ e^{-u(\vec{q})}/Q }[/math], where [math]\displaystyle{ u(\vec{q}) }[/math] is the reduced potential. It can also refer to simply [math]\displaystyle{ e^{-u(\vec{q})} }[/math], although this is not a full weight because it is not normalized by the partition function.
State (variable [math]\displaystyle{ K }[/math]) A unique set of conditions and parameters that completely describe a thermodynamic system. For all examples here, the upper case variable [math]\displaystyle{ K }[/math] will denote all states (or the count thereof), and the lower case variable [math]\displaystyle{ k }[/math] will refer to a specific state.
[math]\displaystyle{ \lambda }[/math] Alchemical variable, thermodynamic path coordinate The most common symbol for denoting/tracking the progress along some alchemical pathway. Since most free energy calculations involve an alchemical thermodynamic path, this site was named accordingly.
Annihilating/Annihilated The act of removing ALL of an atoms interactions with the system, this includes the molecule's interactions with itself. Intermolecular force = OFF, Intramolecular forces = OFF. Also known as disappeared or destroyed. Not to be confused with "decoupling."
Decoupling The act of removing an atom or molecules interactions with its surroundings ONLY; the interactions between the molecule itself remain active Intermolecular force = OFF, Intramolecular forces = ON. Not to be confused with "annihilation."

Facts from the Free Energy Difference Definition

There are a number of key notes we can learn from the definition of free energy differences. Each of these can be important in interpreting simulation results.

  1. Only free energy differences are ever calculated. There is never a calculation where absolute free energies are needed (and rarely can they be calculated at all) as all of the biological or thermodynamical quantities of interest are based on a free energy difference. As such, there must always be a minimum of two defined thermodynamic states.
    • Even absolute free energies of binding are still free energy differences between two states: the ligand restricted to the binding site, and the ligand free to explore all other configurations.
  2. Free energy differences between states at different temperatures are usually not what you want to be calculating for problems of interest. If it did, you would get [math]\displaystyle{ \Delta A_{ij} = -k_B T_j \ln Q_j + k_B T_i \ln Q_i }[/math], which is no longer a ratio calculation and not needed for biological systems of interest. Temperature dependence on free energy is more likely to be "what is [math]\displaystyle{ \Delta A_{ij} }[/math] change at two different temperatures?"
  3. There are two different phase space volumes. [math]\displaystyle{ \Gamma_i }[/math] and [math]\displaystyle{ \Gamma_j }[/math] are often the same, but they are not required to be. The methods presented here almost always assumes that the two phase spaces overlap. However, when the spaces do not overlap, these methods break down and it is difficult to identify this problem without in depth knowledge of your system. Consider the example of a hard sphere solute with radius [math]\displaystyle{ \sigma }[/math] at state [math]\displaystyle{ i }[/math] and a Lennard Jones repulsion/dispersion potential, with the same [math]\displaystyle{ \sigma }[/math] at state [math]\displaystyle{ j }[/math]. Since [math]\displaystyle{ \Gamma _i }[/math] will not have molecules at a distance less than [math]\displaystyle{ \sigma }[/math], but [math]\displaystyle{ \Gamma_j }[/math] will, the two phase spaces are not the same and these methods will either break down or return very error-prone results.
    • The degree to which the phase spaces are shared is called the "phase space overlap". Efficient free energy calculations require significant phase space overlap. There are a number of strategies to address lack of overlap between target spaces, but determining the best way for any given situation is still a research question.
    • It should also be noted that "near zero probability" and "always zero probability" are two distinct things when considering phase space. So long as there is a chance for an observation to be made, no matter how small, it is considered part of the phase space.

References

  1. >Woods, C. J., Manby, F. R., and Mulholland, A. J. (2008) An efficient method for the calculation of quantum mechanics/molecular mechanics free energies. J. Chem. Phys. 128, 014109. - Find at Cite-U-Like