3. 热导率计算的理论框架

分子动力学（MD）计算热导率，本质上是在“平衡态涨落”与“非平衡态响应”两条路线之间做选择。理解这两条路线的理论根基，才能判断在给定势函数、体系规模和计算预算下，应该用哪种方法、设置哪些参数、如何判断结果是否可信。本节把 Green-Kubo 线性响应理论、Einstein 关系、非平衡分子动力学（NEMD）以及直接方法放在同一框架下梳理，并给出适用场景的对比。

一、线性响应与 Green-Kubo 公式

Green-Kubo（GK）方法来自线性响应理论。其核心思想是：体系在平衡态下自发存在微观热流涨落，而热导率正是热流自相关函数的时间积分。对三维各向同性体系，GK 公式为

κ = (1 / (3 V k_B T²)) ∫₀^∞ ⟨ J(0) · J(t) ⟩ dt

其中 V 是体积，k_B 是玻尔兹曼常数，T 是温度，J 是微观热流矢量，⟨·⟩ 表示平衡系综平均。J 的表达式依赖势函数形式：对两体势，J 包含粒子动能输运项与粒子间相互作用势能输运项；对多体势（如 Stillinger-Weber、Tersoff、EAM），必须使用对应势函数导出的热流表达式，否则会系统性低估或高估 κ。这是 GK 实现中最容易出错的地方。

实际操作中，GK 计算流程为：先在 NVT 或 NPT 系综下把体系平衡到目标温度，再切换到 NVE 系综采样；每隔若干步记录一次瞬时热流 J(t)，采样步长应小于热流自相关时间；然后计算热流自相关函数（HCACF）并做时间积分。判断收敛的关键是 HCACF 是否衰减到零附近并保持噪声水平，以及 κ(t) 的积分曲线是否出现平台。若积分曲线持续上升或剧烈振荡，说明采样不足或热流表达式有误。

二、Einstein 关系：从自相关到均方位移

Einstein 关系是 GK 公式的等价形式，它把热导率与“热流积分”的均方位移联系起来：

κ = (1 / (3 V k_B T²)) lim_{t→∞} (1 / (2t)) ⟨ |ΔR(t)|² ⟩
ΔR(t) = ∫₀^t J(τ) dτ

与 GK 相比，Einstein 关系在长时极限下线性增长更稳定，适合处理 HCACF 长尾衰减慢的体系，例如低热导率或强声子散射体系。但它的收敛同样依赖足够长的采样时间，且对热流积分漂移敏感。实践中常把 GK 与 Einstein 关系作为交叉验证：若两者给出的 κ 在误差范围内一致，结果可信度较高。

三、非平衡分子动力学（NEMD）

NEMD 不走平衡涨落路线，而是主动在体系中建立温度梯度或热流，再根据傅里叶定律 κ = -J / ∇T 反推热导率。常见实现有两种：

• Müller-Plathe 交换法：把模拟盒沿 z 方向分成若干层，定期将冷层中动能最大的粒子与热层中动能最小的粒子交换速度，从而在两端形成温差。测量各层温度后拟合温度梯度，再结合交换能量与时间得到热流。
• 热浴法：在体系两端分别施加 Nosé-Hoover 热浴或 Langevin 热浴，直接维持高温端与低温端，中间区域自然形成温度梯度。

NEMD 的优点是信号强、计算量相对小，适合大体系或需要快速估算的场景。但必须注意有限尺寸效应：当体系长度 L 与声子平均自由程可比时，κ 会随 L 变化。通常需要计算多个 L 下的 κ，再用 1/κ(L) = 1/κ∞ + a/L 外推到无限长极限。此外，热浴强度、交换频率、分层厚度都会影响结果，需要做参数收敛测试。

四、直接方法与非平衡态方法的适用场景对比

“直接方法”通常指直接对平衡态轨迹做 GK 或 Einstein 分析，不额外施加扰动；而“非平衡态方法”指 NEMD 这类主动建立梯度的方案。两者并非互斥，而是互补。选择时主要看体系尺寸、势函数复杂度、目标精度和计算资源。

一个实用的判断原则是：如果体系较小、势函数复杂、且希望获得与实验或第一性原理可比的高精度结果，优先用 GK 或 Einstein 关系，并做多独立轨迹平均；如果体系较大、需要快速扫描参数或研究界面热阻，NEMD 更合适，但必须完成有限尺寸外推和热浴参数收敛。对于各向异性材料，还要分别计算不同方向的热流分量，不能简单取平均。

最后需要强调：无论选择哪条路线，热流表达式、系综选择、采样长度和尺寸外推都是决定结果可靠性的关键。GK 方法中，若 HCACF 在 1 ps 内就衰减到零，说明采样步长可能过大；若 20 ps 后仍有明显振荡，则需延长采样并检查势函数截断。NEMD 中，若温度分布在靠近热浴区域出现非线性弯曲，应剔除该区域后再拟合梯度。把这些判断标准落实到每一次计算中，才能让热导率结果具备可复现性和物理意义。