How to Analyze MD Simulation Results: RMSD, RMSF, Radius of Gyration, Hydrogen Bonds and SASA

Step one

Handle periodic boundaries first: order, group selection and what goes wrong if you skip it

Periodic boundary handling is the first step of trajectory analysis. When a protein crosses the box boundary it is cut into two pieces or jumps to the other side as a whole, and RMSD shows spikes of several nm.

The GROMACS manual recommends this order: make molecules whole (-pbc whole); cluster if needed (-pbc cluster); remove jumps using the first frame as reference (-pbc nojump); center on some group (centering shifts the system, so nojump cannot be used afterwards); if needed, put molecules back in the box with -pbc or -ur; fit last (-fit) and use no PBC-related option after fitting. The trjconv manual also states that -pbc, -fit, -ur and -center cannot always be combined in one call to get the desired result, and that multiple calls should be used.

The fitted trajectory is only for visualization and MDAnalysis. gmx rms, gmx rmsf and gmx covar fit internally. gmx hbond, gmx sasa, gmx rdf and gmx dssp use periodic boundaries for distances, while fitting rotates the coordinates and leaves the box vectors unchanged, so these tools should read the centered, unfitted trajectory. All commands below use md_center.xtc.

A forum case shows the effect of a single step: a GROMACS 2024 DNA system processed only with -pbc mol -ur compact -center gave an RMSD jumping between 0 and 6 nm; after the three steps whole, nojump and center, the RMSD was 0.2 to 0.6 nm.

# Single-chain soluble protein (e.g. lysozyme): two calls are enough
# Call 1: center on Protein, put every molecule back in the box; output System
printf "Protein\nSystem\n" | gmx trjconv -s md.tpr -f md.xtc -o md_center.xtc -pbc mol -center -ur compact
# Call 2 (for visualization and MDAnalysis only): fit on Backbone, output System; no PBC option after this step
printf "Backbone\nSystem\n" | gmx trjconv -s md.tpr -f md_center.xtc -o md_fit.xtc -fit rot+trans

# Multi-chain protein or protein-ligand complex: official order whole -> nojump -> center -> fit
# 0) Create a protein+ligand group (ligand residue name LIG as an example)
printf '"Protein" | "LIG"\nq\n' | gmx make_ndx -f md.tpr -o index.ndx
# 1) Make molecules broken by the box whole
printf "System\n" | gmx trjconv -s md.tpr -f md.xtc -o md_whole.xtc -pbc whole
# 2) Build a reference frame with protein and ligand on the same side (check it in VMD/PyMOL), then use it for nojump
printf "Protein_LIG\nSystem\n" | gmx trjconv -s md.tpr -f md_whole.xtc -n index.ndx -o frame0.gro -dump 0 -pbc mol -center -ur compact
printf "System\n" | gmx trjconv -s frame0.gro -f md_whole.xtc -o md_nojump.xtc -pbc nojump
# 3) Center on protein+ligand; do not use -pbc nojump after centering
printf "Protein_LIG\nSystem\n" | gmx trjconv -s md.tpr -f md_nojump.xtc -n index.ndx -o md_center.xtc -center
# 4) For visualization and MDAnalysis only: fit on the protein backbone; no PBC option after this step
printf "Backbone\nSystem\n" | gmx trjconv -s md.tpr -f md_center.xtc -n index.ndx -o md_fit.xtc -fit rot+trans

# Real dissociation or a display artifact: mindist uses PBC and does not depend on trjconv processing
printf "Protein\nLIG\n" | gmx mindist -s md.tpr -f md.xtc -n index.ndx -od mindist_pl.xvg -tu ns
# Minimum distance between the protein and its periodic image; it should stay above the non-bonded cutoff (rcoulomb, rvdw)
printf "Protein\n" | gmx mindist -s md.tpr -f md.xtc -pi -od mindist_pi.xvg -tu ns

Version differences

Command changes in GROMACS 2023–2026: where copying old tutorials breaks

Most tutorials are based on GROMACS 2018–2022. The table lists commands whose behavior has changed in 2026.

RMSD

RMSD: group selection, ligand RMSD and how to read the curve

RMSD (root mean square deviation) is, for each frame after fitting to the reference structure, the root mean square of the distances between the selected atoms and the corresponding reference atoms. It measures the overall deviation and does not show which part is moving.

gmx rms asks first for the fit group and then for the calculation group. For a protein, choose Backbone or C-alpha both times. For ligand RMSD, choose the protein Backbone as fit group and the ligand as calculation group: the value then includes the ligand's translation and rotation in the pocket as well as its internal conformational change. If the ligand is chosen both times, even a ligand that leaves the pocket gives a small RMSD; reviewers often ask about this.

With md.tpr as reference, the reference is the start of the production run; with em.tpr it is the crystal structure after energy minimization. The difference at the start of the two curves reflects the deviation that already occurred during equilibration.

Read three things from the curve: how long the initial rise lasts; around which value and with what amplitude it fluctuates afterwards; and whether there are step-like jumps. For a step-like jump, rule out a periodic boundary problem first, then decide whether it is a real conformational transition; steps in ligand RMSD often correspond to a change of binding mode or dissociation.

# Backbone RMSD: first prompt = fit group, second prompt = calculation group
printf "Backbone\nBackbone\n" | gmx rms -s md.tpr -f md_center.xtc -o rmsd_bb.xvg -tu ns
# RMSD relative to the crystal structure (em.tpr)
printf "Backbone\nBackbone\n" | gmx rms -s em.tpr -f md_center.xtc -o rmsd_xtal.xvg -tu ns
# Ligand RMSD: fit on protein backbone, compute on ligand, no fit on the ligand itself
printf "Backbone\nLIG\n" | gmx rms -s md.tpr -f md_center.xtc -n index.ndx -o rmsd_lig.xvg -tu ns
# Ligand internal conformational change only (a different quantity; report it separately)
printf "LIG\nLIG\n" | gmx rms -s md.tpr -f md_center.xtc -n index.ndx -o rmsd_lig_self.xvg -tu ns
# All-to-all RMSD matrix, to see whether the system revisits sampled states
printf "Backbone\nBackbone\n" | gmx rms -s md.tpr -f md_center.xtc -m rmsd_matrix.xpm -dt 100

RMSF

RMSF: per-residue fluctuation, and why a peak is not automatically an active site

RMSF (root mean square fluctuation) is the standard deviation of each atom's displacement from its average position over the trajectory. gmx rmsf fits by least squares first by default, and -res gives per-residue averages.

# Drop the first 10 ns (-b in ps because -tu is not set); per-residue C-alpha RMSF
printf "C-alpha\n" | gmx rmsf -s md.tpr -f md_center.xtc -b 10000 -res -o rmsf.xvg -oq bfac.pdb

• Remove the equilibration part (-b) before computing RMSF. Including heating and equilibration raises the RMSF of every residue.
• RMSF peaks usually appear at the N- and C-termini and surface loops. Interpreting a peak directly as an active site has no basis; compare with known functional sites, crystallographic B-factors or the RMSF of the apo system.
• For multi-chain proteins compute each chain separately, or fit on a single chain first; otherwise relative motion between chains is added to every residue.
• When comparing two systems, average RMSF over replicas and plot the interval between replicas. Differences of about 0.05 nm on a single trajectory are often within the replica-to-replica spread.
• -oq converts RMSF into B-factors in a PDB file; color by B-factor in PyMOL and compare side by side with the crystal B-factors.

Radius of gyration and SASA

Radius of gyration and SASA: commands, groups and two common misreadings

The radius of gyration (Rg) is the mass-weighted root mean square distance of the selected atoms from their center of mass and measures how compact the protein is. SASA (solvent accessible surface area) is the area of the surface traced by the center of a 0.14 nm probe sphere rolling over the molecule.

# Radius of gyration (new implementation since GROMACS 2024, selected with -sel)
gmx gyrate -s md.tpr -f md_center.xtc -sel Protein -o gyrate.xvg -tu ns
# SASA: -surface contains all non-solvent atoms, -output takes subsets of it
gmx sasa -s md.tpr -f md_center.xtc -surface Protein -o sasa.xvg -or resarea.xvg -tu ns
# Protein-ligand system
gmx sasa -s md.tpr -f md_center.xtc -surface 'group "Protein" or resname LIG' \
         -output 'group "Protein"' 'resname LIG' -o sasa_pl.xvg -tu ns

Hydrogen bonds

Hydrogen bonds: number, occupancy and lifetime

Since GROMACS 2024, gmx hbond uses a geometric criterion: donor–acceptor distance no more than 0.35 nm (-hbr) and hydrogen–donor–acceptor angle no more than 30° (-hba), with N and O as default donor and acceptor elements.

The new gmx hbond requires the -r and -t selections to be either identical or non-overlapping. The forum has a report of the 2024 tool returning “has no donors AND has no acceptors” for a POPC selection while hbond-legacy worked; GitLab issue 4985 reports that small-molecule carboxylate and sulfonate groups were not recognized as acceptors. If your ligand has such groups, compute again with hbond-legacy or MDAnalysis.

Tested on this page (GROMACS 2026.3): with a ligand topology from ACPYPE, gmx hbond -t 'resname LIG' stops with Selection 'resname LIG' has no donors AND has no acceptors! Nothing to be done., although the ligand has a hydroxyl group and hbond-legacy finds its hydrogen bond to Gln102. The cause is that the [ atomtypes ] written by ACPYPE have no atomic-number column, so the ligand atoms have atomnumber -1 in the tpr, and the new gmx hbond identifies donors and acceptors by element (-de and -ae default to N O). Add the atomic numbers to the atomtypes and rerun grompp; see step 3 of the GROMACS protein-ligand simulation page. Check the atomic numbers in a tpr with gmx dump -s md.tpr | grep atomnumber.

# GROMACS 2024+ new implementation: both -r and -t are required; -num must be given explicitly
gmx hbond -s md.tpr -f md_center.xtc -r Protein -t Protein -num hbnum.xvg -tu ns
# Protein-ligand H-bonds; -o writes the atom indices of each H-bond pair
gmx hbond -s md.tpr -f md_center.xtc -r Protein -t 'resname LIG' -num hbnum_pl.xvg -o hbond_pl.ndx -tu ns
# For lifetimes or an existence matrix use the legacy tool: -ac gives Luzar-Chandler rate constants, -hbm gives per-pair existence per frame
printf "Protein\nLIG\n" | gmx hbond-legacy -s md.tpr -f md_center.xtc -n index.ndx \
         -num hbnum_legacy.xvg -ac hbac.xvg -hbn hbond.ndx -hbm hbmap.xpm

Secondary structure

Secondary structure: gmx dssp replaces do_dssp

gmx dssp assigns secondary structure from the hydrogen-bond pattern between residues and outputs one-letter codes such as H (α-helix), E (β-strand), G (3₁₀-helix), I (π-helix), P (κ-helix, i.e. polyproline II), T, S, B and ~ (loop).

# GROMACS 2023+ has DSSP v4 built in; no external dssp binary or DSSP variable needed
gmx dssp -s md.tpr -f md_center.xtc -sel Protein -o dssp.dat -num dssp_num.xvg -tu ns
# Structures without hydrogens: build pseudo-hydrogens from C and O
gmx dssp -s protein_noH.pdb -f protein_noH.pdb -sel Protein -hmode dssp -clear -o dssp_noH.dat

• -num outputs the number of residues in each secondary structure type per frame, suitable for a stacked plot over time; each line of dssp.dat is one frame and can be plotted in Python as a residue × time map.
• The default -hmode gromacs uses hydrogens present in the structure; for structures without hydrogens use -hmode dssp together with -clear to remove residues missing critical atoms.
• Tested on this page: on the same lysozyme structure with hydrogens removed, gmx dssp with the default -hmode gromacs gives no error but labels nearly all helices as S (bend); with -hmode dssp the normal H and E assignments come back. Forgetting -hmode dssp on a structure without hydrogens produces no warning.
• If results differ from old do_dssp or VMD, first check that the same DSSP version is compared: gmx dssp is equivalent to DSSP v4 by default, and -nopolypro corresponds to DSSP v2.
• A single residue switching from coil to β-bridge in a few frames is a common fluctuation. Only secondary structure changes that recur across replicas and persist for a substantial fraction of the time are suitable as conclusions.

Equilibration and convergence

How to tell whether a simulation has equilibrated and converged: practical checks

Equilibration means the system is no longer influenced by the starting structure; convergence means the quantity of interest is reliably estimated from the available sampling. Neither can be judged from whether the RMSD curve has flattened.

Basis for the number of replicas: the Communications Biology 2023 reliability checklist asks for at least 3 simulations per condition with statistical analysis, and for evidence that results are independent of the initial configuration. Knapp et al. (JCTC 2018) ran 100 replicas for each of two systems; their rule of thumb is a minimum of 5 to 10 replicas, and conclusions from several shorter replicas are more reliable than from a single long trajectory.

Limits of the cosine content: Berk Hess noted on the gmx-users list that convergence of the subspace spanned by the first 10 eigenvectors does not mean the first few eigenvectors themselves have converged. A low cosine content only rules out “close to random diffusion”; it does not prove sufficient sampling.

# Block-averaging error estimate (Hess 2002), production part only
gmx analyze -f gyrate.xvg -b 10 -ee gyrate_errest.xvg
# Cosine content: project one PC at a time, then run gmx analyze -cc
printf "C-alpha\nC-alpha\n" | gmx covar -s md.tpr -f md_center.xtc -b 10000 -o eigenval.xvg -v eigenvec.trr -av average.pdb
printf "C-alpha\nC-alpha\n" | gmx anaeig -s md.tpr -f md_center.xtc -b 10000 -v eigenvec.trr -eig eigenval.xvg -first 1 -last 1 -proj pc1.xvg
gmx analyze -f pc1.xvg -cc pc1_cosine.xvg
# Run covar on each half (production 10-60 ns here), then compare the subspaces spanned by the first 10 eigenvectors
printf "C-alpha\nC-alpha\n" | gmx covar -s md.tpr -f md_center.xtc -b 10000 -e 35000 -o eigval_h1.xvg -v eigvec_h1.trr
printf "C-alpha\nC-alpha\n" | gmx covar -s md.tpr -f md_center.xtc -b 35000 -e 60000 -o eigval_h2.xvg -v eigvec_h2.trr
gmx anaeig -v eigvec_h1.trr -v2 eigvec_h2.trr -first 1 -last 10 -over overlap.xvg

# Extend a finished run (ps; +100 ns here)
gmx convert-tpr -s md.tpr -extend 100000 -o md_ext.tpr
gmx mdrun -deffnm md -s md_ext.tpr -cpi md.cpt   # appends to the existing files by default
# With -noappend you get md.part0002.xtc etc. (the number is the simulation part and increases with each part); concatenate before analysis
gmx trjcat -f md.xtc md.part0002.xtc -o md_all.xtc

PCA and FEL

Principal component analysis and free energy landscapes: gmx covar, anaeig, sham

PCA diagonalizes the covariance matrix of the fitted atomic coordinates; the eigenvectors with the largest eigenvalues (principal components) describe the collective motions with the largest amplitude. A free energy landscape (FEL) converts the 2D distribution P of the trajectory on two principal components into G = −kT ln P.

# 1) Covariance matrix on C-alpha only, equilibration removed (-b in ps)
printf "C-alpha\nC-alpha\n" | gmx covar -s md.tpr -f md_center.xtc -b 10000 -o eigenval.xvg -v eigenvec.trr -av average.pdb
# 2) Project on PC1 and PC2; -2d writes two columns (PC1 PC2) without a time column
printf "C-alpha\nC-alpha\n" | gmx anaeig -s md.tpr -f md_center.xtc -b 10000 -v eigenvec.trr -eig eigenval.xvg \
       -first 1 -last 2 -2d 2dproj.xvg
# 3) Extreme structures along PC1, to see the direction of motion
printf "C-alpha\nC-alpha\n" | gmx anaeig -s md.tpr -f md_center.xtc -v eigenvec.trr -first 1 -last 1 -extr pc1_extreme.pdb -nframes 10
# 4) Free energy landscape: -notime is required; set -tsham to the simulation temperature
gmx sham -f 2dproj.xvg -notime -tsham 300 -ngrid 40 40 40 -nlevels 50 -ls fel.xpm -lp prob.xpm

• Memory and time for gmx covar grow at least with the square of the number of atoms; when memory runs out it ends with a Segmentation fault. Use C-alpha or backbone atoms, not all atoms.
• The two columns written by gmx anaeig -2d are PC1 and PC2, without a time column. gmx sham treats the first column as time by default (-time is on), so without -notime you get a 1D distribution.
• To check that gmx sham treats the input as 2D: the screen output should read Read 2 sets of N points; without -notime it reads Read 1 sets of N points and reports There are 40 bins in the 1-dimensional histogram (tested on this page).
• -tsham defaults to 298.15 K; set it to the simulation temperature. Without -xmin/-xmax the axis ranges are taken from the data. gmx sham outputs kJ/mol with the minimum at 0.
• To compare FELs of two systems (for example with and without ligand), concatenate both trajectories and run covar once, then project each; only then are both maps in the same coordinates. PCA done separately gives different PC1 directions that cannot be compared directly.
• eigenval.xvg gives the eigenvalues; compute the fraction of total variance for PC1 and PC2 and put it in the figure caption.
• Minimum depths in an FEL depend on the amount of sampling. The FEL of a single short trajectory shows where that trajectory stayed, not the equilibrium distribution; an FEL from concatenated replicas is more reliable, and the number of frames used should be reported.
• Getting representative structures from a minimum (an experience-based approach from gmx-users): find the frames whose PC1 and PC2 fall in the minimum bin in 2dproj.xvg, then export them with gmx trjconv -dump. fel_from_2dproj.py below prints the frame indices in the lowest bin; indices start from the anaeig -b start time, and multiplying by the output interval gives the time.
• To build an FEL from two other quantities (for example RMSD and Rg), paste the value columns of two xvg files into a two-column file and pass it to gmx sham -notime; gmx sham only converts a histogram into energies, and you decide which quantities go in.

Other metrics

MSD and RDF

MSD (mean square displacement) is used to compute diffusion coefficients; the RDF (radial distribution function) describes the density of one type of particle around reference particles.

gmx msd fits a straight line on 10% to 90% of the MSD curve by default; if the ends of the curve are not linear, set the range with -beginfit and -endfit. Use a trajectory that has been processed with nojump for MSD (see step 2 of the complex workflow above). gmx rdf uses a bin width of 0.002 nm by default, and -rmax 0 means half the box length.

# Water self-diffusion (fit on 10%-90% of the MSD curve by default)
gmx msd -s md.tpr -f md_nojump.xtc -sel 'resname SOL and name OW' -o msd.xvg
# RDF of water oxygens around the ligand
gmx rdf -s md.tpr -f md_center.xtc -ref 'resname LIG' -sel 'resname SOL and name OW' -selrpos res_com -o rdf.xvg

Python

Ready-to-use Python scripts: multi-replica plots, block averaging, FEL and MDAnalysis

The first script reads the .xvg files of several replicas, plots each replica, the mean and a ±1 SD band, and prints the replica-to-replica standard error and the block-averaged standard error of the production part. The second plots an FEL directly from 2dproj.xvg. The third computes the same metrics with MDAnalysis and MDTraj.

MDAnalysis 2.10.0 fails on a tpr written by GROMACS 2026.3 with Your tpx version is 138, which this parser does not support, yet (tested on this page). Export frame0.pdb first with gmx trjconv -s md.tpr -f md_fit.xtc -o frame0.pdb -dump 0 (output group System) and use it as topology; the MDTraj part of the script also reads this file. A PDB has no charges, so guess_hydrogens does not work; when the script finds no charges it switches to explicit selections by atom name.

count_by_ids() in the script returns atom ids (numbered from 1 in PDB and tpr files), not AtomGroup indices; writing u.atoms[d] directly is off by one atom and prints a neighboring hydrogen as the donor. Tested on this page (MDAnalysis 2.10.0, 50 ps trajectory of the 3HTB complex): after the fix the output is LIG164:O -> GLN102:OE1, consistent with the hydrogen bond between the 2-propylphenol hydroxyl and Gln102 in the crystal structure.

The MDAnalysis RMSF class does not fit, so the trajectory must first be aligned to the average structure. AlignTraj in the script uses in_memory=True and rewrites the coordinates in memory, so RMSF is computed after RMSD and hydrogen bonds.

# Usage: python plot_replicas.py "rep*/rmsd.xvg" 10 RMSD_nm
# Arguments: xvg glob, equilibration time to discard (same unit as the xvg time column), y-axis label
import glob
import sys

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np


def read_xvg(path):
    """Read a GROMACS .xvg: skip # and @ lines, stop at & (first data set only)"""
    rows = []
    with open(path) as fh:
        for line in fh:
            s = line.strip()
            if not s or s[0] in "#@":
                continue
            if s.startswith("&"):
                break
            rows.append([float(v) for v in s.split()])
    return np.array(rows)


def block_se(x, n_blocks):
    """Block-average standard error: split into n_blocks equal blocks, return the SE of the block means"""
    m = len(x) // n_blocks
    blocks = x[: m * n_blocks].reshape(n_blocks, m).mean(axis=1)
    return blocks.std(ddof=1) / np.sqrt(n_blocks)


pattern, t_equil, ylabel = sys.argv[1], float(sys.argv[2]), sys.argv[3]
files = sorted(glob.glob(pattern))
if len(files) < 2:
    sys.exit("need at least 2 replicas, found %d" % len(files))
data = [read_xvg(f) for f in files]
n = min(len(d) for d in data)
t = data[0][:n, 0]
for f, d in zip(files, data):
    if not np.allclose(d[:n, 0], t):
        sys.exit("time axis differs: %s" % f)
Y = np.stack([d[:n, 1] for d in data])  # column 2; total Rg is also column 2
mean, sd = Y.mean(axis=0), Y.std(axis=0, ddof=1)

fig, ax = plt.subplots(figsize=(6, 3.5), dpi=200)
for f, y in zip(files, Y):
    ax.plot(t, y, lw=0.6, alpha=0.5, label=f.split("/")[0])
ax.plot(t, mean, lw=1.5, color="black", label="mean")
ax.fill_between(t, mean - sd, mean + sd, color="gray", alpha=0.3, label="±1 SD")
ax.axvline(t_equil, ls="--", lw=0.8, color="gray")
ax.set_xlabel("Time (same unit as xvg)")
ax.set_ylabel(ylabel)
ax.legend(fontsize=6, frameon=False)
fig.tight_layout()
fig.savefig("replicas.png")

# Production statistics: mean per replica, then SE = SD across replicas / sqrt(n)
prod = Y[:, t >= t_equil]
rep_means = prod.mean(axis=1)
se = rep_means.std(ddof=1) / np.sqrt(len(rep_means))
print("replica means:", np.round(rep_means, 4))
print("mean = %.4f, SE across replicas = %.4f (n=%d)" % (rep_means.mean(), se, len(rep_means)))
# Block averaging: the SE should plateau as blocks grow; no plateau means insufficient sampling
for nb in (40, 20, 10, 5):
    if prod.shape[1] < 2 * nb:  # at least 2 frames per block; skip when there are too few frames
        continue
    print("rep1 blocks=%2d  SE=%.4f" % (nb, block_se(prod[0], nb)))

# Usage: python fel_from_2dproj.py 2dproj.xvg 300
# Same definition as gmx sham: G = -kT ln P, shifted so the minimum is 0
import sys

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np

path, T = sys.argv[1], float(sys.argv[2])
xy = np.array([[float(v) for v in l.split()[:2]] for l in open(path)
               if l.strip() and l.strip()[0] not in "#@&"])
H, xe, ye = np.histogram2d(xy[:, 0], xy[:, 1], bins=40)
P = H / H.sum()
kT = 0.0083144626 * T  # kJ/mol
with np.errstate(divide="ignore"):
    G = -kT * np.log(P)
G -= G[np.isfinite(G)].min()
fig, ax = plt.subplots(figsize=(4.5, 3.8), dpi=200)
im = ax.pcolormesh(xe, ye, np.ma.masked_invalid(G).T, cmap="viridis", shading="flat")
fig.colorbar(im, label="G (kJ/mol)")
ax.set_xlabel("PC1 (nm)")
ax.set_ylabel("PC2 (nm)")
fig.tight_layout()
fig.savefig("fel.png")
print("frames: %d, empty bins: %d / %d" % (len(xy), int((H == 0).sum()), H.size))
# Frame indices in the lowest free-energy bin (counted from the anaeig -b start); export representative structures with gmx trjconv -dump
i, j = np.unravel_index(np.nanargmin(np.where(np.isfinite(G), G, np.nan)), G.shape)
ix = np.clip(np.digitize(xy[:, 0], xe) - 1, 0, len(xe) - 2)
iy = np.clip(np.digitize(xy[:, 1], ye) - 1, 0, len(ye) - 2)
print("frames in the minimum bin:", np.where((ix == i) & (iy == j))[0][:20])

# Usage: python mda_analysis.py frame0.pdb md_fit.xtc
# MDAnalysis 2.10 cannot read GROMACS 2026 tpr files (Your tpx version is 138, which this parser does not support),
# so export frame0.pdb with gmx trjconv -dump 0 as topology; pass md.tpr directly where the tpr can be read
import sys

import MDAnalysis as mda
import mdtraj as md
import numpy as np
from MDAnalysis.analysis import align, rms
from MDAnalysis.analysis.hydrogenbonds import HydrogenBondAnalysis

top, traj = sys.argv[1], sys.argv[2]
u = mda.Universe(top, traj)
protein = u.select_atoms("protein")
has_lig = len(u.select_atoms("resname LIG")) > 0

# 1) RMSD: fit on backbone; ligand RMSD via groupselections, computed after the backbone fit with no extra fit
R = rms.RMSD(u, u, select="backbone",
             groupselections=["resname LIG"] if has_lig else None, ref_frame=0).run()
np.savetxt("rmsd_mda.dat", R.results.rmsd, fmt="%.4f",
           header="frame time_ps backbone_A" + (" ligand_A" if has_lig else ""))

# 2) Radius of gyration (Å, mass-weighted)
rg = np.array([(ts.time, protein.radius_of_gyration()) for ts in u.trajectory])
np.savetxt("rg_mda.dat", rg, fmt="%.4f", header="time_ps rg_A")

# 3) H-bonds: 3.5 Å to stay close to the GROMACS default; state the criteria in the paper
sel = "protein or resname LIG" if has_lig else "protein"
if hasattr(u.atoms, "charges"):  # tpr topology: guess hydrogens and acceptors from charges
    hb = HydrogenBondAnalysis(u, between=["protein", "resname LIG"] if has_lig else None,
                              d_a_cutoff=3.5, d_h_a_angle_cutoff=150)
    hb.hydrogens_sel = hb.guess_hydrogens(sel)
    hb.acceptors_sel = hb.guess_acceptors(sel)
else:  # PDB topology has no charges: select N/O donors, acceptors and H by atom name
    hb = HydrogenBondAnalysis(u, between=["protein", "resname LIG"] if has_lig else None,
                              donors_sel="(%s) and (name N* or name O*)" % sel,
                              hydrogens_sel="(%s) and name H*" % sel,
                              acceptors_sel="(%s) and (name O* or (resname HI* and name ND1 NE2))" % sel,
                              d_a_cutoff=3.5, d_h_a_angle_cutoff=150)
hb.run()
np.savetxt("hbnum_mda.dat", np.column_stack([hb.times, hb.count_by_time()]),
           fmt="%.3f", header="time_ps n_hbonds")
id2ix = {i: k for k, i in enumerate(u.atoms.ids)}  # count_by_ids returns atom ids (1-based), not indices
for d, h, a, n in hb.count_by_ids()[:15]:  # 15 pairs with the highest occupancy
    D, A = u.atoms[id2ix[d]], u.atoms[id2ix[a]]
    print("%s%d:%s -> %s%d:%s  occupancy=%.2f" % (
        D.resname, D.resid, D.name, A.resname, A.resid, A.name, n / u.trajectory.n_frames))
tau, acf = hb.lifetime(tau_max=50)  # continuous autocorrelation; tau in frames
np.savetxt("hb_lifetime_acf.dat", np.column_stack([tau, acf]), fmt="%.4f")

# 4) RMSF: MDAnalysis RMSF does not fit; align to the C-alpha average structure first (in_memory rewrites coordinates in u)
avg = align.AverageStructure(u, u, select="protein and name CA", ref_frame=0).run()
align.AlignTraj(u, avg.results.universe, select="protein and name CA", in_memory=True).run()
ca = u.select_atoms("protein and name CA")
F = rms.RMSF(ca).run()
np.savetxt("rmsf_mda.dat", np.column_stack([ca.resids, F.results.rmsf]), fmt="%.4f",
           header="resid rmsf_A")

# 5) SASA: MDTraj Shrake-Rupley (probe 0.14 nm, result in nm^2)
t = md.load(traj, top="frame0.pdb")
t = t.atom_slice(t.topology.select("protein"))
sasa = md.shrake_rupley(t, probe_radius=0.14, mode="residue")
np.savetxt("sasa_total_nm2.dat", np.column_stack([t.time, sasa.sum(axis=1)]), fmt="%.3f")

Papers

Minimum requirements for reporting these results, and common reviewer questions

The items below come from the Communications Biology 2023 MD reliability and reproducibility checklist and from Grossfield et al.'s best practices for uncertainty. There is no accepted minimum simulation length; it must match the timescale of the process studied.

• At least 3 independent replicas per condition, with a statement of how they were generated (different initial velocities or different starting configurations).
• State how equilibration and production were split, and which part and how many frames were analyzed.
• State the simulation and analysis software and versions, for example GROMACS 2026.3 and the MDAnalysis version; for H-bonds, SASA and similar metrics state the criteria and parameters.
• Provide a system composition table: box dimensions, total number of atoms, ion concentration, protonation states.
• Provide initial coordinates, input files (mdp, top, itp) and final coordinates as supplementary material or in a public repository.
• Give an uncertainty for each mean and state whether it is a standard deviation, a standard error or a 95% confidence interval; with few replicas, plot every replica's value instead of only a mean with error bars.
• Report only significant figures, for example write 1.23456 ± 0.1 as 1.2 ± 0.1.
• Connect simulation results with experimental data where possible, for example crystallographic B-factors, NMR chemical shifts or SAXS curves.

Chinese community

Analysis advice from Chinese sources: what is outdated and what is worth following

The following comes from Sobereva's blog, Jerkwin's blog, answers on the Keinsci forum and CSDN tutorials, checked against the GROMACS 2026 manual and the current versions of the tools.

# DuIvyTools 0.6.0 (install with pip, command name dit)
pip install DuIvyTools
dit xvg_show -f rmsd.xvg                      # view a single curve
dit xvg_compare -f rep1/rmsd.xvg rep2/rmsd.xvg rep3/rmsd.xvg -c 1 1 1   # replicas on one plot
dit xpm_show -f fel.xpm                       # view the free energy landscape from gmx sham
dit xpm2csv -f fel.xpm -o fel.csv             # x, y, z columns for Origin or similar software
dit dssp -f dssp.dat -o dssp.xpm -x "Frame"     # turn dssp.dat from gmx dssp into an old-style xpm and two xvg files

Hand it to an agent

Hand this workflow to an agent

Using the “Lysozyme molecular dynamics simulation” case in Scientify as an example, the analysis task can be described in one sentence.

Example instruction: “Run 3 independent replicas of lysozyme 1AKI at 300 K, 50 ns each; handle periodic boundaries with whole, nojump and centering; compute backbone RMSD, C-alpha RMSF, radius of gyration, SASA, intra-protein hydrogen bonds and DSSP; assess convergence with block averaging and cosine content; plot multi-replica mean ±SD figures and a PC1/PC2 free energy landscape.”

The agent uses the preinstalled GROMACS 2026.3 GPU build on a cloud computer for system setup, simulation and analysis: it generates mdp and tpr files, rents GPUs as needed to run the 3 replicas, executes the trjconv and analysis commands on this page, and summarizes and plots with Python. The workspace keeps the mdp, tpr, logs, xvg files, PNG figures and analysis scripts, so the work can be reproduced. The agent runs an adversarial review of the results, checking that every replica finished and that no jumps remain after PBC handling.

You still need to check yourself: whether the force field and water model suit your question, whether the simulation length covers the process you care about, whether the group selections match what the paper says, and whether the conclusions are supported given the differences between replicas.