OpenMM 教程：分子动力学模拟入门路线

OpenMM 的 Python API 让分子动力学模拟的入门门槛大幅降低——不需要拼接命令行与配置文件，整个流程就是一个 Python 脚本。这篇给出从零到跑通一个蛋白-配体复合物模拟的学习路线。

学习路线

阶段	目标	时间
1. 纯蛋白在水中	跑通最简单的流程	1～2 天
2. 理解每一步	最小化、平衡、生产的含义	2～3 天
3. 加入小分子	力场参数化（最难的一步）	3～5 天
4. 轨迹分析	RMSD、相互作用（见 124《RMSD 与 RMSF》、129《MDAnalysis 分析轨迹》）	2～3 天
5. 自定义力	约束、伞形采样	按需
6. 自由能计算	MM/GBSA、FEP（见 126《MM/GBSA》、127《FEP、RBFE 与 ABFE》）	按需

阶段 1：最简单的可运行例子

pip install openmm pdbfixer mdtraj
# 或 conda install -c conda-forge openmm pdbfixer mdtraj

# 验证安装与可用平台
python -m openmm.testInstallation
# 【应该看到 CUDA 或 OpenCL 平台】
# 只有 CPU 的话，模拟会慢几十倍

from openmm import app, unit, LangevinMiddleIntegrator
import openmm as mm
import sys

# 1) 读入蛋白
pdb = app.PDBFile("protein.pdb")

# 2) 力场
forcefield = app.ForceField("amber14-all.xml", "amber14/tip3p.xml")

# 3) 加氢与溶剂化
modeller = app.Modeller(pdb.topology, pdb.positions)
modeller.addHydrogens(forcefield, pH=7.4)
modeller.addSolvent(forcefield, model="tip3p",
                    padding=1.0 * unit.nanometer,
                    ionicStrength=0.15 * unit.molar)
print(f"体系原子数: {modeller.topology.getNumAtoms()}")

# 4) 创建 System
system = forcefield.createSystem(
    modeller.topology,
    nonbondedMethod=app.PME,              # 【长程静电必须用 PME】
    nonbondedCutoff=1.0 * unit.nanometer,
    constraints=app.HBonds,               # 约束含氢键长 → 2 fs 步长
)

# 5) 积分器
integrator = LangevinMiddleIntegrator(
    300 * unit.kelvin, 1 / unit.picosecond, 2 * unit.femtoseconds)

# 6) 模拟对象
simulation = app.Simulation(modeller.topology, system, integrator)
simulation.context.setPositions(modeller.positions)

# 7) 最小化
print("最小化中...")
simulation.minimizeEnergy()

# 8) 平衡
simulation.context.setVelocitiesToTemperature(300 * unit.kelvin)
simulation.step(50000)      # 100 ps

# 9) 生产
simulation.reporters.append(app.DCDReporter("traj.dcd", 5000))
simulation.reporters.append(app.StateDataReporter(
    sys.stdout, 5000, step=True, potentialEnergy=True,
    temperature=True, speed=True))
simulation.step(500000)     # 1 ns

# 【跑通这个，你就完成了第一个 MD 模拟】

阶段 2：理解每一步为什么

# ---- 为什么要最小化 ----
#   初始结构（尤其加了氢和水之后）可能有原子重叠
#   → 直接跑 MD 会因为巨大的力而「爆炸」
#   → 最小化消除这些不合理接触
#
#   【检查】：最小化后能量应该显著下降且为合理的负值

# ---- 为什么要分阶段平衡 ----
#   NVT（固定体积升温）：
#     让体系达到目标温度
#   NPT（固定压力）：
#     让密度达到平衡（约 1 g/cm³）
#
#   【带位置约束逐步释放】：
#     一开始约束蛋白重原子，只让水弛豫
#     然后逐步减小约束力常数
#     → 【避免蛋白结构在水未平衡时被破坏】
#
#   【常见错误】：跳过平衡直接生产
#     → 前面的轨迹不可用，甚至模拟不稳定

# ---- 关键参数的含义 ----
#
#   nonbondedMethod=PME
#     长程静电用 Ewald 求和
#     → 【必须用】：简单截断会造成严重伪影
#
#   constraints=HBonds
#     约束含氢的键长（它们振动最快）
#     → 让时间步长从 1 fs 提高到 2 fs
#
#   hydrogenMass=1.5*amu（氢质量重分配）
#     把重原子的一部分质量转给氢
#     → 【时间步长可到 4 fs，速度快近一倍】
#     → 不影响热力学性质
#
#   Langevin 积分器
#     加入摩擦与随机力，控温
#     摩擦系数 1/ps 是常用值
#
#   Barostat 频率 25 步
#     太频繁会影响动力学，太稀疏平衡慢

# ---- 怎么判断模拟是否正常 ----
#   □ 温度稳定在 300 K 附近（±3 K）
#   □ 密度稳定在约 1.0 g/cm³
#   □ 势能稳定（不持续下降或上升）
#   □ 【蛋白主链 RMSD 稳定在 1~3 Å】（见 124）
#   □ 没有 NaN 或异常大的能量

阶段 3：加入小分子（最难的一步）

# 【为什么难】：蛋白有现成力场，小分子必须参数化（见 014）

from openff.toolkit import Molecule
from openmmforcefields.generators import SMIRNOFFTemplateGenerator
from openmm import app

# 1) 读入配体（【必须是含正确键级与三维结构的 SDF】）
ligand = Molecule.from_file("ligand.sdf")

# 2) 创建模板生成器
smirnoff = SMIRNOFFTemplateGenerator(
    molecules=ligand, forcefield="openff-2.2.0.offxml")

# 3) 注册到力场
forcefield = app.ForceField("amber14-all.xml", "amber14/tip3p.xml")
forcefield.registerTemplateGenerator(smirnoff.generator)

# 4) 合并蛋白与配体
protein_pdb = app.PDBFile("protein_fixed.pdb")
modeller = app.Modeller(protein_pdb.topology, protein_pdb.positions)
modeller.add(ligand.to_topology().to_openmm(),
             ligand.conformers[0].to_openmm())

# 5) 之后与纯蛋白流程相同

# 【最常见的报错】：
#   "No template found for residue XXX"
#   → 说明力场找不到某个残基的参数
#   → 检查：
#     - 配体的模板生成器注册了吗？
#     - 有没有非标准残基、辅因子、金属？
#     - 结构是否已用 PDBFixer 修复？（见 105）
#
# 【第二常见的问题】：
#   模拟一开始就爆炸（NaN）
#   → 原因：初始结构有严重重叠
#   → 应对：
#     1) 增加最小化步数
#     2) 【用位置约束逐步释放的平衡流程】
#     3) 检查配体是否与蛋白原子重叠
#     4) 先用较小的时间步长（1 fs）跑一段

OpenMM 的独特能力：自定义力

# 【这是 OpenMM 相对其它引擎最实用的优势】
#   可以用字符串表达式定义任意势能项

import openmm as mm
from openmm import unit

# ---- 位置约束（平衡时常用）----
restraint = mm.CustomExternalForce(
    "k*periodicdistance(x, y, z, x0, y0, z0)^2")
restraint.addGlobalParameter(
    "k", 100.0 * unit.kilojoules_per_mole / unit.nanometer**2)
for p in ["x0", "y0", "z0"]:
    restraint.addPerParticleParameter(p)
for atom in topology.atoms():
    if atom.element.symbol != "H" and atom.residue.name not in ("HOH",):
        restraint.addParticle(atom.index, positions[atom.index])
system.addForce(restraint)

# 运行时调整强度
# simulation.context.setParameter("k", 10.0 * ...)

# ---- 距离约束（保持配体在口袋中）----
dist_restraint = mm.CustomCentroidBondForce(
    2, "0.5*k*max(0, distance(g1,g2) - d0)^2")
dist_restraint.addGlobalParameter(
    "k", 1000.0 * unit.kilojoules_per_mole / unit.nanometer**2)
dist_restraint.addGlobalParameter("d0", 0.5 * unit.nanometer)
dist_restraint.addGroup(ligand_atom_indices)
dist_restraint.addGroup(pocket_atom_indices)
dist_restraint.addBond([0, 1], [])
system.addForce(dist_restraint)

# ---- 用途 ----
#   - 平衡时保护结构
#   - 伞形采样（沿反应坐标加偏置）
#   - 保持配体不漂走（研究特定构象时）
#   - 实现自定义的增强采样方法

关键要点

学习顺序：纯蛋白 → 理解每步 → 加小分子，配体参数化是最难的一步；
「No template found」是最常见的报错——检查模板生成器是否注册、结构是否修复；
带位置约束逐步释放的平衡流程能避免模拟爆炸；
氢质量重分配 + 4 fs 步长让速度快近一倍；长程静电必须用 PME。

延伸资源

OpenMM 实战：128《OpenMM 跑 10–100 ns MD》；OpenFF 教程：014《OpenFF 教程》；MD 入门：123《分子动力学 MD 入门》；
轨迹分析：124《RMSD 与 RMSF》、129《MDAnalysis 分析轨迹》；结构准备：105《PDBFixer》。