fix(axes): take the moment of inertia tensor after making the molecule whole

Author: harryswift01

CodeEntropy/levels/axes.py

@@ -819,7 +819,10 @@ def get_vanilla_axes(self, molecule):
         moment of inertia tensor.
 
         The original description is preserved:
-        - The molecule is made whole to ensure correct handling of PBC.
+
+        - The molecule is made whole to ensure correct handling of PBC. The tensor
+          is taken after that, so it does not depend on where the molecule sits
+          relative to the periodic boundary.
         - The moments are obtained by diagonalising the moment of inertia tensor.
         - Eigenvalues are returned sorted from largest to smallest magnitude.
 
@@ -832,8 +835,8 @@ def get_vanilla_axes(self, molecule):
                 - principal_axes: (3, 3) axes.
                 - moment_of_inertia: (3,) moments sorted descending by absolute value.
         """
-        moment_of_inertia_tensor = molecule.moment_of_inertia(unwrap=True)
         make_whole(molecule.atoms)
+        moment_of_inertia_tensor = molecule.moment_of_inertia()
         principal_axes = self.get_principal_axes(molecule)
 
         eigenvalues, _ = np.linalg.eigh(moment_of_inertia_tensor)

tests/unit/CodeEntropy/levels/test_axes.py

@@ -1,5 +1,6 @@
 from unittest.mock import MagicMock
 
+import MDAnalysis as mda
 import numpy as np
 import pytest
 
@@ -392,6 +393,26 @@ def test_get_principal_axes_passes_wrap_to_moment_of_inertia():
     group.atoms.moment_of_inertia.assert_called_once_with(wrap=True)
 
 
+def test_get_vanilla_axes_moments_ignore_position_relative_to_box():
+    """Shifting a molecule by a box vector must not change its moments."""
+    universe = mda.Universe.empty(
+        n_atoms=3, n_residues=1, atom_resindex=[0, 0, 0], trajectory=True
+    )
+    universe.add_TopologyAttr("masses", [16.0, 1.0, 1.0])
+    universe.add_TopologyAttr("bonds", [(0, 1), (0, 2)])
+    universe.dimensions = [50.0, 50.0, 50.0, 90.0, 90.0, 90.0]
+    positions = np.array([[0.0, 0.0, 0.0], [0.96, 0.0, 0.0], [-0.24, 0.93, 0.0]])
+    ax = AxesCalculator()
+
+    universe.atoms.positions = positions + 25.0
+    _, moments = ax.get_vanilla_axes(universe.atoms)
+    universe.atoms.positions = positions + 25.0 + np.array([50.0, 0.0, 0.0])
+    _, shifted_moments = ax.get_vanilla_axes(universe.atoms)
+
+    assert np.allclose(shifted_moments, moments, rtol=1e-3)
+    assert moments.max() < 10.0
+
+
 def test_find_bonded_atoms_selects_heavy_and_hydrogen_groups():
     ax = AxesCalculator()