Added option to use shift-invert mode for the Lanczos eigenvalue solver

moldga/config.py

@@ -101,6 +101,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:

moldga/config_parser.py

@@ -187,6 +187,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
 

moldga/dga_config.yaml

@@ -58,6 +58,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:

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 ...")

moldga/eliashberg_solver.py

@@ -242,7 +242,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).
@@ -294,20 +294,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
@@ -316,8 +327,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 = [
@@ -326,8 +338,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
@@ -424,20 +437,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,
+    )

moldga/n_point_base.py

@@ -208,7 +208,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

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
 
 
@@ -52,7 +51,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)
@@ -81,8 +80,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