GPUMD 谱热流与频率窗口实操
用Si216运行compute_shc，从原生相关函数得到频率谱，核对角频率、窗函数和GPUMDkit的带符号数据导出。
阅读本篇目录
谱热流分析用来观察不同频率对热输运的贡献。本篇先用Si216的HNEMD演示 compute_shc，重点说明“相关时间”“角频率”“频谱网格”和输出单位；非均匀冷热浴的操作放在下一篇。
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 提交一个HNEMD谱热流任务#
python3 prepare.py
cd cases/p100-hnemd-x-Nc1000-F5e-05-s13579
cat job.slurm
sbatch job.slurm
squeue -u "$USER"
Si216立方盒长16.293 Å，Tersoff-mini势，300 K、1 fs时间步，20 ps平衡后生产200 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 nvt_nhc 300 300 100
compute_hnemd 1000 5e-05 0 0
compute_shc 5 1000 0 160 120
compute 0 1 1000 jp
dump_thermo 1000
dump_xyz 20000 trajectory.xyz velocity potential precision double
run 200000
compute_hnemd产生驱动，compute_shc测量维里–速度相关函数和谱热流，不会再施加第二个驱动力。参数方向0/1/2对应x/y/z；附件有y、z对照及两档驱动力、两组初速度。
3 分清三个时间与频率参数#
| 参数 | 本例 | 如何理解 |
|---|---|---|
| MD时间步 | 1 fs | 积分步长 |
| SHC采样间隔 | 5步=5 fs | 相邻相关采样间隔 |
Nc |
1000 | 正负相关延迟为−4.995至4.995 ps |
num_omega |
160 | 输出160个频率点 |
max_omega |
120 | 角频率上限，换成普通频率为19.0986 THz |
| 生产轨迹 | 200 ps | 用于统计相关函数的总生产时长 |
原生频率网格是0.75、1.50、…、120 rad/ps。普通频率 ν=ω/(2π)，因此普通频率间隔约0.119366 THz。手册将角频率也用THz标示，画通常的频率横轴时仍要除以2π。SHC参数与分组
4 读取原生输出#
wc -l shc.out kappa.out compute.out
head -5 shc.out
cd ../..
module load ase/3.29.0
python3 analyze.py --cases cases --case p100-hnemd-x-Nc1000-F5e-05-s13579 --output analysis-new
ls analysis-new/p100-hnemd-x-Nc1000-F5e-05-s13579
shc.out前1999个数据行是相关函数，后160行是角频率谱，每行三列。原生相关函数已经乘过Hann窗，二次加窗会改变谱。附件用原生相关行独立余弦变换，并保留文本有效位数带来的差异。
用 Nc=200 对照时，正负延迟约−0.995至0.995 ps。长短窗口都输出160点，但独立分辨率主要由相关窗口决定，约分别对应0.1和0.5 THz的尺度，另有窗函数展宽；增加网格点数不会凭空产生更高分辨率。
5 从谱热流换算谱热导率#
附件CSV提供普通频率、谱热流与谱热导率，后者单位W/(m·K·THz)。零频率点由相关函数求和推导，明确标在结果说明中；原生160点没有输出零点。积分时同时记录频率上限和相关窗，不能把有限带宽积分强行归一到总热导率。
| 本次200 ps实例 | 数值/W/(m·K) |
|---|---|
| 独立总热流的HNEMD平均 | 113.3734 |
| SHC零延迟对应总量 | 107.6342 |
| 0至19.0986 THz有限带宽积分 | 107.5714 |
平均温度约299.537 K。几种估计器在短轨迹上存在差别；它们用于理解处理流程，不代表已经获得硅的收敛材料参数。
6 用GPUMDkit查看，并保留原始负值#
在参考结果对应目录执行：
module load gpumdkit/1.5.8
cd reference-results/p100-hnemd-x-Nc1000-F5e-05-s13579
gpumdkit.sh -plt hnemd 1 20 --save --save-data
ls hnemd.png data_shc.npz data_shc.txt
本例1 fs时间步、1000步HNEMD输出，匹配工具默认1 ps横轴。需要定量积分时，注意GPUMDkit1.5.8的默认 k_g_wt 总谱会把负值替换为0.001。原始有符号分量 k_g_wi、k_g_wo仍可使用，应将这两个分量相加，另存为带符号总谱。
资料包保留原工具导出，不覆盖它们，另写 *-toolkit-signed-spectrum.csv。在附件根目录查看：
python3 read_toolkit_data.py
head -5 analysis/p100-hnemd-x-Nc1000-F5e-05-s13579-toolkit-signed-spectrum.csv
这一步应在 view.sh 已生成的 analysis 上执行；自定义目录用 --analysis 目录名。只提取数值数组，不需要反序列化NPZ中的对象型Results。工具TXT中Kt的单位标注省略了Å，原生相关函数此处实际为Å·eV/ps。不要把用于显示的负值处理带入定量积分。