王永偉,郭永亮,,張偉,柯學(xué)志

(1.華東師范大學(xué)物理與材料科學(xué)學(xué)院,上海200241; 2.中國科學(xué)院上海應(yīng)用物理研究所,上海201800)

基于第一原理計(jì)算擬合Ag、Si和C原子間作用勢的研究

王永偉1,郭永亮1,2,張偉2,柯學(xué)志1

(1.華東師范大學(xué)物理與材料科學(xué)學(xué)院,上海200241; 2.中國科學(xué)院上海應(yīng)用物理研究所,上海201800)

為了實(shí)現(xiàn)Ag在SiC晶體中擴(kuò)散的分子動力學(xué)模擬,根據(jù)第一性原理計(jì)算結(jié)果,采用“力匹配”方法對Ag、Si和C原子間相互作用勢進(jìn)行擬合,并用晶格常數(shù)、內(nèi)聚能、體彈性模量、彈性常數(shù)和缺陷形成能等進(jìn)行了勢函數(shù)的驗(yàn)證.結(jié)果表明,擬合所得的勢對Si、C和SiC晶體的內(nèi)聚能、晶格常數(shù)和體彈性模量的計(jì)算非常準(zhǔn)確,最大誤差不超過0.6%;并且,該勢對Si和C晶體中空位與間隙形成能及SiC晶體中Si與C空位形成能的計(jì)算值比采用J.Tersoff給出的勢的計(jì)算值更精確;此外,該勢對16種AgSiC三原子體系缺陷形成能的計(jì)算也比較精確.

Tersoff勢;力匹配;SiC;第一性原理計(jì)算

0 引言

彌散在石墨集體中的TRISO(Tristructural-Isotropic)包覆顆粒燃料是高溫氣冷反應(yīng)堆的基本元件,從內(nèi)到外它的基本組成單元依次是:燃料內(nèi)核(通常是UO2)、多孔石墨緩沖層、內(nèi)致密熱解炭層、SiC層、外熱解炭層,如圖1所示.在TRISO顆粒中SiC層為其提供了結(jié)構(gòu)支持并且起到阻礙主要核裂變產(chǎn)物擴(kuò)散的作用[1].雖然SiC層能夠有效抑制大部分金屬裂變產(chǎn)物的擴(kuò)散,但是也有實(shí)驗(yàn)發(fā)現(xiàn),在高溫條件下仍然有Ag等裂變產(chǎn)物透過SiC層擴(kuò)散到反應(yīng)堆冷卻劑中[1-2].這會導(dǎo)致Ag在反應(yīng)堆中的積淀,這給核電系統(tǒng)的維護(hù)工作造成了很大的困難,同時還存在一定的安全隱患,因?yàn)榉派湫援a(chǎn)物110Ag的半衰期比較長,為253 d.

圖1 TRISO包覆顆粒燃料結(jié)構(gòu)示意圖Fig.1 The structure of Tristructural-Isotropic-coated fuel particle

為了控制裂變產(chǎn)物Ag在SiC層中的擴(kuò)散,明確其擴(kuò)散機(jī)制是至關(guān)重要的.針對Ag在SiC層中的擴(kuò)散現(xiàn)象,很多擴(kuò)散機(jī)制被提出,其中包括SiC的降解[3-4]、通過裂縫和納米孔的擴(kuò)散[5-7]、沿著晶粒邊界和錯位層的擴(kuò)散[8-9]等.但是時至今日,究竟是哪一種擴(kuò)散機(jī)制起作用仍然是不明確的,有待于實(shí)驗(yàn)和理論工作的進(jìn)一步驗(yàn)證.

關(guān)于Ag在SiC中的擴(kuò)散,在實(shí)驗(yàn)[10-14]和基于密度泛函理論(Density Functional Theory, DFT)的第一性原理[15-16]方面都有相應(yīng)的研究,其中對Ag擴(kuò)散系數(shù)的計(jì)算都是根據(jù)計(jì)算出來的主要幾種點(diǎn)缺陷的形成能和遷移勢壘按照公式D=D0exp(?Q/kbT)估算獲得;然而沒有分子動力學(xué)模擬方面的研究,這是因?yàn)槿狈ο鄳?yīng)的原子間相互作用勢.眾所周知,由于計(jì)算條件的限制,第一性原理方法只能在小尺度范圍內(nèi)進(jìn)行計(jì)算研究;而分子動力學(xué)方法可以在大尺度下進(jìn)行長時間的動態(tài)模擬,并且可以通過分子動力學(xué)模擬對擴(kuò)散機(jī)制的微觀過程進(jìn)行直接觀察.分子動力學(xué)方法可以用均方位移法或速度關(guān)聯(lián)函數(shù)法直接計(jì)算擴(kuò)散系數(shù).對于缺陷遷移這類微觀上低概率事件,可以借助于加速分子動力學(xué)方法[17]來進(jìn)行模擬.

進(jìn)行分子動力學(xué)模擬研究需要相應(yīng)的勢函數(shù)作為支持,在AgSiC三原子體系中存在Si-Si、C-C、Si-C、Ag-Ag、Ag-Si和Ag-C等6種相互作用關(guān)系.但是,現(xiàn)有的勢中沒有Ag-Si和Ag-C勢.況且,已有的Si-Si勢、C-C勢、Si-C勢[18-22]和Ag-Ag勢[23]并不是針對我們研究的體系擬合的,需用進(jìn)一步修正甚至重新擬合.因此,需要建立一個AgSiC三原子體系的相互作用勢,為分子動力學(xué)方法研究Ag在SiC中擴(kuò)散提供支持.

1 理論方法

1.1 計(jì)算原理及公式

本文中勢函數(shù)的驗(yàn)證過程中所涉及的體系內(nèi)聚能、缺陷形成能和體彈性模量等的計(jì)算分別采用如下公式.

內(nèi)聚能的計(jì)算公式[24]為

其中,Ebulk和Eatomi分別表示單個化學(xué)勢晶體能量和單個原子能量(其中包括自旋能和零點(diǎn)振動能兩部分).

缺陷形成能的計(jì)算公式[15,17]為

其中,Edef表示包含缺陷的晶胞能量,Eundef表示完整晶胞能量,I表示元素種類,?nI表示完整晶胞中I類原子數(shù)減去缺陷晶胞中I類原子數(shù),μI表示對應(yīng)化合物構(gòu)型(SiC晶體)下I類原子的化學(xué)勢,EI是相應(yīng)單質(zhì)構(gòu)型(在這里就是Ag、Si和C晶體的)下單個I原子的能量,γI表示元素I在相應(yīng)單質(zhì)構(gòu)型下的化學(xué)勢,可根據(jù)文獻(xiàn)[25]確定它的取值范圍,即, 0≤γI≤ESiC?ESi?EC.這里分3種情況分別取[15]:γSi=γC=(ESiC?ESi?EC)/2; γSi=ESiC?ESi?EC,γC=0;γSi=0,γC=ESiC?ESi?EC.本文計(jì)算SiC晶體缺陷形成能時只考慮γC=(ESiC?ESi?EC[15])/2的情況.

晶格常數(shù)a、內(nèi)聚能E、體彈性模量B0等值是通過Murnaghan方程[26]對體積–能量關(guān)系的擬合結(jié)果獲得,Murnaghan方程的一般表達(dá)式為

本文采用Tersoff勢[21]來描述原子間相互作用,其一般形式為

其中,排斥(Repulsive)勢fR和吸引(Attractive)勢fA如式(5)所示,截?cái)?Cutoff)函數(shù)如式(6)所示,即

以及表示原子周圍環(huán)境的三體作用項(xiàng)

式(4)至式(7)中i,j,k表示原子序號,rik表示原子i和k之間的距離,θijk表示ij和ik原子間化學(xué)鍵之間的角度.為了簡單起見,令χ=1,則以上各式共有11個待定參數(shù).其中A,B,λ和μ這4個參數(shù)為兩體相互作用部分,β,n,c,d和h為三體相互作用部分,R,S的取值跟截?cái)喟霃接嘘P(guān).

在擬合過程中POTFIT采用最小二乘法、共軛梯度法、模擬退火法等,來優(yōu)化參考構(gòu)型的力和能量等物理量的第一性原理計(jì)算值和勢函數(shù)計(jì)算值之間的差值函數(shù),使之取最小值,從而獲得勢函數(shù)各參量的值.差值函數(shù)的一般形式為

其中,

分別代表來自內(nèi)聚能、力和壓強(qiáng)的貢獻(xiàn).在以上各式中Wi表示權(quán)重,NC表示所有構(gòu)型的數(shù)目,NA表示所有原子的數(shù)目,Fi,xj表示i原子在xj方向受到的力,σi,j表示i原子應(yīng)力張量的第j個數(shù)值.本文擬合時只考慮了式(8)中的能量和力兩項(xiàng)對差值函數(shù)的影響.

1.2 計(jì)算參數(shù)設(shè)置

本文采用Tersoff勢作為Ag、Si和C原子間作用的力場函數(shù),其合理性已經(jīng)在參考文獻(xiàn)[18-22,27]中得到了很好驗(yàn)證.擬合勢函數(shù)所需的第一性分子動力學(xué)、靜態(tài)優(yōu)化和電子步優(yōu)化等計(jì)算及驗(yàn)證勢函數(shù)所需的彈性常數(shù)、缺陷形成能、晶格常數(shù)和內(nèi)聚能等第一性原理的計(jì)算采用基于密度泛函理論(DFT)的VASP[28]軟件包.勢函數(shù)的擬合采用基于“力匹配”方法的POTFIT[29]軟件包.勢函數(shù)驗(yàn)證過程中用擬合的勢計(jì)算彈性常數(shù)、晶格常數(shù)、內(nèi)聚能、體彈性模量和SiC晶體缺陷形成能等值時采用基于“牛頓力學(xué)”的Lammps軟件包.

本文用VASP軟件包的投影綴加波方法進(jìn)行密度泛函計(jì)算,采用PBE[30]形式的GGA[31]交換關(guān)聯(lián)泛函求解Kohn-Sham方程.平面波截?cái)嗄苋?70 eV.Monkhorst-Pack K-point[32]網(wǎng)格:8原子體系取13×13×13;大于60原子體系取7×7×7.在SiC和AgSiC缺陷的第一性分子動力學(xué)計(jì)算中,步長取1.55 fs,溫度參考文獻(xiàn)[10-14]中的實(shí)驗(yàn)溫度取1 073～1 773 K.

Lammps計(jì)算過程中涉及的靜態(tài)優(yōu)化過程都采用最陡下降算法,步長取1.55 fs.POTFIT優(yōu)化過程中,體系能量權(quán)重取100,退火溫度取1 500 K.

2 結(jié)果與討論

2.1 構(gòu)型選取

在勢函數(shù)的擬合過程中要選取包含所研究的物理過程或物理?xiàng)l件的晶體結(jié)構(gòu)作為第一性原理計(jì)算的構(gòu)型.為了描述Ag在SiC晶體中的擴(kuò)散,AgSiC三原子缺陷體系構(gòu)型的選取David Shrader在文獻(xiàn)[15]中提到的16種缺陷作為基本構(gòu)型,它們分別是AgC、Ag TC、AgCVC、AgCVSi、AgC2VSi、AgCVSiVC、AgCVSiSiC、AgCVSiCSi、AgSi、Ag TSi、AgSiVSi、AgSiVC、AgSi2VC、AgSiVCVSi、AgSiVCSiC和AgSiVCCSi.在這里AgC和AgSi分別表示Ag原子替換了C和Si原子,VSi和VC分別表示在Ag原子最近的Si空位和C空位,它們前面的數(shù)字表示空位的數(shù)目;SiC和CSi則分別表示1個Si原子替換了離Ag原子最近的C原子和C原子替換了離Ag原子最近的Si原子, Ag TC和Ag TSi則分別表示Ag填隙在4個C的中心和4個Si的中心.為了便于理解,我們畫出了其中5種缺陷構(gòu)型,如圖2所示,其中,圖2(a)AgSiVC表示Ag原子取代了Si原子的位置,并且在Ag原子的4個最近鄰位置之一出現(xiàn)了1個C空位;圖2(b)AgSiVCVSi則表示取代Si的Ag原子的四個最近鄰位置上出現(xiàn)1個C空位,同時在Ag次近鄰位置上出現(xiàn)了一個Si空位;圖2(c)AgSi2VC表示取代Si的Ag原子的4個最近鄰位置上出現(xiàn)2個C空位;圖2(d)AgSiVCCSi表示取代Si的Ag最近鄰位置上出現(xiàn)1個C空位,同時次近鄰位置的Si被C取代;圖2(e)AgSiVCSiC則表示取代Si的Ag最近鄰位置上出現(xiàn)1個C空位,同時最近鄰位置的另一個C被Si原子取代.其他11種構(gòu)型表示方法以此類推.

圖2 5種AgSiC三原子體系缺陷示意圖Fig.2 Five kinds of AgSiC three-atom system defects

本文研究的AgSiC三原子體系的主體是SiC材料,所以在勢的擬合過程中要兼顧SiC晶體的一些性質(zhì).因此我們首先對Si-Si,C-C和Si-C原子間作用勢進(jìn)行了擬合.在本文中Si-Si、C-C和Si-C原子間Tersoff勢的截?cái)喟霃絒21]分別取3.0?A,2.1?A和2.5?A;而在完全優(yōu)化后的SiC晶格中Si-Si,C-C和Si-C原子間距分別約為3.1?A、3.1?A和1.9?A.Si-Si和C-C間距都為3.1?A大于對應(yīng)的勢函數(shù)截?cái)喟霃?.0?A和2.1?A,因此在計(jì)算完美SiC晶體性質(zhì)時, Si-Si和C-C勢對其結(jié)果幾乎沒有影響;而Si-C間距1.9?A在Si-C勢的截?cái)喟霃?.5?A之內(nèi),因此Si-C作用在SiC晶體性質(zhì)計(jì)算時有很大影響.我們用文獻(xiàn)[21]中給出的Tersoff勢進(jìn)行了驗(yàn)證,結(jié)果表明,只有在SiC晶體中存在間隙等缺陷時,Si-Si和C-C作用勢才會對SiC晶體性質(zhì)的計(jì)算結(jié)果產(chǎn)生影響,否則沒有影響.這是因?yàn)殚g隙的存在使得Si-Si和C-C原子間距分別在3.0?A和2.1?A以下產(chǎn)生了分布.因此,我們認(rèn)為Si-Si和C-C作用勢對SiC晶體性質(zhì)直接產(chǎn)生的影響很小,可以最先單獨(dú)進(jìn)行擬合,這樣使得擬合過程快速和簡單,并且得到的SiC體系的勢還可以兼顧C(jī)和Si晶體各自的性質(zhì).

基于以上分析,具體擬合過程分3步進(jìn)行:首先,分別擬合Si-Si和C-C勢;然后,固定Si-Si和C-C勢的各參數(shù)不變進(jìn)一步擬合Si-C勢;最后,固定前兩步中的擬合結(jié)果不變,進(jìn)行Ag-Si和Ag-C作用勢的擬合.相應(yīng)構(gòu)型的選取也分為如下3類.

1.Si-Si和C-C作用勢構(gòu)型

(1)首先,完全優(yōu)化64原子Si(C)晶格;固定晶格常數(shù)分別對65(間隙)、63(空位)原子的Si(C)構(gòu)型進(jìn)行等體積優(yōu)化,并計(jì)算相應(yīng)構(gòu)型Si(C)原子間距的分布及所對應(yīng)權(quán)重(這里的權(quán)重與相同間距出現(xiàn)的次數(shù)成正比).Si和C單質(zhì)體系的間隙、空位結(jié)構(gòu)如圖3所示.

(2)根據(jù)(1)中原子間距分布,放大或縮小8原子完全優(yōu)化的晶格,使其原子間距剛好滿足(1)中得到的原子間距分布.一般會產(chǎn)生10至20個構(gòu)型.

圖3 Si(C)晶體中空位和間隙示意圖Fig.3 Vacancy and interstitial of Si(C)crystal

2.Si-C作用勢構(gòu)型

(1)狀態(tài)方程.為了確定SiC晶體的晶格常數(shù)、形成能、體彈性模量和彈性常數(shù),我們對SiC完整晶格進(jìn)行了完全優(yōu)化(平衡體積為V0),并選取了0.90V0至1.15V0的15個構(gòu)型.

(2)缺陷.為了初步考慮AgSiC體系中的缺陷,我們對AgC等16種缺陷固定體積進(jìn)行了優(yōu)化(采用64原子SiC晶體完全優(yōu)化的晶格常數(shù)),然后去掉缺陷構(gòu)型中Ag原子并保持原有Si原子、C原子位置不變,進(jìn)行DFT靜態(tài)計(jì)算(只做電子步計(jì)算).缺陷結(jié)構(gòu)參照圖2.此外,我們還考慮了SiC中的本征空位VC和VSi兩種構(gòu)型,其結(jié)構(gòu)如圖4所示.

圖4 SiC晶體中Si和C空位示意圖Fig.4 Si and C vacancy of SiC crystal

(3)分子動力學(xué)過程.為了獲得SiC晶體在相應(yīng)溫度下可能出現(xiàn)的構(gòu)型及其力和能量等信息,我們?nèi)?4原子的SiC超胞,采用NVT系綜在1 073～1 773 K下進(jìn)行第一性原理分子動力學(xué)計(jì)算,每個分子動力學(xué)過程進(jìn)行500步,步長1.55 fs,從軌跡中選取20個不同溫度下的構(gòu)型.

3.Si-Ag和C-Ag作用勢構(gòu)型

(1)缺陷靜態(tài)優(yōu)化.對AgC等16種缺陷進(jìn)行固定體積優(yōu)化,其晶格常數(shù)及體積取完全優(yōu)化的SiC晶體的值.

(2)缺陷分子動力學(xué)過程.為了獲得各缺陷體系在相應(yīng)溫度下可能出現(xiàn)的構(gòu)型及其力和能量等信息,我們用第一性原理分子動力學(xué)方法對AgC等16種缺陷在NVT系綜下運(yùn)行500步,步長1.55 fs,溫度1 073～1 773 K.每種缺陷在不同溫度下選取10個構(gòu)型.

2.2 擬合過程

擬合過程中,首先,根據(jù)第2.1節(jié)中所述的3類構(gòu)型選取構(gòu)型,并用第一性原理計(jì)算各構(gòu)型的力和能量等信息.然后,用POTFIT最優(yōu)化Tersoff勢函數(shù)中的11個待定參數(shù),調(diào)整各參數(shù)的變化范圍直到結(jié)果收斂.最后,用所得的Tersoff勢計(jì)算各體系晶格常數(shù)a、內(nèi)聚能E、彈性常數(shù)Cij、體彈性模量B0、缺陷形成能和16種缺陷的內(nèi)聚能,并跟第一性原理計(jì)算結(jié)果進(jìn)行對比,根據(jù)對比結(jié)果調(diào)整各個構(gòu)型的權(quán)重,反復(fù)擬合,直到對比結(jié)果相吻合.

在Si-Si和C-C勢的擬合過程中,我們只對B、λ和μ這3個參數(shù)進(jìn)行了改動,其他參數(shù)借鑒文獻(xiàn)[21]中的數(shù)值保持不變.之所以這么做是由于本文擬合基于第一性原理計(jì)算,而文獻(xiàn)中[21]中的擬合基于實(shí)驗(yàn)數(shù)據(jù),兩者能量不匹配,而A、B、λ和μ等參數(shù)主要確定體系能量的高低,這里B、λ和μ參數(shù)足以達(dá)到能量匹配的效果;R和S參數(shù)確定相互作用勢的截?cái)喾秶?β、n、c、d和h跟原子周圍環(huán)境因素有關(guān),都可以不用改變.

在擬合Si-C勢的擬合過程中,除了表示截?cái)嗟腞和S這2個參數(shù)與文獻(xiàn)[21]保持一致外,其他參數(shù)都進(jìn)行了改變.這是由于除了能量不匹配外,本文考慮的缺陷與文獻(xiàn)[21]中考慮的缺陷不完全同,從而部分原子的周圍環(huán)境因素不完全相同,跟周圍環(huán)境因素有關(guān)的參數(shù)需用重新擬合.

在Ag-Si和Ag-C作用勢的擬合過程中,β和n參數(shù)的選取參考了文獻(xiàn)[33]中的結(jié)果直接取1,R和S的取值結(jié)合SiC晶體中原子最近鄰和次近鄰間距進(jìn)行調(diào)整.

2.3勢的驗(yàn)證及結(jié)果

由第一性原理計(jì)算的Si、C和SiC晶體的晶格常數(shù)、內(nèi)聚能、體彈性模量、彈性常數(shù)、缺陷形成能等物理量見表1 DFT所示.這些數(shù)據(jù)作為勢函數(shù)驗(yàn)證的主要依據(jù),其準(zhǔn)確性決定驗(yàn)證環(huán)節(jié)的可靠程度.因此我們拿這些數(shù)據(jù)跟其他第一性原理計(jì)算結(jié)果(表1 Refa和Refb)進(jìn)行了對比,由表1 DFT、Refa和Refb的數(shù)據(jù)可以看出,本文的第一性原理計(jì)算值跟其他第一性理論計(jì)算值非常接近,其誤差跟計(jì)算條件及精度有關(guān),可以接受.

圖5中(a)、(b)和(c)分別給出了Si、C和SiC晶體的內(nèi)聚能-體積關(guān)系及用Murnaghan方程擬合的結(jié)果.圖中,縱坐標(biāo)為平均每個原子的內(nèi)聚能,橫坐標(biāo)為平均每個原子的體積,黑色“△”、藍(lán)色“□”和紅色“○”分別表示文獻(xiàn)[21]中所給出的勢(Tersoff)、第一性原理(DFT)和本文擬合的勢(This work)的計(jì)算值,相應(yīng)的同種顏色曲線代表Murnaghan方程擬合結(jié)果.由圖5可以看出,由第一性原理數(shù)據(jù)的Murnaghan曲線(藍(lán)色)和本文給出的勢的數(shù)據(jù)的Murnaghan曲線(黑色)基本重合,兩者符合得很好.由Murnaghan方程擬合得到的晶格常數(shù)、內(nèi)聚能和體彈性模量的值見表1所示.

圖5 Si、C和SiC晶體在不同體積下內(nèi)聚能的計(jì)算值及用Murnaghan方程擬合的結(jié)果Fig.5 The calculated values of cohesive energies at different volumes in Si,C and SiC crystals and the Murnaghan equation fitting results of these values

表1給出了用本文擬合得到的勢計(jì)算獲得的Si、C和SiC晶體中晶格常數(shù)a、內(nèi)聚能E、體彈性模量B0、彈性常數(shù)Cij、空位形成能Vac、間隙形成能Int、Si空位形成能VSi和C空位形成能VC等物理量數(shù)值,并跟第一性原理計(jì)算值、實(shí)驗(yàn)值及別人第一性原理計(jì)算值等對比.其中,exp、Tersoff和Ref分別表示實(shí)驗(yàn)值、文獻(xiàn)[21]中勢函數(shù)的計(jì)算值和參考文獻(xiàn)的第一性原理計(jì)算值,相應(yīng)的數(shù)據(jù)來源已在表格中標(biāo)出;This work和DFT分別表示本文擬合所得勢的計(jì)算值和第一性原理計(jì)算值;空位Vac、間隙Int、Si空位VSi和C空位VC的結(jié)構(gòu)如圖3和圖4所示.由表1(a)和表1(b)可以看出,本文擬合的勢對Si和C晶體的晶格常數(shù)、內(nèi)聚能、彈性常數(shù)和體彈性模量的計(jì)算值跟第一性原理計(jì)算值非常吻合;而對Si和C晶體中空位和間隙形成能的計(jì)算值跟文獻(xiàn)[21]給出勢的計(jì)算值比較接近,但與其相比更接近第一性原理計(jì)算值,其結(jié)果有所改善;這是由于我們在Si-Si和C-C勢的擬合過程中只擬合了B、λ和μ這3個參數(shù),同時在構(gòu)型選取時考慮了缺陷體系的原子間距分布的緣故.由表1(c)SiC可以看出,本文擬合的勢對其他物理量的計(jì)算值跟第一性原理計(jì)算值吻合得很好.

表1 Si、C和SiC晶體中各物理量的計(jì)算值Tab.1 The calculated values of various physical quantities in Si,C and SiC crystals

表2給出了計(jì)算AgSiC缺陷形成能所需的化學(xué)勢的值,并與文獻(xiàn)[15]中的數(shù)據(jù)進(jìn)行了對比,其計(jì)算方法及物理意義已在第1.1節(jié)中介紹.由第1.1節(jié)中的介紹可知,本文和文獻(xiàn)[15]中給出的各化學(xué)勢的差異歸根結(jié)底是ESi、EC和ESiC值的不同.文獻(xiàn)[15]中取塊體能量–5.44 eV、–9.20 eV和–15.08 eV作為它們的值,本文取內(nèi)聚能–4.545 eV、–7.773 eV和–12.862 eV作為它們的值.考慮Si和C的原子能–0.880 eV和–1.316 eV后,根據(jù)式(1)可由內(nèi)聚能計(jì)算塊體能,其值分別為–5.425 eV、–9.089 eV和–15.059 eV.這跟文獻(xiàn)[15]中給出的值仍然有微小差別,這是由于計(jì)算精度及方法不同所致.

表2 Ag、Si和C的化學(xué)勢Tab.2 Chemical potentials for Ag、Si and C

表3給出了本文擬合的勢對16種AgSiC三原子體系缺陷形成能的計(jì)算值跟第一性原理計(jì)算值及文獻(xiàn)[15]中給出的值的對比結(jié)果及相對誤差.可以看出,本文擬合的勢計(jì)算16種AgSiC體系缺陷形成能的值與第一性原理計(jì)算的值符合得很好.除了AgCVSiVC、AgCVSiSiC、AgCVC和AgSiVCSiC4這種缺陷形成能誤差分別為2.12 eV、2.31 eV、4.83 eV和2.31 eV外,其他缺陷形成能的誤差都不超過0.84 eV.由于我們考慮的缺陷種類較多,缺陷結(jié)構(gòu)差異較大,擬合得到的勢很難顧及到每一個缺陷的性質(zhì),這是經(jīng)驗(yàn)勢本身的局限性所致.個別構(gòu)型的缺陷形成能顯得大,但這也是可以接受的.如參考文獻(xiàn)[20-21]中的勢對int(S)和VSi+VC缺陷形成的計(jì)算值跟第一性原理計(jì)算值的最大誤差分別為6.7 eV和5.3 eV,相對誤差分別為40%和42%.

表3 16種缺陷體系缺陷形成能的計(jì)算值Tab.3 The calculated formation energies of 16 kinds of defect systems

最終得到的Si、C和Ag原子間相互作用勢參數(shù)如表4所示,本文所研究的缺陷體系中只有一個Ag原子,故暫不考慮Ag原子間相互作用,Ag-Ag勢參數(shù)手動設(shè)置使之為零.我們用表4中的勢分別計(jì)算了Si、C和SiC的多個物理量并跟實(shí)驗(yàn)值(exp)進(jìn)行了對比,如表1所示,兩者有一定的差別,這是由于我們的擬合是基于第一性原理計(jì)算進(jìn)行的.這個誤差的根源是第一性原理計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果之間的差別.總體上來說我們的擬合是成功的.

表4 Si、C和Ag原子間相互作用勢參數(shù)表Tab.4 Interatomic potential parameters of Si,C and Ag

3 結(jié)論

本文采用“力匹配”方法,基于第一性原理計(jì)算擬合了存在輻照缺陷的SiC晶體結(jié)構(gòu)中Ag、Si和C原子間相互作用勢.這種方法避免了傳統(tǒng)的基于實(shí)驗(yàn)數(shù)據(jù)擬合勢的方法中數(shù)據(jù)源不足的問題,為復(fù)雜體系材料原子間作用勢的擬合提供了可能.擬合結(jié)果表明,通過該方法獲得的原子間相互作用勢能夠很好地描述Si、C和SiC晶體材料的晶格常數(shù)、內(nèi)聚能、彈性模量和缺陷形成能等.采用Tersoff形式的勢來描述Ag摻雜的SiC晶體結(jié)構(gòu)中Ag、Si和C原子間相互作用行之有效,不但可以兼顧Si、C和SiC晶體各自的性質(zhì),而且還能準(zhǔn)確計(jì)算AgSiC體系的缺陷形成能.當(dāng)然,由于經(jīng)驗(yàn)勢本身的局限性,我們對缺陷形成能的擬合中,個別構(gòu)型的誤差還是較大,這方面還有一定的提升空間.為擬合更為精確的勢,需要更加全面地去考慮各個因素的影響,我們的現(xiàn)有結(jié)果為進(jìn)一步的研究奠定了初步基礎(chǔ).

[1]RUBIN S D.TRISO-coated particle fuel phenomenon identifi cation and ranking tables(PIRTs)for fi ssion product transport due to manufacturing,operations and accidents[R].USA:US-NRC,2004.

[2]VERFONDERN K.Fuel performance and fi ssion product behavior in gas-cooled reactors No.TECDOC-978[R]. Vienna:IAEA,1997.

[3]MINATO K,SAWA K,KOYA T,et al.Fission product release behavior of individual coated fuel particles for high-temperature gas-cooled reactors[J].Nucl Technol,2000,131:36-47.

[4]SCHENK W,POTT G,NABIELEK H.Fuel accident performance testing for small HTRs[J].J Nucl Mater, 1990,175:19-30.

[5]MINATO K,OGAWA T,FUKUDA K,et al.Release behavior of metallic fi ssion products from HTGR fuel particles at 1 600 to 1 900°C[J].J Nucl Mater,1993,202:47-53.

[6]FRIEDLAND E,MALHERBE J B,VANDERBERG N G,et al.Study of silver diff usion in silicon carbide[J].J Nucl Mater,2009,389:326-331.

[7]MACLEAN H,BALLINGER R,KOLAYA L,et al.The eff ect of annealing at 1500?C on migration and release of ion implanted silver in CVD silicon carbide[J].J Nucl Mater,2006,357:31-47.

[8]BULLOCK R E.Fission-product release during postirradiation annealing of several types of coated fuel particles[J].J Nucl Mater,1984,125:304-319.

[9]PETTI D,BUONGIORNO J,MAKI J,et al.Key diff erences in the fabrication,irradiation and high temperature accident testing of US and German TRISO-coated particle fuel,and their implications on fuel performance[J]. Nucl Eng Des,2003,222:281-297.

[10]NABIELEK H,BROWN P E,OFFERMAN P.Silver release from coated particle fuel[J].Nucl Technol,1977, 35:483-493.

[11]VERFONDERN K,MARTIN R C,MOORMANNN R.Methods and data for HTGR fuel performance and radionuclide release modeling during normal operation and accidents for safety analyses No.JUEL-2721[R]. Germany:Forschungszentrum J¨ulich GmbH,1993.

[12]AMIAN W,STOVER D.Diff usion of silver and cesium in silicon-carbide coatings of fuel particles for hightemperature gas-cooled reactors[J].Nucl Technol,1983,61:475-486.

[13]FRIEDLAND E,MALHERBE J B,VANDERBERG N G,et al.Study of silver diff usion in silicon carbide[J]. J Nucl Mater,2009,389:326-331.

[14]MACLEAN H J.Silver transport in CVD silicon carbide[D].Boston:MIT,2004.

[15]SHRADER D,KHALIL S,GERCZAK T,et al.Ag diff usion in cubic silicon carbide[J].J Nucl Mater,2010,408: 257-271.

[16]KHALIL S,SWAMINATHAN N,SHRADER D,et al.Diff usion of Ag alongΣ3 grain boundaries in 3C-SiC[J]. Phys Rev B,2011,84:214104.

[17]VOTER A F.Hyperdynamics:Accelerated molecular dynamics of infrequent events[J].Phys Rev Lett,1997,78: 3908-3911.

[18]TERSOFF J.New empirical approach for the structure and energy of covalent systems[J].Phys Rev B,1988, 37:6991-7000.

[19]TERSOFF J.New empirical model for the structural properties of silicon[J].Phys Rev Lett,1986,56:632-635.

[20]TERSOFF J.Empirical interatomic potential for carbon,with applications to amorphous carbon[J].Phys Rev Lett,1988,61:2879-2882.

[21]TERSOFF J.Modeling solid-state chemistry:Interatomic potentials for multicomponent systems[J].Phys Rev B,1989,39:5566-5568.

[22]TERSOFF J.Empirical interatomic potential for silicon with improved elastic properties[J].Phys Rev B,1988, 38:9902-9905.

[23]FOILES S M,BASKES M I,DAW M S.Embedded-atom-method functions for the fcc metals Cu,Ag,Au,Ni, Pd,Pt,and their alloys[J].Phys Rev B,1986,33:7983-7991.

[24]LI X P,CEPERLEY D M,MARTIN R M.Cohesive energy of silicon by the Green’s-function monte carlo method[J].Phys Rev B,1991,44:10929-10932.

[25]KOHAN A F,CEDER G,MORGAN D,et al.First-principles study of native point defects in ZnO[J].Phys Rev B,2000,61:15019-15027.

[26]MURNAGHAN F D.The compressibility of media under extreme pressures[J].Proceeding of the National Academy of Sciences of the United States of America,1944,30(9):244-247.

[27]BUTLER K T,VULLUM P E,MUGGERUD A M,et al.Structural and electronic properties of silver/silicon interfaces and implications for solar cell performance[J].Phys Rev B,2011,83(23):2155-2161.

[28]KRESSE G,FURTHMULLER J.Effi cient iterative schemes for ab initio total-energy calculations using a planewave basis set[J].Phys Rev B,1996,54(16):11169-11186.

[29]BROMMER P,G¨AHLER F.Potfi t:Eff ective potentials from ab-initio data[J].Simul Mater Sci Eng,2007,15: 295-304.

[30]PERDEW J P,BURKE K,ERNZERHOF M.Generalized gradient approximation made simple[J].Phys Rev Lett,1996,77:3865-3868.

[31]MARTIN R M.Electronic Structure:Basic Theory and Practical Methods[M].Cambridge:Cambridge University Press,2004.

[32]MONKHORST H J,PACK J D.Special points for Brillouin-zone integrations[J].Phys Rev B,1976,13:5188-5200.

[33]ERHART P,ALBE K.Analytical potential for atomistic simulations of silicon,carbon,and silicon carbide[J]. Phys Rev B,2005,71(3):035211.

[34]MOORE C E.Atomic Energy Levels Volumer[M].Washington D C:NBS,1949.

[35]AADERSON O L.The use of ultrasonic measurements under modest pressure to estimate compression at high pressure[J].J Phys Chem Solids,1966,27:547-565.

[36]DONOHUE J.The structures of the elements[J].Diamond and Related Materials,1974,24(4):436.DOI: 10.1016/j.diamond.2011.01.035.

[37]YIN M T,COHEN M L.Microscopic theory of the phase transformation and lattice dynamics of Si[J].Phys Rev Lett,1980,45:1004-1007.

[38]CAR R,KELLY P J,OSHIYAMA A,et al.Microscopic theory of atomic diff usion mechanisms in silicon[J]. Phys Rev Lett,1984,52:1814-1817.

[39]BARAFF G A,SCHLUTER M.Migration of interstitials in silicon[J].Phys Rev B,1984,30:3460-3469.

[40]BREWER L.Lawrence berkeley laboratory report No.LB-3720[R].California:Lawrence Berkeley Laboratory, 1977.

[41]MCSKIMIN H J,ANDREATCH P.The elastic stiff ness moduli of diamond[J].J Appl Phys,1972,43:985-987.

[42]BERNHOLC J,ANTONELLI A,DELSOLE T M,et al.Mechanism of self-diff usion in diamond[J].Phys Rev Lett,1988,61:2689-2692.

[43]LEE D H,JOANNOPOULOS J D.Simple scheme for deriving atomic force constants:Application to SiC[J]. Phys Rev Lett,1982,48:1846-1849.

(責(zé)任編輯：李藝)

Interatomic potential fi tting study of Ag,Si and C based on
fi rst-principle calculations

WANG Yong-wei1,GUO Yong-liang1,2,ZHANG Wei2,KE Xue-zhi1
(1.School of Physics and Materials Science,East China Normal University, Shanghai 200241,China; 2.Shanghai Institute of Applied Physics,Chinese Academy of Sciences, Shanghai 201800,China)

In order to perform molecular dynamics simulations of Ag diff usion in SiC crystals,we use”force-matching”method to fi t the interatomic interaction potentials of Ag, Si and C based on fi rst-principle calculations.The eff ectiveness of our obtained potential functions are verifi ed by the calculations of the lattice constants,cohesive energies,bulk modulus,elastic constants and defect formation energies,etc.The results show that the values of cohesive energies,lattice constants and bulk modulus can be calculated precisely by our potentials,and the maximum error does not exceed 0.6%.The values of vacanciesand interstitial formation energies of Si and C crystals and Si and C vacancy formation of SiC crystals calculated by our potentials are more accurate than that calculated by J.Tersoff’s potentials.In addition,the formation energies of the 16 kinds of AgSiC three-atom defect systems can also be calculated accurately by our potentials.

Tersoff potential;force-matching;SiC;first-principle calculations

O411.2;O483;O414.2

A

10.3969/j.issn.1000-5641.2017.04.010

1000-5641(2017)04-0114-12

2016-07-22

王永偉,男,碩士研究生,研究方向?yàn)榻y(tǒng)計(jì)和凝聚態(tài)理論.E-mail：wangyongwe13@163.com.

柯學(xué)志,男,教授,博士生導(dǎo)師,研究方向?yàn)榻y(tǒng)計(jì)和凝聚態(tài)理論. E-mail：xzke@phy.ecnu.edu.cn.郭永亮,男,博士研究生,研究方向?yàn)榻y(tǒng)計(jì)和凝聚態(tài)理論.E-mail：ylguo@stu.ecnu.edu.cn.