机器学习势函数与量子声子方法结合：第一性原理预测SrTiO3氧同位素效应

1. 项目概述：当机器学习势函数遇见量子顺电体

在凝聚态物理和材料计算领域，钛酸锶（SrTiO3）一直是一个充满魅力的“问题儿童”。它在低温下展现出量子顺电性——一种介电常数极高但自发极化被量子涨落抑制的奇特状态，而非像传统铁电体那样发生相变。更令人着迷的是，仅仅将氧原子从常见的16O替换为稍重的18O，就能“唤醒”其内在的铁电性，诱导出一个清晰的低温铁电相变。这个现象，即氧同位素效应，自被发现以来就吸引了大量实验和理论研究，因为它直接触及了量子力学、晶格动力学和电子结构之间微妙的相互作用核心。

然而，从第一性原理出发定量预测这一现象，长期以来都是一个巨大的挑战。问题的根源在于，要准确描述SrTiO3，必须同时处理两大难题：核的量子效应和晶格振动的强非谐性。传统的密度泛函理论（DFT）分子动力学模拟，其计算成本在考虑量子核效应（如路径积分）时变得难以承受。而忽略这些效应，又会错误地预测出铁电基态，与实验观察到的量子顺电性背道而驰。

近年来，计算材料学领域的两股技术浪潮为解决这一困境带来了曙光。一是机器学习原子间势的发展，它能够以接近经典力场的计算成本，复现出第一性原理级别的精度，使得对数千原子体系进行大规模采样成为可能。二是随机自洽谐波近似等方法的成熟，这类方法能够在自洽框架下，同时处理温度、非谐性和量子核效应，直接计算体系的自由能及其对结构的影响。

我们这个工作的核心，就是将MLIPs与SSCHA方法强力结合，构建了一个全新的计算框架。我们不再依赖于拟合实验数据的唯象模型，而是试图从“第一性原理”出发，仅输入原子种类和位置，去预测SrTiO3在16O和18O情况下的完整相图，特别是其软模频率随温度、体积和晶格参数的变化。我们的目标很明确：验证这个结合了前沿机器学习与量子声子方法的计算方案，能否定性甚至定量地复现出那个著名的氧同位素效应，并深入理解量子涨落如何在原子尺度上“操控”材料的宏观相变行为。

2. 核心方法解析：MLIPs与SSCHA的强强联合

要理解我们如何攻克SrTiO3的难题，首先得拆解我们手中的两件核心“工具”：机器学习原子间势和随机自洽谐波近似。它们各自解决了传统计算中的不同瓶颈，组合起来才形成了我们这次研究的完整技术栈。

2.1 机器学习原子间势：从“计算黑盒”到“高效代理”

传统的第一性原理计算，如DFT，通过求解电子薛定谔方程来获得体系的能量和原子受力。这个过程虽然精度高，但计算量随原子数呈三次方或更高次方增长，对于需要大量采样以统计非谐效应和量子涨落的有限温度模拟来说，简直是天文数字。

机器学习原子间势的思路则是一种“曲线救国”。我们不再每次都从头求解复杂的量子力学方程，而是先用DFT计算一个相对较小但具有代表性的数据集（包含不同原子构型的能量、力和应力），然后用一个灵活的机器学习模型（如神经网络、高斯过程或矩张量势）去学习从原子构型到这些物理量的映射关系。一旦模型训练完成，它就像一个高度精确的“代理模型”或“力场”，可以在毫秒级别内给出新构型的能量和受力，而其精度却非常接近背后的第一性原理方法。

在我们的工作中，我们直接沿用了Verdi等人开发的一个关键MLIP。这个势函数并非普通的DFT势，而是一个基于随机相位近似训练的势函数。这里就涉及到一个重要的技术选择：为什么是RPA？因为早期的研究发现，对于SrTiO3这种电子关联效应重要的体系，常见的DFT泛函（如PBE、PBEsol甚至SCAN）要么严重高估晶格常数，要么错误地预测铁电基态。RPA作为一种更高级的交换关联近似，虽然计算极其昂贵，但被证明能更好地描述SrTiO3的量子顺电软模行为。Verdi等人采用了一种“Delta-learning”策略：先训练一个基于PBEsol的便宜MLIP，再训练第二个MLIP来学习RPA与PBEsol结果之间的差异，最后将两者结合，得到一个高效的“RPA级”MLIP。我们为了兼容更新的VASP 6.4.0版本，基于他们的数据和参数重新拟合了这个势函数。

注意：MLIP的成败极度依赖于训练数据的质量和广度。对于SrTiO3，其势能面在软模方向非常平坦，这意味着训练数据必须充分覆盖可能的大振幅原子位移，否则MLIP在预测这些区域的能量和受力时会产生巨大误差，直接影响相变行为的判断。

2.2 随机自洽谐波近似：让量子核与非谐性“自洽”起来

有了高效的力场，我们还需要一个能处理量子核效应和非谐性的理论框架。这就是随机自洽谐波近似大显身手的地方。

传统的晶格动力学计算基于谐波近似，即假设原子在平衡位置附近做简谐振动。这显然无法描述SrTiO3中导致相变的强非谐势能面（典型的双阱势）。SSCHA的核心思想是变分法。它不再寻找原子的静态平衡位置，而是寻找一个高斯型核密度矩阵，该密度矩阵能使系统的有限温度自由能最小化。这个高斯分布的均值和协方差矩阵（对应于有效位置和有效力常数）是通过自洽迭代确定的。

具体流程可以概括为：

1. 初始化：从一个谐波力常数矩阵和平衡位置开始。
2. 采样：根据当前的高斯分布，随机生成大量（成千上万）的原子构型（即一个“系综”）。
3. 力计算：用MLIP快速计算这些采样构型中每个原子所受的力。
4. 更新：利用这些力和构型，通过统计方法更新高斯分布的均值和协方差矩阵，目标是使系统的自由能估计值降低。
5. 迭代：重复步骤2-4，直到均值和力常数矩阵收敛。

SSCHA的妙处在于，它通过这个自洽过程，自动包含了量子涨落（通过核的有限宽度）和非谐效应（通过在高斯分布下对真实势能面的平均）。最终收敛时，我们得到的是考虑了所有量子与热效应的“有效”平衡结构和“有效”声子谱（通过计算Hessian矩阵得到）。如果某个振动模式的频率变为虚数，就预示着结构失稳，可能发生相变。

2.3 技术栈整合与工作流程

我们的完整计算流程是上述两个方法的深度集成：

1. 体系构建：我们针对四方相SrTiO3（考虑了低于105K的反倾转相变）构建了3x3x3的超胞，包含270个原子。分别设置氧原子质量为16和18原子质量单位。
2. SSCHA迭代计算：
   • 我们从谐波解出发，进行SSCHA几何弛豫，但固定了非极性对称性，即不允许体系自发极化。八面体的倾转角在模拟中被允许弛豫。
   • 为了平衡精度和效率，我们采用了分阶段增加采样量的策略：先使用1000个结构的系综进行零温计算，收敛后重启，使用10000个结构，最后使用40000个结构。这种策略能快速接近收敛区域，再用大系综精细调整。
   • 收敛判据设定为：辅助SSCHA力常数矩阵和 centroid 位置相关的能量梯度标准差的1e-4倍。
3. 频率与相图计算：收敛后，我们计算Hessian矩阵以获得声子频率，特别是关注面内极化的Eu模和面外极化的A2u模这两个软模。通过系统性地改变晶胞体积和c/a轴比（四方畸变度），我们为不同温度和同位素质量���制了完整的Eu软模频率相图。
4. 双势阱与极化计算：对于靠近相变边界的结构，我们利用SSCHA得到的声子模式本征矢量，沿着软模坐标方向“冻结”该模式，计算体系能量随模式振幅（可转换为极化强度）的变化曲线，从而直接得到双势阱的深度和平衡极化强度。

这套组合拳的优势在于，MLIP解决了SSCHA所需海量力评估的计算瓶颈，而SSCHA则为MLIP提供了一个严格处理量子与热涨落的物理框架，使得大规模、有限温度、包含量子效应的第一性原理品质模拟成为可能。

3. 计算细节与参数选择：在精度与成本的钢丝上行走

进行这样一项研究，每一个参数的选择背后都是一次精度与计算成本的权衡。下面我详细拆解几个关键环节中的实操要点和背后的考量。

3.1 超胞与收敛性测试：规模与采样量的博弈

我们选择3x3x3的超胞（270原子）并非随意。对于SrTiO3的软模，其波矢位于布里渊区中心（Γ点）。使用足够大的超胞可以更好地模拟Γ点的动力学，并容纳可能的长程相互作用。同时，270个原子对于MLIP+SSCHA的计算来说是一个在现有算力下可管理的规模。我们曾测试过2x2x2的超胞，发现某些声子频率，特别是低频软模，对超胞尺寸仍有一定敏感性，3x3x3则能给出更稳定的结果。

系综大小的选择是SSCHA计算的核心。SSCHA的精度严重依赖于蒙特卡洛采样的质量。我们采用的分阶段策略（1000 -> 10000 -> 40000）是一种实用技巧。初始阶段用较小系综快速逼近收敛区域，避免在错误的方向上浪费计算资源。最终阶段使用4万个结构，是为了确保Hessian矩阵计算的准确性。我们引入了一个关键指标：Kong-Liu比率。这个比率评估了采样构型在重加权计算中的有效性，比率越高（我们要求>95%），说明采样越充分，估计的Hessian越可靠。如果达不到这个标准，我们就需要用新的随机种子生成一个新的4万构型系综重新计算，以克服随机噪声。

实操心得：SSCHA的随机性是一个不容忽视的误差来源。我们重复计算了相图中大部分点（使用不同的随机种子），发现软模频率的中位数差异约为0.5 cm⁻¹，但90%分位数可达1.5 cm⁻¹，个别 outlier 的差异更大。这意味着，对于频率在0 cm⁻¹附近的点（即相变边界），其“铁电”或“顺电”的判定可能因随机噪声而摇摆。因此，在分析相图边界时，必须将这种不确定性考虑在内，结论应基于统计趋势而非单个数据点。

3.2 结构选择与实验对标：在理论的“偏差”中寻找“正确”

一个关键的挑战是：我们的计算基于RPA-MLIP，而RPA本身已知会系统性高估晶格体积。如果我们直接使用实验测得的晶格参数进行计算，得到的结果很可能因为体积偏差而失真。因此，我们不能简单地用实验晶格常数作为输入。

我们的策略是“以终为始”，通过结果来反推合理的计算结构。具体步骤如下：

1. 目标锁定：我们以实验测得的低温下ST16O的Eu和A2u软模频率（分别约7.1 cm⁻¹和16.5 cm⁻¹）为黄金标准。
2. 参数扫描：我们在一个合理的体积和c/a比值范围内进行网格化扫描，对每个点进行SSCHA计算，得到其软模频率。
3. 最优匹配：计算每个扫描点的频率与实验目标频率的欧氏距离（平方和），选择距离最小的那个点所对应的体积和c/a比，作为我们“最匹配实验声子行为”的理论结构。
4. 应用与验证：将这个“最优结构”的参数，同时应用于ST16O和ST18O的计算中。我们验证了在这个结构下，ST16O的软模频率在0-10K之间基本饱和，符合量子顺电体的特征。这就为我们研究ST18O的相变提供了一个自洽的基准。

这个被选中的结构具有体积V=59.49 ų和c/a=1.0027。请注意，这与1.5K下的实验结构（V≈59.38 ų, c/a≈1.0005）确实存在偏差，但重要的是，它复现了核心的物理现象——软模行为。这体现了计算中的一种常见思路：当绝对精度因方法局限难以达到时，确保相对趋势和物理机制的正确性更为关键。

3.3 双势阱拟合与极化计算：从能量曲线到宏观物理量

为了定量研究铁电相变的特征，我们需要从SSCHA计算出的能量-模式振幅曲线中提取双势阱的深度、位置以及对应的自发极化。

1. 模式冻结：利用收敛后的SSCHA计算得到的Eu软模本征矢量，我们沿着该模式的方向，对原子位置进行集体位移，构造一系列具有不同极化强度的赝结构。
2. 单点能计算：对这些冻结了软模的赝结构，用MLIP计算其总能量。这里需要注意的是，这些计算是静态计算，不再进行SSCHA迭代，因为我们已经固定了原子的集体位移模式。
3. 曲线拟合：将能量相对于极化强度（或模式振幅）作图，通常会看到一个典型的双阱形状。我们用一个四次多项式（或对称的双阱势函数）去拟合能量最低点附近的区域。拟合的范围需要谨慎选择：范围太小，无法捕捉势阱全貌；范围太大，高振幅区域可能涉及势能面中其他未被考虑的非谐项，导致拟合失真。我们的策略是优先保证双阱区域（即两个极小值点及中间的势垒）的拟合质量。
4. 极化转换：模式振幅是一个无量纲的量。为了得到宏观的自发极化强度（单位：μC/cm²），我们使用了VASP的密度泛函微扰理论（DFPT）模块，计算了高对称结构下的玻恩有效电荷。然后，通过玻恩有效电荷张量和原子的位移量，可以计算出整个超胞的极化矢量。对于沿着特定软模方向的位移，其极化大小与模式振幅成正比。

通过这套流程，我们得到了ST18O在0K下的自发极化约为6.2 μC/cm²。作为对比，Filipič等人的实验测量值（对94.7% 18O富集的晶体外推至0K）约为3.37 μC/cm²。我们的计算值在数量级上与实验一致，但存在约两倍的偏差。这个偏差可能来源于RPA对体积的高估、SSCHA高斯近似的局限性，以及极化计算本身的不确定性。

4. 结果分析与讨论：一幅由计算绘制的量子相变图景

基于上述严谨的计算框架，我们得到了一系列揭示SrTiO3氧同位素效应微观物理的图像。这些结果不仅定性地复现了实验，更提供了相空间中的全局视角。

4.1 核心发现：量子涨落是开关，同位素质量是扳机

我们最核心的发现直观地展示在能量-极化曲线中。对于选定的“最优结构”：

• 纯RPA-MLIP（无SSCHA）：计算显示Eu软模频率为很大的虚数，能量曲线呈现一个很深的双势阱，错误地预测ST16O为铁电体。这印证了忽略量子核效应会导致灾难性错误。
• RPA-MLIP + SSCHA (ST16O)：引入SSCHA后，量子涨落“撑起”了原子，使其不再局限于双势阱的某一个最小值，而是在两个阱之间通过量子隧穿或零点振动进行涨落。这使得双势阱被有效填平，Eu模频率变为一个小的正实数（~10.9 cm⁻¹），体系稳定在顺电态。尽管能量曲线仍显示出强烈的非谐性（非抛物线形），但已没有稳定的极化最小值。
• RPA-MLIP + SSCHA (ST18O)：当氧原子质量从16增加到18，量子涨落的幅度减弱（因为更重的粒子波动性降低）。这使得系统无法再抵抗双势阱的吸引，Eu模频率重新变为虚数（i 17 cm⁻¹），并在能量曲线上出现一个更浅的双势阱，对应一个有限的自发极化。这清晰地表明，18O的替换通过抑制量子涨落，诱导出了铁电序。

这个结果从第一性原理上直接证实了长期以来对SrTiO3量子顺电性起源的解释：是氧原子核的量子零点运动（涨落）压制了本应发生的铁电失稳。同位素效应之所以显著，正是因为质量的改变直接调制了量子涨落的强度。

4.2 温度依赖性与相变温度预测

我们进一步计算了ST18O的Eu软模频率和双势阱深度随温度的变化。结果显示，随着温度从0K升高，软模频率的虚部绝对值逐渐减小（即不稳定性减弱），双势阱深度也逐渐变浅。大约在40-45K时，双势阱的深度已经小于SSCHA计算本身的误差范围。到50K时，二次项变为正，双势阱消失，转变为单势阱，标志着顺电相的恢复。

由此，我们预测ST18O的铁电-顺电相变温度Tc约为40K。这与实验报道的25-26K相比有所高估。这个偏差是系统性的，主要可能源于两个层面：一是RPA对晶格体积的高估，这通常会增强铁电不稳定性（因为更大的体积往往有利于铁电性）；二是SSCHA方法本身的近似，它用高斯分布来近似核密度矩阵，对于涉及势阱间隧穿的过程可能描述不足。尽管如此，40K的预测值在定性上是正确的，它确认了相变的存在并给出了合理的温度尺度。

4.3 全局相图：结构、温度与同位素的协同效应

为了更全面地理解相变条件，我们绘制了Eu软模频率随体积和c/a比变化的相图，涵盖了0K、30K、50K三个温度以及两种同位素。这幅相图蕴含了丰富的信息：

1. 结构效应：在任何给定的温度和同位素下，增大体积或减小c/a比（即增大a轴长度）都会促进铁电相的出现（蓝色区域扩大）。这与实验上通过化学掺杂（如Ca掺杂）或外延应变来诱导SrTiO3铁电性的观察完全一致。体积膨胀减轻了原子间的“拥挤”，使得极化的原子位移更容易发生。
2. 温度效应：升温会抑制铁电性，将相变边界向更大体积和更大c/a比的方向推移。这是典型位移型相变的特征：热涨落破坏了有序的极化排列。
3. 同位素效应：这是相图最精彩的部分。在每一个温度下，ST18O的铁电相稳定区域都比ST16O更广阔。我们将ST18O的相边界（黄色虚线）叠加在ST16O的相图上，可以清晰地看到一条“沟壑”：在相同的结构参数下，ST16O可能是顺电的（红色），而ST18O却已进入铁电区（蓝色）。这定量地展示了同位素替换如何有效地拓宽了铁电相的稳定范围。

为了更直观地展示同位素效应的强度，我们计算了ST18O与ST16O Eu模频率差异的绝对值。结果显示，频率差异在相变边界附近最为显著。这符合临界现象的一般规律：在相变点附近，系统对外部扰动（如质量变化）异常敏感，微小的扰动就能引发巨大的响应。相比之下，如果在计算中不使用SSCHA（即忽略量子与非谐效应），最大的频率差异也只有5 cm⁻¹，远小于SSCHA给出的结果。这再次凸显了包含量子与非谐效应对于正确描述同位素效应的决定性作用。

4.4 与实验结构的对比及计算局限性

我们将实验测得的晶体结构（不同温度下的体积和c/a比）标记在相图上（白色十字）。有趣的是，对于ST16O，其实验结构点恰好落在顺电相区域内；而对于ST18O，其实验点则位于铁电相的边界附近。尽管我们的理论结构（黑十字）与实验结构（白十字）在具体数值上有偏移（主要源于RPA的体积误差），但两者在定性上完全一致：ST16O是顺电的，ST18O是（或非常接近）铁电的。更重要的是，实验结构随温度变化的趋势（体积收缩，c/a比变化）与我们为匹配频率而选择的理论结构的变化趋势是一致的。

这项工作也清晰地揭示了当前计算方法的局限性：

1. 电子结构方法的精度瓶颈：RPA高估体积，SCAN泛函虽能给出更好的结构但描述软模行为不佳。目前尚无一种普适的电子结构方法能同时完美解决这两个问题。这是当前材料计算领域的一个核心挑战。
2. SSCHA的近似与随机性：高斯近似无法完美描述双势阱间的隧穿效应。计算的随机噪声在相变边界附近引入了不可忽视的不确定性。
3. 应力训练的缺失：我们的RPA-MLIP没有在训练集中包含应力数据，这使得在SSCHA中完全弛豫晶胞自由度变得困难，我们只能进行固定晶胞形状的弛豫。

这些局限性共同导致了我们对Tc和自发极化值的定量预测存在偏差。它们提醒我们，尽管MLIP+SSCHA的组合带来了革命性的能力，但对于像SrTiO3这样极其“娇气”的体系，做出完全无需实验校准的定量预测依然非常困难。

5. 方法论的启示与未来展望

这项研究不仅仅是对SrTiO3氧同位素效应的一次成功模拟，更是一次对“计算驱动发现”新范式的有力展示。MLIP与SSCHA的结合，为我们打开了一扇窗，得以窥见那些同时涉及电子、声子、量子与热涨落的复杂物理过程。

从方法论角度看，我们的工作验证了这种技术路线在处理“量子顺电性”这类难题上的可行性。它表明，通过从第一性原理出发，结合先进的机器学习与统计物理方法，我们可以在没有预设序参量或唯象模型的情况下，预测出复杂的量子相变行为。这对于探索其他量子顺电体（如KTaO3）或更广泛的关联量子材料具有重要的借鉴意义。

这项研究也为理解SrTiO3中与铁电涨落相关的超导同位素效应提供了新的计算基础。既然我们能够模拟氧同位素对铁电软模的直接影响，那么未来就有可能进一步研究这些软模涨落如何影响电子配对，从而为理解该体系的超导机制提供微观见解。

展望未来，这个领域的发展方向是清晰的：

1. 发展更精确的电子结构方法：需要能够同时准确描述结构参数和电子关联（尤其是范德瓦尔斯相互作用和激发态效应）的新一代泛函或多体方法。
2. 超越SSCHA的高斯近似：将路径积分分子动力学与MLIP结合，可以更严格地处理核的量子效应，特别是隧穿过程，这对于厘清ST18O中铁电相变是纯粹的位移型还是包含有序-无序成分至关重要。
3. 构建更通用、更稳健的MLIP：扩大训练数据集，包含应力、不同电荷态、缺陷构型等，并发展主动学习策略，使MLIP能在模拟过程中自动发现并学习势能面的未知区域。
4. 研究更复杂的调控手段：利用此框架，可以系统研究应变、电场、掺杂浓度等多种外场对量子顺电和铁电相变的耦合影响，为功能器件的设计提供理论蓝图。

尽管存在定量预测的挑战，但我们的工作无疑迈出了关键的一步。它证明，通过精心整合最前沿的计算工具，我们已经有能力对材料中最微妙、最核心的量子现象进行基于第一性原理的探索和预测。这不仅是计算能力的胜利，更是计算思想的一次进阶——从解释已知，走向预测未知。