基于氧化锌原胞的晶格优化方法研究*

牙莉荀，黄宝丹，劳妃玲，许钟华，房 慧，陈春燕

(广西民族师范学院物理与电子工程学院，广西崇左 532200)

0 引言

氧化锌(ZnO)是一种宽禁带的半导体材料，因其优良的光电性质而备受关注，有着较长的研究历史。许多研究组曾利用第一性原理对ZnO材料进行理论研究，比如关于ZnO基本电子结构性质、光学性质和磁学性质的研究[1-4]，关于元素掺杂ZnO的性质研究[5-8]，以及关于ZnO本征缺陷的研究[9-11]等,但目前还没有关于利用不同晶格优化方法优化ZnO原胞结构的研究报道。

准确的结构优化是获取可信材料性质的前提。VASP是一款基于第一性原理计算方法的材料性质计算程序包，在材料计算研究中有着广泛的应用[12-14]。VASP程序包中也提供了多种晶体优化的方法。一般来说，状态方程拟合法可以得到更准确的结果，但过程非常繁琐；自动优化法的结果会有所偏差，但过程方便简单。2种优化方法的对比研究未见有前人报道。

对于很多研究者，特别是初学者而言，常因如何选择合适的晶格优化方法而困扰。本文以ZnO原胞作为晶格优化对象，采用VASP计算程序包，分别使用状态方程拟合法和自动优化法对ZnO原胞进行晶格优化，对比研究2种优化方法所得到晶格体积和晶格参数的差异，并研究在不同的初始原胞体积下，自动优化法所得结果的差异，为一般材料的晶格体积优化方法的选择提供一定指导和参考。

1 计算方法

本文使用基于密度泛函理论(DFT)[15]的第一性原理方法来研究ZnO晶体体积的优化方法，其中使用VASP计算程序包进行晶格优化，使用VESTA[16]软件进行ZnO原胞模型构建以及晶格参数的读取。在具体优化过程中，采用超软赝势(USPP)来描述电子与离子的交换作用，使用PW91形式的交互关联泛函来处理关联作用，系统总能量的收敛判据为1×10-5eV/atom，力收敛判据为1×10-3eV/Å，布里渊区k点网格取样采用Gamma方法。k网格数目(KPOINTS)以及平面波截断能(ENCUT)在不同优化方法下进行具体的优化后选取。所有计算都进行2次计算：一次结构优化计算，一次静态计算。读取静态计算得到的能量作为体系能量。

KPOINTS取值的优化方法如下：设置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取值的关系图，然后根据关系图选取合适的KPOINTS值。

ENCUT取值的优化方法如下：设置KPOINTS为9×9×9，ENCUT取值分别为350，400，450，500，550，600 eV，计算并绘制体系能量和ENCUT取值的关系图，然后根据关系图选取合适的ENCUT值。

ZnO原胞模型的构建方法如下：使用VESTA软件构建纤锌矿ZnO的原胞模型，其空间群为P63mc (NO.186)，每个纤锌矿ZnO原胞包含4个原子(其中2个Zn原子，2个O原子)，Zn原子与O原子各自以密堆积的方式排列， Zn原子位于相邻4个O原子所形成的四面体中心， O原子位于相邻4个Zn原子所形成的四面体中心。ZnO原胞的晶格参数参考多篇实验研究数据[17-21]后选取，采用的晶格参数为a=b=3.249 3 Å，c=5.205 4 Å,α=β=90°,γ=120°。

状态方程拟合法优化ZnO原胞的步骤如下：设置参数ISIF=2，即在优化过程中保持晶格体积不变，只优化离子位置。对ENCUT和KPOINTS参数进行优化测试后选定ENCUT和KPOINTS的取值。在保持原始晶格c/a=1.602不变的情况下，对原胞体积进行缩放，缩放系数分别为0.80，0.85，0.90，0.95，1.00，1.05，1.10，1.15，1.20。经过缩放后，初始体积发生了变化，计算并绘制体系能量与体积之间的关系图(E-V图)。然后根据Brich-Murnaghan状态方程[22](BM状态方程)对E-V图进行拟合，得到平衡体积V0，其中BM状态方程如下：

(1)

其中V0，E0，B0,B′0分别为平衡体积、平衡态能量、平衡态弹性模量和弹性模量的一阶导数。最后将得到的平衡体积再次进行优化计算，得到晶格参数。

自动优化法优化ZnO原胞的步骤如下：设置参数ISIF=3，即在优化过程中同时优化晶格的体积和离子位置。对ENCUT和KPOINTS参数进行优化测试后，选定ENCUT和KPOINTS的取值，直接对ZnO原胞进行晶格优化，即可得到优化后的晶格体积和晶格参数。

测试不同初始体积下自动优化法优化结果的步骤如下：设置ISIF=3，ENCUT和KPOINTS参数选取测试后的取值，在保持原始晶格c/a=1.602不变的情况下，对初始体积进行缩放，缩放系数分别为0.80，0.85，0.90，0.95，1.00，1.05，1.10，1.15，1.20，得到不同初始体积的ZnO原胞。直接对不同初始体积的ZnO原胞进行晶格优化计算，得到优化的结果。

2 结果与分析

2.1 ZnO原胞模型的构建

图1是在VESTA软件中构建的纤锌矿ZnO原胞模型，其中灰色大球代表锌(Zn)原子，红色小球代表氧(O)原子。

图1 纤锌矿ZnO原胞模型

2.2 状态方程拟合法

2.2.1 KPOINTS和ENCUT参数优化

在设定ISIF=2的情况下，优化KPOINTS和ENCUT参数，以保证计算过程能够稳定地收敛。图2a为ISIF=2时，体系能量(Energy)随KPOINTS取值的变化，其中横坐标数字n代表KPOINTS取值为n×n×n(n从3变化到15)。从图2a可以知道，随着KPOINTS取值的增大，整个体系的能量逐渐降低。在KPOINTS取值为6×6×6以上，整个系统的能量基本趋于稳定，说明KPOINTS取值大于6×6×6，系统计算就可以达到很好的收敛效果。KPOINTS取值越大，计算结果越精确，但同时所需要的计算量也越大。综合考虑，在此部分计算中KPOINTS取值为9×9×9。

图2b为ISIF=2时，体系能量随ENCUT取值的变化。从图2b可知，随着ENCUT取值增大，体系能量逐渐降低，最后达到稳定平衡状态。当ENCUT取值由500 eV变到550 eV时，能量变化只有1.4×10-3eV，说明当ENCUT取值为500 eV时，系统计算已经达到很好的收敛效果。随着ENCUT取值变大，计算量也会随之增大。综合考虑， 在此部分的计算中，ENCUT的取值为500 eV。

图2 ISIF=2时,体系能量随KPOINTS和ENCUT取值的变化

2.2.2E-V图拟合

选取优化后的计算参数(KPOINTS取值为9×9×9，ENCUT取值为500 eV),在保持晶体c/a=1.602不变的情况下，对体积进行缩放，初始原胞的体积为47.595 3 Å3。缩放后对应的晶体体积为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可以明显看出，体系能量随着晶体体积的增大呈现先减少后增加的“碗状”变化:体积从24.368 8 Å3增大到47.595 3 Å3的过程中，整个体系能量逐渐降低，而体积从47.595 3 Å3增大到82.244 6 Å3的过程中，整个体系能量是逐渐增加的。根据BM状态方程的拟合，得到平衡体积V0=49.206 2 Å3，在平衡体积的条件下再次进行晶格优化,得到晶格常数a=b=3.285 6 Å，c=5.263 5 Å。状态方程拟合法优化结束后，ZnO原胞仍然保持纤锌矿的结构。

图3 体系能量随晶体体积的变化

实验过程中，合成方法环境温度以及形貌等因素都会对ZnO晶格参数产生影响，不同实验组测量得到的ZnO晶格常数也会有所区别[17-21]。部分实验组测量的实验数据见表1。由于ZnO原胞建模时所采用的晶格参数是参考多篇实验数值获得，故采用建模时所使用的晶格参数a=b=3.249 3 Å，c=5.205 4 Å，体积V=47.595 3 Å3作为ZnO实验参考值来对计算优化的数值进行分析比较。相比于实验参考值，状态方程拟合法得到的晶格常数a，b，c均偏大，误差都为1.12%，晶体体积的误差为3.38%。综上可见，晶格常数的误差在2%内，晶体体积的误差在4%内，属于合理的误差范围，证明计算结果的可靠性。

表1 不同实验组得到的ZnO晶格常数

2.3 自动优化法

2.3.1 计算参数的优化

图4a为ISIF=3时,体系能量随KPOINTS取值的变化，其中横坐标数字n代表KPOINTS取值为n×n×n(n从3变化到15)。从图4a可以知道，随着KPOINTS取值增大，整个体系能量逐渐降低。在KPOINTS取值为6×6×6以上，整个系统的能量基本趋于稳定，说明KPOINTS取值大于6×6×6，系统计算就可以达到很好的收敛效果。KPOINTS取值越大，计算结果越精确，同时所需要的计算量也越大。综合考虑，在此部分计算中KPOINTS取值为9×9×9。

图4b为ISIF=3时，体系能量随ENCUT取值的变化。从图4b可知，随着ENCUT取值增大，体系能量逐渐降低，最后达到稳定平衡状态。当ENCUT取值由500 eV变到550 eV时，能量变化只有1.4×10-3eV，说明当ENCUT取值为500 eV时，系统的能量基本达到最小值并趋于稳定，说明此时系统计算可以达到很好的收敛效果。随着ENCUT的取值变大，计算量也会随之增大。综合考虑，在此部分的计算中，ENCUT的取值为500 eV。

图4 ISIF=3时,体系能量随KPOINTS和ENCUT取值的变化

2.3.2 自动优化的结果

在VASP的计算中，设置ISIF=3，即在晶格优化过程中，计算系统会同时优化离子位置和晶格体积。自动优化结束后，ZnO原胞仍保持纤锌矿结构，晶格常数a=b=3.282 0 Å，c=5.260 1 Å，晶格体积为49.068 4 Å3。与实验参考值相比，晶格常数a，b的误差为1.01%，晶格常数c的误差为1.05%，晶体体积的误差为3.1%，均属于合理的误差范围，证明计算结果的可靠性。与状态方程拟合法得到的参数相比，自动优化法晶格常数a，b的误差为0.11%，晶格常数c的误差为0.06%，晶体体积的误差为0.28%。2种优化方法的结果相近，所以对于一般晶格的优化，直接采用方便快捷的自动优化法即可。

2.3.3 不同初始体积下自动优化的结果

为进一步研究不同的初始体积下使用自动优化法对最后计算结果是否有影响，本文首先通过不同的缩放系数得到不同的晶格初始体积，再通过系统自动优化得到晶格参数。初始的晶体体积分别为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为自动优化得到的晶体体积随不同初始体积的变化情况。由图5可以看出，整体而言，随着初始体积的增加，优化后体积一直在波动，无明显的函数关系。

图5 自动优化后的晶格体积随初始体积的变化

不同初始体积条件下自动优化得到的晶格参数见表2。对于不同的晶格初始体积，优化后的体积为48.826 0-49.108 4 Å3；晶格常数a，b为3.267 4-3.286 9 Å；晶格常数c为5.239 8-5.296 8 Å。

表2 不同初始体积下自动优化法得到的晶格参数

不同初始体积下自动优化法得到的晶格参数与实验参考值的误差见表3。从表3可以看出，对比实验参考值，不同初始体积自动优化得到的晶格体积的误差为2.59%- 3.18%，晶格常数a，b的误差为0.56%-1.16%，晶格常数c的误差为0.66%-1.76%，都在合理的误差范围内，证明在不同的初始体积下，自动优化法仍能给出合理的优化结果。

表3 不同初始体积下自动优化法得到的晶格参数与实验参考值的误差

不同初始体积下自动优化法得到的晶格参数与状态方程拟合值的误差见表4。从表4可以看出，对比状态方程拟合值，不同初始体积自动优化得到的晶格体积误差为0.20%-0.77%，晶格常数a，b的误差为0.04%-0.55%，晶格常数c的误差为0.02%-0.63%。整体误差不超过1%，证明在不同晶格初始体积的条件下，自动优化法仍然能给出与状态方程拟合法相似的结果。由此进一步说明，对于类似ZnO原胞晶格体系的优化，采用方便简单的自动优化法取代状态方程拟合法是可行的。

表4 不同初始体积下自动优化法得到的晶格参数与状态方程拟合值的误差

3 结论

与实验参考值相比，2种优化方法得到的晶格常数误差不超过2%，晶体体积误差不超过4%，属于合理的误差范围，2种优化方法都能得到比较合理的晶体参数结果。使用状态方程拟合法和自动优化法优化得到的晶格体积之间的误差为0.28%，晶格常数a，b的误差为0.11%，晶格常数c的误差为0.06%。在不同的初始体积下，自动优化法的结果随着初始体积的增加无明显变化关系。与状态方程拟合值相比，在不同的初始体积下使用自动优化法得到的晶格体积误差为0.20%-0.77%，晶格常数a，b的误差为0.04%-0.55%，晶格常数c的误差为0.02%-0.63%，误差都不超过1%。可见，对于类似ZnO原胞简单晶格体系的优化，完全可以使用系统自动优化法代替状态方程拟合法。本研究结果为晶格体积优化方法的选择提供了一定的理论参考价值。