From ff820b9775e165675049885b8409efe563ba98e7 Mon Sep 17 00:00:00 2001 From: alberto-carta <49793269+alberto-carta@users.noreply.github.com> Date: Wed, 27 Sep 2023 19:46:58 +0200 Subject: [PATCH 1/2] Added all changes needed to broyden mix sigma + deprecation fix in mixer --- python/solid_dmft/csc_flow.py | 13 + python/solid_dmft/dmft_cycle.py | 43 ++- .../dmft_tools/greens_functions_mixer.py | 325 +++++++++++++++++- python/solid_dmft/read_config.py | 31 +- 4 files changed, 387 insertions(+), 25 deletions(-) diff --git a/python/solid_dmft/csc_flow.py b/python/solid_dmft/csc_flow.py index 507b7090..b4850fdd 100644 --- a/python/solid_dmft/csc_flow.py +++ b/python/solid_dmft/csc_flow.py @@ -256,10 +256,13 @@ def csc_flow_control(general_params, solver_params, dft_params, advanced_params) # Reads in iteration offset if restarting iteration_offset = 0 + iter_csc = 0 if mpi.is_master_node() and os.path.isfile(general_params['seedname']+'.h5'): with HDFArchive(general_params['seedname']+'.h5', 'r') as archive: if 'DMFT_results' in archive and 'iteration_count' in archive['DMFT_results']: iteration_offset = archive['DMFT_results']['iteration_count'] + if 'DMFT_results' in archive and 'iter_csc' in archive['DMFT_results']: + iter_csc = archive['DMFT_results']['iter_csc'] iteration_offset = mpi.bcast(iteration_offset) iter_dmft = iteration_offset+1 @@ -333,6 +336,15 @@ def csc_flow_control(general_params, solver_params, dft_params, advanced_params) print('DMFT cycle took {:10.3f} seconds'.format(end_time_dmft-start_time_dmft)) print('='*80 + '\n') + + iter_csc += 1 + + if mpi.is_master_node() and os.path.isfile(general_params['seedname']+'.h5'): + + with HDFArchive(general_params['seedname']+'.h5', 'a') as archive: + if not 'iter_csc' in archive['DMFT_results']: + archive['DMFT_results'].create_group('iter_csc') + archive['DMFT_results']['iter_csc'] = iter_csc # If all steps are executed or calculation is converged, finish DFT+DMFT loop if is_converged or iter_dmft > general_params['n_iter_dmft'] + iteration_offset: break @@ -340,6 +352,7 @@ def csc_flow_control(general_params, solver_params, dft_params, advanced_params) # Restarts DFT mpi.barrier() start_time_dft = timer() + mpi.report(f' solid_dmft: Starting CSC iteration {iter_csc}') mpi.report(' solid_dmft: Running {}...'.format(dft_params['dft_code'].upper())) # Runs DFT and converter diff --git a/python/solid_dmft/dmft_cycle.py b/python/solid_dmft/dmft_cycle.py index f8c11ad0..ee542431 100755 --- a/python/solid_dmft/dmft_cycle.py +++ b/python/solid_dmft/dmft_cycle.py @@ -56,6 +56,7 @@ from solid_dmft.dmft_tools import manipulate_chemical_potential as manipulate_mu from solid_dmft.dmft_tools import initial_self_energies as initial_sigma from solid_dmft.dmft_tools import greens_functions_mixer as gf_mixer +from collections import deque def _determine_block_structure(sum_k, general_params, advanced_params): @@ -342,10 +343,27 @@ def dmft_cycle(general_params, solver_params, advanced_params, dft_params, E_kin_dft = None # check for previous broyden data oterhwise initialize it: - if mpi.is_master_node() and general_params['g0_mix_type'] == 'broyden': + #if mpi.is_master_node() and general_params['g0_mix_type'] == 'broyden': + # if not 'broyler' in archive['DMFT_results']: + # archive['DMFT_results']['broyler'] = [{'mu' : [],'V': [], 'dV': [], 'F': [], 'dF': []} + # for _ in range(sum_k.n_inequiv_shells)] + if mpi.is_master_node() and general_params['mix_type']=='broyden': if not 'broyler' in archive['DMFT_results']: - archive['DMFT_results']['broyler'] = [{'mu' : [],'V': [], 'dV': [], 'F': [], 'dF': []} - for _ in range(sum_k.n_inequiv_shells)] + #archive['DMFT_results']['broyler'] = [{'mu' : [],'V': [], 'dV': [], 'F': [], 'dF': []} + # for _ in range(sum_k.n_inequiv_shells)] + + broyler_list = [{'mu' : [],'V': [], 'dV': [], 'F': [], 'dF': [], + 'broy_max_it': general_params['broy_max_it'], + 'deg_shell':sum_k.deg_shells[_icrsh1] } + for _icrsh1 in range(sum_k.n_inequiv_shells)] + archive['DMFT_results']['broyler'] = broyler_list + + mpi.report("Testing mixing data initialization") + broyler_obj = gf_mixer.BroylerClass(broyler_list) + mpi.report(broyler_obj.__dict__) + mpi.report("Comparing with broyler list") + mpi.report(broyler_list) + # Generates a rotation matrix to change the basis if general_params['set_rot'] != 'none': @@ -631,8 +649,12 @@ def _dmft_step(sum_k, solvers, it, general_params, # mixing of G0 if wanted from the second iteration on if it > 1: - solvers[icrsh] = gf_mixer.mix_g0(solvers[icrsh], general_params, icrsh, archive, - G0_freq_previous[icrsh], it, sum_k.deg_shells[icrsh]) + if 'G0'in general_params['mix_quantity']: + mpi.report(f"XXXXXXX Calling {general_params['mix_type']} mixing on {general_params['mix_quantity']}") + solvers[icrsh] = gf_mixer.mix_general(general_params, icrsh, solvers[icrsh], G0_freq_previous[icrsh], it = it, F_type='G0', + deg_shell=sum_k.deg_shells[icrsh], archive=archive) + # solvers[icrsh] = gf_mixer.mix_g0(solvers[icrsh], general_params, icrsh, archive, + # G0_freq_previous[icrsh], it, sum_k.deg_shells[icrsh]) if general_params['solver_type'] in ['cthyb', 'ctint', 'hubbardI', 'inchworm']: solvers[icrsh].G0_freq << make_hermitian(solvers[icrsh].G0_freq) @@ -690,7 +712,16 @@ def _dmft_step(sum_k, solvers, it, general_params, # if CPA average Sigma over impurities before mixing if general_params['dc'] and general_params['dc_type'] == 4: solvers = gf_mixer.mix_cpa(cpa_G0_freq, sum_k.n_inequiv_shells, solvers) - solvers = gf_mixer.mix_sigma(general_params, sum_k.n_inequiv_shells, solvers, Sigma_freq_previous) + #solvers = gf_mixer.mix_sigma(general_params, sum_k.n_inequiv_shells, solvers, Sigma_freq_previous) + + if 'Sigma'in general_params['mix_quantity']: + if it > 1: + mpi.report(f"XXXXXXX Calling {general_params['mix_type']} mixing on {general_params['mix_quantity']}") + for icrsh in range(sum_k.n_inequiv_shells): + solvers[icrsh] = gf_mixer.mix_general(general_params, icrsh, solvers[icrsh], Sigma_freq_previous[icrsh], F_type='Sigma', it=it, + deg_shell=sum_k.deg_shells[icrsh], archive=archive) + + # calculate new DC # for the hartree solver the DC potential will be formally set to zero as it is already present in the Sigma diff --git a/python/solid_dmft/dmft_tools/greens_functions_mixer.py b/python/solid_dmft/dmft_tools/greens_functions_mixer.py index 2ac2ca25..d68423e7 100644 --- a/python/solid_dmft/dmft_tools/greens_functions_mixer.py +++ b/python/solid_dmft/dmft_tools/greens_functions_mixer.py @@ -33,7 +33,7 @@ from triqs.gf.tools import inverse from . import legendre_filter - +from collections import deque def _bgf_to_vec(bgf, deg_shells=None): """ @@ -44,9 +44,9 @@ def _bgf_to_vec(bgf, deg_shells=None): deg_vecs = [] for deg_block in deg_shells: deg_vecs.append(bgf[deg_block[0]]) - vec = np.hstack(gf.data.flatten() for gf in deg_vecs) + vec = np.hstack([gf.data.flatten() for gf in deg_vecs]) else: - vec = np.hstack(bgf[bl].data.flatten() for bl in bgf.indices) + vec = np.hstack([bgf[bl].data.flatten() for bl in bgf.indices]) return vec @@ -72,6 +72,117 @@ def _vec_to_bgf(vec, bgf, deg_shells=None): start = end return G +def _broyden_update_general(it, broyler, general_params, deg_shells, F_freq, F_freq_previous): + """ + calculates the broyden update of G0 using the algorithm presented in: + doi.org/10.1103/PhysRevB.80.125125 (Rok Zitko) & + doi.org/10.1103/PhysRevB.38.12807 (D.D. Johnson) + + optimize function to find roots: F(G0)=dmft_step(G0) - G0 = 0 + where dmft_step gives back the normal defined G0=inverse(Gloc^-1 - Sigma) + + Parameters + ---------- + it : int + iteration number + broyler: dict + dict containing the broyden variables + general_paramters: general params dict + + Returns + ------- + """ + #alpha = general_params['g0_mix'] + alpha = general_params['mix_greed'] + mpi.report(f"Calling broyden step") + # hard code w_0 + w_0 = 0.01 + + F_broyden_update = F_freq.copy() + + # build here leg compression + if type(F_freq.mesh) is MeshImFreq and 'n_l' in general_params: + F = legendre_filter.apply(make_gf_from_fourier(F_freq), order=general_params['n_l'] ) + F_prev = legendre_filter.apply(make_gf_from_fourier(F_freq_previous), order=general_params['n_l'] ) + else: + F = F_freq.copy() + F_prev = F_freq_previous.copy() + + # if this is the first time broyden update is called (it==2) + # store init F first + if it == 2: + # V holds the F as 1d numpy vectors + broyler['V'].append(_bgf_to_vec(F_prev, deg_shells)) + # F = dmft_step_F - V[-1] + broyler['F'].append(_bgf_to_vec(F, deg_shells)-broyler['V'][-1]) + # for the second iteration the broyden vec corresponds to simple linear mixing with alpha + broyden_vec = np.zeros(broyler['V'][0].shape[0],dtype='complex') + + # for all other iteration calc V_mp1 = V[-1] + alpha * F[-1] - broyden_vec + # broyden vec is sum(n=1 to it-1) gamma_m_n U^n + # gamma_m_n = sum(k=1 to it-1) c_k^m beta_k_n^m + # gamma_m_n is a number for each iteration and U^n is of dim F for each iteration + else: + # get size of broyden update vector + dim = broyler['V'][0].shape[0] + broyden_vec = np.zeros(dim,dtype='complex') + + # define m as it-1-1 (last minus for idx starting at 0) + m = it-2 + + # update broyden input params + broyler['F'].append(_bgf_to_vec(F, deg_shells)-broyler['V'][-1]) + norm_F = np.linalg.norm(broyler['F'][-1]-broyler['F'][-2]) + print('broyden norm: {:.3f}'.format(norm_F)) + # dF is normed differences between last two functions + broyler['dF'].append( (broyler['F'][-1]-broyler['F'][-2])/norm_F ) + # dV the normed difference between the two last input F + broyler['dV'].append( (broyler['V'][-1] - broyler['V'][-2])/norm_F ) + + # now build broyden correction vector + + #only keep last m iterations + mem_iter = general_params['broy_max_it']-1 + start = 0 + mat_dim = m + if m > mem_iter and not general_params['broy_max_it'] == -1: + for n in range(0, m-mem_iter): + start += 1 + mat_dim -= 1 + + # first create beta_m + A = np.zeros((mat_dim,mat_dim),dtype='complex') + gamma = np.zeros(mat_dim,dtype='complex') + for k_i in range(start,m): + for n_i in range(start,m): + A[k_i-start,n_i-start] = np.dot(broyler['dF'][n_i].conjugate().transpose(),broyler['dF'][k_i]) + if n_i == k_i: + A[k_i-start,n_i-start] += w_0**2 + beta_m = np.linalg.inv(A) + for n_i in range(start,m): + U_n = alpha * broyler['dF'][n_i] + broyler['dV'][n_i] + + # construct gamma + for k_i in range(start,m): + c_k = np.dot(broyler['dF'][k_i].conjugate().transpose(),broyler['F'][-1]) + gamma[n_i-start] += c_k*beta_m[k_i-start,n_i-start] + broyden_vec += gamma[n_i-start] * U_n + + # update input F + V_mp1 = broyler['V'][-1] + alpha * broyler['F'][-1] - broyden_vec + broyler['V'].append(V_mp1) + + # convert Gl back to F_freq + if type(F_freq.mesh) is MeshImFreq and 'n_l' in general_params: + F_l_update = _vec_to_bgf(broyler['V'][-1], F, deg_shells) + for block, g in F_l_update: + g.enforce_discontinuity(np.identity(g.target_shape[0])) + F_broyden_update[block].set_from_legendre(g) + + else: + F_broyden_update << _vec_to_bgf(broyler['V'][-1], F, deg_shells) + + return F_broyden_update, broyler def _broyden_update(it, broyler, general_params, deg_shells, G0_freq, G0_freq_previous): """ @@ -207,20 +318,214 @@ def mix_g0(solver, general_params, icrsh, archive, G0_freq_previous, it, deg_she return solver +def mix_general(general_params, icrsh, solver, F_freq_previous, F_type, deg_shell, it, archive): + mix_greed = general_params['mix_greed'] + mix_type = general_params['mix_type'] + + if F_type not in ['Gimp', 'G0', 'Sigma']: + raise ValueError(f"You are trying to mix the quantity '{F_type}', allowed values are: 'Gimp', 'G0', 'Sigma'") + + def _load_F_solver(): + if F_type == 'G0': + F_freq = solver.G0_freq + elif F_type == 'Gimp': + F_freq = solver.G_freq + elif F_type == 'Sigma': + F_freq = solver.Sigma_freq + + return F_freq + + def _save_to_solver(): + if F_type == 'G0': + solver.G0_freq << mpi.bcast(solver.G0_freq) + elif F_type == 'Gimp': + solver.G_freq << mpi.bcast(solver.G_freq) + elif F_type == 'Sigma': + solver.Sigma_freq << mpi.bcast(solver.Sigma_freq) + + + + if mix_greed < 1.0 and mpi.is_master_node(): + if mix_type == 'linear': + print(f'Linear mixing {F_type} with previous iteration by factor {mix_greed:.3f}\n') + F_freq = _load_F_solver() + F_freq << (mix_greed * F_freq + (1-mix_greed) * F_freq_previous) + + elif mix_type == 'broyden': + print(f'Broyden mixing {F_type} with previous iteration by factor {mix_greed:.3f}\n') + mpi.report('\n############\n!!!! WARNING !!!! broyden mixing is still in early testing stage ! Use with caution.\n############\n') + # TODO implement broyden mixing for Sigma + + broyler = archive['DMFT_results']['broyler'] + F_freq = _load_F_solver() + + BroylerObj = BroylerClass(broyler) + F_update = BroylerObj.broyden_step(it, icrsh, general_params, F_freq, F_freq_previous) + + F_freq << F_update + archive['DMFT_results']['broyler'] = BroylerObj.broyler_to_dict() + + + + + # IMPORTANT!!!: the broadcast must be called outside of the MPI master node to avoid deadlocks! + _save_to_solver() + + + return solver + +class BroylerClass: + + def __init__(self, broyler_list): + self.n_inequiv_shells = len(broyler_list) + self.broy_max_it = broyler_list[0]['broy_max_it'] + self.deg_shells = [broyler_list[icrsh]['deg_shell'] for icrsh in range(self.n_inequiv_shells)] + self.mu = [deque( broyler_list[icrsh]['mu'], maxlen = self.broy_max_it ) for icrsh in range(self.n_inequiv_shells)] + self.V = [deque( broyler_list[icrsh]['V'] , maxlen = self.broy_max_it ) for icrsh in range(self.n_inequiv_shells)] + self.dV = [deque( broyler_list[icrsh]['dV'], maxlen = self.broy_max_it ) for icrsh in range(self.n_inequiv_shells)] + self.F = [deque( broyler_list[icrsh]['F'] , maxlen = self.broy_max_it ) for icrsh in range(self.n_inequiv_shells)] + self.dF = [deque( broyler_list[icrsh]['dF'], maxlen = self.broy_max_it ) for icrsh in range(self.n_inequiv_shells)] + + def broyler_to_dict(self): + broyler_list_output = [] + for icrsh in range(self.n_inequiv_shells): + broyler_list_output.append({ + 'mu': list(self.mu[icrsh]), + 'V': list(self.V[icrsh]), + 'dV': list(self.dV[icrsh]), + 'F': list(self.F[icrsh]), + 'dF': list(self.dF[icrsh]), + 'broy_max_it': self.broy_max_it, + 'deg_shell': self.deg_shells[icrsh], + }) + + return broyler_list_output + + def broyden_step(self, it, icrsh, general_params, F_freq, F_freq_previous): + """ + calculates the broyden update of G0 using the algorithm presented in: + doi.org/10.1103/PhysRevB.80.125125 (Rok Zitko) & + doi.org/10.1103/PhysRevB.38.12807 (D.D. Johnson) + + optimize function to find roots: F(G0)=dmft_step(G0) - G0 = 0 + where dmft_step gives back the normal defined G0=inverse(Gloc^-1 - Sigma) + + Parameters + ---------- + it : int + iteration number + broyler: dict + dict containing the broyden variables + general_paramters: general params dict + + Returns + ------- + """ + alpha = general_params['mix_greed'] + deg_shells = self.deg_shells[icrsh] + + + mpi.report(f"Object method: broyden step") + # hard code w_0 + w_0 = 0.01 + + F_broyden_update = F_freq.copy() + + # build here leg compression + if type(F_freq.mesh) is MeshImFreq and 'n_l' in general_params: + F = legendre_filter.apply(make_gf_from_fourier(F_freq), order=general_params['n_l'] ) + F_prev = legendre_filter.apply(make_gf_from_fourier(F_freq_previous), order=general_params['n_l'] ) + else: + F = F_freq.copy() + F_prev = F_freq_previous.copy() + + # if this is the first time broyden update is called (it==2) + # store init F first + if it == 2: + # V holds the F as 1d numpy vectors + self.V[icrsh].append(_bgf_to_vec(F_prev, deg_shells)) + # F = dmft_step_F - V[-1] + self.F[icrsh].append(_bgf_to_vec(F, deg_shells)-self.V[icrsh][-1]) + # for the second iteration the broyden vec corresponds to simple linear mixing with alpha + broyden_vec = np.zeros(self.V[icrsh][0].shape[0],dtype='complex') + + # for all other iteration calc V_mp1 = V[-1] + alpha * F[-1] - broyden_vec + # broyden vec is sum(n=1 to it-1) gamma_m_n U^n + # gamma_m_n = sum(k=1 to it-1) c_k^m beta_k_n^m + # gamma_m_n is a number for each iteration and U^n is of dim F for each iteration + else: + # get size of broyden update vector + dim = self.V[icrsh][0].shape[0] + broyden_vec = np.zeros(dim,dtype='complex') + + # define m as it-1-1 (last minus for idx starting at 0) + m = it-2 + + # update broyden input params + self.F[icrsh].append(_bgf_to_vec(F, deg_shells)-self.V[icrsh][-1]) + norm_F = np.linalg.norm(self.F[icrsh][-1]-self.F[icrsh][-2]) + print('broyden norm: {:.3f}'.format(norm_F)) + # dF is normed differences between last two functions + self.dF[icrsh].append( (self.F[icrsh][-1]-self.F[icrsh][-2])/norm_F ) + # dV the normed difference between the two last input F + self.dV[icrsh].append( (self.V[icrsh][-1] - self.V[icrsh][-2])/norm_F ) + + # now build broyden correction vector + + #only keep last m iterations + mem_iter = self.broy_max_it-1 + start = 0 + mat_dim = min(self.broy_max_it, m) + + # first create beta_m + A = np.zeros((mat_dim,mat_dim),dtype='complex') + gamma = np.zeros(mat_dim,dtype='complex') + for k_i in range(start,mat_dim): + for n_i in range(start,mat_dim): + A[k_i-start,n_i-start] = np.dot(self.dF[icrsh][n_i].conjugate().transpose(),self.dF[icrsh][k_i]) + if n_i == k_i: + A[k_i-start,n_i-start] += w_0**2 + beta_m = np.linalg.inv(A) + for n_i in range(start,mat_dim): + U_n = alpha * self.dF[icrsh][n_i] + self.dV[icrsh][n_i] + + # construct gamma + for k_i in range(start,mat_dim): + c_k = np.dot(self.dF[icrsh][k_i].conjugate().transpose(),self.F[icrsh][-1]) + gamma[n_i-start] += c_k*beta_m[k_i-start,n_i-start] + broyden_vec += gamma[n_i-start] * U_n + + # update input F + V_mp1 = self.V[icrsh][-1] + alpha * self.F[icrsh][-1] - broyden_vec + self.V[icrsh].append(V_mp1) + + # convert Gl back to F_freq + if type(F_freq.mesh) is MeshImFreq and 'n_l' in general_params: + F_l_update = _vec_to_bgf(self.V[icrsh][-1], F, deg_shells) + for block, g in F_l_update: + g.enforce_discontinuity(np.identity(g.target_shape[0])) + F_broyden_update[block].set_from_legendre(g) + + else: + F_broyden_update << _vec_to_bgf(self.V[icrsh][-1], F, deg_shells) + + return F_broyden_update + + + + + def mix_sigma(general_params, n_inequiv_shells, solvers, Sigma_freq_previous): - if general_params['sigma_mix'] < 1.0 and mpi.is_master_node(): - print('mixing sigma with previous iteration by factor {:.3f}\n'.format(general_params['sigma_mix'])) + if general_params['mix_greed'] < 1.0 and mpi.is_master_node(): + print('mixing sigma with previous iteration by factor {:.3f}\n'.format(general_params['mix_greed'])) for icrsh in range(n_inequiv_shells): - solvers[icrsh].Sigma_freq << (general_params['sigma_mix'] * solvers[icrsh].Sigma_freq - + (1-general_params['sigma_mix']) * Sigma_freq_previous[icrsh]) - + solvers[icrsh].Sigma_freq << (general_params['mix_greed'] * solvers[icrsh].Sigma_freq + + (1-general_params['mix_greed']) * Sigma_freq_previous[icrsh]) for icrsh in range(n_inequiv_shells): solvers[icrsh].Sigma_freq << mpi.bcast(solvers[icrsh].Sigma_freq) return solvers - - def init_cpa(sum_k, solvers, general_params): # initialize cpa_G_time, cpa_G0_freq; extract cpa_G_loc cpa_G_time = sum_k.block_structure.create_gf(ish=0, gf_function=Gf, space='solver', diff --git a/python/solid_dmft/read_config.py b/python/solid_dmft/read_config.py index a4d88e89..ea9e35b3 100755 --- a/python/solid_dmft/read_config.py +++ b/python/solid_dmft/read_config.py @@ -435,6 +435,13 @@ # Workaround to get the default configparser boolean converter BOOL_PARSER = lambda b: ConfigParser()._convert_to_boolean(b) +# Converts to an integer if possible, implemented to keep legacy mode in DC_type +def TRY_INT_PARSER(s): + try: + result = int(s) + except: + result = str(s) + return result # TODO: it might be nicer to not have optional parameters at all and instead use # explicit default values @@ -485,7 +492,7 @@ 'dc': {'converter': BOOL_PARSER, 'used': True, 'default': True}, - 'dc_type': {'converter': int, 'valid for': lambda x, _: x in (0, 1, 2, 3, 4), + 'dc_type': {'converter': TRY_INT_PARSER, 'valid for': lambda x, _: x in [0, 1, 2, 3, 4, 'cFLL','sFLL', 'cAMF', 'sAMF', 'cHeld'], 'used': lambda params: params['general']['dc']}, 'prec_mu': {'converter': float, 'valid for': lambda x, _: x > 0, 'used': True}, @@ -557,7 +564,20 @@ 'afm_order': {'converter': BOOL_PARSER, 'used': lambda params: params['general']['magnetic'], 'default': False}, - + + 'mix_quantity': {'converter': lambda s: list(map(str, s.split(','))), + 'used': True, + 'default': []}, + + 'mix_type': {'valid for': lambda x, _: x in ('linear', 'broyden'), + 'used': True, 'default': 'linear'}, + + 'mix_greed': {'converter': float, 'valid for': lambda x, _: x > 0, + 'used': True, 'default': 1.0}, + + 'broy_max_it': {'converter': int, 'valid for': lambda x, _: x >= 1 or x==-1 , + 'used': True, 'default': 2}, + 'sigma_mix': {'converter': float, 'valid for': lambda x, params: x >= 0 and (np.isclose(params['general']['g0_mix'], 1) or np.isclose(x, 1)), @@ -569,8 +589,6 @@ 'g0_mix_type': {'valid for': lambda x, _: x in ('linear', 'broyden'), 'used': True, 'default': 'linear'}, - 'broy_max_it': {'converter': int, 'valid for': lambda x, _: x >= 1 or x==-1 , - 'used': lambda params: params['general']['g0_mix_type'] == 'broyden', 'default': -1}, 'calc_energies': {'converter': BOOL_PARSER, 'used': True, 'default': False}, @@ -1229,11 +1247,6 @@ def read_config(config_file): raise ValueError('The following parameters are not valid:' + invalid_error_string) - # warning if sigma mixing is used, remove in future versions - if parameters['general']['sigma_mix'] < 1.0: - if parameters['general']['g0_mix'] < 1.0: - raise ValueError('You shall not use Sigma and G0 mixing together!') - # Workarounds for some parameters # make sure that pick_solver_struct and map_solver_struct are a list of dict From aaffb03925b778e67c8a5aa096eec5a03684a3ae Mon Sep 17 00:00:00 2001 From: alb-carta <49793269+alberto-carta@users.noreply.github.com> Date: Thu, 30 Nov 2023 11:37:15 +0100 Subject: [PATCH 2/2] added possibility of updating mu at every impurity iteration --- python/solid_dmft/dmft_cycle.py | 17 +++++++++++++++-- python/solid_dmft/read_config.py | 23 ++++++++++++----------- 2 files changed, 27 insertions(+), 13 deletions(-) diff --git a/python/solid_dmft/dmft_cycle.py b/python/solid_dmft/dmft_cycle.py index ee542431..ec4e74b8 100755 --- a/python/solid_dmft/dmft_cycle.py +++ b/python/solid_dmft/dmft_cycle.py @@ -628,6 +628,11 @@ def _dmft_step(sum_k, solvers, it, general_params, for icrsh in range(sum_k.n_inequiv_shells): # copy the block of G_loc into the corresponding instance of the impurity solver # TODO: why do we set solvers.G_freq? Isn't that simply an output of the solver? + ### + + if general_params['update_mu_each_imp']: + G_loc_all = sum_k.extract_G_loc(iw_or_w='iw') + solvers[icrsh].G_freq << G_loc_all[icrsh] density_shell_pre[icrsh] = np.real(solvers[icrsh].G_freq.total_density()) @@ -687,6 +692,14 @@ def _dmft_step(sum_k, solvers, it, general_params, mpi.report('Actual time for solver: {:.2f} s'.format(timer() - start_time)) # some printout of the obtained density matrices and some basic checks from the unsymmetrized solver output + + ###new + if general_params['update_mu_each_imp']: + sum_k.put_Sigma([solvers[i].Sigma_freq for i in range(sum_k.n_inequiv_shells)]) + sum_k = manipulate_mu.update_mu(general_params, sum_k, it, archive) + + + density_shell[icrsh] = np.real(solvers[icrsh].G_freq_unsym.total_density()) density_tot += density_shell[icrsh]*shell_multiplicity[icrsh] density_mat_unsym[icrsh] = solvers[icrsh].G_freq_unsym.density() @@ -715,10 +728,10 @@ def _dmft_step(sum_k, solvers, it, general_params, #solvers = gf_mixer.mix_sigma(general_params, sum_k.n_inequiv_shells, solvers, Sigma_freq_previous) if 'Sigma'in general_params['mix_quantity']: - if it > 1: + if it > 0: mpi.report(f"XXXXXXX Calling {general_params['mix_type']} mixing on {general_params['mix_quantity']}") for icrsh in range(sum_k.n_inequiv_shells): - solvers[icrsh] = gf_mixer.mix_general(general_params, icrsh, solvers[icrsh], Sigma_freq_previous[icrsh], F_type='Sigma', it=it, + solvers[icrsh] = gf_mixer.mix_general(general_params, icrsh, solvers[icrsh], Sigma_freq_previous[icrsh], F_type='Sigma', it=it+1, deg_shell=sum_k.deg_shells[icrsh], archive=archive) diff --git a/python/solid_dmft/read_config.py b/python/solid_dmft/read_config.py index ea9e35b3..d3ddc5fe 100755 --- a/python/solid_dmft/read_config.py +++ b/python/solid_dmft/read_config.py @@ -564,7 +564,10 @@ def TRY_INT_PARSER(s): 'afm_order': {'converter': BOOL_PARSER, 'used': lambda params: params['general']['magnetic'], 'default': False}, - + + 'update_mu_each_imp': {'converter': BOOL_PARSER, + 'used': True, 'default': False}, + 'mix_quantity': {'converter': lambda s: list(map(str, s.split(','))), 'used': True, 'default': []}, @@ -578,16 +581,14 @@ def TRY_INT_PARSER(s): 'broy_max_it': {'converter': int, 'valid for': lambda x, _: x >= 1 or x==-1 , 'used': True, 'default': 2}, - 'sigma_mix': {'converter': float, - 'valid for': lambda x, params: x >= 0 and (np.isclose(params['general']['g0_mix'], 1) - or np.isclose(x, 1)), - 'used': True, 'default': 1.0}, - - 'g0_mix': {'converter': float, 'valid for': lambda x, _: x >= 0, - 'used': True, 'default': 1.0}, - - 'g0_mix_type': {'valid for': lambda x, _: x in ('linear', 'broyden'), - 'used': True, 'default': 'linear'}, + #'sigma_mix': {'converter': float, + # 'valid for': lambda x, params: x > 0 and (np.isclose(params['general']['g0_mix'], 1) + # or np.isclose(x, 1)), + # 'used': True, 'default': 1.0}, + #'g0_mix': {'converter': float, 'valid for': lambda x, _: x > 0, + # 'used': True, 'default': 1.0}, + #'g0_mix_type': {'valid for': lambda x, _: x in ('linear', 'broyden'), + # 'used': True, 'default': 'linear'}, 'calc_energies': {'converter': BOOL_PARSER, 'used': True, 'default': False},