LAMMPS计算热导率：两种方法的实操经验与对比

热导率是材料热管理性能的核心参数，结合自己做过的多个体系，聊聊LAMMPS计算热导率的两种主流方法——平衡MD与非平衡MD的实操经验。

【1. 为什么用MD算热导率】

热导率(κ)描述了材料传递热量的能力，在芯片散热、热界面材料、航空航天热防护等场景里是核心设计参数。实验测量热导率成本高、条件严苛，MD模拟可以给出原子尺度的传热机制，而且计算成本比第一性原理低几个数量级——一个几百原子的体系，用EMD算热导率在我的工作站上跑几天就能拿到结果。

我在做氮化硼(BN)纳米片热导率的时候，第一次用MD算出来发现：单层BN的热导率高达600 W/m·K，比很多金属还高——这个结果让我非常意外。后来查了文献，发现这其实和石墨烯(~3000 W/m·K)的高热导率是同一个物理来源：强共价键网络和低声子散射率。MD把这个现象复现出来了，让我对这个方法多了一层信任。

【2. 平衡分子动力学(EMD)：Green-Kubo方法】

Green-Kubo方法的思路是：根据涨落-耗散定理，热流自相关函数的时间积分等于热导率。公式是 κ = (1/(k_B T² V)) ∫₀^∞ ⟨J(t)·J(0)⟩ dt，其中J是热流矢量。这个方法在均匀体系里效果最好，不需要构造温度梯度，统计上更稳定。

“`

# LAMMPS Green-Kubo EMD 示例

fix J all ave/correlate $v &

type map command file J.dat &

auto $1

# 计算热流自相关函数并积分

# 输出热导率

variable k equal trap(f_J)/${T_const}/${T_target}/${V})

“`

我一般在做完能量最小化和NPT平衡后，跑50-100 ps的NVE系综采集热流数据。热流自相关函数的收敛性是关键——如果曲线没有衰减到零就停止积分，得到的热导率会严重高估。

【3. 非平衡分子动力学(NEMD)：我的标准做法】

NEMD的思路更直观：构造一个温度梯度(一边加热、一边冷却)，测量稳态热流，用傅里叶定律 κ = q / ∇T 求出热导率。这个方法更接近实验的真实测量方式，也更容易理解。

“`

# LAMMPS NEMD 示例

# 设置热源和热汇

fix 1 hot heat 1 10.0 ${T_hot}

fix 2 cold heat 1 10.0 ${T_cold}

# 每N步计算一次热流

compute myKE all ke/atom

compute myPE all pe/atom

compute flux all heat/flux myKE myPE

# 热导率输出

fix evflux all ave/time 10 100 1000 c_flux[1] c_flux[2] c_flux[3] file flux.out

“`

我的NEMD标准设置：热源和热汇各10层原子，温差控制在10-20 K以内(太大会有非线性效应)，中间插5-10层的”缓冲层”来平滑温度分布。这些参数不是固定的，需要根据体系做调整——比如二维材料和三维体系的热阻结构完全不同。

【4. 尺寸效应：有限尺寸修正】

MD模拟的体系尺寸(尤其是模拟盒子在热流方向的长度)对热导率结果有显著影响。NEMD的热导率随盒子长度增加而增大，最终趋于一个极限值;EMD也有类似问题，但尺寸效应相对小一些。

我的做法是：用至少三个不同长度(50 Å、100 Å、200 Å)的盒子分别算，然后做 extrapolation 到无穷大。这个修正在低热导率材料(聚合物、液体)里尤其重要——有限尺寸效应可能造成50%以上的偏差。

【5. 常见问题与误差来源】

势函数的选择是热导率计算里最大的不确定性来源。同一体系用不同势函数，热导率可能差出2-3倍。我在重要项目里会对比至少两种势函数的结果，差异太大的话需要进一步分析原因——是势函数本身对声子传输的描述有问题，还是参数化时的条件与当前体系不匹配?

另一个常见问题是NEMD的温度梯度不线性——这通常说明热源和热汇的能量注入/抽取速率太高，导致局部非平衡效应。降低fix heat的频率或者增加缓冲层可以缓解。