CASTEP計算理論總結+實例分析

1、CASTE就算理論總結XBAPRSCASTEP寺點是適合于計算周期性結構，對于非周期性結構一般要將特定的部分作為周期性結構，建 立單位晶胞后方可進行計算。CASTEP十算步驟可以概括為三步：首先建立周期性的目標物質的晶體；其次對建立的結構進行優(yōu)化，這包括體系電子能量的最小化和幾何結構穩(wěn)定化。最后是計算要求的性質， 如電子密度分布(Electron density distribution) ,能帶結構(Band structure)、狀態(tài)密度分布(Density of states)、聲子能譜(Phonon spectrum)、聲子狀態(tài)密度分布 (DOS of phonon),軌道群分布(Or

2、bital populations) 以及光學性質(Optical properties) 等。本文主要將就各個步驟中的計算原理進行闡述， 并結合作者對計算實踐經驗，在文章最后給出了幾個計算事例，以備參考。CASTEP十算總體上是基于 DFT,但實現運算具體理論有：離子實與價電子之間相互作用采用鷹勢來表示；超晶胞的周期性邊界條件；平面波基組描述體系電子波函數；廣泛采用快速 fast Fourier transform (FFT)對體系哈密頓量進行數值化計算；體系電子自恰能量最小化采用迭帶計算的方式；采用最普遍使用的交換-相關泛函實現 DFT的計算，泛函含概了精確形式和屏蔽形式。一，CASTE沖

3、周期性結構計算優(yōu)點與MS中其他計算包不同，非周期性結構在CASTE沖不能進行計算。將晶面或非周期性結構置于一個有限長度空間方盒中，按照周期性結構來處理，周期性空間方盒形狀沒有限制。之所以采用周期性結構原 因在于：依據 Bloch定理，周期性結構中每個電子波函數可以表示為一個波函數與晶體周期部分乘積的 形式。他們可以用以晶體倒易點陣矢量為波矢一系列分離平面波函數來展開。這樣每個電子波函數就是 平面波和，但最主要的是可以極大簡化 Kohn-Sham方程。這樣動能是對角化的，與各種勢函數可以表示 為相應Fourier形式。 、 、 、工 2|k G、V (G-G) V(G-G) V (G-G)C 、

4、= ；.C. «GG ionHxci,k+Gi i,k+G采用周期性結構的另一個優(yōu)點是可以方便計算出原子位移引起的整體能量的變化，在CASTE即引入外力或壓強進行計算是很方便的，可以有效實施幾何結構優(yōu)化和分子動力學的模擬。平面波基組可以直接達 到有效的收斂。計算采用超晶胞結構的一個缺點是對于某些有單點限缺陷結構建立模型時，體系中的單個缺陷將以無限缺陷陣列形式出現，因此在建立人為缺陷時，它們之間的相互距離應該足夠的遠，避免缺陷之間相 互作用影響計算結果。在計算表面結構時，切片模型應當足夠的薄，減小切片間的人為相互作用。CASTE即采用白交換-相關泛函有局域密度近似(LDQ (LDA)、

5、廣義梯度近似(GGA和非定域交換-相關 泛函。CASTE沖提供的唯一定域泛函是 CA-P乙Perdew and Zunger 將Ceperley and Alder 數值化結果 進行了參數擬和。交換-相關泛函的定域表示形式是目前較為準確的一種描述。NameDescriptionReferencePW91 Perdew-Wang generalized-gradient approximation, PW91Perdew and WangPBE Perdew-Burke-Ernzerhof functional, PBEPerdew et al.RPBE Revised Perdew-Burke

6、-Ernzerhof functional, RPBEHammeret al.采用梯度校正的非定域或廣義梯度近似泛函與電子密度梯度d % 和電子密度P都有關，這樣可以同時提高能量和結構預測的準確性，但計算耗時。CASTE即提供的非定域?函有三種：PBE泛函與PW91泛函計算在本質上實際是相同的，但在電子密度變化迅速體系中PBE泛函實用性更好；RPBN特別用來提高DFT描述金屬表面吸附分子能量的泛函，White and Bird 描述了各種梯度校正泛函計算方法，利用廣義梯度近似計算總能量使用平面波基組與定域泛函相比并不直接。包含梯度近似的交換-相關泛函計算時對電子密度數據的精度要求較高，對計算機

7、內存占用會增大。通過采用與平面波基組總能量計算中分 裂交換-相關能量采用一系列空間網格相一致的方法來定義交換-相關勢。平面波基組(Plane wave basis set )Bloch理論表明每個k點處電子波函數都可以展開成離散的平面波基組形式，理論上講這種展開形式包含的平面波數量是無限多的。然而相對于動能較大的情況，動能| k+G|2很小時平面波系數 G+g更重要。Figure 1. Schematic represeintationi of the cutoff energy concept調節(jié)平面波基組，其中包含的平面波 動能小于某個設定的截止能量，如圖所示(球體半徑與截止能量平方根成

8、比例)：總能量計算會因為平面波特定能量 截止而產生誤差，通過增加體系能量 截止數值就可以減小誤差幅度。理論Egt= -G2° fnax上截止能量必須提高到總能量計算 結果達到設定的精確度為止，如果你 在進行關于相穩(wěn)定性的研究，而需要 對比每個相能量的絕對值時，這是一 種推薦計算方法。不過，同一個結構 在低的截止能量下收斂引起的差別 要小于總體能量本身。因此可以選用 合適的平面波基組對幾何結構進行 優(yōu)化或進行分子動力學研究。以上的方法對Brillouin區(qū)取樣收斂測試同樣成立。有限平面波基組的校正采用平面波基組的一個問題是截止能量與基組數量的變化是間斷的，一般而言在k點基組(k-poi

9、ntset)中不同k點對應不同能量截止( cutoffs )時就會產生這種不連續(xù)性。此外，在截止能量不變時， 晶胞形狀和尺寸的變化都會引起平面波基組的間斷。通過采用更加致密的k點基組就可以解決這個問題， 與特定平面波基組相關的加權性也會消除。然而即使在k基組取樣很致密的情況下，這個問題依然存在， 對其近似的解決方法就是引入一個校正因子( correction factor ),利用某個狀態(tài)基組計算使用了無限 數量的k點與實際采用的數量之間的差別來確定。晶體結構在進行幾何優(yōu)化時如果基組不能真正的達到 絕對的收斂，有限基組的糾正就很重要。比如硅的規(guī)范-保守潢勢很“軟”，在平面波基組截止能量是2 0

10、 0 c V時就已經可以得到準確的計算結果了。但如果計算狀態(tài)方程時使用上述截止能量(比如體積 與總能量和壓強都有關系)，能量最小時對應的體積與體系內壓為零時對應體積是不同的。在提高截止 能量和增加k點取樣基礎上重復對狀態(tài)方程的計算，這兩個體積之間的差別會越來越小。此外截止能量低時計算得到的E-V曲線呈現鋸齒狀，提高截止能量計算的曲線連續(xù)而平滑。E-V曲線中出現鋸齒狀的原因在于平面波基組在相同的截止能量時由于晶體點陣常數不同引起的平面波基組數量的間斷。對總能量 進行有限基組的校正，使得我們可以在一個恒定數量基組狀態(tài)下進行計算，即使采用了恒定的截止能量 這個更強制條件也可以糾正計算結果。Milma

11、n等詳細的討論了這種計算方法的細節(jié)。進行這種校正所需要的唯一的參數就是 dEtot/d lnE cut, Etot是體系總能量，Ecut是截止能量。dEtot/d lnE cut的值很好的表示了 能量截止和k點取樣計算收斂性質。當它的數值（每個原子）小于 0.01 eV/atom 時，計算就達到了良好 的收斂精度，對于大多數計算0.1 eV/atom就足夠了。非定域交換-相關泛函基于LDA或GGA勺泛函的Kohn-Sham方程在計算能帶帶隙上存在低估。這對晶體或分子相關性質以及 能量的描述是沒有影響的。然而要理解半導體和絕緣體性質，就必須得到關于電子能帶結構的準確的描 述。DFT能帶帶隙計算誤

12、差可以通過引入經驗“剪刀”校正，相對于價帶而言導帶產生了一個剛性的變化。當實驗提供的能帶帶隙準確時，光學性質計算得到了較為準確的結果。電子結構實驗數據缺乏時采 用“剪刀”工具進行預測性研究或對能帶帶隙調整是不可靠的。關于DFT計算中能帶帶隙問題已經發(fā)展許多技術，但這些技術大多復雜而且很耗時，實際計算中最常用的是屏蔽交換（Sx-LDA）,建立在廣義Kohn-Sham方法基礎上。廣義Kohn-Sham泛函允許我們將總能量交換分布泛函分離為非定域、定域以及屏蔽密度組元。在 CASTEPf算中采用的廣義 Kohn-Sham方法有： HF: exact exchange, no correlation

13、HF-LDA: exact exchange, LDA correlation, sX: screened exchange, no correlation sX-LDA: screened exchange, LDA correlation與LDA和GGAff比No local functionals 也有一些缺陷。在屏蔽交換泛函中不存在已知形式應力張量表 達方式，因此沒有完全的非定域勢可以用于單位晶胞結構優(yōu)化或進行NPT/NPH力學。這樣利用這些泛函計算的光學性質很有可能是不準確的。在哈密頓量中引入一個完全非定域組元就可以解決這個問題， 這個額外的矩陣元破壞了光學矩陣元素由位置算符轉換為動

14、量算符常用表達形式，使得哈密頓量對易很 復雜。規(guī)范保守潢勢和超軟潢勢震勢是利用平面波基組計算體系總能量中關鍵的一個概念，價電子與離子實之間強烈的庫侖勢用全 勢表示時由于力的長程作用很難準確的用少量的Fourier變換組元表示。解決這個問題的另一種方法從體系電子的波函數入手，我們將固體看作價電子和離子實的集合體。離子實部分由原子核和緊密結合的 芯電子組成。價電子波函數與離子實波函數滿足正交化條件，全電子DFT理論處理價電子和芯電子時采取等同對待，而在震勢中離子芯電子是被凍結的，因此采用震勢計算固體或分子性質時認為芯電子是不 參與化學成鍵的，在體系結構進行調整時也不涉及到離子的芯電子。為了滿足正交

15、化條件全電子波函數 中的價電子波函數在芯區(qū)劇烈的振蕩，這樣的波函數很難采用一個合適的波矢來表達。在震勢近似中芯 電子和強烈?guī)靵鰟萏娲鸀橐粋€較弱的震勢作用于一系列震波函數。震勢可以用少量的Fourier變換系數來表示。理想的鷹勢在芯電子區(qū)域是沒有駐點的，因此需要平面波矢數量很少。眾所周知的是現在將震 勢與平面波矢相結合對描述化學鍵是很有用的。全離子勢的散射性質可以通過構筑震勢得到重現，價電 子波函數相位變化與芯電子角動量成分有關，因此震勢的散射性質就與軌道角動量是相關的。震勢最普 遍表達方式是：Vnl =三 | lm> Vi <lm|where | lm> are the sp

16、herical harmonics and V i is the Pseudopotential for angular momentuml .在不同角動量通道均采用同一個震勢值稱為定域震勢(Local Pseudopotential ),定域震勢計算效率更高，一些元素采用定域震勢就可以達到準確描述。震勢的硬度( hardness)在震勢的應用中是一個 重要的概念，當一個震勢可以用很少的Fourier變換組元就可以準確描述時稱為“軟鷹勢”，硬震勢與此相反。早期發(fā)展的準確規(guī)范保守潢勢很快就發(fā)現在過渡元素和第一周期元素(C、N、O,等)中的描述十分“硬”，提高規(guī)范保守潢勢收斂性質的各種方法都已經被

17、提出，在CASTE即采用了由Lin等提出的動能優(yōu)化而來的規(guī)范保守潢勢。Vanderbilt提出了另一種更基本的方法，放寬規(guī)范保守潢勢的要求，從而生成更軟的震勢。在超軟震勢方法中，芯電子區(qū)的震平面波函數可以盡可能的“軟”，這樣截止能量 就可以大幅度的減少。超軟潢勢與規(guī)范保守潢勢相比除了 “更軟”以外還有其它的優(yōu)點，在一系列預先 設定的能量范圍內遺傳算法確保了良好的散射性質，從而使震勢獲得更好變換性和準確性。超軟震勢通 常將外部芯區(qū)按照價層處理，每個角動量通道中的占據態(tài)都包含了復合矢。這樣就增加了震勢的變換性 和準確性，但同時是以消耗計算效率為代價的?？赊D移性是震勢的主要優(yōu)點。震勢是通過孤立的原子

18、或 離子特定的電子排部狀態(tài)下構建的，因此可以準確的描述原子在那些特定排部下芯區(qū)的散射性質。在相 應條件下產生的震勢可以用于各種原子電子排部狀態(tài)以及各種各樣的固體中，同樣也確保了在不同的能 量范圍內具有正確的散射狀態(tài)。Milman給出了不同化學環(huán)境和一系列結構中采用震勢描述準確性事例。非定域震勢即使在最有效離散表示情況下，體系能量震勢計算依然占用了大量計算時間。此外，在倒易 點陣空間采用非定域震勢會因原子數目增多而耗時以原子立方數增大，因此對于大體系是很適用的。M 勢非定域性是指只有在超過原子芯區(qū)時它才會擴展，由于芯區(qū)是很小的，特別是當體系包含有許多的真 空腔體時，在實空間采用震勢來計算就有很大

19、的優(yōu)勢。這時計算量隨體系中原子數目平方增長，因此是 很適合大體系計算的。將電子劃分為芯電子和價電子在處理交換-相關相互作用時會產生新問題，在原子芯區(qū)兩個亞體系疊加在震勢產生過程中很難完全去屏蔽。在震勢能量算符中與電子密度存在非線形關系 的項就是交換-相關能。Louie等采用了一種簡單的方法來處理芯電子和價電子密度之間非線性的交換-相關能。這種方法在很大程度上提高了震勢的可變換性，特別是自旋極化的計算更為準確。當準芯區(qū)電 子不能簡單處理為價電子時非線性核校正就很重要。另一方面將他們簡單地包含在價層亞體系中從本質 上可以避免 NLCCM理的必要性。規(guī)范保守潢勢采用鷹勢計算關鍵在于可以有效的對化學鍵

20、的價電子進行可再現的近似，震勢與全勢在超過離子實半徑 以后具有完全相同的函數形式。of the all-electron and pseudized wave functionsFigure 1. Schematic representation and potentials兩個函數平方幅度的積分數值應該 是相同的，這等同于要求震勢波函數具 有規(guī)范-保守性，比如每個震波函數只能 描述一個電子的行為。這樣的條件就確 保了潢勢可以再現正確的散射 (Scattering Properties) 性質。生成震勢的典型方法如下所述：選 擇某個特定的電子排部狀態(tài)(不一定就是基態(tài))全部電子計算在一個孤立的原

21、子中進行。從而得到原子價電子能量本征值和價電子波函數。選 擇一個離子震勢或震波函數參數形式，通過對參數的調節(jié)，使得鷹原子計算和全電子原子震勢計算采用 相同的交換-相關勢，在超過截止半徑后與價電子波函數形式相同，震勢的本征值等于價電子的本征值。如果電子波函數和震勢波函數滿足正交歸一，兩者在截止半徑以外的匹配性決定了規(guī)范-保守條件自動成立。離子震勢的截止半徑是實際物理芯區(qū)的二到三倍。截止半徑越小，震勢越“硬”而適用性(transferability )好。計算精度和效率決定了實際中采用的截止半徑的大小。規(guī)范-保守潢勢優(yōu)化在固體計算中依據能量的截止存在一系列優(yōu)化震勢的方法，Lin基于Rappe早期工

22、作提出了下列震勢產生方法：在截止半徑(cutoff radius )內，震勢波函數可以表達為：PS:(r)4、a.jl(q.r) jl(qiRc), '(Rc) i=1 i l i ,jc(qiRC) l(RC)jl(qir)是球形Bessel函數，在r=0和r=Rc之間有(i-1 )個零點。為保證潢勢的實用性，截止半徑越大越好。超過截止矢量qc對動能最小化可得到系數c(i。A匚J Ek(qi ,qcPS(qqc八- dr PS, (r)' 2 1Ps(r)- dqq20 ll 0在第一個方程中讓 qc等于q4。其他的三個限制條件使得震波函數在進行Lagrange連乘(Lagr

23、angemultipliers )時保持正交化(normalization ),并且使震波函數在Rc處的第一個二介偏微分是連續(xù)的。半徑相關Kohn-Sham方程反轉標準步驟產生的一個具有理想收斂性質的平滑震勢函數。Lee提出了進一步改進的方法，在CASTE啜據庫中固體規(guī)范保守潢勢就是采用他的思想設計的。這種通用的方法消除了在特定的截止半徑處震波函數的二介偏微分必須是連續(xù)的條件，因為它是自動滿足這個條件的。這樣對于特定截止半徑Rc允許我們通過調節(jié) qc提高震勢的精度和計算效率。超軟潢勢(ULTRASOFT PSEDUPOTENTIAL為了能夠使平面波基組計算中所采用的截止能量盡可能的小，Vand

24、erbilt提出了超軟震勢方法。 眾所周知規(guī)范-保守鷹勢在收斂優(yōu)化中存在本身缺陷，所以就設計了另一種方法。超軟震勢基礎是在大多數情 況下只有當緊密結合原子價軌道加權性分數大部分在芯區(qū)時，利用平面波基組計算才要求較高的截止能量。在這種情況下，減少平面波基組的唯一方法就是解除( violate )規(guī)范-保守潢勢成立條件，將這些 軌道中的電子從芯區(qū)移去。芯區(qū)的震勢就可以盡可能的“軟”，從而使截止能量降低達到要求。從技術上講，通過引入一個廣義的正交歸一化條件就可以完成。為了覆蓋全部電子電荷，在芯區(qū)對由電子波函數模平方產生的電子密度進行適度放大(augmented)。電子密度劃分成兩部分：擴展在整個晶體

25、中“軟”部分和定域在芯區(qū)的“硬”部分。固體中超軟震勢公式超軟震勢中總能量與采用其他震勢平面波方法時相同，非定域勢Vnl表達如下:Nlnm , IDnm)投影算符3和系數D9分別表征鷹勢和原子種類的差別, 度可以表不為：指數I對應于一個原子位置?？偰芰坑秒娮用躰(r)=工卜i(r)|2+ 工 QIm(r)bi nn(mm)i 1 nm,I ' 八 1O超軟震勢完全由定域部分,是波函數，Q(r)是嚴格位于芯區(qū)的附加函數(Augment function)Voci0n (r)和系數D(0), Q and ，這些變量計算方法在下文中將做介紹。引入一個廣義正交歸一條件來解除規(guī)范-保守潢勢的限制條

26、件:S 是哈密頓重疊算符(Hermitian overlap operator )znm,Iqnm系數q是通過對Q (r)積分得到，超軟震勢的Kohn-Sham方程可以寫為:H代表了動能和定域勢能之和，如下所示:H=T+Vefr+nmj在Vf中包含離子定域勢 Vocion(r), Hartree勢和交換-相關勢等項。通過定義一些新參數就可以將因附加 (augmented)電子密度而產生所有項全部包含在潢勢的非定域部分。IDnm(0)D nmI drVeff(r)Qnm(r)S,波函數與D有關而且事實上投影算符函數 3 (projector functionVadr)1.n(n-T-Voc)*以

27、及它們矩陣內積(inner productsnm= 4n 'm /這樣就可以定義用于固體計算的變量(V iocion(r), D ,Q and )Qnm(r): n(r) mn(r) m(r);qnm= drCnm(r)與規(guī)范-保守潢勢對比，不同之處在于在超軟震勢中存在重疊算符)數量要比規(guī)范-保守潢勢中大兩倍多。與附加(augmented)電荷相關的一系列計算可以在實空間(real space) 中進行，這與函數中定域勢的性質有關。多余的步驟不會對計算效率產生較大的影響。在Laasonen文獻中提供了超軟震勢計算的詳細方法以及總能量微分表達式。震勢生成：與規(guī)范-保守震勢情況一樣，在自由

28、原子上對所有的電子進行計算，得到屏蔽原子勢(screened atomic potential )。每個角動量選擇一系列的參考能量專，一般兩個能量參考點就足夠了。這些能量參考范圍必須包含良好散射性質，在每個參考能量處求解與半徑相關的Kohn-Sham方程，得到規(guī)則初始點。選擇截止半徑，對上面產生的每個全電子波函數構筑一個鷹勢4看/捫融、,蚓(R處與颯陽I M逗1:3罰3罰蛇盲R最后就形成了定域軌道(超過R時就消失)：。門)=(B (B bnm，m) m采用去屏蔽(descreening procedure )方法計算 Vocion(r),D系數：ionV.( r ) = V ( r ) - V

29、u ( r ) - V ( r )loclocHxcD nm)= B nm m q nm 一 drV loc ( r ) n ( r )在去屏蔽方法中可以引入非線性核校正方法(The nonlinear core correction (NLCC) ),這與規(guī)范-保守潢勢中所采用的方法完全一致。在以下情況下超軟潢勢是很適用的： 潢本征值與所有電子本征值相同，在芯半徑截止區(qū)以外震軌道波函數與價電子波函數匹配一致；對于每 個參考能量散射性質都是正確的，這樣通過增加參考能量點數目就可以系統(tǒng)的提高震勢的適用性；在參 考電子排部情況下，震勢價電子密度與全價電子密度相同。關于非線性核校正Louie等人第一

30、次提出了非線性芯校正，使得震勢對磁系統(tǒng)的描述更準確。然而，對于非自旋極化體系中準芯區(qū)電子,NLCC1具有同樣的作用。 DFT總能量準確表達需要 NLCC如下:Etot E. » E 匕)- Eioneexc在震勢的計算式中，電子密度分別來自于芯區(qū)電子和價電子。將芯區(qū)能量假設為一常數并切不計入計算。用一個價電子密度和由震勢計算得到的離子定域勢Eon來代替總電子密度， 這樣芯區(qū)電子與價電子之間所有的相互作用全部轉移到震勢上。由此可以推斷電子密度線性化只是對動能和簡單非線性交換-相關能的一個近似，很明顯當芯區(qū)電子和價電子在空間很好分離時是一個良好的近似。但如果兩個區(qū)域電子密度 的疊加密切時

31、，計算體系本身就會產生錯誤，進一步減弱震勢實用性。解決NLCC問題的方法就是調節(jié)震勢生成方法以及在固體中計算方法。在產生震勢時每個角動量通道對應一個屏蔽勢，并且滿足一定的條件，比如規(guī)范-保守，震波函數本征值與全電子波函數本征值相同等。這些屏蔽勢(screened potentials )對應的原子震波函數(atomic pseudowavefunctions )僅表示價電子。從這些波函數可以得到價層震電 子密度(pseudo charge density ),通過對勢的屏蔽得到"光禿”離子勢( bare ):Vion(r)=Vl(r)-Vee(r)-Vxc <(r)由于交換-相

32、關勢泛函是電子密度的非線性函數，對自旋極化體系采用這種方法產生的離子勢與價電子排 列有關。Louie等提出了將上面方程替換為如下表達：ViOn(r)=Vl(r)-Vee (r)-Vxc1 (r)(r)在屏蔽原子勢中減去總交換-相關勢。此外，在計算交換-相關勢時芯區(qū)電荷必須加到價電子中去，這個 額外原子狀態(tài)信息傳遞給 CASTEP在所有計算中芯區(qū)電荷認為(deemed)是相同的，這種做法的一個缺 點是在利用潢勢計算時芯區(qū)電荷很難準確的用Fourier網格表示。而且通常芯區(qū)電子密度比價電子密度大，這很容易將與價電子密度有關的影響掩蓋掉。以下部分將對部分芯區(qū)校正方程建立做介紹，該方法 充分的認識到價

33、電子與芯電子密度重疊的區(qū)域才是我們感興趣的?？拷雍说男倦娮用芏炔粫a生物理結果，雖然有如上所述的一些問題。部分NLCC采用一個在特定半徑以外與 口一致的函數替代全芯電子密度，在原子核周圍這個函數起伏是平滑的。在CASTE呻對一些特定元素在震勢中采用的部分芯區(qū)校正使用了數值化的芯區(qū)電子密度。在規(guī)范保守震勢中雖然有相關的內容，但在計算中并沒有采用這個方法。A Introduction to DFT第一性原理(The first principle )計算也稱為從頭算起(ab-initial calculation ),由于固體的許多基本的物理性質是由其微觀的電子結構決定的，因此通過求解多粒子系

34、統(tǒng)的Schodinger方程，來獲取固體全部的微觀信息從而預測宏觀的性質。利用這個思想建立的能量的哈密頓量非相對論形式可以表 示如下：匕14reirn1n2,rn，-2 j mjk%,*,.e,in1,n2,付(d22zjeZ1Z25+ 2 jp于“團邊監(jiān).國屁但；.rn)12 renej i,j rernj Ju2 rl i1i2j1j2書網(ie1,e2圖.eirn11n2.rnj)考慮到原子核與核外電子質量差別以及電子馳豫時間比原子核馳豫時間要小三個數量級，因此利用 Born-Oppenheimer近似將原子核運動和電子的運動分離，從而將體系波函數劃分為電子波函數和原子核波函數兩個部分，

35、分別用中 和4表示：%(一,r c, rc,rr J c L) = V (L,r c,rc,r .),(rJ C r,), , ,ik el e2 e3 e, n1,n2, njel e2 e3 ei n1,n2, nj能量的哈密頓量可以分解為如下的兩個方程:(一口甲+ 工-2L2me i 12 卜eTej 1 i2+ £i,jzje2rei-rnj廣(re"2，re3,，r、ei)= Eke%r r e2, e3,ei)n,r .) 二 E| (r . r _nj kn n1, n2,，rnj)：2222(-$"+£三+£j"jre

36、f|情人電n1”第一性原理嚴格求解僅在氫分子中實現了，對于多粒子體系的計算幾乎是不可能的。目前均采用不同的近似方法來實現計算，主要方法有Hartree-Fock 近似和DFT近似。在Hartree-Fock 近似中體系的哈密頓量表示如下：E TotaliHF - j、i' (Jij Kij )箕 為第i個電子的Hartree-Fock的軌道能，Jij是庫侖積分，表示電子靜電互斥能，Kij為交換積分。交換積分所代表的交換能指電子由于自旋平行而引起的電子軌道庫侖能量減少的部分。密度泛函理論(Density Functional Theory)建立了將多電子體系化為單電子方程的理論基礎，并且

37、給出了有效勢計算方法，是目前研究多粒子體系性質的一種普遍使用的重要方法。該理論認為對于處于外勢場V (r)中相互作用的多電子系統(tǒng)，電子密度分布函數p (r)是決定該系統(tǒng)基態(tài)物理性質的基本變量。密度泛函理論中體系的能量泛函表示如下：Et(：) = T(：) U(：)Excenergy; E xc(：):exchangeande2Et，(r)(>()力口(仁：(r)：()r - rdrdrExc：(r)T (丹 ： Kinetic energy; U ( D :classical electrostatic correlation energy由上表達式可見體系能量是電子密度的泛函，因此可以

38、進一步將上式表達為:在上式中第一項為電子在外場中的勢能，第二項為電子的動能，第三項為電子相互之間的庫侖能，第四 項為交換相關能，最后一項形式是未知的。系統(tǒng)的電子密度分布可以表達如下：利用上式可以將動能項表示為:：()i (r)US表達為:n NUP(r)=乙£ %(r)i a-zFa-ri(r)1 _Z2i' i(r1) j(r2i(1) j(r2)j ；Tr2N zaz a ： aRa-RN1r2a a"aRrR=-” ：(口) a= _：(r)VN：()# vnnExc( 形式確定有兩種方法： 局域密度近似(LDA,Local Density Approxima

39、tion )和廣義梯度近似(GGA,General Gradient Approximation )。在局域密度近似(LDA)中采用了均勻電子氣的分布函數推倒出了非均勻電子氣的交換-相關能泛函，從而彳#到Exc(P)的具體形式。從近期計算結果相關報道來看采用局域密度近似(LDQ計算在絕緣體中會產生較大的誤差，而且對帶隙寬的半導體等得到不正確的結果。采用局域密度近似(LDA)主要的缺陷現歸納如下：1 .對光學躍遷帶隙預測很差(一般是過低估計帶隙寬度)。這雖然對基態(tài)性質如電荷密度，總能量以及力影響不大，但在導帶狀態(tài)計算中卻是個大問題，如關于光學性質，運輸性質等的計算。在諸如光伏 裝置等領域的研究中

40、，帶隙就是個很重要的問題。采用“剪刀”(Scissors)工具在固體帶隙計算中很有用，但對我們未獲得實驗結果的物質，是不能采用這個方法的。2 .對類似于二氧化硅這樣的電子氣分布極不均勻體系，基本假設中關于電子密度分布在空間是緩慢變化的條件是不滿足的，這樣的體系采用LDA處理就存在難題。3 . LDA簡單的認為計算體系是順磁性(Paramagnetic)的，對于包含未配對(Unpaired)自旋體系采用局域自旋密度近似(LSDQ (對自旋向上(spin up)和向下(spin down)的電子分別采用密度泛函計算)是 很有用的，比如費米能級 (Fermi level)處半填充的系統(tǒng)。4 .最后一

41、個很少關注的領域就是玻璃陶瓷工業(yè)，LDA對弱的結合鍵(如偶極漲落)很難描述，氫鍵(Hydrogen bond)在LDA中也無法獲得準確的計算結果。GGAM以則改進了 L(S)DA,將相關交換能確定為電子密度極其梯度的函數，在GGA派中以Perdew等人Kohn-Sham 方程:認為交換相關能的泛函形式應該以一定的物理規(guī)律為基礎，構造了著名的PBE泛函。將電子密度分布函數帶入體系能量電子密度泛函中，對泛函變分求極小值，可以得到Ei”i2F 二 一' V ； (r)KST% ：()= V ( rKST：(r ) d ri ExcT Ti6P ( r )r - r交換-相關能可以按照下式計算

42、:E = (r)E ：(r)drxcxcnumberof particles;xcr):exchange-correlationenergy perparticles in an uniformelectron gas ;：(r):distribution function of electron density.Exc-56P( r )稱為交換-相關勢和，表示為:口 xC ：?(r) xcExc' '(r )在Castep計算中采用了周期性邊界條件，單電子的軌道波函數滿足Bloch定理，采用平面波展開式有:、"(r)二 e""i(r)周期性邊界條

43、件下的波函數擴展為一系列分離的平面波波矢，這些波矢與晶體的倒易點陣矢量相聯系。iG *RiOCi,Ge2.2晶體光學性質的計算基于以下原理：電磁波在真空以及某種材料介質中傳播時差別可以用一個復數式的折射指數來表示:n ik在真空中N為實數，而且其大小為1;在其他介質中時若材料對于光是透明的則是一個純實數，虛部對應材料的吸收系數(Adsorption Coefficient )。它們之間的關系方程 2所示：2k吸收系數表示的是電磁波通過單位厚度的材料時能量的衰減分數，通?？梢杂貌牧辖苟鸁岬漠a生來衡量。反射系數(Reflection Coefficient )可以簡單通過將垂直光束照射材料的表面引

44、起1N(n-12 k2(n 12 k2在計算光學性質時一般先計算虛部的介電常數，其他的性質與介電常數之間建立關系。虛部介電常數計 算式由下方程確定：這樣折射指數的實部和虛部以及介電常數之間的關系可以寫為:2nk光導率(Optical conductivity )也是一個普遍用來描述材料光學性質的物理量。光導率的表達式為方 程7:這個參數用來描述金屬的光學性質，但在CASTE即將計算范圍擴大到了絕緣體和半導體。計算過程的主要的區(qū)別在于前者的光學譜中IR部分與內部能帶之間的轉變密切相關，而者則在計算內時并沒有完全考慮到這些因素。從虛部介電常數可以進一步得到材料電子的能量損失函數(Energy Lo

45、ss Function) ,它描述了電子通過均勻的電介質時能量的損失情況，計算式如下所示：-1Im（“）8在實驗中我們可以測定的光學性質參數有吸收系數”（0）和反射系數R（o ）。從理論上而言，得到這些參數以后可以將方程 2、3、4表示為復數的形式之后得到表達式1中的實數部和虛數部。但在實際情況下由于入射光源的復雜性，而且晶體結構中極化效應使得材料介電常數并非是各向同性的。此外材料表面幾何結構也不是理想的平滑表面。這些因素就限制了對其光學參數的預測。在CASTE呻提供的光學性質的計算支持體系極化，但狀態(tài)只能在同種自旋間相互轉化。晶體中聲子和電子之間的相互作用可以用電子基態(tài)波函數中包含的含時微擾

46、項來表示，聲子電場擾動引起了電子函數占據態(tài)和未占據態(tài)間的轉變（磁場引起的效應要弱一個因數V/C）,這些激發(fā)態(tài)（激子）聚集態(tài)稱為等離波子。單獨的態(tài)激發(fā)稱為單粒子激子，這些激子對光譜產生的結果是導帶和價帶的狀態(tài)密度之間的連接可以通過選擇合適的加權性矩陣元素來實現。在CASTER部介電常數的方t算按照方程9進u矢量定義光束電場的極化性質。這個表達類似于含時微擾的Fermi-Golden定理，比）可看作真實占據態(tài)與未占據態(tài)之間轉換的細節(jié)。介電常數就描述了一種因果效應，它的實數部和虛數部之間由Kramers-Kronig 變換相聯系。利用這個變換就可以得到介電常數的實數部?（切）。用于描述電子態(tài)轉變的位

47、置算符矩陣元素通常用動量算符矩陣元素來表示，這樣可以在倒易點陣空 間直接的進行計算。局域勢函數會影響計算，在CASTEP十算中一般采用非定域勢函數。本文在進行BFGS晶體結構幾何優(yōu)化時就選擇了非局域勢函數。經過矯正后的矩陣元素可以描述如下：平沖力=上Qk中力+ "限%屈*”10利用超軟鷹勢（Ultra soft Pseudopotential ）計算時會增加額外矩陣元素，在目前CASTEP十算中這部分矩陣元素并沒有涉及。采用規(guī)范保守勢計算結果發(fā)現與采用超軟震勢計算符合的很好，因此額外 的那部分矩陣元素對于計算結果的影響不大。晶體光學性質IR部分受能帶內部的影響較大，采用經驗Drude

48、表達形式就可以精確地描述這個影響。%（十寸1-TDDrude校正的光導率。和Drude限制系數 為與材料許多實際參數有關，一般這些參數可以通過實驗得到。結合上式和式 7就可以了解介電函數中Drude的貢獻，同樣可以得到在其它光學常數中的分布。Drude限制參數描述了計算過程中未涉及因素引起光譜寬化現象，比如電子間的散射效應(包括Auger效應)、電子與聲子之間的散射效應以及電子與晶體結構缺陷之間的散射效應等。在CASTE沖光學性質計算結果的準確性與下列因素有關：1 .導帶數量(Number of conduction bands ):直接決定了 Kramers-Kronig 變換的準確性。2

49、.截止能量(Energy cutoff):體系能量進行迭代計算過程中，電子基態(tài)能量本征值精度直接影響能帶 結構以及光學性質，提高截止能量的數值可以提高計算精度，可以得到更準確未占據態(tài)的自恰電荷密度 和震動自由度。3 .迭代計算中 K點數量(Number of k-points in the SCF calculation):與截止能量對體系基態(tài)能量計算影響一樣，K點數量越多，迭代計算能量越準確。4 .積分 Brillouin zone K 點數量(Number of k-points for Brillouin zone integration): 在計算光學性質矩陣元素時Brillouin

50、zone 選取的K點數量應當是合適的，與電子能量相比，矩陣元素在Brillouin zone變化更快，因此必須選取足夠數量的K點來提高矩陣元素計算結果的準確性。從目前計算結果對比來看，提高上述參數的準確性時，光譜中特征峰可以快速地達到實際的要求。當然CASTE呻對光學性質的計算還有不少的局限性，電介質極化引起的局域場效應在現在計算中被忽略了，這對光譜計算有一定的影響，但在目前計算方式下將是無法進行的。準粒子和DFT能帶帶隙以及激子等都會影響計算結果。 狀態(tài)密度在Brillouin zone 區(qū)的表示：給定能帶n對應的狀態(tài)密度 Nn (E)定義為：Nn(k) =dk3,4(E- En(k)En(

51、k)描述了特定的能帶分布情況，積分在整個Brillouin zone 進行。另外一種表示狀態(tài)密度的方法基于N ( E)d E與第N級能帶在能量E到E+dE范圍內允許波矢量數成比例。總體狀態(tài)密度N (E)就是對所有的能帶允許電子波矢量求和，從能帶極小值積分到費米能級就得到了晶體中包含的所有的電子數。在 自旋極化體系中狀態(tài)密度可以用向上自旋(多數自旋(majority spin)和向下自旋(少數自旋(minorityspin)分別進行計算，他們的和就是整體狀態(tài)密度分布，它們的差值稱為自旋狀態(tài)密度分布。借助于狀 態(tài)密度這個數學概念可以直接對電子能量分布進行積分而避免了對整個Brillouin zon

52、e 積分。狀態(tài)密度分布經常用于快速直觀的分析晶體的電子能帶結構，比如價帶寬度、絕緣體中能隙以及主要特征譜峰強 度分析，這對于解釋實驗各種譜數據有很大的幫助。狀態(tài)密度還可以了解當晶體外部環(huán)境如壓力等發(fā)生 變化時電子能帶的變化情況。狀態(tài)密度數值化計算方法很多，最簡單的方法是對各個能帶電子能級進行采用柱狀圖取樣Gaussian擬和。用這種方法繪制的狀態(tài)密度分布圖不存在類似于van-Hove奇點尖銳分布，但只需要少量的K點即可。其他的準確方法基于對Brillouin zone參考點之間采用線形或二次方內叉法。目前最可靠和普遍使用的方法是四面體叉入法，但這種方法與Brillouin zone網格特殊點是

53、不融合的。因此CASTEP!用了由Ackland發(fā)展的簡單的線性內叉法，對 Monkhorst-Pack倒易基組平行六面體采用線性內叉法，能帶能 量組合基組進行柱狀取樣。2.4 偏態(tài)密度(PDOS和局域X態(tài)密度(LDOS)偏態(tài)密度(PDOS和局域狀態(tài)密度是一種分析電子能帶結構有效的半經驗方法。局域狀態(tài)密度表示了 體系中不同原子在各個能譜范圍內電子狀態(tài)分布情況。偏態(tài)密度(PDOS進一步將上述分布以角動量貢獻進行量化分析。了解狀態(tài)密度分布峰值中S、P和D軌道貢獻是很有用的。LDO解口 PDOS供了一種定量分析電子雜化狀態(tài)的方法，對于解釋XPS和光譜峰值的起源很有幫助。PDOS計算基于 Mullik

54、en population 分析，每個給定原子軌道在能帶各個能量范圍內分布均表示出來，特定原子所有軌道的狀態(tài)密度分布和以LDOSI示出來。與整體態(tài)密度計算相似，采用了高斯混合算法或線形內叉法。 Brillouin zone 積分取樣大快固體中電子狀態(tài)只允許存在于由邊界條件確定一系列k點中，固體周期性結構中包含了無限數量的電子，這對應于無限數量的k點。無限數目的電子波函數計算利用Bloch定理轉變?yōu)橛糜邢迶盗縦點計算有限數量的波函數。每個k點處電子占據態(tài)都會對電子勢有貢獻，因此在理論上要進行無限數量的計算。對于十分臨近的k點，它們的電子波函數幾乎是完全相同的，因此在DFT表達中對所有k點求和(等

55、價于對整個Brillouin zone 積分)可以采用有效的離散化數值計算，即在Brillouin zone 選取有限數量的特殊點。進一步考慮到對稱性，只對 Brillouin zone 無法簡并的部分才計入計算過程。Payne以及Srivastava and Weaire等人的文獻提供特殊k點選擇方法以及求和加權的評論。采用上述方法以后，選用很少的k點對絕緣體電子狀態(tài)計算就可以獲得對電子勢和總能量準確的近似。對于金屬體系而言為了得到費米能級準確性，需要更致密的k點數量。采用更多k點數量就可以減小因K點數量限制而產生的對總能量計算的誤差，與獲得基組數量方程收斂方法類似。當對對稱性不同的兩個體系

56、的能量進行對比時，與k點取樣相關的計算收斂精度要更高，例如比較FCC或HCP結構相對穩(wěn)定性。在這種情況下計算誤差是不可避免的，因此能量必須達到絕對收斂精度。要注意的是，體系總能量不會因k點數量的不同而發(fā)生變化，因此即使收斂精度很低時能量計算也一樣，這就與平面波基組截止能量的收斂計算不同，后者平面基組增大時總能量會減少。Monkhorst-Pack 特殊點(special points )Monkhorst-Pack發(fā)展了一種目前普遍采用的特殊k點產生方法，最初只在立方體系中使用，后來Monkhorst-Pack將其進一步擴展到了六方晶格中，在倒易空間沿著坐標軸生成均勻規(guī)則分布的k點網絡。Monkhorst-Pack 網絡采用三個積分來定義，qi where i=1,2,3 ，確定了與主坐標軸之間的偏差。這些積分得到了下面的一些數字：ur=(2 r-qi -1)/2 qiwhere r varies fro