diff --git a/README.md b/README.md index 3cc5569..1bbcb40 100644 --- a/README.md +++ b/README.md @@ -159,12 +159,27 @@ the PSI4 VQM24 dataset: python scripts/validate_vqm24.py ``` +Validate zero-point energies across 255 model chemistries: + +```bash +python scripts/validate_model_chemistry_zpe.py +``` + +Validate complete DFTB+ and GFN2-xTB entropy workflows against an experimental +gas-phase compilation: + +```bash +python scripts/validate_experimental_entropy.py --engine both +``` + These commands verify their downloaded archives and keep the external data out of the repository. See the [anthraquinone reference benchmark](docs/benchmarks/anthraquinone_workflow.rst) and the thermochemistry benchmarks for [QM9](docs/benchmarks/qm9_thermochemistry.rst) and -[VQM24 transition states](docs/benchmarks/vqm24_transition_states.rst) for +[VQM24 transition states](docs/benchmarks/vqm24_transition_states.rst), plus +the [experimental entropy](docs/benchmarks/entropy_accuracy.rst) and +[model-chemistry ZPE](docs/benchmarks/model_chemistry_zpe.rst) comparisons, for their scope and interpretation. Run linting: diff --git a/docs/benchmarks/entropy_accuracy.rst b/docs/benchmarks/entropy_accuracy.rst index 4c33e3c..23cc705 100644 --- a/docs/benchmarks/entropy_accuracy.rst +++ b/docs/benchmarks/entropy_accuracy.rst @@ -1,94 +1,87 @@ -Gas-phase entropy accuracy -========================== +Experimental gas-phase entropy +============================== -Standard molar entropies S°(298.15 K, 1 bar) from ThermoScreening compared with -experiment, for both engines. These are method accuracies (GFN2-xTB and -DFTB3/3ob), not tool errors. The thermochemistry implementation itself matches -ASE ``IdealGasThermo`` to within 0.0003 cal/(mol K). +ThermoScreening is checked against the experimental entropy compilation +published by Bains, Petkowski, Zhan, and Seager. The source combines 1,403 +reported entropy values in eight source columns, including NIST, JANAF, and +Yaws. -.. list-table:: - :header-rows: 1 +The benchmark downloads the original workbook from +`Zenodo record 4661783 `__, verifies its +checksum, and keeps it outside the repository. + +Scope +----- + +The comparison covers 29 neutral gas-phase molecules containing H, C, N, O, F, +P, S, and Cl. It includes linear and nonlinear molecules, triplet oxygen, +doublet nitric oxide, high-symmetry species, and flexible molecules with +low-frequency torsions. + +For each molecule, the reference value is the median of the available +experimental columns in the workbook. This avoids selecting a single +transcription when the compilation contains conflicting entries, as it does +for methane. Nitric oxide is the documented exception: its two entries differ +by almost 10 J/(mol K), so the JANAF value is used. Values are converted from +J/(mol K) to cal/(mol K). + +Both engines start from the corresponding ASE G2 geometry and run a complete +optimization, Hessian, symmetry analysis, and quasi-RRHO thermochemistry +calculation at 298.15 K and 1 bar. GFN2-xTB uses a force threshold of +0.001 eV/A because methyl torsions remain sensitive at the normal screening +tolerance. + +Run the benchmark +----------------- + +DFTB+ and the 3ob parameter set are required for the default run: + +.. code-block:: bash - * - Molecule - - Symmetry number - - xTB (GFN2) - - DFTB+ 3ob - - Experiment - * - H₂O - - 2 - - 45.08 - - 45.07 - - 45.1 - * - CH₄ - - 12 - - 44.44 - - 44.47 - - 44.5 - * - NH₃ - - 3 - - 45.96 - - 46.00 - - 46.0 - * - N₂ - - 2 - - 45.79 - - 45.76 - - 45.8 - * - CO₂ - - 2 - - 51.23 - - 51.53 - - 51.1 - * - CO - - 1 - - 47.23 - - 47.14 - - 47.2 - * - CH₃OH - - 1 - - 56.63 - - 57.00 - - 57.3 - * - C₂H₆ - - 6 - - 54.38 - - 54.54 - - 54.8 - -Values are in cal/(mol K). - -Deviation from experiment -------------------------- + python scripts/validate_experimental_entropy.py + +Run both independently implemented engines with: + +.. code-block:: bash + + python scripts/validate_experimental_entropy.py --engine both + +The xTB comparison requires the optional ``tblite`` package. The script uses a +separate working directory for every molecule and engine so fixed-name Hessian +and vibration files cannot leak between cases. + +Results +------- .. list-table:: :header-rows: 1 * - Engine - - MAD + - MAE - RMSD - Maximum absolute deviation - - Mean signed deviation - * - xTB (GFN2) - - 0.17 - - 0.28 - - 0.67 (CH₃OH) - - -0.13 - * - DFTB+ 3ob - - 0.14 - - 0.21 - - 0.43 (CO₂) - - -0.04 - -Both engines reproduce experimental gas-phase entropies to well under -1 cal/(mol K). - -Notes ------ + * - DFTB3/3ob + - 0.377 + - 0.635 + - 2.317 (ethanol) + * - GFN2-xTB + - 0.457 + - 0.777 + - 2.387 (ethanol) + +All deviations are in cal/(mol K). The regression gate requires an MAE no +greater than 0.75 and a maximum deviation no greater than 3.0 cal/(mol K). + +Interpretation +-------------- + +This is an end-to-end method benchmark: it exercises geometry optimization, +frequencies, automatic rotational symmetry numbers, open-shell handling, +quasi-RRHO, and entropy assembly. The remaining deviations include the +accuracy limits of DFTB3/3ob and GFN2-xTB and are not numerical errors in the +thermochemistry formulas. Those formulas are separately tested against ASE +``IdealGasThermo``. -- Symmetry numbers are detected correctly for all eight molecules. -- The values use ``pressure=100000`` (1 bar). The default pressure of 101325 Pa - raises the entropy by 0.026 cal/(mol K). -- Residual error for methanol and ethane comes from treating low-frequency - torsions as harmonic oscillators; ``quasi_rrho=True`` softens this. -- The Linux release gate runs a mandatory DFTB+ calculation. Broader xTB and - solvent checks remain optional integration tests. +The source workbook is not redistributed. See the +`dataset publication `__ for provenance, +experimental-source details, and limitations. diff --git a/docs/benchmarks/model_chemistry_zpe.rst b/docs/benchmarks/model_chemistry_zpe.rst new file mode 100644 index 0000000..28dde31 --- /dev/null +++ b/docs/benchmarks/model_chemistry_zpe.rst @@ -0,0 +1,62 @@ +Model-chemistry zero-point energies +=================================== + +The zero-point energy implementation is checked against the quantum-chemistry +outputs published with the model-chemistry benchmark by Green, Wu, Doner, and +Pang. The source spans density-functional, double-hybrid, perturbation, and +composite methods across many basis sets. + +The validation script downloads ``qm_outputs.json`` from +`Zenodo record 14733712 `__, verifies its +checksum, and keeps the 55 MB source file outside the repository. + +Run the benchmark +----------------- + +.. code-block:: bash + + python scripts/validate_model_chemistry_zpe.py + +No external electronic-structure engine is needed. ThermoScreening recomputes +the zero-point energy directly from each source frequency spectrum. + +Coverage +-------- + +The script scans all 41,492 frequency records across 255 model chemistries. +Eleven source rows have fewer modes than their geometries require; their exact +identities are pinned and excluded. The remaining 41,481 records must be +structurally valid. + +A deterministic 1,000-record sample covers every one of the 255 model +chemistries, evenly spaced records, and a seeded random remainder. The sample +contains molecules with H, C, N, O, F, S, Cl, and Br and up to 26 atoms. + +Results +------- + +.. list-table:: + :header-rows: 1 + + * - Quantity + - Mean absolute deviation + - Maximum absolute deviation + * - Zero-point energy + - 0.345971 microhartree + - 1.459415 microhartree + +The regression gate is 0.5 microhartree for the mean and 2.0 microhartree for +the maximum. The residual is consistent with constants and rounding in the +source outputs. + +Interpretation +-------------- + +This benchmark validates frequency counting, unit conversion, and zero-point +energy assembly over a much broader method and element matrix than QM9. It +does not assess whether a particular electronic-structure method predicts an +experimentally accurate frequency. + +See the +`dataset publication `__ +for the model-chemistry definitions and scientific scope. diff --git a/docs/index.rst b/docs/index.rst index 5da0853..a79b387 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -108,6 +108,7 @@ Features configuration api benchmarks/entropy_accuracy + benchmarks/model_chemistry_zpe benchmarks/anthraquinone_workflow benchmarks/qm9_thermochemistry benchmarks/vqm24_transition_states diff --git a/scripts/validate_experimental_entropy.py b/scripts/validate_experimental_entropy.py new file mode 100644 index 0000000..042773c --- /dev/null +++ b/scripts/validate_experimental_entropy.py @@ -0,0 +1,335 @@ +#!/usr/bin/env python3 +"""Validate gas-phase entropies against an experimental compilation.""" + +from __future__ import annotations + +import argparse +import hashlib +import os +import posixpath +import re +import shutil +import statistics +import urllib.request +import zipfile +from pathlib import Path +from xml.etree import ElementTree + +import numpy as np +from ase.build import molecule + +from ThermoScreening.calculator.dftbplus import dftb_3ob_parameters +from ThermoScreening.thermo.api import dftbplus_thermo, xtb_thermo + + +RECORD_ID = "4661783" +ARCHIVE_NAME = "Measured_entropy_V2.1.xlsx" +DOWNLOAD_URL = ( + f"https://zenodo.org/records/{RECORD_ID}/files/{ARCHIVE_NAME}?download=1" +) +ARCHIVE_MD5 = "aa275812f3b010c9474ab8804154600b" +SHEET_NAME = "Measured Entropy" +TEMPERATURE = 298.15 +PRESSURE = 100_000.0 +JOULE_PER_CALORIE = 4.184 +MAE_LIMIT = 0.75 +MAXIMUM_LIMIT = 3.0 + +CASES = ( + ("H2", "hydrogen", 0.0), + ("CH4", "methane", 0.0), + ("H2O", "water", 0.0), + ("N2", "nitrogen", 0.0), + ("CO", "carbon monoxide", 0.0), + ("C2H2", "ethyne", 0.0), + ("HCN", "hydrogen cyanide", 0.0), + ("O2", "oxygen", 1.0), + ("NO", "nitric oxide", 0.5), + ("SH2", "hydrogen sulfide", 0.0), + ("PH3", "phosphane", 0.0), + ("NH3", "ammonia", 0.0), + ("CO2", "carbon dioxide", 0.0), + ("H2CO", "formaldehyde", 0.0), + ("C2H4", "ethene", 0.0), + ("OCS", "carbon oxide sulfide", 0.0), + ("CS2", "carbon disulfide", 0.0), + ("SO2", "sulfur dioxide", 0.0), + ("HCOOH", "formic acid", 0.0), + ("CF4", "perfluoromethane", 0.0), + ("CH3Cl", "chloromethane", 0.0), + ("H2CCl2", "dichloromethane", 0.0), + ("CCl4", "perchloromethane", 0.0), + ("CH3OH", "methanol", 0.0), + ("C2H6", "ethane", 0.0), + ("CH3CHO", "acetaldehyde", 0.0), + ("CH3CH2OH", "ethanol", 0.0), + ("CH3COCH3", "acetone", 0.0), + ("C6H6", "benzene", 0.0), +) + +REFERENCE_COLUMN_OVERRIDES = {"nitric oxide": "JANAF"} + +MAIN_NS = "http://schemas.openxmlformats.org/spreadsheetml/2006/main" +REL_NS = "http://schemas.openxmlformats.org/package/2006/relationships" +OFFICE_REL_NS = ( + "http://schemas.openxmlformats.org/officeDocument/2006/relationships" +) + + +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 _shared_strings(workbook: zipfile.ZipFile) -> list[str]: + try: + root = ElementTree.fromstring(workbook.read("xl/sharedStrings.xml")) + except KeyError: + return [] + return [ + "".join(node.text or "" for node in item.iter(f"{{{MAIN_NS}}}t")) + for item in root.findall(f"{{{MAIN_NS}}}si") + ] + + +def _sheet_path(workbook: zipfile.ZipFile, sheet_name: str) -> str: + root = ElementTree.fromstring(workbook.read("xl/workbook.xml")) + relationship_id = None + for sheet in root.findall(f".//{{{MAIN_NS}}}sheet"): + if sheet.attrib.get("name") == sheet_name: + relationship_id = sheet.attrib[ + f"{{{OFFICE_REL_NS}}}id" + ] + break + if relationship_id is None: + raise RuntimeError(f"Workbook has no {sheet_name!r} sheet.") + + relationships = ElementTree.fromstring( + workbook.read("xl/_rels/workbook.xml.rels") + ) + for relationship in relationships.findall(f"{{{REL_NS}}}Relationship"): + if relationship.attrib.get("Id") == relationship_id: + target = relationship.attrib["Target"].lstrip("/") + if target.startswith("xl/"): + return target + return posixpath.normpath(posixpath.join("xl", target)) + raise RuntimeError(f"Cannot resolve the {sheet_name!r} sheet.") + + +def _column_index(reference: str) -> int: + letters = re.match(r"[A-Z]+", reference) + if letters is None: + raise RuntimeError(f"Invalid cell reference {reference!r}.") + index = 0 + for letter in letters.group(): + index = index * 26 + ord(letter) - ord("A") + 1 + return index - 1 + + +def _cell_value(cell, shared_strings: list[str]): + value = cell.find(f"{{{MAIN_NS}}}v") + cell_type = cell.attrib.get("t") + if cell_type == "inlineStr": + return "".join( + node.text or "" + for node in cell.iter(f"{{{MAIN_NS}}}t") + ) + if value is None or value.text is None: + return None + if cell_type == "s": + return shared_strings[int(value.text)] + if cell_type in {"str", "e"}: + return value.text + return float(value.text) + + +def _worksheet_rows(path: Path) -> list[list[object | None]]: + with zipfile.ZipFile(path) as workbook: + shared_strings = _shared_strings(workbook) + sheet = ElementTree.fromstring( + workbook.read(_sheet_path(workbook, SHEET_NAME)) + ) + + rows = [] + for row in sheet.findall(f".//{{{MAIN_NS}}}row"): + values: list[object | None] = [] + for cell in row.findall(f"{{{MAIN_NS}}}c"): + index = _column_index(cell.attrib["r"]) + if index >= len(values): + values.extend([None] * (index + 1 - len(values))) + values[index] = _cell_value(cell, shared_strings) + rows.append(values) + return rows + + +def _experimental_entropies(path: Path) -> dict[str, float]: + rows = _worksheet_rows(path) + expected_header = [ + "Name", + "SMILES", + "NIST", + "NIST-12", + "Yaws", + "Sander et al", + "Sarangam", + "Mu and He", + "Rihani", + "JANAF", + ] + if rows[2][:10] != expected_header: + raise RuntimeError("Experimental entropy workbook columns changed.") + + source_rows = { + str(row[0]).strip().lower(): row + for row in rows[3:] + if row and row[0] is not None + } + entropies = {} + for _, source_name, _ in CASES: + row = source_rows.get(source_name) + if row is None: + raise RuntimeError(f"Missing experimental row {source_name!r}.") + values = [ + float(value) + for value in row[2:10] + if isinstance(value, (int, float)) + ] + if not values: + raise RuntimeError(f"No experimental values for {source_name!r}.") + if source_name in REFERENCE_COLUMN_OVERRIDES: + column = expected_header.index( + REFERENCE_COLUMN_OVERRIDES[source_name] + ) + reference = row[column] + if not isinstance(reference, (int, float)): + raise RuntimeError( + f"Missing {REFERENCE_COLUMN_OVERRIDES[source_name]} " + f"value for {source_name!r}." + ) + else: + reference = statistics.median(values) + entropies[source_name] = float(reference) / JOULE_PER_CALORIE + return entropies + + +def _run_engine( + engine: str, + references: dict[str, float], + work_dir: Path, +) -> tuple[np.ndarray, list[tuple[str, float]]]: + errors = [] + details = [] + for ase_name, source_name, spin in CASES: + job_dir = work_dir / engine.replace("+", "plus") / ase_name + shutil.rmtree(job_dir, ignore_errors=True) + job_dir.mkdir(parents=True) + common = { + "atoms": molecule(ase_name), + "temperature": TEMPERATURE, + "pressure": PRESSURE, + "spin": spin, + "directory": job_dir, + "quasi_rrho": True, + } + if engine == "dftb+": + thermo = dftbplus_thermo( + **common, + **dftb_3ob_parameters, + ) + else: + thermo = xtb_thermo( + **common, + method="GFN2-xTB", + fmax=0.001, + ) + actual = thermo.total_entropy("cal/(mol*K)") + error = actual - references[source_name] + errors.append(error) + details.append((ase_name, error)) + return np.asarray(errors), details + + +def _default_cache_dir() -> Path: + root = Path(os.environ.get("XDG_CACHE_HOME", Path.home() / ".cache")) + return root / "thermoscreening" / f"experimental-entropy-{RECORD_ID}" + + +def main() -> int: + """Download the dataset, run selected engines, and report status.""" + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--engine", + choices=("dftb+", "xtb", "both"), + default="dftb+", + help="Engine to validate (default: dftb+).", + ) + parser.add_argument( + "--cache-dir", + type=Path, + default=_default_cache_dir(), + help="Dataset and working-directory cache.", + ) + args = parser.parse_args() + + cache_dir = args.cache_dir.expanduser().resolve() + archive = _download_archive(cache_dir) + references = _experimental_entropies(archive) + engines = ("dftb+", "xtb") if args.engine == "both" else (args.engine,) + + print(f"Zenodo record: {RECORD_ID}") + print(f"Experimental molecules: {len(CASES)}") + print("Reference: median reported S; JANAF value for nitric oxide") + print("Treatment: quasi-RRHO; xTB optimization fmax=0.001 eV/A") + + passed = True + for engine in engines: + errors, details = _run_engine(engine, references, cache_dir / "work") + absolute = np.abs(errors) + mean = float(np.mean(absolute)) + rmsd = float(np.sqrt(np.mean(np.square(errors)))) + worst_index = int(np.argmax(absolute)) + worst_name, worst_error = details[worst_index] + print( + f"{engine}: MAE {mean:.3f}, RMSD {rmsd:.3f}, " + f"maximum {abs(worst_error):.3f} cal/(mol*K) ({worst_name})" + ) + passed = ( + passed + and mean <= MAE_LIMIT + and abs(worst_error) <= MAXIMUM_LIMIT + ) + + print("Result: PASS" if passed else "Result: FAIL") + return 0 if passed else 1 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/validate_model_chemistry_zpe.py b/scripts/validate_model_chemistry_zpe.py new file mode 100644 index 0000000..8bddc86 --- /dev/null +++ b/scripts/validate_model_chemistry_zpe.py @@ -0,0 +1,276 @@ +#!/usr/bin/env python3 +"""Validate zero-point energies against a broad model-chemistry dataset.""" + +from __future__ import annotations + +import argparse +import hashlib +import json +import os +import shutil +import urllib.request +from pathlib import Path + +import numpy as np +from ase import Atoms + +from ThermoScreening.thermo.api import run_thermo +from ThermoScreening.thermo.atoms import Atom +from ThermoScreening.thermo.system import resolved_dof + + +RECORD_ID = "14733712" +ARCHIVE_NAME = "qm_outputs.json" +DOWNLOAD_URL = ( + f"https://zenodo.org/records/{RECORD_ID}/files/{ARCHIVE_NAME}?download=1" +) +ARCHIVE_MD5 = "58e758a09e54cbaea7cc675c59b2ec06" +MODEL_CHEMISTRY_COUNT = 334 +FREQUENCY_MODEL_CHEMISTRY_COUNT = 255 +FREQUENCY_RECORD_COUNT = 41_492 +VALID_RECORD_COUNT = 41_481 +SAMPLE_SIZE = 1_000 +RANDOM_SEED = 20_260_724 +MAE_LIMIT_MICROHARTREE = 0.5 +MAX_LIMIT_MICROHARTREE = 2.0 + +EXPECTED_MODE_MISMATCHES = { + ("B2GP-PLYP_def2tzvpp", "241"), + ("M06-2X_cbsb7", "291"), + ("M06L_def2-tzvppd", "291"), + ("M06L_def2tzvpp", "286"), + ("M06L_def2tzvpp", "291"), + ("M06_def2-tzvppd", "291"), + ("PBE_def2tzvpp", "241"), + ("wB97X-D_6-31g-2df-p", "231"), + ("wB97X-D_aug-cc-pV-Q-d-Z", "177"), + ("wB97X-D_aug-cc-pV-T-d-Z", "231"), + ("wB97X-D_def2tzvpp", "291"), +} + + +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 _parse_geometry(text: str) -> Atoms: + rows = [line.split() for line in text.splitlines() if line.strip()] + return Atoms( + numbers=[int(row[0]) for row in rows], + positions=[ + [float(value) for value in row[1:4]] + for row in rows + ], + ) + + +def _atom_list(atoms: Atoms) -> list[Atom]: + return [ + Atom(symbol=atom.symbol, position=atom.position) + for atom in atoms + ] + + +def _species_sort_key(species_id: str) -> tuple[int, int | str]: + if species_id.isdigit(): + return 0, int(species_id) + return 1, species_id + + +def _scan_records( + data: dict, +) -> tuple[list[tuple[str, str, Atoms, np.ndarray, float]], set[tuple[str, str]]]: + if len(data) != MODEL_CHEMISTRY_COUNT: + raise RuntimeError( + f"Expected {MODEL_CHEMISTRY_COUNT} model chemistries, " + f"found {len(data)}." + ) + + records = [] + mismatches = set() + frequency_records = 0 + frequency_methods = 0 + + for method in sorted(data): + method_records = data[method] + has_frequencies = False + for species_id in sorted(method_records, key=_species_sort_key): + record = method_records[species_id] + if not isinstance(record, dict) or record.get("frequencies") is None: + continue + + has_frequencies = True + frequency_records += 1 + required = ("geometry", "energy", "frequencies", "zpe") + if any(record.get(field) is None for field in required): + raise RuntimeError( + f"Incomplete frequency record {method}/{species_id}." + ) + + atoms = _parse_geometry(record["geometry"]) + frequencies = np.asarray(record["frequencies"], dtype=float) + expected_dof = resolved_dof(_atom_list(atoms), len(frequencies)) + if len(frequencies) != expected_dof: + mismatches.add((method, species_id)) + continue + if np.any(frequencies <= 0): + raise RuntimeError( + f"Non-positive frequency in {method}/{species_id}." + ) + records.append( + ( + method, + species_id, + atoms, + frequencies, + float(record["zpe"]), + ) + ) + frequency_methods += int(has_frequencies) + + if frequency_methods != FREQUENCY_MODEL_CHEMISTRY_COUNT: + raise RuntimeError( + f"Expected {FREQUENCY_MODEL_CHEMISTRY_COUNT} frequency methods, " + f"found {frequency_methods}." + ) + if frequency_records != FREQUENCY_RECORD_COUNT: + raise RuntimeError( + f"Expected {FREQUENCY_RECORD_COUNT} frequency records, " + f"found {frequency_records}." + ) + if mismatches != EXPECTED_MODE_MISMATCHES: + raise RuntimeError( + "Model-chemistry geometry/mode mismatches changed: " + f"{sorted(mismatches)}." + ) + if len(records) != VALID_RECORD_COUNT: + raise RuntimeError( + f"Expected {VALID_RECORD_COUNT} valid records, found {len(records)}." + ) + return records, mismatches + + +def _sample_indices( + records: list[tuple[str, str, Atoms, np.ndarray, float]], +) -> list[int]: + selected = {} + for index, record in enumerate(records): + selected.setdefault(record[0], index) + + indices = set(selected.values()) + indices.update( + np.linspace(0, len(records) - 1, 450, dtype=int).tolist() + ) + remaining = np.asarray(sorted(set(range(len(records))) - indices)) + rng = np.random.default_rng(RANDOM_SEED) + indices.update( + rng.choice( + remaining, + size=SAMPLE_SIZE - len(indices), + replace=False, + ).tolist() + ) + if len(indices) != SAMPLE_SIZE: + raise RuntimeError(f"Expected {SAMPLE_SIZE} sampled records.") + if len({records[index][0] for index in indices}) != len(selected): + raise RuntimeError("The sample does not cover every frequency method.") + return sorted(indices) + + +def _calculate_errors( + records: list[tuple[str, str, Atoms, np.ndarray, float]], + indices: list[int], +) -> np.ndarray: + errors = [] + for index in indices: + _, _, atoms, frequencies, expected_zpe = records[index] + thermo = run_thermo( + frequencies, + atoms=atoms, + energy=0.0, + spin=0.0, + symmetry_number=1, + ) + errors.append((thermo.total_EeZPE() - expected_zpe) * 1e6) + return np.asarray(errors) + + +def _default_cache_dir() -> Path: + root = Path(os.environ.get("XDG_CACHE_HOME", Path.home() / ".cache")) + return root / "thermoscreening" / f"model-chemistry-{RECORD_ID}" + + +def main() -> int: + """Download the dataset, run the comparison, and report status.""" + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--cache-dir", + type=Path, + default=_default_cache_dir(), + help="Dataset cache directory.", + ) + args = parser.parse_args() + + archive = _download_archive(args.cache_dir.expanduser().resolve()) + with archive.open(encoding="utf-8") as handle: + data = json.load(handle) + + records, mismatches = _scan_records(data) + indices = _sample_indices(records) + errors = _calculate_errors(records, indices) + absolute = np.abs(errors) + mean = float(np.mean(absolute)) + maximum = float(np.max(absolute)) + + print(f"Zenodo record: {RECORD_ID}") + print(f"Model chemistries scanned: {MODEL_CHEMISTRY_COUNT}") + print( + "Model chemistries with frequencies: " + f"{FREQUENCY_MODEL_CHEMISTRY_COUNT}" + ) + print(f"Frequency records scanned: {FREQUENCY_RECORD_COUNT}") + print(f"Source mode-count mismatches excluded: {len(mismatches)}") + print(f"Valid records: {len(records)}") + print(f"Records sampled: {len(indices)}") + print(f"ZPE: MAE {mean:.6f} microhartree") + print(f"ZPE: maximum {maximum:.6f} microhartree") + + passed = ( + mean <= MAE_LIMIT_MICROHARTREE + and maximum <= MAX_LIMIT_MICROHARTREE + ) + print("Result: PASS" if passed else "Result: FAIL") + return 0 if passed else 1 + + +if __name__ == "__main__": + raise SystemExit(main())