机器学习势函数：构建通用模型加速非晶合金材料设计与性能预测

1. 项目概述：用机器学习“势函数”为无序世界建模

在材料科学，尤其是非晶合金（或称金属玻璃）的研究中，我们长期面临一个核心矛盾：我们渴望从原子层面理解并预测材料的性能，但计算工具要么太“笨”，要么太“贵”。

传统的经验势函数，比如嵌入原子法（EAM），计算速度飞快，可以模拟上百万个原子，但它的精度天花板很低，本质上是对原子间相互作用的简化拟合，对于成分复杂、结构无序的非晶合金，其预测结果常常与实验或高精度计算相去甚远。而作为“金标准”的第一性原理计算，比如基于密度泛函理论（DFT）的方法，精度极高，能揭示电子层面的相互作用，但其计算成本是天文数字。模拟一个几百个原子的体系可能就需要数天甚至数周的计算资源，对于需要统计意义、包含数千原子且结构无序的非晶体系，DFT几乎寸步难行。这就好比，你想研究一座城市的交通流，EAM给了你一张模糊的卫星照片，而DFT则要求你精确追踪每一辆车的发动机状态——前者看不清细节，后者根本算不过来。

机器学习原子间势函数（MLIP）的出现，正是为了弥合这道鸿沟。它的核心思想非常巧妙：我们不直接去解复杂的量子力学方程，而是用机器学习模型（通常是深度神经网络）去“学习”DFT计算产生的大量“标准答案”——即不同原子构型下的能量、原子受力等信息。一旦模型训练完成，它就能像一个经验丰富的“插值专家”和“外推智者”，在面对新的、未见过的原子排列时，以接近DFT的精度，瞬间给出能量和力的预测，而计算成本仅比经验势函数略高。

这个项目的目标，就是为非晶合金这个特殊的材料家族，打造一个“通用型”的MLIP。为什么强调“通用”？因为过去的研究大多针对某个特定体系（如Cu-Zr）开发专用势函数，换个成分就得从头再来，费时费力。我们希望能建立一个覆盖多种元素、多种成分的“大一统”模型，让研究人员能像使用一个标准工具包一样，快速对各类非晶合金进行可靠模拟。

2. 核心思路与方案设计：如何构建一个“通才”势函数

2.1 目标定义：从“专用工具”到“通用平台”

我们的核心目标不是为某一个明星非晶合金（比如Pd40Ni40P20）做一个完美的势函数，而是构建一个能覆盖19个代表性二元、三元非晶合金体系，涉及25种化学元素的通用模型。这一定位决定了我们整个技术路线的基石：

1. 广度优先的数据策略：我们选择的19个体系（如Ca-Al, Pd-Si, Cu-Zr, Mg-Cu-Y等）并非随意挑选。它们是非晶合金领域的“经典家族”，涵盖了从轻质（Mg基）到高强（Co基、Ir基），从低玻璃转变温度到高玻璃转变温度的广泛性能谱。这确保了模型学习到的是非晶态中普适的原子相互作用模式，而非某个体系的特例。
2. 预训练-微调范式：从头训练一个覆盖如此多元素的深度势函数模型，需要海量的数据和巨大的算力。我们采用了更高效的策略：基于一个已在庞大晶体材料数据库（Materials Project）上预训练好的模型——CHGNet——进行微调。CHGNet本身已经学习了89种元素的“化学直觉”，我们只需要用非晶合金的专用数据去“校正”它对于无序结构的认知偏差。这好比一位精通多种语言的翻译（预训练模型），我们只需要对他进行特定领域（非晶合金）的术语和文体培训（微调），他就能出色地完成该领域的翻译工作，事半功倍。

2.2 数据集构建：为无序结构“拍照取样”

数据集的质量直接决定了MLIP的成败。对于非晶合金，其原子排列长程无序、短程有序，如何构建一个能充分反映其结构多样性的数据集是关键。

我们的构建流程可以概括为“熔融-快冷-采样”：

1. 成分选择：在每个二元/三元体系中，按原子百分比梯度选取成分点（如A20B80, A40B60…）。这确保了模型能学习到成分变化对原子环境的影响。
2. 第一性原理分子动力学（AIMD）模拟：对每个选定的成分，使用VASP软件进行高温（2500 K）下的AIMD模拟。高温使得原子剧烈运动，能有效遍历更多的可能构型。
3. 构型采样：从AIMD模拟的轨迹中，按固定时间间隔“快照”出原子的瞬时位置。每个成分我们采集数百至上千个这样的“快照”。每个快照都对应一个由DFT计算出的精确能量和每个原子所受的力，这就是我们模型的“标准答案”。
4. 数据增强：为了进一步提升模型的鲁棒性，我们对每个采样构型的晶胞施加均匀的体积缩放（如0.95, 0.98, 1.02, 1.05倍）。这相当于让模型学习同一原子排列在不同压力（或密度）下的响应，极大地增强了其对于体积变化的预测能力。

最终，我们构建了一个包含20,400个构型的训练集和一个独立的5,100个构型的测试集。这个数据集的规模和多样性，是支撑通用型MLIP的基础。

实操心得：数据集的“代表性”陷阱构建数据集时，最容易犯的错误是只采集平衡态或能量最低的构型。对于非晶合金，我们必须刻意包含高温、高能态的构型，因为快冷形成玻璃的过程本身就穿越了这些非平衡态。如果数据集只包含完美的晶体或松弛后的玻璃态，训练出的模型将无法正确模拟熔融、扩散、形变等动力学过程。我们的策略是在高温AIMD中采样，正是为了捕获这些至关重要的“高能态”信息。

2.3 模型训练与评估：让“通才”胜过“专才”

我们基于CHGNet预训练模型，使用我们的非晶合金数据集进行微调。损失函数采用预测能量、原子受力与DFT计算值之间的均方误差（MSE）。

评估结果令人振奋：

• 精度：在独立测试集上，模型对原子能量的预测平均绝对误差（MAE）仅为5.06 meV/atom，对原子受力的MAE为128.51 meV/Å。这个精度已经非常接近DFT计算本身的误差范围，足以支撑可靠的分子动力学模拟。
• 通用性验证：我们做了一个关键对比：将我们用一个模型统一训练所有19个体系数据得到的“通用模型”，与为每个体系单独训练一个“专用模型”进行比较。结果显示，通用模型在大多数体系上的预测误差低于或等同于专用模型。这彻底颠覆了“一个模型专攻一个体系”的传统认知，证明了通过足够广泛和高质量的数据，完全可以训练出一个性能优异的“通才”势函数。

这背后的逻辑是，不同非晶合金体系间共享着许多相似的局部原子团簇（如二十面体）和键合模式。通用模型通过接触更多样化的数据，反而能更好地归纳出这些底层物理规律，其泛化能力更强。而专用模型如果数据稍有不足，就容易过拟合到特定体系的噪声上。

3. 势函数应用：从原子模拟到宏观性能预测

训练出一个高精度的势函数只是第一步，更重要的是用它来解决实际问题。我们利用训练好的MLIP，通过LAMMPS软件进行大规模分子动力学模拟，预测非晶合金的三个关键宏观性能：密度、杨氏模量和玻璃转变温度（Tg）。

3.1 模拟流程详解

1. 非晶结构生成（密度计算）：

   • 步骤：模拟实验的铜模铸造快冷过程。将体系在2500 K下熔化10 ps，然后在22 ps内快速冷���至300 K（室温），最后在300 K下弛豫10 ps，得到玻璃态结构。
   • 关键：通过分析得到的径向分布函数（RDF），确认其具有非晶态典型的“劈裂的第二峰”，以验证生成的确实是玻璃态而非晶体。
   • 密度获取：在300 K下进行等温等压（NPT）模拟50 ps，取最后10 ps的平均体积计算密度。

2. 杨氏模量计算：

   • 步骤：对弛豫好的非晶结构，分别沿X、Y、Z三个方向施加单轴拉伸应变，应变速率需足够慢以保证准静态条件。
   • 计算：记录应力-应变曲线，在弹性变形阶段（通常取应变0.5%-2%区间）计算曲线的斜率，即为该方向的杨氏模量。
   • 结果：取三个方向模量的平均值作为该材料的杨氏模量，以消除各向异性的影响。

3. 玻璃转变温度（Tg）确定：

   • 原理：玻璃在升温过程中，体积-温度曲线会在Tg处发生斜率变化。
   • 步骤：在NPT系综下，从远低于预估Tg的温度开始，以一定升温速率（在模拟中，这个速率远高于实验）进行一系列模拟。在每个温度点平衡后，记录体系的平均体积。
   • 拟合：分别对低温区（玻璃态）和高温区（过冷液体态）的体积-温度数据进行线性拟合，两条拟合直线的交点对应的温度即为模拟得到的Tg。

3.2 预测结果与实验对比

我们将模拟结果与文献中报道的实验数据进行了系统对比：

• 密度：预测值与实验值吻合极好，平均相对误差仅约3%。这说明MLIP能够非常准确地描述非晶合金在室温下的原子堆积状态。
• 杨氏模量：模拟值普遍高于实验值约16%。这是一个系统性偏差，主要源于两方面：
  1. 尺度效应：尽管我们模拟了1000个原子，远超DFT的尺度，但与实验样品（包含约10^23个原子）相比仍微不足道。小尺寸模拟体系缺少宏观材料中的缺陷、自由体积分布等，可能导致模量偏高。
  2. 测量方法差异：实验测量杨氏模量有超声波共振、纳米压痕、单轴拉伸等多种方法，不同方法本身就会给出略有差异的结果。
• 玻璃转变温度（Tg）：模拟值系统性地高于实验值约125 K。这主要归因于模拟的升温速率远高于实验。实验中常用0.33 K/s的慢速升温以保证均匀加热，而分子动力学模拟受限于时间尺度，升温速率往往在10^10 K/s量级。更高的冷却/升温速率会导致测得的Tg向高温方向移动，这是玻璃物理中的已知现象。

重要提示：系统误差的校正对于杨氏模量和Tg的系统性偏差，我们不必视其为模型的失败。相反，一旦通过少量实验数据标定出这个偏差（如模量乘以0.86，Tg减去125 K），校正后的MLIP预测就能与实验值高度吻合。这实际上为高通量筛选提供了完美工具：我们可以用MLIP快速计算成千上万个成分的“模拟模量”和“模拟Tg”，然后通过校正公式映射到“预测实验值”，从而高效锁定性能最优的成分区间，极大缩小实验试错的范围。

3.3 方法学对比：MLIP的独特优势

为了凸显MLIP的价值，我们以经典的非晶合金Cu50Zr50为例，对比了三种方法：

1. MLIP（本工作）：1000原子体系，精度高，计算成本适中。
2. EAM经验势：可模拟百万原子，计算速度最快，但精度最低，无法准确反映成分变化对性能的细微影响。
3. 第一性原理计算（DFT）：精度最高，但只能算100原子左右的小体系。对于非晶这种无序体系，100个原子的统计代表性严重不足，导致其预测的宏观性能（如模量）反而可能不准确。

结论非常清晰：MLIP在计算精度和模拟尺度/效率之间取得了最佳平衡。它既能达到接近DFT的精度，又能以接近经验势函数的效率进行数千原子规模、纳秒时间尺度的模拟，这正是研究非晶合金所必需的。

4. 微观机理解析：破解实验异常之谜

MLIP不仅是一个预测工具，更是一把打开非晶合金“结构-性能”黑箱的钥匙。我们通过一个典型案例展示了这一点。

在优化Co-Ta-B非晶合金性能时，根据传统的“混合法则”经验，用杨氏模量更高的Ir元素部分替代Co，预期会提升合金的整体模量。但实验结果却显示，随着Ir含量增加，模量和Tg均下降，与经验预测相反。

我们的MLIP模拟完美复现了这一反常趋势！不仅如此，通过分析模拟得到的原子结构，我们揭示了其微观机理：

1. 径向分布函数（RDF）分析：发现随着Ir加入，第一近邻峰逐渐向更大的原子间距方向移动。
2. 原子堆积效率（APE）计算：显示Ir的加入使原子堆积变得松散。
3. 局域配位环境分析：发现Ir的加入并未显著改变Co、Ta、B等主要元素的局域配位数。

这些微观结构分析表明：Ir的加入（与Co同族，尺寸较大）导致整体原子结构变得更松散，从而降低了材料的刚度和热稳定性（模量和Tg下降）。但同时，由于Ir与Co化学性质相似，这种替代没有严重破坏玻璃形成能力。这解释了为什么性能变化趋势与基于纯元素性质的简单混合法则预测不符——微观结构的集体效应起了主导作用。

这一案例的强大之处在于：MLIP不仅预测了与实验一致的反常趋势，还提供了原子尺度的“现场观测”，让研究人员能够“看到”并理解变化背后的物理根源。这是单纯的数据拟合模型或经验势函数无法做到的。

5. 实操指南与经验总结

5.1 如何复现与使用此MLIP

我们的模型、代码和数据集已在GitHub上开源。如果你想在自己的研究中使用或借鉴此方法，可以参考以下步骤：

1. 环境准备：

   • 安装Python科学计算环境（推荐Anaconda）。
   • 安装机器学习框架（如PyTorch）和科学计算软件（VASP用于生成数据，LAMMPS用于模拟）。
   • 安装CHGNet库 (pip install chgnet) 以及其依赖。

2. 数据准备与训练：

   • 如果你有自己的体系，需要按照我们描述的方法（2.2节）生成DFT数据集。关键是确保构型的多样性和代表性。
   • 数据格式需转换为CHGNet支持的格式（如Atoms对象并包含energy,forces,stress等属性）。
   • 加载我们提供的预训练模型或CHGNet官方模型作为起点。
   • 划分训练集、验证集和测试集（例如8:1:1）。
   • 定义优化器（如Adam）、学习率（初始可设为1e-3）和损失函数（MSE），进行微调训练。注意监控验证集误差，防止过拟合。

3. 模拟与预测：

   • 将训练好的模型转换为LAMMPS可调用的格式（CHGNet支持导出为CHGNetCalculator，并可进一步封装）。
   • 在LAMMPS输入脚本中，通过pair_style和pair_coeff命令调用该势函数。
   • 按照3.1节描述的流程编写LAMMPS模拟脚本，进行结构弛豫、退火、力学性能测试等。

5.2 常见问题与避坑指南

1. 问题：训练时损失震荡不收敛，或验证集误差远大于训练集。

   • 排查：首先检查数据集质量。构型是否足够多样？能量和力的DFT计算是否收敛？数据��是否有异常值（如未完全弛豫的高能构型）？
   • 解决：尝试降低学习率，或使用学习率调度器（如ReduceLROnPlateau）。增加验证集的比例，确保其能代表数据分布。对于非晶体系，确保高温采样充分。

2. 问题：模拟过程中能量爆炸或原子飞散。

   • 排查：这是MLIP使用中最常见的问题。根本原因是模型遇到了训练数据分布之外的“陌生”原子构型（即外推失效）。
   • 解决：
     • 回溯训练数据：检查当前模拟的原子间距、配位环境是否远超训练集的范围。尝试在模拟中减小时间步长。
     • 增强数据：在问题构型附近进行额外的DFT计算，将新数据加入训练集重新微调模型。
     • 使用集成或不确定性估计：一些先进的MLIP框架可以提供预测的不确定性。当不确定性过高时，提示此处可能不可靠。

3. 问题：模拟得到的性能趋势正确，但绝对值与实验有系统偏差（如我们的模量和Tg）。

   • 排查与解决：这通常是尺度效应和模拟条件与实验条件不匹配导致的，不一定是模型本身的错误。如我们所述，需要通过少量实验数据对模拟结果进行系统校正。将校正后的关系作为高通量筛选的标尺，而非追求绝对值的完全一致。

4. 问题：想将模型用于训练集之外的新元素或新体系。

   • 策略：我们的通用模型已涵盖25种元素，对新体系有一定的外推能力，尤其是通过“相似元素替代”生成的成分。但如果引入全新元素（如训练集中完全没有的氧、氮等），模型性能会急剧下降。
   • 建议：采用“主动学习”或“增量学习”策略。先用现有模型对新区间进行初步模拟，识别出模型不确定性高的构型，对这些构型进行DFT计算，然后将新数据加入训练集，对模型进行进一步微调。这是一个迭代优化过程。

5.3 未来拓展与个人体会

这项工作构建的不仅仅是一个势函数，更是一个数据驱动的非晶合金研究新范式。我个人最深的一点体会是：在材料科学中，一个兼具精度与效率的计算工具，能从根本上改变研究范式。它让“计算指导实验”从一句口号变成了日常实践。

基于这个通用MLIP，后续可以开展许多有趣的工作：

• 高通量成分筛选：自动生成成千上万个候选成分，用MLIP快速预测其密度、模量、Tg甚至断裂韧性（通过模拟拉伸），从中筛选出最有潜力的少数几个进行实验验证，将研发周期从“年”缩短到“月”。
• 复杂性能预测：进一步模拟疲劳、腐蚀、磨损等更复杂的性能，这些在实验中耗时耗力，但在模拟中可以通过设计合理的“计算实验”来探究。
• 动力学过程研究：研究非晶合金的变形机制（剪切带萌生与扩展）、晶化过程、扩散行为等，这些都需要在足够大的时间和空间尺度上观察，MLIP提供了可能。
• 构建非晶合金专用数据库：将不同成分-结构-性能的模拟结果系统化，构建非晶合金的“材料基因”数据库，为更高级的机器学习（如生成模型、强化学习）提供燃料。

最后，对于想进入这个领域的研究者，我的建议是：从理解物理和化学直觉开始，而不是从调参开始。一个好的MLIP，其灵魂在于高质量、物理意义明确的数据集。你需要清楚你想要模型学习什么（是平衡性质还是反应路径？），然后有针对性地去生成数据。机器学习是强大的放大器，但它无法弥补对基础科学问题理解的缺失。将扎实的材料物理知识、清晰的科学问题，与强大的MLIP工具相结合，才能真正加速新材料的发现与设计。