From 182cda5b4233e7ddce628d861a63e393cdd5c162 Mon Sep 17 00:00:00 2001 From: harryswift01 Date: Mon, 5 Oct 2026 11:38:33 +0100 Subject: [PATCH] fix: use scale-relative tolerance for null-mode filtering in vibrational entropy --- CodeEntropy/entropy/vibrational.py | 12 +++++++--- .../baselines/benzaldehyde/grouping_each.json | 22 +++++++++---------- .../entropy/test_vibrational_entropy_math.py | 10 ++++++--- 3 files changed, 27 insertions(+), 17 deletions(-) diff --git a/CodeEntropy/entropy/vibrational.py b/CodeEntropy/entropy/vibrational.py index 3e02a01..53de6d3 100644 --- a/CodeEntropy/entropy/vibrational.py +++ b/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)) diff --git a/tests/regression/baselines/benzaldehyde/grouping_each.json b/tests/regression/baselines/benzaldehyde/grouping_each.json index 0a21847..deacf42 100644 --- a/tests/regression/baselines/benzaldehyde/grouping_each.json +++ b/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 } } } diff --git a/tests/unit/CodeEntropy/entropy/test_vibrational_entropy_math.py b/tests/unit/CodeEntropy/entropy/test_vibrational_entropy_math.py index 2554571..07d14f7 100644 --- a/tests/unit/CodeEntropy/entropy/test_vibrational_entropy_math.py +++ b/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