分子动力学模拟揭秘SiC高压相变：机器学习势函数与缺陷效应研究

1. 项目概述与核心思路

碳化硅（SiC）作为第三代半导体的核心材料，以其优异的宽带隙、高热导率和高击穿场强，在高温、高频、高功率电子器件领域展现出巨大潜力。然而，当这类器件在极端工况（如航天器再入大气层、核反应堆内壁、高超音速飞行器前缘）下运行时，材料会承受极高的压力和温度，其微观结构可能发生根本性的相变，进而导致宏观性能的急剧退化甚至失效。传统实验手段，如金刚石对顶砧（DAC）结合激光加热，虽然能创造极端条件，但难以实时、原位地观测原子尺度的动态相变过程和缺陷的微观作用机制。

这正是分子动力学（MD）模拟大显身手的地方。简单来说，MD模拟就像给成千上万个原子拍一部“电影”。我们给每个原子赋予初始位置和速度（对应温度），并规定它们之间如何相互作用（通过“力场”或“势函数”描述），然后基于牛顿第二定律，一步步计算每个原子在极短时间（飞秒级）内的受力、加速度和新的位置。通过追踪这部“原子电影”，我们就能直观地看到材料在高温高压下是如何熔化、如何从一种晶体结构转变为另一种、以及缺陷如何像“催化剂”一样加速或改变这一过程。

本次分享的核心，正是基于这样一套MD模拟流程，深入探究SiC在超高压（最高90 GPa，相当于地核深处压力）下的相变行为。我们特别关注两个核心问题：第一，具有岩盐结构（B1相）的SiC是如何熔化成液体的？这个过程是突然的还是渐进的？第二，材料中不可避免存在的原子空位或间隙原子等缺陷，究竟会将相变温度“拉低”多少？这直接关系到材料在实际应用中的可靠性边界。为了高效且准确地回答这些问题，我们采用了“两相模拟”这一经典而强大的技术方案。

2. 模拟方案设计与技术选型解析

2.1 为什么选择“两相模拟”方法？

在确定相变温度时，传统方法是在一系列温度下进行单相模拟，观察序参量（如径向分布函数、平均平方位移）的突变。但这种方法效率较低，且对相变点的判断存在主观性。两相模拟则巧妙地将两种相（这里是固相B1-SiC和液相SiC）放在同一个模拟盒子中，中间以界面隔开。然后，我们在恒温恒压（NPT系综）下运行模拟，观察界面随时间推移的移动方向。

其背后的物理逻辑非常直观：如果当前温度高于真实相变温度，液相更稳定，固-液界面会向固相一侧移动（固体熔化）；反之，界面会向液相一侧移动（液体凝固）。当界面基本保持不动时，系统温度就非常接近真实的固-液两相平衡温度。这种方法相当于直接“询问”系统在给定压力下的平衡态，结果更为直接和可靠。从计算效率上看，一次两相模拟就能提供一个压力点下相变温度的强约束，远优于进行大量单相模拟进行外推。

2.2 势函数的选择：机器学习势函数的优势

MD模拟的准确性，极度依赖于描述原子间相互作用的势函数。对于SiC这样的多元体系，传统经验势（如Tersoff、Vashishta势）虽然计算快，但在极端条件（高压、高温）下的可靠性往往存疑。而第一性原理分子动力学（AIMD）虽然精度高，但计算成本使其无法用于我们所需的大体系（数万原子）和长时程（纳秒级）模拟。

因此，我们选择了基于机器学习（ML）训练的势函数，具体是采用高斯过程回归或神经网络拟合大量第一性原理计算数据得到的势。这类ML势函数能在保持接近AIMD精度的同时，将计算速度提升数个数量级。在本次工作中，我们引用了近期发表的高质量Si-C体系ML势函数，它能够准确描述从低压锌矿（B3）结构到高压岩盐（B1）结构，再到液态乃至缺陷结构的能量和力。这是整个模拟工作的基石，确保了后续所有现象观察和数据分析的物理可信度。

2.3 模拟体系构建与参数设置要点

我们构建了包含16000个原子的初始体系。其中一半原子按B1结构（类似NaCl的立方结构，Si和C各自形成面心立方格子，相互穿插）排列，代表固相；另一半原子则通过将B1结构加热至远高于其熔点后淬火得到无序结构，代表液相。两相之间形成一个清晰的平面界面。

注意：初始液相的制备至关重要。必须确保淬火后的结构是真正均匀、各向同性的液体，而非残留了晶格序的非晶固体。一个检验方法是计算其径向分布函数（RDF），液体的RDF第二峰应是分裂的，而非像非晶固体那样只有一个宽峰。我们通过将体系加热至5000 K以上并保持一段时间，再快速冷却至目标温度区间来制备。

模拟在NPT系综下进行，使用Nosé-Hoover热浴和压浴来控制温度和压力。压力从10 GPa到90 GPa，选取了多个点。温度则根据前期试探性模拟，围绕预估的相变点设置一个区间。时间步长设为1飞秒（fs），每个模拟通常运行3-4纳秒（ns），以确保界面有足够的时间发生可观测的移动或达到稳态。所有模拟均使用LAMMPS这一开源MD软件包完成，其强大的并行计算能力和丰富的系综、fix命令为我们处理复杂边界条件和输出各类物理量提供了便利。

3. 高压下B1-SiC熔化过程的原子尺度解读

3.1 相变过程的动态监测指标

为了定量刻画熔化过程，我们实时监测了多个物理量，它们共同构成了一幅完整的相变动力学图像：

1. B1相分数：通过共近邻分析（CNA）或多面体模板匹配（PTM）等结构识别算法，实时计算体系中仍保持B1结构的原子比例。这是最直接的序参量，其从1下降到0的过程描述了固相的消失。
2. 最大碳团簇尺寸：在SiC熔化或分解过程中，碳原子可能发生偏聚，形成碳团簇。我们通过团簇分析算法追踪最大碳团簇中包含的原子数。这个指标能敏感地反映是否发生了碳的析出或分相，这是某些条件下SiC分解而非熔化的标志。
3. 系统温度与压力张量：监测系统整体温度（T）和压力张量的六个独立分量（Pxx, Pyy, Pzz, Pxy, Pxz, Pyz）。在NPT系综下，温度是设定值，但实际模拟温度会有微小涨落。压力张量则能反映体系是否处于均匀、平衡的力学状态。对角分量（Pxx等）的平均值应接近目标压力，而非对角分量（剪切应力）应接近于零，表明体系没有大的内应力。

3.2 从85到100 GPa：熔化机制的演变

通过分析85、90、100 GPa三个高压点下的模拟数据（对应原文Figure S13），我们观察到了有趣的趋势：

在85 GPa下，B1相分数随时间的下降相对平缓，最大碳团簇尺寸始终很小（几个原子），且压力张量各分量保持稳定。这表明在此压力下，B1-SiC通过均匀成核熔化直接转变为均匀的液态SiC，没有发生碳的显著偏聚。

在90 GPa下，B1相分数的下降过程中出现了一个短暂的“平台期”，同时最大碳团簇尺寸出现了一个明显的峰值，随后又减小。压力张量的非对角分量在此期间出现了短暂的波动。这暗示熔化过程可能更加复杂，界面处可能瞬时形成了亚稳态的富碳区域或团簇，但最终体系仍趋于均匀的液体。

到了100 GPa，B1相分数的下降曲线与85 GPa类似，但最大碳团簇的基线尺寸略有增大。更重要的是，压力张量的非对角分量在整个模拟期间都显示出比低压下更大的涨落。这说明在极高��压力下，熔化的液相其原子堆垛更加紧密，局部剪切变形更容易发生，液体的“粘性”或结构弛豫特性可能发生了变化。

实操心得：监测压力张量的非对角分量是一个常被忽视但极其有用的诊断工具。如果Pxy, Pxz, Pyz持续偏离零值且数值较大，往往意味着模拟盒子存在较大的内应力，或者体系（特别是固体部分）发生了塑性变形或滑移。在两相模拟中，这可能是界面结构不合理或模拟未充分平衡的信号。一个稳定的两相模拟，其非对角压力分量应在零附近对称涨落。

3.3 如何准确判定相变温度？

从两相模拟中提取相变温度，需要综合判断。我们的判据是：当B1相分数在至少1-2 ns的模拟时间内保持基本恒定（变化幅度<5%），且体系温度、压力、团簇尺寸等指标均达到稳态时，即认为当前温度即为该压力下的固-液平衡温度（Tm）。如果界面缓慢向固相移动，则当前温度高于Tm；反之则低于Tm。通过在不同温度点进行几次模拟，我们可以将Tm定位在一个很窄的区间内（例如±50 K）。

这种方法得到的相变温度，其物理图像非常清晰，避免了从单相模拟外推时对“突变点”定义的主观性。我们将不同压力下的Tm连接起来，就得到了SiC在高压下的熔化相界线。

4. 缺陷浓度对相变温度影响的定量分析

4.1 缺陷模型的构建与引入

实际材料中总是存在缺陷。为了研究其影响，我们在完美的B1-SiC晶体中引入了不同浓度的点缺陷。我们主要研究了两种典型缺陷：肖特基缺陷对（同时移除一个Si原子和一个C原子，形成一对空位，保持电中性）和反位缺陷（一个Si原子和相邻的C原子交换位置）。缺陷浓度设置为0（完美晶体）、1/1000（每1000个原子对中有1个缺陷对）和1/100（高缺陷浓度）。

注意：引入缺陷时，需要特别注意体系的弛豫。创建缺陷后，必须先在0 K下进行能量最小化，然后在目标温度和压力下进行足够长时间的NPT平衡模拟（通常数百皮秒），让缺陷周围的晶格畸变得以弛豫，并让体系体积调整到平衡状态，然后再开始两相模拟。否则，初始的应力集中会严重影响相变行为。

4.2 缺陷如何“催化”熔化？

模拟结果（对应原文Figure S14）清晰地展示了缺陷的显著影响。在相同的压力下（例如10， 30， 60， 90 GPa），随着缺陷浓度的增加，固-液界面向固相一侧移动的速度明显加快。换句话说，在相同的模拟温度和时间内，高缺陷浓度的体系其B1相分数下降得更快。

这意味着，缺陷的存在有效降低了SiC的熔化温度。对于1/100的高缺陷浓度，估算的相变温度Tm比完美晶体下降了约100-200 K（具体数值取决于压力）。这个效应在较高压力下似乎更为明显。

其微观机制可以理解为：缺陷（特别是空位和反位缺陷）破坏了晶格的周期性，引入了局部的应变场和化学键的畸变。这些区域原子的结合能较低，振动更剧烈，相当于在晶体中预先埋下了“薄弱点”或“非均匀形核中心”。当温度升高时，熔化优先从这些缺陷处开始，并向周围扩展，从而降低了整体熔化所需的热力学驱动力，表现为表观熔化温度的下降。

4.3 对材料设计的启示

这一发现对极端环境用SiC材料的设计具有直接指导意义。它表明，材料的纯度、结晶质量（缺陷密度）是其高温高压稳定性的关键控制因素之一。在制备用于超高压环境的SiC部件时，必须通过优化生长工艺（如化学气相沉积CVD的参数）、进行后退火处理等手段，尽可能降低点缺陷浓度。同时，这也提示我们，在利用MD模拟预测材料性能极限时，如果不考虑缺陷，可能会得到过于乐观的（即偏高的）相变温度或强度值。将缺陷效应纳入计算材料设计的考量范畴，是走向“现实材料”模拟的重要一步。

5. 模拟尺寸效应与不确定性评估

5.1 体系大小是否足够？——尺寸收敛性测试

一个常见的质疑是：你用16000个原子模拟的结果，能代表宏观材料吗？有限尺寸效应会不会带来误差？为了回答这个问题，我们进行了严格的尺寸收敛性测试。在30 GPa和60 GPa两个压力点，我们分别构建了64000原子和128000原子的更大两相体系，重复了相同的模拟流程。

结果令人鼓舞（对应原文Figure S15）。对比不同体系尺寸计算得到的相变温度Tm，其差异小于5 K。这个差异远小于我们由温度间隔和界面移动速度估算的Tm不确定度（约±50 K）。因此，我们可以确信，对于研究固-液相变这一现象，16000原子的体系在热力学上已经足够大，其界面能、体积涨落等有限尺寸效应的影响已可忽略不计。

实操心得：进行尺寸收敛性测试是计算工作中体现严谨性的重要环节。选择的关键压力/温度点应有代表性（如相界线中点、拐点等）。更大的体系虽然更“保险”，但计算成本呈立方增长。我们的测试表明，对于两相模拟，确保每相有足够厚的区域（通常每相厚度大于2-3倍原子间相互作用截断半径）比单纯追求总原子数更重要。16000原子的体系在计算精度和效率之间取得了良好平衡，使得我们能够系统地扫描多个压力点。

5.2 如何量化与监控模拟的不确定性？

MD模拟的结果并非确定值，而是存在统计涨落和系统误差。我们采用了一种基于贝叶斯思想的主动学习框架来评估和监控不确定性。简单说，我们不仅使用一个ML势函数，而是使用一组（或一个能给出不确定性估计的）势函数。在模拟过程中，实时监测模型对当前原子构型预测力的不确定性。

我们定义了一个“不确定性指标”，例如原子受力预测的标准差。当这个指标超过一个预设的阈值时（对应原文Figure S16中的Threshold），就意味着模拟进入了当前ML势函数训练数据覆盖不足的构型空间（可能是全新的原子排布方式）。这时，模拟会自动触发第一性原理计算，对这个“陌生”构型进行精确计算，并将新数据加入训练集，更新（或重新训练）势函数。更新后的势函数再继续驱动MD模拟。

这种方法将MD模拟从“开环”变成了“闭环”。结果显示，在我们研究的压力和温度范围内，不确定性指标绝大部分时间都远低于阈值，仅在极少数高温高压的瞬时涨落中触及阈值。这双重验证了：第一，我们使用的ML势函数在整个模拟路径上是可靠的；第二，我们最终报道的相变温度是基于一个高置信度的势函数模型得到的，结果稳健。

6. 相图其他区域的探索与现象讨论

6.1 高压下的固-固相变：B3到B1

除了固-液熔化线，我们也简要探索了SiC在高压下的另一种重要相变：从常压稳定的闪锌矿结构（B3，又称立方SiC）向高压稳定的岩盐结构（B1）的转变。我们采用了压力扫描的方法：在固定温度下（如300K， 1000K， 2000K， 3000K），对B3-SiC体系逐步施加压力，并监测一个反映B1相含量的序参量（例如占据简单立方格点的原子数）。

模拟结果（对应原文Figure S17）��示，随着压力升高，B1相的含量从0开始增加，在某个压力区间内完成转变。相变压力随温度升高而略有增加，这与克莱珀龙-克劳修斯方程定性一致。这个固-固相变过程通常比熔化更快，在几十皮秒内即可完成。通过这种方法，我们可以勾勒出B3-B1的相界线，它与熔化线共同构成了SiC高压相图的主体框架。

6.2 低压下的升华与高压液相的稳定性

一个有趣的问题是：在较低压力下，SiC是直接熔化，还是先升华成气相？为了探究这一点，我们在0-15 GPa的压力范围内，对B3-SiC进行了系列NPT模拟，并大幅提高温度。

观察发现（对应原文Figure S18）：在0， 1， 5 GPa下，当温度超过约4000-4200 K时，模拟盒子的体积开始持续、单调地膨胀，且无法稳定。这表明体系发生了持续的升华，原子逐渐脱离凝聚态，向气相扩散，无法形成一个稳定的液相。而在10和15 GPa下，体系体积经历一个初始的瞬态膨胀后，便稳定在一个新的平衡值。这说明在此压力下，SiC能够形成稳定的凝聚态液体，升华被抑制。

这定性地指出了一个压力阈值（大约在5-10 GPa之间）。低于此阈值，高温下的SiC倾向于升华；高于此阈值，则能形成稳定的液相。这个发现对于理解SiC在低压高温环境（如某些气相沉积过程或烧蚀环境）下的行为非常重要。它提醒我们，材料的相变路径并非一成不变，压力可以根本性地改变其高温下的稳定相。

7. 常见问题、排查技巧与实操建议

7.1 模拟不收敛或结果异常？自检清单

1. 能量和温度持续漂移：

   • 检查：首先检查势函数文件是否加载正确，原子类型是否匹配。然后，在NVE（微正则）系综下运行一个短时间（如10 ps）的模拟，看总能量是否守恒（波动应小于1e-5/原子）。如果能量不守恒，问题很可能出在势函数本身或力/能量的计算设置上。
   • 检查：在NPT/NVT系综下，确保热浴和压浴的弛豫时间参数设置合理。弛豫时间太短会导致系统被过度“拉扯”，产生振荡；太长则控温控压效果差。通常设为时间步长的100-1000倍是一个好的起点。

2. 两相界面移动过快或过慢：

   • 过快：可能当前温度离相变点太远，或者初始液相/固相制备不当（如液相中有晶核）。检查初始液相的RDF，确保是无序液体。可以尝试在更接近预估Tm的温度下进行模拟。
   • 过慢：模拟时间可能不够。固-液界面的移动速度在接近Tm时很慢，可能需要数纳秒才能观察到明显移动。确保模拟时间足够长（通常需要界面移动超过其自身厚度数倍的距离）。同时，检查压力是否准确施加，高压下原子扩散慢，界面移动也可能更慢。

3. 压力张量非对角分量过大：

   • 检查：模拟盒子形状是否允许充分弛豫。在NPT系综下，确保盒子的三个方向以及倾斜度（xy， xz， yz）都是可变的（fix npt命令中设置tri关键词）。运行一段时间的平衡模拟，让盒子形状调整到应力最小的状态。
   • 检查：界面本身可能引入了剪切应力。确保两相界面是低指数晶面（如(100)面），并且界面两侧的晶体取向匹配，以减少晶格失配应力。

7.2 数据分析中的陷阱与技巧

1. 结构识别算法的选择：

   • CNA对完美晶体非常有效，但对高度畸变的结构或液体可能误判。PTM算法更稳健，能区分更多结构类型，但计算量稍大。对于SiC，特别是高压B1相，建议使用PTM或结合多种算法（如CNA+中心对称参数）进行判断，并通过可视化软件（如OVITO）手动抽查确认。

2. 团簇分析中的截断半径：

   • 碳团簇分析中，定义两个碳原子“相连”的截断半径非常关键。设得太小，会漏掉实际成键的原子；设得太大，会把本不连接的原子算成一个团簇。一个可靠的方法是先计算液态或非晶SiC中C-C对的径向分布函数（RDF），将第一个谷值的位置作为截断半径的参考。

3. 相变温度的统计误差：

   • 两相模拟给出的Tm本身存在不确定性，主要来源于界面移动速度的统计涨落和有限的模拟时间。可以通过以下方式估算：在接近Tm的温度进行多次独立模拟（改变随机数种子），统计界面位置随时间变化的斜率，通过线性拟合外推到界面位置为零的点，并计算拟合的误差棒。通常，我们报道的Tm会附带一个±50 K左右的误差范围。

7.3 计算资源与效率优化建议

1. 并行策略：LAMMPS在空间分解上并行效率很高。对于我们的两相体系，由于原子分布均匀，使用-sf opt或-sf gpu（如果有GPU）能获得很好的加速比。确保MPI进程数或GPU卡数与体系在三个方向上的网格划分匹配，避免通信开销过大。

2. 输出频率与磁盘I/O：输出轨迹（dump文件）是主要的I/O瓶颈和存储消耗者。对于长达数纳秒的模拟，不需要每一帧都保存。根据界面移动速度，设置合适的输出间隔（例如每10-50 ps输出一帧）。同时，只输出必要的物理量，如坐标、速度、元素类型。其他衍生量可以在模拟结束后通过LAMMPS的fix ave/time或自行编写作后处理分析。

3. 势函数评估优化：ML势函数的评估通常是计算中最耗时的部分。如果使用LAMMPS的pair_style mliap或通过库（如libtorch）调用PyTorch模型，确保已启用Intel MKL或CUDA加速。对于固定的研究体系，可以考虑将势函数转换为计算更快的表格势或解析形式，但需谨慎测试转换精度。

通过这套从模型构建、模拟运行到数据分析的完整流程，我们不仅获得了SiC在极端条件下的具体相变数据，更建立了一套研究此类问题的可靠方法学。它告诉我们，面对复杂的高温高压相变问题，精心设计的模拟方案、稳健的势函数、严谨的收敛性测试以及对不确定性的量化，是得到可信物理结论的基石。而将缺陷等实际因素纳入考量，则让我们的计算模型向真实的材料世界又迈进了一步。