P6.X.4Structure factors for perfect crystals

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*} $$

Show solution

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

python

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)

The output of this program confirms agreement between the explicit and analytical formulas for calculating the structure factor:

### fcc ###
[0 1 0]: -0j 0 True
[0 2 0]: (4+0j) 4 True
[0 3 0]: -0j 0 True
[1 0 0]: -0j 0 True
[1 1 0]: 0j 0 True
[1 2 0]: -0j 0 True
[1 3 0]: 0j 0 True
[2 0 0]: (4+0j) 4 True
[2 1 0]: -0j 0 True
[2 2 0]: (4+0j) 4 True
[2 3 0]: -0j 0 True
[3 0 0]: -0j 0 True
[3 1 0]: 0j 0 True
[3 2 0]: -0j 0 True
[3 3 0]: 0j 0 True
[0 0 1]: -0j 0 True
[0 1 1]: 0j 0 True
[0 2 1]: -0j 0 True
[0 3 1]: 0j 0 True
[1 0 1]: 0j 0 True
[1 1 1]: (4+0j) 4 True
[1 2 1]: -0j 0 True
[1 3 1]: (4+0j) 4 True
[2 0 1]: -0j 0 True
[2 1 1]: -0j 0 True
[2 2 1]: -0j 0 True
[2 3 1]: -0j 0 True
[3 0 1]: 0j 0 True
[3 1 1]: (4+0j) 4 True
[3 2 1]: -0j 0 True
[3 3 1]: (4+0j) 4 True
[0 0 2]: (4+0j) 4 True
[0 1 2]: -0j 0 True
[0 2 2]: (4+0j) 4 True
[0 3 2]: -0j 0 True
[1 0 2]: -0j 0 True
[1 1 2]: -0j 0 True
[1 2 2]: -0j 0 True
[1 3 2]: -0j 0 True
[2 0 2]: (4+0j) 4 True
[2 1 2]: -0j 0 True
[2 2 2]: (4+0j) 4 True
[2 3 2]: -0j 0 True
[3 0 2]: -0j 0 True
[3 1 2]: -0j 0 True
[3 2 2]: -0j 0 True
[3 3 2]: -0j 0 True
[0 0 3]: -0j 0 True
[0 1 3]: 0j 0 True
[0 2 3]: -0j 0 True
[0 3 3]: 0j 0 True
[1 0 3]: 0j 0 True
[1 1 3]: (4+0j) 4 True
[1 2 3]: -0j 0 True
[1 3 3]: (4+0j) 4 True
[2 0 3]: -0j 0 True
[2 1 3]: -0j 0 True
[2 2 3]: -0j 0 True
[2 3 3]: -0j 0 True
[3 0 3]: 0j 0 True
[3 1 3]: (4+0j) 4 True
[3 2 3]: -0j 0 True
[3 3 3]: (4+0j) 4 True
OK
### bcc ###
[0 1 0]: -0j 0 True
[0 2 0]: (2+0j) 2 True
[0 3 0]: -0j 0 True
[1 0 0]: -0j 0 True
[1 1 0]: (2+0j) 2 True
[1 2 0]: -0j 0 True
[1 3 0]: (2+0j) 2 True
[2 0 0]: (2+0j) 2 True
[2 1 0]: -0j 0 True
[2 2 0]: (2+0j) 2 True
[2 3 0]: -0j 0 True
[3 0 0]: -0j 0 True
[3 1 0]: (2+0j) 2 True
[3 2 0]: -0j 0 True
[3 3 0]: (2+0j) 2 True
[0 0 1]: -0j 0 True
[0 1 1]: (2+0j) 2 True
[0 2 1]: -0j 0 True
[0 3 1]: (2+0j) 2 True
[1 0 1]: (2+0j) 2 True
[1 1 1]: -0j 0 True
[1 2 1]: (2+0j) 2 True
[1 3 1]: -0j 0 True
[2 0 1]: -0j 0 True
[2 1 1]: (2+0j) 2 True
[2 2 1]: -0j 0 True
[2 3 1]: (2+0j) 2 True
[3 0 1]: (2+0j) 2 True
[3 1 1]: -0j 0 True
[3 2 1]: (2+0j) 2 True
[3 3 1]: -0j 0 True
[0 0 2]: (2+0j) 2 True
[0 1 2]: -0j 0 True
[0 2 2]: (2+0j) 2 True
[0 3 2]: -0j 0 True
[1 0 2]: -0j 0 True
[1 1 2]: (2+0j) 2 True
[1 2 2]: -0j 0 True
[1 3 2]: (2+0j) 2 True
[2 0 2]: (2+0j) 2 True
[2 1 2]: -0j 0 True
[2 2 2]: (2+0j) 2 True
[2 3 2]: -0j 0 True
[3 0 2]: -0j 0 True
[3 1 2]: (2+0j) 2 True
[3 2 2]: -0j 0 True
[3 3 2]: (2+0j) 2 True
[0 0 3]: -0j 0 True
[0 1 3]: (2+0j) 2 True
[0 2 3]: -0j 0 True
[0 3 3]: (2+0j) 2 True
[1 0 3]: (2+0j) 2 True
[1 1 3]: -0j 0 True
[1 2 3]: (2+0j) 2 True
[1 3 3]: -0j 0 True
[2 0 3]: -0j 0 True
[2 1 3]: (2+0j) 2 True
[2 2 3]: -0j 0 True
[2 3 3]: (2+0j) 2 True
[3 0 3]: (2+0j) 2 True
[3 1 3]: -0j 0 True
[3 2 3]: (2+0j) 2 True
[3 3 3]: -0j 0 True
OK
### dia ###
[0 1 0]: (-0-0j) 0j True
[0 2 0]: -0j 0j True
[0 3 0]: -0j 0j True
[1 0 0]: (-0-0j) 0j True
[1 1 0]: 0j 0j True
[1 2 0]: -0j 0j True
[1 3 0]: 0j 0j True
[2 0 0]: -0j 0j True
[2 1 0]: -0j 0j True
[2 2 0]: (8+0j) (8+0j) True
[2 3 0]: (-0-0j) 0j True
[3 0 0]: -0j 0j True
[3 1 0]: 0j 0j True
[3 2 0]: (-0-0j) 0j True
[3 3 0]: -0j 0j True
[0 0 1]: (-0-0j) 0j True
[0 1 1]: 0j 0j True
[0 2 1]: -0j 0j True
[0 3 1]: -0j 0j True
[1 0 1]: 0j 0j True
[1 1 1]: (4+4j) (4+4j) True
[1 2 1]: -0j 0j True
[1 3 1]: (4-4j) (4-4j) True
[2 0 1]: -0j 0j True
[2 1 1]: -0j 0j True
[2 2 1]: (-0-0j) 0j True
[2 3 1]: -0j 0j True
[3 0 1]: 0j 0j True
[3 1 1]: (4-4j) (4-4j) True
[3 2 1]: -0j 0j True
[3 3 1]: (4+4j) (4+4j) True
[0 0 2]: -0j 0j True
[0 1 2]: -0j 0j True
[0 2 2]: (8+0j) (8+0j) True
[0 3 2]: (-0-0j) 0j True
[1 0 2]: -0j 0j True
[1 1 2]: -0j 0j True
[1 2 2]: (-0-0j) 0j True
[1 3 2]: -0j 0j True
[2 0 2]: (8+0j) (8+0j) True
[2 1 2]: (-0-0j) 0j True
[2 2 2]: -0j 0j True
[2 3 2]: -0j 0j True
[3 0 2]: (-0-0j) 0j True
[3 1 2]: -0j 0j True
[3 2 2]: -0j 0j True
[3 3 2]: -0j 0j True
[0 0 3]: -0j 0j True
[0 1 3]: -0j 0j True
[0 2 3]: (-0-0j) 0j True
[0 3 3]: -0j 0j True
[1 0 3]: -0j 0j True
[1 1 3]: (4-4j) (4-4j) True
[1 2 3]: 0j 0j True
[1 3 3]: (4+4j) (4+4j) True
[2 0 3]: (-0-0j) 0j True
[2 1 3]: 0j 0j True
[2 2 3]: (-0-0j) 0j True
[2 3 3]: -0j 0j True
[3 0 3]: 0j 0j True
[3 1 3]: (4+4j) (4+4j) True
[3 2 3]: -0j 0j True
[3 3 3]: (4-4j) (4-4j) True
OK