LAMMPS eFF电子波包与半径优化

eFF使用可移动、可改变宽度的Gaussian电子波包，并让它与核相互作用。它与普通原子力场、DFT及经典自旋方向模型都不同。本例从一个核加一个电子开始，使用可以独立写出解析解的氢原子Gaussian模型。

下载两篇共用的完整输入、实算数据和绘图包，保存为 example.zip，在自己的工作目录解压：

python3 -m zipfile -e example.zip .
cd lammps-eff-hydrogen

sbatch job.slurm

CPU分区按56核整节点申请，本例顺序实际用1核。CPU分区维护时等待恢复，不在登录节点自行跑生产计算；在已分配CPU作业内可 bash run.sh。现有环境为LAMMPS22Jul2025 Update6 CPU、oneAPI2024.2/Intel MPI2021.13，分析绘图用sisl0.16.4。

两篇共用一次运行，包含14个小条件。输出 VALIDATION.json、完整 analysis-fields.npz 和 plots/；reference-results/保留本文原始科学结果。默认输入带校验，修改设置请复制至新的计算目录、使用自己的CPU提交脚本，不覆盖已完成数据。

1. 先认清电子数据格式#

cat inputs/radius-1p5/hydrogen.data
cat inputs/radius-1p5/input.lammps

完整结构：

Owned H nucleus and one Gaussian electronic packet

2 atoms
2 atom types

-50 50 xlo xhi
-50 50 ylo yhi
-50 50 zlo zhi

Masses

1 1.000794
2 0.00054857990943

Atoms # electron

1 1 1 0 0 0 0 0
2 2 0 1 1.5 0 0 0

atom_style electron各列是 id type q espin eradius x y z。本例核的espin=0、q=1；电子的espin=+1、q列为0，角色本身表示一个电子。不能将这个q列直接求和当成物理净电荷：本例物理净电荷是一个核电荷减一个电子，即0。

电子espin=+1/−1是eFF的离散自旋标签，不是SPIN进动中的三分量方向。只有一个电子时没有电子间Pauli项，所以本次正负标签对照能量一致；不能据此推断多电子体系也与自旋无关。

单位是Bohr、Hartree和fs。核的示例输入质量1.000794u，电子质量0.00054857990943u；核与电子中心均保持在原点，质量不影响这里的静态能量扫描。

2. 完整单点输入与输出分项#

units electron
atom_style electron
atom_modify map array
boundary p p p
newton on
read_data hydrogen.data
pair_style eff/cut 20.0 
pair_coeff * *
comm_modify vel yes
compute components all pair eff/cut
compute packet all property/atom espin eradius ervel erforce
compute effectiveT all temp/eff
compute effectiveP all pressure effectiveT
compute kinetic all ke/eff
thermo 1
thermo_style custom step time pe ke etotal temp press c_components[1] c_components[2] c_components[3] c_components[4] c_kinetic
thermo_modify norm no temp effectiveT press effectiveP format float %.17g
run 0
write_dump all custom final.dump id type q x y z vx vy vz fx fy fz c_packet[*] modify sort id format float %.17g
write_data final.data nocoeff

pair_style eff/cut 20指定20 Bohr相互作用截断；100 Bohr盒子满足本例的最小镜像要求。compute components all pair eff/cut依次给出Gaussian内部动能、Pauli、电荷相互作用和半径限幅势。本例单电子Pauli为0，未启用限幅时第四项为0。

compute packet ...输出离散spin、电子半径、径向速度和径向力，后者是宽度自由度的力，不是某一条x/y/z机械力。

3. 实际能量和径向力#

中心宽度s的解析有效势是 U(s)=3/(2s²)−2√(2/π)/s，对应径向力 −dU/ds=3/s³−2√(2/π)/s²。

脚本还在s=1.5两侧实际计算±1e−4 Bohr的两个单点，用能量中心差分检查径向力；不是只比较同一份文件的重复字段。完整原子编号、角色、中心位置和四个能量分项也同时核对。

4. 让原生最小化优化宽度#

cat inputs/radius-minimize/input.lammps
cat inputs/radius-minimize/final.dump

这个条件从s=1.5开始，使用 min_style cg、min_modify dmax 0.1 和 minimize 0 1e-10 1000 10000。实际31次迭代后，得到宽度 1.879971203176 Bohr、有效势 -0.424413181577 Hartree；解析最优值是1.87997120597325 Bohr和−0.424413181578388 Hartree。

日志的停止原因是 linesearch alpha is zero，输入的1e−10不能直接当成已经达到的残余力。本例另用最终宽度、能量和径向力与解析值比较；检查自己的任务时，也应同时看停止原因与实际输出。

图中密度由Gaussian模型推导，左图取实际优化宽度，右图取解析最优宽度，二者在本图显示精度内一致；它不是DFT密度网格。这个单Gaussian变分模型不等于精确氢波函数，也不提供氢原子的真实光谱。

若要看宽度随时间的实际变化，继续径向动力学与能量读取。参考：eFF模型说明。