Packmol接LAMMPS：显式拓扑与刚性水MD

下载完整建模、拓扑、运行和绘图脚本。解压后进入water-md子目录。

1 使用同一份Packmol坐标#

本篇承接Packmol接GROMACS，使用相同216水分子装箱，但为LAMMPS明确生成电荷、键和角拓扑。LAMMPS 30Sep2026的8173m CPU模块采用Intel oneAPI 2024.2编译。

先运行上一例的prepare.py、Packmol和convert.py。然后：

module purge
module load gpumdkit/1.5.8
python lammps_data.py
head -25 water.data

生成648个原子、432条O–H键、216个H–O–H角。Atoms # full每行依次为原子ID、分子ID、类型、电荷、x、y、z；类型1为O，类型2为H。一个分子的三个原子使用相同分子ID。

本例不靠坐标推断分子身份：输入按OHH逐分子排列，生成器据此写两条键和一个角。换成其他分子必须重新定义拓扑，不能沿用该生成器。

2 完整LAMMPS输入#

in.water：

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 0.1554 3.16557
pair_coeff 1 2 0.0 1.0
pair_coeff 2 2 0.0 1.0
bond_style zero
bond_coeff 1 1.0
angle_style zero
angle_coeff 1 109.47
special_bonds lj/coul 0.0 0.0 0.0
kspace_style pppm 1.0e-5
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
fix rigid all shake 0.0001 20 0 b 1 a 1
thermo_style custom step temp press pe etotal density
thermo 100
minimize 0.0 1.0e-4 500 5000
reset_timestep 0
velocity all create 300.0 27129 mom yes rot yes dist gaussian
fix thermostat all nvt temp 300.0 300.0 100.0
timestep 1.0
variable temp_now equal temp
variable press_now equal press
variable pe_now equal pe
fix output all ave/time 10 1 10 v_temp_now v_press_now v_pe_now file lammps-thermo.dat
dump traj all custom 100 water.dump id mol type x y z
dump_modify traj sort id
run 20000
write_data water-final.data

units real使用Å、fs、kcal/mol和atm。LJ参数明确给在输入中；长程静电使用PPPM，实空间截断8 Å。fix shake维持两个O–H键和H–O–H角，避免把水内振动作为自由度积分。bond_style zero、angle_style zero提供约束几何，不能删掉SHAKE后仍视为SPC/E刚性水。

3 实际运行#

module purge
module load lammps/30Sep2026-8173m
export OMP_NUM_THREADS=1
lmp -in in.water -log lammps-water.log

先最小化，再生成300 K初速度，NVT运行20000步，步长1 fs，共20 ps。lammps-thermo.dat每10步输出温度、压力和势能；water.dump保存坐标，water-final.data保存最终结构和拓扑。

4 核对结果和单位#

本次10–20 ps温度平均 299.01 K，时间序列标准差11.68 K。压力换算为bar后平均 397.45 bar、标准差850.42 bar；最终O–H距离偏离1 Å的最大值约 9.6×10⁻⁷ Å。

图中压力仍用LAMMPS原生atm、势能用kcal/mol。换算到GROMACS常用单位时，1 atm=1.01325 bar，1 kcal/mol=4.184 kJ/mol。两套例子采用不同恒温器、随机初速度和部分数值处理，不能逐时刻比较轨迹，也没有把平均势能差当作两个程序精度的验证。

5 一起生成两套曲线#

运行完整附件后：

module purge
module load gpumdkit/1.5.8
python analyze.py
cat VALIDATION.json

脚本读取两套真实能量输出，检查20 ps结束、温度数据有限、最小化记录和O–H约束，并生成两张曲线与初始结构图。温度接近目标只是初步平衡的一个检查；此例没有完成NPT，也未给出独立统计误差、输运性质或可推广的平衡密度。

参考：LAMMPS官方SPC/SPC/E模型与SHAKE示例。