I am writing to clear a doubt regarding NPT and NVT simulation in lammps, which is as follows:

I take an initial configuration of randomly packed water molecules (TIP4P/Ice) in a certain volume. The density that I use for packing is calculated from an EOS derived from MC simulations of liquid water. The density is 0.9888 g/cm^3. (for T = 280K and P=25Bar)

I Use this configuration to run an NVT simulation and calculate average pressure of the system. But the pressure shows to be in the negative.

I am using: compute P all pressure T
fix PAVG all ave/time 500 1000 500000 c_P file Pavg.dat
to get the average pressure values.

Can someone please explain why?

The input file looks like this:

boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic ### section1 ok
# remove hybrid if not necessary
pair_style hybrid lj/cut/tip4p/long 1 2 1 1 0.1577 8.5
pair_modify tail yes
kspace_style pppm/tip4p 1.0e-5 ### section2 ok
read_data data_D0.9888_random.lmp ### correct
pair_coeff 1 1 lj/cut/tip4p/long 0.210831 3.166800 # O O
pair_coeff 1 2 lj/cut/tip4p/long 0.000000 1.000000 # O H
pair_coeff 2 2 lj/cut/tip4p/long 0.000000 1.000000 # H H ##correct
minimize 1.0e-4 1.0e-6 100 1000 ## Yes required since starting from a random configuration
reset_timestep 0
fix SHAKE all shake 0.0001 20 0 b 1 a 1 ## bonds and angles constrained
neighbor 2.0 bin ## Okay
neigh_modify delay 0 every 1 check yes
timestep 1.0 ##Correct
variable TK equal 280
variable P equal 24.673 ## T and P correct, P only for knowledge
velocity all create ${TK} 12345 ## correct
fix TVSTAT all nvt temp ${TK} ${TK} 100 ## Thermostat
thermo 1000
thermo_style custom step etotal ke pe evdwl ecoul elong temp press vol density
run 1000000 ## Equilibration run
reset_timestep 0 ## Timestep = 0
compute T all temp
**compute P all pressure T ## Defining compute to get Average Pressures**
**fix PAVG all ave/time 500 1000 500000 c_P file Pavg.dat ## Averaging Pressure**
restart 1000000 Restart*.dat ## Writing restarts after every 1e6 timesteps
run 5000000 ## Production run 5ns

Each water model has a unique pressure-density relation. If the density you found was not calibrated using the exact same water model, then I would say that what you observe was expected. For liquid water, a tiny change in volume can lead to huge pressure difference as it is almost incompressible.

If you want to set a particular pressure instead of a density, just use NPT instead of NVT.

I understand this. But the EOS was indeed calibrated for TIP4P/Ice.

Fundamentally then, I should get the same average pressure right? I am asking this from a conceptual point of view, rather than a way out. I want to understand why am I getting negative pressure?

This is definitely a small system, you are at risk of finite size effect. You can try different box size and different cut-off. It could be that for larger box and larger cut off the (average) pressure will tend toward a different value.

Yes exactly. Its a convergence test. For bulk water using tip4p/epsilon, I see a clear convergence only for box size larger than 1000 molecules, and cut-off larger than 12 A.

Also please note that long-range treatment – both the choice of LJ cutoff and the use / non-use of long range electrostatics (and LJ tail corrections) – are part of a force field, not just the “obvious” parameters like \sigma and \epsilon, and especially if you’re interested in pressure. Different cutoffs give different results (unless proven otherwise).

Having said that, I just noticed that in your script you use tail corrections, which mean in theory that your results shouldn’t depend too much on cutoff. It’s still worth checking, and it’s also worth checking if the model’s original paper used tail corrections.

Oh okay I see. For the LJ interactions, the cut-off does not seem to matter as much since the LJ interactions do infact give good convergence at 2.5 times the sigma value.