Free Energy Calculations

Monte Carlo simulations can be used to create large numbers of configurations of a molecular system, with each configuration appearing with its correct (relative) Boltzmann probability. This collection of configurations is called an ensemble. You can do many things with a correct Boltzmann weighted ensemble of configurations. One possibility is to use it to calculate a free energy.

Change to the protoms_tutorial/ethane_methanol/free_energy directory. In here, you will see the ProtoMS input command file input.cmd. This script performs a simulation that calculates the relative hydration free energy of ethane and methanol. The simulation works by simulating both ethane and methanol simultaneously in the centre of a box of water. A lambda coordinate is used to turn off the interactions between ethane and water as the interactions between methanol and water are switched on. This is a "dual topology" method of calculating a relative free energy. The input file simulates the system at lambda=0.5. This means that 50% of the methanol and 50% of the ethane interactions with solvent are calculated. The energy difference for each configuration and lambda=0 and lambda=1 is calculated. This energy difference is averaged using the equation;

dG = -kT ln < exp( -dE / kT ) >_lambda

where dE is the energy difference between lambda=1 or lambda=0 and lambda=0.5. kT is Boltzmann's constant multiplied by temperature. The angle brackets < ... > mean "average over configurations in an ensemble", with _lambda meaning the ensemble generated for a specific value of lambda. The resulting value, dG is the difference in free energy between lambda=1 and lambda=0.5, or lambda=0 and lambda=0.5. The sum of these two free energy differences is the relative hydration free energy.

Run the free energy simulation by typing;

$ /path/protoms2 input.cmd

The final lines of the output contain the results. Look for the lines;

RESULTS Backwards Free Energy =      1.9251908496670806503     ( 2.013836560     )
RESULTS Forwards  Free Energy =     0.72142800990255362414     ( 1.045066552     )

These give the free energy average (and error in brackets) for the 0.5->0.0 perturbation and the 0.5->1.0 perturbation. As you can see, the error is very large, because the simulation is too short.

---

---

Question
Increase the number of steps of the simulation and see if you can reduce the size of the errors. What is the predicted relative hydration free energy of ethane and methanol from your simulation? How does this compare to experiment? What sources of error are there in the simulation?

---

Question
Another way to reduce the size of the errors is to reduce the gap between the reference state (lambda=0.5) and the perturbed states (lambda=0 and lambda=1). Run a series of simulations that break the relative free energy into a series of steps, e.g. 0.1->0.0 and 0.1->0.2, then 0.3->0.2 and 0.3->0.4, then ... 0.9->0.8 and 0.9->1.0.

Run the simulations and see whether the resulting relative hydration free energy agrees more closely with experiment.

---