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

    纖鋅礦結(jié)構(gòu)AlN、GaN及ZnO自發(fā)極化的第一性原理研究

    2017-12-02 03:55:42牛海波竹有章李冠強(qiáng)
    關(guān)鍵詞:偶極矩鋅礦閃鋅礦

    牛海波, 竹有章, 李冠強(qiáng)

    (1.西安交通大學(xué)城市學(xué)院 物理教學(xué)部, 陜西 西安 710018; 2.陜西科技大學(xué) 文理學(xué)院, 陜西 西安 710021)

    纖鋅礦結(jié)構(gòu)AlN、GaN及ZnO自發(fā)極化的第一性原理研究

    牛海波1, 竹有章1, 李冠強(qiáng)2

    (1.西安交通大學(xué)城市學(xué)院 物理教學(xué)部, 陜西 西安 710018; 2.陜西科技大學(xué) 文理學(xué)院, 陜西 西安 710021)

    根據(jù)現(xiàn)代極化理論,分別利用Berry phase方法和最大局域化Wannier函數(shù)方法系統(tǒng)計(jì)算了纖鋅礦結(jié)構(gòu)AlN、GaN及ZnO中的自發(fā)極化,并從電子項(xiàng)和離子項(xiàng)引起的極化變化具體分析了自發(fā)極化的起源.研究表明AlN自發(fā)極化中電子項(xiàng)的貢獻(xiàn)占據(jù)主導(dǎo)地位,而GaN及ZnO自發(fā)極化中以離子項(xiàng)的貢獻(xiàn)為主.研究發(fā)現(xiàn)其他文獻(xiàn)計(jì)算自發(fā)極化時(shí),由于計(jì)算模型和參考模型使用相同的體積,導(dǎo)致計(jì)算結(jié)果偏小.利用Wannier中心,從結(jié)構(gòu)中最小重復(fù)單元的電偶極矩出發(fā)對原有計(jì)算公式進(jìn)行了修正,使得自發(fā)極化的計(jì)算結(jié)果更趨合理.研究發(fā)現(xiàn)在特定建模條件下,可以不用計(jì)算參照模型中的極化,利用Wannier中心確定纖鋅礦結(jié)構(gòu)中正負(fù)電荷的重心,通過經(jīng)典的靜電學(xué)理論直接計(jì)算出自發(fā)極化,直觀解釋了自發(fā)極化的形成.

    纖鋅礦結(jié)構(gòu); 自發(fā)極化; Wannier中心; 第一性原理

    0 引言

    AlN、GaN及ZnO作為典型的第三代半導(dǎo)體材料,具有直接帶隙且?guī)秾挾却蟮娘@著優(yōu)點(diǎn)[1,2],因而呈現(xiàn)出優(yōu)異的發(fā)光特性,在藍(lán)光、深紫外發(fā)光器件及光電探測器件等領(lǐng)域具有重要的應(yīng)用價(jià)值.如隨著高質(zhì)量P型GaN(禁帶寬度3.4 eV)制造工藝的突破[3],GaN基藍(lán)綠光發(fā)光二極管已經(jīng)研制成功并商品化,占據(jù)了發(fā)光二極管市場的絕大多數(shù)份額.2006年,Taniyasu等[4]成功制作出直接基于AlN的發(fā)射波長為210 nm的深紫外發(fā)光二極管,成為迄今為止在氮化物半導(dǎo)體領(lǐng)域獲得的最短發(fā)射波長.通過將AlN、GaN及InN組成三元或四元合金,發(fā)光譜更是覆蓋了從紅外到紫外的整個(gè)區(qū)域[5].這些優(yōu)點(diǎn)使得AlN、GaN及ZnO成為目前最有發(fā)展前景的半導(dǎo)體發(fā)光材料,吸引了眾多研究者的廣泛關(guān)注,并得到迅猛發(fā)展.

    Ⅲ族氮化物及ZnO等寬禁帶半導(dǎo)體材料的突出特征是晶體結(jié)構(gòu)中存在較大的自發(fā)極化.由于這些晶體自然狀態(tài)下的穩(wěn)定結(jié)構(gòu)多為六方對稱的纖鋅礦結(jié)構(gòu),這種結(jié)構(gòu)的特點(diǎn)是每一個(gè)原子與其最近鄰的四個(gè)原子所成的化學(xué)鍵中,沿[0001]方向的鍵比其他三個(gè)稍長,使得原胞中正負(fù)電荷中心不重合,產(chǎn)生電偶極矩,形成了內(nèi)建電場,即存在自發(fā)極化.極化效應(yīng)對材料的光電性質(zhì)存在重要影響.

    近年來受到廣泛關(guān)注的AlGaN/GaN異質(zhì)結(jié)材料,正是因?yàn)榻Y(jié)構(gòu)中內(nèi)生的極化場,導(dǎo)致即使在非故意摻雜情況下也存在超高密度的二維電子氣(載流子濃度高達(dá)1013/cm2)[6,7],非常適合制作一些特殊性能的傳感器件.另一方面,極化場對摻雜元素的濃度及摻雜位置也有較大的影響[8,9],如在異質(zhì)結(jié)中極化場使得氫元素向表面運(yùn)動(dòng);強(qiáng)極化場還會(huì)造成能帶彎曲[10],阻礙電子和空穴對的復(fù)合,降低器件的發(fā)光效率[11].在太陽能電池中,極化場降低了吸收層的總電場,減小了光生載流子的收集效果,從而降低了太陽能電池的性能[12].極化場是研究寬禁帶半導(dǎo)體性質(zhì)時(shí)不可回避的本質(zhì)問題,研究極化場性質(zhì),并對其加以調(diào)控,避短揚(yáng)長,將加快寬禁帶半導(dǎo)體的發(fā)展與應(yīng)用.

    本文根據(jù)現(xiàn)代極化理論,通過構(gòu)建參照模型和計(jì)算模型,分別利用Berry phase及最大局域化Wannier函數(shù)(MLWF)方法[13-15]對纖鋅礦AlN、 GaN及ZnO的自發(fā)極化進(jìn)行了系統(tǒng)的計(jì)算,研究了影響晶體自發(fā)極化大小的各種因素,同時(shí)對其他文獻(xiàn)報(bào)道的結(jié)果進(jìn)行了分析,指出文獻(xiàn)中所用方法的不足并進(jìn)行了公式修正,使得計(jì)算結(jié)果更加合理.此外根據(jù)計(jì)算得到的Wannier 中心,從三種晶體結(jié)構(gòu)的最小重復(fù)單元中電荷的重心分布出發(fā),重新計(jì)算了三種結(jié)構(gòu)的自發(fā)極化,對自發(fā)極化的形成給出了更直觀的解釋.

    1 計(jì)算方法及模型

    根據(jù)現(xiàn)代極化理論,為了計(jì)算纖鋅礦結(jié)構(gòu)中的自發(fā)極化值,首先需要建立一種正負(fù)電荷中心重合的對稱結(jié)構(gòu)作為初始結(jié)構(gòu)(也稱參照結(jié)構(gòu)),這樣才能保證參照結(jié)構(gòu)中的自發(fā)極化為0.纖鋅礦結(jié)構(gòu)與該參照結(jié)構(gòu)的極化值之差才是纖鋅礦結(jié)構(gòu)中的自發(fā)極化.立方對稱的閃鋅礦結(jié)構(gòu)滿足上述建模要求,一般被用作參照模型.這種思想可以表述為,

    ΔP=Pw-Pz

    (1)

    式(1)中:ΔP為纖鋅礦結(jié)構(gòu)中的自發(fā)極化,Pw、Pz分別為纖鋅礦和閃鋅礦形式上的極化值,該值可用式(2)表示成離子項(xiàng)Pion和電子項(xiàng)Pel之和,

    P=Pion+Pel

    (2)

    由于介質(zhì)中的離子分布非常局域,因此離子項(xiàng)的貢獻(xiàn)可以用經(jīng)典的點(diǎn)電荷模型進(jìn)行計(jì)算,即

    (3)

    對Pel的計(jì)算是現(xiàn)代極化理論的精華,Berry phase方法給出的結(jié)果為[15],

    (4)

    n代表了電子占據(jù)的能帶,利用最大局域化Wannier函數(shù)方法,Pel可表示為[13],

    (5)

    式(5)中:ωj和qj為第j個(gè)Wannier中心的位置和所帶的電荷數(shù),V為晶胞體積.根據(jù)現(xiàn)代極化理論,在系統(tǒng)絕熱緩慢地從初始狀態(tài)到末了狀態(tài)的演變過程中,由于每個(gè)狀態(tài)的極化值在計(jì)算時(shí)量子幾何相位(貝利相)是對2π取模計(jì)算出來的,導(dǎo)致會(huì)產(chǎn)生極化量子數(shù)eR/V[15],其中R為晶格參數(shù),V為單胞體積.如果每個(gè)狀態(tài)的極化量子數(shù)不同,將對極化值的改變量ΔP產(chǎn)生干擾.而在實(shí)際計(jì)算時(shí)如本文所研究的自發(fā)極化,由于|ΔP|=|eR/V|,即演變過程中每個(gè)狀態(tài)具有相同的量子數(shù),因此ΔP可避開極化量子數(shù)的干擾,使得自發(fā)極化具有確定值.

    由于纖鋅礦結(jié)構(gòu)沿[0001]方向按照ABAB順序堆垛而成,而閃鋅礦結(jié)構(gòu)沿[111]方向按照ABCABC順序堆垛,因此為了不僅保證兩種結(jié)構(gòu)中原子數(shù)目一致,而且要求計(jì)算量最小,構(gòu)建了1×1×3的纖鋅礦型晶胞及1×1×2的閃鋅礦型晶胞,每種模型共有12個(gè)原子.以AlN為例,分別包含了6個(gè)N原子及6個(gè)Al原子.建模時(shí)首先構(gòu)建閃鋅礦參照模型,然后按照ABAB的堆垛順序?qū)⒄漳P椭械腘、Al原子進(jìn)行移動(dòng),初步構(gòu)建成為纖鋅礦模型.最后對閃鋅礦和纖鋅礦模型進(jìn)行充分的結(jié)構(gòu)優(yōu)化(晶格參數(shù)、原子位置都允許變化).所建模型如圖1所示.

    (a) 閃鋅礦結(jié)構(gòu) (b) 纖鋅礦結(jié)構(gòu)

    與其他文獻(xiàn)計(jì)算中所用到的大的超胞模型相比,本文所建模型原子數(shù)目少,因此計(jì)算量減小,此外本文對閃鋅礦參照模型和纖鋅礦模型都進(jìn)行了完全的結(jié)構(gòu)優(yōu)化,使得兩種模型各自呈現(xiàn)出典型的閃鋅礦及纖鋅礦結(jié)構(gòu)特點(diǎn),而不必保持相同的體積,這是提高計(jì)算合理性及準(zhǔn)確性的重要保證.

    采用第一性原理計(jì)算軟件Quantum-Espresso[16]進(jìn)行結(jié)構(gòu)優(yōu)化、靜態(tài)自洽及非自洽計(jì)算,MLWF由Wannier90軟件包進(jìn)行計(jì)算[17,18].利用廣義梯度近似(GGA)的PBE來處理電子之間的交換關(guān)聯(lián)能,選擇的贗勢為Vanderbilt超軟贗勢.平面波截?cái)嗄転?0 Ry,選取Monkorst-Park特殊K點(diǎn)對全布里淵區(qū)求和.ZnO、GaN自發(fā)極化計(jì)算中Zn、Ga的3d態(tài)電子也作為價(jià)電子,相應(yīng)的利用Wannier90計(jì)算MLWF時(shí)也包含了Zn、Ga的3d能帶.利用Berry phase方法計(jì)算[0001]方向上的極化時(shí),該方向上需要密集K點(diǎn).由于過密的K點(diǎn)并不能帶來更高的準(zhǔn)確性,反而會(huì)使計(jì)算量大幅增加,本文采用9×9×7對K點(diǎn)進(jìn)行了加密,計(jì)算中總能變化收斂的標(biāo)準(zhǔn)為1.0×10-6eV,原子間相互作用力的收斂標(biāo)準(zhǔn)為0.05 eV/nm.

    2 結(jié)果與討論

    極化性質(zhì)對結(jié)構(gòu)的晶格參數(shù)非常敏感,因此在計(jì)算自發(fā)極化前首先對AlN、GaN、ZnO的原胞進(jìn)行了仔細(xì)的結(jié)構(gòu)馳豫,然后再對以原胞為基礎(chǔ)構(gòu)建的閃鋅礦及纖鋅礦模型進(jìn)行優(yōu)化,得到的三種晶體原胞的晶格參數(shù)如表1所示.其中c/a為晶格參數(shù)c與垂直方向的晶格參數(shù)a的比值,參數(shù)u為纖鋅礦結(jié)構(gòu)中平行于[0001]方向的鍵長與晶格參數(shù)c的比值,反映了[0001]方向原子層的間距.三種原胞的晶格參數(shù)實(shí)驗(yàn)值由X射線衍射測量得到[1,2],圖2給出了優(yōu)化后的AlN閃鋅礦及纖鋅礦結(jié)構(gòu)的鍵角.

    表1 纖鋅礦AlN、GaN及ZnO原胞晶格參數(shù)及兩種結(jié)構(gòu)鍵長

    (a) 閃鋅礦 (b) 釬鋅礦

    從表1可以看出,計(jì)算得到的纖鋅礦結(jié)構(gòu)原胞的晶格參數(shù)與實(shí)驗(yàn)值非常接近,即使在計(jì)算值與實(shí)驗(yàn)值相差最大的ZnO中,相對變化也只有0.89%.三種晶體的晶格參數(shù)計(jì)算值普遍稍高于實(shí)驗(yàn)值,這與第一性原理計(jì)算中交換關(guān)聯(lián)能使用GGA近似有關(guān).GGA近似比較適合電子密度非均勻的體系,一般情況下采用GGA近似計(jì)算得到的晶格參數(shù)均會(huì)輕微增加,這并不影響結(jié)構(gòu)的對比分析.觀察兩種原胞中的鍵長,在纖鋅礦結(jié)構(gòu)中,垂直于[0001]方向的三個(gè)鍵長相等,且均小于[0001]方向上的鍵長,其鍵角為108.20 °、110.71 °,呈現(xiàn)出典型的纖鋅礦結(jié)構(gòu)特點(diǎn);在閃鋅礦參照結(jié)構(gòu)中,四個(gè)鍵長一致,鍵角(如圖2(a)所示)均為正四面體的標(biāo)準(zhǔn)鍵角109.47 °,符合閃鋅礦的正四面體結(jié)構(gòu)特征,說明了結(jié)構(gòu)優(yōu)化的合理性以及所建模型的準(zhǔn)確,也保證了參照結(jié)構(gòu)不存在自發(fā)極化.

    2.1 利用Berry phase方法計(jì)算自發(fā)極化

    首先分別對建立的閃鋅礦參照結(jié)構(gòu)和纖鋅礦結(jié)構(gòu)進(jìn)行靜態(tài)自洽計(jì)算,然后設(shè)置[0001]方向?yàn)闃O化方向進(jìn)行非靜態(tài)計(jì)算,可得到兩種結(jié)構(gòu)下的貝利相及極化值.需要注意的是,根據(jù)現(xiàn)代極化理論,晶體在某一狀態(tài)時(shí)的極化值具有多值性,此時(shí)所求得的只是形式上的極化值,并不是該結(jié)構(gòu)中的自發(fā)極化,它也沒有確切的物理意義.只有纖鋅礦結(jié)構(gòu)與閃鋅礦結(jié)構(gòu)中形式上極化的差值才是纖鋅礦結(jié)構(gòu)的自發(fā)極化.表2中列出了AlN、GaN及ZnO 閃鋅礦及纖鋅礦結(jié)構(gòu)中的貝利相(離子項(xiàng)及電子項(xiàng)分別列出)、對應(yīng)的形式上的極化值P及纖鋅礦結(jié)構(gòu)中的自發(fā)極化ΔP,并與其它文獻(xiàn)進(jìn)行了對比.三種晶體中自發(fā)極化均為負(fù)值,說明它們的極化方向與[0001]方向相反.表2中貝利相φ與極化值P之間的關(guān)系為,

    (6)

    式(6)中:V為超胞的體積,R為超胞[0001]方向上的晶格矢量,q為電荷電量.由表2可以發(fā)現(xiàn),AlN、GaN的自發(fā)極化計(jì)算值與其他文獻(xiàn)的結(jié)果非常接近,Laehnemann等[21]利用實(shí)驗(yàn)測得纖鋅礦GaN的自發(fā)極化為-0.022±0.007 C/m2,也與本文的計(jì)算結(jié)果符合較好,說明了所用計(jì)算模型及Berry phase方法的正確性.ZnO的計(jì)算值接近文獻(xiàn)[19]的報(bào)道而與文獻(xiàn)[20]差別較大,主要是因?yàn)橛?jì)算中使用的交換關(guān)聯(lián)能及晶格參數(shù)與文獻(xiàn)中不同所導(dǎo)致,特別是參數(shù)u值不同,因?yàn)樽园l(fā)極化對u值更敏感.我們利用文獻(xiàn)[20]給出的u值代替表1中ZnO的u值重新進(jìn)行了計(jì)算,最后得出ZnO 的自發(fā)極化值為-0.036 C/m2,接近文獻(xiàn)[20]的結(jié)果.

    表2 AlN、GaN及ZnO的貝利相及自發(fā)極化

    注意到計(jì)算出來的AlN、 GaN的自發(fā)極化值(取極化的絕對值)均稍大于文獻(xiàn)[19]報(bào)道的結(jié)果,我們認(rèn)為這是由于文獻(xiàn)[19]中沒有考慮到閃鋅礦參照結(jié)構(gòu)及纖鋅礦兩種超胞具有不同的體積,而是統(tǒng)一使用了纖鋅礦超胞的體積所導(dǎo)致的.具體來說,如圖1所示,兩種模型盡管原子數(shù)目保持一致,但是由于結(jié)構(gòu)對稱性的要求,他們不可能保持相同的晶格參數(shù)與體積.而極化強(qiáng)度是用單位體積中的偶極矩來表征的,如果直接用計(jì)算出的兩種結(jié)構(gòu)的極化強(qiáng)度做差來計(jì)算自發(fā)極化,就表示纖鋅礦和閃鋅礦模型具有相同的體積,這正是引起誤差的原因.準(zhǔn)確的做法是應(yīng)該用兩種結(jié)構(gòu)中的偶極矩做差再除以對應(yīng)的纖鋅礦模型的體積,這是很多之前的報(bào)道中所忽略的.進(jìn)一步分析表3列出的兩種結(jié)構(gòu)的超胞體積,可以看到兩種超胞的體積相差很小,AlN的相對體積變化在0.28%以內(nèi),但是這種微小的體積差異卻是由超胞的晶格參數(shù)c變化較明顯引起的,分析表3中給出的晶格參數(shù)的相對變化,AlN達(dá)到了1.3%, GaN達(dá)到了0.25%,這已經(jīng)足以影響到最終的自發(fā)極化.而且可以看到正是因?yàn)锳lN晶格參數(shù)的c相對變化明顯大于GaN,導(dǎo)致了本文所計(jì)算的AlN自發(fā)極化數(shù)據(jù)與文獻(xiàn)報(bào)道結(jié)果差異更顯著.

    表3 AlN、GaN及ZnO閃鋅礦與纖鋅礦結(jié)構(gòu)晶格參數(shù)、超胞體積的變化

    對于上述計(jì)算得到的自發(fā)極化大于文獻(xiàn)報(bào)道的情況,也可以通過壓電極化進(jìn)行定量的解釋,當(dāng)晶體受到外應(yīng)力引起結(jié)構(gòu)變形時(shí),[0001]晶向上極化的變化為[19],

    δP=2e31ε1+e33ε3

    (7)

    式(7)中:ε1=(a-a0)/a0為晶體a軸受到的應(yīng)力,ε3=(c-c0)/c0為晶體c軸受到的應(yīng)力,e31,e33為壓電系數(shù).現(xiàn)在以AlN為例來進(jìn)行說明,令閃鋅礦結(jié)構(gòu)參數(shù)使用兩組數(shù)據(jù),一組為其他文獻(xiàn)直接使用的纖鋅礦數(shù)據(jù)(記為狀態(tài)1),一組為本章所優(yōu)化過的閃鋅礦結(jié)構(gòu)數(shù)據(jù)(記為狀態(tài)2),由表3可以看到,AlN閃鋅礦結(jié)構(gòu)的晶格參數(shù)a小于纖鋅礦的晶格參數(shù)a0,而參數(shù)c明顯大于纖鋅礦結(jié)構(gòu)中的c0參數(shù),因此當(dāng)閃鋅礦結(jié)構(gòu)從狀態(tài)1變化到狀態(tài)2時(shí),ε1lt;0,ε3gt;0,而e33gt;0,e31lt;0,因此δP=P2-P1gt;0,由于P1、P2均為負(fù)值,因此|P2|lt;|P1|,再利用纖鋅礦結(jié)構(gòu)中的極化值與閃鋅礦中的極化值做差求自發(fā)極化時(shí),使用狀態(tài)1的數(shù)據(jù)明顯會(huì)使最后的自發(fā)極化值減小.

    表4 結(jié)構(gòu)從閃鋅礦變化到纖鋅礦

    2.2 利用MLWF方法分析自發(fā)極化

    上節(jié)利用Berry phase方法計(jì)算了晶體的自發(fā)極化,Vanderbilt等也提出了最大局域化Wannier函數(shù)方法來計(jì)算極化,通過引進(jìn)局域標(biāo)準(zhǔn)求得最局域的Wannier函數(shù),解決了Wannier函數(shù)長期以來存在不唯一性的缺點(diǎn),使得MLWF方法在諸多方面得到了迅速運(yùn)用.在MLWF方法中,電荷可以形象地認(rèn)為局域在Wannier函數(shù)空間分布的幾何中心上,該幾何中心即為Wannier中心,如此一來,可將Wannier中心看做點(diǎn)電荷,從而整個(gè)晶體就可以分為帶正電的離子和帶負(fù)電的Wannier中心,晶體自發(fā)極化可用下式來計(jì)算.

    (8)

    式(8)中:ωj和qj為第j個(gè)Wannier中心的位置和所帶的電荷數(shù);ri和qi分別為第i個(gè)離子的位置與電荷量,V為晶胞體積.對于AlN晶體,可以認(rèn)為是由+3價(jià)的Al離子、+5價(jià)的N離子以及-2價(jià)的Wannier中心組成.在計(jì)算纖鋅礦GaN、ZnO中的自發(fā)極化時(shí),由于將Ga、Zn原子的3d態(tài)電子也作為價(jià)電子,因此最后計(jì)算時(shí)需要將GaN晶體分為+13價(jià)的Ga離子、+5價(jià)的N離子及-2價(jià)的Wannier中心,ZnO晶體則分成+12價(jià)的Zn離子、+6價(jià)的O離子及-2價(jià)的Wannier中心.在構(gòu)造Wannier函數(shù)時(shí),初始空間分布選擇sp3雜化軌道與d軌道,Wannier初始中心則選擇N、O位置及Ga、Zn位置.圖3給出了計(jì)算得到的AlN及ZnO的Wannier中心分布,可以看到ZnO中一部分Wannier中心幾乎與Zn原子位置重合,這是從3d態(tài)能帶得到的Wannier中心,反映了3d態(tài)電子的強(qiáng)局域性,其空間延展函數(shù)值也只有0.26 ?2.事實(shí)上在計(jì)算中我們發(fā)現(xiàn),由于d電子局域在Ga、 Zn原子位置,導(dǎo)致它對GaN、ZnO的自發(fā)極化貢獻(xiàn)很弱,在GaN中,d電子對自發(fā)極化的貢獻(xiàn)僅占1%,在ZnO中,d電子對自發(fā)極化的貢獻(xiàn)也僅有0.81%,因此主要的貢獻(xiàn)來自于最外層的價(jià)電子.

    (a) 纖鋅礦AIN Wannier中心及Wannier函數(shù)空間分布

    (b) 纖鋅礦ZnO Wannier中心及Wannier函數(shù)空間分布

    根據(jù)Wannier中心位置計(jì)算出閃鋅礦及纖鋅礦結(jié)構(gòu)中的形式上的極化值,二者之差即為纖鋅礦AlN的自發(fā)極化.為了進(jìn)行比較,分別計(jì)算了纖鋅礦結(jié)構(gòu)AlN[0001]晶向,以及與它垂直的兩個(gè)晶向上的自發(fā)極化以及離子部分ΔPion和電子部分ΔPel對極化的貢獻(xiàn),結(jié)果如表5所示.同時(shí)也將Berry phase方法計(jì)算得到的三種晶體的自發(fā)極化也列于表中進(jìn)行對比.

    式(8)是諸多文獻(xiàn)中所使用的,在分析晶體保持體積不變時(shí)的極化性質(zhì)時(shí),它是完全成立的,但是我們發(fā)現(xiàn)在計(jì)算諸如本文所涉及的自發(fā)極化時(shí),直接利用該式會(huì)使計(jì)算結(jié)果明顯偏離公認(rèn)值,如在AlN自發(fā)極化計(jì)算中,直接使用式(8),最后的結(jié)果為-0.033 C/m2,遠(yuǎn)小于利用Berry phase方法得到的結(jié)果-0.086 C/m2,造成這一差異的原因在于閃鋅礦參照結(jié)構(gòu)與纖鋅礦結(jié)構(gòu)超胞的體積并不一致,晶格常數(shù)如c參數(shù)也不相同,這種不一致使得直接利用式(8)會(huì)帶來明顯的差異.在其他文獻(xiàn)中,計(jì)算自發(fā)極化時(shí)由于閃鋅礦與纖鋅礦超胞都采用同一個(gè)體積,因此不存在這個(gè)情況或者說忽略了,但是兩種超胞采用同一個(gè)體積計(jì)算自發(fā)極化是不合適的,因?yàn)槿绱藰?gòu)建的閃鋅礦超胞并不是完全的正面體結(jié)構(gòu),導(dǎo)致參照結(jié)構(gòu)自發(fā)極化不為零.那么為什么晶格參數(shù)不一致會(huì)導(dǎo)致這樣大的差異,我們認(rèn)為具體還在于對利用MLWF計(jì)算電偶極矩的內(nèi)涵沒有理解透徹,事實(shí)上,計(jì)算閃鋅礦及纖鋅礦超胞中的電偶極矩應(yīng)該從結(jié)構(gòu)中的最小重復(fù)單元如N-Al對中的偶極矩入手進(jìn)行分析.圖4以纖鋅礦AlN為例給出了晶體中最小重復(fù)單元中的偶極矩示意圖.

    表5 AlN、GaN及ZnO三個(gè)晶向上的自發(fā)極化

    圖4 纖鋅礦AlN結(jié)構(gòu)最小重復(fù)單元中 的偶極矩(黑色箭頭)示意圖

    現(xiàn)在以圖4方框所示的最小重復(fù)單元即N-Al對的偶極矩為例進(jìn)行分析,WF表示W(wǎng)annier中心.MLWF方法將AlN晶體看作由-2價(jià)的Wannier中心、+5價(jià)的N離子以及+3價(jià)的Al離子組成的點(diǎn)電荷集合.為了計(jì)算電偶極矩,可以進(jìn)一步將N、 Al離子各分成帶電量相同的四等份,即將+5價(jià)的N離子看作由4個(gè)帶電量分別為+1.25的離子組成,將+3價(jià)的Al離子看作由4個(gè)帶電量分別為+0.75的離子組成,而每個(gè)Wannier中心則相應(yīng)看作由帶電量為-1.25及-0.75的兩個(gè)離子組成,這樣就可以將N-Al對中的偶極矩表示為圖4所示的偶極矩(由箭頭表示)的合成.以上述分析為基礎(chǔ),AlN超胞中總的電偶極矩即為結(jié)構(gòu)中所有N-Al 對中的偶極矩之和.值得注意的是,在超胞中,處理最上方或者最下方的N-Al的偶極矩時(shí),根據(jù)所建超胞的實(shí)際情況需要增加或者減少一個(gè)晶格常數(shù)項(xiàng),而該項(xiàng)在兩種超胞體積不同的情形下非常重要,在其他文獻(xiàn)中閃鋅礦與纖鋅礦采用同一體積,會(huì)使得這一項(xiàng)兩兩相減而正好消掉,因此被式(8)所忽視,我們用圖5來說明.

    (a) N面 (b) Al面

    在計(jì)算圖5(a)所示的黑色方框中N-Al對的偶極矩時(shí),由于N原子已經(jīng)位于晶胞的最上方,因此計(jì)算偶極矩時(shí)所需的三個(gè)Al原子就必然位于晶胞之外,根據(jù)對稱性,超胞外Al原子的Z坐標(biāo)就需要由晶胞中最下方的Al原子的Z坐標(biāo)加上晶格常數(shù)c得來,而且由于需要三個(gè)晶胞外的Al原子來計(jì)算偶極矩,因此相比直接利用式(8)計(jì)算,黑框中的N-Al對的偶極矩就增加了ΔPdipol,可寫為,

    (9)

    式(9)中:ZAl為最下方的Al原子的Z坐標(biāo), 為纖鋅礦超胞[0001]方向的晶格參數(shù),對于整個(gè)纖鋅礦超胞中的偶極矩,則表示為,

    (10)

    因此最終纖鋅礦結(jié)構(gòu)中的自發(fā)極化計(jì)算修正為,

    (11)

    cZB為閃鋅礦超胞[0001]方向的晶格參數(shù),GaN、ZnO的自發(fā)極化計(jì)算的修正與上述情況相似,只不過在ZnO中由于O為+6價(jià),此時(shí)需要把式(11)中的3/4更換為2/4,本文計(jì)算中所建閃鋅礦及纖鋅礦超胞模型均為圖5(a)的情形,即N、O原子位于最上方.因此根據(jù)修正后的式(11),計(jì)算得到的結(jié)果如表5所示,實(shí)現(xiàn)了與Berry phase方法得到的結(jié)果一致.對于圖5(b)中所示的Al面結(jié)構(gòu),其修正與N面相似,不再贅述.

    2.3 利用電荷重心自發(fā)極化

    前面利用Berry phase及MLWF方法計(jì)算自發(fā)極化時(shí),都需要計(jì)算閃鋅礦參照結(jié)構(gòu)中的極化值,然后纖鋅礦與閃鋅礦形式上的極化值之差為自發(fā)極化,這樣計(jì)算起來稍顯復(fù)雜,分析自發(fā)極化的形成也不直觀.MLWF方法可形象地將整個(gè)晶體分為局域的帶正電的離子與帶負(fù)電的Wannier中心,我們發(fā)現(xiàn)在當(dāng)前建模條件下,利用纖鋅礦結(jié)構(gòu)的對稱性及Wannier中心也是幾乎對稱分布的特點(diǎn),可以方便地確定正負(fù)離子的重心,然后將其看成是經(jīng)典靜電學(xué)中的正負(fù)中心點(diǎn)電荷,帶電量分別為+8e、-8e,然后根據(jù)

    (12)

    可以直接得到纖鋅礦結(jié)構(gòu)中的自發(fā)極化,而不用再計(jì)算閃鋅礦參照結(jié)構(gòu)中的極化,式(12)中q=-8e,為負(fù)中心點(diǎn)電荷帶電量.這樣更能直觀地解釋自發(fā)極化的形成.下面以AlN結(jié)構(gòu)中一個(gè)重復(fù)單元為例進(jìn)行說明.圖6給出了AlN閃鋅礦及纖鋅礦結(jié)構(gòu)中的鍵長及Wannier中心與N原子之間的距離.在AlN結(jié)構(gòu)中,N原子軌道sp3雜化,與周圍四個(gè)Al原子成鍵,根據(jù)MLWF的特點(diǎn),可以認(rèn)為圖6中每個(gè)Al帶電量為0.75.在閃鋅礦結(jié)構(gòu)中,由于四個(gè)N-Al鍵的鍵長都相同,N原子位于正面體的中心位置,因此正電荷的重心位于N原子上;而四個(gè)Wannier中心也是對稱的分布在N原子周圍,其電荷重心也位于N原子上,因此整體上正負(fù)電荷重心重合,不存在電偶極矩,也就沒有自發(fā)極化.在纖鋅礦結(jié)構(gòu)中,沿[0001]方向的鍵長大于其余三個(gè)鍵長,因此四個(gè)Al原子的正電荷重心偏向N原子上方,導(dǎo)致最后總的正電荷重心處在N原子上方;此外四個(gè)Wannier中心的電荷重心位置也偏離了N原子,這是由于[0001]方向上Wannier中心比其余三個(gè)Wannier中心更靠近N原子,使得整個(gè)負(fù)電荷重心處于N原子的下方.因此結(jié)構(gòu)中整體上正負(fù)電荷重心明顯分離,產(chǎn)生了電偶極矩,方向豎直向上,與[0001]方向相反,這與表4得到的分析結(jié)果完全一致.

    圖6 AlN兩種結(jié)構(gòu)中鍵長及Wannier中心

    為了分析因正負(fù)電荷中心分離產(chǎn)生的偶極矩大小,根據(jù)結(jié)構(gòu)的對稱性計(jì)算了圖6(b)所示纖鋅礦結(jié)構(gòu)的正負(fù)電荷的重心位置(僅考慮Z坐標(biāo)),以及據(jù)此計(jì)算得來的自發(fā)極化,如表6所示,由于Ga、Zn原子的d電子局域在Ga、Zn原子的位置,對極化的貢獻(xiàn)很小,因此計(jì)算時(shí)不考慮d電子.

    表6 纖鋅礦結(jié)構(gòu)正負(fù)電荷重心位置、自發(fā)極化

    首先從自發(fā)極化的計(jì)算結(jié)果可以看到與利用MLWF方法的結(jié)果幾乎一致,說明了分析正負(fù)電荷重心過程的合理性,也說明利用這種通過正負(fù)電荷重心計(jì)算自發(fā)極化的方法可以不用再構(gòu)建閃鋅礦參照結(jié)構(gòu),符合人們傳統(tǒng)直觀的認(rèn)識.從正負(fù)電荷重心之間的間距Δr可以看出,AlN在三種晶體中正負(fù)電荷重心間距最大,達(dá)到GaN的2.6倍,表明了AlN較強(qiáng)的極性.

    3 結(jié)論

    根據(jù)現(xiàn)代極化理論,通過構(gòu)建閃鋅礦參照模型及纖鋅礦計(jì)算模型,利用Berry phase方法及MLWF方法系統(tǒng)地計(jì)算了纖鋅礦結(jié)構(gòu)的AlN、 GaN及ZnO三種晶體的自發(fā)極化,與其它報(bào)道進(jìn)行了比較.從電子部分及離子部分的貢獻(xiàn)分析了自發(fā)極化的形成,發(fā)現(xiàn)在AlN中電子部分的貢獻(xiàn)占主導(dǎo)地位,而GaN及ZnO中則已離子部分的貢獻(xiàn)為主.同時(shí)發(fā)現(xiàn)其他文獻(xiàn)報(bào)道中所用計(jì)算方法的不足,會(huì)導(dǎo)致計(jì)算結(jié)果偏小,指出這是由于將閃鋅礦及纖鋅礦超胞采用同一體積所導(dǎo)致的,從電偶極矩的角度出發(fā)分析了原因,并給出了修正公式.

    [1] Strite S,Morkoc H.GaN,AlN and InN:A review[J].Journal of Vacuum Science amp; Technology B,1992,10(4):1 237-1 266.

    [2] Ozgur U,Alivov Y I,Liu C,et al.A comprehensive review of ZnO materials and devices[J].Journal of Applied Physics,2005,98(4):1-103.

    [3] Nakamura S,Iwasa N,Senoh M,et al.Hole compensation mechanism of p-type GaN films[J].Japanese Journal of Applied Physics,1992,31(5A):1 258-1 266.

    [4] Taniyasu Y,Kasu M,Makimoto T.An aluminium nitride light-emitting diode with a wavelength of 210 nanometres[J].Nature,2006,441(7091):325-328.

    [5] Nakamura S.The roles of structural imperfections in InGaN-based blue light-emitting diodes and laser diodes[J].Science,1998,281(5379):956-961.

    [6] 孔月嬋,鄭有炓.Ⅲ族氮化物異質(zhì)結(jié)構(gòu)二維電子氣研究進(jìn)展[J].物理學(xué)進(jìn)展,2006,26(2):127-145.

    [7] 薛麗君,劉 明,王 燕,等.AlGaN/GaN異質(zhì)結(jié)極化行為與二維電子氣[J].半導(dǎo)體技術(shù),2004,29(7):63-65.

    [8] 申 曄,邢懷中,俞建國,等.極化誘導(dǎo)的內(nèi)建電場對Mnδ摻雜的GaN/AlGaN量子阱居里溫度的調(diào)制[J].物理學(xué)報(bào),2007,56(6):3 453-3 457.

    [9] Kozodoy P,Smorchkova Y P,Hansen M,et al.Polarization-enhanced Mg doping of AlGaN/GaN superlattices[J].Applied Physics Letters,1999,75(16):2 444-2 446.

    [10] Allen M W,Miller P,Reeves R J,et al.Influence of spontaneous polarization on the electrical and optical properties of bulk,single crystal ZnO[J].Applied Physics Letters,2007,90(6):301-373.

    [11] 馬繼昭.極化效應(yīng)對紫外發(fā)光二極管光電特性影響的研究[D].南京:南京大學(xué),2013.

    [12] 董可秀.Ⅲ族氮化物的極化效應(yīng)及其在光電子器件中的應(yīng)用[D].南京:南京大學(xué),2013.

    [13] Marzari N,Vanderbilt D.Maximally localized generalized Wannier functions for composite energy bands[J].Physical Review B,1997,56(20):12 847-12 865.

    [14] Souza I,Marzari N,Vanderbilt D.Maximally localized Wannier functions for entangled energy bands[J].Physical Review B,2002,65(3):321-325.

    [15] Vanderbilt D,King-smith R D.Electric polarization as a bulk quantity and its relation to surface charge[J].Physical Review B Condensed Matter,1993,48(7):4 442-4 455.

    [16] Giannozzi P,Baroni S,Bonini N,et al.Quantum Espresso:A modular and open-source software project for quantum simulations of materials[J].Journal of Physics Condensed Matter An Institute of Physics Journal,2009,21(39):1-19.

    [17] Mostofi A A,Yates J R,Lee Y S,et al.wannier90:A tool for obtaining maximally-localised Wannier functions[J].Computer Physics Communications,2008,178(9):685-699.

    [18] Mostofi A A,Yates J R,Pizzi G,et al.An updated version of wannier90:A tool for obtaining maximallylocalised Wannier functions[J].Computer Physics Communications,2014,185(8):2 309-2 310.

    [19] Bernardini F,Fiorentini V,Vanderbilt D.Spontaneous polarization and piezoelectric constants of III-V nitrides[J].1997,56(16):10 024-10 027.

    [20] Gopal P,Spaldin N A.Polarization,piezoelectric constants,and elastic constants of ZnO,MgO,and CdO[J].Journal of Electronic Materials,2006,35(4):538-542.

    [21] Laehnemann J,Brandt O,Jahn U,et al.Direct experimental determination of the spontaneous polarization of GaN[J].Physical Review B,2012,86(8):1-5.

    【責(zé)任編輯:陳佳】

    First-principlesstudyofspontaneouspolarizationofwurtziteAlN,GaNandZnO

    NIU Hai-bo1, ZHU You-zhang1, LI Guan-qiang2

    (1.Department of Physics, Xi′an Jiaotong University City College, Xi′an 710018, China; 2. School of Arts and Sciences, Shaanxi University of Science amp; Technology, Xi′an 710021, China)

    By using Berry phase and maximally localized Wannier functions methods based on the modern polarization theory,the spontaneous polarization of wurtzite AlN,GaN and ZnO are calculated respectively.The origin of spontaneous polarization is studied from the polarization variation of ionic and electronic part.It is found that the spontaneous polarization of AlN is primarily from the electronic part,while the spontaneous polarization of GaN and ZnO is mainly from the ionic part.Since the computational model and reference model have a same volume in the literature,it is found that the result is reduced compared with our result.On the basis of the discussion of local dipole,a correction formula is given to make the calculation more reasonable.Using the charge depth of the formula unit,we also proposed a new method to directly calculate the spontaneous polarization of wurtzite structure without the computation of reference structure.The production of the spontaneous polarization can be intuitively explanation in terms of this method.

    wurtzite structure; spontaneous polarization; Wannier center; first principles

    2017-08-15

    陜西省教育廳專項(xiàng)科研計(jì)劃項(xiàng)目(16JK2099)

    牛海波(1981-),男,山西運(yùn)城人,副教授,博士,研究方向:寬禁帶半導(dǎo)體中的極化場性質(zhì)

    2096-398X(2017)06-0171-08

    O469

    A

    猜你喜歡
    偶極矩鋅礦閃鋅礦
    偶極矩及其排列構(gòu)型
    物理與工程(2024年4期)2024-01-01 00:00:00
    氨基三亞甲基膦酸在閃鋅礦和方鉛礦浮選分離中的應(yīng)用
    金屬礦山(2023年8期)2023-09-19 00:41:10
    西藏甲瑪斑巖成礦系統(tǒng)閃鋅礦礦物學(xué)特征及其地質(zhì)意義*
    鈣(鎂)離子在菱鋅礦表面吸附的量子化學(xué)研究
    對稱和不對稱分子諧波輻射與其結(jié)構(gòu)的內(nèi)在關(guān)系
    青海北祁連陰凹槽塞浦路斯型銅鋅礦特征及找礦標(biāo)志
    Cu-X(X=C,Si,Ge,Sn,Pb)摻雜對閃鋅礦ZnS 可見光吸收的影響研究
    電子是什么形狀?
    貴州五指山特大型鉛鋅礦床閃鋅礦的Rb-Sr定年及其地質(zhì)意義
    澳大利亞杜加爾河鋅礦實(shí)現(xiàn)商業(yè)化生產(chǎn)
    精品国产超薄肉色丝袜足j| 一本色道久久久久久精品综合| 久久精品国产亚洲av涩爱| tube8黄色片| 999精品在线视频| av女优亚洲男人天堂| 啦啦啦啦在线视频资源| 亚洲精品日本国产第一区| 丝袜脚勾引网站| 性色av一级| 伊人久久国产一区二区| 午夜免费观看性视频| 老司机靠b影院| 欧美精品亚洲一区二区| 青春草国产在线视频| 亚洲伊人色综图| 另类亚洲欧美激情| 欧美日韩一区二区视频在线观看视频在线| 久久99精品国语久久久| 国产精品秋霞免费鲁丝片| 久久久久久久精品精品| 午夜激情av网站| 亚洲国产欧美日韩在线播放| 男女无遮挡免费网站观看| 欧美人与善性xxx| 久久精品人人爽人人爽视色| 日本一区二区免费在线视频| 精品一区二区三区av网在线观看 | 成人漫画全彩无遮挡| 黑人欧美特级aaaaaa片| 国产精品久久久久成人av| 另类精品久久| h视频一区二区三区| 最近中文字幕高清免费大全6| 最近手机中文字幕大全| 久久国产精品男人的天堂亚洲| 免费在线观看完整版高清| 国产成人精品在线电影| 亚洲熟女精品中文字幕| 美女扒开内裤让男人捅视频| 免费在线观看视频国产中文字幕亚洲 | 中文字幕制服av| 1024香蕉在线观看| 国产精品 国内视频| 欧美日韩成人在线一区二区| 久久久久视频综合| 少妇猛男粗大的猛烈进出视频| 久热爱精品视频在线9| 日韩精品免费视频一区二区三区| 性高湖久久久久久久久免费观看| 国产精品成人在线| 丰满乱子伦码专区| 老司机影院成人| 天天躁狠狠躁夜夜躁狠狠躁| 性高湖久久久久久久久免费观看| 又大又爽又粗| 熟女少妇亚洲综合色aaa.| 久久综合国产亚洲精品| 中文字幕制服av| 亚洲精华国产精华液的使用体验| 国产成人精品无人区| 观看av在线不卡| 国产亚洲欧美精品永久| 欧美日韩av久久| 国产一级毛片在线| 国产成人欧美| 亚洲一区中文字幕在线| 天堂8中文在线网| 欧美精品av麻豆av| 午夜91福利影院| 涩涩av久久男人的天堂| 亚洲中文av在线| 色精品久久人妻99蜜桃| www.精华液| 国产av国产精品国产| 一边摸一边抽搐一进一出视频| 久久久久久久大尺度免费视频| 欧美黑人精品巨大| 国产男女超爽视频在线观看| 黄片无遮挡物在线观看| av网站在线播放免费| 老司机在亚洲福利影院| 国产乱人偷精品视频| 尾随美女入室| 亚洲av国产av综合av卡| 亚洲国产av影院在线观看| 观看美女的网站| 青青草视频在线视频观看| 国产又爽黄色视频| 亚洲精品av麻豆狂野| xxx大片免费视频| 99精品久久久久人妻精品| 亚洲av中文av极速乱| 日本欧美视频一区| 国产成人欧美在线观看 | 欧美日韩一级在线毛片| 在线精品无人区一区二区三| 国产色婷婷99| 精品亚洲乱码少妇综合久久| 欧美成人精品欧美一级黄| 下体分泌物呈黄色| av不卡在线播放| 日日撸夜夜添| 操出白浆在线播放| 免费观看a级毛片全部| 国产成人精品久久久久久| 一本一本久久a久久精品综合妖精| 一区二区日韩欧美中文字幕| 秋霞伦理黄片| 欧美黑人精品巨大| 亚洲精品一区蜜桃| 国产乱人偷精品视频| 精品亚洲乱码少妇综合久久| 青春草视频在线免费观看| 黄片小视频在线播放| 国产伦理片在线播放av一区| 亚洲av成人不卡在线观看播放网 | 女人高潮潮喷娇喘18禁视频| 91国产中文字幕| 亚洲国产日韩一区二区| 亚洲国产精品999| 国产亚洲最大av| 亚洲一级一片aⅴ在线观看| 国产一区二区在线观看av| 又大又黄又爽视频免费| 亚洲五月色婷婷综合| 最近中文字幕高清免费大全6| 亚洲av男天堂| 国产高清不卡午夜福利| 满18在线观看网站| 狂野欧美激情性xxxx| 热re99久久国产66热| 国产无遮挡羞羞视频在线观看| 国产淫语在线视频| 这个男人来自地球电影免费观看 | 大片免费播放器 马上看| 一区二区三区激情视频| 久久毛片免费看一区二区三区| 丝袜人妻中文字幕| 精品一区在线观看国产| 国产亚洲一区二区精品| 亚洲国产精品999| 1024视频免费在线观看| 午夜激情久久久久久久| 国产乱来视频区| 精品人妻在线不人妻| 午夜福利网站1000一区二区三区| 各种免费的搞黄视频| 亚洲国产最新在线播放| 无遮挡黄片免费观看| 国产免费现黄频在线看| 别揉我奶头~嗯~啊~动态视频 | 一本一本久久a久久精品综合妖精| 国产亚洲av片在线观看秒播厂| 久久久久久人人人人人| 热re99久久精品国产66热6| 18禁动态无遮挡网站| 亚洲欧美清纯卡通| xxx大片免费视频| 伦理电影免费视频| 国产精品免费大片| 一本一本久久a久久精品综合妖精| 男女无遮挡免费网站观看| 男女无遮挡免费网站观看| 少妇被粗大的猛进出69影院| 免费少妇av软件| 国产又爽黄色视频| 久久婷婷青草| 人妻一区二区av| 国产不卡av网站在线观看| 精品亚洲成a人片在线观看| 国产成人欧美| 嫩草影视91久久| 精品亚洲乱码少妇综合久久| 久久精品国产a三级三级三级| 无遮挡黄片免费观看| 欧美变态另类bdsm刘玥| 精品国产一区二区三区久久久樱花| 美女主播在线视频| netflix在线观看网站| 免费黄网站久久成人精品| 波野结衣二区三区在线| av又黄又爽大尺度在线免费看| 免费在线观看视频国产中文字幕亚洲 | 另类亚洲欧美激情| 久久性视频一级片| 丰满乱子伦码专区| 一级黄片播放器| 欧美在线黄色| 纵有疾风起免费观看全集完整版| 在线观看三级黄色| 9热在线视频观看99| 欧美精品高潮呻吟av久久| 人体艺术视频欧美日本| 精品国产一区二区久久| 亚洲欧美一区二区三区久久| 国产日韩欧美亚洲二区| 国产探花极品一区二区| 一区二区日韩欧美中文字幕| 亚洲精品乱久久久久久| 亚洲国产精品一区二区三区在线| 久久精品久久精品一区二区三区| 免费在线观看完整版高清| 国产欧美日韩一区二区三区在线| 亚洲欧美一区二区三区国产| 777久久人妻少妇嫩草av网站| 97人妻天天添夜夜摸| 另类精品久久| 国产三级在线视频| 中文字幕人成人乱码亚洲影| 久久久精品欧美日韩精品| 老司机午夜十八禁免费视频| 岛国在线观看网站| 91精品三级在线观看| 麻豆国产av国片精品| 欧美最黄视频在线播放免费| 女人高潮潮喷娇喘18禁视频| 欧美老熟妇乱子伦牲交| 色av中文字幕| 欧美日韩黄片免| 久久久水蜜桃国产精品网| 久99久视频精品免费| 日日爽夜夜爽网站| 成人av一区二区三区在线看| 亚洲av成人av| 日韩欧美国产在线观看| 亚洲专区中文字幕在线| 美女高潮到喷水免费观看| 在线观看免费午夜福利视频| 熟妇人妻久久中文字幕3abv| 欧美+亚洲+日韩+国产| 国产欧美日韩一区二区三| 女同久久另类99精品国产91| 国产极品粉嫩免费观看在线| 欧美久久黑人一区二区| 欧美性长视频在线观看| 亚洲成人久久性| 国产精品免费一区二区三区在线| 成在线人永久免费视频| 一级a爱片免费观看的视频| 国产xxxxx性猛交| 搡老岳熟女国产| 最近最新免费中文字幕在线| 操美女的视频在线观看| 久久国产精品人妻蜜桃| 777久久人妻少妇嫩草av网站| 村上凉子中文字幕在线| 亚洲人成77777在线视频| 久9热在线精品视频| 法律面前人人平等表现在哪些方面| 一本综合久久免费| 91字幕亚洲| 88av欧美| 又大又爽又粗| 成年女人毛片免费观看观看9| 两个人看的免费小视频| 国语自产精品视频在线第100页| 最好的美女福利视频网| 久久精品国产清高在天天线| 午夜亚洲福利在线播放| 国产极品粉嫩免费观看在线| 日日夜夜操网爽| 精品人妻1区二区| 欧美成人性av电影在线观看| 久久久久精品国产欧美久久久| 悠悠久久av| 精品久久久久久久毛片微露脸| 别揉我奶头~嗯~啊~动态视频| 成人永久免费在线观看视频| 国产私拍福利视频在线观看| bbb黄色大片| 欧美绝顶高潮抽搐喷水| 丰满的人妻完整版| bbb黄色大片| 天天一区二区日本电影三级 | 亚洲激情在线av| 午夜精品久久久久久毛片777| 窝窝影院91人妻| 午夜影院日韩av| 亚洲aⅴ乱码一区二区在线播放 | 少妇被粗大的猛进出69影院| 这个男人来自地球电影免费观看| 中文字幕高清在线视频| 在线观看一区二区三区| 女人精品久久久久毛片| 999精品在线视频| 日韩欧美三级三区| 国产精品久久电影中文字幕| 18禁观看日本| 亚洲av熟女| 在线国产一区二区在线| 美女午夜性视频免费| 丝袜人妻中文字幕| 亚洲自偷自拍图片 自拍| 欧美黑人精品巨大| 精品欧美一区二区三区在线| 欧美日韩黄片免| 大陆偷拍与自拍| 日本三级黄在线观看| 亚洲色图av天堂| 久久亚洲精品不卡| 国产成人啪精品午夜网站| 亚洲国产高清在线一区二区三 | 纯流量卡能插随身wifi吗| 俄罗斯特黄特色一大片| 久久伊人香网站| 在线观看www视频免费| 久久久国产精品麻豆| 国产成人欧美在线观看| 成人特级黄色片久久久久久久| 成人免费观看视频高清| 久久伊人香网站| 一进一出抽搐gif免费好疼| 国内毛片毛片毛片毛片毛片| 久久久国产成人精品二区| 欧美成人午夜精品| 国产精品自产拍在线观看55亚洲| www.精华液| 亚洲精品在线美女| 亚洲av美国av| 亚洲成a人片在线一区二区| 少妇 在线观看| 国产一区二区激情短视频| 中出人妻视频一区二区| 精品午夜福利视频在线观看一区| 日韩成人在线观看一区二区三区| 一区在线观看完整版| 91九色精品人成在线观看| 精品欧美国产一区二区三| 国产成人精品在线电影| 成人手机av| 免费一级毛片在线播放高清视频 | 91九色精品人成在线观看| 少妇 在线观看| 在线视频色国产色| 美女高潮喷水抽搐中文字幕| 国产成人影院久久av| √禁漫天堂资源中文www| 久久精品亚洲精品国产色婷小说| 黑人巨大精品欧美一区二区蜜桃| 十八禁网站免费在线| www.www免费av| 脱女人内裤的视频| 在线视频色国产色| 一本久久中文字幕| 欧美日本视频| 亚洲第一青青草原| 亚洲成a人片在线一区二区| 黄网站色视频无遮挡免费观看| 啦啦啦 在线观看视频| 国产麻豆69| 欧美不卡视频在线免费观看 | 一边摸一边抽搐一进一小说| 国产成人精品久久二区二区免费| 久久亚洲精品不卡| 精品无人区乱码1区二区| 女性被躁到高潮视频| 亚洲成av片中文字幕在线观看| 性少妇av在线| 两个人视频免费观看高清| 国产亚洲欧美在线一区二区| 免费看美女性在线毛片视频| 亚洲av熟女| 亚洲成a人片在线一区二区| 一区二区三区国产精品乱码| 国语自产精品视频在线第100页| 精品人妻1区二区| 美女高潮到喷水免费观看| 桃红色精品国产亚洲av| svipshipincom国产片| 9热在线视频观看99| 在线观看免费视频网站a站| 18禁美女被吸乳视频| aaaaa片日本免费| 亚洲天堂国产精品一区在线| 久热爱精品视频在线9| 国产精品,欧美在线| 国产亚洲av嫩草精品影院| 国产精品一区二区在线不卡| 久久久久亚洲av毛片大全| 国产精品 国内视频| 黄色毛片三级朝国网站| 两个人免费观看高清视频| 日本一区二区免费在线视频| 亚洲专区国产一区二区| 91成年电影在线观看| 国产片内射在线| 每晚都被弄得嗷嗷叫到高潮| 免费女性裸体啪啪无遮挡网站| 两个人看的免费小视频| 国产成人系列免费观看| 青草久久国产| 成人国产综合亚洲| 成人亚洲精品一区在线观看| 91大片在线观看| 日韩精品青青久久久久久| 丝袜人妻中文字幕| 色婷婷久久久亚洲欧美| 久久香蕉国产精品| 国产片内射在线| 91成人精品电影| 后天国语完整版免费观看| 夜夜躁狠狠躁天天躁| 99国产极品粉嫩在线观看| 久久欧美精品欧美久久欧美| 久久国产乱子伦精品免费另类| 操美女的视频在线观看| 国产精品亚洲一级av第二区| 制服丝袜大香蕉在线| 操美女的视频在线观看| 久久久久国内视频| 亚洲第一欧美日韩一区二区三区| 精品久久久久久久人妻蜜臀av | 成人亚洲精品av一区二区| 亚洲精品一卡2卡三卡4卡5卡| 亚洲激情在线av| 可以在线观看的亚洲视频| 欧美乱码精品一区二区三区| 妹子高潮喷水视频| 十八禁网站免费在线| 一进一出抽搐动态| 久久久久精品国产欧美久久久| 在线观看66精品国产| 99re在线观看精品视频| 一个人免费在线观看的高清视频| bbb黄色大片| 一二三四在线观看免费中文在| 91成年电影在线观看| 久久国产精品人妻蜜桃| 国产免费av片在线观看野外av| 国产成人av教育| 亚洲 欧美 日韩 在线 免费| 啪啪无遮挡十八禁网站| 女人爽到高潮嗷嗷叫在线视频| 91精品三级在线观看| 国产精品 欧美亚洲| 男女午夜视频在线观看| 欧美日韩一级在线毛片| 久久国产精品影院| 精品免费久久久久久久清纯| 人妻丰满熟妇av一区二区三区| 亚洲片人在线观看| 女人精品久久久久毛片| 国产精品一区二区精品视频观看| 黄色视频,在线免费观看| 精品国产乱子伦一区二区三区| 亚洲国产看品久久| bbb黄色大片| 精品国产亚洲在线| av网站免费在线观看视频| 国产成人精品久久二区二区免费| 波多野结衣一区麻豆| 国产精品99久久99久久久不卡| 黄网站色视频无遮挡免费观看| 一区二区三区激情视频| avwww免费| 国产精品一区二区精品视频观看| av视频免费观看在线观看| 久久草成人影院| 纯流量卡能插随身wifi吗| 久久午夜亚洲精品久久| 一边摸一边做爽爽视频免费| 两个人视频免费观看高清| 国产亚洲精品一区二区www| 免费av毛片视频| 在线观看日韩欧美| 国产熟女xx| 欧美在线黄色| 色精品久久人妻99蜜桃| 美女扒开内裤让男人捅视频| 久久久久久免费高清国产稀缺| 色综合站精品国产| 日本vs欧美在线观看视频| 又黄又爽又免费观看的视频| 亚洲欧美日韩无卡精品| 亚洲狠狠婷婷综合久久图片| bbb黄色大片| 亚洲一区二区三区不卡视频| 不卡一级毛片| 高清黄色对白视频在线免费看| 9色porny在线观看| 亚洲av成人不卡在线观看播放网| 人成视频在线观看免费观看| 制服人妻中文乱码| 人人妻人人澡欧美一区二区 | 久久精品国产亚洲av高清一级| 一卡2卡三卡四卡精品乱码亚洲| 亚洲专区国产一区二区| 天堂√8在线中文| 中文字幕色久视频| 一区二区三区激情视频| videosex国产| 大码成人一级视频| 久久久久久免费高清国产稀缺| 免费观看人在逋| 久久热在线av| 黄片播放在线免费| 伊人久久大香线蕉亚洲五| 国产熟女午夜一区二区三区| 不卡一级毛片| 好男人在线观看高清免费视频 | 丁香欧美五月| 最近最新中文字幕大全免费视频| 国产亚洲欧美精品永久| 日本欧美视频一区| 欧美一级毛片孕妇| 午夜福利高清视频| 亚洲成av人片免费观看| 国产视频一区二区在线看| 狠狠狠狠99中文字幕| 一边摸一边抽搐一进一出视频| 欧美日韩亚洲综合一区二区三区_| 黑人操中国人逼视频| videosex国产| 成人国产综合亚洲| 波多野结衣一区麻豆| 亚洲国产精品sss在线观看| 国产日韩一区二区三区精品不卡| 俄罗斯特黄特色一大片| 久久香蕉精品热| 国产精品秋霞免费鲁丝片| 亚洲精品久久成人aⅴ小说| 日日夜夜操网爽| 高清黄色对白视频在线免费看| 欧美亚洲日本最大视频资源| 婷婷丁香在线五月| 变态另类丝袜制服| 手机成人av网站| 久99久视频精品免费| 日韩精品免费视频一区二区三区| 丁香六月欧美| 日本 欧美在线| 成人国产一区最新在线观看| 国产激情久久老熟女| av有码第一页| av免费在线观看网站| 热99re8久久精品国产| 精品人妻在线不人妻| 中文字幕高清在线视频| 视频在线观看一区二区三区| 亚洲少妇的诱惑av| 欧美久久黑人一区二区| 久久人妻熟女aⅴ| 亚洲精品中文字幕在线视频| 欧美成人免费av一区二区三区| 国产一区二区三区综合在线观看| 激情在线观看视频在线高清| bbb黄色大片| 国产免费av片在线观看野外av| 在线国产一区二区在线| 日韩高清综合在线| 久久久国产欧美日韩av| 久久人人爽av亚洲精品天堂| 久久精品国产亚洲av高清一级| 少妇粗大呻吟视频| √禁漫天堂资源中文www| 露出奶头的视频| 亚洲欧美日韩另类电影网站| 国产97色在线日韩免费| 国产精品一区二区三区四区久久 | 精品少妇一区二区三区视频日本电影| 操美女的视频在线观看| 人人妻人人爽人人添夜夜欢视频| 天天一区二区日本电影三级 | 一级a爱片免费观看的视频| 啦啦啦观看免费观看视频高清 | 欧美亚洲日本最大视频资源| 成人亚洲精品av一区二区| av有码第一页| 在线观看免费视频日本深夜| 国产精品亚洲av一区麻豆| 人人妻人人爽人人添夜夜欢视频| 岛国在线观看网站| 国产成人av教育| 宅男免费午夜| 少妇熟女aⅴ在线视频| 欧美日韩福利视频一区二区| 国产主播在线观看一区二区| 久久精品国产综合久久久| 91字幕亚洲| 男女之事视频高清在线观看| 亚洲国产精品合色在线| 极品教师在线免费播放| 黑人操中国人逼视频| а√天堂www在线а√下载| 中文亚洲av片在线观看爽| 午夜精品久久久久久毛片777| 老熟妇仑乱视频hdxx| 亚洲第一欧美日韩一区二区三区| 国产欧美日韩一区二区三| 国产精品一区二区免费欧美| 国产黄a三级三级三级人| 亚洲专区字幕在线| 别揉我奶头~嗯~啊~动态视频| 久久精品国产亚洲av香蕉五月| 咕卡用的链子| 亚洲成av人片免费观看| 精品一品国产午夜福利视频| 丝袜美足系列| 午夜精品在线福利| 久久精品aⅴ一区二区三区四区| 欧美日本视频| 久9热在线精品视频| 日韩精品免费视频一区二区三区| 99久久99久久久精品蜜桃| 国产av又大| av视频免费观看在线观看| 久久精品亚洲精品国产色婷小说| 99精品久久久久人妻精品| 欧美午夜高清在线| 在线十欧美十亚洲十日本专区| 一进一出抽搐动态| 在线免费观看的www视频| 亚洲最大成人中文| 国产成人av教育|