基于實時路徑積分的雙分子化學反應動力學計算方法

李永樂范文斌任偉

(1.上海大學理學院量子與分子結構國際中心,上海 200444;2.上海大學理學院上海市高溫超導重點實驗室,上海 200444)

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

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

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

1 方法

1.1 路徑積分分子動力學

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

式中:h為普朗克常數;p和q為串珠的經典動量和坐標;為系統哈密頓算符;βP的定義為

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

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

從化學反應的量子理論[36]

出發,計算得到熱速率常數.式中:Qr為反應物的配分函數;cfs(t)為流關聯函數,

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

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

式中:R∞為設定的2個反應物相互作用可忽略的距離;μR為反應物的折合質量.

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

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

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

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

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

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

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

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

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

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

1.2 運用Cayley傳播子加速RPMD

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

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

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

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

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

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

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

滿足

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

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

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

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

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

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

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

表1 Cayley傳播子RPMD(BCOCB)和標準RPMD(BAB+O)計算H+H2在300 K的速率結果對比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微正則系綜速率常數計算

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

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

式中:

直接應用上述公式會有較大的統計誤差,而應用最大熵方法可以得到有效的精確結果,因為計算是基于離散的數據點進行的.

首先,將熱速率常數和微正則系綜速率常數的關系式改寫為矩陣形式:

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

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

1.4 應用非平衡RPMD計算反應幾率和反應截面

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

1.5 應用機器學習電子結構計算輔助RPMD計算

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

2 結論

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

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

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

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

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

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

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