CO分子四个电子态的振转谱:两种效应修正方法的比较∗

2018-12-02徐慧颖刘勇李仲缘杨玉军闫冰

物理学报 2018年21期

徐慧颖刘勇李仲缘杨玉军闫冰

1)(吉林大学原子与分子物理研究所,长春 130012)2)(吉林大学大数据和网络管理中心,长春 130012)(2018年8月1日收到;2018年8月29日收到修改稿)

基于完全活性空间自洽场方法和多参考组态相互作用(multi-reference configuration interaction method,MRCI)方法,采用MRCI+Q/CBS(TQ5)+CV+SR(方法A)和aug-cc-pwCVnZ-DK(n=T,Q,5)(方法B)方案,分别计算了包含Davidson修正(+Q)、芯-价电子关联(core-valence correlation correction,CV)效应以及标量相对论(scalar relativistic,SR)效应的CO分子的基态X1Σ+和激发态a3Π,a′3Σ+和A1Π的势能曲线.在此基础上,获得了这些电子态的振-转谱.通过与实验结果比较发现:方法A适合a′3Σ+和A1Π等较高激发态的振-转谱的计算,方法B更适合基态X1Σ+和第一激发态a3Π的振-转谱的精细计算.该研究可以为其他小分子高精度振-转谱快速计算方案选择提供参考.

1 引言

分子内部精细的电子结构研究对探索物质的性质、结构和功能具有重要意义.对分子的微观结构与光谱性质的计算一直是原子分子物理学科的重要内容之一.目前,多电子分子体系的准确计算,采用的主要方案是多参考组态相互作用方案.由于需要较多的计算资源,对于简单的分子体系准确的全电子计算也存在较大困难.分子的外壳层性质对于诸多物理化学过程具有更重要的意义,为了实现对其准确快速计算,计算中通常将靠近原子核的电子作用采用一个有效势来代替,计算后再将价电子与芯结构之间的关联作用加到整体能量计算中.此外,为了得到更准确的能量,需要考虑相对论效应.为了准确计算该效应需要求解Dirac方程,其计算更为困难,相对论效应按性质可分为标量相对论(scalar relativistic,SR)效应(达尔文项和质量-速度项)和矢量相对论效应(自旋-轨道耦合项),本文仅考虑了SR效应.

对于SR效应的修正,人们在研究中采用了许多全电子近似的相对论方法.例如,Jong等[1]计算了Hartree-Fock(HF)水平的CF4,Si4和Br2CO的SR效应对平衡核间距Re及总解离能ΣDe的改变.对于芯-价电子关联(core-valence correlation correction,CV)效应,人们针对第二行的原子Al-Ar的研究[2]使用了关联效应的相关一致基组.对这两个效应的修正可以通过如下两种方案进行.

1)方法A

针对这两个修正效应贡献于体系的相对能量特点,采用不包含这两个相互作用的基组较大基函数,计算体系中不包含这两个相互作用的能量,然后通过选择新的较简单的基函数快速计算其贡献.2014年,Abbiche等[3]通过多参考组态相互作用(multi-reference configuration interaction,MRCI)+CV+SR方法对CP分子的X2Σ+和A2Π等7个较低电子态进行了计算.

2)方法B

在基函数选择时选择包含这两个相互作用的基函数,计算体系哈密顿中包含这两个效应的总能量.2014年,Li等[4]通过aug-cc-pwCVnZ-DK基组计算了GeH+的8条电子态的势能曲线(potential energy curves,PECs)和光谱常数等.

这两种方案均可以给出较准确的计算结果,并能实现较快速的计算.但对于其适应的条件还需要深入的研究,即使对于较简单分子体系的不同电子态的计算也是如此.

为了探索这两个效应修正方案的差异,本文选择CO分子加以研究.CO的研究在等离子体、星际物质、生物医学以及环境检测与分析具有重要意义.过去的数十年间,科研人员已对CO分子进行了大量的研究[5−21].2013年Lu等[22]采用aug-cc-pV5Z基组计算了CO分子较低的单重和三重态(包括X1Σ+和A1Π等6个电子态)的PECs和基态X1Σ+在振动量子数ν=20之前的振-转能级结构,其计算中并没考虑到CV和SR效应修正.同年,Shi等[23]选用基组aug-cc-pV5Z计算了MRCI水平的CO分子8条较低电子态的PECs和振-转谱,同时分别采用基组为cc-pCVQZ和cc-pV5Z考虑了CV效应和SR效应对于PECs的影响(aug-ccpV5Z+CV+DK).通过两种方案,计算了两种效应对CO分子能量的影响,进而分析了分子的振-转谱.研究发现方法A适合较高激发态的振-转谱的计算,方法B适合CO分子基态和第一激发态的振-转谱的精细计算.

2 方法和计算细节

为了得到精确的CO分子的PECs和振-转谱,使用了molpro2012[24]软件包计算了分子的电子结构.为了实现效率和精确度的平衡,首先选择HF方法获得分子体系的波函数;然后在此基础上利用态平均的完全活性空间自洽场(complete active space seif-consistent field,CASSCF)方法[25,26]来描述电子之间的静力学相关效应;最后采用MRCI方法[27,28]将电子间的动态相关效应也包含进来.

考虑CO分子的基态X1Σ+和激发态a3Π,a′3Σ+和A1Π,根据其对称性,分子轨道从头计算选择在C∞v的子群C2v中进行.C∞v和C2v的不可约表示的对应关系分别为:Σ+-A1,Π-B1+B2,∆-A1+A2,Σ−-A2.在计算这4个电子态的PECs时,选取了一系列单点能量计算,计算的核间距区域R和间隔∆R分别为R=0.8—2.5 Å,∆R=0.01 Å;R=2.5—5.2 Å,∆R=0.1 Å (1 Å=0.1 nm). 每步的具体计算如下.首先,对每个单点进行HF自洽场计算,在HF计算中,采用了单组态Slater行列式描述CO的电子基态,此处只考虑了自旋反平行电子的相关作用.然后,以HF方法产生的分子轨道作为初始轨道,进行CASSCF计算,分别获得CO分子的基态X1Σ+和激发态a3Π,a′3Σ+和A1Π的波函数.在CASSCF计算中选择CO的10个分子轨道和10个n=2电子作为活性空间,包括4个a1,两个b1和两个b2对称性的分子轨道,它们对应C原子的原子轨道2s2p和O原子的轨道2s2p.C的外层电子2s22p2和O的外层电子2s22p4被放置在活性空间内,剩下的4个电子被冻结而不进行相关能计算.最后,进行MRCI计算,但由于大小一致性误差的存在,在MRCI中加入了Davidson修正(+Q)[29],估计了电子的四重激发修正,这样就保证了在无穷远分子的总能量近似等于单独计算C和O原子的能量之和.

对于CV效应和SR效应的修正通过两种方案进行.方法A首先采用基组aug-cc-pVnZ(n=T,Q,5)[30,31]通过HF,CASSCF和MRCI方法计算CO分子4个束缚态的PECs.为了讨论内层电子相关效应的影响,在MRCI计算中,考虑了C和O原子的1s2内壳层的单、双电子激发产生的CV效应.在不考虑CV效应的计算中,将上述的单、双电子激发关联轨道冻结.在上述所有计算中,为了考查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个束缚态的PECs.选择的基组包含了这两种效应修正的基函数,计算体系哈密顿中包含这两个效应的总能量.为了减小基组带来的误差,对两种修正方案的基组(n=T,Q,5)外推到完全基组(complete basis set,CBS).

最后,根据方法A和方法B计算获得的X1Σ+,a3Π,a′3Σ+和A1Π态的PECs,通过LEVEL[34]程序拟合求解CO分子的一维径向Schrödinger方程,得到各个电子态的振-转谱.

3 结果和讨论

利用方法A和方法B计算方案得到了CO分子基态X1Σ+和激发态a3Π,a′3Σ+和A1Π的PECs.它们对应的解离极限是C(3P)+O(3P).根据方法A计算的PECs(图1),利用了LEVEL程序对4个电子态进行拟合得到光谱参数,通过比较发现其值与实验值十分接近,特别是对于转动常数Be和平衡核间距Re的计算,与实验值的相对误差均小于0.0001 Å和0.001 cm−1.

图1 方法A计算的基态X1Σ+和激发态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计算方案得到的基态X1Σ+和激发态a3Π,a′3Σ+和A1Π 的PECs,可以计算其振动能级Gν. 表1根据两种计算方案列出了基态X1Σ+和第一激发态a3Π的前21个Gν,从表1可以看出,根据方法A,在ν=0—10时(转动量子数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%.根据方法B,在ν=0—10时(J=0),振动能级Gν的计算结果范围为1082.11—21344.55 cm−1,相对偏差为0.33—13.41 cm−1,相对误差范围为0.03%—0.06%.在较低振动态时方法B比方法A的计算结果更接近实验测量值[35]. 在ν>10时,方法B的计算结果也更精确,Gν计算结果范围为23228.85—39143.40 cm−1,对应的相对偏差为15.67—144.53 cm−1. 通过方法B计算的Gν的RMSE为0.13%.

对于第一激发态a3Π,实验[36]上只给出了a3Π态的前8个振动能级(ν=0—7)的测量值. 根据方法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%.通过数据对比发现,对于基态X1Σ+和第一激发态a3Π,方法B比方法A计算的振-转谱更符合实验测量值.

通过方法A和方法B计算的激发态a′3Σ+和A1Π的PECs拟合出前21个振动能级列于表2.表2同时列出了相应的实验测量值[36].对于a′3Σ+态,方法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Π态, 方法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%.通过分析上述结果可知,对于激发态a′3Σ+和A1Π,方法A方案更接近于实验的测量.

为了更清楚地对比两种计算方案得到的振-转谱的计算精度,计算了模拟结果与实验结果的相对误差,如图2所示.从图2可以看出,两种方法计算的振转能级与实验的相对整体都小于1%.但是对于不同电子态还存在一定差距,对于基态X1Σ+和第一激发态a3Π,方法A的计算结果相对误差大于方法B,而对于激发态A1Π和a′3Σ+,方法A的计算结果优于方法B.这一计算进而可以推广到更高激发态的计算.通过对上述结果分析可知,方法B在波函数中包含修正效应,具有更好的普适性,通常对于小分子计算可以采用该方案.而对于较重分子,可以应用方法A,可以在相对小的计算需求条件下得到较准确的结果.本文采用不同的方案计算了SR效应与CV修正,并且后者对精密计算是不可或缺的.对于最低的两个电子态,考虑了两种效应对计算高斯基组的依赖性(方法B),发现振动-转动的计算精度提高了,可见,对于这两个电子态对电子关联计算的要求较高;对于更高的电子态,电子云分布相对松散,反而简单地采用单一高斯基组获得的电子关联修正即能达到相应的计算准确性;当然,由于振-转谱的计算实质上只有相对能量起作用,这里也包含了不同电子态电子关联效应的抵消效应.

表1 通过方法A和方法B计算得到的CO分子X1Σ+和a3Π态的振动能级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计算基态X1Σ+和激发态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 ν.