MM/GBSA：快速自由能估算的优势和局限

MM/GBSA 用广义玻恩（GB）近似替代 MM/PBSA 中的泊松-玻尔兹曼方程求解。速度快 10～100 倍，精度略低——这个权衡让它成为实践中用得最多的自由能估算方法。

GB 与 PB 的区别

# PB（泊松-玻尔兹曼）：
#   把溶剂当作连续介质，
#   数值求解泊松-玻尔兹曼方程得到静电溶剂化能
#   → 需要在三维格点上迭代求解
#   → 精确但慢
#
# GB（广义玻恩）：
#   用解析公式近似溶剂化能：
#
#     ΔG_pol ≈ -½ (1 - 1/ε) Σ_ij  q_i q_j / f_GB(r_ij, α_i, α_j)
#
#   其中 α_i 是原子 i 的「有效玻恩半径」——
#   表示该原子被埋藏的程度
#
#   核心思想：
#     深埋的原子（大 α）溶剂化能小
#     暴露的原子（小 α）溶剂化能大
#
#   → 只需算原子间距离与有效半径，无需解方程
#   → 快得多
#
# 精度差异的主要来源：
#   有效玻恩半径的计算是近似的
#   → 对带电基团、深埋原子的误差较大
#   → 【净电荷不同的分子之间比较时误差最明显】

GB 模型的选择

igb	模型	特点
1	HCT	最早的版本，现在少用
2	OBC1	常用，平衡较好
5	OBC2	最常用的默认选择
7	GBn	改进的有效半径计算
8	GBn2	较新，对某些体系更好

不同 GB 模型给出的绝对值可能相差数 kcal/mol——因此比较不同来源的 MM/GBSA 结果时必须确认用了相同的 igb。实践建议：在自家体系上用已知活性数据测试 igb=2、5、8，选排序能力最好的，然后固定不变。

用 OpenMM + 简化流程计算

# 除了 gmx_MMPBSA，也可以用 OpenMM 直接算
# 思路：对轨迹的每一帧，分别算复合物/蛋白/配体的能量

from openmm import app, unit
import openmm as mm
import mdtraj as md
import numpy as np

def gbsa_energy(topology, system, positions):
    """单点能计算（含 GB 溶剂化）"""
    integrator = mm.VerletIntegrator(1.0 * unit.femtosecond)
    context = mm.Context(system, integrator,
                         mm.Platform.getPlatformByName("CPU"))
    context.setPositions(positions)
    state = context.getState(getEnergy=True)
    e = state.getPotentialEnergy().value_in_unit(unit.kilocalorie_per_mole)
    del context, integrator
    return e

# 完整流程（示意）：
#   1) 用 implicitSolvent=app.OBC2 创建三个 system：
#      复合物、蛋白、配体
#   2) 遍历轨迹帧，各算一次能量
#   3) ΔG = 平均(E_复合物) - 平均(E_蛋白) - 平均(E_配体)
#
# 实际项目建议直接用成熟工具：
#   gmx_MMPBSA   （功能最全，支持多种 MD 软件）
#   MMPBSA.py    （AmberTools 自带）
#   OpenMM 脚本  （需要自己实现，适合定制需求）

# 快速版本：不跑 MD，只用对接姿势做单点
#   → 「最小化 + 单点 GBSA」
#   → 速度接近对接，精度略优于对接分数
#   → 适合作为对接后的第一道重打分
#   → 但可靠性明显低于基于 MD 系综的计算

典型流程与时间预算

步骤	时间（单卡 GPU，中等体系）
体系准备 + 参数化	数十分钟（人力为主）
最小化 + 平衡	约 1 小时
生产 MD（20 ns）	约 2～4 小时
MM/GBSA 计算（500 帧）	数分钟
（对比）MM/PBSA	数小时
（对比）FEP 单对	数小时～十几小时

瓶颈在 MD 而非 GBSA 计算本身。因此如果要评估很多分子，主要成本是每个分子都要跑 MD。

什么时候用 MM/GBSA

• 对接后的重打分：把对接的前几百个分子缩到几十个（见 095《Docking Score》 的分级策略）；
• 同系列类似物的排序：这是最合适的场景——骨架相同、电荷相同、大小相近；
• 热点残基识别：逐残基分解（见 125《MM/PBSA》）；
• 不适合：
  • 跨骨架比较；
  • 净电荷不同的分子之间比较——GB 对电荷特别敏感；
  • 需要定量的绝对亲和力；
  • 差异很小的类似物（需要 FEP 的精度）。

常见错误

# 错误一：把结果当作绝对结合自由能
#   MM/GBSA 的绝对值通常严重高估结合强度
#   （数 kcal/mol 到十几 kcal/mol 的偏差很常见）
#   → 只用相对值

# 错误二：不跳过平衡阶段
#   前 10~20 ns 的能量还在漂移
#   → startframe 必须设在平衡之后
#   → 判断方法：看能量随时间的曲线是否平稳

# 错误三：只跑一个副本就下结论
#   MM/GBSA 结果的涨落可能很大
#   → 至少 3 个独立副本
#   → 报告均值 ± 标准误

# 错误四：比较净电荷不同的分子
#   如比较中性分子与带负电的羧酸类似物
#   → 溶剂化能的误差会主导结果
#   → 这类比较应该用 FEP

# 错误五：忽略配体的构象重组能
#   单轨迹方法假设配体结合前后构象不变
#   → 对柔性配体，重组能可能很大
#   → 可用三轨迹方法（分别模拟复合物/蛋白/配体）
#     但误差不再抵消，结果反而可能更差

# 错误六：混用不同的 igb 值比较
#   不同 GB 模型的绝对值不可比
#   → 一个项目内固定 igb

# 正确的报告方式：
#   「用 MM/GBSA（igb=5，3 个独立 20 ns 副本，
#     取后 10 ns，每 20 ps 一帧，未含熵项）
#     得到的相对排序为：A < B < C，
#     其中 A 与 B 的差异在误差范围内」

在方法谱系中的位置

# 精度递增，成本递增：
#
#   对接打分       毫秒   相关性低
#      ↓
#   CNN 重打分     秒     略好
#      ↓
#   【MM/GBSA】    小时   中等 ← 性价比合适的中间层
#      ↓
#   MM/PBSA        小时+  略好
#      ↓
#   Boltz-2        分钟   中到高（见 121）
#      ↓
#   FEP/RBFE       小时++ 高（见 127）
#      ↓
#   实验测定       天~周  金标准
#
# 实际项目中的典型用法：
#   对接 10 万 → MM/GBSA 前 500 → FEP 前 30 → 合成 10
#
# 注意：
#   Boltz-2 的出现让这个谱系有了新选择 ——
#   速度接近 MM/GBSA 而准确度可能更好
#   → 但需要在自家体系上验证后再决定

延伸资源

• MM/PBSA：125《MM/PBSA》；FEP：127《FEP、RBFE 与 ABFE》；MD 入门：123《分子动力学 MD 入门》；
• OpenMM：128《OpenMM 跑 10–100 ns MD》；Boltz-2：121《Boltz-2 技术报告精读》。