Free Energy Calculation

User discussions
1

GROMACS version: 5.1.4
GROMACS modification: No
Hi all,

I have run a free-energy perturbation (FEP) with the slow-growth approach. So, my lambda increases at each MD step till reaching the value of 1 at the last step of my MD.
So, for each step, I get force on a given lambda values.
Do you know how to integrate it along the lambda growth to get the free energy? gmx bar does not be able to handle such a structured xvg file. I think it expects a constant lambda value and different xvg files instead.
Considering I am also running reverse alchemical transformation.

Thanks everybody for your contribution.

Cheers

VG

2

Hi vg87,

Slightly indirect answer - if you can do slow-growth in both ways, you can also do fast-growth and that will be more robust in almost all cases. For a neat set of tools to analyse that type of data see

http://pmx.mpibpc.mpg.de/

You can even use some of the analysis scripts there to integrate your slow-growth data.

1 Like
3

Thanks a lot!

I will make a try!

V

4

Actually what I have done so far:

I am transforming one of the non-bridging oxigen of a phosphate group in a TA dinucleotide into an S.
So I am passing from Rp to Sp stereochemistry and asses the difference in terms of energetics.

Here the ‘free-energy’ setting I have intensely discussed with Justin Lemkul along these days:

== SETTING ==

; Free energy control stuff
free_energy = yes
init_lambda_state = 0
delta_lambda = 0
calc_lambda_neighbors = 1 ; only immediate neighboring windows
couple-intramol = yes
; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
vdw_lambdas = 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
coul_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
bonded_lambdas = 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
mass_lambdas = 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
restraint_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Options for the decoupling
sc-alpha = 0.5
sc-coul = no ; linear interpolation of Coulomb (none in this case)
sc-power = 1
sc-sigma = 0.3
nstdhdl = 10

while here the results:

point 0 - 1, DG 6.36 +/- 0.04
point 1 - 2, DG 4.99 +/- 0.07
point 2 - 3, DG 3.48 +/- 0.04
point 3 - 4, DG 2.72 +/- 0.05
point 4 - 5, DG 1.87 +/- 0.04
point 5 - 6, DG 2.33 +/- 0.06
point 6 - 7, DG 1.67 +/- 0.06
point 7 - 8, DG 0.86 +/- 0.03
point 8 - 9, DG 0.09 +/- 0.08
point 9 - 10, DG -0.81 +/- 0.03
point 10 - 11, DG 1.82 +/- 0.01
point 11 - 12, DG 1.54 +/- 0.02
point 12 - 13, DG 1.29 +/- 0.01
point 13 - 14, DG 1.04 +/- 0.01
point 14 - 15, DG 0.75 +/- 0.01
point 15 - 16, DG 0.40 +/- 0.07
point 16 - 17, DG 0.10 +/- 0.08
point 17 - 18, DG -0.11 +/- 0.01
point 18 - 19, DG -0.38 +/- 0.01
point 19 - 20, DG -0.70 +/- 0.03

total 0 - 20, DG 29.30 +/- 0.30

Now, experimentally I know that passing from a phosphate group (in TA) to a phosphorothioate group (of course in TA) has a dG = +2.5 kcal mol

As you may note, values are really really far.
Each window is 2 ns long, but even with 5 ns thing do not change.

Have you got any idea guys?

Thanks in advance for any kind contribution.

Regards

VG