[lammps-users] Triclinic box specification


I have been trying to set up a triclinic box in LAMMPS. To get some help
I searched the mailing list and found the following message by Pieter
in which he discussed the relationship between the conventional specification
(a,b,c,alpha,beta,gamma) and the LAMMPS representation:

The lattice of atoms and the shape of the simulation box are independent
issues. So the lattice with a1,a2,a3 vectors and a basis that you define
has nothing to do with the region prism you use to define the box.
You can make the unit cell the same shape as the box if you like, but
that's up to you. If you viz the simulation box it will be due to the region
prism command, not the atom lattice.

When Pieter writes: h = {xx, yy, zz, yz, xz, xy};
that's how the tilt is stored in the code (Voigt notation),
but the doc pages tell how you to specify it in the input file,
either for the region prism or read_data commands.

We could look at extra input options for the alternate conventions
you mention.


I am aware that the lattice and the "region prism" to define the box are independent. However,
in many cases, you would like the simulation box to be an integer number of lattice unit cells.

The question that I have regarding the order of the tilt distances does concern what is written
in the LAMMPS doc pages.

For the the "region prism" command it says :

prism args = xlo xhi ylo yhi zlo zhi xy xz yz

In my particular case I wish to create a periodic simulation box that has precisely the dimensions of 10 lattice repeat vectors
in each of the 3 directions. The geometry of the lattice, containing 4 basis atoms, is defined by (a,b,c,90,beta,90). This means
that only the angle between the a (x) and c (z) axes is different from 90. This means that the only nonzero
tilt shift would be the xz. According to the order of the tilt displacements in the doc pages, xz should appear as
the second tilt shift. But if I put this number in the second place, the other two being zero, the number of atoms in the
region is NOT 4000, as it should. But if I put the tilt displacement xz in the first place, with the second and third equal to
zero, the region contains exactly 4000 atoms, as expected. Probably I�m doing something wrong, but I have not been able
to figure out what.

Steve Plimpton wrote:

Hmmm. Can you viz the atoms you've created and verify
that the box shape is indeed tilted in xz when you make
the 2nd parameter non-zero? That's what should happen.
Myabe you're not getting 4000 exactly due to round-offs
near the faces of the tilted box. But if you're system is
3d periodic, that shouldn't happen.


I believe I know what is going on. Indeed it appears to be a matter of numerical precision. The
order of the tilt displacements is correct.

If I have the following commands, the program generates a box containing 3840 atoms (wrong).
If you compare the elements of the vector a3 and the zz and xz values in the "region" command,
the latter are exactly 10 times the former.

The attached input script creates 4000 atoms for
me with either set of lattice/region commands commented
in/out. In serial or parallel. This is with the
current, fully patched version of LAMMPS. This
is on a Dell dual quad-core Pentium 64-bit Linux box.
So I'm not seeing any round-off issues.

Are you running the most current LAMMPS?


in.bug.2 (560 Bytes)


I found out what was going on. There was an incorrect command after the "lattice" and "region"

Instead of

create_box 1 firstr
create_atoms 1 box

I had

create_box 1 firstr
create_atoms 1 region firstr

If you use these last commands there appears to be this numerical round-off problem.
If you run the attached script and use the first set of lattice and region commands, 3810 atoms are created.
If you use the second set, where I have added a small number to two numbers in the region command it creates 4000 atoms.

If instead one uses the create_atoms 1 box command, both create 4000 atoms.



Steve Plimpton wrote:

in.bug2 (613 Bytes)

I see the problem now.

When you create atoms, there are two tests that are made.
One if the atom is in the region (if you use the region option).
One if the atom is in the simulation box (always done).

The former test is absolute, which may have round-off issues.
The latter test includes a fudge factor to get things right
across periodic boundaries.

If you want to fill the entire simluation box, don't use the
region option. Normally this wouldn't be a problem, but
with triclinic boxes the round-off issue is unavoidable I think.

To avoid round-off with a region option, you could make
the region slightly larger than needed. But you can't do
that if the region is also used to define the box.


This section from the doc page for "create_atoms" is more consise:

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.