[gmx-users] validation of g_hbond

[David]

On Tue, 2005-10-25 at 11:51 -0500, Moore, Jonathan (J) wrote:
> 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?
Xu et al., J. Phys. Chem. B 2002, 106, 2054-2060
containing TIP4P and SPC/E simulation. It's important to read the small
print, e.g. about whether NVT or NPT was used etc.
> 
> 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 is all correct. The original implementation was based on the
definition used by Kabsch and Sander, Biopolymers 22 (1983) p. 2577,
which used the 60 deg. angle.

> 
> 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?
-ac -map


> 
> 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?
I don't know.

> 
> 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
There are still problems with the -sel option. There is an open bugzilla
at: http://bugzilla.gromacs.org/show_bug.cgi?id=25



-- 
David.
________________________________________________________________________
David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group,
Dept. of Cell and Molecular Biology, Uppsala University.
Husargatan 3, Box 596,          75124 Uppsala, Sweden
phone:  46 18 471 4205          fax: 46 18 511 755
spoel at xray.bmc.uu.se    spoel at gromacs.org   http://xray.bmc.uu.se/~spoel
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++