From 28875fb5c1f07be21c6412d28e5073eefb7fdf08 Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Fri, 24 Jul 2026 10:38:00 +0200 Subject: [PATCH] Fix extreme near-linear classification --- ThermoScreening/thermo/system.py | 2 +- docs/benchmarks/qm9_thermochemistry.rst | 36 +++++++------ scripts/validate_qm9.py | 71 ++++++++++++++++++++++--- tests/thermo/test_system.py | 49 ++++++++++++++--- 4 files changed, 125 insertions(+), 33 deletions(-) diff --git a/ThermoScreening/thermo/system.py b/ThermoScreening/thermo/system.py index 29b1d3e..5f4f74a 100644 --- a/ThermoScreening/thermo/system.py +++ b/ThermoScreening/thermo/system.py @@ -125,7 +125,7 @@ def linearity(atoms: List[Atom]) -> bool: eigenvalues = np.linalg.eigvalsh(inertia_tensor) # a linear molecule has exactly one vanishing principal moment of inertia # (the molecular axis); use a relative tolerance rather than an exact zero. - return bool(eigenvalues[0] <= 1e-8 * eigenvalues[-1]) + return bool(eigenvalues[0] <= 1e-9 * eigenvalues[-1]) def dimensionality(atoms: List[Atom]) -> int: diff --git a/docs/benchmarks/qm9_thermochemistry.rst b/docs/benchmarks/qm9_thermochemistry.rst index 3d07954..5193b0d 100644 --- a/docs/benchmarks/qm9_thermochemistry.rst +++ b/docs/benchmarks/qm9_thermochemistry.rst @@ -15,9 +15,10 @@ Run the benchmark from a source checkout: python scripts/validate_qm9.py The command downloads the original 86 MB archive, verifies its MD5 checksum, -and streams a pinned 1,000-molecule sample without extracting the full dataset. -The sample contains the first 100 molecules, 450 evenly spaced records, and 450 -seeded random records. +and streams the dataset without extracting it. It checks the linear/nonlinear +classification of all 133,885 geometries and recomputes thermochemistry for a +pinned 1,000-molecule sample. The sample contains the first 100 molecules, +record 14,564, 450 evenly spaced records, and a seeded random remainder. Method ------ @@ -44,20 +45,20 @@ printed to six decimal places and ``Cv`` to three decimal places. - Mean absolute error - Maximum absolute error * - ``U0`` - - 0.246147 microhartree - - 0.497836 microhartree + - 0.243309 microhartree + - 0.499333 microhartree * - ``U`` - - 0.401046 microhartree - - 1.336198 microhartree + - 0.399690 microhartree + - 1.389329 microhartree * - ``H`` - - 0.404191 microhartree + - 0.394658 microhartree - 1.382253 microhartree * - ``G`` - - 1.622078 microhartree - - 3.211252 microhartree + - 1.612468 microhartree + - 3.207543 microhartree * - ``Cv`` - 0.000259 cal/(mol K) - - 0.000585 cal/(mol K) + - 0.000575 cal/(mol K) The maximum Gibbs-energy difference is approximately 0.0020 kcal/mol. The benchmark therefore reproduces the published quantities to the precision @@ -66,12 +67,13 @@ available from QM9's rounded coordinates, frequencies, and energies. Near-linear geometry -------------------- -QM9 record 25 is a slightly bent cyanogen geometry with six vibrational modes. -Its smallest-to-largest principal-moment ratio is about ``9.7e-8``. A looser -linearity tolerance treated it as exactly linear and expected seven modes. -The external benchmark identified this edge case; the regression test now -requires the geometry to be handled as nonlinear, matching the QM9 -thermochemistry. +QM9 records 25 and 14,564 are slightly bent geometries with six and eighteen +vibrational modes, respectively. Their smallest-to-largest principal-moment +ratios are about ``9.7e-8`` and ``2.9e-9``. Looser linearity tolerances treated +them as exactly linear and expected an extra vibrational mode. The regression +tests require both geometries to be handled as nonlinear. The benchmark also +checks the mode count for every QM9 geometry, including 433 records whose +frequency rows contain two complete mode sets. Interpretation -------------- diff --git a/scripts/validate_qm9.py b/scripts/validate_qm9.py index e067707..3647cf8 100644 --- a/scripts/validate_qm9.py +++ b/scripts/validate_qm9.py @@ -1,5 +1,5 @@ #!/usr/bin/env python3 -"""Validate ThermoScreening against a deterministic sample of QM9.""" +"""Validate ThermoScreening against the QM9 dataset.""" from __future__ import annotations @@ -15,6 +15,8 @@ from ase import Atoms from ThermoScreening.thermo.api import run_thermo +from ThermoScreening.thermo.atoms import Atom +from ThermoScreening.thermo.system import dof FIGSHARE_ARTICLE_ID = "1057646" @@ -24,9 +26,10 @@ MOLECULE_COUNT = 133_885 SAMPLE_SIZE = 1_000 RANDOM_SEED = 20_260_724 +EDGE_CASE_IDS = {14_564} # QM9 does not store Gaussian's rotational symmetry number. Its published -# thermochemistry identifies these three exact small geometries as sigma=2; +# thermochemistry identifies these exact geometries as sigma=2; # the remaining records in the pinned sample use sigma=1. SYMMETRY_NUMBER_OVERRIDES = {3: 2, 4: 2, 23: 2} @@ -75,11 +78,18 @@ def _download_archive(cache_dir: Path) -> Path: def _sample_ids() -> set[int]: rng = np.random.default_rng(RANDOM_SEED) selected = set(range(1, 101)) + selected.update(EDGE_CASE_IDS) selected.update(np.linspace(101, MOLECULE_COUNT, 450, dtype=int).tolist()) remaining = np.asarray( - sorted(set(range(101, MOLECULE_COUNT + 1)) - selected) + sorted(set(range(1, MOLECULE_COUNT + 1)) - selected) + ) + selected.update( + rng.choice( + remaining, + size=SAMPLE_SIZE - len(selected), + replace=False, + ).tolist() ) - selected.update(rng.choice(remaining, size=450, replace=False).tolist()) if len(selected) != SAMPLE_SIZE: raise RuntimeError(f"Expected {SAMPLE_SIZE} sample identifiers.") return selected @@ -115,17 +125,33 @@ def _parse_record(text: str): return int(properties[1]), atoms, frequencies, expected -def _calculate_errors(archive: Path) -> dict[str, np.ndarray]: +def _source_dof(atom_count: int, frequency_count: int) -> tuple[int, bool]: + nonlinear = 3 * atom_count - 6 + linear = 3 * atom_count - 5 + if frequency_count in (nonlinear, 2 * nonlinear): + return nonlinear, frequency_count == 2 * nonlinear + if frequency_count in (linear, 2 * linear): + return linear, frequency_count == 2 * linear + raise RuntimeError( + f"Unexpected frequency count {frequency_count} for {atom_count} atoms." + ) + + +def _calculate_errors( + archive: Path, +) -> tuple[dict[str, np.ndarray], int]: selected = _sample_ids() errors = {name: [] for name in LIMITS} processed = set() + scanned = 0 + duplicated_frequency_rows = 0 with tarfile.open(archive, mode="r:bz2") as tar: for member in tar: stem = ( member.name.removeprefix("dsgdb9nsd_").removesuffix(".xyz") ) - if not stem.isdigit() or int(stem) not in selected: + if not stem.isdigit(): continue handle = tar.extractfile(member) @@ -137,6 +163,26 @@ def _calculate_errors(archive: Path) -> dict[str, np.ndarray]: if index != int(stem): raise RuntimeError(f"Index mismatch in {member.name}.") + source_dof, duplicated = _source_dof( + len(atoms), len(frequencies) + ) + actual_dof = dof( + [ + Atom(symbol=atom.symbol, position=atom.position) + for atom in atoms + ] + ) + if actual_dof != source_dof: + raise RuntimeError( + f"Geometry/mode mismatch for QM9 record {index}: " + f"calculated {actual_dof}, source {source_dof}." + ) + scanned += 1 + duplicated_frequency_rows += int(duplicated) + + if index not in selected: + continue + thermo = run_thermo( frequencies, atoms=atoms, @@ -161,7 +207,14 @@ def _calculate_errors(archive: Path) -> dict[str, np.ndarray]: missing = selected - processed if missing: raise RuntimeError(f"Missing {len(missing)} sampled QM9 records.") - return {name: np.asarray(values) for name, values in errors.items()} + if scanned != MOLECULE_COUNT: + raise RuntimeError( + f"Expected {MOLECULE_COUNT} QM9 records, scanned {scanned}." + ) + return ( + {name: np.asarray(values) for name, values in errors.items()}, + duplicated_frequency_rows, + ) def _default_cache_dir() -> Path: @@ -181,9 +234,11 @@ def main() -> int: args = parser.parse_args() archive = _download_archive(args.cache_dir.expanduser().resolve()) - errors = _calculate_errors(archive) + errors, duplicated_frequency_rows = _calculate_errors(archive) print(f"Figshare article: {FIGSHARE_ARTICLE_ID}") + print(f"QM9 geometries scanned: {MOLECULE_COUNT}") + print(f"Duplicated frequency rows: {duplicated_frequency_rows}") print(f"QM9 molecules sampled: {SAMPLE_SIZE}") passed = True for name, values in errors.items(): diff --git a/tests/thermo/test_system.py b/tests/thermo/test_system.py index a15dcb8..d3fefef 100644 --- a/tests/thermo/test_system.py +++ b/tests/thermo/test_system.py @@ -304,18 +304,53 @@ def test_linearity_detects_off_axis_linear_molecule(): assert linearity(_atoms_at(hcn, np.zeros(3))) is True -def test_linearity_rejects_slightly_bent_near_linear_geometry(): +@pytest.mark.parametrize( + "spec,expected_dof", + [ + ( + [ + ("N", (0.0174573422, -1.1613421749, -0.0041534236)), + ("C", (0.0025321236, -0.0034427793, 0.0017993915)), + ("C", (-0.0161141932, 1.3722098653, 0.0093487294)), + ("N", (-0.0326430226, 2.5300827589, 0.0160909027)), + ], + 6, + ), + ( + [ + ("N", (0.0871074894, -6.2768885493, -0.0324934680)), + ("C", (0.0711275515, -5.1137528555, -0.0261250526)), + ("C", (0.0525342088, -3.7552163012, -0.0186749511)), + ("C", (0.0360166993, -2.5388252685, -0.0119835707)), + ("C", (0.0177779288, -1.1939113440, -0.0045786311)), + ("C", (0.0013190793, 0.0224804809, 0.0021258345)), + ("C", (-0.0169207904, 1.3810213128, 0.0096207384)), + ("N", (-0.0325296666, 2.5441627248, 0.0160279504)), + ], + 18, + ), + ], + ids=["qm9-25", "qm9-14564"], +) +def test_linearity_rejects_slightly_bent_near_linear_geometry( + spec, expected_dof +): + atoms = _atoms(spec) + + assert linearity(atoms) is False + assert dof(atoms) == expected_dof + + +def test_linearity_tolerates_small_coordinate_noise(): atoms = _atoms( [ - ("N", (0.0174573422, -1.1613421749, -0.0041534236)), - ("C", (0.0025321236, -0.0034427793, 0.0017993915)), - ("C", (-0.0161141932, 1.3722098653, 0.0093487294)), - ("N", (-0.0326430226, 2.5300827589, 0.0160909027)), + ("O", (1e-6, 0.0, -1.16)), + ("C", (0.0, -1e-6, 0.0)), + ("O", (-1e-6, 0.0, 1.16)), ] ) - assert linearity(atoms) is False - assert dof(atoms) == 6 + assert linearity(atoms) is True @pytest.mark.parametrize(