徐又捷 郭迎春? 王兵兵

1)(華東師范大學(xué)物理與電子科學(xué)學(xué)院,上海 200241)

2)(中國科學(xué)院物理研究所光物理實驗室,北京凝聚態(tài)物理國家研究中心,北京 100190)

3)(中國科學(xué)院大學(xué),北京 100049)

針對較大分子振動頻率的量化計算,提出了一個節(jié)省計算成本的方法.含N 個原子的分子的振動頻率的計算通常需要計算3N 維勢能超曲面及其二階導(dǎo)數(shù)構(gòu)成的Hessian 矩陣,然后解其特征方程得到全部簡正振動模式的振動頻率.N 越大,計算成本越大.本文提出,針對那些由平衡結(jié)構(gòu)和對稱性就能完全確定的振動模式,可以逐個計算其振動頻率.當(dāng)僅考慮一個振動模式時,3N 維的Hessian 矩陣的計算轉(zhuǎn)化為一維的勢能曲線的計算.基于簡諧振子近似推導(dǎo)單一振動模式下分子勢能曲線的表達(dá)式,接著量化計算勢能曲線,將勢能曲線擬合到表達(dá)式中以獲得振動頻率.相比計算3N 維勢能超曲面及其二階導(dǎo)數(shù)的Hessian 矩陣,僅計算一維勢能曲線而節(jié)省下來的計算資源可以允許選擇更高級別的計算方法和采用更為完備的基組,提高計算的精度.本文首先以計算水分子的B2 振動模式的振動頻率為例,說明了這種方法的可行性.接著將這種方法應(yīng)用到SF6 分子中.多參考組態(tài)相互作用(MRCI)方法是計算電子相關(guān)能的有效方法,本文采用MRCI/6-311G*基組分別計算了SF6 的A1g,Eg,T2g 和T2u 四個振動模式的振動頻率,通過與其他方法的結(jié)果以及實驗結(jié)果相比較,本文計算的四個頻率的相對誤差最小.

1 引言

量化計算是研究分子結(jié)構(gòu)和光譜的重要工具,振動頻率的研究對紅外光譜探測具有重要意義.對于N原子分子,量化計算振動頻率需要計算3N維的分子勢能超曲面及其二階導(dǎo)數(shù)的Hessian 矩陣.多參考組態(tài)相互作用(multi-reference configuration interaction,MRCI)等組態(tài)相互作用方法因使用多個Slater 行列式近似分子軌道,能夠計算更多的相關(guān)能[1,2].預(yù)期能夠得到較精確的振動頻率.對于較大分子,由于電子多、基組大,對計算機(jī)的內(nèi)存和速度都有很高的要求,計算成本較大.本文針對較大分子的簡正頻率計算,提出對3N–6 個簡正振動模式分別計算,將3N維的超曲面及其3N維的二階導(dǎo)數(shù)的Hessian 矩陣的計算轉(zhuǎn)化成一維勢能曲線的計算.這種計算方法大大降低了計算成本.本文將首先采用較低成本的量化計算方法或?qū)ΨQ分析,獲得分子振動的簡正模式,然后針對單一的簡正振動模式,將分子勢能表達(dá)成單一變量且含有參數(shù)振動頻率的函數(shù),進(jìn)而采用MRCI 方法以及較大的基組計算此一維勢能曲線,最后擬合勢能曲線到勢能函數(shù)中獲得振動頻率.

本文在第2 部分將上述方案首先應(yīng)用到研究比較成熟的水分子的振動頻率[3]計算中,來探討這種方法的可靠性.首先量化計算獲得水分子的平衡結(jié)構(gòu).然后在此平衡結(jié)構(gòu)的基礎(chǔ)上,探討3 個簡正振動模式中哪個簡正振動模式可以完全確定,哪個簡正振動模式除考慮對稱性外,還必須考慮力常數(shù).針對可以完全確定的簡正振動模式,進(jìn)行高級別的勢能曲線的計算,推導(dǎo)振動頻率和勢能曲線的關(guān)系,進(jìn)而由勢能曲線擬合獲得振動頻率.最后采用相同的方法,計算基于勢能曲面二階導(dǎo)數(shù)構(gòu)成的Hessian 矩陣,進(jìn)而得到振動頻率,將上面的兩個頻率結(jié)果做比較,確定這種方法的可靠性.

本文在第3 部分將利用此方案計算SF6的振動頻率.SF6是被廣泛使用的高壓電氣設(shè)備絕緣氣體[4?7].SF6分子的振動對SF6分子的高次諧波[8,9]及光電離截面[10,11]等都有直接的影響.人們在實驗上[12?16]對SF6的振動頻率進(jìn)行了大量的測量.理論上,采用LDA(local density approximation)/GGA(generalized gradient approximation)泛函密度、MP2 (second-order Moller-Plesset)方法、RCCSD(T)(restricted coupled cluster theory with singles,doubles,and perturbative triple excitations)等方法都對SF6的振動頻率進(jìn)行了計算[17?19].到目前為止,還沒有看到采用MRCI 方法計算SF6振動頻率的報道.

本文第3 部分安排如下:第1 步,利用分子的對稱性以及高級別的計算,獲得平衡結(jié)構(gòu);第2 步,挑選由對稱性和平衡結(jié)構(gòu)完全確定的簡正振動模式;第3 步,針對第2 步挑選出來的每一個簡正振動模式,一方面推導(dǎo)出勢能與單一變量之間的函數(shù)關(guān)系,另一方面,由MRCI/6-311G*對其進(jìn)行計算,得到勢能曲線,然后將曲線擬合到前面推導(dǎo)的公式中,獲得每個簡正振動模式的振動頻率.

2 水分子振動頻率的計算

2.1 水分子平衡結(jié)構(gòu)的確定

運(yùn)用Molpro[20,21]量化計算軟件,采用MRCI/def2-TVZPP 對水分子做了幾何優(yōu)化計算.獲得水分子的平衡結(jié)構(gòu)為:O—H 鍵鍵長l=1.819 a.u.,鍵角θ=102.755°.

2.2 由分子對稱性和平衡結(jié)構(gòu)確定簡正振動模式的坐標(biāo)

由分子的對稱性和平衡結(jié)構(gòu)可完全確定分子的一部分簡正振動模式的簡正坐標(biāo).其他的簡正振動模式的坐標(biāo)還需要分子力常數(shù).文獻(xiàn)[22]對H2O分子的振動模式的獲得做了詳細(xì)的探討.這里簡述此過程.

H2O 分子屬于C2v對稱群,C2v對稱群有A1,A2,B1和B2四種表示.圖1 給出了水分子的3 個核相對其平衡位置的位移構(gòu)成的九維位移基矢空間,將C2v群的每個操作作用到9 個位移單位基矢上,得到群的一個可約的矩陣表示;進(jìn)一步將其約化成不可約表示,去掉平動和轉(zhuǎn)動的不可約表示,得到代表振動的不可約表示為:A1,A1和B2.下面將借助投影算符來獲得每個振動模式的簡正坐標(biāo).以B2振動模式為例,將不可約表示B2的特征標(biāo)投影算符

作用在9 個位移基矢單位矢量上(位移坐標(biāo)系如圖1 所示)得到

圖1 水分子的位移坐標(biāo)系示意圖,其中1 為氧核,2,3 為兩個氫核Fig.1.Schematic diagram of the displacement coordinate system of water molecule.1 represents oxygen nuclei and 2,3 represent the two hydrogen nucleus.

(1)式等號右邊的四個對稱操作算符依次對應(yīng)C2v群中的變換操作:恒等變換、繞z軸旋轉(zhuǎn)180°、xz面鏡面反射和yz面鏡面反射.可見對稱性是B2的正則振動矢量可以寫作:

當(dāng)中b,c是待定系數(shù),α是任一常數(shù).由質(zhì)心守恒,y方向上有

由角動量守恒條件,x方向上的角動量有

這里mH和mO分別是氫原子和氧原子的質(zhì)量,z1為 O 原子到質(zhì)心的距離.于是c=?2mH/mO,b=cos(θ/2),進(jìn)而可得到對稱性是B2的簡正振動模式為

同樣可以得另外兩個A1簡正振動模式為

其中α′和α′′為任一常數(shù);b1,b2必須由H2O 分子的分子力常數(shù)確定.

由上面求解振動模式的過程可見,分子的振動模式的簡正坐標(biāo)能否由平衡結(jié)構(gòu)和對稱性完全確定,與力常數(shù)無關(guān),取決于是否存在與其具有不同的頻率和相同對稱性的其他振動模式:如果存在,則不能完全確定;不存在,則能.我們將計算可以由平衡結(jié)構(gòu)完全確定的振動模式的本征振動頻率.對于H2O,將計算B2振動模式的本征頻率.

2.3 H2O 分子B2 振動模式頻率的計算

本小節(jié)中,將針對B2簡正振動模式,推導(dǎo)分子勢能曲線的表達(dá)式.然后結(jié)合由Molpro 軟件計算出的勢能曲線,擬合出此振動模式的振動頻率.

2.3.1B2振動模式下的勢能表達(dá)式

根據(jù)波恩-奧本海默(Bonn-Oppenheimer approximation)近似,核在電子的平均勢能場中運(yùn)動.將核的運(yùn)動近似為平衡位置的簡諧振動,則平衡位置附近的勢能ε可以表示為

其中ε0代表分子處于平衡位置的能量;i=1,2,3,分別表示氧核和兩個氫核;Δri表示第i個核相對其平衡位置的偏離.當(dāng)只激發(fā)一個簡正振動模式時,所有核以相同的角速度ω同相振動,振動幅度呈正比關(guān)系,比例系數(shù)由振動坐標(biāo)(如(6)式)給出.由此可見,勢能的表達(dá)式可由一個參量表達(dá).針對水分子的B2振動模式,選取的參量為圖1 中2 號H 核的x方向的位移量,記作Δx2,則勢能表達(dá)為

2.3.2B2振動模式下勢能曲線的量化計算

采用MRCI 方法和/def2-TVZPP 基組,針對B2振動模式,即輸入核坐標(biāo)時,核坐標(biāo)之間的關(guān)系滿足B2振動模式的簡正坐標(biāo)關(guān)系,即滿足(6)式.由此得到水分子B2振動模式的勢能曲線如圖2 所示.圖2 中的橫坐標(biāo)是圖1 中2 號氫核x方向的相對平衡位置的偏離.

圖2 水分子B2 振動模式的勢能曲線圖Fig.2.Potential curve of of H2O molecule for B2 vibrational mode.

2.3.3 振動頻率的擬合計算

將圖2 的數(shù)據(jù)擬合到(8)式中,得到水分子B2振動模式的頻率,為3866.18 cm–1.同時利用MRCI方法,采用def2-TVZPP 基組,直接進(jìn)行頻率計算,即計算整個超勢能面,計算勢能的二階導(dǎo)數(shù)構(gòu)成的Hessian 矩陣,進(jìn)而求得H2O 的3 個振動頻率,其中包括B2模式的振動頻率,為3869.74 cm–1.可見二者之差在幾個波數(shù)的量級上,是基本一致的.通常理論值和實驗值的差別一般在幾十個波數(shù)以上,所以這個差值是可以忽略不計的.

另外,說明一下擬合過程中偏移量的選擇.原則上,偏移量越大,諧振近似越偏離,計算的振動頻率誤差將越大,我們的目標(biāo)是計算分子的平衡點處能量的二次導(dǎo)數(shù),即勢能曲線底部最小值點的曲率半徑.所以將偏移量控制在鍵長的1%的范圍內(nèi),計算表明,在此范圍內(nèi),由于偏移量取值所帶來的偏差會在一兩個波數(shù)的范圍內(nèi).

關(guān)于基組的選擇,由于本文擬合的方法需要計算準(zhǔn)確的平衡位置(這是獲得精確的簡正振動模式的振動坐標(biāo)的要求)和勢能曲線,所以要采用盡可能的完備的基組.那么是否可以由基組外推技術(shù)來進(jìn)一步減小誤差呢? 采用MRCI/ccPVDZ,MRCI/ccPVTZ 和MRCI/ccPVQZ,針對水分子的B2振動模式,進(jìn)行了擬合計算,得到相應(yīng)的振動頻率分別為3877 cm–1,3859 cm–1和3882 cm–1,他們并沒有隨著基組的完備而單調(diào)地趨近實驗值3943 cm–1.我們沒有看到,基組外推能縮小誤差.另外,通過考察計算化學(xué)數(shù)據(jù)庫[3]所給出的修正因子隨基組的完備而變化的情況來探討這個問題.修正因子(scaled factor)是為了使從頭算得到的理論頻率更好地和實驗值相符合,而將理論值乘以一個在0.8 至1.0 之間的因子.如果基組外推,能夠減小頻率的偏差,那么,隨著基組完備性的增強(qiáng),修正因子就應(yīng)該有向1 趨近的趨勢,考察由計算化學(xué)數(shù)據(jù)庫所提供的修正因子隨基組的變化(如QCISD 方法下,基組為3-21G,3-21G*,6-31G,6-31G*,6-31G**,6-31+G**,6-311G*,6-311G**,6-31G(2df,p),6-311+G(3df,2p),對應(yīng)的修正因子依次為:0.969,0.961,0.964,0.952,0.941,0.945,0.957,0.954,0.947,0.954),沒有看到單調(diào)的向1的方向趨近的趨勢.綜上可見,基組外推不能縮小誤差.

3 SF6 分子的振動頻率計算

接下來把上面的擬合計算振動頻率的方法應(yīng)用于計算SF6分子的振動頻率.SF6較H2O 分子復(fù)雜得多.它含有70 個電子,具有15 個簡正振動模式,所以直接采用MRCI/6-311G* 計算21 階的Hessian 矩陣而得到振動頻率是非常困難的.這里將首先采用MRCI/6-311G*確定分子的平衡結(jié)構(gòu),然后針對由對稱性和平衡結(jié)構(gòu)能完全確定的簡正振動模式分別進(jìn)行頻率的擬合計算.

3.1 SF6 的平衡結(jié)構(gòu)

SF6分子是有超高對稱性的正六面體分子,如圖3 所示.我們唯一需確定的構(gòu)型參數(shù)是S—F 鍵的平衡鍵長l0.所以在保持對稱不變的前提下,對6 個鍵長l進(jìn)行相同的改變,逐點進(jìn)行MRCI 水平的單點能計算.首先采用HF/6-311G*獲得電子組態(tài)信息:(core)22(4a1g)2(3t1u)6(2eg)4(5a1g)2(4t1u)6(1t2g)6(3eg)4(5t1u)6(1t2u)6(1t1g)6(6a1g)0(6t1u)0,然后采用完全活性空間自洽場(complete active space self-consistent field)方法確定參考組態(tài).采用的活性空間由1t1g,6a1g和6t1g軌道構(gòu)成,活性電子數(shù)為6 個.最后由MRCI/6-311G*方法獲得勢能隨S—F 鍵長l的變化曲線如圖4 所示,由勢能曲線得到SF6的平衡鍵長為l0=1.558 ?.

圖3 SF6 分子結(jié)構(gòu)示意圖,其中 1 為S 核,2—7 為6 個F 核Fig.3.Schematic diagram of SF6 molecular structure. 1 represents Sulfur nuclei and 2–7 represent Fluorine nucleus.

圖4 SF6 的勢能隨S—F 鍵長l 的變化曲線Fig.4.Potential curve of SF6 as the function of S—F bond length l.

3.2 選擇由平衡結(jié)構(gòu)和對稱性可以完全確定的簡正振動模式

SF6分子所屬Oh群.與水分子類似,在SF6分子的7 個原子核相對平衡位置的位移構(gòu)成的21 維位移空間,約化Oh群的矩陣表示成不可約表示為A1g,Eg,T1u,T1u,T1u,T1g,T2g和T2u,去掉純轉(zhuǎn)動的T1u和平動的T1g表示,獲得SF6分子的15 個簡正振動的對稱性為:A1g,Eg,T1u,T1u,T2g和T2u.A1g是沒有頻率簡并的,Eg是頻率二重簡并的.兩個T1u,T2g和T2u都是頻率三重簡并的.對于兩個T1u振動模式,它們的頻率不同但對稱性相同,其簡正坐標(biāo)是依賴于分子常數(shù)的.所以本文僅針對A1g,Eg,T2g和T2u四種振動模式進(jìn)行頻率計算.

直接使用投影算符獲得SF6分子的簡正振動坐標(biāo)是可行的.但是由于維數(shù)多,所以過程繁雜.這里采用量化計算軟件最經(jīng)濟(jì)的算法HF 方法計算振動頻率.雖然這樣得到的振動頻率的值不是很精確,但是它給出的A1g,Eg,T2g和T2u四種振動模式的簡正坐標(biāo)是準(zhǔn)確的,因為它們是由平衡結(jié)構(gòu)和對稱性確定的(由于SF6的高度對稱性,這四種振動模式的簡正坐標(biāo)與平衡結(jié)構(gòu)也無關(guān)),與量化計算的方法和基組的選取無關(guān).A1g,Eg,T2g和T2u四種振動模式的簡正坐標(biāo)列于表1 中,其中下標(biāo)1 表示S 核,2—7 表示6 個F 核,如圖3 所示.

3.3 SF6 的A1g,Eg,T2g 和T2u 振動模式的頻率ω1,ω2,ω3 和ω4 的計算

3.3.1 各振動模式勢能曲線的表達(dá)式

和水分子情況相似,當(dāng)僅激發(fā)A1g振動模式時,由表1 第1 列的振動坐標(biāo),SF6的勢能為

當(dāng)僅激發(fā)Eg振動模式時,由表1 第2 列的振動坐標(biāo),SF6的勢能為

由表1 第3 列的振動坐標(biāo),SF6的勢能為

當(dāng)僅激發(fā)T2g振動模式時,由表1 第4 或5 或6 列的振動坐標(biāo),SF6的勢能為

當(dāng)僅激發(fā)T2u振動模式時,由表1 第7 或8 或9 列的振動坐標(biāo),SF6的勢能為

表1 SF6 的A1g,Eg,T2g 和T2u 振動模式Table 1.A1g,Eg,T2g and T2u vibrational modes of SF6 molecule.

3.3.2 各振動模式的勢能曲線的從頭計算以及振動頻率的擬合計算

和2.1 節(jié)中獲得勢能曲線的方法相同,采用MRCI/6-311G*獲得各振動模式的勢能曲線,如圖5所示.將勢能曲線數(shù)據(jù)擬合到(11)式、(12)式、(14)式和(15)式,得到A1g,Eg,T2g和T2u振動模式的頻率,具體結(jié)果見表2.表2 還給出了利用HF(Hartree-Fock)、CISD (configuration interaction wavefunction including all single and double excitations),MP2,CCSD(T)和B3LYP (Becke,Lee-Yang-Parr)方法由Hessian 矩陣計算的結(jié)果,同時還給出了相應(yīng)方法計算出的平衡S—F 鍵長l0.我們的目的是比較各個方法的結(jié)果,所以采用的基組都是6-311G*.為方便比較,圖6 給出了上述不同方法計算的4 種振動頻率相對于實驗值的相對誤差.可見,MRCI/6-311G*方法計算的4 個頻率總體來說,相對誤差是最小的;平衡鍵長也最接近實驗值.另外,從圖6 可見,MRCI 方法獲得的4 個頻率的連線與其他方法的結(jié)果相比,最接近平行于橫軸的直線,4 個模式的相對誤差起伏最小,MRCI的結(jié)果與實驗值是最匹配的.

表2 SF6 的A1g,Eg,T2g 和T2u 振動模式的振動頻率ω1,ω2,ω3 和ω4Table 2. Vibrational frequencies ω1,ω2,ω3 and ω4 of the vibrational modes A1g,Eg,T2g and T2u of SF6.

圖5 SF6 不同振動模式的勢能曲線 (a),(b),(c)和(d)分別對應(yīng)A1g,Eg,T2g 和T2u 振動模式Fig.5.Potential curves of different vibrational modes of SF6 molecule:(a),(b),(c) and (d) correspond to A1g,Eg,T2g and T2u vibrational modes respectively.

圖6 不同方法下使用相同基組6-311G*計算的4 種振動頻率相對于實驗值的相對誤差Fig.6.Relative errors to the experimental values of the four vibrational frequencies obtained by different methods and same basis set 6-311G*.

4 總結(jié)

本文給出了一種計算多原子分子振動頻率的節(jié)約成本的計算方法.即,將計算3N維的勢能的二階導(dǎo)數(shù)的Hessian 矩陣分解成計算針對單一振動模式的一維勢能曲線,然后推導(dǎo)勢能曲線的表達(dá)式,通過用表達(dá)式擬合勢能曲線而得到振動頻率.本文的方法適用的振動模式要滿足:與其對稱性相同而頻率不同的振動模式不存在,因為這樣的振動模式的簡正坐標(biāo)可以由平衡結(jié)構(gòu)和對稱性完全確定,與分子力常數(shù)無關(guān).本文計算了水分子的B2振動模式的頻率,說明了這種方法的正確性.進(jìn)而采用這種方法,對SF6的A1g,Eg,T2g,T2u四種振動模式的頻率進(jìn)行了MRCI/6-311G*級別的量化計算.通過和HF,MP2,CISD,CCSD(T)和B3LYP方法以及相同的基組運(yùn)算的結(jié)果以及實驗值相比較,MRCI 方法的結(jié)果的相對誤差是最小的.