Problems with TIP4Pice

Hi everyone,

I want to do NVT and NPT simulations of liquid water using the TIP4P/ice water model. But the simulations crash unless I use a very small (<0.05fs) timestep. I use all the tricks in the book i.e. energy minimization before the actual simulation etc, but irrespective, the simulations crash after a while unless I choose a very small [unpractical] timestep stating that one hydrogen is missing. Below I include my input script and data file for one such simulation.

PS, I am using the most recent release of LAMMPS released on July 4, 2012. Also the problem persists even if I don’t use ‘replicate’ and use energy minimization on a configuration that I have obtained otherwise.

Regards

Amir

DATA FILE single1.dat

One tip4p/ice water molecule

0 3.9199999 xlo xhi
0 3.9199999 ylo yhi
0 3.9199999 zlo zhi
2 atom types
3 atoms
1 bond types
2 bonds
1 angle types
1 angles

Masses

1 15.9994
2 1.0079

Atoms

1 1 1 -1.1794 2.0000 2.0000 2.0000
2 1 2 0.5897 1.243049672736339 1.414117723381705 2.0000
3 1 2 0.5897 2.756950327263661 1.414117723381705 2.0000

Velocities

1 0 0 0
2 0 0 0
3 0 0 0

Bond Coeffs

1 100 0.9572

Bonds

1 1 1 2
2 1 1 3

Angle Coeffs

1 300 104.52

Angles

1 1 2 1 3

==== INPUT SCRIPT ====

Simulation of tip4p/ice water NPT seed 1

units real
dimension 3
boundary p p p
atom_style full
atom_modify map array

bond_style harmonic
angle_style harmonic
pair_style lj/cut/coul/long/tip4p 1 2 1 1 0.1577 8.5 8.5

read_data single1.dat
replicate 16 16 16

group oxy type 1 # Name atom types
group hyd type 2
group water union oxy hyd

pair_coeff 1 1 0.210701615058703 3.1668
pair_coeff 2 2 0.0 0.0

kspace_style pppm/tip4p 1e-4 # Electrostatics
dielectric 1.0
pair_modify tail yes # long-range correction for truncated Lennard-Jones interactions

fix 5 all shake 1e-12 200000 0 b 1 a 1 # Make molecules rigid using shake

thermo_style custom step temp press vol ke pe emol etotal

neighbor 1.0 bin
neigh_modify delay 0 every 1 check yes

thermo 100
timestep 1.0 #nothing greater than 0.05 works
dump 1 all xyz 100 position.xyz
fix 6 all npt temp 220.0 220.0 200.0 iso 1.0 1.0 1000.0

restart 100 tip4piceNPT1a.out tip4piceNPT1b.out
run 1000000

But the simulations crash

How does it crash? Do you get an error message? Does it run
for a while? Does the thermodynamic output go bad? Running
NPT with a single molecule of water is an odd model. If you are not
careful, the box will collapse around it rapidly and the molecule will
interact with itself in an odd manner.

Steve