QUESTION about lammps simualtion: WARNING: Bond/angle/dihedral extent > half of periodic box length

Hi, lammps-users

i’m tying to use lammps to test a zeolite potential model for different types of zeolites based on molecular dynamics calculation, and the initial zeolite structure were dumped with another package GULP and has been optimized to have minimum energy, but i got problem with bond length identification in lammps, here is the warnings

WARNING: Inconsistent image flags (…/domain.cpp:774)
WARNING: Bond/angle/dihedral extent > half of periodic box length (…/domain.cpp:895)
in the dump file from GULP (i have attached), all the bonds in the configuration is defined, some of the bonds connect atoms on the opposite sides of the unit cell, but it make sense, under periodic boundary condition. The problem is that it seems lammps directly evaluates those bonds length across the unit cell, like the periodic boundary condition is not implemented, which gives huge amount of positive total energy of the initial configuration. However, the same configuration gives negative energy in GULP, which is much more reasonable, in the case of zeolite. Further checking on the visualization of minimisation in lammps shows that the unit cell shrink dramatically. That is why i am convinced that the issue here is the bond length identificatio_n._
the magic thing is that for sodalite, quartz and other zeolite types, the warning associated with bond length won’t appear when increase my simulation box to be 3 by 3 by 3, but for faujasite, which is much more complicated, simply increase the box size can work any more. So i am seeking for another well- grounded solution to this pr_oblem. I have red the lammps mailing list, people are asking the same question, but it does no help. The input script is post as follows**._ Any piece of help will be appreaciated_.**_

lammps input script for Faujasite with MZHB

units metal
dimension 3
boundary p p p
atom_style full

lattice custom 34.500900 a1 0 1 1 a2 1 0 1 a3 1 1 0

read_data MZHB_SOD_original.lmp

potential

bond_style harmonic
bond_coeff 1 11.650000 1.6200000

angle_style harmonic
angle_coeff 1 3.4000000 109.47000
angle_coeff 2 1.1100000 149.80000

pair_style lj/cut/coul/long 12.000
pair_coeff 1 1 0.0085997000 1.9599828638 12.0000
pair_coeff 1 2 0.0052661900 1.7580342186 12.0000
pair_coeff 2 2 0.0032249000 1.5768920229 12.0000

special_bonds lj 1.0 1.0 1.0 coul 1.0 1.0 1.0

kspace_style ewald 1.0e-4

MZHB_FAU_original.lmp (255 KB)

LAMMPS does not care if bonds cross periodic boundaries,
it should compute the energy/forces in the bond
correctly regardless. The warning you are getting
has to do with the image flags for the 2 atoms in
the bond that crosses the boundary. They are probably
the same, when they should be different by +/- 1.
But that is just a bookkeeping issue that does not
affect dynamics, but things like computing the
center-of-mass of a molecule (as a diagnostic).

So something else must be causing the huge initial
energy. Possibly your periodic box size is not
commensurate with the lattice of atoms, i.e. 2
atoms on opposite sides of the box are highly
overlapped in a periodic sense.

Steve

Hi Steve,

Thanks for you reply!! But i have check my initial structure for lammps, and it looks normal, actually the same as the structure optimized by GULP, since it won’t affect the dynamics, i will leave it for now, However, i got another problem, I am running a NVE simulation for a shell model, in order to calculate the dynamic properties like IR spectrum from the trajectory. Since in lammps, the core and shell of the atoms are specified as different atoms, i actually got an xyz file with the coordinate of cores and shells respectively, i am looking for some help as how i can use fixed partial charges method to do the IR calculation based on the xyz file. Any piece of help would be greatly apreciated. Thanks

Best regards

Jiasen Guo

Hi Steve,

Thanks for you reply!! But i have check my initial structure for lammps, and
it looks normal, actually the same as the structure optimized by GULP, since
it won't affect the dynamics, i will leave it for now, However, i got

you haven't looked at your structure correctly. you've neglected too
look at the dimensions and shape of your simulation cell.
you data file has a structure for a triclinic cell, but your box
definition is for a cubic box that is the min/max of the extend of
your triclinic coordinates.
that *is* incorrect, because of that your bonds that would be "normal"
and acceptable for periodic boundaries, are no longer acceptable, and
stretched very long.

so you must not ignore this and correct your data file.

axel.

Hi Steve,

Thanks for you reply!! But i have check my initial structure for lammps, and
it looks normal, actually the same as the structure optimized by GULP, since
it won't affect the dynamics, i will leave it for now, However, i got

you haven't looked at your structure correctly. you've neglected too
look at the dimensions and shape of your simulation cell.
you data file has a structure for a triclinic cell, but your box
definition is for a cubic box that is the min/max of the extend of
your triclinic coordinates.
that *is* incorrect, because of that your bonds that would be "normal"
and acceptable for periodic boundaries, are no longer acceptable, and
stretched very long.

so you must not ignore this and correct your data file.

as an add-on remark. part of your problem is, that your geometry of
your unit cell does not follow the requirements that LAMMPS imposes.
you will have to reorient your crystal so that the a vector is on the
x-axis, the b vector in the xy-plane and the c vector so you get a
right handed coordinate system.
you can see the difference, when visualizing your geometry with
showing periodic replica.

axel.