cp2k/tests/SE
Dynamics of Condensed Matter cccd2f3934
Fix TB gradients, virials, and k-point symmetry
# 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.
2026-05-06 00:47:46 +02:00
..
regtest Add tests for force constant 2025-09-12 09:43:02 +02:00
regtest-2-1 Adjust tolerance of SE/regtest-2-1/hcl_ot.inp 2025-07-20 10:20:55 +02:00
regtest-2-2 Fix TB gradients, virials, and k-point symmetry 2026-05-06 00:47:46 +02:00
regtest-3-1 Tests: Convert TEST_TYPES to matchers.py 2025-03-01 19:07:38 +01:00
regtest-3-2 Tests: Convert TEST_TYPES to matchers.py 2025-03-01 19:07:38 +01:00
regtest-3-3 CG/Linesearch (#4141) 2025-04-15 20:13:30 +02:00
regtest-3-4 CG/Linesearch (#4141) 2025-04-15 20:13:30 +02:00
regtest-4 Add support for functions with multiple arguments to fparser 2026-02-04 10:40:02 +01:00
sample_cif Unset executable flag [ci skip] 2019-11-20 14:47:56 +01:00
sample_xtl support for XTL file format (MSI) 2009-05-19 21:10:32 +00:00
sample_xyz Unset executable flag [ci skip] 2019-11-20 14:47:56 +01:00
untested_inputs Format all input files 2023-12-18 22:41:19 +01:00