Number of HBonds in TIP3P water

[laraujo-1]

Hi all,

I am looking for some user guidance on how to obtain reference hydrogen bond (H-bond) values in water using MDAnalysis.

The number of H-bonds per molecule in TIP3P water is well-known to be about 3.3, as mentioned in the MDAnalysis user guide. I performed a 10 ns simulation of TIP3P water with GROMACS, and using gmx hbond with a donor–acceptor cutoff (r_DA = 3.5 Å) and angle cutoff (a_HAD = 30°), I obtained:

HBonds/water = 3.32 ± 0.02

However, I am having difficulty reproducing this value using MDAnalysis.

I understand that MDAnalysis uses the D–H–A angle (rather than H–D–A). I have seen some people using a 120° cutoff for this angle, but even then, the results are not consistent with the reference value. For the same donor–acceptor cutoff of 3.5 Å:

• Using an angle cutoff of 120° → value too high (~4.10)
• Using an angle cutoff of 150° → value too low (~2.56)

Here is the code I am using:

import MDAnalysis as mda
from MDAnalysis.analysis.hydrogenbonds import HydrogenBondAnalysis as HBA

# --- Load trajectory ---
topology = "./1-water_TIP3p/md.tpr"
trajectory = "./1-water_TIP3p/md.xtc"
u = mda.Universe(topology, trajectory)

# --- H-bond parameters ---
r_cut = 3.5       # Å
angle_cut = 150.0 # degrees
step_size = 100

# --- H-bond analysis ---
h = HBA(
    universe=u,
    donors_sel="name OW",             
    hydrogens_sel="name HW1 HW2",
    acceptors_sel="name OW",
    d_a_cutoff=r_cut,
    d_h_a_angle_cutoff=angle_cut,
    update_selections=True)

n_waters = len(u.select_atoms("resname SOL and name OW"))
h.run(step=step_size)

# --- Compute number of H-bonds per water ---
hb_per_mol = 2 * h.count_by_time() / n_waters
hb_average = np.mean(hb_per_mol)
print("Number of hydrogen bonds per water molecule over time", hb_average)

Could someone provide guidance on how to get the appropriate parameters to reproduce literature values for TIP3P water.
Also any additional considerations for MDAnalysis H-bond detection compared to GROMACS gmx hbond I think would be valuable for the community.

Thanks in advance and best regards,

Lucas

[p-j-smith]

Hi @laraujo-1, have you checked the indices of hydrogen bonds found by gmx hbond and MDAnalysis? It would be good to know which hydrogen bonds were not found by MDAnalysis when setting d_h_a_angle_cutoff=150, and what the distributions of D-H-A and H-D-A angles look like for these missing hydrogen bonds.

Also, as you're selecting all water molecules at all frames, you can set update_selections=False, and as your topology contains bonds you can omit the donor_sel as donor atoms directly from hydrogen_sel (any atom bonded to a hydrogen in hydrogen_sel will be a donor). So you can change your code to:

h = HBA(
    universe=u,
    hydrogens_sel="name HW1 HW2",
    acceptors_sel="name OW",
    d_a_cutoff=r_cut,
    d_h_a_angle_cutoff=angle_cut,
    update_selections=False,
)

[laraujo-1]

Hi @p-j-smith,

I followed your suggestion. Basically, GROMACS identifies the same hydrogen bonds as MDAnalysis plus some additional ones when using d_h_a_angle_cutoff=150. When using d_h_a_angle_cutoff=120, the opposite happens. What I did was collect the indexes into files and compare them. Below I am using d_h_a_angle_cutoff=150 and d_a_cutoff=3.5:

I also analyzed the distribution of angles and distances for the hydrogen bonds that MDAnalysis misses. As expected, all distance criteria were satisfied, the donor-acceptor distances for the missing hydrogen bonds are always below the cutoff. The same is also true of the H-D-A angle cutoff of gromacs.

Regarding the angle distributions, it’s a bit trickier. The missing hydrogen bonds have DHA angles ranging from 135° to just below 150°:

From a quick check, using a DHA angle cutoff of ~135° one misses almost nothing. However, I would be cautious because this analysis was done only for the last frame of the trajectory of a box of water. Nevertheless, this suggests that the default value of 150° for the DHA angle might deserve some review.

[p-j-smith]

Nice analysis @laraujo-1!

From a quick check, using a DHA angle cutoff of ~135° one misses almost nothing. However, I would be cautious because this analysis was done only for the last frame of the trajectory of a box of water. Nevertheless, this suggests that the default value of 150° for the DHA angle might deserve some review.

Changing the default value would be a breaking change, which we could do for the 3.0 release. We should at least make a note of the of this in the docs or user guide, as the difference between results from gmx hbond and MDAnalysis is a fairly common question. Or we could change HydrogenBondAnalysis to use the H-D-A angle.

[laraujo-1]

Thanks @p-j-smith.

Changing the default value would be a breaking change, which we could do for the 3.0 release. We should at least make a note of the of this in the docs or user guide, as the difference between results from gmx hbond and MDAnalysis is a fairly common question. Or we could change HydrogenBondAnalysis to use the H-D-A angle.

I think that's a great idea and If I may suggest something, I think it would be nice to have the H–D–A angle as the default, with a cutoff of 30°, following some reference studies (such as the 3.3 value for TIP3P water). This would increase confidence that the code behaves as expected. But it is also a good opportunity to add extra flexibility in the code by allowing the user to select either the H–D–A or D–H–A angle (or both) and set their own cutoff. Since there are so many MD codes with different implementations of hydrogen bond analysis, it would be great if MDAnalysis could be the place where all of them meet.

Another point is that many people use MDAnalysis to analyze hydrogen bonds from AIMD simulations or trajectories generated with machine learning potentials. Providing this flexibility would make it easier to compare results across different studies and improve reproducibility.

Looking forward to the 3.0 release!