# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
# SPDX-License-Identifier: MIT
#
# DGAmore — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
# Eliashberg Equation Solver for Strongly Correlated Electron Systems
r"""
Local Schwinger-Dyson step. Given the two-particle DMFT Green's functions and the bare interaction, the functions
here build the local vertex hierarchy per spin channel — the generalized susceptibility :math:`\chi_{r}`, the
irreducible vertex :math:`\Gamma_{r}` (with the Kitatani shell asymptotics), the auxiliary susceptibility
:math:`\chi^{*}_{r}`, the three-leg vertex :math:`\gamma_{r}` (``vrg``), the full vertex :math:`F_{r}`, and the
physical susceptibility — and recompute the local self-energy via the Schwinger-Dyson equation as a sanity check
against the DMFT input. Equation numbers refer to the author's master's thesis (Chapter 3). A second set of
functions implements the alternative ab-initio DGA formulation.
"""
import numpy as np
import dgamore.config as config
from dgamore.bubble_gen import BubbleGenerator
from dgamore.greens_function import GreensFunction
from dgamore.interaction import LocalInteraction
from dgamore.local_four_point import LocalFourPoint
from dgamore.matsubara_frequencies import MFHelper
from dgamore.n_point_base import SpinChannel, FrequencyNotation
from dgamore.self_energy import SelfEnergy
[docs]
def create_generalized_chi(g2: LocalFourPoint, g_dmft: GreensFunction) -> LocalFourPoint:
r"""
Returns the generalized susceptibility, see also Eq. (3.41) in my master's thesis,
:math:`\chi_{r;abcd}^{\omega\nu\nu'} = \beta (G_{r;abcd}^{(2);\omega\nu\nu'} - 2 \delta_{r,\mathrm{dens}}
\delta_{\omega 0} G_{ab}^{\nu} G_{cd}^{\nu'})`. The disconnected term is subtracted only in the density
(ph) channel at :math:`\omega = 0`.
:param g2: Two-particle (DMFT) Green's function :math:`G^{(2)}_{r}` as a :class:`LocalFourPoint`.
:param g_dmft: The local (DMFT) :class:`GreensFunction`.
:return: The generalized susceptibility :math:`\chi_{r}` as a :class:`LocalFourPoint` (half niw range).
"""
chi = config.sys.beta * g2.to_half_niw_range()
if g2.channel == SpinChannel.DENS and g2.frequency_notation == FrequencyNotation.PH:
g_loc_slice_mat = g_dmft.mat[..., g_dmft.niv - config.box.niv_core : g_dmft.niv + config.box.niv_core]
ggv_mat = g_loc_slice_mat[:, :, None, None, :, None] * g_loc_slice_mat[None, None, :, :, None, :]
chi[:, :, :, :, 0, ...] -= 2.0 * config.sys.beta * ggv_mat
return chi
[docs]
def create_gamma_r(gchi_r: LocalFourPoint, gchi0_inv: LocalFourPoint, beta: float) -> LocalFourPoint:
r"""
Returns the ph-irreducible vertex
:math:`\Gamma_{r;abcd}^{\omega\nu\nu'} = \beta^2 [(\chi_{r;abcd}^{\omega\nu\nu'})^{-1} -
(\delta_{\nu\nu'}\chi_{0;abcd}^{\omega\nu\nu'})^{-1}]`.
:param gchi_r: The generalized susceptibility :math:`\chi_{r}`.
:param gchi0_inv: The inverse bare bubble :math:`\chi_0^{-1}` (diagonal in :math:`\nu\nu'`).
:param beta: Inverse temperature :math:`\beta`.
:return: The irreducible vertex :math:`\Gamma_{r}` as a :class:`LocalFourPoint`.
"""
return beta**2 * (gchi_r.invert() - gchi0_inv)
[docs]
def create_gamma_r_with_shell_correction(
gchi_r: LocalFourPoint, gchi0: LocalFourPoint, u_loc: LocalInteraction
) -> LocalFourPoint:
r"""
Calculates the irreducible vertex with the shell correction as described by Motoharu Kitatani
et al. 2022 J. Phys. Mater. 5 034005; DOI 10.1088/2515-7639/ac7e6d. More specifically, see equations A.4 and A.8.
The irreducible vertex has an additional factor of :math:`1/\beta^2` compared to DGApy. This is also described in
my master's thesis, Sec. 3.7.2.
:param gchi_r: The generalized susceptibility :math:`\chi_{r}` (core frequency box).
:param gchi0: The bare bubble :math:`\chi_0` over the full frequency box.
:param u_loc: The bare local interaction :math:`U`.
:return: The shell-corrected irreducible vertex :math:`\Gamma_{r}` as a :class:`LocalFourPoint`.
"""
chi_tilde_shell = (gchi0.invert() + 1.0 / config.sys.beta**2 * u_loc.as_channel(gchi_r.channel)).invert()
chi_tilde_core_inv = chi_tilde_shell.cut_niv(config.box.niv_core).invert()
return config.sys.beta**2 * (gchi_r.invert() - chi_tilde_core_inv) + u_loc.as_channel(gchi_r.channel)
[docs]
def create_auxiliary_chi(gamma_r: LocalFourPoint, gchi0_inv: LocalFourPoint, u_loc: LocalInteraction) -> LocalFourPoint:
r"""
Returns the auxiliary susceptibility, see Eq. (3.60) in my master's thesis,
:math:`\chi^{*;\omega\nu\nu'}_{r;abcd} = ((\chi_{0;abcd}^{\omega\nu})^{-1} +
(\Gamma_{r;abcd}^{\omega\nu\nu'} - U_{r;abcd})/\beta^2)^{-1}`.
:param gamma_r: The irreducible vertex :math:`\Gamma_{r}`.
:param gchi0_inv: The inverse bare bubble :math:`\chi_0^{-1}` (core box).
:param u_loc: The bare local interaction :math:`U`.
:return: The auxiliary susceptibility :math:`\chi^{*}_{r}` as a :class:`LocalFourPoint`.
"""
return (gchi0_inv + (gamma_r - u_loc.as_channel(gamma_r.channel)) / config.sys.beta**2).invert()
[docs]
def create_generalized_chi_with_shell_correction(
gchi_aux_sum: LocalFourPoint, gchi0: LocalFourPoint, u_loc: LocalInteraction
) -> LocalFourPoint:
r"""
Calculates the generalized susceptibility with the shell correction as described by
Motoharu Kitatani et al. 2022 J. Phys. Mater. 5 034005; DOI 10.1088/2515-7639/ac7e6d. Eq. A.15. This is also
described in my master's thesis, Sec. 3.7.2.
:param gchi_aux_sum: The frequency-summed auxiliary susceptibility :math:`\sum_{\nu\nu'} \chi^{*}_{r}`.
:param gchi0: The bare bubble :math:`\chi_0` over the full frequency box.
:param u_loc: The bare local interaction :math:`U`.
:return: The shell-corrected physical susceptibility :math:`\chi_{r}^{\omega}` as a :class:`LocalFourPoint`.
"""
gchi0_full_sum = 1.0 / config.sys.beta * gchi0.sum_over_all_vn(config.sys.beta)
gchi0_core_sum = 1.0 / config.sys.beta * gchi0.cut_niv(config.box.niv_core).sum_over_all_vn(config.sys.beta)
return ((gchi_aux_sum + gchi0_full_sum - gchi0_core_sum).invert() + u_loc.as_channel(gchi_aux_sum.channel)).invert()
[docs]
def create_full_vertex_from_gamma(gamma_r, gchi0, u_loc):
r"""
Returns the local full vertex in the ``niv_full`` region from the irreducible vertex,
:math:`F = \Gamma [1 + \frac{1}{\beta^2} \chi_0 \Gamma]^{-1}` (with :math:`\Gamma` padded with :math:`U` beyond the core box).
:param gamma_r: The irreducible vertex :math:`\Gamma_{r}` (core box).
:param gchi0: The bare bubble :math:`\chi_0` (with its fermionic axis taken on the diagonal).
:param u_loc: The bare local interaction :math:`U`, used to pad the shell.
:return: The full vertex :math:`F_{r}` as a :class:`LocalFourPoint`.
"""
gamma_urange = gamma_r.pad_with_u(u_loc.as_channel(gamma_r.channel), config.box.niv_full)
return gamma_urange @ (
LocalFourPoint.identity_like(gamma_urange) + 1.0 / config.sys.beta**2 * gchi0 @ gamma_urange
).invert(False)
[docs]
def create_full_vertex(gchi_r: LocalFourPoint, gchi0_inv: LocalFourPoint) -> LocalFourPoint:
r"""
Returns the local full vertex in the ``niv_core`` region, see Eq. (3.58) in my master's thesis,
:math:`F_{r;abcd}^{\omega\nu\nu'} = \beta^2 (\chi_{0;abcd}^{-1} - \chi_{0;abef}^{-1} \chi_{r;fehg}
\chi_{0;ghcd}^{-1})`.
:param gchi_r: The generalized susceptibility :math:`\chi_{r}`.
:param gchi0_inv: The inverse bare bubble :math:`\chi_0^{-1}`.
:return: The full vertex :math:`F_{r}` as a :class:`LocalFourPoint`.
"""
return config.sys.beta**2 * (gchi0_inv - gchi0_inv @ gchi_r @ gchi0_inv)
[docs]
def create_vrg(gchi_aux: LocalFourPoint, gchi0_inv: LocalFourPoint) -> LocalFourPoint:
r"""
Returns the three-leg vertex, see Eq. (3.63) in my master's thesis,
:math:`\gamma_{r;abcd}^{\omega\nu} = \beta (\chi^{\omega\nu\nu}_{0;ablm})^{-1} (\sum_{\nu'}
\chi^{*;\omega\nu\nu'}_{r;mlcd})`.
:param gchi_aux: The auxiliary susceptibility :math:`\chi^{*}_{r}`.
:param gchi0_inv: The inverse bare bubble :math:`\chi_0^{-1}` (core box).
:return: The three-leg vertex :math:`\gamma_{r}` (``vrg``) as a :class:`LocalFourPoint`.
"""
gchi_aux_sum = gchi_aux.sum_over_vn(config.sys.beta, axis=(-1,))
return config.sys.beta * (gchi0_inv @ gchi_aux_sum)
[docs]
def create_vertex_functions(
g2_r: LocalFourPoint,
gchi0: LocalFourPoint,
gchi0_inv_core: LocalFourPoint,
g_dmft: GreensFunction,
u_loc: LocalInteraction,
) -> tuple[LocalFourPoint, LocalFourPoint, LocalFourPoint, LocalFourPoint, LocalFourPoint]:
r"""
Builds the full local vertex hierarchy for a single spin channel: the irreducible vertex :math:`\Gamma_{r}`
(with shell correction), the frequency-summed (shell-corrected) physical susceptibility, the three-leg vertex
:math:`\gamma_{r}`, the full vertex :math:`F_{r}`, and the generalized susceptibility :math:`\chi_{r}`. Employs
explicit asymptotics as proposed by Motoharu Kitatani et al. 2022 J. Phys. Mater. 5 034005;
DOI 10.1088/2515-7639/ac7e6d for the local irreducible vertex.
:param g2_r: The two-particle (DMFT) Green's function :math:`G^{(2)}_{r}` for this channel.
:param gchi0: The bare bubble :math:`\chi_0` over the full frequency box.
:param gchi0_inv_core: The inverse bare bubble :math:`\chi_0^{-1}` over the core box (diagonal in :math:`\nu`).
:param g_dmft: The local (DMFT) :class:`GreensFunction`.
:param u_loc: The bare local interaction :math:`U`.
:return: The tuple ``(gamma_r, gchi_r_sum, vrg_r, f_r, gchi_r)`` of local vertex functions.
"""
logger = config.logger
gchi_r = create_generalized_chi(g2_r, g_dmft)
logger.info(f"Local generalized susceptibility chi^wvv' ({gchi_r.channel.value}) calculated.")
gamma_r = create_gamma_r_with_shell_correction(gchi_r, gchi0, u_loc)
gchi0 = gchi0.take_vn_diagonal()
logger.info(
f"Local irreducible vertex Gamma^wvv' ({gamma_r.channel.value})"
f"{" with asymptotic correction" if config.box.niv_shell > 0 else ""} calculated."
)
f_r = create_full_vertex_from_gamma(gamma_r, gchi0, u_loc)
logger.info(f"Local full vertex F^wvv' ({f_r.channel.value}) calculated.")
gchi_r_aux = create_auxiliary_chi(gamma_r, gchi0_inv_core, u_loc)
logger.info(f"Local auxiliary susceptibility chi^*wvv' ({gchi_r_aux.channel.value}) calculated.")
vrg_r = create_vrg(gchi_r_aux, gchi0_inv_core)
logger.info(f"Local three-leg vertex gamma^wv ({vrg_r.channel.value}) calculated.")
gchi_r_aux_sum = gchi_r_aux.sum_over_all_vn(config.sys.beta)
del gchi_r_aux
gchi_r_aux_sum = create_generalized_chi_with_shell_correction(gchi_r_aux_sum, gchi0, u_loc)
logger.info(
f"Updated local susceptibility chi^w ({gchi_r_aux_sum.channel.value})"
f"{" with asymptotic correction" if config.box.niv_shell > 0 else ""}."
)
return gamma_r, gchi_r_aux_sum, vrg_r, f_r, gchi_r
[docs]
def get_local_hartree_fock(u_loc: LocalInteraction, occ: np.ndarray) -> np.ndarray:
r"""
Returns the local Hartree-Fock (static, frequency-independent) self-energy
:math:`\Sigma^{HF}_{ab}` from the bare interaction and the local occupation, i.e. the density-channel
interaction contracted with the occupation, see Eq. (3.55) in my master's thesis.
The interaction tensor is stored with the inter-orbital density :math:`U'` at :math:`U_{abab}` (the convention
of :meth:`Hamiltonian.kanamori_interaction_dp` and the w2dynamics ``umatrix`` files), whereas the density-channel
projection contracted as ``"abcd,dc->ab"`` picks up :math:`U_{aabb}`. The middle two orbital indices are
therefore swapped (``"abcd->acbd"``) before the projection, so that the Hartree term uses :math:`U'` while the
Fock term still uses :math:`U_{adcb}`. This only affects multi-orbital systems with off-diagonal interactions;
single-orbital and purely orbital-diagonal interactions are unchanged.
:param u_loc: The bare local interaction :math:`U`.
:param occ: The local occupation matrix :math:`n_{ab}`, shape ``[n_bands, n_bands]``.
:return: The Hartree-Fock self-energy, shape ``[n_bands, n_bands]``.
"""
return u_loc.permute_orbitals("abcd->acbd").as_channel(SpinChannel.DENS).times("abcd,dc->ab", occ)
[docs]
def get_loc_self_energy_vrg(
vrg_dens: LocalFourPoint,
vrg_magn: LocalFourPoint,
gchi_dens_sum: LocalFourPoint,
gchi_magn_sum: LocalFourPoint,
g_dmft: GreensFunction,
u_loc: LocalInteraction,
) -> SelfEnergy:
r"""
Performs the local self-energy calculation using the Schwinger-Dyson equation, i.e. the local variant of Eq. (3.64)
in my master's thesis. This is done to verify the implementation of the Schwinger-Dyson equation with the three-leg
vertex and the local susceptibility against the DMFT self-energy. Note that there will never be a perfect match due
to the sampling method of w2dynamics and the stochastic nature of the CTQMC solver. Nevertheless, the results should
be very close. For more details, see also Paul Worm's PhD thesis, Eq. (3.70) and Anna Galler's PhD Thesis, P. 76 ff.
:param vrg_dens: The density three-leg vertex :math:`\gamma_{\mathrm{dens}}`.
:param vrg_magn: The magnetic three-leg vertex :math:`\gamma_{\mathrm{magn}}`.
:param gchi_dens_sum: The frequency-summed density susceptibility :math:`\chi_{\mathrm{dens}}^{\omega}`.
:param gchi_magn_sum: The frequency-summed magnetic susceptibility :math:`\chi_{\mathrm{magn}}^{\omega}`.
:param g_dmft: The local (DMFT) :class:`GreensFunction`.
:param u_loc: The bare local interaction :math:`U`.
:return: The local :class:`SelfEnergy` (including the Hartree-Fock term).
"""
# 1=i, 2=j, 3=k, 4=l, 7=o, 8=p
g_wv = g_dmft.get_g_wv(MFHelper.wn(config.box.niw_core), config.box.niv_core)
inner = vrg_dens - vrg_dens @ u_loc.as_channel(SpinChannel.DENS) @ gchi_dens_sum
inner -= vrg_magn - vrg_magn @ u_loc.as_channel(SpinChannel.MAGN) @ gchi_magn_sum
inner = 0.5 * inner.to_full_niw_range()
sigma_sum = -1.0 / config.sys.beta * u_loc.times("kjop,ilpowv,lkwv->ijv", inner, g_wv)
hartree_fock = get_local_hartree_fock(u_loc, config.sys.occ_dmft)[..., None]
return SelfEnergy((hartree_fock + sigma_sum)[None, None, None, ...], beta=config.sys.beta)
# ----------------------------------------------- AbinitioDGA algorithms -----------------------------------------------
#
# DEVELOPMENT / TESTING ONLY. The functions below are an alternative ("ab-initio DGA") rewriting of the local
# Schwinger-Dyson cross-check: they build the local self-energy from the full vertex F and the density three-leg
# vertex derived from it (density channel only, by design), as opposed to the auxiliary-susceptibility (Hedin)
# form used by the production routine above. They are NOT part of the production routine (``DGAmore.py`` calls
# :func:`perform_local_schwinger_dyson`) and exist only as an internal consistency check, hence are left untested.
[docs]
def get_loc_self_energy_gamma_abinitio_dga(
gamma_dens: LocalFourPoint, u_loc: LocalInteraction, g_loc: GreensFunction
) -> SelfEnergy:
r"""
DEVELOPMENT / TESTING ONLY. Returns the local self-energy with the density three-leg :math:`\gamma` from
ab-initio DGA (density channel only, by design),
.. math:: \Sigma_{ij}^{\nu} = -\frac{1}{\beta} \sum_\omega U_{iabc}\, \gamma_{cbdj}^{\omega\nu}\, G_{ad}^{\omega-\nu}.
:param gamma_dens: The density three-leg vertex :math:`\gamma_{\mathrm{dens}}` (ab-initio convention).
:param u_loc: The bare local interaction :math:`U`.
:param g_loc: The local :class:`GreensFunction`.
:return: The local :class:`SelfEnergy` (including the Hartree term).
"""
g_wv = g_loc.get_g_wv(MFHelper.wn(config.box.niw_core), config.box.niv_core)
sigma = -1.0 / config.sys.beta * u_loc.times("iabc,cbdjwv,adwv->ijv", gamma_dens.to_full_niw_range(), g_wv)
hartree = u_loc.as_channel(SpinChannel.DENS).times("abcd,dc->ab", config.sys.occ_dmft)[..., None]
return SelfEnergy((sigma + hartree)[None, None, None, ...], beta=config.sys.beta)
