lammps计算热导率

这张图是本项目用 LAMMPS 算某种材料热导率的两种方法对比：左边是 NEMD（非平衡分子动力学）的温度梯度曲线，右边是 Green-Kubo 的 running integral。NEMD 通过给体系两端加温差，测量热流 J = −κ·∂T/∂x；Green-Kubo 则从热流自相关函数积分得到 κ。本案例 GK 积分收敛到约 9.5 W/(m·K)，两种方法互相验证，结果才可信。

这张图在说什么

左图温度从热端 380 K 线性降到冷端 280 K，中间主体区斜率稳定，说明达到了稳态。线性区必须远离热/冷端，否则边界效应会污染 κ。右图 running integral 前期振荡大，10 ps 后逐渐收敛到 9.5 W/(m·K) 平台。如果积分区间太短，κ 会被显著低估或高估。

NEMD vs Green-Kubo 怎么选

| 方法 | 原理 | 优点 | 缺点 | 适用 |

|——|——|——|——|——|

| NEMD | 施加温差，测稳态热流 | 直观、对大乘员快 | 尺寸效应大，需外推 | 大体系、半导体 |

| Green-Kubo | 平衡态热流自相关 | 无边界扰动、尺寸效应小 | 统计噪声大、收敛慢 | 小体系、液体 |

本项目对晶体材料一般两种都跑：NEMD 给工程级初值，GK 做交叉验证，偏差 < 20% 才采信。

LAMMPS 关键设置

| 参数 | 取值/方法 | 说明 |

|——|——|——|

| 势函数 | Tersoff / SW / ReaxFF / MLIP | 晶格动力学/热导对势函数极敏感 |

| 系综 | NVT（两端）+ NVE（主体） | NEMD 常用；GK 用 NVE |

| 时间步长 | 0.5–1 fs | 含氢体系取 0.5 fs |

| 模拟时长 | GK ≥ 1 ns，NEMD ≥ 5 ns | 保证收敛 |

| 长程静电 | pppm | 离子/极性体系必须 |

| 采样间隔 | 5–10 fs | 热流自相关函数需要高频采样 |

NEMD 计算步骤

1. 先 NPT 平衡到目标温度、零压，得到平衡结构。
2. 沿热导方向建长条模型，两端各取 1–2 nm 做热浴。
3. 一端升温、一端降温，维持恒定温差 ΔT。
4. 主体区温度线性分布稳定后，统计热流 J 和温度梯度 ∂T/∂x。
5. κ = J / (A·|∂T/∂x|)，A 为截面面积。

尺寸效应与外推

NEMD 的 κ 随模型长度 L 增加而增大，因为长波声子散射受尺寸限制。本项目通常做 3–5 个不同长度（20、40、80、160 nm）的外推，按 1/κ ∝ 1/L 外推到体相。不做外推的 NEMD 结果可能比实验低 30–50%。

Green-Kubo 收敛判断

GK 公式：κ = (1/3Vk_BT²) ∫₀^∞ ⟨J(0)·J(t)⟩ dt

– 积分上限要覆盖热流自相关函数衰减到 0 的时间；对晶体通常 20–50 ps。

– 需要多个独立初始构象取平均，降低统计噪声。

– 本案例 running integral 在 10 ps 后进入平台，20 ps 基本稳定，取 10–20 ps 区间平均得 κ = 9.5 W/(m·K)。

这套结果怎么用

热导率是散热材料、热电材料、电池热管理的核心参数。低 κ 适合做热电转换（大温差），高 κ 适合做散热/封装。本项目算完会把 κ 与温度关系（300 K、400 K、500 K）一起给出，指导材料选型。

更多分子动力学模拟与热导率计算流程可参考科研学术网 [lammps计算热导率](https://www.keyanxueshu.com/category/md/) 栏目，完整案例见首页 https://www.keyanxueshu.com/。

局限

经典势函数不能描述声子-电子相互作用、量子效应（低温声子）；缺陷、晶界、界面热阻需单独建模；NEMD 尺寸效应必须外推。