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

In C+n(n=1—10)團簇的密度泛函理論研究?

2018-06-19 10:03:48張陳俊王養麗陳朝康
物理學報 2018年11期
關鍵詞:結構

張陳俊 王養麗 陳朝康

1 引 言

隨著摻雜碳團簇在星際空間被發現[1?6],這一領域的研究迅速引起了理論和實驗物理學家們的極大興趣.因為把別的原子摻到碳原子團簇中后,團簇一系列的結構性質(如結構特征、電子特性、能級結構、穩定性等)都會發生改變,對這些類型團簇的研究有助于引導人們對新物質的發現,也有助于指導新型固體材料的制備與研發.研究者們用激光濺射固體樣品實驗可以觀察到一系列含雜原子的碳團簇[7?10].人們也期望從理論上進一步揭示各種富碳摻雜團簇的形成機理,故開展了大量對摻雜碳團簇的研究.關于摻第IIIA族元素碳基團簇的理論和實驗研究也引起了人們極大的關注.Chuchev和Belbruno[11]以及Liang和Zhang[12]研究了CnB,CnB2(n=4—10)和C2+nB10?n(n=0—10)團簇的結構和穩定性;Wang等[13],Nakajima[14]以及Pascoli和Lavendy[15]對CnB?(n<13),CnB?(n 6 7)和CnB+(n=9—15)團簇的結構進行了報道;Largo等[16]以及Li和Tang[17]用密度泛函理論對CnAl(n=0—10)和CnAl/CnAl+/CnAl?/CnAl2(n=0—10)團簇進行了理論研究;Chertihin等[18]用Laser-Ablated實驗方法得到了CAl2和C2Al2團簇的紅外譜圖.然而將重元素摻雜到碳團簇中的理論與實驗研究很困難,因為重原子電子數多,而且有較強的相對論效應.作為第IIIA族中的重金屬元素的Ga,In,Tl摻雜到碳團簇的理論研究僅見于我們之前關于GaCn團簇的結構和穩定性的理論研究的報道[19].因此,對摻雜第IIIA族元素的碳團簇的認識還不夠,再加上此前得到的構型數量也相對較少,難以總結出這一系列二元碳團簇的結構演化規律.鑒于這種情形,本文以原子團簇的幾何結構和電子性質及其隨團簇尺寸演變為主線,對In摻雜碳的陽離子團簇InC+n的結構、穩定性、磁性和極化率進行系統的理論研究,以期能分析總結出含第IIIA族元素的碳團簇的結構規律及所具有的物理化學性質.研究結果可以成為相關實驗的理論依據,對實驗研究產生導向性,也將為研究更大尺寸、更多元化的團簇提供理論基礎.

2 計算方法

采用密度泛函(DFT)理論中的雜化密度泛函B3LYP方法(源于Becke,Lee,Yang和Parr的姓氏),是將Hartree-Fock的交換能和DFT中的交換能有效組合而得到的,對本系列團簇的計算有實用性和有效性.采用由Hay和Wadt開發的LANL2DZ基組[20],優勢是在基組中選用了雙ζ價電子,由于較重元素原子的核較大,所以需要將這些原子核附近的電子用一種包含了相對論效應的有效核電勢(贗勢場ECP)近似來處理,在含有重元素的計算中是較好的基組.在前期計算中性InCn團簇時,同時使用了兩種基組:一種是C原子用6-311+G*,In原子用LANL2DZ基組;而另一種是C原子和In原子都用LANL2DZ基組.比較發現,兩種方案得到的有關各團簇基態結構、成鍵情況、穩定性的奇偶振蕩規律以及極化率等規律都相同.另一方面,一些有關重元素摻雜團簇的計算文獻中[21?24]使用的基組也是LANL2DZ.基于以上原因,綜合考慮計算的經濟性和有效性,在計算陽離子團簇InC+n時,直接使用了LANL2DZ基組.在計算過程中,首先用窮舉法對團簇的各種可能構型進行了結構設計和點群確定;優化過程對所有可能存在的各種構型都考慮了3種自旋多重態,在計算結果中篩選出能量較低的一些構型;其次是將初次優化所得能量較低的穩定結構進行更為精確的優化和計算,最終確定了團簇的基態結構.之前大量的研究已表明,這種采用分步優化計算的方法既可以保證計算結果的精確性,又可以節省大量的計算時間.所有的計算工作均在Gaussian 03程序包下進行.

3 結果與分析

經過大量的計算,最終得到了各團簇的各種可能存在的構型,在此只給出基態結構和另兩種能量相對較低的結構.對于每一InC+n團簇都考慮了自旋多重度為1,3,5的情形.在本文中各個態都是真實的平衡態,它們的振動頻率都沒有虛頻;所列舉的總能量中都包含零點振動能.

3.1 幾何結構

在圖1和圖2中展示出了直線型的基態構型和另兩種能量相對較低的結構.表1中給出了基態構型的電子態、總能量、最低振動頻率和相對能.表2中列出了團簇基態結構的電子態、最低振動頻率、總能量、自旋污染期望值、偶極距、相對能.

InC+的基態是In原子處于碳鏈一端的直線型結構,它的基態(3Σ)能量比單態和五重態分別低43.1 kcal/mol(1 cal=4.184 J)和40.4 kcal/mol.團簇的In—C鍵長為2.3187?.

InC+2的能量最低結構也是直線型,另一種能量較低的是三角構型,僅比直線型高17.3 kcal/mol,所以三角構型也被看作是穩定存在態.由表2可以看出,團簇的基態是三重態(3Σ),團簇的In—C鍵和C—C鍵長分別為2.8495?和1.2726?.

InC+3為In原子位于一端的線狀結構(In-CCC),這與BC+3[25]的基態為圓環有所不同.它的1Σ態能量最低.第二和第三穩定結構分別是平面四邊形結構和In原子位于C3單圓環一端的結構,它們的能量分別比基態的高58.4 kcal/mol和39.2 kcal/mol.優化得到的結果顯示團簇的In—C鍵長為2.7432?,C1—C2鍵長為1.2903?,C2—C3鍵長為1.3427?.可以看出In—C鍵明顯比C—C鍵要長,處于末端的C2—C3鍵比C1—C2鍵長一點.

圖2 InC+n團簇的非線型穩定結構Fig.2.Optimized geometries of In C+n non-linear isomers.

表1 InC+n團簇最穩定結構的電子態、最低振動頻率f lv、總能量E、相對能E relTable 1.Electronic state,the lowest vibrational frequencies f lv,total energies E,and relative energies E r el for InC+n clusters.

InC+4的能量最低結構是線狀結構(InCCCC),基態為3A.其扇形結構比基態能量高出42.2 kcal/mol.另一個能量較低的結構是C鏈一端連Ga原子,另一端連C4平面四邊形結構.In—C鍵長為2.6212 ?,C1—C2/C2—C3/C3—C4鍵長分別為1.3094/1.3091/1.3402?,可以看出處于末端的C3—C4鍵明顯要比其他兩個長,說明失去電子使得C鏈末端的C—C鍵相互作用減弱.

InC+5中線狀結構依然是最穩定結構.亞穩定結構是C鏈一端連In原子另一端連C3單圓環結構,它的能量比基態的高55.5 kcal/mol.第三穩定結構是圓環或類扇形結構,相比線性結構能量高出67.5 kcal/mol.經過優化得出團簇的In—C鍵長為2.5697 ?,C1—C2/C2—C3/C3—C4/C4—C5的鍵長分別為1.2760/1.3191/1.2796/1.3349?,可以看出C—C鍵存在有規律的振蕩,Codd—Ceven鍵比Ceven—Codd鍵要短.

表2 InC+n團簇線型結構的電子態、最低振動頻率f lv、總能量E,自旋污染期望值?s2?、偶極距μ、相對能E relTable 2.Electronic state,the lowest vibrational frequencies f lv total energies E?s2?expectation values,dipole momentsμ,and relative energies E rel for the linear isomers InC+n clusters.

InC+6中最穩定的結構是線狀結構(InCCCCCC),基態為3Σ態.第二穩定結構是圓環狀結構,比基態能量高出36.5 kcal/mol.團簇的第三穩定結構是C6單圓環上連一個Ga原子結構.In—C鍵長為2.5244 ?,C1—C2/C2—C3/C3—C4/C4—C5/C5—C6鍵長分別為1.2894/1.3124/1.2893/1.2980/1.3449?,可以看出,處于C鏈末端的C5—C6鍵明顯要比其他的C—C鍵長,說明失去電子使得C鏈末端的C—C鍵相互作用減弱,另外還可以看出除末端的C5—C6鍵外,其他的C—C鍵存在奇偶振蕩規律.

InC+7中In原子位于一端的線狀結構(InCCCCCCC)在它的1A態能量最低.另外兩個能量相對較低的結構分別是圓環狀結構和C7單圓環上連一個In原子結構,它們的能量分別比基態高50.4 kcal/mol和45.1 kcal/mol. 團簇的In—C鍵長為2.4824 ?,C1—C2/C2—C3/C3—C4/C4—C5/C5—C6/C6—C7鍵長分別為1.2712/1.3263/1.2658/1.3050/1.2841/1.3317?,C—C鍵存在奇偶振蕩特性.

InC+8中最穩定結構也是In原子位于一端的線狀結構,具有,3Σ態.第二穩定結構是圓環狀結構,比最穩定結構高39.0 kcal/mol.第三穩定結構是C8單圓環上連一個In原子結構,相對線狀結構高出50.7 kcal/mol.團簇的In—C鍵長2.4535?,C1—C2/C2—C3/C3—C4/C4—C5/C5—C6/C6—C7/C7—C8鍵 長 分 別 為1.2828/1.3174/1.2809/1.2938/1.2915/1.2940/1.3395?.可以看出,處于C鏈末端的C7—C8鍵明顯要比其他的C—C鍵長,說明失去電子使得C鏈末端的C—C鍵相互作用減弱.

InC+9的基態(1A)結構也是線狀結構(C∞v).另外兩個能量相對較低的結構分別是圓環狀結構和C9單圓環上連一個In原子結構,它們的能量分別比基態的高46.6 kcal/mol和47.9 kcal/mol.團簇的In—C 鍵長為2.4302?,C1—C2/C2—C3/C3—C4/C4—C5/C5—C6/C6—C7/C7—C8/C8—C9鍵長分別為1.2680/1.3317/1.2621/1.3111/1.2704/1.3012/1.2872/1.3293?,C—C鍵也存在奇偶振蕩特性.

In中最穩定結構也是In原子位于一端的線狀結構,與團簇AlC+10[26](基態結構則為圓環狀)不同.它的基態為3A.第二穩定結構是圓環狀結構,相對差16.3 kcal/mol.第三穩定結構是C10單圓環上連一個In原子結構,相比基態高16.5 kcal/mol.團簇的In—C鍵長為2.4236?,C1—C2/C2—C3/C3—C4/C4—C5/C5—C6/C6—C7/C7—C8/C8—C9/C9—C10鍵長分別為1.2743/1.3258/1.2723/1.3020/1.2840/1.2915/1.2917/1.2961/1.3327?,可以看出,處于C鏈末端的C9—C10鍵要比其他的C—C鍵長.

分析以上結果可以看出:團簇的全局能量最低構型是In原子位于一端的直線型或準直線型結構,而之前已經報道的BC+n和AlC+n團簇的基態結構中有部分為圓環狀.對于InC+n團簇基態而言,n為偶數的團簇是三重態,除InC+(三重態)外,n為奇數的團簇基態是單態.團簇基態這種單態和三重態交替現象與中性團簇InCn的電子結構有關[27,28],線性的InCn團簇有4n+3個價電子,核外電子排布為:

(除InC:(core)1σ22σ21π23σ1).由以上電子排布情況可以看出,n為奇數的中性團簇基態的電子結構為σ2π1,相應的失去一個電子的陽離子InC+n團簇基態電子結構為σ2(單態),而不是σ1π1,所以這種類型的團簇基態為單態;n為偶數的中性團簇基態電子結構為σ2π3,相對應的失去一個電子的陽離子InC+團簇基態電子結構為σ1π3(三重態).所以,n為偶數的團簇基態為三重態.

團簇的所有的電子基態自旋污染都不嚴重,很接近純自旋值.由表2還可以看出團簇基態的偶極距都不為零,而且團簇的偶極距是逐漸增大的(除個別n=2團簇),說明團簇的極性在逐漸增大.

團簇基態結構的In—C鍵都相對較長,而且比相應的中性團簇的In—C鍵都要長.而C—C鍵則都相對較短,形成的是較強的C=C雙鍵;n為奇數的團簇的C—C鍵存在明顯的奇偶振蕩規律(如InC+7團簇中的C—C鍵長分別為1.271/1.326/1.266/1.305/1.284/1.332?).C—C成鍵軌道上凈電荷分布存在的差異性,導致C—C鍵長出現明顯的長短振蕩規律,產生類聚乙炔化合物的結構特征,這一點與純團簇Cn類似.由自然鍵軌道法(NBO)分析基態結構各原子上的凈電荷分布[29,30]可以看出,在In和Cn相互作用形成InCn團簇基態結構的過程中,In原子的電荷發生了向C原子的轉移,這種電荷轉移的作用使得C—In鍵上的C原子呈負電性,In原子顯正電性,C—In鍵離子化.分析InCn到InC+n結構的電荷分布可知,電離1個電子的InC+n結構中原子得電子傾向比相應的中性結構明顯增大,C原子電負性略有降低.

通過觀察圖2,可以看出團簇的亞穩定或第三穩定結構中一定有圓環型或扇形結構.在得到的所有的前三穩定結構中,共有5種類型:直線型;圓環或類扇形型;C鏈一端連In原子另一端連C3的單圓環型;C鏈一端連In原子另一端連C4的平面四邊形結構;C10單圓環上連一個In原子的結構.

3.2 基態結構的穩定性

增量結合能(EIB)對團簇的相對穩定性分析起著十分重要的作用.為此計算了陽離子團簇InC+n的EIB.增量結合能的值越大表示穩定性越好,其定義為:InC+n→ InC+n?1+C(DN1),能量變化過程為:EIB=n E(C)+E(InC+n?1)? E(InC+n).團簇基態的增量結合能顯示在圖3,可以看出,陽離子團簇InC+n穩定性具有明顯的奇偶振蕩規律,n為1,3,5,7,9的團簇相對要穩定,是奇強偶弱. 為了進一步驗證團簇的相對穩定性,計算了能量的二階差分D2(En).其定義為:D2(En)=En+1+En?1?2En,其中n為團簇中所含C原子的個數,En表示InC+n團簇的總能量,能量的二階差分值越大表示相應的結構越穩定.從圖4中可以明顯看出,團簇基態的穩定性存在明顯的奇偶振蕩規律,穩定性隨尺寸的增加表現出的是奇強偶弱特性,這與由計算增量結合能得到的結果完全符合.而團簇基態穩定性的這種奇偶振蕩規律也是由電子結構所決定的,n為奇數的團簇基態(單態)電子填滿了π軌道(π4),而n為偶數的團簇基態(三重態)電子只填充了一半的π軌道(π2).故團簇穩定性產生了奇強偶弱的振蕩.此外,n為奇數團簇的C—C鍵產生類聚乙炔化合物的結構特征,這也是這類團簇穩定性更好的原因.

本文還計算了InCn團簇的絕熱電離能(Ead),定義為Ead=Ei?En(Ei為優化好的陽離子團簇基態能量,En為優化好的中性團簇基態能量).圖5繪出了團簇絕熱電離勢隨原子數目變化的趨勢.由圖5可以看出,Ead的值表現出明顯的奇偶振蕩效應,n為偶數的比n為奇數的團簇的值要大,說明n為偶數時中性團簇InCn較穩定,而陽離子團簇InC+n剛好相反,n為偶數時較不穩定.這進一步證實了InC+n團簇奇強偶弱的穩定特性.

圖3 陽離子團簇InC+n基態的增量結合能隨團簇尺寸的變化規律Fig.3.Incremental binding energies for InC+n ground state structures vs.the number of carbon atoms.

圖4 In C+n團簇基態的能量二階差分隨團簇尺寸的變化規律Fig.4. The second diff erence in energy for In C+n ground state structures vs. the number of carbon atoms.

圖5 團簇基態的絕熱電離能隨團簇尺寸的變化規律Fig.5.The adiabatic ionization potential for ground state structures vs.the number of carbon atoms.

3.3 磁性分析

摻雜團簇的磁性一直是人們重點研究對象之一,所以對InC+n團簇基態的磁矩做了相應的研究.首先計算出各混合團簇基態總磁矩的值是μB或2μB,還計算了系列團簇的平均每原子磁矩,計算結果如圖6所示.可以看出,團簇的磁矩隨團族尺寸的增加呈現出明顯的奇弱偶強振蕩趨勢,n=1,3,5,7,9的團簇出現了“淬滅”現象,這種現象主要原因在于這些團簇的自旋多重度為1,它們的能隙中最高占據軌道是非簡并的,而且它們的電子為滿殼層排布.

圖6 In C+n團簇平均每原子磁矩隨碳原子數的演化規律Fig.6.The magnetic moment of average atom for the InC+n clusters vs.the number of carbon atoms.

3.4 團簇的極化率

極化率也是團簇的眾多的物理化學性質之一,對極化率的研究可以幫助獲得其受外界影響的程度以及變化規律.這一物理量可以反映外電場對偶極距的影響程度,可以用來描述與物質的非線性作用,同時它也能對分子間色散力和碰撞散射截面等產生影響.基于上述重要用途,研究者也經常把團簇作為重點研究對象.本文計算了團簇基態線性結構的極化率張量?α?,極化率張量的平均值?α?/n,以及反映極化程度的極化率各向異性不變量?α,計算公式如下:

結果列于表3中.通過分析表3中的數據可以得出:極化率主要存在于xx,yy,zz方向,其他xy,xz,yz方向的分量基本都為零.通過詳細分析表中數據可知,張量的平均值?α?/n隨團簇尺寸的增大逐漸增大(除在InC+2處出現小幅振蕩).由此可得團簇的大小對極化率極為敏感,大尺寸的團簇很容易被外場極化,這主要是因為較強外場的作用破壞了電子的最初分布,再加上原子和原子之間的相互作用,整個團簇受極化的可能性增大,最終非線形光學效應凸顯.各向異性不變量?α在n 6 3的小范圍內做了微小波動,而后是單調上升(見圖7),這進一步證實團簇越大越容易被極化.圖7還顯示InC+的各向異性不變量最小,而且各方向的極化率相差不大,因此可知它最不易被外場極化.

圖7 In C+n團簇的極化率的各向異性不變量隨團簇尺寸的變化規律Fig.7.The anisotropic invariants of of the InC+n clusters polarizability vs.the size.

表3 In C+n團簇線型結構的極化率Table 3.Polarizability of the linear isomers InC+n clusters.

4 結 論

本文采用B3LYP/LANL2DZ方法對InC+n團簇進行了詳細研究,得到以下主要結論.1)In原子位于一端的直線型或準直線型結構是團簇的最穩定構型;n為偶數的團簇基態是三重態,n為奇數的團簇基態是單態(除InC+外);團簇基態的In—C鍵都相對較長,而C—C鍵存在明顯的奇偶振蕩規律;直線型團簇的所有電子基態自旋污染都很小,很接近純自旋值;直線型、圓環或類扇形型是相對較穩定的構型.2)通過對增量結合能和能量二階差分的分析可以得出,隨著團簇尺寸的增加,團簇表現出強烈的奇強偶弱振蕩規律,電離能Ead的計算結果進一步證實了這種振蕩規律的正確性.3)團簇的磁矩隨團族尺寸的增加呈現出明顯的奇弱偶強振蕩趨勢,n=1,3,5,7,9的團簇出現了“淬滅”現象.4)團簇基態結構的極化率張量主要分布在xx,yy,zz方向,xy,xz,yz方向的分量值基本上都為零,極化率張量的平均值及各向異性不變量都隨n的增大而增大,但在小尺寸范圍內出現微小振蕩.

感謝西北大學現代物理研究所提供的計算資源.

[1]Pauzat F,Ellinger Y 1989 Astron.Astrophys.216 305

[2]Maccarthy M T,Kalmus P,Gottlieb C A 1996 Astrophysics 467 125

[3]Cernicharo J,Guelin M 1996 Astron.Astrophysics 309 27

[4]Guelin M,Cernicharo J,Travers M J 1997 Astrophysics 37 1

[5]Schermann G,Grosser T,Hampel F 1997 Chem.Eur.J.3 1105

[6]Dembinski R,Bartik T,Bartik B 2000 J.Am.Chem.Soc.122 10

[7]Becker S,Dietze H 1988 Int.J.Mass Spectrom.82 287

[8]Consalvo D,Mele A,Stranges D 1989 Int.J.Mass Spectrom.91 319

[9]Liu Z Y,Wang C R,Huang R B 1995 Int.J.Mass Spectrom.141 201

[10]Liu Z Y,Huang R B,Tang Z C 1998 J.Chem.Phys.229 335

[11]Chuchev K,BelBruno J J 2004 J.Phys.Chem.108 5226

[12]Liang J X,Zhang C 2010 Acta Chim.Sin.68 7

[13]Wang C R,Huang R B,Liu Z Y 1995 Chem.Phys.Lett.242 55

[14]Nakajima A,Taguwa T,Nakao K 1995 J.Chem.Phys.103 2050

[15]Pascoli G,Lavendy H 2002 Opt.Plasma Phys.19 339

[16]Largo A,Redondo P,Barriento S 2002 J.Phys.Chem.A 106 4217

[17]Li G L,Tang Z C 2003 J.Phys.Chem.A 107 5317

[18]Chertihin G V,Andrews L,Taylor P R 1994 J.Am.Chem.Soc.116 3513

[19]Zhang C J,Jiang Z Y,Wang Y L 2013 Comput.Theor.Chem.1004 12

[20]Wadt W R,Hay P J 1985 J.Chem.Phys.82 284

[21]Jia L C,Zhao R N,Han J G,Sheng L S,Cai W P 2008 J.Phys.Chem.A 112 4375

[22]Li G L,Xing X P,Tang Z C 2003 J.Chem.Phys.118 6884

[23]Qi J Y,Dang L,Chen M D,Wu W,Zhang Q E,Au C T 2008 J.Phys.Chem.A 112 12456

[24]Li G L,Wang C Y 2007 J.Mol.Struct.824 48

[25]Wang L J,Zhang C J,Wu H S 2005 Acta Phys.Chim.Sin.21 244(in Chinese)[王利江,張聰杰,武海順2005物理化學學報21 244]

[26]Ma W J,Song X,Zhang X M,Wu H S 2010 Acta Phys.Chim.Sin.26 1396(in Chinese)[馬文瑾,宋翔,張獻明,武海順2010物理化學學報26 1396]

[27]Cheng L J 2012 J.Chem.Phys.136 104301

[28]Cheng L J,Yang J L 2013 J.Chem.Phys.138 141101

[29]Li L F,Cheng L J 2013 J.Chem.Phys.138 094312

[30]Feng Y Q,Cheng L J 2015 RSC Adv.5 62543

猜你喜歡
結構
DNA結構的發現
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
循環結構謹防“死循環”
論《日出》的結構
縱向結構
縱向結構
我國社會結構的重建
人間(2015年21期)2015-03-11 15:23:21
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
主站蜘蛛池模板: 日本欧美中文字幕精品亚洲| 精品综合久久久久久97超人该| 狠狠色噜噜狠狠狠狠色综合久 | 波多野结衣中文字幕久久| 一级毛片a女人刺激视频免费| 97在线免费| 国产精品浪潮Av| 亚洲国产91人成在线| 久久精品这里只有精99品| 欧美视频二区| 国产午夜一级毛片| 女人18毛片水真多国产| 亚洲无码91视频| 成人伊人色一区二区三区| 手机在线看片不卡中文字幕| 91毛片网| 无码区日韩专区免费系列| 欧美精品v| 国产美女自慰在线观看| 久视频免费精品6| 免费啪啪网址| 呦视频在线一区二区三区| 亚洲第一区在线| 四虎成人在线视频| 大陆精大陆国产国语精品1024| 青青草原国产av福利网站| 久久精品最新免费国产成人| 亚洲最大看欧美片网站地址| 日韩精品成人在线| 不卡午夜视频| 精品久久高清| 日韩大片免费观看视频播放| 国产一区成人| 91青青草视频在线观看的| 亚洲午夜天堂| 国产大片黄在线观看| 午夜国产理论| 国产九九精品视频| a毛片在线| 亚洲最猛黑人xxxx黑人猛交| 国产精品亚洲一区二区三区z| 亚洲无码视频一区二区三区 | julia中文字幕久久亚洲| 亚洲一区二区三区国产精华液| 污视频日本| 波多野结衣中文字幕一区二区| 人妻一本久道久久综合久久鬼色| 首页亚洲国产丝袜长腿综合| 久久黄色视频影| 操国产美女| 亚洲精品777| 国产亚洲视频中文字幕视频| 欧美精品v欧洲精品| 国产亚洲欧美日韩在线一区二区三区| 国产精品黄色片| 久久黄色小视频| 麻豆国产在线不卡一区二区| 婷婷久久综合九色综合88| 色天堂无毒不卡| 国产人妖视频一区在线观看| 欧美啪啪精品| 91麻豆精品国产91久久久久| 91青青在线视频| 国产福利一区视频| 亚洲欧美自拍一区| 美女视频黄又黄又免费高清| 刘亦菲一区二区在线观看| 国产91特黄特色A级毛片| 又猛又黄又爽无遮挡的视频网站| 日本免费福利视频| 成年人国产视频| 免费毛片视频| jizz国产在线| 一边摸一边做爽的视频17国产| 久久亚洲综合伊人| 国产精品yjizz视频网一二区| 国产一在线观看| 九九九国产| 欧美第二区| 亚洲丝袜中文字幕| 日韩天堂视频| 欧美日韩中文字幕在线|