[gmx-users] What does the -don switch calculate in g_hbond?

[Andrew DeYoung]

Hi,

In g_hbond, there is an output switch -don, which by default gives an output
file named donor.xvg.  What does this -don switch calculate?

For example, my system contains 22 donors (i.e., oxygens) and 22 hydrogens
-- in the form of 22 hydroxyl (OH) groups.  My system also contains 1024
acceptors -- in the form of 256 PO4^3- ions (each of which possesses 4
acceptors, i.e., oxygens).  So the two groups I select via an index file and
pass to g_hbond are "GroupOH" and "GroupPO4".  Thus, oxygen acts as a donor
in OH groups and as an acceptor in PO4 groups. 

The donor.xvg output from the -don switch gives this sort of output:

@    title "Donor properties"
@    xaxis  label "Time (ps)"
@    yaxis  label "Number"
@ s0 legend "Nbound"
@ s1 legend "Nfree"
 0.000e+00      20      2
 5.000e-02      20      2
 1.000e-01      20      2
 1.500e-01      20      2
 2.000e-01      20      2
 2.500e-01      21      1
 3.000e-01      21      1
 3.500e-01      21      1
 4.000e-01      21      1
 4.500e-01      21      1
 5.000e-01      21      1
 5.500e-01      21      1
 6.000e-01      21      1
...

What are "Nbound" and "Nfree"?  I guess that they are simply the number of
donors involved in hydrogen bonds and the number of donors NOT involved in
hydrogen bonds, respectively.  I have 22 donors (i.e., 22 oxygens on 22 OH
groups), so it makes sense that the columns add up to 22 at each time frame
in my trajectory (I saved configurations every 0.050 ps).  Evidently, the
vast majority of OH groups are involved in hydrogen bonds, since Nbound >>
Nfree at every frame.

However, if OH groups are the only source of donors (O) and PO4 groups are
the only source of acceptors (O), then shouldn't Nbound be equal to the
number of hydrogen bonds at every frame?  This does not seem to be the case.
If I look at the number of hydrogen bonds over time (in hbnum.xvg, from the
-num switch in g_hbond), I get this output: 

@    title "Hydrogen Bonds"
@    xaxis  label "Time (ps)"
@    yaxis  label "Number"
@ s0 legend "Hydrogen bonds"
@ s1 legend "Pairs within 0.227905 nm"
         0          20           0
      0.05          20           1
       0.1          20           0
      0.15          19           0
       0.2          20           0
      0.25          20           0
       0.3          20           0
      0.35          20           0
       0.4          19           1
      0.45          18           0
       0.5          19           2
      0.55          18           1
       0.6          17           1
...

So, the "Hydrogen Bonds" column in hbnum.xvg does not match the "Nbound" in
donor.xvg.  Is this possible, given the "GroupOH" and "GroupPO4" groups
(which are non-identical and non-overlapping) I have passed to g_hbond?

Thank you very much for your time!

Andrew DeYoung
Carnegie Mellon University