金屬Nb的Finnis-Sinclair勢開發(fā)及勢函數(shù)形式對材料性能的影響*

高靜怡 孫嘉興 王遜 周剛 王皞 劉艷俠? 徐東生

1) (遼寧大學(xué)物理學(xué)院， 沈陽 110036)

2) (沈陽建筑大學(xué)理學(xué)院， 沈陽 110168)

3) (中國科學(xué)院金屬研究所， 沈陽 110016)

對于計(jì)算材料科學(xué)的研究者來說， 經(jīng)常由于找不到合適的原子間勢而工作受阻.本文將在Finnis-Sinclair勢的框架下， 通過開發(fā)金屬Nb的Finnis-Sinclair勢而給出較詳細(xì)的原子間勢擬合、檢驗(yàn)、修正的過程.首先建立原子間勢與材料宏觀性能之間的關(guān)系， 然后通過再現(xiàn)金屬Nb的結(jié)合能、體模量、表面能、空位形成能及平衡點(diǎn)陣常數(shù)的實(shí)驗(yàn)數(shù)據(jù)的方法擬合金屬Nb的Finnis-Sinclair勢.利用所構(gòu)建的原子間勢計(jì)算金屬Nb的彈性常數(shù)、剪切模量及柯西壓力來檢驗(yàn)勢函數(shù).討論勢函數(shù)曲線形狀對間隙形成能的影響， 進(jìn)而根據(jù)間隙能的計(jì)算數(shù)據(jù)修正已構(gòu)建的原子間勢.討論截?cái)嗑嚯x的處理方法.本文的結(jié)果一方面為構(gòu)建原子間勢函數(shù)庫提供資料， 為構(gòu)建與Nb相關(guān)的合金原子間勢奠定基礎(chǔ); 另一方面， 為開發(fā)和改善原子間勢質(zhì)量提供方法和依據(jù).

1 引 言

原子尺度模擬是微觀層次研究材料性能行之有效的手段和途徑， 但是計(jì)算機(jī)模擬結(jié)果的可靠性直接來源于原子間勢的質(zhì)量[1]， 而且模擬的計(jì)算量強(qiáng)烈依賴于原子間勢的復(fù)雜度.因此開發(fā)形式簡單、質(zhì)量可靠的原子間勢及建立易于操作的原子間勢構(gòu)建方法極為重要.目前， 原子間勢的開發(fā)大多通過有限的勢函數(shù)參數(shù)來再現(xiàn)材料的部分性質(zhì)而獲得[2-4]， 因此， 根據(jù)部分性能獲得的原子間勢不可能描述材料的全部性能.由于研究目的不同， 有時(shí)即使是同一金屬[5-10]或合金[11，12]也需要開發(fā)不同的原子間勢.

Nb基合金及含Nb合金在航空航天、醫(yī)學(xué)、核工業(yè)等領(lǐng)域有著重要的應(yīng)用[13-15].為了從原子尺度研究Nb及其合金的性能， 首先需要構(gòu)建其原子間勢.在金屬及合金中應(yīng)用廣泛的原子間勢模型主要是1983—1984年間由Daw和Baskes[16，17]提出的Embedded-atom method (EAM)勢模型， 以及1984年由Finnis和 Sinclair[2]提出的Finnis-Sinclair (FS)勢.二者數(shù)學(xué)表達(dá)式都由兩項(xiàng)組成，第一項(xiàng)都為排斥勢， 解釋相同， 都是原子芯之間的靜電勢; 第二項(xiàng)都為吸引勢.二者的主要區(qū)別在于第二項(xiàng)， 體現(xiàn)在兩個(gè)方面: 一是來源不同， EAM勢來源于密度泛函理論， FS勢來源于緊束縛理論的二階矩近似; 二是對第二項(xiàng)即電荷密度函數(shù)和嵌入能函數(shù)的解釋不同， EAM勢中， 把金屬中的每個(gè)原子看成是鑲嵌在其他原子在該位置形成的電子氣中的一個(gè)雜質(zhì)， 并假設(shè)嵌入能是局域電子密度及其高階導(dǎo)數(shù)的函數(shù)， 且電子密度由原子的電荷密度疊加而成; 而FS勢中， 每個(gè)原子內(nèi)聚能的大小與原子間最近鄰的鍵合數(shù)的平方成正比， 在緊束縛理論中， 金屬中電子能帶的帶能是其中占據(jù)態(tài)單電子能之和， 而且?guī)苡蓱B(tài)密度的二階矩來表征， 因此密度函數(shù)是態(tài)密度的二階矩.雖然這兩種原子間勢的來源和解釋不同， 但由于二者的數(shù)學(xué)形式相近，為此， 人們習(xí)慣上將FS勢歸入為EAM模型.本文為了敘述上的方便， 仍然稱作FS勢.到目前為止，在金屬及合金中廣泛使用的原子間勢形式依然是EAM勢及FS勢， 以及對這兩種勢的各種修正勢[18-20].

1984年Finnis和 Sinclair[2]開發(fā)了7種bcc結(jié)構(gòu)金屬Nb， V， Ta， Cr， Mo， W及Fe的FS勢， 首次考慮了原子之間的多體效應(yīng)， 克服了對勢的缺點(diǎn).但該勢的缺點(diǎn)是在近距離下原子之間的排斥力不足， 特別是V和Nb， 在短距離內(nèi)甚至表現(xiàn)出有吸引力， 導(dǎo)致分子動(dòng)力學(xué)模擬時(shí)出現(xiàn)原子落在一起的異常行為[21].為此， Ackland[21]和Rebonato[22]分別對該問題進(jìn)行了修正， 對于小于第一近鄰的對勢曲線進(jìn)行了硬化.另一方面， 該勢的原子間相互作用距離考慮到第二近鄰， 相互作用距離比較短.本文在該勢的基礎(chǔ)上， 首先， 擴(kuò)展相互作用距離到第三近鄰， 并修正近距離相互作用偏軟的現(xiàn)象; 其次， 建立原子間勢與材料宏觀性質(zhì)之間的關(guān)系， 討論函數(shù)形式及截?cái)嗑嚯x對宏觀性能的影響， 并根據(jù)材料的宏觀性質(zhì)修正原子間勢函數(shù)形式.

構(gòu)建原子間勢是一項(xiàng)復(fù)雜的工程， 本文將從模型的選擇， 截?cái)嗑嚯x處理， 勢函數(shù)擬合， 勢函數(shù)檢驗(yàn)及應(yīng)用， 以及根據(jù)檢驗(yàn)和應(yīng)用結(jié)果修正勢函數(shù)幾方面詳細(xì)介紹勢函數(shù)的開發(fā)過程， 并就截?cái)嗑嚯x和勢函數(shù)形式對原子間勢的質(zhì)量及材料性能影響方面進(jìn)行了討論， 期望為開發(fā)和改善原子間勢質(zhì)量提供參考.

2 模 型

本文在FS勢框架下[2]開發(fā)金屬Nb的原子間勢， 將原子間相互作用擴(kuò)展到第三近鄰.對于理想單質(zhì)金屬， 所有原子都是等價(jià)的， 因此體系的總能可以表示為單原子能量之和， 單個(gè)原子的相互作用能表示為:

其中，uP表示原子對之間的相互作用能，uN表示多體相互作用對總能的貢獻(xiàn).它們可用(2)式和(3)式表示:

其中，A為勢參數(shù)， 是大于零的常數(shù).根據(jù)緊束縛理論， 取ρ為參考原子i處的局域電子密度， 表示如下:

其中，φ(ri)為重疊積分的平方和[2]或?yàn)樵与姾擅芏萚17].

原始模型中的φ(ri) 為拋物線型函數(shù)， 當(dāng)將截?cái)嗑嚯x擴(kuò)展到第三近鄰時(shí)， 使用原函數(shù)形式給出了不合理的原子間相互作用， 為此本文修正該函數(shù)為四次冪形式:

原始模型中的對勢函數(shù)V(r)為四次多項(xiàng)式， 本文修正為六次多項(xiàng)式:

函數(shù)形式的選擇和c與d的取值在后面詳細(xì)討論.c與d是截?cái)鄥?shù)， 選取在第三近鄰與第四近鄰之間，c和d的值分為:

其中，a為Nb的晶格常數(shù).

3 截?cái)嗵幚矸椒?/h2>

為了減少分子動(dòng)力學(xué)模擬計(jì)算原子受力所用時(shí)間， 通常不采用相互作用減少到零的自然截?cái)嗟姆绞剑?而是在一個(gè)方便的力程上將原子間相互作用截?cái)?有研究表明， 如果預(yù)先不做截?cái)嗵幚恚?計(jì)算一個(gè)分子動(dòng)力學(xué)步的99%的時(shí)間都將用于計(jì)算使粒子運(yùn)動(dòng)所需的力.

如果勢能不是光滑截?cái)啵?那么在截?cái)帱c(diǎn)上的力就會(huì)出現(xiàn)δ函數(shù)形式的奇異性， 因此， 必須考慮勢能的截?cái)鄬ο到y(tǒng)特性的影響.首先需要保證勢能及其一階、二階導(dǎo)數(shù)在截?cái)帱c(diǎn)連續(xù).勢能對距離的一階導(dǎo)數(shù)對力有影響， 勢能對距離的二階導(dǎo)數(shù)對材料的彈性常數(shù)、體模量、剪切模量等有影響.這些函數(shù)是否連續(xù)、力程的大小對材料的結(jié)構(gòu)穩(wěn)定性、力學(xué)穩(wěn)定性、相變等都有影響.

可以采取多種方法對勢函數(shù)進(jìn)行截?cái)啵?如:1)簡單地將函數(shù)平移; 2)將函數(shù)在適當(dāng)位置截?cái)啵僭黾咏匚埠瘮?shù); 3)在適當(dāng)位置迭加一個(gè)截尾函數(shù); 4) 將截?cái)嗪瘮?shù)直接包含到勢函數(shù)表達(dá)式中.

其中方法1)比較粗糙， 且不能保證勢函數(shù)的一、二階導(dǎo)數(shù)連續(xù)， 只適合做靜態(tài)的與能量相關(guān)的計(jì)算.方法2)和方法3)需要增加參數(shù)的數(shù)量， 以保證在函數(shù)連接點(diǎn)的各階導(dǎo)數(shù)連續(xù)， 且要慎重考慮連接點(diǎn)的選取.方法4)是目前使用比較多的方法，優(yōu)點(diǎn)是與勢函數(shù)同時(shí)進(jìn)行參數(shù)擬合.目前常用的主要有兩種形式， 一是Finnis和Sinclair使用的階躍函數(shù)的形式

即在選取的勢函數(shù)上直接乘以H(x); 另一種是Mishin[23]使用的截?cái)喾绞剑?/p>

也是在勢函數(shù)上直接乘以ψ(x).

4 原子間勢與晶體性質(zhì)

不論是擬合原子間勢還是檢驗(yàn)原子間勢， 均需建立原子間勢與擬合晶體性質(zhì)參數(shù)及檢驗(yàn)晶體性質(zhì)參數(shù)之間的關(guān)系.本文建立了內(nèi)聚能、平衡點(diǎn)陣常數(shù)、體模量、表面能、空位形成能、間隙能及切變模量之間的關(guān)系.

對于體心立方晶體， 設(shè)晶格常數(shù)為a， 平衡態(tài)下的原子體積，Ωe=a3/2.本文構(gòu)建的原子間勢的相互作用在原子的第三和第四近鄰之間截?cái)?bcc晶體第一近鄰有8個(gè)原子，距離為第二近鄰有6個(gè)原子， 距離為a; 第三近鄰有12個(gè)原子，距離為

采用(10)式的符號來標(biāo)記不同近鄰下的對相互作用及電子密度，

4.1 壓力和體模量

其中，PN和PP分別表示為N體相互作用和兩體相互作用部分的力.由

得

當(dāng)r=a時(shí)， 可得到平衡狀態(tài)下的壓力:

由

得

當(dāng)r=a時(shí)， 可得到平衡狀態(tài)下的體模量:

4.2 表面能

表面能Esurf的計(jì)算公式為:

其中，Es為具有表面時(shí)體系的總能，Et為不包含表面時(shí)體系的總能，S為計(jì)算體系表面的面積.考慮(100)面截?cái)嗑嚯x內(nèi)原子能量的變化， (100)面的表面能為:

其中，

4.3 未弛豫空位形成能

單空位形成能為:

其中，En—1為具有n個(gè)原子的晶體在其內(nèi)部產(chǎn)生一個(gè)空位時(shí)體系的總能，En為具有n個(gè)原子的理想晶體體系的總能.考慮空位周圍截?cái)嗑嚯x內(nèi)原子能量的變化， 即為單空位形成能.

4.4 內(nèi)聚能

在FS勢下， 每個(gè)原子的內(nèi)聚能為:

4.5 平衡條件

平衡狀態(tài)下， 晶體內(nèi)各原子受力為零， 總能量最低， 即utot關(guān)于r的一階導(dǎo)數(shù)為0.因此， 當(dāng)r=a時(shí)的平衡條件為:

其中，

4.6 自間隙形成能

自間隙形成能的計(jì)算公式為:

其中，En表示具有n個(gè)粒子的理想晶體的總能量，En+1表示存在一個(gè)間隙原子時(shí)體系的總能量.本文考慮圖1所示的6種間隙構(gòu)型.

本文采用計(jì)算缺陷范圍內(nèi)的原子間相互作用能量變化的方法計(jì)算每種間隙的形成能.八面體、四面體及擠列子3種間隙構(gòu)型中， 在截?cái)嗑嚯x內(nèi)的每個(gè)原子除了距離的原子有8個(gè)， 距離為a的原子有6個(gè)， 距離為的原子有12個(gè)外， 還增加了一個(gè)間隙原子的相互作用， 截?cái)嗑嚯x內(nèi)的這些原子與間隙原子的距離及等價(jià)原子數(shù)見表1.啞鈴構(gòu)型中， 截?cái)嗑嚯x內(nèi)的原子分布要復(fù)雜一些， 可通過Materials Studio軟件建立相應(yīng)的構(gòu)型， 再分析近鄰原子分布情況.

圖1 bcc結(jié)構(gòu)6種間隙構(gòu)型Fig.1.Six interstitial configurations of bcc structure.

表1 截?cái)嗑嚯x內(nèi)的各間隙原子的距離及等價(jià)原子數(shù)Table 1.Distance and equivalent atomic number of each atom withinthe cutoff distance from the interstitial atom.

根據(jù)間隙原子的近鄰分布可計(jì)算間隙原子的形成能.

八面體間隙形成能為:

其余構(gòu)型的間隙形成能可采用同樣方法計(jì)算.

4.7 立方晶體的切變模量C44和C'

我們知道， 應(yīng)變能密度為

其中，σij=Cijklεij.

按照Voigt記法， 可得應(yīng)變能密度為

立方晶體有兩個(gè)切變模量C44和C′(≡(C11—C12)/2).

對于C′， 考慮立方晶體的平面應(yīng)變，X方向伸長時(shí)應(yīng)變?yōu)棣?1，Y方向縮短時(shí)應(yīng)變?yōu)棣?2，Z方向不變時(shí)應(yīng)變?yōu)?.保持體積不變， 則ε11= —ε22.對于C44， 相當(dāng)于在[011]方向的拉伸， 只有應(yīng)變?yōu)棣?2， 將應(yīng)變分量代入(24)式， 得到兩種應(yīng)變下的應(yīng)變能密度分別為

可得

應(yīng)變能密度W為單位體積的內(nèi)能， 則

經(jīng)推導(dǎo)(見附錄A)得

根據(jù)(31)式給出的體模量、切變模量， 柯西壓與彈性常數(shù)之間的關(guān)系， 可在這些量之間相互計(jì)算，

5 原子間勢擬合

勢函數(shù)表達(dá)式(3)式和(6)式中的c0，c1，c2，A是勢函數(shù)參數(shù)， 本文通過擬合金屬Nb的內(nèi)聚能、平衡點(diǎn)陣常數(shù)、體模量、表面能及空位形成能的實(shí)驗(yàn)數(shù)據(jù)獲得， 實(shí)驗(yàn)數(shù)據(jù)[2]見表2中第1行數(shù)據(jù).

根據(jù)第4節(jié)中建立的晶體性能與勢函數(shù)的關(guān)系， 然后采用均方差最小方法， 調(diào)整勢函數(shù)參數(shù)的大小， 使得勢函數(shù)能再現(xiàn)這些晶體性能， 擬合得到的原子間勢參數(shù)見表3中第1行數(shù)據(jù)， 使用所得的原子間勢計(jì)算的晶體性能列于表2中第2行.

從表3可以看出， 擬合的均方差非常小， 因此，利用所建原子間勢計(jì)算的晶格常數(shù)、內(nèi)聚能、未弛豫空位形成能、 體模量的結(jié)果與實(shí)驗(yàn)值完全相同，只有表面能有微小差別.

表2 擬合用金屬Nb的實(shí)驗(yàn)數(shù)據(jù)及計(jì)算結(jié)果Table 2.Experimental and calculation data of metal Nb for fitting interatomic potential.

表3 金屬Nb的FS勢參數(shù)及擬合均方差Table 3.FS potential parameters of metal Nb and fitting mean square error.

6 勢函數(shù)的檢驗(yàn)與應(yīng)用

6.1 彈性常數(shù)

根據(jù)第3節(jié)建立的勢函數(shù)與晶體性能的關(guān)系，利用所得勢函數(shù)計(jì)算了金屬Nb的3個(gè)獨(dú)立的彈性常數(shù)及柯西壓力.結(jié)果見表4中第2行數(shù)據(jù).從表4中數(shù)據(jù)可以看出， 除了C44偏差較大外， 其余結(jié)果均較接近實(shí)驗(yàn)值.

6.2 間隙形成能

根據(jù)第3節(jié)建立的勢函數(shù)與晶體性能的關(guān)系，利用所得勢函數(shù)計(jì)算了Nb的不同間隙構(gòu)型的間隙形成能， 結(jié)果見表5.其中， FS(87)一列數(shù)據(jù)是文獻(xiàn)[22]作者使用修正的FS勢計(jì)算的弛豫的間隙形成能.為了便于與本文的數(shù)據(jù)進(jìn)行比較， 本文使用文獻(xiàn)[22]的勢函數(shù)計(jì)算了各種間隙構(gòu)型的未馳豫間隙形成能， 列于表4中“FS(87)未弛豫”一列.

從表5可以看出， 各種方法計(jì)算的間隙形成能中， 或者〈110〉啞鈴構(gòu)型的間隙能最低， 或者是〈111〉啞鈴構(gòu)型的最低， 本文的計(jì)算結(jié)果也是〈110〉啞鈴構(gòu)型最低.但是本文的計(jì)算結(jié)果普遍偏大， 分析其原因， 一方面我們的計(jì)算結(jié)果均未進(jìn)行弛豫，另一方面， 說明我們所建的勢函數(shù)在近距離處排斥力過大， 曲線偏硬.

表4 金屬Nb的彈性常數(shù)(單位為1011 Pa)Table 4.Elastic constants of metal Nb (in 1011 Pa).

7 勢函數(shù)的修正

為了使所建立的勢函數(shù)能更好地描述間隙形成能， 我們對勢函數(shù)在近距離的行為進(jìn)行修正， 嘗試加入修正項(xiàng)進(jìn)行軟化處理， 即將對勢函數(shù)表達(dá)式(6)式在第一近鄰之內(nèi)減去一項(xiàng)，

其中修正項(xiàng)形式如下， 且僅在第一近鄰內(nèi)有效，

按照前述方法重新擬合勢參數(shù)并進(jìn)行檢驗(yàn).發(fā)現(xiàn)加入修正項(xiàng)后擬合得到的勢參數(shù)幾乎沒有變化，利用修正后的勢函數(shù)計(jì)算的各構(gòu)型的間隙形成能均減小， 但與文獻(xiàn)數(shù)據(jù)仍有一定的偏差， 為此繼續(xù)對勢函數(shù)進(jìn)行修正.

將間隙形成能的計(jì)算值[24]作為目標(biāo)值， 調(diào)整修正項(xiàng)的冪次， 當(dāng)取修正項(xiàng)形式為(34)式時(shí)， 計(jì)算結(jié)果接近目標(biāo)值.

修正前后的對勢曲線、有效對勢曲線及電子密度曲線如圖2所示， 為了便于比較， 同時(shí)繪出了文獻(xiàn)[2]， 及按文獻(xiàn)[2]函數(shù)形式將截?cái)嗑嚯x擴(kuò)展到第三近鄰的函數(shù)曲線.

表5 金屬Nb的間隙形成能Table 5.Interstitial formation energy of metal Nb.

由圖2(b)和圖2(c)可以看出， 原始勢函數(shù)曲線(文獻(xiàn)[2])在近距離處表現(xiàn)出吸引行為， 文獻(xiàn)[21]通過在小于第一近鄰范圍疊加一項(xiàng)B(b0—rij)3exp(—αrij)函數(shù)來增加排斥力， 修正了這種偏軟現(xiàn)象， 其中B，b0，α為勢參數(shù); 文獻(xiàn)[22]通過對比分析由原子間勢計(jì)算的力和實(shí)驗(yàn)數(shù)據(jù)之間的差別， 在小于第一近鄰范圍疊加了一項(xiàng)K(r1e—rij)n函數(shù)， 修正了原函數(shù)在近距離偏軟的現(xiàn)象， 其中K和n為勢參數(shù).

圖2 勢函數(shù)曲線 (a) 電子密度曲線; (b) 對勢曲線; (c) 有效對勢曲線Fig.2.Potential function curve: (a) Electron density curve;(b) potential curve; (c) effective pair potential curve.

我們知道， 從能量對距離的導(dǎo)數(shù)可以獲得力，從曲線形狀來看， 能量曲線的斜率越大， 排斥力越大.從數(shù)學(xué)上考慮， 對于函數(shù)(x0—x)n來說，x＜x0時(shí)， 隨著n的增加， 曲線的斜率將越來越大， 因此，通過幾項(xiàng)冪函數(shù)的迭加形成的多項(xiàng)式可以獲得比較滿意的曲線形式.由于近距離處的原子之間相互作用對間隙能的影響較大， 因此本文將間隙形成能作為目標(biāo)值， 通過調(diào)整多項(xiàng)式的冪次及系數(shù)， 最終得到(34)式的修正項(xiàng)， 使用獲得的加修正項(xiàng)的Nb勢函數(shù)計(jì)算的間隙形成能接近理論計(jì)算結(jié)果.

8 討 論

8.1 函數(shù)形式對勢函數(shù)性能的影響

從勢函數(shù)與晶體性能的關(guān)系表達(dá)式可以看出，勢函數(shù)值及勢函數(shù)的一階、二階導(dǎo)數(shù)值對材料性能的各個(gè)物理量的影響.如: 從(11) 式、(20)式可以看出對勢函數(shù)、電子密度函數(shù)在第一、二、三近鄰處的斜率對壓力P、平衡條件有影響; 從(13)式、(27)式—(30)式可以看出對勢函數(shù)、電子密度函數(shù)、以及多體相互作用項(xiàng)f(ρ)在第一、二、三近鄰處的斜率、曲率對體模量B及彈性常數(shù)有影響; 從(15) 式—(18)式及(21)式可以看出對勢函數(shù)、電子密度函數(shù)， 以及多體相互作用項(xiàng)f(ρ)在第一、二、三近鄰處的值對與能量相關(guān)的物理量， 如表面能、空位能等有影響.比較擬合所得勢函數(shù)計(jì)算的各物理量與實(shí)驗(yàn)值或理論值來修正勢函數(shù)曲線.如本文中根據(jù)間隙能的結(jié)果調(diào)整對勢曲線在近距離的行為， 進(jìn)而達(dá)到修正勢函數(shù)的目的.

本文也嘗試以修改勢函數(shù)曲線形式來考察擬合結(jié)果的變化.原始文獻(xiàn)[2]電子密度函數(shù)及對勢函數(shù)形式為

為方便討論， 幾種函數(shù)形式的截?cái)嗑嚯x均選取如下形式:使用不同形式的電子密度函數(shù)及對勢函數(shù)形式利用表2中的實(shí)驗(yàn)數(shù)據(jù)擬合獲得的勢函數(shù)參數(shù)見表6.

由表6可以看出， 前3組數(shù)據(jù)的擬合參數(shù)不同， 是因?yàn)楹瘮?shù)形式不同.后3組的勢參數(shù)完全相同， 是由于本文在擬合勢參數(shù)時(shí)使用的擬合數(shù)據(jù)a，及都是平衡態(tài)下的數(shù)據(jù)， 這些數(shù)據(jù)只與平衡態(tài)下曲線在第一、二、三近鄰處的值、斜率及曲率有關(guān)， 而后三組勢函數(shù)只對小于第一近鄰部分的曲線進(jìn)行修正， 這種修正并沒有改變曲線在第一、二、三近鄰處的值、斜率及曲率， 因此擬合獲得的勢函數(shù)參數(shù)沒有改變.也就是說， 當(dāng)使用平衡態(tài)下的數(shù)據(jù)作為擬合數(shù)據(jù)時(shí)， 調(diào)整勢函數(shù)曲線只要不改變曲線在各個(gè)近鄰處的值、斜率及曲率， 擬合的勢參數(shù)就不變.

使用獲得的各個(gè)勢函數(shù)計(jì)算表1中的晶格常數(shù)、表面能、空位能及體模量在均方差范圍內(nèi)與實(shí)驗(yàn)數(shù)據(jù)一致.除擬合數(shù)據(jù)之外的各物理量的計(jì)算結(jié)果見表7.

從表7可以看出， 第1列數(shù)據(jù)與表4的彈性常數(shù)及表5的間隙形成的計(jì)算值相差懸殊， 特別是C'出現(xiàn)了負(fù)值的情況， 表明原始的勢函數(shù)形式不適合于擴(kuò)展到第三近鄰.第2列數(shù)據(jù)修改了原始的對勢函數(shù)， 將四次多項(xiàng)式修改為六次多項(xiàng)式， 保持電子密度為原始函數(shù)形式， 結(jié)果有所改善， 但與實(shí)驗(yàn)值相比仍然相差很大.第3列為本文提出的函數(shù)形式， 將電子密度修改為四次冪的形式， 將對勢修改為六次多項(xiàng)式的形式， 彈性常數(shù)的結(jié)果有了明顯的改善， 但是間隙形成能的結(jié)果仍然遠(yuǎn)高于實(shí)驗(yàn)值.第4列和第5列數(shù)據(jù)為分別使用(33)式和(34)式對本文的對勢函數(shù)進(jìn)行了修正， 也就是對小于第一近鄰部分的對勢曲線進(jìn)行了軟化， 結(jié)果表明對彈性常數(shù)幾乎沒有影響， 原因與前述對勢參數(shù)的影響相同， 本文這種對對勢函數(shù)在近距離處的軟化處理并沒有改變對勢曲線在第一、二、三近鄰處的值、斜率及曲率， 而平衡態(tài)下的彈性常數(shù)只受平衡態(tài)下原子間勢函數(shù)在各個(gè)近鄰的值、斜率及曲率的影響，因此這種修正對彈性常數(shù)的計(jì)算結(jié)果沒有影響.第4列數(shù)據(jù)的間隙形成能與前3列數(shù)據(jù)相比明顯降低， 而第五列數(shù)據(jù)為根據(jù)間隙能的實(shí)驗(yàn)結(jié)果對本文的對勢函數(shù)進(jìn)行的進(jìn)一步修正的結(jié)果， 可以看出間隙形成能稍高于實(shí)驗(yàn)數(shù)據(jù)， 這是可以理解的， 因?yàn)楸疚挠?jì)算的間隙能是未弛豫的.表明本文的這種軟化修正影響間隙能的結(jié)果.這是由于間隙原子與近鄰原子的距離往往小于bcc結(jié)構(gòu)的第一近鄰， 而近距離排斥力大， 表明晶體比較硬， 這時(shí)向晶體中放入一個(gè)間隙原子將需要大的能量; 當(dāng)排斥力減小時(shí)， 表明晶體變軟， 這時(shí)向晶體中放入一個(gè)間隙原子所需要的能量將減小， 因此間隙能降低.

表6 不同函數(shù)形式的勢參數(shù)Table 6.Potential parameters of different functional forms.

表7 不同函數(shù)形式的各物理量計(jì)算結(jié)果Table 7.Calculation results of each physical quantity in different function forms.

8.2 截?cái)嗑嚯x對勢函數(shù)性能的影響

為了考察截?cái)嗑嚯x對勢函數(shù)性能的影響， 本文使用修正后的勢函數(shù)形式即(5)式， (32)式和(34)式， 計(jì)算不同截?cái)嗑嚯x下各物理量的變化.先保持電子密度函數(shù)的截?cái)嗑嚯x不變， 選擇不同的對勢截?cái)嗑嚯x， 分別擬合勢函數(shù)， 再利用所得的勢函數(shù)計(jì)算材料的各種物理量， 對勢函數(shù)的截?cái)嗑嚯x形式為:

當(dāng)x分別取0.55， 0.7， 0.80時(shí)， 各種物理量的計(jì)算結(jié)果如表8所示.

表8 不同對勢截?cái)嗑嚯x下的各物理量計(jì)算結(jié)果Table 8.Calculation results of each physical quantity under different pair potential cutoff distance.

表8中前4行數(shù)據(jù)為擬合數(shù)據(jù)， 在均方差范圍內(nèi)與實(shí)驗(yàn)結(jié)果一致.第5—9行為與力學(xué)性質(zhì)相關(guān)的數(shù)據(jù)， 幾種截?cái)嗑嚯x下， 結(jié)果差別不明顯.后6行數(shù)據(jù)為間隙能的計(jì)算結(jié)果， 可以看出， 間隙形成能變化不明顯， 但x= 0.8時(shí)的均方差最小.

隨后， 本文先保持對勢函數(shù)的截?cái)嗑嚯xc=不變， 選擇不同的電子密度函數(shù)的截?cái)嗑嚯x， 分別擬合勢函數(shù)， 再利用所得的勢函數(shù)計(jì)算材料的各種性能， 電子密度函數(shù)的截?cái)嗑嚯x形式為

當(dāng)y分別取0.45， 0.50， 0.60時(shí)， 各種物理量的計(jì)算結(jié)果如表9所示.

表9 不同電子密度截?cái)嗑嚯x下的各物理量計(jì)算結(jié)果Table 9.Calculation results of each physical quantity under different electron density cutoff distance.

從表9可以看出， 前4行數(shù)據(jù)為擬合數(shù)據(jù)， 在均方差范圍內(nèi)與實(shí)驗(yàn)結(jié)果一致， 其余彈性常數(shù)及間隙形成能的結(jié)果均變化不大， 在本文的情況下， 電子密度的截?cái)嗑嚯x對勢函數(shù)的性能影響不大.本文的截?cái)喾绞經(jīng)]有改變曲線在各個(gè)近鄰處的值、斜率及曲率， 并且本文沒有計(jì)算晶體被壓縮和拉伸等情況， 由于在各個(gè)近鄰原子范圍內(nèi)沒有出現(xiàn)原子數(shù)的變化， 因此截?cái)喾绞綄τ?jì)算結(jié)果幾乎沒有影響.但是如果采取其他方式對函數(shù)進(jìn)行截?cái)啵?如(9)式的形式， 則不同的截?cái)嗑嚯x對曲線的斜率和曲率有影響， 進(jìn)而對各物理量的結(jié)果將會(huì)產(chǎn)生影響.另外，如果計(jì)算拉伸等體積或形狀有變化的情況， 如果截?cái)嗑嚯x使相同近鄰原子內(nèi)的原子數(shù)發(fā)生變化時(shí)， 對計(jì)算結(jié)果也將產(chǎn)生影響， 為此進(jìn)行這類計(jì)算時(shí)， 截?cái)嗑嚯x的選取需要慎重.

9 結(jié) 論

本文從文獻(xiàn)[2]的FS勢出發(fā)， 將原子間的相互作用擴(kuò)展到第三近鄰，構(gòu)造了過渡金屬Nb的原子間勢，并較詳細(xì)地?cái)⑹隽碎_發(fā)原子間勢的方法.研究了勢函數(shù)曲線形式對晶體性質(zhì)的影響， 同時(shí)研究了函數(shù)形式、截?cái)嗑嚯x的選擇對材料性質(zhì)的描述及勢函數(shù)質(zhì)量的影響.得到如下結(jié)論.

1) 原始的FS勢函數(shù)形式不適合原子間相互作用擴(kuò)展到第三近鄰.經(jīng)過分析及嘗試發(fā)現(xiàn)， 修正電子密度函數(shù)為四次冪的形式， 對勢函數(shù)形式為六次多項(xiàng)式的形式時(shí)， 原子間勢能較好地描述原子間的相互作用.

2)以間隙形成能的結(jié)果作為目標(biāo)值修正了對勢函數(shù)在近距離的行為， 修正后的勢函數(shù)給出了接近DFT計(jì)算結(jié)果的間隙形成能.

3)當(dāng)使用平衡態(tài)下的物理量作為擬合數(shù)據(jù)時(shí)，調(diào)整勢函數(shù)曲線形式， 只要不改變函數(shù)曲線在各個(gè)近鄰處的函數(shù)值、函數(shù)斜率及曲率， 對擬合的勢參數(shù)沒有影響， 對于彈性常數(shù)的結(jié)果也沒有影響， 改變近距離處曲線的形狀影響間隙能的大小.

4)在本文的截?cái)喾绞较拢?改變截?cái)嗑嚯x對勢參數(shù)和晶體性能的計(jì)算結(jié)果均沒有太大影響.

附錄A 立方晶體剪切模量與原子間勢的關(guān)系

晶體中任意一點(diǎn)(xi，yi，zi)距原點(diǎn)的距離為ri=

對于C44， 晶體中任意一點(diǎn)(xi，yi，zi)的彈性位移為δxi=ε12yi， δyi=ε12xi， δzi=0， 由于應(yīng)變分量是小量， 故:

對于C′， 晶體中任意一點(diǎn)((xi，yi，zi)的彈性位移為δxi=ε11xi， δyi= — ε11yi， δzi= 0， 由于應(yīng)變分量是小量， 故

對于bcc結(jié)構(gòu)金屬Nb， 設(shè)晶格常數(shù)為a，Ωe=a3/2 ， 考慮前三層原子， 第一層有8個(gè)原子， 坐標(biāo)為距原點(diǎn)距離為第二層有6個(gè)原子， 坐標(biāo)為 (±a，0，0) ， ( 0，±a，0) ， ( 0，0，±a) ， 距原點(diǎn)距離為a; 第三層有12個(gè)原子， 坐標(biāo)為(±a，±a，0) ， (±a，0，±a) ， ( 0，±a，±a) ， 距原點(diǎn)距離為

先計(jì)算C44， 根據(jù)(25)式、(26)式及(1)式—(4)式可得

再計(jì)算C′， 同樣根據(jù)(25)式、(26)式及(1)式—(4)式可得

也可以采用另外一種方法計(jì)算彈性常數(shù)C11，C12和C'.

晶體的彈性變形通過彈性常數(shù)矩陣描述，

其中，E是晶體內(nèi)能，εij是應(yīng)變分量.

在絕對零度下， 晶體內(nèi)能即是原子間相互作用勢能.

有研究表明， 在體系處于平衡態(tài)時(shí)， 依據(jù)彈性理論可以得到原子間勢與彈性常數(shù)之間的關(guān)系[25]為

其中，

(A1)式利用了平衡條件

其中，

利用Voigt標(biāo)記法， 可得立方晶體的彈性常數(shù)C11，C12為:

考慮原子分布情況， 對于金屬Nb， 相互作用在第三和第四近鄰之間截?cái)啵?彈性常數(shù)C11，C12涉及的各項(xiàng)結(jié)果如下:

可得

由(A2)式可得