亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        氮化鎵相圖預(yù)測(cè)及其高壓熔化特性研究*

        2022-10-16 09:23:40雷振帥孫小偉劉子江宋婷田俊紅
        物理學(xué)報(bào) 2022年19期
        關(guān)鍵詞:鋅礦巖鹽勢(shì)函數(shù)

        雷振帥 孫小偉 劉子江 宋婷 田俊紅

        (蘭州交通大學(xué)數(shù)理學(xué)院,蘭州 730070)

        采用經(jīng)典分子動(dòng)力學(xué)模擬,結(jié)合第一性原理計(jì)算及晶格動(dòng)力學(xué)方法對(duì)氮化鎵(GaN)纖鋅礦結(jié)構(gòu)與巖鹽結(jié)構(gòu)在 0—80 GPa 壓力范圍內(nèi)的相圖進(jìn)行了預(yù)測(cè).第一性原理計(jì)算與分子動(dòng)力學(xué)模擬得到的零溫下 GaN 纖鋅礦到巖鹽結(jié)構(gòu)的相變壓力分別為44.3 GPa 和 45.9 GPa,與已有研究的實(shí)驗(yàn)結(jié)果吻合;通過(guò)外推纖鋅礦結(jié)構(gòu) GaN的熔化曲線得到其零壓下的熔化溫度為2295 K,當(dāng)壓力增大到33.3 GPa 時(shí),纖鋅礦結(jié)構(gòu)熔化曲線與巖鹽結(jié)構(gòu)熔化曲線相交,兩種結(jié)構(gòu)的熔化溫度均隨壓力的增大而上升;GaN 還可能存在超離子相,纖鋅礦結(jié)構(gòu)在壓力大于2.0 GPa 且溫度高于 2550 K 時(shí)發(fā)生超離子相轉(zhuǎn)變,巖鹽結(jié)構(gòu)在壓力溫度大于 33.1 GPa 和4182 K 后發(fā)生超離子相轉(zhuǎn)變,二者的相轉(zhuǎn)變溫度均會(huì)隨著壓力的增大而升高;GaN 纖鋅礦和巖鹽結(jié)構(gòu)的相界線并非為一直線,在高溫下相界線斜率為正,隨著溫度的降低逐漸變?yōu)橐粭l具有負(fù)斜率的曲線.

        1 引言

        Ⅲ-Ⅴ族化合物半導(dǎo)體氮化鎵(GaN)具有禁帶寬度大、熱導(dǎo)率高、電子飽和速度快等優(yōu)異特性,是發(fā)展高頻、高功率電子器件的優(yōu)良半導(dǎo)體材料,在發(fā)光二極管、晶體管與軍事電子設(shè)備中均具有廣泛的應(yīng)用前景[1].然而,相比于微電子領(lǐng)域,人們從基礎(chǔ)性質(zhì)的角度給予 GaN的研究還不夠充分,如高溫下固-固相變的克勞修斯-克拉珀龍斜率與熔化溫度仍存在很大爭(zhēng)議[2-6],究其原因,主要是因?yàn)镚aN的分解溫度低于熔化溫度,在熔化時(shí)具有非常高的氮?dú)馄胶庹魵鈮篬7],在進(jìn)行GaN的高溫實(shí)驗(yàn)時(shí),需要保持高的氮?dú)鈮毫σ苑乐蛊浞纸?這使得實(shí)驗(yàn)研究變得尤為困難,在一定程度上限制了GaN的開發(fā)和應(yīng)用.因此,GaN的固-固相變和熔化溫度等基礎(chǔ)性質(zhì)成為了該材料發(fā)展過(guò)程中的重要研究?jī)?nèi)容.

        理論和實(shí)驗(yàn)研究表明,GaN 具有3 種典型結(jié)構(gòu),即纖鋅礦結(jié)構(gòu)、閃鋅礦結(jié)構(gòu)和巖鹽結(jié)構(gòu)[3,4,8].在環(huán)境條件下GaN 一般以纖鋅礦結(jié)構(gòu)穩(wěn)定存在,閃鋅礦結(jié)構(gòu)可通過(guò)氣相外延和分子束外延沉積在立方(001)襯底上,以薄膜形式穩(wěn)定存在于環(huán)境條件中[9],巖鹽結(jié)構(gòu)為穩(wěn)定的高壓相.GaN 纖鋅礦到巖鹽結(jié)構(gòu)的相變發(fā)生在 37—52 GPa 之間[8,10],而閃鋅礦到巖鹽結(jié)構(gòu)的相變發(fā)生在 36—42 GPa 之間[11].GaN 纖鋅礦和閃鋅礦結(jié)構(gòu)均為Ga-N 正四面體,二者的結(jié)構(gòu)相近,實(shí)驗(yàn)上還未觀察到該相變,但也有研究表明該相變發(fā)生在 10.87 GPa[11].GaN的高溫高壓實(shí)驗(yàn)條件苛刻,而可以計(jì)算高溫相變的分子模擬又缺乏準(zhǔn)確的原子力場(chǎng),因此,對(duì) GaN 進(jìn)行高溫高壓相圖的研究是一項(xiàng)困難的工作.Van Vechten[2]提出的化學(xué)鍵模型預(yù)測(cè)了 GaN 在高溫下的固-固相變壓力,但該模型預(yù)測(cè)的是纖鋅礦結(jié)構(gòu)與β-tin 結(jié)構(gòu)的相變.關(guān)于 GaN 纖鋅礦到巖鹽結(jié)構(gòu)的相界線研究較少,僅有 Zhou 等[3]的數(shù)值模擬和 Sadovyi 等[4]的高溫高壓實(shí)驗(yàn)結(jié)果,后者對(duì)GaN的固-固相變行為進(jìn)行了較為全面的分析.

        GaN 熔化溫度的研究結(jié)果存在較大分歧,Van Vechten[2]提出的模型預(yù)測(cè)了常壓下GaN的熔化溫度為2791 K,與之后的實(shí)驗(yàn)結(jié)果比較接近,但該模型預(yù)測(cè)的熔化曲線呈一次函數(shù)下降趨勢(shì),這與之后的許多研究結(jié)果不符;Karpiński 等[7]在高壓碳化鎢砧槽中進(jìn)行了高溫實(shí)驗(yàn),發(fā)現(xiàn)壓力為6.0 GPa時(shí) GaN的熔點(diǎn)高于 2573 K;Utsumi 等[5]的X 射線衍射實(shí)驗(yàn)給出了壓力大于6.0 GPa 時(shí),GaN的熔化溫度為2493 K,隨后 Sokol 等[12]與Saitoh 等[13]的研究也支持了Utsumi 等[5]的結(jié)果;Porowski 等[14]為提高高溫高壓下對(duì)樣品特征的檢測(cè)和分析質(zhì)量,使用了比早期實(shí)驗(yàn)材料體積大了100 倍的GaN 樣品,該研究發(fā)現(xiàn)在 6—9 GPa的壓力范圍內(nèi),溫度達(dá)到 3400 K 時(shí)也只是觀察到 GaN的分解而不是熔化,這與 Harafuji 等[6]與Nord 等[15]的分子動(dòng)力學(xué)模擬結(jié)果相近.值得注意的是,Harafuji 等[6]的研究中出現(xiàn)了GaN 熔化溫度的不連續(xù)現(xiàn)象,作者將其解釋為帶部分電荷的Ga 原子的固體電解質(zhì)行為,在本工作中也發(fā)現(xiàn)了Ga 原子的類似行為(超離子相),但并沒(méi)有出現(xiàn)熔化溫度的不連續(xù)現(xiàn)象.

        本文采用Born-Mayer 與Morse 形式結(jié)合的相互作用勢(shì)進(jìn)行經(jīng)典分子動(dòng)力學(xué)模擬,并通過(guò)基于密度泛函理論的第一性原理計(jì)算方法與晶格動(dòng)力學(xué)方法,對(duì)GaN 晶體纖鋅礦與巖鹽結(jié)構(gòu)的彈性常數(shù)、體積模量與晶格能進(jìn)行了計(jì)算.在此基礎(chǔ)上對(duì)GaN 在高壓下的體積壓縮行為進(jìn)行了研究,同時(shí)采用兩相法和單相法預(yù)測(cè)了不同結(jié)構(gòu)GaN的熔化溫度與超離子相,最后給出了GaN 在寬廣的溫度和壓力范圍內(nèi)的P-T相圖.

        2 模型及方法

        2.1 GaN 勢(shì)函數(shù)及其驗(yàn)證

        Alder 和Wainwright[16]于 1957 年提出分子動(dòng)力學(xué)方法,該方法是一種原子層次的計(jì)算機(jī)模擬方法,由于不涉及電子而降低了計(jì)算量,可用于模擬龐大的原子體系.分子動(dòng)力學(xué)模擬的核心是準(zhǔn)確描述原子間相互作用的勢(shì)函數(shù),現(xiàn)有的勢(shì)函數(shù)種類可以分為描述離子鍵的二體勢(shì)、描述共價(jià)鍵的三體勢(shì)以及描述金屬鍵的多體勢(shì).GaN 已有的勢(shì)函數(shù)非常豐富,如三體 Tersoff 勢(shì)[15]與 Stillinger-Weber勢(shì)[17,18]、多體修正嵌入原子勢(shì)[19]以及對(duì)勢(shì)類型的殼層模型Buckingham 勢(shì)[20]與Morse 勢(shì)[21]等,這些勢(shì)函數(shù)主要用于 GaN的缺陷、力學(xué)性質(zhì)等溫度較低的計(jì)算,在高溫高壓條件下可能不適用,因此需要對(duì)勢(shì)函數(shù)進(jìn)行篩選,本文最終確定使用2005 年Zhang 等[21]通過(guò)晶格反演方法得到的對(duì)勢(shì).GaN具有混合共價(jià)-離子鍵[15],使用二體勢(shì)或三體勢(shì)均可進(jìn)行描述.圖1為利用基于密度泛函理論的第一性原理計(jì)算方法得到的GaN 纖鋅礦和巖鹽結(jié)構(gòu)的電子局域函數(shù)圖,由圖1 可以看出兩種結(jié)構(gòu)的電子幾乎都局域在 N 原子周圍,離子鍵更加明顯,因此對(duì)勢(shì)形式的勢(shì)函數(shù)可能會(huì)更好地描述GaN 各原子對(duì)間的相互作用:

        圖1 GaN 電子局域函數(shù)圖 (a) 纖鋅礦結(jié)構(gòu)(110)晶面;(b) 巖鹽結(jié)構(gòu)(100)晶面Fig.1.Electron localization function of GaN (a) wurtzite structure (110) and (b) rocksalt structure (100) crystal planes.

        (1)式中第1 項(xiàng)為庫(kù)侖項(xiàng),Zi,Zj分別為i原子與j原子的有效電荷,ε0為真空介電常數(shù),r為原子間距,第2 項(xiàng)為i原子與j原子的短程排斥項(xiàng).(2)式為N-N,Ga-Ga 間短程排斥相互作用,使用排斥指數(shù)形式的勢(shì)函數(shù)描述,其中D,γ,R分別為勢(shì)函數(shù)參數(shù);(3)式為Ga-N 間短程排斥相互作用,使用 Morse 形式的勢(shì)函數(shù)描述.排斥指數(shù)形式的勢(shì)函數(shù)在分子動(dòng)力學(xué)模擬中使用較少,而 Born-Mayer 勢(shì)函數(shù)[22]則使用的更加廣泛,這是因?yàn)槠湫问胶?jiǎn)單且參數(shù)更少,在計(jì)算材料熱力學(xué)性質(zhì)時(shí)往往也能給出可靠的計(jì)算結(jié)果,因此,本文將 Zhang等[21]的排斥指數(shù)形式的勢(shì)函數(shù)擬合到了 Born-Mayer 形式:

        式中A,η為勢(shì)函數(shù)參數(shù).使用該形式描述 N-N,Ga-Ga 間的短程相互作用,Ga-N 間的相互作用仍使用 Morse 形式描述,最終的勢(shì)函數(shù)參數(shù)見表1.

        表1 GaN的Morse 勢(shì)參數(shù)與本文擬合得到的Born-Mayer 勢(shì)參數(shù)Table 1.The Morse potential parameters and fitted Born-Mayer potential parameters of GaN.

        為驗(yàn)證勢(shì)函數(shù)的準(zhǔn)確性,本文分別使用第一性原理計(jì)算方法與晶格動(dòng)力學(xué)方法對(duì) GaN 不同結(jié)構(gòu)的彈性性質(zhì)與晶格能進(jìn)行了計(jì)算.第一性原理計(jì)算中,為使體系能量在完備平面波基矢水平上足夠收斂,采用了 BFGS 算法[23]對(duì)晶體結(jié)構(gòu)進(jìn)行幾何優(yōu)化,以獲得局域能量最低結(jié)構(gòu),選取由 Perdew 等修訂的PBE 形式的廣義梯度近似方法[24]描述電子間的交換關(guān)聯(lián)能,選取非局域超軟贗勢(shì)[25]描述離子實(shí)和價(jià)電子間的相互作用,計(jì)算過(guò)程中 Ga 和N 原子的外層電子組態(tài)分別為3d104s24p1和 2s22p3,平面波截?cái)嗄芫鶠?50 eV,GaN 纖鋅礦與巖鹽結(jié)構(gòu)的倒易空間布里淵區(qū)k點(diǎn)分別取值13×13×7和9×9×9,在結(jié)構(gòu)優(yōu)化過(guò)程中體系能量收斂標(biāo)準(zhǔn)為5×10—6eV/atom,優(yōu)化后作用在晶胞中每個(gè)原子上的力小于 0.01 eV/?,晶胞應(yīng)力偏差低于0.02 GPa,公差偏移小于 5×10—4?.計(jì)算彈性常數(shù)時(shí),最大應(yīng)變振幅設(shè)置為0.003,每個(gè)應(yīng)變循環(huán)6 步,即產(chǎn)生 6 個(gè)扭曲結(jié)構(gòu),所有的第一性原理計(jì)算均利用 CASTEP 軟件包[26]完成.晶格動(dòng)力學(xué)計(jì)算使用表1 中的勢(shì)函數(shù)參數(shù),截?cái)喟霃竭x擇 10 ?,所有的晶格動(dòng)力學(xué)計(jì)算均利用 GULP 軟件包[27]完成.得到的GaN 纖鋅礦與巖鹽結(jié)構(gòu)的相關(guān)參數(shù)見表2,通過(guò)晶格動(dòng)力學(xué)方法計(jì)算的纖鋅礦結(jié)構(gòu)的晶格常數(shù)、彈性常數(shù)與體積模量接近第一性原理計(jì)算和實(shí)驗(yàn)結(jié)果,這表明該勢(shì)函數(shù)能夠準(zhǔn)確描述GaN 纖鋅礦結(jié)構(gòu)的力學(xué)性質(zhì),但巖鹽結(jié)構(gòu)的彈性常數(shù)與第一性原理計(jì)算數(shù)據(jù)存在差距,晶格動(dòng)力學(xué)計(jì)算的兩種結(jié)構(gòu)的晶格能與已有數(shù)據(jù)非常接近,這表明該勢(shì)函數(shù)能夠很好地描述 GaN 在高溫下尤其是熔化后的性質(zhì).為驗(yàn)證勢(shì)函數(shù)在高壓下的有效性,使用基于該勢(shì)函數(shù)的分子動(dòng)力學(xué)方法計(jì)算了GaN 纖鋅礦與巖鹽結(jié)構(gòu)在 300 K 時(shí)的體積比率隨壓力的變化情況,如圖2 所示,GaN 纖鋅礦結(jié)構(gòu)的體積比率隨壓力的變化與已有的實(shí)驗(yàn)與計(jì)算結(jié)果均吻合,在 45.9 GPa 相變到巖鹽結(jié)構(gòu)時(shí),體積變化率為14.4%,比Xia 等[8]的實(shí)驗(yàn)結(jié)果(17.9%)低,但接近于Pandey 等[28]的從頭算結(jié)果,說(shuō)明該勢(shì)函數(shù)能夠描述 GaN 在高壓條件下的熱力學(xué)狀態(tài).

        圖2 GaN 纖鋅礦和巖鹽結(jié)構(gòu)在 300 K 下的體積比率與已有結(jié)果對(duì)比,所有數(shù)據(jù)均歸一化至纖鋅礦結(jié)構(gòu)初始體積,紅色與黑色圓點(diǎn)為分子動(dòng)力學(xué)模擬結(jié)果(實(shí)線為擬合曲線);藍(lán)色正三角形與綠色菱形分別為Xia 等[8]與 Ueno 等[10]的X 射線衍射實(shí)驗(yàn)結(jié)果;洋紅色虛線為Pandey 等[28]的從頭算結(jié)果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 纖鋅礦與巖鹽結(jié)構(gòu)在零壓下的晶格常數(shù),彈性常數(shù) 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的熔化溫度

        為保證原子數(shù)目相同,本文分別使用8×8×16與6×8×16的正交體系進(jìn)行 GaN 纖鋅礦與巖鹽結(jié)構(gòu)的分子動(dòng)力學(xué)模擬.兩個(gè)體系均包含6144個(gè)原子,采用周期性邊界條件,時(shí)間步長(zhǎng)1 fs,所有的分子動(dòng)力學(xué)計(jì)算均利用 LAMMPS 軟件[32]完成.目前分子動(dòng)力學(xué)模擬中計(jì)算熔化溫度的方法主要有單相法、兩相法、孔洞法、Z 方法以及改進(jìn)的Z 方法[33-36],本文選擇兩相法進(jìn)行 GaN的熔化溫度計(jì)算,該方法準(zhǔn)確度很高,而且也能在一定程度上控制壓力,在進(jìn)行熔化曲線的計(jì)算時(shí),具有較好的操作空間.兩相法計(jì)算熔化溫度的關(guān)鍵在于構(gòu)建能夠長(zhǎng)時(shí)間存在的兩相共存體系,為達(dá)到這個(gè)目的,首先將整個(gè)模擬體系在等溫等壓系綜(NPT)中以略低于熔點(diǎn)的溫度弛豫30 ps,保證整個(gè)體系的應(yīng)力平衡,然后將體系沿z方向平均分為兩部分,一部分設(shè)置受力為0,將其固定,另一部分在遠(yuǎn)高于熔點(diǎn)的溫度下弛豫 50 ps,保證其完全熔化,再將完全熔化的這部分原子降溫到略高于熔點(diǎn)的溫度,降溫時(shí)間為10 ps,最后釋放固定的原子,將整個(gè)體系置于微正則系綜(NVE)中自由弛豫,形成固液共存體系,若整個(gè)體系能夠以固液共存狀態(tài)長(zhǎng)時(shí)間存在(本文中為200 ps),且固液部分體積無(wú)顯著變化,則可以認(rèn)為此時(shí)體系的溫度為最終熔化溫度.圖3為構(gòu)建的GaN 兩相共存體系,圖4為對(duì)應(yīng)的原子數(shù)密度,可以清楚地看出體系明顯分為兩部分,數(shù)密度波動(dòng)大的部分處于固態(tài),數(shù)密度波動(dòng)較小的部分處于液態(tài).

        圖3 處于兩相共存狀態(tài)的GaN 體系Fig.3.GaN system in the two-phase coexistence state.

        圖4 兩相共存狀態(tài)下 GaN 體系沿 z 軸方向的原子數(shù)密度Fig.4.Atomic number density along the z-direction of the GaN system in the two-phase coexistence state.

        2.3 GaN的超離子相

        超離子相的確定不涉及過(guò)熱問(wèn)題,因此本文使用單相法對(duì) GaN的超離子相進(jìn)行確定.為準(zhǔn)確控制壓力和溫度,選擇使用NPT系綜進(jìn)行模擬.首先對(duì)整個(gè)體系在某一溫度下進(jìn)行 20 ps的弛豫,確保原子間受力平衡,然后進(jìn)行 50 ps的弛豫對(duì)數(shù)據(jù)進(jìn)行提取,最后以 25 K的升溫間隔,對(duì)不同溫度下的體系進(jìn)行模擬.判斷超離子相的方法是觀察體系中原子的擴(kuò)散情況,若體系中一種原子始終位于平衡位置,而另一種原子相對(duì)于參考原點(diǎn)存在較大位移,即原子離開平衡位置,在晶格間跳躍或在體系中自由移動(dòng),則可認(rèn)為該體系進(jìn)入了超離子相.Cazorla 等[37]通過(guò)計(jì)算均方位移判斷了CaF2的超離子相,均方位移可以定義為

        式中,ri(t)為原子i在t時(shí)刻的位置,t0為任意的時(shí)間原點(diǎn),〈〉表示時(shí)間平均.分別計(jì)算每種原子的均方位移,并進(jìn)行時(shí)間平均,即可得到整個(gè)體系中某一類原子的擴(kuò)散情況.對(duì)于固體,體系的均方位移一般會(huì)隨著時(shí)間的推移在最大值附近波動(dòng),若體系為液體,則均方位移會(huì)出現(xiàn)隨時(shí)間均勻增大的現(xiàn)象.對(duì)GaN 體系中的Ga 原子與N 原子的均方位移分別進(jìn)行了計(jì)算,以求找出某一溫度下,一種原子在平衡位置附近振動(dòng),使得均方位移在最大值附近波動(dòng),而另一種原子突破晶格平衡位置的限制,在整個(gè)晶體內(nèi)自由移動(dòng),均方位移呈隨時(shí)間增大而上升的狀態(tài).

        為了對(duì)比明顯,本文在10 GPa的壓力下分別計(jì)算了溫度為2900 K、2925 K、3500 K 和4000 K時(shí)GaN的均方位移,如圖5 所示.由于單相法存在過(guò)熱現(xiàn)象,在10 GPa的壓力下,溫度達(dá)到 4000 K時(shí)GaN 仍沒(méi)有熔化,因此才能觀察到明顯的超離子現(xiàn)象,但這并不代表GaN的熔化溫度在 10 GPa時(shí)大于4000 K,本文中GaN的準(zhǔn)確熔化溫度為通過(guò)兩相法計(jì)算的結(jié)果(見第3 部分).從圖5 中可以看出,N 原子從 2900 K 到 4000 K的均方位移都沒(méi)有出現(xiàn)隨時(shí)間增大的趨勢(shì),2900 K 與2925 K 溫度相近,曲線幾乎重合,隨著溫度的升高,3500 K與4000 K的均方位移均在最大值附近波動(dòng).Ga 原子在2900 K時(shí)均方位移在最大值附近波動(dòng),2925 K 時(shí)出現(xiàn)了隨時(shí)間增大的現(xiàn)象,這說(shuō)明在該溫度下,GaN 已經(jīng)進(jìn)入了超離子相,Ga 原子離開平衡位置,在晶格間跳躍或在晶體中自由移動(dòng),隨著溫度的增大,Ga 原子的均方位移斜率也出現(xiàn)了較大變化,說(shuō)明 Ga 原子的擴(kuò)散速率增大.圖6為10 GPa 下 4000 K時(shí)的GaN 模擬體系,從原始體系中將 Ga 原子亞晶格與 N 原子亞晶格分離后,可以很明顯看出,Ga 原子亞晶格出現(xiàn)了局部熔化,而 N 原子亞晶格仍具有穩(wěn)定的結(jié)構(gòu),可以認(rèn)為Ga原子在 N 原子構(gòu)成的剛性框架內(nèi),以近乎液體的狀態(tài)進(jìn)行擴(kuò)散,GaN 亞晶格的局部熔化使得 Ga原子有能力游離于整個(gè)晶體內(nèi)部并起到導(dǎo)電作用,因此從固態(tài)相變到超離子相可能會(huì)使 GaN 從半導(dǎo)體變?yōu)閷?dǎo)體.

        圖5 壓力為10 GPa 時(shí),GaN 在不同溫度下的Ga 原子與 N 原子的均方位移,內(nèi)插圖為總覽圖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 時(shí)的GaN 原子構(gòu)型,體系已經(jīng)處于超離子相 (a) 原始構(gòu)型;(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的固-固相界線

        克勞修斯-克拉珀龍方程用于描述單組分系統(tǒng)在相平衡時(shí)壓力隨溫度的變化率,可以用于確定兩相的相界線斜率,適用于本工作中確定 GaN的固-固相界線斜率,其定義為

        式中T為溫度,P為壓力,ΔV為相變時(shí)的體積差,ΔH為相變潛熱.本文使用克勞修斯-克拉珀龍方程確定了GaN 纖鋅礦與巖鹽結(jié)構(gòu)的相界線、纖鋅礦超離子相與巖鹽結(jié)構(gòu)的相界線以及纖鋅礦超離子相與巖鹽超離子相的相界線.由于不同固相的熔化曲線交點(diǎn)必然處于兩相的固-固相界線上,為一個(gè)固-固-液三相共存點(diǎn),在該點(diǎn)的溫度和壓力狀態(tài)下,分別對(duì) GaN的兩種結(jié)構(gòu)進(jìn)行單相法分子動(dòng)力學(xué)模擬,計(jì)算兩種結(jié)構(gòu)的體積差與焓差,代入(6)式即可得到固-固相界線在該三相共存點(diǎn)處的斜率.假設(shè)GaN的相界線在近距離上為一條直線,從該點(diǎn)出發(fā),沿著斜率方向?qū)ο乱粋€(gè)相平衡點(diǎn)進(jìn)行短距離追蹤,并在得到的新平衡點(diǎn)處再次進(jìn)行兩種結(jié)構(gòu)的單相法模擬并計(jì)算克勞修斯-克拉珀龍斜率,以此類推即可得到 GaN 纖鋅礦與巖鹽結(jié)構(gòu)的完整相界線.

        3 結(jié)果與討論

        3.1 GaN 纖鋅礦結(jié)構(gòu)熔化曲線的“異?!?/h3>

        圖7為纖鋅礦結(jié)構(gòu) GaN 熔化溫度的已有研究結(jié)果與本文計(jì)算結(jié)果的對(duì)比,對(duì) GaN 纖鋅礦結(jié)構(gòu)的熔化溫度研究結(jié)果差異顯著,本文的分子動(dòng)力學(xué)模擬結(jié)果處于已有的計(jì)算機(jī)模擬與實(shí)驗(yàn)結(jié)果之間,零壓下GaN的熔化溫度為2295 K,在 20 GPa 時(shí)為4170 K,這與 Porowski 等[14]的計(jì)算結(jié)果吻合,但與Van Vechten[2]的熔化模型存在顯著區(qū)別.在實(shí)驗(yàn)研究方面,Utsumi 等[5]通過(guò)施加 6.0 GPa 以上的壓力防止 GaN 分解,得到的熔化溫度為2493 K,他們將壓力提高到 6.4 GPa 和 6.8 GPa時(shí),熔化溫度也未出現(xiàn)明顯變化;Sokol 等[12]在7.5 GPa的壓力下對(duì)GaN 進(jìn)行了熔化實(shí)驗(yàn),研究發(fā)現(xiàn)當(dāng)壓力大于7.5 GPa 后 GaN 才會(huì)發(fā)生完全熔化;Porowski 等[14]的實(shí)驗(yàn)中使用了較大體積的GaN 樣品,這能夠降低實(shí)驗(yàn)結(jié)果的不確定性,他們發(fā)現(xiàn)在 6—9 GPa的壓力范圍內(nèi),溫度達(dá)到 3400 K時(shí)也只能觀察到 GaN的分解而不是熔化,這樣的結(jié)果與 Utsumi 等[5]和 Sokol 等[12]的研究結(jié)果不符,但明顯支持了Harafuji 等[6]和本工作的分子動(dòng)力學(xué)模擬結(jié)果;Harafuji 等[6]的分子動(dòng)力學(xué)模擬所使用的勢(shì)函數(shù)為庫(kù)侖項(xiàng)、短程排斥項(xiàng)、范德瓦耳斯項(xiàng)和共價(jià)項(xiàng)結(jié)合的形式,并采用兩相法對(duì) GaN的熔化曲線進(jìn)行了計(jì)算,原子總數(shù)為9408 個(gè),與本文構(gòu)建兩相共存體系的方法不同,Harafuji 等[6]先使GaN 體系在20 ps 內(nèi)產(chǎn)生一個(gè)固液等分結(jié)構(gòu),隨后將體系在不同溫度下弛豫 60 ps,通過(guò)觀察固液界面的移動(dòng)來(lái)確定熔化溫度,相比于Harafuji 等[6]的分子動(dòng)力學(xué)模擬,本文延長(zhǎng)了模擬中兩相共存階段的弛豫時(shí)間,使之達(dá)到了200 ps,更長(zhǎng)的弛豫時(shí)間可以使整個(gè)GaN 體系更加接近平衡態(tài),對(duì)提高模擬結(jié)果的可靠性有利,本文的勢(shì)函數(shù)選擇和兩相法操作過(guò)程與 Harafuji 等[6]不同而導(dǎo)致結(jié)果出現(xiàn)了差異,但熔化曲線的趨勢(shì)是一致的;Van Vechten[2]提出的GaN熔化曲線始終呈下降趨勢(shì),該結(jié)果是通過(guò)參考 Si的熔化溫度隨著壓力升高而減小這一事實(shí),結(jié)合經(jīng)典電負(fù)性理論框架提出的,但之后的研究均發(fā)現(xiàn) GaN的熔化曲線趨勢(shì)與 Si 并不相同.大多數(shù) Ⅲ-Ⅴ 族與 Ⅱ-Ⅵ 族半導(dǎo)體的熔化機(jī)制極為復(fù)雜,在熔化時(shí)常常伴隨著配位數(shù)的改變,在環(huán)境條件下這些半導(dǎo)體均具有與 Si 晶體相似的正四面體鍵合方式,稱為正四面體鍵半導(dǎo)體(Tetrahedrally Bonded Semiconductors,TBSs)[38],這種結(jié)構(gòu)中的每個(gè)原子都與 4 個(gè)最近鄰原子結(jié)合,形成四面體配位網(wǎng)格.大部分 TBSs的熔化溫度會(huì)隨著壓力的增大而降低[2],這是因?yàn)門BSs的熔化機(jī)制由兩個(gè)單獨(dú)的過(guò)程組成: 一個(gè)是溫度的增大導(dǎo)致 TBSs的配位數(shù)從 4 增大到 6,而這個(gè)轉(zhuǎn)變溫度會(huì)隨著壓力的增大而降低;另一個(gè)是在材料的自然熔化過(guò)程中,熔化溫度會(huì)隨著壓力的增大而增大[38].兩個(gè)過(guò)程的結(jié)合導(dǎo)致 Si,Ge 等 TBSs 直接從 4 配位固態(tài)熔化為密度更高的6 配位液態(tài)(短程有序,仍存在部分共價(jià)鍵),使得最終TBSs的熔化溫度會(huì)隨著壓力的增大而降低,但 Harafuji 等[6]與 Porowski等[14]的計(jì)算結(jié)果表明 GaN的熔化溫度會(huì)隨著壓力的增大而升高,這與其他 TBSs的情況并不相同.Porowski 等[38]在 2019 年對(duì)該問(wèn)題進(jìn)行了深入研究,他們對(duì) Drozd-Rzoska 模型[14]進(jìn)行線性外推,得到 GaN 在零壓下 4 到 6 配位的理論轉(zhuǎn)變溫度為6000 K,遠(yuǎn)高于實(shí)際的熔化溫度,因此會(huì)在4 配位固態(tài)下直接熔化,形成低密度的4 配位液態(tài),4 配位固態(tài)的熔化溫度會(huì)隨著壓力的增大而升高,因此出現(xiàn)了與其余 TBSs 不同的熔化情況.壓力的增大使GaN的4 到 6 配位轉(zhuǎn)變溫度降低,當(dāng)轉(zhuǎn)變溫度低于 4 配位固態(tài)的熔化溫度時(shí),GaN的熔化溫度才會(huì)隨著壓力的增大而降低.如果僅從表象來(lái)看,GaN的這種現(xiàn)象與其他 TBSs相比是“異常”的,但實(shí)際上,Si,Ge 等 TBSs 在負(fù)壓力下也可能會(huì)出現(xiàn)這種情況,從本質(zhì)上來(lái)說(shuō),GaN的這種“異?!眱H僅是 TBSs的熔化本質(zhì)在正壓下的體現(xiàn).由于實(shí)驗(yàn)上實(shí)現(xiàn)負(fù)壓非常困難,想要探索 Si,Ge 等 TBSs的熔化機(jī)制就必須克服負(fù)壓?jiǎn)栴},而 GaN的這種在正壓下的“異?!比刍赡苁翘剿?TBSs 熔化機(jī)制的一個(gè)捷徑.

        圖7 GaN 纖鋅礦結(jié)構(gòu)熔化溫度的已有研究結(jié)果與本文計(jì)算結(jié)果的對(duì)比Fig.7.Comparison between existing results on the melting temperature of GaN wurtzite structures and the results calculated in this paper.

        在本文的分子動(dòng)力學(xué)模擬中,GaN 纖鋅礦結(jié)構(gòu)在相變到巖鹽結(jié)構(gòu)之前都沒(méi)有出現(xiàn)熔化溫度的降低,這說(shuō)明 GaN的4 到 6 配位轉(zhuǎn)變溫度在相變前均高于4 配位固態(tài)的熔化溫度.另外,在分子動(dòng)力學(xué)模擬中出現(xiàn)了與 Harafuji 等[6]相似的情況,即超離子相,最初認(rèn)為超離子相可能是 GaN的6 配位液態(tài),但結(jié)果中出現(xiàn)的超離子相仍具有固態(tài)結(jié)構(gòu),同時(shí)在GaN的巖鹽結(jié)構(gòu)中也出現(xiàn)了超離子相,這說(shuō)明 GaN超離子相與 6 配位液態(tài)并不一樣,是一種未被發(fā)現(xiàn)的新相,該相是在 GaN 熔化前出現(xiàn)的一個(gè)高溫高壓相,但目前為止并沒(méi)有相關(guān)的實(shí)驗(yàn)證據(jù)表明 GaN 存在該相,這需要相關(guān)工作者繼續(xù)進(jìn)行研究.

        3.2 GaN 纖鋅礦-巖鹽結(jié)構(gòu)的相界線

        圖8為使用第一性原理計(jì)算得到的GaN 3 種結(jié)構(gòu)(纖鋅礦、閃鋅礦和巖鹽結(jié)構(gòu))相對(duì)焓隨壓力的變化關(guān)系,在 43.6 GPa 之前,3 種結(jié)構(gòu)的穩(wěn)定性順序?yàn)? 纖鋅礦結(jié)構(gòu) > 閃鋅礦結(jié)構(gòu) > 巖鹽結(jié)構(gòu).GaN 閃鋅礦與巖鹽結(jié)構(gòu)的相對(duì)焓相交于43.6 GPa,與 Sun 等[39]的計(jì)算結(jié)果 48 GPa 相近,在這之后巖鹽結(jié)構(gòu)的穩(wěn)定性大于閃鋅礦結(jié)構(gòu),纖鋅礦與巖鹽結(jié)構(gòu)的相對(duì)焓相交于 44.3 GPa,之后巖鹽結(jié)構(gòu)成為最穩(wěn)定的高壓相,在 0—80 GPa 壓力范圍內(nèi)GaN 纖鋅礦與閃鋅礦結(jié)構(gòu)的相對(duì)焓均無(wú)交點(diǎn),這表明在零溫下不會(huì)發(fā)生纖鋅礦到閃鋅礦結(jié)構(gòu)的相變.圖9為使用分子動(dòng)力學(xué)方法對(duì)GaN 纖鋅礦與巖鹽結(jié)構(gòu)固-固相界線的追蹤結(jié)果與已知研究結(jié)果的對(duì)比,目前對(duì) GaN 固-固相變壓力的研究結(jié)果有限且存在很大分歧,Zhou 等[3]使用了第一性原理計(jì)算結(jié)合準(zhǔn)諧近似的方法,得到零溫下GaN 纖鋅礦到巖鹽結(jié)構(gòu)的相變壓力為32.4 GPa,準(zhǔn)諧近似方法彌補(bǔ)了第一性原理計(jì)算無(wú)法得到非零溫度下材料性質(zhì)的缺陷,但由于該方法僅考慮低階的非諧效應(yīng),因而在涉及到更高溫度的相變壓力計(jì)算中難以得到可靠的結(jié)果,這使得 Zhou 等[3]在高溫下得到的相變壓力比已有結(jié)果低,但零溫下的相變壓力與 Xia 等[8]的實(shí)驗(yàn)結(jié)果和 Saitta 等[40]的結(jié)果相近;Sadovyi 等[4]采用第一性原理計(jì)算結(jié)合自洽聲子方法對(duì)GaN 纖鋅礦結(jié)構(gòu)和巖鹽結(jié)構(gòu)的相界線進(jìn)行了計(jì)算,該自洽聲子方法引入了 6 階非諧力,這對(duì)計(jì)算結(jié)果的準(zhǔn)確性有較大提升,Sadovyi 等[4]在零溫下得到的GaN 纖鋅礦到巖鹽結(jié)構(gòu)的相變壓力為45.0 GPa,與已有結(jié)果吻合[11,41];Christensen等[42]與等[43]對(duì)GaN 纖鋅礦和巖鹽結(jié)構(gòu)相變壓力的第一性原理計(jì)算結(jié)果分別為51.8 GPa和56.06 GPa,與Ueno 等[10]的實(shí)驗(yàn)值 52.2 GPa吻合,但略高于本文的計(jì)算結(jié)果.本文在使用分子動(dòng)力學(xué)模擬對(duì) GaN 相變壓力進(jìn)行計(jì)算時(shí),采用了經(jīng)過(guò)驗(yàn)證的有效勢(shì)函數(shù),結(jié)合克勞修斯-克拉珀龍方程,并采取盡可能小的壓力增幅對(duì)GaN的相界線進(jìn)行逐步追蹤計(jì)算,得到較低溫度下纖鋅礦到巖鹽結(jié)構(gòu)的相變壓力為45.9 GPa,這與已有的實(shí)驗(yàn)結(jié)果[41]、計(jì)算結(jié)果[4,11]以及本工作的第一性原理計(jì)算結(jié)果均吻合,但在相同壓力下的相轉(zhuǎn)變溫度比Zhou 等[3]和 Sadovyi 等[4]的結(jié)果高.Sadovyi 等[4]推斷 GaN 巖鹽結(jié)構(gòu)的熔化溫度在 37 GPa 時(shí)應(yīng)該略高于 2100 K,而本文的計(jì)算結(jié)果為4900 K,高的熔化溫度使固-固-液三相共存點(diǎn)位置更高,因此從該點(diǎn)追蹤到的相界線必然會(huì)有高的相轉(zhuǎn)變溫度.綜合來(lái)看,研究中計(jì)算方法的選擇不同是導(dǎo)致GaN 結(jié)構(gòu)相變壓力分歧較大的主要原因.

        圖8 第一性原理計(jì)算得到的GaN 纖鋅礦、閃鋅礦與巖鹽三種結(jié)構(gòu)的相對(duì)焓,內(nèi)插圖為放大的相對(duì)焓交點(diǎn)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 纖鋅礦與巖鹽結(jié)構(gòu)固-固相界線的已有數(shù)據(jù)與本文計(jì)算結(jié)果的對(duì)比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 纖鋅礦結(jié)構(gòu)與巖鹽結(jié)構(gòu)在0—80 GPa 壓力范圍內(nèi)的P-T相圖,圖中給出了GaN的5 種存在狀態(tài): 固態(tài)纖鋅礦結(jié)構(gòu)、纖鋅礦超離子相、固態(tài)巖鹽結(jié)構(gòu)、巖鹽超離子相和液態(tài),其中固態(tài)纖鋅礦結(jié)構(gòu)、固態(tài)巖鹽結(jié)構(gòu)與液態(tài)存在的溫度范圍較寬,兩個(gè)超離子相存在的溫度范圍較窄.通過(guò)外推纖鋅礦結(jié)構(gòu)的熔化曲線,可以得到GaN 在環(huán)境條件下的熔化溫度為2295 K,隨著壓力的增大,GaN 纖鋅礦結(jié)構(gòu)的熔化溫度逐漸升高,并于 P1 點(diǎn)與巖鹽結(jié)構(gòu)的熔化曲線相交,該點(diǎn)為GaN的液態(tài)、纖鋅礦超離子相、巖鹽超離子相的三相共存點(diǎn);從 P1 點(diǎn)開始,GaN 巖鹽結(jié)構(gòu)的熔化溫度隨著壓力的增大而逐漸升高,相比于纖鋅礦結(jié)構(gòu)的熔化曲線,巖鹽結(jié)構(gòu)的熔化曲線更加陡峭,在 80 GPa 時(shí),熔化溫度已經(jīng)達(dá)到了 7646 K;GaN 纖鋅礦結(jié)構(gòu)與其超離子相的相界線為一直線,于 P2 點(diǎn)與纖鋅礦結(jié)構(gòu)的熔化曲線相交,該點(diǎn)為液態(tài)、纖鋅礦超離子相、固態(tài)纖鋅礦結(jié)構(gòu)的三相共存點(diǎn),纖鋅礦結(jié)構(gòu)與超離子相的相界線并沒(méi)有與縱坐標(biāo)相交,說(shuō)明GaN 在常壓下不存在超離子相,該相為一高溫高壓相,只有在 2.0 GPa 之后,且在高溫的驅(qū)動(dòng)下,才進(jìn)入超離子相;GaN 巖鹽結(jié)構(gòu)的超離子相轉(zhuǎn)變溫度非常接近熔化溫度,但隨著壓力的增大而逐漸遠(yuǎn)離熔化溫度,在 80 GPa 下,巖鹽結(jié)構(gòu)超離子相的轉(zhuǎn)變溫度為6162 K.對(duì)GaN 固-固相界線的追蹤從 P1 點(diǎn)開始,從 P1 點(diǎn)得到的克勞修斯-克拉珀龍斜率為331.9 K/GPa,在該斜率方向上壓力溫度條件為32.6 GPa,4321 K 處作為新起點(diǎn),得到第2 個(gè)斜率: —275.8 K/GPa,此時(shí)的克勞修斯-克拉珀龍斜率由正變負(fù),沿該方向的GaN 相界線于 P3 點(diǎn)與巖鹽超離子相的相界線相交,P3 點(diǎn)為纖鋅礦超離子相、巖鹽超離子相與固態(tài)巖鹽結(jié)構(gòu)的三相共存點(diǎn);以 P3 點(diǎn)作為新起點(diǎn)繼續(xù)進(jìn)行固-固相界線的追蹤,在點(diǎn) P4 處GaN的固-固相界線與纖鋅礦超離子相的相界線相交,該點(diǎn)為固態(tài)纖鋅礦結(jié)構(gòu)、纖鋅礦超離子相、固態(tài)巖鹽結(jié)構(gòu)的三相共存點(diǎn),克勞修斯-克拉珀龍斜率為—81.2 K/GPa,為相界線中克勞修斯-克拉珀龍斜率絕對(duì)值的最小值;從 P4 點(diǎn)開始再無(wú)特殊的三相共存點(diǎn),只需一步步進(jìn)行相界線的追蹤,并適當(dāng)調(diào)整每次追蹤的距離,即可得到完整的GaN 固態(tài)纖鋅礦與固態(tài)巖鹽結(jié)構(gòu)的相界線.在本文的計(jì)算結(jié)果中,GaN的相界線不是一條直線,這在固態(tài)纖鋅礦結(jié)構(gòu)與固態(tài)巖鹽結(jié)構(gòu)的相界線中尤為明顯,在壓力大于 40 GPa 后,相界線斜率開始突然增大,并于 45.9 GPa 與x軸相交.相圖中特殊的三相共存點(diǎn)位置見表3.

        圖10 GaN 0—80 GPa 壓力范圍內(nèi)的P-T 相圖,圖中的所有點(diǎn)均為分子動(dòng)力學(xué)模擬值,黑色和紅色線為GaN纖鋅礦結(jié)構(gòu)、巖鹽結(jié)構(gòu)熔化溫度的擬合曲線;海藍(lán)色與紫色線為GaN 纖鋅礦超離子相、巖鹽超離子相的轉(zhuǎn)變邊界擬合曲線;藍(lán)色、綠色與橙色點(diǎn)線分別為GaN 纖鋅礦超離子相-巖鹽超離子相、纖鋅礦超離子相-固態(tài)巖鹽結(jié)構(gòu)、固態(tài)纖鋅礦結(jié)構(gòu)-固態(tài)巖鹽結(jié)構(gòu)的相轉(zhuǎn)變邊界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 相圖中三相共存點(diǎn)的位置Table 3.Location of three-phase coexistence points in GaN P-T phase diagram.

        4 結(jié)論

        本文采用基于 Morse 與 Born-Mayer 組合勢(shì)的分子動(dòng)力學(xué)模擬方法,對(duì) GaN 纖鋅礦與巖鹽結(jié)構(gòu)在寬廣的溫度壓力范圍內(nèi)的P-T相圖進(jìn)行了可靠預(yù)測(cè).首先使用第一性原理計(jì)算方法對(duì) GaN 纖鋅礦結(jié)構(gòu)和巖鹽結(jié)構(gòu)的晶格常數(shù)、彈性常數(shù)、體積模量與晶格能進(jìn)行了計(jì)算,然后與使用基于已選勢(shì)函數(shù)的晶格動(dòng)力學(xué)方法計(jì)算的結(jié)果進(jìn)行了對(duì)比,該勢(shì)函數(shù)能夠給出與第一性原理計(jì)算和實(shí)驗(yàn)結(jié)果相近的力學(xué)、熱力學(xué)性質(zhì)與晶格能,相近的晶格能說(shuō)明該勢(shì)函數(shù)能夠很好地再現(xiàn) GaN 在高溫下尤其是熔化后的性質(zhì).使用分子動(dòng)力學(xué)模擬方法計(jì)算了300 K 下 GaN 纖鋅礦結(jié)構(gòu)與巖鹽結(jié)構(gòu)體積比率隨壓力的變化情況,在纖鋅礦結(jié)構(gòu)下,體積比率隨壓力的變化情況與實(shí)驗(yàn)數(shù)據(jù)符合良好,相變?yōu)閹r鹽結(jié)構(gòu)時(shí)的體積變化率為14.4%;計(jì)算了GaN 纖鋅礦結(jié)構(gòu)與巖鹽結(jié)構(gòu)的熔化曲線,零壓下纖鋅礦結(jié)構(gòu)的熔化溫度為2295 K,未發(fā)現(xiàn)熔化溫度隨壓力增大而降低的情況;對(duì) GaN 纖鋅礦與巖鹽結(jié)構(gòu)的固-固相界線進(jìn)行了追蹤,零溫下纖鋅礦到巖鹽結(jié)構(gòu)的相變壓力為45.9 GPa,與第一性原理計(jì)算結(jié)果44.3 GPa 符合良好;對(duì) GaN 可能存在的超離子相進(jìn)行了確定,在該狀態(tài)下,GaN的Ga 原子亞晶格局部熔化,Ga 原子能夠在晶體內(nèi)部自由移動(dòng)并起導(dǎo)電作用,纖鋅礦結(jié)構(gòu)在壓力大于2.0 GPa 且溫度高于 2550 K 時(shí)進(jìn)入超離子相,而巖鹽結(jié)構(gòu)在壓力溫度大于 33.1 GPa 和 4182 K 后發(fā)生超離子相轉(zhuǎn)變.

        猜你喜歡
        鋅礦巖鹽勢(shì)函數(shù)
        航天器姿態(tài)受限的協(xié)同勢(shì)函數(shù)族設(shè)計(jì)方法
        鈣(鎂)離子在菱鋅礦表面吸附的量子化學(xué)研究
        數(shù)學(xué)理論與應(yīng)用(2022年1期)2022-04-15 09:03:32
        金屬鎢級(jí)聯(lián)碰撞中勢(shì)函數(shù)的影響
        青海北祁連陰凹槽塞浦路斯型銅鋅礦特征及找礦標(biāo)志
        澳大利亞杜加爾河鋅礦實(shí)現(xiàn)商業(yè)化生產(chǎn)
        SOME RESULTS OF WEAKLY f-STATIONARY MAPS WITH POTENTIAL
        大蒜
        中國(guó)東部巖鹽礦區(qū)建造鹽穴儲(chǔ)氣庫(kù)地質(zhì)條件分析
        察爾汗鹽湖巖鹽路基應(yīng)力作用下溶蝕特性試驗(yàn)研究
        一区二区三区放荡人妻| 97久久精品人妻人人搡人人玩| 人妻激情另类乱人伦人妻| 青草福利在线| 亚洲色图在线视频免费观看| 国产精品一区二区偷拍| 久久婷婷五月综合97色一本一本| 4444亚洲人成无码网在线观看 | 国产精品国产三级国产an不卡| 久久久久亚洲av成人人电影| 国产av精国产传媒| 99久久国产亚洲综合精品| 青青草好吊色在线视频| 国内永久福利在线视频图片| 午夜成人无码福利免费视频 | 一本久久综合亚洲鲁鲁五月夫| 美女人妻中出日本人妻| 又色又爽又黄还免费毛片96下载| 国产人成午夜免电影观看| 国产后入内射在线观看| 亚洲精品av一区二区| 天天燥日日燥| 亚洲丁香五月激情综合| 亚洲天堂线上免费av| 日本边添边摸边做边爱喷水| 久久午夜无码鲁丝片直播午夜精品 | 成年人视频在线播放麻豆| 一个色综合中文字幕人妻激情视频 | 毛片一级精油按摩无码| 亚洲黄色av一区二区三区| 性xxxx18免费观看视频| 亚洲中文无码久久精品1| 亚洲国产精品二区三区| 午夜久久久久久禁播电影| 无码国产色欲xxxxx视频| 蜜臀aⅴ永久无码一区二区| 国产女同舌吻1区2区| 国产精品一区二区久久乐下载 | 日本av亚洲中文字幕| 中国老熟妇自拍hd发布| 精品一区二区三区久久久|