fix: correct atomic thermochemistry and radical symmetry detection - #233
fix: correct atomic thermochemistry and radical symmetry detection#233tdpham2 wants to merge 1 commit into
Conversation
Codex Review SummaryThis comment shows the latest Codex review activity on this pull request.
ℹ️ About Codex in GitHubYour team has set up Codex to review pull requests in this repo. Reviews are triggered when you
Codex reacts with 👀 while any review is running, comments if it has suggestions, and reacts with 👍 once all reviews finish with no findings. |
keceli
left a comment
keceli
left a comment
There was a problem hiding this comment.
The physics here is right and the diagnosis is correct. Monatomic thermo returning entropy = 0.0 with H = G = E_pot at every T and P was straightforwardly wrong, and routing it through IdealGasThermo with geometry="monatomic" gives the proper Sackur-Tetrode translational term. The charge_spin_check=False reasoning is right too — PointGroupAnalyzer reads geometry only, so the check was rejecting OH and NO over a property it doesn't use.
I checked the restructuring specifically, since dedenting thermo out of the vibration block is where this could have broken: all_energies is initialized before either branch, final_structure is in scope from line 706, the relocated block is still inside the outer try ending at 984, and vibrational_frequencies output is byte-identical to main. 15 of the 17 new tests fail on main, so the suite does gate the fix.
Requesting changes on one regression, plus a set of things worth folding in while this is open.
The blocker
temperature is Optional[float] = None with no validator, and the monatomic short-circuit you removed was the only thermo path that never read it. Details inline at ase_core.py:892 — short version, driver="thermo" on a single atom with no temperature returns status: "success" on main and TypeError: must be real number, not NoneType on this branch. That's an agent pulling an atomic reference for an atomization energy, which is the workflow this fix exists to enable, and temperature is the field an LLM is most likely to omit since the schema gives it no default. All 17 new tests pass an explicit temperature.
Already ruled out, so you don't re-spend the time
data-ev-k->cell.dataset.evKis the correct camelCase mapping, and the converter reads from the dataset so toggling is idempotent.symmetrynumber=1for monatomic is required, not merely harmless —get_entropyraises ifsigmaorspinis None.Vibrations.run()andwrite_mode()leaveatoms.positionsuntouched, so the relocated thermo block sees the optimized geometry.get_symmetry_numberreturns the correct sigma for every non-monatomic species tested.run_ase_coreis the only thermochemistry producer; parsl/MACE/academy all delegate to it.
One environment note unrelated to the diff: ruff check . passes on the long new lines only because E501 isn't in the selected rule set (select = ["E4","E7","E9","F"]); ruff format --check would reformat them.
🤖 Findings generated with Claude Code (Opus 5) and posted by @keceli.
| spin=(multiplicity - 1) / 2.0, | ||
| ) | ||
| thermo_data = { | ||
| "enthalpy": float(thermo.get_enthalpy(temperature=temperature)), |
There was a problem hiding this comment.
Regression. temperature is Optional[float] = Field(default=None) (ase_input.py:326) with no validator requiring it for driver="thermo". The monatomic branch you removed was the only thermo path that never read it — on main it returned potential_energy and 0.0 directly (ase_core.py:867-873). Now every monatomic run reaches ASE, which formats '... T = %.2f K:' % temperature and raises.
run_ase_core(ASEInputSchema(input_structure_file='cu.xyz',
calculator={'calculator_type': 'emt'}, driver='thermo'))
# main: {'status': 'success', 'thermochemistry': {'enthalpy': 3.51, ...}}
# here: {'status': 'failure', 'error_type': 'TypeError',
# 'message': 'must be real number, not NoneType'}That shape — an atomic reference with no explicit temperature — is exactly what an agent computing an atomization energy produces, and the schema gives the model no default to anchor on. All 17 new tests pass an explicit temperature, so none cover it.
A 298.15 default on the field, or a model_validator requiring temperature when driver == "thermo", would close it.
Separately, on these same three getters: ASE defaults them to verbose=True, writing ~2 KB / 53 formatted lines to raw stdout per call (get_gibbs_energy re-invokes the other two internally, so enthalpy and entropy each print and compute twice). ThermoChem._vprint is a bare sys.stdout.write, which the stderr logging setup doesn't capture. mcp_tools.py guards its calls with redirect_stdout, but mcp/ase_mcp_hpc.py calls run_ase_core three times with no guard, and server_utils.py:31 defaults to --transport stdio — where stdout is the JSON-RPC channel. This is pre-existing for polyatomic thermo, but the diff newly routes monatomic runs, previously silent, through it. verbose=False on all three plus G = H - T*S fixes both the noise and the double computation.
| # IdealGasThermo expects total spin S; calculators expose 2S+1. | ||
| multiplicity = getattr(calc_model, "get_multiplicity", lambda: None)() or 1 | ||
| thermo = IdealGasThermo( | ||
| vib_energies=all_energies, |
There was a problem hiding this comment.
This hands IdealGasThermo the full 3N complex array rather than the mode_indices subset the same function computed 135 lines above.
_vibrational_mode_indices selects range(6, 3N) in eigenvalue order. ASE 3.25 instead does vib_energies.sort(key=np.abs) then takes [-(3N-6):] — magnitude order — and then _clean_vib_energies(ignore_imag_modes=False) raises on anything imaginary. The two selections disagree whenever a genuine imaginary mode is smaller in magnitude than a trans/rot residual.
Reproduced: square-planar Cu4 + EMT + driver='thermo' returns {'status':'failure','error_type':'ValueError','message':'Imaginary vibrational energies are present.'}, discarding the whole run including the optimized structure and the frequencies that would otherwise have been written. When the imaginary mode is the smaller one it silently succeeds but integrates a different mode set than the one reported in vibrational_frequencies.
Pre-existing verbatim, but this diff rewrote and relocated the call, and the correct subset is already in scope:
vib_energies=[all_energies[i] for i in mode_indices],| symmetrynumber = get_symmetry_number(final_structure) | ||
|
|
||
| # IdealGasThermo expects total spin S; calculators expose 2S+1. | ||
| multiplicity = getattr(calc_model, "get_multiplicity", lambda: None)() or 1 |
There was a problem hiding this comment.
or 1 silently forces a singlet, so the open-shell atoms and radicals this PR targets get zero electronic entropy on most calculators.
Only mace/tblite/orca/nwchem/fairchem define get_multiplicity. EMTCalc and AIMNET2Calc have no multiplicity field at all, and mace_parsl_schema.mace_input_schema has neither multiplicity nor charge (parsl_tools.py:61 hardcodes calculator_type='mace_mp'), so the entire MACE parsl/MCP surface can never set it.
Measured: a Cu atom at 298.15 K gives S = 0.00166373 eV/K = 160.5 J/(mol K). Literature Cu(g) is 166.4; adding the omitted kB*ln 2 gives 166.29. Per-atom error in G runs 0.018 eV (H, Cl, Cu doublets) to 0.036 eV (N quartet) — one-directional, so it accumulates in atomization and BDE workflows with no warning.
Worth either defaulting from the element's ground-state term symbol, or logging when a thermo run falls back to the singlet so the omission is visible.
There's also a consistency hazard one layer down, now that the new docs paragraph formalizes this coupling: nwchem_calc.py:116 reports the multiplicity here even when get_calculator drops it from the input deck (theory='scf', or any theory outside {dft,mp2,ccsd,tce,tddft}), and mace_calc.py:211 returns {} from get_atoms_properties() for every type except mace_polar. Both pair a spin-agnostic enthalpy with a spin-shifted entropy.
| "frequencies": [], | ||
| "frequency_unit": "cm-1", | ||
| } | ||
| if driver in {"vib", "thermo", "ir"} and not ( |
There was a problem hiding this comment.
Two things about this guard.
It special-cases one driver instead of the physical fact. driver="vib" and driver="ir" on a single atom still run the full finite-difference loop. Verified on a Cu atom: vib returns success with energies: [] after 7 calculator evaluations that _vibrational_mode_indices then discards at line 243 (if num_atoms == 1: return []), and ir fails with PropertyNotImplementedError: dipole property not implemented after another 7. With ORCA/NWChem that's 14 wasted SCFs or subprocess launches. docs/calculators.md:42 states single atoms skip finite-difference vibrations, which only holds for thermo. The general condition is len(atoms) > 1, the same shape already used for the optimizer at line 684.
Skipping the block also skips its cleanup. freq_file.unlink() and the rewrite at 764-770, and for f in glob.glob(f'{mol_stem}_vib.*.traj'): os.unlink(f) at 783-787, live inside this block and nowhere else.
Reproduced: run thermo on a Cu2 dimer via mol.xyz in a log dir, overwrite mol.xyz with a single atom, run thermo again — frequencies_mol.csv and mol_vib.5.traj from the dimer are still there, where main unlinked both. main_interface.py:1602 falls back to _latest_artifact_path(log_dir, 'frequencies*.csv') and artifacts.py:126 buckets *_vib.*.traj, so the artifact panel shows the previous molecule's modes beside the atom's result. Hoisting the cleanup above the guard would fix it.
| molecule = aaa.get_molecule(atoms) | ||
| # Rotational symmetry depends on geometry, not the calculator's electronic | ||
| # state. Reconstructed Atoms lack magnetic moments for radical species. | ||
| molecule = aaa.get_molecule(atoms, charge_spin_check=False) |
There was a problem hiding this comment.
The fix is correct — on pymatgen 2025.3.10, OH/NO/CH3 raise ValueError: Charge of 0 and spin multiplicity of 1 is not possible without this flag and return sigma 1/1/6 with it.
Two durability notes on the same function:
charge_spin_check isn't a declared parameter of get_molecule. The signature is get_molecule(atoms, cls=Molecule, **cls_kwargs), so this rides undeclared kwargs through into Molecule.__init__. Nothing at the call site validates it, and pyproject.toml pins a bare "pymatgen" with no floor — so a release that stops forwarding unknown kwargs turns every polyatomic driver='thermo' run and every get_symmetry_number tool call into a TypeError several frames deep in pymatgen. A pymatgen>= floor matching environment.yml, or constructing the Molecule directly, would pin this down.
Monatomic input raises here. get_symmetry_number(AtomsData(numbers=[8], positions=[[0,0,0]])) raises AttributeError: 'PointGroupAnalyzer' object has no attribute 'symmops'. run_ase_core now avoids it via the geometry, symmetrynumber = "monatomic", 1 branch, but this function is also exported as a standalone agent tool (registry/tools.py:133 -> ase_tools.py:118), so an agent reasoning about an atomic reference gets an unhandled AttributeError instead of the obvious 1. The new _MOLECULES table has no single-atom case, and tests/test_tools.py:106 only asserts isinstance(..., int).
Minor, same function: cell and pbc are forwarded into the rebuilt Atoms, but get_molecule discards both — Molecule takes no Lattice. An H2O whose second H sits at x = 10 - 0.76 in a 10 A periodic cell returns sigma = 1 instead of 2, which is kB*ln 2 = 0.0178 eV in T*S at 298 K. Since this is an exported tool, an agent handing it a trajectory frame gets a confidently wrong integer rather than an error.
| label.textContent = unit === 'ev' ? 'eV/K' : | ||
| unit === 'kjmol' ? 'kJ/(mol K)' : 'kcal/(mol K)'; | ||
| }); | ||
| document.querySelectorAll('.entropy-value').forEach(cell => { |
There was a problem hiding this comment.
This branch hardcodes .toFixed(6) while the .energy-value loop it was copied from sets precision = 2 for kJ/mol and kcal/mol (lines 816/819).
So one click on kJ/mol renders the three rows as -1264.63, 0.188869, -1320.93 in the same block — entropy at six decimals beside energies at two.
Entropy genuinely needs more digits than energy here (it's ~1e-3 eV/K), so matching precision exactly would lose information; something like 6 for eV/K and 4 for the molar units would keep both readable.
| <button onclick="toggleEnergyUnit('kcalmol')" data-unit="kcalmol">kcal/mol</button> | ||
| </div> | ||
| <strong>Thermochemistry Values</strong> (<span class="energy-unit">eV</span>):<br> | ||
| <strong>Thermochemistry Values</strong>:<br> |
There was a problem hiding this comment.
The toggle is still labelled "Energy Unit:" with eV / kJ/mol / kcal/mol buttons, but it now switches entropy too — clicking kJ/mol relabels the entropy row to kJ/(mol K).
Button text and produced value disagree, which is the ambiguity this PR set out to remove. Someone clicking kcal/mol to read enthalpy silently rewrites entropy as well; screenshotting the panel while remembering only the button pressed carries off an entropy in the wrong unit.
docs/calculators.md:41 compounds it — "HTML reports also support kJ/(mol K) and kcal/(mol K) for entropy" sends a reader looking for a separate selector that was never added. The doc sentence and this label are the only two places the unit control is described, and they disagree.
Relabelling to just "Units:" would make it accurate.
Also, now that the block-header .energy-unit span is gone: a payload carrying entropy but neither enthalpy nor Gibbs renders the full toggle with no unit label anywhere.
| ): | ||
| vibrations = Mock(side_effect=AssertionError("an atom has no vibrational modes")) | ||
| monkeypatch.setattr(ase.vibrations, "Vibrations", vibrations) | ||
| monkeypatch.setattr(EMTCalc, "get_multiplicity", lambda self: multiplicity, raising=False) |
There was a problem hiding this comment.
raising=False here fabricates a method that doesn't exist — hasattr(EMTCalc, 'get_multiplicity') is False.
So six of the nine parametrized cases exercise a spin EMT can never report. Renaming get_multiplicity across the five calculator schemas that genuinely implement it, or a typo in the getattr key at ase_core.py:882, leaves all nine green while every real EMT run silently falls back to the singlet default.
MaceCalc(multiplicity=...) actually implements it, and the molecular test in this same file already uses that — worth doing the same here so the test exercises the real lookup path.
Two more assertions in this test are weaker than they read:
vibrations.assert_not_called()(line 70) andassert not list(tmp_path.glob("*.traj"))(line 71) can't fail in the regression they guard. Run this file against main and it dies earlier, at line 32 —run_ase_core's blanketexcept Exceptionswallows the Mock'sAssertionErrorand returnsstatus='failure'. Even if reached, main produced zero.trajfor a monatomic thermo run; its only stale artifact wasfrequencies_cu.csv, which a*.csvglob would catch and this one doesn't.result["result"]["thermochemistry"] == output["thermochemistry"](line 34) compares the same dict round-tripped through exact JSON float serialization, so it can only trip on NaN.
The _assert_reference comparison and the analytic potential_energy + 2.5 * kB * T check are both real, so the test does its main job — these are the surrounding assertions.
| tmp_path, monkeypatch, symbols, positions, symmetry, multiplicity, geometry, | ||
| ): | ||
| atoms = Atoms(symbols, positions=positions) | ||
| energies = np.zeros(3 * len(atoms)) |
There was a problem hiding this comment.
This fixture can't catch the mode-selection bug it looks like it covers.
np.zeros(3 * len(atoms)) with only the last mode_count entries nonzero is already in sort(key=np.abs) order and contains no imaginary values. Real Vibrations.get_energies() returns dtype complex128 with the near-zero trans/rot modes stored as pure imaginary (verified on EMT Cu2) — which is precisely what makes ASE's magnitude-based reselection diverge from mode_indices.
Line 128 then builds the reference IdealGasThermo(vib_energies=energies, ...) from that same array, so ASE applies identical slicing to both sides. Changing vib_energies=all_energies to [all_energies[i] for i in mode_indices] — or the reverse — leaves all 17 tests green.
A fixture returning complex energies with imaginary trans/rot entries would make this test discriminating.
Also mode_count = 1 if geometry == "linear" else 3 hardcodes the current atom counts rather than 3N-5 / 3N-6, so adding CO2 to _MOLECULES silently breaks the arithmetic.
| Ideal-gas thermochemistry uses the requested temperature and pressure. Enthalpy | ||
| and Gibbs energy are reported in eV; entropy is reported in eV/K with a separate | ||
| `entropy_unit` field. HTML reports also support kJ/(mol K) and kcal/(mol K) for | ||
| entropy. Single atoms include translational and electronic-spin contributions |
There was a problem hiding this comment.
Two claims in this paragraph don't hold as written.
"Single atoms include translational and electronic-spin contributions without running finite-difference vibrations" is true only for driver="thermo". The guard at ase_core.py:727 is driver == "thermo" and len(atoms) == 1, so vib and ir on a single atom still run the full displacement loop (see that comment).
"Calculator multiplicity determines the electronic-spin contribution" is unreachable on EMT, AIMNet2, and the entire MACE parsl/MCP surface — none of them can report a multiplicity, so those runs always get S = 0 regardless of the species. Worth saying that the contribution is included when the calculator reports a multiplicity, and that it defaults to a singlet otherwise, so readers computing atomization energies with MACE know the term is missing.
2 participants
Summary
Monatomic thermochemistry previously returned zero entropy and set enthalpy/Gibbs energy equal to potential energy at every temperature and pressure. Use ASE IdealGasThermo with monatomic geometry, no vibrations, symmetry 1, and multiplicity-derived total spin. Skip finite-difference vibrations for atomic thermochemistry while retaining empty vibrational output with units.
Geometric symmetry detection now disables pymatgen's charge/spin consistency check, allowing OH and NO radicals without changing the calculator's electronic state. Add
entropy_unit="eV/K"and distinguish entropy from energy in HTML reports, including conversions and legacy result rendering.Related issues
Closes #141. Release prerequisite for #227.
Type of change
How was this tested?
ruff check .andgit diff --checkpassed.Checklist
mainand targetsmainruff check .passespytest tests/ -k "not tblite"passes