From f7d7b6f5c36e47d34d639c1d93bc2a54146924e5 Mon Sep 17 00:00:00 2001 From: julpe Date: Tue, 10 Mar 2026 12:05:40 +0100 Subject: [PATCH 1/2] Added option to use shift-invert mode for the Lanczos eigenvalue solver. --- moldga/config.py | 1 + moldga/config_parser.py | 1 + moldga/dga_config.yaml | 2 + moldga/dga_main.py | 33 ++++++++- moldga/eliashberg_solver.py | 106 ++++++++++++++++++++-------- moldga/n_point_base.py | 5 +- tests/test_eliashberg_end_to_end.py | 75 ++++++++++++++++++-- 7 files changed, 186 insertions(+), 37 deletions(-) diff --git a/moldga/config.py b/moldga/config.py index b791bf1..5cec87a 100644 --- a/moldga/config.py +++ b/moldga/config.py @@ -99,6 +99,7 @@ def __init__(self): self.symmetry: str = "random" self.subfolder_name: str = "Eliashberg" self.include_local_part: bool = True + self.use_shift_invert_mode: bool = False class LambdaCorrectionConfig: diff --git a/moldga/config_parser.py b/moldga/config_parser.py index 3e11715..7f77ec7 100644 --- a/moldga/config_parser.py +++ b/moldga/config_parser.py @@ -183,6 +183,7 @@ def _build_eliashberg_config(self, conf_file) -> EliashbergConfig: conf.symmetry = self._try_parse(section, "symmetry", "random") conf.subfolder_name = self._try_parse(section, "subfolder_name", "Eliashberg") conf.include_local_part = self._try_parse(section, "include_local_part", True) + conf.use_shift_invert_mode = self._try_parse(section, "use_shift_invert_mode", False) return conf diff --git a/moldga/dga_config.yaml b/moldga/dga_config.yaml index 3e72a93..70590e8 100644 --- a/moldga/dga_config.yaml +++ b/moldga/dga_config.yaml @@ -53,6 +53,8 @@ eliashberg: symmetry: "random" # Symmetry of the gap function that will be used as a starting point for the power iteration. Available: p-wave-x, p-wave-y, d-wave, random include_local_part: False # If True, the local part of the pairing vertex is included in the Eliashberg calculation. This is only relevant if one might expect s-wave symmetry. # setting this to True will reduce the frequency box of the pairing vertex slightly due to an additional frequency shift that is necessary. + use_shift_invert_mode: False # If True, the code will use shift-invert mode for the power iteration additionally to the regular Lanczos, + # which picks out the eigenvalues closest to 1. These eigenvalues are then saved to a different file. subfolder_name: "Eliashberg" # name for the subfolder where certain outputs are generated self_energy_interpolation: diff --git a/moldga/dga_main.py b/moldga/dga_main.py index 725e5a7..1725151 100644 --- a/moldga/dga_main.py +++ b/moldga/dga_main.py @@ -244,9 +244,16 @@ def execute_dga_routine(): if not np.allclose(config.lattice.q_grid.nk, config.lattice.k_grid.nk): raise ValueError("Eliashberg equation can only be solved when nq = nk.") logger.info("Starting with Eliashberg equation.") - lambdas_sing, lambdas_trip, gaps_sing, gaps_trip = eliashberg_solver.solve( - giwk_dga, g_loc, u_loc, v_nonloc, comm - ) + ( + lambdas_sing, + lambdas_trip, + gaps_sing, + gaps_trip, + lambdas_sing_si, + lambdas_trip_si, + gaps_sing_si, + gaps_trip_si, + ) = eliashberg_solver.solve(giwk_dga, g_loc, u_loc, v_nonloc, comm) if config.output.save_quantities and comm.rank == 0: np.savetxt( @@ -256,9 +263,21 @@ def execute_dga_routine(): fmt="%.9f", ) + if lambdas_sing_si is not None and lambdas_trip_si is not None: + np.savetxt( + os.path.join(config.output.eliashberg_path, "eigenvalues_si.txt"), + [lambdas_sing_si.real, lambdas_trip_si.real], + delimiter=",", + fmt="%.9f", + ) + for i in range(len(gaps_sing)): gaps_sing[i].save(name=f"gap_sing_{i+1}", output_dir=config.output.eliashberg_path) gaps_trip[i].save(name=f"gap_trip_{i+1}", output_dir=config.output.eliashberg_path) + + if gaps_sing_si is not None and gaps_trip_si is not None: + gaps_sing_si[i].save(name=f"gap_sing_si_{i+1}", output_dir=config.output.eliashberg_path) + gaps_trip_si[i].save(name=f"gap_trip_si_{i+1}", output_dir=config.output.eliashberg_path) logger.info("Saved singlet and triplet gap functions to files.") if config.output.do_plotting and comm.rank == 0: @@ -270,6 +289,14 @@ def execute_dga_routine(): plotting.plot_gap_function( gaps_trip[i], kx, ky, name=f"gap_trip_{i+1}", output_dir=config.output.eliashberg_path ) + + if gaps_sing_si is not None and gaps_trip_si is not None: + plotting.plot_gap_function( + gaps_sing_si[i], kx, ky, name=f"gap_sing_si_{i+1}", output_dir=config.output.eliashberg_path + ) + plotting.plot_gap_function( + gaps_trip_si[i], kx, ky, name=f"gap_trip_si_{i+1}", output_dir=config.output.eliashberg_path + ) logger.info("Plotted singlet and triplet gap functions.") logger.info("Exiting ...") diff --git a/moldga/eliashberg_solver.py b/moldga/eliashberg_solver.py index fd4bf70..b78b00c 100644 --- a/moldga/eliashberg_solver.py +++ b/moldga/eliashberg_solver.py @@ -231,7 +231,7 @@ def get_initial_gap_function(shape: tuple, channel: SpinChannel) -> np.ndarray: # --- Eliashberg eigensolver (Lanczos / ARPACK) --- -def solve_eliashberg_lanczos(gamma_q_r_pp: FourPoint, gchi0_q0_pp: FourPoint): +def solve_eliashberg_lanczos(gamma_q_r_pp: FourPoint, gchi0_q0_pp: FourPoint, use_shift_invert_mode: bool = False): """ Solves the Eliashberg equation for the superconducting eigenvalue and gap function using ARPACK. Returns (lambdas, gaps). @@ -243,9 +243,11 @@ def solve_eliashberg_lanczos(gamma_q_r_pp: FourPoint, gchi0_q0_pp: FourPoint): allowed_ranks=(0, 1), ) - gamma_q_r_pp = gamma_q_r_pp.map_to_full_bz( - config.lattice.q_grid.irrk_inv, config.lattice.q_grid.nk - ).decompress_q_dimension() + gamma_q_r_pp = ( + gamma_q_r_pp.compress_q_dimension() + .map_to_full_bz(config.lattice.q_grid.irrk_inv, config.lattice.q_grid.nk) + .decompress_q_dimension() + ) logger.log_memory_usage(f"Gamma_pp_{gamma_q_r_pp.channel.value}", gamma_q_r_pp, 1, allowed_ranks=(0, 1)) sign = 1 if gamma_q_r_pp.channel == SpinChannel.SING else -1 @@ -285,20 +287,31 @@ def mv(gap: np.ndarray): n_eig = config.eliashberg.n_eig eig_label = "" if n_eig > 1 else f" {n_eig}" plural = "" if n_eig == 1 else "s" - logger.info( - f"Starting Lanczos method to retrieve largest{eig_label} eigenvalue{plural} and eigenvector{plural} " - f"for the {gamma_q_r_pp.channel.value}let channel.", - allowed_ranks=(0, 1), - ) - lambdas, gaps = sp.sparse.linalg.eigsh( - mat, k=n_eig, tol=config.eliashberg.epsilon, v0=gap0, which="LA", maxiter=10000 - ) + if not use_shift_invert_mode: + logger.info( + f"Starting Lanczos method to retrieve largest{eig_label} eigenvalue{plural} and eigenvector{plural} " + f"for the {gamma_q_r_pp.channel.value}let channel.", + allowed_ranks=(0, 1), + ) + + lambdas, gaps = sp.sparse.linalg.eigsh( + mat, k=n_eig, tol=config.eliashberg.epsilon, v0=gap0, which="LA", maxiter=10000 + ) + else: + logger.info( + f"Starting Lanczos method to retrieve {eig_label} eigenvalue{plural} and eigenvector{plural} " + f"for the {gamma_q_r_pp.channel.value}let channel that are closest to one using shift-invert mode.", + allowed_ranks=(0, 2, 3), + ) + + lambdas, gaps = sp.sparse.linalg.eigsh( + mat, k=n_eig, tol=config.eliashberg.epsilon, v0=gap0, sigma=1.0, which="LM", maxiter=10000 + ) logger.info( - f"Finished Lanczos method for the largest{eig_label} eigenvalue{plural} and eigenvector{plural} " - f"for the {gamma_q_r_pp.channel.value}let channel.", - allowed_ranks=(0, 1), + f"Finished Lanczos method for the {gamma_q_r_pp.channel.value}let channel{" in shift-invert mode" if config.eliashberg.use_shift_invert_mode else ""}.", + allowed_ranks=(0, 1, 2, 3), ) order = lambdas.argsort()[::-1] # sort eigenvalues in descending order @@ -307,8 +320,9 @@ def mv(gap: np.ndarray): logger.info( f"Largest{eig_label} eigenvalue{plural} for the {gamma_q_r_pp.channel.value}let " - f"channel {"is" if n_eig == 1 else "are"}: " + ", ".join(f"{lam:.6f}" for lam in lambdas), - allowed_ranks=(0, 1), + f"channel {"using shift-invert mode" if use_shift_invert_mode else ""} " + f"{"is" if n_eig == 1 else "are"}: " + ", ".join(f"{lam:.6f}" for lam in lambdas), + allowed_ranks=(0, 1, 2, 3), ) gaps = [ @@ -317,8 +331,9 @@ def mv(gap: np.ndarray): ] logger.info( - f"Finished solving the Eliashberg equation for the {gamma_q_r_pp.channel.value}let channel.", - allowed_ranks=(0, 1), + f"Finished solving the Eliashberg equation for the {gamma_q_r_pp.channel.value}let " + f"channel{" in shift-invert mode" if use_shift_invert_mode else ""}.", + allowed_ranks=(0, 1, 2, 3), ) return lambdas, gaps @@ -411,20 +426,53 @@ def solve( ) gchi0_q0_pp = mpi_dist_irrk.bcast(gchi0_q0_pp, root=0) + gamma_sing_pp = mpi_dist_irrk.bcast(gamma_sing_pp, root=0) gamma_trip_pp = mpi_dist_irrk.bcast(gamma_trip_pp, root=0) lambdas_sing, lambdas_trip, gaps_sing, gaps_trip = (None,) * 4 - if mpi_dist_irrk.my_rank == 0: - lambdas_sing, gaps_sing = solve_eliashberg_lanczos(gamma_sing_pp, gchi0_q0_pp) - if mpi_dist_irrk.mpi_size == 1 or mpi_dist_irrk.my_rank == 1: - lambdas_trip, gaps_trip = solve_eliashberg_lanczos(gamma_trip_pp, gchi0_q0_pp) + lambdas_sing_si, lambdas_trip_si, gaps_sing_si, gaps_trip_si = (None,) * 4 - mpi_dist_irrk.delete_file() + size = mpi_dist_irrk.mpi_size + + root_sing = 0 % size + root_trip = 1 % size + root_sing_si = 2 % size + root_trip_si = 3 % size + + if mpi_dist_irrk.my_rank == root_sing: + lambdas_sing, gaps_sing = solve_eliashberg_lanczos(gamma_sing_pp, gchi0_q0_pp, False) + + if mpi_dist_irrk.my_rank == root_trip: + lambdas_trip, gaps_trip = solve_eliashberg_lanczos(gamma_trip_pp, gchi0_q0_pp, False) + + if config.eliashberg.use_shift_invert_mode: + if mpi_dist_irrk.my_rank == root_sing_si: + lambdas_sing_si, gaps_sing_si = solve_eliashberg_lanczos(gamma_sing_pp, gchi0_q0_pp, True) - lambdas_sing = mpi_dist_irrk.bcast(lambdas_sing, root=0) - lambdas_trip = mpi_dist_irrk.bcast(lambdas_trip, root=1 if mpi_dist_irrk.mpi_size > 1 else 0) + if mpi_dist_irrk.my_rank == root_trip_si: + lambdas_trip_si, gaps_trip_si = solve_eliashberg_lanczos(gamma_trip_pp, gchi0_q0_pp, True) - gaps_sing = mpi_dist_irrk.bcast(gaps_sing, root=0) - gaps_trip = mpi_dist_irrk.bcast(gaps_trip, root=1 if mpi_dist_irrk.mpi_size > 1 else 0) + lambdas_sing = mpi_dist_irrk.bcast(lambdas_sing, root=root_sing) + gaps_sing = mpi_dist_irrk.bcast(gaps_sing, root=root_sing) - return lambdas_sing, lambdas_trip, gaps_sing, gaps_trip + lambdas_trip = mpi_dist_irrk.bcast(lambdas_trip, root=root_trip) + gaps_trip = mpi_dist_irrk.bcast(gaps_trip, root=root_trip) + + lambdas_sing_si = mpi_dist_irrk.bcast(lambdas_sing_si, root=root_sing_si) + gaps_sing_si = mpi_dist_irrk.bcast(gaps_sing_si, root=root_sing_si) + + lambdas_trip_si = mpi_dist_irrk.bcast(lambdas_trip_si, root=root_trip_si) + gaps_trip_si = mpi_dist_irrk.bcast(gaps_trip_si, root=root_trip_si) + + mpi_dist_irrk.delete_file() + + return ( + lambdas_sing, + lambdas_trip, + gaps_sing, + gaps_trip, + lambdas_sing_si, + lambdas_trip_si, + gaps_sing_si, + gaps_trip_si, + ) diff --git a/moldga/n_point_base.py b/moldga/n_point_base.py index e6512a2..94e214b 100644 --- a/moldga/n_point_base.py +++ b/moldga/n_point_base.py @@ -204,7 +204,10 @@ def _malloc_trim(cls): Returns unused heap memory to the OS using glibc malloc_trim. Only available on Linux systems. """ - import os + try: + import os + except ImportError: + return if cls._malloc_trim_available is False: return diff --git a/tests/test_eliashberg_end_to_end.py b/tests/test_eliashberg_end_to_end.py index c90348f..dfecffc 100644 --- a/tests/test_eliashberg_end_to_end.py +++ b/tests/test_eliashberg_end_to_end.py @@ -3,12 +3,11 @@ import numpy as np import pytest +import scipy.sparse.linalg from moldga import config, eliashberg_solver, dga_io from moldga.dga_logger import DgaLogger from moldga.greens_function import GreensFunction -from moldga.local_four_point import LocalFourPoint -from moldga.n_point_base import SpinChannel from tests import conftest @@ -51,7 +50,7 @@ def test_eliashberg_equation_without_local_part(setup, niw_core, niv_core, niv_s g_dga = GreensFunction(np.load(f"{folder}/giwk_dga.npy")) - lambdas_sing, lambdas_trip, gaps_sing, gaps_trip = eliashberg_solver.solve( + 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([3.85828144, 3.70361068, 3.65005429, 3.5992988]), atol=1e-4) @@ -79,8 +78,76 @@ def test_eliashberg_equation_with_local_part(setup, niw_core, niv_core, niv_shel g_dga = GreensFunction(np.load(f"{folder}/giwk_dga.npy")) - lambdas_sing, lambdas_trip, gaps_sing, gaps_trip = eliashberg_solver.solve( + 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([3.7036108, 3.5992989, 3.32485204, 3.32485072]), atol=1e-4) assert np.allclose(lambdas_trip, np.array([2.72114656, 2.72114542, 2.69452022, 2.69451905]), atol=1e-4) + + +@pytest.mark.parametrize( + "niw_core, niv_core, niv_shell, include_local_part, use_shift_invert_mode", + [(20, 20, 10, True, True), (20, 20, 10, True, False), (20, 20, 10, False, True), (20, 20, 10, False, False)], +) +def test_eliashberg_equation_with_shift_invert_mode( + setup, niw_core, niv_core, niv_shell, include_local_part, use_shift_invert_mode +): + folder, comm_mock = setup + + config.box.niw_core = niw_core + config.box.niv_core = niv_core + config.box.niv_shell = niv_shell + + g_dmft, s_dmft, g2_dens, g2_magn = dga_io.load_from_w2dyn_file_and_update_config() + + config.eliashberg.perform_eliashberg = True + config.output.output_path = folder + config.output.eliashberg_path = config.output.output_path + config.eliashberg.include_local_part = include_local_part + config.eliashberg.use_shift_invert_mode = use_shift_invert_mode + + u_loc = config.lattice.hamiltonian.get_local_u() + v_nonloc = config.lattice.hamiltonian.get_vq(config.lattice.q_grid) + + g_dga = GreensFunction(np.load(f"{folder}/giwk_dga.npy")) + + _real_eigsh = scipy.sparse.linalg.eigsh + + def mock_eigsh(mat, k=1, tol=None, v0=None, sigma=None, which="LM", maxiter=None): + if sigma == 1.0 and which.upper() == "LM": + return np.array([1.0, 0.8, 0.9, 0.7]), np.random.rand(4, 4, 1, 2, 2, 20, 4) + else: + try: + return _real_eigsh(mat, k=k, tol=tol, v0=v0, sigma=sigma, which=which, maxiter=maxiter) + except Exception: + raise RuntimeError() + + with patch.object(scipy.sparse.linalg, "eigsh", side_effect=mock_eigsh) as mock: + (l_sing, l_trip, gaps_sing, gaps_trip, l_sing_si, l_trip_si, gaps_sing_si, gaps_trip_si) = ( + eliashberg_solver.solve(g_dga, g_dmft, u_loc, v_nonloc, comm_mock) + ) + + if use_shift_invert_mode: + assert mock.call_count == 4 + assert any( + call.kwargs.get("sigma") == 1.0 and call.kwargs.get("which") == "LM" for call in mock.call_args_list + ) + else: + assert mock.call_count == 2 + + if include_local_part: + assert np.allclose(l_sing, np.array([3.7036108, 3.5992989, 3.32485204, 3.32485072]), atol=1e-4) + assert np.allclose(l_trip, np.array([2.72114656, 2.72114542, 2.69452022, 2.69451905]), atol=1e-4) + else: + assert np.allclose(l_sing, np.array([3.85828144, 3.70361068, 3.65005429, 3.5992988]), atol=1e-4) + assert np.allclose(l_trip, np.array([3.34166718, 2.9909934, 2.72114652, 2.72114537]), atol=1e-4) + + if use_shift_invert_mode: + assert not np.allclose(l_sing_si, l_sing, atol=1e-4) + assert not np.allclose(l_trip_si, l_trip, atol=1e-4) + + if not use_shift_invert_mode: + assert l_sing_si is None + assert l_trip_si is None + assert gaps_sing_si is None + assert gaps_trip_si is None From ca59069e24869de85123633c797410788a9acd27 Mon Sep 17 00:00:00 2001 From: julpe Date: Tue, 10 Mar 2026 12:05:40 +0100 Subject: [PATCH 2/2] Added option to use shift-invert mode for the Lanczos eigenvalue solver. --- moldga/config.py | 1 + moldga/config_parser.py | 1 + moldga/dga_config.yaml | 2 + moldga/dga_main.py | 33 ++++++++- moldga/eliashberg_solver.py | 106 ++++++++++++++++++++-------- moldga/n_point_base.py | 5 +- tests/test_eliashberg_end_to_end.py | 75 ++++++++++++++++++-- 7 files changed, 186 insertions(+), 37 deletions(-) diff --git a/moldga/config.py b/moldga/config.py index b791bf1..5cec87a 100644 --- a/moldga/config.py +++ b/moldga/config.py @@ -99,6 +99,7 @@ def __init__(self): self.symmetry: str = "random" self.subfolder_name: str = "Eliashberg" self.include_local_part: bool = True + self.use_shift_invert_mode: bool = False class LambdaCorrectionConfig: diff --git a/moldga/config_parser.py b/moldga/config_parser.py index 3e11715..7f77ec7 100644 --- a/moldga/config_parser.py +++ b/moldga/config_parser.py @@ -183,6 +183,7 @@ def _build_eliashberg_config(self, conf_file) -> EliashbergConfig: conf.symmetry = self._try_parse(section, "symmetry", "random") conf.subfolder_name = self._try_parse(section, "subfolder_name", "Eliashberg") conf.include_local_part = self._try_parse(section, "include_local_part", True) + conf.use_shift_invert_mode = self._try_parse(section, "use_shift_invert_mode", False) return conf diff --git a/moldga/dga_config.yaml b/moldga/dga_config.yaml index 3e72a93..70590e8 100644 --- a/moldga/dga_config.yaml +++ b/moldga/dga_config.yaml @@ -53,6 +53,8 @@ eliashberg: symmetry: "random" # Symmetry of the gap function that will be used as a starting point for the power iteration. Available: p-wave-x, p-wave-y, d-wave, random include_local_part: False # If True, the local part of the pairing vertex is included in the Eliashberg calculation. This is only relevant if one might expect s-wave symmetry. # setting this to True will reduce the frequency box of the pairing vertex slightly due to an additional frequency shift that is necessary. + use_shift_invert_mode: False # If True, the code will use shift-invert mode for the power iteration additionally to the regular Lanczos, + # which picks out the eigenvalues closest to 1. These eigenvalues are then saved to a different file. subfolder_name: "Eliashberg" # name for the subfolder where certain outputs are generated self_energy_interpolation: diff --git a/moldga/dga_main.py b/moldga/dga_main.py index 725e5a7..1725151 100644 --- a/moldga/dga_main.py +++ b/moldga/dga_main.py @@ -244,9 +244,16 @@ def execute_dga_routine(): if not np.allclose(config.lattice.q_grid.nk, config.lattice.k_grid.nk): raise ValueError("Eliashberg equation can only be solved when nq = nk.") logger.info("Starting with Eliashberg equation.") - lambdas_sing, lambdas_trip, gaps_sing, gaps_trip = eliashberg_solver.solve( - giwk_dga, g_loc, u_loc, v_nonloc, comm - ) + ( + lambdas_sing, + lambdas_trip, + gaps_sing, + gaps_trip, + lambdas_sing_si, + lambdas_trip_si, + gaps_sing_si, + gaps_trip_si, + ) = eliashberg_solver.solve(giwk_dga, g_loc, u_loc, v_nonloc, comm) if config.output.save_quantities and comm.rank == 0: np.savetxt( @@ -256,9 +263,21 @@ def execute_dga_routine(): fmt="%.9f", ) + if lambdas_sing_si is not None and lambdas_trip_si is not None: + np.savetxt( + os.path.join(config.output.eliashberg_path, "eigenvalues_si.txt"), + [lambdas_sing_si.real, lambdas_trip_si.real], + delimiter=",", + fmt="%.9f", + ) + for i in range(len(gaps_sing)): gaps_sing[i].save(name=f"gap_sing_{i+1}", output_dir=config.output.eliashberg_path) gaps_trip[i].save(name=f"gap_trip_{i+1}", output_dir=config.output.eliashberg_path) + + if gaps_sing_si is not None and gaps_trip_si is not None: + gaps_sing_si[i].save(name=f"gap_sing_si_{i+1}", output_dir=config.output.eliashberg_path) + gaps_trip_si[i].save(name=f"gap_trip_si_{i+1}", output_dir=config.output.eliashberg_path) logger.info("Saved singlet and triplet gap functions to files.") if config.output.do_plotting and comm.rank == 0: @@ -270,6 +289,14 @@ def execute_dga_routine(): plotting.plot_gap_function( gaps_trip[i], kx, ky, name=f"gap_trip_{i+1}", output_dir=config.output.eliashberg_path ) + + if gaps_sing_si is not None and gaps_trip_si is not None: + plotting.plot_gap_function( + gaps_sing_si[i], kx, ky, name=f"gap_sing_si_{i+1}", output_dir=config.output.eliashberg_path + ) + plotting.plot_gap_function( + gaps_trip_si[i], kx, ky, name=f"gap_trip_si_{i+1}", output_dir=config.output.eliashberg_path + ) logger.info("Plotted singlet and triplet gap functions.") logger.info("Exiting ...") diff --git a/moldga/eliashberg_solver.py b/moldga/eliashberg_solver.py index fd4bf70..b78b00c 100644 --- a/moldga/eliashberg_solver.py +++ b/moldga/eliashberg_solver.py @@ -231,7 +231,7 @@ def get_initial_gap_function(shape: tuple, channel: SpinChannel) -> np.ndarray: # --- Eliashberg eigensolver (Lanczos / ARPACK) --- -def solve_eliashberg_lanczos(gamma_q_r_pp: FourPoint, gchi0_q0_pp: FourPoint): +def solve_eliashberg_lanczos(gamma_q_r_pp: FourPoint, gchi0_q0_pp: FourPoint, use_shift_invert_mode: bool = False): """ Solves the Eliashberg equation for the superconducting eigenvalue and gap function using ARPACK. Returns (lambdas, gaps). @@ -243,9 +243,11 @@ def solve_eliashberg_lanczos(gamma_q_r_pp: FourPoint, gchi0_q0_pp: FourPoint): allowed_ranks=(0, 1), ) - gamma_q_r_pp = gamma_q_r_pp.map_to_full_bz( - config.lattice.q_grid.irrk_inv, config.lattice.q_grid.nk - ).decompress_q_dimension() + gamma_q_r_pp = ( + gamma_q_r_pp.compress_q_dimension() + .map_to_full_bz(config.lattice.q_grid.irrk_inv, config.lattice.q_grid.nk) + .decompress_q_dimension() + ) logger.log_memory_usage(f"Gamma_pp_{gamma_q_r_pp.channel.value}", gamma_q_r_pp, 1, allowed_ranks=(0, 1)) sign = 1 if gamma_q_r_pp.channel == SpinChannel.SING else -1 @@ -285,20 +287,31 @@ def mv(gap: np.ndarray): n_eig = config.eliashberg.n_eig eig_label = "" if n_eig > 1 else f" {n_eig}" plural = "" if n_eig == 1 else "s" - logger.info( - f"Starting Lanczos method to retrieve largest{eig_label} eigenvalue{plural} and eigenvector{plural} " - f"for the {gamma_q_r_pp.channel.value}let channel.", - allowed_ranks=(0, 1), - ) - lambdas, gaps = sp.sparse.linalg.eigsh( - mat, k=n_eig, tol=config.eliashberg.epsilon, v0=gap0, which="LA", maxiter=10000 - ) + if not use_shift_invert_mode: + logger.info( + f"Starting Lanczos method to retrieve largest{eig_label} eigenvalue{plural} and eigenvector{plural} " + f"for the {gamma_q_r_pp.channel.value}let channel.", + allowed_ranks=(0, 1), + ) + + lambdas, gaps = sp.sparse.linalg.eigsh( + mat, k=n_eig, tol=config.eliashberg.epsilon, v0=gap0, which="LA", maxiter=10000 + ) + else: + logger.info( + f"Starting Lanczos method to retrieve {eig_label} eigenvalue{plural} and eigenvector{plural} " + f"for the {gamma_q_r_pp.channel.value}let channel that are closest to one using shift-invert mode.", + allowed_ranks=(0, 2, 3), + ) + + lambdas, gaps = sp.sparse.linalg.eigsh( + mat, k=n_eig, tol=config.eliashberg.epsilon, v0=gap0, sigma=1.0, which="LM", maxiter=10000 + ) logger.info( - f"Finished Lanczos method for the largest{eig_label} eigenvalue{plural} and eigenvector{plural} " - f"for the {gamma_q_r_pp.channel.value}let channel.", - allowed_ranks=(0, 1), + f"Finished Lanczos method for the {gamma_q_r_pp.channel.value}let channel{" in shift-invert mode" if config.eliashberg.use_shift_invert_mode else ""}.", + allowed_ranks=(0, 1, 2, 3), ) order = lambdas.argsort()[::-1] # sort eigenvalues in descending order @@ -307,8 +320,9 @@ def mv(gap: np.ndarray): logger.info( f"Largest{eig_label} eigenvalue{plural} for the {gamma_q_r_pp.channel.value}let " - f"channel {"is" if n_eig == 1 else "are"}: " + ", ".join(f"{lam:.6f}" for lam in lambdas), - allowed_ranks=(0, 1), + f"channel {"using shift-invert mode" if use_shift_invert_mode else ""} " + f"{"is" if n_eig == 1 else "are"}: " + ", ".join(f"{lam:.6f}" for lam in lambdas), + allowed_ranks=(0, 1, 2, 3), ) gaps = [ @@ -317,8 +331,9 @@ def mv(gap: np.ndarray): ] logger.info( - f"Finished solving the Eliashberg equation for the {gamma_q_r_pp.channel.value}let channel.", - allowed_ranks=(0, 1), + f"Finished solving the Eliashberg equation for the {gamma_q_r_pp.channel.value}let " + f"channel{" in shift-invert mode" if use_shift_invert_mode else ""}.", + allowed_ranks=(0, 1, 2, 3), ) return lambdas, gaps @@ -411,20 +426,53 @@ def solve( ) gchi0_q0_pp = mpi_dist_irrk.bcast(gchi0_q0_pp, root=0) + gamma_sing_pp = mpi_dist_irrk.bcast(gamma_sing_pp, root=0) gamma_trip_pp = mpi_dist_irrk.bcast(gamma_trip_pp, root=0) lambdas_sing, lambdas_trip, gaps_sing, gaps_trip = (None,) * 4 - if mpi_dist_irrk.my_rank == 0: - lambdas_sing, gaps_sing = solve_eliashberg_lanczos(gamma_sing_pp, gchi0_q0_pp) - if mpi_dist_irrk.mpi_size == 1 or mpi_dist_irrk.my_rank == 1: - lambdas_trip, gaps_trip = solve_eliashberg_lanczos(gamma_trip_pp, gchi0_q0_pp) + lambdas_sing_si, lambdas_trip_si, gaps_sing_si, gaps_trip_si = (None,) * 4 - mpi_dist_irrk.delete_file() + size = mpi_dist_irrk.mpi_size + + root_sing = 0 % size + root_trip = 1 % size + root_sing_si = 2 % size + root_trip_si = 3 % size + + if mpi_dist_irrk.my_rank == root_sing: + lambdas_sing, gaps_sing = solve_eliashberg_lanczos(gamma_sing_pp, gchi0_q0_pp, False) + + if mpi_dist_irrk.my_rank == root_trip: + lambdas_trip, gaps_trip = solve_eliashberg_lanczos(gamma_trip_pp, gchi0_q0_pp, False) + + if config.eliashberg.use_shift_invert_mode: + if mpi_dist_irrk.my_rank == root_sing_si: + lambdas_sing_si, gaps_sing_si = solve_eliashberg_lanczos(gamma_sing_pp, gchi0_q0_pp, True) - lambdas_sing = mpi_dist_irrk.bcast(lambdas_sing, root=0) - lambdas_trip = mpi_dist_irrk.bcast(lambdas_trip, root=1 if mpi_dist_irrk.mpi_size > 1 else 0) + if mpi_dist_irrk.my_rank == root_trip_si: + lambdas_trip_si, gaps_trip_si = solve_eliashberg_lanczos(gamma_trip_pp, gchi0_q0_pp, True) - gaps_sing = mpi_dist_irrk.bcast(gaps_sing, root=0) - gaps_trip = mpi_dist_irrk.bcast(gaps_trip, root=1 if mpi_dist_irrk.mpi_size > 1 else 0) + lambdas_sing = mpi_dist_irrk.bcast(lambdas_sing, root=root_sing) + gaps_sing = mpi_dist_irrk.bcast(gaps_sing, root=root_sing) - return lambdas_sing, lambdas_trip, gaps_sing, gaps_trip + lambdas_trip = mpi_dist_irrk.bcast(lambdas_trip, root=root_trip) + gaps_trip = mpi_dist_irrk.bcast(gaps_trip, root=root_trip) + + lambdas_sing_si = mpi_dist_irrk.bcast(lambdas_sing_si, root=root_sing_si) + gaps_sing_si = mpi_dist_irrk.bcast(gaps_sing_si, root=root_sing_si) + + lambdas_trip_si = mpi_dist_irrk.bcast(lambdas_trip_si, root=root_trip_si) + gaps_trip_si = mpi_dist_irrk.bcast(gaps_trip_si, root=root_trip_si) + + mpi_dist_irrk.delete_file() + + return ( + lambdas_sing, + lambdas_trip, + gaps_sing, + gaps_trip, + lambdas_sing_si, + lambdas_trip_si, + gaps_sing_si, + gaps_trip_si, + ) diff --git a/moldga/n_point_base.py b/moldga/n_point_base.py index e6512a2..94e214b 100644 --- a/moldga/n_point_base.py +++ b/moldga/n_point_base.py @@ -204,7 +204,10 @@ def _malloc_trim(cls): Returns unused heap memory to the OS using glibc malloc_trim. Only available on Linux systems. """ - import os + try: + import os + except ImportError: + return if cls._malloc_trim_available is False: return diff --git a/tests/test_eliashberg_end_to_end.py b/tests/test_eliashberg_end_to_end.py index c90348f..dfecffc 100644 --- a/tests/test_eliashberg_end_to_end.py +++ b/tests/test_eliashberg_end_to_end.py @@ -3,12 +3,11 @@ import numpy as np import pytest +import scipy.sparse.linalg from moldga import config, eliashberg_solver, dga_io from moldga.dga_logger import DgaLogger from moldga.greens_function import GreensFunction -from moldga.local_four_point import LocalFourPoint -from moldga.n_point_base import SpinChannel from tests import conftest @@ -51,7 +50,7 @@ def test_eliashberg_equation_without_local_part(setup, niw_core, niv_core, niv_s g_dga = GreensFunction(np.load(f"{folder}/giwk_dga.npy")) - lambdas_sing, lambdas_trip, gaps_sing, gaps_trip = eliashberg_solver.solve( + 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([3.85828144, 3.70361068, 3.65005429, 3.5992988]), atol=1e-4) @@ -79,8 +78,76 @@ def test_eliashberg_equation_with_local_part(setup, niw_core, niv_core, niv_shel g_dga = GreensFunction(np.load(f"{folder}/giwk_dga.npy")) - lambdas_sing, lambdas_trip, gaps_sing, gaps_trip = eliashberg_solver.solve( + 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([3.7036108, 3.5992989, 3.32485204, 3.32485072]), atol=1e-4) assert np.allclose(lambdas_trip, np.array([2.72114656, 2.72114542, 2.69452022, 2.69451905]), atol=1e-4) + + +@pytest.mark.parametrize( + "niw_core, niv_core, niv_shell, include_local_part, use_shift_invert_mode", + [(20, 20, 10, True, True), (20, 20, 10, True, False), (20, 20, 10, False, True), (20, 20, 10, False, False)], +) +def test_eliashberg_equation_with_shift_invert_mode( + setup, niw_core, niv_core, niv_shell, include_local_part, use_shift_invert_mode +): + folder, comm_mock = setup + + config.box.niw_core = niw_core + config.box.niv_core = niv_core + config.box.niv_shell = niv_shell + + g_dmft, s_dmft, g2_dens, g2_magn = dga_io.load_from_w2dyn_file_and_update_config() + + config.eliashberg.perform_eliashberg = True + config.output.output_path = folder + config.output.eliashberg_path = config.output.output_path + config.eliashberg.include_local_part = include_local_part + config.eliashberg.use_shift_invert_mode = use_shift_invert_mode + + u_loc = config.lattice.hamiltonian.get_local_u() + v_nonloc = config.lattice.hamiltonian.get_vq(config.lattice.q_grid) + + g_dga = GreensFunction(np.load(f"{folder}/giwk_dga.npy")) + + _real_eigsh = scipy.sparse.linalg.eigsh + + def mock_eigsh(mat, k=1, tol=None, v0=None, sigma=None, which="LM", maxiter=None): + if sigma == 1.0 and which.upper() == "LM": + return np.array([1.0, 0.8, 0.9, 0.7]), np.random.rand(4, 4, 1, 2, 2, 20, 4) + else: + try: + return _real_eigsh(mat, k=k, tol=tol, v0=v0, sigma=sigma, which=which, maxiter=maxiter) + except Exception: + raise RuntimeError() + + with patch.object(scipy.sparse.linalg, "eigsh", side_effect=mock_eigsh) as mock: + (l_sing, l_trip, gaps_sing, gaps_trip, l_sing_si, l_trip_si, gaps_sing_si, gaps_trip_si) = ( + eliashberg_solver.solve(g_dga, g_dmft, u_loc, v_nonloc, comm_mock) + ) + + if use_shift_invert_mode: + assert mock.call_count == 4 + assert any( + call.kwargs.get("sigma") == 1.0 and call.kwargs.get("which") == "LM" for call in mock.call_args_list + ) + else: + assert mock.call_count == 2 + + if include_local_part: + assert np.allclose(l_sing, np.array([3.7036108, 3.5992989, 3.32485204, 3.32485072]), atol=1e-4) + assert np.allclose(l_trip, np.array([2.72114656, 2.72114542, 2.69452022, 2.69451905]), atol=1e-4) + else: + assert np.allclose(l_sing, np.array([3.85828144, 3.70361068, 3.65005429, 3.5992988]), atol=1e-4) + assert np.allclose(l_trip, np.array([3.34166718, 2.9909934, 2.72114652, 2.72114537]), atol=1e-4) + + if use_shift_invert_mode: + assert not np.allclose(l_sing_si, l_sing, atol=1e-4) + assert not np.allclose(l_trip_si, l_trip, atol=1e-4) + + if not use_shift_invert_mode: + assert l_sing_si is None + assert l_trip_si is None + assert gaps_sing_si is None + assert gaps_trip_si is None