鄭興榮， 付 云， 付文羽， 李向富， 李繼弘

(1.隴東學(xué)院物理系， 慶陽 745000； 2. 貴州民族大學(xué)理學(xué)院， 貴陽 550025)

基于原子晶體構(gòu)型的高壓固氬中的多體相互作用

鄭興榮1， 付 云2， 付文羽1， 李向富1， 李繼弘1

(1.隴東學(xué)院物理系， 慶陽 745000； 2. 貴州民族大學(xué)理學(xué)院， 貴陽 550025)

X射線衍射實(shí)驗(yàn)顯示固氬是面心立方(fcc)晶格結(jié)構(gòu)，目前對晶體氬的研究只限于兩體，三體以及四體相互作用勢. 本文利用多體展開方法和超分子單、雙 (三)重激發(fā)耦合簇理論(CCSD(T))對固氬fcc晶格結(jié)構(gòu)的三體和四體的幾何構(gòu)型、幾何參數(shù)、不同體積下所有三體和四體構(gòu)型的勢能以及各構(gòu)型所占比例等幾個方面進(jìn)行了準(zhǔn)確的量子化學(xué)計算. 結(jié)果表明:所有三體構(gòu)型中對總的三體勢能貢獻(xiàn)最大的是構(gòu)型1、構(gòu)型6、構(gòu)型12和構(gòu)型23；三體勢及其交換部分和色散部分的計算結(jié)果與現(xiàn)有解析經(jīng)驗(yàn)勢在長程部分符合得非常好，但在短程部分有較小差異. 所有的四體構(gòu)形中對總的四體勢能貢獻(xiàn)最大的是構(gòu)型1，構(gòu)型2，構(gòu)型4，構(gòu)型5，構(gòu)型7和構(gòu)型8；四體勢及其交換勢部分和色散部分的計算結(jié)果尚無解析經(jīng)驗(yàn)勢可比較.利用這些特殊構(gòu)型的相關(guān)數(shù)據(jù)并結(jié)合其它構(gòu)型，可擬合出更準(zhǔn)確的三體經(jīng)驗(yàn)勢函數(shù)及其參數(shù)，也為擬合四體經(jīng)驗(yàn)勢函數(shù)及其參數(shù)提供了重要的參考價值.

固氬； fcc晶格結(jié)構(gòu)； CCSD(T)； 三體勢； 四體勢

1 引 言

惰性元素具有球?qū)ΨQ的滿殼層電子結(jié)構(gòu)，一般情況下很難與其他物質(zhì)發(fā)生相互作用，但是在太陽、地球和其他行星內(nèi)部的惰性元素均以高壓狀態(tài)存在[1,2]，且能夠發(fā)生化學(xué)反應(yīng). 氬在實(shí)驗(yàn)上常用作靜態(tài)高壓實(shí)驗(yàn)的傳壓介質(zhì)[3]，在理論上是用量子力學(xué)方法研究高壓原子晶體結(jié)構(gòu)的理想體系[4-6]. 人們利用金剛石對頂砧技術(shù)(DAC)測得了氬的物態(tài)方程，在室溫下氬的最大壓強(qiáng)已達(dá)到114 GPa，體積壓縮比為3.06[7,8]. 常溫低壓下實(shí)驗(yàn)測得固氬是fcc晶格結(jié)構(gòu)[7]. Aziz等人利用Hartree-Fock方法提出了Hartree-Fock Dispersion(HFD)-B[9,10], HFDID[11]以及Stawomir[12]等經(jīng)驗(yàn)勢，這些勢函數(shù)能較好的描述兩體相互作用,但當(dāng)原子間距較小時與實(shí)驗(yàn)值有一定的差異. 隨后Petr等人提出的兩體勢函數(shù)[13]能較準(zhǔn)確地描述固氬的兩體相互作用, 也能模擬出氬的多組實(shí)驗(yàn)數(shù)據(jù)及固氬的零溫零壓性質(zhì). Loubeyre[14,15]和Freiman[16]等人提出了氬的三體經(jīng)驗(yàn)勢函數(shù)，壓強(qiáng)高于某一值時計算出的面心立方(fcc)氬的壓縮線和結(jié)合能對應(yīng)的理論壓強(qiáng)值低于實(shí)驗(yàn)值；而當(dāng)晶格體積小于某一值的高壓區(qū)，由于三體勢貢獻(xiàn)過多的吸引效應(yīng)，出現(xiàn)了氬的結(jié)合能隨體積的減小而下降的不合理趨勢. 因此，我們一方面有必要通過量子化學(xué)計算來檢測該三體經(jīng)驗(yàn)勢的準(zhǔn)確性，另一方面有必要探索更高階多體相互作用在固氬中的貢獻(xiàn). 本文采用多體展開方法和超分子單、雙 (三)重激發(fā)耦合簇理論(CCSD(T))方法和augmented correlation consistent quadruple zeta(aug-cc-pVQZ)基矢[17]計算了高壓下固氬的fcc晶體中不同最近鄰距離下任意三個原子間的三體構(gòu)型和任意四個原子間的四體構(gòu)型，對三體和四體的幾何構(gòu)型、幾何參數(shù)、不同體積下所有三體和四體構(gòu)型的勢能以及各構(gòu)型所占比例等幾個方面對固氬的多體構(gòu)型作了系統(tǒng)的分析. 由此得出不同最近鄰距離下不同多體勢的貢獻(xiàn)程度，對于驗(yàn)證現(xiàn)有三體勢函數(shù)的解析式以及參數(shù)的準(zhǔn)確性和擬合四體勢函數(shù)表達(dá)式可能具有比較重要的理論價值.

2 理論模型

固氬晶體中，任一原子被其不同鄰近原子包圍. 由于屏蔽效應(yīng)，任一原子o對結(jié)合排斥勢的貢獻(xiàn)主要來自于o與周圍近鄰殼層(i=1,2,3,…)的有限的(n-1)個近鄰原子的相互作用. 原子o和這(n-1)個原子構(gòu)成一n原子的團(tuán)簇Arn，該團(tuán)簇的勢能等于這n個原子從無窮遠(yuǎn)處移到當(dāng)前位置處團(tuán)簇總能量的變化

Vn(0,1,2,3,…,n-1)=E(r0,r1,r2,…,rn-1)-nE0,

(1)

其中E0代表孤立原子的基態(tài)能量，E(r0,r1,r2,…,rn-1) 代表n個原子處于r0,r1,r2,…,rn-1位置時團(tuán)簇的基態(tài)總能量，它可以通過求解體系的薛定諤方程得到

φ=Eφ.

(2)

對于處在位置坐標(biāo)r0和ri下的兩氬原子o和i，它們所構(gòu)成團(tuán)簇的相互作用能，可用一個兩體勢描述

u2(o,i)=E(ro,ri)-2E0,

(3)

其中u2(o,i)代表兩氬原子o，i間的兩體相互作用，E(r0,ri)是兩個氬原子分別位于r0，r0位置時它們所構(gòu)成原子團(tuán)簇Ar2的基態(tài)能量，E0代表一個孤立氬原子的基態(tài)能量. 當(dāng)?shù)谌齻€氬原子j靠近這兩個氬原子時，它們之間的的電子云分布彼此會發(fā)生變化，從而使得三個氬原子間存在三體相互作用，即三體勢

u3(o,i,j)=E(ro,ri,rj)-3E0

-u2(o,i)-u2(o,j)-u2(i,j),

(4)

上式中E(r0,ri,rj)是三個氬原子分別位于r0,ri,rj位置時三個氬原子所構(gòu)成團(tuán)簇Ar3的基態(tài)能量. 同理，當(dāng)?shù)谒膫€原子k接近這三個相互作用的原子時，三個原子之間的電子云分布會發(fā)生變化，從而使得四個氬原子間存在四體相互作用，即四體勢

u4(o,i,j,k)=E(ro,ri,rj,rk)-4E0-[u3(o,i,j)+u3(o,i,k)+u3(o,j,k)+u3(i,j,k)]-[u2(o,i)+u2(o,j)+u2(o,k)+u2(i,j)+u2(i,k)+u2(j,k)]

(5)

同理，我們可以得到五體勢，六體勢……

在高壓條件下，晶格中原子之間的距離減少.任一中心原子o與其近鄰原子間的相互作用Vn可展開為原子o與(n-1)個近鄰原子的兩體式、三體勢、四體勢、五體勢……之和，即

=U2(O)+U3(O)+U4(O)+…

(6)

其中U2(O)，U3(O)，U4(O)分別為中心原子o與近鄰的(n-1)個原子中任意一個、兩個、三個原子的兩體，三體，四體相互作用，且U2(O)，U3(O)，U4(O)分別為各項(xiàng)對所有的近鄰原子求和. 在實(shí)際求解過程中我們運(yùn)用Gamess計算程序來求出氬原子團(tuán)簇的多體勢能，并運(yùn)用CCSD(T)的從頭計算方法和多體展開理論，對fcc固氬的性質(zhì)作了研究.

3 結(jié)果與討論

3.1 三原子構(gòu)型的相互作用勢能

3.1.1 Ar3三原子團(tuán)簇的構(gòu)型

圖1 fcc晶格中Ar3團(tuán)簇的幾何圖形Fig. 1 Geometries of Ar3 clusters in fcc lattice

表1 fcc晶格中Ar3團(tuán)簇的幾何參數(shù)，θ1，θ2，θ3是三角形的三個內(nèi)角(從小到大的順序排列)，NP是三邊長求和后的三分之一，NG表示權(quán)重因子

Table 1 Geometry parameters of Ar3clusters in fcc lattice,θ1，θ2，θ3are the three angles (in degree and the value from smallest to largest) of the triangle,NPis a third of the sum of three size lengths,NGrepresents the weight factor

Geometryθ1θ2θ3NPNG形狀126060601.0001.7322424等邊三角形34554.749.848.254.749.865.970.580.465.91.8211.9001.626124872銳角等腰三角形674545901.1381.6093612等腰直角三角形840.240.299.62.03748鈍角等腰三角形935.354.7901.38272直角三角形1035.335.3109.52.09812鈍角等腰三角形1133.673.273.21.48872銳角等腰三角形121330301201.2442.1557224鈍角等腰三角形143060901.57748直角三角形151629.225.436.747.9114.1106.81.9311.6564848鈍角三角形1724.124.1131.82.20924鈍角等腰三角形181919.518.435.326.6125.31351.7271.5502424鈍角三角形2016.816.8146.42.26024鈍角等腰三角形212215.810.919.519.1144.71502.0491.7932448鈍角三角形232425001801.3331.8862.3096312平角三角形

3.1.2 不同體積下所有三體構(gòu)型的能量

固氬的三體勢構(gòu)型共有25種，我們比較了這25種三原子構(gòu)型在不同體積下三體勢能所占的比例，如圖2所示. 任何一組三原子構(gòu)型的能量不僅與三體勢φ3有關(guān)，還決定于相應(yīng)的權(quán)重因子NG(見表1). 圖2給出了體積分別為4.534 cm3/mol，10.385 cm3/mol，19.867 cm3/mol，27.253 cm3/mol時四種體積下所有三體構(gòu)型下的能量相對大小，由圖可以看出，在所有的構(gòu)型中對總的三體勢能貢獻(xiàn)最大的是構(gòu)型1、構(gòu)型6、構(gòu)型12和構(gòu)型23. 構(gòu)型1是一個權(quán)重因子為24的等邊三角形，構(gòu)型6是一個權(quán)重因子為36的等腰直角三角形，構(gòu)型12是一個權(quán)重因子為72的銳角等腰三角形，構(gòu)型23是一個權(quán)重因子為6的平角. 在體積較小時由于構(gòu)型1的比例占據(jù)較大，所以其他幾個重要的構(gòu)型(構(gòu)型6、構(gòu)型12、構(gòu)型23)在圖中作用不太明顯，但隨著體積的增大，構(gòu)型1的比例減小，這幾個重要的構(gòu)型也變得越來越重要. 因此，我們可以利用這幾個構(gòu)型結(jié)合其他所有構(gòu)型的數(shù)值進(jìn)行數(shù)據(jù)擬合.

圖2 不同體積下所有三體構(gòu)型的三體勢能Fig. 2 Potential energies of all three-body configurations at different volumes

3.1.3 氬三原子團(tuán)簇的特殊構(gòu)型的量子化學(xué)計算

由圖2可知，固氬三原子的幾個特殊構(gòu)型--構(gòu)型1(等邊三角形，頂角θ=60°)，構(gòu)型6(等腰直角三角形，頂角θ=90°)，構(gòu)型12(鈍角等腰三角形，頂角θ=120°)，構(gòu)型23(平角θ=180°)在總的三體勢能中所占比例很大，起到了重要作用. 因此，我們給出了由Pauli排斥效應(yīng)引起的交換勢部分U3S和長程色散勢部分U3C(通常指Axilrod-Teller項(xiàng)，簡稱AT項(xiàng))隨三角構(gòu)型的底邊長度的變化關(guān)系，如圖3所示.此外，與已有解析公式給出的交換勢部分和色散勢[14,16]部分進(jìn)行了比較. 結(jié)果表明，本文理論計算結(jié)果與較早文獻(xiàn)[16]有較大的偏差，但與經(jīng)過修正后的理論結(jié)果[14]在長程色散部分和短程交換部分都有很好的一致性. 若以這幾個特殊構(gòu)型的相關(guān)數(shù)據(jù)為基礎(chǔ)并結(jié)合其它構(gòu)型，則可擬合出更準(zhǔn)確的經(jīng)驗(yàn)勢函數(shù)及其參數(shù).

3.2 四原子構(gòu)型的相互作用勢能

3.2.1 Ar4四原子團(tuán)簇的構(gòu)型

3.2.2 不同體積下所有四體構(gòu)型的能量

由上面分析我們知道，四體勢的構(gòu)型共有28種，其中第8個構(gòu)型有兩種情況(如表2)，所以共有27種不同的構(gòu)型，我們比較了這28種四原子構(gòu)型在不同體積下四體勢能中的比例，如圖5所示. 任何一組四原子構(gòu)型的能量不僅與φ4有關(guān)，還決定于相應(yīng)的權(quán)重因子NG(見表2). 圖5表明構(gòu)型1到構(gòu)型8非常重要，對四體勢做出最大貢獻(xiàn)的是權(quán)重因子為8的第一種構(gòu)型，隨著體積的增大它的貢獻(xiàn)也在減少.

3.2.3 氬四原子團(tuán)簇的特殊構(gòu)型的量子化學(xué)計算

由上圖5可以看出，前幾個構(gòu)型在所有構(gòu)型中占的比例較大，起的作用較大，尤其是構(gòu)型1(正四面體)，其次是構(gòu)型2，構(gòu)型4，構(gòu)型5，構(gòu)型7和構(gòu)型8. 從兩體勢和三體勢的理論計算與較早文獻(xiàn)[13,14]比較可知，CCSD方法對于計算固氬的abintio數(shù)據(jù)具有較高的可靠性. 因此，在圖6中給出了貢獻(xiàn)較大的構(gòu)型1，構(gòu)型2，構(gòu)型4和構(gòu)型7的交換勢部分U4S、色散勢部分U4C以及二者之和四體勢能U4的結(jié)果隨最近原子間距R的變化關(guān)系，為擬合更準(zhǔn)確的四體經(jīng)驗(yàn)勢函數(shù)及其參數(shù)提供了重要的參考價值.

圖3 氬三體勢的CCSD(T)結(jié)果與經(jīng)驗(yàn)三體勢結(jié)果的比較Fig. 3 Comparison of CCSD results with the three-body experience potential results

表2 fcc晶格中Ar4團(tuán)簇的幾何參數(shù).三邊長r12、r13、r14以最近鄰原子間距R為單位，θ213是r12、r13間的夾角，θ214是r12、r14的夾角，θ314是r13、r14的夾角，NG表示權(quán)重因子

Table 2 Geometry parameters of Ar4clusters in fcc lattice. The lengths of three edgesr12,r13,r14are in unit of the nearest interatomic distanceR,θ213is the angle betweenr12andr13,θ214is the angle betweenr12andr14,θ314is the angle betweenr13andr14,NGrepresents the weight factor

Geometryr12r13r14θ213θ214θ314NG11116060608211160906048311160601202441129045451251116018012048611160901209671116012012024812121290901802412911218045135241011290451352411111901201202412112120459048131121204513548141111201201208151126090135481611260135135241711218090901218112909090241911290135135122012245909048211224513590242212245135180242311212090135482412290901801225122901359048261221351359012272229090908

4 結(jié) 論

本文運(yùn)用CCSD(T)的從頭計算方法和多體展開理論，從不同最近鄰殼層距離下任意三個原子間的三體構(gòu)型和任意四個原子間的四體構(gòu)型，對三體和四體的幾何構(gòu)型、幾何參數(shù)、不同體積下所有三體和四體構(gòu)型的能量以及各構(gòu)型所占比例等幾個方面對固氬fcc晶體的多體構(gòu)型作了系統(tǒng)的分析. 通過計算，我們可以看出，所有25組三體構(gòu)型中對總的三體勢能貢獻(xiàn)最大的是構(gòu)型1、構(gòu)型6、構(gòu)型12和構(gòu)型23；對三體勢及其交換部分和色散部分進(jìn)行了量子化學(xué)計算，計算結(jié)果與現(xiàn)有解析經(jīng)驗(yàn)勢在長程部分符合得非常好，但在短程部分有較小差異. 所有28組四體構(gòu)形中對總的四體勢能貢獻(xiàn)最大的是構(gòu)型1，構(gòu)型2，構(gòu)型4，構(gòu)型5，構(gòu)型7和構(gòu)型8；對四體勢及其交換勢部分和色散部分也進(jìn)行了量子化學(xué)計算，計算結(jié)果尚無解析經(jīng)驗(yàn)勢可比較. 若以這幾個特殊構(gòu)型的相關(guān)數(shù)據(jù)為基礎(chǔ)并結(jié)合其它構(gòu)型，可擬合出更準(zhǔn)確的三體經(jīng)驗(yàn)勢函數(shù)及其參數(shù)，也為擬合四體經(jīng)驗(yàn)勢函數(shù)及其參數(shù)提供了重要的參考價值.

圖4 fcc晶格中Ar4團(tuán)簇的幾何形狀Fig. 4 Geometries of Ar4 clusters in fcc lattice

圖5 不同體積下所有四體構(gòu)型的勢能Fig.5 Potential energies of all four-body potential configurations at different volumes

圖6 固氬的四體勢的量子化學(xué)計算Fig.6 Quantum chemistry calculation of four-body potential for solid argon

[1] Pepin R O. On the origin and early evolution of terrestrial planet atmospheres and meteoritic volatiles[J].Icarus, 1991, 92(2).

[2] Allegre C J, Hofmann A, Nions R O,etal. The argon constraints on mantle structure[J].Geophys.Res.Lett., 1996, 23: 3555.

[3] Mao H K, Xu J A, Bell P M,etal. Calibration of the ruby pressure gauge to 800 kbar under quasi-hydrostatic condition[J].Geophys.Res., 1986, 91: 4673.

[4] Pechenik E, Kelson I, Makov G,etal. Many-body model of rare gases at high pressures[J].Phys.Rev. B, 2008, 78(13): 134109.

[5] Cazorla C, Errandonea D, Sola E,etal. High-pressure phases, vibrational properties, and electronic structure of Ne(He)2and Ar(He)2: A First-principles study[J].Phys.Rev. B, 2009, 80(6): 064105-1.

[6] Tse J S, Klug D D, Shpakov V，etal. High pressure elastic properties of solid argon from first-principles density functional and quasi-harmonic lattice dynamic calculations[J].SolidStateCommunications, 2002, 1098(02): 557.

[7] Errandonea D, Boehler R, Japel S,etal. Structural transformation of compressed solid Ar: An x-ray diffraction study to 114 GPa[J].Phys.Rev. B, 2006, 73(9): 092106.

[8] Ross M, Mao H K, Bell P M,etal. The equation of state of dense argon: a comparison of shock and static studies[J].J.Chem.Phys., 1986, 85: 1028.

[9] Aziz R A, Slaman M J. The argon and krypton interatomic potentials revisited[J].Mol.Phys., 1986, 58(4): 679.

[10] Aziz R A, Slaman M J. The repulsive wall of the Ar-Ar interatomic potential reexamined[J].J.Chem.Phys., 1990, 92(2): 1030.

[11] Aziz R A. A highly accurate interatomic potential for argon[J].J.Chem.Phys., 1993, 99: 4518.

[12] Cybulski S M, Toczylowski R R. Ground state potential energy curves for He2, Ne2, Ar2, He-He, He-Ar, and Ne-Ar: A coupled-cluster study[J].J.Chem.Phys., 1999, 111(23): 10520.

[13] Slavicek P, Kalus R, Paska P,etal., State-of-the-art correlated ab initio potential energy curves for heavy rare gas dimers: Ar2, Kr2, and Xe2[J].J.Chem.Phys., 2003, 119(4): 2101.

[14] Freiman Y A, Tretyak S M. Many-body interactions and high-pressure equations of state in rare-gas solids[J].LowTemperaturePhysics, 2007, 33(6-7):545.

[15] Loubeyre P. Three-body exchange interaction in dense helium[J].Phys.Rev.Lett., 1987, 58: 1857.

[16] Loubeyre P. Three-body-exchange interaction in dense rare gases[J].Phys.Rev. B, 1988, 37(10):5432.

[17] Angela K W, David E W, Kirk A P,etal. Gaussian basis sets for use in correlated molecular calculations. IX. The atoms gallium through krypton[J].J.Chem.Phys., 1999, 110(16): 7667.

Many-body interaction of highly compressed solid argon based on atomic crystal configuration

ZHENG Xing-Rong1, FU Yun2, FU Wen-Yu1, LI Xiang-Fu1, LI Ji-Hong1

(1. Department of Physics, Longdong University, Qingyang 745000, China；2. College of Science, Guizhou Mizu University, Guiyang 550025, China)

X-ray diffraction experiments shows that solid argon is a face-centered cubic crystal structure (fcc structure). At present, crystal argon is studied involving two-, three- and four-body potentials. Using the many-body expansion method and the double cluster with full single and double excitations plus perturbative treatment of triples (CCSD(T)), the properties of fcc structure for solid argon which include geometrical configuration, geometrical parameters, potential energy of three- and four-body potentials and the proportion of all configurations at different volumes were accurately calculated. It is concluded that the configuration 1, 6, 12 and 23 play most important role in all three-body potential configurations. The calculation results of three-body potential, exchange potential and dispersion potential are in good agreement with the analytic experience potential in long-range part, but have small differences in short-range part. The configuration 1, 2, 4, 5, 7 and 8 play very important role in all four-body potential configurations, and there is not analytic experience potential to compare between the calculated results of four-body potential, exchange potential and dispersion potential. Using the data of these special configurations in combination with other configurations can fit a more accurate three-body experience potential function and its parameters, and also provide important reference value about fit four-body experience potential function and its parameters.

Solid argon; Face-centered cubic crystal structure; CCSD(T); Three-body potential; Four-body potential

2014-08-03

慶陽市自然科學(xué)基金(ZJ201306)

鄭興榮(1986—), 男，甘肅天水人，助教，主要從事材料計算研究.E-mail:zhengxingrong2006@163.com

103969/j.issn.1000-0364.2015.10.030

O641; O561.1

A

1000-0364(2015)05-0896-07