Salt Bridge Analysis in GROMACS: gmx saltbr vs pairdist

To analyse a salt bridge in a GROMACS trajectory, track the minimum distance between the acidic side-chain oxygens (Asp OD1/OD2, Glu OE1/OE2) and the basic side-chain nitrogens (Lys NZ, Arg NE/NH1/NH2) with gmx pairdist, then count the frames under 0.4 nm to get occupancy. Treat gmx saltbr as a rough scan on a protein-only system, not as the answer.

After hydrogen bonds, the next question a supervisor asks about a simulation is usually about one specific charged pair: does the Asp-Lys bridge seen in the crystal structure survive, or does the ligand’s carboxylate keep its contact with an arginine in the pocket? GROMACS ships a tool with the obvious name, gmx saltbr, and most students run it first. It rarely gives them what they wanted. This guide from the StemSkills Lab team (10+ years in structural bioinformatics, drug design and molecular modeling) explains what gmx saltbr actually computes, then walks through the distance-based method that gives a clean, defensible occupancy number.

All flags below were checked against the GROMACS 2026.4 manual. If you are new to trajectory analysis, start with the complete guide to molecular dynamics with GROMACS and come back here.

What counts as a salt bridge in an MD simulation?

A salt bridge is a close-range electrostatic contact between two oppositely charged side chains, usually Asp or Glu on one side and Lys, Arg or protonated His on the other. Most analyses use a geometric criterion: a carboxylate oxygen and a basic nitrogen within about 4 Å (0.4 nm).

That 4 Å figure comes from structural surveys. In Barlow and Thornton’s 1983 study “Ion-pairs in proteins” (J. Mol. Biol. 168:867-885), the authors derived a working definition of 4 Å or less between charged groups from distance distributions in 38 proteins. Under that definition, on average one-third of charged residues took part in ion pairs, and only 17% of those pairs were buried. Most salt bridges sit on or near the surface, where water competes for them, which is exactly why they break and re-form during a simulation.

Kumar and Nussinov’s review, “Close-range electrostatic interactions in proteins” (ChemBioChem, 2002), makes the same point from the dynamics side: “salt bridges and their stabilities fluctuate in proteins.” The same review notes that a salt bridge can be stabilizing or destabilizing, so a bridge that breaks in your trajectory is a result to report, not automatically a failed simulation.

Different tools use different cutoffs, so always state yours:

What does gmx saltbr actually compute?

gmx saltbr plots the distance between charged groups as a function of time. Its whole option list is short: -f, -s, -b, -e, -dt, -t (default 1000) and -[no]sep. There is no -n index option and no selection option, so you cannot point it at one residue pair or at a ligand.

The critical detail is how it defines a “charged group”. In the current gmx_saltbr.cpp source, every atom whose partial charge in the .tpr is non-zero becomes its own group, and the tool stores the distance between every pair of such atoms for every frame. In an all-atom force field almost every atom carries a partial charge. The consequences:

• It is not a salt-bridge detector. plus-min.xvg holds every oppositely charged atom pair that came closer than -t, including backbone carbonyl oxygens and amide nitrogens. Those are not salt bridges.
• Water and ions are included if you pass the full solvated .tpr.
• Memory grows with the square of the atom count. Rough arithmetic: 3,000 charged atoms give about 4.5 million pairs, and at 4 bytes per single-precision distance that is about 18 MB per frame, or about 18 GB for 1,000 frames.
• -t defaults to 1000, and GROMACS distances are in nm, so the default filters nothing.

Outputs go to fixed names: plus-plus.xvg, min-min.xvg and plus-min.xvg. With -sep you get one file per pair named from residue name, residue number and atom number (the manual pattern is sb-(Resname)(Resnr)-(Atomnr)), and the manual itself warns there “may be MANY”.

Can I still use gmx saltbr for a first scan?

Yes, if you shrink the problem first. Make a protein-only run input file and trajectory, thin the frames, and set a real -t:

# 1. protein-only .tpr (pick "Protein" at the prompt)
gmx convert-tpr -s md.tpr -o protein.tpr

# 2. matching protein-only trajectory, molecules made whole
gmx trjconv -s md.tpr -f md.xtc -o protein.xtc -pbc mol -dt 100
#    select "Protein" for output

# 3. scan: only pairs that ever come within 0.4 nm are written
gmx saltbr -s protein.tpr -f protein.xtc -t 0.4

The manual flags a subset .tpr from gmx convert-tpr as “not fully functional”, which is fine for analysis but not for running mdrun. Open plus-min.xvg, read the legend lines (each one names the two atoms as residue name, number and atom number), and note the side-chain pairs you care about. Ignore backbone pairs. You now have a candidate list, and the real measurement starts in the next step. If you are unsure how -pbc mol and -dt behave, see the gmx trjconv tutorial.

How do I track one salt bridge with gmx pairdist?

gmx pairdist computes the minimum distance (the default -type min) between a reference selection (-ref) and one or more other selections (-sel), with periodic boundaries on by default. That is exactly the salt-bridge question: how close is the nearest acidic oxygen to the nearest basic nitrogen in this frame?

Suppose your scan pointed at Asp 45 and Lys 120 (replace these with your own residue numbers):

gmx pairdist -s md.tpr -f md.xtc \
  -ref 'resid 45 and name OD1 OD2' \
  -sel 'resid 120 and name NZ' \
  -o sb_asp45_lys120.xvg -tu ns

Include both carboxylate oxygens. OD1 and OD2 are chemically equivalent and swap roles constantly, so tracking only one of them makes a stable bridge look like it breaks every few picoseconds. The atom names for the common residues in AMBER, CHARMM and OPLS-AA topologies:

• Asp: OD1 OD2; Glu: OE1 OE2
• Lys: NZ; Arg: NE NH1 NH2
• His: ND1 NE2, but only when the residue is protonated on both (see troubleshooting)

resid is a synonym for resnr in the GROMACS selection syntax and uses the numbering in your input file. In a multi-chain system where numbers repeat, add a chain term or use resindex. If you prefer named groups, build them with gmx make_ndx and pass -n index.ndx, as shown in our make_ndx index groups tutorial.

Why not gmx distance? It measures fixed atom pairs (1-2, 3-4 and so on) within one selection, so it cannot take the minimum over two oxygens and three arginine nitrogens. The manual itself points you to gmx pairdist for minimum distances between two selections.

How do I check a salt bridge between the protein and a charged ligand?

Same command, with the ligand as one side. Take the charged atom names from your ligand topology (.itp or .mol2), because they are whatever your parameterization tool assigned:

gmx pairdist -s md.tpr -f md.xtc \
  -ref 'resname LIG and name O1 O2' \
  -sel 'resid 210 and name NE NH1 NH2' \
  -o sb_lig_arg210.xvg -tu ns

Here O1 O2 stand in for a ligand carboxylate and Arg 210 for a pocket arginine. A ligand only forms a salt bridge if its topology really carries a formal charge, so check the net charge in the .itp before you look for one. For a fingerprint of every interaction type at once, including ionic ones, ProLIF interaction analysis is the better tool, and hydrogen bonds to the same ligand are covered in our gmx hbond protein-ligand guide.

How do I calculate salt bridge occupancy in Python?

Occupancy is the percentage of frames in which the minimum O-N distance is at or below your cutoff. This script reads the .xvg, skips the # and @ header lines, and reports occupancy, the mean distance and the longest unbroken stretch:

import numpy as np

def load_xvg(path):
    rows = [line.split() for line in open(path)
            if line.strip() and not line.startswith(("#", "@"))]
    return np.array(rows, dtype=float)

data = load_xvg("sb_asp45_lys120.xvg")
time, dist = data[:, 0], data[:, 1]   # gmx pairdist writes nm
cutoff = 0.40                          # 4.0 Angstrom O-N criterion

formed = dist <= cutoff
print(f"Frames analysed: {len(dist)}")
print(f"Occupancy: {100 * formed.mean():.1f} %")
print(f"Mean O-N distance: {dist.mean():.3f} nm")

# longest unbroken stretch with the bridge formed
runs, current = [], 0
for f in formed:
    current = current + 1 if f else 0
    runs.append(current)
dt = time[1] - time[0]
print(f"Longest continuous contact: {max(runs) * dt:.1f} (time units of the file)")

Report occupancy together with the cutoff and the time window, for example “Asp45-Lys120 salt bridge, 0.40 nm O-N cutoff, 78% occupancy over the last 80 ns”. Drop the equilibration period first with -b on the gmx pairdist command. To plot the distance trace, follow plotting GROMACS xvg files with Python and draw a horizontal line at 0.4 nm.

How do I read the salt bridge distance plot?

• Flat and well below 0.4 nm: a tight, persistent bridge. This is the classic direct contact.
• Switching between a low plateau and a second plateau just above the cutoff: the bridge breaks and re-forms. Often a water molecule slips between the two groups, giving a water-mediated contact. Report occupancy rather than calling it either stable or lost.
• Starts low, rises and stays well above the cutoff: the bridge broke for good. Check whether the side chains rotated away or whether a larger conformational change is happening, using RMSD and the radius of gyration and SASA analysis.
• Isolated spikes to several nm for one frame: almost always a periodic boundary artifact, not chemistry (see troubleshooting).

gmx saltbr vs gmx pairdist vs VMD vs MDAnalysis: which should I use?

A practical split: use the VMD plugin or a thinned gmx saltbr run to find candidates, then measure each candidate with gmx pairdist and the occupancy script. Note that the VMD plugin’s written distance (centre of mass of the side-chain oxygens to centre of mass of the nitrogens) differs from its detection criterion, as its documentation says, so do not mix its numbers with gmx pairdist minimum distances in one table.

What are the common salt bridge analysis problems and how do I fix them?

• gmx saltbr runs for hours or runs out of memory. It is storing every pair of partially charged atoms, solvent included. Use a protein-only .tpr and trajectory from gmx convert-tpr and gmx trjconv, and thin frames with -dt.
• Thousands of sb- files appear. You used -sep with the default -t 1000, so every pair was written. Set -t 0.4 or drop -sep.
• plus-min.xvg is full of backbone atoms. Expected: the tool works on atoms with partial charges, not on charged residues. Filter the legend for side-chain atom names or switch to gmx pairdist.
• A bridge from the crystal structure never forms. Check protonation. A histidine built as neutral (HID or HIE in AMBER naming, HSD or HSE in CHARMM) cannot form a salt bridge, and an Asp or Glu you protonated is no longer charged. gmx pdb2gmx has -his, -asp, -glu, -lys and -arg flags for interactive choices; see protein protonation states for docking and MD.
• The distance jumps to several nm for single frames. A periodic boundary artifact. Keep the default -pbc on in gmx pairdist; if you also ran -nopbc or analysed a trajectory that was centred on something else, re-process it with gmx trjconv -pbc mol.
• The bridge looks broken half the time but the structure has not changed. You probably tracked one carboxylate oxygen. Select both OD1 OD2 (or OE1 OE2) and let the minimum handle the swap.
• Occupancy is 100% or 0% for everything. Check units: GROMACS writes nm, so the cutoff is 0.4, not 4.

Where does salt bridge analysis fit in a full MD analysis?

It is a second-order check. Run it after you have confirmed the system is equilibrated and the protein is stable (RMSD, RMSF, radius of gyration), and alongside hydrogen bond analysis, since both answer “which specific contacts hold this structure or this complex together”. The order of skills, from setup to analysis, is laid out in the computational biology skills roadmap.

Frequently asked questions

What cutoff should I use for a salt bridge in GROMACS?

0.4 nm (4 Å) between a side-chain carboxylate oxygen and a basic nitrogen is the most widely used choice and follows Barlow and Thornton’s 4 Å definition. VMD’s plugin defaults to a stricter 3.2 Å. Whichever you pick, state it in your methods.

Does gmx saltbr accept an index file?

No. Its only options are -f, -s, -b, -e, -dt, -t and -sep. To restrict it, make a reduced .tpr with gmx convert-tpr and a matching trajectory with gmx trjconv.

What is the difference between a salt bridge and a hydrogen bond in MD?

A salt bridge needs two formally charged groups; a hydrogen bond needs a donor, a hydrogen and an acceptor within distance and angle limits. Many salt bridges also satisfy hydrogen bond geometry, so the same contact can appear in both analyses.

Can a salt bridge be water-mediated?

Yes. If a water molecule bridges the two charged groups, the direct O-N distance sits above 0.4 nm and a 0.4 nm cutoff counts it as broken. Report direct and water-mediated contacts separately if the distinction matters to your conclusion.

How long should the simulation be to judge salt bridge stability?

Long enough to see the bridge break and re-form several times, or to show it does not. A single short run gives one sample; replicate simulations with different starting velocities give a much more convincing occupancy.