GPUMDkit NEMD温度与热流绘图
用实际冷热端Si计算检查分组、能量交换和有限长度热输运。
阅读本篇目录
下载输入、实际结果与绘图脚本。以下命令在解压后的目录执行。
module load gpumdkit/1.5.8
bash run.sh
cat LAST-RUN.txt
run.sh在新的work目录中复跑工具，不会重新启动DFT或分子动力学。下面也列出每一步的单独命令。GPUMDkit为1.5.8、Python3.11；GPUMD/NEP为5.9.1，实际GPU计算使用P100单卡。
1. 模型分组要与工具一致#
本例1152原子Si，16×3×3常规胞。x为开边界，两端各一个常规胞属于固定组0；中间由左到右为组1–14，热端组1、冷端组14，y和z周期。使用Tersoff Si_Fan_2019，种子81、82，各平衡20 ps、生产100 ps。
本版nemd绘图入口按“热端组1、冷端为最后一组”的布局读取数据。不能把周期盒内组0加热、组6冷却的旧布局直接套入；周期盒还存在两条传热路径，热流分配也不一样。更改绘图输入不能替代重新建立物理模型。
potential potential.txt
velocity 300 seed 81
time_step 1
ensemble nvt_nhc 300 300 100
fix 0
run 20000
ensemble heat_lan 300 100 30 1 14
fix 0
compute 0 10 100 temperature
dump_thermo 1000
dump_xyz 10000 trajectory.xyz
run 100000
fix 0在每段重新声明。compute 0 10 100 temperature每1000步输出一次，即每1 ps一行；每行15组温度，再加冷热端交换能量两列，总计17列。
2. 先看原始数据，再调用绘图#
cd reference/md/nemd-81
head -n 2 compute.out
tail gpumd.log
gpumdkit.sh -plt nemd Auto 1 20 --save --save-data
ls nemd.png data_nemd.txt data_nemd.npz
cd ../../..
Auto根据模型分组确定冷热端中心距离，本例7.0603 nm；第二参数1表示使用三维实体截面积，无真空厚度修正。第三参数20仅用于存在SHC文件时的频率截断，本例没有shc.out，因此不会产生频谱热流曲线。
3. 怎么读工具结果#
| 种子 | 后半段冷热端温差 / K | 原生端点法κ / W·m⁻¹·K⁻¹ | 后半段能量线性拟合法κ / W·m⁻¹·K⁻¹ |
|---|---|---|---|
| 81 | 47.639 | 9.441 | 9.896 |
| 82 | 51.036 | 8.649 | 8.851 |
| 工具先计算热导G，再乘冷热端距离得到这里的有效κ。冷热端存在温跃，小模型中间温度并非理想直线，因此不能将它叫作已收敛的体相热导率。 |
本版原生端点法用第51至100行的能量差，却以50 ps作除数，实际两端时刻相差49 ps；它与直接线性拟合还存在统计取样差别。上表分别保留原生数值和按真实时间线性拟合的结果，不混成一组数据。下面可直接计算后半段冷热端功率：
import numpy as np
x = np.loadtxt('compute.out')
t = np.arange(1, len(x)+1) # 每行1 ps，本例第一行是1 ps
i = len(t)//2
qin = -np.polyfit(t[i:], x[i:, -2], 1)[0]
qout = np.polyfit(t[i:], x[i:, -1], 1)[0]
print('hot/cold, eV/ps:', qin, qout)
print('mean heat power:', (qin+qout)/2)
种子81和82冷热端功率差分别约2.64%和1.50%。只有一次运行时原生输出的±0来自重复样本数为1，不意味着误差为零。换成自己的体系，应增加独立重复、延长稳态统计、比较长度，并检查固定层、温跃和截面积定义。