自由能模拟计算

我第一次做 自由能模拟计算是给一个酶做单点突变的 ΔΔG 预测，只设了 6 个 λ 窗口，算出来的数值误差棒比效应还大，根本没法下结论。后来把窗口加到 16 个、每窗跑 5 ns，误差棒收窄到 0.3 kcal/mol 以内，结论才立得住。那次让我懂，自由能模拟计算的精度不是软件给的，是采样给的。

一、为什么 自由能模拟计算先做方法选型

自由能差有两条主流路线：炼金法（FEP/TI）算的是两态之间的差，精度高但要对每对状态做 λ 变换；终点法（MM/PBSA）只算初态和终态各自的能量再相减，便宜但丢掉了炼金路径的cancel误差。我做配体优化默认走 FEP，做大体系粗筛才用 MM/PBSA。曾经为了省算力全用 MM/PBSA 报绝对结合自由能，误差 2~3 kcal/mol，和实验对不上。做 自由能模拟计算时，方法选型决定了你结论的置信上限。

二、核心原理：炼金变换怎么收自由能

FEP 把两态用 λ 连续连接，每个 λ 窗口算 ΔU 的 Zwanzig 指数平均；TI 则直接对 ⟨∂U/∂λ⟩ 积分。BAR 方法在两相邻窗口间用 Bennett 接受比，比 Zwanzig 更高效。我用 gmx bar 汇总时，会先看每个窗口的自洽误差，确认相邻窗口的 ΔG 重叠充分。关键认知：自由能差可以算得很准（误差<1 kcal/mol），但绝对自由能因为熵项难采永远算不准——所以永远报差，不报绝对值。

三、关键技术要点：窗口重叠与平衡

自由能模拟计算最致命的是 λ 窗口之间没重叠，forward 和 reverse 的 ΔG 对不上。我一般要求相邻窗口的 dA/dλ 分布重叠 > 70%，否则就加密窗口或加采样。另一个坑是平衡时间不够，前 1 ns 还在弛豫就被算进平均，会系统性偏。我默认每窗口至少 2 ns 平衡 + 3 ns 生产，软核势参数（soft-core）也要调，避免 λ 接近 0/1 时能量发散。

四、实操流程：从拓扑到 ΔG

我的标准动作：构建双拓扑或单拓扑 → 规划 λ 窗口（突变区软核）→ 每窗平衡+生产 → gmx bar 汇总 → 做误差棒和重叠检查 → 和实验或文献比对。最容易被砍的是窗口重叠检查，但 自由能模拟计算 的可信度全靠它，我默认保留。

补一个水合自由能的实算：苯分子水合自由能，8 窗口 FEP 算 −1.1 kcal/mol，和实验 −0.9 接近；只设 4 窗口、采样 2 ns，结果飘到 −1.8，误差翻倍。如果拿这个去指导亲水性改造，方向就反了。这之后我做 FEP，窗口数和每窗时长都是先按重叠质量倒推，不再拍脑袋定。

五、常见踩坑：溶剂模型和力场

自由能模拟计算常见错是 TIP3P 和 OPC 水混用、力场版本不统一，导致终点法能量不可比。另一个坑是忽视长程静电（PME 截断），自由能对静电极其敏感，截断差一点 ΔG 就飘 1 kcal/mol。还有人用 MM/PBSA 时把熵（NMODE 或准谐）省略，结果在柔性大的配体上误差被放大到 3 kcal/mol 以上。

六、复盘：自由能是相对量

回过头看，自由能模拟计算最值钱的是用可控的采样把”哪个突变/哪个配体更好”量化到误差棒以内，而不是给你一个绝对数。我接项目默认先问要算差还是算绝对、要什么精度，再定窗口和时长。把方法、重叠、力场三件事串好，ΔG 才敢拿去指导分子设计。被证明有用的，是愿意为重叠质量多开四个窗口的耐心。下次有人甩给我一个”结合能 −10 kcal/mol”，我会先问——你算的是差还是绝对值，误差棒多大。