periodic boundary conditions for triclinic box

Question

I am trying to do analysis on a data set of atomic coordinates generated form lammps. I simulated an alpha glycine crystal in a triclinic box. The box vectors look like the following, where xy,xz,and yz, are the tilt factors.

    xlo = -0.34; xhigh = 40.466
    ylo = -1.6283; yhigh = 49.50861
    zlo = 0.77; zhigh = 39.79
    xy  = -8.62; xz = 0.0; yz = 0.0

In the case of a cubic box, the way to implement periodic boundary conditions is trivial.

if(x.gt.xhigh)x=x-xbox
if(x.lt.xlo)x=x+box

etc.

However, I cannot seem to figure out how to handle the tilt factors, for the triclinic box, to implement the PBC's correctly . Could someone please provide a code snippet outlining this for me?

Answer 1

Well there's always the general approach. Let's say we have a particle at a position $\bf r$ and the lattice vectors, i.e. the vectors that forming the sides of the simulation cell, are $\bf a_i$. We can then express $\bf r$ in the basis of $\bf a_i$:

$${\bf r}=\sum_i f_i{\bf a_i}$$

You can express this as a matrix vector multiply $${\bf r}={\bf Af}$$ where A is a matrix whose columns are the lattice vectors, and it can be see that f are the fractional coordinates. So defining ${\bf B}={\bf A}^{-1}$ you can find the factional coordinates simply by

$${\bf f}={\bf Br}$$

Now we want all our images to be in the same cell ... One way of doing this is to ensure $0\le f_i < 1$, and to achieve this we can simply modify $f_i$ to be $f_i - floor(f_i)$ as this is simply translates the particles to equivalent position in the reference cell due to the periodic boundary conditions. So for arbitrary cell (in arbitrary dimensionality) the recipe is

• Form ${\bf A}$
• Invert ${\bf A}$ to form ${\bf B}$, this should always be well determined numerically as the determinant of the matrix is the volume of the cell, and so should be far from zero. You also probably want these vectors anyway for the Ewald sum ...
• From the position of a particle at $\bf r$ form $\bf f$ by ${\bf f}={\bf Br}$
• Translate into the reference cell by $g_i=f_i-floor(f_i)$
• Form the new position in real space ${\bf t}$ by ${\bf t}={\bf Ag}$

Answer 2

I have never implemented periodicity for this symmetry class, although I have done it for hexagonal crystals before.

Let us start with a (2D) box like the one depicted in the image

In this case, we have one axis aligned with $x$ and the other one with the vector $u_2$. The sides of the box are given by $y_\text{low}$, $y_\text{high}$, and the oblique line $ax + by=c_1$, $ax + by=c_2$. And make the sides $A$ and $B$.

Then, when a particle goes above $y_\text{high}$ you reassign it with

$$(x,y) = (x,y) - B u_2$$

in the other case, when $y < y_\text{low}$

$$(x,y) = (x,y) + B u_2 \enspace .$$

In the $x$ direction the new position is easier is just $$(x,y) = (x,y) - A u_1$$ and $$(x,y) = (x,y) + A u_1 \enspace .$$

The difficulty in this case is knowing when the particle goes out of the boundary. In the example above, the particle goes out to the right if

$$ax + by < c_2$$

and goes out to the left if

$$ax + by \geq c_1 \enspace .$$

Thus, the outside lies below the line $ax + by = c_2$ if the shear of the box is to the right and above if it is to the left. This reasoning can be extended to more dimensions (although in your particular example there is just one shear in the box).