Mayer能量分解及APEX4程序使用简介

1 前言

目前已有一大批能量分解方法，用于将相互作用能拆成各个成分，如Morokuma方法、Ziegler-Rauk方法、对称匹配微扰(SAPT)方法、自然能量分解方法(NEDA)等等。这些方法在原理、算法上有所差异，共同特点是都需要将体系拆成两个片段，将其在孤立状态作为“参考态”，然后在此基础上经过一系列处理得到各个能量项。从这个角度看，Mayer提出的能量分解方法CECA (Chemical Energy Component Analysis，见CPL,332,381)与上述方法有很大差异。CECA不需要人为定义片段，它可以直接获得体系中每个原子的能量以及每一对儿原子间的能量（如果想获得片段间作用能将相应原子对儿作用能加和即可），每一对儿原子间的能量可以被分解为"动能"+"静电"+"交换"+"重叠"项(见Theor. Chem. Acc.,109,91)。另外CECA也不需要设定参考态，它只利用量子化学计算的整个体系的单粒子密度矩阵信息。最初由于CECA对三、四中心电子积分只是近似地考虑，导致全部原子能量加上全部原子对儿作用能与总能量略有偏差，后来分解算法经过修改（见CPL,382,265），使得总和与总能量精确一致，在下文将统称为CECA，仅当需要进行区别时会特意称为“精确CECA”。由于在Mayer的文章中对CECA算法描述得比较复杂混乱，建议参考JCP,129,144111中对Mayer方法的介绍，原理基本等价但思路不同，简单清晰得多，下面也以这种方式对CECA进行描述。

2 原理

体系的密度函数可以这样写为原子成分的加和
ρ(r)   =∑[a]ρ_a(r)     其中ρ_a(r)   =∑[i∈a]∑[j]P_i,j*χ_j(r)*χ_i(r)，χ代表基函数，P是密度矩阵
体系的密度矩阵函数也可以拆解为原子的加和
ρ(r|r')=∑[a]ρ_a(r|r')  其中ρ_a(r|r')=∑[i∈a]∑[j]P_i,j*χ_j(r)*χ_i(r')
体系Hartree-Fock总能量可写为
E=核互斥能+电子动能+核吸引势能+电子互斥能+交换能
 =∑[a]∑[b>a]Z_a*Z_b/R_a,b + ∑[i]∑[j]P_i,j*<i|T|j> - ∑[a]∫Z_a*ρ(r)/|R_a-r|dr
  + (1/2)*∫∫ρ(r)*ρ(r')/|r-r'|drdr' - (1/2)*∫∫ρ(r|r')*ρ(r'|r)/|r-r'|drdr'
将其中ρ(r)和ρ(r|r')都替换为原子成分加和的形式，就会看到上式后三项都可以写为∑[a]∑[b]这样原子对儿加和的形式。核互斥能已经是这种形式，而电子动能项也可以写成这种形式，即∑[a]∑[b]∑[i∈a]∑[j∈b]P_i,j*<i|T|j>。因此，HF总能量就可以写为原子对儿加和形式
E=∑[a]∑[b>a]E_a,b+∑[a]E_a,a
这里E_a,b代表a、b原子间的作用能，E_a,a是原子a在实际体系中自身的能量。E_a,b表达式为
E_a,b=2∑[i∈a]∑[j∈b]P_i,j*<i|T|j> + Z_a*Z_b/R_a,b - ∫Z_a*ρ_b(r)/|R_a-r|dr - ∫Z_b*ρ_a(r)/|R_b-r|dr
      + ∫∫ρ_a(r)*ρ_b(r')/|r-r'|drdr' - ∫∫ρ_a(r|r')*ρ_b(r'|r)/|r-r'|drdr'
第一项是原子间相互作用的动能贡献，第2~5项是静电贡献，最后一项是交换贡献。

虽然严格来说，KS-DFT计算的体系的能量分解还得额外考虑怎么把交换相关泛函积分划分为原子对儿加和形式，但是如果只将KS-DFT认为是用于提供比较好的密度（矩阵）函数的话，那么CECA的形式也可以直接用在KS-DFT波函数上。尽管此时显然各项加和的结果不等于量子化学计算输出的总能量，因为忽略了相关能，且计算交换能的方式不同。

原子a在实际体系中的能量E_a,a相对于它在自由状态的能量的变化可正可负，称为preparation能量，一方面来自于电子密度的极化，另一方面来自于电子在原子间的净转移。前者总是使原子能量升高，而后者，对于失电子的原子会使其能量升高，而使得电子的原子能量会降低。E_a,b和键的强度、键的解离能有很直接的关联，在某种意义上类似于键级。由于CECA的结果是依赖于几何结构的，研究键的强度的时候应当用稳定的几何结构。注意E_a,b并不精确等于键的解离能，因为键的解离能是以原子在自由状态为参考的，而E_a,b是以原子在实际体系中的状态为参考，但CECA仍不失为一种十分省事的估算体系内各化学键键能的方法。

CECA对密度矩阵的划分实际上用的就是Mulliken布居分析用的平均划分形式。正如同Mulliken布居对基组敏感，CECA的结果也有一定基组依赖性，且不随基组质量提高而收敛。不过由于E_a,b、E_a,a的数量级较大，基组造成的影响从相对值上来看并不大。然而由于preparation能量数量级较小，基组的改变可能使其正负号发生变化。哪种基组下结果是最合理的，这没有确切答案。但至少含弥散函数的基组切不可用，因为弥散函数会明显破坏基组的平衡性，这也是为什么用弥散函数后Mulliken电荷往往很烂。笔者认为，CECA可以尝试改用以自然原子轨道(NAO)为基的形式，这样基组依赖性问题会小得多，结果会更可靠。

3 APEX4程序的使用

APEX4是实现CECA的程序，使用简单，运行速度快，仅需要Gaussian（G92~G09皆兼容）的fch文件即可，执行中不需要再调用量子化学程序，因为程序包里面本身就带了电子积分代码。缺点是不支持f及以上角动量函数。RHF、UHF波函数都支持，也支持KS-DFT波函数（注意上文对KS-DFT情况的讨论）。APEX4源代码可以从这里免费下载http://occam.ttk.hu/programs/。APEX4具体的执行方法和Mayer的原文，以及上面对CECA的原理介绍有一定出入，但最根本的原理都是差不多的。

编译过程很简单。
对于Linux下的编译，下载并解压后，将Makefile里的f77改成自己喜欢的编译器。然后在此目录下执行make命令即可。我这里用ifort 11.0可以顺利编译。
对于Windows下的编译，比如用CVF6.5，将所有文件加入工程里，然后build即可。我编译好的1.02版可由此下载http://pan.baidu.com/s/1hq3fLHE
 
APEX4执行也很容易。将体系的fch文件改名为Test.FChk后，放到与APEX4可执行文件相同路径下，在命令行下面（对于Windows就是MSDOS环境）启动APEX4即可。程序就会分析Test.FChk并输出能量分解矩阵。对于小体系几秒钟就能出结果。我提供的压缩包里的Test.FChk是乙烷在HF/6-31G**下面优化任务对应的fch文件。（注：可以在route section直接写上formcheck关键词，%chk那行可以不写，对于windows版Gaussian，任务完毕后就会直接在scratch目录下生成Test.FChk文件。这就免得先生成chk再用formchk转换再改名了。）

程序最先输出的是Mayer键级矩阵和原子价，随后输出各种各样的能量矩阵。其中包括以下四种总能量矩阵，等号右边是它们的非对角元和其它矩阵的运算关系。这四种矩阵的对角元就是各个E_a,a，非对角元就是E_a,b（注：E_a,b=E_b,a）。全部对角元的加和以及全部非对角元的加和的一半（即全部原子间作用能）也会在程序中输出。

"CECA"   = ELECTROSTATIC + EXCHANGE + OVERLAP + FINITE-BASIS CORRECTIONS
"CECA + 3- and 4-center contributions" = CECA + Remaining 3- and 4-center
"CECA/T" = ELECTROSTATIC + EXCHANGE + OVERLAP + KINETIC ENERGY
"Exact"  = "CECA/T" + Remaining 3- and 4-center （注意这只是形式上的关系，"Exact"计算时并不是靠这样加和来获得的）

它们的差异可以从两个角度来看：
(1)精确与近似分解的角度："CECA/T"和"CECA"都是对总能量的近似分解，它们对三、四中心积分都通过近似方式处理，导致上三角（或下三角，后同）矩阵元加和不精确等于SCF能量（尽管十分接近）。这差异可以通过加上"Remaining 3- and 4-center"项来修正（其数值很小），这样上三角矩阵元加和就精确等于SCF能量了。而"Exact"，也就是之前提到的"精确CECA"，用了不同分解算法，没有对三、四中心积分进行近似处理，所以"Exact"的上三角矩阵元加和直接就是SCF能量。
(2)对动能的考虑："CECA"和"CECA + 3- and 4-center contributions"都只将电子动能分解到各个单中心项（E_a,a）里，而"CECA/T"和"Exact"则将电子动能同时分解到单中心项和双中心项（正如前面的原理介绍所述），这也使得它们不需要"FINITE-BASIS CORRECTIONS"这一项。由于这个原因，"CECA/T"和"Exact"的非对角元明显没有"CECA"和"CECA + 3- and 4-center contributions"的那么负，而对角元则更负。

从结果的意义上来讲，不要用"CECA"和"CECA + 3- and 4-center contributions"矩阵，原子间作用项过大，化学意义差。"CECA/T"和"Exact"数值很接近，如果想精确地将总能量分解为单中心和双中心项，就用"Exact"。由于"Remaining 3- and 4-center"这个东西物理意义很不清楚，所以如果想将单、双中心项再分解为各种类型能量成分，就用"CECA/T"及构成它的各个矩阵。换句话说，APEX4里面具体用的方法没法像第二节介绍的那样，既保持总能量精确，又使能量各个成分都有清楚的物理解释。

程序还会输出这样的矩阵"ELECTROSTATICS: POINT-CHARGE APPROXIMATION"，它代表将体系电荷的连续分布用原子电荷模型代替后计算出的静电作用矩阵。"ELECTROSTATICS: DEVIATION FROM P.-CH. APPROIMATION"矩阵就是"ELECTROSTATIC CONTRIBUTIONS"矩阵与点电荷近似矩阵的差值。当原子间距离越远，静电相互作用能就可以越好地通过点电荷模型来表达，因此这个偏差矩阵数值就会越小。

4 例子

这个是HF/6-31G**优化下的交错式乙烷的"Exact"能量矩阵，可见C-C键作用能是-0.1895 a.u.。氢原子间虽然不直接相连，但也有弱相互作用。同一个甲基的氢，如H6-H7，呈现0.0136 a.u.互斥作用，体现了位阻效应。而一个甲基氢与另一个甲基的位置相反的氢（中心反演对称），比如H2-H6，呈现-0.001 a.u.（-0.627kcal/mol）稳定化作用。其中静电项为0.0015 a.u.，为正值这不难理解，因为两个氢电性相同；使能量降低的主要因素来自于动能项，为-0.0024 a.u.。H2-H6氢相互作用能为负值在某种程度上也可以解释为C1-H2成键向另一个C5-H6反键空轨道的离域作用（NBO E2分析为-3.56kcal/mol）。
   "Exact" energy decomposition matrix
             1  C     2  H     3  H     4  H     5  C     6  H     7  H     8  H

   1  C  -37.6605  -0.1722  -0.1722  -0.1722  -0.1895   0.0071   0.0071   0.0071
   2  H   -0.1722  -0.4738   0.0136   0.0136   0.0071  -0.0010   0.0045   0.0045
   3  H   -0.1722   0.0136  -0.4738   0.0136   0.0071   0.0045  -0.0010   0.0045
   4  H   -0.1722   0.0136   0.0136  -0.4738   0.0071   0.0045   0.0045  -0.0010
   5  C   -0.1895   0.0071   0.0071   0.0071 -37.6605  -0.1722  -0.1722  -0.1722
   6  H    0.0071  -0.0010   0.0045   0.0045  -0.1722  -0.4738   0.0136   0.0136
   7  H    0.0071   0.0045  -0.0010   0.0045  -0.1722   0.0136  -0.4738   0.0136
   8  H    0.0071   0.0045   0.0045  -0.0010  -0.1722   0.0136   0.0136  -0.4738
对于乙烯，C-C作用能是-0.266 a.u.，对于乙炔是-0.485 a.u.，随键强增加作用能逐渐加大，比较符合化学直觉。

不过，CECA方法的可靠性目前被检验得还不够广泛，对于一些体系不敢说可靠。比如HF/6-31G**下的甲酰胺，C-O是双键，Mayer键级为1.84，而C-N单键键级为1.07，可是相应的作用能却是C-N比C=O要大，违背了化学直觉，所以使用CECA的时候需要慎重。
   "Exact" energy decomposition matrix
             1  N     2  H     3  H     4  C     5  H     6  O

   1  N  -54.2605  -0.2359  -0.2402  -0.3997   0.0210   0.1095
   2  H   -0.2359  -0.4219   0.0375   0.0659   0.0061  -0.0301
   3  H   -0.2402   0.0375  -0.4166   0.0700   0.0015  -0.0390
   4  C   -0.3997   0.0659   0.0700 -37.4528  -0.1686  -0.2831
   5  H    0.0210   0.0061   0.0015  -0.1686  -0.4638   0.0029
   6  O    0.1095  -0.0301  -0.0390  -0.2831   0.0029 -74.8428

5 改进的化学能量成份分析

在PCCP,14,337 (2012)里，Mayer分析了前面介绍的能量分解方法的缺点，并提出了一个改进的分解方法，可以让原子间能量项更有化学意义，尽管从绝对值上看还是偏离实际键能比较远。具体原理就不谈了，这里直接介绍怎么计算。目前此方法只支持闭壳层HF/DFT波函数。

到http://occam.ttk.hu/programs/下载NEWENPART，然后解压，进入executables文件夹，用文本编辑器将futs文件里的effao-x改为./effao-x，apost4-x改为./apost4-x。然后将要算的fch文件改名为Test.FChk后放到此文件夹里。然后运行诸如./futs ltwd，结果就输出到了ltwd.out里。
注：如果fch文件是在windows下生成的，记得先用dos2unix把文件转换一下

此程序除了输出前面提到和没提到的各种能量分解方式的结果外，最后还有如下输出，对应于PCCP,14,337 (2012)介绍的方法（仍以甲酰胺为例）
       Two-center components in kcal/mol
 
             1  N     2  H     3  H     4  C     5  H     6  O

   1  N      0.00  -240.10  -238.48  -183.34    12.09   116.89
   2  H   -240.10     0.00    21.15    32.27     4.79   -12.65
   3  H   -238.48    21.15     0.00    36.59     4.54   -19.93
   4  C   -183.34    32.27    36.59     0.00  -131.04  -305.70
   5  H     12.09     4.79     4.54  -131.04     0.00    16.10
   6  O    116.89   -12.65   -19.93  -305.70    16.10     0.00
数值越负，代表相应两个原子间键的强度越大。这对于讨论分子在质谱过程的裂解顺序是有帮助的，比Mayer键级更可靠。
这里的结果和上一节的明显不同，这里的结果正确地反映出了C=O键的能量（305.70kcal/mol）是明显大于C-N键的（183.34kcal/mol）。

接着还输出体系中各个原子相对于ROHF下计算的自由原子能量的差值（自由原子能量是在什么基组下得到的不清楚，可能是6-31G**）
   Energy differences with respect to free atoms (kcal/mol)
 
  1  N -121.42    2  H  153.67    3  H  168.02    4  C  430.20    5  H   68.29
  6  O -179.17

虽然不敢说这种改进的能量分解的可靠性有多高，但如果只是讨论原子间相互作用的总能的话，比起前面几节讨论的方法更好。