李名銳,馮　娜，周　剛，初　哲，趙　南，陳春林

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

?

H2-H2各向異性作用勢的ab-initio法計算

李名銳,馮娜，周剛，初哲，趙南，陳春林

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

摘要：選用aug-cc-pVTZ及aug-cc-pVQZ基組對不同H2-H2位形進行了量子化學計算，外推出了氫分子間各向異性ab-initio作用勢。依據實驗數據確定出氫分子相互作用的5-體模型，在全局最優(yōu)算法基礎上擬合出了包含電四極矩項、多體極化作用項及三體關聯(lián)項等在內的氫分子間各向異性作用勢。結果驗證了擬合勢與ab-initio勢能曲線幾乎重合，兩者偏差落在(-1.35 K,1.29 K)范圍內；得出的氫分子第二維里系數與實驗結果一致。

關鍵詞：氫分子；各向異性作用勢；從頭計算

原子與分子相互作用勢是分子計算的前提與基礎[1]。長期以來，人們對H2-H2分子間相互作用勢做了大量分子束散射、維里系數與輸運系數測量等實驗和理論研究[2]，給出了形式多樣的各向同性勢能表達式，如SG勢[3]等，然而這些研究僅考慮分子質心間相互作用的各向同性作用勢，不大適用于描述分子發(fā)生劇烈振動與轉動的高溫稠密液氫系統(tǒng)，此時需從分子層面上構建出各向異性作用勢。

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。由于氫分子極化率張量僅與電荷的坐標及其極化率參數有關，與分子環(huán)境無關，為避免位于同條直線且方向相同的兩誘導偶極矩間的極化作用出現無窮大值，本文采取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認為前者可提高稠密液氫系統(tǒng)狀態(tài)方程的精度。

圖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取向最?。慌c多體極化作用勢不同，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]，以考慮稠密液氫系統(tǒng)的三體關聯(lián)作用。為降低擬合計算的復雜度，利用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)含有的待優(yōu)化參數較多且復雜，若選用常用的最小二乘法進行擬合，很可能得到的是局部最優(yōu)勢能曲面。為此，本文將在通用全局優(yōu)化算法基礎上[12]，采用Simplex法[13]進行優(yōu)化求解，各參數的擬合結果列于表2。本文獲得了精度較高的全局最優(yōu)勢能面，其擬合相關系數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個數值，代入擬合出的氫分子間各向異性勢能函數進行統(tǒng)計平均，可以得到對應的各向同性作用勢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結論

本文從量子計算、相互作用勢模型以及擬合算法等方面綜合優(yōu)化提高了各向異性作用勢的精度。在CCSD(T)水平下選取aug-cc-pVTZ及aug-cc-pVQZ基組對不同H2-H2位形進行大量量子化學計算，利用外推法進一步提高了氫分子間各向異性ab-initio作用勢精度。依據實驗數據確定出氫分子間各向異性相互作用的5-體模型，并考慮氫分子間的三體關聯(lián)作用項，利用存在多個局部極小值點的復雜ab-initio勢減去電四極矩作用勢和多體極化作用勢，將相減后勢阱變化相對平緩的勢能曲面作為擬合目標函數，以降低擬合計算的復雜度，在此基礎上采取全局最優(yōu)算法進行擬合。結果表明：本文擬合勢與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