The (symmetric) matrix representing the inertia tensor of a collection of masses, $m_i$, with positions $(x_i, y_i, z_i)$ relative to their centre of mass is \begin{align*} \mathbf{I} = \left( \begin{array}{lll} I_{xx} & I_{xy} & I_{xz}\\ I_{xy} & I_{yy} & I_{yz}\\ I_{xz} & I_{yz} & I_{zz} \end{array}\right), \end{align*} where \begin{align*} I_{xx} &= \sum_i m_i(y_i^2 + z_i^2), & \quad I_{yy} &= \sum_i m_i(x_i^2 + z_i^2), & \quad I_{zz} &= \sum_i m_i(x_i^2 + y_i^2),\\ I_{xy} &= -\sum_i m_ix_iy_i, & \quad I_{yz} &= -\sum_i m_iy_iz_i, & \quad I_{xz} &= -\sum_i m_ix_iz_i. \end{align*}
There exists a transformation of the coordinate frame such that this matrix is diagonal: the axes of this transformed frame are called the principal axes and the diagonal inertia matrix elements, $I_a \le I_b \le I_c$ are the principal moments of inertia.
Write a program to calculate the principal moments of inertia of a molecule, given the position and masses of its atoms relative to some arbitrary origin. Your program should first relocate the atom coordinates relative to its centre of mass and then determine the principal moments of inertia as the eigenvalues of the matrix $\mathbf{I}$.
A molecule may be classified as follows according to the relative values of $I_a$, $I_b$ and $I_c$:
$I_a = I_b = I_c$: spherical top;
$I_a = I_b < I_c$: oblate symmetric top;
$I_a < I_b = I_c$: prolate symmetric top;
$I_a < I_b < I_c$: asymmetric top.
Determine the principal moments of inertia of and classify the molecules $\mathrm{NH_3}$, $\mathrm{CH_4}$, $\mathrm{CH_3Cl}$ and $\mathrm{O_3}$ given the data available in the file molecule-data.zip. Also determine the rotational constants, $A$, $B$ and $C$, related to the moments of inertia through $Q = h/(8\pi^2cI_q)$ ($Q=A,B,C; q=a,b,c$) and usually expressed in $\mathrm{cm^{-1}}$.
One solution is presented here.
import sys
import numpy as np
# Atomic mass unit (kg), Planck constant (J s), speed of light (m s-1)
u, h, c = 1.66053886e-27, 6.62606957e-34, 2.99792458e8
def translate_to_cofm(masses, xyz):
# Position of centre of mass in original coordinates
cofm = sum(masses[:,np.newaxis] * xyz) / np.sum(masses)
# Transform to CofM coordinates and return
xyz -= cofm
return xyz
def get_inertia_matrix(masses, xyz):
# Moment of intertia tensor
xyz = translate_to_cofm(masses, xyz)
x, y, z = xyz.T
Ixx = np.sum(masses * (y**2 + z**2))
Iyy = np.sum(masses * (x**2 + z**2))
Izz = np.sum(masses * (x**2 + y**2))
Ixy = -np.sum(masses * x * y)
Iyz = -np.sum(masses * y * z)
Ixz = -np.sum(masses * x * z)
I = np.array([[Ixx, Ixy, Ixz], [Ixy, Iyy, Iyz], [Ixz, Iyz, Izz]])
return I
def get_principal_moi(I):
Ip = np.linalg.eigvals(I)
# Sort and convert principal moments of inertia to SI (kg.m2)
Ip.sort()
return Ip
def classify_molecule(A, B, C):
if np.isclose(A, B):
if np.isclose(B, C):
return 'Spherical top'
return 'Oblate symmetric top'
if np.isclose(B, C):
return 'Prolate symmetric top'
return 'Asymmetric top'
def read_xyz(filename):
try:
data = np.loadtxt(filename, skiprows=2)
except FileNotFoundError:
print('No such file:', filename)
sys.exit(1)
except ValueError as e:
print('Malformed data in {}: {}'.format(filename, e))
sys.exit(1)
return data[:,0], data[:,1:]
try:
masses, xyz = read_xyz(sys.argv[1])
except IndexError:
print('Usage: {} <xyz filename>'.format(sys.argv[0]))
sys.exit(1)
I = get_inertia_matrix(masses, xyz)
Ip = get_principal_moi(I)
Ip *= u / 1.e20
A, B, C = h / 8 / np.pi**2 / c / 100 / Ip
rotor_type = classify_molecule(A, B, C)
print('{}: A={:.6f}, B={:.6f}, C={:.6f} cm-1'.format(rotor_type, A, B, C))
The key function is get_principal_moi which returns the sorted eigenvalues of the inertia matrix. When run from the command line, the output is:
$ python ex6-9-e-molecule-mofi.py CH4.dat
Spherical top: A=5.307906, B=5.307906, C=5.307906 cm-1
$ python ex6-9-e-molecule-mofi.py CH3Cl.dat
Prolate symmetric top: A=5.131101, B=0.439604, C=0.439602 cm-1
$ python ex6-9-e-molecule-mofi.py NH3.dat
Oblate symmetric top: A=9.965800, B=9.965783, C=6.340151 cm-1
$ python ex6-9-e-molecule-mofi.py O3.dat
Asymmetric top: A=3.524891, B=0.444762, C=0.394931 cm-1