结合自由能计算在药物设计、蛋白-配体相互作用和主客体化学中承担着从虚拟筛选到先导化合物优化的关键使命。但在遇到一个具体体系时,第一个需要回答的问题不是”怎么算”而是”用哪种方法算”——MM/PBSA、自由能微扰(FEP)和热力学积分(TI),三者的精度跨度和计算成本差出两个数量级。

MM/PBSA的理论基础是把结合自由能拆成气相分子力学项、极性溶剂化能(PB)和非极性溶剂化能(SA)三部分,分别用MD采样的构象系综做后处理平均。不需要在结合/解离态之间做显式的路径采样,这是它跟FEP/TI最大的不同——也因此计算量低了不止一个数量级。
一个MD轨迹跑完蛋白-配体复合物(通常20-50 ns),抽出200-500帧等间距snapshot,每帧分别算复合物、蛋白、配体的MM能量+PB极性项+SA非极性项,差值平均即为结合自由能。整个计算机时约等于一次MD的长度——几十到几百核时。
但MM/PBSA有结构性的精度上限。隐式水模型对特定水分子介导的氢键网络完全盲视(蛋白-配体结合界面若存在有序水分子,PB的连续介质近似直接失效),而且不包含构象熵的完整贡献。正常情况下MM/PBSA的误差在2-5 kcal/mol范围——对区分nM和μM级别的亲和力不够,但对μM到mM级别的rank-ordering基本胜任。
MM/PBSA的正确使用场景是:已有一批对接pose需要快速排序;或者蛋白-配体界面水分子较少(疏水结合占主导)。不适合的场景:金属离子介导的配位结合、核酸结合(静电主导,PB的线性近似失效)、以及追求<1 kcal/mol级别精度的亲和力预测。
自由能微扰的核心哲学是”在两个化学态之间炼金术般地逐渐变化”——不是直接算A→B的结合自由能差,而是在参数λ空间沿着一条渐变的路径采样。每个λ窗口内,Hamiltonian介于A和B之间:H(λ) = (1-λ)H_A + λ·H_B,体系在每个λ窗口独立平衡,相邻窗口之间计算自由能差。
FEP的精度取决于λ窗口的密度和每个窗口的采样时长。对典型的蛋白-配体结合自由能,需要12-24个λ窗口(每个窗口2-5 ns),总模拟时长可达100-500 ns——计算成本是MM/PBSA的50-100倍。
FEP能做什么MM/PBSA做不到的?区分—CH₃和—CF₃取代基对结合自由能的影响。这类微小化学修饰的结合自由能差通常在0.5-2 kcal/mol级别——MM/PBSA的噪声淹没了信号,FEP的信噪比可以让这个差异从噪声中浮现出来。
但FEP的成败高度依赖三个前提条件:扰动路径不走环(即λ路径上不出现不可逆的非平衡过程)、力场对两个末端态的构象空间都有足够代表性、以及采样循环足够。尤其是扰动基团较大时(苯环→萘环),构象空间急剧膨胀,每个λ窗口的采样需求直接翻倍。
热力学积分与FEP共享同样的炼金术框架,但用积分形式代替差分:ΔG = ∫₀¹ ⟨∂H(λ)/∂λ⟩_λ dλ——沿λ路径对Hamiltonian导数做系综平均再积分。理论上TI比FEP更数序化,因为在连续极限下ΔG的表达式是否定的。
实际经验中,TI在积分端点(λ=0和λ=1)附近容易出问题——∂H/∂λ在这两个端点可能出现奇点或大梯度,数值积分对端点行为极敏感。解决方案是在端点加密λ窗口(端点的λ间距从0.05缩小到0.02),但这增加了计算成本。
TI的核心优势在于一次模拟可以同时评估多个λ窗口——计算哈密顿量导数的成本远低于FEP的显式能量差。但当扰动涉及电荷变化时(配体从电中性变为带电),周期性边界条件的人工镜像电荷效应必须用适当的静电校正方案处理,否则TI积出来的自由能会有系统性偏移。
一个可操作的经验规则:如果实验结合自由能在-5到-10 kcal/mol范围且只需要rank-ordering(排序而非定量预测),MM/PBSA足够;如果需要预测<1 kcal/mol的修饰效应,FEP是必选。
在这个计算集群资源日益廉价、但方法选择仍需要物理直觉的时代,结合自由能计算的真正门槛不在算力的多寡,而在于对方法学近似层次的理解——知道MM/PBSA的2 kcal/mol误差从哪来,比盲目上FEP更有价值。
更多内容请访问 https://www.keyanxueshu.com/
吸附自由能计算:从吸附能到表面反应自由能的完整路径
吉布斯自由能计算:溶剂化自由能与热力学循环
化学反应吉布斯自由能计算:从OK电子能到反应ΔG的修正链路
结合自由能计算:MM/PBSA、FEP与TI的方法选择
CASTEP计算形成能:Materials Studio环境下的缺陷热力学建模
迁移能垒理论计算:NEB/CI-NEB的初猜策略与鞍点定位
能带模拟计算:三大主流软件的赝势路线与收敛策略
能带结构计算:色散关系、带隙类型与有效质量的正确打开方式
Gaussian计算结合能:BSSE校正与超分子方法的精度博弈
Gaussian计算:基组、相关能与溶剂模型构成的参数决策链
高斯算活化能:过渡态搜索、IRC验证与热力学校正的完整实践
二维材料DFT计算中的范德华修正与层数依赖性问题