This page was generated from doc/tutorials/crystal_directions.ipynb. Interactive online version: .

Crystal directions

In this tutorial we will perform operations with and plot directions with resepct to the crystal reference frame using Miller indices with the orix.vector.Miller class.

Many of the initial examples, explaining basic crystallographic computations with crystal lattice directions \([uvw]\) and crystal lattice planes \((hkl)\), are taken from the textbook Introduction to Conventional Transmission Electron Microscopy ([De Graef, 2003]). Some of the later examples are also inspired by MTEX’ excellent documentation on Miller indices and operations with them.

%matplotlib inline

from diffpy.structure import Lattice, Structure
import matplotlib.pyplot as plt
import numpy as np
from orix.crystal_map import Phase
from orix.quaternion import Orientation, Rotation, symmetry
from orix.vector import Miller, Vector3d


plt.rcParams.update(
    {
        "figure.figsize": (7, 7),
        "font.size": 20,
        "axes.grid": True,
        "lines.markersize": 10,
        "lines.linewidth": 2,
    }
)

To start with, let’s create a tetragonal crystal with lattice parameters \(a\) = \(b\) = 0.5 nm and \(c\) = 1 nm and symmetry elements described by point group 4

tetragonal = Phase(
    point_group="4",
    structure=Structure(lattice=Lattice(0.5, 0.5, 1, 90, 90, 90)),
)
print(tetragonal)
print(tetragonal.structure.lattice)

<name: . space group: None. point group: 4. proper point group: 4. color: tab:blue>
Lattice(a=0.5, b=0.5, c=1, alpha=90, beta=90, gamma=90)

Here, the Phase class attaches a point group symmetry, Symmetry, to a Structure containing a Lattice (where the crystal axes are defined), and possibly some Atoms.

Directions \([uvw]\)

A crystal lattice direction \(\mathbf{m} = u \cdot \mathbf{a} + v \cdot \mathbf{b} + w \cdot \mathbf{c}\) is a vector with coordinates \(u, v, w\) with respect to the crystal axes \(\mathbf{a}, \mathbf{b}, \mathbf{c}\), and is denoted \([uvw]\). We can create a set of these Miller indices by passing a list/array/tuple to uvw in Miller

m1 = Miller(uvw=[[1, 2, 0], [3, 1, 1]], phase=tetragonal)
m1.scatter(axes_labels=["x", "y", None], c=["b", "r"])
m1

Miller (2,), point group 4, uvw
[[1. 2. 0.]
 [3. 1. 1.]]

Here, we plotted the lattice directions in the stereographic projection using the Vector3d.scatter() plotting method, which also works for Miller.scatter() since the Miller class inherits most of the functionality in the Vector3d class.

Let’s compute the dot product in nanometres and the angle in degrees between the vectors

m1[0].dot(m1[1])

array([1.25])

m1[0].angle_with(m1[1], degrees=True)

array([53.3007748])

The length of a direct lattice vector is available via the Miller.length property, given in lattice parameter units (nm in this case)

Miller(uvw=[0, -0.5, 0.5], phase=tetragonal).length

array([0.55901699])

Planes \((hkl)\)

A crystal lattice plane \((hkl)\) is described by its normal vector \(\mathbf{n} = h\cdot\mathbf{a^*} + k\cdot\mathbf{b^*} + l\cdot\mathbf{c^*}\), where \(\mathbf{a^*}, \mathbf{b^*}, \mathbf{c^*}\) are the reciprocal crystal axes. As for crystal directions \([uvw]\), we can create a set of these Miller indices by passing hkl instead of uvw to Miller

m2 = Miller(hkl=m1.uvw, phase=tetragonal)
m2.scatter(c=["y", "g"], marker="*")
m2

Miller (2,), point group 4, hkl
[[1. 2. 0.]
 [3. 1. 1.]]

Let’s compute the angle in degrees between the lattice plane normals

m2[0].angle_with(m2[1], degrees=True)

array([45.69879831])

Note that the lattice plane normals \((hkl)\) are not always parallel to the lattice directions \([uvw]\), even if the indices are the same

fig = m1.scatter(return_figure=True, c=["b", "r"])
m2.scatter(figure=fig, c=["y", "g"], marker="*")

We can get the reciprocal components of the lattice vector [114] (i.e. which lattice plane the [114] direction is perpendicular to) by accessing Miller.hkl for a Miller object with crystal directions (or Miller.uvw for a Miller object with crystal plane normals)

Miller(uvw=[1, 1, 4], phase=tetragonal).hkl

array([[0.25, 0.25, 4.  ]])

We see that the lattice vector [114] is perpendicular to the lattice plane (1 1 16).

The length of reciprocal lattice vectors can also be accessed via the Miller.length property, and equals \(\left|\mathbf{g}_{\mathrm{hkl}}\right| = \frac{1}{d_{\mathrm{hkl}}}\) in reciprocal lattice parameter units (nm^-1 in this case), meaning we can obtain the interplanar spacing \(d_{\mathrm{hkl}}\)

m2.length

array([4.47213595, 6.40312424])

1 / m2.length

array([0.2236068 , 0.15617376])

Zone axes

The cross product of the lattice vectors [110] and [111] in the tetragonal crystal, in direct space, is described by a vector expressed in reciprocal space, \((hkl)\)

m3 = Miller(uvw=[[1, 1, 0], [1, 1, 1]], phase=tetragonal)
m3perp = m3[0].cross(m3[1])
m3perp

Miller (1,), point group 4, hkl
[[ 0.25 -0.25  0.  ]]

The exact “indices” are returned, however, we can get a new Miller instance with the indices rounded down or up to the closest smallest integer via the Miller.round() method

m3perp.round()

Miller (1,), point group 4, hkl
[[ 1. -1.  0.]]

The maximum index that Miller.round() might return can be set by passing the max_index parameter.

We can plot these direct lattice vectors \([uvw]\) and the vectors normal to the cross product vector, i.e. normal to the reciprocal lattice vector \((hkl)\)

fig = m3.scatter(return_figure=True, c="r")
m3perp.draw_circle(figure=fig, color="b", linestyle="--")

Likewise, the cross product of reciprocal lattice vectors (110) and (111) results in a direct lattice vector

m4 = Miller(hkl=[[1, 1, 0], [1, 1, 1]], phase=tetragonal)
m4perp = m4[0].cross(m4[1]).round()
m4perp

Miller (1,), point group 4, uvw
[[ 1. -1.  0.]]

fig = m4.scatter(return_figure=True, c="b")
m4perp.draw_circle(figure=fig, color="r", linestyle="--")

Trigonal and hexagonal index conventions

Crystal lattice vectors and planes in lattices with trigonal and hexagonal crystal symmetry are typically expressed in Weber symbols \([UVTW]\) and Miller-Bravais indices \((hkil)\). The definition of \([UVTW]\) used in orix follows Introduction to Conventional Transmission Electron Microscopy (DeGraef, 2003)

\begin{align} U &= \frac{2u - v}{3},\\ V &= \frac{2v - u}{3},\\ T &= -\frac{u + v}{3},\\ W &= w. \end{align}

For a plane, the \(h\), \(k\) and \(l\) indices are the same in \((hkl)\) and \((hkil)\), and \(i = -(h + k)\).

The first three Miller indices always sum up to zero, i.e.

\begin{align} U + V + T &= 0,\\ h + k + i &= 0. \end{align}

trigonal = Phase(
    point_group="321",
    structure=Structure(lattice=Lattice(4.9, 4.9, 5.4, 90, 90, 120)),
)
trigonal

<name: . space group: None. point group: 321. proper point group: 32. color: tab:blue>

m5 = Miller(UVTW=[2, 1, -3, 1], phase=trigonal)
m5

Miller (1,), point group 321, UVTW
[[ 2.  1. -3.  1.]]

m6 = Miller(hkil=[1, 1, -2, 3], phase=trigonal)
m6

Miller (1,), point group 321, hkil
[[ 1.  1. -2.  3.]]

fig = m5.scatter(return_figure=True, c="C0", grid_resolution=(30, 30))
m6.scatter(figure=fig, c="C1")

One can change the coordinate format of the Miller class, however this does not change the vector, since all vectors are stored with respect to the cartesian coordinate system internally

print(m6, "\n\n", m6.data)

Miller (1,), point group 321, hkil
[[ 1.  1. -2.  3.]]

 [[0.20408163 0.35347976 0.55555556]]

m6.coordinate_format = "UVTW"
print(m6, "\n\n", m6.data)

Miller (1,), point group 321, UVTW
[[ 0.0278  0.0278 -0.0555  0.1029]]

 [[0.20408163 0.35347976 0.55555556]]

Getting the closest integer indices, however, changes the vector

m6round = m6.round()
print(m6round, "\n\n", m6round.data)

Miller (1,), point group 321, UVTW
[[ 3.  3. -6. 11.]]

 [[22.05       38.19172031 59.4       ]]

Symmetrically equivalent directions and planes

The point group symmetry elements of the crystal lattice can be applied to to describe symmetrically equivalent crystal directions and planes. This applies to crystals in all seven systems, but we’ll use the cubic crystal as an example because of its high symmetry

cubic = Phase(point_group="m-3m", structure=Structure())
print(cubic, "\n", cubic.structure.lattice.abcABG())

<name: . space group: None. point group: m-3m. proper point group: 432. color: tab:blue>
 (1.0, 1.0, 1.0, 90.0, 90.0, 90.0)

The directions parallel to the crystal axes (\(\mathbf{a}\), \(\mathbf{b}\), \(\mathbf{c}\)) given by \([100]\), \([\bar{1}00]\), \([010]\), \([0\bar{1}0]\), \([001]\), and \([00\bar{1}]\) (\(\bar{1}\) means “-1”) are symmetrically equivalent, and can be obtained using Miller.symmetrise()

m7 = Miller(uvw=[1, 0, 0], phase=cubic)
m7.symmetrise(unique=True)

Miller (6,), point group m-3m, uvw
[[ 1.  0.  0.]
 [ 0.  1.  0.]
 [-1.  0.  0.]
 [ 0. -1.  0.]
 [ 0.  0.  1.]
 [ 0.  0. -1.]]

Without passing unique=True, since the cubic crystal symmetry is described by 48 symmetry operations (or elements), 48 directions would have been returned.

These six directions, known as a family, may be expressed collectively as \(\left<100\right>\), the brackets implying all six permutations or variants of 1, 0, 0. This particular family is said to have a multiplicity of 6

m7.multiplicity

array([6])

m8 = Miller(uvw=[[1, 0, 0], [1, 1, 0], [1, 1, 1]], phase=cubic)
m8

Miller (3,), point group m-3m, uvw
[[1. 0. 0.]
 [1. 1. 0.]
 [1. 1. 1.]]

m8.multiplicity

array([ 6, 12,  8])

Let’s plot the symmetrically equivalent directions from the direction families \(\left<100\right>\), \(\left<110\right>\), and \(\left<111\right>\) impinging on the upper hemisphere. By also returning the indices of which family each symmetrically equivalent direction belongs to from Miller.symmetrise(), we can give a unique color per family

m9, idx = m8.symmetrise(unique=True, return_index=True)

fig = m9[idx == 0].scatter(c="C0", return_figure=True)
for i in range(1, m9.size):
    m9[idx == i].scatter(c=f"C{i}", figure=fig)

Similarly, symmetrically equivalent planes \((hkl)\) can be collectively expressed as planes of the form \(\{hkl\}\)

m10 = Miller(hkl=[[1, 0, 0], [1, 1, 0], [1, 1, 1]], phase=cubic)
m10.multiplicity

array([ 6, 12,  8])

fig = m10[0].symmetrise(unique=True).scatter(c="C0", return_figure=True)
for i in range(1, m10.size):
    m10[i].symmetrise(unique=True).scatter(c=f"C{i}", figure=fig)

We computed the angles between directions and plane normals earlier in this tutorial. In general, Miller.angle_with() does not consider symmetrically equivalent directions, unless use_symmetry=True is passed. Consider \((100)\) and \((\bar{1}00)\) and \((111)\) and \((\bar{1}11)\) in the stereographic plot above

m11 = Miller(hkl=[[1, 0, 0], [1, 1, 1]], phase=cubic)
m12 = Miller(hkl=[[-1, 0, 0], [-1, 1, 1]], phase=cubic)

m11.angle_with(m12, degrees=True)

array([180.        ,  70.52877937])

m11.angle_with(m12, use_symmetry=True, degrees=True)

array([0., 0.])

Thus, passing use_symmetry=True ensures that the smallest angles between m11 and the symmetrically equivalent directions to m12 are found.## Directions and planes in rotated crystals

Directions and planes in rotated crystals

Let’s consider the orientation of a cubic crystal rotated \(45^{\circ}\) about the sample \(\mathbf{Z}\) axis. Orientations in orix are defined as rotations from the sample to the crystal (so-called “lab2crystal”), so we have to invert the orientation to get the vector expressed in the sample reference frame. Furthermore, as per the discussion by [Rowenhorst et al., 2015], the axis-angle representation of a \(45^{\circ}\) rotation about the sample \(\mathbf{Z}\) is written as \((\mathbf{\hat{n}}, \omega) = ([00\bar{1},90^{\circ})\)

o = Orientation.from_axes_angles(
    (0, 0, -1), 45, symmetry=cubic.point_group, degrees=True
)
o

Orientation (1,) m-3m
[[ 0.9239  0.      0.     -0.3827]]

We can apply this orientation to a crystal direction \([uvw]\) to find this direction in this particular crystal with respect to the sample coordinate system

m12 = Miller(uvw=[1, 1, 1], phase=cubic)
~o * m12

Miller (1,), point group m-3m, uvw
[[-0.      1.4142  1.    ]]

# [uvw] in unrotated crystal
fig = m12.scatter(c="r", return_figure=True, axes_labels=["X", "Y", "Z"])

# [uvw] in rotated crystal
(~o * m12).scatter(figure=fig, c="b", marker="x")

We see that the \([111]\) vector in the crystal is aligned with the sample \(\mathbf{Y}\) direction.

We can apply the crystal symmetry to obtain the coordinates with respect to the sample reference frame for all crystallographically equivalent directions

(~o * m12.symmetrise(unique=True)).scatter(c="b", marker="x")

o2 = Orientation.from_euler(
    [10, 20, 30], symmetry=trigonal.point_group, degrees=True
)
o2

Orientation (1,) 321
[[ 0.9254 -0.171   0.0302 -0.3368]]

m13 = Miller(hkil=[1, -1, 0, 0], phase=trigonal)
m13

Miller (1,), point group 321, hkil
[[ 1. -1. -0.  0.]]

p = ~o2 * m13.symmetrise(unique=True)
p.scatter(
    hemisphere="both",
    grid_resolution=(30, 30),
    figure_kwargs=dict(figsize=(10, 5)),
    axes_labels=["X", "Y", "Z"],
)

The stereographic plots above are essentially the \(\{1\bar{1}00\}\) pole figure representation of the orientation \(O_2\).

---

A diamond [111] pole figure

Let’s make a pole figure in the [111] direction of the diamond structure, as seen in this figure from Wikipedia.

The figure caption reads as follows

The spots in the stereographic projection show the orientation of lattice planes with the 111 in the center. Only poles for a non-forbidden Bragg reflection are shown between Miller indices -10 <= (h,k,l) <= 10. The green spots contain Miller indices up to 3, for example 111, 113, 133, 200 etc in its fundamental order. Red are those raising to 5, ex. 115, 135, 335 etc, while blue are all remaining until 10, such as 119, 779, 10.10.00 etc.

diamond = Phase(space_group=227)
md = Miller.from_highest_indices(phase=diamond, uvw=[10, 10, 10])
md

Miller (9260,), point group m-3m, uvw
[[ 10.  10.  10.]
 [ 10.  10.   9.]
 [ 10.  10.   8.]
 ...
 [-10. -10.  -8.]
 [-10. -10.  -9.]
 [-10. -10. -10.]]

Remove duplicates under symmetry using Miller.unique()

md2 = md.unique(use_symmetry=True)
md2.size

285

Symmetrise to get all symmetrically equivalent directions

md3 = md2.symmetrise(unique=True)
md3

Miller (9260,), point group m-3m, uvw
[[  1.   0.   0.]
 [  0.   1.   0.]
 [ -1.   0.   0.]
 ...
 [ 10.  10. -10.]
 [-10.  10. -10.]
 [-10. -10. -10.]]

Remove forbidden reflections in face-centered cubic structures (all hkl must be all even or all odd)

selection = np.sum(np.mod(md3.hkl, 2), axis=1)
allowed = np.array([i not in [1, 2] for i in selection], dtype=bool)
md4 = md3[allowed]
md4

Miller (2330,), point group m-3m, uvw
[[  2.   0.   0.]
 [  0.   2.   0.]
 [ -2.   0.   0.]
 ...
 [ 10.  10. -10.]
 [-10.  10. -10.]
 [-10. -10. -10.]]

Assign colors to each class of vectors as per the description on Wikipedia

uvw = np.abs(md4.uvw)
green = np.all(uvw <= 3, axis=-1)
red = np.any(uvw > 3, axis=-1) * np.all(uvw <= 5, axis=-1)
blue = np.any(uvw > 5, axis=-1)
rgb_mask = np.column_stack([red, green, blue])

# Sanity check
print(np.count_nonzero(rgb_mask) == md4.size)

True

Rotate directions so that [111] impinges the unit sphere in the north pole (out of plane direction)

vz = Vector3d.zvector()
v111 = Vector3d([1, 1, 1])
r1 = Rotation.from_axes_angles(v111.cross(vz), v111.angle_with(vz).data)
r2 = Rotation.from_axes_angles(vz, -15, degrees=True)
md5 = r2 * r1 * md4

Restrict to upper hemisphere and remove duplicates

is_upper = md5.z > 0
md6 = md5[is_upper]
rgb_mask2 = rgb_mask[is_upper]

_, idx = md6.unit.unique(return_index=True)
md7 = md6[idx]
rgb_mask2 = rgb_mask2[idx]

Finally, plot the vectors

rgb = np.zeros_like(md7.uvw)
rgb[rgb_mask2] = 1

md7.scatter(c=rgb, s=10, grid=False, figure_kwargs=dict(figsize=(12, 12)))