diff --git a/README.md b/README.md index c8db4a1..3cc5569 100644 --- a/README.md +++ b/README.md @@ -152,11 +152,20 @@ sample from QM9: python scripts/validate_qm9.py ``` -Both commands verify their downloaded archives and keep the external data out +Validate transition-state thermochemistry against a deterministic sample from +the PSI4 VQM24 dataset: + +```bash +python scripts/validate_vqm24.py +``` + +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 [QM9 thermochemistry benchmark](docs/benchmarks/qm9_thermochemistry.rst) -for their scope and interpretation. +and the thermochemistry benchmarks for +[QM9](docs/benchmarks/qm9_thermochemistry.rst) and +[VQM24 transition states](docs/benchmarks/vqm24_transition_states.rst) for +their scope and interpretation. Run linting: diff --git a/ThermoScreening/thermo/api.py b/ThermoScreening/thermo/api.py index 37c9619..3660adf 100644 --- a/ThermoScreening/thermo/api.py +++ b/ThermoScreening/thermo/api.py @@ -14,7 +14,7 @@ from ThermoScreening import __package_name__ from .inputFileReader import InputFileReader -from .system import System, dof +from .system import System, resolved_dof from .thermo import Thermo from .atoms import Atom from ..calculator import Geoopt, Hessian, Modes @@ -416,7 +416,7 @@ def run_thermo( for i, symbol in enumerate(atom_symbol) ] - expected_dof = dof(atom_list) + expected_dof = resolved_dof(atom_list, len(vibrational_frequencies)) if len(vibrational_frequencies) < expected_dof: raise TSValueError( "The number of vibrational frequencies does not match with the degree of freedom." diff --git a/ThermoScreening/thermo/system.py b/ThermoScreening/thermo/system.py index 5f4f74a..7c2a444 100644 --- a/ThermoScreening/thermo/system.py +++ b/ThermoScreening/thermo/system.py @@ -21,6 +21,10 @@ from .cell import Cell +_LINEARITY_TOLERANCE = 1e-9 +_STRICT_LINEARITY_TOLERANCE = 1e-12 + + def _real_position(position: np.ndarray) -> np.ndarray: """ Return a real floating-point atom position for symmetry analysis. @@ -84,6 +88,24 @@ def _pq_atomic_system(atoms: List[Atom]) -> PQAtomicSystem: ) +def _principal_moment_ratio(atoms: List[Atom]) -> float: + masses = np.array([atom.mass for atom in atoms], dtype=float) + positions = np.array([_real_position(atom.position) for atom in atoms]) + positions = positions - np.average(positions, axis=0, weights=masses) + + x, y, z = positions[:, 0], positions[:, 1], positions[:, 2] + inertia_tensor = np.zeros((3, 3)) + inertia_tensor[0, 0] = np.sum(masses * (y**2 + z**2)) + inertia_tensor[1, 1] = np.sum(masses * (x**2 + z**2)) + inertia_tensor[2, 2] = np.sum(masses * (x**2 + y**2)) + inertia_tensor[0, 1] = inertia_tensor[1, 0] = -np.sum(masses * x * y) + inertia_tensor[0, 2] = inertia_tensor[2, 0] = -np.sum(masses * x * z) + inertia_tensor[1, 2] = inertia_tensor[2, 1] = -np.sum(masses * y * z) + + eigenvalues = np.linalg.eigvalsh(inertia_tensor) + return float(eigenvalues[0] / eigenvalues[-1]) + + def linearity(atoms: List[Atom]) -> bool: """ Checks if the system is linear or non-linear. @@ -108,24 +130,9 @@ def linearity(atoms: List[Atom]) -> bool: if len(atoms) == 2: return True - masses = np.array([atom.mass for atom in atoms], dtype=float) - positions = np.array([_real_position(atom.position) for atom in atoms]) - # relocate to the center of mass so linearity is translation-invariant - positions = positions - np.average(positions, axis=0, weights=masses) - - x, y, z = positions[:, 0], positions[:, 1], positions[:, 2] - inertia_tensor = np.zeros((3, 3)) - inertia_tensor[0, 0] = np.sum(masses * (y**2 + z**2)) - inertia_tensor[1, 1] = np.sum(masses * (x**2 + z**2)) - inertia_tensor[2, 2] = np.sum(masses * (x**2 + y**2)) - inertia_tensor[0, 1] = inertia_tensor[1, 0] = -np.sum(masses * x * y) - inertia_tensor[0, 2] = inertia_tensor[2, 0] = -np.sum(masses * x * z) - inertia_tensor[1, 2] = inertia_tensor[2, 1] = -np.sum(masses * y * z) - - 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-9 * eigenvalues[-1]) + return _principal_moment_ratio(atoms) <= _LINEARITY_TOLERANCE def dimensionality(atoms: List[Atom]) -> int: @@ -231,6 +238,27 @@ def dof(atoms: List[Atom]) -> int: raise TSValueError("The number of atoms must be greater than 0.") +def resolved_dof(atoms: List[Atom], frequency_count: int) -> int: + """ + Resolve a false linear classification from an exact nonlinear mode count. + + A projected ``3N-6`` spectrum resolves a geometry inside the near-linear + ambiguity band as nonlinear. Strictly linear geometries remain linear, so + an incomplete spectrum is still rejected. A ``3N-5`` spectrum is not + sufficient to override a nonlinear geometry because it can contain a + residual rotational mode. + """ + number_of_atoms = len(atoms) + nonlinear_dof = 3 * number_of_atoms - 6 + if ( + number_of_atoms >= 3 + and frequency_count == nonlinear_dof + and _principal_moment_ratio(atoms) > _STRICT_LINEARITY_TOLERANCE + ): + return nonlinear_dof + return dof(atoms) + + def default_spin(atoms: List[Atom], charge: float) -> float: """ Guess the spin quantum number S as the minimum-spin ground state from the @@ -669,7 +697,7 @@ def __init__( self._solvent = solvent self._electronic_energy = electronic_energy self._vibrational_frequencies = vibrational_frequencies - self._dof = dof(atoms) + self._dof = resolved_dof(atoms, len(vibrational_frequencies)) self._dim = dim(atoms) self._mass = mass(atoms) self._center_of_mass = center_of_mass(atoms, self._mass) diff --git a/ThermoScreening/thermo/thermo.py b/ThermoScreening/thermo/thermo.py index be42686..22e6c3e 100644 --- a/ThermoScreening/thermo/thermo.py +++ b/ThermoScreening/thermo/thermo.py @@ -9,7 +9,7 @@ from ThermoScreening.utils.physicalConstants import PhysicalConstants from ThermoScreening import __package_name__ -from .system import System, linearity +from .system import System def _real_scalar(value) -> float: @@ -387,7 +387,8 @@ def _rotational_dof(self) -> int: atoms = self._system.atoms if len(atoms) == 1: return 0 - return 2 if linearity(atoms) else 3 + linear_dof = 3 * len(atoms) - 5 + return 2 if self._system.dof == linear_dof else 3 def _compute_vibrational_partition_function(self): diff --git a/docs/benchmarks/vqm24_transition_states.rst b/docs/benchmarks/vqm24_transition_states.rst new file mode 100644 index 0000000..fa0519b --- /dev/null +++ b/docs/benchmarks/vqm24_transition_states.rst @@ -0,0 +1,93 @@ +VQM24 transition-state benchmark +================================ + +The `VQM24 dataset `_ +contains 835,947 converged structures calculated with PSI4 1.7 at the +omegaB97X-D3/cc-pVDZ level. Its chemical space includes C, N, O, F, Si, P, S, +Cl, and Br. The associated +`Scientific Data paper `_ +provides geometries, harmonic frequencies, and thermochemical properties for +minima and saddle structures. + +Run the transition-state benchmark from a source checkout: + +.. code-block:: bash + + python scripts/validate_vqm24.py + +The command downloads the 112 MB ``DFT_saddles.npz`` archive, verifies its MD5 +checksum, and keeps it outside the repository. It scans all 51,072 saddle +records and recomputes thermochemistry for a deterministic sample of 1,000 +first-order transition states. + +Method +------ + +A first-order transition state must contain exactly one imaginary frequency. +The full scan finds 26,445 such records. It excludes 21 records with ``3N-3`` +stored modes because they do not represent a projected vibrational spectrum. + +One near-linear C2NPSi record contains ``3N-5`` stored modes while its +published thermal corrections use three rotational degrees of freedom. The +extra 5.0625 cm\ :sup:`-1` mode makes the source representation internally +inconsistent, so the benchmark reports and excludes the record instead of +changing the thermochemistry model to reproduce it. + +The sample combines the first 100 comparable records, 450 evenly spaced +records, a pinned near-linear C-Si-C regression, and a seeded random remainder. The +electronic baseline is reconstructed as ``U0 - ZPVE`` so the comparison +isolates thermochemistry from inconsistencies between independent source +fields. Calculations use 298.15 K, 101325 Pa, a singlet electronic state, and +rotational symmetry number one. + +Reference result +---------------- + +Energy differences are reported in microhartree. + +.. list-table:: + :header-rows: 1 + + * - Property + - Mean absolute error + - Maximum absolute error + * - ``U0`` + - 0.002560 microhartree + - 0.006189 microhartree + * - ``U298`` + - 0.005651 microhartree + - 0.020289 microhartree + * - ``H`` + - 0.005880 microhartree + - 0.019318 microhartree + * - ``G`` + - 7.164475 microhartree + - 25.527762 microhartree + * - ``S`` + - 0.015089 cal/(mol K) + - 0.053738 cal/(mol K) + * - ``Cv`` + - 0.000256 cal/(mol K) + - 0.000517 cal/(mol K) + +The maximum Gibbs-energy difference is approximately 0.016 kcal/mol. + +Near-linear geometry +-------------------- + +VQM24 array index 2 is a nonlinear C-Si-C transition state with three +projected vibrational modes. Its smallest principal moment is only about +``4.7e-10`` of the largest, below a practical geometry-only linearity cutoff. +An exact ``3N-6`` supplied mode count therefore resolves the geometry as +nonlinear within the near-linear ambiguity band. Strictly linear geometries, +full Cartesian spectra, and ``3N-5`` spectra continue to use the geometry-based +classification. + +Interpretation +-------------- + +This benchmark independently validates transition-state mode handling and +ideal-gas rigid-rotor harmonic-oscillator thermochemistry against PSI4 over a +broader elemental range than QM9. It does not evaluate the predictive accuracy +of the density-functional method or the physical validity of the underlying +saddle structures. diff --git a/docs/index.rst b/docs/index.rst index a2dcc2a..5da0853 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -110,6 +110,7 @@ Features benchmarks/entropy_accuracy benchmarks/anthraquinone_workflow benchmarks/qm9_thermochemistry + benchmarks/vqm24_transition_states Indices ------- diff --git a/scripts/validate_vqm24.py b/scripts/validate_vqm24.py new file mode 100644 index 0000000..c7307d1 --- /dev/null +++ b/scripts/validate_vqm24.py @@ -0,0 +1,290 @@ +#!/usr/bin/env python3 +"""Validate transition-state thermochemistry against the VQM24 dataset.""" + +from __future__ import annotations + +import argparse +import hashlib +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 +from ThermoScreening.utils.physicalConstants import PhysicalConstants + + +RECORD_ID = "15442257" +ARCHIVE_NAME = "DFT_saddles.npz" +DOWNLOAD_URL = ( + f"https://zenodo.org/records/{RECORD_ID}/files/{ARCHIVE_NAME}?download=1" +) +ARCHIVE_MD5 = "82bfaf515f720d45cc5fe03e401b73f4" +RECORD_COUNT = 51_072 +FIRST_ORDER_COUNT = 26_445 +UNEXPECTED_MODE_COUNT = 21 +INCONSISTENT_MODE_ROWS = {38_387} +SAMPLE_SIZE = 1_000 +RANDOM_SEED = 20_260_724 +EDGE_CASE_IDS = {2} +IMAGINARY_TOLERANCE = 1e-8 + +MILLIHARTREE_TO_CAL_PER_MOL = ( + PhysicalConstants["H"] + * PhysicalConstants["N_A"] + / PhysicalConstants["cal"] + / 1000.0 +) + +LIMITS = { + "U0": 0.01, + "U298": 0.03, + "H": 0.03, + "G": 30.0, + "S": 0.06, + "Cv": 0.0006, +} + + +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 _frequencies(raw: np.ndarray) -> np.ndarray: + values = np.asarray(raw, dtype=complex) + return np.where( + np.abs(values.imag) > IMAGINARY_TOLERANCE, + -np.abs(values.imag), + values.real, + ).astype(float) + + +def _atom_list(numbers: np.ndarray, coordinates: np.ndarray) -> list[Atom]: + atoms = Atoms(numbers=numbers, positions=coordinates) + return [ + Atom(symbol=atom.symbol, position=atom.position) + for atom in atoms + ] + + +def _classify_records(data) -> tuple[list[int], list[int], set[int]]: + comparable = [] + unexpected = [] + inconsistent = set() + first_order = 0 + + for index, (numbers, raw) in enumerate(zip(data["atoms"], data["freqs"])): + frequencies = _frequencies(raw) + if np.count_nonzero(frequencies < 0) != 1: + continue + + first_order += 1 + atom_count = len(numbers) + if len(frequencies) not in ( + 3 * atom_count - 6, + 3 * atom_count - 5, + ): + unexpected.append(index) + continue + + dof = resolved_dof( + _atom_list(numbers, data["coordinates"][index]), + len(frequencies), + ) + if len(frequencies) != dof: + inconsistent.add(index) + continue + comparable.append(index) + + if len(data["atoms"]) != RECORD_COUNT: + raise RuntimeError( + f"Expected {RECORD_COUNT} VQM24 saddle records, " + f"found {len(data['atoms'])}." + ) + if first_order != FIRST_ORDER_COUNT: + raise RuntimeError( + f"Expected {FIRST_ORDER_COUNT} first-order saddles, " + f"found {first_order}." + ) + if len(unexpected) != UNEXPECTED_MODE_COUNT: + raise RuntimeError( + f"Expected {UNEXPECTED_MODE_COUNT} records with unexpected " + f"mode counts, found {len(unexpected)}." + ) + if inconsistent != INCONSISTENT_MODE_ROWS: + raise RuntimeError( + "VQM24 geometry/mode inconsistencies changed: " + f"{sorted(inconsistent)}." + ) + return comparable, unexpected, inconsistent + + +def _sample_ids(comparable: list[int]) -> list[int]: + selected = set(comparable[:100]) + selected.update( + np.asarray(comparable)[ + np.linspace(0, len(comparable) - 1, 450, dtype=int) + ].tolist() + ) + selected.update(EDGE_CASE_IDS) + + remaining = np.asarray(sorted(set(comparable) - selected)) + rng = np.random.default_rng(RANDOM_SEED) + selected.update( + rng.choice( + remaining, + size=SAMPLE_SIZE - len(selected), + replace=False, + ).tolist() + ) + if len(selected) != SAMPLE_SIZE: + raise RuntimeError(f"Expected {SAMPLE_SIZE} sampled records.") + return sorted(selected) + + +def _calculate_errors(data, selected: list[int]) -> dict[str, np.ndarray]: + errors = {name: [] for name in LIMITS} + + for index in selected: + atoms = Atoms( + numbers=data["atoms"][index], + positions=data["coordinates"][index], + ) + thermo = run_thermo( + _frequencies(data["freqs"][index]), + atoms=atoms, + energy=( + float(data["U0"][index]) + - float(data["zpves"][index]) + ), + temperature=298.15, + pressure=101325, + spin=0.0, + symmetry_number=1, + transition_state=True, + ) + actual = { + "U0": thermo.total_EeZPE(), + "U298": thermo.total_EeEtot(), + "H": thermo.total_EeHtot(), + "G": thermo.total_EeGtot(), + "S": thermo.total_entropy("cal/(mol*K)"), + "Cv": thermo.total_heat_capacity("cal/(mol*K)"), + } + expected = { + "U0": float(data["U0"][index]), + "U298": float(data["U298"][index]), + "H": float(data["H"][index]), + "G": float(data["G"][index]), + "S": ( + float(data["S"][index]) + * MILLIHARTREE_TO_CAL_PER_MOL + ), + "Cv": float(data["Cv"][index]), + } + for name in LIMITS: + factor = 1e6 if name in {"U0", "U298", "H", "G"} else 1.0 + errors[name].append((actual[name] - expected[name]) * factor) + + 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" / f"vqm24-{RECORD_ID}" + + +def main() -> int: + """Download VQM24, 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()) + required_arrays = ( + "atoms", + "coordinates", + "freqs", + "zpves", + "U0", + "U298", + "H", + "G", + "S", + "Cv", + ) + with np.load(archive, allow_pickle=True) as archive_data: + data = {name: archive_data[name] for name in required_arrays} + + comparable, unexpected, inconsistent = _classify_records(data) + selected = _sample_ids(comparable) + errors = _calculate_errors(data, selected) + + print(f"Zenodo record: {RECORD_ID}") + print(f"VQM24 saddle records scanned: {RECORD_COUNT}") + print(f"First-order saddles: {FIRST_ORDER_COUNT}") + print(f"Unexpected mode-count records excluded: {len(unexpected)}") + print(f"Inconsistent near-linear records excluded: {len(inconsistent)}") + print(f"Transition states sampled: {SAMPLE_SIZE}") + + passed = True + for name, values in errors.items(): + absolute = np.abs(values) + unit = ( + "microhartree" + if name in {"U0", "U298", "H", "G"} + else "cal/(mol*K)" + ) + 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_api.py b/tests/thermo/test_api.py index 48784dd..e1f7b99 100644 --- a/tests/thermo/test_api.py +++ b/tests/thermo/test_api.py @@ -281,6 +281,49 @@ def test_run_thermo_transition_state(self): assert thermo.imaginary_mode_wavenumber() == pytest.approx(-300.0) + def test_run_thermo_resolves_near_linear_transition_state(self): + atoms = Atoms( + numbers=[6, 14, 6], + positions=[ + [-2.7178537e-05, -6.4067e-08, 0.615695131465], + [-2.7178537e-05, -6.4067e-08, -1.080363120435], + [9.0542863e-05, 2.13433e-07, 1.903074839465], + ], + ) + frequencies = np.array([-38.4307, 821.6903, 1973.3601]) + + thermo = run_thermo( + frequencies, + atoms=atoms, + engine="dftb+", + energy=-327.0, + symmetry_number=1, + transition_state=True, + ) + + assert thermo._system.dof == 3 + assert thermo._n_rot == 3 + assert thermo.imaginary_mode_wavenumber() == pytest.approx(-38.4307) + + def test_run_thermo_rejects_incomplete_strictly_linear_spectrum(self): + carbon_dioxide = Atoms( + "OCO", + positions=[ + [0.0, 0.0, -1.16], + [0.0, 0.0, 0.0], + [0.0, 0.0, 1.16], + ], + ) + + with pytest.raises( + TSValueError, + match="number of vibrational frequencies", + ): + run_thermo( + np.array([667.0, 1333.0, 2349.0]), + atoms=carbon_dioxide, + ) + def test_run_thermo_non_transition_state_rejects_imaginary_frequency(self): water = Atoms( "H2O", diff --git a/tests/thermo/test_system.py b/tests/thermo/test_system.py index d3fefef..168d7cf 100644 --- a/tests/thermo/test_system.py +++ b/tests/thermo/test_system.py @@ -4,7 +4,17 @@ import numpy as np from ase.atoms import Atoms import ThermoScreening.thermo.system as system_module -from ThermoScreening.thermo.system import System, dim, dof, linearity, rotational_symmetry_number, default_spin, frequency_dof, check_frequency_length +from ThermoScreening.thermo.system import ( + System, + check_frequency_length, + default_spin, + dim, + dof, + frequency_dof, + linearity, + resolved_dof, + rotational_symmetry_number, +) from ThermoScreening.thermo.atoms import Atom from ThermoScreening.thermo.cell import Cell from ThermoScreening.exceptions import TSValueError @@ -351,6 +361,47 @@ def test_linearity_tolerates_small_coordinate_noise(): ) assert linearity(atoms) is True + assert resolved_dof(atoms, 3) == 4 + + +def test_projected_mode_count_resolves_near_linear_nonlinear_geometry(): + atoms = _atoms( + [ + ("C", (-2.7178537e-05, -6.4067e-08, 0.615695131465)), + ("Si", (-2.7178537e-05, -6.4067e-08, -1.080363120435)), + ("C", (9.0542863e-05, 2.13433e-07, 1.903074839465)), + ] + ) + + assert linearity(atoms) is True + assert resolved_dof(atoms, 3) == 3 + + +def test_extra_projected_mode_does_not_override_nonlinear_geometry(): + atoms = _atoms( + [ + ("C", (1.3501048715, 0.000276513571, -5.1988239e-05)), + ("C", (0.056074661843, 0.000176868465, -3.4554166e-05)), + ("P", (-2.766004783046, -0.000109263102, 2.2827209e-05)), + ("N", (-1.16828267501, 0.000118312157, 1.8129501e-05)), + ("Si", (3.043900054112, -0.000132717544, 2.773555e-06)), + ] + ) + + assert linearity(atoms) is False + assert resolved_dof(atoms, 10) == 9 + + +def test_full_cartesian_mode_count_uses_geometry_classification(): + atoms = _atoms( + [ + ("O", (0.0, 0.0, -1.16)), + ("C", (0.0, 0.0, 0.0)), + ("O", (0.0, 0.0, 1.16)), + ] + ) + + assert resolved_dof(atoms, 9) == dof(atoms) == 4 @pytest.mark.parametrize(