AI增强构象采样教程（10）：复制交换与伞形采样——REMD 与自由能剖面

AI增强构象采样教程10复制交换与伞形采样——REMD 与自由能剖面版本声明块工具/软件OpenMM 8.x温度副本ReplicaExchange、自定义CustomCVForce/restraint、Python 3.9、PyMBARpymbar或 WHAM 脚本坐标 km/h、能量 kJ/mol。语言/环境Pythonopenmm、numpy、pymbar。本文参照 openmm-cookbook umbrella_sampling 教程代码思路。本文目标补齐元动力学的另两条路径——REMD 搞定高能量垒未知 CV伞形采样搞定沿已知 CV 精确算一维自由能剖面并给三方法选型表。一句话结论复制交换REMD跑一组温度由低到高的副本周期性按 Metropolis 判据尝试交换相邻副本温度让高温副本翻越高垒、低温副本保住精细构象实现无 CV、无需先验的全局采样伞形采样则先验选定一个 CV沿它切分 N 个窗口各加DISTANCE偏置谐振势中心 s0、力常数 k跑完各窗直方图再用 WHAM 或 PyMBARMultiStateSampler把 N 个窗口的重叠分布汇兑成一维自由能剖面 G(s)——REMD 适合不知道该偏置哪个坐标时兜底伞形采样适合明确了反应坐标、要精确 ΔG时实测。〇、本篇要解决的认知问题复制交换到底交换什么副本温度副本与 Metropolis 判据怎么连起来的REMD 的优势是什么——为什么它不需要预先知道 CV伞形采样为什么能算出精确的一维自由能剖面跟元动力学跑熟了再看有什么分工多窗口的偏置怎么加、直方图怎么汇兑成 G(s)WHAM 与 PyMBAR 差别在哪元动力学、REMD、伞形采样三者该在什么场景分别选用一、机制解析1.1 复制交换REMD无 CV 的全局翻垒REMD 跑一列副本replica 通常温度不同T0 T1 … T_{N−1}。相邻副本 T_i 与 T_{i1} 之间周期性尝试交换体系温度按 Metropolis 判据接受P_acc min(1, e^{−(Δβ/2)(ΔU)}), Δβ β_{i1} − β_i, ΔU U_{i1} − U_i高温副本更有能量翻越自由能垒、低温副本保留高分辨构象——通过不断换高温副本探索到的构象能下传到低温副本。普通温度 REMD 的关键约束是局部曲率一致的交换接受率通常要求温度梯度和几何尺度匹配接受率过低则交换失效。OpenMM 的ReplicaExchange把这一整套编排多个上下文、同步交换封装你在Simulation加偏置并跑replicaExchange()。1.2 分界与局限REMD 的优势是无 CV、无先验就能提升整体遍历适合盆地结构完全未知代价是交换可能因几何不匹配而失效、且比元动力学贵得多多副本并行。这也是为什么AI 增强采样里 REMD 常被拿来和 metaD 互补。1.3 伞形采样沿已知 CV 分窗偏置拉到位移伞形采样umbrella sampling沿一个先验选定的反应坐标 s按窗口中心 s0 依序切分每个窗口给体系加一个偏置势按住它典型为谐振V_umb(s) (k/2)·(s − s0)²窗口之间直方图必须重叠相邻窗口交叠是 WHAM/PyMBAR 能拼接的前提把所有窗的直方图 n_k(s) 与偏置 V_umb,k 合成无偏估计即可得 P(s)再由 F(s):−kBT·ln P(s)C 得到一维自由能剖面。1.4 WHAM 与 PyMBAR 的本质与差别两者都是多窗口统计拼接器WHAM加权直方图分析。用每窗的偏置势做迭代权重重建无偏概率把线性直达窗口数据用梯形/样条拼接简单、快但要求足够重叠且窗口间热力学充分衔接。PyMBAR多态 Bennett 接受率pymbar.MultiStateSampler把每个窗口视作一个热力学态用统计优化做整体热力学计量对重叠不完美更鲁棒、误差更小是大体系的正式版替代。一句话选 WHAM when windows overlap well code simplicity选 PyMBAR when need rigorous underestimate low overlap robustness。本文 python 用pymbar演示汇兑。二、完整代码与逐行剖析参照 openmm-cookbook umbrella 教程思路2.1 伞形采样构建体系、加沿距离 CV 的谐振偏置一个窗口importopenmmfromopenmmimportapp,unitfromopenmmimportCustomBondForceimportnumpyasnp# —— 用 OpenMM cookbook umbrella_sampling 的思路每窗口一个偏置 ——topapp.PDBFile(pocket.pdb)# 清洗好的蛋白-口袋结构forcefieldapp.ForceField(amber14-all.xml,amber14/tip3pfb.xml)systemforcefield.createSystem(top.topology,nonbondedMethodapp.PME)# 原子对配体 1 号与口袋锚定原子OpenMM 原子索引 0 起始atom_i,atom_j10,50k_fc1500.0*unit.kilojoule_per_mole/unit.nanometer**2# 偏置力常数s00.35*unit.nanometer# 当前窗口中心每个窗口一个值# CustomBondForce沿两原子距离 s 的谐振偏置 V k/2 (s - s0)^2# distance(s1,s2) 是 OpenMM 内建量perBondParameter 存每个键的 (k, s0)engCustomBondForce(0.5*k*(distance(s1,s2)-r0)^2)eng.addPerBondParameter(unit.kilojoule_per_mole/unit.nanometer**2)# keng.addPerBondParameter(unit.nanometer)# r0eng.addBond(atom_i,atom_j,[k_fc.value_in_unit(unit.kilojoule_per_mole/unit.nanometer**2),s0.value_in_unit(unit.nanometer)])system.addForce(eng)逐要点CustomBondForce的表达式里distance(s1,s2)返回 nm、沿 1 起始的 perBondParameter 索引对应k与r0addBond传入该键的(k, s0)数值。每个窗口只需把s0换成该窗口中心0.35→0.75 nm 逐窗增其余代码不变——这正是伞形采样的一次偏置、多窗口换心。2.2 全窗口循环跑 收集直方图可运行骨架逐窗口采样defrun_window(s0_nm,nsteps20000):对给定窗口中心 s0_nm 建 Simulation 跑一段返回该窗口的距离采样数组。systemforcefield.createSystem(top.topology,nonbondedMethodapp.PME)# 每窗新 systemengCustomBondForce(0.5*k*(distance(s1,s2)-r0)^2)eng.addPerBondParameter(1.0/unit.nanometer**2)eng.addPerBondParameter(1.0)eng.addBond(atom_i,atom_j,[k_fc.value_in_unit(unit.kilojoule_per_mole/unit.nanometer**2),s0_nm])system.addForce(eng)simapp.Simulation(top.topology,system,openmm.LangevinMiddleIntegrator(300*unit.kelvin,1.0/unit.picosecond,2.0*unit.femtosecond))sim.context.setPositions(top.positions)sim.minimizeEnergy()# 记录每步距离用 CustomBondForce 计算的 distance或直接用内建 distance()dist[]for_inrange(nsteps):sim.step(1)# 完整直方图提取见 openmm cookbook此处示意把窗口采样返回给 2.4 汇兑returnnp.array(dist)windowsnp.arange(0.35,0.80,0.05)# 9 个窗口中心 (nm)all_dist[run_window(w)forwinwindows]# 各窗口的距离采样列表逐要点每窗口新开一个Simulation、重复最小化采样保证偏置独立窗口中心 0.05 nm 间距需要较小的力常数才能保证相邻直方图重叠。相邻窗口直方图必须重叠——若某窗距离越界/空窗应提高窗口中心密度或调低k_fc铁律 5 的现实体现。完整、逐行最严谨的参照是 openmm cookbook 的umbrella_sampling.ipynbhttps://openmm.github.io/openmm-cookbook/dev/notebooks/tutorials/umbrella_sampling.html它把多窗口构型、偏置力常数、Egb汇兑与误差一次讲透。2.3 REMD 温度副本OpenMM 代码思路供改造# —— REMD建立 N 个温度副本的 Simulation ——Temps[300,330,360,400]*unit.kelvin sims[]forTinTemps:system_cforcefield.createSystem(top.topology,nonbondedMethodapp.PME)integratoropenmm.LangevinMiddleIntegrator(T,1.0/unit.picosecond,2.0*unit.femtosecond)simapp.Simulation(top.topology,system_c,integrator)sim.context.setPositions(top.positions)sims.append(sim)# 封装成 replica-exchangechunkSize/attemptFrequency 参数以官方文档为准# rep ReplicaExchange(sims, replica_idx, chunkSize100, attemptFrequency100)# rep.run(1000)要点各副本用不同LangevinMiddleIntegrator温度ReplicaExchange负责同步交换与温度更新交换接受判据由库内部按 Metropolis 实现细节以官方文档为准。REMD 不需定义 CV——这正是它与 metaD 的本质分工。2.4 PyMBAR 汇兑各窗口直方图 → 自由能剖面frompymbarimportMultiStateSampler,MultiStateSamplerOptionsimportnumpyasnp# samples_k: 每窗口的 (距离采样, 偏置势)供 pymbar 统计拼接samples[np.column_stack([h,v])forh,vinzip(all_dist,all_bias)]optionsMultiStateSamplerOptions()samplerMultiStateSampler(options)sampler.fit(samples)# 拟合各态自由能# 输出的 ΔF 即各窗口间自由能差拼成 FESF(s)ΔF[窗口]C详见 pymbar 文档要点PyMBAR 把每窗口视为一个热力学态用整体统计估计重叠分布的自由能差比 WHAM 对弱重叠更稳健若要一维连续剖面常先拼各窗 ΔF 再作样条方法接口以官方文档为准。三、常见报错与排查现象根因解法REMD 交换接受率极低温度梯度过大/几何变化快减小相邻副本温度差、增加温度副本数降低交换尝试频率伞形采样某窗口直方图为空偏置力常数过大、窗口跨过强垒不重叠增大窗口重叠度或略减小力常数先检验各窗口距离分布连续WHAM/PyMBAR 抛数据不重叠/欠充分相邻窗口偏置过强无重叠区增大重叠或改采更强鲁棒的 PyMBAR只能用元动力学无 CV得到的总出彩伞形剖面却算不出反应坐标选得不对不是真正的自由能隔栅用 driver/投影检验 CV 是否横穿垒必要时回到 metaD 找 CV四、动手练习练习 1REMD 判据接 2.3 建 4–8 个副本随机温度梯度跑一段。判据replicaExchange成功完成不报错且记录中交换接受率落在可接受区间如 10%以官方文档阈值为准。练习 2伞形窗口判据沿 0.35–0.75 nm 开 9 窗口跑伞形采样。判据每个窗口的距离直方图非空且相邻窗口有重叠区重叠 bin ≥1 个。练习 3剖面一致性判据用 PyMBAR 拼出 G(s)。判据自由能剖面主谷位置与第 09 篇 metaD 重建 FES 的谷位置一致在窗口分辨率内两方法交叉验证同一反应坐标。五、小结与下一篇预告本篇补齐了与元动力学互补的两条路线REMD 用温度副本 Metropolis 交换实现无先验全局翻垒OpenMMReplicaExchange伞形采样沿先验 CV 分窗加谐振偏置、再用 WHAM/PyMBAR 把多窗口直方图汇兑成一维自由能剖面。三者的选型逻辑metaD 效率高但需 CV 且可能被 CV 错选误导REMD 无需 CV 但贵且难收敛伞形采样精确但要求先验反应坐标与足够重叠。实战阶段将按先用 metaD 找 CV、再用伞形确证 ΔG的组合调用。第 11 篇预告《抗体 CDR loop 的构象集合采样》把增强采样落地到抗体 CDR H3 高热塑性环上用 metaD/gREST 采 CDR loop ensemble 并提取构象。三种增强采样方法选型对比表方法是否需要先验 CV主作用对自由能熵相对成本典型场合元动力学metaD是偏置 CV引导探索、重建 FES高效、易被 CV 错选误导低-中结合/变构、构象集合初探第 08-09 篇复制交换REMD否全局翻垒、温度巡检稳健但贵高盆地结构未知、体系小到可多副本伞形采样umbrella是反应坐标更精确精确算一维自由能剖面精确但要求重叠充分中要可靠 ΔG 的反应路径断点实测本篇认知问题回显FAQQ1复制交换交换什么Metropolis 判据怎么连A1一组温度递增的副本之间按 P_accmin(1,e^{−(Δβ/2)ΔU}) 周期性交换温度让高温副本翻越高垒、低温副本保留精细构象实现全局遍历。Q2REMD 为什么不需要预先知道 CVA2它用温度差而非坐标偏置提升遍历只要求副本温度梯度合理无需定义任何集合变量适合盆地结构未知的兜底。Q3伞形采样和元动力学的分工A3伞形采样沿先验 CV 分窗加偏置、需重叠充分给出精确一维自由能剖面元动力学不先定窗、效率高但依赖 CV 选得对。前者适合确证 ΔG后者适合初探。Q4多窗口偏置怎么加、怎么汇兑A4每窗口给体系加谐振偏置 V(k/2)(s−s0)²跑完收集直方图WHAM 用加权直方图拼接要求充分重叠PyMBAR 用多态统计优化合并、对弱重叠更鲁棒。Q5三方法各适合什么场合A5metaD 适合结合/构象初探且效率高REMD 适合未知盆地结构的全局翻垒但昂贵伞形采样适合沿清晰反应坐标精确算 ΔG。