Crystallography

Anri’s crystallography is a set of plain functions on arrays. A structure with atoms is a Dans_Diffraction.Crystal, which also reads the CIF; everything else needs only the lattice parameters and the space-group number.

import anri.utils
anri.utils.setup()  # before any JAX computation: don't grab most of the GPU memory at start-up

import Dans_Diffraction
import numpy as np

import anri.crystal

xtl = Dans_Diffraction.Crystal("../../../tests/data/cif/Fe.cif")
lpars = anri.crystal.lattice_parameters(xtl)
sg = anri.crystal.space_group(xtl)
lpars, sg

(array([ 2.8630355,  2.8630355,  2.8630355, 90.       , 90.       ,
        90.       ]),
 229)

The B matrix (with \(g = U B h\), no \(2\pi\)). B_matrix evaluates the JAX definition, anri.crystal.lpars_to_B, in float64, so there is one definition that can also be jitted and differentiated:

B = anri.crystal.B_matrix(lpars)
B

array([[ 3.49279637e-01, -5.61684495e-17, -5.61684495e-17],
       [ 0.00000000e+00,  3.49279637e-01, -2.13872095e-17],
       [ 0.00000000e+00,  0.00000000e+00,  3.49279637e-01]])

Now the reflections up to a \(d^*\) limit, with the space group’s systematic absences removed. The wavelength gives their \(2\theta\), and caps \(d^*\) at \(2/\lambda\). They come sorted by \(d^*\):

dsmax = 1.0
wavelength = 0.18
refl = anri.crystal.reflections(lpars, sg, wavelength, dsmax)
refl["hkl"][:8], refl["ds"][:8], refl["tth"][:8]

(array([[-1, -1,  0],
        [-1,  0, -1],
        [-1,  0,  1],
        [-1,  1,  0],
        [ 0, -1, -1],
        [ 0, -1,  1],
        [ 0,  1, -1],
        [ 0,  1,  1]]),
 array([0.493956, 0.493956, 0.493956, 0.493956, 0.493956, 0.493956,
        0.493956, 0.493956]),
 array([5.09596643, 5.09596643, 5.09596643, 5.09596643, 5.09596643,
        5.09596643, 5.09596643, 5.09596643]))

Grouping them into rings of equal \(d^*\) gives each reflection a ring index, and each ring its \(d^*\); counting the members gives the multiplicities:

ring, ring_ds = anri.crystal.rings(refl["ds"], tol=1e-4)
ring_tth = np.degrees(2 * np.arcsin(ring_ds * wavelength / 2))
ring_mult = np.bincount(ring)
ring_ds, ring_tth, ring_mult

(array([0.493956  , 0.69855927, 0.85555689, 0.987912  ]),
 array([ 5.09596643,  7.20916425,  8.83230631, 10.20204589]),
 array([12,  6, 24, 12]))

refl["hkl"][ring == 1]  # the {200} ring

array([[-2,  0,  0],
       [ 0, -2,  0],
       [ 0,  0, -2],
       [ 0,  0,  2],
       [ 0,  2,  0],
       [ 2,  0,  0]])

Because we have a structure with atoms, we can compute \(|F|^2\) for each reflection, with Debye-Waller factors and anomalous dispersion at this wavelength. This is where reflections that the space group allows but the atom positions extinguish (e.g. {222} in diamond) show up, as \(|F|^2 \approx 0\):

F2 = anri.crystal.structure_factors(xtl, refl["hkl"], wavelength)
F2_ring = np.array([F2[ring == r].mean() for r in range(len(ring_ds))])
F2_ring

/tmp/ipykernel_4785/689639548.py:1: UserWarning: No isotropic thermal factors (U_iso or B_iso) for Fe: intensities have no Debye-Waller attenuation.
  F2 = anri.crystal.structure_factors(xtl, refl["hkl"], wavelength)

array([1393.51724348,  955.07679029,  710.99911354,  560.11383862])

The space group’s operations, as 4 x 4 matrices on fractional coordinates, and the rotations of its Laue group in the Cartesian crystal frame (for misorientations and orientation grids):

S = anri.crystal.symmetry_matrices(sg)
ops = anri.crystal.laue_rotations(S, B)
S.shape, ops.shape

((96, 4, 4), (24, 3, 3))

We can confirm the rings against ImageD11 too:

from ImageD11.unitcell import unitcell

uc = unitcell(lpars, sg)
uc.makerings(dsmax, tol=1e-4)
print(np.abs(ring_ds - np.array(uc.ringds)).max())
print(all(sorted(map(tuple, refl["hkl"][ring == i].tolist())) == sorted(uc.ringhkls[ds]) for i, ds in enumerate(uc.ringds)))

2.220446049250313e-16
True