CO分子四個電子態(tài)的振轉譜:兩種效應修正方法的比較?

徐慧穎 劉勇 李仲緣 楊玉軍 閆冰

1)(吉林大學原子與分子物理研究所,長春 130012)2)(吉林大學大數(shù)據(jù)和網(wǎng)絡管理中心,長春 130012)(2018年8月1日收到;2018年8月29日收到修改稿)

基于完全活性空間自洽場方法和多參考組態(tài)相互作用(multi-reference configuration interaction method,MRCI)方法,采用MRCI+Q/CBS(TQ5)+CV+SR(方法A)和aug-cc-pwCVnZ-DK(n=T,Q,5)(方法B)方案,分別計算了包含Davidson修正(+Q)、芯-價電子關聯(lián)(core-valence correlation correction,CV)效應以及標量相對論(scalar relativistic,SR)效應的CO分子的基態(tài)X1Σ+和激發(fā)態(tài)a3Π,a′3Σ+和A1Π的勢能曲線.在此基礎上,獲得了這些電子態(tài)的振-轉譜.通過與實驗結果比較發(fā)現(xiàn):方法A適合a′3Σ+和A1Π等較高激發(fā)態(tài)的振-轉譜的計算,方法B更適合基態(tài)X1Σ+和第一激發(fā)態(tài)a3Π的振-轉譜的精細計算.該研究可以為其他小分子高精度振-轉譜快速計算方案選擇提供參考.

1 引 言

分子內部精細的電子結構研究對探索物質的性質、結構和功能具有重要意義.對分子的微觀結構與光譜性質的計算一直是原子分子物理學科的重要內容之一.目前,多電子分子體系的準確計算,采用的主要方案是多參考組態(tài)相互作用方案.由于需要較多的計算資源,對于簡單的分子體系準確的全電子計算也存在較大困難.分子的外殼層性質對于諸多物理化學過程具有更重要的意義,為了實現(xiàn)對其準確快速計算,計算中通常將靠近原子核的電子作用采用一個有效勢來代替,計算后再將價電子與芯結構之間的關聯(lián)作用加到整體能量計算中.此外,為了得到更準確的能量,需要考慮相對論效應.為了準確計算該效應需要求解Dirac方程,其計算更為困難,相對論效應按性質可分為標量相對論(scalar relativistic,SR)效應(達爾文項和質量-速度項)和矢量相對論效應(自旋-軌道耦合項),本文僅考慮了SR效應.

對于SR效應的修正,人們在研究中采用了許多全電子近似的相對論方法.例如,Jong等[1]計算了Hartree-Fock(HF)水平的CF4,Si4和Br2CO的SR效應對平衡核間距Re及總解離能ΣDe的改變.對于芯-價電子關聯(lián)(core-valence correlation correction,CV)效應,人們針對第二行的原子Al-Ar的研究[2]使用了關聯(lián)效應的相關一致基組.對這兩個效應的修正可以通過如下兩種方案進行.

1)方法A

針對這兩個修正效應貢獻于體系的相對能量特點,采用不包含這兩個相互作用的基組較大基函數(shù),計算體系中不包含這兩個相互作用的能量,然后通過選擇新的較簡單的基函數(shù)快速計算其貢獻.2014年,Abbiche等[3]通過多參考組態(tài)相互作用(multi-reference configuration interaction,MRCI)+CV+SR方法對CP分子的X2Σ+和A2Π等7個較低電子態(tài)進行了計算.

2)方法B

在基函數(shù)選擇時選擇包含這兩個相互作用的基函數(shù),計算體系哈密頓中包含這兩個效應的總能量.2014年,Li等[4]通過aug-cc-pwCVnZ-DK基組計算了GeH+的8條電子態(tài)的勢能曲線(potential energy curves,PECs)和光譜常數(shù)等.

這兩種方案均可以給出較準確的計算結果,并能實現(xiàn)較快速的計算.但對于其適應的條件還需要深入的研究,即使對于較簡單分子體系的不同電子態(tài)的計算也是如此.

為了探索這兩個效應修正方案的差異,本文選擇CO分子加以研究.CO的研究在等離子體、星際物質、生物醫(yī)學以及環(huán)境檢測與分析具有重要意義.過去的數(shù)十年間,科研人員已對CO分子進行了大量的研究[5?21].2013年Lu等[22]采用aug-cc-pV5Z基組計算了CO分子較低的單重和三重態(tài)(包括X1Σ+和A1Π等6個電子態(tài))的PECs和基態(tài)X1Σ+在振動量子數(shù)ν=20之前的振-轉能級結構,其計算中并沒考慮到CV和SR效應修正.同年,Shi等[23]選用基組aug-cc-pV5Z計算了MRCI水平的CO分子8條較低電子態(tài)的PECs和振-轉譜,同時分別采用基組為cc-pCVQZ和cc-pV5Z考慮了CV效應和SR效應對于PECs的影響(aug-ccpV5Z+CV+DK).通過兩種方案,計算了兩種效應對CO分子能量的影響,進而分析了分子的振-轉譜.研究發(fā)現(xiàn)方法A適合較高激發(fā)態(tài)的振-轉譜的計算,方法B適合CO分子基態(tài)和第一激發(fā)態(tài)的振-轉譜的精細計算.

2 方法和計算細節(jié)

為了得到精確的CO分子的PECs和振-轉譜,使用了molpro2012[24]軟件包計算了分子的電子結構.為了實現(xiàn)效率和精確度的平衡,首先選擇HF方法獲得分子體系的波函數(shù);然后在此基礎上利用態(tài)平均的完全活性空間自洽場(complete active space seif-consistent field,CASSCF)方法[25,26]來描述電子之間的靜力學相關效應;最后采用MRCI方法[27,28]將電子間的動態(tài)相關效應也包含進來.

考慮CO分子的基態(tài)X1Σ+和激發(fā)態(tài)a3Π,a′3Σ+和A1Π,根據(jù)其對稱性,分子軌道從頭計算選擇在C∞v的子群C2v中進行.C∞v和C2v的不可約表示的對應關系分別為:Σ+-A1,Π-B1+B2,?-A1+A2,Σ?-A2.在計算這4個電子態(tài)的PECs時,選取了一系列單點能量計算,計算的核間距區(qū)域R和間隔?R分別為R=0.8—2.5 ?,?R=0.01 ?;R=2.5—5.2 ?,?R=0.1 ? (1 ?=0.1 nm). 每步的具體計算如下.首先,對每個單點進行HF自洽場計算,在HF計算中,采用了單組態(tài)Slater行列式描述CO的電子基態(tài),此處只考慮了自旋反平行電子的相關作用.然后,以HF方法產(chǎn)生的分子軌道作為初始軌道,進行CASSCF計算,分別獲得CO分子的基態(tài)X1Σ+和激發(fā)態(tài)a3Π,a′3Σ+和A1Π的波函數(shù).在CASSCF計算中選擇CO的10個分子軌道和10個n=2電子作為活性空間,包括4個a1,兩個b1和兩個b2對稱性的分子軌道,它們對應C原子的原子軌道2s2p和O原子的軌道2s2p.C的外層電子2s22p2和O的外層電子2s22p4被放置在活性空間內,剩下的4個電子被凍結而不進行相關能計算.最后,進行MRCI計算,但由于大小一致性誤差的存在,在MRCI中加入了Davidson修正(+Q)[29],估計了電子的四重激發(fā)修正,這樣就保證了在無窮遠分子的總能量近似等于單獨計算C和O原子的能量之和.

對于CV效應和SR效應的修正通過兩種方案進行.方法A首先采用基組aug-cc-pVnZ(n=T,Q,5)[30,31]通過HF,CASSCF和MRCI方法計算CO分子4個束縛態(tài)的PECs.為了討論內層電子相關效應的影響,在MRCI計算中,考慮了C和O原子的1s2內殼層的單、雙電子激發(fā)產(chǎn)生的CV效應.在不考慮CV效應的計算中,將上述的單、雙電子激發(fā)關聯(lián)軌道凍結.在上述所有計算中,為了考查CV效應的影響,基組均采用cc-pCVQZ.利用MRCI方法和非收縮aug-cc-pVQZ基組.通過計算三階Douglas-Kroll[32]和Hess[33]單電子積分獲得SR效應,即包含了質量速度項和Darwin項兩種相對論效應,然后將兩個修正效應貢獻于體系的相對能量加到不包含這兩種效應的能量上.在方法B中,選用了新的基組aug-cc-pwCVnZ-DK(n=T,Q,5)通過HF,CASSCF和MRCI方法計算CO分子4個束縛態(tài)的PECs.選擇的基組包含了這兩種效應修正的基函數(shù),計算體系哈密頓中包含這兩個效應的總能量.為了減小基組帶來的誤差,對兩種修正方案的基組(n=T,Q,5)外推到完全基組(complete basis set,CBS).

最后,根據(jù)方法A和方法B計算獲得的X1Σ+,a3Π,a′3Σ+和A1Π態(tài)的PECs,通過LEVEL[34]程序擬合求解CO分子的一維徑向Schr?dinger方程,得到各個電子態(tài)的振-轉譜.

3 結果和討論

利用方法A和方法B計算方案得到了CO分子基態(tài)X1Σ+和激發(fā)態(tài)a3Π,a′3Σ+和A1Π的PECs.它們對應的解離極限是C(3P)+O(3P).根據(jù)方法A計算的PECs(圖1),利用了LEVEL程序對4個電子態(tài)進行擬合得到光譜參數(shù),通過比較發(fā)現(xiàn)其值與實驗值十分接近,特別是對于轉動常數(shù)Be和平衡核間距Re的計算,與實驗值的相對誤差均小于0.0001 ?和0.001 cm?1.

圖1 方法A計算的基態(tài)X1Σ+和激發(fā)態(tài)a3Π,A1Π和a′3Σ+ 的PECs(1 hartreee=27.2114 eV)Fig.1.PECs of ground state X1Σ+and excited states a3Π,A1Π and a′3Σ+calculated by method A.

基于方法A和方法B計算方案得到的基態(tài)X1Σ+和激發(fā)態(tài)a3Π,a′3Σ+和A1Π 的PECs,可以計算其振動能級Gν. 表1根據(jù)兩種計算方案列出了基態(tài)X1Σ+和第一激發(fā)態(tài)a3Π的前21個Gν,從表1可以看出,根據(jù)方法A,在ν=0—10時(轉動量子數(shù)J=0),振動能級Gν的計算結果范圍為1082.75—21354.78 cm?1,相對偏差為0.97—23.63 cm?1,相對偏差的變化范圍為0.09%—0.11%. 但在ν>10時,方法A的誤差較大,振動能級Gν的誤差的變化范圍為26.56—157.15 cm?1,并且由方法A計算的振動能級Gν的均方根誤差(root mean square error,RMSE)為0.17%.根據(jù)方法B,在ν=0—10時(J=0),振動能級Gν的計算結果范圍為1082.11—21344.55 cm?1,相對偏差為0.33—13.41 cm?1,相對誤差范圍為0.03%—0.06%.在較低振動態(tài)時方法B比方法A的計算結果更接近實驗測量值[35]. 在ν>10時,方法B的計算結果也更精確,Gν計算結果范圍為23228.85—39143.40 cm?1,對應的相對偏差為15.67—144.53 cm?1. 通過方法B計算的Gν的RMSE為0.13%.

對于第一激發(fā)態(tài)a3Π,實驗[36]上只給出了a3Π態(tài)的前8個振動能級(ν=0—7)的測量值. 根據(jù)方法A計算方案的振動能級結構的范圍為870.5872—12302.4342 cm?1,相對偏差為2.43—56.43 cm?1,RMSE為0.37%.方法B計算方案的振動能級的范圍為870.06—12296.18 cm?1,相對偏差為1.9—50.1 cm?1,RMSE為0.32%.通過數(shù)據(jù)對比發(fā)現(xiàn),對于基態(tài)X1Σ+和第一激發(fā)態(tài)a3Π,方法B比方法A計算的振-轉譜更符合實驗測量值.

通過方法A和方法B計算的激發(fā)態(tài)a′3Σ+和A1Π的PECs擬合出前21個振動能級列于表2.表2同時列出了相應的實驗測量值[36].對于a′3Σ+態(tài),方法A方案的計算結果與實驗值[36]符合得較好,計算結果范圍為615.4681—21161.94 cm?1,相對偏差為2.89—129.16 cm?1,RMSE為0.54%.方法B計算的結果范圍為615.63—21172.16 cm?1,相對偏差為3.05—139.38 cm?1,RMSE為0.57%.同 樣, 對 于A1Π態(tài), 方 法A計 算 結 果 范圍 為755.3099—23926.31 cm?1, 相 對 偏 差 為1.82—164.31 cm?1,RMSE為0.49%;方法B的計算結果范圍為755.2961—23938.30 cm?1,相對偏差為1.81—182.27 cm?1,RMSE為0.53%.通過分析上述結果可知,對于激發(fā)態(tài)a′3Σ+和A1Π,方法A方案更接近于實驗的測量.

為了更清楚地對比兩種計算方案得到的振-轉譜的計算精度,計算了模擬結果與實驗結果的相對誤差,如圖2所示.從圖2可以看出,兩種方法計算的振轉能級與實驗的相對整體都小于1%.但是對于不同電子態(tài)還存在一定差距,對于基態(tài)X1Σ+和第一激發(fā)態(tài)a3Π,方法A的計算結果相對誤差大于方法B,而對于激發(fā)態(tài)A1Π和a′3Σ+,方法A的計算結果優(yōu)于方法B.這一計算進而可以推廣到更高激發(fā)態(tài)的計算.通過對上述結果分析可知,方法B在波函數(shù)中包含修正效應,具有更好的普適性,通常對于小分子計算可以采用該方案.而對于較重分子,可以應用方法A,可以在相對小的計算需求條件下得到較準確的結果.本文采用不同的方案計算了SR效應與CV修正,并且后者對精密計算是不可或缺的.對于最低的兩個電子態(tài),考慮了兩種效應對計算高斯基組的依賴性(方法B),發(fā)現(xiàn)振動-轉動的計算精度提高了,可見,對于這兩個電子態(tài)對電子關聯(lián)計算的要求較高;對于更高的電子態(tài),電子云分布相對松散,反而簡單地采用單一高斯基組獲得的電子關聯(lián)修正即能達到相應的計算準確性;當然,由于振-轉譜的計算實質上只有相對能量起作用,這里也包含了不同電子態(tài)電子關聯(lián)效應的抵消效應.

表1 通過方法A和方法B計算得到的CO分子X1Σ+和a3Π態(tài)的振動能級Gν(ν=0—20,J=0)Table 1.Vibrational energy levels Gν of X1Σ+and a3Π state of CO molecule calculated by method A and method B(ν =0–20,J=0).

圖2 通過方法A和方法B計算基態(tài)X1Σ+和激發(fā)態(tài)a3Π,A1Π和a′3Σ+振動能級Gν的相對誤差隨ν的變化Fig.2.Relative error of the vibrational energy level Gν of ground state X1Σ+and excited states a3Π,A1Π and a′3Σ+calculated by method A and method B varies with ν.

4 結 論

本文采用了兩種方案計算包含SR和CV效應的CO基態(tài)X1Σ+和激發(fā)態(tài)a3Π,a′3Σ+和A1Π的PECs和振-轉譜,通過和實驗結果的比較發(fā)現(xiàn):

1)計算基態(tài)X1Σ+和第一激發(fā)態(tài)a3Π的振-轉譜時,基組為aug-cc-pwCVnZ-DK(n=T,Q,5)方法(方法B)的計算結果更符合實驗測量值,方法A的計算方案擬合得到的振-轉譜誤差較大;

2)計算較高的激發(fā)態(tài)a′3Σ+和A1Π或更高的激發(fā)態(tài)時,方法A的計算結果更加準確,相對誤差較小.

得到的CO分子基態(tài)X1Σ+和激發(fā)態(tài)a3Π,a′3Σ+和A1Π的較完整的振動能級信息,對今后的實驗研究具有參考價值,同時為獲得高精度分子的振-轉譜計算方案的選擇提供了依據(jù).

感謝吉林大學超算中心的計算支持.