From acd71b31b2cebe767c5ab5590c089996751bc957 Mon Sep 17 00:00:00 2001 From: julpe Date: Mon, 8 Jun 2026 10:43:28 +0200 Subject: [PATCH] Optimized Pulay mixing and ordering of quantities to calculate to stabilize convergence. --- dgamore/greens_function.py | 6 +++++- dgamore/nonlocal_sde.py | 37 +++++++++++++++++++++++++++---------- 2 files changed, 32 insertions(+), 11 deletions(-) diff --git a/dgamore/greens_function.py b/dgamore/greens_function.py index 04f4c236..003aadbb 100644 --- a/dgamore/greens_function.py +++ b/dgamore/greens_function.py @@ -39,7 +39,11 @@ def get_total_fill(mu: float, ek: np.ndarray, sigma_mat: np.ndarray, beta: float g_loc_mat = np.mean(g_full_mat, axis=0) eigenvals, eigenvecs = np.linalg.eig(beta * (hloc.real + smom0 - mu_bands)) - rho_diag = np.where(eigenvals > 0, np.exp(-eigenvals) / (1 + np.exp(-eigenvals)), 1 / (1 + np.exp(eigenvals))) + + rho_diag = np.empty_like(eigenvals) + mask = eigenvals > 0 + rho_diag[mask] = np.exp(-eigenvals[mask]) / (1 + np.exp(-eigenvals[mask])) + rho_diag[~mask] = 1 / (1 + np.exp(eigenvals[~mask])) rho_diag = np.einsum("...i,ij->...ij", rho_diag, np.eye(n_bands)) rho_loc = eigenvecs @ rho_diag @ np.linalg.inv(eigenvecs) diff --git a/dgamore/nonlocal_sde.py b/dgamore/nonlocal_sde.py index 93e6df22..e8ba1a6f 100644 --- a/dgamore/nonlocal_sde.py +++ b/dgamore/nonlocal_sde.py @@ -889,6 +889,9 @@ def calculate_self_energy_q( sigma_new += delta_sigma sigma_new = sigma_new.concatenate_self_energies(sigma_dmft) + logger.info("Applying mixing strategy to the self-energy.") + sigma_new = apply_mixing_strategy(sigma_new, sigma_old, sigma_dmft, current_iter) + old_mu = mu_history[-1] if comm.rank == 0: mu_finding_failed = False @@ -917,9 +920,6 @@ def calculate_self_energy_q( if config.self_consistency.max_iter > 1: logger.info("Updated occupation matrix with new Green's function.") - logger.info("Applying mixing strategy to the self-energy.") - sigma_new = apply_mixing_strategy(sigma_new, sigma_old, sigma_dmft, current_iter) - if comm.rank == 0: sigma_new.save(name=f"sigma_dga_iteration_{current_iter}", output_dir=config.output.output_path) logger.info(f"Saved sigma for iteration {current_iter} as numpy array.") @@ -1028,20 +1028,32 @@ def get_result(idx: int): iter_diff = get_result(-1) - get_proposal(-1) f_i[:n_total] = iter_diff.real f_i[n_total:] = iter_diff.imag - - # Solve min||F @ c - f_i|| via least squares (more stable than explicit inverse) - coeffs, _, _, _ = np.linalg.lstsq(f_matrix, f_i, rcond=None) + 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) + cutoff = 1e-5 * (s[0] if s.size else 1.0) + mask = s > cutoff + if not np.any(mask): + logger.warning("Pulay SVD ill-conditioned - falling back to linear mixing.") + return alpha * sigma_new + (1 - alpha) * sigma_old + coeffs = vh[mask].T @ ((u[:, mask].T @ f_i) / s[mask]) # Pulay update: x_{n+1} = x_n + alpha*f_i - (R + alpha*F) @ c update = alpha * f_i - (r_matrix + alpha * f_matrix) @ coeffs + norm_u = np.linalg.norm(update) + if norm_f > 0 and norm_u > 10.0 * norm_f: + update *= 10.0 * norm_f / norm_u + logger.warning(f"Pulay step clamped (norm_u={norm_u:.3e}, norm_f={norm_f:.3e}).") update = update[:n_total] + 1j * update[n_total:] # 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 " - f"a mixing parameter of {config.self_consistency.mixing}." + f"Pulay mixing applied (m={n_hist}, alpha={alpha:.3f}, norm_f={norm_f:.3e}, rel_res={rel_res:.3e})." ) return sigma_new @@ -1069,6 +1081,9 @@ 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) @@ -1102,7 +1117,7 @@ def get_result(idx: int): coeffs = vh[mask].T @ (s_reg * (u[:, mask].T @ f_vec)) except np.linalg.LinAlgError: - logger.warning("Anderson SVD failed — falling back to linear mixing.") + logger.warning("Anderson SVD failed - falling back to linear mixing.") return alpha * sigma_new + (1 - alpha) * sigma_old # Undamped Anderson proposal: x_n + f_n - (dX + dF) @ c @@ -1123,7 +1138,9 @@ 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}).") + logger.info( + f"Anderson acceleration applied (m={n_hist}, alpha={alpha:.3f}, norm_f={norm_f:.3e}, rel_res={rel_res:.3e})." + ) return sigma_new