TIP4P/Ice (2005) force field for simulation of ice 1h

Ankit_Mishra · May 12, 2017, 8:47am

Dear lammps Users,

I am trying to simulate ice 1h structure which is hexagonal in nature using TIP4P/Ice (2005) force field. My objective is to obtain phonon density of states of water at T = 100 K (though temperature does not effect DOS, but still I am doing at low temperature to avoid diffusion effects) using a force field which properly describes this phase. And compare it with the behaviour of water at similar T & P using a reaxFF implementation for a system, which also has water in it. This way I can compare how good the reaxFF implementation is.

Since ice 1h is the most stable phase at my required condition of temperature and pressure, hence I am trying to obtain its phonon DOS.

I have constructed the input script based on the information related to charges, and LJ parameters on the lammps web site (attached with this email).

But after a trying a lot of things, I am still unable to obtain a stable structure at even 10K.

So here I list all the problems I encountered and solutions I tried:

I encountered the first problem in kspace_style, since TIP4P/Ice require kspace_style pppm/tip4p for coloumb interaction calculations, and this style only works with orthogonal strcutures.

So I changed my hexagonal structure to a orthogonal one, by some basic transformations.

Then I tried to simulate this orthogonal structure and minimize it using style CG/hftn but both did not work, and the system blew up in less than 50 MD steps with a time step of 1 fs.
Then I tried other ways to minimize it, using suggestions on lammps page, i.e. putting viscous damping using fix viscous, and fix nve/limit. I tried both of them alone and combined initially, to obtain a stable system. But even that did not work, as system continuously increases in temperature and results in blow up.
I have also tried making initial time step extremely small, so as to make integration very fine, but even that does not control system temperature.

I really do not know now what to do now.

I think there can be some problem with my data file but I do not know, what specific problem. As structure is okay on visualization. There are no free H or O atoms in the system. As everything is bonded. There are 1152 water molecules in my system, and the total number of atoms is 3456. Also I have maintained the order of O atom id < H atom id and there arrangement in the data file.

So I do not know if there is anything wrong with it, as I have tried my best to avoid any mistakes in the data file which I am sending to everybody’s scrutiny, to save your time and effort. Though any suggestion related to possible problems in it will be really helpful.

Can someone recommend me what to do based on my objective and problems that I encountered.

I am attaching my input and data file for clarity.

I will be really grateful for suggestions.

Thanks,
Ankit

in.H2O (2.72 KB)

data.water (203 KB)

akohlmey · May 12, 2017, 12:01pm

Dear lammps Users,

I am trying to simulate ice 1h structure which is hexagonal in nature
using TIP4P/Ice (2005) force field. My objective is to obtain phonon
density of states of water at T = 100 K (though temperature does not effect
DOS, but still I am doing at low temperature to avoid diffusion effects)
using a force field which properly describes this phase. And compare it
with the behaviour of water at similar T & P using a reaxFF implementation
for a system, which also has water in it. This way I can compare how good
the reaxFF implementation is.

Since ice 1h is the most stable phase at my required condition of
temperature and pressure, hence I am trying to obtain its phonon DOS.

I have constructed the input script based on the information related to
charges, and LJ parameters on the lammps web site (attached with this
email).

But after a trying a lot of things, I am still unable to obtain a stable
structure at even 10K.

So here I list all the problems I encountered and solutions I tried:

1. I encountered the first problem in kspace_style, since TIP4P/Ice
require kspace_style pppm/tip4p for coloumb interaction calculations, and
this style only works with orthogonal strcutures.

as of LAMMPS version 6 Jan 2017 this statement is incorrect.

So I changed my hexagonal structure to a orthogonal one, by some basic
transformations.

it is generally more convenient and more efficient to run simulations in
orthogonal cells.

2. Then I tried to simulate this orthogonal structure and minimize it
using style CG/hftn but both did not work, and the system blew up in less
than 50 MD steps with a time step of 1 fs.

you define bonds and angle with force constants of 0.0. no surprise you
are getting garbage. you must use either fix shake or have suitable values
there (or larger, if you want to emulate fix shake) for scenarios in which
shake is not supported.

3. Then I tried other ways to minimize it, using suggestions on lammps
page, i.e. putting viscous damping using fix viscous, and fix nve/limit. I
tried both of them alone and combined initially, to obtain a stable system.
But even that did not work, as system continuously increases in temperature
and results in blow up.

see above, that is just nonsense.

4. I have also tried making initial time step extremely small, so as to
make integration very fine, but even that does not control system
temperature.

besides, why should the temperature be "constant" when your system isn't
equilibrated yet? this looks a lot like you need to spend more time
studying MD basics.

i've pared down your massive and needlessly complex and convoluted input
and it seems to be running fine in the NVT ensemble from the get go.

units real
boundary p p p
atom_style full
read_data data.water

mass 1 15.9994 # O
mass 2 1.008 # H

pair_style lj/cut/tip4p/long 1 2 1 1 0.1577 12.0
kspace_style pppm/tip4p 1.0e-5
pair_modify tail yes

pair_coeff 1 1 0.21084 3.1668
pair_coeff 1 2 0.0 0.0
pair_coeff 2 2 0.0 0.0

bond_style harmonic
bond_coeff 1 0.0 0.9572

angle_style harmonic
angle_coeff 1 0.0 104.52

neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes

timestep 1.0
velocity all create 10 34455 dist gaussian mom yes rot yes
fix 1 all nvt temp 10 10 100
fix 4 all shake 0.0001 10 1000 b 1 a 1

thermo_style custom step temp pe ke etotal press
thermo 100
run 1000

$ mpirun -np 4 ~/compile/lammps-icms/src/lmp_omp -in in.H2O -sf omp
LAMMPS (4 May 2017)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread.
(../comm.cpp:90)
  using 1 OpenMP thread(s) per MPI task
using multi-threaded neighbor list subroutines
Reading data file ...
  orthogonal box = (0 0 0) to (31.28 40.663 29.44)
  2 by 2 by 1 MPI processor grid
  reading atoms ...
  3456 atoms
  scanning bonds ...
  2 = max bonds/atom
  scanning angles ...
  1 = max angles/atom
  reading bonds ...
  2304 bonds
  reading angles ...
  1152 angles
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj: 0 0 0
  special bond factors coul: 0 0 0
  2 = max # of 1-2 neighbors
  1 = max # of 1-3 neighbors
  1 = max # of 1-4 neighbors
  2 = max # of special neighbors
Finding SHAKE clusters ...
  0 = # of size 2 clusters
  0 = # of size 3 clusters
  0 = # of size 4 clusters
  1152 = # of frozen angles
PPPM initialization ...
  extracting TIP4P info from pair style
WARNING: Using 12-bit tables for long-range coulomb (../kspace.cpp:321)
  G vector (1/distance) = 0.258867
  grid = 25 30 24
  stencil order = 5
  estimated absolute RMS force accuracy = 0.00369145
  estimated relative force accuracy = 1.11167e-05
  using double precision FFTs
  3d grid and FFT values/proc = 12958 4680
Last active /omp style is kspace_style pppm/tip4p/omp
Neighbor list info ...
  update every 1 steps, delay 0 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 14.3154
  ghost atom cutoff = 14.3154
  binsize = 7.1577, bins = 5 6 5
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair lj/cut/tip4p/long/omp, perpetual
      attributes: half, newton on, omp
      pair build: half/bin/newton/omp
      stencil: half/bin/3d/newton
      bin: standard
Setting up Verlet run ...
  Unit style : real
  Current step : 0
  Time step : 1
SHAKE stats (type/ave/delta) on step 0
  1 0.854612 0.0716265 6912
  1 109.171 3.05582
Per MPI rank memory allocation (min/avg/max) = 10.68 | 10.88 | 11.08 Mbytes
Step Temp PotEng KinEng TotEng Press
       0 15.002171 -12092.961 102.98699 -11989.974 9961.9129
     100 7.8817572 -18880.145 54.106731 -18826.039 -3929.4866
     200 10.183447 -18895.078 69.907381 -18825.171 -2630.8048
     300 12.622467 -18908.331 86.650779 -18821.68 -3014.7697
     400 12.09286 -18903.123 83.015134 -18820.107 -3450.2393
     500 7.54285 -18873.918 51.780199 -18822.138 -2649.3121
     600 6.8884116 -18866.929 47.287606 -18819.641 -3977.6456
     700 7.7839339 -18875.136 53.435192 -18821.701 -2479.4775
     800 11.398194 -18895.435 78.246384 -18817.189 -3450.9473
     900 14.024945 -18912.645 96.278521 -18816.366 -3011.2215
SHAKE stats (type/ave/delta) on step 1000
  1 0.9572 5.07075e-08 6912
  1 104.52 4.743e-06
    1000 11.626875 -18897.209 79.816235 -18817.392 -2855.3987
Loop time of 27.7428 on 4 procs for 1000 steps with 3456 atoms

Performance: 3.114 ns/day, 7.706 hours/ns, 36.045 timesteps/s
77.9% CPU use with 4 MPI tasks x 1 OpenMP threads