Add GBSA/ALPB implicit solvation to the DFTB+ pipeline - #66
Merged
Conversation
Add a `solvent=` knob so geometry optimisation, energy, and the Hessian all
run in implicit solvent. DFTB+ includes the GBSA/SASA term in the second
derivatives, so the optimised geometry and the frequencies stay consistent
with the solvated energy.
- calculator/dftbplus.py: `_solvation_kwargs` builds the
`Solvation = GeneralizedBorn { ParamFile = ... }` block (empty for the
gas-phase default), passing an absolute parameter-file path since DFTB+
resolves it relative to the run directory.
- cli/dftb_setup.py: download GBSA parameter files (grimme-lab/gbsa-parameters)
by solvent name; `setup-dftb --solvent <name>` fetches them.
- thermo/api.py, thermo/screening.py, cli/thermo.py: thread `solvent`
(and an explicit `solvation_param_file` override) through `dftbplus_thermo`,
`screen`, and the `screen`/`setup-dftb` CLI.
- thermo/screening.py: report `Eelec_hartree` (the electronic energy, which
carries the solvation term) and `G_total_hartree` (electronic + Gibbs
correction), so the solvation free energy is visible in the results.
The published GBSA parameters are fit for GFN-xTB; used with 3ob/mio they are
an approximation, so an explicit method-consistent parameter file can be
supplied instead. Validated against real DFTB+: solvated water is stabilised
by ~14 kcal/mol and the effect appears in the reported free energy.
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #66 +/- ##
==========================================
+ Coverage 96.31% 96.43% +0.12%
==========================================
Files 22 22
Lines 1358 1404 +46
==========================================
+ Hits 1308 1354 +46
Misses 50 50
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
galjos
added a commit
that referenced
this pull request
Jul 5, 2026
Adds a `quasi_rrho` option (CLI `--quasi-rrho`) for the vibrational entropy. ## Why The rigid-rotor-harmonic-oscillator entropy of a mode **diverges as its frequency → 0**, so low-frequency modes (floppy torsions, weak/hindered rotors, non-covalent complexes) get a spuriously large entropy — a well-known RRHO weakness. This is also the dominant error left in implicit-solvent free energies (#66). ## What Grimme's quasi-RRHO ([Chem. Eur. J. 2012, 18, 9955](https://doi.org/10.1002/chem.201200497)): each mode's entropy is interpolated between the harmonic-oscillator and free-rotor values, weighted by `w = 1 / (1 + (100/ν)^4)`. High-frequency modes stay harmonic; low-frequency modes approach the finite free-rotor limit. - **`thermo/thermo.py`** — `Thermo(..., quasi_rrho=False)`; `_compute_vibrational_entropy` uses the new `_quasi_rrho_entropy` interpolation when enabled. The harmonic path is numerically unchanged, so the **default preserves existing results**. - **`thermo/api.py`, `thermo/screening.py`, `cli/thermo.py`** — thread `quasi_rrho` through `run_thermo`, `dftbplus_thermo`, `screen`, and `screen --quasi-rrho`. Engine-independent (it acts on the frequencies), so it also applies to a future xTB engine. ## Usage ``` thermo screen mols/ --solvent water --quasi-rrho ``` ## Validation - **Physics**: a molecule with 25/40 cm⁻¹ modes has its entropy reduced by ~3.4 cal/mol/K (≈ +1 kcal/mol in G) — the spurious low-mode entropy is tamed; a high-frequency-only case is unchanged (ΔS ≈ 0.01). - **Real DFTB+**: ethane (lowest real vibration 277 cm⁻¹, all > 100) gives qRRHO ≈ harmonic, and the harmonic S = 54.5 cal/mol/K matches the experimental ~54.8 — confirming the correct mode set is used. - Offline suite green (206 passed); real-DFTB+ integration green (54 passed, 0 skips); codecov patch fully covered.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Adds a
solvent=knob so a screen can run in implicit solvent. The geometry optimisation, energy, and Hessian all run in solution — DFTB+ folds the GBSA/SASA term into the second derivatives, so the optimised geometry and the frequencies stay consistent with the solvated energy (verified: gas vs. solvated Hessians differ).What changed
calculator/dftbplus.py—_solvation_kwargsbuilds theSolvation = GeneralizedBorn { ParamFile = ... }block (empty for the gas-phase default). The parameter-file path is made absolute because DFTB+ resolves it relative to the run directory (jobs run in per-molecule dirs).cli/dftb_setup.py— download GBSA parameter files (grimme-lab/gbsa-parameters) by solvent name;thermo setup-dftb --solvent <name>fetches them. 14 solvents (water, acetonitrile, dmso, thf, chcl3, …).thermo/api.py,thermo/screening.py,cli/thermo.py— threadsolvent(and an explicitsolvation_param_fileoverride) throughdftbplus_thermo,screen, and thescreen/setup-dftbCLI.thermo/screening.py— reportEelec_hartree(the electronic energy, which carries the solvation term) andG_total_hartree(electronic + Gibbs correction). Without this the solvation free energy would be invisible, since the existing columns are thermal corrections only.Usage
Scientific note
The published GBSA/ALPB parameters are fit for GFN-xTB; used with the DFTB (3ob/mio) Hamiltonians they are an approximation. Pass
solvation_param_file=to use a method-consistent set instead. (On a future xTB engine this caveat disappears —--alpbuses the natively-fit parameters.)Validation
screen(parameter_set="3ob", solvent="water"): solvated water is stabilised by ~14 kcal/mol (ΔG_total = -13.89 kcal/mol), now surfaced in the CSV.DFTB_PREFIX(33 passed, 0 skips); codecov patch fully covered.