GPUMD GKMA模态热流实操
先生成本征矢，再提交Si216的GKMA，导出模态与总热流互相关积分并对照原生HAC结果。
阅读本篇目录
本篇演示GKMA如何把总热流分解到振动模态，再计算各模态组对热导率的带符号贡献。先按本征矢准备完成匹配结构；这里使用Si216、300 K、1 fs时间步和50 ps生产轨迹。
1 下载与查看#
下载Si模态分析输入、真实结果和绘图脚本，解压后先查看已有结果：
unzip gpumd-Si216-GKMA-HNEMA-wave229*.zip -d modal-example
cd modal-example
sha256sum -c MANIFEST.sha256
bash view.sh
ls analysis/figures
view.sh只处理附件中的参考结果，不提交新计算。新计算使用普通账号，一任务一张GPU、一个CPU核。本例环境为GPUMD5.9.1、GCC13.3、CUDA12.9；绘图使用ASE3.29.0、NumPy和Matplotlib。附件提供P100与V100 16G模板，套餐用户可用 yskit 取得匹配账户的分区模板，保留其中的程序启动命令。
2 按顺序提交两个任务#
在附件根目录执行，第一次先生成案例目录：
python3 prepare.py
cd cases/basis-p100-d0.005
basis_job=$(sbatch --parsable job.slurm)
cd ../p100-gkma-hac-control
sbatch --dependency=afterok:$basis_job job.slurm
squeue -u "$USER"
如果已执行过 prepare.py，不要重复创建案例目录。afterok让GKMA在本征矢成功生成后启动。模板会复制对应本征矢，再生成下面的实际 run.in：
potential potential.txt
time_step 1
velocity 300 seed 13579
ensemble nvt_nhc 300 300 100
run 20000
ensemble nve
compute_gkma 10 1 648 bin_size 24
compute 0 10 1 jp
compute_hac 10 500 1
dump_thermo 1000
dump_xyz 5000 trajectory.xyz velocity force potential precision double
run 50000
先用NHC平衡20 ps，再用NVE采样50 ps。compute_gkma 10 1 648 bin_size 24 表示每10步采样，包含第1到648个模态，每24个模态为一组，共27组。compute 0 10 1 jp 独立输出总热流，采样时刻与GKMA对齐；本例额外启用 compute_hac 10 500 1，方便在同一条轨迹核对积分。
每个run只能启用一个模态分析动作，GKMA和HNEMA应放在不同任务目录。GKMA参数
3 查看输出，再提取与画图#
wc -l heatmode.out compute.out hac.out
head -3 heatmode.out
head -3 compute.out
cd ../..
module load ase/3.29.0
python3 analyze.py --cases cases --case p100-gkma-hac-control --output analysis-new
ls analysis-new/p100-gkma-hac-control
ls analysis-new/figures
实际生产轨迹有5000个采样时刻。heatmode.out 每个时刻输出27行，每行五个分量，依次为x-in、x-out、y-in、y-out、z；先按“时刻×模态组×分量”还原，再对模态组求和。不能把所有行直接按时间平均。
4 积分时必须保留交叉项#
每个模态组的相关函数应为 C_b(τ)=<J_b(t) J_total(t+τ)>。将每组仅与自身做自相关再相加，会遗漏组间交叉项。附件脚本用FFT计算组与总热流的互相关，并以直接点积检查部分延迟。
本例采样间隔10 fs；500个相关点对应0至4.99 ps。4.99 ps是相关延迟，不是50 ps的MD生产时长。积分采用300 K分母以直接对照原生HAC；NVE实际平均温度约298.3865 K，另列在结果中。
| 实算检查量 | P100结果 |
|---|---|
| 模态求和与独立JP最大差 | 1.20×10⁻⁶，GPUMD原生热流单位 |
| NVE最大能量偏移 | 7.34×10⁻⁶ eV/原子 |
| 模态相关求和与原生HAC最大差 | 3.36×10⁻⁹，原生相关函数单位 |
| 模态积分求和与原生HAC最大差 | 5.49×10⁻⁷ W/(m·K) |
V100 16G同类对照的积分最大差为3.15×10⁻⁷ W/(m·K)。这说明两条数据处理路径吻合，不能据此断言50 ps足以获得材料收敛热导率。
5 读取各组贡献#
find analysis-new/p100-gkma-hac-control -maxdepth 1 -type f
cat analysis-new/SUMMARY.json | less
附件CSV保留频率中心、每组模态数，以及不同积分延迟处的带符号贡献。例如第1组24个模态、平均频率约3.0141 THz，在1、3、4.99 ps的贡献分别约0.08873、0.26273、0.32104 W/(m·K)。部分其他组为负，代表有限采样下的互相关贡献；不要截断为零后再求和。
更换模态分组可使用 f_bin_size 1.0 按1 THz分箱。先核对模态求和，再比较频率组贡献；扩大结构、延长生产、增加初速度重复后，才能讨论材料结论。需要均匀驱动的模态方法时，继续HNEMA操作。