LAMMPS 硅直接NEMD：热浴、热流与温度剖面

下载完整输入、原生温度剖面、热浴功与绘图脚本，保存为 example.zip：

python3 -m zipfile -e example.zip .
cd lammps-direct-thermal-examples
bash results.sh
bash plot.sh

环境为LAMMPS 30Sep2026、Kokkos 5.2.1双精度、C++20 / CUDA 12.8构建，MPI与模块说明见版本环境。本批纯Si在P100、Si/Ge在V100 16G实际运行，每例一张GPU、一个MPI进程。下面的下载脚本可按自己的资源选择一种GPU入口。

reference-results是实际结果，inputs是可继续使用的输入，reference-data是提取后的CSV。查看与作图不启动新的MD。

1. 建什么体系#

Si采用金刚石结构、5.43 Å晶格常数，横截面4×4个常规胞。三种盒长为8、12、16个常规胞，分别1024、1536、2048原子。采用LAMMPS随附 SiCGe.tersoff 中的 Si(D) 参数，原文件保留Tersoff 1989来源说明。

全周期边界，热区位于z/L=0.20—0.30，冷区0.70—0.80，目标温度330、270 K。两区间存在两个周期传热方向。各长度使用两个独立初速度种子；相应目录名和参数在PLAN.json。

2. 先平衡，再小跑，再追加采样#

每例流程为：固定晶胞最小化，300 K NVT平衡10 ps，NVE配合冷热Langevin热浴采样50 ps；随后从该末态的DATA继续200 ps。时间步1 fs。

追加阶段保留位置与速度，重新按空间区域建立热浴，并重置热浴随机种子。这是继续采样的分段操作，不保证完整积分器与随机数状态的逐步重现。需要严格状态续算时，应另行使用与你方法相匹配的restart设置。

以下是Si最小例子的关键输入；完整输入见 inputs/Si-L8-s74011/short.in：

units metal
atom_style atomic
boundary p p p
lattice diamond 5.43
region box block 0 4 0 4 0 8
create_box 2 box
create_atoms 1 box
mass 1 28.085
mass 2 72.63
pair_style tersoff
pair_coeff * * SiCGe.tersoff Si(D) Ge
neighbor 0.5 bin
min_style cg
minimize 1e-8 0.02 1000 10000
velocity all create 300 74011 mom yes rot no dist gaussian
timestep .001
fix eq all nvt temp 300 300 .1
run 10000
unfix eq
reset_timestep 0
region hotzone block INF INF INF INF 8.688 13.032 units box
region coldzone block INF INF INF INF 30.408 34.752 units box
group hot region hotzone
group cold region coldzone
fix integrate all nve
fix hotbath hot langevin 330 330 .1 74111 tally yes zero yes
fix coldbath cold langevin 270 270 .1 74211 tally yes zero yes

纯Si例中没有类型2原子；类型映射为同一份Si/Ge工作流保留。Kokkos启动参数使用 neigh half newton on，与该Tersoff构建的要求配套。

3. 保存可直接提取的输出#

compute layers all chunk/atom bin/1d z lower .05 units reduced
fix profile all ave/chunk 10 100 1000 layers temp norm sample ave one file temperature.profile format " %.16g"
variable energy equal etotal
variable wh equal f_hotbath
variable wc equal f_coldbath
fix work all ave/time 100 1 100 v_energy v_wh v_wc file bath-work.dat format " %.16g"
run 50000
write_data final.data nocoeff

Si的8胞长度用20层，12/16胞用40层，避免分层比原子层间距还细而出现空层。每1 ps输出一个时间平均温度剖面，每0.1 ps保存热浴能量。自定义format中的前导空格保留列间隔，提取时才能直接分列。

提交与提取#

P100资源使用：

sbatch p100.slurm

V100 16G资源使用：

sbatch v100.slurm

选择与你权限匹配的一种即可。每个脚本包含8个数组任务，默认最多同时2个；按自己的并发权限调整。想先跑第0个Si小例子，可用 sbatch --array=0 p100.slurm。完整批次结束后：

module load ase/3.29.0
python3 collect.py --root work --output results
bash plot.sh --data-dir results

work和results保存自己的计算，已有目录时脚本拒绝覆盖。只跑一个例子时可先检查其目录；整批提取脚本要求8例两阶段都完成。

4. 从累计热浴功得到热流#

本例 hotbath 的累计量为负，对应向原子输入能量；coldbath为正，对应移除能量。先检查二者随时间的斜率，再取正的输入/输出率。

Qh = -slope(hotbath tally)
Qc =  slope(coldbath tally)
J = (Qh + Qc) / (4 A)

分母同时包含两浴斜率的平均，以及周期体系的两个传热方向。A为横截面积，不是体积。采用metal单位时，1 eV/(ps·Å²)=1.602176634×10¹³ W/m²。

这批追加段中，E_total+hotbath+coldbath的漂移保持在约百分之几eV量级，热源输入与冷源输出斜率接近。温度和能量收支都应看，单凭目标温度330/270 K不能判断已经获得可靠热流。

方法说明见 LAMMPS热导计算指南 与 Langevin热浴及能量累计。长度效应和表观热导提取继续看尺寸与采样比较。