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

1. 项目概述当机器学习遇见缺陷工程如何精准预测材料热导率在热电材料研发的前沿有一个长期困扰理论家和实验家的“经典难题”基于完美晶体模型的第一性原理计算预测出的晶格热导率κL常常远高于实验测量值。以明星热电材料半赫斯勒Half-Heusler, HH化合物TaFeSb为例理论预测的室温κL值在18到30 Wm⁻¹K⁻¹之间而实验测得的数值却只有9 Wm⁻¹K⁻¹左右。这个巨大的鸿沟不仅让理论模型的可靠性受到质疑更阻碍了我们基于计算进行高效材料设计的步伐。问题的根源往往隐藏在那些“不完美”之中。在实际合成的材料样品里本征点缺陷——比如原子空位、间隙原子或原子错位——几乎不可避免。这些微观尺度的“瑕疵”传统上被认为是性能的杀手但在热电领域它们却可能是提升性能的“功臣”。缺陷能强烈散射传热的声子从而有效降低晶格热导率这正是提升热电优值ZT的关键。然而精确计算含缺陷大体系的声子谱和热导率对计算资源的要求是惊人的。传统的密度泛函微扰理论DFPT或有限位移法在面对包含缺陷、需要构建数十甚至上百个原子的大超胞时计算成本变得难以承受。这正是我们这项工作的切入点。我们不再将缺陷视为需要规避的麻烦而是将其作为核心的设计变量。我们采用了一种创新的“DFT机器学习”混合计算策略以TaFeSb为模型体系系统探究了本征缺陷对其电子结构、声子性质乃至最终晶格热导率的微观影响机制。简单来说我们先用高精度的DFT计算训练出一个能够“学会”原子间相互作用的机器学习势函数MLIP再用这个势函数去驱动大规模的分子动力学模拟高效计算含缺陷体系的声子和热导率。这套方法的核心价值在于它既保持了接近第一性原理的精度又将计算规模扩展到了传统方法难以企及的程度为在复杂材料中开展系统的缺陷工程研究打开了一扇新的大门。2. 核心思路与技术路线拆解为什么是“DFTMLIP”要理解我们方案的优势首先得看清传统方法的瓶颈。研究缺陷对声子输运的影响理想流程是构建包含缺陷的大超胞 - 计算体系的二阶谐波和三阶非谐力常数 - 求解玻尔兹曼输运方程BTE得到κL。其中力常数的计算是核心也是最耗资源的一步。2.1 传统方法的“不可能三角”精度、尺度与成本在材料计算领域长期存在一个“不可能三角”计算精度、系统尺度和计算成本三者难以兼得。纯DFT如DFPT精度最高是金标准。但对于包含缺陷的3x3x3超胞约135个原子计算三阶力常数所需的微扰数量巨大对计算资源和时间的要求是天文数字几乎不可行。经典分子动力学MD可以轻松处理数千甚至上万个原子尺度最大。但问题在于精度常用的经验势如Lennard-Jones难以准确描述半赫斯勒化合物中复杂的金属-共价键合作用计算结果可靠性存疑。纯机器学习势MLIP训练理论上可以兼顾精度和尺度。但如果训练数据完全来自含缺陷的大超胞DFT计算那么数据生成本身就会陷入DFT的成本泥潭。我们的策略正是为了打破这个三角。2.2 我们的“DFT引导的MLIP”策略分而治之高效融合我们设计的计算框架是一个精巧的“分阶段”流程核心思想是用高精度但小尺度的DFT计算来“教导”一个可以用于大尺度模拟的MLIP模型。第一阶段高精度DFT计算奠定基础缺陷形成能与电子结构分析对完美的TaFeSb原胞以及包含各种可能缺陷空位、反位、间隙原子、原子交换的3x3x3超胞进行严格的DFT计算。这一步虽然涉及超胞但只是单点能计算成本尚可接受。我们通过形成能计算筛选出最可能稳定存在的缺陷类型Fe间隙原子IFe和SbTa反位缺陷并分析它们对电子能带结构、态密度的影响特别是带隙的变化。生成训练数据我们不直接对含缺陷的大超胞做声子计算而是用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.1we设为1这表示我们更看重能量的精确拟合这在后续计算结合能、形成能时至关重要。势函数验证训练完成后必须验证MTP的可靠性。我们用训练好的MTP和纯DFTDFPT和有限位移法分别计算完美TaFeSb的声子色散关系和态密度。当两者结果高度吻合时如图4a,b所示我们才有信心用MTP进行后续更昂贵的计算。第三阶段大尺度声子与热导率计算一旦MTP通过验证它就从“学生”变成了“工具”。计算力常数对包含目标缺陷如IFe, SbTa的大超胞我们使用基于MTP的经典分子动力学进行长时间模拟。从模拟轨迹中通过统计方法提取出体系的二阶和三阶力常数。这一步完全绕开了DFT计算速度提升了数个数量级。求解热导率将得到的力常数输入到求解玻尔兹曼输运方程BTE的软件包如KALDO中在弛豫时间近似RTA下计算从低温到高温的晶格热导率κL。这套方法的核心优势在于它将DFT“贵而精”的计算能力浓缩到了一个可以快速评估的MLIP模型中从而实现了对含缺陷大体系声子性质的“既快又准”的研究。这不仅仅是计算工具的组合更是一种研究范式的转变从回避缺陷的理想化计算转向直面并利用缺陷的真实材料模拟。3. 缺陷的筛选与形成机制谁才是“最可能出现的客人”在材料生长过程中哪种缺陷会自发形成是由其形成能决定的。形成能越低缺陷在热力学上越稳定出现的概率就越高。我们的计算首先聚焦于此为后续研究锁定目标。3.1 缺陷类型与形成能计算我们系统考察了四类本征点缺陷如图1示意空位Vα从α原子位点移除一个原子。反位缺陷βα用β类型的原子替换α位点上的原子。间隙缺陷Iα在HH晶体结构固有的空位4d位置添加一个α原子。原子交换β⇔α交换两个不同位点上的α和β原子。形成能计算公式如下空位Δ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和能带结构发生了显著变化图3Fe间隙缺陷IFe在费米能级下方、非常靠近价带顶的位置引入了一个狭窄的杂质能带。该能带的态密度投影分析显示它主要来源于间隙Fe原子的d电子态。这个杂质能带就像在禁带中打入了一个“钉子”将理论带隙从0.85 eV显著降低至0.40 eV。SbTa反位缺陷同样引入了一个杂质能带但其位置在费米能级上方、导带底附近。这个能带的电子态更加离域不仅来源于占据Ta位的Sb原子也涉及其周围的Fe原子。其效果是将带隙减小至0.80 eV。4.3 对热电性能的潜在影响这种电子结构的改变具有双重效应带隙减小使理论值更接近实验测量值0.53 eV部分解释了理论与实验的差异。引入杂质态在带隙中引入的局域化或半局域化态会成为额外的载流子散射中心可能影响载流子迁移率。同时它们也可能提供额外的载流子输运通道。具体是提升还是降低电导率强烈依赖于费米能级的位置即掺杂水平。在我们的中性缺陷计算中未额外掺杂费米能级钉扎在杂质能带附近这通常会增加电子热导率并可能降低塞贝克系数对热电优值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 局域模如何散射声子这些缺陷引入的局域振动模式是降低热导率的“功臣”其物理机制主要体现在两方面共振散射当传播的声子行波的频率与缺陷局域模的频率接近时会发生强烈的共振散射。局域模就像一个“陷阱”吸收声子的能量极大地缩短了声子的寿命和平均自由程。IFe缺陷引入的局域模正好位于声学支顶部这个频率区域的声子对热导贡献很大因此散射效果尤为显著。破坏周期性增强非谐性缺陷破坏了晶体完美的平移对称性使得原子间的力常数发生局域扰动。这不仅增加了声子-缺陷散射与频率的平方成正比对高频声子散射更强还可能增强三阶非谐力常数从而加剧声子-声子Umklapp散射这是高温下热阻的主要来源。实操心得MTP对声子计算精度的验证。图4(a,b)的对比至关重要。它显示了用MTP有限位移法计算的声子态密度与纯DFTDFPT和有限位移法结果高度一致。这直接证明了我们训练的MTP势函数在描述原子间相互作用特别是二阶力常数上具有DFT级别的精度。只有在完成这一步验证后我们才有足够的底气使用MTP去计算更复杂的含缺陷体系以及用于提取三阶力常数的大规模MD模拟。这是整个工作可信度的基石。6. 晶格热导率的定量计算与实验吻合所有分析的最终落脚点是定量预测缺陷对宏观可观测物理量——晶格热导率κL的影响。我们通过求解玻尔兹曼输运方程BTE实现了这一目标。6.1 计算方法与参数设置我们使用KALDO软件包在弛豫时间近似RTA下求解BTE。所需的输入是体系的二阶谐波和三阶非谐力常数。力常数获取对于完美和含缺陷体系我们使用基于训练好的MTP进行长时间纳秒级的NPT系综分子动力学模拟。从平衡后的轨迹中采用相关函数法提取出所有原子对之间的二阶和三阶力常数。为了对比完美TaFeSb的三阶力常数也使用了更精确但昂贵的DFTDFPT方法通过D3Q代码计算。计算细节采用足够密的q点网格如24x24x24对布里渊区进行采样以确保声子谱和群速度的收敛。计算温度范围从100K到800K覆盖了热电应用的主要温区。6.2 计算结果分析与讨论图5展示了我们的核心计算结果可以清晰地看到几个关键现象完美晶体理论值高估无论是基于DFTDFPT最精确还是MTPMD计算完美TaFeSb的理论κL曲线在室温附近都比实验值高出约2-4倍。这直接印证了引言中提到的理论与实验的巨大差距。MTPMD的结果比DFTDFPT还要高约20%说明我们的MLIP在非谐相互作用上可能存在轻微的高估但这不影响定性趋势。缺陷显著降低热导率引入3.7%浓度的IFe或SbTa缺陷后计算出的κL曲线大幅下降向实验数据靠拢。其中IFe缺陷的降低效果比SbTa更明显。这与声子谱的分析一致IFe引入的局域模位于声学支顶部对主要热载流子低频声学声子的散射更强。浓度效应我们计算了更高缺陷浓度12.5%的情况。结果显示κL进一步降低甚至略低于部分实验点。结合缺陷形成能均大于1 eV的事实可以推断实际样品中的缺陷浓度可能接近我们计算的较低值~3.7%并且可能是多种缺陷共存的混合状态。与实验的吻合在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) 训练并验证该材料的MLIP4) 用MLIP-MD研究缺陷效应。关键在于步骤2)的采样要充分且MLIP的验证如声子、弹性常数必须严格。7. 方法论的普适性、局限性与未来展望我们的研究不仅揭示了TaFeSb中特定缺陷的物理效应更验证了“DFTMLIP”这一计算范式在研究复杂材料缺陷问题上的强大能力。7.1 “DFTMLIP”范式的优势与挑战优势精度与效率的平衡在保持接近DFT精度的前提下将声子及热导率计算的能力从数十原子提升到数百甚至上千原子体系使得研究真实缺陷浓度、缺陷团簇乃至晶界成为可能。可扩展性一旦MLIP训练完成它可以被用于各种基于MD的模拟如扩散系数计算、力学性能测试、高温相变研究等性价比极高。高通量筛选潜力结合主动学习算法可以自动化地生成训练数据并迭代改进势函数为大规模筛选有潜力的热电材料缺陷构型提供了工具。挑战与注意事项训练数据的质量与数量MLIP的精度完全依赖于训练数据。对于含有强电子关联性或复杂化学键的材料生成具有代表性的训练数据集本身可能是一项挑战。需要精心设计AIMD的初始条件和采样策略。外推风险MLIP在训练数据所覆盖的原子构型和能量范围内是可靠的但对于远离训练集的极端构型如高缺陷浓度、表面、断裂其预测可能不准确。使用时需明确其适用范围。计算流程的复杂性整个流程涉及DFT计算、AIMD模拟、MLIP训练、力常数提取和BTE求解等多个环节需要熟练使用多种计算软件如Quantum ESPRESSO, LAMMPS, phonopy, KALDO等并对每一步的物理含义和参数设置有深刻理解。7.2 未来研究方向基于本工作有几个方向值得深入探索多缺陷协同效应实际材料中往往多种缺陷共存。可以研究Fe间隙和SbTa反位缺陷之间的相互作用以及它们与掺杂原子的协同效应如何进一步调控声子谱和热导率。有限温度与动力学效应我们的形成能计算是0K下的静态值。在实际生长温度下缺陷的形成熵和振动自由能也会影响其平衡浓度。可以结合MLIP进行高温MD模拟直接观察缺陷的动态形成与湮灭过程。超越弛豫时间近似RTARTA是BTE的一种简化在低温或某些材料中可能不准确。未来可以使用MLIP结合更精确的迭代法求解BTE或直接使用非平衡分子动力学NEMD计算热导率进行交叉验证。面向工程的材料设计将这套计算框架集成到材料设计平台中。针对目标热电性能如高ZT值逆向设计最佳的缺陷类型、浓度和分布为实验合成提供明确的理论指导。回过头看这项工作的最大价值在于它提供了一个强有力的案例通过融合第一性原理的精度与机器学习的效率我们能够以前所未有的清晰度窥见微观缺陷如何左右宏观性能。在热电材料乃至更广阔的功材料领域这种“从缺陷中要性能”的思路辅以先进的计算手段正引领我们从经验摸索走向理性设计。对于计算材料学的研究者而言掌握并发展此类多尺度模拟方法无疑是应对未来复杂材料体系挑战的关键技能。