复杂体系热导率计算的新途径：机器学习势函数与噪声抑制技术的协同加速

《Journal of Chemical Theory and Computation》：Accelerating First-Principles Molecular-Dynamics Thermal Conductivity Calculations for Complex Systems

本文系统评述了基于第一性原理分子动力学（MD）的热导率计算新策略，重点探讨了机器学习势函数（MLIPs）与多种噪声抑制技术（如cepstral分析、KUTE）在复杂体系（如InAs纳米线）中的协同应用。作者通过构建高精度MACE势函数，有效解决了传统方法计算成本高、收敛难的问题，为低热导和高热导材料提供了定制化分析方案，显著提升了模拟效率与可靠性。

　　

Accelerating First-Principles Molecular-Dynamics Thermal Conductivity Calculations for Complex Systems

MACE Surrogate Models

MACE是一种消息传递神经网络力场，其模型相对于三维欧几里得群是等变的，这意味着其内部特征和输出能够以标量或张量的形式保持一致地变换。在本研究中，球形谐波的最大阶数被设置为2，因为增加到3只会带来精度的边际提升，但会显著增加计算速度的成本。相互作用的截断半径设置为5 ?，使用了两层具有64个非耦合特征通道的网络，消息传递层的数量设置为2。径向贝塞尔函数的数量设置为8，用于实现平滑截断的多项式阶数设置为5。每层的关联阶数设置为3，用于构建消息的多层感知器的三层宽度设置为64，并使用SiLU作为非线性激活函数。

在训练过程中，总数据集被随机分割，其中90%分配给训练集，10%分配给验证集。训练模型800个周期通常被认为足以实现适当的收敛。在前500个周期或直到损失函数在50次迭代中没有改善为止，训练以较高的力权重和较低的能量权重进行。在随后的300个周期中，能量权重增加，力权重减少。使用了AMSGrad算法，批量大小为10。

Data Set Acquisition and Model Validation

本研究聚焦于图1所示的两个特定InAs纳米线的热传输模拟。它们是六边形形状的锌矿（ZB）[111]和纤锌矿（WZ）[001]纳米线。选择这些结构是因为它们在体相水平上总体相似，且WZ纳米线表面比ZB纳米线光滑得多，这可能导致热导率的巨大差异。

为了创建一个有用且可转移的势函数，使其能够模拟超出密度泛函理论（DFT）可及晶胞尺寸的各种纳米线以及其他InAs结构，在创建机器学习势函数时使用了多样化的数据集。事实上，图1中的特定纳米线并未直接包含在ab initio数据集中，而是用于验证模型。所述数据集包括对应于许多不同表面构型的环境，并包含体相、表面和纳米线结构。训练使用监督式主动学习方法进行，应用了四种不同的策略生成数据：（i）从初始结构进行随机位移；（ii）从初始结构进行随机位移并附加随机双轴应变；（iii）使用模型委员会进行分子动力学模拟，然后选择不确定性最高的结构；（iv）基于模型集合评估的对抗性损失进行优化。最终数据集包含7357个不同的DFT计算结构，每个结构最多包含180个原子。

作为数据集的起点，对所有单个原子坐标应用了从标准差范围在0.01到0.1 ?的高斯分布中抽取的随机位移。对于体相结构和部分纳米线结构，在随机位移之上还施加了应变。对于纳米线，这些应变始终平行于周期性方向，而对于体相结构，整个施加的应变张量根据标准差为0.01的正态分布进行随机化。基于随机位移，总共计算了3115个不同的结构。

分子动力学模拟主要作为补充参考数据来源添加，以实现感兴趣结构的热稳定性，因为当数据库中仅包含少量初始随机位移结构进行训练时，模型会迅速发散到非物理区域。该过程迭代进行：这些分子动力学模拟使用朗之万恒温器在模型早期版本的300 K温度下进行，一旦模型达到改进的热稳定性程度，则在700 K下进行。模拟使用五个独立的MACE模型组成的委员会进行。一旦分子动力学温度达到目标温度的五倍，则认为模拟不稳定并中止。不确定性通过聚合委员会在给定构型中单个力预测的标准差来获得。

进行长时间的分子动力学模拟来采样多样化的参考结构集可能效率低下，进行更传统的不确定性最大化可能更有意义。这可以通过基于对抗性攻击的方法来实现。核心概念是基于最大化对抗性损失，该损失表述为模型集合的力不确定性σ2与在所需温度下评估的玻尔兹曼因子的乘积。该玻尔兹曼因子用于惩罚即使具有高水平不确定性但不物理的结构。在最大化过程中使用了BFGS算法。在每次优化之前，使用标准差为0.05 ?的高斯分布对初始结构的原子位置进行扰动，100个不同的随机种子导致几个独特的局部最小值。基于原子位置之间的绝对差异，丢弃等效结构，并将独特结构附加到数据集，在几次迭代过程中通过此方法总共添加了2097个独特结构。

关于数据集中体相相和纳米结构的混合，起点是通过随机施加位移和双轴应变生成的500个体相ZB和WZ结构的集合（300个无应变和200个应变结构）。通过声子能带结构的优异一致性判断，这足以准确描述两个体相相。接下来添加了一些初始ZB表面构型，从150个随机位移开始。取向选择为最有可能在目标纳米线中作为刻面出现的候选方向：ZB [0 1?1]、[112]和[111]以及WZ [100]。数据集中总共包含1332个不同的表面板结构。

数据集中最大的部分，包含总共5525个不同的结构，由薄纳米线结构组成，这些结构主要需要确保复杂表面环境的热稳定性。主要组成部分是基于ZB的[111]取向纳米线，具有六边形、圆形和正方形形状，具有[0 1?1]和[112]刻面。最后，还包含了包含表面缺陷的表面和纳米线结构，以进一步增加数据集的多样性。

图2b显示了基于整个数据集的 parity 图。总体一致性非常好，力的均方根误差（RMSE）为59 meV ?^–1，清楚地表明模型能够学习不同纳米线结构的不同构型。

振动特性也表现出优异的一致性，使用phonopy中实现的有限差分方法获得，位移距离为0.01 ?。图2c显示了从DFT和MACE模型评估的WZ体相相的声子能带结构比较。其一致性质量与ZB相相似，这证实了非体相结构的添加不会降低体相性能。

为了展示目标纳米线的性能，在图2d中，我们比较了每个纳米线100个随机位移结构上的力误差。这些结构中没有一个直接包含在参考数据集中。尽管如此，RMSE值与训练集的值范围相似，ZB和WZ纳米线分别为60和36 meV ?^–1。这展示了该模型对相关结构的可转移性。

Ab-Initio Data Generation Procedure

DFT计算使用维也纳Ab-initio模拟软件包（VASP）进行，采用投影增强波（PAW）方法，使用针对固体修订的Perdew-Burke-Ernzerhof（PBEsol）交换相关泛函，并考虑In的s、p和d作为价态，As仅考虑s和p作为价态。平面波能量截断设置为450 eV，高斯展宽宽度设置为0.05 eV。对于锌矿结构，k点网格设置为常规立方晶胞的5 × 5 × 5，所有其他结构的k点网格根据各自的晶格参数进行调整以符合等效密度，总是向上取整。能量截断和k点密度是基于仔细的收敛测试选择的，以实现1 meV/原子的精度。自洽过程的能量收敛标准设置为10^–8eV，结构优化的收敛标准设置为最大力为10^–3eV ?^–1。投影算子在倒易空间中评估以提高精度。

Green–Kubo Thermal-Conductivity Calculations

基于涨落-耗散定理，GK方法从平衡分子动力学模拟期间热流J的涨落中提取热导率。具体地，热导率κ由热流自相关函数（HFACF）的时间积分给出。

对于质量输运可忽略的固体（如离子扩散），可以合理地排除热流中的对流分量，以显著减少HFACF中的噪声，同时不影响产生的热导率。应该强调，这对于晶体也是一种近似，并且与考虑全热流相比会导致微小差异。更严格的方法将涉及考虑相对于平衡位置的有界位移。然而，我们将专注于一种更普遍适用的减少噪声的方法，即利用热流的规范不变性来移除不贡献于热导率的部分。这是通过从总热流中减去原子维里S_i的时间平均的贡献来实现的。

对于准一维亚周期系统，晶格热导率可以视为一维标量，只有沿周期轴的热流J的投影与计算相关。同样，方程3中的外积简化为该投影在两个不同时间评估的普通乘积。因此，当只有一个周期方向时，所需“展开”的额外计算工作量相当小，因为额外原子的数量有限。

我们将参考文献19中创建的热流实现适配到MACE模型，并使用方程4的直接评估验证了我们的实现。为了对MACE进行基于势流的标准固定，我们重用了为热流计算的量。

所有GK模拟的分子动力学运行使用大规模原子/分子并行模拟器（LAMMPS）进行，而用于构建数据集的分子动力学模拟通过原子模拟环境（ASE）中可用的积分器进行。在每种情况下，使用1 fs的积分时间步长，并且在每个生产运行之前，在NVT系综中至少有20 ps的平衡期。在平衡期之后的进一步准备步骤中，通过缩放原子力同时保持系统的总能量，将纳米线的角动量强制为零。在所有模拟中，我们专注于300 K的温度。GK模拟中纳米线的长度为42.2 ?，为了获得体积，我们基于能量优化结构的元素特异性范德华球计算其所有原子周围的非重叠体积。

Cepstral Analysis

Cepstral分析是一种评估Green-Kubo积分的有效技术。它基于热流的功率谱S(ν)，该功率谱由热流自相关函数的傅里叶变换给出。除了预因子，频率ν = 0处的功率谱与热导率方程3具有相同的形式。因此，重点在于分析频谱的低频部分。需要选择一个合适的截止频率，对频谱进行重新采样，并形成其对数逆离散傅里叶变换以创建所谓的“倒谱”。

该倒谱用于获得S(ν = 0)的值，这一步我们通过仅使用前几个倒谱系数C_n来有效地应用低通滤波器。最优倒谱系数数量P*通过最小化二阶Akaike信息准则（AIC_c）客观地选择。

为了在AIC_c存在多个局部最小值的情况下提高方法的一致性，我们额外建议使用模型平均，这是模型选择中常用的技术。为此，我们对每个模型（每个倒谱系数数量）赋予权重w_i。

然后使用模型平均（MA）从cepstral分析预测总热导率κ_cepstral,MA。应该注意，使用模型平均时的不确定性估计总是高于不使用它的情况。

Uncertainty-Based Analysis Based on KUTE

KUTE是Otero-Lema等人最近提出的一种基于不确定性传播直接评估GK积分的方法。为了描述该方法，我们首先从单个分子动力学运行中引入HFACF元素（R_k）的离散公式，这是数值积分方程3所必需的。

为了进行统计分析，整个热流轨迹被折叠成M个等长的片段，并计算HFACF的平均元素。

使用梯形方案数值获得HFACF的累积积分。

为了评估该表达式的不确定性，KUTE忽略了R?_k之间的相关性，并执行标准传播。

这种传播的不确定性然后用于创建积分在相关时间结束时的加权平均值。

这导致了一次GK运行的电导率和不确定性。为了在KUTE方法中获得最终值，进行了几次独立模拟，然后使用来自方程22的不确定性再次使用加权平均值评估各个独立模拟的热导率。

KUTE方法的一个关键近似是在传播电导率积分的不确定性时忽略了协方差。然而，相邻的J_i+k^(A)J_i^(A)值将是高度相关的，这也会在R?_k之间引入相关性。因此，我们提出了一个扩展，基于M个片段的平均值，为HFACF的元素引入完整的协方差矩阵。

并且我们使用更简单的欧拉积分方法来替换方程18。

方程27的不确定性然后被计算，不忽略相关性。

Results and Discussion

Low-Conductivity ZB-Phase Nanowire

正如引言中提到的，已经开发了多种方法来有效地处理噪声相关的函数，以在可接受的时间尺度内评估热导率。为了演示传统应用的方法，我们在图3a,b中分别显示了ZB纳米线系统的热流自相关函数及其积分。这里，相关函数在11次独立运行上取平均，每次模拟时间为1 ns，相关时间为0.5 ns。热流每5个时间步（每个时间步为1 fs）评估一次以节省计算时间，因为我们的测试表明，对于本工作中考虑的系统，较低的采样率不会影响结果。从图3可以清楚地看出，即使经过大量的模拟时间平均，热导率基于相关长度仍然显示出很大的波动。

在没有平均的情况下，波动变得更加明显，使得难以估计热导率。在文献中，这通常通过手动选择主要平台区域来解决。这种选择基于这样一个事实：经过一定