Molecular Dynamics Guide

Force Field Selection

Common Force Fields

OpenMM System Setup

from openmm.app import PDBFile, ForceField, Modeller, Simulation
from openmm.app import PME, HBonds, NoCutoff
from openmm import LangevinMiddleIntegrator, MonteCarloBarostat
from openmm.unit import kelvin, atmospheres, nanometers, picoseconds

def setup_openmm_simulation(pdb_path: str,
                              temperature: float = 300,
                              pressure: float = 1.0,
                              timestep: float = 0.002) -> Simulation:
    """
    Set up an OpenMM molecular dynamics simulation.
    pdb_path: path to prepared PDB file
    temperature: simulation temperature in Kelvin
    pressure: pressure in atmospheres
    timestep: integration timestep in picoseconds
    """
    pdb = PDBFile(pdb_path)
    forcefield = ForceField("amber14-all.xml", "amber14/tip3pfb.xml")

    modeller = Modeller(pdb.topology, pdb.positions)
    modeller.addSolvent(forcefield, padding=1.0 * nanometers,
                        ionicStrength=0.15)

    system = forcefield.createSystem(
        modeller.topology,
        nonbondedMethod=PME,
        nonbondedCutoff=1.0 * nanometers,
        constraints=HBonds,
    )

    # Barostat for NPT ensemble
    system.addForce(
        MonteCarloBarostat(pressure * atmospheres, temperature * kelvin)
    )

    integrator = LangevinMiddleIntegrator(
        temperature * kelvin,
        1.0 / picoseconds,
        timestep * picoseconds,
    )

    simulation = Simulation(modeller.topology, system, integrator)
    simulation.context.setPositions(modeller.positions)

    # Energy minimization
    simulation.minimizeEnergy()

    return simulation

Trajectory Analysis

Structural Analysis with MDAnalysis

import MDAnalysis as mda
from MDAnalysis.analysis import rms, align, diffusionmap
import numpy as np

def analyze_trajectory(topology: str, trajectory: str) -> dict:
    """
    Comprehensive trajectory analysis: RMSD, RMSF, radius of gyration.
    topology: topology file (GRO, PDB, PSF)
    trajectory: trajectory file (XTC, TRR, DCD)
    """
    u = mda.Universe(topology, trajectory)
    protein = u.select_atoms("protein and name CA")

    # RMSD over time (C-alpha atoms)
    ref = mda.Universe(topology)
    rmsd_analysis = rms.RMSD(u, ref, select="protein and name CA")
    rmsd_analysis.run()
    rmsd_data = rmsd_analysis.results.rmsd  # shape: (n_frames, 3)

    # RMSF per residue
    align.AlignTraj(u, ref, select="protein and name CA", in_memory=True).run()
    rmsf = rms.RMSF(protein).run()

    # Radius of gyration
    rg_values = []
    for ts in u.trajectory:
        rg_values.append(protein.radius_of_gyration())

    return {
        "n_frames": len(u.trajectory),
        "rmsd_mean_nm": np.mean(rmsd_data[:, 2]) / 10,  # A to nm
        "rmsd_final_nm": rmsd_data[-1, 2] / 10,
        "rmsf_mean_nm": np.mean(rmsf.results.rmsf) / 10,
        "rg_mean_nm": np.mean(rg_values) / 10,
        "rg_std_nm": np.std(rg_values) / 10,
        "simulation_time_ns": u.trajectory[-1].time / 1000,
    }

Hydrogen Bond Analysis

from MDAnalysis.analysis.hydrogenbonds import HydrogenBondAnalysis

def analyze_hbonds(universe: mda.Universe,
                    donor_sel: str = "protein",
                    acceptor_sel: str = "protein") -> dict:
    """Analyze hydrogen bonds over the trajectory."""
    hbonds = HydrogenBondAnalysis(
        universe,
        donors_sel=f"({donor_sel}) and (name N* or name O*)",
        acceptors_sel=f"({acceptor_sel}) and (name O* or name N*)",
        d_a_cutoff=3.5,
        d_h_a_angle_cutoff=150,
    )
    hbonds.run()

    return {
        "total_hbonds_detected": len(hbonds.results.hbonds),
        "mean_per_frame": len(hbonds.results.hbonds) / hbonds.n_frames,
        "unique_pairs": len(set(
            (int(r[1]), int(r[3])) for r in hbonds.results.hbonds
        )),
    }

Free Energy Methods

Umbrella Sampling

Umbrella sampling computes the potential of mean force (PMF) along a reaction coordinate:

1. Generate windows along the reaction coordinate (e.g., distance between two groups)
2. Run restrained simulations at each window with a harmonic bias
3. Combine windows using WHAM (Weighted Histogram Analysis Method)
4. Report the free energy profile (PMF)

Alchemical Free Energy Perturbation

Used for computing binding free energies and solvation free energies:

Lambda schedule: 0.0, 0.1, 0.2, ..., 0.9, 1.0
At lambda=0: full interaction (bound state)
At lambda=1: no interaction (unbound state)

Each lambda window: independent MD simulation
Analysis: MBAR or TI to combine lambda windows

Tools and Software

• GROMACS: High-performance MD engine (free, GPU-accelerated)
• OpenMM: Python-native MD with GPU support
• AMBER: Comprehensive MD package (academic license)
• NAMD: Scalable MD for large biomolecular systems
• MDAnalysis: Python trajectory analysis library
• MDTraj: Lightweight trajectory analysis
• PyMOL / VMD: Molecular visualization and movie generation
• PLUMED: Free energy and enhanced sampling methods plugin