机器学习势函数驱动的高压材料全局结构搜索：以铋基化合物为例

1. 项目概述：当机器学习势函数遇上高压材料探索

在材料科学的前沿，尤其是在高压物理和新型功能材料设计领域，我们面临着一个核心挑战：如何在浩瀚的原子构型空间中，高效且准确地找到能量最低、最稳定的晶体结构？这个问题，就是晶体结构预测。传统上，我们依赖第一性原理计算，比如基于密度泛函理论的方法，来精确计算给定原子排列下的能量。但问题在于，对于一个稍微复杂点的二元或三元体系，可能的原子排列方式几乎是天文数字。用DFT去逐一计算，其计算成本是任何超级计算机都难以承受的。这就好比要在整个太平洋里找一颗特定的珍珠，用最精密的仪器（DFT）去检测每一滴水，理论上可行，但现实中不可能。

这就是机器学习势函数大显身手的地方。简单来说，我们可以把它理解为一个“超级速算员”。我们先让DFT这位“精密测量师”计算出成千上万个不同原子构型的能量和受力，然后用这些数据去训练一个机器学习模型。一旦训练完成，这个MLP模型就能在瞬间（通常是DFT计算时间的百万分之一甚至更短）给出一个新构型的能量预测，而且精度可以非常接近DFT。这样一来，我们就能用这个“速算员”去快速筛选太平洋里的海水，迅速定位到那些可能藏有珍珠的区域，然后再请“精密测量师”DFT去做最终确认。

我这次分享的工作，正是将这套“MLP速算员+DFT精算师”的组合拳，应用到了一个非常具体且有趣的领域：高压下的铋基二元化合物体系。为什么是铋？铋本身是一种具有独特电子结构和物理性质的后过渡金属，其高压相变非常丰富。而它与碱金属、碱土金属、稀土金属形成的化合物，常常展现出超导、拓扑绝缘体、热电等诱人的性质。在高压下，这些性质可能被进一步调控甚至激发出全新的物态。然而，这些体系在高压下的相图，尤其是那些尚未被实验发现的亚稳相，很大程度上仍是未知领域。

我们的目标很明确：利用基于机器学习势函数的随机结构搜索方法，系统性地对十个铋基二元体系（K-Bi, Mg-Bi, Ca-Bi, Sr-Bi, Ba-Bi, Sc-Bi, Y-Bi, La-Bi, Eu-Bi, Gd-Bi）在0到20 GPa的压力范围内，进行地毯式的全局结构搜索。我们不仅想验证已知的实验结构，更希望能像探矿一样，挖掘出大量隐藏在相图深处、尚未被报道的热力学稳定化合物，为后续的高压实验合成和物性研究绘制一份详尽的“藏宝图”。这项工作对于从事计算材料学、高压物理、以及新型功能材料设计的同行来说，是一次方法论上的完整实践，也提供了丰富的数据参考。

2. 核心方法解析：MLP-RSS工作流与关键技术细节

要把这个“大海捞针”的任务系统化，需要一个严谨、高效且可复现的工作流程。我们采用的核心方法是“机器学习势函数驱动的随机结构搜索”，其完整流程可以拆解为几个关键环节，每个环节都有需要注意的“坑”和技巧。

2.1 第一性原理数据集构建：质量的基石

一切始于高质量的训练数据。MLP的性能上限，很大程度上由训练数据的质量和广度决定。我们的策略是针对每一个A-Bi二元体系，分别训练一个专用的MLP。

数据生成策略：我们并非漫无目的地生成结构。初始数据集来源于几个方面：

1. 已知晶体结构数据库：从ICSD等数据库中提取该二元体系及纯组元A和Bi在常压和高压下已知的晶体结构。
2. 主动学习迭代：这是核心。我们先用一个小的初始数据集训练一个“粗糙”的MLP，然后用它驱动随机结构搜索。在搜索过程中，会生成大量新的、能量较低的候选结构。我们从中选取一部分代表性结构（例如，能量接近凸包、或结构多样性高的），用DFT进行精确的单点能或弛豫计算，将这些新的（输入结构，DFT能量/受力）对加入训练集。
3. 扰动与畸变：对已知稳定结构进行小幅度的晶格畸变、原子位置扰动，以增强MLP对势能面局部变化的敏感性。

实操心得：数据平衡的艺术训练数据不仅要“广”，还要“平衡”。如果数据集中全是能量很低的结构，MLP对高能区域的预测会非常不准，导致在搜索中无法正确翻越能量势垒。我们通常会确保数据覆盖从低能到较高能（比如高于凸包50-100 meV/atom）的广泛区域。同时，各化学成分比例的数据量也要大致均衡，避免模型偏向某个组分。

DFT计算参数：所有DFT计算我们使用VASP软件，采用PAW赝势和PBE泛函。平面波截断能设置为520 eV，k点网格密度保证倒易空间采样间距约0.03 Å⁻¹。对于含Eu和Gd的体系，由于4f电子的强关联效应，我们采用了DFT+U方法（U值根据文献设置为7.0 eV）。这些参数是在精度和计算成本间的折衷，经验证足以提供可靠的训练数据。

2.2 机器学习势函数模型选择与训练

我们选择了多项式机器学习势函数。与现在流行的神经网络势函数相比，多项式势函数虽然在表达极端复杂势能面时能力稍弱，但其具有物理可解释性强、训练稳定、推理速度极快（对超大规模搜索至关重要）的优点。对于金属和金属间化合物体系，其势能面相对“平滑”，多项式模型通常能取得很好的效果。

描述符是关键：我们使用原子中心对称函数作为描述符，它将每个原子周围的局部化学环境转化为一组数学向量。对于二元体系，我们考虑了A-A, A-Bi, Bi-Bi三种原子对的径向对称函数，以及三种组合的角向对称函数。截断半径通常设为5-6 Å，以包含足够的近邻原子信息。

训练与验证：数据集按8:1:1的比例随机划分为训练集、验证集和测试集。训练目标是最小化模型预测的能量、受力与DFT值之间的损失函数。我们密切监控验证集上的误差随训练周期的变化，一旦验证误差开始上升（过拟合），就提前停止训练。最终，我们要求测试集上每个原子的能量均方根误差低于2-3 meV，受力的RMSE低于0.1 eV/Å。这个精度足以可靠地区分能量差在5-10 meV/atom以上的结构。

2.3 随机结构搜索与凸包构建

拥有训练好的MLP后，就可以进行大规模的RSS了。对于每个A-Bi体系，在每个目标压力点（0, 5, 10, 15, 20 GPa），我们进行以下操作：

1. 结构生成：使用随机晶格向量和随机原子位置，生成数千至上万个初始结构。每个结构的原胞原子数（Z）覆盖一个范围（例如4到32个原子），以探索不同化学计量比。
2. MLP预筛选与弛豫：用训练好的MLP快速计算每个初始结构的能量，并基于MLP进行结构弛豫（最小化原子受力）。由于MLP评估一次仅需毫秒级，这步可以快速淘汰掉大量明显不合理的极高能结构。
3. 聚类与去重：对弛豫后的结构，通过比较其晶格参数、原子位置和指纹函数（如SOAP）进行聚类，去除重复或极其相似的结构，保留能量最低的几百个独特候选结构。
4. DFT精修与凸包计算：将MLP筛选出的低能候选结构，用DFT进行严格的几何优化和能量计算。最后，利用这些DFT精修后的能量，绘制该压力下该体系的形成焓-成分凸包图。落在凸包上的点，即对应热力学稳定的化合物；凸包上方的点，则为亚稳相。

注意事项：压力处理在MLP训练和RSS中��压力是一个外部变量。我们的做法是将压力作为描述符的一个额外输入特征，或者更常见的是，在每次弛豫时固定目标压强（通过设置VASP的PSTRESS参数），让MLP/DFT在恒压条件下优化晶胞体积。这样训练出的MLP或进行的DFT计算，直接对应特定压力下的能量。

3. 结果深度解读：十个铋基二元体系的高压相图宝藏

经过上述系统性的搜索，我们得到了十个体系在0-20 GPa压力下丰富的成果。这里我挑几个有代表性的体系，深入聊聊我们发现的有趣结构和背后的物理。

3.1 碱土金属体系（Mg-Bi, Ca-Bi, Sr-Bi, Ba-Bi）：从简单到复杂

Mg-Bi体系相对简单，但仍有新发现。在常压下，已知的稳定相是Mg3Bi2（反La2O3型结构）。我们的搜索不仅确认了这一点，还在高压下预测了一系列新的富Mg或富Bi的稳定化合物，例如Mg5Bi、Mg7Bi3、MgBi2等，它们具有诸如P3m1、I4/mmm、C2/m等空间群。这些结构大多由Mg和Bi原子层交替堆叠或形成特定的多面体配位构成。一个有趣的趋势是：随着压力增加，高对称性的结构（如立方相）往往更倾向于向低对称性（单斜、正交）转变，因为压力可以稳定更复杂的配位多面体。

Ca-Bi, Sr-Bi, Ba-Bi体系则展现出丰富的结构化学。这些碱土金属原子半径较大，与Bi形成的化合物结构类型更多样。我们预测了大量实验未报道的稳定相。例如，在Ca-Bi体系中，除了已知的Ca3Bi2，我们还发现了Ca11Bi3 (Cmmm)、Ca5Bi3 (Mn5Si3型)、Ca4Bi3 (Th3P4型) 等一系列新相。特别值得注意的是，许多这些预测的稳定结构，都属于典型的“Zintl相”或极性金属间化合物结构原型，如CuAu型、Th3P4型、Mn5Si3型等。这说明在高压下，电负性差异带来的电荷转移效应依然显著，并驱动形成了这些具有特定价电子规则的相。

实操心得：如何判断一个“新结构”是否合理？当RSS吐出一堆空间群和晶体结构时，不能只看形成焓。我们有一套“合理性检查清单”：

1. 几何优化稳定性：DFT弛豫后，所有原子受力是否真的接近零（<0.01 eV/Å）？晶胞应力是否平衡？
2. 动力学稳定性：计算声子谱，确认在整个布里渊区没有虚频（软模）。有虚频意味着结构在某个振动模式下不稳定，可能会畸变成其他结构。
3. 热力学稳定性：不仅对纯组元稳定，还要对分解成其他预测的二元相稳定（即落在凸包上）。
4. 化学直观性：检查键长、配位数是否合理。例如，Bi-Bi键长不应过短，金属-Bi键长应在典型范围内。不合理的配位（如Bi的配位数异常低）可能预示着计算问题。

3.2 稀土金属体系（Sc-Bi, Y-Bi, La-Bi, Eu-Bi, Gd-Bi）：f电子的挑战与机遇

稀土体系因其f电子的强关联和可能存在的磁性而更具挑战性，但也因此可能蕴含更奇特的物性。

Sc-Bi与Y-Bi体系：Sc和Y没有f电子，相对“干净”。我们预测了众多新相，如Sc6Bi、Sc5Bi3、ScBi2等。一个有趣的发现是，在Sc-Bi中，除了常见的NiAs型ScBi，我们还预测了R3m和FeSi型的ScBi相，它们能量非常接近，可能是竞争相。在Y-Bi体系中，YBi以NaCl型结构稳定存在，这与许多稀土-主族化合物一致。

La-Bi体系：La的4f轨道为空，但其体积较大。我们预测的稳定相包括La3Bi (Cu3Au型)、La2Bi (多种变体)、LaBi (CuAu型)等。这些结构中，La与Bi之间显示出强烈的离子-共价混合键合特征。

Eu-Bi与Gd-Bi体系：这是本次研究的难点和亮点。Eu和Gd具有半满和近半满的4f电子壳层，磁性效应显著。为了平衡计算成本与精度，我们在主搜索中采用了将4f电子处理为核电子（冻结在芯态）的赝势。令人惊喜的是，即使在这个近似下，我们的MLP-RSS方法仍然成功地预测出了许多已知的含Eu/Gd的铋化物结构，例如Ho11Ge10型的Eu11Bi10、Th3P4型的Eu4Bi3/Gd4Bi3、Mn5Si3型的Eu5Bi3/Gd5Bi3。这说明对于这些化合物的结构稳定性，几何堆积和静电相互作用可能起着主导作用，而f电子的精细关联效应对其晶体结构的影响是第二位的。

为了验证这一点，我们专门对Eu3Bi、EuBi、EuBi3、Gd3Bi、GdBi等几个关键相，用包含了4f价电子的赝势并开启DFT+U进行了额外的电子结构计算。如图21和22所示，对于富Bi的化合物（如EuBi3、EuBi、GdBi），两种赝势计算得到的费米能级附近的态密度（DOS）非常接近，晶格参数变化也很小（<1%）。然而，对于富稀土相的Eu3Bi，DOS在Eu的p、d轨道投影上显示出明显差异，晶格参数变化更大（从5.16 Å到5.31 Å）。这证实了我们的判断：在富稀土相中，f电子的离域性和关联效应确实对电子结构和晶格常数有更明显的影响，但对整体结构类型的预测影响相对有限。

3.3 纯组元高压相：方法的试金石

在搜索二元化合物之前，我们对所有纯金属（Na, K, Mg, Ca, Sr, Ba, Sc, Y, La, Eu, Gd, Bi）在0-20 GPa下的稳定相也进行了预测，这相当于对我们MLP在单质体系精度的检验。

结果总体令人满意。例如，我们成功预测了Na和K从BCC到FCC的相变，Ca和Sr从FCC到BCC的相变，Ba从BCC到HCP的相变，以及Sc、Y、La等稀土金属的一系列复杂相变序列（如HCP -> α-Sm型 -> α-La型）。这些预测与已知的高压实验数据基本吻合，证明了我们训练的MLP在描述元素金属高压行为方面的可靠性。

关于Bi的一个小插曲：对于纯Bi，已知其在高压下有一系列复杂相变（Bi-I, Bi-II, Bi-III, Bi-V）。我们的MLP搜索到了Bi-III（宿主-客体结构）和Bi-V（BCC）相，但未能直接找到最稳定的Bi-I（A7）相，而是找到了一个能量极其接近的Imma空间群结构。这与之前专门针对Bi训练的MLP工作结果一致。这提醒我们，MLP的预测能力存在边界，对于某些具有特殊电子关联或非常局域化键合的结构，可能需要更专门的数据集或更复杂的模型架构。不过，由于Bi-I和Imma结构能量几乎简并，这个微小偏差对后续二元化合物稳定性的判断影响甚微。

4. 技术价值与实操启示：MLP-RSS的威力与边界

通过这十个体系的研究，我们可以更深刻地体会到MLP驱动的高通量全局结构搜索在材料发现中的巨大价值，同时也看清了其当前的局限性。

4.1 方法论的威力：效率与覆盖度的革命

1. 搜索空间的海量覆盖：传统的基于DFT的直接搜索，受限于计算成本，通常只能探索几百到几千个结构。而我们的MLP-RSS方法，在每个成分-压力点都能轻松生成并评估数万甚至数十万个结构。图11-20中展示的每个子图里密密麻麻的点（局部极小值结构），就是这种覆盖度的直观体现。许多能量略高于凸包但结构新颖的亚稳相也被挖掘出来，为研究动力学形成路径提供了线索。
2. 凸包预测的高精度：尽管MLP存在微小误差，但通过“MLP粗筛+DFT精修”的两步策略，最终基于DFT能量绘制的凸包（图11-20中的(b)子图）非常清晰可靠。MLP预测的凸包（(a)子图）与DFT结果高度一致，这说明MLP出色地学习了势能面的整体拓扑，能够准确识别出真正的低能区域。
3. 新材料的“预测-验证”闭环：本研究预测了上百个实验未报道的、在热力学上稳定的铋化物新相。这些预测结果本身，就是给实验学家的一份“采购清单”。例如，预测的某些化合物具有已知超导或拓扑材料的原型结构（如NaCl型、CuAu型），这强烈暗示它们可能具备类似的优异性质，值得优先开展高压合成实验进行验证。

4.2 实操中的挑战与应对策略

尽管成果丰硕，但在实际运行这样一套流程时，会遇到不少挑战：

1. 数据生成与MLP训练的循环成本：主动学习迭代虽然高效，但每一轮都需要进行新的DFT计算，整个过程是串行的，总耗时依然可观。我们的策略是并行化与自动化：编写脚本自动管理结构生成、MLP训练、RSS运行、DFT任务提交和结果分析流程。同时，利用计算集群并行计算多个不同化学成分或压力的点。
2. “维度灾难”与搜索效率：即使有MLP，在超高维的构型空间中盲目随机搜索效率仍可能低下。我们结合了对称性约束（在生成随机结构时施加可能的空间群）和基于已知结构原型的衍生（对已知稳定结构进行元素替换、晶格畸变），作为RSS的补充，以更有针对性地探索潜在的低能区域。
3. MLP的“外推风险”：MLP在训练数据覆盖的区域内（内插）预测很准，但对于完全陌生的原子环境（外推）可能给出荒谬的结果。在RSS中，我们设置了一个能量阈值过滤器，如果MLP弛豫后的结构能量异常高，我们会用DFT快速计算一个单点能进行校验，如果DFT能量与MLP预测严重不符，则将该结构加入训练集，重新训练MLP。这是一个重要的在线学习环节。
4. 磁性、强关联与自旋轨道耦合：如Eu/Gd体系所示，对于磁性体系，标准的非自旋极化MLP会失效。我们的解决方案是使用非磁性或顺磁性的DFT数据训练MLP进行结构搜索，因为晶体结构主要由几何和静电因素决定。在找到候选稳定结构后，再使用包含磁性的DFT进行最终的能量排序和性质计算。对于自旋轨道耦合效应显著的重元素（如Bi本身），在最终DFT精修时需要考虑进去，但这通常不纳入MLP训练，因其对总能量的贡献相对较小且计算昂贵。

4.3 对未来工作的展望

基于这次的经验，我认为这个领域有几个明确的下一步方向：

1. 从二元走向多元：目前的成功证明了方法在二元体系的威力。自然的延伸是三元甚至四元体系。挑战在于组合爆炸，需要发展更高效的描述符和主动学习策略，以及可能考虑使用图神经网络等更能捕捉复杂化学环境的MLP模型。
2. 集成性质预测：不仅预测稳定结构，还能同步预测其基本电子性质（如带隙、态密度、声子谱）。可以训练“多任务”MLP模型，或在MLP弛豫后，用超快但近似的方法（如紧束缚模型、预训练的神经网络）进行性质初筛。
3. 动力学过程模拟：利用MLP进行纳秒甚至微秒级的分子动力学模拟，研究高压下的相变路径、成核动力学，而不仅仅是零温下的热力学稳定性。这对于理解实验上观测到的亚稳相形成至关重要。
4. 与实验的深度互动：将计算预测的X射线衍射图谱、拉曼光谱等与高压原位实验数据直接对比，实现“计算引导实验、实验验证计算”的实时互动模式，加速新材料的确证过程。

回过头看，这项工作的核心价值在于它不仅仅是一份关于铋基化合物的高压相图报告，更是一次对**“机器学习势函数+全局结构搜索”** 这一现代计算材料学研究范式的完整演示和压力测试。它证明了，通过精心设计的流程和对其局限性的清醒认识，我们确实有能力系统性地探索未知的材料空间，将新材料发现从“偶然邂逅”逐步推向“理性设计”的新阶段。对于想进入这个领域的研究者，我的建议是：从一个小而具体的二元体系开始，亲手走通“数据生成-MLP训练-RSS搜索-DFT验证-结果分析”的全流程，踩过里面所有的坑，你对这个强大工具的理解才会真正深入骨髓。