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.

Relative hydration free energy of p-cresol with respect to toluene


In this exercise, we will study the mutation of toluene into p-cresol in order to calculate their relative hydration free energy. For that we will use the following thermodynamic cycle:

thermodynamic cycle

ΔGhyd(T): hydration free energy of the toluene.
ΔGhyd(C): hydration free energy of the cresol.
ΔG1: free energy associated to the mutation of toluene into cresol in vacuo.
ΔG2: same as ΔG1 in water.

From the cycle we get the following equation:
ΔGhyd(C) - ΔGhyd(T) = - ΔG1 + ΔG2

The procedure is very similar to the one used for the hydration free energy of toluene since you still have to perform the mutation both in water and vacuo using the same protocol (first equilibration and then free energy calculation with equilibration and data collection for each lambda point).

The only difference comes from the fact that only a few atoms will be perturbed this time (and that we will make one atom "appear" instead of making some "disappear").

During this exercise you will:

1. Edit the structure file and the topology file to describe the mutation

2. Perform an equilibration of the system

3. Perform the free energy calculation

4. Analyze the results

First create a directory for this new exercise and then enter it:

mkdir p-cresol

cd p-cresol

Create the directories for the simulations in water and in vacuo (mkdir water vacuo) and copy the files you need from the previous exercise:

cp ../toluene/water/EQ/confout.gro   water/conf_ini.gro

cp ../toluene/water/   water/

cp ../toluene/vacuo/EQ/confout.gro   vacuo/conf_ini.gro

Also this time we will focus on the simulations in water but the steps are exactly the same in vacuo.
Go to the water directory (cd water)

Step 1. Edit the topology file and the structure file.

Make the following changes in the topology file

- Under the directive [ atoms ], change lines corresponding to the atoms 11 12 13 into the following:

    11      C      1    PHE     CZ      6       -0.1     12.011     C     0.15        12.011
    12     HC      1    PHE     HZ      6        0.1      1.008    OA   -0.548         1.008
    13    DUM      1    PHE    DUM      6          0      1.008     H    0.398         1.008

- Under the directive [ bonds ], change corresponding lines into the following:

    7    11     2    gb_15   gb_15
    9    11     2    gb_15   gb_15
   11    12     2    gb_3    gb_12
   12    13     2    gb_1    gb_1

- Under the directive [ pairs ], add the following:

    7    13     1
    9    13     1 

- Under the directive [ angles ], change corresponding lines intothe following:

    3     7    11     2    ga_26   ga_26
    8     7    11     2    ga_24   ga_24
    5     9    11     2    ga_26   ga_26
   10     9    11     2    ga_24   ga_24
    7    11     9     2    ga_26   ga_26
    7    11    12     2    ga_24   ga_26
    9    11    12     2    ga_24   ga_26
   11    12    13     2    ga_11   ga_11

- Under the first directive [ dihedrals ], add the following:

    7    11    12    13     1    gd_2    gd_2

- Under the second directive [ dihedrals ], change corresponding lines into the following:

    2     3     7    11     2    gi_1    gi_1
    2     5     9    11     2    gi_1    gi_1
    3     7    11     9     2    gi_1    gi_1
    5     9    11     7     2    gi_1    gi_1
    7    11     3     8     2    gi_1    gi_1
    9    11     5    10     2    gi_1    gi_1
   11     7     9    12     2    gi_1    gi_1

You can use this modified topology to perform the simulation in vacuo but like in the previous exercise, do not forget to remove the last line relative to to water molecules.

Structure file conf_ini.gro:

On the second line there is the number of atoms that you should increase by 1. You also have to add the dummy atom that will be turned into a H atom and give it coordinates very close to the ones of the previous atom:

    1PHE     HZ   12   1.243   1.843   1.171  0.0722 -0.2555  1.1745
    1PHE    DUM   13   1.240   1.840   1.170  0.0000  0.0000  0.0000

Do not worry about the numbering of the atoms as Gromacs does not use this information but read the atoms sequentially instead (only the order is important).

Simulation in water

Step 2. Equilibration of the system.

The force field will take care of positioning the dummy atom correctly during this stage while the temperature coupling will give the correct velocity distribution for all the atoms.

Create a new directory EQ (mkdir EQ), enter it (cd EQ) and copy the parameter file equilibration.mdp.

Check the contents of the file:

less equilibration.mdp

You can see that there is no new velocity generation as we use the set of velocity already present in the structure file.

Create the run input file (topol.tpr) and start the simulation:

grompp  -f  equilibration.mdp  -c  ../conf_ini.gro  -p  ../

mdrun  -s  topol.tpr

Step 3. Free energy calculation.

The procedure is exactly the same as in the Step 6 of the previous exercise.

Parameter files:

Step 4. Analysis.

The procedure is exactly the same as in the Step 7 of the previous exercise.

The result in this case should look like the file dHdl.xvg.
What value do you get for the integral ?

You still need to make the calculation in vacuo to be able to compare the results with the experimental value.

Simulation in vacuo.
For the simulation in vacuo, you can use the topology created for the simulations in water but you have to delete the last line relative to the solvent. For the structure file, proceed as indicated in Step 1 and for the simulations follow the procedure starting from Step 2 with the following parameter files:

Integrate the plot obtained by gathering the results for the different λ values in order to get the mutation free energy in vacuo.


With the results of the simulations performed both in water and vacuo, you can calculate the relative hydration free energy of p-cresol with respect to toluene.

ΔΔGhydT -> C = ΔGhydp-cresol -  ΔGhydtoluene =  ΔGmutationT -> C(water)  -  ΔGmutationT -> C(vacuo).

Compare it to the experimental value:  ΔGhydp-cresol  -  ΔGhydtoluene =  - 22.1 kJ.mol-1

<<- Back

This page was last modified on Jan. 19th 2004