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.



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?


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

enter image description here

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).

| cite | improve this answer | |

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}$
| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.