机器学习势函数与量子热浴结合：精准模拟钛酸钡相变中的核量子效应

1. 项目概述：当机器学习势函数遇上量子热浴

在计算材料科学领域，我们一直面临着一个核心矛盾：精度与效率的权衡。研究像钛酸钡（BaTiO₃）这样的经典铁电材料相变，我们需要在原子尺度上追踪成千上万个原子在温度和压力下的动态演化。第一性原理分子动力学（AIMD）虽然能提供高精度的势能面描述，但其计算成本是天文数字，通常只能处理几百个原子、几个皮秒的模拟，这对于研究相变这种涉及长程关联和缓慢动力学的过程，无异于杯水车薪。而传统经验势函数虽然快，但其精度和可迁移性常常令人担忧，尤其是在描述复杂的多体相互作用和相变临界区域时。

这个矛盾，正是机器学习势函数（MLP）大显身手的地方。它本质上是一个“超级插值器”，通过神经网络等模型，从相对少量的高精度第一性原理计算数据中，学习到原子构型与系统能量、原子受力之间的复杂映射关系。一旦训练完成，MLP在预测时的计算代价仅相当于一个经验势函数，但其精度却可以无限逼近第一性原理计算。这就像培养了一位“计算专家”，它熟读了所有高深的理论典籍（第一性原理数据），并能以极快的速度解答新的问题（预测新构型的能量和受力）。

然而，当我们把目光投向低温或涉及轻元素的体系时，另一个物理效应——核量子效应（NQE）——就变得不可忽视了。在经典分子动力学（MD）中，原子核被当作遵循牛顿力学的经典粒子处理。但在量子力学中，原子核具有波粒二象性，存在零点能和量子隧穿等效应。对于钛酸钡，其中的氧原子质量较轻，在低温下其量子涨落会显著影响晶格振动、铁电极化乃至相变温度。忽略NQE，可能会导致对相变边界、有序-无序转变温度的预测出现系统性偏差。

量子热浴（QTB）方法，正是为了在经典MD框架内“植入”量子效应而生的巧妙工具。它通过在朗之万方程的热浴项中引入一个与频率相关的有色噪声，使得体系能量的涨落满足量子涨落-耗散定理，从而让经典振子的能量分布逼近量子谐振子的结果。简单说，QTB让经典MD模拟“听”起来像量子体系。将高精度的MLP与能描述NQE的QTB结合（即QTB-MLPMD），我们便获得了一把利器，既能高效处理大体系，又能相对准确地捕捉量子效应，从而定量评估像压力-温度相图中核量子效应的具体影响。这对于深入理解铁电材料的微观机理、预测其在外场下的响应行为，具有关键意义。

2. 核心思路与技术选型解析

2.1 为什么选择机器学习势函数作为基石？

在钛酸钡相变研究中，选择MLP而非其他势函数，是基于几个深层次的考量。首先，钛酸钡的相变涉及复杂的序参量耦合，包括铁电位移（Ti原子相对于氧八面体中心的偏移）、氧八面体的旋转以及晶格应变。这些序参量相互耦合，其对应的势能面存在多个浅的极小值（对应不同相），能垒高度对温度和压力极其敏感。传统的刚性离子模型或壳层模型势函数，其函数形式是预先设定的，难以灵活捕捉这种复杂、多尺度的耦合关系。而MLP，特别是基于原子局部环境描述符（如SOAP、ACE或Behler-Parrinello对称函数）的神经网络势，具有强大的泛化能力和高维拟合能力，可以从数据中自动学习到这些复杂的相互作用。

其次，研究的核心目标之一是构建精确的压力-温度（P-T）相图。这要求势函数在从立方相到四方相、正交相乃至三方相的整个构型空间内都保持高精度，尤其是在相边界附近，能量的微小差异就会导致相变预测的巨变。MLP的训练策略可以针对性地在相变路径附近进行主动学习或增强采样，密集采集第一性原理数据，确保在这些关键区域有出色的预测能力。从补充材料图S1-S3可以看出，研究者对比了MLP预测与第一性原理计算（FP）在能量、原子力（特别是Ti离子受力）和应力上的一致性。这些图显示数据点紧密分布在y=x对角线两侧，相关系数（R²）极高，这直接验证了所构建的MLP在整个模拟采样空间（包括不同压力下的4×4×4和5×5×5超胞）都具有接近第一性原理的精度，为后续的大规模模拟奠定了可信的基础。

注意：评估MLP精度时，不能只看总能量的误差，原子力和应力的预测精度往往更为关键。因为MD模拟的动力学轨迹由受力直接驱动，力的误差会直接导致错误的原子运动轨迹和相变动力学。应力误差则直接影响压力计算和晶格常数的弛豫。因此，像文中那样系统性地验证能量、力、应力三项指标，是评估MLP是否适用于相变研究的黄金标准。

2.2 量子热浴方法：在经典框架内引入量子灵魂

理解了MLP作为“准确引擎”的价值后，我们需要一个能模拟量子效应的“变速箱”，这就是量子热浴。经典MD使用经典谐振子模型，其平均能量在温度T下为k_B T，在低温下趋于零。但量子谐振子有零点能，即使在绝对零度，能量也不为零。QTB方法的核心思想是修改经典朗之万方程中的随机力项。

经典的朗之万方程是：m_i dv_i/dt = F_i - γ m_i v_i + R_i(t)。其中R_i(t)是高斯白噪声，满足 <R_i(0)R_j(t)> = 2γ m_i k_B T δ_ij δ(t)。为了引入量子效应，QTB将白噪声替换为有色噪声，其自相关函数与频率相关噪声谱密度相连：<R_i(0)R_j(t)> = (m_i γ / π) ∫_0^∞ dω J(ω) cos(ωt) δ_ij。这里的关键是谱密度J(ω)。为了重现量子谐振子的能量期望，需要令J(ω)满足量子涨落-耗散定理：J(ω) = (ħω/2) coth(ħω/(2k_B T))。

实际操作中，我们通过傅里叶变换在频率空间生成满足该谱密度的随机力，再变换回时域加到原子上。这样，每个振动模式所感受到的热涨落就具有了量子特征，体系的总动能和势能分布会更接近量子结果。QTB是一种近似方法，它对于描述非简谐性不强的体系中的核量子效应（如零点能、量子涨落对热容的影响）非常有效，且计算开销仅比经典MD略高，远比路径积分分子动力学（PIMD）等全量子方法高效。

2.3 模拟规模与收敛性验证：尺寸效应至关重要

从补充材料的图S4-S6和S7-S9可以看出，研究者系统地进行了不同超胞尺寸（12×12×12, 16×16×16, 20×20×20）的模拟。这绝非随意为之，而是研究铁电相变时必须考虑的有限尺寸效应。铁电相变与长程的偶极-偶极相互作用密切相关，模拟盒子太小，会人为地增强周期性边界条件带来的镜像相互作用，可能抑制或扭曲铁电畴的形成和演化，导致相变温度、晶格参数乃至相变序的预测出现偏差。

通过比较不同尺寸超胞下晶格参数和Ti-O离子位移随时间演化的曲线（图S4-S9），可以评估模拟结果是否收敛于体材料性质。通常，当超胞尺寸增大到一定程度后，这些宏观观测量的统计平均值和涨落不再发生显著变化，就认为基本消除了尺寸效应。文中采用高达20×20×20的超胞（对应数万个原子），正是为了确保模拟能够捕捉到真实的体相行为，使得基于QTB-MLPMD计算得到的相变压力（图S10-S11中的黄色竖线）和晶格参数更加可靠。这种对模拟收敛性的严谨验证，是确保研究结论物理可信度的基石。

3. 实操流程与关键步骤拆解

3.1 机器学习势函数的构建与训练流程

构建一个适用于钛酸钡相变研究���MLP，是一个系统性的工程，绝非简单调用一个现成工具包就能完成。其核心流程可以拆解为以下几步：

第一步：第一性原理数据集的生成与准备。这是整个工作的“燃料”。数据集的质量和覆盖面直接决定了MLP的上限。针对钛酸钡，我们需要采集覆盖其所有相关相（立方C、四方T、正交O、三方R）以及相变路径中间态的原子构型。采样方法包括：

1. 常温常压下的分子动力学采样：在不同初始相下运行短时间的AIMD，采集动态轨迹中的构型。
2. 势能面扫描：沿着特定的反应坐标（如某个铁电软模的振幅、晶格常数）手动构造一系列静态畸变构型。
3. 主动学习/迭代训练：这是更高效的策略。先用一个初步训练的MLP跑MD，当MLP遇到预测不确定性高（如神经网络委员会预测方差大）的构型时，将这些构型提交给第一性原理计算，并将新数据加入训练集，重新训练MLP。如此循环，直至MLP在感兴趣的温度和压力范围内都保持高置信度。

第二步：原子结构描述符的选择与计算。这是将原子坐标转化为机器学习模型可读输入的“翻译官”。对于钛酸钡这种钙钛矿结构，常用的有平滑重叠原子位置（SOAP）描述符或原子中心对称函数（ACSF）。它们能将每个原子周围的局部化学环境（邻居原子的种类、距离、角度）转化为一个固定长度的、旋转和平移不变的向量。这个步骤需要仔细调节描述符的超参数，如截断半径（应大于相互作用的范围）、高斯函数的宽度等，以确保能充分捕捉到关键的原子间相互作用。

第三步：机器学习模型的训练与验证。将数据集（描述符作为输入，第一性原理计算的总能量、原子力、应力（可选）作为输出）按比例（如8:1:1）划分为训练集、验证集和测试集。常用的模型有神经网络（如Behler-Parrinello型）、高斯过程回归或近年来流行的图神经网络（如MACE、Allegro）。训练时，损失函数通常是能量、力和应力的加权均方误差之和。这里有一个关键技巧：力的权重通常要设得比能量大很多（例如100-1000倍），因为力的数据点数量远多于能量（每个构型有3N个力分量，但只有一个能量），且如前所述，力的精度对MD模拟至关重要。训练过程中要密切关注验证集误差，防止过拟合。

第四步：势函数的部署与MD引擎集成。训练好的模型需要封装成LAMMPS、ASE或i-PI等主流MD程序可以调用的势函数格式。这通常需要编写一个包装器，实现计算给定原子位置下的能量、力和应力的接口。集成后，必须进行严格的基准测试，如图S1-S3所示，在独立的测试集上全面评估其精度。

3.2 量子热浴耦合分子动力学的实现细节

将QTB与基于MLP的MD模拟结合，在技术实现上需要注意以下几个要点：

1. 积分器的选择与参数设置：QTB引入的频率相关噪声要求积分器能够稳定地处理有色随机力。通常采用速度Verlet等时间可逆的积分算法。关键参数是朗之万方程中的摩擦系数γ。γ太小，系统与热浴耦合弱，达到平衡慢；γ太大，动力学会被过度阻尼，可能影响相变动力学。需要针对具体体系进行测试，选择一个能使系统在合理时间内达到平衡且不显著改变系统本征动力学的γ值。

2. 噪声谱的生成与截断：理论上，量子涨落-耗散定理定义的谱密度J(ω)频率范围是0到无穷大。实际计算中，需要设定一个合理的截止频率ω_c。这个频率应高于体系中所有重要的振动模式的频率（可以通过先运行一段经典MD并计算振动谱来估计）。对于钛酸钡，其最高光学支频率通常在几百个波数（cm⁻¹）量级，据此可以设定ω_c。生成有色噪声时，通常在频率空间生成符合指定谱密度的随机数，再通过逆傅里叶变换得到时域噪声。需要确保噪声的随机种子设置正确，以便结果可重复。

3. 模拟流程的编排：一个典型的QTB-MLPMD模拟流程如下：

• 初始化：构建特定相（如立方相）的超胞，根据目标温度初始化原子速度（麦克斯韦-玻尔兹曼分布）。
• 平衡阶段：在NPT或NVT系综下运行足够长的QTB-MD，使体系能量、温度、压力、晶格参数等物理量达到稳定。图S4-S6中时间序列的起始部分就是平衡过程。
• 生产阶段：继续运行QTB-MD，采集用于统计分析的数据。此阶段需要足够长，以获取良好的统计平均并观察可能的相变涨落。
• 观测量的计算：从轨迹中提取晶格常数（a, b, c）、晶胞体积、Ti原子相对于其最近邻氧八面体中心的平均位移（作为铁电极化序参量的代理）、以及它们的统计平均值和标准偏差（如图S10-S12中的μ和σ）。

4. 压力与温度的控制：研究P-T相图需要精确控制压力。通常使用 Parrinello-Rahman 等压控温方法，将目标压力作为输入。通过改变目标压力值，在一系列压力点下进行模拟，从而勾勒出相边界。温度则由QTB热浴本身控制，其设定的温度参数T决定了噪声谱的强度。

3.3 相变识别与数据分析方法

从模拟的海量轨迹数据中准确识别相变点和相结构，是另一个技术关键。文中主要采用了两种互补的分析方法：

1. 晶格参数与序参量监测：这是最直观的方法。如图S10所示，在不同恒定温度下，逐步改变压力进行模拟，记录平衡后的晶格参数a, b, c。当体系发生相变时，晶格对称性改变，会反映在晶格参数的关系上。例如：

• 立方相 (C)：a = b = c。
• 四方相 (T)：a = b ≠ c（通常c轴略长）。
• 正交相 (O)：a ≠ b ≠ c，但特定关系取决于原胞选取。
• 三方相 (R)：a = b = c，但夹角α = β = γ ≠ 90°（在伪立方晶胞描述下，表现为a, b, c有微小差异且伴随原子位移）。 通过观察晶格参数随压力的突变或交叉点，可以初步判断相变压力。同时，如图S11所示，监测Ti离子相对于O离子的位移（δTi-O），在顺电立方相中，这个位移的平均值应为零；在铁电相中，它会沿某个或某几个方向出现非零的平均值，直接指示了极化方向。

2. 统计分析与误差估计：文中提到晶格参数和位移的数值是作为“μ1, μ2, μ3, with an error bar as σ1, σ2, σ3 in Eq. (4)”给出的。这里的μ likely是生产阶段轨迹的时间平均值，而σ是其标准误差（标准差除以平方根观测次数），用于衡量统计的不确定性。在相变点附近，这些序参量的涨落（σ）通常会显著增大，这是相变临界现象的一个特征。因此，误差棒的大小本身也提供了相变区域的信息。

3. 相变压力的确定：图S10和S11中的黄色竖线标识了估计的相变压力。这通常是通过寻找序参量（如c/a比值，或某个方向的δTi-O）随压力变化曲线的拐点，或者观察不同相的自由能（可通过热力学积分或直接比较相共存模拟获得）交叉点来确定的。在MLP-MD模拟中，有时也会直接观察在某个压力下，体系是否会从一种相的自发涨落演变为另一种相。

4. 结果解读与核量子效应量化分析

4.1 压力-温度相图的重构与对比

图S10和S11是整个研究的成果核心，它们直观地展示了结合QTB（考虑NQE）和经典热浴（CTB，忽略NQE）两种情况下，钛酸钡的相图如何随压力变化。

观察图S10（晶格参数）和图S11（Ti-O位移），一个清晰的模式是：在所有研究的温度下（100K-300K），使用QTB（左列子图）与使用CTB（右列子图）所得到的相变压力（黄色竖线位置）存在系统性差异。具体来说，考虑核量子效应（QTB）后，相变压力线整体向高压方向移动。例如，在150K（图S10c,d），从正交相(O)到四方相(T)的转变，QTB模拟给出的相变压力高于CTB模拟。

这背后的物理图像是：核量子效应，特别是轻原子（氧）的零点振动，对晶格有“软化”或“膨胀”效应。零点能使得原子即使在绝对零度也在平衡位置附近振动，这种额外的动能贡献等效于一种“量子压力”，抵抗外部的机械压缩。因此，当考虑NQE时，需要施加更大的外部压力才能迫使晶体从一种相（通常原子位置更不对称）转变为另一种相（通常更对称）。换句话说，NQE稳定了低温下的铁电相，使其在更高的压力下才转变为顺电相。

4.2 体积-压力关系与实验验证

图S12提供了另一个关键验证：晶胞体积随压力的变化。图中将QTB和CTB在不同温度下模拟得到的体积-压力曲线，与实验测得的常温常压体积范围（黄色区域，64.2-64.4 ų）进行了对比。

这里可以解读出两点重要信息：

1. 绝对体积的准确性：在常压附近（P≈0），无论是QTB还是CTB模拟，预测的晶胞体积与实验值吻合得相当好。这首先验证了MLP本身在描述钛酸钡基态和近平衡结构方面的准确性。
2. 压缩行为的差异：随着压力增加，体积减小。对比QTB和CTB的曲线，可以发现考虑NQE（QTB）的体系，在相同压力下体积略大于经典（CTB）体系。这再次印证了量子零点振动对晶格的“支撑”作用，使得量子体系更难被压缩。这种差异在低温下（如100K）应该更为明显，因为此时量子效应占比更大。

实操心得：将模拟结果（如晶格常数、体积）与已知的高质量实验数据在某个点（通常是常温常压）进行对标，是验证整个计算方案（包括MLP和模拟方法）可靠性的黄金标准。如果在这个基准点上都有较大偏差，那么对相变压力等更精细物理量的预测就需要格外谨慎。图S12的这种对比，极大地增强了整个研究结论的可信度。

4.3 超胞尺寸对模拟结果的影响评估

图S4-S9展示了在150K下，不同超胞尺寸（12^3, 16^3, 20^3）的模拟中，晶格参数和Ti-O位移随时间演化的轨迹。仔细分析这些图，我们可以获得关于模拟收敛性的信息：

• 晶格参数（图S4-S6）：对于较小的12^3超胞，晶格参数a, b, c随时间有可见的波动，且三条曲线分离不够清晰（在铁电相中，它们应该分开）。随着超胞增大到16^3和20^3，波动减小，晶格参数趋于稳定，并且不同轴向的参数差异（各向异性）更加明确和稳定地呈现出来。这表明20^3的超胞已经能够较好地代表体材料，有限尺寸效应得到了有效抑制。
• 离子位移（图S7-S9）：Ti-O位移的演化图传达了类似的信息。在小超胞中，位移的涨落较大，且可能由于周期性边界条件强加的虚假关联，使得极化方向不够稳定。在大超胞中，位移序参量表现出更稳定的非零平台（对应铁电相）或零平台（对应顺电相），涨落减小。

结论是：对于钛酸钡这种长程相互作用重要的铁电体，使用足够大的超胞（文中20×20×20）对于获得收敛的、可靠的相变性质和序参量至关重要。基于小超胞的模拟可能会低估相变能垒，或者得到被尺寸效应扭曲的相变动力学。

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

在实际操作QTB-MLPMD模拟研究相变的过程中，会遇到一系列典型问题。以下是一些“踩坑”经验总结：

问题1：MLP在相变点附近预测不准，导致模拟“卡住”或得到错误相。

• 表现：模拟中体系能量异常升高，原子受力突然变大，甚至出现原子“飞散”的非物理情况。
• 原因：训练数据集缺乏相变路径附近（即势能面鞍点区域）的构型样本。MLP在数据稀疏的外推区域预测不确定性高，给出错误能量/力。
• 解决方案：实施主动学习（Active Learning）。具体步骤：a) 用初始MLP运行一段高温MD或施加应变，故意让体系探索构型空间；b) 实时监控MLP的预测不确定性（如委员会模型的方差）；c) 当不确定性超过阈值时，保存当前构型，中断模拟；d) 用第一性原理计算这些“不确定构型”的能量和力；e) 将新数据加入训练集，重新训练并更新MLP；f) 用新MLP继续模拟。如此迭代，直至MLP在目标相变区域表现稳健。

问题2：QTB模拟中，体系温度或能量持续漂移，无法达到稳定。

• 表现：系统总动能或温度偏离设定值，并随时间单调上升或下降。
• 原因：可能是积分时间步长Δt设置过大，导致数值不稳定；或者摩擦系数γ与噪声谱的离散化参数不匹配；也可能是QTB实现中存在数值误差累积。
• 解决方案：首先检查并减小时间步长（如从1 fs减至0.5 fs）。其次，验证有色噪声生成算法的正确性，确保生成的随机力满足预期的自相关函数和谱密度。可以单独测试噪声生成模块，计算其功率谱是否与理论的J(ω)一致。最后，确保QTB与所选系综（NVT/NPT）控温控压方法的兼容性。

问题3：如何确定一个相变是真实的，而不是由模拟参数或有限尺寸效应引起的假象？

• 表现：在某个压力/温度点，序参量发生突变，但重复模拟结果不一致，或者超胞尺寸改变后突变点移动很大。
• 解决方案：进行系统的有限尺寸缩放分析。就像文中做的那样，至少在3个不同尺寸的超胞下重复关键点的模拟。如果相变点（如晶格参数突变对应的压力）随着超胞增大而收敛到一个稳定值，那么可以认为是体相变。此外，可以计算序参量（如极化矢量）的分布函数，在真的一级相变点附近，通常会观察到双峰分布（两相共存），而连续相变则表现为单峰分布的展宽。

问题4：计算资源消耗巨大，如何优化？

• 挑战：QTB-MLPMD模拟，特别是大超胞、长时间模拟，对计算资源要求很高。
• 优化策略：
  1. MLP推理优化：选择计算效率高的MLP架构，如使用线性模型或经过高度优化的神经网络（如使用LAMMPS的pair_style pace或pair_style nnp）。利用GPU加速MLP的能量和力计算。
  2. 并行计算：将MD模拟在数百甚至数千个CPU核心上进行空间分解并行。MLP评估部分也应支持并行。
  3. 采样策略：并非所有压力点都需要从头跑长时模拟。可以采用热力学积分或复本交换（平行回火）方法，更高效地扫描相图。
  4. 数据存储：只保存必要的观测量的时间序列或统计量，而非每一帧的全体原子坐标，以节省I/O和存储空间。

问题5：如何区分核量子效应与其他效应（如非谐性）的影响？

• 深度分析：QTB方法本身已经包含了简谐近似下的量子效应。要更精确地评估非谐性对NQE的贡献，可以将QTB-MD的结果与更精确但昂贵得多的路径积分分子动力学（PIMD）结果进行对比。如果两者在关注的性质（如相变温度偏移）上吻合良好，说明QTB的近似在此体系中是足够的。如果差异显著，则可能需要考虑更高���的方法，如基于PIMD的机器学习势函数模拟，但这会带来巨大的计算负担。在大多数情况下，对于像钛酸钡这样的氧化物，QTB已被证明是捕捉NQE对相变影响的一种有效且高效的近似。