手机版

自由能计算 — 从伞形采样到热力学积分的四种方法对比

发布时间：2026-07-16 来源：科研学术网

小中大

字号：

自由能计算是分子模拟中最高频的需求之一——离子通道的透过自由能、配体-蛋白的结合自由能、化学反应的活化自由能——这些问题的答案最终都归结为一个ΔG。但计算ΔG的方法多达四五种，每种方法的适用范围、计算成本和精度差异巨大。选错了方法，轻则多花几倍的计算时间，重则算出来的ΔG和实验差出量级上的差距。这篇文章把自由能计算的四条主流路线并排对比，给出明确的选型逻辑。

一、自由能计算的四条主流路线

伞形采样（Umbrella Sampling）是求沿反应坐标的平均力势（PMF）。把一个反应坐标（如离子通道轴向位置）切成若干窗口，每个窗口加一个谐振势约束在窗口中心附近采样，最后用WHAM（加权直方图分析法）拼接出完整的PMF曲线。适用场景：离子/小分子沿通道的透过自由能、蛋白-蛋白解离PMF、配体从蛋白解离的路径自由能。

热力学积分（Thermodynamic Integration）是通过λ参数连接初态和末态，积分dU/dλ对λ的系综平均得到ΔG。适用场景：小分子在水和有机相之间分配的溶剂化自由能、氨基酸突变对蛋白稳定性的影响。

FEP（自由能微扰）和TI在原理上类似（都是λ路径法），但FEP用相邻窗口之间的自由能差递推累加，而TI直接积分连续函数。适用场景：配体结合自由能精确排名、先导化合物优化中的R基团修饰效果预测。

MM/PBSA是事后分析法，对一条MD轨迹做后处理，分别计算气相相互作用能、极性溶剂化能（PB方程）、非极性溶剂化能（SASA）的平均值。适用场景：大量配体的快速筛选排名。

二、伞形采样——窗口密度和力常数的选择

伞形采样是自由能计算中最”手工”的方法——窗口数量、每个窗口的位置和力常数都需要自己设定，设定的质量直接影响PMF的精度。

窗口密度：对于沿z轴的离子通道PMF计算（通道长度约40 Å），通常设置40-60个窗口，间距0.5-1.0 Å。间距不能太大——相邻窗口的采样分布必须有足够重叠，让WHAM能可靠拼接。重叠度的最低要求是：相邻两个窗口的质心分布直方图的半高全宽区域有交集，重叠区域的采样量≥10%的总采样量。

力常数：谐振约束的力常数k决定了每个窗口的采样分布宽度。k太小，离子容易漂出窗口，采样不充分；k太大，分布太窄，需要的窗口数量暴增。对单价离子（Na⁺、K⁺），k取2-5 kcal/(mol·Å²)为宜，对应的分布半宽约0.5-0.8 Å。

力常数和窗口间距的耦合关系：k越大→分布越窄→窗口间距应该越小。如果你用k=5、间距取1.0 Å，相邻窗口重叠度不够，PMF曲线会出现”锯齿”——这不是物理信号，是WHAM拼接失败的人为噪声。

三、热力学积分——λ路径上的收敛判据

热力学积分的核心操作是在λ=0到λ=1之间取若干采样点，每个λ值跑一段MD，计算dU/dλ的系综平均，然后数值积分。λ取值点的分布不是等间距最优——在dU/dλ变化剧烈的区域（通常靠近λ=0和λ=1）需要加密。

TI收敛的判据有两条。第一：正向（λ从0到1）和反向（λ从1到0）算出来的ΔG应该一致，偏差<1 kcal/mol。如果正向-反向偏差>2 kcal/mol，说明λ采样点不够或每个点的模拟时间不足。

第二：每个λ点的dU/dλ时间序列应该是白噪声——如果存在明显的长期漂移趋势，说明该λ点的体系没有达到平衡。对于涉及大构象变化的突变（如色氨酸→丙氨酸，蛋白核心出现空穴），λ靠近1的点需要更长的平衡时间（5-10 ns vs常规的1-2 ns）。

四、FEP——小分子结合自由能计算的精度极限

FEP（自由能微扰）理论上的精度是四种方法中最高的——因为它对构象空间的采样是通过MD在每个λ窗口充分平衡和采样的，而不是事后分析。但FEP的”精度极限”取决于两个因素：力场的精度和采样的充分性。

力场精度是硬天花板。用GAFF通用小分子力场做FEP，即使采样完全充分，结合自由能和实验的RMSD也在1.5-2.0 kcal/mol——这不是FEP方法的问题，是力场本身的误差。用专门针对类药物分子优化的力场（如OpenFF或拟合了QM数据的自定义参数），RMSD可以压缩到1.0 kcal/mol以内。

采样的充分性则取决于体系。如果配体结合伴随着蛋白的构象变化（如loop从开变关），而MD的时间尺度不足以覆盖这种变化，FEP的结果就是偏的——它只采样了单一的蛋白构象子态，没有捕捉到结合过程中的构象选择效应。

五、MM/PBSA——从轨迹中提取自由能的效率权衡

MM/PBSA的位置很微妙——精度不如FEP和TI，速度远快于伞形采样。它最适合的场景是：你需要对30-50个候选配体做快速排名，不在乎绝对值的精度，只需要知道A比B大概率更强。

MM/PBSA的核心效率来源于它不额外增加MD模拟——它对已有的MD轨迹做后处理。一条20 ns的蛋白-配体复合物轨迹，取最后10 ns做200帧均匀采样，PB计算每帧约30秒，总后处理时间约2小时——比FEP快两个数量级。

但MM/PBSA有一个结构性的缺陷：熵项的估算。结合自由能中的熵变（-TΔS）在MM/PBSA中通常用简正模分析估算，这个过程的不确定度约3-5 kcal/mol。很多MM/PBSA研究干脆忽略熵项——不做-TΔS的MM/PBSA叫MM/PBSA（少了ΔS），实际上变成了”焓主导的结合能排名”，而非真正的结合自由能。