The scattering of X-rays by the atoms in a crystal is described mathematically by the structure factor, $F_{hkl}$. In particular, the intensity of a diffracted beam scattered by a crystal plane with Miller indices $(hkl)$ with a vector $\boldsymbol{K} = h\boldsymbol{\hat{x}}^* + k\boldsymbol{\hat{y}}^* +l\boldsymbol{\hat{z}}^* = 2\pi\left(\frac{h}{a} \boldsymbol{\hat{x}} + \frac{k}{b} \boldsymbol{\hat{y}} + \frac{l}{c} \boldsymbol{\hat{z}}\right)$ is given by the square of the quantity:
$$ F_{hkl} = \sum_j f_j \exp\left( -2\pi i \boldsymbol{K}\cdot \boldsymbol{r_j} \right), $$ where $f_j$ is the form factor of atom $j$ and $\boldsymbol{r_j}$ is its position in the crystal lattice unit cell.
In this exercise, we will focus on crystals that are cubic (so $a=b=c$ and $\boldsymbol{K} = \frac{2\pi}{a}(\boldsymbol{\hat{x}} + \boldsymbol{\hat{y}} + \boldsymbol{\hat{z}})$) and monatomic (so that $f_j$ is the same for all atoms and we will set it equal to unity).
Write a function to calculate $F_{hkl}$ using the above formula for a given crystal structure. Your functions should take two arguments: a sequence of $m$ scattering vectors, hkl
, provided as a NumPy array with dimensions (m,3)
and an array of atom positions within the unit cell with dimensions (n,3)
.
Test your function against the following analytical results:
Body-centred cubic crystals
Atom positions: $$ \textstyle r_1 = (0,0,0) \quad r_2 = (\frac{a}{2},\frac{a}{2},\frac{a}{2}) $$
Structure factor: \begin{align*} F_{hkl} = \left\{ \begin{array}{ll} 2f, & h+k + l \; \mathrm{even}\\ 0, & h+k + l \; \mathrm{odd} \end{array}\right. \end{align*}
For example, this function could be defined by the code:
def F_bcc(hkl):
"""Structure factor for a perfect bcc crystal."""
# F = 2 if h+k+l is even; otherwise F = 0.
return 2 * (~np.sum(hkl, axis=1) % 2)
Here, the (m,3)
array hkl
is summed over axis=1
(the three $h$, $k$, $l$ values) to give an array of shape (m,)
. Each $m$ value is taken modulo 2 to leave the remainder 0 for even values and 1 for odd values. The complement (~
) of this is then 1 for even values and 0 for odd values. We return twice this value as the structure factor for a bcc crystal.
Face-centred cubic crystals
Atom positions: $$ \textstyle r_1 = (0,0,0) \quad r_2 = (\frac{a}{2},\frac{a}{2},0) \quad r_3 = (0, \frac{a}{2},\frac{a}{2}) \quad r_4 = (\frac{a}{2}, 0, \frac{a}{2}) $$
Structure factor: \begin{align*} F_{hkl} = \left\{ \begin{array}{ll} 4f, & h, k, l \; \mathrm{all\;even\;or\;all\;odd}\\ 0, & \mathrm{otherwise} \end{array}\right. \end{align*}
Diamond crystal structure (Harder)
Atom positions $$ \textstyle r_1 = (0,0,0) \quad r_2 = (\frac{a}{4},\frac{a}{4},\frac{a}{4}) \quad r_3 = (\frac{a}{2},\frac{a}{2},0) \quad r_4 = (\frac{3a}{4},\frac{3a}{4},\frac{a}{4}) $$ $$ \textstyle r_5 = (\frac{a}{2},0,\frac{a}{2}) \quad r_6 = (0,\frac{a}{2},\frac{a}{2}) \quad r_7 = (\frac{3a}{4},\frac{a}{4},\frac{3a}{4}) \quad r_8 = (\frac{a}{4},\frac{3a}{4},\frac{3a}{4}) $$
Structure factor: \begin{align*} F_{hkl} = \left\{ \begin{array}{ll} 4f [1 + (-i)^{h+k+l}], & h+k+l \; \mathrm{odd}\\ 8f, & h+k+l = 4n \; \mathrm{(i.e.\;divisible\;by\;4)}\\ 0 & \mathrm{otherwise} \end{array}\right. \end{align*}
The following code tests the NumPy vectorized implementation of the general structure factor summation, calc_F
against the three analytical results for bcc, fcc and diamond crystal structures for $hkl$ vectors with $h,k,l \le 3$.
import numpy as np
# An array of hkl vectors for all combinations of h,k,l up to 3 except (0,0,0)
hkl = np.array(np.meshgrid(*[[0,1,2,3],]*3)).T.reshape(-1,3)[1:]
def calc_F(hkl, atoms):
"""Return the structure factor for K = hkl = (h,k,l) and lattice atoms."""
return np.sum(np.exp(-1j * np.inner(2*np.pi*hkl, atoms)), axis=1)
def calc_absF(hkl, atoms):
"""Return the absolute value of the structure factor, F.
Given a sequence of vectors K = hkl, and lattice positions atoms, return
the absolute value of the structure factor, Fhkl.
"""
return np.abs(calc_F(hkl, atoms))
# Here are the "analytical" structure factors for some cubic crystals.
def F_bcc(hkl):
"""Structure factor for a perfect bcc crystal."""
# F = 2 if h+k+l is even; otherwise F = 0.
return 2 * (~np.sum(hkl, axis=1) % 2)
def F_fcc(hkl):
"""Structure factor for a perfect fcc crystal."""
# F = 4 if h,k,l all odd or all even; otherwise F = 0.
return (np.all(np.isclose(hkl % 2, (0,0,0)), axis=1) |
np.all(np.isclose(hkl % 2, (1,1,1)), axis=1)
) * 4
def F_dia(hkl):
"""Structure factor for a perfect diamond crystal."""
# Unless h,k,l are all odd or all even, F = 0. Otherwise:
# F = 4(1 +(-1)^(h+k+l)) for h+k+l = odd;
# F = 8 for h+k+l = 4n; otherwise F = 0.
all_odd_or_even = (np.all(np.isclose(hkl % 2, (0,0,0)), axis=1) |
np.all(np.isclose(hkl % 2, (1,1,1)), axis=1))
sum_mod_4 = (np.sum(hkl, axis=1) % 4)
sgn = -1 * (sum_mod_4==1) + 1 * (sum_mod_4==3)
return all_odd_or_even * ( (np.isclose(sum_mod_4, 0))*8 +
(sum_mod_4 % 2)*(1 + sgn*1j)*4 )
def check_F(crystal_type, atoms, Ffunc):
# An array of hkl vectors for all combinations of h,k,l up to 3
# except (0,0,0): (1,0,0), (0,1,0), ..., (3,3,2), (3,3,3).
hkl = np.array(np.meshgrid(*[[0,1,2,3],]*3)).T.reshape(-1,3)[1:]
print('### {} ###'.format(crystal_type))
Fhkl = calc_F(hkl, atoms)
Fcheck = Ffunc(hkl)
check = np.isclose(Fhkl, Fcheck)
print('\n'.join( '{}: {} {} {}'.format(*fields) for fields in
zip(hkl, np.round(Fhkl,3), np.round(Fcheck, 3), check) ))
print('OK' if np.all(check) else 'Error')
fcc_atoms = np.array([[0,0,0], [0.5,0.5,0], [0,0.5,0.5], [0.5,0,0.5]])
check_F('fcc', fcc_atoms, F_fcc)
bcc_atoms = np.array([[0,0,0],[0.5,0.5,0.5]])
check_F('bcc', bcc_atoms, F_bcc)
dia_atoms = np.array([[0,0,0], [0.25, 0.25, 0.25], [0.5, 0.5, 0],
[0.75, 0.75, 0.25], [0.5, 0, 0.5], [0, 0.5, 0.5],
[0.75, 0.25, 0.75], [0.25, 0.75, 0.75]])
check_F('dia', dia_atoms, F_dia)