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
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
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.
Thank you for your answer !
I did the same calculations with a corrected lambda vector taking into account the change in the mass and bonded parameters, but I still get wrong values (-29.48 +/- 0.03 kcal/mol). Since it’s really close to the previous result, the cause of the error seems to be the same.
Folllowing your idea, and after reading again the .mdp section of the gromacs user guide, I also tried to separate all the contribution changes in the lambda vector to see if I could get a hint about where the issue is coming from. I first changed the Coulombic, then VdW, then the mass and finally the bonded components (see below the lambda vectors).
I’m not really sure if the results I obtain make sense and are reliable, but by computing the free energy differences for each window with respect to the first one (ethane), I can see that changing the Coulombic components leads to an absolute difference of around 46 kcal/mol, which then goes down to ~32 kcal/mol when changing the VdW interactions. Changing the mass and bonded components then leads to the ~29,5 kcal/mol final value. So my guess is that the problem comes from the Coulombic interactions.
Does anyone have any thoughts on that ?
; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
coul_lambdas = 0.00 0.20 0.40 0.60 0.80 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
vdw_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.20 0.40 0.60 0.80 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.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.20 0.40 0.60 0.80 1.00 1.00 1.00 1.00 1.00 1.00
bonded_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.20 0.40 0.60 0.80 1.00
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 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
temperature_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
Try with couple-intramol = no, otherwise you should do the same perturbations in vacuum as well.
EDIT: Now that I think of it, I must admit I don’t know how couple-intramol works with mutation perturbations. At least it’s worth testing if it makes a difference.
Well that’s a good point, I initially assumed you do the transformation in both solvent and vacuum and then subtract the results to close the thermodynamic cycle, but is this the case? If not, then that would be an obvious source of error.
1 Like
The more I’ve been thinking about this (I’m personally almost never running relative (mutation) free energy calculations) I’m quite certain you should always run a reference state in vacuum. I’m quite sure that’s what milosz mean as well.
1 Like
Thank you very much for your assistance !
By doing the same transformation in vacuum, I did obtain a value of -22.86 +/- 0.03 kcal/mol that, when susbtracted from my original value, gives me the expected result of ~6 kcal/mol. I do admit that, even though it seems logical to me it is required to close the thermodynamic cycle, I thought this vacuum transformation would lead to a value of 0, since there is no solvent to interact with. I did not expect such a value, that seems high for a molecule interacting with itself, but I think I’m mixing some things up…
Just as a side note, I also tried with couple_intramol=no and the values are almost unchanged ( -29.30 +/- 0.03 kcal/mol and -29.34 +/- 0.03 kcal/mol). I looked again on the forum for topics about this parameter and found this clear answer by J.Lemkul : Relationship between couple-moltype and couple-intramol - #5 by jalemkul
I’m quoting him directly (I don’t know if there is a better way to do it):
couple-intramol is integral to the calculation, it specifies whether intramolecular force field terms are coupled as a function of λ. For an absolute free energy calculation, it needs to be set to no but for a transformation like you’re describing, in which the molecular structure changes, it should be set to yes.
Anyway, thanks again !
Yeah, the value obtained for the transformation in vacuum is just the free energy associated with adding/removing classical intramolecular interaction terms - internal electrostatics/vdW as well as the bonded terms energy, it’s almost always going to be non-zero so has to be corrected for. Happy there were no further errors!
1 Like