计算化学 · 自由能微扰（FEP）计算实例

AI 计算化学实例：Scientify 用自由能微扰（FEP）计算 6 个分子的水合自由能

Scientify 自行搭建分子体系、配置 GROMACS 自由能微扰参数，在租用的 GPU 上运行 120 个 λ 状态的分子动力学，再用 MBAR 求出水合自由能。

每个点是一个分子。虚线为计算值等于实验值，灰带为 ±1 kcal/mol。

与实验值的 RMSE

0.83 kcal/mol

计算量

120 个 λ 状态 · 88 ns

算力

1 × NVIDIA T4 · 约 11 小时

模拟体系

溶质、水分子和周期性盒子

下图是 Scientify 用 PyMOL 在云电脑上渲染的体系构象，坐标取自 λ = 0 状态生产模拟的最后一帧。溶质居中，每个水分子取离溶质最近的周期性镜像。

苯分子位于 547 个 TIP3P 水分子构成的菱形十二面体盒子中心的分子渲染图 — 苯与 547 个 TIP3P 水分子，共 1,653 个原子。红色为水的氧原子，白色为氢原子，深灰色为碳原子。
盒子是菱形十二面体，三维周期性边界。图中水分子已展开到溶质周围，因此呈球形分布。
正交投影，CPU 光线追踪并开启环境光遮蔽。

甲烷第一水合层 35 个水分子，无溶质与水的氢键
苯第一水合层 50 个水分子，无溶质与水的氢键
甲苯第一水合层 52 个水分子，无溶质与水的氢键
乙醇第一水合层 43 个水分子，溶质与水形成 4 条氢键（灰色虚线）
苯酚第一水合层 48 个水分子，溶质与水形成 2 条氢键（灰色虚线）
乙酰胺第一水合层 35 个水分子，溶质与水形成 3 条氢键（灰色虚线）

模型与参数

从力场到采样，每一项参数由 Scientify 配置

下列数值取自实际运行的 GROMACS 输入文件。每组参数附有原始文件片段，苯在 λ = 10 状态的完整生产参数文件可在页尾下载。

力场: GAFF
原子电荷: AM1-BCC（FreeSolv 官方参数文件）
水模型: TIP3P

原始文件mobley_3053621/system.top

[ atomtypes ]
;name   bond_type     mass     charge   ptype   sigma         epsilon       Amb
 ca       ca          0.00000  0.00000   A     3.39967e-01   3.59824e-01 ; 1.91  0.0860
 ha       ha          0.00000  0.00000   A     2.59964e-01   6.27600e-02 ; 1.46  0.0150
OW OW 15.9994 0.0 A 0.315061 0.636386
HW HW 1.008 0.0 A 0.0 0.0
 
[ atoms ]
;   nr  type  resi  res  atom  cgnr     charge      mass
     1   ca     1   MOL    C1    1     -0.13010     12.01000 ; qtot -0.130
     7   ha     1   MOL    H7    7      0.13010      1.00800 ; qtot -0.650

计算流程

建模、平衡、采样、分析，四个阶段依次完成

每个分子的 20 个 λ 状态各自走完整的平衡与采样流程，结果经检查后再进入 MBAR。

1120 个状态
建模
读取 6 个分子的 GAFF/AM1-BCC 拓扑，用 TIP3P 水溶剂化，生成 120 个 λ 状态的输入文件，并通过 grompp 检查。
218 ns 平衡
平衡
每个状态依次做能量最小化、50 ps NVT 和两段 50 ps NPT，使温度、压力和密度稳定。
370 ns 生产
采样
在 T4 上运行生产模拟，短程非键相互作用在 GPU 上计算，PME 在 CPU 上计算，每 0.2 ps 记录全部 20 个状态的能量差。
46 个分子
分析与检查
去相关后做 20 态 MBAR，检查相邻状态重叠和随采样长度的稳定性。乙酰胺的估计值仍在移动，Scientify 将其采样延长到每状态 1 ns。

λ 路径

120 个状态的能量导数

曲线是每个 λ 状态生产采样的 ⟨∂U/∂λ⟩ 平均值。关闭静电阶段，曲线起点按分子极性排列：甲烷接近 0，乙酰胺约 24 kcal/mol。关闭范德华阶段，曲线在 λ≈0.75–0.8 处出现软核路径的典型低谷，到 λ = 1 回到 0。相邻状态的 MBAR 重叠最低为 0.07。

数据来自 windows.csv。沿解耦方向计算，单位为 kcal/mol 每单位 λ。

计算结果

6 个分子的水合自由能

计算值与实验值的均方根误差为 0.83 kcal/mol，在常用的化学精度（1 kcal/mol）以内。苯和甲苯与实验值相差 0.02 kcal/mol。

分子	计算值	实验值	计算 − 实验
甲烷	2.42 ± 0.04	2.00	+0.42
苯	-0.88 ± 0.07	-0.90	+0.02
甲苯	-0.92 ± 0.08	-0.90	-0.02
乙醇	-3.55 ± 0.07	-5.00	+1.45
苯酚	-5.56 ± 0.10	-6.60	+1.04
乙酰胺	-8.84 ± 0.05	-9.71	+0.87

RMSE

0.83

MUE

0.64

R²

0.96

单位为 kcal/mol，计算值后为 MBAR 统计误差。乙醇、苯酚和乙酰胺的计算值偏正，这是 GAFF/AM1-BCC 力场对极性分子的已知系统偏差，原论文的同方法结果呈现相同趋势。

方法来源

FreeSolv（Mobley & Guthrie, 2014）汇编了小分子水合自由能的实验值与自由能微扰计算值，并公开了 GAFF/AM1-BCC 参数和 GROMACS 输入文件。Scientify 阅读论文和官方输入文件，提取解耦顺序、λ 路径、软核参数、截断与 PME 设置。

官方输入面向旧版 GROMACS。Scientify 将其映射到 GROMACS 2024.4 的参数写法，例如范德华势函数切换和长程色散校正，并用 grompp 逐个状态检查。在同一协议下，甲烷、苯、甲苯、乙醇和乙酰胺的计算值与论文计算值相差 0.1–0.3 kcal/mol。

输入与结果文件

可以下载核对的文件

以下文件来自本次运行的工作区，未经改写。

适用范围

读数时需要知道的条件

生产采样为每状态 0.5 ns（乙酰胺 1 ns），短于原协议的 5 ns。
统计误差来自 MBAR，不包含力场误差。
原子电荷使用 FreeSolv 官方 AM1-BCC 参数，未重新计算。
样本为 6 个小分子，用于展示完整计算流程，不用于比较不同力场的精度。

参考资料

方法与数据来源

把计算课题交给 Scientify

描述体系和要算的性质，Scientify 配置模型和参数、按需租用 GPU 运行模拟，并交付结果与全部输入文件。

开始计算