机器学习势函数在碳化硅极端环境模拟中的应用与验证 1. 项目概述当机器学习“学会”了原子间的“对话”在材料科学的前沿尤其是在核能、航空航天和极端环境电子器件领域碳化硅SiC正扮演着越来越关键的角色。这种材料以其卓越的硬度、高热导率、高击穿电压和优异的抗辐照性能而闻名。然而一个核心的挑战始终横亘在材料研究者面前我们如何精准预测SiC在高温、高压乃至强粒子辐照等极端服役环境下的微观结构演化与性能退化实验手段固然直接但在极端条件下例如数千开尔文、上百吉帕斯卡的压力实时、原位地观测原子尺度的动态过程其难度和成本都高得惊人。这时计算模拟特别是分子动力学MD模拟就成为了我们窥探原子世界的一双“眼睛”。MD模拟通过求解牛顿运动方程追踪每一个原子在给定势函数描述原子间相互作用力下的轨迹从而模拟材料的熔化、相变、缺陷产生与演化等过程。可以说模拟结果的可靠性几乎完全取决于势函数的准确性。传统上我们依赖经验势函数如Tersoff势、ZBL势等。这些势函数基于物理直觉构建了简洁的解析形式计算效率极高能够模拟数百万原子、微秒尺度的系统。但它们的“软肋”也同样明显其函数形式简单参数有限往往是在特定条件下如常温常压的完美晶体拟合的。一旦将模拟场景推向极端条件——比如高压下SiC从闪锌矿结构3C向岩盐结构RS的相变或者高能粒子撞击产生的非平衡态缺陷团簇——传统势函数的预测精度就会急剧下降甚至给出完全错误的物理图像。这就好比用一个在平地上校准的指南针去指引深海或高空航行其可靠性难以保证。近年来机器学习势函数ML-IAP的兴起为这一困境带来了革命性的突破。其核心思想非常巧妙我们不再试图用一个人为设计的简单公式去“猜测”复杂的原子间相互作用而是让机器学习模型直接从高精度的第一性原理如密度泛函理论DFT计算数据中“学习”这种相互作用。DFT虽然精度高但计算量巨大通常只能处理几百个原子、几个皮秒的体系。ML-IAP则充当了一个“桥梁”和“加速器”它通过训练内化了DFT所描述的复杂、高维的能量曲面在保持接近DFT精度的同时将计算速度提升了数个量级使得百万原子规模、微秒时长的模拟成为可能。本项工作的核心正是开发并验证了一个针对Si-C体系的、基于表格化高斯近似势tabGAP框架的机器学习势函数。我们不仅仅满足于“造出一个好用的势函数”更是要用它作为探索未知领域的“探针”去系统性地回答两个关键的科学与工程问题SiC在从常压到极端高压~110 GPa、从室温到极端高温~6500 K的广阔范围内其稳定的物相是什么相变的边界在哪里这直接关系到SiC在超高压设备如金刚石压砧、行星内部碳系外行星或高超音速飞行器热防护系统中的行为。当高能粒子如中子、离子辐照SiC时不同晶体结构如2H和3C的“抗损伤”能力有何差异初始的缺陷是如何产生和演化的这决定了SiC在核反应堆包壳、空间辐射环境电子器件中的使用寿命和可靠性。通过构建一个覆盖了从完美晶体、点缺陷、非晶态到熔融态、甚至纯硅和纯碳相的庞大训练数据库我们训练出的tabGAP势函数在能量和力预测上达到了与DFT高度一致的水平。基于此我们首次通过大规模分子动力学模拟绘制了SiC迄今为止最完整、最细致的压力-温度相图并深入揭示了2H和3C两种重要晶型在辐照下的阈值位移能分布与碰撞级联损伤的微观机制。这些发现不仅澄清了以往实验和模拟中的一些争议更为SiC材料在极端环境下的应用提供了原子尺度的设计依据和理论预言。2. 势函数构建从数据海洋到精准模型的炼金术构建一个优秀的机器学习势函数远非简单地收集数据、跑个训练脚本那么简单。它更像是一场精心策划的“炼金术”目标是将离散的量子力学计算结果炼化成能够连续、快速、准确预测原子间作用的“万能公式”。这里数据库的构建、模型的选择与训练策略每一个环节都至关重要。2.1 训练数据库覆盖“可能性”的战术设计一个ML-IAP的泛化能力即处理训练数据未见过的原子构型的能力直接取决于其训练数据库的广度与质量。对于SiC这样一个多晶型丰富、且在极端条件下可能发生分解、非晶化、熔化等复杂相变的体系我们的数据库设计必须具有前瞻性和系统性。我们的训练数据库包含了3460个原子构型对应185,542个局部原子环境。关键在于这些构型并非随机生成而是被策略性地分为四个互补的“域”以确保模型能应对各种挑战分散构型域基础相互作用校准这个域包含了二聚体Si-Si, C-C, Si-C和三聚体的小分子。它们的作用非常“基础”但关键二聚体主要用于精确描述原子在极短距离下的强排斥作用。这在模拟高能碰撞事件如辐射损伤的初始阶段时至关重要因为原子可能被挤压到非常近的距离。三聚体以较低的计算成本覆盖了各种可能的键角组合。这确保了模型在预测三体相互作用键角弯曲能时不会外推到不合理的能量值避免了在模拟中出现“原子穿模”或异常稳定的畸形结构。晶体SiC多晶型域核心相空间这是数据库的“主力军”涵盖了五种主要的SiC晶型2H、3C、4H、6H和高压岩盐RS结构。但我们没有止步于完美的晶体晶格畸变我们对原胞施加了从-5%到5%的各向同性应变模拟从高拉伸到高压缩的力学环境。有限温度采样通过第一性原理分子动力学AIMD在高温下采样获取原子热振动导致的瞬时构型让模型“见识”到温度带来的无序性。点缺陷引入特别是在3C-SiC中我们引入了弗伦克尔对空位-间隙原子对等点缺陷。这是模拟辐射损伤的基础因为辐照的本质就是大量点缺陷的产生和演化。纯元素相域分解与偏析的预言基础SiC在高温下可能发生分解生成液态硅和固态碳石墨或金刚石。为了准确模拟这一“ incongruent melting”非共晶熔化过程我们必须让模型“认识”纯硅和纯碳的各种形态。因此我们纳入了碳金刚石、石墨、石墨烯、非晶碳。硅金刚石结构硅、液态硅、非晶硅。 这样当模拟中Si-C键断裂后模型才能正确地驱动硅原子聚集成液滴碳原子形成石墨片或金刚石团簇而不是表现出其他非物理行为。无序体相域高温与无序态的考验我们通过将3C-SiC在AIMD中快速加热至3000 K使其熔化再以极高的淬火速率10^14 K/s冷却生成了非晶态SiC。同时在淬火过程中我们对模拟盒子施加了±5%的体积变化对应约±45 GPa的压力范围得到了不同密度下的非晶结构。这迫使模型学习在完全失去长程有序后原子间如何重新成键、配位这对于模拟材料在辐照或冲击加载下的非晶过程至关重要。实操心得数据库的“质量”与“数量”博弈在构建数据库时一个常见的误区是盲目追求构型数量。实际上构型的“代表性”和“多样性”远比单纯的数量更重要。我们采用SOAP平滑原子位置重叠描述符将高维的原子结构映射到二维空间进行可视化如图1a可以直观地检查数据库是否覆盖了所有我们关心的相空间区域各晶型、非晶、液态等。如果发现某些区域空白就需要有针对性地补充采样。此外对于高能构型如二聚体、剧烈畸变的结构我们在训练时会赋予较低的权重因为虽然它们的力很大、对拟合误差贡献显著但它们在平衡态附近出现的概率极低。这种“分层加权”的策略能防止模型为了拟合这些极端点而牺牲了在常见构型上的精度。2.2 模型选择与验证为何是tabGAP在众多ML-IAP框架中如DeePMD、GAP、NequIP等我们选择了表格化高斯近似势tabGAP。这背后有深刻的工程与效率考量计算效率tabGAP的核心思想是将学习到的原子间相互作用“表格化”。在模拟运行时不需要进行复杂的神经网络前向传播而是通过查表结合插值来获取能量和力。这使得其计算速度比许多基于神经网络的势函数快一个数量级以上特别适合我们计划开展的百万原子、微秒尺度的超大规模分子动力学模拟。精度与泛化平衡tabGAP基于高斯过程回归本质上是一种核方法。它在小到中等规模的数据集上就能表现出良好的拟合能力并且通过核函数的设计能提供对预测不确定性的估计。这对于我们判断模拟结果在哪些区域可能外推、可靠性如何是一个很有用的内置诊断工具。与长程作用兼容tabGAP框架天然支持与描述长程库仑作用的屏蔽库仑势相结合。对于SiC这种具有一定离子性的材料这一特性有助于更准确地描述其远距离相互作用。模型的验证是信心的来源。我们进行了多层次、多维度的验证基础拟合精度如图1c, d所示tabGAP预测的能量和力与DFT参考值在非常宽的范围内能量-10到10 eV/原子力-400到400 eV/Å高度一致均方根误差RMSE分别低至0.013 eV/原子和0.679 eV/Å。对于晶体块体误差更是低至0.2 meV/原子和0.007 eV/Å。状态方程图1b展示了五种SiC晶型在零温下的能量-体积曲线。tabGAP的预测实线与DFT数据点方块几乎完全重合并且正确复现了各相在常压下的稳定性顺序4H最稳定其次是6H、3C、2HRS最不稳定。这是势函数能正确预测压力驱动相变的基础。声子谱声子谱直接反映了晶体的动力学稳定性与热学性质。如图2所示tabGAP计算的3C、2H和RS相的声子色散关系与DFT计算结果吻合良好。对于高压下的RS相其声子频率随压力增加而升高的趋势也被正确捕捉。虽然在光学支和RS结构上存在微小差异但这在覆盖极端能量范围的势函数中是预期之内且可接受的。非晶态结构这是对势函数“描述无序能力”的终极考验之一。我们比较了在压缩、无应变、拉伸条件下tabGAP与AIMD模拟得到的非晶SiC的径向分布函数RDF和键角分布图3。两者在关键特征上表现出高度一致性例如在约1.58 Å处出现的C-C键峰清晰地指示了从sp3到sp2杂化的石墨化倾向Si-C键的键长和配位环境也得到了准确描述。这些验证表明tabGAP有能力处理高度无序、偏离理想化学计量比的复杂原子环境。3. 极端条件相图绘制SiC的“生存地图”有了经过严格验证的tabGAP势函数我们便拥有了探索未知P-T疆域的“罗盘”。绘制SiC的压力-温度相图不仅是对势函数预测能力的终极测试更能揭示其在极端环境下丰富的相变与分解行为为理解和设计相关应用提供关键地图。3.1 模拟方法与相边界确定为了构建可靠的相图我们结合了多种热力学计算方法固-液界面法用于直接模拟熔化过程。我们构建一个包含固相和液相的体系在NPT系综下运行MD观察界面随时间的移动方向从而确定在给定压力下的熔点。自由能计算对于固-固相变如3C ↔ RS直接模拟相变动力学极其缓慢。我们采用热力学积分或晶格反演等方法如CALPHAD型计算计算各相在不同P-T条件下的吉布斯自由能自由能最低的相即为稳定相。直接MD模拟对于分解过程我们直接将晶体加热至目标温度观察其是否分解为Si和C以及分解产物的形态石墨还是金刚石。通过大量此类模拟我们将结果绘制于图4a并与已有的实验和DFT数据对比。整个相图可以根据压力划分为四个特征区域其相变序列从低温到高温依次展开。3.2 相图深度解析从升华到高压分解1. 低压区 0.01 MPa升华与表面石墨化在极低压和高温下如10 Pa, 2000 KSiC不再熔化而是直接升华。我们的模拟清晰地捕捉到了这一过程Si原子从表面优先挥发留下富碳的表面层并自组织形成石墨烯结构图4b-i。这与实验上通过SiC表面Si原子升华来制备石墨烯的观测完全一致。当压力稍高1-100 MPa温度极高5000-7000 K时块体SiC会直接解离成单个原子或小分子进入气相区。这个区域定义了SiC在近真空或低压高温环境如某些太空或等离子体环境下的稳定性边界。2. 中压区0.01 MPa – 1 GPa非共晶熔化与石墨析出这是最贴近许多实际应用如常压高温烧结的区域。模拟发现SiC在此区域发生非共晶熔化即它并非直接熔化为均匀的SiC液体而是分解为液态硅和固态碳的混合物。碳以石墨蜂窝状晶格的形式析出图4b-ii。随着温度升高石墨在约3500 K时开始溶解体系逐渐转变为均匀的SiC熔体图4b-iii。此外我们预测在约100 MPa、2000 K时3C-SiC会向6H-SiC发生相变这与化学气相沉积CVD制备SiC薄膜的高温实验观察相符。3. 高压区1 – 60 GPa三阶段分解与金刚石的形成随着压力突破1 GPaSiC的分解行为变得更加复杂呈现三个阶段亚稳态分解阶段在温度刚超过固相线时碳原子开始聚集但体系并未达到全局平衡SiC晶体、液态硅和固态碳短暂共存。稳定固-液混合相温度继续升高SiC完全分解形成液态硅 金刚石的混合物图4b-vi。这是因为在高压下金刚石比石墨更稳定。我们的模拟给出了这一分解边界的预测。均匀液相或相分离液相温度进一步升高金刚石熔化。在压力低于30 GPa时最终形成均匀的SiC熔体图4b-vii。然而当压力高于30 GPa时一个有趣的现象出现了体系难以形成均匀液体而是倾向于发生硅-碳原子的液相分离形成富碳区和富硅区图4b-ix。这表明在超高压下Si-C混合的焓可能为正倾向于分相。4. 超高区 60 GPa岩盐相与极端相分离在低于3000 K的温区当压力超过60 GPa3C-SiC会转变为更致密的岩盐RS结构。我们的自由能计算给出的相边界与实验和DFT结果吻合良好。在更高温度下RS相也会分解其行为与高压区类似但液相分离的趋势更为显著。在一个极端点200 GPa, 8000 K的模拟中我们甚至观察到液态硅开始结晶而碳仍保持液态的奇特景象图4b-这揭示了在极端行星内部条件下可能存在的物态。注意事项模拟与实验的差异解读我们的模拟预测的分解起始温度在高于10 GPa时比部分实验值高出约200 K。这并非一定是势函数的误差更可能源于实验测量的挑战在超过2000 K的高温下红外测温或热电偶测量会因SiC表面的热辐射损失而低估真实温度。此外实验样品中的杂质如氧、金属夹杂物通常会降低观测到的分解温度。因此高精度的模拟有时能为理解实验现象提供更“纯净”的理论参照。4. 辐射损伤模拟从单点缺陷到级联碰撞辐射损伤是核材料和空间应用材料性能退化的核心机制。一个高能粒子初级撞出原子PKA撞入晶格会像保龄球击倒球瓶一样引发一连串的原子碰撞形成所谓的“碰撞级联”最终产生大量空位、间隙原子及其团簇。评估材料抗辐射能力的第一个关键参数就是阈值位移能。4.1 阈值位移能材料抗辐照的“起跑线”阈值位移能TDE定义为将一个晶格原子永久击出其格点位置所需的最小动能。它是一个具有强烈方向依赖性的物理量。我们通过两种方法进行验证和计算准静态拖动计算如图5所示我们选择一个原子Si或C沿着特定的晶体学方向如[111]缓慢拖动同时用DFT和tabGAP计算体系总能的变化。曲线上的峰值代表需要克服的能量势垒。结果显示tabGAP的能量曲线与DFT高度吻合远优于传统的Tersoff/ZBL混合势。大规模方向采样为了获得全面的TDE分布我们对每个PKA类型Si或C在2H和3C-SiC中随机选取了约7000-8000个方向进行完整的MD模拟统计产生永久位移所需的最小能量。计算得到的TDE分布图图6揭示了深刻的信息总体趋势无论是tabGAP还是Tersoff/ZBL都预测Si原子的TDE普遍高于C原子。这是因为Si原子更重与周围原子的结合更强。势函数差异Tersoff/ZBL势普遍给出了更宽、平均值更高的TDE分布。尤其在3C-SiC中它对Si和C的TDE差值预测是tabGAP的两倍之大。这种差异会显著影响后续碰撞级联中缺陷产生的数量和类型。方向各向异性TDE地图图6a,b中的热图清晰展示了强烈的方向依赖性。例如在2H-SiC中C原子的高TDE60 eV仅出现在如[111]等少数特定方向。在3C-SiC中对于Si的PKATersoff/ZBL预测沿开方向如[1-1-1]和闭方向如[111]的TDE存在明显不对称性而tabGAP预测两者都在20 eV左右。文献中基于AIMD的研究对此也有不同报道这可能与TDE的具体定义有关是要求PKA本身永久位移还是只要有任何原子永久位移。我们的tabGAP结果支持后一种定义下不对称性减弱的观点。4.2 碰撞级联模拟揭示多晶型依赖的损伤机制TDE是入门券而真实的辐射损伤是一个动态的、包含数百万原子、皮秒到纳秒尺度的非平衡过程。我们利用tabGAP的高效性模拟了不同初始能量的PKA在2H和3C-SiC中引发的完整碰撞级联。模拟揭示了传统势函数难以捕捉的关键现象缺陷团簇的多晶型依赖性。在3C-SiC中碰撞产生的空位和间隙原子倾向于快速扩散并湮灭或者形成小而分散的团簇。这得益于其立方对称的结构为缺陷迁移提供了更多等效的路径。而在2H-SiC中由于其六方结构带来的各向异性缺陷特别是间隙原子的迁移受到更多限制。它们更容易在碰撞区域附近聚集形成更大、更稳定的缺陷团簇如双空位、三空位或者间隙原子环。这种差异至关重要。大而稳定的缺陷团簇是材料中位错环、空洞等扩展缺陷的“种子”会严重劣化材料的力学性能和尺寸稳定性。因此2H-SiC在抗辐照肿胀方面可能天生就比3C-SiC更脆弱。传统的Tersoff势由于无法准确描述不同晶型下缺陷形成的细微能量差异很可能低估或错误预测这种团簇行为。实操心得大规模辐射模拟的参数设置进行百万原子规模的碰撞级联模拟时计算参数的选择直接影响结果的物理意义和计算效率盒子尺寸与边界条件模拟盒子必须足够大以确保级联产生的冲击波在到达边界前充分衰减避免周期性边界条件带来的自相互作用。通常盒子边长至少是级联区域估计尺寸的3-5倍。我们使用自由表面或固定边界层来吸收能量。系综选择级联初期~0.1 ps能量高度局域应采用微正则系综NVE以保持能量守恒。在热峰冷却和缺陷演化阶段~10-100 ps应切换到等温系综NVT或等温等压系综NPT使用朗之万或Nosé-Hoover热浴将多余的热量耗散掉模拟真实材料中的热传导。时间步长由于涉及高能碰撞初期必须使用非常小的时间步长如0.1 fs。在热峰冷却后可以逐步增大到0.5-1.0 fs以提高效率。需要编写脚本自动监测体系温度或最大原子速度动态调整步长。结果分析缺陷识别常用Wigner-Seitz元胞法或邻域分析法。不仅要统计点缺陷数量更要分析它们的空间分布、团簇尺寸分布和类型。可视化软件如OVITO是理解复杂级联形态不可或缺的工具。5. 工具链、复现指南与未来展望5.1 核心工具与工作流本项目的研究高度依赖于一套成熟的开源计算软件生态复现或开展类似工作可遵循以下路径第一性原理计算与数据生成VASP / Quantum ESPRESSO / ABINIT用于生成训练数据库所需的DFT能量、力和应力数据。需要精心设计采样策略覆盖目标P-T空间和结构空间。LAMMPS (AIMD)也可用于生成经典势函数下的初始构型再通过DFT进行单点能校正这是一种高效的主动学习策略。机器学习势函数训练QUIP/GAP本工作使用的tabGAP框架集成在QUIPQuantum Mechanics and Interatomic Potentials软件包中。它提供了完整的GAP模型训练和部署接口。DeePMD-kit另一个广泛使用的深度机器学习势框架基于深度学习在精度和效率上也有出色表现尤其适合非常复杂的多元素体系。训练流程1) 准备DFT计算数据2) 提取或计算原子环境描述符如SOAP、ACSF3) 使用上述软件包进行模型训练与超参数优化4) 在独立的测试集上验证模型精度。大规模分子动力学模拟LAMMPS业界标准的MD模拟软件支持多种ML-IAP的接口如通过pair_style quip或pair_style deepmd调用。我们所有的大规模相变和辐射损伤模拟均在LAMMPS中完成。关键输入需要编写详细的LAMMPS输入脚本定义盒子、创建晶体、设置势函数、定义系综、施加压力/温度、引入PKA等。数据分析与可视化OVITO原子模拟可视化的首选工具用于观察原子构型、相变过程、缺陷分布等。Python (NumPy, SciPy, Matplotlib, ASE)用于自动化数据处理、计算RDF、键角分布、分析缺陷、绘制相图和TDE地图等。ASE原子模拟环境库提供了与多种计算代码交互的通用接口。5.2 常见问题与排查技巧训练误差低但模拟结果发散可能原因数据库缺乏对高能非平衡态如碰撞级联初期的覆盖描述符的截断半径设置过小未能捕捉到足够的原子环境信息训练时对高能量构型的正则化权重设置不当。排查检查模拟发散瞬间的原子构型计算其SOAP描述符并与训练数据库的SOAP分布进行对比看是否处于外推区域。可以针对这些“问题构型”补充DFT计算加入训练集重新训练主动学习。相变温度/压力与实验偏差大可能原因DFT泛函本身存在系统误差如GGA通常低估相变压力训练数据未充分包含相变路径附近的过渡态构型模拟尺寸有限存在有限尺寸效应。排查使用更高级的DFT方法如杂化泛函、RPA对关键点进行校验。采用更大的超胞进行相变模拟或使用专门计算相界的方法如固-液界面法、热力学积分。辐射模拟中缺陷数量异常可能原因TDE计算不准确导致PKA初始动能设置不合理模拟时间不够长缺陷未充分演化势函数对点缺陷的形成能描述有误。排查务必先进行全面的TDE方向采样获得可靠的阈值。延长模拟时间观察缺陷数量是否趋于稳定。计算空位、间隙原子的形成能与DFT或实验值对比。计算速度慢可能原因ML-IAP本身比经典势慢描述符计算是瓶颈未使用GPU加速或并行效率低。优化对于tabGAP确保使用了其高效的表格查找模式。对于DeePMD启用GPU支持可带来数十倍加速。检查LAMMPS的并行设置确保域分解与通信开销合理。5.3 未来拓展与应用前景这项工作开发的tabGAP势函数及其揭示的物理规律为SiC及相关领域的研究打开了新的大门势函数的扩展与迁移当前的势函数专注于Si-C二元体系。一个自然的扩展是加入氧O、氢H等元素用于研究SiC在氧化、腐蚀环境下的行为或作为聚变堆中面向等离子体材料时的氢滞留问题。这需要构建包含Si-C-O/H的更大训练数据库。多尺度模拟桥梁ML-IAP的精度和效率使其成为连接第一性原理与宏观连续介质模型如有限元分析的理想桥梁。我们可以用ML-IAP模拟纳米尺度裂纹扩展、位错运动提取关键参数如表面能、位错芯结构、Peierls应力传递给更高尺度的模型实现从原子到器件的跨尺度预测。材料设计与筛选这套方法论可以推广到其他重要的极端环境材料如氮化硼BN、碳化锆ZrC、钨W等。通过高通量生成训练数据、自动化训练流程可以快速构建一系列材料的精准ML-IAP用于筛选在特定辐照、高温、高压下性能最优的材料。动力学过程深入研究利用此势函数可以更深入地研究相变和损伤的微观动力学路径。例如高压下3C到RS相变的原子集体运动模式是什么碰撞级联中是哪些特定的原子碰撞序列导致了2H-SiC中更大的缺陷团簇这些原子尺度的机理研究是设计抗辐照、相变增韧新材料的基础。机器学习势函数正在从根本上改变计算材料科学的研究范式。它不再仅仅是一个计算工具而是一个能够从数据中学习并发现新物理的“合作研究者”。对于像SiC这样在国计民生关键领域扮演核心角色的材料拥有一个能在极端条件下可靠工作的原子尺度“数字替身”其价值不言而喻。它让我们能够在计算机中安全、经济、反复地“测试”材料的极限从而加速高性能SiC材料的研发、优化与工程应用。