BeC分子基态和低激发态光谱性质和解析势能函数*

张计才 孙金锋 施德恒 朱遵略

(河南师范大学物理与材料科学学院, 新乡 453007)

BeC是一个具有丰富低激发电子态的分子, 本文基于动态权重完全活性空间自冾场方法获得的参考波函数, 采用多参考组态相互作用方法对BeC分子进行高精度的从头计算, 获得了BeC分子 , , ,和 共5个电子态的势能曲线. 为了获得精确的光谱结果, 在计算中考虑了标量相对论效应修正, 并把相互作用能外推至完全基组极限. 在此基础上获得了这些态的光谱常数和偶极距, 以及一些允许跃迁的跃迁偶极距、弗兰克-康登因子和辐射寿命. 最后, 通过扩展的Rydberg函数拟合获得了基态势能曲线精确的解析表达式.

1 引 言

作为在银河系中质量丰度位列第四的元素, 碳在地壳中的含量仅为0.2%, 质量丰度排名为第十五, 而铍是最轻的碱土金属元素, 被广泛应用于制造轻便硬质合金. 铍以及含铍合金具有独特的物理化学特性, 在电子、原子能和航空航天等不同的工业领域中得到广泛应用. 碳化铍在等离子体物理[1,2]、宇宙化学[3-5]以及新材料的开发等领域有着非常大的潜在应用. 但是到目前为止, 仅有少数几个实验和理论工作来研究铍-碳二元团簇[6-15], 因对Be与C相互作用的研究能够为碳化铍二元团簇的研究提供可靠的理论依据, 所以, 对双原子分子BeC的理论研究引起了广泛关注[6-19].

在分子轨道理论中, 强束缚表明成键分子轨道的高占有和较大的键级. 在研究分子的某个激发态时, 该激发态可能与其他电子态相互作用, 这个态的电子组态随着核间距离的改变可能发生剧烈的变化. 由于电子态的相互作用, 离解能将不再适合用于表征键的强度, 而电子态的平衡位置和振动频率则可以. 通常把具有比基态平衡位置小的和谐振频率大的态称之为强束缚态 (strongly-bound state).如果分子电子态的能量低于相应态的离解极限, 则称该电子态为束缚态. 如果该电子态的总能比所有可能的离解通道的总能低, 则称为热力学稳定态.1993 年, Wright和 Kolbuszewski[7]通过研究发现BeC分子具有强束缚双激发态的特性, 并给出了这些低激发电子态的势能曲线和光谱常数; 同年,Borin和Omellas[8]使用多参考单双激发组态相互作用方法研究了BeC分子的电子结构和光谱性质,计算了13个低激发电子态的势能曲线, 给出了这些态的光谱常数、激发能、振动能级以及偶极距等数据. 在 2004 年, Pelegrini等[10]利用 6-311+G基组并使用多参考组态相互作用 (multi-reference configuration interaction, MRCI) 方 法 研 究 了BeC分子第一离解极限下的基态和激发态两个电子态, 获得了平衡位置、谐振频率以及离解能等相关光谱常数. Teberekidis等[19]使用单双激发的多参考布里渊-维格纳耦合簇 (multireference Brillouin-Wigner coupled cluster method restricted to single and double excitations,MRBW-CCSD) 方法, 通过外推至完全基组极限,获得了基态的平衡位置、离解能等光谱常数. BeC分子有非常低的电子激发态, 该激发态大约高于基态 2300 cm–1, 因此, 利用单参考方法研究该体系将是一个挑战. 实验上已经获得BeC分子的基态转动常数 B0= 1.162 cm–1, 但由于势的复杂形状和实验数据的稀缺, 没有获得平衡位置等其他光谱常数的信息[6].

完全活性空间自洽场方法在确定分子激发态的电子结构, 尤其是具有多参考特征的激发态等方面是一种非常重要的方法. 然而, 当一个避免交叉态未包含态平均时, 计算将导致另一个态的势能曲线存在不连续行为. 由于分子系统随着能量的增加态密度也增加, 这时在进行态平均计算时采取增加越来越多态的方法有时是无效的. 为了处理这样的问题, Deskevich 等[20]和 Dawes等[21]提出了一种动态权重方法, 即在完全活性空间自洽场 (complete active space self-consistent feld, CASSCF)

计算中通过自动改变态的平均权重来获得光滑的势能曲线. 本文采用动态权重态平均完全活性空间自洽场 (dynamically weighted state averaged complete active space self-consistent feld, DW-SACASSCF) 方法得到的函数作为利用内收缩多参考组态相互作用 (internally contracted multi-reference configuration interaction, icMRCI)[22,23]方法进行计算的参考波函数, 研究了BeC分子第一和第二离解极限下 5 个电子态 (,,,和) 的光谱和跃迁性质, 并对基态的解析势能函数进行了有益探索.

2 计算方法

由 Deskevich等[20]和 Dawes等[21]所提出的DW-SA-CASSCF能量表达式为

式中 ECAS,i是第 i个态的 CAS能量,是这个态的权重, 该参数是所有CAS态能量的显函数, 求和遍布所有态平均中的N个态. 能量相关权重函数定义为[20,21]:

对权重参数的选择需要避免分子态以及分子态的权重发生快速变化, Deskevich等[20]和Dawes等[21]推荐使用的双曲正割函数则很好地满足了要求. 函数中衰减系数 β 依靠经验进行优化, 通常取值为 3 eV. 程序开始运行时, 所有的态具有相同的权重, 然后进行第二次迭代计算, 此时各个态的权重由(2)式来确定, 反过来也会稍微改变各个态的能量值, 直到SA-CASSCF能量达到自洽时迭代停止 (一般进行3—6次迭代). 当误差改变小于指定的值 (大多数情况下,10–6a.u.是一个合理的选择, 此时随后的 MRCI计算收敛将好于 10–7a.u., 1 a.u. = 2625.4984 kJ/mol)时, 从DW-MCSCF计算中获得的轨道被当作随后icMRCI计算的参考轨道. 目前, 尽管动态权重程序已经内置于Molpro中, 但在输入文件中写入一个简单的循环结构将更便于应用.

所以, 在应用Werner等[24]开发的量子化学从头计算Molpro2010程序包计算BeC分子的电子结构时, 本文采用动态权重方法来获得MRCI的参考态. 为了更精确地预测分子结构(0.001 Å)或能量 (小于 1 kcal/mol, 1 kcal/mol = 4.184 kJ/mol),对于第一、第二主族的原子, 必须以适当的方式考虑内层的非价电子对体系的影响[25]. 所以本文直接使用考虑了核价相关效应的aug-cc-pCVQZ和aug-cc-pCV5Z基组[26,27], 并在此基础上将能量外推至完全基组极限. 计算时取核间距从1.18—11.8 Å之间的60个点做单点能计算.

分子电子态的总能由参考能和相关能分别进行外推得到. 参考能采用Karton-Martin外推公式[28]

表1 BeC 分子 6 个态的离解极限关系Table 1. Dissociation relationship of six electronic states of BeC molecule.

相关能使用Halkier等[29]获得的两点外推公式

在进行标量相对论修正计算时, 常采用二阶Douglas-Kroll-Hess (DKH2) 哈 密 顿 近 似 . 由 于Hirao及其合作者[30-32]使用三阶DKH3 获得了与实验更接近的结果, 因此考虑BeC分子的相对论效应修正时, 本文使用三阶DKH3 哈密顿近似进行相对论修正(记为+DK).

3 结果和讨论

3.1 势能曲线和光谱参数

Be原子的基态和第一激发态光谱项为1Sg和3Pu. 如表1所列, BeC分子前三个离解通道是由基态Be原子和C原子的三个最低价态3Pg,1Dg和1Sg构成. Barker等[6]通过第一原理计算获得Be原子第一激发态与基态的能级差是21846 cm–1(加 Davidson 修正) 或 22033 cm–1(不加 Davidson修正), 发现不加修正的能量更接近于实验值21978.3 cm–1[33], 更高离解通道 Be(3Pu)+C(1Dg)的计算也表明未加修正的能量更接近实验值. 本文通过进一步计算发现 (见表1), 第三离解极限Be(1Sg)+C(1Sg) 也存在同样的情况, 但对于第二离解极限Be(1Sg)+C(1Dg) 包含Davidson修正的能量 (10169.70 cm–1)比不包含 Davidson修正的能量(10124.12 cm–1)更接近于实验值(10192.66 cm–1)[34].所以, 在本文中仅研究了前两个离解通道对应的5个电子态.

图1是利用 icMRCI+Q/Q5+DK 方法计算获得的BeC分子5个态的势能曲线. 由于态与更高的态 (也即态) 存在避免交叉, 使态具有了双势阱结构的势能曲线, 同时也表明这2个对称性相同的电子态之间存在较强的相互作用, 为了得到准确的结果, 在计算中必须同时对这两个电子态进行优化. 从图1可以看出, 这5个电子态均为束缚态, 都具有非常深的势阱. 使用LEVEL8.0[35]程序对这些势进行拟合可以获得相应的光谱常数, 具体数值见表2. 为了便于对比, 已有的实验和理论计算结果也列在表2中. 对于BeC分子的基态, 其特征组态为. 我们预测的势阱深度为 2.1873 eV, 势阱能够支持34个振动束缚态. 实验报道了基态振动谱项= 902 cm–1, 我们计算的结果为 911 cm–1,与实验符合得很好. 计算的 B0= 1.165 cm–1也与相应的实验值1.162 cm–1非常接近, 同时我们的光谱结果也与其他理论结果相符合.

图 1 BeC 分 子 5 个 电 子 态 的 势 能 曲 线 (1 Hartree =2625.4984 kJ/mol)Fig. 1. Potential energy curves of five states of BeC molecule (1 Hartree = 2625.4984 kJ/mol).

表2 BeC 分子 , , , 和 等 5 个态的光谱参数Table 2. Spectroscopic constants of the five states of BeC molecule.

表2 BeC 分子 , , , 和 等 5 个态的光谱参数Table 2. Spectroscopic constants of the five states of BeC molecule.

注: Cal. 为理论计算值.

0 0.1673 918.08 7.350 17.87 1.1783 15.644 2.1873 Cal. [6] 0 0.1661 937.9 9.6 — 1.19 — —Cal. [7] 0 0.1693 905 — — — — 2.04 Cal. [8] 0 0.1667 951 8.42 — 1.183 — 2.39 Cal. [10] 0 0.1683 925 11.25 — — — 2.14 Cal. [19] 0 0.1680 — — — — — 2.04/cm–1 De/eV Te/cm–1 Re/nm /cm–1 /cm–1102 уe/cm–1 Be/cm–18916.35 0.1771 772.74 8.692 199.89 1.0518 28.367 1.0777 Cal. [7]9033.41 0.1799 — — — — — 0.92 Cal. [8] 9466 0.1756 764 14.69 — 1.0752 — 1.16 Cal. [10] 8961 0.1791 874 26.26 — — — 1.037823.39 0.1675 933.50 8.301 7.56 1.1714 16.700 2.4408 Cal. [7]8872.100.1693 904 — — — — 2.27 Cal. [8] 8732 0.1668 956 7.6 — 1.1757 — 2.6312139.14 0.16698 936.21 17.553 126.98 1.1764 16.000 1.9060 Cal. [7]12582.240.1693 905 — — — — 1.81 Cal. [8] 13579 0.167 955 7.3 — 1.1732 — 2.02 10909.21 0.1760 834.80 7.138 9.20 1.0551 14.365 2.0933 Cal. [7]11291.760.1778 818 — — — — 1.97 Cal. [8] 11618 0.1758 847 6.94 — 1.0612 — 2.24

第二个解离通道 Be(1Sg)+C(1Dg) 对应三个电子态:,,. 从图 1 可以看出, 这三个态也都有很深的势阱. 这三个电子态的势阱能够容纳的振动态数目分别为37, 34 和28. 通过对它们的波函数分析发现在平衡位置附近与态相似, 也具有单组态特征, 而和两个电子态则具有显著的多组态特征.

3.2 跃迁特性

图2 BeC 分子 5 个电子态的偶极矩 (1 Debye = 3.336 ×10–30 C·m)Fig. 2. Dipole moments of five electronic states of BeC(1 Debye = 3.336 × 10–30 C·m).

采用icMRCI方法, 在cc-pCVQZ基组的基础上计算了BeC分子5个电子态在核间距R = 0.1—0.8 nm的永久偶极距, 结果如图2所示. 从图2中可以看出, 当核间距离较小时, 5个态的电偶极距均为负值, 表明分子极性为, 即Be失去电子, 而C得到电子. 这从图3基态的Mulliken电荷布局可以更加明显的展示出来: 当核间距较小时, Be带正电, C 带负电, 而当核间距比较大时,Be和C两个原子将均不带净电荷, 电偶极距趋于零值, 同时, Be和C原子的Mulliken电荷也分别为零.

BeC分子激发态到基态的跃迁特性, 对于开展该分子光谱实验研究具有非常重要的指导意义.本文计算了BeC分子核间距在0.1—0.8 nm的—,—以及—等三个允许跃迁的跃迁偶极矩, 计算结果如图4所示.—有一个突变源自与一个更高的激发态之间的避免交叉引起的. Wright和Kolbuszewski[7]以及 Borin和Ornellas[8]在 1993年通过计算均发现了这个现象. BeC分子的低激发态都有较深的势阱, 使得该分子具有非常丰富的光谱信息 .—以 及—之 间 的 跃 迁 偶 极矩随着核间距的变化发生跃变, 这是由于和更高的激发态的避免交叉导致波函数产生组态交换引起的.

图3 BeC 分子基态的 Mulliken 电荷分布Fig. 3. Mulliken population for the ground state BeC molecule.

图4 跃迁偶极距曲线Fig. 4. Transition dipole moment curves.

根据前面的弗兰克-康登因子, 通过(5)式可以计算这些振动态的自发辐射寿命

3.3 解析势能函数

目前, 获得分子势能数据主要是通过第一原理计算、Rydbery-Klein-Rees(RKR)反演法和逆向微扰法等三种途径, 但得到的都是离散的数据点. 本文通过研究已有的各种形式的解析势能函数, 其中下面这种扩展的Rydberg函数[37]在相同参数的情况下给出了最优的拟合结果

式中De是离解能, Re为分子的平衡位置, Ck为可调参数. 为了使均方根偏差达到较小的值, 拟合曲线时需要平衡考虑精度和最少可调参数

表3 — , — , — 跃迁的弗兰克-康登因子Table 3. Franck-condon factor for − , − , − transitions.

表3 — , — , — 跃迁的弗兰克-康登因子Table 3. Franck-condon factor for − , − , − transitions.

—FC—FC —FC — FC —FC —FC —0—0 0.5439 0—1 0.3653 1—0 0.3049 1—2 0.4143 1—3 0.1670 2—0 0.1080 2—1 0.2660 2—3 0.3601 2—4 0.2276 3—2 0.1654 3—4 0.2908 3—5 0.2691 4—2 0.1763 4—5 0.2255 4—6 0.2993 5—3 0.1561 5—6 0.1639 5—7 0.3216 6—3 0.1070 6—4 0.1273 6—6 0.1339 6—7 0.1052 6—8 0.3348 6—9 0.1183 7—4 0.1095 7—5 0.1000 7—7 0.1618 7—9 0.3342 7—10 0.1531 —0—0 0.6237 0—1 0.3153 1—0 0.2701 1—1 0.1740 1—2 0.3999 1—3 0.1368 2—1 0.2689 2—3 0.3470 2—4 0.2322 3—1 0.1518 3—2 0.1803 3—4 0.2122 3—5 0.2964 4—2 0.1894 4—4 0.1027 4—5 0.1045 4—6 0.3042 4—7 0.1238 5—3 0.1634 5—5 0.1435 5—7 0.3138 6—3 0.1236 6—4 0.1268 6—6 0.1567 6—8 0.2968 7—4 0.1289 7—7 0.1495 7—9 0.2690 8—5 0.1219 8—8 0.1212 8—10 0.2380 9—11 0.2025 10—12 0.1699—0—0 0.6164 0—1 0.2717 1—0 0.3232 1—1 0.1693 1—2 0.2617 2—1 0.4016 2—3 0.1699 2—4 0.1850 3—1 0.3901 3—2 0.3513 3—5 0.1510 4—2 0.2362 4—3 0.2156 4—4 0.1073 4—6 0.1068 5—3 0.3014 5—4 0.1011 5—5 0.1501 6—4 0.3109 6—6 0.1613 7—4 0.1292 7—5 0.3175 7—7 0.1442 8—5 0.1801 8—6 0.2945 8—8 0.1054 9—6 0.2192 9—7 0.2612 10—7 0.2563 10—8 0.2265

表4 , , 和 态几个振动能级的辐射寿命Table 4. Radiative lifetime of the vibrational energy levels for , , and states.

表4 , , 和 态几个振动能级的辐射寿命Table 4. Radiative lifetime of the vibrational energy levels for , , and states.

= 7跃迁辐射寿命 = 0 = 1 = 2 = 3 = 4 = 5 = 6 — 191.10 166.59 161.38 159.48 151.88 176.09 197.32 239.85—14.66 14.64 14.46 13.82 12.99 12.16 11.65 11.31 Cal.[10] 14.0 14.5 15.6—2730.00 1010.00 630.00 460.00 380.00 330 280.00 230.00

表5 BeC 分子基态解析势能函数参数值Table 5. The values for the analytic parameters for the ground state of BeC molecule.

式中Ee(Ri) 表示的是核间距Ri处从头计算的势能值, V(Ri) 表示的是核间距Ri处解析表达式得到的数值, N是需要拟合的从头计算势能数据点的数目. 进行拟合前, De和 Re取自于从本文的理论计算值.

为了精确地获得拟合参数, 本文利用1stOpt 6.0软件[38]独特的稳健全局优化算法对BeC分子基态的势能数据点进行了曲线拟合. 当取10个可调参数时, 如表5所列, 得到的均方根偏差已小于 5.03 × 10–5a.u., 表明该解析势能函数可以精确描述基态势能曲线.

4 结 论

利用高精度的从头计算MRCI+Q方法, 本文研究了BeC分子第一和第二离解极限对应的基态和 4 个 低激发 态 (,,和)的光谱性质, 计算得到这五个态的势能曲线, 并使用LEVEL程序拟合得到了相应的光谱常数. 计算的结果与最新报道的实验值符合得很好. 当前的计算还给出了5个电子态的永久偶极矩, 结果表明在这几个激发态的Re处BeC的极性均为.本文还讨论了偶极允许跃迁—,—,—的跃迁特性, 其中包括相应的跃迁偶极矩、Franck-Condon因子以及辐射寿命. 结果显示:—和—自发跃迁的辐射寿命处在量级,—自发跃迁的辐射寿命 则 为 ms量 级 . 最 后 , 我 们 还 通 过 扩 展 的Rydberg函数拟合获得了基态势能曲线精确的解析表达式.