LAMMPS SPC/E水、长程静电与SHAKE

本文使用鸢算科技集群 LAMMPS 30 Sep 2026 的真实计算。下载输入、参考轨迹、分析脚本与数据图（约 72.5 MB），解压后进入 lammps-examples。cases 保存参考结果，inputs 保存干净输入，使用 new-case.py 创建自己的计算副本。

要完成什么#

本例自行建立 216 个 SPC/E 水分子，共 648 原子。water.data 包含分子 ID、电荷、原子、键和角；atom_style full 用于保存这些信息。采用 real 单位、PPPM 长程静电和 8 Å 短程截断，先最小化，后 NVT 10 ps、NPT 30 ps。

氧/氢电荷分别为 −0.8476/+0.4238 e，总电荷为零。fix shake 约束 O–H 键和 H–O–H 角，动力学步长 1 fs。最小化阶段未施加 SHAKE，采用输入中的谐和键角稳定初态；进入动力学后约束恢复刚性几何。

这是长程静电版本的 SPC/E 教学模型，和原始有限截断模型并不完全等价。密度、压力和 RDF 可以用来检查建模，但 30 ps 不能证明扩散或介电性质收敛。更换水模型时要同时替换电荷、几何、LJ 和约束设定。

1 准备并提交#

python new-case.py water-spce
cd work/water-spce
sbatch run.sh

run.sh 是本算例匹配的核时 P100 或 V10016G 模板。套餐用户按自己的权限使用 yskit 生成提交脚本，保留下面的计算命令。完整模拟在计算节点执行。其他例子回到资料包根目录，替换目录名后分别提交。

2 实际计算命令与输入#

water-spce#

mpirun -np 1 lmp -sf gpu -pk gpu 1 -in in.lammps > stdout.log 2>&1

下面是 water-spce/in.lammps 的完整输入；其他目录的完整输入也在资料包中。

units real
atom_style full
boundary p p p
read_data water.data
pair_style lj/cut/coul/long 8.0
pair_coeff 1 1 .1553 3.166
pair_coeff 1 2 0 1
pair_coeff 2 2 0 1
bond_style harmonic
bond_coeff 1 1000 1
angle_style harmonic
angle_coeff 1 100 109.47
special_bonds lj/coul 0 0 0
kspace_style pppm 1e-4
neighbor 2 bin
neigh_modify every 1 delay 0 check yes
minimize 1e-8 1e-5 2000 10000
fix rigid all shake 1e-5 20 0 b 1 a 1
velocity all create 300 712341 mom yes rot yes dist gaussian
timestep 1
fix eq all nvt temp 300 300 100
thermo 100
thermo_style custom step time atoms temp pe ke etotal press vol density
thermo_modify norm no format float %.12g
run 10000
unfix eq
reset_timestep 0
fix int all npt temp 300 300 100 iso 1 1 1000
compute rdf all rdf 100 1 1
fix rdfout all ave/time 10 100 1000 c_rdf[*] file rdf.dat mode vector
dump tr all custom 200 trajectory.lammpstrj id type x y z vx vy vz
dump_modify tr sort id
run 30000
write_data final.data
write_restart final.restart

数据文件头和原子部分节选如下，完整拓扑使用附件，不要把节选当成完整文件运行。

Water topology

648 atoms
432 bonds
216 angles

2 atom types
1 bond types
1 angle types

0 18.64 xlo xhi
0 18.64 ylo yhi
0 18.64 zlo zhi

Masses

1 15.9994
2 1.008

Atoms # full

1 1 1 -0.8476 1.5533333333333335 1.5533333333333335 1.5533333333333335
2 1 2 0.4238 2.5533333333333337 1.5533333333333335 1.5533333333333335
3 1 2 0.4238 1.2200200857650962 2.4961494760650513 1.5533333333333335
4 2 1 -0.8476 4.66 1.5533333333333335 1.5533333333333335
5 2 2 0.4238 5.66 1.5533333333333335 1.5533333333333335
6 2 2 0.4238 4.326686752431763 2.4961494760650513 1.5533333333333335
7 3 1 -0.8476 7.766666666666667 1.5533333333333335 1.5533333333333335
8 3 2 0.4238 8.766666666666666 1.5533333333333335 1.5533333333333335
9 3 2 0.4238 7.433353419098429 2.4961494760650513 1.5533333333333335
10 4 1 -0.8476 10.873333333333335 1.5533333333333335 1.5533333333333335
11 4 2 0.4238 11.873333333333335 1.5533333333333335 1.5533333333333335
12 4 2 0.4238 10.540020085765098 2.4961494760650513 1.5533333333333335
13 5 1 -0.8476 13.979999999999999 1.5533333333333335 1.5533333333333335
14 5 2 0.4238 14.979999999999999 1.5533333333333335 1.5533333333333335

3 查看输出与实际数据#

tail -n 30 stdout.log
ls -lh

stdout.log 是计算输出；轨迹保存在 trajectory.lammpstrj（适用的算例），结构保存在 final.data。静态扫描和 run 0 对照不一定生成动力学轨迹。以下动态均值取最后一个热力学输出块的后半段；静态结果取对应输出，单位随输入确定。

LJ/珠簧例子使用约化单位；metal 为 K、g/cm³、eV；水的 real 单位为 K、g/cm³、kcal/mol。势能列是整个体系的总量，跨不同原子数不可直接相比。

4 提取数据与重画#

回到附件根目录，使用已经安装的 ASE Python 环境执行真实数据分析脚本：

module load ase/3.29.0
python analyze.py --root .
ls analysis
ls figures

脚本使用 NumPy/Matplotlib，读取参考日志、轨迹和相关函数，生成 CSV 与 PNG。查看 analyze.py 中本篇对应目录的处理段，可以把路径换成你的新结果。完整脚本按本专题的全部参考目录组织，单个新算例不需要假造其他目录来运行；提取对应段即可。

进一步查阅#

• LAMMPS专题目录
• 官方30 Sep 2026手册
• GPU与KOKKOS包的详细使用范围