Issue with TIP4P Flexible Potential for water

Hi,

I am trying to run a simulation with an already equilibrated structure in a box that is enough space for additional atoms, but when I try to extend the box dimensions a get the following error for LAMMPS:

Setting up run …

WARNING: Inconsistent image flags (…/domain.cpp:645)

Memory usage per processor = 8.72688 Mbytes

Step Temp E_pair E_mol TotEng Press

0 100 -2834.819 1099.2163 -1542.7443 -2045.3261

ERROR on proc 10: Bond atoms 484 485 missing on proc 10 at step 290 (…/neigh_bond.cpp:65)

ERROR on proc 9: Bond atoms 394 395 missing on proc 9 at step 290 (…/neigh_bond.cpp:65)

APPLICATION TERMINATED WITH THE EXIT STRING: Terminated (signal 15)

The only change in my box dimensions are shown below

-0.0203835 18.6204 xlo xhi

-0.0203835 18.6204 ylo yhi

-0.0203835 18.6204 zlo zhi

to

-0.0203835 18.6204 xlo xhi

-0.0203835 18.6204 ylo yhi

-0.0203835 37.2408 zlo zhi

I can run a simulation with the first set of dimensions but the second yields the error. How could I run a simulation with water in a box larger than the equilibrated structure? My input file is below…My ultimate goal is to equilibrate H2O with SiC to calculate the equilibrium properties of the interface between the two structures. Thanks

INPUT FILE

#TIP4P/2005f model: M. A. González and J. L. F. Abascal, A flexible

#model for water based on TIP4P/2005, J. Chem. Phys. 135 (2011) 224516.

#variables-------------------------------------------------------------------------

variable Text equal 100.0

variable Pext equal 0.0

variable Nrun equal 100000

variable Nrun1 equal 50000

variable Nrun2 equal 150000

variable ts equal 0.2

variable Tdamp equal ${ts}*100

variable Pdamp equal ${ts}*1000

variable Nf equal ${Nrun}/100

variable Ne equal 10

variable Nr equal {Nf}/{Ne}

variable Ndump equal ${Nrun}/2

variable Nr_rdf equal 0.5*{Nrun}/{Ne}

variable watMoleMass equal 18.0153 # /(g/mol)

variable nAvog equal 6.0221415e23 # Avogadro’s number

variable watMoleculeMass equal ({watMoleMass}/{nAvog}) # /(g/molecule)

variable A3_in_cm3 equal 1e-24 # Angstrom^3 in cm^3

variable nAtoms equal atoms

variable nMolecules equal v_nAtoms/3

units real

dimension 3

boundary p p p

atom_style full

read_data H2O_648N_Coord.txt

#create groups ###---------------------------------------------------------

mass 1 1.00794 # H

mass 2 15.9994 # O

#include forcefield.TIP4P-2005.txt

TIP4P/2005 potential parameters-------------------------

group O type 2

group H type 1

group H2O type 1 2

TIP4P/2005 flexible potenrial parameters---------------------------------

pair_style lj/cut/tip4p/long 2 1 1 1 0.1546 14.0

kspace_style pppm/tip4p 1e-5

pair_coeff 1 1 0.0 0.0

pair_coeff 1 2 1 0.0 1.5795

pair_coeff 2 2 0.1852573718 3.1644

bond_style morse

bond_coeff 1 103.38934 2.287 0.9419

angle_style harmonic

angle_coeff 1 43.95435 107.4

neighbor-----------------------------------------------------------------

neighbor 2.0 bin

neigh_modify delay 0 every 1 check yes

initial condition------------------------------------------------------

velocity all create 100 1234546 dist gaussian mom yes rot yes

#dumps and restart-----------------------------------------------------------

dump waterdump all custom 1000 flex.lammpstrj id type x y z vx vy vz

dump trj all atom ${Nf} WaterTrj.data

dump_modify trj scale no sort id

dump deqlb all atom ${Nrun} FinalCoord.data

dump_modify deqlb scale no sort id

dump coord all custom ${Ndump} DumpCoord.xyz element x y z

#dump_modify coord element H O

#dump_modify deqlb element H O

restart 200000 t.restart1 t.restart2

timestep ${ts}

fix 2 all nvt temp 1 ${Text} 250

run ${Nrun1}

unfix 2

fix 3 all nvt temp {Text} {Text} 250

run ${Nrun2}

unfix 3

fix 4 all npt iso {Pext} {Pext} {Pdamp} temp {Text} {Text} {Tdamp}

variable Dens equal v_nMolecules*{watMoleculeMass}/(vol*{A3_in_cm3})

fix DensAve all ave/time {Ne} {Nr} ${Nf} v_Dens file wat.dens.data

compute bonds all property/local btype batom1 batom2

compute angles all property/local atype aatom1 aatom2 aatom3

compute vel H2O property/atom vx vy vz

#dumps

dump bonds all local ${Ndump} dump.bond.txt c_bonds[1] c_bonds[2] c_bonds[3]

dump angles all local ${Ndump} dump.angle.txt c_angles[1] c_angles[2] c_angles[3] c_angles[4]

output-------------------------------------------------------------------

thermo_style custom step temp pe ke etotal press vol lx

thermo_modify flush yes norm yes

thermo ${Nf}

run ${Nrun}

dump velocities H2O custom ${Nf} Velocities.txt c_vel[2]

write_restart restart.TIP4P_1procs_0.2fs_648N_100K

Hi,

I am trying to run a simulation with an already equilibrated structure in a
box that is enough space for additional atoms, but when I try to extend the
box dimensions a get the following error for LAMMPS:

Setting up run ...

WARNING: Inconsistent image flags (../domain.cpp:645)

Memory usage per processor = 8.72688 Mbytes

Step Temp E_pair E_mol TotEng Press

       0 100 -2834.819 1099.2163 -1542.7443 -2045.3261

ERROR on proc 10: Bond atoms 484 485 missing on proc 10 at step 290
(../neigh_bond.cpp:65)

ERROR on proc 9: Bond atoms 394 395 missing on proc 9 at step 290
(../neigh_bond.cpp:65)

APPLICATION TERMINATED WITH THE EXIT STRING: Terminated (signal 15)

The only change in my box dimensions are shown below

  -0.0203835 18.6204 xlo xhi

  -0.0203835 18.6204 ylo yhi

  -0.0203835 18.6204 zlo zhi

to

  -0.0203835 18.6204 xlo xhi

  -0.0203835 18.6204 ylo yhi

  -0.0203835 37.2408 zlo zhi

I can run a simulation with the first set of dimensions but the second
yields the error. How could I run a simulation with water in a box larger
than the equilibrated structure? My input file is below......My ultimate
goal is to equilibrate H2O with SiC to calculate the equilibrium properties
of the interface between the two structures. Thanks

you have water molecules in your system that straddle the periodic
boundary. atom positions are recorded in LAMMPS as positions inside
the simulation cell and counters for how many times an atoms has
passed through the box in that direction. for a molecule that
straddles the box boundary, one of those counters is larger by one for
some atoms than for the rest. by increasing the box length, you propel
those atoms far apart resulting in the warning and the errors you see.

try using the change_box command instead of editing the data file to
increase the box length.

axel..