GPUMD 周期非平衡热流实操

直接非平衡分子动力学（NEMD）在结构中设置热源和冷源，读取温度剖面与热浴交换的能量。本篇先把一个周期算例完整跑通；下一篇继续计算界面温差和热阻。

1 下载、画出结构#

下载输入、16组原生结果和绘图脚本，解压进入目录：

unzip gpumd-SiGe-periodic-NEMD-wave227*.zip -d nemd-example
cd nemd-example
sha256sum -c MANIFEST.sha256
bash view.sh
ls analysis/figures

view.sh 使用 ASE、NumPy 和 Matplotlib 处理已有结果，不提交新计算。图片来自真实轨迹末帧和原生输出。结构、温度分区及两条热流路径要一起看。

体系沿 x 分成 Si、Ge 两半，晶格常数统一为 5.5 Å。这是人为指定晶胞的共格受应变模型，不代表两种材料各自的零压晶格。晶胞三个方向均为周期边界，化学界面有两处：x=Lx/2，以及周期边界 x=0。纯 Si 采用同一几何分区作对照。

坐标沿 x 平移了 a/8，避免把原子平面直接压在界面分界上。各原子带有一个整数 group；热源为第6组，冷源为第19组，编号从0开始。它们分别位于 Lx/4 和 3Lx/4，两条周期热流路径等长。

2 准备新的计算目录#

python3 prepare.py
cd cases/p100-SiGe-Nx16-s13579
cat config.json
cat job.slurm

新目录与 reference-results 分开，避免原生输出追加到旧结果。config.json 中的 nx、system、seed 规定结构与初速度；先按附件原值跑通。提交脚本执行 run_native.py，生成 model.xyz、potential.txt 和完整 run.in，随后启动 GPUMD。

本例采用 GPUMD 5.9.1、GCC 13.3、CUDA 12.9，P100和V100 16G分别使用对应模块；单任务一张卡、一个CPU核。附件的 job.slurm 用于对应GPU核时分区，套餐用户通过 yskit 取得与账户及GPU型号匹配的模板，保留程序运行命令即可。不要同时套用多个 mpirun 进程启动同一个目录。

3 读懂完整输入#

SiGe 使用原版 SiGe_Tersoff_1989.txt，纯 Si 使用其中 Si 参数相同的 Si_Tersoff_1989.txt，不是 NEP 模型。势文件已包含在包内，复现时不要随意换势。

potential potential.txt
time_step 1
minimize fire 0.001 5000
velocity 300 seed 13579
ensemble nvt_bdp 300 300 100
run 50000
ensemble heat_bdp 300 100 30 6 19
run 100000
ensemble heat_bdp 300 100 30 6 19
compute 0 1000 1 temperature potential
dump_thermo 1000
dump_xyz 20000 trajectory.xyz velocity force precision double
run 400000

依次做固定晶胞的原子位置优化、50 ps整体恒温、100 ps热流预热、400 ps正式采样。时间步为1 fs。这里的优化不改变晶胞，不能理解为晶格已松弛到零压。

heat_bdp 300 100 30 6 19 中，300 K是中间温度，30 K是两侧各自偏移，因此目标温度为 330/270 K，总温差60 K；100为以时间步计的耦合参数。局部热浴的实际均温需要从数据读取，不能直接把目标值写成实测值。

compute 0 1000 1 temperature potential 表示用第0种分组、每1000步采样、每1次采样输出。因而每行是同一时刻的瞬时温度、势能与热浴累计能量，每1 ps一行。这里不做输出内时间平均，便于核对同步能量收支。正式阶段共400行，分析舍去前100 ps，平均101–400 ps的300个采样点。

4 提交、检查结果#

sbatch job.slurm
squeue -u "$USER"
tail -n 20 native.out
wc -l compute.out

任务完成后应有 model.xyz、run.in、compute.out、thermo.out、trajectory.xyz 和 final.xyz。正式采样每20 ps保存一帧，共20帧；final.xyz 是最后一帧的原始字节副本。

本例使用的 heat-BDP 执行分支中，thermo.out 的动能、势能以及 XYZ 全局 energy/stress 字段没有随每步刷新。读取当前能量应使用 compute.out；XYZ可用于坐标、速度和力的展示，不要把其中的全局energy字段当作当前训练标签。这一点也说明，只有作业结束还不够，必须看列定义与输出是否实际变化。

5 提取温度和热流#

cd ../..
module load ase/3.29.0
python3 analyze.py --cases cases --output analysis-new
head -n 6 analysis-new/p100-SiGe-Nx16-s13579/profile.csv
head -n 6 analysis-new/p100-SiGe-Nx16-s13579/bath-energy.csv

compute.out 的54列依次为：26个组温度（K）、26个组的总势能（eV）、热源与冷源热浴累计能量（eV）。组势能已经是该组总量，不要再乘原子数；组温度转换为动能时需要乘各组原子数。

热浴能量采用“从体系转移到热浴”的符号：热源向体系供能，斜率为负；冷源带走能量，斜率为正。用101–400 ps线性拟合两者斜率，取其绝对值平均作为总传热功率。

本例横截面积为 A=272.25 Ų。周期结构有两条等价传热路径，单条路径的热流密度为：

Qdot = (abs(hot_slope) + abs(cold_slope)) / 2       # eV/ps
J = Qdot * 1.602176634e-7 / (2*A*1e-20)           # W/m²

斜率绝对值平均中的 /2 和周期两条热流路径中的 /2 是两个不同原因。单路径固定端模型不能照搬第二个因子。

例如P100短盒、种子13579的SiGe模型，热源和冷源斜率分别为 −0.051227和+0.057012 eV/ps，得到单条路径 J=1.59245×10⁹ W/m²；实际热源/冷源均温约 332.680/274.525 K。

6 短采样波动较大时延长计算#

独立新轨迹的冷热浴热流在400 ps内仍可能相差较大。一次V100短采样约为22%；不要改种子去挑好看的数据。保留原目录，另外准备工作副本，将正式阶段延长至120万步（1.2 ns），再比较热浴斜率与能量收支。附件程序支持在 config.json 中加入 production_steps：

cd ..
unzip gpumd-SiGe-periodic-NEMD-wave227*.zip -d nemd-long
cd nemd-long
python3 prepare.py
python3 - <<'PY'
import json
from pathlib import Path
p=Path('cases/v100-SiGe-Nx16-s13579/config.json')
a=json.loads(p.read_text())
a['production_steps']=1200000
p.write_text(json.dumps(a,indent=2)+'\n')
PY

随后进入 cases/v100-SiGe-Nx16-s13579，检查模板并 sbatch job.slurm。请在新的、尚未运行的工作目录使用此设置；不能直接在已有结果目录重跑。前面的50 ps恒温和100 ps热流预热保持不变。最后使用同一分析命令，脚本按实际步数读取数据，仍舍去正式阶段前100 ps。BDP热浴本身还使用按时钟初始化的随机序列，velocity ... seed固定的是初速度，不保证新MD轨迹逐字节相同。参考16组保留原400 ps数据，新增长轨迹属于采样长度检查，不是又一个独立初速度样本。

7 同步核对能量#

python3 - <<'PY'
import numpy as np
from ase.io import read
c='reference-results/p100-SiGe-Nx16-s13579/'
a=read(c+'model.xyz')
n=np.bincount(a.arrays['group'].reshape(-1))
x=np.loadtxt(c+'compute.out')
ke=1.5*8.617343e-5*(x[:,:26]*n).sum(axis=1)
pe=x[:,26:52].sum(axis=1)
e=ke+pe+x[:,52]+x[:,53]
print('max drift (eV/atom):',np.max(abs(e-e[0]))/len(a))
PY

这个参考例的最大漂移约 1.33×10⁻⁵ eV/atom。应比较同一时刻的体系能量与热浴能量，不要把时间平均的体系能量和区间末端热浴累计值直接相加。

温度剖面能看见界面附近的跃变，但把全长温差直接代入傅里叶定律，会混入界面与热浴效应。下一步按SiGe界面温差与热阻分别拟合两侧。

参数语义见GPUMD热输运系综与compute手册。上述短纳米模型用于学习操作，不作为收敛材料热导或界面热阻。