[gmx-users] Hydrogen bonding energy calculation

leoneverlasting leoneverlasting at gmail.com
Tue Oct 6 14:19:46 CEST 2015

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!


More information about the gromacs.org_gmx-users mailing list