《Small》：Large-Scale Cooperative Sulfur Vacancy Dynamics in Two-Dimensional Mos2 From Machine Learning Interatomic Potentials
编辑推荐：
这篇研究采用机器学习原子间势函数（MLIP），实现了对二维二硫化钼（MoS2）单层中硫（S）空位在纳秒尺度上动态迁移的高精度分子动力学（MD）模拟。研究揭示了空位辅助的协同跳跃机制，阐明了空位聚集成簇并最终形成长达数十纳米线缺陷的原子级过程，为理解其在催化、传感及新型忆阻器（memristor）等应用中的功能行为提供了关键的动力学见解。
1. 引言
硫空位迁移在二维MoS2的多种现象中扮演核心角色，涉及催化、化学传感、光致发光、电子输运以及下一代计算器件（如基于二维MoS2的忆阻晶体管和忆阻元件）等领域。理解并最终控制此类器件功能，需要对支配空位迁移和聚集的动态机制有定量的认识。这些机制强烈依赖于局部原子环境（如缺陷结构、电荷态和迁移路径的变化），并表现出协同特性。例如，室温下足够的迁移率需要相邻空位间的相互作用（即空位辅助传输），这可将迁移能垒从2-3 eV降低至0.6-0.8 eV。传统的分子动力学（MD）方法受限于时间尺度，难以捕捉相对缓慢的硫空位排列动态演化。虽然基于密度泛函理论（DFT）的能量势垒进行动力学蒙特卡洛（kMC）模拟可以扩展可及的时间和长度尺度，但它依赖于预定义的迁移事件，对于超越简单缺陷构型的复杂情况，构建精确模型极具挑战性。
本研究通过使用专门训练的机器学习原子间势函数（MLIPs）进行MoS2的MD模拟，来应对这一挑战。MLIPs是一种新兴方法，可在纳秒时间尺度上以接近DFT的精度建模大规模空位迁移和缺陷演化。具体而言，我们比较了两种MLIP框架：采用高斯近似势（GAP）的即时学习方法，以及对等变基础模型（MACE MP-0）进行微调的方法。通过训练得到的MLIPs，我们模拟了硫空位在多个纳秒时间尺度上的迁移，捕捉了传统原子模拟技术无法触及的结构演化过程，揭示了MoS2单层中硫空位动力学的底层机制，并为实验中观察到的由小缺陷团簇聚集成横跨数十纳米的扩展线缺陷现象提供了原子级别的解释。
2. 结果与讨论
2.1 MLIP模型的训练
我们应用了两种机器学习架构来构建MLIPs：基于原子簇扩展描述结构环境的等变图神经网络MACE，以及使用原子描述符并实现即时学习的GAP。对于每种架构，我们测试了两种不同的训练数据生成策略。所有模型都针对DFT计算的力和能量进行微调，但在训练集的结构数量以及这些结构的获取方式上有所不同。模型在测试集上的力和能量误差大多在预期范围内。我们强调，这些测试误差是基于从使用待测MLIP的MD模拟中采样的平衡和非平衡构型进行严格评估得到的，这通常会导致比主要关注近平衡结构的MLIP基准测试报告的误差更高。
2.2 对协同空位辅助跳跃MLIP预测势能曲线的评估
为了更有针对性地评估MLIPs的准确性，我们使用微动弹性带（NEB）方法计算并比较了特定硫缺陷跳跃的势能曲线。我们重点研究了一个硫原子跳入一对相邻空位旁间隙位点的过程（即空位辅助跳跃）。将不同MLIPs版本的NEB势能曲线与DFT参考计算进行比较发现，MACE基础模型甚至给出了定性错误的结果。所有经过微调的MACE模型表现均显著优于基础模型，但其与DFT NEB计算的一致性很大程度上取决于训练集的选择。使用从基础模型的MD轨迹中采样的结构进行微调的模型，在势能曲线最小值附近准确，但在过渡区域显著偏离DFT参考值。相比之下，使用从从头算分子动力学（AIMD）模拟中采样的大量快照（如50,000个）进行微调的MACE模型，能够非常准确地预测势能曲线，其预测的激活能与DFT参考值相符。
对于GAP MLIPs，在包含双空位或三空位团簇的结构上通过即时学习训练的模型，在能量最小值附近与DFT势能曲线吻合良好，在最大值附近仅有微小偏差，预测的前向跳跃能垒（0.74和0.78 eV）接近DFT值（0.83 eV）。而在单一空位结构上训练的模型，或在多种结构不同的空位团簇上训练的模型，则显示出较大偏差。总结来说，使用从AIMD模拟中采样的大型训练集，我们能够微调出一个非常准确地描述这些硫跳跃的模型。而基础模型及所有基于基础模型MD结构微调的模型都严重低估了过渡态的势能。
2.3 其他跳跃过程的讨论
先前的研究已对硫原子跳入空位的所有可能构型的能量势垒进行了全面分析，并将我们评估的跳跃识别为具有最低能量势垒的跳跃之一。相对较低的势垒归因于目标位点旁存在一个“辅助”空位，这使得硫原子能够进入间隙位点。相比之下，硫原子跳入单一空位的能量势垒要高得多（约1.81 eV）。根据指数阿伦尼乌斯定律，在任何相关技术温度下，硫原子跳入两个相邻空位之一的概率远高于跳入单一空位。此外，我们还计算了硫原子跳入由三个空位组成的直线中的势能曲线，发现其能量势垒与跳入两个相邻空位相似（约0.8 eV）。
2.4 MLIP预测的径向分布函数评估
硫-硫径向分布函数（RDF）提供了另一种评估MLIPs在MD模拟中质量的方法。结果显示，作为参考的AIMD模拟（分别使用CP2K和VASP）与最准确的GAP和MACE MLIPs预测的RDF之间具有良好的一致性，只有MACE基础模型无法准确预测RDF。
2.5 硫空位动力学
获得能够以令人满意的精度预测势能曲线和RDF的MLIPs后，我们将其应用于确定不同缺陷结构在单层MoS2中的均方位移（MSD）和缺陷跳跃速率，时间尺度远超AIMD能力范围。扩展的可及模拟时间也允许识别协同效应。对包含三个硫空位的结构进行的MD模拟显示，我们最好的经过微调的MACE和GAP模型预测的MSD彼此合理吻合，平均每个移动硫原子在1 ns时间跨度内移动约3 ?。而MACE基础模型则预测了更大的位移和几乎线性的MSD增长，表明硫原子在该时间尺度上不受限制地扩散，这归因于其系统性的“软化”，导致能垒被错误降低。
当结构中存在两个或三个空位时，可观察到硫原子的显著位移，MSD最终达到约5 ?2的平台，表明硫原子运动在这个有限系统中在纳秒时间尺度上受到空间约束。具有两个空位的系统MSD收敛更快，而具有三个空位的系统在趋于平稳前MSD持续增加数纳秒。钼原子在我们的MD模拟中表现出非常有限的动力学，即使存在硫空位时也是如此。
通过在800至1200 K温度下进行长时间的MD模拟，我们分析了相对罕见的硫原子跳跃的统计规律。跳跃速率与温度的关系符合阿伦尼乌斯定律，使我们能够从斜率估计硫跳跃的“有效”能量势垒。对于两个空位的结构，从MD模拟得出的“有效”势垒（0.81 eV）与NEB计算的“静态”值（0.83 eV）非常吻合。而对于三个空位的结构，由MD模拟阿伦尼乌斯拟合得出的能量势垒（0.63 eV）与NEB计算值（0.85 eV）相差0.22 eV。这种差异凸显了动力学效应的重要性，这些效应在静态NEB计算中未被捕捉，例如邻近原子的振动贡献和涉及同时跳跃的协同机制。
通过分析10 ns长的MACE MD模拟中不同初始空位构型附近的位点依赖跳跃概率，我们发现硫原子迁移性主要由“空位辅助跳跃”机制支配。该机制要求存在至少两个相邻的硫空位，且跳跃的硫原子必须是这两个空位的第一近邻，同时跳跃路径上的间隙位点不能被下方原子层中的钼原子阻挡。我们将这些条件称为“跳跃规则”。在小的空位团簇情况下，这些跳跃规则将硫原子迁移性限制在围绕空位团簇的局部区域，这些区域是等边三角形。根据跳跃规则并得到MD模拟证实，所有硫原子跳跃都发生在一个涉及3个或6个硫晶格位点的三角形图案内。这解释了两种情况下MSD的差异。值得注意的是，在MLIP可达的较长时间尺度上，空位三角形内的所有跳跃都以相当的频率发生。
观察到的跳跃规则与动力学约束模型非常相似。分析表明，对于任何有限的空位浓度和无限大的系统，几乎肯定存在一个可以跨越整个系统增长的可动区域。换句话说，一旦可动区域达到临界尺寸，进一步的生长几乎是必然的。例如，在实验可达到的空位密度下（如1%），一个长度为5的空位线，其可动区域有15个相邻硫晶格位点，至少有一个位点是空位的概率高达86%。
2.6 硫空位聚集的讨论
“跳跃规则”概念解释了具有可动硫原子和硫空位的区域的形成，但不足以解释实验事实，例如使用透射电子显微镜（TEM）频繁观察到的扩展直线空位线。MLIP MD模拟则同时考虑了能量差异和构型依赖的跃迁势垒高度差异。
我们观察到，空位团簇可以通过捕获位于其可动区域外但与之相邻的额外空位来实现生长，这导致团簇可动区域的扩大。这种空位捕获机制在之前的经典反应力场MD研究中未被观察到，很可能是由于其不准确的能量势垒预测所致。而本研究开发的MACE MLIP模型能可靠地重现DFT结果，从而揭示了导致实验观察到扩展硫空位团簇的集体空位动力学机制。
空位团簇倾向于沿着其可动区域的边缘排列成直线，以最小化其能量。通过DFT计算比较直线和锯齿形构型的总能量差（可解释为额外空位与已有空位线的结合能）发现，空位结合能随着空位线长度的增加而单调递减。因此，形成能是一条凹曲线，更长的空位线在能量上比总空位数相同的多个较短空位线更有利。同时，锯齿形构型的能量差始终大于其直线对应构型。这两个观察结果共同解释了空位团簇倾向于组织成沿其可动区域边缘的长直线的现象。这种扩展空位线的形成已在实验中通过TEM观察到，并被认为在忆阻元件的功能行为中起关键作用。
此外，在有限温度的MD模拟中，我们通常观察到空位团簇在其可动区域边缘附近处于动态平衡，在转变回直线之前，会继续形成弯曲构型。在更长时间尺度上（对于小的N为数纳秒），空位线甚至可以通过一系列单独的跳跃迁移到其可动区域的另一条边。在某些MD运行中，我们观察到连续跳跃在很短的时间跨度内发生，一个硫原子的迁移刚好在或与第二个跳跃同时发生，使得两个跳跃以相关的方式进行。这凸显了基于MLIP的长时间尺度MD模拟的价值，因为其他方法（如kMC模型）需要事先了解这种相关的空位动力学才能准确捕捉协同空位跳跃。
关于直线缺陷在有限温度下是否会分裂以增加熵的问题，在我们的MLIP MD模拟中（N ≤ 5），在线缺陷初始形成后没有观察到其分裂。偶尔有单个硫原子跳入空位线中心附近的空位，暂时将其分成两段，但它们会在进一步跳跃稳定新构型并导致永久性线重排之前，返回到其原始位点。静态DFT计算证实，存在相当强的能量驱动力抵消熵驱动的空位线瓦解，这解释并证实了我们在纳秒尺度MD模拟中观察到的稳定空位线。
不同空位构型之间形成能和能量势垒的不可忽视的变异，凸显了本研究开发的MLIPs相对于先前用于研究MoS2中硫动力学的kMC模型的一个关键优势。
3. 结论
MoS2中的大型硫空位团簇即使在高温下也会在AIMD计算无法触及的时间尺度上形成并发生动态变化。底层过程依赖于团簇构型的细节，其复杂程度使kMC模拟难以处理。这使得基于MLIP的MD模拟成为唯一可行的替代方案。我们发现，MLIPs只需极少的微调即可准确描述MoS2稳定构型附近的势能面，但为了对过渡态进行定性和定量正确的描述，需要精心构建和选择适当的训练数据集。
我们比较了GAP和等变GNN MACE两种MLIP架构。GAP框架固有地支持即时学习策略，无需“手动”选择代表性训练结构即可获得准确结果。相比之下，MACE的准确性对训练数据集的大小和代表性更为敏感。我们在从AIMD模拟中采样的大型数据集上成功训练了MACE MLIP，该模型对重要的“空位辅助跳跃”的过渡态能量和RDF都做出了高度准确的预测。最终，两种方法都可靠地再现了支配空位迁移的底层势能面。这表明，在确定未来哪种MLIP方法将被最广泛采用时，可用性（理解为应用的简单性和数据需求的效率）可能超过微小的精度差异。
重要的是，训练准确的MLIP模型使我们能够进行远超以往在从头算精度下可及时间尺度的MD模拟，从而获得了对MoS2中硫空位局域动力学的新见解。这些动力学由特定的低能硫跳跃机制（“空位辅助跳跃”）驱动，该机制发生在至少两个空位相邻处。通过分析纳秒级MLIP-MD模拟的MSD、跳跃速率和位点依赖跳跃概率，我们在原子水平上观察到了硫空位集体迁移性的新机制。这些原子过程是静态DFT计算或kMC模型无法触及的，凸显了长时间尺度、高精度模拟的附加价值。
在低空位密度下，只形成小团簇。此类团簇的迁移性被限制在一个三角形区域内，由于支配空位迁移的特定硫跳跃，其边长大约小于或等于团簇内的空位数。通过阿伦尼乌斯关系计算“有效”能量势垒，我们发现这些团簇内的迁移性具有强烈的尺寸依赖性。对不同温度的MD模拟分析也强调了动力学效应的重要性——例如邻近原子的振动贡献和涉及同时跳跃的协同机制——这些是先前用于研究MoS2中硫动力学的静态NEB计算和kMC模拟组合所无法捕捉的。
此外，我们证明小团簇可以捕获位于其受限运动区域相邻的空位，从而扩大其运动范围。原则上，这种机制可以在实验相关的空位密度下导致大甚至无限的可动空位团簇，但需要很长时间才能实现。
MD模拟和静态DFT计算表明，在这些可动区域内，空位优先形成直线以最小化其能量。综上所述，这些发现解释了在室温下观察到大量中等长度空位线，以及在高温下观察到数量较少但显著更长（跨越数十纳米）的空位线的实验现象。此处观察到的小空位团簇的迁移性表明，空位可能优先在晶界处聚集，那里更高的缺陷密度将允许更长程的扩散。这凸显了对多晶MoS2系统中空位动力学进行进一步理论建模的必要性。
我们期望未来的研究能更详细地探索浓度依赖性、历史依赖的玻璃态行为以及大尺度忆阻性之间相互交织的方面。这些都源于任何空位辅助传播所固有的非线性。为此，很可能基于从MLIP MD推导出的浓度和历史依赖的迁移率，开发出有效的连续介质模型。
