TIP4P-2005 water model

Hello guys,
I’m trying my water model with TIP4P-2005. Here is my script.
However, when I start, it just stays running while nothing happens on the screen no matter how long I wait.
So, is there anything wrong about my script?
Thanks.

#Pure water

units real
dimension 3
boundary p p p
atom_style full

read_data 267watercube.data

#Force field parameters
#H=1, O=2

set type 1 charge 0.5564
set type 2 charge -1.1128

#TIP4P-2005 water
pair_style lj/cut/tip4p/long 2 1 1 1 0.1546 13.0 8.5
pair_modify tail yes mix arithmetic
bond_style harmonic
angle_style harmonic
kspace_style pppm/tip4p 1.0e-5

pair_coeff 1 1 0.0 0.0 #H-H
pair_coeff 2 2 0.1852 3.1589 #O-O
pair_coeff 1 2 0.0 0.0 #O-H
bond_coeff 1 0.0 0.9572 #H-O
angle_coeff 1 0.0 104.52 #H-O-H

neighbor 3.0 bin
neigh_modify every 1 delay 1 check yes

#Group
group water type 1 2
group oxgen type 2

#Calculations

compute tp water temp
compute peratom water stress/atom NULL
compute pr water reduce sum c_peratom[1] c_peratom[2] c_peratom[3]
compute pr3 water reduce sum c_peratom[4] c_peratom[5] c_peratom[6]

variable T equal 300.0
variable P equal 1.0
variable dt equal 0.5
variable ci equal 5 #correlation interval
variable si equal 2 #sample interval
variable ti equal {ci}*{si} #total interval
variable th equal 1000 #thermo interval
variable di equal 1000 #dump interval
variable Tdamp equal {dt}*100 variable Pdamp equal {dt}*1000
variable den equal density

#Run

timestep {dt} thermo {th}
thermo_style custom step c_tp ke pe etotal vol press v_den

minimize 1.0e-4 1.0e-6 100 1000

fix shakeh2o water shake 0.0001 200 0 b 1 a 1

fix l all langevin {T} {T} ${Tdamp} 19930830
fix e all nve

dump 1 all atom ${di} langevin-nve.lammpstrj

run 10000

unfix l
unfix e
undump 1

fix 1 all npt temp {T} {T} {Tdamp} iso {P} {P} {Pdamp}

dump 1 all atom ${di} npt.lammpstrj

run 10000

undump 1

which stage of the script is stalling? your output

should tell you that. If its the minimize, try running

with larger tolerance and or variable th=1

Steve

Hello guys,
    I'm trying my water model with TIP4P-2005. Here is my script.
    However, when I start, it just stays running while nothing happens on
the screen no matter how long I wait.
    So, is there anything wrong about my script?

​you have a force constant of 0 for both bond and angle interaction. that
would be fine for as long as you *also* use fix shake, but that is not the
case during minimization (because fix shake only works for MD). so there is
nothing keeping your hydrogens from getting extremely close to the oxygens
and thus leading to the explosion of potential energy and ultimately the
divergence you see. GI-GO.

axel.​