NOTE: This tutorial was intended to be used with Gromacs version 3.1.4 and may not be compatible with newer versions. The tutorial will be revised in the near future. Note, a brief tutorial on Free Energy calculation with Gromacs 3.3.x can be found at the DillGroup Wiki.
Experimentally, a difference in free energy is determined either from the relative probability of finding the system in a given state, or from the reversible work required to transform the system from one initial state into another. Computationally, both approaches can be used. However, counting how often the system is in one state, for instance how often a complex is formed when simulating a mixture of two compounds, is extremely inefficient. In practice, a perturbation is applied to force the transition from one state to another. Then, statistical procedures are used to calculate the work done on the system by the perturbation.
The free energy is a state function, which means that the free energy difference is only depending on the initial and end state, no matter what path is taken to go from one to the other. As a consequence, you can choose any nonphysical path you can think of to perform your calculations as long as you can relate them through thermodynamic cycles to the physical process you are interested in.
The most widely used methods to calculate free energy differences are:
The relative free energy between two state A and B is expressed as an integral from λA to λB over the ensemble average of the derivative of the Hamiltonian with respect to the coupling parameter λ.
The integration can be performed continuously while slowly changing the coupling parameter λ from λA to λB during the simulation (the so-called slow growth method ). However this approach can be problematic when the system lag behind the changing hamiltonian and never equilibrate appropriately. A more controlled approach is to simulate the system at a number of fixed λ points and to evaluate the integral numerically. This way the convergence of the simulations at each λ point can be checked independently. This is the approach we will use during the practical.
The free energy difference is calculated as an ensemble average over the state A but the equation can equally be written as an ensemble average over the state B:
In finite sampling, the ensemble average only converges to the correct answer if configurations sampled in state A also have a high probability in state B. The end state B must therefore not be too different from the reference state A.
A simplified method was proposed that reduces the required computational time by avoiding the simulation of uninteresting intermediate state. The method estimates relative free energy by extrapolating in one single step from a well chosen reference state.
Free energy calculations can be performed using either Monte-Carlo or Molecular Dynamics simulations, but during this practical only Molecular Dynamics will be used. As always with this type of methods, the extent of the sampling reached during the simulation will be one of the two primarily limiting factors concerning the accuracy of the calculations. The other one being the underlying model, or force field, used to describe the system.
The following exercises show some of the practical applications that can be made of the free energy calculations.
solvation free energy:
relative solvation free energy: