Principle Axes documentation is confusing (in tutorial) and missing (in the docs)

[kain88-de]

• question https://stackoverflow.com/q/49239475/2207958
• answer https://stackoverflow.com/a/49268247

Expected behaviour

The principal axes coordinates are not transposed.

Actual behaviour

principal axes coordinates are transposed.

Code to reproduce the behaviour

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())
Lambda =  U.T.dot(I.dot(U))

print("center of mass", CA.center_of_mass())
print("moment of inertia", I)
print("principal axes", U)
print("Lambda = U'IU", Lambda)

Currently version of MDAnalysis:

(run python -c "import MDAnalysis as mda; print(mda.__version__)")
0.17.0

EDIT: fixed matrix multiplication

[orbeckst]

Primarily this concerns the tutorial, see https://github.com/MDAnalysis/MDAnalysisTutorial/issues/18, but we should also properly document the AtomGroup.principal_axes() and moment_of_inertia() methods.

[orbeckst]

import numpy as np

import MDAnalysis as mda
from MDAnalysis.tests.datafiles import PSF, DCD

u = mda.Universe(PSF, DCD)
ca = u.select_atoms("name CA")

U = ca.principal_axes()
I = ca.moment_of_inertia()

# NOTE: U.I.UT instead of UT.I.U because of #33
Lambda = U.dot(I.dot(U.T))
print(Lambda)

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

gives

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

[orbeckst]

Or we can simply say that principal_axes() returns UT instead of U:

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

U = UT.T

Lambda = UT.dot(I.dot(U))

[orbeckst]

Btw, we should also use np.eigh() instead of np.eig() in principal_axes() because I is symmetric by construction.

[orbeckst]

Where is the documentation for principal_axes() and for moment_of_inertia()???? It is not part of AtomGroup – is this because it is transplanted from core.topologyattrs?

@richardjgowers how do we get this back into the docs? (center_of_mass() is also missing – this is not good!)

EDIT: I raised #1845

[orbeckst]

The docs are not really wrong but it is confusing that the principal axes are returned as row-vectors and not as column vectors. See my answer https://stackoverflow.com/a/49268247 for details.

[richardjgowers]

@lilyminium does this ring a bell for you?

[lilyminium]

This should probably go into the list of AtomGroup methods I've been procrastinating on. For now, I think it would be good to explicitly say that it returns the transpose of the eigenvector matrix in the docstring. I guess you could also convert it to a tuple or list or iterator to make it more obvious, but a) that would mess up anyone already using U.T and b) numpy.dot will accept tuples and lists anyway (although not iterators).

FYI Mathematica does the same thing with Eigenvectors[] by returning a list of column vectors that looks like an array of row vectors to matrix operations.

[orbeckst]

Let's just state clearly that U.T is returned.

EDIT: see https://www.mdanalysis.org/MDAnalysisTutorial/atomgroups.html#principalaxes