Mayer能量分解方法及其在Amesp中的使用
能量分解方法能够将分子之间的相互作用分解成各种具有物理意义的成分，比如键能、静电相互作用、排斥作用、极化、诱导、电荷转移以及色散作用等。常用的能量分解方法有如Kitaura-Morokuma方法1、Ziegler-Rauk方法2、对称匹配微扰（SAPT）方法3、LMO-EDA方法4、GKS-EDA方法5、自然能量分解（NEDA）方法6等。而本文将介绍可以获得分子中原子的能量以及原子对之间的相互作用的Mayer能量分解方法7及其在Amesp中的使用。
1 理论方法
本小节将介绍Mayer能量分解的原理，体系的Hartree-Fock总能量为：
其中D为总的密度矩阵：
h为单电子哈密顿项：
在Mayer能量分解中，分子中原子A的能量EA为：
而原子对A和B之间的相互作用EAB为：
公式（5）中包含原子A和B之间的动能（第一项）、电子-核吸引能（第二项和第三项）、电子-电子之间的静电相互作用（第四项）以及排斥作用（第五项），而体系的总能量写为：
其中EAB可以用来表述键能，只需要执行一次计算便可估算所有化学键之间的键能，从而能够使用键能判断分子的稳定性。
使用(3)式和(4)式可以将Hartree-Fock的能量进行分解，也可以对DFT计算得到的波函数来使用(3)式和(4)式进行能量分解计算，但是总能量(6)式与DFT总能量是不一致的，因为其中少了交换相关项。在Vyboishchikov等人8的工作中，他们将交换相关项加入进Mayer能量分解的框架中，交换相关能的表达式为：
将其分解到
其中分解到原子A的电子密度为：
由于每个电子的交换相关能密度（the exchange-correlation energy density per electron）εxc(r)关于电子密度是非线性的，因此εAxc(r)≠εxc[ρA(r)]。在Vyboishchikov等人的工作中，εxc(r)使用一组以原子为中心的辅助基函数进行展开，而εAxc(r)则以原子A为中心的辅助基函数表示：
在(11)式中，ξk为待定的拟合系数，使用最小二乘法求得。在Amesp中，为保证总能量在拟合过程中不变，添加了以下约束条件：
求解如下线性方程组，即可得到拟合系数ξk：
上式中：
值得注意的是，在εxc(r)中乘以一个权重函数w(r)不影响总能量的结果：
根据(17)式将(15)式变为：
这里w(r)的选取具有一定随意性，(15)式中选择的是w(r)=1，但是这种情况下拟合得到的ξk均匀的覆盖到了包括原子核附近和远离原子核的部分，这使得原子核附近和远离原子核的部分区分度降低，因此w(r)=1会导致结果完全错误。为了区分出近核和远核的区域，Vyboishchikov等人选取w(r)=ρ(r)1/2，求解出拟合系数ξk后，εAxc(r)写为：
最后整理得到DFT的Mayer能量分解的表达式为：
2 Mayer能量分解在Amesp中的使用
这里介绍一个简单的使用Amesp计算NH3分子Mayer能量分解的例子，其输入为：
% npara 4
! b3lyp def2-TZVPP
>method
eda mayer
end
>xyz 0 1
N 0.80283598 1.70165446 -0.00000000
H 1.17826884 0.76097872 0.00000000
H 1.17830444 2.17209355 0.81462077
H 1.17830444 2.17209355 -0.81462077
end上述计算默认的辅助基组为def2/J，若想使用其他的辅助基组，需要在>ope模块做如下修改：
>ope
dfset def2-JK
end其中dfset后面为具体的辅助基组，也可以使用自定义基组，关键词为define，做法可以参考《Amesp中基组的使用》。输出的结果部分为：
Total energy part of Mayer's EDA:
1 2 3 4
1 N -54.67098761
2 H -0.17093239 -0.47193843
3 H -0.17088434 0.00361248 -0.47189135
4 H -0.17088434 0.00361248 0.00357257 -0.47189136其中对角线部分为EA，非对角线部分是EAB，如EN-H = -0.17093239 a.u. = -107.26 kcal/mol。这是对总能量进行分解，若想获得每一部分（如动能，核吸引能，交换相关能等）的能量分解，可以在>ope模块中添加out 1关键词即可：
% npara 4
! b3lyp def2-TZVPP
>method
eda mayer
end
>ope
out 1
end
>xyz 0 1
N 0.80283598 1.70165446 -0.00000000
H 1.17826884 0.76097872 0.00000000
H 1.17830444 2.17209355 0.81462077
H 1.17830444 2.17209355 -0.81462077
end此时的输出结果为：
Kinetic energy part of Mayer's EDA:
1 2 3 4
1 N 53.60060679
2 H 0.47320842 0.43047686
3 H 0.47338246 -0.00174079 0.43058376
4 H 0.47338246 -0.00174079 -0.00175097 0.43058375
Nuc-Ele energy part of Mayer's EDA:
1 2 3 4
1 N -129.34760060
2 H -7.28639164 -1.03297637
3 H -7.28750256 -0.56952335 -1.03326085
4 H -7.28750257 -0.56952335 -0.56967776 -1.03326084
Nuc-Nuc energy part of Mayer's EDA:
1 2 3 4
1 N 0.00000000
2 H 3.65732408 0.00000000
3 H 3.65718893 0.32477379 0.00000000
4 H 3.65718893 0.32477379 0.32479973 0.00000000
EleCoul energy part of Mayer's EDA:
1 2 3 4
1 N 27.41172240
2 H 3.37101696 0.30539460
3 H 3.37220054 0.25633131 0.30560249
4 H 3.37220055 0.25633130 0.25644353 0.30560248
HF-Exch energy part of Mayer's EDA:
1 2 3 4
1 N -1.21581671
2 H -0.07748585 -0.02922896
3 H -0.07749506 0.00114666 -0.02924850
4 H -0.07749506 0.00114666 0.00114723 -0.02924849
DFT-Exc energy part of Mayer's EDA:
1 2 3 4
1 N -5.11989949
2 H -0.30860437 -0.14560456
3 H -0.30865864 -0.00737513 -0.14556826
4 H -0.30865865 -0.00737513 -0.00738918 -0.14556826
Total energy part of Mayer's EDA:
1 2 3 4
1 N -54.67098761
2 H -0.17093239 -0.47193843
3 H -0.17088434 0.00361248 -0.47189135
4 H -0.17088434 0.00361248 0.00357257 -0.47189136需要注意的是EAB并不严格等于键能，但是其可以用于方便估算键能。若只想使用DFT的波函数来使用(3)式和(4)式(Hartree-Fock)进行能量分解计算，只需要在>ope模块中添加mayerdft off关键词即可，值得注意的是，此时的分解后相加得到的总能量和DFT能量是不一致的。对于计算片段和片段之间的能量分解，只需要把与该片段相关的原子求和即可。
3 参考文献
1. Int. J. Quantum Chem. 10, 325 (1976). 2. Inorg. Chem. 18, 1558 (1979). 3. J. Mole. Struct. 547, 293 (2001). 4. J. Chem. Phys. 131, 014102 (2009). 5. J. Phys. Chem. A 118, 2531 (2014). 6. J. Chem. Phys. 100, 2900 (1994). 7. Chem. Phys. Lett. 382, 265 (2003). 8. J. Chem. Phys. 122, 244110 (2005).