运用深度势方法连接Mg-La合金的原子尺度能量学与宏观腐蚀动力学

《Mechanics of Materials》：Bridging atomic-scale energetics and macroscopic corrosion kinetics in Mg-La alloys using a deep potential approach

　　张浩祥|陈佳辉|高若桐|张向阳|何志杰|荣菊|余晓华中国昆明科技大学材料科学与工程学院，昆明650093摘要镁镧合金是一种具有潜力的可生物降解金属材料，但其原子尺度上的结构-性能-腐蚀关系目前仍不够清晰。本研究开发了一种基于人工智能的科学深度势模型，并对其进行了严格验证，该模型能够实现对镁镧合金结构、力学及电化学性能的多尺度精确模拟。该深度势模型达到了密度泛函理论级别的精度，能量预测的均方根误差为5.79×10^-3电子伏/原子，力预测的均方根误差为1.63×10^-3电子伏/埃，且具备很强的通用性。它能够准确再现能量-体积关系、晶格常数、声子色散、弹性性能以及缺陷、表面和内聚能等参数，相对误差均在1%以下。该模型预测的体模量和杨氏模量与第一性原理计算及实验数据高度吻合，证明了材料的力学和动力学稳定性。通过将深度势模型得到的表面能与巴特勒-沃尔默框架相结合，可直接在原子尺度上预测阳极极化曲线。计算得到的腐蚀电位和电流密度（约10^-4.9安·厘米^-2）与实验结果趋势一致，且显示出明显的晶体学取向依赖性。这项工作建立了一个基于物理原理、具有通用性的框架，将复杂合金中的原子能学与宏观腐蚀行为联系起来。1.引言镁基合金因其极低的密度以及优异的力学性能组合，属于最轻的结构金属材料类别，在追求轻量化与高性能设计的现代工业中受到了广泛关注。特别是在生物医学领域，镁基合金凭借其出色的生物相容性以及与人体骨组织极为相近的力学性能，相比传统金属生物材料具有不可替代的价值。在各种镁合金体系中，稀土强化的镁镧合金由于镧元素的微合金化效应而展现出独特优势。这类合金既能实现晶粒细化强化，又能提升耐腐蚀性，因此被视为突破传统镁合金性能局限的有希望的候选材料。然而，镁镧合金的实际应用仍面临两大挑战：首先，难以实现对强度、延展性和韧性等关键力学性能的精准协同调控；其次，在复杂的生理环境或动态腐蚀环境中，仍缺乏对降解速率和腐蚀形貌的有效稳定控制。这些问题严重制约了镁镧合金在工程中的应用以及临床转化。因此，为解决这些难题，有必要探索新的调控策略，深入理解镁合金的材料特性，优化其微观结构设计及加工工艺，从而在降低成本的同时大幅拓展其应用潜力。为系统建立镁合金微观结构与宏观性能之间的内在关联，进而实现其合理的性能设计，当前研究主要依靠计算建模与实验表征这两种互补的方法。这两种方法相互印证，共同构成了该领域的主流研究范式。在计算方面，分子动力学模拟作为原子尺度的核心理论工具，能够动态展示材料在热机械耦合作用下的结构演化、缺陷行为及相变过程，从而为理解性能形成的机制提供动态视角。例如，李等人利用这种方法研究了稀土元素在镁基体中的扩散与偏聚行为，为揭示稀土元素提升镁合金整体性能的微观机制提供了重要理论见解。不过，这类研究大多依赖于传统的经验势函数，如嵌入原子法。这类势函数在描述含稀土体系如镁镧合金中的多体相互作用、电荷转移以及局部化学环境时精度有限，导致预测的能量景观或动力学行为与实际物理情况存在偏差。在实验方面，包括原位电化学测试、表面形貌分析以及成分分析在内的高精度表征技术被广泛用于实时捕捉镁合金在服役条件下的腐蚀动态及界面演变过程。这些实验为理论模型和模拟结果提供了重要验证。然而，实验观测也面临诸多挑战：一方面，复杂腐蚀界面的高时空异质性往往导致获得的机制信息较为局部且零散；另一方面，来自不同技术、不同尺度的数据难以整合与关联。因此，当前研究面临着双重限制：计算模拟受限于经验势函数的精度，而实验观测则受到零散、多尺度数据的制约。所以，克服这两大障碍，整合多尺度信息以建立普适性的设计原则，仍是该领域的核心挑战。近年来，机器学习原子间势函数作为一种“人工智能用于科学”理念下的关键技术，为突破传统经验势函数的精度限制提供了新途径。尤其是基于深度学习的势函数模型，由于其在描述多组分和强非线性原子相互作用方面的卓越能力，被视为实现镁合金高精度原子尺度建模的关键解决方案。与传统依赖预定义函数形式和有限拟合参数的经验势函数不同，机器学习原子间势函数利用深度神经网络自适应学习高维势能面。结合主动学习策略进行数据增强与优化后，它们能够在大幅降低计算成本的同时实现接近第一性原理的精度，从而在精度与可扩展性之间取得良好平衡。大量研究表明，机器学习原子间势函数在捕捉局部化学环境的变化、非线性多体相互作用以及复杂缺陷结构方面具有显著优势。它们能够准确识别传统势函数无法描述的原子间相互作用的细微差异，为预测材料的力学行为、扩散过程及腐蚀动力学等性能提供更可靠的原子尺度依据。例如，孙等人对镁基材料的腐蚀行为进行了系统的第一性原理研究，揭示了成分调控与局部电子结构演变之间的内在关系。他们进一步指出，将高精度的第一性原理计算与分子动力学模拟相结合，是开发性能可预测、行为可机制解释的材料的重要趋势。在力学行为研究方面，程等人研究了塑性变形与缺陷演变，证明基于深度势函数的分子动力学模拟能够准确再现位错运动、孪晶形成及应变局域化等关键动态过程，从而为阐明镁合金的复杂变形机制提供了强大的原子尺度工具。同时，赵等人通过将电化学实验与深度势函数分子动力学模拟相结合，系统研究了镁合金腐蚀过程中的界面反应与结构演变。他们的研究结果表明，机器学习原子间势函数能够与实验观测结果高度一致，大大提升了在真实服役条件下行建模的可靠性与适用性。这些进展共同表明，尤其是深度势函数方法，其在提高镁基合金原子尺度建模精度、阐明复杂服役行为背后的机制以及实现以性能为导向的合理设计方面具有巨大潜力，为镁合金在生物医学工程等关键领域的进一步发展奠定了坚实的理论和方法基础。基于上述研究背景与理论需求，本研究建立了一个高精度、具有物理可解释性的多尺度计算框架，用于系统研究镁镧合金的力学响应及电化学腐蚀行为。具体而言，该研究采用了基于深度神经网络的深度势方法，并将其与经典的巴特勒-沃尔默电化学动力学模型相结合，实现了对原子尺度能量学与宏观阳极极化行为的统一描述。通过引入主动学习策略对构型空间进行自适应采样，所构建的深度势模型在能量预测和原子力预测方面的均方根误差分别仅为5.79×10^-3电子伏/原子和1.63×10^-3电子伏/埃，对于未见过的新构型的决定系数R2可达0.98。在此基础上，该深度势模型能够准确再现镁镧合金的能量-体积关系及平衡晶格参数，预测的晶格常数偏差极小。此外，关键弹性常数的预测误差约为0.2吉帕，为描述体相变形及界面原子相互作用提供了定量可靠的力学与结构参考。进一步地，从深度势函数计算得到的表面能被作为具有物理意义的原子尺度参数引入巴特勒-沃尔默模型，建立了原子尺度表面稳定性与电化学阳极溶解动力学之间的直接定量关联。利用这一多尺度耦合框架，可以直接预测镁镧合金的阳极极化曲线，计算得到的腐蚀电流密度与先前报道的实验结果吻合良好。这些发现表明，将深度势方法与电化学动力学建模相结合，为基于物理原理对复杂镁基合金的腐蚀行为进行建模与定量预测提供了一种实用且有效的方法。2.计算方法DP-Gen算法借助DeepMD-Kit工具包构建了一个包含数据集获取、神经网络势函数训练以及模型逐步优化的自洽流程。这种迭代方式有助于高效构建高精度的深度势模型。此外，该框架还能分析从头算分子动力学中出现的复杂能量变化，因此能够生成能够可靠模拟材料在不同热力学状态下行为的原子间势函数，展现出较强的通用性。最后，它适用于多种多组分合金体系。虽然目前的势函数是为镁镧体系定制的，但背后的深度势/DP-GEN框架具有很高的通用性。若要将其扩展到其他镁-稀土二元或镁镧-X三元体系，只需在DFT训练数据集中添加目标元素特有的构型（如体相、缺陷、表面以及多组分局部环境），然后使用相同的主动学习流程和模型架构对特定体系的势函数进行重新训练，即可为更广泛的镁基合金开发高精度的机器学习原子间势函数，实现数据驱动的精准建模。2.1.DFT计算所有第一性原理计算均采用维也纳从头算模拟软件包VASP完成。交换-相关效应通过广义梯度近似下的佩尔杜-伯克-恩策霍夫函数处理，而电子-离子相互作用则采用投影增强波方法描述，从而实现了对电子结构的精确建模。为确保布里渊区内的均匀采样，通过系统的收敛性测试确定了相关参数，包括500电子伏的平面波能量截断值，以及间距为0.15埃^-1的以Γ点为中心的k点网格。计算精度还通过双重收敛标准进行控制：连续电子自洽循环之间的总能量变化小于10^-6电子伏，最大原子力分量低于0.02电子伏/埃。随后，从候选体系的DFT计算中获得了大约5000个原子构型构成的训练数据集，用于训练深度势模型。此外，还对能量-体积关系及相关平面缺陷能（如层错能）进行了系统的定量分析。同时，还利用Phonopy软件结合VASP计算了声子色散谱。为确保深度势模型针对镁镧体系具有高精度和强通用性，最初构建的训练数据集涵盖了多种局部原子环境，共计5000个构型（见表S1），包括平衡体相结构、低配位表面环境、点缺陷结构以及受应变扰动的非平衡构型。如此全面的覆盖使得模型能够系统学习影响镁镧体系结构稳定性、力学响应及表面相关性质的局部原子环境。所有构型的能量、原子力和应力张量均通过第一性原理计算得出，随后被用于训练深度势函数。值得注意的是，明确纳入表面和缺陷结构有助于提升模型准确捕捉低配位及局部畸变环境的能力，这对于后续的表面能分析及腐蚀相关性质的研究至关重要。2.2.MD模拟结构优化及分子动力学模拟则是使用LAMMPS软件包进行的。基于体心立方结构的镁镧合金，构建了一个5×5×5的超胞，其中包含总共250个原子。在三个空间方向上都采用了周期性边界条件，以消除表面效应。分子动力学模拟的时间步长为1飞秒，每个时间步长内原子位置、速度和加速度均通过速度维勒特积分器进行更新。在开展分子动力学模拟之前，先使用最速下降算法进行能量最小化处理，直到系统达到规定的收敛标准（能量变化低于10^-7电子伏，最大原子力低于10^-3电子伏/埃）。为实现热力学平衡，在300 K和零外压的等温等压（NPT）系综中进行了30 ps（30,000步）的模拟。温度和压力通过带有适当阻尼参数的Nosé-Hoover恒温器和压力控制器进行调控。平衡达成后，又在正则（NVT）系综中进行了50 ps（50,000步）的产生模拟，以收集统计数据。为了进行状态方程（E-V）计算，在NPT系综下进行模拟，以便原子位置和晶胞体积都能充分松弛。原子间相互作用的截断半径设定为6 ?。为验证所开发的深度势（DP）模型的预测精度，将其结果与包括MACE（补充材料S1）和NequIP（补充材料S2）在内的现有原子间势模型的结果进行了比较(41)、(42)。

2.3 训练过程
图1示意性地展示了构建MgLa合金深度势的工作流程。关键步骤包括：(i)生成初始训练数据集，(ii)进行迭代参数优化，(iii)通过主动学习循环逐步扩展配置集。初始训练集涵盖了多种元素、合金以及晶体结构类型，以确保数据的多样性和代表性。这一基础数据集的构成详细说明见补充材料S3。此外，2.1节（“DFT计算”）明确了用于初始分子动力学（AIMD）模拟的具体参数。随后，使用dpdata工具将AIMD轨迹转换为DeePMD-kit可用的格式。在该框架中，系统的总能量表现为原子能量的总和，每个原子能量由一个深度神经网络描述，该网络可将局部原子环境映射到相应的能量贡献值。每个原子的局部环境是通过在7.0 ?的截断半径范围内的平滑嵌入描述符（se_e2_a）来表示的。与传统将相互作用明确分解为二体或三体项的经验势不同，DP模型不依赖于预定义的成对或角度函数形式，而是利用深度神经网络隐式捕捉高阶多体相互作用。因此，DP模型可被视为完全的多体势。嵌入网络包含三个隐藏层，神经元数量分别为20、40和80，而拟合网络则包含三个全连接层，每层有200个神经元。所得数据随后按8:2的比例分为训练集和测试集。这一精心设计的工作流程最终产生了庞大的原始数据集，为后续模型开发提供了坚实基础。训练深度势模型的关键超参数记录在S4中。补充材料S5则阐述了潜在的优化算法，并详细说明了主动学习采样循环的理论依据。本研究采用DP-GEN框架结合随机梯度下降方法，通过最小化包含总能量和原子间力项的复合损失函数，从而获得高精度的原子势能面。借助预训练模型的特征空间，一种创新的主动学习策略能够从新生成的MD轨迹中筛选出最具信息量的结构快照，这种方法显著扩大了训练数据的相空间覆盖范围。详细方法见S6。经过DFT验证后，扩展后的配置集会经历多轮迭代训练，最终得到的模型具有很高的预测精度以及出色的跨系统适用性。这一集成框架为研究MgLa合金的微观力学行为提供了可靠、原子级精度的计算解决方案。更广泛地说，所提出的方法作为一种变革性的工具，有望用于开发各种多元合金体系的原子间势，从而为计算材料科学开辟新的方法路径。

下载：下载高分辨率图像（693KB）
下载：下载完整尺寸图像
图1. 基于深度学习和DFT参考数据，用于训练和拟合Mg-La合金深度势模型的示意图。

2.4 电化学模型
本研究聚焦于Mg-La合金体系的阳极极化行为。根据Butler-Volmer理论公式，描述阳极溶解的电流-电位关系由以下方程给出(43)、(44)：
(1)
其中，阳极反应的交换电流密度用以下表达式确定(45)：
(2)
在此公式中，n表示参与阳极半反应的电子数。对于金属体系，单原子溶解过程中形成的单价阳离子（M+）通常被视为短暂中间体。电荷转移系数α通常是通过将模型与实验测量值拟合来获得的。参数ΔG0描述了与原子溶解相关的活化自由能，大约占断裂原子间键并产生离子物种所需总能量的30%。该模型涉及的物理常数包括F、k、T、h和R。此处，Ue表示平衡电位，而/ρ则表示表面能量密度。在高极化条件下，施加的电位会偏离其平衡值，此时可忽略方程1中的第二项，从而得到简化表达式：
(3)
而在等温等压条件下，形成单位面积界面所需的最小不可逆功决定了该界面的表面能量，其数学表达式如下：
(4)
在此表达式中，表示板状结构的总能量，而则表示相应块体Mg合金的能量。

3 结果与讨论
3.1 晶体性质
3.1.1 深度学习模型的整体性能
为评估所构建DP模型的准确性和可靠性，我们首先在训练数据分布范围内，系统地比较了DP模型的预测结果与参考密度泛函理论（DFT）计算结果。图2展示了预测的总能量和原子力与DFT结果的散点图。所有数据点都紧密分布在参考线y=x周围，直观显示出两者之间存在高度的线性一致性。具体而言，能量和原子力预测的均方根误差分别仅为5.79 × 10-3 eV/原子和1.63 × 10-3 eV/?。这些关键指标共同表明，DP模型的预测精度超过了大多数已报道的传统机器学习势模型。为进一步评估该模型在更接近实际模拟条件下的泛化能力，我们通过从0 K、300 K和1200 K下的分子动力学轨迹中随机抽取200组配置，进而进行单点DFT计算，构建了一个独立的测试集。如图S1所示，该“未见过的”数据集上的DP模型预测结果与DFT结果仍保持强相关性（R2=0.98），其能量均方根误差与训练集中的数值相当。这直接证明了该模型出色的外推稳健性，以及其在有限温度下进行动力学模拟时的可靠性。总之，系统的基准测试验证了所构建DP模型的高预测精度和出色的泛化能力，为后续研究Mg-La合金的力学和腐蚀行为提供了可靠的原子尺度计算基础。

下载：下载高分辨率图像（545KB）
下载：下载完整尺寸图像
图2. Mg-La合金系统的DFT参考数据与深度势（DP）模型预测结果的比较：(a)总能量，(b-d)沿x、y和z方向的原子力分量。总能量的均方根误差为5.79 × 10-3 eV/原子，而原子力的均方根误差分别为1.65 × 10-3（x方向）、1.67 × 10-3（y方向）和1.57 × 10-3 eV/?（z方向）。

3.1.2 能量-体积（E-V）曲线
为评估所构建DP模型在连续变形条件下的可靠性，我们系统地计算了纯组分（Mg和La）以及MgLa合金的能量-体积（E-V）关系。如图3所示，DP模型预测的E-V曲线在整个应变范围内都与参考DFT结果高度吻合，整体均方根误差低于0.022 eV/原子。具体而言，MgLa的预测平衡体积为31.5 ?3/原子（图3c），与DFT参考值的相对偏差仅为0.23%。这些关键结构参数的高度一致性体现了DP模型在描述MgLa体系原子相互作用方面的出色可靠性和预测精度。此外，为量化这一微小偏差对力学性能的实际影响，我们采用了经典的Birch-Murnaghan状态方程结合Voigt-Reuss-Hill平均方案，进行了严格的敏感性分析，以研究“体积偏差”如何导致“模量误差”。分析表明，体积与弹性模量之间存在近似线性关系，即便体积偏差为0.23%，所带来的模量不确定性也仍在可接受范围内（<2%）(46)、(47)。这一结论与以往关于高熵合金和过渡金属合金的第一性原理及分子动力学研究结果一致，这些研究表明微观结构的微小波动会对宏观力学响应造成可控且可接受的扰动(48)。基于此，0.23%的体积偏差并不会影响本研究的根本稳健性。

下载：下载高分辨率图像（372KB）
下载：下载完整尺寸图像
图3. 通过DFT计算和深度势（DP）模型预测得到的Mg（a）、La（b）和MgLa（c）的能源-体积曲线，插图显示了相应的晶体结构。

3.1.3 晶格常数
基于已确定的E-V曲线，我们运用Birch-Murnaghan状态方程对MgLa体系的平衡晶格常数进行了拟合和测定。如表1所示，DP模型预测的MgLa合金晶格常数为4.016 ?，与DFT参考值3.978 ?非常接近，相对偏差仅为0.96%。这一精度水平超过了MACE和NequIP等传统势模型。这一结果进一步证明了所构建DP模型在关键晶体结构参数方面的高精度和可靠性。更重要的是，对平衡结构的精准预测为后续所有关于力学性能和腐蚀行为的原子尺度模拟提供了可靠的起点。DP模型在这一阶段的出色表现，为其在后续章节中模拟复杂使用场景的应用奠定了坚实的结构基础。

表1. 不同方法计算得到的Mg、La和MgLa的晶胞参数对比。
金属 | 晶格常数（?） | 来源
---|---|---
Mg | a=3.209，c=5.211 | expa=3.205，c=5.214
DFT | a=3.206，c=5.215 | DPa=3.201，c=5.209
MACE | a=3.202，c=5.194 | NEQUIPa=3.200，c=5.190
MEAM | la=5.323，c=5.323 | expa=5.323，c=5.323
DFT | a=5.248，c=5.248 | DPa=5.112，c=5.112
MACE | a=5.104，c=5.104 | NEQUIP
MgLa | expa=3.978，c=3.978 | DFTa=4.016，c=4.016
DPa=4.136，c=4.136 | MACEa=4.055，c=4.055
NEQUIP

3.1.4 声子谱
声子谱反映了材料中原子的内在振动特性，其谱特征决定了材料在热扰动下的晶格响应。为评估DP模型在描述晶格动力学方面的可靠性，我们对MgLa合金体系，使用了DP势和第一性原理密度DFT计算得到的声子色散关系进行了系统比较。声子色散曲线中出现虚频是判断动力学稳定性的重要标准，因为这类频率通常意味着可能存在结构弛豫或相变。如图4所示，在所研究的各个体系的整个布里渊区内均未观察到虚频。所有的声子分支都是连续的且始终为正，这表明它们的基态结构具有动力学稳定性。定量比较显示，对于纯Mg和纯La，DP势预测的声子频率与DFT结果非常吻合，整体平均偏差低于2.0%。即使在以复杂振动模式为特征的高频区域，最大频率差异也保持在4.9 cm-1以下。对于更为复杂的MgLa合金体系，DP模型同样能够在低频和高频范围内精确再现DFT的声子特征，其标准差不超过1.3 cm-1。这些结果表明，所采用的DP势在描述晶格动力学性质方面达到了近乎DFT的精度，为复杂合金的声子模拟提供了准确高效的计算工具。此外，这一能力还为先进生物医学镁基合金的高温稳定性进行原子尺度评估和设计奠定了可靠的基础。

下载：下载高分辨率图像（660KB）
下载：下载完整尺寸图像
图4. 在布里渊区的高对称路径上，Mg、La和MgLa的声子色散关系，对比了密度泛函理论（DFT）计算结果与深度势（DP）模型预测结果。图(a)-(c)分别对应Mg、La和MgLa。

3.2 热物理性质
3.2.1 弹性刚度矩阵
弹性常数是连接微观原子间相互作用与宏观力学响应的重要物理量。它们直接反映了材料内部相互作用对外部载荷的响应方式，因而决定了材料的力学行为和稳定性。为评估DP模型在预测多组分合金体系弹性性能方面的可靠性，研究人员根据DP模型生成的原子力计算出了弹性常数。如表2所示，使用DP势能计算出的所有弹性常数均满足与其相应晶体结构相对应的机械稳定性准则。具体而言，对于六方密排结构的纯镁，需满足以下条件：C11 > 0、C33 > 0、C44 > 0、C11 - C12 > 0，以及(C11 + C12)C33 - 2C132 > 0；对于立方结构的纯镧和MgLa合金，则需满足C11 > 0、C44 > 0、C11 - C12 > 0以及C11 + 2C12 > 0这些稳定性准则。这些结果表明，所有研究的体系在其基态下都具有机械稳定性。尤其对于作为结构材料的MgLa合金，其计算得到的剪切弹性常数C44为33.8 GPa，与DFT参考值（约34.0 GPa）极为接近，偏差仅约为0.2 GPa。这些结果充分体现了DP模型在预测多组分合金体系方面的高精度，为其在生物医学领域，如承重型镁基植入材料中的应用提供了坚实的理论支撑。

表2. 不同方法计算得到的镁、镧及MgLa合金的弹性常数对比

金属 | 弹性性质 | Exp (Gpa) | DFT (Gpa) | DP (Gpa) | MACE (GPa) | NEQUIP (Gpa)
---|---|---|---|---|---|---
镁 | C11 | 63.4 | 63.8 | 64.0 | 64.4 | 61.3
| C12 | 25.9 | 26.1 | 26.1 | 24.1 | 27.1
| C13 | 21.7 | 21.8 | 21.9 | 22.9 | 23.3
| C33 | 66.4 | 67.3 | 67.5 | 69.3 | 70.6
| C44 | 18.4 | 18.6 | 18.7 | 20.2 | 21.6

镧 | C11 | 33.0 | 35.0 | 33.7 | 38.5 | 38.6
| C12 | 14.0 | 14.0 | 14.8 | 12.4 | 12.4
| C44 | 18.0 | 15.0 | 16.4 | 13.1 | 13.1

MgLa合金 | C11 | 47.0 | 48.0 | 49.2 | 47.7 | 47.5
| C12 | 27.0 | 28.0 | 25.3 | 25.7 | 26.0
| C44 | 35.0 | 34.0 | 33.8 | 35.8 | 36.4

2. 体积模量与剪切模量
在得到弹性刚度张量后，研究人员采用经典的Voigt-Reuss-Hill（VRH）平均方法计算出了体积模量（B）和剪切模量（G）。VRH方法通过对上界（Voigt）和下界（Reuss）数值进行平均，能够可靠地估算出多晶材料的有效弹性模量，这一方法在第一性原理计算和分子动力学研究中都被广泛应用。体积模量和剪切模量的Voigt与Reuss上下界分别由公式（5）和（6）给出，其中分别代表在Reuss和Voigt边界条件下计算得到的体积模量，而则对应着相应的剪切模量。如表3所示，对于MgLa合金，DP势能预测的体积模量为B=33.27 GPa，剪切模量为G=22.29 GPa，这两项数值都与DFT计算结果极为接近，偏差均小于1.42 GPa。相比之下，传统的MEAM势能在预测这些宏观弹性模量时会出现更大的偏差。这些结果表明，所构建的DP势能在预测MgLa合金的宏观弹性模量（B和G）方面，其精确度可与第一性原理计算相媲美，这也证明了DP模型能够可靠地捕捉该合金体系的本质机械响应。

表3. 不同方法计算得到的镁、镧及Mg-La合金的弹性模量、剪切模量、杨氏模量及泊松比

金属 | 方法 | 模量 | B | G | ν
---|---|---|---|---|---
镁 | Exp | 36.9 | 0 | 23.2 | 30 | 54.2 | 0 | 0 | 0.255
| DFT | 36.3 | 6 | 0 | 21.6 | 50 | 54.3 | 0 | 0 | 0.255
| DP | 36.8 | 0 | 21.7 | 0 | 54.4 | 0 | 0 | 0.254
| MACE | 38.0 | 0 | 22.5 | 0 | 56.2 | 0 | 0 | 0.253
| Nequlp | 36.0 | 0 | 24.5 | 0 | 59.9 | 0 | 0 | 0.223
镁 | MEAM | 37.5 | 0 | 19.7 | 0 | 55.8 | 0 | 0 | 0.257

镧 | Exp | 20.3 | 33 | 0 | 13.9 | 30 | 5 | 0 | 34.0 | 50 | 0 | 0.221
| DFT | 21.0 | 0 | 36 | 0 | 13.0 | 30 | 5 | 0 | 32.3 | 30 | 0 | 0.243
| DP | 21.1 | 0 | 36 | 0 | 13.1 | 30 | 5 | 0 | 32.6 | 60 | 0 | 0.242
| MACE | 21.1 | 0 | 36 | 0 | 13.0 | 30 | 5 | 0 | 32.5 | 40 | 0 | 0.243
| Nequlp | 21.1 | 0 | 36 | 0 | 13.1 | 30 | 5 | 0 | 32.6 | 20 | 0 | 0.243

MgLa合金 | Exp | 33.6 | 70 | 21.2 | 50 | 30 | 5 | 0 | 52.6 | 80 | 0 | 0.239
| DFT | 34.6 | 70 | 20.8 | 50 | 30 | 5 | 0 | 52.1 | 70 | 0 | 0.249
| DP | 33.2 | 70 | 22.2 | 50 | 30 | 5 | 0 | 54.6 | 50 | 0 | 0.226
| MACE | 33.0 | 70 | 22.3 | 50 | 30 | 5 | 0 | 54.7 | 30 | 0 | 0.224
| Nequlp | 33.1 | 70 | 22.3 | 50 | 30 | 5 | 0 | 54.8 | 20 | 0 | 0.225

3.2. 杨氏模量与泊松比的验证
杨氏模量是描述材料轴向刚度的核心力学参数，而泊松比则反映了材料在单轴应力作用下的横向应变与轴向应变的比值。二者共同构成了判断材料弹性行为的重要指标。为验证DP模型在计算弹性性能方面的准确性，研究人员采用了以下公式，对包括纯镁、纯镧以及MgLa合金在内的多种体系的杨氏模量与泊松比进行了测算：

（7）（8）

测试结果显示，DP模型预测的MgLa合金杨氏模量与DFT参考值极为接近，相对误差较低。相比之下，MACE和NequIP模型预测的杨氏模量与DFT基准值的偏差较大，其中NequIP的相对误差甚至高达5.2%。至于泊松比，DP模型为MgLa合金算出的值为0.226，与DFT的0.249相比也相当吻合。同样，传统的MEAM势在预测杨氏模量和泊松比时也会出现较大的偏差。综合来看，所构建的DP势在再现DFT级别的宏观弹性性能方面具有很高的可靠性，其预测精度明显优于MACE和NequIP模型，这充分体现了DP势在高效、精准地分析复杂合金体系弹性性能方面的优势。

3.3. 能量特性
3.3.1 缺陷形成能
为评估DP模型描述局部结构畸变和键合重排能力，研究人员利用该模型计算了关键的缺陷形成能。为系统验证其预测精度，研究人员在相同的计算条件下，将DP模型与包括采用PBE泛函的第一性原理DFT计算、MACE模型以及NequIP模型在内的其他代表性方法进行了对比。如图5a和5b所示，DP模型预测的间隙形成能与DFT-PBE参考值非常接近，误差低于0.78%，显示出极高的准确性。同样，DP模型预测的空位形成能也表现优异，偏差小于0.99%。相比之下，MACE模型的预测相对误差更大，可达2.07%，而NequIP模型的最大偏差甚至达到2.88%，表明其预测存在更为明显的偏差。作为对比，传统MEAM势为纯镁计算的间隙形成能和空位形成能分别为2.075电子伏特和1.266电子伏特（详见表S3）。这些结果表明，与传统势函数以及一些现有的机器学习势函数相比，DP模型在预测缺陷能量方面具有更出色的精度，尤其是在处理复杂的缺陷构型时。

下载：下载高分辨率图像（406KB）
下载：下载完整尺寸图像
图5. 镁、镧及MgLa体系缺陷能和表面能的对比：(a)间隙形成能；(b)空位形成能；(c)表面能；(d)内聚能。图中柱状图上的数值为DFT计算结果与深度势模型、MACE模型以及NequIP模型预测值的对比结果。

3.3.2 表面能与结合能
此外，研究人员还进一步研究了MgLa合金的表面能和结合能。如图5c和5d所示，DP模型预测的值与DFT-PBE参考值非常接近，相对误差分别仅为0.3%和0.12%。相比之下，传统MEAM势为纯镁的表面能和结合能预测出的数值偏差较大（详见表S3）。这些结果表明，DP模型能够准确描述体相（内聚能）和表面（表面能）环境中的能量状态，而这两种环境代表了截然不同的原子配位环境，这也体现了该模型在描述复杂合金体系中的势能分布时的稳定性和可靠性。相反，MACE和NequIP等其他基准模型则存在较大的误差，其预测值在不同体系或不同构型下波动较大，反映出它们在捕捉势能面的细微特征以及实现跨构型推广方面的局限性。这一现象与计算材料科学领域的普遍认知一致：要可靠地预测包含表面、界面和缺陷的复杂异质体系的能量特性，势函数模型需要包含足够且严谨的物理约束，以确保其具备出色的泛化能力。为提升DP模型在Mg-La体系中的适用性，研究团队在训练数据集中不仅包含了平衡态晶体结构，还纳入了多种非平衡态的原子构型，比如存在晶格应变的结构、受热扰动的结构以及其他局部畸变的体系。此外，本研究还采用了DP-GEN主动学习迭代框架。通过在不同热力学条件下探索构型空间，该框架可以识别出预测不确定性较高的候选结构，再通过DFT计算对其准确标注，随后将其重新纳入训练集。这种迭代式的处理方式能够不断扩展训练数据的范围，使其涵盖更多远离平衡态的原子环境，从而显著提升模型在实际模拟中的稳定性和泛化能力。总体而言，DP势在描述MgLa合金中与界面和缺陷相关的微观能量特性方面，展现了与DFT相当的准确度和可靠性。这为模拟表面反应和腐蚀的初始阶段提供了可靠的计算基础，同时也凸显了该模型在先进合金设计及性能预测方面的巨大潜力。

3.4 电化学特性
3.4.1 计算方法
为建立从原子尺度特性到宏观腐蚀性能之间的预测桥梁，研究人员开发了一种集成多尺度计算框架。该框架的核心思路如下：首先，利用经过DFT基准测试验证且结合分子静力学方法的DP模型，计算出不同晶面的表面能等关键原子尺度参数，从而实现高通量评估；之后将这些参数引入经典的Butler-Volmer电化学动力学模型中，为高通量筛选耐腐蚀成分提供可行方法。该框架被应用于MgLa合金体系，用于高效筛选其耐腐蚀组分。在该框架中，研究人员引入了表面能密度这一概念，用以量化不同晶面在腐蚀过程中的相对敏感程度。如表S3和表S4所示，随着成分和晶面取向的变化而变化的、通过DP模型计算得到的表面能和表面活化能，与实验观测到的腐蚀趋势高度吻合，这证明了这一表征方法的有效性。需要指出的是，在电化学腐蚀的理论描述中，材料的功函数（或电子发射势）通常被视为一个重要影响因素。然而，最新研究表明，在决定阳极极化曲线Tafel段的高过电位区域这一关键动力学区域内，功函数的影响可能并非最为重要。本研究开发的框架侧重于将微观表面稳定性与宏观腐蚀速率联系起来，这一思路更强调界面反应动力学，而非体相的电子特性。

3.4.2 具有特定成分的三种Mg-La合金的极化曲线
本研究将深度势方法应用于阳极溶解过程，从而能够预测MgLa合金在电化学条件下的极化曲线。图6展示了几种MgLa合金体系在中性溶液（pH = 7）中，其最密排晶面（Mg2La的(110)晶面、Mg7La的(001)晶面、MgLa2的(010)晶面）的阳极极化曲线，同时还给出了相应的析氢反应动力学预测结果（详见补充信息S7和图S2）。通过结合分析阳极极化曲线和析氢反应动力学，研究人员计算出了纯镁、纯镧以及其他合金体系的腐蚀电位和腐蚀电流密度。系统的研究结果表明，合金的成分比例会影响到这些关键的腐蚀参数。值得注意的是，模拟得到的腐蚀参数与实验测量值高度吻合。例如，在3.5重量百分比的NaCl溶液中，镁的(001)晶面的模拟腐蚀电流密度为?3.408，腐蚀电位为?1.277伏特，这两个数值都在实验观测的范围之内（log[io]：-4.265至约-4.012；Ecorr：-1.266至-1.172）（49）、（50）。计算结果与以往实验研究报道的趋势总体相符，进一步证明了基于DP方法的电化学预测框架的准确性。而对于Mg7La合金，其计算得到的腐蚀电位为?1.256伏特，腐蚀电流密度为10^-3.526安培/平方厘米；相比之下，MgLa2的腐蚀电位更高，为?1.112伏特，腐蚀电流密度则更低，为10^-4.782安培/平方厘米。这些结果表明，纯镁的腐蚀倾向明显强于纯镧，这是因为镧具有更强的表面钝化能力。需要说明的是，我们的计算是基于最密排的镁晶面进行的，而实验测量则是针对含有多种晶面取向的多晶样品进行的，这也是导致计算值与实验值之间存在较小偏差的原因。在所研究的MgLa合金体系中，MgLa2的耐腐蚀性最强，而Mg7La的腐蚀倾向最为严重。其背后的微观机制如下：较高的镧含量有助于形成稳定且致密的钝化膜，从而有效抑制阳极溶解；而较高的镁含量则会提升合金表面的电化学活性，加速金属的溶解速度。总而言之，这一将阳极溶解模型与深度势方法相结合的框架，在合理设计耐腐蚀合金以及预测其性能方面展现出卓越的效果和巨大潜力。

下载：下载高分辨率图像（362KB）
下载：下载完整尺寸图像
图6. 中性溶液（pH = 7）中Mg-La合金的阳极极化曲线：(a) Mg2La，(b) Mg7La，(c) MgLa2。所有电位值均以标准氢电极为参考。图中的标记分别表示腐蚀电位（E(corr)）以及对应的腐蚀电流密度（i(corr)）。

3.4.3 验证Mg-La合金特定晶面的极化曲线
为更深入地了解MgLa合金腐蚀的动力学特性，本研究采用了大规模计算模型来研究其宏观电化学极化行为。如图所示S3中，针对两种不同的MgLa合金成分，构建了三个低指数晶面的原子尺度模型。需要指出的是，为简化模型并聚焦于合金的固有腐蚀特性，这些表面被假设为理想的、完整的且无氧化物存在。通过将表面能等固有参数纳入改进后的界面反应模型，得到了MgLa合金的极化曲线（见图7）。结果表明，Mg3La合金的腐蚀行为显著依赖于晶面取向。具体而言，(100)和(110)晶面的阳极电流密度分别为10^-3.708和10^-3.694 A/cm2，腐蚀电位分别为-1.243 V和-1.245 V；而(111)晶面的阳极电流密度为10^-3.724 A/cm2，腐蚀电位为-1.241 V，表明其在所研究的晶面中具有最佳的耐腐蚀性。在该建模框架下，不同晶面之间的阳极电流密度和腐蚀电位差异主要由表面能密度决定，而表面能密度又由晶面的活化能和原子排列方式决定。较高的表面能密度会显著降低阳极溶解的能量障碍，加快反应速率，提高腐蚀电流密度，并降低腐蚀电位。因此，Mg3La的腐蚀活性顺序为(110) > (100) > (111)，对应的阳极电流密度比值为icorr(111)/icorr(100) = 0.964，icorr(111)/icorr(110) = 0.933，这与表面能密度的理论排序完全一致。这一趋势在MgLa3中也得到了进一步验证：(010)、(100)和(001)晶面的腐蚀电流密度分别为10^-4.925、10^-4.775和10^-4.715 A/cm2，相应的腐蚀电位则分别为-1.089、-1.108和-1.115 V，呈单调上升趋势。在MgLa体系中，(010)晶面的计算阳极电流密度仅为(100)和(111)晶面的71%和62%。此外，Liu等人测得的MgLa3的阳极电流密度（约10^-5 A/cm2）也与这些预测值十分接近。总之，深度势模型与Butler-Volmer框架的结合，为理解和预测多组分合金的腐蚀行为提供了一种强大的新方法。图7显示了Mg3La与MgLa3的腐蚀行为相似，说明两者之间不会发生电偶腐蚀。因此，将这两种合金结合使用可兼顾性能与成本。基于大规模原子模型的精确极化曲线预测，标志着腐蚀研究正从传统的实验表征向以计算驱动、以设计为导向的方向发展。

下载：下载高分辨率图像（317KB）
下载：下载全尺寸图像
图7. 中性溶液（pH=7）中Mg-La金属间化合物的动电势极化曲线。(a) Mg3La：(111)、(100)、(110)；(b) MgLa3：(010)、(100)、(111)。插图标出了E(corr)和i(corr)，用于比较不同晶面下的初始腐蚀情况。虚线表示H+/H2平衡电位；所有电位均以标准氢电极为参考并经过pH校正。

4. 结论
在本研究中，通过第一性原理计算与深度学习方法相结合，系统地开发并验证了一种高精度的Mg-La合金深度势模型，该模型能够实现晶体结构、力学性能以及电化学腐蚀行为的跨尺度研究。基准测试结果表明，该深度势模型能够高度准确地预测包括能量、原子力及弹性响应在内的关键物理量，其能量和原子力预测误差分别控制在10^-3 eV/原子和10^-3 eV/埃的数量级，为后续的大规模分子动力学模拟提供了可靠的势能面。在结构与力学性能方面，该模型能够准确再现Mg-La合金的能量-体积关系、平衡晶格常数以及声子色散特性。同时，它还能成功预测弹性常数、体模量、剪切模量和杨氏模量，所有结果均满足严格的力学和动态稳定性要求，并与第一性原理计算及实验数据高度吻合。在能量特性方面，该模型对于缺陷形成能、表面能和内聚能等关键指标的相对误差均低于1%，远优于传统的经验势函数以及一些现有的机器学习原子间势函数。更为重要的是，从该深度势模型中得到的表面能等原子尺度能量描述符被引入Butler-Volmer电化学动力学框架中，从而建立了原子尺度表面稳定性与宏观阳极溶解动力学之间的明确定量关联。基于这种多尺度耦合，直接预测出了Mg-La合金的阳极极化曲线，其腐蚀电流密度约为10^-4.9 A/cm2，与现有实验数据在数值和趋势上均十分吻合。总体而言，结合深度势模型与电化学动力学的多尺度仿真方法，为模拟复杂镁基合金的腐蚀行为提供了一种具有物理意义且可定量验证的途径。该方法不仅深化了对Mg-La合金中原子尺度能量特征与电化学响应之间内在关系的理解，还为其他多组分金属体系的结构-性能-腐蚀关系的计算研究提供了通用性的参考框架。

CRediT作者贡献声明：
张祥阳：方法学。
荣菊：项目管理。
何志杰：写作——审阅与编辑。
余晓华：写作——审阅与编辑、项目管理、资金获取。
张浩翔：写作——审阅与编辑、写作——初稿撰写、方法学、研究设计、正式分析、概念构建。
高若桐：软件应用。
陈佳辉：数据整理。

数据可用性：
所有数据均已包含在补充信息中。