高峰 张红 张常哲 赵文丽† 孟庆田‡

1) (山东农业大学信息科学与工程学院, 泰安 271018)

2) (山东师范大学物理与电子科学学院, 济南 250358)

基于Molpro 2012程序包, 应用包含Davidson修正的多参考组态相互作用方法, 使用AVX Z和AVX dZ(X = T, Q, 5, 6)基组进行单点能从头算, 然后采用Aguado-Paniagua函数进行拟合, 得到了SiH+(X1Σ+)离子在不同基组、不同方法和是否考虑自旋-轨道耦合(SOC)情况下的解析势能函数(APEFs).以APEFs为基础,计算了SiH+(X1Σ+)离子的解离能De, 平衡键长Re, 振动频率ωe, 光谱常数Be, αe和ωeχe, 同时讨论了SOC对该体系的影响.本文的计算结果与其他理论计算符合得较好, 与实验数值也基本吻合.基于SOC-AV6dZ方法下的APEF, 通过求解径向薛定谔方程, 给出了SiH+(X1Σ+)离子的前23个振动能级(j = 0), 并详细列出了每1个振动能级及其相应的经典拐点, 每个振动态的转动常数和6个离心畸变常数, 且提供了振动能级图.该工作对于实验和后续的理论工作有参考和指导作用.

1 引 言

分子离子是自然界中普遍存在的物质形态.在星际化学研究中, 分子离子的电子结构和光谱特性对了解在恒星大气和星际空间中发生的物理和化学过程有重要意义[1-6].例如, CH+离子, 作为星云扩散过程中形成的较大碳氢化合物, 研究者们在理论和实验上对其进行了大量的研究[7-9].由于在星际空间中含有丰富的Si和H元素[10], 从理论上来讲, SiH+离子应广泛存在于星际介质和恒星大气中, 实际上, 在太阳光球层的光谱中已观测到了SiH+离子的存在[1].考虑到SiH+离子在星际化学和等离子体物理中重要性[1-4], 不论是在实验还是在理论上对其都予以了足够的关注.

实验上, Douglas和Lutz[2]首先在含硅烷的氦的空心阴极放电发射谱中观察到了SiH+离子的不同电子态 A1Π—X1Σ+跃迁光谱带, 并测量了(1, 0),(2, 1)和(2, 0)的谱线位置.1970年, Grevesse和Sauval[1,3]在太阳光谱中观测到SiH+离子A1Π—X1Σ+的吸收光谱, 确认了太阳光球层SiH+的存在, 同时计算出了SiH+离子 A1Π—X1Σ+光谱带中(0, 0)和(0, 1)跃迁以及SiH+离子 A2Δ—X2Π 光谱带中(0, 0)跃迁的振子强度.1977年, Singh和Vanlandingham[11]在实验室中测量了SiH+基态振动跃迁(1, 0), (2, 1)和(2, 0)的谱线位置, 并计算了其振动能级.Carlson等[12]确认并分析了SiH+离子 A1Π—X1Σ+的(2, 0)和(3, 0)振动跃迁光谱带, 测定了 A1Π 态的ν = 0, 1, 2和3时的寿命, 并把他们转换成A—X跃迁的振子强度的表达形式.Hishikawa和Karawajczyk[13]测定了SiH+离子的(0, 3)和(0, 4)的A—X跃迁光谱带, 并结合(0, 0),(0, 1), (0, 2), (1, 0), (1, 1)[2], (2, 0), (3, 0)[12]光谱带的跃迁数据计算出了SiH+( X1Σ+)的光谱常数.Davies和Martineau[14]使用激光光谱学的方法在硅烷的等离子体中测量了SiH+基态 X1Σ+态的基本振转光谱带.2016年, Mosnier等[15]测量了SiH+离子L-内壳层光谱中主要共振的能量和光解离截面, 并通过组态相互作用的方法对光谱特征进行了理论分析.

理论上, 人们也对SiH+离子开展了一系列高水平的从头算研究.1986年, Hirst等[16]用多参考组态相互作用(MRCI)的方法[17]计算了SiH+离子 X1Σ+, A1Π , a3Π 和 13Σ+态的势能曲线, 并给出了 X1Σ+和 a3Π 的光谱常数.1998年, Matos等[18]采用完全活性空间的自洽场(CASSCF)和多体微扰理论的方法计算了SiH+离子的 X1Σ+态的势能曲线和光谱常数, 同时使用完全活性空间的耦合簇方法和MRCI的方法研究了SiH+离子的 A1Π 态.Sannigrahi等[19]使用大高斯基组, 应用组态相互作用的方法计算了SiH+离 子 X1Σ+和 A1Π 态的光谱常数, 同时计算了激发态 A1Π 的寿命和A1Π—X1Σ+的跃迁振子强度 f00, 并且与实验进了对比.2018年, 应用Molpro 2015程序包[20],Zhang等[21]采用包含Davidson修正的多参考组态相互作用方法(MRCI(Q))计算了SiH+离子最低的7个 Λ —S 态的势能曲线和光谱常数, 计算中使用了AWCV5Z-DK基组.同年, 应用ORCA2.9.0程序[22], Biglari等[23]采用包含标量相对论效应的MRCI方法计算了SiH+离子的基态和低激发态的势能函数, 同时计算了自旋允许的Einstein跃迁系数A和激发态的束缚振动能级的寿命.

尽管已经对SiH+离子的各个电子态进行了大量的理论和实验研究, 但是精确的SiH+离子X1Σ+态的解析势能函数(APEF)、光谱常数、振转能级和自旋-轨道耦合(SOC)等方面仍然缺乏系统的研究.本文基于Molpro 2012程序包[24], 应用MRCI(Q)的方法, 使用aug-cc-pVXZ(AVXZ)和aug-ccpV(X+d)Z(AVXdZ)基组进行能量点的从头计算,其中X = Q, 5, 6.然后应用Aguado-Paniagua函数[25,26]拟合势能曲线函数, 并对AVX Z和AVX dZ(X = Q, 5)基组的能量点应用外插到完全基组集极限(CBS)方法[27,28]进行修正.基于该势能函数,计算了SiH+离子基态X1Σ+态的光谱常数和振动能级, 并研究了SOC效应对该体系的影响.该工作对于实验和后续的理论工作有参考价值和指导作用.

2 理 论

2.1 单点能计算

SiH+( X1Σ+) 离子所有的单点能都采用MRCI(Q)方法进行从头算, 计算采用AVX Z和AVX dZ(X = Q, 5, 6)基组.SiH+离子是双原子离子, 对称点群是 C∞v, 在Molpro 2012程序包中使用其阿贝尔点群C2v子群.在本文的计算中将5个分子轨道确定为活化空间, 分别是3个a1, 1个b1和1个b2分子轨道.Si+原子在3S3P轨道上的3个电子,H原子1S轨道上的1个电子被放在活化空间, 即这4个电子分布在上述5个分子轨道上, 剩下的10个电子放在5个闭壳层轨道, 包括3个a1,1个b1, 1个b2.除此之外, 还有310个外部轨道,分别是108个a1, 76个b1, 76个b2和50个a2分子轨道.SiH+离子的核间距都选择在0.8a0—29a0之间.

2.2 CBS方法

计算电子结构时通常采用1个高斯型轨道展开基组, 基组的不完备性产生基组重叠误差(BSSE)[29,30], 影响计算的准确性.Varandas[27,28]提出, BSSE误差可以通过有限基组单点能计算外推到完备基组(CBS)极限来纠正, 其基于Dunning的相关一致基组外推电子能量的方案如下.

MRCI计算水平下的能量可写成

其中下标X表示计算中使用的AVXZ基组的类型, R是空间坐标的向量, 上标CAS和dc分别表示完备活性空间能量及动态相关能量.

应用Karton和Martin[31]提出的KM方案可得CAS能量

根据Varandas[27,28]提出的统一的单重态和三重态外推方案, 动态相关(dc)能量为

其中系数A5满足下式关系

2.3 势能曲线函数

本文中的势能函数曲线采用的是Aguado-Paniagua函数, 即双原子势能函数 VAB可以表示为短程势和长程势之和[25,26]:

这种函数形式可以保证双原子势能 VAB在RAB→0时为无限大, 同时保证原子势能在解离极限时( RAB→∞ )趋于0.(5)式中的线性参量ai(i=0,1,2,···,n) 和非线性参量 βi(i=1,2) 可以从拟合过程中得到, 本文中取n = 12, 这是综合考虑后拟合的最优结果.

3 结果和讨论

3.1 势能函数

分别使用几种不同的基组和方法(AVQZ, AV-5Z, AV6Z, CBS(Q, 5), AVQdZ, AV5dZ, AV6dZ,CBS(Qd, 5d)), 获得了SiH+(X1Σ+)离子的0.8a0—29a0范围内的112个从头算能量点.基于这些能量点, 应用Aguado-Paniagua函数拟合得到APEF.为了提高精确度, 使用了12个参数, 得到了几种不同基组和不同方法下的APEFs, 函数的具体参数列于表1 (表中只列出了主要讨论的AV6Z, CBS(Q, 5), AV6dZ和CBS(Qd, 5d)).

图1是使用AV6Z基组和CBS(Q, 5)方法得到的SiH+(X1Σ+)离子的从头算能量点(蓝色圆点)和拟合给出的势能曲线(红色实线).两幅图中的上部区域是短程作用范围, 由于在该范围内能量变化范围较大, 因此纵坐标采用对数标度, 横坐标的取值范围是1a0—2a0, 纵坐标的取值范围为10—2Eh—101Eh.中间区域横坐标取值范围是2a0—10a0, 纵坐标的取值范围为0至—0.14Eh, 其内插图为长程作用范围, 相应的横坐标取值范围是6a0—26a0和纵坐标的取值范围为0至—0.005Eh.通过将曲线纵坐标标度的减小和横坐标标度的增大, 清晰地展示出渐近区域内势能曲线和从头算能量点相吻合的特征.底部是拟合的势能曲线与从头算能量点之间的误差, 单位为cm—1.从整体来看, 两种方案拟合得到的APEFs不论是在短程区域还是长程区域都表现出平滑的行为特征, 从头算能量点和拟合曲线之间符合得非常好, 误差均小于10 cm—1.从表1可以看出, 使用CBS(Q, 5)和AV6Z基组的方均根误差(RMSD)为1.61755859 × 10—2和1.60176420 × 10—2kcal/mol.

图1 SiH+(X1Σ+)在CBS(Q, 5)和AV6Z基组下的势能曲线和从头算能量点Fig.1.Potential energy curves and ab initio points at CBS(Q, 5) and AV6Z results.

表1 SiH+(X1Σ+) APEFs的拟合参数Table 1.Parameters of APEFs for SiH+(X1Σ+).

为了提高计算精度, 本文还使用了AVQdZ,AV5dZ, AV6dZ基组及CBS(Qd, 5d)方法计算和拟合了SiH+(X1Σ+)离子的势能曲线.AVX Z和AVX dZ基组二者都属于Dunning的关联一致基组, 而后者在AVX Z基础上加入了扩散函数, 得到了紧凑的扩充集.需要说明的是, 在使用AVQdZ和AV5dZ基组计算时, 对于Si和H元素都是使用AVQdZ和AV5dZ基组; 当使用AV6dZ基组计算时, 由于Molpro程序基组使用的限定, 对于Si和H元素分别使用AV6dZ基组和AV6Z基组计算.

图2给出了分别使用CBS(Qd, 5d)方法和AV6dZ基组所得势能曲线和从头算能量点, 蓝色圆点是从头计算得到的能量点, 红色实线是拟合得到的势能曲线函数.从图2可以看出, 不论是在长程区还是在短程区, 使用CBS(Qd, 5d)方法和AV6dZ基组所得势能曲线光滑且与从头算能量点都符合较好, RMSD较小, 分别为1.60042370 ×10—2kcal/mol和1.627559848 × 10—2kcal/mol.在短程区域, 应用AV6dZ基组的拟合曲线和从头算能量点的误差上限为15 cm—1, CBS(Qd, 5d)方法的误差上限为10 cm—1, CBS(Qd, 5d)方法相对于应用较耗费机时的AV6dZ基组计算更有优势.

图2 SiH+(X1Σ+)应用CBS(Qd, 5d)和AV6dZ基组的势能曲线和从头算能量点Fig.2.Potential energy curves and ab initio points at CBS(Qd, 5d) and AV6dZ results.

3.2 自旋-轨道耦合计算

关于SiH+(X1Σ+)离子的SOC效应的讨论涉及到对应同一解离极限 S i+(2Pu)+H(2Sg) 的 X1Σ+,A1Π , b3Σ+和 a3Π 分子态[19,21].应用AV6dZ基组,在限定这4个分子态的态平均多组态自洽场(SACASSCF)和包含Davidson修正的多参考组态相互作用(SA-AV6dZ)的计算基础上, 再考虑SOC的影响, 所得从头算能量点(SOC-AV6dZ)如图3所示.

图3 SiH+离子的 X 1Σ+ , A 1Π , b 3Σ+ 和 a 3Π 态在SOCAV6dZ基组下的从头算能量点Fig.3.The ab initio points of X 1Σ+ , A 1Π , b 3Σ+ and a3Πstates for SiH+ cation at SOC-AV6dZ results.

基于SA-AV6dZ和SOC-AV6dZ两种方法得到的单点能, 应用Aguado-Paniagua函数进行拟合, 拟合的参数如表1所列.从表1可以看出, SAAV6dZ和SOC-AV6dZ两种方法的RMSD分别为9.45767662 × 10—3和1.11170443 × 10—2kcal/mol,是略小于使用基组AV6Z, AV6dZ以及应用CBS(Q, 5), CBS(Qd, 5d)方法的结果.图4为相应的从头算单点能和拟合势能曲线, 可以看出, 在长程区和短程区域, 使用SA-AV6dZ和SOC-AV6dZ方法所得的势能曲线光滑且与单点能量点符合都比较好, 误差上限分别为10 cm—1和20 cm—1, 其中单纯的态平均计算结果(SA-AV6dZ)更优.

图4 SiH+(X1Σ+)在SA-AV6dZ和SOC-AV6dZ基组下的从头算能量点与拟合势能曲线Fig.4.Potential energy curves and ab initio points at SA-AV6dZ and SOC-AV6dZ results.

3.3 光谱常数

基于不同的基组和拟合方法(AVQZ, AV5Z,AV6Z, CBS(Q, 5), AVQdZ, AV5dZ, AV6dZ, CBS(Qd, 5d), SA-AV6dZ和SOC-AV6dZ)得到的APEFs, 本文计算了SiH+(X1Σ+)离子的平衡键长Re,解离能De, 振动频率ωe, 光谱常数ωeχe, αe, βe.各个光谱常数的公式如下:

其中, f2, f3和f4分别是二阶、三阶和四阶力常数[32,33].

SiH+(X1Σ+)离子的光谱常数列于表2, 表中同时列出了实验[12,34]和其他一些理论研究结果[16,18,21,23].通过对比发现, 本文的计算结果与其他理论计算[16,18,21,23]符合得较好, 与实验数值[12,34]也基本吻合.需要说明的是, 就解离能De而言, 与实验最接近的是使用AVQZ基组得到的结果, 实验值[12,34]De为0.123203Eh, 用AVQZ基组得到的De为0.124640Eh, 与实验值的误差为1.166367%,理论计算值和实验值差别较大.Zhang等[21]计算结果是0.125317Eh, 与实验值的误差是1.715867%;Biglari等[23]的数值是0.124980Eh, 与实验值的误差为1.442335%, 与实验同样偏差较大, 而与本工作的计算结果符合较好.另外, 应用更巨大的基组、更高级的方法和考虑SOC效应(AV5Z, AV6Z,CBS(Q, 5), AVQdZ, AV5dZ, AV6dZ, CBS(Qd,5d), SA-AV6dZ和SOC-AV6dZ), 本应得到更精确的结果, 却发现其理论计算值与实验值的偏差反而增大.基于此, 有理由对实验结果提出疑问, 建议实验工作者重新测定SiH+离子的包含解离能在内的各项光谱常数.

表2 SiH+(X1Σ+)的平衡键长Re, 解离能De, 振动频率ωe, 光谱常数ωeχe, αe和βeTable 2.Spectroscopic constants compared with the experimental values and other theoretical results for SiH+(X1Σ+).

从表2还可知, 基于态平均技术的SA-AV6dZ方法的平衡键长Re、解离能De分别是2.848531a0和0.127264Eh, 考虑SOC效应后的(SOC-AV6dZ)的Re和De分别是2.848382a0和0.126533Eh, 平衡位置基本吻合, 误差为0.005231%, 解离能误差为0.577715%, 偏差较小.基于SA-AV6dZ方法的振动频率ωe, 光谱常数ωeχe, αe, βe分别为2163.448 cm—1, 7.625581 cm—1, 0.216725 cm—1和41.893 cm—1, 考虑SOC效应后的数值分别为2164.033 cm—1, 7.626378 cm—1, 0.217885 cm—1和42.158 cm—1, 误差分别为0.027033%, 0.010451%,0.532391%和0.628588%, 偏差均较小.根据以上数值对比, 发现SOC效应对SiH+离子光谱常数各项的影响较小.

3.4 振动能级

基于SiH+(X1Σ+)离子的APEF, 核运动的径向薛定谔方程[35]为

其中, μ是约化质量, j是转动量子数, v是振动量子数, r是核间距, V(r)是该体系的APEF.在给定振动能级的基础上, 转动能级表达形式如下:

式中, G(v)为体系的振动能级; Bv为体系惯性转动常数; Dv, Hv, Lv, Mv, Nv和Ov分别是为6个离心畸变常数.

表3和表4列出了SiH+(X1Σ+)离子在SOCAV6dZ方法下, 转动量子数j = 0时的前23个振动能级G(v), 经典拐点, 转动常数Bv和6个离心畸变常数Dv, Hv, Lv, Mv, Nv和Ov.理论上说, SOCAV6dZ方法的计算结果由于采用更巨大的基组(AV6dZ)、更高级的方法(SA-MRCI(Q))和考虑SOC效应可以得到最高精度的APEF、光谱常数和振动能级.图5给出了SiH+(X1Σ+)离子在SOCAV6dZ方法下, j = 0时的前23个振动能级图.

表3 SiH+(X1Σ+)离子在SOC-AV6dZ方法下,j = 0时的前23个振动能级G(v)、经典拐点和惯性转动常数BvTable 3.Vibrational levels G(v), classical turn point androtational constant Bv for SiH+(X1Σ+)when j = 0 at SOC-AV6dZ result.

表4 SiH+(X1Σ+)离子在SOC-AV6dZ方法下, j = 0时的前23个振动能级的6个离心畸变常数Dv, Hv, Lv, Mv, Nv和OvTable 4.Six centrifugal distortion constants Dv, Hv, Lv, Mv, Nv和Ov for the top 23 vibrational states of SiH+(X1Σ+) when j = 0 at SOC-AV6dZ result.

从表3可以看出, 对于j = 0的前23个能级,其经典拐点的取值范围较大, 从v = 0的Rmin=2.62981a0, Rmax= 3.11141a0到v = 22的Rmin=2.00734a0, Rmax= 16.60576a0, 这与能级较多相关.v = 22的振动能级(27769.256 cm—1)与解离极限(27770.8 cm—1)非常接近, 二者差值仅为1.544 cm—1.并且v = 22和v = 21的振动能级差值非常小, 为34.794 cm—1, 仅为v = 22的振动能级的0.125%, 以至于在图5的能级图中几乎区别不出这两个振动能级.

图5 SiH+(X1Σ+)离子在SOC-AV6 dZ方法下, j = 0时的前23个振动能级Fig.5.Top 23 vibrational energy levels of SiH+(X1Σ+) when j = 0 at SOC-AV6 dZ result.

由表4可以看出, 在这6个离心畸变常数中, 除Dv范围变化较小外(—3.7712102 × 10—4——5.721315 × 10—3), 其余的5个数值变化范围非常大, 例如, Ov的范围为6.9631 × 10—26— —2.8187 ×10—8, 横跨了18个数量级.所以在表4中只有Dv和Hv把数量级标注在第1行, 其余各常数的数量级都附在数值后.

4 结 论

本文基于Molpro 2012程序包, 应用MRCI(Q)的方法进行单点能从头算, 然后采用Aguado-Paniagua函数进行拟合, 得到了SiH+(X1Σ+)离子的不同基组、不同方法和是否考虑SOC情况

(AVQZ, AV5Z, AV6Z, CBS(Q, 5), AVQdZ,AV5dZ, AV6dZ, CBS(Qd, 5d), SA-AV6dZ和SOC-AV6dZ)下的APEFs.基于APEFs计算了光谱常数De, Re, ωe, Be, αe和ωeχe, 本文的计算结果与其他理论计算[16,18,21,23]符合得较好, 与实验数值[12,34]也基本吻合, 同时讨论了SOC效应对光谱常数的影响.基于SOC-AV6dZ方法下的APEF,通过求解径向薛定谔方程, 首次计算了SiH+(X1Σ+)离子前23个振动态(j = 0), 并详细列出了每一个振动能级及其相应的经典拐点, 每个振动态的转动常数和6个离心畸变常数并提供了振动能级图.该工作对于实验和后续的理论工作有参考和指导作用.