Add parameters to enable improved estimators

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,21 +366,29 @@ 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:
 
             if self.complex:
                 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

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",

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