[gmx-users] validation of g_hbond

[Moore, Jonathan (J)]

There seems to be much confusion in the gmx_users archive about how H bonds are defined GROMACS, as it has changed between 3.1 and 3.2 and isn't always accurately reflected in the documentation (3.2.1 manual still lists the cut-off angle as 60 and g_hbond -h refers to the Donor - Hydrogen - Acceptor angle).

In a message to gmx-users on Dec 5, 2003 (http://www.gromacs.org/pipermail/gmx-users/2003-December/008182.html), David wrote:
"When I quite recently tried to reproduce literature results of water simulation and found some oddities, I looked into all this stuff. The number of hbonds can now be reproduced exactly compared to literature values."

David, could you give us the details of this validation exercise?  For example, to exactly which literature results did you compare?  Also, can you provide details (or references that do) concerning the justification of the acut and rcut values?  I think I saw somewhere the 0.35 nm is the location of the first minimum in the O-O rdf in SPC water.  What about the choice of 30 for the angle?

I'd like to know more details about how g_hbond was validated (especially given the uncertainty in the email archives).

Looking in the source code:

In g_hbond, the existence of an H bond is defined by a distance being less than or equal to a certain value (rcut) and an angle being less than or equal to a certain value (acut).  For H=hydrogen, D=donor, A=acceptor, looking at the source code, it looks to me like the current default definition is:

acut=30, rcut=0.35, a is the angle between the DH vector and the DA vector, r is the distance between D and A

That definition is changed significantly from v 3.1, apparently to bring the results into agreement with literature values.

As a command line option, you can change the value of acut, change the value of rcut, and change the distance to be that between H and A instead of between D and A.  However, the angle that the documentation refers to (e.g., in g_hbond -h from v 3.3rc3), is "Donor - Hydrogen - Acceptor".  I don't think that is the angle used by the code.  Instead, it appears to me to be "Hydrogen - Donor - Acceptor", i.e. the angle between the DH vector and the DA vector.  

This message from October 2005 refers to a serious error in v 3.3
http://www.gromacs.org/pipermail/gmx-users/2005-October/017634.html
Exactly which g_hbond results are affected by this error?

This message from September 2005 reports a discrepancy between GROMACS and VMD regarding the calculation of the number of H bonds:
http://www.gromacs.org/pipermail/gmx-users/2005-September/017293.html
I didn't see a response to the message.  Do GROMACS and VMD use the same H bond definition?

This message from February 2005 reports a discrepancy between GROMACS and a homegrown code regarding the calculation of the number of H bonds:
http://www.gromacs.org/pipermail/gmx-users/2005-February/014060.html

Thanks,
Jonathan

____________________________
Jonathan Moore, Ph.D.
Research and Engineering Sciences - New Products
Core R&D
The Dow Chemical Company
1702 Building, Office 4E
Midland, MI 48674  USA
Phone:  (989) 636-9765 
Fax: (989) 636-4019
E Mail: jmoore2 at dow.com