AI增强构象采样教程（9）：PLUMED 深度——CV 定义与 FES 重建

AI增强构象采样教程9PLUMED 深度——CV 定义与 FES 重建版本声明块工具/软件PLUMED 2.9METAD/COORDINATIONNUMBER等动作、命令行plumed sum_hills/plumed driver、Python 3.9numpy、OpenMM 8.x承接第 08 篇的 HILLS/COLVAR 产物。语言/环境PLUMED 输入.dat Bash 命令行 Python。单位PLUMED 内部距离 nm、角度弧度能量 kJ/mol。本文目标从能跑一个 metaD升级到能设计合理 CV 并自己重建 FES——读 HILLS、用 sum_hills/driver 还原自由能面、找局部极小。一句话结论PLUMED 的多种 CV 都由一行标签: 动作 参数算出来——DISTANCE ATOMS1,2、ANGLE/TORSION、COORDINATIONNUMBER SPECIES… SPECIES2… R_0…、SASA ATOMS…——再被METAD ARG… SIGMA… HEIGHT… PACE… BIASFACTOR… TEMP…偏置重构时plumed sum_hills --hills HILLS --kt kBT --mintozero --bin n --outfile fes.dat把累积偏置的负像变成 FES或用plumed driver--plumed plumed.dat --mf_* 轨迹--traj_stride对任意构象批量算 CV最后 Pythonloadtxt读fes.dat沿 CV 找局部极小凹度判定。〇、本篇要解决的认知问题距离/角/二面角之外为什么还要配位数COORDINATIONNUMBER和 SASA 这两种非几何 CVMETAD ARG SIGMA HEIGHT PACE BIASFACTOR TEMP每个关键字到底管什么不写会怎样HILLS 与 COLVAR 两个文件的列结构到底长什么样怎么区分偏置台账和轨迹账本plumed sum_hills与plumed driver各在什么场合用--kt/--mintozero各有什么用拿到fes.datCV, 自由能两列后怎么用 Python 找局部极小并得自由能差一、机制解析1.1 CV 类型与适用场景PLUMED 的 CV 是一行标签: 动作 参数的声明动作输出一个或多个标量供后续偏置/打印引用。除几何类外配位数与 SASA 更贴近结合/解开这类事件动作语法要点场景/单位DISTANCEATOMS1,2两原子欧氏距离距离单位 nmANGLEATOMS1,2,3弧度TORSIONATOMS1,2,3,4弧度二面角COORDINATIONNUMBERSPECIES… SPECIES2… R_0…键基函数平滑0~N 无单位判断是否接触SASAATOMS…溶剂可及表面积ŲCOORDINATIONNUMBER 特别适合蛋白-配体体系描述口袋里发生了几处接触它把近与远用 R_0 处的平滑阶梯函数过渡天然是连续可微的反应坐标SASA 则适合看溶剂暴露/去溶剂化——直接对应结合过程中的熵/溶剂重排。选 CV 的铁律仍见第 07 篇CV 必须横穿关心的垒。1.2 metad 与 well-tempered 的关键字# 普通 metadmetad:METAD ARGd SIGMA0.05 HEIGHT1.2 PACE500 LABELmetad# 良温 metadmetad:METAD ARGd SIGMA0.05 HEIGHT1.2 PACE500 BIASFACTOR10 TEMP300 LABELmetadARG传给它加大势的 CV 标签可多个多维 metaD。SIGMA高斯宽度与 CV 同单位HEIGHT单次高斯高度PACE沉积周期步。BIASFACTOR良温 γ必须同时给TEMP。二者缺一会报需良温温度类错误以官方文档为准。未给BIASFACTOR即普通版本高度恒定、会持续增长第 08 篇已论证。1.3 HILLS/COLVAR 的列格式记账实物HILLSMETAD 自动写每 PACE 一行#! FIELDS time metad.d metad.height metad.sigma 0.00 0.3110 1.2000 0.0500 1.00 3.4123 1.1922 0.0500 ...即时间、CV 中心第 N 个 ARG 则 N1 列、高度、宽度。良温下高度列呈衰减趋势。COLVARPRINT 写每 STRIDE 一行#! FIELDS time d metad.bias 100.00 0.3221 -14.112 ...即时间、各 CV 值、当前总偏置势。1.4 重建 FESsum_hills 与 driver 分工sum_hills累积偏置的负像恰是 FES。正确用法关键在两个选项--kt kBT良温版必须传精确的 kBT单位 kJ/molkTkB×T否则 sum_hills 默认按普通版本求极限体系温度无论如何也要对--mintozero把 FES 整体上移使最小值为 0便于读局部极小与相对深。其他--bin n网格网格点数、--stride分段重建、--outfile fes.dat。driver不跑动力学的计算器——读入固定轨迹traj.pdb/traj.xtc对每一帧批量算可选 CV 并输出。典型用途是把已采样的构象投影到一个独立 CV 上验证如 band/contact而不必重跑偏置。fes.dat 输出两列CV 值 自由能假定已经 mintozero。找极小即找此一维函数的局部凹点。二、完整代码与逐行剖析2.1 PLUMED 输入文件全文input.dat含多 CV 定义与良温 metad# 一、几何 CV配体原子 1-2 的距离与二面角索引 1 起始d:DISTANCE ATOMS10,18a:ANGLE ATOMS10,18,25t:TORSION ATOMS7,10,18,25# 二、配位 CV配体原子 1-20 与口袋残基原子 21-60 的接触数R_00.35 nm 阈值c:COORDINATIONNUMBER SPECIES1-20 SPECIES221-60 R_00.35# 三、SASA CV口袋侧链原子集s:SASA ATOMS31-40# 四、良温 metad 偏置在这一个 CVc上mt:METAD ARGc SIGMA0.2 HEIGHT1.2 PACE500 BIASFACTOR10 TEMP300 LABELmt# 五、打印COLVAR 每 100 步记一次各 CV 值与偏置势PRINT ARGd,a,t,c,s,mt.bias STRIDE100 FILECOLVAR要点SPECIES1-20是原子序号区间COORDINATIONNUMBER 的 R_0 单位 nm0.35 nm≈3.5 Å最近接触距离。此文件可在 OpenMM 里用PlumedForce(open(input.dat).read())直接载入承接第 08 篇。2.2 Bashsum_hills 重建 FES driver 批量算 CV# sum_hills良温重建kT(300K)2.494 kJ/mol--mintozero 使最小值为 0plumed sum_hills--hillsHILLS--kt2.494--bin100--mintozero--outfilefes.dat# 查看 fes.dat 头 5 行CV, F 两列head-5fes.dat# driver对轨迹批量算 CV不跑动力学plumed driver--plumeddriver_input.dat--mf_xtctraj.xtc--traj_stride1# driver_input.dat 只含 CV 定义与 PRINT去掉 METAD 偏置如# d: DISTANCE ATOMS10,18# PRINT ARGd STRIDE1 FILEproj.dat逐要点--kt 2.494必须与模拟温度一致300 K--mintozero让全局极小对齐到 0driver用--mf_xtc/--mf_pdb指定轨迹、--traj_stride指定步距、--plumed指向只算 CV 的输入文件。2.3 Python读取 fes.dat 找局部极小并算自由能差importnumpyasnp# 读 FES两列CV, F已 mintozero 时 F 最小为 0datanp.loadtxt(fes.dat)cv,Fdata[:,0],data[:,1]# 找局部极小某点 F 同时小于左右邻居用连续三点差分判别dFnp.gradient(F,cv)minima_idx[]foriinrange(1,len(F)-1):ifdF[i-1]0anddF[i]0:# 斜率由负转正 局部极小凹谷minima_idx.append(i)foriinminima_idx:print(f局部极小 CV{cv[i]:.3f}, F{F[i]:.3f}kJ/mol)ifminima_idx:Fm[F[i]foriinminima_idx]# 全局最小最低自由能差 次小与最小之差铁律 5需 FES 已收敛才有意义print(全局极小 F%.3f, 相邻谷落差 ΔF≈%.2f kJ/mol%(min(Fm),sorted(Fm)[1]-min(Fm)iflen(Fm)1else0.0))逐要点np.gradient对一维 FES 求导用差分符号由负转正判局部极小ΔF 即构象间稳定差但前提是 sum_hills 用的 --kt 与重建网格正确、且 HILLS 已衰减收敛——否则 ΔF 值无意义呼应铁律 5。三、常见报错与排查现象根因解法COORDINATIONNUMBER报 SPECIES/SPECIES2 索引越界原子序号从 0 起写PLUMED 原子索引 1 起始区间区间写法1-20表示 1…20核查拓扑原子数sum_hills 输出警告温度与 --kt 不一致或 FES 全负未给正确--kt良温或未--mintozero良温版必须--kt kBT--mintozero上移全局最小到 0BIASFACTOR 给了但没给 TEMP → 报错良温需要温度计算衰减补TEMP300或去掉 BIASFACTOR 用普通版本driver 无输出/第一列全 NaN--plumed文件里还残留 METAD 偏置不该跑动力学driver 输入只保留 CV 定义 PRINT去掉 METADfes.dat 第一列不分bucket网格异常--bin过小/网格粒度不足增大--bin或用--outgrid自定义四、动手练习练习 1重建判据对第 08 篇的 HILLS 跑plumed sum_hills --hills HILLS --kt 2.494 --bin 100 --mintozero --outfile fes.dat。判据fes.dat第二列最小接近 0且曲线在 CV 出现 1 个以上凹谷。练习 2极小定位判据用 2.3 的 Python 读 fes.dat。判据脚本打印的局部极小数与直方图目测谷数一致且 ΔF 落在合理区间如 0–15 kJ/mol过大说明 FES 未收敛。练习 3driver 复核判据用driver_input.dat对traj.xtc前 100 帧批量算 d。判据proj.dat的 d 值分布与 COLVAR 中 d 列分布一致两套笔终给同一 CV 坐标验证笔误/标签错配。五、小结与下一篇预告本篇把 PLUMED 用透五种代表性的 CVDISTANCE/ANGLE/TORSION/COORDINATIONNUMBER/SASA、metad/良温关键字与记账格式HILLS/COLVAR、以及两条重建 FES 的路——sum_hills偏置负像转 FES需 --kt/–mintozero和driver固定轨迹批量算 CV并用 Python 定位局部极小、读相对深度。至此你已具备完整闭环跑良温 metaD → 重建 FES → 找稳定构象。第 10 篇预告《复制交换与伞形采样》将给出另两条互补的增强/自由能策略——REMD温度副本交换与伞形采样沿 CV 分窗 WHAM/PyMBAR 汇兑并给出三方法适用场景对比表。本篇认知问题回显FAQQ1为什么还要配位数和 SASA 这类非几何CVA1它们直接编码接触数/溶剂暴露这类对结合/去溶剂化事件本质的反应坐标比单纯距离更接近物理过程且更平滑可微。Q2METAD 各关键字管什么A2ARG 指定偏置 CVSIGMA 高斯宽度、HEIGHT 高度、PACE 沉积周期BIASFACTOR 良温因子配 TEMP经衰减使 HILLS 收敛不给 BIASFACTOR 即普通版本、高度恒定增长。Q3HILLS/COLVAR 列结构A3HILLS 每 PACE 一行时间、各 CV 中心、高度、宽度COLVAR 每 STRIDE 一行时间、各 CV 值、当前偏置势。前者是偏置台账后者是轨迹账本。Q4sum_hills 与 driver 各何时用A4sum_hills 把累积偏置负像变成 FES良温须 --kt、–mintozero 对齐最小driver 不跑动力学、对固定轨迹批量算 CV用于独立投影验证。Q5怎么用 Python 找局部极小A5np.loadtxt读到两列 CV/F用np.gradient求导、找斜率由负转正的点ΔF 为谷间落差但须确认 FES 已收敛且重建参数正确。