From 2c9630dd538efe83d943a67a25ce09db8f17c4b8 Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Mon, 6 Jul 2026 00:05:49 +0200 Subject: [PATCH] Add a native-xtb engine with implicit solvation for open-shell species The tblite xTB engine (in-process) can't do solvation (tblite's ASE calculator exposes no ALPB/GBSA), and running xTB through DFTB+ can't do open-shell (a known DFTB+ limitation). The native xtb binary does all of it, so add it as the "xtb-cli" engine -- covering charged radicals in solution with method-consistent parameters. - calculator/xtb_cli.py: run `xtb --ohess --gfn --chrg --uhf <2S> [--alpb ]`, then parse the energy (xtbopt.xyz), optimised geometry, and frequencies (vibspectrum), feeding them into run_thermo like the other engines. The xtb executable is found via $XTB_COMMAND or PATH; nothing is imported from xtb, so the package still imports without it. - thermo/api.py: `xtb_cli_thermo`, resolving spin to --uhf so radicals run open-shell automatically; solvent -> --alpb. - thermo/screening.py, cli/thermo.py: `screen(engine="xtb-cli")` / `--engine xtb-cli`; solvent now applies to dftb+ and xtb-cli. Validated against the real xtb 6.7.1: gas-phase water S = 45.05 cal/mol/K (experiment ~45.1); an OH radical runs open-shell and is stabilised by ~6.5 kcal/mol in water (--alpb). Frequencies/energy match the tblite path. --- ThermoScreening/calculator/xtb_cli.py | 153 ++++++++++++++++++++++++++ ThermoScreening/cli/thermo.py | 8 +- ThermoScreening/thermo/api.py | 77 +++++++++++++ ThermoScreening/thermo/screening.py | 28 +++-- tests/calculator/test_xtb_cli.py | 127 +++++++++++++++++++++ tests/thermo/test_api.py | 38 +++++++ tests/thermo/test_screening.py | 25 +++++ 7 files changed, 446 insertions(+), 10 deletions(-) create mode 100644 ThermoScreening/calculator/xtb_cli.py create mode 100644 tests/calculator/test_xtb_cli.py diff --git a/ThermoScreening/calculator/xtb_cli.py b/ThermoScreening/calculator/xtb_cli.py new file mode 100644 index 0000000..f79d269 --- /dev/null +++ b/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 diff --git a/ThermoScreening/cli/thermo.py b/ThermoScreening/cli/thermo.py index 6ed6195..31f6441 100644 --- a/ThermoScreening/cli/thermo.py +++ b/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 diff --git a/ThermoScreening/thermo/api.py b/ThermoScreening/thermo/api.py index 7732e2c..74b40d3 100644 --- a/ThermoScreening/thermo/api.py +++ b/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 diff --git a/ThermoScreening/thermo/screening.py b/ThermoScreening/thermo/screening.py index 3ea444c..25b5fdb 100644 --- a/ThermoScreening/thermo/screening.py +++ b/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, diff --git a/tests/calculator/test_xtb_cli.py b/tests/calculator/test_xtb_cli.py new file mode 100644 index 0000000..e9afc03 --- /dev/null +++ b/tests/calculator/test_xtb_cli.py @@ -0,0 +1,127 @@ +import os +import shutil + +import numpy as np +import pytest +from ase import Atoms + +from ThermoScreening.calculator import xtb_cli + + +# --- executable resolution --------------------------------------------------- # + +def test_resolve_xtb_prefers_explicit_then_env(monkeypatch): + monkeypatch.setenv("XTB_COMMAND", "/opt/xtb") + assert xtb_cli.resolve_xtb("/usr/bin/xtb") == "/usr/bin/xtb" # explicit wins + assert xtb_cli.resolve_xtb() == "/opt/xtb" # then env + + +def test_resolve_xtb_missing_raises(monkeypatch): + monkeypatch.delenv("XTB_COMMAND", raising=False) + monkeypatch.setattr(xtb_cli.shutil, "which", lambda command: None) + with pytest.raises(FileNotFoundError, match="xtb"): + xtb_cli.resolve_xtb() + + +# --- output parsing ---------------------------------------------------------- # + +def test_parse_vibspectrum_takes_wavenumber_and_sorts(tmp_path): + vib = tmp_path / "vibspectrum" + vib.write_text( + "$vibrational spectrum\n" + "# mode symmetry wave number IR intensity selection rules\n" + " 1 -0.00 0.00000 - \n" + " 6 0.00 0.00000 - \n" + " 7 a 1541.16 133.03481 YES\n" + " 8 a 3635.99 6.71206 YES\n" + "$end\n", + encoding="utf-8", + ) + out = xtb_cli._parse_vibspectrum(str(vib)) + + assert len(out) == 4 + assert np.all(np.diff(out) >= 0) # sorted ascending + assert list(out[-2:]) == [1541.16, 3635.99] # real modes on top + + +def test_parse_optimised_energy(tmp_path): + xyz = tmp_path / "xtbopt.xyz" + xyz.write_text( + "2\n energy: -4.438237080538 gnorm: 0.0004 xtb: 6.7.1\n" + "O 0.0 0.0 0.0\nH 0.0 0.0 0.97\n", + encoding="utf-8", + ) + assert xtb_cli._parse_optimised_energy(str(xyz)) == pytest.approx(-4.438237080538) + + +# --- run_xtb (mocked subprocess, no binary needed) --------------------------- # + +def test_run_xtb_builds_command_and_parses(monkeypatch, tmp_path): + monkeypatch.chdir(tmp_path) + monkeypatch.setenv("XTB_COMMAND", "/fake/xtb") + captured = {} + + def fake_run(argv, **kwargs): + captured["argv"] = argv + (tmp_path / "xtbopt.xyz").write_text( + "2\n energy: -4.4 gnorm: 0 xtb: x\nO 0 0 0\nH 0 0 0.97\n", encoding="utf-8" + ) + (tmp_path / "vibspectrum").write_text( + "$vibrational spectrum\n 1 -0.00 0.0 -\n" + " 6 a 3500.0 1.0 YES\n$end\n", encoding="utf-8" + ) + return type("R", (), {"returncode": 0, "stdout": "", "stderr": ""})() + + monkeypatch.setattr(xtb_cli.subprocess, "run", fake_run) + + atoms, energy, freqs = xtb_cli.run_xtb( + Atoms("OH", positions=[[0, 0, 0], [0, 0, 0.97]]), + charge=-1, unpaired=1, method="GFN2-xTB", solvent="water", + ) + + assert energy == pytest.approx(-4.4) + assert freqs[-1] == pytest.approx(3500.0) + argv = captured["argv"] + assert "--ohess" in argv + assert argv[argv.index("--gfn") + 1] == "2" + assert argv[argv.index("--uhf") + 1] == "1" + assert argv[argv.index("--chrg") + 1] == "-1" + assert argv[argv.index("--alpb") + 1] == "water" + + +def test_run_xtb_rejects_unknown_method(monkeypatch): + # method check happens before the executable is even resolved + with pytest.raises(ValueError, match="Unknown xTB method"): + xtb_cli.run_xtb(Atoms("H2", positions=[[0, 0, 0], [0, 0, 0.7]]), method="B3LYP") + + +def test_run_xtb_raises_on_failure(monkeypatch, tmp_path): + monkeypatch.chdir(tmp_path) + monkeypatch.setenv("XTB_COMMAND", "/fake/xtb") + monkeypatch.setattr( + xtb_cli.subprocess, "run", + lambda argv, **kw: type("R", (), {"returncode": 1, "stdout": "boom", "stderr": "e"})(), + ) + with pytest.raises(RuntimeError, match="xtb failed"): + xtb_cli.run_xtb(Atoms("H2", positions=[[0, 0, 0], [0, 0, 0.7]])) + + +# --- real integration (needs the xtb binary) -------------------------------- # + +xtb_available = shutil.which("xtb") is not None or "XTB_COMMAND" in os.environ + + +@pytest.mark.skipif(not xtb_available, reason="the native xtb binary is not available.") +def test_xtb_cli_radical_solvation_end_to_end(tmp_path): + from ThermoScreening.thermo.api import xtb_cli_thermo + + def oh(): + return Atoms("OH", positions=[[0, 0, 0], [0, 0, 0.97]]) + + gas = xtb_cli_thermo(oh(), directory=str(tmp_path / "gas")) + solvated = xtb_cli_thermo(oh(), directory=str(tmp_path / "sol"), solvent="water") + + # open-shell radical (doublet: S ~ 42-43 cal/mol/K incl. R ln2) + assert gas.total_entropy("cal/(mol*K)") == pytest.approx(42.5, abs=2.0) + # solvation stabilises the radical + assert solvated.total_EeGtot() < gas.total_EeGtot() diff --git a/tests/thermo/test_api.py b/tests/thermo/test_api.py index 8f3291d..bd4f6fd 100644 --- a/tests/thermo/test_api.py +++ b/tests/thermo/test_api.py @@ -554,5 +554,43 @@ def fake_run_thermo(frequencies, atoms=None, engine=None, spin=None, assert seen["energy"] == -5.0 +def test_xtb_cli_thermo_pipeline(monkeypatch, tmp_path): + import ThermoScreening.thermo.api as api + + seen = {} + + def fake_run_xtb(atoms, charge=0.0, unpaired=0, method="GFN2-xTB", solvent=None): + seen.update(charge=charge, unpaired=unpaired, method=method, solvent=solvent) + return "optimized-atoms", -4.5, np.array([1500.0, 3600.0, 3700.0]) + + def fake_run_thermo(frequencies, atoms=None, engine=None, spin=None, + quasi_rrho=False, **kwargs): + seen.update(engine=engine, spin=spin, quasi_rrho=quasi_rrho, + energy=kwargs.get("energy")) + return "thermo-result" + + monkeypatch.setattr(api, "run_xtb", fake_run_xtb) + monkeypatch.setattr(api, "run_thermo", fake_run_thermo) + + # neutral OH radical -> auto doublet -> 1 unpaired electron; water solvation + result = api.xtb_cli_thermo( + Atoms("OH", positions=[[0, 0, 0], [0, 0, 0.97]]), + directory=str(tmp_path / "j"), + solvent="water", + method="GFN1-xTB", + quasi_rrho=True, + ) + + assert result == "thermo-result" + assert seen["engine"] == "xtb" + assert seen["spin"] == 0.5 # auto electron-count guess (9 electrons) + assert seen["unpaired"] == 1 # round(2*S) passed to xtb --uhf + assert seen["charge"] == 0.0 + assert seen["solvent"] == "water" + assert seen["method"] == "GFN1-xTB" + assert seen["quasi_rrho"] is True + assert seen["energy"] == -4.5 + + if __name__ == "__main__": unittest.main() diff --git a/tests/thermo/test_screening.py b/tests/thermo/test_screening.py index 7f1de08..5d62281 100644 --- a/tests/thermo/test_screening.py +++ b/tests/thermo/test_screening.py @@ -174,6 +174,28 @@ def fake_dftb(*args, **kwargs): assert captured["method"] == "GFN1-xTB" +def test_screen_dispatches_to_xtb_cli_engine(monkeypatch, tmp_path): + _write_xyz(tmp_path / "mol.xyz") + + captured = {} + + def fake_xtb_cli(atoms, method="GFN2-xTB", solvent=None, **kwargs): + captured.update(called="xtb-cli", method=method, solvent=solvent) + return _FakeThermo() + + monkeypatch.setattr(screening, "xtb_cli_thermo", fake_xtb_cli) + screening.screen( + str(tmp_path), + out=str(tmp_path / "r"), + directory=str(tmp_path / "runs"), + engine="xtb-cli", + solvent="water", + ) + + assert captured["called"] == "xtb-cli" + assert captured["solvent"] == "water" # xtb-cli supports solvation + + def test_screen_rejects_unknown_engine(tmp_path): _write_xyz(tmp_path / "mol.xyz") with pytest.raises(TSValueError, match="Unknown engine"): @@ -320,6 +342,9 @@ def test_cli_parse_args_routes_screen(): assert xtb_args.engine == "xtb" assert xtb_args.method == "GFN1-xTB" + cli_args = cli.parse_args(["screen", "molecules.csv", "--engine", "xtb-cli"]) + assert cli_args.engine == "xtb-cli" + def test_cli_run_screen_returns_failure_count(monkeypatch): import ThermoScreening.cli.thermo as cli