做药物设计或者先导化合物优化的人迟早都会碰到同一个问题手里有十几个看起来都不错的小分子到底哪个跟靶标结合得更牢实验上可以跑SPR、ITC但成本高、周期长没法天天筛。计算上可以打分函数粗筛但精度有限真正要拿得出手的排名还是得靠自由能计算。而在自由能计算方法里GROMACS的FEP自由能微扰是目前学术界和工业界都用得最广泛、资料最全、踩坑教程相对最多的一条路。这篇文章就把我用GROMACS做FEP计算小分子结合自由能的完整流程包括热力学循环怎么搭、为什么要这么搭、mdp配置文件里每个关键参数到底在干什么、跑完以后数据怎么分析全部拆开讲清楚。文章最后会附带可以直接拿来改的完整mdp模板以及我实际跑模拟时踩过的坑和排查思路。适合有一定分子动力学基础、想上手自由能计算但又不想去啃一堆晦涩理论文献的读者。1. FEP到底在算什么一张图看懂热力学循环先明确一件事FEP算的不是“结合能”算的是“结合自由能差”。结合能是单个构象的能量差结合自由能则是体系在热力学系综下所有微观状态的统计平均所以它天然包含了熵效应、溶剂效应、构象涨落。这也是为什么FEP算出来的数值比单纯用MM/PBSA或者打分函数更接近实验值。1.1 自由能微扰的数学本质FEP的理论基础其实不复杂就是用统计力学里的指数平均公式ΔG -kT ln ⟨exp(-(U_B - U_A)/kT)⟩_A这个公式说的是如果我想知道从状态A变到状态B的自由能差那么我在A状态里做大量采样把每个构象下A和B两个状态的能量差拿出来做一个带指数的Boltzmann平均就能得到ΔG。道理不难但直接做会有一个致命问题如果A和B差别很大两个状态的构象空间重叠很少指数平均会严重发散算出来的自由能完全不可信。这就好比你想比较北京和上海两地的房价差异但只在北京随机拍了几间房然后又拿这几间房的“北京-上海价差”去推整体差异——样本根本覆盖不了上海的房价分布。所以实际做FEP时一定要引入中间态。1.2 耦合参数λ的引入GROMACS的做法是在哈密顿量里塞进一个耦合参数λ让体系的势能函数跟λ挂钩U(λ) (1-λ)U_0 λU_1当λ0时体系是完整的物理态比如配体与蛋白正常相互作用当λ1时体系是非物理态比如配体与环境的相互作用被完全关闭。中间λ从0渐变到1的过程就是配体从“完整存在”慢慢变成“不存在”的假想路径。每跑一个λ窗口就能得到该窗口下∂U/∂λ的系综平均然后对所有窗口积分就得到完整的自由能差。这个过程叫热力学积分TI而GROMACS里默认推荐的BAR/MBAR分析则基于每个窗口的能量差分布。不管哪种本质都是把巨大变化拆成一系列微小的可逆变化。1.3 为什么结合自由能要算两个“消失过程”单独把配体在蛋白里“变没”得到的自由能变化叫相对结合自由能中的一个分量但它不是最终答案。因为配体在结合前后溶剂环境完全不同要算结合自由能ΔG_bind需要构建一个热力学循环。具体做法是同时跑两个独立的FEP过程第一个过程配体在纯水溶剂中从完整状态“消失”decoupling第二个过程配体在蛋白-配体复合物中从完整状态“消失”把这两个过程的自由能差相减ΔG_bind ΔG_complex - ΔG_water这样做的巧妙之处在于热力学循环把难算的结合过程转换成了两个相对容易收敛的“消除”过程。体系状态函数只取决于始末态所以路径怎么走不影响理论上的结果但会影响实际模拟的收敛速度和误差大小。这个方法常常被称为“double decoupling”或“alchemical binding free energy”是当前基于分子动力学计算结合自由能的主流方案。1.4 什么场景适合用FEPFEP不是万能药。它特别适合骨架相同、取代基不同的一系列类似物做活性排序预测单个位点突变对结合的影响丙氨酸扫描评估几个候选分子与同一靶标的相对结合强弱但它不适合蛋白构象变化剧烈的体系、结合和解离涉及大尺度构象重排的体系、共价抑制剂体系。这类体系用FEP会非常难收敛算出来的自由能可能跟实验值差好几kcal/mol。2. 计算流程总览与方案选型跑之前先把路线定死FEP计算最怕的不是算得慢而是方案没定好就开跑跑完发现路径设计有误所有窗口全部白费。所以动工之前先把整体方案捋清楚。2.1 单拓扑还是双拓扑GROMACS支持两种自由能计算策略单拓扑single topology和双拓扑dual topology。单拓扑的意思是用一套原子坐标、一套成键参数通过键连原子的质量或电荷插值来实现两种状态的转换。它适用于两个状态结构非常相似的情况比如同一个配体上某个取代基从H变成CH3或者从CH3变成CF3。这种方案优势是计算量小、收敛性好因为大多数原子始终存在坐标基本不变。双拓扑则是给两个状态各准备一套完整的原子和参数两套原子在空间上可能重叠通过耦合参数逐步关闭一组、打开另一组。它适用于两个分子结构差异比较大的情况比如完全不同的骨架替换。结合自由能计算中的“配体消失”属于单拓扑里的特殊情况——初态是完整的配体末态是配体与环境完全没有非键相互作用。我一般直接用单拓扑路线也就是mdp里couple-lambda0 vdw-qcouple-lambda1 none配体所有非键相互作用从完整逐渐关闭到零。2.2 慢生长还是分窗采样FEP有两种采样策略慢生长slow growth和分窗采样ladder。慢生长是让λ在一条轨迹里从0持续变到1理论上是绝热过程但实际模拟中λ变化速率必须非常慢否则体系跟不上变化自由能严重偏差。这个方法在早期研究中常用但效率太低现在几乎没人用。分窗采样是目前的主流把[0,1]的λ区间离散成十几个甚至更多窗口每个窗口在固定的λ值下做独立MD模拟最后用BAR/MBAR方法整合所有窗口的能量数据。每个窗口的模拟之间没有严格要求路径一致性只要每个窗口采样充分整体结果就可靠。我习惯用17个窗口λ值分布如下0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.50, 0.60, 0.70, 0.80, 0.85, 0.90, 0.95, 1.00有经验的读者会发现我在0到0.4之间放得密0.4到0.8之间放得疏0.8到1.0又加密。原因很简单λ接近0和接近1的时候能量变化梯度非常大容易出现采样不充分和端点发散所以必须加密窗口中间区域能量变化平缓窗口可以稀疏一些。如果你第一次跑没有把握直接上21个窗口λ从0到1每隔0.05取一个稳妥优先。2.3 力场与参数化的统一FEP结果可信度的上限由力场精度决定。我强烈建议整个体系使用同一个力场家族的参数不要混搭。蛋白用AMBER ff14SB或ff19SB的话配体就用GAFF2或者更好的广义力场参数配体电荷用AM1-BCC计算溶剂用TIP3P水。这些都是经过大量自由能计算验证的组合。具体操作上配体参数化可以用antechamber加acpype。先用antechamber给配体分配GAFF2原子类型跑AM1-BCC电荷然后acpype把prmtop转成GROMACS拓扑和坐标。蛋白拓扑则直接由pdb2gmx生成。特别提醒一个坑配体残基名要统一且唯一。比如把配体残基名设为LIG那么在mdp里couple-moltype LIG蛋白拓扑里不能有其他残基重名。我曾经因为配体残基名叫LIG但体系里恰好有个晶体水也用了类似命名导致GROMACS完全无法定位要耦合的原子报错找了大半天。2.4 体系准备的基本功体系准备阶段虽然听起来跟普通MD一样但有几个自由能计算特有的点需要注意水盒子推荐用dodecahedron十二面体比立方体省体积但四面镜像距离更合理蛋白质加配体后离盒子边缘至少1.0 nm必须加离子中和体系电荷推荐0.15 M NaCl模拟生理盐环境配体和蛋白初始接触姿态要合理最好先用对接或实验构象确定结合模式FEP只能评估给定结合模式下的自由能不能帮你找结合模式能量最小化时配体位置要稳住不要让配体在优化过程中滑出结合口袋3. 核心mdp参数逐行拆解每行都在干什么mdp配置是FEP计算的核心之一网上流传的模板很多但很多人只是一股脑复制出了错不知道改哪。这一节我把关键参数全部拆开讲附带可直接用于生产的完整模板。3.1 必须理解的关键参数先说几个FEP特有的参数这些参数看不懂就贸然开跑基本等于盲飞。free-energy yes开启自由能计算告诉GROMACS在MD过程中额外计算耦合参数相关的能量数据。couple-moltype LIG指定哪一组原子参与“耦合/退耦合”过程。这里填残基名GROMACS会把这个残基的所有原子标记为需要随λ变化的组。couple-lambda0 vdw-q初态的耦合状态表示配体与环境的范德华作用和静电相互作用都完整存在。这里还有一个细节couple-lambda0可以用vdw、vdw-q、none等不同组合表示不同物理状态。couple-lambda1 none末态的耦合状态表示配体与环境的所有非键相互作用全部关闭。配体变成一个“幽灵分子”它存在的意义只是为了维持蛋白-配体复合物的构象不至于因配体消失而崩溃。这里有个容易混淆的地方none指的只是配体与环境之间的非键相互作用关闭配体自身的键长、键角、二面角等成键相互作用以及配体内部原子之间的非键相互作用默认还是存在的。在接受自由能计算时不会出问题因为分子内非键相互作用在结合态和溶液态之间的贡献在设计单拓扑退耦合时已经通过热力学循环抵消掉了。如果你想连分子内相互作用一起关就需要设置couple-intramol yes。不过我要坦白说大多数标准FEP流程直接用默认的no因为热力学循环设计已经让分子内项在相减时消掉了。init-lambda-state 对应的窗口编号GROMACS里λ窗口的编号从0开始每个mdp文件里只跑一个λ窗口。所以你要准备十几份mdp文件每个文件里的init-lambda-state改成对应的序号这也是FEP计算“管理上最烦人”的地方。nstdhdl 20每隔20步写一次dhdl能量数据。这个是分析自由能的关键文件频率太疏会导致数据点不够太密会让文件巨大。20步对应40 fs一个数据点对于300 K下的MD模拟来说足够。separate-dhdl-file yes每个λ窗口的dhdl数据单独存文件避免不同窗口数据混在一起分析时更清晰。dhdl-derivative yes额外输出∂U/∂λ数据。对TI分析有用对BAR分析也有辅助参考价值建议打开。calc-lambda-neighbors -1让GROMACS计算所有相邻λ窗口之间的能量差这样BAR/MBAR分析可以直接使用。如果设成-1就等价于计算了所有λ状态之间的重叠非常方便。3.2 软核势soft-core到底在解决什么问题FEP里最经典的问题出现在λ接近0或1的端点区。想象一下配体的原子和环境原子之间本来有正常的排斥力与吸引力当你逐渐关闭这些相互作用时原子之间可能出现严重的“重叠”或“走近”现象此时如果仍然用正常的LJ势函数计算能量会瞬间爆炸到几百万kJ/mol整个模拟直接崩溃。软核势soft-core就是为了解决这个问题设计的。它把LJ势函数做了数学变换让势能项在相互作用很弱时不再发散而是缓慢趋于有限值。这样即使在端点区粒子之间也可以“穿过”彼此而不产生不可接受的能量尖峰。mdp里相关的参数是sc-alpha 0.5软核强度这个值在GROMACS官方教程和大量文献里都被验证是可靠的。不要随意调大太大了会让势能面变形严重影响结果准确性。sc-sigma 0.3软核作用的特征距离。理论上应该跟体系中典型原子的σ相当GROMACS默认0.3 nm对于大多数有机小分子体系是合适的。如果体系里有特别大的重原子可以适当调大。sc-power 1软核变换的幂次GROMACS推荐1在多数情况下表现稳定。3.3 完整的生产mdp模板下面这套模板是我在当前项目里实际使用的GROMACS版本是2021.x如果是2018或之前版本个别参数名可能有差异但大框架通用。title FEP production ; 运行参数 integrator md dt 0.002 ; 2 fs时间步长 nsteps 5000000 ; 总步数5 ns生产模拟 ; 输出频率 nstxout-compressed 5000 ; 每10 ps存一帧压缩轨迹 nstlog 5000 nstcalcenergy 1 nstenergy 100 ; 键约束 constraints h-bonds constraint-algorithm lincs continuation yes ; 继续前一阶段模拟不做初始化 ; 邻居搜索 cutoff-scheme Verlet nstlist 20 rlist 1.0 ; 静电 coulombtype PME rcoulomb 1.0 fourierspacing 0.12 ; 范德华 vdwtype Cut-off rvdw 1.0 ; 自由能核心区 free-energy yes couple-moltype LIG couple-lambda0 vdw-q couple-lambda1 none couple-intramol no init-lambda-state 0 nstdhdl 20 separate-dhdl-file yes dhdl-derivative yes calc-lambda-neighbors -1 ; 软核 sc-alpha 0.5 sc-sigma 0.3 sc-power 1 sc-coul yes ; 温度耦合 tcoupl v-rescale tc-grps Protein_LIG SOL tau_t 0.1 0.1 ref_t 300 300 ; 压力耦合 pcoupl Parrinello-Rahman pcoupl-type isotropic tau_p 2.0 ref_p 1.0 compressibility 4.5e-5 ; 周期性边界 pbc xyz ; 色散校正 DispCorr EnerPres这套mdp的运行逻辑很简单配体与环境之间的静电和LJ相互作用随λ从完整逐步关闭其他一切维持常规MD。3.4 温度与压力耦合需要注意的细节温度耦合我习惯用v-rescale而不是Berendsen。原因很简单v-rescale满足正则系综的要求能给出正确的涨落Berendsen虽然也常用但它只是“弱耦合”恒温器会抑制涨落算自由能时不如v-rescale严谨。如果模拟时间足够长且体系构象变化不大用no恒温即直接NVT不加恒温器也可以但对大多数体系不推荐。压力耦合在平衡阶段用Berendsen生产阶段我改成Parrinello-Rahman。原因不复杂Parrinello-Rahman在长时间模拟中能正确产生NPT系综的盒涨落但对初始压力的波动敏感如果一上来就用它盒子可能剧烈振荡平衡阶段用Berendsen先让体系稳住体积再切到Parrinello-Rahman是我实测下来最稳的流程。有一点必须注意自由能计算对外界条件比较敏感前后两个FEP腿水相和复合物相的温度、压力耦合参数必须保持一致不然由系综不同引入的系统误差会直接影响最终ΔG_bind。4. 完整实操流程从配体参数化到跑完所有λ窗口理论说得再多不如照着走一遍。下面是我实际操作时的步骤记录每一步都经过验证。4.1 配体参数化与体系搭建第一步给配体做参数。我的流程是从ChemDraw或者PubChem拿到SMILESOpenBabel转成3D结构然后antechamber分配GAFF2力场类型跑AM1-BCC电荷。antechamber -i lig.mol2 -fi mol2 -o lig.am1bcc.mol2 -fo mol2 -c bcc -nc 0 -at gaff2 -rn LIG这里-nc 0表示配体净电荷为0如果你的配体带电要改成实际的净电荷千万别填错。填错电荷会直接导致整个体系的静电参数错误算出来的自由能完全没意义。然后acpype可以把AMBER格式转成GROMACS可用的拓扑acpype -i lig.am1bcc.mol2 -o gromacs -b LIG生成LIG_GMX.gro和LIG_GMX.top。蛋白部分直接用GROMACS自带的pdb2gmx就行力场选AMBER ff14SB。蛋白与配体组合成一个复合物坐标文件后我习惯先用gmx editconf定义盒子gmx editconf -f complex.pdb -o complex_box.gro -d 1.2 -bt dodecahedron然后gmx solvate加水、gmx genion加离子中和。到这里体系准备完成。4.2 能量最小化与平衡能量最小化我通常做两轮先用最速下降法再用共轭梯度法。这一步的目的不是跳出局部极小而是消除体系中可能存在的原子碰撞和过高应力。平衡分两个阶段NVT然后NPT。NVT平衡300 K下跑200 ps把体系温度拉起来NPT平衡再加1 bar压力跑500 ps让盒子体积稳定。注意NVT和NPT平衡阶段都要加上位置约束把配体和蛋白重原子约束住避免溶剂还没平衡好结构就跑飞了。平衡阶段mdp和生产的区别在于free-energy参数必须和保持一致否则整个热力学路径在平衡与生产交接处断裂。很多人只跑普通MD平衡忘了平衡阶段也需要free-energy yes导致生产的时候体系还没有在耦合哈密顿量下达到平衡前几个ns全部作废。4.3 准备并运行所有λ窗口生产模拟前把mdp文件复制17份每个文件里的init-lambda-state改为0到16。然后写一个简单的bash循环脚本批量提交所有窗口任务。for i in 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 do mkdir lambda_${i} cd lambda_${i} gmx grompp -f fep_${i}.mdp -c ../npt.gro -p ../topol.top -o fep.tpr -maxwarn 3 gmx mdrun -deffnm fep -v cd .. done这里有一点要提醒grompp阶段如果报错先检查maxwarn不该随便加。maxwarn3只是跳过警告但如果警告里涉及到自由能参数错误跳过警告等于埋雷。我只有在遇到“原子类型匹配警告”这类无关紧要的提示时才用maxwarn。每个窗口生产模拟跑多长这个没有绝对标准但经验是先跑5 ns试水看所有窗口的∂U/∂λ随时间是否稳定、相邻窗口能量差分布是否重叠好。如果重叠很差加长时间到10 ns甚至20 ns。对于药物分子结合这类体系5 ns每个窗口通常偏短最终结果误差会比较大建议直接上10 ns起。4.4 构建Lambda窗口表和并行运行技巧如果机器核数多另一个高效方式是使用mdrun -multi选项一次提交所有λ窗口但需要准备一个lambda窗口表。GROMACS里这个表是通过-fepw参数指定的内容格式大概是0 0.0000 0.0000 1 0.0500 0.0500 ...或者你直接在每个tpr里设好init-lambda-state用mdrun -multi -lambda 指定文件。两种方法都行我个人更喜欢前者——每个窗口独立目录独立mdp出了问题可以单独重跑某个窗口不需要整批重来。4.5 水相“腿”的准备配体在纯水中的FEP计算本质上和复合物体系一样只是把蛋白去掉只保留一个配体放在水盒子里。配体初始坐标可以从复合物里抠出来也可以重新放到盒子中心但一定要保证配体不被自己的周期镜像干扰盒子里配体到边缘的距离至少1.5 nm。复合物相和水相的模拟参数、窗口划分、λ分布应当完全一致这是自由能相减能够相互抵消系统误差的前提。我两次跑完对比过如果一边用17个窗口一边用13个窗口即使总自由能结果接近误差也会明显增大。5. 数据分析从dhdl.xvg到最终的ΔG_bind跑完模拟只是第一步真正麻烦的是数据分析这一步。很多新手跑完几十个ns的模拟却不知道该怎么从一堆dhdl.xvg文件里拿到最终答案我详细讲一下。5.1 用gmx bar估算ΔGGROMACS自带的gmx bar命令可以直接读取多个窗口的dhdl.xvg文件然后用BAR方法估算相邻窗口之间的自由能差最后加和成总的ΔG同时给误差估计。标准命令如下gmx bar -f lambda_0/fep_dhdl.xvg lambda_1/fep_dhdl.xvg ... lambda_16/fep_dhdl.xvg -o bar.xvg -g bar.log需要注意的是dhdl.xvg文件里记录了所有λ状态的能量差得益于calc-lambda-neighbors -1所以gmx bar可以一次性整合所有窗口的数据不需要你手动算两两窗口的ΔG。输出结果里会有一张表列出每对相邻λ窗口的ΔG和误差最后还有总的ΔG和累计误差。我个人习惯把每一次生产的轨迹分成两半分别算前半段和后半段的ΔG如果两条半段的自由能差超过0.5 kcal/mol就说明模拟时间不够长、采样没有收敛结果不可采信。5.2 MBAR与更精细的分析gmx bar本质上只是BAR的重复应用更严谨的做法是MBARMultistate Bennett Acceptance Ratio它能同时利用所有窗口的所有数据进行全局优化估计理论误差更小。GROMACS新版可以输出dhdl数据后用pymbar或alchemical-analysis Python脚本做MBAR分析。alchemical-analysis工具是我的首选用起来也简单alchemical-analysis -d lambda_*/fep_dhdl.xvg -u kcal -m BAR --overlap它会自动生成收敛性分析图、重叠矩阵图还能画自由能累积曲线非常直观。唯一的麻烦是要装Python环境和依赖但现在的conda已经能一键装好不算大问题。5.3 计算ΔG_bind以复合物相和水相分别跑完FEP得到ΔG_complex和ΔG_water最终ΔG_bind ΔG_complex - ΔG_water从热力学循环也可以校验符号的合理性配体结合得越强结合自由能负值越大。如果算出来的结果是正的先把数据复查一遍检查有没有窗口跑飞、有没有初始结构不合理。另外要说明的是FEP的一次运行结果是一个样本严格的做法是跑2到3次独立重复初始速度不同取平均并计算标准误。这样最终误差才能反映统计不确定性而不是单次模拟的偶然波动。5.4 收敛性判断三板斧判断FEP计算是否收敛我一般看三个指标第一是每个λ窗口的∂U/∂λ时间序列是否自洽。把每个窗口的dhdl数据按时间画折线如果曲线在模拟后期是一条围绕平均值稳定波动的水平线说明该窗口采样充分如果持续漂移说明没有收敛。第二是相邻窗口的能量差分布重叠。BAR方法要求相邻窗口的能量差分布有明显重叠区域分布完全分离意味着两个窗口之间缺少协同采样结果不可靠。alchemical-analysis里的重叠矩阵图就是干这个的重叠比例最好都大于0.03。第三是正反路径一致性。理论上从λ0往λ1跑和从λ1往λ0跑得到的ΔG应该严格相同。实际中如果两条路径相差超过1 kcal/mol说明存在明显的采样滞后模拟时间要加倍。6. 常见问题速查表与避坑技巧实录这部分内容是我踩坑之后总结的。FEP计算的坑非常多我把最常见的问题列成表方便随时翻阅。问题表现可能原因排查与解决mdrun报错“No atom (LIG) in topology”couple-moltype里的残基名与拓扑中不一致检查拓扑文件中配体残基名确保与mdp一致λ接近0或1时能量暴涨软核参数未开启或设置不当确认sc-alpha0.5、sc-power1检查软核段是否有被注释dhdl.xvg文件为空或没有数据nstdhdl未设置或free-energyyes没有打开检查mdp中的free-energy、nstdhdl、separate-dhdl-file所有窗口的ΔG波动巨大每个窗口采样不足或相邻λ窗口间隔过大增加窗口密度延长模拟时间检查状态重叠复合物相配体从口袋中跑出来初始结构未做约束平衡或结合模式不合理重新准备初始结构NVT/NPT平衡阶段对配体加位置约束前后两条FEP腿结果差异极大两条腿的参数不一致或盒子大小差异太大统一两边mdp所有参数确保盒子尺寸和离子浓度一致gmx bar输出总误差超过1 kcal/mol采样时间不足或某些窗口分布重叠差延长模拟加密不良窗口的λ分布MD过程中蛋白结构严重漂移力场不兼容或约束参数错误检查力场搭配确认约束组设置正确水相FEP中配体自相互作用干扰盒子太小配体与周期镜像接触增大盒子到配体边缘至少1.5 nm用十二面体盒子我额外补充几个容易在“看不太出来”的地方出问题的经验。第一个是关于配体的质子化状态。配体在生理pH下的质子化状态一定要在参数化之前确认好不然带电基团的质子化状态错了静电相互作用整体偏移算出来的自由能跟实验值对不上。最简单的方式是用计算pKa的工具先预测一下或者在文献里查化合物在不同pH下的解离常数。第二个是二面角参数。GAFF2对普通有机分子一般表现良好但遇到杂环、卤素、含有未共用电子对的氮原子时二面角参数有时候会不太理想。如果做完FEP发现某个系列分子的排名和实验趋势有系统性偏差可以优先怀疑配体的二面角参数而不是自由能方法本身。解决办法是重新参数化二面角或者换用像CGenFF这样专门针对小分子的参数集。第三个是体系里是否含有金属离子。许多靶标蛋白含有金属离子比如锌指蛋白、金属酶这类体系的FEP计算必须格外小心因为金属离子的配位作用对配体结合影响极大但AMBER的默认参数对金属离子的处理常常不够精确需要额外使用阳离子虚拟模型cation dummy model或者专门的金属中心参数。新手第一次做FEP千万别选金属蛋白体系练手会被虐得很惨。第四个是水分子的处理。结合口袋里如果有结晶水或者结构水它们对结合自由能的影响可能很大。为此有的严格流程会把这些水分子保留在体系中并让它们也参与自由能计算。但这样会显著增加计算复杂性新手阶段先不必上这个强度只要注意配体周围不要有异常靠近的水分子就行。7. FEP计算的成本控制与算力规划最后一个实用话题FEP很贵怎么规划算力才不浪费。做一个完整的结合自由能计算包括两条腿、17个窗口、每个窗口10 ns生产总模拟时间是340 ns。复合物体系如果有3万到5万个原子用一张主流GPU卡跑每天大约能完成100到200 ns取决于GPU型号也就是说整体跑完需要2到4天。如果做17个分子的一系列排名就要做好计算量乘以分子数的心理准备。成本控制的核心思路是先用短时间初筛再用长时间精算。第一次跑的时候每个窗口跑2到3 ns快速扫一遍所有分子的能量学趋势。看趋势是否合理不合理就直接淘汰或修正不要浪费时间在错误体系上精算。等初筛通过后再把最终要用的分子加上长时间模拟每个窗口至少10 ns条件允许跑到20 ns保证最终结果的统计精度。另外一定要并行利用多个GPU。GROMACS的mdrun能够同时管理多个窗口如果机器上有8张卡开8个mdrun进程每个进程处理几个窗口整体效率会显著提升。如果只有一台单卡机器那就老老实实排队跑但可以通过缩短初筛时间来控制总成本。还有一点建议跑之前先把所有需要分析的中间文件都设置好别在跑完之后才发现轨迹没输出或者dhdl数据没保存。跑完一次FEP再重新跑不管是时间还是电费都是很大的浪费。8. 结语FEP是一条值得走的路但要带着脑子走做FEP跟做普通MD最大的区别在于普通MD跑出轨迹看一眼构象跑完就完事FEP整个流程里每一个环节的失误都会累积成最终的误差而且这个误差往往是隐藏的不像能量爆炸、原子穿模那样肉眼可见。所以在跑FEP之前先花时间想清楚体系、力场、耦合路径、采样长度这些关键决策比急着把窗口跑出来重要得多。我个人这几年做FEP最深的体会是别指望第一次跑就能拿到和实验完美吻合的数字更别指望FEP能自动给你答案。FEP的意义在于当你有一系列结构高度相关的分子时它能给出比打分函数可靠得多、比实验便宜得多的相对排名而这个排名结果足以支撑你先导化合物优化的下一步决策。如果你正在为某个药物设计项目做FEP计算我的建议是先用一两个已知活性数据的分子做验证把整条流程跑通并确认误差可控然后再批量扩展。这样上算力之前你至少知道自己的计算方案在这套体系里是能出活儿的。祝每位做自由能计算的同行都能少踩几个坑多出几个准得让实验组惊讶的漂亮数据。