Test regression suite against MDAnalysis pre-release

CodeEntropy/entropy/vibrational.py

@@ -148,7 +148,7 @@ def _matrix_eigenvalues(matrix: np.ndarray) -> np.ndarray:
             Eigenvalues as a NumPy array.
         """
         matrix = np.asarray(matrix, dtype=float)
-        return la.eigvals(matrix)
+        return la.eigvalsh(matrix)
 
     def _convert_lambda_units(self, lambdas: np.ndarray) -> np.ndarray:
         """Convert eigenvalues into SI units using run_manager.
@@ -185,7 +185,10 @@ def _frequencies_from_lambdas(self, lambdas: np.ndarray, temp: float) -> np.ndar
         """Convert eigenvalues to frequencies with robust filtering.
 
         Filters out eigenvalues that are complex, non-positive, or near-zero to
-        avoid invalid frequencies and unstable entropies.
+        avoid invalid frequencies and unstable entropies. The near-zero
+        tolerance scales with the spectrum's own magnitude (as in
+        ``numpy.linalg.matrix_rank``) rather than using a fixed cutoff, since
+        the eigensolver noise floor around a true null mode is not constant.
 
         Args:
             lambdas: Eigenvalues (post unit conversion).
@@ -197,8 +200,11 @@ def _frequencies_from_lambdas(self, lambdas: np.ndarray, temp: float) -> np.ndar
         lambdas = np.asarray(lambdas)
         lambdas = np.real_if_close(lambdas, tol=1000)
 
+        scale = float(np.max(np.abs(lambdas))) if lambdas.size else 0.0
+        atol = lambdas.size * np.finfo(float).eps * scale
+
         valid_mask = (
-            np.isreal(lambdas) & (lambdas > 0) & (~np.isclose(lambdas, 0, atol=1e-7))
+            np.isreal(lambdas) & (lambdas > 0) & (~np.isclose(lambdas, 0, atol=atol))
         )
 
         removed = int(len(lambdas) - np.count_nonzero(valid_mask))

CodeEntropy/levels/axes.py

@@ -149,7 +149,7 @@ def get_residue_axes(
             rot_center = np.array(residue.center_of_mass())
         else:
             make_whole(data_container.atoms)
-            trans_axes = data_container.atoms.principal_axes()
+            trans_axes = self.get_principal_axes(data_container.atoms)
             if len(edge_atom_set) == 1:
                 edge_atom = edge_atom_set[0]
                 rot_center, rot_axes = self.get_terminal_axes(
@@ -223,7 +223,7 @@ def get_residue_axes_from_topology(
             trans_axes = rot_axes
         else:
             make_whole(mol.atoms)
-            trans_axes = mol.atoms.principal_axes()
+            trans_axes = self.get_principal_axes(mol.atoms)
             rot_axes, moment_of_inertia = self.get_vanilla_axes(residue_atoms)
             center = residue_atoms.center_of_mass(unwrap=True)
 
@@ -281,7 +281,7 @@ def get_UA_axes(self, data_container, index: int, res_position):
                 # only the one residue => use principal axes
                 residue = data_container
                 trans_center = data_container.atoms.center_of_mass(unwrap=True)
-                trans_axes = data_container.atoms.principal_axes()
+                trans_axes = self.get_principal_axes(data_container.atoms)
             else:
                 # residue of interest has at least one neighbour
                 if res_position == -1 or res_position == 1:
@@ -421,7 +421,7 @@ def get_UA_axes_from_topology(
             )
         else:
             make_whole(residue_atoms)
-            trans_axes = residue_atoms.principal_axes()
+            trans_axes = self.get_principal_axes(residue_atoms)
 
         center = heavy_atom.position
         rot_axes, moment_of_inertia = self.get_bonded_axes_from_topology(
@@ -611,7 +611,7 @@ def get_terminal_axes(self, residue, edge, dimensions):
         if len(bonded_atoms) == 0:
             # there is only one heavy atom in the residue
             rot_center = edge.position
-            rot_axes = residue.atoms.principal_axes()
+            rot_axes = self.get_principal_axes(residue.atoms)
         else:
             average_bonded = np.zeros(3)
             for bonded_atom in bonded_atoms:
@@ -834,9 +834,9 @@ def get_vanilla_axes(self, molecule):
         """
         moment_of_inertia_tensor = molecule.moment_of_inertia(unwrap=True)
         make_whole(molecule.atoms)
-        principal_axes = molecule.principal_axes()
+        principal_axes = self.get_principal_axes(molecule)
 
-        eigenvalues, _ = np.linalg.eig(moment_of_inertia_tensor)
+        eigenvalues, _ = np.linalg.eigh(moment_of_inertia_tensor)
         order = np.argsort(np.abs(eigenvalues))[::-1]
         moment_of_inertia = eigenvalues[order]
 
@@ -1099,7 +1099,7 @@ def get_custom_principal_axes(
                 - principal_axes: (3, 3) principal axes (rows).
                 - moment_of_inertia: (3,) principal moments.
         """
-        eigenvalues, eigenvectors = np.linalg.eig(moment_of_inertia_tensor)
+        eigenvalues, eigenvectors = np.linalg.eigh(moment_of_inertia_tensor)
         order = np.abs(eigenvalues).argsort()[::-1]  # descending order
         transposed = np.transpose(eigenvectors)  # columns -> rows
         moment_of_inertia = eigenvalues[order]
@@ -1113,6 +1113,66 @@ def get_custom_principal_axes(
 
         return principal_axes, moment_of_inertia
 
+    @staticmethod
+    def _positive_largest_component(
+        vector: np.ndarray, tol: float = 1e-6
+    ) -> np.ndarray:
+        """Return ``vector`` signed so that its largest component is positive.
+
+        Components within ``tol`` of the largest magnitude count as ties, and the
+        first of them decides the sign.
+        """
+        i = np.flatnonzero(np.abs(vector) >= np.abs(vector).max() - tol)[0]
+        return vector if vector[i] > 0 else -vector
+
+    def get_principal_axes(
+        self, group, wrap: bool = False, rel_tol: float = 1e-8
+    ) -> np.ndarray:
+        """Compute reproducible principal axes of an atom group.
+
+        Like MDAnalysis's ``AtomGroup.principal_axes()``, but using ``np.linalg.eigh``
+        and a canonical frame, so the result does not depend on which eigenvectors
+        the LAPACK build happens to return.
+
+        - Axes are sorted by descending eigenvalue and returned as rows.
+        - Non-degenerate axes take the sign that makes their largest component
+          positive, and the last axis completes a right-handed frame.
+        - If two moments agree to within ``rel_tol`` (e.g. a linear bead), the
+          in-plane axes are built from the lab axis most perpendicular to the
+          unique axis. If all three agree, the lab axes are returned.
+
+        Args:
+            group: MDAnalysis atom group (or anything with ``.atoms``).
+            wrap: Whether to wrap atoms into the primary unit cell first.
+            rel_tol: Relative tolerance, scaled by the largest moment, below which
+                two moments are treated as degenerate.
+
+        Returns:
+            np.ndarray: (3, 3) principal axes (rows), right-handed.
+        """
+        e_val, e_vec = np.linalg.eigh(group.atoms.moment_of_inertia(wrap=wrap))
+        order = np.argsort(e_val)[::-1]
+        vals, axes = e_val[order], e_vec[:, order].T
+
+        gap = rel_tol * max(abs(vals[0]), np.finfo(float).tiny)
+        d01 = abs(vals[0] - vals[1]) <= gap
+        d12 = abs(vals[1] - vals[2]) <= gap
+
+        if d01 and d12:
+            return np.eye(3)
+
+        if d01 or d12:
+            n = self._positive_largest_component(axes[2] if d01 else axes[0])
+            i = np.flatnonzero(np.abs(n) <= np.abs(n).min() + 1e-6)[0]
+            u = np.eye(3)[i] - n[i] * n
+            u /= np.linalg.norm(u)
+            v = np.cross(n, u)
+            return np.array([u, v, n] if d01 else [n, u, v])
+
+        a0 = self._positive_largest_component(axes[0])
+        a1 = self._positive_largest_component(axes[1])
+        return np.array([a0, a1, np.cross(a0, a1)])
+
     def get_UA_masses(self, molecule) -> list[float]:
         """Return united-atom (UA) masses for a molecule.
 

CodeEntropy/levels/nodes/covariance.py

@@ -506,7 +506,7 @@ def _build_ua_vectors(
                 # principal axes
                 make_whole(residue.atoms)
                 make_whole(bead)
-                trans_axes = residue.atoms.principal_axes()
+                trans_axes = axes_manager.get_principal_axes(residue.atoms)
                 rot_axes, moi = axes_manager.get_vanilla_axes(bead)
                 center = bead.center_of_mass(unwrap=True)
 
@@ -658,7 +658,7 @@ def _get_residue_axes(
         make_whole(mol.atoms)
         make_whole(bead)
 
-        trans_axes = mol.atoms.principal_axes()
+        trans_axes = axes_manager.get_principal_axes(mol.atoms)
         rot_axes, moi = axes_manager.get_vanilla_axes(bead)
         center = bead.center_of_mass(unwrap=True)
         return (
@@ -688,7 +688,7 @@ def _get_polymer_axes(
         make_whole(mol.atoms)
         make_whole(bead)
 
-        trans_axes = mol.atoms.principal_axes()
+        trans_axes = axes_manager.get_principal_axes(mol.atoms)
         rot_axes, moi = axes_manager.get_vanilla_axes(bead)
         center = bead.center_of_mass(unwrap=True)
 

tests/regression/baselines/benzaldehyde/axes_off.json

@@ -2,15 +2,15 @@
   "groups": {
     "0": {
       "components": {
-        "united_atom:Transvibrational": 0.08982962903796131,
-        "united_atom:Rovibrational": 32.16018134884085,
-        "residue:FTmat-Transvibrational": 88.7671666695003,
-        "residue:FTmat-Rovibrational": 61.61036267672132,
+        "united_atom:Transvibrational": 0.08580917390710577,
+        "united_atom:Rovibrational": 32.7616905666739,
+        "residue:FTmat-Transvibrational": 87.38657089591656,
+        "residue:FTmat-Rovibrational": 62.30339378754512,
         "united_atom:Conformational": 0.0,
         "residue:Conformational": 0.0,
         "residue:Orientational": 20.481571492615355
       },
-      "total": 203.1091118167158
+      "total": 203.01903591665805
     }
   }
 }

tests/regression/baselines/benzaldehyde/combined_forcetorque_false.json

@@ -2,15 +2,15 @@
   "groups": {
     "0": {
       "components": {
-        "united_atom:Transvibrational": 0.07119323721997475,
-        "united_atom:Rovibrational": 49.68669738152346,
-        "residue:Transvibrational": 69.48692941204929,
-        "residue:Rovibrational": 68.46147102540942,
+        "united_atom:Transvibrational": 0.05129761410051273,
+        "united_atom:Rovibrational": 49.686697381523466,
+        "residue:Transvibrational": 69.21750344995603,
+        "residue:Rovibrational": 67.99888878209293,
         "united_atom:Conformational": 0.0,
         "residue:Conformational": 0.0,
         "residue:Orientational": 20.481571492615355
       },
-      "total": 208.1878625488175
+      "total": 207.4359587202883
     }
   }
 }

tests/regression/baselines/benzaldehyde/frame_window.json

@@ -2,15 +2,15 @@
   "groups": {
     "0": {
       "components": {
-        "united_atom:Transvibrational": 40.30267601961045,
-        "united_atom:Rovibrational": 38.21906858443615,
-        "residue:FTmat-Transvibrational": 73.41098578352612,
-        "residue:FTmat-Rovibrational": 57.881504393660364,
+        "united_atom:Transvibrational": 42.23402688550062,
+        "united_atom:Rovibrational": 38.219068584436165,
+        "residue:FTmat-Transvibrational": 73.4310301096772,
+        "residue:FTmat-Rovibrational": 57.470843728209715,
         "united_atom:Conformational": 0.0,
         "residue:Conformational": 0.0,
         "residue:Orientational": 20.481571492615355
       },
-      "total": 230.29580627384846
+      "total": 231.83654080043905
     }
   }
 }

tests/regression/baselines/benzaldehyde/grouping_each.json

@@ -3,7 +3,7 @@
     "0": {
       "components": {
         "united_atom:Transvibrational": 0.0,
-        "united_atom:Rovibrational": 0.062499436482478445,
+        "united_atom:Rovibrational": 0.062499436482478535,
         "residue:FTmat-Transvibrational": 16.83857036769833,
         "residue:FTmat-Rovibrational": 0.0,
         "united_atom:Conformational": 0.0,
@@ -27,79 +27,79 @@
     "2": {
       "components": {
         "united_atom:Transvibrational": 0.0,
-        "united_atom:Rovibrational": 0.0022147884278620462,
-        "residue:FTmat-Transvibrational": 17.449372254125862,
+        "united_atom:Rovibrational": 0.0022147884278620493,
+        "residue:FTmat-Transvibrational": 17.449372254125855,
         "residue:FTmat-Rovibrational": 0.0,
         "united_atom:Conformational": 0.0,
         "residue:Conformational": 0.0,
         "residue:Orientational": 31.909265242174982
       },
-      "total": 49.36085228472871
+      "total": 49.3608522847287
     },
     "3": {
       "components": {
         "united_atom:Transvibrational": 0.0,
-        "united_atom:Rovibrational": 138.65701144028458,
-        "residue:FTmat-Transvibrational": 15.803976930845373,
+        "united_atom:Rovibrational": 0.001983835543157824,
+        "residue:FTmat-Transvibrational": 15.80397693084538,
         "residue:FTmat-Rovibrational": 0.0,
         "united_atom:Conformational": 0.0,
         "residue:Conformational": 0.0,
         "residue:Orientational": 0.0
       },
-      "total": 154.46098837112996
+      "total": 15.805960766388537
     },
     "4": {
       "components": {
         "united_atom:Transvibrational": 0.0,
         "united_atom:Rovibrational": 0.003874883925216773,
-        "residue:FTmat-Transvibrational": 23.899916844406114,
+        "residue:FTmat-Transvibrational": 23.899916844406135,
         "residue:FTmat-Rovibrational": 0.0,
         "united_atom:Conformational": 0.0,
         "residue:Conformational": 0.0,
         "residue:Orientational": 0.0
       },
-      "total": 23.90379172833133
+      "total": 23.903791728331353
     },
     "5": {
       "components": {
         "united_atom:Transvibrational": 0.0,
-        "united_atom:Rovibrational": 0.010824819910336835,
-        "residue:FTmat-Transvibrational": 11.05479779559611,
+        "united_atom:Rovibrational": 0.010824819910336868,
+        "residue:FTmat-Transvibrational": 11.054797795596109,
         "residue:FTmat-Rovibrational": 0.0,
         "united_atom:Conformational": 0.0,
         "residue:Conformational": 0.0,
         "residue:Orientational": 31.909265242174982
       },
-      "total": 42.97488785768143
+      "total": 42.97488785768142
     },
     "6": {
       "components": {
         "united_atom:Transvibrational": 0.0,
-        "united_atom:Rovibrational": 0.00018982532508837594,
-        "residue:FTmat-Transvibrational": 19.723925127490933,
+        "united_atom:Rovibrational": 0.00018982532508837537,
+        "residue:FTmat-Transvibrational": 19.72392512749092,
         "residue:FTmat-Rovibrational": 0.0,
         "united_atom:Conformational": 0.0,
         "residue:Conformational": 0.0,
         "residue:Orientational": 0.0
       },
-      "total": 19.72411495281602
+      "total": 19.724114952816006
     },
     "7": {
       "components": {
         "united_atom:Transvibrational": 0.0,
         "united_atom:Rovibrational": 0.004033941970593011,
-        "residue:FTmat-Transvibrational": 18.90676882142135,
+        "residue:FTmat-Transvibrational": 18.906768821421363,
         "residue:FTmat-Rovibrational": 0.0,
         "united_atom:Conformational": 0.0,
         "residue:Conformational": 0.0,
         "residue:Orientational": 0.0
       },
-      "total": 18.910802763391942
+      "total": 18.910802763391956
     },
     "8": {
       "components": {
         "united_atom:Transvibrational": 0.0,
-        "united_atom:Rovibrational": 0.020362757918507005,
+        "united_atom:Rovibrational": 0.020362757918507036,
         "residue:FTmat-Transvibrational": 16.13204209116967,
         "residue:FTmat-Rovibrational": 0.0,
         "united_atom:Conformational": 0.0,
@@ -111,14 +111,14 @@
     "9": {
       "components": {
         "united_atom:Transvibrational": 0.0,
-        "united_atom:Rovibrational": 138.65525900218074,
-        "residue:FTmat-Transvibrational": 21.42548529236268,
+        "united_atom:Rovibrational": 0.00023139743934636996,
+        "residue:FTmat-Transvibrational": 21.42548529236267,
         "residue:FTmat-Rovibrational": 0.0,
         "united_atom:Conformational": 0.0,
         "residue:Conformational": 0.0,
         "residue:Orientational": 0.0
       },
-      "total": 160.08074429454342
+      "total": 21.425716689802016
     }
   }
 }