Finding Sample Tilt for Desired Zone axis in TEM

[anderscmathisen]

Hi,

This is not really an issue, but something I have been working on and thought would be productive to share. Did not know if this is too specific to be included in Orix in any way, but might be useful as a tutorial to highlight how orientations may be combined.

The main idea here is to find the x- and y-tilt of the double tilt gonio system in a TEM, which orients the crystal to be in a chosen zone axis. To solve this problem, the orientation, $g_0$, of the crystal in a reference sample tilt $(x_0, y_0)$ needs to be known (e.g from template matching using Pyxem on SPED data), and the position of the gonio tilt axes with respect to the sample reference frame needs to be known. The position of the gonio tilt axes can be found in a tilt series.

The rotation matrix which describes the rotation of the double tilt holder in a TEM, when the x-tilt is changed from $\alpha_0$ to $\alpha_0 + \alpha$, and the y-tilt is changed from $\beta_0$ to $\beta_0 + \beta$ is given by the product $T_1 * T_2$ [cite], where

$$ \begin{align} T_1(\alpha) & = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos{\alpha} & -\sin{\alpha}\\ 0 & \sin{\alpha} & \cos{\alpha} \end{bmatrix} \\ T_2(\beta; \alpha_0) & = \begin{bmatrix} \cos{\beta} & -\sin{\beta}\sin{\alpha_0} & \sin{\beta}\cos{\alpha_0} \\ \sin{\beta}\sin{\alpha_0} & \cos^2\alpha_0 + \sin^2\alpha_0\cos{\beta} & \sin{\alpha_0}\cos{\alpha_0}(1-\cos{\beta})\\ -\sin{\beta}\cos{\alpha_0} & \sin{\alpha_0}\cos{\alpha_0}(1-\cos{\beta}) & \sin^2\alpha_0 + \cos^2\alpha_0\cos{\beta} \end{bmatrix} \end{align} $$

The product $T_1 * T_2$ I'll refer to as $T(\alpha,\beta;\alpha_0, \beta_0)$.

$T(\alpha,\beta;\alpha_0, \beta_0)$ of course operates in the reference frame of the gonio tilt axes, while the orientation $g_0$ is defined with respect to the crystal and sample reference frame as lab-to-crystal, where the sample reference frame usually is defined by the directions if the scan in SPED, and will not necessarily coincides with the gonio reference frame. $T(\alpha,\beta;\alpha_0, \beta_0)$ can therefore not simply be multiplied with $g_0$ to find the transformed orientation. The gonio and sample reference frame do, however, share the optical axis of the TEM as the z-axis, and so the two are related by a rotation about the z-axis, which I'll call $R(\theta)$.

Thus, to meaningfully transform the orientation $g_0$ with $T(\alpha,\beta;\alpha_0, \beta_0)$ to find the new orientation $g$ when the sample tilt is $(x, y)$, the transformation used is

$$ g = (R(-\theta)T(x-x_0, y-y_0; x_0, y_0)R(\theta)g_0^{-1})^{-1} $$

where the first inversion of $g_0$ converts to crystal-to-lab definition of orientations (describing the sample reference frame), the multiplication with $R(\theta)$ converts the sample reference frame to the gonio reference frame, where the transformation $T(x-x_0, y-y_0; x_0, y_0)$ is applied, after which $R(-\theta)$ converts back to the sample reference frame, and the final inversion converts back to the standard lab-to-crystal definition of orientations.

To use the transformation above to find what sample tilt leaves the crystal in a desired zone axis, one can simply search for values of $x$ and $y$ which minimizes the angle between the desired zone axis in the new, transformed orientation (or any symmetrically equivalent zone axis) and the optical axis of the TEM. This minimization can be done using SciPys Minimize functions.

Below, I have included some code implementing the procedure. From my testing, the approach appears to work well, and the predicted x- and y-tilt is usually within 1° of the actual zone axis.

import numpy as np
from orix.quaternion import Rotation, Orientation
from orix.vector import Miller
from scipy.optimize import minimize

class GonioPosition:
    """This class represents a sample tilt position"""
    def __init__(self, x_tilt, y_tilt, degrees=True) -> None:
        if degrees:
            self.x_tilt = np.radians(x_tilt)
            self.y_tilt = np.radians(y_tilt)
        else:
            self.x_tilt = x_tilt
            self.y_tilt = y_tilt


    def __repr__(self) -> str:
        return f"x_tilt: {np.rad2deg(self.x_tilt) :.2f}, y_tilt: {np.rad2deg(self.y_tilt) :.2f}"



class RotationGenerator:
    """This class is used to generate the rotation matrix T(alpha,beta; alpha_0, beta_0"""
    def __init__(
        self,
        new_gonio_pos: GonioPosition,
        old_gonio_pos: GonioPosition = GonioPosition(0, 0),
    ) -> None:
        self.new_gonio_pos = new_gonio_pos
        self.old_gonio_pos = old_gonio_pos

        self.alpha_0 = self.old_gonio_pos.x_tilt
        self.beta_0 = self.old_gonio_pos.y_tilt

        self.alpha = self.new_gonio_pos.x_tilt - self.old_gonio_pos.x_tilt
        self.beta = self.new_gonio_pos.y_tilt - self.old_gonio_pos.y_tilt

    def get_T1(self):
        a11 = 1
        a12 = 0
        a13 = 0
        a21 = 0
        a22 = np.cos(self.alpha)
        a23 = -np.sin(self.alpha)
        a31 = 0
        a32 = np.sin(self.alpha)
        a33 = np.cos(self.alpha)

        return np.array([[a11, a12, a13], 
                         [a21, a22, a23], 
                         [a31, a32, a33]])

    def get_T2(self):
        a11 = np.cos(self.beta)
        a12 = -np.sin(self.beta) * np.sin(self.alpha_0)
        a13 = np.sin(self.beta) * np.cos(self.alpha_0)
        a21 = np.sin(self.beta) * np.sin(self.alpha_0)
        a22 = np.cos(self.alpha_0) ** 2 + np.sin(self.alpha_0) ** 2 * np.cos(self.beta)
        a23 = np.sin(self.alpha_0) * np.cos(self.alpha_0) * (1 - np.cos(self.beta))
        a31 = -np.sin(self.beta) * np.cos(self.alpha_0)
        a32 = np.sin(self.alpha_0) * np.cos(self.alpha_0) * (1 - np.cos(self.beta))
        a33 = np.sin(self.alpha_0) ** 2 + np.cos(self.alpha_0) ** 2 * np.cos(self.beta)

        return np.array([[a11, a12, a13], 
                         [a21, a22, a23], 
                         [a31, a32, a33]])

    def get_full_rotation_matrix(self):
        return self.get_T1() @ self.get_T2()


def rotate_inplane(ori, angle ):
    """This function is used to get the rotation R(theta) converting 
    between the sample and gonio reference frames."""

    rotatior = Rotation.from_axes_angles([0,0,-1], np.deg2rad(angle))

    return Orientation(rotatior * ori, symmetry=ori.symmetry)


def transform_to_gonio(ori, new_gonio, old_gonio, tilt_axes_align_angle = 9.61):
    """This function describes the full transformation of a crystal orientation, as described above
         The tilt_axes_align_angle parameter is the theta parameter in R(theta)."""
    rotgen = RotationGenerator(new_gonio, old_gonio)
    rotator = Rotation.from_matrix(rotgen.get_full_rotation_matrix())
    
    ori = ~ori
    ori = rotate_inplane(ori, tilt_axes_align_angle)
    ori = Orientation(rotator * ori, symmetry = ori.symmetry)
    ori = rotate_inplane(ori, -tilt_axes_align_angle)
    ori = ~ori

    return ori

And applying transform_to_gonio() to find a desired zone axis by Scipy Minimization

def solve_Gonio_position(current_orientation : Orientation,
                         current_goniopos : GonioPosition, 
                         desired_zoneaxis : Miller):
    
    equivalent_zoneaxes = desired_zoneaxis.symmetrise()
    def smallest_angle_for_tilt(gonio : np.ndarray, current_orientation : Orientation):
        
        new_gonio = GonioPosition(*gonio)
        
        new_orientation = transform_to_gonio(current_orientation, new_gonio, current_goniopos)
        
        angles = np.rad2deg((~new_orientation * equivalent_zoneaxes).angle_with(Miller(uvw=[0,0,1], phase = desired_zoneaxis.phase)))

            return np.min(angles)
    
    #Include bounds for the minimization to only give results within the max tilt-range of the sample holder
    res = minimize(smallest_angle_for_tilt, np.array([0,0]), args=current_orientation, bounds=((-30,30), (-30,30)), method = "Nelder-Mead")

    x_opt =res.x

    return GonioPosition(*x_opt)

The application of the transformation described above to find the x- and y-tilt of specific zone axes is very useful on its own (especially for polycrystalline sample with small grains), and has saved me many hours on the TEM. I do think, however, that the transformation could be used for many other problems. For instance, if one is studying slanted grain boundaries in TEM, and somehow manage to describe the normal vector of the grain boundary, $\vec{n}$, then it should be possible to use the transformation described above in some way to find the set of all x- and y- tilts (probably a line through the space of x- and y-tilts) which keep the grain boundary edge on (e.g. $\vec{n}$ perpendicular to optical axis), while rotating the sample around $\vec{n}$. Allowing the grain boundary to be studied edge on, from multiple different perspectives.

I would be happy to assist with implementing some tutorial for the docs if that is desired, but don't really know where to start with that (code would definitively need a redo, above is only slapped together without any real concern for readability). Also please let me know if you want me to provide any data to show the application of the functions above, or in general have any questions or remarks.

[hakonanes]

Hi @anderscmathisen,

thank you for sharing your procedure with us here. As you say, it can be a real time saver, so it is definitely worth creating a tutorial showing an example of how to use it. If you think such a tutorial benefits from showing diffraction patterns, it might be best to add it to pyxem's documentation. If you think plotting unit cells and pole figures is enough, then it should be a perfect fit for orix. Perhaps a simple helper function to show a unit cell embedded within a tilted sample under an incoming electron beam could help understand what's going on as well...

Question: Do you combine your set of rotations $g$ from right-to-left (first-to-last)? orix uses passive rotations, and according to the Rowenhorst et al. (2015) paper, they are combined from left-to-right (first-to-last), like shown in a new example in the dev docs.

[anderscmathisen]

Hmm yes I suppose I have accidentally used the right-to-left convention in combining the rotations.

In the passive rotation left-to-right convention, I guess the transformation would be

$$ g = (g_0^{-1}R(\theta)T(x-x_0, y-y_0; x_0, y_0)R(-\theta))^{-1}? $$

Feel free to correct me if I still have missunderstood the passive nature of rotations.

[...] If you think such a tutorial benefits from showing diffraction patterns, it might be best to add it to pyxem's documentation. If you think plotting unit cells and pole figures is enough, then it should be a perfect fit for orix. Perhaps a simple helper function to show a unit cell embedded within a tilted sample under an incoming electron beam could help understand what's going on as well...

I don't necessarily think the tutorial would need to show diffraction patterns. I think the core of the idea here is rooted in orientations and rotations, and how to combine them. So I believe this would be a better fit for the docs of Orix. A helper function showing unit cells would really be great here.