6.1 自由能计算 TI 与 FEP

6.1 自由能计算 TI 与 FEP

本节摘要：自由能差是化学与材料问题的通用货币（溶解、分配、结合、相变），但普通模拟给不出它——因为直接的"硬切换"（把一个分子瞬间变成另一个）概率为零。解法是炼金术：把变换拆成 λ 从 0 到 1 的一串中间态，每个态跑一次普通 MD，把 dH/dλ 积分起来。本节以甲苯的溶剂化自由能为演算主线，讲清热力学循环、软偶合、λ 方案设计与梯度积分的完整流程。

5.3 的统计量都在"同一个物理体系内部"打转，本节开始问跨态问题：这个分子在水里比在辛醇里舒服多少（分配系数的来源）？这个官能团换成另一个，结合会强多少（药物改造的核心问题）？这类量的定义就是两态自由能差 ΔG = −k_B T ln Z_B/Z_A，而两态之间如果直接切换，重叠概率小到宇宙年龄都采不到一个样本——这就是普通 MD 够不着自由能的原因。

热力学循环：差值可以绕道算

直接算不了的差值，可以借一个假想的中间态绕道。算相对结合自由能的经典循环：配体 A 与 B 竞争同一个结合位点，我们要的 ΔΔG = ΔG_bind^B − ΔG_bind^A（真实过程：两个解结合实验），可以用"在结合态把 A 变成 B、在溶液态把 A 变成 B"这两条非物理但可模拟的路径拼出来：

ΔΔG = ΔG_alch(bound) − ΔG_alch(solvated)

绕道的合法性由状态函数性质保证：自由能只看首末态、不看路径，物理路径与非物理路径的首末态相同，差值必然相等。这条"绕道自由"是全部炼金术方法的理论地基，而工程问题只剩一个：非物理路径怎么拆成可采样的台阶。

答案是 λ 耦合参数：把 B 的相互作用（电荷与 LJ）按 (1−λ)·A + λ·B 插值，λ 从 0（纯 A）走到 1（纯 B）。拆成 M 个 λ 点（典型 11 到 21 个），每点一次平衡加采样。拆分的精度陷阱在端点：λ=0 附近 B 的粒子几乎无相互作用，可能与 A 的原子重叠——软偶合（soft-core）势在 λ 中段修改 LJ 的距离依赖，把 r⁻¹² 的无穷排斥压成有限值，专门治这个"端点重叠"病。GROMACS 里 sc-alpha、sc-sigma 两个参数即为此设，默认值（0.5 与 0.3）通常不必动。

图 6-1 溶剂化自由能的热力学循环与 λ 阶梯

mdp 方案与会话

甲苯水溶剂化是自由能方法的"果蝇实验"（实验值约 −3.3 kcal/mol，即约 −13.8 kJ/mol，适合验流程）。关键 mdp 段：

free-energy              = yes
init-lambda              = 0          ; 每个窗口改它，0, 0.1, ..., 1.0
delta-lambda             = 0          ; 慢增长模式才非零；窗口法保持 0
couple-moltype           = TOL        ; 被耦合的分子名
couple-lambda0           = vdw-q      ; lambda=0 时范德华与电荷全开
couple-lambda1           = none       ; lambda=1 时全部关掉（消失）
sc-alpha                 = 0.5
couple-intramol          = no         ; 分子内项不随 lambda 变，省一半麻烦
nstdhdl                  = 10         ; 每 20 fs 记一次 dH/dlambda
calc-lambda-neighbors    = -1         ; 记录全部相邻态的能量差，事后 BAR 可用
​

每个 λ 窗口的完整流程与第3、4章完全一致（EM、NVT、NPT、生产各 1 到 5 ns）——λ 窗口数 × 上述时长就是预算，这就是 6.3 存在的理由。采样完导出梯度文件（每窗口一个 dhdl.xvg），分析交给 gmx bar 或本节的积分演算。

演算：把 11 个梯度积成自由能

TI（热力学积分）的公式：ΔG = ∫₀¹ ⟨∂H/∂λ⟩ dλ。数值上就是梯形法。下面脚本对一套甲苯溶剂化的典型梯度数据积分，并演示"窗口少了会怎样"：

import numpy as np

# 11 个 lambda 窗口的 <dH/dlambda>（kJ/mol），甲苯水溶剂化的典型形态：
# 电荷项贡献集中在 lambda 前段（水对偶极响应快），色散腔项拖到后段
lam = np.linspace(0, 1, 11)
dhdl = np.array([-9.2, -11.8, -14.1, -15.6, -15.2, -13.1,
                 -9.8, -6.4, -4.1, -2.6, -1.9])

dg_trapz = np.trapz(dhdl, lam)
print(f"11 窗 TI 梯形积分: ΔG = {dg_trapz:.2f} kJ/mol（实验参考约 -13.8）")

# 稀疏化实验：只用 6 个窗口（步长 0.2）再积一遍
lam6, dhdl6 = lam[::2], dhdl[::2]
print(f" 6 窗 TI 梯形积分: ΔG = {np.trapz(dhdl6, lam6):.2f} kJ/mol")
# 曲线弯曲处丢点，误差立刻显现 -> 窗口密度要沿曲率分布，不是均匀分布

# 不确定度的朴素估计：每窗口梯度标准差 1 kJ/mol、独立样本 500 个
sigma_w = 1.0 / np.sqrt(500)
print(f"单窗口梯度标准误约 {sigma_w:.3f} kJ/mol，11 窗独立误差合成 "
      f"约 {sigma_w * np.sqrt(11) * 0.5:.2f} kJ/mol")
​

两点判读：稀疏化实验是窗口设计的自检手段——把 11 点抽成 6 点重积，差值超过目标精度就说明曲率大的区段窗口不够；误差合成给出"要 1 kJ/mol 精度需要多长采样"的量化依据。BAR/FEP 走的是另一条数学路（两态重叠的双侧指数平均，gmx bar 一条命令），它利用了相邻窗口的双向信息，同剂量采样下精度普遍优于 TI——实践建议：跑窗口法采样，BAR 与 TI 都出数，两者一致是流程健康的最强证据，差得多则先怀疑窗口重叠不足。

四种估计方法的家族对照

TI 与 FEP 之外还有两个名字高频出现，四种方法一张表分清：

实操推荐路径：grompp 里 calc-lambda-neighbors = -1 让每个窗口记录到所有相邻态的能量差，采样结束后 BAR 直接出全表；想再榨精度，把同一批数据丢给 MBAR 类后处理脚本重算，零额外模拟成本。TI 保留为对照——BAR 与 TI 相差超过各自误差棒时，第一嫌疑是窗口重叠不足，第二嫌疑是某窗口没跑平衡。

演算：从能量差直估一个窗口对的 BAR

BAR 的名字唬人，单窗口对的含义朴素：两边重叠区的能量差分布互为镜像，自洽解出 ΔF。用两个正态样本演示数量级（真实分析交给 gmx bar，这里建立直觉）：

import numpy as np

rng = np.random.default_rng(9)
n = 5000
# 两个相邻 lambda 窗口：能量差近似正态，均值相反、方差相同（对称近邻）
dA = rng.normal(-6.0, 2.0, n)    # 从窗口 A 看去"变到 B"的能量差
dB = rng.normal(6.0, 2.0, n)     # 从窗口 B 看去"变回 A"的能量差

# BAR 方程的朴素定点迭代：f = ln<exp(-beta(dA-f))> - ln<exp(-beta(dB+f))>
KB = 0.0083144621; T = 300; beta = 1 / (KB * T)
f = 0.0
for _ in range(200):
    lhs = -np.log(np.mean(np.exp(-beta * (dA - f)))) / beta
    rhs =  np.log(np.mean(np.exp(-beta * (dB + f)))) / beta
    f = f + 0.5 * (lhs - rhs)
print(f"BAR 估计 ΔF ≈ {f:.2f} kJ/mol（设计真值 6）")

# FEP 单侧估计对照：只用 A 侧指数平均（重叠差时方差爆炸）
f_fep = -np.log(np.mean(np.exp(-beta * dA))) / beta
print(f"单侧 FEP 估计 ≈ {f_fep:.2f} kJ/mol")
​

跑几次换种子，你会看到 BAR 稳在真值附近而单侧 FEP 时有震荡——这就是"双侧信息更稳"的数值证据，也是 BAR 成为默认的原因。

两个高频疑问

λ 步长怎么定？ 沿曲线曲率定：电荷消去段（λ 前半）梯度变化陡，窗口密一档；色散腔段（后半）平缓，可以疏一档。GROMACS 支持非均匀格点（couple-lambdas 或 fep-lambdas 列表），比均匀 11 点更省。判断够不够的硬办法就是 6.1 正文的稀疏化重积分——数据在手，答案自见。

相对自由能与绝对自由能哪个先做？ 会算绝对溶剂化（本节的甲苯）是一切的地基：流程短、有实验值可对、能验证你的 λ 方案与 BAR 管线全通。直接上相对结合自由能（多窗口加多配体）是常见翻车姿势——管线里的坑在简单体系暴露是学费，在课题体系暴露是事故。

速记：直接切换采不到，热力学循环把差值拆成两条非物理路径；λ 阶梯加软偶合防端点重叠；couple-intramol = no 省掉分子内校正；TI 是梯形积分、BAR 用双侧重叠信息，两算对照是健康检查；窗口密度沿曲率铺，预算 = 窗口数 × 单窗时长。整条循环的思路有了，下一节换一条路：沿着反应坐标一段一段攻。