Fix extreme near-linear classification

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:

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
 --------------

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():

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(