From 7766c35597550e4d73de992d3c778295a54d9083 Mon Sep 17 00:00:00 2001 From: julpe Date: Fri, 27 Mar 2026 13:28:52 +0100 Subject: [PATCH 1/2] Added support for La3Ni2O7 --- moldga/config.py | 2 + moldga/dga_io.py | 94 +++++++++++++++-- moldga/hamiltonian.py | 73 +++++++++++-- moldga/symmetrize_new.py | 191 +++++++++++++++++++++++++-------- moldga/w2dyn_aux.py | 88 ++++++++-------- tests/test_dga_io.py | 215 ++++++++++++++++++++++++++++++++++++++ tests/test_hamiltonian.py | 123 +++++++++++++++++++++- tests/test_symmetrize.py | 177 +++++++++++++++++++++++++++++-- 8 files changed, 844 insertions(+), 119 deletions(-) create mode 100644 tests/test_dga_io.py diff --git a/moldga/config.py b/moldga/config.py index e90e808..1a0b0d8 100644 --- a/moldga/config.py +++ b/moldga/config.py @@ -148,6 +148,8 @@ def __init__(self): self.mu: float = 0.0 self.n: float = 0.0 self.n_bands: int = 1 + self.nd_bands: int = 1 + self.np_bands: int = 0 self.occ: np.ndarray = np.ndarray(0) self.occ_k: np.ndarray = np.ndarray(0) self.occ_dmft: np.ndarray = np.ndarray(0) diff --git a/moldga/dga_io.py b/moldga/dga_io.py index 6863f9e..7736b6a 100644 --- a/moldga/dga_io.py +++ b/moldga/dga_io.py @@ -31,9 +31,13 @@ def uniquify_path(path: str = None): return path -def load_from_w2dyn_file_and_update_config() -> tuple[GreensFunction, SelfEnergy, LocalFourPoint, LocalFourPoint]: +def load_from_w2dyn_file_and_update_config( + combine_two_atoms_to_one_obj: bool = False, +) -> tuple[GreensFunction, SelfEnergy, LocalFourPoint, LocalFourPoint]: """ Loads data from the w2dyn file and updates the config file. + If combine_atoms_to_one_obj is True, we are doubling the orbital space and putting the data from the + inequivalent atoms into the diagonal blocks after we average over them. """ file = w2dyn_aux.W2dynFile(fname=str(os.path.join(config.dmft.input_path, config.dmft.fname_1p))) @@ -51,22 +55,60 @@ def load_from_w2dyn_file_and_update_config() -> tuple[GreensFunction, SelfEnergy config.lattice.interaction.vpp = file.get_vpp() config.sys.mu = file.get_mu() - config.sys.n_bands = file.get_nd() + file.get_np() + config.sys.nd_bands = file.get_nd() + config.sys.np_bands = file.get_np() + config.sys.n_bands = config.sys.nd_bands + config.sys.np_bands config.sys.n = file.get_totdens() config.sys.occ_dmft = 2 * np.mean(file.get_rho1(), axis=(1, 3)) + if combine_two_atoms_to_one_obj: + config.sys.n_bands *= 2 + config.sys.nd_bands *= 2 + config.sys.np_bands *= 2 + config.sys.n *= 2 + + config.sys.occ_dmft = np.zeros((config.sys.n_bands, config.sys.n_bands)) + rho_1_mean = 0.5 * (np.mean(file.get_rho1(), axis=(1, 3)) + np.mean(file.get_rho1(ineq=2), axis=(1, 3))) + config.sys.occ_dmft[0:2, 0:2] = 2 * rho_1_mean + config.sys.occ_dmft[2:4, 2:4] = 2 * rho_1_mean + if config.sys.n == 0: config.sys.n = 2 * np.sum(config.sys.occ_dmft) file2 = w2dyn_aux.W2dynG4iwFile(fname=str(os.path.join(config.dmft.input_path, config.dmft.fname_2p))) - g2_dens = LocalFourPoint(file2.read_g2_full_multiband(config.sys.n_bands, name="dens"), channel=SpinChannel.DENS) - g2_magn = LocalFourPoint(file2.read_g2_full_multiband(config.sys.n_bands, name="magn"), channel=SpinChannel.MAGN) + g2_dens = LocalFourPoint( + file2.read_g2_full_multiband(file.get_nd() + file.get_np(), name="dens"), channel=SpinChannel.DENS + ) + g2_magn = LocalFourPoint( + file2.read_g2_full_multiband(file.get_nd() + file.get_np(), name="magn"), channel=SpinChannel.MAGN + ) + + if combine_two_atoms_to_one_obj: + g2_dens_2 = LocalFourPoint( + file2.read_g2_full_multiband(file.get_nd() + file.get_np(), ineq=2, name="dens"), channel=SpinChannel.DENS + ) + g2_magn_2 = LocalFourPoint( + file2.read_g2_full_multiband(file.get_nd() + file.get_np(), ineq=2, name="magn"), channel=SpinChannel.MAGN + ) + + def construct_large_g2(g2_1: LocalFourPoint, g2_2: LocalFourPoint) -> LocalFourPoint: + g2_mean = 0.5 * (g2_1.mat + g2_2.mat) + del g2_2 + g2_1.mat = np.zeros((4, 4, 4, 4, *g2_mean.shape[4:])) + g2_1.mat[0:2, 0:2, 0:2, 0:2] = g2_mean + g2_1.mat[2:4, 2:4, 2:4, 2:4] = g2_mean + del g2_mean + return g2_1 + + g2_dens = construct_large_g2(g2_dens, g2_dens_2) + g2_magn = construct_large_g2(g2_magn, g2_magn_2) + file2.close() update_frequency_boxes(g2_dens.niw, g2_dens.niv) def extend_orbital(arr: np.ndarray) -> np.ndarray: - return np.einsum("i...,ij->ij...", arr, np.eye(config.sys.n_bands)) + return np.einsum("i...,ij->ij...", arr, np.eye(arr.shape[0])) giw_spin_mean = np.mean(file.get_giw(), axis=1) # [band,spin,niv] niv_dmft = giw_spin_mean.shape[-1] // 2 @@ -74,6 +116,16 @@ def extend_orbital(arr: np.ndarray) -> np.ndarray: giw_spin_mean = giw_spin_mean[..., niv_dmft - niv_cut : niv_dmft + niv_cut] g_dmft = GreensFunction(extend_orbital(giw_spin_mean)) + if combine_two_atoms_to_one_obj: + giw_spin_mean_2 = np.mean(file.get_giw(ineq=2), axis=1)[..., niv_dmft - niv_cut : niv_dmft + niv_cut] + giw_spin_mean = 0.5 * (giw_spin_mean + giw_spin_mean_2) + del giw_spin_mean_2 + giw_spin_mean_large = np.zeros((2 * giw_spin_mean.shape[0], *giw_spin_mean.shape[1:])) + giw_spin_mean_large[0:2] = giw_spin_mean + giw_spin_mean_large[2:4] = giw_spin_mean + g_dmft = GreensFunction(extend_orbital(giw_spin_mean_large)) + del giw_spin_mean_large + siw_spin_mean = np.mean(file.get_siw(), axis=1) # [band,spin,niv] siw_spin_mean = extend_orbital(siw_spin_mean)[None, None, None, ...] siw_dc_spin_mean = np.mean(file.get_dc(), axis=-1) # [band,spin] @@ -81,6 +133,28 @@ def extend_orbital(arr: np.ndarray) -> np.ndarray: siw_spin_mean = siw_spin_mean[..., niv_dmft - niv_cut : niv_dmft + niv_cut] sigma_dmft = SelfEnergy(siw_spin_mean, estimate_niv_core=True) + siw_dc_spin_mean + if combine_two_atoms_to_one_obj: + siw_spin_mean_2 = np.mean(file.get_siw(ineq=2), axis=1) + siw_spin_mean_2 = extend_orbital(siw_spin_mean_2)[None, None, None, ...] + siw_dc_spin_mean_2 = np.mean(file.get_dc(ineq=2), axis=-1) + siw_dc_spin_mean_2 = extend_orbital(siw_dc_spin_mean_2)[None, None, None, ..., None] + siw_spin_mean_2 = siw_spin_mean_2[..., niv_dmft - niv_cut : niv_dmft + niv_cut] + siw_spin_mean = 0.5 * (siw_spin_mean + siw_spin_mean_2) + del siw_spin_mean_2 + siw_dc_spin_mean = 0.5 * (siw_dc_spin_mean + siw_dc_spin_mean_2) + del siw_dc_spin_mean_2 + + siw_spin_mean_large = np.zeros( + (1, 1, 1, 2 * siw_spin_mean.shape[3], 2 * siw_spin_mean.shape[3], siw_spin_mean.shape[-1]) + ) + siw_spin_mean_large[:, :, :, 0:2, 0:2, ...] = siw_spin_mean + siw_spin_mean_large[:, :, :, 2:4, 2:4, ...] = siw_spin_mean + siw_dc_spin_mean_large = np.zeros_like(siw_spin_mean_large) + siw_dc_spin_mean_large[:, :, :, 0:2, 0:2, ...] = siw_dc_spin_mean + siw_dc_spin_mean_large[:, :, :, 2:4, 2:4, ...] = siw_dc_spin_mean + sigma_dmft = SelfEnergy(siw_spin_mean_large, estimate_niv_core=True) + siw_dc_spin_mean_large + del siw_spin_mean_large, siw_dc_spin_mean_large + del giw_spin_mean, siw_spin_mean, siw_dc_spin_mean file.close() @@ -198,11 +272,17 @@ def set_hamiltonian(er_type: str, er_input: str | list, int_type: str, int_input if int_type == "one_band_from_dmft" or int_type == "" or int_type is None: return ham.single_band_interaction(config.lattice.interaction.udd) elif int_type == "kanamori_from_dmft": - return ham.kanamori_interaction( - config.sys.n_bands, + return ham.kanamori_interaction_dp( + config.sys.nd_bands, + config.sys.np_bands, config.lattice.interaction.udd, + config.lattice.interaction.upp, + config.lattice.interaction.udp, config.lattice.interaction.jdd, + config.lattice.interaction.jpp, + config.lattice.interaction.jdp, config.lattice.interaction.vdd, + config.lattice.interaction.vpp, ) elif int_type == "custom": if not isinstance(int_input, str): diff --git a/moldga/hamiltonian.py b/moldga/hamiltonian.py index 7505107..1d838bd 100644 --- a/moldga/hamiltonian.py +++ b/moldga/hamiltonian.py @@ -107,26 +107,81 @@ def interaction_orbital_diagonal(self, u: float, n_bands: int = 1) -> "Hamiltoni return self._add_interaction_term(interaction_elements) - def kanamori_interaction(self, n_bands: int, udd: float, jdd: float, vdd: float = None) -> "Hamiltonian": + def kanamori_interaction_d(self, n_bands: int, udd: float, jdd: float, vdd: float = None) -> "Hamiltonian": """ - Adds the Kanamori interaction terms to the Hamiltonian. The interaction terms are defined by the Hubbard 'udd' (U), - the exchange interaction 'jdd' (J) and the pair hopping 'vdd' (V or sometimes U'). 'vdd' is an optional parameter, if left empty, - it is set to U-2J. + Adds the Kanamori interaction terms ONLY for d orbitals to the Hamiltonian. + The interaction terms are defined by the Hubbard 'udd' (U), + the exchange interaction 'jdd' (J) and the pair hopping 'vdd' (V or sometimes U'). + 'vdd' is an optional parameter, if left empty, it is set to V=U-2J. """ + return self.kanamori_interaction_dp(nd_bands=n_bands, udd=udd, jdd=jdd, vdd=vdd) + + def kanamori_interaction_p(self, n_bands: int, upp: float, jpp: float, vpp: float = None) -> "Hamiltonian": + """ + Adds the Kanamori interaction terms ONLY for p orbitals to the Hamiltonian. + The interaction terms are defined by the Hubbard 'udd' (U), + the exchange interaction 'jpp' (J) and the pair hopping 'vpp' (V or sometimes U'). + 'vpp' is an optional parameter, if left empty, it is set to V=U-2J. + """ + return self.kanamori_interaction_dp(np_bands=n_bands, upp=upp, jpp=jpp, vpp=vpp) + + def kanamori_interaction_dp( + self, + nd_bands: int = 0, + np_bands: int = 0, + udd: float = 0.0, + upp: float = 0.0, + udp: float = 0.0, + jdd: float = 0.0, + jpp: float = 0.0, + jdp: float = 0.0, + vdd: float = None, + vpp: float = None, + ) -> "Hamiltonian": + """ + Adds the full Kanamori interaction terms for d and p orbitals to the Hamiltonian. + The interaction terms are defined by the local interaction Hubbard U, + the exchange interaction J and the pair hopping V or sometimes U'. + vdd (vpp) (vdp) are optional parameters, if left empty, they are set to V=U-2J. + """ + + # Default V=U-2J if vdd is None: vdd = udd - 2 * jdd + if vpp is None: + vpp = upp - 2 * jpp + + def orb_type(i: int) -> str: + return "d" if i < nd_bands else "p" # d bands come first + + def get_params(o1: int, o2: int) -> tuple[float, float, float]: + # Return correct (U, J, V) depending on orbital types of o1, o2 + t1, t2 = orb_type(o1), orb_type(o2) + + if t1 == "d" and t2 == "d": + return udd, jdd, vdd + elif t1 == "p" and t2 == "p": + return upp, jpp, vpp + else: + return 0, jdp, udp # udp is used as "Vdp", since there is no Udp possible that fits U_llll r_loc = [0, 0, 0] + n_tot = nd_bands + np_bands interaction_elements = [] - for a, b, c, d in it.product(range(n_bands), repeat=4): + + for a, b, c, d in it.product(range(n_tot), repeat=4): bands = [a + 1, b + 1, c + 1, d + 1] + + # choose parameters based on (a,b) pair (equivalently (c,d)) + u, j, v = get_params(a, b) + if a == b == c == d: # U_{llll} - interaction_elements.append(InteractionElement(r_loc, bands, udd)) - elif (a == d and b == c) or (a == b and c == d): # U_{lmml} or U_{llmm} - interaction_elements.append(InteractionElement(r_loc, bands, jdd)) + interaction_elements.append(InteractionElement(r_loc, bands, u)) + elif (a == d and b == c) or (a == b and c == d): # U_{lmml}, U_{llmm} + interaction_elements.append(InteractionElement(r_loc, bands, j)) elif a == c and b == d: # U_{lmlm} - interaction_elements.append(InteractionElement(r_loc, bands, vdd)) + interaction_elements.append(InteractionElement(r_loc, bands, v)) return self._add_interaction_term(interaction_elements) diff --git a/moldga/symmetrize_new.py b/moldga/symmetrize_new.py index 42d37fb..53eac82 100644 --- a/moldga/symmetrize_new.py +++ b/moldga/symmetrize_new.py @@ -8,15 +8,17 @@ import glob import itertools as it import os +import re import readline +from pathlib import Path import h5py import numpy as np -def index2component_general(num_bands: int, n: int, ind: int) -> tuple[np.ndarray, np.ndarray, np.ndarray]: +def index2component_general_4(num_bands: int, n: int, ind: int) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """ - Returns the band and spin components corresponding to a compound index. + Returns the band and spin components corresponding to a compound index for four-legged objects. """ bandspin = np.zeros(n, dtype=np.int_) spin = np.zeros(n, dtype=np.int_) @@ -33,9 +35,23 @@ def index2component_general(num_bands: int, n: int, ind: int) -> tuple[np.ndarra return bandspin, band, spin -def component2index_general(num_bands: int, bands: list, spins: list) -> int: +def component2index_general_2(num_bands: int, bands: list, spins: list) -> int: """ - Computes a compound index from band and spin indices. + Computes a compound index from band and spin indices for a two-legged object. + """ + assert num_bands > 0, "Number of bands has to be set to non-zero positive integers." + + n_spins = 2 + dims_bs = 2 * (num_bands * n_spins,) + dims_1 = (num_bands, n_spins) + + bandspin = np.ravel_multi_index((bands, spins), dims_1) + return np.ravel_multi_index(bandspin, dims_bs) + 1 + + +def component2index_general_4(num_bands: int, bands: list, spins: list) -> int: + """ + Computes a compound index from band and spin indices for a four-legged object. """ assert num_bands > 0, "Number of bands has to be set to non-zero positive integers." @@ -47,9 +63,9 @@ def component2index_general(num_bands: int, bands: list, spins: list) -> int: return np.ravel_multi_index(bandspin, dims_bs) + 1 -def index2component_band(num_bands: int, n: int, ind: int) -> list: +def index2component_band_4(num_bands: int, n: int, ind: int) -> list: """ - Computes only orbital indices from a compound index. + Computes only orbital indices from a compound index for four-legged objects. """ b = [] ind_tmp = ind - 1 @@ -59,9 +75,9 @@ def index2component_band(num_bands: int, n: int, ind: int) -> list: return b -def component2index_band(num_bands: int, n: int, b: list) -> int: +def component2index_band_4(num_bands: int, n: int, b: list) -> int: """ - Computes a compound index from orbital indices only. + Computes a compound index from only orbital indices for four-legged objects. """ ind = 1 for i in range(n): @@ -69,22 +85,81 @@ def component2index_band(num_bands: int, n: int, b: list) -> int: return ind -def get_worm_components(num_bands: int) -> list[int]: +def _get_worm_components_2(num_bands: int, orbs: list[list[int]]) -> list[int]: """ - Returns the list of worm components for a given number of bands, where only relevant spin combinations for the - density and magnetic channels in the case of SU(2) symmetry are picked. If one wants to speed up the w2dynamics - simulation, one can furthermore restrict the worm components to only allow spins = [0, 0, 0, 0], [0, 0, 1, 1] at - the cost of more stochastic noise. + Returns the worm components for two-legged objects. + """ + spins = [0, 0], [1, 1] + component_indices = [] + for o in orbs: + for s in spins: + component_indices.append(int(component2index_general_2(num_bands, o, s))) + return sorted(component_indices) + + +def get_worm_components_2(num_bands: int) -> list[int]: + """ + Returns the list of worm components for a given number of bands for two-legged objects, + where only relevant spin combinations for the + density and magnetic channels in the case of SU(2) symmetry are picked. + """ + orbs = [list(orb) for orb in it.product(range(num_bands), repeat=2)] + return _get_worm_components_2(num_bands, orbs) + + +def get_worm_components_partial_2(num_bands: int) -> list[int]: + """ + Returns the list of worm components for a given number of bands for two-legged objects, + where only relevant spin combinations for the + density and magnetic channels in the case of SU(2) symmetry are picked. + It only lists orbital-diagonal components. + """ + orbs = [[orb, orb] for orb in range(num_bands)] + return _get_worm_components_2(num_bands, orbs) + + +def _get_worm_components_4(num_bands: int, orbs: list[list[int]]) -> list[int]: + """ + Returns the worm components for 4-legged objects. """ - orbs = [list(orb) for orb in it.product(range(num_bands), repeat=4)] spins = [0, 0, 0, 0], [1, 1, 1, 1], [0, 0, 1, 1], [1, 1, 0, 0], [1, 0, 0, 1], [0, 1, 1, 0] component_indices = [] for o in orbs: for s in spins: - component_indices.append(int(component2index_general(num_bands, o, s))) + component_indices.append(int(component2index_general_4(num_bands, o, s))) return sorted(component_indices) +def get_worm_components_4(num_bands: int) -> list[int]: + """ + Returns the list of worm components for a given number of bands for four-legged objecst, + where only relevant spin combinations for the + density and magnetic channels in the case of SU(2) symmetry are picked. + """ + orbs = [list(orb) for orb in it.product(range(num_bands), repeat=4)] + return _get_worm_components_4(num_bands, orbs) + + +def get_worm_components_partial_4(num_bands: int) -> list[int]: + """ + Returns the list of worm components for a given number of bands for four-legged objects, + where only relevant spin combinations for the + density and magnetic channels in the case of SU(2) symmetry are picked. + It only lists worm components where the orbitals are not of type ijjj, jijj, jjij or jjji. + """ + orbs = [ + list(orb) + for orb in it.product(range(num_bands), repeat=4) + if not ( + orb[1] == orb[2] == orb[3] != orb[0] # ijjj + or orb[0] == orb[2] == orb[3] != orb[1] # jijj + or orb[0] == orb[1] == orb[3] != orb[2] # jjij + or orb[0] == orb[1] == orb[2] != orb[3] # jjji + ) + ] + return _get_worm_components_4(num_bands, orbs) + + def extract_g2_general(group_string: str, indices: list, file: h5py.File, niw: int, niv: int) -> tuple: r""" Extracts the two-particle Green's function components from the vertex file for given indices and group string. @@ -102,7 +177,7 @@ def extract_g2_general(group_string: str, indices: list, file: h5py.File, niw: i elements = elements.transpose(0, 1, 3, 2) # construct G2dens and G2magn the output file - bands, spins = zip(*(index2component_general(n_bands, 4, int(i))[1:3] for i in indices)) + bands, spins = zip(*(index2component_general_4(n_bands, 4, int(i))[1:3] for i in indices)) # since we are SU(2) symmetric, we only have to pick out the elements where the spin is either # [0,0,0,0] or [1,1,1,1] for uu component, [0,0,1,1] or [1,1,0,0] for ud component and [0,1,1,0] or [1,0,0,1] for ud_bar component @@ -143,16 +218,16 @@ def extract_g2_general(group_string: str, indices: list, file: h5py.File, niw: i return g2_uuuu, g2_dddd, g2_dduu, g2_uudd, g2_uddu, g2_duud -def save_to_file(g2_list: list[np.ndarray], names: list[str], niw: int, nb: int): +def save_to_file(g2_list: list[np.ndarray], names: list[str], niw: int, nb: int, ineq: int): """ Saves the given g2 to the output file. """ assert len(g2_list) == len(names) for wn in range(2 * niw + 1): for i, j, k, l in it.product(range(nb), repeat=4): - idx = component2index_band(nb, 4, [i, j, k, l]) + idx = component2index_band_4(nb, 4, [i, j, k, l]) for g2, name in zip(g2_list, names): - output_file[f"ineq-001/{name}/{wn:05}/{idx:05}/value"] = g2[i, j, k, l, wn].transpose() + output_file[f"ineq-{ineq:03}/{name}/{wn:05}/{idx:05}/value"] = g2[i, j, k, l, wn].transpose() def get_niw_niv(vertex_file, g4iw_groupstring, indices): @@ -179,8 +254,6 @@ def complete(text, state): default_filename = "Vertex.hdf5" default_output_filename = "g4iw_sym.hdf5" - g4iw_groupstring = "worm-last/ineq-001/g4iw-worm" - readline.parse_and_bind("tab: complete") readline.set_completer_delims(" \t\n;") readline.set_completer(complete) @@ -191,42 +264,68 @@ def complete(text, state): input_filename = input_filename.strip() if input_filename else default_filename output_filename = output_filename.strip() if output_filename else default_output_filename + input_filename = str(Path(input_filename).with_suffix(".hdf5")) + output_filename = str(Path(output_filename).with_suffix(".hdf5")) + vertex_file = h5py.File(input_filename, "r") output_file = h5py.File(output_filename, "w") - n_bands = int(vertex_file[".config"].attrs[f"atoms.1.nd"]) + int(vertex_file[".config"].attrs[f"atoms.1.np"]) + group = vertex_file["worm-last"] - indices = None - try: - indices = list(vertex_file[g4iw_groupstring].keys()) - except KeyError: - print("WARNING: No g4iw-worm group found in the input file. Aborting.") - exit() + ineq_numbers = [] + for key in group.keys(): + match = re.match(r"ineq-(\d+)", key) + if match: + ineq_numbers.append(int(match.group(1))) - niw, niv = get_niw_niv(vertex_file, g4iw_groupstring, indices) + ineq_numbers.sort() - print("Number of bands:", n_bands) - print("Number of fermionic Matsubara frequencies:", niv) - print("Number of bosonic Matsubara frequencies:", niw) + n_bands = int(vertex_file[".config"].attrs[f"atoms.1.nd"]) + int(vertex_file[".config"].attrs[f"atoms.1.np"]) - print("Extracting G2 ...") - g2_uuuu, g2_dddd, g2_dduu, g2_uudd, g2_uddu, g2_duud = extract_g2_general( - g4iw_groupstring, indices, vertex_file, niw, niv - ) - print("G2 extracted. Calculating G2_dens and G2_magn ...") - g2_dens = 0.5 * (g2_uuuu + g2_dddd + g2_uudd + g2_dduu) - g2_magn = 0.5 * (g2_uddu + g2_duud) + for ineq in ineq_numbers: + print("-----------------------------------------") + print("Processing inequivalent atom number:", ineq) + print("-----------------------------------------") + g4iw_groupstring = f"worm-last/ineq-{ineq:03}/g4iw-worm" + + indices = None + try: + indices = list(vertex_file[g4iw_groupstring].keys()) + except KeyError: + if ineq == max(ineq_numbers): + print(f"WARNING: No g4iw-worm group found for atom {ineq} in the input file. Aborting.") + exit() + else: + print( + f"WARNING: No g4iw-worm group found for atom {ineq} in the input file. Continuing with next atom." + ) + continue + + niw, niv = get_niw_niv(vertex_file, g4iw_groupstring, indices) + + print("Number of bands:", n_bands) + print("Number of fermionic Matsubara frequencies:", niv) + print("Number of bosonic Matsubara frequencies:", niw) + + print(f"Extracting G2 for atom {ineq} ...") + g2_uuuu, g2_dddd, g2_dduu, g2_uudd, g2_uddu, g2_duud = extract_g2_general( + g4iw_groupstring, indices, vertex_file, niw, niv + ) + print(f"G2 extracted. Calculating G2_dens and G2_magn for atom {ineq} ...") + g2_dens = 0.5 * (g2_uuuu + g2_dddd + g2_uudd + g2_dduu) + g2_magn = 0.5 * (g2_uddu + g2_duud) - del g2_uuuu, g2_dddd, g2_dduu, g2_uudd, g2_uddu, g2_duud - gc.collect() - print("G2_dens and G2_magn calculated. Writing to file ...") + del g2_uuuu, g2_dddd, g2_dduu, g2_uudd, g2_uddu, g2_duud + gc.collect() + print(f"G2_dens and G2_magn for atom {ineq} calculated. Writing to file ...") - save_to_file([g2_dens, g2_magn], ["dens", "magn"], niw, n_bands) - del g2_dens, g2_magn - gc.collect() - print("G2_dens and G2_magn successfully written to file.") + save_to_file([g2_dens, g2_magn], ["dens", "magn"], niw, n_bands, ineq) + del g2_dens, g2_magn + gc.collect() + print(f"G2_dens and G2_magn for atom {ineq} successfully written to file.") output_file.close() vertex_file.close() + print(f"{len(ineq_numbers)} inequivalent atom(s) written to {output_filename}.") print("Done!") exit() diff --git a/moldga/w2dyn_aux.py b/moldga/w2dyn_aux.py index ac59aff..8be86a4 100644 --- a/moldga/w2dyn_aux.py +++ b/moldga/w2dyn_aux.py @@ -34,23 +34,23 @@ def open(self): """ self._file = h5py.File(self._fname, "r") - def atom_group(self, dmft_iter="dmft-last", atom=1): + def ineq_group(self, dmft_iter="dmft-last", ineq=1): """ - Returns the group string for a given DMFT iteration and atom. + Returns the group string for a given DMFT iteration and ineq. """ - return dmft_iter + f"/ineq-{atom:03}" + return dmft_iter + f"/ineq-{ineq:03}" - def get_nd(self, atom: int = 1) -> int: + def get_nd(self, ineq: int = 1) -> int: """ Returns the number of d orbitals of the DMFT calculation. """ - return self._from_atom_config("nd", atom=atom) + return self._from_ineq_config("nd", ineq=ineq) - def get_np(self, atom: int = 1) -> int: + def get_np(self, ineq: int = 1) -> int: """ Returns the number of p orbitals of the DMFT calculation. """ - return self._from_atom_config("np", atom=atom) + return self._from_ineq_config("np", ineq=ineq) def get_beta(self) -> float: r""" @@ -70,107 +70,107 @@ def get_totdens(self) -> float: """ return self._file[".config"].attrs["general.totdens"] - def get_jdd(self, atom: int = 1) -> float: + def get_jdd(self, ineq: int = 1) -> float: """ Extracts the Hund's coupling for d orbitals. """ - return self._from_atom_config("jdd", atom=atom) + return self._from_ineq_config("jdd", ineq=ineq) - def get_jdp(self, atom: int = 1) -> float: + def get_jdp(self, ineq: int = 1) -> float: """ Extracts the Hund's coupling between d and p orbitals. """ - return self._from_atom_config("jdp", atom=atom) + return self._from_ineq_config("jdp", ineq=ineq) - def get_jpp(self, atom: int = 1) -> float: + def get_jpp(self, ineq: int = 1) -> float: """ Extracts the Hund's coupling for p orbitals. """ - return self._from_atom_config("jpp", atom=atom) + return self._from_ineq_config("jpp", ineq=ineq) - def get_jppod(self, atom: int = 1) -> float: + def get_jppod(self, ineq: int = 1) -> float: """ Extracts the offdiagonal terms for jpp. """ - return self._from_atom_config("jppod", atom=atom) + return self._from_ineq_config("jppod", ineq=ineq) - def get_udd(self, atom: int = 1) -> float: + def get_udd(self, ineq: int = 1) -> float: """ Extracts the Hubbard U for d orbitals. """ - return self._from_atom_config("udd", atom=atom) + return self._from_ineq_config("udd", ineq=ineq) - def get_udp(self, atom: int = 1) -> float: + def get_udp(self, ineq: int = 1) -> float: """ Extracts the Hubbard U between d and p orbitals. """ - return self._from_atom_config("udp", atom=atom) + return self._from_ineq_config("udp", ineq=ineq) - def get_upp(self, atom: int = 1) -> float: + def get_upp(self, ineq: int = 1) -> float: """ Extracts the Hubbard U for p orbitals. """ - return self._from_atom_config("upp", atom=atom) + return self._from_ineq_config("upp", ineq=ineq) - def get_uppod(self, atom: int = 1): + def get_uppod(self, ineq: int = 1): """ Extracts the offdiagonal terms for upp. """ - return self._from_atom_config("uppod", atom=atom) + return self._from_ineq_config("uppod", ineq=ineq) - def get_vdd(self, atom: int = 1) -> float: + def get_vdd(self, ineq: int = 1) -> float: """ Extracts the intersite interaction between d orbitals. """ - return self._from_atom_config("vdd", atom=atom) + return self._from_ineq_config("vdd", ineq=ineq) - def get_vpp(self, atom: int = 1) -> float: + def get_vpp(self, ineq: int = 1) -> float: """ Extracts the intersite interaction between p orbitals. """ - return self._from_atom_config("vpp", atom=atom) + return self._from_ineq_config("vpp", ineq=ineq) - def get_siw(self, dmft_iter: str = "dmft-last", atom: int = 1) -> list: + def get_siw(self, dmft_iter: str = "dmft-last", ineq: int = 1) -> list: """ Extracts the DMFT self-energy in Matsubara frequency space as [band, spin, iv]. """ - return self._file[self.atom_group(dmft_iter=dmft_iter, atom=atom) + "/siw/value"][()] + return self._file[self.ineq_group(dmft_iter=dmft_iter, ineq=ineq) + "/siw/value"][()] - def get_giw(self, dmft_iter: str = "dmft-last", atom: int = 1) -> list: + def get_giw(self, dmft_iter: str = "dmft-last", ineq: int = 1) -> list: """ Extracts the DMFT Green's function in Matsubara frequency space as [band, spin, iv]. """ - return self._file[self.atom_group(dmft_iter=dmft_iter, atom=atom) + "/giw/value"][()] + return self._file[self.ineq_group(dmft_iter=dmft_iter, ineq=ineq) + "/giw/value"][()] - def get_occ(self, dmft_iter: str = "dmft-last", atom: int = 1) -> list: + def get_occ(self, dmft_iter: str = "dmft-last", ineq: int = 1) -> list: """ Extracts the occupation matrix as [band1, spin1, band2, spin2]. """ - return self._file[self.atom_group(dmft_iter=dmft_iter, atom=atom) + "/occ/value"][()] + return self._file[self.ineq_group(dmft_iter=dmft_iter, ineq=ineq) + "/occ/value"][()] - def get_rho1(self, dmft_iter: str = "dmft-last", atom: int = 1) -> list: + def get_rho1(self, dmft_iter: str = "dmft-last", ineq: int = 1) -> list: """ Extracts the 1-particle density matrix as [band1, spin1, band2, spin2]. """ - return self._file[self.atom_group(dmft_iter=dmft_iter, atom=atom) + "/rho1/value"][()] + return self._file[self.ineq_group(dmft_iter=dmft_iter, ineq=ineq) + "/rho1/value"][()] - def get_rho2(self, dmft_iter: str = "dmft-last", atom: int = 1) -> list: + def get_rho2(self, dmft_iter: str = "dmft-last", ineq: int = 1) -> list: """ Extracts the 2-particle density matrix as [band1, spin1, band2, spin2, band3, spin3, band4, spin4]. """ - return self._file[self.atom_group(dmft_iter=dmft_iter, atom=atom) + "/rho2/value"][()] + return self._file[self.ineq_group(dmft_iter=dmft_iter, ineq=ineq) + "/rho2/value"][()] - def _from_atom_config(self, key: str, atom: int = 1): + def _from_ineq_config(self, key: str, ineq: int = 1): """ - Extracts a value from the .config group for a given atom. + Extracts a value from the .config group for a given ineq. """ - return self._file[".config"].attrs[f"atoms.{atom:1}.{key}"] + return self._file[".config"].attrs[f"atoms.{ineq:1}.{key}"] - def get_dc(self, dmft_iter: str = "dmft-last", atom: int = 1) -> list: + def get_dc(self, dmft_iter: str = "dmft-last", ineq: int = 1) -> list: """ Extracts the DMFT double-counting correction as [band, spin]. """ - return self._file[self.atom_group(dmft_iter=dmft_iter, atom=atom) + "/dc/value"][()] + return self._file[self.ineq_group(dmft_iter=dmft_iter, ineq=ineq) + "/dc/value"][()] class W2dynG4iwFile: @@ -210,11 +210,11 @@ def read_g2_full_multiband(self, n_bands: int, ineq: int = 1, name: str = "dens" self._file[f"{channel_group_string}/00000/{first_index:05}/value"][()] ) # extract niv from the size of the array - g2 = np.zeros((n_bands,) * 4 + (niw_full,) + 2 * (niv_full,), dtype=np.complex64) + g2 = np.zeros((n_bands,) * 4 + (niw_full,) + 2 * (niv_full,), dtype=np.complex128) for wn in range(niw_full): wn_group_string = f"{channel_group_string}/{wn:05}" for ind in self._file[wn_group_string].keys(): - bands = sym.index2component_band(n_bands, 4, int(ind)) + bands = sym.index2component_band_4(n_bands, 4, int(ind)) val = self._file[f"{wn_group_string}/{ind}/value"][()].T g2[bands[0], bands[1], bands[2], bands[3], wn, ...] = val diff --git a/tests/test_dga_io.py b/tests/test_dga_io.py new file mode 100644 index 0000000..4b4e731 --- /dev/null +++ b/tests/test_dga_io.py @@ -0,0 +1,215 @@ +import numpy as np +import pytest +import moldga.dga_io as dga_io + + +class DummyW2dynFile: + def __init__(self, *args, **kwargs): + pass + + def get_beta(self): + return 10.0 + + def get_mu(self): + return 0.5 + + def get_nd(self): + return 1 + + def get_np(self): + return 1 + + def get_totdens(self): + return 1.0 + + # interactions (not relevant but required) + def get_udd(self): + return 1.0 + + def get_udp(self): + return 1.0 + + def get_upp(self): + return 1.0 + + def get_uppod(self): + return 1.0 + + def get_jdd(self): + return 1.0 + + def get_jdp(self): + return 1.0 + + def get_jpp(self): + return 1.0 + + def get_jppod(self): + return 1.0 + + def get_vdd(self): + return 1.0 + + def get_vpp(self): + return 1.0 + + # distinguish ineq=1 and ineq=2 clearly + def get_rho1(self, ineq=1, **kwargs): + val = 1.0 if ineq == 1 else 3.0 + return val * np.ones((2, 2, 2, 2)) + + def get_giw(self, ineq=1, **kwargs): + val = 2.0 if ineq == 1 else 6.0 + return val * np.ones((2, 2, 40)) + + def get_siw(self, ineq=1, **kwargs): + val = 5.0 if ineq == 1 else 7.0 + return val * np.ones((2, 2, 40)) + + def get_dc(self, ineq=1, **kwargs): + val = 1.0 if ineq == 1 else 3.0 + return val * np.ones((2, 2)) + + def close(self): + pass + + +class DummyW2dynG4iwFile: + def __init__(self, *args, **kwargs): + pass + + def read_g2_full_multiband(self, n_bands, ineq=1, **kwargs): + val = 10.0 if ineq == 1 else 20.0 + return val * np.ones((n_bands, n_bands, n_bands, n_bands, 5, 6, 6)) + + def close(self): + pass + + +class DummyLocalFourPoint: + def __init__(self, mat, channel=None): + self.mat = mat + self.niw = mat.shape[4] + self.niv = mat.shape[5] + + def cut_niw_and_niv(self, *args): + return self + + def symmetrize_orbitals(self, *args): + return self + + def symmetrize_v_vp(self): + return self + + +class DummyGreensFunction: + def __init__(self, data): + self.data = data + + def symmetrize_orbitals(self, *args): + return self + + +class DummySelfEnergy: + def __init__(self, data, estimate_niv_core=False): + self.data = data + + def __add__(self, other): + return self + + def symmetrize_orbitals(self, *args): + return self + + +@pytest.fixture +def cfg(tmp_path): + class Dummy: + pass + + cfg = Dummy() + + cfg.sys = Dummy() + cfg.sys.n = 1 + cfg.sys.nd_bands = 1 + cfg.sys.np_bands = 1 + cfg.sys.n_bands = 2 + cfg.sys.occ_dmft = None + + cfg.box = Dummy() + cfg.box.niw_core = 2 + cfg.box.niv_core = 2 + cfg.box.niv_shell = 1 + cfg.box.niv_full = None + + cfg.dmft = Dummy() + cfg.dmft.input_path = "" + cfg.dmft.fname_1p = "f1" + cfg.dmft.fname_2p = "f2" + cfg.dmft.symmetrize_orbitals = None + cfg.dmft.do_sym_v_vp = False + + cfg.lattice = Dummy() + cfg.lattice.interaction = Dummy() + cfg.lattice.type = "t_tp_tpp" + cfg.lattice.er_input = [1, 0, 0] + cfg.lattice.interaction_type = "" + cfg.lattice.interaction_input = None + cfg.lattice.orbital_basis = None + + cfg.lattice.k_grid = Dummy() + cfg.lattice.k_grid.nk_tot = 1 + cfg.lattice.k_grid.specify_orbital_basis = lambda *a: None + + cfg.lattice.q_grid = Dummy() + cfg.lattice.q_grid.nk_tot = 1 + cfg.lattice.q_grid.specify_orbital_basis = lambda *a: None + + cfg.output = Dummy() + cfg.output.output_path = str(tmp_path) + cfg.output.plotting_subfolder_name = "Plots" + cfg.output.save_quantities = False + cfg.output.do_plotting = False + + cfg.eliashberg = Dummy() + cfg.eliashberg.subfolder_name = "Eliashberg" + cfg.eliashberg.perform_eliashberg = False + + cfg.logger = Dummy() + cfg.logger.info = lambda *a: None + + return cfg + + +def test_combine_two_atoms_branch(monkeypatch, cfg): + monkeypatch.setattr(dga_io, "config", cfg) + monkeypatch.setattr(dga_io.w2dyn_aux, "W2dynFile", DummyW2dynFile) + monkeypatch.setattr(dga_io.w2dyn_aux, "W2dynG4iwFile", DummyW2dynG4iwFile) + monkeypatch.setattr(dga_io, "LocalFourPoint", DummyLocalFourPoint) + monkeypatch.setattr(dga_io, "GreensFunction", DummyGreensFunction) + monkeypatch.setattr(dga_io, "SelfEnergy", DummySelfEnergy) + monkeypatch.setattr(dga_io, "set_hamiltonian", lambda *a, **k: "ham") + + g, sigma, g2_dens, g2_magn = dga_io.load_from_w2dyn_file_and_update_config(combine_two_atoms_to_one_obj=True) + + assert cfg.sys.n_bands == 4 + + mat = g2_dens.mat + + # diagonal blocks + assert np.allclose(mat[0:2, 0:2, 0:2, 0:2], 15.0) + assert np.allclose(mat[2:4, 2:4, 2:4, 2:4], 15.0) + + # off-diagonal blocks should be zero + assert np.allclose(mat[0:2, 0:2, 2:4, 2:4], 0.0) + + g_data = g.data + assert np.allclose(g_data[0, 0], g_data[2, 2]) # duplicated blocks + + sigma_data = sigma.data + assert sigma_data.shape[3] == 4 # orbital doubled + + occ = cfg.sys.occ_dmft + assert np.allclose(occ[0:2, 0:2], 4.0) + assert np.allclose(occ[2:4, 2:4], 4.0) + assert np.allclose(occ[0:2, 2:4], 0.0) + assert np.allclose(occ[2:4, 0:2], 0.0) diff --git a/tests/test_hamiltonian.py b/tests/test_hamiltonian.py index 70b0a78..7a99e60 100644 --- a/tests/test_hamiltonian.py +++ b/tests/test_hamiltonian.py @@ -76,12 +76,12 @@ def test_single_band_interaction_sets_correct_u(): def test_kanamori_interaction_defaults_1_band(): - h = Hamiltonian().kanamori_interaction(n_bands=1, udd=5.0, jdd=1.0) + h = Hamiltonian().kanamori_interaction_d(n_bands=1, udd=5.0, jdd=1.0) assert np.isclose(h._ur_local[0, 0, 0, 0], 5.0) def test_kanamori_interaction_with_vdd_1_band(): - h = Hamiltonian().kanamori_interaction(n_bands=1, udd=5.0, jdd=1.0, vdd=2.0) + h = Hamiltonian().kanamori_interaction_d(n_bands=1, udd=5.0, jdd=1.0, vdd=2.0) assert np.isclose(h._ur_local[0, 0, 0, 0], 5.0) @@ -92,7 +92,7 @@ def test_kanamori_interaction_with_vdd_2_band(): "vdd": np.random.rand(), } - h = Hamiltonian().kanamori_interaction(n_bands=2, **params) + h = Hamiltonian().kanamori_interaction_d(n_bands=2, **params) assert np.isclose(h._ur_local[0, 0, 0, 0], params["udd"]) assert np.isclose(h._ur_local[1, 1, 1, 1], params["udd"]) @@ -206,3 +206,120 @@ def test_read_write_hr_hk_files(): wannier_hk_twoband, _ = Hamiltonian().read_hk_w2k(f"{folder}/wannier_twoband_24x24.hk") ek_ref = wannier_hk_twoband.get_ek(k_grid).reshape(ek.shape) assert np.allclose(ek, ek_ref) + + +def test_kanamori_d_basic(): + ham = Hamiltonian() + n = 3 + u_val = 4.0 + j = 1.0 + + ham.kanamori_interaction_d(n_bands=n, udd=u_val, jdd=j) + u = ham.get_local_u() + + v = u_val - 2 * j + + for a in range(n): + for b in range(n): + for c in range(n): + for d in range(n): + if a == b == c == d: + assert np.isclose(u[a, b, c, d], u_val) + elif (a == d and b == c) or (a == b and c == d): + assert np.isclose(u[a, b, c, d], j) + elif a == c and b == d: + assert np.isclose(u[a, b, c, d], v) + else: + assert np.isclose(u[a, b, c, d], 0.0) + + +def test_kanamori_p_basic(): + ham = Hamiltonian() + n = 2 + u_val = 3.0 + j = 0.5 + + ham.kanamori_interaction_p(n_bands=n, upp=u_val, jpp=j) + u = ham.get_local_u() + + v = u_val - 2 * j + + for a in range(n): + for b in range(n): + for c in range(n): + for d in range(n): + if a == b == c == d: + assert np.isclose(u[a, b, c, d], u_val) + elif (a == d and b == c) or (a == b and c == d): + assert np.isclose(u[a, b, c, d], j) + elif a == c and b == d: + assert np.isclose(u[a, b, c, d], v) + else: + assert np.isclose(u[a, b, c, d], 0.0) + + +def test_kanamori_dp_block_structure(): + ham = Hamiltonian() + + nd, npb = 2, 2 + + udd, jdd = 8.0, 1.0 + upp, jpp = 4.0, 0.5 + udp, jdp = 2.0, 0.2 + + ham.kanamori_interaction_dp(nd_bands=nd, np_bands=npb, udd=udd, upp=upp, udp=udp, jdd=jdd, jpp=jpp, jdp=jdp) + + u = ham.get_local_u().mat + vdd = udd - 2 * jdd + vpp = upp - 2 * jpp + + def is_d(i): + return i < nd + + for a in range(nd + npb): + for b in range(nd + npb): + for c in range(nd + npb): + for d in range(nd + npb): + + if is_d(a) and is_d(b): + uu, jj, vv = udd, jdd, vdd + elif (not is_d(a)) and (not is_d(b)): + uu, jj, vv = upp, jpp, vpp + else: + uu, jj, vv = 0, jdp, udp + + if a == b == c == d: + assert np.isclose(u[a, b, c, d], uu) + elif (a == d and b == c) or (a == b and c == d): + assert np.isclose(u[a, b, c, d], jj) + elif a == c and b == d: + assert np.isclose(u[a, b, c, d], vv) + else: + assert np.isclose(u[a, b, c, d], 0.0) + + +def test_kanamori_dp_index_split(): + ham = Hamiltonian() + + nd, npb = 1, 1 + + ham.kanamori_interaction_dp(nd_bands=nd, np_bands=npb, udd=10.0, upp=5.0, udp=2.0, jdd=1.0, jpp=0.5, jdp=0.1) + + u = ham.get_local_u() + + assert np.isclose(u[0, 0, 0, 0], 10.0) + assert np.isclose(u[1, 1, 1, 1], 5.0) + + assert np.isclose(u[0, 1, 0, 1], 2.0) + assert np.isclose(u[0, 1, 1, 0], 0.1) + + +def test_kanamori_dp_no_unexpected_terms(): + ham = Hamiltonian() + + ham.kanamori_interaction_dp(nd_bands=2, np_bands=1, udd=6.0, upp=3.0, udp=1.0, jdd=1.0, jpp=0.5, jdp=0.2) + + u = ham.get_local_u() + + assert np.isclose(u[0, 1, 2, 0], 0.0) + assert np.isclose(u[2, 0, 1, 2], 0.0) diff --git a/tests/test_symmetrize.py b/tests/test_symmetrize.py index 0bd3fc5..513c5b4 100644 --- a/tests/test_symmetrize.py +++ b/tests/test_symmetrize.py @@ -12,25 +12,25 @@ @pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) def test_index2component_general_and_back(num_bands): for ind in range(1, 16 * num_bands**4 + 1): - bandspin, band, spin = index2component_general(num_bands, 4, ind) - ind_back = component2index_general(num_bands, list(band), list(spin)) + bandspin, band, spin = index2component_general_4(num_bands, 4, ind) + ind_back = component2index_general_4(num_bands, list(band), list(spin)) assert ind_back == ind @pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) def test_index2component_general_and_back_raises_if_index_too_large_or_too_small(num_bands): with pytest.raises(ValueError): - bandspin, band, spin = index2component_general(num_bands, 4, 16 * num_bands**4 + 1) - _ = component2index_general(num_bands, list(band), list(spin)) + bandspin, band, spin = index2component_general_4(num_bands, 4, 16 * num_bands**4 + 1) + _ = component2index_general_4(num_bands, list(band), list(spin)) with pytest.raises(ValueError): - bandspin, band, spin = index2component_general(num_bands, 4, 0) - _ = component2index_general(num_bands, list(band), list(spin)) + bandspin, band, spin = index2component_general_4(num_bands, 4, 0) + _ = component2index_general_4(num_bands, list(band), list(spin)) def test_component2index_general_invalid_num_bands(): with pytest.raises(AssertionError): - component2index_general(0, [0], [0]) + component2index_general_4(0, [0], [0]) @pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) @@ -38,14 +38,171 @@ def test_index2component_band_and_back(num_bands): orbs = list(it.product(range(num_bands), repeat=4)) for orb in orbs: - ind = component2index_band(num_bands, 4, list(orb)) - indices_back = index2component_band(num_bands, 4, ind) + ind = component2index_band_4(num_bands, 4, list(orb)) + indices_back = index2component_band_4(num_bands, 4, ind) assert indices_back == list(orb) @pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) def test_get_worm_components(num_bands): - result = get_worm_components(num_bands) + result = get_worm_components_4(num_bands) if num_bands == 1: assert result == [1, 4, 7, 10, 13, 16] assert len(result) == 6 * num_bands**4 + + +@pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) +def test_component2index_general_2_returns_int(num_bands): + result = component2index_general_2(num_bands, [0, 0], [0, 0]) + assert isinstance(result, (int, np.integer)) + + +@pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) +def test_component2index_general_2_index_within_range(num_bands): + max_index = (2 * num_bands) ** 2 + for b0 in range(num_bands): + for b1 in range(num_bands): + for s0, s1 in it.product([0, 1], repeat=2): + idx = component2index_general_2(num_bands, [b0, b1], [s0, s1]) + assert 1 <= idx <= max_index + + +@pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) +def test_component2index_general_2_all_indices_unique(num_bands): + indices = [ + component2index_general_2(num_bands, [b0, b1], [s0, s1]) + for b0 in range(num_bands) + for b1 in range(num_bands) + for s0, s1 in it.product([0, 1], repeat=2) + ] + assert len(indices) == len(set(indices)) + + +def test_component2index_general_2_invalid_num_bands_raises(): + with pytest.raises(AssertionError): + component2index_general_2(0, [0, 0], [0, 0]) + + +def test_component2index_general_2_single_band_covers_all_four_indices(): + results = {component2index_general_2(1, [0, 0], [s0, s1]) for s0, s1 in it.product([0, 1], repeat=2)} + assert results == {1, 2, 3, 4} + + +@pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) +def test_get_worm_components_2_is_sorted(num_bands): + result = get_worm_components_2(num_bands) + assert result == sorted(result) + + +@pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) +def test_get_worm_components_2_length(num_bands): + assert len(get_worm_components_2(num_bands)) == 2 * num_bands**2 + + +@pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) +def test_get_worm_components_2_no_duplicates(num_bands): + result = get_worm_components_2(num_bands) + assert len(result) == len(set(result)) + + +@pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) +def test_get_worm_components_2_all_indices_positive(num_bands): + assert all(idx >= 1 for idx in get_worm_components_2(num_bands)) + + +@pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) +def test_get_worm_components_partial_2_is_sorted(num_bands): + result = get_worm_components_partial_2(num_bands) + assert result == sorted(result) + + +@pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) +def test_get_worm_components_partial_2_length(num_bands): + assert len(get_worm_components_partial_2(num_bands)) == 2 * num_bands + + +@pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) +def test_get_worm_components_partial_2_no_duplicates(num_bands): + result = get_worm_components_partial_2(num_bands) + assert len(result) == len(set(result)) + + +@pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) +def test_get_worm_components_partial_2_is_subset_of_full(num_bands): + assert set(get_worm_components_partial_2(num_bands)).issubset(set(get_worm_components_2(num_bands))) + + +def test_get_worm_components_partial_2_single_band_matches_full(): + assert get_worm_components_partial_2(1) == get_worm_components_2(1) + + +@pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) +def test_get_worm_components_partial_4_is_sorted(num_bands): + result = get_worm_components_partial_4(num_bands) + assert result == sorted(result) + + +@pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) +def test_get_worm_components_partial_4_no_duplicates(num_bands): + result = get_worm_components_partial_4(num_bands) + assert len(result) == len(set(result)) + + +@pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) +def test_get_worm_components_partial_4_is_subset_of_full(num_bands): + assert set(get_worm_components_partial_4(num_bands)).issubset(set(get_worm_components_4(num_bands))) + + +def test_get_worm_components_partial_4_single_band_matches_full(): + assert get_worm_components_partial_4(1) == get_worm_components_4(1) + + +@pytest.mark.parametrize("num_bands", [2, 3, 4]) +def test_get_worm_components_partial_4_excludes_ijjj_type_orbitals(num_bands): + partial_indices = set(get_worm_components_partial_4(num_bands)) + spins = [[0, 0, 0, 0], [1, 1, 1, 1], [0, 0, 1, 1], [1, 1, 0, 0], [1, 0, 0, 1], [0, 1, 1, 0]] + for i, j in it.permutations(range(num_bands), 2): + for excluded_orb in [[i, j, j, j], [j, i, j, j], [j, j, i, j], [j, j, j, i]]: + for s in spins: + idx = int(component2index_general_4(num_bands, excluded_orb, s)) + assert idx not in partial_indices + + +@pytest.mark.parametrize("num_bands", [2, 3, 4]) +def test_get_worm_components_partial_4_strictly_smaller_than_full(num_bands): + assert len(get_worm_components_partial_4(num_bands)) < len(get_worm_components_4(num_bands)) + + +@pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) +def test_component2index_band_4_roundtrip(num_bands): + for orb in it.product(range(num_bands), repeat=4): + idx = component2index_band_4(num_bands, 4, list(orb)) + assert index2component_band_4(num_bands, 4, idx) == list(orb) + + +@pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) +def test_component2index_band_4_all_indices_unique(num_bands): + indices = [component2index_band_4(num_bands, 4, list(orb)) for orb in it.product(range(num_bands), repeat=4)] + assert len(indices) == len(set(indices)) + + +@pytest.mark.parametrize("num_bands", [1, 2, 3, 4]) +def test_component2index_band_4_index_range(num_bands): + indices = [component2index_band_4(num_bands, 4, list(orb)) for orb in it.product(range(num_bands), repeat=4)] + assert min(indices) == 1 + assert max(indices) == num_bands**4 + + +def test_component2index_band_4_single_band(): + assert component2index_band_4(1, 4, [0, 0, 0, 0]) == 1 + assert index2component_band_4(1, 4, 1) == [0, 0, 0, 0] + + +@pytest.mark.parametrize("num_bands", [1, 2, 3]) +def test_both_indexing_schemes_produce_correct_number_of_unique_indices(num_bands): + band_indices = {component2index_band_4(num_bands, 4, list(orb)) for orb in it.product(range(num_bands), repeat=4)} + spin_indices = { + component2index_general_4(num_bands, list(orb), [0, 0, 0, 0]) for orb in it.product(range(num_bands), repeat=4) + } + assert len(band_indices) == num_bands**4 + assert len(spin_indices) == num_bands**4 From c0ada977152ba5df9b41a98a206f96ee0e78fc4f Mon Sep 17 00:00:00 2001 From: julpe Date: Mon, 30 Mar 2026 08:17:45 +0200 Subject: [PATCH 2/2] Added anderson mixing alongside pulay mixing to try out other convergence. --- moldga/config.py | 1 + moldga/nonlocal_sde.py | 117 ++++++++++++++++++++++++++++++++++++++++- 2 files changed, 116 insertions(+), 2 deletions(-) diff --git a/moldga/config.py b/moldga/config.py index 1a0b0d8..deff208 100644 --- a/moldga/config.py +++ b/moldga/config.py @@ -88,6 +88,7 @@ def __init__(self): self.previous_sc_path: str = "./" self.use_lambda_correction: bool = False self.restrict_chi_phys: bool = False + self.anderson_prev_res: float | None = None class EliashbergConfig: diff --git a/moldga/nonlocal_sde.py b/moldga/nonlocal_sde.py index 83c269c..dc79b6d 100644 --- a/moldga/nonlocal_sde.py +++ b/moldga/nonlocal_sde.py @@ -536,7 +536,7 @@ def calculate_self_energy_q( sigma_old, starting_iter = get_starting_sigma(config.self_consistency.previous_sc_path, sigma_dmft) if starting_iter > 0: logger.info( - f"Using previous calculation and starting the self-consistency loop at iteration {starting_iter+1}." + f"Using previous calculation and starting the self-consistency loop at iteration {starting_iter + 1}." ) sigma_old = sigma_old.cut_niv(config.box.niw_core + config.box.niv_full + 10) @@ -771,9 +771,122 @@ def get_result(idx: int): 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}." + f"Pulay mixing applied with {n_hist} previous iterations and " + f"a mixing parameter of {config.self_consistency.mixing}." ) + return sigma_new + if ( + config.self_consistency.mixing_strategy == "anderson" + and current_iter > n_hist + and config.self_consistency.save_iter + and config.output.save_quantities + ): + last_results = read_last_n_sigmas_from_files( + n_hist, config.output.output_path, config.self_consistency.previous_sc_path + ) + + sigma_dmft_stacked = np.tile(sigma_dmft.mat, (config.lattice.k_grid.nk_tot, 1, 1, 1)) + last_proposals = [sigma_dmft_stacked] + last_results + last_results = last_results + [sigma_new.mat] + + niv = sigma_new.current_shape[-1] // 2 + niv_core = config.box.niv_core + + 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] + + shape = last_results[-1].shape + n_total = int(np.prod(shape)) + + r_matrix = np.zeros((2 * n_total, n_hist), dtype=np.float64) + f_matrix = np.zeros_like(r_matrix) + f_i = np.zeros((2 * n_total), dtype=np.float64) + + flat = lambda x: x.reshape(-1) + + # Build dx and df + for i in range(n_hist): + x_i = flat(last_proposals[-2 - i]) + x_ip1 = flat(last_proposals[-1 - i]) + + f_i_ = flat(last_results[-2 - i]) - x_i + f_ip1_ = flat(last_results[-1 - i]) - x_ip1 + + dx = x_ip1 - x_i + df = f_ip1_ - f_i_ + + r_matrix[:n_total, i] = dx.real + r_matrix[n_total:, i] = dx.imag + + f_matrix[:n_total, i] = df.real + f_matrix[n_total:, i] = df.imag + + # Current residual + f_curr = flat(last_results[-1]) - flat(last_proposals[-1]) + f_i[:n_total] = f_curr.real + f_i[n_total:] = f_curr.imag + + # SVD-based regularized solve + try: + u, s, vh = np.linalg.svd(f_matrix, full_matrices=False) + + # truncate small singular values + s_max = s[0] if len(s) > 0 else 1.0 + cutoff = 1e-10 * s_max + mask = s > cutoff + + if not np.any(mask): + raise np.linalg.LinAlgError("All singular values truncated") + + s_trunc = s[mask] + u_trunc = u[:, mask] + vh_trunc = vh[mask, :] + + # Tikhonov regularization + lam = 1e-8 + s_inv = s_trunc / (s_trunc**2 + lam) + coeffs = vh_trunc.T @ (s_inv * (u_trunc.T @ f_i)) + except np.linalg.LinAlgError: + logger.warning("Anderson SVD failed -> fallback to linear.") + return alpha * sigma_new + (1 - alpha) * sigma_old + + # 4) adaptive damping + norm_f = np.linalg.norm(f_i) + if config.self_consistency.anderson_prev_res is not None: + prev = config.self_consistency.anderson_prev_res + if norm_f < prev: + alpha_eff = min(1.2 * alpha, 1.0) + else: + alpha_eff = 0.5 * alpha + else: + alpha_eff = alpha + + config.self_consistency.anderson_prev_res = float(norm_f) + + # Anderson update + update = alpha_eff * f_i - (r_matrix + alpha_eff * f_matrix) @ coeffs + + # safety clipping + norm_u = np.linalg.norm(update) + if norm_f > 0 and norm_u > 5 * norm_f: + update *= 5 * norm_f / norm_u + + update = update[:n_total] + 1j * update[n_total:] + + candidate = last_proposals[-1] + update.reshape(shape) + + # fallback if step is bad + new_residual = np.linalg.norm(flat(last_results[-1]) - flat(candidate)) + if new_residual > norm_f: + logger.warning("Anderson step rejected -> fallback to linear.") + return alpha * sigma_new + (1 - alpha) * sigma_old + + # Update the new self-energy + sigma_new.mat[..., niv - niv_core : niv + niv_core] = candidate + + logger.info(f"Anderson mixing (robust) applied (m={n_hist}, alpha_eff={alpha_eff:.3f}).") + return sigma_new sigma_new = config.self_consistency.mixing * sigma_new + (1 - config.self_consistency.mixing) * sigma_old