机器学习势函数结合DFT：揭示缺陷如何降低半赫斯勒化合物晶格热导率

1. 项目概述与核心问题在热电材料的研究领域半赫斯勒化合物一直是个“明星选手”它们拥有不错的电学性能但一个长期困扰研究者的难题是理论计算出的晶格热导率总是比实验测量值高出一大截。这可不是个小问题晶格热导率是决定热电转换效率的关键参数之一理论预测不准后续的材料设计和性能优化就无从谈起。以TaFeSb为例基于完美晶体模型的第一性原理计算给出的室温晶格热导率预测值在18到30 Wm⁻¹K⁻¹之间而实验测得的数值却只有9 Wm⁻¹K⁻¹左右。这个巨大的鸿沟让很多基于理想模型的理论筛选和预测工作变得意义有限。过去大家把目光更多地放在了缺陷对电子结构的影响上因为计算相对容易。但对于缺陷如何影响声子输运、进而改变晶格热导率却研究甚少。原因很简单用传统的第一性原理方法比如密度泛函微扰理论DFPT来计算包含缺陷的大超胞的声子谱计算量是天文数字普通课题组根本玩不转。这篇工作正是要啃下这块硬骨头。我们不再回避缺陷而是主动拥抱它利用机器学习原子间势这把“利器”结合密度泛函理论系统探究了本征点缺陷如空位、反位、间隙原子在半赫斯勒化合物TaFeSb中究竟扮演了什么样的角色。我们的目标很明确第一找出最可能自然形成的缺陷类型第二量化这些缺陷对电子能带结构的影响第三也是最重要的揭示它们如何扰动晶格振动、散射声子最终导致晶格热导率的大幅降低。这不仅是为了解释一个具体的实验差异更是为了提供一套可推广的计算框架让缺陷工程从“经验摸索”走向“理性设计”。2. 方法论详解当DFT遇见机器学习势函数要研究缺陷首先得把模型建对。传统的第一性原理方法精度高但算不动大体系经典分子动力学速度快但势函数精度不够难以描述半赫斯勒化合物中复杂的键合作用。我们的策略是取二者之长用DFT保证电子结构计算的精度用机器学习原子间势来高效处理包含缺陷的大体系晶格动力学计算。2.1 计算框架与软件栈整个研究的技术栈可以概括为“DFT打底MLIP攻坚专用工具分析”。密度泛函理论计算全部在Quantum Espresso软件包中完成。我们采用PBE泛函处理交换关联能使用了投影缀加波方法。为了平衡精度和效率平面波截断能设置为100 Ry波函数和1000 Ry电荷密度。布里渊区采样使用了6×6×6的Monkhorst-Pack k点网格。这些参数对于半赫斯勒化合物这类金属间化合物来说是经过验证的可靠选择能够在可接受的计算成本下获得收敛良好的总能和电子结构。缺陷建模是关键一步。我们构建了3×3×3的超胞来模拟缺陷这大约对应3.7%的缺陷浓度。有同行可能会问这个浓度是不是太高了实际上前期研究表明对于点缺陷形成能的计算这个尺寸的超胞已经能够较好地逼近稀薄极限继续增大超胞对形成能结果的影响很小但计算量会呈立方级增长。在超胞中我们系统考察了四种缺陷空位、反位缺陷、间隙原子和原子交换。机器学习原子间势的核心我们选择了矩张量势。MTP是一种基于局部原子环境描述符的势函数它将体系的总能表达为每个原子在其截断半径内邻居环境贡献的求和。其强大之处在于通过足够的DFT数据训练后它能以接近DFT的精度预测能量、力和应力但计算成本却低了好几个数量级。我们通过LAMMPS进行从头算分子动力学模拟生成了约3000个构型作为训练集其中既包含完美的TaFeSb也包含含有缺陷的构型如TaFe₁.₁₂₅Sb。声子与热导率计算的流程是用训练好的MTP通过PHONOPY软件的有限位移法计算包含缺陷体系的声子谱。获得力常数后再利用KALDO软件包求解玻尔兹曼输运方程在弛豫时间近似下最终得到晶格热导率。作为基准对比完美TaFeSb的声子谱和热导率我们也用DFTDFPT和DFT有限位移法算了一遍以确保MTP的可靠性。注意机器学习势的训练是成败的关键。训练集必须足够“丰富”要涵盖你关心的所有原子构型和温度范围。我们的策略是运行从300K到1100K的NPT系综AIMD模拟以采样不同的原子位置和晶格振动模式。训练时能量、力和应力的权重设置1 0.1 0.1也需要根据体系特性微调过高的力权重可能导致势函数对局部结构过于敏感而忽略长程相互作用。2.2 关键参数与收敛性测试在实际操作中有几个参数需要格外小心MTP截断半径我们设置为5 Å。这个值需要大于所研究相互作用的主要范围。对于TaFeSb5 Å足以包含次近邻原子这对于准确描述声子谱是必要的。截断半径太小会丢失重要相互作用太大会增加计算量并可能引入噪声。超胞尺寸与k点网格对于缺陷计算超胞要足够大以避免缺陷之间的周期性镜像相互作用。我们的测试表明3×3×3超胞对于点缺陷形成能已基本收敛。同时超胞变大会导致k点可采样数减少因此需要测试k点网格的收敛性。对于2×2×2的训练用超胞我们使用了2×2×2的k点网格这在对训练集进行大规模AIMD采样时是效率与精度的折中。声子计算网格用有限位移法计算力常数时需要构建足够大的超胞我们用了4×4×4的超胞来准确描述长波声子。同时后续用BTE计算热导率时需要对倒易空间进行非常密集的q点采样例如30×30×30以准确积分声子态密度和弛豫时间。踩坑心得初期我们曾尝试用更小的训练集或更短的AIMD模拟时间来生成训练数据结果得到的MTP在预测某些高频光学声子支时出现明显偏差。后来发现必须让AIMD模拟充分遍历构型空间特别是在高温下如1100K以采样到更大的原子位移这对于训练势函数准确描述非谐相互作用至关重要。这部分的计算量投入是省不得的。3. 缺陷稳定性与电子结构影响材料中的缺陷不是凭空出现的哪种缺陷更容易形成由它的形成能决定。形成能越低缺陷在热力学上越稳定在实验制备过程中出现的概率就越高。3.1 缺陷形成能排序与物理图像我们计算了TaFeSb中多种本征缺陷的形成能结果一目了然Fe间隙原子的形成能最低约为1.28 eV。紧随其后的是Sb占据Ta位的反位缺陷形成能约为1.50 eV。其他如Ta空位、Fe反位等缺陷的形成能则要高得多普遍在2 eV以上。这个排序背后有清晰的晶体化学图像。半赫斯勒化合物的晶体结构空间群F-43m本身就在(0.75, 0.75, 0.75)位置存在一个本征的空位。这个空位就像一个“预留座位”尺寸和化学环境使得Fe原子坐进去的能量代价相对较小。而SbTa反位缺陷的形成可能与Ta和Sb原子在某些条件下的互扩散倾向有关。这些低形成能的缺陷在材料生长过程中很容易自发形成因此在实际的TaFeSb样品中它们的存在是大概率事件而不是特例。提示计算形成能时化学势的取值非常关键。我们采用了元素固态单质的基态能量作为化学势参考Ta和Fe用bcc结构Sb用三结构。这种取值方式假设材料与元素单质处于平衡是处理这类三元化合物缺陷的常用方法。但如果考虑材料处于非平衡生长条件如某种元素过饱和化学势的取值范围会变化从而影响形成能的绝对值但通常不会改变不同缺陷之间的相对稳定性排序。3.2 缺陷诱导的电子态与能带调控缺陷不仅影响结构更会深刻改变电子结构。我们对完美晶体以及含有3.7%浓度的Fe间隙和SbTa反位缺陷的体系进行了电子态密度和能带计算。Fe间隙缺陷它在费米能级下方、价带顶之上引入了一个窄的杂质能带。这个杂质能带主要来自间隙Fe原子的d电子态局域性非常强。它的出现相当于在原本的带隙中“塞”进了一些电子态导致理论计算的带隙从完美晶体的0.85 eV显著减小到0.40 eV。SbTa反位缺陷同样引入了一个杂质能带但这个能带位于费米能级上方导带底之下。通过投影态密度分析发现这个杂质态的电子云更加离域不仅来自占据Ta位的Sb原子也涉及其周围Fe原子的贡献。它使带隙减小到0.80 eV。这两个结果具有重大意义实验上通过红外光谱测得的TaFeSb带隙约为0.53 eV。无论是完美晶体的0.85 eV还是使用更高级泛函如mBJ计算出的~0.9 eV都与实验值相差甚远。而引入缺陷后理论值特别是Fe间隙缺陷的0.40 eV与实验值的差距大大缩小。这强有力地表明实验中观测到的“更窄的带隙”很可能正是这些本征缺陷存在的直接证据。这些缺陷引入的局域或半局域态会成为载流子的散射中心或额外的输运通道必然会对电导率、塞贝克系数等热电输运性质产生复杂影响。4. 缺陷对声子谱与晶格动力学的扰动晶格热导率的高低本质上取决于声子晶格振动的量子如何被散射。缺陷作为晶格中的“异物”是散射声子的高手。但要定量理解其机制必须从声子谱的微观变化入手。4.1 声子态密度与局域模分析我们首先对比了用不同方法计算完美TaFeSb声子态密度的结果DFPT、有限位移法DFT-FD、以及基于MTP的有限位移法MTP-FD。三者在中低频区域吻合得非常好只在高于5 THz的部分有些微差异。这验证了我们训练的MTP在描述晶格动力学方面具有接近DFT的精度为后续研究缺陷体系打下了信任基础。分析含有缺陷的体系的声子态密度发现了非常有趣的现象声子带隙变窄完美TaFeSb的声子谱在约7 THz处存在一个明显的带隙宽度约0.95 THz。引入Fe间隙缺陷后这个带隙宽度缩小到0.68 THz引入SbTa反位缺陷后带隙缩小到0.74 THz。出现局域模在带隙缩小的区域出现了新的声子模式。对于Fe间隙缺陷在~6.5 THz处出现了一组尖锐的峰。通过投影态密度分析发现这些模式几乎完全局域在间隙Fe原子上周围其他原子的贡献可以忽略不计。对于SbTa反位缺陷则在高于7.25 THz的光学支区域产生了一些新的模式这些模式主要局域在缺陷周围最近邻的Fe原子上被替换的Sb原子本身贡献很小。这里的物理图像很直观间隙Fe原子质量较轻且所处的局域环境与晶格格点原子不同其振动频率会落在原本晶体声子模式的“空白区”带隙内形成局域振动模。而SbTa反位缺陷由于Sb和Ta原子质量、键强差异扰动了周围Fe原子的受力情况从而在其振动频率附近诱发出新的局域化模式。4.2 局域模如何散射声子这些缺陷诱导的局域声子模式是降低晶格热导率的“元凶”之一主要通过两种机制起作用共振散射当传播的声子尤其是能量与局域模接近的声子经过缺陷附近时会与这些局域模发生强烈的相互作用共振导致声子被强烈散射寿命急剧缩短。这类似于光学中的共振吸收。改变声子群速度与相空间局域模的出现可能“压扁”或改变原有声子带的色散关系降低声子的群速度即热传播速度。同时它们为三声子散射过程特别是倒逆过程提供了更多的通道和相空间进一步增强了声子-声子散射。我们的计算清晰地显示Fe间隙缺陷引入的、位于声学支顶端的局域平带对声子的散射作用尤为显著。因为中低频的声学支声子通常是热输运的主要贡献者这些声子被强烈散射直接导致热导率大幅下降。5. 晶格热导率计算与实验差异弥合理论计算的最终检验标准是能否解释或预测实验。我们分别计算了完美TaFeSb、含有3.7%和12.5%浓度缺陷的TaFeSb的晶格热导率随温度的变化并将其与实验数据对比。5.1 完美晶体的理论高估作为基准用最精确的DFTDFPT方法计算完美TaFeSb的晶格热导率在室温附近的理论值大约是实验值的4倍。而使用MTPMD方法虽然效率高但精度略逊计算的结果甚至更高比DFT结果还高出约20%。这个对比首先确认了两件事第一完美晶体模型确实严重高估热导率第二我们的MTP方法在趋势上是可靠的虽然绝对值有偏差但用于研究缺陷引起的相对变化是足够精准的。5.2 缺陷引入带来的显著降低引入缺陷后晶格热导率出现了戏剧性的下降对于3.7%的Fe间隙缺陷室温下热导率理论值下降至接近实验值的两倍左右。对于更高浓度12.5%的Fe间隙缺陷理论曲线在400K至800K的温度区间内与实验数据惊人地接近。SbTa反位缺陷也能降低热导率但效果弱于Fe间隙缺陷这与它引入的局域模能量较高、对主要热输运声子散射较弱有关。温度区间的解读为什么在400-800K区间符合得最好在更低温度下如300K声子输运中“正常散射”过程占主导而弛豫时间近似在处理这类三声子散射过程时本身就有局限。在更高温度下多晶样品中的晶界散射、甚至可能存在的极化子效应等机制开始凸显而这些在我们的单晶模型中没有考虑。因此中温区间最能反映体材料本征的点缺陷散射机制。5.3 对实际样品缺陷浓度的启示计算表明需要相对较高的缺陷浓度如12.5%才能使理论热导率曲线与实验数据在较宽温区重合。然而我们的形成能计算显示即使是最稳定的Fe间隙缺陷其形成能也在1.28 eV以上这意味着在热平衡条件下如此高的缺陷浓度是很难达到的。这看似矛盾实则指向一个更现实的图景实际样品中很可能同时存在多种类型的缺陷Fe间隙、SbTa反位可能还有其他类型。虽然每种单一缺陷的浓度不高但它们的散射效应会叠加。此外样品中可能还存在位错、晶界等扩展缺陷以及制备过程中引入的非平衡缺陷。所有这些散射机制共同作用才能将热导率从理想晶体的高值拉低到实验观测的低值。我们的计算聚焦于最可能存在的两种点缺陷已经能解释大部分差异这证明了点缺陷散射是导致理论实验不符的心物理机制。6. 方法优势、局限与实操建议6.1 机器学习势在此类研究中的优势这项研究成功的关键在于采用了“DFT MLIP”的混合策略。其优势是压倒性的突破计算尺度限制直接用DFTDFPT计算包含缺陷的3×3×3超胞的声子谱和二阶、三阶力常数计算量极其庞大。而用训练好的MTP进行分子动力学模拟来提取力常数计算成本降低了数个数量级使得研究成为可能。保持高精度MTP通过拟合DFT数据能够复现复杂的多体相互作用其精度远高于传统的经验势如Lennard-Jones, MEAM等特别适用于键合性质复杂的半赫斯勒化合物。路径可推广这套“DFT生成训练数据 - 训练MLIP - MLIP驱动大尺度MD/晶格动力学计算”的流程可以无缝应用到其他含有缺陷、界面或非晶态的材料体系中是计算材料学中一个非常强大的范式。6.2 当前方法的局限性与改进方向当然没有方法是万能的我们的研究也存在一些局限也是未来可以改进的方向缺陷浓度与分布我们只考虑了低浓度下孤立的点缺陷。实际材料中缺陷可能存在簇聚、或沿特定方向有序排列这些都会改变散射机制。用更大的超胞和更长的MD模拟来研究缺陷-缺陷相互作用是下一步的工作。散射机制建模我们使用了玻尔兹曼输运方程下的弛豫时间近似。RTA在处理强非谐性或某些正常散射过程占主导的情况时不够精确。未来可以考虑使用迭代法求解BTE或直接采用非平衡分子动力学模拟来计算热导率进行交叉验证。电子-声子耦合本研究聚焦于晶格热导率。但缺陷同时影响电子和声子电子-声子耦合散射也可能改变热输运。要全面评估热电性能需要将电导率、塞贝克系数与热导率放在同一框架下计算考虑其相互关联。机器学习势的泛化能力MTP在训练数据覆盖的构型空间内表现良好但对于训练集未充分采样的极端畸变或全新的原子排列其外推能力存疑。在研究中需要密切关注MD模拟中是否出现了能量或力异常大的构型必要时需将其加入训练集进行主动学习迭代优化。6.3 给同行研究者的实操建议如果你打算在自己的体系中应用类似方法以下几点经验可能对你有帮助训练集构建是重中之重不要只对完美晶体进行AIMD采样。务必包含你计划研究的缺陷构型并在不同的温度下运行模拟以确保训练集能覆盖缺陷可能引起的局部晶格畸变和振动模式。构型数量不是唯一标准质量多样性比数量更重要。先做小规模测试在投入大量资源训练全尺寸势函数之前先用一个小体系如2×2×2超胞和较短的时间步长测试不同的MTP复杂度等级和训练参数快速评估方法的可行性。始终设置DFT基准对于完美晶体或某个关键性质如弹性常数、某个声子频率一定要有纯DFT的计算结果作为基准用于量化MLIP的误差。我们的研究显示MTP对热导率的绝对预测可能有~20%的偏差但对缺陷引起的相对变化捕捉得很准。理解物理而不仅仅是拟合数据当MLIP给出一个出乎意料的结果时比如某个缺陷导致热导率异常升高要回头去分析声子谱、态密度、甚至动画观察振动模式。机器学习是工具背后的物理机制需要研究者自己去解读和洞察。7. 结论与展望缺陷工程的新视角这项研究通过结合密度泛函理论和机器学习原子间势清晰地揭示了本征点缺陷在半赫斯勒化合物TaFeSb的热输运性质中扮演的关键角色。我们不仅从热力学上确认了Fe间隙和SbTa反位是最可能存在的缺陷更重要的是从微观机理上阐明了它们如何通过引入局域的电子态和声子模来显著降低材料的晶格热导率从而弥合了理论与实验之间长期存在的巨大差距。这项工作最大的启示在于它为我们提供了一套强大的计算工具和清晰的物理图像使得“缺陷工程”不再是一个黑箱。过去我们通过实验掺杂、合金化来引入缺陷以降低热导率多少有些“试错”的成分。现在我们可以先通过计算预测哪种缺陷最有效形成能低、散射强甚至可以在设计新材料时就有意引入特定的缺陷结构来优化热电性能。例如是否有可能通过控制生长条件来选择性促进Fe间隙缺陷的形成从而最大化其对热导率的降低效果从更广阔的视角看DFTMLIP的方法组合为研究复杂材料体系中的各类缺陷点缺陷、位错、晶界、非谐效应、以及有限温度下的性质打开了新的大门。它使得在保持第一性原理精度的前提下模拟包含数千个原子、在皮秒至纳秒时间尺度演化的体系成为可能。这对于理解热电材料、超导材料、电池电极材料等众多功能材料在实际工作条件下的微观过程具有不可估量的价值。最后我想分享一点个人在调试计算中的体会机器学习势的训练看似自动化实则非常依赖研究者的物理直觉和经验。最初我们只用完美晶体的数据训练MTP然后直接去算缺陷体系的声子结果高频部分完全失真。后来才意识到缺陷周围的局部力场与完美晶体截然不同训练集必须“喂给”势函数足够多这类畸变环境的例子。这就像教一个AI识别动物不能只给它看站着的猫还得有趴着的、跑着的、甚至炸毛的。对于材料模拟而言构型空间的充分采样是任何机器学习方法取得成功的基础。