晶格热导率计算：从声子BTE到非谐效应的数值实现

晶格热导率计算在热电材料筛选和热管理设计中是不可回避的一环。这个项目曾围绕ZrSi2的热导率预测展开了长达三周的反复计算——初始结果比实验值高出一倍，排查非谐力常数和q点采样的过程揭示了这类计算中暗藏的多重数值陷阱。

晶格热导率的微观表述基于声子玻尔兹曼输运方程（BTE），在弛豫时间近似下：

κ_L = (1/V) Σ_{q,ν} C_{q,ν} v_{q,ν}² τ_{q,ν}

其中C为声子比热，v为群速度，τ为弛豫时间。前三项依赖谐和近似下的声子色散，而τ的核心来源是三声子散射过程，需要三阶非谐力常数（IFC）的支撑。

项目最初低估了三阶力常数计算的复杂度。二阶力常数（谐和）通过有限位移法获取相对直观——位移幅度0.03 Å，对称性约简后约几十个超胞计算即可。但三阶力常数的独立项数目急剧膨胀：ZrSi2原胞含6个原子，三阶力常数经对称性约简后仍有超过4000个独立位移构型需要计算。

项目采用Phonopy-Phono3py工具链完成三阶力常数的提取。超胞尺寸的选择是一个关键决策点——2×2×2超胞（48原子）与3×3×3超胞（162原子）的结果差异显著：

– 2×2×2超胞：κ = 18.3 W/mK（300K）

– 3×3×3超胞：κ = 14.1 W/mK（300K）

– 实验值：约12.5 W/mK

2×2×2超胞的结果偏高30%，根源在于长程声子散射通道未被充分采样。根据Physical Review B相关研究的建议，三阶力常数计算的超胞需满足最短截断半径至少覆盖第三近邻。项目最终采用2×2×2超胞配合8 Å截断半径，而非更大的超胞——在精度提升和计算可行性之间做出了这个取舍。

位移幅度也需要谨慎选取。项目测试了0.01、0.03、0.05 Å三个幅度，发现0.01 Å时力常数噪声偏大（受SCF收敛精度限制），0.05 Å则引入了非简谐效应的高阶项污染。0.03 Å是多数文献推荐的标准值，项目采纳了这个经验。

弛豫时间τ的计算涉及对所有三声子散射过程的求和，这依赖q点网格的密度。项目从11×11×11逐步加密到31×31×31：

– 11×11×11：κ = 19.2 W/mK

– 19×19×19：κ = 15.8 W/mK

– 25×25×25：κ = 14.5 W/mK

– 31×31×31：κ = 14.2 W/mK

差距不会说谎——稀疏q点严重高估热导率，因为它遗漏了大量小q值散射通道。项目最终采用31×31×31网格，此时结果收敛至1%以内。

值得一提的是，声子寿命的计算中存在一个容易被忽略的细节：等权重近似与简并处理的差异。某些高对称方向上存在声子简并，如果不做特殊处理，散射相空间会被人为放大。项目在处理Γ点附近的声学支简并时遇到了数值发散，最终通过引入小宽化参数（σ = 0.1 THz）将发散抑制。

天然材料中同位素质量涨落引入额外散射。项目采用Tamura模型计算同位素散射率，对ZrSi2的影响约为5-8%——天然丰度下Zr有五种稳定同位素，质量涨落参数g = 1.7×10⁻⁴。根据Nature报道的热导率测量研究，同位素纯化可使某些半导体的热导率提升数倍，但在金属间化合物中影响相对温和。

加入同位素散射后，计算值降至13.2 W/mK，与实验值12.5 W/mK的偏差缩小至5.6%。

回过头审视，这套计算框架的可靠性高度依赖两个假设：玻恩-奥本海默近似和三声子散射主导。对于高温下四声子散射不可忽略的体系（如金刚石、MgO），仅考虑三声子过程会高估热导率。项目在计算ZrSi2时温度限于800 K以下，恰好规避了这一限制。

另一个值得警醒的局限在于电子热导的剥离。ZrSi2是金属间化合物，实测热导率包含电子和晶格两部分。项目通过Wiedemann-Franz定律从总热导率中扣除电子贡献（Lorenz数L = 2.44×10⁻⁸ WΩ/K²），剩余部分才与纯晶格热导率计算对比。如果忽略这一步，理论与实验的”偏差”将毫无物理意义。

项目最终建立的流程是：谐和力常数收敛验证 → 三阶力常数超胞测试 → q点密度扫描 → 同位素修正 → 电子贡献剥离。每一步的数值决策都有据可查，结果才经得起同行审视。晶格热导率计算的价值不在于得到一个数字，而在于理解这个数字背后每一步近似引入的不确定度。