Add in-memory thermochemistry from PySCF results

ThermoScreening/thermo/api.py

@@ -568,6 +568,109 @@ def cclib_thermo(
     )
 
 
+def _pyscf_frequencies_from_hessian(mol, hessian):
+    """Vibrational frequencies (cm^-1) from a PySCF Hessian, via pyscf."""
+    try:
+        from pyscf.hessian import thermo as pyscf_thermo_module
+    except ImportError as exc:
+        raise TSValueError(
+            "pyscf is required to compute frequencies from a Hessian; install it "
+            "with 'pip install thermoscreening[pyscf]'."
+        ) from exc
+    # imaginary_freq=False stores imaginary modes as negative real wavenumbers
+    # (the convention the other engines use), so a saddle point isn't silently
+    # flattened into a minimum by dropping a complex frequency.
+    freq_info = pyscf_thermo_module.harmonic_analysis(  # pragma: no cover
+        mol, hessian, imaginary_freq=False
+    )
+    return np.asarray(freq_info["freq_wavenumber"], dtype=float)  # pragma: no cover
+
+
+def pyscf_thermo(
+    mean_field,
+    frequencies=None,
+    hessian=None,
+    energy=None,
+    temperature=298.15,
+    pressure=101325,
+    charge=0.0,
+    spin=None,
+    quasi_rrho=False,
+):
+    """
+    Run the thermochemistry from an in-memory PySCF calculation.
+
+    PySCF is a Python library, so its results live in memory rather than a file.
+    Pass a converged mean-field object; its geometry (Bohr -> Angstrom) and
+    energy (``mean_field.e_tot``) are used, and the vibrational frequencies are
+    taken from ``frequencies`` or computed from ``hessian`` via
+    ``pyscf.hessian.thermo``.
+
+    Parameters
+    ----------
+    mean_field : object
+        A converged PySCF mean-field object (has ``.mol`` and ``.e_tot``).
+    frequencies : array_like, optional
+        Vibrational frequencies in cm^-1. If omitted, ``hessian`` is used.
+    hessian : array_like, optional
+        A PySCF Hessian (from ``mf.Hessian().kernel()``); frequencies are derived
+        from it. Requires the ``pyscf`` extra. Ignored if ``frequencies`` is given.
+    energy : float, optional
+        Electronic energy in Hartree. Defaults to ``mean_field.e_tot``.
+    temperature : float
+        Temperature in K. Default 298.15.
+    pressure : float
+        Pressure in Pa. Default 101325.
+    charge : float
+        System charge. Default 0.0.
+    spin : float, optional
+        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.
+
+    Returns
+    -------
+    Thermo
+        The thermo calculation object.
+
+    Raises
+    ------
+    TSValueError
+        If neither ``frequencies`` nor ``hessian`` is given (or pyscf is missing
+        when a ``hessian`` is used).
+    """
+    from ase import Atoms
+    from ase.units import Bohr
+
+    mol = mean_field.mol
+    symbols = [mol.atom_pure_symbol(index) for index in range(mol.natm)]
+    positions = np.asarray(mol.atom_coords()) * Bohr  # Bohr -> Angstrom
+    atoms = Atoms(symbols=symbols, positions=positions)
+
+    if energy is None:
+        energy = float(mean_field.e_tot)
+
+    if frequencies is None:
+        if hessian is None:
+            raise TSValueError(
+                "Provide frequencies=... or hessian=... to pyscf_thermo."
+            )
+        frequencies = _pyscf_frequencies_from_hessian(mol, hessian)
+    frequencies = np.asarray(frequencies, dtype=float)
+
+    return run_thermo(
+        vibrational_frequencies=frequencies,
+        atoms=atoms,
+        energy=energy,
+        temperature=temperature,
+        pressure=pressure,
+        charge=charge,
+        spin=spin,
+        engine="dftb+",
+        quasi_rrho=quasi_rrho,
+    )
+
+
 def execute(input_file: str) -> Thermo:
     """
     Execute the thermo calculation. Returns thermo calculation object.

docs/api.rst

@@ -9,6 +9,7 @@ High-level pipeline
 .. autofunction:: ThermoScreening.thermo.api.xtb_cli_thermo
 .. autofunction:: ThermoScreening.thermo.api.orca_thermo
 .. autofunction:: ThermoScreening.thermo.api.cclib_thermo
+.. autofunction:: ThermoScreening.thermo.api.pyscf_thermo
 .. autofunction:: ThermoScreening.thermo.api.run_thermo
 .. autofunction:: ThermoScreening.thermo.api.execute
 

docs/usage.rst

@@ -135,6 +135,20 @@ which auto-detects the program via `cclib <https://cclib.github.io/>`_:
 ``read_cclib`` returns the parsed ``(atoms, frequencies, energy)`` if you want
 them directly; ``energy=`` overrides the parsed energy.
 
+**PySCF** results live in memory, so pass the mean-field object directly (with
+the ``pyscf`` extra, ``pip install thermoscreening[pyscf]``):
+
+.. code-block:: python
+
+   from ThermoScreening.thermo.api import pyscf_thermo
+
+   mf = mol.RKS(xc="b3lyp").run()
+   hessian = mf.Hessian().kernel()
+   thermo = pyscf_thermo(mf, hessian=hessian)   # or frequencies=<cm^-1 array>
+
+The geometry and energy (``mf.e_tot``) come from the mean field; frequencies are
+taken from ``frequencies=`` or derived from ``hessian=`` via ``pyscf.hessian``.
+
 Temperature scans
 -----------------
 

pyproject.toml

@@ -49,6 +49,10 @@ Issues = "https://github.com/MolarVerse/ThermoScreening/issues"
 qm = [
     "cclib >= 1.8, < 2.0",
 ]
+# In-memory thermochemistry from PySCF results
+pyscf = [
+    "pyscf",
+]
 test = [
     "pytest",
     "coverage",

tests/calculator/test_pyscf.py

@@ -0,0 +1,88 @@
+import importlib.util
+import math
+
+import numpy as np
+import pytest
+
+_HAS_PYSCF = importlib.util.find_spec("pyscf") is not None
+
+import ThermoScreening.thermo.api as api
+from ThermoScreening.thermo.api import pyscf_thermo, _pyscf_frequencies_from_hessian
+from ThermoScreening.exceptions import TSValueError
+
+
+class _FakeMol:
+    natm = 3
+
+    def atom_coords(self):  # Bohr
+        return np.array([[0.0, 0.0, 0.221722],
+                         [0.0, 1.430901, -0.886569],
+                         [0.0, -1.430901, -0.886569]])
+
+    def atom_pure_symbol(self, index):
+        return ["O", "H", "H"][index]
+
+
+class _FakeMeanField:
+    def __init__(self, e_tot=-76.4):
+        self.mol = _FakeMol()
+        self.e_tot = e_tot
+
+
+def test_pyscf_thermo_from_frequencies():
+    thermo = pyscf_thermo(_FakeMeanField(), frequencies=[1600.0, 3700.0, 3800.0])
+    assert thermo.electronic_energy() == pytest.approx(-76.4)  # from mf.e_tot
+    assert math.isfinite(thermo.total_EeGtot())
+
+
+def test_pyscf_thermo_energy_override():
+    thermo = pyscf_thermo(_FakeMeanField(), frequencies=[1600.0, 3700.0, 3800.0], energy=-77.0)
+    assert thermo.electronic_energy() == pytest.approx(-77.0)
+
+
+def test_pyscf_thermo_requires_frequencies_or_hessian():
+    with pytest.raises(TSValueError, match="Provide frequencies"):
+        pyscf_thermo(_FakeMeanField())
+
+
+def test_pyscf_thermo_hessian_path(monkeypatch):
+    monkeypatch.setattr(
+        api, "_pyscf_frequencies_from_hessian",
+        lambda mol, hessian: np.array([1600.0, 3700.0, 3800.0]),
+    )
+    thermo = pyscf_thermo(_FakeMeanField(), hessian=object())
+    assert math.isfinite(thermo.total_EeGtot())
+
+
+def test_pyscf_frequencies_from_hessian_missing_dependency(monkeypatch):
+    # force the pyscf import inside the helper to fail, robustly -- independent
+    # of whether pyscf is installed or already imported (which would defeat a
+    # sys.modules sentinel) -- by intercepting the import itself
+    import builtins
+
+    real_import = builtins.__import__
+
+    def fake_import(name, *args, **kwargs):
+        if name.startswith("pyscf"):
+            raise ImportError("no pyscf")
+        return real_import(name, *args, **kwargs)
+
+    monkeypatch.setattr(builtins, "__import__", fake_import)
+    with pytest.raises(TSValueError, match="pyscf is required"):
+        _pyscf_frequencies_from_hessian(object(), object())
+
+
+@pytest.mark.skipif(not _HAS_PYSCF, reason="pyscf is not installed")
+def test_pyscf_frequencies_from_hessian_real():
+    # real pyscf: harmonic_analysis on a water HF/STO-3G Hessian returns the
+    # 3N-6 vibrational modes as a real (cm^-1) array (imaginary_freq=False)
+    from pyscf import gto, scf
+
+    mol = gto.M(atom="O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469",
+                basis="sto-3g", verbose=0)
+    mf = scf.RHF(mol).run()
+    hessian = mf.Hessian().kernel()
+
+    frequencies = _pyscf_frequencies_from_hessian(mol, hessian)
+    assert frequencies.dtype == float  # real (negative for any imaginary mode)
+    assert len(frequencies) == 3       # 3N - 6 for a bent triatomic