• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

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

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

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

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

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

    1 引言

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

    理論和實驗研究表明,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)可通過氣相外延和分子束外延沉積在立方(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)相近,實驗上還未觀察到該相變,但也有研究表明該相變發(fā)生在 10.87 GPa[11].GaN的高溫高壓實驗條件苛刻,而可以計算高溫相變的分子模擬又缺乏準確的原子力場,因此,對 GaN 進行高溫高壓相圖的研究是一項困難的工作.Van Vechten[2]提出的化學(xué)鍵模型預(yù)測了 GaN 在高溫下的固-固相變壓力,但該模型預(yù)測的是纖鋅礦結(jié)構(gòu)與β-tin 結(jié)構(gòu)的相變.關(guān)于 GaN 纖鋅礦到巖鹽結(jié)構(gòu)的相界線研究較少,僅有 Zhou 等[3]的數(shù)值模擬和 Sadovyi 等[4]的高溫高壓實驗結(jié)果,后者對GaN的固-固相變行為進行了較為全面的分析.

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

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

    2 模型及方法

    2.1 GaN 勢函數(shù)及其驗證

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

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

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

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

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

    圖2 GaN 纖鋅礦和巖鹽結(jié)構(gòu)在 300 K 下的體積比率與已有結(jié)果對比,所有數(shù)據(jù)均歸一化至纖鋅礦結(jié)構(gòu)初始體積,紅色與黑色圓點為分子動力學(xué)模擬結(jié)果(實線為擬合曲線);藍色正三角形與綠色菱形分別為Xia 等[8]與 Ueno 等[10]的X 射線衍射實驗結(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的正交體系進行 GaN 纖鋅礦與巖鹽結(jié)構(gòu)的分子動力學(xué)模擬.兩個體系均包含6144個原子,采用周期性邊界條件,時間步長1 fs,所有的分子動力學(xué)計算均利用 LAMMPS 軟件[32]完成.目前分子動力學(xué)模擬中計算熔化溫度的方法主要有單相法、兩相法、孔洞法、Z 方法以及改進的Z 方法[33-36],本文選擇兩相法進行 GaN的熔化溫度計算,該方法準確度很高,而且也能在一定程度上控制壓力,在進行熔化曲線的計算時,具有較好的操作空間.兩相法計算熔化溫度的關(guān)鍵在于構(gòu)建能夠長時間存在的兩相共存體系,為達到這個目的,首先將整個模擬體系在等溫等壓系綜(NPT)中以略低于熔點的溫度弛豫30 ps,保證整個體系的應(yīng)力平衡,然后將體系沿z方向平均分為兩部分,一部分設(shè)置受力為0,將其固定,另一部分在遠高于熔點的溫度下弛豫 50 ps,保證其完全熔化,再將完全熔化的這部分原子降溫到略高于熔點的溫度,降溫時間為10 ps,最后釋放固定的原子,將整個體系置于微正則系綜(NVE)中自由弛豫,形成固液共存體系,若整個體系能夠以固液共存狀態(tài)長時間存在(本文中為200 ps),且固液部分體積無顯著變化,則可以認為此時體系的溫度為最終熔化溫度.圖3為構(gòu)建的GaN 兩相共存體系,圖4為對應(yīng)的原子數(shù)密度,可以清楚地看出體系明顯分為兩部分,數(shù)密度波動大的部分處于固態(tài),數(shù)密度波動較小的部分處于液態(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的超離子相

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

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

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

    圖5 壓力為10 GPa 時,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 時的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)在相平衡時壓力隨溫度的變化率,可以用于確定兩相的相界線斜率,適用于本工作中確定 GaN的固-固相界線斜率,其定義為

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

    3 結(jié)果與討論

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

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

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

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

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

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

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

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

    4 結(jié)論

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

    猜你喜歡
    鋅礦巖鹽勢函數(shù)
    航天器姿態(tài)受限的協(xié)同勢函數(shù)族設(shè)計方法
    鈣(鎂)離子在菱鋅礦表面吸附的量子化學(xué)研究
    次可加勢函數(shù)拓撲壓及因子映射
    金屬鎢級聯(lián)碰撞中勢函數(shù)的影響
    青海北祁連陰凹槽塞浦路斯型銅鋅礦特征及找礦標(biāo)志
    澳大利亞杜加爾河鋅礦實現(xiàn)商業(yè)化生產(chǎn)
    SOME RESULTS OF WEAKLY f-STATIONARY MAPS WITH POTENTIAL
    大蒜
    中國東部巖鹽礦區(qū)建造鹽穴儲氣庫地質(zhì)條件分析
    察爾汗鹽湖巖鹽路基應(yīng)力作用下溶蝕特性試驗研究
    美女大奶头黄色视频| 高清在线视频一区二区三区| 亚洲精品国产av蜜桃| 一个人免费看片子| 欧美黄色片欧美黄色片| 美女高潮到喷水免费观看| 国产一区二区三区综合在线观看| 日本黄色日本黄色录像| 成年人免费黄色播放视频| 欧美+日韩+精品| 亚洲,欧美精品.| 黄色一级大片看看| 深夜精品福利| 亚洲久久久国产精品| 人妻 亚洲 视频| a 毛片基地| 一本大道久久a久久精品| 日韩欧美精品免费久久| 老司机亚洲免费影院| 日本欧美国产在线视频| 精品国产超薄肉色丝袜足j| xxx大片免费视频| 少妇人妻 视频| 日韩av在线免费看完整版不卡| 91aial.com中文字幕在线观看| 色婷婷av一区二区三区视频| 亚洲,一卡二卡三卡| 亚洲精品在线美女| 久久亚洲国产成人精品v| 91午夜精品亚洲一区二区三区| 一边亲一边摸免费视频| 国产av国产精品国产| 天堂俺去俺来也www色官网| 亚洲人成网站在线观看播放| 久久99一区二区三区| 中国三级夫妇交换| 性色avwww在线观看| 伊人久久国产一区二区| 亚洲成色77777| 又大又黄又爽视频免费| 黄色毛片三级朝国网站| 飞空精品影院首页| 国产一区二区三区综合在线观看| 国产精品一国产av| 久久热在线av| 另类精品久久| 亚洲成av片中文字幕在线观看 | 亚洲精品乱久久久久久| 性色avwww在线观看| 精品亚洲成a人片在线观看| 国产高清不卡午夜福利| 日韩不卡一区二区三区视频在线| 午夜免费鲁丝| 黄色 视频免费看| 欧美人与善性xxx| 午夜日本视频在线| av天堂久久9| 麻豆精品久久久久久蜜桃| 午夜免费男女啪啪视频观看| 1024视频免费在线观看| 男女国产视频网站| 欧美日韩av久久| 午夜免费观看性视频| 亚洲一码二码三码区别大吗| 高清av免费在线| 下体分泌物呈黄色| 秋霞伦理黄片| a级毛片在线看网站| 国产精品 欧美亚洲| 亚洲欧美精品综合一区二区三区 | 日本vs欧美在线观看视频| 国产精品久久久久久av不卡| 九九爱精品视频在线观看| 看十八女毛片水多多多| 久久这里只有精品19| 制服丝袜香蕉在线| 久久久久久久亚洲中文字幕| 人妻 亚洲 视频| 亚洲精品视频女| 色哟哟·www| 亚洲成人手机| av免费观看日本| 女性生殖器流出的白浆| 国产成人午夜福利电影在线观看| 免费大片黄手机在线观看| 最近最新中文字幕大全免费视频 | 久久久久久久久久人人人人人人| 日本vs欧美在线观看视频| 一级毛片我不卡| 久久久国产欧美日韩av| 久久人妻熟女aⅴ| 国产成人精品婷婷| 黄色视频在线播放观看不卡| 另类亚洲欧美激情| 汤姆久久久久久久影院中文字幕| 我的亚洲天堂| 男人舔女人的私密视频| 美女脱内裤让男人舔精品视频| 少妇的丰满在线观看| 国产精品熟女久久久久浪| 在线观看免费视频网站a站| 午夜日本视频在线| 国产综合精华液| 欧美日韩亚洲高清精品| 成人免费观看视频高清| 七月丁香在线播放| 成年人午夜在线观看视频| 日韩欧美一区视频在线观看| 丝袜人妻中文字幕| 人人妻人人添人人爽欧美一区卜| 狠狠精品人妻久久久久久综合| 久久鲁丝午夜福利片| 五月伊人婷婷丁香| 好男人视频免费观看在线| 欧美av亚洲av综合av国产av | 一级片'在线观看视频| 日韩 亚洲 欧美在线| 蜜桃在线观看..| 男人操女人黄网站| 热re99久久国产66热| av在线老鸭窝| 男女边摸边吃奶| 久久综合国产亚洲精品| 久久久久久人人人人人| 老司机影院毛片| 亚洲国产欧美在线一区| 少妇人妻精品综合一区二区| 亚洲久久久国产精品| 亚洲欧洲国产日韩| 最近手机中文字幕大全| 成年女人在线观看亚洲视频| 在线观看免费日韩欧美大片| 日本爱情动作片www.在线观看| 国产深夜福利视频在线观看| 亚洲国产精品一区二区三区在线| 久久久精品94久久精品| 国产精品久久久久久精品古装| 五月开心婷婷网| 国产精品亚洲av一区麻豆 | 亚洲视频免费观看视频| 日本午夜av视频| 亚洲国产最新在线播放| 久久狼人影院| 只有这里有精品99| 黄色毛片三级朝国网站| 在线观看一区二区三区激情| 亚洲一区二区三区欧美精品| 久久青草综合色| 欧美成人午夜免费资源| 寂寞人妻少妇视频99o| 中文字幕最新亚洲高清| 晚上一个人看的免费电影| 久久国产精品大桥未久av| xxx大片免费视频| 人妻系列 视频| 日韩,欧美,国产一区二区三区| 亚洲国产av影院在线观看| 老司机影院成人| 熟女av电影| 国产国语露脸激情在线看| 午夜福利影视在线免费观看| 日韩欧美一区视频在线观看| 观看av在线不卡| 黄色怎么调成土黄色| 韩国av在线不卡| 女性生殖器流出的白浆| 熟女少妇亚洲综合色aaa.| www日本在线高清视频| 久久久a久久爽久久v久久| 日本免费在线观看一区| 色吧在线观看| 欧美国产精品va在线观看不卡| 香蕉精品网在线| 人妻一区二区av| 大香蕉久久成人网| 中文字幕最新亚洲高清| 欧美 亚洲 国产 日韩一| 男人舔女人的私密视频| 老汉色av国产亚洲站长工具| 免费观看a级毛片全部| 91精品国产国语对白视频| 在线亚洲精品国产二区图片欧美| 国产成人精品久久久久久| 中文字幕最新亚洲高清| av线在线观看网站| 日韩一本色道免费dvd| 国精品久久久久久国模美| 亚洲少妇的诱惑av| 18在线观看网站| 亚洲婷婷狠狠爱综合网| 80岁老熟妇乱子伦牲交| 最新中文字幕久久久久| 曰老女人黄片| 黄色配什么色好看| 一级片'在线观看视频| 国产一级毛片在线| 精品国产乱码久久久久久男人| 亚洲伊人久久精品综合| 天天躁日日躁夜夜躁夜夜| 婷婷色综合www| 久久久久久久久免费视频了| 99久国产av精品国产电影| 久久精品国产亚洲av天美| 久久精品国产亚洲av涩爱| 在线观看一区二区三区激情| 亚洲国产日韩一区二区| 亚洲伊人久久精品综合| 一级毛片我不卡| 国产亚洲av片在线观看秒播厂| 久久久久精品人妻al黑| 国产爽快片一区二区三区| 香蕉精品网在线| 一区二区日韩欧美中文字幕| 久久久国产精品麻豆| 一区在线观看完整版| av国产精品久久久久影院| 在线观看www视频免费| 国产伦理片在线播放av一区| 桃花免费在线播放| 国产精品一二三区在线看| 欧美激情高清一区二区三区 | 啦啦啦中文免费视频观看日本| 一级爰片在线观看| 日本色播在线视频| 日韩欧美精品免费久久| 免费在线观看完整版高清| av免费在线看不卡| 亚洲成人一二三区av| 免费高清在线观看日韩| 十分钟在线观看高清视频www| av不卡在线播放| 日韩电影二区| 中文字幕人妻熟女乱码| a 毛片基地| 亚洲成色77777| 久久久久人妻精品一区果冻| 999精品在线视频| 看免费成人av毛片| av网站免费在线观看视频| 成年人免费黄色播放视频| 午夜av观看不卡| 丰满迷人的少妇在线观看| 亚洲一区中文字幕在线| 亚洲,欧美,日韩| 精品久久蜜臀av无| 国产欧美日韩一区二区三区在线| a级毛片在线看网站| videossex国产| 午夜福利,免费看| 国产精品人妻久久久影院| 女的被弄到高潮叫床怎么办| 母亲3免费完整高清在线观看 | 伊人久久国产一区二区| 午夜av观看不卡| 午夜激情av网站| 成人亚洲欧美一区二区av| 欧美日韩视频高清一区二区三区二| 日韩视频在线欧美| 一区二区av电影网| 超碰97精品在线观看| 亚洲精品aⅴ在线观看| 91精品伊人久久大香线蕉| 热99久久久久精品小说推荐| 国产一区亚洲一区在线观看| 久久久久久人人人人人| 国产极品天堂在线| 亚洲一区中文字幕在线| 欧美最新免费一区二区三区| 久久人人爽av亚洲精品天堂| 成年动漫av网址| 性少妇av在线| 亚洲国产欧美日韩在线播放| 天天躁夜夜躁狠狠久久av| 丰满少妇做爰视频| 黄色配什么色好看| 国产一区二区 视频在线| 精品视频人人做人人爽| 色吧在线观看| 三上悠亚av全集在线观看| 国产精品国产三级专区第一集| 老司机亚洲免费影院| 人人妻人人澡人人爽人人夜夜| 两个人看的免费小视频| 久久精品aⅴ一区二区三区四区 | 飞空精品影院首页| 亚洲人成电影观看| 亚洲欧美中文字幕日韩二区| 夫妻性生交免费视频一级片| 黄色怎么调成土黄色| 久久毛片免费看一区二区三区| 成年女人毛片免费观看观看9 | 免费在线观看视频国产中文字幕亚洲 | 午夜日本视频在线| a级毛片黄视频| 99国产综合亚洲精品| 啦啦啦在线免费观看视频4| 欧美老熟妇乱子伦牲交| 国产精品 欧美亚洲| tube8黄色片| 国产麻豆69| 蜜桃国产av成人99| 在线观看免费日韩欧美大片| 一二三四在线观看免费中文在| 成年美女黄网站色视频大全免费| 好男人视频免费观看在线| 亚洲国产日韩一区二区| 国产成人av激情在线播放| 少妇被粗大的猛进出69影院| 国产激情久久老熟女| tube8黄色片| 亚洲伊人久久精品综合| 熟妇人妻不卡中文字幕| 美女主播在线视频| 婷婷色av中文字幕| 免费人妻精品一区二区三区视频| 97精品久久久久久久久久精品| 永久网站在线| 夜夜骑夜夜射夜夜干| 看十八女毛片水多多多| 国精品久久久久久国模美| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 国产欧美日韩一区二区三区在线| 制服诱惑二区| 一级黄片播放器| 美女大奶头黄色视频| www日本在线高清视频| 精品国产乱码久久久久久男人| 嫩草影院入口| videosex国产| 久久 成人 亚洲| 久久精品国产亚洲av涩爱| 国产一区二区 视频在线| 国产又色又爽无遮挡免| 欧美国产精品一级二级三级| 日韩电影二区| 久久久国产一区二区| 久久久久久久亚洲中文字幕| 国产97色在线日韩免费| 亚洲少妇的诱惑av| 日韩不卡一区二区三区视频在线| av女优亚洲男人天堂| 妹子高潮喷水视频| 国产欧美日韩一区二区三区在线| 岛国毛片在线播放| 精品一区二区免费观看| 高清欧美精品videossex| 久久毛片免费看一区二区三区| 国产熟女午夜一区二区三区| 啦啦啦啦在线视频资源| 久久精品亚洲av国产电影网| 亚洲欧美日韩另类电影网站| 男的添女的下面高潮视频| 中文乱码字字幕精品一区二区三区| 中文字幕人妻丝袜制服| 国产精品久久久久成人av| 精品国产一区二区三区久久久樱花| 精品卡一卡二卡四卡免费| 久久精品人人爽人人爽视色| 夜夜骑夜夜射夜夜干| 国产激情久久老熟女| 9热在线视频观看99| 捣出白浆h1v1| 最近中文字幕高清免费大全6| 国产精品二区激情视频| 美女xxoo啪啪120秒动态图| 午夜免费观看性视频| 国产男女内射视频| 日韩视频在线欧美| 日韩av免费高清视频| 国产淫语在线视频| 欧美中文综合在线视频| 亚洲人成网站在线观看播放| 国产片内射在线| 亚洲人成网站在线观看播放| av一本久久久久| 18+在线观看网站| 国产精品一区二区在线不卡| 国产成人午夜福利电影在线观看| 一本—道久久a久久精品蜜桃钙片| 91成人精品电影| 色婷婷av一区二区三区视频| 新久久久久国产一级毛片| 99久久中文字幕三级久久日本| av线在线观看网站| 国产成人免费观看mmmm| 亚洲欧美色中文字幕在线| 久久久久视频综合| 精品国产一区二区三区四区第35| 啦啦啦在线观看免费高清www| 巨乳人妻的诱惑在线观看| 又粗又硬又长又爽又黄的视频| 国产xxxxx性猛交| 精品久久蜜臀av无| 一区福利在线观看| 97在线视频观看| 少妇被粗大的猛进出69影院| 成人影院久久| 国产伦理片在线播放av一区| 色视频在线一区二区三区| 不卡视频在线观看欧美| 卡戴珊不雅视频在线播放| av有码第一页| 9191精品国产免费久久| 欧美激情高清一区二区三区 | 久久久久久人人人人人| 日韩制服骚丝袜av| 久久久a久久爽久久v久久| 久久久久国产一级毛片高清牌| 亚洲欧美色中文字幕在线| 一级片'在线观看视频| 精品人妻在线不人妻| 国产高清国产精品国产三级| 国产成人一区二区在线| 久久国产亚洲av麻豆专区| 国产男女超爽视频在线观看| 一本大道久久a久久精品| 国产福利在线免费观看视频| 亚洲欧美一区二区三区国产| 超色免费av| 亚洲精品成人av观看孕妇| 91国产中文字幕| 妹子高潮喷水视频| 人人妻人人澡人人看| 国产在线一区二区三区精| 久久久久视频综合| 80岁老熟妇乱子伦牲交| 色视频在线一区二区三区| 女人久久www免费人成看片| 丰满少妇做爰视频| 人成视频在线观看免费观看| 精品少妇一区二区三区视频日本电影 | 日本免费在线观看一区| 亚洲欧美中文字幕日韩二区| 这个男人来自地球电影免费观看 | 九草在线视频观看| 一区二区三区乱码不卡18| 亚洲欧美一区二区三区黑人 | 免费播放大片免费观看视频在线观看| 日本91视频免费播放| 777米奇影视久久| 考比视频在线观看| 久久99热这里只频精品6学生| 久久精品人人爽人人爽视色| 中文字幕精品免费在线观看视频| 人成视频在线观看免费观看| 麻豆乱淫一区二区| 观看av在线不卡| 亚洲少妇的诱惑av| 国产激情久久老熟女| 亚洲综合精品二区| av网站免费在线观看视频| 国产亚洲一区二区精品| 日韩制服骚丝袜av| h视频一区二区三区| av天堂久久9| 岛国毛片在线播放| 国产精品熟女久久久久浪| 这个男人来自地球电影免费观看 | 女性生殖器流出的白浆| www日本在线高清视频| 这个男人来自地球电影免费观看 | 亚洲经典国产精华液单| 2021少妇久久久久久久久久久| 老鸭窝网址在线观看| 日韩中文字幕欧美一区二区 | 日本色播在线视频| 午夜福利视频在线观看免费| 精品一区在线观看国产| www.av在线官网国产| av电影中文网址| 日韩一本色道免费dvd| 婷婷成人精品国产| 一本—道久久a久久精品蜜桃钙片| 两个人免费观看高清视频| 久久99热这里只频精品6学生| 最近中文字幕2019免费版| 午夜激情av网站| 日本色播在线视频| 一级爰片在线观看| 黑人猛操日本美女一级片| 青春草视频在线免费观看| 啦啦啦在线观看免费高清www| 亚洲,欧美精品.| 日本av手机在线免费观看| 国产成人一区二区在线| 人人妻人人添人人爽欧美一区卜| 叶爱在线成人免费视频播放| 少妇熟女欧美另类| 精品国产一区二区久久| 欧美在线黄色| 精品少妇久久久久久888优播| 熟女电影av网| 久久久久国产精品人妻一区二区| 视频区图区小说| 精品视频人人做人人爽| 亚洲色图综合在线观看| 午夜福利在线观看免费完整高清在| 亚洲精品成人av观看孕妇| 国产国语露脸激情在线看| 青春草视频在线免费观看| 观看美女的网站| 熟女电影av网| 一区二区三区激情视频| 亚洲av.av天堂| 18+在线观看网站| 精品少妇内射三级| 一级,二级,三级黄色视频| 精品99又大又爽又粗少妇毛片| 在线观看三级黄色| 伦精品一区二区三区| 亚洲精品国产av成人精品| 男女边摸边吃奶| 韩国精品一区二区三区| 97在线视频观看| 女人久久www免费人成看片| av女优亚洲男人天堂| 69精品国产乱码久久久| 天天影视国产精品| 精品少妇久久久久久888优播| 天堂中文最新版在线下载| 另类精品久久| 大码成人一级视频| 亚洲精品美女久久av网站| 免费日韩欧美在线观看| 观看av在线不卡| 欧美97在线视频| 精品少妇久久久久久888优播| 午夜福利在线观看免费完整高清在| 妹子高潮喷水视频| 在线观看三级黄色| 免费大片黄手机在线观看| 男人爽女人下面视频在线观看| 波多野结衣一区麻豆| 欧美精品高潮呻吟av久久| 国产综合精华液| 亚洲欧美精品自产自拍| 国产精品久久久久久久久免| 国产亚洲欧美精品永久| 美女高潮到喷水免费观看| 欧美日韩精品网址| 亚洲美女视频黄频| 亚洲视频免费观看视频| av电影中文网址| 两性夫妻黄色片| 国产成人精品一,二区| 人人妻人人添人人爽欧美一区卜| 男女免费视频国产| 国产欧美日韩一区二区三区在线| 午夜久久久在线观看| videosex国产| 久久人人爽人人片av| 日日撸夜夜添| 国产精品二区激情视频| 国产片特级美女逼逼视频| 中文字幕人妻熟女乱码| 久热这里只有精品99| 卡戴珊不雅视频在线播放| 麻豆av在线久日| 男女午夜视频在线观看| 综合色丁香网| 纵有疾风起免费观看全集完整版| 亚洲美女搞黄在线观看| 伦理电影免费视频| 一级毛片黄色毛片免费观看视频| 久久人人97超碰香蕉20202| 精品国产一区二区久久| av在线app专区| 中国三级夫妇交换| 国产精品麻豆人妻色哟哟久久| 97在线视频观看| 亚洲成人av在线免费| 久久久精品区二区三区| 自拍欧美九色日韩亚洲蝌蚪91| 久久久久国产网址| 亚洲欧美精品综合一区二区三区 | 香蕉精品网在线| 日本猛色少妇xxxxx猛交久久| 人人妻人人澡人人爽人人夜夜| 国产乱来视频区| xxx大片免费视频| 五月开心婷婷网| 一边摸一边做爽爽视频免费| 亚洲成国产人片在线观看| 久久久久精品人妻al黑| 亚洲av免费高清在线观看| 美女主播在线视频| 亚洲人成网站在线观看播放| 婷婷色综合www| 看非洲黑人一级黄片| 少妇的逼水好多| 最近手机中文字幕大全| 五月伊人婷婷丁香| 亚洲人成网站在线观看播放| 一本大道久久a久久精品| 少妇被粗大猛烈的视频| 国产成人午夜福利电影在线观看| 免费在线观看视频国产中文字幕亚洲 | 欧美日韩精品网址| 亚洲av男天堂| 十八禁网站网址无遮挡| 最近2019中文字幕mv第一页| 成人国产av品久久久| 国产av码专区亚洲av| 啦啦啦啦在线视频资源| 国产激情久久老熟女| 最近中文字幕2019免费版| 亚洲国产欧美日韩在线播放| 激情五月婷婷亚洲| 永久网站在线| 黑人巨大精品欧美一区二区蜜桃| 久久热在线av| 中文天堂在线官网| 精品亚洲成国产av| 久久 成人 亚洲| 成人手机av|