/mdtraj-trajectory-analysis

name: "mdtraj-trajectory-analysis"
description: "mdtraj molecular dynamics trajectory analysis (Python). Reads DCD/XTC/TRR/NetCDF/H5/PDB topologies and trajectories; computes RMSD vs time, radius of gyration, per-residue RMSF, residue-residue contact frequency maps, phi/psi torsions for Ramachandran plots (general + Gly/Pro), and 8-state DSSP secondary structure. Modules: trajectory I/O, geometry (distances/angles/dihedrals), structural analysis (RMSD/Rg/RMSF/SASA), contacts, hydrogen bonds, secondary structure (DSSP), NMR observables. For broader atom-selection grammar use mdanalysis-trajectory; for running MD simulations use OpenMM/GROMACS."
license: "LGPL-2.1"

mdtraj Trajectory Analysis

Overview

mdtraj is a dependency-light Python library for analyzing MD trajectories. Reads DCD/XTC/TRR/NetCDF/H5/AMBER/GROMACS/CHARMM/OpenMM into a `Trajectory` object backed by NumPy arrays, then exposes geometry, RMSD/Rg/RMSF/SASA, contacts, hydrogen bonds, torsions, and 8-state DSSP as pure-Python functions.

> **Units**: mdtraj uses **nm** and **ps** internally. Multiply distances by 10 for Å, divide time by 1000 for ns. Torsions are in **radians** — `np.degrees()`.

When to Use

• RMSD vs time, Rg, per-residue RMSF for stability across MD replicates
• Residue-residue contact frequency maps from a trajectory ensemble
• Backbone phi/psi for Ramachandran (general, Gly, Pro)
• 8-state DSSP per residue per frame for secondary-structure time series
• Lightweight ad-hoc analyses where MDAnalysis's full framework is overkill
• NMR observables (J-couplings, chemical shifts via SHIFTX2)
• Use **mdanalysis-trajectory** instead when you need MDAnalysis's selection grammar, AnalysisBase parallelism, or LAMMPS/NAMD-specific readers
• Use **OpenMM/GROMACS** to *run* the simulation — this skill is post-simulation only

Prerequisites

• `mdtraj`, `numpy`, `pandas`, `matplotlib`
• Topology (PDB/GRO/PRMTOP/PSF) + trajectory (DCD/XTC/TRR/NetCDF/H5)
• Python 3.9+, conda-forge build recommended

Check before installing — inside a pixi/conda env mdtraj is usually present:

python3 -c "import mdtraj" 2>/dev/null || conda install -c conda-forge mdtraj numpy pandas matplotlib

Quick Start

import mdtraj as md

traj = md.load("traj.xtc", top="topology.pdb")
ca = traj.topology.select("name CA")
traj.superpose(traj, frame=0, atom_indices=ca)

rmsd_ang = md.rmsd(traj, traj, frame=0, atom_indices=ca) * 10.0  # nm -> Å
print(f"Frames: {traj.n_frames}  RMSD: {rmsd_ang.min():.2f}–{rmsd_ang.max():.2f} Å")

Core API

Module 1: Trajectory I/O

Load whole or streamed. Format auto-detected from extension.

import mdtraj as md

traj = md.load("rep1.xtc", top="protein.pdb")

# Stream large trajectories — avoids OOM
for chunk in md.iterload("rep1.xtc", top="protein.pdb", chunk=500):
    rmsd_chunk = md.rmsd(chunk, chunk, frame=0)

# Save subset
ca = traj.topology.select("name CA")
traj.atom_slice(ca).save_dcd("ca_only.dcd")

Selecting atoms and slicing frames:

backbone = traj.topology.select("backbone")
chain_a  = traj.topology.select("chainid 0")

first_ns      = traj[:1000]
every_10th    = traj[::10]
last_half_bb  = traj[traj.n_frames // 2:].atom_slice(backbone)

Module 2: RMSD, Rg, RMSF

`md.rmsd` superposes internally; RMSF you compute manually after explicit superpose.

import mdtraj as md, numpy as np

traj = md.load("rep1.xtc", top="protein.pdb")
ca = traj.topology.select("name CA")

rmsd_ang = md.rmsd(traj, traj, frame=0, atom_indices=ca) * 10.0    # Å
rg_ang   = md.compute_rg(traj) * 10.0                              # Å
time_ns  = traj.time / 1000.0
print(f"<RMSD>={rmsd_ang.mean():.2f} Å,  <Rg>={rg_ang.mean():.2f} Å")

# Per-CA RMSF — average-structure reference
ca_traj = traj.atom_slice(ca)
ca_traj.superpose(ca_traj, frame=0)
diff     = ca_traj.xyz - ca_traj.xyz.mean(axis=0)
rmsf_ang = np.sqrt((diff ** 2).sum(axis=2).mean(axis=0)) * 10.0
res_ids  = [a.residue.resSeq for a in ca_traj.topology.atoms]
print(f"Max RMSF: residue {res_ids[np.argmax(rmsf_ang)]} = {rmsf_ang.max():.2f} Å")

Module 3: Contacts and Distance Maps

Threshold distances to get contact frequency.

import mdtraj as md, numpy as np

traj = md.load("rep1.xtc", top="protein.pdb")
distances_nm, pairs = md.compute_contacts(traj, contacts="all", scheme="closest-heavy")
# distances_nm: (n_frames, n_pairs);  pairs: (n_pairs, 2) of residue indices

contact_freq = (distances_nm < 0.5).mean(axis=0)                  # 5 Å cutoff
n_res    = traj.n_residues
freq_map = np.zeros((n_res, n_res))
for (i, j), f in zip(pairs, contact_freq):
    freq_map[i, j] = freq_map[j, i] = f
print(f"Persistent contacts (>0.8): {(contact_freq > 0.8).sum()}")

Module 4: Torsions (Ramachandran)

phi/psi returned in radians; intersect on residue since first residue has no phi and last has no psi.

import mdtraj as md, numpy as np

traj = md.load("rep1.xtc", top="protein.pdb")
phi_ix, phi_rad = md.compute_phi(traj)
psi_ix, psi_rad = md.compute_psi(traj)

def res_of(indices): return np.array([traj.topology.atom(ix[1]).residue.index for ix in indices])
phi_res, psi_res = res_of(phi_ix), res_of(psi_ix)
common  = np.intersect1d(phi_res, psi_res)
phi_deg = np.degrees(phi_rad[:, np.isin(phi_res, common)])
psi_deg = np.degrees(psi_rad[:, np.isin(psi_res, common)])

Filter to Gly / Pro residues:

res_names = [traj.topology.residue(r).name for r in common]
gly_cols  = [i for i, n in enumerate(res_names) if n == "GLY"]
pro_cols  = [i for i, n in enumerate(res_names) if n == "PRO"]
phi_gly, psi_gly = phi_deg[:, gly_cols].ravel(), psi_deg[:, gly_cols].ravel()
phi_pro, psi_pro = phi_deg[:, pro_cols].ravel(), psi_deg[:, pro_cols].ravel()

Module 5: DSSP Secondary Structure (8-state)

`md.compute_dssp(traj, simplified=False)` returns `(n_frames, n_residues)` of one-character codes:

| Code | Meaning | |------|--------------------| |