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

1)(吉林大學(xué)原子與分子物理研究所,長春 130012)2)(吉林大學(xué)大數(shù)據(jù)和網(wǎng)絡(luò)管理中心,長春 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)、芯-價電子關(guān)聯(lián)(core-valence correlation correction,CV)效應(yīng)以及標(biāo)量相對論(scalar relativistic,SR)效應(yīng)的CO分子的基態(tài)X1Σ+和激發(fā)態(tài)a3Π,a′3Σ+和A1Π的勢能曲線.在此基礎(chǔ)上,獲得了這些電子態(tài)的振-轉(zhuǎn)譜.通過與實(shí)驗(yàn)結(jié)果比較發(fā)現(xiàn):方法A適合a′3Σ+和A1Π等較高激發(fā)態(tài)的振-轉(zhuǎn)譜的計算,方法B更適合基態(tài)X1Σ+和第一激發(fā)態(tài)a3Π的振-轉(zhuǎn)譜的精細(xì)計算.該研究可以為其他小分子高精度振-轉(zhuǎn)譜快速計算方案選擇提供參考.

1 引 言

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

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

1)方法A

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

2)方法B

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

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

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

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

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

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

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

3 結(jié)果和討論

利用方法A和方法B計算方案得到了CO分子基態(tài)X1Σ+和激發(fā)態(tài)a3Π,a′3Σ+和A1Π的PECs.它們對應(yīng)的解離極限是C(3P)+O(3P).根據(jù)方法A計算的PECs(圖1),利用了LEVEL程序?qū)?個電子態(tài)進(jìn)行擬合得到光譜參數(shù),通過比較發(fā)現(xiàn)其值與實(shí)驗(yàn)值十分接近,特別是對于轉(zhuǎn)動常數(shù)Be和平衡核間距Re的計算,與實(shí)驗(yàn)值的相對誤差均小于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時(轉(zhuǎn)動量子數(shù)J=0),振動能級Gν的計算結(jié)果范圍為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ν的計算結(jié)果范圍為1082.11—21344.55 cm?1,相對偏差為0.33—13.41 cm?1,相對誤差范圍為0.03%—0.06%.在較低振動態(tài)時方法B比方法A的計算結(jié)果更接近實(shí)驗(yàn)測量值[35]. 在ν>10時,方法B的計算結(jié)果也更精確,Gν計算結(jié)果范圍為23228.85—39143.40 cm?1,對應(yīng)的相對偏差為15.67—144.53 cm?1. 通過方法B計算的Gν的RMSE為0.13%.

對于第一激發(fā)態(tài)a3Π,實(shí)驗(yàn)[36]上只給出了a3Π態(tài)的前8個振動能級(ν=0—7)的測量值. 根據(jù)方法A計算方案的振動能級結(jié)構(gòu)的范圍為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計算的振-轉(zhuǎn)譜更符合實(shí)驗(yàn)測量值.

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

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

表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 結(jié) 論

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

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

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

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

感謝吉林大學(xué)超算中心的計算支持.