How to use MDAnalysis to principal_axes and moment_of_inertia with a group of atoms?

Question

I am trying to use MDAnalysis (MDAnalysis.__version__ == 0.17.0) API functions principal_axes() and moment_of_inertia() to calculate these matrices for a group of selected atoms as described in the doc:

import MDAnalysis
from MDAnalysis.tests.datafiles import PSF, DCD
import numpy as np

u = MDAnalysis.Universe(PSF, DCD)

CA = u.select_atoms("protein and name CA")

I = np.matrix(CA.moment_of_inertia())
U = np.matrix(CA.principal_axes())
print("center of mass", CA.center_of_mass())
print("moment of inertia", I)
print("principal axes", U)
print("Lambda = U'IU", np.transpose(U)*I*U)

The output:

center of mass [ 0.06873595 -0.04605918 -0.24643682]
moment of inertia [[ 393842.2070687     -963.01376005   -6050.68541811]
 [   -963.01376005  474434.9289629    -3902.61617054]
 [  -6050.68541811   -3902.61617054  520207.91703069]]
principal axes [[-0.04680878 -0.08278738  0.99546732]
 [ 0.01813292 -0.9964659  -0.08201778]
 [-0.99873927 -0.01421157 -0.04814453]]
Lambda = U'IU [[ 519493.24344558   -4093.3268841    11620.96444297]
 [  -4093.3268841   473608.1536763     7491.56715845]
 [  11620.96444297    7491.56715845  395383.6559404 ]]

This looks wrong, one of the reason is that U'IU isn't diagonal as mentioned in the doc:

Maybe I need to align the protein to the center of mass to calculate the moment of inertia with respect to that.

The last line to calculate lamdba fails because `U*I` in numpy is not the matrix product but a element wise product. To calculate the matrix product you need to use `np.dot(U, I)`. — Max Linke, Mar 13 '18 at 11:20
@MaxLinke true, but also when using `np.matmul` that what happens. I didn't know I can use `np.do`t for matrix multiplication. — 0x90, Mar 13 '18 at 12:50

0x90 · Answer 1 · 2019-08-12T04:38:20.467

1

The thing is that in the documentation they say U'IU, but U is the transpose of the return value from CA.principal_axes() (see the source code):

    # Sort
    indices = np.argsort(e_val)[::-1]
    # Return transposed in more logical form. See Issue 33.
    return e_vec[:, indices].T

Matlab to confirm:

>> I=[ 393842.2070687     -963.01376005   -6050.68541811 ;  -963.01376005  474434.9289629    -3902.61617054;  -6050.68541811   -3902.61617054  520207.91703069];
>> U=[-0.04680878 -0.08278738  0.99546732; 0.01813292 -0.9964659  -0.08201778;-0.99873927 -0.01421157 -0.04814453];
>> U*I*U'

ans =

   1.0e+05 *

    5.2082    0.0000   -0.0000
    0.0000    4.7413   -0.0000
   -0.0000   -0.0000    3.9354

edited Aug 12 '19 at 04:38

answered Mar 13 '18 at 05:51

0x90

39,472
36
165
245

This is indeed the issue. I opened an [issue](https://github.com/MDAnalysis/mdanalysis/issues/1828) to fix this. Thanks for diagnosing this yourself. – Max Linke Mar 13 '18 at 11:21
@MaxLinke not sure it's a real issue, because in the source they say it's done intentionally due to issue #33. Can you fix the code to reproduce to use `np.matmul` please? – 0x90 Mar 13 '18 at 12:53
Sorry, the [tutorial docs on AtomGroup](https://www.mdanalysis.org/MDAnalysisTutorial/atomgroups.html#important-methods-and-attributes-of-atomgroup) must not have been updated since before issue [#33](https://github.com/MDAnalysis/mdanalysis/issues/33). We'll update the docs. I also raise [MDAnalysisTutorial #18](https://github.com/MDAnalysis/MDAnalysisTutorial/issues/18). Thanks for the excellent feedback. – orbeckst Mar 13 '18 at 23:55
@orbeckst thank you. I am now a bit confused. Can you please add an answer how to find the moment of inertia and principal axes and eigenvalues? Please note in the issue that we need to use `np.matmul` or `np.dot` and not the `*` operator in issue [#18](https://github.com/MDAnalysis/MDAnalysisTutorial/issues/18). – 0x90 Mar 14 '18 at 00:08
1

@0x90, thanks for noting the `*`, I replaced it with `.dot()` in the two issues. – orbeckst Mar 14 '18 at 00:34

score 1 · Accepted Answer · answered Mar 14 '18 at 00:59

The docs in the tutorial on AtomGroup.principal_axes() are in principle correct but it is confusing that the return value of AtomGroup.principal_axes() is not the matrix U but its transpose, U.T:

The AtomGroup.principal_axes() method returns an array [p1, p2, p3] where the principal axes p1, p2, p3 are arrays of length 3; this layout as row vectors was chosen for convenience (so that one can extract the vectors with p1, p2, p3 = ag.principal_axes()). To form a matrix U where the principal axes are the column vectors as in the usual treatment of the principal axes one has to transpose. For example:

import MDAnalysis
from MDAnalysis.tests.datafiles import PSF, DCD
import numpy as np

u = MDAnalysis.Universe(PSF, DCD)

CA = u.select_atoms("protein and name CA")

I = CA.moment_of_inertia()
UT = CA.principal_axes()

# transpose the row-vector layout UT = [p1, p2, p3]
U = UT.T

# test that U diagonalizes I
Lambda = U.T.dot(I.dot(U))

print(Lambda)

# check that it is diagonal (to machine precision)
print(np.allclose(Lambda - np.diag(np.diagonal(Lambda)), 0))

The matrix Lambda should be diagonal (the last print should show True):

[[ 5.20816990e+05 -6.56706349e-10 -2.83491351e-12]
[-6.62283524e-10  4.74131234e+05 -2.06979926e-11]
[-6.56687024e-12 -2.07159142e-11  3.93536829e+05]]
True

Finally, if you want to calculate "by hand":

values, evecs = np.linalg.eigh(I)
indices = np.argsort(values)
U = evecs[:, indices]

This gives U with the principal axes as column vectors.

How to use MDAnalysis to principal_axes and moment_of_inertia with a group of atoms?

2 Answers2