GPCR靶点分子动力学与自由能计算 - 第2节

2.1 经典力学到底在算什么？

做自由能计算，你得先理解分子模拟在干嘛。说白了，就是用牛顿第二定律去推原子怎么动。每个原子带一个质量，受力了就会加速。你在跑 MD 的时候，每步都在解：

F = ma
a = F / m

这个 F 怎么来的？不是量子力学里的电子云，而是靠力场（Force Field）算出来的。力场就像一套「经验公式」，告诉你原子间怎么推怎么拉。

我刚开始做 GPCR 模拟那会儿，总觉得力场太粗糙。直到有一次对比 QM 和 FF 的势能面，发现对大多数构象，力场的误差控制在 1 kcal/mol 以内。嗯，足够了。

2.1.1 势能函数长什么样？

常见力场的总势能可以拆成这几块：

E_total = E_bond + E_angle + E_dihedral + E_vdw + E_elec

E_bond：两原子成键，拉伸谐振
E_angle：三个原子形成的键角弯曲
E_dihedral：四个原子形成的二面角扭曲
E_vdw：范德华作用（Lennard-Jones 12-6）
E_elec：库仑静电作用

我个人习惯：拿到一个新力场，先看一眼它的二面角项参数矩阵。那里面藏着很多 torsion barrier 的信息，直接影响自由能差的计算。

2.2 力场到底是什么？

你想想看，力场就是一台「标尺」。它把分子中每个原子的类型、电荷、平衡键长、力常数都定好了。你给出一组坐标，它立刻告诉你体系的总能量和每个原子上的力。

我记得有一次调试一个配体参数，发现键长力常数给得太软，导致模拟中键疯狂伸缩，能量爆炸。后来改成 AMBER 推荐值，才稳住。

2.2.1 力场的组成模块

模块	描述	典型参数
原子类型	区分不同化学环境（如 sp3 碳 vs 芳香碳）	CT, CA, N3, O2
键参数	K_r, r_eq	AMBER: 310 kcal/(mol·Å²)
角参数	K_θ, θ_eq	CHARMM: 60 kcal/(mol·rad²)
二面角参数	V_n, γ, n	OPLS: 0.15 kcal/mol
非键参数	ε, σ（vdW）; q（电荷）	AMBER: ε=0.1094, σ=3.4

2.3 主流力场：AMBER、CHARMM、OPLS

做 GPCR 自由能计算，碰得最多的就是 AMBER 和 CHARMM。两个都是生物分子力场里的老大哥，但脾气不太一样。

2.3.1 AMBER 力场

AMBER 的设计哲学是「参数化尽量拟合小分子气相的 QM 数据」。它的非键参数尤其注重静电匹配。

适用体系：蛋白质、核酸、小分子配体（配合 GAFF/GAFF2）
典型版本：ff14SB, ff19SB, Lipid21
特点：对蛋白主链二面角描述精准，适合长时间尺度的 GPCR 模拟

我在项目里常用 ff14SB 加 GAFF2 处理配体。你注意一下，AMBER 的 1-4 静电和 vdW 项会有缩放因子，默认是 1/1.2。这个细节容易被人忽略。

# AMBER 力场参数文件片段 (prmtop 里实际存的是二进制，但简要格式如下)
%FLAG POINTERS
%FORMAT(10I8)
    ...
%FLAG BOND_FORCE_CONSTANT
%FORMAT(5E16.8)
  310.00000000E+00 ...
%FLAG BOND_EQUIL_VALUE
  1.52500000E+00 ...

2.3.2 CHARMM 力场

CHARMM 更强调对凝聚相性质的拟合，尤其是液体密度和水合自由能。常见的有 CHARMM36, C36m。

适用体系：蛋白质、膜、糖类、离子
典型版本：CHARMM36m（对结构柔性的 GPCR 优化）
特点：对脂质膜描述很细致，膜蛋白模拟里的首选

我曾经踩过的坑：CHARMM 力场的 CMAP 修正项专为蛋白主链设计。如果你做 GPCR 的 loop 区采样，不要自作主张关掉 CMAP，否则二面角分布会偏。

2.3.3 OPLS 及其他

OPLS（Optimized Potentials for Liquid Simulations）偏重有机小分子和结合自由能预测。它的 vdW 参数往往比 AMBER/CHARMM 更「软」一些。

力场	典型应用	我推荐的场景
AMBER ff14SB + GAFF2	蛋白-配体结合自由能	GPCR 胞外域配体结合
CHARMM36m	膜蛋白模拟	GPCR 跨膜区结构稳定
OPLS-AA/M	小分子溶剂化自由能	配体相对结合能力快速筛选

2.4 一张图看懂力场能量项

2.5 自由能计算中如何选力场？

GPCR 靶点自由能计算最怕力场选错。我遇到过不少人，用 AMBER 跑膜蛋白，结果脂质参数不兼容，模拟跑着跑着膜就破了。

我个人建议：

若你的 GPCR 在膜环境中且需要长时间采样 → CHARMM36m（脂质参数最成熟）
若你主要关注胞外区配体结合且配体种类多 → AMBER ff19SB + GAFF2（小分子参数覆盖广）
若你想快速扫一遍相对结合自由能 → OPLS-AA/M (计算快，收敛好)

2.5.1 一个典型的工作流

以 AMBER 体系为例，做 GPCR-配体自由能计算的大致步骤：

用 tleap 生成蛋白拓扑（ff19SB）
用 GAFF2 生成配体力场参数（antechamber + parmchk2）
嵌入膜（用 CHARMM-GUI 或 Membrane Builder 给 AMBER 转格式）
平衡模拟（NVT → NPT）
自由能计算（TI / FEP / MBAR）

# 用 antechamber 生成配体参数（示例）
antechamber -i ligand.mol2 -fi mol2 -o ligand.prepc -fo prepc -c bcc -at gaff2
parmchk2 -i ligand.prepc -f prepc -o ligand.frcmod -a gaff2

曾经有个项目，我用 CHARMM36 跑 GPCR，配体参数用 CGenFF 生成后直接拿来用。结果发现二面角势垒参数是缺失的，系统提示自动补了默认值。跑出来的自由能差了 2.5 kcal/mol。后来手工补了 QM 扫描的参数才修正。

2.6 小结

力场是自由能计算的基石。你花点时间理解各能量项的物理含义，选对力场版本，配好脂质模型，后面就顺了。别指望一个力场打天下——AMBER 和 CHARMM 各有所长，按场景来选就行。