Chemistry Asked by crystallography on October 5, 2021
So let’s say for a given molecule I have the volume of a unit cell, the shape of the unit cell (monoclinic, cubic, etc), the cell angles alpha, beta, and gamma, the cell lengths a,b, and c, and a list of atoms with their corresponding fractional units. How do I convert these fractional coordinates of the atoms into cartesian coordinates for any given cell shape? Orthogonal or not?
Without doubt, the previously provided answers will work well if the crystallographic model is a description in space group triclinic $P1$. However on occasion, the symmetry of the unit cell is higher and the molecule of interest is on a special position, for example a centre of inversion. As a result, the underlying .cif
file need not to contain as many explicit coordinate tuples as atoms in the molecule. On the other hand, on first glance by visual inspection, the model might appear incomplete. This answer aims to complement the previous answers from the perspective not wanting to write new programs.
Entry COD 2102215 of the COD about 2,4,6-trimethoxy-1,3,5-triazine serves as an example. The initial display (of the asymmetric unit) by the COD and the complete molecular structure are provided below:
The .cif
describes only coordinate tuples for seven atom as such, which is exactly a third of the molecule. This is because the threefold symmetry axis of the molecule aligns with a threefold axis of the unit cell. A script like codcif2sdf
, part of the cod-tools collection, however is able to recognize the described symmetries of the unit cell acting on the molecule, and to join the molecular fragments into a .sdf
:
codcif2sdf 2102215.cif > 2102215.sdf
The later format is one many other programs accept as input. Equally, its content may be rewritten in other formats like the most elementary structure format .xyz
, or to yield visualizations e.g., by OpenBabel
obabel -isdf 2102215.sdf -oxyz -O 2102215.xyz
obabel -isdf 2102215.sdf -opng -O 2102215.png
Further details about OpenBabel are available here.
Answered by Buttonwood on October 5, 2021
The answer by porphyrin is correct and arguably the proper way to perform the transformation. But there is another / simpler way to convert between fractional and cartesian coordinates that requires a little less programming.
Provided you have the lattice vectors of the cell $vec{a}_1$, $vec{a}_2$, and $vec{a}_3$ (in cartesian coordinates), you can organize these into a matrix with the lattice vectors stored column-wise:
$$A = begin{bmatrix}vec{a}_1 & vec{a}_2 & vec{a}_3end{bmatrix} = begin{bmatrix} a_{11} & a_{21} & a_{31} \ a_{12} & a_{22} & a_{32} \ a_{13} & a_{23} & a_{33} end{bmatrix}$$
If we denote an atomic position in cartesian space as $vec{x}$ and the corresponding point in fractional (or crystal) coordinates as $vec{y}$, we can convert between the two using the following equation:
$$vec{x} = Avec{y}$$
To convert from cartesian coordinates to fractional coordinates, you can just trivially invert this equation to get:
$$vec{y} = A^{-1}vec{x}$$
So, now if you have a set of $m$ atomic positions in fractional coordinates, you can organize them into an $mtimes 3$ matrix denoted as $Y$. To get the corresponding set of atomic positions in cartesian coordinates denoted as $X$ (also an $mtimes 3$ matrix), you can use the following equation:
$$X = (AY^T)^T$$
and similarly to convert a set of cartesian coordinates to fractional coordinates you can use the inverse relation:
$$Y = (A^{-1}X^T)^T$$
In Python you can do the following:
import numpy as np
# Pick your favorite lattice vectors
a1 = np.array([...])
a2 = np.array([...])
a3 = np.array([...])
# Build the matrix of lattice vectors stored column-wise
# and get its inverse
A = np.vstack([a1, a2, a3]).T
A_inv = np.linalg.inv(A)
# A random set of fractional coordinates stored row-wise
Y = np.array([[0.25, 0.75, 0.1],
[0.75, 0.25, 0.3],
[0.25, 0.75, 0.5]]
)
# Compute the cartesian coordinates
X = np.matmul(A, Y.T).T
# Perform the inverse operation to get fractional coordinates
Y_ = np.matmul(A_inv, X.T).T
You should be able to confirm that the arrays Y
and Y_
are the same.
Answered by Stephen Weitzner on October 5, 2021
You have to use a matrix to convert. This is derived in most textbooks on crystallography, such as McKie & McKie 'Essentials of Crystallography'.
The matrix is
$$ M=begin{bmatrix} a & 0 & \bcos(gamma)& bsin(gamma) & 0\ccos(beta) & cn_2 & csqrt{sin^2(beta)-n_2^2} end{bmatrix}$$
where $$n_2=frac{cos(alpha)-cos(gamma)cos(beta)}{sin(gamma)}$$
and $a,b,c $ are the unit cell dimensions and $alpha,beta,gamma$ the angles in radians. The determinant of $M$ is the cell volume.
If $V_1$ is the vector of fractional cell coordinates $x/a,y/b,z/c$ etc. the matrix multiplication calculation $hat d_{12}=(V_1-V_2)cdot M $ and then $L_{12}=sqrt{ hat d_{12}cdot hat d_{12}}$ gives the bond length atom 1 to 2. The angle $theta$ between two bonds 12 and 23 is given by a dot product $cos(theta)=hat d_{12}cdot hat d_{23}/(L_{12}L_{23})$
An example calculation in python is shown below for a triclinic crystal. The @ is matrix multiplication.
import numpy as np
a = 7.55 # cell parameters
b = 4.99
c = 12.50
alpha = 122.5*np.pi/180
beta = (95+18/60)*np.pi/180
gama = (118+54/60)*np.pi/180
V1 = np.array( [-0.2812 , -0.0628 , 0.1928 ] ) # x/a, y/a, z/a
V2 = np.array( [-0.2308 , -0.0972 , 0.2931 ] )
V3 = np.array( [-0.3639 , -0.1913 , 0.3521 ] )
n2 = (np.cos(alpha)-np.cos(gama)*np.cos(beta))/np.sin(gama)
M = np.array([[a,0,0],[b*np.cos(gama),b*np.sin(gama),0],
[c*np.cos(beta),c*n2,c*np.sqrt(np.sin(beta)**2-n2**2)]])
dcm1 = (V1-V2) @ M # row x matrix
L12 = np.sqrt(dcm1 @ dcm1) # sqrt(dot product)
L12
Answered by porphyrin on October 5, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP