机器学习势函数在辐射损伤模拟中的性能评估与优化策略

1. 项目概述：为什么势函数是辐射损伤模拟的“命门”

搞材料模拟，尤其是辐射损伤这块的同行，估计都深有体会：分子动力学模拟的结果准不准，八成得看你的原子间势函数靠不靠谱。这玩意儿说白了就是一套描述原子之间怎么“推拉拽扯”的数学公式，是模拟的物理基础。你用它来算晶格振动、扩散这些近平衡态的性质，可能感觉还行；但一旦涉及到高能粒子轰击、原子被撞飞、产生大量缺陷这种极端非平衡过程，很多势函数就开始“露怯”了。最近几年，机器学习势函数火得不行，号称能逼近第一性原理计算的精度，又比直接做第一性原理分子动力学快好几个数量级，听起来简直是辐射损伤模拟的“梦中情势”。但问题是，它真的那么神吗？在不同场景下表现是否一致？这就是我们这次要深扒的核心。

我们这次就拿金属镍开刀。镍是核反应堆结构材料、聚变堆面向等离子体材料的关键组分，理解它在辐射下的行为至关重要。手头有几套势函数：两套基于不同赝势（Ni和Ni pv）训练的表格化高斯近似势，也就是tabGAP，这是机器学习势的代表；还有两套经典的嵌入原子法势，分别是Stoller和Bonny势。我们的目标很明确，就是把这些势函数扔到辐射损伤模拟的“高压锅”里煮一煮，看看它们在阈值位移能、单次级联碰撞、以及最要命的重叠级联累积损伤这些关键测试中，到底表现如何。结果发现，水还挺深，不是所有机器学习势都能在所有环节“通吃”，短程排斥相互作用这个细节，成了决定模拟结果走向的关键胜负手。

2. 核心思路：从静态测试到动态“实战”的评估体系

评估一个势函数在辐射损伤模拟中的适用性，不能只看它算出来的晶格常数、弹性模量准不准。那只是“体检”的常规项目。真正的考验在于它能否准确描述从两个原子无限接近时的剧烈排斥（短程），到正常晶格距离的平衡相互作用（中程），再到缺陷形成、迁移等涉及能量变化的复杂过程。我们的评估逻辑是层层递进的，从简单到复杂，从静态到动态，构建一个完整的“压力测试”链条。

2.1 评估维度的逻辑拆解

首先，最基础的检查是看势函数对原子间相互作用的描述是否“健康”。我们画了二聚体曲线，就是计算两个镍原子在不同距离下的总能量。这能直观看出短程排斥部分的“硬度”。如果这里就和第一性原理计算的结果差很远，那后续的高能碰撞模拟基本可以不用看了，结果肯定跑偏。紧接着，我们做了准静态拖拽计算，模拟一个原子在晶格中沿着特定方向被缓慢拖开时，周围晶格原子的弛豫和能量变化。这个测试比二聚体曲线更进了一步，它包含了晶格环境的影响，能反映势函数对点缺陷形成能的预测趋势。

通过了这些“静态”或“准静态”测试，才能进入真正的动态模拟环节。第一步是计算阈值位移能。这是辐射损伤领域的一个核心参数，指的是将一个晶格原子永久撞离其位置所需的最小动能。我们不是只算几个高对称方向就完事，而是系统性地计算了不同入射方向上的阈值位移能，生成一个三维的“能量球面”。这能全面评估势函数对各向异性效应的预测能力。

然后就是模拟真实的辐射事件——初级离位级联。用一个初始动能为几千电子伏特的初级撞出原子去轰击晶格，用分子动力学模拟整个级联过程：能量沉积、热峰形成、冷却、以及最终稳定下来的缺陷（空位和间隙原子对，即Frenkel对）数量。这一步看的是势函数在极端非平衡、高能条件下的瞬时响应能力。

最后，也是最接近真实辐照环境的测试：重叠级联模拟。在现实中，材料是持续受到辐照的，新的碰撞级联会不断在已经存在损伤的区域产生。我们模拟了成千上万次级联事件在同一个模拟盒子中依次发生，观察缺陷浓度是如何随着辐照剂量累积的，以及最终会形成什么样的缺陷结构（比如位错环、层错四面体等）。这个测试考察的是势函数在长时间尺度、复杂缺陷相互作用下的表现，是对其综合预测能力的终极考验。

2.2 势函数选型背后的考量

为什么选这四套势函数？这里头有对比的巧思。两套tabGAP势函数（Ni和Ni pv）的核心区别在于训练它们所用的第一性原理数据，采用了不同的赝势来处理镍的电子结构。Ni pv赝势 explicitly 处理了镍的3p半芯电子，而标准的Ni赝势没有。这个看似技术性的选择，会直接影响对电子云排布，尤其是内层电子在原子非常接近时的描述，从而显著影响短程排斥力。用它们俩对比，就是想剥离出“机器学习框架相同，但底层量子力学输入不同”所带来的影响。

而Stoller和Bonny这两套EAM势，是经过多年发展、在辐射损伤领域被广泛使用和引用的经验势。它们代表了在机器学习势流行之前，我们所能获得的、经过一定校验的“传统智慧”。尤其是Bonny势，它在描述镍及其合金的辐照肿胀等方面有一定声誉。将这组传统势函数与新的机器学习势进行对比，不仅能判断机器学习势是否真的超越了前辈，更能揭示不同势函数在物理描述上的根本差异点在哪里。

注意：势函数的评估一定要有参照系。单纯说“我的势函数预测缺陷浓度是X”没有意义，必须和可靠的实验数据、更高级别的第一性原理计算结果，或者其他经过广泛验证的势函数进行对比。本次工作中，我们尽可能地将模拟结果与已有的实验阈值位移能数据、以及第一性原理分子动力学的结果进行交叉验证。

3. 关键测试一：短程相互作用与阈值位移能

阈值位移能是辐射损伤研究的“敲门砖”，它直接决定了有多少初始动能会转化为稳定的原子离位。如果势函数连这个基本参数都预测不准，后续的级联模拟就像建立在沙滩上的城堡。

3.1 二聚体曲线与准静态拖拽：短程排斥的“体检报告”

我们先从最基础的二聚体曲线看起。把两个镍原子从很远拉到几乎重叠，计算总能量变化。理想情况下，这条曲线在短程部分应该与第一性原理计算的结果高度吻合。我们的结果显示，两套tabGAP势函数在短程部分出现了肉眼可见的分歧。基于Ni pv赝势训练的tabGAP势，其排斥“墙”更“硬”，即随着距离减小，能量上升得更陡峭；而基于标准Ni赝势的tabGAP势则相对“软”一些。这个差异在二聚体距离小于0.5埃时变得非常明显。相比之下，Stoller的EAM势与Ni pv tabGAP的短程行为更接近，而Bonny势则展现出截然不同的、更“硬”的排斥特征。

准静态拖拽测试进一步放大了这种差异。当我们模拟一个原子沿<100>、<110>、<111>等高对称方向被缓慢拖离时，需要克服的能垒高度直接反映了势函数对晶格恢复力的描述。从提供的能量差图可以清晰看到，不同势函数给出的能量曲线形状和峰值差异显著。例如，在<111>方向，某些势函数预测的拖拽能峰可高达120 eV以上，而另一些则低于80 eV。这直接预示着它们在阈值位移能上会有巨大差别。

实操心得：不要只看势函数论文里报告的平衡态性质拟合误差。一定要亲自画一下关键方向的二聚体曲线和拖拽曲线，与可靠的DFT或更高精度的计算结果对比。短程部分的微小差异，在级联模拟中会被指数级放大。很多现成的势函数文件可能未充分优化这一区域，需要使用者保持警惕。

3.2 阈值位移能的方向各向异性与平均值

我们系统计算了数千个不同入射方向上的阈值位移能，绘制成三维极坐标图。图像非常直观地展示了镍中TDE的强烈各向异性：沿着密排方向（如<110>）通常需要更低的能量就能造成离位，而沿着开放方向（如<111>）则需要更高的能量。这是由晶格通道效应决定的——原子在密排方向运动时，更容易找到“缝隙”穿过而不引起巨大的晶格畸变。

具体到数值，不同势函数的预测可谓五花八门。对于<100>方向，实验值大约在38 eV左右。我们的tabGAP势预测值在21-23 eV，明显偏低；Stoller势预测27 eV，Bonny势高达49 eV。在<110>方向，实验值约21 eV，tabGAP势在19-25 eV区间，Stoller势23 eV，Bonny势又是49 eV。最夸张的是<111>方向，实验值认为大于60 eV，tabGAP势预测29-39 eV，Stoller势35 eV，而Bonny势达到了惊人的69 eV，其平均值更是接近100 eV，远超其他势函数和实验认知。

这些差异说明了什么？首先，Bonny势在几乎所有方向都给出了异常高的TDE，这可能源于其势函数形式中人为加强了原子间的束缚力，以提高某些平衡态性质（如空位形成能）的拟合精度，但却牺牲了非平衡碰撞过程的准确性。其次，两套tabGAP势虽然基于相同的机器学习框架，但由于训练数据源（赝势）不同，在TDE预测上也出现了不可忽视的差异，尤其是在<111>方向。这凸显了训练数据的量子力学计算精度对最终势函数在极端条件下性能的决定性影响。

表1：不同势函数预测的镍阈值位移能对比（单位：eV）

4. 关键测试二：单次初级离位级联模拟

阈值位移能是“门槛”，而单次级联模拟则是看“进门之后”会发生什么。我们模拟了从1 keV到10 keV不同初始动能的初级撞出原子引发的级联过程。

4.1 缺陷产生数量与PKA能量的关系

模拟50次取平均后，我们绘制了稳定Frenkel对数量随PKA能量变化的曲线。一个有趣的发现是：缺陷产生数量并不与平均阈值位移能简单相关。例如，平均TDE最低的Ni tabGAP势，在几乎所有PKA能量下，预测产生的空位数都是最高的。而平均TDE最高的Bonny势，产生的缺陷数量却并非最低，在某些能量区间甚至超过了Stoller势。

这打破了“TDE越低，产生缺陷越多”的简单直觉。原因在于，缺陷最终数量不仅取决于初始离位是否容易发生（TDE），更取决于级联发展过程中，热峰的能量耗散、原子混合、以及后续的空位-间隙原子复合效率。一个“较软”的势函数（如Ni tabGAP），其原子间排斥力较弱，可能导致级联区域原子的碰撞截面更大，能量沉积更分散，形成更弥散、非球形的热峰。这种热峰冷却较慢，原子有更多时间扩散和重组，反而可能抑制了近距离的空位-间隙原子复合，导致最终保留下来的缺陷更多。

4.2 级联时空演化与热峰特征

通过分析5 keV PKA下缺陷数量随时间演化的曲线，我们获得了更动态的图景。所有势函数都展示了典型的级联发展三阶段：快速损伤产生期（~0.1 ps内）、热峰冷却与复合期（~0.1-1 ps）、以及缺陷稳定期（>1 ps）。然而，不同势函数在峰值损伤（热峰最大规模）和最终稳定缺陷数量上的排序并不一致。

Ni tabGAP势的峰值损伤反而是最低的，但最终缺陷残留最多。这印证了上面的推测：它的级联更“软”、更弥散，热峰温度可能相对较低，但范围广，在冷却过程中复合不充分。而Ni pv tabGAP和Stoller势则表现出更高的峰值损伤和相对紧凑的热峰，但最终残留缺陷却较少，说明其热峰冷却更快，或者原子在热峰内更容易复合。

排查技巧：当比较不同势函数的级联结果时，不要只看最终缺陷数。一定要分析级联的时空演化过程，观察热峰的大小、形状、寿命。可以使用原子温度云图、径向分布函数随时间变化等工具。有时，最终缺陷数量相近的两个势函数，其级联的微观物理过程可能截然不同，这会影响对后续现象（如缺陷迁移、聚集）的预测。

5. 关键测试三：重叠级联与长期缺陷演化

单次级联像是“闪电战”，而重叠级联模拟则是“持久战”，它更接近材料在反应堆内经受长期辐照的真实情况。我们在一个模拟晶胞中，依次注入数千个5 keV的PKA，让新产生的级联与之前残留的缺陷结构发生相互作用。

5.1 缺陷浓度的饱和现象

随着注入级联数量（对应辐照剂量）的增加，所有势函数模拟的体系其缺陷浓度都逐渐趋于饱和。这是一个典型的辐照损伤现象：新缺陷的产生与旧缺陷的恢复（复合、湮灭、聚集形成大缺陷）达到了动态平衡。然而，令人震惊的是，达到饱和的缺陷浓度水平，在不同势函数间存在巨大差异。

Ni pv tabGAP势预测的饱和缺陷浓度显著低于其他三种势函数。即使我们将模拟剂量加倍（从2000次级联增加到4000次），它的缺陷浓度依然稳定在低水平，没有向其他势函数的结果靠拢。而其他三种势函数（Ni tabGAP， Stoller， Bonny）的饱和浓度则彼此接近。这个结果非常反直觉，因为在前面的单次级联测试中，Ni pv tabGAP和Stoller势的行为颇为相似，但在累积损伤的长期竞争中，它们却分道扬镳。

5.2 最终缺陷结构的微观分析

观察2000次级联后的原子构型，所有势函数都预测形成了由肖克利不全位错链构成的大尺寸间隙型位错环，以及由空位聚集而成的层错四面体。这是面心立方金属中典型的辐照缺陷结构。定性上看，缺陷类型相似。但定量上，Ni pv tabGAP势体系中的缺陷密度明显更低，大尺寸的位错环结构似乎也更“稀疏”。

这说明，决定长期辐照损伤演化的，不仅仅是单次级联产生缺陷的“数量”，更是这些缺陷的“命运”。缺陷的稳定性、迁移能力、以及它们之间相互作用的势垒，这些都由势函数在近平衡区域（即缺陷周围弛豫后的原子构型）的能量面精确性所决定。Ni pv tabGAP势可能更准确地描述了间隙原子团簇的迁移势垒或空位团簇的稳定性，使得缺陷更容易在后续级联的扰动下复合或发生改变，从而抑制了缺陷的累积。

重要提示：重叠级联模拟是检验势函数预测长期辐照行为的“试金石”，但计算成本极高。在进行此类模拟前，务必用单次级联、阈值位移能、点缺陷性质等测试对势函数进行初步筛选。如果势函数在这些基础测试中表现怪异，那么它在重叠级联中给出错误结果的可能性极大。不要盲目相信一个势函数仅仅因为它标榜“机器学习”或“基于DFT”。

6. 关键���试四：与实验对比的RBS/C模拟

模拟的终极目标是解释和预测实验。卢瑟福背散射/沟道谱技术是实验上表征近表面区域晶格损伤的常用手段。我们通过模拟计算了经过重叠级联损伤后的镍晶体在不同深度处的RBS/C谱，并与实验数据直接对比。

6.1 模拟与实验谱的匹配度

我们将不同势函数产生的最终缺陷结构作为输入，计算了相应的模拟RBS/C谱。对比发现，在采用了更合理的核能量沉积剖面进行采样（而非简单假设一个统一的TDE值）后，Ni pv tabGAP势模拟出的谱线与实验测量结果最为接近。其他势函数，包括另一套tabGAP和两套EAM势，模拟出的背散射产额均显著高于实验值。

这意味着，Ni pv tabGAP势所预测的缺陷类型、浓度及其空间分布，综合起来更符合实验观测到的晶格无序度。这是一个强有力的证据，表明在镍的案例中，采用更精确的赝势（Ni pv）来生成训练数据，所获得的机器学习势函数，在预测真实辐照损伤累积方面，可能具有更高的保真度。

6.2 采样方法的影响

我们对比了两种从模拟中提取结构用于RBS/C计算的方法：一种是传统的基于NRT-dpa模型并假设固定TDE（如40 eV）的方法；另一种是直接基于模拟中实际沉积的核能量剖面进行采样。结果表明，后一种方法虽然增加了计算复杂度，但能更真实地反映不同势函数在TDE上的差异对最终损伤剂量估算的影响，从而使模拟与实验的对比更具物理意义。忽略这种差异，可能会错误地归因或掩盖势函数之间的性能差别。

7. 势函数优化尝试：调整短程排斥相互作用

既然发现了短程排斥相互作用是导致不同势函数行为差异的关键，一个很自然的想法是：能否对训练好的机器学习势函数进行“微调”，专门优化其短程部分？由于tabGAP势函数形式明确分离了对相互作用项和多体项，这为我们进行此类操作提供了便利。

我们尝试了三种修改方案：1) 对“较软”的Ni tabGAP势，用其原始外部排斥势重新拟合其两体项；2) 将Ni tabGAP势的两体项在短程部分平滑过渡到Ni pv tabGAP势的“较硬”两体项；3) 反向操作，将Ni pv tabGAP势的两体项在短程部分过渡到Ni tabGAP的“较软”形式。

7.1 修改后的效果验证

将这些修改后的势函数重新投入重叠级联模拟，结果颇具启发性。对于原版Ni tabGAP势（缺陷产生多），通过方法1和方法2使其短程部分“硬化”后，其预测的饱和缺陷浓度显著下降，变得与未修改的Ni pv tabGAP势结果非常接近。反之，将Ni pv tabGAP势“软化”后，其缺陷浓度有所上升，但仍远低于原版Ni tabGAP势。

这系列实验有力地证明了两点：首先，短程排斥相互作用确实是控制级联损伤演化，特别是最终缺陷累积浓度的关键杠杆之一。其次，对于表格化的机器学习势，在训练后对其两体排斥部分进行有针对性的修正或替换，在技术上是可行的，并且能显著改变其在非平衡模拟中的行为。

7.2 优化策略的局限与展望

然而，这种“外科手术式”的修改必须极其谨慎。我们只是简单地基于二聚体曲线的视觉匹配进行平滑过渡，缺乏系统的优化。修改后的势函数虽然在重叠级联测试中表现改变，但其在点缺陷形成能、表面能、弹性常数等平衡性质上是否依然保持精度，需要重新全面验证。粗暴的修改可能会破坏势函数在训练数据域内的内在一致性。

一个更稳健的思路是将短程排斥的准确性作为目标，从一开始就纳入机器学习势的训练过程中。例如，在训练数据集里，不仅包含平衡构型、弹性变形、点缺陷、声子等信息，还应 explicitly 加入高能二聚体、小间距晶格压缩、甚至是从第一性原理分子动力学级联模拟中提取的高能非平衡原子构型。让机器学习模型在训练阶段就“见识”并学习到这些极端情况下的相互作用，从而生成一个在宽能量范围内都可靠的势函数。

8. 总结与对从业者的建议

通过这一系列从简到繁、从静态到动态的测试，我们对机器学习势函数在镍辐射损伤模拟中的性能有了更立体和深刻的认识。核心结论是：一个适用于辐射损伤模拟的优秀势函数，必须同时在两个“战场”上表现优异——既要能精确描述平衡态和近平衡态的物理性质（如晶格常数、弹性模量、缺陷形成能），也必须能准确刻画短程高能排斥相互作用。

对于从事相关模拟的研究者，我个人的体会是：

1. 切勿迷信“黑箱”：不要因为一个势函数标榜“机器学习”或“第一性原理精度”就全盘接受。务必进行系统的基准测试，特别是针对你的研究问题所涉及的关键物理过程。对于辐射损伤，二聚体曲线、阈值位移能计算和单次级联模拟是最低限度的验证。
2. 理解训练数据的“基因”：仔细阅读势函数原始文献，弄清楚其训练数据集是如何构建的。使用了哪种赝势？是否包含了高能数据？训练集对你想研究的性质（如沿特定方向的离位）是否有足够的覆盖？本次工作中Ni和Ni pv赝势带来的差异就是活生生的例子。
3. 长期行为是终极裁判：对于研究辐照损伤累积、微结构演化等长期效应，单次级联的结果可能具有误导性。只要计算资源允许，尽量进行重叠级联模拟，或者至少研究缺陷团簇的稳定性与迁移率。势函数在长期模拟中表现出的差异，往往揭示了其在描述复杂缺陷相互作用能垒上的根本不同。
4. 实验对比是金标准：尽可能地将模拟结果与实验观测进行定量对比，如RBS/C谱、透射电镜观测的缺陷密度与类型、电阻率变化等。即使无法完全吻合，趋势上的一致性能极大增强模拟结果的可信度。
5. 审慎进行势函数修改：如果确实需要对现有势函数进行修改以改进其短程行为，必须采用系统的方法，并在修改后对势函数的全面性能进行重新评估。记录下所有的修改步骤和参数，确保工作的可重复性。

机器学习势函数无疑为材料辐照模拟带来了新的强大工具，但它并非“银弹”。它放大了我们对高精度训练数据的需求，也要求我们以更严谨、更系统的方式去验证和解读模拟结果。这条路没有捷径，唯有通过细致扎实的基准测试和物理分析，才能让这些强大的计算工具真正服务于对材料辐照行为更深刻、更准确的理解。