Rounding errors with special EDGE key word

I have a suggestion regarding the EDGE specifier. The following
demonstrates a problem, at least it does on my machine:

units metal
lattice bcc 2 orient x -1 -1 2 orient y 1 -1 0 orient z 1 1 1
region simbox block 0 6 0 7 0 17
create_box 1 simbox
region mbox block EDGE EDGE EDGE EDGE 0 17
region mbox2 block EDGE EDGE 0 7 0 17
region mbox3 block -0.00000001 6.00000001 EDGE EDGE -0.00000001 17
mass 1 1
pair_style lj/cut 2.5
pair_coeff * * 1.0 1.0

# Work correctly
#create_atoms 1 box
#create_atoms 1 region mbox3
# Work incorrectly
#create_atoms 1 region simbox
#create_atoms 1 region mbox
#create_atoms 1 region mbox2

Uncommenting any of the "create_atoms" commands will demonstrate what
happens. All of these commands should generate 5712 atoms, but due to
rounding errors some of them generate fewer. The command,
   create_atoms 1 region simbox
generates 5693, with some atoms missing in the x and z directions. This
happens because there are (very tiny) rounding errors at the edges, which
means we get things like -1.0E-14 < 0.00000, and it declares it outside
the box and an atom is then missing. This happens most often with
rotated lattices, of course.

Does it make sense to add/subtract a very small number (the machine
epsilon, for example) to any boundary defined by the EDGE command, to make
sure the boundary of the region lies outside the box boundary? Otherwise,
this causes problems such as the above that EDGE was designed to solve in
the first place, I suspect.

Karl D. Hammond
University of Tennessee, Knoxville
[email protected]...

"You can never know everything, and part of what you know is always
   wrong. A portion of wisdom lies in knowing that. A portion of courage
   lies in going on anyway."
"Nothing ever goes as you expect. Expect nothing, and you will not be
   surprised."

karl,

Does it make sense to add/subtract a very small number (the machine
epsilon, for example) to any boundary defined by the EDGE command, to make
sure the boundary of the region lies outside the box boundary? Otherwise,
this causes problems such as the above that EDGE was designed to solve in
the first place, I suspect.

in my opinion, the root of the problem that you are observing
is that the atom positions of the lattice are ambiguous, since
they are defined to be *exactly* on top of the boundary.

just shifting your box by some small amount (or the
corresponding region definitions) in each direction will
remove this ambiguity and you get a consistent number
of atoms.

ciao,
   axel.

The problem is more fundamental, and it has
to do with create_atoms, not with region definitions.

From the create_atoms doc page:

For the {box} style, the create_atoms command fills the entire
simulation box with atoms on the lattice. If your simulation box is
periodic, you should insure its size is a multiple of the lattice
spacings, to avoid unwanted atom overlaps at the box boundaries. If
your box is periodic and a multiple of the lattice spacing in a
particular dimension, LAMMPS is careful to put exactly one atom at the
boundary (on either side of the box), not zero or two.

For the {region} style, the geometric volume is filled that is inside
the simulation box and is also consistent with the region volume. See
the "region"_region.html command for details. Note that a region can
be specified so that its "volume" is either inside or outside a
geometric boundary. Also note that if your region is the same size as
a periodic simulation box (in some dimension), LAMMPS does not
implement the same logic as with the {box} style, to insure exactly
one atom at the boundary. if this is what you desire, you should
either use the {box} style, or tweak the region size to get precisely
the atoms you want.

If you are using a region with create_atoms (instead of box),
then as Axel said, if you have atoms exactly (within epsilon)
on the boundaries, it is up to you to detect and fix that by adjusting
your region extent or your box.

Steve