Add VQM24 transition-state benchmark

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:
 

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

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)

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

docs/benchmarks/vqm24_transition_states.rst

@@ -0,0 +1,93 @@
+VQM24 transition-state benchmark
+================================
+
+The `VQM24 dataset <https://doi.org/10.5281/zenodo.15442257>`_
+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 <https://doi.org/10.1038/s41597-025-05428-4>`_
+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.

docs/index.rst

@@ -110,6 +110,7 @@ Features
    benchmarks/entropy_accuracy
    benchmarks/anthraquinone_workflow
    benchmarks/qm9_thermochemistry
+   benchmarks/vqm24_transition_states
 
 Indices
 -------