H2-H2各向异性作用势的ab-initio法计算

李名锐,冯　娜，周　刚，初　哲，赵　南，陈春林

(西北核技术研究所,西安　710024)



H2-H2各向异性作用势的ab-initio法计算

李名锐,冯娜，周刚，初哲，赵南，陈春林

(西北核技术研究所,西安710024)

摘要：选用aug-cc-pVTZ及aug-cc-pVQZ基组对不同H2-H2位形进行了量子化学计算，外推出了氢分子间各向异性ab-initio作用势。依据实验数据确定出氢分子相互作用的5-体模型，在全局最优算法基础上拟合出了包含电四极矩项、多体极化作用项及三体关联项等在内的氢分子间各向异性作用势。结果验证了拟合势与ab-initio势能曲线几乎重合，两者偏差落在(-1.35 K,1.29 K)范围内；得出的氢分子第二维里系数与实验结果一致。

关键词：氢分子；各向异性作用势；从头计算

原子与分子相互作用势是分子计算的前提与基础[1]。长期以来，人们对H2-H2分子间相互作用势做了大量分子束散射、维里系数与输运系数测量等实验和理论研究[2]，给出了形式多样的各向同性势能表达式，如SG势[3]等，然而这些研究仅考虑分子质心间相互作用的各向同性作用势，不大适用于描述分子发生剧烈振动与转动的高温稠密液氢系统，此时需从分子层面上构建出各向异性作用势。

1各向异性作用势构建

1.1分子取向设定

1.2ab-initio作用势计算

基组是构建高精度ab-initio作用势的关键因素之一。本文利用GAUSSIAN03程序包[7]在MPn(n=2,3,4)和CCSD(T)水平下计算了单个H2分子的振动频率、键长及离解能，以考核各基组(6-31G、6-311G、aug-cc-pVDZ、aug-cc-pVTZ、aug-cc-pVQZ等)的精度，如表1所列。在相同基组下，MP2水平得出的振动频率值大于MP4相应的值，而MP2的键长及离解能却小于MP4的相应结果。在CCSD(T)水平下，利用aug-cc-pVTZ和aug-cc-pVQZ结果，并根据式(1)外推出的aug-cc-pV34Z值与实验值最为吻合。

(1)

图1　氢分子H2-H2取向设定Fig.1　Configuration settings for the H2-H2 dimer

ParameterExperimentalvalueMethod6-31G6-311GpVDZpVTZpVQZpV34ZVibrationalfrequency/cm-14401.2MP24533.64458.14463.74517.74515.04513.0MP34459.74366.44406.04464.54463.24462.3MP44414.94317.54380.04435.64433.74432.3CCSD(T)4370.34270.64347.64404.04402.84401.9Bondlength/bohr1.4009MP21.39371.39311.42661.39351.39141.3899MP31.40181.40181.43371.39861.39651.3950MP41.40621.40631.43621.40101.39921.3979CCSD(T)1.40971.41031.43941.40411.40181.4001Dissociationenergy/Ha0.16467MP20.137340.136100.147380.155090.156550.15762MP30.142710.141440.153900.160740.161700.16240MP40.144430.143220.155610.162280.163240.16394CCSD(T)0.145270.144150.156320.162840.163940.16474

图2　氢分子ab-initio作用势Fig.2　Ab-initio potential curves of the H2-H2 dimer

1.3多体极化作用势

假设在静电场F中有N个偶极子极化点，那么m点的诱导偶极矩

(2)

(3)

式中，αm是m点的极化率张量;Tmn为偶极子场张量;I为3×3单位矩阵;rmn=|rmn|为点m与点n的间距，令xk(k=1,2,3)代表矢量rmn在笛卡儿坐标系中k轴的分量;rr′为张量xkxl。那么氢分子的极化率张量

(4)

其中, 矩阵A=[α-1+T]-1。由于氢分子极化率张量仅与电荷的坐标及其极化率参数有关，与分子环境无关，为避免位于同条直线且方向相同的两诱导偶极矩间的极化作用出现无穷大值，本文采取Thole提出的方法进行修正[9]，取

(5)

式中,λ1(r)=1-(b2r2/2+br+1)exp(-br)，λ2(r)=1-(b3r3/6+b2r2/2+br+1)exp(-br)，对氢原子而言,b=2.130 4。

(6)

(7)

图3是不同取向下两氢分子间多体极化作用势的计算结果。可以看出，4种取向的多体极化作用势均为负值，L取向的极化能绝对值最大，H取向与X取向的势能曲线几乎重合。尽管氢分子间多体极化作用势值远小于ab-initio势值，但Belof认为前者可提高稠密液氢系统状态方程的精度。

图3　氢分子间多体极化作用势Fig.3　The multi-body polarization potentialcurves for the H2-H2 dimer

1.4电四极矩相互作用势

电四极矩相互作用势为[10]

(8)

式中，rij为两氢分子的质心间距; Q为氢分子的电四极矩;cosγ=cosαcosβ+sinαsinβcosφ。

图4是不同取向下氢分子电四极矩相互作用势的计算结果。4种取向中仅T取向的电四极矩作用势为负值；L取向的电四极矩作用势绝对值最大，T取向次之，X取向最小；与多体极化作用势不同，H取向与X取向的电四极矩作用势差别明显。

图5是将ab-initio势减去电四极矩作用势和多体极化作用势的结果。相减后不同分子取向的作用势势阱深度间的差异已明显减小，说明电四极矩作用及多体极化作用是导致分子间作用势为各向异性的主要因素之一。考虑到复杂的ab-initio势能面上存在多个局部极小值点，直接拟合该势能面的精度可能较低，若利用ab-initio势减去电四极矩作用势和多体极化作用势，并将此相对平缓势能面作为拟合目标函数，会使拟合变得简单些。

图4　氢分子电四极矩相互作用势Fig.4　The electric quadrupole potential curvesfor the H2-H2 dimer

图5　相减后的作用势Fig.5　The fitting target potential curves of fourdistinct relative H2 orientations

1.5拟合函数

通常认为多中心势能函数比较适用于描述分子间的相互作用，但2-体或3-体分子模型还不能较好地重构出ab-initio势能面，需采用5-体模型来获得足够的拟合精度[11]。本文将在Van构造的5-体作用模型基础上，附加SG势函数中的球状Axilrod-Teller贡献C9项[3]，以考虑稠密液氢系统的三体关联作用。为降低拟合计算的复杂度，利用aug-cc-pV34Z势减去电四极矩作用势及多体极化作用势，将得到的较为平缓势能面作为拟合目标势能面。构建的拟合势能函数为

(9)

式中，f1(rij)=(1+exp(-2(δijrij-2)))-15，f2(rij)=1-exp(-βijrij)；依据氢分子电四极矩实验值设定qH=0.393 27e，qN=-2qH，qM=0；假定令βij=1bohr-1，不同位置点具有相同参数δij。由氢分子对称性可知，共需拟合37个参数。

2结果分析

由于拟合势能函数式(9)含有的待优化参数较多且复杂，若选用常用的最小二乘法进行拟合，很可能得到的是局部最优势能曲面。为此，本文将在通用全局优化算法基础上[12]，采用Simplex法[13]进行优化求解，各参数的拟合结果列于表2。本文获得了精度较高的全局最优势能面，其拟合相关系数R=0.999 9；拟合势能值与ab-initio计算值的偏差落在(-4.28 μHa,4.08 μHa)或(-1.35 K,1.29 K) 范围内, 如图6所示。

图7是不同拟合作用势与ab-initio势的对比结果，图中实心正方形为本文计算出的ab-initio作用势，空心圆为Belof[5]的3-体LJ势拟合结果，空心正三角为本文拟合结果。T取向时，Belof的3-体拟合势与ab-initio势存在着明显差异，且无论Belof如何增加LJ作用势的参数均不能改善该拟合结果[5]。相比而言，不同取向下本文拟合出的势能曲线与ab-initio势能曲线均几乎重合。

固定分子间距r，利用MC法分别在方位角α∈[0°,180°]，β∈[0°,180°]，φ∈[0°,180°]上随机抽取106个数值，代入拟合出的氢分子间各向异性势能函数进行统计平均，可以得到对应的各向同性作用势Uiso(r)。图8是Uiso(r)与SG作用势的对比结果。各向同性作用势Uiso(r)的势阱位置与SG作用势的势阱位置相同，但前者势阱较后者要深出4 K左右。整体上看，本文拟合出的各向同性作用势Uiso(r)较SG势略显偏硬。获得各向同性作用势后可依据式(10)计算出氢的第二维里系数：

表2　采用Simplex法的拟合结果

图7　拟合势与ab-initio势比较Fig.7　The fitting potential vs. ab-initio curves

(10)

图8　平均后的各向同性势Fig.8　Average results of the isotropic potential

图9　第二维里系数B2(T)Fig.9　Results of the second virial coefficient

3结论

本文从量子计算、相互作用势模型以及拟合算法等方面综合优化提高了各向异性作用势的精度。在CCSD(T)水平下选取aug-cc-pVTZ及aug-cc-pVQZ基组对不同H2-H2位形进行大量量子化学计算，利用外推法进一步提高了氢分子间各向异性ab-initio作用势精度。依据实验数据确定出氢分子间各向异性相互作用的5-体模型，并考虑氢分子间的三体关联作用项，利用存在多个局部极小值点的复杂ab-initio势减去电四极矩作用势和多体极化作用势，将相减后势阱变化相对平缓的势能曲面作为拟合目标函数，以降低拟合计算的复杂度，在此基础上采取全局最优算法进行拟合。结果表明：本文拟合势与ab-initio势能曲线几乎重合，两者偏差落在(-4.28 μHa,4.08 μHa)或(-1.35 K,1.29 K)范围内；与Van及Belof作用势相比，得出的氢分子第二维里系数结果与实验值更为吻合。

参考文献

[1]令狐荣锋, 徐梅, 吕兵, 等. He原子与O2分子相互作用势及碰撞微分截面的研究[J]. 原子与分子物理学报, 2013, 30(4): 591-596. (LINGHU Rong-feng, XU Mei, LYU Bing, et al. A theoretical study on interaction potential energy surface of He-O2[J]. Journal of Atomic and Molecular Physics, 2013, 30(4): 591-596.)

[2]MATSUISHI K, GREGORYANZ E, MAO H K, et al. Equation of state and intermolecular interactions in fluid hydrogen from Brillouin scattering at high pressures and temperatures[J]. J Chem Phys, 2003, 118(23): 10 683-10 695.

[3]SILVERA I F, GOLDMAN V V. The isotropic intermolecular potential for h2and d2in the solid and gas phase[J]. J Chem Phys, 1978, 69(9): 4 209-4 213.

[4]DIEP P, JOHNSON J K. An accurate H2-H2interaction potential from first principles[J]. J Chem Phys, 2000, 112(10): 4 465-4 473.

[5]BELOF J L, STERN A C, SPACE B. An accurate and transferable intermolecular diatomic hydrogen potential for condensed phase simulation[J]. J Chem Theory Comput, 2008, 4(8): 1 332-1 337.

[6]VAN T P. Ab-initio calculation of intermolecular potentials, prediction of second virial coefficients for dimers H2-H2, H2-O2, F2-F2and H2-F2, and Monte Carlo simulations of the vapor-liquid equilibria for hydrogen and fluorine[D]. Cologne: University of Cologne, 2006.

[7]FRISCH M J, TRUCKS G W, SCHLEHGEL H B, et al. GAUSSIAN 03, Revision B.02, Gaussian[M]. USA: Wallingford, 2003.

[8]BOYS S F, BERNARDI F. The calculations of small molecular interactions by differences of separate total energy: some procedures with reduced errors[J]. Mol Phys, 1970, 19(4): 553-556.

[9]THOLE B T. Molecular polarizabilities calculated with a modified dipole interaction[J]. Chem Phys, 1981, 59(3): 341-350.

[10]ALLEN M P, TILDESLEY D J. Computer Simulation of Liquids[M]. Oxford: Oxford University Press, 1987.

[11]BOCK S, BICH E, VOGEL E. A new intermolecular potential energy surface for carbon dioxide from ab initio calculations[J]. J Chem Phys, 2000, 257(2/3): 147-156.

[12]WILLIAM F C, RONCARATTI L F, GARGANO R, et al. Fitting potential energy surface of reactive system via genetic algorithm[J]. Int J Quantum Chem, 2006, 106(3): 2 650-2 657.

[13]YEN J, LIAO J C, LEE BOGIU . A hybrid approach to modeling metabolic systems using a genetic algorithm and simplex method[J]. IEEE Trans Syst Man Cybem B, 1998, 28(2): 173-191.

[14]LIDE D R. Handbook of Chemistry and Physics[M]. 84th ed. Boca Raton: CRC Press, 2004.

Ab-Initio Calculations for the Anisotropic Interaction Potential of H2-H2Dimer

LI Ming-rui,FENG Na,ZHOU Gang,CHU Zhe,ZHAO Nan,CHEN Chun-lin

(Northwest Institute of Nuclear Technology,Xi’an710024,China)

Abstract:The new ab-initio anisotropic intermolecular potential of H2-H2 system is obtained by using the extrapolation method and the data from a large number of quantum chemical calculations with the basis sets of aug-cc-pVTZ and aug-cc-pVQZ.Based on the ab-initio simulation results, an analytical 5-site potential function comprising the electric quadrupole interaction, multi-body polarization interaction and three-body correlation effects has been constructed. The adjustable parameters of the potential function are determined by the global optimization method and some experimental data. The results show that the fitting potential curves are almost coincident with the ab-initio potential curves mentioned above, with a narrow residual range from -1.35 K to 1.29 K; and the second virial coefficients of hydrogen dimers calculated by the fitting potential are also in good agreement with those of experiments.

Key words:hydrogen molecule;anisotropic potential;ab-initio calculation

中图分类号：O562;TB33

文献标志码:A

文章编号：2095-6223(2016)010901(7)

作者简介：李名锐(1983- )，男，湖北黄石人，助理研究员，博士，主要从事冲击波物理研究。E-mail:limingrui@nint.ac.cn

收稿日期：2015-05-04；修回日期：2015-10-27