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

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

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

?

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

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

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

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

關(guān)鍵詞：氫分子；各向異性作用勢；從頭計算

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

1各向異性作用勢構(gòu)建

1.1分子取向設(shè)定

1.2ab-initio作用勢計算

基組是構(gòu)建高精度ab-initio作用勢的關(guān)鍵因素之一。本文利用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相應(yīng)的值，而MP2的鍵長及離解能卻小于MP4的相應(yīng)結(jié)果。在CCSD(T)水平下，利用aug-cc-pVTZ和aug-cc-pVQZ結(jié)果，并根據(jù)式(1)外推出的aug-cc-pV34Z值與實驗值最為吻合。

(1)

圖1　氫分子H2-H2取向設(shè)定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多體極化作用勢

假設(shè)在靜電場F中有N個偶極子極化點，那么m點的誘導(dǎo)偶極矩

(2)

(3)

式中，αm是m點的極化率張量;Tmn為偶極子場張量;I為3×3單位矩陣;rmn=|rmn|為點m與點n的間距，令xk(k=1,2,3)代表矢量rmn在笛卡兒坐標(biāo)系中k軸的分量;rr′為張量xkxl。那么氫分子的極化率張量

(4)

其中, 矩陣A=[α-1+T]-1。由于氫分子極化率張量僅與電荷的坐標(biāo)及其極化率參數(shù)有關(guān)，與分子環(huán)境無關(guān)，為避免位于同條直線且方向相同的兩誘導(dǎo)偶極矩間的極化作用出現(xiàn)無窮大值，本文采取Thole提出的方法進(jìn)行修正[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是不同取向下兩氫分子間多體極化作用勢的計算結(jié)果?？梢钥闯觯?種取向的多體極化作用勢均為負(fù)值，L取向的極化能絕對值最大，H取向與X取向的勢能曲線幾乎重合。盡管氫分子間多體極化作用勢值遠(yuǎn)小于ab-initio勢值，但Belof認(rèn)為前者可提高稠密液氫系統(tǒng)狀態(tài)方程的精度。

圖3　氫分子間多體極化作用勢Fig.3　The multi-body polarization potentialcurves for the H2-H2 dimer

1.4電四極矩相互作用勢

電四極矩相互作用勢為[10]

(8)

式中，rij為兩氫分子的質(zhì)心間距; Q為氫分子的電四極矩;cosγ=cosαcosβ+sinαsinβcosφ。

圖4是不同取向下氫分子電四極矩相互作用勢的計算結(jié)果。4種取向中僅T取向的電四極矩作用勢為負(fù)值；L取向的電四極矩作用勢絕對值最大，T取向次之，X取向最??；與多體極化作用勢不同，H取向與X取向的電四極矩作用勢差別明顯。

圖5是將ab-initio勢減去電四極矩作用勢和多體極化作用勢的結(jié)果。相減后不同分子取向的作用勢勢阱深度間的差異已明顯減小，說明電四極矩作用及多體極化作用是導(dǎo)致分子間作用勢為各向異性的主要因素之一?？紤]到復(fù)雜的ab-initio勢能面上存在多個局部極小值點，直接擬合該勢能面的精度可能較低，若利用ab-initio勢減去電四極矩作用勢和多體極化作用勢，并將此相對平緩勢能面作為擬合目標(biāo)函數(shù)，會使擬合變得簡單些。

圖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擬合函數(shù)

通常認(rèn)為多中心勢能函數(shù)比較適用于描述分子間的相互作用，但2-體或3-體分子模型還不能較好地重構(gòu)出ab-initio勢能面，需采用5-體模型來獲得足夠的擬合精度[11]。本文將在Van構(gòu)造的5-體作用模型基礎(chǔ)上，附加SG勢函數(shù)中的球狀A(yù)xilrod-Teller貢獻(xiàn)C9項[3]，以考慮稠密液氫系統(tǒng)的三體關(guān)聯(lián)作用。為降低擬合計算的復(fù)雜度，利用aug-cc-pV34Z勢減去電四極矩作用勢及多體極化作用勢，將得到的較為平緩勢能面作為擬合目標(biāo)勢能面。構(gòu)建的擬合勢能函數(shù)為

(9)

式中，f1(rij)=(1+exp(-2(δijrij-2)))-15，f2(rij)=1-exp(-βijrij)；依據(jù)氫分子電四極矩實驗值設(shè)定qH=0.393 27e，qN=-2qH，qM=0；假定令βij=1bohr-1，不同位置點具有相同參數(shù)δij。由氫分子對稱性可知，共需擬合37個參數(shù)。

2結(jié)果分析

由于擬合勢能函數(shù)式(9)含有的待優(yōu)化參數(shù)較多且復(fù)雜，若選用常用的最小二乘法進(jìn)行擬合，很可能得到的是局部最優(yōu)勢能曲面。為此，本文將在通用全局優(yōu)化算法基礎(chǔ)上[12]，采用Simplex法[13]進(jìn)行優(yōu)化求解，各參數(shù)的擬合結(jié)果列于表2。本文獲得了精度較高的全局最優(yōu)勢能面，其擬合相關(guān)系數(shù)R=0.999 9；擬合勢能值與ab-initio計算值的偏差落在(-4.28 μHa,4.08 μHa)或(-1.35 K,1.29 K) 范圍內(nèi), 如圖6所示。

圖7是不同擬合作用勢與ab-initio勢的對比結(jié)果，圖中實心正方形為本文計算出的ab-initio作用勢，空心圓為Belof[5]的3-體LJ勢擬合結(jié)果，空心正三角為本文擬合結(jié)果。T取向時，Belof的3-體擬合勢與ab-initio勢存在著明顯差異，且無論Belof如何增加LJ作用勢的參數(shù)均不能改善該擬合結(jié)果[5]。相比而言，不同取向下本文擬合出的勢能曲線與ab-initio勢能曲線均幾乎重合。

固定分子間距r，利用MC法分別在方位角α∈[0°,180°]，β∈[0°,180°]，φ∈[0°,180°]上隨機抽取106個數(shù)值，代入擬合出的氫分子間各向異性勢能函數(shù)進(jìn)行統(tǒng)計平均，可以得到對應(yīng)的各向同性作用勢Uiso(r)。圖8是Uiso(r)與SG作用勢的對比結(jié)果。各向同性作用勢Uiso(r)的勢阱位置與SG作用勢的勢阱位置相同，但前者勢阱較后者要深出4 K左右。整體上看，本文擬合出的各向同性作用勢Uiso(r)較SG勢略顯偏硬。獲得各向同性作用勢后可依據(jù)式(10)計算出氫的第二維里系數(shù)：

表2　采用Simplex法的擬合結(jié)果

圖7　擬合勢與ab-initio勢比較Fig.7　The fitting potential vs. ab-initio curves

(10)

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

圖9　第二維里系數(shù)B2(T)Fig.9　Results of the second virial coefficient

3結(jié)論

本文從量子計算、相互作用勢模型以及擬合算法等方面綜合優(yōu)化提高了各向異性作用勢的精度。在CCSD(T)水平下選取aug-cc-pVTZ及aug-cc-pVQZ基組對不同H2-H2位形進(jìn)行大量量子化學(xué)計算，利用外推法進(jìn)一步提高了氫分子間各向異性ab-initio作用勢精度。依據(jù)實驗數(shù)據(jù)確定出氫分子間各向異性相互作用的5-體模型，并考慮氫分子間的三體關(guān)聯(lián)作用項，利用存在多個局部極小值點的復(fù)雜ab-initio勢減去電四極矩作用勢和多體極化作用勢，將相減后勢阱變化相對平緩的勢能曲面作為擬合目標(biāo)函數(shù)，以降低擬合計算的復(fù)雜度，在此基礎(chǔ)上采取全局最優(yōu)算法進(jìn)行擬合。結(jié)果表明：本文擬合勢與ab-initio勢能曲線幾乎重合，兩者偏差落在(-4.28 μHa,4.08 μHa)或(-1.35 K,1.29 K)范圍內(nèi)；與Van及Belof作用勢相比，得出的氫分子第二維里系數(shù)結(jié)果與實驗值更為吻合。

參考文獻(xiàn)

[1]令狐榮鋒, 徐梅, 呂兵, 等. He原子與O2分子相互作用勢及碰撞微分截面的研究[J]. 原子與分子物理學(xué)報, 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

文獻(xiàn)標(biāo)志碼:A

文章編號：2095-6223(2016)010901(7)

作者簡介：李名銳(1983- )，男，湖北黃石人，助理研究員，博士，主要從事沖擊波物理研究。E-mail:limingrui@nint.ac.cn

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