第3章 电荷计算：RESP与HF/6-31G*方法的参数设置、电荷约束的陷阱、多重构象电荷平均的注意事项
电荷分配，说白了就是决定分子里每个原子带多少电。力场里90%的静电相互作用都靠它，一旦弄错，后面跑出的结果基本没法看。我刚开始接触力场参数化时，觉得电荷嘛，随便找个半经验算算就行……后来被狠狠坑过一次——模拟的蛋白-配体结合自由能差了8 kcal/mol，调了俩月才发现是电荷的问题。嗯，从那以后我再也不敢在电荷上偷懒了。
这一章，我们专门聊RESP电荷配合HF/6-31G*的计算套路。不讲理论推导，只讲怎么设参数、怎么避开常见坑、怎么处理多重构象的平均。
3.1 RESP方法与HF/6-31G*基组——为什么是黄金组合？
RESP（Restrained Electrostatic Potential）电荷，本质上是通过拟合QM算出来的静电势，得到一套与静电势匹配最好的点电荷。它自带一个抛物线约束，防止电荷跑偏。而HF/6-31G*这个基组，我个人的理解是它“故意”高估了一点偶极矩（大约高10%~20%），这样在力场里用的时候反而能补偿无极化水的环境——说白了，它就是为凝聚相模拟优化的。
所以，如果你做的是生物分子在小分子配体或者蛋白-配体体系的力场模拟，RESP + HF/6-31G* 几乎是业界的默认选择。
实际操作时，大部分人用Gaussian做ESP计算，然后用Antechamber或者手动脚本跑RESP。我习惯在Gaussian里加入以下关键词：
#P HF/6-31G* Opt POP=CHELPG SCF=Tight IOp(6/33=2)解释一下：
- POP=CHELPG：输出基于Merz-Kollman格点采样的ESP数据（默认是CHelpG，但CHELPG更稳）
- IOp(6/33=2)：将格点密度设为2（每Bohr²约2个点），够用
- SCF=Tight：收敛阈值高一点，减少静电势的噪声
避坑：我见过有人用B3LYP/6-31G*算ESP，然后套RESP……其实不是不能用，但这样出来的电荷和GAFF/AMBER的默认组合不匹配。一旦后续兼容性出了问题，定位起来很蛋疼。建议统一用HF/6-31G*。
3.2 电荷约束的陷阱——别乱加枷锁
RESP分为两个阶段（两步）：
- 阶段1：不加任何原子等效约束，自由拟合出电荷
- 阶段2：对等效原子（比如甲基上的三个氢）施加等价约束，重新拟合
等效约束的目的是让对称性相同的原子带相同电荷，这对保持分子旋转对称性很重要。但问题来了——很多人走到第二步就把所有碳氢键都约束成一样，结果电荷极度失真。
举个例子，苯环的六个碳在RESP阶段2应该被约束为等价吗？答案是：如果苯环是平面的、没有取代基破坏对称性，就可以。但一旦有一个取代基（比如对位有-OH），六个碳的电荷其实不一样，再强行约束成等价，你拟合下的RMSD会飙高。
我常用的约束规则（以AMBER的resp程序为例）：
输入文件：chg.in
&files
-i // 第一步无约束作业
-o // 输出文件
/
&resp
iqopt=2 ! 使用两阶段流程
iqwt=0.001 ! 约束强度
/
二阶段约束写法（tleap格式）：
Residue: LIG
Atom1 1 6 0.0 ! 等效组1：三个H
Atom2 2 7 0.0
Atom3 3 8 0.0
.....小技巧：约束强度 iqwt 默认0.001（对应AMBER），数值越大约束越紧。如果你发现等效原子电荷差太多，又非约束不可，可以稍微调高到0.01试试，别超过0.1。
核心原则：只有经过结构对称性验证的原子才能加等效约束。拿不准时，先不加约束跑阶段1，看RMSD和电荷值再决定。
3.3 多重构象电荷平均——不能简单平均坐标
一个柔性分子有多个低能构象，你用哪一个构象的电荷？答案是：全部都要，然后做ESP平均后的RESP拟合。
具体流程我总结为三步：
- 对每个构象分别做 HF/6-31G* 结构优化（固定成该构象的局部极小点）
- 分别计算 每个构象的ESP格点（格点坐标完全相同，但用的是该构象的波函数）
- 将所有构象的ESP格点合并成一个大数据集，然后对这个合并数据集做RESP拟合
这里最坑的是第二步——格点坐标必须保持一致。我曾在项目中遇到一个分子有8个构象，分别算ESP时默认格点中心的原点设置不同，导致合并后每个构象的格点不重叠。最后拟合出的电荷完全无法用。后来我用了固定格点原点的方法：
在Gaussian的输入文件里添加：
IOp(6/50=1) ! 使用用户指定的格点原点
IOp(6/51=0, 0, 0) ! 将原点固定在大分子质心（一般设置分子质心处）或者在Antechamber里用 antechamber -fi gout -fo prepi -c resp -at gaff -nc 0 ... 时，如果分子有柔性，可以先生成一个包含所有构象ESP的“残基ESP文件”，再统一运行resp程序。
避坑指南：我曾经直接把8个构象的坐标平均成一个伪结构，然后算ESP……结果电荷完全不对，因为分子内氢键位置变化导致电荷分布完全变了。千万不要为了省事去平均坐标！必须平均ESP格点。
多重构象的数量多少合适？我个人的习惯是：分子内可旋转键数量≤3时，取3~5个构象；可旋转键≥4时，至少取8~10个构象。超过15个构象的收益开始递减，而且计算成本飙升，建议用自动构象搜索工具（如OpenEye Omega, RDKit ETKDG）生成。
3.4 实战参数速查表
| 参数/步骤 | 推荐值/操作 | 说明 |
|---|---|---|
| QM方法/基组 | HF/6-31G* | 与GAFF/AMBER配套 |
| ESP采点方法 | CHELPG (Gaussian) | 比CHelpG更可靠 |
| 格点密度 | 2~4 points/Bohr² | IOp(6/33=2) |
| RESP约束强度 | 0.001 (默认) | 阶段2使用 |
| 等效约束对象 | 甲基H、苯环对称位 | 其他位置慎用 |
| 多重构象数 | 5~10个 | 覆盖低能构象空间 |
| 电荷验证 | RMSD<0.1 kcal/mol | 偶极矩偏差<10% |
3.5 便捷脚本示例
如果你用AmberTools，以下命令可以一键走完单构象RESP流程（假设你已经有了Gaussian的.fchk和.log文件）：
antechamber -i ligand.log -fi gout -o ligand.prep -fo prepi \
-c resp -at gaff -nc 0 -j 4 -rn LIG对于多重构象，则需手动准备一个ESP合并文件。我写过一个简单的bash脚本片段（仅供参考）：
#!/bin/bash
# 合并多个构象的ESP（已固定格点原点）
for conf in conf*.log; do
antechamber -i $conf -fi gout -o ${conf%.log}.esp -fo esp
done
# 用resp程序手动合并
cat conf1.esp conf2.esp conf3.esp > all.esp
resp -i all.esp -o charge.out -p ligand.pdb -q个人建议：第一次跑之前，先用一个构象验证整个流程通不通。再看RMSD和偶极矩，没问题后再上多重构象——否则中间出bug排查起来真的很崩溃。
对了，还有一点需要注意：电荷总和的整数性。RESP拟合出来的电荷总和可能不是精确的整数（比如中性分子不是0.000），需要做一个微调。AmberTools的resp程序会自动修正，但如果你自己写脚本，记得加上最后一步的电荷归整。
以上就是RESP电荷计算的核心内容。参数设置其实就那么几个关键点，但每个点里都藏着无数坑。我踩过的那些、同行踩过的那些，上面基本都提到了。电荷搞定了，力场参数的骨架才算真正立起来——剩下的就是二面角、范德华那些事了。