C++H2 反應(yīng)的動(dòng)力學(xué)研究:基于一個(gè)新構(gòu)建的勢(shì)能面*

李文濤 袁美玲 王杰敏

1)(洛陽(yáng)師范學(xué)院物理與電子信息學(xué)院,洛陽(yáng) 471934)

2)(濰坊科技學(xué)院通識(shí)學(xué)院,壽光 262700)

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

采用一種明顯依賴于電子-電子距離的多參考組態(tài)相互作用方法,結(jié)合ACVQZ 基組進(jìn)行了體系的從頭算計(jì)算.采用置換不變多項(xiàng)式神經(jīng)網(wǎng)絡(luò)方法對(duì)從頭算得到的18222 個(gè)單點(diǎn)能進(jìn)行了擬合.此外,將勢(shì)能面的特征與可得到的實(shí)驗(yàn)結(jié)果和理論結(jié)果進(jìn)行了比較.結(jié)果表明新構(gòu)建的勢(shì)能面比以往的勢(shì)能面更加精確.基于新構(gòu)建的勢(shì)能面,采用準(zhǔn)經(jīng)典軌線方法進(jìn)行了C++H2 反應(yīng)的動(dòng)力學(xué)計(jì)算,報(bào)道了積分截面和微分截面等動(dòng)力學(xué)信息,并與之前的理論結(jié)果進(jìn)行了比較.動(dòng)力學(xué)結(jié)果表明在反應(yīng)過(guò)程中插入反應(yīng)機(jī)理占據(jù)主導(dǎo)地位.

1 引言

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

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

2 勢(shì)能面

2.1 從頭算計(jì)算

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

2.2 勢(shì)能面擬合

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

其中,Ntot是從頭算能量點(diǎn)的個(gè)數(shù);是連接權(quán)重;是相位;f1和f2是訓(xùn)練函數(shù),其數(shù)學(xué)表達(dá)式為

過(guò)擬合問(wèn)題一直是神經(jīng)網(wǎng)絡(luò)方法中的一個(gè)棘手問(wèn)題,為了避免過(guò)擬合的發(fā)生,將輸入數(shù)據(jù)隨機(jī)分為3 個(gè)部分,即訓(xùn)練部分(90%)、測(cè)試部分(5%)和驗(yàn)證部分(5%).采用均方根誤差(root mean square error,RMSE)來(lái)測(cè)試勢(shì)能面的精確度,它的表達(dá)式為

其中,N是總的能量點(diǎn)個(gè)數(shù),是擬合值,是從頭算能量.對(duì)于本文的體系,擬合RMSE僅為1.27 meV.

2.3 勢(shì)能面特征

測(cè)試勢(shì)能面的精確程度的一個(gè)有效途徑就是將計(jì)算得到的光譜常數(shù)與實(shí)驗(yàn)進(jìn)行比對(duì).當(dāng)?shù)? 個(gè)原子遠(yuǎn)離另外兩個(gè)原子中心距離為30 Bohr 時(shí),改變另外兩個(gè)原子之間的間距,就可以在勢(shì)能面上得到CH+和H2分子的勢(shì)能曲線.將該曲線用最小二乘法擬合得到光譜常數(shù),與實(shí)驗(yàn)結(jié)果和其他理論結(jié)果一并列入到了表1 中.

如表1 所列,本文擬合得到的CH+和H2分子的光譜數(shù)據(jù)與實(shí)驗(yàn)值十分吻合.對(duì)于CH+離子,平衡核間距re、諧振頻率ωe、離解能De與實(shí)驗(yàn)值的差距僅為0.001 Bohr,2.75 cm–1,3 meV.對(duì)于H2分子,平衡核間距re與實(shí)驗(yàn)結(jié)果一致,諧振頻率ωe和離解能De與實(shí)驗(yàn)值的差距為6.42 cm–1和5 meV.實(shí)驗(yàn)與理論的微小差距說(shuō)明新構(gòu)建的勢(shì)能面很好地描述了當(dāng)?shù)? 個(gè)原子遠(yuǎn)離另外兩個(gè)原子時(shí)的勢(shì)能曲線.

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

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

圖1 (a)鍵角在136.68°時(shí),化學(xué)鍵伸縮的等勢(shì)線圖(等勢(shì)線的起點(diǎn)為–9.1 eV,間隔為0.4 eV);(b) C+離子以C2v 對(duì)稱性接近H2 分子中心的等勢(shì)線圖(等勢(shì)線的起點(diǎn)為–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) 當(dāng) RCH+=2.136 Bohr 時(shí),H 原子繞CH+離子運(yùn)動(dòng)的等勢(shì)線;(b) 當(dāng) RHH=1.401 Bohr 時(shí),C+離子圍繞H2 分子運(yùn)動(dòng)的等勢(shì)線.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 勢(shì)能面不同角度的最小能量路徑以及來(lái)自文獻(xiàn)[11]的理論結(jié)果Fig.3.The minimum energy paths of the potential energy surface at different angles along with the results obtained from Ref.[11].

3 動(dòng)力學(xué)

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

其中,Rtra是反應(yīng)軌線數(shù)目,Ntra是總軌線數(shù)目,Δθ是角度步長(zhǎng),Nθ是角方向格點(diǎn)數(shù)目.

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

圖4 C++H2 反應(yīng)碰撞在0.4—1.0 eV 范圍內(nèi)的積分截面以及文獻(xiàn)[12,18]的結(jié)果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+反應(yīng)在不同碰撞能下的微分截面.從圖5 不難看出,對(duì)于所有能量點(diǎn),散射信號(hào)都集中在0°和180°附近,并且前后向?qū)ΨQ.這說(shuō)明,在反應(yīng)過(guò)程中插入反應(yīng)機(jī)理占據(jù)主導(dǎo)地位.反應(yīng)路徑上的深勢(shì)阱(圖3)支持長(zhǎng)壽命的束縛態(tài)和準(zhǔn)束縛態(tài),從而使得產(chǎn)物可以到達(dá)任何可接觸到的通道.

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

4 結(jié)論