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

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

基于原子晶體構型的高壓固氬中的多體相互作用

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

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

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

固氬； fcc晶格結構； CCSD(T)； 三體勢； 四體勢

1 引 言

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

2 理論模型

固氬晶體中，任一原子被其不同鄰近原子包圍. 由于屏蔽效應，任一原子o對結合排斥勢的貢獻主要來自于o與周圍近鄰殼層(i=1,2,3,…)的有限的(n-1)個近鄰原子的相互作用. 原子o和這(n-1)個原子構成一n原子的團簇Arn，該團簇的勢能等于這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位置時團簇的基態(tài)總能量，它可以通過求解體系的薛定諤方程得到

φ=Eφ.

(2)

對于處在位置坐標r0和ri下的兩氬原子o和i，它們所構成團簇的相互作用能，可用一個兩體勢描述

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

(3)

其中u2(o,i)代表兩氬原子o，i間的兩體相互作用，E(r0,ri)是兩個氬原子分別位于r0，r0位置時它們所構成原子團簇Ar2的基態(tài)能量，E0代表一個孤立氬原子的基態(tài)能量. 當?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位置時三個氬原子所構成團簇Ar3的基態(tài)能量. 同理，當?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)分別為各項對所有的近鄰原子求和. 在實際求解過程中我們運用Gamess計算程序來求出氬原子團簇的多體勢能，并運用CCSD(T)的從頭計算方法和多體展開理論，對fcc固氬的性質(zhì)作了研究.

3 結果與討論

3.1 三原子構型的相互作用勢能

3.1.1 Ar3三原子團簇的構型

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

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

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 不同體積下所有三體構型的能量

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

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

3.1.3 氬三原子團簇的特殊構型的量子化學計算

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

3.2 四原子構型的相互作用勢能

3.2.1 Ar4四原子團簇的構型

3.2.2 不同體積下所有四體構型的能量

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

3.2.3 氬四原子團簇的特殊構型的量子化學計算

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

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

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

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 結 論

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

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

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

圖6 固氬的四體勢的量子化學計算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

慶陽市自然科學基金(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