这个主题我拖了很久才动手写。FEP这东西嘴上说“原理简单”真到了自己动手跑的时候配置项一个接一个报错一条接一条GROMACS版本更新后参数名又变网上教程还经常互相矛盾。我前前后后跑了三个体系才把整套流程摸顺。这篇我把我自己的完整操作记录整理出来从体系构建到mdp逐行解释再到批处理脚本和数据分析全部写清楚你按顺序执行就能复现。文章默认你用的是GROMACS 2020及以上版本我用2021.5测试过配置语法在2019到2023之间通用。体系用你自己手上的“蛋白配体”复合物结构就行不需要特定的靶标只要你有pdb文件和一个配体mol2/结构文件就能跟下来。1. FEP到底在算什么结合自由能的物理意义与算法定位1.1 结合自由能不是结合能很多刚转计算化学的朋友容易把“结合自由能”和“相互作用能”混为一谈。蛋白和小分子之间的范德华相互作用、静电相互作用算出来是负的值越大看着结合越强但这并不是热力学上的结合自由能。结合自由能ΔG_bind描述的是这样一个热力学过程配体从水溶液中转移到蛋白结合口袋两者形成复合物过程的自由能变化。它天然包含三个关键贡献配体脱溶剂化的代价、配体与蛋白在结合态的相互作用增益、以及结合前后整个体系的熵变化。一个配体跟口袋相互作用很强但如果它被水包围得天衣无缝要它“脱离水”进入口袋需要先付出巨大的脱溶剂能量净结果可能反而不结合。FEP方法的价值就在于它不直接猜“相互作用有多强”而是通过一系列热力学状态插值把脱溶剂、结合、熵变这些因素全部隐式地统计进来。这也是为什么FEP的精度通常比打分函数高一个量级实验相关性R²能做到0.6到0.9而普通打分函数往往不到0.3。1.2 从Zwanzig公式到lambda窗口FEP的数学骨架FEP的底层公式是Zwanzig方程ΔG -kBT ln ⟨exp[-(V1 - V0)/kBT]⟩0其中V0和V1分别是体系在态0和态1下的势能⟨⟩0表示在态0的系综下取平均。如果直接把配体从“有相互作用”瞬间变成“无相互作用”两个状态势能差巨大指数项会发散采样根本收敛不了。解决办法就是在两个端态之间插入若干个中间态每个中间态用一个耦合参数λ来定义势能V(λ) (1-λ)V0 λV1λ从0变到1体系从完整的配体状态平滑过渡到配体“消失”的状态。相邻窗口之间的自由能差用BAR或MBAR统计估算然后逐段相加得到总的自由能变化。这就是GROMACS中FEP生产模拟的核心逻辑。1.3 绝对结合自由能与相对结合自由能先选对方向开始配mdp之前一定要搞清楚自己算的是绝对还是相对结合自由能。标题里说的“小分子结合自由能”通常走的是绝对结合自由能路线分别计算配体在蛋白口袋环境中“被去耦合”从有vdW和静电相互作用变成无相互作用的自由能变化ΔG_complex以及配体在水溶液中“被去耦合”的自由能变化ΔG_solvent两者相减得到结合自由能ΔG_bind ΔG_complex - ΔG_solvent这里两个ΔG都是正值表示“拿走”配体与周围环境相互作用需要做的功。相减之后如果ΔG_bind为负说明配体在口袋中比在水中更“舒服”结合有利。相对结合自由能则是让两个相似配体在同一个结合口袋里互相转化计算的是两个配体结合自由能的差值ΔΔG_bind。相对计算对力场误差抵消更充分精度通常比绝对计算更高但一次只能回答“化合物A比化合物B亲和力高多少”这种相对排序问题题目要做的是绝对数值所以我下面全部按双体系绝对结合自由能的流程来写。2. 体系准备决定FEP成败的前处理工程2.1 蛋白结构处理与质子化状态FEP对体系的敏感性远高于普通平衡模拟。蛋白结构里的原子命名、残基命名、缺失原子、非标准残基任何一个问题都会在pdb2gmx阶段直接报错。我的标准流程是先清理pdb文件只保留蛋白和配体相关链去掉水分子、金属离子和其他配体如果这些不是研究目标。然后运行gmx pdb2gmx -f protein.pdb -o protein_processed.gro -ff gromos54a7 -water tip3p -ignh这里我推荐用-ignh让pdb2gmx重新按力场规则加氢避免原pdb中氢位置不合理导致后续能量爆炸。力场选gromos54a7还是amber99sb-ildn主要看配体参数怎么生成。我用GAFF参数生成配体时蛋白端配amber99sb-ildn或gromos54a7都可以混用但要保证后续配体topology里的原子类型和电荷来自同一套逻辑。如果蛋白端选了amber力场配体建议用GAFF2如果蛋白端用charmm36配体最好用CGenFF。混搭不是不行但会增加隐性误差来源新手不要自己给自己增加难度。蛋白中的组氨酸要特别注意质子化。pdb2gmx默认按中性pH处理大部分残基但组氨酸有HID、HIE、HIP三种状态会影响氢键网络和电荷分布。建议用gmx pdb2gmx -his交互模式逐个确认或者用PDB2PQR先算一遍质子化状态再转换。2.2 配体力场参数生成GAFF/AM1-BCC路径配体topology是FEP流程里最容易被忽视的坑。GROMACS自带的力场参数库里没有你的配体你必须自己生成一份完整的itp文件包括原子类型、电荷、键合参数。我用的是GAFF2力场加AM1-BCC电荷这是学术界做FEP最常见的组合之一。准备配体的mol2文件然后安装AmberTools用antechamber生成GAFF2类型的mol2和frcmod文件再用ACPYPE转换成GROMACS格式antechamber -i ligand.mol2 -fi mol2 -o ligand_bcc.mol2 -fo mol2 -c bcc -nc 0 -at gaff2 acpype -i ligand_bcc.mol2 -o ligand-nc 0要根据配体真实净电荷改。AM1-BCC电荷计算对初始构象敏感建议先把配体做一个简单的DFT优化或者至少用MMFF94力场优化一下避免只能量化的奇怪构象被带进电荷计算。acpype输出里面有几个关键文件ligand_GMX.top、ligand_GMX.gro、posre_ligand.itp。注意检查[ moleculetype ]的名称通常是LIG后面mdp里的couple-moltype要跟它严格一致。还有一点很多人不知道acpype生成的top文件里可能包含多个#include层级复合物topology合并时只需要把配体的itp部分include进去不要把acpype的完整top文件原样塞进蛋白的topol.top。2.3 双体系组装溶剂化、加离子与拓扑文件核对FEP需要两个独立的模拟体系一个叫复合物体系蛋白配体水离子一个叫溶剂体系配体水离子。两个体系必须分别构建、分别平衡、分别跑FEP最后两个自由能相减。复合物体系组装# 蛋白和配体结构合并 gmx editconf -f ligand_GMX.gro -o ligand_GMX.pdb cat protein_processed.gro ligand_GMX.pdb complex_pre.gro # 手动合并topol.top包括蛋白、配体、水 gmx editconf -f complex_pre.gro -o complex_box.gro -d 1.0 -bt cubic gmx solvate -cp complex_box.gro -cs spc216.gro -o complex_solv.gro -p topol.top gmx grompp -f ions.mdp -c complex_solv.gro -p topol.top -o ions.tpr gmx genion -s ions.tpr -p topol.top -o complex_ions.gro -neutral -pname NA -nname CL-d 1.0表示蛋白或配体原子到盒子边缘的最小距离是1.0 nm这是计算量、边界效应和周期性映像相互作用的折中。FEP体系建议盒子稍微大一点1.2 nm更保险因为配体“消失”过程中不排除会有构象漂移盒子太小会让配体跟自己的周期性映像相互作用。溶剂体系组装基本同理只不过把蛋白坐标去掉只留配体。注意配体可能带电genion -neutral加反离子时要留意配体本身电荷的正负不要加反了。两个体系的topol.top核对非常关键。gmx solvate和gmx genion会自动更新分子数量但前提是top文件里的#include顺序和[ molecules ]部分结构正确。我的经验是每次运行后都手动打开topol.top检查一眼确认配体在[ molecules ]里占据了一行且#include ligand.itp出现在[ system ]之前。3. mdp配置逐行拆解从平衡到FEP生产3.1 最小化与平衡阶段的mdpFEP生产之前体系必须先做能量最小化和NVT/NPT平衡。平衡阶段的mdp跟普通MD没有本质区别。我用的最小化mdp文件minimization.mdpintegrator steep nsteps 5000 emtol 1000 emstep 0.01 nstlist 10 cutoff-scheme Verlet coulombtype PME rvdw 1.0 rcoulomb 1.0平衡阶段分两步。NVT平衡300 psNPT平衡500 ps。NVT的mdpnvt.mdpintegrator md dt 0.002 nsteps 150000 nstxout-compressed 1000 nstlog 1000 tcoupl v-rescale tc-grps System tau-t 0.1 ref-t 300 nstlist 20 cutoff-scheme Verlet coulombtype PME rvdw 1.0 rcoulomb 1.0 constraints h-bondsNPT平衡在NVT基础上加压浴pcoupl berendsen tau-p 1.0 ref-p 1.0 compressibility 4.5e-5这里我用tc-grps System单组耦合。虽然有些教程建议把蛋白、配体、溶剂分开热浴但FEP生产阶段配体处于“半消失”状态自身的温度定义本来就模糊分开耦合反而容易出问题。单组耦合简单稳定实测没有问题。3.2 FEP生产阶段的核心mdp全文这是整篇文章的重头戏。FEP生产阶段我用free_energy.mdp文件内容如下integrator md dt 0.002 nsteps 2000000 nstxout-compressed 5000 nstlog 1000 nstcalcenergy 100 dhdl-frequency 100 separate-dhdl-file yes tcoupl v-rescale tc-grps System tau-t 0.1 ref-t 300 pcoupl berendsen tau-p 1.0 ref-p 1.0 compressibility 4.5e-5 nstlist 20 cutoff-scheme Verlet coulombtype PME rcoulomb 1.0 rvdw 1.0 DispCorr no constraints h-bonds free-energy yes init-lambda-state 0 couple-moltype LIG couple-lambda0 vdw-q couple-lambda1 none couple-intramol no sc-alpha 0.5 sc-coul yes sc-power 1 sc-sigma 0.3 calc-lambda-neighbors 2nsteps 2000000配合dt 0.002是4 ns这是每个lambda窗口的生产长度。lambda窗口数我后面会详细讲如果窗口少、步数短可以把nsteps加到5000000也就是10 ns。nstcalcenergy 100和dhdl-frequency 100配合每0.2 ps记录一次能量和dhdl数据。这个频率对BAR计算足够了不需要更密。3.3 软核参数、lambda与耦合方式为什么这样设free-energy yes开启FEP模块init-lambda-state 0表示grompp生成tpr时默认从第0个lambda窗口开始。couple-moltype LIG指定要对哪个分子做耦合这里的LIG必须与配体itp里的[ moleculetype ]名字一致。couple-lambda0 vdw-q和couple-lambda1 none的意思是λ0时配体有完整的范德华和静电相互作用λ1时配体与环境的相互作用完全关闭。couple-intramol no很关键它表示配体内部的键合和非键相互作用不随λ变化。如果设成yes等于把配体自身也“溶解”掉了配体内部化学键都要拆不仅毫无物理意义还会让自由能结果彻底没法解释。sc-alpha 0.5是软核势的alpha参数。不加软核时配体的vdW球在λ接近1时会因为两个原子靠得太近而产生无穷大的排斥能这就是著名的“端点奇异性”问题。软核势将范德华相互作用改为V_sc (1-λ) * α * VvdW 形式避免原子间排斥能发散。α取0.5是经验最优值既不会过度软化导致采样失真又能有效消除发散。sc-coul yes表示静电相互作用也使用软核处理。有人会设置sc-coul no只软化vdW但对带电配体来说在λ接近0或1时静电项的奇异性一样会造成能量抖动我建议带电配体一律开sc-coul yes。sc-power 1是软核势的幂次sc-sigma 0.3是软核sigma值单位是nm对应原子之间的“等效半径”。这些参数在不同教程里有细微差别但0.3是GROMACS官方测试集里的默认推荐值新手上手别乱改。4. 跑起来双体系多窗口FEP的完整运行流程4.1 复合物体系的平衡与FEP生产两个体系都按“最小化 → NVT平衡 → NPT平衡 → FEP生产”四步跑。复合物体系命令示例gmx grompp -f minimization.mdp -c complex_ions.gro -p topol.top -o em.tpr gmx mdrun -deffnm em gmx grompp -f nvt.mdp -c em.gro -p topol.top -o nvt.tpr gmx mdrun -deffnm nvt gmx grompp -f npt.mdp -c nvt.gro -p topol.top -o npt.tpr gmx mdrun -deffnm npt平衡完成后检查npt.gro的结构完整性。一个常见问题是配体在平衡过程中漂出结合口袋特别是一开始配体放的位置不好、口袋开口较大时。如果发现配体跑出口袋建议重新检查初始结构或者加一个极弱的位置约束力常数为10 kJ/mol/nm²做额外约束平衡。FEP本身不需要位置约束因为配体必须自由采样才能得到正确的结合自由能。然后进入FEP生产。我习惯为每个lambda窗口单独建目录各自跑各自的tpr。以11个lambda窗口为例mkdir -p win_{0..10} # 每个窗口生成tpr for i in $(seq 0 10); do sed s/^init-lambda-state .*/init-lambda-state $i/ free_energy.mdp free_energy_lambda$i.mdp gmx grompp -f free_energy_lambda$i.mdp -c npt.gro -p topol.top -o win_$i/tpr.tpr done # 提交运行 for i in $(seq 0 10); do cd win_$i gmx mdrun -s tpr.tpr -deffnm md md.log 21 cd .. done wait不要用sed替换多个变量时把mdp其他行弄坏init-lambda-state在整个mdp文件里只出现一次安全。这里还有一种更高效的方法只生成一个tpr用mdrun -lambda覆盖初始lambda窗口编号每个窗口读取同一个tpr文件。实测可行但新手容易把-lambda和mdp里的init-lambda-state关系搞混我宁可用上面这种老实方法多花一点磁盘空间换省心。4.2 配体水溶液的FEP生产溶剂体系配体水的FEP生产流程与复合物体系完全一样只是没有蛋白tpr生成和运行命令几乎相同。注意平衡阶段同样要先NVT后NPT不能因为体系小就跳过平衡。配体-水体系平衡很快NPT 300 ps就够了但FEP生产阶段的步数建议跟复合物体系保持一致这样两个方向的采样精度才匹配误差不会由某一方主导。溶剂体系平衡之后重点检查水的密度。GROMACS里SPC水模型在300 K、1 bar下密度应该接近0.995 g/cm³如果偏离超过2%多半是NPT平衡没跑充分或者盒子里残余应力没释放强行跑FEP会让自由能产生系统偏差。4.3 多窗口并行提交与运行状态检查FEP多窗口最适合并行因为每个lambda窗口之间相互独立可以放心地分配在不同GPU或节点上。每个窗口跑4 ns11个窗口就是44 ns的MD单个GPU大约需要10到20小时取决于体系大小。如果要加速优先减少窗口数而不是减少每个窗口的模拟时间因为每个窗口采样不充分会导致BAR误差直接变大。运行过程中要定期检查日志tail -n 30 md.log grep Step, Time md.log grep Average md.log重点关注温度是否稳定在300 K附近压力是否有剧烈波动超过几百bar能量有没有出现NaN或inf。FEP模拟中能量NaN常见原因包括初始构象有原子重叠、NPT密度异常、或者软核参数没开导致vdW发散。5. 数据分析从dhdl.xvg到结合自由能5.1 用gmx bar做BAR/MBAR统计跑完所有窗口之后复合物体系每个窗口输出一个dhdl.xvg文件。把这些文件按lambda窗口编号顺序合并给gmx barcd complex_directory gmx bar -f win_*/md_dhdl.xvg -o bar.xvg -b 2000-b 2000是指从第2000 ps开始统计跳过前2 ns的“非平衡”阶段。FEP生产模拟开始时虽然坐标来自NPT平衡但配体的耦合状态从0突变到某一个中间λ值体系需要一段时间的再平衡这部分轨迹不应该纳入自由能计算。具体跳过多少取决于收敛速度我至少跳过10%到20%的轨迹复杂体系甚至会跳过一半。gmx bar输出的结果里会给出每个相邻窗口之间的自由能差、总自由能、以及误差估计。对复合物体系记录输出结果中“Total”那一行的值记为ΔG_complex对溶剂体系同样操作得到ΔG_solvent。最终结合自由能ΔG_bind ΔG_complex - ΔG_solvent单位默认是kJ/mol论文里通常换算成kcal/mol除以4.184即可。5.2 收敛性判断什么时候结果可信一个很常见的问题是“我跑完了怎么知道结果可不可信”第一看gmx bar报告的误差。GROMACS的BAR误差是基于每个窗口的涨落和窗口间重叠度计算的如果总误差超过1.5 kJ/mol说明采样不足或窗口间隔过大结果需要谨慎对待。第二做分块分析。把轨迹分成前后两半分别算ΔG如果两半差异超过2 kJ/mol说明前半段和后半段的相空间采样不一致模拟时间不够长需要延长nsteps或者增加独立重复。第三检查相邻窗口的lambda状态重叠。BAR方法要求相邻λ窗口的势能分布有足够重叠。如果某个窗口对之间出现“无重叠”警告gmx bar会提示说明λ间距太大需要插入更多中间窗口。5.3 结果验证与物理合理性检查一个正常的绝对结合自由能计算结果通常在-20到-60 kJ/mol范围内对应Kd在微摩尔到纳摩尔级别。如果算出50 kJ/mol或者-200 kJ/mol这种极端值先别急着怀疑实验数据大概率是流程哪里错了。把自由能换算成解离常数时用Kd exp(ΔG_bind / RT)其中R 8.314 J/(mol·K)T 300 K时RT 2.494 kJ/mol。如果ΔG_bind -30 kJ/molKd exp(-12.03) ≈ 6.0×10⁻⁶ M也就是6 μM这个量级的小分子结合算比较典型。如果ΔG_bind是-8 kJ/molKd约4 mM基本是不结合或很弱的程度。另一个物理合理性检查是看单个窗口的能量分布。打开某个中间λ窗口的dhdl.xvg用gmx analyze或者直接脚本统计如果能量分布出现双峰说明体系在该λ状态存在多个稳定构象亚态单纯的BAR计算可能低估自由能垒的影响需要考虑增强采样或增加模拟长度。6. 实操中的坑与我的调参经验6.1 lambda窗口数目、步数与计算量第一次跑FEP我用过5个窗口结果误差极大bar都报告警告。后来加密到11个改进明显。对大多数小分子11到21个窗口是比较合理的区间。窗口太多会浪费算力窗口太少BAR统计的方差会爆炸。带电配体建议用17到21个窗口因为静电场相空间比纯vdW更复杂中性疏水配体用11个窗口就够。单个窗口4 ns是最低标准。如果gmx bar的误差还是超过1.5 kJ/mol先把每个窗口的步数加到6到8 ns再考虑加窗口。虽然增加窗口数也能降低误差但增加窗口数对采样充分性的改善不如延长模拟时间来得直接。6.2 端点奇异性为什么必须开软核这个我相信每个做FEP的人都踩过。第一次跑的时候我没开sc-alpha直接让λ从0变到1结果fep模拟跑到一半能量直接爆掉日志里全是NaN。原理是λ接近1时配体的vdW相互作用几乎消失但残余的微弱排斥力如果遇到两个原子距离极近势能会飙升到天文数字积分直接发散。软核势解决的就是这个问题。α参数越大势能函数被软化得越厉害常规取值0.5对于原子半径特别小的体系可以试试0.3到0.7的扫描。sc-power 1和sc-power 2的区别在于软核势在原点附近的形状默认1即可2偶尔会带来更好的数值稳定性但对普通小分子差别不大。还有一个小技巧如果你在FEP生产阶段发现某个λ窗口的模拟总是不稳定可以单独把这个窗口的init-lambda-state参数调成相邻窗口然后从相邻窗口的gro坐标继续跑。本质上是让体系先适应一个新的λ值再切回去。6.3 力场参数和体系电荷带来的隐性误差FEP对体系总电荷极度敏感。PME静电算法下体系如果有净电荷长程静电的自由能贡献会依赖于盒子大小导致结果不再具有热力学一致性。所以加离子时务必保证体系净电荷为零。复合物体系里蛋白加配体的总净电荷可能不是整数尤其配体带1电荷、蛋白某个组氨酸质子化状态变动都会改变总电荷。遇到这种情况先检查配体电荷通常已知再检查蛋白的净电荷最后确认反离子数量。另一个隐性误差来自配体力场参数选择不一致。我曾经遇到一个体系配体用GAFF生成但蛋白端用了charmm36结果ΔG_bind明显偏正后来全部改成amber力场之后结果才合理。交叉力场匹配虽然可以跑但非键参数中LJ参数来自不同力场混合规则会产生系统性偏差。对你的第一套FEP强烈建议“蛋白力场配体力场”从同一逻辑体系里选。6.4 版本差异与常见报错速查GROMACS各版本间FEP参数变化比较大很多老教程的mdp直接扔到2023版会直接报错。常遇到的就是nstdhdl在2020之后被dhdl-frequency替代写老参数会警告但不一定报错新版本干脆不认。DispCorr EnerPres如果用在新版部分测试场景下会跟软核势冲突。FEP生产mdp里我建议直接DispCorr no省得版本底层实现差异带来坑。calc-lambda-neighbors -1在老版本里表示自动新版建议写正整数2这样输出窗口前后各两个邻居的dhdl数据MBAR更稳定。常见报错“Some atoms are missing in molecule type”通常是因为couple-moltype名称跟配体itp里的分子名对不上。这个好排查打开配体的itp文件看[ moleculetype ]下面的名字改到mdp里就行。还有报错“ERROR: The group System has no atoms”这类说明拓扑文件或结构文件里分子数量为零或者index文件没对上多半是solvate或genion步骤出错回头检查topol.top。遇到“NaN in bond, angle”这类报错先确认是不是初始结构原子重叠太多也检查平衡阶段是否真的跑完了NPT没完成就直接在FEP生产里跑很容易出问题。根据我个人经验FEP整套流程跑顺一次之后换体系再跑就轻松多了主要工作量会转移到配体参数生成和体系平衡上。这篇教程覆盖的是绝对结合自由能的最基础协议如果你后续要做相对结合自由能只需要把couple-lambda0和couple-lambda1的设计从“到none”改成“配体A到配体B”的双拓扑结构原理和mdp配置完全同源。先用这套标准流程拿到第一个可信的ΔG_bind再考虑更复杂的方案路会稳很多。 GROMACS FEP绝对结合自由能计算全流程：从体系构建到数据分析
这个主题我拖了很久才动手写。FEP这东西嘴上说“原理简单”真到了自己动手跑的时候配置项一个接一个报错一条接一条GROMACS版本更新后参数名又变网上教程还经常互相矛盾。我前前后后跑了三个体系才把整套流程摸顺。这篇我把我自己的完整操作记录整理出来从体系构建到mdp逐行解释再到批处理脚本和数据分析全部写清楚你按顺序执行就能复现。文章默认你用的是GROMACS 2020及以上版本我用2021.5测试过配置语法在2019到2023之间通用。体系用你自己手上的“蛋白配体”复合物结构就行不需要特定的靶标只要你有pdb文件和一个配体mol2/结构文件就能跟下来。1. FEP到底在算什么结合自由能的物理意义与算法定位1.1 结合自由能不是结合能很多刚转计算化学的朋友容易把“结合自由能”和“相互作用能”混为一谈。蛋白和小分子之间的范德华相互作用、静电相互作用算出来是负的值越大看着结合越强但这并不是热力学上的结合自由能。结合自由能ΔG_bind描述的是这样一个热力学过程配体从水溶液中转移到蛋白结合口袋两者形成复合物过程的自由能变化。它天然包含三个关键贡献配体脱溶剂化的代价、配体与蛋白在结合态的相互作用增益、以及结合前后整个体系的熵变化。一个配体跟口袋相互作用很强但如果它被水包围得天衣无缝要它“脱离水”进入口袋需要先付出巨大的脱溶剂能量净结果可能反而不结合。FEP方法的价值就在于它不直接猜“相互作用有多强”而是通过一系列热力学状态插值把脱溶剂、结合、熵变这些因素全部隐式地统计进来。这也是为什么FEP的精度通常比打分函数高一个量级实验相关性R²能做到0.6到0.9而普通打分函数往往不到0.3。1.2 从Zwanzig公式到lambda窗口FEP的数学骨架FEP的底层公式是Zwanzig方程ΔG -kBT ln ⟨exp[-(V1 - V0)/kBT]⟩0其中V0和V1分别是体系在态0和态1下的势能⟨⟩0表示在态0的系综下取平均。如果直接把配体从“有相互作用”瞬间变成“无相互作用”两个状态势能差巨大指数项会发散采样根本收敛不了。解决办法就是在两个端态之间插入若干个中间态每个中间态用一个耦合参数λ来定义势能V(λ) (1-λ)V0 λV1λ从0变到1体系从完整的配体状态平滑过渡到配体“消失”的状态。相邻窗口之间的自由能差用BAR或MBAR统计估算然后逐段相加得到总的自由能变化。这就是GROMACS中FEP生产模拟的核心逻辑。1.3 绝对结合自由能与相对结合自由能先选对方向开始配mdp之前一定要搞清楚自己算的是绝对还是相对结合自由能。标题里说的“小分子结合自由能”通常走的是绝对结合自由能路线分别计算配体在蛋白口袋环境中“被去耦合”从有vdW和静电相互作用变成无相互作用的自由能变化ΔG_complex以及配体在水溶液中“被去耦合”的自由能变化ΔG_solvent两者相减得到结合自由能ΔG_bind ΔG_complex - ΔG_solvent这里两个ΔG都是正值表示“拿走”配体与周围环境相互作用需要做的功。相减之后如果ΔG_bind为负说明配体在口袋中比在水中更“舒服”结合有利。相对结合自由能则是让两个相似配体在同一个结合口袋里互相转化计算的是两个配体结合自由能的差值ΔΔG_bind。相对计算对力场误差抵消更充分精度通常比绝对计算更高但一次只能回答“化合物A比化合物B亲和力高多少”这种相对排序问题题目要做的是绝对数值所以我下面全部按双体系绝对结合自由能的流程来写。2. 体系准备决定FEP成败的前处理工程2.1 蛋白结构处理与质子化状态FEP对体系的敏感性远高于普通平衡模拟。蛋白结构里的原子命名、残基命名、缺失原子、非标准残基任何一个问题都会在pdb2gmx阶段直接报错。我的标准流程是先清理pdb文件只保留蛋白和配体相关链去掉水分子、金属离子和其他配体如果这些不是研究目标。然后运行gmx pdb2gmx -f protein.pdb -o protein_processed.gro -ff gromos54a7 -water tip3p -ignh这里我推荐用-ignh让pdb2gmx重新按力场规则加氢避免原pdb中氢位置不合理导致后续能量爆炸。力场选gromos54a7还是amber99sb-ildn主要看配体参数怎么生成。我用GAFF参数生成配体时蛋白端配amber99sb-ildn或gromos54a7都可以混用但要保证后续配体topology里的原子类型和电荷来自同一套逻辑。如果蛋白端选了amber力场配体建议用GAFF2如果蛋白端用charmm36配体最好用CGenFF。混搭不是不行但会增加隐性误差来源新手不要自己给自己增加难度。蛋白中的组氨酸要特别注意质子化。pdb2gmx默认按中性pH处理大部分残基但组氨酸有HID、HIE、HIP三种状态会影响氢键网络和电荷分布。建议用gmx pdb2gmx -his交互模式逐个确认或者用PDB2PQR先算一遍质子化状态再转换。2.2 配体力场参数生成GAFF/AM1-BCC路径配体topology是FEP流程里最容易被忽视的坑。GROMACS自带的力场参数库里没有你的配体你必须自己生成一份完整的itp文件包括原子类型、电荷、键合参数。我用的是GAFF2力场加AM1-BCC电荷这是学术界做FEP最常见的组合之一。准备配体的mol2文件然后安装AmberTools用antechamber生成GAFF2类型的mol2和frcmod文件再用ACPYPE转换成GROMACS格式antechamber -i ligand.mol2 -fi mol2 -o ligand_bcc.mol2 -fo mol2 -c bcc -nc 0 -at gaff2 acpype -i ligand_bcc.mol2 -o ligand-nc 0要根据配体真实净电荷改。AM1-BCC电荷计算对初始构象敏感建议先把配体做一个简单的DFT优化或者至少用MMFF94力场优化一下避免只能量化的奇怪构象被带进电荷计算。acpype输出里面有几个关键文件ligand_GMX.top、ligand_GMX.gro、posre_ligand.itp。注意检查[ moleculetype ]的名称通常是LIG后面mdp里的couple-moltype要跟它严格一致。还有一点很多人不知道acpype生成的top文件里可能包含多个#include层级复合物topology合并时只需要把配体的itp部分include进去不要把acpype的完整top文件原样塞进蛋白的topol.top。2.3 双体系组装溶剂化、加离子与拓扑文件核对FEP需要两个独立的模拟体系一个叫复合物体系蛋白配体水离子一个叫溶剂体系配体水离子。两个体系必须分别构建、分别平衡、分别跑FEP最后两个自由能相减。复合物体系组装# 蛋白和配体结构合并 gmx editconf -f ligand_GMX.gro -o ligand_GMX.pdb cat protein_processed.gro ligand_GMX.pdb complex_pre.gro # 手动合并topol.top包括蛋白、配体、水 gmx editconf -f complex_pre.gro -o complex_box.gro -d 1.0 -bt cubic gmx solvate -cp complex_box.gro -cs spc216.gro -o complex_solv.gro -p topol.top gmx grompp -f ions.mdp -c complex_solv.gro -p topol.top -o ions.tpr gmx genion -s ions.tpr -p topol.top -o complex_ions.gro -neutral -pname NA -nname CL-d 1.0表示蛋白或配体原子到盒子边缘的最小距离是1.0 nm这是计算量、边界效应和周期性映像相互作用的折中。FEP体系建议盒子稍微大一点1.2 nm更保险因为配体“消失”过程中不排除会有构象漂移盒子太小会让配体跟自己的周期性映像相互作用。溶剂体系组装基本同理只不过把蛋白坐标去掉只留配体。注意配体可能带电genion -neutral加反离子时要留意配体本身电荷的正负不要加反了。两个体系的topol.top核对非常关键。gmx solvate和gmx genion会自动更新分子数量但前提是top文件里的#include顺序和[ molecules ]部分结构正确。我的经验是每次运行后都手动打开topol.top检查一眼确认配体在[ molecules ]里占据了一行且#include ligand.itp出现在[ system ]之前。3. mdp配置逐行拆解从平衡到FEP生产3.1 最小化与平衡阶段的mdpFEP生产之前体系必须先做能量最小化和NVT/NPT平衡。平衡阶段的mdp跟普通MD没有本质区别。我用的最小化mdp文件minimization.mdpintegrator steep nsteps 5000 emtol 1000 emstep 0.01 nstlist 10 cutoff-scheme Verlet coulombtype PME rvdw 1.0 rcoulomb 1.0平衡阶段分两步。NVT平衡300 psNPT平衡500 ps。NVT的mdpnvt.mdpintegrator md dt 0.002 nsteps 150000 nstxout-compressed 1000 nstlog 1000 tcoupl v-rescale tc-grps System tau-t 0.1 ref-t 300 nstlist 20 cutoff-scheme Verlet coulombtype PME rvdw 1.0 rcoulomb 1.0 constraints h-bondsNPT平衡在NVT基础上加压浴pcoupl berendsen tau-p 1.0 ref-p 1.0 compressibility 4.5e-5这里我用tc-grps System单组耦合。虽然有些教程建议把蛋白、配体、溶剂分开热浴但FEP生产阶段配体处于“半消失”状态自身的温度定义本来就模糊分开耦合反而容易出问题。单组耦合简单稳定实测没有问题。3.2 FEP生产阶段的核心mdp全文这是整篇文章的重头戏。FEP生产阶段我用free_energy.mdp文件内容如下integrator md dt 0.002 nsteps 2000000 nstxout-compressed 5000 nstlog 1000 nstcalcenergy 100 dhdl-frequency 100 separate-dhdl-file yes tcoupl v-rescale tc-grps System tau-t 0.1 ref-t 300 pcoupl berendsen tau-p 1.0 ref-p 1.0 compressibility 4.5e-5 nstlist 20 cutoff-scheme Verlet coulombtype PME rcoulomb 1.0 rvdw 1.0 DispCorr no constraints h-bonds free-energy yes init-lambda-state 0 couple-moltype LIG couple-lambda0 vdw-q couple-lambda1 none couple-intramol no sc-alpha 0.5 sc-coul yes sc-power 1 sc-sigma 0.3 calc-lambda-neighbors 2nsteps 2000000配合dt 0.002是4 ns这是每个lambda窗口的生产长度。lambda窗口数我后面会详细讲如果窗口少、步数短可以把nsteps加到5000000也就是10 ns。nstcalcenergy 100和dhdl-frequency 100配合每0.2 ps记录一次能量和dhdl数据。这个频率对BAR计算足够了不需要更密。3.3 软核参数、lambda与耦合方式为什么这样设free-energy yes开启FEP模块init-lambda-state 0表示grompp生成tpr时默认从第0个lambda窗口开始。couple-moltype LIG指定要对哪个分子做耦合这里的LIG必须与配体itp里的[ moleculetype ]名字一致。couple-lambda0 vdw-q和couple-lambda1 none的意思是λ0时配体有完整的范德华和静电相互作用λ1时配体与环境的相互作用完全关闭。couple-intramol no很关键它表示配体内部的键合和非键相互作用不随λ变化。如果设成yes等于把配体自身也“溶解”掉了配体内部化学键都要拆不仅毫无物理意义还会让自由能结果彻底没法解释。sc-alpha 0.5是软核势的alpha参数。不加软核时配体的vdW球在λ接近1时会因为两个原子靠得太近而产生无穷大的排斥能这就是著名的“端点奇异性”问题。软核势将范德华相互作用改为V_sc (1-λ) * α * VvdW 形式避免原子间排斥能发散。α取0.5是经验最优值既不会过度软化导致采样失真又能有效消除发散。sc-coul yes表示静电相互作用也使用软核处理。有人会设置sc-coul no只软化vdW但对带电配体来说在λ接近0或1时静电项的奇异性一样会造成能量抖动我建议带电配体一律开sc-coul yes。sc-power 1是软核势的幂次sc-sigma 0.3是软核sigma值单位是nm对应原子之间的“等效半径”。这些参数在不同教程里有细微差别但0.3是GROMACS官方测试集里的默认推荐值新手上手别乱改。4. 跑起来双体系多窗口FEP的完整运行流程4.1 复合物体系的平衡与FEP生产两个体系都按“最小化 → NVT平衡 → NPT平衡 → FEP生产”四步跑。复合物体系命令示例gmx grompp -f minimization.mdp -c complex_ions.gro -p topol.top -o em.tpr gmx mdrun -deffnm em gmx grompp -f nvt.mdp -c em.gro -p topol.top -o nvt.tpr gmx mdrun -deffnm nvt gmx grompp -f npt.mdp -c nvt.gro -p topol.top -o npt.tpr gmx mdrun -deffnm npt平衡完成后检查npt.gro的结构完整性。一个常见问题是配体在平衡过程中漂出结合口袋特别是一开始配体放的位置不好、口袋开口较大时。如果发现配体跑出口袋建议重新检查初始结构或者加一个极弱的位置约束力常数为10 kJ/mol/nm²做额外约束平衡。FEP本身不需要位置约束因为配体必须自由采样才能得到正确的结合自由能。然后进入FEP生产。我习惯为每个lambda窗口单独建目录各自跑各自的tpr。以11个lambda窗口为例mkdir -p win_{0..10} # 每个窗口生成tpr for i in $(seq 0 10); do sed s/^init-lambda-state .*/init-lambda-state $i/ free_energy.mdp free_energy_lambda$i.mdp gmx grompp -f free_energy_lambda$i.mdp -c npt.gro -p topol.top -o win_$i/tpr.tpr done # 提交运行 for i in $(seq 0 10); do cd win_$i gmx mdrun -s tpr.tpr -deffnm md md.log 21 cd .. done wait不要用sed替换多个变量时把mdp其他行弄坏init-lambda-state在整个mdp文件里只出现一次安全。这里还有一种更高效的方法只生成一个tpr用mdrun -lambda覆盖初始lambda窗口编号每个窗口读取同一个tpr文件。实测可行但新手容易把-lambda和mdp里的init-lambda-state关系搞混我宁可用上面这种老实方法多花一点磁盘空间换省心。4.2 配体水溶液的FEP生产溶剂体系配体水的FEP生产流程与复合物体系完全一样只是没有蛋白tpr生成和运行命令几乎相同。注意平衡阶段同样要先NVT后NPT不能因为体系小就跳过平衡。配体-水体系平衡很快NPT 300 ps就够了但FEP生产阶段的步数建议跟复合物体系保持一致这样两个方向的采样精度才匹配误差不会由某一方主导。溶剂体系平衡之后重点检查水的密度。GROMACS里SPC水模型在300 K、1 bar下密度应该接近0.995 g/cm³如果偏离超过2%多半是NPT平衡没跑充分或者盒子里残余应力没释放强行跑FEP会让自由能产生系统偏差。4.3 多窗口并行提交与运行状态检查FEP多窗口最适合并行因为每个lambda窗口之间相互独立可以放心地分配在不同GPU或节点上。每个窗口跑4 ns11个窗口就是44 ns的MD单个GPU大约需要10到20小时取决于体系大小。如果要加速优先减少窗口数而不是减少每个窗口的模拟时间因为每个窗口采样不充分会导致BAR误差直接变大。运行过程中要定期检查日志tail -n 30 md.log grep Step, Time md.log grep Average md.log重点关注温度是否稳定在300 K附近压力是否有剧烈波动超过几百bar能量有没有出现NaN或inf。FEP模拟中能量NaN常见原因包括初始构象有原子重叠、NPT密度异常、或者软核参数没开导致vdW发散。5. 数据分析从dhdl.xvg到结合自由能5.1 用gmx bar做BAR/MBAR统计跑完所有窗口之后复合物体系每个窗口输出一个dhdl.xvg文件。把这些文件按lambda窗口编号顺序合并给gmx barcd complex_directory gmx bar -f win_*/md_dhdl.xvg -o bar.xvg -b 2000-b 2000是指从第2000 ps开始统计跳过前2 ns的“非平衡”阶段。FEP生产模拟开始时虽然坐标来自NPT平衡但配体的耦合状态从0突变到某一个中间λ值体系需要一段时间的再平衡这部分轨迹不应该纳入自由能计算。具体跳过多少取决于收敛速度我至少跳过10%到20%的轨迹复杂体系甚至会跳过一半。gmx bar输出的结果里会给出每个相邻窗口之间的自由能差、总自由能、以及误差估计。对复合物体系记录输出结果中“Total”那一行的值记为ΔG_complex对溶剂体系同样操作得到ΔG_solvent。最终结合自由能ΔG_bind ΔG_complex - ΔG_solvent单位默认是kJ/mol论文里通常换算成kcal/mol除以4.184即可。5.2 收敛性判断什么时候结果可信一个很常见的问题是“我跑完了怎么知道结果可不可信”第一看gmx bar报告的误差。GROMACS的BAR误差是基于每个窗口的涨落和窗口间重叠度计算的如果总误差超过1.5 kJ/mol说明采样不足或窗口间隔过大结果需要谨慎对待。第二做分块分析。把轨迹分成前后两半分别算ΔG如果两半差异超过2 kJ/mol说明前半段和后半段的相空间采样不一致模拟时间不够长需要延长nsteps或者增加独立重复。第三检查相邻窗口的lambda状态重叠。BAR方法要求相邻λ窗口的势能分布有足够重叠。如果某个窗口对之间出现“无重叠”警告gmx bar会提示说明λ间距太大需要插入更多中间窗口。5.3 结果验证与物理合理性检查一个正常的绝对结合自由能计算结果通常在-20到-60 kJ/mol范围内对应Kd在微摩尔到纳摩尔级别。如果算出50 kJ/mol或者-200 kJ/mol这种极端值先别急着怀疑实验数据大概率是流程哪里错了。把自由能换算成解离常数时用Kd exp(ΔG_bind / RT)其中R 8.314 J/(mol·K)T 300 K时RT 2.494 kJ/mol。如果ΔG_bind -30 kJ/molKd exp(-12.03) ≈ 6.0×10⁻⁶ M也就是6 μM这个量级的小分子结合算比较典型。如果ΔG_bind是-8 kJ/molKd约4 mM基本是不结合或很弱的程度。另一个物理合理性检查是看单个窗口的能量分布。打开某个中间λ窗口的dhdl.xvg用gmx analyze或者直接脚本统计如果能量分布出现双峰说明体系在该λ状态存在多个稳定构象亚态单纯的BAR计算可能低估自由能垒的影响需要考虑增强采样或增加模拟长度。6. 实操中的坑与我的调参经验6.1 lambda窗口数目、步数与计算量第一次跑FEP我用过5个窗口结果误差极大bar都报告警告。后来加密到11个改进明显。对大多数小分子11到21个窗口是比较合理的区间。窗口太多会浪费算力窗口太少BAR统计的方差会爆炸。带电配体建议用17到21个窗口因为静电场相空间比纯vdW更复杂中性疏水配体用11个窗口就够。单个窗口4 ns是最低标准。如果gmx bar的误差还是超过1.5 kJ/mol先把每个窗口的步数加到6到8 ns再考虑加窗口。虽然增加窗口数也能降低误差但增加窗口数对采样充分性的改善不如延长模拟时间来得直接。6.2 端点奇异性为什么必须开软核这个我相信每个做FEP的人都踩过。第一次跑的时候我没开sc-alpha直接让λ从0变到1结果fep模拟跑到一半能量直接爆掉日志里全是NaN。原理是λ接近1时配体的vdW相互作用几乎消失但残余的微弱排斥力如果遇到两个原子距离极近势能会飙升到天文数字积分直接发散。软核势解决的就是这个问题。α参数越大势能函数被软化得越厉害常规取值0.5对于原子半径特别小的体系可以试试0.3到0.7的扫描。sc-power 1和sc-power 2的区别在于软核势在原点附近的形状默认1即可2偶尔会带来更好的数值稳定性但对普通小分子差别不大。还有一个小技巧如果你在FEP生产阶段发现某个λ窗口的模拟总是不稳定可以单独把这个窗口的init-lambda-state参数调成相邻窗口然后从相邻窗口的gro坐标继续跑。本质上是让体系先适应一个新的λ值再切回去。6.3 力场参数和体系电荷带来的隐性误差FEP对体系总电荷极度敏感。PME静电算法下体系如果有净电荷长程静电的自由能贡献会依赖于盒子大小导致结果不再具有热力学一致性。所以加离子时务必保证体系净电荷为零。复合物体系里蛋白加配体的总净电荷可能不是整数尤其配体带1电荷、蛋白某个组氨酸质子化状态变动都会改变总电荷。遇到这种情况先检查配体电荷通常已知再检查蛋白的净电荷最后确认反离子数量。另一个隐性误差来自配体力场参数选择不一致。我曾经遇到一个体系配体用GAFF生成但蛋白端用了charmm36结果ΔG_bind明显偏正后来全部改成amber力场之后结果才合理。交叉力场匹配虽然可以跑但非键参数中LJ参数来自不同力场混合规则会产生系统性偏差。对你的第一套FEP强烈建议“蛋白力场配体力场”从同一逻辑体系里选。6.4 版本差异与常见报错速查GROMACS各版本间FEP参数变化比较大很多老教程的mdp直接扔到2023版会直接报错。常遇到的就是nstdhdl在2020之后被dhdl-frequency替代写老参数会警告但不一定报错新版本干脆不认。DispCorr EnerPres如果用在新版部分测试场景下会跟软核势冲突。FEP生产mdp里我建议直接DispCorr no省得版本底层实现差异带来坑。calc-lambda-neighbors -1在老版本里表示自动新版建议写正整数2这样输出窗口前后各两个邻居的dhdl数据MBAR更稳定。常见报错“Some atoms are missing in molecule type”通常是因为couple-moltype名称跟配体itp里的分子名对不上。这个好排查打开配体的itp文件看[ moleculetype ]下面的名字改到mdp里就行。还有报错“ERROR: The group System has no atoms”这类说明拓扑文件或结构文件里分子数量为零或者index文件没对上多半是solvate或genion步骤出错回头检查topol.top。遇到“NaN in bond, angle”这类报错先确认是不是初始结构原子重叠太多也检查平衡阶段是否真的跑完了NPT没完成就直接在FEP生产里跑很容易出问题。根据我个人经验FEP整套流程跑顺一次之后换体系再跑就轻松多了主要工作量会转移到配体参数生成和体系平衡上。这篇教程覆盖的是绝对结合自由能的最基础协议如果你后续要做相对结合自由能只需要把couple-lambda0和couple-lambda1的设计从“到none”改成“配体A到配体B”的双拓扑结构原理和mdp配置完全同源。先用这套标准流程拿到第一个可信的ΔG_bind再考虑更复杂的方案路会稳很多。