From c6d71ab78873cc3860ec671c4aa4366c32e0eb81 Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Thu, 23 Jul 2026 22:59:27 +0200 Subject: [PATCH 1/2] Add anthraquinone reference benchmark --- README.md | 12 ++ ThermoScreening/thermo/api.py | 37 ++++ ThermoScreening/thermo/system.py | 20 +- docs/benchmarks/anthraquinone_workflow.rst | 72 +++++-- scripts/validate_anthraquinone_zenodo.py | 240 +++++++++++++++++++++ tests/thermo/test_api.py | 25 +++ tests/thermo/test_system.py | 34 +++ 7 files changed, 426 insertions(+), 14 deletions(-) create mode 100644 scripts/validate_anthraquinone_zenodo.py diff --git a/README.md b/README.md index 02aeb25..27199ef 100644 --- a/README.md +++ b/README.md @@ -138,6 +138,18 @@ Run the full test suite: python -m pytest -q ``` +Reproduce the published DFTB/3ob/COSMO anthraquinone potentials from the raw +Zenodo data: + +```bash +python scripts/validate_anthraquinone_zenodo.py +``` + +The command verifies the downloaded archives, processes all 39 charge-state +calculations, and compares 25 calibrated potentials. See the +[anthraquinone reference benchmark](docs/benchmarks/anthraquinone_workflow.rst) +for the data source, historical settings, and interpretation. + Run linting: ```bash diff --git a/ThermoScreening/thermo/api.py b/ThermoScreening/thermo/api.py index 318ccfe..37c9619 100644 --- a/ThermoScreening/thermo/api.py +++ b/ThermoScreening/thermo/api.py @@ -350,6 +350,7 @@ def run_thermo( charge=0.0, atoms=None, spin=None, + symmetry_number=None, quasi_rrho=False, transition_state=False, ): @@ -375,6 +376,11 @@ def run_thermo( The system charge. atoms : ase.Atoms, optional An optimized geometry to use directly instead of reading ``coord_file``. + spin : float, optional + Spin quantum number S. Defaults to the minimum-spin electron-count guess. + symmetry_number : int, optional + Rotational symmetry number used in the rotational entropy. Defaults to + automatic point-group analysis. quasi_rrho : bool If True, use Grimme's quasi-RRHO treatment for the vibrational entropy instead of the pure harmonic oscillator. Default False. @@ -424,6 +430,7 @@ def run_thermo( pbc=pbc, charge=charge, spin=spin, + symmetry_number=symmetry_number, ) thermo_setup = Thermo( @@ -447,6 +454,7 @@ def orca_thermo( pressure=101325, charge=0.0, spin=None, + symmetry_number=None, quasi_rrho=False, transition_state=False, ): @@ -475,6 +483,9 @@ def orca_thermo( System charge. Default 0.0. spin : float, optional Spin quantum number S. Defaults to the minimum-spin electron-count guess. + symmetry_number : int, optional + Rotational symmetry number used in the rotational entropy. Defaults to + automatic point-group analysis. quasi_rrho : bool If True, use Grimme's quasi-RRHO vibrational entropy. Default False. transition_state : bool @@ -514,6 +525,7 @@ def orca_thermo( pressure=pressure, charge=charge, spin=spin, + symmetry_number=symmetry_number, engine="dftb+", quasi_rrho=quasi_rrho, transition_state=transition_state, @@ -527,6 +539,7 @@ def cclib_thermo( pressure=101325, charge=0.0, spin=None, + symmetry_number=None, quasi_rrho=False, transition_state=False, ): @@ -557,6 +570,9 @@ def cclib_thermo( System charge. Default 0.0. spin : float, optional Spin quantum number S. Defaults to the minimum-spin electron-count guess. + symmetry_number : int, optional + Rotational symmetry number used in the rotational entropy. Defaults to + automatic point-group analysis. quasi_rrho : bool If True, use Grimme's quasi-RRHO vibrational entropy. Default False. transition_state : bool @@ -593,6 +609,7 @@ def cclib_thermo( pressure=pressure, charge=charge, spin=spin, + symmetry_number=symmetry_number, engine="dftb+", quasi_rrho=quasi_rrho, transition_state=transition_state, @@ -626,6 +643,7 @@ def pyscf_thermo( pressure=101325, charge=0.0, spin=None, + symmetry_number=None, quasi_rrho=False, transition_state=False, ): @@ -657,6 +675,9 @@ def pyscf_thermo( System charge. Default 0.0. spin : float, optional Spin quantum number S. Defaults to the minimum-spin electron-count guess. + symmetry_number : int, optional + Rotational symmetry number used in the rotational entropy. Defaults to + automatic point-group analysis. quasi_rrho : bool If True, use Grimme's quasi-RRHO vibrational entropy. Default False. transition_state : bool @@ -704,6 +725,7 @@ def pyscf_thermo( pressure=pressure, charge=charge, spin=spin, + symmetry_number=symmetry_number, engine="dftb+", quasi_rrho=quasi_rrho, transition_state=transition_state, @@ -796,6 +818,7 @@ def dftbplus_thermo( directory=None, parameter_set="3ob", spin=None, + symmetry_number=None, spin_constants=None, solvent=None, solvation_param_file=None, @@ -829,6 +852,9 @@ def dftbplus_thermo( (even -> 0, odd -> 0.5). When S > 0 the DFTB+ steps run colinear spin-polarised (so radicals are treated open-shell automatically); S = 0 keeps the restricted closed-shell calculation. + symmetry_number : int, optional + Rotational symmetry number used in the rotational entropy. Defaults to + automatic point-group analysis. spin_constants : dict, optional Element -> spin-constant mapping matching the Slater-Koster set in use. Defaults to the 3ob constants; pass ``SPIN_CONSTANTS_MIO`` (or the value @@ -922,6 +948,7 @@ def dftbplus_thermo( engine='dftb+', charge=charge, spin=spin, + symmetry_number=symmetry_number, quasi_rrho=quasi_rrho, ) @@ -934,6 +961,7 @@ def xtb_thermo( pressure=101325, charge=0.0, spin=None, + symmetry_number=None, method="GFN2-xTB", directory=None, quasi_rrho=False, @@ -961,6 +989,9 @@ def xtb_thermo( Spin quantum number S. Defaults to the minimum-spin electron-count guess (even -> 0, odd -> 0.5). ``round(2*S)`` unpaired electrons are passed to xTB, so radicals run open-shell automatically. + symmetry_number : int, optional + Rotational symmetry number used in the rotational entropy. Defaults to + automatic point-group analysis. method : str GFN-xTB parametrisation, ``"GFN2-xTB"`` (default) or ``"GFN1-xTB"``. directory : str, optional @@ -1012,6 +1043,7 @@ def xtb_thermo( engine="xtb", charge=charge, spin=spin, + symmetry_number=symmetry_number, quasi_rrho=quasi_rrho, ) @@ -1024,6 +1056,7 @@ def xtb_cli_thermo( pressure=101325, charge=0.0, spin=None, + symmetry_number=None, method="GFN2-xTB", solvent=None, directory=None, @@ -1052,6 +1085,9 @@ def xtb_cli_thermo( Spin quantum number S. Defaults to the minimum-spin electron-count guess (even -> 0, odd -> 0.5); ``round(2*S)`` unpaired electrons are passed to xtb, so radicals run open-shell automatically. + symmetry_number : int, optional + Rotational symmetry number used in the rotational entropy. Defaults to + automatic point-group analysis. method : str GFN parametrisation: ``"GFN2-xTB"`` (default), ``"GFN1-xTB"`` or ``"GFN0-xTB"``. @@ -1088,6 +1124,7 @@ def xtb_cli_thermo( engine="xtb", charge=charge, spin=spin, + symmetry_number=symmetry_number, quasi_rrho=quasi_rrho, ) diff --git a/ThermoScreening/thermo/system.py b/ThermoScreening/thermo/system.py index 80196e0..49d13ac 100644 --- a/ThermoScreening/thermo/system.py +++ b/ThermoScreening/thermo/system.py @@ -602,6 +602,7 @@ def __init__( electronic_energy: float | None = None, vibrational_frequencies: np.ndarray | None = None, spin: float | None = None, + symmetry_number: int | None = None, ) -> None: """ Initializes the System with the given parameters. @@ -626,6 +627,12 @@ def __init__( The electronic energy of the system. vibrational_frequencies : np.ndarray, optional, default=None The vibrational frequencies of the system. + spin : float, optional, default=None + Spin quantum number. Defaults to the minimum-spin electron-count + guess. + symmetry_number : int, optional, default=None + Rotational symmetry number used in the rotational entropy. Defaults + to automatic point-group analysis. Raises ------ @@ -644,6 +651,13 @@ def __init__( if vibrational_frequencies is None: raise TSValueError("Vibrational frequencies must be provided.") + if symmetry_number is not None and ( + isinstance(symmetry_number, bool) or symmetry_number < 1 + ): + raise TSValueError( + "The rotational symmetry number must be a positive integer." + ) + self._atoms = atoms self._charge = charge self._periodicity = periodicity @@ -659,7 +673,11 @@ def __init__( self._dim = dim(atoms) self._mass = mass(atoms) self._center_of_mass = center_of_mass(atoms, self._mass) - self._symmetry_number = rotational_symmetry_number(atoms) + self._symmetry_number = ( + rotational_symmetry_number(atoms) + if symmetry_number is None + else symmetry_number + ) self._rotational_group = rotational_group_calc(atoms) if isinstance(self._cell, Cell): self._spacegroup_number = spacegroup_number(atoms, self._cell) diff --git a/docs/benchmarks/anthraquinone_workflow.rst b/docs/benchmarks/anthraquinone_workflow.rst index beab1ec..1c93703 100644 --- a/docs/benchmarks/anthraquinone_workflow.rst +++ b/docs/benchmarks/anthraquinone_workflow.rst @@ -1,17 +1,63 @@ -Anthraquinone workflow regression +Anthraquinone reference benchmark ================================= -The test suite contains a compact three-molecule anthraquinone fixture covering -the parent compound and two hydroxy substitution positions. It validates: +ThermoScreening can reproduce the DFTB/3ob/COSMO redox-potential tables from +Kröll *et al.*, `Phys. Chem. Chem. Phys. 24, 2022 +`_. The raw calculations and processed +tables are published in the immutable `Zenodo record 20796838 +`_. -- deterministic SMILES embedding and molecular formulas; -- construction of the oxidized, singly reduced and doubly reduced states; -- separate reference calibration of both reduction steps; -- potential-inversion classification; and -- reproducible input-set provenance. +Run the benchmark from a source checkout: -The fixture tests workflow invariants with controlled state energies. It is not -an absolute-potential accuracy benchmark and does not copy results from an -external dataset. Quantitative comparisons require the same structures, -electronic-structure method, solvent model, reference electrode and standard -states. +.. code-block:: bash + + python scripts/validate_anthraquinone_zenodo.py + +The command: + +- downloads the raw DFTB/3ob/COSMO-DMF outputs and published plot tables; +- verifies both archives against pinned SHA-256 checksums; +- processes 13 molecules in three charge states with :func:`run_thermo`; +- uses the archived analysis settings, including ``spin=1`` for monoanions and + ``symmetry_number=1`` because the historical workflow disabled symmetry + analysis; +- calibrates the parent anthraquinone result against each published table; and +- compares 13 first-reduction and 12 overall two-electron potentials. + +Reference result +---------------- + +.. list-table:: + :header-rows: 1 + + * - Metric + - Result + - Acceptance limit + * - Mean absolute error + - 0.021320 mV + - 0.05 mV + * - Maximum absolute error + - 0.334662 mV + - 0.5 mV + +The largest difference is the first reduction of molecule 7 +(1,4-diaminoanthraquinone). The remaining discrepancy is below the precision +relevant to the published millivolt-scale table. + +Interpretation +-------------- + +This benchmark verifies that ThermoScreening reproduces the archived +thermochemistry and redox post-processing when given the same raw electronic +energies, frequencies, geometries, spin convention, symmetry number, and +reference calibration. It does not establish experimental predictive accuracy. + +The archived ``DFTB_2e`` column is the calibrated neutral-to-dianion +two-electron average. It is reproduced with ``n_electrons=2`` and must not be +interpreted as the stepwise anion-to-dianion ``E2`` reported by the current +three-state redox workflow. The current workflow calibrates ``E1`` and ``E2`` +separately and reports ``E2e`` as their arithmetic mean. + +The ordinary test suite also retains a compact synthetic anthraquinone fixture. +That fixture checks workflow invariants without network access; the Zenodo +command is the quantitative external-data benchmark. diff --git a/scripts/validate_anthraquinone_zenodo.py b/scripts/validate_anthraquinone_zenodo.py new file mode 100644 index 0000000..04f9ab5 --- /dev/null +++ b/scripts/validate_anthraquinone_zenodo.py @@ -0,0 +1,240 @@ +#!/usr/bin/env python3 +"""Reproduce the published DFTB/3ob/COSMO anthraquinone potentials.""" + +from __future__ import annotations + +import argparse +import csv +import hashlib +import os +import shutil +import tarfile +import urllib.request +from pathlib import Path + +import numpy as np + +from ThermoScreening.thermo.api import read_vibrational, run_thermo +from ThermoScreening.thermo.reactions import ( + calibrate_reduction_reference, + reduction_potential, +) + + +RECORD_ID = "20796838" +BASE_URL = f"https://zenodo.org/records/{RECORD_ID}/files" +ARCHIVES = { + "04_dftb_3ob_cosmo_raw_outputs.tar.xz": ( + "c0ef53aa96547e10cf8478192abcd0dc0d09a0464ff41877016aa9fe461436dd" + ), + "08_processed_tables.tar.xz": ( + "4a2a083fd1c702e882844bfca210848b7dcf8afbec09ffb87e4e31db4121ab47" + ), +} +MAX_ERROR_MV = 0.5 +MAX_MEAN_ERROR_MV = 0.05 + +LABEL_TO_ORDER = { + "AQ": 1, + "1,2-OH": 2, + "1,4-OH": 3, + "1,5-OH": 4, + "1,8-OH": 5, + "1,2-NH2": 6, + "1,4-NH2": 7, + "2,6-NH2": 8, + "1-OH": 9, + "2-OH": 10, + "1-NH2": 11, + "2-NH2": 12, + "1-NH-4-OH": 13, + "1-NH2-4-OH": 13, +} + + +def _sha256(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as handle: + for block in iter(lambda: handle.read(1024 * 1024), b""): + digest.update(block) + return digest.hexdigest() + + +def _download_archive(cache_dir: Path, name: str, checksum: str) -> Path: + destination = cache_dir / name + if destination.exists() and _sha256(destination) == checksum: + return destination + + temporary = destination.with_suffix(destination.suffix + ".part") + request = urllib.request.Request( + f"{BASE_URL}/{name}?download=1", + headers={"User-Agent": "ThermoScreening-reference-validation"}, + ) + try: + with ( + urllib.request.urlopen(request, timeout=60) as response, + temporary.open("wb") as output, + ): + shutil.copyfileobj(response, output) + if _sha256(temporary) != checksum: + raise RuntimeError(f"Checksum mismatch for {name}.") + temporary.replace(destination) + finally: + temporary.unlink(missing_ok=True) + return destination + + +def _prepare_data(cache_dir: Path) -> Path: + cache_dir.mkdir(parents=True, exist_ok=True) + for name, checksum in ARCHIVES.items(): + archive = _download_archive(cache_dir, name, checksum) + with tarfile.open(archive, mode="r:xz") as tar: + tar.extractall(cache_dir, filter="data") + return cache_dir + + +def _read_published_potentials(path: Path) -> dict[int, float]: + rows = {} + with path.open(newline="", encoding="utf-8") as handle: + for row in csv.reader(handle, delimiter=";"): + if not row or row[0].startswith("#"): + continue + try: + order = LABEL_TO_ORDER[row[0]] + except KeyError as exc: + raise RuntimeError( + f"Unknown molecule label {row[0]!r} in {path}." + ) from exc + if order in rows: + raise RuntimeError(f"Duplicate molecule {order} in {path}.") + rows[order] = float(row[1]) / 1000.0 + return rows + + +def _calculate_states(data_dir: Path): + calculations = ( + data_dir / "dftb_3ob_cosmo_dmf" / "calculations" + ) + states = {} + for order in range(1, 14): + matches = list(calculations.glob(f"{order:02d}_*")) + if len(matches) != 1: + raise RuntimeError( + f"Expected one calculation directory for molecule {order}, " + f"found {len(matches)}." + ) + for charge in (0, -1, -2): + thermo_dir = matches[0] / f"run_{charge}" / "thermo" + frequencies = read_vibrational( + str(thermo_dir / "frequency.txt"), "dftb+" + ) + energy = float( + (thermo_dir / "electronic_energy.txt").read_text(encoding="utf-8") + ) + states[order, charge] = run_thermo( + frequencies, + coord_file=str(thermo_dir / "geo_opt.xyz"), + energy=energy, + charge=charge, + spin=1.0 if charge == -1 else 0.0, + symmetry_number=1, + ) + return states + + +def _compare(data_dir: Path) -> tuple[np.ndarray, tuple[str, int, float]]: + table_dir = ( + data_dir + / "processed_results" + / "tables" + / "dft_and_dftb_paper_plot_tables" + / "dftb_legacy_tables_from_plot_workflow" + ) + published_e1 = _read_published_potentials( + table_dir / "calc_exp_1_DFTB_COSMO_3ob.csv" + ) + published_e2e = _read_published_potentials( + table_dir / "calc_exp_2_DFTB_COSMO_3ob.csv" + ) + if set(published_e1) != set(range(1, 14)): + raise RuntimeError("The first-reduction table does not contain 13 molecules.") + if set(published_e2e) != {*range(1, 12), 13}: + raise RuntimeError("The two-electron table does not contain 12 molecules.") + states = _calculate_states(data_dir) + + e1_reference = calibrate_reduction_reference( + states[1, 0], + states[1, -1], + experimental_potential=published_e1[1], + ) + e2e_reference = calibrate_reduction_reference( + states[1, 0], + states[1, -2], + experimental_potential=published_e2e[1], + n_electrons=2, + ) + + errors = [] + labels = [] + for order, published in published_e1.items(): + calculated = reduction_potential( + states[order, 0], + states[order, -1], + reference_potential=e1_reference, + ) + errors.append(calculated - published) + labels.append(("E1", order)) + for order, published in published_e2e.items(): + calculated = reduction_potential( + states[order, 0], + states[order, -2], + n_electrons=2, + reference_potential=e2e_reference, + ) + errors.append(calculated - published) + labels.append(("E2e", order)) + + errors_mv = np.asarray(errors) * 1000.0 + worst_index = int(np.argmax(np.abs(errors_mv))) + worst = (*labels[worst_index], float(errors_mv[worst_index])) + return errors_mv, worst + + +def _default_cache_dir() -> Path: + root = Path(os.environ.get("XDG_CACHE_HOME", Path.home() / ".cache")) + return root / "thermoscreening" / f"anthraquinone-{RECORD_ID}" + + +def main() -> int: + """Download the reference data, run the comparison, and report status.""" + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--cache-dir", + type=Path, + default=_default_cache_dir(), + help="Archive and extraction cache directory.", + ) + args = parser.parse_args() + + data_dir = _prepare_data(args.cache_dir.expanduser().resolve()) + errors_mv, worst = _compare(data_dir) + mean_error = float(np.mean(np.abs(errors_mv))) + max_error = float(np.max(np.abs(errors_mv))) + + print(f"Zenodo record: {RECORD_ID}") + print("Charge-state calculations: 39") + print(f"Published potentials: {len(errors_mv)}") + print(f"Mean absolute error: {mean_error:.6f} mV") + print(f"Maximum absolute error: {max_error:.6f} mV") + print( + f"Worst case: {worst[0]}, molecule {worst[1]}, " + f"signed error {worst[2]:+.6f} mV" + ) + + passed = max_error <= MAX_ERROR_MV and mean_error <= MAX_MEAN_ERROR_MV + print("Result: PASS" if passed else "Result: FAIL") + return 0 if passed else 1 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tests/thermo/test_api.py b/tests/thermo/test_api.py index badba10..48784dd 100644 --- a/tests/thermo/test_api.py +++ b/tests/thermo/test_api.py @@ -326,6 +326,31 @@ def test_run_thermo_ase_atoms_matches_file_path(self): from_file.total_gibbs_free_energy("H") ) + def test_run_thermo_uses_explicit_symmetry_number(self): + hydrogen = Atoms( + "H2", + positions=[[0.0, 0.0, 0.0], [0.0, 0.0, 0.74]], + ) + frequencies = np.array([4400.0]) + + automatic = run_thermo(frequencies, atoms=hydrogen, energy=-1.0) + explicit = run_thermo( + frequencies, + atoms=hydrogen, + energy=-1.0, + symmetry_number=1, + ) + + assert automatic._system.rotational_symmetry_number == 2 + assert explicit._system.rotational_symmetry_number == 1 + entropy_difference = ( + explicit.total_entropy("cal/(mol*K)") + - automatic.total_entropy("cal/(mol*K)") + ) + assert entropy_difference == pytest.approx( + 8.314462618 / 4.184 * np.log(2.0) + ) + def test_run_in_directory_isolates_and_restores(tmp_path): from ThermoScreening.thermo.api import _run_in_directory diff --git a/tests/thermo/test_system.py b/tests/thermo/test_system.py index 145bb3f..778aa54 100644 --- a/tests/thermo/test_system.py +++ b/tests/thermo/test_system.py @@ -340,6 +340,40 @@ def test_system_rejects_negative_spin(): ) +def test_system_accepts_explicit_symmetry_number(): + atoms = [ + Atom(symbol="H", position=np.array([0.0, 0.0, 0.0])), + Atom(symbol="H", position=np.array([0.0, 0.0, 0.74])), + ] + + system = System( + atoms, + charge=0, + electronic_energy=-1.0, + vibrational_frequencies=np.array([4400.0]), + symmetry_number=1, + ) + + assert system.rotational_symmetry_number == 1 + + +@pytest.mark.parametrize("symmetry_number", [0, -1, True]) +def test_system_rejects_invalid_symmetry_number(symmetry_number): + atoms = [ + Atom(symbol="H", position=np.array([0.0, 0.0, 0.0])), + Atom(symbol="H", position=np.array([0.0, 0.0, 0.74])), + ] + + with pytest.raises(TSValueError, match="positive integer"): + System( + atoms, + charge=0, + electronic_energy=-1.0, + vibrational_frequencies=np.array([4400.0]), + symmetry_number=symmetry_number, + ) + + def test_rotational_symmetry_number_accepts_property(monkeypatch): class FakePointGroupAnalyzer: get_rotational_symmetry_number = 7 From 478bca74bc264c680283afa7eee78d46f1f561f2 Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Fri, 24 Jul 2026 08:41:36 +0200 Subject: [PATCH 2/2] Add QM9 thermochemistry benchmark --- README.md | 14 +- ThermoScreening/thermo/system.py | 2 +- docs/benchmarks/qm9_thermochemistry.rst | 82 ++++++++++ docs/index.rst | 1 + scripts/validate_qm9.py | 202 ++++++++++++++++++++++++ tests/thermo/test_system.py | 14 ++ 6 files changed, 311 insertions(+), 4 deletions(-) create mode 100644 docs/benchmarks/qm9_thermochemistry.rst create mode 100644 scripts/validate_qm9.py diff --git a/README.md b/README.md index 27199ef..c8db4a1 100644 --- a/README.md +++ b/README.md @@ -145,10 +145,18 @@ Zenodo data: python scripts/validate_anthraquinone_zenodo.py ``` -The command verifies the downloaded archives, processes all 39 charge-state -calculations, and compares 25 calibrated potentials. See the +Validate ideal-gas thermochemistry against a deterministic 1,000-molecule +sample from QM9: + +```bash +python scripts/validate_qm9.py +``` + +Both commands verify their downloaded archives and keep the external data out +of the repository. See the [anthraquinone reference benchmark](docs/benchmarks/anthraquinone_workflow.rst) -for the data source, historical settings, and interpretation. +and [QM9 thermochemistry benchmark](docs/benchmarks/qm9_thermochemistry.rst) +for their scope and interpretation. Run linting: diff --git a/ThermoScreening/thermo/system.py b/ThermoScreening/thermo/system.py index 49d13ac..29b1d3e 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-6 * eigenvalues[-1]) + return bool(eigenvalues[0] <= 1e-8 * eigenvalues[-1]) def dimensionality(atoms: List[Atom]) -> int: diff --git a/docs/benchmarks/qm9_thermochemistry.rst b/docs/benchmarks/qm9_thermochemistry.rst new file mode 100644 index 0000000..3d07954 --- /dev/null +++ b/docs/benchmarks/qm9_thermochemistry.rst @@ -0,0 +1,82 @@ +QM9 thermochemistry benchmark +============================= + +The `QM9 dataset `_ +contains optimized geometries, harmonic frequencies, and thermochemical +properties for 133,885 neutral CHONF molecules. The associated +`Scientific Data paper `_ reports that +the calculations used B3LYP/6-31G(2df,p) and provides ``U0``, ``U``, ``H``, +``G``, and ``Cv`` at 298.15 K. + +Run the benchmark from a source checkout: + +.. code-block:: bash + + 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. + +Method +------ + +QM9 stores the zero-point energy and internal energy at 0 K, so the electronic +energy is reconstructed as ``U0 - ZPVE``. ThermoScreening then recomputes all +thermal quantities from the published geometry and frequencies. + +The XYZ records do not store the rotational symmetry number used by Gaussian. +The pinned source settings use ``sigma=2`` for records 3, 4, and 23 and +``sigma=1`` for the other sampled records. This avoids treating a symmetry +metadata difference as a thermochemistry error. + +Reference result +---------------- + +Energy differences are reported in microhartree. The source energies are +printed to six decimal places and ``Cv`` to three decimal places. + +.. list-table:: + :header-rows: 1 + + * - Property + - Mean absolute error + - Maximum absolute error + * - ``U0`` + - 0.246147 microhartree + - 0.497836 microhartree + * - ``U`` + - 0.401046 microhartree + - 1.336198 microhartree + * - ``H`` + - 0.404191 microhartree + - 1.382253 microhartree + * - ``G`` + - 1.622078 microhartree + - 3.211252 microhartree + * - ``Cv`` + - 0.000259 cal/(mol K) + - 0.000585 cal/(mol K) + +The maximum Gibbs-energy difference is approximately 0.0020 kcal/mol. The +benchmark therefore reproduces the published quantities to the precision +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. + +Interpretation +-------------- + +This is an independent numerical validation of the ideal-gas rigid-rotor +harmonic-oscillator implementation over varied molecular sizes and structures. +It does not measure the predictive accuracy of B3LYP or compare computed +thermochemistry with experiment. diff --git a/docs/index.rst b/docs/index.rst index 0ae2ead..a2dcc2a 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -109,6 +109,7 @@ Features api benchmarks/entropy_accuracy benchmarks/anthraquinone_workflow + benchmarks/qm9_thermochemistry Indices ------- diff --git a/scripts/validate_qm9.py b/scripts/validate_qm9.py new file mode 100644 index 0000000..e067707 --- /dev/null +++ b/scripts/validate_qm9.py @@ -0,0 +1,202 @@ +#!/usr/bin/env python3 +"""Validate ThermoScreening against a deterministic sample of QM9.""" + +from __future__ import annotations + +import argparse +import hashlib +import os +import shutil +import tarfile +import urllib.request +from pathlib import Path + +import numpy as np +from ase import Atoms + +from ThermoScreening.thermo.api import run_thermo + + +FIGSHARE_ARTICLE_ID = "1057646" +ARCHIVE_NAME = "dsgdb9nsd.xyz.tar.bz2" +DOWNLOAD_URL = "https://ndownloader.figshare.com/files/3195389" +ARCHIVE_MD5 = "ad1ebd51ee7f5b3a6e32e974e5d54012" +MOLECULE_COUNT = 133_885 +SAMPLE_SIZE = 1_000 +RANDOM_SEED = 20_260_724 + +# QM9 does not store Gaussian's rotational symmetry number. Its published +# thermochemistry identifies these three exact small geometries as sigma=2; +# the remaining records in the pinned sample use sigma=1. +SYMMETRY_NUMBER_OVERRIDES = {3: 2, 4: 2, 23: 2} + +LIMITS = { + "U0": 0.6, + "U": 1.5, + "H": 1.5, + "G": 3.5, + "Cv": 0.0007, +} + + +def _md5(path: Path) -> str: + digest = hashlib.md5(usedforsecurity=False) + with path.open("rb") as handle: + for block in iter(lambda: handle.read(1024 * 1024), b""): + digest.update(block) + return digest.hexdigest() + + +def _download_archive(cache_dir: Path) -> Path: + cache_dir.mkdir(parents=True, exist_ok=True) + destination = cache_dir / ARCHIVE_NAME + if destination.exists() and _md5(destination) == ARCHIVE_MD5: + return destination + + temporary = destination.with_suffix(destination.suffix + ".part") + request = urllib.request.Request( + DOWNLOAD_URL, + headers={"User-Agent": "ThermoScreening-reference-validation"}, + ) + try: + with ( + urllib.request.urlopen(request, timeout=120) as response, + temporary.open("wb") as output, + ): + shutil.copyfileobj(response, output) + if _md5(temporary) != ARCHIVE_MD5: + raise RuntimeError(f"Checksum mismatch for {ARCHIVE_NAME}.") + temporary.replace(destination) + finally: + temporary.unlink(missing_ok=True) + return destination + + +def _sample_ids() -> set[int]: + rng = np.random.default_rng(RANDOM_SEED) + selected = set(range(1, 101)) + selected.update(np.linspace(101, MOLECULE_COUNT, 450, dtype=int).tolist()) + remaining = np.asarray( + sorted(set(range(101, MOLECULE_COUNT + 1)) - selected) + ) + 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 + + +def _number(value: str) -> float: + return float(value.replace("*^", "e")) + + +def _parse_record(text: str): + lines = text.splitlines() + atom_count = int(lines[0]) + properties = lines[1].split() + atom_rows = [line.split() for line in lines[2 : 2 + atom_count]] + atoms = Atoms( + [row[0] for row in atom_rows], + positions=[ + [_number(value) for value in row[1:4]] + for row in atom_rows + ], + ) + frequencies = np.asarray( + [_number(value) for value in lines[2 + atom_count].split()] + ) + expected = { + "zpve": _number(properties[11]), + "U0": _number(properties[12]), + "U": _number(properties[13]), + "H": _number(properties[14]), + "G": _number(properties[15]), + "Cv": _number(properties[16]), + } + return int(properties[1]), atoms, frequencies, expected + + +def _calculate_errors(archive: Path) -> dict[str, np.ndarray]: + selected = _sample_ids() + errors = {name: [] for name in LIMITS} + processed = set() + + 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: + continue + + handle = tar.extractfile(member) + if handle is None: + raise RuntimeError(f"Cannot read {member.name}.") + index, atoms, frequencies, expected = _parse_record( + handle.read().decode("utf-8") + ) + if index != int(stem): + raise RuntimeError(f"Index mismatch in {member.name}.") + + thermo = run_thermo( + frequencies, + atoms=atoms, + energy=expected["U0"] - expected["zpve"], + temperature=298.15, + pressure=101325, + spin=0.0, + symmetry_number=SYMMETRY_NUMBER_OVERRIDES.get(index, 1), + ) + actual = { + "U0": thermo.total_EeZPE(), + "U": thermo.total_EeEtot(), + "H": thermo.total_EeHtot(), + "G": thermo.total_EeGtot(), + "Cv": thermo.total_heat_capacity("cal/(mol*K)"), + } + for name in LIMITS: + factor = 1.0 if name == "Cv" else 1e6 + errors[name].append((actual[name] - expected[name]) * factor) + processed.add(index) + + 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()} + + +def _default_cache_dir() -> Path: + root = Path(os.environ.get("XDG_CACHE_HOME", Path.home() / ".cache")) + return root / "thermoscreening" / "qm9" + + +def main() -> int: + """Download QM9, run the comparison, and report status.""" + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--cache-dir", + type=Path, + default=_default_cache_dir(), + help="Archive cache directory.", + ) + args = parser.parse_args() + + archive = _download_archive(args.cache_dir.expanduser().resolve()) + errors = _calculate_errors(archive) + + print(f"Figshare article: {FIGSHARE_ARTICLE_ID}") + print(f"QM9 molecules sampled: {SAMPLE_SIZE}") + passed = True + for name, values in errors.items(): + absolute = np.abs(values) + unit = "cal/(mol*K)" if name == "Cv" else "microhartree" + maximum = float(np.max(absolute)) + mean = float(np.mean(absolute)) + print(f"{name}: MAE {mean:.6f}, maximum {maximum:.6f} {unit}") + passed = passed and maximum <= LIMITS[name] + + print("Result: PASS" if passed else "Result: FAIL") + return 0 if passed else 1 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tests/thermo/test_system.py b/tests/thermo/test_system.py index 778aa54..a15dcb8 100644 --- a/tests/thermo/test_system.py +++ b/tests/thermo/test_system.py @@ -304,6 +304,20 @@ 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(): + 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)), + ] + ) + + assert linearity(atoms) is False + assert dof(atoms) == 6 + + @pytest.mark.parametrize( "symbols,charge,expected", [