Add optional D3(BJ) dispersion to the DFTB+ engine

ThermoScreening/calculator/dftbplus.py

@@ -197,6 +197,52 @@ def _solvation_kwargs(solvent=None, param_file=None, install_root=None):
     }
 
 
+# Grimme D3 with Becke-Johnson damping, using the parameters the DFTB+ manual
+# recommends for the DFTB3/3ob Hamiltonian (Brandenburg et al., J. Chem. Phys.
+# 143, 054110 (2015)).
+_D3_BJ_3OB = {"a1": 0.5719, "a2": 3.6017, "s6": 1.0, "s8": 0.5883}
+
+
+def _dispersion_kwargs(dispersion=None):
+    """
+    ASE ``Dftb`` keyword arguments adding a Grimme dispersion correction.
+
+    Parameters
+    ----------
+    dispersion : str, optional
+        Dispersion model. ``"d3-bj"`` adds Grimme D3 with Becke-Johnson damping
+        using the 3ob-recommended parameters. ``None`` (default) adds nothing,
+        so existing gas-phase results are unchanged.
+
+    Returns
+    -------
+    dict
+        The ``Hamiltonian_Dispersion*`` kwargs (empty when ``dispersion`` is None).
+
+    Raises
+    ------
+    TSValueError
+        If ``dispersion`` is not a supported model.
+    """
+    if dispersion is None:
+        return {}
+
+    if dispersion.lower() != "d3-bj":
+        raise ValueError(
+            f"Unknown dispersion model {dispersion!r}; supported: 'd3-bj'."
+        )
+
+    p = _D3_BJ_3OB
+    return {
+        "Hamiltonian_Dispersion": "DftD3 {",
+        "Hamiltonian_Dispersion_Damping": "BeckeJohnson {",
+        "Hamiltonian_Dispersion_Damping_a1": p["a1"],
+        "Hamiltonian_Dispersion_Damping_a2": p["a2"],
+        "Hamiltonian_Dispersion_s6": p["s6"],
+        "Hamiltonian_Dispersion_s8": p["s8"],
+    }
+
+
 class Geoopt(Dftb):
     """
     Custom DFTB+ calculator to optimize the system with the 'GeometryOptimisation' driver (Rational).

ThermoScreening/cli/thermo.py

@@ -107,6 +107,11 @@ def _command_parser():
         help="GBSA/ALPB implicit-solvation solvent (e.g. 'water') applied to every "
         "molecule. Default gas phase. Install with 'setup-dftb --solvent <name>'.",
     )
+    screen_parser.add_argument(
+        "--dispersion", default=None, choices=["d3-bj"],
+        help="Add a Grimme dispersion correction to the DFTB+ Hamiltonian "
+        "('d3-bj', 3ob-recommended parameters). Default none.",
+    )
     screen_parser.add_argument(
         "--quasi-rrho", action="store_true",
         help="Use Grimme's quasi-RRHO vibrational entropy (better for low-frequency "
@@ -233,6 +238,7 @@ def run_screen(parser_args):
         directory=parser_args.directory,
         parameter_set=parser_args.parameter_set,
         solvent=parser_args.solvent,
+        dispersion=parser_args.dispersion,
         quasi_rrho=parser_args.quasi_rrho,
         engine=parser_args.engine,
         method=parser_args.method,

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, _solvation_kwargs, SPIN_CONSTANTS_3OB
+from ..calculator.dftbplus import _spin_kwargs, _solvation_kwargs, _dispersion_kwargs, SPIN_CONSTANTS_3OB
 from ..calculator.xtb import optimise_and_frequencies, xtb_calculator
 from ..calculator.xtb_cli import run_xtb
 
@@ -517,6 +517,7 @@ def dftbplus_thermo(
         spin_constants=None,
         solvent=None,
         solvation_param_file=None,
+        dispersion=None,
         quasi_rrho=False,
         **kwargs
     ):
@@ -555,6 +556,10 @@ def dftbplus_thermo(
     solvation_param_file : str, optional
         Explicit path to a GBSA parameter file, overriding ``solvent`` (use a
         method-consistent set instead of the default GFN-fit one).
+    dispersion : str, optional
+        Dispersion correction to add to the Hamiltonian. ``"d3-bj"`` uses Grimme
+        D3 with Becke-Johnson damping (3ob-recommended parameters). Defaults to
+        no dispersion. Requires a DFTB+ build with D3 support.
     quasi_rrho : bool
         If True, use Grimme's quasi-RRHO treatment for the vibrational entropy,
         which tames the entropy of low-frequency modes. Default False.
@@ -595,7 +600,8 @@ def dftbplus_thermo(
     # Implicit solvation (empty for the gas-phase default). Applied to both the
     # optimisation and the Hessian so the geometry and frequencies are consistent.
     solvation_kwargs = _solvation_kwargs(solvent, solvation_param_file)
-    engine_kwargs = {**spin_kwargs, **solvation_kwargs, **kwargs}
+    dispersion_kwargs = _dispersion_kwargs(dispersion)
+    engine_kwargs = {**spin_kwargs, **solvation_kwargs, **dispersion_kwargs, **kwargs}
 
     with _run_in_directory(directory):
         # run geometry optimization

ThermoScreening/thermo/screening.py

@@ -180,6 +180,7 @@ def screen(
     spin=None,
     parameter_set="3ob",
     solvent=None,
+    dispersion=None,
     quasi_rrho=False,
     engine="dftb+",
     method="GFN2-xTB",
@@ -219,6 +220,9 @@ def screen(
         DFTB+ engine the GBSA parameters are GFN-xTB-fit, so solvation free
         energies are only qualitative (and can be non-monotonic in the dielectric
         for small neutral solutes).
+    dispersion : str, optional
+        Dispersion correction for the DFTB+ engine. ``"d3-bj"`` adds Grimme
+        D3(BJ) with the 3ob-recommended parameters. Defaults to none.
     quasi_rrho : bool
         If True, use Grimme's quasi-RRHO treatment for the vibrational entropy
         (recommended for flexible molecules with low-frequency modes). Default
@@ -306,6 +310,7 @@ def screen(
                     spin=job.spin,
                     spin_constants=spin_constants,
                     solvent=solvent,
+                    dispersion=dispersion,
                     quasi_rrho=quasi_rrho,
                     **parameters,
                 )

docs/usage.rst

@@ -29,8 +29,8 @@ The ``thermo`` command has four subcommands.
 
     .. code-block:: bash
 
-       # DFTB+ (3ob) with water solvation and quasi-RRHO entropy
-       thermo screen molecules/ --parameter-set 3ob --solvent water --quasi-rrho
+       # DFTB+ (3ob) with water solvation, D3(BJ) dispersion and quasi-RRHO entropy
+       thermo screen molecules/ --parameter-set 3ob --solvent water --dispersion d3-bj --quasi-rrho
 
        # native xtb, radical anions in solution
        thermo screen molecules/ --engine xtb-cli --solvent water

tests/calculator/test_dftbplus.py

@@ -395,6 +395,30 @@ def get_global_number_of_atoms(self):
     assert hessian.read().shape == (9, 9)
 
 
+def test_dispersion_kwargs_none_is_empty():
+    assert dftbplus_module._dispersion_kwargs() == {}
+    assert dftbplus_module._dispersion_kwargs(None) == {}
+
+
+def test_dispersion_kwargs_d3_bj_parameters():
+    kw = dftbplus_module._dispersion_kwargs("d3-bj")
+    assert kw["Hamiltonian_Dispersion"] == "DftD3 {"
+    assert kw["Hamiltonian_Dispersion_Damping"] == "BeckeJohnson {"
+    assert kw["Hamiltonian_Dispersion_Damping_a1"] == 0.5719
+    assert kw["Hamiltonian_Dispersion_Damping_a2"] == 3.6017
+    assert kw["Hamiltonian_Dispersion_s6"] == 1.0
+    assert kw["Hamiltonian_Dispersion_s8"] == 0.5883
+
+
+def test_dispersion_kwargs_is_case_insensitive():
+    assert dftbplus_module._dispersion_kwargs("D3-BJ") == dftbplus_module._dispersion_kwargs("d3-bj")
+
+
+def test_dispersion_kwargs_rejects_unknown_model():
+    with pytest.raises(ValueError, match="Unknown dispersion model"):
+        dftbplus_module._dispersion_kwargs("d4")
+
+
 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) == {}

tests/thermo/test_api.py

@@ -489,6 +489,30 @@ def test_dftbplus_thermo_injects_solvation_into_both_steps(monkeypatch, tmp_path
     assert seen["hessian_kwargs"]["Hamiltonian_Solvation"] == "GeneralizedBorn {"
 
 
+def test_dftbplus_thermo_injects_dispersion_into_both_steps(monkeypatch, tmp_path):
+    api, seen = _mock_pipeline(monkeypatch)
+
+    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"),
+        dispersion="d3-bj",
+    )
+
+    for step in ("geoopt_kwargs", "hessian_kwargs"):
+        assert seen[step]["Hamiltonian_Dispersion"] == "DftD3 {"
+        assert seen[step]["Hamiltonian_Dispersion_Damping_a1"] == 0.5719
+
+
+def test_dftbplus_thermo_no_dispersion_by_default(monkeypatch, tmp_path):
+    api, seen = _mock_pipeline(monkeypatch)
+    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 "Hamiltonian_Dispersion" not in seen["geoopt_kwargs"]
+
+
 def test_dftbplus_thermo_gas_phase_has_no_solvation(monkeypatch, tmp_path):
     api, seen = _mock_pipeline(monkeypatch)
     api.dftbplus_thermo(

tests/thermo/test_main.py

@@ -168,6 +168,29 @@ def test_main_runs_doctor(monkeypatch, capsys):
     assert capsys.readouterr().out == "not ready\n"
 
 
+def test_run_screen_forwards_dispersion(monkeypatch):
+    captured = {}
+
+    def fake_screen(source, **kwargs):
+        captured.update(kwargs)
+        return []
+
+    monkeypatch.setattr(thermo, "screen", fake_screen)
+    args = argparse.Namespace(
+        source="mols", out="results", charge=0.0, temperature=298.15,
+        pressure=101325, directory="screening", parameter_set="3ob",
+        solvent=None, dispersion="d3-bj", quasi_rrho=False, engine="dftb+",
+        method="GFN2-xTB", resume=False,
+    )
+    assert thermo.run_screen(args) == 0
+    assert captured["dispersion"] == "d3-bj"
+
+
+def test_parse_args_screen_dispersion_choice():
+    assert thermo.parse_args(["screen", "mols"]).dispersion is None
+    assert thermo.parse_args(["screen", "mols", "--dispersion", "d3-bj"]).dispersion == "d3-bj"
+
+
 def test_main_runs_conformers(monkeypatch):
     monkeypatch.setattr(
         thermo, "parse_args", lambda: argparse.Namespace(command="conformers")

tests/thermo/test_screening.py

@@ -146,6 +146,26 @@ def fake_thermo(atoms, solvent=None, **kwargs):
     assert captured["solvent"] == "water"
 
 
+def test_screen_passes_dispersion_to_dftbplus_thermo(monkeypatch, tmp_path):
+    _write_xyz(tmp_path / "mol.xyz")
+
+    captured = {}
+
+    def fake_thermo(atoms, dispersion=None, **kwargs):
+        captured["dispersion"] = dispersion
+        return _FakeThermo()
+
+    monkeypatch.setattr(screening, "dftbplus_thermo", fake_thermo)
+    screening.screen(
+        str(tmp_path),
+        out=str(tmp_path / "r"),
+        directory=str(tmp_path / "runs"),
+        dispersion="d3-bj",
+    )
+
+    assert captured["dispersion"] == "d3-bj"
+
+
 def test_screen_dispatches_to_xtb_engine(monkeypatch, tmp_path):
     _write_xyz(tmp_path / "mol.xyz")
 
@@ -423,8 +443,8 @@ 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, engine="dftb+", method="GFN2-xTB",
-        resume=False,
+        solvent=None, dispersion=None, quasi_rrho=False, engine="dftb+",
+        method="GFN2-xTB", resume=False,
     )
 
     assert cli.run_screen(args) == 1  # one molecule failed