Add a native-xtb engine with implicit solvation for open-shell species

ThermoScreening/calculator/xtb_cli.py

@@ -0,0 +1,153 @@
+"""Native xTB (GFN) engine via the ``xtb`` command-line program.
+
+Unlike the in-process tblite route (:mod:`ThermoScreening.calculator.xtb`), the
+native ``xtb`` binary exposes open-shell (``--uhf``), charge (``--chrg``) and
+implicit solvation (``--alpb``), so it covers charged radicals in solution. A
+single ``xtb --ohess`` call optimises the geometry and computes the Hessian; the
+energy, optimised geometry and frequencies are then read back and fed into
+``run_thermo`` exactly like the other engines.
+"""
+
+import os
+import shutil
+import subprocess
+
+import numpy as np
+import ase.io
+
+# GFN method name -> ``--gfn`` argument
+_GFN_METHODS = {"GFN2-xTB": "2", "GFN1-xTB": "1", "GFN0-xTB": "0"}
+
+
+def resolve_xtb(command=None):
+    """
+    Locate the ``xtb`` executable.
+
+    Order: explicit ``command``, then ``$XTB_COMMAND``, then ``xtb`` on PATH.
+
+    Raises
+    ------
+    FileNotFoundError
+        If no ``xtb`` executable can be found.
+    """
+    candidate = command or os.environ.get("XTB_COMMAND") or shutil.which("xtb")
+    if not candidate:
+        raise FileNotFoundError(
+            "The 'xtb' executable was not found. Install it (e.g. "
+            "`conda install -c conda-forge xtb`) or set the XTB_COMMAND "
+            "environment variable to its path."
+        )
+    return candidate
+
+
+def _parse_optimised_energy(path="xtbopt.xyz"):
+    """
+    Read the optimised total energy (Hartree) from an ``xtbopt.xyz`` comment.
+
+    The second line looks like ``energy: -5.070544 gnorm: ... xtb: 6.7.1``.
+    """
+    with open(path, "r", encoding="utf-8") as handle:
+        handle.readline()  # atom count
+        comment = handle.readline()
+
+    tokens = comment.split()
+    for index, token in enumerate(tokens):
+        if token == "energy:":
+            return float(tokens[index + 1])
+    raise ValueError(f"No energy found in {path!r} comment line: {comment!r}")
+
+
+def _parse_vibspectrum(path="vibspectrum"):
+    """
+    Parse frequencies (cm^-1) from a Turbomole ``vibspectrum`` file.
+
+    Each data line is ``mode [symmetry] wavenumber IR-intensity selection``; the
+    wavenumber is the third-from-last token whether or not a symmetry label is
+    present. Imaginary modes are negative; translational/rotational modes are
+    ~0. The result is sorted ascending so ``System`` keeps the highest ``dof``
+    real vibrations (as with the other engines).
+    """
+    frequencies = []
+    in_block = False
+    with open(path, "r", encoding="utf-8") as handle:
+        for line in handle:
+            stripped = line.strip()
+            if stripped.startswith("$vibrational"):
+                in_block = True
+                continue
+            if stripped.startswith("$end"):
+                break
+            if not in_block or stripped.startswith("#") or not stripped:
+                continue
+            tokens = stripped.split()
+            if len(tokens) >= 4 and tokens[0].isdigit():
+                frequencies.append(float(tokens[-3]))
+    return np.sort(np.array(frequencies))
+
+
+def run_xtb(
+    atoms,
+    charge=0.0,
+    unpaired=0,
+    method="GFN2-xTB",
+    solvent=None,
+    command=None,
+):
+    """
+    Run ``xtb --ohess`` and return the optimised geometry, energy, frequencies.
+
+    Runs in the current working directory (xtb writes several files there), so
+    call it inside a per-job directory.
+
+    Parameters
+    ----------
+    atoms : ase.Atoms
+        The initial geometry.
+    charge : float
+        Total charge (``--chrg``).
+    unpaired : int
+        Number of unpaired electrons, i.e. round(2*S) (``--uhf``).
+    method : str
+        ``"GFN2-xTB"`` (default), ``"GFN1-xTB"`` or ``"GFN0-xTB"``.
+    solvent : str, optional
+        ALPB implicit-solvation solvent (``--alpb``), e.g. ``"water"``.
+    command : str, optional
+        Path to the ``xtb`` executable (defaults to PATH / ``$XTB_COMMAND``).
+
+    Returns
+    -------
+    tuple(ase.Atoms, float, np.ndarray)
+        Optimised atoms, energy in Hartree, sorted real frequencies in cm^-1.
+
+    Raises
+    ------
+    ValueError
+        If ``method`` is unknown.
+    RuntimeError
+        If the xtb run fails.
+    """
+    try:
+        gfn = _GFN_METHODS[method]
+    except KeyError:
+        known = ", ".join(sorted(_GFN_METHODS))
+        raise ValueError(f"Unknown xTB method {method!r}; choose one of: {known}.")
+
+    executable = resolve_xtb(command)
+    ase.io.write("xtb_input.xyz", atoms)
+
+    argv = [
+        executable, "xtb_input.xyz", "--ohess", "--gfn", gfn,
+        "--chrg", str(int(round(charge))), "--uhf", str(int(unpaired)),
+    ]
+    if solvent is not None:
+        argv += ["--alpb", solvent]
+
+    result = subprocess.run(argv, capture_output=True, text=True, check=False)
+    if result.returncode != 0 or not os.path.isfile("xtbopt.xyz"):
+        tail = (result.stdout or "")[-800:] + (result.stderr or "")[-800:]
+        raise RuntimeError(f"xtb failed (exit {result.returncode}):\n{tail}")
+
+    optimized = ase.io.read("xtbopt.xyz")
+    energy_hartree = _parse_optimised_energy("xtbopt.xyz")
+    frequencies = _parse_vibspectrum("vibspectrum")
+    return optimized, energy_hartree, frequencies

ThermoScreening/cli/thermo.py

@@ -112,13 +112,15 @@ def _command_parser():
         "modes) instead of the pure harmonic oscillator.",
     )
     screen_parser.add_argument(
-        "--engine", default="dftb+", choices=["dftb+", "xtb"],
+        "--engine", default="dftb+", choices=["dftb+", "xtb", "xtb-cli"],
         help="Calculation engine (default 'dftb+'). 'xtb' uses GFN-xTB via tblite "
-        "(no Slater-Koster files); parameter-set/solvent apply only to dftb+.",
+        "(gas-phase); 'xtb-cli' uses the native xtb binary, which also does "
+        "implicit solvation of charged radicals. parameter-set applies to dftb+; "
+        "solvent to dftb+ and xtb-cli.",
     )
     screen_parser.add_argument(
         "--method", default="GFN2-xTB", choices=["GFN2-xTB", "GFN1-xTB"],
-        help="GFN-xTB parametrisation when --engine xtb (default 'GFN2-xTB').",
+        help="GFN-xTB parametrisation for the xtb engines (default 'GFN2-xTB').",
     )
 
     return parser

ThermoScreening/thermo/api.py

@@ -20,6 +20,7 @@
 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
+from ..calculator.xtb_cli import run_xtb
 
 
 logger = logging.getLogger(__package_name__).getChild("api")
@@ -692,3 +693,79 @@ def xtb_thermo(
         )
 
     return thermo
+
+
+def xtb_cli_thermo(
+        atoms,
+        temperature=298.15,
+        pressure=101325,
+        charge=0.0,
+        spin=None,
+        method="GFN2-xTB",
+        solvent=None,
+        directory=None,
+        quasi_rrho=False,
+    ):
+    """
+    Run the thermo pipeline with the native ``xtb`` program (GFN-xTB).
+
+    Unlike :func:`xtb_thermo` (in-process tblite, gas-phase), this uses the
+    ``xtb`` binary, which natively supports open-shell (``--uhf``), charge
+    (``--chrg``) and implicit solvation (``--alpb``) -- so charged radicals in
+    solution are handled with method-consistent parameters. A single
+    ``xtb --ohess`` optimises and computes the Hessian.
+
+    Parameters
+    ----------
+    atoms : ase.Atoms
+        The initial geometry.
+    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
+        (even -> 0, odd -> 0.5); ``round(2*S)`` unpaired electrons are passed to
+        xtb, so radicals run open-shell automatically.
+    method : str
+        GFN parametrisation: ``"GFN2-xTB"`` (default), ``"GFN1-xTB"`` or
+        ``"GFN0-xTB"``.
+    solvent : str, optional
+        ALPB implicit-solvation solvent (e.g. ``"water"``). Native to xtb, so the
+        parameters are method-consistent. Defaults to gas phase.
+    directory : str, optional
+        Working directory (created if needed). Defaults to the current directory.
+    quasi_rrho : bool
+        If True, use Grimme's quasi-RRHO vibrational entropy. Default False.
+
+    Returns
+    -------
+    Thermo
+        The thermo calculation object.
+    """
+
+    if spin is None:
+        electrons = round(float(sum(atoms.get_atomic_numbers())) - charge)
+        spin = 0.0 if electrons % 2 == 0 else 0.5
+    unpaired = int(round(2.0 * float(spin)))
+
+    with _run_in_directory(directory):
+        optimized_atoms, potential_energy, frequencies = run_xtb(
+            atoms, charge=charge, unpaired=unpaired, method=method, solvent=solvent,
+        )
+
+        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, xtb_thermo
+from .api import dftbplus_thermo, xtb_thermo, xtb_cli_thermo
 
 logger = logging.getLogger(__package_name__).getChild("screening")
 logger = setup_logger(logger)
@@ -183,19 +183,21 @@ def screen(
         (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+.
+        Calculation engine: ``"dftb+"`` (default), ``"xtb"`` (GFN-xTB via tblite,
+        gas-phase) or ``"xtb-cli"`` (native ``xtb`` binary, which additionally
+        supports implicit solvation of charged radicals). ``parameter_set``
+        applies only to DFTB+; ``solvent`` applies to DFTB+ and xtb-cli.
     method : str
-        GFN-xTB parametrisation when ``engine="xtb"`` (``"GFN2-xTB"`` default).
+        GFN-xTB parametrisation for the xtb engines (``"GFN2-xTB"`` default).
 
     Returns
     -------
     list of dict
         One result record per molecule, including failed ones (status="error").
     """
-    if engine not in ("dftb+", "xtb"):
+    if engine not in ("dftb+", "xtb", "xtb-cli"):
         raise TSValueError(
-            f"Unknown engine {engine!r}; choose 'dftb+' or 'xtb'."
+            f"Unknown engine {engine!r}; choose 'dftb+', 'xtb' or 'xtb-cli'."
         )
 
     if engine == "dftb+":
@@ -217,7 +219,19 @@ def screen(
         logger.info(f"Screening {job.name} (charge {job.charge})")
         try:
             atoms = ase.io.read(str(job.path))
-            if engine == "xtb":
+            if engine == "xtb-cli":
+                thermo = xtb_cli_thermo(
+                    atoms,
+                    temperature=temperature,
+                    pressure=pressure,
+                    charge=job.charge,
+                    directory=str(root / job.name),
+                    spin=job.spin,
+                    method=method,
+                    solvent=solvent,
+                    quasi_rrho=quasi_rrho,
+                )
+            elif engine == "xtb":
                 thermo = xtb_thermo(
                     atoms,
                     temperature=temperature,