Add a GFN-xTB engine (tblite) alongside DFTB+

ThermoScreening/calculator/xtb.py

@@ -0,0 +1,101 @@
+"""xTB (GFN) calculator via tblite: geometry optimisation and vibrations.
+
+Runs entirely in-process through the tblite ASE calculator (no external binary
+and no Slater-Koster files -- GFN-xTB parameters are built in for H-Rn). The
+geometry is optimised with an ASE optimiser and the frequencies come from ASE's
+finite-difference vibrational analysis, so the result plugs into ``run_thermo``
+exactly like the DFTB+ path.
+"""
+
+import numpy as np
+
+from ..utils.physicalConstants import PhysicalConstants
+
+
+def _eV_to_hartree(energy_ev):
+    """
+    Convert an energy in eV (ASE's unit) to Hartree (the thermo unit).
+    """
+
+    return energy_ev * PhysicalConstants["eV"] / PhysicalConstants["H"]
+
+
+def _real_frequencies_cm(frequencies):
+    """
+    Convert ASE's complex vibrational frequencies to a sorted real array.
+
+    ASE returns frequencies in cm^-1 as a complex array where imaginary modes
+    carry the value in the imaginary part. They are mapped to negative real
+    frequencies (the convention the ``System`` uses) and sorted ascending, so
+    the translational/rotational and any imaginary modes sit at the bottom and
+    ``System`` keeps the highest ``dof`` real vibrations.
+
+    Parameters
+    ----------
+    frequencies : np.ndarray
+        Complex vibrational frequencies in cm^-1 (from ``Vibrations``).
+
+    Returns
+    -------
+    np.ndarray
+        Sorted real frequencies in cm^-1 (imaginary modes negative).
+    """
+    frequencies = np.asarray(frequencies)
+    real = np.where(
+        np.abs(frequencies.imag) > 1e-6,
+        -np.abs(frequencies.imag),
+        frequencies.real,
+    )
+    return np.sort(real)
+
+
+def optimise_and_frequencies(atoms, calc, fmax=0.01):
+    """
+    Optimise ``atoms`` with ``calc`` and return geometry, energy, frequencies.
+
+    Parameters
+    ----------
+    atoms : ase.Atoms
+        Initial geometry (with any ``info['charge']`` / ``info['spin']`` already
+        set for the calculator).
+    calc : ase.calculators.calculator.Calculator
+        The ASE calculator to attach (e.g. a tblite ``TBLite`` instance).
+    fmax : float
+        Force convergence threshold for the optimisation (eV/A).
+
+    Returns
+    -------
+    tuple(ase.Atoms, float, np.ndarray)
+        The optimised atoms, the energy in Hartree, and the sorted real
+        vibrational frequencies in cm^-1.
+    """
+    from ase.optimize import BFGS
+    from ase.vibrations import Vibrations
+
+    atoms = atoms.copy()
+    atoms.calc = calc
+
+    BFGS(atoms, logfile=None).run(fmax=fmax)
+    energy_hartree = _eV_to_hartree(atoms.get_potential_energy())
+
+    vibrations = Vibrations(atoms, name="xtb_vib")
+    vibrations.run()
+    frequencies = _real_frequencies_cm(vibrations.get_frequencies())
+    vibrations.clean()
+
+    return atoms, energy_hartree, frequencies
+
+
+def xtb_calculator(method="GFN2-xTB"):
+    """
+    Build a tblite GFN-xTB ASE calculator (charge and spin are read from
+    ``atoms.info``).
+
+    Parameters
+    ----------
+    method : str
+        ``"GFN2-xTB"`` or ``"GFN1-xTB"``.
+    """
+    from tblite.ase import TBLite  # pragma: no cover - optional heavy dependency
+
+    return TBLite(method=method, verbosity=0)  # pragma: no cover

ThermoScreening/cli/thermo.py

@@ -111,6 +111,15 @@ def _command_parser():
         help="Use Grimme's quasi-RRHO vibrational entropy (better for low-frequency "
         "modes) instead of the pure harmonic oscillator.",
     )
+    screen_parser.add_argument(
+        "--engine", default="dftb+", choices=["dftb+", "xtb"],
+        help="Calculation engine (default 'dftb+'). 'xtb' uses GFN-xTB via tblite "
+        "(no Slater-Koster files); parameter-set/solvent apply only to dftb+.",
+    )
+    screen_parser.add_argument(
+        "--method", default="GFN2-xTB", choices=["GFN2-xTB", "GFN1-xTB"],
+        help="GFN-xTB parametrisation when --engine xtb (default 'GFN2-xTB').",
+    )
 
     return parser
 
@@ -193,6 +202,8 @@ def run_screen(parser_args):
         parameter_set=parser_args.parameter_set,
         solvent=parser_args.solvent,
         quasi_rrho=parser_args.quasi_rrho,
+        engine=parser_args.engine,
+        method=parser_args.method,
     )
 
     failed = sum(1 for record in results if record["status"] != "ok")

ThermoScreening/thermo/api.py

@@ -19,6 +19,7 @@
 from .atoms import Atom
 from ..calculator import Geoopt, Hessian, Modes
 from ..calculator.dftbplus import _spin_kwargs, _solvation_kwargs, SPIN_CONSTANTS_3OB
+from ..calculator.xtb import optimise_and_frequencies, xtb_calculator
 
 
 logger = logging.getLogger(__package_name__).getChild("api")
@@ -612,3 +613,82 @@ def dftbplus_thermo(
         )
 
     return thermo
+
+
+def xtb_thermo(
+        atoms,
+        temperature=298.15,
+        pressure=101325,
+        charge=0.0,
+        spin=None,
+        method="GFN2-xTB",
+        directory=None,
+        quasi_rrho=False,
+        fmax=0.01,
+    ):
+    """
+    Run the thermo pipeline with GFN-xTB (tblite): optimise, compute the
+    Hessian/frequencies, and evaluate the thermochemistry.
+
+    Needs no Slater-Koster files or spin constants -- GFN-xTB parameters are
+    built in for H-Rn and open-shell systems are handled natively via the number
+    of unpaired electrons.
+
+    Parameters
+    ----------
+    atoms : ase.Atoms
+        Atoms object. The initial geometry.
+    temperature : float
+        The temperature in K. Default is 298.15.
+    pressure : float
+        The pressure in Pa. Default is 101325.
+    charge : float
+        The system charge. Default is 0.0.
+    spin : float, optional
+        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.
+    method : str
+        GFN-xTB parametrisation, ``"GFN2-xTB"`` (default) or ``"GFN1-xTB"``.
+    directory : str, optional
+        Working directory to run in (created if needed). Defaults to the current
+        directory.
+    quasi_rrho : bool
+        If True, use Grimme's quasi-RRHO vibrational entropy. Default False.
+    fmax : float
+        Force convergence threshold for the optimisation (eV/A). Default 0.01.
+
+    Returns
+    -------
+    Thermo
+        The thermo calculation object.
+    """
+
+    # Resolve the spin the same way as the DFTB+ path so the calculation and the
+    # analysis use the same multiplicity.
+    if spin is None:
+        electrons = round(float(sum(atoms.get_atomic_numbers())) - charge)
+        spin = 0.0 if electrons % 2 == 0 else 0.5
+
+    prepared = atoms.copy()
+    prepared.info["charge"] = int(round(charge))
+    prepared.info["spin"] = int(round(2.0 * float(spin)))
+
+    with _run_in_directory(directory):
+        optimized_atoms, potential_energy, frequencies = optimise_and_frequencies(
+            prepared, xtb_calculator(method), fmax=fmax
+        )
+
+        thermo = run_thermo(
+            frequencies,
+            atoms=optimized_atoms,
+            temperature=temperature,
+            pressure=pressure,
+            energy=potential_energy,
+            engine="xtb",
+            charge=charge,
+            spin=spin,
+            quasi_rrho=quasi_rrho,
+        )
+
+    return thermo

ThermoScreening/thermo/screening.py

@@ -20,7 +20,7 @@
 from ThermoScreening.exceptions import TSValueError
 from ThermoScreening.utils.custom_logging import setup_logger
 
-from .api import dftbplus_thermo
+from .api import dftbplus_thermo, xtb_thermo
 
 logger = logging.getLogger(__package_name__).getChild("screening")
 logger = setup_logger(logger)
@@ -144,6 +144,8 @@ def screen(
     parameter_set="3ob",
     solvent=None,
     quasi_rrho=False,
+    engine="dftb+",
+    method="GFN2-xTB",
 ):
     """
     Run a thermochemistry screen over a set of molecules.
@@ -180,14 +182,26 @@ def screen(
         If True, use Grimme's quasi-RRHO treatment for the vibrational entropy
         (recommended for flexible molecules with low-frequency modes). Default
         False (pure harmonic oscillator).
+    engine : str
+        Calculation engine: ``"dftb+"`` (default) or ``"xtb"`` (GFN-xTB via
+        tblite). The ``parameter_set``/``solvent`` options apply only to DFTB+.
+    method : str
+        GFN-xTB parametrisation when ``engine="xtb"`` (``"GFN2-xTB"`` default).
 
     Returns
     -------
     list of dict
         One result record per molecule, including failed ones (status="error").
     """
-    default_parameters, spin_constants = resolve_parameter_set(parameter_set)
-    parameters = default_parameters if parameters is None else parameters
+    if engine not in ("dftb+", "xtb"):
+        raise TSValueError(
+            f"Unknown engine {engine!r}; choose 'dftb+' or 'xtb'."
+        )
+
+    if engine == "dftb+":
+        default_parameters, spin_constants = resolve_parameter_set(parameter_set)
+        parameters = default_parameters if parameters is None else parameters
+
     jobs = _load_jobs(source, charge, spin)
     root = Path(directory)
 
@@ -203,18 +217,30 @@ def screen(
         logger.info(f"Screening {job.name} (charge {job.charge})")
         try:
             atoms = ase.io.read(str(job.path))
-            thermo = dftbplus_thermo(
-                atoms,
-                temperature=temperature,
-                pressure=pressure,
-                charge=job.charge,
-                directory=str(root / job.name),
-                spin=job.spin,
-                spin_constants=spin_constants,
-                solvent=solvent,
-                quasi_rrho=quasi_rrho,
-                **parameters,
-            )
+            if engine == "xtb":
+                thermo = xtb_thermo(
+                    atoms,
+                    temperature=temperature,
+                    pressure=pressure,
+                    charge=job.charge,
+                    directory=str(root / job.name),
+                    spin=job.spin,
+                    method=method,
+                    quasi_rrho=quasi_rrho,
+                )
+            else:
+                thermo = dftbplus_thermo(
+                    atoms,
+                    temperature=temperature,
+                    pressure=pressure,
+                    charge=job.charge,
+                    directory=str(root / job.name),
+                    spin=job.spin,
+                    spin_constants=spin_constants,
+                    solvent=solvent,
+                    quasi_rrho=quasi_rrho,
+                    **parameters,
+                )
             record.update(_thermo_summary(thermo))
         except Exception as exc:  # pylint: disable=broad-except
             # isolate failures so one bad molecule does not abort the screen

ThermoScreening/thermo/thermo.py

@@ -113,7 +113,7 @@ def __init__(
         self._engine = engine
         self._quasi_rrho = quasi_rrho
 
-        if self._engine != "dftb+":
+        if self._engine not in ("dftb+", "xtb"):
             raise TSValueError("The engine is not supported.")
         if self._temperature < 0:
             raise TSValueError("The temperature is negative.")

tests/calculator/test_xtb.py

@@ -0,0 +1,60 @@
+import importlib.util
+
+import numpy as np
+import pytest
+
+from ThermoScreening.calculator.xtb import (
+    _eV_to_hartree,
+    _real_frequencies_cm,
+    optimise_and_frequencies,
+)
+
+
+def test_eV_to_hartree_conversion():
+    # 1 Hartree is ~27.2114 eV
+    assert _eV_to_hartree(27.211386) == pytest.approx(1.0, rel=1e-4)
+    assert _eV_to_hartree(0.0) == 0.0
+
+
+def test_real_frequencies_cm_maps_imaginary_negative_and_sorts():
+    # imaginary modes (in the imaginary part) -> negative; then sorted ascending
+    freqs = np.array([1000 + 0j, 0 + 50j, 200 + 0j, 0 + 5j])
+    out = _real_frequencies_cm(freqs)
+
+    assert list(out) == [-50.0, -5.0, 200.0, 1000.0]
+
+
+def test_optimise_and_frequencies_with_emt(monkeypatch, tmp_path):
+    # exercise the real ASE optimiser + finite-difference vibrations with a
+    # dependency-free calculator (EMT), so no tblite is required
+    from ase import Atoms
+    from ase.calculators.emt import EMT
+
+    monkeypatch.chdir(tmp_path)
+    cu2 = Atoms("Cu2", positions=[[0, 0, 0], [0, 0, 2.4]])
+
+    optimized, energy_hartree, frequencies = optimise_and_frequencies(
+        cu2, EMT(), fmax=0.05
+    )
+
+    assert isinstance(energy_hartree, float)
+    assert len(frequencies) == 6  # 3N modes for 2 atoms
+    assert np.all(np.diff(frequencies) >= 0)  # sorted ascending
+    # the single real stretch (top mode) is a positive vibration
+    assert frequencies[-1] > 0
+
+
+tblite_available = importlib.util.find_spec("tblite") is not None
+
+
+@pytest.mark.skipif(not tblite_available, reason="tblite (GFN-xTB) is not installed.")
+def test_xtb_thermo_runs_real_gfn2(tmp_path):
+    # real GFN2-xTB end-to-end: gas-phase water entropy is close to experiment
+    from ase.build import molecule
+    from ThermoScreening.thermo.api import xtb_thermo
+
+    thermo = xtb_thermo(molecule("H2O"), directory=str(tmp_path / "w"))
+    entropy = thermo.total_entropy("cal/(mol*K)")
+
+    # experimental standard molar entropy of gaseous water ~ 45.1 cal/mol/K
+    assert entropy == pytest.approx(45.1, abs=2.0)