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

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

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

1 引 言

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

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

1)方法A

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

2)方法B

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

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

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

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

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

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

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

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

3 結(jié)果和討論

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

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

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

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

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

表1 通過方法A和方法B計(jì)算得到的CO分子X1Σ+和a3Π態(tài)的振動(dòng)能級(jí)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計(jì)算基態(tài)X1Σ+和激發(fā)態(tài)a3Π,A1Π和a′3Σ+振動(dòng)能級(jí)Gν的相對(duì)誤差隨ν的變化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é) 論

本文采用了兩種方案計(jì)算包含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)計(jì)算基態(tài)X1Σ+和第一激發(fā)態(tài)a3Π的振-轉(zhuǎn)譜時(shí),基組為aug-cc-pwCVnZ-DK(n=T,Q,5)方法(方法B)的計(jì)算結(jié)果更符合實(shí)驗(yàn)測(cè)量值,方法A的計(jì)算方案擬合得到的振-轉(zhuǎn)譜誤差較大;

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

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

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