Water density computation using TIP4P

LAMMPS LAMMPS Beginners
1

Hi,
I am simulating bulk water using TIP4P/2005 model in NPT to obtain water density. The problem is that the density values are very low (around 0.007) in gr/cm3. Could you help me to solve the problem? input files are attached.
Thank you
input.txt (1.8 KB)
data.lmpdat (560.9 KB)

2

The group function count() is the number of atoms in the group.

Therefore, your variable dens equal count(all)/vol will print the density in atoms/Angs3. The canonical mass density can be accessed in the thermo output directly with the density keyword, as in

thermo_style   custom step temp press vol etotal enthalpy v_dens density
3

Thanks for the reply,
But I tried it first and the results are the same…

4

Very unlikely. But feel free to submit an input that reproduce the behavior you are describing together with your LAMMPS version.

5

You are starting from a very low density. Why not build an initial system that is closer to the target density?

When using ambient pressure it will take almost and eternity until you have an equilibrium. In general, it is often easier to expand a system than to shrink it.

You can accelerate the process by creating a more dense system from the beginning (see for example 8.4.4. TIP4P water model — LAMMPS documentation). You can also enforce a faster compression, either through initially using fix npt with a higher target pressure until your system is closer to the target density or through using fix nvt in combination with fix deform.

Please also note the “nreset” keyword of fix npt.

6

Thanks for your valuable reply.
When I use more water molecules with NPT, the initial pressure of the system will be larger and the system tends to expand (till the pressure reaches to 1 atm). So, the same density of 0.0007 g/cm3 will be obtained. I dont know what is the problem.

7

0.007 is the initial density of your sample.
With your settingsEDIT, I can run an NPT simulation of the attached, pre-equilibrated water box and the density is just fine:
water02.data.gz (59.7 KB)

# Atom Definition
read_data water02.data
mass 1 1.0079
mass 2 15.9994

If you want to stick to your setup, the only way to get liquid water is to actually compress the initial sample with fix deform. Here I use a loop as the PPPM will fail if you change the box vectors too much:

# Time integration
timestep        1.0
velocity        all create 300.0 4928459 mom yes rot yes dist gaussian #check gaussian or uniform
fix             1        all      nvt temp 300.0 300.0 100
fix             2        all      deform 1 x scale .95 y scale .95 z scale .95 units box remap x

# Shrink the sample to 1 g/cm3.
variable        i loop 32
label           cycle1
 run             2000
 next            i
jump            SELF cycle1
unfix           1
unfix           2

and then continue with an NPT simulation.

EDIT: changing the coupling constant of the barostat to $(1000*ts).

8

Thanks for your valuable replies.
But still a question which is elementary. The problem is to find the right density of liquid water in ambient condition. Is it possible with a NPT run? As you noted, the number of water molecules seems to be small to represents liquid water. But, with running NPT, the initial density is about 0.007 g/cm3 and at the same time the pressure is higher than 6 atm. Therefore, by proceeding the simulation, the pressure tends to reduce to 1.0 atm (as sets in the script) and consequently, the density being lower and it will not result the right value of 1.0 g/cm3. I mean, if the number of water molecules is considered small, the pressure should be lower than the set pressure of 1.0 atm and in the NPT run, the volume of box would shrink. But, we have the opposite trend now. Am I wrong with units of the density or the simulation setup?

9

The density is fine but the pressure fluctuates in an extremely wide range. This could be due to the fact that the density of liquids is not a strong function of the pressure.

10

There is no right density. The density you obtain is a function of:

  1. The force field model and parameters.
  2. The simulation settings (as you have discovered with your starting structure).
  3. The thermodynamic variables N, V, T, P.

No one said that. I have provided a smaller sample of 1000 molecules to make the point that your settings are fine to simulate liquid water at 300 K and 1 atm. A lot of properties of liquid water show a little dependence from the system size, as highlighted in this excellent presentation by @akohlmey.

As already pointed out, this happens because your initial guess is too far away from equilibrium conditions. What you are observing is a metastable gas which won’t condense into a liquid unless you provide a condensation seed (look for literature on how to determine melting/boiling temperatures with the two-phase method).

The settings and units are fine.

You have to quote the relevant commands to make any useful statement. In your original script, you used a coupling constant for the barostat which is way too small:

fix        1 all npt temp 300.0 300.0 100 iso 1.0 1.0 100

The recommended value is $(1000*ts).

1 Like
11

I got it, thanks for the time you spent for the reply.