Fix rotational thermochemistry for linear and monatomic species - #59
Merged
Conversation
The rotational block in Thermo was hardwired to the nonlinear rigid rotor: it multiplied all three rotational temperatures and used (3/2)RT energy and (3/2)R heat capacity unconditionally. For a linear molecule one principal moment of inertia is ~zero, so a rotational temperature went to infinity and the partition function to zero, giving rotational entropy = -inf (and thus -inf total entropy and +inf Gibbs) for every linear molecule (CO2, N2, O2, CO, acetylene, ...). Monatomic species were likewise broken. Branch the rotational contribution on geometry via the existing linearity() classifier: n_rot = 0 (monatomic), 2 (linear), 3 (nonlinear), with S_rot=R(ln q+n_rot/2), E_rot=(n_rot/2)RT, Cv_rot=(n_rot/2)R, and the correct linear partition function q_rot=T/(sigma*theta). Also use eigvalsh (not eig) on the symmetric inertia tensor, and make dof() return 0 for a monatomic species (no vibrational modes) so the vibrational formulas no longer divide by zero. Validated against ASE IdealGasThermo: H2O (nonlinear) 45.135, CO2 (linear) 51.064 vs 51.063, N2 (linear) 45.768, Ar (monatomic) 36.983 cal/(mol*K) — all match; the nonlinear anthraquinone regression is unchanged.
This was referenced Jun 29, 2026
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #59 +/- ##
==========================================
+ Coverage 95.93% 96.15% +0.21%
==========================================
Files 22 22
Lines 1280 1301 +21
==========================================
+ Hits 1228 1251 +23
+ Misses 52 50 -2
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:
|
This was referenced Jun 29, 2026
galjos
added a commit
that referenced
this pull request
Jun 30, 2026
Closes #60. An imaginary (non-positive) vibrational mode in the kept `dof` set made the harmonic-oscillator formulas return `NaN` entropy/Gibbs **silently** (`np.log(1 - exp(-theta/T))` of a negative argument). For the screening framework that means a `NaN` row instead of a clean `error`. ## Fix `_vibrational_contribution` raises `TSValueError` when `real_vibrational_frequencies` contains a non-positive value — a non-minimum geometry is now reported as an error (and isolated by the screening loop) instead of producing `NaN`. ## Why the kept set, not `has_imaginary_frequencies` Real DFTB+ output carries small **negative** translation/rotation frequencies in the *input* (e.g. `frequency.txt` starts `-35.7, -14.34, -3.23`) that `frequency_dof` correctly **drops**. Guarding on `System.has_imaginary_frequencies` (the full input) would wrongly reject every real molecule; the kept set of a true minimum is strictly positive (verified: the anthraquinone fixture keeps 66 modes, min 40.6 cm⁻¹, and still computes S = 103.743). ## Tests A 3-atom system with an imaginary mode in the kept set now raises `TSValueError`; the existing nonlinear/anthraquinone regressions (negative input, positive kept) are unchanged. 163 passed; pylint 8.90. Found by the thermochemistry audit alongside #59; #61 (spin/multiplicity) remains.
galjos
added a commit
that referenced
this pull request
Jul 2, 2026
Closes #61. Electronic entropy (`q_elec = 2S+1`, `S_elec = R ln q`) took its spin from `spin(charge)` — **charge parity** — which is wrong for open-shell species. ## Change - Replace it with `default_spin(atoms, charge)`: the **minimum-spin ground state** from the electron-count parity — `even → singlet` (S=0), `odd → doublet` (S=0.5). Even electron count is a singlet ~99% of the time, and an odd count is essentially always a doublet, so this auto-handles closed-shell molecules **and** simple radicals (e.g. every 1-electron-reduced redox species) correctly. - Add an explicit, validated (non-negative) `spin` argument to `System`, threaded through `run_thermo`, `dftbplus_thermo`, and `screen` (optional `spin` manifest column) for the cases the guess cannot know: **even-electron high-spin ground states** (triplet O2 → `spin=1`, → `S_elec = R ln 3`). ```python System(atoms, charge=0, spin=1.0, ...) # triplet O2 (even electrons, user decides) # radicals (odd electrons) are auto-guessed as doublets ``` ## Tests - `default_spin` electron-count parametrized (water/O2/hydroxide → 0; methyl radical → 0.5); explicit-spin override; negative-spin rejection. - Triplet O2 (`spin=1`) → `S_elec = R ln 3`. - Manifest `spin` column → `ScreeningJob.spin`; `screen` passes it to `dftbplus_thermo`. 172 passed; pylint 8.86. Follow-up (separate PR, once this lands): opt-in **spin-polarised DFTB+** so a declared/guessed open-shell spin also drives the calculation (verified 3ob spin constants ready). Closes out the thermochemistry audit (#59 rotational, #62 imaginary modes, #61 spin).
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.
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.pywas hardwired to the nonlinear rigid rotor — it multiplied all three rotational temperatures (theta_x*theta_y*theta_z) and used(3/2)RTenergy and(3/2)Rheat capacity unconditionally. A linear molecule has one ~zero principal moment of inertia, sotheta -> 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_rotwere 50% too high for linear (should beRT/R) and nonzero for monatomic (should be 0).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
linearity()classifier:n_rot = 0(monatomic),2(linear),3(nonlinear).S_rot = R(ln q_rot + n_rot/2),E_rot = (n_rot/2)RT,Cv_rot = (n_rot/2)R; linear partition functionq_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.np.linalg.eigvalsh(noteig) on the symmetric inertia tensor — real, sorted eigenvalues.Validation (vs ASE IdealGasThermo, total entropy cal/(mol*K))
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_frequenciesis computed but ignored), and spin/multiplicity is inferred from charge parity only (open-shell neutrals like O2/OH get the wrong electronic entropy).