Project

General

Profile

Bug #2573

Different mdp files describing the same change yield different free energy on the same trajectory

Added by zhiyi wu about 1 year ago. Updated about 1 year ago.

Status:
New
Priority:
High
Assignee:
-
Category:
-
Target version:
-
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

I was trying to calculate the free energy of decharging n-phenylglycinonitrile in a box of water. A topology file is made where the A-state is the original molecule and B-state is the same atom types with all partial charge scaled to 0.
The simulation was run with five windows (fep-lambdas = 0.0 0.25 0.5 0.75 1.0) and yeilds 12.327 +- 0.109, which is the consistent with other studies.
The mdp files were in the fep folder and the result can be reproduced by

for i in {0..4};
do
gmx grompp -f ./fep/prod.0.mdp -c ligand.gro -p ligand.top -o ./fep/prod_$i.tpr -maxwarn 1
gmx mdrun -s ./fep/prod_$i.tpr -rerun ./simulations/$i/prod.trr -dhdl ./fep/dhdl_$i.xvg
done

However, if the fep-lambdas were split to vdw and coul contributions. The result became 8.153 +- 0.030.
The mdp file can be found in coul_vdw folder, where the

fep-lambdas = 0.0 0.25 0.5 0.75 1.0

is changed into

coul-lambdas = 0.0 0.25 0.5 0.75 1.0
vdw-lambdas = 0.0 0.0 0.0 0.0 0.0

The result can be reproduced by

for i in {0..4};
do
gmx grompp -f ./coul_vdw/prod.0.mdp -c ligand.gro -p ligand.top -o ./coul_vdw/prod_$i.tpr -maxwarn 1
gmx mdrun -s ./coul_vdw/prod_$i.tpr -rerun ./simulations/$i/prod.trr -dhdl ./coul_vdw/dhdl_$i.xvg
done

Since the trajectories are exactly the same, I would expect the result wouldn't change due to this change in mdp files.

This incorrect result can also be reproduced by changing coul-lambdas to fep-lambdas.

fep-lambdas = 0.0 0.25 0.5 0.75 1.0
vdw-lambdas = 0.0 0.0 0.0 0.0 0.0

Which yeilded the same answer as is shown in the fep_vdw folder.

I initially suspect that there were some factors other than the ones controlled by coul-lambdas, contributed free energy to the result.
Thus, I test this hypothesis by changing the mdp file as is shown in fep_coul_0, to compute the free energy conributation other than coul.

fep-lambdas = 0.0 0.25 0.5 0.75 1.0
coul-lambdas = 0.0 0.0  0.0 0.0  0.0

Which yeild 0 kcal/mol.
I have also tested the case where coul-lambdas is 1.0
fep-lambdas = 0.0 0.25 0.5 0.75 1.0
coul-lambdas = 1.0 1.0  1.0 1.0  1.0

Which yeild the same 0 kcal/mol.

Since the file is larger than 50 MB, I will share the file through google drive.
https://drive.google.com/file/d/18ELypM_TaINYJzx0kwWizfliSAxfxR7F/view?usp=sharing

History

#1 Updated by zhiyi wu about 1 year ago

The number 8.153 +- 0.030 kcal/mol is very close to the free energy of decoupling the charge, which is where all the interaction between the surrounding and the ligand is scaled to the zero but intra ligand coul interactions were still present.

#2 Updated by zhiyi wu about 1 year ago

The problem seems to due to the sc-coul. I have resubmitted the bug in issue 2580 with more clear wording.

Also available in: Atom PDF