2. 分子动力学基础与热导率计算关联

热导率不是分子动力学（MD）模拟直接输出的物理量，而是从原子轨迹中“反推”出来的输运系数。理解这一点的关键在于：MD 每一步只做一件事——在给定势函数和边界条件下，用牛顿运动方程推进原子位置与速度。热导率计算的全部参数选择，本质上都是在回答同一个问题：这条轨迹能否真实反映材料内部的热输运过程？本节围绕牛顿运动方程、势函数、系综、边界条件、时间步长与截断半径六个要素，说明它们如何共同决定热导率模拟的可信度。

牛顿运动方程：轨迹的生成引擎

MD 的核心是求解每个原子的牛顿第二定律：

m_i * d²r_i/dt² = F_i = -∂U/∂r_i

其中 m_i 为原子质量，r_i 为位置，U 为体系总势能，F_i 为原子所受合力。数值实现上普遍采用 Velocity-Verlet 算法，它同时更新位置与速度，时间可逆且能量漂移小。热导率计算依赖速度自相关函数（Green-Kubo 方法）或稳态温度梯度（非平衡方法），二者都要求轨迹在长时间内保持能量守恒和动量守恒。若积分器步长过大，高频声子模式被截断，热流自相关函数会失真，热导率往往被系统性低估。

势函数：决定热输运的物理内核

势函数描述原子间相互作用，直接决定声子色散、群速度和寿命，因此是热导率模拟中影响最大的输入。常见类型包括：

• Lennard-Jones（LJ）势：适合惰性气体与简单液体，参数少、计算快，但无法描述共价键的方向性。
• EAM/MEAM 势：适合金属，能较好再现弹性常数与声子谱。
• Tersoff/Stillinger-Weber 势：适合硅、碳等半导体，可描述键角与配位。
• 机器学习势：精度接近第一性原理，但训练集必须覆盖目标温度与结构。

选择势函数时，应优先核对其声子谱、热膨胀系数与实验值是否吻合。同一材料用不同势函数算出的热导率可能相差数倍，这不是数值误差，而是物理模型差异。

系综与边界条件：控制变量与消除尺寸效应

系综决定哪些宏观量被固定。热导率模拟常用：

• NVE：用于 Green-Kubo 方法的生产段，保证能量守恒，热流自相关函数才有明确物理意义。
• NVT：用 Nosé-Hoover 或 Langevin 热浴控温，用于平衡段。
• NPT：用于获得目标温度下的平衡晶格常数，避免因体积不准引入应力。

边界条件方面，体材料通常采用周期性边界条件（PBC）消除表面效应；一维纳米结构（如纳米线、纳米管）可在输运方向用 PBC，横向留真空。非平衡方法中，常沿输运方向设置热源与冷源，形成稳态温度梯度，再按傅里叶定律 J = -κ ∇T 提取热导率。此时必须确认温度分布在统计上是线性的，否则说明热源功率或尺寸不合适。

时间步长与截断半径：数值稳定性的两条红线

时间步长 Δt 必须小于最快振动周期的约 1/10。对含氢体系，O-H 或 C-H 伸缩振动周期约 10 fs，Δt 应取 0.5 fs 甚至更小；对硅、铜等重元素体系，1～2 fs 通常可行。判断标准是：NVE 段总能量漂移应小于 10⁻⁴～10⁻⁵ eV/atom/ps。

截断半径 r_c 决定势能计算范围。对 LJ 势，r_c 一般取 2.5σ～3.5σ；对 EAM 势，需覆盖电子密度衰减到可忽略的距离。截断过小会人为削弱长程相互作用，导致声子软化、热导率偏低；截断过大则计算量急剧上升。工程上常配合长程校正或切换函数，使势能在 r_c 处平滑归零。

参数如何共同影响热导率结果

六个要素并非独立。势函数决定声子寿命，时间步长决定高频模式是否被正确积分，截断半径决定长程作用是否完整，系综决定热流涨落是否被正确采样，边界条件决定尺寸效应是否被消除，牛顿方程积分器则保证轨迹本身可信。任何一项设置不当，都会在热导率上表现为系统性偏差，而非随机噪声。

一个可复现的参数检查清单

下面给出硅体系 Green-Kubo 热导率模拟的示例参数组合，用于说明各参数如何协同。以下数值为教学示例，实际使用前应做收敛性测试。

最后需要强调，热导率模拟的收敛性必须逐项验证：先固定其他参数，单独改变 Δt、r_c、体系尺寸和关联时间，观察 κ 是否进入平台区。只有当所有参数都收敛，得到的 κ 才具有物理可比性。这也是 MD 热导率计算区别于简单“跑一个模拟”的核心工程要求。