李文濤于文濤3)姚明海

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

2)(中國科學院大連化學物理研究所,分子反應動力學國家重點實驗室,大連 116023)

3)(東風縣第三中學,遼源 136300)

1 引 言

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

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

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

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

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

2 動力學計算方法

量子含時波包量子方法在動力學計算中已廣泛使用[23,24],本文只做一些簡單的介紹,詳細的信息可以參考文獻[25—28].對于給定的總角動量J,H+Li2反應體系的哈密頓量在體坐標系下可以寫為

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

其中Ca,Cb是吸收勢的強度;?t是時間步長;ra,rb,rend是吸收勢對應的位置.

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

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

其中是散射矩陣,d?是方位角,?是散射角,是魏格納轉動矩陣.

3 結果與討論

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

3.1 反應概率

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

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

圖2 在0—0.4 eV的碰撞能范圍內,H+Li2反應總角動量J=0,10,20,30,40,50的反應概率以及文獻[14,17]的結果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的反應概率見圖2.如圖2所示,在J=0的反應概率上沒有閾值出現(xiàn)并且在研究的碰撞能范圍內有一些共振峰,這些共振峰的產生是因為在反應路徑上有一個能產生大量束縛態(tài)和準束縛態(tài)的深勢阱;伴隨著碰撞能的增加,共振峰的強度變弱,共振峰之間的間距變寬.這表明隨著碰撞能的增加,碰撞絡合物的壽命開始變短.當總角動量J變大以后,反應閾值開始出現(xiàn)并且隨著J的增大而增大.反應閾值的出現(xiàn)是因為離心勢的變大,離心勢隨著J的增大而變大.

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

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

從圖3中不難看出,當k=15時,得到的結果是不收斂的.從圖3(b)可以看出,反應概率在高碰撞能范圍的巨大差距主要是因為文獻[17]沒有把所有的k-block值考慮在計算中.

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

圖3 當總角動量J=30時,不同k-block值的反應概率(a)和文獻[17]結果的比較(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反應在總角動量J=0,60和70時的反應概率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反應的積分截面在0—0.2 eV的范圍內隨著碰撞能的增加而逐漸降低,這是一個典型的放熱反應特征.由圖5可知,同位素取代增大了反應的積分截面,尤其是在低碰撞能的時候.這是可以預料的,因為同位素取代也導致了反應概率的增加.

圖5 H/D+Li2反應在碰撞能0—0.2 eV范圍內的積分截面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反應在0.1和0.15 eV時的微分截面Fig.6. The differential cross sections of H/D+Li2 reaction at 0.1 and 0.15 eV,respectively.

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

4 結 論

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

[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物理學報63 063402]

[24]Zhang J,Wei W,Gao S B,Meng Q T 2015Acta Phys.Sin.64 063101(in Chinese)[張靜，魏巍，高守寶，孟慶田2015物理學報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