2. 自由能计算理论基础归纳总结
做 MM/GBSA 之前，我建议你先搞明白自由能到底在算什么。说白了，自由能就是体系在给定条件下能做多少“有用功”。
我在项目中遇到过不少同学，上来就对着 Amber 脚本一通跑，结果算出来的结合自由能正负号都不对。嗯，根本原因就是理论基础没打牢。这一节我就带你过一遍核心概念，不搞花架子，全是干活用的东西。
2.1 自由能是什么？
热力学第二定律告诉我们，系统自发过程的方向是熵增。但实际生物分子体系往往在恒温恒压下，这时候判断自发用的就是吉布斯自由能 G。
你想想看，两个分子结合，如果 ΔG < 0，那结合就是自发的；如果 ΔG 很大负值，结合就很强。反过来，ΔG > 0 则结合困难。
对于结合过程：ΔG_bind = G_complex - G_protein - G_ligand
我刚开始学的时候总把焓和熵搞混。后来自己搭了个小模型走了一遍，才明白：焓变反映键能、非键作用的变化，熵变反映自由度变化。两个都重要，但 MM/GBSA 里熵常常被简化处理——后面我会专门说这个坑。
2.2 MM/GBSA的数学骨架
MM/GBSA 本质上是一种端到端的能量分解方法。它把结合自由能拆成几个物理上可解释的部分：
ΔG_bind = ΔE_MM + ΔG_GB + ΔG_SA - TΔS这里每个项都有物理含义：
| 能量项 | 含义 | 我踩过的坑 |
|---|---|---|
| ΔE_MM | 分子力学气相相互作用能（静电 + vdW + 内能） | 内能变化在刚性受体时可以直接忽略，但柔性对接时不能偷懒 |
| ΔG_GB | 广义Born溶剂化能（极性溶剂贡献） | Born半径的选择很敏感，不同力场参数换一套，结果能差 2-3 kcal/mol |
| ΔG_SA | 非极性溶剂化能（空穴 + 疏水作用） | 通常用 SASA 乘以经验系数，但这个系数0.005~0.007随便调，你得对着实验标定 |
| -TΔS | 熵贡献（平动+转动+振动） | 计算成本高，我很多时候直接设为零，只做相对排序 |
注意：ΔE_MM 的计算需要从轨迹中提取平均结构，或者对多帧取平均——我个人习惯用单轨迹法（single trajectory），因为受体和配体构象变化完全耦合，误差抵消好。
2.3 溶剂化模型：GB到底怎么算？
GB 模型将溶剂视为连续电介质，用 Born 公式近似描述原子在溶剂中的极化能。核心是每个原子的 Born 半径：
ΔG_pol ≈ - (1/2) * (1/ε_in - 1/ε_out) * Σ_i q_i^2 / R_i这里 R_i 就是有效 Born 半径，它受周围原子的屏蔽影响。我曾经碰到过一个坑：同一个配体，用 GBOBC 和 GBNeck 参数跑出来的结果，排名完全反过来。后来查文献才发现，OBC 对大疏水腔体 overestimate，Neck 更准但计算慢。
2.4 熵的陷阱
熵项 (-TΔS) 在 MM/GBSA 里是个烫手山芋。理论上要用正则模态分析或准谐波模型计算，但你想一下：一个复杂蛋白配体复合物有几千个原子，算一次振动熵得跑对角化——我试过，48 核跑了两天才出结果。
而且实际项目中，我们更关心不同配体之间的相对结合强弱。如果所有配体的熵变差不多，那忽略 -TΔS 对排序几乎没影响。我总结了一个经验规则：如果配体是刚性类似物，放心忽略熵；如果配体差别很大（比如链状 vs 环状），那必须算熵。
2.5 一张图理清思路
下面这张 SVG 展示了 MM/GBSA 自由能计算的标准流程。我从轨迹处理到最终打分，每一步的关键决策都标出来了。你以后写论文或者做报告，直接拿这个图改改就能用。
2.6 实战中的常见坑
最后我把这些年碰到的坑总结一下，你直接对照自查：
- 轨迹平衡判断——别用前几纳秒的轨迹，RMSD 没稳定之前的数据全是噪声。我习惯取最后 5 ns 的 100 帧。
- 截断半径设置——LJ 和静电的 cut 至少 12 Å，GB 模型则对长程效应更敏感。我常用 20 Å 的 cut 保证收敛。
- 熵的计算代价——如果非要加熵，只对非简谐振动部分用正则模态，平动转动直接解析公式就够了。
- 结果单位混淆——Amber 默认单位是 kcal/mol，很多新手拿着 kJ/mol 的数据对比文献，直接差 4 倍。注意换算。
- 忽略离子效应——如果你体系有带电残基且溶剂中有盐，GB 模型一定要设 istrng > 0，否则静电贡献会被严重高估。
嗯，理论部分就这些。记住一条核心：MM/GBSA 不是第一性原理，它的价值在于相对排序而不是绝对数值。你只要对标实验数据调整参数，就能得到靠谱的结论。