From 3fdf18100d35507aa4b6bfee5c3c45e466455005 Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Wed, 8 Jul 2026 08:16:06 +0200 Subject: [PATCH] Add transition-state thermochemistry and Eyring rate constants Thermo.run() previously raised on any imaginary vibrational frequency, so a transition state (a first-order saddle point) could not be evaluated at all. Add Thermo(..., transition_state=True): requires exactly one imaginary mode among the kept frequencies, excludes it from the vibrational partition function/entropy/energy/heat-capacity sums (the reaction coordinate contributes no thermal vibrational term), and exposes its wavenumber via Thermo.imaginary_mode_wavenumber(). Non-TS behavior (the default) is unchanged; temperature_scan propagates the flag. Thread transition_state=False through run_thermo, orca_thermo, cclib_thermo and pyscf_thermo -- the "import an externally-computed structure" engines, since a DFTB+/xtb geometry optimization is a minimizer and cannot land on a saddle point. Add ThermoScreening.thermo.kinetics: - eyring_rate_constant(reactants, ts, temperature, kappa=1.0): TST rate constant from the activation free energy, using the same Thermo/(coefficient, Thermo) stoichiometry convention as reaction_free_energy. - wigner_tunneling_correction(imaginary_wavenumber, temperature): a simple first-order tunneling estimate from the TS's imaginary frequency. Closes #104 --- ThermoScreening/thermo/__init__.py | 1 + ThermoScreening/thermo/api.py | 37 ++++++++- ThermoScreening/thermo/kinetics.py | 92 ++++++++++++++++++++ ThermoScreening/thermo/thermo.py | 74 ++++++++++++++--- docs/api.rst | 8 +- docs/usage.rst | 33 ++++++++ tests/calculator/test_orca.py | 7 ++ tests/calculator/test_pyscf.py | 7 ++ tests/calculator/test_qm.py | 7 ++ tests/thermo/test_api.py | 26 ++++++ tests/thermo/test_kinetics.py | 129 +++++++++++++++++++++++++++++ tests/thermo/test_thermo.py | 94 +++++++++++++++++++++ 12 files changed, 501 insertions(+), 14 deletions(-) create mode 100644 ThermoScreening/thermo/kinetics.py create mode 100644 tests/thermo/test_kinetics.py diff --git a/ThermoScreening/thermo/__init__.py b/ThermoScreening/thermo/__init__.py index 8c0af7d..e6a6f89 100644 --- a/ThermoScreening/thermo/__init__.py +++ b/ThermoScreening/thermo/__init__.py @@ -10,3 +10,4 @@ from .conformers import generate as generate_conformers, write_conformers from .reactions import reaction_free_energy, reduction_potential from .ensemble import boltzmann_weights, ensemble_free_energy, lowest_gibbs +from .kinetics import eyring_rate_constant, wigner_tunneling_correction diff --git a/ThermoScreening/thermo/api.py b/ThermoScreening/thermo/api.py index 79518f9..b196caa 100644 --- a/ThermoScreening/thermo/api.py +++ b/ThermoScreening/thermo/api.py @@ -351,6 +351,7 @@ def run_thermo( atoms=None, spin=None, quasi_rrho=False, + transition_state=False, ): """ Run the thermo calculation. Returns thermo calculation object. @@ -377,6 +378,11 @@ def run_thermo( quasi_rrho : bool If True, use Grimme's quasi-RRHO treatment for the vibrational entropy instead of the pure harmonic oscillator. Default False. + transition_state : bool + If True, treat the geometry as a first-order saddle point: exactly one + imaginary vibrational frequency is required and excluded from the + vibrational thermochemistry (see ``Thermo.imaginary_mode_wavenumber``), + instead of raising. Default False (a minimum is expected). The coordinate file should be in xyz format. @@ -385,11 +391,13 @@ def run_thermo( ------- Thermo The thermo calculation object. - + Raises ------ TSValueError - If the number of vibrational frequencies does not match with the degree of freedom. + If the number of vibrational frequencies does not match with the degree + of freedom, or if the imaginary-mode count doesn't match + ``transition_state`` (see ``Thermo``). """ # Build the atom list either from an optimized ASE Atoms object or by # reading the coordinate file. @@ -424,6 +432,7 @@ def run_thermo( pressure=pressure, engine=engine, quasi_rrho=quasi_rrho, + transition_state=transition_state, ) thermo_setup.run() @@ -439,6 +448,7 @@ def orca_thermo( charge=0.0, spin=None, quasi_rrho=False, + transition_state=False, ): """ Run the thermochemistry from an ORCA ``.hess`` file (DFT-quality data). @@ -465,6 +475,12 @@ def orca_thermo( Spin quantum number S. Defaults to the minimum-spin electron-count guess. quasi_rrho : bool If True, use Grimme's quasi-RRHO vibrational entropy. Default False. + transition_state : bool + If True, treat the file as a first-order saddle point (e.g. an ORCA + OptTS + freq run): its one imaginary frequency is required and excluded + from the vibrational thermochemistry instead of raising. See + ``Thermo.imaginary_mode_wavenumber``. Default False (a minimum is + expected). Returns ------- @@ -497,6 +513,7 @@ def orca_thermo( spin=spin, engine="dftb+", quasi_rrho=quasi_rrho, + transition_state=transition_state, ) @@ -508,6 +525,7 @@ def cclib_thermo( charge=0.0, spin=None, quasi_rrho=False, + transition_state=False, ): """ Run the thermochemistry from a QM output file via cclib. @@ -535,6 +553,12 @@ def cclib_thermo( Spin quantum number S. Defaults to the minimum-spin electron-count guess. quasi_rrho : bool If True, use Grimme's quasi-RRHO vibrational entropy. Default False. + transition_state : bool + If True, treat the file as a first-order saddle point (e.g. a Gaussian + TS optimization + freq run): its one imaginary frequency is required and + excluded from the vibrational thermochemistry instead of raising. See + ``Thermo.imaginary_mode_wavenumber``. Default False (a minimum is + expected). Returns ------- @@ -565,6 +589,7 @@ def cclib_thermo( spin=spin, engine="dftb+", quasi_rrho=quasi_rrho, + transition_state=transition_state, ) @@ -596,6 +621,7 @@ def pyscf_thermo( charge=0.0, spin=None, quasi_rrho=False, + transition_state=False, ): """ Run the thermochemistry from an in-memory PySCF calculation. @@ -627,6 +653,12 @@ def pyscf_thermo( Spin quantum number S. Defaults to the minimum-spin electron-count guess. quasi_rrho : bool If True, use Grimme's quasi-RRHO vibrational entropy. Default False. + transition_state : bool + If True, treat the geometry as a first-order saddle point: its one + imaginary frequency is required and excluded from the vibrational + thermochemistry instead of raising. See + ``Thermo.imaginary_mode_wavenumber``. Default False (a minimum is + expected). Returns ------- @@ -668,6 +700,7 @@ def pyscf_thermo( spin=spin, engine="dftb+", quasi_rrho=quasi_rrho, + transition_state=transition_state, ) diff --git a/ThermoScreening/thermo/kinetics.py b/ThermoScreening/thermo/kinetics.py new file mode 100644 index 0000000..172bd16 --- /dev/null +++ b/ThermoScreening/thermo/kinetics.py @@ -0,0 +1,92 @@ +"""Transition-state-theory rate constants from computed ``Thermo`` objects. + +Pure post-processing helpers that combine the absolute Gibbs free energies +(``Thermo.total_EeGtot()``) of a reactant and a transition state +(``Thermo(..., transition_state=True)``, see :mod:`ThermoScreening.thermo.thermo`) +into an Eyring rate constant, plus a simple tunneling correction. +""" + +import math + +from ..utils.physicalConstants import PhysicalConstants +from .reactions import _total_gibbs + + +def eyring_rate_constant(reactants, ts, temperature=298.15, kappa=1.0): + """ + Eyring transition-state-theory rate constant. + + ``k = kappa * (kB T / h) * exp(-dG-double-dagger / (R T))``, with the + activation free energy ``dG-double-dagger = G(ts) - sum(G(reactants))``. + + Parameters + ---------- + reactants : iterable + The reactant species, in the same convention as + :func:`ThermoScreening.thermo.reactions.reaction_free_energy`: each + entry is a ``Thermo`` (stoichiometric coefficient 1) or a + ``(coefficient, Thermo)`` tuple. Pass a single-entry list for a + unimolecular reaction. + ts : Thermo + The transition state, computed with ``transition_state=True`` (see + :class:`ThermoScreening.thermo.thermo.Thermo`). + temperature : float + Temperature in K, matching the temperature ``reactants`` and ``ts`` + were computed at. Default 298.15. + kappa : float + Transmission coefficient (a tunneling/recrossing correction). Default + 1.0 (no correction); see :func:`wigner_tunneling_correction` for a + simple estimate from the transition state's imaginary frequency. + + Returns + ------- + float + The TST rate constant. For a single reactant this is a first-order + rate constant in s^-1. For multiple reactants (a bimolecular or higher + reaction) this is the pseudo-first-order rate assuming each reactant is + at its computed standard state; converting to a true bimolecular+ rate + constant (e.g. L mol^-1 s^-1) requires an additional standard-state + correction that is not applied here. + """ + delta_g_ddagger = ts.total_EeGtot() - sum( + _total_gibbs(reactant) for reactant in reactants + ) # Hartree per particle + delta_g_j_per_mol = delta_g_ddagger * PhysicalConstants["H"] * PhysicalConstants["N_A"] + return ( + kappa + * (PhysicalConstants["kB"] * temperature / PhysicalConstants["h"]) + * math.exp(-delta_g_j_per_mol / (PhysicalConstants["R"] * temperature)) + ) + + +def wigner_tunneling_correction(imaginary_wavenumber, temperature=298.15): + """ + Wigner's first-order tunneling correction to a TST rate constant. + + ``kappa = 1 + (1/24) * (h c |nu_imag| / (kB T))^2`` (Wigner, Z. Phys. Chem. + B 19, 203 (1932)), a small quantum correction for tunneling along the + reaction coordinate. Valid only for small corrections; for deep tunneling + (large kappa) use a more complete treatment (e.g. Eckart). + + Parameters + ---------- + imaginary_wavenumber : float + The transition state's imaginary-mode wavenumber in cm^-1 (as returned + by ``Thermo.imaginary_mode_wavenumber()``; the sign is ignored). + temperature : float + Temperature in K. Default 298.15. + + Returns + ------- + float + The dimensionless transmission coefficient kappa (>= 1), for use as + ``eyring_rate_constant(..., kappa=...)``. + """ + u = ( + PhysicalConstants["h"] + * PhysicalConstants["c"] + * abs(imaginary_wavenumber) + * 10**2 # cm^-1 -> m^-1 + / (PhysicalConstants["kB"] * temperature) + ) + return 1.0 + u**2 / 24.0 diff --git a/ThermoScreening/thermo/thermo.py b/ThermoScreening/thermo/thermo.py index 710471a..be42686 100644 --- a/ThermoScreening/thermo/thermo.py +++ b/ThermoScreening/thermo/thermo.py @@ -62,6 +62,7 @@ def __init__( system: System, engine: str, quasi_rrho: bool = False, + transition_state: bool = False, ): """ Initializes the Thermo class with the temperature, pressure, system information @@ -81,6 +82,14 @@ def __init__( If True, use Grimme's quasi-RRHO treatment for the vibrational entropy (interpolating low-frequency modes towards a free rotor) instead of the pure rigid-rotor-harmonic-oscillator model. + transition_state : bool + If True, treat the geometry as a first-order saddle point: exactly + one imaginary (negative) vibrational frequency is required and + excluded from the vibrational partition function (the reaction + coordinate contributes no thermal vibrational term), instead of + raising. Its wavenumber is exposed via + :meth:`imaginary_mode_wavenumber`. Default False (a minimum is + expected; any imaginary frequency raises). Raises ------ @@ -112,6 +121,8 @@ def __init__( self._system = system self._engine = engine self._quasi_rrho = quasi_rrho + self._transition_state = transition_state + self._imaginary_mode_wavenumber = None if self._engine not in ("dftb+", "xtb"): raise TSValueError("The engine is not supported.") @@ -151,9 +162,9 @@ def temperature_scan(self, temperatures): The electronic energy, geometry, and vibrational frequencies are temperature-independent, so this reuses the same :class:`System` (and - this object's pressure, engine, and quasi-RRHO setting) and only the - thermal terms are recomputed -- i.e. a temperature scan from a single - Hessian. This object is left unchanged. + this object's pressure, engine, quasi-RRHO and transition-state setting) + and only the thermal terms are recomputed -- i.e. a temperature scan from + a single Hessian. This object is left unchanged. Parameters ---------- @@ -174,6 +185,7 @@ def temperature_scan(self, temperatures): system=self._system, engine=self._engine, quasi_rrho=self._quasi_rrho, + transition_state=self._transition_state, ) thermo.run() scan.append(thermo) @@ -388,7 +400,7 @@ def _compute_vibrational_partition_function(self): """ self._vib_temp_K = ( PhysicalConstants["h"] - * self._system.real_vibrational_frequencies + * self._real_vibrational_frequencies * PhysicalConstants["c"] * 10**2 / (PhysicalConstants["kB"]) @@ -468,7 +480,7 @@ def _quasi_rrho_entropy(self, harmonic): weight = 1.0 / ( 1.0 - + (self._QRRHO_FREQ_CM / self._system.real_vibrational_frequencies) ** 4 + + (self._QRRHO_FREQ_CM / self._real_vibrational_frequencies) ** 4 ) return weight * harmonic + (1.0 - weight) * free_rotor @@ -537,6 +549,11 @@ def _vibrational_contribution(self): """ Computes the vibrational contribution of the system. + For a transition state (``transition_state=True``), the one required + imaginary mode is set aside (stored via :meth:`imaginary_mode_wavenumber`) + and excluded from the vibrational sums below, which then run over the + remaining real modes only. + Returns ------- None @@ -544,16 +561,35 @@ def _vibrational_contribution(self): Raises ------ TSValueError - If a kept vibrational frequency is imaginary (non-positive), which - would otherwise make the harmonic formulas return NaN. - """ - - if np.any(self._system.real_vibrational_frequencies <= 0): + If ``transition_state`` is False and a kept vibrational frequency is + imaginary (non-positive), which would otherwise make the harmonic + formulas return NaN. + If ``transition_state`` is True and the kept frequencies do not have + exactly one imaginary (non-positive) mode, i.e. the geometry is not a + first-order saddle point. + """ + + frequencies = self._system.real_vibrational_frequencies + imaginary = frequencies[frequencies <= 0] + + if self._transition_state: + if len(imaginary) != 1: + raise TSValueError( + "A transition state must have exactly one imaginary " + f"vibrational frequency (a first-order saddle point); found " + f"{len(imaginary)}: {list(imaginary)}." + ) + self._imaginary_mode_wavenumber = float(imaginary[0]) + frequencies = frequencies[frequencies > 0] + elif len(imaginary): raise TSValueError( "Imaginary (non-positive) vibrational frequencies are present; " - "the geometry is not a minimum." + "the geometry is not a minimum. Pass transition_state=True if " + "this is intentionally a first-order saddle point." ) + self._real_vibrational_frequencies = frequencies + self._compute_vibrational_partition_function() self._compute_vibrational_entropy() self._compute_vibrational_energy() @@ -1018,3 +1054,19 @@ def electronic_energy(self) -> float: """ return _real_scalar(self._system.electronic_energy) + + def imaginary_mode_wavenumber(self): + """ + The transition state's imaginary-mode wavenumber (cm^-1, negative). + + Set by :meth:`run` when this object was constructed with + ``transition_state=True``; useful for a tunneling correction (see + :func:`ThermoScreening.thermo.kinetics.wigner_tunneling_correction`). + + Returns + ------- + float or None + The imaginary mode's wavenumber, or ``None`` for a non-transition- + state ``Thermo`` (the default) or before ``run()`` has been called. + """ + return self._imaginary_mode_wavenumber diff --git a/docs/api.rst b/docs/api.rst index e266f58..8140c1b 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -33,6 +33,12 @@ Reactions and redox .. autofunction:: ThermoScreening.thermo.reactions.reaction_free_energy .. autofunction:: ThermoScreening.thermo.reactions.reduction_potential +Transition states and kinetics +------------------------------- + +.. autofunction:: ThermoScreening.thermo.kinetics.eyring_rate_constant +.. autofunction:: ThermoScreening.thermo.kinetics.wigner_tunneling_correction + Conformer ensembles ------------------- @@ -48,7 +54,7 @@ Thermochemistry core .. autoclass:: ThermoScreening.thermo.thermo.Thermo :members: total_energy, total_enthalpy, total_gibbs_free_energy, total_entropy, total_heat_capacity, total_EeGtot, electronic_energy, - temperature_scan + temperature_scan, imaginary_mode_wavenumber Coordinate and frequency readers -------------------------------- diff --git a/docs/usage.rst b/docs/usage.rst index 74b24b5..ea7ec19 100644 --- a/docs/usage.rst +++ b/docs/usage.rst @@ -200,6 +200,39 @@ conditions), combine them into reaction free energies and reduction potentials: across similar species, with a higher-accuracy method, or with a ``reference_potential`` calibrated against experiment. +Transition states and rate constants +------------------------------------- + +By default ``Thermo`` requires a minimum (no imaginary frequencies). Pass +``transition_state=True`` to the ``*_thermo`` functions to evaluate a +first-order saddle point instead: its one required imaginary (reaction +coordinate) mode is excluded from the vibrational thermochemistry rather than +raising, and exposed via ``Thermo.imaginary_mode_wavenumber()``. This works +with any engine that imports an externally-computed structure -- ``orca_thermo``, +``cclib_thermo`` and ``pyscf_thermo`` -- since a DFTB+/xtb geometry +*optimization* is a minimizer and cannot itself locate a saddle point: + +.. code-block:: python + + from ThermoScreening.thermo.api import orca_thermo + from ThermoScreening.thermo import eyring_rate_constant, wigner_tunneling_correction + + reactant = orca_thermo("reactant.hess") + ts = orca_thermo("ts.hess", transition_state=True) + + nu_imag = ts.imaginary_mode_wavenumber() # cm^-1, negative + kappa = wigner_tunneling_correction(nu_imag, temperature=298.15) + k = eyring_rate_constant([reactant], ts, temperature=298.15, kappa=kappa) + +``eyring_rate_constant`` uses the same reactant convention as +``reaction_free_energy`` (a list of ``Thermo`` or ``(coefficient, Thermo)`` +entries), so a bimolecular reaction is ``eyring_rate_constant([a, b], ts)``. +For a single reactant the result is a first-order rate constant in s\ :sup:`-1`; +for multiple reactants it is the pseudo-first-order TST rate (no standard-state +correction is applied). ``wigner_tunneling_correction`` is a small, first-order +tunneling estimate -- valid only when it stays close to 1; for deep tunneling +use a more complete treatment. + End-to-end example ------------------ diff --git a/tests/calculator/test_orca.py b/tests/calculator/test_orca.py index 0b6511c..1f51578 100644 --- a/tests/calculator/test_orca.py +++ b/tests/calculator/test_orca.py @@ -110,6 +110,13 @@ def test_orca_thermo_energy_override(tmp_path): assert thermo.electronic_energy() == pytest.approx(-77.0) +def test_orca_thermo_transition_state(tmp_path): + # mode 6 (the softest "vibrational" mode) is imaginary -> a TS + ts_hess = _HESS.replace(" 6 1600.000000", " 6 -300.000000") + thermo = orca_thermo(_write_hess(tmp_path, ts_hess, "ts.hess"), transition_state=True) + assert thermo.imaginary_mode_wavenumber() == pytest.approx(-300.0) + + _CO2_HESS = """\ $act_energy -188.500000 diff --git a/tests/calculator/test_pyscf.py b/tests/calculator/test_pyscf.py index 674c4d2..855b090 100644 --- a/tests/calculator/test_pyscf.py +++ b/tests/calculator/test_pyscf.py @@ -40,6 +40,13 @@ def test_pyscf_thermo_energy_override(): assert thermo.electronic_energy() == pytest.approx(-77.0) +def test_pyscf_thermo_transition_state(): + thermo = pyscf_thermo( + _FakeMeanField(), frequencies=[-300.0, 1600.0, 3800.0], transition_state=True, + ) + assert thermo.imaginary_mode_wavenumber() == pytest.approx(-300.0) + + def test_pyscf_thermo_requires_frequencies_or_hessian(): with pytest.raises(TSValueError, match="Provide frequencies"): pyscf_thermo(_FakeMeanField()) diff --git a/tests/calculator/test_qm.py b/tests/calculator/test_qm.py index 83727e8..27e054f 100644 --- a/tests/calculator/test_qm.py +++ b/tests/calculator/test_qm.py @@ -87,6 +87,13 @@ def test_cclib_thermo_uses_file_energy(monkeypatch, tmp_path): assert math.isfinite(thermo.total_EeGtot()) +def test_cclib_thermo_transition_state(monkeypatch, tmp_path): + ts_data = dict(_WATER, vibfreqs=np.array([-300.0, 1600.0, 3800.0])) + _fake_ccread(monkeypatch, **ts_data) + thermo = cclib_thermo(str(tmp_path / "ts.log"), transition_state=True) + assert thermo.imaginary_mode_wavenumber() == pytest.approx(-300.0) + + def test_cclib_thermo_energy_override(monkeypatch, tmp_path): _fake_ccread(monkeypatch, **_WATER) thermo = cclib_thermo(str(tmp_path / "water.log"), energy=-77.0) diff --git a/tests/thermo/test_api.py b/tests/thermo/test_api.py index e9203aa..5a2ad2c 100644 --- a/tests/thermo/test_api.py +++ b/tests/thermo/test_api.py @@ -265,6 +265,32 @@ def test_read_gen_uses_pqanalysis(self, read_gen_file_mock): assert cell is None assert pbc is False + def test_run_thermo_transition_state(self): + # bent (non-linear) triatomic, dof=3: 6 near-zero trans/rot modes then + # one imaginary (reaction-coordinate) mode and two real vibrations + water = Atoms( + "H2O", + positions=[[0.0, 0.76, -0.48], [0.0, -0.76, -0.48], [0.0, 0.0, 0.12]], + ) + frequencies = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -300.0, 1600.0, 3800.0]) + + thermo = run_thermo( + frequencies, atoms=water, engine="dftb+", energy=-76.0, + transition_state=True, + ) + + assert thermo.imaginary_mode_wavenumber() == pytest.approx(-300.0) + + def test_run_thermo_non_transition_state_rejects_imaginary_frequency(self): + water = Atoms( + "H2O", + positions=[[0.0, 0.76, -0.48], [0.0, -0.76, -0.48], [0.0, 0.0, 0.12]], + ) + frequencies = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -300.0, 1600.0, 3800.0]) + + with pytest.raises(TSValueError, match="geometry is not a minimum"): + run_thermo(frequencies, atoms=water, engine="dftb+", energy=-76.0) + def test_run_thermo_rejects_wrong_frequency_count(self): coord_file = Path(__file__).resolve().parents[1] / "data/thermo/geo_opt.xyz" diff --git a/tests/thermo/test_kinetics.py b/tests/thermo/test_kinetics.py new file mode 100644 index 0000000..67fb56b --- /dev/null +++ b/tests/thermo/test_kinetics.py @@ -0,0 +1,129 @@ +import math + +import pytest + +from ThermoScreening.thermo.kinetics import ( + eyring_rate_constant, + wigner_tunneling_correction, +) + +_KB_OVER_H = 1.380649e-23 / 6.62607015e-34 # kB/h, s^-1 K^-1 + + +class _FakeThermo: + """A stand-in exposing only what the kinetics helpers use.""" + + def __init__(self, eegtot): + self._eegtot = eegtot + + def total_EeGtot(self): + return self._eegtot + + +def test_eyring_rate_constant_zero_barrier_is_kT_over_h(): + reactant = _FakeThermo(-10.0) + ts = _FakeThermo(-10.0) # dG-ddagger = 0 + temperature = 298.15 + + k = eyring_rate_constant([reactant], ts, temperature=temperature) + + assert k == pytest.approx(_KB_OVER_H * temperature) + + +def test_eyring_rate_constant_known_barrier(): + # solve k = kB T / h * exp(-dG/RT) = 1 analytically for dG, independently + # of eyring_rate_constant, and check it reproduces k = 1 s^-1 + R, N_A, H = 8.314462618, 6.02214076e23, 4.359744722207101e-18 + temperature = 298.15 + dG_j_per_mol = R * temperature * math.log(_KB_OVER_H * temperature) + dG_hartree = dG_j_per_mol / (H * N_A) + + reactant = _FakeThermo(-10.0) + ts = _FakeThermo(-10.0 + dG_hartree) + + k = eyring_rate_constant([reactant], ts, temperature=temperature) + + assert k == pytest.approx(1.0) + + +def test_eyring_rate_constant_stoichiometric_reactants(): + # bimolecular convention: reactants is a list of Thermo/(coeff, Thermo); + # only the total G(reactants) matters, split however + ts = _FakeThermo(-9.0) + single = [_FakeThermo(-10.0)] + split = [_FakeThermo(-6.0), _FakeThermo(-4.0)] + + k_single = eyring_rate_constant(single, ts, temperature=300.0) + k_split = eyring_rate_constant(split, ts, temperature=300.0) + + assert k_single == pytest.approx(k_split) + + +def test_eyring_rate_constant_stoichiometry_tuple(): + ts = _FakeThermo(-9.0) + reactant = _FakeThermo(-5.0) + # 2 A -> TS, same total G as two separate -5.0 reactants + k_tuple = eyring_rate_constant([(2.0, reactant)], ts, temperature=300.0) + k_split = eyring_rate_constant( + [_FakeThermo(-5.0), _FakeThermo(-5.0)], ts, temperature=300.0 + ) + + assert k_tuple == pytest.approx(k_split) + + +def test_eyring_rate_constant_kappa_scales_linearly(): + reactant = _FakeThermo(-10.0) + ts = _FakeThermo(-9.9) + + k1 = eyring_rate_constant([reactant], ts, temperature=298.15, kappa=1.0) + k2 = eyring_rate_constant([reactant], ts, temperature=298.15, kappa=2.5) + + assert k2 == pytest.approx(2.5 * k1) + + +def test_eyring_rate_constant_higher_barrier_is_slower(): + reactant = _FakeThermo(-10.0) + low_ts = _FakeThermo(-9.98) + high_ts = _FakeThermo(-9.90) + + k_low = eyring_rate_constant([reactant], low_ts, temperature=298.15) + k_high = eyring_rate_constant([reactant], high_ts, temperature=298.15) + + assert k_low > k_high > 0 + + +def test_wigner_tunneling_correction_small_mode_near_one(): + kappa = wigner_tunneling_correction(-50.0, temperature=298.15) + assert kappa == pytest.approx(1.0, abs=0.01) + assert kappa > 1.0 + + +def test_wigner_tunneling_correction_ignores_sign(): + assert wigner_tunneling_correction(-500.0, 298.15) == pytest.approx( + wigner_tunneling_correction(500.0, 298.15) + ) + + +def test_wigner_tunneling_correction_increases_with_wavenumber(): + small = wigner_tunneling_correction(-100.0, 298.15) + large = wigner_tunneling_correction(-1000.0, 298.15) + assert large > small > 1.0 + + +def test_wigner_tunneling_correction_matches_closed_form(): + # kappa = 1 + (1/24) * (h c |nu| / (kB T))^2, computed independently here + h, c, kB = 6.62607015e-34, 299792458.0, 1.380649e-23 + nu, temperature = 1200.0, 310.0 + u = h * c * nu * 100.0 / (kB * temperature) + expected = 1.0 + u**2 / 24.0 + + assert wigner_tunneling_correction(-nu, temperature) == pytest.approx(expected) + + +def test_public_api_exported(): + from ThermoScreening.thermo import eyring_rate_constant as erc + from ThermoScreening.thermo import wigner_tunneling_correction as wtc + from ThermoScreening.thermo import kinetics + + assert erc is kinetics.eyring_rate_constant + assert wtc is kinetics.wigner_tunneling_correction diff --git a/tests/thermo/test_thermo.py b/tests/thermo/test_thermo.py index 9d1666f..58ec640 100644 --- a/tests/thermo/test_thermo.py +++ b/tests/thermo/test_thermo.py @@ -1,5 +1,6 @@ import pytest import os +import math import unittest import numpy as np @@ -111,6 +112,99 @@ def _valid_system(): ) +def _bent_system(vibrational_frequencies): + # a non-linear (right-angle) triatomic, so dof = 3*3 - 6 = 3 -- unlike + # _valid_system's collinear atoms, this lets the kept (last-dof) frequencies + # include an imaginary mode without it being a trans/rot mode + atoms = [ + Atom(symbol="H", position=np.array([0.0, 0.0, 0.0])), + Atom(symbol="H", position=np.array([1.0, 0.0, 0.0])), + Atom(symbol="H", position=np.array([0.0, 1.0, 0.0])), + ] + return System( + atoms, + periodicity=False, + cell=None, + charge=0, + electronic_energy=-1.0, + vibrational_frequencies=vibrational_frequencies, + ) + + +def _valid_ts_system(): + # 6 near-zero trans/rot modes, then one imaginary (reaction-coordinate) mode + # and two real vibrations -- a first-order saddle point + return _bent_system(np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -50.0, 1500.0, 3000.0])) + + +def test_thermo_transition_state_excludes_imaginary_mode(): + thermo = Thermo( + temperature=298.15, pressure=101325, system=_valid_ts_system(), + engine="dftb+", transition_state=True, + ) + thermo.run() + + assert thermo.imaginary_mode_wavenumber() == pytest.approx(-50.0) + # the vibrational sums ran over only the two real modes + assert list(thermo._real_vibrational_frequencies) == [1500.0, 3000.0] + assert math.isfinite(thermo.total_EeGtot()) + + +def test_thermo_transition_state_temperature_scan_propagates_flag(): + thermo = Thermo( + temperature=298.15, pressure=101325, system=_valid_ts_system(), + engine="dftb+", transition_state=True, + ) + thermo.run() + + scan = thermo.temperature_scan([280.0, 320.0]) + assert all(t.imaginary_mode_wavenumber() == pytest.approx(-50.0) for t in scan) + + +def test_thermo_non_transition_state_rejects_imaginary_mode(): + with pytest.raises(TSValueError, match="geometry is not a minimum"): + Thermo( + temperature=298.15, pressure=101325, system=_valid_ts_system(), + engine="dftb+", # transition_state defaults to False + ).run() + + +def test_thermo_transition_state_rejects_wrong_imaginary_count(): + # two imaginary modes: not a valid first-order saddle point + two_imaginary = _bent_system( + np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -50.0, -60.0, 3000.0]) + ) + with pytest.raises(TSValueError, match="exactly one imaginary"): + Thermo( + temperature=298.15, pressure=101325, system=two_imaginary, + engine="dftb+", transition_state=True, + ).run() + + # zero imaginary modes: also not a valid saddle point + zero_imaginary = _bent_system( + np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 100.0, 1500.0, 3000.0]) + ) + with pytest.raises(TSValueError, match="exactly one imaginary"): + Thermo( + temperature=298.15, pressure=101325, system=zero_imaginary, + engine="dftb+", transition_state=True, + ).run() + + +def test_thermo_imaginary_mode_wavenumber_none_before_run_or_for_minimum(): + ts_thermo = Thermo( + temperature=298.15, pressure=101325, system=_valid_ts_system(), + engine="dftb+", transition_state=True, + ) + assert ts_thermo.imaginary_mode_wavenumber() is None # before run() + + minimum = Thermo( + temperature=298.15, pressure=101325, system=_valid_system(), engine="dftb+" + ) + minimum.run() + assert minimum.imaginary_mode_wavenumber() is None # not a transition state + + def test_thermo_rejects_unsupported_engine(): with pytest.raises(TSValueError, match="engine is not supported"): Thermo(temperature=298.15, pressure=101325, system=_valid_system(), engine="orca")