自由能扰动(FEP)原理：指数平均与Bennett接受比

1. 从Zwanzig公式说起

自由能扰动（Free Energy Perturbation, FEP）的核心，说白了就是拿一个态的采样去估计另一个态的自由能差。最简单的表达式来自Zwanzig：

ΔF = -k_B T ln ⟨exp(-β ΔU)⟩_A

其中 ΔU = U_B - U_A，⟨·⟩_A 表示在态 A 的系综中求平均。这个公式看起来干净利落，但实际用起来坑不少。我最早做FEP时，以为只要跑一条A态的轨迹，把能量差指数平均一下就能交差了。结果算出来的自由能差反复震荡，换了几个随机种子也不稳定。后来才意识到：指数平均的收敛性，比想象中差得多。

2. 指数平均法（EXP）——最直接的方案

指数平均，英文叫Exponential Averaging，缩写EXP。它的思路很简单：

• 采样N个构型，来自态A的平衡分布
• 计算每个构型上的势能差 ΔU
• 计算 exp(-β ΔU) 的算术平均，再取对数

但这里有个关键问题：指数函数会把负的ΔU放大，而正的ΔU几乎被压缩到零。一旦分布尾部出现几个很大的负ΔU（即B态能量显著低于A态），它们会对平均值产生巨大影响，甚至主导结果。你想想看，几十个正常采样抵不过一个“幸运”构型。我在一个蛋白-配体结合自由能项目中就遇到过：用EXP算出的结果比实验值偏低了3 kcal/mol，因为采样跑进了一个配体偶然卡进疏水口袋的罕见态，能量差异常高。

3. Bennett接受比法（BAR）——更稳健的改进

Bennett接受比法（BAR）是EXP的升级版。它不只用A态的采样，还把B态的采样也拿进来，通过求解一个非线性方程得到最大似然的自由能差。公式长这样：

ΔF = -k_B T ln( ⟨f(-β(ΔU - C))⟩_A / ⟨f(β(ΔU - C))⟩_B ) + C

其中 f(x) = 1/(1+e^x) 是费米-狄拉克函数，C 是一个需要自洽求解的偏移量。迭代过程一般收敛很快，三五步就能稳定。

BAR相比EXP的优势很明显：

• 它把两个方向的采样都考虑进来，相当于“双向桥接”
• 方差更小，尤其在重叠不足时仍能保持偏差可控
• 对分布尾部不那么敏感

我用BAR接过好几个之前EXP跑崩的案例。有一次做自由能微扰的窗口设置（15个λ窗口），EXP算出的总自由能差曲线发“毛”，标准差0.8 kcal/mol。换成BAR后，标准差降到0.2 kcal/mol，曲线也光滑了。嗯，这里要注意：BAR在低重叠时仍然有偏差，但远小于EXP。

4. 实战对比：EXP vs BAR

下面这张表总结了我个人在项目中的对比经验：

为什么BAR对尾部的容忍度更高？你可以这么理解：EXP只问“从A看B，B有多好”，而BAR问“从两边看，两边有多像”。费米函数有天然的截断能力，不会让单个极端值把整个平均值带偏。

5. 我的建议与避坑

讲几个我亲自踩过的坑：

1. 别迷信“采样越长越准”。 EXP的方差不会随采样量平方根衰减——它衰减得慢多了。我曾经跑过100 ns的A态轨迹，用EXP算一个窗口的ΔF，结果跑10个不同种子得到1.5 kcal/mol的散布。换成BAR，同样数据量，散布只有0.3 kcal/mol。
2. 迭代求解BAR时要监控收敛。 我习惯写一个小循环，要求前后两次ΔF差异小于0.001 k_B T。如果20步还不收敛，很可能是初始猜测太差，或者两个态几乎不重叠。
3. 不要跳过双向采样，只为了省一半计算量。 很多教程说EXP只需要A态采样就够了，但代价是结果可能完全错误。我一般会跑A→B和B→A各一半的总步数，然后用BAR合并，效果远好于双倍步数只跑一个方向再用EXP。

# BAR 迭代的 Python 伪代码（实际项目中会用 pymbar）
C = 0.0
for i in range(20):
    numerator   = mean( fermi( -beta * (dUA_B - C) ) )  # A 态采样
    denominator = mean( fermi(  beta * (dUB_A - C) ) )  # B 态采样
    C_new = - (1/beta) * log(numerator / denominator) + C
    if abs(C_new - C) < 1e-6:
        break
    C = C_new
deltaF = C

上面这段代码虽然简单，但足够对付单窗口的BAR计算。实际项目中我会用 pymbar 库，它处理了自动分组、统计误差估计和多个态同时分析。但理解这段逻辑，你就能明白BAR到底在干什么。