# Efficient way to generate a crystal lattice in a box

I'm looking for a fast way to generate a crystal lattice. I know you can do a bunch of linear combinations of lattice vectors, but this can be costly. Here's my current code for doing so in Python, for a FCC lattice:

import numpy as np

a = 5.472 # lattice parameter

# Lattice vector matrix
A = a*np.array([[0.0, 0.5, 0.5],
[0.5, 0.0, 0.5],
[0.5, 0.5, 0.0]] )

# Basis in direct coordinates
basisd = np.array ( [[0.0, 0.0, 0.0]] )

# Define the orthohombic box
bx = np.array([ [a, 0., 0.],
[0., a, 0.],
[0., 0., a] ] )

# Define the supercell dimensions
supercell = np.array([3., 3., 3.])

# Define maximum scope of n values to check with r= n1*a1 + n2*a2 + n3*a3 (use a larger number for a larger lattice)
nvals = 10;

# Convert to cartesian coordinates
basisc = np.matmul(basisd, A)
bx = np.multiply(bx, supercell)
bx_inv = np.linalg.inv(bx)

a1 = A[0]
a2 = A[1]
a3 = A[2]

basis_shape = np.shape(basisc)
batoms = basis_shape[0]
pos = []
nums = []
for b in range(0,batoms):
for n1 in range(-nvals, nvals+1):
for n2 in range(-nvals, nvals+1):
for n3 in range(-nvals, nvals+1):
r = n1*a1 + n2*a2 + n3*a3 + basisc[b]

if 0 <= r[0] and r[0] < bx[0][0] and 0 <= r[1] and r[1] < bx[1][1] and 0 <= r[2] and r[2] < bx[2][2]:

pos.append(r)


I had to throw an if-statement in that loop to ensure that I'm only capturing positions that are inside the box, since I'm sampling many positions outside the box. The final array pos has all the lattice positions in the box. This code is fine for a small 3X3X3 unit cell, but it is incredibly slow for larger cells, where I have to sample more lattice vector linear combinations.

Instead of looping over many integers $n_1, n_2, n_3$ to calculate positions $\textbf{r} = n_1\textbf{a}_1+n_2\textbf{a}_2 + n_3\textbf{a}_3$, is there a faster way to do this?

• This is a question about coding, not physics. Computational Science SE would be a better place to ask such questions. May 11, 2018 at 19:50
• You may also want to look at Code Review for some Python-oriented pointers. May 11, 2018 at 19:56
• I know people use Python for coding but honestly it's dreadfully inefficient when you have to use Numpy. At least try PyPy or better yet Java or C++. May 11, 2018 at 20:39
• @StephenG That nested Python loop is the killer. Generally, if you're using explicit Python loops on Numpy arrays you're doing it wrong, although sometimes such loops are unavoidable. Numpy is fine when it can handle the looping itself at compiled speed, and it's performing arithmetic with machine numeric types. FWIW, converting machine numbers to Python objects is relatively slow, extracting the machine numeric value from an object is quick. May 12, 2018 at 21:46

You don't seem to be using the symmetry of the crystal to your advantage. In fact you could get away to simply compute coordinates for a single $n_1n_2$ plane, and then make N copies along the $n_3$ direction while adding $\bf{a}_3$ in every iteration to the whole plane