Question about the .principal_axes() outputs. Conflicts with document?

[XuyangTsin]

Hello,

I am using The MDAnalysis v2.7, trying to calculate the local membrane norm.
I searched the lipid heads in local area and use the principal_axes to achieve.

def local_memb_norm(universe, component_resid, cutoff):
    local_surface = universe.select_atoms(f"same residue as name PO4 CNC2 COC and around {str(cutoff)} resid {component_resid}")
    local_norm = local_surface.principal_axes()[2]
    return local_norm

According to the theory I understand, the eigenvector which has the smallest eigenvalue should be the one pointing at direction that shows least variations.
Considering the shape of the selected lipid heads is like a plate, the membrane norm should be the one with smallest eigenvalue.
As it is Documented for .principal_axes(), the returned eigenvectors e1,e2,e3 are listed with largest to smallest eigenvalues. That is why I tried to use the e3 as membrane norm.
However, when I check the e3 vector, it is not the one pointing to the positive Z like a membrane norm. In contrast, the e1 is.
Here is a example I took:

U0.select_atoms(f"name PO4 CNC2 COC and around 15 resid 1").principal_axes()

Here is the returned vectors:

array([[-0.11293756, -0.10384797,  0.98816026],
       [-0.11303746, -0.98672347, -0.11661611],
       [ 0.98715127, -0.12486947,  0.09969943]])

You can find the 1st is the one pointing to positive Z, like a membrane norm, instead of the 3rd one.
Here is the coordinates of the lipid heads selection:

U0.select_atoms(f"name PO4 CNC2 COC and around 15 resid 1").positions

array([[ 71.2     ,  37.8     , 149.29999 ],
       [ 49.8     ,  28.5     , 148.      ],
       [ 56.1     ,  21.499998, 144.1     ],
       [ 63.1     ,  18.4     , 146.7     ],
       [ 64.7     ,  23.5     , 147.59999 ],
       [ 82.9     ,  35.      , 148.      ],
       [ 50.2     ,  40.      , 146.79999 ],
       [ 82.      ,  21.      , 147.8     ],
       [ 58.199997,  37.5     , 147.      ],
       [ 72.      ,  20.699999, 147.29999 ],
       [ 69.      ,  24.3     , 150.9     ],
       [ 55.899998,  16.5     , 142.2     ]], dtype=float32)

I did a brief trial to calculate the eigenvalues and eigenvectors using these coordinates:

# Coordinates of the heads
coords = np.array([
    [71.2, 37.8, 149.29999],
    [49.8, 28.5, 148.0],
    [56.1, 21.499998, 144.1],
    [63.1, 18.4, 146.7],
    [64.7, 23.5, 147.59999],
    [82.9, 35.0, 148.0],
    [50.2, 40.0, 146.79999],
    [82.0, 21.0, 147.8],
    [58.199997, 37.5, 147.0],
    [72.0, 20.699999, 147.29999],
    [69.0, 24.3, 150.9],
    [55.899998, 16.5, 142.2]
], dtype=np.float32)

# Compute the covariance matrix
cov_matrix = np.cov(coords, rowvar=False)

# Calculate eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)

# sort from largest to smallest
sorted_indices = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[sorted_indices]
eigenvectors = eigenvectors[:, sorted_indices]

# results
print("Eigenvalues:", eigenvalues)
print("Eigenvectors (principal axes):")
for i, eigenvector in enumerate(eigenvectors.T):
    print(f"Principal Axis a{i+1}: {eigenvector}")

And the output:

Eigenvalues: [126.41967259  70.51007948   3.03457085]
Eigenvectors (principal axes):
Principal Axis a1: [-0.99112536  0.10293388 -0.08411389]
Principal Axis a2: [0.09180748 0.98766212 0.12686575]
Principal Axis a3: [ 0.09613489  0.11801757 -0.98834707]

Regardless of handedness, the vectors a1,a2,a3 are pretty close to ones from .principal_axes(), but in opposite sequence.

'From .principal_axes()'
array([[-0.11293756, -0.10384797,  0.98816026],
       [-0.11303746, -0.98672347, -0.11661611],
       [ 0.98715127, -0.12486947,  0.09969943]])

So what can cause such problem, feels like .principal_axes() did not sort them from largest to smallest eigenvalues?
Or did I misunderstand the math?

Thank you.

[orbeckst]

The code is for principal_axes() is straightforward

    def principal_axes(group, wrap=False):
        atomgroup = group.atoms
        e_val, e_vec = np.linalg.eig(atomgroup.moment_of_inertia(wrap=wrap))

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

        # Make sure the right hand convention is followed
        if np.dot(np.cross(e_vec[0], e_vec[1]), e_vec[2]) < 0:
            e_vec *= -1

        return e_vec

and does exactly what its docs say: calculate eigenvalues and eigenvectors, then sort from highest to lowest eval.

If you look at wikipedia's List_of_3D_inertia_tensors you can see that the one for the solid cuboid of sides $w = d = a$ and height $h$ (i.e., your membrane) is, when diagonalized,

$$ \mathbf{I} = \begin{pmatrix} \frac{1}{12} M (h^2 + a^2) & 0 & 0 \\ 0 & \frac{1}{12} M (h^2 + a^2) & 0 \\ 0 & 0 & \frac{1}{6} M a^2 \end{pmatrix} $$

with the entries representing the eigenvalues for the principal axes along x, y, and z. The last eigenvalue $I_z$ is greater than $I_x = I_y$ when $a > h$, i.e., as long as the x-y extension of your membrane patch is larger than the thickness $h$, the membrane normal should correspond approximately to the principal axes with the largest eigenvalue.

Thus, you made a mistake in your assumption Considering the shape of the selected lipid heads is like a plate, the membrane norm should be the one with smallest eigenvalue.