📚 课程目录第1章 FEP自由能计算简介与原理

第1章：FEP自由能计算简介与原理

我知道很多刚接触自由能计算的同学，第一反应就是「这玩意到底准不准？值得花那么大功夫搭建pipeline吗？」。我在2018年第一次做FEP的时候，也带着这个疑问。后来跑完几个项目才摸清楚——FEP不是玩具，它是药物设计里少数能给出定量结合自由能差值的方法。

1.1 什么是自由能，为什么我们需要它？

自由能，说白了就是系统在特定温度压力下能做「有用功」的那部分能量。在分子模拟里，我们最关心的是结合自由能（ΔG_bind）和相对结合自由能（ΔΔG）。

你想想看，一个药物分子和靶点结合得牢不牢，最终判断标准就是ΔG。实验上测的是IC₅₀或K_d，但计算能直接给出ΔG，省去合成和测试的时间。我做过的案例里，FEP预测的ΔΔG和实验值偏差通常在1 kcal/mol以内——这个精度足以指导先导化合物优化。

  关键概念： ΔG = –RT ln(Keq)。1 kcal/mol的自由能差异，对应结合常数差异约5倍。所以哪怕只算准0.5 kcal/mol，也能淘汰掉一半以上的候选分子。

1.2 FEP方法的核心原理

FEP的全称是Free Energy Perturbation，中文叫「自由能微扰」。它的数学基础其实很简单——用统计力学里的Zwanzig公式：

ΔG(A→B) = –k_BT ln⟨exp[–(U_B – U_A)/k_BT]⟩_A

公式里⟨⟩_A表示在状态A的系综下取平均。说白了，我们只需要在状态A的模拟轨迹里，分别算出状态A和状态B的能量差，然后取指数平均，就能得到A到B的自由能变化。

但这个公式有个致命问题——只有当A和B的相空间重叠足够大时，计算结果才可靠。我最初做FEP时，一上来就试了苯环上取代基从‑Cl直接变‑CF₃，结果收敛性差得要命。后来才学会用多步微扰（也叫λ窗口）。

1.3 耦合参数λ——核心工具

为了让相空间重叠足够好，我们把A到B的转变拆成若干中间态，每个中间态由一个耦合参数λ控制：

U(λ) = (1 – λ) U_A + λ U_B

λ从0到1平滑变化。通常我们设置11～21个λ窗口，每个窗口跑几千步的平衡+采样。这种方式就叫热力学积分（TI）或逐步微扰（STEP）。

我记得有一次，有个学生问我：「为什么不用单步FEP直接算？省时间啊。」我告诉他：单步FEP只适用于非常相似的分子（比如‑H变‑CH₃）。对大部分改造，窗口数少于10个就会产生不可忽视的滞后误差。这不是理论问题，是实践里的坑——我踩过。

ΔG计算方法	特点	典型窗口数	适用场景
单步FEP (Zwanzig)	直接指数平均，快但依赖相空间重叠	1	极相似结构（‑H/‑CH₃）
逐步FEP (BAR/MBAR)	多窗口加 Bennett 接受率，精度高	11～21	中等差异（‑Cl/‑CF₃、‑OH/‑OCH₃）
热力学积分 (TI)	对λ梯度积分，适合连续变化	15～31	带电变化、成键/断键

个人经验： 做药物改造时，我通常先跑7个窗口（λ = 0, 0.17, 0.33, 0.50, 0.67, 0.83, 1.0）做粗扫，看到ΔΔG超过1.5 kcal/mol再加密到15窗口。这样能省50%的计算量。

1.4 从原理到Pipeline——一张图说清楚

下面这张图概括了FEP计算的核心逻辑链：

_A 和 U_B 到文件 BAR / MBAR 多权重→ΔG(A→B) 误差分析滞回、统计误差 ⚠ 关键注意事项 • λ窗口重叠: 相邻窗口的相空间重叠系数应 ≥ 0.3 • 平衡时间: 每个窗口至少500 ps 删除以消除初始偏置 • 静电力处理: 使用PME + 长程校正，否则ΔG偏差可达2 kcal/mol

我的习惯： 每个λ窗口的模拟长度我一般不固定，而是跑一个自适应收敛检测——每100 ps计算一次ΔG的滑动平均，直到连续3次标准差小于0.1 kcal/mol才停止。这比固定跑2 ns更高效。

1.5 常见的「虚功」误区

做FEP最容易犯的错误是把「自由能」和「内能」混为一谈。我在项目评审会上见过好几次，有人提交的ΔG结果其实只算了焓变（ΔH）。自由能里包含熵项（–TΔS），而熵的变化在结合过程中占比很大——尤其是配体从溶液进入蛋白口袋时，水分子的释放贡献了可观的熵增。

另一个我踩过的坑是静电力截断。早期我用8 Å的截断半径做FEP，结果算出来的ΔΔG全是正的——预测全部失败。后来改成PME（粒子网格Ewald）加12 Å实空间截断，RMSD才降到0.8 kcal/mol。

避坑提醒： 如果你用AMBER或GROMACS跑FEP，务必确认静电算法是「PME」而不是「cut-off」。前者在GPU上性能不错，别舍不得这点算力。

1.6 什么时候应该用FEP？

不是所有问题都适合FEP。我个人判断标准有3条：

相似性：待比较的两个分子骨架重叠度高（RMSD < 0.5 Å），否则相空间重叠不够
力场覆盖：所用力场（GAFF、CGenFF、OPLS）能覆盖所有原子类型，缺失参数需要手动拟合
计算资源：单个窗口至少能跑1 ns以上的平衡+采样，GPU集群是标配

举个例子，我之前做过的EGFR抑制剂优化，从母核噻唑啉换成吡啶并酮，核心骨架RMSD只有0.3 Å，11窗口的FEP就给出了ΔΔG=–1.2 kcal/mol，后续合成验证的IC₅₀确实提升了3倍。反过来，如果两个配体结合构象完全不同，FEP就爱莫能助了——这种情况我建议用MM/GBSA或Metadynamics。

1.7 为什么选FEP作为Pipeline核心？

市面上有很多自由能方法：TI、FEP、BAR、MBAR、WHAM、MM/GBSA、MM/PBSA……但FEP（结合BAR或MBAR后处理）在相对结合自由能计算上精度最高、收敛最快。原因是FEP的微扰范围小，每个窗口的构象变化不大，所以相空间采样充分。

我自己的pipeline里，后处理统一用MBAR替代原来的BAR——因为MBAR能一次性利用所有窗口的能量数据，得到全局最优的无偏估计。虽然计算量稍大，但误差条带更窄，而且能直接给出统计不确定性。