Thermodynamic Integration

Free Energy Fundamentals
Theory Definitions Exponential Averaging / Zwanzig Thermodynamic Integration Bennett Acceptance Ratio WHAM MBAR
Methods of Free Energy Simulations Thermodynamic Cycle Intermediate States (Constructing) Soft Core Potentials Running Simulations Extracting Simulation Information Analysis / Bootstrap Sampling
Improving Sampling and Accelerating Simulations Umbrella Sampling Hamiltonian Replica Exchange [math]\displaystyle{ \lambda }[/math]-dynamics / Expanded Ensemble

Free Energy How-to's
Example Calculations Absolute Solvation Free Energy Relative Binding Affinity Absolute Binding Affinity
Tutorials Ethanol Solvation - GROMACS Expanded Ensemble - GROMACS Absolute Binding Free Energy - GROMACS

Thermodynamic Integration (TI) is one of the most common methods for calculating free energy differences, and the easiest to understand. The basic relationship can be calculated by taking the derivative of the free energy difference with respect to [math]\displaystyle{ \lambda }[/math]. It is one of the few methods that require calculation of [math]\displaystyle{ \frac{\partial U(\lambda,\vec{q})}{\partial\lambda} }[/math]. The need to calculate the derivative is also one of its limitations as many simulation packages will not calculate this natively, and can cause problems when numerically evaluating it at [math]\displaystyle{ r = 0 }[/math]. Despite this, it is still one of the most accurate methods if used correctly and it is generally recommended that those new to the field because of its simplicity and ease of use if available in your code if interest.

Derivation

Starting with the identity of free energy

[math]\displaystyle{ \displaystyle A = -\beta^{-1} \ln Q }[/math]

taking the derivative with respect to [math]\displaystyle{ \lambda }[/math] yields

[math]\displaystyle{ \displaystyle dA/d\lambda = -\beta^{-1}\frac{d}{d\lambda} \ln \int e^{-\beta U(\lambda,\vec{q})}d\vec{q} = -\beta^{-1} \frac{\frac{d}{d\lambda}\int e^{-\beta U(\lambda,\vec{q})}d\vec{q}}{Q} }[/math]

which can then be written as

[math]\displaystyle{ \displaystyle dA/d\lambda = -\beta^{-1}\frac{-\beta\int \frac{dU(\lambda,\vec{q})}{d\lambda} e^{-\beta U(\lambda,\vec{q})}d\vec{q}}{Q} = \left\langle \frac{dU(\lambda,\vec{q})}{d\lambda}\right\rangle_\lambda }[/math].

Finally, one can do integration over the whole range of [math]\displaystyle{ \lambda }[/math] to get the final TI equation

[math]\displaystyle{ \displaystyle \Delta A = \int_0^1 \left\langle \frac{dU(\lambda,\vec{q})}{d\lambda}\right\rangle_\lambda d\lambda }[/math].

Estimating Free Energies with TI

The above derivation makes it rather simple to estimate free energies from TI as there is no iterative solution needed, and information from only a single state is needed to calculate the derivative. Since we can only perform simulations at a number of [math]\displaystyle{ \lambda }[/math] states, numeric integration schemes are required.

All numeric integration schemes have the form

[math]\displaystyle{ \displaystyle \Delta A \approx \sum_{k=1}^K w_k \left\langle \frac{dU(\lambda,\vec{q})}{d\lambda}\right\rangle_k }[/math]

where the weights, [math]\displaystyle{ w_k }[/math] will depend on which numeric integration style is chosen. We recommend starting users in free energy calculations use the trapezoid rule as it is very straightforward and maximizes flexibility in choice of [math]\displaystyle{ \lambda }[/math] spacing. Under the trapezoid rule, even lambda spacing weights are [math]\displaystyle{ \displaystyle w_1 = w_k = 1/[2(K-1)] }[/math] and [math]\displaystyle{ w_{k \ne 1,K} = 1/(K-1) }[/math]. Other weighting schemes have also been tried but are not recommended for beginners as they each have their own benefits and disadvantages are are often system specific. The worse case scenario is that you have to run a few extra intermediate states to get accurate results.

Calculating the Statistical Variance of TI

The variance of TI is one of the simplest, but has one catch experimenters should be aware of. Because the information from calculating [math]\displaystyle{ \left\langle \frac{dU}{d\lambda} \right\rangle }[/math] requires information from neighboring states, the variances from each interval, which each take information from 2 states, do not add independently.

Let's look at this will all the mathematical details. Since each of the averages is generated from different simulations, the variance for TI is

[math]\displaystyle{ \mathrm{var}\left(\Delta A\right) = \sum_{k=1}^{K}w_k^2 \mathrm{var}\left(\frac{dU}{d\lambda}\right)_k }[/math].

To illustrate how this is different from independent adding, consider the trapezoid rule example. Taking into account the correct equation the variance would be

[math]\displaystyle{ \mathrm{var}\left(\Delta A_{1,K}\right) = \frac{1}{4}\mathrm{var}\left(\frac{dU}{d\lambda}\right)_1 + \mathrm{var}\left(\frac{dU}{d\lambda}\right)_2 + \cdots + \mathrm{var}\left(\frac{dU}{d\lambda}\right)_{K-1} +\frac{1}{4}\mathrm{var}\left(\frac{dU}{d\lambda}\right)_K }[/math].

Compare this to the incorrect method shown below.

[math]\displaystyle{ \mathrm{var}\left(\Delta A_{i,i+1}\right) = \frac{1}{4}\mathrm{var}\left(\frac{dU}{d\lambda}\right)_i + \frac{1}{4}\mathrm{var}\left(\frac{dU}{d\lambda}\right)_{i+1} }[/math]

[math]\displaystyle{ \begin{alignat}{2} \mathrm{var}\left(\Delta A_{1,K}\right) &=\sum_{i=1}^{K-1}\mathrm{var}\left(\Delta A_{i,i+1}\right) \\ &= \frac{1}{4}\mathrm{var}\left(\frac{dU}{d\lambda}\right)_1 + \frac{1}{2}\mathrm{var}\left(\frac{dU}{d\lambda}\right)_2 + \cdots + \frac{1}{2}\mathrm{var}\left(\frac{dU}{d\lambda}\right)_{K-1} +\frac{1}{4}\mathrm{var}\left(\frac{dU}{d\lambda}\right)_K \\ \end{alignat} }[/math]

The second set of equations is clearly quite different from the correct first set. To get the error, each individual average should be multiplied by [math]\displaystyle{ \sqrt{2\tau} }[/math] to correct for the correlation time at each state; error is then [math]\displaystyle{ \sqrt{\mathrm{var}\left(\Delta A_{i,K}\right)} }[/math].

Problems with TI

Although TI is one of the simplest methods, it also suffers from some of the more difficult challenges that must be overcome. For instance, if the curvature of [math]\displaystyle{ \frac{dU}{d\lambda} }[/math] is large, the bias introduced by discrete [math]\displaystyle{ \lambda }[/math] states becomes significant. This said, when using TI, it is of the utmost importance that researchers verify that sufficient states are present such that the free energy is independent of number of states.

Infinite [math]\displaystyle{ dU/d\lambda }[/math]

One of the largest problems with TI is evaluation of [math]\displaystyle{ \frac{dU}{d\lambda} }[/math] at [math]\displaystyle{ r = 0 }[/math] when standard Lennard-Jones potentials are used. The simplest thermodynamic pathway one can chose is the linear pathway of

[math]\displaystyle{ U(\lambda,\vec{q}) = (1-\lambda)U_0(\vec{q}) + \lambda U_1(\vec{q}) + U_{unaffected}(\vec{q}) }[/math]

which is acceptable for changes in parameters, but not for annihilating or decoupling a site because of the [math]\displaystyle{ r^{-12} }[/math] term in the Lennard-Jones potentials in [math]\displaystyle{ U_i }[/math]. The linear transformation will always have an infinite potential at [math]\displaystyle{ r=0 }[/math] leading to numeric instabilities for evaluating [math]\displaystyle{ \frac{dU}{d\lambda} }[/math] in TI. Although there are ways to get around this, they do not converge very will with any function of [math]\displaystyle{ \lambda }[/math] taking the form shown above. However, if one were to use a soft core potential, this problem can be mostly avoided.

Modifying Simulations

Because most other free energy methods do not need to evaluate [math]\displaystyle{ \frac{dU}{d\lambda} }[/math], they do not natively support evaluating this in the code. If the thermodynamic path is constructed with a linear transformation, then the derivative can be evaluated in post-processing knowing the energy of the system. However, if the transformation is done with soft core potentials, then the derivative will need to be evaluated in code, and it will often be necessary for researchers to modify the code in-house as many simulation packages do not evaluate this quantity.