• <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)力作用下溶蝕特性試驗研究
    日韩三级伦理在线观看| 国产一区亚洲一区在线观看| 国产精品一区二区在线观看99| 国产色婷婷99| 蜜桃国产av成人99| 高清在线视频一区二区三区| 天天躁狠狠躁夜夜躁狠狠躁| 成年动漫av网址| 国产又色又爽无遮挡免| 少妇精品久久久久久久| 亚洲国产毛片av蜜桃av| 韩国精品一区二区三区| av线在线观看网站| 啦啦啦啦在线视频资源| 精品99又大又爽又粗少妇毛片| 激情视频va一区二区三区| 国产精品久久久久久久久免| 三上悠亚av全集在线观看| 狠狠精品人妻久久久久久综合| 秋霞在线观看毛片| 亚洲精品成人av观看孕妇| 97在线视频观看| 欧美少妇被猛烈插入视频| 欧美日韩视频高清一区二区三区二| av一本久久久久| 久久99热这里只频精品6学生| 一边亲一边摸免费视频| 久久久国产一区二区| 欧美精品国产亚洲| 日韩中文字幕视频在线看片| 成人影院久久| 亚洲国产最新在线播放| 国产精品熟女久久久久浪| 亚洲欧美日韩另类电影网站| 麻豆精品久久久久久蜜桃| 黄色 视频免费看| 男女午夜视频在线观看| 国产精品蜜桃在线观看| 久久精品国产a三级三级三级| 午夜精品国产一区二区电影| 乱人伦中国视频| 欧美成人午夜免费资源| 国产乱人偷精品视频| 欧美精品高潮呻吟av久久| 少妇被粗大的猛进出69影院| 国产成人午夜福利电影在线观看| 少妇被粗大的猛进出69影院| 久久毛片免费看一区二区三区| 日韩视频在线欧美| 91成人精品电影| 国产av一区二区精品久久| 亚洲欧美成人精品一区二区| 欧美激情极品国产一区二区三区| 久久这里有精品视频免费| 欧美激情 高清一区二区三区| 午夜福利,免费看| 精品人妻熟女毛片av久久网站| 丝袜喷水一区| 波多野结衣av一区二区av| 久久婷婷青草| 18+在线观看网站| 成人午夜精彩视频在线观看| 一区福利在线观看| h视频一区二区三区| 一区二区日韩欧美中文字幕| 考比视频在线观看| 在线观看免费视频网站a站| av视频免费观看在线观看| 熟女av电影| 高清视频免费观看一区二区| 黄色一级大片看看| 亚洲国产日韩一区二区| 看十八女毛片水多多多| 成年美女黄网站色视频大全免费| 久久精品国产亚洲av高清一级| 中文字幕另类日韩欧美亚洲嫩草| 色网站视频免费| 侵犯人妻中文字幕一二三四区| 另类精品久久| 波多野结衣av一区二区av| 国产成人精品在线电影| 汤姆久久久久久久影院中文字幕| 中国三级夫妇交换| 亚洲色图综合在线观看| 侵犯人妻中文字幕一二三四区| 中文字幕制服av| 秋霞伦理黄片| 欧美bdsm另类| 免费av中文字幕在线| 中文字幕最新亚洲高清| 青春草亚洲视频在线观看| 在线观看免费高清a一片| 黄色视频在线播放观看不卡| 女的被弄到高潮叫床怎么办| 欧美精品高潮呻吟av久久| 午夜福利在线免费观看网站| 精品福利永久在线观看| 电影成人av| 国产成人一区二区在线| 美女脱内裤让男人舔精品视频| 国产精品.久久久| 日本-黄色视频高清免费观看| 亚洲色图综合在线观看| 久久久久视频综合| 久久久久久久久免费视频了| 在线 av 中文字幕| 香蕉精品网在线| 午夜免费男女啪啪视频观看| av有码第一页| 亚洲人成电影观看| 精品一区二区三区四区五区乱码 | 性色av一级| 天堂俺去俺来也www色官网| 天天躁夜夜躁狠狠躁躁| 免费少妇av软件| 免费看不卡的av| 看十八女毛片水多多多| 边亲边吃奶的免费视频| 中文字幕人妻熟女乱码| xxx大片免费视频| 亚洲婷婷狠狠爱综合网| 9热在线视频观看99| 国产日韩欧美视频二区| 国产亚洲精品第一综合不卡| 18+在线观看网站| 在线天堂最新版资源| 男人操女人黄网站| 国产精品亚洲av一区麻豆 | 99久久中文字幕三级久久日本| 免费不卡的大黄色大毛片视频在线观看| 制服人妻中文乱码| 日本欧美国产在线视频| www.自偷自拍.com| 国产一区有黄有色的免费视频| 日韩精品免费视频一区二区三区| 免费在线观看完整版高清| 国产av码专区亚洲av| 亚洲精品,欧美精品| 人人妻人人澡人人看| 国产男女内射视频| 男女边吃奶边做爰视频| 最近中文字幕高清免费大全6| 一本大道久久a久久精品| 国产成人欧美| 国产精品国产av在线观看| 叶爱在线成人免费视频播放| 成人手机av| 大片电影免费在线观看免费| 中文字幕av电影在线播放| 水蜜桃什么品种好| 女人被躁到高潮嗷嗷叫费观| 久久婷婷青草| 精品国产一区二区三区久久久樱花| freevideosex欧美| 2018国产大陆天天弄谢| 国产1区2区3区精品| 少妇人妻精品综合一区二区| 精品亚洲成国产av| 久久国产精品大桥未久av| 国产精品一国产av| 大片电影免费在线观看免费| 王馨瑶露胸无遮挡在线观看| 免费看av在线观看网站| 免费黄频网站在线观看国产| 亚洲男人天堂网一区| 最新的欧美精品一区二区| 可以免费在线观看a视频的电影网站 | 狂野欧美激情性bbbbbb| 免费少妇av软件| 色94色欧美一区二区| 国产黄频视频在线观看| 男女高潮啪啪啪动态图| 热99国产精品久久久久久7| 高清在线视频一区二区三区| 国产一区二区激情短视频 | 国产精品成人在线| 青春草亚洲视频在线观看| 亚洲五月色婷婷综合| 国产色婷婷99| av.在线天堂| 久久99蜜桃精品久久| 国产亚洲午夜精品一区二区久久| 国产1区2区3区精品| 精品国产超薄肉色丝袜足j| 久久精品aⅴ一区二区三区四区 | 成人黄色视频免费在线看| 2018国产大陆天天弄谢| 欧美国产精品va在线观看不卡| 亚洲国产成人一精品久久久| 欧美变态另类bdsm刘玥| 丝袜美腿诱惑在线| 成人手机av| 亚洲欧洲精品一区二区精品久久久 | av视频免费观看在线观看| 久久久久久久久久人人人人人人| 纵有疾风起免费观看全集完整版| 日本av免费视频播放| 99久久中文字幕三级久久日本| 国产精品一国产av| 搡女人真爽免费视频火全软件| 街头女战士在线观看网站| 在线观看三级黄色| 亚洲国产欧美网| 中文字幕制服av| 国产精品不卡视频一区二区| 国产av精品麻豆| 免费观看a级毛片全部| 激情五月婷婷亚洲| 亚洲一区中文字幕在线| 男女午夜视频在线观看| 国产成人午夜福利电影在线观看| 免费高清在线观看日韩| 成人国产麻豆网| 亚洲欧洲精品一区二区精品久久久 | 成年人免费黄色播放视频| 欧美精品人与动牲交sv欧美| tube8黄色片| 人妻 亚洲 视频| 高清黄色对白视频在线免费看| 亚洲精品,欧美精品| 久久ye,这里只有精品| 亚洲欧美一区二区三区国产| 国产精品久久久av美女十八| 熟女电影av网| 国产欧美日韩综合在线一区二区| 国产日韩一区二区三区精品不卡| a 毛片基地| 美女午夜性视频免费| 91午夜精品亚洲一区二区三区| 菩萨蛮人人尽说江南好唐韦庄| 夫妻午夜视频| 免费日韩欧美在线观看| 九九爱精品视频在线观看| 免费大片黄手机在线观看| 中文字幕人妻熟女乱码| 卡戴珊不雅视频在线播放| 国产一级毛片在线| 黄色配什么色好看| 两性夫妻黄色片| 啦啦啦视频在线资源免费观看| 另类精品久久| 精品酒店卫生间| freevideosex欧美| 国产片内射在线| 人妻少妇偷人精品九色| 日本欧美视频一区| 午夜福利在线观看免费完整高清在| 国产欧美日韩综合在线一区二区| 观看美女的网站| 老汉色∧v一级毛片| 国产在线视频一区二区| 亚洲美女视频黄频| 亚洲综合精品二区| 少妇精品久久久久久久| 大陆偷拍与自拍| 菩萨蛮人人尽说江南好唐韦庄| 午夜免费鲁丝| 老汉色∧v一级毛片| 黄色配什么色好看| 国产精品秋霞免费鲁丝片| 日日摸夜夜添夜夜爱| 女人被躁到高潮嗷嗷叫费观| 国产探花极品一区二区| 只有这里有精品99| 亚洲欧洲精品一区二区精品久久久 | 国产成人精品福利久久| 少妇被粗大猛烈的视频| 中文字幕精品免费在线观看视频| 免费黄网站久久成人精品| 国产av精品麻豆| 爱豆传媒免费全集在线观看| 女性被躁到高潮视频| 久久国产精品男人的天堂亚洲| 老司机影院毛片| 国产老妇伦熟女老妇高清| 一区二区三区激情视频| 捣出白浆h1v1| 国产有黄有色有爽视频| tube8黄色片| 18+在线观看网站| 国产免费福利视频在线观看| 久久久久久久久久久久大奶| 亚洲国产av影院在线观看| www.熟女人妻精品国产| 日韩欧美精品免费久久| 曰老女人黄片| 欧美日韩亚洲高清精品| 一级爰片在线观看| 九草在线视频观看| 男女啪啪激烈高潮av片| 最近最新中文字幕免费大全7| 久久国产精品男人的天堂亚洲| 国产xxxxx性猛交| 国产女主播在线喷水免费视频网站| 免费人妻精品一区二区三区视频| 国产精品偷伦视频观看了| 黑人猛操日本美女一级片| 亚洲欧美清纯卡通| 亚洲av综合色区一区| 美女中出高潮动态图| 久久国产精品大桥未久av| 国产精品成人在线| av又黄又爽大尺度在线免费看| 青青草视频在线视频观看| 亚洲欧美精品综合一区二区三区 | 亚洲欧洲日产国产| 欧美精品av麻豆av| 热re99久久国产66热| 成年人午夜在线观看视频| 少妇 在线观看| 午夜免费观看性视频| 黄频高清免费视频| 亚洲国产av新网站| 黑人欧美特级aaaaaa片| 亚洲综合精品二区| 美女视频免费永久观看网站| 久久久精品免费免费高清| 黄色配什么色好看| 日韩av不卡免费在线播放| 一本色道久久久久久精品综合| 黄片播放在线免费| 亚洲欧美日韩另类电影网站| 十八禁高潮呻吟视频| 不卡视频在线观看欧美| 久久久亚洲精品成人影院| 男女国产视频网站| 王馨瑶露胸无遮挡在线观看| 一级a爱视频在线免费观看| 少妇的逼水好多| 亚洲欧美精品自产自拍| 一级毛片 在线播放| 91精品国产国语对白视频| 欧美亚洲日本最大视频资源| 日韩在线高清观看一区二区三区| 丝袜美腿诱惑在线| 桃花免费在线播放| 最近中文字幕高清免费大全6| 超碰成人久久| 久久国产精品男人的天堂亚洲| 欧美日韩一区二区视频在线观看视频在线| 人人妻人人澡人人看| 观看av在线不卡| 极品人妻少妇av视频| 午夜福利视频精品| 国产一区有黄有色的免费视频| 久久国产精品大桥未久av| 精品午夜福利在线看| 水蜜桃什么品种好| 久久久精品免费免费高清| 国产精品国产三级专区第一集| 久久久精品94久久精品| 欧美成人午夜精品| 国产成人一区二区在线| 国产极品粉嫩免费观看在线| av电影中文网址| 男人添女人高潮全过程视频| 亚洲一区二区三区欧美精品| 亚洲,一卡二卡三卡| 激情五月婷婷亚洲| 免费日韩欧美在线观看| 两个人免费观看高清视频| 久久99蜜桃精品久久| 91精品国产国语对白视频| 国产精品久久久av美女十八| 亚洲av日韩在线播放| 天堂8中文在线网| 欧美日韩亚洲高清精品| 亚洲一码二码三码区别大吗| 亚洲精品国产色婷婷电影| 久久久久久人人人人人| 黄频高清免费视频| 搡老乐熟女国产| 超碰97精品在线观看| 99精国产麻豆久久婷婷| 欧美国产精品一级二级三级| 丝袜在线中文字幕| 国产av码专区亚洲av| 可以免费在线观看a视频的电影网站 | 巨乳人妻的诱惑在线观看| 超色免费av| 一区二区三区乱码不卡18| 青春草视频在线免费观看| 国产精品成人在线| 欧美日韩国产mv在线观看视频| 一级爰片在线观看| 国产极品天堂在线| 91在线精品国自产拍蜜月| tube8黄色片| 欧美精品一区二区免费开放| 亚洲国产最新在线播放| 男女无遮挡免费网站观看| 超色免费av| 国产免费一区二区三区四区乱码| 国产乱人偷精品视频| 女性生殖器流出的白浆| 各种免费的搞黄视频| 一级黄片播放器| 国产成人精品福利久久| 一本色道久久久久久精品综合| 欧美日韩av久久| 90打野战视频偷拍视频| 中文精品一卡2卡3卡4更新| 国产亚洲一区二区精品| 成人亚洲欧美一区二区av| av免费在线看不卡| 秋霞在线观看毛片| 午夜福利在线免费观看网站| 一级片'在线观看视频| 婷婷色av中文字幕| 精品国产乱码久久久久久男人| 国产成人精品婷婷| 国产精品久久久久久精品古装| 欧美日本中文国产一区发布| 18禁裸乳无遮挡动漫免费视频| 最近的中文字幕免费完整| 久久鲁丝午夜福利片| av又黄又爽大尺度在线免费看| 国产黄色视频一区二区在线观看| 不卡视频在线观看欧美| 高清av免费在线| www.自偷自拍.com| 国产女主播在线喷水免费视频网站| 成人漫画全彩无遮挡| 免费黄频网站在线观看国产| 久久99一区二区三区| 久久ye,这里只有精品| 国产爽快片一区二区三区| 欧美日韩一区二区视频在线观看视频在线| 亚洲内射少妇av| 纯流量卡能插随身wifi吗| 人成视频在线观看免费观看| 啦啦啦中文免费视频观看日本| 老熟女久久久| 老鸭窝网址在线观看| 国产成人免费无遮挡视频| 久久精品aⅴ一区二区三区四区 | 国产高清国产精品国产三级| 丝瓜视频免费看黄片| 亚洲精品自拍成人| 久久 成人 亚洲| 成人午夜精彩视频在线观看| 日韩中字成人| 亚洲国产色片| 黄色一级大片看看| 男女高潮啪啪啪动态图| 97在线人人人人妻| tube8黄色片| 在线看a的网站| 久久国产亚洲av麻豆专区| 国产精品香港三级国产av潘金莲 | 国产精品久久久久久久久免| 男女免费视频国产| 丰满少妇做爰视频| 国产精品香港三级国产av潘金莲 | 午夜福利视频在线观看免费| 亚洲精品自拍成人| 久久国产亚洲av麻豆专区| kizo精华| 亚洲激情五月婷婷啪啪| 欧美日韩国产mv在线观看视频| 日本91视频免费播放| 久久精品亚洲av国产电影网| av电影中文网址| 岛国毛片在线播放| 一区二区av电影网| 在线免费观看不下载黄p国产| 老汉色av国产亚洲站长工具| 男女边吃奶边做爰视频| 伊人久久国产一区二区| 国产精品一国产av| 一二三四在线观看免费中文在| 在线精品无人区一区二区三| 亚洲精品国产一区二区精华液| 欧美亚洲日本最大视频资源| 91午夜精品亚洲一区二区三区| videossex国产| 老汉色av国产亚洲站长工具| 天堂8中文在线网| 欧美人与善性xxx| 99热网站在线观看| 久久久精品区二区三区| 高清黄色对白视频在线免费看| 国产高清国产精品国产三级| 波多野结衣av一区二区av| 成人手机av| 亚洲国产精品成人久久小说| 免费大片黄手机在线观看| 熟女少妇亚洲综合色aaa.| 久久久久精品性色| 国产探花极品一区二区| 久久久精品94久久精品| 波多野结衣av一区二区av| 日韩不卡一区二区三区视频在线| 男人爽女人下面视频在线观看| 精品久久久精品久久久| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | av又黄又爽大尺度在线免费看| 亚洲成人一二三区av| 麻豆av在线久日| 青春草国产在线视频| 国产精品成人在线| 性色av一级| 丝袜人妻中文字幕| 2021少妇久久久久久久久久久| 中文字幕人妻丝袜一区二区 | 国产精品三级大全| 美女午夜性视频免费| 可以免费在线观看a视频的电影网站 | 美女xxoo啪啪120秒动态图| 国产成人精品一,二区| 亚洲av成人精品一二三区| 欧美 日韩 精品 国产| 日本午夜av视频| 狠狠精品人妻久久久久久综合| 国产一区二区三区av在线| 侵犯人妻中文字幕一二三四区| 男女免费视频国产| 春色校园在线视频观看| 久久国产精品大桥未久av| 亚洲综合精品二区| 狠狠婷婷综合久久久久久88av| 日本爱情动作片www.在线观看| h视频一区二区三区| 中文精品一卡2卡3卡4更新| 美女午夜性视频免费| 高清av免费在线| 午夜免费观看性视频| 国产精品三级大全| 精品一区二区三卡| 久久久久视频综合| 久久国产精品男人的天堂亚洲| 99热网站在线观看| www.自偷自拍.com| 精品一区二区三区四区五区乱码 | 男女下面插进去视频免费观看| 丰满少妇做爰视频| 国产av精品麻豆| 美国免费a级毛片| 肉色欧美久久久久久久蜜桃| 不卡av一区二区三区| 亚洲第一区二区三区不卡| 亚洲欧美清纯卡通| 国产女主播在线喷水免费视频网站| 人妻人人澡人人爽人人| 精品国产乱码久久久久久男人| 男男h啪啪无遮挡| 久久久国产精品麻豆| 亚洲国产最新在线播放| 免费女性裸体啪啪无遮挡网站| 久久99热这里只频精品6学生| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 精品国产国语对白av| 叶爱在线成人免费视频播放| 亚洲欧美一区二区三区久久| 国产在线视频一区二区| 久久久久久久精品精品| 免费大片黄手机在线观看| av卡一久久| 国产av精品麻豆| 久久国内精品自在自线图片| 制服丝袜香蕉在线| 国产精品熟女久久久久浪| 精品福利永久在线观看| 亚洲精品第二区| 2022亚洲国产成人精品| 日韩一区二区三区影片| 香蕉国产在线看| 欧美亚洲日本最大视频资源| 久久99精品国语久久久| 欧美人与性动交α欧美软件| 99精国产麻豆久久婷婷| av不卡在线播放| 亚洲一码二码三码区别大吗| 亚洲欧美中文字幕日韩二区| 大话2 男鬼变身卡| 国产黄色视频一区二区在线观看| av网站免费在线观看视频| 一边摸一边做爽爽视频免费| 久久免费观看电影| 亚洲欧美日韩另类电影网站| 如日韩欧美国产精品一区二区三区| 韩国av在线不卡| 国产一区二区三区av在线| 韩国精品一区二区三区| 亚洲精品乱久久久久久| 一区二区日韩欧美中文字幕| 9色porny在线观看| 久久久国产精品麻豆| 亚洲四区av| 汤姆久久久久久久影院中文字幕| 黄网站色视频无遮挡免费观看| 久久99精品国语久久久| 久热久热在线精品观看| 日韩不卡一区二区三区视频在线| 成年av动漫网址| 一本色道久久久久久精品综合| 亚洲欧美色中文字幕在线| 大片电影免费在线观看免费| 欧美黄色片欧美黄色片| 国产免费视频播放在线视频| 精品国产露脸久久av麻豆| 免费日韩欧美在线观看| 天天操日日干夜夜撸| 秋霞在线观看毛片| 制服诱惑二区| 亚洲伊人色综图| 中文精品一卡2卡3卡4更新| 午夜福利,免费看| 国产成人免费无遮挡视频| 日本免费在线观看一区| 国产黄色视频一区二区在线观看| 日韩一本色道免费dvd| 亚洲五月色婷婷综合| 男人操女人黄网站| 乱人伦中国视频|