Problems with bicrystal lattice configuration

Dear LAMMPS users,

I've been studying grain boundaries' lowest energy configuration through LAMMPS,
but I've found few errors concerning on it. I'm not sure whether these
problems comes from my script error, or inherent bugs in lattice
command in LAMMPS.

<scripts>

variable Cu equal 3.615
lattice fcc ${Cu}

variable Xsize equal (16^2+1)^(0.5)*5
variable Zsize equal 3
variable Ysize equal 20
variable Ymsize equal -${Ysize}

region Rbox prism 0 \{Xsize\} {Ymsize} \{Ysize\} 0 {Zsize} 0.0 0.0 0.0
region Rugrain prism 0 \{Xsize\} 0 {Ysize} 0 \{Zsize\} 0\.0 0\.0 0\.0 region Rdgrain prism 0 {Xsize} \{Ymsize\} 0 0 {Zsize} 0.0 0.0 0.0

create_box 2 Rbox

lattice fcc 3.615 orient x 16 1 0 orient y -1 16 0 orient z 0 0 1
create_atoms 1 region Rugrain
lattice fcc 3.615 orient x 16 -1 0 orient y 1 16 0 orient z 0 0 1
create_atoms 2 region Rdgrain

group Gugrain region Rugrain
group Gdgrain region Rdgrain

region Rinterface block INF INF -10 10 INF INF units box
group Ginterface region Rinterface

problems are:

1. Unexpected voids are found in the boundary of a cell. To be detail,
this voids are not observable Inside the cell, but at the very
boundary between the cell and its translated image by x direction.
After few days of struggling, I've found that this problem is
sometimes solved through changing number of processors (from 16 to 1),
but on other times, it is solved by changing the x-direction size of
the cell.
(I'm running multiple simulations by changing interface misorientation
angle. in the above script, it is 7.15 deg. So, the other time means,
literally, another misorientation angle)
I cannot grasp any consistency to solve this problem.

2. While the number of atoms in the interface region has to be
constant w/ changing y-direction size of the cell, (as shown in the
script) i've found it does not.
for instance, when Ysize = 21, number of atoms in Ginterface = 1053,
even though when Ysize = 20, number of atoms in Ginterface = 1044.
it is very crucial issue in grain boundary energy because slight
change can change overall configuration of a interface.

it anyone have idea on these problems, please leave any comments.

thanks in advance.

Hayoung Chung.

The way the lattice command works is simple
and documented in great detail on the lattice doc
page. If you are using it to fill 2 separate regions
with 2 different orientations of atoms, it is
simply doing the simple thing twice. As the doc
page states, when filling a region, you may or may
not get atoms exactly on the boundaries of the region
due to round-off issues. So it is up to you to check
that and compensate for it if necessary. The best
way to check the lattice and create_atoms command
are producing what you want is to dump the result
to a file and check it manually or by visualization.
You can also use the delete_atoms command after
creating the GB to remove (nearly) overlapping atoms.

I can't think of any case where changing the number
of processors will change what the create_atoms
command produces. It may change the numbering/ordering
of the atoms, but not the number created I don't think.

You are also free to build a grain boundary geometry
exactly how you wish using an external program
and simply input the atom coords into LAMMPS, e.g.
via a read_data file.

Steve

The way the lattice command works is simple
and documented in great detail on the lattice doc
page. If you are using it to fill 2 separate regions
with 2 different orientations of atoms, it is
simply doing the simple thing twice. As the doc
page states, when filling a region, you may or may
not get atoms exactly on the boundaries of the region
due to round-off issues. So it is up to you to check
that and compensate for it if necessary. The best
way to check the lattice and create_atoms command
are producing what you want is to dump the result
to a file and check it manually or by visualization.
You can also use the delete_atoms command after
creating the GB to remove (nearly) overlapping atoms.

I can't think of any case where changing the number
of processors will change what the create_atoms
command produces. It may change the numbering/ordering
of the atoms, but not the number created I don't think.

yes, steve, it can, but it requires two condidtions:
you have to have lattice positions _exactly_ on
the region/box boundaries, and you have to use
an x86 process with the floating point unit in 80-bit
mode (the default). in this case it matters, if numbers
are computed using variables that can remain in
registers or not, since having to store a number from
a register to the stack enforces rounding to 64-bit
and then suddenly what is "in" or "out" can change.
this particularly happens when using higher
optimization levels, which are allowed to violate
some IEEE-754 rounding rules.

this effect can be limited by requiring the FP unit
to be in 64-bit mode using -pc64 (intel compilers)
or -mpc64 (recent GNU compilers).

cheers,
    axel.

Dear Axel and Steve.

Thanks for such quick and candid response.
However, I'm having a hard time understanding Axel's comment. Does it
means that I can circumvent the existing problem by increasing a cell
slightly? like, from "3 lattice" to "3.01 lattice". Or, does executing
simulation by single-core processor is the only solution if I do not
have any other choices in CPU structure?

And again, (sorry) does Steve's comment on "lattice command" applies
to my second question too? I started looking into the source code on
lattice, but i still cannot understand why altering y direction
influences the number of atoms belongs to interface region.

Thanks again.

Hayong.

Dear Axel and Steve.

Thanks for such quick and candid response.
However, I'm having a hard time understanding Axel's comment. Does it
means that I can circumvent the existing problem by increasing a cell
slightly? like, from "3 lattice" to "3.01 lattice". Or, does executing

no that will just change the positions where
rounding/truncation issues can occur. you
will _still_ have lattice positions that are
located _exactly_ on the box boundary, but
floating point numbers are not "exact enough".

simulation by single-core processor is the only solution if I do not
have any other choices in CPU structure?

the simplest method to avoid having lattice positions
located exactly on the simulation cell boundaries
is to slightly translate the box. e.g. in the melt example use:

region box block 0.1 10.1 0.1 10.1 0.1 10.1

instead of:

region box block 0 10 0 10 0 10

And again, (sorry) does Steve's comment on "lattice command" applies
to my second question too? I started looking into the source code on
lattice, but i still cannot understand why altering y direction
influences the number of atoms belongs to interface region.

again, this is due to rounding/truncation issues with floating point
numbers. for example 1/3 cannot be represented and 1/3+1/3+1/3
will result in a number that is smaller than 1.0.

it is almost always better to choose any region boundaries
so that there is no potential for overlap with lattice positions.

axel.