Strange phenomena during simulating water freezing in model TIP4P/2005

LAMMPS LAMMPS Beginners
1

Dear all

I am a beginner of lammps and I am trying to simulate the freezing process of water with TIP4P/2005 model.

Firstly, I tried to build a system with half water and half ice. After that I used nvt to control the temperature under its melting point. However the ice started to collapse and began to act weirdly since the nvt process have started. I am trying to find out what happened but have no clue. Maybe I have made a stupid mistake. Please help me with it.

here is my input file:

dimension 3
boundary p p p
units real
atom_style full

neighbor 3.0 bin
neigh_modify every 1 delay 0 check yes

pair_style lj/cut/tip4p/long 1 2 1 1 0.1546 10.0 8.5
bond_style harmonic
angle_style harmonic
kspace_style pppm/tip4p 0.0001

read_data tip3p.data

pair_coeff 1 1 0.1852 3.1589
pair_coeff 2 2 0.0 0.0
pair_coeff 1 2 0.0 0.0
bond_coeff 1 0.0 0.9572
angle_coeff 1 55.0 104.52

region WATER block INF INF INF INF INF 0.0 units box
group water region WATER
group ice subtract all water #all - water

timestep 1.0
velocity all create 100.0 123456 dist gaussian

thermo 1000
thermo_style custom step temp press etotal

dump 1 all custom 10000 testing-1.xyz id type x y z

fix 1 all shake 0.0001 20 0 b 1 a 1
fix 2 all npt temp 100.0 100.0 100.0 aniso 0.0 0.0 100.0
run 10000000

unfix 2

fix 3 water nvt temp 400.0 400.0 100.0
fix 4 ice spring/self 10.0
run 20000

unfix 3
unfix 4

fix 5 all recenter 0.0 0.0 0.0

fix ice6 ice nvt temp 100.0 100.0 100.0
fix water6 water nvt temp 100.0 100.0 100.0
run 15000000

here is part of my data file

576 atoms
384 bonds
192 angles
  0 dihedrals
  0 impropers

2 atom types
1 bond types
1 angle types

-7.688405550     7.751594450 xlo xhi
-7.762002600     6.737997400 ylo yhi

-13.373225400 13.316774600 zlo zhi

Masses

1 15.999400 # o*
2 1.007970 # h*

Atoms # full

  1      1   1 -0.820000    -6.298405550    -5.362002600    -9.863225400   0   0   0 # o*
  2      1   2  0.410000    -7.188405550    -5.712002600    -9.893225400   0   0   0 # h*
  3      1   2  0.410000    -5.878405550    -5.712002600   -10.653225400   0   0   0 # h*
  4      1   1 -0.820000    -2.478405550    -2.802002600   -12.153225400   0   0   0 # o*
  5      1   2  0.410000    -3.368405550    -2.462002600   -12.123225400   0   0   0 # h*
  6      1   2  0.410000    -2.058405550    -2.462002600   -11.363225400   0   0   0 # h*

Thanks and best wishes to you all

RX

2

If you want to ask questions about LAMMPS, you should post in the LAMMPS category, not “Materials Project”.

3

I doubt anybody will be able to help you if you remain vague on the issue.

Also your protocol is very exotic to say the least. You cool down the system to 100 K, then heat up the water for a very short time at 400K while not thermalizing the ice, and then make a long run at 100 K. Is this intended?

4

Why increase the neighbor list skin?

Why the choice of cutoffs (10.0 for lj/cut and 8.5 for coul/long)?

This is too coarse a kspace convergence for variable cell. Pressure is not converged at 1.0e-4, you probably need 1.0e-6. Only forces are acceptable at 1.0e-4 due to error cancellation.

What is the purpose of such a long equilibration at 100K? A simpler approach would be to start from a proper Ice geometry as data file directly and just gently bring that close to the melting point. You may not even need to adjust the box dimensions for that.

This is pretty pointless, if those atoms are not subject to time integration.

Unlike the previous equilibration, this seems far too short to get anywhere near a properly equilibrated liquid water system after it had been cooled.

This should not be needed. If you have a drift of the system (a so-called flying ice cube), then fix recenter will not cure that problem, only hide it and your thermostat will actually condition the system to the wrong temperature. If you follow the proper protocol for setting up a coexistence simulation by starting from the properly equilibrated frozen structure, there should be no drift forming. There can be a small net momentum due to the finite accuracy of assigned initial velocities, but that can be removed with the velocity command, once the system has been equilibrated to a proper bulk ice.

This makes no sense at all. First, you want a single, global thermostat and second you want to be closer to the melting point. Initially you want to equilibrate the system to the exact melting point according to the water potential used and only later change it to see the frozen part growing.

This is far too small to yield any meaningful results. You want a much larger system to be free from most finite size effects.

I and others have discussed how to properly set up coexistence simulations multiple times in the past and thus you can look up the complete protocol in the archives section. Since you are new to LAMMPS, my recommendation is to forget about melting and freezing and first try to reproduce bulk liquid and frozen water systems individually with the chosen potential and match the published results from the corresponding publication(s). That will also allow you to determine the volume at the melting point (so no more need to use the difficult to control NPT integrator, but you can jump ahead to use the simpler NVT setup at the expected density) and confirm the most efficient and suitable settings for simulation parameters like timestep, neighborlist skin and updates, shake convergence etc.

Depending on your level of experience, it may be also worth doing a coexistence simulation for a metal or LJ system first, since those are easier to set up and control than molecular systems, especially water with its strong directional interactions due to the hydrogen bonds.

5

I am not sure if it is my problem of description. I intended to heat part of the ice up to water and remain the rest of the ice. As a result this is a system with half water and half ice in my opinion.
Anyway, I really appreciate with your help and I will double check my protocol again.

6

Thank you for you prompt response!
I will try to start over with the simple ones.
Further more, I really appreciate your useful and pertinent suggestions.

7

Have you solved this problem？I meet the similar problem of collapsed ice. if you solved this problem,can you show me your input script as a reference.
thank you

8

Sorry，i didn’t use TIP4P/2005 model anymore. Hope you can find the solution.
Best wishes!