机器学习势函数微调：精准预测卤化物固态电解质离子电导率

1. 项目概述：当机器学习势函数遇上卤化物固态电解质

在固态电池的研发竞赛中，卤化物固态电解质（如Li3YCl6-xBrx）因其与高电压正极的良好兼容性、较高的室温离子电导率（>1 mS cm⁻¹）和优异的机械稳定性，成为了明星候选材料。然而，要真正理解并优化这类材料的离子传输性能，我们需要在原子尺度上“看清”锂离子的迁移路径和动力学过程。传统的第一性原理分子动力学（AIMD）虽然精度高，但其巨大的计算成本（模拟纳秒尺度可能需要数周甚至数月）严重限制了我们对大体系、长时标扩散行为的探索。

这正是机器学习势函数（Machine Learning Interatomic Potential, MLIP）大显身手的舞台。它就像一个“学霸”，通过学习海量的密度泛函理论（DFT）数据，掌握了原子间相互作用的“规律”，从而能以接近DFT的精度，实现比DFT快成千上万倍的分子动力学（MD）模拟。通用机器学习势函数（uMLIP），如CHGNet、M3GNet，更是“博学多才”，在庞大的材料数据库上进行了预训练，理论上能泛化到多种化学体系。但问题是，当面对像卤化物电解质这样具有特殊局部结构和无序性的新材料家族时，这些“通才”模型还能否保持高精度？答案往往是否定的，直接应用常会出现能量预测偏差、高温模拟不稳定等问题。

本文要分享的，正是我们针对这一挑战所构建的一套完整解决方案：如何通过一套系统、高效的微调（Fine-tuning）流程，将一个通用的CHGNet势函数，精准地“调教”成专属于Li3YCl6-xBrx卤化物电解质家族的“专家”模型。我们不仅解决了uMLIP在特定体系应用时的精度和稳定性难题，更关键的是，建立了一个从“无序实验结构”到“可靠离子电导率预测”的标准化工作流。这个工作流的核心价值在于，它用极低的数据成本（仅需数百个DFT计算），就将模拟的预测精度提升到了接近从头算的水平，同时将计算速度提升了四个数量级，使得在普通计算集群上开展纳秒级的、可预测相变和离子传输的模拟成为现实。无论你是计算材料学的研究者，还是对固态电池材料设计感兴趣的工程师，这套方法都能为你提供一条从原子模拟到性能预测的清晰、高效的路径。

2. 核心思路拆解：从无序到有序，从通用到专用

我们的目标很明确：对Li3YCl6-xBrx（LYCB）这一混合卤化物电解质家族，实现其晶体结构、相稳定性和锂离子电导率的高精度、高效率预测。面临的初始挑战也是材料模拟中的经典难题：实验给出的晶体结构往往是“无序”的。以Li3YCl6和Li3YBr6为例，其精修结构中的锂（Li）和钇（Y）位点存在部分占位，这意味着在模拟中，我们需要从海量可能的原子排列方式中，找到一个能量最低的、有物理意义的“有序”模型作为起点。

2.1 结构枚举与能量排序：为模拟奠定物理基础

直接对包含部分占位的超大超胞进行DFT能量计算是不现实的，组合数会随着超胞尺寸指数增长。我们的策略是“分而治之，高效筛选”：

1. 构建超胞与对称性枚举：首先，我们基于实验结构构建了1×1×2的超胞。对于Li3YCl6（三方相，P-3m1）和Li3YBr6（单斜相，C2/m），这分别得到了包含60和40个原子的模型。接着，利用enumlib这类对称性枚举工具，系统地生成所有对称性不等价的原子有序排列方式。这一步产生了约6万个LYC和20万个LYB的候选结构。

2. 预筛选与精修：如果直接用DFT优化这数十万个结构，计算量依然惊人。因此，我们引入了一个“快速裁判”——预训练的M3GNet通用势函数。用它来快速优化所有枚举结构并计算其总能量。虽然M3GNet的绝对精度可能不如DFT，但其相对能量排序通常是可靠的。我们从M3GNet排序的结果中，挑选出能量最低的20个结构。

3. DFT验证与基态确定：对这20个最有希望的候选结构，我们动用“金标准”DFT进行完整的结构优化和能量计算。这一步确保了最终确定的基态有序结构（我们称之为LYC0和LYB0）的可靠性。LYC0的空间群为P-3c1，LYB0为C2，它们为后续的所有模拟提供了物理上合理的初始构型。

实操心得：这一步的关键在于信任uMLIP在相对能量预测上的能力，而非绝对能量。用M3GNet做初筛，再用DFT做最终裁决，是一种在计算效率和精度之间取得完美平衡的策略。对于混合卤化物LYCB，我们通过对LYC0和LYB0这两个端点结构进行卤素（Cl/Br）替换，再结合类似的枚举排序流程，就能高效地得到整个成分范围内（x=0到6）的有序结构模型。

2.2 迭代式微调：让通用模型成为领域专家

获得了可靠的初始结构后，下一个挑战是：预训练的CHGNet模型在LYCB体系上，特别是在有限温度下，表现不佳。如图2所示，它在高温（T≥600 K）下能量预测误差增大，体积预测严重高估，甚至在500 K以上的NpT（恒温恒压）模拟中会失稳。这说明预训练模型未能充分捕捉该体系在有限温度下所经历的结构和热力学多样性。

我们的解决方案是一个迭代式、主动学习启发的微调流程（图1蓝色部分），其核心思想是：让模型在它即将出错的“边界”附近学习。

1. 启动模拟与数据采集：从LYC0和LYB0结构出发，使用当前版本的CHGNet模型（初始为预训练模型）在某个起始温度（如T₀=200 K）下进行NpT-MD模拟，使体系在目标温度和压力（1 bar）下充分平衡。随后进行NVT（恒温恒容）生产模拟，产生包含数十万帧原子构型的轨迹。

2. 代表性结构采样：存储每一帧构型既不必要，计算成本也高。我们采用最远点采样（Furthest Point Sampling）算法，基于原子位置的均方根误差，从整个MD轨迹中挑选出约250个在几何结构上最具多样性的“代表性”构型。这确保了训练数据能最大程度地覆盖相空间。

3. DFT标注：对这250个左右的结构，执行静态DFT计算，精确得到每个构型的总能、原子受力和应力张量。这些数据构成了本轮微调的“黄金标准”训练集。

4. 模型微调：以当前的CHGNet模型为起点，用新生成的（LYC+LYB）DFT数据对其进行微调。这里有一个关键技巧：我们并非冻结所有层，而是有选择地解冻部分网络层进行训练，同时优化批次大小、学习率等超参数。这样既能让模型快速适应新体系，又避免了完全重新训练导致遗忘原有知识（灾难性遗忘）。

5. 温度爬升与迭代：将微调后得到的新模型（如CHGNet_200K）用于下一个更高温度（T₁=400 K）的MD模拟，重复步骤1-4。如此循环，逐步将训练温度提升至800K。通过这种方式，模型的学习数据从低温平衡结构，逐步扩展到包含高温热涨落、甚至可能接近相变的复杂构型，其预测能力和稳定性也随之层层递进。

注意事项：温度增量ΔT的选择需要权衡。ΔT太小（如50K），迭代次数过多，计算成本增加；ΔT太大（如300K），相邻温度区间的结构差异可能过大，导致在新温度下模拟一开始就失稳，无法采集到有效数据。我们选择200K作为一个平衡点，既能保证平稳过渡，又控制了迭代成本。

3. 微调效果深度评估：精度、稳定性与效率的飞跃

经过从200K到800K的迭代微调，我们得到了一系列模型（CHGNet_200K, _400K, _600K, _800K）。那么，这些“专家”模型究竟比原来的“通才”强在哪里？

3.1 静态与热力学性质的显著提升

我们首先在一个独立的测试集（包含从0K到800K的LYCB结构）上评估了模型的预测精度（图2b, c）。

• 能量、力与应力：预训练CHGNet对LYCB体系的总能预测平均绝对误差（MAE）在~11 meV/atom。而经过微调后，最佳模型CHGNet_600K将能量MAE降低至3.5 meV/atom以下，力的MAE低于50 meV/Å，应力的MAE低于85 MPa。这个精度水平已经与当前最先进的专用MLIP相当，意味着在原子相互作用细节的刻画上，微调后的模型已经达到了基于DFT数据的MLIP所能期望的精度天花板。
• 体积预测与热稳定性：这是微调带来的最直观改进之一。如图2d,e所示，预训练CHGNet在高温下严重高估晶胞体积（误差达25-40%），且模拟易失稳。而CHGNet_600K的预测与另一种先进势函数SevenNet以及直接DFT计算的结果高度吻合，偏差在2%以内。更重要的是，CHGNet_600K能够稳定地进行高达800K的长时间NpT-MD模拟，这为研究材料的热膨胀、高温相变等性质奠定了基础。

值得注意的是，即使微调只使用了端点成分LYC和LYB的数据，模型对中间成分LYCB的预测依然保持高精度，这证明了该方法良好的跨成分迁移能力。

3.2 与实验的对比及“体积效应”的深入探讨

将模拟体积与实验测量值对比，我们发现即使经过微调，预测体积仍比实验值大5-10%。这并非模型缺陷，而是需要理解的物理：

1. DFT泛函的系统误差：我们训练数据所用的PBE泛函本身就有轻微高估晶格常数的倾向。模型学到了DFT的“世界观”，其预测自然与DFT趋势一致。
2. 理想晶体 vs. 真实材料：模拟用的是完美、无缺陷的晶体超胞。而真实样品存在晶界、点缺陷、微应变等，这些都会导致实验测得的平均体积偏小。
3. 压力调控验证：为了验证体积差异对关键性质（离子电导率）的影响，我们做了一个巧妙的对照实验：在10 kbar的外压下进行NpT平衡，获得一个体积更接近实验值的压缩晶胞。后续的离子扩散模拟将分别在1 bar（DFT精度体积）和10 kbar（近实验体积）平衡的两种晶胞上进行。这帮助我们区分了“本征离子输运能力”和“晶格体积效应”。

4. 锂离子输运性质的高通量预测

有了稳定可靠的CHGNet_600K势函数和经过NpT平衡的初始结构，我们就可以大规模、低成本地探究LYCB家族的离子输运行为。我们在400K到750K的温度范围内，对多种成分进行了2纳秒的NVT-MD模拟，并从中提取了锂离子的自扩散系数（D）、活化能（Eₐ）和通过Nernst-Einstein关系换算得到的室温离子电导率（σ₃₀₀K）。

4.1 LYC与LYB：各向异性与各向同性的鲜明对比

模拟结果揭示了LYC和LYB截然不同的锂离子迁移机制（图4）：

• Li₃YCl₆ (LYC)：扩散具有强烈的各向异性。在500K时，锂离子沿c轴方向的扩散系数比在ab面内高2-3倍。从锂离子概率密度图可以清晰看到，锂离子主要沿着c轴方向，通过共享面的八面体-四面体-八面体路径进行快速迁移，而在ab面内通过共享边的路径迁移则困难得多。这与其层状晶体结构提供的沿c轴开放通道是吻合的。
• Li₃YBr₆ (LYB)：扩散几乎是各向同性的。锂离子在所有方向上的迁移能力相近。概率密度图显示，迁移主要通过共享边的八面体-四面体-八面体路径进行，没有像LYC那样明显的优先通道。

这种差异直接影响了它们的宏观离子电导率。对于体积压缩到近实验值的结构（10 kbar平衡），CHGNet_600K预测LYC的室温电导率约为15.7 mS cm⁻¹，活化能0.24 eV；而LYB的室温电导率约为0.58 mS cm⁻¹，活化能0.36 eV。LYC更高的电导率和更低的活化能，与其提供的快速一维通道密切相关。

4.2 混合卤化物LYCB：成分调控与性能优化

我们进一步将研究扩展到整个LYCB固溶体系列（x=3到6，即富溴端）。如图5所示，随着Cl含量增加（x减小），室温离子电导率呈现先相对平稳、后单调上升的趋势（在x=3到5之间），同时活化能相应降低。

• 机理分析：这种性能提升可归因于Cl部分取代Br后引入的晶格无序和应变。较小的Cl⁻离子取代较大的Br⁻离子，会导致局部晶格畸变，可能拓宽了锂离子的迁移通道，或者创造了更多低能垒的跃迁路径，从而降低了整体扩散活化能。
• 与实验的对比与启示：我们的预测显示，通过Br→Cl取代，可以将LYB的离子电导率提升一个数量级。虽然模拟没有完全复现实验中在x≈4.5处出现的尖锐电导率峰值，但趋势一致。这暗示实验中的峰值可能还与更复杂的因素有关，如有序-无序相变、缺陷浓度变化等，这些是未来模拟可以深入的方向。

实操心得：电导率计算的误差处理。从MD轨迹计算扩散系数D时，需要拟合均方位移（MSD）曲线。一个关键细节是拟合区间的选择。我们排除了短时标的弹道扩散区域（< 几ps）和长时标因统计性变差而波动大的区域（> 总模拟时间的70%）。对于2 ns的轨迹，我们将最大拟合窗口限制在500 ps，以确保线性拟合的可靠性。误差棒则通过将轨迹分成多个区块进行统计得到，这比单纯看拟合的R²值更能反映结果的统计不确定性。

5. 完整工作流再现与关键实操指南

如果你想在自己的研究体系中复现或借鉴这套方法，以下是基于我们实践总结的核心步骤和工具链：

5.1 软件环境与数据准备

• 核心软件包：
  • DFT计算：VASP。需要准备好相应的赝势和INCAR参数（我们使用PBE泛函，平面波截断能520 eV，k点网格保证能量收敛至0.5 meV/atom以内）。
  • 结构处理与枚举：Pymatgen, ASE。用于构建超胞、结构操作和IO。
  • 对称性枚举：enumlib。集成在ASE中，用于生成有序构型。
  • MLIP与快速筛选：matgl(包含M3GNet模型)。用于初始结构的快速能量评估和排序。
  • MLIP微调与推理：chgnet(PyTorch后端)。用于加载预训练模型、进行微调训练和作为ASE的Calculator进行单点计算或结构优化。
  • 分子动力学模拟：LAMMPS（需编译支持CHGNet等MLIP的版本）。我们使用GPU加速的版本进行大规模MD。
• 初始数据：你需要准备目标体系的实验晶体结构文件（CIF格式），明确其中的部分占位原子。

5.2 分步实操流程

1. 结构建模与基态搜索：

   # 示例：使用Pymatgen和ASE生成有序构型并用M3GNet排序 from pymatgen.core import Structure from pymatgen.analysis.ewald import EwaldSummation from ase import Atoms from matgl import load_model import numpy as np # 1. 加载实验结构，构建超胞 exp_structure = Structure.from_file("Li3YCl6_disordered.cif") supercell = exp_structure * (1, 1, 2) # 构建1x1x2超胞 # 2. 使用enumlib生成有序排列 (此处为概念代码，实际调用ASE的enum函数) # from ase.build import make_supercell # 调用枚举算法，生成所有对称性不等价的Li/Y有序排列列表: ordered_structures_list # 3. 加载M3GNet模型 potential = load_model("M3GNet-MP-2021.2.8-PES") # 4. 循环优化并计算能量 energies = [] for idx, struct in enumerate(ordered_structures_list): ase_atoms = AseAtomsAdaptor.get_atoms(struct) ase_atoms.calc = potential opt = BFGS(ase_atoms) opt.run(fmax=0.01) # 优化到受力小于0.01 eV/A energy = ase_atoms.get_potential_energy() energies.append(energy) # 可选：每100个结构保存一次进度 # 5. 排序并选取能量最低的20个 sorted_indices = np.argsort(energies) top_20_structures = [ordered_structures_list[i] for i in sorted_indices[:20]]

   将这20个结构用DFT进行精确优化，确定能量最低的作为后续模拟的基态有序模型。

2. 迭代微调流程实现： 这是一个循环过程，自动化脚本是关键。核心循环如下：

   # 伪代码展示迭代逻辑 current_model = load_pretrained_chgnet() current_temperature = 200 # K max_temperature = 800 # K temperature_step = 200 # K while current_temperature <= max_temperature: # 步骤1: 用当前模型进行NpT-MD，平衡体积 # 调用LAMMPS，输入：当前模型，基态结构， current_temperature, 1 bar # 产出：平衡后的结构和NVT轨迹文件 # 步骤2: 从轨迹中采样代表性结构 (例如使用最远点采样) sampled_structures = farthest_point_sampling(trajectory_file, n_structures=250) # 步骤3: 对采样结构进行DFT单点计算，获取E, F, Stress标签 dft_data = run_vasp_calculations(sampled_structures) # 步骤4: 微调CHGNet模型 # 加载当前模型，解冻最后几层，用dft_data训练若干epoch tuned_model = fine_tune_chgnet(current_model, dft_data, frozen_layers=10, learning_rate=1e-4, epochs=50) current_model = tuned_model save_model(current_model, f"CHGNet_{current_temperature}K.pth") # 步骤5: 升温，进入下一轮 current_temperature += temperature_step

3. 离子电导率计算： MD模拟完成后，使用自编脚本或工具（如pymatgen的DiffusionAnalyzer）分析轨迹：

   from pymatgen.analysis.diffusion.analyzer import DiffusionAnalyzer # 假设已经将LAMMPS轨迹处理成了pymatgen的Structure对象列表 `structures` 和时间列表 `times` analyzer = DiffusionAnalyzer.from_structures(structures, times, specie="Li", temperature=500) # 单位K diffusivity = analyzer.diffusivity # 单位 cm^2/s conductivity = analyzer.conductivity # 单位 mS/cm activation_energy = analyzer.activation_energy # 单位 eV

5.3 常见问题与排查技巧

• 问题1：微调过程中，模型在高温MD模拟中突然崩溃（能量/体积发散）。

  • 排查：这通常是因为训练数据未能覆盖当前高温下的某些高能构型。检查上一步采样出的训练结构是否足够“多样”。最远点采样的距离度量是否合理？可以尝试增加采样数量，或手动检查轨迹，看看崩溃前是否出现了异常的键长或键角。
  • 解决：降低温度增量ΔT，比如从200K改为100K，让模型更平缓地学习。或者在当前温度下，延长MD平衡时间，采集更充分的数据后再进行微调。

• 问题2：微调后的模型在能量预测上表现很好，但体积预测依然与实验有较大差距。

  • 排查：首先确认DFT计算所用的交换关联泛函。如果使用的是PBE，其本身就有约1%-2%的体积高估。其次，检查NpT模拟的压力设置是否正确（通常为1 bar），以及模拟时间是否足够长以达到充分平衡（我们通常平衡>100 ps）。
  • 解决：如果追求与实验体积一致，可以采用“压力校正”策略，即像我们一样，在较高的外压（如10 kbar）下进行NpT平衡，以获得与实验匹配的晶胞参数，再进行NVT扩散模拟。这能更真实地反映实际材料中的离子输运环境。

• 问题3：计算得到的离子电导率数量级与文献或实验差异很大。

  • 排查：
    1. 统计性：模拟时间是否足够长？锂离子的均方位移（MSD）是否达到了至少其离子间平均距离的平方（~ Ų）？通常需要模拟到MSD曲线进入清晰的线性扩散区。
    2. 体系尺寸：超胞是否足够大？避免因周期性边界条件导致的锂离子自相关。一般建议超胞边长至少是扩散粒子平均自由程的2倍以上。
    3. 体积效应：确认模拟所用的晶胞体积是来自NpT平衡后的，而不是0K优化的体积。体积对电导率影响极大。
    4. 公式应用：使用Nernst-Einstein关系计算电导率时，是否考虑了所有可移动离子（N）和电荷（q）？对于Li3YCl6，每个晶胞有3个可移动的Li⁺。
  • 解决：确保MD模拟时间足够（我们通常用2 ns），在分析时仔细选择MSD的线性拟合区间，并正确计算误差。将你的计算流程在一个已知结果的标准体系（如Li10GeP2S12）上测试一遍，可以快速定位问题。

6. 方法论的普适性与未来展望

我们建立的这套“实验结构枚举 → uMLIP快速筛选 → DFT精确定标 → 迭代式主动学习微调 → 高通量性质预测”工作流，其威力远不止于LYCB这一个体系。它为解决含有无序、部分占位或复杂缺陷的材料体系的原子模拟提供了一个通用模板。

• 普适性：任何从实验上已知存在无序，且其有序化模型对性能有关键影响的材料，如超离子导体、无序岩盐正极材料、高熵合金等，都可以套用此流程。核心是将难以处理的“无序”问题，转化为对有限个“有序”候选模型的筛选和优化问题。
• 效率与精度的平衡：该方法最突出的优势是数据效率。我们仅用了约250个/温度的DFT计算（主要花在LYC和LYB两个端元上），就使一个通用模型在复杂卤化物家族上达到了专用模型的精度。相比从头训练一个MLIP需要成千上万个数据点，这种方法节省了1-2个数量级的DFT计算成本。
• 对未来研究的启示：
  1. 主动学习的自动化：目前的迭代温度爬升是一种简单的主动学习策略。未来可以集成更智能的采样策略，例如基于模型预测不确定性（Uncertainty）来选择对模型改进最有价值的构型进行DFT计算，从而进一步提升数据效率。
  2. 多目标性质预测：当前微调主要针对能量、力和应力。可以进一步引入带电荷的MLIP（如CHGNet本身支持），或同时用DFT数据训练电子态密度、能带等电子性质预测头，实现一个模型预测多种性质。
  3. 面向更复杂的动力学过程：当前工作流能出色地处理平衡态扩散。未来可探索将其应用于非平衡过程，如电极/电解质界面的反应、枝晶生长初期过程模拟等，这可能需要采集包含化学反应断键/成键的构型数据。

回过头看，这项工作的最大价值在于它架起了一座桥梁：一边是追求高精度但计算昂贵的量子力学方法，另一边是追求高效率但精度有限的经典力场或早期MLIP。通过有针对性的、数据高效的微调，我们让强大的通用基础模型在特定的专业领域内发挥了极致性能。这不仅仅是解决了一个具体的材料科学问题，更是为计算驱动的新材料设计与发现，提供了一套可复制、可推广的方法论工具。在固态电��乃至更广阔的电化学储能材料探索中，这类能够深度融合先验知识、高效利用计算资源、并产出可靠微观机理的研究范式，将会变得越来越重要。