NVT Simulations not giving the correct pressure?

Hey All,

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

  1. 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)
  2. 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?

Then some parameters could explain the difference you obtain:

  • Cut-off : 8.5 A is quite small, are you sure about the value?
  • Box size : how many molecules are you using?

I am using 360 molecules in a cubic box of size 22.157 A.

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.

Okay. But how do we then decide the sufficiently large box size?

Is it like I have to try different systems and find a large enough one after which the average pressure values do not change?

For now I was following the thumb rule that the cut-off radius be ~2.5sigma and the box length >= 2cut-off radius.

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.

Thank you so much for the help.


Have you tried with the epsilon value from the LAMMPS manual?
See my comment to your post in the other thread.

I did see your comment. I am trying with that as well.

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.

Thank you so much @srtee e