Ramachandran Plot from a GROMACS Trajectory (gmx rama)

Quick answer: To make a Ramachandran plot from a GROMACS trajectory, run gmx rama -s md.tpr -f md_noPBC.xtc -o rama.xvg. It reads the backbone dihedrals from the run input file and writes one phi, psi and residue-label line per residue for every frame. Plot all lines as a 2D density to see the conformations your protein visited during the run.

A Ramachandran plot of a single model tells you whether that model’s backbone is plausible. A Ramachandran plot of a trajectory tells you something different: whether the backbone stays plausible while the protein moves, which residues explore new regions, and whether a loop you modelled by hand settles or keeps flipping. This guide from the StemSkills Lab team (10+ years in structural bioinformatics, drug design and molecular modeling) walks through the full workflow on a GROMACS run, from the trajectory you already have to a publication-quality density plot and per-residue time series.

The worked example uses hen egg-white lysozyme (PDB 1AKI), the system from Justin Lemkul’s widely used GROMACS lysozyme tutorial. Every command works the same on your own protein.

What does a Ramachandran plot show for an MD trajectory?

It shows the distribution of backbone dihedral angles phi (rotation about N to CA) and psi (rotation about CA to C) for every residue across every frame you analyse. The idea goes back to G. N. Ramachandran, C. Ramakrishnan and V. Sasisekharan, who showed in 1963 that steric clashes restrict which phi/psi combinations a peptide can adopt (J. Mol. Biol. 7:95, 1963).

For a static structure each residue gives one point. For a trajectory each residue gives a cloud of points, one per frame. The shape of that cloud carries information a single model cannot:

• A tight cloud inside the alpha-helix or beta-sheet region means stable secondary structure.
• A cloud that spreads between two regions means the residue switches conformation, common in loops and termini.
• A cloud that sits in a disallowed region for most of the run points to a modelling or topology problem, not thermal motion.

This is the natural follow-up to validating a protein structure with MolProbity: you checked the starting model, now you check whether the simulation keeps it healthy.

What inputs does gmx rama need?

Two files: a run input file (.tpr) and a trajectory. According to the official gmx rama manual page, the tool selects the phi/psi combinations from your topology and “computes these as a function of time.”

• -s (default topol.tpr): the portable run input file. It supplies atom names and residue information, so a plain .gro will not do here.
• -f (default traj.xtc): the trajectory. Accepted formats include xtc, trr, cpt, gro, g96, pdb, tng and h5md.
• -o (default rama.xvg): the output file.
• -b, -e, -dt: first frame, last frame and frame stride, all in ps by default.
• -xvg: header formatting, one of xmgrace (default), xmgr or none.

Unlike most gmx analysis tools, gmx rama does not ask you to pick an index group. It finds every backbone C, N, CA, C, N pattern in the topology and computes phi and psi for each residue that completes it.

Most groups process the raw trajectory once with gmx trjconv and reuse that file for every analysis. If you have not done this yet, the standard command is:

gmx trjconv -s md.tpr -f md.xtc -o md_noPBC.xtc -pbc mol -center

Choose Protein for centering and System for output. Our gmx trjconv tutorial explains each option. Strictly, gmx rama already makes molecules whole internally using the .tpr (the source code calls the GROMACS remove-PBC routine on every frame), so periodic jumps rarely corrupt phi and psi. Using the same processed trajectory across all your analyses still keeps frame counts and time stamps consistent.

How do I run gmx rama on my trajectory?

The basic command:

gmx rama -s md.tpr -f md_noPBC.xtc -o rama.xvg

For long runs, thin the frames with -dt. This analyses one frame every 100 ps:

gmx rama -s md.tpr -f md_noPBC.xtc -o rama.xvg -dt 100

To analyse only the equilibrated part, skip the first 10 ns (10000 ps):

gmx rama -s md.tpr -f md_noPBC.xtc -o rama_prod.xvg -b 10000 -dt 100

To get the starting model’s plot for comparison, point -f at a single structure such as the energy-minimised coordinates:

gmx rama -s md.tpr -f em.gro -o rama_start.xvg

What is inside the rama.xvg file?

After the # comment and @ formatting lines, every data line has three fields. In the GROMACS source (gmx_rama.cpp) each line is written as phi, psi and a label:

-62.1  -41.7  ALA-10
-70.4  -35.2  ALA-11
...

(The numbers above illustrate the layout only.) Three things to know:

1. There is no time column. Lines for all residues of frame 1 come first, then all residues of frame 2, and so on. Frame index equals line number divided by the number of residues.
2. Angles are in degrees, from -180 to 180.
3. The label is RESNAME-NUMBER, built from the residue name and the residue number stored in the topology (the source uses the format %s-%d).

List the residue labels present in your file with:

awk '!/^[#@]/{print $3}' rama.xvg | sort -u | head

Pull out one residue’s frames with a whole-word match, so ALA-1 does not also catch ALA-10:

grep -w "ALA-10" rama.xvg > ala10.dat

How do I plot all frames as a 2D density in Python?

Scatter plots of hundreds of thousands of points turn into a solid blob. A 2D histogram shows where the protein actually spends its time. This script reads rama.xvg and draws a density map with numpy and matplotlib:

import numpy as np
import matplotlib.pyplot as plt

phi, psi, labels = [], [], []
with open("rama.xvg") as fh:
    for line in fh:
        if line.startswith(("#", "@")):
            continue
        a, b, lab = line.split()
        phi.append(float(a))
        psi.append(float(b))
        labels.append(lab)

phi = np.array(phi)
psi = np.array(psi)

fig, ax = plt.subplots(figsize=(6, 6))
h = ax.hist2d(phi, psi, bins=180, range=[[-180, 180], [-180, 180]],
              cmap="viridis", cmin=1)
fig.colorbar(h[3], ax=ax, label="count")
ax.set_xlabel("phi (degrees)")
ax.set_ylabel("psi (degrees)")
ax.set_xticks(range(-180, 181, 60))
ax.set_yticks(range(-180, 181, 60))
ax.set_aspect("equal")
fig.savefig("rama_density.png", dpi=300)

bins=180 gives 2 degree bins, and cmin=1 leaves empty bins blank so the background stays white. For a log colour scale, pass norm=matplotlib.colors.LogNorm() to hist2d; it makes sparsely visited regions visible next to the crowded helix region. More general plotting tips are in our guide to plotting GROMACS xvg files in Python.

How do I read the main regions?

• Lower left (phi near -60, psi near -45): right-handed alpha helix.
• Upper left (phi between about -150 and -60, psi between about 120 and 180): beta strands and polyproline II.
• Upper right (phi positive, psi near 45): left-handed helix, mostly glycine and some asparagine and aspartate.
• Right half generally: sparsely populated for non-glycine residues.

Can I do the same analysis with MDAnalysis?

Yes. The Ramachandran class in MDAnalysis.analysis.dihedrals reads the GROMACS files directly and can draw reference contours for you:

import MDAnalysis as mda
from MDAnalysis.analysis.dihedrals import Ramachandran

u = mda.Universe("md.tpr", "md_noPBC.xtc")
rama = Ramachandran(u.select_atoms("protein")).run(step=10)

print(rama.results.angles.shape)   # (n_frames, n_residues, 2)
ax = rama.plot(ref=True, s=2, alpha=0.2)
ax.figure.savefig("rama_mda.png", dpi=300)

Per the MDAnalysis documentation, results.angles is an array of shape n_frames by n_residues by 2, and ref=True adds a background of allowed and marginally allowed regions. Those reference regions were built from 500 PDB files from the Lovell et al. (2003) set, with cutoffs enclosing 90% and 99% of the data points. The class warns about and drops the first and last residues, since they lack phi or psi. If you already use MDAnalysis for trajectory analysis in Python, this keeps everything in one script. Cite it as Michaud-Agrawal et al., J. Comput. Chem. 32:2319 (2011).

How is a trajectory Ramachandran plot different from static validation?

Static tools judge one model against reference statistics. Trajectory tools describe the ensemble. Use both: validate the starting model first, then check the ensemble.

How do I track one residue that leaves the favoured region?

When the density map shows points in an unusual region, find which residue put them there. Extract it with grep -w as above, then plot phi and psi against frame:

import numpy as np
import matplotlib.pyplot as plt

d = np.loadtxt("ala10.dat", usecols=(0, 1))
dt_ps = 100   # the -dt you used, or your output interval
t_ns = np.arange(len(d)) * dt_ps / 1000

fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(t_ns, d[:, 0], ".", ms=2, label="phi")
ax.plot(t_ns, d[:, 1], ".", ms=2, label="psi")
ax.set_xlabel("time (ns)")
ax.set_ylabel("angle (degrees)")
ax.set_ylim(-180, 180)
ax.legend()
fig.savefig("ala10_timeseries.png", dpi=300)

Two patterns matter. A residue that visits an odd region briefly and returns is thermal motion. A residue that jumps and stays there for the rest of the run has found a new state; check it against your RMSF profile and secondary structure from gmx dssp to see whether a helix end frayed or a loop rearranged.

How many outliers are normal during MD?

There is no agreed threshold for trajectories, so do not borrow a number blindly. The well-known MolProbity goals, more than 98% of residues in favoured regions and fewer than 0.05% outliers, come from high-quality experimental structures and apply to single models. The favoured contours themselves derive from the analysis by Lovell et al. (Proteins 50:437, 2003).

Every MD frame is a snapshot with thermal noise, so individual frames will show occasional points outside the favoured region. A practical, defensible approach:

1. Compare the ensemble density with the starting model’s plot (rama_start.xvg).
2. Look for residues that sit in disallowed regions for a large share of frames, not for a few.
3. Extract representative frames, for example gmx trjconv -s md.tpr -f md_noPBC.xtc -o frame_50ns.pdb -dump 50000, and run them through MolProbity to apply the real contours.

If a homology or AlphaFold model drifts into persistent outliers in a low-confidence loop, that is a hint the region needs refinement. See our guides on reading AlphaFold pLDDT and PAE and refining a homology model.

Why do glycine and proline look different?

Glycine has no side chain, so it reaches regions that are crowded for other residues, including positive phi values. Proline’s ring locks phi near -65 degrees, so its points form a narrow vertical band. Lovell et al. (2003) built separate contours for general residues, glycine, proline and pre-proline for exactly this reason. When you plot a trajectory, split glycine and proline out:

is_gly = np.array([l.startswith("GLY") for l in labels])
is_pro = np.array([l.startswith("PRO") for l in labels])
general = ~(is_gly | is_pro)

Then call hist2d on phi[general], phi[is_gly] and phi[is_pro] (with matching psi) in three panels. Glycine points at positive phi are expected, not a warning.

How do I compare wild type and mutant ensembles?

Run gmx rama on both trajectories with the same -b and -dt, then plot the two density maps side by side with the same bins and colour scale, or subtract the normalised histograms to get a difference map. Focus on residues near the mutation site: a substitution to proline or glycine often changes the local backbone distribution most clearly. The full simulation setup is covered in our wild type versus mutant MD guide.

Troubleshooting: what goes wrong and how do I fix it?

gmx rama refuses my -s file

Cause: you passed a .gro or .pdb to -s. Fix: -s takes a .tpr only. Use the .tpr from gmx grompp for the production run. A .gro is fine for -f.

rama.xvg is huge and Python runs out of memory

Cause: one line per residue per frame. A 100 ns run saved every 10 ps has 10,001 frames; for a protein with about 125 residues that is well over a million lines. Fix: use -dt 100 or larger, restrict to a window with -b/-e, or load with np.loadtxt(..., usecols=(0, 1), comments=("#", "@")) and keep labels only when you need them.

Residue numbers do not match my PDB file

Cause: labels use the residue numbers stored in the topology. If you ran gmx pdb2gmx with -renum, residues were renumbered consecutively. With several chains, the same number can appear once per chain, giving duplicate labels. Fix: map numbers with the processed .gro, and for multi-chain systems rely on line order within each frame rather than the label alone, or use MDAnalysis, which keeps segment information.

My ligand or cofactor is missing from the output

Cause: gmx rama only computes dihedrals for atoms named C, N, CA, C, N in sequence. Ligands, waters, ions and most cofactors never match. Fix: this is expected. For ligand torsions use gmx angle with a dihedral index file. Note that matching is by atom name, so a chain break or missing backbone atom can produce a meaningless angle at the gap; check residues next to breaks.

Points pile up on the edges at plus and minus 180

Cause: angles are periodic, so a residue in the beta region with psi near 180 flips between 179 and -179. Time series show vertical jumps that are not real. Fix: for time series, unwrap with np.degrees(np.unwrap(np.radians(d[:, 1]))), or shift psi with np.where(psi < -120, psi + 360, psi) and adjust the axis.

MDAnalysis raises an error about the selection

Cause: with the default check_protein=True, a selection containing non-protein atoms raises an error, and the docs warn that duplicate or missing atoms (altlocs, incomplete residues) can cause errors. Fix: select "protein" only, and for files with alternate locations select one, for example "protein and not altloc B".

What should I report in a thesis or paper?

State the tool and version, the frames analysed (time window and stride), and show the ensemble density next to the starting model’s plot. Cite GROMACS as Abraham et al., SoftwareX 1-2:19 (2015), plus MDAnalysis if you used it. Name any residues that left the favoured region persistently and explain why. Reviewers trust a short, specific explanation more than a claim that “all residues were in allowed regions”.

Ramachandran analysis is one step in the trajectory-analysis stage of our GROMACS molecular dynamics learning path. To see where it fits among all the skills a computational biologist needs, follow the computational biology skills roadmap.

Frequently asked questions

Does gmx rama need an index file?

No. It has no -n option and does not prompt for a group. It computes phi and psi for every residue that has the backbone C, N, CA, C, N pattern in the topology.

Why are the first and last residues missing?

Phi needs the C atom of the previous residue and psi needs the N atom of the next one. Terminal residues lack one of them, so they have no complete phi/psi pair (unless capping groups supply the missing atoms).

Can I make a Ramachandran plot of a single GROMACS structure?

Yes. Pass a single-frame file to -f, for example gmx rama -s md.tpr -f em.gro -o rama_start.xvg. For formal validation of one model, MolProbity is the better choice because it applies curated contours.

Should I analyse the whole trajectory or only production?

Analyse the equilibrated part. Use -b to skip the initial stretch where RMSD is still rising, so early relaxation does not blur the ensemble.

Is gmx rama or MDAnalysis better?

They compute the same angles. gmx rama is a one-line command that writes plain text; MDAnalysis gives a numpy array and built-in reference contours. Pick the one that fits the rest of your analysis pipeline.