[gmx-users] Free energy and dhdl calculation in Gromacs

[Aho, Noora]

Hi all,

I am wondering, how is the dh/dl calculated when doing free energy runs in Gromacs.

I performed a free energy calculation for ASP deprotonation, from lambda=0 (protonated) to lambda=1.0 (deprotonated) with points 0.0, 0.1, 0.2,…,0.9,1.0, to get dhdl for each lambda. Then I calculated dhdl values “by hand". The following were the steps:

1) run a short simulation (only charges changed) for one lambda point -> get dhdl
2) perform -rerun using a topology with the a-state charges
3) perform -rerun using a topology with the b-state charges
4) for both reruns, get coulomb energies (or total potential energies) and calculate V_b - V_a

For each lambda, I get different values for dhdl and V_b - V_a from the reruns. But by definition, is this not exactly how dhdl is calculated (linear interpolation, then dhdl=Vb-Va)?  My results were, for example for lambda=0.7, dhdl=-430.9 kJ/mol and V_b - V_a=-408.1 kJ/mol.

Also, for lambda=0 and lambda=1, I get huge (~10^3) dhdl values from the free energy runs, however V_b - V_a from reruns gives reasonable values. I know that for lambda=1 this might be because no soft-core potentials were used and the charge there disappears, but for lambda=0 I don’t see the reason.

Lastly, if I change the cut-off scheme from Group to Verlet, the dhdl values (and V_b - V_a) are really different. For example for lambda=0.7, with Verlet I get dhdl=-513.3 kJ/mol, V_b - V_a=-485.7 kJ/mol, and with Group the above mentioned dhdl=-430.9 kJ/mol and V_b - V_a=-408.1 kJ/mol. 


So altogether my questions are:
Question 1: how is dh/dl calculated for a lambda point in Gromacs?
Question 2: why dhdl can be huge for lambdas 0,1 even though rerun with V_b - V_a gives a correct result?
Question 3: why is the dhdl value so different for Group and Verlet cut-off schemes?

Best regards,
Noora Aho