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.
Using the context of two-dimensional plane groups, in this assignment you will practice the following skills:
Point
and Triangle
classes I have defined below. Practice creating these objects and accessing their attributes.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.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.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 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:
If your notebook is just computer code your assignment will be marked incomplete.
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.
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:
sympy
as defined rotates clockwiseRather 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)
Compute the following rotation matrices:
In [ ]:
# Put your code here.
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
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]))
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
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
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
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
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.
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.
In [ ]:
vectorOne = np.array([1,2,3])
vectorTwo = np.array([1,1,2])
# Put your code here.
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:
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))