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

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?

-

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 Xi, Yi and Zi denote coordinates of i-th atom in an arbitrary Cartesian coordinate frame O and xi, yi and zi denote respective coordinates in the frame Ocom translated to the center of mass of the molecule.

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

$$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$$

Transformation from O to Ocom is simply translation by [-Xcom, -Ycom, -Zcom]:

$$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$$

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

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

(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 yi and zi above we have

$$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$$

which is exactly the formula for Ixx 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.

-
My program is returning numbers far larger than expected (I'm getting 30000 for the eigenvalues when I am expecting something between 0 and 40). What units should I use for atomic mass? I'm using standard atomic weight given on Wikipedia, so for example nitrogen is just 14.0067. – Gbean Jul 26 '12 at 21:24
As always, you can use whatever units you prefer as long as you remain consistent. It is easier if you annotate/comment potential variable declarations and code and with the units you use at any stage of input, computation and output. In particular, when you compare the eigenvalues you get with what you expect the units must also be the same. – Adam Zalcman Jul 26 '12 at 23:19