蔣元祺

1) (南昌師范學(xué)院物理系, 南昌 330032)

2) (吉林大學(xué)物理學(xué)院, 長春 130012)

1 引 言

凝固與熔化是自然界中廣泛存在的兩類自然現(xiàn)象[1,2], 在工業(yè)生產(chǎn)和日常生活中發(fā)揮著重要作用, 是研究物質(zhì)相變和結(jié)構(gòu)演化的有效途徑. 通常認(rèn)為結(jié)晶應(yīng)該只發(fā)生在晶化溫度附近, 熔化應(yīng)該只發(fā)生在熔點(diǎn)溫度處, 但現(xiàn)實(shí)中存在“過冷”[1]與“預(yù)熔”[2]現(xiàn)象, 且在凝固的形核機(jī)制與晶體的熔化機(jī)制方面仍存在諸多未解之謎[1,2]. 近年針對(duì)液態(tài)過渡金屬 (transition metal, TM)的快凝過程, 發(fā)現(xiàn)其中存在著大量的五次對(duì)稱性原子團(tuán)與二十面體團(tuán)簇[1], 它們不但可以抑制形核[1?4], 影響金屬玻璃的形成能力[5,6], 而且還與部分畸變的晶體型原子團(tuán)簇之間存在競爭與演變現(xiàn)象[7,8]. 2014年Zhong等[9]利用1014K/s的超高冷卻速率, 首次在實(shí)驗(yàn)上將具有體心立方(BCC)結(jié)構(gòu)的液態(tài)難熔單質(zhì)金屬 Ta, V, Mo 以及 W 制備成了金屬玻璃. 這極大地鼓舞了從事計(jì)算模擬的研究人員[10,11]. 他們以此為基礎(chǔ), 先后從特征原子團(tuán)簇的競爭與演變[7]、化學(xué)硬度的高低[8]、近黃金比序[11]的新型原子結(jié)構(gòu)以及化學(xué)序參數(shù)的變化[12]等角度, 相繼揭示了液態(tài)Cu-Zr合金與單質(zhì)金屬Ag, Ta等在快凝過程中其原子與電子層面的物理機(jī)制. 此外還有大量工作深入分析了團(tuán)簇的遺傳[13,14]、自組織[15]及其枝狀連接[16], 闡述了團(tuán)簇結(jié)構(gòu)與玻璃形成能力之間的關(guān)系[17,18], 促進(jìn)了學(xué)界對(duì)凝固機(jī)理的理解, 進(jìn)一步深化了對(duì)非晶形成的認(rèn)識(shí)[19,20].

不過遺憾的是, 近二十年來對(duì)同樣是以溫度為主導(dǎo)的自然現(xiàn)象—金屬熔化的研究力度遠(yuǎn)不及對(duì)凝固的研究, 絕大部分研究集中在晶體結(jié)構(gòu)的識(shí)別[21]、晶粒生長[22]、納米顆粒的表面預(yù)熔[23]、晶格振動(dòng)[24]及之后液滴的形成[25]等, 對(duì)其中團(tuán)簇的研究也主要側(cè)重于結(jié)構(gòu)穩(wěn)定性[26,27]、化學(xué)序[28,29]、熱力學(xué)[30,31]與動(dòng)力學(xué)行為[32]、團(tuán)簇結(jié)構(gòu)競爭[33]與團(tuán)簇相轉(zhuǎn)變[33,34]、不同初始構(gòu)型[35]及不同幾何尺寸[36,37]對(duì)團(tuán)簇熔點(diǎn)的影響[38,39]、金屬小團(tuán)簇熔點(diǎn)與其塊體材料熔點(diǎn)的異同[40,41](包括團(tuán)簇原子數(shù)太少而難以熔化的機(jī)理解釋[42])、核-殼[43]二元反對(duì)稱Mackay二十面體與扶手型結(jié)構(gòu)之間的轉(zhuǎn)變[44,45], 以及近兩年新興起的用機(jī)器學(xué)習(xí)[46]來研究和識(shí)別納米顆粒熔化過程中的原子分布與競爭關(guān)系等方面. 而對(duì)熔化過程中特征原子團(tuán)簇的結(jié)構(gòu)演變、分布以及不同熔化速率對(duì)微觀原子結(jié)構(gòu)的影響涉及較少, 一些具體科學(xué)問題仍不清楚[2,47], 且針對(duì)原子團(tuán)簇的相關(guān)學(xué)術(shù)報(bào)道呈逐年下降趨勢. 當(dāng)然這與當(dāng)前的研究熱點(diǎn)有關(guān), 但對(duì)熔化現(xiàn)象, 尤其是對(duì)其中典型原子團(tuán)簇的深入研究會(huì)更好地促進(jìn)對(duì)一級(jí)相變過程的理解. 本文以含有16000個(gè)V原子的難熔金屬為研究對(duì)象, 以300 K為熔化的起始溫度, 模擬了其在5種速率下的熔化過程, 初步分析了BCC、面心立方(FCC)、六角密堆(HCP)、二十面體(ICO)以及簡單立方(SC)原子團(tuán)隨溫度及熔化速率的演變關(guān)系, 基本澄清了隱藏在其背后的物理機(jī)制.

2 計(jì)算方法與參數(shù)設(shè)置

將含有16000個(gè)釩原子的BCC晶體結(jié)構(gòu)置于一個(gè)立方盒子中, 在常溫常壓(normal pressure and temperature, NPT)系綜下, 按周期性邊界條件, 采用分子動(dòng)力學(xué) (molecular dynamics, MD)程序LAMMPS[48]來模擬其熔化過程. 原子間相互作用勢采用嵌入原子勢(embedded-atom model,EAM)[49]. 起始熔化溫度為 300 K, 熔化速率分別為g1= 1 × 1011K/s,g2= 1 × 1012K/s,g3= 1 ×1013K/s,g4= 1 × 1014K/s 與g5= 1 × 1015K/s,時(shí)間步長為1 fs, 實(shí)時(shí)輸出溫度T、勢能P、體積V以及原子的位置坐標(biāo)信息. 在模擬完成后, 利用Ovito程序[50]與多面體模板匹配法(polyhedral template matching method, PTMM)[51]進(jìn)一步可視化分析熔化過程中體系特征原子團(tuán)簇的構(gòu)型. 文中涉及標(biāo)準(zhǔn)BCC, FCC, HCP以及ICO的從頭算分子動(dòng)力學(xué)(ab initioMD)模擬時(shí), 采用的是正則系綜 (NVT), 模擬的時(shí)間步長為 5 fs, 模擬步數(shù)為500 步, 總模擬時(shí)間為 2.5 ps.

3 結(jié)果分析與討論

3.1 體系勢能、體積與熔化速率的關(guān)系

針對(duì)含有16000個(gè)釩原子的BCC結(jié)構(gòu), 圖1中首先給出了5種不同熔化速率下勢能曲線隨溫度的演變過程. 從圖 1(a)與(b)中可以看出, 4種熔化速率 (g1= 1 × 1011K/s,g2= 1 × 1012K/s,g3= 1 × 1013K/s,g4= 1 × 1014K/s)的能量曲線存在明顯的突變, 其熔化轉(zhuǎn)變溫度分別為Tm1=3246 K,Tm2= 3309 K,Tm3= 3384 K 以及Tm4=3553 K. 而隨著熔化速率進(jìn)一步升高至 1×1015K/s,體系勢能曲線的突變消失, 在T= 3250 K 附近呈現(xiàn)漸變現(xiàn)象, 直至4236 K與另外4種熔化曲線重合. 即熔化速率越高,Tm也越高. 導(dǎo)致這一現(xiàn)象產(chǎn)生的原因是過高的熔化速率使得原子體系在極短的時(shí)間內(nèi)來不及達(dá)到穩(wěn)定的動(dòng)力學(xué)與熱力學(xué)平衡.與此同時(shí), 圖2中也同時(shí)給出了5種不同熔化速率下體積隨溫度的演變過程. 從圖2知, 在熔化速率g≤g4(1 × 1014K/s)時(shí), 不同熔化速率下體系的體積存在明顯的突變, 說明有一級(jí)相變產(chǎn)生, 其突變所對(duì)應(yīng)的溫度點(diǎn)與圖1中勢能曲線的轉(zhuǎn)變點(diǎn)呈現(xiàn)完美的一一對(duì)應(yīng)關(guān)系, 反映出固態(tài)釩的體積明顯小于液態(tài)釩的體積. 當(dāng)考慮不同熔化速率下體系的壓強(qiáng)隨溫度的波動(dòng)變化時(shí), 圖3中明確顯示熔化速率最低時(shí) (g1= 1 × 1011K/s), 壓強(qiáng)波動(dòng)最大, 而當(dāng)熔化速率越來越高時(shí), 壓強(qiáng)的波動(dòng)幅度則逐漸減小, 當(dāng)速率升高至 1 × 1015K/s時(shí), 壓強(qiáng)的變化幅度在圖3中是最小的. 該現(xiàn)象說明包含16000個(gè)釩原子的周期性模擬體系在熔化速率較低時(shí), 研究對(duì)象有相對(duì)充分的時(shí)間來達(dá)到相對(duì)穩(wěn)定的亞穩(wěn)或者穩(wěn)定平衡狀態(tài), 而速率較高時(shí), 體系則經(jīng)歷劇烈的熱力學(xué)過程.

圖 1 5種不同熔化速率下金屬釩體系的勢能變化與溫度之間的演變關(guān)系, (b)是(a)的局部放大圖Fig. 1. The potential energy of the vanadium metal as a function of temperature at various melting rates g, (b) is a partial enlarged view of (a).

圖 2 5種不同熔化速率下金屬釩體系的體積變化與溫度之間的演變關(guān)系, (b)是(a)的局部放大圖Fig. 2. The volume of the vanadium metal as a function of temperature at various melting rates g, (b) is a partial enlarged view of (a).

圖 3 不同熔化速率下體系壓強(qiáng)隨溫度的演變關(guān)系Fig. 3. The pressure of the vanadium metal as a function of temperature at various melting rates g.

3.2 徑向分布函數(shù)與微結(jié)構(gòu)的演化

針對(duì)熔化過程中系統(tǒng)微結(jié)構(gòu)的演變, 徑向分布函數(shù) (radial distribution function, RDF)是一種簡單有效的方法, 與實(shí)驗(yàn)上經(jīng)常使用的結(jié)構(gòu)因子互為傅里葉變換, 是表征結(jié)構(gòu)的一個(gè)重要參數(shù), 數(shù)學(xué)定義如下[1]:

其gαβ(r) 表示以所有α粒子為中心, 在距中心原子r處的球殼內(nèi)的單位體積中找到β粒子的平均幾率.N是粒子總數(shù),Nα和Nβ分別是α類和β類原子 的 數(shù) 目 ,|rij|是 第i個(gè) 原 子 (α類 )和 第j(β類)個(gè)原子間的距離,ρ是原子數(shù)密度. 對(duì)于不區(qū)分原子種類的單原子系統(tǒng), 其表達(dá)式簡化如下[1]:

結(jié)合上述RDF的定義以及文獻(xiàn)[1]中的示意圖, 圖 4 中給出了金屬釩以 1 × 1012K/s的速率熔化時(shí)的RDF曲線圖. 可知, RDF曲線中峰和谷的變化過程可以分為如下四個(gè)階段: 首先在低溫300 K時(shí), 明銳分立的衍射峰反映出完美的晶體結(jié)構(gòu)特征; 在隨后的 400—3200 K (溫度間隔 DT=100 K)的溫度區(qū)間內(nèi), 隨著溫度的逐漸升高, 分立尖銳的峰逐漸變得模糊, 尤其是在r= 5 ?處的兩個(gè)較窄的、分立的峰逐漸合并為一個(gè)較寬的峰, 即300 K時(shí)的第二峰逐漸減小并融入第三峰; 當(dāng)溫度繼續(xù)升高, 在 3200 K ≤T≤ 3320 K 的溫度區(qū)間內(nèi), 除第一峰外, RDF曲線其余峰變得更加平緩,但仍未完全消失, 呈現(xiàn)出典型的固態(tài)向液態(tài)熔化演變的過渡特征; 當(dāng)溫度繼續(xù)升高至3320 K以上,在r≥ 5 ?的區(qū)域內(nèi) RDF 系列峰逐漸消失, 呈現(xiàn)液態(tài)金屬的典型特征. RDF曲線的峰與谷的演變規(guī)律與勢能曲線與溫度之間的演變規(guī)律具有很好的對(duì)應(yīng)關(guān)系.

圖 4 熔化速率為 g2 = 1 × 1012 K/s時(shí) RDF 在不同溫度區(qū)間的演變關(guān)系Fig. 4. The RDF for vanadium metal at different temperature obtained by rapid heating rate g2 = 1 × 1012 K/s.

在此基礎(chǔ)上, 圖5更加直觀地給出了釩晶體熔化過程中不同類型原子團(tuán)簇(FCC, HCP, BCC,ICO 以及 SC)的演變示意圖. 可知 300 K 時(shí), 體系的確呈現(xiàn)完美的晶體結(jié)構(gòu), 原子排列整齊, 然而當(dāng)溫度繼續(xù)升高至3300 K時(shí), 原子已經(jīng)明顯偏離了其平衡位置, BCC顯著減少, 之后當(dāng)溫度每升高10 K其結(jié)構(gòu)變化都是異常劇烈的. 當(dāng)T升高至4500 K時(shí), 體系的原子排布已經(jīng)是完全“無序”的排列, 直觀地印證了圖4中RDF反映的體系結(jié)構(gòu)特征. 與此同時(shí), 也要承認(rèn)以速率 1 × 1012K/s進(jìn)行的熔化所得的熔點(diǎn)溫度與實(shí)驗(yàn)熔化溫度有一定差距, 這主要是由于熔化速率過高所致.

3.3 不同熔化速率下微結(jié)構(gòu)的分布與演化

在完成5種不同熔化速率的MD模擬之后,要進(jìn)一步獲悉熔化過程中晶體型原子團(tuán)簇之間的演變和彼此之間的相互轉(zhuǎn)化, 以及要發(fā)生這種轉(zhuǎn)化所需的基本條件等, 首先就需要對(duì)其中承載基本信息的原子團(tuán)簇進(jìn)行識(shí)別. 本文將運(yùn)用多面體模板匹配法 (polyhedral template matching method, PTMM)[51]來完成這一任務(wù). 簡單來說, PTMM 方法是基于均方根偏差(root-mean-square deviation,RMSD)的數(shù)學(xué)方法來判斷兩點(diǎn)之間結(jié)構(gòu)的相似性.對(duì)于給定的兩個(gè)原子A與B, 其RMSD可定義為

圖 5 晶體釩熔化過程不同類型團(tuán)簇的演變示意圖 (g2 =1 × 1012 K/s)Fig. 5. Different type atomic schematic diagraph in vanadium system at heating rate g2 = 1 × 1012 K/s.

而對(duì)于兩個(gè)以上的多原子體系的疊加問題, 主要是要找到一個(gè)旋轉(zhuǎn)矢量A與一個(gè)縮放矢量B, 使得RMSD最小. 這種情況下RMSD的定義式轉(zhuǎn)變?nèi)缦耓51]:

其中N是原子數(shù),s是A的最佳縮放因子,Q是右手 正 交 矩 陣 ,分 別是A與B原子的重心,ai與bi分別代表A與B原子的任意變化矢量. 通?？梢岳谜痪仃嚺c縮放因子的共同調(diào)節(jié)使得RMSD可以保持一個(gè)測定的不變量, 從而獲得變換后各點(diǎn)距原點(diǎn)的距離[51]:

圖 6 利用多面體模板匹配法分析得到的5種不同熔化速率下FCC, HCP, BCC, ICO以及SC類型團(tuán)簇原子隨溫度的演變關(guān)系Fig. 6. The population of the FCC, HCP, BCC, ICO and SC types atoms in V system as a function of temperature obtained from the polyhedral template matching at various rates, respectively.

眾所周知, 金屬釩在常溫下是BCC結(jié)構(gòu), 圖5中已經(jīng)初步給出了 FCC, HCP, BCC, ICO 以及 SC團(tuán)簇原子隨溫度的變化示意圖, 在此基礎(chǔ)上, 圖6給出了 5 種不同熔化速率 (g1= 1 × 1011K/s,g2=1 × 1012K/s,g3= 1 × 1013K/s,g4= 1 × 1014K/s以及g5= 1 × 1015K/s) 下 FCC, HCP, BCC, ICO以及SC團(tuán)簇原子數(shù)隨溫度變化的具體情況. 從圖6總的來看, 在低溫區(qū)間內(nèi), 不論是高熔化速率還是低熔化速率, 當(dāng)溫度由室溫逐漸升高時(shí), BCC團(tuán)簇在起始階段呈現(xiàn)緩慢減少的趨勢, 但當(dāng)溫度逐漸逼近Tm時(shí), BCC原子急劇減少, HCP則迅速增加,之后 ICO, FCC以及 SC等原子團(tuán)簇也在T=Tm這一瞬間陡然增加. 而一旦當(dāng)溫度T>Tm時(shí),HCP與ICO則由最大值開始急劇下降, 同時(shí)FCC也開始減少. 與此同時(shí)雖然SC的數(shù)量在高溫階段較少, 但其數(shù)目并未隨溫度的增加而發(fā)生明顯變化. 簡言之, 圖 6中隱藏著如下的物理規(guī)律: 即在BCC急劇消失的時(shí)候, 在熔化溫度附近有部分原子團(tuán)簇快速轉(zhuǎn)化為HCP, ICO, FCC以及SC團(tuán)簇, 這一熔化規(guī)律并不隨熔化速率的升高而發(fā)生改變. 之后為了進(jìn)一步了解晶體釩在熔化過程中, 不同熔化速率對(duì)同一類型團(tuán)簇原子數(shù)目的影響, 圖7分別給出了 FCC, HCP, BCC, ICO, SC 以及其他類型團(tuán)簇的原子分布隨溫度的演變關(guān)系. 從圖7可見, 不論是 FCC, HCP, ICO 還是 SC 類型的原子團(tuán)簇, 其數(shù)目都在熔化溫度Tm附近明顯增加, 與此同時(shí)伴隨著BCC的急劇減少(圖7(c)). 對(duì)同一團(tuán)簇而言, 考慮熔化速率效應(yīng), 發(fā)現(xiàn)在速率g≤ 1 ×1012K/s時(shí), 速率的改變對(duì)團(tuán)簇?cái)?shù)目的影響并不顯著. 當(dāng)g大于 1 × 1012K/s時(shí), 熔化速率越高, BCC類型團(tuán)簇減少得越慢, 而 FCC, HCP, ICO 以及SC則增長得越慢. 之后隨著溫度逐漸升高, FCC,HCP以及ICO的數(shù)目在達(dá)到頂峰后才逐漸減少.

圖 7 利用多面體模板匹配法分析得到的5種不同熔化速率下各種類型團(tuán)簇原子分布隨溫度的演變關(guān)系 (a) FCC; (b) HCP;(c) BCC; (d) ICO; (e) SC; (f) 其他類型團(tuán)簇Fig. 7. The fraction of the various types atoms in V system as a function of temperature obtained from the polyhedral template matching at five various rates: (a) FCC; (b) HCP; (c) BCC; (d) ICO; (e) SC; (f) other types atoms counts.

圖 8 FCC, HCP, BCC 以及 ICO 團(tuán)簇平均每原子勢能隨模擬步長的演變趨勢 (a) 300 K; (b) 500 K; (c) 1100 KFig. 8. The variation of the potential energy of per atom of FCC, HCP, BCC and ICO cluster as a function of time step, respectively:(a) 300 K; (b) 500 K; (c) 1100 K.

針對(duì)圖6與圖7中 ICO原子團(tuán)簇在高溫階段持續(xù)存在的物理現(xiàn)象, 圖8進(jìn)一步利用ab initioMD模擬了 BCC(15個(gè)原子)、FCC(13個(gè)原子)、HCP(13個(gè)原子)以及ICO(13個(gè)原子)分別在300,500 K和 1100 K時(shí)的結(jié)構(gòu)穩(wěn)定性與團(tuán)簇壽命. 結(jié)果顯示在 300 K和 500 K時(shí), FCC與 HCP均在模擬結(jié)束后轉(zhuǎn)變?yōu)榱薎CO, BCC也在模擬結(jié)束后轉(zhuǎn)變?yōu)榱撕蠭CO團(tuán)簇核心的幾何結(jié)構(gòu), 而ICO則幾乎保持其初始結(jié)構(gòu)而未發(fā)生變化. 當(dāng)溫度升高至1100 K時(shí), BCC演變?yōu)橐粋€(gè)含有畸變ICO的結(jié)構(gòu), 而FCC與HCP則經(jīng)過一系列不同的亞穩(wěn)結(jié)構(gòu)最終演變?yōu)榉且?guī)則構(gòu)型. 唯獨(dú)ICO仍然保持其基本構(gòu)型, 且未發(fā)生明顯改變, 這說明在高溫階段ICO的穩(wěn)定性的確高于晶體型原子團(tuán)簇, 這也為ICO在溫度高于熔化溫度后的大量出現(xiàn)提供了團(tuán)簇層面的解釋.

此外, 考慮到熔化過程中的MD模擬[35]與實(shí)驗(yàn)[37,40]均是在有限溫度下進(jìn)行的, 所以本文通過計(jì)算團(tuán)簇的原子振動(dòng)頻率[52]后, 進(jìn)一步從熵與吉布斯自由能的角度來研究各團(tuán)簇的熱力學(xué)穩(wěn)定性[52],期望以此找到具有5次對(duì)稱性而缺乏平移對(duì)稱性的ICO能夠在高溫區(qū)穩(wěn)定存在的原始因素. 因?yàn)殪嘏c吉布斯自由能均是描述熱力學(xué)系統(tǒng)的重要態(tài)函數(shù), 其大小均可清晰反映系統(tǒng)所處狀態(tài)的穩(wěn)定情況, 指明熱力學(xué)過程進(jìn)行的方向, 可以為研究對(duì)象提供定量表述, 即熵越高, 吉布斯自由能越低, 孤立體系越穩(wěn)定. 基于此, 圖9中分別給出了FCC,HCP, ICO及BCC(為了與其他三種標(biāo)準(zhǔn)團(tuán)簇具有可比性，此處是僅含13原子的標(biāo)準(zhǔn)BCC碎片)原子團(tuán)簇的熵與吉布斯自由能隨溫度的演變關(guān)系. 可知, 雖然各團(tuán)簇的熵與自由能隨溫度的演變趨勢恰好相反, 但其所指示的物理意義[52]卻完全一致, 即與 FCC, HCP, BCC 相比, ICO 的熵最高, 吉布斯自由能最低, 而且隨著溫度逐漸升高, 該現(xiàn)象愈加顯著. 這充分表明高溫液態(tài)金屬區(qū)域ICO的熱力學(xué)穩(wěn)定性遠(yuǎn)遠(yuǎn)優(yōu)于各晶體型原子團(tuán), 為液態(tài)金屬中ICO原子團(tuán)的大量存在提供了熱力學(xué)層面的理論解釋.

圖 9 FCC, HCP, BCC 及 ICO 原子團(tuán)簇的熵與吉布斯自由能隨溫度的演變關(guān)系 (a) 熵; (b) 吉布斯自由能Fig. 9. The variation of the entropy and Gibbs free energy of FCC, HCP, BCC and ICO cluster as a function of temperature, respectively: (a) Entropy; (b) Gibbs free energy.

4 結(jié) 論

本文采用MD方法, 通過LAMMPS模擬了5 種不同熔化速率 (g1= 1 × 1011K/s,g2= 1 ×1012K/s,g3= 1 × 1013K/s,g4= 1 × 1014K/s與g5= 1 × 1015K/s)下 16000 個(gè)釩原子在同一初始幾何結(jié)構(gòu)下局域特征原子結(jié)構(gòu)的熔化行為. 研究表明: 在不同速率的熔化過程中, 釩的熔點(diǎn)隨著速率的升高而依次升高, 熔化速率對(duì)不同類型原子團(tuán)簇的絕對(duì)數(shù)目影響顯著, 但對(duì)不同類型原子團(tuán)簇?cái)?shù)目分布的相對(duì)次序無顯著影響. 針對(duì)特征原子團(tuán)簇進(jìn)行的ab initioMD模擬和熱力學(xué)分析表明, ICO原子團(tuán)簇在高溫區(qū)之所以能夠持續(xù)存在, 不僅是因?yàn)镮CO原子團(tuán)簇的穩(wěn)定性和團(tuán)簇壽命明顯優(yōu)于FCC, HCP以及BCC, 而且是因?yàn)镮CO原子團(tuán)簇在高溫區(qū)具有相對(duì)較高的團(tuán)簇熵與相對(duì)較低的吉布斯自由能.