[gmx-developers] XVG Energy Discrepancies Between GROMACS Versions

[Miroslav Suruzhon]

Hello,

I was running some free energy perturbation simulations in GROMACS, and 
I had problems reproducing data from other people's simulations in 
GROMACS. After several days of trying to find out the issue, I realised 
that different versions of GROMACS give very different results in the 
.xvg file for the same input .tpr file and .trr trajectory (using gmx 
mdrun -rerun). For example, showing only the first few timesteps for one 
of the perturbations for all versions of GROMACS since 2018.1 (except 
2018.2):

XVG:

@ s0 legend "dH/d\xl\f{} fep-lambda = 0.0000"
@ s1 legend "dH/d\xl\f{} mass-lambda = 0.0000"
@ s2 legend "dH/d\xl\f{} vdw-lambda = 1.0000"
@ s3 legend "\xD\f{}H \xl\f{} to (0.0000, 0.0000, 1.0000)"
@ s4 legend "\xD\f{}H \xl\f{} to (0.2000, 0.0000, 1.0000)"
@ s5 legend "pV (kJ/mol)"


2018.1:

0.0000 -174.60474 20.225147 277.66620 1.9073486e-05 -3.2027580 3.9291420
0.5000 -82.147018 8.3129530 236.27664 1.5258789e-05 -1.2502190 3.9439571
1.0000 -156.04514 11.252308 204.58133 7.6293945e-06 -2.5048264 3.9501503


2018.3:

0.0000 -174.61128 20.225662 277.66617 0.0000000 -3.2027481 3.9291420
0.5000 -82.146118 8.3130722 236.27664 1.5258789e-05 -1.2502289 3.9439571
1.0000 -156.04543 11.252201 204.58139 1.5258789e-05 -2.5047165 3.9501503


2018.4:

0.0000 -174.61128 20.225662 277.66617 2.6702881e-05 -34.366931 3.9291420
0.5000 -82.146118 8.3130722 236.27664 -2.2888184e-05 -16.014955 3.9439571
1.0000 -156.04543 11.252201 204.58139 -1.5258789e-05 -30.665316 3.9501503


2018.5:

0.0000 -174.61128 20.225662 277.66617 0.0000000 -34.366950 3.9291420
0.5000 -82.146126 8.3130722 236.27664 0.0000000 -16.014955 3.9439571
1.0000 -156.04543 11.252201 204.58139 0.0000000 -30.665323 3.9501503


2019:

0.0000 -18.344168 0.0000000 277.66617 1.1444092e-05 -3.1135101
0.5000 -7.9936657 0.0000000 236.27664 0.0000000 -1.1844406
1.0000 -15.156097 0.0000000 204.58139 3.0517578e-05 -2.4874344


Ignoring the floating point errors, we can see that 2018.1 and 2018.3 
give more or less the same results, which are very different from those 
produced by 2018.4 and 2018.5, resulting in extremely large 
discrepancies in free energies - more than 20 kJ/mol (not shown here, 
but the other .xvg files exhibit the same behaviour).  However, these 
discrepancies only seem to hold for the delta H values and not dH/dl. 
Finally, free energies calculated from 2019 using gmx bar seem to 
reproduce 2018.1 and 2018.3 but we can see that there is a missing 
column of data compared to other versions and most of the dH/dl 
derivatives do not agree with other versions of GROMACS. Note that this 
issue only seems to persist when perturbing the Coulomb interactions. 
Van der Waals perturbation energies seem consistent between versions 
(not shown here).


Since the only variable for all of these simulations is the GROMACS 
version (which was compiled on the same computer with the same 
libraries), I highly suspect that there is (at least one) bug with the 
free energy implementation of the code (or the output thereof). However, 
I am not sure which version gives the "correct" results. Any feedback 
would be highly appreciated. I have not uploaded any files but if any 
more information is required, I am happy to provide it.


Thank you very much.


Kind regards,

Miroslav