GPUMD SiGe界面温差与热阻
对真实SiGe周期轨迹分别拟合界面两侧温度，用两条周期路径的热流求热阻，并比较纯Si对照、盒长和随机种子。
阅读本篇目录
本篇接着周期NEMD实操的16组计算，演示如何选拟合区间、外推界面温差和计算界面热阻。先保留每条轨迹、每个界面的带符号结果，再比较尺寸与种子，不取绝对值强行制造正热阻。
1 重现温度与拟合#
下载SiGe NEMD资料包，解压后执行：
cd nemd-example
bash view.sh
ls analysis/p100-SiGe-Nx16-s13579
每个目录有 profile.csv 和 bath-energy.csv。前者给出分区编号、名义中心、原子数、实际平均温度及镜像分区温度；后者保留每1 ps的热浴累计能量、同步体系能量和总能量变化。所有拟合值位于 analysis/SUMMARY.json。
体系有两个化学界面。热源在Si半区，冷源在Ge半区；热流从热源分别沿正、负x方向传向冷源。图中的一组虚线是中间界面两侧的线性拟合，不能把穿过热浴和界面的一条全长直线当成均匀材料温度梯度。
2 选择两侧拟合区间#
这里固定26个分区。s=x/Lx，第g组的名义中心是 (g+0.5)/26。从101–400 ps取平均，避开热浴和界面紧邻区。
| 分支 | 左侧拟合窗口 | 右侧拟合窗口 | 外推位置 |
|---|---|---|---|
| A：穿过中间界面 | Si：0.29<s<0.46 | Ge：0.54<s<0.71 | x=Lx/2 |
| B：穿过周期边界 | Ge：0.79<u<0.96 | Si：1.04<u<1.21 | x=Lx |
第二分支把 s<0.25 的坐标加1，得到展开坐标u。原子和数据没有移动，只是为了对周期边界两侧作连续的坐标表示。纯Si对照的两侧均为Si。
热浴、强非线性界面层和只有少量原子的分区可能影响拟合。替换自己的结构时，先画图，再按同一规则给各组选择足够的拟合点；不要为得到期望结果逐个调整窗口。本例每侧4个分区，附件保留原窗口。
module load ase/3.29.0
python3 - <<'PY'
import json
s=json.load(open('analysis/SUMMARY.json'))
m=next(x for x in s['metrics'] if x['name']=='p100-SiGe-Nx16-s13579')
for k in ['fit_branchA_left_KperA','fit_branchA_right_KperA',
'interface_A_signed_deltaT_K','interface_B_signed_deltaT_K',
'two_path_flux_Wm2','interface_mean_signed_R_m2K_W']:
print(k,m[k])
PY
拟合数组依次为斜率（K/Å）和截距（K）。B分支表中先列Ge再列Si，而温差仍按“热端Si外推值减冷端Ge外推值”定义；不要仅按横坐标左右顺序机械相减。
3 外推温差、计算热阻#
对每个界面，在同一个界面位置分别代入两侧拟合线：
DeltaT_A = T_Si_fit(Lx/2) - T_Ge_fit(Lx/2)
DeltaT_B = T_Si_fit(Lx) - T_Ge_fit(Lx)
R_A = DeltaT_A / J
R_B = DeltaT_B / J
R_mean = (R_A + R_B) / 2
J 是上篇已经除过两条周期路径的热流密度。此处不要再除2，也不能对热阻再乘横截面积。
以P100短盒、种子13579为例：
| 量 | 实际结果 |
|---|---|
| J | 1.59245×10⁹ W/m² |
| 界面A温差 | 30.9233 K |
| 界面B温差 | 34.5377 K |
| 两界面热阻平均 | 2.05536×10⁻⁸ m²·K/W |
两个界面处在同一条轨迹内，不能把它们当成两次独立重复，也不能把300个时间采样点当成300个独立初速度样本。
4 比较长度、种子与纯Si对照#
下表每行使用两条独立初速度轨迹。先在每条轨迹内平均两个界面，再给出两种子的平均及最小—最大值。单位统一为 10⁻⁹ m²·K/W。
| GPU | Nx | 体系 | 两种子均值 | 两种子范围 |
|---|---|---|---|---|
| P100 | 16 | 纯Si对照 | -0.3007 | -0.3459—-0.2555 |
| P100 | 16 | Si/Ge | 17.1740 | 13.7944—20.5536 |
| P100 | 24 | 纯Si对照 | -0.2111 | -0.2290—-0.1931 |
| P100 | 24 | Si/Ge | 18.0597 | 16.3974—19.7221 |
| V100 16G | 16 | 纯Si对照 | -0.2370 | -0.2883—-0.1858 |
| V100 16G | 16 | Si/Ge | 23.7699 | 20.1971—27.3427 |
| V100 16G | 24 | 纯Si对照 | -0.0128 | -0.0153—-0.0103 |
| V100 16G | 24 | Si/Ge | 14.0071 | 13.5619—14.4523 |
纯Si没有Si/Ge化学界面。对它强行在同样几何位置拟合，也会得到小的非零甚至负的“跃变”，反映有限长度、统计波动和拟合模型的偏差；不要称为负的物理界面热阻。
SiGe例子在这两种长度和短采样下的结果仍有明显种子及平台差异。这张表是在展示真实数据如何比较，不能据此宣布热阻已收敛，也不能从两个长度外推可靠的无限尺寸值。5.5 Å共格应变、有限横截面、固定势函数与330/270 K热浴共同限定了这些结果。
需要材料定量结果时，应增加独立种子、延长稳态采样、扩大横截面并检查热流驱动强度、拟合窗口及有限尺寸；每次只改少数因素，保留原始数据，不用平滑或取绝对值代替统计。
5 处理自己新跑的数据#
module load ase/3.29.0
python3 analyze.py --cases cases --output analysis-new
python3 - <<'PY'
import json
s=json.load(open('analysis-new/SUMMARY.json'))
for m in s['metrics']:
print(m['name'],m['hot_mean_K'],m['cold_mean_K'],
m['two_path_flux_Wm2'],m['interface_mean_signed_R_m2K_W'])
PY
结果目录必须使用新名称，以免把不同计算混在一起。脚本包含固定的本例几何、时间步和分区检查；换成自己的体系，需要同步修改这些定义，而不是删除提示继续计算。
本例的绘图使用ASE读取结构、NumPy线性拟合、Matplotlib出图。GPUMDkit 1.5.8的现有NEMD绘图入口采用另一种分区排列，不能直接套到本例“26组、热源6、冷源19”的双周期路径布局。对于已有工具，先核对分区、边界和列定义，不能只换文件名就套用。
有关温度和热浴能量列，见GPUMD compute输出。完整数据、固定窗口、带符号结果和脚本均保留在下载包内。