Problem regarding TIP4P/SPC water model

Hello everyone,

I built a TIP4P water model with packmol with 357 molecules in 22x22x22 angstrom box. Then I used Topotools to create the bonds and angles among the molecules. When I ran the simulation in LAMMPS it was giving me non-numeric box dimension when I used NPT ensmble. When I used NVT, it was giving me out of range atoms-can’t compute pppm.

After consulting the previous emails regarding these sorts of problems from LAMMPS user mailing list, I understood that my simulation might be unstable. So, I changed the box dimensions to 5x5x5 angstrom with only 3 molecules of water just to see what I was doing wrong. Still I have the same problems. I used the input parameters from LAMMPS documentation.

I tried SPC water model with 3 molecules in 5x5x5 and 20x20x20 box dimensions. The pressure I got is of the order of 10^34 or NAN. I checked with different timesteps. So, I am stuck with both TIP4P and SPC water model. I am sending you both my input file and screenshot of the errors. I will greatly appreciate someone’s help on this.

TIP4P (3 molecules):

tip4p water

atom_style full
units real
dimension 3
boundary p p p
read_data data.h20_small

group hydrogen type 1
group oxygen type 2
#group all type 1 type 2

set group oxygen charge -1.11
set group hydrogen charge 0.5564

velocity all create 300.0 123456789 dist gaussian mom yes rot yes

neighbor 2.0 bin

neigh_modify delay 0

Define interaction parameters

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

pair_coeff 1 1 0.0 0.0
pair_coeff 1 2 0.0 0.0
pair_coeff 2 2 0.16275 3.16435

bond_style harmonic
bond_coeff 1 0.0 0.9572

angle_style harmonic
angle_coeff 1 0.0 104.52

fix 1 all shake 0.0001 100 1 b 1 a 1

initialize

thermo 10
timestep 0.001

------------- Equilibration and thermalisation ----------------

fix 2 all press/berendsen iso 1.0 1.0 0.1

fix 3 all nvt temp 300.0 300.0 0.1
#fix 3 all npt temp 300.0 300.0 0.1 iso 1.01325 1.01325 1000

dump 2 all xyz 1 water.xyz
run 500

DUMP file:
9
Atoms. Timestep: 0
2 -0.000576 1.52239 3.71945
1 0.758776 0.941048 3.42646
1 0.357216 2.31842 4.20638
2 1.5924 3.46099 3.95979
1 1.82369 2.91375 4.76444
1 0.71911 3.15121 3.58537
2 0.754134 1.1164 4.88936
1 0.456617 0.234294 4.52359
1 1.21872 1.63526 4.17261
9
Atoms. Timestep: 1
2 -377.653 -275.741 -273.393
1 -20.29 -526.907 -398.175
1 -196.759 79.1415 -54.9923
2 379.068 644.728 -159.406
1 467.321 386.103 203.482
1 -18.7762 490.423 -315.206
2 -19.1083 -468.184 282.362
1 -146.438 -863.734 98.8119
1 188.883 -247.404 -55.8088
9
Atoms. Timestep: 2
2 -nan -nan -nan
1 -nan -nan -nan
1 -nan -nan -nan
2 -nan -nan -nan
1 -nan -nan -nan
1 -nan -nan -nan
2 -nan -nan -nan
1 -nan -nan -nan
1 -nan -nan -nan

SPC file 3 molecules (22x22x22 angstrom):

SPC water

atom_style full
units real
dimension 3
boundary p p p
read_data data.lowdensity

group hydrogen type 1
group oxygen type 2
#group all type 1 type 2

set group oxygen charge -0.82
set group hydrogen charge 0.41

#velocity all create 200.0 123456789 dist gaussian mom yes rot yes
velocity all create 298.0 1234

neighbor 2.0 bin

neigh_modify every 1 delay 0 check yes
timestep 2.0

Define interaction parameters

pair_style lj/cut/coul/long 9.8
pair_modify tail yes
kspace_style pppm 1.0e-5

pair_coeff 1 1 0.0 0.0
pair_coeff 1 2 0.0 0.0
pair_coeff 2 2 0.1553 3.166

bond_style harmonic
bond_coeff 1 0.0 1.0

angle_style harmonic
angle_coeff 1 0.0 109.47

fix 1 all shake 0.0001 10 1 b 1 a 1
fix 2 all npt temp 298.0 298.0 1.0 iso 1.01325 1.01325 100
fix 3 all momentum 1 linear 1 1 1

dump 2 all xyz 1 spc.xyz
thermo 100
run 1000

#unfix 2

#reset_timestep 0
#fix 2 all nvt temp 300.0 300.0 0.1
#thermo 10
#run 500

initialize

#thermo_style custom step temp pe ke etotal press vol density cpu

------------- Equilibration and thermalisation ----------------

#fix 2 all press/berendsen iso 1.0 1.0 0.1

Regards,
Baig

Hello everyone,

I built a TIP4P water model with packmol with 357 molecules in 22x22x22
angstrom box. Then I used Topotools to create the bonds and angles among the
molecules. When I ran the simulation in LAMMPS it was giving me non-numeric
box dimension when I used NPT ensmble. When I used NVT, it was giving me out
of range atoms-can't compute pppm.
After consulting the previous emails regarding these sorts of problems from
LAMMPS user mailing list, I understood that my simulation might be unstable.
So, I changed the box dimensions to 5x5x5 angstrom with only 3 molecules of
water just to see what I was doing wrong. Still I have the same problems. I
used the input parameters from LAMMPS documentation.

please note, that you cannot simply copy parameters without knowing
what units they are in.
in the cases below, your time constants for the thermostat/barostat
fixes seem awfully small for your choices of units and timestep.

in general, it is a bad idea to run a variable cell simulation on a
non-equilibrated system with high potential energy.

in addition, if you have difficulties to get a stable simulation from
a generated initial configuration, you should try and do the simplest
method first, i.e. using fix nve for time integration and then perhaps
use a dissipative thermostat (langevin, temp/csld, etc.) to better
handle transfer of energy.

axel.

Hello everyone,

Thank you Axel for replying to me. I changed the input file for SPC model to a low density water box. Now, I am having numeric index out of bound errors. I defined types of atoms, bonds, angles, impropers and dihedrals in the data file, and the corresponding coefficients are in the input file. I don’t understand where my atom IDs might be wrong because there’s only 9 atoms in the box. Can someone please help me? I have copied both input file and data file for you to see.

SPC water box (3 molecules)

SPC water

echo screen
log debug.log

units real

neighbor 2.0 bin
#neigh_modify delay 0 every 1 one 10 check yes
neigh_modify delay 5

boundary p p p

atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style none
improper_style none

pair_style lj/cut/coul/cut 9.0 9.0

read_data data.lowdensity
pair_modify mix arithmetic tail yes

group hydrogen type 1
group oxygen type 2
group water type 1 2

set group oxygen charge -0.82
set group hydrogen charge 0.41

pair_coeff 1 1 0.0 0.0
pair_coeff 1 2 0.0 0.0
pair_coeff 2 2 0.1553 3.166

bond_coeff 2 0.0 1.0
angle_coeff 2 0.0 109.47

velocity water create 300.0 12345 dist gaussian

fix 1 water shake 1.0e-8 100 0 b 1 a 1
fix 2 water npt temp 300.0 300.0 100.0 iso 1.01325 1.01325 500

timestep 1.0
thermo 50
thermo_style custom step temp press vol etotal ke pe evdwl ecoul elong

dump 1 all xyz 1 spc.xyz
run 10000

DATA file:
LAMMPS data file. CGCMM style. atom_style full generated by VMD/TopoTools v1.5 on Tue Dec 01 16:55:25 CST 2015
9 atoms
6 bonds
6 angles
0 dihedrals
0 impropers
2 atom types
1 bond types
1 angle types
0 dihedral types
0 improper types

-10.00000 30.00000 xlo xhi
-10.00000 30.00000 ylo yhi
-10.00000 30.00000 zlo zhi

# Pair Coeffs

Hello everyone,

Thank you Axel for replying to me. I changed the input file for SPC model to
a low density water box. Now, I am having numeric index out of bound errors.
I defined types of atoms, bonds, angles, impropers and dihedrals in the data
file, and the corresponding coefficients are in the input file. I don't
understand where my atom IDs might be wrong because there's only 9 atoms in
the box. Can someone please help me? I have copied both input file and data
file for you to see.

look closely at the output to the screen and note the line in the
input where the error occurs.
it has nothing to do with atom indices, but bond (and angle) type.
your data file has only one type (1), but you are setting parameters
for type 2. that clearly means the type index is out of bounds.

axel.

SPC water box (3 molecules)
# SPC water

echo screen
log debug.log

units real

neighbor 2.0 bin
#neigh_modify delay 0 every 1 one 10 check yes
neigh_modify delay 5

boundary p p p

atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style none
improper_style none

pair_style lj/cut/coul/cut 9.0 9.0

read_data data.lowdensity
pair_modify mix arithmetic tail yes

group hydrogen type 1
group oxygen type 2
group water type 1 2

set group oxygen charge -0.82
set group hydrogen charge 0.41

pair_coeff 1 1 0.0 0.0
pair_coeff 1 2 0.0 0.0
pair_coeff 2 2 0.1553 3.166

bond_coeff 2 0.0 1.0
angle_coeff 2 0.0 109.47

velocity water create 300.0 12345 dist gaussian

fix 1 water shake 1.0e-8 100 0 b 1 a 1
fix 2 water npt temp 300.0 300.0 100.0 iso 1.01325 1.01325 500

timestep 1.0
thermo 50
thermo_style custom step temp press vol etotal ke pe evdwl ecoul elong

dump 1 all xyz 1 spc.xyz
run 10000

DATA file:
LAMMPS data file. CGCMM style. atom_style full generated by VMD/TopoTools
v1.5 on Tue Dec 01 16:55:25 CST 2015
9 atoms
6 bonds
6 angles
0 dihedrals
0 impropers
2 atom types
1 bond types
1 angle types
0 dihedral types
0 improper types

[...]