Wyckoff Positions
Contents
Wyckoff Positions¶
The equivalent sites in a structure are linked to each other via the group’s symmetry operators. If \(B\) site can be reached by applying a symmetry element \(g_i\) of the group \(\mathcal{G}\) to a site \(A\), such that:
then site \(A\) is equivalent to site \(B\) and they belong to the same orbit. The number of sites in an orbit gives the multiplicity of the said Wyckoff Position. Since all the sites in an orbit are equivalent to each other and all the rest can be derived from any one in the orbit, each site in the orbit can act as the representative of the orbit.
Crystals¶
Thanks to the periodicity of the crystals, it is ensured that each one of the unit cells building up the crystal is identical to each other. In other words any translation along a crystal axis equal to an integer number of the corresponding crystal lattice vector in that direction will bring us to an equivalent site. This is indicated as the de facto crystal translation operators identified by:
As their inverses and products (combinations) are also symmetry operators according to the group theory axioms, they span whole space.
Using this fact, for every site lying outside our unit cell of interest, there is an equivalent site in our unit cell. For conveniency, we work on the unit cell confined within \(\{x_i\}\in[0,1)\). Thus, for example, a site with coordinates \(x'=(-0.2,1.3,0.7)\) is equivalent to a site within the unit cell with coordinates \(x=(0.8,0.3,0.7)\).
We can either add/subtract 1 from the coordinates until arriving in the reference unit cell or more practically, subtract the “floor” of the positions from the site’s coordinates, i.e.,
where the \(\lfloor \dots \rfloor\) indicates the floor function.
import numpy as np
x_p = np.array([-0.2,1.3,0.7])
x = x_p - np.floor(x_p)
print(x)
[0.8 0.3 0.7]
Note
Any crystal lattice in 3D can be classified into one of the 14 Bravais lattice types. Here, the second letter of the Bravais lattice symbol (the capital one) indicates the centering with the corresponding additional crystal translation operators:
Primitive (P): no additional translation operator
Base-centered (S): can be one of the A,B,C axes, where:
A: \(t(0,0.5,0.5)\)
B: \(t(0.5,0,0.5)\)
C: \(t(0.5,0.5,0)\)
Body-centered (I): \(t(0.5,0.5,0.5)\)
Face-centered (F): \(t(0,0.5,0.5)\),\(t(0.5,0,0.5)\),\(t(0.5,0.5,0)\)
Obtaining the orbits¶
An orbit containing the equivalent sites of a given site is constructed by applying every symmetry operator of the group. Some of these sites may coincide and at most, the number of distinct sites will be equal to the number of symmetry operators of the group (meaning that, every operator applied will correspond to a different site in the reference unit cell). This kind of sites where the number of orbit elements is equal to the number of symmetry operators of the group are called “general Wyckoff positions”. A Wyckoff position is denoted by a symbol containing a number and a letter, e.g., 4d. The number denotes the number of elements in the corresponding orbit and called the “multiplicity” of the site. The letter is reserved for further classification which will be explained shortly.
Example: Construction of orbits¶
For the space group P4 (#75), construct the orbits for the following sites:
\(x_1 = (0.1,0.7,0.4)\)
\(x_2 = (0,0.5,0.3)\)
\(x_3 = (0.5,0.5,0.4)\)
\(x_4 = (0,0,0.9)\)
The symmetry operators for the P4 (#75) space group can be obtained from the GENPOS tool:
Defining them in the augmented matrix form, we have:
o1 = np.eye(4,4)
o2 = np.diag([-1,-1,1,1])
o4p = np.array([[0,-1,0,0],[1,0,0,0],[0,0,1,0],[0,0,0,1]])
o4m = np.dot(o4p,o2)
Note
Notice how we defined the \(4^-\) operator from the multiplication of \(4^+\) and \(2\). We could actually derive all the operators by first defining \(4^+\) (or \(4^-\) for that matter), and by taking its powers, although this is not a generally valid method – the groups that have this property are called cyclic groups.
o4m
array([[ 0, 1, 0, 0],
[-1, 0, 0, 0],
[ 0, 0, 1, 0],
[ 0, 0, 0, 1]])
Collecting them in a list and iterating over them, we construct the orbits:
x1 = np.array([[0.1],[0.7],[0.4],[1]])
x2 = np.array([[0],[0.5],[0.3],[1]])
x3 = np.array([[0.5],[0.5],[0.4],[1]])
x4 = np.array([[0],[0],[0.9],[1]])
sites = [x1,x2,x3,x4]
operators = [o1,o4p,o2,o4m]
op_names = ["o1","o4p","o2","o4m"]
orbits = []
for site in sites:
orbit = []
print("Calculating the orbit of:\n ({:5.2f},{:5.2f},{:5.2f})"\
.format(*site.flatten()))
for i in range(len(operators)):
operator = operators[i]
op_site = np.dot(operator, site)
print("Applying {:>4s} yields: ({:5.2f},{:5.2f},{:5.2f})"\
.format(op_names[i],*op_site.flatten()),end="")
if(np.any(op_site<0) or np.any(op_site>=1)):
op_site = op_site - np.floor(op_site)
print(" -> ({:.2f},{:.2f},{:.2f})"\
.format(*op_site.flatten()),end="")
orbit.append(op_site)
print("")
# Remove the duplicate sites in the orbit
orbit = np.unique(orbit,axis=0)
orbits.append(orbit)
print("-"*30)
Calculating the orbit of:
( 0.10, 0.70, 0.40)
Applying o1 yields: ( 0.10, 0.70, 0.40) -> (0.10,0.70,0.40)
Applying o4p yields: (-0.70, 0.10, 0.40) -> (0.30,0.10,0.40)
Applying o2 yields: (-0.10,-0.70, 0.40) -> (0.90,0.30,0.40)
Applying o4m yields: ( 0.70,-0.10, 0.40) -> (0.70,0.90,0.40)
------------------------------
Calculating the orbit of:
( 0.00, 0.50, 0.30)
Applying o1 yields: ( 0.00, 0.50, 0.30) -> (0.00,0.50,0.30)
Applying o4p yields: (-0.50, 0.00, 0.30) -> (0.50,0.00,0.30)
Applying o2 yields: ( 0.00,-0.50, 0.30) -> (0.00,0.50,0.30)
Applying o4m yields: ( 0.50, 0.00, 0.30) -> (0.50,0.00,0.30)
------------------------------
Calculating the orbit of:
( 0.50, 0.50, 0.40)
Applying o1 yields: ( 0.50, 0.50, 0.40) -> (0.50,0.50,0.40)
Applying o4p yields: (-0.50, 0.50, 0.40) -> (0.50,0.50,0.40)
Applying o2 yields: (-0.50,-0.50, 0.40) -> (0.50,0.50,0.40)
Applying o4m yields: ( 0.50,-0.50, 0.40) -> (0.50,0.50,0.40)
------------------------------
Calculating the orbit of:
( 0.00, 0.00, 0.90)
Applying o1 yields: ( 0.00, 0.00, 0.90) -> (0.00,0.00,0.90)
Applying o4p yields: ( 0.00, 0.00, 0.90) -> (0.00,0.00,0.90)
Applying o2 yields: ( 0.00, 0.00, 0.90) -> (0.00,0.00,0.90)
Applying o4m yields: ( 0.00, 0.00, 0.90) -> (0.00,0.00,0.90)
------------------------------
for i in range(len(orbits)):
print("Orbit of x{:} ({:.2f},{:.2f},{:.2f}):"\
.format(i+1,*sites[i].flatten()))
for j in range(len(orbits[i])):
orb_element = orbits[i][j]
print(" {:d}. ({:.2f},{:.2f},{:.2f})"\
.format(j+1,*orb_element.flatten()))
print("Hence, the multiplicity of this Wyckoff position is: {:d}"\
.format(j+1))
print("-"*54)
Orbit of x1 (0.10,0.70,0.40):
1. (0.10,0.70,0.40)
2. (0.30,0.10,0.40)
3. (0.70,0.90,0.40)
4. (0.90,0.30,0.40)
Hence, the multiplicity of this Wyckoff position is: 4
------------------------------------------------------
Orbit of x2 (0.00,0.50,0.30):
1. (0.00,0.50,0.30)
2. (0.50,0.00,0.30)
Hence, the multiplicity of this Wyckoff position is: 2
------------------------------------------------------
Orbit of x3 (0.50,0.50,0.40):
1. (0.50,0.50,0.40)
Hence, the multiplicity of this Wyckoff position is: 1
------------------------------------------------------
Orbit of x4 (0.00,0.00,0.90):
1. (0.00,0.00,0.90)
Hence, the multiplicity of this Wyckoff position is: 1
------------------------------------------------------
Compare these with the Wyckoff positions list obtained from the WYCKPOS tool:
One can also directly use the WPASSIGN tool to obtain the orbit elements. For our case, entering the corresponding space group and site information (lattice parameters are not used for these calculations, so they can be arbitrary, as well as the atomic symbols):
# Space Group ITA number
75
# Lattice parameters
5.0 5.0 7.0 90 90 90
# Number of independent atoms in the asymmetric unit
3
# [atom type] [number] [WP] [x] [y] [z]
X 1 - 0.10 0.70 0.40
X 2 - 0.00 0.50 0.30
X 3 - 0.50 0.50 0.40
X 4 - 0.00 0.00 0.90
Upon clicking the “Show” button, the following result page appears:
Checking the listing of WYCKPOS we encountered above, we see that the ‘coordinates’ columns list the form for fix and free parameters for each Wyckoff position. For our case, we see that, corresponding to the ‘4d’ position, we have:
(x,y,z) (-x,-y,z) (-y,x,z) (y,-x,z)
so, for the \(x_1 = (0.1,0.7,0.4)\), where \(x=0.1\), \(y=0.7\) and \(z=0.4\), the orbit elements are constructed accordingly:
A “connection” with the symmetry operators¶
Consider the \((-x,-y,z)\) entry: let’s try to construct a \(3\times3\) “matrix” from this by interpreting each component as a row. If we have an ‘x’, let’s interpret it as ‘[1 0 0]’; ‘-y’ corresponds to ‘[0 -1 0]’; a component like ‘y-z’ is interpreted as [0 1 -1], etc. (More information on the xyz notation can be found in the Symmetry Operators chapter. So, our \((-x,-y,z)\) entry will be interpreted as:
Checking the list of the symmetry operators of our space group P4(#75), we see that, this matrix actually corresponds to the 2-fold rotation operator, \(\{2_{001}|0\}\):
In other words, the matrix represantation of a symmetry operator directly tells how it’s going to map a general position of \((x,y,z)\)!
Site Symmetry¶
Consider the two sites \(x_3 = (0.5,0.5,0.4)\) and \(x_4 = (0,0,0.9)\) we had in our case above: Both had orbits containing only themselves (hence a multiplicity of 1) but they were assigned different Wyckoff positions (\(x_3\) was assigned to ‘1b’, whereas \(x_4\) was assigned to ‘1a’). Let’s investigate two other sites: \(x_5 = (0.5,0.5,0.7)\) and \(x_6 = (0,0,0.4)\). Doing similar calculations we had before, we’ll see that \(x_5\) will be assigned to “1b” and \(x_6\) to “1a”. Where does this distinction come from? The answer lies in the site symmetry concept which incorporates the lattice translations.
Applying the 4 symmetry operators to \(x_4\) via the code above we saw that all of them mapped the site to itself:
Calculating the orbit of:
( 0.00, 0.00, 0.90)
Applying o1 yields: ( 0.00, 0.00, 0.90) -> (0.00,0.00,0.90)
Applying o4p yields: ( 0.00, 0.00, 0.90) -> (0.00,0.00,0.90)
Applying o2 yields: ( 0.00, 0.00, 0.90) -> (0.00,0.00,0.90)
Applying o4m yields: ( 0.00, 0.00, 0.90) -> (0.00,0.00,0.90)
However, for \(x_3\), further translations to the reference unit cell were necessary:
Calculating the orbit of:
( 0.50, 0.50, 0.40)
Applying o1 yields: ( 0.50, 0.50, 0.40) -> (0.50,0.50,0.40)
Applying o4p yields: (-0.50, 0.50, 0.40) -> (0.50,0.50,0.40)
Applying o2 yields: (-0.50,-0.50, 0.40) -> (0.50,0.50,0.40)
Applying o4m yields: ( 0.50,-0.50, 0.40) -> (0.50,0.50,0.40)
For the \(4^+\) we actually acquired \((-0.50, 0.50, 0.40)\) and from there applied a translation of \((1,0,0)\) for its equivalent site in the reference unit cell.
For the \(2\) we actually acquired \((-0.50, -0.50, 0.40)\) and from there applied a translation of \((1,1,0)\) for its equivalent site in the reference unit cell.
For the \(4^-\) we actually acquired \((0.50, -0.50, 0.40)\) and from there applied a translation of \((0,1,0)\) for its equivalent site in the reference unit cell.
Site symmetry group contains the operators that map a site onto itself. In our case, the site symmetry group for \(x_4\) consists of: \(\{1,4^+,2,4^-\}\) whereas \(x_3\)’s site symmetry group is made of: \(\{1,\{4^+|100\},\{2|110\},\{4^-|010\}\}\). Notice that, both sets satisfy the group axioms, hence they are groups by themselves. As they have different set of elements, they are assigned to different Wyckoff positions. In other words: positions that are assigned to the same Wyckoff positions belong to the same symmetry group.
Caution
Two sites having the same Wyckoff position doesn’t necessarily mean that they belong to the same orbit (i.e., one can be accessed from the other by means of applying symmetry operations of the group). In our case, even though \(x_3\) and \(x_5\) are both “1b” (likewise, both \(x_4\) and \(x_6\) being “1a”) Wyckoff positions, they belong to different orbits (of size 1).
Fractional Coordinates¶
As you have seen, the positions of the sites throughout this chapter ranged from 0 to 1, with 0 corresponding to the origin and 1 marking the edge of our unit cell. This interpretation is also valid in unit cells that are not necessarily cubic. Consider a triclinic unit cell with the given lattice parameters:
a = 5.44700 Å
b = 7.35076 Å
c = 9.56728 Å
α = 82.1683°
β = 79.4250°
γ = 46.3873°
α denotes the angle between the b and c lattice vectors; β the angle between the a and c; γ between the a and c.
Going between the parametric form and the vectorial form¶
One can go between the parametric form \((a,b,c,\alpha,\beta,\gamma)\) and the vectorial form of:
Note
The conversion between the two forms is basic geometry: the lattice parameters \((a,b,c)\) are the lengths of the vectors \(\left(\vec a,\vec b,\vec c\right)\), and the angles can be derived using the definition of dot product, e.g.,
Going from the parametric form to the vectorial is a little bit complicated: one first defines the \(\vec a\) vector as \((a_1,0,0)\) with \(|\vec a| = a = a_1\), then derives the components of the \(\vec b\) \((b_1,b_2)\) in relation to the \(\vec a\) vector, with the information on the length of \(\vec b\) being b from these two equations:
and for \(\vec c = (c_1,c_2,c_3)\) a similar procedure is followed:
Here are the sample functions that does these conversions:
def par2vec(par):
# Converts from parametric form to vectorial form
# Lattice parameters are input as a 6 dimensional vector
# Returns the lattice vectors as a 3x3 matrix, with each
# vector as a row-vector
a = par[0]
b = par[1]
c = par[2]
alpha = np.deg2rad(par[3])
beta = np.deg2rad(par[4])
gamma = np.deg2rad(par[5])
a1 = a
vec_a = np.array([a1,0,0])
b1 = b*np.cos(gamma)
b2 = np.sqrt(b**2 - b1**2)
vec_b = np.array([b1,b2,0])
c1 = c*np.cos(beta)
c2 = ((b*c*np.cos(alpha))- b1*c1)/b2
c3 = np.sqrt(c**2 - c1**2 - c2**2)
vec_c = np.array([c1,c2,c3])
abc_vec = np.array([vec_a,vec_b,vec_c])
return abc_vec
par1 = np.array([5.44700,7.35076,9.56728,82.1683,79.4250,64.3873])
vec1 = par2vec(par1)
print(vec1)
[[5.447 0. 0. ]
[3.17762796 6.62845028 0. ]
[1.75581063 0.60401363 9.38536857]]
def vec2par(vec):
# Converts from vectorial form to parametric form
# Vectors are input as a 3x3, with each lattice
# as a row-vector
# Returns the lattice parameters as a
# 6 dimensional vector
vec_a = vec[0,:]
vec_b = vec[1,:]
vec_c = vec[2,:]
a = np.linalg.norm(vec_a)
b = np.linalg.norm(vec_b)
c = np.linalg.norm(vec_c)
alpha = np.rad2deg(np.arccos(np.dot(vec_b,vec_c)/(b*c)))
beta = np.rad2deg(np.arccos(np.dot(vec_a,vec_c)/(a*c)))
gamma = np.rad2deg(np.arccos(np.dot(vec_a,vec_b)/(a*b)))
par = np.array([a,b,c,alpha,beta,gamma])
return par
par1 = vec2par(vec1)
print(par1)
[ 5.447 7.35076 9.56728 82.1683 79.425 64.3873 ]
Caution
When deriving the vectorial notation, we first fixed \(\vec{a}\) on a coordinate axis, then defined \(\vec b\) in relation with \(\vec{a}\), and then \(\vec{c}\) wrt \(\vec a\) and \(\vec b\). But in reality, we could choose any starting orientation we wished – so, any 3 vectors obtained by arbitrarily rotating the constructed \(\vec a,\vec b,\vec c\) will also be compatible. This means that there is an infinite number of possible lattice vectors that correspond to the given lattice parameters.
Metric Tensor¶
Fractional coordinates are very useful in simplifying the calculations and visualizing the relations. For example, the site with the coordinates \(\left(\tfrac{1}{2},\tfrac{1}{2},\tfrac{1}{2}\right)\) is always at the center of the unit cell and two sites with coordinates like \(\left(0.3,0.4,0.7\right)\) and \(\left(2.3,1.4,0.7\right)\) are in the unit cells that are 2 unit cells forward along the \(a\) direction and 1 unit cell along the \(b\) direction. If we had used Cartesian coordinates, it wouldn’t be this easy to identify them.
Using the lattice cell parameters defined above (a = 5.44700 Å, b = 7.35076 Å, c = 9.56728 Å, α = 82.1683°, β = 79.4250°, γ = 46.3873°), the center site location in Cartesian coordinates would be: \((2.7235, 3.6754, 4.7836)\), with the units given in Å. And the two sites’ coordinates would be: \((1.6341, 2.9403, 6.6971)\) and \((12.5281, 10.2911, 6.6971)\), respectively (in Å).
As you can see, working with fractional coordinates makes much more sense and they are easier. However, when we want to compare observational data, we need a way to easily convert our results to real life results. We define metric tensor to act as a bridge between the results obtained using fractional coordinates and the real life.
Metric tensor (\(G\)) is a \(3\times3\) symmetrical matrix defined as:
and it is used to convert vectors in the fractional system to real life system, with the dot product between two vectors re-defined as:
A sample function to construct the metric tensor from lattice vectors is given below:
def metric_from_vector(vec):
# Constructs the metric tensor from the input lattice vectors
# Vectors are input as a 3x3, with each lattice vector given
# as a row-vector
# Returns the metric tensor as a 3x3 matrix
G = np.empty((3,3))
for i in range(3):
for j in range(i,3):
G[i,j] = G[j,i] = np.dot(vec[i,:],vec[j,:])
return G
and the accompanying parameters_from_metric() function:
def parameters_from_metric(G):
a = np.sqrt(G[0,0])
b = np.sqrt(G[1,1])
c = np.sqrt(G[2,2])
alpha = np.rad2deg(np.arccos(G[1,2]/(b*c)))
beta = np.rad2deg(np.arccos(G[0,2]/(a*c)))
gamma = np.rad2deg(np.arccos(G[0,1]/(a*b)))
return np.array([a,b,c,alpha,beta,gamma])
Here are two applications of the metric tensor:
Calculation of distance between two points:¶
Suppose we want to calculate the actual distance between the \((0.1,0.2,0.3)\) and \((0.6, 0.1,0.4)\) sites (given obviously in fractional coordinate system) for the triclinic unit cell of parameters \((5.4470, 7.3508, 9.5673, 82.1683, 79.4250, 64.3873)\)
# Step 1: Convert the lattice parameters to vectorial form
# and construct the metric tensor
par1 = np.array([5.4470,7.3508,9.5673,82.1683,79.4250,64.3873])
vec_abc = par2vec(par1)
G = metric_from_vector(vec_abc)
# Step 2: Calculate the displacement vector (in fractional system)
# for the two sites
r1 = np.array([0.1,0.2,0.3])
r2 = np.array([0.6,0.1,0.4])
displacement_vector = r2 - r1
# We are being asked the length of this displacement vector
# in real life. Normally, its length would be the square root
# of the dot product of the vector with itself, but to connect
# it to the real world, we insert the metric tensor in between
# the dot multiplication
distance = np.sqrt(np.linalg.multi_dot(\
(displacement_vector,G,displacement_vector.T)))
print("The distance between the two sites is: {:.5f} Å"\
.format(distance))
The distance between the two sites is: 2.81194 Å
and here is the same data obtained using VESTA:
Calculation of the angle between the two position vectors:¶
This time let’s calculate the actual angle between the two position vectors connecting the two sites in the above example to the origin (using the same lattice parameters):
From the formula, we know that, the angle between two vectors can be found via the definition of dot product:
where, the length of a vector is defined as:
So, inserting the metric tensor into each of the dot products, we have:
(\(\theta'\) indicates the actual, real life angle as opposed to the angle in fractional system)
par1 = np.array([5.4470,7.3508,9.5673,82.1683,79.4250,64.3873])
vec_abc = par2vec(par1)
G = metric_from_vector(vec_abc)
r1 = np.array([0.1,0.2,0.3])
r2 = np.array([0.6,0.1,0.4])
cos_theta_p = np.linalg.multi_dot((r1,G,r2.T))\
/np.sqrt(np.linalg.multi_dot((r1,G,r1.T))\
*np.linalg.multi_dot((r2,G,r2.T)))
theta_p = np.rad2deg(np.arccos(cos_theta_p))
print("The angle between the two vectors is: {:.2f}°"\
.format(theta_p))
The angle between the two vectors is: 22.87°
Same calculation done with VESTA:
You can download and further experiment on this hypothetical triclinic case in VESTA: here (If you don’t readily have VESTA, you can download it from: http://jp-minerals.org/vesta/en/download.html)
Another practical application of the metric tensor is, the square root of its determinant directly being the volume of the cell.