[gmx-users] Hydrogen bonding energy calculation

[leoneverlasting]

Dear gromacs users,

I have a question regarding the hydrogen bonding energy calculation. I am trying to reproduce hydrogen bonding energy in David’s paper from J. Phys. Chem. B 2006, 110, 4393-4398.

I set up my systems using OPLS force field, with TIP3P, SPC and SPCE water models. My systems just contain water molecules. The results I got for tip3p are as follows:

Doing autocorrelation according to the theory of Luzar and Chandler.
Normalization for c(t) = 0.000585267 for gh(t) = 2.92622e-08
Hydrogen bond thermodynamics at T = 298.15 K
Fitting parameters chi^2 =  0.0139988
Q =          0
--------------------------------------------------
Type      Rate (1/ps) Time (ps)  DG (kJ/mol)  Chi^2
Forward         0.534      1.873       6.083   0.0139988
Backward        0.632      1.583       5.667
One-way         0.605      1.652       5.773
Integral        0.326      3.070       7.309
Relaxation      0.534      1.871       6.081

where I think 6.083 kJ/mol is the activation energy for hydrogen bonding. However this value is much larger than 3.95 kJ/mol as in the 2006 JPCB paper. 

A few questions that I have are:
1. Does the version of gromacs matters? In the paper, GROMACS 3.3.1 was used for simulation and GRIMACES 3.3 was used for analysis. I am using GROMACS 5.0.1.
2. All the other systems parameters are set up as described in Erik’s 2003 paper in J. Chem. Phys. 119, No 14, 8, except that I use SHAKE algorithm to constrain the water molecules, they use SETTLE. Could anyone elaborate on the difference that might cause?

Thanks for all your help in advance!

Ben