氫分子間色散相互作用能的自然軌道泛函

盛曉偉

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

1 引 言

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

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

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

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

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

圖1 兩相互平行的氫分子結構圖Fig. 1 The geometry of two parallel hydrogen molecules

該體系的基態行列式為ag，b3u兩分子軌道雙占據，即：

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

(1)

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

(2)

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

(3)

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

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

(4)

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

(5)

ΨE(a)所包含的6個行列式與Ψ0均只有兩個軌道不同，利用Slater-Condon規則可計算得到Ψ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)類激發組態所貢獻的電子間相互作用能，其具體形式如下：

〈1ag1b3u|nb3umag〉]

(12)

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

〈1ag1b3u|nb2umb1g〉]

(13)

〈1ag1b3u|nb1umb2g〉]

(14)

〈1ag1b3u|naumb3g〉]

(15)

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

(16)

3 計算結果與討論

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

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

(17)

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

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

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

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

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

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

4 結 論

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

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

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

(19)

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