# Fix TB gradients, virials, tblite reference checks, and k-point symmetry
## Summary
This PR continues the CP2K/tblite interface work originally started by Johann Pototschnig
(`jpoto`) and relies in part on Sebastian Ehlert's `tblite` library. It makes the CP2K tight
binding paths substantially more consistent for nuclear gradients, virials/stress tensors, and
k-point symmetry, especially for DFTB, CP2K-native GFN0/GFN1, and `&XTB &TBLITE` GFN1/GFN2/IPEA1.
Main changes:
- fix several CP2K/tblite force and virial terms for periodic systems,
- add a native `tblite` CLI reference path under `&XTB &TBLITE &REFERENCE_CLI`,
- enable and test k-point symmetry reduction for CP2K-native TB and CP2K/tblite TB paths,
- improve DFTB SCC/Ewald virials and update affected DFTB references,
- add/extend regtests for molecular and truly periodic systems, including force/stress DEBUG checks,
- update the CP2K toolchain patches for `tblite`/embedded DFTD4 and standalone DFTD4.
## User-visible input changes
- `&DFT &XTB` can now select the library implementation via
`GFN_TYPE TBLITE`; in that case an `&XTB &TBLITE` section is required.
- `&XTB &TBLITE` supports `METHOD GFN1`, `METHOD GFN2`, and `METHOD IPEA1`.
- `&XTB &TBLITE` has an optional `REFERENCE_CLI` keyword and a lone
`&REFERENCE_CLI` subsection. If present/enabled, CP2K writes a `.gen` geometry, runs native
`tblite run --grad`, and compares energy, gradient, and virial.
- `SCC_MIXER TBLITE` is documented as a modified Broyden mixer; the associated controls are named
`TBLITE_MIXER_DAMPING` and `TBLITE_MIXER_MAX_SCF`.
- `DO_EWALD` remains available as a deprecated/debug override; the manual now documents the required
periodic `POISSON` setup for CP2K-native DFTB/xTB and the `TBLITE` exception.
- CP2K warns if `CELL%PERIODIC` and `POISSON%PERIODIC` differ, but still runs.
## What was fixed
### CP2K/tblite GFN1/GFN2/IPEA1
- Periodic image-cell matrix construction now fills all image cells instead of only image cell 1.
- H0, overlap/potential, GFN2 multipole, CN-chain-rule, and D4-related virial paths were checked
term by term against finite differences.
- Periodic self-image virial signs/factors were corrected in the CP2K/tblite interface.
- The previously bad strongly distorted Si8/GFN2 Gamma stress case was traced to a hard D4 two-body
cutoff discontinuity and fixed by adding a smooth switching region in the D4 path.
- `RUN_TYPE DEBUG` force/stress checks now behave consistently for molecular and periodic
CP2K/tblite GFN1/GFN2 cases; isolated stress is allowed for tblite with a warning because it is a
diagnostic virial, not a physical isolated-system stress.
### Native tblite CLI reference
The `REFERENCE_CLI` path initially showed good molecular agreement but poor periodic virials for
Ar-fcc:
- GFN1 periodic virial diff was about `1.15e-2`,
- GFN2 periodic virial diff was about `2.95e-2`.
The parser, units, row order, and sign convention were checked and were correct. The remaining error
was in native `tblite`: periodic H0 self-image pairs (`iat == jat`, non-zero translation) used the
normal pair sigma factor `0.5`. The correct symmetrized cell-derivative factor for these self-image
pairs is `0.25`.
Updated `REFERENCE_CLI` comparison after the native tblite H0 sigma fix:
| Case | SCF | Energy diff | Gradient diff sum/max | Virial diff sum/max |
| --- | ---: | ---: | ---: | ---: |
| CH2O/GFN1, molecular | 12 | `6.775398020409e-10` | `3.296429618975e-08` / `1.577689361792e-08` | `2.884140131919e-08` / `2.625531908551e-08` |
| CH2O/GFN2, molecular | 15 | `6.651852402229e-10` | `1.437665769393e-07` / `6.050656990025e-08` | `1.240845506761e-07` / `8.358997352631e-08` |
| Ar-fcc/GFN1, periodic | 7 | `3.595097552989e-09` | `1.125627343711e-09` / `3.688615600028e-10` | `5.441065574157e-09` / `1.653417298474e-09` |
| Ar-fcc/GFN2, periodic | 20 | `8.975906240494e-09` | `2.140900031075e-08` / `4.424293700931e-09` | `3.094086449720e-07` / `8.057020835922e-08` |
`REFERENCE_CLI` is intentionally skipped/guarded for CP2K multi-k-point runs: native `tblite` CLI is
a useful Gamma reference, but it is not a CP2K-equivalent k-point/symmetry implementation.
### DFTB and CP2K-native xTB/GFN0/GFN1
- DFTB SCC/Ewald background handling was fixed so the homogeneous background contributes to energy
and atomic properties without entering the SCC overlap virial as a constant potential shift.
- The TB real-space Ewald self-image virial prefactor was fixed.
- CP2K-native DFTB/GFN0/GFN1 k-point symmetry support was extended and tested with the CP2K atomic
backend and SPGLIB backend.
- Full atomic k-point symmetry is used for high-symmetry production geometries; for displaced DEBUG
geometries the reduced symmetry is recomputed instead of reusing an invalid high-symmetry map.
- DFTB3 forces, Gamma stress, k-point forces, and k-point stress were tested and covered by new
regtests.
## What was tried
The debugging was deliberately systematic rather than a single-code-path patch:
- molecular and truly periodic crystals were checked separately; large vacuum cells were not used as
periodic substitutes,
- GFN1 and GFN2 were tested with tight `EPS_DEFAULT`/`EPS_SCF`, DX scans, and high `MAX_SCF`,
- CP2K/tblite results were compared to native `tblite` CLI for energy, gradient, and virial,
- H0, overlap/potential, Coulomb/multipole, CN-chain-rule, D4, D3/s-D3, DFTD4 and SCC/Ewald pieces
were isolated where practical and checked against finite differences,
- CP2K-native DFTB, DFTB3, GFN0, GFN1, CP2K/tblite GFN1/GFN2/IPEA1, and native `tblite` CLI were
tested separately,
- k-point full grids, CP2K atomic symmetry, SPGLIB symmetry, displaced atoms, skew cells, and
MPI-rank reproducibility were checked,
- 1D/2D periodicity, charged systems, and open-shell/spin behavior were probed.
What helped:
- restoring periodic H0/CN self-image contributions,
- fixing self-image virial signs/factors,
- fixing the native tblite H0 self-image sigma factor,
- fixing the DFTB SCC/Ewald background gauge contribution,
- recomputing k-point symmetry for lowered-symmetry DEBUG geometries,
- smoothing the hard D4 two-body cutoff in the stress-sensitive GFN2/Si8 case,
- updating the external library toolchain patches and DFTB reference values.
What did not help or was not the root cause:
- simply tightening `EPS_SCF`, `EPS_DEFAULT`, `DX`, or `MAX_SCF`,
- changing only SCC damping/mixer parameters for the hard GFN2 stress cases,
- treating the native `tblite` CLI periodic virial as authoritative before fixing its own H0 sigma
self-image term,
- using large vacuum boxes as stand-ins for true periodic crystals,
- assuming all heteronuclear crystals are good regtests: several are dominated by overlap
conditioning/SCF-state sensitivity rather than a clean force/stress bug.
## What works now
- CP2K/tblite GFN1/GFN2 molecular force DEBUG checks are at `~1e-8`.
- CP2K/tblite GFN2 isolated diagnostic stress is `~4e-9`.
- CP2K/tblite GFN1/GFN2 Ar-fcc k-point force/stress DEBUG runs with `64 -> 4` k-point reduction.
- CP2K/tblite GFN2 `PERIODIC XZ` Ar-layer k-point stress DEBUG is reduced from the reproduced `5.83e-4` remainder to `2.66e-5`; forces are `2.29e-6`.
- CP2K/tblite GFN1/GFN2 skew/asym Si8 Gamma stress is now in the `~1e-6` range rather than
`~1e-3`.
- CP2K/tblite IPEA1 runs with k-point reduction and force/stress DEBUG coverage.
- CP2K-native DFTB non-SCC/SCC and GFN0/GFN1 k-point symmetry paths pass the new tests.
- DFTB3 force, Gamma stress, k-point force, and k-point stress tests pass.
- DFTD4/DFTD3/D3(BJ) CP2K tests pass with the updated build; external DFTD4/tblite patches are ready
for the toolchain.
Representative final regtest values:
| Path | Representative result |
| --- | ---: |
| CP2K/tblite GFN2 CH2O force DEBUG | `M042 = 1e-8` |
| CP2K/tblite GFN2 CH2O isolated stress DEBUG | `M042 = 4.067e-9` |
| CP2K/tblite GFN2 Si8 Gamma stress DEBUG | `M042 = 4.76966e-7` |
| CP2K/tblite GFN2 Ar-fcc k-point stress DEBUG | `M042 = 3.33e-7`, `N_special_kpoints = 4` |
| CP2K/tblite GFN2 Ar-layer XZ k-point stress DEBUG | `M042 = 2.66e-5`, `N_special_kpoints = 2` |
| CP2K/tblite GFN1 Ar-fcc k-point stress DEBUG | `M042 = 1.6417e-8`, `N_special_kpoints = 4` |
| CP2K/tblite IPEA1 Ar-fcc k-point stress DEBUG | `M042 = 2.4656e-8`, `N_special_kpoints = 4` |
| CP2K-native GFN0 Si k-point symmetry | `N_special_kpoints = 4`; displaced geometry `12` |
| DFTB non-SCC Si k-point stress | `M042 = 3.807e-9`, `N_special_kpoints = 4` |
| DFTB SCC C k-point stress | `M042 = 5.409e-9` |
| DFTB SCC ZnS k-point stress | `M042 = 2.18113e-7`, `N_special_kpoints = 4` |
| DFTB3 H2O k-point stress | `M042 = 6.89e-10`, `N_special_kpoints = 4` |
## Remaining limitations
- The reproduced `PERIODIC XZ` CP2K/tblite GFN2 Ar-layer k-point virial remainder is now much smaller (`2.66e-5` instead of `5.83e-4`), but not yet in the `1e-6` range; forces are good.
- `PERIODIC X` full stress DEBUG is limited by the general CP2K cell/Poisson/stress machinery, not by
a TB-specific quick fix.
- CP2K/tblite open-shell/LSD still aborts intentionally:
`LSD option not compatible with tblite library`.
- Native `tblite` CLI remains a serial external Gamma-reference path; it is not a k-point/symmetry
reference implementation.
- Some heteronuclear CP2K/tblite crystals with k-points are sensitive to overlap conditioning and
SCF state; Ar/Si remain the robust PR regtest choices.
## Verification
- formatting:
- `tools/precommit/format_fortran.py` on changed Fortran files,
- `tools/precommit/format_input_file.py` on changed/new CP2K inputs,
- `tools/precommit/precommit.py -m` on non-Fortran/non-input files.
- patch checks:
- `tblite-0.5.0-gradient-fixes.patch`: `patch -p1 --dry-run` against fresh `tblite-0.5.0`,
- `dftd4-4.1.0-gradient-fixes.patch`: `patch -p1 --dry-run` against fresh `dftd4-4.1.0`.
- build:
- `cmake --build build-tblite-mpi-check --target cp2k-bin -- -j 8`
- regtests with `mpiranks=2`, `ompthreads=1`:
- `xTB/regtest-tblite-gfn2`, `xTB/regtest-tblite-gfn1-grad`,
`xTB/regtest-tblite-ipea1`, `DFTB/regtest-scc-2`: `49 / 49` correct,
- `DFTB/regtest-scc`, `DFTB/regtest-nonscc`: after updating affected references,
`89 / 89` correct,
- additionally checked in the same PR-ready run: `xTB/regtest-1`, `xTB/regtest-3`,
`xTB/regtest-gfn0`, `QS/regtest-dft-vdw-corr-4`; those were green, with only the old DFTB
references requiring the follow-up update above.
This PR shows an attempt to implement Gaussian smearing method.
The broading of it is much narrower than Fermi-Dirac distribution, and in most cases it should give a very resonable result. However, although Gaussian smearing is useful in practice, it does not have a clear physical meaning and therefore cannot be used in situations where electronic temperature or thermal excitations are of interest.
Most of the logic is adapted from how the Fermi–Dirac method works; the main differences lie in the mathematical formulas. In addition, an analytical form of df_i/dε_j is attempted here.
