# Free-energy perturbation theory

We begin our treatment of free energy differences by examining the problem of transforming a system from one thermodynamic state to another. Let these states be denoted generically as \( {\cal A} \) and \({\cal B}\). At the microscopic level, these two states are characterized by potential energy functions \(U_{\cal A}({\textbf r}_1,...,{\textbf r}_N)\) and \(U_{\cal B}({\textbf r}_1,...,{\textbf r}_N)\). For example, in a drug-binding study, the state \( {\cal A} \) might correspond to the unbound ligand and enzyme, while \({\cal B}\) would correspond to the bound complex. In this case, the potential \(U_{\cal A}\) would exclude all interactions between the enzyme and the ligand and the enzyme, whereas they would be included in the potential \(U_{\cal B } \).

The Helmholtz free energy difference between the states \( {\cal A}\) and \({\cal B }\) is simply \( A_{\cal A \cal B}= A_{\cal B}-A_{\cal A}\). The two free energies \(A_{\cal A}\) and \(A_{\cal B}\) are given in terms of their respective canonical partition functions \(Q_{\cal A}\) and \(Q_{\cal B}\), respectively by \(A_{\cal A} = -kT\ln Q_{\cal A}\) and \(A_{\cal B} = -kT \ln Q_{\cal B}\), where

\( Q_{\cal A}(N,V,T)\) | \( C_N\int\;d^N{\textbf p}\;d^N{\textbf r}\;exp\left\{-\beta\left[\sum_{i=1}^N{\textbf p_i^2 \over 2m_i} +U_{\cal A}({\textbf r}_1,...,{\textbf r}_N)\right]\right\}\) | ||

\(\underline { {Z_{\cal A}(N,V,T) \over N!\lambda^{3N}} }\) | |||

\( Q_{\cal B}(N,V,T)\) | \( C_N\int\;d^N{\textbf p}\;d^N{\textbf r}\;exp\left\{-\beta\left[\sum_{i=1}^N{\textbf p_i^2 \over 2m_i} +U_{\cal B}({\textbf r}_1,...,{\textbf r}_N)\right]\right\}\) | ||

\(\underline { {Z_{\cal B}(N,V,T) \over N!\lambda^{3N}} }\) | (1) |

The free energy difference is, therefore,

\[A_{\cal A\cal B}= A_{\cal B} - A_{\cal A} = -kT\ln\left ( {Q_{\cal A} \over Q_{\cal A} }\right)= -kT\ln\left({Z_{\cal B} \over Z_{\cal A}}\right)\] | (2) |

where \(Z_{\cal A}\) and \(Z_{\cal B}\) are the configurational partition functions for states \(\cal A\) and \(\cal B \), respectively,

\(Z_{\cal A}\) | \(\underline { \int\;d^N{\textbf r}\;e^{-\beta U_{\cal A}({\textbf r}_1,...,{\textbf r}_N)} }\) | ||

\(Z_{\cal B}\) | \(\underline { \int\;d^N{\textbf r}\;e^{-\beta U_{\cal B}({\textbf r}_1,...,{\textbf r}_N)} }\) | (3) |

The ratio of full partition functions \(Q_{\cal B}/Q_{\cal A}\) reduces to the ratio of configurational partition functions \(Z_{\cal B}/Z_{\cal A}\) because the momentum integrations in the former cancel out of the ratio.

Eqn. (2) is difficult to implement in practice because in any numerical calculation via either molecular dynamics or Monte Carlo, we do not have direct access to the partition function only averages of phase-space functions corresponding to physical observables. However, if we are willing to extend the class of phase-space functions whose averages we seek to functions that do not necessarily correspond to direct observables, then the ratio of configurational partition functions can be manipulated to be in the form of such an average. Consider inserting unity into the expression for \(Z_{\cal B} \) as follows:

\(Z_{\cal B} \) | \(\underline { \int\;d^N{\textbf r}\;e^{-\beta U_{\cal B}({\textbf r}_1,...,{\textbf r}_N)} }\) | ||

\(\underline { \int\;d^N{\textbf r}\;e^{-\beta U_{\cal B}({\textbf r}_1,...,{\textbf r}_N)}e^{-\beta U_{\cal A}({\textbf r}_1,...,{\textbf r}_N)}e^{\beta U_{\cal A}({\textbf r}_1,...,{\textbf r}_N)} }\) | |||

\( \underline { \int\;d^N{\textbf r}\;e^{-\beta U_{\cal A}({\textbf r}_1,...,{\textbf r}_N)}e^{- \beta (U_{\cal B}({\textbf r}_1,...,{\textbf r}_N)-U_{\cal A}({\textbf r}_1,...,{\textbf r}_N))} }\) | (4) |

If we now take the ratio \({Z_{\cal B} \over Z_{\cal A}}\), we find

\(\underline { {Z_{\cal B} \over Z_{\cal A}} }\) | \(\underline { {1 \over Z_{\cal A}}\int\;d^N{\textbf r}\;e^{-\beta U_{\cal A}({\textbf r}_1,...,{\textbf r}_N)}e^{-\beta(U_{\cal B}({\textbf r}_1,...,{\textbf r}_N)-U_{\cal A}({\textbf r}_1,...,{\textbf r}_N))}}\) | ||

\( \left<e^{-\beta (U_{\cal B}({\textbf r}_1,...,{\textbf r}_N)-U_{\cal A}({\textbf r}_1,...,{\textbf r}_N))}\right>_{\cal A}\) | (5) |

where the notation \(\left<\cdots\right>_{\cal A}\) indicates an average taken with respect to the canonical configurational distribution of the state \(\cal A\). Substituting eqn. (5) into eqn. (2), we find

\[A_{\cal A \cal B}= -kT\ln\left<e^{-\beta (U_{\cal B}-U_{\cal A})}\right>_{\cal A}\] | (6) |

Eqn. (6) is known as the *free-energy perturbation* formula; it should be reminiscent of the thermodynamic perturbation formula used to derive the van der Waals equation. Eqn. (6) can be interpreted as follows: We start with microstates \(\{{\textbf r}_1,...,{\textbf r}_N\}\) selected from the canonical ensemble of state \({\cal A} \) and use these to compute \(Z_{\cal B} \) by placing them in the state \( \cal B \) by simply changing the potential energy from \(U_{\cal A} \) to \(U_{\cal B} \). In so doing, we need to ``unbias'' our choice to sample the configurations from the canonical distribution of state \(\cal A\) by removing the weight factor \(exp (- \beta U_{\cal A}) \) from which the microstates are sample and reweighting the states by the factor \(exp (- \beta U_{\cal B}) \) corresponding to state \(\cal B\). This leads to eqn. (5). The difficulty with this approach is that the microstates corresponding to the canonical distribution of state \(\cal A\) may not be states of high probability in the canonical distribution of state \(\cal B\). If this is the case, then the potential enegy difference \(U_{\cal B} - U_{\cal A} \) will be large, he exponential factor \( exp \left [ -\beta (U_{\cal B} - U_{\cal A}) \right ] \) will be negligibly small, and the free energy difference will be very slow to converge in an actual simulation. For this reason, it is clear that the free-energy perturbation formula is only useful for cases in which the two states \(\cal A\) and \(\cal B\) are not that different from each other.

If \(U_{\cal B} \) is not a small perturbation to \(U_{\cal A} \), then the free-energy perturbation formula can still be salvaged by introducing a set of \(M - 2 \) intermediate states with potentials \(U_{\alpha} (r_1, \cdots, r_N ) \), where \(\alpha = 1, \cdots, M \), \(\alpha = 1 \) corresponds to the state \(\cal A\) and \(\alpha = M \) corresponds to the state \(\cal B\). Let \(\Delta U_{\alpha,\alpha+1} = U_{\alpha+1}-U_{\alpha}\). We can now imagine transforming the system from state \(\cal A\) to state \(\cal B\) by passing through these intermediate states and computing the average of \(\Delta U_{\alpha,\alpha+1}\) in the state \(\alpha\). Applying the free-energy perturbation formula to this protocol yields the free-energy difference as

\[A_{\cal A\cal B}= -kT\sum_{\alpha=1}^{M-1}\ln\left<e^{-\beta\Delta U_{\alpha,\alpha+1}}\right>_{\alpha} \] | (7) |

where \(\langle\cdots\rangle_{\alpha}\) means an average taken over the distribution \(exp (-\beta U_a) \). The key to applying eqn. (7) is choosing the intermediate states so as to achieve sufficient overlap between the intermediate states without requiring a large number of them, i.e. choosing the thermodynamic path between states \(\cal A\) and \(\cal B\) effectively.