Bonds Not Crossing the Boundary/Image Plane

Dear Sir/Madam,

I’ve been trying to run a 2D simulation on a simple hexagonal system (resembling a graphene sheet).
The potential that we are using is very simple and consists of only harmonic bonds and angles; there are no non-bonding interactions or dihedrals etc…

I’ve generated the structure on material studio v6.0 and have converted the structure into a readable format via msi2lmp and have made some slight modifications to the file to suit my needs.
I’m running a simple minimization simulation however the program is returning a number of warnings which are:

  • WARNING: Inconsistent image flags (…/domain.cpp:594)
  • WARNING: Bond/angle/dihedral extent > half of periodic box length (…/domain.cpp:667)

I’ve checked the system using VMD and it appears that the atoms at the boundary are not bonding with their image atoms but with the atoms in the same cell (as can be seen in this image: http://tinypic.com/view.php?pic=2j29euc&s=5#.UrgQ0vRDtl0).

I’ve defined all of the bonds in the structure file (even bonds crossing the image plane) however the problem still persist. From the diagnostic tests that I’ve carried out it appears that the bonds crossing the periodic boundaries are being bonded with atoms within the same cell. I’ve tried searching for a possible solution and have tried tampering about with the script as well as structure file however I have not managed to address this problem.

For reference, here’s a copy of the script that I’m using: http://pastebin.com/raw.php?i=9LjwJmSR

Kind Regards & Best Wishes,
John Smith

Dear Sir/Madam,

I've been trying to run a 2D simulation on a simple hexagonal system
(resembling a graphene sheet).
The potential that we are using is very simple and consists of only harmonic
bonds and angles; there are no non-bonding interactions or dihedrals etc...

I've generated the structure on material studio v6.0 and have converted the
structure into a readable format via msi2lmp and have made some slight
modifications to the file to suit my needs.
I'm running a simple minimization simulation however the program is
returning a number of warnings which are:

- WARNING: Inconsistent image flags (../domain.cpp:594)
- WARNING: Bond/angle/dihedral extent > half of periodic box length
(../domain.cpp:667)

I've checked the system using VMD and it appears that the atoms at the
boundary are not bonding with their image atoms but with the atoms in the
same cell (as can be seen in this image:
http://tinypic.com/view.php?pic=2j29euc&s=5#.UrgQ0vRDtl0).

i don't see any bonds at all in this image.

please note that a) vmd's periodic display is just taking the regular
rendered cell and replicating it, so it cannot directly render bonds
between periodic images and b) that LAMMPS will set up bonds to the
"closest image", thus the warning about a bond longer than half the
periodic box is harmless in your case. the warning is printed because
LAMMPS cannot distinguish between your case and the case where a
molecule is "broken" due to inconsistent image flags (which can spell
problems down the road, e.g. when trying to use the replicate command.

I've defined all of the bonds in the structure file (even bonds crossing the
image plane) however the problem still persist. From the diagnostic tests
that I've carried out it appears that the bonds crossing the periodic
boundaries are being bonded with atoms within the same cell. I've tried
searching for a possible solution and have tried tampering about with the
script as well as structure file however I have not managed to address this
problem.

For reference, here's a copy of the script that I'm using:
http://pastebin.com/raw.php?i=9LjwJmSR

this is useless without the data file.

axel.

Dear Dr. Kohlmeyer,

First of all I’d like to thank you for your quick reply.

We’ve managed to solve the warning about the bond longer than half the periodic box length by enlarging the unit cell size used.

The only issue I have is that, whenever I’m giving very strong force constants to the bonds spanning the periodic boundaries (and weaker force constants to the rest), the atoms next to the boundaries move inwards. Thus we are suspecting that the atoms are still bonding with the atoms in the same cell however this should give a warning but does not if the size is big enough.

Please find the data file in the following url: http://pastebin.com/raw.php?i=f4tvLgEK

Regards,
John Smith

Dear Sir/Madam,

I’ve been trying to run a 2D simulation on a simple hexagonal system
(resembling a graphene sheet).
The potential that we are using is very simple and consists of only harmonic
bonds and angles; there are no non-bonding interactions or dihedrals etc…

I’ve generated the structure on material studio v6.0 and have converted the
structure into a readable format via msi2lmp and have made some slight
modifications to the file to suit my needs.
I’m running a simple minimization simulation however the program is
returning a number of warnings which are:

  • WARNING: Inconsistent image flags (…/domain.cpp:594)
  • WARNING: Bond/angle/dihedral extent > half of periodic box length
    (…/domain.cpp:667)

I’ve checked the system using VMD and it appears that the atoms at the
boundary are not bonding with their image atoms but with the atoms in the
same cell (as can be seen in this image:
http://tinypic.com/view.php?pic=2j29euc&s=5#.UrgQ0vRDtl0).

i don’t see any bonds at all in this image.

please note that a) vmd’s periodic display is just taking the regular
rendered cell and replicating it, so it cannot directly render bonds
between periodic images and b) that LAMMPS will set up bonds to the
“closest image”, thus the warning about a bond longer than half the
periodic box is harmless in your case. the warning is printed because
LAMMPS cannot distinguish between your case and the case where a
molecule is “broken” due to inconsistent image flags (which can spell
problems down the road, e.g. when trying to use the replicate command.

I’ve defined all of the bonds in the structure file (even bonds crossing the
image plane) however the problem still persist. From the diagnostic tests
that I’ve carried out it appears that the bonds crossing the periodic
boundaries are being bonded with atoms within the same cell. I’ve tried
searching for a possible solution and have tried tampering about with the
script as well as structure file however I have not managed to address this
problem.

For reference, here’s a copy of the script that I’m using:
http://pastebin.com/raw.php?i=9LjwJmSR

this is useless without the data file.

axel.

Dear Dr. Kohlmeyer,

First of all I'd like to thank you for your quick reply.

We've managed to solve the warning about the bond longer than half the
periodic box length by enlarging the unit cell size used.

The only issue I have is that, whenever I'm giving very strong force
constants to the bonds spanning the periodic boundaries (and weaker force
constants to the rest), the atoms next to the boundaries move inwards. Thus
we are suspecting that the atoms are still bonding with the atoms in the
same cell however this should give a warning but does not if the size is big
enough.

you are not making much sense here. the cell dimensions have to be set
properly to so that all bonds are of the same length. you cannot just
make the box arbitrarily larger and expect the system to behave is if
it is not.

axel.

Dear Dr. Kohlmeyer,

Firstly, thank you for your rapid reply.

The cell dimensions and atomic positions are actually calculated so that the system is already in its ideal minimum. As can be seen from the structure file and potentials.

The problem is that the atoms forming the row having a minimum or maximum y-value for their atomic positions, move inwards. This movement is not expected, since the system should be in the optimum position already. Thus the problem is that I cannot explain this movement of atoms after a minimization.

Please find the image that I tried to send earlier: http://tinypic.com/r/2rw1d3b/5 (apologies for sending the incorrect URL) depicting the atom positions (without bonds) of the original (left) and minimized (right) structure.

Regards,

Dear Dr. Kohlmeyer,

First of all I’d like to thank you for your quick reply.

We’ve managed to solve the warning about the bond longer than half the
periodic box length by enlarging the unit cell size used.

The only issue I have is that, whenever I’m giving very strong force
constants to the bonds spanning the periodic boundaries (and weaker force
constants to the rest), the atoms next to the boundaries move inwards. Thus
we are suspecting that the atoms are still bonding with the atoms in the
same cell however this should give a warning but does not if the size is big
enough.

you are not making much sense here. the cell dimensions have to be set
properly to so that all bonds are of the same length. you cannot just
make the box arbitrarily larger and expect the system to behave is if
it is not.

axel.

Dear Dr. Kohlmeyer,

Firstly, thank you for your rapid reply.

The cell dimensions and atomic positions are actually calculated so that the
system is already in its ideal minimum. As can be seen from the structure
file and potentials.

The problem is that the atoms forming the row having a minimum or maximum
y-value for their atomic positions, move inwards. This movement is not
expected, since the system should be in the optimum position already. Thus
the problem is that I cannot explain this movement of atoms after a
minimization.

ok. i think i understand now what is going on. a dead giveaway is also
the insanely high energy. the problem is that you don't have the
periodic atom that is closest to the other atom forming the bond is
among the ghost atoms. the amount of ghost atoms is normally inferred
from the non-bonded cutoff. since you don't have that set, you have
to manually tell LAMMPS how far it has to communicate ghost atoms.

if you add:

communicate single cutoff 11.0

it should magically work (well, it does for me).

axel.

From the communicate doc page:

The cutoff option allows you to set a ghost cutoff distance, which is the distance from the borders of a processor’s sub-domain at which ghost atoms are acquired from other processors. By default the ghost cutoff = neighbor cutoff = pairwise force cutoff + neighbor skin. See the neighbor command for more information about the skin distance. If the specified Rcut is greater than the neighbor cutoff, then extra ghost atoms will be acquired. If it is smaller, the ghost cutoff is set to the neighbor cutoff.

These are simulation scenarios in which it may be useful or even necessary to set a ghost cutoff > neighbor cutoff:

  • a single polymer chain with bond interactions, but no pairwise interactions
  • bonded interactions (e.g. dihedrals) extend further than the pairwise cutoff
  • ghost atoms beyond the pairwise cutoff are needed for some computation

In the first scenario, a pairwise potential is not defined. Thus the pairwise neighbor cutoff will be 0.0. But ghost atoms are still needed for computing bond, angle, etc interactions between atoms on different processors, or when the interaction straddles a periodic boundary.

Steve