氮化镓相图预测及其高压熔化特性研究*

雷振帅 孙小伟 刘子江 宋婷 田俊红

(兰州交通大学数理学院,兰州 730070)

采用经典分子动力学模拟,结合第一性原理计算及晶格动力学方法对氮化镓(GaN)纤锌矿结构与岩盐结构在 0—80 GPa 压力范围内的相图进行了预测.第一性原理计算与分子动力学模拟得到的零温下 GaN 纤锌矿到岩盐结构的相变压力分别为44.3 GPa 和 45.9 GPa,与已有研究的实验结果吻合;通过外推纤锌矿结构 GaN的熔化曲线得到其零压下的熔化温度为2295 K,当压力增大到33.3 GPa 时,纤锌矿结构熔化曲线与岩盐结构熔化曲线相交,两种结构的熔化温度均随压力的增大而上升;GaN 还可能存在超离子相,纤锌矿结构在压力大于2.0 GPa 且温度高于 2550 K 时发生超离子相转变,岩盐结构在压力温度大于 33.1 GPa 和4182 K 后发生超离子相转变,二者的相转变温度均会随着压力的增大而升高;GaN 纤锌矿和岩盐结构的相界线并非为一直线,在高温下相界线斜率为正,随着温度的降低逐渐变为一条具有负斜率的曲线.

1 引言

Ⅲ-Ⅴ族化合物半导体氮化镓(GaN)具有禁带宽度大、热导率高、电子饱和速度快等优异特性,是发展高频、高功率电子器件的优良半导体材料,在发光二极管、晶体管与军事电子设备中均具有广泛的应用前景[1].然而,相比于微电子领域,人们从基础性质的角度给予 GaN的研究还不够充分,如高温下固-固相变的克劳修斯-克拉珀龙斜率与熔化温度仍存在很大争议[2-6],究其原因,主要是因为GaN的分解温度低于熔化温度,在熔化时具有非常高的氮气平衡蒸气压[7],在进行GaN的高温实验时,需要保持高的氮气压力以防止其分解,这使得实验研究变得尤为困难,在一定程度上限制了GaN的开发和应用.因此,GaN的固-固相变和熔化温度等基础性质成为了该材料发展过程中的重要研究内容.

理论和实验研究表明,GaN 具有3 种典型结构,即纤锌矿结构、闪锌矿结构和岩盐结构[3,4,8].在环境条件下GaN 一般以纤锌矿结构稳定存在,闪锌矿结构可通过气相外延和分子束外延沉积在立方(001)衬底上,以薄膜形式稳定存在于环境条件中[9],岩盐结构为稳定的高压相.GaN 纤锌矿到岩盐结构的相变发生在 37—52 GPa 之间[8,10],而闪锌矿到岩盐结构的相变发生在 36—42 GPa 之间[11].GaN 纤锌矿和闪锌矿结构均为Ga-N 正四面体,二者的结构相近,实验上还未观察到该相变,但也有研究表明该相变发生在 10.87 GPa[11].GaN的高温高压实验条件苛刻,而可以计算高温相变的分子模拟又缺乏准确的原子力场,因此,对 GaN 进行高温高压相图的研究是一项困难的工作.Van Vechten[2]提出的化学键模型预测了 GaN 在高温下的固-固相变压力,但该模型预测的是纤锌矿结构与β-tin 结构的相变.关于 GaN 纤锌矿到岩盐结构的相界线研究较少,仅有 Zhou 等[3]的数值模拟和 Sadovyi 等[4]的高温高压实验结果,后者对GaN的固-固相变行为进行了较为全面的分析.

GaN 熔化温度的研究结果存在较大分歧,Van Vechten[2]提出的模型预测了常压下GaN的熔化温度为2791 K,与之后的实验结果比较接近,但该模型预测的熔化曲线呈一次函数下降趋势,这与之后的许多研究结果不符;Karpiński 等[7]在高压碳化钨砧槽中进行了高温实验,发现压力为6.0 GPa时 GaN的熔点高于 2573 K;Utsumi 等[5]的X 射线衍射实验给出了压力大于6.0 GPa 时,GaN的熔化温度为2493 K,随后 Sokol 等[12]与Saitoh 等[13]的研究也支持了Utsumi 等[5]的结果;Porowski 等[14]为提高高温高压下对样品特征的检测和分析质量,使用了比早期实验材料体积大了100 倍的GaN 样品,该研究发现在 6—9 GPa的压力范围内,温度达到 3400 K 时也只是观察到 GaN的分解而不是熔化,这与 Harafuji 等[6]与Nord 等[15]的分子动力学模拟结果相近.值得注意的是,Harafuji 等[6]的研究中出现了GaN 熔化温度的不连续现象,作者将其解释为带部分电荷的Ga 原子的固体电解质行为,在本工作中也发现了Ga 原子的类似行为(超离子相),但并没有出现熔化温度的不连续现象.

本文采用Born-Mayer 与Morse 形式结合的相互作用势进行经典分子动力学模拟,并通过基于密度泛函理论的第一性原理计算方法与晶格动力学方法,对GaN 晶体纤锌矿与岩盐结构的弹性常数、体积模量与晶格能进行了计算.在此基础上对GaN 在高压下的体积压缩行为进行了研究,同时采用两相法和单相法预测了不同结构GaN的熔化温度与超离子相,最后给出了GaN 在宽广的温度和压力范围内的P-T相图.

2 模型及方法

2.1 GaN 势函数及其验证

Alder 和Wainwright[16]于 1957 年提出分子动力学方法,该方法是一种原子层次的计算机模拟方法,由于不涉及电子而降低了计算量,可用于模拟庞大的原子体系.分子动力学模拟的核心是准确描述原子间相互作用的势函数,现有的势函数种类可以分为描述离子键的二体势、描述共价键的三体势以及描述金属键的多体势.GaN 已有的势函数非常丰富,如三体 Tersoff 势[15]与 Stillinger-Weber势[17,18]、多体修正嵌入原子势[19]以及对势类型的壳层模型Buckingham 势[20]与Morse 势[21]等,这些势函数主要用于 GaN的缺陷、力学性质等温度较低的计算,在高温高压条件下可能不适用,因此需要对势函数进行筛选,本文最终确定使用2005 年Zhang 等[21]通过晶格反演方法得到的对势.GaN具有混合共价-离子键[15],使用二体势或三体势均可进行描述.图1为利用基于密度泛函理论的第一性原理计算方法得到的GaN 纤锌矿和岩盐结构的电子局域函数图,由图1 可以看出两种结构的电子几乎都局域在 N 原子周围,离子键更加明显,因此对势形式的势函数可能会更好地描述GaN 各原子对间的相互作用:

图1 GaN 电子局域函数图 (a) 纤锌矿结构(110)晶面;(b) 岩盐结构(100)晶面Fig.1.Electron localization function of GaN (a) wurtzite structure (110) and (b) rocksalt structure (100) crystal planes.

(1)式中第1 项为库仑项,Zi,Zj分别为i原子与j原子的有效电荷,ε0为真空介电常数,r为原子间距,第2 项为i原子与j原子的短程排斥项.(2)式为N-N,Ga-Ga 间短程排斥相互作用,使用排斥指数形式的势函数描述,其中D,γ,R分别为势函数参数;(3)式为Ga-N 间短程排斥相互作用,使用 Morse 形式的势函数描述.排斥指数形式的势函数在分子动力学模拟中使用较少,而 Born-Mayer 势函数[22]则使用的更加广泛,这是因为其形式简单且参数更少,在计算材料热力学性质时往往也能给出可靠的计算结果,因此,本文将 Zhang等[21]的排斥指数形式的势函数拟合到了 Born-Mayer 形式:

式中A,η为势函数参数.使用该形式描述 N-N,Ga-Ga 间的短程相互作用,Ga-N 间的相互作用仍使用 Morse 形式描述,最终的势函数参数见表1.

表1 GaN的Morse 势参数与本文拟合得到的Born-Mayer 势参数Table 1.The Morse potential parameters and fitted Born-Mayer potential parameters of GaN.

为验证势函数的准确性,本文分别使用第一性原理计算方法与晶格动力学方法对 GaN 不同结构的弹性性质与晶格能进行了计算.第一性原理计算中,为使体系能量在完备平面波基矢水平上足够收敛,采用了 BFGS 算法[23]对晶体结构进行几何优化,以获得局域能量最低结构,选取由 Perdew 等修订的PBE 形式的广义梯度近似方法[24]描述电子间的交换关联能,选取非局域超软赝势[25]描述离子实和价电子间的相互作用,计算过程中 Ga 和N 原子的外层电子组态分别为3d104s24p1和 2s22p3,平面波截断能均为650 eV,GaN 纤锌矿与岩盐结构的倒易空间布里渊区k点分别取值13×13×7和9×9×9,在结构优化过程中体系能量收敛标准为5×10—6eV/atom,优化后作用在晶胞中每个原子上的力小于 0.01 eV/Å,晶胞应力偏差低于0.02 GPa,公差偏移小于 5×10—4Å.计算弹性常数时,最大应变振幅设置为0.003,每个应变循环6 步,即产生 6 个扭曲结构,所有的第一性原理计算均利用 CASTEP 软件包[26]完成.晶格动力学计算使用表1 中的势函数参数,截断半径选择 10 Å,所有的晶格动力学计算均利用 GULP 软件包[27]完成.得到的GaN 纤锌矿与岩盐结构的相关参数见表2,通过晶格动力学方法计算的纤锌矿结构的晶格常数、弹性常数与体积模量接近第一性原理计算和实验结果,这表明该势函数能够准确描述GaN 纤锌矿结构的力学性质,但岩盐结构的弹性常数与第一性原理计算数据存在差距,晶格动力学计算的两种结构的晶格能与已有数据非常接近,这表明该势函数能够很好地描述 GaN 在高温下尤其是熔化后的性质.为验证势函数在高压下的有效性,使用基于该势函数的分子动力学方法计算了GaN 纤锌矿与岩盐结构在 300 K 时的体积比率随压力的变化情况,如图2 所示,GaN 纤锌矿结构的体积比率随压力的变化与已有的实验与计算结果均吻合,在 45.9 GPa 相变到岩盐结构时,体积变化率为14.4%,比Xia 等[8]的实验结果(17.9%)低,但接近于Pandey 等[28]的从头算结果,说明该势函数能够描述 GaN 在高压条件下的热力学状态.

图2 GaN 纤锌矿和岩盐结构在 300 K 下的体积比率与已有结果对比,所有数据均归一化至纤锌矿结构初始体积,红色与黑色圆点为分子动力学模拟结果(实线为拟合曲线);蓝色正三角形与绿色菱形分别为Xia 等[8]与 Ueno 等[10]的X 射线衍射实验结果;洋红色虚线为Pandey 等[28]的从头算结果Fig.2.The volume ratios of GaN with wurtzite and rocksalt structures at 300 K are compared with existing results,all data are normalized to the initial volume of the wurtzite GaN,the red and black dots are the molecular dynamics simulations results (solid line is the fitted curve).The green diamond and blue square triangle are the X-ray diffraction experimental results by Ueno et al.[10] and Xia et al.[8],respectively.The magenta dashed line is the ab initio result by Pandey et al.[28].

表2 GaN 纤锌矿与岩盐结构在零压下的晶格常数,弹性常数 Cij,体积模量 B 及晶格能 ETable 2.The lattice constants,elastic constants Cij,bulk modulus B and lattice energy E of GaN with wurtzite and rocksalt structures at zero pressure.

2.2 GaN的熔化温度

为保证原子数目相同,本文分别使用8×8×16与6×8×16的正交体系进行 GaN 纤锌矿与岩盐结构的分子动力学模拟.两个体系均包含6144个原子,采用周期性边界条件,时间步长1 fs,所有的分子动力学计算均利用 LAMMPS 软件[32]完成.目前分子动力学模拟中计算熔化温度的方法主要有单相法、两相法、孔洞法、Z 方法以及改进的Z 方法[33-36],本文选择两相法进行 GaN的熔化温度计算,该方法准确度很高,而且也能在一定程度上控制压力,在进行熔化曲线的计算时,具有较好的操作空间.两相法计算熔化温度的关键在于构建能够长时间存在的两相共存体系,为达到这个目的,首先将整个模拟体系在等温等压系综(NPT)中以略低于熔点的温度弛豫30 ps,保证整个体系的应力平衡,然后将体系沿z方向平均分为两部分,一部分设置受力为0,将其固定,另一部分在远高于熔点的温度下弛豫 50 ps,保证其完全熔化,再将完全熔化的这部分原子降温到略高于熔点的温度,降温时间为10 ps,最后释放固定的原子,将整个体系置于微正则系综(NVE)中自由弛豫,形成固液共存体系,若整个体系能够以固液共存状态长时间存在(本文中为200 ps),且固液部分体积无显著变化,则可以认为此时体系的温度为最终熔化温度.图3为构建的GaN 两相共存体系,图4为对应的原子数密度,可以清楚地看出体系明显分为两部分,数密度波动大的部分处于固态,数密度波动较小的部分处于液态.

图3 处于两相共存状态的GaN 体系Fig.3.GaN system in the two-phase coexistence state.

图4 两相共存状态下 GaN 体系沿 z 轴方向的原子数密度Fig.4.Atomic number density along the z-direction of the GaN system in the two-phase coexistence state.

2.3 GaN的超离子相

超离子相的确定不涉及过热问题,因此本文使用单相法对 GaN的超离子相进行确定.为准确控制压力和温度,选择使用NPT系综进行模拟.首先对整个体系在某一温度下进行 20 ps的弛豫,确保原子间受力平衡,然后进行 50 ps的弛豫对数据进行提取,最后以 25 K的升温间隔,对不同温度下的体系进行模拟.判断超离子相的方法是观察体系中原子的扩散情况,若体系中一种原子始终位于平衡位置,而另一种原子相对于参考原点存在较大位移,即原子离开平衡位置,在晶格间跳跃或在体系中自由移动,则可认为该体系进入了超离子相.Cazorla 等[37]通过计算均方位移判断了CaF2的超离子相,均方位移可以定义为

式中,ri(t)为原子i在t时刻的位置,t0为任意的时间原点,〈〉表示时间平均.分别计算每种原子的均方位移,并进行时间平均,即可得到整个体系中某一类原子的扩散情况.对于固体,体系的均方位移一般会随着时间的推移在最大值附近波动,若体系为液体,则均方位移会出现随时间均匀增大的现象.对GaN 体系中的Ga 原子与N 原子的均方位移分别进行了计算,以求找出某一温度下,一种原子在平衡位置附近振动,使得均方位移在最大值附近波动,而另一种原子突破晶格平衡位置的限制,在整个晶体内自由移动,均方位移呈随时间增大而上升的状态.

为了对比明显,本文在10 GPa的压力下分别计算了温度为2900 K、2925 K、3500 K 和4000 K时GaN的均方位移，如图5 所示.由于单相法存在过热现象,在10 GPa的压力下,温度达到 4000 K时GaN 仍没有熔化,因此才能观察到明显的超离子现象,但这并不代表GaN的熔化温度在 10 GPa时大于4000 K,本文中GaN的准确熔化温度为通过两相法计算的结果(见第3 部分).从图5 中可以看出,N 原子从 2900 K 到 4000 K的均方位移都没有出现随时间增大的趋势,2900 K 与2925 K 温度相近,曲线几乎重合,随着温度的升高,3500 K与4000 K的均方位移均在最大值附近波动.Ga 原子在2900 K时均方位移在最大值附近波动,2925 K 时出现了随时间增大的现象,这说明在该温度下,GaN 已经进入了超离子相,Ga 原子离开平衡位置,在晶格间跳跃或在晶体中自由移动,随着温度的增大,Ga 原子的均方位移斜率也出现了较大变化,说明 Ga 原子的扩散速率增大.图6为10 GPa 下 4000 K时的GaN 模拟体系,从原始体系中将 Ga 原子亚晶格与 N 原子亚晶格分离后,可以很明显看出,Ga 原子亚晶格出现了局部熔化,而 N 原子亚晶格仍具有稳定的结构,可以认为Ga原子在 N 原子构成的刚性框架内,以近乎液体的状态进行扩散,GaN 亚晶格的局部熔化使得 Ga原子有能力游离于整个晶体内部并起到导电作用,因此从固态相变到超离子相可能会使 GaN 从半导体变为导体.

图5 压力为10 GPa 时,GaN 在不同温度下的Ga 原子与 N 原子的均方位移,内插图为总览图Fig.5.Mean square displacement of the Ga and N atoms with different temperatures for GaN at 10 GPa,in which the inset is a general view.

图6 压力为10 GPa,温度为4000 K 时的GaN 原子构型,体系已经处于超离子相 (a) 原始构型;(b) Ga 原子亚晶格;(c) N 原子亚晶格Fig.6.GaN in superionic state at 10 GPa and 4000 K (a)original configuration,(b) Ga and (c) N atomic sublattices.

2.4 GaN的固-固相界线

克劳修斯-克拉珀龙方程用于描述单组分系统在相平衡时压力随温度的变化率,可以用于确定两相的相界线斜率,适用于本工作中确定 GaN的固-固相界线斜率,其定义为

式中T为温度,P为压力,ΔV为相变时的体积差,ΔH为相变潜热.本文使用克劳修斯-克拉珀龙方程确定了GaN 纤锌矿与岩盐结构的相界线、纤锌矿超离子相与岩盐结构的相界线以及纤锌矿超离子相与岩盐超离子相的相界线.由于不同固相的熔化曲线交点必然处于两相的固-固相界线上,为一个固-固-液三相共存点,在该点的温度和压力状态下,分别对 GaN的两种结构进行单相法分子动力学模拟,计算两种结构的体积差与焓差,代入(6)式即可得到固-固相界线在该三相共存点处的斜率.假设GaN的相界线在近距离上为一条直线,从该点出发,沿着斜率方向对下一个相平衡点进行短距离追踪,并在得到的新平衡点处再次进行两种结构的单相法模拟并计算克劳修斯-克拉珀龙斜率,以此类推即可得到 GaN 纤锌矿与岩盐结构的完整相界线.

3 结果与讨论

3.1 GaN 纤锌矿结构熔化曲线的“异常”

图7为纤锌矿结构 GaN 熔化温度的已有研究结果与本文计算结果的对比,对 GaN 纤锌矿结构的熔化温度研究结果差异显著,本文的分子动力学模拟结果处于已有的计算机模拟与实验结果之间,零压下GaN的熔化温度为2295 K,在 20 GPa 时为4170 K,这与 Porowski 等[14]的计算结果吻合,但与Van Vechten[2]的熔化模型存在显著区别.在实验研究方面,Utsumi 等[5]通过施加 6.0 GPa 以上的压力防止 GaN 分解,得到的熔化温度为2493 K,他们将压力提高到 6.4 GPa 和 6.8 GPa时,熔化温度也未出现明显变化;Sokol 等[12]在7.5 GPa的压力下对GaN 进行了熔化实验,研究发现当压力大于7.5 GPa 后 GaN 才会发生完全熔化;Porowski 等[14]的实验中使用了较大体积的GaN 样品,这能够降低实验结果的不确定性,他们发现在 6—9 GPa的压力范围内,温度达到 3400 K时也只能观察到 GaN的分解而不是熔化,这样的结果与 Utsumi 等[5]和 Sokol 等[12]的研究结果不符,但明显支持了Harafuji 等[6]和本工作的分子动力学模拟结果;Harafuji 等[6]的分子动力学模拟所使用的势函数为库仑项、短程排斥项、范德瓦耳斯项和共价项结合的形式,并采用两相法对 GaN的熔化曲线进行了计算,原子总数为9408 个,与本文构建两相共存体系的方法不同,Harafuji 等[6]先使GaN 体系在20 ps 内产生一个固液等分结构,随后将体系在不同温度下弛豫 60 ps,通过观察固液界面的移动来确定熔化温度,相比于Harafuji 等[6]的分子动力学模拟,本文延长了模拟中两相共存阶段的弛豫时间,使之达到了200 ps,更长的弛豫时间可以使整个GaN 体系更加接近平衡态,对提高模拟结果的可靠性有利,本文的势函数选择和两相法操作过程与 Harafuji 等[6]不同而导致结果出现了差异,但熔化曲线的趋势是一致的;Van Vechten[2]提出的GaN熔化曲线始终呈下降趋势,该结果是通过参考 Si的熔化温度随着压力升高而减小这一事实,结合经典电负性理论框架提出的,但之后的研究均发现 GaN的熔化曲线趋势与 Si 并不相同.大多数 Ⅲ-Ⅴ 族与 Ⅱ-Ⅵ 族半导体的熔化机制极为复杂,在熔化时常常伴随着配位数的改变,在环境条件下这些半导体均具有与 Si 晶体相似的正四面体键合方式,称为正四面体键半导体(Tetrahedrally Bonded Semiconductors,TBSs)[38],这种结构中的每个原子都与 4 个最近邻原子结合,形成四面体配位网格.大部分 TBSs的熔化温度会随着压力的增大而降低[2],这是因为TBSs的熔化机制由两个单独的过程组成: 一个是温度的增大导致 TBSs的配位数从 4 增大到 6,而这个转变温度会随着压力的增大而降低;另一个是在材料的自然熔化过程中,熔化温度会随着压力的增大而增大[38].两个过程的结合导致 Si,Ge 等 TBSs 直接从 4 配位固态熔化为密度更高的6 配位液态(短程有序,仍存在部分共价键),使得最终TBSs的熔化温度会随着压力的增大而降低,但 Harafuji 等[6]与 Porowski等[14]的计算结果表明 GaN的熔化温度会随着压力的增大而升高,这与其他 TBSs的情况并不相同.Porowski 等[38]在 2019 年对该问题进行了深入研究,他们对 Drozd-Rzoska 模型[14]进行线性外推,得到 GaN 在零压下 4 到 6 配位的理论转变温度为6000 K,远高于实际的熔化温度,因此会在4 配位固态下直接熔化,形成低密度的4 配位液态,4 配位固态的熔化温度会随着压力的增大而升高,因此出现了与其余 TBSs 不同的熔化情况.压力的增大使GaN的4 到 6 配位转变温度降低,当转变温度低于 4 配位固态的熔化温度时,GaN的熔化温度才会随着压力的增大而降低.如果仅从表象来看,GaN的这种现象与其他 TBSs相比是“异常”的,但实际上,Si,Ge 等 TBSs 在负压力下也可能会出现这种情况,从本质上来说,GaN的这种“异常”仅仅是 TBSs的熔化本质在正压下的体现.由于实验上实现负压非常困难,想要探索 Si,Ge 等 TBSs的熔化机制就必须克服负压问题,而 GaN的这种在正压下的“异常”熔化可能是探索 TBSs 熔化机制的一个捷径.

图7 GaN 纤锌矿结构熔化温度的已有研究结果与本文计算结果的对比Fig.7.Comparison between existing results on the melting temperature of GaN wurtzite structures and the results calculated in this paper.

在本文的分子动力学模拟中,GaN 纤锌矿结构在相变到岩盐结构之前都没有出现熔化温度的降低,这说明 GaN的4 到 6 配位转变温度在相变前均高于4 配位固态的熔化温度.另外,在分子动力学模拟中出现了与 Harafuji 等[6]相似的情况,即超离子相,最初认为超离子相可能是 GaN的6 配位液态,但结果中出现的超离子相仍具有固态结构,同时在GaN的岩盐结构中也出现了超离子相,这说明 GaN超离子相与 6 配位液态并不一样,是一种未被发现的新相,该相是在 GaN 熔化前出现的一个高温高压相,但目前为止并没有相关的实验证据表明 GaN 存在该相,这需要相关工作者继续进行研究.

3.2 GaN 纤锌矿-岩盐结构的相界线

图8为使用第一性原理计算得到的GaN 3 种结构(纤锌矿、闪锌矿和岩盐结构)相对焓随压力的变化关系,在 43.6 GPa 之前,3 种结构的稳定性顺序为: 纤锌矿结构 ＞ 闪锌矿结构 ＞ 岩盐结构.GaN 闪锌矿与岩盐结构的相对焓相交于43.6 GPa,与 Sun 等[39]的计算结果 48 GPa 相近,在这之后岩盐结构的稳定性大于闪锌矿结构,纤锌矿与岩盐结构的相对焓相交于 44.3 GPa,之后岩盐结构成为最稳定的高压相,在 0—80 GPa 压力范围内GaN 纤锌矿与闪锌矿结构的相对焓均无交点,这表明在零温下不会发生纤锌矿到闪锌矿结构的相变.图9为使用分子动力学方法对GaN 纤锌矿与岩盐结构固-固相界线的追踪结果与已知研究结果的对比,目前对 GaN 固-固相变压力的研究结果有限且存在很大分歧,Zhou 等[3]使用了第一性原理计算结合准谐近似的方法,得到零温下GaN 纤锌矿到岩盐结构的相变压力为32.4 GPa,准谐近似方法弥补了第一性原理计算无法得到非零温度下材料性质的缺陷,但由于该方法仅考虑低阶的非谐效应,因而在涉及到更高温度的相变压力计算中难以得到可靠的结果,这使得 Zhou 等[3]在高温下得到的相变压力比已有结果低,但零温下的相变压力与 Xia 等[8]的实验结果和 Saitta 等[40]的结果相近;Sadovyi 等[4]采用第一性原理计算结合自洽声子方法对GaN 纤锌矿结构和岩盐结构的相界线进行了计算,该自洽声子方法引入了 6 阶非谐力,这对计算结果的准确性有较大提升,Sadovyi 等[4]在零温下得到的GaN 纤锌矿到岩盐结构的相变压力为45.0 GPa,与已有结果吻合[11,41];Christensen等[42]与等[43]对GaN 纤锌矿和岩盐结构相变压力的第一性原理计算结果分别为51.8 GPa和56.06 GPa,与Ueno 等[10]的实验值 52.2 GPa吻合,但略高于本文的计算结果.本文在使用分子动力学模拟对 GaN 相变压力进行计算时,采用了经过验证的有效势函数,结合克劳修斯-克拉珀龙方程,并采取尽可能小的压力增幅对GaN的相界线进行逐步追踪计算,得到较低温度下纤锌矿到岩盐结构的相变压力为45.9 GPa,这与已有的实验结果[41]、计算结果[4,11]以及本工作的第一性原理计算结果均吻合,但在相同压力下的相转变温度比Zhou 等[3]和 Sadovyi 等[4]的结果高.Sadovyi 等[4]推断 GaN 岩盐结构的熔化温度在 37 GPa 时应该略高于 2100 K,而本文的计算结果为4900 K,高的熔化温度使固-固-液三相共存点位置更高,因此从该点追踪到的相界线必然会有高的相转变温度.综合来看,研究中计算方法的选择不同是导致GaN 结构相变压力分歧较大的主要原因.

图8 第一性原理计算得到的GaN 纤锌矿、闪锌矿与岩盐三种结构的相对焓,内插图为放大的相对焓交点Fig.8.Relative enthalpies of GaN with wurtzite,zincblende,and rocksalt structures are calculated by first-principles calculations,and the interpolation shows the enlarged relative enthalpies intersection.

图9 GaN 纤锌矿与岩盐结构固-固相界线的已有数据与本文计算结果的对比Fig.9.Comparison of the existing solid-solid phase boundary results of GaN with wurtzite and rocksalt structures with the results calculated in this paper.

3.3 GaN的P-T 相图

图10为GaN 纤锌矿结构与岩盐结构在0—80 GPa 压力范围内的P-T相图,图中给出了GaN的5 种存在状态: 固态纤锌矿结构、纤锌矿超离子相、固态岩盐结构、岩盐超离子相和液态,其中固态纤锌矿结构、固态岩盐结构与液态存在的温度范围较宽,两个超离子相存在的温度范围较窄.通过外推纤锌矿结构的熔化曲线,可以得到GaN 在环境条件下的熔化温度为2295 K,随着压力的增大,GaN 纤锌矿结构的熔化温度逐渐升高,并于 P1 点与岩盐结构的熔化曲线相交,该点为GaN的液态、纤锌矿超离子相、岩盐超离子相的三相共存点;从 P1 点开始,GaN 岩盐结构的熔化温度随着压力的增大而逐渐升高,相比于纤锌矿结构的熔化曲线,岩盐结构的熔化曲线更加陡峭,在 80 GPa 时,熔化温度已经达到了 7646 K;GaN 纤锌矿结构与其超离子相的相界线为一直线,于 P2 点与纤锌矿结构的熔化曲线相交,该点为液态、纤锌矿超离子相、固态纤锌矿结构的三相共存点,纤锌矿结构与超离子相的相界线并没有与纵坐标相交,说明GaN 在常压下不存在超离子相,该相为一高温高压相,只有在 2.0 GPa 之后,且在高温的驱动下,才进入超离子相;GaN 岩盐结构的超离子相转变温度非常接近熔化温度,但随着压力的增大而逐渐远离熔化温度,在 80 GPa 下,岩盐结构超离子相的转变温度为6162 K.对GaN 固-固相界线的追踪从 P1 点开始,从 P1 点得到的克劳修斯-克拉珀龙斜率为331.9 K/GPa,在该斜率方向上压力温度条件为32.6 GPa,4321 K 处作为新起点,得到第2 个斜率: —275.8 K/GPa,此时的克劳修斯-克拉珀龙斜率由正变负,沿该方向的GaN 相界线于 P3 点与岩盐超离子相的相界线相交,P3 点为纤锌矿超离子相、岩盐超离子相与固态岩盐结构的三相共存点;以 P3 点作为新起点继续进行固-固相界线的追踪,在点 P4 处GaN的固-固相界线与纤锌矿超离子相的相界线相交,该点为固态纤锌矿结构、纤锌矿超离子相、固态岩盐结构的三相共存点,克劳修斯-克拉珀龙斜率为—81.2 K/GPa,为相界线中克劳修斯-克拉珀龙斜率绝对值的最小值;从 P4 点开始再无特殊的三相共存点,只需一步步进行相界线的追踪,并适当调整每次追踪的距离,即可得到完整的GaN 固态纤锌矿与固态岩盐结构的相界线.在本文的计算结果中,GaN的相界线不是一条直线,这在固态纤锌矿结构与固态岩盐结构的相界线中尤为明显,在压力大于 40 GPa 后,相界线斜率开始突然增大,并于 45.9 GPa 与x轴相交.相图中特殊的三相共存点位置见表3.

图10 GaN 0—80 GPa 压力范围内的P-T 相图,图中的所有点均为分子动力学模拟值,黑色和红色线为GaN纤锌矿结构、岩盐结构熔化温度的拟合曲线;海蓝色与紫色线为GaN 纤锌矿超离子相、岩盐超离子相的转变边界拟合曲线;蓝色、绿色与橙色点线分别为GaN 纤锌矿超离子相-岩盐超离子相、纤锌矿超离子相-固态岩盐结构、固态纤锌矿结构-固态岩盐结构的相转变边界Fig.10.P-T phase diagram of GaN in the pressure range of 0—80 GPa,WZ and RS are used to denote the wurtzite and rocksalt structures of GaN in the diagram.All points are molecular dynamics simulations values in the diagram,the black and red lines are the fitted curves of melting temperature of the GaN with wurtzite and rocksalt structures.The navy blue and purple lines are the fitted curves of phase transition boundary of wurtzite superionic and rocksalt superionic for GaN.The blue,green and orange dotted lines are the phase transition boundary of GaN with wurtzite phase superionic-rocksalt phase superionic,wurtzite phase superionic-rocksalt solid and wurtzite solid-rocksalt solid,respectively.

表3 GaN P-T 相图中三相共存点的位置Table 3.Location of three-phase coexistence points in GaN P-T phase diagram.

4 结论

本文采用基于 Morse 与 Born-Mayer 组合势的分子动力学模拟方法,对 GaN 纤锌矿与岩盐结构在宽广的温度压力范围内的P-T相图进行了可靠预测.首先使用第一性原理计算方法对 GaN 纤锌矿结构和岩盐结构的晶格常数、弹性常数、体积模量与晶格能进行了计算,然后与使用基于已选势函数的晶格动力学方法计算的结果进行了对比,该势函数能够给出与第一性原理计算和实验结果相近的力学、热力学性质与晶格能,相近的晶格能说明该势函数能够很好地再现 GaN 在高温下尤其是熔化后的性质.使用分子动力学模拟方法计算了300 K 下 GaN 纤锌矿结构与岩盐结构体积比率随压力的变化情况,在纤锌矿结构下,体积比率随压力的变化情况与实验数据符合良好,相变为岩盐结构时的体积变化率为14.4%;计算了GaN 纤锌矿结构与岩盐结构的熔化曲线,零压下纤锌矿结构的熔化温度为2295 K,未发现熔化温度随压力增大而降低的情况;对 GaN 纤锌矿与岩盐结构的固-固相界线进行了追踪,零温下纤锌矿到岩盐结构的相变压力为45.9 GPa,与第一性原理计算结果44.3 GPa 符合良好;对 GaN 可能存在的超离子相进行了确定,在该状态下,GaN的Ga 原子亚晶格局部熔化,Ga 原子能够在晶体内部自由移动并起导电作用,纤锌矿结构在压力大于2.0 GPa 且温度高于 2550 K 时进入超离子相,而岩盐结构在压力温度大于 33.1 GPa 和 4182 K 后发生超离子相转变.