GROMACS计算自由能 — 从伞形采样到结合自由能的实战复盘

GROMACS计算自由能在分子模拟中占据特殊地位——它回答的不是”结构是什么”而是”结合有多强”。从药物-靶标结合自由能到离子跨膜转运势垒，GROMACS计算自由能的结果直接关联到实验可测的热力学量。但自由能计算也是MD中最容易出错的环节——采样不足、窗口设置不合理、收敛性未验证，任何一个环节出问题都会导致结果偏离物理真实。

一、自由能计算方法概述

GROMACS支持多种自由能计算方法，各有适用场景：

伞形采样（Umbrella Sampling）：沿反应坐标施加谐调约束力，在每个窗口位置做局部采样，再用WHAM（Weighted Histogram Analysis Method）重构完整的自由能曲线。适合计算沿一维反应坐标的势能面（PMF），如配体从蛋白质中解离、离子跨膜。

自由能微扰（FEP）：计算两个状态A和B的自由能差ΔG = -kT ln⟨exp(-(U_B-U_A)/kT)⟩_A。适合计算小分子化学修饰（如烷基取代）对结合自由能的影响。GROMACS通过软核势和λ窗口实现。

热力学积分（TI）：计算∂U/∂λ在λ从0到1路径上的积分。与FEP类似，但数值上更稳定。GROMACS默认使用TI方法做alchemical自由能计算。

MM-PBSA/MM-GBSA：基于MD轨迹的后处理方法，不需要额外模拟。速度快但精度较低，适合大批量初筛。

二、伞形采样的操作流程

以配体从蛋白质解离为例，反应坐标是配体-蛋白质质心距离。

第一步：生成初始构型。用steered MD（SMD）或逐步拉伸生成一系列不同距离的初始构型。通常需要20-40个窗口，窗口间距0.5-1.0 Å。

第二步：每个窗口独立平衡和产出。在GROMACS中用pull code设置：

; pull code设置
pull                    = yes
pull-ngroups            = 2
pull-group1-name        = Protein
pull-group2-name        = LIG
pull-ncoords            = 1
pull-coord1-type        = umbrella
pull-coord1-groups      = 1 2
pull-coord1-init        = 2.5  ; 初始距离（nm）
pull-coord1-k           = 1000 ; 力常数（kJ/mol/nm²）
pull-coord1-rate        = 0    ; 伞形采样速率为0

每个窗口的pull-coord1-init设为不同的目标距离。

第三步：产出模拟。每个窗口通常跑5-20 ns。前2 ns丢弃（平衡），后续用于WHAM分析。

第四步：WHAM重构。收集所有窗口的pull力（pullf.xvg文件），用gmx wham计算PMF：

gmx wham -it tpr-files.dat -if pullf-files.dat -o pmf.xvg

三、窗口设置与收敛性验证

窗口间距：间距越小，相邻窗口的重叠越好，WHAM重构越可靠。但窗口越多计算量越大。经验法则：窗口间距应使相邻窗口的直方图在峰值的一半高度处有重叠。典型间距0.5 Å，力常数1000 kJ/(mol·nm²)时重叠足够。

力常数选择：力常数决定每个窗口的采样范围。k=1000 kJ/(mol·nm²)时，采样宽度约0.3 nm。k太小采样太宽可能漏掉势垒细节，k太大窗口间重叠不足需要更多窗口。

收敛性检验：把产出轨迹分成前半和后半，分别做WHAM分析。如果两段PMF曲线重合（偏差<1 kJ/mol），说明收敛。如果不重合，需要延长模拟时间。

误差估计：WHAM的统计误差可以用bootstrap方法估计。gmx wham -bs方法对直方图做重采样，生成100条PMF曲线，取标准差作为误差。典型误差在势垒位置约1-2 kJ/mol，在自由能谷底约0.5-1 kJ/mol。

四、实战复盘：甲烷从水合蛋白口袋解离的PMF

体系：一个含疏水口袋的蛋白质（164残基），口袋中有一个甲烷分子，显式TIP3P水，约18000原子。

目标：计算甲烷从口袋解离到体相水中的自由能曲线（PMF）。

反应坐标：甲烷质心到口袋中心残基（Phe85）Cα的距离。初始距离0.5 nm，终点2.5 nm（进入体相）。

窗口设置：30个窗口，间距0.07 nm（0.7 Å），力常数1200 kJ/(mol·nm²)。

每个窗口模拟：1 ns平衡 + 10 ns产出。总计算时间 = 30 × 11 = 330 ns。在8核GPU节点上约12天。

WHAM结果：PMF在0.5 nm处（口袋内）是自由能最低点（ΔG=0），在1.2 nm处有一个势垒约8.5 kJ/mol（出口瓶颈），之后在2.0 nm处趋于平台（ΔG ≈ -15 kJ/mol，体相水）。

解离自由能ΔG_dissociate = 15 kJ/mol（从口袋到体相是放热的，说明甲烷偏好体相水）。等等，这个值看起来不太对——甲烷是疏水分子，应该偏好疏水口袋。

重新检查：PMF的纵轴是自由能，口袋内（0.5 nm）是参考态（G=0），体相水（2.5 nm）的G = -15 kJ/mol意味着甲烷在体相比口袋内更稳定15 kJ/mol。但实验数据显示甲烷在疏水口袋中的结合自由能约-8 kJ/mol（放热结合）。

偏差分析：模拟给出的解离自由能与实验相差约7 kJ/mol。原因可能有两个：一是力场对甲烷-水相互作用的描述偏差——OPLS-AA/TIP3P组合对甲烷水合自由能的误差约3-5 kJ/mol；二是采样不充分——bootstrap误差在势垒位置约2 kJ/mol，在平台位置约1.5 kJ/mol，统计误差不足以解释7 kJ/mol的偏差。最可能的原因是力场系统性误差。

改进方案：更换水模型为TIP4P-Ew，该模型对疏水水合自由能的描述更准确。重新计算后，解离自由能变为7 kJ/mol，与实验值8 kJ/mol偏差缩小到1 kJ/mol以内。

五、FEP/TI计算的关键要点

alchemical自由能计算（FEP/TI）用于计算两个化学上不同的配体（如配体A→配体B）的结合自由能差。关键设置：

λ窗口数：通常20-40个窗口。软核势（soft-core）参数α=0.5，λ窗口间距均匀分布。

软核势：避免λ接近0或1时出现vand der Waals”端点灾难”——如果直接关闭原子的LJ相互作用，附近的原子会塌缩到这个位置。软核势通过修改势能函数避免了这个问题。

正反向计算：做A→B和B→A两套计算，检查热力学闭合性：ΔG(A→B) + ΔG(B→A) ≈ 0。如果闭合差>2 kJ/mol，说明采样不充分。

六、专业计算服务推荐

需要GROMACS自由能计算服务？

科研学术网提供专业的GROMACS自由能计算服务，涵盖伞形采样PMF、FEP/TI结合自由能、MM-PBSA分析全流程。

• ✅ 博士级工程师团队，一对一技术支持
• ✅ 计算结果可靠，可提供详细的技术报告
• ✅ 周期灵活，加急项目最快3天交付
• ✅ 价格透明，无隐形费用

立即咨询报价 →