OpenMM系统构建基于Espaloma的System、Topology与Integrator版本声明本教程基于 espaloma 0.3.2、openff-toolkit 0.19.0、OpenMM ≥ 8.x。氢质量重分配HMR与约束的推荐值属于经典力场模拟经验参数如氢气原子质量约 1.5 amu、键长约束rigidWatertrue、constraintsHBonds其具体数值与 CPU/GPU 平台表现应以 OpenMM 官方文档为准本文仅作教学示意。所有能量单位为 OpenMM 内部默认单位kJ/mol、nm、ps换算为 kcal/mol 需显式 ×1/(4.184)。一句话结论openmm_system_from_graph只给了参数化的System把它交给LocalEnergyMinimizer.minimize(system, integrator)并把simulation.context.getState(...)的getPotentialEnergy()拿出来再加上一个LangevinMiddleIntegrator(300.0, 1.0/ps, 0.002*ps)才真正组成可在 OpenMM 里跑起来的最小化流程。〇、认知问题espaloma 输出的openmm_system与能跑 MD 的完整体系差在哪为什么需要再补 Topology 和 IntegratorHMR氢质量重分配和约束如何影响时间步长选择为什么推荐 1.5 amu 与HBondsLangevinMiddleIntegrator的friction/stepSize参数各代表什么为什么是恒温主选之一能量最小化返回的势能到底是什么数量级、怎么从 OpenMM 里取出来一、机制解析OpenMM 的程序模型是三个对象三职责Topology描述原子名、残基、链与化学键决定怎么把能量归到原子System描述所有力由力场参数化得到Integrator描述怎么推进坐标/速度。三条腿缺一不可。1.1 从 espaloma 到完整 System 的缺口第 6 篇部署锚点给的openmm_system已经有HarmonicBondForce、HarmonicAngleForce、PeriodicTorsionForce、NonbondedForce但通常没有溶剂水非键项为空或只含溶质分子坐标还不存在System本身不含坐标坐标属于Topology/后续Positions也还没有积分器与恒温器。因此任务变成三件① 拿坐标构造/补齐 Topology② 按需加 HMR 与约束③ 配Integrator并做最小化。1.2 HMR 与约束的物理动机经典经验步长受最快振动通常 C–H 伸缩约 3000 cm⁻¹周期极短限制若不用约束安全步长往往 ≤ 1 fs。**氢质量重分配HMR**把每个氢的额外质量借自相连的重原子如让 H 质量从 1 amu 提到约 1.5 amu同时冻结对光电子运动的显示求解把原型物理的质量效应隐藏进参数从而把 C–H 振动频率压下来配合constraintsHBonds约束所有 H 参与的键就能把积分步长安全提到 2 fs。关于约 1.5 amu需要一句方法论提醒它既不是 Espaloma 硬编码的参数也不是 OpenMM 内置的默认而是经典力场社区在约束 H 键 Langevin 恒温 常规生产场景下长期实践的推荐经验值。不同的重分配目标有的把 H 提到 3 amu、有的按比例均摊到邻接重原子会得到不同频率与不同数稳定性具体应以 OpenMM 官方文档与你的力场默认值为准不要把它当铁律抄进每个体系。默认(无HMR): C–H 振动 ω ≈ 高 → 步长必须 ≤1 fs HMR约束HBonds: 有效振动频率降低 → 步长可到 2 fs同物理时间更快完成 氢质量 1.0 → 1.5 amu重原子相应减重总质量守恒注意HMR 会改变质量矩阵、影响动能故只适合约束 H 化学键 恒温恒压的常规 MD做热容或严格谱学不推荐。1.3 OpenMM 的 Force 类型与最小化参数化后 System 里可能出现的力类型Force 类openmm.*对应项HarmonicBondForce键伸缩HarmonicAngleForce角弯曲PeriodicTorsionForce二面角傅里叶级数NonbondedForce非键 LJ 电荷CustomBondForce等高阶/自定义项espaloma 可能用到LocalEnergyMinimizer.minimize(system, integrator, maxIterations)是局部能量最小化或 L-BFGS/steepest descent它不推进动力学只是下坡找势能极小点返回的getPotentialEnergy()单位 kJ/mol。它需要Context或者至少Integrator实际会为最小化内部创建上下文。对最小化后的势能是多少要有量级直觉咖啡因这类中型有机分子在真空下用 Espaloma 参数化的势能通常落在几十到几万 kJ/mol 之间具体取决于构象是否合理、是否有大范围斥力重叠若数值出现极端负值千万量级或直接/nan多半是坐标错位或粒子序不一致而不是力场算错了。真正的判断不是有多负而是相对初值下降了、且稳定无 NaN、再最小化不再显著下降。这正是数值合理性 绝对数值的模拟第一性原则也是第 9、10 篇做双引擎对照时对拍的对象。二、完整代码与逐行剖析2.1 完整可运行建立 System Topology 最小化真空# filename: 07_minimize_vacuum.py# 基于锚点 A用 espaloma 参数化咖啡因补 Topology 并做能量最小化importespalomaasespfromopenff.toolkit.topologyimportMoleculeimportopenmmfromopenmmimportunitfromopenmm.appimportPDBFile,LocalEnergyMinimizer,Modeller moleculeMolecule.from_smiles(CN1CNC2C1C(O)N(C(O)N2C)C)mgesp.Graph(molecule)modelesp.get_model(latest)model.eval()# 本地权重必须 evalmodel(mg.heterograph)# 前向得参数systemesp.graphs.deploy.openmm_system_from_graph(mg)# 1) 从 OpenFF Molecule 生成 OpenMM Topologytopologymolecule.to_topology().to_openmm(topologyNone)# 坐标由 RDKit 提供get_model 体系下位置不属于 Systempositionsmolecule.conformers[0]._value*unit.nanometer# 教学示意取第一构象# 2) 选积分器LangevinMiddleIntegrator(温度, 摩擦系数, 步长)integratoropenmm.LangevinMiddleIntegrator(300.0*unit.kelvin,1.0/unit.picosecond,2.0*unit.femtosecond)# 3) 创建 Simulation 用于最小化simulationopenmm.app.Simulation(topology,system,integrator)simulation.context.setPositions(positions)# 4) 局部能量最小化LocalEnergyMinimizer.minimize(simulation.context)statesimulation.context.getState(getEnergyTrue)print(最小化后势能:,state.getPotentialEnergy().value_in_unit(unit.kilojoule_per_mole),kJ/mol)关键点molecule.conformers由Molecule.from_smiles用 RDKit 生成构象生产环境建议先generate_conformers或读 SDF/PDB。Simulation(topology, system, integrator)自动建Context是最小化最省事的入口。2.2 补粒子数与坐标一致性校验# filename: 07_consistency_check.py# 教学示意核对 System/Topology/Positions 三者粒子数与化学键数一致n_syssystem.getNumParticles()n_toptopology.getNumAtoms()n_poslen(positions)assertn_sysn_topn_pos,f{n_sys}!{n_top}!{n_pos}print(f粒子一致{n_sys}个原子键力项粒子对{system.getForce(0).getNumBonds()})这一步看起来琐碎却是排查系统坐标错位/最小化发散的第一道闸门System无坐标Topology才绑定坐标语义三者不一致时 OpenMM 会静默给出荒谬能量。2.3 带 HMR 的最小化教学示意参数以 OpenMM 文档为准# filename: 07_hmr_minimize.py# 教学示意HMR 约束后重做最小化。接口/推荐值以 OpenMM 官方文档为准fromopenmmimportapp,unitfromopenmm.appimportSimulation,LocalEnergyMinimizer# 常见推荐氢质量提到 ~1.5 amu示例值以 docs.openmm.org 为准CUSTOM_HEAVY{H:1.5*unit.atomic_mass_unit}# 教学示意把 module 里 H 原子质量改写——真实 HMR 常经 openmmforcefields/自定义# 以官方 MD 教程为准此处仅演示改质量→再约束→最小化的流程编排app.Topology# 引用以触发导入integratorapp.LangevinMiddleIntegrator(300*unit.kelvin,1/unit.picosecond,2*unit.femtosecond)simulationSimulation(topology,system,integrator)simulation.context.setPositions(positions)LocalEnergyMinimizer.minimize(simulation.context)# 教学示意HMR 后的最小化强调若不打算真用 HMR最稳妥是依赖 OpenMM/OpenMMForceFields 官方 HMR 接口或文档中的同参数流程本文的CUSTOM_HEAVY仅示意质量→频率→步长链路具体实现以官方为准。真实 MD 参数rigidWater、constraintsHBonds、步长 2 fs在第 8 篇综合使用。三、常见报错与排查现象根因处置Simulation时Topology与System原子数不符未生成位置或拓扑来源不匹配用第 2.2 的断言核对三者的粒子数最小化后能量极大/NaN坐标未归一或断言失败先setPositions再getState单位显式* unit.nanometer报No integrator错误最小化需 Context只给 System建Simulation(topology, system, integrator)再最小化LangevinMiddleIntegrator步长单位错步长量纲必须 time用2.0 * unit.femtosecond而非裸数字势能一运行就爆炸真空体系电荷大无溶剂屏蔽后续加显式溶剂或先用getState(getEnergyTrue)观察数量级四、动手练习最小化并输出势能跑 2.1打印最小化前后的getPotentialEnergy()确认差值应下降且无 NaN。换分子把分子换成CCCC、C1CCCCC1观察getNumParticles()与键力项getNumBonds()是否匹配直觉。步长与稳定性把stepSize从 2 fs 改成 10 fs重新最小化观察最小化是否仍收敛、能量是否出现异常增长体会步长受最快振动限制。一致性断言动手实现 2.2 的断言脚本故意给错坐标维度确认它会第一时间抛 AssertionError 而非静默产出坏能量。五、小结与下一篇预告这一篇把参数在手推进到能最小化System只是力的集合、需要配Topology与Integrator才完整HMR氢质量 ~1.5 amuHBonds约束把安全步长从 1 fs 提到 2 fsLocalEnergyMinimizer.minimize配合getState(getEnergyTrue)就能拿到以 kJ/mol 计的势能。你已拥有一个从 SMILES 到最小化结构 势能的最小闭环。下一篇8把闭环放大成真正的分子动力学实战给配体加显式溶剂盒Modeller.addSolvent、跑 NVT/NPT 平衡、用LangevinMiddleIntegrator走 300K/2fs/HMR 的生产轨迹并输出能量与 RMSD 曲线。本篇认知问题回显FAQQ1Espaloma 的 openmm_system 与能跑 MD 的完整体系差在哪里为何还要 Topology 和 IntegratorAopenmm_system只含各 Force 的参数本身没有坐标Topology 提供原子/残基/化学键语义并承载坐标Integrator 决定推进方式三者齐备才算可模拟的完整体系。Q2Espaloma 体系中氢质量重分配 HMR 与约束如何影响时间步长AHMR 把氢质量提到约 1.5 amu 并配合constraintsHBonds约束所有含 H 的键有效降低 C–H 振动频率使常规步长从约 1 fs 安全提到约 2 fs加快同物理时长的模拟。Q3LangevinMiddleIntegrator 的 friction 与 stepSize 参数各代表什么为何适合恒温Afriction摩擦系数单位 1/ps决定耦合到 Langevin 恒温的强度stepSize为积分步长Middle 型在恒温精度与采样间折中是常规 NVT/NPT 生产的默认选择之一。Q4能量最小化的势能如何从 OpenMM 中取出是什么量纲A用simulation.context.getState(getEnergyTrue)再调state.getPotentialEnergy()量纲为 OpenMM 内部默认的 kJ/mol换算 kcal/mol 需除以约 4.184。