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

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

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

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

1 計算過程

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

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

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

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

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

對圖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 晶格原胞中計算原子距離從0.2～0.7 nm，步長為0.01～0.05 nm 共計40 個格點的內(nèi)聚能。并減去基態(tài)能的值即得到晶格內(nèi)聚能曲線。如圖2 所示。根據(jù)以上計算結(jié)果用自編程序計算得到原子距離從0.2～1.2 nm，步長為0.02 nm 的晶格反演勢曲線，如圖3 所示。

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

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

1.2 反演勢曲線的擬合

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 應(yīng)用和分析

2.1 銠的聲子譜的計算

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

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

布里淵區(qū)對稱點坐標(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é)果表明，三種方法對應(yīng)的聲子譜曲線的變化趨勢是相似的，表明反演對勢可以有效的反應(yīng)原子間的相互作用。同時，采用Sutton-Chen 方法計算聲子譜時所需的時間比對勢方法所需的時間要多40倍，說明對勢函數(shù)在計算量上有明顯的優(yōu)勢。

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

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

采用式(8)雙指數(shù)型函數(shù)對圖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

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

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

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

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

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

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

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

由式(12)：

計算得到273～373 K 溫度段的銠的線膨脹系數(shù)的平均值為αLRh=0.71978×10-5K-1。與文獻(xiàn)[17]實驗數(shù)據(jù)對比αLRh=0.85×10-5K-1，計算結(jié)果的相對誤差為15.3%，說明計算數(shù)據(jù)與實驗數(shù)據(jù)相比基本符合，但值較小。這主要是因為忽略了金屬自由電子氣的影響。根據(jù)格林乃森方程：

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

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

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

2.2.2 體彈性模量的計算

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

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

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

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

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

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

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

3 結(jié)論

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

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

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