diff --git a/ThermoScreening/calculator/xtb.py b/ThermoScreening/calculator/xtb.py new file mode 100644 index 0000000..b72a77f --- /dev/null +++ b/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 diff --git a/ThermoScreening/cli/thermo.py b/ThermoScreening/cli/thermo.py index 5afe918..6ed6195 100644 --- a/ThermoScreening/cli/thermo.py +++ b/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") diff --git a/ThermoScreening/thermo/api.py b/ThermoScreening/thermo/api.py index 5ff315e..7732e2c 100644 --- a/ThermoScreening/thermo/api.py +++ b/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 diff --git a/ThermoScreening/thermo/screening.py b/ThermoScreening/thermo/screening.py index 6b06d4b..3ea444c 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 +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 diff --git a/ThermoScreening/thermo/thermo.py b/ThermoScreening/thermo/thermo.py index 35fa479..3565da0 100644 --- a/ThermoScreening/thermo/thermo.py +++ b/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.") diff --git a/tests/calculator/test_xtb.py b/tests/calculator/test_xtb.py new file mode 100644 index 0000000..d186772 --- /dev/null +++ b/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) diff --git a/tests/thermo/test_api.py b/tests/thermo/test_api.py index 418b43b..8f3291d 100644 --- a/tests/thermo/test_api.py +++ b/tests/thermo/test_api.py @@ -510,5 +510,49 @@ def test_dftbplus_thermo_forwards_quasi_rrho(monkeypatch, tmp_path): assert seen["run_thermo_quasi_rrho"] is True +def test_xtb_thermo_pipeline(monkeypatch, tmp_path): + import ThermoScreening.thermo.api as api + + seen = {} + + def fake_calculator(method): + seen["method"] = method + return "xtb-calc" + + def fake_optimise(atoms, calc, fmax=0.01): + seen["info"] = dict(atoms.info) + seen["calc"] = calc + return "optimized-atoms", -5.0, np.array([1500.0, 3600.0, 3700.0]) + + def fake_run_thermo(frequencies, atoms=None, engine=None, spin=None, + quasi_rrho=False, **kwargs): + seen["engine"] = engine + seen["spin"] = spin + seen["quasi_rrho"] = quasi_rrho + seen["energy"] = kwargs.get("energy") + return "thermo-result" + + monkeypatch.setattr(api, "xtb_calculator", fake_calculator) + monkeypatch.setattr(api, "optimise_and_frequencies", fake_optimise) + monkeypatch.setattr(api, "run_thermo", fake_run_thermo) + + # OH radical: 9 electrons -> auto doublet -> 1 unpaired electron for xTB + result = api.xtb_thermo( + Atoms("OH", positions=[[0, 0, 0], [0, 0, 0.97]]), + directory=str(tmp_path / "j"), + method="GFN1-xTB", + quasi_rrho=True, + ) + + assert result == "thermo-result" + assert seen["engine"] == "xtb" + assert seen["spin"] == 0.5 # auto electron-count guess + assert seen["info"]["spin"] == 1 # unpaired electrons passed to xTB + assert seen["info"]["charge"] == 0 + assert seen["method"] == "GFN1-xTB" + assert seen["quasi_rrho"] is True + assert seen["energy"] == -5.0 + + if __name__ == "__main__": unittest.main() diff --git a/tests/thermo/test_screening.py b/tests/thermo/test_screening.py index 5a9f189..7f1de08 100644 --- a/tests/thermo/test_screening.py +++ b/tests/thermo/test_screening.py @@ -146,6 +146,40 @@ def fake_thermo(atoms, solvent=None, **kwargs): assert captured["solvent"] == "water" +def test_screen_dispatches_to_xtb_engine(monkeypatch, tmp_path): + _write_xyz(tmp_path / "mol.xyz") + + captured = {} + + def fake_xtb(atoms, method="GFN2-xTB", **kwargs): + captured["method"] = method + captured["called"] = "xtb" + return _FakeThermo() + + def fake_dftb(*args, **kwargs): + captured["called"] = "dftb+" + return _FakeThermo() + + monkeypatch.setattr(screening, "xtb_thermo", fake_xtb) + monkeypatch.setattr(screening, "dftbplus_thermo", fake_dftb) + screening.screen( + str(tmp_path), + out=str(tmp_path / "r"), + directory=str(tmp_path / "runs"), + engine="xtb", + method="GFN1-xTB", + ) + + assert captured["called"] == "xtb" + assert captured["method"] == "GFN1-xTB" + + +def test_screen_rejects_unknown_engine(tmp_path): + _write_xyz(tmp_path / "mol.xyz") + with pytest.raises(TSValueError, match="Unknown engine"): + screening.screen(str(tmp_path), out=str(tmp_path / "r"), engine="orca") + + def test_screen_passes_quasi_rrho_to_dftbplus_thermo(monkeypatch, tmp_path): _write_xyz(tmp_path / "mol.xyz") @@ -277,6 +311,14 @@ def test_cli_parse_args_routes_screen(): qrrho_args = cli.parse_args(["screen", "molecules.csv", "--quasi-rrho"]) assert qrrho_args.quasi_rrho is True + assert args.engine == "dftb+" # default engine + assert args.method == "GFN2-xTB" + + xtb_args = cli.parse_args( + ["screen", "molecules.csv", "--engine", "xtb", "--method", "GFN1-xTB"] + ) + assert xtb_args.engine == "xtb" + assert xtb_args.method == "GFN1-xTB" def test_cli_run_screen_returns_failure_count(monkeypatch): @@ -289,7 +331,7 @@ def test_cli_run_screen_returns_failure_count(monkeypatch): args = Namespace( source="x", out="res", charge=0.0, temperature=298.15, pressure=101325.0, directory="screening", parameter_set="3ob", - solvent=None, quasi_rrho=False, + solvent=None, quasi_rrho=False, engine="dftb+", method="GFN2-xTB", ) assert cli.run_screen(args) == 1 # one molecule failed diff --git a/tests/thermo/test_thermo.py b/tests/thermo/test_thermo.py index 191ee1a..6ce591a 100644 --- a/tests/thermo/test_thermo.py +++ b/tests/thermo/test_thermo.py @@ -102,7 +102,15 @@ def _valid_system(): def test_thermo_rejects_unsupported_engine(): with pytest.raises(TSValueError, match="engine is not supported"): - Thermo(temperature=298.15, pressure=101325, system=_valid_system(), engine="xtb") + Thermo(temperature=298.15, pressure=101325, system=_valid_system(), engine="orca") + + +def test_thermo_accepts_xtb_engine(): + # xtb is a supported engine (GFN-xTB via tblite) + thermo = Thermo( + temperature=298.15, pressure=101325, system=_valid_system(), engine="xtb" + ) + assert thermo._engine == "xtb" def test_thermo_rejects_negative_temperature():