999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

有機電光晶體4-(4-二甲基氨基苯乙烯基)甲基吡啶對甲基苯磺酸鹽的太赫茲光譜研究?

2018-01-18 19:01:24連宇翔戴澤林許向東谷雨李欣榮王福楊春成曉夢周華新
物理學報 2017年24期
關鍵詞:振動優化

連宇翔 戴澤林 許向東 谷雨 李欣榮 王福楊春 成曉夢 周華新

1)(電子科技大學光電信息學院,電子薄膜與集成器件國家重點實驗室,成都 610054)2)(四川師范大學化學與材料科學學院,成都 610068)

(2017年6月28日收到;2017年7月25日收到修改稿)

1 引 言

太赫茲(THz)波是指頻率從0.1—10 THz的電磁波,介于毫米波與遠紅外之間(3 mm—30μm).THz波具有良好的介質穿透性、低電離能和相干性等優異特性.而且,許多大分子化合物的特征譜也在THz頻段,因此THz技術在物理、化學、生物、天文和醫藥科學等基礎研究領域以及醫學成像、環境監測、安全檢查等應用研究領域均具有巨大的科研價值和廣闊的應用前景[1].近年來,利用有機電光晶體的光整流效應進行太赫茲波的產生和探測已成為太赫茲光子學領域重要的研究課題.與傳統的無機電光晶體相比,性能優良的有機電光晶體更易設計和合成.為了滿足太赫茲源和太赫茲探測的需要,大量具有較高的二階非線性系數的新的有機電光晶體已相繼被合成出來,其中最具代表性的是一種有機吡啶磺酸鹽,即4-(4-二甲基氨基苯乙烯基)甲基吡啶對甲基苯磺酸鹽(DAST).Zhang等[2]開展的DAST的光整流產生太赫茲輻射的實驗表明,DAST的整流電場強度是無機電光晶體LiTaO3的185倍、GaAs和InP晶體的42倍.這主要得益于DAST較高的二階非線性系數(在1318 nm波長激光作用下,DAST的二階非線性系數d111=(1010±110)pm/V)和電光系數(在1535 nm波長激光作用下,DAST的電光系數r111=(47±8)pm/V)[3].然而,DAST在太赫茲波段還表現出強烈的吸收和色散[4,5],這對獲得連續寬帶的太赫茲光源極為不利.強烈的吸收帶會嚴重地降低太赫茲產生效率,而強烈色散也會使得相位匹配不容易實現.類似的情況也在其他有機電光分子晶體如PNP,OH1[6]中出現.因此,研究有機電光晶體在太赫茲波段的吸收和色散性質,對于理解其太赫茲波產生和探測的特點以及優化和設計高效的用于太赫茲源和探測的新材料等,具有重要的意義.

Walther等[4]對DAST單晶在0—3 THz波段a軸和b軸的吸收和折射率色散性質進行了系統的測量,結合簡單的橫向光學聲子計算,認為DAST在此波段的振動主要是源于陰陽離子間的光學聲子振動.然而,Walther等并未對這些吸收峰對應的振動方式進行具體的歸屬和驗證.Glavcheva等[5]測試了DAST粉末在0.63—10.6 THz的吸收光譜,但缺乏深入的分析.顯然,詳細的模式分析不僅有助于理解聲子振動的起源,為今后合成性能更優的新的DAST衍生物以及其他有機電光晶體提供有價值的指導,同時也有助于理解DAST作為太赫茲發射和探測的特點及注意環節.遺憾的是,到目前為止,此項研究在國內外文獻中十分罕見.雖然Saito等[7]使用第一性原理對DAST進行了理論計算,然而其計算結果與實驗測量的結果出入較大,這可能是由于計算參數和初始結構不合理.而且,作者僅對103.2 cm?1處的吸收峰做了模式分析[7],缺乏對低頻的聲子振動進行詳細的討論.

詳細地分析和歸屬太赫茲振動譜帶有利于更加深入地了解低頻振動動力學的物理機理,這一任務可通過基于密度泛函理論(DFT)的量子力學計算實現.傳統的中紅外光譜振動模式主要取決于分子內的共價鍵作用.與之不同,太赫茲光譜與分子間的相互作用,例如氫鍵和范德瓦耳斯力密切相關.因此,理論計算的一個關鍵問題是如何合理地描述這些非共價鍵作用.最近,我們通過對簡單磺酸類有機物牛磺酸進行DFT光譜分析發現[8],對于氫鍵分子晶體,色散校正在模擬分子在晶體中所處的環境時十分重要[9,10].與簡單的牛磺酸小分子相比較,DAST的分子量更大、結構更加復雜,因此需要通過引入色散校正的方法對DAST進行結構優化及太赫茲光譜計算.為此,本文首次采用周期性固態色散校正的密度泛函理論(dispersion-corrected density functional theory,DFT-D2),對有機電光晶體DAST在0—4 THz的聲子振動模式進行詳細的分析和歸屬.在收斂測試的基礎上,使用逐步提高收斂精度優化法進行紅外光譜計算.根據DFT計算結果,探究DAST主要振動模式的起源.本文不僅探索出使用一個從頭算量子力學程序Cambridge sequential total energy package(CASTEP)預測復雜分子紅外光譜的合適方法,還揭示了DFT計算在太赫茲光子學中的重要應用價值.對探究有機電光晶體的太赫茲響應物理原理具有重要的指導價值.

2 理論計算

2.1 DFT計算

DFT是用來研究多電子體系電子結構的一種量子力學的方法,在物理、化學、材料等領域都有廣泛應用.研究表明單個分子的DFT計算,不能很好地模擬范德瓦耳斯力在分子間的作用,導致模擬計算結果與實驗有一定差別[11,12].為了解決這一問題,采用Grimme[13]的方法將范德瓦耳斯作用中的色散力校正項引入到DFT計算中(DFT-D2法)[8].在該方法中,體系中各原子對之間的范德瓦耳斯相互作用通過一個簡單的成對力場來描述.本文使用CASTEP對有機電光晶體DAST進行結構優化和太赫茲光譜計算.與液相或氣相的獨立分子建模計算相比較,CASTEP考慮了固體的周期性晶格排布,更為準確地描述了分子吸收光譜以及集體振動模式.同時,在太赫茲波段的分子振動模式識別中,尤其是研究分子間或分子內作用力時,周期性DFT計算已經被廣泛采用[14,15].據此,本文采用周期性固態DFT-D2方法進行DAST的太赫茲光譜計算.

2.2 計算參數設置及優化

選取Marder等建立的DAST晶格參數和原子坐標作為計算的初始結構[16].在幾何優化之前,分別對截斷能和k點網格進行測試,以便選取合適的計算參數用于結構優化和光譜計算.所有周期性DFT計算均在CASTEP[17]上完成.優化后的平面波截斷能設為800 eV,選用廣義梯度[18]近似下的Perdew-Burke-Ernzerhof(PBE)泛函作為交換關聯函數的近似.幾何優化和光譜計算均采用模守恒贗勢和Monkhorst-Packk點網格(5×5×5).使用Grimme的PBE-D2方法[19]將色散校正項引入PBE泛函.該方法中庫侖作用通過成對的力場加以描述.在PBE-D2方法的計算中,晶格常數和原子坐標均完全弛豫.幾何優化的總能量收斂標準均設為10?6Hartree.計算得到的光譜采用5 cm?1半高全寬的洛倫茲線型繪制.

3 結果與討論

有機電光晶體DAST作為太赫茲發射源,主要利用其較高的二階非線性系數.DAST分子由一個帶正電的有機吡啶鹽和一個帶負電的磺酸鹽組成,如圖1(a)所示.其中,有機吡啶鹽由電子給體(二甲基氨基)、共軛π橋(苯乙烯基)和電子受體(甲基吡啶)組成.帶正電荷的甲基吡啶具有很強的吸電子能力,并與二甲基氨基通過苯乙烯基形成推-拉電子結構,因此該結構具有很大的一階超極化率.DAST的有機吡啶鹽鏈可看作一個偶極子(如圖1(a)所示),其極化率p可表示為

其中,q為電荷電量,l為電荷間位移矢量,ε0為真空介電常數,E為外加電場;α,β和γ分別為分子的一階、二階和三階非線性系數.在DAST晶體中,磺酸鹽(負)與有機吡啶鹽(正)電荷間的庫侖作用大于有機吡啶鹽間的偶極子-偶極子作用,使得有機吡啶鹽沿著與a軸夾角約為20°的方向排列(如圖1(d)所示).而且,整個DAST在三維空間中以非中心對稱的形式排列,以獲得最大的宏觀極化率.宏觀極化率P可表示為

其中,χ(1),χ(2)和χ(3)分別為一階、二階和三階非線性系數.χ(2)與β成正比關系[20],即分子的超極化率越大,晶體的二階非線性系數越大.

圖1 (網刊彩色)DAST分子結構及其原胞沿各軸方向上的分子排列 (a)DAST的分子結構(由右側一個帶正電的有機吡啶鹽和左側一個帶負電的磺酸鹽組成);(b)DAST原胞沿a軸(晶格夾角A)的分子排列;(c)DAST原胞沿b軸(晶格夾角B)的分子排列;(d)DAST原胞沿c軸(晶格夾角Γ)的分子排列Fig.1.(color online)DAST molecular structure and the molecular arrangement of the cells along the axis:(a)DAST molecular structure consisting of a positively-charged organic pyridine salt and a negatively-charged sulfonate;(b)the molecular arrangement of the cells along a axis(lattice angle A);(c)the molecular arrangement of the cells along b axis(lattice angle B);(d)the molecular arrangement of the cells along c axis(lattice angle Γ).

從以上討論可以看出,DAST之所以具有非常高的二階非線性系數,主要是源于有機吡啶鹽較大的分子極化率及其非中心對稱的分子排列方式.磺酸鹽的存在主要是為了使有機吡啶鹽能夠在晶體內有序地定向排列.這種成鹽的方式是設計和合成二階非線性材料的重要方法[16].陰陽離子間的庫侖作用為DAST穩定的非中心對稱結構提供了保障.雖然DAST極高的二階非線性系數使得其通過光整流技術在太赫茲源的應用上具有優勢,但是DAST的這種結構也使其在太赫茲波段存在本征吸收帶,影響太赫茲波發射效率.下面重點分析DAST在太赫茲波段強的聲子吸收帶的起源.

在進行幾何優化時,發現直接對初始結構進行高精度計算得到的最終結構與實驗值差別較大(優化后的晶格常數與實驗值相差超過10%),這可能是由于實驗所測的DAST初始結構不太合理,與基態相差較大.為此,本文采用逐步增加收斂精度的方法.首先,進行低精度的計算,然后逐步過渡到高精度.圖2顯示了從低精度到高精度(分別對應CASTEP幾何優化中的coarse,medium, fi ne和ultra fi ne)優化過程中晶格常數a,b,c和晶格夾角A,B,Γ(圖2)的變化.

從圖2可以看出,在逐步提高收斂精度優化過程中,晶格常數、晶格夾角與實驗值相比波動很小(小于1%),而且優化過程中晶格對稱性保持不變,晶格常數始終保持a=b,晶格夾角A=B.這表明,逐步優化的方法是可行的.為了方便比較,表1列出了ultra fi ne后的晶格常數與晶格夾角同初始結構實驗數據的對比.圖2中,隨著收斂精度的提高,晶格常數a,b,c的變化較小(Δa= Δb=0.0095 ?,0.0012%;Δc=0.1154 ?,0.0065%);晶格夾角A,B,Γ雖然在medium和 fi ne相對初始結構出現了超過1°的變化,但是最終的夾角參數(ultra fi ne)與初始實測結構相比,變化也很小(ΔA= ΔB=05012°,0.0055%;ΔΓ=0.089°,0.0009%).幾何優化本質上是尋找體系的最小能量狀態,低精度的優化有助于將優化的方向逐步推向正確的收斂方向.逐步優化的方法對于大體系的結構優化十分有效,既能節省時間和計算成本,又能找到體系的基態.優化結果表明,采用逐步收斂的方法,能夠較好地得到DAST的基態穩定結構,使得理論計算得到的太赫茲光譜與實測結果更好地相符.

圖2 (網刊彩色)逐步優化過程中DAST晶格參數變化(晶格常數a,b,c;晶格夾角A,B,Γ;對應收斂精度從低精度initial到高精度ultra fi ne)Fig.2.(color online)Changes of the DAST lattice parameters during gradual optimization(lattice constants of a,b and c;lattice angles of A,B and Γ;convergence accuracy from initial to ultra fi ne).

表1 DAST初始結構實驗數據同結構收斂優化后對比(晶格常數a,b,c;晶格夾角A,B,Γ)Table 1.Comparison of the experimental data for the initial structure of DAST and those obtained by structure optimization(lattice constants of a,b and c;lattice angles of A,B and Γ).

最近,我們采用DFT-D2方法,對牛磺酸小分子的太赫茲光譜進行計算,得到與實驗較為相符的仿真結果[8].據此,本文采用DFT-D2方法,對復雜程度更高、分子量更大的DAST在0—4 THz波段的太赫茲光譜進行了計算,結果如圖3所示.作為比較,圖3還給出了Glavcheva等[5]的DAST多晶粉末的太赫茲光譜測量結果.圖3顯示,DFT計算很好地重現了實測的5個吸收峰M1—M5的位置,這為進一步分析DAST在這一波段的振動起源奠定了基礎.注意到雖然Saito等[7]基于第一性原理的光譜計算也重現出了M1和M5,與本文的計算結果類似.但是,Saito等[7]的計算沒有很好地重現M2,M3,M4三個振動模式,而且其計算得到的峰形與實驗結果差別較大.與之不同,本文對M1—M5振動模式都很好地進行了重現.據我們所知,本文是第一次通過DFT理論計算的方法完全預測出DAST在0—4 THz范圍的所有振動模式.而且,在本文仿真過程中,紅外光譜計算沒有虛頻產生,說明采用逐步優化的方法得到的基態結構是穩定的.同時說明采用DFT-D2方法,仿真獲得的DAST的太赫茲光譜與實測數據的符合程度高,說明色散校正的引入確實可以較好地模擬分子在晶體中的環境.

圖3 (網刊彩色)在0—4 THz波段,通過DFT仿真計算的DAST太赫茲光譜與實驗測量的DAST粉末太赫茲光譜(文獻[5])的比較Fig.3.(color online)Comparison of the DAST THz spectrum in the frequency band of 0–4 THz simulated by DFT and that for the DAST powder experimentally measured(Ref.[5])

為便于比較,表2列出了理論計算和實測吸收峰的位置,它們之間的平均誤差僅為0.088 THz(約2.9 cm?1).值得注意的是,M1和M3兩個較強的吸收帶也在以往DAST單晶偏振測試中出現[5],是DAST在太赫茲低頻波段的主要吸收峰.

表2 IR計算的頻率位置與實測結果的比較Table 2.Comparison of the frequency positions obtained by IR calculation and the experimentallymeasured results.

為了探究DAST的物理特性,本文繼續對DAST在太赫茲波段的聲子模式歸屬進行研究.圖4顯示了振動模式M1—M5的振動位移矢量圖.下面結合CASTEP的動畫模擬功能,討論DAST實測的太赫茲吸收峰的起源.單分子振轉模式表現分子內各基團的振動和轉動,而晶胞分子的振轉模式則比較復雜,表現為晶胞各分子的對稱和非對稱彎曲的協同振動以及分子內基團振動.在1.12 THz(對應計算的M1)處,DAST陽離子(有機吡啶鹽)和陰離子(磺酸鹽)分別在各自的(苯環)平面內發生平移振動.陰陽離子的振動方向如圖中紅色箭頭所示.該模式可歸屬為陰陽離子的相對運動,即光學聲子模式.從圖4可以看出,陰陽離子的相對振動方向基本沿著a軸,這解釋了當使用沿著a軸偏振的太赫茲入射波進行光譜測試時,DAST單晶在M1處的吸收峰強度高于沿b軸偏振測試的吸收強度的現象[4].1.46 THz(對應計算的M2)和1.54 THz(對應計算的M3)主要源于有機磺酸鹽的振動.其中,1.46 THz的振動是由于磺酸鹽沿a軸的轉動,而1.54 THz則是由于整個磺酸離子沿著c軸的平動.這兩處的吸收峰與陽離子基本無關.在文獻[5]中,Glavcheva等[5]還比較測量了DAST和DASC的太赫茲光譜.DASC的結晶類型和化學結構與DAST相似,不同之處是磺酸鹽上的—CH3被—Cl取代.Glavcheva等[5]的實測結果顯示,兩者的低頻振動模式M1—M5位置完全相同.這暗示著DAST的振動模式M1—M5與單個—CH3官能團無關,而更有可能與陰陽離子間的集體振動有關.這不僅體現了太赫茲振動模式本質上的復雜性,也有助于人們更加直觀地認識到中紅外光譜與太赫茲光譜的不同:中紅外光譜對于物質的化學成分(化學鍵)更加敏感,而太赫茲光譜則對于物質的結構排列更為敏感.由于有機電光晶體DAST二階非線性的貢獻主要來源于有機吡啶鹽陽離子大的分子極化率及其非中心對稱的分子排列,這為設計具有二階非線性的有機電光晶體DAST新型衍生物提供了一個方向,即通過設計或改變陰離子基團來消除DAST在1.46 THz(對應計算的M2)和1.54 THz(對應計算的M3)處的吸收,從而有望提高太赫茲波的產生效率.仔細觀察還發現,M2和M3的振動中陰離子的苯環振動較強,因此可以通過改變陰離子的結構,例如將苯環替代成烷烴或其他官能團,或將苯環的H原子用其他官能團取代,牽制苯環的振動,據此達到減弱或消除M2和M3處吸收的目的.反之,如果需要利用DAST在太赫茲波段的吸收特性,則可以考慮通過取代陰陽離子基團來增強DAST在太赫茲波段的吸收.2.63 THz(對應M4)處主要是有機吡啶鹽兩個苯環的扭轉振動以及陰離子繞a軸的轉動.而3.16 THz處最強的吸收峰(對應M5)與M4的振動類型類似,但振動的方向不同(如圖中紅色箭頭所指示).綜合這幾處的振動模式可以發現,有機磺酸鹽陰離子在DAST的低頻振動中起了關鍵的作用.本文首次對DAST在0—4 THz范圍太赫茲光譜的特征峰進行了系統歸屬.

圖4 (網刊彩色)通過CASTEP動畫模擬功能得到的DAST分子振動模式M1—M5對應的振動位移矢量圖,其中,M1源于陽離子(有機吡啶鹽)、陰離子(磺酸鹽)分別在各自的(苯環)平面內發生平移振動;M2,M3分別源于陰離子沿a軸、c軸轉動;M4源于陽離子兩個苯環扭轉振動、陰離子繞a軸的轉動;M5源于陰陽離子苯環扭轉轉動Fig.4.(color online)Vibration displacement vector diagrams for the DAST molecular structure obtained by CASTEP animation simulation function.M1 is originated from the DAST cation(organic pyridinium salt)and anion(sulfonate)undergo translational vibrations in their respective(benzene ring)plane;M2 and M3 are assigned to the anion rotates along the a axis and c axis,respectively;M4 is related to the torsional vibrations of two benzene rings on cation and the anion rotates around the a axis;M5 is attributed to the torsional vibrations of the benzene rings in the cation and anion.

4 結 論

首次將DFT-D2用于計算結構復雜的有機磺酸鹽DAST的太赫茲光譜.在收斂測試(k點和截斷能測試)的基礎上,采用逐步提高精度進行幾何優化的方法尋找DAST收斂的基態穩定結構.結果表明,該方法不僅很好地重現了DAST的結構參數(例如晶格常數和夾角),而且在此結構基礎上計算的太赫茲光譜與實測光譜較好地相符.在此基礎上,首次對DAST在0—4 THz范圍太赫茲光譜的特征峰進行系統歸屬.模式分析結果表明,1.12 THz處的振動是DAST的陰陽離子的光學聲子模式,1.46 THz和1.54 THz兩處的振動主要與其磺酸鹽相關,而2.63 THz和3.16 THz兩處的振動則源于陽離子的扭轉振動和陰離子的轉動.這說明在不改變DAST二階非線性效應的基礎上,能夠通過設計、改變磺酸鹽陰陽離子基團的方式,減弱DAST在太赫茲波段的吸收,提高材料在太赫茲波段的發射效率.研究有機電光晶體DAST在太赫茲波段的振動模式歸屬不僅有助于全面掌握DAST在太赫茲的吸收特性,對研制性能更加優越的新型太赫茲有機敏感材料,探索DAST基材料在太赫茲探測器中的新應用等都具有重要的意義,還為今后合成具有更優二階非線性的DAST新型衍生物提供了重要的參考和指導.本文結果揭示了DFT在太赫茲光子學中的重要應用價值,為使用DFT計算預測其他用于有機電光晶體的太赫茲吸收性能提供了可資的借鑒方法,有助于深入探究有機電光晶體的太赫茲響應的物理原理.

[1]Ferguson B,Zhang X C 2002Nat.Mater.1 26

[2]Zhang X C,Ma F,Jin Y,Lu T M,Boden E P,Phelps P D,Stewart K R,Yakymyshyn C P 1992Appl.Phys.Lett.61 3080

[3]Bosshard C,Spreiter R,Degiorgi L,Gunter P 2002Phys.Rev.B66 205107

[4]Walther M,Jensby K,Keiding S R 2000Opt.Lett.25 911

[5]Glavcheva Z,Umezawa H,Mineno Y,Odani T,Okada S,Ikeda S,Taniuchi T,Nakanish H 2005Jpn.J.Appl.Phys.44 5231

[6]Kim J,Kwon O P,Brunner F D J,Jazbinsek M,Lee S H 2015J.Phys.Chem.C119 10031

[7]Saito S,Inerbaev T M,Mizuseki H,Igarashi N,Note R,Kawazoe Y 2006Chem.Phys.Lett.432 157

[8]Dai Z L,Xu X D,Gu Y,Li X R,Wang F,Lian Y X,Fan K,Chen X M,Chen Z G,Sun M H,Jiang Y D,Yang C,Xu J 2017J.Chem.Phys.146 124119

[9]Kim J,Kwon O P,Jazbinsek M,Park Y C,Lee Y S 2015J.Phys.Chem.C119 12598

[10]King M D,Buchanan W D,Korter T M 2011Phys.Chem.Chem.Phys.13 4250

[11]Takahashi M 2014Crystals4 74

[12]Grimme S,Ehrlich S,Goerigk L 2011J.Comput.Chem.32 1456

[13]Grimme S 2004J.Comput.Chem.25 1463

[14]Miles R E,Zhang X C,Eisele H,Krotkus A 2007Terahertz Frequency Detection and Identi fi cation of Materials and Objects(Netherlands:Springer)pp147–163

[15]Zhang Y,Peng X H,Chen Y,Chen J,Curioni A,Andreoni W,Nayak S K,Zhang X C 2008Chem.Phys.Lett452 59

[16]Seidler T,Stadnicka K,Champagne B 2014J.Chem.Phys.141 104109

[17]Clark S J,Segall M D,Pickard C J,Hasnip P J,Probert M I J,Refson K,Payne M C 2005Zeitschrift für Kristallographie220 567

[18]Perdew J P,Burke K,Ernzerhof M 1996Phys.Rev.Lett.77 3865

[19]Grimme S 2006J.Comput.Chem.27 1787

[20]Marder S R,Perry J W,Yakymyshyn C P 1994Chem.Mater.6 1137

猜你喜歡
振動優化
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
噴水推進高速艇尾部振動響應分析
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
This “Singing Highway”plays music
振動攪拌 震動創新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
主站蜘蛛池模板: 国产精品制服| 国产青榴视频| lhav亚洲精品| 免费网站成人亚洲| 亚洲欧洲日产国产无码AV| 亚洲人成电影在线播放| 国产一级裸网站| 国产 在线视频无码| 国产成人综合亚洲网址| 成人在线天堂| 国产在线观看第二页| 国产h视频免费观看| 欧美精品成人| 欧美中文一区| 国产激情国语对白普通话| 97成人在线视频| 无码精品一区二区久久久| 女人18毛片一级毛片在线 | 特级aaaaaaaaa毛片免费视频 | 亚洲精品欧美日韩在线| 亚洲人视频在线观看| 91国内视频在线观看| 久久亚洲国产最新网站| 麻豆国产精品视频| 国产不卡网| 午夜日韩久久影院| 四虎国产精品永久一区| 毛片手机在线看| 久草视频精品| 国产在线精品人成导航| 亚洲天堂网视频| 青青青草国产| 国产丝袜91| 全部免费特黄特色大片视频| 国产色图在线观看| 高清无码不卡视频| 国产亚洲精久久久久久久91| 国产精品视频系列专区| 久久这里只精品国产99热8| 欧美精品亚洲精品日韩专区| 9丨情侣偷在线精品国产| 亚洲av成人无码网站在线观看| 亚洲日本一本dvd高清| 六月婷婷精品视频在线观看 | 国产91丝袜在线播放动漫| 国产国产人成免费视频77777| 午夜国产不卡在线观看视频| 高潮毛片免费观看| 日韩性网站| 亚洲一区二区黄色| 亚洲天堂久久| av大片在线无码免费| yjizz视频最新网站在线| 国产对白刺激真实精品91| 国产精品欧美日本韩免费一区二区三区不卡| 国产精品30p| 国产1区2区在线观看| 欧美劲爆第一页| 在线视频97| 人妻91无码色偷偷色噜噜噜| 中文字幕在线看| 亚洲成aⅴ人在线观看| 亚洲天堂.com| 亚洲精品午夜无码电影网| 国禁国产you女视频网站| 91精品国产一区自在线拍| 久久综合五月婷婷| 四虎永久在线| 无码国内精品人妻少妇蜜桃视频| 亚洲国产天堂久久九九九| 中文字幕不卡免费高清视频| 久久精品人人做人人爽电影蜜月| 萌白酱国产一区二区| 在线毛片网站| 中文字幕乱码中文乱码51精品| 亚欧成人无码AV在线播放| 97精品伊人久久大香线蕉| 思思热在线视频精品| 亚洲欧美激情小说另类| 免费看av在线网站网址| 日本久久网站| 91丝袜在线观看|