Incorrect relative solvation free energy values using single topology alchemical transformation

User discussions
, ,
1

GROMACS version: 2022.5

GROMACS modification: No

Dear gromacs users,

I am currently trying to perform FEP calculations. More precisely, I want to calculate the relative solvation free energy between ethane and ethanol using a single topology approach. I generated an hybrid topology for the ethane → ethanol alchemical transformation using the alchemical-setup tool developped by D.Mobley and P.Klimovich:
Klimovich, P.V., Mobley, D.L. A Python tool to set up relative free energy calculations in GROMACS | Journal of Computer-Aided Molecular Design J Comput Aided Mol Des 29, 1007–1014 (2015).

However, the results I am getting (using the MBAR method and the alchelyb library for python) are far from the expected values: (-29.37 +/- 0.02) kcal/mol when I am supposed to obtain a value around -6 kcal/mol according to Mobley’s Hydration free energy datbase (you can find it here : FreeSolv/database.txt at master · MobleyLab/FreeSolv · GitHub )

When doing absolute solvation free energy calculations, I do get the correct values : (2.18 +/- 0.02 kcal/mol) kcal/mol for ethane and (-4.36 +/- 0.02 kcal/mol) for ethanol, which makes me think that I am not giving gromacs the correct instructions for the alchemical transformation. You can find here the hybrid topology .itp and the .mdp file I am using to perform my relative solvation free energy calculations.

ethane_ethanol.itp (6.2 KB)

06-NPT-SFE_ethane_ethanol.mdp (4.7 KB)

I’ve tried to do the same relative calculations by generating the hybrid topology using another tool (pmx) and also doing a simpler alchemical transformation : ethanol → ethanethiol which does not involve any dummy atom, and I still get wrong values (again around 30 kcal/mol when i’m supposed to get around 3 kcal/mol )

Do you have any insight on what I am doing wrong ?

Just in case, I am also proving the .mdp file as well as the .itp files I used for the absolute calculations (here for ethanol but the ethane ones are the same). Do not hesitate to ask for other information should you need them.

06-NPT_ethanol_FEP-Coul-VdW.mdp (5.8 KB)
ethanol.top (5.0 KB)

Thank you in advance.
Best regards,
NicolasL

2

Hello everyone,
Since my post was automatically hidden when I first published it (and I’m sorry if I did not comply with the community guidelines), I would like to reiterate my question. Does anyone have any insight on why I do not obtain the correct values for my FEP calculations?

Thank you in advance !
NicolasL

3

One thing I’ve noticed is that in the .mdp settings, you’re transforming the Coulombic and VdW components of the Hamiltonian, but your bonded and mass terms remain in state A, or lambda=0. That means you’re not changing, for example, the C-H/C-O bond length, nor adjusting the dihedrals/angles that are supposed to change between ethane and ethanol. Not sure if this can produce such a large discrepancy in the results, but I’d look into this.