摘要:采用经典分子动力学(molecular dynamics, MD)方法, 模拟了16000个钒原子在5种不同熔化速率(γ₁ = 1 × 10¹¹ K/s, γ₂ = 1 × 10¹² K/s, γ₃ = 1 × 10¹³ K/s, γ₄ = 1 × 10¹⁴ K/s 与γ₅ = 1 × 10¹⁵ K/s)下原子结构的熔化行为. 结果表明: 不同熔化速率对难熔金属钒的熔点影响明显, 不过随着温度升高, 体系特征原子结构诸如体心立方(BCC)、六角密堆(HCP)、面心立方(FCC)、简单立方(SC)以及二十面体(ICO)的相对分布次序并不随熔化速率的改变而改变, 温度仍然是影响原子结构分布的主要因素. 通过从头算分子动力学(ab initio MD)与热力学分析发现, ICO能够在液态金属区域稳定存在, 一方面是因为其孤立团簇的相对稳定性和团簇寿命要优于晶体型原子团簇, 另一方面是因为其拥有相对较高的团簇熵与相对较低的自由能.
关键词: 分子动力学/
钒/
熔化/
团簇

English Abstract

全文HTML

--> --> -->

凝固与熔化是自然界中广泛存在的两类自然现象^[1,2], 在工业生产和日常生活中发挥着重要作用, 是研究物质相变和结构演化的有效途径. 通常认为结晶应该只发生在晶化温度附近, 熔化应该只发生在熔点温度处, 但现实中存在“过冷”^[1]与“预熔”^[2]现象, 且在凝固的形核机制与晶体的熔化机制方面仍存在诸多未解之谜^[1,2]. 近年针对液态过渡金属(transition metal, TM)的快凝过程, 发现其中存在着大量的五次对称性原子团与二十面体团簇^[1], 它们不但可以抑制形核^[1-4], 影响金属玻璃的形成能力^[5,6], 而且还与部分畸变的晶体型原子团簇之间存在竞争与演变现象^[7,8]. 2014年Zhong等^[9]利用10¹⁴ K/s的超高冷却速率, 首次在实验上将具有体心立方(BCC)结构的液态难熔单质金属Ta, V, Mo以及W制备成了金属玻璃. 这极大地鼓舞了从事计算模拟的研究人员^[10,11]. 他们以此为基础, 先后从特征原子团簇的竞争与演变^[7]、化学硬度的高低^[8]、近黄金比序^[11]的新型原子结构以及化学序参数的变化^[12]等角度, 相继揭示了液态Cu-Zr合金与单质金属Ag, Ta等在快凝过程中其原子与电子层面的物理机制. 此外还有大量工作深入分析了团簇的遗传^[13,14]、自组织^[15]及其枝状连接^[16], 阐述了团簇结构与玻璃形成能力之间的关系^[17,18], 促进了学界对凝固机理的理解, 进一步深化了对非晶形成的认识^[19,20].
不过遗憾的是, 近二十年来对同样是以温度为主导的自然现象—金属熔化的研究力度远不及对凝固的研究, 绝大部分研究集中在晶体结构的识别^[21]、晶粒生长^[22]、纳米颗粒的表面预熔^[23]、晶格振动^[24]及之后液滴的形成^[25]等, 对其中团簇的研究也主要侧重于结构稳定性^[26,27]、化学序^[28,29]、热力学^[30,31]与动力学行为^[32]、团簇结构竞争^[33]与团簇相转变^[33,34]、不同初始构型^[35]及不同几何尺寸^[36,37]对团簇熔点的影响^[38,39]、金属小团簇熔点与其块体材料熔点的异同^[40,41](包括团簇原子数太少而难以熔化的机理解释^[42])、核-壳^[43]二元反对称Mackay二十面体与扶手型结构之间的转变^[44,45], 以及近两年新兴起的用机器学习^[46]来研究和识别纳米颗粒熔化过程中的原子分布与竞争关系等方面. 而对熔化过程中特征原子团簇的结构演变、分布以及不同熔化速率对微观原子结构的影响涉及较少, 一些具体科学问题仍不清楚^[2,47], 且针对原子团簇的相关学术报道呈逐年下降趋势. 当然这与当前的研究热点有关, 但对熔化现象, 尤其是对其中典型原子团簇的深入研究会更好地促进对一级相变过程的理解. 本文以含有16000个V原子的难熔金属为研究对象, 以300 K为熔化的起始温度, 模拟了其在5种速率下的熔化过程, 初步分析了BCC、面心立方(FCC)、六角密堆(HCP)、二十面体(ICO)以及简单立方(SC)原子团随温度及熔化速率的演变关系, 基本澄清了隐藏在其背后的物理机制.

将含有16000个钒原子的BCC晶体结构置于一个立方盒子中, 在常温常压(normal pressure and temperature, NPT)系综下, 按周期性边界条件, 采用分子动力学(molecular dynamics, MD)程序LAMMPS^[48]来模拟其熔化过程. 原子间相互作用势采用嵌入原子势(embedded-atom model, EAM)^[49]. 起始熔化温度为300 K, 熔化速率分别为γ₁ = 1 × 10¹¹ K/s, γ₂ = 1 × 10¹² K/s, γ₃ = 1 × 10¹³ K/s, γ₄ = 1 × 10¹⁴ K/s 与γ₅ = 1 × 10¹⁵ K/s, 时间步长为1 fs, 实时输出温度T、势能P、体积V以及原子的位置坐标信息. 在模拟完成后, 利用Ovito程序^[50]与多面体模板匹配法(polyhedral template matching method, PTMM)^[51]进一步可视化分析熔化过程中体系特征原子团簇的构型. 文中涉及标准BCC, FCC, HCP以及ICO的从头算分子动力学(ab initio MD)模拟时, 采用的是正则系综(NVT), 模拟的时间步长为5 fs, 模拟步数为500步, 总模拟时间为2.5 ps.

2

-->

3.1.体系势能、体积与熔化速率的关系

针对含有16000个钒原子的BCC结构, 图1中首先给出了5种不同熔化速率下势能曲线随温度的演变过程. 从图1(a)与(b)中可以看出, 4种熔化速率(γ₁ = 1 × 10¹¹ K/s, γ₂ = 1 × 10¹² K/s, γ₃ = 1 × 10¹³ K/s, γ₄ = 1 × 10¹⁴ K/s)的能量曲线存在明显的突变, 其熔化转变温度分别为T_m1 = 3246 K, T_m2 = 3309 K, T_m3 = 3384 K以及T_m4 = 3553 K. 而随着熔化速率进一步升高至$1\!\times \!10^{15}$ K/s, 体系势能曲线的突变消失, 在T = 3250 K附近呈现渐变现象, 直至4236 K与另外4种熔化曲线重合. 即熔化速率越高, T_m也越高. 导致这一现象产生的原因是过高的熔化速率使得原子体系在极短的时间内来不及达到稳定的动力学与热力学平衡. 与此同时, 图2中也同时给出了5种不同熔化速率下体积随温度的演变过程. 从图2知, 在熔化速率γ ≤ γ₄ (1 × 10¹⁴ K/s)时, 不同熔化速率下体系的体积存在明显的突变, 说明有一级相变产生, 其突变所对应的温度点与图1中势能曲线的转变点呈现完美的一一对应关系, 反映出固态钒的体积明显小于液态钒的体积. 当考虑不同熔化速率下体系的压强随温度的波动变化时, 图3中明确显示熔化速率最低时(γ₁ = 1 × 10¹¹ K/s), 压强波动最大, 而当熔化速率越来越高时, 压强的波动幅度则逐渐减小, 当速率升高至1 × 10¹⁵ K/s时, 压强的变化幅度在图3中是最小的. 该现象说明包含16000个钒原子的周期性模拟体系在熔化速率较低时, 研究对象有相对充分的时间来达到相对稳定的亚稳或者稳定平衡状态, 而速率较高时, 体系则经历剧烈的热力学过程.
图 1 5种不同熔化速率下金属钒体系的势能变化与温度之间的演变关系, (b)是(a)的局部放大图
Figure1. The potential energy of the vanadium metal as a function of temperature at various melting rates γ, (b) is a partial enlarged view of (a).

图 2 5种不同熔化速率下金属钒体系的体积变化与温度之间的演变关系, (b)是(a)的局部放大图
Figure2. The volume of the vanadium metal as a function of temperature at various melting rates γ, (b) is a partial enlarged view of (a).

图 3 不同熔化速率下体系压强随温度的演变关系
Figure3. The pressure of the vanadium metal as a function of temperature at various melting rates γ.

2 -->

3.2.径向分布函数与微结构的演化

针对熔化过程中系统微结构的演变, 径向分布函数(radial distribution function, RDF)是一种简单有效的方法, 与实验上经常使用的结构因子互为傅里叶变换, 是表征结构的一个重要参数, 数学定义如下^[1]:

${g_{{\rm{\alpha \beta }}}}(r) = \dfrac{N}{{4{\rm{\pi }}\rho {N_{\rm{\alpha }}}{N_{\rm{\beta }}}}}\sum\limits_{i = 1}^{{N_{\rm{a}}}} {\sum\limits_{j = 1}^{{N_{\rm{\beta }}}} {\delta \left( {r - \left| {{{{r}}_{ij}}} \right|} \right)} }, $

其${g_{{\rm{\alpha \beta }}}}(r)$表示以所有${\rm{\alpha }}$粒子为中心, 在距中心原子r处的球壳内的单位体积中找到${\rm{\beta }}$粒子的平均几率. N是粒子总数, ${N_{\rm{\alpha }}}$和${N_{\rm{\beta }}}$分别是${\rm{\alpha }}$类和${\rm{\beta }}$类原子的数目, $\left| {{{{r}}_{ij}}} \right|$是第i个原子(${\rm{\alpha }}$类)和第j (${\rm{\beta }}$类)个原子间的距离, $\rho $是原子数密度. 对于不区分原子种类的单原子系统, 其表达式简化如下^[1]:

$g(r) = \dfrac{1}{{4{\rm{\pi }}\rho N}}\sum\limits_{i = 1}^N {\sum\limits_{j = 1,j \ne i}^N {\delta (r - \left| {{{{r}}_{ij}}} \right|)} } .$

结合上述RDF的定义以及文献[1]中的示意图, 图4中给出了金属钒以1 × 10¹² K/s的速率熔化时的RDF曲线图. 可知, RDF曲线中峰和谷的变化过程可以分为如下四个阶段: 首先在低温300 K时, 明锐分立的衍射峰反映出完美的晶体结构特征; 在随后的400—3200 K (温度间隔ΔT = 100 K)的温度区间内, 随着温度的逐渐升高, 分立尖锐的峰逐渐变得模糊, 尤其是在r = 5 ?处的两个较窄的、分立的峰逐渐合并为一个较宽的峰, 即300 K时的第二峰逐渐减小并融入第三峰; 当温度继续升高, 在3200 K ≤ T ≤ 3320 K的温度区间内, 除第一峰外, RDF曲线其余峰变得更加平缓, 但仍未完全消失, 呈现出典型的固态向液态熔化演变的过渡特征; 当温度继续升高至3320 K以上, 在r ≥ 5 ?的区域内RDF系列峰逐渐消失, 呈现液态金属的典型特征. RDF曲线的峰与谷的演变规律与势能曲线与温度之间的演变规律具有很好的对应关系.
图 4 熔化速率为γ₂ = 1 × 10¹² K/s时RDF在不同温度区间的演变关系
Figure4. The RDF for vanadium metal at different temperature obtained by rapid heating rate γ₂ = 1 × 10¹² K/s.

在此基础上, 图5更加直观地给出了钒晶体熔化过程中不同类型原子团簇(FCC, HCP, BCC, ICO以及SC)的演变示意图. 可知300 K时, 体系的确呈现完美的晶体结构, 原子排列整齐, 然而当温度继续升高至3300 K时, 原子已经明显偏离了其平衡位置, BCC显著减少, 之后当温度每升高10 K其结构变化都是异常剧烈的. 当T升高至4500 K时, 体系的原子排布已经是完全“无序”的排列, 直观地印证了图4中RDF反映的体系结构特征. 与此同时, 也要承认以速率1 × 10¹² K/s进行的熔化所得的熔点温度与实验熔化温度有一定差距, 这主要是由于熔化速率过高所致.
图 5 晶体钒熔化过程不同类型团簇的演变示意图(γ₂ = 1 × 10¹² K/s)
Figure5. Different type atomic schematic diagraph in vanadium system at heating rate γ₂ = 1 × 10¹² K/s.

2 -->

3.3.不同熔化速率下微结构的分布与演化

在完成5种不同熔化速率的MD模拟之后, 要进一步获悉熔化过程中晶体型原子团簇之间的演变和彼此之间的相互转化, 以及要发生这种转化所需的基本条件等, 首先就需要对其中承载基本信息的原子团簇进行识别. 本文将运用多面体模板匹配法(polyhedral template matching method, PTMM)^[51]来完成这一任务. 简单来说, PTMM方法是基于均方根偏差(root-mean-square deviation, RMSD)的数学方法来判断两点之间结构的相似性. 对于给定的两个原子A与B, 其RMSD可定义为

${\rm{RMSD}}({{A}},{{B}}) = \sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {\left\| {{{{A}}_i} - {{{B}}_j}} \right\|_{}^2} }, $

而对于两个以上的多原子体系的叠加问题, 主要是要找到一个旋转矢量A与一个缩放矢量B, 使得RMSD最小. 这种情况下RMSD的定义式转变如下^[51]:

$\begin{split}&{\rm{RMSD}}({{A}},{{B}}) = \\& \mathop {{\rm{min}}}\limits_{s {{Q}}}\sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {{{\Big\| {s\left[ {{{{a}}_i} - {{A}}} \right] - {{\big( {{{Q}}{{[{{{b}}_i}-{{B}}]}^{\rm{T}}}} \big)}^{\rm{T}}}} \Big\|}^2}} },\end{split}$

其中N是原子数, s是A的最佳缩放因子, Q是右手正交矩阵, ${{A}} = \dfrac{1}{N}\displaystyle\sum\limits_{i = 1}^N {{{ {{a}}}_i}}$与${{B}} = \dfrac{1}{N}\displaystyle\sum\limits_{i = 1}^N {{{ {{b}}}_i}}$分别是A与B原子的重心, ${ {{a}}_i}$与${ {{b}}_i}$分别代表A与B原子的任意变化矢量. 通常可以利用正交矩阵与缩放因子的共同调节使得RMSD可以保持一个测定的不变量, 从而获得变换后各点距原点的距离^[51]:

$\begin{split}&{\rm{RMSD}}({{A}},{{B}}) = \\&\mathop {{\rm{min}}}\limits_{s{{Q}}} \sqrt {\dfrac{1}{N}\sum\limits_{i = 1}^N {{{\Big\| {s\left[ {{{{a}}_i} - {{A}}} \right] - \dfrac{1}{{l({{B}})}}{{\big( {{{Q}}{{[{b_i} - {{B}}]}^{\rm{T}}}} \big)}^{\rm{T}}}} \Big\|}^2}} } ,\end{split}$

其中$l( {{B}}) = \dfrac{1}{N}\displaystyle\sum\nolimits_{i = 1}^N {\left\| {{b_i} - {{B}}} \right\|}$. 之后确定一个中心原子和其近邻原子之间的结构, 就只需先给出标准模板团簇构型(如图5中的ICO, BCC, HCP, FCC及SC多面体), 然后通过RMSD来判断该标准团簇模板与目标团簇构型之间的结构相似性. RMSD数值越小, 条件越发苛刻, 目标结构与模板的相似度也就越高, 所表征的团簇结构也就越接近于模板结构, 该方法称之为PTMM, 目前已被广泛应用于对非晶与玻璃态结构的表征, 接下来将用其分析不同类型团簇的原子结构.
众所周知, 金属钒在常温下是BCC结构, 图5中已经初步给出了FCC, HCP, BCC, ICO以及SC团簇原子随温度的变化示意图, 在此基础上, 图6给出了5种不同熔化速率(γ₁ = 1 × 10¹¹ K/s, γ₂ = 1 × 10¹² K/s, γ₃ = 1 × 10¹³ K/s, γ₄ = 1 × 10¹⁴ K/s以及γ₅ = 1 × 10¹⁵ K/s)下FCC, HCP, BCC, ICO以及SC团簇原子数随温度变化的具体情况. 从图6总的来看, 在低温区间内, 不论是高熔化速率还是低熔化速率, 当温度由室温逐渐升高时, BCC团簇在起始阶段呈现缓慢减少的趋势, 但当温度逐渐逼近T_m时, BCC原子急剧减少, HCP则迅速增加, 之后ICO, FCC以及SC等原子团簇也在T = T_m这一瞬间陡然增加. 而一旦当温度T > T_m时, HCP与ICO则由最大值开始急剧下降, 同时FCC也开始减少. 与此同时虽然SC的数量在高温阶段较少, 但其数目并未随温度的增加而发生明显变化. 简言之, 图6中隐藏着如下的物理规律: 即在BCC急剧消失的时候, 在熔化温度附近有部分原子团簇快速转化为HCP, ICO, FCC以及SC团簇, 这一熔化规律并不随熔化速率的升高而发生改变. 之后为了进一步了解晶体钒在熔化过程中, 不同熔化速率对同一类型团簇原子数目的影响, 图7分别给出了FCC, HCP, BCC, ICO, SC以及其他类型团簇的原子分布随温度的演变关系. 从图7可见, 不论是FCC, HCP, ICO还是SC类型的原子团簇, 其数目都在熔化温度T_m附近明显增加, 与此同时伴随着BCC的急剧减少(图7(c)). 对同一团簇而言, 考虑熔化速率效应, 发现在速率γ ≤ 1 × 10¹² K/s时, 速率的改变对团簇数目的影响并不显著. 当γ大于1 × 10¹² K/s时, 熔化速率越高, BCC类型团簇减少得越慢, 而FCC, HCP, ICO以及SC则增长得越慢. 之后随着温度逐渐升高, FCC, HCP以及ICO的数目在达到顶峰后才逐渐减少.
图 6 利用多面体模板匹配法分析得到的5种不同熔化速率下FCC, HCP, BCC, ICO以及SC类型团簇原子随温度的演变关系
Figure6. 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.

图 7 利用多面体模板匹配法分析得到的5种不同熔化速率下各种类型团簇原子分布随温度的演变关系　(a) FCC; (b) HCP; (c) BCC; (d) ICO; (e) SC; (f)其他类型团簇
Figure7. 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.

针对图6与图7中 ICO原子团簇在高温阶段持续存在的物理现象, 图8进一步利用ab initio MD模拟了BCC(15个原子)、FCC(13个原子)、HCP(13个原子)以及ICO(13个原子)分别在300, 500 K和1100 K时的结构稳定性与团簇寿命. 结果显示在300 K和500 K时, FCC与HCP均在模拟结束后转变为了ICO, BCC也在模拟结束后转变为了含有ICO团簇核心的几何结构, 而ICO则几乎保持其初始结构而未发生变化. 当温度升高至1100 K时, BCC演变为一个含有畸变ICO的结构, 而FCC与HCP则经过一系列不同的亚稳结构最终演变为非规则构型. 唯独ICO仍然保持其基本构型, 且未发生明显改变, 这说明在高温阶段ICO的稳定性的确高于晶体型原子团簇, 这也为ICO在温度高于熔化温度后的大量出现提供了团簇层面的解释.
图 8 FCC, HCP, BCC以及ICO 团簇平均每原子势能随模拟步长的演变趋势　(a) 300 K; (b) 500 K; (c) 1100 K
Figure8. 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

此外, 考虑到熔化过程中的MD模拟^[35]与实验^[37,40]均是在有限温度下进行的, 所以本文通过计算团簇的原子振动频率^[52]后, 进一步从熵与吉布斯自由能的角度来研究各团簇的热力学稳定性^[52], 期望以此找到具有5次对称性而缺乏平移对称性的ICO能够在高温区稳定存在的原始因素. 因为熵与吉布斯自由能均是描述热力学系统的重要态函数, 其大小均可清晰反映系统所处状态的稳定情况, 指明热力学过程进行的方向, 可以为研究对象提供定量表述, 即熵越高, 吉布斯自由能越低, 孤立体系越稳定. 基于此, 图9中分别给出了FCC, HCP, ICO及BCC(为了与其他三种标准团簇具有可比性，此处是仅含13原子的标准BCC碎片)原子团簇的熵与吉布斯自由能随温度的演变关系. 可知, 虽然各团簇的熵与自由能随温度的演变趋势恰好相反, 但其所指示的物理意义^[52]却完全一致, 即与FCC, HCP, BCC相比, ICO的熵最高, 吉布斯自由能最低, 而且随着温度逐渐升高, 该现象愈加显著. 这充分表明高温液态金属区域ICO的热力学稳定性远远优于各晶体型原子团, 为液态金属中ICO原子团的大量存在提供了热力学层面的理论解释.
图 9 FCC, HCP, BCC及ICO 原子团簇的熵与吉布斯自由能随温度的演变关系　(a) 熵; (b) 吉布斯自由能
Figure9. 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.

本文采用MD方法, 通过LAMMPS模拟了5种不同熔化速率(γ₁ = 1 × 10¹¹ K/s, γ₂ = 1 × 10¹² K/s, γ₃ = 1 × 10¹³ K/s, γ₄ = 1 × 10¹⁴ K/s与γ₅ = 1 × 10¹⁵ K/s)下16000个钒原子在同一初始几何结构下局域特征原子结构的熔化行为. 研究表明: 在不同速率的熔化过程中, 钒的熔点随着速率的升高而依次升高, 熔化速率对不同类型原子团簇的绝对数目影响显著, 但对不同类型原子团簇数目分布的相对次序无显著影响. 针对特征原子团簇进行的ab initio MD模拟和热力学分析表明, ICO原子团簇在高温区之所以能够持续存在, 不仅是因为ICO原子团簇的稳定性和团簇寿命明显优于FCC, HCP以及BCC, 而且是因为ICO原子团簇在高温区具有相对较高的团簇熵与相对较低的吉布斯自由能.