📚 课程目录第4章自由能计算方法概览

自由能计算方法概览：选对工具，事半功倍

自由能计算的方法有很多，说实话，刚入行那会儿我也被这些缩写搞得晕头转向。FEP、TI、BAR、MBAR、MM/GBSA……每个都号称能算自由能，但到底什么时候该用哪个？这章我就把这些方法掰开揉碎，用我踩过的坑帮你理清楚。

核心思路：没有万能方法，只有最合适的工具。精度和成本永远是交易关系。

1. 自由能计算的几个基本流派

我个人习惯把自由能计算方法分成两类：一类是严格的统计力学方法，比如FEP、TI、BAR、MBAR；另一类是近似方法，比如MM/GBSA和MM/PBSA。前者精度高但计算量巨大，后者速度快但精度有天花板。

方法	全称	核心原理	精度等级	计算成本
FEP	自由能微扰	两个状态间直接微扰	★★★★★	极高
TI	热力学积分	沿路径积分能量导数	★★★★☆	高
BAR/MBAR	Bennett接受比	统计最优的比值估计	★★★★★	高
MM/GBSA	分子力学/连续溶剂	隐式溶剂模型+简正模熵	★★★☆☆	低

2. FEP：最直观，也最容易翻车

自由能微扰（FEP）的原理很直白：把两个状态A和B的配分函数比一下，然后取负对数。公式长这样：

ΔG(A→B) = -kT ln ⟨exp[-(V_B - V_A)/kT]⟩_A

嗯，这里有个坑。我曾经在一个激酶项目里直接跑FEP，从R基到R'基变换。结果算出来的ΔG和实验值差了3 kcal/mol。排查了半天，发现是λ窗口没覆盖好——两个态的能量重叠太小，采样根本不够。

避坑指南：FEP的成败高度依赖微扰前后状态的能量重叠。如果⟨exp[-(V_B-V_A)/kT]⟩严重偏离1，结果绝对不可信。我曾经犯过这个错误，白白浪费了两周算力。

实际操作时，你需要把λ从0到1切分成多个窗口。比如做9个λ窗口，每个窗口跑2 ns平衡 + 2 ns采样。脚本大致是这样：

# FEP lambda窗口配置文件 (AMBER格式)
&cntrl
  imin=0, irest=0, ntx=1,
  nstlim=500000, dt=0.002,
  ntpr=5000, ntwx=5000, ntwr=5000,
  icfe=1, clambda=0.0,
  ifsc=1,
/

clambda从0.0跑到1.0，步长0.1的话就是11个窗口。我习惯再增加两个端点窗口（0.01和0.99），防止奇点问题。

3. TI：路径积分，更稳健但更慢

热力学积分（TI）不直接算自由能差，而是对耦合参数λ求能量导数，然后沿路径积分：

ΔG = ∫₀¹ ⟨∂V/∂λ⟩_λ dλ

说白了，就是从态A到态B走一条连续路径，在每个λ点上采样平均力。这个方法的好处是——不容易出现重叠不足的问题。只要路径平滑，结果一般很稳。

但我得说，TI的计算量确实大。每个λ点都要跑独立的MD模拟。如果λ点设得太少，积分误差会很大；设得太多，Rigid-body算到天荒地老。

我的习惯做法：先跑一个5点的粗扫描（λ = 0.0, 0.25, 0.5, 0.75, 1.0），看看∂V/∂λ的曲线形状。如果曲线有剧烈拐点，再在拐点附近加密到11点或17点。这比你盲目跑20个点更有效率。

4. BAR和MBAR：统计学的力量

BAR（Bennett接受比方法）是我个人最喜欢的方法。它用最大似然估计来优化接受比，理论上在有限采样下比FEP和TI更高效。多状态版本叫MBAR，可以同时处理多个λ窗口的数据。

判断收敛的指标就是BAR的标准差。我曾经的项目里，BAR给出的误差在0.2 kcal/mol以内，而同样采样量的TI误差到了0.6 kcal/mol。差距还是很明显的。

为什么？因为BAR自动利用了正反两个方向的接受信息，而TI只是单方向积分。你想想看，双向信息肯定比单向更可靠。

用pymbar跑BAR只需要几行代码：

import pymbar

# u_kln: 每个lambda窗口每个构象的能量
mbar = pymbar.MBAR(u_kln, N_k, verbose=True)
Deltaf_ij, dDeltaf_ij = mbar.getFreeEnergyDifferences()

print(f"ΔG(A→B) = {Deltaf_ij[0, -1]:.2f} ± {dDeltaf_ij[0, -1]:.2f} kcal/mol")

注意：输入的能量矩阵必须包含所有窗口对所有构象的评估。也就是说，每个λ窗口的每一帧，你都要计算其他窗口的能量。计算量一下子就上去了。

5. MM/GBSA：快速粗筛，别当精算工具

MM/GBSA（分子力学/广义生表面溶剂）是实战中的快速筛选利器。它的思路是先对复合物、受体、配体做MD采样，然后提取轨迹帧，再用隐式溶剂模型算能量，最后加上熵校正：

ΔG_binding = ΔE_MM + ΔG_solv - TΔS

我个人习惯的用法是：先跑2ns的模拟，取最后1ns的100个构象做MM/GBSA计算。

MM/GBSA的真正战场：不是计算精确的结合自由能，而是排序——在几十个类似分子里找出活性相对较高的那几个。排序比绝对值重要得多。

但你要警惕它的缺陷。熵的贡献（TΔS）非常难算，多数人直接忽略。然而在某些体系里（比如配体构象刚柔性差异大），熵项能改变排序结果。我曾经因为省计算量没算熵，筛选出的top分子在活性实验里翻了车。

一个典型MM/GBSA输入脚本：

# MMPBSA.py输入 (AMBER)
&general
  startframe=500, endframe=1000, interval=5,
  verbose=2, keep_files=0,
/
&gb
  igb=5, saltcon=0.150, surften=0.0072,
/

igb=5是使用GB(OBC)模型，盐浓度0.15 M。我个人更喜欢igb=2（ALPB模型），在某些极性空腔体系中表现更好。两个都试试，看哪个与实验趋势一致。

6. 一张图看懂方法选型

下面这张图梳理了核心方法与应用的对应关系，帮你快速做选择。

图里有个潜规则：MM/GBSA和MM/PBSA几乎只用于结合自由能排序，别拿它去算配体构象变化或绝对溶剂化能——那是个大坑。

7. 选型决策清单

最后给你一个实战决策清单。我每次开工前都会快速过一遍：

我是做先导优化（摸清构效关系）？ → FEP 或 BAR，相对自由能。
我要判断一个新靶点能不能结合？ → 先用MM/GBSA粗筛，再用FEP验证。
我要进行500个以上化合物排序？ → MM/GBSA，因为FEP完全跑不过来。
我需要精确到0.5 kcal/mol以内的ΔΔG？ → BAR/MBAR，多次重复取统计。
我只看趋势，不需要精确值？ → MM/GBSA或分子对接打分，别浪费算力。

再次提醒：MM/GBSA的绝对自由能数值基本不可信，但它对类似化合物的相对排序有时还不错。我的经验是——先做个小型验证集（比如5-10个已知活性的分子），如果MM/GBSA的排序和实验趋势对不上，就别用了，换成FEP或BAR。

方法选对了，后续工作就顺了。这章的内容就这么多，下次你打开计算脚本时，想想自己到底想要什么——是精确的ΔG，还是快速的排序？选对了工具，你才能真正发挥自由能计算在药物优化中的价值。