mdtraj-trajectory-analysis · diff
git:20260528.6b09274 to git:20260528.63fd25c
294 added, 175 removed. Audit A to A.
---
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.
+ mdtraj is a fast, dependency-light Python library for reading, writing, and analyzing molecular dynamics (MD) trajectories. It handles every common MD format (DCD, XTC, TRR, NetCDF, H5, AMBER, GROMACS, CHARMM, OpenMM) into a single `Trajectory` object backed by NumPy arrays, then exposes geometry, RMSD/Rg/RMSF/SASA, contacts, hydrogen bonds, torsions, and 8-state DSSP secondary structure as pure-Python functions. This skill covers the analysis-side workflows: structural stability metrics, contact frequency maps, Ramachandran plots (general + Gly/Pro), and DSSP time series across multi-replicate ensembles.
- > **Units**: mdtraj uses **nm** and **ps** internally. Multiply distances by 10 for Å, divide time by 1000 for ns. Torsions are in **radians** — `np.degrees()`.
+ > **Unit note**: mdtraj works in **nanometers** and **picoseconds** internally. Multiply distances by 10 for Ångström output. Torsions are in **radians** — convert with `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
+ - Computing RMSD vs time, radius of gyration, and per-residue RMSF to assess structural stability across MD replicates
+ - Building residue-residue contact frequency maps from a trajectory ensemble
+ - Extracting backbone phi/psi torsions to plot Ramachandran diagrams (general, glycine-specific, proline-specific)
+ - Assigning 8-state DSSP secondary structure per residue per frame for time-series visualization
+ - Lightweight ad-hoc analyses where MDAnalysis's atom selection grammar or AnalysisBase framework is overkill
+ - NMR observables (J-couplings, chemical shifts via SHIFTX2 wrapper) directly on a trajectory
+ - Use `mdanalysis-trajectory` instead when you need MDAnalysis's richer selection language, AnalysisBase parallelism, or LAMMPS/NAMD-specific readers
+ - For running the MD simulation itself, use OpenMM or GROMACS — this skill is post-simulation analysis only
## Prerequisites
- - `mdtraj`, `numpy`, `pandas`, `matplotlib`
- - Topology (PDB/GRO/PRMTOP/PSF) + trajectory (DCD/XTC/TRR/NetCDF/H5)
- - Python 3.9+, conda-forge build recommended
+ - **Python packages**: `mdtraj`, `numpy`, `pandas`, `matplotlib`, `seaborn`
+ - **Data**: topology (PDB / GRO / PRMTOP / PSF) + trajectory (DCD / XTC / TRR / NetCDF / H5)
+ - **Environment**: Python 3.9+, conda-forge build recommended for binary wheels
- Check before installing — inside a pixi/conda env mdtraj is usually present:
+ Check before installing — inside a pixi/conda env `mdtraj` is usually already present. Inside a pixi project invoke as `pixi run python ...`:
```bash
- python3 -c "import mdtraj" 2>/dev/null || conda install -c conda-forge mdtraj numpy pandas matplotlib
+ command -v python3 && python3 -c "import mdtraj" 2>/dev/null \
+ || conda install -c conda-forge mdtraj numpy pandas matplotlib seaborn
```
## Quick Start
+ Load a trajectory, align to first frame, compute RMSD on alpha-carbons.
+
```python
import mdtraj as md
+ import numpy as np
traj = md.load("traj.xtc", top="topology.pdb")
ca = traj.topology.select("name CA")
- traj.superpose(traj, frame=0, atom_indices=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} Å")
+
+ print(f"Frames: {traj.n_frames} Residues: {traj.n_residues}")
+ print(f"RMSD range: {rmsd_ang.min():.2f}–{rmsd_ang.max():.2f} Å")
```
## Core API
### Module 1: Trajectory I/O
- Load whole or streamed. Format auto-detected from extension.
+ Load and write trajectories. mdtraj auto-detects format from extension.
```python
import mdtraj as md
+ # Load whole trajectory
traj = md.load("rep1.xtc", top="protein.pdb")
- # Stream large trajectories — avoids OOM
+ # Stream large trajectories chunk-by-chunk (avoids OOM)
+ chunk_rmsd = []
for chunk in md.iterload("rep1.xtc", top="protein.pdb", chunk=500):
- rmsd_chunk = md.rmsd(chunk, chunk, frame=0)
+ chunk_rmsd.append(md.rmsd(chunk, chunk, frame=0))
+ print(f"Total chunks: {len(chunk_rmsd)}")
- # Save subset
+ # Save aligned subset
ca = traj.topology.select("name CA")
traj.atom_slice(ca).save_dcd("ca_only.dcd")
```
Selecting atoms and slicing frames:
```python
+ # Topology selection language (subset of VMD/MDAnalysis grammar)
backbone = traj.topology.select("backbone")
+ protein = traj.topology.select("protein")
chain_a = traj.topology.select("chainid 0")
- first_ns = traj[:1000]
+ # Frame slicing
+ first_ns = traj[:1000] # first 1000 frames
every_10th = traj[::10]
- last_half_bb = traj[traj.n_frames // 2:].atom_slice(backbone)
+ last_half_ca = traj[traj.n_frames // 2:].atom_slice(backbone)
+ print(f"Backbone-only, second half: {last_half_ca}")
```
- ### Module 2: RMSD, Rg, RMSF
+ ### Module 2: RMSD, Rg, RMSF (Structural Stability)
- `md.rmsd` superposes internally; RMSF you compute manually after explicit superpose.
+ mdtraj returns distances in nm — multiply by 10 for Å.
```python
- import mdtraj as md, numpy as np
+ import mdtraj as md
+ import numpy as np
traj = md.load("rep1.xtc", top="protein.pdb")
ca = traj.topology.select("name CA")
+ # RMSD vs frame 0 (after internal superposition by md.rmsd)
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} Å")
+ time_ns = traj.time / 1000.0 # ps -> ns
+
+ # Radius of gyration (uses ALL atoms by default; pass masses=None for unweighted)
+ rg_ang = md.compute_rg(traj) * 10.0 # Å
+ print(f"<RMSD>={rmsd_ang.mean():.2f}±{rmsd_ang.std():.2f} Å, <Rg>={rg_ang.mean():.2f} Å")
```
+ RMSF (per-atom or per-residue):
+
```python
- # Per-CA RMSF — average-structure reference
+ # RMSF on CA atoms — 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
+ mean_xyz = ca_traj.xyz.mean(axis=0) # (n_atoms, 3)
+ diff = ca_traj.xyz - mean_xyz
+ rmsf_nm = np.sqrt((diff ** 2).sum(axis=2).mean(axis=0)) # per CA
+ rmsf_ang = rmsf_nm * 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} Å")
+ print(f"Most flexible residue: {res_ids[np.argmax(rmsf_ang)]} RMSF={rmsf_ang.max():.2f} Å")
```
### Module 3: Contacts and Distance Maps
- Threshold distances to get contact frequency.
+ Residue-residue contact distances; threshold to get contact frequency.
```python
- import mdtraj as md, numpy as np
+ import mdtraj as md
+ import 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
+ # Distances between every residue pair (closest heavy atom)
+ distances_nm, residue_pairs = md.compute_contacts(traj, contacts="all", scheme="closest-heavy")
+ # distances_nm: (n_frames, n_pairs); residue_pairs: (n_pairs, 2)
+
+ # Contact frequency: fraction of frames within 0.5 nm (5 Å)
+ cutoff = 0.5
+ contact_freq = (distances_nm < cutoff).mean(axis=0) # (n_pairs,)
+
+ # Reshape to square matrix
+ n_res = traj.n_residues
freq_map = np.zeros((n_res, n_res))
- for (i, j), f in zip(pairs, contact_freq):
+ for (i, j), f in zip(residue_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.
+ Compute phi/psi per residue per frame; mdtraj returns radians.
```python
- import mdtraj as md, numpy as np
+ import mdtraj as md
+ import 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)
+ # md.compute_phi / compute_psi return (atom_indices, angles_rad)
+ phi_indices, phi_rad = md.compute_phi(traj)
+ psi_indices, psi_rad = md.compute_psi(traj)
+
+ # Match phi and psi by residue (phi[i] and psi[i] don't always align; intersect on residue)
+ 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_indices), res_of(psi_indices)
+ 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)])
+ print(f"Residues with both phi and psi: {len(common)}")
```
- Filter to Gly / Pro residues:
+ Filtering Ramachandran samples to Gly / Pro residues:
```python
+ # Residue identity for the columns above
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()
+ print(f"Gly samples: {len(phi_gly)}, Pro samples: {len(phi_pro)}")
```
### Module 5: DSSP Secondary Structure (8-state)
- `md.compute_dssp(traj, simplified=False)` returns `(n_frames, n_residues)` of one-character codes:
+ `md.compute_dssp(traj, simplified=False)` returns an `(n_frames, n_residues)` array of single-character DSSP codes:
| Code | Meaning |
|------|--------------------|
| `H` | alpha-helix |
| `B` | beta-bridge |
| `E` | beta-sheet |
| `G` | 3_10-helix |
| `I` | pi-helix |
| `T` | turn |
| `S` | bend |
| ` ` (space) or `C` | coil |
```python
import mdtraj as md
traj = md.load("rep1.xtc", top="protein.pdb")
- dssp = md.compute_dssp(traj, simplified=False)
+ dssp = md.compute_dssp(traj, simplified=False) # 8-state
+ # Replace whitespace with explicit 'C' for plotting
dssp[dssp == " "] = "C"
print(f"DSSP grid: {dssp.shape} unique codes: {set(dssp.flatten())}")
```
- ### Module 6: Hydrogen Bonds and SASA
+ ### Module 6: Visualization Helpers (matplotlib)
- ```python
- import mdtraj as md
- import pandas as pd
+ mdtraj returns NumPy arrays — plotting is plain matplotlib/seaborn. Convention used across this skill:
- traj = md.load("rep1.xtc", top="protein.pdb")
+ ```python
+ import matplotlib.pyplot as plt
- # Baker–Hubbard: returns (n_hbonds, 3) of [donor, H, acceptor] atom indices
- hbonds = md.baker_hubbard(traj, freq=0.5, periodic=False)
- label = lambda a: f"{traj.topology.atom(a).residue}-{traj.topology.atom(a).name}"
- hbonds_df = pd.DataFrame({"donor": [label(d) for d, _, _ in hbonds],
- "acceptor": [label(a) for _, _, a in hbonds]})
+ # Consistent style for stability metrics
+ plt.rcParams.update({"figure.dpi": 120, "savefig.bbox": "tight",
+ "axes.spines.top": False, "axes.spines.right": False})
- # SASA in nm^2 per residue per frame
- sasa_nm2 = md.shrake_rupley(traj, mode="residue")
- total_ang2 = sasa_nm2.sum(axis=1) * 100.0 # nm^2 -> Å^2
- print(f"Persistent H-bonds: {len(hbonds_df)}, <SASA>: {total_ang2.mean():.0f} Ų")
+ fig, ax = plt.subplots(figsize=(8, 4))
+ ax.plot(time_ns, rmsd_ang, lw=1, label="RMSD")
+ ax.set_xlabel("Time (ns)"); ax.set_ylabel("RMSD (Å)")
+ ax.legend(); fig.savefig("rmsd.png")
```
## Key Concepts
- ### Units (nm, ps, radians)
+ ### Concept 1: Trajectory units (nm, ps, radians)
- `traj.xyz` is in nm, `traj.time` in ps, torsions in radians. Convert before plotting:
+ mdtraj stores `traj.xyz` in **nm**, `traj.time` in **ps**, and torsions in **radians**. The conversion is your responsibility before plotting:
```python
+ xyz_ang = traj.xyz * 10.0
+ time_ns = traj.time / 1000.0
import numpy as np
- xyz_ang = traj.xyz * 10.0
- time_ns = traj.time / 1000.0
- phi_deg = np.degrees(phi_rad)
+ phi_deg = np.degrees(phi_rad)
```
- ### Pair indices vs `resSeq`
+ ### Concept 2: Pair indices vs residue indices
- `md.compute_contacts(..., contacts="all")` returns **0-based residue indices** (contiguous across chains), not PDB `resSeq`. Reshape into `(n_res, n_res)` before heatmaps. Map back via `traj.topology.residue(i).resSeq` only when labeling.
+ `md.compute_contacts(..., contacts="all")` returns `(n_pairs, 2)` of **residue indices** (0-based, contiguous across chains). Reshape into a square `(n_res, n_res)` matrix before plotting heatmaps. Don't confuse with `resSeq` (PDB numbering, may have gaps and chain repeats).
- ### 8-state vs simplified DSSP
+ ### Concept 3: 8-state vs simplified DSSP
- `simplified=False` → 8 states. `simplified=True` → 3 states (H/E/C), which loses 3_10 (G), pi-helix (I), beta-bridge (B), bend (S), turn (T).
+ `simplified=False` → 8 states (H/B/E/G/I/T/S/' '). `simplified=True` → 3 states (H/E/C). The 8-state output is required to distinguish 3_10 (G), pi-helix (I), beta-bridge (B), bend (S), and turn (T) — collapsing to 3 states loses these signals.
## Common Workflows
### Workflow 1: Multi-Replicate Stability + Flexibility + Contacts + Ramachandran + DSSP
- **Goal**: structural-analysis report for three MD replicates — RMSD/Rg over time, per-residue RMSF, residue contact frequency map (rep1), Ramachandran general + Gly + Pro (rep1), DSSP 8-state time series (rep1). Plots use matplotlib defaults; the consumer picks palette downstream.
+ **Goal**: full structural-analysis report for three MD replicates — RMSD/Rg over time, per-residue RMSF, residue contact frequency map (replica 1), Ramachandran (general + Gly + Pro for replica 1), and DSSP 8-state time series (replica 1). Plots use matplotlib defaults so the consumer picks their own palette downstream.
```python
- import mdtraj as md, numpy as np, matplotlib.pyplot as plt
+ import mdtraj as md
+ import numpy as np
+ import matplotlib.pyplot as plt
from pathlib import Path
- replicas = {"rep1": "rep1.xtc", "rep2": "rep2.xtc", "rep3": "rep3.xtc"}
+ # --- Inputs -----------------------------------------------------------------
topology = "protein.pdb"
+ replicas = {
+ "rep1": "rep1.xtc",
+ "rep2": "rep2.xtc",
+ "rep3": "rep3.xtc",
+ }
outdir = Path("figures"); outdir.mkdir(exist_ok=True)
- trajs = {name: md.load(p, top=topology) for name, p in replicas.items()}
+ # --- Load all three replicates ---------------------------------------------
+ trajs = {name: md.load(path, top=topology) for name, path in replicas.items()}
ca = next(iter(trajs.values())).topology.select("name CA")
- # 1. RMSD vs time + Rg (Å) for all replicates
+ # ============================================================================
+ # 1. RMSD vs time + Radius of gyration (Å) for all three replicates
+ # ============================================================================
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
+
for name, traj in trajs.items():
traj.superpose(traj, frame=0, atom_indices=ca)
- ax1.plot(traj.time / 1000.0, md.rmsd(traj, traj, frame=0, atom_indices=ca) * 10.0, lw=1, label=name)
- ax2.plot(traj.time / 1000.0, md.compute_rg(traj) * 10.0, lw=1, label=name)
- ax1.set(xlabel="Time (ns)", ylabel="RMSD (Å)", title="Backbone RMSD")
+ rmsd_ang = md.rmsd(traj, traj, frame=0, atom_indices=ca) * 10.0
+ rg_ang = md.compute_rg(traj) * 10.0
+ t_ns = traj.time / 1000.0
+ ax1.plot(t_ns, rmsd_ang, lw=1, label=name)
+ ax2.plot(t_ns, rg_ang, lw=1, label=name)
+
+ ax1.set(xlabel="Time (ns)", ylabel="RMSD (Å)", title="Backbone RMSD vs reference")
ax2.set(xlabel="Time (ns)", ylabel="Rg (Å)", title="Radius of gyration")
for ax in (ax1, ax2): ax.legend(frameon=False)
- fig.tight_layout(); fig.savefig(outdir / "01_rmsd_rg.png"); plt.close(fig)
+ fig.tight_layout(); fig.savefig(outdir / "01_rmsd_rg.png", dpi=200)
+ plt.close(fig)
- # 2. Per-residue RMSF (Å) across replicates
+ # ============================================================================
+ # 2. Per-residue RMSF (Å) across three replicates
+ # ============================================================================
fig, ax = plt.subplots(figsize=(10, 4))
for name, traj in trajs.items():
- ct = traj.atom_slice(ca); ct.superpose(ct, frame=0)
- rmsf = np.sqrt(((ct.xyz - ct.xyz.mean(axis=0)) ** 2).sum(axis=2).mean(axis=0)) * 10.0
- ax.plot([a.residue.resSeq for a in ct.topology.atoms], rmsf, lw=1, label=name)
+ 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]
+ ax.plot(res_ids, rmsf_ang, lw=1, label=name)
ax.set(xlabel="Residue", ylabel="RMSF (Å)", title="Per-residue flexibility")
ax.legend(frameon=False); fig.tight_layout()
- fig.savefig(outdir / "02_rmsf.png"); plt.close(fig)
+ fig.savefig(outdir / "02_rmsf.png", dpi=200); plt.close(fig)
- # 3. Contact frequency map (rep1, 5 Å cutoff)
+ # ============================================================================
+ # 3. Contact frequency map for replica 1
+ # ============================================================================
rep1 = trajs["rep1"]
dist_nm, pairs = md.compute_contacts(rep1, contacts="all", scheme="closest-heavy")
- freq = (dist_nm < 0.5).mean(axis=0)
- freq_map = np.zeros((rep1.n_residues, rep1.n_residues))
+ cutoff_nm = 0.5 # 5 Å
+ freq = (dist_nm < cutoff_nm).mean(axis=0)
+ n_res = rep1.n_residues
+ freq_map = np.zeros((n_res, n_res))
for (i, j), f in zip(pairs, freq):
freq_map[i, j] = freq_map[j, i] = f
np.fill_diagonal(freq_map, 1.0)
+
fig, ax = plt.subplots(figsize=(6, 5))
im = ax.imshow(freq_map, vmin=0, vmax=1, origin="lower")
- ax.set(xlabel="Residue", ylabel="Residue", title="Contact frequency (rep1)")
+ ax.set(xlabel="Residue", ylabel="Residue", title="Contact frequency (rep1, 5 Å cutoff)")
fig.colorbar(im, ax=ax, label="Fraction of frames in contact")
- fig.tight_layout(); fig.savefig(outdir / "03_contact_map.png"); plt.close(fig)
+ fig.tight_layout(); fig.savefig(outdir / "03_contact_map.png", dpi=200)
+ plt.close(fig)
- # 4. Ramachandran (rep1) — phi/psi 2D density
+ # ============================================================================
+ # 4. Ramachandran (replica 1) — 2D density of phi/psi
+ # ============================================================================
phi_ix, phi_rad = md.compute_phi(rep1)
psi_ix, psi_rad = md.compute_psi(rep1)
- def res_of(ix, t): return np.array([t.topology.atom(i[1]).residue.index for i in ix])
- phi_res, psi_res = res_of(phi_ix, rep1), res_of(psi_ix, rep1)
- common = np.intersect1d(phi_res, psi_res)
+
+ def res_of(indices, traj):
+ return np.array([traj.topology.atom(ix[1]).residue.index for ix in indices])
+
+ phi_res = res_of(phi_ix, rep1); psi_res = res_of(psi_ix, rep1)
+ common = np.intersect1d(phi_res, psi_res)
phi_deg = np.degrees(phi_rad[:, np.isin(phi_res, common)]).ravel()
psi_deg = np.degrees(psi_rad[:, np.isin(psi_res, common)]).ravel()
- H, _, _ = np.histogram2d(phi_deg, psi_deg, bins=72, range=[[-180, 180], [-180, 180]])
+
fig, ax = plt.subplots(figsize=(5, 5))
+ H, xe, ye = np.histogram2d(phi_deg, psi_deg, bins=72, range=[[-180, 180], [-180, 180]])
ax.imshow(H.T, origin="lower", extent=[-180, 180, -180, 180], aspect="equal")
- ax.set(xlabel=r"$\phi$ (°)", ylabel=r"$\psi$ (°)", title="Ramachandran (rep1)")
- fig.tight_layout(); fig.savefig(outdir / "04_ramachandran_all.png"); plt.close(fig)
+ ax.set(xlabel=r"$\phi$ (°)", ylabel=r"$\psi$ (°)",
+ title="Ramachandran (rep1) — density")
+ fig.tight_layout(); fig.savefig(outdir / "04_ramachandran_all.png", dpi=200)
+ plt.close(fig)
+ # ============================================================================
# 5. Glycine + proline Ramachandran
+ # ============================================================================
res_names = [rep1.topology.residue(r).name for r in common]
- for label, cols, fname in [("Glycine", [i for i, n in enumerate(res_names) if n == "GLY"], "05_rama_gly.png"),
- ("Proline", [i for i, n in enumerate(res_names) if n == "PRO"], "06_rama_pro.png")]:
- if not cols: continue
+ 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"]
+
+ for label, cols, fname in [("Glycine", gly_cols, "05_rama_gly.png"),
+ ("Proline", pro_cols, "06_rama_pro.png")]:
+ if not cols:
+ continue
phi = np.degrees(phi_rad[:, cols]).ravel()
psi = np.degrees(psi_rad[:, cols]).ravel()
H, _, _ = np.histogram2d(phi, psi, bins=60, range=[[-180, 180], [-180, 180]])
fig, ax = plt.subplots(figsize=(5, 5))
ax.imshow(H.T, origin="lower", extent=[-180, 180, -180, 180], aspect="equal")
- ax.set(xlabel=r"$\phi$ (°)", ylabel=r"$\psi$ (°)", title=f"{label} Ramachandran (rep1)")
- fig.tight_layout(); fig.savefig(outdir / fname); plt.close(fig)
+ ax.set(xlabel=r"$\phi$ (°)", ylabel=r"$\psi$ (°)",
+ title=f"{label} Ramachandran (rep1, frequency)")
+ fig.tight_layout(); fig.savefig(outdir / fname, dpi=200); plt.close(fig)
- # 6. DSSP 8-state time series (rep1)
- dssp = md.compute_dssp(rep1, simplified=False); dssp[dssp == " "] = "C"
+ # ============================================================================
+ # 6. DSSP 8-state time series for replica 1
+ # ============================================================================
+ dssp = md.compute_dssp(rep1, simplified=False)
+ dssp[dssp == " "] = "C"
+
codes_order = ["C", "E", "B", "S", "T", "H", "I", "G"]
code_to_int = {c: i for i, c in enumerate(codes_order)}
- dssp_int = np.vectorize(lambda c: code_to_int.get(c, 0))(dssp)
+ dssp_int = np.vectorize(lambda c: code_to_int.get(c, 0))(dssp) # (n_frames, n_res)
+
fig, ax = plt.subplots(figsize=(10, 5))
ax.imshow(dssp_int.T, aspect="auto", origin="lower", interpolation="nearest",
extent=[rep1.time[0] / 1000.0, rep1.time[-1] / 1000.0, 0, rep1.n_residues])
- ax.set(xlabel="Time (ns)", ylabel="Residue", title="DSSP 8-state (rep1)")
- fig.tight_layout(); fig.savefig(outdir / "07_dssp_8state.png"); plt.close(fig)
+ ax.set(xlabel="Time (ns)", ylabel="Residue", title="DSSP 8-state evolution (rep1)")
+ fig.tight_layout(); fig.savefig(outdir / "07_dssp_8state.png", dpi=200)
+ plt.close(fig)
- # Dump DSSP code mapping so downstream code can pick its own legend/palette
+ # Optional: dump the DSSP code→index mapping alongside the figure so a downstream
+ # script can pick its own color scheme and legend
np.savez(outdir / "07_dssp_8state.npz", dssp=dssp, dssp_int=dssp_int, codes_order=codes_order)
+
print("All figures written to figures/")
```
- ### Workflow 2: H-Bond Persistence + SASA Across Replicates
+ ### Workflow 2: Hydrogen-Bond Persistence + SASA Time Series
- **Goal**: persistent hydrogen bonds (≥50% of frames) and total SASA time series per replicate.
+ **Goal**: classify hydrogen bonds present in >50% of frames (Baker–Hubbard) and track solvent-accessible surface area over time.
```python
- import mdtraj as md, numpy as np, pandas as pd
+ import mdtraj as md
+ import numpy as np
+ import pandas as pd
- replicas = {"rep1": "rep1.xtc", "rep2": "rep2.xtc", "rep3": "rep3.xtc"}
- top = "protein.pdb"
- rows = []
- for name, p in replicas.items():
- t = md.load(p, top=top)
- hbonds = md.baker_hubbard(t, freq=0.5, periodic=False)
- sasa_ang2 = (md.shrake_rupley(t, mode="residue").sum(axis=1) * 100.0)
- rows.append({"replica": name, "n_persistent_hbonds": len(hbonds),
- "sasa_mean_ang2": sasa_ang2.mean(), "sasa_std_ang2": sasa_ang2.std()})
- print(pd.DataFrame(rows).to_string(index=False))
+ traj = md.load("rep1.xtc", top="protein.pdb")
+
+ # Baker–Hubbard H-bonds: returns (n_hbonds, 3) of [donor, H, acceptor] atom indices
+ hbonds = md.baker_hubbard(traj, freq=0.5, periodic=False)
+ label = lambda a: f"{traj.topology.atom(a).residue}-{traj.topology.atom(a).name}"
+ hbonds_df = pd.DataFrame({
+ "donor": [label(d) for d, _, _ in hbonds],
+ "acceptor": [label(a) for _, _, a in hbonds],
+ })
+ print(f"Persistent H-bonds (>=50% frames): {len(hbonds_df)}")
+
+ # SASA in nm^2 per residue per frame; sum residues for total
+ sasa_nm2 = md.shrake_rupley(traj, mode="residue") # (n_frames, n_residues)
+ total_ang2 = sasa_nm2.sum(axis=1) * 100.0 # nm^2 -> Å^2
+ print(f"<Total SASA> = {total_ang2.mean():.0f} ± {total_ang2.std():.0f} Ų")
```
## Key Parameters
| Parameter | Module | Default | Range / Options | Effect |
|-----------|--------|---------|-----------------|--------|
- | `chunk` (`md.iterload`) | I/O | `100` | `100`-`5000` | Frames per stream chunk; trade speed vs RAM |
- | `atom_indices` (`md.rmsd`) | Stability | `None` | array of ints | Restrict superposition + RMSD (e.g., CA only) |
- | `frame` (`md.rmsd`) | Stability | `0` | `0`-`n_frames-1` | Reference frame |
- | `contacts` (`md.compute_contacts`) | Contacts | `"all"` | `"all"`, `[(i,j),...]` | All residue pairs or custom list |
+ | `chunk` (`md.iterload`) | I/O | `100` | `100`-`5000` | Frames per streamed chunk; raise for speed, lower for RAM |
+ | `atom_indices` (`md.rmsd`) | Stability | `None` (all) | array of ints | Restrict superposition + RMSD to subset (e.g., CA only) |
+ | `frame` (`md.rmsd`) | Stability | `0` | `0`-`n_frames-1` | Reference frame for RMSD |
+ | `contacts` (`md.compute_contacts`) | Contacts | `"all"` | `"all"`, `[(i,j),...]` | Compute every residue pair, or a custom pair list |
| `scheme` (`md.compute_contacts`) | Contacts | `"closest-heavy"` | `"ca"`, `"closest"`, `"closest-heavy"`, `"sidechain"`, `"sidechain-heavy"` | Atom subset for inter-residue distance |
- | `cutoff` (contact threshold, user-set) | Contacts | `0.5` nm | `0.4`-`0.8` nm | In-contact distance (5 Å for CA, 4 Å for closest-heavy) |
- | `simplified` (`md.compute_dssp`) | DSSP | `True` | `True`/`False` | `False` → 8-state codes |
- | `freq` (`md.baker_hubbard`) | H-bonds | `0.1` | `0.0`-`1.0` | Min frame fraction to report |
+ | `cutoff` (contact threshold, user-defined) | Contacts | `0.5` nm | `0.4`-`0.8` nm | Distance below which pair counts as in contact (5 Å is standard for CA, 4 Å for closest-heavy) |
+ | `simplified` (`md.compute_dssp`) | Secondary structure | `True` | `True`/`False` | `False` returns 8-state codes; `True` returns 3-state H/E/C |
+ | `freq` (`md.baker_hubbard`) | H-bonds | `0.1` | `0.0`-`1.0` | Min fraction of frames an H-bond must appear to be reported |
| `mode` (`md.shrake_rupley`) | SASA | `"atom"` | `"atom"`, `"residue"` | Per-atom or summed per residue |
- | `probe_radius` (`md.shrake_rupley`) | SASA | `0.14` nm | `0.12`-`0.20` nm | Solvent probe (1.4 Å standard) |
- | `periodic` (geometry funcs) | Geometry | `True` | bool | Minimum-image convention if box info present |
+ | `probe_radius` (`md.shrake_rupley`) | SASA | `0.14` nm | `0.12`-`0.20` nm | Solvent probe radius (1.4 Å is the convention) |
+ | `periodic` (geometry funcs) | Geometry | `True` | bool | Apply minimum-image convention if box info is present |
## Best Practices
- 1. **Convert units before plotting** — nm → Å (×10), ps → ns (÷1000), rad → deg (`np.degrees`). #1 cause of "my RMSD is 10× too small".
- 2. **Superpose before RMSF** — `md.rmsd` superposes internally; RMSF doesn't. Without explicit `traj.superpose(...)`, RMSF picks up rigid-body motion.
- 3. **`md.iterload` for big trajectories** — 1 µs at 1 ps/frame is millions of frames. Streaming avoids OOM.
- 4. **`simplified=False` for fine-grained DSSP** — 3-state collapses 3_10 / pi-helix / bend transitions.
- 5. **Pair indices, not `resSeq`** — `md.compute_contacts` returns 0-based contiguous indices. Map to `resSeq` only for labeling.
- 6. **Strip waters early** — `traj.atom_slice(traj.topology.select("protein"))` cuts memory/time 5–20× on solvated systems.
+ 1. **Always convert units before plotting**: mdtraj returns nm / ps / radians. Multiply by 10 (Å), divide by 1000 (ns), wrap in `np.degrees()` (deg). Forgetting this is the #1 cause of "my RMSD looks 10× too small".
+ ```python
+ rmsd_ang = md.rmsd(traj, traj, 0) * 10.0
+ ```
+
+ 2. **Superpose before RMSD/RMSF**: `md.rmsd` superposes internally each call. For RMSF you compute manually — call `traj.superpose(...)` first or RMSF will include rigid-body motion.
+
+ 3. **Use `md.iterload` for trajectories that don't fit in RAM**: a 1 µs all-atom trajectory at 1 ps/frame is millions of frames. Streaming avoids OOM and is often faster than `md.load` for selective metrics.
+
+ 4. **Don't trust `simplified=True` DSSP for fine-grained analysis**: collapsing 8 states to 3 hides 3_10 / pi-helix / bend transitions. Use `simplified=False` whenever the secondary structure signal matters.
+
+ 5. **Pair indices, not `resSeq`**: `md.compute_contacts` returns **0-based residue indices** (contiguous across chains), not the PDB `resSeq`. Map back via `traj.topology.residue(i).resSeq` only when labeling.
+
+ 6. **Strip waters early**: `traj = traj.atom_slice(traj.topology.select("protein"))` before downstream analysis cuts memory and speed by 5–20× on solvated systems.
+
## Common Recipes
### Recipe: Convert Trajectory Format
+ When to use: convert AMBER NetCDF to GROMACS XTC, or DCD to multi-frame PDB.
+
```python
import mdtraj as md
traj = md.load("input.nc", top="topology.parm7")
traj.save_xtc("output.xtc")
- traj.save_pdb("output.pdb")
- print(f"Converted {traj.n_frames} frames -> output.xtc")
+ traj.save_pdb("output.pdb") # whole trajectory in one multi-MODEL PDB
+ print(f"Converted {traj.n_frames} frames → output.xtc")
```
- ### Recipe: Per-Residue Mean ± SD RMSF Across Replicates
+ ### Recipe: Compute Per-Residue Average RMSF Across Replicates
- Ensemble flexibility input for flexibility-aware docking or homology refinement.
+ When to use: ensemble-averaged flexibility, used as input to flexibility-aware docking or homology refinement.
```python
- import mdtraj as md, numpy as np, pandas as pd
+ import mdtraj as md
+ import numpy as np
+ import pandas as pd
paths = {"rep1": "rep1.xtc", "rep2": "rep2.xtc", "rep3": "rep3.xtc"}
- top = "protein.pdb"
- rmsf_by_rep = {}
+ top = "protein.pdb"
+ rmsf_per_rep = {}
+
for name, p in paths.items():
t = md.load(p, top=top)
- ct = t.atom_slice(t.topology.select("name CA")); ct.superpose(ct, frame=0)
- rmsf_by_rep[name] = np.sqrt(((ct.xyz - ct.xyz.mean(axis=0)) ** 2).sum(axis=2).mean(axis=0)) * 10.0
+ ca = t.topology.select("name CA")
+ ct = t.atom_slice(ca); ct.superpose(ct, frame=0)
+ diff = ct.xyz - ct.xyz.mean(axis=0)
+ rmsf_per_rep[name] = np.sqrt((diff ** 2).sum(axis=2).mean(axis=0)) * 10.0
- df = pd.DataFrame(rmsf_by_rep)
- df["mean"], df["std"] = df.mean(axis=1), df.std(axis=1)
+ df = pd.DataFrame(rmsf_per_rep)
+ df["mean"] = df.mean(axis=1)
+ df["std"] = df.std(axis=1)
df.to_csv("rmsf_replicates.csv", index_label="residue_index")
+ print(df.head())
```
- ### Recipe: Save Frames Above RMSD Threshold
+ ### Recipe: Save Only Frames Where RMSD > Threshold
- Curate "open" conformations for clustering or downstream analysis.
+ When to use: cluster late conformational states or build a curated set of "open" conformations.
```python
import mdtraj as md
traj = md.load("rep1.xtc", top="protein.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
- open_state = traj[rmsd_ang > 4.0]
+ open_state = traj[rmsd_ang > 4.0] # frames with RMSD > 4 Å
open_state.save_dcd("open_state.dcd")
- print(f"Saved {open_state.n_frames}/{traj.n_frames} frames")
+ print(f"Saved {open_state.n_frames}/{traj.n_frames} frames as open state")
```
- ### Recipe: Cross-Replicate Helix Fraction
+ ### Recipe: Cross-Replicate DSSP Helix Fraction
+ When to use: quantitative summary of secondary structure stability across replicates.
+
```python
- import mdtraj as md, numpy as np
+ import mdtraj as md
+ import numpy as np
paths = ["rep1.xtc", "rep2.xtc", "rep3.xtc"]
- helix_fracs = [(md.compute_dssp(md.load(p, top="protein.pdb"), simplified=False) == "H").mean() for p in paths]
- print(f"α-helix per replicate: {[f'{h:.2f}' for h in helix_fracs]}")
+ top = "protein.pdb"
+ helix_fracs = []
+ for p in paths:
+ t = md.load(p, top=top)
+ dssp = md.compute_dssp(t, simplified=False)
+ helix_fracs.append((dssp == "H").mean())
+ print(f"α-helix fraction per replicate: {[f'{h:.2f}' for h in helix_fracs]}")
print(f"Mean ± std: {np.mean(helix_fracs):.2f} ± {np.std(helix_fracs):.2f}")
```
## Troubleshooting
| Problem | Cause | Solution |
|---------|-------|----------|
- | `RuntimeError: No topology` | XTC/DCD/TRR carry no topology | Pass `top="protein.pdb"` (or `.prmtop`) to `md.load` |
- | RMSD 10× too small | mdtraj returns nm | Multiply by 10.0 |
- | RMSF spikes at termini | Intrinsic terminal flexibility, not a bug | Slice termini off or report explicitly |
- | DSSP returns only `H/E/C` | `simplified=True` default | Pass `simplified=False` |
- | Empty/`' '` DSSP cells | Space = coil, not missing | `dssp[dssp == " "] = "C"` |
- | MemoryError on big trajectory | Whole trajectory in RAM | Switch to `md.iterload(..., chunk=500)` |
- | phi/psi arrays differ in residue count | First residue has no phi; last has no psi | Intersect residue index sets before pairing |
- | Ramachandran off-center | Angles still in radians | Wrap with `np.degrees`; set `range=[[-180, 180], [-180, 180]]` |
- | `md.compute_contacts` returns NaN | Some pairs lack heavy atoms for `closest-heavy` (e.g., glycine sidechain) | Use `scheme="closest"` or `"ca"`; drop NaN columns |
- | Contact map asymmetric | Only filled upper triangle | Mirror: `freq_map[j, i] = freq_map[i, j]` |
- | `traj.time` all zeros | Some writers leave time blank | Reconstruct: `traj.time = np.arange(traj.n_frames) * dt_ps` |
+ | `RuntimeError: No topology` | Trajectory format (XTC/DCD/TRR) carries no topology | Pass `top="protein.pdb"` or `top="topology.prmtop"` to `md.load` |
+ | RMSD values 10× too small | mdtraj returns nm, not Å | Multiply by `10.0` before plotting/reporting |
+ | RMSF spikes at termini | Termini have higher intrinsic flexibility; not a bug | Slice them off (`ct.atom_slice(ct.topology.select("residue 2 to -2"))`) or report explicitly |
+ | `md.compute_dssp` returns only `H/E/C` | `simplified=True` is the default | Pass `simplified=False` for 8-state output |
+ | Empty DSSP cells / `' '` codes | Spaces represent coil — they aren't missing | Replace with `'C'` before plotting: `dssp[dssp == " "] = "C"` |
+ | Memory error on large trajectory | Whole trajectory loaded into RAM | Switch to `md.iterload(path, top=top, chunk=500)` |
+ | Phi/psi arrays have different residue counts | First residue has no phi; last has no psi | Intersect the residue index sets before pairing |
+ | Ramachandran plot off-center | Forgot `np.degrees(...)`; angles are radians | Wrap with `np.degrees` and set `range=[[-180, 180], [-180, 180]]` in `histogram2d` |
+ | `md.compute_contacts` returns NaNs | Some pairs lack heavy atoms for `closest-heavy` (e.g., glycine sidechain) | Use `scheme="closest"` or `"ca"`; drop NaN columns before plotting |
+ | Contact-map heatmap is asymmetric | Only filled the upper triangle | Mirror with `freq_map[j, i] = freq_map[i, j]` after filling |
+ | Trajectory has no `time` array | Some writers store every step as `time=0` | Reconstruct: `traj.time = np.arange(traj.n_frames) * dt_ps` |
## Related Skills
- - **mdanalysis-trajectory** — richer atom-selection grammar and AnalysisBase framework
- - **autodock-vina-docking** / **smina-molecular-docking** — feed representative MD frames as docking targets
- - **chembl-database-bioactivity** — overlay ligand bioactivity on MD-derived flexibility maps
+ - **mdanalysis-trajectory** — alternative trajectory library with richer atom-selection grammar and the AnalysisBase framework
+ - **autodock-vina-docking** / **smina-molecular-docking** — feed representative frames from MD as docking targets
+ - **opentargets-database** / **chembl-database-bioactivity** — fetch ligand bioactivity to overlay on MD-derived flexibility maps
## References
- [mdtraj documentation](https://mdtraj.org/) — official docs, API reference, tutorials
- [McGibbon et al. (2015)](https://doi.org/10.1016/j.bpj.2015.08.015) — "MDTraj: A Modern Open Library for the Analysis of Molecular Dynamics Trajectories", *Biophys J*
- - [DSSP algorithm](https://doi.org/10.1002/bip.360221211) — Kabsch & Sander (1983)
- - [Baker–Hubbard H-bond criterion](https://doi.org/10.1016/S0006-3495(99)77244-7) — Baker & Hubbard (1984)
- - [Shrake–Rupley SASA](https://doi.org/10.1016/0022-2836(73)90011-9) — Shrake & Rupley (1973)
- - [mdtraj GitHub](https://github.com/mdtraj/mdtraj) — source, issues, examples
+ - [DSSP algorithm](https://doi.org/10.1002/bip.360221211) — Kabsch & Sander (1983) "Dictionary of protein secondary structure"
+ - [Baker–Hubbard H-bond criterion](https://doi.org/10.1016/S0006-3495(99)77244-7) — Baker & Hubbard (1984) "Hydrogen bonding in globular proteins"
+ - [Shrake–Rupley SASA](https://doi.org/10.1016/0022-2836(73)90011-9) — Shrake & Rupley (1973) "Environment and exposure to solvent of protein atoms"
+ - [mdtraj on GitHub](https://github.com/mdtraj/mdtraj) — source, issue tracker, examples