How to Analyze RMSD and RMSF from a GROMACS MD Trajectory (Step-by-Step)

To analyze a GROMACS trajectory, run gmx rms to compute backbone RMSD (how far the structure drifts from its starting pose over time) and gmx rmsf to compute per-residue RMSF (which residues move the most). A flat RMSD plateau means your simulation has equilibrated; tall RMSF peaks mark flexible loops and chain termini.

You finished your first molecular dynamics run and now have a trajectory file sitting on disk. A trajectory by itself is not a result. The two analyses every reviewer and thesis committee expects first are RMSD and RMSF, because together they answer the question “did this simulation actually behave like a stable protein?” This guide shows the exact GROMACS commands, the group selections that trip people up, and how to read the output, using the classic lysozyme-in-water system as the worked example.

This is a spoke in our learn molecular dynamics with GROMACS series. If you are mapping out where MD fits in a computational biology career, see the full computational biology skills roadmap.

What do RMSD and RMSF actually measure in an MD trajectory?

RMSD, the root mean square deviation, measures how much the whole structure has moved away from a reference frame, averaged over all selected atoms, and reported for every saved time point. It is a single number per frame, so plotted against time it tells you whether the structure has settled down or is still drifting.

RMSF, the root mean square fluctuation, flips the calculation. Instead of one value per frame, it gives one value per atom or per residue, averaged over the whole trajectory. It tells you where the motion is: a rigid helix core sits low, while flexible loops and the N and C termini spike upward.

GROMACS reports both quantities in nanometers, not angstroms. One nm equals 10 angstroms, so an RMSD of 0.2 nm is 2 angstroms. Keep that in mind when you compare your numbers against papers that quote angstroms.

How do you prepare a GROMACS trajectory before analysis?

Before any RMSD or RMSF calculation, fix periodic boundary conditions. During the run, the protein can drift across the edge of the simulation box and reappear on the opposite side, which creates huge artificial jumps in RMSD that have nothing to do with real motion. Correct this once with gmx trjconv:

gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_nopbc.xtc -pbc mol -center

When prompted, choose Protein as the group to center, then System as the group to output. The -pbc mol flag keeps each molecule whole, and -center recenters the protein in the box. From here on, run every analysis on the cleaned md_0_1_nopbc.xtc file, not the raw trajectory. If jumps still appear, rerun with -pbc nojump instead, which prevents atoms from ever crossing a boundary.

How do you calculate RMSD in GROMACS with gmx rms?

The gmx rms tool does a least-squares fit to remove overall translation and rotation, then computes the deviation. The canonical command for the lysozyme tutorial is:

gmx rms -s md_0_1.tpr -f md_0_1_nopbc.xtc -o rmsd.xvg -tu ns

You will be asked to pick a group twice. The first selection is the group used for least-squares fitting; the second is the group whose RMSD is calculated. For a standard protein, choose Backbone for both. Backbone (N, C-alpha, C atoms) is the standard choice because it tracks the fold without being dominated by wagging side chains. The -tu ns flag sets the time axis to nanoseconds so your plot reads cleanly.

The reference frame here comes from the structure inside md_0_1.tpr, which is the start of the production run. To instead measure deviation from the energy-minimized crystal structure, point the -s flag at the minimization run input:

gmx rms -s em.tpr -f md_0_1_nopbc.xtc -o rmsd_xtal.xvg -tu ns

Comparing both curves is a common reviewer request: deviation from the start shows internal stability, while deviation from the crystal shows how far the model wandered from the experimental structure. The official flag list lives in the GROMACS gmx rms documentation.

How do you read an RMSD plot?

Open rmsd.xvg and look at the shape, not just the final value. A healthy trajectory rises during the first part of the run as the protein relaxes from its starting coordinates, then flattens into a plateau and stays there. That plateau is the signature of equilibration: the structure has found a stable basin and is fluctuating around it rather than drifting.

A curve that keeps climbing for the entire run has not equilibrated. Either the protein is genuinely unfolding, or your production run is too short and needs to be extended. A curve with a sudden vertical step almost always means a periodic boundary artifact slipped through, so go back and rerun gmx trjconv. As a rough rule of thumb, a stable, well-folded globular protein often plateaus at a backbone RMSD on the order of 0.1 to 0.3 nm, though the right value depends on protein size and flexibility, so treat it as orientation rather than a hard pass mark.

How do you calculate per-residue RMSF with gmx rmsf?

RMSF answers a different question: which parts of the protein move. Use the -res flag to average the fluctuation per residue, which produces a far more readable plot than raw per-atom output:

gmx rmsf -s md_0_1.tpr -f md_0_1_nopbc.xtc -o rmsf.xvg -res

When prompted, choose Protein (or C-alpha for a cleaner backbone view). The output rmsf.xvg plots residue number on the x-axis against fluctuation in nm on the y-axis. You can also write the fluctuation values straight into the B-factor column of a PDB file, which lets you color the structure by flexibility in a viewer like PyMOL or VMD:

gmx rmsf -s md_0_1.tpr -f md_0_1_nopbc.xtc -oq bfactor.pdb -res

The full option set is in the GROMACS gmx rmsf documentation.

How do you read an RMSF plot?

Read RMSF as a map of the protein’s mechanics. Low, flat regions correspond to secondary structure, helices and beta sheets that hold their shape. Tall peaks correspond to loops, linkers, and the flexible N and C termini, which almost always show the highest values because their ends are unconstrained. If you docked a ligand or are studying an active site, check the RMSF of those residues specifically: a rigid binding pocket and a flexible recognition loop tell a real structural story you can write up.

RMSD vs RMSF: which one answers which question?

How do you plot the .xvg files GROMACS produces?

GROMACS writes results in .xvg format, designed for the Grace plotting program. The fastest preview is:

xmgrace rmsd.xvg

For a thesis figure you control, load the file in Python instead. Every line starting with # or @ is Grace metadata, so skip it:

import numpy as np
import matplotlib.pyplot as plt

data = np.loadtxt("rmsd.xvg", comments=["#", "@"])
plt.plot(data[:, 0], data[:, 1])
plt.xlabel("Time (ns)")
plt.ylabel("RMSD (nm)")
plt.savefig("rmsd.png", dpi=300)

The same snippet works for rmsf.xvg by relabeling the axes to residue number and RMSF.

Troubleshooting common RMSD and RMSF errors

• RMSD shows a huge sudden jump. A molecule crossed the periodic boundary. Rerun gmx trjconv with -pbc nojump, or -pbc mol -center, before analysis.
• “Index group Backbone not found.” Your system has non-standard residues, so the default protein groups were not built. Create a custom group with gmx make_ndx -f md_0_1.tpr -o index.ndx and pass it with -n index.ndx.
• “Selected group does not contain the same number of atoms.” The -s reference and the -f trajectory do not match. Make sure the .tpr you pass was used to generate that exact .xtc.
• Python cannot read the .xvg. You forgot to skip the metadata. Pass comments=["#", "@"] to np.loadtxt so the header lines are ignored.
• RMSD never plateaus. The system has not equilibrated. Extend the production run, and confirm your NVT and NPT equilibration steps actually converged before you trust the production trajectory.

Frequently asked questions

What is a good RMSD value for a protein simulation?

There is no universal cutoff. What matters is that the curve plateaus and stays flat rather than climbing. For a small, well-folded globular protein, backbone RMSD often settles in the region of 0.1 to 0.3 nm, but larger or more flexible proteins legitimately sit higher. Judge stability by the shape of the curve, not a single threshold.

Should I use Backbone or C-alpha for RMSD?

Backbone (N, C-alpha, C) is the standard choice and is slightly more stringent than C-alpha alone. Either is defensible as long as you state which you used. Avoid fitting on all atoms, because side-chain motion inflates the deviation and hides the fold’s real behavior.

Why is RMSF highest at the ends of my protein?

The N and C termini are tethered on only one side, so they fluctuate more than residues buried in secondary structure. High terminal RMSF is normal and expected; it is not a sign of a broken simulation.

Do I need to correct PBC before every analysis?

Yes. Run gmx trjconv with -pbc mol -center once and save a cleaned trajectory, then use that file for RMSD, RMSF, radius of gyration, and any other analysis. Skipping this step is the most common cause of nonsense plots.

Can I compute RMSD against the experimental crystal structure?

Yes. Point gmx rms -s at the energy-minimized input (for example em.tpr) instead of the production input. This measures how far your model drifted from the starting experimental coordinates rather than from the first production frame.

Where this fits in your MD workflow

RMSD and RMSF are the first two analyses, not the last. Once you trust that your trajectory equilibrated, you move on to radius of gyration (gmx gyrate) for compactness, hydrogen-bond counts, and, for drug-design work, ligand stability inside the binding pocket. Getting RMSD and RMSF right first means everything downstream rests on a trajectory you have actually checked.

This walkthrough follows Justin A. Lemkul’s widely used lysozyme tutorial, published as part of “From Proteins to Perturbed Hamiltonians: A Suite of Tutorials for the GROMACS-2018 Molecular Simulation Package” in the Living Journal of Computational Molecular Science (DOI: 10.33011/livecoms.1.1.5068). For the engine itself, the reference is Abraham et al., “GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers,” SoftwareX 1-2 (2015) 19-25.

Written by the StemSkills Lab team, with 10+ years in sequence and structural bioinformatics, drug discovery and design, and multiscale molecular modeling.