AI增强构象采样教程（7）：构象采样与自由能景观的心智模型——Boltzmann 分布、集合变量与元动力学机制

AI增强构象采样教程7构象采样与自由能景观的心智模型——Boltzmann 分布、集合变量与元动力学机制版本声明块工具/软件MDAnalysis依赖于 openmm 轨迹读入、PLUMED 2.9元动力学、Python 3.9。语言/环境Pythonmdanalysis、numpy。单位距离 Å、温度 K、能量 kJ/mol势能读自轨迹。本文目标在写任何增强采样代码之前把构象概率、自由能、集合变量、偏置势四件事想通透这套心智模型决定你第 08-10 篇会不会用错 CV 和污染 FES。一句话结论构象在温度 T 下的概率服从 Boltzmann 分布 P≈e^{−U/(kBT)}自由能景观沿集合变量CV例如蛋白配体的某个二面角或 RMSD投影 F(CV)−kBT·ln P(CV)所以采样够不足的 CV 区间就等价于补全 FES 的低能谷元动力学metadynamics通过在历史访问过的区域持续累积高斯偏置势把体系从已探明的谷里顶出去这正是它比普通 MD 采样更快、也更容易制造假收敛的原因——MDAnalysis 能先让你对这条单轨迹的 RMSD/势能分布有直觉再谈偏置。〇、本篇要解决的认知问题为什么说负责测构象的频率不是概率均匀而是 Boltzmann 权重什么是集合变量CV为什么要把它从全体系坐标里降维出来构象采样任务到底难在哪——为什么普通 MD 在毫秒/微秒尺度上不够元动力学如何工作偏置势是什么、为什么它依赖历史历史依赖自由能景观FES怎么从采样得到的 CV 直方图里还原出来一、机制解析1.1 Boltzmann 分布构象不是均匀出现的系统在温度 T 处于微态 x 的概率正比于 g(x)·e^{−U(x)/(kBT)}其中 U(x) 是势能kB 是玻尔兹曼常数g(x) 是简并度。这意味着低能态占比指数级占优一个 5 kJ/mol 的能量差在 300 KkBT≈2.49 kJ/mol下对应的概率比约为 e^{5/2.49}≈7.4 倍。于是采样职能的本质是——在状态数量随自由度指数爆炸的高维空间里让马尔可夫链在可接受时间内覆盖那些概率非零的低能构象。这就是经典 MD 的瓶颈时间步长受制于最快的振动典型 1–2 fs只能走到微秒量级而蛋白配体复合物的构象重排、蛋白-配体结合/解离常出现在毫秒到秒尺度见第 01 篇你让普通 MD 跑 1 μs可能只在一个自由能盆地里打转其他低能盆地根本探访不到。1.2 集合变量CV全体系坐标的降维坐标全体系有 3N 个笛卡尔坐标但真正决定构象是否不同的往往只是少数几个慢自由度的组合。集合变量集体变量collective variableCV就是一组从体系坐标映射到低维的标量/向量函数例如CV定义典型取值/单位DISTANCE两原子/两基团距离nm 或 ÅANGLE三原子夹角弧度/度TORSION四原子二面角弧度/度RMSD相对参考结构的主链/口袋匹配ÅCOORDINATIONNUMBER对某组的配位数接触数无单位0~NSASA溶剂可及表面积Ų选择 CV 是一把双刃剑CV 必须能横穿你关心的自由能垒否则偏置白费但也绝不能把体系中本会关联的自由度全塞进一个数而失去物理辨信息。心智模型是好的 CV ≈ 反应坐标reaction coordinate——它单调穿过起始态与目标态之间的过渡态。1.3 FES从直方图还原自由能若我们沿 CVs 把构象空间分成 bin分箱则直方图 n(s)∝ 该箱内采样数当采样充分后F(s) −kBT·ln n(s) C即 FES 是采样直方图的负对数差一个常数。这条公式是全文的引擎第 09 篇plumed sum_hills、第 10 篇 WHAM/PyMBAR 最终都在回归这条式子。而普通 MD 采样不足的箱n(s)≈0ln 发散的——所以高自由能区根本没法被可靠估计这正是要引入增强采样的根本动机。1.4 元动力学历史依赖的偏置势元动力学metadynamics维护一个随历史累积的偏置势 V_bias(s,t)每当体系在某 CV 值 s 附近停留时间达到 PACE 步就叠加一个高斯数V_bias(s,t) W·exp[−((s−s(t))²)/(2·σ²)]参数与单位数值以官方文档为准HEIGHTW单次高斯高度单位 kJ/mol元动力学能量默认 kcal/mol 也可配。SIGMAσ高斯宽度单位同 CVÅ/度/弧度决定偏置沿 CV 的分辨率。PACE每多少步沉积一个高斯单位模拟步。良温元动力学well-tempered加BIASFACTOR γ高度随已累积势指数衰减确保最终收敛避免元动力学过度填充。机制已探过的区域被持续垫高的势能顶出去 → 体系被迫翻越旧谷去新盆地 → 历史累积的 −V_bias 恰是 −F(s) 的负像于是 F(s)−V_bias(∞)普通版本到极限或 F(s)∝V_bias良温到极限。这就是历史依赖两字的核心偏置来自过去的访问记录也因此收敛必须监控 bias 衰减详见铁律 5。ASCII 概念图F(CV) 叠加偏置势历史累积 高斯随时间沉积 高 │╲ ╱ ╲ ╱ ┊ ┊ ┊ ┃ │ ╲ ┌┐ ╱ → ╲╱ ···· ▓▓ (V_bias 垫高旧谷) │ ╲ ┌┘└┐ ╱ 偏向 体系被顶出→翻越右谷 低 │ └──┘ └─┘ F₀ 左谷被填 └──────────────────────── └──────────────────────→ s 两谷均深普通MD困在左谷 元动力学达到右谷得到全局 FES1.5 MDAnalysis先看轨迹再谈偏置在把元动力学涂上体系之前最该做的是用 MDAnalysis 对自己的一条 MD 轨迹做体检读轨迹、算 RMSD相对参考构象、提取势能看采样是否真的在该盆地内游荡。这既不依赖 PLUMED又能让你对偏置到底改变了什么产生直觉参照。二、完整代码与逐行剖析2.1 最小分析脚本MDAnalysis 加载单条轨迹算 RMSD 与势能importnumpyasnpimportMDAnalysisasmdafromMDAnalysis.analysisimportrms,align# 加载拓扑与轨迹top: 平衡结构 pdb/grotraj: 一条 MD 轨迹umda.Universe(npt.gro,md.xtc)refmda.Universe(npt.gro)# 参考构象可另取第一帧# 选主链原子组做对齐与 RMSD减掉整体平动/转动影响selu.select_atoms(name CA)# Cα 原子组蛋白骨架代表# 逐帧把轨迹对齐到参考再计算 RMSDrmsd_ca[]fortsinu.trajectory:align.alignto(u,ref,selectname CA)# 旋转平移对齐rrms.rmsd(u.select_atoms(name CA).positions,ref.select_atoms(name CA).positions)rmsd_ca.append(r)rmsd_canp.array(rmsd_ca)# 势能从 GROMACS 能量文件读-f md.edrMDAnalysis 需 reader# 此处演示用 numpy 直接打印 RMA 分布势能读取见下方注释print(RMSD (Å): min%.2f max%.2f mean%.2f std%.2f%(rmsd_ca.min(),rmsd_ca.max(),rmsd_ca.mean(),rmsd_ca.std()))逐行要点u mda.Universe(top, traj)第一参数拓扑含原子名第二参数轨迹支持pdb/gro xtc/dcd。alignto对齐后再rms计算否则刚体平动/转动会污染 RMSD。势能可从md.edr用MDAnalysis的能量读取模块或gmx energy -f md.edr -o energy.xvg导出再用numpy.loadtxt读列。2.2 势能直方图的 numpy 快速还原ΔG vs 计数importnumpyasnp kB0.0083144621# kJ/(mol·K)T300.0kTkB*T# ≈2.494 kJ/mol300 K# 直接把一条轨迹的势能列已从 gmx energy 的 Potential 导出 energy.xvg# 读成数组转成 -kBT·ln P 的相对自由能示意bin 计数归一成概率potnp.loadtxt(energy.xvg,comments[#,])[:,1]hist,edgesnp.histogram(pot,bins50)Phist/hist.sum()F-kT*np.log(P[P0])# 只对非零箱取对数print(自由能盆地取样点约,int(P[P0].sum()*len(pot)),帧)print(相对 F 跨度 (kJ/mol):,F.max()-F.min())三、常见报错与排查现象根因解法rms.rmsd返回空/报尺寸不符两组原子数不一致原子名选择错误确认sel与ref用同一select_atoms表达式长度相等轨迹/拓扑对不上读入原子数 mismatchtop 与 traj 分别对应不同体系统一用平衡结构做 top轨迹对应的体系必须一致直方图大量P0导致 F 发散采样不足的 bin这是普通 MD 覆盖不足的信号正是引增强采样的动机增大 bins 或轨迹长度势能范围异常巨大上亿未选取平衡段NPT 之前的过渡段污染用trajectory切片跳过前 N 帧如u.trajectory[N:]四、动手练习练习 1RMSD 判据对第 06 篇平衡后的md.xtc加载计算 RMA。判据std 1.5 Å说明采样在一个盆地内未翻越若 std 达 3 Å 以上说明体系已跨过势垒——两种不同答案对应不同的元动力学设计。练习 2FES 直觉判据用 2.2 脚本把势能转成相对自由能。判据F.max()-F.min()在 3–8 kJ/mol 区间约相当 1–3 kBT说明该盆地内的热涨落可被普通 MD 覆盖若跨越更大则需引入偏置。练习 3CV 可视化把第 1 步 Cα RMSD 存成ison数组并画直方图。判据直方图是单峰——它是单一盆地的客观佐证与练习 1 的 std 结论一致。五、小结与下一篇预告本篇建立了增强采样的心智模型Boltzmann 分布决定构象概率、CV 把高维空间降维、FES 是采样直方图的负对数、元动力学用历史累积的高斯偏置势填谷翻越、良温版本保证收敛并用 MDAnalysis 对一个实际轨迹做了 RMSD 与势能的最小体检先具备偏置改变了什么的直觉。这是后三篇首个元动力学实验、PLUMED 深度、REMD/伞形的共同坐标系。第 08 篇预告《第一个增强采样实验OpenMMPLUMED 元动力学》将用 openmm-plumed 的PlumedForce定义 1 个 CV 跑 well-tempered metaD并解读 HILLS 输出。本篇认知问题回显FAQQ1为什么构象不是均匀出现的要服从 Boltzmann 权重A1系统在热平衡下取状态 x 的概率正比于 e^{−U(x)/(kBT)}低能态以指数级占优5 kJ/mol 的能量差在 300 K 下使高能态概率约为低能态的 1/7.4故采样本质上覆盖的是少量低能盆地。Q2什么力集合变量CV为何要降维A2CV 是从 3N 维坐标映射到低维的标量函数距离、二面角、RMSD 等只保留决定构象差异的慢自由度使自由能可以被投影成可理解的低维 FES。Q3构象采样难在哪A3普通 MD 时间步受 1–2 fs 限制只能走微秒量级而蛋白配体重排常出现在毫秒到秒尺度体系被卡在单个低能盆地跨不过需穿越的自由能垒高自由能区间根本无法估计。Q4元动力学如何工作A4它随历史在已访问的 CV 区域持续叠加高斯偏置势HEIGHT/SIGMA/PACE 控制强度把已探明盆地垫高、让体系翻越到未探区域良温版本按 BIASFACTOR 衰减高度保证收敛。Q5FES 如何从直方图还原A5由 F(s)−kBT·ln P(s)C 还原即采样直方图的负对数加常数散列到 P0 的箱 F 发散这正是普通 MD 覆盖不足、必须引增强采样的原因。