F+CHD3→HF+CD3 反應C—H 伸縮振動激發的量子動力學研究*

周勇

(安徽師范大學物理與電子信息學院,蕪湖 241000)

對多原子化學反應進行精確的理論模擬研究,是當前量子動力學領域的重要挑戰之一.本文運用七維的量子含時波包方法,對反應物振動基態和C—H 伸縮振動第一激發態的F+CHD3 反應進行動力學研究.分析了不同振動態條件下的反應概率,發現當碰撞能低于0.06 eV 時反應概率曲線均表現出許多劇烈的振蕩峰,支持實驗上所推測的動力學共振現象.當碰撞能介于0.06—0.3 eV 之間時,振動激發態反應中HF 產物通道的反應概率低于振動基態反應,與實驗觀測結果一致.模擬發現: 相比于振動基態的情況,低能碰撞時C—H 伸縮振動激發態反應的非含時波函數在過渡態區域更傾向于靠近D 原子一側,并將這一現象歸結為振動激發態勢能面在大碰撞角區域更顯著的能量優勢,解釋了伸縮振動激發對HF 產物通道的抑制作用.本研究為實驗結果的解釋提供了重要的理論支持,有助于更深入地理解多原子反應中振動模式激發對動力學過程的影響.

1 引言

過去幾十年中,氣相化學反應的精確理論研究取得了很大的進展.尤其是隨著含時波包方法的發展,現在已經可以在沒有動力學近似的情況下,對四原子反應進行初始態選擇的研究[1,2].理論工作者報道的關于H+H2O → H2+OH 反應[3]和HD+OH → H2O+D 反應[4]微分截面的研究結果,標志著對四原子反應問題的理論解決已經接近尾聲.目前,在量子動力學領域的一個主要挑戰是發展用于研究包含多于4 個原子化學反應的精確方法.

作為重要的多原子反應體系,F+CH4反應及其同位素取代反應近年來吸引了來自實驗和理論方面的廣泛關注[5-11].Liu等[12]利用“離子速度成像”(velocity map ion imaging)技術獲得了F+CH4反應的配對相關的(pair-correlated)態分辨激發函數;通過研究反應能量閾值附近出現的“臺階”現象,他們推測該反應中可能存在反應共振現象[13].Liu 及其合作者們[14]還使用類似的實驗技術研究了F+CHD3反應HF 通道產物對之間的模式相關性(mode correlation),并依據在反應積分截面的低能區域所發現的明顯臺階特征,推測該反應通道也存在反應共振[15].理論研究方面,Troya等[16]將CH3看作偽原子構造了三維的基態解析勢能面,并基于此勢能面的準經典軌線研究復現了實驗上測得的HF 內態分布,證實了實驗中發現的立體反應動力學現象.Espinosa-Garcia等[17]構造了F(2P)+CH4及其同位素反應的勢能面,并將其用于反應速率常數計算和動力學研究.Czakó等[18]構建了F+CH4體系的全維全域勢能面,在此基礎上進行的準經典軌線計算得到了和實驗結果吻合的HF 產物振轉分布.Nyman 和Espinosa-Garcia[19]基于三維量子散射模型的計算,為F+CH4反應中的共振現象提供了理論證據.Wang 和Czakó[20]在四維量子動力學計算得到的總反應概率和積分截面中發現了反應共振的相關特征.Westermann等[21]基于全維的多組態含時Hartree(MCTDH)方法計算了F+CH4反應的過渡態光譜,發現了反應入口通道范德瓦耳斯(vdW)勢阱里的反應共振.von Horsten和Clary[22]對F+CHD3反應進行的兩維量子動力學研究,計算了兩個產物通道末態分辨的反應概率和截面,并分析了低能碰撞區域處的共振態.Qi等[23]和Zhao等[24]分別開展了F+CHD3反應的八維和七維波包動力學計算,均報道了低碰撞能條件下總反應概率曲線上的快速振蕩結構,并將之歸結為反應前的vdW 勢阱.最近,Zhang 課題組構建了一個新的F+CH4 反應勢能面[25]并報道了對F+CHD3→HF+CD3振轉基態反應的六維態-態量子動力學研究結果,在總的和態分辨的反應概率曲線中觀察到顯著的振蕩結構并將其歸結為由HF 化學鍵軟化所導致的HF(v'=3)-CD3振動絕熱勢阱[26].

Zhang等[27]開展了對F+CHD3反應中的C—H鍵伸縮振動激發的實驗研究.對于這一早期勢壘(early barrier)反應,根據經驗推測反應物的振動應該不會和反應坐標有效耦合,因此對反應產出的影響很小;然而,實驗中發現C—H 鍵的伸縮振動激發反而抑制了C—H 鍵的斷裂,即傾向于DF+CHD2產物通道.Czakó等[28]通過對F+CHD3(vC—H=0,1)反應準經典軌線的統計分析發現: C—H鍵的伸縮振動激發會導致F 原子在反應鞍點區域傾向于靠近C—D 鍵.雖然這一發現可以在一定程度上解釋C—H 鍵伸縮振動激發對C—H 鍵斷裂的抑制作用和實驗上所觀察到的鍵選擇性現象,但是對于F 原子接近CHD3分子時會偏向于D 原子一側的原因,Czakó等[28]并未給出更進一步的解釋.此外,受準經典軌線方法本身的限制,其僅能對F+CHD3反應進行比較粗糙的描述,例如其很難精確捕捉F+CHD3振動激發態反應中可能存在的量子反應共振現象.

量子理論可以為化學反應動力學過程提供最精確的描述.然而,與H+CH4,H+CHD3,H+CD4等被廣泛研究的六原子反應體系[29-33]相比,對F+CHD3反應的量子波包動力學研究更加具有挑戰性.一方面,入射原子由質量較小的H 原子變為質量較大的F 原子,對于可能發生共振現象的低能碰撞,波包需要傳播的時間要長得多;其次,由于F 原子和CHD3分子之間存在長程相互作用,為獲得更準確的結果,需要將初始波包的位置放在很遠的位置(大于15 Bohr),這將進一步延長波包的傳播時間.本文采用七維的含時波包方法,研究了振動基態和C—H 伸縮振動第一激發態F+CHD3反應當總角動量Jtotal=0 的情況,發現低碰撞能情況下反應概率曲線上均出現了暗示反應共振的強烈振蕩現象,重現了實驗上觀測到的C—H 鍵伸縮振動激發對HF 產物通道的抑制現象,并基于振動態平均的勢能面對反應中的鍵選擇性進行了解釋.

2 理論方法

2.1 X+YCZ3 型反應的量子含時波包方法

圖1 定義了用于研究X+YCZ3→XY+CZ3型反應的八維Jacobi 坐標系統.R為從YCZ3分子的質心到X 原子的向量,r為從CZ3基團的質心到Y 原子的向量.s為C—Z 鍵的鍵長,在7D 計算中該鍵長被固定為反應物C—Z 鍵的平衡鍵長.χ為C—Z 鍵和CZ3基團的C3V對稱軸的夾角.為描述角度坐標和系統的轉動,這里引入了4 個框架,分別為空間固定(space-fixed)框架、XYCZ3固定框架、YCZ3固定框架和CZ3固定框架.YCZ3框架的z軸沿著向量r,向量s總是在該框架的xz平面上.CZ3固定框架的z軸沿著其對稱軸,同時第一個Z 原子始終在此框架的x-z平面中.定義R和r之間的彎曲角為θ1,r和s之間的彎曲角為θ2,YCZ3圍繞向量r的扭轉角為φ1,CZ3基團圍繞其對稱軸的扭轉角為φ2.

圖1 用于研究X+YCZ3 反應的8 維Jacobi 坐標系統Fig.1.8D Jacobi coordinate system for X+YCZ3 reactions.

基于八維模型的X+YCZ3反應的哈密頓量用上述Jacobi 坐標表示為[29,34]

含時波函數用宇稱匹配(parity-adapted)的轉動基組展開為

對應于某個態(Jtot,M,ε) 的系統初始波函數寫成局域波包G0(R) 和對應于(n0,J0,K0,p0) 狀態下YCZ3的本征函數這兩部分的直積,其中n0,J0,K0,p0分別表示初始振動態、總角動量及其在XYCZ3固定框架z軸上的投影和YCZ3分子的宇稱.

本工作使用“分裂算符”方法進行波函數的傳播:

其中參考哈密頓量H0定義為

而參考勢能U定義為

對應于某個特定初始態的在整個能量范圍內的總反應概率,可以首先通過對含時波函數進行傅里葉變換得到非含時波函數ψiE,然后在分割面r=rs上對非含時波函數及其一階導數積分得到:

2.2 計算參數

本文使用2.1 節中所介紹的方法對F+CHD3反應進行研究.為了減小計算中所用基組的大小,對R和r坐標進行“L 形”的波函數展開.R坐標的取值范圍為3—21 Bohr,共使用302 個正弦基組,其中在相互作用區域共使用98 個格點.對于r方向,反應漸進區和相互作用區分別使用7 個和50 個基函數.C—D 鍵的長度固定為反應物CHD3分子的平衡鍵長(2.06 Bohr);同時,使用了9 個基函數來描述CD3基團的傘形振動.約束轉動基函數大小的參數分別為Jmax=114,lmax=90,jmax=24和kmax=24.在考慮了CD3基團的C3V對稱性的情況下,展開波函數時共使用了284927 個轉動基函數,在對轉動基函數進行數值積分的過程中共使用了2394875 個角度格點.高斯波包的初始中心位于R0=18.5 Bohr,波包寬度δ=0.2 Bohr,波包中心能量為E0=0.2 eV.本工作中采用的是Czakó等[18]構造的全維全域勢能面.

3 結果與討論

圖2 顯示了總角動量Jtotal=0 情況下,反應物振動基態和C—H 鍵伸縮振動第一激發態的總反應概率隨碰撞能的變化.由圖2 可知,無論反應物是處于振動基態還是振動激發態,反應的能量閾值均非常小.Czakó等[18]在UCCSD(T)/aVTZ 水平上進行的電子結構計算表明,該反應的反應能壘很低,僅為0.5 kcal/mol 左右.考慮到該H 原子抽取反應通道可能存在的隧穿(tunneling)效應,這里看到的低反應能量閾值情況與Czakó等[18]的計算結果吻合.在當前所考慮的碰撞能范圍內,總反應概率的最大值可以達到15%左右,遠遠高于H+CH4及其同位素反應的相關結果(約0.3%),這表明了F 原子具有更強的反應活性.

圖2 總角動量Jtotal=0 時,F+CHD3 (vC—Hstr=0,1) →HF+CD3 反應的總反應概率Fig.2.Total reaction probabilities of the F+CHD3 (vC-Hstr=0,1) → HF+CD3 reaction with total angular momentum Jtotal=0.

對于反應物處于振動基態的情況,當碰撞能Ec＜ 0.06 eV 時,反應概率隨碰撞能劇烈變化并形成了許多尖銳的振蕩峰.這一現象與Liu等[26]基于神經網絡勢能面所進行的六維量子波包結果吻合得很好,說明反應過程中確實存在共振現象.從Ec=0.06 eV 附近開始,隨著碰撞能量的增大,反應概率隨碰撞能的振蕩現象減弱,同時反應概率緩慢增大.Zhou等[15]在實驗中發現,當碰撞能小于1.2 kcal/mol(0.052 eV)時,反應激發函數曲線中會出現明顯的“臺階”特征,并據此推測反應中可能存在共振現象.這一點現在也被本工作的理論模擬結果所支持.此外,實驗結果還顯示,當碰撞能升高至0.052 eV 附近時激發函數會明顯增大,本工作模擬得到的總反應概率中也發現了類似的現象.此外,von Horsten 和Clary[22]的二維量子動力學研究也發現: 在低能(Ec＜ 6 kJ/mol,合0.062 eV)碰撞區域,共振結構起主導作用;在5—6 kJ/mol的碰撞能附近,由于HF (v'=3)產物通道的打開,總反應概率也會突然增大并且隨著碰撞能的增大而上升.這些結果也都與本工作的模擬結果一致.

對于反應物C—H 鍵伸縮振動第一激發態的情況,反應概率隨碰撞能的變化曲線與基態反應情況類似: 當碰撞能Ec＜ 0.06 eV 時,可以看到許多明顯的振蕩峰;隨著碰撞能的增大,反應概率逐漸增大.不同的是,當碰撞能Ec在0.06—0.3 eV 范圍內時,振動激發態的反應概率要小于基態反應;這一現象與Zhang等[27]的實驗結果在定性上吻合得很好.定量而言,理論得到的基態反應概率最大僅為激發態反應概率的1.5 倍左右,這與實驗結果[27]測得的反應激發函數比值(5—10 倍)差別仍比較大.為了實現理論與實驗之間的一對一比較,后續有必要進一步地計算振動激發態的反應積分截面.

Zhou等[15]進行的實驗研究表明,對于振動基態反應,當碰撞能小于1.2 kcal/mol(合0.052 eV)時,生成的HF 產物基本上完全處于v′=2 振動態.為從理論上證明這一點,圖3 給出了低能碰撞情況下產物通道的非含時波函數.這里采用了與三原子態-態反應研究中“反應物→產物”坐標變換類似的方法,以便更清楚地觀察產物的振動態特征.圖3中豎直方向為產物反應坐標,其中下方為反應過渡態區域,上方為產物漸進區;水平方向對應于HF產物坐標.圖3(a)顯示了振動基態反應、碰撞能Ec=0.018 eV 的情況.由圖3(a)可知,在產物漸進區域波函數僅能分辨出2 個節點,此時波函數對應于HF 產物伸縮振動的v′=2 態,這一點與實驗觀察到的結果[15]吻合.在反應過渡態區域波函數沿著水平方向出現了3 個明顯的節點,說明在該區域HF 伸縮振動對應于v′=3 激發態;Liu等[26]基于神經網絡勢能面所進行的六維量子波包模擬也發現了類似的現象,并將過渡態區HF(v′=3)態到產物漸進區HF(v′=2)態的轉變歸結為Feshbach共振.

圖3 產物通道的非含時波函數(a)振動基態,碰撞能Ec=0.018 eV;(b) C—H 鍵伸縮振動第一激發態,碰撞能Ec=0.025 eV;圖中豎直方向和水平方向分別對應散射坐標和HF 產物坐標Fig.3.Time-independent wave function in product channel:(a) Ground vibrational state with collision energy Ec=0.018 eV;(b) the first C—H stretching vibrational excited state with collision energy Ec=0.025 eV.The vertical and horizontal directions correspond to the scattering coordinate and the HF product coordinate,respectively.

量子化學計算結果[15]表明,對于振動基態反應,當反應碰撞能Ec＞1.14 kcal/mol(0.05 eV),HF(v′=3)產物通道將打開.我們計算得到這一振動激發態的激發能為0.369 eV,這一數值與Czakó和Bowman[35]使用變分組態相關(variational configuration interaction)方法計算得到的結果2985 cm-1(0.370 eV) 吻合得非常好.因此,從可資用能量的角度推測,對于C—H 鍵伸縮振動第一激發態,即使對于很低的碰撞能,也可能會生成HF(v′=3)產物.圖3(b)顯示了CHD3分子C—H鍵伸縮振動第一激發態、碰撞能Ec=0.025 eV 時的產物通道非含時波函數.由圖3(b)可知,在反應過渡態區域附近,與圖3(a)中振動基態的情況類似,波函數沿著水平方向也表現為3 個節點,對應HF(v′=3)態.然而,與振動基態情況不同的是,隨著產物反應坐標從過渡態區域變化到產物漸進區域,振動激發態反應的產物波函數仍然保留著3 個節點的特征,對應于HF(v′=3)產物.此點與上述基于能量分析的推測結果一致.

為了研究伸縮振動激發抑制C—H 鍵斷裂的原因,選擇兩個有代表性的碰撞能,分別為0.06 eV和0.2 eV.由圖2 可知,當碰撞能Ec=0.06 eV時,反應物伸縮振動第一激發態所對應的反應概率明顯低于振動基態情況(相對差別～35%);而當Ec=0.2 eV 時,二者反應概率則非常接近(相對差別約4%).

采用與之前研究H+CD4反應[29]時類似的做法,圖4 顯示了碰撞能Ec=0.06 eV 時,反應物區域的非含時波函數隨反應坐標R、入射F 原子傾斜角θ(即圖1 中定義的θ1)變化的情況.對于反應物處于振動基態的情況(圖4(a)),反應坐標R=9 Bohr 處的波函數在大傾斜角和小傾斜角區域的分布之間差別不大.隨著趨近于過渡態區域,可發現波函數在大傾斜角區域的分布逐漸增多,即H 原子逐漸偏離C—H 鍵并向C—D 鍵靠攏.類似地,Czakó和Bowman[36]對準經典軌線的統計分析也表明,振動基態情況下入射F 原子較靠近C—D鍵和較靠近C—H 鍵這兩種情況的比例大約為3∶1,與CHD3分子中對應原子數目的比值相近.圖4(b)顯示了反應物處于C—H 鍵伸縮振動第一激發態的情況.與振動基態情況類似,隨著波函數向過渡態區域逐漸靠近,其更加傾向于分布在大傾斜角一端.不過,與振動基態情況不同的是,當C—H 伸縮振動模式被激發后,過渡態區域附近的波函數在小傾斜角區域的分布要明顯變少;絕大部分波函數分布在大傾斜角區域,并在傾斜角θ=140°—160°附近積聚成很高的峰.這一區別將直接導致C—H 鍵伸縮振動激發所對應的HF 產物生成概率大大減小.之前Czakó和Bowman[36]所進行的準經典軌線分析表明: C—H 伸縮振動激發情況下入射F 原子較靠近C—D 鍵和較靠近C—H鍵這兩種情況的比例遠遠高于D/H 原子數目之比(最高可達11 左右).這一點與圖4(b)中所展示的結果在定性上也是一致的.

圖4 碰撞能Ec=0.06 eV 時,反應物區的非含時波函數隨反應坐標R 和入射F 原子傾斜角θ 的變化(a)振動基態;(b) C—H鍵伸縮振動第一激發態Fig.4.Time-independent wave functions in the reactant region as a function of reaction coordinate R and the inclination angle θ of incoming F atom at a collision energy of Ec=0.06 eV: (a) Ground vibrational state;(b) the first C—H stretching vibrational state.

圖5 展示了碰撞能Ec=0.2 eV 時,對應于反應物振動基態和C—H 鍵伸縮振動第一激發態的非含時波函數.對于振動基態的情況(圖5(a)),與低碰撞能結果(圖4(a))類似,反應過渡態處的波函數同樣具有較大的比例分布在小碰撞角區域.類似地,對于振動激發的情況(圖5(b)),即便是當非常靠近過渡態時,在小傾斜角區域也仍然分布有較多的波函數,這一點與圖4(b)中所示的低碰撞能情況有明顯不同.通過比較圖4(b)和圖5(b)可知,的確可以將C—H 鍵伸縮振動激發對HF 產物通道的抑制作用歸結為入射F 原子相對于C—H鍵的偏離.

圖5 碰撞能Ec=0.2 eV 時,反應物區的非含時波函數隨反應坐標R 和入射F 原子傾斜角θ 的變化(a)振動基態;(b) C—H鍵伸縮振動第一激發態Fig.5.Time-independent wave functions in the reactant region as a function of reaction coordinate R and the inclination angle θ of incoming F atom at a collision energy of Ec=0.2 eV: (a) Ground vibrational state;(b) the first C—H stretching vibrational state.

為了解釋C—H 鍵伸縮振動模式被激發時入射F 原子偏離C—H 鍵的原因,接下來進一步分析不同振動模式下反應的有效勢能(effective potential).這里借鑒了前人研究H+H2O 反應時所采用的分析方法[37],將勢能對C—H 鍵的鍵長r進行平均,即:

其中,?v1(r) 是C—H 鍵的振動本征函數,Vref(r)是r方向的參考勢.如圖6 所示,當反應坐標R值較大,即遠離反應過渡態區域時,反應物振動基態和C—H 鍵伸縮振動第一激發態的平均勢能之間差別不大.隨著逐漸靠近強相互作用區域,兩個不同振動態對應的勢能面均開始向大傾斜角傾斜.對于振動基態勢能面,小傾斜角和大傾斜角區域之間的差別相對不那么顯著;而當C—H 鍵伸縮振動模式被激發之后,與基態勢能面相比,小傾斜角區域的勢能有小幅度的升高,同時大傾斜角區域的勢能則進一步的降低.從能量優勢的角度而言,振動激發情況下,波包將更傾向于以類似于“坐滑梯”的方式達到大傾斜角區域,從而導致了入射F 原子偏離C—H 鍵而向C—D 鍵靠近.這一點同時也印證了Zhang等[27]之前為解釋實驗結果所推測的、F 原子與伸縮振動激發的C—H 鍵之間相互作用的各向異性改變.

圖6 振動態平均的反應有效勢能隨反應坐標R 和入射F 原子傾斜角θ 的變化(a) 振動基態;(b) C—H 鍵伸縮振動第一激發態Fig.6.Averaged effective potential energy over vibrational states as a function of the reaction coordinate R and the inclination angle θ of incoming F atom: (a) Ground vibrational state;(b) the first C—H stretching vibrational excited state.

4 結論

本文采用七維量子含時波包方法研究了反應物振動基態和C—H 伸縮振動第一激發態的F+CHD3反應動力學.在總角動量為0 情況下,當碰撞能低于0.06 eV 時振動基態的反應概率曲線上出現許多尖銳的振蕩峰,與前人基于不同減維模型或勢能面所得到的理論模擬結果吻合得很好,支持實驗上所推測的動力學共振現象.當碰撞能低于0.06 eV 時,振動激發態的反應概率曲線也表現出類似于振動基態反應的劇烈振蕩峰,暗示此時反應過程中也可能存在動力學共振現象.計算發現,對介于0.06—0.3 eV 之間的碰撞能,振動激發態反應HF 產物通道的反應概率低于振動基態反應,與實驗觀測結果一致.與前人的準經典軌線模擬結果類似,低碰撞能情況下振動激發態反應的非含時波函數在過渡態區域更傾向于分布在大碰撞角附近(D 原子側);而對于高碰撞能情況,不同振動態下過渡態區域非含時波函數的碰撞角分布之間則差別不大.對反應有效勢能的分析表明,相對于振動基態的情況,C—H 伸縮振動激發態勢能面在大碰撞角區域的能量優勢更加明顯,由此導致伸縮振動激發對HF 產物通道的抑制作用.