企业官网建设流程全解析

做药物设计或者先导化合物优化的人，迟早都会碰到同一个问题：手里有十几个看起来都不错的小分子，到底哪个跟靶标结合得更牢？实验上可以跑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 什么场景适合用FEP

FEP不是万能药。它特别适合：

• 骨架相同、取代基不同的一系列类似物做活性排序
• 预测单个位点突变对结合的影响（丙氨酸扫描）
• 评估几个候选分子与同一靶标的相对结合强弱

但它不适合：蛋白构象变化剧烈的体系、结合和解离涉及大尺度构象重排的体系、共价抑制剂体系。这类体系用FEP会非常难收敛，算出来的自由能可能跟实验值差好几kcal/mol。

2. 计算流程总览与方案选型：跑之前先把路线定死

FEP计算最怕的不是算得慢，而是方案没定好就开跑，跑完发现路径设计有误，所有窗口全部白费。所以动工之前，先把整体方案捋清楚。

2.1 单拓扑还是双拓扑

GROMACS支持两种自由能计算策略：单拓扑（single topology）和双拓扑（dual topology）。

单拓扑的意思是用一套原子坐标、一套成键参数，通过键连原子的质量或电荷插值来实现两种状态的转换。它适用于两个状态结构非常相似的情况，比如同一个配体上某个取代基从H变成CH3，或者从CH3变成CF3。这种方案优势是计算量小、收敛性好，因为大多数原子始终存在，坐标基本不变。

双拓扑则是给两个状态各准备一套完整的原子和参数，两套原子在空间上可能重叠，通过耦合参数逐步关闭一组、打开另一组。它适用于两个分子结构差异比较大的情况，比如完全不同的骨架替换。

结合自由能计算中的“配体消失”属于单拓扑里的特殊情况——初态是完整的配体，末态是配体与环境完全没有非键相互作用。我一般直接用单拓扑路线，也就是mdp里couple-lambda0 = vdw-q，couple-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拿到SMILES，OpenBabel转成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不该随便加。maxwarn=3只是跳过警告，但如果警告里涉及到自由能参数错误，跳过警告等于埋雷。我只有在遇到“原子类型匹配警告”这类无关紧要的提示时才用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估算ΔG

GROMACS自带的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的重复应用，更严谨的做法是MBAR（Multistate 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计算的坑非常多，我把最常见的问题列成表，方便随时翻阅。

我额外补充几个容易在“看不太出来”的地方出问题的经验。

第一个是关于配体的质子化状态。配体在生理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计算，我的建议是：先用一两个已知活性数据的分子做验证，把整条流程跑通并确认误差可控，然后再批量扩展。这样上算力之前，你至少知道自己的计算方案在这套体系里是能出活儿的。祝每位做自由能计算的同行都能少踩几个坑，多出几个准得让实验组惊讶的漂亮数据。