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

        ?

        利用密度泛函理論研究α-聯(lián)噻吩體系H(C4H2S)nH的結(jié)構(gòu)和電子光譜

        2013-11-21 01:39:30趙建想
        化學(xué)研究 2013年5期
        關(guān)鍵詞:體系優(yōu)化實(shí)驗(yàn)

        寧 攀,趙建想

        (河南大學(xué) 化學(xué)化工學(xué)院 環(huán)境和分析科學(xué)研究所,河南 開封 475004)

        圖1 α-聯(lián)噻吩Fig.1 α-Oligothiophenes

        在過去的幾十年里,以α-位的C-C單鍵相連接的聯(lián)噻吩(α-聯(lián)噻吩)(如圖1所示)由于其具有平面的或是類似線性的共軛π體系結(jié)構(gòu)而被廣泛研究[1-5]. 它們突出的光學(xué)、電學(xué)、氧化還原、光致發(fā)光和電荷傳輸?shù)忍匦宰屗鼈冊谟袡C(jī)場效應(yīng)晶體管[6]、發(fā)光二極管[7]以及光伏打電池[8-9]等方面表現(xiàn)出廣泛的應(yīng)用和發(fā)展前景. 早在1974年,KUHN等人[10]就報(bào)道了α-五聯(lián)噻吩在Langmuir-Blodgett薄膜中的電流測量值;1999年,DENIS F等研究了α-聯(lián)噻吩的晶體和薄膜結(jié)構(gòu)[11];2011年,SIEGERT等人利用光致光電子能譜法研究了在氣相中將自由基陰離子接入中性α-聯(lián)噻吩的最低電子激發(fā)態(tài)[12]. 但這些研究的對象都是八個(gè)噻吩環(huán)以內(nèi)的α-聯(lián)噻吩,對于含有八個(gè)以上噻吩環(huán)的α-聯(lián)噻吩還未見報(bào)道.

        為更好地了解噻吩環(huán)個(gè)數(shù)與其結(jié)構(gòu)的關(guān)系,并預(yù)測較長的聯(lián)噻吩體系的光譜性質(zhì),本文作者采用密度泛函理論優(yōu)化了H(C4H2S)nH (n= 2~13)的基態(tài)和最低激發(fā)態(tài)的幾何構(gòu)型,在平衡構(gòu)型的基礎(chǔ)上進(jìn)一步計(jì)算了它們的垂直激發(fā)能和垂直發(fā)射能.

        1 計(jì)算方法

        在B3LYP/6-31G(d)和CAM-B3LYP/ 6-31G(d)水平下優(yōu)化了α-聯(lián)噻吩 H(C4H2S)nH(n= 2~13)的基態(tài)幾何構(gòu)型,并在相同水平下計(jì)算了它們的振動(dòng)頻率. 采用TD-B3LYP和TD-CAM-B3LYP方法優(yōu)化了它們的激發(fā)態(tài). 用含時(shí)密度泛函理論的TD-B3LYP方法,分別結(jié)合cc-pVDZ, cc-pVTZ以及6-31G(d)基組計(jì)算了α-聯(lián)噻吩H(C4H2S)nH (n= 2~13)的垂直躍遷能. 考慮到溶劑對有機(jī)分子光譜的影響,在計(jì)算的過程中使用了極化連續(xù)介質(zhì)(PCM)模型[13-14]. 最后,在TD-CAM-B3LYP/cc-pVDZ和TD-CAM-B3LYP/cc-pVTZ水平下計(jì)算了它們的垂直發(fā)射能. 以上所有計(jì)算都是利用Gaussian 09程序完成的.

        2 結(jié)果和討論

        2.1 基態(tài)幾何構(gòu)型和穩(wěn)定性

        本文中要討論的α-聯(lián)噻吩(以下簡稱α-nT,n代表噻吩環(huán)的個(gè)數(shù))有三種構(gòu)型,包括對稱性為C1和C2的螺旋上升型和環(huán)型的兩種順式結(jié)構(gòu)(以下分別簡稱cis-C1和cis-C2)以及對稱性為C1的稍有彎曲的帶狀結(jié)構(gòu)的反式構(gòu)型(以下簡稱trans),如圖2,圖3和圖4所示. 利用B3LYP和CAM-B3LYP方法優(yōu)化了α-nT的基態(tài)構(gòu)型,這兩種方法計(jì)算的結(jié)果非常接近,優(yōu)化得到的C-C鍵的鍵長和二面角相差很小.

        圖2 α-聯(lián)噻吩H(C4H2S)nH (n = 2~13) 的幾何結(jié)構(gòu)Fig.2 Structures of α-oligothiophenes H(C4H2S)nH (n = 2-13)

        圖3 B3LYP/6-31G(d) 方法下優(yōu)化得到的cis-13T的基態(tài)構(gòu)型Fig.3 Optimized structures of cis-13T at B3LYP/6-31G(d) level

        圖4 B3LYP/6-31G(d) 方法下優(yōu)化得到的tran-13T的基態(tài)構(gòu)型Fig.4 Optimized structures of tran-13T (S0) at B3LYP/6-31G(d) level

        B3LYP和CAM-B3LYP兩種方法計(jì)算得到的振動(dòng)頻率都是正值,說明優(yōu)化得到的結(jié)構(gòu)對應(yīng)勢能面上的能量極小點(diǎn),均是穩(wěn)定構(gòu)型. 而且它們的最小彎曲振動(dòng)頻率都小于100 cm-1,可見此類化合物是“非剛性的”,而且隨著n的增加,彎曲振動(dòng)頻率逐漸減小,化合物會(huì)越來越“軟”.

        2.2 基態(tài)和激發(fā)態(tài)的比較

        TD-B3LYP方法優(yōu)化的α-nT (n= 10)的第一激發(fā)態(tài)21A (21B)的構(gòu)型繪制于圖5.cis-C1和cis-C2構(gòu)型由基態(tài)的螺旋上升型和環(huán)型結(jié)構(gòu)變成了激發(fā)態(tài)的環(huán)狀結(jié)構(gòu). 第一激發(fā)態(tài)的相鄰噻吩環(huán)之間的二面角接近于0°或180°,即近似于平面;而C-C鍵的鍵長跟基態(tài)相比均變小,這種縮短是由于激發(fā)態(tài)分子的平面性增強(qiáng),電子離域程度增大導(dǎo)致的. TD-CAM-B3LYP方法優(yōu)化的結(jié)果也是如此.

        圖5 TD-B3LYP方法下優(yōu)化得到的α-10T的第一激發(fā)態(tài)的構(gòu)型Fig.5 Optimized structures of α-10T (S1) at TD-B3LYP/6-31G(d) level

        2.3 能隙

        圖6是體系在氣相中運(yùn)用B3LYP方法得到的HOMO,LUMO能量及其能隙ε與n的關(guān)系圖. 從圖中可以看出,隨著n的增加,體系的HOMO能量逐漸升高,而LUMO能量則逐漸降低,因此LUMO與HOMO的能量差逐漸減小,我們推測是由于體系的共軛程度隨著n的增大而增強(qiáng),導(dǎo)致前線軌道的能隙ε減小,從HOMO到LUMO的電子躍遷逐漸變得容易,從而導(dǎo)致吸收波長增長.

        圖6 B3LYP方法下得到的氣相中的α-nT的HOMO能量,LUMO能量和能隙εFig.6 The energy levels of HOMO, LUMO and ε of α-nT by B3LYP method in gas phase

        2.4 垂直躍遷能

        表1-表3列出了TD-B3LYP方法分別結(jié)合cc-pVDZ, cc-pVTZ 以及 6-31G(d)基組計(jì)算得到的α-nT (n= 2~13)的21A (21B)←X1A電子躍遷的垂直躍遷能和相應(yīng)的實(shí)驗(yàn)值.

        表1 cis-α-nT (C1)的S0→S1的垂直躍遷能(ΔE,eV)和振子強(qiáng)度(f)Table 1 The vertical transition energies (ΔE, eV) and oscillator strengths (f) of cis-α-nT (C1) of S0→S1

        aBy TD-B3LYP/cc-pVDZ in gas;bBy TD-B3LYP/cc-pVDZ in dichloromethane

        表2 cis-α-nT (C2)S0→S1的垂直躍遷能(ΔE, eV)和振子強(qiáng)度(f)Table 2 The vertical transition energies (ΔE, eV) and oscillator strengths (f) of cis-α-nT (C2) of S0→S1

        aBy TD-B3LYP/cc-pVDZ in gas;bBy TD-B3LYP/cc-pVTZ in gas;

        cBy TD-B3LYP/6-31G(d) in dichloromethane;

        dBy TD-B3LYP/cc-pVTZ in dichloromethane

        表3 trans-α-nT (n = 2~13)S0→S1的垂直躍遷能(ΔE, eV)和振子強(qiáng)度(f)和實(shí)驗(yàn)值Table 3 The vertical transition energies (ΔE, eV) and oscillator strengths (f) of trans-α-nT(n =2-13) of S0→S1 and the experimental value

        aBy TD-B3LYP/cc-pVDZ in gas;bBy TD-B3LYP/cc-pVTZ in gas;

        cBy TD-B3LYP/cc-pVDZ in dichloromethane;d實(shí)驗(yàn)值摘自參考文獻(xiàn)[12]

        2.4.1 cis-α-nT (C1)

        表1列出了氣相和二氯甲烷溶劑中在TD-B3LYP/cc-pVDZ水平下計(jì)算得到的cis-α-nT (C1) (n= 2~13)的21A (21B)←X1A電子躍遷的垂直躍遷能. 表中數(shù)據(jù)顯示,隨著n的增加,除了7T以外,垂直躍遷能逐漸減小,這表明21A (21B)←X1A電子躍遷隨著n的增大變得越來越容易,圖6很直觀的描述了這一現(xiàn)象.

        因?yàn)榧ぐl(fā)態(tài)21A (21B)都是由HOMO→LUMO的電子躍遷引起的,所以cis-α-nT (C1)的能隙與21A (21B)←X1A躍遷的躍遷能遞減的趨勢密切相關(guān). 由圖7可見,隨著噻吩環(huán)個(gè)數(shù)n的增加,21A (21B)←X1A躍遷的垂直躍遷能逐漸減小,并趨于平緩,這與HOMO-LUMO能隙變化趨勢是一致的.

        2.4.2 cis-α-nT (C2)

        在不同水平下計(jì)算的cis-α-nT (C2) (n= 2~13)的21A (21B)←X1A電子躍遷的垂直躍遷能結(jié)果列于表2中. 數(shù)據(jù)表明,不同的基組和溶劑對此體系的躍遷能影響不大.

        在TD-B3LYP/cc-pVDZ方法下計(jì)算的相鄰躍遷能之間的差值為:0.73, 0.43, 0.27, 0.19, 0.14, 0.09, 0.08, 0.06, 0.05, 0.06 和 0.04 eV. 從這些數(shù)據(jù)可以看出,隨著n的增大,躍遷能呈現(xiàn)遞減的趨勢. 在此方法下計(jì)算的躍遷能與體系大小n的依賴關(guān)系繪于圖8.

        圖7 cis-α-nT (C1)的21A (21B)←X1A的垂直躍遷能與n的關(guān)系曲線Fig.7 Plot of the vertical transition energies of the 21A (21B)←X1A transition for cis-α-nT (C1) vs n

        圖8 cis-α-nT (C2)的21A (21B) ←X1A的垂直躍遷能與n的關(guān)系曲線Fig.8 Plot of the vertical transition energies of the 21A (21B)←X1A transition for cis-α-nT (C2) vs n

        從圖8中可知體系躍遷能與體系大小n呈現(xiàn)非線性關(guān)系. 擬合表2中在TD-B3LYP/cc-pVDZ水平下得到的躍遷能數(shù)據(jù),得到cis-α-nT (C2) (n= 2~13)體系氣相中的21A (21B)←X1A躍遷能和體系大小n的關(guān)系式:

        ΔE=A+B1e-n/D1+B2e-n/D2

        (1)

        式中,A= 1.685 39,B1= 1.451 12,B2= 4.486 21,D1= 9.411 89, 和D2= 1.644 97. 擬合誤差和相關(guān)系數(shù)分別為0.000 03 eV和 0.999 95. 用式(1)可以預(yù)測,當(dāng)體系足夠大時(shí),體系的躍遷能ΔE將收斂至1.68 eV (735.57 nm). 為了更好的說明這種非線性關(guān)系,其他水平下計(jì)算的躍遷能與體系大小n的依賴關(guān)系也繪于圖7中.

        2.4.3 trans-α-nT

        表3列出了在不同理論水平下計(jì)算得到的trans-α-nT (n= 2~13)的21A←X1A電子躍遷的垂直躍遷能(ΔE),振子強(qiáng)度(f)和相應(yīng)的實(shí)驗(yàn)值. 由表3數(shù)據(jù)可知,理論計(jì)算的2T-4T的垂直躍遷能與實(shí)驗(yàn)值符合得很好,但是隨著鏈的增長誤差越來越大. 表4列出了以前的理論研究和實(shí)驗(yàn)上得到的trans-α-nT的第一激發(fā)態(tài)的躍遷能數(shù)據(jù)以及本文中的理論值. 從表中數(shù)據(jù)可以發(fā)現(xiàn),在TD-B3LYP/cc-pVDZ和TD-B3LYP/cc-pVTZ水平下計(jì)算得到的數(shù)據(jù)與先前報(bào)道的DFT/MRCI和 CASPT2 方法下計(jì)算結(jié)果都很接近. 考慮到計(jì)算成本,我們在TD-B3LYP/cc-pVTZ水平下計(jì)算了n= 2~4的聯(lián)噻吩體系.

        表4 氣相中的第一單重激發(fā)態(tài)能量(eV)Table 4 The S1 singlet state vertical transition energies in gas phase (all energies in eV)

        aBy TD-B3LYP/cc-pVDZ;bBy TD-B3LYP/cc-pVTZ;

        cCASPT2// MP2, 摘自參考文獻(xiàn) [15-16];dtrans-α-2T氣相中的LIF光譜, 摘自參考文獻(xiàn) [17];

        eDFT/MRCI, 摘自參考文獻(xiàn) [12];f摘自參考文獻(xiàn) [12]

        實(shí)驗(yàn)值摘自參考文獻(xiàn)[12]圖9 trans-α-nT的21A←X1A的垂直躍遷能與n的關(guān)系曲線Fig.9 Plot of the vertical transition energies of the 21A←X1A transition for trans-α-nT vs n

        圖9trans-α-nT是在TD-B3LYP/cc-pVDZ水平下氣相和二氯甲烷溶劑中的21A←X1A垂直躍遷能與體系大小n的關(guān)系. 從預(yù)測的垂直躍遷能可以看出,隨著n的增大能量在逐漸減小,說明隨著n的增加21A←X1A電子躍遷越來越容易.

        將氣相和二氯甲烷溶劑中TD-B3LYP/cc-pVDZ水平下得到的垂直躍遷能數(shù)據(jù)用公式(1)進(jìn)行擬合,可以得到: 在氣相中,A=1.830 79,B1=1.930 39,B2=4.309 28,D1=4.800 47, 和D2=1.365 59;擬合誤差和相關(guān)系數(shù)分別為0.000 05 eV和 0.999 93. 當(dāng)體系足夠大時(shí),ΔE將收斂至1.83 eV (677.45 nm);在二氯甲烷溶劑中,A=1.756 45,B1=1.917 29,B2=4.694 34,D1=4.337 55 和D2=1.297 1;擬合誤差和相關(guān)系數(shù)分別為0.000 05 eV和0.999 94. 當(dāng)體系足夠大時(shí),ΔE將收斂至1.75 eV (708.42 nm). 這說明溶劑對體系的垂直躍遷能有一定的影響.

        從圖8可以看出,不同水平下計(jì)算的垂直躍遷能在n比較小的時(shí)候與實(shí)驗(yàn)值符合得很好,但隨著鏈的增長,與相應(yīng)的實(shí)驗(yàn)結(jié)果之間的誤差在增加.

        2.5 垂直發(fā)射能

        在TD-CAM-B3LYP方法優(yōu)化的第一激發(fā)態(tài)構(gòu)型的基礎(chǔ)上,計(jì)算了α-聯(lián)噻吩 H(C4H2S)nH (n= 2~10)的21A (21B)→X1A的電子發(fā)射能. 利用TD-CAM-B3LYP/cc-pVDZ 和TD-CAM-B3LYP/cc-pVTZ方法計(jì)算的發(fā)射波長λ(nm)和振子強(qiáng)度f以及相應(yīng)的實(shí)驗(yàn)值列在表5中. 不同水平下計(jì)算的發(fā)射波長相互之間很接近,而且都跟實(shí)驗(yàn)值符合得較好,說明TD-CAM-B3LYP方法是計(jì)算發(fā)射光譜的一種行之有效的方法. 計(jì)算的發(fā)射光譜與相應(yīng)的吸收光譜相比呈現(xiàn)明顯的紅移,是由于激發(fā)態(tài)弛豫造成的.

        表5 TD-CAM-B3LYP計(jì)算的α-nT 的垂直發(fā)射能(nm)和振子強(qiáng)度以及相應(yīng)的實(shí)驗(yàn)值Table 5 The vertical emission energies (λ in nm) and oscillator strengths (f) of α-nT at TD-CAM-B3LYP level, along with the experimental value

        aBy TD-CAM-B3LYP/cc-pVDZ;bBy TD-CAM-B3LYP/cc-pVTZ;c摘自參考文獻(xiàn)[18]

        3 結(jié)論

        分別運(yùn)用B3LYP/6-31G(d) 和 CAM-B3LYP/6-31G(d)兩種密度泛函理論優(yōu)化了α-聯(lián)噻吩H(C4H2S)nH (n= 2~13)的基態(tài)構(gòu)型. 振動(dòng)頻率分析結(jié)果表明得到的三種結(jié)構(gòu)都是穩(wěn)定構(gòu)型. 用TD-DFT方法優(yōu)化了對應(yīng)的第一激發(fā)態(tài)的平衡構(gòu)型,與基態(tài)相比,激發(fā)態(tài)構(gòu)型中C-C鍵的長度縮短了,這是由于激發(fā)態(tài)分子的平面性增強(qiáng),電子的離域程度增大導(dǎo)致的.

        在TD-B3LYP/cc-pVDZ、TD-B3LYP/cc-pVTZ和TD-B3LYP/6-31G(d)三種理論水平下,預(yù)測的α-聯(lián)噻吩體系三種構(gòu)型的21A (21B) ←X1A躍遷的垂直躍遷能都具有相似的非線性尺寸依賴. 根據(jù)計(jì)算結(jié)果,導(dǎo)出了這個(gè)體系的垂直躍遷能與體系大小n的解析關(guān)系式.

        體系的垂直發(fā)射能是在TD-CAM-B3LYP/cc-pVDZ 和 TD-CAM-B3LYP/cc-pVTZ水平下計(jì)算得到的,預(yù)測的最大發(fā)射波長與相應(yīng)的實(shí)驗(yàn)值符合得較好. 上述結(jié)果將為進(jìn)一步研究更長的聯(lián)噻吩體系的發(fā)光性質(zhì)提供必要的理論信息.

        參考文獻(xiàn):

        [1] CASADO J, PAPPENFUS T M, MILLER L L, et al. Nitro-functionalized oligothiophenes as a novel type of electroactive molecular material: spectroscopic, electrochemical, and computational study[J]. J Am Chem Soc, 2003, 125(9): 2524-2534.

        [2] ZHANG Xin Nian, MATZGER A J. Effect of ring fusion on the electronic absorption and emission properties of oligothiophenes[J]. J Org Chem, 2003, 68(25): 9813-9815.

        [3] CARAS-QUINTERO D, BAUERLE P. Synthesis of the first enantiomerically pure and chiral, disubstituted 3,4-ethylenedioxythiophenes (EDOTs) and corresponding stereo-and regioregular PEDOTs[J]. Chem Commun, 2004 (8): 926-927.

        [4] MURPHY A R, FRECHET J M J, CHANG P, et al. Organic thin film transistors from a soluble oligothiophene derivative containing thermally removable solubilizing groups[J]. J Am Chem Soc, 2004, 126(6): 1596-1597.

        [5] ZHANG Xin Nian, CTéA P, MATZGER A J. Synthesis and structure of fusedα-oligothiophenes with up to seven rings[J]. J Am Chem Soc, 2005, 127(30): 10502-10503.

        [6] BRUSSO J L, HIRST O D, DADVAND A, et al. Two-dimensional structural motif in thienoacene semiconductors: synthesis, structure, and properties of tetrathienoanthracene isomers[J]. Chem Mater, 2008, 20(7): 2484-2494.

        [7] PEREPICHKA I F, PEREPICHKA D F, MENG H. Light-emitting polythiophenes[J]. Adv Mater, 2005, 17(19): 2281-2305.

        [8] KWON T H, ARMEL V, NATTESTAD A, et al. Dithienothiophene (DTT)-based dyes for dye-sensitized solar cells: synthesis of 2,6-dibromo-DTT[J]. J Org Chem, 2011, 76(10): 4088-4093.

        [9] BOUDREAULT P L T, NAJARI A, LECLERC M. Processable low-bandgap polymers for photovoltaic applications[J]. Chem Mater, 2011, 23: 456-469.

        [10] SCHOELER U, TEWS K H, KUHN H. Potential model of dye molecule from measurements of the photocurrent in monolayer assemblies[J]. J Chem Phys, 1974, 61: 5009.

        [11] FICHOU D. Structural order in conjugated oligothiophenes and its implications on opto-electronic devices[J]. J Mater Chem, 2000, 10: 571-588.

        [12] SIEGERT S, VOGELER F, MARIAN C M, et al. Throwing light on dark states ofα-oligothiophenes of chain lengths 2 to 6: radical anion photoelectron spectroscopy and excited-state theory[J]. Phys Chem Chem Phys, 2011, 13: 10350-10363.

        [14] TOMASI J, PERSICO M. Molecular interactions in solution: an overview of methods based on continuous distributions of the solvent[J]. Chem Rev, 1994, 94(7): 2027-2094.

        [15] RUBIO M, MERCHN M, POU-AMéRIGO R, et al. The low-lying excited states of 2,2′-bithiophene: a theoretical analysis[J]. Chem Phys Chem, 2003, 4(12): 1308-1315.

        [16] RUBIO M, MERCHN M, ORTI E. A theoretical study on the low-lying excited states of 2,2′:5′,2″-terthiophene and 2,2′:5′,2″:5″,2?-quaterthiophene[J]. Chem Phys Chem, 2005, 6(7): 1357-1368.

        [17] BIRNBAUM D, FICHOU D, KOHLER B E. The lowest energy singlet state of tetrathiophene, an oligomer of polythiophene[J]. J Chem Phys, 1992, 96: 165-169.

        [18] BECKER R S, MACANITA A L, ELISEI F, et al. Comprehensive evaluation of the absorption, photophysical, energy transfer, structural, and theoretical properties ofα-oligothiophenes with one to seven rings[J]. J Chem Phys, 1996, 100(48): 18683-1869.

        猜你喜歡
        體系優(yōu)化實(shí)驗(yàn)
        記一次有趣的實(shí)驗(yàn)
        超限高層建筑結(jié)構(gòu)設(shè)計(jì)與優(yōu)化思考
        民用建筑防煙排煙設(shè)計(jì)優(yōu)化探討
        關(guān)于優(yōu)化消防安全告知承諾的一些思考
        一道優(yōu)化題的幾何解法
        構(gòu)建體系,舉一反三
        做個(gè)怪怪長實(shí)驗(yàn)
        NO與NO2相互轉(zhuǎn)化實(shí)驗(yàn)的改進(jìn)
        實(shí)踐十號上的19項(xiàng)實(shí)驗(yàn)
        太空探索(2016年5期)2016-07-12 15:17:55
        “曲線運(yùn)動(dòng)”知識體系和方法指導(dǎo)
        国产精品亚洲一区二区三区妖精| 亚洲高潮喷水无码av电影| 精品午夜福利1000在线观看| 无码高清视频在线播放十区| 亚洲精品成人一区二区三区| 精品亚洲麻豆1区2区3区| 免费国产黄网站在线观看| 国产精品一区二区资源| 一区二区三区在线观看视频| 亚洲午夜狼人综合影院| 首页 综合国产 亚洲 丝袜| 亚洲欧美中文在线观看4| 亚洲女同恋中文一区二区| 久久久99精品免费视频| 男女裸交无遮挡啪啪激情试看| 国产91网| 久久国产女同一区二区| 亚洲大尺度无码无码专区| 成 人 免费 黄 色 视频 | 亚洲精品有码在线观看| 蜜桃在线观看视频在线观看| 手机在线亚洲精品网站| 中文字幕精品一区二区2021年| 中文字幕一区二区三区四区在线| 日本特殊按摩在线观看| 人妻 丝袜美腿 中文字幕| 日日碰狠狠丁香久燥| 久久久精品中文无码字幕| 三级日本理论在线观看| 精品国产av色一区二区深夜久久| 欧美日本亚洲国产一区二区| 国产一区二区三区经典| 久久久国产精品无码免费专区| 亚洲精品aa片在线观看国产| 无码超乳爆乳中文字幕| 亚洲1区第2区第3区在线播放 | 久久精品亚洲成在人线av乱码| 成人国产精品一区二区视频| 国产一级毛片AV不卡尤物| 白丝美女扒开内露出内裤视频| 色偷偷色噜噜狠狠网站30根|