談談軌道成份的計算方法
注:本文屬于深入、詳細的討論。如果只是著急想馬上得到軌道成分計算結果,直接看本文第6節或Multiwfn手冊的4.8節的例子即可。
談談軌道成份的計算方法
文/Sobereva @北京科音
First release: 2012-Feb-24 Last update: 2020-Jun-3
1 前言
經常有人發帖問怎么計算軌道成份,文獻中也經常給出軌道成份,用于討論電子結構、分析電荷轉移、根據前線軌道理論分析反應位點等等。然而軌道成份的計算方法卻極少被討論,時常看到錯誤的說法在網上以訛傳訛。在《分子軌道成份的計算》(化學學報, 69, 2393)一文中作者深入詳細地討論了軌道成份計算方法的原理并通過實際體系進行了比較,此文可以在此下載:/usr/uploads/file/20150610/20150610022849_33146.pdf。但是對于初學者來說可能一些討論并不很容易理解,而且程序的具體使用步驟未在其中涉及。本文算是那篇文章的補充,本文對各種計算方法的原理介紹比較粗略,忽略了不少細節討論,但是對初學者來說更為易懂,并會詳細說明實際計算時的過程,重點說明了Multiwfn中分析軌道成份功能的基本用法。希望對于懷有“軌道成份怎么計算”問題的人在看完此文以及《分子軌道成份的計算》后都能找到答案。
筆者的另一篇文章《通過軌道離域指數(ODI)衡量軌道的空間離域程度》(http://www.shanxitv.org/525)和此文內容關系較密切,其中筆者提出了一種名為ODI的指數,基于軌道中的原子所占成份來計算,可以很可靠地衡量軌道離域程度,感興趣的讀者建議一看。
2 軌道成份的概念
計算軌道成份也就是計算構成軌道的各個組成部分的貢獻。層次從小到大可以分為:(1)計算基函數的貢獻 (2)計算原子軌道的貢獻 (3)計算原子的貢獻 (4)計算某個片段的貢獻。
正如同蛋糕有不同切法,計算軌道成份的方法也有很多種,沒有一種方法是絕對正確的,畢竟軌道成份不是可觀測量,定義完全是人為的。很多文獻中給出軌道成份時不注明是怎么計算的,這是不應該的,這導致數據難以被后人重現,也沒法得知作者所用的軌道成份計算方法及結果是否合理。
不管什么成份計算方法,最基本的要求是能夠歸一,即所有部分的貢獻和為100%。另一個基本要求是成份必須是正值,負值是沒有意義的。在滿足這兩個要求下,還應當盡量滿足其它一些要求:(1)滿足旋轉不變性。因為分子的朝向是任意的,不能因為分子在空間中旋轉了,成份的計算結果產生過大變化。(2)結果應當能隨基組增大而收斂,否則相當于沒法得到一個確定的結果 (3)計算速度應當快,算法應當簡單,易于實現 (4)物理意義應當清晰 (5)結果應當和分子軌道圖形定性上相符 (6)適用體系應當廣泛 (7)數值穩定性好,避免出現異常值。
接下來談談計算基函數、原子軌道、原子以及片段的成份的基本原理和具體計算步驟。本文所指的“軌道”可以是任何類型軌道,比如分子軌道、定域化軌道、自然軌道等等,但舉例都用分子軌道舉例,對于其它類型軌道的處理方法是一樣的。
3 計算基函數所占的成份
基函數本身沒有直接的化學意義,它只是求解量子化學問題的數學工具。之所以要提到它是因為它可以認為是構成軌道的最小單元,在計算軌道中其它組成部分的成份時會涉及到它。另外,對于諸如STO-3G這種極小基來說,一個基函數就對應于一個原子軌道,計算基函數的成份此時就是計算原子軌道的成份,是有化學意義的。
這里先來考慮一個最簡單的異核雙原子體系HeH+,用Gaussian在HF/STO-3G計算后的pop=full輸出是
Molecular Orbital Coefficients:
1 2
O V
Eigenvalues -- -1.52378 -0.26764
1 1 He 1S 0.90013 -0.64981
2 2 H 1S 0.19438 1.09302
可知占據的那條分子軌道波函數可寫為ψ1=0.9*χ_1 + 0.1944*χ_2。其中χ是基函數符號。根據軌道的歸一化條件可寫為
1=<ψ1|ψ1>=0.9^2*<χ_1|χ_1>+2*0.9*0.1944*<χ_1|χ_2>+0.1944^2*<χ_2|χ_2>=0.9^2*S11+2*0.9*0.1944*S12+0.1944^2*S22。
上式中,0.9^2和0.1944^2分別是χ_1和χ_2的“定域項”,表明了它們獨自對ψ1的貢獻。2*0.9*0.1944*S12是交叉項,代表了χ_1和χ_2對ψ1的聯合貢獻。
由于基函數是歸一化了的,所以S11和S22都為1,因此<ψ1|ψ1>=0.8478+2*0.9*0.1944*S12。很顯然,很多人用軌道系數的平方來計算軌道成份是大錯特錯,0.9^2+0.1944^2等于0.8478而明顯不得1!即這么算出來的軌道成份和不是100%。這是因為基函數間不是正交的,重疊積分S12不為零,所以計算軌道成份必須考慮怎么把交叉項劃分給各個基函數。
計算基函數成份的方法有三種,分別是Mulliken、Stout-Politzer和SCPA方法。本質上,這三種方式的差異僅在于怎么劃分交叉項。
Mulliken方法在計算χ_n對軌道的貢獻時,將χ_n的定域項全部算作χ_n的貢獻,而將χ_n與其它基函數間的交叉項的一半劃給χ_n。對于上例,He對軌道的貢獻是(0.9^2+0.9*0.1944*S12)*100%。利用iop(3/33=1)關鍵詞可輸出重疊矩陣可知S12=0.4343,因此He的貢獻為88.6%。實際上,計算軌道成份的Mulliken方法就對應于Mulliken布居用到的方法。Mulliken布居方法中,某條軌道上某基函數的布居數就是這個基函數的貢獻值乘上軌道占據數,因此χ_1在ψ1這條雙占據軌道上的Mulliken布居數就是88.6%*2=1.772。由于以He為中心的只有χ_1這么一個基函數,而且此體系只有ψ1是占據軌道,因此He原子的Mulliken布居數就是1.772。由于He的核電荷為2,因此它的Mulliken電荷是2-1.772=0.228,和Gaussian輸出的一致。
Mulliken方法對交叉項的平分做法被一些人批評,因為這顯得沒物理意義。于是有些人提出了改進方法。Stout-Politzer方法考慮了基函數的不均等性,它將交叉項按照相應兩個基函數的系數平方比來劃分。對于上例,χ_1所劃得的交叉項就是2*0.9*0.1944*S12* 0.9^2/(0.9^2+0.1944^2)。雖然Stout-Politzer方法看起來比Mulliken方法考慮得更周到,但是實際結果反倒往往有所惡化,數值不穩定性增加,出現負值的傾向加劇,還破壞了Mulliken方法的旋轉不變性。此方法不建議使用。
很多人在用系數平方來計算軌道成份時也注意到了必須滿足歸一條件,于是就乘了個歸一化因子,因此χ_1的貢獻就是0.9^2/(0.9^2+0.1944^2)*100%。這個做法實際上對應于SCPA布居方法所用到的軌道劃分方法。這個方法的計算公式經過簡單推導(見《分子軌道成份的計算》),可知本質上和Mulliken、Stout-Politzer方法的差異也僅是劃分交叉項的做法不同而已。Mulliken和Stout-Politzer時常出現沒意義的負的成份值,對于虛軌道尤為嚴重,而SCPA方法的好處是結果一定為正。另外,SCPA方法計算簡單,不需要考慮重疊積分,通過比較發現數值合理性也比Mulliken和Stout-Politzer方法都略強,是推薦使用的計算基函數貢獻的方法。
4 計算原子軌道所占的成份
在分子體系中實際上沒有清晰的原子軌道的概念,但是以某些方法將原子軌道概念還原出來,并計算其在各軌道中的成份,有助于分析電子結構特征、成鍵本質。
4.1 基函數加和方法
第一種計算方法是利用基組與原子軌道的對應關系。對于極小基,每個基函數描述一個原子軌道,因此就按上一節討論的方法就得到了原子軌道的貢獻。而對于擴展基,則是以一個或多個基函數來描述一個原子軌道體現的效應,因此可以把原子軌道對應的基函數的貢獻相加作為原子軌道的貢獻。
例如對于HF/6-31G*下計算水分子,Gaussian輸出的占據軌道系數是這樣的:
1 2 3 4 5
(A1)--O (A1)--O (B2)--O (A1)--O (B1)--O
Eigenvalues -- -20.55790 -1.34610 -0.71418 -0.57083 -0.49821
1 1 O 1S 0.99462 -0.20953 0.00000 -0.07310 0.00000
2 2S 0.02117 0.47576 0.00000 0.16367 0.00000
3 2PX 0.00000 0.00000 0.00000 0.00000 0.63927
4 2PY 0.00000 0.00000 0.50891 0.00000 0.00000
5 2PZ -0.00134 -0.09475 0.00000 0.55774 0.00000
6 3S 0.00415 0.43535 0.00000 0.32546 0.00000
7 3PX 0.00000 0.00000 0.00000 0.00000 0.51183
8 3PY 0.00000 0.00000 0.30390 0.00000 0.00000
9 3PZ 0.00046 -0.04982 0.00000 0.40482 0.00000
10 4XX -0.00394 -0.00103 0.00000 0.01187 0.00000
11 4YY -0.00421 0.02692 0.00000 0.00079 0.00000
12 4ZZ -0.00409 0.02132 0.00000 -0.04630 0.00000
13 4XY 0.00000 0.00000 0.00000 0.00000 0.00000
14 4XZ 0.00000 0.00000 0.00000 0.00000 -0.03417
15 4YZ 0.00000 0.00000 -0.05088 0.00000 0.00000
16 2 H 1S 0.00032 0.13302 0.23243 -0.14007 0.00000
17 2S -0.00021 0.00173 0.10729 -0.08280 0.00000
18 3 H 1S 0.00032 0.13302 -0.23243 -0.14007 0.00000
19 2S -0.00021 0.00173 -0.10729 -0.08280 0.00000
注意這樣的表里每行對應一個基函數,而基函數前的數字編號沒有物理意義,它只是為了區分以相同原子為中心的類型相同的基函數,編號越靠前的通常分布區域越接近原子核。這些行并非直接對應原子軌道,數字也不是原子軌道的主量子數。很多初學者容易在這點上混淆。
對于6-31G*基組,每個內層原子軌道用一個收縮度為6的基函數來體現,故這個基函數的貢獻就是相應內層原子軌道的貢獻。也就是說,氧的1s原子軌道成份就是這里的1S基函數的成份(本文中原子軌道符號用小寫字母,基函數符號用大寫字母)。而6-31G*中每個價層原子軌道用一個收縮度為3和一個收縮度為1的相應類型的基函數一起來表現,這兩個基函數的貢獻之和可認為是對應的原子軌道的貢獻。也就是說,氧的2px原子軌道的貢獻可認為是這里2PX和3PX基函數的貢獻之和,而氫的1s原子軌道的貢獻就是其1S和2S基函數貢獻之和。表中的4XX、4YY等等這些D型基函數沒有物理意義,顯然不可能對應于本來就不存在的2d原子軌道,這些D型基函數存在的意義只是為了在自洽場計算中更充分地描述出實際的電子密度的極化效果
按照基函數與原子軌道對應關系來計算原子軌道貢獻的方法有很大局限性。有些基組基函數很多,不容易搞清楚對應關系(但可以巧妙利用Multiwfn的布居分析獲得對應關系,見《利用布居分析判斷基函數與原子軌道的對應關系》http://www.shanxitv.org/418)。對于彌散、極化函數也沒法處理,考慮它們的話沒法確定應該將它們歸在哪個原子軌道里,而不考慮的話則原子軌道貢獻和不為100%(但可以對已考慮了的原子軌道貢獻值做歸一化處理)。而且由于基函數之間的混合、缺乏對基函數在原子軌道中權重的考慮,導致這種方法在原理上不甚嚴格。雖然缺點能列舉出很多,但實際應用中這種計算方法還是比較便利的,對于大多數情況來說沒什么問題。但并不建議在含有彌散基函數的情況下,或者基組不平衡較強(即某些原子的基組質量和相鄰原子基組質量相差懸殊。如某原子用6-311G(2d),而旁邊的原子只用3-21G)時使用本方法。由于用這種方法計算原子軌道貢獻時先要計算基函數的貢獻,這就牽扯到用什么方法來計算基函數的貢獻,我仍建議用SCPA方法。
用本節的方法對于計算d及以上角動量原子軌道的貢獻時有一點需要注意,也就是應當盡量用球諧型高斯函數,而不建議用笛卡爾型高斯函數。這點不打算在這里多談,請參閱《談談5d、6d型d殼層基函數與它們在Gaussian中的標識》一文(http://www.shanxitv.org/51)以及《分子軌道成份的計算》中的討論。
4.2 NAO方法
另一種計算原子軌道成份的方法是借助于NBO分析。NBO分析可以從體系的密度矩陣信息中得到自然原子軌道(NAO),NAO的數目與原先基函數相同。NAO當中被稱為“極小集”的那部分就對應于一般意義上的原子軌道。并且NBO分析還可以將分子軌道以NAO來展開。由于NAO彼此間都是正交的,不用考慮怎么劃分交叉項,因此NAO在相應軌道中的系數平方就是它在軌道中的成份了。
在Gaussian中,用了pop=nboread,輸入文件末尾空一行寫上$NBO NAOMO $END,在NBO模塊中就會輸出以NAO為基的MO的展開系數,還是以HF/6-31G*的水為例:
MOs in the NAO basis:
NAO 1 2 3 4 5 6 7 8
---------- ------- ------- ------- ------- ------- ------- ------- -------
1. O 1 (S) 0.9958 -0.0896 0.0000 -0.0177 0.0000 0.0046 0.0000 0.0000
2. O 1 (S) 0.0847 0.8834 0.0000 0.2928 0.0000 0.2288 0.0000 0.0000
3. O 1 (S) 0.0003 0.0081 0.0000 -0.0246 0.0000 -0.2297 0.0000 0.0000
4. O 1 (S) 0.0000 0.0000 0.0000 0.0000 0.0000 0.0554 0.0000 0.0000
5. O 1 (px) 0.0000 0.0000 0.0000 0.0000 0.9993 0.0000 0.0000 0.0000
6. O 1 (px) 0.0000 0.0000 0.0000 0.0000 -0.0173 0.0000 0.0000 0.0000
7. O 1 (py) 0.0000 0.0000 0.8554 0.0000 0.0000 0.0000 -0.3288 -0.3600
8. O 1 (py) 0.0000 0.0000 0.0334 0.0000 0.0000 0.0000 0.0250 0.4768
9. O 1 (pz) 0.0023 -0.1585 0.0000 0.9170 0.0000 -0.2598 0.0000 0.0000
10. O 1 (pz) -0.0011 -0.0006 0.0000 -0.0124 0.0000 0.0254 0.0000 0.0000
11. O 1 (d1) 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
12. O 1 (d2) 0.0000 0.0000 0.0000 0.0000 -0.0342 0.0000 0.0000 0.0000
13. O 1 (d3) 0.0000 0.0000 -0.0388 0.0000 0.0000 0.0000 0.0177 -0.1517
14. O 1 (d4) 0.0001 -0.0176 0.0000 0.0075 0.0000 0.0212 0.0000 0.0000
15. O 1 (d5) 0.0000 0.0057 0.0000 -0.0352 0.0000 0.0106 0.0000 0.0000
16. H 2 (S) 0.0244 0.3046 0.3644 -0.1882 0.0000 -0.4212 0.3712 0.4027
17. H 2 (S) 0.0028 -0.0165 0.0095 -0.0079 0.0000 -0.4839 0.5547 -0.3845
18. H 3 (S) 0.0244 0.3046 -0.3644 -0.1882 0.0000 -0.4212 -0.3712 -0.4027
19. H 3 (S) 0.0028 -0.0165 -0.0095 -0.0079 0.0000 -0.4839 -0.5547 0.3845
下面是在一開始輸出的NAO布居分析
NAO Atom No lang Type(AO) Occupancy Energy
----------------------------------------------------------
1 O 1 S Cor( 1S) 1.99992 -20.39645
2 O 1 S Val( 2S) 1.74644 -1.14691
3 O 1 S Ryd( 3S) 0.00135 1.42530
4 O 1 S Ryd( 4S) 0.00000 3.93946
5 O 1 px Val( 2p) 1.99707 -0.49473
6 O 1 px Ryd( 3p) 0.00060 1.16772
7 O 1 py Val( 2p) 1.46327 -0.31260
8 O 1 py Ryd( 3p) 0.00223 1.30664
9 O 1 pz Val( 2p) 1.73207 -0.41806
10 O 1 pz Ryd( 3p) 0.00031 1.19893
11 O 1 dxy Ryd( 3d) 0.00000 2.03061
12 O 1 dxz Ryd( 3d) 0.00234 2.06405
13 O 1 dyz Ryd( 3d) 0.00301 2.91411
14 O 1 dx2y2 Ryd( 3d) 0.00073 2.47969
15 O 1 dz2 Ryd( 3d) 0.00254 2.02361
16 H 2 S Val( 1S) 0.52321 0.33308
17 H 2 S Ryd( 2S) 0.00086 0.70497
18 H 3 S Val( 1S) 0.52321 0.33308
19 H 3 S Ryd( 2S) 0.00086 0.70497
表中Cor對應于內核NAO,Val對應于價層NAO,標著這兩類標簽的NAO就是極小集部分,可以簡單認為它們就是一般意義的原子軌道。而Ryd型NAO是用來描述那些沒有被Cor和Val型NAO所描述的信息,它們對占據的分子軌道貢獻極小而一般可以忽略,但是對于高能態的虛軌道則可能有很大貢獻。通常我們分析的都是占據軌道和能量最低的幾個虛軌道,所以一般不必考慮Ryd型NAO。
將上面兩個表對照,就能知道怎么計算各個原子軌道成份了。例如,要計算第三個MO中氧的2py軌道的貢獻,從上面表中看到py Val( 2p)對應的NAO標號是7,因此就將系數表中第三列第七行的數值0.8554取平方乘100%即可,即73.2%。
4.3 AOIM方法
AOIM方法有兩個版本,第一個版本是劉文劍和黎樂民提出的,而這里涉及的是后來由Haigang Lu、黎樂民等人提出的版本。相應的分析程序AOIM可以在http://faculty.sxu.cn/luhg/aoim.html下載(已失效)。
AOIM是通過正交投影的方法在分子體系中還原出原子軌道,每個原子軌道以一個STO函數表達,這個過程中令誤差函數最小化來盡可能地避免波函數信息的損失。同時也得到了每個分子軌道向各個原子軌道的展開系數。AOIM相當于將原始的在擴展基下表達的波函數轉化為了在極小基下表達,因此按照本文第三節的方法,就可以計算各軌道中每個原子軌道的成份了,在AOIM中這一步用的是Mulliken方法。
AOIM程序需要讀入Gaussian的fch文件來進行分析。在DOS下執行比如aoim < ltwd.fch就可以對ltwd.fch進行AOIM分析。還是那個水分子的例子,AOIM輸出的軌道成份信息部分是
MO components in correctted STOs(>0.05):
MO: 1 E= -20.55786 A.U.
1.00 1 O1S
MO: 2 E= -1.34606 A.U.
0.81 1 O2S
0.08 2 H1S
0.08 3 H1S
MO: 3 E= -0.71414 A.U.
0.60 1 O2Py
0.20 2 H1S
0.20 3 H1S
MO: 4 E= -0.57083 A.U.
0.09 1 O2S
0.81 1 O2Pz
...(略)
可見第三條MO的氧的2py軌道的貢獻是60%。和之前用NAO方法的結果73.2%有一定差異,這是正常的,不同計算方法不可能總在定量上很一致,尤其是軌道涉及多個原子間交疊的情況差異會更明顯,但多數情況下本節介紹的三種方法結果還是定性一致的。
AOIM程序只會輸出貢獻值大于5%的原子軌道,輸出的分子軌道只包含占據軌道和LUMO及LUMO+1軌道,如果想考察貢獻值更低的基函數以及更高階的虛軌道,可以自行修改代碼,這很容易做到。
NAO和AOIM方法的基組穩定性都明顯好于加和基函數貢獻的方法,能隨著基組質量增加而很快收斂,而且原理都比較嚴格。可惜AOIM程序年久失修,對Gaussian09的fch文件不兼容,對大體系分析速度慢,不支持f及以上的角動量基函數,也不支持球諧型基函數,且在個別體系比如氧化鈹上有嚴重問題。對于計算原子軌道成份的方法,整體來說我最推薦NAO方法。
5 計算原子以及片段所占的成份
如果只是想定性了解各個原子對軌道貢獻的相對大小,最簡單的辦法就是顯示軌道圖形,將isovalue調到一個能夠區分不同位置軌道波函數數值大小的值,然后圖形觀察各個原子附近的等值面涵蓋范圍的大小即可,范圍越大表明貢獻越大。
而定量計算原子所占成份最直接的辦法就是將原子所擁有的全部原子軌道所占成份加和,這里我也是建議用NAO方法來計算原子軌道所占成份。
還有一類基于實空間劃分的方法可用于計算原子所占成份。這類方法都是將整個分子空間劃分為一個個原子區域,然后積分某個軌道在相應原子空間內的模的平方來作為此原子在此軌道中的成份。一種廣為熟知的劃分方法就是Bader提出的QTAIM劃分方法,是根據電子密度零通量面來劃分,但這種劃分一方面是積分困難,另一方面也沒有很強的化學意義。比較適合用于計算原子貢獻的劃分方式是Hirshfeld劃分,這種劃分是模糊的,原子的空間是相互交疊的,沒有明確的邊界,具體來說就是空間中每個位置上各個原子都占有一定權重。設ρ_pro(r)為分子中每個原子都處在自由狀態時電子密度的疊加(也叫Promolecule密度),r代表空間坐標,而ρ_A(r)是A原子在自由狀態下的密度,那么A原子的權重函數就是ρ_A(r)/ρ_pro(r)。在離A原子核越近的區域中A原子的權重會越大。這種Hirshfeld劃分的化學意義比較直觀明確,積分也比較容易。在計算結果上與NAO方法定性一致,但計算量還是要大于NAO方法。對于一些虛軌道,尤其是高能態時,Ryd型NAO所占成份不可忽視,也就是說加和極小集NAO的成份會明顯偏離100%,此時就不建議用NAO方法了,而建議用Hirshfeld方法。
利用Becke函數也可以進行模糊式劃分,形式上和Hirshfeld劃分類似,結果也比較相似。相對于Hirshfeld劃分的好處是原子權重不需要通過原子的密度來構建,而只需要提供原子半徑即可。Hirshfeld還有種重要的變體叫Hirshfeld-I,它通過迭代方式逐漸調整原子權重函數,使之可以響應原子實際所處的化學環境,基于Hirshfeld-I劃分也可以計算原子對軌道的貢獻。
計算分子中某個片段的方法沒什么可說的,就是將其中原子所占成份加和即可。
6 在Multiwfn中計算軌道成份
在Multiwfn(http://www.shanxitv.org/multiwfn)中的主功能8里提供了分析軌道成份的功能,支持Mulliken、Stout-Politzer、SCPA方法,還獨家支持Hirshfeld、Hirshfeld-I、Becke、NAO、AIM方法。功能強大、操作簡單,輸出信息清晰易讀,比其它程序好用得多。
使用Mulliken、Stout-Politzer、SCPA方法進行分析的話必須用包含基函數信息的文件作為輸入文件,如.mwfn、.fch、.molden、GAMESS-US輸出文件等,不清楚的話詳見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》(http://www.shanxitv.org/379)。以Becke、Hirshfeld、Hirshfeld-I、AIM方法分析原子貢獻的話任何帶有波函數信息的文件都可以用(.fch/.molden/.wfn/.wfx等)。以NAO方法分析的話需要提供NBO程序的輸出文件。這里還是以分析HF/6-31G*分子軌道成份為例進行說明。選項的詳細介紹請閱讀手冊的3.10節。
6.1 在Multiwfn中使用Mulliken、Stout-Politzer、SCPA方法進行分析
啟動Multiwfn,輸入.fch或.molden等含有基函數信息的文件路徑后,選8就進入了成份分析界面。如果你的經驗不豐富的話,可以先選-1進入分子片段設定界面,然后輸入all來查看一下當前體系中的所有基函數信息(之后輸入q退回):
Basis: 1 Shell: 1 Center: 1(O ) Type:S
Basis: 2 Shell: 2 Center: 1(O ) Type:S
Basis: 3 Shell: 3 Center: 1(O ) Type:X
Basis: 4 Shell: 3 Center: 1(O ) Type:Y
Basis: 5 Shell: 3 Center: 1(O ) Type:Z
Basis: 6 Shell: 4 Center: 1(O ) Type:S
Basis: 7 Shell: 5 Center: 1(O ) Type:X
Basis: 8 Shell: 5 Center: 1(O ) Type:Y
Basis: 9 Shell: 5 Center: 1(O ) Type:Z
Basis: 10 Shell: 6 Center: 1(O ) Type:XX
Basis: 11 Shell: 6 Center: 1(O ) Type:YY
Basis: 12 Shell: 6 Center: 1(O ) Type:ZZ
Basis: 13 Shell: 6 Center: 1(O ) Type:XY
Basis: 14 Shell: 6 Center: 1(O ) Type:XZ
Basis: 15 Shell: 6 Center: 1(O ) Type:YZ
Basis: 16 Shell: 7 Center: 2(H ) Type:S
Basis: 17 Shell: 8 Center: 2(H ) Type:S
Basis: 18 Shell: 9 Center: 3(H ) Type:S
Basis: 19 Shell: 10 Center: 3(H ) Type:S
每個基函數對應的殼層,所屬中心,基函數類型都一目了然。每個殼層包含一組基函數,即S殼層包含一個S基函數,P殼層包含X,Y,Z三個基函數,6D型殼層包含XX,YY,ZZ,XY,XZ,YZ共六個基函數...
比如想用Mulliken方法分析第3條分子軌道的成份,就選擇1,然后輸入軌道編號3,就得到了以下結果:
Threshold of absolute value: > 0.500000%
Orbital: 3 Energy(a.u.): -0.71414263 Occ: 2.00000000 Type: Alpha&Beta
Basis Type Atom Shell Local Cross term Total
4 Y 1(O ) 3 25.89789% 14.57170% 40.46960%
8 Y 1(O ) 5 9.23602% 16.82836% 26.06437%
15 YZ 1(O ) 6 0.25896% 0.83797% 1.09693%
16 S 2(H ) 7 5.40298% 7.23541% 12.63839%
17 S 2(H ) 8 1.15080% 2.39536% 3.54616%
18 S 3(H ) 9 5.40298% 7.23541% 12.63839%
19 S 3(H ) 10 1.15080% 2.39536% 3.54616%
Sum up listed above: 48.50044% 51.49956% 100.00000%
Sum up all basis functions: 48.50044% 51.49956% 100.00000%
Composition of each shell, threshold of absolute value: > 0.500000%
Shell 3 Type: P in atom 1(O ) : 40.46960%
Shell 5 Type: P in atom 1(O ) : 26.06437%
Shell 6 Type: D in atom 1(O ) : 1.09693%
Shell 7 Type: S in atom 2(H ) : 12.63839%
Shell 8 Type: S in atom 2(H ) : 3.54616%
Shell 9 Type: S in atom 3(H ) : 12.63839%
Shell 10 Type: S in atom 3(H ) : 3.54616%
Composition of each atoms:
Atom 1(O ) : 67.630897%
Atom 2(H ) : 16.184551%
Atom 3(H ) : 16.184551%
默認情況下只有貢獻值大于0.5%的項被輸出,這個閾值由Multiwfn目錄下的settings.ini里的compthres參數控制。由于這條軌道中各個基函數的貢獻要么為0要么就大于0.5%,因此Sum up listed above(顯示出的基函數的貢獻和)的值與Sum up all basis functions(全部基函數的貢獻和)相同。
Total項是Local和Cross term的總和。Local就是之前的說的定域項的貢獻,Cross term是指這個基函數與其它所有基函數間交叉項中歸屬于這個基函數的貢獻。
Multiwfn還輸出了每個殼層對軌道的貢獻,即殼層內基函數的貢獻之和。注意除非分子的朝向已經清楚指定了,才可以在文中討論某個基函數的貢獻,否則討論基函數的貢獻其實是沒有意義的,不少文獻都犯了這個錯誤。因為分子在旋轉后,各個基函數的貢獻也會發生變化(除了S型以外,因為它是球對稱的)。在不指明朝向時,只能對某個殼層的貢獻進行討論,因為Mulliken方法得出的殼層的貢獻值是不會因為分子旋轉而改變的。對于原子軌道的貢獻來說也是一樣,必須指明了朝向才能說某個原子軌道的貢獻,否則只能說原子軌道殼層的貢獻。注意在Stout-Politzer方法下,即便是殼層的貢獻值也存在旋轉依賴性。而SCPA方法下,必須使用的是球諧型高斯函數才能嚴格滿足旋轉不變性,不過使用笛卡爾型高斯函數時旋轉依賴性并不大。這個問題在這里不做繼續探討,可參見《分子軌道成份的計算》的討論。
由于4、8號基函數對應于氧的2py軌道,因此40.46960%+26.06437%=66.53%就是2py對第三條分子軌道的貢獻,和前文用AOIM和NAO方法算出來的60%和73.2%都在比較相近。由于第三個和第五個殼層的貢獻分別和第4號、第8號基函數貢獻相同,也就是說在此MO當中氧的2px和2pz的成份都為0。
輸出信息最后給出了通過基函數貢獻加和得到的每個原子的貢獻。
輸入0退回上一級菜單。下面我們要進行片段分析。Multiwfn可以十分靈活、方便地自定義片段,片段不僅可以包含幾個原子,也可以將特定的基函數、殼層、原子混合在一個片段里,片段的貢獻值就是通過其中包含的基函數的貢獻加和得到。這里我們先把氧的2p原子軌道殼層定義為第一個片段。選擇-1進入第一個片段的定義界面,從all命令給出的列表以及6-31G*的基組定義中可以得知第三個和第五個殼層對應于氧的2p原子軌道殼層,于是就輸入s 3,5。此時如果輸入list命令來查看當前片段內已有的基函數列表,就會看到這兩個殼層內包含的3,4,5,7,8,9號基函數都已經被加入了。輸入q保存并退出片段定義界面。
然后選-2進入第二個片段定義界面。這里將兩個氫原子(原子編號為2和3)定義為第二個片段,于是就輸入a 2,3。輸入q保存并退回。
選擇選項4,Multiwfn就會使用Mulliken方法將第一個片段對所有軌道的貢獻一次性都輸出出來:
Orb# Type Energy Occ c^2 Int.cross Ext.cross Total
1 AB -20.558 2.000 0.000% 0.000% 0.000% 0.000%
2 AB -1.346 2.000 1.146% 0.474% 0.926% 2.546%
3 AB -0.714 2.000 35.134% 15.513% 15.887% 66.534%
4 AB -0.571 2.000 47.503% 22.651% 9.568% 79.722%
5 AB -0.498 2.000 67.063% 32.820% 0.000% 99.883%
6 AB 0.213 0.000 29.745% 10.680% -36.963% 3.462%
7 AB 0.307 0.000 81.284% 27.464% -105.805% 2.942%
8 AB 1.032 0.000 50.954% 6.230% -35.251% 21.933%
9 AB 1.133 0.000 79.565% -35.898% -3.568% 40.099%
10 AB 1.168 0.000 200.037% -100.056% 0.000% 99.980%
11 AB 1.178 0.000 63.463% -23.135% -4.472% 35.855%
12 AB 1.385 0.000 345.645% -160.782% -76.700% 108.163%
13 AB 1.431 0.000 159.364% -58.247% -66.716% 34.400%
14 AB 2.021 0.000 1.705% -0.037% -1.019% 0.650%
15 AB 2.031 0.000 0.000% 0.000% 0.000% 0.000%
16 AB 2.067 0.000 0.110% 0.027% 0.000% 0.136%
17 AB 2.636 0.000 55.219% -2.688% -50.347% 2.184%
18 AB 2.966 0.000 80.748% -0.574% -79.747% 0.427%
19 AB 3.978 0.000 11.297% -3.714% -6.502% 1.082%
片段1的貢獻值Total是由c^2、Int.cross、Ext.cross三項加和得到的。c^2指的是片段內基函數的系數平方和。Int.cross(內交叉項)指的是片段1包含的基函數之間的交叉項之和。Ext.cross(外交叉項)指的是片段1的基函數與體系中其余基函數間的交叉項中被劃歸到片段1的部分。可以看到虛軌道當中出現了大量負值,這是Mulliken方法的弊病之一,因此如果需要分析虛軌道,建議用保證結果為正的SCPA,但此時將不會輸出交叉項信息。軌道的Type顯示的AB代表這個軌道可以同時占alpha電子和Beta電子,也即閉殼層軌道。
由于已經定義了片段2,所以下面還會輸出片段1、2間的交叉項(如果只定義了片段1,則只輸出上面那些內容)
Cross term between fragment 1 and 2 and their individual parts:
Orb# Type Energy Occ Frag1 part Frag2 part Total
1 AB -20.558 2.000 0.0000% 0.0000% 0.0000%
2 AB -1.346 2.000 0.9263% 0.9263% 1.8525%
3 AB -0.714 2.000 15.8871% 15.8871% 31.7743%
4 AB -0.571 2.000 9.5681% 9.5681% 19.1363%
5 AB -0.498 2.000 0.0000% 0.0000% 0.0000%
6 AB 0.213 0.000 -36.9633% -36.9633% -73.9266%
7 AB 0.307 0.000 -105.8051% -105.8051% -211.6102%
8 AB 1.032 0.000 -35.2511% -35.2511% -70.5022%
9 AB 1.133 0.000 -3.5678% -3.5678% -7.1357%
10 AB 1.168 0.000 0.0000% 0.0000% 0.0000%
11 AB 1.178 0.000 -4.4724% -4.4724% -8.9447%
12 AB 1.385 0.000 -76.6996% -76.6996% -153.3991%
13 AB 1.431 0.000 -66.7162% -66.7162% -133.4324%
14 AB 2.021 0.000 -1.0185% -1.0185% -2.0370%
15 AB 2.031 0.000 0.0000% 0.0000% 0.0000%
16 AB 2.067 0.000 0.0000% 0.0000% 0.0000%
17 AB 2.636 0.000 -50.3467% -50.3467% -100.6933%
18 AB 2.966 0.000 -79.7466% -79.7466% -159.4933%
19 AB 3.978 0.000 -6.5016% -6.5016% -13.0032%
由于Mulliken方法是將交叉項平分,所以片段1和片段2所占的交叉項是相同的。兩個片段間總交叉項描述了兩個片段間的相互作用程度,這與Mulliken鍵級直接相關。比如第三條MO,0.3177乘上軌道的占據數2得到的0.6354實際上就是此MO對片段1、2間的Mulliken鍵級的貢獻。
6.2 在Multiwfn中使用Hirshfeld方法分析原子的貢獻
下面介紹一下用Hirshfeld方法計算原子的貢獻。生成各個原子的Hirshfeld權重時要用到原子在自由狀態下的密度。Multiwfn提供了兩種方法計算:(1)直接用程序內置的原子密度,詳見手冊附錄3,這種做法對用戶來說最方便 (2)基于原子波函數文件計算原子密度,詳見手冊3.7.3節或《使用Multiwfn作電子密度差圖(http://www.shanxitv.org/113)的第6~8節。兩種做法結果雖然有異,但差別很小。為了簡便,本文例子使用方法(1)。
在軌道成份分析功能里選8進入Hirshfeld分析,然后選1用程序內置的原子密度,之后Multiwfn會先進行初始化。然后比如要分析第3條分子軌道就輸入3,將看到以下結果:
After normalization:
Atom 1(O ) : 64.926881%
Atom 2(H ) : 17.536559%
Atom 3(H ) : 17.536559%
由于數值積分精度原因,直接算出來的結果的加和值可能偏離100%一些,所以程序會自動做一下歸一化。
如果你想獲得某個原子在指定分子軌道范圍內的貢獻的話,就輸入-2,之后輸入原子編號和軌道范圍即可,比如輸入1和1,5就得到了氧原子在1至5號MO中的貢獻:
Orb# Type Energy Occ Composition Population
1 Alpha&Beta -20.558 2.000 99.676867% 1.993537
2 Alpha&Beta -1.346 2.000 80.133155% 1.602663
3 Alpha&Beta -0.714 2.000 64.926646% 1.298533
4 Alpha&Beta -0.571 2.000 82.437155% 1.648743
5 Alpha&Beta -0.498 2.000 88.761277% 1.775226
Population of this atom in these orbitals: 8.318702
這里還給出了Population值,也就是原子對軌道的貢獻值乘以軌道的占據數。最后給出了它在這些軌道的布居數加和值。由于1至5號就是全部的占據軌道,因此氧原子在水分子中的總電子數就是8.318702,也因此它的Hirshfeld電荷就是8-8.318702=-0.3187。
還可以選-9來定義片段,然后比如輸入2,3就把兩個氫都加進片段里了。之后輸入軌道序號后,就會看到片段的貢獻(Fragment contribution)。也可以選-3,然后輸入軌道范圍,來得到此片段對這些軌道各自的貢獻。
6.3 在Multiwfn中使用Becke、Hirshfeld-I、AIM方法分析原子的貢獻
在軌道成份分析功能里選9就能進入Becke分析。由于使用Becke方法和使用Hirshfeld方法分析原子的貢獻的操作完全相同,因此這里不再舉例。一般而言,這兩種方法所得結果是比較接近的,用哪個都可以。對于大體系的話Hirshfeld方法更便宜,而且由于Hirshfeld方法的物理意義更強,我更傾向于用Hirshfeld方法。
使用Multiwfn基于Hirshfeld-I劃分計算軌道成份在本文也不舉例了,因為這種做法雖然比Hirshfeld方法更嚴格,但需要額外多花不少時間來通過迭代改進原子空間,而從結果來看也沒顯著優勢。如果你對結果的物理意義要求特別高的話倒是可以考慮用這個。一般的用法是把Multiwfn自帶的examples目錄下的atmrad目錄拷到當前目錄下,然后進入軌道成分分析模塊后,選擇Hirshfeld-I方法,然后選相應選項開始計算即可。
Multiwfn還支持基于AIM劃分計算軌道中原子的成份,由于不常用,這里也不給出例子了。這種做法計算的耗時很高,因為產生AIM盆的步驟非常昂貴,而結果并不比基于Hirshfeld或Becke方法算的好。如果確實想用這種做法,看Multiwfn手冊4.8.6節。
6.4 在Multiwfn中使用NAO方法分析原子軌道、殼層、原子的貢獻
如4.2節的例子所示,計算原子軌道的貢獻就是將以NAO為基的分子軌道系數矩陣的相應矩陣元取平方乘100%。Multiwfn自身并不會生成這個矩陣,用戶需要將含有這個矩陣的NBO程序的輸出文件作為Multiwfn的輸入文件。對于Gaussian,將route section寫上pop=nboread,將分子坐標末尾空一行寫上$NBO NAOMO $END,則這個Gaussian輸出文件就能作為Multiwfn進行NAO分析的輸入文件了(不要誤用.fch文件作為輸入文件!)。如果你用的是獨立的NBO程序,即GENNBO,則應在.47文件的$NBO和$END中間添上NAOMO關鍵詞,產生的輸出文件也將可以作為Multiwfn的輸入文件。
進入Multiwfn后選8,再選7,就進入了NAO分析模塊。選0再輸入分子軌道序號就能輸出這個分子軌道的詳細成份信息。輸出哪些項可以通過2 Select output mode選項來控制,默認是只輸出貢獻大于0.5%的Cor和Val型的NAO和殼層。比如分析第3條軌道會輸出:
NAO# Center Label Type Composition
7 1(O ) py Val( 2p) 73.154%
16 2(H ) S Val( 1S) 13.279%
18 3(H ) S Val( 1S) 13.279%
Condensed NAO terms to shells:
Atom: 1(O ) Shell: 5( 2p Val) 73.154%
Atom: 2(H ) Shell: 8( 1S Val) 13.279%
Atom: 3(H ) Shell: 10( 1S Val) 13.279%
Condensed NAO terms to atoms:
Center Composition
1(O ) 73.416%
2(H ) 13.288%
3(H ) 13.288%
Core composition: 0.000%
Valence composition: 99.711%
Rydberg composition: 0.280%
輸出信息簡單易懂。比如可見主要貢獻來自于O的2p殼層,更具體來說是2py。當前軌道的Rydberg composition只有0.28%,比較小,說明對這個軌道做NAO軌道成份分析是有意義的。
在NAO分析模塊中也同樣可以定義片段來計算它對一批軌道的貢獻,我們這里計算兩個氫上的NAO對1至8號分子軌道上的貢獻。因此輸入
-1 //進入片段設定界面
a 2,3 //將2、3原子上的NAO都加入片段
q //保存片段
1 //計算片段對一批軌道的貢獻
1-8 //要考察的軌道序號范圍
結果如下所示
Orb.# Core Valence Rydberg Total
1 0.000% 0.119% 0.002% 0.121%
2 0.000% 18.556% 0.054% 18.611%
3 0.000% 26.557% 0.018% 26.576%
4 0.000% 7.076% 0.012% 7.089%
5 0.000% 0.000% 0.000% 0.000%
6 0.000% 35.482% 46.832% 82.314%
7 0.000% 27.543% 61.561% 89.104%
8 0.000% 32.433% 29.553% 61.986%
由于氫沒有內核軌道,所以Core項全部為零。Valence、Rydberg項分別是指這個片段內所有Val型、Ryd型NAO的貢獻。Total就是Core、Valence、Rydberg項的總和。前五個分子軌道是占據軌道,與之前討論的一致,Ryd的成份非常低,這時可以說這兩個氫對MO3的貢獻是26.557%。而對于后面的虛軌道,則Ryd成份很大,不能忽略。由于Ryd型軌道物理意義不清,與原子軌道沒有直接對應關系,難以討論,所以此時就不建議再用NAO方法來分析虛軌道成份了,而更適合Hirshfeld方法分析原子的貢獻。
7 總結
本文將軌道成份的計算原理、要注意的問題和具體操作過程都做了討論。計算軌道成份不能稀里糊涂,一定要思路清楚。很多人急于算出軌道成份,在網上隨便找一個來路不明的號稱能計算軌道成份的軟件就用于自己的計算,或者網上道聽途說一個計算方法,比如什么“軌道成份就是系數的平方”,就按照這個方法算,這都是很危險的。
對于計算基函數的貢獻,最建議用SCPA方法;對于計算原子軌道的貢獻,最建議用NAO方法;對于計算原子的貢獻,Hirshfeld、Hirshfeld-I、Becke方法以及將NAO方法得到的原子軌道的貢獻進行加和的方法都推薦使用。注意NAO方法不適合分析虛軌道。Multiwfn也支持AIM方法計算原子的貢獻,但由于此方法耗時高,所以一般不用。
加和SCPA方法得到的基函數的貢獻來得到原子軌道的貢獻,以及進一步加和得到原子軌道的貢獻,雖然原理上不如上面建議的方法,但是只要沒有彌散函數,且基組不平衡性較小時,且分析的不是較高能態虛軌道時,結果一般沒什么問題。在使用SCPA方法分析原子軌道成份時最好(但不是必須)使用球諧型高斯函數以避免旋轉依賴性,也避免笛卡爾型高斯函數引入的額外的基函數影響分析結果。
最后再次提醒一下,討論某個基函數、某個原子軌道貢獻時必須指明分子朝向,否則沒有意義。不指定朝向時只能討論殼層的貢獻。