Add xtb_fukui_indices: per-atom reactivity via xtb --vfukui

ThermoScreening/calculator/xtb_cli.py

@@ -151,3 +151,109 @@ def run_xtb(
     energy_hartree = _parse_optimised_energy("xtbopt.xyz")
     frequencies = _parse_vibspectrum("vibspectrum")
     return optimized, energy_hartree, frequencies
+
+
+def _parse_fukui(output):
+    """
+    Parse the per-atom ``Fukui functions:`` table from ``xtb --vfukui`` stdout.
+
+    The table looks like::
+
+        Fukui functions:
+             #        f(+)     f(-)     f(0)
+             1C       0.024    0.024    0.024
+             8O       0.131    0.129    0.130
+
+    Returns a list of ``(symbol, f_plus, f_minus, f_zero)`` tuples, in the
+    input geometry's atom order.
+
+    Raises
+    ------
+    ValueError
+        If no ``Fukui functions:`` table is found in ``output``.
+    """
+    lines = output.splitlines()
+    for index, line in enumerate(lines):
+        if line.strip() == "Fukui functions:":
+            break
+    else:
+        raise ValueError("No 'Fukui functions:' table found in xtb output.")
+
+    rows = []
+    for line in lines[index + 2:]:
+        tokens = line.split()
+        if len(tokens) != 4:
+            break
+        label, f_plus, f_minus, f_zero = tokens
+        symbol = label.lstrip("0123456789")
+        rows.append((symbol, float(f_plus), float(f_minus), float(f_zero)))
+    return rows
+
+
+def run_xtb_fukui(
+    atoms,
+    charge=0.0,
+    unpaired=0,
+    method="GFN2-xTB",
+    solvent=None,
+    command=None,
+):
+    """
+    Run ``xtb --vfukui`` and return per-atom Fukui reactivity indices.
+
+    A single-point (non-geometry-optimising) calculation: pass an already
+    optimised geometry (e.g. from :func:`run_xtb`). Runs in the current
+    working directory, so call it inside a per-job directory.
+
+    Parameters
+    ----------
+    atoms : ase.Atoms
+        The geometry to analyse.
+    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
+    -------
+    list of tuple(str, float, float, float)
+        ``(symbol, f_plus, f_minus, f_zero)`` per atom, in input geometry order.
+        ``f_plus`` (susceptibility to nucleophilic attack / electron gain),
+        ``f_minus`` (electrophilic attack / electron loss), ``f_zero`` (radical
+        attack, the average of the two).
+
+    Raises
+    ------
+    ValueError
+        If ``method`` is unknown, or xtb's output has no Fukui table.
+    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", "--vfukui", "--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:
+        tail = (result.stdout or "")[-800:] + (result.stderr or "")[-800:]
+        raise RuntimeError(f"xtb failed (exit {result.returncode}):\n{tail}")
+
+    return _parse_fukui(result.stdout)

ThermoScreening/thermo/api.py

@@ -22,7 +22,7 @@
 from ..calculator.orca import read_orca_hess
 from ..calculator.qm import read_cclib, _describe_source
 from ..calculator.xtb import optimise_and_frequencies, xtb_calculator
-from ..calculator.xtb_cli import run_xtb
+from ..calculator.xtb_cli import run_xtb, run_xtb_fukui
 
 
 logger = logging.getLogger(__package_name__).getChild("api")
@@ -1074,3 +1074,59 @@ def xtb_cli_thermo(
         )
 
     return thermo
+
+
+def xtb_fukui_indices(atoms, charge=0.0, spin=None, method="GFN2-xTB", solvent=None, directory=None):
+    """
+    Per-atom Fukui reactivity indices via the native ``xtb`` program.
+
+    A single-point analysis (no geometry optimisation) -- pass an already
+    optimised geometry, e.g. from :func:`xtb_cli_thermo`. Each index estimates
+    how much electron density a site would gain or lose if the whole molecule
+    were reduced or oxidised by one electron, which is a direct, cheap way to
+    localise *where* a computed reduction/oxidation is likely to happen (or
+    where an unexpected side reaction might), rather than relying on whole-
+    molecule quantities like ``total_EeGtot()`` or frontier orbital energies.
+
+    Parameters
+    ----------
+    atoms : ase.Atoms
+        The (already optimised) geometry to analyse.
+    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), matching :func:`xtb_cli_thermo`.
+    method : str
+        GFN parametrisation: ``"GFN2-xTB"`` (default), ``"GFN1-xTB"`` or
+        ``"GFN0-xTB"``.
+    solvent : str, optional
+        ALPB implicit-solvation solvent (e.g. ``"water"``). Defaults to gas phase.
+    directory : str, optional
+        Working directory (created if needed). Defaults to the current directory.
+
+    Returns
+    -------
+    list of tuple(str, float, float, float)
+        ``(symbol, f_plus, f_minus, f_zero)`` per atom, in ``atoms`` order.
+        ``f_plus`` is the susceptibility to nucleophilic attack (electron gain,
+        i.e. reduction) at that site, ``f_minus`` to electrophilic attack
+        (electron loss, oxidation), and ``f_zero`` (their average) to radical
+        attack.
+
+    Raises
+    ------
+    ValueError
+        If ``method`` is unknown, or xtb's output has no Fukui table.
+    RuntimeError
+        If the xtb run fails.
+    """
+    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):
+        return run_xtb_fukui(
+            atoms, charge=charge, unpaired=unpaired, method=method, solvent=solvent,
+        )

docs/api.rst

@@ -7,6 +7,7 @@ High-level pipeline
 .. autofunction:: ThermoScreening.thermo.api.dftbplus_thermo
 .. autofunction:: ThermoScreening.thermo.api.xtb_thermo
 .. autofunction:: ThermoScreening.thermo.api.xtb_cli_thermo
+.. autofunction:: ThermoScreening.thermo.api.xtb_fukui_indices
 .. autofunction:: ThermoScreening.thermo.api.orca_thermo
 .. autofunction:: ThermoScreening.thermo.api.cclib_thermo
 .. autofunction:: ThermoScreening.thermo.api.pyscf_thermo

docs/usage.rst

@@ -259,6 +259,35 @@ conditions), combine them into reaction free energies and reduction potentials:
    across similar species, with a higher-accuracy method, or with a
    ``reference_potential`` calibrated against experiment.
 
+Reactivity site prediction (Fukui indices)
+-------------------------------------------
+
+``reaction_free_energy`` and ``reduction_potential`` tell you *whether* a
+reaction is favourable for the molecule as a whole; they say nothing about
+*where* on the molecule it happens, or whether some other site might react
+first (an unintended side reaction competing with the one you're modelling).
+``xtb_fukui_indices`` runs a single ``xtb --vfukui`` calculation on an
+already-optimised geometry and reports, per atom, how susceptible that site
+is to gaining electron density (``f_plus``, i.e. reduction), losing it
+(``f_minus``, oxidation), or either (``f_zero``, radical attack):
+
+.. code-block:: python
+
+   import ase.io
+   from ThermoScreening.thermo.api import xtb_cli_thermo, xtb_fukui_indices
+   from ThermoScreening.thermo.conformers import generate
+
+   mol = generate("O=C1C=CC(=O)C=C1", max_conformers=1)[0]   # benzoquinone
+   xtb_cli_thermo(mol, directory="opt")   # optimise; xtbopt.xyz is written to "opt"
+
+   optimized = ase.io.read("opt/xtbopt.xyz")
+   rows = xtb_fukui_indices(optimized)
+   most_oxidisable = max(rows, key=lambda row: row[2])   # (symbol, f_plus, f_minus, f_zero)
+
+Each of ``f_plus``/``f_minus``/``f_zero`` sums to ~1 over all atoms, so a
+value well above the 1/n_atoms baseline flags that site as the dominant one
+for that kind of reactivity.
+
 Acid dissociation (pKa)
 ------------------------
 

tests/calculator/test_xtb_cli.py

@@ -89,6 +89,70 @@ def fake_run(argv, **kwargs):
     assert argv[argv.index("--alpb") + 1] == "water"
 
 
+def test_parse_fukui_extracts_table():
+    output = (
+        "some preceding SCF chatter\n"
+        "\n"
+        "Fukui functions:\n"
+        "     #        f(+)     f(-)     f(0)\n"
+        "     1C       0.024    0.024    0.024\n"
+        "     8O       0.131    0.129    0.130\n"
+        "\n"
+        "trailing text after the table\n"
+    )
+    rows = xtb_cli._parse_fukui(output)
+
+    assert rows == [("C", 0.024, 0.024, 0.024), ("O", 0.131, 0.129, 0.130)]
+
+
+def test_parse_fukui_raises_without_table():
+    with pytest.raises(ValueError, match="No 'Fukui functions:' table"):
+        xtb_cli._parse_fukui("no such table here")
+
+
+def test_run_xtb_fukui_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
+        stdout = "Fukui functions:\n     #        f(+)     f(-)     f(0)\n     1O       0.5      0.5      0.5\n"
+        return type("R", (), {"returncode": 0, "stdout": stdout, "stderr": ""})()
+
+    monkeypatch.setattr(xtb_cli.subprocess, "run", fake_run)
+
+    rows = xtb_cli.run_xtb_fukui(
+        Atoms("OH", positions=[[0, 0, 0], [0, 0, 0.97]]),
+        charge=-1, unpaired=1, method="GFN2-xTB", solvent="water",
+    )
+
+    assert rows == [("O", 0.5, 0.5, 0.5)]
+    argv = captured["argv"]
+    assert "--vfukui" in argv
+    assert "--ohess" not 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_fukui_rejects_unknown_method(monkeypatch):
+    with pytest.raises(ValueError, match="Unknown xTB method"):
+        xtb_cli.run_xtb_fukui(Atoms("H2", positions=[[0, 0, 0], [0, 0, 0.7]]), method="B3LYP")
+
+
+def test_run_xtb_fukui_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_fukui(Atoms("H2", positions=[[0, 0, 0], [0, 0, 0.7]]))
+
+
 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"):
@@ -125,3 +189,29 @@ def oh():
     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()
+
+
+@pytest.mark.skipif(not xtb_available, reason="the native xtb binary is not available.")
+def test_xtb_fukui_indices_end_to_end(tmp_path):
+    import ase.io
+    from ThermoScreening.thermo.api import xtb_cli_thermo, xtb_fukui_indices
+
+    water = Atoms("OH2", positions=[[0, 0, 0.119], [0, 0.763, -0.477], [0, -0.763, -0.477]])
+    opt_dir = tmp_path / "opt"
+    xtb_cli_thermo(water, directory=str(opt_dir))
+    optimized = ase.io.read(str(opt_dir / "xtbopt.xyz"))
+
+    rows = xtb_fukui_indices(optimized, directory=str(tmp_path / "fukui"))
+
+    assert len(rows) == 3
+    symbols = [symbol for symbol, *_ in rows]
+    assert symbols == ["O", "H", "H"]
+    # Fukui functions are normalised: each column sums to ~1 over all atoms
+    assert sum(f_plus for _, f_plus, _, _ in rows) == pytest.approx(1.0, abs=0.05)
+    assert sum(f_minus for _, _, f_minus, _ in rows) == pytest.approx(1.0, abs=0.05)
+    # the oxygen lone pairs make O the dominant site for electron loss (f_minus,
+    # i.e. oxidation); the two symmetric H atoms dominate electron gain (f_plus)
+    oxygen_row = rows[0]
+    hydrogen_rows = rows[1:]
+    assert oxygen_row[2] > max(row[2] for row in hydrogen_rows)
+    assert oxygen_row[1] < min(row[1] for row in hydrogen_rows)

tests/thermo/test_api.py

@@ -642,5 +642,47 @@ def fake_run_thermo(frequencies, atoms=None, engine=None, spin=None,
     assert seen["energy"] == -4.5
 
 
+def test_xtb_fukui_indices_pipeline(monkeypatch, tmp_path):
+    import ThermoScreening.thermo.api as api
+
+    seen = {}
+
+    def fake_run_xtb_fukui(atoms, charge=0.0, unpaired=0, method="GFN2-xTB", solvent=None):
+        seen.update(charge=charge, unpaired=unpaired, method=method, solvent=solvent)
+        return [("O", 0.5, 0.5, 0.5), ("H", 0.25, 0.25, 0.25), ("H", 0.25, 0.25, 0.25)]
+
+    monkeypatch.setattr(api, "run_xtb_fukui", fake_run_xtb_fukui)
+
+    # neutral OH radical -> auto doublet -> 1 unpaired electron
+    result = api.xtb_fukui_indices(
+        Atoms("OH", positions=[[0, 0, 0], [0, 0, 0.97]]),
+        directory=str(tmp_path / "j"),
+        solvent="water",
+        method="GFN1-xTB",
+    )
+
+    assert result == [("O", 0.5, 0.5, 0.5), ("H", 0.25, 0.25, 0.25), ("H", 0.25, 0.25, 0.25)]
+    assert seen["unpaired"] == 1  # auto electron-count guess -> doublet -> 1 unpaired
+    assert seen["charge"] == 0.0
+    assert seen["solvent"] == "water"
+    assert seen["method"] == "GFN1-xTB"
+
+
+def test_xtb_fukui_indices_explicit_spin(monkeypatch, tmp_path):
+    import ThermoScreening.thermo.api as api
+
+    seen = {}
+
+    def fake_run_xtb_fukui(atoms, charge=0.0, unpaired=0, method="GFN2-xTB", solvent=None):
+        seen["unpaired"] = unpaired
+        return []
+
+    monkeypatch.setattr(api, "run_xtb_fukui", fake_run_xtb_fukui)
+
+    api.xtb_fukui_indices(Atoms("H2", positions=[[0, 0, 0], [0, 0, 0.74]]), charge=-2.0, spin=1.0)
+
+    assert seen["unpaired"] == 2
+
+
 if __name__ == "__main__":
     unittest.main()