Hydrogen bond calculation

[horlust]

Hi everyone!

I am computing the hydrogen bonds within a system containing 510 water molecules. I do something like this:

hbonds = HBA(universe=u, d_a_cutoff=3.5, d_h_a_angle_cutoff=30)

hydrogens_sel = hbonds.guess_hydrogens("resname SOL")
acceptors_sel = hbonds.guess_acceptors("resname SOL")
hydrogen_donors_sel = "resname SOL and name OW"

hbonds.hydrogens_sel = f"({hydrogens_sel})"
hbonds.acceptors_sel = f"({acceptors_sel})"
hbonds.donors_sel = f"({hydrogen_donors_sel})"

hbonds.run(verbose=True)

However, the results are:

The amount of hydrogen bonds per molecule of water is of around 9. It is very unlikely, and different from the output I get using GROMACS:

Do you know what could be the problem?

Thanks!

[remaker789]

The definitions of angle are different between MDAnalysis and gmx hbond. MDAnalysis uses the D-H-A angle(see this), but gmx hbond use the H-D-A angle(see this). You can use 'd_h_a_angle_cutoff=150' in MDAnalysis, however, the results of hydrogen bonds from MDAnalysis will still be less than gmx bond due to the definition of angle.