牙莉荀，黃寶丹，勞妃玲，許鐘華，房 慧，陳春燕

(廣西民族師范學院物理與電子工程學院，廣西崇左 532200)

0 引言

氧化鋅(ZnO)是一種寬禁帶的半導體材料，因其優(yōu)良的光電性質(zhì)而備受關(guān)注，有著較長的研究歷史。許多研究組曾利用第一性原理對ZnO材料進行理論研究，比如關(guān)于ZnO基本電子結(jié)構(gòu)性質(zhì)、光學性質(zhì)和磁學性質(zhì)的研究[1-4]，關(guān)于元素摻雜ZnO的性質(zhì)研究[5-8]，以及關(guān)于ZnO本征缺陷的研究[9-11]等,但目前還沒有關(guān)于利用不同晶格優(yōu)化方法優(yōu)化ZnO原胞結(jié)構(gòu)的研究報道。

準確的結(jié)構(gòu)優(yōu)化是獲取可信材料性質(zhì)的前提。VASP是一款基于第一性原理計算方法的材料性質(zhì)計算程序包，在材料計算研究中有著廣泛的應用[12-14]。VASP程序包中也提供了多種晶體優(yōu)化的方法。一般來說，狀態(tài)方程擬合法可以得到更準確的結(jié)果，但過程非常繁瑣；自動優(yōu)化法的結(jié)果會有所偏差，但過程方便簡單。2種優(yōu)化方法的對比研究未見有前人報道。

對于很多研究者，特別是初學者而言，常因如何選擇合適的晶格優(yōu)化方法而困擾。本文以ZnO原胞作為晶格優(yōu)化對象，采用VASP計算程序包，分別使用狀態(tài)方程擬合法和自動優(yōu)化法對ZnO原胞進行晶格優(yōu)化，對比研究2種優(yōu)化方法所得到晶格體積和晶格參數(shù)的差異，并研究在不同的初始原胞體積下，自動優(yōu)化法所得結(jié)果的差異，為一般材料的晶格體積優(yōu)化方法的選擇提供一定指導和參考。

1 計算方法

本文使用基于密度泛函理論(DFT)[15]的第一性原理方法來研究ZnO晶體體積的優(yōu)化方法，其中使用VASP計算程序包進行晶格優(yōu)化，使用VESTA[16]軟件進行ZnO原胞模型構(gòu)建以及晶格參數(shù)的讀取。在具體優(yōu)化過程中，采用超軟贗勢(USPP)來描述電子與離子的交換作用，使用PW91形式的交互關(guān)聯(lián)泛函來處理關(guān)聯(lián)作用，系統(tǒng)總能量的收斂判據(jù)為1×10-5eV/atom，力收斂判據(jù)為1×10-3eV/?，布里淵區(qū)k點網(wǎng)格取樣采用Gamma方法。k網(wǎng)格數(shù)目(KPOINTS)以及平面波截斷能(ENCUT)在不同優(yōu)化方法下進行具體的優(yōu)化后選取。所有計算都進行2次計算：一次結(jié)構(gòu)優(yōu)化計算，一次靜態(tài)計算。讀取靜態(tài)計算得到的能量作為體系能量。

KPOINTS取值的優(yōu)化方法如下：設(shè)置ENCUT為500 eV，KPOINTS分別取點3×3×3，4×4×4，5×5×5，6×6×6，7×7×7，8×8×8，9×9×9，10×10×10，11×11×11，12×12×12，13×13×13，14×14×14，15×15×15，計算并繪制體系能量和KPOINTS取值的關(guān)系圖，然后根據(jù)關(guān)系圖選取合適的KPOINTS值。

ENCUT取值的優(yōu)化方法如下：設(shè)置KPOINTS為9×9×9，ENCUT取值分別為350，400，450，500，550，600 eV，計算并繪制體系能量和ENCUT取值的關(guān)系圖，然后根據(jù)關(guān)系圖選取合適的ENCUT值。

ZnO原胞模型的構(gòu)建方法如下：使用VESTA軟件構(gòu)建纖鋅礦ZnO的原胞模型，其空間群為P63mc (NO.186)，每個纖鋅礦ZnO原胞包含4個原子(其中2個Zn原子，2個O原子)，Zn原子與O原子各自以密堆積的方式排列， Zn原子位于相鄰4個O原子所形成的四面體中心， O原子位于相鄰4個Zn原子所形成的四面體中心。ZnO原胞的晶格參數(shù)參考多篇實驗研究數(shù)據(jù)[17-21]后選取，采用的晶格參數(shù)為a=b=3.249 3 ?，c=5.205 4 ?,α=β=90°,γ=120°。

狀態(tài)方程擬合法優(yōu)化ZnO原胞的步驟如下：設(shè)置參數(shù)ISIF=2，即在優(yōu)化過程中保持晶格體積不變，只優(yōu)化離子位置。對ENCUT和KPOINTS參數(shù)進行優(yōu)化測試后選定ENCUT和KPOINTS的取值。在保持原始晶格c/a=1.602不變的情況下，對原胞體積進行縮放，縮放系數(shù)分別為0.80，0.85，0.90，0.95，1.00，1.05，1.10，1.15，1.20。經(jīng)過縮放后，初始體積發(fā)生了變化，計算并繪制體系能量與體積之間的關(guān)系圖(E-V圖)。然后根據(jù)Brich-Murnaghan狀態(tài)方程[22](BM狀態(tài)方程)對E-V圖進行擬合，得到平衡體積V0，其中BM狀態(tài)方程如下：

(1)

其中V0，E0，B0,B′0分別為平衡體積、平衡態(tài)能量、平衡態(tài)彈性模量和彈性模量的一階導數(shù)。最后將得到的平衡體積再次進行優(yōu)化計算，得到晶格參數(shù)。

自動優(yōu)化法優(yōu)化ZnO原胞的步驟如下：設(shè)置參數(shù)ISIF=3，即在優(yōu)化過程中同時優(yōu)化晶格的體積和離子位置。對ENCUT和KPOINTS參數(shù)進行優(yōu)化測試后，選定ENCUT和KPOINTS的取值，直接對ZnO原胞進行晶格優(yōu)化，即可得到優(yōu)化后的晶格體積和晶格參數(shù)。

測試不同初始體積下自動優(yōu)化法優(yōu)化結(jié)果的步驟如下：設(shè)置ISIF=3，ENCUT和KPOINTS參數(shù)選取測試后的取值，在保持原始晶格c/a=1.602不變的情況下，對初始體積進行縮放，縮放系數(shù)分別為0.80，0.85，0.90，0.95，1.00，1.05，1.10，1.15，1.20，得到不同初始體積的ZnO原胞。直接對不同初始體積的ZnO原胞進行晶格優(yōu)化計算，得到優(yōu)化的結(jié)果。

2 結(jié)果與分析

2.1 ZnO原胞模型的構(gòu)建

圖1是在VESTA軟件中構(gòu)建的纖鋅礦ZnO原胞模型，其中灰色大球代表鋅(Zn)原子，紅色小球代表氧(O)原子。

圖1 纖鋅礦ZnO原胞模型

2.2 狀態(tài)方程擬合法

2.2.1 KPOINTS和ENCUT參數(shù)優(yōu)化

在設(shè)定ISIF=2的情況下，優(yōu)化KPOINTS和ENCUT參數(shù)，以保證計算過程能夠穩(wěn)定地收斂。圖2a為ISIF=2時，體系能量(Energy)隨KPOINTS取值的變化，其中橫坐標數(shù)字n代表KPOINTS取值為n×n×n(n從3變化到15)。從圖2a可以知道，隨著KPOINTS取值的增大，整個體系的能量逐漸降低。在KPOINTS取值為6×6×6以上，整個系統(tǒng)的能量基本趨于穩(wěn)定，說明KPOINTS取值大于6×6×6，系統(tǒng)計算就可以達到很好的收斂效果。KPOINTS取值越大，計算結(jié)果越精確，但同時所需要的計算量也越大。綜合考慮，在此部分計算中KPOINTS取值為9×9×9。

圖2b為ISIF=2時，體系能量隨ENCUT取值的變化。從圖2b可知，隨著ENCUT取值增大，體系能量逐漸降低，最后達到穩(wěn)定平衡狀態(tài)。當ENCUT取值由500 eV變到550 eV時，能量變化只有1.4×10-3eV，說明當ENCUT取值為500 eV時，系統(tǒng)計算已經(jīng)達到很好的收斂效果。隨著ENCUT取值變大，計算量也會隨之增大。綜合考慮， 在此部分的計算中，ENCUT的取值為500 eV。

圖2 ISIF=2時,體系能量隨KPOINTS和ENCUT取值的變化

2.2.2E-V圖擬合

選取優(yōu)化后的計算參數(shù)(KPOINTS取值為9×9×9，ENCUT取值為500 eV),在保持晶體c/a=1.602不變的情況下，對體積進行縮放，初始原胞的體積為47.595 3 ?3?？s放后對應的晶體體積為24.368 8，29.229 5，34.697 0，40.807 0，47.595 3，55.097 5，63.349 3，72.386 5，82.244 6 ?3。

圖3為體系能量隨晶體體積的變化情況，即E-V圖。從圖3可以明顯看出，體系能量隨著晶體體積的增大呈現(xiàn)先減少后增加的“碗狀”變化:體積從24.368 8 ?3增大到47.595 3 ?3的過程中，整個體系能量逐漸降低，而體積從47.595 3 ?3增大到82.244 6 ?3的過程中，整個體系能量是逐漸增加的。根據(jù)BM狀態(tài)方程的擬合，得到平衡體積V0=49.206 2 ?3，在平衡體積的條件下再次進行晶格優(yōu)化,得到晶格常數(shù)a=b=3.285 6 ?，c=5.263 5 ?。狀態(tài)方程擬合法優(yōu)化結(jié)束后，ZnO原胞仍然保持纖鋅礦的結(jié)構(gòu)。

圖3 體系能量隨晶體體積的變化

實驗過程中，合成方法環(huán)境溫度以及形貌等因素都會對ZnO晶格參數(shù)產(chǎn)生影響，不同實驗組測量得到的ZnO晶格常數(shù)也會有所區(qū)別[17-21]。部分實驗組測量的實驗數(shù)據(jù)見表1。由于ZnO原胞建模時所采用的晶格參數(shù)是參考多篇實驗數(shù)值獲得，故采用建模時所使用的晶格參數(shù)a=b=3.249 3 ?，c=5.205 4 ?，體積V=47.595 3 ?3作為ZnO實驗參考值來對計算優(yōu)化的數(shù)值進行分析比較。相比于實驗參考值，狀態(tài)方程擬合法得到的晶格常數(shù)a，b，c均偏大，誤差都為1.12%，晶體體積的誤差為3.38%。綜上可見，晶格常數(shù)的誤差在2%內(nèi)，晶體體積的誤差在4%內(nèi)，屬于合理的誤差范圍，證明計算結(jié)果的可靠性。

表1 不同實驗組得到的ZnO晶格常數(shù)

2.3 自動優(yōu)化法

2.3.1 計算參數(shù)的優(yōu)化

圖4a為ISIF=3時,體系能量隨KPOINTS取值的變化，其中橫坐標數(shù)字n代表KPOINTS取值為n×n×n(n從3變化到15)。從圖4a可以知道，隨著KPOINTS取值增大，整個體系能量逐漸降低。在KPOINTS取值為6×6×6以上，整個系統(tǒng)的能量基本趨于穩(wěn)定，說明KPOINTS取值大于6×6×6，系統(tǒng)計算就可以達到很好的收斂效果。KPOINTS取值越大，計算結(jié)果越精確，同時所需要的計算量也越大。綜合考慮，在此部分計算中KPOINTS取值為9×9×9。

圖4b為ISIF=3時，體系能量隨ENCUT取值的變化。從圖4b可知，隨著ENCUT取值增大，體系能量逐漸降低，最后達到穩(wěn)定平衡狀態(tài)。當ENCUT取值由500 eV變到550 eV時，能量變化只有1.4×10-3eV，說明當ENCUT取值為500 eV時，系統(tǒng)的能量基本達到最小值并趨于穩(wěn)定，說明此時系統(tǒng)計算可以達到很好的收斂效果。隨著ENCUT的取值變大，計算量也會隨之增大。綜合考慮，在此部分的計算中，ENCUT的取值為500 eV。

圖4 ISIF=3時,體系能量隨KPOINTS和ENCUT取值的變化

2.3.2 自動優(yōu)化的結(jié)果

在VASP的計算中，設(shè)置ISIF=3，即在晶格優(yōu)化過程中，計算系統(tǒng)會同時優(yōu)化離子位置和晶格體積。自動優(yōu)化結(jié)束后，ZnO原胞仍保持纖鋅礦結(jié)構(gòu)，晶格常數(shù)a=b=3.282 0 ?，c=5.260 1 ?，晶格體積為49.068 4 ?3。與實驗參考值相比，晶格常數(shù)a，b的誤差為1.01%，晶格常數(shù)c的誤差為1.05%，晶體體積的誤差為3.1%，均屬于合理的誤差范圍，證明計算結(jié)果的可靠性。與狀態(tài)方程擬合法得到的參數(shù)相比，自動優(yōu)化法晶格常數(shù)a，b的誤差為0.11%，晶格常數(shù)c的誤差為0.06%，晶體體積的誤差為0.28%。2種優(yōu)化方法的結(jié)果相近，所以對于一般晶格的優(yōu)化，直接采用方便快捷的自動優(yōu)化法即可。

2.3.3 不同初始體積下自動優(yōu)化的結(jié)果

為進一步研究不同的初始體積下使用自動優(yōu)化法對最后計算結(jié)果是否有影響，本文首先通過不同的縮放系數(shù)得到不同的晶格初始體積，再通過系統(tǒng)自動優(yōu)化得到晶格參數(shù)。初始的晶體體積分別為24.368 8，29.229 5，34.697 0，40.807 0，47.595 3，55.097 5，63.349 3，72.386 5，82.244 6 ?3。圖5為自動優(yōu)化得到的晶體體積隨不同初始體積的變化情況。由圖5可以看出，整體而言，隨著初始體積的增加，優(yōu)化后體積一直在波動，無明顯的函數(shù)關(guān)系。

圖5 自動優(yōu)化后的晶格體積隨初始體積的變化

不同初始體積條件下自動優(yōu)化得到的晶格參數(shù)見表2。對于不同的晶格初始體積，優(yōu)化后的體積為48.826 0-49.108 4 ?3；晶格常數(shù)a，b為3.267 4-3.286 9 ?；晶格常數(shù)c為5.239 8-5.296 8 ?。

表2 不同初始體積下自動優(yōu)化法得到的晶格參數(shù)

不同初始體積下自動優(yōu)化法得到的晶格參數(shù)與實驗參考值的誤差見表3。從表3可以看出，對比實驗參考值，不同初始體積自動優(yōu)化得到的晶格體積的誤差為2.59%- 3.18%，晶格常數(shù)a，b的誤差為0.56%-1.16%，晶格常數(shù)c的誤差為0.66%-1.76%，都在合理的誤差范圍內(nèi)，證明在不同的初始體積下，自動優(yōu)化法仍能給出合理的優(yōu)化結(jié)果。

表3 不同初始體積下自動優(yōu)化法得到的晶格參數(shù)與實驗參考值的誤差

不同初始體積下自動優(yōu)化法得到的晶格參數(shù)與狀態(tài)方程擬合值的誤差見表4。從表4可以看出，對比狀態(tài)方程擬合值，不同初始體積自動優(yōu)化得到的晶格體積誤差為0.20%-0.77%，晶格常數(shù)a，b的誤差為0.04%-0.55%，晶格常數(shù)c的誤差為0.02%-0.63%。整體誤差不超過1%，證明在不同晶格初始體積的條件下，自動優(yōu)化法仍然能給出與狀態(tài)方程擬合法相似的結(jié)果。由此進一步說明，對于類似ZnO原胞晶格體系的優(yōu)化，采用方便簡單的自動優(yōu)化法取代狀態(tài)方程擬合法是可行的。

表4 不同初始體積下自動優(yōu)化法得到的晶格參數(shù)與狀態(tài)方程擬合值的誤差

3 結(jié)論

與實驗參考值相比，2種優(yōu)化方法得到的晶格常數(shù)誤差不超過2%，晶體體積誤差不超過4%，屬于合理的誤差范圍，2種優(yōu)化方法都能得到比較合理的晶體參數(shù)結(jié)果。使用狀態(tài)方程擬合法和自動優(yōu)化法優(yōu)化得到的晶格體積之間的誤差為0.28%，晶格常數(shù)a，b的誤差為0.11%，晶格常數(shù)c的誤差為0.06%。在不同的初始體積下，自動優(yōu)化法的結(jié)果隨著初始體積的增加無明顯變化關(guān)系。與狀態(tài)方程擬合值相比，在不同的初始體積下使用自動優(yōu)化法得到的晶格體積誤差為0.20%-0.77%，晶格常數(shù)a，b的誤差為0.04%-0.55%，晶格常數(shù)c的誤差為0.02%-0.63%，誤差都不超過1%?？梢?，對于類似ZnO原胞簡單晶格體系的優(yōu)化，完全可以使用系統(tǒng)自動優(yōu)化法代替狀態(tài)方程擬合法。本研究結(jié)果為晶格體積優(yōu)化方法的選擇提供了一定的理論參考價值。