Run and analyze molecular dynamics simulations with OpenMM and MDAnalysis. Set up protein and protein-ligand systems with PDBFixer, choose force fields and water models (AMBER14, CHARMM36m, ff19SB, GAFF2, TIP3P), solvate and add ions, run energy minimization, NVT/NPT equilibration and production MD on GPU, then analyze trajectories for RMSD, RMSF, radius of gyration, hydrogen bonds, native contacts, PCA and free energy surfaces. Use this skill for protein stability under mutation, ligand binding-mode and residence-time questions, conformational sampling, membrane proteins, and disordered ensembles. Also trigger on OpenMM, MDAnalysis, mdtraj, Simulation.step, LangevinMiddleIntegrator, PDBFixer, DCD or XTC trajectory, RMSD analysis, or production MD.
Run and analyze molecular dynamics simulations with OpenMM and MDAnalysis. Set up protein and protein-ligand systems with PDBFixer, choose force fields and water models (AMBER14, CHARMM36m, ff19SB, GAFF2, TIP3P), solvate and add ions, run energy minimization, NVT/NPT equilibration and production MD on GPU, then analyze trajectories for RMSD, RMSF, radius of gyration, hydrogen bonds, native contacts, PCA and free energy surfaces. Use this skill for protein stability under mutation, ligand binding-mode and residence-time questions, conformational sampling, membrane proteins, and disordered ensembles. Also trigger on OpenMM, MDAnalysis, mdtraj, Simulation.step, LangevinMiddleIntegrator, PDBFixer, DCD or XTC trajectory, RMSD analysis, or production MD.
license
MIT
compatibility
Requires Python 3.11+ with openmm and mdanalysis, best installed from conda-forge; PDBFixer and nglview are optional extras. A CUDA or OpenCL GPU is effectively required — production MD on CPU is 10-100x slower, so nanoseconds become days. Trajectory analysis alone runs fine on CPU.
allowed-tools
Read Write Edit Bash
metadata
{"version":"1.2","skill-author":"Kuan-lin Huang"}
Molecular Dynamics
Overview
Molecular dynamics (MD) simulation computationally models the time evolution of molecular systems by integrating Newton's equations of motion. This skill covers two complementary tools:
OpenMM (https://openmm.org/): High-performance MD simulation engine with GPU support, Python API, and flexible force field support
MDAnalysis (https://mdanalysis.org/): Python library for reading, writing, and analyzing MD trajectories from all major simulation packages
Checked against: OpenMM 8.5.2 and MDAnalysis 2.10.0, August 2026.
Read references/mdanalysis_analysis.md for the trajectory
analysis catalogue — selection language, RMSD/RMSF, contacts, hydrogen bonds, PCA, clustering and
free energy surfaces — and load it when the question is about analysing a finished run rather than
producing one.
"""
Prepare an OpenMM system from a PDB file.
Args:
pdb_file: Path to cleaned PDB file (use PDBFixer for raw PDB files)
forcefield_name: Force field XML file
water_model: Water model XML file
Returns:
pdb, forcefield, system, topology
"""
# Particle Mesh Ewald for long-range electrostatics
1.0
# Constrain hydrogen bonds (allows 2 fs timestep)
True
0.0005
return
2. Energy Minimization
from openmm.app import *
from openmm import *
from openmm.unit import *
defminimize_energy(modeller, system, output_pdb="minimized.pdb",
max_iterations=1000, tolerance=10.0):
"""
Energy minimize the system to remove steric clashes.
Args:
modeller: Modeller object with topology and positions
system: OpenMM System
output_pdb: Path to save minimized structure
max_iterations: Maximum minimization steps
tolerance: Convergence criterion in kJ/mol/nm
Returns:
simulation object with minimized positions
"""# Set up integrator (doesn't matter for minimization)
integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
# Create simulation# Use GPU if available (CUDA or OpenCL), fall back to CPUtry:
platform = Platform.getPlatformByName('CUDA')
properties = {'DeviceIndex': '0', 'Precision': 'mixed'}
except Exception:
try:
platform = Platform.getPlatformByName('OpenCL')
properties = {}
except Exception:
platform = Platform.getPlatformByName('CPU')
properties = {}
simulation = Simulation(
modeller.topology, system, integrator,
platform, properties
)
simulation.context.setPositions(modeller.positions)
# Check initial energy
state = simulation.context.getState(getEnergy=True)
print(f"Initial energy: {state.getPotentialEnergy()}")
# Minimize
simulation.minimizeEnergy(
tolerance=tolerance*kilojoules_per_mole/nanometer,
maxIterations=max_iterations
)
state = simulation.context.getState(getEnergy=True, getPositions=True)
print(f"Minimized energy: {state.getPotentialEnergy()}")
# Save minimized structurewithopen(output_pdb, 'w') as f:
PDBFile.writeFile(simulation.topology, state.getPositions(), f)
return simulation
3. NVT Equilibration
from openmm.app import *
from openmm import *
from openmm.unit import *
defrun_nvt_equilibration(simulation, n_steps=50000, temperature=300,
report_interval=1000, output_prefix="nvt"):
"""
NVT equilibration: constant N, V, T.
Equilibrate velocities to target temperature.
Args:
simulation: OpenMM Simulation (after minimization)
n_steps: Number of MD steps (50000 × 2fs = 100 ps)
temperature: Temperature in Kelvin
report_interval: Steps between data reports
output_prefix: File prefix for trajectory and log
"""# Add position restraints for backbone during NVT# (Optional: restraint heavy atoms)# Set temperature
simulation.context.setVelocitiesToTemperature(temperature*kelvin)
# Add reporters
simulation.reporters = []
# Log file
simulation.reporters.append(
StateDataReporter(
f"{output_prefix}_log.txt",
report_interval,
step=True,
potentialEnergy=True,
kineticEnergy=True,
temperature=True,
volume=True,
speed=True
)
)
# DCD trajectory (compact binary format)
simulation.reporters.append(
DCDReporter(f"{output_prefix}_traj.dcd", report_interval)
)
print(f"Running NVT equilibration: {n_steps} steps ({n_steps*2/1000:.1f} ps)")
simulation.step(n_steps)
print("NVT equilibration complete")
return simulation
4. NPT Equilibration and Production
defrun_npt_production(simulation, n_steps=500000, temperature=300, pressure=1.0,
report_interval=5000, output_prefix="npt"):
"""
NPT production run: constant N, P, T.
Args:
n_steps: Production steps (500000 × 2fs = 1 ns)
temperature: Temperature in Kelvin
pressure: Pressure in bar
report_interval: Steps between reports
"""# Add Monte Carlo barostat for pressure control
system = simulation.context.getSystem()
system.addForce(MonteCarloBarostat(pressure*bar, temperature*kelvin, 25))
simulation.context.reinitialize(preserveState=True)
# Update reporters
simulation.reporters = []
simulation.reporters.append(
StateDataReporter(
f"{output_prefix}_log.txt",
report_interval,
step=True,
potentialEnergy=True,
temperature=True,
density=True,
speed=True
)
)
simulation.reporters.append(
DCDReporter(f"{output_prefix}_traj.dcd", report_interval)
)
# Save checkpoints
simulation.reporters.append(
CheckpointReporter(f"{output_prefix}_checkpoint.chk", 50000)
)
print(f"Running NPT production: {n_steps} steps ({n_steps*2/1000000:.2f} ns)")
simulation.step(n_steps)
print("Production MD complete")
return simulation
Trajectory Analysis with MDAnalysis
1. Load Trajectory
import MDAnalysis as mda
from MDAnalysis.analysis import rms, align, contacts
import numpy as np
import matplotlib.pyplot as plt
defload_trajectory(topology_file, trajectory_file):
"""
Load an MD trajectory with MDAnalysis.
Args:
topology_file: PDB, PSF, or other topology file
trajectory_file: DCD, XTC, TRR, or other trajectory
"""
u = mda.Universe(topology_file, trajectory_file)
print(f"Universe: {u.atoms.n_atoms} atoms, {u.trajectory.n_frames} frames")
print(f"Time range: 0 to {u.trajectory.totaltime:.0f} ps")
return u
2. RMSD Analysis
defcompute_rmsd(u, selection="backbone", reference_frame=0):
"""
Compute RMSD of selected atoms relative to reference frame.
Args:
u: MDAnalysis Universe
selection: Atom selection string (MDAnalysis syntax)
reference_frame: Frame index for reference structure
Returns:
numpy array of (time, rmsd) values
"""# Align trajectory to minimize RMSD
aligner = align.AlignTraj(u, u, select=selection, in_memory=True)
aligner.run()
# Compute RMSD
R = rms.RMSD(u, select=selection, ref_frame=reference_frame)
R.run()
rmsd_data = R.results.rmsd # columns: frame, time, RMSDreturn rmsd_data
defplot_rmsd(rmsd_data, title="RMSD over time", output_file="rmsd.png"):
"""Plot RMSD over simulation time."""
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(rmsd_data[:, 1] / 1000, rmsd_data[:, 2], 'b-', linewidth=0.5)
ax.set_xlabel("Time (ns)")
ax.set_ylabel("RMSD (Å)")
ax.set_title(title)
ax.axhline(rmsd_data[:, 2].mean(), color='r', linestyle='--',
label=f'Mean: {rmsd_data[:, 2].mean():.2f} Å')
ax.legend()
plt.tight_layout()
plt.savefig(output_file, dpi=150)
return fig
3. RMSF Analysis (Per-Residue Flexibility)
defcompute_rmsf(u, selection="backbone", start_frame=0):
"""
Compute per-residue RMSF (flexibility).
Returns:
resids, rmsf_values arrays
"""# Select atoms
atoms = u.select_atoms(selection)
# Compute RMSF
R = rms.RMSF(atoms)
R.run(start=start_frame)
# Average by residue
resids = []
rmsf_per_res = []
for res in u.select_atoms(selection).residues:
res_atoms = res.atoms.intersection(atoms)
iflen(res_atoms) > 0:
resids.append(res.resid)
rmsf_per_res.append(R.results.rmsf[res_atoms.indices].mean())
return np.array(resids), np.array(rmsf_per_res)
4. Protein-Ligand Contacts
defanalyze_contacts(u, protein_sel="protein", ligand_sel="resname LIG",
radius=4.5, start_frame=0):
"""
Track protein-ligand contacts over trajectory.
Args:
radius: Contact distance cutoff in Angstroms
"""
protein = u.select_atoms(protein_sel)
ligand = u.select_atoms(ligand_sel)
contact_frames = []
for ts in u.trajectory[start_frame:]:
# Find protein atoms within radius of ligand
distances = contacts.contact_matrix(
protein.positions, ligand.positions, radius
)
contact_residues = set()
for i inrange(distances.shape[0]):
if distances[i].any():
contact_residues.add(protein.atoms[i].resid)
contact_frames.append(contact_residues)
return contact_frames
# For ligand parameterization, use OpenFF toolkit or ACPYPE# uv pip install openff-toolkitfrom openff.toolkit import Molecule, ForceField as OFFForceField
from openff.interchange import Interchange
defparameterize_ligand(smiles, ff_name="openff-2.0.0.offxml"):
"""Generate GAFF2/OpenFF parameters for a small molecule."""
mol = Molecule.from_smiles(smiles)
mol.generate_conformers(n_conformers=1)
off_ff = OFFForceField(ff_name)
interchange = off_ff.create_interchange(mol.to_topology())
return interchange
Best Practices
Always minimize before MD: Raw PDB structures have steric clashes
Equilibrate before production: NVT (50–100 ps) → NPT (100–500 ps) → Production
Use GPU: Simulations are 10–100× faster on GPU (CUDA/OpenCL)
2 fs timestep with HBonds constraints: Standard; use 4 fs with HMR (hydrogen mass repartitioning)
Analyze only equilibrated trajectory: Discard first 20–50% as equilibration
Save checkpoints: MD runs can fail; checkpoints allow restart
Periodic boundary conditions: Required for solvated systems
PME for electrostatics: More accurate than cutoff methods for charged systems
Composing with the rest of the bundle
uniprot-rcsb → before: the structure, plus the check that the residues you care about are
actually resolved. Minimising a model with a 12-residue gap through your binding loop wastes the
whole run.
binding-site-analysis → before: whether the pocket is real and whether it is cryptic. MD is the
standard way to open a cryptic site, but only if you know that is the question.
autodock-vina / diffdock → before: MD is how a docked pose is tested. A pose that leaves
the site in 10 ns was not a pose. Docking scores rank; MD tells you whether the ranking
survives contact with dynamics.
boltz → before: a predicted complex, when no experimental structure exists — but treat a
predicted holo structure as a hypothesis and check its confidence before spending GPU-days.
free-energy-perturbation → after: rigorous ΔΔG. FEP is MD with an alchemical schedule and
proper convergence checks; reach for it when you need numbers rather than a movie.
degraders → after: ternary complex stability, which is exactly a question about persistence
over time rather than a static pose.
Report the ensemble, not the frame. A single snapshot from a trajectory is a screenshot of a
distribution — quote the mean and spread over the equilibrated portion, and say how much you
discarded as equilibration.
OpenMM paper: Eastman P et al. (2017) OpenMM 7: Rapid development of high performance
algorithms for molecular dynamics. PLOS Computational Biology. PMID: 28746339
MDAnalysis paper: Michaud-Agrawal N et al. (2011) J Computational Chemistry. PMID: 21500218
