BUG: Pulay mixing is possibly not behaving as intended

moldga/__init__.py

@@ -0,0 +1,6 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+

moldga/brillouin_zone.py

@@ -1,6 +1,11 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 """
-Module to handle operations within the (irreducible) Brillouin zone. Copied over from Paul Worm's code.
-Only modified the constant arrays and made enums out of them for type hinting.
+Module to handle operations within the (irreducible) Brillouin zone. Heavily inspired by DGApy.
 """
 
 import warnings

moldga/bubble_gen.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import moldga.config as config
 from moldga.greens_function import GreensFunction
 from moldga.matsubara_frequencies import MFHelper

moldga/config.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import numpy as np
 
 import moldga.brillouin_zone as bz

moldga/config_parser.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import argparse
 import os
 
@@ -207,8 +213,8 @@ def _build_self_energy_interpolation_config(self, conf_file):
         section = conf_file["self_energy_interpolation"]
 
         conf.do_interpolation = self._try_parse(section, "do_interpolation", False)
-        conf.beta_target = self._try_parse(section, "beta_target", 1.0)
-        conf.niv_target = self._try_parse(section, "niv_target", 10)
+        conf.beta_target = self._try_parse(section, "target_beta", 1.0)
+        conf.niv_target = self._try_parse(section, "target_niv", 10)
 
         return conf
 

moldga/dga_io.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import os
 
 import moldga.brillouin_zone as bz

moldga/dga_logger.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import logging
 import os
 from datetime import datetime

moldga/dga_main.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import itertools as it
 import logging
 import os

moldga/eliashberg_solver.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import os
 from typing import Tuple
 

moldga/four_point.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import gc
 from copy import deepcopy
 

moldga/gap_function.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import numpy as np
 
 from moldga.brillouin_zone import KGrid

moldga/greens_function.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 from copy import deepcopy
 
 import numpy as np

moldga/hamiltonian.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import itertools as it
 import logging
 

moldga/interaction.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 from copy import deepcopy
 
 import numpy as np

moldga/lambda_correction.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import numpy as np
 
 from moldga import config

moldga/local_four_point.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 from moldga.interaction import LocalInteraction, Interaction
 from moldga.local_n_point import LocalNPoint
 from moldga.matsubara_frequencies import MFHelper

moldga/local_n_point.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import itertools
 import os
 

moldga/local_sde.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import moldga.config as config
 from moldga.bubble_gen import BubbleGenerator
 from moldga.greens_function import GreensFunction

moldga/matsubara_frequencies.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 from enum import Enum
 
 import numpy as np

moldga/mpi_distributor.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import gc
 import os
 

moldga/n_point_base.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import gc
 from abc import ABC
 from copy import deepcopy
@@ -490,29 +496,54 @@ def _map_to_full_bz(self, k_grid: KGrid, num_orbital_dimensions: int, nq: tuple
                 if np.allclose(u_ref, identity):
                     continue
 
-                uc_ref = u_ref.conj()
                 idx = np.array(indices)
 
-                if num_orbital_dimensions == 2:
-                    if path_2 is None:
-                        path_2 = np.einsum_path(
-                            "ap,bq,kpq...->kab...", u_ref, uc_ref, self.mat[idx], optimize="optimal"
-                        )[0]
-                    self.mat[idx] = np.einsum("ap,bq,kpq...->kab...", u_ref, uc_ref, self.mat[idx], optimize=path_2)
-                elif num_orbital_dimensions == 4:
-                    if path_4 is None:
-                        path_4 = np.einsum_path(
+                def _is_permutation_matrix(u: np.ndarray) -> bool:
+                    return (
+                        np.allclose(np.abs(u), np.abs(u).astype(int))  # only 0s and 1s
+                        and np.allclose(u.sum(axis=0), 1)  # one 1 per column
+                        and np.allclose(u.sum(axis=1), 1)  # one 1 per row
+                    )
+
+                if _is_permutation_matrix(u_ref):
+                    # For real permutation matrices, U @ M @ U^T is just index reordering.
+                    # Find the permutation: perm[i] = j means orbital i gets content from orbital j.
+                    perm = np.argmax(u_ref.real, axis=1)
+
+                    if num_orbital_dimensions == 2:
+                        self.mat[idx] = self.mat[idx][:, perm, ...][:, :, perm, ...]
+                    elif num_orbital_dimensions == 4:
+                        self.mat[idx] = self.mat[idx][:, perm, ...][:, :, perm, ...][:, :, :, perm, ...][
+                            :, :, :, :, perm, ...
+                        ]
+                else:
+                    uc_ref = u_ref.conj()
+                    if num_orbital_dimensions == 2:
+                        if path_2 is None:
+                            path_2 = np.einsum_path(
+                                "ap,bq,kpq...->kab...", u_ref, uc_ref, self.mat[idx], optimize="optimal"
+                            )[0]
+                        self.mat[idx] = np.einsum("ap,bq,kpq...->kab...", u_ref, uc_ref, self.mat[idx], optimize=path_2)
+                    elif num_orbital_dimensions == 4:
+                        if path_4 is None:
+                            path_4 = np.einsum_path(
+                                "ap,bq,cr,ds,kpqrs...->kabcd...",
+                                u_ref,
+                                uc_ref,
+                                u_ref,
+                                uc_ref,
+                                self.mat[idx],
+                                optimize="optimal",
+                            )[0]
+                        self.mat[idx] = np.einsum(
                             "ap,bq,cr,ds,kpqrs...->kabcd...",
                             u_ref,
                             uc_ref,
                             u_ref,
                             uc_ref,
                             self.mat[idx],
-                            optimize="optimal",
-                        )[0]
-                    self.mat[idx] = np.einsum(
-                        "ap,bq,cr,ds,kpqrs...->kabcd...", u_ref, uc_ref, u_ref, uc_ref, self.mat[idx], optimize=path_4
-                    )
+                            optimize=path_4,
+                        )
 
         self.mat = self.mat.reshape((np.prod(self.nq), *self.original_shape[1:]))
         self.update_original_shape()

moldga/nonlocal_sde.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import gc
 import glob
 import os
@@ -704,11 +710,11 @@ def apply_mixing_strategy(
     is either 'linear' or 'pulay'.
     """
     logger = config.logger
-    strategy = config.self_consistency.mixing_strategy
     n_hist = config.self_consistency.mixing_history_length
+    alpha = config.self_consistency.mixing
 
     if (
-        strategy == "pulay"
+        config.self_consistency.mixing_strategy == "pulay"
         and current_iter > n_hist
         and config.self_consistency.save_iter
         and config.output.save_quantities
@@ -720,10 +726,10 @@ def apply_mixing_strategy(
         last_proposals = [sigma_dmft_stacked] + last_results
         last_results = last_results + [sigma_new.mat]
 
-        niv_dmft = sigma_new.current_shape[-1] // 2
+        niv = sigma_new.current_shape[-1] // 2
         niv_core = config.box.niv_core
-        last_proposals = [sigma[..., niv_dmft - niv_core : niv_dmft + niv_core] for sigma in last_proposals]
-        last_results = [sigma[..., niv_dmft - niv_core : niv_dmft + niv_core] for sigma in last_results]
+        last_proposals = [sigma[..., niv - niv_core : niv + niv_core] for sigma in last_proposals]
+        last_results = [sigma[..., niv - niv_core : niv + niv_core] for sigma in last_results]
         logger.info(f"Loaded last {n_hist} self-energies from files.")
 
         shape = last_results[-1].shape
@@ -749,17 +755,20 @@ def get_result(idx: int):
 
             f_matrix[:, i] -= r_matrix[:, i]
 
+        # Residual: F(x_n) - x_n, where x_n = last_proposals[-1] = sigma_old (core window)
         iter_diff = get_result(-1) - get_proposal(-1)
         f_i[:n_total] = iter_diff.real
         f_i[n_total:] = iter_diff.imag
 
-        update = config.self_consistency.mixing * f_i
-        fact1 = (r_matrix + config.self_consistency.mixing * f_matrix) @ np.linalg.inv(f_matrix.T @ f_matrix)
-        update -= fact1 @ (f_matrix.T @ f_i)
+        # Solve min||F @ c - f_i|| via least squares (more stable than explicit inverse)
+        coeffs, _, _, _ = np.linalg.lstsq(f_matrix, f_i, rcond=None)
+
+        # Pulay update: x_{n+1} = x_n + alpha*f_i - (R + alpha*F) @ c
+        update = alpha * f_i - (r_matrix + alpha * f_matrix) @ coeffs
         update = update[:n_total] + 1j * update[n_total:]
-        sigma_new.mat[..., niv_dmft - niv_core : niv_dmft + niv_core] = sigma_old.compress_q_dimension().mat[
-            ..., niv_dmft - niv_core : niv_dmft + niv_core
-        ] + update.reshape(shape)
+
+        # Update the new self energy
+        sigma_new.mat[..., niv - niv_core : niv + niv_core] = get_proposal(-1).reshape(shape) + update.reshape(shape)
 
         logger.info(
             f"Pulay mixing applied with {n_hist} previous iterations and a mixing parameter of {config.self_consistency.mixing}."

moldga/plotting.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import os
 
 import numpy as np

moldga/self_energy.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import itertools
 import itertools as it
 from copy import deepcopy

moldga/symmetrize_new.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import gc
 import glob
 import itertools as it

moldga/w2dyn_aux.py

@@ -1,3 +1,9 @@
+# SPDX-FileCopyrightText: 2025-2026 Julian Peil <julian.peil@tuwien.ac.at>
+# SPDX-License-Identifier: MIT
+#
+# moLDGA — Multi-Orbital Ladder Dynamical Vertex Approximation (LDGA) &
+#          Eliashberg Equation Solver for Strongly Correlated Electron Systems
+
 import h5py
 
 import moldga.symmetrize_new as sym