Rigorous alchemical methods such as Thermodynamic Integration (TI，热力学积分) and Free Energy Perturbation (FEP，自由能微扰) have shown considerable success in making accurate predictions.

# Free energy perturbation

Free energy perturbation (FEP) is a method based on statistical mechanics that is used in computational chemistry for computing free energy differences from molecular dynamics or Metropolis Monte Carlo simulations. The FEP method was introduced by Robert W. Zwanzig in 1954.[1] According to the free-energy perturbation method, the free energy difference for going from state A to state B is obtained from the following equation, known as the Zwanzig equation:

where T is the temperature, kB is Boltzmann’s constant, and the triangular brackets denote an average over a simulation run for state A. In practice, one runs a normal simulation for state A, but each time a new configuration is accepted, the energy for state B is also computed. The difference between states A and B may be in the atom types involved, in which case the ΔF obtained is for “mutating” one molecule onto another, or it may be a difference of geometry, in which case one obtains a free energy map along one or more reaction coordinates. This free energy map is also known as a potential of mean force or PMF. Free energy perturbation calculations only converge properly when the difference between the two states is small enough; therefore it is usually necessary to divide a perturbation into a series of smaller “windows”, which are computed independently. Since there is no need for constant communication between the simulation for one window and the next, the process can be trivially parallelized by running each window in a different CPU, in what is known as an “embarrassingly parallel” setup.

FEP calculations have been used for studying host-guest binding energetics, pKa predictions, solvent effects on reactions, and enzymatic reactions. For the study of reactions it is often necessary to involve a quantum-mechanical representation of the reaction center because the molecular mechanics force fields used for FEP simulations can’t handle breaking bonds. A hybrid method that has the advantages of both QM and MM calculations is called QM/MM.

Umbrella sampling is another free-energy calculation technique that is typically used for calculating the free-energy change associated with a change in “position” coordinates as opposed to “chemical” coordinates, although Umbrella sampling can also be used for a chemical transformation when the “chemical” coordinate is treated as a dynamic variable (as in the case of the Lambda dynamics approach of Kong and Brooks). An alternative to free energy perturbation for computing potentials of mean force in chemical space is thermodynamic integration. Another alternative, which is probably more efficient, is the Bennett acceptance ratio method.

# Thermodynamic integration

Thermodynamic integration is a method used to compare the difference in free energy between two given states (e.g., A and B) whose potential energies $U_{A}$ and $U_{B}$ have different dependences on the spatial coordinates. Because the free energy of a system is not simply a function of the phase space coordinates of the system, but is instead a function of the Boltzmann-weighted integral over phase space (i.e. partition function), the free energy difference between two states cannot be calculated directly. In thermodynamic integration, the free energy difference is calculated by defining a thermodynamic path between the states and integrating over ensemble-averaged enthalpy changes along the path. Such paths can either be real chemical processes or alchemical processes. An example alchemical process is the Kirkwood’s coupling parameter method.[1]

## Derivation

Consider two systems, A and B, with potential energies $U_{A}$ and $U_{B}$. The potential energy in either system can be calculated as an ensemble average over configurations sampled from a molecular dynamics or Monte Carlo simulation with proper Boltzmann weighting. Now consider a new potential energy function defined as:

Here, $\lambda$ is defined as a coupling parameter with a value between 0 and 1, and thus the potential energy as a function of $\lambda$ varies from the energy of system A for $\lambda =0$ and system B for $\lambda =1$. In the canonical ensemble, the partition function of the system can be written as:

In this notation, $U_{s}(\lambda )$ is the potential energy of state $s$ in the ensemble with potential energy function $U(\lambda )$ as defined above. The free energy of this system is defined as:

If we take the derivative of F with respect to λ, we will get that it equals the ensemble average of the derivative of potential energy with respect to λ.

The change in free energy between states A and B can thus be computed from the integral of the ensemble averaged derivatives of potential energy over the coupling parameter $\lambda$ .[2] In practice, this is performed by defining a potential energy function $U(\lambda )$, sampling the ensemble of equilibrium configurations at a series of $\lambda$ values, calculating the ensemble-averaged derivative of $U(\lambda )$ with respect to $\lambda$ at each $\lambda$ value, and finally computing the integral over the ensemble-averaged derivatives.

Umbrella sampling is a related free energy method. It adds a bias to the potential energy. In the limit of an infinite strong bias it is equivalent to thermodynamic integration.[3]

ref: