企业官网建设流程全解析

1. 项目概述：当机器学习势函数遇上嵌套采样，如何绘制极端高压下的镁相图？

在计算材料科学的前沿，我们常常面临一个核心矛盾：精度与效率的权衡。第一性原理计算，如密度泛函理论（DFT），能提供近乎“金标准”的精度，但其巨大的计算成本将我们牢牢限制在数百个原子、皮秒时间尺度的模拟中。而当我们想探究材料在极端条件（例如高达600 GPa的超高压）下的相变行为、熔化曲线时，这无异于杯水车薪。传统经验势函数（如EAM）速度虽快，但其精度和可迁移性往往在远离训练集的区域急剧下降，预测结果可能与实验和第一性原理结果相去甚远。

近年来，机器学习势函数（Machine Learning Interatomic Potentials, MLIPs）的崛起，为打破这一僵局提供了利器。它本质上是一个“数据驱动的拟合器”，通过学习海量DFT计算产生的能量、力、应力数据，构建出一个既能保持量子力学精度、又能以经典力场速度进行计算的代理模型。其中，原子簇展开（Atomic Cluster Expansion, ACE）框架因其严格的数学完备性、高精度和出色的计算效率，成为了构建通用型机器学习势函数的明星方法。

然而，有了一个准确的势函数，只是解决了“力”的问题。要绘制一张可靠的压力-温度（P-T）相图，我们需要计算每个相在不同（P, T）条件下的自由能，并找到自由能最低的稳定相。这涉及到对高维构型空间的彻底探索，传统分子动力学模拟在跨越较高的自由能势垒（如固-固相变）时效率低下。这时，嵌套采样（Nested Sampling, NS）这项源自天体物理和贝叶斯推断的算法，展现了其独特价值。它通过一种“剥洋葱”式的采样策略，系统地探索整个势能面，能够直接、高效地计算配分函数和自由能，特别擅长处理多稳态和相变问题。

本项工作，正是将ACE机器学习势函数与嵌套采样这两大尖端工具深度融合，系统攻克了金属镁在0-600 GPa超高压范围内的相图预测难题。我们不仅验证了已知的HCP到BCC的相变，更在约467 GPa处预测了一个新的BCC到FCC的相变，并绘制了直至熔化温度的热力学曲线。这个过程，是一套完整的“从第一性原理数据到宏观相图”的计算范式，对于研究其他材料在极端条件下的行为具有普遍的参考意义。如果你正在从事高压物理、相变理论或计算材料设计，那么下文对方法细节、实操陷阱和结果分析的拆解，或许能为你扫清一些障碍。

2. 核心方法解析：为什么是ACE+NS？

2.1 原子簇展开（ACE）：构建高精度、可迁移的机器学习势函数

原子簇展开并非一个具体的神经网络模型，而是一个构建原子间描述符的数学框架。它的核心思想是将整个原子系统的势能，表达为所有原子局域环境基函数的对称化乘积之和。听起来抽象，我们可以用一个类比来理解：想象你要描述一个人的面部特征（局域环境），你可以用眼睛、鼻子、嘴巴等基本单元（基函数）的组合来描述。ACE做的就是系统地定义这些“基本单元”（如原子间距、角度等），并确保无论人脸如何旋转（体系发生对称操作），描述符都保持不变（旋转、平移、置换不变性）。

具体到公式，原子i的能量E_i可以写为：

E_i = ∑_B c_B * B(A_i)

其中，A_i是原子i的局域环境，B是一组精心设计的基函数（通常是球谐函数与径向函数的耦合），c_B是需要通过拟合DFT数据来学习的系数。ACE框架的“展开”指的是，这些基函数B是通过对原子邻居的坐标进行多项式展开得到的，其阶数（如C4, O4, D14）控制了描述的复杂度和精度。

为什么选择ACE而不是其他MLIP（如GAP、Neural Network Potentials）？

1. 严格的完备性与收敛性：ACE在数学上被证明是完备的，意味着只要展开阶数足够高，它可以以任意精度逼近任何光滑的势能面。这为系统性地提高精度提供了明确路径。
2. 卓越的计算效率：ACE描述符的计算和能量/力的评估可以高度向量化和优化，相比许多深度神经网络势函数，它在分子动力学模拟中通常更快。
3. 优秀的可迁移性：通过对大量不同构型（晶体、缺陷、液相等）的DFT数据进行训练，ACE势函数能够可靠地外推到训练集未直接覆盖的相空间区域，这对于高压相图探索至关重要。

在本研究中，我们训练了一个C4 O4 D14参数的ACE势函数。这里的参数选择是经验与测试的平衡：C4（4-body correlation order）和 O4（4-body polynomial order）足以捕捉镁在高压下复杂的多体相互作用，而D14（最大截断半径14埃）确保了足够的相互作用范围。训练数据来源于主动学习（Active Learning）循环，初始数据来自一个已知的EAM势函数产生的错误相图（见附录A），这引导采样去探索HCP、BCC、FCC以及液相区域，确保了训练数据的多样性和代表性。

2.2 嵌套采样（NS）：高效遍历势能面，直击自由能计算

自由能是决定相稳定性的终极判据。但对于复杂的凝聚态体系，自由能无法直接测量，计算也极其困难。嵌套采样提供了一条优雅的路径。

NS的核心思想是，将高维构型空间的积分（计算配分函数Z）转化为一维的积分问题。它通过引入一个“约束能量”E的概念来工作：

1. 初始化：在给定的体积（或压力）下，随机生成大量（如数千个）原子构型（称为“活点”）。
2. 迭代“淘汰”：
   • 找出当前所有活点中势能最高的那个构型。
   • 记录其能量E，并将其从活点集合中移除。
   • 从剩余的活点中随机复制一个，并对其施加一个马尔可夫链蒙特卡洛（MCMC）随机行走，但有一个关键约束：新生成的构型势能必须低于刚才淘汰的E。
   • 将这个新构型加入活点集合，保持总活点数不变。
3. 结果：重复步骤2，能量界限E会逐渐降低，活点集合会向势能更低的区域（即更稳定的构型）收缩。在这个过程中，NS自动地、按权重探索了从高能无序态（如液相）到低能有序态（如晶体）的所有可能构型。

最终，我们得到一条能量E与构型空间体积X(E)（其对数正比于熵）的关系曲线。通过统计力学公式，可以直接从中推导出所有热力学量，如亥姆霍兹自由能F、熵S、恒压热容C_P等。C_P在相变点会出现尖峰，这为我们精确定位相变温度提供了清晰信号。*

NS用于相图计算的独特优势：

• 同时处理多相：一次NS运行可以同时采样到固体、液体甚至亚稳态，无需预先指定相。
• 高效跨越势垒：MCMC步长和约束条件使其能够有效跨越固-固相变能垒，而传统MD可能需要过长的模拟时间等待涨落。
• 直接输出自由能：避免了热力学积分等复杂且容易累积误差的后处理步骤。

2.3 方法联动的完整工作流

将ACE和NS结合，形成了一个强大的闭环工作流：

1. 第一性原理数据生成：使用CASTEP软件包，在0-600 GPa选取多个压力点，对HCP、BCC、FCC等候选结构进行结构优化和弹性常数计算，生成高精度的能量、力、应力数据。收敛性测试（平面波截断能、k点网格）至关重要，附录E详细展示了如何确保每个原子能量误差<1 meV。
2. ACE��函数训练与验证：利用QUIP或ACEpotentials.jl等软件，用步骤1的数据训练ACE势。关键验证包括：在训练集和测试集上的能量/力误差；预测已知晶体结构的晶格常数、弹性常数；计算Bain路径（BCC到FCC的变形路径，见附录D）并与DFT结果对比。
3. 嵌套采样计算：使用pymatnest软件包，在目标压力点下，对64个原子的超胞进行NS模拟。活点数（如2000）、MCMC步数等参数需要测试以确保收敛。NS输出每个温度下的平均能量、体积等。
4. 热力学量与相图构建：从NS输出的能量-温度曲线，通过数值微分得到热容C_P，其峰值对应相变温度。同时，可以计算焓H = U + PV和吉布斯自由能G。通过比较不同相在相同（P, T）下的G，即可确定稳定相，绘制相图。
5. 缺陷性质计算（补充验证）：使用训练好的ACE势，在ASE环境中计算空位形成焓、自间隙子形成焓和层错能。这些性质是势函数在点缺陷尺度上精度的试金石，也与材料的力学行为密切相关。

3. 实操要点与避坑指南

3.1 第一性原理计算：精度基石与参数收敛

所有机器学习势函数的“天花板”就是其训练数据的质量。对于高压研究，DFT计算设置必须格外小心。

关键参数与收敛标准：

• 赝势（Pseudopotential）：在超高压下（如600 GPa），电子波函数被强烈压缩，需要验证所用赝势（通常是超软赝势）的可靠性。我们进行了Δ测试（Delta Test，附录F），比较了超软赝势和更“硬”的赝势在压缩体积下的能量曲线差异。结果显示Δ值仅为0.035 meV/atom，证明我们的选择在目标压力范围内是可靠的。
• k点网格与截断能：对于不同的晶体结构和晶胞，需要分别测试。我们采用了固定k点网格的策略（见表IX），而不是根据晶胞大小自动生成。这是因为在NS的晶格动力学过程中，晶胞形状和体积会变化，自动重新生成k点会导致能量出现不连续的“跳跃”，干扰NS的采样。固定网格基于0.015 Å⁻¹的最大间距确定，确保了所有相关结构k点密度的一致性。
• 收敛判据：几何优化的收敛标准需要设得足够紧。我们使用了：总能量变化<0.02 meV，原子受力<1 meV/Å，应力<0.01 GPa，原子位移<0.001 Å。这些严苛的标准确保了训练数据中能量的微小差异是真实的物理效应，而非数值噪声。

避坑提示：切勿在高压计算中直接套用常压下的DFT参数。务必进行系统的收敛性测试，并特别关注赝势在极端压缩下的表现。Δ测试是一个被广泛接受的验证方法。

3.2 ACE势函数训练：数据、描述符与正则化

训练数据库的构建是成功的一半。我们的策略是“主动学习”：

1. 初始种子：从一个已知的、但不准确的镁EAM势函数产生的相图开始（附录A）。这个“错误”的相图（例如预测了不存在的FCC稳定区）恰恰为我们指明了需要重点采样的相空间区域。
2. 迭代丰富：用初始ACE势运行NS或分子动力学，探测其预测不确定性高的区域（如相边界附近、高温液相），将这些新构型用DFT重新计算，加入训练库。
3. 多样性保障：数据库需包含：不同压力下的完美晶体、施加了随机扰动（声子）的晶体、液态构型、点缺陷（空位、间隙子）、堆垛层错以及Bain路径上的中间结构。

描述符与超参数选择：

• 截断半径（r_cut）：需要大于最近邻原子距离的第二或第三近邻。对于高压镁，我们设为14 Å，以确保在600 GPa的极端压缩下仍有足够的描述范围。
• 展开阶数（C, O）：这本质上是精度与计算成本的权衡。我们通过测试发现，对于镁，C4 O4在预测能量、力和应力上已经达到与DFT的误差接近（~1 meV/atom）。盲目提高阶数会增加过拟合风险，且显著增加计算量。
• 正则化（Regularization）：在损失函数中加入L2正则化项，防止系数c_B过大，提高势函数的数值稳定性和外推能力。

验证不止于误差：训练集和测试集上的低误差（RMSE）是必要的，但非充分的。必须进行物理性质验证：

• 晶格动力学：计算声子谱，确保在所有波矢下没有虚频（负频率），否则结构动力学不稳定。
• 弹性常数：与DFT计算的弹性常数对比，确保弹性张量的各向异性关系正确。
• Bain路径：计算BCC沿[001]方向拉伸至FCC路径上的能量变化，与DFT结果严格对比（图19）。这是检验势函数对剪切变形响应是否准确的关键。

3.3 嵌套采样实施：参数、收敛性与误差分析

系统大小与有限尺寸效应：我们主要使用64个原子的超胞。NS的结果存在有限尺寸效应，表现为相变时热容峰是展宽的，而非理想的尖峰。附录B指出，随着系统增大，峰会变尖锐并向低温移动。因此，相变温度的误差棒（error bar）通常取为热容峰的半高全宽（FWHM），这比单纯重复运行NS的统计波动更能反映系统尺寸带来的不确定性。

NS关键参数设置：

• 活点数（N_live）：决定了采样分辨率。一般需要几百到几千。我们使用了2000个活点，在计算成本和结果平滑度之间取得了平衡。
• MCMC步数（N_steps）：每个淘汰迭代中，用于产生新构型的随机行走步数。步数太少，采样不充分，无法跨越能垒；步数太多，计算浪费。通常需要测试以确保构型空间被充分探索。一个经验法则是，让步数足够大，使得新构型与母本构型的均方位移（RMSD）达到晶格常数的量级。
• MCMC步长：需要动态调整，以保持一个合理的接受率（如40%-60%）。pymatnest通常内置了自适应调整算法。

处理势能面“空洞”（PES Holes）：在主动学习初期，或当势函数在某个高能区域拟合不佳时，可能会产生非物理的极低能量点（“空洞”）。NS采样一旦陷入空洞，就会卡住。我们尝试过两种解决方案：

1. 最小键长限制（附录C）：直接拒绝任何出现小于某个阈值的原子间距的构型。这在中等压力下有效，但在极高压力下，液相本身的原子间距也很小，此方法会错误地排除真实液相构型。
2. 委员会模型（Committee Model）：训练多个（如5个）不同的ACE势函数（通过数据子集或不同随机种子）。在NS采样时，如果一个构型被大部分模型预测为高能，但被某一个模型预测为极低能，这很可能是一个“空洞”。我们可以拒绝或惩罚这样的构型。这是更稳健、更物理的方法，也是我们最终采用的方案。

3.4 相变判定与后处理

NS直接输出的是每个温度下的平均势能和平均体积 。焓H =+ P。相变点通过以下特征判定：

• 固-液相变：在焓-温度曲线上出现明显的斜率变化（因为熔化潜热），同时体积发生跃变。热容C_P会出现一个尖锐的峰。
• 固-固相变：焓-温度曲线斜率也会变化（因为热容不同），但通常没有潜热，表现为一个连续的转变。热容C_P会出现一个相对宽泛的峰。BCC<->FCC这类扩散型相变，在有限尺寸体系下，其热容峰可能非常宽，需要仔细分析。

熔化潜热与体积变化计算（附录G）：由于NS数据有噪声，直接取相变温度前后两点的差值不准确。我们的做法是，在相变温度前后各取一段温度区间（如±100-2000 K），分别拟合一条���线到焓-温度数据上。两条直线在相变温度T_m处的差值即为熔化潜热ΔH_m。体积变化ΔV_m用同样方法计算。这种方法平滑了噪声，给出了更可靠的值。

4. 镁高压相图计算结果与讨论

基于上述ACE+NS框架，我们绘制了镁在0-600 GPa范围内的完整P-T相图。

4.1 固-固相变序列

我们的计算清晰地揭示了镁在高压下的相变序列：

1. HCP -> BCC：在约56 GPa（0 K下）发生。这与之前的多数第一性原理研究和部分高压实验观测相符。HCP是镁的常温常压稳定相，在压力下逐渐失去稳定性，让位于更紧密堆积的BCC相。
2. BCC -> FCC：在约468 GPa（0 K下）发生。这是一个新的理论预测。在如此高的压力下，电子结构发生重大重构，FCC结构因其更高的配位数和对称性，在焓上变得比BCC更有利。Bain路径分析（附录D）显示了从BCC到FCC连续变形的能量路径，确认了该相变的过渡机制。

深度分析：为什么是FCC而不是其他结构？在超高压下，材料的性质由电子费米面的拓扑和离子实-电子相互作用主导。FCC结构具有更高的对称性和配位数（12），在某些电子填充情况下，其能带结构可能更稳定。我们的计算通过直接比较HCP、BCC、FCC以及大量堆垛变体（见下文）在目标压力下的焓，确认了FCC的稳定性。

4.2 熔化曲线与热力学性质

我们计算的熔化曲线从常压延伸至600 GPa。与之前的实验和EAM势结果相比，我们的ACE-NS结果：

• 在低压段（< 20 GPa）与实验数据吻合良好。
• 纠正了EAM势在35-40 GPa错误预测的FCC相以及不准确的熔化温度。
• 在高压下，熔化温度持续上升，但上升斜率逐渐趋于平缓。熔化潜热ΔH_m和熔化体积变化ΔV_m随压力增加而减小（图21），这与高压下原子间相互作用增强、液相与固相结构差异减小的物理图像一致。

4.3 堆垛变体与结构搜索

为了排除其他复杂结构（如长周期堆垛结构）在高压下稳定的可能性，我们进行了一项系统性搜索（第H节，表X）。在六方晶胞中，生成了最多12层原子的所有可能堆垛序列（共11,676种唯一结构），在1 GPa下进行松弛并计算其能量。结果表明，在研究的压力范围内，标准的HCP、BCC、FCC结构是最稳定的，没有发现更复杂的堆垛变体能量更低。这增强了我们对所预测相图简洁性的信心。

4.4 与已有研究的对比及意义

我们将结果与此前主要的实验和计算研究进行了全面对比（图16）：

• Stinton et al. (2014)：实验研究，提供了中压段的相边界。
• Moriarty & Althoff (1995), Mehta et al. (2006)：早期的第一性原理研究，预测了HCP-BCC相变压强。
• Ibrahim et al. (2023)：另一项基于ACE势函数的工作，其结果与我们的高度一致，形成了独立的交叉验证。
• Gorman et al. (2022)：最新的实验工作，在太帕（TPa）压力下观测到了镁的开放结构，这超出了我们当前600 GPa的研究范围，但指明了下一步的探索方向。

本工作的核心价值在于，它首次采用统一的、从第一性原理精度出发的机器学习势函数，结合严格统计力学的嵌套采样方法，系统性地预测了镁在超宽压力范围内的完整相图。这套方法学框架具有高度的可移植性，可以应用于其他元素或合金体系，特别是那些实验难以触及的极端条件区域。

5. 常见问题、挑战与解决方案实录

在实际操作这套流程时，会遇到许多在标准教程中不会提及的“坑”。这里记录下我们踩过的一些雷和解决方案。

Q1：DFT训练数据已经非常准确，但训练出的ACE势在NS模拟中还是出现了非物理的“空洞”，导致采样崩溃。怎么办？A1：这是MLIP训练中常见的“外推”问题。即使训练集误差很低，势函数在从未见过的、高能的、非局域构型（如严重扭曲的键、原子穿透）区域行为可能是不可预测的。解决方案：

• 强化主动学习：在生成初始训练数据时，不仅要采样平衡结构，还要有意识地引入非平衡构型。例如，运行高温MD（甚至让部分原子熔化）、对晶体施加大的剪切应变、随机删除或插入原子等。
• 使用委员会模型进行在线过滤：正如我们采用的，在NS采样过程中，用多个模型进行“投票”，拒绝被委员会中多数模型判定为异常低能的构型。这是目前最有效的防护措施。
• 引入物理约束：在势函数形式中内置一些硬约束，如短程强排斥项。但这可能会影响拟合灵活性。

Q2：NS计算非常耗时，尤其是为了收敛需要较大的系统尺寸（如256原子）。有什么加速策略？A2：NS的耗时主要在每一步的MCMC能量评估上。

• 并行化：NS算法天生易于并行。每个活点的MCMC行走是独立的，可以在数百个CPU核心上同时进行。pymatnest支持MPI并行。
• 混合精度与GPU加速：确保你的ACE势评估代码支持GPU计算。像QUIP、ACEpotentials.jl等库在这方面优化得很好。同时，在MCMC过程中，可以使用单精度浮点数进行能量和力的评估，以换取速度，只要最终收集统计数据时用双精度即可。
• 系统尺寸缩放：先在小体系（如32、64原子）上摸清相变的大致位置和NS参数，然后有针对性地对大体系进行关键压力点的计算。并非所有压力点都需要大体系模拟。

Q3：如何确定NS已经“收敛”？活点数N_live和MCMC步数N_steps到底设多少合适？A3：没有绝对标准，但有以下判断方法：

• 重复性：用不同的随机种子运行两次NS，观察得到的焓-温度曲线和热容峰位置是否在误差范围内一致。
• 证据值（Evidence）的稳定性：NS会计算配分函数的对数（ln Z）。增加活点数，如果ln Z的变化小于某个阈值（如0.1），则可以认为收敛。
• 能量分布的探索：观察被淘汰构型的能量序列。如果能量下降曲线在后期出现很长的平台，说明低能区域已被充分探索。如果能量突然快速下降后很久没有新构型被淘汰，可能采样不足。
• 经验法则：N_live至少是体系自由度的几倍。对于64原子的三维体系，2000是一个合理的起点。N_steps需要测试：设置一个值运行，检查MCMC的接受率和产生的构型与母本的RMSD。如果接受率太低或RMSD太小，就增加步数。

Q4：从NS数据中提取相变温度时，热容峰很宽，如何精确定位？A4：对于固-固相变，有限尺寸效应会导致热容峰非常宽。这时，峰顶位置并不总是最可靠的指标，尤其是对于一级相变。更好的方法是：

1. 计算自由能差：直接利用NS输出的配分函数，计算两相在目标温度附近的吉布斯自由能差ΔG(T)。ΔG(T) = 0对应的温度即为相变温度。
2. 使用序参量（Order Parameter）：定义能区分两相的序参量（如BCC的键取向序、FCC的四面体序）。在NS过程中同时计算该序参量的平均值。序参量随温度的突变点，与自由能差为零的点应对应。
3. 多峰拟合：对热容曲线进行多峰（如两个高斯峰）拟合，将两相共存区的贡献分开，取两峰之间自由能相等的温度。

Q5：对于合金或多组分体系，这套方法还能用吗？复杂度会增加多少？A5：完全可以，但复杂度呈指数级增长。

• ACE势函数：ACE框架天然支持多组分。描述符需要包含化学物种信息。训练数据需要覆盖所有可能的成分和有序/无序构型，数据量需求巨大。主动学习变得更加关键。
• 嵌套采样：构型空间不仅包括原子坐标，还包括化学序（Chemical Ordering）。这需要扩展采样变量，例如在MCMC步中不仅要移动原子，还要以一定的概率交换不同种类的原子。这大大增加了采样难度和计算量。可能需要结合复制品交换（Replica Exchange）等技术来提高采样效率。
• 相图类型：从单组分的P-T相图，变为多组分的成分-温度相图或等压截面图。计算量从一条线扩展为一个面或一个体。

最后，分享一点个人体会：机器学习势函数与先进采样方法的结合，正在从根本上改变计算相图学的研究范式。它不再仅仅是DFT的“快速替代品”，而是成为一个能够自主探索未知相空间、发现新现象的强大发现工具。这个领域的门槛正在从“会不会写DFT输入文件”，转向“能否设计高效的数据生成策略、构建稳健的ML模型、并理解统计采样算法的深层原理”。掌握这套从电子结构到宏观性质的完整计算链条，将成为未来计算材料学家的一项核心技能。所有的代码（pymatnest, NS_database_builder）和数据都已开源，希望这份详细的技术复盘，能帮助你更快地上手和拓展自己的工作。