LAMMPS反向NEMD的温度剖面与误差

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

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

1 先建立两档热流强度#

选择 3：LAMMPS反向NEMD。体系为2048个Lennard–Jones原子，约化密度0.8、约化温度1，使用周期边界。每50或100步交换一次动能，每档四个初速度种子。本例采用LAMMPS 30Sep2026的P100 GPU包构建，编译环境见版本环境。

cat work/lammps-nemd-50-s27121/in.lammps
bash submit.sh

第一个种子的完整输入如下。所有温度、距离、时间、热流和热导率均为 LJ约化单位，不能直接把表中数字写成W/(m·K)。

units lj
atom_style atomic
boundary p p p
lattice fcc 0.8
region box block 0 8 0 8 0 8
create_box 1 box
create_atoms 1 box
mass 1 1
pair_style lj/cut 2.5
pair_coeff * * 1 1 2.5
pair_modify shift yes
neighbor 0.3 bin
neigh_modify every 1 delay 0 check yes
thermo 100
thermo_style custom step time atoms temp pe ke etotal press vol density
thermo_modify norm no format float %.12g
velocity all create 1 27121 mom yes dist gaussian
timestep .004
fix eq all nvt temp 1 1 .5
run 8000
unfix eq
reset_timestep 0
fix int all nve
compute keatom all ke/atom
variable T atom c_keatom/1.5
compute layers all chunk/atom bin/1d z lower .05 units reduced
fix profile all ave/chunk 10 100 1000 layers v_T file temperature.profile
fix flux all thermal/conductivity 50 z 20
variable transfer equal f_flux
fix exchange all ave/time 100 1 100 v_transfer file exchange.dat
dump tr all custom 10000 trajectory.lammpstrj id type x y z vx vy vz
dump_modify tr sort id
run 200000
write_data final.data

fix thermal/conductivity 50 z 20 将z方向划为20层，在第1层和第11层之间交换动能。把50改成100即得到另一档强度。该命令记录累计交换能量；它不同于每步直接添加固定热量。

2 输出有哪些、先看什么#

head -n 8 work/lammps-nemd-50-s27121/temperature.profile
tail work/lammps-nemd-50-s27121/exchange.dat
python analyze.py
python plot.py

自己的8个任务完成后，运行 python analyze.py --root work。temperature.profile 每1000步写一个20层剖面，本例共200个剖面；exchange.dat 是时间步和累计交换能量。先看后半段剖面是否稳定，再拟合热流和温度梯度。

3 周期边界为什么要除以2#

周期体系的热量从热层沿两个方向传向冷层。单个方向的热流为：

J = 累计交换能量对时间的斜率 / (2 × 横截面积)
κ = J / |温度梯度|

本例对后100000步平均，分别拟合 z/Lz=0.125–0.425 和 0.625–0.925 两段，避开交换层。两个梯度应方向相反；提取脚本对其绝对值取平均。

4 实际结果与解读#

交换更频繁时，温度梯度更明显。两档结果的区间存在重叠，但这只是一项小尺寸、有限时间的教学比较。自己的体系仍需增加尺寸、采样长度和交换间隔档位；若剖面明显弯曲，应减弱热流并重新检查线性响应。

官方命令与周期热流说明；若希望从平衡态相关函数开始，阅读Green–Kubo教程。