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?


2 Answers 2


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}$
  • $\begingroup$ Great explanation! This should be the accepted answer. $\endgroup$
    – Jan Joswig
    Commented Dec 4, 2020 at 9:38
  • $\begingroup$ Hi! Can I know the source so I can read more about this general approach for PBCs? $\endgroup$ Commented Apr 12, 2021 at 8:31
  • $\begingroup$ Hmmmm, I'm not sure if I have seen this written out formally anywhere outside computer code, but given the definition of fractional coordinates it is just simple manipulation, the above notation is very much my own. Something like Allen and Tildesley "Computer Simulation of Liquids" may well cover it but I can't check, some wonderful person has stolen my copy, grrrrr .... $\endgroup$
    – Ian Bush
    Commented Apr 12, 2021 at 8:40
  • $\begingroup$ @IanBush Alright, not a problem. Will you please be so kind to explain to me how periodic boundary conditions are implemented here simply by wrapping the atoms back into the box. If I had to, say, calculate the distance between each atom, should I not also incorporate the Minimum Image Convention here? If yes, is there a simple way to do it? $\endgroup$ Commented Apr 12, 2021 at 9:36
  • $\begingroup$ That's a separate question - please ask it as such. Comments are not the place for that. $\endgroup$
    – Ian Bush
    Commented Apr 12, 2021 at 10:05

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

  • $\begingroup$ Hi, I know this is really old but can you explain how you implemented PBC for a hexagonal crystal? $\endgroup$ Commented Apr 8, 2021 at 13:18
  • $\begingroup$ @pauling0304, you should ask another question with the details for your problem. $\endgroup$
    – nicoguaro
    Commented Apr 8, 2021 at 14:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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