任縣利，陳 松，周建仁，謝 明，王塞北，陳靜洪，畢亞男

(昆明貴金屬研究所 稀貴金屬綜合利用新技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，昆明 650106)

原子間相互作用勢(shì)是凝聚態(tài)物質(zhì)在原子尺度上進(jìn)行計(jì)算模擬的基礎(chǔ)，特別是使用分子動(dòng)力學(xué)等方法對(duì)物質(zhì)的結(jié)構(gòu)和性質(zhì)進(jìn)行模擬研究時(shí)。獲得精確的原子間勢(shì)函數(shù)一直是模擬計(jì)算的重點(diǎn)和前提，并直接決定了模擬結(jié)果的準(zhǔn)確性和有效性。陳氏晶格反演勢(shì)理論，作為一種對(duì)勢(shì)理論，已經(jīng)被廣泛的應(yīng)用于各種材料的研究中，包括金屬間化合物，金屬陶瓷，離子晶體，半導(dǎo)體和金屬氫化物[1-5]等材料領(lǐng)域。陳氏晶格理論以數(shù)論中的莫比烏斯理論為基礎(chǔ)，可以獲得中心原子與任意近鄰下的原子之間的相互作用，從而得到精確的原子間相互作用勢(shì)，反演過程以數(shù)論為基礎(chǔ)，不含經(jīng)驗(yàn)成分，得到的結(jié)果是準(zhǔn)確有效的。而傳統(tǒng)的多體勢(shì)模型如EAM 勢(shì)(嵌入原子勢(shì))，勢(shì)函數(shù)中包含的參數(shù)往往通過擬合實(shí)驗(yàn)數(shù)據(jù)獲得，含有較大經(jīng)驗(yàn)性，并且在不同條件下，參數(shù)的值不同。同時(shí)，對(duì)于不同種材料，對(duì)應(yīng)有不同的勢(shì)函數(shù)形式。所以多體勢(shì)模型往往普適性不強(qiáng)[6-8]。

本文以貴金屬銠為研究對(duì)象，基于第一性原理和陳氏晶格反演理論，構(gòu)建銠的精確反演對(duì)勢(shì)。并利用對(duì)勢(shì)的計(jì)算結(jié)果計(jì)算銠的聲子譜以驗(yàn)證其有效性和可靠性。計(jì)算銠的線膨脹系數(shù)等物理量，為銠的研究尤其是在工程應(yīng)用上的研究提供參考。

1 計(jì)算過程

1.1 反演勢(shì)的構(gòu)建

基于莫比烏斯理論，陳氏晶格反演方法[9]將得到的晶格內(nèi)聚能反演，從而得到原子間的對(duì)勢(shì)，可以總結(jié)為以下兩個(gè)函數(shù)之間的轉(zhuǎn)換。如果：

其中E(x)為原子內(nèi)聚能，x為最近鄰距離，r(n)為第n近鄰原子的配位數(shù)，b(n)為第n近鄰原子的相對(duì)距離，那么：

函數(shù)(2)中的I(n)為反演系數(shù)，可以通過式(3)計(jì)算：

集合{b(m)}滿足乘法半群，并且通過計(jì)算機(jī)編程使得b(1),b(2)…b(n)按從小到大排序，并且b(1)=1。b-1[b(m)/b(n)]是一種數(shù)學(xué)運(yùn)算，表示當(dāng)b(m)/b(n)的值屬于集合{b(m)}并等于b(k)時(shí)，那么該運(yùn)算得到的值即為k。計(jì)算使用 Material Studio 軟件中的CASTEP 模塊，采用基于局域密度近似(LDA)的贗勢(shì)，計(jì)算了fcc(面心立方)金屬銠在不同原子距離下的孤立原子基態(tài)能。計(jì)算過程中布里淵區(qū)k空間網(wǎng)格劃分為16*16*16，截?cái)嗄茉O(shè)置為340 eV，迭代過程收斂精度選擇1*10-5eV。計(jì)算選取的價(jià)電子為4d85s1，計(jì)算過程不考慮自旋極化[10-12]。孤立原子基態(tài)能計(jì)算結(jié)果如圖1 所示。采用式(4)函數(shù)：

對(duì)圖1 所示曲線進(jìn)行擬合，其中E0表示孤立原子的基態(tài)能且有E0Rh= ?601.483 eV。

圖1 銠的孤立原子基態(tài)能曲線Fig.1 The isolated atomic ground state energy curve of Rh

在相同的參數(shù)設(shè)置條件下，在銠的fcc 晶格原胞中計(jì)算原子距離從0.2～0.7 nm，步長為0.01～0.05 nm 共計(jì)40 個(gè)格點(diǎn)的內(nèi)聚能。并減去基態(tài)能的值即得到晶格內(nèi)聚能曲線。如圖2 所示。根據(jù)以上計(jì)算結(jié)果用自編程序計(jì)算得到原子距離從0.2～1.2 nm，步長為0.02 nm 的晶格反演勢(shì)曲線，如圖3 所示。

圖2 銠的晶格內(nèi)聚能曲線Fig.2 The lattice cohesive energy curve of Rh

圖3 銠的晶格反演勢(shì)曲線Fig.3 The lattice inversion potential curve of Rh

1.2 反演勢(shì)曲線的擬合

精確的擬合函數(shù)，尤其是全局性符合程度高的函數(shù)，是下一步精確計(jì)算的基礎(chǔ)。采用Origin 軟件對(duì)圖3 所示曲線進(jìn)行擬合，擬合質(zhì)量通過軟件給出的相關(guān)系數(shù)進(jìn)行評(píng)估。相關(guān)系數(shù)的值在0～1 之間，值越接近于1，表明擬合的效果越好。本文分別采用Rose 函數(shù)，Morse 函數(shù)及本文提出的新型的雙指數(shù)型函數(shù)對(duì)反演勢(shì)曲線進(jìn)行擬合，并對(duì)擬合結(jié)果進(jìn)行對(duì)比和分析。

1.2.1 Rose 函數(shù)和Morse 函數(shù)的擬合

Rose 函數(shù)常被用來做對(duì)勢(shì)函數(shù)的擬合[13]，其函數(shù)形式如下：

其中D、R0和α為擬合得到的參數(shù)，φ(r)表示對(duì)勢(shì)勢(shì)能，r為最近鄰原子距離。擬合得到的相關(guān)系數(shù)值為0.99781，可見擬合質(zhì)量符合要求。

Morse 函數(shù)已被廣泛運(yùn)用于fcc 金屬的對(duì)勢(shì)函數(shù)的擬合中，其函數(shù)形式如式(6)。

擬合得到的參數(shù)值如表1 所列，表1 中同時(shí)將計(jì)算得到的數(shù)據(jù)與Flahive 等[14]計(jì)算的數(shù)據(jù)對(duì)比。

表1 利用Morse 函數(shù)擬合得到的銠參數(shù)對(duì)比Tab.1 Fitting parameters in Morse potential function of Rh

對(duì)比發(fā)現(xiàn)：R0和α的值相對(duì)差別較小，而D的值相差較大。這主要是兩者選擇的計(jì)算溫度不同，本文計(jì)算選取0 K 為基準(zhǔn)態(tài)；而文獻(xiàn)[14]的數(shù)據(jù)計(jì)算基準(zhǔn)態(tài)為氣態(tài)，并且文獻(xiàn)[14]中R0由平衡態(tài)最近鄰原子距離決定，α與勢(shì)函數(shù)的二階求導(dǎo)有關(guān)，盡管隨著溫度的變化，原子間距離會(huì)以平衡態(tài)最近鄰距離為中心做熱振動(dòng)，導(dǎo)致晶格常數(shù)發(fā)生改變，但根本上平衡態(tài)最近鄰原子距離始終保持不變，所以當(dāng)擬合選取的溫度不同時(shí)，這兩個(gè)參數(shù)的值差別較小。而D值與勢(shì)函數(shù)對(duì)應(yīng)的最低勢(shì)能有關(guān)，由勢(shì)函數(shù)選取的基準(zhǔn)態(tài)決定，并且氣態(tài)條件下，基準(zhǔn)態(tài)能量要小于0 K 時(shí)孤立原子的基態(tài)能，所以導(dǎo)致了D值有0.0896 eV 的偏差，且文獻(xiàn)[14]對(duì)應(yīng)的值較小。所以本研究得到的勢(shì)函數(shù)是準(zhǔn)確有效的。

1.2.2 雙指數(shù)型勢(shì)函數(shù)的擬合

為了進(jìn)一步提高函數(shù)擬合精度，本文提出了包含有5 個(gè)參數(shù)的雙指數(shù)型勢(shì)函數(shù)，其形式如下：

擬合得到的相關(guān)系數(shù)為1，表明擬合精度高，得到的參數(shù)如表2 所列。

表2 利用雙指數(shù)型函數(shù)擬合得到銠的參數(shù)Tab.2 Fitting parameters in double exponential potential function of rhodium

為了對(duì)比和分析3 種函數(shù)擬合對(duì)勢(shì)曲線的效果，將其擬合結(jié)果進(jìn)行放大觀察和分析。如圖4 所示。結(jié)果表明，Rose 函數(shù)和Morse 函數(shù)雖然在對(duì)勢(shì)函數(shù)曲線的長程段和短程段擬合效果較好，但是在曲線的拐點(diǎn)附近與曲線的一致性符合較差。所以兩個(gè)函數(shù)的全局性精度相對(duì)不高。而本文提出的雙指數(shù)型勢(shì)函數(shù)，在擬合的全局上表現(xiàn)非常優(yōu)異，且相關(guān)系數(shù)為1，表明擬合精度很高，為下一步的精確計(jì)算，提供有力的基礎(chǔ)。

圖4 三種勢(shì)函數(shù)擬合結(jié)果的比較Fig.4 The comparison of the fitting results of these three potential functions

2 應(yīng)用和分析

2.1 銠的聲子譜的計(jì)算

為了驗(yàn)證對(duì)勢(shì)函數(shù)的有效性，通過Material Studio 軟件中的GULP 模塊分別采用反演對(duì)勢(shì)方法和EAM 勢(shì)中的Sutton-Chen 多體勢(shì)方法以及通過軟件中的CASTEP 模塊(采用有限位移法)計(jì)算了銠的聲子譜，如圖5 所示。

圖5 銠的聲子譜Fig.5 The phonon spectra of Rh

布里淵區(qū)對(duì)稱點(diǎn)坐標(biāo)為G(0, 0, 0)、X(0.5, 0, 0.5)、W(0.5, 0.25, 0.75)、K(0.75, 0.375, 0.375)、L(0.5, 0.5, 0.5)。結(jié)果表明，三種方法對(duì)應(yīng)的聲子譜曲線的變化趨勢(shì)是相似的，表明反演對(duì)勢(shì)可以有效的反應(yīng)原子間的相互作用。同時(shí)，采用Sutton-Chen 方法計(jì)算聲子譜時(shí)所需的時(shí)間比對(duì)勢(shì)方法所需的時(shí)間要多40倍，說明對(duì)勢(shì)函數(shù)在計(jì)算量上有明顯的優(yōu)勢(shì)。

對(duì)于不同的材料體系，EAM 勢(shì)的函數(shù)模型不同且均包含有經(jīng)驗(yàn)成分[15]。并且其函數(shù)推導(dǎo)過程往往以Rose 函數(shù)為基礎(chǔ)，導(dǎo)致在后續(xù)的計(jì)算精度不夠。而反演對(duì)勢(shì)函數(shù)基于數(shù)學(xué)理論，通過嚴(yán)格的數(shù)學(xué)證明，不含有任何經(jīng)驗(yàn)成分，是一種全局精確的勢(shì)函數(shù)。對(duì)于聲子譜的計(jì)算結(jié)果也是有效和可靠的。

2.2 勢(shì)函數(shù)的應(yīng)用和計(jì)算

采用式(8)雙指數(shù)型函數(shù)對(duì)圖2 所示的銠的晶格內(nèi)聚能曲線進(jìn)行擬合，擬合得到的參數(shù)如表3，擬合相關(guān)系數(shù)為0.99996。

表3 利用雙指數(shù)型函數(shù)擬合內(nèi)聚能曲線得到的銠參數(shù) Tab.3 Fitting parameters of cohesive curve using double exponential potential function for rhodium

對(duì)該函數(shù)求一階導(dǎo)數(shù)，并令其值為0，得到0 K下，平衡態(tài)原子最近鄰距離r0Rh= 0.2686139 nm。在此基礎(chǔ)上，分別計(jì)算了銠的線膨脹系數(shù)，體彈性模量以及格林乃森常數(shù)等物理量，并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比和分析。

2.2.1 線膨脹系數(shù)的計(jì)算

本方法基于計(jì)算得到的勢(shì)函數(shù)的解析式以及玻爾茲曼分布函數(shù)[16]，計(jì)算在不同溫度下，銠的熱振動(dòng)平均位移。

其中δ=r?r0，r0為0 K 下平衡態(tài)原子最近鄰距離。考慮到原子距離不可能為負(fù)值，且從a到∞變化，其中a為銠的離子半徑。

由于線膨脹系數(shù)是一維物理量，V(r)應(yīng)為一維原子鏈上的勢(shì)能，而公式(8)中的u(r)為原子間總能，是三維的量。所以V(r)與u(r)之間的關(guān)系應(yīng)為

將式(8)及表3 中的數(shù)據(jù)代入到式(10)中，用自編程序計(jì)算得到273～373 K 范圍內(nèi)不同溫度下的原子平均熱振動(dòng)位移的數(shù)值解。結(jié)果如圖6 所示。

圖6 銠的原子熱振動(dòng)位移-溫度曲線Fig.6 Relationship between average atomic thermal vibration displacement and temperature curve of Rh

由式(12)：

計(jì)算得到273～373 K 溫度段的銠的線膨脹系數(shù)的平均值為αLRh=0.71978×10-5K-1。與文獻(xiàn)[17]實(shí)驗(yàn)數(shù)據(jù)對(duì)比αLRh=0.85×10-5K-1，計(jì)算結(jié)果的相對(duì)誤差為15.3%，說明計(jì)算數(shù)據(jù)與實(shí)驗(yàn)數(shù)據(jù)相比基本符合，但值較小。這主要是因?yàn)楹雎粤私饘僮杂呻娮託獾挠绊?。根?jù)格林乃森方程：

其中κ為體彈性模量，αV為體膨脹系數(shù)且與線膨脹系數(shù)之間的關(guān)系為αV=3αL，γ為格林乃森常數(shù)，CV為定容比熱容。金屬的比熱容包括晶格比熱容和電子比熱容兩部分[16]，其表達(dá)式如下：

當(dāng)溫度高于德拜溫度時(shí)，晶格比熱容起主導(dǎo)作用，但是當(dāng)溫度較低時(shí)，電子對(duì)金屬的比熱容會(huì)有顯著貢獻(xiàn)，所以在273～373K 范圍內(nèi)，金屬自由電子對(duì)比熱容的影響不能忽略。所以導(dǎo)致計(jì)算結(jié)果值偏小。

但該計(jì)算方法是有效且十分有意義的，傳統(tǒng)的計(jì)算線膨脹系數(shù)的理論方法，主要依靠近似原子間作用勢(shì)和分子動(dòng)力學(xué)方法等，每種方法都有自己的假設(shè)，許多假設(shè)只有在有限的溫度范圍內(nèi)才有效[18]，且計(jì)算量大，流程復(fù)雜。本方法結(jié)合原子間相互作用勢(shì)和玻爾茲曼分布函數(shù)，計(jì)算量小，步驟簡單，具有一定精度，適合實(shí)際使用。

2.2.2 體彈性模量的計(jì)算

體彈性模量和原子間相互作用勢(shì)的關(guān)系[19]可以用下式表示：

其中U為原子內(nèi)聚能，V為單個(gè)原子在fcc 晶格結(jié)構(gòu)中所占的體積，所以有：

將公式(8)帶入到上式得：

由公式(10)計(jì)算得到在室溫293 K 下原子的平均熱振動(dòng)位移為δ-Rh=5.51×10-4nm，所以室溫下原子平衡態(tài)最近鄰原子距離為r0'=0.2719819 nm，r0'由公式r0'=r0+δ-計(jì)算得到。帶入到公式(16)中，得室溫(293K)下銠的體彈性模量為κRh=3.132099×1011N·m-2，與實(shí)驗(yàn)數(shù)據(jù)[17]κRh=2.76×1011N·m-2對(duì)比，誤差為13.5%，說明計(jì)算結(jié)果是準(zhǔn)確有效的。

2.2.3 格林乃森常數(shù)的計(jì)算

基于公式(13)及銠的線膨脹系數(shù)和體彈性模量的計(jì)算值，其中體積V 通過下式計(jì)算：

NA為阿伏伽德羅常數(shù)，r0為平衡態(tài)最近鄰原子距離。且CVRh=24.98 J/mol·K，計(jì)算得到293 K 下，銠的格林乃森常數(shù)為γRh=2.25，與實(shí)驗(yàn)數(shù)據(jù)[17]對(duì)比γRh=2.26，相對(duì)誤差為0.8%，說明得到的勢(shì)函數(shù)勢(shì)準(zhǔn)確有效[20-23]。

3 結(jié)論

1) 以第一性原理計(jì)算為基礎(chǔ)，通過陳氏晶格反演方法，得到了fcc 金屬銠的精確的晶格反演對(duì)勢(shì)曲線，并采用本研究提出的雙指數(shù)型勢(shì)函數(shù)對(duì)曲線進(jìn)行擬合，得到對(duì)勢(shì)函數(shù)的精確解析式。

2) 分別采用反演對(duì)勢(shì)函數(shù)，Sutton-Chen 函數(shù)及有限位移法計(jì)算了銠的聲子譜，對(duì)比發(fā)現(xiàn)本文得到的勢(shì)函數(shù)是有效的。

3) 提出了基于反演勢(shì)和玻爾茲曼分布函數(shù)計(jì)算線膨脹系數(shù)的方法，計(jì)算了銠的線膨脹系數(shù)、體彈性模量和格林乃森常數(shù)，與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比相對(duì)誤差分別為15.3%、13.5%和0.8%。表明本研究提出的計(jì)算方法是準(zhǔn)確且有效的。