# 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.
The new keyword IGNORE_CONVERGENCE_FAILURE can be set to obtain the old permissive behaviour.
It replaces the keyword NOTCONV_STOPALL which enforced until now the strict behaviour which is now the default.
* Fix ks_csr_write and s_csr_write not implemented for k-points
Printing the Kohn-Sham and overlap matrices is not implemented for
k-point calculations. Both matrices are now printed at the
k grid used in the scf procedure.
* Added csr_write to TB
Added calls to ks_scr_write and s_csr_write to
qs_scf_post_tb.F. Before these print sections were silently
ignored.
* Move ks_csr_write and s_csr_write to new module
The ks_csr_write and s_csr_write move to the new module
qs_post_scf_csr_write. Also moved the bulk of the printing
function to a helper function to avoid copying of code.
* Add ks_csr_write and s_csr_write to regtests
Add the ks_csr_write and s_csr_write sections to three regtests from
DFTB, xTB and QS. These cover gamma-point, k-point with and without
real wave fucntions.
This PR enables ELPA by default when incorporated at compile-time:
1. Enables ELPA by default (if compiled into the binary),
2. Selects ELPA-kernel according to CPUID.
The code uses m_cpuid (machine module) and introduced get_elpa_kernel_index (cp_fm_elpa module) to select the ELPA-kernel according to the CPUID.
The GLOBAL entry `PREFERRED_DIAG_LIBRARY SL` is given for a number of tests that would otherwise exceed tolerances as currently used (mostly CDFT tests). Slightly adjusted some tolerances of non-CDFT tests when reasonable, and one for the PSMP case.
Other (unrelated) changes: fixed some typos (source code comments).
