Add mio as a selectable DFTB parameter set

ThermoScreening/calculator/dftbplus.py

@@ -51,15 +51,16 @@ 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 = {
+# Atomic spin constants (Hartree): 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 the atom-resolved SCC. These
+# are parameters tied to the Slater-Koster set + functional, so each parameter set
+# has its own.
+#
+# 3ob-3-1 (PBE): taken from the authoritative ``spinw.hsd`` shipped with the set
+# (calculated with PBE/slateratom) and match it exactly. Verified end-to-end: an
+# OH radical runs spin-polarised (0.40 eV below restricted, S_elec = R ln 2).
+SPIN_CONSTANTS_3OB = {
     "H": "{ -0.07174 }",
     "C": "{ -0.02265 }",
     "N": "{ -0.02545 }",
@@ -77,27 +78,41 @@ def _slako_dir(slako_dir=None):
     "I": "{ -0.01144 }",
 }
 
+# mio-1-1 reuses the 3ob spin constants. mio is the LDA-based Elstner-1998 set
+# (PRB 58, 7260) and ships no spin constants of its own; the only well-documented
+# organic spin constants (3ob's ``spinw.hsd`` and the DFTB+ manual's own H2O
+# example, H = -0.072 / O Wpp = -0.028) are PBE values. The atomic spin constant
+# is only weakly functional-dependent for H/C/N/O/S, so the authoritative 3ob
+# values are the best available choice for mio too. Verified end-to-end: a mio OH
+# radical runs spin-polarised and is stabilised relative to the restricted run.
+SPIN_CONSTANTS_MIO = SPIN_CONSTANTS_3OB
 
-def _spin_kwargs(atoms, spin):
+# Default (3ob) spin constants.
+SPIN_CONSTANTS = SPIN_CONSTANTS_3OB
+
+
+def _spin_kwargs(atoms, spin, spin_constants=SPIN_CONSTANTS_3OB):
     """
     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.
+    injects the ``spin_constants`` for the elements present.
 
     Parameters
     ----------
     atoms : ase.Atoms
         The atoms whose elements need spin constants.
     spin : float or None
         Spin quantum number S.
+    spin_constants : dict
+        Element -> spin-constant brace string, matching the Slater-Koster set.
 
     Raises
     ------
     ValueError
-        If an element has no tabulated 3ob spin constant.
+        If an element has no tabulated spin constant.
     """
     if spin is None or float(spin) <= 0.0:
         return {}
@@ -107,11 +122,12 @@ def _spin_kwargs(atoms, spin):
         return {}
 
     elements = sorted(set(atoms.get_chemical_symbols()))
-    missing = [element for element in elements if element not in SPIN_CONSTANTS]
+    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."
+            f"{', '.join(missing)}: no spin constant is tabulated for this "
+            "parameter set."
         )
 
     kwargs = {
@@ -121,7 +137,7 @@ def _spin_kwargs(atoms, spin):
         "Hamiltonian_SpinConstants_ShellResolvedSpin": "No",
     }
     for element in elements:
-        kwargs[f"Hamiltonian_SpinConstants_{element}"] = SPIN_CONSTANTS[element]
+        kwargs[f"Hamiltonian_SpinConstants_{element}"] = spin_constants[element]
     return kwargs
 
 
@@ -477,3 +493,75 @@ def read(self):
     # Parser options
     ParserOptions_ParserVersion=12,
 )
+
+
+# mio-1-1 (DFTB2): the original mio set is a second-order model, so it has no
+# ThirdOrderFull / Hubbard derivatives. Same SCC + Fermi machinery as 3ob;
+# spin-polarised runs use SPIN_CONSTANTS_MIO. Verified end-to-end against real
+# DFTB+ on an OH radical.
+dftb_mio_parameters = dict(
+    # SCC
+    Hamiltonian_SCC="Yes",
+    Hamiltonian_MaxSCCIterations=250,
+    Hamiltonian_SCCTolerance='1.0e-7',
+    Hamiltonian_ReadInitialCharges="No",
+
+    # Fermi smearing
+    Hamiltonian_Filling="Fermi {",
+    Hamiltonian_Filling_empty="Temperature [Kelvin] = 300",
+
+    # Convergence helper
+    Hamiltonian_Mixer="DIIS{}",
+
+    # Are guessed by ase
+    Hamiltonian_MaxAngularMomentum_="",
+
+    # Analysis
+    Analysis_="",
+    Analysis_CalculateForces="Yes",
+    Analysis_MullikenAnalysis="Yes",
+
+    # Parser options
+    ParserOptions_ParserVersion=12,
+)
+
+
+# Selectable DFTB parameter sets: name -> (Hamiltonian parameters, spin constants).
+# ``dftbplus_thermo`` / ``screen`` pick a set by name so a run needs only a
+# structure + charge (+ optional spin); everything else follows from the set.
+DFTB_PARAMETER_SETS = {
+    "3ob": (dftb_3ob_parameters, SPIN_CONSTANTS_3OB),
+    "3ob-3-1": (dftb_3ob_parameters, SPIN_CONSTANTS_3OB),
+    "mio": (dftb_mio_parameters, SPIN_CONSTANTS_MIO),
+    "mio-1-1": (dftb_mio_parameters, SPIN_CONSTANTS_MIO),
+}
+
+
+def resolve_parameter_set(parameter_set):
+    """
+    Resolve a parameter-set name to ``(hamiltonian_parameters, spin_constants)``.
+
+    Parameters
+    ----------
+    parameter_set : str
+        One of the keys of :data:`DFTB_PARAMETER_SETS` (e.g. ``"3ob"``, ``"mio"``).
+
+    Returns
+    -------
+    tuple(dict, dict)
+        A fresh copy of the Hamiltonian parameters and the matching spin
+        constants for the set.
+
+    Raises
+    ------
+    ValueError
+        If ``parameter_set`` is not a known set.
+    """
+    try:
+        parameters, spin_constants = DFTB_PARAMETER_SETS[parameter_set]
+    except KeyError:
+        known = ", ".join(sorted(DFTB_PARAMETER_SETS))
+        raise ValueError(
+            f"Unknown DFTB parameter set {parameter_set!r}; choose one of: {known}."
+        )
+    return dict(parameters), spin_constants

ThermoScreening/cli/dftb_setup.py

@@ -14,12 +14,52 @@
 
 
 DEFAULT_PARAMETER_SET = "3ob-3-1"
-DEFAULT_SLAKO_URL = (
-    "https://github.com/dftbparams/3ob/releases/latest/download/"
-    f"{DEFAULT_PARAMETER_SET}.tar.xz"
-)
 REQUIRED_PARAMETER_FILE = "C-C.skf"
 
+# Download URLs for the supported Slater-Koster sets (dftbparams GitHub releases).
+PARAMETER_SET_URLS = {
+    "3ob-3-1": (
+        "https://github.com/dftbparams/3ob/releases/latest/download/3ob-3-1.tar.xz"
+    ),
+    "mio-1-1": (
+        "https://github.com/dftbparams/mio/releases/latest/download/mio-1-1.tar.xz"
+    ),
+}
+
+# Short names accepted by the CLI (mapped to the canonical archive/directory name).
+PARAMETER_SET_ALIASES = {"3ob": "3ob-3-1", "mio": "mio-1-1"}
+
+
+def _canonical_set_name(parameter_set: str) -> str:
+    """
+    Map a short set name (e.g. ``"mio"``) to its canonical name (``"mio-1-1"``).
+    """
+
+    return PARAMETER_SET_ALIASES.get(parameter_set, parameter_set)
+
+
+def slako_url(parameter_set: str = DEFAULT_PARAMETER_SET) -> str:
+    """
+    Return the download URL for a Slater-Koster parameter set.
+
+    Raises
+    ------
+    ValueError
+        If ``parameter_set`` is not a known set.
+    """
+
+    name = _canonical_set_name(parameter_set)
+    try:
+        return PARAMETER_SET_URLS[name]
+    except KeyError:
+        known = ", ".join(sorted(PARAMETER_SET_URLS))
+        raise ValueError(
+            f"Unknown parameter set {parameter_set!r}; choose one of: {known}."
+        )
+
+
+DEFAULT_SLAKO_URL = PARAMETER_SET_URLS[DEFAULT_PARAMETER_SET]
+
 
 @dataclass(frozen=True)
 class Diagnostic:
@@ -40,13 +80,16 @@ def default_install_root() -> Path:
     return Path.home() / ".local" / "share" / "thermoscreening" / "slakos"
 
 
-def default_parameter_dir(install_root: str | Path | None = None) -> Path:
+def default_parameter_dir(
+    install_root: str | Path | None = None,
+    parameter_set: str = DEFAULT_PARAMETER_SET,
+) -> Path:
     """
-    Return the default 3ob parameter directory under an install root.
+    Return the parameter directory for ``parameter_set`` under an install root.
     """
 
     root = Path(install_root).expanduser() if install_root is not None else default_install_root()
-    return root / DEFAULT_PARAMETER_SET
+    return root / _canonical_set_name(parameter_set)
 
 
 def _download_file(url: str, destination: Path, timeout: int = 60) -> None:
@@ -77,15 +120,32 @@ def _safe_extract_tar(archive_path: Path, destination: Path) -> None:
 
 def install_slakos(
     install_root: str | Path | None = None,
-    url: str = DEFAULT_SLAKO_URL,
+    url: str | None = None,
     force: bool = False,
+    parameter_set: str = DEFAULT_PARAMETER_SET,
 ) -> Path:
     """
-    Download and extract the default Slater-Koster parameter set.
+    Download and extract a Slater-Koster parameter set.
+
+    Parameters
+    ----------
+    install_root : str or Path, optional
+        Directory where parameter sets are installed. Defaults to the user-local
+        share directory.
+    url : str, optional
+        Archive URL. Defaults to the release URL for ``parameter_set``.
+    force : bool
+        Re-download even when the set already exists.
+    parameter_set : str
+        Set to install (``"3ob"``/``"3ob-3-1"`` or ``"mio"``/``"mio-1-1"``).
     """
 
+    name = _canonical_set_name(parameter_set)
+    if url is None:
+        url = slako_url(name)
+
     root = Path(install_root).expanduser() if install_root is not None else default_install_root()
-    parameter_dir = default_parameter_dir(root)
+    parameter_dir = default_parameter_dir(root, name)
     marker_file = parameter_dir / REQUIRED_PARAMETER_FILE
 
     if marker_file.exists() and not force:
@@ -94,7 +154,7 @@ def install_slakos(
     root.mkdir(parents=True, exist_ok=True)
 
     with tempfile.TemporaryDirectory(prefix="thermoscreening-dftb-") as tmp_dir:
-        archive_path = Path(tmp_dir) / f"{DEFAULT_PARAMETER_SET}.tar.xz"
+        archive_path = Path(tmp_dir) / f"{name}.tar.xz"
         _download_file(url, archive_path)
         _safe_extract_tar(archive_path, root)
 

ThermoScreening/cli/thermo.py

@@ -5,7 +5,6 @@
 import time
 
 from ThermoScreening.cli.dftb_setup import (
-    DEFAULT_SLAKO_URL,
     check_dftb_setup,
     dftb_prefix_export,
     format_diagnostics,
@@ -34,17 +33,23 @@ def _command_parser():
 
     setup_parser = subparsers.add_parser(
         "setup-dftb",
-        help="Download the default DFTB+ Slater-Koster parameter set.",
+        help="Download a DFTB+ Slater-Koster parameter set.",
     )
     setup_parser.add_argument(
         "--install-root",
         default=None,
         help="Directory where Slater-Koster parameter sets are installed.",
     )
+    setup_parser.add_argument(
+        "--parameter-set",
+        default="3ob",
+        choices=["3ob", "mio"],
+        help="Parameter set to download (default '3ob').",
+    )
     setup_parser.add_argument(
         "--url",
-        default=DEFAULT_SLAKO_URL,
-        help="Archive URL for the default Slater-Koster parameter set.",
+        default=None,
+        help="Archive URL override (defaults to the release URL for the set).",
     )
     setup_parser.add_argument(
         "--force",
@@ -84,6 +89,11 @@ def _command_parser():
         "--directory", default="screening",
         help="Root working directory; each molecule runs in <directory>/<name>.",
     )
+    screen_parser.add_argument(
+        "--parameter-set", default="3ob", choices=["3ob", "mio"],
+        help="Slater-Koster parameter set (default '3ob'). Selects the Hamiltonian "
+        "parameters and matching spin constants.",
+    )
 
     return parser
 
@@ -120,6 +130,7 @@ def run_setup_dftb(parser_args):
         install_root=parser_args.install_root,
         url=parser_args.url,
         force=parser_args.force,
+        parameter_set=parser_args.parameter_set,
     )
 
     print("Slater-Koster files: ", parameter_dir)
@@ -152,6 +163,7 @@ def run_screen(parser_args):
         temperature=parser_args.temperature,
         pressure=parser_args.pressure,
         directory=parser_args.directory,
+        parameter_set=parser_args.parameter_set,
     )
 
     failed = sum(1 for record in results if record["status"] != "ok")

ThermoScreening/thermo/api.py

@@ -18,7 +18,7 @@
 from .thermo import Thermo
 from .atoms import Atom
 from ..calculator import Geoopt, Hessian, Modes
-from ..calculator.dftbplus import _spin_kwargs
+from ..calculator.dftbplus import _spin_kwargs, SPIN_CONSTANTS_3OB
 
 
 logger = logging.getLogger(__package_name__).getChild("api")
@@ -504,6 +504,7 @@ def dftbplus_thermo(
         charge=0.0,
         directory=None,
         spin=None,
+        spin_constants=None,
         **kwargs
     ):
     """
@@ -529,6 +530,10 @@ def dftbplus_thermo(
         (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.
+    spin_constants : dict, optional
+        Element -> spin-constant mapping matching the Slater-Koster set in use.
+        Defaults to the 3ob constants; pass ``SPIN_CONSTANTS_MIO`` (or the value
+        from :func:`resolve_parameter_set`) when running the mio set.
 
     Other Parameters
     ----------------
@@ -547,7 +552,10 @@ def dftbplus_thermo(
         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)
+    if spin_constants is None:
+        spin_constants = SPIN_CONSTANTS_3OB
+
+    spin_kwargs = _spin_kwargs(atoms, spin, spin_constants)
 
     with _run_in_directory(directory):
         # run geometry optimization