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. 核心方法解析为什么是ACENS2.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严格的完备性与收敛性ACE在数学上被证明是完备的意味着只要展开阶数足够高它可以以任意精度逼近任何光滑的势能面。这为系统性地提高精度提供了明确路径。卓越的计算效率ACE描述符的计算和能量/力的评估可以高度向量化和优化相比许多深度神经网络势函数它在分子动力学模拟中通常更快。优秀的可迁移性通过对大量不同构型晶体、缺陷、液相等的DFT数据进行训练ACE势函数能够可靠地外推到训练集未直接覆盖的相空间区域这对于高压相图探索至关重要。在本研究中我们训练了一个C4 O4 D14参数的ACE势函数。这里的参数选择是经验与测试的平衡C44-body correlation order和 O44-body polynomial order足以捕捉镁在高压下复杂的多体相互作用而D14最大截断半径14埃确保了足够的相互作用范围。训练数据来源于主动学习Active Learning循环初始数据来自一个已知的EAM势函数产生的错误相图见附录A这引导采样去探索HCP、BCC、FCC以及液相区域确保了训练数据的多样性和代表性。2.2 嵌套采样NS高效遍历势能面直击自由能计算自由能是决定相稳定性的终极判据。但对于复杂的凝聚态体系自由能无法直接测量计算也极其困难。嵌套采样提供了一条优雅的路径。NS的核心思想是将高维构型空间的积分计算配分函数Z转化为一维的积分问题。它通过引入一个“约束能量”E的概念来工作初始化在给定的体积或压力下随机生成大量如数千个原子构型称为“活点”。迭代“淘汰”找出当前所有活点中势能最高的那个构型。记录其能量E并将其从活点集合中移除。从剩余的活点中随机复制一个并对其施加一个马尔可夫链蒙特卡洛MCMC随机行走但有一个关键约束新生成的构型势能必须低于刚才淘汰的E。将这个新构型加入活点集合保持总活点数不变。结果重复步骤2能量界限E会逐渐降低活点集合会向势能更低的区域即更稳定的构型收缩。在这个过程中NS自动地、按权重探索了从高能无序态如液相到低能有序态如晶体的所有可能构型。最终我们得到一条能量E与构型空间体积X(E)其对数正比于熵的关系曲线。通过统计力学公式可以直接从中推导出所有热力学量如亥姆霍兹自由能F、熵S、恒压热容C_P等。C_P在相变点会出现尖峰这为我们精确定位相变温度提供了清晰信号。*NS用于相图计算的独特优势同时处理多相一次NS运行可以同时采样到固体、液体甚至亚稳态无需预先指定相。高效跨越势垒MCMC步长和约束条件使其能够有效跨越固-固相变能垒而传统MD可能需要过长的模拟时间等待涨落。直接输出自由能避免了热力学积分等复杂且容易累积误差的后处理步骤。2.3 方法联动的完整工作流将ACE和NS结合形成了一个强大的闭环工作流第一性原理数据生成使用CASTEP软件包在0-600 GPa选取多个压力点对HCP、BCC、FCC等候选结构进行结构优化和弹性常数计算生成高精度的能量、力、应力数据。收敛性测试平面波截断能、k点网格至关重要附录E详细展示了如何确保每个原子能量误差1 meV。ACE函数训练与验证利用QUIP或ACEpotentials.jl等软件用步骤1的数据训练ACE势。关键验证包括在训练集和测试集上的能量/力误差预测已知晶体结构的晶格常数、弹性常数计算Bain路径BCC到FCC的变形路径见附录D并与DFT结果对比。嵌套采样计算使用pymatnest软件包在目标压力点下对64个原子的超胞进行NS模拟。活点数如2000、MCMC步数等参数需要测试以确保收敛。NS输出每个温度下的平均能量、体积等。热力学量与相图构建从NS输出的能量-温度曲线通过数值微分得到热容C_P其峰值对应相变温度。同时可以计算焓H U PV和吉布斯自由能G。通过比较不同相在相同P, T下的G即可确定稳定相绘制相图。缺陷性质计算补充验证使用训练好的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势函数训练数据、描述符与正则化训练数据库的构建是成功的一半。我们的策略是“主动学习”初始种子从一个已知的、但不准确的镁EAM势函数产生的相图开始附录A。这个“错误”的相图例如预测了不存在的FCC稳定区恰恰为我们指明了需要重点采样的相空间区域。迭代丰富用初始ACE势运行NS或分子动力学探测其预测不确定性高的区域如相边界附近、高温液相将这些新构型用DFT重新计算加入训练库。多样性保障数据库需包含不同压力下的完美晶体、施加了随机扰动声子的晶体、液态构型、点缺陷空位、间隙子、堆垛层错以及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采样一旦陷入空洞就会卡住。我们尝试过两种解决方案最小键长限制附录C直接拒绝任何出现小于某个阈值的原子间距的构型。这在中等压力下有效但在极高压力下液相本身的原子间距也很小此方法会错误地排除真实液相构型。委员会模型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. 镁高压相图计算结果与讨论基于上述ACENS框架我们绘制了镁在0-600 GPa范围内的完整P-T相图。4.1 固-固相变序列我们的计算清晰地揭示了镁在高压下的相变序列HCP - BCC在约56 GPa0 K下发生。这与之前的多数第一性原理研究和部分高压实验观测相符。HCP是镁的常温常压稳定相在压力下逐渐失去稳定性让位于更紧密堆积的BCC相。BCC - FCC在约468 GPa0 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 与已有研究的对比及意义我们将结果与此前主要的实验和计算研究进行了全面对比图16Stinton 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. 常见问题、挑战与解决方案实录在实际操作这套流程时会遇到许多在标准教程中不会提及的“坑”。这里记录下我们踩过的一些雷和解决方案。Q1DFT训练数据已经非常准确但训练出的ACE势在NS模拟中还是出现了非物理的“空洞”导致采样崩溃。怎么办A1这是MLIP训练中常见的“外推”问题。即使训练集误差很低势函数在从未见过的、高能的、非局域构型如严重扭曲的键、原子穿透区域行为可能是不可预测的。解决方案强化主动学习在生成初始训练数据时不仅要采样平衡结构还要有意识地引入非平衡构型。例如运行高温MD甚至让部分原子熔化、对晶体施加大的剪切应变、随机删除或插入原子等。使用委员会模型进行在线过滤正如我们采用的在NS采样过程中用多个模型进行“投票”拒绝被委员会中多数模型判定为异常低能的构型。这是目前最有效的防护措施。引入物理约束在势函数形式中内置一些硬约束如短程强排斥项。但这可能会影响拟合灵活性。Q2NS计算非常耗时尤其是为了收敛需要较大的系统尺寸如256原子。有什么加速策略A2NS的耗时主要在每一步的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对于固-固相变有限尺寸效应会导致热容峰非常宽。这时峰顶位置并不总是最可靠的指标尤其是对于一级相变。更好的方法是计算自由能差直接利用NS输出的配分函数计算两相在目标温度附近的吉布斯自由能差ΔG(T)。ΔG(T) 0对应的温度即为相变温度。使用序参量Order Parameter定义能区分两相的序参量如BCC的键取向序、FCC的四面体序。在NS过程中同时计算该序参量的平均值。序参量随温度的突变点与自由能差为零的点应对应。多峰拟合对热容曲线进行多峰如两个高斯峰拟合将两相共存区的贡献分开取两峰之间自由能相等的温度。Q5对于合金或多组分体系这套方法还能用吗复杂度会增加多少A5完全可以但复杂度呈指数级增长。ACE势函数ACE框架天然支持多组分。描述符需要包含化学物种信息。训练数据需要覆盖所有可能的成分和有序/无序构型数据量需求巨大。主动学习变得更加关键。嵌套采样构型空间不仅包括原子坐标还包括化学序Chemical Ordering。这需要扩展采样变量例如在MCMC步中不仅要移动原子还要以一定的概率交换不同种类的原子。这大大增加了采样难度和计算量。可能需要结合复制品交换Replica Exchange等技术来提高采样效率。相图类型从单组分的P-T相图变为多组分的成分-温度相图或等压截面图。计算量从一条线扩展为一个面或一个体。最后分享一点个人体会机器学习势函数与先进采样方法的结合正在从根本上改变计算相图学的研究范式。它不再仅仅是DFT的“快速替代品”而是成为一个能够自主探索未知相空间、发现新现象的强大发现工具。这个领域的门槛正在从“会不会写DFT输入文件”转向“能否设计高效的数据生成策略、构建稳健的ML模型、并理解统计采样算法的深层原理”。掌握这套从电子结构到宏观性质的完整计算链条将成为未来计算材料学家的一项核心技能。所有的代码pymatnest, NS_database_builder和数据都已开源希望这份详细的技术复盘能帮助你更快地上手和拓展自己的工作。
