molecular-dynamics

分子动力学

概述

分子动力学 (MD) 模拟通过积分牛顿运动方程，以计算方式对分子系统随时间的演化进行建模。本技能涵盖两个互补的工具：

• OpenMM (https://openmm.org/)：高性能 MD 模拟引擎，支持 GPU、Python API 以及灵活的力场支持
• MDAnalysis (https://mdanalysis.org/)：用于读取、写入和分析来自所有主流模拟软件包的 MD 轨迹的 Python 库

已审查版本： OpenMM 8.6.1 和 MDAnalysis 2.10.0。参考平台冒烟 检查使用小型合成系统；GPU 性能和科学收敛性 并未通过这些检查得到确立。OpenMM 最新 API 页面标识的是一个开发构建， 因此已安装的 8.6.1 API 才是可执行的参考。

安装到专用环境：

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

何时使用本技能

在以下情况下使用分子动力学：

• 蛋白质稳定性分析：突变如何影响蛋白质动力学？
• 药物结合模拟：表征配体的结合模式和驻留时间
• 构象采样：探索蛋白质柔性和构象变化
• 蛋白质-蛋白质相互作用：对界面动力学和结合能量学进行建模
• RMSD/RMSF 分析：量化相对于参考结构的结构涨落
• 自由能估算：计算结合自由能或构象自由能
• 膜模拟：在脂质双分子层中对蛋白质进行建模
• 内在无序蛋白：研究 IDR 构象集合

核心工作流：OpenMM 模拟

示例需要针对具体系统进行准备和验证。在参数化之前，请审查生物 组装体、可变位置、缺失环、末端、质子化、二硫键、 配体和辅因子。仅凭一条 MD 轨迹并不能 确立结合自由能或驻留时间。

以下函数构成一个 Python 模块：请一并执行它们的导入块。 为每个阶段使用新的输出前缀，以避免覆盖先前的结果。

1. 系统准备

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

    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. 能量最小化

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

    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 平衡

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
    """
    # This example is unrestrained; add validated restraints before Context creation.
    if any("Barostat" in type(f).__name__ and f.getFrequency() > 0
           for f in simulation.system.getForces()):
        raise ValueError("NVT requires all barostats disabled or absent")

    # Set both the thermostat target and initial velocities.
    simulation.integrator.setTemperature(temperature*kelvin)
    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)
    )

    duration = (n_steps * simulation.integrator.getStepSize()).value_in_unit(picoseconds)
    print(f"Running NVT equilibration: {n_steps} steps ({duration:.1f} ps)")
    simulation.step(n_steps)
    print("NVT equilibration complete")

    return simulation

4. NPT 平衡与生产运行

首先使用专用的 NPT 平衡前缀调用此函数。在单独调用生产运行之前，检查密度、 能量、结构和重复稳定性；默认时长只是一个示例，并非 平衡或收敛的判据。

def run_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
    """
    # Keep the Langevin thermostat and barostat at the same temperature.
    simulation.integrator.setTemperature(temperature*kelvin)
    system = simulation.system
    barostats = [f for f in system.getForces() if "Barostat" in type(f).__name__]
    if not system.usesPeriodicBoundaryConditions():
        raise ValueError("NPT requires a periodic system")
    if not barostats:
        system.addForce(MonteCarloBarostat(pressure*bar, temperature*kelvin, 25))
        simulation.context.reinitialize(preserveState=True)
    elif len(barostats) != 1 or not isinstance(barostats[0], MonteCarloBarostat):
        raise ValueError("This example supports one isotropic MonteCarloBarostat")
    else:
        barostats[0].setDefaultPressure(pressure*bar)
        barostats[0].setDefaultTemperature(temperature*kelvin)
        barostats[0].setFrequency(25)
    # Existing Context parameters must also be updated on repeated calls.
    simulation.context.setParameter(MonteCarloBarostat.Pressure(), pressure)
    simulation.context.setParameter(MonteCarloBarostat.Temperature(), temperature)

    # 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)
    )

    duration = (n_steps * simulation.integrator.getStepSize()).value_in_unit(nanoseconds)
    print(f"Running NPT stage: {n_steps} steps ({duration:.3f} ns)")
    simulation.step(n_steps)
    simulation.saveCheckpoint(f"{output_prefix}_checkpoint.chk")
    simulation.saveState(f"{output_prefix}_state.xml")
    print("NPT stage complete")
    return simulation

使用 MDAnalysis 进行轨迹分析

1. 加载轨迹

使用与 DCD 原子数和顺序完全一致的溶剂化拓扑（例如上面的 minimized.pdb），而不是原始的未溶剂化输入。MDAnalysis 默认将长度转换为 Å，将时间转换为 ps；OpenMM 的裸坐标使用 nm。请保留 帧的盒子矢量和实际时间戳；不要根据帧编号推断时间。

import MDAnalysis as mda
from MDAnalysis.analysis import rms, align, contacts
import numpy as np
import matplotlib.pyplot as plt

def load_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")
    first, last = u.trajectory[0].time, u.trajectory[-1].time
    print(f"Time range: {first:g} to {last:g} ps")
    u.trajectory[0]
    return u

2. RMSD 分析

def compute_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 with columns [frame index, time in ps, RMSD in Angstroms]
    """
    # RMSD performs its own fit; avoid rotating the stored trajectory here.
    # Make the selected molecule whole before analysis of periodic trajectories.
    R = rms.RMSD(u, select=selection, ref_frame=reference_frame)
    R.run()

    rmsd_data = R.results.rmsd  # columns: frame, time, RMSD
    return rmsd_data

def plot_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 分析（逐残基柔性）

在进行 RMSF 之前，先将分子补全为完整分子，并对一个单独的结构分析副本进行对齐； RMSF 不会进行对齐。不要对盒子未被旋转的已旋转轨迹运行周期性接触分析。有关对齐、PCA、DSSP、氢键和群体导出自由能，请参见分析参考。

def compute_rmsf(u, selection="protein and name CA", start_frame=0):
    """
    Compute per-residue RMSF (flexibility).

    Returns:
        residue_keys, rmsf_values; keys are (resindex, segid, resid, resname)
    """
    # Select atoms
    atoms = u.select_atoms(selection)

    if not len(atoms):
        raise ValueError("RMSF selection is empty")
    # A multi-atom selection returns mean atomic RMSF, not COM 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)
        if len(res_atoms) > 0:
            resids.append((res.ix, res.segid, res.resid, res.resname))
            # RMSF is indexed within the selected AtomGroup, not the Universe.
            rmsf_per_res.append(R.results.rmsf[atoms.resindices == res.ix].mean())

    return resids, np.array(rmsf_per_res)

4. 蛋白质-配体接触

def analyze_contacts(u, protein_sel="protein", ligand_sel="resname LIG",
                      radius=4.5, start_frame=0, periodic=True):
    """
    Track protein-ligand contacts over trajectory.

    Args:
        radius: Atom-pair cutoff in Angstroms; any pair defines a residue contact.
        periodic: Require box data and use minimum-image distances.
    Returns sets of (resindex, segid, resid, resname), one per analyzed frame.
    """
    from MDAnalysis.lib.distances import distance_array

    protein = u.select_atoms(protein_sel)
    ligand = u.select_atoms(ligand_sel)

    if not len(protein) or not len(ligand):
        raise ValueError("Protein and ligand selections must both be nonempty")
    if not np.isfinite(radius) or radius <= 0:
        raise ValueError("radius must be finite and positive")
    contact_frames = []
    for ts in u.trajectory[start_frame:]:
        if periodic and (ts.dimensions is None or not np.all(np.isfinite(ts.dimensions))
                         or np.any(ts.dimensions[:3] <= 0)):
            raise ValueError("Periodic contacts require valid frame box dimensions")
        # This includes hydrogen atoms unless the caller selects heavy atoms.
        distances = contacts.contact_matrix(
            distance_array(protein.positions, ligand.positions, box=ts.dimensions if periodic else None),
            radius,
        )
        contact_residues = set()
        for i in range(distances.shape[0]):
            if distances[i].any():
                res = protein[i].residue
                contact_residues.add((res.ix, res.segid, res.resid, res.resname))
        contact_frames.append(contact_residues)

    return contact_frames

力场与准备

为科学体系选择经过验证的蛋白质/配体/水/离子组合。 工作示例保留 AMBER14/TIP3P-FB；这并不是普适推荐。 当前 OpenMM 还内置了 amber19-all.xml（ff19SB、DNA OL21、RNA OL3、lipid21）， 以及 charmm36_2024.xml 及其自带的 charmm36_2024/water.xml。通用 TIP3P 与 CHARMM 修饰的 TIP3P 不可互换。IDP 系综对蛋白质-水平衡尤为 敏感。四点位水需要额外的粒子以及与之匹配的 Modeller 几何结构，这与上述三点位示例不同。

AMBER 蛋白质 XML 无法为任意配体参数化。请使用经过审查的 GAFF2 或 OpenFF 路线，并明确立体化学、质子化、键级、原子 映射和电荷方法。仅含配体的 Interchange 并不是蛋白质-配体体系。 有关保守的 PDBFixer 和 OpenFF 示例、安装要求及当前发行版详情，请参见准备参考。

验证与溯源

• 从最小化开始；能量有限并不能证明结构有效。
• 使用分阶段的 NVT → NPT 平衡 → 生产。此示例没有约束。
• 对含 H 键施加约束并采用 2 fs 步长是一个起始选择。更大的 步长/HMR 需要针对该体系进行积分器、稳定性和可观测量验证。
• 使用逐帧盒子计算周期性距离。使用可信的键拓扑使分子对结构 可观测量保持完整；解缠无法推断缺失的键。
• 评估预热、自相关、有效样本、副本和不确定性。 短暂稳定的温度/密度轨迹不能确立构象采样。
• 保存拓扑/原子顺序、System/Integrator XML、力场文件/版本、 种子、平台/精度、参数和分析选择。检查点与 平台/硬件/版本相关；XML 状态更易移植，但不能 保留所有内部随机数生成器状态。测试实际设置的重新启动。

其他资源

• OpenMM 文档：https://openmm.org/documentation.html
• MDAnalysis 用户指南：https://docs.mdanalysis.org/
• GROMACS（替代 MD 引擎）：https://manual.gromacs.org/
• NAMD（替代方案）：https://www.ks.uiuc.edu/Research/namd/
• CHARMM-GUI（基于 Web 的体系构建器）：https://charmm-gui.org/
• AmberTools（免费 Amber 工具）：https://ambermd.org/AmberTools.php
• OpenMM 论文：Eastman P et al. (2017) PLOS Computational Biology. PMID: 28278240
• MDAnalysis 论文：Michaud-Agrawal N et al. (2011) J Computational Chemistry. PMID: 21500218