EDA-FF：基于分子力场的能量分解
在研究分子间弱相互作用时，单看一个总结合能往往不够。比如两个分子形成复合物，到底是静电作用占主导，还是色散作用占主导？哪些原子促进结合，哪些原子反而产生排斥？这些问题通常需要借助能量分解分析（energy decomposition analysis, EDA）回答。
常见的 SAPT、ALMO-EDA、sobEDA 等方法从电子结构计算出发，物理图像较完整，但计算量也明显高于普通的单点能计算。对于较大的体系，这些方法并不适用。为解决此问题，sobereva提出了 EDA-FF（energy decomposition analysis based on forcefield），直接使用经典分子力场的非键相互作用表达式，把片段间作用拆成静电、排斥和色散三项，用于迅速得到弱相互作用成分。
EDA-FF 分解的不是完整的量子化学相互作用能，而是所选原子电荷和力场参数定义出的经典非键相互作用能。
原理
对于一个体系，假设其被划分为片段 $A$ 和片段 $B$。EDA-FF 遍历所有跨片段原子对，即 $i\in A$、$j\in B$，计算每一对原子的非键相互作用：
\[E_{ij}^{\mathrm{nb}} = E_{ij}^{\mathrm{ele}} + E_{ij}^{\mathrm{rep}} + E_{ij}^{\mathrm{disp}}\]随后对所有跨片段原子对求和：
\[E_{AB}^{X} =\sum_{i\in A}\sum_{j\in B}E_{ij}^{X}, \qquad X\in\{\mathrm{ele},\mathrm{rep},\mathrm{disp}\}\]因此总的力场相互作用能为：
\[E_{AB}^{\mathrm{FF}} = E_{AB}^{\mathrm{ele}} + E_{AB}^{\mathrm{rep}} + E_{AB}^{\mathrm{disp}}\]如果定义了多个片段，只需再对所有不同片段对求和：
\[E_{\mathrm{int}}^{X} = \sum_{a<b}\sum_{i\in F_a}\sum_{j\in F_b}E_{ij}^{X}\]因此，EDA-FF 算法本质上是一个跨片段原子对势的求和。只要有几何结构、原子电荷和力场原子类型，就能完成分析。且得益于力场势函数的简单形式，EDA-FF 计算可以以非常小的代价完成。
能量项
静电作用
两个原子的静电相互作用使用点电荷库仑势表示：
\[E_{ij}^{\mathrm{ele}} = k_{\mathrm e}\frac{q_iq_j}{r_{ij}}\]其中:
- $q_i$ 和 $q_j$ 是原子电荷
- $r_{ij}$ 是原子间距离
- 在原子单位制中 $k_{\mathrm e}=1$
显然，异号电荷给出负值，产生吸引；同号电荷给出正值，产生排斥；作用按 $1/r$ 衰减。但这里的原子电荷并不是一个唯一的量子力学可观测量。不同电荷划分方法可能给出明显不同的数值，因此 EDA-FF 的静电项会直接依赖所选电荷模型。静电项实际上想近似的是一个片段在另一个片段所在空间位置产生的静电势。因此，所选原子电荷应当尽量重现分子范德华表面附近的静电势，如MK电荷、CHELPG电荷等。对于体系大到难以使用波函数计算的情况，可以考虑使用EEM电荷。
这里并不推荐使用RESP电荷。RESP电荷的优势是可以解决柔性体系等价原子电荷不一致的问题，但在EDA-FF场景下，体系是静态的，实际上并不需要RESP电荷进行的约束处理；相反，RESP的约束还会略微削弱MK电荷对静电势的重现性。
范德华作用（交换互斥项与色散项）
Multiwfn的EDA-FF 使用常见的 12–6 Lennard-Jones 形式描述范德华作用：
\[E_{ij}^{\mathrm{vdW}} = E_{ij}^{\mathrm{rep}} + E_{ij}^{\mathrm{disp}}\]其中排斥部分为：
\[E_{ij}^{\mathrm{rep}} = \varepsilon_{ij} \left(\frac{R_{ij}^{0}}{r_{ij}}\right)^{12}\]色散吸引部分为：
\[E_{ij}^{\mathrm{disp}} = -2\varepsilon_{ij} \left(\frac{R_{ij}^{0}}{r_{ij}}\right)^6\]合并后得到：
\[E_{ij}^{\mathrm{vdW}} = \varepsilon_{ij}\left[\left(\frac{R_{ij}^{0}}{r_{ij}}\right)^{12} - 2\left(\frac{R_{ij}^{0}}{r_{ij}}\right)^6\right]\]这里 $\varepsilon_{ij}$ 是势阱深度，$R_{ij}^{0}$ 是势能最低点对应的原子间距离。当 $r_{ij}=R_{ij}^{0}$ 时：
\[E_{ij}^{\mathrm{rep}} = +\varepsilon_{ij}\] \[E_{ij}^{\mathrm{disp}} = -2\varepsilon_{ij}\] \[E_{ij}^{\mathrm{vdW}} = -\varepsilon_{ij}\]
12–6 势中排斥项、色散项和总范德华作用的关系。
这张图也解释了为什么几何结构对 EDA-FF 极其重要：排斥项按 $r^{-12}$ 变化，两个原子稍微靠近一点，正的排斥能就可能迅速爆炸；色散项按 $r^{-6}$ 衰减，变化相对平缓。
EDA-FF 输出中的 repulsion 只是用 $r^{-12}$ 势模拟的短程排斥。它可以近似对应 Pauli/交换互斥的物理效果，但并不等同于 SAPT 或其它量子化学 EDA 中严格定义的交换项。
不同力场计算参数的方式不同。UFF 的规则为：
\[\varepsilon_{ij} = \sqrt{\varepsilon_i\varepsilon_j}\] \[R_{ij}^{0} = \sqrt{R_i^{0}R_j^{0}}\]AMBER/GAFF 则使用：
\[\varepsilon_{ij} = \sqrt{\varepsilon_i\varepsilon_j}\] \[R_{ij}^{0}=R_i^{*}+R_j^{*}\]其中 $R_i^{*}$ 是力场为单个原子类型定义的非键半径。
对于普通有机体系，AMBER 或 GAFF 通常比 UFF 更适合做这类分析。UFF 最大的优点是元素覆盖广，但其通用参数有时会在量子化学优化的结构上严重高估短程排斥。若体系含有 AMBER/GAFF 不支持的元素，必须谨慎补充参数，不能因为程序成功输出了数字就默认结果合理。
原子贡献
EDA-FF 的另一个特点是可以输出每个原子对片段间相互作用的贡献。对于能量分量 $X$，原子 $i$ 的贡献可理解为：
\[C_i^X = \frac{1}{2}\sum_{j,\,\mathrm{frag}(j)\neq\mathrm{frag}(i)}E_{ij}^{X}\]也就是把每个跨片段原子对能量平均分给参与该原子对的两个原子。于是：
\[\sum_i C_i^X=E_{\mathrm{int}}^X\]应用
multiwfn
multiwfn支持基于uff、gaff、amber的EDA-FF。详见sobereva老师的博文使用Multiwfn做基于分子力场的能量分解分析
baneda
作为懒狗，笔者觉得multiwfn的输入太复杂了，因此基于笔者先前写好的一些库做了一个baneda程序进行便捷的EDA-FF计算。
baneda推荐的输入是片段的chg文件。用户可以提前将体系的片段计算好MK电荷，产生chg文件。如果嫌这一步也麻烦，可以利用banetask一键完成：
1
2
3
4
5
6
7
8
9
10
$dimeropt
%keywords opt freq b3lyp 6-31g* em=gd3bj
@foreach? frag
$[frag.name]
%keywords nosymm b3lyp 6-31g* em=gd3bj
%process
banewfn mk.bwc [inputname].fchk
%artifact charge.[frag.name] mk.chg
@end
随后，将交给baneda：
1
baneda frag1.chg frag2.chg --base-types gaff
baneda的–base-types可以基于内置的API自动分辨原子类型，如此即可完成EDA-FF计算：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
EDA-FF energy decomposition (kJ/mol)
Inputs: 2
[1] frag1.chg CHG frame 1/1 atoms 3 topology none
[2] frag2.chg CHG frame 1/1 atoms 3 topology none
Atoms: 6
Fragments: 2
Electrostatic model: 1/r
LJ model: typed Rmin/2 additive
Fragment definitions:
1 frag1.chg (3 atoms): 1-3
2 frag2.chg (3 atoms): 4-6
Charge assignment:
Fragment 1 frag1.chg:
source: input charges
calculated total: 0.000000 e
Fragment 2 frag2.chg:
source: input charges
calculated total: 0.000000 e
Type sources: base force field=6
Fragment-pair energies:
Frag A Frag B Electrostatic Repulsion Dispersion vdW Total
1 2 -0.405080 0.001119 -0.053343 -0.052224 -0.457304
Total over unique fragment pairs: Elec=-0.405080 Rep=0.001119 Disp=-0.053343 vdW=-0.052224 Total=-0.457304
对于体系很大，难以计算原子电荷的情况，baneda内置了一套EEM电荷参数集，此时输入文件需要包含拓扑，推荐使用obabel进行转换：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
[warning] baneda.charge.input-overridden: fragment 'C3' specifies net charge 0; its input atomic charges are ignored and replaced by EEM charges
[warning] baneda.charge.input-overridden: fragment 'C' specifies net charge 0; its input atomic charges are ignored and replaced by EEM charges
[warning] baneda.charge.input-overridden: fragment 'G' specifies net charge 0; its input atomic charges are ignored and replaced by EEM charges
EDA-FF energy decomposition (kJ/mol)
Input: C3GC.mol2
Format: MOL2
Frame: 1/1
Topology: input bond table
Atoms: 101
Fragments: 3
Electrostatic model: 1/r
LJ model: typed Rmin/2 additive
Fragment definitions:
1 C3 (72 atoms): 30-101
2 C (13 atoms): 1-4,10-12,18,20-24
3 G (16 atoms): 5-9,13-17,19,25-29
Charge assignment:
Fragment 1 C3:
source: EEM
target charge: 0
parameter set: B3LYP/6-31G* MK EEM (2009)
topology: input bond table
calculated total: -0.000000 e
reciprocal condition estimate: 2.024448e-03
relative residual: 7.486468e-18
Fragment 2 C:
source: EEM
target charge: 0
parameter set: B3LYP/6-31G* MK EEM (2009)
topology: input bond table
calculated total: 0.000000 e
reciprocal condition estimate: 2.318199e-02
relative residual: 1.264575e-17
Fragment 3 G:
source: EEM
target charge: 0
parameter set: B3LYP/6-31G* MK EEM (2009)
topology: input bond table
calculated total: -0.000000 e
reciprocal condition estimate: 2.197031e-02
relative residual: 1.167326e-17
Type sources: base force field=101
Fragment-pair energies:
Frag A Frag B Electrostatic Repulsion Dispersion vdW Total
1 2 9.084336 44.859021 -94.691387 -49.832366 -40.748030
1 3 4.971214 62.080949 -132.623501 -70.542552 -65.571338
2 3 -75.274858 60.255979 -45.541447 14.714532 -60.560326
Total over unique fragment pairs: Elec=-61.219308 Rep=167.195949 Disp=-272.856334 vdW=-105.660385 Total=-166.879693
不要用obabel产生的，那个是照着NPA电荷拟合的，不适合用于能量分解计算。
EDA-GFN-FF
EDA-FF还可以基于更优秀的力场进行计算。GFN-FF无需指定原子类型，也无需提供原子电荷（xTB有一套自己的EEQ电荷拟合算法），非常适合简易计算。笔者基于GFN-FF开源库简单做了一个EDA-GFN-FF程序，可以在GitHub查看。
除了前述片段定义方式，EDA-GFN-FF也可以接受Gaussian格式的输入文件，预先定义好片段，直接传给程序：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
Input: C3GC.gjf
Format: gaussian
Atoms: 101
Total charge: 0
Multiplicity: 1
Fragments: 3
Fragment definition
Frag Atoms Charge
1 13 0
2 16 0
3 72 0
GFN-FF total energy: -2.41834422568402E+01 Eh
GFN-FF gradient norm: 1.15814031067106E-01 Eh/a0
EDA-GFNFF interfragment decomposition (kcal/mol)
FragA FragB Electrostatic Repulsion Dispersion H-bond X-bond Total NCI
1 2 -23.820722 7.630888 -2.853420 -15.216343 0.000000 -34.259597
1 3 0.000372 1.532240 -11.865242 -0.155619 0.000000 -10.488249
2 3 0.027754 1.881684 -16.177708 -0.202090 0.000000 -14.470360
------------------------------------------------------------------------------------------------------------
All pairs -23.792595 11.044812 -30.896370 -15.574053 0.000000 -59.218207
Total interfragment EDA-GFNFF
Electrostatic: -3.79159140728572E-02 Eh
Repulsion: 1.76010278386348E-02 Eh
Dispersion: -4.92365000879241E-02 Eh
H-bond: -2.48188334544019E-02 Eh
X-bond: 0.00000000000000E+00 Eh
Three-term: -6.95513863221465E-02 Eh
Total NCI: -9.43702197765484E-02 Eh
Total NCI: -59.218207 kcal/mol
Total NCI: -247.768978 kJ/mol
GFN-FF具有特色的氢键项和卤键项，因此在计算相关体系会比普通力场好得多。例如当前测试体系中，利用B3LYP-D3(BJ)/6-311+G**结合counterpoise校正计算出，有氢键的Frag1和Frag2之间的结合能是-34.409656 kcal/mol，跟GFN-FF的计算值非常接近，说明H-bond校正项相当有效。（当然，其他两个片段计算值有一定差异）
共 1 个文件，位于 /assets/posts/65/。