面心立方Ce同構相變的分子動力學模擬

第伍旻杰 胡曉棉

1) (中國工程物理研究院研究生院,北京 100088)

2) (北京應用物理與計算數(shù)學研究所,計算物理國家重點實驗室,北京 100088)

基于嵌入原子法,本文給出了一個金屬Ce原子間的相互作用勢.利用該勢分別計算了γ-Ce和α-Ce的晶格常數(shù)、結合能、彈性常數(shù),計算結果與實驗或第一性原理研究中得出的數(shù)值符合得較好.給出了兩相Ce中如點缺陷形成能、表面能、層錯能以及孿晶能等晶體缺陷形成能.通過分析兩相Ce的聲子譜,得出了不同溫度下兩相的晶格振動熵差,其中在室溫條件下約為0.67kB/atom.還利用分子動力學模擬得出了該相變的等溫線,并且利用徑向分布函數(shù)分析了相變前后兩相的晶體結構,確認了該相變?yōu)槊嫘牧⒎酵瑯嬒嘧?即Ce的α-γ相變.由此表明,本文的嵌入原子法勢,不僅可以分別合理地描述γ-Ce和α-Ce,還可以反映γ-Ce和α-Ce兩相之間的相變.

1 引 言

鈰(cerium,Ce)是一種稀土金屬元素,其f電子處于局域與非局域的轉變點上,具有低熔點、非對稱晶體結構、多種同素異形體共存、發(fā)生相變伴隨體積變化等物理、化學性質[1,2].Ce及其稀土合金已廣泛應用于鋼鐵、有色金屬及其合金、發(fā)火合金、稀土永磁合金、儲氫材料等諸多領域[3].

在較低的溫度和壓強條件下,金屬Ce有兩種面心立方(FCC)相(α相和γ相)以及一種雙密排六方相(β相)[4].在室溫和0.7GPa壓縮條件下Ce會發(fā)生γ → α同構相變[4-6],晶格常數(shù)從5.16 ?突變到4.84 ?,同時體積減小14%-17%.α-γ相變的臨界點(C.P.)壓強PC=1.44-1.96 GPa[4,7,8],溫度TC=460-600 K[4,7,8].

作為Ce的α-γ相變的機制,從電子結構來看,Mott相變和Kondo體積坍縮(KVC)兩種理論模型長期爭論不下[4,9,10].也有研究認為,這兩種模型分別反映了Ce的α-γ相變的不同側面[11].除此之外,Amadon等[12]認為在室溫條件下熵的效應不可忽略,Ce的α-γ相變是由熵驅動的.而熵的貢獻可分為電子的和晶格振動的兩個自由度.利用高壓高分辨率同步X射線粉末衍射實驗,Jeong等[13]測得Ce的α-γ相變過程中,晶格振動的熵變(ΔSvibα-γ≈(0.75 ± 0.15)kB/atom)約占總熵變(約1.54kB/atom)的一半.Wang等[6]利用第一性原理計算得出的總熵變約為2.86kB/atom,其中晶格振動的熵變約為0.94kB/atom.這些都表明晶格振動的熵變在總熵變中占據(jù)相當可觀的一部分.

對于Ce的α-γ相變,以往的模擬計算工作中主要分宏觀數(shù)值模擬和第一性原理電子結構計算兩方面.基于多相物態(tài)方程的宏觀數(shù)值模擬[14-16]無法給出相變過程中微尺度結構的細節(jié); 而相關的第一性原理計算[6,11,17-19]受制于其時間和空間尺度,只能對此相變進行均勻、平衡態(tài)的分析.而非平衡態(tài)分子動力學模擬因其特點,可以在微、納米尺度提供更多的信息,如相變過程中微結構的演化,從而可以將離散原子尺度與宏觀尺度溝通起來.

在分子動力學模擬中,分子間相互作用勢必不可少.已有工作中針對FCC結構Ce的經(jīng)驗勢不多：例如Singh和Singh[20]將兩體作用的雜化近自由電子緊束縛成鍵模型用在一些FCC結構的稀土金屬(包括Ce),得出的結果在聲子譜和彈性常數(shù)等方面與實驗符合較好; Hachiya和Ito[21]在此基礎上利用近自由電子鍵角相關緊束縛成鍵模型(最早用于過渡金屬及其合金)進行了改進; Sheng等[22]利用勢能面擬合法,得出了包括Ce的14中FCC結構金屬的嵌入原子法(EAM)勢; Fu和Zhao[23]擬合了FCC鑭(La)和Ce的Gupta勢.雖然這些經(jīng)驗作用勢都與Ce的實驗和第一性原理的數(shù)據(jù)符合得很好,但只涉及了γ-Ce,而與α-Ce無關,因此無法用于Ce的α-γ相變的分子動力學模擬.Voter和Chen[24]通過對Voter-Chen EAM進行修正,提高兩相之間的勢壘(0.142 eV/atom)來維持各相的穩(wěn)定性,提出過一種能夠描述Ce的α-γ相變的經(jīng)驗勢[25,26],不過該勢的具體描述至今未見發(fā)表.為了能夠對Ce的α-γ相變進行大規(guī)模分子動力學模擬,本文展示了一種新的FCC結構Ce的EAM勢,該勢可對γ-Ce,α-Ce以及Ce的α-γ相變進行合理描述.

2 EAM模型勢的構建

α-Ce和γ-Ce的冷能曲線可以分別用Rose物態(tài)方程來表示[25,27],也就是說FCC結構Ce的冷能曲線應當是一個分段函數(shù)：晶格常數(shù)較大的部分為γ-Ce的Rose曲線,晶格常數(shù)較小的部分則為α-Ce的Rose曲線.兩段各有一個極小值點對應于各相平衡位置的晶格常數(shù),兩段曲線交叉于兩極小值點中間形成一個勢壘.

在嵌入原子法中,單類型原子系統(tǒng)的總勢能Etot可以表示為

其中Φ(rij)為原子i與j之間的對勢,F為原子i處于電子云密度ρi的嵌入能,f(rij)為原子j在原子i處的電子密度分布函數(shù).

在文獻[28]中,嵌入能F和電子密度ρ的關系表示為

式中a,b,c為待定參數(shù).Φ(rij)和f(rij)的解析形式為立方項求和：

可以根據(jù)精度要求選取n的上限,an,rn,bn,sn為待定參數(shù); Θ(x)為Heaviside函數(shù),若x ≤ 0,Θ(x)=0; 若x＞0,Θ(x)=1.已經(jīng)有研究將其用于擬合面心立方金屬的原子間相互作用勢.在本文工作中先對Φ(rij)和f(rij)分別只取一項,即

且令bγ=1 ?-3,sγ=6.6 ? (此即相互作用勢的截斷距離).利用最小二乘法,以γ-Ce的Rose物態(tài)方程曲線和彈性常數(shù)為目標對(2)式-(4)式中的參數(shù)a,b,c,aγ,rγ進行擬合.由此得出的γ-Ce的冷能曲線與α-Ce的Rose物態(tài)方程曲線相交于r=3.498 ?,于是對Φ(r),f(r)和F(ρ)函數(shù)應當采用分段函數(shù)的形式(分段點為r=3.498 ?,ρ=389.094).由于勢函數(shù)的長程部分(r＞3.498 ?)已經(jīng)確定,接下來只能調整短程的部分(r ＜ 3.498 ?).在r ＜ 3.498 ?的范圍內仍采用相似的形式,同時為了保證勢函數(shù)的連續(xù)性,則有

嵌入函數(shù)F(ρ)則通過嵌入能F (即總勢能與對勢能部分之差)與電子云密度ρ的對應關系利用立方樣條插值得出.(5)式和(6)式中的參數(shù)aα,rα,Φ0,bα,sα,f0參照α-Ce的Rose物態(tài)方程曲線和彈性常數(shù),并考慮α-γ相變的溫度和壓強條件來確定.最終的對勢函數(shù)Φ(r)、電子密度函數(shù)f(r)、嵌入能函數(shù)F(ρ)如圖1(a)和圖1(b)所示.

圖1 Ce的EAM勢 (a)對勢函數(shù)和電子密度分布函數(shù); (b)嵌入能函數(shù)Fig.1.EAM potential for Ce：(a) The pair function Φ(rij) and the density function f(rij); (b) the embedded function F(ρ).

3 結果與討論

為了驗證第2節(jié)中所得的EAM勢,本文采用了LAMMPS[29]對面心立方Ce晶體的一些性質以及α-γ相變過程進行了分子動力學模擬計算.

3.1 基本性質

圖2 由本文EAM勢得出的面心立方Ce的冷能曲線Fig.2.The cold energy of fcc Ce calculated with the newly developed EAM potential.

圖2展示了由EAM勢計算得到的Ce的冷能曲線.該曲線分為兩段,每段各有一個極小點,分別對應于α-Ce (a=4.81 ?)和γ-Ce (a=5.14 ?).可以看出,在兩個極小值點中間存在一個高度為15.3-20.8 meV/atom的勢壘,遠小于Chen等[25]的經(jīng)驗勢中的對應值(0.142 eV/atom).表1中展示了由EAM計算得出的α-Ce和γ-Ce的一些基本性質,包括晶格常數(shù)、結合能、彈性常數(shù)等,與實驗或第一性原理計算的結果符合得很好.

3.2 晶體缺陷

3.2.1 點缺陷

為了定量地研究在Ce中的點缺陷,采用本文的EAM勢在10×10×10 的FCC超胞中計算了α-Ce和γ-Ce中的空位和填隙原子的形成能.晶體點缺陷的形成能(Ef)的定義為

其中Edef為包含單個點缺陷超胞的最小化總能量,N為無缺陷超胞的原子數(shù)(在此N=4000).于是正號對應于填隙原子,負號對應于空位.Ecoh為無缺陷晶體中單個原子的結合能.

表1 面心立方Ce的基本性質的EAM計算值與實驗和第一性原理結果的比較Table 1.EAM predicted properties of Ce lattice in comparison with experimental and ab initiodata.

對于γ-Ce,填隙原子和空位的形成能分別為1.93和0.85 eV,分別對應于無缺陷單晶中單個原子結合能的44.7%和19.7%; α-Ce中的填隙原子和空位的形成能分別為2.97和1.15 eV,分別對應于無缺陷單晶中單個原子結合能的68.7%和26.6%.以往的研究中所得的FCC過渡金屬的空位能約為1-3 eV[35],這里所得的填隙原子和空位的形成能處在與之可比的范圍內,可以認為是合理的.

3.2.2 表面能

采用本文的EAM勢同時通過Polak-Ribiere共軛梯度法[36]對表面構型進行充分弛豫,以此計算了α-Ce和γ-Ce的三種低密勒指數(shù)(即(100),(110)和(111))晶面的表面形成能,計算結果列在了表2中.對于γ-Ce本文中計算所得的表面能明顯低于以往經(jīng)驗勢的計算結果[22,23].不過兩相Ce中的表面能的大小順序都符合γ(111)＜ γ(100)＜γ(110),表明其中(111)表面最穩(wěn)定,這方面與Sheng等[22]以及Fu和Zhao[23]對γ-Ce的計算相結果相一致.比較α-Ce和γ-Ce相同晶向的表面能,則有γα＜ γγ,如對于(100)晶面γα(100)＜ γγ(100),另外對于γ(110)和γ(111)亦同.到目前尚未見關于金屬Ce表面能實驗測量結果的報道,本文的表面能數(shù)值可作為未來相關測量實驗的參考.

3.2.3 面缺陷

材料的穩(wěn)定和非穩(wěn)定層錯能被認為是決定其力學行為(如金屬的塑性變形)的重要物理量.利用本文的EAM勢計算得出了α-Ce和γ-Ce中堆垛層錯和(111)孿晶缺陷的廣義面缺陷能曲線.廣義層錯能曲線可通過對無缺陷單晶沿(111)晶面進行剛性剪切得到,而廣義孿晶缺陷能曲線則可通過對預置層錯的無缺陷單晶進行剛性剪切得出.最終α-Ce和γ-Ce中的廣義層錯能曲線和廣義孿晶能曲線分別如圖3(a)和圖3(b)所示.

由此可得γ-Ce的非穩(wěn)定層錯能γusf=0.543 J·m-2,穩(wěn)定層錯能γssf=0.457 J·m-2(與Sheng等[22]利用EAM計算的結果相近); 而在α-Ce中γusf=0.821 J·m-2,γssf=0.734 J·m-2.雖然此處的結果與?stlin等[37]的第一性原理計算結果相差較大,但是對比α-Ce和γ-Ce兩相的穩(wěn)定層錯能和非穩(wěn)定層錯能均滿足γγsf＜ γαsf.Swygenhoven等[38]認為要了解FCC晶體變形的機制僅有非穩(wěn)定層錯能或者穩(wěn)定層錯能都是不全面的,還應當比較穩(wěn)定層錯能與非穩(wěn)定層錯能的比值γssf/γusf.該比值接近于1的時候傾向于越過整個勢壘產(chǎn)生全位錯; 反之若該比值較小則傾向于產(chǎn)生偏位錯(作為參考,在Al中該比值為0.97,在Ni中為0.55,在Cu中為0.13).在γ-Ce中γssf/γusf=0.842,α-Ce中γssf/γusf=0.875,數(shù)值較高,與Al中的比值比較接近.因此可以推測本文的EAM勢在MD模擬的FCC結構Ce塑性變形中有可能觀察到全位錯產(chǎn)生.

表2 γ-Ce and α-Ce中晶體缺陷的形成能Table 2.Calculated formation energy of lattice defectsin γ-Ce and α-Ce.

圖3 面心立方Ce中的廣義層錯能和廣義孿晶能 (a) α-Ce; (b) γ-CeFig.3.Generalized stacking fault energy and generalized twining fault energy curve for (a) α-Ce and (b) γ-Ce.

除位錯激活與滑移之外,形變孿晶是FCC結構晶體的另一種形變機制.由本文EAM勢計算得出的非穩(wěn)定孿晶能γutf,在γ-Ce中γutf=0.768 J·m-2,在α-Ce中γutf=1.167 J·m-2.Tadmor和Hai[39,40]發(fā)現(xiàn)形成形變孿晶的趨勢與比值γutf/γusf相關,該比值越高越不容易產(chǎn)生形變孿晶.對于γ-Ce,γutf/γusf=1.414; 對于α-Ce,γutf/γusf=1.420,均接近且高于其在Al中的比值(約1.3).由此可推測利用本文EAM勢進行MD模擬很難出現(xiàn)形變孿晶.

3.3 晶格動力學

3.3.1 聲子態(tài)密度和晶格振動熵

圖4所示為溫度T=300 K條件下晶格常數(shù)分別等于4.89 ? (α-Ce,壓強約0.7 GPa)和5.16 ? (γ-Ce,零壓)狀態(tài)的聲子態(tài)密度.兩相的態(tài)密度有顯著差異：1)“肩部”由1.21 THz (γ-Ce) 移動到了 1.43 THz (α-Ce); 2)后續(xù)增長的尖峰由1.81 THz (γ-Ce) 移動到了 2.31 THz (α-Ce);3)再之后的峰由2.58 THz (γ-Ce)移動到了3.36 THz (α-Ce).

Ce的α-γ相變被認為是由熵驅動的,也就是說熵在其中扮演重要位置.利用圖4中的態(tài)密度,可以分別計算出兩種狀態(tài)的晶格振動的熵(如圖5).圖5中右下角虛線圖為兩種狀態(tài)的晶格振動熵的差值,溫度T=300 K時的熵差為ΔSvib=0.67kB/atom,此與Jeong等[13]的結果(0.75 ±0.15)kB/atom相近.然而由于經(jīng)典分子動力學模擬的特性,無法給出電子結構自由度對熵的貢獻,在此對Ce的α-γ相變過程的經(jīng)典分子動力學模擬中只能將固體系統(tǒng)的熵全部歸于晶格振動自由度.于是本文工作中分子動力學模擬的Ce的α-γ相變過程的總熵變等于晶格振動熵變,即ΔS=ΔSvib=0.67kB/atom (此值較實驗值1.54kB/atom少約0.87 kB/atom).

圖4 面心立方Ce的聲子態(tài)密度Fig.4.Phonon density of states for FCC Ce.

圖5 α-Ce和γ-Ce兩相的晶格振動熵Svib以及熵差ΔSvib(右下插圖)隨溫度T的變化Fig.5.The vibrational entropy Svib and its change ΔSvib between two phases (the inset) plotted as functions of temperature.

3.3.2 聲子色散譜

利用聲子色散曲線可以更嚴格地檢驗經(jīng)驗作用勢.圖6所示為300 K溫度條件下計算利用本文EAM計算所得α-Ce (a=4.89 ?,黑實線) 和γ-Ce(a=5.16 ?,藍色點虛線) 的聲子色散關系譜.兩相均有三支聲學波,其中一支為縱波,另兩支為橫波.之前的實驗和理論研究中發(fā)現(xiàn)對于γ-Ce在[ζζζ]方向橫波的L點處存在反常下沉(軟模)[6,19,41].而本文計算中無論α-Ce還是γ-Ce均未找到此反常下沉.

圖6 面心立方Ce的聲子色散譜Fig.6.Phonon dispersion relations for FCC Ce.

3.4 α-γ相變

利用本文中的EAM勢并通過分子動力學模擬研究了Ce的α-γ相變.對初始為FCC結構的Ce進行靜水壓(各向同性)加載,計算所得的壓強-體積等溫線如圖7(a)-圖7(c)所示,其中圖7(a)和圖7(b)的等溫線由NPT系綜分別增、減壓強得出,圖7(c)的等溫線由NVT系綜弛豫得出.如在溫度為T=300 K的條件下,增大壓強到0.80 GPa時發(fā)生相變,體積突變減小2.68 ?3/atom; 反之,減小壓強到0.59 GPa時發(fā)生逆相變,體積突變增大2.83 ?3/atom.由此看出300 K溫度條件下存在明顯的(0.21 GPa)遲滯.當溫度升高至約550 K時,圖7(a)和圖7(b)中的等溫線變?yōu)楣饣€,體積隨壓強連續(xù)變化.對應圖7(c)中的極小值點與極大值點重合于拐點處,此時拐點的斜率為零.當達到更高的溫度,如600,700和800 K時,拐點的斜率變?yōu)樨撝?由此得出臨界點溫度Tc≈550 K,壓強Pc≈1.21 GPa.根據(jù)圖7(a)和圖7(b)繪出的靜水壓加載下該相變的溫度-壓強相變路徑如圖8所示,此可視為該相變的P-T相圖.

為了確認該相變否是為Ce的α-γ同構相變,利用徑向分布函數(shù)對相變前后的晶格結構進行了分析.圖9中比較了溫度為300 K時不同壓強加載下的徑向分布函數(shù)：零壓加載下的前三個峰值對應的間距值分別為3.62,5.16 和6.30 ?; 增壓加載當壓強增大至0.7 GPa (圖7中A點所指狀態(tài)) 時對應的這三個間距值分別為3.56,5.06 和6.18 ?;減壓加載當壓強降至0.7 GPa時 (圖7中B點所指狀態(tài)) 對應的這三個間距值分別為3.42,4.89和5.98 ?.這三種不同加載情況下的徑向分布函數(shù)均符合FCC結構的最近三階近鄰間距之間的關系.于是可以認為相變前后的兩相均為FCC結構,以上模擬的相變正是Ce的α-γ同構相變,其中晶格常數(shù)較大的一相對應于γ-Ce,而晶格常數(shù)較小的一相對應于α-Ce.也就是說,本文的EAM勢可以用于FCC結構Ce的α-γ同構相變的分子動力學模擬,對α-γ同構相變進行多尺度的研究.

圖7 面心立方Ce的等溫線 (a) NPT增壓; (b) NPT減壓; (c) NVTFig.7.Isotherms of FCC Ce：(a) NPT,pressure increased; (b) NPT,pressure decreased; (c) NVT.

圖8 Ce的γ→α和α→γ相變路徑Fig.8.The path of Ce γ→α and α→γ phase transition.

圖9 零壓以及圖7中A,B兩點狀態(tài)的徑向分布函數(shù)Fig.9.Radial distribution function for zero pressure(black),the A state (red),and the B state (blue) pointed in Fig.7.

4 結 論

根據(jù)Rose物態(tài)方程還有彈性常數(shù)等實驗以及第一性原理計算數(shù)據(jù),本文給出了FCC結構Ce的EAM勢.利用該EAM勢對兩種FCC結構Ce進行了定量的討論.計算得出的兩相Ce的結合能、彈性常數(shù)等一些基本性質,以及晶體缺陷和聲子譜,與實驗和第一性原理計算的結果符合較好.此外該EAM勢可以用于FCC結構Ce的α-γ同構相變的分子動力學模擬,并且通過MD模擬得出的P-T相圖也與實驗相符合.

需要注意的是,以往的研究中表明Ce的α-γ相變是熵驅動的(特別是當處在室溫條件下),其中包括了電子自由度和晶格振動自由度兩個方面.但由于分子動力學模擬的局限性,無法體現(xiàn)電子自由度在其中的貢獻,只能給出晶格振動方面的影響.本文中計算得出的常溫下兩相的晶格振動熵差約為0.67kB/atom.

本文的EAM勢主要針對的是FCC結構的α-Ce和γ-Ce以及Ce的α-γ同構相變,未考慮到其他的因素.如同樣處于較低的溫度和壓強條件下的β相,以及溫度較高時的δ相(體心立方結構)和液態(tài)相.對此需要在后續(xù)的工作中考慮并改進.