密度矩陣泛函理論及其計算方法

IRIMIA Marinela,王 堅

(1.湖州師范學院 國際學院,浙江 湖州 313000; 2.湖州師范學院 理學院,浙江 湖州 313000)

0 引 言

客觀世界是由原子和分子組成的,它們遵從量子力學的規律.因此,原則上只要求解量子力學方程,就可以知道材料的各種物理化學性能.與實驗比較,求解量子力學方程的計算成本輕,因此其已成為材料設計的重要手段.

一般的原子、分子和固體系統,其物理性質可通過求解下面的量子力學方程,即薛定諤方程來求解：

HΨ=EΨ.

采用原子單位,哈密頓算符H由以下幾部分組成：

其中,第一項是電子的動能算符,第二項是原子核與電子間的庫侖吸引勢能,第三項是電子間的庫侖排斥勢能,第四項是原子核間的庫侖排斥勢能.假定波恩-奧本海默近似成立,即在考慮電子運動時,原子核可以看作是靜止的,其原因是核比電子的質量大得多.對只有一個電子的氫原子,這一方程可以嚴格求解,且解的結果已被光譜實驗完全證實,由此奠定了量子力學的基礎.

對由兩個電子或更多電子組成的原子、分子和固體系統,這個方程無法嚴格求解,只能通過近似方法求解.因為電子是費米子,波函數對于任意兩個電子交換要具有反對稱性.為滿足這一性質,波函數Ψ可采用行列式的形式,因為行列式的兩行或兩列交換后,行列式值就變符號,正好符合反對稱性要求.最簡單的波函數形式是用一個行列式來表達的,這種方法也稱為Hartree-Fock方法.為提高精確度,可采用多個行列式的線性組合來近似波函數,這種方法稱為組態相互作用(簡稱CI).一個行列式對應一個組態,因此這是一種多組態的波函數方法.這在理論上是理想的,但隨著電子數目的增加,行列式數目成指數式增加,因此這種方法很難推廣到多電子系統.

20世紀50年代,人們認識到計算系統的能量不需要波函數的全部信息,只需要密度矩陣.例如,一階密度矩陣的定義是[1]：

它是除了第一個電子的坐標外,對其余電子的坐標進行積分后得到的結果.每個電子有3個空間坐標和1個自旋坐標,積分包含空間坐標的積分和自旋坐標的求和.根據定義,它是厄米矩陣,因此可以對角化,本征值為實數,

其中,本征函數χi稱為“自然軌道”;本征值λi稱為自然軌道χi的占有數,0≤λi≤1.

因為電子是全同粒子,只要知道一個電子的動能,總動能就等于它的N倍.利用一階密度矩陣,電子的動能可以簡化為：

同樣,原子核與電子的庫侖吸引勢能也可簡化為單個電子的積分：

這里的γ(1,1)是γ(1′,1)的對角元,稱為電子密度.

我們還可以定義二階密度矩陣(或電子對密度)：

它是波函數和它的復共軛中除了第一和第二個電子的坐標外,對其余電子的坐標進行積分后剩下的結果.其中,N(N-1)是電子對的總數目.電子間的庫侖勢能是兩體算符,利用電子的全同性,只要考慮一對電子的庫侖勢能,然后乘上電子對總數即可,因此可以用Γ(1,2)來表示,即：

綜上所述,與電子相關的能量可以用密度矩陣表示為：

因此,只要知道相當于密度矩陣γ(1′,1)和Γ(1,2)的信息,而不需要波函數的全部信息,就可以獲得電子系統的能量.

1964年,Kohn和Hohenberg提出了密度泛函理論(簡稱DFT)[2].他們用反證法證明電子系統基態的能量是電子密度的泛函.然而這一定理只是個存在性證明,這個密度泛函的具體表達形式仍然未知.按照上面的分析,得出表1所示的結果.只有原子核與電子的庫侖吸引勢能可以用密度表達,電子能量的其余兩項都不能用密度嚴格表達,所以只能尋找近似模型.Kohn和Sham最早通過均勻電子氣體模型給出了一個稱為“局域密度近似(LDA)”的泛函模型[3].固體中的電子運動近似均勻電子氣體,因此這個模型在固體物理領域最先得到應用.但原子或分子中的電子密度一般是不均勻的,因此20世紀90年代,Perdew等對密度不均勻體系提出了“廣義梯度近似”模型[4].DFT才開始在原子分子物理和化學領域得到應用.

表1 密度泛函方法與密度矩陣泛函方法的比較

1975年,Gilbert提出了密度矩陣泛函理論[5],這種理論用一階密度矩陣作為泛函變量.與DFT方法比較,這一方法的優點是電子動能,以及原子核與電子的庫侖勢能都可以用一階密度矩陣嚴格表達,只有電子間的庫侖勢能需要尋找泛函模型.密度矩陣泛函方法與密度泛函方法的比較見表1.從泛函模擬的角度看,這種理論方法比DFT更有希望.因此不少DFT學者,如M Levy、E K U Gross、E J Baerends等,都成為密度矩陣泛函理論的倡導者.

1 密度矩陣泛函

用CI波函數展開Γ(1,2),電子間的庫侖勢能可表示為：

其中,〈ij|kl〉為雙電子積分.在波函數方法中,Γij,kl與軌道指標ij和kl所在的行列式的波函數展開系數有關.在密度矩陣泛函理論中,自變量是一階密度矩陣,包括自然軌道χi及其占有數λi,因此原則上Γij,kl是自然軌道和占有數的函數.Γij,kl有4個下標,每個下標有M個軌道選項,因此Γij,kl共有M4項.如果軌道數M=10,那么就有1萬個關于Γij,kl的泛函,這么多泛函逐個落實是不現實的.為使計算可行,泛函設計時常作近似處理,一般只保留指標為Γij,ij和Γij,ji的項,而忽略其他指標項,文獻[6]中把這種泛函模型稱為“JK”型.原本Γij,kl有M4項,實際模擬的只有M2項.

同時,由于雙電子積分〈ij|kl〉已包含軌道,所以在設計Γij,kl的泛函模型時往往忽略軌道變量,而是假定Γij,kl只與占有數λi有關,這種簡化的泛函模型稱為“代數型”[7].

對實際的分子,我們用波函數推導的電子對密度對文獻中發表的密度矩陣泛函進行逐項比較.結果發現,“JK”型和“代數型”泛函都不能重復從波函數得到Γij,kl,它們不具有一一對應的關系.這說明這些泛函模型都存在唯一性、對稱性和普適性問題[7-8].

為了克服這些問題,經過長期的摸索,我們對以下形式的Γ(1,2)分解專門進行研究：

Γ(1,2)=γ(1,1)γ(2,2)-γ(1,2)γ(2,1)+λ(1,2),

其中,第一項是獨立粒子引起的電子對密度,第二項是交換效應.這兩項都可用一階密度矩陣表達.后一項λ(1,2)稱為“cumulant”,其對應的電子能量為：

著名DFT學者M Levy把Ecum稱為密度矩陣泛函理論的關聯能,它是能量泛函中唯一需要用泛函模擬的部分.因為關聯能使能量降低,所以Ecum總是取負數.對Ecum進行模擬可避免模擬Γij,kl時符號的不確定性.這個符號的不確定性問題,文獻中稱為“phase dillema”[9].通過積分,Ecum隱含了所有可能指標ij,kl引起的貢獻,避免了“JK”型和“代數型”泛函的近似性,也避免了模擬Γij,kl時有M4項之多的麻煩.

信息論之父Shannon曾提出,“信息熵”函數S=-∑ipilogpi可用來描述統計分布pi.我們將占有數分布λi看作一個統計分布,用S=-∑iλilogλi作為與這個分布相聯系的信息熵,發現這個熵函數與關聯能Ecum有較好的線性關系[10]：

Ecum=-κS-b,

其中,κ和b是與系統有關的常數,可以通過結合能、鍵長等理論或實驗數據來擬合.按照這一表達式,當熵取最大值時,關聯能達到最小值.關聯能越低,意味著系統越穩定,這正是我們期望的,也是符合物理要求的.考慮到電子是費米子,應符合狄拉克-費米統計,信息熵的嚴格表達形式應為：

S=-∑i[λilogλi+(1-λi)log(1-λi)].

2 自洽場方程

Gilbert在提出密度矩陣泛函理論時,發現當占有數λi為分數時,軌道能級都是簡并的[5].這一結果與實驗或DFT計算的結果不符,他在文章中表示意外.因為能級簡并,無法獲得DFT那樣的Kohn-Sham方程,計算時只能依靠非線性優化方法,計算效率很低.這一問題成為密度矩陣泛函理論的一個歷史性問題.

用信息熵表達關聯能后,系統的電子能量可表示為：

E=∑iλihii+Y-κS-b,

其中,與單電子有關的能量,hii=〈χi|h|χi〉,包括動能以及原子核對電子的庫侖吸引勢能,單電子算符h的定義為：

第二項Y包含雙電子能量中可以用一階密度矩陣表達的部分：

其中,

稱為雙電子積分.

考慮到自然軌道的正交歸一性,以及占有數的求和性質∑iλi=N,我們用拉格朗日函數方法構造一個無條件優化的變分函數：

其中,μ和Λij是拉格朗日因子.

其中,算子Jj和Kj與Hartree-Fock方法中定義的庫侖算子和交換算子完全一樣,

為簡化書寫,定義算子：

這個算子與Hartree-Fock方法中的Fock算子比較,只多了系數λj.

(λi-λj)〈j|f|i〉=0.

若使這個方程成立,則f|i>應具有本征方程的形式,即：

fχi(1)=iχi(1).

〈i|f|i〉+κ[logλi-log(1-λi)]=μ.

利用〈i|f|i〉=i,可以解得：

編寫量子力學計算軟件是一項長期積累的工作.因為本文的本征值方程同Hartree-Fock方程非常接近,所以我們可以把Hartree-Fock程序作為出發點來開展研究.Hartree-Fock程序已研發了90多年,文獻資料相當豐富,程序代碼也很多,如原子的高斯基函數庫、雙電子積分、矩陣對角化等方面的內容.利用這些文獻資料,可以大大加快新方法計算程序的編寫.

Hartree-Fock方法是一個非常有效、穩定的計算方法.由于本征值方程的相似性,本方法也具有相當的穩定性、可靠性和可行性.在自洽場收斂方面,本文的方法比Hartree-Fock方法或DFT的Kohn-Sham方程更有優勢,其原因是軌道占有數可連續變化,而在Hartree-Fock方法或Kohn-Sham方法中,占有數是跳躍性的,從1變成0或從0變成1,這容易使自洽場在收斂過程中,在近乎簡并的軌道上發生跳躍,從而出現振蕩現象而難于收斂.

3 計算舉例：氫分子的分解

氫分子的分解問題是密度泛函理論的一個弱點,楊偉濤等曾專門在Science雜志上討論這個問題[12].在DFT計算中,密度用軌道來展開,軌道上填充電子的方式按能級由低到高的次序來填充,有電子占據的軌道,占有數為1,否則為0.這一填充規則與Hartree-Fock方法一致,相當于用單個行列式表達波函數.在波函數計算領域,人們早就知道,單個行列式不能描述具有靜態關聯的系統,比如遠離平衡態的氫分子.

圖1為利用不同方法計算的氫分子能量隨原子間距變化的結果.當氫分子分解成兩個氫原子后,能量等于兩個氫原子的能量之和,每個氫原子的能量是-0.5個原子單位,所以總能量應趨向于-1.由圖1可知,Hartree-Fock(HF)方法、DFT的LDA泛函方法都沒有收斂到正確的結果(-1),而本文的方法和CI方法都給出了正確的結果.因為在密度矩陣泛函理論中,軌道占有數λi可以為分數,這與CI波函數方法一樣,所以它們的結果一致.這一結果說明,占有數可以為分數,密度矩陣泛函理論是一個比DFT方法更靈活、更完美、更精確的理論方法.對于四個氫原子組成的正方形結構和50個氫原子組成的一維分子鏈,我們也獲得了類似的結果[13].

圖1 氫分子分解過程的能量曲線

4 結 語

電子結構計算的核心問題是電子關聯.在波函數計算方法中,電子關聯是通過增加行列式數目來計算的,這種方法收斂很慢.而在密度矩陣泛函理論方法中,電子關聯是通過軌道占有數來引入的.波函數方法目前只能處理10個左右電子的關聯.在沒有發現本征值方程前,密度矩陣泛函理論的計算采用非線性優化方法,計算量遠比DFT的Kohn-Sham方法復雜.本文的本征值方程方法,使密度矩陣泛函理論的計算工作量相當于單個行列式的Hartree-Fock方法;在同樣的計算機條件下,關聯的電子數目可成千上萬倍地增加.這為復雜分子的研究創造了條件,是一個歷史性突破.