DFT+机器学习势函数精准预测材料热导率：以TaFeSb缺陷工程为例

1. 项目概述：当机器学习遇见缺陷工程，如何精准预测材料热导率？

在热电材料研发的前沿，有一个长期困扰理论家和实验家的“经典难题”：基于完美晶体模型的第一性原理计算，预测出的晶格热导率（κL）常常远高于实验测量值。以明星热电材料半赫斯勒（Half-Heusler, HH）化合物TaFeSb为例，理论预测的室温κL值在18到30 Wm⁻¹K⁻¹之间，而实验测得的数值却只有9 Wm⁻¹K⁻¹左右。这个巨大的鸿沟，不仅让理论模型的可靠性受到质疑，更阻碍了我们基于计算进行高效材料设计的步伐。

问题的根源，往往隐藏在那些“不完美”之中。在实际合成的材料样品里，本征点缺陷——比如原子空位、间隙原子或原子错位——几乎不可避免。这些微观尺度的“瑕疵”，传统上被认为是性能的杀手，但在热电领域，它们却可能是提升性能的“功臣”。缺陷能强烈散射传热的声子，从而有效降低晶格热导率，这正是提升热电优值ZT的关键。然而，精确计算含缺陷大体系的声子谱和热导率，对计算资源的要求是惊人的。传统的密度泛函微扰理论（DFPT）或有限位移法，在面对包含缺陷、需要构建数十甚至上百个原子的大超胞时，计算成本变得难以承受。

这正是我们这项工作的切入点。我们不再将缺陷视为需要规避的麻烦，而是将其作为核心的设计变量。我们采用了一种创新的“DFT+机器学习”混合计算策略，以TaFeSb为模型体系，系统探究了本征缺陷对其电子结构、声子性质乃至最终晶格热导率的微观影响机制。简单来说，我们先用高精度的DFT计算训练出一个能够“学会”原子间相互作用的机器学习势函数（MLIP），再用这个势函数去驱动大规模的分子动力学模拟，高效计算含缺陷体系的声子和热导率。这套方法的核心价值在于，它既保持了接近第一性原理的精度，又将计算规模扩展到了传统方法难以企及的程度，为在复杂材料中开展系统的缺陷工程研究，打开了一扇新的大门。

2. 核心思路与技术路线拆解：为什么是“DFT+MLIP”？

要理解我们方案的优势，首先得看清传统方法的瓶颈。研究缺陷对声子输运的影响，理想流程是：构建包含缺陷的大超胞 -> 计算体系的二阶（谐波）和三阶（非谐）力常数 -> 求解玻尔兹曼输运方程（BTE）得到κL。其中，力常数的计算是核心，也是最耗资源的一步。

2.1 传统方法的“不可能三角”：精度、尺度与成本

在材料计算领域，长期存在一个“不可能三角”：计算精度、系统尺度和计算成本，三者难以兼得。

• 纯DFT（如DFPT）：精度最高，是金标准。但对于包含缺陷的3x3x3超胞（约135个原子），计算三阶力常数所需的微扰数量巨大，对计算资源和时间的要求是天文数字，几乎不可行。
• 经典分子动力学（MD）：可以轻松处理数千甚至上万个原子，尺度最大。但问题在于精度：常用的经验势（如Lennard-Jones）难以准确描述半赫斯勒化合物中复杂的金属-共价键合作用，计算结果可靠性存疑。
• 纯机器学习势（MLIP）训练：理论上可以兼顾精度和尺度。但如果训练数据完全来自含缺陷的大超胞DFT计算，那么数据生成本身就会陷入DFT的成本泥潭。

我们的策略，正是为了打破这个三角。

2.2 我们的“DFT引导的MLIP”策略：分而治之，高效融合

我们设计的计算框架是一个精巧的“分阶段”流程，核心思想是用高精度但小尺度的DFT计算，来“教导”一个可以用于大尺度模拟的MLIP模型。

第一阶段：高精度DFT计算奠定基础

1. 缺陷形成能与电子结构分析：对完美的TaFeSb原胞以及包含各种可能缺陷（空位、反位、间隙原子、原子交换）的3x3x3超胞，进行严格的DFT计算。这一步虽然涉及超胞，但只是单点能计算，成本尚可接受。我们通过形成能计算，筛选出最可能稳定存在的缺陷类型（Fe间隙原子IFe和SbTa反位缺陷），并分析它们对电子能带结构、态密度的影响，特别是带隙的变化。
2. 生成训练数据：我们不直接对含缺陷的大超胞做声子计算，而是用DFT进行从头算分子动力学（AIMD）模拟。模拟对象包括完美的2x2x2超胞，以及含有低浓度缺陷的类似体系。在300K到1100K的不同温度下运行AIMD，采集体系演化过程中数千个原子构型（快照）。对于每一个构型，用DFT精确计算其总能量、每个原子所受的力以及体系的应力张量。这些（构型，能量，力，应力）数据对，就构成了训练机器学习势的“教科书”。

注意：训练数据的质量是关键。AIMD的温度范围要足够宽，以覆盖材料在应用温区附近的原子振动；构型采样要足够充分，以确保MLIP能学习到各种可能的原子局部环境。我们使用了约3000个构型，这是一个在精度和成本间取得平衡的经验值。

第二阶段：机器学习势的构建与验证我们选择矩张量势（Moment Tensor Potential, MTP）作为我们的MLIP框架。MTP将体系的总能量表达为每个原子在其截断半径内局部环境贡献的求和。通过优化一组参数，使MTP对训练集构型预测的能量、力和应力，与DFT计算值的差异（即损失函数）最小化。

• 训练目标函数：最小化 Σ [we*(E_DFT - E_MTP)² + wf*Σ|F_DFT - F_MTP|² + ws*Σ|σ_DFT - σ_MTP|²]这里，we, wf, ws是分配给能量、力和应力的权重。我们将wf和ws设为0.1，we设为1，这表示我们更看重能量的精确拟合，这在后续计算结合能、形成能时至关重要。
• 势函数验证：训练完成后，必须验证MTP的可靠性。我们用训练好的MTP和纯DFT（DFPT和有限位移法）分别计算完美TaFeSb的声子色散关系和态密度。当两者结果高度吻合时（如图4a,b所示），我们才有信心用MTP进行后续更昂贵的计算。

第三阶段：大尺度声子与热导率计算一旦MTP通过验证，它就从“学生”变成了“工具”。

1. 计算力常数：对包含目标缺陷（如IFe, SbTa）的大超胞，我们使用基于MTP的经典分子动力学进行长时间模拟。从模拟轨迹中，通过统计方法提取出体系的二阶和三阶力常数。这一步完全绕开了DFT，计算速度提升了数个数量级。
2. 求解热导率：将得到的力常数输入到求解玻尔兹曼输运方程（BTE）的软件包（如KALDO）中，在弛豫时间近似（RTA）下，计算从低温到高温的晶格热导率κL。

这套方法的核心优势在于，它将DFT“贵而精”的计算能力，浓缩到了一个可以快速评估的MLIP模型中，从而实现了对含缺陷大体系声子性质的“既快又准”的研究。这不仅仅是计算工具的组合，更是一种研究范式的转变：从回避缺陷的理想化计算，转向直面并利用缺陷的真实材料模拟。

3. 缺陷的筛选与形成机制：谁才是“最可能出现的客人”？

在材料生长过程中，哪种缺陷会自发形成，是由其形成能决定的。形成能越低，缺陷在热力学上越稳定，出现的概率就越高。我们的计算首先聚焦于此，为后续研究锁定目标。

3.1 缺陷类型与形成能计算

我们系统考察了四类本征点缺陷（如图1示意）：

1. 空位（Vα）：从α原子位点移除一个原子。
2. 反位缺陷（βα）：用β类型的原子替换α位点上的原子。
3. 间隙缺陷（Iα）：在HH晶体结构固有的空位（4d位置）添加一个α原子。
4. 原子交换（β⇔α）：交换两个不同位点上的α和β原子。

形成能计算公式如下：

• 空位：ΔE(Vα) = E_defect - E_perfect + μ_α
• 反位缺陷：ΔE(βα) = E_defect - E_perfect + μ_α - μ_β
• 间隙缺陷：ΔE(Iα) = E_defect - E_perfect - μ_α
• 原子交换：ΔE(β⇔α) = E_defect - E_perfect

其中，E是超胞总能量，μ是相应元素在其稳定固态相中的化学势（我们采用Materials Project数据库中的值：Ta和Fe为体心立方，Sb为三角晶系）。

3.2 计算结果与物理意义

计算得到的形成能（图2）揭示了一个清晰的趋势：

• 最稳定的缺陷：Fe间隙原子（IFe），形成能最低，约为1.28 eV。这意味着在TaFeSb的生长过程中，额外的Fe原子占据晶体结构中本来的空位，是能量上最有利的过程。
• 次稳定缺陷：SbTa反位缺陷，形成能约为1.50 eV，略高于IFe。即一个Sb原子占据了本应是Ta原子的位置。
• 其他缺陷：如TaSb反位、FeTa反位等，形成能显著更高（>2 eV），在热平衡条件下浓度会很低。

这个结果与半赫斯勒体系的化学特性高度吻合。HH化合物的通式为XYZ，其晶体结构可以看作三个面心立方（fcc）子晶格相互穿插而成，并天然存在一个空位。这个空位为小原子（如Fe）的填入提供了空间。同时，Sb和Ta的电负性、原子半径差异可能导致一定程度的原子互换倾向。

实操心得：化学势的选择。形成能计算严重依赖于化学势μ的取值。在实际材料生长中（如熔体生长、溅射），原子的化学势并非孤立固态的化学势，而是受生长环境（富Ta、富Fe或富Sb）影响的变量。我们的计算假设了生长条件接近热力学平衡，且各元素化学势由其单质标准态决定。若需精确预测特定工艺下的缺陷浓度，需要构建更复杂的相图并确定在该生长条件下的实际化学势，这是一个重要的深化方向。

4. 缺陷对电子结构的重构：带隙收缩与杂质态引入

明确了最可能存在的缺陷后，我们紧接着分析它们如何改变材料的“电子骨架”。电子结构是决定材料电导率、塞贝克系数等电输运性质的基础。

4.1 完美TaFeSb的电子结构基线

完美的TaFeSb是一种窄带隙半导体，我们的DFT-PBE计算给出的带隙约为0.85 eV。需要注意的是，DFT普遍存在带隙低估的问题，但有趣的是，对于许多HH化合物，理论值反而高于实验值。例如，实验通过红外光谱测得TaFeSb的带隙约为0.53 eV，而更高级的杂化泛函计算甚至给出~0.9 eV的更大值。这暗示，除了交换关联泛函的局限性，材料本身的缺陷可能是导致差异的重要原因。

4.2 缺陷诱导的电子态局域化

引入3.7%浓度的IFe和SbTa缺陷后，电子态密度（DOS）和能带结构发生了显著变化（图3）：

• Fe间隙缺陷（IFe）：在费米能级下方、非常靠近价带顶的位置，引入了一个狭窄的杂质能带。该能带的态密度投影分析显示，它主要来源于间隙Fe原子的d电子态。这个杂质能带就像在禁带中打入了一个“钉子”，将理论带隙从0.85 eV显著降低至0.40 eV。
• SbTa反位缺陷：同样引入了一个杂质能带，但其位置在费米能级上方、导带底附近。这个能带的电子态更加离域，不仅来源于占据Ta位的Sb原子，也涉及其周围的Fe原子。其效果是将带隙减小至0.80 eV。

4.3 对热电性能的潜在影响

这种电子结构的改变具有双重效应：

1. 带隙减小：使理论值更接近实验测量值（0.53 eV），部分解释了理论与实验的差异。
2. 引入杂质态：在带隙中引入的局域化或半局域化态，会成为额外的载流子散射中心，可能影响载流子迁移率。同时，它们也可能提供额外的载流子输运通道。具体是提升还是降低电导率，强烈依赖于费米能级的位置（即掺杂水平）。在我们的中性缺陷计算中（未额外掺杂），费米能级钉扎在杂质能带附近，这通常会增加电子热导率，并可能降低塞贝克系数，对热电优值ZT的影响需要综合权衡。

注意事项：超胞折叠与能带展开。计算含缺陷超胞的电子结构时，由于平移对称性被破坏，传统的原胞布里渊区不再适用，计算出的能带会在超胞的小布里渊区内折叠，显得非常复杂。为了与完美晶体的能带对比，我们使用了BandUP等代码进行“能带展开”，将超胞计算的结果映射回原始原胞的布里渊区，从而清晰地观察缺陷引起的能带扰动和杂质态的出现。这是分析缺陷电子结构的关键后处理步骤。

5. 声子谱的扰动：局域模与声子工程

热电材料中，声子（晶格振动）是热量的主要载体。降低晶格热导率的核心，就在于增加对声子的散射。缺陷正是最有效的声子散射源之一。我们通过训练好的MTP势函数，计算了含缺陷体系的声子态密度，揭示了缺陷影响声子输运的微观机制。

5.1 完美TaFeSb的声子谱特征

完美TaFeSb的声子态密度显示，在约7 THz频率附近，存在一个明显的“声学支-光学支”带隙，宽度约0.95 THz（图4a,b）。这个带隙的存在，意味着在这个频率范围内的声子模式是禁止的。声学支声子（低频）主要负责长波传热，而高频光学支声子对热导的贡献相对较小。

5.2 缺陷诱导的局域化振动模式

引入缺陷后，声子谱发生了关键性变化：

• Fe间隙缺陷（IFe）：在约6.5 THz处，产生了一系列新的、尖锐的声子态，正好位于原本的带隙之中，并将带隙宽度压缩至0.68 THz（图4c）。通过分波态密度（PDOS）分析发现，这些新出现的态几乎完全局域在间隙Fe原子上，周围其他原子的贡献微乎其微（图4e）。这是一种典型的局域振动模。
• SbTa反位缺陷：同样减小了带隙（至0.74 THz），但其产生的新声子态主要出现在高于7.25 THz的光学支区域（图4d）。PDOS分析表明，这些态主要局域在被替换原子（Sb）周围的Fe原子壳层上，而非缺陷原子本身（图4f）。

5.3 局域模如何散射声子？

这些缺陷引入的局域振动模式，是降低热导率的“功臣”，其物理机制主要体现在两方面：

1. 共振散射：当传播的声子（行波）的频率与缺陷局域模的频率接近时，会发生强烈的共振散射。局域模就像一个“陷阱”，吸收声子的能量，极大地缩短了声子的寿命和平均自由程。IFe缺陷引入的局域模正好位于声学支顶部，这个频率区域的声子对热导贡献很大，因此散射效果尤为显著。
2. 破坏周期性，增强非谐性：缺陷破坏了晶体完美的平移对称性，使得原子间的力常数发生局域扰动。这不仅增加了声子-缺陷散射（与频率的平方成正比，对高频声子散射更强），还可能增强三阶非谐力常数，从而加剧声子-声子Umklapp散射，这是高温下热阻的主要来源。

实操心得：MTP对声子计算精度的验证。图4(a,b)的对比至关重要。它显示了用MTP+有限位移法计算的声子态密度，与纯DFT（DFPT和有限位移法）结果高度一致。这直接证明了我们训练的MTP势函数在描述原子间相互作用（特别是二阶力常数）上具有DFT级别的精度。只有在完成这一步验证后，我们才有足够的底气使用MTP去计算更复杂的含缺陷体系，以及用于提取三阶力常数的大规模MD模拟。这是整个工作可信度的基石。

6. 晶格热导率的定量计算与实验吻合

所有分析的最终落脚点，是定量预测缺陷对宏观可观测物理量——晶格热导率κL的影响。我们通过求解玻尔兹曼输运方程（BTE）实现了这一目标。

6.1 计算方法与参数设置

我们使用KALDO软件包，在弛豫时间近似（RTA）下求解BTE。所需的输入是体系的二阶（谐波）和三阶（非谐）力常数。

• 力常数获取：对于完美和含缺陷体系，我们使用基于训练好的MTP进行长时间（纳秒级）的NPT系综分子动力学模拟。从平衡后的轨迹中，采用相关函数法提取出所有原子对之间的二阶和三阶力常数。为了对比，完美TaFeSb的三阶力常数也使用了更精确但昂贵的DFT+DFPT方法（通过D3Q代码）计算。
• 计算细节：采用足够密的q点网格（如24x24x24）对布里渊区进行采样，以确保声子谱和群速度的收敛。计算温度范围从100K到800K，覆盖了热电应用的主要温区。

6.2 计算结果分析与讨论

图5展示了我们的核心计算结果，可以清晰地看到几个关键现象：

1. 完美晶体理论值高估：无论是基于DFT+DFPT（最精确）还是MTP+MD计算，完美TaFeSb的理论κL曲线在室温附近都比实验值高出约2-4倍。这直接印证了引言中提到的理论与实验的巨大差距。MTP+MD的结果比DFT+DFPT还要高约20%，说明我们的MLIP在非谐相互作用上可能存在轻微的高估，但这不影响定性趋势。
2. 缺陷显著降低热导率：引入3.7%浓度的IFe或SbTa缺陷后，计算出的κL曲线大幅下降，向实验数据靠拢。其中，IFe缺陷的降低效果比SbTa更明显。这与声子谱的分析一致：IFe引入的局域模位于声学支顶部，对主要热载流子（低频声学声子）的散射更强。
3. 浓度效应：我们计算了更高缺陷浓度（12.5%）的情况。结果显示κL进一步降低，甚至略低于部分实验点。结合缺陷形成能均大于1 eV的事实，可以推断实际样品中的缺陷浓度可能接近我们计算的较低值（~3.7%），并且可能是多种缺陷共存的混合状态。
4. 与实验的吻合：在400K至800K的中温区域，考虑缺陷后的理论计算曲线与实验数据吻合得非常好。这一温区是热电发电最具应用价值的范围。在更低温度（<300K），理论值仍高于实验，这可能是因为RTA近似在低温下失效（此时Normal散射过程占主导）。在更高温度，多晶样品的晶界散射、可能存在的极化子效应等其他机制开始凸显，这些在我们的单晶模型中没有考虑。

6.3 对热电优值ZT的启示

晶格热导率κL的降低，直接有利于提升热电优值ZT = (S²σ/κ)T，其中S是塞贝克系数，σ是电导率，κ是总热导率（κ = κ_e + κ_L）。我们的工作表明，在TaFeSb中，本征的Fe间隙和SbTa反位缺陷，通过强烈的声子散射，能将κL降低至接近实验值的水平。这意味着，实际测得的高ZT值，部分正是源于这些“不可避免”的缺陷。这为材料设计提供了新思路：与其追求难以实现的完美晶体，不如主动引入或调控这类低形成能的特定缺陷，作为优化热电性能的“预设计”手段。

常见问题与排查：

• 计算的热导率仍然偏高？首先检查力常数的收敛性，包括截断半径、MD模拟时长和采样精度。其次，确认BTE计算中q点网格是否足够密。最后，考虑是否忽略了其他散射机制，如晶界散射（对多晶样品至关重要）、同位素散射等。可以尝试在BTE中显式加入点缺陷散射率进行计算。
• MTP训练误差大，导致声子谱有虚频？虚频（负频率）的出现通常意味着势函数在某个原子构型下预测的力常数矩阵不正定，即结构不稳定。这往往是因为训练数据未能充分覆盖该缺陷构型附近的势能面。解决方法是在AIMD采样时，针对可能产生虚频的缺陷区域进行增强采样，或在该区域手动构建微扰构型加入训练集。
• 如何将这种方法应用于其他材料？本框架具有普适性。对于任何新材料，步骤是：1) DFT计算完美与小缺陷超胞的电子结构；2) 基于完美和微扰结构进行AIMD采样生成训练数据；3) 训练并验证该材料的MLIP；4) 用MLIP-MD研究缺陷效应。关键在于步骤2)的采样要充分，且MLIP的验证（如声子、弹性常数）必须严格。

7. 方法论的普适性、局限性与未来展望

我们的研究不仅揭示了TaFeSb中特定缺陷的物理效应，更验证了“DFT+MLIP”这一计算范式在研究复杂材料缺陷问题上的强大能力。

7.1 “DFT+MLIP”范式的优势与挑战

优势：

• 精度与效率的平衡：在保持接近DFT精度的前提下，将声子及热导率计算的能力从数十原子提升到数百甚至上千原子体系，使得研究真实缺陷浓度、缺陷团簇乃至晶界成为可能。
• 可扩展性：一旦MLIP训练完成，它可以被用于各种基于MD的模拟，如扩散系数计算、力学性能测试、高温相变研究等，性价比极高。
• 高通量筛选潜力：结合主动学习算法，可以自动化地生成训练数据并迭代改进势函数，为大规模筛选有潜力的热电材料缺陷构型提供了工具。

挑战与注意事项：

• 训练数据的质量与数量：MLIP的精度完全依赖于训练数据。对于含有强电子关联性或复杂化学键的材料，生成具有代表性的训练数据集本身可能是一项挑战。需要精心设计AIMD的初始条件和采样策略。
• 外推风险：MLIP在训练数据所覆盖的原子构型和能量范围内是可靠的，但对于远离训练集的极端构型（如高缺陷浓度、表面、断裂），其预测可能不准确。使用时需明确其适用范围。
• 计算流程的复杂性：整个流程涉及DFT计算、AIMD模拟、MLIP训练、力常数提取和BTE求解等多个环节，需要熟练使用多种计算软件（如Quantum ESPRESSO, LAMMPS, phonopy, KALDO等），并对每一步的物理含义和参数设置有深刻理解。

7.2 未来研究方向

基于本工作，有几个方向值得深入探索：

1. 多缺陷协同效应：实际材料中往往多种缺陷共存。可以研究Fe间隙和SbTa反位缺陷之间的相互作用，以及它们与掺杂原子的协同效应，如何进一步调控声子谱和热导率。
2. 有限温度与动力学效应：我们的形成能计算是0K下的静态值。在实际生长温度下，缺陷的形成熵和振动自由能也会影响其平衡浓度。可以结合MLIP进行高温MD模拟，直接观察缺陷的动态形成与湮灭过程。
3. 超越弛豫时间近似（RTA）：RTA是BTE的一种简化，在低温或某些材料中可能不准确。未来可以使用MLIP结合更精确的迭代法求解BTE，或直接使用非平衡分子动力学（NEMD）计算热导率，进行交叉验证。
4. 面向工程的材料设计：将这套计算框架集成到材料设计平台中。针对目标热电性能（如高ZT值），逆向设计最佳的缺陷类型、浓度和分布，为实验合成提供明确的理论指导。

回过头看，这项工作的最大价值在于它提供了一个强有力的案例：通过融合第一性原理的精度与机器学习的效率，我们能够以前所未有的清晰度，窥见微观缺陷如何左右宏观性能。在热电材料乃至更广阔的功���材料领域，这种“从缺陷中要性能”的思路，辅以先进的计算手段，正引领我们从经验摸索走向理性设计。对于计算材料学的研究者而言，掌握并发展此类多尺度模拟方法，无疑是应对未来复杂材料体系挑战的关键技能。