金屬銠的反演勢構建和應用

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

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

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

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

1 計算過程

1.1 反演勢的構建

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

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

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

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

對圖1 所示曲線進行擬合，其中E0表示孤立原子的基態能且有E0Rh= ?601.483 eV。

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

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

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

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

1.2 反演勢曲線的擬合

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

1.2.1 Rose 函數和Morse 函數的擬合

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

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

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

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

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

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

1.2.2 雙指數型勢函數的擬合

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

擬合得到的相關系數為1，表明擬合精度高，得到的參數如表2 所列。

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

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

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

2 應用和分析

2.1 銠的聲子譜的計算

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

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

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

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

2.2 勢函數的應用和計算

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

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

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

2.2.1 線膨脹系數的計算

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

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

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

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

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

由式(12)：

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

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

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

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

2.2.2 體彈性模量的計算

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

其中U為原子內聚能，V為單個原子在fcc 晶格結構中所占的體積，所以有：

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

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

2.2.3 格林乃森常數的計算

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

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

3 結論

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

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

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