domains / science / molecular-dynamics-simulation / molecular-dynamics

molecular-dynamics

Skill card

combined

The two copies describe the same OpenMM/MDAnalysis workflow, but they disagree on environment setup, version pinning, platform selection, and minimization parameters. The K-Dense copy is more rigorous (dedicated venv, pinned versions, explicit platform argument, 2 fs integrator, neutralization flag, and a warning that trajectories alone do not establish binding free energy), while the FreedomIntelligence copy adds a conda/pip install path and a CUDA/OpenCL/CPU fallback. Neither is a strict superset of the other, so crowning one would hide a real alternative.

Repositories

Other copies

Combined skill

---
name: molecular-dynamics
description: Runs and analyzes molecular dynamics simulations with OpenMM and MDAnalysis. Sets up protein/small molecule systems, defines force fields, runs energy minimization and production MD, and analyzes trajectories (RMSD, RMSF, contact maps, free energy surfaces). For structural biology, drug binding, and biophysics.
license: MIT
compatibility: Requires Python 3.11+ with OpenMM and MDAnalysis; matplotlib for plots. Optional PDBFixer and OpenFF need separate installation. Network access for installation; local simulation and analysis run offline.
metadata:
  version: "1.3"
  last-reviewed: "2026-10-01"
  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

**Reviewed versions:** OpenMM 8.6.1 and MDAnalysis 2.10.0. Reference-platform smoke checks use small synthetic systems; GPU performance and scientific convergence are not established by them. OpenMM latest API pages identify a development build, so the installed 8.6.1 API is the executable reference.

**Installation into a dedicated environment:**
```bash
uv venv --python 3.12 .venv-md
uv pip install --python .venv-md/bin/python "openmm==8.6.1" "MDAnalysis==2.10.0" matplotlib pandas
# Windows: use .venv-md/Scripts/python.exe instead
```

Alternative installation (less controlled):
```bash
conda install -c conda-forge openmm mdanalysis nglview
# or
pip install openmm mdanalysis
```

## When to Use This Skill

Use molecular dynamics when:

- **Protein stability analysis**: How does a mutation affect protein dynamics?
- **Drug binding simulations**: Characterize binding mode and residence time of a ligand
- **Conformational sampling**: Explore protein flexibility and conformational changes
- **Protein-protein interaction**: Model interface dynamics and binding energetics
- **RMSD/RMSF analysis**: Quantify structural fluctuations from a reference structure
- **Free energy estimation**: Compute binding free energy or conformational free energy
- **Membrane simulations**: Model proteins in lipid bilayers
- **Intrinsically disordered proteins**: Study IDR conformational ensembles

## Core Workflow: OpenMM Simulation

Examples require system-specific preparation and validation. Review biological assembly, alternate locations, missing loops, termini, protonation, disulfides, ligands and cofactors before parameterization. An MD trajectory alone does not establish binding free energy or residence time.

The functions below form one Python module: execute their import blocks together. Use a new output prefix for every stage to avoid overwriting earlier results.

### 1. System Preparation

```python
from openmm.app import *
from openmm import *
from openmm.unit import *

def prepare_system_from_pdb(pdb_file, forcefield_name="amber14-all.xml",
                              water_model="amber14/tip3pfb.xml", ph=7.0):
    """
    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
        ph: pH for hydrogen addition

    Returns:
        modeller, system
    """
    # Load PDB
    pdb = PDBFile(pdb_file)

    # Load force field
    forcefield = ForceField(forcefield_name, water_model)

    # Add hydrogens and solvate
    modeller = Modeller(pdb.topology, pdb.positions)
    modeller.addHydrogens(forcefield, pH=ph)

    # TIP3P geometry also serves TIP3P-FB; XML supplies the parameters.
    # 10 Å padding; added salt excludes neutralizing counterions.
    modeller.addSolvent(
        forcefield,
        model='tip3p',
        padding=10*angstroms,
        ionicStrength=0.15*molar, neutralize=True
    )

    print(f"System: {modeller.topology.getNumAtoms()} atoms, "
          f"{modeller.topology.getNumResidues()} residues")

    # Create system
    system = forcefield.createSystem(
        modeller.topology,
        nonbondedMethod=PME,         # Particle Mesh Ewald for long-range electrostatics
        nonbondedCutoff=1.0*nanometer,
        constraints=HBonds,           # Constrain bonds involving H (not intermolecular H-bonds)
        rigidWater=True,
        ewaldErrorTolerance=0.0005
    )

    return modeller, system
```

### 2. Energy Minimization

```python
from openmm.app import *
from openmm import *
from openmm.unit import *

def minimize_energy(modeller, system, output_pdb="minimized.pdb",
                     max_iterations=1000, tolerance=10.0, platform_name=None):
    """
    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
        platform_name: Optional explicit platform (e.g. 'CPU', 'CUDA', 'OpenCL')

    Returns:
        simulation object with minimized positions
    """
    # Keep the integrator at 2 fs for the subsequent NVT/NPT examples.
    integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)

    # Automatic platform selection, or an explicit tested platform (e.g. CPU).
    # A registered GPU plugin does not prove a working device or driver.
    platform = Platform.getPlatformByName(platform_name) if platform_name else None
    simulation = Simulation(modeller.topology, system, integrator, platform)
    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 structure
    with open(output_pdb, 'w') as f:
        PDBFile.writeFile(simulation.topology, state.getPositions(), f)

    return simulation
```

### 3. NVT Equilibration

```python
from openmm.app import *
from openmm import *
from openmm.unit import *

def run_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: restrain heavy atoms to reference structure)
    simulation.context.setTemperature(temperature*kelvin)
    simulation.step(n_steps)
    
    # Save trajectory
    with open(f"{output_prefix}.dcd", 'wb') as f:
        DCDFile.writeFile(simulation.topology, simulation.context.getState(getPositions=True).getPositions(), f)
    
    return simulation
```

### 4. NPT Equilibration

```python
from openmm.app import *
from openmm import *
from openmm.unit import *

def run_npt_equilibration(simulation, n_steps=50000, temperature=300,
                            pressure=1.0, report_interval=1000, output_prefix="npt"):
    """
    NPT equilibration: constant N, P, T.
    Equilibrate pressure and box dimensions.

    Args:
        simulation: OpenMM Simulation (after NVT)
        n_steps: Number of MD steps
        temperature: Temperature in Kelvin
        pressure: Pressure in bar
        report_interval: Steps between data reports
        output_prefix: File prefix for trajectory and log
    """
    # Add barostat
    simulation.context.setPressure(pressure*bar)
    simulation.step(n_steps)
    
    # Save trajectory
    with open(f"{output_prefix}.dcd", 'wb') as f:
        DCDFile.writeFile(simulation.topology, simulation.context.getState(getPositions=True).getPositions(), f)
    
    return simulation
```

### 5. Production MD

```python
from openmm.app import *
from openmm import *
from openmm.unit import *

def run_production_md(simulation, n_steps=500000, temperature=300,
                       pressure=1.0, report_interval=10000, output_prefix="production"):
    """
    Production MD: long simulation for sampling.

    Args:
        simulation: OpenMM Simulation (after NPT)
        n_steps: Number of MD steps
        temperature: Temperature in Kelvin
        pressure: Pressure in bar
        report_interval: Steps between data reports
        output_prefix: File prefix for trajectory and log
    """
    simulation.context.setTemperature(temperature*kelvin)
    simulation.context.setPressure(pressure*bar)
    simulation.step(n_steps)
    
    # Save trajectory
    with open(f"{output_prefix}.dcd", 'wb') as f:
        DCDFile.writeFile(simulation.topology, simulation.context.getState(getPositions=True).getPositions(), f)
    
    return simulation
```

### 6. Trajectory Analysis with MDAnalysis

```python
import MDAnalysis as mda
import numpy as np
import matplotlib.pyplot as plt

def analyze_trajectory(reference_pdb, trajectory_dcd, output_prefix="analysis"):
    """
    Analyze MD trajectory: RMSD, RMSF, contact maps, free energy surfaces.

    Args:
        reference_pdb: Path to reference PDB structure
        trajectory_dcd: Path to trajectory DCD file
        output_prefix: File prefix for output plots
    """
    u = mda.Universe(reference_pdb, trajectory_dcd)
    
    # RMSD
    rmsd = mda.analysis.rms.RMSD(u, reference=u.atoms, select='protein')
    rmsd.run()
    
    # RMSF
    rmsf = mda.analysis.rms.RMSF(u, reference=u.atoms, select='protein')
    rmsf.run()
    
    # Plot results
    plt.figure(figsize=(10, 6))
    plt.plot(rmsd.rmsd, label='RMSD')
    plt.xlabel('Time (ps)')
    plt.ylabel('RMSD (Å)')
    plt.title('RMSD Analysis')
    plt.legend()
    plt.savefig(f"{output_prefix}_rmsd.png", dpi=300)
    plt.close()
    
    plt.figure(figsize=(10, 6))
    plt.bar(range(len(rmsf.atoms)), rmsf.rmsf, width=1.0, align='center')
    plt.xlabel('Residue Index')
    plt.ylabel('RMSF (Å)')
    plt.title('RMSF Analysis')
    plt.savefig(f"{output_prefix}_rmsf.png", dpi=300)
    plt.close()
    
    return u
```

## Conflicts

- **K-Dense-AI/scientific-agent-skills** pins OpenMM 8.6.1 and MDAnalysis 2.10.0 in a dedicated uv venv; **FreedomIntelligence/OpenClaw-Medical-Skills** uses conda-forge or pip without version pins.
- **K-Dense-AI/scientific-agent-skills** uses a 2 fs LangevinMiddleIntegrator for minimization and keeps it for NVT/NPT; **FreedomIntelligence/OpenClaw-Medical-Skills** uses a 4 fs integrator for minimization.
- **K-Dense-AI/scientific-agent-skills** passes an explicit platform_name argument and notes that a registered GPU plugin does not prove a working device; **FreedomIntelligence/OpenClaw-Medical-Skills** hard-codes a CUDA/OpenCL/CPU fallback chain.
- **K-Dense-AI/scientific-agent-skills** adds neutralize=True to addSolvent; **FreedomIntelligence/OpenClaw-Medical-Skills** omits it.
- **K-Dense-AI/scientific-agent-skills** includes a pH parameter in addHydrogens; **FreedomIntelligence/OpenClaw-Medical-Skills** does not.