Classify linearity from COM-relocated inertia tensor - #51
Merged
Conversation
linearity() built its tensor from raw atom positions using the second-moment form (sum of mass-weighted outer products) and tested the smallest eigenvalue with an exact == 0 check. That made the result position-dependent and gave a planar molecule lying in a coordinate plane a spurious zero eigenvalue, so it was misclassified as linear; off-origin linear molecules were misclassified as non-linear. Build the true inertia tensor from center-of-mass-relocated coordinates and use a relative tolerance for the vanishing principal moment. This feeds dof() (3N-5 vs 3N-6), so the fix corrects the degrees of freedom for any geometry that is not pre-centered.
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #51 +/- ##
==========================================
+ Coverage 93.80% 93.83% +0.03%
==========================================
Files 21 21
Lines 1194 1200 +6
==========================================
+ Hits 1120 1126 +6
Misses 74 74
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
galjos
added a commit
that referenced
this pull request
Jun 29, 2026
Found by an audit of the RRHO math against ASE `IdealGasThermo`: the nonlinear path is exact, but **linear and monatomic species were broken**. ## Bug The rotational block in `thermo.py` was hardwired to the nonlinear rigid rotor — it multiplied all three rotational temperatures (`theta_x*theta_y*theta_z`) and used `(3/2)RT` energy and `(3/2)R` heat capacity unconditionally. A linear molecule has one ~zero principal moment of inertia, so `theta -> inf`, `q_rot -> 0`, and: - `S_rot = R(ln q_rot + 3/2) = -inf` → **total entropy = -inf, Gibbs = +inf** for every linear molecule (CO2, N2, O2, CO, HCN, acetylene, …). - `E_rot`/`Cv_rot` were 50% too high for linear (should be `RT`/`R`) and nonzero for monatomic (should be 0). - Monatomic species additionally hit a divide-by-zero in the vibrational formula because `dof()` returned 3 (→ three zero-frequency "modes") instead of 0. The earlier `linearity()` fix (#51) only corrected the geometry *classification*; the Thermo rotational math never consulted it. ## Fix - Branch the rotational contribution on geometry using the existing `linearity()` classifier: `n_rot = 0` (monatomic), `2` (linear), `3` (nonlinear). - General formulas: `S_rot = R(ln q_rot + n_rot/2)`, `E_rot = (n_rot/2)RT`, `Cv_rot = (n_rot/2)R`; linear partition function `q_rot = T/(sigma*theta)`. The nonlinear path is unchanged. - `dof()` returns 0 for a monatomic species (no rotational/vibrational DOF), so the vibrational formulas no longer divide by zero. - Use `np.linalg.eigvalsh` (not `eig`) on the symmetric inertia tensor — real, sorted eigenvalues. ## Validation (vs ASE IdealGasThermo, total entropy cal/(mol*K)) | molecule | geometry | TS | ASE | |---|---|---|---| | H2O | nonlinear | 45.135 | 45.135 | | CO2 | linear | 51.064 | 51.063 | | N2 | linear | 45.768 | 45.768 | | Ar | monatomic | 36.983 | 36.983 | New parametrized tests assert TS == ASE for all four geometry classes (ASE is already a dependency, so this runs in CI with no DFTB+ binary), plus a regression check that linear/monatomic rotational entropy is finite. The nonlinear anthraquinone regression (internal `_rotational_*` attributes) is unchanged. ## Not in this PR (filed separately, both LOW severity) The same audit found two minor issues: imaginary/negative vibrational modes silently produce NaN (`has_imaginary_frequencies` is computed but ignored), and spin/multiplicity is inferred from charge parity only (open-shell neutrals like O2/OH get the wrong electronic entropy).
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #49.
linearity()built its tensor from raw atom positions using the second-moment form (Σ mᵢ rᵢ⊗rᵢ) and tested the smallest eigenvalue with an exact== 0. Two consequences:linearity()feedsdof()(3N−5 for linear vs 3N−6 for non-linear), so a wrong classification corrupts the vibrational degrees of freedom and thus the thermochemistry for any geometry that is not pre-centered.Fix
Build the true inertia tensor (
Σ mᵢ(|rᵢ|²I − rᵢ⊗rᵢ)) from center-of-mass-relocated coordinates — mirroringthermo.py::_compute_inertia_tensor— and decide the vanishing principal moment with a relative tolerance instead of an exact zero.Tests
No existing expectations changed (the old exact-zero check already returned non-linear for real off-axis geometries). 144 passed, 0 skipped (suite run with DFTB+ binaries); pylint 7.98/10.