How to Plot GROMACS Results in Python from .xvg Files (RMSD, RMSF and Energy)

To plot GROMACS results in Python, read each .xvg file with a small parser that skips the header lines beginning with @ and #, load the remaining columns into NumPy, then draw them with Matplotlib. The same three-step pattern turns rmsd.xvg, rmsf.xvg and energy.xvg into publication-ready figures without ever opening xmgrace.

Every GROMACS analysis command writes its numbers to an .xvg file: gmx rms for RMSD, gmx rmsf for per-residue fluctuation, gmx energy for energy terms, gmx gyrate for radius of gyration. The files are plain text, but beginners get stuck because the tools were built to hand these files to Grace (xmgrace), which many students never install. The good news: an .xvg file is just a header plus a table of numbers, and Python reads it in a few lines. This guide gives you one reusable parser and the exact plotting code for the three figures almost every thesis needs.

What is an .xvg file and why can’t you just load it?

An .xvg file is the Grace data format that GROMACS analysis tools write by default. It has two kinds of header lines you must skip before the numbers start:

• Lines beginning with # are comments (the command that generated the file, the GROMACS version).
• Lines beginning with @ are Grace formatting instructions (the title, axis labels, and the legend for each data series).

Everything after those headers is a whitespace-separated table: the first column is almost always the x value (time, or residue number), and the remaining columns are the y values. If you point a naive loader at the raw file, it hits the first @ line and fails with an error like could not convert string to float: '@'. The fix is simply to tell the loader that @ and # mark lines to ignore. The format itself is described in the official GROMACS analysis documentation.

How do you read an .xvg file in Python?

You have two solid options. For a quick look, NumPy’s loadtxt can skip both comment characters in one call:

import numpy as np
data = np.loadtxt("rmsd.xvg", comments=["@", "#"])

The comments argument accepts a list of strings, so both header types are stripped and you get a clean 2D array. This is documented behaviour of numpy.loadtxt. It works, but it throws away the axis labels GROMACS already wrote for you. A slightly longer parser keeps those labels so your plots are self-documenting:

import numpy as np

def read_xvg(path):
    """Read a GROMACS .xvg file. Returns (data, meta)."""
    meta = {"title": "", "xaxis": "", "yaxis": "", "legends": []}
    rows = []
    with open(path) as fh:
        for line in fh:
            if line.startswith("#"):
                continue
            if line.startswith("@"):
                if "xaxis" in line and '"' in line:
                    meta["xaxis"] = line.split('"')[1]
                elif "yaxis" in line and '"' in line:
                    meta["yaxis"] = line.split('"')[1]
                elif line.startswith("@ s") and "legend" in line:
                    meta["legends"].append(line.split('"')[1])
            continue
            rows.append([float(v) for v in line.split()])
    return np.array(rows), meta

This one function handles every .xvg file in the rest of this guide. It returns the numeric block as a NumPy array plus a small dictionary of the labels GROMACS stored in the @ lines.

Which reader should you use?

Note the pandas caveat: pandas.read_csv accepts only a single comment character, so you would strip @ lines separately. That is why the two NumPy routes above are the cleaner default for .xvg data.

How do you plot an RMSD curve from gmx rms?

First generate the file. RMSD measures how far the structure has drifted from a reference over time, and it is the standard equilibration check. The RMSD and RMSF analysis workflow produces it like this:

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

GROMACS asks for a least-squares fit group and an RMSD group (choose Backbone for both). The -tu ns flag reports time in nanoseconds instead of the default picoseconds. The output has two columns: time in the first, RMSD in nanometres in the second. Now plot it:

import matplotlib.pyplot as plt

data, meta = read_xvg("rmsd.xvg")
time, rmsd = data[:, 0], data[:, 1]

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(time, rmsd, lw=1.2, color="#1f77b4")
ax.set_xlabel(meta["xaxis"] or "Time (ns)")
ax.set_ylabel(meta["yaxis"] or "RMSD (nm)")
ax.set_title("Backbone RMSD")
fig.tight_layout()
fig.savefig("rmsd.png", dpi=300)

Because GROMACS reports RMSD in nanometres, multiply by 10 if a reviewer wants Ångström (1 nm = 10 Å). A curve that rises then flattens is the sign of an equilibrated run; a curve still climbing at the end means you need more sampling.

How do you plot per-residue RMSF from gmx rmsf?

RMSF measures how much each residue wiggles across the trajectory, so its natural x-axis is residue number, not time. Add the -res flag to average per residue:

gmx rmsf -s md.tpr -f md.xtc -o rmsf.xvg -res

The output columns are residue number and RMSF in nanometres. Plotting is the same pattern with a filled profile that makes flexible loops obvious:

data, meta = read_xvg("rmsf.xvg")
res, rmsf = data[:, 0], data[:, 1]

fig, ax = plt.subplots(figsize=(7, 4))
ax.fill_between(res, rmsf, color="#ff7f0e", alpha=0.5)
ax.plot(res, rmsf, color="#d95f0e", lw=1)
ax.set_xlabel("Residue")
ax.set_ylabel("RMSF (nm)")
fig.savefig("rmsf.png", dpi=300)

The tall peaks are your flexible regions (usually loops and termini); the flat valleys are the rigid core. If your x-axis shows atom numbers instead of residues, you forgot the -res flag.

How do you plot energy terms from gmx energy?

gmx energy reads the binary .edr file and is interactive: it lists the available terms and waits for you to type the ones you want. To script it, pipe the selection in with echo:

echo "Potential" | gmx energy -f md.edr -o potential.xvg

The output columns are time in picoseconds and the energy term in kJ/mol. GROMACS reports energies in kJ/mol throughout, which is worth keeping in mind when you interpret a plot: at 300 K the thermal energy RT is about 2.49 kJ/mol (the gas constant 8.314 J/mol/K times 300 K), so fluctuations of a few kJ/mol around a stable mean are normal thermal noise, not drift.

data, meta = read_xvg("potential.xvg")
time, energy = data[:, 0], data[:, 1]

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(time, energy, lw=0.8, color="#2ca02c")
ax.set_xlabel("Time (ps)")
ax.set_ylabel("Potential energy (kJ/mol)")
fig.savefig("potential.png", dpi=300)

If you select several terms at once (for example Potential and Kinetic-En.), the file gains one column per term after the time column, and the meta["legends"] list tells you which is which so you can label each line.

How do you combine RMSD, RMSF and energy into one figure?

A multi-panel figure is what most theses and papers actually print. Matplotlib’s subplots grid handles it cleanly:

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

for ax, fname, ylabel, colour in [
    (axes[0], "rmsd.xvg", "RMSD (nm)", "#1f77b4"),
    (axes[1], "rmsf.xvg", "RMSF (nm)", "#ff7f0e"),
    (axes[2], "potential.xvg", "Energy (kJ/mol)", "#2ca02c"),
]:
    d, m = read_xvg(fname)
    ax.plot(d[:, 0], d[:, 1], lw=1, color=colour)
    ax.set_ylabel(ylabel)
    ax.set_xlabel(m["xaxis"] or "x")

fig.tight_layout()
fig.savefig("summary.png", dpi=300)

The Matplotlib library was built for exactly this, described by its author as a tool for “publication-quality image generation” across platforms (Hunter, Computing in Science & Engineering, 2007). Once this loop works, adding a fourth panel for radius of gyration from gmx gyrate is a one-line change.

Which GROMACS command writes which file?

Keep this table next to your terminal. It maps each analysis command to the file it writes and the columns you will plot.

Troubleshooting: common .xvg plotting errors

• ValueError: could not convert string to float: ‘@’: your loader is reading the header. Use the read_xvg function, or np.loadtxt(..., comments=["@", "#"]).
• IndexError on data[:, 1]: the file has only one column, usually because gmx energy got no term selected. Pipe the term in with echo "Potential" | gmx energy ... and confirm the file has more than one column.
• Time axis runs to huge numbers: the file is in picoseconds. Add -tu ns to gmx rms, or divide the time column by 1000 in Python.
• RMSF x-axis shows atoms, not residues: rerun gmx rmsf with the -res flag.
• Nothing appears and the script hangs on plt.show(): you are on a headless server with no display. Save the figure with fig.savefig("plot.png") instead, or set the backend first with import matplotlib; matplotlib.use("Agg").

For a complete plan from your first simulation to these analysis plots, the GROMACS molecular dynamics track sequences every step, and the wider computational biology skills roadmap shows where MD analysis fits among the other skills grad schools look for. Once these plots feel routine, the next step is turning them into a free energy landscape.

Frequently asked questions

Do I need xmgrace installed to plot GROMACS output?
No. The .xvg file is plain text. Any tool that reads columns of numbers works, and the Python route with NumPy and Matplotlib avoids installing Grace entirely.

Why does GROMACS report RMSD in nanometres and not Ångström?
Nanometres are the standard length unit throughout GROMACS. To convert, multiply by 10, since 1 nm equals 10 Å. Do the conversion in Python before plotting if your field prefers Ångström.

Can I open an .xvg file directly in pandas?
Yes, but pandas.read_csv accepts only one comment character. Strip the @ lines first, or read with the read_xvg function here and pass the resulting array to a DataFrame.

How do I plot several energy terms on the same axis?
Select all of them in gmx energy, then loop over the columns after the time column and use the meta["legends"] entries as labels.

What does the read_xvg meta dictionary give me?
It captures the title, x-axis label, y-axis label, and each series legend that GROMACS wrote in the @ header, so your plots label themselves without hard-coding text.

Written by the StemSkills Lab team, computational scientists with 10+ years in sequence and structural bioinformatics, drug discovery and design, and multiscale molecular modeling. GROMACS analysis routines follow the official manual; NumPy array handling follows Harris et al., Nature, 2020.