Molecular dynamics (MD) simulations numerically integrate Newton's equations of motion for all atoms in a system, producing a time-resolved trajectory of atomic positions. In drug discovery, MD is used to study protein conformational flexibility, validate and refine docking poses, compute binding free energies more accurately than docking scoring functions, and understand the mechanism of drug dissociation.

MD is computationally expensive — a 100 ns simulation of a protein-ligand complex in explicit water (~50,000 atoms) takes 4–24 hours on a modern GPU — but it provides information about dynamic processes that static crystal structures cannot capture. This guide uses Python with the OpenMM simulation toolkit, which runs on CUDA GPUs through a Python API.

System Preparation

Force Field Selection

Component	Recommended Force Field	Notes
Protein	AMBER ff14SB	Best validated for protein structure and dynamics
Small molecule (GAFF)	GAFF2 + AM1-BCC charges	Widely used, some known issues with rings
Small molecule (OpenFF)	OpenFF Sage 2.1	Better torsional parameters, open source
Water	TIP3P or OPC	TIP3P is fast; OPC better for large proteins
Ions	Joung-Cheatham (AMBER compatible)	Na⁺/Cl⁻ with TIP3P

from openmm.app import *
from openmm import *
from openmm.unit import *
import numpy as np
from pathlib import Path

def prepare_protein(pdb_file: str) -> Modeller:
    """
    Load and prepare a protein structure using OpenMM Modeller:
    - Add missing hydrogens at pH 7.4
    - Add standard residues
    - Solvate in a cubic TIP3P water box

    Returns OpenMM Modeller object ready for system creation.
    """
    # Load PDB
    pdb = PDBFile(pdb_file)

    # Force field: protein ff14SB + TIP3P water
    forcefield = ForceField("amber14-all.xml", "amber14/tip3p.xml")

    # Create Modeller to handle missing atoms and protonation
    modeller = Modeller(pdb.topology, pdb.positions)
    modeller.addHydrogens(forcefield, pH=7.4)

    # Solvate: 10 Å padding on each side, 0.15 M NaCl
    modeller.addSolvent(
        forcefield,
        model="tip3p",
        padding=10.0 * angstroms,
        ionicStrength=0.15 * molar,
    )
    return modeller


def create_openmm_system(modeller: Modeller) -> tuple:
    """
    Create an OpenMM System for GPU simulation.

    Returns (system, simulation) tuple ready for energy minimisation.
    """
    forcefield = ForceField("amber14-all.xml", "amber14/tip3p.xml")

    # Create system with PME electrostatics and SHAKE constraints on H-bonds
    system = forcefield.createSystem(
        modeller.topology,
        nonbondedMethod=PME,
        nonbondedCutoff=10.0 * angstroms,
        constraints=HBonds,
        hydrogenMass=1.5 * amu,  # Allows 4 fs timestep
    )

    # NPT integrator: Langevin thermostat (300K) + Monte Carlo barostat (1 atm)
    integrator = LangevinMiddleIntegrator(300 * kelvin, 1.0 / picosecond, 4.0 * femtoseconds)
    system.addForce(MonteCarloBarostat(1 * bar, 300 * kelvin))

    # Use GPU via CUDA platform
    platform = Platform.getPlatformByName("CUDA")
    simulation = Simulation(modeller.topology, system, integrator, platform)
    simulation.context.setPositions(modeller.positions)

    return system, simulation

Ligand Parametrization with OpenFF

from openff.toolkit.topology import Molecule, Topology
from openff.toolkit.typing.engines.smirnoff import ForceField as OFFForceField
from openmmforcefields.generators import SMIRNOFFTemplateGenerator
from openmm.app import ForceField

def parametrize_ligand_openff(ligand_smiles: str) -> SMIRNOFFTemplateGenerator:
    """
    Generate OpenFF Sage 2.1.0 force field parameters for a small molecule.
    These parameters can be injected into an OpenMM ForceField for protein-ligand simulations.

    Args:
        ligand_smiles: SMILES string of the ligand (with 3D conformer from docking)

    Returns:
        SMIRNOFFTemplateGenerator that can be registered with OpenMM ForceField
    """
    # Load molecule and assign partial charges (AM1-BCC via ELF10)
    molecule = Molecule.from_smiles(ligand_smiles, allow_undefined_stereo=True)
    molecule.generate_conformers(n_conformers=10)
    molecule.assign_partial_charges(partial_charge_method="am1bccelf10")

    # Create SMIRNOFF template generator from OpenFF Sage
    generator = SMIRNOFFTemplateGenerator(
        molecules=[molecule],
        forcefield="openff_unconstrained-2.1.0.offxml",
    )
    return generator


def build_complex_forcefield(ligand_smiles: str) -> ForceField:
    """
    Build a combined protein (AMBER ff14SB) + ligand (OpenFF Sage) ForceField.
    """
    generator = parametrize_ligand_openff(ligand_smiles)
    forcefield = ForceField("amber14-all.xml", "amber14/tip3p.xml")
    forcefield.registerTemplateGenerator(generator.generator)
    return forcefield

Production MD Protocol

from openmm.app import *
from openmm import *
from openmm.unit import *

def run_production_md(
    simulation: Simulation,
    n_steps: int = 25_000_000,   # 25M steps × 4 fs = 100 ns
    report_interval: int = 10_000,  # Save every 40 ps
    trajectory_file: str = "trajectory.dcd",
    log_file: str = "md.log",
) -> None:
    """
    Run production NPT MD simulation.

    Energy minimisation and NVT/NPT equilibration are assumed done prior.

    Args:
        simulation: OpenMM Simulation object (already equilibrated)
        n_steps: Number of MD steps
        report_interval: Steps between trajectory/log writes
        trajectory_file: DCD output file path
        log_file: State data log path
    """
    # Add reporters
    simulation.reporters.append(
        DCDReporter(trajectory_file, report_interval)
    )
    simulation.reporters.append(
        StateDataReporter(
            log_file,
            report_interval,
            step=True, time=True, potentialEnergy=True,
            kineticEnergy=True, totalEnergy=True,
            temperature=True, volume=True, density=True,
            speed=True, progress=True, remainingTime=True,
            totalSteps=n_steps, separator="\t",
        )
    )

    print(f"Running {n_steps} steps ({n_steps * 4e-6:.1f} ns)...")
    simulation.step(n_steps)
    print("Production MD complete.")

Trajectory Analysis: RMSD and RMSF

import numpy as np
import mdtraj as md  # pip install mdtraj

def analyse_trajectory(
    topology_file: str,
    trajectory_file: str,
    ligand_resname: str = "LIG",
    protein_backbone: bool = True,
) -> dict:
    """
    Analyse an MD trajectory for protein stability (RMSD) and flexibility (RMSF).

    Args:
        topology_file: PDB file with topology (first frame reference)
        trajectory_file: DCD trajectory file
        ligand_resname: Residue name of the ligand for pose stability RMSD
        protein_backbone: If True, compute backbone RMSD; otherwise, all heavy atoms

    Returns:
        Dict with numpy arrays: time_ns, protein_rmsd_A, per_residue_rmsf_A, ligand_rmsd_A
    """
    # Load trajectory with mdtraj
    traj = md.load(trajectory_file, top=topology_file)
    reference = traj[0]

    # Protein backbone atom selection
    if protein_backbone:
        protein_atoms = traj.topology.select("backbone")
    else:
        protein_atoms = traj.topology.select("protein and not element H")

    # Align trajectory to reference on protein backbone
    traj.superpose(reference, atom_indices=protein_atoms)

    # Protein RMSD (backbone, vs reference)
    protein_rmsd = md.rmsd(traj, reference, atom_indices=protein_atoms) * 10  # nm → Å

    # Per-residue RMSF
    rmsf_nm = md.rmsf(traj, reference, atom_indices=protein_atoms)
    # Map back to residue index (take mean RMSF over backbone atoms per residue)
    rmsf_per_residue_A = []
    ca_atoms = traj.topology.select("name CA")
    rmsf_ca = md.rmsf(traj, reference, atom_indices=ca_atoms) * 10  # nm → Å

    # Ligand RMSD (pose stability)
    lig_atoms = traj.topology.select(f"resname {ligand_resname} and not element H")
    if len(lig_atoms) > 0:
        ligand_rmsd = md.rmsd(traj, reference, atom_indices=lig_atoms) * 10  # nm → Å
    else:
        ligand_rmsd = np.array([])

    # Time axis: assume 40 ps between frames (10,000 steps × 4 fs)
    time_ns = np.arange(len(traj)) * 0.04  # ns

    return {
        "time_ns": time_ns,
        "protein_rmsd_A": protein_rmsd,
        "rmsf_per_ca_A": rmsf_ca,
        "ligand_rmsd_A": ligand_rmsd,
        "mean_protein_rmsd": float(protein_rmsd.mean()),
        "max_ligand_rmsd": float(ligand_rmsd.max()) if len(ligand_rmsd) > 0 else None,
    }

Binding Free Energy with MBAR

Multistate Bennett Acceptance Ratio (MBAR) can estimate ΔΔG between two ligands or between different thermodynamic states by combining data from multiple MD simulations at different lambda windows. The Python library pymbar implements this efficiently.

import numpy as np
from pymbar import MBAR  # pip install pymbar

def estimate_ddg_mbar(
    u_kln: np.ndarray,
    n_samples_k: list[int],
    temperature: float = 300.0,
) -> dict:
    """
    Estimate free energy differences using MBAR.

    Args:
        u_kln: Reduced potential energy matrix, shape (K, K, N_max)
                u_kln[k, l, n] = reduced potential energy of sample n
                from state k evaluated in state l
        n_samples_k: List of number of samples from each state k
        temperature: Simulation temperature in Kelvin

    Returns:
        Dict with free energy differences (kcal/mol) and uncertainty
    """
    mbar = MBAR(u_kln, n_samples_k, verbose=False)
    results = mbar.getFreeEnergyDifferences(return_dict=True)

    # Convert from kT to kcal/mol
    kT = 0.001987 * temperature  # kcal/mol at T K
    delta_f = results["Delta_f"] * kT
    d_delta_f = results["dDelta_f"] * kT

    # ΔF from state 0 to state K-1
    n_states = len(n_samples_k)
    ddg = float(delta_f[0, n_states - 1])
    ddg_err = float(d_delta_f[0, n_states - 1])

    return {
        "ddg_kcal_mol": round(ddg, 3),
        "ddg_uncertainty_kcal_mol": round(ddg_err, 3),
        "delta_f_matrix_kcal_mol": delta_f.tolist(),
    }

Enhanced Sampling: Metadynamics Overview

Standard MD can get trapped in local energy minima and fail to explore all relevant conformational states within accessible simulation times. Enhanced sampling methods accelerate transitions across energy barriers:

Metadynamics: adds history-dependent bias along user-defined collective variables (CVs) to push the system out of visited states. Implemented in PLUMED (plumed.org), which interfaces with OpenMM via a PLUMED reporter.
Replica Exchange MD (REMD): runs multiple replicas at different temperatures and periodically swaps coordinates between adjacent replicas, using high-temperature replicas to cross barriers.
Steered MD / SMD: applies a harmonic restraint moving over time to force the ligand out of the binding site, measuring the work required to estimate binding free energy via Jarzynski's equality.
Adaptive Sampling (OpenPathSampling, DeepDriveMD): uses ML models to predict which conformations to seed next, maximising conformational coverage with minimal simulation time.

Key Takeaways

Use AMBER ff14SB + OpenFF Sage 2.1 for protein-ligand simulations; OpenFF has better torsional parameters for small molecules than GAFF2, especially for heterocycles.
For a 100 ns simulation: energy minimise → 1 ns NVT → 10 ns NPT equilibration → production run. Always verify RMSD convergence (<2 Å) before trusting analysis.
Ligand RMSD > 2 Å during MD indicates the docked pose is unstable and the binding mode should be revisited before any free energy calculations.
MBAR (via pymbar) is the gold-standard method for combining alchemical windows in free energy calculations — it uses all data from all states optimally, unlike TI or BAR.
For rough binding affinity ranking of 5–50 compounds, metadynamics-based ΔG estimation (with PLUMED) balances accuracy and compute cost better than full FEP.

Try it live

Explore a related model hands-on: open Molecular Dynamics Simulation and change its parameters while it runs. Everything stays in your browser.

▶ Open Molecular Dynamics Simulation simulation

Molecular Dynamics for Drug Discovery