采用量子含時(shí)波包方法研究H/D+Li2→ LiH/LiD+Li反應(yīng)?

李文濤于文濤3)姚明海

1)(渤海大學(xué)基礎(chǔ)教研部,錦州 121000)

2)(中國科學(xué)院大連化學(xué)物理研究所,分子反應(yīng)動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室,大連 116023)

3)(東風(fēng)縣第三中學(xué),遼源 136300)

1 引 言

近數(shù)十年來,計(jì)算機(jī)科學(xué)與數(shù)值方法取得的重大進(jìn)展使得人們可以用量子的方法對復(fù)雜的三原子反應(yīng)、多原子碰撞反應(yīng)以及三原子三體重組反應(yīng)等[1?4]進(jìn)行精確的理論模擬研究.HLi2體系作為一個(gè)典型的且具有深勢阱的復(fù)雜放熱反應(yīng),在以往的理論研究中常采用近似的方法對其進(jìn)行理論研究.計(jì)算機(jī)科學(xué)和數(shù)值方法的進(jìn)展使得對HLi2體系進(jìn)行精確的理論計(jì)算成為可能.此外,由于HLi2體系是含有最少電子的最小混合團(tuán)簇體系并且也是最穩(wěn)定的中性三聚物之一,這些性質(zhì)使得HLi2成為研究金屬的化學(xué)吸附現(xiàn)象和改變合金性質(zhì)的最佳體系之一.這些特性也使得HLi2體系在實(shí)驗(yàn)和理論上都得到了廣泛的研究[5?17].

實(shí)驗(yàn)上,Wu和Ihle[5]通過質(zhì)譜測量液體鋰中稀釋的氫溶液證明了HLi2分子的穩(wěn)定存在并報(bào)道了HLi2分子的離解能(89.7±5.0)kcal/mol.基于一個(gè)優(yōu)化的勢能面,Kim和Herschbach[6]更正了勢能深度的最小值為(87.9±3.0)kcal/mol.1993年,Vezin等[7]報(bào)道了HLi2分子實(shí)驗(yàn)光學(xué)光譜,并通過模擬第一振動帶的轉(zhuǎn)動結(jié)構(gòu)得到了HLi2分子基態(tài)與B激發(fā)態(tài)的幾何結(jié)構(gòu).為了解反應(yīng)機(jī)理,進(jìn)行動力學(xué)散射計(jì)算的理論模擬是十分必要的.而勢能面是動力學(xué)理論計(jì)算的基石,因此在過去的幾十年中也有大量的理論工作集中在該體系的勢能面構(gòu)建上.

采用自洽場波函和全組態(tài)相互作用方法并結(jié)合一個(gè)收縮高斯基組(H(5s1p/3s1p),Li(8s3p/4s3p)),Siegbahn和Schaefer[8]在1975年構(gòu)建了第一個(gè)HLi2分子體系的勢能面.然而計(jì)算中采用的基組過小,使得其構(gòu)建的勢能面有很大的改善空間.1997年,Yan等[9]使用285個(gè)從頭算能量點(diǎn)構(gòu)建了一個(gè)新的勢能面.在他們的工作中,采用多參考組態(tài)相互作用單雙激發(fā)方法和6-311G(2df,2pd)基組進(jìn)行了從頭算計(jì)算并確定了HLi2分子體系的放熱能為32.85 kcal/mol.2010年,Maniero等[10]對394個(gè)從頭算能量點(diǎn)進(jìn)行八階的多項(xiàng)式擬合,從而構(gòu)建了一個(gè)新的勢能面.通過對勢能面的地形特征分析表明,H+Li2反應(yīng)是一個(gè)沒有閾值的放熱反應(yīng),放熱能量約為34 kcal/mol.采用雙重多體展開公式對3726個(gè)從頭算能量點(diǎn)進(jìn)行擬合,Song等[11]構(gòu)建了一個(gè)HLi2分子體系基態(tài)的全維勢能面.在計(jì)算過程中,采用多參考組態(tài)相互作用方法,對于Li原子采用cc-pV5Z基組,對于H原子采用aug-cc-pV5Z基組進(jìn)行了從頭算計(jì)算.與Maniero等得到結(jié)論類似,明確H+Li2反應(yīng)是一個(gè)放熱反應(yīng)并且放熱能為33.668 kcal/mol.最新的勢能面是2017年Yuan等[12]采用置換不變多項(xiàng)式神經(jīng)網(wǎng)絡(luò)方法對大約30000個(gè)從頭算能量點(diǎn)進(jìn)行擬合得到的.在他們的工作中,采用MRCI-F12方法結(jié)合aug-cc-pVTZ基組進(jìn)行了從頭算計(jì)算.由于H+Li2反應(yīng)是一個(gè)沒有閾值的放熱反應(yīng),因此長程相互作用勢可能會對低能的反應(yīng)碰撞產(chǎn)生大的影響.該勢能面采用比較密集的從頭算能量點(diǎn)在反應(yīng)漸進(jìn)區(qū),以便于得到比較準(zhǔn)確的長程相互作用勢.

在動力學(xué)方面,準(zhǔn)經(jīng)典軌線方法(QCT)和量子方法(QM)經(jīng)常被用來進(jìn)行動力學(xué)計(jì)算.Kim等[13]采用QCT方法分別計(jì)算了平動能為0.5,2和10 kcal/mol時(shí)H+Li2反應(yīng)的積分截面.2011年,Vila等[14]采用非含時(shí)的量ABC程序在總角動量J=0時(shí)對H+Li2→LiH+Li反應(yīng)進(jìn)行了動力學(xué)模擬,采用J-shifting方法對于J>0進(jìn)行了計(jì)算.理論模擬結(jié)果表明,反應(yīng)物的振轉(zhuǎn)激發(fā)對反應(yīng)起阻礙作用.此外,在0—4000 K的溫度范圍內(nèi)計(jì)算了H+Li2反應(yīng)體系的速率常數(shù).2014年,Cunha等[15]采用QCT方法研究了H+Li2體系的同位素效應(yīng)并報(bào)道了平動能范圍在0—9 kcal/mol范圍內(nèi)H/D/Mu+Li2反應(yīng)的積分截面.2015年,Gao等[16]采用含時(shí)量子波包方法和二階分裂算符進(jìn)行了H+Li2體系的動力學(xué)計(jì)算,在0—0.4 eV的碰撞能范圍內(nèi)報(bào)道了該體系的反應(yīng)概率和積分截面.此外,他們還發(fā)現(xiàn)Li2分子的振動激發(fā)在低碰撞能范圍內(nèi)阻礙反應(yīng)的進(jìn)行.2016年,為了明確科里奧利耦合在反應(yīng)過程的重要性,Gao等[17]重新對該體系進(jìn)行了動力學(xué)計(jì)算,結(jié)果表明科里奧利耦合在反應(yīng)過程中起著重要的作用.然而,該體系是一個(gè)具有深勢阱的反應(yīng)體系,要得到收斂的結(jié)果比較耗時(shí).因此,他們在計(jì)算中截?cái)嗔丝偨莿恿吭隗w固定坐標(biāo)軸上投影的數(shù)目(k=1[16]和k=6[17]).

綜上所述,對H+Li2反應(yīng)體系已有充分的理論研究.然而,由于該體系是一個(gè)高放熱并且在反應(yīng)路徑上具有一個(gè)深勢阱的體系,要得到收斂的結(jié)果比較耗時(shí),因此以往的理論研究中往往采用QCT方法或者近似的量子方法.此外,同位素效應(yīng)是在不對體系進(jìn)行大的改變的基礎(chǔ)上進(jìn)一步探究反應(yīng)機(jī)理的重要手段,在物理化學(xué)、原子與分子物理領(lǐng)域已有同位素對比的廣泛研究[18?22].而對于H+Li2反應(yīng)體系,到目前為止尚未見在量子理論水平進(jìn)行同位素對比的研究.本文的目標(biāo)主要集中在兩個(gè)方面:1)以往的理論研究中對總角動量在體固定坐標(biāo)軸上投影(k-block)的數(shù)目進(jìn)行了截?cái)?而與科里奧利耦合相關(guān)的反應(yīng)概率、積分截面等動力學(xué)性質(zhì)顯著地依賴于k-block的數(shù)目,尤其是總角動量J比較大時(shí);收斂的動力學(xué)結(jié)果有助于理解動力學(xué)反應(yīng)機(jī)理,本文進(jìn)行包含所有k-block數(shù)目的動力學(xué)計(jì)算并與之前的理論研究進(jìn)行對比;2)將H原子取代為較重的D原子,在量子水平上進(jìn)行同位素變換的研究,進(jìn)一步研究動力學(xué)反應(yīng)機(jī)理.

2 動力學(xué)計(jì)算方法

量子含時(shí)波包量子方法在動力學(xué)計(jì)算中已廣泛使用[23,24],本文只做一些簡單的介紹,詳細(xì)的信息可以參考文獻(xiàn)[25—28].對于給定的總角動量J,H+Li2反應(yīng)體系的哈密頓量在體坐標(biāo)系下可以寫為

其中G(R)是高斯波包,?v0j0是Li2分子的振轉(zhuǎn)方程,(M,j0,l0)是體坐標(biāo)下體系的初始振轉(zhuǎn)態(tài),ε是體系的宇稱.在傳播過程中,采用快速傅里葉變換方法[29]進(jìn)行徑向動能算符的計(jì)算,采用離散表象變換方法[30]進(jìn)行勢能算符的計(jì)算,并采用二階分裂算符傳播構(gòu)建好的指定初始態(tài)的量子波包.為了避免波包在邊界處反彈造成二次傳播,在計(jì)算中采用如下形式的吸收勢:

其中Ca,Cb是吸收勢的強(qiáng)度;?t是時(shí)間步長;ra,rb,rend是吸收勢對應(yīng)的位置.

通過取代旋量表示下的散射矩陣可以得到態(tài)分辨的積分截面

和態(tài)分辨的微分截面

其中是散射矩陣,d?是方位角,?是散射角,是魏格納轉(zhuǎn)動矩陣.

3 結(jié)果與討論

在計(jì)算中采用Song等[11]報(bào)道的勢能面,該勢能面的主要特征如圖1所示.由圖1可知,H+Li2反應(yīng)是一個(gè)放熱反應(yīng),其放熱能約為1.46 eV(考慮零點(diǎn)能),并且在反應(yīng)路徑上有一個(gè)深勢阱,相對于反應(yīng)物漸進(jìn)區(qū),其深度約為2.6 eV.

3.1 反應(yīng)概率

對于H+Li2反應(yīng),在總角動量J=0—60的范圍內(nèi)進(jìn)行了動力學(xué)計(jì)算.對于其同位素D+Li2反應(yīng),計(jì)算的總角動量范圍為0—70.為了得到收斂的動力學(xué)結(jié)果,在總角動量J=0時(shí),本文進(jìn)行了大量的動力學(xué)參數(shù)測試.得到的最優(yōu)參數(shù)列于表1.對于所有總角動量J>0的情況,采用相同的計(jì)算參數(shù).此外,在計(jì)算中本文考慮了與k-block密切相關(guān)的科里奧利耦合效應(yīng).在計(jì)算中對于每一個(gè)總角動量J,采用的k-block數(shù)目為(J+1)和100之間的最小值,這樣可以保證總角動量所有的投影都包含在計(jì)算中,以便得到收斂的結(jié)果.

圖1 H+Li2反應(yīng)5?—180?的最小反應(yīng)能量路徑Fig.1. The minimum energy reaction path from 5?–180? for H+Li2reaction.

圖2 在0—0.4 eV的碰撞能范圍內(nèi),H+Li2反應(yīng)總角動量J=0,10,20,30,40,50的反應(yīng)概率以及文獻(xiàn)[14,17]的結(jié)果Fig.2.The total reaction probabilities for the total angular momentum J=0,10,20,30,40,50 of the H+Li2reaction in the collision energy range from 0–0.4 eV as well as the results from Ref.[14,17].

總角動量J=0,10,20,30,40,50的反應(yīng)概率見圖2.如圖2所示,在J=0的反應(yīng)概率上沒有閾值出現(xiàn)并且在研究的碰撞能范圍內(nèi)有一些共振峰,這些共振峰的產(chǎn)生是因?yàn)樵诜磻?yīng)路徑上有一個(gè)能產(chǎn)生大量束縛態(tài)和準(zhǔn)束縛態(tài)的深勢阱;伴隨著碰撞能的增加,共振峰的強(qiáng)度變?nèi)?共振峰之間的間距變寬.這表明隨著碰撞能的增加,碰撞絡(luò)合物的壽命開始變短.當(dāng)總角動量J變大以后,反應(yīng)閾值開始出現(xiàn)并且隨著J的增大而增大.反應(yīng)閾值的出現(xiàn)是因?yàn)殡x心勢的變大,離心勢隨著J的增大而變大.

為了與之前的理論結(jié)果進(jìn)行比較,圖2同時(shí)給出了文獻(xiàn)[14]和文獻(xiàn)[17]報(bào)道的理論結(jié)果.如圖2所示,對于J=0,本文結(jié)果與文獻(xiàn)[17]的結(jié)果符合得很好,與文獻(xiàn)[14]在低能部分符合得比較好,但是對于高能部分則有較大的差距.這可能是因?yàn)樵谟?jì)算中采用不同勢能面造成的:文獻(xiàn)[14]采用Maniero等[10]構(gòu)建的勢能面僅包含394個(gè)從頭算能量點(diǎn),并且采用的基組也有進(jìn)一步提升的空間;而文獻(xiàn)[17]采用Song等[11]構(gòu)建的勢能面從頭算能量點(diǎn)的數(shù)目幾乎是Maniero等[10]的10倍,此外,勢阱的深度也比Maniero等[10]構(gòu)建的勢能面更深(1 kcal/mol)也更加接近實(shí)驗(yàn)值.勢能面的地形結(jié)構(gòu)往往會對動力學(xué)結(jié)果產(chǎn)生較大的影響.在低碰撞能時(shí),由于平動能低,勢阱深度的微弱差距對中間絡(luò)合物的壽命產(chǎn)生的影響相對較小,所以在低能時(shí),反應(yīng)概率符合得比較好.當(dāng)平動能增加時(shí),勢阱深度的差異會對中間絡(luò)合物的壽命產(chǎn)生較大的影響,所以在圖2中不難看出隨著碰撞能的增加,文獻(xiàn)[14]的反應(yīng)概率與本文反應(yīng)概率的差距在逐漸變大.

對于J>0,本文的結(jié)果與文獻(xiàn)[17]的結(jié)果有較大的差距,有時(shí)甚至不能得到相同的閾值.分析認(rèn)為,隨著碰撞能和離心勢的增加,反應(yīng)概率會逐漸變小.這是一個(gè)典型的放熱反應(yīng)特征.因?yàn)殡S著能量的變大,H原子或D原子的速度逐漸變快,以致于Li2分子來不及與快速運(yùn)動的H原子或者D原子形成HLi鍵或者DLi鍵,從而降低反應(yīng)概率.而文獻(xiàn)[17]的結(jié)果明顯有悖于這樣的趨勢.例如,J=10的反應(yīng)概率小于J=20的反應(yīng)概率,這明顯是不合理的.此外,本文結(jié)果與文獻(xiàn)[17]結(jié)果的差距也隨著碰撞能的增加變得越來越明顯.分析認(rèn)為,這樣的結(jié)果或許是因?yàn)槲墨I(xiàn)[17]為了節(jié)省計(jì)算量對k-block數(shù)目進(jìn)行了截?cái)?k=6).在文獻(xiàn)[17]中,Gao等[17]對J=30,50時(shí)分別對k=1—10,15進(jìn)行了測試,發(fā)現(xiàn)k=6的結(jié)果與k=15的結(jié)果比較接近,從而認(rèn)為k=6對于多數(shù)J應(yīng)該是收斂的,并將k=6應(yīng)用到所有J的計(jì)算當(dāng)中.分析認(rèn)為這樣的做法是值得商榷的,因?yàn)閗=15無論是對J=30或J=50得到的結(jié)果都不是真正意義上的收斂結(jié)果.如果k=15得到的結(jié)果是不收斂的,那么文獻(xiàn)[17]進(jìn)行的比較將變得沒有意義.因此本文以J=30為例對k=1,5,10,15,20,25的結(jié)果進(jìn)行了計(jì)算和比較,結(jié)果如圖3所示.

從圖3中不難看出,當(dāng)k=15時(shí),得到的結(jié)果是不收斂的.從圖3(b)可以看出,反應(yīng)概率在高碰撞能范圍的巨大差距主要是因?yàn)槲墨I(xiàn)[17]沒有把所有的k-block值考慮在計(jì)算中.

為了解同位素效應(yīng)的影響,圖4給出了碰撞能0—0.4 eV的范圍內(nèi)J=0,60和70的反應(yīng)概率.由圖4可知,當(dāng)H原子被D原子取代之后,反應(yīng)概率得到了提升,這或許是因?yàn)镈Li鍵的鍵能弱于HLi鍵的鍵能,DLi鍵比HLi鍵更加容易形成與斷裂.圖4還表明,要獲得收斂的積分截面和微分截面需要大量的分波計(jì)算.本文對于H+Li2反應(yīng)的最大總角動量J=60,表明當(dāng)碰撞能低于0.23 eV時(shí).獲得的積分截面和微分截面是收斂的.對于D+Li2反應(yīng),由于H原子被D原子取代,導(dǎo)致離心勢變小.本文計(jì)算的最大總角動量J=70在碰撞能低于0.15 eV時(shí),獲得的動力學(xué)結(jié)果是收斂的.

圖3 當(dāng)總角動量J=30時(shí),不同k-block值的反應(yīng)概率(a)和文獻(xiàn)[17]結(jié)果的比較(b)Fig.3.The reaction probabilities of the different k-block values(a)and compared with the values obtained from Ref.[17](b).

圖4 H/D+Li2反應(yīng)在總角動量J=0,60和70時(shí)的反應(yīng)概率Fig.4.The total reaction probabilities of the H/D+Li2reaction as the total angular momentum J=0,60,and 70,respectively.

3.2 積分截面與微分截面

如圖5所示,H/D+Li2反應(yīng)的積分截面在0—0.2 eV的范圍內(nèi)隨著碰撞能的增加而逐漸降低,這是一個(gè)典型的放熱反應(yīng)特征.由圖5可知,同位素取代增大了反應(yīng)的積分截面,尤其是在低碰撞能的時(shí)候.這是可以預(yù)料的,因?yàn)橥凰厝〈矊?dǎo)致了反應(yīng)概率的增加.

圖5 H/D+Li2反應(yīng)在碰撞能0—0.2 eV范圍內(nèi)的積分截面Fig.5.The total integral cross sections of H/D+Li2 reaction in the collision energy range from 0 to 0.2 eV.

圖6 H/D+Li2反應(yīng)在0.1和0.15 eV時(shí)的微分截面Fig.6. The differential cross sections of H/D+Li2 reaction at 0.1 and 0.15 eV,respectively.

當(dāng)碰撞能為0.1和0.15 eV時(shí),H/D+Li2反應(yīng)體系的微分截面如圖6所示.可見同位素取代同樣導(dǎo)致了微分截面的增加,但并未對反應(yīng)機(jī)理產(chǎn)生影響.無論是在0.1 eV還是在0.15 eV時(shí),微分截面都是前后對稱的,這表明在整個(gè)反應(yīng)過程中插入反應(yīng)機(jī)理占據(jù)主導(dǎo)地位.

4 結(jié) 論

采用含時(shí)量子波包方法對H/D+Li2反應(yīng)體系進(jìn)行了動力學(xué)計(jì)算,計(jì)算過程中考慮了科里奧利耦合效應(yīng),報(bào)道了體系的反應(yīng)概率、積分截面、微分截面并考慮了同位素效應(yīng).與之前報(bào)道的理論結(jié)果的比較表明,科里奧利耦合效應(yīng)在反應(yīng)過程中扮演重要的角色.由于本文考慮了總角動量J在體固定坐標(biāo)軸的所有投影,所以本文的結(jié)果更加準(zhǔn)確、合理.當(dāng)H原子被相對較重的D原子取代以后,體系的反應(yīng)概率、積分截面、微分截面均增加了.然而同位素效應(yīng)并未對反應(yīng)機(jī)理產(chǎn)生影響.從微分截面的結(jié)果可知,在整個(gè)反應(yīng)過程中插入反應(yīng)機(jī)理在反應(yīng)過程中占據(jù)主導(dǎo)地位.

[1]Chu T S,Zhang Y,Han K L 2006Int.Rev.Phys.Chem.25 201

[2]Pérez-Ríos J,Greene C H 2015J.Chem.Phys.143 041105

[3]Wang B B,Han Y C,Cong S L 2016J.Chem.Phys.145 204304

[4]Wang B B,Han Y C,Gao W,Cong S L 2017Phys.Chem.Chem.Phys.19 22926

[5]Wu C H,Ihle H R 1977J.Chem.Phys.66 4356

[6]Kim S K,Herschbach D R 1987Faraday Discuss.Chem.Soc.84 159

[7]Vezin B,Dugourd P,Rayane D,Labastie P,Broyer M 1993Chem.Phys.Lett.202 209

[8]Siegbahn P,Schaefer H F 1975J.Chem.Phys.62 3488.

[9]Yan G S,Xian H,Xie D Q 1997Sci.China Ser.B:Chem.40 342

[10]Maniero A M,Acioli P H,Silva G M,Gargano R 2010Chem.Phys.Lett.490 123

[11]Song Y Z,Li Y Q,Gao S B,Meng Q T 2014Eur.Phys.J.D68 3

[12]Yuan M L,Li W T,Chen M D 2017Int.J.Quant.Chem.e25380

[13]Kim S K,Jeoung S C,Tan A L C,Herschbach D R 1991J.Chem.Phys.95 3854

[14]Vila H V R,Leal L A,Martins J B L,Skouteris D,eSilva G M,Gargano R 2012J.Chem.Phys.136 34319

[15]Cunha W F,Leal L A,Cunha T F,Silva G M,Martins J B L,Gargano R 2014J.Mol.Model20 2315

[16]Gao S B,Zhang J,Song Y Z,Meng Q T 2015Eur.Phys.J.D69 111

[17]Gao S B,Zhang L,Song Y Z,Meng Q T 2016Chem.Phys.Lett.651 233

[18]Fu B N,Zhang D H 2012J.Chem.Phys.136 194301

[19]Shen P R,Han Y C,Li J L,Chen C J,Cong S L 2015Laser Phys.Lett.12 045302

[20]Pang Y H,Wang B B,Han Y C,Cong S L,Niu Y Y 2016Chin.J.Chem.Phys.29 297

[21]Gao W,Wang B B,Hu X J,Chai S,Han Y C,Greenwood J B 2017Phys.Rev.A96 013426

[22]Yuan J C,Cheng D H,Chen M D 2014RSC Adv.4 36189

[23]Duan Z X,Qiu M H,Yao C X 2014Acta Phys.Sin.63 063402(in Chinese)[段志欣,邱明輝,姚翠霞2014物理學(xué)報(bào)63 063402]

[24]Zhang J,Wei W,Gao S B,Meng Q T 2015Acta Phys.Sin.64 063101(in Chinese)[張靜，魏巍，高守寶，孟慶田2015物理學(xué)報(bào)64 063101]

[25]Yuan K J,Cheng Y,Liu X H,Harich S,Yang X M,Zhang D H 2006Phys.Rev.Lett.96 103202

[26]Hankel M,Smith S C,Allan R J,Gray S K,Balint-Kurti G G 2006J.Chem.Phys.125 164303

[27]Fu B,Zhang D H 2007J.Phys.Chem.A111 9516

[28]Zhang D H 2006J.Chem.Phys.125 133102

[29]Koslof fR 1988J.Phys.Chem.92 2087

[30]Light J C,Carrington T 2000Adv.Chem.Phys.114 263