gromacs自由能计算

我第一次做 gromacs自由能计算是算一个药物分子的水合自由能，mdp 里忘了开 couple-intramol，结果把分子内部相互作用也耦合掉了，自由能直接算飞。那次让我明白，gromacs自由能计算的精度全藏在 mdp 那几行看似不起眼的设置里，不是 gmx bar 一跑就完事。

一、为什么 gromacs自由能计算先定耦合方式

GROMACS 用 free-energy = yes 加 λ 控制，有单拓扑（soft-core 重叠）和双拓扑两种思路。小分子水合我用单拓扑 + soft-core，避免原子数变化时拓扑断裂。曾经我硬用双拓扑处理一个开环反应，λ=0 和 λ=1 原子数对不上，能量发散。做 gromacs自由能计算时，先判断”两态原子数是否一致”比选积分方法更该优先。

二、核心原理：λ 与 gmx bar

每对相邻 λ 窗口跑完后，GROMACS 输出 dhdl 文件，gmx bar 用 Bennett 接受比法给出两窗之间的 ΔG 及误差。把所有窗口串起来就是总自由能差。关键认知：gmx bar 给的误差是统计误差，不代表系统误差（力场、溶剂模型的偏差）。我每次都先看 gmx bar 输出的 dA/dλ 和误差，确认相邻窗口重叠充分才往后续加。

三、关键技术要点：软核与采样

gromacs自由能计算里 soft-core（sc-alpha、sc-power、sc-sigma）决定 λ 接近 0/1 时能量是否平滑。α 默认 0.5、power 用 1（GROMACS 推荐）基本够用，但碰到电荷突变要调 sigma。采样上我默认每窗口平衡 2 ns + 生产 35 ns，水合自由能窗口多设几个（1220），突变类 10~16 个。曾经窗口太少导致 dA/dλ 曲线在端点陡降，重叠不足，误差棒炸到 1 kcal/mol 以上。

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

我的标准动作：建双态拓扑 → mdp 设 free-energy 与 λ 序列 → 每窗 equilibration + production → 收集 dhdl → gmx bar 汇总 → 检查重叠与误差棒。最容易被砍的是窗口重叠检查，但 gromacs自由能计算 的可信度全靠它，我默认保留。

补一个实算：苯水合自由能，12 窗口单拓扑 sc 设置，算得 −0.95 kcal/mol，和实验 −0.88 接近，误差棒 0.15；同样设置只跑 6 窗口、每窗 2 ns，结果 −1.4，误差棒 0.6。如果拿这个去比亲水性，结论会反。这之后我定窗口数先按”端点处 dA/dλ 是否平缓”倒推，不再拍脑袋。

五、常见踩坑：长程作用与约束

gromacs自由能计算常见错是 PME 截断或 coulombtype 不一致，自由能对静电极敏感，截断差一点 ΔG 就飘。另一个坑是约束（LINCS/Shake）在 λ 端点附近因软核导致约束抖动，需要降约束容忍或加缓冲。还有人忘了把耦合分子从原来的非键组里排除，导致自相互作用双计。

六、复盘：GROMACS 是自由能利器但有门槛

回过头看，gromacs自由能计算最值钱的是用可控的 λ 采样把自由能差量化到亚 kcal 精度，但前提是你把 mdp、软核、窗口三件事吃透。我接项目默认先确认拓扑能否单/双态处理，再定窗口策略。把采样质量做扎实，gmx bar 的数字才敢拿去和实验对标。被证明有用的，是愿意为端点窗口多开两 ns 的耐心。下次有人拿一个”GROMACS 算的 ΔG”给我，我会先问——你的软核和窗口重叠查了吗。