From 62fdee17e37ecca865469019becac2d90b4bfff1 Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Thu, 2 Jul 2026 12:48:22 +0200 Subject: [PATCH 1/2] Run open-shell species spin-polarised in DFTB+ automatically Follow-up to #63: make the spin actually drive the DFTB+ calculation, not just the electronic-degeneracy entropy. dftbplus_thermo now resolves the spin (the electron-count guess when not given) and, when S > 0, runs Geoopt/Hessian colinear spin-polarised via _spin_kwargs: SpinPolarisation = Colinear with UnpairedElectrons = round(2S) and the 3ob-PBE SpinConstants injected for the elements present. Closed-shell (S = 0) keeps the restricted calculation unchanged. So a radical (odd electrons) is auto-detected and treated open-shell with no extra input; only even-electron triplets (O2) need an explicit spin. The 3ob spin constants are cross-verified (3ob Elstner list vs QUASINANO2015 vs DFTB+ manual/mio); elements without one (Br, I, metals) raise a clear error. Energies should be validated with a real spin-polarised run. --- ThermoScreening/calculator/dftbplus.py | 66 +++++++++++++++++++++ ThermoScreening/thermo/api.py | 19 ++++++- tests/calculator/test_dftbplus.py | 25 ++++++++ tests/thermo/test_api.py | 79 +++++++++++++++++++++++++- 4 files changed, 186 insertions(+), 3 deletions(-) diff --git a/ThermoScreening/calculator/dftbplus.py b/ThermoScreening/calculator/dftbplus.py index e08f22f..97665ac 100644 --- a/ThermoScreening/calculator/dftbplus.py +++ b/ThermoScreening/calculator/dftbplus.py @@ -51,6 +51,72 @@ def _slako_dir(slako_dir=None): return selected_dir + os.sep +# Atomic spin constants (Hartree) for the 3ob-3-1 PBE Slater-Koster set, i.e. the +# single spin constant of the highest occupied shell per element (Wss for H, Wpp +# otherwise), used with ShellResolvedSpin = No to match this tool's atom-resolved +# 3ob SCC. These are parameters that must match the SK set's functional (PBE) and +# parametrisation; using a different SK set requires its own spin constants. +# Cross-verified: 3ob (Elstner list) vs QUASINANO2015 (arXiv:1605.01360) vs the +# DFTB+ manual/mio (Koehler). Validate energies with a real spin-polarised run. +SPIN_CONSTANTS = { + "H": "{ -0.07174 }", + "C": "{ -0.02265 }", + "N": "{ -0.02545 }", + "O": "{ -0.02785 }", + "F": "{ -0.02990 }", + "P": "{ -0.01490 }", + "S": "{ -0.01549 }", + "Cl": "{ -0.01606 }", +} + + +def _spin_kwargs(atoms, spin): + """ + ASE ``Dftb`` keyword arguments enabling colinear spin polarisation. + + Returns an empty dict for ``spin`` in (None, 0), so the closed-shell + (restricted) calculation is left exactly as before. For spin S > 0 it enables + ``SpinPolarisation = Colinear`` with ``UnpairedElectrons = round(2*S)`` and + injects the 3ob ``SpinConstants`` for the elements present. + + Parameters + ---------- + atoms : ase.Atoms + The atoms whose elements need spin constants. + spin : float or None + Spin quantum number S. + + Raises + ------ + ValueError + If an element has no tabulated 3ob spin constant. + """ + if spin is None or float(spin) <= 0.0: + return {} + + unpaired = int(round(2.0 * float(spin))) + if unpaired <= 0: + return {} + + elements = sorted(set(atoms.get_chemical_symbols())) + missing = [element for element in elements if element not in SPIN_CONSTANTS] + if missing: + raise ValueError( + "Spin-polarised DFTB+ is not available for element(s) " + f"{', '.join(missing)}: no 3ob spin constant is tabulated." + ) + + kwargs = { + "Hamiltonian_SpinPolarisation": "Colinear {", + "Hamiltonian_SpinPolarisation_UnpairedElectrons": unpaired, + "Hamiltonian_SpinConstants_": "", + "Hamiltonian_SpinConstants_ShellResolvedSpin": "No", + } + for element in elements: + kwargs[f"Hamiltonian_SpinConstants_{element}"] = SPIN_CONSTANTS[element] + return kwargs + + class Geoopt(Dftb): """ Custom DFTB+ calculator to optimize the system with the 'GeometryOptimisation' driver (Rational). diff --git a/ThermoScreening/thermo/api.py b/ThermoScreening/thermo/api.py index e4f896b..f015192 100644 --- a/ThermoScreening/thermo/api.py +++ b/ThermoScreening/thermo/api.py @@ -18,6 +18,7 @@ from .thermo import Thermo from .atoms import Atom from ..calculator import Geoopt, Hessian, Modes +from ..calculator.dftbplus import _spin_kwargs logger = logging.getLogger(__package_name__).getChild("api") @@ -523,6 +524,12 @@ def dftbplus_thermo( geometry/Hessian/modes files are written here, so separate jobs can run without clobbering each other. Defaults to the current directory. + spin : float, optional + Spin quantum number S. Defaults to the minimum-spin electron-count guess + (even -> 0, odd -> 0.5). When S > 0 the DFTB+ steps run colinear + spin-polarised (so radicals are treated open-shell automatically); S = 0 + keeps the restricted closed-shell calculation. + Other Parameters ---------------- **kwargs : dict @@ -534,14 +541,22 @@ def dftbplus_thermo( The thermo calculation object. """ + # Resolve the spin (electron-count guess when not given) so the calculation + # and the analysis use the same multiplicity; S > 0 enables spin polarisation. + if spin is None: + electrons = round(float(sum(atoms.get_atomic_numbers())) - charge) + spin = 0.0 if electrons % 2 == 0 else 0.5 + + spin_kwargs = _spin_kwargs(atoms, spin) + with _run_in_directory(directory): # run geometry optimization - geoopt = Geoopt(atoms=atoms, charge=charge, **kwargs) + geoopt = Geoopt(atoms=atoms, charge=charge, **spin_kwargs, **kwargs) potential_energy = geoopt.potential_energy() optimized_atoms = geoopt.read() # run hessian calculation - Hessian(atoms=optimized_atoms, charge=charge, **kwargs) + Hessian(atoms=optimized_atoms, charge=charge, **spin_kwargs, **kwargs) # run normal mode calculation modes = Modes() diff --git a/tests/calculator/test_dftbplus.py b/tests/calculator/test_dftbplus.py index 77acb2b..778756f 100644 --- a/tests/calculator/test_dftbplus.py +++ b/tests/calculator/test_dftbplus.py @@ -390,6 +390,31 @@ def get_global_number_of_atoms(self): assert hessian.read().shape == (9, 9) +def test_spin_kwargs_restricted_and_fractional(): + h2 = Atoms("H2", positions=[[0, 0, 0], [0, 0, 0.74]]) + assert dftbplus_module._spin_kwargs(h2, None) == {} + assert dftbplus_module._spin_kwargs(h2, 0.0) == {} + assert dftbplus_module._spin_kwargs(h2, 0.1) == {} # rounds to 0 unpaired + + +def test_spin_kwargs_builds_colinear_block(): + kw = dftbplus_module._spin_kwargs( + Atoms("OH", positions=[[0, 0, 0], [0.97, 0, 0]]), 0.5 + ) + assert kw["Hamiltonian_SpinPolarisation"] == "Colinear {" + assert kw["Hamiltonian_SpinPolarisation_UnpairedElectrons"] == 1 + assert kw["Hamiltonian_SpinConstants_ShellResolvedSpin"] == "No" + assert kw["Hamiltonian_SpinConstants_H"] == "{ -0.07174 }" + assert kw["Hamiltonian_SpinConstants_O"] == "{ -0.02785 }" + + +def test_spin_kwargs_rejects_element_without_constant(): + with pytest.raises(ValueError, match="not available for element"): + dftbplus_module._spin_kwargs( + Atoms("Br2", positions=[[0, 0, 0], [0, 0, 2.3]]), 0.5 + ) + + @pytest.mark.skipif( not dftbplus_ready, reason="DFTB+ executables are not installed or cannot start.", diff --git a/tests/thermo/test_api.py b/tests/thermo/test_api.py index bb56157..b8cbd45 100644 --- a/tests/thermo/test_api.py +++ b/tests/thermo/test_api.py @@ -21,6 +21,7 @@ unit_frequency, ) from ase.build import molecule +from ase import Atoms from ThermoScreening.exceptions import TSNotImplementedError, TSValueError class TestApi(unittest.TestCase): @@ -345,6 +346,7 @@ def test_dftbplus_thermo_runs_pipeline_in_directory(monkeypatch, tmp_path): class FakeGeoopt: def __init__(self, atoms, charge, **kwargs): seen["cwd_during"] = Path.cwd().resolve() + seen["geoopt_kwargs"] = kwargs def potential_energy(self): return -1.0 @@ -371,13 +373,88 @@ def fake_run_thermo(frequencies, atoms=None, **kwargs): job = tmp_path / "job" start = Path.cwd() - result = api.dftbplus_thermo("initial-atoms", directory=str(job)) + water = Atoms("OH2", positions=[[0, 0, 0.12], [0, 0.76, -0.48], [0, -0.76, -0.48]]) + result = api.dftbplus_thermo(water, directory=str(job)) assert result == "thermo-result" assert seen["run_thermo_atoms"] == "optimized-atoms" assert seen["hessian_atoms"] == "optimized-atoms" assert seen["cwd_during"] == job.resolve() # pipeline ran inside the job dir assert Path.cwd() == start # working directory restored afterwards + # closed-shell water -> restricted, no spin polarisation in the DFTB+ kwargs + assert "Hamiltonian_SpinPolarisation" not in seen["geoopt_kwargs"] + + +def _mock_pipeline(monkeypatch): + import ThermoScreening.thermo.api as api + + seen = {} + + class FakeGeoopt: + def __init__(self, atoms, charge, **kwargs): + seen["geoopt_kwargs"] = kwargs + + def potential_energy(self): + return -1.0 + + def read(self): + return "optimized-atoms" + + class FakeHessian: + def __init__(self, atoms, charge, **kwargs): + seen["hessian_kwargs"] = kwargs + + class FakeModes: + def __init__(self): + self.wave_numbers = np.array([1.0, 2.0, 3.0]) + + def fake_run_thermo(frequencies, atoms=None, spin=None, **kwargs): + seen["run_thermo_spin"] = spin + return "thermo-result" + + monkeypatch.setattr(api, "Geoopt", FakeGeoopt) + monkeypatch.setattr(api, "Hessian", FakeHessian) + monkeypatch.setattr(api, "Modes", FakeModes) + monkeypatch.setattr(api, "run_thermo", fake_run_thermo) + return api, seen + + +def test_dftbplus_thermo_spin_polarises_radical_automatically(monkeypatch, tmp_path): + api, seen = _mock_pipeline(monkeypatch) + # OH: 9 electrons (odd) -> auto doublet -> spin-polarised, no user input + api.dftbplus_thermo( + Atoms("OH", positions=[[0, 0, 0], [0.97, 0, 0]]), directory=str(tmp_path / "j") + ) + + assert seen["run_thermo_spin"] == 0.5 # analysis multiplicity matches + assert seen["geoopt_kwargs"]["Hamiltonian_SpinPolarisation"] == "Colinear {" + assert seen["geoopt_kwargs"]["Hamiltonian_SpinPolarisation_UnpairedElectrons"] == 1 + assert seen["hessian_kwargs"]["Hamiltonian_SpinConstants_O"] == "{ -0.02785 }" + + +def test_dftbplus_thermo_restricted_for_closed_shell(monkeypatch, tmp_path): + api, seen = _mock_pipeline(monkeypatch) + # water: 10 electrons (even) -> singlet -> restricted, no spin kwargs + api.dftbplus_thermo( + Atoms("OH2", positions=[[0, 0, 0.12], [0, 0.76, -0.48], [0, -0.76, -0.48]]), + directory=str(tmp_path / "j"), + ) + + assert seen["run_thermo_spin"] == 0.0 + assert "Hamiltonian_SpinPolarisation" not in seen["geoopt_kwargs"] + + +def test_dftbplus_thermo_explicit_triplet(monkeypatch, tmp_path): + api, seen = _mock_pipeline(monkeypatch) + # O2 is even-electron but a triplet -> user declares spin=1 -> 2 unpaired + api.dftbplus_thermo( + Atoms("O2", positions=[[0, 0, 0], [0, 0, 1.2]]), + directory=str(tmp_path / "j"), + spin=1.0, + ) + + assert seen["run_thermo_spin"] == 1.0 + assert seen["geoopt_kwargs"]["Hamiltonian_SpinPolarisation_UnpairedElectrons"] == 2 if __name__ == "__main__": From 6a9f7faaf628150c03276c66339a4cb26cc49b38 Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Thu, 2 Jul 2026 15:12:16 +0200 Subject: [PATCH 2/2] Cover the full 3ob element set from the authoritative spin constants The 3ob-3-1 set ships its own spin constants in extras/spinw.hsd (PBE, slateratom). All eight organic values already matched it exactly; extend the table with the remaining 3ob elements from that authoritative file (Na, Mg, K, Ca, Zn, Br, I) so spin-polarised calculations work for any element the tool's 3ob parameters already support. Verified a Br-containing radical runs spin-polarised and converges. --- ThermoScreening/calculator/dftbplus.py | 22 +++++++++++++++------- tests/calculator/test_dftbplus.py | 3 ++- 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/ThermoScreening/calculator/dftbplus.py b/ThermoScreening/calculator/dftbplus.py index 97665ac..d428ca5 100644 --- a/ThermoScreening/calculator/dftbplus.py +++ b/ThermoScreening/calculator/dftbplus.py @@ -51,22 +51,30 @@ def _slako_dir(slako_dir=None): return selected_dir + os.sep -# Atomic spin constants (Hartree) for the 3ob-3-1 PBE Slater-Koster set, i.e. the -# single spin constant of the highest occupied shell per element (Wss for H, Wpp -# otherwise), used with ShellResolvedSpin = No to match this tool's atom-resolved -# 3ob SCC. These are parameters that must match the SK set's functional (PBE) and -# parametrisation; using a different SK set requires its own spin constants. -# Cross-verified: 3ob (Elstner list) vs QUASINANO2015 (arXiv:1605.01360) vs the -# DFTB+ manual/mio (Koehler). Validate energies with a real spin-polarised run. +# Atomic spin constants (Hartree) for the 3ob-3-1 PBE Slater-Koster set: the spin +# constant of the highest occupied shell per element (Wss for H, Wpp for the +# p-block, valence Wss for the s-block and Zn), used with ShellResolvedSpin = No +# to match this tool's atom-resolved 3ob SCC. Taken from the authoritative +# ``spinw.hsd`` shipped with the 3ob-3-1 set (PBE, slateratom) and match it +# exactly. These are parameters tied to the SK set + functional; a different SK +# set (e.g. mio) needs its own constants. Verified end-to-end: an OH radical runs +# spin-polarised (0.40 eV below restricted, S_elec = R ln 2). SPIN_CONSTANTS = { "H": "{ -0.07174 }", "C": "{ -0.02265 }", "N": "{ -0.02545 }", "O": "{ -0.02785 }", "F": "{ -0.02990 }", + "Na": "{ -0.01528 }", + "Mg": "{ -0.01667 }", "P": "{ -0.01490 }", "S": "{ -0.01549 }", "Cl": "{ -0.01606 }", + "K": "{ -0.01075 }", + "Ca": "{ -0.01196 }", + "Zn": "{ -0.01680 }", + "Br": "{ -0.01377 }", + "I": "{ -0.01144 }", } diff --git a/tests/calculator/test_dftbplus.py b/tests/calculator/test_dftbplus.py index 778756f..abed4af 100644 --- a/tests/calculator/test_dftbplus.py +++ b/tests/calculator/test_dftbplus.py @@ -409,9 +409,10 @@ def test_spin_kwargs_builds_colinear_block(): def test_spin_kwargs_rejects_element_without_constant(): + # Fe is not in the 3ob spin-constant set with pytest.raises(ValueError, match="not available for element"): dftbplus_module._spin_kwargs( - Atoms("Br2", positions=[[0, 0, 0], [0, 0, 2.3]]), 0.5 + Atoms("Fe2", positions=[[0, 0, 0], [0, 0, 2.3]]), 0.5 )