朴正贤（JongHyeon Park）|小田卓二（Takuji Oda）
韩国首尔国立大学核工程系，冠岳区冠岳路1号，邮编08826
**摘要**
多尺度建模被广泛用于研究不同时间和长度尺度下材料的辐射损伤过程和效应。对于初级损伤模拟，通常采用二元碰撞近似（BCA）和分子动力学（MD）方法。尽管BCA在计算上效率高，但在预测详细的辐射诱导缺陷结构方面缺乏MD的准确性。在这项工作中，我们开发了一种名为PRDgen的机器学习方法/代码，该代码基于MD模拟数据进行训练，能够利用反冲能量和方向作为输入来生成初级损伤结构。PRDgen由两部分组成：一部分是基于高斯过程回归的缺陷数量和空间分布统计的准随机缺陷结构生成器；另一部分是基于卷积神经网络的辐射损伤结构分类器。这两部分都使用MD反冲模拟结果进行训练。我们系统地评估了不同的反冲能量和方向采样策略，以确定在使用有限MD数据时仍能获得足够准确性的训练数据集。结果表明，PRDgen的速度大约是MD的75-200倍，同时能够再现MD模拟中观察到的几个关键结构特征，包括点缺陷的数量和空间分布。然而，由于缺陷结构生成器中缺乏描述缺陷间关联的统计模型，因此在缺陷聚类方面的再现效果仍然不够理想。初步测试表明，通过考虑缺陷间的关联，效率和保真度都可以显著提高。这些结果表明，PRDgen框架为多尺度建模工作流程提供了实用且高效的方法，以生成初级辐射损伤结构。
**1. 引言**
核反应堆中使用的材料会受到中子辐照的损伤，并可能出现诸如辐照脆化和孔洞膨胀等不良现象[1] [2]。这些材料退化会影响核能系统中结构和功能材料的性能[3] [4] [5]。因此，了解辐射缺陷的形成和发展方式及其对材料性能的影响至关重要。
在晶体材料中，初级撞击原子（PKAs）通过与中子的碰撞获得额外的动能而产生。如果PKA接收到足够高的能量，由于一系列碰撞中的能量传递，会生成许多缺陷，这种现象称为碰撞级联[6] [7]。在碰撞级联开始后的几十皮秒内，所有反冲离子的速度都会降至热振动能量的水平，只留下一些形成的空位和自间隙原子（SIA）[6] [7] [8] [9]。随后，剩余的缺陷会通过扩散、聚集、解离和沉降吸收等热过程在更长的时间尺度上进一步演化。最终，微观结构的变化（通常会导致材料性能的变化）会在更长的时间尺度上发生，例如几天或几年。由于这种广泛的时间和长度尺度范围，有必要适当地结合不同的理论方法，即多尺度建模，以准确预测辐射损伤现象和效应[9]。
多尺度建模的第一步通常是估计碰撞级联后剩余的辐射缺陷的数量和结构。这项任务可以通过几种计算和理论方法完成[9]，例如二元碰撞近似（BCA）[10]、使用分子动力学的反冲模拟[6] [9]，以及使用分子动力学的Frenkel对积累（FPA）方法[11] [12]。反冲MD模拟是一种计算密集型方法，适用于模拟原子的多体碰撞动力学，并能提供损伤过程和形成的缺陷结构的所有细节。大量的研究利用MD模拟来研究辐射诱导的缺陷[13] [14]。另一方面，在BCA中，材料中原子之间的碰撞被近似视为一系列独立的二元碰撞，从而大大降低了计算成本。对于每次碰撞，计算经典散射积分以确定散射角和能量损失[9] [15] [16]。BCA被认为是一种比MD更快的计算离子范围、能量沉积剖面和反冲分布的方法。然而，BCA无法准确提供缺陷结构，因为它只考虑了两个原子之间的短程排斥作用。此外，BCA还被认为会高估产生的缺陷数量[17]。与此同时，FPA通过在MD模拟中以一定的距离约束原子位移来随机生成Frenkel对，也可以显著降低计算成本，允许模拟高达高剂量的损伤积累。然而，由于缺乏描述缺陷结构中缺陷间关联的统计模型，缺陷聚集的再现效果仍然不足。初步测试表明，通过考虑缺陷间关联，效率和保真度都可以显著提高。这些结果表明，PRDgen框架为多尺度建模工作流程提供了一种实用且高效的方法，用于生成初级辐射损伤结构。
**2. 方法**
**2.1. 所提出的方法概述**
如图1所示，所提出的方法使用了两种机器：一种使用高斯过程回归（GPR）的候选缺陷结构生成器（称为机器A），另一种使用卷积神经网络（CNN）的辐射损伤结构分类器（称为机器B）。首先，机器A根据给定的反冲能量和方向生成辐射损伤的候选结构。然后，将生成的结构提供给机器B，机器B会赋予一个分类分数，表明该结构是否类似于MD得到的初级损伤数据。接着应用一个接受规则对该分数进行判断。如果候选结构被接受，它将在给定的反冲条件下作为PRDgen输出返回；否则，机器A会生成另一个候选结构，重复这一过程。
**下载：** 下载高分辨率图像（49KB）
**下载：** 下载完整尺寸图像
**图1. 生成给定反冲能量和方向下的辐射缺陷结构的机器构造和使用概述。** 机器A和B分别使用GPR和CNN深度学习， 并使用MD反冲模拟结果进行训练。通过结合这两种机器来生成辐射缺陷结构。
机器A使用GPR学习初级损伤结构的统计特性，如缺陷数量和位置分布；而机器B使用CNN学习初级损伤结构的形状。对于这两者，训练数据都是通过MD反冲模拟准备的。如果仅使用机器B随机生成候选结构，产生物理上有效的缺陷结构的概率会极低，导致计算效率低下。另一方面，仅使用机器A则需要构建许多统计模型来捕捉缺陷结构的所有相关方面，如缺陷簇大小分布、缺陷间的空间相关性等，这是一个复杂且艰巨的任务。因此，我们设计了一种混合方法，结合了机器A和机器B：机器A根据MD结果提供统计指导的候选结构，而机器B则选择与MD模拟学习到的模式相匹配的结构。我们期望这种方法能够以计算高效的方式生成具有MD级别准确性的真实缺陷结构。
本章的后续部分将详细解释本研究中使用的计算方法：第2.2节介绍MD模拟设置；第2.3节介绍机器A的GPR模型；第2.4节介绍机器B的CNN模型。
**2.2. MD模拟**
**2.2.1. 反冲模拟设置**
为了生成机器A和B构建和测试所需的数据，使用了LAMMPS代码[23]进行MD反冲模拟。对于[100]位移的模拟，使用了144×72×72的超胞（1492992个W原子）；对于其他位移方向的模拟，使用了110×100×90的超胞（990000个W原子）。系统首先在30 K和1巴的压力下使用NPT系综和Nose-Hoover恒温器[24] [25]以及Shinoda等人提出的压力控制方案[26]进行平衡，该方案结合了Martyna等人的静水压力方程[27]和Parrinello与Rahman的应变能[28]。随后，在给原子提供相应的PKA额外动能后，使用NVE系综和变时间步长方法模拟时间演化。具体来说，LAMMPS代码的fix dt/reset命令设置了每个时间步长内原子的最大允许位移为0.025 Å（Xmax），最大允许时间步长为4 fs，每个时间步长内原子的最大允许动能变化为15 eV，以及变量时间步长的更新频率为5。每次反冲模拟一直进行到缺陷数量充分收敛，这通常在6~15 ps内实现，具体取决于反冲能量。例如，32 keV反冲的MD模拟中缺陷数量的时间演化在补充材料（第T节）中提供。在这些设置下，NVE模拟过程中的能量得到了合理的保守。例如，在32 keV反冲的MD模拟中，总能量的典型波动约为1 eV，这一点通过补充材料中的两个测试MD模拟案例得到了验证（第J节）。我们认为这种能量波动水平对于本研究的目的来说是足够小的。
对于势能模型，我们使用了Derlet等人[30]最初参数化、随后Bjorkas等人[31]修改的嵌入原子方法（EAM）[29]势能，该方法能够与Ziegler-Biersack-Littmark（ZBL）势能在短距离内平滑连接，适用于反冲模拟。这种势能模型不仅给出了与实验[32] [33]相当的阈值位移能量，还准确再现了在高能量辐照下钨中的大位错循环和缺陷聚集行为[34]。
在当前的MD反冲模拟中，没有考虑电子效应，即电子阻止（ES）[35]和电子-声子耦合[36]。虽然可以使用摩擦力[37]有效地结合ES，这种方法已在以往的研究中被广泛采用，但电子-声子耦合的处理更为复杂，且尚未完全建立[9]。为了评估电子散射（ES）的潜在影响，我们使用了基于SRIM代码[38]中的阻止能力表的摩擦力进行测试模拟。忽略ES会导致缺陷产生的数量被高估多达20%，这表明对于本研究中考虑的反弹能量（低于32 keV），ES虽然影响有限，但不可忽视。这种差异可以通过使用损伤能量来解决，损伤能量被定义为对缺陷形成有贡献的PKA能量的一部分。因此，机器A和B被训练来仅再现由原子碰撞过程产生的辐射损伤，而不考虑电子效应。我们仅在没有ES效应起显著作用的情况下使用损伤函数。此外，在模拟单元的外围区域对原子应用了恒温器，这种方法在之前的研究中被频繁用来减轻能量和动量的重新进入，但在当前的模拟中没有实施。测试模拟确认，应用外围恒温器对产生的缺陷数量和形态几乎没有影响。关于ES和外围恒温器应用的详细结果以及用于生成分子动力学（MD）数据的LAMMPS输入脚本之一，在补充材料（第A节）中提供。
2.2.2. 反弹能量和方向
由于缺陷结构受到PKA的反弹能量和方向的极大影响[4]，我们需要为不同的反弹能量和方向条件准备辐射损伤数据。一种简单的方法是使用MD模拟创建所有可能的反弹条件来生成训练数据。然而，反弹条件的多样化增加了生成训练数据所需的时间。相反，使用较少的反弹条件可能会降低机器的回归和分类性能。因此，我们研究了能够最佳平衡机器性能和时间效率的适当反弹能量和方向区间。具体来说，我们通过改变训练数据的大小来测试构建的机器的性能，分别称为能量测试和方向测试。
在能量测试中，使用了9个能量：0.5、1、2、4、6、8、16、24和32 keV。参考文献[39]显示，在DEMO聚变反应堆环境下的W的PKA能量谱中，大于10 eV的平均PKA能量为3.2 keV，而该谱仍以低概率延伸到数百keV。大约98%和96%的反弹能量大于1 eV和90 eV的PKA（分别对应于W的阈值位移能量）的反弹能量低于32 keV[39]。因此，我们将本研究的最大PKA能量设定为32 keV，以便在保持MD数据生成成本适中的同时，覆盖大部分相关的PKA群体。尽管将能量范围扩展到数百keV是未来工作的一个课题，但在本研究中，我们采用了一种混合方法，将我们的方法与BCA代码结合来模拟超过32 keV的反弹，如第4.3节所述。在方向测试中，准备了15个方向，使用方法如下：首先，使用之前的研究[40]、[41]中的方法，在单位球面上生成了500个准均匀点。随后，考虑到面心立方（bcc）晶格的对称性，只选择了满足x≥y≥z≥0的点，从而得到11个方向。此外，我们包括了[100]、[110]和[111]方向作为主要方向，以及[321]方向作为bcc金属（包括W）中具有高阈值位移能量的典型方向[7]、[32]、[41]。能量和方向测试的结果在第3.1节中描述。由于有限温度下的碰撞模拟结果受到原子坐标因热振动而波动的显著影响[41]，因此对于每个135种设置（9个能量×15个方向），收集了100个具有不同反弹事件开始时间的样本。这15个方向在补充材料（第D节）中列出。
2.2.3. 分子动力学（MD）结果的处理
MD模拟结果是原子的三维（3D）笛卡尔坐标。由于两种机器的关键信息都是缺陷位置信息，因此使用Wigner-Seitz（WS）缺陷分析[8]从原子位置数据中提取了空位和间隙位错（SIAs）的笛卡尔坐标。然后，将MD模拟单元中的所有空位和SIAs的（x, y, z）坐标（在本研究中为144×72×72或110×100×90的超胞）除以晶格常数，并在移动坐标后转换为125×125×125的单元格，使得PKA的初始位置为（50, 50, 50）。注意，我们没有区分每个单元格的角位点和体心位点。因此，如果给定单元格的两个位点都是空位的，例如，空位的数量被计为2个。接着，如图2中的流程图所示，缺陷位置信息被转换成用于在机器A上进行生成过程回归（GPR）和在机器B上进行卷积神经网络（CNN）训练的形式。
2.4. 机器B使用卷积神经网络（CNN）
对于机器A上的GPR，需要准备和学习损伤结构的统计数据。在本研究中，我们使用了缺陷数量分布和缺陷位置分布。缺陷数量分布是从每个135个反弹能量和方向设置的MD结果中提取的，因此创建了135个概率质量函数（PMFs）。对于缺陷位置分布，提取了每个反弹能量和方向的空位和SIAs的x/y/z坐标，因此总共生成了810个（3个坐标×2种缺陷类型×135个反弹设置）PMFs。这些PMFs被分为低能量（0-8 keV）和高能量（8-32 keV）两组，然后分别用作机器A上GPR的训练和测试数据集。更多细节在补充材料（第B-1节）中提供。
对于机器B上的CNN，为了减小数据的大小和维度，将缺陷的3D坐标投影到三个平面上，得到了六组125×125的2D单元格数据：（x, y）和（y, z）以及（x, z），分别对应空位和SIAs。例如，将（5.5, 30.0, 23.5）的标准化坐标投影为（x, y）=（5, 30）、（y, z）=（30, 23）和（x, z）=（5, 23）。然后，为了将这些信息当作2D图像数据进行处理，計算了每个投影单元格中空位和SIAs的数量，并通过将缺陷数量设置为每个像素的值来构建125×125像素的数据。最后，将所有6个数组（每个方向3个投影）中的数据合并成一个（125, 125, 6）维度的集成数组，宽度 and 长度各为125，深度为6，用于输入到机器B的CNN中。尽管使用完整的3D数据可以保留更多的结构细节，但由于数据的极端稀疏性和3D卷积的大量内存需求，在计算上不切实际。例如，单个30 keV的MD结果会产生一个由超过780万个值组成的125×125×125×4的数组，而其中只有大约200个值对应于实际缺陷。在本研究中，我们采用了2D投影方法，因为它仍然可以保留重要的结构特征，同时显著提高了计算的简单性和效率。如果通过适当的预处理能够提取有意义的数据，使用3D数据将是一个未来的研究课题。
2.3. 机器A使用生成过程回归（GPR）模型
机器A的目标是生成满足MD获得的辐射缺陷重要结构统计特性的候选缺陷结构。这些统计特性可以用任何反弹能量的平滑概率分布函数（PDFs）来表示。然而，由于我们只在135个反弹设置中进行了MD反弹模拟，缺陷数量和位置分布是离散的，因此MD结果的缺陷统计特性是PMFs。为了从MD统计中获得平滑的PDF，GPR使用第2.2.3节中获得的PMFs来推断PDF。输入有四个类别：（1）反弹能量，（2,3）两个对应于球坐标theta和phi的反弹方向参数，以及（4）用于数量分布的缺陷数量和用于位置分布的缺陷位置。输出是相应的概率。使用了径向基函数（RBF）核，初始方差为1.0，初始长度尺度为（1, 1, 1, 1），然后通过最大化边际似然来优化方差和长度尺度。由于每个类别的输入值都是标准化的，我们分别乘以一个缩放系数来为每个输入类别分配权重。如果缩放因子太小，两个相邻变量将被视为几乎相等，或在训练中的影响较小。例如，在低能量情况下，每个样本的缺陷位置变化相当小，因此需要一个相对较大的缩放因子。噪声常数是我们GPR模型中最有影响力的超参数。较小的噪声常数可以提高GP回归的准确性，但也增加了过拟合的风险。因此，我们将噪声常数设置为不会导致过拟合的最小值。因此，对于每个低（0-8 keV）和高（8-32 keV）反弹能量集合，分别获得了1个缺陷数量的PDF和6个缺陷位置的PDF，包括空位和SIAs的x/y/z坐标。GPR模型参数优化的更多细节在补充材料（第F节）中提供。
使用GPR模型的PDF，可以根据给定的反弹能量和方向条件生成候选缺陷结构。首先，生成一个在[0,nmax)范围内的随机整数n和一个在[0,1)范围内的随机分数。如果形成n个缺陷的概率密度大于这个分数，n就被设置为Frenkel对的数量。否则，重复这个过程直到确定Frenkel对的数量。我们设置nmax = 100，这高于在MD模拟中发现的Frenkel对的最大数量。其次，生成一个在[0,mmax)范围内的整数m和一个在[0,1)范围内的随机分数，然后如果在一个坐标m处有缺陷（SIA或空位）的概率密度大于这个分数，m就被设置为x坐标。mmax应等于机器B中使用的单元格大小，因此mmax = 125。同样，y和z坐标也被确定。随后，缺陷位置被设置为（x, y, z）或（x+0.5, y+0.5, z+0.5），分别代表bcc晶格的角位点和体心位点，概率相等。对于每个空位和SIAs重复这个过程n次，生成由n个空位和n个SIAs组成的候选缺陷结构。这个过程保证了生成的候选缺陷结构满足从MD结果生成的GPR模型估计的缺陷数量分布和缺陷位置分布。为了避免同一个bcc晶格位点被多个SIAs或空位重复占据，机器A在生成候选结构时检查位点占据情况并在已经分配给另一个缺陷的位点进行重新采样。
2.4. 机器B使用卷积神经网络（CNN）
2.4.1. 数据准备
机器B的目的是将给定的缺陷配置分类为真数据（True-Data），它有合理的概率是由MD反弹模拟生成的，或者假数据（False-Data），它由MD反弹模拟生成的概率非常低。为此，我们使用了CNN。为了训练CNN，我们准备了真数据和假数据的参考数据。参考真数据是从MD模拟中获得的。数据被格式化为（125, 125, 6）维度的数组数据，如第2.2.3节所解释的。参考假数据通过两种方式生成：简单随机结构和准随机结构。简单随机结构数据是通过随机分布缺陷生成的。首先，Frenkel对的数量被随机确定为4到80之间的一个整数。然后，在125×125×125的超胞中以空位和SIA的均匀分布形式生成随机缺陷结构，并将其转换为(125, 125, 6)维度的数组数据。由于简单随机结构数据缺乏辐射缺陷的特性（如真实的缺陷分布），因此分类变得非常简单，不能单独用于训练分类器。为了提高分类器的性能，我们添加了准随机结构数据，这些数据在某些方面类似于 molecular dynamics (MD) 模拟的结果，但在实际的 MD 模拟中很少生成。为此，我们使用了机器A，它可以准随机生成满足 MD 结果中观察到的缺陷数量和位置分布的候选缺陷结构，如第2.3节所解释的。我们假设机器A生成真实数据的概率足够低。然而，应该注意的是，机器A偶尔可能会生成真实数据，这种情况发生的可能性取决于 recoil 能量。机器A 对机器B训练的影响在第3.3节中讨论。总共准备了1500个简单随机结构和54000个准随机缺陷结构，作为机器B的训练和测试用的虚假数据。对于真实数据，通过MD模拟准备了13500个数据，每种recoil设置对应100个MD数据。在评估分类器性能时，训练和测试集是分别为每个recoil能量类别或recoil方向类别独立准备的，而不是从所有recoil条件的样本中随机选取的。
2.4.2. 卷积神经网络
CNN架构在图3中示意性地展示。有三个卷积层用于提取特征，并使用了ReLU激活函数，三个最大池化层用于下采样，以及一个批量归一化层。输入数据通过这些层，CNN提取特征。特征提取后，卷积层的最终输出被送入全连接层。最后，使用2个ReLU层和1个SoftMax层将输出分类为真实数据或虚假数据。CNN为真实数据类别输出一个softmax概率。在当前的PRDgen实现中，根据这个输出，接受或拒绝机器A生成的候选结构，使用的决策阈值为0.5。
3. 结果
3.1. 机器A的GPR性能
机器A配备了14个GPR模型，用于两种能量类别（高和低）的7个PDF（1个用于Frenkel对的数量，6个用于SIA和空位的3维位置）。每个GPR模型通过调整噪声常数进行优化，其回归性能通过验证均方根误差（RMSE）与真实数据（即MD结果）进行比较来评估。此外，我们检查了缺陷位置和数量的拟合趋势作为PKA能量和PKA方向的函数，以验证没有过拟合，并且得到的拟合结果是平滑和合理的。基于对不同噪声常数的这些评估，对于每个PDF类别，我们选择了在避免过拟合的同时实现最低验证RMSE的模型。
例如，图4a展示了GPR模型为[110]方向和32 keV recoil能量预测的SIA y坐标的PDF，以及相应的训练和测试数据。图4b和图4c分别展示了GPR模型根据recoil能量和recoil方向生成的彩色图。如图4a中SIA y位置分布所示，GPR模型不仅以合理的精度再现了MD结果（图4a），而且还生成了在能量和方向域上平滑插值的PDF（图4b和图4c）。关于构建的GPR模型性能的更多细节在补充材料（第B-2节）中提供。
3.2. 机器B的测试条件优化
3.2.1. 能量测试
为了确定必要的和足够的recoil能量和方向类别数量，我们测试了机器B的插值能力。在能量测试中使用了9种能量：0.5、1、2、4、6、8、16、24和32 keV。我们从仅使用一种recoil能量开始训练，并增加recoil能量类别的数量来评估训练数据数量和类型对分类性能的影响。
表1显示了用于recoil能量测试的数据列表。对于每个能量类别，使用1200个MD模拟数据（总共1500个MD数据，每个方向150个数据）作为真实数据的训练数据，以及660个准随机数据和140个简单随机数据作为虚假数据的训练数据。测试数据由1 keV到40 keV之间以1 keV间隔计算的1000个MD模拟数据组成（每种能量25个MD数据；为了简化，仅在[100]方向附近）和相当数量的（990个）准随机结构作为虚假数据。
3.3. 机器B的分类性能
基于第3.2节中关于必要和足够数据的测试结果，我们决定在训练数据集中包括(recoil能量类别为(1, 2, 4, 8, 16, 32) keV和15个选定方向的recoil方向。尽管15个方向的数据数量超过了第3.2.2节中确定的必要数据数量，但由于15个方向的案例在图7中的性能略优于3个方向的案例，我们使用了所有15个方向的数据以最佳利用本研究中积累的MD数据。超参数被优化以最大化整体分类精度，详细信息在补充材料（第C节）中提供。随后，使用5400个测试数据评估了创建的机器B的分类性能，获得了96.5%的准确率。测试数据涵盖了9个能量类别（训练中使用了6个类别：0.5、6和24 keV），15个方向类别，2种数据类型（真实数据或虚假数据），以及每种情况下20个样本。这相当于5400个数据点（9×15×2×20）。图8显示了每个recoil能量类别的错误率，即每个recoil能量的分类失败率。
根据这些测试结果，我们认为测试9方案是一个最佳选择，该方案包括从1 keV开始以加倍步骤递增的recoil能量数据，例如1 keV、2 keV、4 keV、8 keV、16 keV和32 keV。它在1-40 keV范围内的所有能量上都表现良好，整体准确率超过97%。
3.2.2. 方向测试
如第2.2.2节所述，准备了15个方向。类似于能量测试，我们首先只使用一个recoil方向进行测试，逐渐增加更多的recoil方向类别，并分析了性能变化。表2展示了4个测试案例中使用的recoil方向。对于每个方向类别，分别为真实数据和虚假数据准备了450个数据。每个测试案例的训练数据集大小与recoil方向的数量成正比，保持真实数据与虚假数据的比例为1。测试结果显示在图7中。只有bcc晶格中满足x≥y≥z≥0的不可约方向被展示。在包含所有15个方向的测试4中，每个测试方向的分类错误率低于6%。测试3仅包括三个低指数方向([100], [110], [111])，在[621]方向附近的精度损失较小。相比之下，仅使用一个或两个方向无法获得高精度，如测试1和2所示。因此，我们认为使用[100], [110], [111]方向的数据是一个最佳选择。它实现了所有方向平均超过95%的精度，并且没有方向的错误率显著较大。
表2. 方向测试的详细信息。
方向测试指标 训练数据集（recoil方向）
[111] 900
[110] [111] 1800
[100] [110] [111] 2700
4 13500
3.3. 机器B的分类性能
基于第3.2节关于必要和足够数据的测试结果，我们决定在训练数据集中包括(recoil能量类别为(1, 2, 4, 8, 16, 32) keV和15个选定方向的recoil方向。虽然15个方向的数据数量超过了第3.2.2节中确定的必要数据数量，但由于15个方向的案例在图7中的性能略优于3个方向的案例，我们使用了所有15个方向的数据以充分利用本研究中的MD数据。超参数被优化以最大化整体分类精度，详细信息在补充材料（第C节）中提供。随后，使用5400个测试数据评估了创建的机器B的分类性能，获得了96.5%的准确率。测试数据涵盖了9个能量类别（训练中使用了6个类别：0.5, 6和24 keV），15个方向类别，2种数据类型（真实数据或虚假数据），以及每种情况下20个样本。这总共是5400个数据点（9×15×2×20）。图8显示了每个recoil能量类别的错误率。为了评估机器B在recoil能量类别方面的性能变化，准备了另一个测试数据集，固定recoil方向为[100]，recoil能量间隔为1 keV。这个额外的数据集总共包含3200个样本，包括每个能量类别从1 keV到32 keV的25个真实数据样本和75个虚假数据样本。尽管在高recoil能量下的错误率足够低，但在6 keV以下，虚假数据的错误率增加到超过15%。特别是在8 keV以下的能量范围内，错误率相对较高。这一现象的原因及实际影响在第4.1节中讨论。
为了进一步测试机器B的性能，我们使用了训练中使用的最终能量类别集（1, 2, 4, 8, 16, 32 keV）进行了留一法评估。结果显示，即使在移除一个能量类别后，分类器仍然具有广泛的鲁棒性，尽管在低能量范围内的精确度对于移除某个低能量类别相对敏感。详细信息在补充材料（第P.3.4节）中提供。最后，通过将机器A与机器B连接起来，建立了一个基于机器学习的辐射缺陷生成器PRDgen，用于生成候选缺陷结构，并评估生成的结构是否合理，如图1所示。当PRDgen接收到反冲能量和方向作为输入时，机器A首先根据使用分子动力学（MD）结果通过GPR学到的概率来确定弗伦克尔对（Frenkel pairs）的数量。如果弗伦克尔对的数量为零，则直接返回“无缺陷”的结果，而无需进一步处理。否则，机器A根据训练好的缺陷位置分布生成候选缺陷结构，然后机器B验证该缺陷结构是否合理。这个过程重复进行，直到机器B将机器A生成的缺陷结构分类为真实数据。机器A、机器B和PRDgen的代码都在补充材料（第G节）中提供，同时为了方便用户，还提供了PRDgen的Windows版本。Windows版本无需安装额外的库即可运行。
我们评估了PRDgen与分子动力学（MD）模拟在以下方面的性能：（1）生成辐射缺陷结构所需的计算时间；（2）生成缺陷结构中包含的平均弗伦克尔对数量；（3）生成缺陷结构的外观；（4）从初始接触点（PKA）到生成点缺陷的距离统计；（5）两种缺陷之间的相关性，即最近邻壳层分布和簇尺寸分布；（6）生成的缺陷结构的多样性。评估结果在以下小节中描述。
3.4.1 计算速度
我们测试了三种反冲能量（10、20和30 keV），这些能量没有包含在训练数据集中，反冲方向是随机选择的。每个计算都在同一台计算机上执行，该计算机是一部具有20个CPU核心的Linux服务器（2 × Intel Xeon Silver 4210；2.2 GHz；每CPU 10个核心）。使用LAMMPS的MD模拟是通过Open MPI在20个核心上并行执行的，而PRDgen则在单个CPU核心上运行，因为当前版本的PRDgen没有实现并行化。对于LAMMPS和PRDgen，都从计算开始到生成辐射损伤结构的墙钟时间进行了评估。对于每种反冲能量情况，准备了100个样本以获得平均值和标准差。计算时间估计方法的更多详细信息在补充材料（第H节）中提供。构建PRDgen的一次性成本，包括训练数据集的生成和两台机器的训练，在表3中没有考虑。这一成本在补充材料（第O节）中单独总结。
表3. MD与PRDgen的计算速度比较。括号中的值表示样本标准差。PRDgen在单个CPU核心上运行，而使用LAMMPS代码的MD则在20个CPU核心上并行运行。
| 反冲能量 | 计算类型 | 平均时间（秒） | 平均弗伦克尔对数量 |
|--------|---------|-----------|-------------------|
| 10 keV | MD | 660 | 10.5 |
| | PRDgen | 2.31 | 3.73 |
| 20 keV | MD | 1535 | 20.3 |
| | PRDgen | 8.61 | 15.30 |
| 30 keV | MD | 3499 | 32.2 |
| | PRDgen | 31.97 | 31.1 |
表3表明，PRDgen生成的弗伦克尔对数量与MD一致，同时计算速度在30 keV时加快了75倍，在10 keV时加快了200倍。由于单核PRDgen作业可以同时运行，因此生成许多损伤结构的有效吞吐量可以大致与可用CPU核心的数量成比例增加，假设竞争（I/O和内存带宽）可以忽略不计：例如，在20个CPU核心时，计算速度可以加快1500倍。
PRDgen的计算标准差远大于表3中的MD计算标准差。实际上，PRDgen的运行时间并不遵循简单的高斯分布，由于其依赖于重复的提议采样和基于卷积神经网络（CNN）的过滤，因此呈现出长尾分布。因此，尽管平均值提供了生成许多辐射损伤结构成本的有用估计，但它并不能完全捕捉运行时间统计。为了补充表3中的基于平均值的比较，补充材料（第I节）提供了（i）提议级别的接受率（接受的类别预测数除以CNN评估的提议总数）和（ii）每种反冲能量和方向的端到端墙钟百分比（P50/P90/P99）。
虽然这不是直接的比较，因为速度是使用不同的计算机资源评估的，但DDPM会比PRDgen更高效。具体来说，使用单个GPU（Nvidia RTX 4090）的DDPM大约比使用112个CPU核心（8 × Intel Xeon E5-2680v4；2.4 GHz；每CPU 14个核心）的MD快105倍。
3.4.2 视觉评估
图9显示了使用BCA（Iradina代码[42]）、PRDgen和MD（LAMMPS代码）生成的30 keV反冲能量下的主要损伤结构。在BCA计算中，有效（阈值）位移能量被设置为90 eV，这是ASTM标准[43]推荐的，并通过之前的MD研究[32]验证。在图9中，为了确保三种方法之间的一致性，不仅在BCA模拟中考虑了ES，还在MD和PRDgen计算中考虑了ES。在MD中，通过能量依赖的摩擦力来实现ES，其低能量截止值为10 eV，这是使用SRIM代码实现的ES功率表计算的。在PRDgen中，ES效应通过NRT损伤函数来考虑，如补充材料（第A节）中所述。
下载：下载高分辨率图像（93KB）
下载：下载全尺寸图像
图9. 由(a) BCA、(b) PRDgen和(c) MD在30 keV反冲能量下生成的缺陷结构。模拟结果在一个30nm×30nm×30nm的单元格中进行了比较。红色、蓝色和绿色点分别表示硅间隙（SIAs）、空位（vacancies）和PKA起点。如正文中所述，所有三种计算都考虑了ES效应。
在图9中，BCA模拟的损伤数量显然高估了与MD相比的缺陷数量，而PRDgen提供的结果与MD相当。PRDgen生成的损伤结构大小也与MD生成的类似。然而，高反冲能量下形成的缺陷簇并没有充分再现，PRDgen产生的缺陷在空间上倾向于均匀分布。这些特性在接下来的章节中进行了定量评估。补充材料（第E节）提供了MD和PRDgen生成的主要损伤结构的额外视觉比较。
3.4.3 缺陷数量统计
我们分析了PRDgen与MD生成的缺陷数量的统计。首先，提取了每个样本的弗伦克尔对总数，并计算了样本平均值。这里的样本平均值定义为第2.2.2节描述的15种反冲方向和8种反冲能量（1/2/4/6/8/16/24/32 keV）每组的100个样本的平均弗伦克尔对数量。此外，还计算了每个能量下所有方向的样本平均值的平均值和标准差。这个标准差表明了每个反冲能量的缺陷数量对反冲方向的依赖性。图10比较了PRDgen统计和MD统计中8种反冲能量的样本平均值的平均值和标准差。
下载：下载高分辨率图像（103KB）
下载：下载全尺寸图像
图10. MD与PRDgen在生成的弗伦克尔对数量方面的比较，作为反冲能量的函数。误差条表示每个反冲能量类别的样本平均值的标准差，表明了生成的弗伦克尔对数量对反冲方向的依赖性。
为了评估超出基于平均值分析的一致性，我们使用两个互补的距离比较了从MD和PRDgen结果中提取的弗伦克尔对计数的完整分布：双样本Kolmogorov–Smirnov统计量（KS）和Wasserstein-1距离（W1）。对于每种反冲条件，我们收集了n = 100个独立的MD级联和n = 100个被PRDgen接受的输出，并直接从经验缺陷计数样本NFP计算KS和W1。KS定义为经验累积分布函数之间的最大绝对偏差，而W1是以弗伦克尔对为单位的地球移动距离。如补充材料（第K节）中所述，这些分析确认PRDgen即使在完整分布级别上也高保真地再现了MD的缺陷计数统计。
3.4.4 缺陷位置分布
通过测量点缺陷从初始PKA位置的均值和最大距离来量化点缺陷的全球范围，针对4种反冲方向（[100]、[110]、[111]和[321]）。PRDgen的结果在图11中与MD结果进行了比较。PRDgen生成的均值和最大距离与MD模拟得到的结果非常接近。正如MD结果所示，已知SIAs位于损伤中心相对较远的位置。这一特性也被PRDgen准确再现了。尽管这种基于距离的分析没有捕捉到缺陷结构的所有复杂结构特征，但它至少表明了缺陷从初始PKA位置的范围被准确表示。
下载：下载高分辨率图像（186KB）
下载：下载全尺寸图像
图11. 比较了MD和PRDgen在32 keV PKA能量下，从初始PKA位置到生成的点缺陷位置的统计分布：(a) 平均距离和(b) 最大距离。误差条表示每个数据点的平均值的标准误差。
3.4.5 缺陷之间的相关性
我们检查了最近邻壳层分布（NND）和簇尺寸分布作为多个缺陷之间的两个关键相关性。对于SIA-SIA和空位-空位对的NND分析，首先考虑周期性边界条件得到每个缺陷与其他缺陷的最近邻距离。随后，定义了壳层k的NND分数为fk=1/N∑i=1NIk(diNN)，其中N是选定类型的缺陷数量，diNN是缺陷i的最近邻距离，Ik(d)是一个指示函数，定义为Ik(d)=1，如果dbelongstothek−thnearestneighborshell，否则为0。
图12显示了1 keV和32 keV下[100]反冲方向的NND分数分析结果。绘制的值和误差条分别代表从每种条件下生成的100个损伤结构计算得到的fk的均值和标准差。对于SIA-SIA对，MD结果显示，在较高能量下，最近邻壳层（k = 1, 2, 3）的分数随着反冲能量的增加而增加，表明SIAs在更高能量下更加聚集。PRDgen并没有清楚地再现这一趋势，这表明当前模型没有充分捕捉到MD中观察到的短程SIA-SIA相关性的能量依赖性增加。相比之下，对于VAC-VAC对，无论是MD还是PRDgen都没有显示出随着反冲能量增加而显著增加的低壳层分数。这表明短程空位-空位接近性不如SIA明显，因此MD和PRDgen之间的差异对于空位来说较小。
下载：下载高分辨率图像（207KB）
下载：下载全尺寸图像
图12. 比较了1 keV和32 keV下SIAs和空位的最近邻壳层分数分布，从左到右分别显示了1 keV和32 keV下的SIA-SIA NND分数和VAC-VAC NND分数。符号和误差条分别表示100个结构的均值和标准差。
图13显示了同一反冲能量和方向条件下SIAs和空位的簇尺寸分布。缺陷根据固定的bcc壳层连接标准被分配到同一簇中，其中对于SIAs，两个缺陷如果在3个最近邻壳层内相连，则被认为是相连的；对于空位，则是2个最近邻壳层。这些连接标准是基于之前在bcc-W [6]、[44]、[45]中的簇识别研究选择的。特别是，早期与密度泛函理论计算一致的结果表明，最稳定的di-间隙配置由两个平行的<111>哑铃组成，它们由第三最近邻距离分隔，支持了对SIAs的3个最近邻标准。团簇大小1对应于一个孤立点缺陷。下载：下载高清图片（214KB）下载：下载完整尺寸图片图13. 在1 keV和32 keV下，SiAs和空位的团簇大小分布。从左到右分别是1 keV和32 keV下的SiAs团簇大小分布，接着是1 keV和32 keV下的空位团簇大小分布。符号和误差条分别表示100个结构的平均值和标准差。在图13中，对于SiAs，分子动力学（MD）结果显示，小团簇大小的概率随反冲能量的增加而不会显著增加，而大团簇大小的尾部逐渐出现。相比之下，PRDgen在所有能量下都主要表现为团簇大小1，大于1的团簇大小出现的概率很低。这证实了当前的PRDgen框架无法充分再现低反冲能量下已经存在的小SiAs团簇的概率，也无法再现高反冲能量下较大SiAs团簇的出现。对于空位，MD和PRDgen都主要集中于团簇大小1，表明空位聚集相对较弱。对于NND和团簇大小分布分析，[110]、[111]和[321]反冲方向以及8 keV反冲能量的相应结果在补充材料（第M节）中提供。这些结果与本节中描述的[100]反冲方向在1 keV和32 keV下的趋势相似。3.4.6. 生成的损伤结构的多样性使用基于CNN的损伤结构分类器，即Machine B，分析了PRDgen生成的损伤结构与MD生成的损伤结构的多样性。为此，我们对CNN从损伤结构中提取的特征应用了主成分分析（PCA），并在由第一和第二主成分（PC1和PC2）跨越的二维空间中比较了它们的多样性。图14比较了MD、Machine A和PRDgen生成的损伤结构的CNN特征图。Machine A生成的数据对应于准随机结构。Machine A和PRDgen的输出之间的区别在于生成的结构是否通过Machine B分类器。这一比较清楚地表明，虽然MD和PRDgen在低能量下生成的缺陷多样性相当吻合，但随着反冲能量的增加，这种一致性显著下降。这种行为可能与PRDgen无法再现缺陷聚集有关，而在高能量下聚集变得更加明显。下载：下载高清图片（193KB）下载：下载完整尺寸图片图14. 由MD、Machine A和PRDgen生成的在1、8和32 keV反冲能量下的辐射损伤结构的CNN特征图。每种方法绘制了400个结构，每个反冲方向（[100]、[110]、[111]和[321]）各100个结构。在1 keV下，准随机结构与MD结果的重叠相对较大，表明在某些情况下，低能量下生成的缺陷位置可能不会表现出强烈的缺陷间相关性。同时，即使在低能量（例如1 keV）下，MD和PRDgen在特征空间中的结构分布也存在可检测的差异，这表明在其他情况下缺陷间相关性已经起作用。尽管Machine B滤波器在整个能量范围内提高了生成缺陷结构的多样性和保真度，但这种改进程度仍不足以与MD结果达到良好对应，特别是在高反冲能量下。4. 讨论4.1. 低反冲能量下分类器性能的局限性如图6和图8所示，在低能量范围内出现了相对较高的错误。为了进一步分析其原因，我们调整了训练数据集中真实数据与虚假数据的比例，并评估了分类性能，如图15所示。下载：下载高清图片（107KB）下载：下载完整尺寸图片图15. 不同训练数据比例下的分类性能比较。训练数据比例为真数据:虚假数据为1:1.5（黑色）、2:1（红色）、4:1（黑色）和8:1（绿色）。左侧和右侧面板分别显示了真数据错误率和虚假数据错误率。在图15中观察到8 keV以下的性能有两个趋势。首先，错误率随着能量的降低而迅速增加。其次，分类性能受到用于训练的真数据与虚假数据比例的显著影响。具体来说，增加虚假数据的比例往往会减少虚假数据的错误（较低的假阳性率），但会增加真数据的错误（较高的假阴性率），而增加真数据的比例则产生相反的效果。与此一致的是，正确识别的真数据比例（即召回率）随着真数据比例的增加而提高。这一观察表明，用于生成准随机缺陷结构的Machine A在低能量下以相对较高的概率生成真数据缺陷结构，因为低能量下的缺陷空间分布较窄，缺陷没有大量聚集。图14中显示的CNN特征图也支持这一解释：Machine A生成的缺陷结构与MD生成的缺陷结构在所跨越的特征空间中的重叠在低能量下更大。在这种情况下，如果对真数据赋予较高权重（通过在训练中包含更多真数据来实现），Machine A生成的类真数据结构倾向于被分类为真数据，以最小化总错误。另一方面，如果对虚假数据赋予较高权重，Machine A可以高概率生成的一些真数据将被分类为虚假数据。因此，前者的虚假数据错误率增加，后者的真数据错误率增加。简而言之，在低能量下，真数据和虚假数据的分类准确性之间存在权衡。应该注意的是，即使在高能量下，Machine A偶尔也会生成真数据；然而，这种情况的发生频率较低，通常约为1%，这是根据Machine A生成被Machine B接受的损伤结构的频率估算的。因此，我们得出结论，低反冲能量下的相对较高错误率并不表示分类器的性能低下，而主要是反映了缺陷结构的固有特性。在这种情况下，需要将Machine A生成的类真数据结构判定为真数据，以便所有合理的缺陷结构都被正确分类为真数据。这意味着我们需要在一定程度上牺牲虚假数据的分类准确性。因此，我们调整了用于训练的真数据与虚假数据的比例，以便将真数据错误率降低到大约10%以下，同时尽可能抑制虚假数据错误的增加。4.2. 高反冲能量下PRDgen性能的局限性通过高能量碰撞级联MD模拟得到的辐射损伤结构可以分为三类：（1）缺陷分布广泛且聚集较少的结构，（2）缺陷密集且显著聚集的结构，以及（3）反冲原子在其初始动能下传播较远距离的结构，导致缺陷分布延长。图16展示了每种类型典型的MD生成缺陷结构。PRDgen几乎总是生成类型（1）的结构，而无法生成类型（2）和（3）的结构。下载：下载高清图片（84KB）下载：下载完整尺寸图片图16. MD模拟得到的典型缺陷结构。（a）缺陷分布广泛且聚集较少。（b）大部分缺陷显著聚集。（c）缺陷沿反冲原子的延长轨迹分布。MD模拟是在32 keV反冲下进行的，每个缺陷结构都被放在30nm×30nm×30nm的单元中进行了比较。红色和蓝色点分别表示SIA和空位。为了澄清为什么不生成类型（2）和（3）的结构是由于生成器（Machine A）还是分类器（Machine B）的问题，我们分别为10、20和30 keV反冲额外生成了40个MD数据，并用它们来测试Machine B是否能够将其分类为真数据。分类性能结果如图17所示，同时给出了Machine A生成的准随机缺陷结构和PRDgen生成的真数据结构的结果。在30 keV反冲的40个MD数据中，包含了所有三种类型的结构：类型（1）占75%，类型（2）占20%，类型（3）占约5%。如图17所示，Machine B正确分类了大多数MD数据，除了两个10 keV反冲的数据，这表明不仅类型（1）的缺陷结构，类型（2）和（3）的缺陷结构也可以被识别为可能的辐射缺陷结构。还应注意，在图17中，PRDgen获得的缺陷结构的分类得分大多在50-80%之间，略高于虚假数据和真数据之间的界限50%，远低于MD数据对于大多数30 keV反冲数据的90%以上的分类得分。这一结果表明，Machine B正确识别了高能量反冲下缺陷聚集的重要性。下载：下载高清图片（191KB）下载：下载完整尺寸图片图17. Machine B计算的分类得分的标准化直方图。上图显示了Machine A生成的准随机结构的分类得分，下图显示了MD模拟和PRDgen得到的结构的分类得分。对于MD和PRDgen，每种反冲能量使用了40个数据。对于准随机缺陷结构，每种反冲能量使用了300个数据。因此，我们得出结论，PRDgen无法生成类型（2）和（3）的缺陷结构，因为Machine A生成这些类型的概率非常低，这是因为Machine A没有考虑增强聚集或延长的统计模型。此外，x、y和z位置分布的概率密度函数（PDFs）是独立使用的，也没有考虑不同缺陷坐标之间的相关性，这也使得Machine A难以形成类型（2）和（3）的缺陷结构。总之，PRDgen可以合理地生成类型（1）的缺陷结构，而类型（2）和（3）的缺陷结构几乎无法形成，这是由于Machine A的性能有限。在低反冲能量下，类型（1）占主导地位，因此PRDgen的表现尚可。然而，在高反冲能量下，由于类型（2）和（3）的缺陷结构在现实中（即MD中）以不可忽视的频率形成，但Machine A难以生成，所以PRDgen的性能下降。同时，Machine B可以正确分类所有类型的缺陷结构，包括类型（2）和（3）。4.3. 与BCA集成以预测更高反冲能量下的损伤由于PRDgen是在32 keV以下的MD模拟结果上训练的，它只能可靠地生成这个能量范围内或略高范围内的缺陷结构。为了将预测扩展到更高反冲能量，我们将PRDgen与BCA代码Iradina [42]结合使用。在这个框架中，Iradina将高能量反冲分解为能量低于32 keV的次级 knock-on 原子（SKAs）。每个能量低于32 keV的SKA被传递给PRDgen以生成缺陷结构，而能量较高的SKAs则在Iradina中进一步模拟，直到它们的能量降至32 keV以下。然后，所有SKA生成的缺陷结构被组合起来形成总体损伤配置。Iradina在原子碰撞主要由两个核之间的短程排斥相互作用主导的高反冲能量下提供了足够的准确性，这是BCA有效的范围。然而，在较低能量下，多体碰撞和多体原子间相互作用变得显著，降低了BCA的准确性。此外，BCA倾向于高估缺陷的产生，通常产生与NRT模式 [46] 一致的结果，但远高于MD模拟 [6]、[17]、[47] 的结果。相比之下，PRDgen通过基于MD数据派生的统计分布生成缺陷来避免这个问题，尽管它仍然局限于最多32 keV的能量范围。我们结合了这两种方法的方法，称为BCA+PRDgen，利用了两种方法的优势来描述广泛的能量范围内的初级损伤结构。我们对75 keV反冲进行了BCA+PRDgen的测试。图18将其损伤结构与MD和完整BCA的结果进行了比较。对于每种方法，分别为PKA方向[111]和[321]生成了30个缺陷结构。所有情况下都考虑了ES效应。在BCA+PRDgen中，通过在BCA阶段应用ES，然后使用NRT损伤函数对SKAs进行事后校正来模拟这一效应。对于[321]反冲，形成的Frenkel对数量分别为：BCA为401.4（15.4），BCA+PRDgen为66.7（15.1），MD为78.1（8.4）。括号中的值表示标准偏差。对于[111]反冲，生成的Frenkel对数量分别为402.7（19.0），68.1（14.6）和67.9（16.4）。这些结果表明，BCA+PRDgen大大减少了BCA对缺陷的过度估计，产生了接近MD的结果。为了评估点缺陷的空间分布，图19比较了从初始PKA位置出发的SIA和空位的归一化直方图，确认BCA+PRDgen在缺陷空间分布的再现上与MD非常吻合。为了进一步检验BCA+PRDgen框架的性能，我们对50 keV和60 keV的反冲进行了额外的测试，结果再次证实了其与MD的良好一致性。详细结果在补充材料（第N节）中有描述。
下载：下载高分辨率图像（95KB）
下载：下载全尺寸图像
图18. 由75 keV的321反冲引起的缺陷结构，在（a）BCA、（b）BCA + PRDgen和（c）MD模拟中的表现。所有模拟结果都是在30nm×30nm×30nm的单元中进行的比较。红色、蓝色和绿色点分别代表SIA、空位和PKA的起点。如正文中所述，所有三种计算中都考虑了ES效应。
下载：下载高分辨率图像（161KB）
下载：下载全尺寸图像
图19. 由75 keV PKA引起的初级辐射损伤的点缺陷距离的归一化直方图，该数据来自MD模拟和BCA+PRDgen。x轴值以W晶格参数（Åa0=3.165Å）进行归一化。
关于计算效率，使用第3.4.1节中相同的计算资源，BCA+PRDgen平均需要80秒，而MD则需要21,788秒，速度提高了270倍。由于BCA+PRDgen不使用模拟单元，因此其运行时间随PKA能量的增加而增加的程度较小，这与MD不同。这些结果表明，BCA+PRDgen能够在广泛的能量范围内高效且有效地生成初级损伤结构。
尽管有这些优势，未来仍有几个问题需要研究。在缺陷形态方面，MD结果显示随着反冲能量的增加，尤其是SIA的缺陷会聚集[6]，而BCA+PRDgen未能再现大的缺陷簇，如图18所示。这一限制源于两个因素：（1）机器A生成聚集缺陷的概率过低，如4.2节关于高反冲能量下PRDgen性能的限制所讨论的；（2）PRDgen仅针对32 keV以下的反冲能量进行了训练，因此无法再现更高反冲能量下的聚集行为。关于因素（2），先前的研究[34]、[47]表明，在bcc-W的MD模拟中，缺陷簇的数量和大小在大约40 keV以上迅速增加，这是当前模型未能学到的现象。因此，为了成功地将PRDgen与BCA代码结合，确定两种方法之间的适当切换能量至关重要。较低的切换能量可以降低总计算成本，包括生成机器所需的成本，但可能会因为BCA在低能量下的可靠性降低而影响准确性。例如，在研究硅中的沟道效应时，BCA已经在大约100 eV以上的反冲能量下与MD进行了验证[48]。然而，据我们所知，目前还没有系统性地评估BCA在模拟与SKAs形成相关的碰撞动力学方面的有效性。需要进一步的研究来确定一个平衡准确性和效率的最佳切换能量。
4.4. 与DDPM的比较
本研究中开发的两种基于ML的辐射缺陷生成器PRDgen和Liao等人[22]开发的DDPM，在架构、学习策略和缺陷结构生成方法上存在本质差异。首先，DDPM使用生成式AI模型，而PRDgen结合了GPR和CNN模型。DDPM在一个统一模型中整合了训练和生成，并直接将PKA能量张量插入到噪声样本中。相比之下，PRDgen采取了模块化的方法。缺陷结构对物理条件（如PKA能量和方向）的依赖性是使用GPR（即机器A）明确且统计性地建模的。此外，独立的基于CNN的分类器机器B被训练用来独立于其PKA条件评估候选结构的结构真实性，从而可以间接考虑机器A未能隐式捕获的统计特征。通过机器A直接且明确地考虑缺陷结构的统计特性是PRDgen的一个重要特点，这使得它能够持续改进和扩展。
然而，由于当前版本中只实现了两个统计模型（一个用于缺陷数量，另一个用于缺陷位置），PRDgen在三个方面（1）速度、（2）缺陷聚类统计的准确性以及（3）考虑的反冲能量范围（PRDgen最多32 keV，DDPM最多100 keV）上未能达到DDPM的性能。相比之下，PRDgen中只考虑了反冲方向依赖性，而DDPM中没有考虑。第（3）点主要是由于为本研究准备的训练数据有限。对于更本质的问题（即（1）和（2），需要进一步的研究来提高PRDgen的性能和可用性。
实际上，计算速度和缺陷聚类统计的准确性这两个问题密切相关。如第4.2节所讨论的，机器A无法再现MD模拟在高反冲能量下常见的缺陷聚集现象，因为它没有包含促进聚集形成的统计模型。因此，如图17所示，即使机器A生成的结构被机器B分类为真实数据，其评分仍然相对较低。此外，如表3所示，随着反冲能量的增加，PRDgen的效率也会下降，因为机器A在更高能量下生成真实数据的能力大大降低。这些发现表明，在机器A中添加统计模型以增强缺陷聚集效果将使其生成物理上合理的缺陷结构的速度更快，从而提高PRDgen的整体速度。开发这样的模型，包括那些用于聚类形成概率、聚类大小分布以及x、y、z方向上缺陷坐标之间的相关性的模型，是未来工作的一个重要方向。
4.5. 缺陷间相关性的初步评估
迄今为止的结果和讨论一致表明，PRDgen的一个主要限制是缺乏缺陷间相关性。尽管识别和实现关键的缺陷间相关性需要更全面的研究，但我们进行了一个初步测试，将一些简单的相关性纳入机器A中，以检验引入缺陷间相关性是否如预期那样能够提高PRDgen的保真度和效率。在这个初步测试中，我们考虑了两种缺陷间相关性：（1）聚类大小分布和（2）最近邻距离（NND）。机器A被扩展以包括这两个统计量，除了GPR模型实现的缺陷数量和位置统计量。针对SIA和空位，分别构建了这两个相关性的拟合方程作为反冲能量的函数。为了简化起见，没有考虑反冲方向依赖性。通过这种扩展，缺陷聚集得到了增强，再现了MD模拟得到的聚类大小分布，而整体缺陷数量和缺陷位置的整体趋势仍然遵循原始机器A中实现的概率分布。此后，将这种包含两种缺陷间相关性的扩展版本的PRDgen称为PRDgenEX。用于此初步测试的聚类大小分布和最近邻壳层分数的统计模型的准备和实现是在补充材料（第U节）中描述的。
图20展示了PRDgenEX的性能。图20仅展示了一组原始PRDgen表现相对有限的代表性指标。包括NND分数和不同反冲能量下的聚类大小分布在内的更多详细结果、端到端墙时比较、生成缺陷结构的视觉比较、详细的分类分数分布以及它们的分布都在补充材料（第R节）中提供。
下载：下载高分辨率图像（443KB）
下载：下载全尺寸图像
图20. 基线PRDgen与扩展模型PRDgenEX在（a）SIA–SIA NND分数、（b）SIA聚类大小分布、（c,d）1 keV和32 keV损伤结构的PCA特征图、（e）分类分数以及（f）计算时间方面的比较。
对于最近邻和聚类大小分布，图20以SIA为例展示了结果。与PRDgen相比，与MD结果的一致性显著提高。聚类大小分布与MD结果几乎完全一致，而NND分数仍然存在一些限制，特别是第二和第三邻居分数的低估。这是因为当一个缺陷同时具有第一个和第三个邻居原子时，根据最近邻分布的定义，只有第一个邻居被计算在内，从而降低了第三邻居分数。虽然纠正这个问题很简单，但在这个初步测试中我们没有尝试纠正它，因为这个测试的目的是为了概念验证。
PCA特征图也显示出在1 keV和32 keV下再现MD结果的显著改进。同时，PRDgenEX仍然无法完全再现MD结果所覆盖的特征空间。这表明，额外的缺陷间相关性仍然很重要，可能包括空位-SIA位置相关性。
正如预期的那样，改进不仅体现在生成损伤结构的保真度上，还体现在生成效率上。由于机器A没有包含促进聚集形成的统计模型，因此无法再现MD模拟中在高反冲能量下常见的缺陷聚集。结果如图17所示，即使机器A生成的结构被机器B分类为真实数据，其分数也相对较低。此外，如表3所示，随着反冲能量的增加，PRDgen的效率也会下降，因为机器A在更高能量下生成真实数据的能力大大降低。这些发现表明，在机器A中添加统计模型以增强缺陷聚集效果将使其生成物理上合理的缺陷结构的速度更快，从而提高PRDgen的整体速度。开发这些模型，包括那些用于聚类形成概率、聚类大小分布以及x、y、z方向上缺陷坐标之间的相关性的模型，是未来工作的一个重要方向。
4.5. 缺陷间相关性的初步评估
到目前为止的结果和讨论一致表明，PRDgen的一个主要限制是缺乏缺陷间相关性。尽管识别和实现关键的缺陷间相关性需要更全面的研究，但我们进行了一个初步测试，将一些简单的相关性纳入机器A中，以检验引入缺陷间相关性是否能够如预期那样提高PRDgen的保真度和效率。在这个初步测试中，我们考虑了两种缺陷间相关性：（1）聚类大小分布和（2）最近邻距离（NND）。机器A被扩展以包括这两个统计量，除了GPR模型实现的缺陷数量和位置统计量。分别为SIA和空位构建了这两个相关性的拟合方程。为了简单起见，没有考虑反冲方向依赖性。通过这种扩展，缺陷聚集得到了增强，再现了MD模拟得到的聚类大小分布，而整体缺陷数量和缺陷位置的全局趋势仍然遵循原始机器A中实现的概率分布。此后，这种包含两种缺陷间相关性的扩展版本的PRDgen被称为PRDgenEX。
图20说明了PRDgenEX的性能。图20仅展示了一组原始PRDgen表现相对有限的代表性指标。补充材料（第R节）提供了包括NND分数和不同反冲能量下的聚类大小分布在内的更多详细结果、端到端墙时比较、生成缺陷结构的视觉比较、详细的分类分数分布以及PCA特征图及其分布。
下载：下载高分辨率图像（443KB）
下载：下载全尺寸图像
图20. 基线PRDgen与扩展模型PRDgenEX在（a）SIA–SIA NND分数、（b）SIA聚类大小分布、（c,d）1 keV和32 keV损伤结构的PCA特征图、（e）分类分数以及（f）计算时间方面的比较。
对于最近邻和聚类大小分布，图20以SIA为例展示了结果。与PRDgen相比，与MD结果的一致性显著提高。聚类大小分布与MD结果几乎完全一致，而NND分数仍有一些限制，特别是第二和第三邻居分数的低估。这是因为当一个缺陷同时具有第一和第三邻居原子时，根据最近邻分布的定义，只有第一邻居被计算在内，从而降低了第三邻居分数。尽管纠正这个问题很简单，但在这个初步测试中我们没有尝试纠正它，因为这个测试的目的是为了概念验证。
PCA特征图也显示出在1 keV和32 keV下再现MD结果的显著改进。同时，PRDgenEX仍然无法完全再现MD结果覆盖的特征空间。这表明，额外的缺陷间相关性仍然很重要，可能包括空位-SIA位置相关性。
正如预期的那样，改进不仅体现在生成损伤结构的保真度上，还体现在生成效率上。由于扩展后的机器A的softmax分数分布显示出显著的上升，因此接受率显著提高。例如，在32 keV时，墙时减少了大约80倍。结合第3.4.1节中提到的原始PRDgen的加速效果，PRDgenEX的整体速度大约是MD模拟的6000倍。考虑到用于估计速度的计算架构差异（MD使用32个CPU核心，DDPM使用GPU，PRDgenEX使用单个CPU核心），PRDgenEX的有效生成速度很可能接近DDPM。
总之，尽管需要进一步的研究来建立适当的缺陷间相关函数并将它们系统地纳入PRDgen，但这个初步测试表明PRDgen框架是一种有前途的方法，用于生成辐射损伤结构，并且有可能实现与当前最先进的辐射损伤生成器DDPM相当的效率。
4.6. 机器B中True/False边界的敏感性
另一种提高当前PRDgen保真度的方法是使机器B中的True/False边界更加严格。在当前实现中，我们将机器B中softmax输出的True/False边界（即决策阈值）设置为0.5，相当于50%的分类分数。这个值可以区分大多数准随机损伤结构和MD损伤结构，因为前者数据接近0，后者数据接近1，如图17所示。然而，由于MD损伤结构的值接近1，我们可以将这个边界设置为一个更高的值，比如0.9，以在机器B和MD损伤结构之间获得更好的softmax分数一致性。同时，由于计算成本在很大程度上受到机器B接受率的影响，使用更高的决策阈值会增加计算成本。
图21显示了在保持机器A设置不变的情况下，仅改变机器B中CNN过滤器的决策阈值到softmax值0.3、0.5、0.7和0.9时获得的结果。正如预期的那样，计算成本对阈值非常敏感。另一方面，原始PRDgen的一个局限性在于无法再现缺陷间的相关性，即使提高决策阈值也几乎无法改善这种情况。这表明，仅仅通过改进CNN分类器是无法再现局部结构特征的，必须结合第4.5节中介绍的缺陷生成逻辑的修改。关于机器B中使用的CNN过滤器的接受规则以及决策阈值对SIA–VAC最近邻分布影响的更多细节，请参见补充材料（第Q节）。
图21. 机器B中CNN决策阈值对（a）端到端计算时间、（b）SIA-SIA NND比例以及（c）32 keV下[100]反冲方向的SIA簇大小分布的影响。
4.7. Wigner-Seitz分析的局限性及数据处理替代方法
本研究中的所有缺陷结构均指的是通过WS缺陷分析识别出的空位和SIA的位置。WS缺陷分析根据完美晶体格子构建Voronoi单元，计算每个单元中的原子数量（占据度），并将占据度≥2的单元视为SIA，占据度=0的单元视为空位。尽管WS方法在辐射损伤分析中被广泛采用，但它有两个明显的局限性。首先，占据度=2没有问题，但在高PKA能量或高温下，可能会出现占据度≥3的单元，从而导致缺陷计数精度的不准确。在目前用于训练和评估的13500个MD数据集中，只有57个案例（约0.4%）的Voronoi单元的占据度为3或更高。为了更详细地探讨这个问题，我们在图22中总结了按PKA能量划分的案例发生率。这些结果表明，尽管在整个数据集中这类情况仍然很少见，但它们的频率随着反冲能量的增加而增加。因此，在本研究中，由于使用的反冲能量相对较低，WS方法的这一局限性并未造成问题。然而，如果我们将模型扩展到非常高的PKA能量，例如超过100 keV，WS方法的缺陷计数错误可能会产生显著影响。其次，WS方法仅指示了点缺陷占据的晶格位置，并未保留缺陷配置的原子几何结构。因此，在预处理过程中会丢失SIA哑铃形或群集状结构的局部结构信息，这些特征也因此被排除在学习阶段之外。
图22. 在所有用于训练/评估的MD数据集中，每个PKA能量下包含至少一个占据度≥3的Voronoi单元的结构比例。
对于未来的发展，采用比WS预处理更能保留局部结构信息的缺陷分析方法可能有助于解决PRDgen目前的局限性。例如，在参考文献[49]中提出的方法中，每个原子都通过其局部原子环境的旋转不变描述符向量来表示，并通过与参考描述符进行比较（使用基于距离和概率的标准）来识别与缺陷相关的环境。这种表示方法可以保留比WS分析中使用的标量占据度更丰富的信息，从而在ML训练开始前更好地保留SIA哑铃形或群集状结构的局部几何信息。在参考文献[50]中提出的方法中，通过从SIA哑铃形/群集状结构构建图并分类连接组件和循环，可以进一步分析缺陷形态，从而能够解决方向、环状结构和混合形态的问题。这些方法可能为未来的PRDgen模型提供更好的基础，特别是对于高能量级联和聚集缺陷。
总体而言，WS分析并不是识别和表示缺陷结构的最佳方法，尤其是对于高能量级联或几何形状复杂的缺陷簇，尽管在当前研究的能量范围内它仍然是一个实用的折中方案。未来对PRDgen的改进将从采用更丰富的局部环境表示和/或结合形态感知的分析方法来增强基于WS的缺陷识别中受益。
4.8. 机器B中改进的投影策略
在当前实现中，机器B使用了WS处理后的缺陷结构的三个正交2D投影。为了检验投影策略本身是否影响分类器的性能，我们额外测试了增加了三个45度投影（分别对应<110>、<011>和<101>的扩展六视图表示。同样的扩展也应用于MD、准随机和简单随机数据集，并重新训练和评估了机器B。然而，真数据和假数据的误差率趋势与使用原始三个视图表示得到的结果大致相似，差异并不显著或系统化，不足以支持明显的改进。这一结果表明，在基于图像识别的当前框架内，仅仅增加投影视图的数量并不能显著恢复在2D投影过程中丢失的结构信息。有关此额外测试和结果的详细信息，请参见补充材料（第L节）。尽管如此，这一结果可能反映了由于缺乏缺陷间相关性所带来的局限性过于显著，以至于额外的投影所带来的相对较小的改进不太明显。因此，一旦在未来的工作中更全面地纳入缺陷间相关性，重新审视减少与2D投影相关的信息损失的好处将值得考虑。
4.9. 应用于其他材料及提高训练效率的策略
尽管PRDgen框架是通用的，原则上可以扩展到其他材料系统，但这种扩展不太可能通过直接重用当前的bcc-W实现来实现，而是需要对每个目标材料进行重新训练。例如，之前的MD模拟表明，在fcc金属的碰撞级联过程中可以直接形成堆垛缺陷四面体[51]，这表明学习材料特定特征的重要性。由于当前PRDgen的实现需要大量的训练数据（如补充材料[O节]中所述），因此对每种材料进行重新训练是一个实际挑战。
有几个可能的方向可以提高训练过程的效率。首先，对于机器A中实现的统计模型，它们对反冲能量和方向的依赖性在不同材料之间可能存在共同特征。在分析了相对较少数量的代表性材料后，可以识别出辐射损伤结构统计中的这些一般特征，从而减少为其他材料构建统计模型所需的努力。此外，一旦明确了控制这些统计模型的基本物理机制，通过适当的核设计在GPR中进行一定程度的外推可能是可能的，这可以减少机器A所需的训练数据量。对于机器B，补充材料[S节]中已经观察到其具有外推能力。
其次，基于bcc-W训练的CNN可以通过几种为基于CNN的图像分类开发的迁移学习方案扩展到其他材料，以高效地为具有相似特征的目标构建分类器。例如，可以重用预训练的CNN作为特征提取器，同时仅重新训练目标材料的最终分类层，或者使用新材料的有限数据集部分微调较高层[52],[53]。此外，采用共享特征提取层和任务特定输出头的多任务学习框架也可能有助于跨相关材料学习通用表示[54]，明确学习任务间特征共享程度的架构（例如交叉编织网络）可能进一步促进这种迁移[55]。
第三，也可以利用之前的MD模拟结果。在这方面，国际原子能机构（IAEA）开发的CascadesDB[56]数据库，其中包含了碰撞级联的MD模拟结果，可能为未来扩展当前框架提供有用的支持。探索有效利用这些现有数据库的方法，可能还需要通过额外的模拟来弥补现有数据集的局限性，这将是一个有趣的研究课题。
5. 结论
在这项研究中，我们开发了一个基于ML的框架，用于生成具有给定反冲能量和方向的bcc-W中的初级辐射损伤结构，并将其实现为一个名为PRDgen的代码。PRDgen由两台机器组成：一个候选缺陷结构生成器和一个辐射损伤结构分类器，分别称为机器A和机器B。机器A使用基于MD反冲模拟结果的GPR生成满足合理缺陷数量和位置分布的准随机缺陷结构。机器B使用基于相同MD结果的损伤结构进行CNN训练。通过系统地测试训练的反冲能量和方向采样条件，我们确定了当前框架的有效数据集组成。特别是，发现从1到32 keV的加倍步骤的反冲能量采样是合适的，而方向测试显示了[100]、[110]和[111]方向的重要性。基于这些结果，最终实现采用了六个反冲能量类别和15个选定的方向来构建当前版本的PRDgen。
PRDgen能够生成合理再现MD损伤结构的几个重要全局特征的缺陷结构，同时将计算成本降低了大约两个数量级。特别是，PRDgen能够以合理的准确性再现缺陷数量和缺陷位置分布的统计特性。此外，混合BCA+PRDgen方案将方法的适用范围扩展到了当前训练范围之外的反冲能量，并且计算成本显著低于MD。同时，当前研究揭示了当前PRDgen实现的一个局限性，即在更高反冲能量下无法再现缺陷的最近邻和簇大小分布。这一局限性主要源于机器A，因为当前版本中实现的统计模型没有包括缺陷间的相关性。
为了验证当前PRDgen实现这一局限性的解释，我们进行了一个初步测试，将缺陷邻接和聚类统计纳入机器A。扩展后的模型命名为PRDgenEX，显著提高了准确性和效率。这些结果表明，建立和纳入考虑这些局部特征的额外统计模型是未来发展的重要方向。此外，扩展能量范围并能够预测更精细的缺陷结构（如SIA哑铃形方向）将是改进框架的重要步骤。
尽管当前存在局限性并且需要进一步开发，但PRDgen框架具有重要的优势，这激励了对其进行持续改进。尽管机器A不依赖于第一性原理物理模型，但它明确使用了MD观察到的缺陷数量和缺陷位置分布，并将其纳入缺陷结构的生成中。在初步测试中检验缺陷间相关性的影响时，从MD模拟构建的统计模型显著提高了机器A的性能。这样，虽然由于深度学习的性质，机器B的可解释性较低，但通过引入额外的具有物理意义的统计模型，可以系统地改进机器A。因此，确定机器A内再现辐射损伤结构所需的统计成分可能会揭示控制初级损伤结构的关键特征。在之前的研究中，只检查了有限的结构统计信息，如形成的缺陷数量和缺陷簇大小分布。因此，PRDgen框架不仅提供了一个实用的替代模型来生成辐射损伤结构，还为揭示初级辐射损伤背后的统计特性提供了一个系统的方法。
CRediT作者贡献声明
Takuji Oda：写作 – 审稿与编辑、监督、软件、项目管理、方法论、调查、形式分析、数据管理。
JongHyeon Park：写作 – 原始草稿、可视化、验证、软件、方法论、调查、形式分析、数据管理。
