diff --git a/dgamore/bubble_gen.py b/dgamore/bubble_gen.py index 8e96aa2..4b30bda 100644 --- a/dgamore/bubble_gen.py +++ b/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 diff --git a/dgamore/eliashberg_solver.py b/dgamore/eliashberg_solver.py index cd25dd0..c7d9978 100644 --- a/dgamore/eliashberg_solver.py +++ b/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]): diff --git a/dgamore/four_point.py b/dgamore/four_point.py index 7a709d7..0b3728f 100644 --- a/dgamore/four_point.py +++ b/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": """ diff --git a/dgamore/greens_function.py b/dgamore/greens_function.py index 7875dc2..233d15c 100644 --- a/dgamore/greens_function.py +++ b/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)) diff --git a/dgamore/n_point_base.py b/dgamore/n_point_base.py index 960bc8a..684fb1a 100644 --- a/dgamore/n_point_base.py +++ b/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, diff --git a/dgamore/nonlocal_sde.py b/dgamore/nonlocal_sde.py index b759617..4bb64a0 100644 --- a/dgamore/nonlocal_sde.py +++ b/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 diff --git a/dgamore/self_energy.py b/dgamore/self_energy.py index e435eab..bc0ad86 100644 --- a/dgamore/self_energy.py +++ b/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]: diff --git a/tests/test_bubble_gen.py b/tests/test_bubble_gen.py index 55e6275..9c565ae 100644 --- a/tests/test_bubble_gen.py +++ b/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() diff --git a/tests/test_eliashberg_end_to_end.py b/tests/test_eliashberg_end_to_end.py index 3cbcd06..66f4553 100644 --- a/tests/test_eliashberg_end_to_end.py +++ b/tests/test_eliashberg_end_to_end.py @@ -213,6 +213,8 @@ def test_eliashberg_gap_functions_carry_sector_parity_and_match_reference(setup) config.output.eliashberg_path = folder config.eliashberg.n_eig = 4 config.eliashberg.resolve_frequency_parity = True + # at the default 1e-6 the admixture of the nearest gap (~1e-3) depends on the BLAS platform + config.eliashberg.epsilon = 1e-10 u_loc = config.lattice.hamiltonian.get_local_u() v_nonloc = config.lattice.hamiltonian.get_vq(config.lattice.k_grid) diff --git a/tests/test_eliashberg_solver.py b/tests/test_eliashberg_solver.py index d50969a..e363568 100644 --- a/tests/test_eliashberg_solver.py +++ b/tests/test_eliashberg_solver.py @@ -2733,6 +2733,24 @@ def test_slice_constructor_is_bit_invariant_under_the_chunk_budget(setup, monkey assert np.array_equal(one_chunk.mat, per_w.mat) +def test_pp_band_reads_the_negative_bosonic_half_by_hermiticity(): + """The w < 0 pp band is -conj(F^{|w|}_{2341}(v', -v)), the same-momentum Hermitian partner of the ladder vertex.""" + no, niv_pp, nq = 2, 2, 3 + n2 = 2 * niv_pp + shape = (nq,) + (no,) * 4 + (n2, n2, n2) + rng = np.random.default_rng(8) + f = rng.standard_normal(shape) + 1j * rng.standard_normal(shape) + _, omega = es._pp_w0_band(niv_pp, n2 - 1) + out = np.zeros((nq,) + (no,) * 4 + (n2, n2), dtype=np.complex128) + es._write_pp_band(out, FourPoint(f.copy(), SpinChannel.DENS, (nq, 1, 1), 1, 2, False, True, True), niv_pp, omega, 0) + f_1432, f_2341 = np.einsum("qadcbwxy->qabcdwxy", f), np.einsum("qbcdawxy->qabcdwxy", f) + ref = np.empty_like(out) + for i, j in np.ndindex(n2, n2): + w = omega[i, j] + ref[..., i, j] = -f_1432[..., w, i, n2 - 1 - j] if w >= 0 else -np.conj(f_2341[..., -w, j, n2 - 1 - i]) + assert np.allclose(out, ref, atol=1e-5) + + def test_streaming_fq_file_has_gather_layout_and_matches_pp_band(setup): """The streamed f_irrq file carries the historical layout and the w' = 0 band map at a spot-checked entry.""" no = 2 diff --git a/tests/test_greens_function.py b/tests/test_greens_function.py index 9c1111e..e0f13a9 100644 --- a/tests/test_greens_function.py +++ b/tests/test_greens_function.py @@ -196,6 +196,47 @@ def test_energies_use_injected_state_not_config(monkeypatch): assert np.isfinite(g.get_epot()) +def _make_complex_hopping_ek(nk: tuple) -> np.ndarray: + """Builds a two-site dispersion whose inter-site hopping ``0.3 + 0.7 e^{ik_x}`` is complex, with H(-k) = H(k)^*.""" + kx = 2 * np.pi * np.arange(nk[0]) / nk[0] + ek = np.zeros((*nk, 2, 2), dtype=np.complex128) + ek[..., 0, 0] = -2 * np.cos(kx)[:, None, None] + ek[..., 1, 1] = 0.5 + ek[..., 0, 1] = (0.3 + 0.7 * np.exp(1j * kx))[:, None, None] + ek[..., 1, 0] = np.conj(ek[..., 0, 1]) + return ek + + +def test_get_fill_nonlocal_is_the_fermi_dirac_density_of_a_complex_hopping_dispersion(): + """Without a self-energy the k-resolved occupation is f(beta (H(k) - mu)) of the complex H(k), for any box.""" + nk, niv, beta, mu = (4, 1, 1), 8, 5.0, 0.2 + ek = _make_complex_hopping_ek(nk) + sig = SelfEnergy(np.zeros((int(np.prod(nk)), 2, 2, 2 * niv)), nk=nk, has_compressed_q_dimension=True, beta=beta) + _, _, occ_k = GreensFunction.get_g_full(sig, mu, ek, beta).get_fill_nonlocal() + eps, vecs = np.linalg.eigh(beta * (ek - mu * np.eye(2))) + exact = (vecs / (1 + np.exp(eps))[..., None, :]) @ np.conj(np.swapaxes(vecs, -1, -2)) + assert np.allclose(occ_k, exact, atol=1e-5) + + +def test_get_fill_nonlocal_keeps_the_real_arithmetic_for_a_real_dispersion(): + """A real dispersion is occupied with the real-part arithmetic f(Re h) + Re(G - G_model), bit for bit.""" + from dgamore.greens_function import _fermi_dirac_density + + nk, nb, niv, beta, mu = (2, 2, 1), 2, 6, 5.0, 0.3 + rng = np.random.default_rng(3) + ek = rng.standard_normal((*nk, nb, nb)) + ek = 0.5 * (ek + ek.swapaxes(-1, -2)) + 0j + shape = (int(np.prod(nk)), nb, nb, 2 * niv) + sig_mat = rng.standard_normal(shape) + 1j * rng.standard_normal(shape) + sig = SelfEnergy(sig_mat, nk=nk, has_compressed_q_dimension=True, beta=beta) + g = GreensFunction.get_g_full(sig, mu, ek, beta) + _, _, occ_k = g.get_fill_nonlocal() + box = np.sum(g._get_gfull_mat().real - g._get_g_model_k_mat().real, axis=-1) / beta + ref = _fermi_dirac_density(ek.real + sig.smom[0][None, None, None] - mu * np.eye(nb), beta) + box + ref.real[np.abs(ref) < 1e-12] = 0.0 + assert np.array_equal(occ_k, ref) + + def test_fermi_dirac_density_matches_diagonal_matmul_reference(): """_fermi_dirac_density (column-scaling) equals the explicit V diag(f) V^-1 construction, bit-for-bit.""" from dgamore.greens_function import _fermi_dirac_density diff --git a/tests/test_n_point_base.py b/tests/test_n_point_base.py index 4437acf..91fff46 100644 --- a/tests/test_n_point_base.py +++ b/tests/test_n_point_base.py @@ -95,6 +95,25 @@ def test_raises_error_when_dividing_by_invalid_type(): obj / "invalid" +def test_conj_in_place_conjugates_the_stored_buffer_and_returns_self(): + """conj(copy=False) writes the conjugate into the existing buffer and returns the same object.""" + mat = np.array([[1 + 2j, 3 - 4j], [-5j, 6]]) + obj = IHaveMat(mat) + buffer = obj.mat + result = obj.conj(copy=False) + assert result is obj and result.mat is buffer + assert np.array_equal(result.mat, np.conj(mat)) + + +def test_conj_with_copy_returns_the_conjugate_and_leaves_the_original_untouched(): + """conj() returns a conjugated copy and leaves the original matrix unchanged.""" + mat = np.array([[1 + 2j, 3 - 4j], [-5j, 6]]) + obj = IHaveMat(mat) + result = obj.conj() + assert result is not obj + assert np.array_equal(result.mat, np.conj(mat)) and np.array_equal(obj.mat, mat) + + def test_reshapes_matrix_and_updates_original_shape(): """Reshaping updates the matrix and tracks the original shape.""" mat = np.array([[1, 2], [3, 4]]) @@ -1386,24 +1405,6 @@ def test_map_to_full_bz_plain_symmetry_kgrid_does_not_call_orbital_transform(c12 assert spy.call_count == 0 -@pytest.mark.parametrize("antiunitary", [False, True]) -def test_map_to_full_bz_of_a_conjugated_object_with_conjugate_rotation_is_the_conjugated_mapping( - c128_storage, antiunitary -): - """Mapping conj(X) with conjugate unitaries equals conj of mapping X bit for bit, anti-unitary members included.""" - grid, _ = _build_auto_kgrid(nx=4, ny=4, nz=2, nb=2) - rng = np.random.default_rng(4) - grid._auto_us = np.array( - [np.linalg.qr(rng.standard_normal((2, 2)) + 1j * rng.standard_normal((2, 2)))[0] for _ in range(32)] - ).reshape(4, 4, 2, 2, 2) - grid._auto_conjs = (rng.random((4, 4, 2)) < 0.5) if antiunitary else np.zeros((4, 4, 2), dtype=bool) - shape = (grid.nk_irr, 2, 2, 2, 2, 3) - x = rng.standard_normal(shape) + 1j * rng.standard_normal(shape) - plain = IAmNonLocal(x.copy(), (4, 4, 2), has_compressed_q_dimension=True)._map_to_full_bz(grid, 4) - conj = IAmNonLocal(np.conj(x), (4, 4, 2), has_compressed_q_dimension=True)._map_to_full_bz(grid, 4, conjugate=True) - assert np.array_equal(conj.mat, np.conj(plain.mat)) - - def test_auto_orbital_groups_are_cached_per_dtype_and_recomputed_for_replaced_rotations(): """KGrid.auto_orbital_groups computes the groups once per dtype and again after the rotations are replaced.""" grid, _ = _build_auto_kgrid(nx=4, ny=4, nz=2, nb=2) diff --git a/tests/test_nonlocal_sde.py b/tests/test_nonlocal_sde.py index 940d21e..2b200fb 100644 --- a/tests/test_nonlocal_sde.py +++ b/tests/test_nonlocal_sde.py @@ -1554,15 +1554,47 @@ def fn(comm, rank): return res[0] +def _time_reversed_qloop_sigma(kernel, giwk): + """Hand-rolled q-loop sigma over the full-BZ kernel, its w < 0 half read at -q by time reversal (conj, v -> -v).""" + k_grid, nk = config.lattice.k_grid, config.lattice.nk + niw, niv, o = config.box.niw_core, config.box.niv_core, config.sys.n_bands + full = FourPoint(kernel.copy(), SpinChannel.NONE, nk, 1, 1, False, True, True).map_to_full_bz(k_grid).mat + q_list = k_grid.get_q_list() + minus_q = np.ravel_multi_index(tuple((-q_list % np.array(nk)).T), nk) + kernel_w = np.concatenate((np.conj(full[minus_q][..., 1:, ::-1])[..., ::-1, :], full), axis=-2) + mat = np.zeros((*nk, o, o, niv), dtype=np.complex128) + for iq, q in enumerate(q_list): + g_shift = np.roll(giwk.mat, tuple(q), axis=(0, 1, 2)) + for iw, w in enumerate(MFHelper.wn(niw)): + g = g_shift[..., giwk.niv - w : giwk.niv + niv - w] + mat += np.einsum("aijdv,xyzadv->xyzijv", kernel_w[iq, ..., iw, niv:], g) + mat *= -0.5 / config.sys.beta / k_grid.nk_tot + return SelfEnergy(mat, nk, False, beta=config.sys.beta).compress_q_dimension().to_full_niv_range() + + +def test_has_time_reversal_accepts_real_hoppings_and_rejects_a_complex_one(): + """Real hoppings obey eps(-k) = eps(k)^T up to hopping-file noise; a complex on-site hopping breaks it.""" + kx = 2 * np.pi * np.arange(4) / 4 + ek = np.zeros((4, 1, 1, 2, 2), dtype=np.complex128) + ek[..., 0, 0] = -2 * np.cos(kx)[:, None, None] + ek[..., 0, 1] = (0.3 + 0.7 * np.exp(1j * kx))[:, None, None] + ek[..., 1, 0] = np.conj(ek[..., 0, 1]) + assert nonlocal_sde._has_time_reversal(ek) + ek[..., 0, 1] += 1e-6j + ek[..., 1, 0] -= 1e-6j + assert nonlocal_sde._has_time_reversal(ek) + ek[..., 0, 1] += 0.2j + ek[..., 1, 0] -= 0.2j + assert not nonlocal_sde._has_time_reversal(ek) + + @pytest.mark.parametrize("auto", [False, True]) @pytest.mark.parametrize("size", [1, 3]) def test_column_sde_matches_the_qloop_reference(auto, size, monkeypatch): - """The column-distributed contraction reproduces the q-loop sum over the kernel mapped to the full BZ.""" + """The column-distributed contraction reproduces the q-loop sum, its w < 0 half read at -q by time reversal.""" monkeypatch.setattr(mpi_utils, "MPI", FAKE_MPI) kernel, giwk = _column_sde_setup(auto) - k_grid = config.lattice.k_grid - full_kernel = FourPoint(kernel.copy(), SpinChannel.NONE, config.lattice.nk, 1, 1, False, True, True) - ref = nonlocal_sde.calculate_sigma_from_kernel(full_kernel.map_to_full_bz(k_grid), giwk, k_grid.get_q_list()) + ref = _time_reversed_qloop_sigma(kernel, giwk) mat = _run_column_sde_parallel(size, kernel, giwk, hostnames=["n0"] * size) sigma = SelfEnergy(mat, config.lattice.nk, False, True, calc_smom=False, beta=config.sys.beta) diff --git a/tests/test_self_energy.py b/tests/test_self_energy.py index 84eaeaf..44faa0c 100644 --- a/tests/test_self_energy.py +++ b/tests/test_self_energy.py @@ -826,3 +826,44 @@ def test_fit_smom_of_an_orbitally_symmetric_self_energy_equals_the_elementwise_f assert mom0.dtype.kind == "f" and mom1.dtype.kind == "f" assert np.array_equal(mom0, np.mean(fitdata.real, axis=-1)) assert np.array_equal(mom1, np.mean(fitdata.imag * vn, axis=-1)) + + +def _time_reversed(mat: np.ndarray, grid: tuple) -> np.ndarray: + """Returns Sigma(-k)^T of a compressed [k, o1, o2, v] array, built from per-axis momentum flips.""" + flipped = np.roll(np.flip(mat.reshape(*grid, *mat.shape[1:]), axis=(0, 1, 2)), 1, axis=(0, 1, 2)) + return np.swapaxes(flipped, 3, 4).reshape(mat.shape) + + +def test_symmetrize_time_reversal_averages_every_momentum_with_its_transposed_partner(): + """symmetrize_time_reversal sets Sigma(k) to (Sigma(k) + Sigma(-k)^T) / 2 and returns the removed asymmetry.""" + rng = np.random.default_rng(5) + grid = (4, 3, 2) + se = _se( + rng.standard_normal((24, 3, 3, 4)) + 1j * rng.standard_normal((24, 3, 3, 4)), + nk=grid, + has_compressed_q_dimension=True, + ) + old = se.mat.copy() + expected = 0.5 * (old + _time_reversed(old, grid)) + removed = se.symmetrize_time_reversal() + assert np.allclose(se.mat, expected, atol=1e-6) + assert np.allclose(se.mat, _time_reversed(se.mat, grid), atol=1e-6) + assert np.isclose(removed, np.abs(expected - old).max(), atol=1e-6) + + +def test_symmetrize_time_reversal_is_bit_invariant_under_the_chunk_size(monkeypatch): + """The frequency-chunked averaging reproduces the single-pass result bit for bit.""" + rng = np.random.default_rng(6) + mat = rng.standard_normal((16, 2, 2, 10)) + 1j * rng.standard_normal((16, 2, 2, 10)) + single = _se(mat.copy(), nk=(4, 4, 1), has_compressed_q_dimension=True) + single.symmetrize_time_reversal() + monkeypatch.setattr(SelfEnergy, "_TR_CHUNK_BYTES", 1) + chunked = _se(mat.copy(), nk=(4, 4, 1), has_compressed_q_dimension=True) + chunked.symmetrize_time_reversal() + assert np.array_equal(chunked.mat, single.mat) + + +def test_symmetrize_time_reversal_rejects_an_uncompressed_momentum_axis(): + """symmetrize_time_reversal needs the compressed momentum axis.""" + with pytest.raises(ValueError): + _se(mat_decompressed.copy(), nk=nk, has_compressed_q_dimension=False).symmetrize_time_reversal()