Energy calculation and bugfix in Eliashberg equation

README.md

@@ -133,7 +133,6 @@ mpirun/mpiexec -np $SLURM_NTASKS DGAmore.py -p "<path to config>" -c "<name of c
 
 The `#SBATCH` options `-o` and `-e` denote the file where the job output and errors, respectively, will be written. In
 this case both will be written into the same file, however, it is possible to use separate files.
-The Python option `-u` tells the program to flush `stdout` and `stderr` in real time.
 
 ---
 
@@ -142,11 +141,11 @@ The Python option `-u` tells the program to flush `stdout` and `stderr` in real
 To configure run parameters, one has to create or modify the existing configuration file, which is a YAML file. The 
 default configuration file is [dga_config.yaml](dgamore/dga_config.yaml), which can be found in the repository directory. 
 Each entry in the configuration file is explained in detail in the file itself, so please refer to the file for more 
-information on how to configure the code. Please note that each section in the configuration file is required, and the 
-code will not run if any section is missing. The default values in the configuration file are set to reasonable 
+information on how to configure the code. The default values in the configuration file are set to reasonable 
 defaults, so if you are unsure about any of the parameters, you can simply use the default values (except for the file 
-paths to the input). If you are unsure about any of the parameters, please refer to the documentation in the 
-configuration file, my [Master's thesis](https://doi.org/10.34726/hss.2025.130528) or contact me via [E-Mail](mailto:julian.peil@tuwien.ac.at).
+paths to the input). If a setting or if a whole section is missing from the file, the code will use the default values. 
+If you are unsure about any of the parameters, please refer to the documentation in the configuration file, my 
+[Master's thesis](https://doi.org/10.34726/hss.2025.130528) or contact me via [E-Mail](mailto:julian.peil@tuwien.ac.at).
 
 --- 
 

dgamore/DGAmore.py

@@ -257,6 +257,7 @@ def write_smom(obj_full: SelfEnergy, obj_ineq: SelfEnergy, sl: slice):
     if config.lambda_correction.perform_lambda_correction and comm.rank == 0:
         chi_d_full.save(name="chi_dens_loc", output_dir=config.output.output_path)
         chi_m_full.save(name="chi_magn_loc", output_dir=config.output.output_path)
+        del chi_d, chi_m
 
     if comm.rank == 0:
         g2_dens_full.save(name="g2_dens_loc", output_dir=config.output.output_path)
@@ -297,19 +298,6 @@ def write_smom(obj_full: SelfEnergy, obj_ineq: SelfEnergy, sl: slice):
         logger.info("Plotted gamma (magn).")
         del gamma_magn_plot, gamma_m_full
 
-        g_dmft_full._ek = ek
-        plotting.chi_checks(
-            [chi_d_full.mat],
-            [chi_m_full.mat],
-            config.sys.beta,
-            ["Loc-tilde"],
-            g_dmft_full.e_kin,
-            name="loc",
-            output_dir=config.output.plotting_path,
-        )
-        del chi_d, chi_m
-        logger.info("Plotted checks of the susceptibility.")
-
         sigma_list = []
         sigma_names = []
         for i, j in it.product(range(config.sys.n_bands), repeat=2):

dgamore/eliashberg_solver.py

@@ -312,12 +312,13 @@ def create_local_ud_diagrams_pp_w0(g_dmft: GreensFunction) -> Tuple[LocalFourPoi
     gchi0_loc_pp_w0 = (
         BubbleGenerator.create_generalized_chi0_pp_w0(g_dmft, gchi_ud_loc_pp_w0.niv, config.sys.beta)
         .extend_vn_to_diagonal()
+        .permute_orbitals("abcd->acbd")
         .flip_frequency_axis(-1, False)
     )
 
     gamma_ud_loc_pp_w0 = config.sys.beta**2 * (
         (gchi_ud_loc_pp_w0 - gchi0_loc_pp_w0).invert() + gchi0_loc_pp_w0.invert()
-    )
+    ).permute_orbitals("abcd->acbd")
 
     gamma_ud_loc_pp_w0.save(output_dir=config.output.eliashberg_path, name="gamma_ud_loc_pp_w0")
 

dgamore/greens_function.py

@@ -105,16 +105,6 @@ def __init__(
             # config.sys.n, config.sys.occ = self._get_fill()
             config.sys.n, config.sys.occ, config.sys.occ_k = self.get_fill_nonlocal()
 
-    @property
-    def e_kin(self):
-        """
-        Returns the kinetic energy of the system, see Eq (22) in G. Rohringer & A. Toschi
-        PHYSICAL REVIEW B 94, 125144 (2016).
-        """
-        ekin = 2 / config.sys.beta * np.sum(np.mean(self._ek[..., None] * self.mat, axis=(0, 1, 2)))
-        assert np.abs(ekin.imag) < 1e-8, "Kinetic energy must be real."
-        return ekin.real
-
     @property
     def ek(self) -> np.ndarray:
         """
@@ -198,7 +188,7 @@ def map_to_full_bz(self, k_grid: KGrid, nq: tuple = None):
 
     def transpose_orbitals(self):
         r"""
-        Transposes the orbitals of the Green's function object like :math:`G_{ab}^k -> G_{ba}^k.
+        Transposes the orbitals of the Green's function object like :math:`G_{ab}^k -> G_{ba}^k`.
         """
         return self.permute_orbitals("ab->ba")
 
@@ -239,6 +229,57 @@ def get_fill_nonlocal(self) -> tuple[float, np.ndarray, np.ndarray]:
         n_el = 2.0 * np.trace(occ_mean).real
         return n_el, occ_mean, occ_k
 
+    def get_ekin(self) -> float:
+        r"""
+        Returns the kinetic energy calculated from the Green's function and the k-dependent occupation.
+        :math:`E_{kin} = \sum_{\sigma\vec{k}ab} \varepsilon(\vec{k})_{ab} n(\vec{k})_{ba}`.
+        """
+        return 2 * np.sum(self._ek * config.sys.occ_k.swapaxes(-1, -2)).real / config.lattice.k_grid.nk_tot
+
+    def get_epot(self, niv_asympt: int = 50000) -> float:
+        r"""
+        Galitskii-Migdal potential energy, moment-corrected.
+
+        E_pot = sum_k Tr[Sigma_inf rho_k]                                  # Hartree (exact, conv. factor)
+              + 1/beta sum_{k,v} Tr[(Sigma - Sigma_inf) G]                 # in-box correlation part
+              + 1/beta [ sum_big - sum_box ] Tr[(Sigma_mod - Sigma_inf) G_mod]   # analytic 1/v^2 tail
+
+        with Sigma_mod - Sigma_inf = -smom1/(iv) and G_mod = [iv + mu - e_k - Sigma_inf]^-1.
+        The model subtraction cancels the 1/v^2 tail of the correlation sum (remainder ~1/v^4),
+        while sum_big supplies the part beyond the stored box.
+        """
+        smom0, smom1 = self._sigma.smom  # Sigma_inf, first tail coeff; both [o1, o2]
+
+        # 1) Hartree: physical (tail-corrected) occupation, convergence factor exact.
+        e_hartree = np.sum(smom0[None, None, None] * config.sys.occ_k.swapaxes(-1, -2)).real
+
+        # 2) In-box correlation part: Tr[(Sigma - Sigma_inf) G], Sigma_inf already counted above.
+        dsigma = self._sigma.decompress_q_dimension().mat - smom0[..., None]
+        g_ba = self.decompress_q_dimension().transpose_orbitals().mat
+        e_corr = (dsigma * g_ba).sum().real / config.sys.beta
+
+        # 3) Analytic 1/v^2 model tail: replace the truncated box value by the large-box one.
+        e_tail = self._model_epot(smom0, smom1, niv_asympt, config.sys.beta) - self._model_epot(
+            smom0, smom1, self.niv, config.sys.beta
+        )
+
+        return (e_hartree + e_corr + e_tail) / config.lattice.k_grid.nk_tot
+
+    def _model_epot(self, smom0, smom1, niv, beta):
+        """
+        Analytic 1/v^2 model tail
+        """
+        h = (self._ek + smom0[None, None, None]).reshape(self.nq_tot, self.n_bands, self.n_bands)
+        lam, u = np.linalg.eig(h)  # once per k
+        u_inv = np.linalg.inv(u)
+        smom1_rot = u_inv @ smom1 @ u  # rotate tail coeff into eigenbasis
+
+        iv = 1j * MFHelper.vn(niv, beta)
+        g_diag = 1.0 / (iv[None, :] + config.sys.mu - lam[:, :, None])  # [k, band, v]
+        # Tr[(-smom1/iv) G_mod] = -sum_i (smom1_rot)_ii * g_diag_i / iv
+        integrand = -np.einsum("kii,kiv->kv", smom1_rot, g_diag) / iv[None, :]
+        return integrand.sum().real / beta
+
     def _get_gfull_mat(self) -> np.ndarray:
         """
         Returns the full Green's function.

dgamore/n_point_base.py

@@ -440,7 +440,7 @@ def find_q(self, q: tuple[int, int, int] = (0, 0, 0)):
     def filter_q_index(self, index: int = 0):
         r"""
         Filters the object to the given index of the momentum dimension and returns a copy. Acts like a filter.
-        Makes it possible to use e.g. only the first component of a non-local object.
+        Makes it possible to use e.g. only the n-th q-component of a non-local object.
         """
         if not self.has_compressed_q_dimension:
             self.compress_q_dimension()
@@ -451,6 +451,19 @@ def filter_q_index(self, index: int = 0):
         copy._nq = (1, 1, 1)
         return copy
 
+    def q_mean(self):
+        r"""
+        Averages over the momentum dimension and returns a copy of the object with nq = (1,1,1).
+        """
+        copy = deepcopy(self)
+        if self.has_compressed_q_dimension:
+            copy.mat = np.mean(copy.mat, axis=0)[None, ...]
+        else:
+            copy.mat = np.mean(copy.mat, axis=(0, 1, 2))[None, None, None, ...]
+        copy.update_original_shape()
+        copy._nq = (1, 1, 1)
+        return copy
+
     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.

dgamore/nonlocal_sde.py

@@ -760,7 +760,7 @@ def calculate_self_energy_q(
         config.sys.n, config.sys.occ, config.sys.occ_k = giwk_full.get_fill_nonlocal()
         giwk_full.cut_niv(niv_cut)
 
-        if sigma_old is sigma_dmft:
+        if np.allclose(sigma_old.cut_niv(config.box.niv_core).mat, sigma_dmft.cut_niv(config.box.niv_core).mat):
             giwk_full.save(output_dir=config.output.output_path, name="g_latt_dmft")
     config.sys.n, config.sys.occ, config.sys.occ_k = comm.bcast(
         (config.sys.n, config.sys.occ, config.sys.occ_k), root=0
@@ -888,10 +888,21 @@ def calculate_self_energy_q(
         # calculated in this code and add the smooth dmft self-energy
         sigma_new += delta_sigma
         sigma_new = sigma_new.concatenate_self_energies(sigma_dmft)
+        delta_sigma = sigma_dmft.cut_niv(config.box.niv_core) - sigma_new.q_mean().cut_niv(config.box.niv_core)
 
         logger.info("Applying mixing strategy to the self-energy.")
         sigma_new = apply_mixing_strategy(sigma_new, sigma_old, sigma_dmft, current_iter)
 
+        # Canonical self-consistency residual. This is the relative (L2) residual used for the convergence check
+        sigma_new = sigma_new.compress_q_dimension()
+        sigma_old = sigma_old.compress_q_dimension()
+
+        sigma_new_test = sigma_new.mat[..., sigma_new.niv : sigma_new.niv + config.box.niv_core]
+        sigma_old_test = sigma_old.mat[..., sigma_new.niv : sigma_old.niv + config.box.niv_core]
+        diff = (sigma_new_test - sigma_old_test).ravel()
+        norm_x = np.linalg.norm(np.concatenate([sigma_old_test.real.ravel(), sigma_old_test.imag.ravel()]))
+        relative_residual = np.linalg.norm(np.concatenate([diff.real, diff.imag])) / norm_x
+
         old_mu = mu_history[-1]
         if comm.rank == 0:
             mu_finding_failed = False
@@ -916,9 +927,17 @@ def calculate_self_energy_q(
             # calculate new occupation matrix from new Green's function (outside asympt region it is the DMFT
             # lattice Green's function)
             _, config.sys.occ, config.sys.occ_k = giwk_occ.get_fill_nonlocal()  # n should not change
+
+            ekin = giwk_occ.get_ekin()
+            logger.info(f"Kinetic energy: {ekin:.4f} [t or eV].")
+
+            epot = giwk_occ.get_epot()
+            logger.info(f"Potential energy: {epot:.4f} [t or eV].")
+            logger.info(f"Total energy: {(ekin + epot):.4f} [t or eV].")
         config.sys.occ, config.sys.occ_k = comm.bcast((config.sys.occ, config.sys.occ_k), root=0)
+
         if config.self_consistency.max_iter > 1:
-            logger.info("Updated occupation matrix with new Green's function.")
+            logger.info("Updated occupation matrix from new Green's function.")
 
         if comm.rank == 0:
             sigma_new.save(name=f"sigma_dga_iteration_{current_iter}", output_dir=config.output.output_path)
@@ -935,15 +954,11 @@ def calculate_self_energy_q(
 
         logger.info("Checking self-consistency convergence.")
         if comm.rank == 0 and current_iter > starting_iter + 1:
-            niv_start = sigma_new.niv
-            niv_end = niv_start + int(np.ceil(config.box.niv_core / 5))
-
-            sigma_converged = np.allclose(
-                sigma_old.compress_q_dimension()[..., niv_start:niv_end],
-                sigma_new.compress_q_dimension()[..., niv_start:niv_end],
-                atol=config.self_consistency.epsilon,
+            sigma_converged = abs(relative_residual) < config.self_consistency.epsilon
+            logger.info(
+                f"Self-energy convergence: {sigma_converged} "
+                f"(relative residual={relative_residual:.3e}, epsilon={config.self_consistency.epsilon:.3e})."
             )
-            logger.info(f"Self-energy convergence: {sigma_converged}.")
 
             mu_converged = (
                 abs(mu_history[-1] - mu_history[-2]) < np.pi / (10 * config.sys.beta)
@@ -1029,8 +1044,6 @@ def get_result(idx: int):
         f_i[:n_total] = iter_diff.real
         f_i[n_total:] = iter_diff.imag
         norm_f = np.linalg.norm(f_i)
-        norm_x = np.linalg.norm(get_proposal(-1))
-        rel_res = norm_f / norm_x if norm_x > 0 else np.inf
 
         # Solve min||F @ c - f_i|| via truncated-SVD pseudoinverse (drops collinear directions)
         u, s, vh = np.linalg.svd(f_matrix, full_matrices=False)
@@ -1052,9 +1065,7 @@ def get_result(idx: int):
         # 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 (m={n_hist}, alpha={alpha:.3f}, norm_f={norm_f:.3e}, rel_res={rel_res:.3e})."
-        )
+        logger.info(f"Pulay mixing applied (m={n_hist}, alpha={alpha:.3f}, norm_f={norm_f:.3e}).")
 
         return sigma_new
     if config.self_consistency.mixing_strategy.lower() == "anderson" and current_iter > n_hist:
@@ -1081,9 +1092,6 @@ def get_result(idx: int):
         f_curr = flat(last_results[-1]) - flat(last_proposals[-1])
         f_vec = np.concatenate([f_curr.real, f_curr.imag])
         norm_f = np.linalg.norm(f_vec)
-        x_curr = flat(last_proposals[-1])
-        norm_x = np.linalg.norm(np.concatenate([x_curr.real, x_curr.imag]))
-        rel_res = norm_f / norm_x if norm_x > 0 else np.inf
 
         # Build dX and dF matrices (n_hist columns)
         # dX[:,i] = x_{n-i} - x_{n-i-1}  (proposal differences)
@@ -1138,14 +1146,10 @@ def get_result(idx: int):
 
         sigma_new.mat[..., sl] = candidate.reshape(shape)
 
-        logger.info(
-            f"Anderson acceleration applied (m={n_hist}, alpha={alpha:.3f}, norm_f={norm_f:.3e}, rel_res={rel_res:.3e})."
-        )
+        logger.info(f"Anderson acceleration applied (m={n_hist}, alpha={alpha:.3f}, norm_f={norm_f:.3e}).")
 
         return sigma_new
 
-    sigma_new = config.self_consistency.mixing * sigma_new + (1 - config.self_consistency.mixing) * sigma_old
-    logger.info(
-        f"Sigma linearly mixed with previous iteration using a mixing parameter of {config.self_consistency.mixing}."
-    )
+    sigma_new = alpha * sigma_new + (1 - alpha) * sigma_old
+    logger.info(f"Sigma linearly mixed (m=1, alpha={alpha}).")
     return sigma_new

tests/test_eliashberg_end_to_end.py

@@ -129,5 +129,5 @@ def test_eliashberg_equation_with_local_part(setup, niw_core, niv_core, niv_shel
     lambdas_sing, lambdas_trip, gaps_sing, gaps_trip = eliashberg_solver.solve(
         g_dga, g_dmft, u_loc, v_nonloc, comm_mock
     )
-    assert np.allclose(lambdas_sing, np.array([15.80373132, 14.68344248, 12.59236782, 10.8154743]), atol=1e-4)
+    assert np.allclose(lambdas_sing, np.array([15.80373154, 14.68344217, 12.59137005, 10.81553018]), atol=1e-4)
     assert np.allclose(lambdas_trip, np.array([6.70800083, 6.70799431, 6.45388305, 6.45387905]), atol=1e-4)

tests/test_mixing.py

@@ -14,7 +14,6 @@
 from dgamore.self_energy import SelfEnergy
 from dgamore.nonlocal_sde import apply_mixing_strategy
 
-
 BETA = 10.0
 NB = 1
 NK = (1, 1, 1)