Add dense AO spectral localizers and quadratic pseudospectra (#6150)

Co-authored-by: Thomas D. Kuehne <tkuehne@cp2k.org>

Author: DCM-Uni-Paderborn

docs/methods/properties/index.md

@@ -12,4 +12,6 @@ raman
 nmr
 stm_images
 inversion_topology
+spectral_localizer
+quadratic_pseudospectrum
 ```

docs/methods/properties/quadratic_pseudospectrum.md

@@ -0,0 +1,62 @@
+# Quadratic pseudospectrum
+
+`FORCE_EVAL / PROPERTIES / QUADRATIC_PSEUDOSPECTRUM` finds states approximately localized in both
+position and energy from a converged GPW/GAPW calculation. It reports no topological index and
+requires no additional library.
+
+For Hermitian covariant AO observables $A_j$, queries $\lambda_j$, positive weights $w_j$ and
+overlap $S$, the generalized eigenproblem is
+
+$$
+Q c = q S c,\qquad
+Q=\sum_j w_j^2(A_j-\lambda_j S)S^{-1}(A_j-\lambda_j S).
+$$
+
+The lowest eigenpairs describe joint approximate states. The output includes $\sqrt{q}$, separate
+energy/position residuals and expectation values. States satisfy $c^\dagger S c=1$. Residuals
+concern projected finite-basis operators, not exact continuum variances.
+
+```text
+&PROPERTIES
+  &QUADRATIC_PSEUDOSPECTRUM
+    FORMULATION FINITE
+    POSITION [bohr] 0 0 0
+    ENERGY [hartree] 0.0
+    KAPPA [hartree*bohr^-1] 0.01
+    NSTATES 2
+    &STATES
+    &END
+  &END
+&END
+```
+
+## Observables
+
+- `FINITE` uses analytic Cartesian moments. It requires matching `CELL/POISSON PERIODIC NONE`
+  without explicit k-points.
+- `PERIODIC` uses sine/cosine coordinates on the full `MP_GRID` torus at the frozen SCF potential.
+  Only periodic directions are localized. Nonperiodic mesh dimensions must be one.
+- `TRANSLATION` uses analytic translated-Gaussian overlaps for $U(a)\psi(r)=\psi(r+a)$. Specify
+  repeated `TRANSLATION_VECTOR` and `WAVE_VECTOR` queries and positive `TRANSLATION_SCALE` instead
+  of `KAPPA`. The target phase is $\exp(i k\cdot a)$ without an additional $2\pi$. This path
+  includes the full continuum translation norm and reports leakage outside the AO span, rather than
+  treating compressed AO translation matrices as unitary.
+
+The periodic paths construct a complete analysis torus, independently of SCF symmetry reduction.
+`SOC T` uses post-SCF GTH spin-orbit coupling after a restricted calculation with suitable
+pseudopotentials.
+
+## Output and limits
+
+`SOLVER DENSE` uses a metric-aware dense reference calculation. `MAX_AO` bounds the complete scalar
+AO space, including periodic cells. `EPS_METRIC` rejects ill-conditioned overlap matrices, and
+`EPS_EIGEN` bounds eigenpair residuals. `NSTATES` selects the number of states.
+
+The `STATES` print section writes complex AO coefficients. Spin-up precedes spin-down for SOC, with
+cell-major ordering inside each spin block on a torus. Individual vectors in a degenerate subspace
+are gauge dependent. Compare subspaces, residuals and observables instead of individual
+coefficients.
+
+Choose energy, position and scales for the physical question and check basis and size dependence.
+Broader examples are available in
+[TopologicalCP2K](https://github.com/DCM-Uni-Paderborn/TopologicalCP2K).

docs/methods/properties/spectral_localizer.md

@@ -0,0 +1,70 @@
+# Spectral localizers
+
+`FORCE_EVAL / PROPERTIES / SPECTRAL_LOCALIZER` evaluates a local topological index and protection
+gap from the converged GPW/GAPW Hamiltonian and AO overlap. It is independent of Wannier90 and Kubo
+transport. No additional library is required.
+
+## Finite systems
+
+For `CELL PERIODIC NONE` and `POISSON PERIODIC NONE`, the operators are the analytic Cartesian AO
+moments. For the ordered plane $(X,Y)$, the covariant localizer is
+
+$$
+L=\begin{pmatrix}
+H-ES & \kappa[(X-xS)-i(Y-yS)]\\
+\kappa[(X-xS)+i(Y-yS)] & -(H-ES)
+\end{pmatrix},\qquad B=\operatorname{diag}(S,S).
+$$
+
+Here $H$, $S$, $X$ and $Y$ are AO matrices, $(x,y)$ and $E$ are the query position and energy, and
+$\kappa$ has units of energy per length. For class A, the index is half the signature of $L$. The
+gap is the smallest absolute generalized eigenvalue of $(L,B)$, not a pivot magnitude. LAPACK
+pivoted LDL factorization is checked against the eigenvalue result.
+
+```text
+&PROPERTIES
+  &SPECTRAL_LOCALIZER
+    FORMULATION FINITE
+    INVARIANT CHERN
+    PLANE XY
+    POSITION [bohr] 0 0 0
+    ENERGY [hartree] 0.0
+    KAPPA [hartree*bohr^-1] 0.01
+  &END
+&END
+```
+
+Choose the energy and scale for the system. Repeated positions and lists of energies and scales
+produce scans. A real scalar Hamiltonian is not expected to have a nonzero class-A Chern index.
+
+## Spin-orbit coupling and periodic systems
+
+`SOC T` adds GTH spin-orbit coupling after a restricted SCF calculation in the complete AO spinor
+space. It needs SOC pseudopotential parameters and is not self-consistent noncollinear DFT.
+`INVARIANT Z2` additionally requires `TIME_REVERSAL T`. The code checks this symmetry and uses an
+orientation-sensitive Pfaffian sign in two dimensions. `DIMENSION 3` uses the class-AII chiral-block
+determinant construction and all three position directions.
+
+`FORMULATION PERIODIC` uses analytic sine/cosine AO integrals on a complete Born-von-Karman torus.
+Set `MP_GRID` to its replication factors and use `ETA` (energy units) instead of `KAPPA`.
+Nonperiodic directions must have size one. For two dimensions, `PLANE` must agree with the two
+periodic directions in both `CELL` and `POISSON`. Three dimensions require `PERIODIC XYZ`. The torus
+is built at the frozen SCF potential, independently of the SCF mesh and its symmetry reduction.
+Finite Cartesian coordinates are not substituted for periodic position operators.
+
+## Numerical controls
+
+`SOLVER DENSE` is the reference implementation. `MAX_AO` bounds the complete scalar AO space
+including torus cells. Spinors double that space, and dense localizers require quadratic storage.
+`EPS_METRIC`, `MATRIX_TOLERANCE` and `EPS_GAP` control metric conditioning, operator checks and gap
+resolution. Unresolved gaps do not yield an index. `SPECTRAL_FLATTENING` optionally replaces the
+shifted Hamiltonian by its metric-covariant matrix sign after checking an electronic gap. Its scale
+and gap are distinct from the final localizer gap.
+
+`DFT / PRINT / AO_MATRICES / POSITION` exports the finite AO moments in the same ordering as overlap
+and Hamiltonian matrices. `SOC` exports the three real antisymmetric GTH components, whose physical
+orbital operators include a factor of $i$. Both exports are restricted to isolated systems.
+
+These are finite-basis local diagnostics. Check stability against basis, system size, query position
+and scale before interpreting a material index. Broader examples are available in
+[TopologicalCP2K](https://github.com/DCM-Uni-Paderborn/TopologicalCP2K).

src/CMakeLists.txt

@@ -678,6 +678,16 @@ list(
   qs_ks_utils.F
   qs_ks_reference.F
   qs_kubo_transport.F
+  qs_spectral_localizer.F
+  qs_localizer_torus.F
+  spectral_localizer.F
+  quadratic_pseudospectrum.F
+  qs_finite_ao.F
+  qs_ao_translations.F
+  qs_quadratic_pseudospectrum.F
+  spectral_localizer_dbcsr.F
+  spectral_localizer_flatten.F
+  spectral_localizer_bloch.F
   qs_linres_atom_current.F
   qs_linres_current.F
   qs_linres_current_utils.F
@@ -1602,6 +1612,9 @@ set(CP2K_PROGS_F
     kpoint_lattice_fft_unittest.F
     kpoint_smearing_unittest.F
     topology_wilson_unittest.F
+    spectral_localizer_unittest.F
+    quadratic_pseudospectrum_unittest.F
+    pw/pw_unittest.F
     topology_symmetry_unittest.F
     topology_curvature_unittest.F
     topology_snapshot_unittest.F
@@ -1615,7 +1628,6 @@ set(CP2K_PROGS_F
     qs_fod_unittest.F
     qs_ot_complex_ref_unittest.F
     subsys/cell_types_unittest.F
-    pw/realspace_grid_cube_unittest.F
     metadyn_tools/graph.F
     motion/dumpdcd.F
     motion/xyz2dcd.F
@@ -2037,6 +2049,9 @@ list(
   "kpoint_lattice_fft_unittest"
   "kpoint_smearing_unittest"
   "topology_wilson_unittest"
+  "spectral_localizer_unittest"
+  "quadratic_pseudospectrum_unittest"
+  "pw_unittest"
   "topology_symmetry_unittest"
   "topology_curvature_unittest"
   "topology_snapshot_unittest"
@@ -2050,7 +2065,6 @@ list(
   "qs_fod_unittest"
   "qs_ot_complex_ref_unittest"
   "cell_types_unittest"
-  "realspace_grid_cube_unittest"
   "graph"
   "dumpdcd"
   "xyz2dcd"
@@ -2081,6 +2095,10 @@ add_executable(gapw_1c_basis_unittest gapw_1c_basis_unittest.F)
 add_executable(kpoint_lattice_fft_unittest kpoint_lattice_fft_unittest.F)
 add_executable(kpoint_smearing_unittest kpoint_smearing_unittest.F)
 add_executable(topology_wilson_unittest topology_wilson_unittest.F)
+add_executable(spectral_localizer_unittest spectral_localizer_unittest.F)
+add_executable(quadratic_pseudospectrum_unittest
+               quadratic_pseudospectrum_unittest.F)
+add_executable(pw_unittest pw/pw_unittest.F)
 add_executable(topology_symmetry_unittest topology_symmetry_unittest.F)
 add_executable(topology_curvature_unittest topology_curvature_unittest.F)
 add_executable(topology_snapshot_unittest topology_snapshot_unittest.F)
@@ -2097,7 +2115,6 @@ add_executable(qs_charge_mixing_unittest qs_charge_mixing_unittest.F)
 add_executable(qs_fod_unittest qs_fod_unittest.F)
 add_executable(qs_ot_complex_ref_unittest qs_ot_complex_ref_unittest.F)
 add_executable(cell_types_unittest subsys/cell_types_unittest.F)
-add_executable(realspace_grid_cube_unittest pw/realspace_grid_cube_unittest.F)
 add_executable(graph metadyn_tools/graph.F)
 add_executable(dumpdcd motion/dumpdcd.F)
 add_executable(xyz2dcd motion/xyz2dcd.F)

src/common/bibliography.F

@@ -98,9 +98,9 @@ MODULE bibliography
                             KuhneHeskeProdan2020, Schreder2021, Schreder2024_1, Schreder2024_2, &
                             Shiga2022, Lindh1995, Chai2024a, Rullan2026, Sundararaman2017, Andreussi2019, &
                             Chai2025a, Neugebauer2023, Chai2024b, Moynihan2017, Cococcioni2005, &
-                            Cances2000, Herbst2022, Zhang2025, Gresch2017, Soluyanov2011, Davidson1975, &
-                            Olsen1988, Stratmann1998, Vecharynski2017, Stathopoulos1998, Bai2012, &
-                            Fukaya2014, Schreder2024_3
+                            Cances2000, Herbst2022, Zhang2025, Gresch2017, Soluyanov2011, Dixon2023, Doll2025, &
+                            Wong2026, Loring2015, Cerjan2023, Bairnsfather2026, Davidson1975, &
+                            Olsen1988, Stratmann1998, Vecharynski2017, Stathopoulos1998, Bai2012, Fukaya2014, Schreder2024_3
 
 CONTAINS
 
@@ -2317,6 +2317,41 @@ SUBROUTINE add_all_references()
                          source="Phys. Rev. B", volume="83", pages="235401", &
                          year=2011, doi="10.1103/PhysRevB.83.235401")
 
+      CALL add_reference(key=Dixon2023, &
+                         authors=s2a("K. Y. Dixon", "T. A. Loring", "A. Cerjan"), &
+                         title="Classifying topology in photonic heterostructures with gapless environments", &
+                         source="Phys. Rev. Lett.", volume="131", pages="213801", &
+                         year=2023, doi="10.1103/PhysRevLett.131.213801")
+
+      CALL add_reference(key=Wong2026, &
+                         authors=s2a("S. Wong", "I. Yamazaki", "C. Siefert", "I. Duff", "T. A. Loring", "A. Cerjan"), &
+                         title="Efficient local classification of parity-based material topology", &
+                         source="Phys. Rev. B", volume="114", pages="084201", &
+                         year=2026, doi="10.1103/6hj9-tgct")
+
+      CALL add_reference(key=Loring2015, &
+                         authors=s2a("T. A. Loring"), &
+                         title="K-theory and pseudospectra for topological insulators", &
+                         source="Ann. Phys.", volume="356", pages="383-416", &
+                         year=2015, doi="10.1016/j.aop.2015.02.031")
+
+      CALL add_reference(key=Doll2025, &
+                         authors=s2a("N. Doll", "T. A. Loring", "H. Schulz-Baldes"), &
+                         title="Topological indices for periodic gapped Hamiltonians and fuzzy tori", &
+                         source="Math. Phys. Anal. Geom.", volume="28", pages="13", &
+                         year=2025, doi="10.1007/s11040-025-09508-0")
+
+      CALL add_reference(key=Cerjan2023, &
+                         authors=s2a("A. Cerjan", "T. A. Loring", "F. Vides"), &
+                         title="Quadratic pseudospectrum for identifying localized states", &
+                         source="J. Math. Phys.", volume="64", pages="023501", &
+                         year=2023, doi="10.1063/5.0098336")
+
+      CALL add_reference(key=Bairnsfather2026, &
+                         authors=s2a("C. A. Bairnsfather", "R. M. Kaufmann", "T. A. Loring", "A. Cerjan"), &
+                         title="Band Unfolding via the Quadratic Pseudospectrum", &
+                         source="arXiv", pages="2605.05423", year=2026, doi="10.48550/arXiv.2605.05423")
+
       CALL add_reference(key=Davidson1975, &
                          authors=s2a("E. R. Davidson"), &
                          title="The iterative calculation of a few of the lowest eigenvalues and corresponding "// &

src/cp_dbcsr_output.F

@@ -579,10 +579,10 @@ SUBROUTINE write_matrix_gen(matrix, matrix_name, before, after, para_env, &
             DO irow = first_row, last_row
                IF (omit_headers) THEN
                   WRITE (UNIT=output_unit, FMT=fmtstr2) &
-                     irow, (matrix(irow, jcol), jcol=from, to)
+                     (matrix(irow, jcol), jcol=from, to)
                ELSE
                   WRITE (UNIT=output_unit, FMT=fmtstr2) &
-                     (matrix(irow, jcol), jcol=from, to)
+                     irow, (matrix(irow, jcol), jcol=from, to)
                END IF
             END DO
          END DO

src/input_cp2k_print_dft.F

@@ -486,6 +486,23 @@ SUBROUTINE create_print_dft_section(section)
       CALL section_add_keyword(print_key, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword, __LOCATION__, name="POSITION", &
+                          description="Print the analytic AO matrices <mu|x,y,z|nu> in Bohr, relative to "// &
+                          "the coordinate origin, in the same AO ordering as OVERLAP and KOHN_SHAM_MATRIX. "// &
+                          "Only isolated GPW/GAPW systems; not a periodic Berry-phase operator. "// &
+                          "Use NDIGITS 16 for numerical postprocessing.", &
+                          default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
+      CALL section_add_keyword(print_key, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, __LOCATION__, name="SOC", &
+                          description="Print the three real antisymmetric AO GTH SOC components Vx, Vy, Vz in hartree. "// &
+                          "The physical orbital operators are i*Vx, i*Vy, i*Vz. "// &
+                          "Requires a restricted, isolated GPW/GAPW calculation with GTH SOC potentials.", &
+                          default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
+      CALL section_add_keyword(print_key, keyword)
+      CALL keyword_release(keyword)
+
       CALL keyword_create(keyword, __LOCATION__, name="FERMI_CONTACT", &
                           description="If the printkey is activated controls the printing of the Fermi contact matrix", &
                           default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)