8
$\begingroup$

The minimum image convention (MIC), see for example a short note of W. Smith, is often used in molecular dynamics or monte carlo simulations of periodic systems with an orthorhombic unit cell. For this special case, it is rather trivial to implement the MIC correctly. How does one apply the MIC to systems with a more general triclinic unit cell? The naive idea to simply generalize the algorithm for orthorhomic cells to triclinic cells, does not seem to work.

I could find this document, which provides some clues, such as the importance of the reduced basis of the lattice vectors, but it does not contain enough details to make a computer implementation.

I will implement the answers below in a test code to check which algorithms work (under given assumptions). The Python test program can be found here: https://gist.github.com/3566972. It tests the lammps implementation and one naive attempt of mine. Both fail. If one finds a bug in this test program, please get in touch.

P.S. For those who get worried about MD or MC codes that use a naive algorithm for triclinic cells: please read the note written by W. Smith. It explains a naive algorithm that works as long as the distance cutoff is shorter than half the shortest distance between the opposite faces of the unit cell. For a given cutoff, one may always construct a triclinic supercell for which naive algorithms will work fine.

$\endgroup$
1
  • $\begingroup$ If this question (or similar) is still of interest to you, you may want to check out Matter Modeling Stack Exchange. $\endgroup$
    – Tyberius
    Apr 16 at 17:36

3 Answers 3

5
$\begingroup$

The key is to take the differences $\Delta x$, $\Delta y$, and $\Delta z$ separately before beginning. Given the edge vectors ${\bf a}$, ${\bf b}$, and ${\bf c}$ that define a unit cell that has one corner at the origin, then your tilt factors are $b_1$, $c_1$, and $c_2$, where $n_m$ defines the $m^{th}$ component of vector ${\bf n}$.

Now, for each component:

  • If $|\Delta x| > 0.5(a_{max} - a_{min})$, then subtract off $a_{max} - a_{min}$ if $\Delta x > 0$, and add it for values less than zero.
  • If $|\Delta y| > 0.5(b_{max} - b_{min})$, then subtract off $b_{max} - b_{min}$ if $\Delta y > 0$, and add it for values less than zero. Also add or subtract $b_1$ from $\Delta x$.
  • If $|\Delta z| > 0.5(c_{max} - c_{min})$, then subtract off $c_{max} - c_{min}$ if $\Delta z > 0$, and add it for values less than zero. Also add or subtract $c_1$ from $\Delta x$ and $c_2$ from $\Delta y$.

My source for this is taken from the code for triclinic domains found in domain.cpp in the LAMMPS molecular dynamics package.

$\endgroup$
2
  • 1
    $\begingroup$ Thanks for the reply. Unfortunately, this algorithms fails for some cases (see test program). I guess it will (also) only work if the cutoff for the relative vectors is smaller than half the distance between two crystal planes. It is still an interesting answer, because it shows that most MD codes do not bother to implemented a more rigorous (and computationally more demanding) solution. $\endgroup$ Sep 1, 2012 at 8:30
  • $\begingroup$ I am not sure I fully understood your answer. Can you clafify what you mean when you say "* If $|\Delta x| > 0.5(a_{max} - a_{min})$, then subtract off $a_{max} - a_{min}$ if $\Delta x > 0$, and add it for values less than zero." The second if statment, is it dependent on the first being true? When you say "and add it" what do you mean by "and" there doesn't seem to be anything to "and" there. $\endgroup$
    – Magpie
    Jan 2, 2013 at 17:37
1
$\begingroup$

The easiest way to do this is to use the fractional coordinate. The algorithm looks like this:

  • Transform coordinate from cartesian to fractional;
  • Round fractional coordinates to [0, 1);
  • Transform back rounded coordinates to cartesian.

To transform from cartesian coordinates to fractional, you have to multiply the coordinate vector by the inverted cell matrix: $$ H^{-1} = \begin{pmatrix} a & b \cos \gamma & c \cos \beta\\ 0 & b \sin \gamma & c (\cos\alpha - \cos\beta \cos\gamma)/\sin\gamma\\ 0 & 0 & \sqrt{c^2 - c^2 \cos^2 \beta - c^2 (\cos\alpha - \cos\beta \cos\gamma)^2/\sin^2\gamma} \end{pmatrix}^{-1}$$

And to transform back to cartesian, you have to multiply the vector by the cell matrix $H$.

The good thing is that this also works for orthorombic cells.

$\endgroup$
1
  • 1
    $\begingroup$ This is what I called "the naive implementation" in my tests, and it does not work in general. This is known to fail sometimes for non-orthorhombic cells when the shortest distance between two atoms exceeds half of the distance between two faces of the unit cell. This algorithm may then produce not the shortest possible relative vector, but a longer one. In the orthorhomic case, it does work properly. $\endgroup$ Jul 11, 2017 at 5:45
0
$\begingroup$

You have box vector as (a, b, c) with angles $(\alpha, \beta, \gamma)$. That can be converted to the box vector as | a :0 :0| | xy:b :0| | xz:yz:c| Where xy- change in x due to y, xz-change in x due to z, yz- change in y due to z. To calculate minimum distance. First calculate rij = $\Delta x + \Delta y + \Delta z$, distance between $i^{th}$ and $j^{th}$ particle.

Now, calculate the number of boxes that needs to be added or deleted. That is dividing the distance with trace of the box vector. N_box=nearest integer(rij/trace(box_vector)). For x, N_box(1)=nearest_integer($\Delta x/a$).

Calculate minimum distance: rij = rij - N_box*box_vector.

That is; $\Delta x = \Delta x$ - (N_box(1) * a + N_box(2) * xy +N_box(3) * xz )

$\Delta y = \Delta y$ - (N_box(2) * b +N_box(3) * yz )

$\Delta z = \Delta z$ - N_box(3) * c.

This should generally work, but in rare cases, when particles are in the corner of the box, then the, updated values of $\Delta z$ or $\Delta y$ may affect $\Delta y$ or $\Delta x$. Therefore for the safer side, I prefer to calculate z to x.

N_box(3) = $\Delta z / c $ ; $\Delta z = \Delta z$ - N_box(3) * c

N_box(2) = ($\Delta y$ + N_box(3)*yz) /b ; $\Delta y = \Delta y$ - N_box(2) * b

N_box(1) = ($\Delta x$ + N_box(3)*xz + N_box(2)*xy) /a ; $\Delta x = \Delta x$ - N_box(1) * a

This will work for all cases. You are welcome to post optimized vector calculations for the steps.

$\endgroup$

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.