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

        ?

        基于實時路徑積分的雙分子化學(xué)反應(yīng)動力學(xué)計算方法

        2022-12-06 14:38:36李永樂范文斌任偉
        關(guān)鍵詞:勢能課題組常數(shù)

        李永樂范文斌任偉

        (1.上海大學(xué)理學(xué)院量子與分子結(jié)構(gòu)國際中心,上海 200444;2.上海大學(xué)理學(xué)院上海市高溫超導(dǎo)重點實驗室,上海 200444)

        化學(xué)反應(yīng)速率常數(shù)的確定在燃燒[1-2]、星際分子[3]以及溫室效應(yīng)探索[4]的化學(xué)反應(yīng)動力學(xué)建模中起決定性作用.基于確定同一反應(yīng)、不同同位素的速率常數(shù),計算得到動力學(xué)同位素效應(yīng)(kinetic isotope effect,KIE),可以進一步研究化學(xué)反應(yīng)中的量子效應(yīng)(如量子隧穿(tunneling)和零點能(zero point energy,ZPE)等),以及動力學(xué)效應(yīng)(如再穿越(recrossing)等).另外,理論計算得到的動力學(xué)同位素效應(yīng)也反映了理論計算中使用的勢能面的性質(zhì),從而幫助進一步優(yōu)化勢能面[5].其他化學(xué)反應(yīng)動力學(xué)的相關(guān)物理量,如微正則系綜速率常數(shù)[6]、反應(yīng)幾率[7]、散射截面[8]等,可以反映化學(xué)反應(yīng)的動態(tài)信息,揭示反應(yīng)的微觀機理.由于化學(xué)動力學(xué)的實驗測量通常是一件困難且昂貴的工作[9],而精確的理論計算可以做出可靠預(yù)測,其結(jié)果還直接對應(yīng)一定的反應(yīng)機理,故有助于促進實驗技術(shù)的發(fā)展.因此,化學(xué)動力學(xué)的理論計算成為研究中不可或缺的重要手段[10-11].

        目前,氣相雙分子反應(yīng)動力學(xué)的理論計算為基于量子散射理論的含時波包動力學(xué)提供了較可靠的方法[12-13].結(jié)合高精度的神經(jīng)網(wǎng)絡(luò)(neural network,NN)勢能面[14-15],其計算結(jié)果可以與當今量子態(tài)分辨高精度測量結(jié)果高度符合.如在最近工作中,Chen等[16]基于新開發(fā)的NN勢能面上的波包動力學(xué)計算,并基于振動態(tài)分辨的反應(yīng)幾率和散射波函數(shù)分析,成功解釋了楊學(xué)明院士團隊在實驗中觀測到的Cl+CH4后向散射的微分截面在0.15 eV處的峰值,是來自于反應(yīng)物的重-輕-重動力學(xué)效應(yīng)而不是共振,是理論與實驗緊密結(jié)合的研究典范.另外,目前已發(fā)表的文獻表明,NN勢能面已經(jīng)能夠擬合多達十個原子的反應(yīng)體系[17],為計算復(fù)雜氣相雙分子反應(yīng)提供了堅實基礎(chǔ).但是,基于高精度勢能面的量子波包動力學(xué)方法,其計算量是根據(jù)原子數(shù)目呈指數(shù)級增長,目前仍不能滿足反應(yīng)動力學(xué)研究中對多原子反應(yīng)的計算需求.在快速計算方面,過渡態(tài)理論(transition state theory,TST)[18]以及基于TST的各種改進方法[19-23]是最普通的選擇.但TST計算忽略了量子力學(xué)的非定域性,以及反應(yīng)中的動力學(xué)效應(yīng),導(dǎo)致TST結(jié)果存在較大誤差,特別是低溫下TST的計算結(jié)果與實驗測量值往往相差若干數(shù)量級.另一種常用計算方法為準經(jīng)典軌線法(quasi-classical trajectory,QCT)[24],該方法根據(jù)量子力學(xué)設(shè)置反應(yīng)物的初態(tài),用經(jīng)典力學(xué)模擬散射過程,可以快速估算反應(yīng)幾率,并進一步求得散射截面、反應(yīng)速率常數(shù)等物理量[11,25-27].不過,由于QCT的軌跡不包含量子效應(yīng),可能會發(fā)生零點能泄漏問題,也沒有標準的指認產(chǎn)物態(tài)量子態(tài)分布的方法[28],故導(dǎo)致目前QCT只能作為研究反應(yīng)動力學(xué)的輔助工具.

        近年來,基于實時路徑積分的珠串分子動力學(xué)(ring-polymer molecular dynamics,RPMD,珠串直譯為環(huán)聚合物)[29-31]被廣泛應(yīng)用于各種氣相雙分子反應(yīng)的熱反應(yīng)速率及KIE的計算中[32].RPMD提供了一套比TST和QCT更精確、比量子散射理論耗費資源更少的計算方案.即使在較低溫度下的深度隧穿區(qū)域,RPMD也能給出足夠精確的結(jié)果[33].雖然標準RPMD方法的計算量比基于波函數(shù)的量子力學(xué)方法小很多,但是受限于目前勢能面計算能量的速度,目前RPMD可計算的系統(tǒng)仍然停留在七原子及以下反應(yīng),無法大幅超越量子波包動力學(xué)的應(yīng)用范圍,更無法達到TST和QCT廣泛適用的要求.如果采用計算速度快但精度低的解析勢能面,則即使應(yīng)用RPMD也無法計算出正確的反應(yīng)速率常數(shù)[34].目前,針對氣相多原子分子反應(yīng)動力學(xué)的理論研究仍然是化學(xué)動力學(xué)領(lǐng)域的前沿和難點之一[35],需要研究人員進一步完善計算方法.近年來,應(yīng)用RPMD計算氣相雙分子化學(xué)反應(yīng)的動力學(xué)又有了新的進展,該研究已經(jīng)擴展到多通道化學(xué)反應(yīng)、微正則系綜速率常數(shù)、表面化學(xué)反應(yīng)幾率和反應(yīng)截面的計算,本工作將對此進行綜述.下面首先簡述標準RPMD的化學(xué)速率常數(shù)計算方法,再簡單介紹反應(yīng)動力學(xué)理論的新進展與新應(yīng)用.

        1 方法

        1.1 路徑積分分子動力學(xué)

        考慮一個質(zhì)量為m的原子,引入P個串珠,有

        式中:h為普朗克常數(shù);p和q為串珠的經(jīng)典動量和坐標;為系統(tǒng)哈密頓算符;βP的定義為

        式中:kB為Boltzmam常數(shù);HP是所謂的珠串哈密頓量,

        且珠串中珠子間首尾相連,有q(P)=q(0),連接相鄰串珠的諧振子頻率為

        從化學(xué)反應(yīng)的量子理論[36]

        出發(fā),計算得到熱速率常數(shù).式中:Qr為反應(yīng)物的配分函數(shù);cfs(t)為流關(guān)聯(lián)函數(shù),

        實際計算中,運用Bennet-Chandler分解[37],將計算歸結(jié)為如下2個步驟.

        步驟1計算自由能變.引入2個分隔面s0和s1,其中s0放在2個反應(yīng)物相互作用可忽略的漸近區(qū)域內(nèi),s1放在過渡態(tài)區(qū)域內(nèi),由此定義出反應(yīng)坐標ξ=s0/(s0?s1).再基于路徑積分分子動力學(xué)(path-integral molecular dynamics,PIMD)沿著ξ的方向,運用增強采樣(enhanced sampling)方法,計算出平均力勢(potential of mean force,PMF)W(ξ),并確定其峰值對應(yīng)的反應(yīng)坐標ξ?,計算得到

        式中:R∞為設(shè)定的2個反應(yīng)物相互作用可忽略的距離;μR為反應(yīng)物的折合質(zhì)量.

        步驟2計算透射系數(shù).在雙dagger號ξ?處進行基于RPMD的的計算:

        式中:δ[ξ]為Dirac delta函數(shù);h[ξ]為Heaviside階躍函數(shù).此處根據(jù)RPMD雙分子反應(yīng)速率理論的基本假設(shè),認為RPMD中得到的對應(yīng)了量子理論中的cfs(t).

        最后,將這2個步驟的結(jié)果相乘得到最終速率.考慮到分子電子態(tài)的簡并有時還需要乘以一項跟電子配分函數(shù)有關(guān)的修正系數(shù)f(T),于是得到

        由于步驟2的計算可以模擬實時動力學(xué),RPMD速率包含再穿越(recrossing)效應(yīng),可以消去人為引入分隔面帶來的影響,因此分隔面的選取有一定的任意性,可在一定范圍內(nèi)任取[38].

        2013年,Suleimanov等[39]開發(fā)了程序RPMDrate,用于計算氣相中雙分子化學(xué)反應(yīng)速率常數(shù).由于使用了Fortran和Python編程語言,故該程序得到了廣泛應(yīng)用.

        在RPMD計算中,每個原子使用P個串珠(bead)來描述,這對計算量的需求呈線性增長.但由于標準RPMD方法中的時間步長通常為0.1 fs,而計算的軌跡需要數(shù)十納秒(ns),且低溫下往往需要多達數(shù)百個串珠來描述一個原子的運動,故導(dǎo)致計算效率又降低數(shù)百倍,難以滿足相關(guān)研究的需求.以下將通過2個例子來說明.

        OH+CH4在高溫下的拔氫過程是天然氣燃燒的重要反應(yīng)步驟,而該反應(yīng)在低溫下是甲烷在大氣層中消耗的關(guān)鍵反應(yīng)步驟[40].劉國平院士課題組于2005年對OH+CH4展開了系列研究工作,系統(tǒng)測量了該反應(yīng)的化學(xué)動力學(xué),如產(chǎn)物之間的振動量子態(tài)的關(guān)聯(lián)[41]及KIE[42]等.2018年,Li等[40]開發(fā)了基于置換不變多項式神經(jīng)網(wǎng)絡(luò)(permutationally invariant polynomial-NN,PIP-NN)的高精度全維勢能面.基于此PIP-NN勢能面,該課題組于2018年運用TST和QCT計算了OH+CH4的速率常數(shù)及OH/OD+CH4、OH+CH4/CD4、OH+12CH4/13CH4這3組KIE.結(jié)果表明,即使基于高精度勢能面,TST和QCT的估算結(jié)果也僅僅在高溫區(qū)可以與實驗結(jié)果符合,而在低于300 K的低溫區(qū)則存在較大誤差.其他課題組針對該體系的工作,大部分也是基于高精度電子結(jié)構(gòu)的計算結(jié)果運用TST或QCT,或者基于低精度的解析勢能面運用高精度的動力學(xué)方法進行計算,如RPMD[43].上述二者都因為引入了誤差更大的計算方法,所以得到的計算結(jié)果與實驗數(shù)據(jù)偏離更大.

        H+C3H8是燃料燃燒研究中重要的典型反應(yīng).早在1992年就有了關(guān)于該反應(yīng)的量子態(tài)分辨的動力學(xué)研究[44].2017年的反應(yīng)動力學(xué)實驗已經(jīng)可以精確測量不同反應(yīng)通道的速率常數(shù)[45],這給理論計算帶來了極大挑戰(zhàn),對于該反應(yīng)目前還沒有RPMD的計算結(jié)果.最近,Laude等[23]在對H及其同位素Mu和C3H8反應(yīng)的瞬子TST計算中,由于只考慮了勢能面上反應(yīng)路徑附近的區(qū)域,以及缺乏對再穿越等動力學(xué)效應(yīng)的計算,故得到的KIE比實驗值小1 000倍,誤差極大.

        由上述2個例子可見,對于復(fù)雜氣相雙分子反應(yīng),必須基于高精度勢能面且運用包含了必要的量子效應(yīng)的動力學(xué)方法(如RPMD),才能計算得到準確的反應(yīng)速率常數(shù)、KIE和反應(yīng)截面等反應(yīng)動力學(xué)中的關(guān)鍵物理量.要達到上述目標,首先必須解決標準RPMD方法計算速度慢的問題,其次必須解決標準RPMD方法本身存在的問題.例如,最近發(fā)現(xiàn)即使是采用足夠多串珠(bead)的標準RPMD計算,也存在數(shù)值不穩(wěn)定性和非遍歷性的問題[46].

        因此,開發(fā)RPMD計算的新方法必須確保其正確模擬量子效應(yīng)的同時提高其計算效率,才能使RPMD被廣泛應(yīng)用.近年來,美國加州理工學(xué)院的Korol等[46-47]進行了初步的嘗試,其提出的Cayley傳播子方法保持了RPMD的保辛性、時間可逆性、強穩(wěn)定性和遍歷性,同時支持較大的時間步長.

        1.2 運用Cayley傳播子加速RPMD

        原始的路徑積分分子動力學(xué)中,體系從t時刻隨時間演化到t+?t的運動方程為

        式中:A、B和O分別為自由珠串(free ring-polymer)隨時間演化、在外加勢能作用下隨時間演化和控溫;a為可調(diào)參數(shù).這里,當a=1時,為標準的對稱化RPMD時間演化方案OBABO;a=0時,時間演化記作BAOAB,即Liu等[48]和Zhang等[49]提出的中間控溫方案(middle thermostat)具有一定的數(shù)值穩(wěn)定性和計算效率優(yōu)勢.而O=0代表不控溫,對應(yīng)的是微正則系綜中的時間演化.

        上述傳播子中,A具有精確解析表達式:

        式中:ωP為串珠之間的諧振子頻率,同式(4).

        在數(shù)值計算中,需要引入近似RPMD傳播子M(?t)來計算exp(A?t).這里M(?t)需要滿足3個條件:

        (2)強穩(wěn)定性:對所有的ωj,P>0,且?t小于某不依賴于ωj,P的常數(shù),M(?t)都是強穩(wěn)定的辛矩陣;

        (3)時間反演對稱性:對于速度反演矩陣

        滿足

        標準的M(?t)使用Verlet積分方法,即

        此傳播子不滿足條件(2),導(dǎo)致標準的RPMD軌跡計算缺乏穩(wěn)定性.而Cayley傳播子考慮了隨機時間步長δt,并對其取統(tǒng)計平均:

        再將其對稱化,得到A?t的Cayley變換:

        這個傳播子可滿足全部的上述3條性質(zhì).Cayley傳播子方法就是利用上述cay(A?t)來代替標準的RPMD傳播子exp(A?t)進行計算,類比上述A、B的記號,記作C.

        考慮北京大學(xué)Liu等[48]提出的具有較高效率中間控溫方案(記作BAOAB),即將自由珠串含時演化拆分成2個部分,在這2個部分中間進行控溫,故Cayley傳播子可以進一步寫成相等的2個部分的乘積:

        并在上述2個半Cayley傳播子中間加入控溫步驟,實現(xiàn)正則系綜的Cayley傳播子PIMD/RPMD模擬,記為BCOCB.這是目前已知數(shù)值穩(wěn)定性最好的正則系綜Cayley傳播子實現(xiàn)方案[50].

        現(xiàn)在,Cayley傳播子可以對所有的串珠間諧振子頻率給出正確的量子分布,且在傳播中可以選取較大的時間步長[50-51].該方法目前已經(jīng)在一維非諧振子、液相水振動光譜計算中獲得了成功應(yīng)用.目前,Cayley傳播子已經(jīng)編寫進原有的RPMD計算程序,在H+H2中進行了測試,結(jié)果如表1所示.由表1可見,應(yīng)用Cayley傳播子后可以通過增大積分的時間步長,將計算速度提高了至少5倍而不改變計算結(jié)果.

        表1 Cayley傳播子RPMD(BCOCB)和標準RPMD(BAB+O)計算H+H2在300 K的速率結(jié)果對比Table 1 Comparion results of Cayley propogation sub-RPMD(BCOCB)and standard RPMD(BAB+O)calculated H+H2 at the rate of 300 K

        1.3 RPMD微正則系綜速率常數(shù)計算

        通過基于RPMD計算出的不同溫度下的大量熱反應(yīng)速率常數(shù),應(yīng)用最近開發(fā)的最大熵方法[53]可以計算得到微正則系綜速率常數(shù).由統(tǒng)計力學(xué)可知,熱速率常數(shù)k(β)和微正則系綜速率常數(shù)N(E)具有如下關(guān)系:

        式中:β=(kBT)?1.運用Laplace逆變換,有

        式中:

        直接應(yīng)用上述公式會有較大的統(tǒng)計誤差,而應(yīng)用最大熵方法可以得到有效的精確結(jié)果,因為計算是基于離散的數(shù)據(jù)點進行的.

        首先,將熱速率常數(shù)和微正則系綜速率常數(shù)的關(guān)系式改寫為矩陣形式:

        式中:S(v)為求得的v矩陣和某個模型λ(Ej)之間的信息熵;χ2(v)為似然函數(shù);Vreg(v)為一個約束函數(shù),用來保證N(Ei)∈[0,1],

        式中:ζ為一個可調(diào)參數(shù),通過測試計算確定.目前這套最大熵估計微正則系綜速率常數(shù)的方法在一維模型體系中得到了驗證,但尚未推廣到雙分子反應(yīng)的計算中.而且,要得到微正則系綜速率常數(shù),必須首先獲得至少數(shù)十個不同溫度下的k(T).標準RPMD的耗時十分巨大,限制了基于RPMD的微正則系綜速率常數(shù)的計算.但是,如果能成功提高RPMD計算效率,則雙分子反應(yīng)的微正則系綜速率常數(shù)的計算也能得以廣泛應(yīng)用.

        1.4 應(yīng)用非平衡RPMD計算反應(yīng)幾率和反應(yīng)截面

        QCT方法計算雙分子碰撞過程的計算量比標準RPMD更小,且相關(guān)物理量的計算較為簡單直接,因此將RPMD與QCT結(jié)合將會是一個計算效率較高的方法.近期研究結(jié)果表明,RPMD模擬不受零點能泄漏問題的影響[7,54-55],計算結(jié)果的準確率將遠高于標準QCT.Welsch等[56]基于非平衡RPMD,將RPMD初步應(yīng)用于QCT的相關(guān)計算,得到了與實驗結(jié)果高度符合的H2+Cu(111)和D2O+Ni(111)的表面散射反應(yīng)幾率[7];Marjollet等[8,57]運用這種思路計算了的Mu/H/D/Cl/F+H2散射截面.這些研究結(jié)果表明,RPMD-QCT計算具有廣泛應(yīng)用前景.但是,標準RPMD在QCT計算中需要很小的時間步長(約0.01 fs)以保證能量守恒,計算效率較低;另外,目前還沒有一套成熟的、基于RPMD的初態(tài)選取和產(chǎn)物態(tài)指認方案.上述2個問題導(dǎo)致了目前結(jié)合標準RPMD的QCT進展緩慢,而新型傳播子RPMD由于其計算效率高、穩(wěn)定性強,有望結(jié)合QCT發(fā)展新型計算方法,并快速獲得復(fù)雜氣相雙分子的熱速率常數(shù)和KIE等化學(xué)動力學(xué)中的關(guān)鍵物理量,擴大RPMD在分子反應(yīng)動力學(xué)理論研究中的應(yīng)用范圍.

        1.5 應(yīng)用機器學(xué)習(xí)電子結(jié)構(gòu)計算輔助RPMD計算

        復(fù)雜氣相雙分子反應(yīng)的勢能面開發(fā)較為困難.首先,應(yīng)用高精度量子化學(xué)方法計算大量高維空間的構(gòu)型要耗費大量機時;其次,構(gòu)型的選取和勢能面的構(gòu)建與完善也需要消耗大量時間.而基于RPMD-QCT,采用基于量子化學(xué)計算得到的梯度進行直接動力學(xué)計算,可以避免構(gòu)建勢能面.但是,氣相雙分子反應(yīng)的直接動力學(xué)需要較為精確的量子化學(xué)計算結(jié)果,目前常見的計算速度較快的量子化學(xué)方法(如密度泛函理論)無法滿足精度需求.近年來,機器學(xué)習(xí)電子結(jié)構(gòu)領(lǐng)域取得了顯著的進展,其中基于分子軌道的機器學(xué)習(xí)(molecular-orbital-based machine learning,MOB-ML)方法[58-59]可以在Hartree-Fock(H-F)的計算速度下,達到耦合簇(coupled cluster single double(triple),CCSD(T))的精度.目前,MOB-ML方法已經(jīng)可以計算閉殼層、開殼層和金屬配合物體系,且具有解析梯度[60],計算效率較高.可見,基于這套機器學(xué)習(xí)量子化學(xué)計算方法直接進行RPMD動力學(xué)模擬,可以擴大RPMD-QCT的應(yīng)用范圍.

        2 結(jié)論

        目前,國內(nèi)外應(yīng)用RPMD計算雙分子化學(xué)反應(yīng)速率的研究人員逐年增多,有大連化學(xué)物理研究所/復(fù)旦大學(xué)張東輝課題組、南京大學(xué)謝代前課題組、重慶大學(xué)李軍課題組、西北工業(yè)大學(xué)孟慶勇課題組、塞浦路斯研究所Yury等.本課題組于2012年首次將RPMD應(yīng)用于具有van der Waals勢阱的重-輕-重氣相雙分子反應(yīng)O+CH4的研究,獲得了與高精度量子波包動力學(xué)一致的結(jié)果.此后,本課題組繼續(xù)致力于RPMD方法在反應(yīng)動力學(xué)中的推廣,取得了一系列的成功,如計算無能壘反應(yīng)O+H2的速率常數(shù)等[61].2017年,Zuo等[62]基于MPI4Py開發(fā)了跨節(jié)點并行的RPMD速率計算版本,應(yīng)用于本課題組及重慶大學(xué)李軍、南京大學(xué)謝代前課題組的研究中.2020年,本課題組還將RPMD的應(yīng)用拓展到多通道漫游反應(yīng)[63].

        近年來,對RPMD直接發(fā)展的有加州理工學(xué)院Miller課題組和中科院北京化學(xué)研究所的史強研究員合作開發(fā)的非平衡RPMD、非絕熱RPMD以及Cayley傳播子RPM;瑞士洛桑聯(lián)邦理工學(xué)院的Ceriotti開發(fā)了適合RPMD的新型控溫方法;美國Stanford大學(xué)的Markland將RPMD應(yīng)用于凝聚相體系;北京大學(xué)的劉劍課題組在RPMD的控溫和非絕熱動力學(xué)方面做了重要工作等.

        綜上所述,針對復(fù)雜氣相雙分子化學(xué)反應(yīng)的反應(yīng)速率常數(shù)、KIE和反應(yīng)截面等反應(yīng)動力學(xué)關(guān)鍵物理量的計算,實時路徑積分RPMD方法體現(xiàn)了強大的應(yīng)用價值.但是,現(xiàn)有的理論研究方法仍然存在明顯不足.如何基于實時路徑積分方法,快速且高精度地計算化學(xué)反應(yīng)動力學(xué)中的速率常數(shù)等關(guān)鍵物理量,目前仍然是研究的熱點.未來的發(fā)展可能包括如下幾個方向:

        (1)運用新型傳播子加速標準的RPMD方法,使計算速度進一步加快;

        (2)開發(fā)RPMD-QCT計算方案,使得對化學(xué)反應(yīng)動力學(xué)關(guān)鍵物理量的計算不需要計算自由能,可以大大提高計算效率;

        (3)結(jié)合機器學(xué)習(xí)量子化學(xué)計算方法的直接動力學(xué),實現(xiàn)快速準確地計算上述物理量.

        上述3個方向的理論進展,將會使研究人員對雙分子氣相化學(xué)反應(yīng)的機理有更深層次的認識,也將擴大RPMD的應(yīng)用范圍,為雙分子反應(yīng)動力學(xué)研究中理論與實踐結(jié)合提供更有力的理論計算工具.

        猜你喜歡
        勢能課題組常數(shù)
        “動能和勢能”知識鞏固
        作 品:景觀設(shè)計
        ——《勢能》
        文化縱橫(2022年3期)2022-09-07 11:43:18
        陽城縣“耕心微寫”課題組
        “動能和勢能”知識鞏固
        “動能和勢能”隨堂練
        關(guān)于Landau常數(shù)和Euler-Mascheroni常數(shù)的漸近展開式以及Stirling級數(shù)的系數(shù)
        原科技大學(xué)新能源開發(fā)與應(yīng)用課題組介紹
        幾個常數(shù)項級數(shù)的和
        萬有引力常數(shù)的測量
        課題組成員
        支點(2015年11期)2015-11-16 10:25:03
        精品国产一区二区三区三级| 一区二区三区内射视频在线观看| 中文字幕一区韩国三级| 麻豆精品国产免费av影片| 在线麻豆精东9制片厂av影现网| 亚洲av无码专区在线播放中文| 亚洲日本天堂| 激情视频在线观看国产中文| 视频在线观看国产自拍| 97精品国产97久久久久久免费 | 亚洲黄色官网在线观看| 日本国产精品久久一线 | 东京无码熟妇人妻av在线网址| 国产精品三级一区二区按摩| 女同视频网站一区二区| 丝袜美腿av在线观看| 国产果冻豆传媒麻婆精东| 中文人妻av久久人妻18| 久久精品熟女不卡av高清| 日韩在线视频专区九区| 精品国品一二三产品区别在线观看| 男女性高爱潮免费观看| 国语憿情少妇无码av| 人妖在线一区二区三区| 中文字幕av一区二区三区人妻少妇| 91热这里只有精品| 日韩精品中文字幕免费人妻| 色视频网站一区二区三区| 日韩高清在线观看永久| 亚洲性69影视| 亚洲天堂色婷婷一区二区 | 精品女同一区二区三区免费播放| 久久午夜av一区二区三区| 色欲aⅴ亚洲情无码av| 国产av影片麻豆精品传媒| 亚洲精品中文字幕尤物综合| 少妇高潮在线精品观看| 无码av免费一区二区三区| 欧美一级欧美一级在线播放| 日本道免费一区日韩精品| 中文字幕一区二区三区视频|