What's the right way to calculate the principal moment of inertia?

Question

I am writing a program that incorporates calculating the principal moment of inertia for a protein residue based on its component atom XYZ coordinates. I am exceedingly confused about which formulas to use in calculating principal moment of inertia for my situation.

Thus far my program does the following:

1. Calculate elements (Ixx, Iyy, Izz, etc.) of symmetric matrix A
2. Find the eigenvalues (and thus the principal moments of inertia)

The formulas I originally used to calculate the elements in step 1 were from pg. 5 of this document. However looking at Wikipedia, the formulas for the symmetric matrix elements are different.

I have also been scouring the internet and found open source code where the center of mass is subtracted from the x, y, and z coordinates before beginning calculation of elements Ixx, Iyy, Izz, etc.

Which formulas/algorithm do I use to calculate principal moment of inertia for my case? I'm not sure where to begin with picking which formulas to use. Is there a source that is accessible to those with a weak physics background that will help me understand which formula/algorithm to use in calculating principal moment of inertia?

Answer 1

The formulae from the NASA document and Wikipedia are simply in different frames: the ones in Wikipedia assume coordinate frame whose axes pass through the center of mass whereas the NASA document assumes arbitrary cartesian coordinate frame.

Let capital X_i, Y_i and Z_i denote coordinates of i-th atom in an arbitrary Cartesian coordinate frame O and x_i, y_i and z_i denote respective coordinates in the frame O_com translated to the center of mass of the molecule.

The coordinates of the center of mass of the molecule in the O frame are:

\begin{equation} X_{com} = \frac{1}{M}\sum_{i} m_i X_i \\ Y_{com} = \frac{1}{M}\sum_{i} m_i Y_i \\ Z_{com} = \frac{1}{M}\sum_{i} m_i Z_i \end{equation}

Transformation from O to O_com is simply translation by [-X_com, -Y_com, -Z_com]:

\begin{equation} x_i = X_i - X_{com} = X_i - \frac{1}{M}\sum_{i} m_i X_i\\ y_i = Y_i - Y_{com} = Y_i - \frac{1}{M}\sum_{i} m_i Y_i\\ z_i = Z_i - Z_{com} = Z_i - \frac{1}{M}\sum_{i} m_i Z_i \end{equation}

Now, take for example the formula for I_xx in the O frame:

\begin{equation} I_{xx} = \sum_{i} m_i (y_i^2 + z_i^2) \end{equation}

(Note that the two-term sum in the parenthesis is simply the Euclidean distance between the i-th atom and the x axis, exactly as one would expect.)

Substituting the formulae for y_i and z_i above we have

\begin{equation} I_{xx} = \sum_{i} m_i [(Y_i - Y_{com})^2 + (Z_i - Z_{com})^2] \\ I_{xx} = \sum_{i} m_i (Y_i^2 + Z_i^2) + \sum_{i} m_i (Y_{com}^2 + Z_{com}^2) - 2 \sum_{i} m_i (Y_i Y_{com} + Z_i Z_{com}) \\ I_{xx} = \sum_{i} m_i (Y_i^2 + Z_i^2) + M (Y_{com}^2 + Z_{com}^2) - 2 M (Y_{com}^2 + Z_{com}^2) \\ I_{xx} = \sum_{i} m_i (Y_i^2 + Z_i^2) - \frac{1}{M}(\sum_{i} m_i Y_i)^2 - \frac{1}{M}(\sum_{i} m_i Z_i)^2 \end{equation}

which is exactly the formula for I_xx in the NASA document. Analogously for the other formulae.

Interestingly, the sum of the last two terms in the final formula is exactly what one would expect from the parallel axis theorem.

The conclusion is that you should use whichever formulae fit your coordinate frame. If the center of your coordinate frame coincides with the center of mass of the molecule, then you can use the simpler formulae from Wikipedia. Otherwise you should use those from the NASA document.