計算材料學-PPT課件

1、計算材料學楊振華第一性原理計算方法 第一性原理方法是一種理想的研究方法，物理學家常稱第一性原理方法，化學家常稱為“從頭算”，但是本質都是一樣的。就是從材料的電子結構出發(fā)，應用量子力學理論，只借助于普朗克常數(shù)h、電子的靜止質量m0、電子電量e、光速c和波爾茲曼常數(shù)k這五個基本的物理常量，以及某些合理的近似而進行計算。這種計算不需要任何其他可調的(經(jīng)驗的或擬合的)參數(shù)就可以如實地求解材料的一些基本物理性能參數(shù)。通過求解多粒子系統(tǒng)總能量的辦法來分析體系的電子結構和原子核構型的關系，從而確定系統(tǒng)的性質 。絕熱近似 波恩(Born M)和奧本海默(Oppenheimer J.E) 提出了絕熱近似，根據(jù)這

2、種近似，可以將原子核運動和電子的運動分開。通過絕熱近似，可以獲得多電子的薛定諤方程 電子作用項 原子核作用項 電子和原子核相互作用項 波恩(Born M)和奧本海默(Oppenheimer J.E) 提出了絕熱近似多電子的薛定諤方程 ，成功地分開了電子的運動與原子核的運動 單粒子算符雙粒子算符哈特利方程 此方程以位于r處的單個電子為研究對象，描述其在晶格勢和其他所有電子的平均勢中的運動規(guī)律 將多電子問題變?yōu)榱藛坞娮訂栴}，但是沒有考慮電子的交換反對稱性 。為了研究電子的交換反對稱性的影響，采用Slater行列式來求能量，經(jīng)過合適的變換，得到了如式所示方程： 單電子的哈特利-?？朔匠?，比哈特利方程

3、多了交換相互作用項。 多電子的薛定諤方程可通過哈利特-?？私坪喕癁閱坞娮佑行Х匠?，如式所示。 包含了電子與電子的交換相互作用，但自旋反平行電子間的排斥相互作用沒有被考慮，即還需考慮電子關聯(lián)相互作用。 為了更加準確地描述多電子系統(tǒng)，Hohenberg P和Kohn W提出了兩個基本的定理：(1) 定理1：不計自旋的全同費密子系統(tǒng)的基態(tài)能量是粒子數(shù)密度函數(shù)的唯一泛函；(2) 定理2：能量泛函在粒子數(shù)不變條件下對正確的粒子數(shù)密度函數(shù)取極小值，并等于基態(tài)能量。定理1的主旨思想是粒子數(shù)密度函數(shù)是一個決定系統(tǒng)基態(tài)物理性質的基本變量；定理2的要點是在粒子數(shù)不變條件下能量泛函對密度函數(shù)的變分就得到系統(tǒng)基態(tài)的

4、能量。密度泛函理論的理論基礎是這兩條基本定理，其基本的思想是原子、分子和固體的基態(tài)物理性質可以用粒子密度函數(shù)來表示。 Hohenberg-Kohn定理說明了粒子數(shù)密度是確定多粒子系統(tǒng)基態(tài)物理性質的基本變量以及能量泛函對粒子數(shù)密度函數(shù)的變分是確定系統(tǒng)基態(tài)的途徑。但是仍然存在三個問題未解決：(1) 如何確定粒子數(shù)密度函數(shù)；(2) 如何確定動能泛函；(3) 如何確定交換關聯(lián)能泛函。為了解決這三個問題，Kohn W與 Sham L.J共同合作，提出了Kohn- Sham方程 。Kohn W和Sham L.J成功地提出了Kohn- Sham方程，用無相互作用的粒子模型代替有相互作用粒子哈密頓量中的相應項

5、，將有相互作用粒子的全部復雜性歸入交換關聯(lián)作用泛函。將多粒子系統(tǒng)的基態(tài)求解轉化為單粒子系統(tǒng)的等效求解，解決第一和第二個問題，對于第三個問題，需要采用局域密度近似來解決。為了求解Kohn- Sham方程，必須構造合適的交換關聯(lián)能。目前比較常用的交換關聯(lián)能主要有以下兩種形式：局域密度近似(LDA)和廣義梯度近似(GGA)。局域密度近似局域密度近似最早是由Kohn W和Sham L.J提出來的，這是一種既簡單可行而又很有效的近似，其基本思想是在局域密度近似中，利用均勻電子氣密度函數(shù)來獲得非均勻電子氣的交換關聯(lián)泛函。 交換關聯(lián)能可以寫為式 Kohn- Sham方程中的交換關聯(lián)勢近似為式 密度為：均勻無

6、相互作用電子氣的交換-關聯(lián)密度，在實際的計算過程中，通常把交換-關聯(lián)密度分成兩部分：交換項和關聯(lián)項。 交換能關聯(lián)能考慮了自旋 Local Density Methods 假設局域電子密度可以被認為是均勻電子氣，或等效地說，電子密度是隨空間緩慢變化的函數(shù)。交換項Local Density Approximation (LDA)Local Spin Density Approximation (LSDA)關聯(lián)項Vosko，Wilk，and Nusair （VWN）GGA （見下）中的 PW91 修改了 VWN 的泛函形式： Gradient Corrected Methods Gradient C

7、orrected or Generalized Gradient Approximation (GGA): 泛函不僅決定于電子密度，還決定于電子密度的梯度。交換項Perdew and Wang (PW86): 修正 LSDA 的泛函形式：加入高階項。Becke (B or B88): 正確的能量密度漸進行為。Becke and Roussel (BR): 加入軌道波函數(shù)的導數(shù)項。Perdew and Wang (PW91)關聯(lián)項Lee, Yang, and Parr （LYP）Perdew：修正 LSDA 的梯度項。Perdew and Wang（PW91 or P91）：。其中 在 LSDA

8、 部分已經(jīng)給出。Becke（B95）：更好地滿足一些基本的物理約束。 混合方法 混合 HF 和 DFT 給出的能量項。 Becke 3 parameter functional (B3)廣義梯度近似 為了對局域密度近似進行提高和改善，引入了電荷密度梯度，即粒子密度的空間分布不僅僅與局域密度有關系，而且與對應點附近的密度有關系。其中最為常用的是廣義梯度近似(GGA)，在GGA近似下，在交換相關能泛函中引入電子密度的梯度來完成。 考慮了電子密度的非局域性，改善了LDA的計算結果。一般GGA的計算結果與實驗結果較為吻合。DFT+U方法簡介 基于密度泛函理論(DFT)的第一性原理計算方法已在材料的晶體

9、結構、磁結構、電子結構以及材料的力學性能計算等方面取得了巨大的成功，但是對于Mott絕緣體(如過渡金屬氧化物和稀土氧化物)，由于其d電子或f電子的強關聯(lián)作用，傳統(tǒng)的第一性原理方法已不能很好地描述其基本性質。 在Mott絕緣體中，當電子從一個一個原子位置跳躍到另外一個原子位置時，如果那個原子位置已經(jīng)擁有一個電子，電子之間就會產(chǎn)生庫倫排斥力作用，這種跳躍需要一定的能量以致能克服這種庫倫斥力作用，如果這個能量大于能帶帶隙，即使能帶沒有全部占滿，電子也很難自由輸運，從而使材料體現(xiàn)絕緣體的特征。 當采用傳統(tǒng)的第一性原理計算Mott絕緣體時，只考慮了交換參數(shù)J，沒有考慮Hubbard參數(shù)U，而在Mott絕

10、緣體中，其決定性的參數(shù)是Hubbard參數(shù)U值，因此采用傳統(tǒng)的計算方法往往會導致失敗。為了解決計算Mott絕緣體的問題，Anisimov等提出了Anisimov 模型，在該模型中，將所研究的電子分為兩個部分：(1) 傳統(tǒng)的DFT算法，在此過程中沒有考慮Hubbard參數(shù)U；(2) 對于d軌道電子或f軌道電子，能帶模型為Hubbard模型，考慮了d軌道或f軌道電子的強關聯(lián)作用。 LDA+U方法為例，電子的總能量計算可以通過下式進行表述： 局域態(tài)的軌道占據(jù)數(shù) 總的局域電子數(shù) 考慮了d軌道或f軌道電子的強關聯(lián)作用，并采用Hartree表達式所計算的能量 原來傳統(tǒng)LDA計算過程所包含的關聯(lián)能，采用LD

11、A+U方法后，此項應該減去 U 為Hubbard參數(shù) 對軌道占據(jù)數(shù)進行微分電子軌道勢當軌道占據(jù)數(shù)分別為1和0時，相應的值表示將采用傳統(tǒng)LDA計算所得的軌道能量分別偏移 VASP計算軟件包簡介VASP，其全稱是Vienna Ab-initio Simulation Package。它基于1989年的CASTEP(1989版)，最早是由Gerorgo Kresse 和Jrgen Furthmller合作，共同開發(fā)出來的，1995年被正式命名為VASP，隨后被開發(fā)者不斷完善。VASP是一種使用贗勢和平面波基組進行從頭量子力學分子動力學計算和第一性原理計算的軟件包，主要用于具有周期性的晶體或表面的計算

12、，可以采用大單胞，也可以用于處理小的分子體系。 與同類的軟件相比，它比較早地實現(xiàn)了超軟贗勢，計算量相對于一般的模守恒贗勢方法大為減少。其對計算領域最大貢獻無疑是在Blchl的基礎上發(fā)展的投影綴加平面波(PAW)方法，這是最重要的。這使得VASP不僅計算速度快，而且精度是abinit和pwscf沒法比的。VASP的精度，比如磁性計算，很多可以跟FLAPW相比，并且計算速度比FLAPW快很多。 在實空間計算勢的非局域部分并保持正交化的數(shù)目減少，使得計算時間小于N3；VASP在電子自洽迭代計算中，采用了RMM-DISS和blocked Davidson等非常有效算法并能自動確定體系的對稱性；此外，V

13、ASP的代碼使用FORTRAN語言編寫，可讀性好，幾乎支持所有的計算機平臺，已廣泛應用于材料科學領域。 VASP基本原理簡介基本知識常用關鍵詞使用說明計算結果處理VASP程序基本原理 VASP是基于贗勢平面波基組的密度泛函程序，其前身是CASTEP 1989版本，其基本原理如下： 根據(jù)Bloch定理，對于周期體系，其電子波函數(shù)可以寫為單胞部分和類波部分的乘積： 其中，單胞部分的波函數(shù)可以用一組在倒易空間的平面波來表示： 這樣，電子波函數(shù)可以寫為平面波的加和： 根據(jù)密度泛函理論，波函數(shù)通過求解KohnSham方程來確定： i：KohnSham本征值Vion：電子與核之間的作用勢VH和VXC：電子

14、的Hartree勢和交換相關勢 基于平面波表示的KohnSham方程： 上式中動能項是對角化的，通過求解上式方括號中的哈密頓矩陣來求解KS方程，該矩陣的大小由截至能(cutoff energy)來決定。嘗試電子密度和嘗試波函數(shù)寫出交換相關勢表達式構造哈密頓量子空間對角化，優(yōu)化迭代自由能的表達式E新電子密度，與嘗試電子密度比較輸出結果，寫波函數(shù)是否程序流程： 與原子軌道基組相比，平面波基組有如下優(yōu)點：無需考慮BSSE校正；平面波基函數(shù)的具體形式不依賴于核的坐標，這樣，一方面，價電子對離子的作用力可以直接用Hellman-Feymann定理得到解析的表達式，計算顯得非常方便，另一方面也使能量的計算

15、在不同的原子構象下具有基本相同的精度；很方便地采用快速傅立葉變換(FFT)技術，使能量、力等的計算在實空間和倒易空間快速轉換，這樣計算盡可能在方便的空間中進行；計算的收斂性和精確性比較容易控制，因為通過截斷能的選擇可以方便控制平面波基組的大小。 平面波基組方法的不足之處：所求得的波函數(shù)很難尋找出一個直觀的物理或化學圖象與化學家習慣的原子軌道的概念相聯(lián)系，即其結果與化學家所感興趣的成鍵和軌道作用圖象很難聯(lián)系出來，這就為我們計算結果的分析帶來了困難；考察某些物理量時，例如原子電荷，涉及到積分范圍的選取，這造成所得物理量的絕對值意義不大；有些方法，例如雜化密度泛函方法不易于采用平面波基組方法實現(xiàn)。V

16、ASP程序基本知識1. VASP程序主要功能：能量計算J. Phys. Chem. C, 2008, 112, 191能帶結構DOS2) 電子結構(能帶結構、DOS、電荷密度分布)電荷密度分布J. Phys. Chem. B, 2005, 109, 192703) 構型優(yōu)化(含過渡態(tài))和反應途徑J. Phys. Chem. B, 2006, 110, 154544) 頻率計算和HREELS能譜模擬J. Phys. Chem. C, 2007, 111, 74375) STM圖像模擬Surf. Sci., 2007, 601, 34886) UPS能譜圖像模擬Surf. Sci., 2007,

17、601, 34887) 材料光學性質計算8) 其它性質計算，包括功函、力學性質等2. 重復平板模型(或層晶模型)： VASP程序采用重復平板模型來模擬零維至三維體系零維分子體系Dv: Vacuum thickness (10 A)二維固體表面說明： 重復平板模型中的平移矢量長度必須合理選擇，以保證：對于分子體系，必須保證相鄰重復單元中最近鄰原子之間的距離必須至少710埃以上；對于一維體系，相鄰兩條鏈最近鄰原子之間的距離必須至少710埃以上；對二維體系，上下兩個平板最近鄰原子之間的距離必須至少710埃以上；4) 嚴格意義上，通過考察體系總能量/能量差值對真空區(qū)大小的收斂情況來確定合理的平移矢量長

18、度。3. K網(wǎng)格大小的選擇： 對于一維至三維體系的計算，需涉及k點數(shù)目的選擇，對于K點的確定，它與布里淵區(qū)的形狀以及對稱性有關。VASP的K點輸入方法有多種，其中最常用的是直接給定K-mesh的大小，然后程序根據(jù)布里淵區(qū)的形狀以及對稱性自動生成各K點的坐標和權重。 對于K-mesh的確定方法，通常通過考察總能量/能量差的收斂程度來確定，能量的收斂標準是1meV/atom。 多數(shù)情況下，對半導體或絕緣體較小的K-mesh能量就可以收斂，對于導體，一般需要較大的K-mesh。硅體相總能量隨K-mesh大小的變化情況4. Cutoff energy大小的選擇： 截至能的大小直接影響到計算結果的精度和

19、計算速度，因此，它是平面波計算方法的一個重要參數(shù)。 理論上截至能越大計算結果也可靠，但截至能大小決定了計算中平面波的數(shù)目，平面波數(shù)目越多計算時間約長、內存開銷越大。 一般根據(jù)所求物理量來確定截至能，例如計算體模量以及彈性系數(shù)時，需要較高的截至能，而通常的構型優(yōu)化只要中等大小的截至能即可，另外動力學模擬時，可選取低的截至能。 不同元素在構造其贗勢時，有各自的截至能，對于VASP，在缺省情況下，選取的是中等大小的截至能，這對于求解多數(shù)物理量是足夠的。嚴格意義上，截至能的確定與K-mesh大小的確定類似，也是通過考察在總能量的收斂情況來確定(即保證總能量收斂至1meV/atom)。硅體相總能量隨cu

20、toff energy大小的變化情況5. VASP輸入和輸出文件：輸入文件(文件名必需大寫) INCAR : 其內容為關鍵詞，確定了計算參數(shù)以及目的；POSCAR : 構型描述文件，主要包括平移矢量、原子類 型和數(shù)目、以及各原子坐標；KPOINTS : K點定義文件，可手動定義和自動產(chǎn)生；POTCAR : 各原子的贗勢定義文件。主要輸出文件 OUTCAR : 最主要的輸出文件，包含了所有重要信息； OSZICAR : 輸出計算過程的能量迭代信息；CONTCAR: 內容為最新一輪的構型(分數(shù)坐標，可用于續(xù)算)；CHGCAR、CHG、PARCHG :用于電荷密度圖繪制；WAVECAR : 波函數(shù)文

21、件；EIGENVAL: 記錄各K點的能量本征值，用于繪制能帶圖；XDATCAR: 構型迭代過程中各輪的構型信息(分數(shù)坐標，用于 動力學模擬)；DOSCAR : 態(tài)密度信息。POSCAR文件內容說明：Silicon bulk (Title) 2.9 (Scaling factor or lattice constant) 0.0 1.0 1.0 (第一個平移矢量的方向) 1.0 0.0 1.0 (第二個平移矢量的方向) 1.0 1.0 0.0 (第三個平移矢量的方向) 2(單胞內原子數(shù)目以及原子種類)Selective dynamics(表示對構型進行部分優(yōu)化，如果沒這行，則表示全優(yōu)化)Dire

22、ct (表示所采用的為分數(shù)坐標，如果內容為Car，則坐標單位為埃) 0.125 0.125 0.125 T T T (各原子坐標以及哪個方向坐標放開優(yōu)化) -0.125 -0.125 -0.125 T T Tsurface of mgo(100) (2*2)Mg 1.00000000000000 5.9459999999999997 0.0000000000000000 0.0000000000000000 0.0000000000000000 5.9459999999999997 0.0000000000000000 0.0000000000000000 0.0000000000000000

23、 20.0000000000000000 20 20 (體系中有2種元素，各自的原子數(shù)目分別為20，20)Selective dynamicsDirect 0.0000000000000000 0.0000000000000000 0.0000000000000000 F F F 0.5000000000000000 0.0000000000000000 0.0000000000000000 F F F 0.5000000000000000 0.5000000000000000 0.0000000000000000 F F F 0.0000000000000000 0.500000000000

24、0000 0.0000000000000000 F F F 0.2500000000000000 0.2500000000000000 0.0000000000000000 F F F 0.7500000000000000 0.2500000000000000 0.0000000000000000 F F F 0.2500000000000000 0.7500000000000000 0.0000000000000000 F F F 0.7500000000000000 0.7500000000000000 0.0000000000000000 F F F POTCAR文件內容說明： VASP

25、程序本身有提供了贗勢庫，只需將體系各類原子的贗勢合并在一起即可，但需注意到：1) 贗勢類型：US型贗勢LDAGGAPW91PBEPAW型贗勢GGAPW91PBELDAUS型贗勢所需截至能較小，計算速度快，PAW贗勢截至能通常較大，而且考慮的電子數(shù)多，計算慢，但精確度高。2) POTCAT中各原子贗勢定義的順序必需與POSCAR中相同：surface of mgo(100) (2*2)Mg 1.00000000000000 5.9459999999999997 0.0000000000000000 0.0000000000000000 0.0000000000000000 5.945999999

26、9999997 0.0000000000000000 0.0000000000000000 0.0000000000000000 20.0000000000000000 20 20Selective dynamicsDirect3) 對各原子的贗勢參數(shù)，我們最關心的是截至能以及電子數(shù);4) POTCAR的泛函類型必需與INCAR中GGA關鍵詞定義的 類型一致；5) 使用zcat命令產(chǎn)生和合并POTCAR文件。對應于中等大小的截至能(構型優(yōu)化時采用)對應于低的截至能(動力學模擬時采用)構造該贗勢時，所采用的泛函類型，這里為PW91電子數(shù)目和組態(tài)KPOINTS文件內容說明： 一般有兩種定義K點的方

27、法：1) 通過定義K-mesh大小，由程序自動產(chǎn)生各K點：Automatic mesh (title)0 (為0時，表示自動產(chǎn)生K點)M (表示采用Monkhorst-Pack方法生成K點坐標)5 5 5(對應于5x5x5網(wǎng)格)0 0 0(原點平移大小)2)手動定義各K點的坐標(一般僅在計算能帶結構時使用)： k-points for MgO(100) (title) 31 (K點數(shù)目)Rec (字母R打頭表示為倒易空間坐標，否則為實空間的坐標) 0.0 0.0 0.0 1.0 (各K點的坐標以及權重) 0.05 0.0 0.0 1.0 0.1 0.0 0.0 1.0 0.15 0.0 0.0

28、 1.0 0.2 0.0 0.0 1.0 0.25 0.0 0.0 1.0 0.3 0.0 0.0 1.0 0.35 0.0 0.0 1.0 0.4 0.0 0.0 1.0 0.45 0.0 0.0 1.0 0.5 0.0 0.0 1.06. VASP安裝和運行： (1) VASP程序安裝： a. 設置編譯環(huán)境：安裝Fortran編譯器，常用為IFC b. 對于并行版本vasp的編譯，還需安裝MPICH c. 編譯vasp自帶的庫文件 d. 對makefile進行修改，包括BLAS和Lapack庫文件所在 目錄，一般可采用IFC所帶的數(shù)學庫 e. 運行make命令進行編譯 (2) 創(chuàng)建輸入文件

29、，包括INCAR，KPOINTS，POSCAR 和POTCAR (3) 運行vasp： 單機版：/bin/vasp.4.5-ifc-mk-sp vasp.out &版本號編譯環(huán)境多個K點Single process并行版本：mpirun np 4 machinefile ./hosts /bin/vasp.4.5-ifc-mk-mp & vasp.out &CPU數(shù)目存放要并行運算的機器名或者IP常用關鍵詞使用說明(部分參考清華大學物理系蘇長榮編寫的VASP安裝和使用說明)(1)(2)一般單胞尺寸大時，選實空間，小單胞選取倒易空間。EDIFF=1e-4ENCUT=數(shù)值用戶手動定義截至能，如果沒

30、有，則由PREC選項確定。(3)(4)(5)(6)EDIFFG=EDIFF10當數(shù)值為負數(shù)時，表示以力作為收斂標準，多數(shù)情況均采用力作為收斂標準。ALGO=38|48該關鍵詞確定能量計算迭代方法38-采用Davidson優(yōu)化方法；(可靠，但速度慢)48-采用RMM-DIIS算法；(常用，速度快)ISYM=0|1|2該關鍵詞確定能量和構型優(yōu)化時是否使用對稱性(將影響到K點數(shù)目和計算量大小)0-不使用對稱性；1-采用對稱性；2-用于PAW型贗勢；(7)(8)NELM=整數(shù)該關鍵詞確定能量自洽場最大迭代輪數(shù)，缺省為60輪；NELMIN=整數(shù) 在構型優(yōu)化中，計算每個構象能量時最少迭代輪數(shù)，一般為34，

31、以保證能量和力的穩(wěn)定性；定義DFT泛函類型，注意要與POTCAR中的贗勢類型一致。(9)(10)ISPIN=1|21-非自旋極化計算(缺省)2-自旋極化計算, 將給出體系磁矩大小(對含有過渡金屬原子體系，一般均要采用自旋極化方法)。(11)(12)(13) ISMEAR選擇：1) 對半導體或絕緣體選取-5，如果單胞較大時，或者所選取k 點數(shù)目少時，用0；2) 對導體，通常用0； SIGMA取值： SIGMA取值的原則是使得計算得到的TS項(OUTCAR中)，分攤到每個原子上時小于1meV，否則得到的總能量不準確，對導體尤其要注意該參數(shù)的選擇。以下為構型優(yōu)化所用關鍵詞：NSW=整數(shù) 構型優(yōu)化的最

32、大輪數(shù)IBRION = -1|0|1|2 構型優(yōu)化方法：-1-構型不變更；0-分子動力學模擬；1-采用準牛頓方法確定新的構型(當初始構型較合理時使用)；2-采用CG方法確定構型(當初始構型離平衡位置較遠時使用)。POTIM=數(shù)值 控制構型優(yōu)化步長，缺省為0.5，對動力學模擬則為時間步長(單位為fs)(14)輸出控制關鍵詞： LCHARG = .FALSE. (輸出電荷密度？) LWAVE = .FALSE. (輸出波函數(shù)？) LVTOT = .FALSE.(輸出靜電勢，求功函時使用)其他關鍵詞： NPAR = 8 (CPU數(shù)目，并行計算時使用) LPLANE = .TRUE.(與并行算法有關)(15)(16) 實例： SYSTEM = Silicon bulk NPAR = 2 LPLANE = .TRUE. Elecronic minimisation ISTART = 0 LREAL= .FALSE. PREC = Medium precission: Mediun/High/Low EDIFF = 1e-4 converge criterion: default = 1e-4 EDIFFG = -0.02 converge criterion for relation loop IALGO =