GPUMD Si/Ge 热接触与固定边界

本例使用鸢算科技集群的 GPUMD 5.9.1 P100 构建（GCC 13.3、CUDA 12.9），用 ASE 3.29、NumPy 和 Matplotlib 建模、提取和画图。每个原生任务只用一张 P100。所用数据是经典势实际计算结果。

下载本篇配套输入、真实结果与脚本，解压进入目录。先查看已有结果：

bash run.sh

analysis/RESULTS.json 为结果汇总，CSV 可继续作图，figures/ 为高清图片。这个步骤不需要 GPU。要重新计算，另解压一份，进入新目录：

sbatch compute.slurm

脚本按小体系静态计算、短热接触的顺序运行，申请一张 P100、一个 CPU 核心。inputs/ 是输入，results/ 是配套参考输出，work/ 是你新运行的结果。不要把新输出追加到原参考目录。完成后自动提取 CSV 和图片。

1. 先区分界面与替位合金#

本例是 Si 与 Ge 两块金刚石结构沿 x 方向相接，使用同一个官方双元素 Tersoff 1989 势。晶格常数统一取 5.5 Å，形成受应变的共格教学模型。它不是实验应变已消除的 Si/Ge 界面，也不是元素随机混合的合金。

13×4×4 常规胞，共 1664 原子，晶胞 71.5×22×22 ų。x 非周期，y/z 周期；x=33 Å 附近为元素分界。按 x 单元分组：最左和最右端共同标作固定组 0，中间 11 个区域为组 1–11。组 1 是热源，组 11 是冷源。

head -n 6 inputs/sige-contact-30/model.xyz
cat inputs/sige-contact-30/CASE.json
cat inputs/sige-contact-30/run.in

关键字段为 pbc="F T T"、group:I:1。组 0 虽然在空间上分为左右两块，仍是一个分组；不能把它当一块连续的温度测量区域。

2. 固定两端，先做预处理与平衡#

potential potential.txt
fix 0
minimize sd 0.0001 3000
velocity 300 seed 160031
time_step 1
ensemble nvt_nhc 300 300 100
fix 0
run 5000

这里用 sd 做固定端的原子预处理。在此版本中，SD 接收固定组参数；不要把它直接替换为不接收固定组参数的 FIRE。fix 0 每一段 run 前都要重新写，防止下一段恢复移动。

velocity 300 seed 160031 初始化可重复的初始速度，本模型未在 XYZ 里预写速度；若 XYZ 已有 vel，这个温度参数不会覆盖原有速度。随后每步 1 fs，平衡 5000 步，即 5 ps。

3. 加冷热浴并输出分组数据#

ensemble heat_nhc 300 100 30 1 11
fix 0
compute 0 10 10 temperature
dump_thermo 100
dump_xyz 1000 snapshot.xyz precision double mass velocity force potential group_labels
run 20000

平均温度参数为 300 K，delta_T=30，所以目标热源 330 K、冷源 270 K，两个浴的温差是 60 K。第二个条件把 delta_T 改为 60，对应 360/240 K。耦合参数 100 指 100 个时间步，本例正好是 100 fs。

compute 0 10 10 temperature 对分组方法 0 每 10 步取样，平均 10 次后输出一行，即每 0.1 ps 一行。生产段共 20 ps，得到 200 行。最后两列是热源、冷源热浴累计能量，前 12 列为组温度，合计 14 列。

4. 提取温度与热浴能量#

head work/sige-contact-30/compute.out
head analysis/sige-contact-30-profile.csv
head analysis/sige-contact-30-history.csv
head analysis/sige-contact-30-energy-budget.csv

Python 读取示例：

import numpy as np
c = np.loadtxt('compute.out')
T_groups = c[:, :12]
Q_hot_thermostat = c[:, -2]
Q_cold_thermostat = c[:, -1]
T_last_10ps = T_groups[100:].mean(axis=0)

组 0 温度为零，分析移动原子的温度时不纳入它。固定位置应按原子编号比较，周期方向用最小镜像；只看坐标从 0 跳到盒长，不能判断固定失败。固定原子的最终输出力可能包含约束反力，主要核对固定坐标及零速度。

5. 看懂能量符号并检查能量收支#

热浴累计能量是“热浴得到的能量”；原子从热源获得能量时，热源列通常向负方向变化。原子获得的净热量是 -(Q_hot+Q_cold)。

此版本 heat_nhc 不刷新全局热力学缓存，因此这个阶段不要用 thermo.out 或 XYZ 头部的旧总势能做随时间的能量判断。包内脚本改读每帧的实际 energy_atom，以及实际 mass/vel 重建动能，再与同一时间的热浴累计能量对应。这些输出是由 dump_xyz ... potential velocity mass 明确请求的。

20 个每 1 ps 的快照中，E_atoms+Q_hot+Q_cold 相对第一帧的最大偏差分别约 5.60、6.02 meV。速度单位是 Å/fs，不能直接与原子质量相乘而忘记动能单位换算；完整转换在 analyze.py 中。

6. 什么时候才能算界面热导#

这 20 ps 主要供学习建模、固定端和数据提取。图中材料内部尚有明显瞬态，不能直接报告收敛界面热导。应先延长平衡与生产，分块比较温度剖面和热流是否稳定，检查热浴大小、尺寸、时间步和势函数。

稳定后分别在界面两侧排除热浴区域拟合温度，再外推到同一界面位置求温差，不能拿冷热浴温差直接充当界面温差。本例 open-x 只有一条传热通道，热流密度不能沿用周期双通道的除 2。横截面积为 22×22=484 Ų，但确定物理模型与稳态以后才有意义。

GPUMDkit NEMD工具可用于已有均匀体系流程；跨界面应分两侧处理温度，不能把整段一条直线的斜率当作本征热导率。本包用 ASE/NumPy/Matplotlib 处理真实分组和单位，没有把不适用的均匀材料拟合套在界面上。

继续看替位合金、层间势双层模型与GPUMD 冷热源。依据：固定组、热输运积分器、compute.out格式。

用ASE检查周期轨迹与MSD#

ASE坐标展开与MSD处理本篇真实20个存储帧，先展开周期方向，再切出1408个移动原子，核对1ps帧间隔与时间单位。短非平衡固体轨迹不报告材料扩散系数。