Ligand RMSD in GROMACS: Fit, Fix PBC, Read the Plot

To calculate ligand RMSD in GROMACS, first re-center the trajectory with gmx trjconv -pbc mol -center so the ligand stops jumping across the box, then run gmx rms choosing Backbone as the fit group and the ligand’s heavy atoms as the RMSD group. The result measures how far the ligand moved relative to the protein pocket.

This is the single most common analysis students get wrong after a protein-ligand run. The command looks trivial, but two choices (what you fit on, and fixing periodic boundaries first) decide whether your plot shows a stable binding pose or a meaningless 4 nm spike. Below is the exact sequence the StemSkills Lab team uses, checked against the official gmx rms and gmx trjconv manual pages (GROMACS 2026.4) and Justin Lemkul’s protein-ligand tutorial. If you are still setting up the simulation itself, start with our protein-ligand MD simulation tutorial in GROMACS and come back once you have a production trajectory.

What does ligand RMSD actually measure?

Ligand RMSD is the root mean square deviation of the ligand’s atom positions in each frame from their positions in a reference structure, after the frames have been superimposed. The question it answers depends entirely on what you superimpose on.

When you fit on the protein backbone and then measure the ligand, the protein’s overall rotation and translation are removed, so the number reports how much the ligand moved inside the pocket. Lemkul’s tutorial puts it plainly: this RMSD is “a good indicator of how well the binding pose was preserved.” That is the quantity reviewers expect in a docking-plus-MD paper.

When you fit on the ligand itself, you remove the ligand’s own translation and rotation, so all you see is internal flexibility (ring puckers, torsion changes). A ligand that drifts completely out of the pocket can still give a flat, low RMSD this way. This is the most frequent silent mistake we see in student reports.

Which fit and RMSD group combination should I use?

gmx rms asks you for two groups, in this order: a group for least-squares fitting, then a group for the RMSD calculation. The table shows what each common combination tells you.

For the general theory of RMSD and RMSF on the protein side, see our guide to analyzing RMSD and RMSF from a GROMACS trajectory. This post covers the ligand-specific steps that guide only mentions in passing.

How do I make an index group for the ligand?

GROMACS usually creates a default group named after the ligand’s residue name (for example JZ4 in Lemkul’s T4 lysozyme system, PDB 3HTB). Hydrogens add noise to the RMSD and many force fields treat them differently, so the convention is to measure heavy atoms only. Build that group with gmx make_ndx:

gmx make_ndx -f em.gro -n index.ndx
 > 13 & ! a H*
 > name 23 JZ4_Heavy
 > q

Here 13 is the number of the ligand group in your own menu and 23 is the number make_ndx assigns to the new group; both will differ on your system, so read the printed list before typing. ! a H* means “not atoms whose name starts with H”. If the ligand has no default group, select it by residue name first with r JZ4 (replace with your residue code). Our gmx make_ndx index groups tutorial covers the selection syntax in more detail.

Tap to explore in 3D

The complexFit on backbone, measure ligandPocket residues around the ligand

Why do I need to fix periodic boundaries before calculating ligand RMSD?

Your simulation runs in a periodic box. When the protein or ligand crosses a box face, GROMACS writes its coordinates wrapped back in on the opposite side. Visually the ligand “teleports”; numerically, its distance from the reference position suddenly becomes several nanometres. Because the ligand is a separate molecule from the protein, it can end up in a different periodic image from the protein even while it sits happily in the pocket.

The fix is to re-center the trajectory on the protein and put whole molecules back in the box, exactly as Lemkul’s tutorial does:

gmx trjconv -s md_0_10.tpr -f md_0_10.xtc -o md_0_10_center.xtc -center -pbc mol -ur compact

Choose Protein for centering and System for output. -pbc mol puts the centre of mass of each molecule inside the box and requires a run input file passed with -s, according to the trjconv manual. -ur compact places atoms at the closest distance from the box centre, which matters for rhombic dodecahedron boxes.

Do not try to fit and re-wrap in the same call. Lemkul’s tutorial warns that “simultaneous PBC re-wrapping and fitting of the coordinates is mathematically incompatible,” and the trjconv manual itself recommends multiple calls when combining -pbc, -fit, -ur and -center. You do not need a separate fitted trajectory for ligand RMSD anyway, because gmx rms fits each frame itself. For every other conversion option, see our gmx trjconv tutorial.

What is the exact gmx rms command for ligand RMSD?

With the centered trajectory and the heavy-atom group ready, run:

gmx rms -s em.tpr -f md_0_10_center.xtc -n index.ndx -tu ns -o rmsd_jz4.xvg

At the first prompt choose Backbone (the fit group). At the second prompt choose JZ4_Heavy (the RMSD group). What each flag does, from the gmx rms manual:

• -s em.tpr: the reference structure. Every frame is compared with these coordinates. Using the energy-minimized structure means the reference is your starting pose. A .tpr also supplies correct masses.
• -f: the re-centered trajectory, not the raw one.
• -n index.ndx: the index file that holds your custom ligand group.
• -tu ns: writes time in nanoseconds instead of the default picoseconds.
• -fit: defaults to rot+trans (rotation and translation), which is what you want. You do not need to type it.
• -mw: mass weighting for superposition, on by default.

To skip the equilibration phase when you report an average, add -b with a start time (in ps by default; for example -b 5000 skips the first 5 ns). To compare several groups against the same fit in one run, use -ng with the number of RMSD groups.

What does a good ligand RMSD value look like?

As a real reference point, Lemkul’s 10 ns tutorial run of 2-propylphenol (JZ4) in T4 lysozyme L99A/M102Q reports a ligand heavy-atom RMSD of about 0.15 ± 0.06 nm (1.5 ± 0.6 Å), which he describes as only a very small change in the ligand’s position. Remember that GROMACS reports RMSD in nanometres; multiply by 10 to get Ångström, the unit most docking papers use.

There is no universal pass mark, but the 2 Å cut-off widely used to judge whether a docking pose reproduces a crystal pose is a sensible yardstick for “the pose was preserved”. Read the shape of the trace, not just the mean:

• Stable: a quick rise in the first nanosecond or so, then a flat plateau with small fluctuations. The pose holds.
• Rearranging: a plateau, a step up, then a new plateau. The ligand found a different binding mode inside the pocket. Look at frames before and after the step in VMD or PyMOL before you call it good or bad.
• Drifting: a steady climb that never levels off. The pose is not stable on this timescale, or the run is too short to tell.
• Unbinding: a climb to large values (often over 1 nm) that never comes back. Confirm with the ligand-to-pocket distance in gmx distance, and by watching the trajectory.

One run is one sample. Before you draw conclusions about stability, read our guide on how many replicate MD simulations you need. To turn the .xvg into a clean figure, use our guide to plotting GROMACS xvg files in Python.

How do I get the average ligand RMSD for my report?

Use gmx analyze on the output file. It prints the average and standard deviation of each data set:

gmx analyze -f rmsd_jz4.xvg

Report the mean and standard deviation over the equilibrated part only (use -b in gmx rms as above), state the fit group and RMSD group explicitly, and say which structure was the reference. Our guide to writing the MD methods section shows how to phrase this.

Why does my ligand RMSD jump or look wrong? (Troubleshooting)

Ligand RMSD of 3 to 6 nm, or sudden vertical spikes. You ran gmx rms on the raw trajectory, so the ligand is crossing periodic boundaries. Re-run the trjconv -pbc mol -center step and use the centered file. If spikes remain in a long run, Lemkul notes that complexes can be hard to center and suggests a custom index group for centering (for example, pocket residues, or a merged protein plus ligand group made in make_ndx with 1 | 13).

Flat, low RMSD even though the ligand visibly left the pocket. You chose the ligand as the fit group. Re-run with Backbone at the first prompt.

Fatal error saying an index “is larger than the number of atoms in the trajectory file”. Your trajectory contains fewer atoms than your index file or .tpr expects, usually because you saved only Protein at the trjconv output prompt. Re-run trjconv and choose System for output, or make the index from a structure with the same atoms as the trajectory.

Warning that masses will be guessed based on residue and atom names. You passed a .gro or .pdb to -s. The gmx rms manual notes that guessed masses are “fine for proteins, but not necessarily for other molecules”. For a ligand, use the .tpr.

A sharp one-frame spike, then back to normal. Often a symmetric group (a phenyl ring or a carboxylate) flipping. gmx rms matches atoms by index, not by chemical symmetry, so a flip counts as deviation even though nothing changed physically. Check the frame visually before reporting it.

The ligand group is missing from the menu. Your ligand residue name was not recognized as a default group. Create it with r LIG in make_ndx (using your residue name) and pass the index with -n. More fixes for setup errors are in our list of common GROMACS errors.

Where does ligand RMSD fit in a full MD analysis?

Ligand RMSD tells you whether the pose survived. It does not tell you why. Pair it with hydrogen bond and contact analysis (see protein-ligand interaction analysis with ProLIF) and, if you need an energy estimate, MM-PBSA binding free energy. For the full sequence of skills from structure preparation to analysis, follow our GROMACS molecular dynamics learning path and the wider computational biology skills roadmap.

For the underlying method, the canonical worked example is Lemkul’s paper, “From Proteins to Perturbed Hamiltonians: A Suite of Tutorials for the GROMACS-2018 Molecular Simulation Package” (Living Journal of Computational Molecular Science, 2018), together with the online protein-ligand complex analysis step.

Frequently asked questions

Should I calculate ligand RMSD with or without hydrogens?

Without. Heavy-atom RMSD is the convention in docking and MD papers, it matches how crystal poses are compared, and it avoids noise from fast-moving hydrogens. Build the group with ! a H* in gmx make_ndx.

What reference structure should gmx rms use for the ligand?

The starting structure of your simulation, typically the energy-minimized em.tpr or the production .tpr. If you want to compare against the crystal pose instead, the reference must contain exactly the same atoms in the same order as the trajectory group.

Is ligand RMSD in GROMACS reported in nm or Å?

Nanometres. Multiply by 10 to convert to Ångström. A value of 0.2 nm is 2 Å.

Can I calculate protein and ligand RMSD in one command?

Yes. Use -ng 2 in gmx rms, choose Backbone as the fit group, then select Backbone and your ligand group as the two RMSD groups. Both columns are written to the same .xvg file.

My ligand RMSD rises to 0.5 nm and stays flat. Is the simulation bad?

Not necessarily. A stable plateau means the ligand settled into a pose, but that pose differs from the starting one by about 5 Å. Compare the two poses visually and check whether the key contacts from docking are still present. It may be the docking pose that was wrong.