Laboratory 09: Matrices and Rotations (Crystallography)

Background

The purpose of this laboratory is to understand and apply the concept of an agebraic group to produce one of the plane group patterns below. In this context, algebraic groups are mathematical objects that represent a collection of geometric transformations and are used in the study of crystallography. For example, when you say that is crystal is "face centered cubic" what you are really saying is that an atomic motif, when transformed by all the algebraic elements of a particular group, will be identical to the face centered cubic structure. The group is the mathematical machinary of crystal structures. The group itself is an abstract idea that is embodied in a crystal structure, a geometric design, a geometric series, a set of matrices, or a set of any kind where the rules that govern the set meet the criteria of a group. The idea of a group exists beyond the example I present here.

What Skills Will I Learn?

Using the context of two-dimensional plane groups, in this assignment you will practice the following skills:

  • Identifying a symmetry operation (rotation, translation, mirror) and determining the associated matrix representation.
  • Applying a symmetry operation to a position vector to transform the vector (or collection of vectors).
  • Seeing the inner workings of a class object and its attributes for the purpose of simplifying a computing task.
  • Using lists to collect and organize objects (such as transformation matrices)

What Steps Should I Take?

  1. Review the idea of "symmetry operations" and familarize yourself with mirror, translational, and rotational symmetry operations.
  2. Review the idea of matrix transformations and how they can be used to represent symmetry operations and how those symmetry operations can be represented as a matrix/vector dot product.
  3. Read about the metric tensor below and practice a few calculations of lengths and angles in unit cell coordinates.
  4. Read and practice writing down 2D, 3D and 3D augmented transformation matrices and learn how Euler angles are defined.
  5. For each case above, use a position vector pointing in a direction you select and apply a transformation using each matrix in turn. Does the resultant vector match your anticipated transformation? (e.g. If you take the position vector pointing at the [0,0,1] position in a cubic lattice and apply a 90 degree rotation about the x-axis, where should the resultant vector be pointing? Repeat this for mirror planes and translations.
  6. Review the concept of a "class" in object oriented programming and familarize yourself with the Point and Triangle classes I have defined below. Practice creating these objects and accessing their attributes.
  7. Review the code that accesses the Python Imaging Library (PIL) and look at the polygon function call and the Triangle points method. You can call the points method from within the polygon function to simplify the drawing of triangles.
  8. Using the Point and Triangle classes and your knowledge of symmetry transformations, reproduce one of the plane groups from the figure below. My suggestion is to start with a single triangle and apply a set of operations (transformation matrices) and collect new triangles into a list. Sketch a flowchart for how you might solve this problem.
  9. Review the very last block of code in this notebook and modify it to complete the assignment. You will need to use combinations of translations, reflections and rotations to accomplish this.
  10. Prepare a new notebook (not just modifications to this one) that describes your approach to reproducing one or more of the plane groups.

You may discover that a unique set of matrices will reproduce the whole pattern. This small set of matrices are an algebraic structure known as a group. So - the "plane group" is the group of symmetry operations that reproduces the structure in the plane starting from a single motif. The plane group representations are reproduced below for reference; this image comes from Hammond's book on crystallography, cited above.

A Sucessful Jupyter Notebook Will

  • Present a description of the essential elements of plane group symmetry, matrix algebra, and group theory to understand how these are related;
  • Identify the audience for which the work is intended;
  • Run the code necessary to draw one of the plane groups;
  • Provide a narrative and equations to explain why your approach is relevant to solving the problem;
  • Provide references and citations to any others' work you use to complete the assignment;
  • Be checked into your GitHub repository by the due date (one week from assignment).

A high quality communication provides an organized, logically progressing, blend of narrative, equations, and code that teaches the reader a particular topic or idea. You will be assessed on:

  • The functionality of the code (i.e. it should perform the task assigned).
  • The narrative you present. I should be able to read and learn from it. Choose your audience wisely.
  • The supporting equations and figures you choose to include.

If your notebook is just computer code your assignment will be marked incomplete.

Reading and Reference

  • M. De Graef and M. McHenry, Structure of Materials, Cambridge University Press, 2nd ed.
  • C. Hammond, The Basics of Crystallography and Diffraction, Oxford Science Publications, 4th ed.

The Plane Groups

Introduction


Operations using vector-like data structures are an essential component of numerical computing, mathematics, science and engineering. In the field of crystallography the use of vectors and rotations in real and reciprocal space helps to simplify the description of and quantitative operations on crystal structures. Matrix operations are used in the solution of partial differential equations. The vector algebra and rotations are most easily performed using matrix tools and representations. The concepts and their mathematical properties will be reviewed and demonstrated using symbolic algebra and numerical methods. A review of matrix algebra will be helpful.

Rotations


A vector can be transformed by translations, rotations, and stretching/shrinking. A matrix multiplication operation can be used to define each individual operation. We can use matrix multiplication to perform combinations of these operations and then this composite operator can be applied to a vector. In general these types of transformations are called Affine Transformations. A rotation in two dimensions is given by:

\begin{equation*} \left[\begin{matrix}\cos{\left (\theta \right )} & - \sin{\left (\theta \right )}\\\sin{\left (\theta \right )} & \cos{\left (\theta \right )}\end{matrix}\right] \end{equation*}

We can rotate a vector $x\mathbf{i} + y\mathbf{j}$ to the position $x'\mathbf{i} + y'\mathbf{j}$ using the following matrix operation:

\begin{equation*} \left( \begin{matrix} x' \\ y' \end{matrix} \right) = \left[\begin{matrix}\cos{\left (\theta \right )} & - \sin{\left (\theta \right )}\\\sin{\left (\theta \right )} & \cos{\left (\theta \right )}\end{matrix}\right] \cdot \left( \begin{matrix} x \\ y \end{matrix} \right) \end{equation*}

In [ ]:
%matplotlib inline

import numpy as np
import sympy as sp
sp.init_session(quiet=True)

In [ ]:
?sp.rot_axis3

In [ ]:
x = sp.symbols('\\theta')
sp.rot_axis3(x)

We can look up definitions, but we can also do some simple tests to see which way things rotate. Let us take a vector pointing in the $\mathbf{x}$ direction and rotate it about $\mathbf{z}$ by 90 degrees and see what happens:


In [ ]:
xUnit = sp.Matrix([1,0,0])
zRotation = sp.rot_axis3(sp.rad(90))

In [ ]:
xUnit

In [ ]:
# Each column can be viewed as where the unit vectors are moved to in the new space.
zRotation

In [ ]:
zRotation*xUnit

In [ ]:
# This should not work.
xUnit*zRotation

In [ ]:
xUnit.T*zRotation

What can we learn from this result? It is now known that:

  • The convention for positive angles is a counterclockwise rotation.
  • The rotation axis function in sympy as defined rotates clockwise
  • There are conventions about active and passive rotations.
  • Don't assume module functions will do what you want - always check.

Rather than rely on module functions, we can define our own rotation function. Using a function called "isclose" or Boolean indexing it is possible to clean up the arrays and remove small numerical values that should be zeros.


In [ ]:
def rotation2D(theta):
    return np.array([[np.cos(theta), np.sin(theta)],
                     [-np.sin(theta),  np.cos(theta)]])

In [ ]:
testArray = rotation2D(0)
testArray

In [ ]:
np.dot(np.array([1,0]),testArray)

DIY: Computing Rotation Matrices


Compute the following rotation matrices:

  • A rotation of 0$^\circ$ about the origin.
  • A rotation of 45$^\circ$ about the origin.
  • A rotation of 60$^\circ$ about the origin.
  • A rotation of 90$^\circ$ about the origin.
  • A rotation of 180$^\circ$ about the origin.
  • A rotation of -270$^\circ$ about the origin.

In [ ]:
# Put your code here.

Cleaning up the Small Values


Sometimes Numpy returns very small numbers instead of zeros.

One strategy is to remove small numbers less than some tolerance and set them equal to zero. Algorithms like these where you compare your data to a tolerance and then operate on the entries that meet certain criteria are not uncommon in numerical methods. This is the tradeoff between symbolic computation and numeric computation.


In [ ]:
testArray[np.abs(testArray) < 1e-5] = 0

testArray

The key is in the Boolean comparision using the < symbol. The expression returns a numpy array of dtype=bool. Let me say here that it is good to check the results of expressions if you are unsure.


In [ ]:
np.abs(testArray) < 1e-5

We can write a function to do this that is a bit more robust. Modifications are done in-place (by reference) so we just return the array passed to the function after some manipulation that we do by Boolean indexing.


In [ ]:
def zeroArray(testArray):
    testArray[np.isclose(testArray, np.zeros_like(testArray))] = 0.0
    return testArray

In [ ]:
modifiedArray = rotation2D(np.pi/2)
modifiedArray

In [ ]:
modifiedArray = zeroArray(rotation2D(np.pi/2))
modifiedArray

Rotations (Continued)


Using the new rotation2D function and zeroArray function we can now write:


In [ ]:
zeroArray(np.dot(np.array([1,0]),rotation2D(np.pi/2)))

A collection of functions for performing transformations is available at http://www.lfd.uci.edu/~gohlke/. This can be imported and the functions used in your code. The fastest way to explore what is available is to import the file and then use autocomplete and docstring viewing functions from the Jupyter notebook.


In [ ]:
import transformations as tfm

In [ ]:
tfm.rotation_matrix(np.pi/2, [0,1,0])

In [ ]:
zeroArray(tfm.rotation_matrix(np.pi/2, [0,1,0]))

Symmetry Operations and Translations in Crystals


A generalized affine transformation in two dimensions can be thought of as an augmented matrix as:

$$\begin{bmatrix} r_1 & r_2 & t_x\\ r_3 & r_4 & t_y\\ 0 & 0 & 1\\ \end{bmatrix}$$

so you could imagine the following:

$$\begin{bmatrix} x'\\ y'\\ 1\\ \end{bmatrix} = \begin{bmatrix} 1 & 0 & t_x\\ 0 & 1 & t_y\\ 0 & 0 & 1\\ \end{bmatrix} \begin{bmatrix}x\\ y\\ 1\\ \end{bmatrix} $$

expanding to:

$$x' = x + t_x $$

and

$$y' = y + t_y $$

We can explicitly write the rotation components as listed earlier:

$$\begin{bmatrix} x'\\ y'\\ 1\\ \end{bmatrix} = \begin{bmatrix} \cos{\theta} & \sin{\theta} & t_x\\ -\sin{\theta} & \cos{\theta} & t_y\\ 0 & 0 & 1\\ \end{bmatrix} \begin{bmatrix}x\\ y\\ 1\\ \end{bmatrix} $$

where the $r_i$ represent the rotation matrix components and the $t_i$ represent the translations components. This can be thought of as a shearing operation in 3D. The Wikipedia article on this topic expands this idea a bit more.

In this format we can use a point description that looks like $(x, y, t)$ and matrix algebra to generate our transformed points. Using SymPy:


In [ ]:
alpha, t_x, t_y, x, y = sp.symbols('alpha t_x t_y x y')

In [ ]:
sa = sp.sin(alpha)
ca = sp.cos(alpha)
M = sp.Matrix([[ca, sa, t_x], [-sa, ca, t_y], [0, 0, 1]])
V = sp.Matrix([x, y, 1])

M*V

Helper Classes, Drawing and Plane Groups


Let us explore a bit of how we can draw an image - and then we have everything we need to start building the plane group representations. To build pictoral representations of plane groups, two classes have been created to simplify the organization of the motif.

Below are the Point and Triangle classes for your use. A Point has storage for an $(x,y)$ position. Rather than building arrays and referencing specific positions within the array a more natural referencing of points is possible. If we define a point p1 = Point(10,20) we can access the points by p1.x and p1.y. The code is more easily read and debugged with this syntax.

Building on the Point class we define a Triangle class. The Triangle permits access of each point and defines an affine() method that will take a Numpy array that represents a transformation matrix. This method returns a new instance of a Triangle and preserves the original points. Writing the code this way avoids explicit handling of the transformation matrices.

These two classes and methods are demonstrated in the second and third code blocks. Building on this demonstration you will be able to complete the homework.


In [ ]:
# Class definitions

from math import sqrt
import numpy as np

class Point:
    """
    A Point object to simplify storage of (x,y) positions.
    p1.x, p1.y, etc...
    """
    def __init__(self,x_init,y_init):
        self.x = x_init
        self.y = y_init

    def shift(self, x, y):
        self.x += x
        self.y += y

    def __repr__(self):
        return "".join(["Point(", str(self.x), ",", str(self.y), ")"])
    
class Triangle:
    """
    A Triangle class constructed from points.  Helps organize information on
    triangles.  Has a points() method for returning points in a form that 
    can be used with polygon drawing from Python Image Library (PIL) and a 
    method, affine(), that applies a matrix transformation to the points.
    """
    def __init__(self, p1_init, p2_init, p3_init):
        self.p1 = p1_init
        self.p2 = p2_init
        self.p3 = p3_init
        
    def points(self):
        x1, y1 = self.p1.x, self.p1.y
        x2, y2 = self.p2.x, self.p2.y
        x3, y3 = self.p3.x, self.p3.y
        
        return [(x1,y1),(x2,y2),(x3,y3)]
        
    def affine(self, affineMatrix):
        """
        Applies an affine transformation to a triangle and changes the 
        points of the triangle.  This code returns a new triangle.  Uses
        Points to simplify augmenting the matrix and dot products.
        """
        x1, y1 = self.p1.x, self.p1.y
        x2, y2 = self.p2.x, self.p2.y
        x3, y3 = self.p3.x, self.p3.y
        
        p1Vector = np.array([[x1, y1, 1]])
        p2Vector = np.array([[x2, y2 , 1]])
        p3Vector = np.array([[x3, y3, 1]])
        
        p1New = np.dot(affineMatrix, p1Vector.T)
        p2New = np.dot(affineMatrix, p2Vector.T)
        p3New = np.dot(affineMatrix, p3Vector.T)
        
        # This line needs to be cleaned up.
        newTriangle = Triangle(Point(p1New[0,0],p1New[1,0]),Point(p2New[0,0],p2New[1,0]),Point(p3New[0,0],p3New[1,0]))
        
        return newTriangle

Using the Python Imaging Library


In order that our transformations can be visualized, we will use the Python Imaging Library and the polygon() method. The polygon() method takes a list of points in the format returned by our Triangle class method points(). The code below is a very simple starting point for the student to begin building a representation of the plane groups.


In [ ]:
# Class demonstrations

%matplotlib inline

import numpy as np
import math
from PIL import Image, ImageDraw

image = Image.new('RGB', (500, 500), 'white')
draw = ImageDraw.Draw(image)

t1 = Triangle(Point(10,10),Point(40,10),Point(10,50))
t2 = Triangle(Point(10,10),Point(40,10),Point(10,50))
translationMatrix = np.array([[1,0,50],[0,1,0],[0,0,1]])
reflectionMatrix = np.array([[-1,0,0],[0,1,0],[0,0,1]])

t2 = t1.affine(reflectionMatrix*translationMatrix)

print(t1.points(), t2.points())

draw.polygon(t1.points(), outline='black', fill='red')
draw.polygon(t2.points(), outline='black', fill='green')

image

In [ ]:
# A slightly more advanced demonstration

%matplotlib inline

import numpy as np
import math
from PIL import Image, ImageDraw

image = Image.new('RGB', (500, 500), 'white')
draw = ImageDraw.Draw(image)

t1 = Triangle(Point(70,10),Point(100,10),Point(70,50))

triangleList = [t1]

translationMatrix = np.array([[1,0,10],[0,1,0],[0,0,1]])
reflectionMatrixY = np.array([[-1,0,0],[0,1,0],[0,0,1]])
reflectionMatrixX = np.array([[1,0,0],[0,-1,0],[0,0,1]])

r90 = np.array([[0,-1,0],[1,0,0],[0,0,1]])

triangleList.append(t1.affine(r90))
triangleList.append((t1.affine(r90)).affine(r90))
triangleList.append((t1.affine(r90)).affine(r90).affine(r90))

tempList = [triangle.affine(reflectionMatrixX) for triangle in triangleList]
triangleList.extend(tempList)

# Using an affine transformation to center the triangles in the drawing
# as canvas coordinates are (0,0) at the top left.
centerMatrix = np.array([[1,0,250],[0,1,250],[0,0,1]])
drawList = [triangle.affine(centerMatrix) for triangle in triangleList]
for triangle in drawList:
    draw.polygon(triangle.points(), outline='black', fill='red')

image

Advanced Topics: The Metric Tensor

Students who learn crystallography for the first time are introduced to the topic through study of cubic crystal structures. This permits an appeal to our intuition about orthonormal (Euclidian) coordinate systems. This approach misses out on a more general method for teaching the topic where the d-spacing and the angle between directions can be worked out for any general crystal system. The method for describing distances in a general reference frame involves the metric tensor. The metric tensor defines how distances are measured in every direction and its components are the dot product of every combination of basis vector in the system of interest. We use the standard lattice parameter designations:

$$ [a, b, c, \alpha, \beta, \gamma] $$

where $\gamma$ is the angle between $\mathbf{a}$ and $\mathbf{b}$, etc. This general system has basis vectors:

$$ \mathbf{a},\mathbf{b}, \mathbf{c} $$

The standard geometric interpretation of an inner product of vectors is:

$$ \mathbf{a} \cdot \mathbf{b} = |a| \; |b| \; \cos{\gamma} $$

so that the metric tensor is:

$$ g_{ij} = \begin{bmatrix} \mathbf{a} \cdot \mathbf{a} & \mathbf{a} \cdot \mathbf{b} & \mathbf{a} \cdot \mathbf{c} \\ \mathbf{b} \cdot \mathbf{a} & \mathbf{b} \cdot \mathbf{b} & \mathbf{b} \cdot \mathbf{c} \\ \mathbf{c} \cdot \mathbf{a} & \mathbf{c} \cdot \mathbf{b} & \mathbf{c} \cdot \mathbf{c} \end{bmatrix} $$

Committing this to memory is simple and it has an intuitive meaning in that each entry measures a different projection of a vector component onto an axis in a general crystal system within lattice coordinates. It is possible to write a small function in SymPy to illustrate this.


In [ ]:
def metricTensor(a=1, b=1, c=1, alpha=sp.pi/2, beta=sp.pi/2, gamma=sp.pi/2):
    return sp.Matrix([[a*a, a*b*sp.cos(gamma), a*c*sp.cos(beta)], \
                     [a*b*sp.cos(gamma), b*b, b*c*sp.cos(alpha)], \
                     [a*c*sp.cos(beta), b*c*sp.cos(alpha), c*c]])

In [ ]:
sp.var('a b c alpha beta gamma u v w h k l')

In [ ]:
M = metricTensor(a=a,
                 b=b,
                 c=c,
                 alpha=alpha,
                 beta=beta,
                 gamma=gamma
                )
M

Two Common Computations


Using the Einstein summation convention, the square of the distance between two points located at the end of vectors $\mathbf{p}$ and $\mathbf{q}$ is given by:

$$ D^2 = (\mathbf{q} - \mathbf{p})_i \; g_{ij} \; (\mathbf{q} - \mathbf{p})_j $$

The dot product between two vectors is given by:

$$ \mathbf{p} \cdot \mathbf{q} = p_i \; g_{ij} \; q_j $$

Note that the vectors $\mathbf{p}$ and $\mathbf{q}$ are in lattice coordinates.

DIY: Angle Between Two Vectors


Find the expression for and compute the angle between two vectors in a general coordinate system. Use the function defined above or write a new one. You are encouraged to use SymPy or Numpy in your calculations. Refer to the earlier lectures that cover these topics for the implementation and technical details.


In [ ]:
# Put your code or markdown here.

DIY: Compute the angle between the $(123)$ plane and the $(112)$ plane in a trigonal crystal system


The trigonal crystal system has the least symmetry. Refer to standard texts for the pattern of lattice parameters.


In [ ]:
vectorOne = np.array([1,2,3])
vectorTwo = np.array([1,1,2])

# Put your code here.

Advanced Topics: Euler Angles

This discussion is derived primarily from Arfken, Chapter 3.3. The figure below is from Arfken:

There are three successive rotations used in the Euler angle formalism - the product of these three rotations is the single operation that transforms one set of coordinates $(x,y,z)$ into another, $(x',y',z')$. The order is important and is the difference between active and passive rotations.

In steps as shown in Figure 3.7 from Arfken:

  1. The first rotation is about $x_3$. In this case the $x'_3$ and $x_3$ axes coincide. The angle $\alpha$ is taken to be positive (counterclockwise). Our new coordinate system is $(x'_1,x'_2,x'_3)$.
  2. The coordinates are now rotated through an angle $\beta$ around the $x'_2$ axis. Our new coordinate system is now $(x''_1,x''_2,x''_3)$.
  3. The final rotation is through the angle $\gamma$ about the $x'''_3$ axis. Our coordinate system is now the $(x'''_1,x'''_2,x'''_3)$. In the case pictured above the $x''_3$ and $x'''_3$ axes coincide.

For example:


In [ ]:
alpha, beta, gamma = sp.symbols('alpha beta gamma')

def rZ(angle):
    sa = sp.sin(angle)
    ca = sp.cos(angle)
    M = sp.Matrix([[ca, sa, 0],
                  [-sa, ca, 0],
                  [0, 0, 1]])
    return M

def rY(angle):
    sb = sp.sin(angle)
    cb = sp.cos(angle)
    M = sp.Matrix([[cb, 0, -sb],
                  [0, 1, 0],
                  [sb, 0, cb]])
    return M

In [ ]:
(rZ(alpha), rY(beta), rZ(gamma))

You'll find that the symbolic triple matrix product matches up with the results in Arfken for the definition of Euler angles $\alpha$, $\beta$, $\gamma$. Note also that this is much easier to compute than by hand! Also - less likely to result in errors.


In [ ]:
rZ(gamma)*rY(beta)*rZ(alpha)

To convert this symbolic expression to a numerical expression, one option is to use lambdafy from SymPy:


In [ ]:
eulerAngles = sp.lambdify((alpha,beta,gamma), rZ(gamma)*rY(beta)*rZ(alpha), "numpy")

In [ ]:
np.array([1,0,0]).dot(eulerAngles(np.pi/2.0,0,0))