Installer avec Codex ou Claude Copiez ce prompt, collez-le dans Codex, Claude ou un autre assistant, puis laissez-le vérifier la page du skill et l'installer pour vous.
Une commande directe contourne le prompt de vérification. Examinez la source avant de l'exécuter.
Parse and analyze molecular dynamics (MD) trajectories from GROMACS, AMBER, LAMMPS,
and NAMD simulations using MDAnalysis. This skill covers RMSD/RMSF computation,
hydrogen bond tracking, contact map generation, and protein-ligand distance analysis
with publication-ready matplotlib visualizations.
When to Use This Skill
You have an MD trajectory file (.xtc, .trr, .dcd, .nc, .lammpstrj) and
want to extract quantitative structural metrics.
You need to compute RMSD (global fold stability) or RMSF (per-residue
flexibility) over a simulation.
You want to enumerate hydrogen bonds between protein and ligand or between
secondary-structure elements.
You need a contact map or residue-residue distance matrix to identify
persistent interactions.
You are measuring protein-ligand binding pocket distances to assess whether a
ligand stays bound throughout the simulation.
You want to export per-frame structural data to a pandas DataFrame for
downstream statistical analysis.
Background & Key Concepts
Universe and AtomGroup
MDAnalysis represents a simulation as a Universe object that couples a topology
file (connectivity, residue names, charges) with one or more trajectory files
(coordinate frames). Selections of atoms are returned as AtomGroup objects, which
support arithmetic, boolean masks, and iteration.
import MDAnalysis as mda
u = mda.Universe("topology.tpr", "trajectory.xtc")
protein = u.select_atoms("protein")
Trajectory Iteration
Iterating over u.trajectory moves the Universe to successive frames. All
AtomGroup.positions arrays are updated automatically — no manual frame loading is
required.
for ts in u.trajectory:
# ts.time is in ps; ts.frame is the 0-based frame index
coords = protein.positions # (N, 3) float32 array, angstroms
RMSD vs RMSF
Metric
Domain
Measures
RMSD
per-frame
global deviation from a reference
RMSF
per-atom
time-averaged fluctuation amplitude
Both metrics require alignment (superposition) to remove overall rotation/translation
before computing distances.
Hydrogen Bond Criteria
The default HBond criterion in MDAnalysis uses donor-acceptor distance ≤ 3.5 Å and
donor-H-acceptor angle ≥ 150°. These thresholds can be tuned for non-standard force
fields.
Contact Maps
A contact is defined when Cα–Cα distance (or heavy-atom distance) drops below a
cutoff (typically 8 Å for Cα, 4.5 Å for heavy atoms). A contact map averaged over
trajectory frames reveals stable structural contacts.
Environment Setup
Installation
# Create a dedicated conda environment (recommended)
conda create -n mdanalysis-env python=3.11 -y
conda activate mdanalysis-env
# Install MDAnalysis and visualization stack
pip install "MDAnalysis>=2.6""matplotlib>=3.7""numpy>=1.24""pandas>=2.0""scipy>=1.11"# Optional: install MDAnalysisTests for sample data
pip install MDAnalysisTests
Verify Installation
import MDAnalysis as mda
print(mda.__version__) # e.g. 2.6.1import MDAnalysis.tests
from MDAnalysisTests.datafiles import PSF, DCD
u = mda.Universe(PSF, DCD)
print(f"Atoms: {u.atoms.n_atoms}, Frames: {u.trajectory.n_frames}")
Supported Formats
Format
Topology
Trajectory
GROMACS
.tpr, .gro
.xtc, .trr
AMBER
.prmtop, .parm7
.nc, .ncdf
NAMD/CHARMM
.psf
.dcd
LAMMPS
.data
.lammpstrj
PDB/mmCIF
.pdb, .cif
—
Core Workflow
Step 1 — Load Universe and Inspect Topology
import MDAnalysis as mda
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# Replace with your actual files
TOPOLOGY = "md_system.tpr"# or .prmtop / .psf / .pdb
TRAJECTORY = "md_traj.xtc"# or .dcd / .nc / .lammpstrj
u = mda.Universe(TOPOLOGY, TRAJECTORY)
print(f"Number of atoms : {u.atoms.n_atoms}")
print(f"Number of residues: {u.residues.n_residues}")
print(f"Number of segments: {u.segments.n_segments}")
print(f"Number of frames : {u.trajectory.n_frames}")
print(f"Time step (ps) : {u.trajectory.dt:.3f}")
print(f"Total time (ns) : {u.trajectory.totaltime / 1000:.2f}")
# Inspect segment / chain namesfor seg in u.segments:
print(f" Segment {seg.segid}: {seg.atoms.n_atoms} atoms")
# Select key atom groups
protein = u.select_atoms("protein")
backbone = u.select_atoms("backbone")
ligand = u.select_atoms("resname LIG") # adjust resname
water = u.select_atoms("resname WAT SOL TIP3")
print(f"\nProtein atoms : {protein.n_atoms}")
print(f"Backbone atoms: {backbone.n_atoms}")
print(f"Ligand atoms : {ligand.n_atoms}")
Step 2 — RMSD Calculation and Visualization
from MDAnalysis.analysis import rms, align
# --- Align trajectory to first frame (in-place) ---
aligner = align.AlignTraj(u, u, select="backbone", in_memory=False)
aligner.run()
# --- Compute RMSD for backbone and C-alpha ---
rmsd_analysis = rms.RMSD(
u,
select="backbone",
groupselections=["backbone", "name CA"],
ref_frame=0,
)
rmsd_analysis.run(verbose=True)
# Results array shape: (n_frames, 3+n_groups)# Columns: frame, time(ps), backbone_rmsd, [group_rmsds...]
df_rmsd = pd.DataFrame(
rmsd_analysis.results.rmsd[:, 1:], # drop frame index column
columns=["Time (ps)", "Backbone RMSD (Å)", "C-alpha RMSD (Å)"],
)
df_rmsd["Time (ns)"] = df_rmsd["Time (ps)"] / 1000.0# --- Plot ---
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(df_rmsd["Time (ns)"], df_rmsd["Backbone RMSD (Å)"], label="Backbone", lw=1.5)
ax.plot(df_rmsd["Time (ns)"], df_rmsd["C-alpha RMSD (Å)"], label="Cα", lw=1.5, ls="--")
ax.set_xlabel("Time (ns)", fontsize=12)
ax.set_ylabel("RMSD (Å)", fontsize=12)
ax.set_title("Backbone and Cα RMSD over Simulation", fontsize=13)
ax.legend()
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig("rmsd_plot.png", dpi=150)
plt.show()
# Save data
df_rmsd.to_csv("rmsd_data.csv", index=False)
print("RMSD data saved to rmsd_data.csv")
pip install "MDAnalysis>=2.6"# On some systems you also need:
pip install "MDAnalysisTests>=2.6"
MemoryError with Large Trajectories
Use in_memory=False in AlignTraj and process frames in chunks:
import MDAnalysis as mda
u = mda.Universe(TOPOLOGY, TRAJECTORY)
chunk_size = 500# frames per chunkfor start inrange(0, u.trajectory.n_frames, chunk_size):
stop = min(start + chunk_size, u.trajectory.n_frames)
for ts in u.trajectory[start:stop]:
pass# process ts here
NoDataError: Universe has no bonds
Many trajectory formats strip bond information. Provide a topology file that contains
bonds (.tpr, .prmtop, .psf) rather than a bare .pdb:
u = mda.Universe("system.tpr", "traj.xtc") # bonds present in .tpr
