LAMMPS两温模型与电子能量交换

两温模型把原子用经典MD描述，把电子用空间网格上的温度场描述，再让两者交换能量。它不等于显式计算电子波函数；电子比热、密度、热导和耦合系数都需要与自己的材料及适用温度匹配。

下载三篇共用的完整输入、实算数组、原生科学输出和绘图包。下载后保存为 example.zip，在自己的目录解压：

python3 -m zipfile -e example.zip .
cd lammps-coupled-models

sbatch job.slurm

本包顺序执行12个小条件，三个主题共用同一次运行。CPU分区按56核整节点规则申请，实际用1核；CPU分区维护时等待恢复，不在登录节点自行跑生产计算。已分配的CPU作业内可 bash run.sh。使用现有LAMMPS22Jul2025 Update6 CPU、oneAPI2024.2、Intel MPI2021.13，分析绘图用sisl0.16.4；不需要重新安装公共软件。

输出 VALIDATION.json、完整 analysis-fields.npz 和 plots/。reference-results/保留本文原始数组与报告，当前目录保存自己的复跑结果。只重画时加载sisl后执行 python analyze.py、python plot.py，不重复计算。

1. 这个小模型包含什么#

32个原子、8 Å立方周期盒，原子质量63.546u；非键相互作用设为零，专门看热交换。电子使用4×4×4网格、比热0.001 eV/(electron·K)、密度1 electron/ų。这些电子参数是演示系数，不能因为原子质量接近Cu就当成Cu材料模型。

五个条件分别是：电子热扩散、无热扩散对照，以及同一0.1 ps耦合过程的三档时间步。后者电子初温1000 K、原子初温100 K，离子电子摩擦系数1000 u/ps，电子热导设0以单独观察交换。

2. 完整耦合输入#

cat inputs/ttm-coupled-fine/input.lammps
head -n 9 inputs/ttm-coupled-fine/initial-electrons.dat

units metal
atom_style atomic
boundary p p p
lattice fcc 4.0
region box block 0 2 0 2 0 2
create_box 1 box
create_atoms 1 box
mass 1 63.546
pair_style zero 3.0
pair_coeff * *
velocity all create 100 46717 mom yes rot no dist gaussian
fix integrator all nve
fix electrons all ttm 46719 0.001 1.0 0 1000 0.0 0.0 4 4 4 infile initial-electrons.dat outfile 1000 electrons.dat
timestep 1e-05
thermo 1
thermo_style custom step time temp pe ke etotal f_electrons[1] f_electrons[2]
thermo_modify norm no format float %.17g
dump atoms all custom 1 atoms.dump id type x y z vx vy vz fx fy fz
dump_modify atoms sort id format float %.17g
run 10000
write_data final.data nocoeff

fix electrons ... ttm的参数依次是随机种子、电子比热、电子密度、热导、电子离子摩擦、电子阻止项、速度阈值以及三维网格尺寸。摩擦单位是质量/时间，不是无量纲阻尼。原生程序要求电子离子摩擦严格大于0；近似解耦的扩散/不扩散对照用 1e-20，并在下面的解析比较中检查其影响足够小。

fix integrator all nve负责原子积分；TTM自己施加相应电子温度的摩擦和随机力。本例没有再加一套原子Langevin热浴，以免把两种控温混在一起。

初温文件每行是从1开始的 ix iy iz temperature。完整4³均匀初温如下：

# grid ix iy iz Te
1 1 1 1000
2 1 1 1000
3 1 1 1000
4 1 1 1000
1 2 1 1000
2 2 1 1000
3 2 1 1000
4 2 1 1000
1 3 1 1000
2 3 1 1000
3 3 1 1000
4 3 1 1000
1 4 1 1000
2 4 1 1000
3 4 1 1000
4 4 1 1000
1 1 2 1000
2 1 2 1000
3 1 2 1000
4 1 2 1000
1 2 2 1000
2 2 2 1000
3 2 2 1000
4 2 2 1000
1 3 2 1000
2 3 2 1000
3 3 2 1000
4 3 2 1000
1 4 2 1000
2 4 2 1000
3 4 2 1000
4 4 2 1000
1 1 3 1000
2 1 3 1000
3 1 3 1000
4 1 3 1000
1 2 3 1000
2 2 3 1000
3 2 3 1000
4 2 3 1000
1 3 3 1000
2 3 3 1000
3 3 3 1000
4 3 3 1000
1 4 3 1000
2 4 3 1000
3 4 3 1000
4 4 3 1000
1 1 4 1000
2 1 4 1000
3 1 4 1000
4 1 4 1000
1 2 4 1000
2 2 4 1000
3 2 4 1000
4 2 4 1000
1 3 4 1000
2 3 4 1000
3 3 4 1000
4 3 4 1000
1 4 4 1000
2 4 4 1000
3 4 4 1000
4 4 4 1000

3. 先用热扩散检查网格和单位#

ttm-diffusion的x向四个温度是1200、1000、800、1000 K，y/z方向重复；平均温度1000 K。热导0.1 eV/(ps·Å·K)，对应这里的热扩散系数100 Ų/ps。4个网格中心位于x=1、3、5、7 Å。

cat inputs/ttm-diffusion/input.lammps
head -n 9 inputs/ttm-diffusion/electrons.dat.100

该单一温度傅里叶模式在本例离散网格和显式步长下，每步振幅乘0.995，因此第n步的预期振幅为 200×0.995ⁿ K。脚本比较全部64个网格点；不是只拟合一条平均温度曲线。ttm-no-diffusion将热导设0，温度模式基本保持。这里比较的是已知离散算子的结果，不代表4³网格已达到连续介质空间收敛。

4. 读取两条能量通道#

原生 f_electrons[1] 是电子总能，单位eV；f_electrons[2] 是当前这一步传给原子的能量，单位eV，不是累计量。本例电子总能也可以由所有网格点独立计算：Ue=Ce×rho×每格体积×ΣTe；平均温度1000 K时为512 eV。

脚本另外从所有32个原子的质量与速度重算动能，从全部受力、速度和时间步重算每步 ΣF·v·dt，再核对电子能量变化。这些步骤帮助发现单位、输出字段或时间索引错误。

5. 不能忽略时间步#

固定同一0.1 ps，对比原子加电子步末总能的最大波动与实际交换能量之比：

本例用1%作为时间步比较标准，因此选择 1e-5 ps，即0.01 fs。前两档波动较大，结果也保留在参考报告中。速度Verlet积分使原子与电子步末能量带有半步错位，二者的和不能被解释为每一步严格恒定；该1%比较也不是所有材料的通用精度承诺。

原子温度沿轨迹上升并有明显涨落，32原子的短轨迹不代表平衡温度或宏观升温曲线。正式使用时需补真实原子势、材料电子参数、空间网格和时间步收敛。激光源、表面烧蚀及温度依赖参数属于其他TTM模式，不是这份常系数输入已经计算的内容。

参考：匹配版本TTM说明。