gmx hbond for Protein-Ligand H-Bonds in GROMACS 2024+

In GROMACS 2024 and later, count protein-ligand hydrogen bonds with gmx hbond -s md.tpr -f md.xtc -r 'group "Protein"' -t 'resname LIG' -num hbnum.xvg -tu ns. The default geometry is a donor-acceptor distance of 0.35 nm and a hydrogen-donor-acceptor angle of 30 degrees. The old interactive tool still exists as gmx hbond-legacy, and there -r means a distance cutoff, not a selection.

Hydrogen bonds are the third question reviewers ask about a protein-ligand simulation, after ligand RMSD and whether the ligand stayed in contact. The trouble is that most forum answers and older tutorials were written for the pre-2024 gmx hbond, which was rewritten in GROMACS 2024. Copy one of those commands into a current install and it either fails or does something you did not intend. This guide from the StemSkills Lab team (10+ years in structural bioinformatics, drug design and molecular modeling) shows the new syntax, the legacy equivalent, how to get per-pair occupancy, and how to read the result.

What changed in gmx hbond in GROMACS 2024?

GROMACS 2024 replaced gmx hbond with a new implementation built on the selection framework, and kept the old program under the name gmx hbond-legacy. The current gmx hbond manual page says so directly: it is “a new implementation of the hbond utility added in GROMACS 2024.” You can confirm the cut-over yourself: the GROMACS 2023 manual has no hbond-legacy page, and the 2024.0 manual has both.

The practical differences:

• No interactive group prompt. The new tool takes two selections on the command line: -r (reference) and -t (target). The legacy tool still asks you to pick two index groups.
• -r changed meaning. In gmx hbond-legacy, -r 0.35 is the distance cutoff in nm. In the new gmx hbond, -r is a selection. An old command such as gmx hbond -r 0.4 -a 40 no longer works on 2024+.
• Different outputs. The new tool writes -num, -dist, -ang, -dan and an index file -o. The legacy tool also offers -hbn, -hbm (existence matrix), -life and -ac (lifetime and autocorrelation analysis).
• The geometry flags differ by version. In GROMACS 2024 the new tool has fixed criteria of 0.35 nm and 30 degrees, plus an -an switch for nitrogen acceptors. From GROMACS 2025 it gained -hbr (distance), -hba (angle), -de and -ae (donor and acceptor elements, both N O by default).

Check your version before copying anything:

gmx --version | head -n 3
gmx help hbond

What does GROMACS count as a hydrogen bond?

Both tools use a purely geometric criterion. A pair counts as hydrogen bonded when the donor-acceptor distance is at most 0.35 nm and the angle between the donor-hydrogen bond and the donor-acceptor vector is at most 30 degrees. These defaults follow the geometric definition used by Luzar and Chandler for water (Nature 379, 55, 1996), the paper the gmx hbond-legacy documentation cites for its kinetics analysis. The 0.35 nm value sits close to the first minimum of the oxygen-oxygen radial distribution function of liquid water.

Donors and acceptors are found from the topology, not from atom names. By default oxygen and nitrogen are both donors and acceptors, and a donor must carry a bonded hydrogen in the topology. Two consequences matter for ligands:

• The tool needs a run input file with bonds (-s md.tpr), not a bare .gro.
• If your ligand topology lacks polar hydrogens, its OH and NH groups cannot act as donors and the count will be too low.

How do I prepare the trajectory before counting H-bonds?

The new tool makes molecules whole and applies periodic boundary conditions by default (-rmpbc and -pbc are both on), so a raw trajectory gives correct distances. Use the same centred, PBC-corrected trajectory you used for RMSD anyway, so the frames you analyse match the frames you look at. If you have not made one yet, follow our gmx trjconv tutorial.

Next, confirm the ligand residue name. Many tutorials use LIG, but CHARMM-GUI, ACPYPE and other tools often write something else:

grep -m 5 LIG md.gro

If nothing prints, open md.gro and read the residue name next to the last protein residue. A wrong residue name is the most common cause of a flat zero line.

How do I count protein-ligand hydrogen bonds with the new gmx hbond?

Run this on GROMACS 2024 or later:

gmx hbond -s md.tpr -f md_center.xtc \
  -r 'group "Protein"' -t 'resname LIG' \
  -num hbnum.xvg -dist hbdist.xvg -ang hbang.xvg \
  -o hbond.ndx -tu ns

What each part does:

• -r 'group "Protein"' picks the default Protein group. You can use any selection, such as a binding-pocket group from an index file passed with -n.
• -t 'resname LIG' picks the ligand. Reference and target must be either identical or completely non-overlapping.
• -num writes the number of hydrogen bonds per frame. This is the plot you will report.
• -dist and -ang write histograms of the distances and angles of the bonds found, a quick check that your criteria are sensible.
• -o writes an index file with the donors, hydrogens, acceptors and the bonded pairs. It is always written, so naming it keeps it out of the way.
• -tu ns puts the time axis in nanoseconds.

The names “reference” and “target” do not limit direction. Both groups are searched for donors and acceptors, so bonds where the ligand donates and bonds where the protein donates are both counted.

How do I change the distance or angle cutoff?

On GROMACS 2025 or later, set -hbr and -hba. The neighbour-search cutoff must be at least as large as the distance cutoff, so raise both together:

gmx hbond -s md.tpr -f md_center.xtc \
  -r 'group "Protein"' -t 'resname LIG' \
  -hbr 0.40 -cutoff 0.40 -hba 30 -num hbnum_040.xvg -tu ns

On GROMACS 2024, the new tool does not expose these criteria. Use gmx hbond-legacy with -r and -a if you need a different definition. Whatever you choose, state it in your methods section; the defaults are the most widely recognised.

How do I run the same analysis with gmx hbond-legacy?

The legacy tool works with index groups. Create a ligand group first if you do not have one (our gmx make_ndx guide covers this), then run:

gmx hbond-legacy -s md.tpr -f md_center.xtc -n index.ndx \
  -num hbnum.xvg -hbn hbond.ndx -hbm hbmap.xpm -tu ns

GROMACS prints the list of groups and asks for two of them. Choose Protein, then your ligand group. To script it, pipe the group numbers from that menu, for example printf "1\n13\n" | gmx hbond-legacy ..., after checking which numbers your index file uses.

Legacy-only options that are still useful:

• -r 0.35 and -a 30: distance (nm) and angle (degrees) cutoffs.
• -noda: measure the hydrogen-acceptor distance instead of donor-acceptor. Leave the default -da on to match the new tool.
• -hbm hbmap.xpm: an existence matrix with one row per hydrogen bond and one column per frame, in the same order as the bonds in the -hbn index file.
• -life and -ac: lifetime distribution and autocorrelation. These need memory proportional to donors times acceptors, so keep the groups small.

On GROMACS 2023 and older, the same command works under the name gmx hbond, because that was the legacy tool.

How do I get per-pair hydrogen bond occupancy?

The count plot tells you how many bonds exist; occupancy tells you which ones matter. Occupancy is the percentage of frames in which a specific donor-acceptor pair is bonded. A 100 ns trajectory saved every 10 ps has 10,001 frames, so a pair present in 6,500 of them has an occupancy of 65%.

With the new tool, add -pf so each frame gets its own section in the index file:

gmx hbond -s md.tpr -f md_center.xtc \
  -r 'group "Protein"' -t 'resname LIG' \
  -num hbnum.xvg -o hbond_pf.ndx -pf -tu ns

Each frame section is named [ hbonds_..._frame_N ], and each line under it holds three atom numbers: donor, hydrogen, acceptor (two numbers if you add -m to merge bonds that differ only in the hydrogen). This short script counts each donor-acceptor pair once per frame and divides by the number of frames in hbnum.xvg:

import numpy as np
from collections import Counter

nframes = len(np.loadtxt("hbnum.xvg", comments=("#", "@")))
counts, seen, in_frame = Counter(), set(), False
for line in open("hbond_pf.ndx"):
    line = line.strip()
    if line.startswith("["):
        counts.update(seen)
        seen = set()
        in_frame = line.startswith("[ hbonds_")
    elif line and in_frame:
        f = line.split()
        seen.add((int(f[0]), int(f[-1])))
counts.update(seen)

for (d, a), n in counts.most_common(10):
    print(f"donor {d:6d}  acceptor {a:6d}  {100 * n / nframes:5.1f}%")

Atom numbers are 1-based and match a PDB written from the same run with gmx trjconv -s md.tpr -f md_center.xtc -dump 0 -o frame0.pdb, so you can look up residue and atom names there. If you prefer residue-level labels straight away, ProLIF reports hydrogen bond occupancy per residue in Python.

With the legacy tool, the -hbm matrix carries the same information: each row is one bond from hbond.ndx, and the share of frames where it is present is its occupancy. Render it with gmx xpm2ps -f hbmap.xpm -o hbmap.eps for a visual check.

How do I read the hydrogen bond count plot?

Plot hbnum.xvg with xmgrace or Python (see our guide to plotting GROMACS .xvg files). The second column is the H-bond count per frame. Raw counts jump between integers every frame, so add a running average over about 1 ns to see the trend. Four patterns cover most cases:

• Stable, non-zero band (for example 2 to 4 throughout): the key interactions persist. Check occupancy to see whether the same pairs are responsible.
• Fluctuating between 0 and 1 or 2: transient bonds. Common for ligands that bind mainly through hydrophobic contact, and not by itself a sign of instability.
• Step down to zero partway through: a key interaction broke. Compare the time with your ligand RMSD to see whether the ligand shifted pose or left the pocket.
• Zero from the first frame: usually a selection or topology problem rather than chemistry. See the troubleshooting section below.

Report the mean and standard deviation of the count over the equilibrated part of the run, together with the occupancy of the top pairs. A count alone, without knowing which residues are involved, says little about binding.

gmx hbond vs hbond-legacy vs VMD vs ProLIF: which should I use?

The defaults are not interchangeable: VMD’s are stricter than GROMACS’s, so counts from the two will differ on the same trajectory. Pick one tool for a project and quote its criteria. VMD defaults are from the VMD HBonds plugin documentation; ProLIF defaults are from its HBAcceptor and HBDonor interaction classes.

What are the common gmx hbond errors and how do I fix them?

• An old command fails on -r 0.35 or an unknown option -a. You are running a legacy command on GROMACS 2024+. Replace gmx hbond with gmx hbond-legacy, or rewrite it with -r and -t selections.
• “Partial overlap between groups”. Your reference and target share some atoms but are not identical, for example -r 'group "System"' with -t 'resname LIG'. Use Protein and the ligand, or two identical selections for an intra-group count.
• “has no donors AND has no acceptors! Nothing to be done.” The selection contains no N or O atoms the tool can use. Check the residue name, and check that the run input file matches this system.
• Flat zero line, no error. The ligand name is wrong, or the ligand topology lacks polar hydrogens so it has no donors. Add -dan hbdan.xvg to see how many donors and acceptors were found per frame.
• Counts that look too high. Water or ions slipped into a selection such as not group "Protein". Select the ligand by residue name only.
• Invalid cutoff error. On 2025+, an -hbr larger than -cutoff is rejected; raise -cutoff too. On 2024, a -cutoff below 0.35 is rejected.

For problems earlier in the pipeline, see common GROMACS errors and how to fix them.

Where does H-bond analysis fit in a full MD analysis?

Hydrogen bonds are one of three checks for whether a ligand stayed bound: ligand RMSD for pose stability, distance and contacts for whether it is still in the pocket, and H-bonds for which specific interactions hold. Our Rg, SASA and hydrogen bond guide covers intra-protein H-bonds for the protein side of the story. The full sequence is laid out in our learn molecular dynamics with GROMACS pillar, and the wider skill path in the computational biology skills roadmap.

Frequently asked questions

Why does my old gmx hbond command not work anymore?

GROMACS 2024 replaced gmx hbond with a selection-based tool. The old program is now gmx hbond-legacy. Rename the command, or rewrite it with -r and -t selections.

What is the default hydrogen bond cutoff in GROMACS?

A donor-acceptor distance of 0.35 nm and a hydrogen-donor-acceptor angle of 30 degrees, in both gmx hbond and gmx hbond-legacy.

Does gmx hbond count bonds in both directions?

Yes. Both selections are searched for donors and acceptors, so bonds where the ligand donates and bonds where the protein donates are both counted.

How do I calculate hydrogen bond occupancy in GROMACS?

With the new tool, write a per-frame index file with -o and -pf, then count how often each donor-acceptor pair appears and divide by the number of frames. With hbond-legacy, read the -hbm existence matrix row by row.

How many hydrogen bonds should a stable ligand have?

There is no universal number. It depends on the ligand’s chemistry and the pocket. What matters is whether the bonds seen in the docked or crystal pose persist with high occupancy across replicate runs.