Add transition-state thermochemistry and Eyring rate constants

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

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

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

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,23 +549,47 @@ 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
 
         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

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