GPUMD NEMD谱热导与有效通道

本篇从周期冷热浴出发，演示如何选局部测量区、取得谱热导，再构造有效经典传输通道函数。先读谱热流窗口与单位，这里沿用其角频率换算与原生窗函数约定。

1 下载与查看#

下载谱热流输入、原始结果、完整NEMD原子能量轨迹与分析脚本：

unzip gpumd-Si-SHC-spectral-transmission-wave230*.zip -d shc-example
cd shc-example
sha256sum -c MANIFEST.sha256
bash view.sh
ls analysis/figures

脚本先处理已有数据，再调用GPUMDkit1.5.8绘图；不会提交新模拟。需要GPUMD5.9.1、ASE3.29.0、GPUMDkit1.5.8及对应的NumPy、Matplotlib环境。实际GPU程序使用GCC13.3和CUDA12.9，一任务一张P100或V100 16G、一个CPU核。模板分区须与账户匹配，套餐用户可从 yskit 取得模板。

2 结构与分组要先画清楚#

Si432沿x为6个常规胞，y/z各3个，晶胞32.586×16.293×16.293 Å。沿x均匀分12组，每组36原子，组2为热浴、组8为冷浴。周期边界让热流沿两条路径流动。

第二套分组把原来的空间组4、5、6合为测量组1，其他原子为组0。因此测量组含108原子，长度8.1465 Å，测量体积约2162.585 ų，是全盒体积的四分之一。这里要区分“空间组编号”和“第二套测量组编号”。

3 提交短任务，再做长采样对照#

python3 prepare.py
cd cases/p100-nemd-delta30-s13579
cat config.json
cat job.slurm
sbatch job.slurm

真实200 ps生产输入如下；时间步1 fs，先20 ps NHC，再100 ps冷热浴预热：

potential potential.txt
minimize fire 0.001 5000
time_step 1
velocity 300 seed 13579
ensemble nvt_nhc 300 300 100
run 20000
ensemble heat_bdp 300 100 30 2 8
run 100000
ensemble heat_bdp 300 100 30 2 8
compute 0 1 1000 temperature potential jp
compute_shc 5 1000 0 160 120 group 1 1
dump_thermo 1000
dump_xyz 2000 trajectory.xyz velocity potential precision double
run 200000

compute_shc ... group 1 1 表示第二套分组、组1。另一种子样例使用 group 1 -1，指该套分组除0外的所有组；本结构只有一个非零组，所以两者选择同样的108个原子。更复杂的结构有多个非零组时，两者不能等同。

需要1 ns对照，可提交已有独立目录：

cd ../p100-nemd-delta30-s13579-long1ns
cat config.json
sbatch job.slurm

1 ns案例把生产步数改为1,000,000，仍单独从建模、平衡开始，并非把200 ps文件接着追加。冷热浴还有额外随机噪声，相同初速度种子的两次作业不保证轨迹前缀相同。

4 导出局部谱与能量数据#

wc -l shc.out compute.out
head -3 compute.out
cd ../..
module load ase/3.29.0
python3 analyze.py --cases cases --case p100-nemd-delta30-s13579 --output analysis-new
ls analysis-new/p100-nemd-delta30-s13579

compute.out每1 ps一行，62列包含12组温度、12组势能、36个分量优先排列的组热流、2个热浴累计能量。200 ps与1 ns分别有200和1000行。分析采用实际热浴温差，不能把输入中的30直接当成最终平均温差。

轨迹每2 ps保存一次原子速度和原子势能，分别有100、500帧。当前GPUMD5.9.1的此 heat_bdp 路径不能直接依靠全局能量头判断能量守恒，本例从原子势能加速度动能重建体系能量，再加同一时刻热浴能量；附件保留完整原子数据，方便独立提取。

5 周期双路径与测量体积#

从两个热浴累计能量的斜率估算总功率后，要除以2，得到单条周期路径的功率，再除截面积和实际冷热浴温差。局部SHC已经只选中其中一条路径上的测量区，换算时不要再除一次2。

局部谱热导 G(ν) 单位为W/(m²·K·THz)。使用附件给出的谱热流时，换算系数为 1.602176634e13/(V_measure·ΔT)，其中V用ų、ΔT用K；不能错用全盒体积，也不能把谱热导单位当成谱热导率单位。

随后定义 M_eff(ν)=G(ν)·A_m²/(k_B·10¹²)。这是面积与经典热容加权的有效通道函数，可超过1，有限采样时也可出现负估计；不是0到1的单通道透射概率，也不代表进行了量子动力学计算。

6 为什么要看多个实例和采样时长#

200 ps短任务适合先跑通输入与分析，但局部维里–速度相关的噪声明显大于热浴功率统计。资料包保留8个短NEMD实例，以及Δ=30 K条件下两种初速度、两个平台的4条独立1 ns对照，全部参与比较：

四条1 ns谱积分均值约0.44098 GW/(m²·K)，热浴估计均值约0.43091 GW/(m²·K)，但单条结果仍相差明显。不能仅因均值接近就宣布收敛；也不能把不同GPU上的统计波动解释成硬件改变材料性质。这里是两组初速度在两个平台的独立运行，不是四种独立初速度种子。

本例为同质硅的有限周期盒演示，不是异质界面热导率测量。研究自己的界面时，还需定义合理的两侧温度外推、界面温降和谱测量区，检查长度、温差、势函数与采样收敛。负值和偏差保留在CSV中，先延长采样与增加重复，不要通过裁剪负值制造平滑的“成功曲线”。

最终得到：原生 shc.out、局部带符号谱CSV、温度剖面、能量功数据以及有效通道图。需要实际材料热导率而不做谱分解时，可先用GPUMD专题中的EMD、NEMD或HNEMD完整流程。