1. 热导率分子动力学计算工程导论

热导率描述的是材料在温度梯度下传导热量的能力。宏观上，它由傅里叶定律定义：热流密度 J 与温度梯度 ∇T 成正比，比例系数即热导率 κ，写作 J = −κ∇T。负号表示热流方向与温度升高方向相反，即热量自发从高温区流向低温区。这一约定在后续所有公式、图示与数据分析中保持一致：温度梯度指向升温方向，热流矢量指向降温方向。热导率的国际单位是 W/(m·K)，它不是一个纯数，而是同时包含能量输运能力与几何尺度的量纲。

从微观角度看，热量在固体中并不由某种“热流体”携带，而是通过原子振动和电子运动传递。对绝缘体和半导体，晶格振动（声子）是主要载体；对金属，自由电子贡献显著。分子动力学（MD）处理的是原子核的运动，因此它天然适合描述声子主导的晶格热导，也能通过经验势或机器学习势部分反映电子—声子耦合的间接影响，但不直接求解电子输运。理解这一点，是判断“什么时候该用 MD、什么时候该用玻尔兹曼输运方程或第一性原理”的前提。

MD 计算热导率的核心思想，是把宏观的傅里叶定律“翻译”成原子尺度的可测量量。常用方法有两类：平衡态分子动力学（EMD）和非平衡态分子动力学（NEMD）。EMD 基于 Green–Kubo 公式，通过热流自相关函数的积分得到热导率：κ = (1/(3Vk_BT²)) ∫⟨J(0)·J(t)⟩dt。NEMD 则人为施加温度梯度或能量交换，直接测量稳态热流与温度梯度之比。两者各有适用边界：EMD 在收敛良好时统计效率高，NEMD 更直观但受有限尺寸效应影响明显。

工程实践中，MD 计算热导率并不是“跑一个脚本就出数”。它是一条从物理建模到数值验证的完整链路。以典型的硅晶体为例：先建立含数千到数万原子的立方超胞，采用 Stillinger–Weber 或 Tersoff 势描述 Si–Si 相互作用；在 NPT 系综下弛豫到目标温度（如 300 K）和零压，得到平衡晶格常数；随后切换到 NVE 系综，用 Velocity-Verlet 积分器以 0.5–1 fs 步长推进。若用 EMD，需要先运行足够长的平衡段（例如 200 ps）消除初始关联，再采集热流自相关函数，相关时间通常取 5–20 ps，并对多个独立初态做系综平均。若用 NEMD，则把超胞沿 x 方向分成若干段，两端分别设置热源和冷源，中间为测量区，运行 1–2 ns 直到温度剖面线性化，再由线性拟合斜率与热流计算 κ。

这里有几个容易踩坑的参数关系值得提前建立直觉。第一，尺寸效应：NEMD 中声子平均自由程可能达到数百纳米，而 MD 超胞通常只有几纳米到几十纳米，导致 κ 被系统性低估。常用做法是计算多个长度 L 下的 κ(L)，再用 1/κ(L) = 1/κ∞ · (1 + λ/L) 外推到无限长。第二，时间步长与截断半径：步长过大会导致能量漂移，过小则浪费算力；截断半径一般取 2.5–3.5 倍势函数特征长度，并配合尾部修正。第三，温度与量子修正：MD 是经典统计，低温下声子热容被高估，常需引入量子修正因子将经典结果映射到真实温度。第四，势函数保真度：同一材料不同势给出的 κ 可能相差数倍，必须与实验或第一性原理数据交叉验证。