Add external thermochemistry benchmarks

README.md

@@ -138,6 +138,26 @@ 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
+```
+
+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)
+and [QM9 thermochemistry benchmark](docs/benchmarks/qm9_thermochemistry.rst)
+for their scope and interpretation.
+
 Run linting:
 
 ```bash

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

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

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
+<https://doi.org/10.1039/D2CP01717B>`_. The raw calculations and processed
+tables are published in the immutable `Zenodo record 20796838
+<https://doi.org/10.5281/zenodo.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.