the effect of Pdamp in fix npt on the equilibration density/energy for Tip4p water

Dear LAMMPS users,

I am doing a testing run on Tip4p water with kspace_style pppm/tip4p at T=298K and P=1bar via fix npt for N=216 water molecules from a initial density of 0.9872 gcm-3 (the orientation of water molecules in the initial configuration are completely artificial, i.e. I have all of them lied up in the same direction, but I assume for liquid water the system should equilibrate itself to natural orientation)

I am using a time step of 2.0fs.

For two different runs, both the one with Pdamp = 50 and the one with Pdamp=1000 has been able to equilibrate the total potential energy to around the same as literature values and with each other (~-10kcal/mol) for the first 400,0000 steps, and the three radial distribution functions also match well with reference radial distribution functions from literature and with each other.

However, strange thing starts to happen after ~400,0000 steps (~10ns), the run with Pdamp = 50 starts to equilibrate to lower potential energy
(~-12kcal/mol) with lower density/large volume and seem to find its new equilibrium there, while the run with Pdamp = 1000 remains the same, with the pressures of both runs are still oscillating around the desired P=1.

Hence I wonder whether the change for the run with Pdamp = 50 is due to the fact that Pdamp is too small? Because from LAMMPS manual, it says that:

“If Pdamp is too small, the pressure and volume can fluctuate wildly; if it is too large, the pressure will take a very long time to equilibrate. A good choice for many models is a Pdamp of around 1000 time steps”

However, in my case the pressure and volumes fluctuation for both runs are around the same, actually the pressure fluctuation for the run with Pdamp = 50 is smaller.

Below is my input if further information on the input setting is needed, otherwise please ignore it:
dimension 3
units real
atom_style full

read_data config-water.dat
replicate 1 1 1

group all type 1 2
group O type 1
group H type 2

neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes one 4000

pair_style lj/cut/tip4p/long 1 2 1 1 0.15 8.5 8.5
kspace_style pppm/tip4p 1.0e-5

pair_modify shift yes
pair_modify tail no

pair_coeff 1 1 0.155 3.1536
pair_coeff 1 2 0.0000 0.0000
pair_coeff 2 2 0.0000 0.0000

bond_style harmonic
angle_style harmonic
dihedral_style none
improper_style none

bond_coeff 1 10000.00 0.9572
angle_coeff 1 10000.0 104.52

velocity all create 298 427 dist gaussian mom yes

restart 1000 old_config_1.dat old_config_2.dat

fix 1 all shake 0.00001 60 0 b 1 a 1
fix 5 all npt temp 298.0 298.0 100 iso 1 1 1000 drag 0.0

thermo 1000

thermo_style custom step temp pe evdwl ecoul press vol etail
dump id all xyz 1000
dump_modify id element O H

timestep 2.0
run 7000000


I don’t know, but you might monitor the

center-of-mass motion of the system. If it

starts to drift over long times, you could
have issues with its velocity contributing
to the temperature of the system and effectively
lowering the kinetic energy (temperature)

of the molecules.