How to Analyze Radius of Gyration, SASA and Hydrogen Bonds from a GROMACS Trajectory (Step-by-Step)

After RMSD and RMSF, run gmx gyrate for radius of gyration (how compact and folded the protein stays), gmx sasa for solvent accessible surface area (how much surface is exposed to water), and gmx hbond for hydrogen bonds (how the internal network holds up). A steady radius of gyration and stable hydrogen-bond count both signal a well-behaved, folded trajectory.

You already checked that your simulation equilibrated with RMSD and mapped its flexible regions with RMSF. Those two answer “is the trajectory stable?” but a thesis committee expects more than stability. The next three analyses describe the protein’s shape, its exposure, and its internal glue, and together they turn a raw trajectory into a defensible structural story. This guide gives the exact GROMACS commands, the group selections that trip beginners, and how to read each plot, using the standard lysozyme-in-water system as the worked example.

This is a spoke in our learn molecular dynamics with GROMACS series and follows directly from how to analyze RMSD and RMSF. If you are mapping where MD fits in a research career, see the full computational biology skills roadmap.

Do you still need to correct PBC before these analyses?

Yes. Every analysis below assumes a trajectory with periodic boundary conditions already fixed. If the protein drifted across a box edge during the run, radius of gyration and surface area will show artificial spikes that have nothing to do with real motion. Clean the trajectory once with gmx trjconv, then run everything on the corrected file:

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

Choose Protein as the group to center and System as the group to output. From here on, use md_0_1_nopbc.xtc, not the raw trajectory. GROMACS reports lengths in nanometers and areas in square nanometers, so remember that 1 nm equals 10 angstroms when you compare against papers that quote angstroms.

How do you calculate radius of gyration in GROMACS with gmx gyrate?

The radius of gyration, Rg, is the mass-weighted root mean square distance of a protein’s atoms from its center of mass. In plain terms it is a single number for how compact the structure is: a tightly folded globule has a low, steady Rg, while a protein that swells or unfolds shows Rg climbing over time. Compute it with:

gmx gyrate -s md_0_1.tpr -f md_0_1_nopbc.xtc -o gyrate.xvg

When prompted, choose the Protein group. The output gyrate.xvg has time on the x-axis and five columns: the total Rg followed by its components around the x, y and z axes (Rg-x, Rg-y, Rg-z). The total is what you plot first; the three components tell you whether the molecule is roughly spherical (all three similar) or elongated along one axis (one component clearly larger). The full option list lives in the GROMACS gmx gyrate documentation.

How do you read a radius of gyration plot?

Read the shape before the number. A folded, stable protein produces an Rg trace that stays flat and fluctuates within a narrow band, mirroring the RMSD plateau you already saw. A trace that drifts steadily upward suggests the structure is expanding or partially unfolding, which is a real finding if you are studying thermal or chemical denaturation, and a red flag if you expected a stable fold. For orientation, hen egg white lysozyme (129 residues) sits near an Rg of about 1.4 nm; your own value depends entirely on protein size, so treat published numbers as context, not a pass mark. Compare the Rg trace side by side with your RMSD curve: when both flatten over the same window, that window is your equilibrated, analysis-ready portion of the trajectory.

How do you calculate solvent accessible surface area with gmx sasa?

Solvent accessible surface area (SASA) measures how much of the protein surface a water-sized probe can touch. It reports whether buried hydrophobic residues stay tucked away or become exposed as the structure moves, which is central to folding and stability arguments. Compute the total surface with:

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

When prompted for the group to compute the surface for, choose Protein. The output sasa.xvg plots total SASA in square nanometers against time. To see which residues are exposed rather than just the total, add the per-residue average output:

gmx sasa -s md_0_1.tpr -f md_0_1_nopbc.xtc -o sasa.xvg -or resarea.xvg -tu ns

The resarea.xvg file gives average area per residue, which pairs naturally with your RMSF plot: flexible, high-RMSF loops are often the same residues that gain surface exposure. Under the hood, gmx sasa uses the double cubic lattice method of Eisenhaber and colleagues (Eisenhaber F, Lijnzaad P, Argos P, Sander C, Scharf M, “The double cubic lattice method,” Journal of Computational Chemistry 16 (1995) 273-284, DOI: 10.1002/jcc.540160303). Full flags are in the GROMACS gmx sasa documentation.

How do you count hydrogen bonds with gmx hbond?

Hydrogen bonds are the internal glue of secondary structure, so tracking their number over time tells you whether the fold’s bonding network survives the simulation. As the GROMACS manual puts it, gmx hbond “computes and analyzes hydrogen bonds.” Count them within the protein with:

gmx hbond -s md_0_1.tpr -f md_0_1_nopbc.xtc -num hbnum.xvg

You are prompted to pick two groups: the first supplies donors, the second acceptors. For intramolecular protein bonds, choose Protein for both. To study how the protein interacts with its environment instead, choose Protein and then Water. The output hbnum.xvg plots the number of hydrogen bonds against time. GROMACS decides what counts as a hydrogen bond using two geometric cutoffs: a donor-acceptor distance of 0.35 nm and a hydrogen-donor-acceptor angle of 30 degrees by default. Those defaults are documented in the GROMACS gmx hbond documentation, and you can tighten or loosen them with the -r and -a flags if a reviewer asks.

How do you read a hydrogen-bond count plot?

For intramolecular bonds, look for a roughly constant count that fluctuates around a stable average. A protein does not hold a fixed number every frame; bonds break and reform constantly, so noise around a flat mean is exactly what a healthy fold looks like. A count that trends steadily downward suggests loss of secondary structure, which should agree with a rising RMSD and rising Rg if the protein is genuinely unfolding. Protein-water hydrogen bonds behave differently: they are far more numerous and fluctuate more, because the solvent shell is in constant exchange, so judge them by their average level rather than expecting a tidy plateau.

Radius of gyration vs SASA vs hydrogen bonds: which answers which question?

How do you plot the .xvg files GROMACS produces?

GROMACS writes results in .xvg format for the Grace plotting program, so the quickest preview is xmgrace gyrate.xvg. For a figure you fully control, load the file in Python and skip the metadata lines that start with # or @:

import numpy as np
import matplotlib.pyplot as plt

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

The same snippet works for sasa.xvg and hbnum.xvg by relabeling the axes. This is the same plotting pattern used in the RMSD and RMSF guide, so once you have it saved you can reuse it across every analysis.

Troubleshooting common gmx gyrate, gmx sasa and gmx hbond errors

• Radius of gyration or SASA shows a sudden jump. A molecule crossed the periodic boundary. Rerun gmx trjconv with -pbc mol -center, or -pbc nojump, before analyzing.
• “Index group Protein not found.” Your system has non-standard residues, so the default groups were not built. Create your own with gmx make_ndx -f md_0_1.tpr -o index.ndx and pass it with -n index.ndx.
• gmx sasa complains about a missing selection. Newer versions use the selection framework; if the interactive prompt does not appear, pass the group explicitly, for example -surface 'group "Protein"'.
• gmx hbond count looks far too high. You likely selected System or Water on both prompts. For intramolecular bonds, both selections must be Protein.
• Python cannot read the .xvg. You forgot to skip the header. Pass comments=["#", "@"] to np.loadtxt so Grace metadata is ignored.

Frequently asked questions

What does radius of gyration tell you in an MD simulation?

It reports how compact the protein is at each frame. A stable, folded protein keeps a flat radius of gyration; a steady rise means the structure is expanding or unfolding. Read it alongside RMSD, because a fold that is both low-RMSD and steady-Rg is genuinely stable.

Is a lower SASA always better?

No. SASA is descriptive, not a score. A drop in surface area means the protein is packing more tightly, and a rise means residues are becoming exposed. Which is “good” depends entirely on your question, such as whether you expect folding, unfolding, or a binding event.

Why does my hydrogen-bond count fluctuate so much?

Hydrogen bonds break and reform on a picosecond timescale, so frame-to-frame variation is normal. Judge the trace by its average and whether it stays roughly constant, not by any single frame. A steady downward trend, not noise, is the sign of structural loss.

What distance and angle define a hydrogen bond in GROMACS?

By default gmx hbond uses a donor-acceptor distance cutoff of 0.35 nm and a hydrogen-donor-acceptor angle cutoff of 30 degrees. State these criteria in your methods section, and adjust them with -r and -a only if you have a specific reason.

Which analysis should I run first after RMSD and RMSF?

Radius of gyration, because it is the fastest way to confirm the fold stays compact over the same window where RMSD plateaued. Then add SASA and hydrogen bonds to describe exposure and the bonding network in more detail.

Where this fits in your MD workflow

Radius of gyration, SASA and hydrogen bonds are the second tier of trajectory analysis. RMSD and RMSF tell you the simulation is trustworthy; these three describe what the protein actually did. Once you can defend all five, the next steps are collective-motion analyses such as principal component analysis, and, for drug-design work, tracking ligand contacts inside the binding pocket. Each analysis you add rests on a trajectory you have already checked, which is what makes the final figures defensible.

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.