GPUMD HNEMD驱动力与统计误差

先下载本批输入、参考数据与绘图脚本，解压进入 full-function-wave1。参考数据在 cases、tools，自己的新计算在 work，避免覆盖或向旧输出追加数据。

unzip full-function-wave1-*.zip
cd full-function-wave1
python start.py

1 同时准备种子和驱动力#

菜单选择 2：HNEMD多种子与驱动力。512原子Si在300 K下，分别使用 1e-5、2e-5、5e-5 Å⁻¹ 的x方向驱动力，每个驱动力四个初速度种子。

cat work/gpumd-hnemd-1e-05-s27101/run.in
bash submit.sh

完整输入：

potential potential.txt
velocity 300 seed 27101
time_step 1
ensemble nvt_nhc 300 300 100
run 50000
ensemble nvt_nhc 300 300 100
compute_hnemd 1000 1e-05 0 0
dump_thermo 1000
run 1000000

HNEMD持续输入能量，生产阶段使用Nose–Hoover链恒温器维持300 K。不要直接照搬EMD的NVE生产阶段，也不要用Langevin恒温器替换这里的NHC。参数和适用条件见GPUMD官方HNEMD说明。程序与编译分支同环境指南。

2 从每个种子提取一个估计#

head work/gpumd-hnemd-1e-05-s27101/kappa.out
python analyze.py
python plot.py

自己的同组计算结束后，用 python analyze.py --root work 提取。第1、2列相加得到本例的κxx；每1000步输出一个1 ps区间平均，本例有1000行，即1 ns。先在每个种子内做时间平均，再对四个种子的估计求平均和区间，不能把1000个相关时间点当作1000次独立重复。

import numpy as np
kappa = np.loadtxt("kappa.out")[:, :2].sum(axis=1)
print("单条轨迹均值：", kappa.mean())
running_mean = np.cumsum(kappa) / np.arange(1, len(kappa) + 1)

3 本次1 ns结果#

误差条很宽，最弱驱动力的一条短轨迹甚至给出负的均值。这是短时间信噪比不足的提示，不能把负值截为零，也不能只选数值最顺眼的一条。图中区间使用四个独立初速度的样本标准差，公式同EMD统计教程。

4 怎样安排自己的计算#

保持结构、势、时间步和温度不变，先延长每个种子的生产时间。例如把 run 1000000 改为 run 20000000，在1 fs时间步下得到20 ns。观察前半段、后半段和整段均值，再比较种子之间的离散。20 ns只是可用的下一档长度，是否足够仍由结果判断。

较强驱动力可能提高信噪比，也可能离开线性响应。应比较多个驱动力下的长时间结果和误差，而不是只追求更平滑的曲线。当前1 ns数据不足以确定线性响应区间，也不能与EMD任意积分窗口数值硬作一致性结论。

列定义：kappa.out。

补充：同一组种子的20 ns采样#

四个初速度种子沿用前文，每个场强分别延长至20 ns，另加入1E-4 Å⁻¹场强。延长轨迹不是增加独立种子数。以下误差仍来自四个种子均值的95% t区间，不能单独证明无限尺寸热导率已收敛。

完整run.in、势文件和kappa.out见补充数据包中的 gpumd。将生产段改为 run 20000000，其他设置保持；做初次流程练习时先使用短采样，再决定是否延长。

强场结果的解释： 1E-4 Å⁻¹得到约2450 W·m⁻¹·K⁻¹，显著高于较弱场强下约210–275的结果，已表现出明显场强依赖；不能将该点作为线性响应热导率或与弱场结果合并平均。温度保持300 K并不保证驱动力足够小。初次复现优先用较弱场强和短轨迹检查流程，材料定量结论仍需场强、尺寸和采样收敛。