GROMACS自由能计算方法对比：FEP自由能微扰与TI热力学积分的实战选择

“这两个分子结合得有多紧？”任何做药物设计或者酶工程的项目，最终都要回答这个问题。实验上可以用ITC（等温滴定量热法）或者SPR（表面等离子体共振）来测，但如果你需要在合成化合物之前知道候选分子的结合能力，就需要计算自由能。

在GROMACS中做自由能计算，两条经典路线是自由能微扰（FEP）和热力学积分（TI）。两套方法在马尔可夫状态采样的框架下数学上等价，具体实施中的差异决定了各自适用场景。

FEP和TI的物理直觉

FEP的逻辑很直接：把体系从状态A（比如蛋白+溶剂）渐变到状态B（蛋白+配体+溶剂），每一步只引入微小扰动，相邻两个状态的自由能差用Zwanzig方程计算。关键在于扰动量必须足够小——如果一步变化太大，状态A的采样覆盖不到状态B的构象空间，方程的分母接近于零，结果的不确定性爆炸。

TI的逻辑是积分：沿着λ路径每隔一段取一个样点，在每个样点上计算哈密顿量对λ的导数（∂H/∂λ），然后沿λ方向积分得到总自由能变化。导数本身是系综平均，采样的质量直接决定积分的精度。如果∂H/∂λ在某个λ区间剧烈变化（典型的是配体从缩小的软核态突然插入蛋白结合口袋），需要在这个区间加密采样点。

两种方法在GROMACS中的实现都比较成熟，选择更多取决于工作流习惯而非精度差异——FEP需要的λ窗口通常更多（20-40个），但每个窗口的模拟时间可以较短；TI对窗口数量要求较低（10-20个），但每个窗口需要更长的模拟以确保∂H/∂λ的统计收敛。

GROMACS中FEP的实操流程

以计算一个小分子配体在水中的溶剂化自由能为例，这是一个典型的单拓扑FEP（配体从完全耦合到完全解耦）：

第一步：准备两个状态的文件。 状态A（配体在水中，完全耦合），状态B（配体消失，只剩水）。用gmx grompp生成两个tpr文件，分别对应λ=0和λ=1的哈密顿量。

第二步：设置λ窗口。 FEP的精度对窗口分布敏感。在耦合/解耦的两端（λ接近0或1），自由能对λ的导数变化剧烈（软核势的非线性区域），窗口需要加密。一个经过验证的分布方案：0.0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0——在λ→1附近拉开到0.05步长。

第三步：软核势参数。 这是决定FEP成败的关键技术细节。当配体在λ→0或λ→1的过程中体积趋近于零时，Lennard-Jones势的非键相互作用会出现奇点（1/r¹²项在r→0时发散）。软核势（soft-core potential）用修正的LJ形式替换了近程排斥项：

sc-alpha = 0.5     # 软核参数，控制"软"的程度
sc-power = 1       # 默认值
sc-sigma = 0.3     # nm，软核半径

sc-alpha设太小（0.1）可能保留过多的硬核排斥，导致λ中间区域采样困难；设太大（1.0）则过度软化，有效体积被高估。0.5是文献中反复验证的对有机小分子配体最稳健的值。

第四步：每个λ窗口的模拟。 每个窗口先做100 ps NPT平衡（配体重原子位置约束），再做1-5 ns生产相MD。如果体系较大（>5万原子）或配体柔性较高，生产相建议延长到10 ns。收敛检查的标准：相邻两个λ窗口之间的∂H/∂λ重叠度应该足够（重叠直方图在采样范围内连续），否则说明窗口间距过大。

第五步：用BAR或MBAR分析。 GROMACS自带的gmx bar工具用Bennett接受比法（BAR）处理FEP数据。BAR的数学比直接的Zwanzig方程更稳健——它利用前向和后向采样的对称性，消除了Zwanzig方程的偏性。如果有多个λ窗口且有重叠的∂H分布，MBAR（多项Bennett接受比法）可以联合所有窗口的信息给出更精确的估计。

GROMACS中TI的实操流程

TI的mdp文件配置比FEP更简洁——不需要处理状态A和B两个tpr文件的耦合，只需在单个mdp中定义λ向量和导数计算方式：

free-energy = yes
couple-lambda0 = vdw-q      # 状态A的耦合类型
couple-lambda1 = none       # 状态B（完全解耦）
couple-moltype = LIG        # 对哪个分子做TI
nstdhdl = 10                # 每10步输出一次∂H/∂λ

; λ窗口定义
init-lambda-state = 0
calc-lambda-neighbors = 1   # 计算相邻窗口的导数，用于误差估计

TI的λ窗口分布和FEP类似——两端加密。但TI的优势在于单窗口的数据独立：如果某个λ窗口的∂H/∂λ方差过大，可以单独延长这个窗口而不影响其他窗口。这个特性在实际操作中很实用——跑完一遍TI后检查误差曲线，对有异常的窗口追加模拟，比FEP的全局重跑灵活。

gmx analyze可以对每个窗口输出的∂H/∂λ做统计分析，gmx bar中的TI模式进行数值积分并估算误差。

FEP vs TI：什么时候选哪个

对于常规的蛋白-配体结合自由能计算，FEP+BAR是目前的主流选择——窗口数虽多但每个窗口可以较短，总时间可控，BAR的误差估计也足够稳健。对于丙氨酸扫描（逐个突变界面残基）和点突变自由能变化，TI更直接——每个突变只需对比WT和突变体的∂H/∂λ积分差值。

收敛判据与验证

自由能计算的陷阱是”看起来收敛了但偏差很大”。三条硬性检查：

正反路径一致性。 前向（耦合→解耦）和后向（解耦→耦合）的自由能值之差不应超过1 kcal/mol（或2 kT）。差距大于这个值说明采样不充分，两个方向没有遍历相同的构象空间。

窗口间重叠度检查。 相邻λ窗口的∂H分布必须有充分的重叠——直观标准是两个分布的主体（均值的±2σ范围）至少重叠50%以上。

独立重复。 从不同的初始构象开始（或不同的随机种子），独立运行两遍计算。两次结果的差异如果超过1.5 kcal/mol，需要增加采样。

自由能计算在药物设计中的应用越来越广泛，但参数细节对结果的影响依然很大。科研学术网（https://www.keyanxueshu.com）-配体结合和突变自由能等多个场景。