GROMACS version: 2018
GROMACS modification: Yes/No
Here post your question
Hello everyone,

We’re currently working with GROMACS using thermodynamic integration as a way to obtain Gibbs free energy of the system.

We used the g_analyze command in GROMACS to integrate the dh/dlambda curve for Coulombic interaction of metal ion, the system that we looked at is Mg2+ ion, solvated in TIP4PEW water model. The results after MD with GROMACS using integration are not in agreement with the BAR method in mdtutorial thought we have used the same output xvg to analyze. (By comparing with literature value, Coulombic contribution calculated by BAR is -1738.89, which is extremely close to literature value whereas integrating gave value of -127.96 kJ/mol. )

My question is, do you guys know what causes the difference between two methods?

If you use the slow-growth method, there are some systematic issues and you would expect this to converge badly. The key issue here is that you would expect the system to be in equilibrium at every lambda point that you integrate over, at the same time you keep moving the lambda, destroying the equilibirum and heavily biasing to the state of the system where you just started from at the new lambda position - introducing a large hystersis. To circumvent this you would have to run this infinitely slow, but even approaching very long time-scales, this method is quite inefficient.

If you are using Crooks theorem with “fast-growth” TI, then I’d be curious to hear about the paramters you used - in my experience this one works very well.