diff --git a/python/w2dyn_cthyb/solver.py b/python/w2dyn_cthyb/solver.py index 7614d7f..a3e38cd 100644 --- a/python/w2dyn_cthyb/solver.py +++ b/python/w2dyn_cthyb/solver.py @@ -108,10 +108,12 @@ def solve(self, **params_kw): n_cycles = params_kw.pop("n_cycles") ### what does the True or False mean? n_warmup_cycles = params_kw.pop("n_warmup_cycles", 5000) ### default max_time = params_kw.pop("max_time", -1) + selfenergy = params_kw.pop("selfenergy", "dyson") worm = params_kw.pop("worm", False) percentageworminsert = params_kw.pop("PercentageWormInsert", 0.20) percentagewormreplace = params_kw.pop("PercentageWormReplace", 0.20) wormcomponents = params_kw.pop("worm_components", None) + wormsampling = worm or (selfenergy in ["improved_worm", "symmetric_improved_worm"]) length_cycle = params_kw.pop("length_cycle", 50) h_int = params_kw.pop("h_int") @@ -124,7 +126,7 @@ def solve(self, **params_kw): random_seed = params_kw.pop("random_seed", 1) move_double = params_kw.pop("move_double", True) - measure_G_l = params_kw.pop("measure_G_l", False) + measure_G_l = params_kw.pop("measure_G_l", True) measure_G_tau = params_kw.pop("measure_G_tau", True) measure_pert_order = params_kw.pop("measure_pert_order", False) statesampling = params_kw.pop("statesampling", False) @@ -215,7 +217,7 @@ def solve(self, **params_kw): cfg["QMC"]["offdiag"] = 1 ### complex worms are not yet existing - if self.complex and worm: + if self.complex and wormsampling: print('complex and worm together not yet implemented') exit() @@ -247,6 +249,8 @@ def solve(self, **params_kw): ### in a more sophisticated way cfg["General"]["beta"] = self.beta + cfg["General"]["siw_moments"] = "estimate" + cfg["General"]["SelfEnergy"] = selfenergy cfg["QMC"]["Niw"] = self.n_iw cfg["QMC"]["Ntau"] = self.n_tau * 2 # use double resolution bins & down sample to Triqs l8r @@ -267,14 +271,20 @@ def solve(self, **params_kw): else: cfg["QMC"]["statesampling"] = 0 - if worm: - + if selfenergy == "improved_worm": + cfg["QMC"]["WormMeasGSigmaiw"] = 1 + cfg["General"]["FTType"] = "none_worm" + elif selfenergy == "symmetric_improved_worm": + cfg["QMC"]["WormMeasQQ"] = 1 + cfg["General"]["FTType"] = "none_worm" + elif worm: # Do not enable measurements if cfg_qmc is supplied in the solve call if not 'cfg_qmc' in params_kw: cfg["QMC"]["WormMeasGiw"] = 1 cfg["QMC"]["WormMeasGtau"] = 1 cfg["QMC"]["WormSearchEta"] = 1 + if wormsampling: ### set worm parameters to some default values if not set by user if percentageworminsert != 0.0: cfg["QMC"]["PercentageWormInsert"] = percentageworminsert @@ -356,6 +366,8 @@ def solve(self, **params_kw): ### solve impurity problem mccfgcontainer = [] + siw_method = cfg["General"]["SelfEnergy"] + smom_method = cfg["General"]["siw_moments"] iter_no = 1 if measure_G_tau or measure_G_l or measure_pert_order: @@ -363,14 +375,20 @@ def solve(self, **params_kw): solver.set_problem(imp_problem) solver.umatrix = U_ijkl result = solver.solve(mccfgcontainer) + result.postprocessing(siw_method, smom_method) gtau = result.other["gtau-full"] + giw = result.giw + siw = result.siw - elif not worm: + elif not wormsampling: solver.set_problem(imp_problem) solver.umatrix = U_ijkl result = solver.solve(iter_no, mccfgcontainer) + result.postprocessing(siw_method, smom_method) gtau = result.other["gtau-full"] + giw = result.giw + siw = result.siw elif worm_get_sector_index(cfg['QMC']) == 2: @@ -413,6 +431,7 @@ def solve(self, **params_kw): solver.set_problem(imp_problem) solver.umatrix = U_ijkl result_aux, result = solver.solve_component(1, 2, comp_ind, mccfgcontainer) + result.postprocessing(siw_method, smom_method) for i in list(result.other.keys()): @@ -437,6 +456,16 @@ def solve(self, **params_kw): gtau[0, b1, s1, b2, s2, :] = result.other[gtau_name] gtau = stat.DistributedSample(gtau, mpi_comm, ntotal=mpi.size) + giw = result.giw + siw = result.siw + + elif worm_get_sector_index(cfg['QMC']) == 3: + + raise NotImplementedError("improved_worm") + + elif worm_get_sector_index(cfg['QMC']) == 10: + + raise NotImplementedError("symmetric_improved_worm") elif cfg["QMC"]["FourPnt"] == 8: # Know that: worm == True and worm_get_sector_index(cfg['QMC']) != 2 @@ -603,6 +632,8 @@ def solve(self, **params_kw): gf_err.data[:] = gf.stderr() self.GF_worm_components.append((component, gf_mean, gf_err)) + else: + raise NotImplementedError("The chosen combination of parameters is not supported") # TRIQS/cthyb like interface to sample all components of the two-particle Green's function @@ -687,16 +718,11 @@ def bs_diagflat(bs_array): self.G_tau, self.G_tau_error = w2dyn_ndarray_to_triqs_BlockGF_tau_beta_ntau( gtau, self.beta, self.gf_struct) - self.G_iw = BlockGf(mesh=self.iw_mesh, gf_struct=self.gf_struct) - - ### I will use the FFT from triqs here... - for name, g in self.G_tau: - bl_size = g.target_shape[0] - known_moments = np.zeros((4, bl_size, bl_size), dtype=complex) - for i in range(bl_size): - known_moments[1,i,i] = 1 + self.G_iw, self.G_iw_error = w2dyn_ndarray_to_triqs_BlockGF_iw_beta_niw( + giw, self.n_iw, self.beta, self.gf_struct) - self.G_iw[name].set_from_fourier(g, known_moments) + self.Sigma_iw, self.Sigma_iw_error = w2dyn_ndarray_to_triqs_BlockGF_iw_beta_niw( + siw, self.n_iw, self.beta, self.gf_struct) ### add perturbation order as observable #print 'measure_pert_order ', measure_pert_order diff --git a/test/python/2orb_Discrete_Bath.py b/test/python/2orb_Discrete_Bath.py index 7af9c4b..eb0e4c5 100755 --- a/test/python/2orb_Discrete_Bath.py +++ b/test/python/2orb_Discrete_Bath.py @@ -115,6 +115,7 @@ with HDFArchive("2orb_Discrete_Bath.out.h5",'w') as results: results["G_iw"] = S.G_iw results["G_tau"] = S.G_tau + results["Sigma_iw"] = S.Sigma_iw from triqs.utility.h5diff import h5diff if args.libcxx: @@ -151,6 +152,7 @@ with HDFArchive("2orb_Discrete_Bath.delta_interface.out.h5",'w') as results: results["G_iw"] = S.G_iw results["G_tau"] = S.G_tau + results["Sigma_iw"] = S.Sigma_iw if args.libcxx: h5diff("2orb_Discrete_Bath.libcxx.ref.h5", "2orb_Discrete_Bath.delta_interface.out.h5", diff --git a/test/python/2orb_Discrete_Bath.ref.h5 b/test/python/2orb_Discrete_Bath.ref.h5 index 1520b33..c471854 100644 Binary files a/test/python/2orb_Discrete_Bath.ref.h5 and b/test/python/2orb_Discrete_Bath.ref.h5 differ diff --git a/test/python/SIAM_Discrete_Bath.py b/test/python/SIAM_Discrete_Bath.py index 3e945fe..2a3bcef 100755 --- a/test/python/SIAM_Discrete_Bath.py +++ b/test/python/SIAM_Discrete_Bath.py @@ -82,6 +82,7 @@ with HDFArchive("SIAM_Discrete_Bath.out.h5",'w') as results: results["G_iw"] = S.G_iw results["G_tau"] = S.G_tau + results["Sigma_iw"] = S.Sigma_iw from triqs.utility.h5diff import h5diff if args.libcxx: @@ -117,6 +118,7 @@ with HDFArchive("SIAM_Discrete_Bath.delta_interface.out.h5",'w') as results: results["G_iw"] = S.G_iw results["G_tau"] = S.G_tau + results["Sigma_iw"] = S.Sigma_iw if args.libcxx: h5diff("SIAM_Discrete_Bath.libcxx.ref.h5","SIAM_Discrete_Bath.delta_interface.out.h5") diff --git a/test/python/SIAM_Discrete_Bath.ref.h5 b/test/python/SIAM_Discrete_Bath.ref.h5 index a3c3cf1..c0b36bf 100644 Binary files a/test/python/SIAM_Discrete_Bath.ref.h5 and b/test/python/SIAM_Discrete_Bath.ref.h5 differ