From 4c8c19327baa5dc01f77215d4a40ae53b5b21c63 Mon Sep 17 00:00:00 2001 From: Henri Menke Date: Fri, 10 Nov 2023 15:40:57 +0100 Subject: [PATCH] Prototype for w2dyn_cthyb --- .../dmft_tools/results_to_archive.py | 3 + python/solid_dmft/dmft_tools/solver.py | 146 +++++++++++++++++- python/solid_dmft/read_config.py | 62 ++++++-- 3 files changed, 200 insertions(+), 11 deletions(-) diff --git a/python/solid_dmft/dmft_tools/results_to_archive.py b/python/solid_dmft/dmft_tools/results_to_archive.py index c89cc5d3..f69c44f7 100644 --- a/python/solid_dmft/dmft_tools/results_to_archive.py +++ b/python/solid_dmft/dmft_tools/results_to_archive.py @@ -98,6 +98,9 @@ def _compile_information(sum_k, general_params, solver_params, solvers, write_to_h5['F_freq_{}'.format(icrsh)] = solvers[icrsh].F_freq write_to_h5['F_time_{}'.format(icrsh)] = solvers[icrsh].F_time + if general_params['solver_type'] == 'w2dyn_cthyb': + pass # TODO + return write_to_h5 def write(archive, sum_k, general_params, solver_params, solvers, it, is_sampling, diff --git a/python/solid_dmft/dmft_tools/solver.py b/python/solid_dmft/dmft_tools/solver.py index fdfd362b..a2152a0b 100755 --- a/python/solid_dmft/dmft_tools/solver.py +++ b/python/solid_dmft/dmft_tools/solver.py @@ -231,6 +231,16 @@ def __init__(self, general_params, solver_params, advanced_params, sum_k, icrsh, self.git_hash = triqs_ctseg_hash self.version = version + elif self.general_params['solver_type'] == 'w2dyn_cthyb': + from w2dyn_cthyb.version import w2dyn_cthyb_hash, version + + # sets up necessary GF objects on ImFreq + self._init_ImFreq_objects() + # sets up solver + self.triqs_solver = self._create_w2dyn_cthyb_solver() + self.git_hash = w2dyn_cthyb_hash + self.version = version + # ******************************************************************** # initialize Freq and Time objects # ******************************************************************** @@ -264,7 +274,8 @@ def _init_ImFreq_objects(self): or self.general_params['solver_type'] == 'cthyb' and self.general_params['legendre_fit'] or self.general_params['solver_type'] == 'ctseg' and self.solver_params['measure_gl'] or self.general_params['solver_type'] == 'ctseg' and self.general_params['legendre_fit'] - or self.general_params['solver_type'] == 'hubbardI' and self.solver_params['measure_G_l']): + or self.general_params['solver_type'] == 'hubbardI' and self.solver_params['measure_G_l'] + or self.general_params['solver_type'] == 'w2dyn_cthyb' and self.solver_params['measure_G_l']): self.n_l = self.general_params['n_l'] self.G_l = self.sum_k.block_structure.create_gf(ish=self.icrsh, gf_function=Gf, space='solver', @@ -610,6 +621,50 @@ def make_positive_definite(G): # call postprocessing self._ctseg_postprocessing() + if self.general_params['solver_type'] == 'w2dyn_cthyb': + + # TODO: factor this block out of cthyb and w2dyn_cthyb + if self.general_params['cthyb_delta_interface']: + mpi.report('\n Using the delta interface for cthyb passing Delta(tau) and Hloc0 directly.') + # prepare solver input + sumk_eal = self.sum_k.eff_atomic_levels()[self.icrsh] + solver_eal = self.sum_k.block_structure.convert_matrix(sumk_eal, space_from='sumk', ish_from=self.sum_k.inequiv_to_corr[self.icrsh]) + # fill Delta_time from Delta_freq sum_k to solver + for name, g0 in self.G0_freq: + self.Delta_freq[name] << iOmega_n - inverse(g0) - solver_eal[name] + known_moments = make_zero_tail(self.Delta_freq[name], 1) + tail, err = fit_hermitian_tail(self.Delta_freq[name], known_moments) + # without SOC delta_tau needs to be real + if not self.sum_k.SO == 1: + self.triqs_solver.Delta_tau[name] << make_gf_from_fourier(self.Delta_freq[name], self.triqs_solver.Delta_tau.mesh, tail).real + else: + self.triqs_solver.Delta_tau[name] << make_gf_from_fourier(self.Delta_freq[name], self.triqs_solver.Delta_tau.mesh, tail) + + + # Make non-interacting operator for Hloc0 + Hloc_0 = Operator() + for spin, spin_block in solver_eal.items(): + for o1 in range(spin_block.shape[0]): + for o2 in range(spin_block.shape[1]): + # check if off-diag element is larger than threshold + if o1 != o2 and abs(spin_block[o1,o2]) < self.solver_params['off_diag_threshold']: + continue + else: + # TODO: adapt for SOC calculations, which should keep the imag part + Hloc_0 += spin_block[o1,o2].real/2 * (c_dag(spin,o1) * c(spin,o2) + c_dag(spin,o2) * c(spin,o1)) + self.solver_params['h_loc0'] = Hloc_0 + else: + # fill G0_freq from sum_k to solver + self.triqs_solver.G0_iw << self.G0_freq + + # Solve the impurity problem for icrsh shell + # ************************************* + self.triqs_solver.solve(h_int=self.h_int, **{ **self.solver_params, **random_seed }) + # ************************************* + + # call postprocessing + self._w2dyn_cthyb_postprocessing() + return # ******************************************************************** @@ -798,6 +853,24 @@ def _create_ctseg_solver(self): return triqs_solver + def _create_w2dyn_cthyb_solver(self): + r''' + Initialize w2dyn_cthyb solver instance + ''' + from triqs_w2dyn_cthyb import Solver as w2dyn_cthyb_solver + + gf_struct = self.sum_k.gf_struct_solver_list[self.icrsh] + # Construct the triqs_solver instances + if self.solver_params['measure_G_l']: + triqs_solver = w2dyn_cthyb_solver(beta=self.general_params['beta'], gf_struct=gf_struct, + n_iw=self.general_params['n_iw'], n_tau=self.general_params['n_tau'], + n_l=self.general_params['n_l']) + else: + triqs_solver = w2dyn_cthyb_solver(beta=self.general_params['beta'], gf_struct=gf_struct, + n_iw=self.general_params['n_iw'], n_tau=self.general_params['n_tau']) + + return triqs_solver + def _create_ftps_solver(self): r''' Initialize ftps solver instance @@ -1173,3 +1246,74 @@ def set_Gs_from_G_l(): self.perturbation_order = self.triqs_solver.histogram return + + def _w2dyn_cthyb_postprocessing(self): + r''' + Organize G_freq, G_time, Sigma_freq and G_l from w2dyn_cthyb solver + ''' + + def set_Gs_from_G_l(): + + # create new G_freq and G_time + for i, g in self.G_l: + g.enforce_discontinuity(np.identity(g.target_shape[0])) + # set G_freq from Legendre and Fouriertransform to get G_time + self.G_freq[i].set_from_legendre(g) + self.G_time[i].set_from_legendre(g) + + # Symmetrize + self.G_freq << make_hermitian(self.G_freq) + self.G_freq_unsym << self.G_freq + self.sum_k.symm_deg_gf(self.G_freq, ish=self.icrsh) + self.sum_k.symm_deg_gf(self.G_time, ish=self.icrsh) + # Dyson equation to get Sigma_freq + self.Sigma_freq << inverse(self.G0_freq) - inverse(self.G_freq) + + return + + # get Delta_time from solver + self.Delta_time << self.triqs_solver.Delta_tau + + # if measured in Legendre basis, get G_l from solver too + if self.solver_params['measure_G_l']: + # store original G_time into G_time_orig + self.G_time_orig << self.triqs_solver.G_tau + self.G_l << self.triqs_solver.G_l + # get G_time, G_freq, Sigma_freq from G_l + set_Gs_from_G_l() + + else: + self.G_freq << make_hermitian(self.triqs_solver.G_iw) + self.G_freq_unsym << self.G_freq + self.sum_k.symm_deg_gf(self.G_freq, ish=self.icrsh) + # set G_time + self.G_time << self.triqs_solver.G_tau + self.sum_k.symm_deg_gf(self.G_time, ish=self.icrsh) + + if self.general_params['legendre_fit']: + self.G_time_orig << self.triqs_solver.G_tau + # run the filter + self.G_l << legendre_filter.apply(self.G_time, self.general_params['n_l']) + # get G_time, G_freq, Sigma_freq from G_l + set_Gs_from_G_l() + elif self.solver_params['perform_tail_fit'] and not self.general_params['legendre_fit']: + # if tailfit has been used replace Sigma with the tail fitted Sigma from cthyb + self.Sigma_freq << self.triqs_solver.Sigma_iw + self.sum_k.symm_deg_gf(self.Sigma_freq, ish=self.icrsh) + else: + # obtain Sigma via dyson from symmetrized G_freq + self.Sigma_freq << inverse(self.G0_freq) - inverse(self.G_freq) + + # if density matrix is measured, get this too + if self.solver_params['measure_density_matrix']: + self.density_matrix = self.triqs_solver.density_matrix + self.h_loc_diagonalization = self.triqs_solver.h_loc_diagonalization + + if self.solver_params['measure_pert_order']: + self.perturbation_order = self.triqs_solver.perturbation_order + self.perturbation_order_total = self.triqs_solver.perturbation_order_total + + if self.general_params['measure_G2_iw_ph']: + pass # TODO + + return diff --git a/python/solid_dmft/read_config.py b/python/solid_dmft/read_config.py index af6f4304..b2833565 100755 --- a/python/solid_dmft/read_config.py +++ b/python/solid_dmft/read_config.py @@ -109,6 +109,7 @@ * 'hubbardI' * 'hartree' * 'ctseg' + * 'w2dyn_cthyb' n_iw : int, optional, default=1025 number of Matsubara frequencies @@ -345,6 +346,29 @@ Sigma_iw will automatically be calculated via http://dx.doi.org/10.1103/PhysRevB.85.205106 +w2dyn_cthyb parameters +====================== +measure_G2_iw_ph : bool, optional, default=False + Measure two-particle Green's function in particle-hole + frequency convention +measure_G2_n_bosonic : int, optional, default=30 + Number of bosonic Matsubara frequencies for two-particle + Green's function measurement +measure_G2_n_fermionic : int, optional, default=30 + Number of fermionic Matsubara frequencies for + two-particle Green's function measurement +worm_components : iterable of int, optional, default=None + Overrides worm components to measure with compound indices + (see ``w2dyn.auxiliaries.compound_index``) in iterable if not + None (Advanced setting: do not change if unsure) +move_global_prob : float, optional, default=0.005 + Overall probability of the global moves +flavourchange_moves : int, optional, default=0 + Use flavourchange-moves instead of 4-operator insertions; + they seemed to be inferior to 4 operator moves... +statesampling : int, optional, default=0 + Activate state-sampling algorithm + hartree parameters ================ with_fock : bool, optional, default=False @@ -499,7 +523,7 @@ 'n_l': {'converter': int, 'valid for': lambda x, _: x > 0, - 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'inchworm', 'hubbardI', 'ctseg'] + 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'inchworm', 'hubbardI', 'ctseg', 'w2dyn_cthyb'] and (params['solver']['measure_G_l'] or params['solver']['legendre_fit'])}, 'n_iw': {'converter': int, 'valid for': lambda x, _: x > 0, @@ -692,19 +716,19 @@ # cthyb parameters # 'length_cycle': {'converter': int, 'valid for': lambda x, _: x > 0, - 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg']}, + 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg', 'w2dyn_cthyb']}, 'n_warmup_cycles': {'converter': lambda s: int(float(s)), 'valid for': lambda x, _: x > 0, - 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg']}, + 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg', 'w2dyn_cthyb']}, 'n_cycles_tot': {'converter': lambda s: int(float(s)), 'valid for': lambda x, _: x >= 0, - 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg']}, + 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg', 'w2dyn_cthyb']}, 'max_time': {'converter': int, 'valid for': lambda x, _: x >= 0, 'default': None, - 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg']}, + 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg', 'w2dyn_cthyb']}, 'imag_threshold': {'converter': float, 'default': None, 'used': lambda params: params['general']['solver_type'] in ['cthyb']}, @@ -716,13 +740,13 @@ 'used': lambda params: params['general']['solver_type'] in ['cthyb']}, 'measure_G_tau': {'converter': BOOL_PARSER, 'default': True, - 'used': lambda params: params['general']['solver_type'] in ['hubbardI', 'ctseg']}, + 'used': lambda params: params['general']['solver_type'] in ['hubbardI', 'ctseg', 'w2dyn_cthyb']}, 'measure_G_iw': {'converter': BOOL_PARSER, 'default': False, 'used': lambda params: params['general']['solver_type'] in ['ctseg']}, 'measure_G_l': {'converter': BOOL_PARSER, 'default': False, - 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'hubbardI', 'ctseg']}, + 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'hubbardI', 'ctseg', 'w2dyn_cthyb']}, 'measure_density_matrix': {'converter': BOOL_PARSER, 'default': False, 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'hubbardI']}, @@ -731,13 +755,13 @@ 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint']}, 'measure_pert_order': {'converter': BOOL_PARSER, 'default': False, - 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg']}, + 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg', 'w2dyn_cthyb']}, 'move_shift': {'converter': BOOL_PARSER, 'default': False, 'used': lambda params: params['general']['solver_type'] in ['cthyb']}, 'random_seed': {'converter': str, 'default': None, - 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg']}, + 'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg', 'w2dyn_cthyb']}, 'perform_tail_fit': {'converter': BOOL_PARSER, 'used': lambda params: params['general']['solver_type'] in ['cthyb'] @@ -863,6 +887,24 @@ 'used': lambda params: params['general']['solver_type'] in ['ftps']}, 'dmrg_tw': {'converter': float, 'valid for': lambda x, _: x > 0, 'default': 1e-9, 'used': lambda params: params['general']['solver_type'] in ['ftps']}, + # + # w2dyn_cthyb parameters + # + 'measure_G2_iw_ph': {'converter': BOOL_PARSER, 'default': False, + 'used': lambda params: params['general']['solver_type'] in ['w2dyn_cthyb']}, + 'measure_G2_n_bosonic': {'converter': int, 'valid for': lambda x, _: x > 0, 'default': 30, + 'used': lambda params: params['general']['solver_type'] in ['w2dyn_cthyb']}, + 'measure_G2_n_fermionic': {'converter': int, 'valid for': lambda x, _: x > 0, 'default': 30, + 'used': lambda params: params['general']['solver_type'] in ['w2dyn_cthyb']}, + 'worm_components': {'converter': lambda s: list(map(int, s.split(','))), 'default': None, + 'used': lambda params: params['general']['solver_type'] in ['w2dyn_cthyb']}, + 'move_global_prob': {'converter': float, 'default': 0.005, + 'used': lambda params: params['general']['solver_type'] in ['w2dyn_cthyb']}, + 'flavourchange_moves': {'converter': int, 'default': 0, + 'used': lambda params: params['general']['solver_type'] in ['w2dyn_cthyb']}, + 'statesampling': {'converter': int, 'default': 0, + 'used': lambda params: params['general']['solver_type'] in ['w2dyn_cthyb']}, + }, 'advanced': {'dc_factor': {'converter': float, 'used': True, 'default': 'none'}, @@ -1235,7 +1277,7 @@ def read_config(config_file): if isinstance(parameters['advanced']['map_solver_struct'], dict): parameters['advanced']['map_solver_struct'] = [parameters['advanced']['map_solver_struct']] - if parameters['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg']: + if parameters['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg', 'w2dyn_cthyb']: parameters['solver']['n_cycles'] = parameters['solver']['n_cycles_tot'] // mpi.size del parameters['solver']['n_cycles_tot']