AI增强构象采样教程（8）：第一个增强采样实验——OpenMM+PLUMED 良温元动力学

AI增强构象采样教程8第一个增强采样实验——OpenMMPLUMED 良温元动力学版本声明块工具/软件OpenMM 8.x、openmm-plumed提供PlumedForce、PLUMED 2.9偏置内核、PDBFixer、Python 3.9力场 ff14SBtip3popenmmforcefields。语言/环境Pythonopenmm、openmm-plumed、numpy。能量单位OpenMM 默认 kJ/molPLUMED 高斯高度 HEIGHT 单位亦 kJ/mol可配置。本文目标完成零到一的良温元动力学运行环境、单 CV、PlumedForce 注入、HILLS 初步解读。为第 09 篇多 CV/FES 重建铺路。一句话结论一条pip install openmm openmm-plumed openmmforcefields补全然后createSystem(..., removeCMMotionFalse)后在system上addForce(PlumedForce(cv_script))其中 PLUMED 脚本用METAD ARGd SIGMA0.05 HEIGHT1.2 PACE500 BIASFACTOR10 TEMP300沉积高斯跑完生成HILLS每 PACE 一步一行时间、CV 中心、高度、宽度与COLVAR每步 CV 值 偏置——先做单 CV 的 well-tempered 实验才能看懂物理体积扩散而不至于被 HILLS 高度误导以为偏置在无脑增长。〇、本篇要解决的认知问题openmm-plumed 和 PLUMED 是什么关系PlumedForce为什么被 push 到system.addForce良温元动力学 vs 普通元动力学代码差别在哪BIASFACTOR 干什么用一个 CVRMSD 或距离PLUMED 脚本怎么写、METAD关键字怎么排SIGMA/HEIGHT 这两个参数分别控制什么选多少合适不瞎填HILLS和COLVAR都是偏置的台账它们每行每列分别代表什么、怎么初步看懂收敛一、机制解析1.1 openmm-plumedPLUMED 是库PlumedForce 是桥PLUMED 本身是一套独立的增强采样库也能以-plumed接入 GROMACS 等分子动力学引擎。openmm-plumed 提供一条openmm.PlumedForcePython 包名openmm-plumed——你只需在createSystem()后的System上加一个PlumedForce(plumed_input)OpenMM 的每一步积分后会自动调用 PLUMED 计算 CV、沉积高斯偏置并写HILLS/COLVAR。它把 PLUMED 的输入脚本作为字符串传来传去因此 CV 与偏置在 PLUMED 侧定义OpenMM 只负责力和积分。这是Python 优先二次开发接口的典型形态。1.2 良温 vs 普通多一行 BIASFACTOR 的差别普通元动力学的高度 W 固定会无界地填平整个 CV 空间长时间后偏置势发散良温元动力学well-tempered加BIASFACTORγ常用 5–20使高斯高度随已累积偏置按 e^{−V_bias/(kB·(γ−1)·T)} 衰减在未探区V_bias 小高度近满快速推进在已探区高度迅速趋零体系松松地停留在每个谷里统计。结果收敛到一个有限、有界且单调逼近真实 −F(s) 的态这是铁律 5监控 bias 收敛能成立的基础。代码层面两者仅差一行BIASFACTOR——但物理机制完全不同切勿把普通版本的 HILLS 直线增长当收敛。1.3 单 CV 的选择与定义对蛋白配体体系最省事的两个初选 CV距离 CV两个关键原子如配体某原子与口袋残基 O的DISTANCE——代码直观、单位 Å。RMSD CV配体相对参考构象的RMSD REFERENCEref.pdb TYPEOPTIMAL——单位 Å适合看整体重排。本篇以距离 CV 为主、RMSD 为注。定义在 PLUMED 输入里先label: 动作 参数算出一个标量再METAD ARGlabel ...把它当偏置参量。二、完整代码与逐行剖析2.1 环境安装命令行首次执行# 补全纯 MD 偏置依赖Python 3.9pipinstallopenmm openmm-plumed openmmforcefields pdbfixer# 校验 PlumedForce 可用python-cimport openmm; from openmm.app import * sSimulation; print(openmm ok)# 若要用 QM/BF 扩展再装 plumed 本体仅跑 metaD 时 openmm-plumed 内置即可2.2 完整可运行脚本单距离 CV 的良温元动力学importopenmmfromopenmmimportapp,unitfromopenmmimportPlumedForce# 由 openmm-plumed 注册进 openmmimportnumpyasnp# 1) 准备体系PDBFixer 清洗 力场ff14SBtip3ppdbapp.PDBFile(protein_lig.pdb)# 已清洗/加氢结构见第 05 篇forcefieldapp.ForceField(amber14-all.xml,amber14/tip3pfb.xml)# 大体系/配体建议用 openmmforcefields 的 GAFF2 参数化此处为可跑演示用标准力场# 2) 建 SystemremoveCMMotionFalse 保证质心不受 PLUMED 偏置判断干扰systemforcefield.createSystem(pdb.topology,nonbondedMethodapp.PME,constraintsapp.HBonds,rigidWaterTrue,removeCMMotionFalse,# 保留质心运动便于偏置按相对量设定)# 3) PLUMED 输入定义 1 个距离 CV 良温 METAD 偏置# 原子对取配体 1 号原子 与 口袋残基参考原子示例索引 10, 50需按你的体系改plumed_script # 一、标量 CV两个原子的距离单位 nmPLUMED 内部用 nm d: DISTANCE ATOMS10,50 # 二、良温元动力学偏置沉积高斯的四大参数 BIASFACTOR 良温因子 metad: METAD ARGd SIGMA0.05 HEIGHT1.2 PACE500 BIASFACTOR10 TEMP300 LABELmetad # 三、输出COLVAR 每 100 步记录一次 CV 值 当前偏置势 PRINT ARGd,metad.bias STRIDE100 FILECOLVAR # 4) 把 PLUMED 包进 OpenMM 的力场并加入 system关键一步plumed_forcePlumedForce(plumed_script)system.addForce(plumed_force)# 5) 建模拟器恒定温度复平衡常用 Langevin 中速阻力integratoropenmm.LangevinMiddleIntegrator(300*unit.kelvin,1.0/unit.picosecond,2*unit.femtosecond)platformopenmm.Platform.getPlatformByName(CUDA)# 无 GPU 可用 CPUsimulationapp.Simulation(pdb.topology,system,integrator,platform)simulation.context.setPositions(pdb.positions)simulation.minimizeEnergy()# 6) 运行localStep 一次 dt2fs跑 10000 步 ≈ 20 ps 演示正式跑长得多simulation.step(10000)# 7) 观察 HILLS偏置沉积记录最后一帧高度见 at 偏置势是否趋零hillsnp.loadtxt(HILLS,skiprows1)print(HILLS 行数%d, 末尾高度(W)列均值%.3f kJ/mol%(len(hills),hills[-20:,3].mean()))colvarnp.loadtxt(COLVAR,skiprows1)print(COLVAR 记录: 最后一帧 d(cv)%.4f nm, bias%.3f kJ/mol%(colvar[-1,1],colvar[-1,2]))逐行剖析PlumedForce来自openmm-plumed正确from openmm import PlumedForce使用。createSystem(... removeCMMotionFalse)偏置势只与内部坐标有关故关闭质心抑制以避免与元动力学冲突以官方文档为准。脚本里d: DISTANCE ATOMS10,50——原子序号是 1 起始的 PLUMED 索引务必对照拓扑核对单位 nmnet 两原子距离。METAD ARGd SIGMA0.05 HEIGHT1.2 PACE500 BIASFACTOR10 TEMP300SIGMA 与 CV 单位一致nm→0.05 nm ≈0.5 ÅHEIGHT 高度 kJ/molPACE 每 500 步沉积一个高斯BIASFACTOR10 良温因子TEMP 用于良温的指数衰减因子。PRINT ... FILECOLVARSTRIDE100控制轨迹记录频率HILLS则由 METAD 每 PACE 自动写。运行 10000 步 ×2fs20ps 只是演示收敛不够正式元动力学通常 ns–μs并按下文铁律 5 监控 HILLS 衰减。2.3 另一个初选 CVRMSD作为完成代码补充plumed_script_rmsd # RMSD CV相对参考结构ref.pdb 取配体/口袋 cv: RMSD REFERENCEref.pdb TYPEOPTIMAL metad: METAD ARGcv SIGMA0.15 HEIGHT1.2 PACE500 BIASFACTOR10 TEMP300 LABELmetad PRINT ARGcv,metad.bias STRIDE100 FILECOLVAR_RMSD plumed_force_rmsdPlumedForce(plumed_script_rmsd)三、常见报错与排查现象根因解法ImportError: No module named openmmplumed未安装openmm-plumed或名字不同pip install openmm-plumed确认导入名以官方 README 为准库注册后from openmm import PlumedForcePLUMED 报atom index X out of range原子序号越界或从 0 起写PLUMED 原子索引 1 起始ATOMS10,50对应拓扑第 10/50 个原子先复核拓扑报System无addForce或偏置无效PlumedForce 未正确加入 system 次序确保在createSystem之后、Simulation之前addForce多力场共存时检查NonbondedForce是否完整HILLS 高度列单调递增不衰减用了普通 METAD无 BIASFACTOR或体系未充分返访旧谷补BIASFACTOR并让返访时长 ≥ 谷内扩散时长偏置势数值异常大/发散、速度失控半时间步与 removeCMMotion 设置冲突或时间步过大检查removeCMMotionFalse与integrator步长必要时降 dt四、动手练习练习 1接线判据把 2.2 脚本换成你的蛋白配体、原子对选好跑 10000 步。判据HILLS文件生成且行数 ≈ 10000/500 20COLVAR行数 ≈ 10000/100 100。练习 2衰减观察判据把模拟加到 50000 步≈100 psnumpy算 HILLS 末尾 50 行的 W 列均值。判据良温下该均值相比前段应明显下降趋进口衰减证明体系返访旧谷而非单向冲开。练习 3对比实验判据删掉BIASFACTOR重新跑 50000 步。判据HILLS高度列不再衰减单调累积与练习 2 形成DID 对比直观凸显良温的价值。五、小结与下一篇预告本篇完成了第一个增强采样实验安装依赖、定义单距离/RMSD CV、把PlumedForce注入 System、跑良温元动力学、从 HILLS/COLVAR 初步判断收敛。关键落地是 SIGMA/HEIGHT/PACE/BIASFACTOR 四个参数的物理含义与单位以及良温保证 HILLS 衰减、普通版本持续增长的判别。这只是开始——单 CV 的 FES 尚未重建。第 09 篇预告《PLUMED 深度CV 定义与 FES 重建》将展开 PLUMED 的多种 CVDISTANCE/ANGLE/TORSION/COORDINATIONNUMBER/SASA、metad/well-tempered语法、HILLS/COLVAR格式并用plumed sum_hills/driver重建 FES、定位局部极小。本篇认知问题回显FAQQ1openmm-plumed 与 PLUMED 什么关系PlumedForce 为何 addForceA1PLUMED 是增强采样库openmm-plumed 提供PlumedForce把 PLUMED 输入脚本包成 OpenMM 的Force在createSystem后system.addForce积分器每一步自动算 CV 与偏置。Q2良温 vs 普通元动力学代码差别A2仅差一行BIASFACTORγ良温按 e^{−V_bias/(kB(γ−1)T)} 衰减高斯高度普通版本高度固定会持续增长。良温保证收敛到有限有界态。Q3单 CV 的 PLUMED 脚本怎么写A3先d: DISTANCE ATOMS10,50或cv: RMSD REFERENCEref.pdb算出 CV再METAD ARGd SIGMA0.05 HEIGHT1.2 PACE500 BIASFACTOR10 TEMP300最后PRINT ... FILECOLVAR。Q4SIGMA/HEIGHT 怎么选A4SIGMA 是高斯沿 CV 的宽度、单位与 CV 一致距离取 ≈0.05nmHEIGHT 是单次高度kJ/mol决定推挤强度PACE 是沉积周期三者共同决定填谷速度合理初值以官方文档为准。Q5HILLS/COLVAR 每列代表什么、怎么看收敛A5HILLS 每沉积周期一行时间、CV 中心、高度、宽度COLVAR 每 STRIDE 行时间、CV 值、偏置势。良温下 HILLS 高度列应衰减不衰减即偏置在单调累积、FES 未收敛。