Skip to content

Rank-reduced Woodbury kernel for the lattice Green's function - #16

Open
harrisonlabollita wants to merge 8 commits into
TRIQS:unstablefrom
harrisonlabollita:pr/woodbury-kernel
Open

Rank-reduced Woodbury kernel for the lattice Green's function#16
harrisonlabollita wants to merge 8 commits into
TRIQS:unstablefrom
harrisonlabollita:pr/woodbury-kernel

Conversation

@harrisonlabollita

Copy link
Copy Markdown
Collaborator

Incremental diff:
pr/obe-printing-h5...pr/woodbury-kernel

Note that this PR shares postprocess.{cpp,hpp} with PR #15 for unrelated
reasons. So we should merge #15 and rebase before merging this PR.

Summary

Every k-summed quantity in ModEST — gloc, the density, the chemical-potential
search, the charge-density correction — inverted the full
(N_\nu \times N_\nu) band-basis matrix
(G_k(\omega) = [\omega + i\delta + \mu - H(k) - P^\dagger \Sigma(\omega) P]^{-1})
at every ((k, \sigma, \omega)). This inversion is unnecessary, because for the fixed k-grid cases (to which
this port applies), (H(k)) is diagonal in the band basis and (\Sigma) is structurally space in the C-space.

This PR applies the Woodbury identity to collapse the inversion. We build up the basic backend primitives that allows us to port all lattice Green's function computations onto one shared Woodbury kernel.

Changes (organized by commit)

02957d3 — extract lattice_gf_helpers.hpp (pure refactor, no numerics change
except as noted). Pulls the direct (N_\nu \times N_\nu) inversion helpers out of
gloc_fixed_grid.hpp into a new shared header so density, chargedensity and
the post-processing routines can use them:

  • upfold_self_energy_at_freq, local_gf_at_k, lattice_gf_at_k move over, losing
    their Mesh template parameter.
  • upfold_self_energy_all_freq moves out of postprocess.cpp.
  • new make_sigma_total, folding Sigma_static into Sigma_dynamic once per call
    rather than at every frequency.
  • removes dead code: the old spectral_function sketch in postprocess.hpp, the
    impurity_levels stub in downfolding.cpp.

One deliberate numerical change here: density_for_matrix_valued_impl now reuses
lattice_gf_at_k and computes trace(G_band) - trace(inv(w + mu - H_k)).

ca59384 — the kernel, and gloc ported. Adds active_subspace_t,
detect_active_subspace, compute_bare_projected, compute_sigma_active,
make_woodbury_setup, apply_K, build_G_band to lattice_gf_helpers.hpp;
gloc_fixed_grid.hpp uses them.

fb552e2 — density and chemical-potential search. trace_G_B_m_G_KS is
replaced by the rank-reduced correction (\mathrm{tr}(K \cdot G^0_{QQ,\mathrm{sq}})),
accumulated per ((k,\sigma,\omega)). The Kohn-Sham term is unchanged, and the
(\Sigma = 0) case now short-circuits to it. Removes detail::G0_C_k_sigma
(downfolding.{hpp,cpp}) and detail::calc_inv_G_G0 (density.hpp).

1b76bf1charge_density_correction. Builds (G_k) from the identity
instead of inverting at every frequency. Two changes here are independent of the
Woodbury port
and worth reviewing separately:

  • adds an explicit matrix_valued (H(k)) fallback branch. The rank reduction
    assumes diagonal (H(k)); previously every case went through lattice_gf_at_k,
    so the matrix-valued path now keeps the direct inversion.
  • parallelises the ((k,\sigma)) loop with OpenMP.

Testing

This is a behavior-preserving change. No new tests
are added, because the existing reference-data suite already covers the kernel
densely and passes unchanged (including the LFS set under
-DEnable_LFS_Tests=ON).

Review notes

  • Read lattice_gf_helpers.hpp first; the block comment at the
    Rank-reduced Woodbury banner states the identity the rest of the header
    implements.
  • G0_QPdag is not dagger(G0_PQ) — the (i\delta) makes (D) complex. Same
    reason Lmat · KR is not Rmat† · … in build_G_band. These are the two easy
    places to introduce a wrong conjugate.
  • LaTeX notes of the implementation are also available on the ModEST paper draft.

Rename the theta projector to "partial projector" throughout the loaders and
post-processing entry points, matching the terminology used in the docs.

- add h5_read / h5_write for spectral_function_kw and spectral_function_w
- add printing for the post-process structs
Add rotate_local_basis to the obe c2py spec so the local-basis rotation used by
the post-processing routines is reachable from Python.
Pull the direct N_v x N_v inversion helpers out of gloc_fixed_grid.hpp into a
new header so density, chargedensity and the post-processing routines can share
them.

- move upfold_self_energy_at_freq, local_gf_at_k and lattice_gf_at_k, dropping
  the Mesh template parameter.
- move upfold_self_energy_all_freq out of postprocess.cpp into the same header
- add make_sigma_total, folding Sigma_static into Sigma_dynamic once per call.
- repoint chargedensity.hpp and density.hpp at the new header.
- remove dead code (the old spectral_function sketch in
  postprocess.hpp, the impurity_levels stub in downfolding.cpp)

density_for_matrix_valued_impl now reuses lattice_gf_at_k and
computes trace(G_band) - trace(inv(w + mu - H_k)) where it
previously computed trace(inv(...) - inv(...)).
For some cases, the rank of Sigma can be less than C, so we use the active
rank instead of M. The equations are:

    M0(w, v) = (w + i*delta + mu - eps_v(k))     [diagonal, H_k diagonal]
    Q        = P_active[active C rows, :]        (rank x N_v)
    K        = (Sigma_active^-1 - Q M0^-1 Q+)^-1 (rank x rank)

    diag(P G_k P+)_m = diag(P M0^-1 P+)_m + ((P M0^-1 Q+) K (Q M0^-1 P+))_mm

We solve (I - Y_aa Sigma_a) X = R and return Sigma_a X.
Sigma_a is singular at large iw_n and for a static-only
self-energy. Guards rank == 0.
Replace trace_G_B_m_G_KS with the rank-reduced correction term
tr(K . G0_QQ_sq), accumulated per (k, sigma, omega). The Kohn-Sham term is
unchanged. When Sigma iszero the density reduces to the Kohn-Sham result,
which is now short-circuited.

Remove detail::G0_C_k_sigma (downfolding.cpp/hpp) and
detail::calc_inv_G_G0 (density.hpp).
Build the band-basis lattice Green's function via

    G_k(v, v') = D(v) delta_vv' + D(v) (Q+ K Q)(v, v') D(v')

instead of inverting the full N_v x N_v matrix at every frequency.

Two changes here are independent of the Woodbury port:

- adds a matrix_valued H(k) fallback branch. Previously every case went through
  lattice_gf_at_k; the rank reduction assumes diagonal H(k), so the matrix-valued
  path now keeps the direct inversion explicitly.
- parallelises the (k, sigma) loop with OpenMP.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant