1. 项目概述：构建一个“聪明”的原子间势能模型

在材料模拟的世界里，我们一直面临着一个核心矛盾： 精度与效率的博弈 。第一性原理计算（如密度泛函理论，DFT）能给出近乎“标准答案”般精确的原子间相互作用，但计算成本高得吓人，通常只能处理几百个原子、几个皮秒的模拟。而传统的经验势函数（如EAM）虽然快，但精度和可迁移性（Transferability）往往堪忧，换个温度或压力条件可能就完全不准了。

机器学习原子间势能（MLIP）的出现，像是一道曙光。它的核心思想很直观：我们不从物理原理去推导势能公式，而是让机器学习模型（比如神经网络）从一个庞大的、由DFT计算生成的数据集中，“学习”原子构型与其总能量、原子受力之间的复杂映射关系。一旦模型训练好，它就能以接近DFT的精度，但快几个数量级的速度，预测新构型的性质。这让我们得以模拟数百万个原子在纳秒甚至微秒尺度上的行为，去研究相变、缺陷演化、力学性能等传统方法难以触及的问题。

然而，这个美好的愿景有个致命前提： 训练数据库的质量 。模型再强大，也是“巧妇难为无米之炊”。如果数据库不能全面、无偏地覆盖你感兴趣的热力学条件（如从室温到几千开尔文，从常压到几百吉帕），那么模型在这些未覆盖区域的外推预测就极不可靠，这就是所谓的“可迁移性”差。

传统构建数据库的方法，比如基于已知晶体结构进行应变、扰动，或者运行昂贵的从头算分子动力学（AIMD），都存在明显短板。前者严重依赖先验知识，可能错过未知的亚稳相；后者采样效率低，且在高能态（如接近熔化）时模拟不稳定，难以系统性地探索整个构型空间。

我这次要分享的，正是一个旨在从根本上解决这个问题的框架： 基于嵌套采样（Nested Sampling, NS）与原子簇展开（Atomic Cluster Expansion, ACE）的自主热力学信息数据库生成方法 。简单说，我们不再“猜”需要什么数据，而是让算法自动地、系统地去“探索”整个原子构型空间，并基于热力学权重来智能地收集训练样本。我们以金属镁（Mg）为例，成功构建了一个能在0-600 GPa压力和0-8000 K温度范围内都保持高精度的MLIP模型。下面，我就来拆解这个框架的每一个环节，分享其中的设计逻辑、实操细节以及我们踩过的坑。

2. 核心思路：为何是嵌套采样与ACE的联姻？

在深入细节前，理解“为什么是这两个技术”至关重要。这决定了整个框架的效率和鲁棒性。

2.1 嵌套采样：一个“公平”的构型空间探险家

嵌套采样本质上是一个用于高效计算配分函数和探索高维空间的统计物理算法。你可以把它想象成一个智能的“探险队”：

1. 初始化 ：在广阔的构型空间（从完全无序的气体到紧密的固体）中，随机放置K个“探险者”（称为walker），每个walker代表一个具体的原子构型。
2. 淘汰与迭代 ：计算每个walker的焓值（在恒定压力下，相当于系统的“成本”）。淘汰掉焓值最高的那个（即当前最“差”的构型）。
3. 繁殖与探索 ：从剩下的walker中随机选一个，复制它，然后让这个副本在构型空间中进行一段随机游走（通过蒙特卡洛或短分子动力学），试图找到一个焓值更低的“新地点”。
4. 接受新成员 ：如果新地点的焓值低于其“父辈”，它就作为新walker加入队伍；否则就重试。
5. 循环 ：重复淘汰最高焓值walker的过程。

这个过程的精妙之处在于， 淘汰的速率直接关联于相空间体积的收缩 。算法会从高焓（无序、气体）区域开始，逐步“挤压”到低焓（有序、固体）区域。最终，它产生的构型样本，其分布 自然地正比于该构型在给定热力学条件下的统计权重 。这意味着：

• 无偏性 ：不需要预先知道镁有HCP、BCC、FCC等相，算法会自动发现它们。
• 热力学代表性 ：在高温熔融态，液相的样本会更多；在低温固态，基态结构的样本会更集中。这正好符合我们构建数据库的需求——样本密度应反映其热力学重要性。
• 全面性 ：能系统地覆盖从气体到固体所有可能的相，包括亚稳相和过渡态。

2.2 原子簇展开：一个高效且严谨的“学习者”

有了高质量数据，我们需要一个强大的模型来学习。原子簇展开是一种构建原子间势能的数学框架，它通过系统地对原子局部环境进行对称性不变的展开来描述势能面。

它的核心优势在于：

• 完备性 ：在理论上，当展开阶数足够高时，ACE可以精确表示任何满足平移、旋转和置换不变性的势能面。
• 高效性 ：其函数形式经过精心设计，能量、力的计算速度极快，且易于求导，非常适合大规模分子动力学模拟。
• 可解释性 ：展开系数与物理上有意义的原子簇关联，模型复杂度（由“阶数”和“度数”控制）可以系统性地增加，以平衡精度和计算成本。

为什么是ACE而不是神经网络（如Neural Network Potentials, NNP）？ 在项目中我们也权衡过。NNP非常强大，但训练更不稳定，需要更精细的超参调整，且推理速度通常慢于ACE。ACE提供了一个更“可控”和“稳定”的框架，对于构建一个需要反复迭代、自动化的主动学习流程来说，确定性更强，调试更简单。

两者的结合点 ：NS负责生成 热力学权重正确 的、 全面 的构型样本；ACE负责高效、准确地 拟合 这些样本数据。NS的探索过程需要调用数百万次势能评估，ACE的快速推理使得这种大规模探索在计算上成为可能。而ACE拟合出的势能模型，又可以作为下一轮NS探索的势函数，形成一个“探索-学习-再探索”的增强回路。

3. 实操流程拆解：从零构建镁的势能模型

我们的完整工作流是一个包含多个“周期”的主动学习循环，下图概括了核心步骤： （ 注：此处原图示意为：1. 嵌套采样生成构型 -> 2. 基于温度选择构型构建数据库 -> 3. DFT从头算评估 -> 4. 拟合ACE势能模型 -> 返回步骤1，使用新模型进行下一轮采样 ）

下面，我分步详解每个环节的实操要点。

3.1 周期0：用经典势函数启动探索

万事开头难。NS需要一个初始的势函数来评估构型的焓值，但我们还没有MLIP。这里采用了一个务实策略： 用一个已知但不完美的经典势函数（如EAM势）来启动 。

• 操作 ：我们选用了一个已有的镁EAM势。它在常压下表现尚可，但已知会错误预测某些高压相变（如将HCP到BCC的相变压力预测错误）。这反而成了我们的优势——一个不完美的起点，正好可以测试后续学习循环的修正能力。
• 采样设置 ：
  • 系统大小 ：使用16个原子的超胞。这是一个权衡，太小可能无法捕捉某些相的特征，太大则NS计算量激增。16原子对于镁的常见晶体结构（HCP, BCC, FCC）的晶胞是兼容的。
  • 压力点 ：在0-45 GPa范围内选择了11个压力点（0, 1, 5, 10, …, 45 GPa）。这个范围覆盖了镁在常温常压到中等压力下的主要相区。
  • NS参数 ：设置K=200个walkers。Walk length（随机游走步数）需要足够长以使walker充分解关联，我们经过测试设为1000步。
• 构型选择与数据库构建 ： NS��行结束后，会输出数十万个构型。我们不可能全部做DFT计算。这里利用了NS的一个副产品： 每个构型都有一个关联的“温度”估计值 （通过配分函数计算得出）。
  1. 我们设定一个温度范围（200-3000 K），这个范围应宽于目标材料在该压力区间的熔点。
  2. 在此温度范围内，按NS迭代步数等间隔选取100个构型。这样能保证在焓值-迭代步数图上均匀取样，从而覆盖从固态到液态的连续变化。
  3. 对11个压力点重复此操作，得到1100个16原子构型（共17,600个原子环境），构成初始数据库。
• DFT计算 ：使用VASP软件包对这些构型进行单点能、力和应力张量的计算。这是整个流程中最耗时的部分，但也是“黄金标准”数据源。
• 首次拟合 ：用这1100个构型的DFT数据，训练第一个ACE模型（记为C1 O4 D14，意为周期1使用的、阶数4、度数14的ACE模型）。

实操心得1：初始势函数的选择 不一定非要EAM。如果你有一个覆盖更广的“基础MLIP”（Foundation MLIP），用它启动会更好。甚至一个非常简单的势（如Lennard-Jones）也可以，只要它能定义一个基本的原子排斥和吸引。关键在于，这个初始势需要能引导NS探索到有物理意义的构型区域，避免一开始就陷入完全不现实的构型空间。

3.2 周期1与2：引入“委员会”不确定性，主动规避陷阱

用C1模型进行新一轮NS采样时，我们遇到了MLIP的典型难题： 势能面“空洞” 。

• 问题 ：MLIP在训练数据未覆盖的区域可能会产生非物理的预测，比如在原子距离极近时给出异常低的能量。NS的 exhaustive 探索特性极易发现这些“空洞”。一旦walker掉入空洞，由于其极低的能量，它会一直被复制，导致整个采样停滞在该非物理区域。
• 解决方案：委员会模型与不确定性量化 我们不是训练一个ACE模型，而是训练一个由5个相同架构但不同随机初始化的模型组成的“委员会”。对于同一个构型，5个模型会给出5个略有不同的能量预测。
  • 标准差（STD）作为不确定性度量 ：计算这5个预测值的标准差。在训练数据密集的区域，所有模型预测一致，STD小；在数据稀疏或未见的区域，模型“各执己见”，STD大。
  • 设置STD阈值 ：我们发现，当构型接近非物理的“空洞”时，其能量预测的STD会异常高（> 1 eV/atom）。而在正常的物理构型区域，即使是在高压高温下，STD通常也小于 0.1 eV/atom。
• 在NS中集成拒绝机制 ： 在NS的“随机游走”步骤（克隆并扰动walker）中，我们增加了一个接受准则： 新构型的委员会STD必须低于一个阈值（我们设为62.5 meV/atom） 。如果超过，则拒绝此次移动，重新尝试。同时，如果超过90%的walker都因高STD被“困住”，则终止本次NS运行。
• 效果 ：如图4所示，在没有STD限制时（左列），NS在高压力下采样到了大量原子间距极短（非物理）且STD极高的构型。启用限制后（右列），这些非物理构型被有效过滤，采样集中在合理的物理构型空间内。

实操心得2：不确定性阈值的设定 62.5 meV/atom这个阈值不是魔法数字，需要通过试验确定。我们的方法是：在已知可靠的构型（如平衡晶格）上计算委员会STD，得到一个基准值（通常很小）。然后，故意构造一些明显非物理的构型（如原子重叠），计算其STD。阈值应设在这两者之间，并留有足够的安全边际。这个机制是保证主动学习流程稳定运行的关键。

3.3 周期3与4：细化采样与加权策略

在周期2之后，模型已经能较好地描述主要相变。周期3的目标是 填充数据库 和 改进对亚稳相的描述 。

• 扩展压力范围 ：将NS采样压力点扩展到0-600 GPa，以覆盖镁的整个高压相图（HCP -> BCC -> FCC）。
• 调整采样策略 ：由于模型已趋完善，“空洞”问题基本解决，我们移除了STD限制，让NS自由探索。
• 样本加权 ：随着数据库增大，我们需要考虑样本权重。对所有样本一视同仁（等权重）可能会让高能、不常见的构型对损失函数产生过大影响。我们采用了基于能量的加权方案： 权重 ∝ exp(-α * (E - E_min)) 其中α是一个小的正数（如0.1），E是构型的每个原子能量，E_min是数据库中能量的最小值。这使得低能量（更稳定）的构型在训练中占有更高权重，符合物理直觉。
• 周期4：专注低焓区与模型增强 为了更精确地捕捉基态和低能激发态的性质（如声子谱、弹性常数），我们在周期4做了两件事：
  1. 减小系统尺寸 ：使用8原子超胞进行NS。更小的系统意味着在相同的计算资源下，我们可以对低焓值（低熵）的固态相进行更密集的采样。
  2. 增加模型复杂度 ：将ACE模型的度数从14提升到18（C5 O4 D18）。更高的度数意味着模型具有更强的拟合能力（更高的“灵活性”），可以刻画更精细的势能面细节，当然计算成本也会增加。

经过5个周期的迭代，我们最终的数据包含 8230个构型，总计100,480个原子环境 。相比传统方法动辄需要数万甚至数十万个DFT计算，这个数据量是相当高效的。

4. 模型性能验证：不仅仅是能量预测

一个MLIP的好坏，不能只看它在训练集上的能量、力拟合误差（虽然我们的模型这些误差都很低，~1 meV/atom）。更重要的是看它能否 正确预测从未在训练中明确出现过的、衍生的物理性质 。

4.1 零温焓-压曲线与相变

我们计算了HCP, BCC, FCC, dHCP（双六方密堆）四种晶体结构在0 K下随压力变化的焓值，并与DFT结果对比。

• 结果 ：如图6所示，即使是仅用初始EAM数据训练的C2模型，已经能正确预测HCP是常压下的稳定相，BCC在高压下更稳定。随着数据增加和模型增强，C5模型预测的HCP->BCC和BCC->FCC的相变压力与DFT结果高度吻合（误差在几个GPa以内）。
• 关键发现 ：亚稳相（如dHCP）的预测误差相对较大。这正是我们方法 热力学导向采样 的自然结果——NS倾向于在热力学权重高的区域（稳定相）采集更多样本。如果你需要极度精确地预测某个亚稳相，可能需要在采样中引入针对性的偏置或额外添加该相的结构。

4.2 声子谱与弹性常数

这些性质对势能面的二阶导数（力常数）非常敏感，是检验MLIP“局部”准确性的试金石。

• 声子谱 ：我们计算了HCP-Mg在0 GPa和50 GPa下的声子色散关系。我们的ACE模型完美复现了DFT的结果，包括所有声子支的频率和走向，没有出现任何虚频（动力学不稳定的标志）。
• 弹性常数 ：计算了零温零压下的弹性常数矩阵（C11, C12, C13, C33, C44）。与DFT和实验值对比，所有独立分量的误差均在5%以内。这证明了模型在预测材料力学响应方面的可靠性。

4.3 压力-温度相图

这是最终的“大考”。我们使用最终版的C5模型，结合NS（及其变体副本交换NS以提升相变分辨率）直接计算了镁在0-100 GPa， 0-4000 K范围内的P-T相图。

• 成���复现 ：模型成功预测了固-固（HCP-BCC）相变线、熔化线以及固-液-气三相点。
• 揭示争议区域 ：在5-20 GPa、高温固态区（熔化线以下），实验上曾观测到可能存在一个类似dHCP的未知相，但结构未定。我们的模拟在这一区域显示了复杂的自由能景观，提示可能存在丰富的相行为，为后续研究提供了线索。

5. 常见问题、挑战与应对策略

在实际操作中，我们遇到了不少典型问题，以下是总结和解决方案：

问题1：NS计算量巨大，如何优化？

• 挑战 ：一次NS运行需要10^8次量级的能量/力评估，即使使用快速的MLIP，也是一笔可观开销。
• 策略 ：
  • 并行化 ：NS算法天生适合并行。K个walker可以分布在多个计算节点上独立评估其能量/力。
  • 系统尺寸 ：在探索阶段，使用尽可能小的超胞（如16或32原子）。在验证最终性质时，再使用大体系。
  • 控制walker数量 ：K值越大，采样分辨率越高，但成本也越高。需要平衡。我们从K=200开始，后续周期可酌情调整。
  • 利用MLIP加速 ：这正是本框架的核心——用便宜的MLIP评估替代昂贵的DFT评估进行探索。

问题2：DFT计算失败（不收敛、SCF振荡）怎么办？

• 挑战 ：NS可能产生一些原子距离非常近或电子结构极其复杂的构型，导致DFT计算无法收敛。
• 策略 ：
  • 预处理过滤 ：在提交DFT计算前，检查构型中最小原子间距。如果小于一个安全阈值（例如，镁的金属半径的0.5倍），则直接丢弃该构型。这通常与高STD的“空洞”构型重合。
  • 灵活的DFT设置 ：对于少数难以收敛的构型，可以尝试放宽收敛标准、使用不同的混合泛函参数、或增加平面波截断能。但需谨慎，避免为个别“怪异”构型引入不一致的计算设置。
  • 委员会STD把关 ：如前所述，高STD是DFT可能失败的强预警信号。将其作为构型是否送入DFT计算的首要过滤器。

问题3：如何确定ACE模型的复杂度（阶数、度数）？

• 挑战 ：模型太简单（欠拟合）精度不够；太复杂（过拟合）计算慢且可能不稳定。
• 策略 ：
  • 从简开始 ：初始周期使用中等复杂度模型（如O4 D12）。
  • 监控验证集误差 ：始终划分一个独立的验证集（不参与训练），观察其能量、力误差随训练周期的变化。如果训练误差持续下降而验证误差上升，可能是过拟合。
  • 通过物理性质检验 ：最终判断标准是衍生性质（如晶格常数、弹性常数、声子谱）的准确性。如果这些性质与DFT有系统偏差，且增加数据无法改善，则可能需要提升模型复杂度（如我们从D14提升到D18）。
  • 计算成本权衡 ：更高阶/度的模型，其单次能量评估耗时呈多项式增长。需要在精度需求和模拟效率间取得平衡。

问题4：数据库样本权重如何设置？

• 挑战 ：不同构型的物理重要性不同，等权重训练可能不是最优。
• 策略 ：
  • 热力学权重 ：我们的方法中，NS产生的样本本身已近似具有正确的玻尔兹曼权重。但在构建离散数据库时，我们采用了基于能量的指数衰减加权（公式见3.3节），这近似于在有限温度下赋予低能构型更高重要性。
  • 力与应力的权重 ：在损失函数中，力和应力项的权重通常需要高于能量项。一个经验法则是，让能量、力、应力对总损失函数的贡献在同一数量级。我们通常设置力权重是能量权重的100-1000倍。
  • 手动调整 ：如果特别关注某个相或某种性质（如表面能），可以手动增加数据库中相关构型的权重。

构建一个鲁棒、可迁移的MLIP数据库，远不止是收集一堆DFT数据那么简单。它需要一套系统的、物理信息驱动的策略，去主动探索未知的构型空间，并智能地分配计算资源。基于嵌套采样和ACE的框架，通过将热力学采样与主动学习、不确定性量化相结合，为我们提供了一条通向这个目标的清晰路径。

镁的案例证明了这条路径的可行性。最终得到的模型，仅用约8000个构型的DFT数据，就在极端宽泛的P-T范围内展现了出色的预测能力。这套方法的美妙之处在于其通用性——它不依赖于对材料的先验知识，可以自动化运行，并且通过委员会不确定性机制保证了流程的稳定性。对于研究复杂相图、高压高温行为或多元合金体系，这种自动化的、热力学导向的数据库构建方法，无疑是一个强大而高效的工具。