From 3b51daf92875d56abe7d29e16581fbcb49201e0d Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Tue, 7 Jul 2026 11:00:54 +0200 Subject: [PATCH 1/2] Add optional D3(BJ) dispersion to the DFTB+ engine Add an opt-in Grimme D3 dispersion correction (Becke-Johnson damping) for the DFTB+ Hamiltonian, using the 3ob-recommended parameters from the DFTB+ manual (a1=0.5719, a2=3.6017, s6=1.0, s8=0.5883). Default is no dispersion so existing results are unchanged. - _dispersion_kwargs helper in calculator/dftbplus.py - dispersion argument on dftbplus_thermo (merged into the Geoopt/Hessian kwargs) - dispersion argument on screen and a --dispersion d3-bj CLI flag Closes #88 --- ThermoScreening/calculator/dftbplus.py | 47 ++++++++++++++++++++++++++ ThermoScreening/cli/thermo.py | 6 ++++ ThermoScreening/thermo/api.py | 10 ++++-- ThermoScreening/thermo/screening.py | 5 +++ docs/usage.rst | 4 +-- tests/calculator/test_dftbplus.py | 25 ++++++++++++++ tests/thermo/test_main.py | 23 +++++++++++++ tests/thermo/test_screening.py | 24 +++++++++++-- 8 files changed, 138 insertions(+), 6 deletions(-) diff --git a/ThermoScreening/calculator/dftbplus.py b/ThermoScreening/calculator/dftbplus.py index f28a896..57a2083 100644 --- a/ThermoScreening/calculator/dftbplus.py +++ b/ThermoScreening/calculator/dftbplus.py @@ -9,6 +9,7 @@ from ase.io import read from ..utils.physicalConstants import PhysicalConstants +from ..exceptions import TSValueError # --------------------------------------------------------------------------- # @@ -197,6 +198,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 TSValueError( + 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). diff --git a/ThermoScreening/cli/thermo.py b/ThermoScreening/cli/thermo.py index 144dadc..18aa2ca 100644 --- a/ThermoScreening/cli/thermo.py +++ b/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 '.", ) + 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, diff --git a/ThermoScreening/thermo/api.py b/ThermoScreening/thermo/api.py index 39b7645..ab06ef5 100644 --- a/ThermoScreening/thermo/api.py +++ b/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 diff --git a/ThermoScreening/thermo/screening.py b/ThermoScreening/thermo/screening.py index f949180..96e958f 100644 --- a/ThermoScreening/thermo/screening.py +++ b/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, ) diff --git a/docs/usage.rst b/docs/usage.rst index 6405f9c..28a4def 100644 --- a/docs/usage.rst +++ b/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 diff --git a/tests/calculator/test_dftbplus.py b/tests/calculator/test_dftbplus.py index 72833bc..41335fe 100644 --- a/tests/calculator/test_dftbplus.py +++ b/tests/calculator/test_dftbplus.py @@ -9,6 +9,7 @@ import ThermoScreening.calculator.dftbplus as dftbplus_module from ThermoScreening.calculator import Geoopt, Hessian, Modes from ThermoScreening.calculator.dftbplus import _slako_dir, dftb_3ob_parameters +from ThermoScreening.exceptions import TSValueError from ThermoScreening.cli.dftb_setup import gbsa_param_path from ThermoScreening.thermo.api import dftbplus_thermo @@ -395,6 +396,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(TSValueError, 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) == {} diff --git a/tests/thermo/test_main.py b/tests/thermo/test_main.py index bdfda5a..d63412f 100644 --- a/tests/thermo/test_main.py +++ b/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") diff --git a/tests/thermo/test_screening.py b/tests/thermo/test_screening.py index 94057ea..4f033f5 100644 --- a/tests/thermo/test_screening.py +++ b/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 From f748c9b111808b5e3d8a25f42ada73eca84b5b33 Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Tue, 7 Jul 2026 11:10:07 +0200 Subject: [PATCH 2/2] Match module ValueError convention and test dispersion injection Raise ValueError (not TSValueError) for an unknown dispersion model, to match _spin_kwargs/resolve_parameter_set in the same module. Add tests asserting dispersion is injected into both the Geoopt and Hessian steps and absent by default. --- ThermoScreening/calculator/dftbplus.py | 3 +-- tests/calculator/test_dftbplus.py | 3 +-- tests/thermo/test_api.py | 24 ++++++++++++++++++++++++ 3 files changed, 26 insertions(+), 4 deletions(-) diff --git a/ThermoScreening/calculator/dftbplus.py b/ThermoScreening/calculator/dftbplus.py index 57a2083..85fd9fd 100644 --- a/ThermoScreening/calculator/dftbplus.py +++ b/ThermoScreening/calculator/dftbplus.py @@ -9,7 +9,6 @@ from ase.io import read from ..utils.physicalConstants import PhysicalConstants -from ..exceptions import TSValueError # --------------------------------------------------------------------------- # @@ -229,7 +228,7 @@ def _dispersion_kwargs(dispersion=None): return {} if dispersion.lower() != "d3-bj": - raise TSValueError( + raise ValueError( f"Unknown dispersion model {dispersion!r}; supported: 'd3-bj'." ) diff --git a/tests/calculator/test_dftbplus.py b/tests/calculator/test_dftbplus.py index 41335fe..1a43b39 100644 --- a/tests/calculator/test_dftbplus.py +++ b/tests/calculator/test_dftbplus.py @@ -9,7 +9,6 @@ import ThermoScreening.calculator.dftbplus as dftbplus_module from ThermoScreening.calculator import Geoopt, Hessian, Modes from ThermoScreening.calculator.dftbplus import _slako_dir, dftb_3ob_parameters -from ThermoScreening.exceptions import TSValueError from ThermoScreening.cli.dftb_setup import gbsa_param_path from ThermoScreening.thermo.api import dftbplus_thermo @@ -416,7 +415,7 @@ def test_dispersion_kwargs_is_case_insensitive(): def test_dispersion_kwargs_rejects_unknown_model(): - with pytest.raises(TSValueError, match="Unknown dispersion model"): + with pytest.raises(ValueError, match="Unknown dispersion model"): dftbplus_module._dispersion_kwargs("d4") diff --git a/tests/thermo/test_api.py b/tests/thermo/test_api.py index bd4f6fd..e9203aa 100644 --- a/tests/thermo/test_api.py +++ b/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(