Exercise 2 - Relative hydration free energy


You will learn to compute free energy changes in the gas phase using several free energy techniques. You will learn the differences between single topology and a dual topology calculations. You will learn to compute free energy changes in a water box. You will learn to analyse the results from the simulations.

Single topology model of ethane to methanol

Go to

cd $HOME/workshop/exercise2/

and setup the example files

cp in/* .

We will start with setting up a single topology template to transform ethane into methanol. This can easily be done with protoms.py, but first we need to consider what ethane atoms will be perturbed into methanol atoms. The sketch below describes the perturbation

    H02  H06               H02  DUM
     |    |                 |    |                     
H03-C01--C05-H07  ==>  H03-C01--O05-DUM
     |    |                 |    |          
    H04  H08               H04  H08
    Ethane (λ=0.0)     Methanol  (λ=1.0)    

in simple terms, the C05 atom of ethane will be perturbed into the O05 atom of methanol and the H06 and H07 atoms of ethanol will be perturbed into dummy atoms (as marked by DUM). The H08 atom will have the same name in both ethane and methanol, but it will change atom type and charge. In fact, H02, H03 and H04 are also being perturbed to different atom types. 

This is sketch is important to keep in mind as we setup the single-topology calculation with

protoms.py -s singletopology -l ethane.pdb methanol.pdb --testrun

the script will now prompt you to type in the corresponding atoms. To your help you have the sketch above and a distance matrix that the program writes out. When prompted for H06 and H07 you leave the corresponding atom blank to indicate that these should be perturbed to dummy atoms.

You should inspect the created template file called ethtmeo_comb.tem. 

This template looks similar to the template you have seen for the molecule of ethane. However, look at the 'atom' section. For instance

atom C05 ETH 3000 3009 DM3 DUM DM2 DUM DM1 DUM

Notice that the first two numbers '3000' and '3009'. These numbers refer to the atomic parameters used for this solute and are listed in the mode clj section of the template. In the previous example the numbers were the same because the forcefield parameters assigned to an atom do no vary as the coupling parameter lambda changes. Here, in the present single topology simulation, the force field parameters must describe a molecule of ethane at lambda = 0 (the first column of numbers) and a molecule of methanol at lambda = 1 (the second column of numbers).

To illustrate the single topology model, a special input file for ProtoMS is provided. It is called run_single.cmd.  Inspect the contents of the file input.cmd. In this command file, notice that there is a new keyword lambda, initially set at 0.000 0.000 0.000. ProtoMS then uses successive chunks to set different values of the coupling parameter lambda and run a one-step MC simulation..

Execute ProtoMS with this file

protoms3 run_single.cmd

and visualize the the all.pdb file in the folder out_single with VMD to see how the solute progressively morphs from ethane to methanol as lambda increases. This is the essence of single topology.

Relative free energy in the gas phase

In this exercise, you will compute the relative free energy of ethane and methanol, in the gas phase. As seen in the lecture, this is one of the two legs required to compute a relative free energy of hydration.

We will be using a single topology model, with a variety of techniques. The ‘correct’ answer as estimated by long calculations and for the parameters used is +2.8 kcal.mol-1.

We will start with computing the free energy using 3 λ-values. 

protoms.py -s singletopology -l ethane.pdb methanol.pdb --singlemap single_cmap.dat --testrun --lambdas 3

This will create a number of command files for ProtoMS. The one we will be using in this exercise is run_comb_gas.cmd, will will run the single-topology free energy calculation in gas phase. For more complicated perturbations it is common to divide it into an electrostatic and a van der Waals perturbation. Such command files are also produced by protoms.py

This job will requires a few processor and will take a couple of minutes to run. Therefore, we will submit the job to the Lyceum queue by executing

./run_gas.bash 3

(to run ProtoMS with the command file one simple execute mpirun -np 3 protoms3 run_comb_gas.cmd, but this is done in run_gas.bash)

You can check your job in the queue with the command qstat.

When the job is finished you can calculate the free energy using

calc_dg.py -d out_comb_gas/

The program will first print the gradient with respect to λ and then it will give you estimate from thermodynamic integration (TI), Bennet's Acceptance Ratio (BAR) and multi BAR (MBAR). The results will look something like this:

lambda     gradient      std
0.000        1.4922   0.0111
0.500        3.0235   0.0010
1.000        2.7565   0.0009

TI (out_comb_gas/) : 2.574 +- 0.003

BAR (out_comb_gas/): 2.812 +- 0.001

MBAR (out_comb_gas/): 2.812 +- 0.000

You will notice that TI gives different results compared to BAR and MBAR. The problem with TI is that the free energy gradients are integrated numerically. Here we have used the trapezium rule, but with only 3 points the true shape of the free energy gradients vs lambda plot (otherwise called the potential of mean force) is likely grossly approximated. Copy the gradient into an Excel sheet and plot it against λ.

Now setup a calculation using 5 λ-values.

protoms.py -s singletopology -l ethane.pdb methanol.pdb --singlemap single_cmap.dat --testrun --lambdas 5

move the old output folder

mv out_comb_gas out_comb_gas_3lam

and then submit the job to the queue. 

./run_gas.bash 5

wait for it to finish and then analyse it as above. Now the TI results is much closer to the BAR and MBAR results. Compared to the previous situation where the integration was made with 3 points, how much has the shape of the curve changed?

Finally, setup, run and analyse a job with 8 λ-values. 

Relative free energy in a box of water

Now that you know how to compute a free energy change with ProtoMS using the TI, BAR and MBAR techniques on a trivial gas phase system, we are going to tackle the more challenging problem of computing free energy changes in a condensed phase. Specifically we will compute the free energy change for perturbing ethane into methanol in a box of water. This time the calculations will take longer since we will have to simulate many more Monte Carlo moves to allow the relaxation of the solvent as the solute is perturbed. Two alternative techniques, single topology and dual topology, will be covered. These calculations will run over lunch.

The experimental relative free energy of hydration between ethane and methanol is -6.90 kcal.mol-1.

Single-topology setup

The water simulations will be run with 100 k equilibration steps and 10 m sampling steps. 8 λ-values will be used. This can be setup using

protoms.py -s singletopology -l ethane.pdb methanol.pdb --singlemap single_cmap.dat --lambdas 8 --nequil 1E5 --nprod 10E6

Look at the command file run_comb_free.cmd. There are number of new keywords compared with the previous gas phase example.

solvent1 ethane_box.pdb
boundary solvent
pressure 1
chunk simulate 10000000 solvent=981 protein=0 solute=16 volume=3

These additional keywords are used because we need to load the coordinates of box of water molecules (solvent1), instruct ProtoMS to use periodic boundary conditions (boundary solvent), set an NPT simulation and therefore specify the pressure.

The last line instructs ProtoMS that we will sample solute degrees of freedom about 1.6% of the time (16 moves out of 1000), solvent degrees of freedom about 98% of the time (981 moves out of 1000), and volume box moves 0.3% of the time (3 moves out of 1000).

Also look at the top of ethane_box.pdb

HEADER box -12.3724 -11.6015 -12.2139 11.7076 12.4785 11.8661

the first line contains a HEADER file which lists the box dimensions. This box was created from a pre-equilibrated box of TIP4P water molecules and extends at least 10 A from the ethane molecule. Copy ethane.pdb and ethane_box.pdb to your workstation to visualize them and make sure you understand what you are going to simulate.

Once you feel ready, submit the job by typing


this job will take about an hour to complete. Notice that this is a rather short run, the default in protoms.py is to run 5 m equilibration moves and 40 m sampling moves.

Dual-topology setup

Now we will setup the dual topology simulations. The command for setting up this type of simulation is very similar to what you already have done. Simply type

protoms.py -s dualtopology -l ethane.pdb methanol.pdb --lambdas 8 --nequil 1E5 --nprod 10E6

Look at the command file run_free.cmd. There are number of new keywords compared with the previous single topology example.

parfile eth-meo.tem
dualtopology1 1 2 synctrans syncrot

softcore1 solute 1
softcore2 solute 2
softcoreparams coul 1 delta 0.2 deltacoul 2.0 power 6 soft66

The template file that is used eth-meo.tem contains templates for both ethane and methanol.

The dualtopology keyword is instructing ProtoMS to pair solute1 with solute2. The synctrans and syncrot keywords are used to synchronize the translations and rotations of the two solutes.

Copy back to the workstation the files ethane_box.pdb, ethane.pdb and methanol.pdb. Visualize the system. Notice how the two solutes are overlapping. The dualtopology keyword caused ProtoMS to ignore intermolecular energies between the two solutes. In the single topology calculations, as lambda varies, the forcefield parameters and geometry of the solute are gradually changed to match methanol. In this dual topology implementation, as lambda increases from 0 to 1, the intermolecular energy of ethane is turned off whereas the intermolecular energy of methanol is turned on.

The softcore keywords are instructing ProtoMS to use a softened intermolecular energy function for the solutes, and the precise functional form can be tuned with the softcoreparams line.

Once you feel ready, submit the job by typing


this job will take about an hour to complete. Notice that this is a rather short run, the default in protoms.py is to run 5 m moves of equilibration and 40 m moves of sampling.

Analysis of the results

The results for the single-topology simulations are located in the folder


and dual-topology simulations are located in the folder


Because these runs are rather short, you are also provided with results files from more extensive simulations in the long_sim folder. This folder thus also contains out_comb_free and out_free.

We will start to look at the convergence of the long simulations. The quantity that we are interested in checking for convergence is the total energy of the system. To check the convergence for the simulation at λ=0 of the single topology simulation, type the following

calc_series.py -f long_sim/out_comb_free/lam-0.000/results -s total -o total_lam-0.000_single

this will produce a PNG file called total_lam-0.000_single.png that you can download and display. The plot will look something like this.

 λ=0, single topology

  λ=1, single topology

The program will also report that point at which the total energy is equilibrated and show it in the plot has a vertical line. This should be around 100 snapshots into the simulation.

Now, perform the same analysis for λ=1. Then you need also need to analyse the dual-topology simulation. Is there a difference in equilibration time? You will probably find that the program the program report a much longer equilibration time. However, looking at the plot we could argue that by visual inspection we don't think this is realistic. It is important to realize that the program gives suggestions that you should to validate by yourself.

To calculate the single-topology free energy type.

calc_dg.py -d long_sim/out_comb_free -s 100

where we remove 100 snapshots from the beginning of the simulation. You should obtain something like this

lambda gradient std
0.000 -3.9473 0.2747
0.143 -0.3176 0.1700
0.286 0.7885 0.1150
0.429 -0.5429 0.1192
0.571 -2.4938 0.1398
0.714 -5.0220 0.1403
0.857 -7.7318 0.1715
1.000 -10.4733 0.1736

TI (out_comb_free/) : -3.220 +- 0.056

BAR (out_comb_free/): -3.198 +- 0.031

MBAR (out_comb_free/): -3.241 +- 0.073

Now repeat the calculation for the dual-topology simulations.

Once you have a free energy change and error estimate for each calculation, draw thermodynamic cycles for the perturbation of ethane in methanol in gas phase and water. You should be able to compute the relative hydration free energy using the results from the single and dual topology calculations.
Compare the results you obtained with single and dual topology. Why do you not need to compute the free energy change in the gas phase with the dual topology calculation?