LAMMPS长程静电与精度选择
实际比较Ewald、PPPM、stagger、cg和MSM，检查容差、显式网格、总能与逐原子力。
阅读本篇目录
带电周期体系通常需要长程静电算法。本篇固定一个16个点电荷的小周期模型，用相同几何比较Ewald、PPPM及相关路径，再检查手动网格与求力容差对实际结果的影响。
下载12组长程静电完整输入、原始能量与力、CSV及绘图工具，解压进入run.sh目录。
1 提交并读取结果#
sbatch submit.slurm
模板为核时P100分区1GPU1CPU，顺序完成12个固定结构单点；使用P100构建中的CPU执行路径，没有GPU后缀，不作GPU性能结论。LAMMPS30Sep2026、GNU C++13.3.0/C++17、Open MPI5.0.10、双精度MKL FFT；分析ASE3.29.0。
套餐用户用yskit生成所属分区脚本，在分配资源中执行bash run.sh；不要在登录节点模拟。inputs为新输入，reference为本文实算输出，work为新结果，已有work请重新解压。脚本完成后自动提取和画图。
2 模型与完整PPPM输入#
8Å立方盒中放入8个+1e和8个−1e电荷，采用CsCl型两子格排布，并把1号原子移动(0.12,0.03,0.04)Å，使受力不再完全由对称性抵消。只计算固定结构的静电能，不包含真实晶体所需的短程排斥或完整材料势。
units real
atom_style charge
atom_modify map array
boundary p p p
region box block 0 8 0 8 0 8 units box
create_box 2 box
mass 1 23.0
mass 2 35.45
create_atoms 1 single 0 0 0 units box
create_atoms 2 single 2 2 2 units box
create_atoms 1 single 0 0 4 units box
create_atoms 2 single 2 2 6 units box
create_atoms 1 single 0 4 0 units box
create_atoms 2 single 2 6 2 units box
create_atoms 1 single 0 4 4 units box
create_atoms 2 single 2 6 6 units box
create_atoms 1 single 4 0 0 units box
create_atoms 2 single 6 2 2 units box
create_atoms 1 single 4 0 4 units box
create_atoms 2 single 6 2 6 units box
create_atoms 1 single 4 4 0 units box
create_atoms 2 single 6 6 2 units box
create_atoms 1 single 4 4 4 units box
create_atoms 2 single 6 6 6 units box
set type 1 charge 1.0
set type 2 charge -1.0
group probe id 1
displace_atoms probe move 0.12 0.03 0.04 units box
pair_style coul/long 3.8
pair_coeff * *
pair_modify table 0
kspace_style pppm 1e-6
neighbor 0.3 bin
thermo_style custom step atoms pe ecoul elong press
thermo_modify norm no format float %.16g
run 0
write_data structure.data
write_dump all custom atoms.dump id type q x y z fx fy fz modify sort id format float %.16g
pair_style coul/long处理与长程算法配套的实空间项，kspace_style处理倒空间部分，两者共同组成静电能。thermo_modify norm no使输出为整盒能量，没有除以16个原子。
pair_modify table 0关闭实空间库仑查表。本例要对照较紧的力容差，先消除默认有限位数表带来的额外精度限制；不能只把kspace容差写得很小，而忽略其他数值设置。
3 本包八个周期条件#
| 目录 | 设置 |
|---|---|
| bulk-ewald-1e10 | Ewald，目标相对力精度1E-10，作为同一模型数值参照 |
| bulk-pppm-1e4/1e6/1e8 | PPPM，三种目标相对力精度，自动选择网格 |
| bulk-pppm-mesh12 | PPPM 1E-6，但显式指定12³网格、5阶映射 |
| bulk-pppm-stagger | 交错网格PPPM/stagger，1E-6 |
| bulk-pppm-cg | PPPM/cg，1E-6，本例所有粒子均带非零电荷 |
| bulk-msm | coul/msm配MSM，1E-6，保留3.8Å实空间截断 |
MSM目录同时使用：
pair_style coul/msm 3.8
pair_coeff * *
pair_modify table 0
kspace_style msm 1e-6
kspace_modify cutoff/adjust no
MSM默认可能为平衡网格工作量而调整实空间截断；对很小的周期盒，增长后的邻居列表会明显变大。本例显式保留既定截断，让网格处理精度要求。更大体系应重新权衡，不把这个设置视为通用性能最优值。
4 在终端查看原生精度和能量#
grep -E 'grid =|estimated relative force accuracy|G vector|Step.*Atoms.*PotEng' work/cases/bulk-pppm-1e6/log.lammps
grep -E 'grid =|estimated relative force accuracy' work/cases/bulk-pppm-mesh12/log.lammps
cat work/analysis/kspace-comparison.csv
| 目录 | 整盒能量/kcal·mol⁻¹ | 相对Ewald能量差/kcal·mol⁻¹ | 力分量RMS差/kcal·mol⁻¹·Å⁻¹ |
|---|---|---|---|
| bulk-ewald-1e10 | -1351.762049818 | 0.000E+00 | 0.000E+00 |
| bulk-pppm-1e4 | -1351.832437089 | -7.039E-02 | 1.910E-03 |
| bulk-pppm-1e6 | -1351.759300522 | 2.749E-03 | 4.321E-05 |
| bulk-pppm-1e8 | -1351.762038436 | 1.138E-05 | 5.164E-07 |
| bulk-pppm-mesh12 | -1351.248183916 | 5.139E-01 | 4.164E-03 |
| bulk-pppm-stagger | -1351.765261001 | -3.211E-03 | 4.625E-05 |
| bulk-pppm-cg | -1351.759300522 | 2.749E-03 | 4.321E-05 |
| bulk-msm | -1351.762021051 | 2.877E-05 | 1.192E-04 |
自动PPPM的1E-4→1E-6→1E-8对照中，力差随精度要求收紧而降低；最紧条件的RMS差约5.16E-7 kcal/mol/Å。不同算法和参数的能量差不一定按同样规律变化，kspace_style的accuracy是目标相对力误差，不是总能收敛阈值。
显式12³网格条件保留在表中作为参数对照。它的原生估计相对力误差约1.57E-4，大于输入的1E-6目标；手动固定过粗的网格会改变自动精度选择，不能只看输入中的容差数字。
PPPM/cg在本模型中与普通PPPM结果相同，不代表任意稀疏电荷体系都无差别。MSM本例总力残余约7.73E-4 kcal/mol/Å，保留这一有限网格误差；用于敏感响应时继续收敛网格与参数，不把输出正常结束当作所有力都已达到需求。
提取与重画#
module load ase/3.29.0
python analyze.py work
python plot_results.py work
python check_results.py work
cat work/analysis/kspace-comparison.csv
ls work/figures
CSV列出方法、目标容差、显式修正、整盒总能、相对同几何Ewald参照的能量差、全部3N个笛卡尔力分量的RMS差、原生估计精度及总力残余。查看本文参考数据时将work改为reference。
参数：kspace_style、kspace_modify、库仑pair样式。薄层边界继续看slab修正实操。