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

雖然不敢說這種改進的能量分解的可靠性有多高，但如果只是討論原子間相互作用的總能的話，比起前面幾節討論的方法更好。