How to Analyze a GROMACS Trajectory in Python with MDAnalysis: A Step-by-Step Guide for Beginners

MDAnalysis is a Python library that reads a GROMACS run directly from its .tpr topology and .xtc trajectory, then hands you the coordinates as NumPy arrays. Three lines get you started: import MDAnalysis as mda, build a Universe from the two files, and select the atoms you care about. From there any per-frame analysis you can write in Python is available to you.

Your trajectory is finished, and the standard figures are done. You ran gmx rms and gmx rmsf, you plotted radius of gyration, SASA and hydrogen bonds, and you turned the .xvg files into publication figures in Python. Then your supervisor asks for the distance between your ligand’s carbonyl oxygen and the catalytic lysine, per frame, for all four systems. There is no gmx tool that does exactly that. This is the point where every computational student either learns Python trajectory analysis or spends a week clicking through VMD.

This guide is written by the StemSkills Lab team from 10+ years in structural bioinformatics, drug design and multiscale molecular modeling. Every function name, keyword argument and attribute below is checked against the MDAnalysis documentation and the MDAnalysis User Guide, which documents version 2.9.0. Check your own version with MDAnalysis.__version__ before you copy anything, because two of the details here changed between major releases.

Why would you leave gmx analysis tools at all?

You should not leave them for the standard figures. For RMSD, RMSF, radius of gyration, hydrogen bonds and secondary structure, the compiled gmx tools are already written, already tested and faster to type than any script. Reaching for Python to reproduce gmx rms is a waste of an afternoon.

Python wins in three specific situations, and it is worth being honest that those are the only three.

• The analysis does not exist as a gmx tool. A per-frame distance between one ligand atom and one named residue, a count of water molecules inside a pocket, an angle between two helices defined by your own atom sets. None of these have a dedicated command, and all of them are a short loop.
• You must repeat the same analysis over many systems. Four mutants times three replicas is 12 trajectories. A loop over 12 file paths is reproducible. Twelve rounds of interactive group selection at a gmx prompt is not.
• The number and the figure must come from one script. Reviewers and examiners ask how a value was computed. A single notebook that reads the .xtc, computes the number and draws the plot answers that question in a way a folder of .xvg files never does.

There is a career argument too. “Python and MDAnalysis” on a grad-school application or a CV names a transferable skill. “Ran gmx rms” does not. If you are working out which skills to build in what order, our computational biology skills roadmap puts scripted analysis exactly here, after you can run a simulation and before you attempt free energy work.

How do you install MDAnalysis?

The User Guide gives two supported routes. The conda-forge route creates an isolated environment, which is what you want if GROMACS, AmberTools or a docking stack are already installed and you would rather not disturb them:

mamba create -n mdanalysis -c conda-forge mdanalysis
mamba activate mdanalysis

The pip route installs into whichever environment is currently active:

pip install --upgrade MDAnalysis

One thing catches people out. MDAnalysisTests is a separate package. The User Guide describes it as optional, roughly 90 MB, used for verifying an installation and for running the User Guide examples. You do not need it for your own data, but you do need it the moment you try to run an example that imports something like MDAnalysis.tests.datafiles. Install it the same way:

pip install --upgrade MDAnalysisTests

Confirm what you actually got before you write any analysis:

import MDAnalysis as mda
print(mda.__version__)

Which files do you feed it, and what happens if you only have the .xtc?

MDAnalysis builds a Universe from a topology plus a trajectory. The User Guide draws the distinction plainly: a topology file “provides static information, such as atom identities (name, mass, etc.), charges, and bond connectivity,” while a trajectory file “provides dynamic information, such as coordinates, velocities, forces, and box dimensions.”

For a GROMACS run, that maps cleanly onto the files you already have:

import MDAnalysis as mda

u = mda.Universe("md.tpr", "md_center.xtc")

print(u)
print("frames:", len(u.trajectory))
print("atoms:", u.atoms.n_atoms)

The .tpr is the better topology by a wide margin. The TPRParser documentation lists the attributes it reads: atom names, types, residue names, residue IDs, segment IDs, chain IDs, masses, charges and elements, plus bonds, angles, dihedrals and impropers. That is the full force-field-level description, which is why selections like type OW or anything charge-weighted work at all.

Pass only an .xtc and you get coordinates with no identities. No residue names, no masses, no bonds, so select_atoms("protein") has nothing to match on. Pass only a .gro and you get names and residues but not charges or bond connectivity. The practical rule: use the .tpr as the topology whenever you have it, and fall back to the .gro only when you do not.

That same TPRParser page also documents its supported range: TPR files from GROMACS 4.0 through 2025.0, with the note that “Files generated by the beta versions of Gromacs 2020 are NOT supported.” This range is the source of the most common first error in this workflow, covered in troubleshooting below.

What are the three objects you actually need to understand?

Nearly all beginner confusion comes from missing one of these three ideas.

The Universe is the whole system. It holds every atom and knows how to reach every frame. You usually create exactly one per trajectory.

An AtomGroup is a subset you made with a selection string. u.select_atoms("protein and name CA") returns one. It behaves like a list of atoms with useful attributes: .positions, .resids, .resnames, .masses, .n_atoms.

The trajectory is an iterator, not an array. This is the idea that trips everyone. The User Guide states it directly: “In order to remain memory-efficient, MDAnalysis does not load every frame of your trajectory into memory at once.” Only one frame exists at a time, held in a Timestep, described as “a data structure that holds information for the current time frame.”

The consequence is that an AtomGroup’s .positions silently change meaning as you move through the trajectory:

ca = u.select_atoms("protein and name CA")

u.trajectory[0]
first = ca.positions.copy()      # copy, or you lose it

u.trajectory[-1]
last = ca.positions.copy()

for ts in u.trajectory:
    # ca.positions is now frame ts.frame, time ts.time ps
    pass

Note the .copy(). Without it you are holding a reference to a buffer that the next frame overwrites, and your “first frame” array quietly becomes the last frame you visited. This single mistake produces more wrong figures than any other in Python trajectory analysis.

How does the atom selection language work?

The selection language replaces the interactive group menu you get from gmx make_ndx. It is parsed left to right, and parentheses group terms. Here are the keywords a docking or MD student uses most, all documented on the MDAnalysis selections page.

Booleans and wildcards compose the way you would hope. The documented example protein and not (name H* or type OW) gives heavy protein atoms with water oxygens excluded, and GL[YN] matches both GLY and GLN.

Two habits worth forming immediately. Distances in a selection string are in ångström, so around 3.5 protein is 3.5 Å, not 3.5 nm. And always check the size of a selection before you run an analysis on it, because a selection that matches nothing returns an empty AtomGroup rather than raising an error:

lig = u.select_atoms("resname LIG")
print(len(lig))       # if this prints 0, fix the name before going further

Why are your MDAnalysis numbers ten times your gmx numbers?

This is the single most useful thing on this page, and it costs people whole afternoons.

The MDAnalysis units page states it without ambiguity: “The units of MDAnalysis trajectories are the Å (ångström) for length and ps (picosecond) for time.” Energies are in kJ/mol, forces in kJ/(mol·Å). The same page confirms that whatever the source format uses, MDAnalysis “converts it to MDAnalysis units when reading the data in, and converts back when writing the data out.”

GROMACS works in nanometres. So a gmx rms .xvg file that reports a backbone RMSD of 0.25 will come back from MDAnalysis as 2.5. Both are correct. They are the same quantity in different units.

The failure mode is plotting both on one axis, or comparing a Python number against an .xvg number and concluding your script is broken. It is not broken. Divide MDAnalysis lengths by 10 to get nanometres, or multiply .xvg lengths by 10 to get ångström, and label the axis with the unit you chose. Time needs no conversion: MDAnalysis reports picoseconds and so does GROMACS.

What do you fix before any of this runs?

Periodic boundary artefacts ruin a Python RMSD exactly as they ruin gmx rms. If molecules jump across the box, your RMSD trace shows sudden jumps that have nothing to do with conformational change, and a ligand distance can leap by a box length between consecutive frames.

Fix it once with GROMACS before you open Python, using the same treatment covered in our gmx trjconv guide:

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

Analyse md_center.xtc, not the raw output. Everything below assumes you did this.

How do you reproduce RMSD in MDAnalysis?

Start here, because you can check the answer against a file you already have. The class is MDAnalysis.analysis.rms.RMSD, with this signature:

RMSD(atomgroup, reference=None, select='all', groupselections=None,
     weights=None, weights_groupselections=False, tol_mass=0.1,
     ref_frame=0, **kwargs)

In practice:

import MDAnalysis as mda
from MDAnalysis.analysis import rms

u = mda.Universe("md.tpr", "md_center.xtc")

R = rms.RMSD(u, u, select="backbone", ref_frame=0).run()

data = R.results.rmsd          # N x 3
frame  = data[:, 0]
time_ps = data[:, 1]
rmsd_A  = data[:, 2]

The documentation describes results.rmsd as an “N×3 numpy.ndarray array with content [[frame, time (ps), RMSD (A)], […], …]”. Passing groupselections adds further RMSD columns for extra domains, computed after the superposition on select, which is how you get a whole-protein fit with a per-domain RMSD in the same run.

Watch the results attribute. The older R.rmsd style has been deprecated since MDAnalysis 2.0.0 and is documented for removal in 3.0.0. Tutorials written before 2021 use it, and copying them is how you get a deprecation warning today and a broken script later. Use R.results.rmsd.

Now the check that makes this section worth doing: divide rmsd_A by 10 and it should land on the same curve as your existing gmx rms .xvg for the same selection and the same reference frame. If it does, your Universe, your selection and your units are all correct, and you can trust the next two analyses.

How do you calculate RMSF, and why does alignment come first?

RMSF is where beginners get numbers that look absurd, and the cause is always the same. The MDAnalysis User Guide is explicit: “rms.RMSF does not allow on-the-fly alignment to a reference, and presumes that you have already aligned the trajectory.” An unaligned trajectory contains rigid-body translation and rotation of the whole molecule, and RMSF happily reports that tumbling as if it were residue flexibility.

The guide also notes that “the RMSF should be calculated to an average structure of the simulation,” not to frame 0. That gives the full recipe:

import MDAnalysis as mda
from MDAnalysis.analysis import align, rms

u = mda.Universe("md.tpr", "md_center.xtc")

average = align.AverageStructure(u, u, select="protein and name CA",
                                 ref_frame=0).run()
ref = average.results.universe

aligner = align.AlignTraj(u, ref, select="protein and name CA",
                          in_memory=True).run()

c_alphas = u.select_atoms("protein and name CA")
R = rms.RMSF(c_alphas).run()

for resid, value in zip(c_alphas.resids, R.results.rmsf):
    print(resid, value)          # value in angstrom

The full alignment signature, for when you need to write the aligned trajectory to disk instead of holding it in RAM, is AlignTraj(mobile, reference, select='all', filename=None, prefix='rmsfit_', weights=None, tol_mass=0.1, match_atoms=True, strict=False, force=True, in_memory=False, writer_kwargs=None, **kwargs). Pass a filename and drop in_memory=True when the trajectory is too large to hold in memory.

As with RMSD, R.results.rmsf is the current attribute and the bare R.rmsf alias is deprecated since 2.0.0. Values are in ångström, so divide by 10 before comparing against a gmx rmsf .xvg.

How do you write an analysis that gmx cannot do?

This is the reason the library exists. Suppose you need the distance between a ligand atom and a named binding-site residue, for every frame, so you can show when the contact breaks. There is no command for that. There is a loop.

import numpy as np
import MDAnalysis as mda
from MDAnalysis.lib.distances import calc_bonds

u = mda.Universe("md.tpr", "md_center.xtc")

lig  = u.select_atoms("resname LIG and name O1")
site = u.select_atoms("protein and resid 145 and name NZ")

print(len(lig), len(site))       # both must be 1 before you continue

times, dists = [], []
for ts in u.trajectory:
    d = calc_bonds(lig.positions[0], site.positions[0], box=ts.dimensions)
    times.append(ts.time)
    dists.append(d)

times = np.array(times)          # ps
dists = np.array(dists)          # angstrom

print("mean:", dists.mean(), "A")
print("frames within 3.5 A:", int((dists < 3.5).sum()))

Three details make this correct rather than approximately correct. calc_bonds(coords1, coords2, box=None, result=None, backend='serial') takes the box argument, and passing ts.dimensions is what makes the distance respect periodic boundaries instead of measuring across the box. The len() check catches the case where name O1 does not exist in your ligand, which would otherwise fail deep inside the loop. And the results are appended as scalars, not as coordinate arrays, which keeps memory flat no matter how many frames you have.

Swap calc_bonds for distance_array(reference, configuration, box=None, result=None, backend='serial') when you want every pairwise distance between two groups, for example every ligand heavy atom against every pocket atom. It returns an n × m array per frame.

That same loop shape covers most custom analyses: select once outside the loop, compute a small number inside it, append. It is how you would count pocket waters with around, or track an inter-helix angle, or run the same measurement across the protein-ligand systems in a mutant series.

How do you plot the result?

Once the numbers are NumPy arrays, plotting is ordinary matplotlib, and everything in our guide to plotting GROMACS results in Python applies without change:

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(times / 1000, dists / 10, lw=1)
ax.set_xlabel("Time (ns)")
ax.set_ylabel("O1 to Lys145 NZ distance (nm)")
fig.tight_layout()
fig.savefig("ligand_contact.png", dpi=300)

Note the two conversions in the plot call: picoseconds to nanoseconds, and ångström to nanometres. Do that arithmetic in the plotting line where the axis label sits beside it, and the unit mistake becomes hard to make.

MDAnalysis vs MDTraj vs the gmx tools: which should you use?

Note the units row: MDTraj works in nanometres, MDAnalysis in ångström. Mixing the two libraries in one project is a unit bug waiting to happen. No speed comparison appears in this table on purpose, because we have no benchmark of our own to cite, and an invented one would be worse than none.

Troubleshooting: real errors and their fixes

How do you cite MDAnalysis in your thesis?

The project asks for both papers, not one. The MDAnalysis citations page says: “To properly credit MDAnalysis, please cite both of the following papers.”

• N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein, “MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations,” Journal of Computational Chemistry, 2011, 32, 2319-2327. doi:10.1002/jcc.21787
• R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein, “MDAnalysis: A Python package for the rapid analysis of molecular dynamics simulations,” Proceedings of the 15th Python in Science Conference, 2016, 98-105. doi:10.25080/majora-629e541a-00e

If your comparison table cites MDTraj, the reference is McGibbon and co-authors, “MDTraj: A Modern Open Library for the Analysis of Molecular Dynamics Trajectories,” Biophysical Journal, 2015, 109(8), 1528-1532, doi:10.1016/j.bpj.2015.08.015. State the MDAnalysis version you used alongside the citation, exactly as you would state your GROMACS version. Our guide to the MD methods section shows where these lines belong in the write-up.

Where does this fit in the rest of your GROMACS work?

Python trajectory analysis sits after the standard analyses, not instead of them. Run the usual set first with gmx: RMSD and RMSF, radius of gyration, SASA and hydrogen bonds, secondary structure with gmx dssp, and PCA for essential dynamics. Open Python for the measurement none of those give you, and for anything you need to repeat across systems. The full sequence, from installing GROMACS to writing up, is laid out in our molecular dynamics with GROMACS pillar guide.

Frequently asked questions

Can MDAnalysis read a GROMACS .xtc without a .tpr?

It can read the coordinates, but you get no atom names, residue names, masses or bonds, so selections like protein match nothing. Pass a topology as the first argument to Universe. The .tpr is best because it also carries charges and connectivity; a .gro works if that is all you have.

Why does my MDAnalysis RMSD not match my gmx rms output?

Check units first. MDAnalysis reports ångström and GROMACS .xvg files report nanometres, so a factor of 10 is expected, not a bug. If the shape of the curve also differs, check that you used the same selection and the same reference frame, and that periodic boundaries were fixed with gmx trjconv -pbc mol.

Do I need to align the trajectory before calculating RMSF?

Yes. The MDAnalysis User Guide states that rms.RMSF does not align on the fly and assumes you have already aligned the trajectory. Compute an average structure with align.AverageStructure, align to it with align.AlignTraj, then run rms.RMSF. Skipping this reports whole-molecule tumbling as residue flexibility.

Should I use MDAnalysis or MDTraj?

Use MDAnalysis when your topology is a GROMACS .tpr, when the trajectory is too large to hold in memory, or when you want its selection language. MDTraj is a reasonable choice for array-style work on trajectories that fit in RAM. Do not mix them in one project, because MDTraj uses nanometres and MDAnalysis uses ångström.

How much Python do I need before starting this?

Loops, lists, functions and basic NumPy indexing are enough for everything in this guide. You do not need object-oriented Python. If you can write a for loop that appends to a list and then plot that list with matplotlib, you can write your own trajectory analysis.

Which MDAnalysis version should I cite and use?

Use a current release and report the exact version from MDAnalysis.__version__ in your methods. The User Guide referenced here documents version 2.9.0. Version matters because the results object replaced the older attribute style in 2.0.0, and those older attributes are documented for removal in 3.0.0.