蔣元祺

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

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

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

將含有16000個釩原子的BCC晶體結(jié)構(gòu)置于一個立方盒子中, 在常溫常壓(normal pressure and temperature, NPT)系綜下, 按周期性邊界條件, 采用分子動力學(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,時間步長為1 fs, 實時輸出溫度T、勢能P、體積V以及原子的位置坐標(biāo)信息. 在模擬完成后, 利用Ovito程序[50]與多面體模板匹配法(polyhedral template matching method, PTMM)[51]進一步可視化分析熔化過程中體系特征原子團簇的構(gòu)型. 文中涉及標(biāo)準(zhǔn)BCC, FCC, HCP以及ICO的從頭算分子動力學(xué)(ab initioMD)模擬時, 采用的是正則系綜 (NVT), 模擬的時間步長為 5 fs, 模擬步數(shù)為500 步, 總模擬時間為 2.5 ps.

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

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

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

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

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

結(jié)合上述RDF的定義以及文獻[1]中的示意圖, 圖 4 中給出了金屬釩以 1 × 1012K/s的速率熔化時的RDF曲線圖. 可知, RDF曲線中峰和谷的變化過程可以分為如下四個階段: 首先在低溫300 K時, 明銳分立的衍射峰反映出完美的晶體結(jié)構(gòu)特征; 在隨后的 400—3200 K (溫度間隔 DT=100 K)的溫度區(qū)間內(nèi), 隨著溫度的逐漸升高, 分立尖銳的峰逐漸變得模糊, 尤其是在r= 5 ?處的兩個較窄的、分立的峰逐漸合并為一個較寬的峰, 即300 K時的第二峰逐漸減小并融入第三峰; 當(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ī)律具有很好的對應(yīng)關(guān)系.

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

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

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

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

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

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

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

圖 7 利用多面體模板匹配法分析得到的5種不同熔化速率下各種類型團簇原子分布隨溫度的演變關(guān)系 (a) FCC; (b) HCP;(c) BCC; (d) ICO; (e) SC; (f) 其他類型團簇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 團簇平均每原子勢能隨模擬步長的演變趨勢 (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.

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

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

圖 9 FCC, HCP, BCC 及 ICO 原子團簇的熵與吉布斯自由能隨溫度的演變關(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 個釩原子在同一初始幾何結(jié)構(gòu)下局域特征原子結(jié)構(gòu)的熔化行為. 研究表明: 在不同速率的熔化過程中, 釩的熔點隨著速率的升高而依次升高, 熔化速率對不同類型原子團簇的絕對數(shù)目影響顯著, 但對不同類型原子團簇數(shù)目分布的相對次序無顯著影響. 針對特征原子團簇進行的ab initioMD模擬和熱力學(xué)分析表明, ICO原子團簇在高溫區(qū)之所以能夠持續(xù)存在, 不僅是因為ICO原子團簇的穩(wěn)定性和團簇壽命明顯優(yōu)于FCC, HCP以及BCC, 而且是因為ICO原子團簇在高溫區(qū)具有相對較高的團簇熵與相對較低的吉布斯自由能.