说明：本文华算科技主要介绍晶格热导率的物理定义、二阶与高阶力常数如何进入声子玻尔兹曼输运方程、RTA 和迭代 BTE 与 Green–Kubo 的差别，以及数值收敛和二维归一化的核查方法。
晶格热导率在计算里是什么量？
晶体处在微小温度梯度中时，声子携带的热流密度满足傅里叶定律。沿笛卡尔方向写成 Jα = −ΣβκLαβ∇βT，其中 κLαβ 是二阶张量。立方晶体的三个对角分量相同，低对称晶体则可能沿不同晶轴给出不同数值。常见单位 W·m−1·K−1 采用三维体积归一化。
线性化声子玻尔兹曼输运方程把宏观张量写成模态求和。令 λ = (q, s) 表示波矢与声子支，原胞体积为 V0，均匀 q 点数为 N，则可写为
κLαβ = 1/(NV0) Σλ CλvλαFλβ。
Cλ 是模态热容，vλα = ∂ωλ/∂qα 是群速度，Fλβ 是碰撞方程求出的广义平均自由位移。弛豫时间近似下有 Fλβ ≈ τλvλβ，熟悉的 Cv2τ 形式便从这里出现。
图1左侧给出一条可追溯的计算链：平衡结构经过原子位移与自洽力计算后，得到二、三、四阶原子间力常数；这些量生成声子频率、散射率与碰撞算符，再由 PBTE 给出 κL。图中还加入自洽声子重整化和非对角热流，说明 常规“谐波声子加三声子散射”只是多种计算近似中的一种，强非谐晶体可能需要更高阶描述。
模态热容由玻色分布和 ℏωλ/kBT 决定，群速度来自色散斜率，Fλ 则由散射算符控制。三者必须在同一结构、温度和 q 点采样下进入求和。电子热导率 κe 不包含在 κL 内；金属和高载流子浓度半导体与实验总热导率比较时，还要处理电子贡献以及可能存在的双极扩散项。
三维体材料的实验值可用于核验整条计算链。图2把 30 种晶体在 300 K 的计算值与实验值放在对数坐标中，蓝色圆点对应显式求解非谐 PBTE 的 AAPL 结果，分布整体靠近等值线；只采用准谐近似的快速估算散点更大。这里检验的是跨材料数量级与排序能否复现，单个样品的晶界、同位素组成和缺陷浓度仍需单独对齐。
原子位移和DFT力怎样生成κL？
平衡结构附近的势能面可按位移展开：E = E0 + 1/2 ΣΦ(2)uu + 1/6 ΣΦ(3)uuu + 1/24 ΣΦ(4)uuuu + ···。DFT 对位移超胞给出的力满足 F = −∂E/∂u，拟合这些力就能恢复展开系数。Φ(2) 决定动力学矩阵及声子本征值，结构质量、电子自洽精度与超胞尺寸会直接传入频率和群速度。
Φ(3) 生成三声子相互作用矩阵元。允许的吸收与发射过程同时满足能量守恒和晶体动量守恒：ωλ ± ωλ′ = ωλ″，q ± q′ = q″ + G。矩阵元、可用相空间和占据数共同给出散射率。τλ 不是单独输入的经验常数，在三声子 PBTE 中它由三阶力常数与守恒条件计算得到。
有限位移通过多组受扰结构的 DFT 力拟合 Φ(n)，密度泛函微扰理论则在响应方程中求谐波量；两种路线获得的物理对象相同，数据组织方式各异。对低对称或大原胞材料，独立参数会迅速增加，稀疏回归可用较少位移结构筛选有效项。图3中，硅的三阶截断半径改变拟合误差和非零参数数目，训练结构从 2 个增加时，300 K 热导率出现明显波动后才落入参考区间。
四阶力常数有两类用途：一类进入四声子散射率，增加可耗散的过程；另一类通过自洽声子理论修正有限温度频率和本征矢。两种效应对 κL 的作用方向可以相反。频率硬化可能抬高群速度或压缩三声子相空间，四声子散射则增加线宽。只在谐波频率上补一个四声子寿命会漏掉频率重整化，高温软模晶体尤其明显。
图4统计了 562 种材料在多种计算近似下的比值。自洽声子修正可使部分材料的热导率上升，也有材料下降；加入四声子散射后，比值整体小于 1。BP、SiC 等高热导材料的高阶修正较温和，若干卤化物的幅度很大。是否加入 Φ(4) 应由温度、软模和三声子结果的偏差迹象决定，不能把同一近似等级无条件套到所有晶体。
RTA/迭代BTE/Green–Kubo算的是？
线性化 PBTE 可写成散射矩阵作用于未知偏离分布的方程。RTA 只保留对角衰减项，每个模态以自身 τλ 回到平衡；迭代 BTE 保留模态间耦合，反复求解 Fλ。Normal 过程守恒总晶体动量，它会重排非平衡分布，却不像 Umklapp 过程那样直接耗散热流。Normal 过程占比高时，对角近似与迭代解的差距往往扩大。
两种 PBTE 解法共用声子频率与非谐力常数，差别来自碰撞算符的求解。若散射矩阵近似对角，RTA 已能接近迭代值；若出现声子流体行为或强集体漂移，应采用迭代解。RTA 和迭代 BTE 的标签必须与结果一同报告，因为同一套 Φ(3) 仍可产生不同的 κL。
晶粒散射常以附加速率写进模态寿命，例如 τ−1B,λ ∼ |vλ|/L。图5a、b显示硅的温度曲线和累积热导率随特征尺度 L 改变；图5c则把四方与单斜 ZrO2 的 AAPL 结果和 Green–Kubo 结果放在一起。本征三声子 PBTE 与有限晶粒样品对应不同散射集合，和实验对照时应明确样品尺度。
平衡分子动力学直接计算体系总热流随时间的涨落。若 J 表示体系总能流，Green–Kubo 关系写成 καβ = 1/(VkBT2) ∫0∞⟨Jα(t)Jβ(0)⟩dt；若代码输出的是已除以体积的热流密度，前因子相应变为 V/(kBT2)。热流变量的归一化约定会改变公式前因子，单位检查可发现这一类错误。
Green–Kubo 结果还受原子间势、体系尺寸、轨迹长度、独立样本数和相关积分窗口控制。经典分子动力学采用经典模态占据，在低于德拜温度时与量子热容有差别。准谐 Green–Kubo（QHGK）把热流相关函数写到振动模态基底中，可处理晶体与无序固体。图6对非晶硅给出 QHGK 和平衡分子动力学的温度曲线，两套结果在中温区相符并带有分子动力学统计误差。
晶体中声子准粒子清楚、线宽远小于频率间隔时，PBTE 的模态图像较直接；强无序、谱线重叠或接近 Ioffe–Regel 区域时，应改用非对角热流和 QHGK/Wigner 描述。选择算法前要先检查频率间隔与线宽的相对大小，软件正常结束只表明数值流程已经运行。
怎样确认晶格热导率已经收敛？
收敛检查应沿数据传递顺序展开。结构弛豫残余力、平面波截断和电子 k 网格先决定 DFT 力的噪声；二阶超胞不足会扭曲声子频率，三阶超胞与相互作用截断不足则会漏掉散射矩阵元。频率看似平滑并不能证明 Φ(3) 已收敛，因为谐波谱主要检验 Φ(2)。
进入输运求和后，还要分别增加 q 网格密度，调整能量守恒的展宽或四面体积分方案，并检查迭代残差。图7给出石墨烯实例：三阶截断从 5.94 Å 增至 6.35 Å 后，纳入 6.19 Å 处的较强相互作用；较短截断下 κL 随 q 网格持续上升，较长截断才形成稳定平台。过短三阶截断会让更密 q 网格延续错误趋势，这种曲线缺少收敛依据。
同位素散射、晶粒散射、点缺陷散射和电子—声子散射要按研究对象选择。理想单晶计算若加入任意晶粒尺度，所得数值已含外禀散射；实验样品若含纳米晶粒，只用本征三声子结果通常会偏高。温度、晶相、同位素组成与样品尺度应和 κL 成对记录，否则数值差异无法归因。
极性晶体还需检查 Born 有效电荷与介电张量产生的非解析修正，尤其当低频光学支参与热输运时。软模结构在谐波近似下出现虚频，基于该谱的常规 PBTE 求和没有稳定的准粒子基础；若目标晶相只在有限温度稳定，应采用相应温度下的声子重整化。虚频来源要在数值误差与有限温度稳相之间分开核查，两种情况对应不同处理。
二维材料的超胞体积含真空层，程序输出的三维化 κL 会随晶胞高度 Lz 成反比变化。把它乘以 Lz 得到面热导 Gs = κLLz，可消除真空厚度的任意选择；若再用物理厚度 h 换回 W·m−1·K−1，应同时写明 h 的定义。二维结果至少报告晶胞高度或选定厚度，文献间才有可复算的换算关系。
可交付的计算结果应包含 κLxx、κLyy、κLzz 或明确的方向平均值，并附温度、PBTE 解法、散射项、力常数阶数、超胞和 q 网格收敛记录。若与实验比较，再补充样品尺度及电子热导率的扣除方式。这样得到的数值对应清楚的哈密顿量、散射集合与归一化约定。
