Complex conjugation fix

dgamore/bubble_gen.py

@@ -324,8 +324,8 @@ def create_generalized_chi0_pp_w0(g_dmft: GreensFunction, niv_pp: int, beta: flo
     def create_generalized_chi0_q_pp_w0(giwk: GreensFunction, niv_pp: int, q_grid: KGrid) -> FourPoint:
         r"""
         Returns the momentum-dependent particle-particle bare bubble at :math:`\omega = 0`,
-        :math:`\chi^{\mathrm{k}}_{0;1234} = G^{\mathrm{k}}_{14}\, G^{-\mathrm{k}}_{23}` with :math:`G^{-\mathrm{k}}_{23}
-        = (G^{\mathrm{k}}_{32})^{*}`.
+        :math:`\chi^{\mathrm{k}}_{0;1234} = G^{\mathrm{k}}_{14}\, G^{-\mathrm{k}}_{23}`, the partner propagator read at
+        the flipped momentum and fermionic frequency, :math:`-\mathrm{k} = (-\mathbf{k}, -\nu)`.
         Note that no factor of :math:`-\beta` is included here.
 
         :param giwk: The momentum-dependent :class:`GreensFunction`.
@@ -334,10 +334,9 @@ def create_generalized_chi0_q_pp_w0(giwk: GreensFunction, niv_pp: int, q_grid: K
         :return: The momentum-dependent pp bubble as a :class:`FourPoint` (no bosonic axis, pp notation, compressed q).
         """
         g = giwk.cut_niv(niv_pp).compress_q_dimension()
-        # transpose_orbitals() returns a fresh private copy, so conjugate it in place to reuse its buffer.
-        g_t = g.transpose_orbitals()
-        np.conj(g_t.mat, out=g_t.mat)
-        gchi0_q_pp_w0 = g.mat[:, :, None, None, :, :] * g_t.mat[:, None, :, :, None, :]
+        # G_23 at (-k, -v); the same-momentum conj(G_32(k, v)) equals it only if G(-k) = G(k)
+        g_minus = g.flip_momentum_axis(copy=True)
+        gchi0_q_pp_w0 = g.mat[:, :, None, None, :, :] * g_minus.mat[:, None, :, :, None, ::-1]
 
         return FourPoint(
             gchi0_q_pp_w0, SpinChannel.NONE, q_grid.nk, 0, 1, True, True, True, FrequencyNotation.PP

dgamore/eliashberg_solver.py

@@ -242,9 +242,10 @@ def _write_pp_band(out: np.ndarray, f_chunk: FourPoint, niv_pp: int, omega: np.n
     Writes the :math:`\omega' = 0` pp band of one ladder-vertex window into the pp accumulator of its momenta.
 
     Each bosonic frequency contributes the anti-diagonal :math:`\omega = \nu - \nu'` (see :func:`_pp_w0_band`), and
-    the negative bosonic half is obtained from the positive one through the complex-conjugation symmetry that
-    :meth:`~dgamore.local_n_point.LocalNPoint.to_negative_niw_range` implements. The orbital permutation, the
-    fermionic flip and the overall minus are the ones of :func:`_transform_vertex_frequencies_w0`.
+    the negative bosonic half is obtained from the positive one at the same momentum through the Hermiticity of the
+    vertex, :math:`F^{(\mathbf{q},-\omega)\nu\nu'}_{1234} = (F^{(\mathbf{q},\omega)(-\nu')(-\nu)}_{4321})^*` (a plain
+    conjugation would need :math:`-\mathbf{q}`). The orbital permutation, the fermionic flip and the overall minus
+    are the ones of :func:`_transform_vertex_frequencies_w0`.
 
     :param out: The pp accumulator of this momentum group, shape ``[nq_group, no, no, no, no, 2 niv_pp, 2 niv_pp]``.
     :param f_chunk: The ladder vertex over a bosonic window, in ph notation and half niw range.
@@ -256,7 +257,13 @@ def _write_pp_band(out: np.ndarray, f_chunk: FourPoint, niv_pp: int, omega: np.n
     # cut to the pp box first so the negative-half copy is pp-sized, not core-sized (the cut is centered, so the
     # fermionic flips inside to_negative_niw_range commute with it); positive then mutates the cut copy in place
     cut = f_chunk.cut_niv(niv_pp)
-    negative = cut.to_negative_niw_range().permute_orbitals("abcd->adcb", copy=False).flip_frequency_axis(-1, False)
+    # Hermitian partner F^{-w}_{1234}(v, v') = conj(F^{w}_{4321}(-v', -v)), then the 1432 slot map: orbitals 2341
+    negative = (
+        cut.to_negative_niw_range()
+        .swap_fermionic_frequency_axes(copy=False)
+        .permute_orbitals("bcda->abcd", copy=False)
+        .flip_frequency_axis(-1, False)
+    )
     positive = cut.permute_orbitals("abcd->adcb", copy=False).flip_frequency_axis(-1, False)
 
     for index in range(positive.current_shape[-3]):

dgamore/four_point.py

@@ -397,18 +397,16 @@ def permute_orbitals(self, permutation: str = "abcd->abcd", copy: bool = True) -
         self.mat = np.einsum(permutation, self.mat, optimize=True)
         return self
 
-    def map_to_full_bz(self, grid: KGrid, nq: tuple = None, conjugate: bool = False):
+    def map_to_full_bz(self, k_grid: KGrid, nq: tuple = None):
         """
         Unfolds the object from the irreducible BZ to the full BZ using the grid's symmetry index map (see
         :meth:`IAmNonLocal._map_to_full_bz`), with four orbital dimensions.
 
-        :param grid: The :class:`KGrid` providing the irreducible-to-full BZ index mapping.
+        :param k_grid: The :class:`KGrid` providing the irreducible-to-full BZ index mapping.
         :param nq: Optional number of momenta per direction for the unfolded grid; defaults to the object's ``nq``.
-        :param conjugate: Whether the object holds the complex conjugate of the quantity the grid's orbital rotations
-            were discovered for; the rotation then runs with the conjugate unitaries.
         :return: ``self`` defined on the full BZ.
         """
-        return self._map_to_full_bz(grid, 4, nq, conjugate)
+        return self._map_to_full_bz(k_grid, 4, nq)
 
     def add(self, other, copy: bool = True) -> "FourPoint":
         """

dgamore/greens_function.py

@@ -358,8 +358,11 @@ def get_fill_nonlocal(self) -> tuple[float, np.ndarray, np.ndarray]:
 
         mu_bands: np.ndarray = self._mu * np.eye(self.n_bands)[None, None, None, ...]
 
-        rho_k = _fermi_dirac_density(self._ek.real + smom0 - mu_bands, self._beta)
-        occ_k = rho_k + np.sum(mat.real - g_model.real, axis=-1) / self._beta
+        # a complex (Hermitian) dispersion keeps its imaginary part; a real one keeps the real-part arithmetic
+        ek = np.real_if_close(self._ek)
+        rho_k = _fermi_dirac_density(ek + smom0 - mu_bands, self._beta)
+        box = mat - g_model if np.iscomplexobj(ek) else mat.real - g_model.real
+        occ_k = rho_k + np.sum(box, axis=-1) / self._beta
         occ_k.real[np.abs(occ_k) < 1e-12] = 0.0
 
         occ_mean = np.mean(occ_k, axis=(0, 1, 2))

dgamore/n_point_base.py

@@ -255,6 +255,19 @@ def scale(self, factor, copy: bool = False):
         self.mat *= factor
         return self
 
+    def conj(self, copy: bool = True):
+        """
+        Complex-conjugates the matrix. The in-place branch (``copy=False``) writes the conjugate into the existing
+        buffer and allocates nothing.
+
+        :param copy: If True, operate on and return a deep copy; if False, mutate and return ``self`` in place.
+        :return: The conjugated object (``self`` when ``copy=False``).
+        """
+        if copy:
+            return self.copy().conj(copy=False)
+        np.conj(self.mat, out=self.mat)
+        return self
+
     def __getitem__(self, item):
         """
         Indexing shortcut: ``obj[item]`` is equivalent to ``obj.mat[item]``.
@@ -806,19 +819,17 @@ def interpolate_q_grid(self, nq_new: tuple[int, int, int], copy: bool = True):
 
         return copy.compress_q_dimension() if compress else copy
 
-    def map_to_full_bz(self, k_grid: "KGrid", nq: tuple = None, conjugate: bool = False):
+    def map_to_full_bz(self, k_grid: "KGrid", nq: tuple = None):
         """
         Maps to full BZ using k_grid's inverse map and precomputed orbital rotation tensors.
 
         :param k_grid: The momentum grid carrying the irreducible-to-full-BZ map and per-k orbital rotations.
         :param nq: Optional override for the number of momenta; if None the object's own ``nq`` is used.
-        :param conjugate: Whether the object holds the complex conjugate of the quantity the grid's orbital rotations
-            were discovered for (see :meth:`_map_to_full_bz`).
         :return: ``self`` expanded to the full BZ (four orbital dimensions transformed).
         """
-        return self._map_to_full_bz(k_grid, 4, nq, conjugate)
+        return self._map_to_full_bz(k_grid, 4, nq)
 
-    def _map_to_full_bz(self, k_grid: "KGrid", num_orbital_dimensions: int, nq: tuple = None, conjugate: bool = False):
+    def _map_to_full_bz(self, k_grid: "KGrid", num_orbital_dimensions: int, nq: tuple = None):
         r"""
         Maps the object from the irreducible to the full Brillouin zone.
 
@@ -841,9 +852,6 @@ def _map_to_full_bz(self, k_grid: "KGrid", num_orbital_dimensions: int, nq: tupl
         :param k_grid: The momentum grid carrying the irreducible-to-full-BZ map and per-k orbital rotations.
         :param num_orbital_dimensions: Number of orbital axes to transform; must be 2 or 4.
         :param nq: Optional override for the number of momenta; if None the object's own ``nq`` is used.
-        :param conjugate: If True, the object holds the complex conjugate of the quantity the rotations were
-            discovered for (e.g. a time-reversed kernel), so the rotation runs with the conjugate unitaries
-            :math:`U^*`; this equals conjugating the mapped un-conjugated object.
         :return: ``self`` expanded to the full BZ.
         :raises ValueError: If the object does not have a compressed momentum dimension.
         """
@@ -871,7 +879,7 @@ def _map_to_full_bz(self, k_grid: "KGrid", num_orbital_dimensions: int, nq: tupl
             us = k_grid._auto_us.reshape(np.prod(k_grid.nk), *k_grid._auto_us.shape[3:])
             self.mat = symmetry_reduction.apply_auto_orbital_transform(
                 self.mat,
-                us=us.conj() if conjugate else us,
+                us=us,
                 sigmas=k_grid._auto_sigmas.reshape(-1),
                 conjs=k_grid._auto_conjs.reshape(-1),
                 num_orbital_dimensions=num_orbital_dimensions,

dgamore/nonlocal_sde.py

@@ -742,16 +742,17 @@ def _run_column_sde(
 
     for :math:`\nu \geq 0` as a real-space product (convolution theorem), distributing the kernel's frequency columns
     over the ranks instead of its momenta. A column is one bosonic frequency :math:`\omega` of one fermionic frequency
-    :math:`\nu`; a negative :math:`\omega` is read from the stored positive one by time reversal,
-    :math:`K^{-\omega,\nu} = (K^{\omega,-\nu})^*`. The tasks - runs of
+    :math:`\nu`; a negative :math:`\omega` is read from the stored positive one by time reversal, which for real
+    hoppings is a plain conjugation in real space, :math:`K^{-\omega,\nu}(\mathbf{R}) = (K^{\omega,-\nu}(\mathbf{R}))^*`
+    (in momentum space it maps :math:`\mathbf{q} \to -\mathbf{q}`). The tasks - runs of
     :data:`~dgamore.memory_estimator.SDE_W_BLOCK` consecutive columns of one :math:`\nu` - are laid out by
     :func:`~dgamore.memory_estimator.column_sde_schedule`, and every rank:
 
     1. receives the irreducible-BZ rows of its tasks' columns in rounds of at most ``chunk_bytes``
        (:func:`~dgamore.mpi_utils.transpose_columns`),
-    2. expands each column to the full BZ with :meth:`FourPoint.map_to_full_bz` (a time-reversed column with the
-       conjugate orbital rotation), Fourier transforms it over the momentum axes in the order z, y, x and contracts
-       it with the frequency slab of ``g_r`` at :math:`\nu - \omega` in bounded real-space row chunks,
+    2. expands each column to the full BZ with :meth:`FourPoint.map_to_full_bz`, Fourier transforms it over the
+       momentum axes in the order z, y, x (a time-reversed column is conjugated afterwards) and contracts it with the
+       frequency slab of ``g_r`` at :math:`\nu - \omega` in bounded real-space row chunks,
     3. folds its tasks' partial sums of each :math:`\nu` in block order; the owner of :math:`\nu` (the rank holding
        its first block) folds the other ranks' sums in rank order, and rank 0 gathers the result.
 
@@ -817,9 +818,10 @@ def columns(tasks: range) -> np.ndarray:
                 for w in ws:
                     col = FourPoint(block[..., c : c + 1], SpinChannel.NONE, config.lattice.nk, 1, 0, False, True, True)
                     c += 1
+                    col = col.map_to_full_bz(k_grid).fft(copy=False, axes=(2, 1, 0))
+                    # time reversal of real hoppings conjugates at fixed R: K^{-w,v}(R) = conj(K^{w,-v}(R))
                     if negative:
-                        col = col.to_negative_niw_range()
-                    col = col.map_to_full_bz(k_grid, conjugate=negative).fft(copy=False, axes=(2, 1, 0))
+                        col = col.conj(copy=False)
                     g_slab = g_r[giwk_niv + (w if negative else -w) + v]
                     k_mat = col.mat[..., 0]
                     for r0 in range(0, n_r, rows):
@@ -1214,6 +1216,19 @@ def _assemble_occupation(
     return 2.0 * np.trace(occ).real, occ, occ_k
 
 
+def _has_time_reversal(ek: np.ndarray) -> bool:
+    r"""
+    Returns whether the dispersion obeys the time reversal of real hoppings,
+    :math:`\varepsilon_{12}(-\mathbf{k}) = \varepsilon_{21}(\mathbf{k})`, to 1e-4 of its largest element: above the
+    six-digit precision of wannier90 hopping files, below any physical complex hopping.
+
+    :param ek: The band dispersion ``[kx, ky, kz, o1, o2]``.
+    :return: Whether the relation holds.
+    """
+    flipped = np.roll(np.flip(ek, axis=(0, 1, 2)), 1, axis=(0, 1, 2))
+    return bool(np.allclose(flipped, np.swapaxes(ek, -1, -2), rtol=0.0, atol=1e-4 * np.abs(ek).max()))
+
+
 def calculate_sigma_proposal(
     sigma_in: SelfEnergy,
     mu: float,
@@ -1234,7 +1249,12 @@ def calculate_sigma_proposal(
     :math:`\mu`: Hartree/Fock, the Dyson Green's function, the bubble, the double-counting, density and magnetic
     kernels, and the FFT Schwinger-Dyson contraction, finished with the noise-removal term and the DMFT tail. The
     tail carries the momentum-dependent part of the Hartree-Fock term, i.e. the contribution of
-    :math:`V^{\mathbf{q}}`, which the impurity self-energy does not contain.
+    :math:`V^{\mathbf{q}}`, which the impurity self-energy does not contain. For several orbitals on a lattice with
+    real hoppings the proposal is finally averaged with its time-reversed partner (see
+    :meth:`~dgamore.self_energy.SelfEnergy.symmetrize_time_reversal`): the Schwinger-Dyson equation is one-sided,
+    with the three-leg vertex on the first external orbital and the interaction on the second, time reversal maps it
+    onto the mirrored form, and the average is the two-sided equation. One band needs no average, its one-sided
+    equation already obeys time reversal.
 
     Single source of truth for the proposal map: it is called once per self-consistency iteration by
     :func:`calculate_self_energy_q`. The local irreducible vertex is frozen, so every
@@ -1408,6 +1428,12 @@ def calculate_sigma_proposal(
     # calculated in this code and add the smooth dmft self-energy
     sigma_prop += delta_sigma
     sigma_prop = sigma_prop.concatenate_self_energies(sigma_dmft, shell_offset=hf_v)
+    # the one-sided SDE breaks Sigma(-k) = Sigma(k)^T for several orbitals; the average is the two-sided SDE
+    if sigma_prop.n_bands > 1 and _has_time_reversal(config.lattice.hamiltonian.get_ek()):
+        removed = sigma_prop.symmetrize_time_reversal()
+        logger.info(f"Time-reversal average of the self-energy proposal: largest change {removed:.3e}.")
+    elif sigma_prop.n_bands > 1:
+        logger.info("The dispersion breaks time reversal (complex hoppings); the self-energy proposal is not averaged.")
     return sigma_prop
 
 

dgamore/self_energy.py

@@ -33,6 +33,10 @@ class SelfEnergy(TwoPoint):
     Matsubara frequencies.
     """
 
+    # Upper bound (in bytes) on the transient of symmetrize_time_reversal: the gathered partners and their average of
+    # one fermionic chunk, so a full-BZ self-energy is never duplicated.
+    _TR_CHUNK_BYTES = 256 * 1024**2
+
     def __init__(
         self,
         mat: np.ndarray,
@@ -289,6 +293,38 @@ def concatenate_self_energies(self, other: "SelfEnergy", shell_offset: np.ndarra
             result_mat, self.nq, self.full_niv_range, self.has_compressed_q_dimension, False, beta=self._beta
         )
 
+    def symmetrize_time_reversal(self) -> float:
+        r"""
+        Enforces the time-reversal relation of real hoppings, :math:`\Sigma^{(\mathbf{k},\nu)}_{12} =
+        \Sigma^{(-\mathbf{k},\nu)}_{21}`, in place by averaging every momentum with its time-reversed partner,
+
+        .. math:: \Sigma^{(\mathbf{k},\nu)}_{12} \to \tfrac{1}{2}\left(\Sigma^{(\mathbf{k},\nu)}_{12}
+            + \Sigma^{(-\mathbf{k},\nu)}_{21}\right).
+
+        The average keeps the Matsubara Hermiticity and the causality of the input. It runs over chunks of the
+        fermionic axis whose transient (the gathered partners and their average) stays below ``_TR_CHUNK_BYTES``.
+
+        :return: The largest change of an element, i.e. half the largest time-reversal asymmetry of the input.
+        :raises ValueError: If the momentum axis is not compressed.
+        """
+        if not self.has_compressed_q_dimension:
+            raise ValueError("The time-reversal average needs a compressed momentum axis.")
+        grid = np.array(self.nq)
+        k = np.array(np.unravel_index(np.arange(grid.prod()), self.nq))
+        minus = np.ravel_multi_index(tuple(-k % grid[:, None]), self.nq)
+
+        mat = self.mat
+        width = max(1, int(self._TR_CHUNK_BYTES // (2 * mat[..., :1].nbytes)))
+        removed = 0.0
+        for start in range(0, mat.shape[-1], width):
+            block = mat[..., start : start + width]
+            average = block[minus].swapaxes(1, 2)
+            average += block
+            average *= 0.5
+            removed = max(removed, float(np.abs(average - block).max()))
+            block[...] = average
+        return removed
+
     def fit_smom_concatenated(
         self, other: "SelfEnergy", shell_offset: np.ndarray | None = None
     ) -> tuple[np.ndarray, np.ndarray]:

tests/test_bubble_gen.py

@@ -147,22 +147,41 @@ def test_create_generalized_chi0_q_pp_w0_matches_reference():
 
     gm = g.cut_niv(niv_pp).compress_q_dimension().mat
     nkt, n = gm.shape[0], 2 * niv_pp
+    minus_k = [np.ravel_multi_index(tuple(-np.array(np.unravel_index(k, nk)) % nk), nk) for k in range(nkt)]
     ref = np.zeros((nkt, nb, nb, nb, nb, n), dtype=np.complex128)
     for k in range(nkt):
         for a in range(nb):
             for b in range(nb):
                 for c in range(nb):
                     for d in range(nb):
-                        # G_14^{kv} * conj(G_32^{kv})   (transpose_orbitals -> g[c,b], conjugated)
-                        ref[k, a, b, c, d, :] = gm[k, a, d, :] * np.conj(gm[k, c, b, :])
+                        # G_14^{k,v} * G_23^{-k,-v}
+                        ref[k, a, b, c, d, :] = gm[k, a, d, :] * gm[minus_k[k], b, c, ::-1]
 
     assert res.mat.shape == ref.shape
     assert res.frequency_notation == FrequencyNotation.PP
     assert np.allclose(res.mat, ref, atol=1e-4)
 
 
+def test_momentum_pp_bubble_transforms_like_a_four_point_object_under_a_unit_cell_relabeling():
+    """Relabeling one orbital's cell (G -> D G D^+, D(-k) = D(k)^*) rephases the pp bubble by D_1 D_2^* D_3 D_4^*."""
+    nk, nb, niv_pp = (4, 1, 1), 2, 3
+    q_grid = bz.KGrid(nk, symmetries=[])
+    g = _make_momentum_g(nk, nb, niv_pp, seed=12)
+    d = np.ones((*nk, nb), dtype=np.complex128)
+    d[..., 1] = np.exp(2j * np.pi * np.arange(nk[0]) / nk[0])[:, None, None]
+    mat = d[..., :, None, None] * g.mat * d[..., None, :, None].conj()
+    shifted = GreensFunction(mat, has_compressed_q_dimension=False, nk=nk)
+    bubble, bubble_shifted = (
+        BubbleGenerator.create_generalized_chi0_q_pp_w0(x, niv_pp, q_grid).decompress_q_dimension().mat
+        for x in (g, shifted)
+    )
+    dc = d.conj()
+    phase = d[..., :, None, None, None] * dc[..., None, :, None, None] * d[..., None, None, :, None]
+    assert np.allclose(bubble_shifted, (phase * dc[..., None, None, None, :])[..., None] * bubble, atol=1e-5)
+
+
 def test_create_generalized_chi0_q_pp_w0_does_not_mutate_input():
-    """create_generalized_chi0_q_pp_w0 conjugates only a private copy and leaves its input untouched."""
+    """create_generalized_chi0_q_pp_w0 flips only a private copy and leaves its input untouched."""
     nk, nb, niv_pp = (2, 2, 1), 2, 2
     g = _make_momentum_g(nk, nb, niv_pp + 1, seed=7)
     g_before = g.mat.copy()