OpenMM 跑 10–100 ns MD：小分子复合物模拟流程

OpenMM 是一个高性能的分子模拟引擎，最大的特点是完整的 Python API——这让它容易嵌入到自动化的计算流程中，而不像传统 MD 软件那样需要拼接命令行与配置文件。

准备：配体参数化

# 蛋白有现成的力场（Amber14、CHARMM36），
# 但【小分子必须单独生成参数】—— 这是 MD 流程中最容易出错的一步

pip install openmm openff-toolkit openmmforcefields pdbfixer mdtraj

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

# 从 SDF 读取配体（必须有正确的三维结构、氢与键级）
ligand = Molecule.from_file("ligand.sdf")

# 生成 OpenFF 力场模板
smirnoff = SMIRNOFFTemplateGenerator(
    molecules=ligand,
    forcefield="openff-2.2.0.offxml",
)

forcefield = app.ForceField("amber14-all.xml", "amber14/tip3p.xml")
forcefield.registerTemplateGenerator(smirnoff.generator)

# 【关键检查】：
#   1) 配体的质子化态是否正确？（见 105）
#   2) 立体化学是否正确？
#   3) 电荷是否正确（净电荷）？
#   4) 是否有力场覆盖不了的基团？
#      → 硼、硅、罕见金属配位等可能失败
#
# 替代方案：
#   GAFF2 + AM1-BCC（用 espaloma 或 antechamber）
#   from openmmforcefields.generators import GAFFTemplateGenerator
#   gaff = GAFFTemplateGenerator(molecules=ligand)

完整的模拟脚本

from openmm import app, unit, LangevinMiddleIntegrator, MonteCarloBarostat
import openmm as mm
from openff.toolkit import Molecule, Topology
from openmmforcefields.generators import SMIRNOFFTemplateGenerator
from pdbfixer import PDBFixer
import sys

# ---- 1) 准备蛋白 ----
fixer = PDBFixer(filename="protein.pdb")
fixer.findMissingResidues()
fixer.findMissingAtoms()
fixer.addMissingAtoms()
fixer.addMissingHydrogens(pH=7.4)
app.PDBFile.writeFile(fixer.topology, fixer.positions,
                      open("protein_fixed.pdb", "w"))

# ---- 2) 合并蛋白与配体 ----
ligand = Molecule.from_file("ligand.sdf")
lig_top = ligand.to_topology().to_openmm()
lig_pos = ligand.conformers[0].to_openmm()

protein_pdb = app.PDBFile("protein_fixed.pdb")
modeller = app.Modeller(protein_pdb.topology, protein_pdb.positions)
modeller.add(lig_top, lig_pos)

# ---- 3) 力场 ----
smirnoff = SMIRNOFFTemplateGenerator(molecules=ligand)
forcefield = app.ForceField("amber14-all.xml", "amber14/tip3p.xml")
forcefield.registerTemplateGenerator(smirnoff.generator)

# ---- 4) 溶剂化 ----
modeller.addSolvent(forcefield,
                    model="tip3p",
                    padding=1.0 * unit.nanometer,   # 边界距离
                    ionicStrength=0.15 * unit.molar,
                    positiveIon="Na+", negativeIon="Cl-")
print(f"体系原子数: {modeller.topology.getNumAtoms()}")

# ---- 5) 创建 System ----
system = forcefield.createSystem(
    modeller.topology,
    nonbondedMethod=app.PME,                    # 长程静电用 PME
    nonbondedCutoff=1.0 * unit.nanometer,
    constraints=app.HBonds,                     # 约束含氢键长 → 可用 2 fs 步长
    rigidWater=True,
    hydrogenMass=1.5 * unit.amu,                # 氢质量重分配 → 可用 4 fs
)
system.addForce(MonteCarloBarostat(1 * unit.bar, 300 * unit.kelvin, 25))

# ---- 6) 积分器与模拟对象 ----
integrator = LangevinMiddleIntegrator(
    300 * unit.kelvin,
    1.0 / unit.picosecond,        # 摩擦系数
    4.0 * unit.femtoseconds,      # 配合氢质量重分配
)
platform = mm.Platform.getPlatformByName("CUDA")
simulation = app.Simulation(modeller.topology, system, integrator, platform)
simulation.context.setPositions(modeller.positions)

# ---- 7) 能量最小化 ----
print("最小化前能量:",
      simulation.context.getState(getEnergy=True).getPotentialEnergy())
simulation.minimizeEnergy(maxIterations=5000)
print("最小化后能量:",
      simulation.context.getState(getEnergy=True).getPotentialEnergy())

# ---- 8) 平衡（带位置约束，逐步释放）----
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 modeller.topology.atoms():
    if atom.residue.name not in ("HOH", "NA", "CL") and atom.element.symbol != "H":
        restraint.addParticle(atom.index, modeller.positions[atom.index])
system.addForce(restraint)
simulation.context.reinitialize(preserveState=True)

simulation.context.setVelocitiesToTemperature(300 * unit.kelvin)
for k in [100.0, 50.0, 10.0, 2.0, 0.0]:
    simulation.context.setParameter("k", k * unit.kilojoules_per_mole / unit.nanometer**2)
    simulation.step(50000)        # 每级 200 ps
    print(f"平衡: k = {k}")

# ---- 9) 生产模拟 ----
simulation.reporters.append(
    app.DCDReporter("trajectory.dcd", 5000))        # 每 20 ps 存一帧
simulation.reporters.append(
    app.StateDataReporter("log.csv", 5000, step=True, time=True,
                          potentialEnergy=True, temperature=True,
                          volume=True, density=True, speed=True))
simulation.reporters.append(
    app.StateDataReporter(sys.stdout, 50000, step=True,
                          temperature=True, speed=True,
                          remainingTime=True, totalSteps=12500000))

simulation.step(12500000)        # 4 fs × 12.5M = 50 ns

simulation.saveState("final.xml")
print("完成")

关键参数的选择

常见问题

• 「No template found for residue」：最常见的错误。力场找不到某个残基的参数——通常是配体、非标准残基或辅因子。检查配体是否正确注册了模板生成器；
• 模拟爆炸（NaN 能量）：初始结构有严重原子重叠。先充分最小化，且平衡时用位置约束逐步释放；
• 温度或密度不稳：平衡时间不足，延长 NPT 平衡；
• 速度太慢：确认用了 CUDA 平台（platform.getName()）；开启氢质量重分配；检查是否用了 mixed 精度；
• 配体参数化失败：检查 SDF 的化学正确性；某些基团（硼、硅、特殊金属配位）OpenFF 可能不支持。

OpenMM 的独特优势

• 自定义力：CustomExternalForce、CustomBondForce 等可以用字符串表达式定义任意的势能项——做伞形采样、位置约束、自定义偏置都很方便；
• 与 Python 生态无缝：可以在模拟循环中插入任意 Python 代码，实时分析或调整；
• 易于嵌入流程：不需要生成中间配置文件，整个流程是一个 Python 脚本；
• OpenMM Setup：官方提供的图形界面工具，可以生成脚本模板——新手的好起点。

延伸资源

• MD 入门：123《分子动力学 MD 入门》；轨迹分析：124《RMSD 与 RMSF》、129《MDAnalysis 分析轨迹》；
• 力场参数化：202《OpenFF Toolkit》；结构准备：105《PDBFixer》；MM/GBSA：126《MM/GBSA》。