盛曉偉

(安徽師范大學(xué) 物理與電子信息學(xué)院， 蕪湖 241000)

1 引 言

色散相互作用源于電子間的非定域關(guān)聯(lián)效應(yīng)，它是原子分子間多種相互作用中的一種. 由于其在氣液相變、液固相變以及冷原子碰撞等重要物理過程中起關(guān)鍵性的作用[1-3]，從而引起了人們的廣泛關(guān)注. 目前對色散相互作用能的理論計算主要有以下三種方法：半經(jīng)驗的勢模型法、電子相關(guān)的從頭計算法以及密度泛函法. 在半經(jīng)驗勢模型方法中，體系的色散相互作用能一般用一個包含若干參數(shù)的解析式進(jìn)行表示. 然而勢模型中參數(shù)值的確定，依賴于實(shí)驗或者精確的從頭計算[4,5]，因此勢模型方法在實(shí)際應(yīng)用中受到一定的限制；電子相關(guān)的從頭計算方法可以很精確地計算出色散相互作用能，并且不依賴于任何參數(shù)[6,7]. 然而這類方法的計算量隨體系的增大而快速增長，很難推廣到大分子體系[8]；密度泛函理論在處理色散相互作用問題上，一直沒有得到很好的解決. 基于局域密度和廣義梯度密度近似的密度泛函，無法描述具有非定域性質(zhì)的色散相互作用[9]. 近些年有很多設(shè)法計算色散相互作用能的半經(jīng)驗密度泛函相繼被報道[10-14]，然而隨著密度泛函對色散相互作用能計算精度的增加，泛函的形式越來越復(fù)雜，計算量也越來越大[15]. Baerends 等[16]在文章中指出：“截止目前，無法構(gòu)造出一個能夠精確處理色散相互作用能的簡單密度泛函”. 可見，發(fā)展一種能夠?qū)υ臃肿娱g色散相互作用實(shí)現(xiàn)高效計算的理論方法仍然是個有待解決的基礎(chǔ)科學(xué)問題.

密度矩陣泛函以非定域量一階密度矩陣為基本變量(一階密度矩陣可以利用其本征函數(shù)自然軌道和本征值自然軌道占居數(shù)展開，因此密度矩陣泛函也稱為自然軌道和自然軌道占居數(shù)泛函)，在計算量上與密度泛函接近，例如：Buijse-Baerends(BB)泛函的計算量按基組的4次方增長，這和一般密度泛函的計算量相當(dāng)[17]. 密度矩陣泛函和密度泛函的最大區(qū)別在于：體系的動能可以用一階密度矩陣進(jìn)行精確表示，無需構(gòu)造Kohn-Sham無相互作用體系，這使其在處理電子相關(guān)問題上較密度泛函更具有優(yōu)勢. 如：BB泛函成功地描述了一些小分子體系中的電荷轉(zhuǎn)移以及化學(xué)鍵的離解過程[18]；Sharma的冪函數(shù)能夠正確計算出半導(dǎo)體以及過渡金屬氧化物絕緣體的帯隙[19]；最近該泛函又被成功的運(yùn)用到同性電子氣和一維以及二維異性Hubbard模型中[20-22].

密度矩陣與色散相互作用都具有非定域的性質(zhì). 因此，利用密度矩陣泛函對原子間色散相互作用進(jìn)行計算引起了廣泛的關(guān)注. Cioslowski和Pernal[23]基于對兩個弱相互作用體系總能量漸進(jìn)形式的研究以及自旋軌道在兩個體系中具有完全定域性質(zhì)的假設(shè)，得到了描述色散相互作用能的密度矩陣泛函所必須滿足的條件，但其并沒有給出該密度矩陣泛函的具體形式. Gritsenko和Baerends[24]通過分析自然軌道和自然軌道占據(jù)數(shù)表象下雙電子體系的三重態(tài)波函數(shù)，得到了三態(tài)氫分子在范德瓦爾斯勢阱位置附近的勢能曲線. 研究結(jié)果表明：該體系中的色散相互作用能在自然軌道表象下收斂速度非?？? 10個自然軌道組合成的組態(tài)波函數(shù)已經(jīng)能夠計算出體系 90% 以上的色散相互作用能，然而在Kohn-Sham和Hartree-Fock軌道基組下，得到同樣比例的色散相互作用能卻都需要60個以上軌道. 由此可見：自然軌道相較于Kohn-Sham軌道和Hartree-Fock軌道更適合描述色散相互作用，密度矩陣泛函較密度泛函在處理色散相互作用上更具有優(yōu)勢. 2007年，Piris等[25]提出了一個完全由交換和庫侖積分構(gòu)成的密度矩陣泛函(PNOF2). 該泛函能夠給出基態(tài) He2的勢阱位置Re，但只能得到50%的勢阱深度De. 2013年，本課題組分析了基態(tài) H2和 He2的對關(guān)聯(lián)函數(shù)，基于色散相互作用的非定域特性，首次辨認(rèn)出對關(guān)聯(lián)函數(shù)中描述色散相互作用的項，在此基礎(chǔ)上找到了基態(tài)H2和He2的色散相互作用能自然軌道泛函[26]. 本文將在此基礎(chǔ)上進(jìn)一步分析兩相互平行的氫分子間色散相互作用能的自然軌道泛函形式.

2 兩相互平行的氫分子間色散相互作用能

兩相互平行的氫分子體系屬于 D2h群， 該群有8種不可約表示:ag，b3u，b1g，b2u，b2g，b1u，b3g和au. 圖1展示了該體系的幾何結(jié)構(gòu).

圖1 兩相互平行的氫分子結(jié)構(gòu)圖Fig. 1 The geometry of two parallel hydrogen molecules

該體系的基態(tài)行列式為ag，b3u兩分子軌道雙占據(jù)，即：

Ψ0=|lagαlagβlb3uαlb3uβ|

(1)

前期研究工作指出， 體系的色散相互作用能可通過計算一些特定雙重激發(fā)組態(tài)所貢獻(xiàn)的電子間相互作用能而得到. 對于基態(tài) He2分子，我們發(fā)現(xiàn)E類雙重激發(fā)組態(tài)為該體系色散相互作用能的主要貢獻(xiàn)項[16]. 本文的研究體系與He2分子具有類似的雙重激發(fā)組態(tài). 因此，E類雙重激發(fā)組態(tài)在兩相互平行的氫分子體系中必然也是其色散相互作用能的主要貢獻(xiàn)項. 該體系中E類雙重激發(fā)組態(tài)為如下四種組合形式：

(2)

電子間相互作用能計算公式如下：

(3)

Γ(1,2)為二階密度矩陣，其定義如下：

Ψ*(1′,2′,3…N)d3…dN

(4)

其中Ψ(1,2,3…N)為該體系的組態(tài)波函數(shù)，為了得到E類激發(fā)組態(tài)所貢獻(xiàn)的電子間相互作用能，我們考慮ΨE(a,b,c,d)激發(fā)組態(tài)和Ψ0的乘積對二階密度矩陣的貢獻(xiàn). 現(xiàn)以ΨE(a)為例來作此分析，

(5)

ΨE(a)所包含的6個行列式與Ψ0均只有兩個軌道不同，利用Slater-Condon規(guī)則可計算得到ΨE(a)與Ψ0乘積所得到如下的六個二階密度矩陣：

[mag(1)nb3u(2)lag(1′)lb3u(2′)-

mag(1)nb3u(2)lb3u(1′)lag(2′)-

nb3u(1)mag(2)lag(1′)lb3u(2′)-

nb3u(1)mag(2)lb3u(1′)lag(2′)]ββββ

(6)

[mag(1)nb3u(2)lag(1′)lb3u(2′)-

mag(1)nb3u(2)lb3u(1′)lag(2′)-

nb3u(1)mag(2)lag(1′)lb3u(2′)+

nb3u(1)mag(2)lb3u(1′)lag(2′)]αααα

(7)

[mag(1)αnb3u(2)βlag(1′)βlb3u(2′)α-

mag(1)αnb3u(2)βlb3u(1′)αlag(2′)β-

nb3u(1)βmag(2)αlag(1′)βlb3u(2′)α+

nb3u(1)βmag(2)αlb3u(1′)αlag(2′)β]

(8)

[mag(1)βnb3u(2)αlag(1′)βlb3u(2′)α-

mag(1)βnb3u(2)αlb3u(1′)αlag(2′)β-

nb3u(1)αmag(2)βlag(1′)βlb3u(2′)α+

nb3u(1)αmag(2)βlb3u(1′)αlag(2′)β]

(9)

[mag(1)αnb3u(2)βlag(1′)αlb3u(2′)β-

mag(1)αnb3u(2)βlb3u(1′)βlag(2′)α-

nb3u(1)βmag(2)αlag(1′)αlb3u(2′)β+

nb3u(1)βmag(2)αlb3u(1′)βlag(2′)α]

(10)

[mag(1)βnb3u(2)αlag(1′)αlb3u(2′)β-

mag(1)βnb3u(2)αlb3u(1′)βlag(2′)α-

nb3u(1)αmag(2)βlag(1′)αlb3u(2′)β+

nb3u(1)αmag(2)βlb3u(1′)βlag(2′)α]

(11)

將式(6-11)代入到式(3)，得E(a)類激發(fā)組態(tài)所貢獻(xiàn)的電子間相互作用能，其具體形式如下：

〈1ag1b3u|nb3umag〉]

(12)

同理可得E(b)，E(c) 和E(d)類激發(fā)組態(tài)所貢獻(xiàn)的電子間色散相互作用能.

〈1ag1b3u|nb2umb1g〉]

(13)

〈1ag1b3u|nb1umb2g〉]

(14)

〈1ag1b3u|naumb3g〉]

(15)

通過分析式(12-15)可知：色散相互作用能來源于自旋相同部分正好是自旋相反的兩倍. 然而，色散相互作用為電子間的電磁相互作用與電子的自旋并無直接關(guān)系，自旋相同與自旋不同電子間色散相互作用能應(yīng)該相等. 對于He2體系，我們前期的工作指出E類激發(fā)組態(tài)波函數(shù)僅能得到自旋不同電子間色散相互作用能的一半. 估算體系的色散相互作用能應(yīng)在式(12-15)求和的基礎(chǔ)上再乘以系數(shù)4/3. 因此，該體系的色散相互作用能可近似為下式：

(16)

3 計算結(jié)果與討論

本文我們討論兩相互平行的氫分子在以下兩種情況下的色散相互作用能.

第一種情況：兩個氫分子都處于平衡位置附近(RH-H=1.4 a.u.)，研究色散相互作用能隨兩氫分子間距的變化(6 a.u.

(17)

由于SAPT理論對處于非平衡體系將失效. 因此，式(16)計算得到的第二種情況下色散相互作用能將和遠(yuǎn)程的色散項C6/R6直接進(jìn)行對比. 該色散項能夠非常精確的估算體系遠(yuǎn)程的色散相互作用能.C6通過直接積分 Casimir-Polder 方程得到[27].

α(iω)為偶級極化率，由CCSD理論水平下的響應(yīng)理論計算得[28]. 表I和表II分表羅列出以上兩種情況下的計算結(jié)果. 以上所有計算均在 aug-cc-pvTZ 基組下進(jìn)行.

對于第一種情況， 式(16)的計算結(jié)果與 SAPT 非常接近. 最大的誤差百分比僅為8.33 %. 當(dāng)R< 8.0 a.u. 時, SAPT 理論計算得到的色散相互作用能比式(16)要大. 但在R> 8.0 a.u. 區(qū)域， 式(16)的計算結(jié)果要比 SAPT大.

表I 兩相互平行的氫分子間色散相互作用能(RH-H=1.4 a.u.). 所有的計算在aug-cc-pvTZ基組下進(jìn)行，數(shù)值單位為原子單位.

在第二種情況下，式(16)同樣和遠(yuǎn)程的色散項C6/R6符合的較好. 當(dāng)RH-H=1.4 a.u. 時有最大誤差百分比 17.94 %. 該誤差百分比將隨RH-H的增大而逐漸減小直到RH-H> 3.5 a.u.. 式(16)的計算結(jié)果在這些位置都比C6/R6大. 這是由于遠(yuǎn)程的色散項C6/R6并不包含高階色散項. 因此，C6/R6一般比精確值要小. 例如：表I顯示在RH-H=1.4 a.u.處，式(16)與 SAPT 的誤差百分比僅有-2.57 %. 當(dāng)RH-H> 3.5 a.u.時，式(16)的結(jié)果比色散項C6/R6的值還要小. 這是因為在這些位置將會有更多的組態(tài)波函數(shù)具有和基態(tài)行列式Ψ0可比擬的展開系數(shù)(氫分子處于離解區(qū)域時非動態(tài)關(guān)聯(lián)效應(yīng)明顯). 這時僅考慮ΨE(a)和Ψ0的乘積對二階密度矩陣的貢獻(xiàn)不再是個很好的近似,多參考組態(tài)效應(yīng)需要考慮，本文對此不作深入討論. 由此可見，式(16)對兩氫分子間色散相互作用能的確能夠給出較準(zhǔn)確的估算.

表II 兩相互平行的氫分子間色散相互作用能(兩氫分子間距R= 7.0 a.u.). 所有的計算在aug-cc-pvTZ基組下進(jìn)行，數(shù)值單位為原子單位.

4 結(jié) 論

本文通過計算來源于E類雙重激發(fā)組態(tài)波函數(shù)的電子間相互作用能來估算兩相互平行的氫分子間色散相互作用能. 通過和高精度的從頭計算(SAPT和CCSD)進(jìn)行比較，表明該方法合理有效.

研究表明,精確描述該體系中的色散相互作用能的自然軌道泛函形式為：

[〈kgku|numg〉-〈kgku|mgnu〉]

(19)

其中，k取遍所有軌道占據(jù)數(shù)nkg/u>1的強(qiáng)占據(jù)軌道，m,n表示所有占據(jù)數(shù)nmg,nnu< 1的弱占據(jù)軌道.fxcdisp為自然軌道占據(jù)數(shù)泛函. 該泛函中的每個積分含有4個軌道，并非是目前密度矩陣泛函廣泛使用的交換和庫倫積分泛函. 因此，現(xiàn)有的密度矩陣泛函并不包含有體系的色散相互作用能. 要得到精確計算色散相互作用能的密度矩陣泛函還需要知道fxcdisp的具體形式，這將是我們今后的工作.