企业官网建设流程全解析

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

版本声明块

• 工具/软件：PLUMED 2.9（METAD/COORDINATIONNUMBER等动作、命令行plumed sum_hills/plumed driver）、Python 3.9+（numpy）、OpenMM 8.x（承接第 08 篇的 HILLS/COLVAR 产物）。
• 语言/环境：PLUMED 输入.dat+ Bash 命令行 + Python。单位：PLUMED 内部距离 nm、角度弧度，能量 kJ/mol。
• 本文目标：从"能跑一个 metaD"升级到"能设计合理 CV 并自己重建 FES"——读 HILLS、用 sum_hills/driver 还原自由能面、找局部极小。

一句话结论：PLUMED 的多种 CV 都由一行"标签: 动作 参数"算出来——DISTANCE ATOMS=1,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 找局部极小（凹度判定）。

〇、本篇要解决的认知问题

1. 距离/角/二面角之外，为什么还要配位数（COORDINATIONNUMBER）和 SASA 这两种"非几何" CV？
2. METAD ARG SIGMA HEIGHT PACE BIASFACTOR TEMP每个关键字到底管什么，不写会怎样？
3. HILLS 与 COLVAR 两个文件的列结构到底长什么样，怎么区分"偏置台账"和"轨迹账本"？
4. plumed sum_hills与plumed driver各在什么场合用，--kt/--mintozero各有什么用？
5. 拿到fes.dat（CV, 自由能两列）后，怎么用 Python 找局部极小并得自由能差？

一、机制解析

1.1 CV 类型与适用场景

PLUMED 的 CV 是一行"标签: 动作 参数"的声明，动作输出一个或多个标量供后续偏置/打印引用。除几何类外，配位数与 SASA 更贴近"结合/解开"这类事件：

COORDINATIONNUMBER 特别适合蛋白-配体体系描述"口袋里发生了几处接触"：它把"近"与"远"用 R_0 处的平滑阶梯函数过渡，天然是连续可微的反应坐标；SASA 则适合看溶剂暴露/去溶剂化——直接对应结合过程中的熵/溶剂重排。选 CV 的铁律仍见第 07 篇：CV 必须横穿关心的垒。

1.2 metad 与 well-tempered 的关键字

# 普通 metad：metad:METAD ARG=d SIGMA=0.05 HEIGHT=1.2 PACE=500 LABEL=metad# 良温 metad：metad:METAD ARG=d SIGMA=0.05 HEIGHT=1.2 PACE=500 BIASFACTOR=10 TEMP=300 LABEL=metad

• ARG：传给它加大势的 CV 标签（可多个，多维 metaD）。
• SIGMA：高斯宽度（与 CV 同单位）；HEIGHT：单次高斯高度；PACE：沉积周期（步）。
• BIASFACTOR：良温 γ；必须同时给TEMP。二者缺一会报"需良温温度"类错误（以官方文档为准）。
• 未给BIASFACTOR即普通版本，高度恒定、会持续增长（第 08 篇已论证）。

1.3 HILLS/COLVAR 的列格式（记账实物）

HILLS（METAD 自动写，每 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 则 N+1 列）、高度、宽度。良温下高度列呈衰减趋势。

COLVAR（PRINT 写，每 STRIDE 一行）：

#! FIELDS time d metad.bias 100.00 0.3221 -14.112 ...

即：时间、各 CV 值、当前总偏置势。

1.4 重建 FES：sum_hills 与 driver 分工

• sum_hills：累积偏置的负像恰是 FES。正确用法关键在两个选项：
  • --kt <kBT>：良温版必须传精确的 kBT（单位 kJ/mol，kT=kB×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 ATOMS=10,18a:ANGLE ATOMS=10,18,25t:TORSION ATOMS=7,10,18,25# 二、配位 CV：配体原子 1-20 与口袋残基原子 21-60 的接触数（R_0=0.35 nm 阈值）c:COORDINATIONNUMBER SPECIES=1-20 SPECIES2=21-60 R_0=0.35# 三、SASA CV：口袋侧链原子集s:SASA ATOMS=31-40# 四、良温 metad 偏置在这一个 CV（c）上：mt:METAD ARG=c SIGMA=0.2 HEIGHT=1.2 PACE=500 BIASFACTOR=10 TEMP=300 LABEL=mt# 五、打印：COLVAR 每 100 步记一次各 CV 值与偏置势PRINT ARG=d,a,t,c,s,mt.bias STRIDE=100 FILE=COLVAR

要点：SPECIES=1-20是原子序号区间；COORDINATIONNUMBER 的 R_0 单位 nm（0.35 nm≈3.5 Å＝最近接触距离）。此文件可在 OpenMM 里用PlumedForce(open("input.dat").read())直接载入（承接第 08 篇）。

2.2 Bash：sum_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 ATOMS=10,18# PRINT ARG=d STRIDE=1 FILE=proj.dat

逐要点：--kt 2.494必须与模拟温度一致（300 K）；--mintozero让全局极小对齐到 0；driver用--mf_xtc/--mf_pdb指定轨迹、--traj_stride指定步距、--plumed指向只算 CV 的输入文件。

2.3 Python：读取 fes.dat 找局部极小并算自由能差

importnumpyasnp# 读 FES：两列（CV, F）；已 mintozero 时 F 最小为 0data=np.loadtxt("fes.dat")cv,F=data[:,0],data[:,1]# 找局部极小：某点 F 同时小于左右邻居（用连续三点差分判别）dF=np.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）。

三、常见报错与排查

四、动手练习

练习 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 未收敛。

练习 3（driver 复核判据）：用driver_input.dat对traj.xtc前 100 帧批量算 d。判据：proj.dat的 d 值分布与 COLVAR 中 d 列分布一致（两套笔终给同一 CV 坐标，验证笔误/标签错配）。

五、小结与下一篇预告

本篇把 PLUMED 用透：五种代表性的 CV（DISTANCE/ANGLE/TORSION/COORDINATIONNUMBER/SASA）、metad/良温关键字与记账格式（HILLS/COLVAR）、以及两条重建 FES 的路——sum_hills（偏置负像转 FES，需 --kt/–mintozero）和driver（固定轨迹批量算 CV）；并用 Python 定位局部极小、读相对深度。至此，你已具备完整闭环：跑良温 metaD → 重建 FES → 找稳定构象。

第 10 篇预告：《复制交换与伞形采样》，将给出另两条互补的增强/自由能策略——REMD（温度副本交换）与伞形采样（沿 CV 分窗 + WHAM/PyMBAR 汇兑），并给出三方法适用场景对比表。

本篇认知问题回显（FAQ）

• Q1：为什么还要配位数和 SASA 这类"非几何"CV？
• A1：它们直接编码"接触数/溶剂暴露"这类对结合/去溶剂化事件本质的反应坐标，比单纯距离更接近物理过程，且更平滑可微。
• Q2：METAD 各关键字管什么？
• A2：ARG 指定偏置 CV；SIGMA 高斯宽度、HEIGHT 高度、PACE 沉积周期；BIASFACTOR 良温因子配 TEMP，经衰减使 HILLS 收敛；不给 BIASFACTOR 即普通版本、高度恒定增长。
• Q3：HILLS/COLVAR 列结构？
• A3：HILLS 每 PACE 一行：时间、各 CV 中心、高度、宽度；COLVAR 每 STRIDE 一行：时间、各 CV 值、当前偏置势。前者是偏置台账，后者是轨迹账本。
• Q4：sum_hills 与 driver 各何时用？
• A4：sum_hills 把累积偏置负像变成 FES（良温须 --kt、–mintozero 对齐最小）；driver 不跑动力学、对固定轨迹批量算 CV，用于独立投影验证。
• Q5：怎么用 Python 找局部极小？
• A5：np.loadtxt读到两列 CV/F，用np.gradient求导、找斜率由负转正的点；ΔF 为谷间落差，但须确认 FES 已收敛且重建参数正确。