Non-numeric atom coords in a simultion of TIP4P water

LAMMPS
1

Dear all,

I want to do a simulation of TIP4P water (W.L Jorgensen’s water model). I went through the LAMMPS page for the TIP4P and I learnt that there are two ways to do this. I chose the explicit method.

I tried to set up the system and I am getting Non-numeric atom coords error. Sheer observation of the input structure did not readily tell me about overlapping atoms. Since the bonds and angles are rigid, I suspect this is due to the electrostatic interaction between the dummy particles. But I don’t seem to understand why and how to fix it.

Any help would be appriciated.

My input script:

# -- Initialise system --
units real
atom_style full

# -- read data --
read_data solvent_tip4p.data
displace_atoms all random 0.5 0.5 0.5 21344

# -- set charges --
# water
set  type  1  charge  0.0
set  type  2  charge  0.52
set  type  3  charge  -1.04

# -- non bonded interactions --
pair_style lj/cut/coul/long 10.0 10.0
pair_coeff 1 1 0.155 3.15365
pair_coeff 2 2 0.0 1.0
pair_coeff 3 3 0.0 1.0
pair_modify mix geometric

# -- kspace --
kspace_style pppm 0.00001

# -- bonded interaction
bond_style zero
bond_coeff 1 0.9572

angle_style zero
angle_coeff 1 104.52

# -- write ff data --
write_data ff.data pair ij

# -- water shake --
fix mymilkshake all shake 0.0001 10 10000 b 1 a 1

# -- neighbour list --
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes

# -- simulation protocol --
# EM
minimize 0.0 0.0 1000 10000
write_data postem.data pair ij

# run a nvt equilibriation
timestep 0.5
velocity all create 300.0 21423623
fix thermo0 all nvt temp 300 300 100.0 
thermo_style custom step temp pe ke press density
thermo 10
dump myDump0 all custom 1 nvt_solvent.lammpstrj id mol type xu yu zu
run 5000
unfix thermo0

# equilibriate pressure
timestep 2
velocity all create 300.0 21423623
fix thermo1 all npt temp 300 300 100.0 iso 0.986923 0.986923 2000.0
thermo_style custom step temp pe ke press density
thermo 10
dump myDump1 all custom 100 npt_solvent.lammpstrj id mol type x y z
run 50000
unfix thermo1

# prod
timestep 2
fix thermo2 all npt temp 300 300 100.0 iso 0.986923 0.986923 2000.0
thermo_style custom step temp pe ke emol enthalpy press vol density 
thermo 10
log solvent.log
dump myDump2 all custom 100 prod_solvent.lammpstrj id mol type x y z
run 200000
write_data solvent_final.data nocoeff

I am attatching the input structures that I have used to test this.
solvent_tip4p.data (148.3 KB)

Thank you.

Regards,
Ved