HCP分子势能面和振动光谱研究

张佳美 徐建刚 张云光 张思雨 郭甜 陈艳南

摘 要    ：本文基于Molpro2019程序包， 采用单双激发耦合簇方法（CCSD）结合基组cc-pVQZ构建HCP分子的一维势能曲线和二维势能曲面， 发现HCP分子面内、面外弯曲振动模式之间的简并现象以及H-C拉伸振动模式对分子势能的重要影响. 本文以势能面为基础， 首次采用振动多组态自洽场方法（VMCSCF）和振动多参考组态相互作用方法（VMRCI）计算HCP分子的基频、倍频、组合频以及振动能量， 频率计算值与实验值吻合较好. 拟合绘制出分子的红外和拉曼振动光谱， 发现振动模式间的费米共振现象. 本文为含磷星际分子的实验和理论研究提供了参考.

关键词 ：HCP分子; 势能面; 振动多参考组态相互作用方法; 振动频率; 振动光谱

中图分类号 ：O561 文献标识码 ：A DOI ：  10.19907/j.0490-6756.2023.044001

The potential energy surface and vibrational spectroscopy of HCP molecule

ZHANG Jia-Mei， XU Jian-Gang， ZHANG Yun-Guang，

ZHANG Si-Yu， GUO Tian， CHEN Yan-Nan

（School of Science， Xian University of Posts and Telecommunications， Xian 710000， China）

In this paper， based on the Molpro2019， the one-dimensional potential energy curves and two-dimensional potential energy surfaces of HCP are constructed by using the single and double excitation coupled cluster method（CCSD）combined with the basis set cc-pVQZ. It is found that the degeneracy between the in-plane and out-of-plane bending vibration modes of HCP molecule and the important influence of the tensile vibration mode on the molecular potential energy. Based on the potential energy surface， the fundamentals， overtones， combination bands and vibration energy of HCP are calculated by the vibrational multiconfiguration self-consistent field method （VMCSCF） and vibrational multireference configuration interaction method （VMRCI）. The calculated vibration frequencies are in good agreement with the experimental values. The infrared and Raman spectrum of HCP are further fitted and plotted， and the Fermi resonance between vibration modes is found. This paper provides  reference for the experimental and theoretical research of phosphorus-related interstellar molecules.

HCP molecule; Potential energy surface; Vibrational multireference configuration interaction method; Vibrational frequency; Vibrational spectrum

1 引 言

HCP分子是星际空间中能探测到的少数含有磷的分子之一. 在PN和CP之后， 它是第三个被发现的含磷星际分子. 研究含磷星际分子对于天体物理和星际磷化学研究具有重要意义， 因此受到专家学者的密切关注.

在实验观测方面， 1961年Gier  ［1］首次成功在实验室条件下制备出了HCP固体聚合物单体， 确定了单体HCP制备的比例， 并在-196°处检测到固体HCP单体的红外光谱与HCP结构完全一致， 为HCP结构的研究提供了进一步支持. 随后， Garneau和Cabana  ［2-5］陆续用不同分辨率的光谱仪实验观察到HCP振转光谱的各个谱带， 并拟合了指定红外光谱的所有波数. 2007年， Agúndez和他的团队  ［6］使用IRAM 30 m望远镜在AGB星周包层IRC+10 216 cm波段处探测到HCP的存在. 这是HCP分子首次在星际介质中被探测到. 在此之前， 研究用到的HCP均基于Gier  ［1］的制备方法制备， 但该方法对制备温度要求严格， 且制备纯度有限， 从而导致实验观察到的HCP分子的红外光谱并不准确.

在计算方面， 随着量子化学领域的发展， 关于HCP分子的量化计算也同步开展  ［7-10］. 上世纪末至今， 对于HCP分子的量化研究内容主要包括几何平衡构型、动力学特性、振动能级、势能面、振动谱等. 1993年， Chen和Chong  ［11］采用deMon密度泛函方法计算了基态HCP分子的势能面， 并使用FORTRAN 77环境下的TRIATOM程序和DVR程序计算出其振转能量， 将理论值与实验值进行了比较. 但该方法存在较大的局限性， 计算精度较差. 1996年， Koput等人  ［12，13］使用耦合聚类方法和微扰法对HCP分子展开计算， 先后发表了两篇论文报道了计算得到的几何结构参数、光谱常数、势能面和振动能级. 2000年， Koput等人  ［14］又使用内收缩多参考态组态相互作用方法对HCP分子展开当时最为精确的量化计算， 提出了HCP分子第一个全局精确势能面， 从而提供了基态HCP分子的完整内部动力学描述. 此外， 围绕HCP分子， 众多学者还广泛开展了对其高激发态振动动力学、同位素分子以及HCP-CPH异构化过程的研究. H ltzl等人  ［15］使用了多种量子化学方法， 包括CCSD（T）、CASSCF、CASPT2、MR-ACPF和MR-ACPF-2， 研究HCP二聚化体系. 2016年和2020年， Wang等人  ［16，17］先后报道了采用动态势法研究得到的HCP及其同位素分子DCP在高激发态下的性能特点.

近二十年来， 多核结构计算水平的提高以及从头算电子结构理论的发展使得在合理的时间范围内确定分子体系的势能面和准确可靠地研究星际分子的动力学特性等相关性质成为可能  ［18-21］. 2014年， Pfeiffer和Rauhut  ［22］开发了基于特定状态的多参考振动相互作用方法（VMRCI）. 该方法基于振动多组态自洽场方法（VMCSCF）实现， 能够实现对体系振动能量的准确计算， 尤其适用于具有强非谐共振和高振动态（倍频或组合频）的系统. 本文首次将VMRCI方法应用于HCP分子振动结构计算， 实现了对分子体系振动特性的精确把握. 在计算过程中， 采用CCSD/cc-pVQZ方法优化HCP分子几何结构， 基于势能面， 采用VSCF、VMCSCF、VMRCI三种方法层层优化计算， 得到HCP分子的基频、倍频和组合频等振动频率以及绝对能量， 绘制出精确的红外和拉曼光谱图， 为实验和理论工作提供参考.

2 理论计算

分子光谱和势能面是研究分子结构的重要工具  ［23-27］. 本文使用Molpro 2019从头算程序包， 以波恩-奥本海默近似下的Hartree-Fock（HF）方法作为计算起点， 在CCSD/cc-PVQZ水平确定HCP分子的平衡几何结构和势能面扩展所需的简谐频率和简正坐标， 计算出的高水平势能面为之后的振动计算提供可靠的输入.

本文采用离散变量表示法求解振动薛定谔方程， 采用VSCF方法计算非谐振动频率. 考虑到HCP是线性闭壳分子， VSCF的计算采用限制性闭壳HF方法， 薛定谔方程中的沃森哈密顿量展开为

H ︿ = 1 2 ∑  αβ    J ︿    α-  π ︿    α  μ   αβ   J ︿    β-  π ︿    β -

1 8 ∑  α  μ   αα- 1 2 ∑  i      2   q  2   i +V  q  1，…， q   3N-6   （1）

其中，  V  q  i  表示势能面， 它由Molpro程序中的SURF模块提供， 在简正坐标下展开. 势能面函数给出的是分子能量随几何形状变化的函数. 本文研究的所有势能面都基于多模展开， 在Molpro从头算程序包中以完全自动化和并行化的方式获得， 并行化算法允许进行多级计算， 从而加快了计算速度. 多级计算指的是在势能的多模展开式中， 不同水平的电子结构方法被用于不同的项， 以获得该项在势能展开中的贡献. 势能面的多级展开具体见式（2）.

由于方程展开中的所有项都是指相对能量而不是总能量， 不影响薛定谔方程的求解结果， 因此这种近似是可行的.

V  q  1，…， q   3N-6 =  ∑  i  V   i  q   i     1D：CCSD  +  ∑  i

本文基于势能面计算了HCP分子的振动结构， 使用VSCF模态作为初始猜测， 再用VMCSCF方法对初始波函数进行优化得到多组态波函数， 考虑相关效应最终展开VMRCI计算. 量化计算基于HF方法展开， 由HF方法产生的误差称为相关能误差. VMCSCF方法通过在波函数中考虑激发组态校正相关能. 考虑的激发组态越多， 校正越大， 计算越精确， 但同时计算时间也会增加. 在VMCSCF计算中， 为了结果的收敛， 激发水平需比势能展开高一阶数， 势能在三阶后被截断，  因此VMCSCF计算多达四重激发  ［22］. 一般量子化学方法只考虑一重和二重这两种主要激发， 其中双重激发最为重要， 对相关能的校正起很大作用， 单激发对能量的校正远远小于双激发， 而VMCSCF计算多达四重激发， 计算精度更高，  VMRCI计算同样多达四重激发， VMCSCF、VMRCI方法考虑了振动状态间的静态和动态相关效应， 最终实现对于体系振动能量的精确计算.

3 结果与讨论

3.1 平衡几何参数和势能面

HCP是闭壳线性分子， 属于 C   ∞   v 点群对称， 非阿贝尔点群， 因此选取其子群 C   2   v 方便展开计算，   C   2   v 点群的不可约表示为A  1、B  1、B  2和A  2， 分别对应对称性为Σ、Σ、П和Π.

量化计算前首先进行几何结构优化， HCP的几何结构参数包括C-H键长、C-P键长以及HCP键角， 最低能量状态下HCP保持线性状态， 键角始终呈180°.在CCSD/cc-pVQZ方法下， 优化后的C-H键长为1.0698 ， C-P键长为1.5367 ， 而CCSD/cc-pVDZ方法下优化后的C-H键长为1.0865 ， C-P键长为1.5615 . 相比于实验值  r （CH）=1.0670 和 r （CP）=1.5420 ［28］， 明显cc-pVQZ基组下的优化结果更加接近实验值. 同基组不同方法下， HF/cc-pVQZ优化后的键长为  r （CH）=1.0611 和 r （CP）=1.5095 ［28］， 对比实验值可以发现相差较大， 仍然是CCSD/cc-pVQZ方法下的优化结果最接近实验值， 其中， C-H键长与实验值差0.0028 ， C-P键长与实验值差0.0053 ， 键角也非常符合.

在CCSD/cc-pVQZ理论水平上优化过几何平衡结构后， 在同一理论水平上确定了势能面的扩展和法坐标. 线性分子的振动自由度为3 n -5（ n 为原子数）， 因此， 作为线性三原子分子， HCP具有四个振动自由度， 即振动模式的数量. 如图1所示， HCP分子有四种正则振动模式， 分别为H-C拉伸振动、C-P拉伸振动、面内弯曲振动和面外弯曲振动. 在法向坐标中， 基本假设每种振动模式是相互独立的.

图2显示了四种振动模式下的一维势能曲线图， 其中简正坐标 q   1、 q   2分别对应HCP分子面内和面外弯曲振动模式，   q   3、 q   4分别对应C-P键拉伸振动模式和H-C键拉伸振动模式. 从图2a中可以看到， HCP面内弯曲振动模式的势能曲线呈“U”型， 图象关于 q   1=0呈现出明显的轴对称特性， 其原因是HCP分子在面内弯曲振动模式下， 简正坐标互为相反数的两个振动状态经过180°面外旋转后重合， 对应的分子形状本质上是相同的， 势能相等， 故而势能曲线对称.  另外， 图2a和图2b的势能曲线图是一样的， 这是因为HCP分子的面内弯曲振动模式和面外弯曲振动模式在90°旋转后实际上是相同的. 而在C-P和H-C拉伸振动模式下， 简正坐标互为相反数时对应着不同的分子形状， 势能不同， 因此势能曲线不再对称， 如图2c和图2d所示. 特别是 q   4坐标下势能曲线明显不对称， 且曲线随坐标变化更加陡峭， 这表明H-C拉伸振动模式在四种振动模式中对HCP分子体系振动势能影响最大. 在分子的相关计算中， 适当利用相应的对称性可以节省大量的计算时间和计算资源.

二维势能面描述了HCP分子在不同的坐标组合下势能变化的图象， 其中二维相对势能面如图3所示， 由于弯曲振动在简正坐标互为相反数时对应的分子形状相同， 因而图中含有 q   1或 q   2坐标的势能面均呈现出对称性. 如图3a所示， 势能面沿 q   1和 q   2坐标延伸， 由于面内、面外弯曲振动模式旋转后重合， 因此图象关于 q   1 = 0和 q   2 = 0同时呈现对称性， 并且在这两条线上势能为零， 分子处于平衡状态， 随着坐标变化势能面均匀向外延伸， 最终势能面四角出现四个势垒， 中心呈洼地状. 类似地， 可以预测以 q   1、 q   3为坐标和以 q   2、 q   3为坐标的势能面也应该是相同的. 结果确实如此， 并且图象分别关于 q   1   = 0和 q   2   = 0对称， 势能面整体稳定， 分子基本处于平衡状态. 图3d与图3e也是相同的， 相较其他坐标下的二维相对势能面更加平稳， 只有当 q   1、 q   2坐标值很大时， 势能才逐渐降低.  q   3、 q   4坐标下的势能面形状较其他势能面相比存在较大差异， 不再呈现对称性， 而是在两坐标值增大、减小时呈现出不同的势能变化规律. 这是因为C-P键拉伸振动模式和C-H键拉伸振动模式振动频率较大， 对分子体系的平衡状态有较大影响.

图4a展示了HCP分子在面内弯曲振动模态（ q   1）与面外弯曲振动模态（ q   2）耦合下的二维完整势能面.  q   1、 q   3坐标下和 q   2、 q   3坐标下的完整势能面图象与之相同， 图上无明显相变点， 整体呈“凹”状， 四角处呈现高势能， 中心处呈低洼状， 在零坐标处出现全局极小值点， 势能对坐标一阶导数为0， 对应着体系的平衡状态， 由极小值点向各个方向势能逐渐升高. 当 q   4作为坐标时， 完整势能面图象发生变化， 如图4b所示， 势能面上出现了一个鞍点. 这是由于H-C伸缩振动模态对分子系统的影响较大， 它通过影响HCP弯曲振动模式与C-P拉伸振动模式之间的耦合而影响体系势能.

3.2 振动能级和振动光谱

在量化计算之前， 首先确定分子振动自由度， 即振动基频的数目. 振动能级一般可表示为

E  vib=∑  i=1   v   i+  d   i 2  h v   i  （3）

其中，  d  i 表示模式 i 的简并度. HCP为线型分子， 具有四个振动自由度， 即四个基频， 其振动动力学标记为（ v   1，  v   2，  v   3，  v   4）. 由于HCP分子是线性分子， 振动能级间存在简并现象， 如表1中所示，  v   1和 v   2模式对应的振动频率始终相同， 存在简并现象. 这一现象也反映在随后列举的的倍频和组合频中.

由表1可以发现， C-P伸缩振动的基频（ v   3 = 1239.79 cm  -1）与HCP弯曲振动的倍频（2 v   2 = 1398.12 cm  -1）频率接近， 产生了2∶1的费米共振， 从而使得对应频率附近的吸收强度大大增强. 因此， 采用VMCSCF和VMRCI方法来计算具有共振现象的HCP分子体系是非常合理的. 通过对比表1中的数据可以发现， 相比于振动组态相互作用方法（VCI）， VMRCI方法计算出的振动频率更接近实验值， 计算精度也更高. 此外， 在四个基频中，  v   4模式对应的H-C拉伸振动频率最大， 对HCP分子的振动状态影响也最大.

VMRCI/cc-pVQZ方法下的HCP分子振动基频、倍频和组合频列于表2中， 同时列出的还有红外和拉曼振动强度. 可以发现， HCP分子红外活性和拉曼活性遵循互斥法则， 另外， 较高的振动强度都出现在基频处， 说明基频跃迁的能量更大. 同时， 本文的计算较精确， 即使一些强度较弱的振动也没有被忽视. 基于强度和频率数据， 可以进一步拟合绘制出HCP分子的光谱图象， 红外光谱图和拉曼光谱图分别如图5和图6所示.

根据表2中的数据， 振动强度最高的频带对应的频率为684.34 cm  -1， 强度值为65.93 km/mol， 对应HCP弯曲振动的基本跃迁. 然而， 在 图5红外光谱图中最强吸收峰处的峰值强度却为131.86 km/mol. 这是由于同一峰位处的两个峰强度发生了叠加现象. 这是光谱图的一个特点. 另一较强的谱带出现在3242.61 cm  -1频率处， 对应C-P拉伸振动的基本跃迁， 对应红外光谱图上次高的峰. 这两个峰是HCP分子红外光谱图的两个主要特征峰. 此外， 一些强度较小的峰位在光谱图上是不明显的， 如3880.41 cm  -1处和4620.81 cm  -1处的两个峰， 其强度分别为2.62 km/mol和1.54 km/mol， 但放大后清晰可见这些峰的存在. 这也可以与计算出的红外强度值进行比对， 从而证明本文计算出的光谱数据是丰富和细致的.

从图6可以看出， 拉曼光谱也有两个主要特征峰， 分别位于3242.61 cm   -1和1239.79 cm  -1处， 对一些弱峰进行放大可以清晰可见其峰位和峰强. 对比图5和图6可以发现， HCP分子的红外和拉曼活性不在同一振动频率下出现， 遵循互斥法则. 这两者计算出来之后相互补充， 共同构成HCP分子的详细光谱. 与红外光谱不同的是， 拉曼光谱信号较弱， 一般在实验中更难测量到. 但同时它又对分子振动分析有着重要意义， 因此本文关于拉曼光谱的计算为HCP分子光谱研究提供了参考数据.

4 结 论

本文基于Molpro程序包， 采用从头算方法对HCP分子势能面和振动光谱进行了研究. 首先， 采用cc-pVQZ/CCSD方法对HCP分子几何结构进行优化. 优化后的C-H键长和C-P键长与实验值非常接近. 计算得到了HCP分子丰富的一维势能曲线、二维相对势能面和二维完整势能面. 其中， H-C伸缩振动模式通过影响HCP弯曲振动模式和CP拉伸振动模式之间的耦合来影响分子体系的势能， 导致完整势能面中出现了鞍点. 基于势能面， 采用VMCSCF、VMRCI从头算方法对具有费米共振现象的CP分子振动能级和振动光谱展开计算. 发现了振动能级中的简并现象以及CP伸缩振动模式与HCP弯曲振动模式之间的费米共振现象. 该方法计算出的振动频率与实验值相吻合， 且更加接近实验值. 进一步绘制出的光谱图象为HCP分子光谱研究提供了参考数据. 本文整体工作对于HCP等含磷分子的实验和后续理论工作具有参考作用.

参考文献：

［1］   Gier T E. HCP， a unique phosphorus compound ［J］. J Am Chem Soc， 1961， 83： 1769.

［2］  Garneau J M， Cabana A. The ν  1 and ν  1 + ν  2 - ν  2 infrared absorption bands of H  12CP ［J］. J Mol Spectrosc， 1978， 69： 319.

［3］  Garneau J M， Cabona A. High-resolution infrared absorption spectrum of H  12CP. The ν  3 band ［J］. J Mol Spectrosc， 1980， 79： 502.

［4］  Garneau J M， Cabana A. The vibration-rotation spectrum of methinophosphide： l-type doubling， l-type resonance and the ν  2， 2ν  2， and 3ν  2 bands of H  12CP; the ν  1 and ν  2 bands of H  13CP ［J］. J Mol Spectrosc， 1981， 87： 490.

［5］  Cabana A， Doucet Y， Garneau J M，  et al.  The vibration-rotation spectrum of methinophosphide：  the overtone bands 2ν  1 and 2ν  3， the summation bands ν  1 +  ν  2 and ν  2 + ν  3， and the difference band ν  1 - ν  2   ［J］. J Mol Spectrosc， 1982， 96： 342.

［6］  Agúndez M， Cernicharo J， Guélin M. Discovery of phosphaethyne（HCP）in space： phosphorus chemistry in circumstellar envelopes ［J］. Astrophys J， 2007， 662： L91.

［7］  Bera P P， Yamaguchi Y， Schaefer III H F. Low-lying quartet electronic states of nitrogen dioxide ［J］. J Chem Phys， 2007， 127： 174303.

［8］  Yu Q， Bowman J. Communication： VSCF/VCI vibrational spectroscopy of H  7O  3  + and H  9O  4  + using high-level， many-body potential energy surface and dipole moment surfaces ［J］. J Chem Phys， 2017， 146： 121102.

［9］  Meier P， Oschetzki D， Pfeiffer F，  et al.  Towards an automated and efficient calculation of resonating vibrational states based on state-averaged multiconfigurational approaches ［J］. J Chem Phys， 2015， 143： 244111.

［10］  Transue W J， Velian A， Nava M，  et al.  A molecular precursor to phosphaethyne and its application in synthesis of the aromatic 1， 2， 3， 4-phosphatriazolate anion ［J］. J Am Chem Soc， 2016， 138： 6731.

［11］ Chen Y T， Chong D P. Comparison of theoretical vibrational and rotational energies of the HCP molecule with experimental values ［J］. J Chem Phys， 1993， 99： 8870.

［12］ Koput J. The equilibrium structure and spectroscopic constants of HCP-an  ab initio  study ［J］. Chem Phys Lett， 1996， 263： 401.

［13］ Koput J， Carter S. The potential energy surface and vibrational-rotational energy levels of HCP ［J］. Spectrochim Acta A， 1997， 53： 1091.

［14］ Beck C， Schinke R， Koput J. Vibrational spectroscopy of phosphaethyne（HCP）. I. Potential energy surface， variational calculations， and comparison with experimental data ［J］. J Chem Phys， 2000， 112： 8446.

［15］ H ltzl T， Szieberth D， Nguyen M，  et al.  Formation of phosphaethyne dimers： a mechanistic study ［J］. Chem Eur J， 2006， 12： 8044.

［16］ Wang A， Sun L， Fang C，  et al.  Study of highly excited vibrational dynamics of HCP integrable system with dynamic potential methods ［J］. Chin Phys B， 2020， 29： 013101.

［17］ Wang A， Sun L， Fang C，  et al.  The study of dynamic potentials of highly excited vibrational states of DCP： from case analysis to comparative study with HCP ［J］. Int J Mol Sci， 2016， 17： 1280.

［18］ Loos P F， Scemama A， Boggio-Pasqua M，  et al.  Mountaineering strategy to excited states： highly accurate energies and benchmarks for exotic molecules and radicals ［J］. J Chem Theory Comput， 2020， 16： 3720.

［19］ 张学富， 吕兵， 宋晓书， 等. BeH分子基态振动能级与光谱常数的理论研究［J］. 四川大学学报： 自然科学版， 2018， 55： 533.

［20］ 万冲， 宋晓书， 吕兵， 等. CaH分子低激发态的势能曲线和光谱性质［J］. 四川大学学报： 自然科学版， 2020， 57： 315.

［21］ 雷良建， 万冲， 王兴炜， 等. NH自由基基态及低激发态从头计算研究［J］. 四川大学学报： 自然科学版， 2020， 57： 953.

［22］ Pfeiffer F， Rauhut G. Multi-reference vibration correlation methods ［J］. J Chem Phys， 2014， 140： 064110.

［23］ Albuquerque J V， Shirsat R N. Prelude to molecular dynamics-II： investigation of potential energy surfaces using gaussian charge models ［J］. Chem Sel， 2020， 5： 11052.

［24］ Barreto P R P， Cruz A C P S， Euclides H O，  et al.  Spherical harmonics representation of the potential energy surface for the H  2…H  2 van der Waals complex ［J］. J Mol Model， 2020， 26： 277.

［25］ Melchakova I， Nikolaeva K M， Kovaleva E A，  et al.  Potential energy surfaces of adsorption and migration of transition metal atoms on nanoporus materials： the case of nanoporus bigraphene and G-C  3N  4［J］. Appl Surf Sci， 2021， 540： 148223.

［26］ Zhang K， Yin L， Liu G. Physically inspired atom-centered symmetry functions for the construction of high dimensional neural network potential energy surfaces ［J］. Comput Mater Sci， 2021， 186： 110071.

［27］ Qian W， Lu B， Tan G，  et al.  Vibrational spectrum and photochemistry of phosphaketene HPCO ［J］. Phys Chem Chem Phys， 2021， 23： 19237.

［28］ NIST. Computational chemistry comparison and benchmark DataBase ［DB/OL］. ［2022-05-22］. https：//cccbdb.nist.gov/.