fix: use scale-relative tolerance for null-mode filtering in vibrational entropy

Author: harryswift01

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

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,26 +27,26 @@
     "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,
+        "united_atom:Rovibrational": 0.001983835543157824,
         "residue:FTmat-Transvibrational": 15.803976930845373,
         "residue:FTmat-Rovibrational": 0.0,
         "united_atom:Conformational": 0.0,
         "residue:Conformational": 0.0,
         "residue:Orientational": 0.0
       },
-      "total": 154.46098837112996
+      "total": 15.80596076638853
     },
     "4": {
       "components": {
@@ -63,7 +63,7 @@
     "5": {
       "components": {
         "united_atom:Transvibrational": 0.0,
-        "united_atom:Rovibrational": 0.010824819910336835,
+        "united_atom:Rovibrational": 0.010824819910336868,
         "residue:FTmat-Transvibrational": 11.05479779559611,
         "residue:FTmat-Rovibrational": 0.0,
         "united_atom:Conformational": 0.0,
@@ -75,7 +75,7 @@
     "6": {
       "components": {
         "united_atom:Transvibrational": 0.0,
-        "united_atom:Rovibrational": 0.00018982532508837594,
+        "united_atom:Rovibrational": 0.00018982532508837537,
         "residue:FTmat-Transvibrational": 19.723925127490933,
         "residue:FTmat-Rovibrational": 0.0,
         "united_atom:Conformational": 0.0,
@@ -99,7 +99,7 @@
     "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,
+        "united_atom:Rovibrational": 0.00023139743934636996,
         "residue:FTmat-Transvibrational": 21.42548529236268,
         "residue:FTmat-Rovibrational": 0.0,
         "united_atom:Conformational": 0.0,
         "residue:Conformational": 0.0,
         "residue:Orientational": 0.0
       },
-      "total": 160.08074429454342
+      "total": 21.425716689802027
     }
   }
 }

tests/unit/CodeEntropy/entropy/test_vibrational_entropy_math.py

@@ -16,17 +16,21 @@ def get_KT2J(self, temp):
     return RM()
 
 
-def test_matrix_eigenvalues_returns_complex_dtype_possible(run_manager):
+def test_matrix_eigenvalues_returns_real_dtype_for_symmetric_matrix(run_manager):
     ve = VibrationalEntropy(run_manager=run_manager)
-    m = np.array([[0.0, -1.0], [1.0, 0.0]])
+    m = np.array([[2.0, -1.0], [-1.0, 2.0]])
     eigs = ve._matrix_eigenvalues(m)
     assert eigs.shape == (2,)
+    assert not np.iscomplexobj(eigs)
 
 
 def test_frequencies_from_lambdas_filters_nonpositive_and_near_zero(run_manager):
+    """Near-zero eigenvalue is scaled to the spectrum's magnitude, not fixed."""
     ve = VibrationalEntropy(run_manager=run_manager)
 
-    lambdas = np.array([-1.0, 0.0, 1e-12, 1.0, 4.0])
+    scale = 4.0
+    noise = np.finfo(float).eps * scale
+    lambdas = np.array([-1.0, 0.0, noise, 1.0, scale])
     freqs = ve._frequencies_from_lambdas(lambdas, temp=298.0)
 
     assert freqs.size == 2