LAMMPS eFF多电子波包与固定芯

eFF把电子表示为Gaussian波包；与单电子氢模型相比，多电子体系还要处理Pauli项和芯电子。本篇从固定结构的H2和一个Li模型入手，重点是数据怎样写、分项怎样读。

下载五篇共用的完整输入、科学数据与绘图包，保存为 example.zip，在自己的新工作目录中：

python3 -m zipfile -e example.zip .
cd lammps-electronic-polar-models

sbatch job.slurm

这五篇共用一次计算，顺序执行40个小条件。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/提供本文参考输出，inputs/提供全部结构、输入和计算后的原生日志。默认包校验输入；修改模型请另建目录、用自己的CPU作业脚本运行，保留原始数据。

1. 先准备核与电子角色#

cat inputs/H2-opposite/system.data
cat inputs/Li-fixedcore/system.data

H2完整结构：

Owned electronic-model input

4 atoms
2 atom types

-1000 1000 xlo xhi
-1000 1000 ylo yhi
-1000 1000 zlo zhi

Masses

1 1.000794
2 0.00054857990943000005

Atoms # electron

1 1 1 0 0 -0.69999999999999996 0 0
2 1 1 0 0 0.69999999999999996 0 0
3 2 0 1 1.5 -0.20000000000000001 0 0
4 2 0 -1 1.5 0.20000000000000001 0 0

atom_style electron的Atoms列是 id type q espin eradius x y z，单位为Bohr/Hartree/fs。espin=0为核，+1/−1为显式电子的离散标签，2为包含两个芯电子的固定Gaussian芯；它们与三分量SPIN方向不是同一种变量。电子的q列为0，物理电子电荷由角色确定，不能把q列简单求和当作物理净电荷。

这里两个核相距1.4Bohr，两个电子中心相距0.4Bohr，宽度均1.5Bohr。相反标签与相同标签两组使用相同固定位置、宽度；未做几何或电子宽度优化。

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

units electron
atom_style electron
atom_modify map array
boundary p p p
read_data system.data
pair_style eff/cut 40.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
thermo 1
thermo_style custom step pe temp press c_components[1] c_components[2] c_components[3] c_components[4]
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 fx fy fz c_packet[*] modify sort id format float %.17g

compute components all pair eff/cut输出四项：Gaussian内部动能、Pauli、电荷相互作用和半径限幅势。内部Gaussian动能已经属于有效势，不应再次额外加到总势能上。compute packet同时输出标签、宽度、径向速度和径向力。

改变标签后，内部动能和电荷项保持一致，Pauli项发生变化。这里比较的是同一未优化Gaussian几何的模型能量，不是H2真实单重态/三重态本征能隙。 相反标签的Pauli项也不必为0。

3. Li固定芯怎么写#

Owned electronic-model input

2 atoms
2 atom types

-1000 1000 xlo xhi
-1000 1000 ylo yhi
-1000 1000 zlo zhi

Masses

1 6.9400000000000004
2 0.00054857990943000005

Atoms # electron

1 1 3 2 0.69999999999999996 0 0 0
2 2 0 1 2 0.20000000000000001 0 0

Li芯的q=3、espin=2、宽度0.7Bohr，另有一个显式价电子，宽度2Bohr、x=0.2Bohr。物理净电荷为 3−2个隐式芯电子−1个显式电子=0，原生q列的和则为3。固定芯宽度是本例输入参数，没有参与优化。

实际有效势−4.947452358974169Ha，四个分项按相同方式提取；固定芯模型的绝对能量参照不可直接与不同全电子或ECP模型混排。图中的密度切片由输入Gaussian宽度推导，不是DFT网格。

4. 径向力与空间力怎样检查#

cat inputs/H2-opposite/final.dump
cat inputs/Li-fixedcore/input.lammps
cat inputs/Li-fixedcore/final.dump

电子宽度的力在 c_packet[4]，位置力在fx/fy/fz。包内对H2两组分别把一个电子宽度改变±1e−5Bohr，对Li价电子分别改变x和宽度，再实际计算两侧单点，以 −[E(+h)−E(−h)]/(2h) 对照对应原生力。它能发现把径向力当作x方向力、或把标签列读错的问题。

已有原生结果时，只重新提取和画图：

module load sisl/0.16.4
python analyze.py
python plot.py

这一步读取实际结果，不重新运行原生模拟。

下一步可看ECP的Si/C两个形式。模型和输入语法参考LAMMPS eFF说明。