Run open-shell species spin-polarised in DFTB+ automatically

ThermoScreening/calculator/dftbplus.py

@@ -51,6 +51,80 @@ def _slako_dir(slako_dir=None):
     return selected_dir + os.sep
 
 
+# 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 }",
+}
+
+
+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).

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()

tests/calculator/test_dftbplus.py

@@ -390,6 +390,32 @@ 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():
+    # Fe is not in the 3ob spin-constant set
+    with pytest.raises(ValueError, match="not available for element"):
+        dftbplus_module._spin_kwargs(
+            Atoms("Fe2", 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.",

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__":