基于無軌道密度泛函理論的金屬鈉板系統能量泛函*

劉潤鳳,胡奇宏,周 芳,程曉迪,王小云,黃勇剛

(吉首大學物理與機電工程學院,湖南 吉首 416000)

密度泛函理論(Density Functional Theory,DFT)能夠準確描述各種材料的基態性質,是計算凝聚態物理電子結構的重要方法[1-2].然而,DFT的計算成本與體系中電子數量的立方成正比,因此往往僅能用于處理具有周期結構或者超小規模的體系,而對于大規模多電子體系,DFT的應用受到極大限制.相比之下,無軌道密度泛函理論(Orbital-Free Density Functional Theory,OF-DFT)的計算成本較低,其實質可以看作是只考慮最小特征值的自恰Kohn-Sham方程[3-6],能直接通過電子密度計算總能量,因此可以突破DFT在模擬尺度上的極限.OF-DFT的主要優點是線性縮放計算成本,可應用于非常大的系統.此外,在有限元方法中實施OF-DFT,進一步擴展了OF-DFT的應用廣泛性[7].將OF-DFT與有限元方法結合,可以實現對各種具有復雜形貌特征的金屬體系的基態性質的精確計算.基于此,筆者擬采用OF-DFT對金屬鈉板基態性質進行深入探討.

使用OF-DFT計算體系的基態密度和功函數時,能量泛函的選擇至關重要,它直接影響方法的準確度.能量泛函由動能泛函和交換關聯泛函組成,目前相關研究工作中,許多近似的動能泛函被用于描述體系的相互作用能[4-5,8-10].最常采用的動能泛函是由Thomas-Fermi (TF) 和von Weizsacker (vW)組合而成的TFvW動能泛函,交換關聯泛函通常采用基于局域密度近似(Local Density Approximation,LDA)的交換關聯勢[7,11-12].雖然TFvW動能泛函被廣泛采用,但vW動能泛函缺乏微觀對應,其比重尚未確定.研究發現,vW動能泛函的比重越大,電子溢出效應越明顯,溢出功也越大[13],因此確定vW動能泛函的比重成為一個亟待解決的問題.另外,雖然不同的LDA交換關聯勢被廣泛采用[7],但其對計算結果的影響尚不明確.

為了確定金屬鈉板獲得準確基態性質所需的vW動能泛函,筆者首先研究了OF-DFT及DFT的相關理論和模型,并介紹了計算基態密度和功函數的方法,進一步研究了vW動能泛函的比重在不同LDA交換關聯泛函中對基態密度誤差與功函數誤差的影響.

1 理論與模型

1.1 無軌道密度泛函理論

在OF-DFT中,通常通過求解歐拉方程來獲得體系的基態性質[4-5,14-15],即

(1)

其中:n表示基態密度;φ0表示靜電勢;μ表示化學勢.能量方程G[n]包含電子系統相互作用的動能(Ts)和交換相關能(EXC),即

G[n]=Ts[n]+EXC[n].

(2)

τ(n,w)=τvW(n,w)+τTF(n),

其中

τvW(n,w)=An-1w,

τTF[n]=Bn5/3.

為了求解(1)式,對方程式兩側做梯度運算,即

(3)

其中E0表示靜電場.(3)式稱為量子流體動力學靜態方程.將(3)式與Maxwell方程耦合,可得

(4)

其中n+表示正電荷密度.在凝膠模型中,正電荷密度在球體內是均勻恒定的,而在球體邊緣外部突然下降至0.將(3)式與(4)式結合,可得

(5)

通過(5)式可以計算體系的基態密度.本研究中,所有的OF-DFT計算都使用有限元(FEM)軟件COMSOL的弱形式模塊來實現.

1.2 OF-DFT中能量泛函對基態性質的影響

圖1 金屬鈉無限大平板示意Fig. 1 Schematic Diagram of an Infinite Plate of Sodium Metal

使用OF-DFT計算金屬鈉無限大平板的基態密度和功函數,并將結果與DFT計算結果進行對比.假定無限大平面在y軸和z軸方向無限大,在x軸方向的厚度為L,那么本研究只需關注無限大平面在x軸方向上的基態性質即可.金屬鈉無限大平板模型如圖1所示.

首先研究OF-DFT中不同能量泛函對基態密度和功函數的影響.圖2和圖3的縱軸使用對數坐標,這可以明顯比較基態密度在金屬邊界處的衰減速度.不同vW動能泛函比重對基態密度的影響如圖2所示.圖2中,交換關聯勢均使用了LDA-WG.從圖2可知,隨著λ的逐漸增大,基態密度的衰減速度逐漸降低.在對數坐標下,基態密度在金屬邊界處的衰減近似為穩定的直線,這說明基態密度在金屬表面的衰減可以看作近似e指數衰減[11].vW動能泛函的比重可以調控這一衰減因子,vW動能泛函比重越大,衰減因子越小,則基態密度在金屬邊界處的衰減速度越慢,體系的溢出效應越明顯.不同LDA交換關聯勢對基態密度的影響如圖3所示.圖3中,參數vW動能泛函比重取值為1.從圖3可知,由LDA-WG,LDA-PZ和LDA-GL計算得到的基態密度幾乎沒有明顯區別,即使使用了對數坐標,差距仍然較小.

以上結果說明,vW動能泛函的比重會影響基態密度的計算結果.vW動能泛函比重的增加,實質上可以理解為增加了體系的電子溢出效應,因此基態密度的衰減速度也會受到影響.相比之下,選擇不同的LDA交換關聯勢對基態密度的影響較小.

圖2 不同vW動能泛函比重對基態密度的影響Fig. 2 Ground State Densities Obtained with Different vW Kinetic Energy Functional Proportions

圖3 不同LDA交換關聯勢對基態密度的影響 Fig. 3 Ground State Densities Obtained with Different LDA Exchange-Correlation Potentials

另一方面,功函數同樣是重要的基態性質.接下來筆者介紹使用OF-DFT求解功函數的方法.OF-DFT可以看作是只考慮最小特征值的自恰Kohn-Sham方程,由于OF-DFT中只考慮最小的特征值,體系中所有電子能量相同,因此(1)式中的化學勢(μ)即為體系功函數,即

W=Ts[n]+EXC[n].

(6)

根據功函數的定義,(6)式需取無窮遠處的值.利用(6)式計算得到的化學式,其在空間各處皆為常數.不同vW動能泛函比重對功函數的影響如圖4所示.圖4中交換關聯勢均使用了LDA-WG.從圖4可知,功函數隨著λ的增大而增大,λ由0.6增大至1.0時,功函數由3.218 eV增大至3.538 eV.不同交換關聯勢對功函數的影響如圖5所示.圖5中,參數vW動能泛函比重值為1.由圖5可知,LDA-WG,LDA-PZ和LDA-GL計算得到的功函數分別為3.538,3.478,3.368 eV.

圖4 不同vW動能泛函比重對功函數的影響Fig. 4 Work Functions Obtained by Different vW Kinetic Energy Functional Proportions

圖5 不同LDA交換關聯勢對功函數的影響Fig. 5 Work Functions Obtained by Different LDA Exchange-Correlation Potentials

計算結果說明,vW動能泛函比重和交換關聯勢的改變都會影響功函數的計算結果.與上述基態密度的變化類似,vW動能泛函比重的增加,實質上可以理解為增加體系的電子溢出效應,從而引起功函數增加.此外,筆者發現,與基態密度不同,交換關聯勢的改變會影響功函數的計算結果,原因在于交換關聯勢在遠場起主要作用.

1.3 密度泛函理論

為了驗證OF-DFT計算結果的準確性,使用開源軟件OCTOPUS計算DFT模塊[16-17],并將DFT計算結果與OF-DFT計算結果進行對比.在堿金屬中,價電子的費米波長遠大于金屬晶格常數,離子的贗勢對電子結構沒有顯著影響.因此,可以用均質的正電荷背景代替離散的離子結構,即使用凝膠模型[18-20].此時,金屬板鈉體系可以看作是Wigner-Seitz半徑(rs)為4的凝膠模型,進一步求解Kohn-Sham方程,這一處理能極大節省計算資源.在OCTOPUS中,可以使用不同的交換關聯勢,為了能與動能泛函在OF-DFT中的作用進行比較,本研究使用與OF-DFT一致的交換關聯勢.

DFT中,不同LDA交換關聯勢下基態密度計算結果如圖6所示.從圖6可知,與OF-DFT類似,DFT中使用不同交換關聯勢計算得到的基態密度并無較大差別.不同LDA交換關聯勢下功函數計算結果如圖7所示.這類計算在文獻[21-23]中已出現過.從圖7可知,功函數呈現振蕩行為,其原因在于隨著平板厚度的逐漸增加,閑帶的能量逐漸變低,在一定的平板寬度上,各自的最低閑帶接觸并最終穿過費米曲面.這些接觸點以λF/2的周期出現(λF≈3.27rs),筆者觀察到,功函數在λF=8時基本達到穩定振蕩.若進一步增加平板厚度,則需要更多的計算資源,同時會導致收斂困難,因此在本研究中,選擇λF=9.5和λF=10時的平均值作為DFT計算得到的功函數,這一方案獲得的研究結果與文獻[14]完全一致.使用LDA-WG,LDA-PZ和LDA-GL交換關聯勢計算得到的功函數分別為3.051,2.910和3.025 eV,這說明交換關聯勢對功函數依然有較大影響.

圖6 DFT中不同LDA交換關聯勢對基態密度的影響Fig. 6 Ground State Density Calculated for Different LDA Exchange-Correlation Potentials in DFT

圖7 DFT中不同LDA交換關聯勢對功函數的影響Fig. 7 Work Function Calculated for Different LDA Exchange-Correlation Potentials in DFT

1.4 功函數誤差和基態密度誤差

筆者計算了不同動能泛函的基態密度誤差[4],計算公式為

(7)

其中S表示任意選擇的有限表面積,方程右側分母表示有限體積中的電子總數.(7)式中的基態密度誤差可以看作是單個電子的平均基態密度誤差.功函數誤差

其中i=1,2,3,代表3種不同的相關勢,即WG,PZ,GL相關勢.需要注意的是,在對比中要保持選擇的交換關聯勢一致.下面筆者將重點比較OF-DFT中不同動能泛函的功函數誤差與基態密度誤差,探尋金屬板鈉體系的最優能量泛函.

2 誤差最小的能量泛函

圖8 TFvW動能泛函中不同vW動能泛函比重對基態 密度誤差及功函數誤差的影響Fig. 8 Ground State Density Error and Work Function Error Obtained by Different vW Kinetic Energy Functional Proportions

不同vW動能泛函比重下的基態密度誤差與功函數誤差計算結果如圖8所示.圖8中包含了LDA-WG,LDA-PZ和LDA-GL 3種交換關聯勢,點線旁的數字代表vW動能泛函的比重.例如:LDA-WG交換關聯勢中,數字1.0代表vW動能泛函的比重為1.0,該點橫坐標值為0.024 9,縱坐標值為0.487,分別對應這一能量泛函下的基態密度誤差和功函數誤差.圖8中,vW動能泛函的比重由0.3增加至1.0,相鄰點之間的差距為0.1.

對于LDA-WG交換關聯勢,隨著vW動能泛函比重的增加,基態密度誤差先減小后增加,由0.017 5(vW動能泛函比重為0.3時)減小至最小值0.016 4(vW動能泛函比重為0.4時),最后上升至0.024 9.使用LDA-PZ和LDA-GL交換關聯勢計算得到的基態密度誤差趨勢與LDA-WG交換關聯勢的基本一致.其中LDA-PZ交換關聯勢的基態密度誤差由0.017 9下降至0.016 8,最后上升至0.025 4;LDA-GL交換關聯勢的基態密度誤差由0.018 3下降至0.017 1,最后上升至0.025 8.3種交換關聯勢的基態密度誤差最小值點均出現在vW動能泛函比重為0.4時.對于功函數,使用3種交換關聯勢計算得到的功函數誤差同樣隨vW動能泛函比重的增加,呈現出先減小后增加的趨勢.使用LDA-WG,LDA-PZ和LDA-GL交換關聯勢計算得到的功函數誤差最小值分別為0.035,0.051,0.043 eV,分別出現在vW動能泛函比重為0.4,0.5,0.5時.當vW動能泛函比重為1時,使用LDA-WG,LDA-PZ和LDA-GL交換關聯勢計算得到的功函數誤差分別上升至0.487,0.467,0.456 eV.

以上結果說明,TFvW動能泛函中,vW動能泛函比重對基態密度誤差和功函數誤差均有較大的影響.原因主要是,隨著vW動能泛函比重增加,OF-DFT計算得到的功函數增加,基態密度的衰減速度降低.此外,盡管交換關聯勢對功函數有較大影響,比如當vW動能泛函比重為1時,LDA-WG,LDA-GL和LAD-PZ交換關聯勢的功函數分別為3.538,3.478,3.368 eV,但對于相同的動能泛函,使用3種交換關聯勢計算得到的功函數誤差卻相差不大,比如當vW動能泛函比重為1.0時,LDA-WG,LDA-PZ和LDA-GL的功函數誤差分別為0.487,0.467,0.456 eV.原因在于,作為參考值的DFT計算結果同樣受交換關聯勢的較大影響.研究結果進一步說明,vW動能泛函比重為0.4的TFvW能夠適配不同的LDA交換關聯勢.

3 結論

本研究以DFT的計算結果作為參考值,計算了OF-DFT在金屬板鈉體系中獲得準確基態密度和功函數所需的能量泛函.在DFT和OF-DFT中,不同交換關聯勢的使用對基態密度的影響較小,而對功函數有較大影響,這是因為交換關聯勢主要在遠場起作用.此外,OF-DFT中,vW動能泛函比重同時影響基態密度和功函數結果,因此找到合適的vW動能泛函比重至關重要.筆者定義了OF-DFT的基態密度誤差與功函數誤差,通過數值實驗發現,vW動能泛函比重取0.4時,利用不同LDA交換關聯勢計算得到的結果具有較短高的準確度.值得注意的是,對于功函數,作為參考值的DFT計算結果受到交換關聯勢的影響,在vW動能泛函比重為0.4時,不同交換關聯勢的基態密度誤差和功函數誤差較小,這說明vW動能泛函比重為0.4的TFvW能夠適配不同的交換關聯勢.因此,vW動能泛函比重為0.4時,OF-DFT能給出與DFT較符合的基態性質.OF-DFT無需計算單粒子波函數,計算高效,能應用到更大尺寸的復雜納米結構中.