C++H2 反應的動力學研究:基于一個新構建的勢能面*

李文濤 袁美玲 王杰敏

1)(洛陽師范學院物理與電子信息學院,洛陽 471934)

2)(濰坊科技學院通識學院,壽光 262700)

3)(棗莊學院人工智能學院,棗莊 277160)

采用一種明顯依賴于電子-電子距離的多參考組態相互作用方法,結合ACVQZ 基組進行了體系的從頭算計算.采用置換不變多項式神經網絡方法對從頭算得到的18222 個單點能進行了擬合.此外,將勢能面的特征與可得到的實驗結果和理論結果進行了比較.結果表明新構建的勢能面比以往的勢能面更加精確.基于新構建的勢能面,采用準經典軌線方法進行了C++H2 反應的動力學計算,報道了積分截面和微分截面等動力學信息,并與之前的理論結果進行了比較.動力學結果表明在反應過程中插入反應機理占據主導地位.

1 引言

自1941 年Douglas 和Herzberg[1]在星際介質中發現CH+離子的高豐度含量后,CH+離子受到了大量的實驗[2?6]和理論[7?29]研究.CH+離子是高活性的,在宇宙環境中很容易與電子、氫原子、氫分子等發生反應,從而被迅速破壞掉.因此,CH+離子在宇宙中的高豐度就變得難以解釋,這也是現代天體物理學中的一個難題.為解釋其在宇宙中高豐度的現象,在接下來的幾十年中,人們對CH+離子的形成和破壞機制進行了廣泛的研究.在實驗方面,1997 年,Hierl 等[2]采用高溫流動余熱裝置測量了C+(2P)和H2分子及其同位素D2分子反應在400—1300 K 內的熱速率常數.實驗結果表明,反應物分子的轉動能和平動能對反應的影響類似,而反應物分子的振動激發對反應的影響遠大于平動能和轉動能.與反應物的振動基態v=0比較,v=1 的振動激發,使H2分子速率常數增大約1000 倍,使D2分子速率常數增大約6000 倍.2011 年,Plasil 等[3]采用溫度可變的22 極離子阱和冷噴氫原子束,研究了CH+離子與H 原子的碰撞反應.他們報道了該反應在10—1000 K 溫度范圍內的速率常數,發現當溫度低于100 K 時,速率常數會出現一個迅速的下降.這說明非轉動的CH+離子在低溫情況下很難與H 原子發生碰撞反應.

在理論方面,大量的理論研究工作者采用構建高精度勢能面并進行動力學計算的方式來獲得指定電子態的速率常數.2005 年,Stoecklin 和Halvick[7]采用多項式擬合的方法構建了體系的基態勢能面,并基于該勢能面采用RIOSA-NIP 方法計算了H+CH+反應的反應概率、積分截面和速率常數.2007 年,Halvick 等[8]基于該勢能面采用準經典軌線方法(quasiclassical trajectory method,QCT)和相空間理論重新研究了H+CH+反應,并報道了該反應的速率常數.2013 年,Zanchet 等[9]基于該勢能面采用含時量子波包方法研究了C++H2反應,并報道了該反應量子態分辨的速率常數.2011 年,Warmbier 和Schneider[10]報道了體系的一個新勢能面.他們采用AVTZ和AVQZ 基組和多參考組態相互作用方法(multireference configuration interaction method,MRCI)進行了從頭算計算.此外,基于該勢能面還采用了QCT 方法和ABC 非含時的量子方法對H+CH+反應進行了動力學計算,并報道了積分截面和速率常數等動力學信息.2015 年,Li 等[11]采用MRCI方法和AV6Z 基組進行了高精度的從頭算計算,并采用多項式擬合的方法構建了體系高精度的勢能面.該勢能面得到廣泛的應用.Guo 等[12],Sundaram 等[13?15]和 Wu 等[16]采用該勢能面對H+CH+反應及其逆反應進行了動力學計算并報道了積分截面、速率常數等動力學信息.最新的勢能面是由Guo 等[17]構建的,在從頭算計算中,他們把基組外推到了完全基組極限,同樣采用多項式擬合的方法構建了新的勢能面.

2 勢能面

2.1 從頭算計算

本文采用Molpro 程序包進行從頭算計算,并且對所有的構型都采用了Cs對稱性.在從頭算計算中,采用了一種明顯依賴于電子-電子距離的內收縮多參考組態相互作用方法(MRCI-F12)[30].此外,對C+離子和H 原子采用了ACVQZ 相關一致基組[31,32].對于所有構型,首先執行包含4 個等權重電子態(3 個2A',1 個2A'')的態平均完全活化空間自洽場(state average-complete active space self-consistent field,SA-CASSCF)[33,34]計算,然后利用SA-CASSCF 計算提供的參考波函數進行了MRCI-F12 計算.為了補償高階截斷誤差,計算中考慮了Davidson 修正.在計算中采用雅克比坐標來排布空間格點(RQ,RHH,α).其中,RQ是C+離子到H2分子的中心距離,在0.6—25 Bohr 的范圍內放置55 個非均勻格點;RHH是H2分子的鍵長,在0.6—25 Bohr 的范圍內放置56 個非均勻格點;α是RQ和RHH之間的夾角,其變化范圍是1°到89°,間隔為15°.在擬合過程中,刪掉了能量高于基態最低能量10 eV 的從頭算能量點.最終,在擬合過程中采用了18222 個從頭算能量點.

2.2 勢能面擬合

置換不變多項式神經網絡(permutive invariant polynomial-neuronal network,PIP-NN)[35,36]方法已經廣泛地應用于勢能面構建,例如H +HBr[37]和K+H2反應[38]等.在此,簡要地介紹一下PIP-NN 方法.眾所周知,神經網絡方法由3 個部分組成,即輸入端、隱藏層、輸出端.對于輸入端,為了避免在邊界處一階導數不連續,采用PIP方法對雅克比坐標進行如下變換:

其中,Ntot是從頭算能量點的個數;是連接權重;是相位;f1和f2是訓練函數,其數學表達式為

過擬合問題一直是神經網絡方法中的一個棘手問題,為了避免過擬合的發生,將輸入數據隨機分為3 個部分,即訓練部分(90%)、測試部分(5%)和驗證部分(5%).采用均方根誤差(root mean square error,RMSE)來測試勢能面的精確度,它的表達式為

其中,N是總的能量點個數,是擬合值,是從頭算能量.對于本文的體系,擬合RMSE僅為1.27 meV.

2.3 勢能面特征

測試勢能面的精確程度的一個有效途徑就是將計算得到的光譜常數與實驗進行比對.當第3 個原子遠離另外兩個原子中心距離為30 Bohr 時,改變另外兩個原子之間的間距,就可以在勢能面上得到CH+和H2分子的勢能曲線.將該曲線用最小二乘法擬合得到光譜常數,與實驗結果和其他理論結果一并列入到了表1 中.

如表1 所列,本文擬合得到的CH+和H2分子的光譜數據與實驗值十分吻合.對于CH+離子,平衡核間距re、諧振頻率ωe、離解能De與實驗值的差距僅為0.001 Bohr,2.75 cm–1,3 meV.對于H2分子,平衡核間距re與實驗結果一致,諧振頻率ωe和離解能De與實驗值的差距為6.42 cm–1和5 meV.實驗與理論的微小差距說明新構建的勢能面很好地描述了當第3 個原子遠離另外兩個原子時的勢能曲線.

表1 CH+和H2 分子的光譜常數Table 1. Spectroscopic constants of the CH+ and H2 molecules.

圖1 和圖2 展示了這個新構建的 CH+2勢能面的主要地形特征.由圖1 可知,在研究涉及的空間范圍內等勢線正確并光滑,其等勢線對應反應過程中一些主要的駐點.圖1(a)描述了鍵角∠HC+H=136.68°時,RCH+鍵伸縮的等勢線圖.很明顯,圖1(a)中有一個很深的勢阱,對應整個勢能面的最低能量點,對應的結構為RCH+=2.08 Bohr,RHH=3.866 Bohr,相比C++H2通道的能量,勢阱的深度約為4.30 eV.圖1(b)描述了C+離子垂直接近H2分子中心點過程的等勢線.可看到兩個駐點分別對應于局域最小值和整體最小值.其中局域最小值的構型為RCH+=2.602 Bohr,RHH=1.695 Bohr.圖2(a)給出了當CH+離子處于平衡核間距的構型下,H原子圍繞CH+離子運動的等勢線圖.在C+離子一側有一個明顯的深勢阱.當H 原子在∠HQC+=0°—60°范圍接近CH+離子時,很容易被該勢阱捕獲,從而形成長壽命的中間產物.圖2(b)給出了當H2分子處于平衡鍵長時,C+離子圍繞H2分子的等勢線圖.在90°附近有一個明顯的深勢阱.可以預見:當C+離子接近H2分子時,很可能被該勢阱吸引,從而進入T 型結構;對于C++H2→CH++H 反應,插入反應機理對該反應起著重要作用.

圖1 (a)鍵角在136.68°時,化學鍵伸縮的等勢線圖(等勢線的起點為–9.1 eV,間隔為0.4 eV);(b) C+離子以C2v 對稱性接近H2 分子中心的等勢線圖(等勢線的起點為–9.1 eV,間隔為0.37 eV)Fig.1.(a) Contour plot for chemical bond stretching,in which the angle is fixed at 136.68° (Contours starting at–9.1 eV and equally spaced by 0.4 eV);(b) contour plot for the C+ ion approach to the midpoint of H2 molecule in the C2v symmetry (Contours starting at–9.1 eV and equally spaced by 0.37 eV).

圖2 (a) 當 RCH+=2.136 Bohr 時,H 原子繞CH+離子運動的等勢線;(b) 當 RHH=1.401 Bohr 時,C+離子圍繞H2 分子運動的等勢線.Fig.2.(a) Contour plot for the H atom moves around CH+ ion at the bond distance RCH+=2.136 Bohr;(b) contour plot for C+ion moves around H2 molecule at its equilibrium geometry RHH=1.401 Bohr.

圖3 勢能面不同角度的最小能量路徑以及來自文獻[11]的理論結果Fig.3.The minimum energy paths of the potential energy surface at different angles along with the results obtained from Ref.[11].

3 動力學

為了進一步驗證勢能面的準確性,基于新構建的勢能面,采用QCT[42]對C++H2反應以0.05 eV的步長,在0.4—1.0 eV 的碰撞能范圍內進行了動力學計算.當入射C+原子與H2分子的質心的距離小于8 ?時,軌線開始計算,并且對于每一個能量隨機采樣1×105條軌線.在軌線的計算中,為了確保數值的穩定性,采用的積分步長為0.02 fs.在正式開始計算之前,對于每一個能量點,需要對最大碰撞參數bmax進行調節.在本文的計算中,采取測試的軌線數目為5×103條,在測試過程中逐漸增大bmax直至沒有反應軌線出現.在積分截面和微分截面的計算中,碰撞參數b在0 和最大碰撞參數bmax之間隨機取值.當產物分子之間的距離大于8 ?時,軌線終止計算.積分截面和微分截面可以通過下列公式得到:

其中,Rtra是反應軌線數目,Ntra是總軌線數目,Δθ是角度步長,Nθ是角方向格點數目.

圖4 給出了C++H2→H+CH+反應的積分截面.為了便于與之前的理論結果進行比較,Guo等[12]和Herráez-Aguilar 等[18]的結果也列入到了圖中.如圖4 所示,積分截面的大小隨著碰撞能的增加而增大,Guo 等[12]和Herráez-Aguilar 等[18]的結果與本文結果具有相同的趨勢.盡管如此,本文的結果與Guo 等[12]和Herráez-Aguilar 等[18]的結果之間還存在著明顯的不同.Herráez-Aguilar等[18]的結果是基于Warmbier 和Schneider[10]的勢能面,采用QCT 得到的.與本文的結果比較,在研究的碰撞能范圍內,本文的結果明顯大于Herráez-Aguilar 等[18]的結果,這是因為在計算中采用了不同的勢能面.Guo 等[12]的結果是基于Li 等[11]的勢能面,采用含時量子波包方法計算得到的.然而在他們的計算中,對角動量J在z軸方向上的投影進行了截斷(Kblock=20).對于固定的J值,完全的科里奧利耦合計算在z軸方向上的投影應為J+1,這表明Guo 等[12]的結果在J≥20 時就可能得不到收斂的動力學結果.因此,Guo 等[12]的結果很可能在低能處是收斂的,而對于較高的碰撞能,他們的結果很可能不收斂.本文采用的QCT方法不能考慮低能處的量子效應,因此,在低能處,本文的結果與Guo 等[12]的結果有一些差距.在高能處,量子效應變得不那么明顯,量子結果應與QCT 結果取得類似的結果.然而,Guo 等[12]的結果與本文的QCT 結果依然存在著較大的差距,這可能是因為Guo 等[12]截斷了J在z軸方向上的投影,從而導致量子結果不收斂.

圖4 C++H2 反應碰撞在0.4—1.0 eV 范圍內的積分截面以及文獻[12,18]的結果Fig.4.Integral cross section of C++H2 reaction in the collision energy range from 0.4 to 1.0 eV along with the values obtained from Ref.[12,18].

圖5 給出了C++H2→ H+CH+反應在不同碰撞能下的微分截面.從圖5 不難看出,對于所有能量點,散射信號都集中在0°和180°附近,并且前后向對稱.這說明,在反應過程中插入反應機理占據主導地位.反應路徑上的深勢阱(圖3)支持長壽命的束縛態和準束縛態,從而使得產物可以到達任何可接觸到的通道.

圖5 C++H2 反應若干碰撞能下的微分截面Fig.5.Differential cross sections of the C++H2 reaction at several collision energies.

4 結論