From ec8da69e347cdc3bf21935e50050031ab46726eb Mon Sep 17 00:00:00 2001
From: Jonathan Karp <jkarp314@gmail.com>
Date: Fri, 27 Aug 2021 17:31:07 -0400
Subject: [PATCH 01/11] edit SumKDFT class to take gf mesh at initialization
 and force Sigma to have that mesh

---
 python/triqs_dft_tools/sumk_dft.py       | 129 +++++++++++------------
 python/triqs_dft_tools/sumk_dft_tools.py |  96 +++++++----------
 2 files changed, 96 insertions(+), 129 deletions(-)

diff --git a/python/triqs_dft_tools/sumk_dft.py b/python/triqs_dft_tools/sumk_dft.py
index e1c474c0..e78d45ae 100644
--- a/python/triqs_dft_tools/sumk_dft.py
+++ b/python/triqs_dft_tools/sumk_dft.py
@@ -42,7 +42,7 @@ from scipy.optimize import minimize
 class SumkDFT(object):
     """This class provides a general SumK method for combining ab-initio code and triqs."""
 
-    def __init__(self, hdf_file, h_field=0.0, use_dft_blocks=False,
+    def __init__(self, hdf_file, h_field=0.0, mesh=None, beta=40, n_iw=1025, use_dft_blocks=False,
                  dft_data='dft_input', symmcorr_data='dft_symmcorr_input', parproj_data='dft_parproj_input',
                  symmpar_data='dft_symmpar_input', bands_data='dft_bands_input', transp_data='dft_transp_input',
                  misc_data='dft_misc_input',bc_data='dft_bandchar_input',fs_data='dft_fs_input'):
@@ -99,6 +99,14 @@ class SumkDFT(object):
             self.fs_data = fs_data
             self.h_field = h_field
 
+            if mesh is None:
+                self.mesh = MeshImFreq(beta=beta, S='Fermion',n_max=n_iw)
+            elif isinstance(mesh, MeshImFreq) or isinstance(mesh, MeshReFreq):
+                self.mesh = mesh
+            else:
+                raise ValueError('mesh must be a triqs mesh of type MeshImFreq or MeshReFreq')
+                
+
             self.block_structure = BlockStructure()
 
             # Read input from HDF:
@@ -467,7 +475,7 @@ class SumkDFT(object):
 
         return gf_rotated
 
-    def lattice_gf(self, ik, mu=None, iw_or_w="iw", beta=40, broadening=None, mesh=None, with_Sigma=True, with_dc=True):
+    def lattice_gf(self, ik, mu=None, broadening=None, mesh=None, with_Sigma=True, with_dc=True):
         r"""
         Calculates the lattice Green function for a given k-point from the DFT Hamiltonian and the self energy.
 
@@ -508,9 +516,7 @@ class SumkDFT(object):
             mu = self.chemical_potential
         ntoi = self.spin_names_to_ind[self.SO]
         spn = self.spin_block_names[self.SO]
-        if (iw_or_w != "iw") and (iw_or_w != "w"):
-            raise ValueError("lattice_gf: Implemented only for Re/Im frequency functions.")
-        if not hasattr(self, "Sigma_imp_" + iw_or_w):
+        if self.Sigma_inp is None:
             with_Sigma = False
         if broadening is None:
             if mesh is None:
@@ -518,45 +524,32 @@ class SumkDFT(object):
             else:  # broadening = 2 * \Delta omega, where \Delta omega is the spacing of omega points
                 broadening = 2.0 * ((mesh[1] - mesh[0]) / (mesh[2] - 1))
 
-        # Are we including Sigma?
-        if with_Sigma:
-            Sigma_imp = getattr(self, "Sigma_imp_" + iw_or_w)
-            sigma_minus_dc = [s.copy() for s in Sigma_imp]
-            if with_dc:
-                sigma_minus_dc = self.add_dc(iw_or_w)
-            if iw_or_w == "iw":
-                # override beta if Sigma_iw is present
-                beta = Sigma_imp[0].mesh.beta
-                mesh = Sigma_imp[0].mesh
-            elif iw_or_w == "w":
-                mesh = Sigma_imp[0].mesh
-                if broadening>0 and mpi.is_master_node():
-                    warn('lattice_gf called with Sigma and broadening > 0 (broadening = {}). You might want to explicitly set the broadening to 0.'.format(broadening))
-        else:
-            if iw_or_w == "iw":
-                if beta is None:
-                    raise ValueError("lattice_gf: Give the beta for the lattice GfReFreq.")
-                # Default number of Matsubara frequencies
-                mesh = MeshImFreq(beta=beta, S='Fermion', n_max=1025)
-            elif iw_or_w == "w":
-                if mesh is None:
-                    raise ValueError("lattice_gf: Give the mesh=(om_min,om_max,n_points) for the lattice GfReFreq.")
-                mesh = MeshReFreq(mesh[0], mesh[1], mesh[2])
-
         # Check if G_latt is present
         set_up_G_latt = False                       # Assume not
-        if not hasattr(self, "G_latt_" + iw_or_w):
+        if not hasattr(self, "G_latt" ):
             # Need to create G_latt_(i)w
             set_up_G_latt = True
         else:                                       # Check that existing GF is consistent
-            G_latt = getattr(self, "G_latt_" + iw_or_w)
+            G_latt = self.G_latt
             GFsize = [gf.target_shape[0] for bname, gf in G_latt]
             unchangedsize = all([self.n_orbitals[ik, ntoi[spn[isp]]] == GFsize[
                                 isp] for isp in range(self.n_spin_blocks[self.SO])])
-            if not unchangedsize:
+            if (not mesh is None) or (not unchangedsize):
                 set_up_G_latt = True
-            if (iw_or_w == "iw") and (self.G_latt_iw.mesh.beta != beta):
-                set_up_G_latt = True  # additional check for ImFreq
+
+        # Are we including Sigma?
+        if with_Sigma:
+            Sigma_imp = self.Sigma_imp
+            sigma_minus_dc = [s.copy() for s in Sigma_imp]
+            if with_dc:
+                sigma_minus_dc = self.add_dc()
+            if isinstance(self.mesh, MeshReFreq) and broadening > 0 and mpi.is_master_node():
+                warn('lattice_gf called with Sigma and broadening > 0 (broadening = {}). You might want to explicitly set the broadening to 0.'.format(broadening))
+        elif not mesh is None:
+            mesh = MeshReFreq(mesh[0], mesh[1], mesh[2])
+
+        if mesh is None:
+            mesh = self.mesh
 
         # Set up G_latt
         if set_up_G_latt:
@@ -565,19 +558,19 @@ class SumkDFT(object):
             gf_struct = [(spn[isp], block_structure[isp])
                          for isp in range(self.n_spin_blocks[self.SO])]
             block_ind_list = [block for block, inner in gf_struct]
-            if iw_or_w == "iw":
+            if isinstance(mesh, MeshImFreq):
                 glist = lambda: [GfImFreq(indices=inner, mesh=mesh)
                                  for block, inner in gf_struct]
-            elif iw_or_w == "w":
+            else:
                 glist = lambda: [GfReFreq(indices=inner, mesh=mesh)
                                  for block, inner in gf_struct]
             G_latt = BlockGf(name_list=block_ind_list,
                              block_list=glist(), make_copies=False)
             G_latt.zero()
 
-        if iw_or_w == "iw":
+        if isinstance(mesh, MeshImFreq):
             G_latt << iOmega_n
-        elif iw_or_w == "w":
+        else:
             G_latt << Omega + 1j * broadening
 
         idmat = [numpy.identity(
@@ -599,7 +592,7 @@ class SumkDFT(object):
                                       sigma_minus_dc[icrsh][bname], gf)
 
         G_latt.invert()
-        setattr(self, "G_latt_" + iw_or_w, G_latt)
+        self.G_latt = G_latt
 
         return G_latt
 
@@ -629,19 +622,19 @@ class SumkDFT(object):
         assert len(Sigma_imp) == self.n_corr_shells,\
             "put_Sigma: give exactly one Sigma for each corr. shell!"
 
-        if all((isinstance(gf, Gf) and isinstance(gf.mesh, MeshImFreq)) for bname, gf in Sigma_imp[0]):
+        if isinstance(self.mesh, MeshImFreq) and all(isinstance(gf, Gf) and gf.mesh == self.mesh for bname, gf in Sigma_imp[0]):
             # Imaginary frequency Sigma:
-            self.Sigma_imp_iw = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[icrsh].mesh, space='sumk')
+            self.Sigma_imp = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[icrsh].mesh, space='sumk')
                                  for icrsh in range(self.n_corr_shells)]
-            SK_Sigma_imp = self.Sigma_imp_iw
-        elif all(isinstance(gf, Gf) and isinstance(gf.mesh, MeshReFreq) for bname, gf in Sigma_imp[0]):
+            SK_Sigma_imp = self.Sigma_imp
+        elif isinstance(self.mesh, MeshReFreq) and all(isinstance(gf, Gf) and gf.mesh == self.mesh for bname, gf in Sigma_imp[0]):
             # Real frequency Sigma:
-            self.Sigma_imp_w = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[icrsh].mesh, gf_function=GfReFreq, space='sumk')
+            self.Sigma_imp = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[icrsh].mesh, gf_function=GfReFreq, space='sumk')
                                 for icrsh in range(self.n_corr_shells)]
-            SK_Sigma_imp = self.Sigma_imp_w
+            SK_Sigma_imp = self.Sigma_imp
 
         else:
-            raise ValueError("put_Sigma: This type of Sigma is not handled, give either BlockGf of GfReFreq or GfImFreq.")
+            raise ValueError("put_Sigma: Sigma_imp must have the same mesh as SumKDFT.mesh.")
 
         # rotation from local to global coordinate system:
         for icrsh in range(self.n_corr_shells):
@@ -654,13 +647,12 @@ class SumkDFT(object):
                     gf << Sigma_imp[icrsh][bname]
 
         #warning if real frequency self energy is within the bounds of the band energies
-        if isinstance(Sigma_imp[0].mesh, MeshReFreq):
+        if isinstance(self.mesh, MeshReFreq):
             if self.min_band_energy is None or self.max_band_energy is None:
                 self.calculate_min_max_band_energies()
-            for gf in Sigma_imp:
-                Sigma_mesh = numpy.array([i for i in gf.mesh.values()])
-                if Sigma_mesh[0] > (self.min_band_energy - self.chemical_potential) or Sigma_mesh[-1] < (self.max_band_energy - self.chemical_potential):
-                    warn('The given Sigma is on a mesh which does not cover the band energy range. The Sigma MeshReFreq runs from %f to %f, while the band energy (minus the chemical potential) runs from %f to %f'%(Sigma_mesh[0], Sigma_mesh[-1], self.min_band_energy, self.max_band_energy))
+            mesh = numpy.array([i for i in self.mesh.values()])
+            if mesh[0] > (self.min_band_energy - self.chemical_potential) or mesh[-1] < (self.max_band_energy - self.chemical_potential):
+                warn('The given Sigma is on a mesh which does not cover the band energy range. The Sigma MeshReFreq runs from %f to %f, while the band energy (minus the chemical potential) runs from %f to %f'%(Sigma_mesh[0], Sigma_mesh[-1], self.min_band_energy, self.max_band_energy))
 
     def transform_to_sumk_blocks(self, Sigma_imp, Sigma_out=None):
         r""" transform Sigma from solver to sumk space
@@ -703,7 +695,7 @@ class SumkDFT(object):
                 G_out=Sigma_out[icrsh])
         return Sigma_out
 
-    def extract_G_loc(self, mu=None, iw_or_w='iw', with_Sigma=True, with_dc=True, broadening=None,
+    def extract_G_loc(self, mu=None, with_Sigma=True, with_dc=True, broadening=None,
                       transform_to_solver_blocks=True, show_warnings=True):
         r"""
         Extracts the local downfolded Green function by the Brillouin-zone integration of the lattice Green's function.
@@ -739,14 +731,14 @@ class SumkDFT(object):
         if mu is None:
             mu = self.chemical_potential
 
-        if iw_or_w == "iw":
-            G_loc = [self.Sigma_imp_iw[icrsh].copy() for icrsh in range(
+        if isinstance(self.mesh, MeshImFreq):
+            G_loc = [self.Sigma_imp[icrsh].copy() for icrsh in range(
                 self.n_corr_shells)]   # this list will be returned
             beta = G_loc[0].mesh.beta
             G_loc_inequiv = [BlockGf(name_block_generator=[(block, GfImFreq(target_shape=(block_dim, block_dim), mesh=G_loc[0].mesh)) for block, block_dim in self.gf_struct_solver[ish].items()],
                                      make_copies=False) for ish in range(self.n_inequiv_shells)]
-        elif iw_or_w == "w":
-            G_loc = [self.Sigma_imp_w[icrsh].copy() for icrsh in range(
+        else:
+            G_loc = [self.Sigma_imp[icrsh].copy() for icrsh in range(
                 self.n_corr_shells)]   # this list will be returned
             mesh = G_loc[0].mesh
             G_loc_inequiv = [BlockGf(name_block_generator=[(block, GfReFreq(target_shape=(block_dim, block_dim), mesh=mesh)) for block, block_dim in self.gf_struct_solver[ish].items()],
@@ -757,13 +749,12 @@ class SumkDFT(object):
 
         ikarray = numpy.array(list(range(self.n_k)))
         for ik in mpi.slice_array(ikarray):
-            if iw_or_w == 'iw':
+            if isinstance(self.mesh, MeshImFreq):
                 G_latt = self.lattice_gf(
-                    ik=ik, mu=mu, iw_or_w=iw_or_w, with_Sigma=with_Sigma, with_dc=with_dc, beta=beta)
-            elif iw_or_w == 'w':
-                mesh_parameters = (G_loc[0].mesh.omega_min,G_loc[0].mesh.omega_max,len(G_loc[0].mesh))
+                    ik=ik, mu=mu, with_Sigma=with_Sigma, with_dc=with_dc)
+            else:
                 G_latt = self.lattice_gf(
-                    ik=ik, mu=mu, iw_or_w=iw_or_w, with_Sigma=with_Sigma, with_dc=with_dc, broadening=broadening, mesh=mesh_parameters)
+                    ik=ik, mu=mu, with_Sigma=with_Sigma, with_dc=with_dc, broadening=broadening)
             G_latt *= self.bz_weights[ik]
 
             for icrsh in range(self.n_corr_shells):
@@ -1446,7 +1437,7 @@ class SumkDFT(object):
         return trans
 
 
-    def density_matrix(self, method='using_gf', beta=40.0):
+    def density_matrix(self, method='using_gf'):
         """Calculate density matrices in one of two ways.
 
         Parameters
@@ -1477,8 +1468,7 @@ class SumkDFT(object):
 
             if method == "using_gf":
 
-                G_latt_iw = self.lattice_gf(
-                    ik=ik, mu=self.chemical_potential, iw_or_w="iw", beta=beta)
+                G_latt_iw = self.lattice_gf(ik=ik, mu=self.chemical_potential)
                 G_latt_iw *= self.bz_weights[ik]
                 dm = G_latt_iw.density()
                 MMat = [dm[sp] for sp in self.spin_block_names[self.SO]]
@@ -1777,7 +1767,7 @@ class SumkDFT(object):
                     self.dc_imp[icrsh][sp] = numpy.dot(T.conjugate().transpose(),
                             numpy.dot(self.dc_imp[icrsh][sp], T))
 
-    def add_dc(self, iw_or_w="iw"):
+    def add_dc(self):
         r"""
         Subtracts the double counting term from the impurity self energy.
 
@@ -1796,8 +1786,7 @@ class SumkDFT(object):
         """
 
         # Be careful: Sigma_imp is already in the global coordinate system!!
-        sigma_minus_dc = [s.copy()
-                          for s in getattr(self, "Sigma_imp_" + iw_or_w)]
+        sigma_minus_dc = [s.copy() for s in self.Sigma_imp]
         for icrsh in range(self.n_corr_shells):
             for bname, gf in sigma_minus_dc[icrsh]:
                 # Transform dc_imp to global coordinate system
@@ -1869,7 +1858,7 @@ class SumkDFT(object):
                 else:
                     gf_to_symm[key].from_L_G_R(v, ss, v.conjugate().transpose())
 
-    def total_density(self, mu=None, iw_or_w="iw", with_Sigma=True, with_dc=True, broadening=None):
+    def total_density(self, mu=None, with_Sigma=True, with_dc=True, broadening=None):
         r"""
         Calculates the total charge within the energy window for a given chemical potential.
         The chemical potential is either given by parameter `mu` or, if it is not specified,
@@ -1922,7 +1911,7 @@ class SumkDFT(object):
         ikarray = numpy.array(list(range(self.n_k)))
         for ik in mpi.slice_array(ikarray):
             G_latt = self.lattice_gf(
-                ik=ik, mu=mu, iw_or_w=iw_or_w, with_Sigma=with_Sigma, with_dc=with_dc, broadening=broadening)
+                ik=ik, mu=mu, with_Sigma=with_Sigma, with_dc=with_dc, broadening=broadening)
             dens += self.bz_weights[ik] * G_latt.total_density()
         # collect data from mpi:
         dens = mpi.all_reduce(mpi.world, dens, lambda x, y: x + y)
diff --git a/python/triqs_dft_tools/sumk_dft_tools.py b/python/triqs_dft_tools/sumk_dft_tools.py
index ff159a9d..21e6cab3 100644
--- a/python/triqs_dft_tools/sumk_dft_tools.py
+++ b/python/triqs_dft_tools/sumk_dft_tools.py
@@ -41,17 +41,17 @@ class SumkDFTTools(SumkDFT):
     Extends the SumkDFT class with some tools for analysing the data.
     """
 
-    def __init__(self, hdf_file, h_field=0.0, use_dft_blocks=False, dft_data='dft_input', symmcorr_data='dft_symmcorr_input',
+    def __init__(self, hdf_file, h_field=0.0, mesh=None, bets=40, n_iw=1025, use_dft_blocks=False, dft_data='dft_input', symmcorr_data='dft_symmcorr_input',
                  parproj_data='dft_parproj_input', symmpar_data='dft_symmpar_input', bands_data='dft_bands_input',
                  transp_data='dft_transp_input', misc_data='dft_misc_input'):
         """
         Initialisation of the class. Parameters are exactly as for SumKDFT.
         """
 
-        SumkDFT.__init__(self, hdf_file=hdf_file, h_field=h_field, use_dft_blocks=use_dft_blocks,
-                         dft_data=dft_data, symmcorr_data=symmcorr_data, parproj_data=parproj_data,
-                         symmpar_data=symmpar_data, bands_data=bands_data, transp_data=transp_data,
-                         misc_data=misc_data)
+        SumkDFT.__init__(self, hdf_file=hdf_file, h_field=h_field, mesh=mesh, beta=beta, n_iw=n_iw,
+                        use_dft_blocks=use_dft_blocks, dft_data=dft_data, symmcorr_data=symmcorr_data,
+                        parproj_data=parproj_data, symmpar_data=symmpar_data, bands_data=bands_data,
+                        transp_data=transp_data, misc_data=misc_data)
 
     # Uses .data of only GfReFreq objects.
     def dos_wannier_basis(self, mu=None, broadening=None, mesh=None, with_Sigma=True, with_dc=True, save_to_file=True):
@@ -82,10 +82,9 @@ class SumkDFTTools(SumkDFT):
         DOSproj_orb : Dict of numpy arrays
                       DOS projected to atoms and resolved into orbital contributions.
         """
-        if (mesh is None) and (not with_Sigma):
-            raise ValueError("lattice_gf: Give the mesh=(om_min,om_max,n_points) for the lattice GfReFreq.")
-        if mesh is None:
-            om_mesh = [x.real for x in self.Sigma_imp_w[0].mesh]
+        if mesh is None or with_Sigma:
+            assert isinstance(self.mesh, MeshReFreq), "mesh must be given if self.mesh is a MeshImFreq"
+            om_mesh = [x.real for x in self.mesh]
             om_min = om_mesh[0]
             om_max = om_mesh[-1]
             n_om = len(om_mesh)
@@ -119,7 +118,7 @@ class SumkDFTTools(SumkDFT):
         for ik in mpi.slice_array(ikarray):
 
             G_latt_w = self.lattice_gf(
-                ik=ik, mu=mu, iw_or_w="w", broadening=broadening, mesh=mesh, with_Sigma=with_Sigma, with_dc=with_dc)
+                ik=ik, mu=mu, broadening=broadening, mesh=mesh, with_Sigma=with_Sigma, with_dc=with_dc)
             G_latt_w *= self.bz_weights[ik]
 
             # Non-projected DOS
@@ -218,10 +217,9 @@ class SumkDFTTools(SumkDFT):
         DOSproj_orb : Dict of numpy arrays
                       DOS projected to atoms and resolved into orbital contributions.
         """
-        if (mesh is None) and (not with_Sigma):
-            raise ValueError("lattice_gf: Give the mesh=(om_min,om_max,n_points) for the lattice GfReFreq.")
-        if mesh is None:
-            om_mesh = [x.real for x in self.Sigma_imp_w[0].mesh]
+        if mesh is None or with_Sigma:
+            assert isinstance(self.mesh, MeshReFreq), "mesh must be given if self.mesh is a MeshImFreq"
+            om_mesh = [x.real for x in self.mesh]
             om_min = om_mesh[0]
             om_max = om_mesh[-1]
             n_om = len(om_mesh)
@@ -255,7 +253,7 @@ class SumkDFTTools(SumkDFT):
         for ik in mpi.slice_array(ikarray):
 
             G_latt_w = self.lattice_gf(
-                ik=ik, mu=mu, iw_or_w="w", broadening=broadening, mesh=mesh, with_Sigma=with_Sigma, with_dc=with_dc)
+                ik=ik, mu=mu, broadening=broadening, mesh=mesh, with_Sigma=with_Sigma, with_dc=with_dc)
             G_latt_w *= self.bz_weights[ik]
 
             # Non-projected DOS
@@ -350,10 +348,9 @@ class SumkDFTTools(SumkDFT):
         if self.symm_op:
             self.symmpar = Symmetry(self.hdf_file, subgroup=self.symmpar_data)
 
-        if (mesh is None) and (not with_Sigma):
-            raise ValueError("lattice_gf: Give the mesh=(om_min,om_max,n_points) for the lattice GfReFreq.")
-        if mesh is None:
-            om_mesh = [x.real for x in self.Sigma_imp_w[0].mesh]
+        if mesh is None or with_Sigma:
+            assert isinstance(self.mesh, MeshReFreq), "mesh must be given if self.mesh is a MeshImFreq"
+            om_mesh = [x.real for x in self.mesh]
             om_min = om_mesh[0]
             om_max = om_mesh[-1]
             n_om = len(om_mesh)
@@ -389,7 +386,7 @@ class SumkDFTTools(SumkDFT):
         for ik in mpi.slice_array(ikarray):
 
             G_latt_w = self.lattice_gf(
-                ik=ik, mu=mu, iw_or_w="w", broadening=broadening, mesh=mesh, with_Sigma=with_Sigma, with_dc=with_dc)
+                ik=ik, mu=mu, broadening=broadening, mesh=mesh, with_Sigma=with_Sigma, with_dc=with_dc)
             G_latt_w *= self.bz_weights[ik]
 
             # Non-projected DOS
@@ -502,10 +499,9 @@ class SumkDFTTools(SumkDFT):
           if not value_read:
             return value_read
 
-        if (mesh is None) and (not with_Sigma):
-            raise ValueError("lattice_gf: Give the mesh=(om_min,om_max,n_points) for the lattice GfReFreq.")
-        if mesh is None:
-            om_mesh = [x.real for x in self.Sigma_imp_w[0].mesh]
+        if mesh is None or with_Sigma:
+            assert isinstance(self.mesh, MeshReFreq), "mesh must be given if self.mesh is a MeshImFreq"
+            om_mesh = [x.real for x in self.mesh]
             om_min = om_mesh[0]
             om_max = om_mesh[-1]
             n_om = len(om_mesh)
@@ -532,7 +528,7 @@ class SumkDFTTools(SumkDFT):
         for ik in mpi.slice_array(ikarray):
 
             G_latt_w = self.lattice_gf(
-                ik=ik, mu=mu, iw_or_w="w", broadening=broadening, mesh=mesh, with_Sigma=with_Sigma, with_dc=with_dc)
+                ik=ik, mu=mu, broadening=broadening, mesh=mesh, with_Sigma=with_Sigma, with_dc=with_dc)
             G_latt_w *= self.bz_weights[ik]
             if(nk!=None):
               for iom in range(n_om):
@@ -667,8 +663,9 @@ class SumkDFTTools(SumkDFT):
         if not value_read:
           return value_read
 
-        if with_Sigma is True:
-            om_mesh = [x.real for x in self.Sigma_imp_w[0].mesh]
+        if with_Sigma is True or mesh is None:
+            assert isinstance(self.mesh, MeshReFreq), "SumkDFT.mesh must be real if with_Sigma is True or mesh is not given"
+            om_mesh = [x.real for x in self.mesh]
             #for Fermi Surface calculations
             if FS:
               jw=[i for i in range(len(om_mesh)) if om_mesh[i] == 0.0]
@@ -690,17 +687,6 @@ class SumkDFTTools(SumkDFT):
             mesh = (om_min, om_max, n_om)
             if broadening is None:
                broadening=0.0
-        elif mesh is None:
-        #default is to set "mesh" to be just for the Fermi surface - omega=0.0
-            om_min = 0.000
-            om_max = 0.001
-            n_om = 3
-            mesh = (om_min, om_max, n_om)
-            om_mesh = numpy.linspace(om_min, om_max, n_om)
-            if broadening is None:
-               broadening=0.01
-            FS=True
-            jw=[i for i in range(len(om_mesh)) if((om_mesh[i]<=om_max)and(om_mesh[i]>=om_min))]
         else:
         #a range of frequencies can be used if desired
             om_min, om_max, n_om = mesh
@@ -732,7 +718,7 @@ class SumkDFTTools(SumkDFT):
             vkc[ik,:] = numpy.matmul(self.bmat,self.vkl[ik,:])
 
             G_latt_w = self.lattice_gf(
-                ik=ik, mu=mu, iw_or_w="w", broadening=broadening, mesh=mesh, with_Sigma=with_Sigma, with_dc=with_dc)
+                ik=ik, mu=mu, broadening=broadening, mesh=mesh, with_Sigma=with_Sigma, with_dc=with_dc)
 
             for iom in range(n_om):
               for bname, gf in G_latt_w:
@@ -854,7 +840,8 @@ class SumkDFTTools(SumkDFT):
         if mu is None:
             mu = self.chemical_potential
         spn = self.spin_block_names[self.SO]
-        mesh = numpy.array([x.real for x in self.Sigma_imp_w[0].mesh])
+        mesh = [x.real for x in self.mesh]
+        n_om = len(mesh)
 
         if plot_range is None:
             om_minplot = mesh[0] - 0.001
@@ -881,8 +868,7 @@ class SumkDFTTools(SumkDFT):
         ikarray = numpy.array(list(range(self.n_k)))
         for ik in mpi.slice_array(ikarray):
 
-            G_latt_w = self.lattice_gf(
-                ik=ik, mu=mu, iw_or_w="w", broadening=broadening)
+            G_latt_w = self.lattice_gf(ik=ik, mu=mu, broadening=broadening)
 
             if ishell is None:
                 # Non-projected A(k,w)
@@ -953,7 +939,7 @@ class SumkDFTTools(SumkDFT):
 
         return Akw
 
-    def partial_charges(self, beta=40, mu=None, with_Sigma=True, with_dc=True):
+    def partial_charges(self, mu=None, with_Sigma=True, with_dc=True):
         """
         Calculates the orbitally-resolved density matrix for all the orbitals considered in the input, consistent with
         the definition of Wien2k. Hence, (possibly non-orthonormal) projectors have to be provided in the partial projectors subgroup of
@@ -995,23 +981,16 @@ class SumkDFTTools(SumkDFT):
         # Set up G_loc
         gf_struct_parproj = [[(sp, self.shells[ish]['dim']) for sp in spn]
                              for ish in range(self.n_shells)]
-        if with_Sigma:
-            G_loc = [BlockGf(name_block_generator=[(block, GfImFreq(target_shape=(block_dim, block_dim), mesh=self.Sigma_imp_iw[0].mesh))
-                                                   for block, block_dim in gf_struct_parproj[ish]], make_copies=False)
-                     for ish in range(self.n_shells)]
-            beta = self.Sigma_imp_iw[0].mesh.beta
-        else:
-            G_loc = [BlockGf(name_block_generator=[(block, GfImFreq(target_shape=(block_dim, block_dim), beta=beta))
-                                                   for block, block_dim in gf_struct_parproj[ish]], make_copies=False)
-                     for ish in range(self.n_shells)]
+        G_loc = [BlockGf(name_block_generator=[(block, GfImFreq(target_shape=(block_dim, block_dim), mesh=self.mesh))
+                                                for block, block_dim in gf_struct_parproj[ish]], make_copies=False)
+                    for ish in range(self.n_shells)]
         for ish in range(self.n_shells):
             G_loc[ish].zero()
 
         ikarray = numpy.array(list(range(self.n_k)))
         for ik in mpi.slice_array(ikarray):
 
-            G_latt_iw = self.lattice_gf(
-                ik=ik, mu=mu, iw_or_w="iw", beta=beta, with_Sigma=with_Sigma, with_dc=with_dc)
+            G_latt_iw = self.lattice_gf(ik=ik, mu=mu, with_Sigma=with_Sigma, with_dc=with_dc)
             G_latt_iw *= self.bz_weights[ik]
             for ish in range(self.n_shells):
                 tmp = G_loc[ish].copy()
@@ -1207,7 +1186,7 @@ class SumkDFTTools(SumkDFT):
         # Define mesh for Green's function and in the specified energy window
         if (with_Sigma == True):
             self.omega = numpy.array([round(x.real, 12)
-                                      for x in self.Sigma_imp_w[0].mesh])
+                                      for x in self.mesh])
             mesh = None
             mu = self.chemical_potential
             n_om = len(self.omega)
@@ -1225,13 +1204,13 @@ class SumkDFTTools(SumkDFT):
                 # In the future there should be an option in gf to manipulate the mesh (e.g. truncate) directly.
                 # For now we stick with this:
                 for icrsh in range(self.n_corr_shells):
-                    Sigma_save = self.Sigma_imp_w[icrsh].copy()
+                    Sigma_save = self.Sigma_imp[icrsh].copy()
                     spn = self.spin_block_names[self.corr_shells[icrsh]['SO']]
                     glist = lambda: [GfReFreq(target_shape=(block_dim, block_dim), window=(self.omega[
                                               0], self.omega[-1]), n_points=n_om) for block, block_dim in self.gf_struct_sumk[icrsh]]
                     self.Sigma_imp_w[icrsh] = BlockGf(
                         name_list=spn, block_list=glist(), make_copies=False)
-                    for i, g in self.Sigma_imp_w[icrsh]:
+                    for i, g in self.Sigma_imp[icrsh]:
                         for iL in g.indices[0]:
                             for iR in g.indices[0]:
                                 for iom in range(n_om):
@@ -1267,8 +1246,7 @@ class SumkDFTTools(SumkDFT):
         ikarray = numpy.array(list(range(self.n_k)))
         for ik in mpi.slice_array(ikarray):
             # Calculate G_w  for ik and initialize A_kw
-            G_w = self.lattice_gf(ik, mu, iw_or_w="w", beta=beta,
-                                  broadening=broadening, mesh=mesh, with_Sigma=with_Sigma)
+            G_w = self.lattice_gf(ik, mu, broadening=broadening, mesh=mesh, with_Sigma=with_Sigma)
             A_kw = [numpy.zeros((self.n_orbitals[ik][isp], self.n_orbitals[ik][isp], n_om), dtype=numpy.complex_)
                     for isp in range(n_inequiv_spin_blocks)]
 

From 6df75d8c3977ee06cfcee7c2422a1f6234c35356 Mon Sep 17 00:00:00 2001
From: Jonathan Karp <jkarp314@gmail.com>
Date: Fri, 27 Aug 2021 17:46:17 -0400
Subject: [PATCH 02/11] bug fixes

---
 python/triqs_dft_tools/sumk_dft.py       | 2 +-
 python/triqs_dft_tools/sumk_dft_tools.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/python/triqs_dft_tools/sumk_dft.py b/python/triqs_dft_tools/sumk_dft.py
index e78d45ae..3e6a9309 100644
--- a/python/triqs_dft_tools/sumk_dft.py
+++ b/python/triqs_dft_tools/sumk_dft.py
@@ -516,7 +516,7 @@ class SumkDFT(object):
             mu = self.chemical_potential
         ntoi = self.spin_names_to_ind[self.SO]
         spn = self.spin_block_names[self.SO]
-        if self.Sigma_inp is None:
+        if self.Sigma_imp is None:
             with_Sigma = False
         if broadening is None:
             if mesh is None:
diff --git a/python/triqs_dft_tools/sumk_dft_tools.py b/python/triqs_dft_tools/sumk_dft_tools.py
index 21e6cab3..ef020eec 100644
--- a/python/triqs_dft_tools/sumk_dft_tools.py
+++ b/python/triqs_dft_tools/sumk_dft_tools.py
@@ -41,7 +41,7 @@ class SumkDFTTools(SumkDFT):
     Extends the SumkDFT class with some tools for analysing the data.
     """
 
-    def __init__(self, hdf_file, h_field=0.0, mesh=None, bets=40, n_iw=1025, use_dft_blocks=False, dft_data='dft_input', symmcorr_data='dft_symmcorr_input',
+    def __init__(self, hdf_file, h_field=0.0, mesh=None, beta=40, n_iw=1025, use_dft_blocks=False, dft_data='dft_input', symmcorr_data='dft_symmcorr_input',
                  parproj_data='dft_parproj_input', symmpar_data='dft_symmpar_input', bands_data='dft_bands_input',
                  transp_data='dft_transp_input', misc_data='dft_misc_input'):
         """

From 1e7c66805f9f0af409520ae9d562ea80ac422b84 Mon Sep 17 00:00:00 2001
From: Jonathan Karp <jkarp314@gmail.com>
Date: Mon, 30 Aug 2021 15:58:30 -0400
Subject: [PATCH 03/11] fix bugs and tests

---
 python/triqs_dft_tools/sumk_dft.py | 4 ++--
 test/python/sigma_from_file.py     | 2 +-
 test/python/srvo3_transp.py        | 3 +--
 test/python/sumkdft_basic.py       | 4 ++--
 4 files changed, 6 insertions(+), 7 deletions(-)

diff --git a/python/triqs_dft_tools/sumk_dft.py b/python/triqs_dft_tools/sumk_dft.py
index 3e6a9309..595c339a 100644
--- a/python/triqs_dft_tools/sumk_dft.py
+++ b/python/triqs_dft_tools/sumk_dft.py
@@ -516,7 +516,7 @@ class SumkDFT(object):
             mu = self.chemical_potential
         ntoi = self.spin_names_to_ind[self.SO]
         spn = self.spin_block_names[self.SO]
-        if self.Sigma_imp is None:
+        if not hasattr(self, "Sigma_imp"):
             with_Sigma = False
         if broadening is None:
             if mesh is None:
@@ -652,7 +652,7 @@ class SumkDFT(object):
                 self.calculate_min_max_band_energies()
             mesh = numpy.array([i for i in self.mesh.values()])
             if mesh[0] > (self.min_band_energy - self.chemical_potential) or mesh[-1] < (self.max_band_energy - self.chemical_potential):
-                warn('The given Sigma is on a mesh which does not cover the band energy range. The Sigma MeshReFreq runs from %f to %f, while the band energy (minus the chemical potential) runs from %f to %f'%(Sigma_mesh[0], Sigma_mesh[-1], self.min_band_energy, self.max_band_energy))
+                warn('The given Sigma is on a mesh which does not cover the band energy range. The Sigma MeshReFreq runs from %f to %f, while the band energy (minus the chemical potential) runs from %f to %f'%(mesh[0], mesh[-1], self.min_band_energy, self.max_band_energy))
 
     def transform_to_sumk_blocks(self, Sigma_imp, Sigma_out=None):
         r""" transform Sigma from solver to sumk space
diff --git a/test/python/sigma_from_file.py b/test/python/sigma_from_file.py
index 7096e14d..30471924 100644
--- a/test/python/sigma_from_file.py
+++ b/test/python/sigma_from_file.py
@@ -40,7 +40,7 @@ for name, s in Sigma_hdf:
     np.savetxt('Sigma_' + name + '.dat', mesh_a_data)
 
 # Read self energy from txt files
-SK = SumkDFTTools(hdf_file =  'SrVO3.ref.h5', use_dft_blocks = True)
+SK = SumkDFTTools(hdf_file='SrVO3.ref.h5', mesh=Sigma_hdf.mesh, use_dft_blocks=True)
 
 # the order in the orig SrVO3 file is not assured, hence order it here
 a_list = sorted([a for a,al in SK.gf_struct_solver[0].items()])
diff --git a/test/python/srvo3_transp.py b/test/python/srvo3_transp.py
index 278457e3..dcbf613b 100644
--- a/test/python/srvo3_transp.py
+++ b/test/python/srvo3_transp.py
@@ -32,10 +32,9 @@ Converter = Wien2kConverter(filename='SrVO3', repacking=True)
 Converter.convert_dft_input()
 Converter.convert_transport_input()
 
-SK = SumkDFTTools(hdf_file='SrVO3.ref.h5', use_dft_blocks=True)
-
 with HDFArchive('SrVO3_Sigma.h5', 'a') as ar:
     Sigma = ar['dmft_transp_input']['Sigma_w']
+    SK = SumkDFTTools(hdf_file='SrVO3.ref.h5', mesh=Sigma.mesh, use_dft_blocks=True)
     SK.set_Sigma([Sigma])
     SK.chemical_potential = ar['dmft_transp_input']['chemical_potential']
     SK.dc_imp = ar['dmft_transp_input']['dc_imp']
diff --git a/test/python/sumkdft_basic.py b/test/python/sumkdft_basic.py
index ebead200..5f6d8f39 100644
--- a/test/python/sumkdft_basic.py
+++ b/test/python/sumkdft_basic.py
@@ -28,8 +28,8 @@ from triqs.utility.h5diff import h5diff
 
 SK = SumkDFTTools(hdf_file = 'SrVO3.ref.h5')
 
-dm = SK.density_matrix(method = 'using_gf', beta = 40)
-dm_pc = SK.partial_charges(beta=40,with_Sigma=False,with_dc=False)
+dm = SK.density_matrix(method = 'using_gf')
+dm_pc = SK.partial_charges(with_Sigma=False, with_dc=False)
 
 with HDFArchive('sumkdft_basic.out.h5','w') as ar:
     ar['dm'] = dm

From fc893a0f489592edd7f7b8f1c30bc9b9bab8d935 Mon Sep 17 00:00:00 2001
From: Jonathan Karp <jkarp314@gmail.com>
Date: Thu, 23 Sep 2021 13:44:50 -0400
Subject: [PATCH 04/11] use mesh from Sigma instead of Sumk in lattice_gf

---
 python/triqs_dft_tools/sumk_dft.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/python/triqs_dft_tools/sumk_dft.py b/python/triqs_dft_tools/sumk_dft.py
index 595c339a..bf5127d4 100644
--- a/python/triqs_dft_tools/sumk_dft.py
+++ b/python/triqs_dft_tools/sumk_dft.py
@@ -543,12 +543,12 @@ class SumkDFT(object):
             sigma_minus_dc = [s.copy() for s in Sigma_imp]
             if with_dc:
                 sigma_minus_dc = self.add_dc()
-            if isinstance(self.mesh, MeshReFreq) and broadening > 0 and mpi.is_master_node():
+            mesh = Sigma_imp[0].mesh
+            if isinstance(mesh, MeshReFreq) and broadening > 0 and mpi.is_master_node():
                 warn('lattice_gf called with Sigma and broadening > 0 (broadening = {}). You might want to explicitly set the broadening to 0.'.format(broadening))
         elif not mesh is None:
             mesh = MeshReFreq(mesh[0], mesh[1], mesh[2])
-
-        if mesh is None:
+        else:
             mesh = self.mesh
 
         # Set up G_latt

From 5c8071a6acc995c96abe49d6cc9811f6c05a73aa Mon Sep 17 00:00:00 2001
From: Jonathan Karp <jkarp314@gmail.com>
Date: Mon, 4 Oct 2021 12:58:44 -0400
Subject: [PATCH 05/11] update docstrings and give warning in lattice_gf if
 with Sigma and mesh is given

---
 python/triqs_dft_tools/sumk_dft.py       | 14 +++++++-------
 python/triqs_dft_tools/sumk_dft_tools.py |  2 --
 2 files changed, 7 insertions(+), 9 deletions(-)

diff --git a/python/triqs_dft_tools/sumk_dft.py b/python/triqs_dft_tools/sumk_dft.py
index bf5127d4..ce1d3f0c 100644
--- a/python/triqs_dft_tools/sumk_dft.py
+++ b/python/triqs_dft_tools/sumk_dft.py
@@ -57,6 +57,11 @@ class SumkDFT(object):
                   The value of magnetic field to add to the DFT Hamiltonian.
                   The contribution -h_field*sigma is added to diagonal elements of the Hamiltonian.
                   It cannot be used with the spin-orbit coupling on; namely h_field is set to 0 if self.SO=True.
+        mesh: MeshImFreq or MeshImFreq, optional. Frequency mesh of Sigma.
+        beta : real, optional
+               Inverse temperature. Used to construct imaginary frequency if mesh is not given. 
+        n_iw : integer, optional
+               Number of Matsubara frequencies. Used to construct imaginary frequency if mesh is not given. 
         use_dft_blocks : boolean, optional
                          If True, the local Green's function matrix for each spin is divided into smaller blocks
                           with the block structure determined from the DFT density matrix of the corresponding correlated shell.
@@ -486,13 +491,6 @@ class SumkDFT(object):
         mu : real, optional
              Chemical potential for which the Green's function is to be calculated.
              If not provided, self.chemical_potential is used for mu.
-        iw_or_w : string, optional
-
-                  - `iw_or_w` = 'iw' for a imaginary-frequency self-energy
-                  - `iw_or_w` = 'w' for a real-frequency self-energy
-
-        beta : real, optional
-               Inverse temperature.
         broadening : real, optional
                      Imaginary shift for the axis along which the real-axis GF is calculated.
                      If not provided, broadening will be set to double of the distance between mesh points in 'mesh'.
@@ -543,6 +541,8 @@ class SumkDFT(object):
             sigma_minus_dc = [s.copy() for s in Sigma_imp]
             if with_dc:
                 sigma_minus_dc = self.add_dc()
+            if not mesh is None:
+                warn('lattice_gf called with Sigma and given mesh. Mesh will be taken from Sigma.')
             mesh = Sigma_imp[0].mesh
             if isinstance(mesh, MeshReFreq) and broadening > 0 and mpi.is_master_node():
                 warn('lattice_gf called with Sigma and broadening > 0 (broadening = {}). You might want to explicitly set the broadening to 0.'.format(broadening))
diff --git a/python/triqs_dft_tools/sumk_dft_tools.py b/python/triqs_dft_tools/sumk_dft_tools.py
index ef020eec..5f47c03a 100644
--- a/python/triqs_dft_tools/sumk_dft_tools.py
+++ b/python/triqs_dft_tools/sumk_dft_tools.py
@@ -950,8 +950,6 @@ class SumkDFTTools(SumkDFT):
 
         with_Sigma : boolean, optional
                      If True, the self energy is used for the calculation. If false, partial charges are calculated without self-energy correction.
-        beta : double, optional
-               In case the self-energy correction is not used, the inverse temperature where the calculation should be done has to be given here.
         mu : double, optional
              Chemical potential, overrides the one stored in the hdf5 archive.
         with_dc : boolean, optional

From 5f2a0c763e6050dd1daa0ede71f75cec5543e62a Mon Sep 17 00:00:00 2001
From: Jonathan Karp <jkarp314@gmail.com>
Date: Tue, 5 Oct 2021 14:13:17 -0400
Subject: [PATCH 06/11] make mesh parameter in lattice_gf a triqs mesh object

---
 python/triqs_dft_tools/sumk_dft.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/python/triqs_dft_tools/sumk_dft.py b/python/triqs_dft_tools/sumk_dft.py
index ce1d3f0c..91b4e2fd 100644
--- a/python/triqs_dft_tools/sumk_dft.py
+++ b/python/triqs_dft_tools/sumk_dft.py
@@ -494,9 +494,8 @@ class SumkDFT(object):
         broadening : real, optional
                      Imaginary shift for the axis along which the real-axis GF is calculated.
                      If not provided, broadening will be set to double of the distance between mesh points in 'mesh'.
-        mesh : list, optional
-               Data defining mesh on which the real-axis GF will be calculated, given in the form
-               (om_min,om_max,n_points), where om_min is the minimum omega, om_max is the maximum omega and n_points is the number of points.
+        mesh : MeshReFreq or MeshImFreq, optional
+                    Mesh to be used if with_Sigma=False. If with Sigma=False and mesh is none then self.mesh is used.
         with_Sigma : boolean, optional
                      If True the GF will be calculated with the self-energy stored in self.Sigmaimp_(w/iw), for real/Matsubara GF, respectively.
                      In this case the mesh is taken from the self.Sigma_imp object.
@@ -547,7 +546,7 @@ class SumkDFT(object):
             if isinstance(mesh, MeshReFreq) and broadening > 0 and mpi.is_master_node():
                 warn('lattice_gf called with Sigma and broadening > 0 (broadening = {}). You might want to explicitly set the broadening to 0.'.format(broadening))
         elif not mesh is None:
-            mesh = MeshReFreq(mesh[0], mesh[1], mesh[2])
+            assert isinstance(mesh, MreshReFreq) or isinstance(mesh, MeshImFreq),  "mesh must be a triqs MeshReFreq or MeshImFreq"
         else:
             mesh = self.mesh
 

From bb4682fafad8e7a6363ff77ed77877e3f34ccb62 Mon Sep 17 00:00:00 2001
From: Jonathan Karp <jkarp314@gmail.com>
Date: Tue, 5 Oct 2021 15:09:52 -0400
Subject: [PATCH 07/11] warning in extract_G_loc if SumkDFT and Sigma have
 different mesh

---
 python/triqs_dft_tools/sumk_dft.py | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/python/triqs_dft_tools/sumk_dft.py b/python/triqs_dft_tools/sumk_dft.py
index 91b4e2fd..2b2890f1 100644
--- a/python/triqs_dft_tools/sumk_dft.py
+++ b/python/triqs_dft_tools/sumk_dft.py
@@ -730,7 +730,15 @@ class SumkDFT(object):
         if mu is None:
             mu = self.chemical_potential
 
-        if isinstance(self.mesh, MeshImFreq):
+        if with_Sigma:
+            mesh = self.Sigma_imp[0].mesh
+            if mesh != self.mesh:
+                warn('self.mesh and self.Sigma_imp[0].mesh are differen! Using mesh from Sigma')
+
+        else:
+            mesh = self.mesh
+
+        if isinstance(mesh, MeshImFreq):
             G_loc = [self.Sigma_imp[icrsh].copy() for icrsh in range(
                 self.n_corr_shells)]   # this list will be returned
             beta = G_loc[0].mesh.beta

From a36f2fdc61d7a77b1264a0482bea303292a57b9d Mon Sep 17 00:00:00 2001
From: Alexander Hampel <ahampel@flatironinstitute.org>
Date: Wed, 15 Dec 2021 09:01:10 -0500
Subject: [PATCH 08/11] fix spaghettis function and update test

---
 python/triqs_dft_tools/sumk_dft.py       | 10 +++++-----
 python/triqs_dft_tools/sumk_dft_tools.py |  9 +++++----
 test/python/srvo3_spectral.py            | 24 ++++++++++--------------
 3 files changed, 20 insertions(+), 23 deletions(-)

diff --git a/python/triqs_dft_tools/sumk_dft.py b/python/triqs_dft_tools/sumk_dft.py
index 2b2890f1..72d8cd17 100644
--- a/python/triqs_dft_tools/sumk_dft.py
+++ b/python/triqs_dft_tools/sumk_dft.py
@@ -59,9 +59,9 @@ class SumkDFT(object):
                   It cannot be used with the spin-orbit coupling on; namely h_field is set to 0 if self.SO=True.
         mesh: MeshImFreq or MeshImFreq, optional. Frequency mesh of Sigma.
         beta : real, optional
-               Inverse temperature. Used to construct imaginary frequency if mesh is not given. 
+               Inverse temperature. Used to construct imaginary frequency if mesh is not given.
         n_iw : integer, optional
-               Number of Matsubara frequencies. Used to construct imaginary frequency if mesh is not given. 
+               Number of Matsubara frequencies. Used to construct imaginary frequency if mesh is not given.
         use_dft_blocks : boolean, optional
                          If True, the local Green's function matrix for each spin is divided into smaller blocks
                           with the block structure determined from the DFT density matrix of the corresponding correlated shell.
@@ -110,7 +110,7 @@ class SumkDFT(object):
                 self.mesh = mesh
             else:
                 raise ValueError('mesh must be a triqs mesh of type MeshImFreq or MeshReFreq')
-                
+
 
             self.block_structure = BlockStructure()
 
@@ -621,12 +621,12 @@ class SumkDFT(object):
         assert len(Sigma_imp) == self.n_corr_shells,\
             "put_Sigma: give exactly one Sigma for each corr. shell!"
 
-        if isinstance(self.mesh, MeshImFreq) and all(isinstance(gf, Gf) and gf.mesh == self.mesh for bname, gf in Sigma_imp[0]):
+        if isinstance(self.mesh, MeshImFreq) and all(isinstance(gf.mesh, MeshImFreq) and isinstance(gf, Gf) and gf.mesh == self.mesh for bname, gf in Sigma_imp[0]):
             # Imaginary frequency Sigma:
             self.Sigma_imp = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[icrsh].mesh, space='sumk')
                                  for icrsh in range(self.n_corr_shells)]
             SK_Sigma_imp = self.Sigma_imp
-        elif isinstance(self.mesh, MeshReFreq) and all(isinstance(gf, Gf) and gf.mesh == self.mesh for bname, gf in Sigma_imp[0]):
+        elif isinstance(self.mesh, MeshReFreq) and all(isinstance(gf, Gf) and isinstance(gf.mesh, MeshReFreq) and gf.mesh == self.mesh for bname, gf in Sigma_imp[0]):
             # Real frequency Sigma:
             self.Sigma_imp = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[icrsh].mesh, gf_function=GfReFreq, space='sumk')
                                 for icrsh in range(self.n_corr_shells)]
diff --git a/python/triqs_dft_tools/sumk_dft_tools.py b/python/triqs_dft_tools/sumk_dft_tools.py
index 5f47c03a..77f12738 100644
--- a/python/triqs_dft_tools/sumk_dft_tools.py
+++ b/python/triqs_dft_tools/sumk_dft_tools.py
@@ -822,8 +822,9 @@ class SumkDFTTools(SumkDFT):
               Data as it is also written to the files.
         """
 
-        assert hasattr(
-            self, "Sigma_imp_w"), "spaghettis: Set Sigma_imp_w first."
+        # check if ReFreqMesh is given
+        assert isinstance(self.mesh, MeshReFreq)
+
         things_to_read = ['n_k', 'n_orbitals', 'proj_mat',
                           'hopping', 'n_parproj', 'proj_mat_all']
         value_read = self.read_input_from_hdf(
@@ -840,7 +841,7 @@ class SumkDFTTools(SumkDFT):
         if mu is None:
             mu = self.chemical_potential
         spn = self.spin_block_names[self.SO]
-        mesh = [x.real for x in self.mesh]
+        mesh = numpy.array([x.value for x in self.mesh])
         n_om = len(mesh)
 
         if plot_range is None:
@@ -858,7 +859,7 @@ class SumkDFTTools(SumkDFT):
             Akw = {sp: numpy.zeros(
                 [self.shells[ishell]['dim'], self.n_k, n_om], numpy.float_) for sp in spn}
 
-        if not ishell is None:
+        if ishell is not None:
             gf_struct_parproj = [
                 (sp, self.shells[ishell]['dim']) for sp in spn]
             G_loc = BlockGf(name_block_generator=[(block, GfReFreq(target_shape=(block_dim, block_dim), mesh=self.Sigma_imp_w[0].mesh))
diff --git a/test/python/srvo3_spectral.py b/test/python/srvo3_spectral.py
index 63b833a7..55a47856 100644
--- a/test/python/srvo3_spectral.py
+++ b/test/python/srvo3_spectral.py
@@ -1,20 +1,16 @@
 from triqs_dft_tools.sumk_dft_tools import *
-from triqs.utility.h5diff import h5diff 
+from triqs.utility.h5diff import h5diff
 from h5 import HDFArchive
 
-beta = 40
-
-SK = SumkDFTTools(hdf_file='SrVO3_spectral.h5', use_dft_blocks=True)
-
-if mpi.is_master_node():
-    with HDFArchive('SrVO3_Sigma.h5', 'a') as ar:
+with HDFArchive('SrVO3_Sigma.h5', 'r') as ar:
         Sigma = ar['dmft_transp_input']['Sigma_w']
-        SK.chemical_potential = ar['dmft_transp_input']['chemical_potential']
-        SK.dc_imp = ar['dmft_transp_input']['dc_imp']
+        chemical_potential = ar['dmft_transp_input']['chemical_potential']
+        dc_imp = ar['dmft_transp_input']['dc_imp']
 
-Sigma = mpi.bcast(Sigma)
-SK.chemical_potential = mpi.bcast(SK.chemical_potential)
-SK.dc_imp = mpi.bcast(SK.dc_imp)
+SK = SumkDFTTools(hdf_file='SrVO3_spectral.h5', use_dft_blocks=True, mesh = Sigma.mesh)
+
+SK.chemical_potential = chemical_potential
+SK.dc_imp = dc_imp
 SK.set_Sigma([Sigma])
 
 dos_wannier = SK.dos_wannier_basis(broadening=0.01, with_Sigma=True, with_dc=True, save_to_file=False)
@@ -24,12 +20,12 @@ spaghetti = SK.spaghettis(broadening=0.01, plot_shift=0.0, plot_range=(-1,1), is
 if mpi.is_master_node():
 
     # with HDFArchive('srvo3_spectral.ref.h5', 'a') as ar:
-    #     ar['dos_wannier'] = dos_wannier 
+    #     ar['dos_wannier'] = dos_wannier
     #     ar['dos_parproj'] = dos_parproj
     #     ar['spaghetti'] = spaghetti
 
     with HDFArchive('srvo3_spectral.out.h5', 'a') as ar:
-        ar['dos_wannier'] = dos_wannier 
+        ar['dos_wannier'] = dos_wannier
         ar['dos_parproj'] = dos_parproj
         ar['spaghetti'] = spaghetti
 

From 66b3719c9b1e818edd09f75c10c39d82f4e6d215 Mon Sep 17 00:00:00 2001
From: Alexander Hampel <ahampel@flatironinstitute.org>
Date: Fri, 4 Mar 2022 12:55:28 -0500
Subject: [PATCH 09/11] improve performance of extract Gloc by factor 5! huge
 cleanup

---
 python/triqs_dft_tools/sumk_dft.py | 232 +++++++++++++++--------------
 1 file changed, 123 insertions(+), 109 deletions(-)

diff --git a/python/triqs_dft_tools/sumk_dft.py b/python/triqs_dft_tools/sumk_dft.py
index 72d8cd17..fa3a8c48 100644
--- a/python/triqs_dft_tools/sumk_dft.py
+++ b/python/triqs_dft_tools/sumk_dft.py
@@ -25,7 +25,7 @@ General SumK class and helper functions for combining ab-initio code and triqs
 """
 
 from types import *
-import numpy
+import numpy as np
 import triqs.utility.dichotomy as dichotomy
 from triqs.gf import *
 import triqs.utility.mpi as mpi
@@ -57,7 +57,7 @@ class SumkDFT(object):
                   The value of magnetic field to add to the DFT Hamiltonian.
                   The contribution -h_field*sigma is added to diagonal elements of the Hamiltonian.
                   It cannot be used with the spin-orbit coupling on; namely h_field is set to 0 if self.SO=True.
-        mesh: MeshImFreq or MeshImFreq, optional. Frequency mesh of Sigma.
+        mesh: MeshImFreq or MeshReFreq, optional. Frequency mesh of Sigma.
         beta : real, optional
                Inverse temperature. Used to construct imaginary frequency if mesh is not given.
         n_iw : integer, optional
@@ -105,13 +105,21 @@ class SumkDFT(object):
             self.h_field = h_field
 
             if mesh is None:
-                self.mesh = MeshImFreq(beta=beta, S='Fermion',n_max=n_iw)
-            elif isinstance(mesh, MeshImFreq) or isinstance(mesh, MeshReFreq):
+                self.mesh = MeshImFreq(beta=beta, S='Fermion', n_max=n_iw)
+                self.mesh_values = np.linspace(self.mesh(self.mesh.first_index()),
+                                               self.mesh(self.mesh.last_index()),
+                                               len(self.mesh))
+            elif isinstance(mesh, MeshImFreq):
                 self.mesh = mesh
+                self.mesh_values = np.linspace(self.mesh(self.mesh.first_index()),
+                                               self.mesh(self.mesh.last_index()),
+                                               len(self.mesh))
+            elif isinstance(mesh, MeshReFreq):
+                self.mesh = mesh
+                self.mesh_values = np.linspace(self.mesh.omega_min, self.mesh.omega_max, len(self.mesh))
             else:
                 raise ValueError('mesh must be a triqs mesh of type MeshImFreq or MeshReFreq')
 
-
             self.block_structure = BlockStructure()
 
             # Read input from HDF:
@@ -537,18 +545,32 @@ class SumkDFT(object):
         # Are we including Sigma?
         if with_Sigma:
             Sigma_imp = self.Sigma_imp
-            sigma_minus_dc = [s.copy() for s in Sigma_imp]
             if with_dc:
                 sigma_minus_dc = self.add_dc()
+            else:
+                sigma_minus_dc = [s for s in Sigma_imp]
             if not mesh is None:
                 warn('lattice_gf called with Sigma and given mesh. Mesh will be taken from Sigma.')
-            mesh = Sigma_imp[0].mesh
-            if isinstance(mesh, MeshReFreq) and broadening > 0 and mpi.is_master_node():
-                warn('lattice_gf called with Sigma and broadening > 0 (broadening = {}). You might want to explicitly set the broadening to 0.'.format(broadening))
+
+            if self.mesh != Sigma_imp[0].mesh:
+                mesh = Sigma_imp[0].mesh
+                if isinstance(mesh, MeshImFreq):
+                    mesh_values = np.linspace(mesh(mesh.first_index()), mesh(mesh.last_index()), len(mesh))
+                else:
+                    mesh_values = np.linspace(mesh.omega_min, mesh.omega_max, len(mesh))
+            else:
+                mesh = self.mesh
+                mesh_values = self.mesh_values
+
         elif not mesh is None:
-            assert isinstance(mesh, MreshReFreq) or isinstance(mesh, MeshImFreq),  "mesh must be a triqs MeshReFreq or MeshImFreq"
+            assert isinstance(mesh, MeshReFreq) or isinstance(mesh, MeshImFreq),  "mesh must be a triqs MeshReFreq or MeshImFreq"
+            if isinstance(mesh, MeshImFreq):
+                mesh_values = np.linspace(mesh(mesh.first_index()), mesh(mesh.last_index()), len(mesh))
+            else:
+                mesh_values = np.linspace(mesh.omega_min, mesh.omega_max, len(mesh))
         else:
             mesh = self.mesh
+            mesh_values = self.mesh_values
 
         # Set up G_latt
         if set_up_G_latt:
@@ -567,28 +589,25 @@ class SumkDFT(object):
                              block_list=glist(), make_copies=False)
             G_latt.zero()
 
-        if isinstance(mesh, MeshImFreq):
-            G_latt << iOmega_n
-        else:
-            G_latt << Omega + 1j * broadening
-
-        idmat = [numpy.identity(
-            self.n_orbitals[ik, ntoi[sp]], numpy.complex_) for sp in spn]
-        M = copy.deepcopy(idmat)
-
-        for ibl in range(self.n_spin_blocks[self.SO]):
+        idmat = [np.identity(
+            self.n_orbitals[ik, ntoi[sp]], np.complex_) for sp in spn]
 
+        # fill Glatt
+        for ibl, (block, gf) in enumerate(G_latt):
             ind = ntoi[spn[ibl]]
             n_orb = self.n_orbitals[ik, ind]
-            M[ibl] = self.hopping[ik, ind, 0:n_orb, 0:n_orb] - \
-                (idmat[ibl] * mu) - (idmat[ibl] * self.h_field * (1 - 2 * ibl))
-        G_latt -= M
+            if isinstance(mesh, MeshImFreq):
+                gf.data[:, :, :] = (idmat[ibl] * (mesh_values[:, None, None] + mu + self.h_field*(1-2*ibl))
+                                    - self.hopping[ik, ind, 0:n_orb, 0:n_orb])
+            else:
+                gf.data[:, :, :] = (idmat[ibl] *
+                                    (mesh_values[:, None, None] + mu + self.h_field*(1-2*ibl) + 1j*broadening)
+                                    - self.hopping[ik, ind, 0:n_orb, 0:n_orb])
 
-        if with_Sigma:
-            for icrsh in range(self.n_corr_shells):
-                for bname, gf in G_latt:
-                    gf -= self.upfold(ik, icrsh, bname,
-                                      sigma_minus_dc[icrsh][bname], gf)
+            if with_Sigma:
+                for icrsh in range(self.n_corr_shells):
+                    gf -= self.upfold(ik, icrsh, block,
+                                      sigma_minus_dc[icrsh][block], gf)
 
         G_latt.invert()
         self.G_latt = G_latt
@@ -649,7 +668,7 @@ class SumkDFT(object):
         if isinstance(self.mesh, MeshReFreq):
             if self.min_band_energy is None or self.max_band_energy is None:
                 self.calculate_min_max_band_energies()
-            mesh = numpy.array([i for i in self.mesh.values()])
+            mesh = np.array([i for i in self.mesh.values()])
             if mesh[0] > (self.min_band_energy - self.chemical_potential) or mesh[-1] < (self.max_band_energy - self.chemical_potential):
                 warn('The given Sigma is on a mesh which does not cover the band energy range. The Sigma MeshReFreq runs from %f to %f, while the band energy (minus the chemical potential) runs from %f to %f'%(mesh[0], mesh[-1], self.min_band_energy, self.max_band_energy))
 
@@ -754,7 +773,7 @@ class SumkDFT(object):
         for icrsh in range(self.n_corr_shells):
             G_loc[icrsh].zero()                          # initialize to zero
 
-        ikarray = numpy.array(list(range(self.n_k)))
+        ikarray = np.array(list(range(self.n_k)))
         for ik in mpi.slice_array(ikarray):
             if isinstance(self.mesh, MeshImFreq):
                 G_latt = self.lattice_gf(
@@ -766,11 +785,8 @@ class SumkDFT(object):
 
             for icrsh in range(self.n_corr_shells):
                 # init temporary storage
-                tmp = G_loc[icrsh].copy()
-                for bname, gf in tmp:
-                    tmp[bname] << self.downfold(
-                        ik, icrsh, bname, G_latt[bname], gf)
-                G_loc[icrsh] += tmp
+                for bname, gf in G_loc[icrsh]:
+                    gf += self.downfold(ik, icrsh, bname, G_latt[bname], gf)
 
         # Collect data from mpi
         for icrsh in range(self.n_corr_shells):
@@ -930,8 +946,8 @@ class SumkDFT(object):
             dm = {}
             for block, block_dim in self.gf_struct_solver[ish].items():
                 # get dm for the blocks:
-                dm[block] = numpy.zeros(
-                    [block_dim, block_dim], numpy.complex_)
+                dm[block] = np.zeros(
+                    [block_dim, block_dim], complex)
                 for ind1 in range(block_dim):
                     for ind2 in range(block_dim):
                         block_sumk, ind1_sumk = self.solver_to_sumk[
@@ -992,7 +1008,7 @@ class SumkDFT(object):
             gf = [g_sh.copy() for g_sh in G]
             for ish in range(len(gf)):
                 for name, g in gf[ish]:
-                    g << 1.0j*(g-g.conjugate().transpose())/2.0/numpy.pi
+                    g << 1.0j*(g-g.conjugate().transpose())/2.0/np.pi
         elif all(isinstance(g_sh.mesh, MeshReTime) for g_sh in G):
             def get_delta_from_mesh(mesh):
                 w0 = None
@@ -1002,15 +1018,15 @@ class SumkDFT(object):
                     else:
                         return w-w0
             gf = [BlockGf(name_block_generator = [(name, GfReFreq(
-                window=(-numpy.pi*(len(block.mesh)-1) / (len(block.mesh)*get_delta_from_mesh(block.mesh)),
-                numpy.pi*(len(block.mesh)-1) / (len(block.mesh)*get_delta_from_mesh(block.mesh))),
+                window=(-np.pi*(len(block.mesh)-1) / (len(block.mesh)*get_delta_from_mesh(block.mesh)),
+                np.pi*(len(block.mesh)-1) / (len(block.mesh)*get_delta_from_mesh(block.mesh))),
                 n_points=len(block.mesh), indices=block.indices)) for name, block in g_sh], make_copies=False)
                 for g_sh in G]
 
             for ish in range(len(gf)):
                 for name, g in gf[ish]:
                     g.set_from_fourier(G[ish][name])
-                    g << 1.0j*(g-g.conjugate().transpose())/2.0/numpy.pi
+                    g << 1.0j*(g-g.conjugate().transpose())/2.0/np.pi
         else:
             raise Exception("G must be a list of BlockGf of either GfImFreq, GfImTime, GfReFreq or GfReTime")
         return gf
@@ -1061,7 +1077,7 @@ class SumkDFT(object):
                                      for ish in range(self.n_inequiv_shells)]
 
         # the maximum value of each matrix element of each block and shell
-        max_gf = [{name:numpy.max(numpy.abs(g.data),0) for name, g in gf[ish]} for ish in range(self.n_inequiv_shells)]
+        max_gf = [{name:np.max(np.abs(g.data),0) for name, g in gf[ish]} for ish in range(self.n_inequiv_shells)]
 
         if include_shells is None:
             # include all shells
@@ -1160,7 +1176,7 @@ class SumkDFT(object):
         # helper function
         def null(A, eps=1e-15):
             """ Calculate the null-space of matrix A """
-            u, s, vh = numpy.linalg.svd(A)
+            u, s, vh = np.linalg.svd(A)
             null_mask = (s <= eps)
             null_space = compress(null_mask, vh, axis=0)
             return null_space.conjugate().transpose()
@@ -1215,9 +1231,9 @@ class SumkDFT(object):
                     # one block to make it equal to the other (at least
                     # for tau=0).
 
-                    e1 = numpy.linalg.eigvalsh(gf1.data[0])
-                    e2 = numpy.linalg.eigvalsh(gf2.data[0])
-                    if numpy.any(abs(e1-e2) > threshold): continue
+                    e1 = np.linalg.eigvalsh(gf1.data[0])
+                    e2 = np.linalg.eigvalsh(gf2.data[0])
+                    if np.any(abs(e1-e2) > threshold): continue
 
                     for conjugate in [False,True]:
                         if conjugate:
@@ -1236,21 +1252,21 @@ class SumkDFT(object):
                         # product to get a linear problem, which consists
                         # of finding the null space of M vec T = 0.
 
-                        M = numpy.kron(numpy.eye(*gf1.target_shape),gf2.data[0])-numpy.kron(gf1.data[0].transpose(),numpy.eye(*gf1.target_shape))
+                        M = np.kron(np.eye(*gf1.target_shape),gf2.data[0])-np.kron(gf1.data[0].transpose(),np.eye(*gf1.target_shape))
                         N = null(M, threshold)
 
                         # now we get the intersection of the null spaces
                         # of all values of tau
                         for i in range(1,len(gf1.data)):
-                            M = numpy.kron(numpy.eye(*gf1.target_shape),gf2.data[i])-numpy.kron(gf1.data[i].transpose(),numpy.eye(*gf1.target_shape))
+                            M = np.kron(np.eye(*gf1.target_shape),gf2.data[i])-np.kron(gf1.data[i].transpose(),np.eye(*gf1.target_shape))
                             # transform M into current null space
-                            M = numpy.dot(M, N)
-                            N = numpy.dot(N, null(M, threshold))
-                            if numpy.size(N) == 0:
+                            M = np.dot(M, N)
+                            N = np.dot(N, null(M, threshold))
+                            if np.size(N) == 0:
                                 break
 
                         # no intersection of the null spaces -> no symmetry
-                        if numpy.size(N) == 0: continue
+                        if np.size(N) == 0: continue
 
                         # reshape N: it then has the indices matrix, matrix, number of basis vectors of the null space
                         N = N.reshape(gf1.target_shape[0], gf1.target_shape[1], -1).transpose([1, 0, 2])
@@ -1268,7 +1284,7 @@ class SumkDFT(object):
                             sum  y[i] N[:,:,i] y[j].conjugate() N[:,:,j].conjugate().transpose() = eye.
                         The object N[:,:,i] N[:,:,j] is a four-index object which we call Z.
                         """
-                        Z = numpy.einsum('aci,bcj->abij', N, N.conjugate())
+                        Z = np.einsum('aci,bcj->abij', N, N.conjugate())
 
                         """
                         function chi2
@@ -1278,16 +1294,16 @@ class SumkDFT(object):
                         """
                         def chi2(y):
                             # reinterpret y as complex number
-                            y = y.view(numpy.complex_)
+                            y = y.view(np.complex_)
                             ret = 0.0
                             for a in range(Z.shape[0]):
                                 for b in range(Z.shape[1]):
-                                    ret += numpy.abs(numpy.dot(y, numpy.dot(Z[a, b], y.conjugate()))
+                                    ret += np.abs(np.dot(y, np.dot(Z[a, b], y.conjugate()))
                                                   - (1.0 if a == b else 0.0))**2
                             return ret
 
                         # use the minimization routine from scipy
-                        res = minimize(chi2, numpy.ones(2 * N.shape[-1]))
+                        res = minimize(chi2, np.ones(2 * N.shape[-1]))
 
                         # if the minimization fails, there is probably no symmetry
                         if not res.success: continue
@@ -1295,10 +1311,10 @@ class SumkDFT(object):
                         if res.fun > threshold: continue
 
                         # reinterpret the solution as a complex number
-                        y = res.x.view(numpy.complex_)
+                        y = res.x.view(np.complex_)
 
                         # reconstruct the T matrix
-                        T = numpy.zeros(N.shape[:-1], dtype=numpy.complex_)
+                        T = np.zeros(N.shape[:-1], dtype=np.complex_)
                         for i in range(len(y)):
                             T += N[:, :, i] * y[i]
 
@@ -1336,9 +1352,9 @@ class SumkDFT(object):
                         # C1(v1^dagger G1 v1) = C2(v2^dagger G2 v2)
                         if (ind1 < 0) and (ind2 >= 0):
                             if conjugate:
-                                self.deg_shells[ish][ind2][block1] = numpy.dot(T.conjugate().transpose(), v2[0].conjugate()), not v2[1]
+                                self.deg_shells[ish][ind2][block1] = np.dot(T.conjugate().transpose(), v2[0].conjugate()), not v2[1]
                             else:
-                                self.deg_shells[ish][ind2][block1] = numpy.dot(T.conjugate().transpose(), v2[0]), v2[1]
+                                self.deg_shells[ish][ind2][block1] = np.dot(T.conjugate().transpose(), v2[0]), v2[1]
                         # the first block is already present
                         # set v2 and C2 so that they are compatible with
                         # C(T gf1 T^dagger) = gf2
@@ -1346,9 +1362,9 @@ class SumkDFT(object):
                         # C1(v1^dagger G1 v1) = C2(v2^dagger G2 v2)
                         elif (ind1 >= 0) and (ind2 < 0):
                             if conjugate:
-                                self.deg_shells[ish][ind1][block2] = numpy.dot(T.conjugate(), v1[0].conjugate()), not v1[1]
+                                self.deg_shells[ish][ind1][block2] = np.dot(T.conjugate(), v1[0].conjugate()), not v1[1]
                             else:
-                                self.deg_shells[ish][ind1][block2] = numpy.dot(T, v1[0]), v1[1]
+                                self.deg_shells[ish][ind1][block2] = np.dot(T, v1[0]), v1[1]
                         # the blocks are not already present
                         # we arbitrarily choose v1=eye and C1=False and
                         # set v2 and C2 so that they are compatible with
@@ -1357,7 +1373,7 @@ class SumkDFT(object):
                         # C1(v1^dagger G1 v1) = C2(v2^dagger G2 v2)
                         elif (ind1 < 0) and (ind2 < 0):
                             d = dict()
-                            d[block1] = numpy.eye(*gf1.target_shape), False
+                            d[block1] = np.eye(*gf1.target_shape), False
                             if conjugate:
                                 d[block2] = T.conjugate(), True
                             else:
@@ -1418,7 +1434,7 @@ class SumkDFT(object):
                 "calculate_diagonalization_matrix: Choices for prop_to_be_diagonal are 'eal' or 'dm'.")
             return 0
 
-        trans = [{block: numpy.eye(block_dim) for block, block_dim in gfs} for gfs in self.gf_struct_sumk]
+        trans = [{block: np.eye(block_dim) for block, block_dim in gfs} for gfs in self.gf_struct_sumk]
 
         for ish in shells:
             trafo = {}
@@ -1427,10 +1443,10 @@ class SumkDFT(object):
                 prop[ish] = self.block_structure.convert_matrix(prop[ish], space_from='sumk', space_to='solver')
             # Get diagonalisation matrix, if not already diagonal
             for name in prop[ish]:
-                if numpy.sum(abs(prop[ish][name]-numpy.diag(numpy.diagonal(prop[ish][name])))) > 1e-13:
-                    trafo[name] = numpy.linalg.eigh(prop[ish][name])[1].conj().T
+                if np.sum(abs(prop[ish][name]-np.diag(np.diagonal(prop[ish][name])))) > 1e-13:
+                    trafo[name] = np.linalg.eigh(prop[ish][name])[1].conj().T
                 else:
-                    trafo[name] = numpy.identity(numpy.shape(prop[ish][name])[0])
+                    trafo[name] = np.identity(np.shape(prop[ish][name])[0])
             # Transform back to sumk if necessay, blocks change in this step!
             if calc_in_solver_blocks:
                 trafo = self.block_structure.convert_matrix(trafo, space_from='solver', space_to='sumk')
@@ -1467,10 +1483,10 @@ class SumkDFT(object):
         dens_mat = [{} for icrsh in range(self.n_corr_shells)]
         for icrsh in range(self.n_corr_shells):
             for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]:
-                dens_mat[icrsh][sp] = numpy.zeros(
-                    [self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], numpy.complex_)
+                dens_mat[icrsh][sp] = np.zeros(
+                    [self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], np.complex_)
 
-        ikarray = numpy.array(list(range(self.n_k)))
+        ikarray = np.array(list(range(self.n_k)))
         for ik in mpi.slice_array(ikarray):
 
             if method == "using_gf":
@@ -1485,7 +1501,7 @@ class SumkDFT(object):
                 ntoi = self.spin_names_to_ind[self.SO]
                 spn = self.spin_block_names[self.SO]
                 dims = {sp:self.n_orbitals[ik, ntoi[sp]] for sp in spn}
-                MMat = [numpy.zeros([dims[sp], dims[sp]], numpy.complex_) for sp in spn]
+                MMat = [np.zeros([dims[sp], dims[sp]], np.complex_) for sp in spn]
 
                 for isp, sp in enumerate(spn):
                     ind = ntoi[sp]
@@ -1508,10 +1524,10 @@ class SumkDFT(object):
                     n_orb = self.n_orbitals[ik, ind]
                     projmat = self.proj_mat[ik, ind, icrsh, 0:dim, 0:n_orb]
                     if method == "using_gf":
-                        dens_mat[icrsh][sp] += numpy.dot(numpy.dot(projmat, MMat[isp]),
+                        dens_mat[icrsh][sp] += np.dot(np.dot(projmat, MMat[isp]),
                                                          projmat.transpose().conjugate())
                     elif method == "using_point_integration":
-                        dens_mat[icrsh][sp] += self.bz_weights[ik] * numpy.dot(numpy.dot(projmat, MMat[isp]),
+                        dens_mat[icrsh][sp] += self.bz_weights[ik] * np.dot(np.dot(projmat, MMat[isp]),
                                                                                projmat.transpose().conjugate())
 
         # get data from nodes:
@@ -1530,7 +1546,7 @@ class SumkDFT(object):
                 for sp in dens_mat[icrsh]:
                     if self.rot_mat_time_inv[icrsh] == 1:
                         dens_mat[icrsh][sp] = dens_mat[icrsh][sp].conjugate()
-                    dens_mat[icrsh][sp] = numpy.dot(numpy.dot(self.rot_mat[icrsh].conjugate().transpose(), dens_mat[icrsh][sp]),
+                    dens_mat[icrsh][sp] = np.dot(np.dot(self.rot_mat[icrsh].conjugate().transpose(), dens_mat[icrsh][sp]),
                                                     self.rot_mat[icrsh])
 
         return dens_mat
@@ -1565,8 +1581,8 @@ class SumkDFT(object):
         eff_atlevels = [{} for ish in range(self.n_inequiv_shells)]
         for ish in range(self.n_inequiv_shells):
             for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]:
-                eff_atlevels[ish][sp] = numpy.identity(
-                    self.corr_shells[self.inequiv_to_corr[ish]]['dim'], numpy.complex_)
+                eff_atlevels[ish][sp] = np.identity(
+                    self.corr_shells[self.inequiv_to_corr[ish]]['dim'], np.complex_)
                 eff_atlevels[ish][sp] *= -self.chemical_potential
                 eff_atlevels[ish][
                     sp] -= self.dc_imp[self.inequiv_to_corr[ish]][sp]
@@ -1579,18 +1595,18 @@ class SumkDFT(object):
             for icrsh in range(self.n_corr_shells):
                 dim = self.corr_shells[icrsh]['dim']
                 for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]:
-                    self.Hsumk[icrsh][sp] = numpy.zeros(
-                        [dim, dim], numpy.complex_)
+                    self.Hsumk[icrsh][sp] = np.zeros(
+                        [dim, dim], np.complex_)
                 for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh]['SO']]):
                     ind = self.spin_names_to_ind[
                         self.corr_shells[icrsh]['SO']][sp]
                     for ik in range(self.n_k):
                         n_orb = self.n_orbitals[ik, ind]
-                        MMat = numpy.identity(n_orb, numpy.complex_)
+                        MMat = np.identity(n_orb, np.complex_)
                         MMat = self.hopping[
                             ik, ind, 0:n_orb, 0:n_orb] - (1 - 2 * isp) * self.h_field * MMat
                         projmat = self.proj_mat[ik, ind, icrsh, 0:dim, 0:n_orb]
-                        self.Hsumk[icrsh][sp] += self.bz_weights[ik] * numpy.dot(numpy.dot(projmat, MMat),
+                        self.Hsumk[icrsh][sp] += self.bz_weights[ik] * np.dot(np.dot(projmat, MMat),
                                                                                  projmat.conjugate().transpose())
             # symmetrisation:
             if self.symm_op != 0:
@@ -1603,7 +1619,7 @@ class SumkDFT(object):
                         if self.rot_mat_time_inv[icrsh] == 1:
                             self.Hsumk[icrsh][sp] = self.Hsumk[
                                 icrsh][sp].conjugate()
-                        self.Hsumk[icrsh][sp] = numpy.dot(numpy.dot(self.rot_mat[icrsh].conjugate().transpose(), self.Hsumk[icrsh][sp]),
+                        self.Hsumk[icrsh][sp] = np.dot(np.dot(self.rot_mat[icrsh].conjugate().transpose(), self.Hsumk[icrsh][sp]),
                                                           self.rot_mat[icrsh])
 
         # add to matrix:
@@ -1628,7 +1644,7 @@ class SumkDFT(object):
             dim = self.corr_shells[icrsh]['dim']
             spn = self.spin_block_names[self.corr_shells[icrsh]['SO']]
             for sp in spn:
-                self.dc_imp[icrsh][sp] = numpy.zeros([dim, dim], numpy.float_)
+                self.dc_imp[icrsh][sp] = np.zeros([dim, dim], np.float_)
         self.dc_energ = [0.0 for icrsh in range(self.n_corr_shells)]
 
     def set_dc(self, dc_imp, dc_energ):
@@ -1706,7 +1722,7 @@ class SumkDFT(object):
                 Ncr[bl] += dens_mat[block].real.trace()
             Ncrtot = sum(Ncr.values())
             for sp in spn:
-                self.dc_imp[icrsh][sp] = numpy.identity(dim, numpy.float_)
+                self.dc_imp[icrsh][sp] = np.identity(dim, np.float_)
                 if self.SP == 0:  # average the densities if there is no SP:
                     Ncr[sp] = Ncrtot / len(spn)
                 # correction for SO: we have only one block in this case, but
@@ -1771,8 +1787,8 @@ class SumkDFT(object):
             if transform:
                 for sp in spn:
                     T = self.block_structure.effective_transformation_sumk[icrsh][sp]
-                    self.dc_imp[icrsh][sp] = numpy.dot(T.conjugate().transpose(),
-                            numpy.dot(self.dc_imp[icrsh][sp], T))
+                    self.dc_imp[icrsh][sp] = np.dot(T.conjugate().transpose(),
+                            np.dot(self.dc_imp[icrsh][sp], T))
 
     def add_dc(self):
         r"""
@@ -1798,12 +1814,10 @@ class SumkDFT(object):
             for bname, gf in sigma_minus_dc[icrsh]:
                 # Transform dc_imp to global coordinate system
                 if self.use_rotations:
-                    dccont = numpy.dot(self.rot_mat[icrsh], numpy.dot(self.dc_imp[icrsh][
-                                    bname], self.rot_mat[icrsh].conjugate().transpose()))
+                    gf -= np.dot(self.rot_mat[icrsh], np.dot(self.dc_imp[icrsh][
+                                       bname], self.rot_mat[icrsh].conjugate().transpose()))
                 else:
-                    dccont = self.dc_imp[icrsh][bname]
-                    
-                sigma_minus_dc[icrsh][bname] -= dccont
+                    gf -= self.dc_imp[icrsh][bname]
 
         return sigma_minus_dc
 
@@ -1844,7 +1858,7 @@ class SumkDFT(object):
                     v, C = degsh[key]
                 else:
                     # for backward compatibility, allow degsh to be a list
-                    v = numpy.eye(*ss.target_shape)
+                    v = np.eye(*ss.target_shape)
                     C = False
                 # the helper is in the basis where the blocks are all equal
                 helper.from_L_G_R(v.conjugate().transpose(), gf_to_symm[key], v)
@@ -1858,7 +1872,7 @@ class SumkDFT(object):
                     v, C = degsh[key]
                 else:
                     # for backward compatibility, allow degsh to be a list
-                    v = numpy.eye(*ss.target_shape)
+                    v = np.eye(*ss.target_shape)
                     C = False
                 if C:
                     gf_to_symm[key].from_L_G_R(v, ss.transpose().copy(), v.conjugate().transpose())
@@ -1915,7 +1929,7 @@ class SumkDFT(object):
         if mu is None:
             mu = self.chemical_potential
         dens = 0.0
-        ikarray = numpy.array(list(range(self.n_k)))
+        ikarray = np.array(list(range(self.n_k)))
         for ik in mpi.slice_array(ikarray):
             G_latt = self.lattice_gf(
                 ik=ik, mu=mu, with_Sigma=with_Sigma, with_dc=with_dc, broadening=broadening)
@@ -2046,16 +2060,16 @@ class SumkDFT(object):
 # Convert Fermi weights to a density matrix
             dens_mat_dft = {}
             for sp in spn:
-                dens_mat_dft[sp] = [fermi_weights[ik, ntoi[sp], :].astype(numpy.complex_) for ik in range(self.n_k)]
+                dens_mat_dft[sp] = [fermi_weights[ik, ntoi[sp], :].astype(np.complex_) for ik in range(self.n_k)]
 
 
         # Set up deltaN:
         deltaN = {}
         for sp in spn:
-            deltaN[sp] = [numpy.zeros([self.n_orbitals[ik, ntoi[sp]], self.n_orbitals[
-                                      ik, ntoi[sp]]], numpy.complex_) for ik in range(self.n_k)]
+            deltaN[sp] = [np.zeros([self.n_orbitals[ik, ntoi[sp]], self.n_orbitals[
+                                      ik, ntoi[sp]]], np.complex_) for ik in range(self.n_k)]
 
-        ikarray = numpy.arange(self.n_k)
+        ikarray = np.arange(self.n_k)
         for ik in mpi.slice_array(ikarray):
             G_latt_iw = self.lattice_gf(
                 ik=ik, mu=self.chemical_potential, iw_or_w="iw")
@@ -2075,11 +2089,11 @@ class SumkDFT(object):
                 if dm_type in ['vasp','qe']:
 # In 'vasp'-mode subtract the DFT density matrix
                     nb = self.n_orbitals[ik, ntoi[bname]]
-                    diag_inds = numpy.diag_indices(nb)
+                    diag_inds = np.diag_indices(nb)
                     deltaN[bname][ik][diag_inds] -= dens_mat_dft[bname][ik][:nb]
 
                     if self.charge_mixing and self.deltaNOld is not None:
-                        G2 = numpy.sum(self.kpts_cart[ik,:]**2)
+                        G2 = np.sum(self.kpts_cart[ik,:]**2)
                         # Kerker mixing
                         mix_fac = self.charge_mixing_alpha * G2 / (G2 + self.charge_mixing_gamma**2)
                         deltaN[bname][ik][diag_inds] = (1.0 - mix_fac) * self.deltaNOld[bname][ik][diag_inds] + mix_fac * deltaN[bname][ik][diag_inds]
@@ -2088,7 +2102,7 @@ class SumkDFT(object):
                     b1, b2 = band_window[isp][ik, :2]
                     nb = b2 - b1 + 1
                     assert nb == self.n_orbitals[ik, ntoi[bname]], "Number of bands is inconsistent at ik = %s"%(ik)
-                    band_en_correction += numpy.dot(deltaN[bname][ik], self.hopping[ik, isp, :nb, :nb]).trace().real * self.bz_weights[ik]
+                    band_en_correction += np.dot(deltaN[bname][ik], self.hopping[ik, isp, :nb, :nb]).trace().real * self.bz_weights[ik]
 
         # mpi reduce:
         for bname in deltaN:
@@ -2157,9 +2171,9 @@ class SumkDFT(object):
                         fout.close()
         elif dm_type == 'vasp':
             if kpts_to_write is None:
-                kpts_to_write = numpy.arange(self.n_k)
+                kpts_to_write = np.arange(self.n_k)
             else:
-                assert numpy.min(kpts_to_write) >= 0 and numpy.max(kpts_to_write) < self.n_k
+                assert np.min(kpts_to_write) >= 0 and np.max(kpts_to_write) < self.n_k
 
             assert self.SP == 0, "Spin-polarized density matrix is not implemented"
 
@@ -2188,7 +2202,7 @@ class SumkDFT(object):
                 with open(filename, 'w') as f:
         # determine the number of spin blocks
                     n_spin_blocks = self.SP + 1 - self.SO
-                    nbmax = numpy.max(self.n_orbitals)
+                    nbmax = np.max(self.n_orbitals)
         # output beta and mu in Hartrees
                     beta = G_latt_iw.mesh.beta * self.energy_unit
                     mu = self.chemical_potential/self.energy_unit
@@ -2223,7 +2237,7 @@ class SumkDFT(object):
             assert self.SP == 0, "Spin-polarized density matrix is not implemented"
 
             subgrp = 'dft_update'
-            delta_N = numpy.zeros([self.n_k, max(self.n_orbitals[:,0]), max(self.n_orbitals[:,0])], dtype=complex)
+            delta_N = np.zeros([self.n_k, max(self.n_orbitals[:,0]), max(self.n_orbitals[:,0])], dtype=complex)
             mpi.report(" %i  -1  ! Number of k-points, default number of bands\n"%(self.n_k))
             for ik in range(self.n_k):
                 ib1 = band_window[0][ik, 0]
@@ -2255,10 +2269,10 @@ class SumkDFT(object):
 
     def calculate_min_max_band_energies(self):
         hop = self.hopping
-        diag_hop = numpy.zeros(hop.shape[:-1])
+        diag_hop = np.zeros(hop.shape[:-1])
         hop_slice = mpi.slice_array(hop)
         diag_hop_slice = mpi.slice_array(diag_hop)
-        diag_hop_slice[:] = numpy.linalg.eigvalsh(hop_slice)
+        diag_hop_slice[:] = np.linalg.eigvalsh(hop_slice)
         diag_hop = mpi.all_reduce(mpi.world, diag_hop, lambda x, y: x + y)
         min_band_energy = diag_hop.min().real
         max_band_energy = diag_hop.max().real
@@ -2297,7 +2311,7 @@ class SumkDFT(object):
     def check_projectors(self):
         """Calculated the density matrix from projectors (DM = P Pdagger) to check that it is correct and
            specifically that it matches DFT."""
-        dens_mat = [numpy.zeros([self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], numpy.complex_)
+        dens_mat = [np.zeros([self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], np.complex_)
                     for icrsh in range(self.n_corr_shells)]
 
         for ik in range(self.n_k):
@@ -2306,7 +2320,7 @@ class SumkDFT(object):
                 n_orb = self.n_orbitals[ik, 0]
                 projmat = self.proj_mat[ik, 0, icrsh, 0:dim, 0:n_orb]
                 dens_mat[icrsh][
-                    :, :] += numpy.dot(projmat, projmat.transpose().conjugate()) * self.bz_weights[ik]
+                    :, :] += np.dot(projmat, projmat.transpose().conjugate()) * self.bz_weights[ik]
 
         if self.symm_op != 0:
             dens_mat = self.symmcorr.symmetrize(dens_mat)
@@ -2316,7 +2330,7 @@ class SumkDFT(object):
             for icrsh in range(self.n_corr_shells):
                 if self.rot_mat_time_inv[icrsh] == 1:
                     dens_mat[icrsh] = dens_mat[icrsh].conjugate()
-                dens_mat[icrsh] = numpy.dot(numpy.dot(self.rot_mat[icrsh].conjugate().transpose(), dens_mat[icrsh]),
+                dens_mat[icrsh] = np.dot(np.dot(self.rot_mat[icrsh].conjugate().transpose(), dens_mat[icrsh]),
                                             self.rot_mat[icrsh])
 
         return dens_mat

From 213016ee1c099b35b922d9bbf5f881d98eb09e33 Mon Sep 17 00:00:00 2001
From: Alexander Hampel <ahampel@flatironinstitute.org>
Date: Fri, 15 Apr 2022 10:31:14 -0400
Subject: [PATCH 10/11] fix transport test, allow for parallel execution of
 tests

---
 python/triqs_dft_tools/sumk_dft_tools.py |   4 ++--
 test/python/SrVO3_Sigma_transport.h5     | Bin 0 -> 318032 bytes
 test/python/srvo3_transp.py              |   7 ++++---
 3 files changed, 6 insertions(+), 5 deletions(-)
 create mode 100644 test/python/SrVO3_Sigma_transport.h5

diff --git a/python/triqs_dft_tools/sumk_dft_tools.py b/python/triqs_dft_tools/sumk_dft_tools.py
index 77f12738..9cf955f2 100644
--- a/python/triqs_dft_tools/sumk_dft_tools.py
+++ b/python/triqs_dft_tools/sumk_dft_tools.py
@@ -862,7 +862,7 @@ class SumkDFTTools(SumkDFT):
         if ishell is not None:
             gf_struct_parproj = [
                 (sp, self.shells[ishell]['dim']) for sp in spn]
-            G_loc = BlockGf(name_block_generator=[(block, GfReFreq(target_shape=(block_dim, block_dim), mesh=self.Sigma_imp_w[0].mesh))
+            G_loc = BlockGf(name_block_generator=[(block, GfReFreq(target_shape=(block_dim, block_dim), mesh=self.Sigma_imp[0].mesh))
                                                   for block, block_dim in gf_struct_parproj], make_copies=False)
             G_loc.zero()
 
@@ -1207,7 +1207,7 @@ class SumkDFTTools(SumkDFT):
                     spn = self.spin_block_names[self.corr_shells[icrsh]['SO']]
                     glist = lambda: [GfReFreq(target_shape=(block_dim, block_dim), window=(self.omega[
                                               0], self.omega[-1]), n_points=n_om) for block, block_dim in self.gf_struct_sumk[icrsh]]
-                    self.Sigma_imp_w[icrsh] = BlockGf(
+                    self.Sigma_imp[icrsh] = BlockGf(
                         name_list=spn, block_list=glist(), make_copies=False)
                     for i, g in self.Sigma_imp[icrsh]:
                         for iL in g.indices[0]:
diff --git a/test/python/SrVO3_Sigma_transport.h5 b/test/python/SrVO3_Sigma_transport.h5
new file mode 100644
index 0000000000000000000000000000000000000000..616aafcc0f68106cc27d9937940efbd5992cdfd3
GIT binary patch
literal 318032
zcmeEP2V7IT_rC=d+^A2*@zc6-S6n!A-MCRvae@QEfmQ{@jT?8J;CkZ5fq)Zd>c))&
z#eokM1Q((rd;D9H=)bAAMISs`2!3pCa&zy^IXNdMIr%18X4$cG(L&`5amqhqV=f<O
zl5t6X((#|wvwRszBu?(z;R_KC-Ert`fG$Lw@*l^U;`9>u{d7L-zMVR8ik#&AB*c$#
zVfoROF-c%yp5lOIr*>9Kf^9I>5*!9%Im1!TdCG{1Q>BxpIOBxTPGg*HQKIeGiBcDm
zPUO{$JSY$f$*-j$dcYCOr8fCBjY9dE{0cBSmnYNmOIlfXEI)+&eq$Y+M%YfzXfu^*
zmgL4ChW&c=?kM66ezd=K^N`>Y`FF>0@A)X5$fp>~H_1oK`yUoZZyytC2OKfOtc(z8
zbUhP(O7VTd0K(UdMvWxq!mt3Q3AK-$pnZ2oJrIuYD;KhYk)(kGBF<2ypJ;;niDD7b
z6k_+a*;%U+=p#&yG%!0C-=f=JOH6A_5mlmHq_*bD6h)ZPl~ja{l<|=hCXE_zJ8^`Q
z{ZwU?8$D_IL}kv5P>~}eR6If2Z(_%hofZ9PzesyV^7-NREr=>d#z)xC4-N>m1D+*O
zmdv&T#($ywq<(&UfL-JX4$$pDhLS5NPoM;1^a$q>N`_1wJJG>q!iXtjoo8h5N1^|Y
z%8RiRM~@w)XbXrQn`PV;{z2`45huo3qE-m!Ro>z8ZO3Qmf_z7Y<9;ttx&_YXcT6>3
z2#zP5K<Pi|T&bE*yb8tnH7F!(aBZua&kM)pJW4l30VlCn&1Zw-t_~=j*jrgK)qDXs
z9zs5d-={4gej@t$A{^p`h%2jTYmGVUzFqtDvsE~;n%XL>(st@7M|&ql*-9xGsVo!W
zXO}U3?K@AgpG>~1tsmjI0y!24UF10yG7f0sfc4yokGq?Sc!&AcTT*=^;bLhtEMq3(
zeQQ2_A?p+gUjXlGx0{Lh_T%4+U7x*$EsIndKBc3GU)6cZ`<D~n!b>YU4Jy({#7he{
z8yRiz7OtGvI&{ZS5r4-ie88Pe5%AkmKWlRb5%0HeLI<xp5%6!`aP$qSh#zxd;=pAW
z!eR3cl9%IVi1?oMk`{RP4Tl!jrk|-eU&LQ&uy;{~$T#r1=x(8ir6T^>O`pWnS#RKq
zenr3CSS8|PB!QF0mUsi7wu|(-zfQzkNISW@?0*e?zpw6S;348YeLJ3O-u5+|v$#y5
zky}Lk(J9^~t^~h=ZKq5x6~9fyi?+&S^(MZ88&cX#z5cg|FLFnGV7}2Scz^QUT2FU~
z_%};Nw0N*H4AyAVs;s$}h<99Svb%1ZFj#cxuEHUT`}fLEyVUsDKXAz%W4|E9^V8oM
zw(1@I5+3zFxT?x_5x-~C(`D!CzJ%vboGcim`2L}K<1JQ>4~5hBP8jt{QNFbCar<uj
zynww1Y&<%6gNWbO#?yTH*XQu$*B)<qnTR(|PZ@Nt#dA2TLfI>xD@1&)=Uc<}EPMt>
zpKlsDW{HTm^qvw^BIqd$*x5?%>n7sIk8IttV&$iBaLY#<rq2}d$u{D@0;M5v@wPMb
zE=^U`-@50UOs_tH?}|K=ZkQnA<@X&=7q9pP7WlNT{+H1resYxt<BzyJhI=|z?(;%X
zf09c9hhy&t!$nI<_8ep_;@cE0x&M_U7@k-gGj@D;5pP!W^u?*`9>Lx{Lavu?FXAr@
z?RCmK?jdZsqOI)^3lZ-Wk?K8X*hARA$4L=rAmRs<50HGm_5k{<zdZPIH4*QBuhW!7
z^9OMFqLVT274r9Oa+kWEJMY6Ct5((sC@SJN1k}HDt>k@p@=|c~`UWC?o<~b-w<SR^
z{ou9^J>ociNJ68<^Ahht_d>lplzq?fkE=JHlm7lLtlDp%i7ZrMrpk6Y=<@0gbnsYh
z9}vXx|9c0&K7Mf<w&n_luD#6h1=cV8pK0hVSovL<*#Dj6`0k0`;q_nNguUJ_C^X!g
z<2y~_jEa7|0o#?nGpgSfj+eD|u{oF&2>)pmTPuA9$B(-<;eJxlKsZm-Y4pfB9PeD9
z`Sx!$uERZNsuY<!k>meP+f!?K`)hDk)ys)}Z8`qnw~*S$M_q-sb-qQv>do<6yNo*N
zu<Qza{#L%RF3<7X-=FWe_1I<D<$C@35(&pQUDW(sr`MO@ZSGP1FO@ld!|vg8hnKzt
z57%8{FtIqtKcAN(4efCeikk)P?3J1hHk%$XI_oZnAJavLTRu+*?nnA<30n(c@S9aG
zUIFRgh{N2FPvtMbQAu5DyKPJd8yB>kyzcNhcnH3`^>lJNXk6s!#4DZ8!qAR}ZG3vB
zgD@jm?^*9o!|IE|nkSj1gLfC&mGRpU02_84Wi2;M2Zh`k4jE>73O=j<z$iU54a{xj
zGjv|QlW?Q?_Pr-gr2)~a_r?FZaU4!qVrutBmIj<YMK-Fi=@{Hou0hwOW7B|fmAy8b
z#vg^XcNOU8)h-QqSM|=ftgAm<E<b2hzg!x4vo!o@LPI~;Eqz;Dy_i(+q-(2}Pum=U
z7E-6@+pnbpafJ7sm-P<AaVNUC&fS&@g6|$kEL;B|ywPiR?K_iGL09XKXLh#rh6haw
zG<5Hr3R+9s*V;Z}AAEegh>uVCRPexKR2#{rJ<xW6{oz*8DPTzD6(yIv*ahX~OI46u
zOab>xS2VcSWG4(<a`m?R`V=rIE%L(rwmV<}{~@zF*r$NHr|Mfbu-FFgtaP|OtXT?3
z?F8Lx4%`fffs+=&1yaDyW0wXcZrucT91m@~|8X){mA{8ed!r4oX!ZMz?(a<oK_$Nh
z9NfGXR;)3%ua$E$s1OXD?(bg>tC~39s?t6gyemCy;R2`S@a=LB|Mta`LGaRY>BW~U
zfzvvM=kt4>1lmPK7v#S!fD;M>|1tZMK(*NT&@WE&V1JLnYer2;0tbf&zge_z25h>@
z-+p)NB)}IcGja6>XLxnC=ew!}lE9$(@v{!MmO}pSIQh%FiQu;7mkVpVjD;T_2fIDo
zln7o1mdV%R?@=&B*6DcautX4=0D9drv4u0wO$!-RJrPtpJYeIAYJ*|(Mdl;I-z9)~
z<44}AwAu#N{LA9Z<|7Hf{_H7-3F8LCZgDnUr%g!!HM<owS^3Bo9(TGUoz^q~lzW-B
z(7<Ih+&=T#+D-BCV8gLP7WpH_!3Kvr6n_{H54N0H<~eHbWZ2K4aUZi8@nA-W*)GA)
zr@{EA;n&w$!~=^pGi(hj&w(2qAAT+t9|w{bO^EK^dI4-!=u!X8$K${ipJy#b%vb_F
zPx@?aFgXr9H*~pDsmBWV{#wmBx9Y`##$}Bvl^!O8{bSb#Pl)^unqDd}<IL#waGvwq
zy~THa2ci<|*Ho5of@8N#s^CB9J1ElbuGgw5TVQIRDjvN`d<S>S+un4Y_cxT?+HoS{
zdMxn16dlRC?1a-!6|`}l6ANa|92s}^`EF>kbfZr~qgc?d_L+h8UG~AQ-AlEd9Ptg*
z8+hWAPlp3=@hAK5FFd}1Z|ypHOS>M1hy0Bst-5{#g`NzZdD-0;uKAyD;EKerAfS4^
zTE+|fp~>1hweER;1>0PTrbG`v28W53TNN1i70kUmV6EfO6Y#dJZ-RNgub|ef(OVzS
zItAY}x7|3<?+YmLVDHw07f!=|pFd2QKI99SBQI$=rsg>)zPjn=ru<*P?9QR@wjRF#
z-&Vi$?97oE(3)FeGyQ}dzA0Qn^nO4L2yJdXXIYI)a6+L<J_A!egPwz4w0Agk8D0~|
zZ!OvNGw4xzQMZ&XS7E#4J>KOuM1!}pd!)2Zxds<(?i0Bv^b=_5>^-OEjX+p_c|^>l
z8lS-Nd7nHk_}+w5fV+9p_>bUl=?g1L`QL`OtJZ7Z#s33%HFD0~R(J2h=YKV_nD_lX
zuuSdzq|TxsxIdz7v3nNpL0T_+rw)Gip`k;4gR!&T0e;_s@YD|v;Hxdk>5I-sfrdp-
z&;Hi@5q!VMyy%GZNHEROpn!ZsFbodzIMc3OB=AhTI6!3h1oo{n+H=anw_u>7#l8Ux
zL*TxyoqalAivYW-4!C7o<{9+3S@g4Mp$HK2EPcqibI;-ZDVNUn=p7D@HBUZJXnH8T
z-gbYVUYp*4*CSSJdE4n9=-zTj$Myfb2J?L)>cnmcgV)1jKW?e}8o1Y-ZhwF6E7&aE
z<3aSySK!olTf_UCU&BWQcNtf^8wRevEa20`=MB8!DxUYaMi_XovEY7-+u_jbiqFC^
z`TqeM7T-BqGCl$(r48;7dNUMsEo2rPQ$G^=9eTHP-S!t?PuM`i1Zfn6Ws26QH19ds
z`MT<<_VRadJ|9u6!Q^M4u4wo8e=2-{6+d>^zS!j{Xi=v|{(B2QLgNnoC#9?k0gWP@
zYZ)edg1cXtgwH<r1ayckFmI^qXP8p2LPZ0^C!pkI>-p9TVqlZMC(PV9>M>~FYe4xv
zy}rP4W!i1&9u^FC#YCO>UimA`Z@K0uKPwoN|Ng9uar{?!ajf@<4b6gqUx%>1?H+!E
zyB7Gq9Q{`?Sa@|>)t<*=VciBp*WD-+3|>Tk=ze4KclfOT_lkEwFlf1dx1`g;IQXRY
z+DBv727|ET=?mLQ;~~G;=JwR0kHO(#of}K+5}^IleY1Q|JO<`brxrW+O@tRx8&*BM
z;t3E3_Lw!Qa}umz=TYRAECjsn8ZxP4t7JHG$o1Dvt~~`SOTH}7SCRs&hr2a+Q13Yy
zc;dmxv(;1Kv5E(KbiDckbW3;j-di>ec0S>p_HxZjU~d6$*Eddw^Dme>_URl3e2;wh
z9``*RS`VnwY1^<@;7Hxh`<qYV<TL%JMo(My8jRRKKCtF1PF}Xxg@8f9Z@`0DV}=&7
z63HV!__SNyBmzX1Y+A?ptVmvNR(xT}!MC7U)9db^sv5`}Ki(BnziAZktTz5ektGK5
z8~Nsp?(ygyC}2@{(}ZvXdHTt<%gmO10PmKTZ{^#@P(JdUXKJgipFs22)cf<+8p`90
zuJnFY;WO|KSY^5XA49p-$A0tc&5Z%k1N~bbsFhD%F0$bq2j?$fqur9NE63)Oo7;RS
z=Qrsq*b=#~_FsGR$$Ry#Z)Gy+8+bf0>TuO(`Q&AXw{eqAjs<JRZ%W-~oL_FRvR%)X
z)4qef72iGHD#<UOy81!?3bW%tld%O|jr!)7+r-<R?Xn;qe3ADzO`ecnZW!NT!i%K|
zU{{MHP1eoNFTeD~w3cXPBDni>(77*b^2?8%IwN_%DhYhBUR}M5XMVZRyp?fPG>S2+
zjzVqYVZRzE9~%#A8vW(od$d1Ij7hb`@<S+}<TQ4o^4e*{)bYxY_}7RZhWN3R_cHS*
z;tC)i4t@dgu#@*0D&~ol*Cq<VL<l8u#%uE{yg$lx0|SGMtj3BNaOH)eD6Doy2AmsR
zb6m*`ZxQkreOhtdn&V{TmQIDpzi9YlKHu0(bMc7VEup67d_=FksQT3Wl>XjTmn+0b
z5fT?9a|-e@gy$FG@5^0yCgo0=BTN->Mu=8Sum~ehMEOtxp1xA!2lA3~oVNQ)iB-O-
zgcx8=<GCh<xhCltXCkbxQJ8ODDtw{!wFoE4xG&<0DBf?(+1iep<RqP7?`mtCkzPvi
zt1K(2&y2~wj5|a_G9o{IgzRP*evrqd5^5Li@O|<^+$<{EMc0g^sbgnmc+g3Ce!pE*
zmQWp6>55<E*%EzOiK8UCP0FIfp&XZuzkx55^R&wUxquFda;`ZS(9uXaSB~bv(X4pJ
zud<dDDyY1d@vE3K*LPYGUFbJg(1L2MQLdmxl2I?H(bwSyjTMb?ASXGXZCr2C6qU|2
zCkx5GonstO9|xh{BI%d@dB5Gg4kAjm_euBLInF)DX#drR$@rR#3T&GFD%a1nAJjGZ
znf9wIqF;{9e|0zarDdO7$@EgLwDaSq9}CG8pCJ^&<6_TD@<FDP!llhFiCd#$6Ft;P
zB3PlajRUk_g|wdm*sno8Ji&fU#*)4%V;vlo7sPJyB*Slp>UxU4EW+y$$oolXsPHu@
zua*8KrT<x2U%jxtA}=*5mVVen1MIU?_}N>LlK(m+`s&!`*Z1>kYWB(19=((+A)LQO
zJ0u}CSWvXt3GenOBayqjDqne(ZrV;rFeAcvRgRE($gINo6MQ<2m(MDEx{!Vt<Mp~k
zKWO?b3H?})2p1H7*(_b}<Vn4se20a)%>ki#Ez^#u+{XP-sIA?VDQ3u(e?P{*D9V2G
z2nU4ZC}ALCfaze(j9}{|*g6TePJ*qIVCy8<ItjK;f~}Ka>m=AZ3ARpxt&_;9brQ56
zMU00pW_UbIMl-sPZ}-lmKUhXc(y5Tlzmxf%r~#UbNBq7a{Wa$!ehN3LJ~cn3zjxK=
z8_Y6@$v%Ah;Y`Y%G)I`KyZClL%r`F;zR>f6S&eVk=HK%gf{H`*P^EAzlIuAjWDhMU
zzh2g_YtXjd>rdB_Sk*<9Ml1n}(R!xqc(kaYnnjOLC!^|#ev0UI!Z7+y?u@J3L&Pa1
zc0>91=Qu|jy_75M{P<lwnq4OOpqu?pq4^96T0X{EjX-ju&2EN_MEOX+6vfzYUf}?3
zH%U7y#qvY=pg*`ZN37TSjav((_)qGAghXG(qe%)2l^jsF6G)=}Jsxs`UdojaE++=>
zgK2}sMw^|m8;dd$xx<gKAdhf>wi6Oew<ne(!iMa^t<n3pW*u%#^F9J(T{r1>lTf)f
zU-_r){-5{LB~IuIGL;>;>G5Ay`&p1bGv5!Kl$9eYh0w1Vz9C<XE(m|3aNO7(rHgUC
zm<g)+{N|!~kQ)j~`{F%LHD3^p%h7hn=%x<;3iZ!vI<&RBma`C%X!NJU1hb-L3kQVc
zh}IjSk2BmF!>uvg8pEwI+#18JG29x%tufph!>uvg8pEyW8g4C&8rSPSZcUm=f3Qqw
zXqKbtG;Ymnj^^SJzl?0Q=6u9YVKGy4zT9$GeZB#C#Ig^!=A22nljaDEbr-khiuvZH
z!k6s!r)u->*(^ZCA@i;(g=3Li&jET~MT{wqaZ@E-e{gH2^-#%)B|tG+&vYHP=Cx3j
z|8$H{XE~p~qZ%iA?YBrZJ-5cx?IGe6nHT!=`7p6w%9VD0{4Q?IER%fD&3>oQd<H#V
zsLgK5mLUxy{nFg3PmHc#bAYy+q@6W$SIw_KxHSpZYyHNpxl#OAtqsL~Y&aogIejN@
z$JOlwlIVYrTWg}1a;1$Q1+ch~@W;M}ijIadO?Y2J3(Sw|DgvkL{Zk<vUjRNI+Xhpk
z&2E~lMHz`6sxW~?ay<uVyD7mG+hI8(49G4VAAPQS?#J;ZY(PBCSsb6KjMVI}Lo=N0
z?t?DKc+C&TgM3lC80QPXRFQnH`_K(fZxoXDMY;~fNj^^;ci)L$q~n3mc!j1zTf1wr
z8xd-(<MguPWe*30<cQWAp|3L>AH(r693R8+F&rPm@i819!|^d3AH(r693R8+=^BpD
zn;O^aK8|l_CjH}D94JBqw`@+QaeR(@HGh!!Wo-6n&iBJlu}gEl+;Ufaz5#i}vJb~+
zmr1#k<_L>*7sodm^UX_zFZBFPR<DcJ=HC+^M8zTV&?<#vkzCIKdR|40DGtE?aKZrn
z!SUs*kIMYR5@bYaJ=1j@pW`9I`Co^GI+FvxzE4t#Ub`M9DY-jTw}*&RWFF_waePJe
zQm(Y~<9Bg<e`S&ny4mj(n$Mu;3$@uz8-Junq+go5wTe+kbAYy+q@9hy^UWl`F5&nz
z<BKX{y(Y*cM1g+T&7bF{TKW;r|2ic4>UH$j_w#D1x}880of#aTX1S`WmvSY9Ta%td
z5@Le|MVp-nJB2b5xy!5al~?Je?Sup~!W7F9;Y)Vm)<TY<FC?5mVeZGR`JF|S<}7Y4
zh^k*w&m=U%7yKi1LB=PJIPU%!rHgSsDW;0#liWc!tnQ(Zv@d3-QJmzn!f`$jzevXe
zq5e5dhqiV%5N~iu#>qNPHY;8hb3jOrXuT2oEW@oa+#18JG29x%tufph!>uvg8pEwI
z+#18JG2EK2;nwV^alP*2*37U!SlfO9GNGYaj;7PNHSulD#UXwflUthe5kEx$RiB!l
z(%-x4^9^Ph#AF|CtzIVOPMRZ3)m_}01oO>Hg)j8H;D3!<3kgESA@i;(g=3Li&jET~
zMT{wyVt+Vcfd1gtg0P+uOMob#^-R}sYvKol^S=%Wb%9hp(Yc9UOCHj9a%WuK9wJVW
zd7MATt%c~NTxsXW@8Z_NGRX(s>~{*yXVCM7+UCVfLXiHDerfL1CPtMI2WY!V+F5Tr
z-%Os<AKY3L)@%L7t#MBX=YJg%eYJV|>-%{%RozY?iT?MvwP?MRD<RyPStyba8!RZ=
z?1bk_l#$3?UX`!BN;hpMB$yFFFH|{_UAQ&bGxQxl{Jq?dTXTGkD9%~jnmcYMRBQQX
z40-7|E+Gk}k@1N*8O8Zj6pB&6nPI9(zNl~LhH)GUN&6D^4~mm~#$hO)5RF1z{}t+=
z({yNScimzTrDU9}>vXg7Wi1DU<cQWAp|7I(3DJ*pkUiU>oJ6k@N}{+NrwhF&ucr(3
zqXHrkXENK5GJA*elXfnz^EI#0FVv4ReXGietnNplcEE@e<Ifn6AT5CEcS(%L7w*0&
zAIa^9<3au?#RBIGa#zilfa9hsP&(-^nvPV>$K$x&Xp~O$!n2fWJ~@tul}71kRL5~v
zEmiZm;<#sPluqoX-y2j=l*(}-IGzxJ3CLtn0`L4zKc7|oQPy5(vLdu(3rRy6NPT-2
zMOOjEQ7FdGTTuBpBXr*mrxSVP;#J>ghU1or_(eJ%F#I~huQU8Q!>=>^I>WCs{5r$0
zGyFQkuQU8Q!>{KQeqG3~Nsp_v`NxCsyaO5M>-aE?6)%f9K>O3g`2CL5xLo(wL2&7)
z=)&tD?B&<3IZl?v9L$>6(d!`0@VuYU+lgPrG({)}N+)_Pr|MJlQ~G;X9lv6hK}`0o
zgUFXjxs&DyQ+0P8gfZrumkMA0KR#QjuUWkgLdb5?eoLWt(GB}GNujeI53ejC#}WTL
z@t4#7*L1}n(QQ%|9S-HVY&_mZhI5hHnk)b30y-oLIwZM(jz$VP8jYSbeWI--;}=~k
zYe~W>>#I@5uQ|u*JFOYinWNuaK?|y>M!A9(Bzp>pV}FbTi~~X(&^E4@6hiZ`q#q&l
z9+sXvIiNlcLcK*K6&&~H>)*U_zg@NWN%z}zeO;w%VZw=D*FT|7Tu^g9qF?SsH0S&E
zSM)D6`{c?`FXc)*KYn*zWk4qRAk#_V(Po$U5=c)>bh7lUP}$7^+OI;|&k&3&Bj4Ah
z^^!`Da()E@Tq(H*1_l`dY^?CmC@%~}a;&e~{Lb!#G@S^Mzxn4r{}L0VmqLXj`s!X%
zb3UQxb8EVK`;R0#GwUTa%heUVlq(^ezfBn=H#S%hwAqQkvM4{1yWA>2jILjEfVLA7
zk`2odiDos<pZMEpynI&S)3wdV5dS!7Wl3m&3D$Q1&&PG170|bY*Gp!{d<?Qf`F^ag
z3ZUwj)H4aq@P%0uR9eCxD;$?LLn&gM&jwRP@|o60R~8LXNZJ>ZawtynS>U*t8A%}}
z(9KY2UW=wfTf6I58_`R~zq&~+D^+%JKuC^gy%G98TPMNRNw9SiY@GyKC&AW9uyqn_
zodjDa!PZHzbrNiyL{6=fp#4Z<JaDn3#`U_7Zx76*KUgL-G|SO+8sBbQPqSZz_+=96
zXwFCc6kb$)<|wh*_;L052ILXTK79M_Ov;@!8>tq4wKjihA^a$jA`%jRs<!)o?oSQE
zeDhK%JdKmdYJ9sk|DIfeI6&rIRSL%<xt;?;_Rt3V!wCbjeqBR>jQMi(!zlvPaua?0
zb3K#dR}xErD4_LB*YRk^=7bZZ22?{rogY<Cs=0oy8!>IH`F-?r`Wjcahgl_8u6ij~
z+WGOjcr>?6@<BKIokH^&^n9T<yP41&=@aRf>Z`u6a{n0zXuC=JPcl5;Ox~(LxHWgI
z*ZPfHi=y~X>VbqrU-=e!!yXd4wT0?-BCBv~EA>*Ygm7ynZIIk#7Pscsme~nyvX`6h
z({@6Fi4DN}sS?_=3%5q^-<tbzYgX-%#^o$-O^VwIQonR)h7-aDqtc>Zjx!#DuFQs^
zP>l1LV5&&Ipnm8orauaaT?+&#PV&XzcvvrzLYG9h;zIp%nhtIKb4M#gFBvE67L}|t
z+06kVIimGOs2mKp#&By4x5jX547bK`YYex>aBB>=#&By4x5jX5x`tchdaI7>bsx7T
z%cMV8CNwn5(R3QOChx1tf4XqQFXPuob3Wpyu%YUceqU~htIszek68BM);45P?xZ=w
zV%^29d1AhKsqlry)BM-CHP?YipUAwcO5s=}*K>fLR}o{1O$Mm8v-*Quvk;+@6H9;y
zqV-JIacgoLRsPe0NT_q9>WR)r^g3ve=J(Oh>1$lw9wJT=>i--M$?K(DY3Il9;?^uP
z$p_u+cM8pC(DQ}b=EeMmBZ(sYQhn7IR_;IJ0Btu(J8OyOo5@@C2e)R0^;*AiYhD!p
zNj;E|=xc&4u?XZ&-VW961d`}~zu%g*Udoj=eiXpsLPA5lkD>B7Xl0u4zJ@lKAJtU^
zPS^XVLO8w<d_J}frbe6Hv>e0irb@-<u?TH9C76CUEJuVL*@fey&vnoJI6l*HNaJ!A
z$7f-$+1?_rCC7O$K>0}SARLcbgi^#<Ktjf<=JT9|;sJ9|NbJ1KK{a0hj(bl-=|tXi
zH53}J&~#{PcO|ZfUg9^<ReD*uvWEjgazyKm(AOD`kKy<jj*sE^7><wO_!y3l;rJMi
zkKy<jj*sE^bPdNBNR8`tAICR4lm2m;(9kSL(`g)^`%KM#72=n1ouN4&@l%MWYtEOW
z#Ao5h)#n?KM=bkrd~TVPJ82eDnfX`gE{<;z=9`xaG4%Y+e~shg=OH~I^Ux}VW073X
z0eW6Vj42Mm{&2zo{lW27H$W0iECHf`)-zql@wvNcwqJxg$z09(h+ccn*PM^kv%ZJw
z_7HK3Q2*yRzIu8oSK9gUyEwikndF0R_B(~<GwAt3ZS!KTOOQm7eyP6e4J-feae%g)
zq@6WcOj;WBK3&4`X~q||#(GVVNvQn#VK;xCo9amMpRQG+ukxj;^+e~(YjJfufh0OJ
zI6lpC)m|^<N(i?mTY==p1`CKbJHf4Db|SCo(*LKIwi6P3(*ny8$(vocwFGx0f2R0*
zxgWO{BtsPEEN(4^s$WviBt%|1j<Y*}N{fCKqf#9AK7>-l_yHNFisZA}g|4J~P)O{W
z<!Tfs`J^~*x1FSr66j_q)IX=`(AMro{f+1)<7C~WmX#{II3OfPwB87PpW)URZjIsA
z7;cT>));P$;no;#jp5c9ZjIsA7;a71aBFVVxL)^hYnIp_tZlylnb6QIN7HHC8oyJs
zUxoN(Bwm{H5kEx;Ri8OZY&L#eeZB#C#Ig^!)+3X0C(TBxg<q|^xHT)xH!qdK)ANG=
zHEu0oFVZJ6@2XNb7RmJ-pyyS@m|_|BhZ6?q4{j|6>lv{Ghyq&AbRD<G@6Vi_i=yg@
z)E5bfUR!x<ejmxN@1eRqL~=l=|8v}0f?mp%c7FUWZjCcUyws*&y4mj(n$Mu;3$@LQ
zNqi6&NWWCy^@f%I_c%b?P14Q=;`wIsHvPe^8DqWHZ`_*rFlp&Y4XB1hUtN!A&ZqiR
z?vJV42_(_~9=BFZFXc)Iw`O@1$&C#b32k;F;25(LxnFyPDWvU$1ZR!$S8ac?3%BOy
zi@p<tznA-QYwo8I#W{;x^TzFjkYpe)9mkmlq0$m*;BnmU0ZI|$7g%DdNIv5m=*sLC
z3W;6gPM|o+XNKdZS4fIbiF6}fsDDn=p{?EZx{Byk9i`|7jjR;e&jBGhqV-0oBtrII
zb{5H1X6-+B9_9a`pF$sG>A8{vg7*KkYI~g3_8-l!8F6CF)u<*08u`oPyN(M{K2kav
zj{7Y}=@vMj-%!<jAvm5e9Ho<fPC`@FeBx#(&Rd|6=ta~+^hNSruSj$ye}_Ue8sRue
zG1d3k;JB*^N+<0-=>?GY`C<DRMDblCQxrE#aHxdz^Ld~j1rhIr`%x<{UANBa`6?v0
z`4qc(qrVv3-U`(d=|{<F`FVxq=M4(Q_ytiHP#pa#`d=^5l?kSn$fG5l?>5RO#py!j
z(Uo)|`)5kgLAOtYPcn+mtj9wrOQ_&KW*E^kx;});C&Bea_?^}BQw+DxaO(`W&T#7t
zx6W|u47bj3>kPNfaO(`W&T#8Fg<IF=@3y>$#;0W5OvXhl%o`j~pT9?vfTV)s3>AJX
z<+X_-8Lb!)XD~+ioHWx%=vm?Q1xo)9$H{b(!Swn9uls}(zb-EE`^fP8+plwRKm1Y_
zcQwCH=y`ohSC=c~am&8-1r;+XchVeTrS7gT5M#c1sqp3h<DRp_U#l)3bh~4k(e9KG
z1FSnd*Q7AlB;61dQ@H(+V|vj02vIaWx`yKx^1FIb{%E23MT_^ST_N*}gzj|3pDHKF
z9i0s2xK#WNeBmOsHCO)61$0OhbVzam9gP%pG#Wi=`b1kv#xJ^7){=zdl=m`z%{flr
zY0aR{9R21BT2M_j$`!OA*;7ay`(qqn91!Aww(-BkOVn>5^Ql7bVd=S(17y5Q#<6by
zXda)EzN8<<nF?`*$dfEg^ZSIJ*SGZS8_xvl>3F6<#`q2WaEgE)-~IF}p*ws=l?#vS
zf-=bmnNEt1wszPu0_h3SL!Bgo6)M{}K;xB2`x%CDXyn7Xv@THT;T6``5Uj7-{J`#n
zG@S^MJ^Z<!y+ycY`$hCs_EvK~ognT%EQGo}L=v5ub%C1Y>V;m)mA3ZZ`vYtL|3mTA
zxlhsUzXX>s4CB~IfwLO#Pw=xeZa(+p)V)3<jmue_x+nJAllrAYGn`<DOG~Vb#2Bqx
zuquQ?G0tZ}<r9n0eclj-#9qaGL~kPRGyR0FBq=0?E{Sf%h2~3WI<&RB0jY>yGzQTv
z8d+(wp94a2MC*-EK3WfjO3c;;uyp}!T>x7bz}5w@bpdQ$09zNp)&;P20c>3WTNjX1
z>jG$hlNb+Ntf=vS*00~r{JbEVIY1>7o@WTlq(4|DG&IZ6bTa===6g&y<fjoT6!G71
z>6-mC#Gm0!)u+&V^*vpEz5#j8vJY<_o=Lfr<_K$b7jGVg`R1j<7aE_P)%$~J^8*Iv
zLwZN_P^EAzlIuAjWDo7I-<&WYtM?00Z?{m(O<MOq-*?<S_*b57TF-QC9R#t5#Lk%%
z&}<KhoeQAqiPRSfiC&A1G{2AJ*Y_~{<VvQOa;2Rgzl%rn%p@Olv)`%B9)=V~dP4M2
z-}Ql&|MxgR+e6ZRO7VO%d7J*=(Y&y}>Ng%Okm5h72NDu}6&KALc8<`kEmUv+kwpJ{
zJlY<;lq;ciq%n9uOdHHlZFa)0B(oFRWG^@0r|pCU6YGiPh|r#0xHWno*W8a=laxUk
zm$SGvs}h>+3-Z2koV)?bN5->JIBwhsrHiqE#FSRe=T{5GgKDFY*fsA`s`-L&TwWEW
z6M56sPzXmv)1j^1wXBZlCF5URrI(c}dpIB@N3`AueVyUa7#@w`(HI_$;n5f#jp5N4
z9*yDA7#@w`(HI_0*YIdz)VN;v@o3(e^asm?hGsdMPUhdqe2-U6&3+Z)myy-boR9b^
zEUIbFm!rgI;m6hI8<0mV`|xOoGbwk{ETl5?uhLySnjhwymkKd7PUgSHquJC&dP4M2
zrEn~g>p4KrtB5hhG1wnY7@$9RG^@wC0gvWYPqY0ZcFwYn=6po2{pxGZN9tMMLv?$I
zI7Q}#{v79MqnC1}ogcr8N3+W$A9S<dDKwu!&lhT&7n7MIeIor*eb*aS{@>#OZ8u3f
zYli2W$=mb?x8{iTTEB5?ZWRAXJ&=&-Ye?g~VGjx2+Cp_Zfh79h<JKnVrCe#_M*+-G
z5*jM@9aZ)*N-2JnX~OFWrI;VpRjPy@2kP2-2_YO`06r(%22-QWZkn}Zc2lL|^H_wo
zn-WaF9hM`)j_ktm(dWD8ejHyy8>DeLi{mqGq1oOduce~Z8-VhW@tPlw2ic%>F&2;j
zOclxJ+8f>X?1MsL=cTPw-{*<r?%hy2kvClpg~ls19opJmlOBj(;y2J$dRe)$hXX=#
zMC*;v*BOqF;rJMikKy<jj*sE^7><wO_!y3l;rJMikKy=q4aetAjq7zE$G0Jq{&AVm
z&@4yOX&j$pFU@`x;+L`MsX5;dKSg)V`Er!_Ed03od;{``Wgm{uGm~;B%|a?O|0><Z
z@omF=^HL#(p1=98aeU%_NKeQ-v`XPvB-e9*o>viLiUY7eoG?ItaD1(wAm9EEOOQbV
zt!KJ6PfF||v2%{rn(Y^%&ZMvAd_=EZ`)kfe>RI1Ib$f_7Mdo|{9LLvDFXc)*KYka-
z*CUgB(9M3Q(0m3xU#M+f%w{msC!(+Vt~adwzsCXEZjyF32G2K>x9Jj&Pcy#A8tXMd
zCLs#+!*2dO57lyzW;;jp)oY06e5C&JN~mrpkVIz&$ER7Y2J59<3E|eHwn%PluwZDj
z6Jd7DPUICG`W2nDosi&6rdW>1Guee(3mJ;!FJU+eb3bm)Z#1GbXK`ymRQ-~ACL!|D
zaU4Gfm6nW89C6%z9!eME2c(!Pl20-X-M4Z@A+c*_BURsLh2#8Wluo{{tD#W;oTfut
zyBjzK(M!h3x=JrASN3o~NRDW|5&AmAtufph!>uvg8pEwI+#18JG29x%tufph!>uvg
zny%s2?5J_Q?&H?%us>MaegQI}p;?ZmlX-eF-y?R>>{lUv853vC`G}t)fU3_NB{myB
zu0G#@JYw00TN|B8xszri)xxjVUEG=@=9`yF;put7{~EUzG6U%o(L<HOu}H4x06niF
z#uQ7jKb$Z?e{gHYA-Ms!CZ3f!I~PdR6R9r}61|qp*8Dz_U*AJ@dx&}jWWMLmacd@e
zDOcM0@w>P+(@gR~H~XDJ^BMGfp|*K3lle%WNWWCy^@f%I_c%b?P14SK<N0RtHvPe^
ziLqYmH*Sq{(`@I6zS=C%oR8FBUJ2Fh1d`}~k6Wv$mvSY9TQggN<i-XIhBiClxs=(7
zyrM(DqLa205}YY$v1<F1UAQ&bLi8O!{Jq?dTXS52D9%~jnmcYMekBviG~}h@xP;xP
zv}AlD-h<-&J`{@a3(PQ8Bwy4vbl-Rf3W;3{b4PKK&v+S%CwQO`zfNZZq5e5dhqiXt
zZ8M^m_=9zpU{<hf;ee1F(Rw5FaWp?+#EEf+u&RbaSu@gkd^adR%187s0mn^^QMv`r
zXV+LYpJ#6r59otJVkZ*9RrNWh8oF;%6NmVH{Cd@Va){z#7g0#YCuSDtIZ{58Aaq~y
z0EI+fg4Uoo$(MlRrt45T{;bwOn?L>regpXq3CUNq-T!lcjw9xnu$}O@e%(5&{qaby
z@+o#JMt^=Q38lYgC9xBxRQ);*M)l}E6oq2^xEob(rp3{HixMa#b|MDXBjL9s0+quU
zKmY5}=uRr#uOjrDQRUD2d0b`7!f`~ON^ppyG@qxAa<q4{NBdDL>zeSn%b33Qou}AO
zCY%%+XUK5<WOke(Uxn(A)Gv8I3-bmCXni2|!-8tBvum6|+lPNU&iL6r{aWAX_Daa!
z)9s{?zL{P?`tw8I=1zRv-BiRo%(vc>>Kh3cOQT^KGZF7w^XUs&r%3n$cwf8SOvJYz
z|6c6+>@93rq|)#y9Yy@A&P(3EocI=ATG44xkv<|`TCmy3XoI(K<-FFRJBEt*J5J#P
z?re&H-<JAWn>&bjzkL%rc-4u3fAfZ;Z%9S_m<tmJF1rv8n|F}B95+M6_pFz+z`Jia
zw753?OvU*k{z`+riz-CEf!9TM3q>py@y~AhB&N=K16TAb`t8Ol5g#K7oHVw?8~C(c
zq}Tm*BHlvU$<<~5Yv}uZbw>ja5%1~S@m%w^ui>1<WeSblBI1uu@h)*C_!VqBWqPUj
zZ6aQ@RVJ%9@fF;V!f@*hx6W|u47bj3>kPNfaO(`W&T#7tx6W|uIfYx-T|b)JWi)ci
zSwC6=o(B;c?Grzm1ecchV@v~3=`GHnP>f#?^bW;IK4X6rH#>$x;%D)u@>%Ud_X94W
zP^dh*lCE_gWiO(G_$zg#+^ih`4F}Y5AxIL4UKlF;pUP_!MKW56BF-q7@cDakZzNH|
z>sFM0caD?&Hw6CXSLk&s3BH>3oU{{xcs^I?b40(Sc>Y%?Uv5cPmn*Xj(z0*e%9TvY
zois;SY(nM{C_f}AOuFKa$OuW(;ZTmt#*_Iw7pbke@_#O%L!zKVk_+f)q@bhG=t<Kj
z+DbBh(Y3OcBpj!_m+@=Par#bc26g7>H&@VtYN}DLpasdELgLsT;{fA;5C^o4{|n&<
zNxz(gq<5<A{-5J$0&%~}OQq=le_d;K;ApgsTm6nB$_Y+GrEn~g>p4K<A;g$sOYFZP
z4A7r-t&RtfM5~ShXg$;Q^B59NkQRy5fNDsn^P=iWHP=u4mGLRf@1vj7*SPw48F7lV
z^MBrcOZ8H&wDaSqUkY)G&!~D~zm#hx`5@Ct5z=N4<!6zenCN8bS)sC<1GGIP{URGY
z??Jw=OY2pY9%5mAb;J6q&5!I(NYjZB`Hz2|v*1VZpDrZPSL5?}!_Lw43&qvje<abF
zS+A;Dt`_N~ToD{Pz3$`R;&2K7X}RTwq4fmhJ{=OdH5slsAN^czk888j#ygpv&h0Xw
zi9Xs+OGsu4Pmo<WetJDvR<$2AU1CgE2u0U_wqA$AeVDvJdMA0Q+0L^4?Af{iwl09J
z3t;O4*t!6=E`Y5IVCw?dx&XE=fUOJ2sdWL`{0RXMnLmLHq4Q!$kGqLK!RmqPI5(^P
z3EJ=~1Q$x?FG)yBrS1No&v%-<N8iZV^(x*J|7qDY!wHseP-zK&94jGzsuZV-aX#X|
zAo(n?Kf?z5Gl)OdirW8&3r6o3KOvl_OQKtGZS9?R2%-bc`RZ2LS!uGL19W>Q#-HRV
z`ObPjV*#|s4^}`C-r<Sp<L|GNanIyDlq~z1`UhX2*U50tHMcLMeTsRiIUmt$PpUru
zpon_$i26JyvJBa`PG)5$<xZL-{L|L2(1vr>cK^@&6*9~>FBQJf_`Iy%PflpxBijCJ
zvzuZ6AiX2F5EUX=B-e95$Zon)^I`gfmz5m+RZo!CM_tFun!U=LJq)7ii56I**A}le
zzmI-SU*qa_6ZIT^_Qy;!`fDY`0Qb~*u1R68NjiqB2;<Hy^ir<0^W%5%vV11_pqu?u
zA-hSBujzJ?=#dPMzlnb64{p*D*PDLfCWZEGbjSN3+u(|z_gB@{@3D*em6ag1%EH{v
zf!}Etk;aHKSR!m6680mB(xzX!*)G!KcH?52@tb7ao=`+{J~D2%ip(qhk6-UwOYaXP
zWQXYehJ^Hx;n5f#jp5N49*yDA7#@w`(HI_$;n5f#jp5N49xbQvXzJtX?4Gaecqnt@
zZ8l%2977r_PA*dVzzs!~m~URn<R5?Z{76>MS8DUeIi{f6B!1Q(1!EU^f&;X_T}-G>
zR`na{&pc!l)-&P(Lj=)!rt9;NeyOVbC)GkiLY)m&PvkxciC)XoG{2AJ*Y{A}9wJWt
z`+3M1y_75M{P^8@$b?MtK{x&F+U%jLA@ZveKfAu`11taUae%gmr2RB8P@QklrTOOq
z3KE6q$GF3&<yYqV`KKTD@aMigjuiiYQi86CzRL5d))SpCuf^5de<abF*|$ftTouqu
zxe~%P$nbg5HdsKk*$J)?vlDqum;OJ!w4IRPn-*A(NZzbo=SB9NqR%=1xBJEr>!1zq
zt!<r{w)=mM_qMb}8ke){$YQAaCG|`~<b6|&E-^lXj5~~HpwgSoLZKKx5LFn(Nj~Wm
z6gR=o6T4<f<qH{)@=0+z`M$1(LbwZ>4qC4WpJWuBS-+m>*Z2{-K7`6Aq3}3a?=Sgv
z`DnTNx9tFtAN76%lI-kT2a}RX9G8?7@jnf<wTs%;xoErp=j|e=ll#HTo1pfIjGz9~
z%3>eOaSjNzlQt!gT#y20bvtQO4ZYvBhU$07_ejWad<@6OaC{8M$8dZM$H#Dd49CZC
zd<@6OaC{8Mms2=CZSxK8Vl<v1{y;L0VPW3j0DT@r6ty0j<^qwEpemZj!Sr*vJ+AK8
z%=~<)=5tJX=;=CddLJ!Sevo+~5*jKHt;)HIQi>nV^TSq{AJtVvO4o4*+S=iSTC5#b
zsn9$Yq5YZ?Ouq#0i$K_sUGu&%*uQClzn69MJx0`idX^cai1MMCNp17G+V21P{6+}&
z80YM~ZxmI(q@GDg=6x+1qxON!=X!QU@qq3q6ypb6C93(na9nPV(uw`%YpdpS#Buk!
zD4l#?S3{wBUz!eW?XF3EM6a=~63WVztsJ1+T`~TeH>C$zj|-stM@RyRe*gYHqt2P^
zC(5#wpMUytdcP9KhML<U(vI0Q(46mwUf0u{?@!CEYtN~VLy!f@K7561Cgo0=BdpS0
ze8p_cH!l^w(Bqh_UbiP?H|c(YHoGb5jP!}&FjT59&qRdmW*F{g5=LY-4ns(O>He?~
z9$V5E$*m18;M7r$_D=T5@1V35#QxIz&#AYow7ve*eJH&!&gzGy$|xD#4(R&2LgQwd
z?JJ?q507Jre9$3jPfVL@ejojuzQ(gpuDtb9uC(*xch?pAWs(oN*)JAq$87L<tE3$h
zY8NwmE(pV|G29x%tufph!>uvg8pEwI+#18JG29x%tufqMPT|(Hj|2bI-%>G?`^>oO
zluLokkC1r}*Djj<BgEe#?xZ;%@wd2JYR;!q%=j0D$v(NNu9tG9ogcsJZ>g6_KIo>u
zMccdw-xKKx(ZhcsitJH6<N!VILE6s{jKd)$>CZk+O|ZTaBqpLjKjt5P?&Bosq1k>B
zeRc1pIUlLNyb`Ln|45?${XR}D^-`|1wf|oISo@z>bmLca((S(lXA0?~+OA|5PS>+H
z`c44;Uhc=~%J50|Ig8VkVcfv4WI~x{hHsb+M)gB*1VQO&m&6zmN*CjN9F<R=fbxYU
zp^(_CC~MXCne<0-$v_m6*Xe2~G%rchp{?BwutD^aajmY>%gU8K9H84>G5(Adr3bo?
z)8%o$D0iGI=ogRLE1_B^?U+}b=6po2y{Y;Xdau5x3$=^19T3V#+ZAGOXxy5%cGdNQ
zs=d*7zOwTF9tY@lmH0EfhpP53bjhEgdA)1^#(QaNU)_I?QQ4gP-;BdG^W&Ue@*_}+
z<S%N9svOeQz;J5}x5jX547bK`YYex>aBB>=#&By4x5jX547Zk3xHWD57UN0Gk4=x^
z@={!VoR2)8#DDku`#u!G=Ot<LcgS>-%>>s-{7!b0HIMg+KiLZV2ZV?s{^SHF&F>R>
zUf<Hy<;pCBgzVe*!6cJ%C(RL7>h8V|e__6PsqlryTV?gW58CXKmkZJpqK7JlW073X
z0rmESwo8ACw?B^MowUe^m9)I;`u+w{(=^)=LY>T6b3US93DY&_%Pj(zeR6e9FXc)*
zKYkZ)e<hQA(9QmzHoN3E8|ew@_vIG7V02`T148rDmiW9w@)lj1pVqv8@olWH`i-~u
zqWDkhfrLa~6XxU%J4fi&7OJ=Zs7LVc@%9h&Qm%fSJ6jTQPTmfMzk_QwQanQZC+#Op
z8a2Ml7)4z9kK;&xi+sm@9+F4z`6$$8f6T6-d={u!>^F~bfVMw8rautN4`D~Yu?|in
zY^QUhN7;^bQn);n+scg{HDZFTbdt0EMCY+1CZL>jJ|ee7PRaXp*MAXs1Zjj9DQZUe
z|HpG%!mnHZ*ZvE;C8$(1A5{v;BDt0W>UK4&{1=YCNV3QXAF+^@cU|{icrQg%2nhh8
z&T6sdd_=zj+%@MD`i{P(XP;b2^-`|1^W%5@7p|G)gKqjS=>CBO3wYErRsG27^?*YC
zTiH-F59l{sbv}^j9qDiU?mB^}1Vj&MC&>F)m^U~ev=6By3CV?x02|tszB8)d7LJ#-
zcCk5_6bS!m6k98O1;>xOHsO9!(Lgv))M@m{IUMg?p!xQ1HLk-wXQ~vLJCWo6PTNy!
zdHZW{cGb&?eQi1Z*0+$_$46a-gX?^Ye$|`f&vY4e(qY*Z_~Na6V_jbH{`coQZasDx
zcE4VKzC^<D#*3Pt>-72(yv{wU|D`g=U)()>?(ou=;NiMU3?>%m_~-Ldq@g`7LUFU8
zoxM`i!DiDVMrYmS@MF5@aLecE!2L+SEn#aR41Tl9#Va5k9C4T%@~Qj<I4Y@2ZMTi-
zVB><8lh++S2M@t_x1LT;2mcm-I=+?9(0Tbz!foc;_ntVF1|EdHFaFn!<8a&(Q@b~^
zG~oFuvQdRi$KalF4Z1EJn+A-l?6ui6{wS=yt3W@mc4@%7s&~F+UH##5`9Z7t<<h{L
zrQt^t8v4O*>D%J!#iW8KU0c0;+U5wfkUBlzek~P<BfRgttaliWJJH2;?zU7AeD^?N
z+4=|Jjb5{B-<g~Wx>|ocv$L%?JZMs&p?l|4&|2EQ*7gzm;ML<re0<8Mg5w^e+DJC-
zfddxUA8r+$0_Ijm_|wznOI46uOaXmMS2VcSWG8&I<mzqr^(jD-8hK%U+a0i=|BzW7
z>{CF8fcn-AEVjXiD;=&6YnB2Ibb{_R2X2O=!AXnY0x2N!*rh><TQ|WS$3xrhf1C_f
z<?rFr-e?0XTK#^b`+Ji?P|0rr2RE;U6>H4xYvr5_Dg;BP`}<eJswU32s<clA?@A9_
zxWH*Se7oGkzkTs!5WKWpdhsPo;Ixk6`TU+Ifp$^R1^I6a;H1L9f6V?Q;PoXw^o!Fx
zIM8G8no(1dK*<roZx-#F0UNIJx8L1532Z4`X5#7%&hXxB&v#V|B!Nfs<7XXiErlKK
zj+4K<n+Td)ez~B;pMH8A?DlX|B3OB?OuiO>kAfkxPRCP+C4$fd(Ce0oEu497TF9X4
ziJ;oy0UJ+L8w{H-G9MBCE&<FNKk`<k)i$u^UlwOJA4vfAXHPjy7&jPpi?iuEZAt>D
z*{z_-%15^FxYHf!w5AE5+{?6u1}>xF_L<k#Zi<fw8;%{a$R9BdHaOg&_``sBu;t7$
z&ry3P!+s8p`<Ts$2Qxa%b_sqy4aPqWzrM~Q9$2iIVQW}<4qWf}@N=>FI1spSLUi}m
z3!p`zNBuV+j{__GpS2h<V+mY)(r0Uf$#KBUz~xG%9xLG6Yc=QGsuu^M${1BDJxm4%
z$F2>Y5cwShUMeu-%;@!SuJhZy#dm%O`HHV!Q(3+VI&PU%!GF+qkOJ;{t(vk0Cibb~
z(W}IF;8f1`rt7@Fp~tNqCnB!Lf_;~wBYBsdaO$aoHqLWmf$z+bac7_JhGmy-^hszG
z3tHAbGqAqPK4{gwRNKiB-@wNKCqDUfH~^P@vj6_V;~O~Nu9LU4>tT4v-$>G`>o-v7
z$<UdX-F@Mj|M>>4Nc;){s@JP!yucrttgTb)p7&R<&828c^zdVFm}t3Gfq`Ga+`9wT
zIu1PnZ`=ANnCJTnYRwwG_3^Az@J(~ujRXC@fD#Y(ZasM6H0<~J!<6YmzJNLMl9pp?
zo`d46n{IB({{_tM9QtnS@e44k@uxF3_laB-`U#jgd(WwPBM??x9uYIC#wXx2?~}&`
z-<xnIa5qmH{}Ie9ePKl@|J(3()q3r__<sP;N6xw1>h4|m^sh!1^S-|a9a1|#sk0~u
z{vA=a*gcE)AibBpQwP8M(Ac5A!Pr^vfcd@y;i(@Uz>qD;>5I-sfmTIN&;Hi@5q!JI
zyy%GZNHEILpn!ZsFuWJ!ai(3nNU%Na;sBB16WFKDXwNAN--3RQ7W)P)41qhlBK+y!
zRR`R%E%OX|+${Rpv`_>Hd6qup+_~rQ{*+5+d-M(m$C@V}C^S73UT?d<Pp?gH!0Qn!
zw!H214|H!iq~rR3UW55Q5p`lWgu(0Ku^+e8eGS}ePPf0m_7!CIB&>OoSvq)kp<Nlj
z4FRxW*HPAT!*o!{t>KVimZ#ve`VWlKL(>4WC(NEOd&2BVyW~CI<u*ivx3hbsv`)DO
znLT0l#HVLE2s4uPp7riDtiFiZ6E;3!;}bSMv7WQ6#w9qRP$i#%DW5^lK`+`n9J&mz
ziQ~7H?D`q>D7~m#N|&pU*%Qs<lbuHs@;O}n0y>|uZM~+>X|?g=tHXyvns5CAj^2*7
z`>X#Mc;{T#mqtszfW3p-oqO=<41D0>aq_zK3vg;%H~NG7S+((IRg-Ns%69$&;=0za
zXL;=$Y*WGD^*+fL(7My+qMnZDVYgum&(*B(1w1b|Svs`n1-0>M`<b_@#($0h3s=<;
z4;v2Q^7Z##ba@s7?#k@!`{kFzq)*dIEV~i|dcHM(3Wv(o#<P<TMf3~xi~+r0ml-st
z_(gcgsL$q@1u?+5Rr5i?<1WJgNlpgFPEf>m=L=hK{UUr*xP<8afEW<k9LZOFY2T4S
zd<>Ye^z4m~i!Q;X)tVRgt`h@Z1g;<B7jX$*?&V*5@?SB)#yTbSNterN^5{2h%yhqd
zpFxSMZ(>H6T!9Uiedj*y`V6l0TXFl=v@0-X=d!Iur+o&;FP7Tv{ql;M{QAVV7@Sh!
zGq8)_xAW%FtMK-+rHyVyM+4vDMzXS%ufg;uc5d@7MuYi7-{+59c}-2;Z+zTce9M?<
za5wac^Z%T#!;XJ9sw8R=4QgHsow4$t>#)DwX})*iXaLN6R$ARJP_4f;=gz&wyL>-^
zgR2wI4&raX-LWs$d(Zy_8nnHdzuxH^aC!AYt33LB0>Qgd*6eJ2Q*FH8_jvIR=J6lF
zW5Y(Hme#)oof`Ei{7>LVaJ|T8!1>&Q*Jt|$PuublbPTDmuzrKvYW{(_ueP5Y+vX#<
zWH}3-Y<dSab+SmlTI?g}JEh*r3g_;?8<(5^T|Ddq82)~OL+6fn)%*!BlblvNE%*Qm
zaWy(~gYQA(GR+z;?DqjQnI;?T8F>$GE_O25z4`}GL*~(YQL7*||AVz*>5>lD-h=rq
zOiN$769j$brK-2{cn>bkYdoxn>3tXz-MUFJ$M?YeXbGd+qwcHuD<;LCU$Cg)d+?&J
zLG<y|`*0EW-Lc`bcOdEW<#$zkJb-oEblV<u@Ev&7e(IGqTOX+TH;$Dl;+EL`9guuE
zc5G6chw!kYQ~re&-+{@arf!d1^AJkj6kS^2QxrJasf~%_n}=%tkm4JTPEKDH1@<MJ
z>%Vf<Bj`H3MfkMgQQ&*|+9plIAHkwgZN_|T5Cs<Alvj;u6|Ckz@i!a)J@{TEm^Ur8
z@&4#w`08%U_C35JLHy^MX5L*N!_7rYRJtmS1Y4fgwY|Fcv6{bSw#9v`%OxU#=iO<m
zhL3mx4{mB&rTgo*V7=As1^>vOK<j1qr#?IS7UUN@pIA~QM9n|bpzoZwe_6i;y&vp)
z)aPvooEy3I{H5w|!4#W5@|oS9Ld#}77SH(}0aku@jC$z%RL!5WsohGS${rD*c-a+i
z#dDv*sXN>Y6mpCJ7h$WD<szTK@AK=mGinh5o{lJcqDP<SYW|<&)fXS+pNE4JE*Gkd
zul)k@QMbyt9S#SHoledi;_(7byIU?SWKKBfkhVT#O@UA~f6?=eeH}x^;owlkMuYE0
zg~Bc6eSB7aeFNMlT(Jrr@)F)XUTfIUD{sL4JLX?n1ie)AFO~B-J887z8&G#`vHcO}
z{(<E;o!#rz@(rl5AmYTf@?p@t{XZY_6@3Ft!zxsqGB8ZdA9cLj%(Fd@z6R-)0}ZP`
z41?w0opfln_%*ng^6ZL1$yd-f>EYPkHm^Z}b8Fn!_Iah|ze4NXe(qar>GA7%Ix?L?
zaFqOh)%iqnKd*-B^Uw&HkQp4MlH>U~;|zSDfq_BBl&G-+DX6?K6rIEPCe?b;b-af*
z9An@LL=Tatyec<&m2QH6wem*oi0lRme+Or)ockBT4`{<HDvsFA_#@Agkn9Sn?Y<$_
z&k|yw&|Kvo&owE`HA%<KqOd+&tVda_Q~>w1K4&#vk;a9HGs>*^IgAU8<5|65DXC{~
zE2Qs%7>Dw2>$#Mo!Zg)CO|K=UFNUJmMQUrVyrl>;pfZZEkup9~iQBfFIKs((sxr!r
zo-}=;GG|7p$dM5$o*;5V+A;dO881~l@93k-=d8xxWE481q@*Mnksm*#9%ScGIs7i5
zFkSJ7<o8Eb9lnTxsBkfW=%=#&l>c)99TKGs<^nnzDRl_x2hx^|Uu7*RltFnf<5!_e
z`c5mN3;pH_T2QSu$`!OoGU|mizdGEev7#{!<QxY`zaYpR_0M9KqwwGM&q)7Qcl+V8
zo2Uo$-}b|iza<|YGjN>EN!5NKEeU2g!S4#nNBWxyrl_7maJm@hli_hS$tQ>Cxv+~U
zB)vmNDxY5fx*v1~g>*S|Ew0VKV1xU)`qQ`7H4f`{dFXMN7?Tr2$#+(dch&uXNCJsI
z{{DVy&VEQEk-n2!q5TCMCnL9XDx~*QBjZk?LX!5$>%8W6i0HNDS<U%`zN2sH>T-oF
zL-y^b=9)>lljaC3we>3s;ZhO3BO&qcX}kaD{=M0lZ(b^Vq4hPZ{d?N%q2pzwCuF^(
zO5s=}*K<I~9&(pd+fS#l6Ek?>G{Q9_oH}-9hO3R}kHv3tO1<5p?cty9<5tlhm7G`t
zR4`i4bbTFkNEubTPP3O#XG7H!osZ~sR9VgMqo32)xX^l2dR#=yJ>gRV9ybyCvwpuO
z!;HBB@;ws{^=wS5p!%`FAJ#UHLG+GZA1k!}QHs~alJ^Vo$*eXiM^p;ZBB6Ye`nbIH
zP)Pc75?p@2tK@l<AX|>&L3o~p&@T=_&l7&rp#@F|!rviun+BujES_k7o_?RM$F<pO
zuVlnW(jV4!dRY0gmIHLV^+VtAJURJVR?i~}wgU^1JbJD}dMK1n+#TIFS&u@Z|HRKI
zi$d@5dyhgw2OY9?b!;63TL;0`L9lfYY#jt!r_J{LVEcZseLvW~A8g+bw(keqM~dwu
z#rBb6`$(~Uq;k_fQf!?zTc^#|X|r|OY@Ie+r!Blr`zLPL_=@KFP=pRs;ETVZc@{FC
z>W1T9-%+|4PrQgfs^+uAamQ$sPUaz_!d3H`M4-4N5{2|Uj;_Up=B;Qtw9U^3;CU>9
z&(JkGS$VRT1N8i?7=O-+(gWSc7Z*8-G*`7fQO8lqvX`NM_;!LvA@j^8A2hc^q#fh5
z%s&&o_NMCd4~wbukE!GIkOj#;e6dL;<xZL-ywhEL@n4v4UMhT{^);*U#X@$I?k8xo
zn=zk}eiD6EA%aD6JqLvBraSIu5(elG4*59NM`AG$1++fuIu6<5i)On?s0*X&iPRSf
ziC$ZN)%-q^U*AJ@yNNhO##?{Bj^vzP%9VD0{O);7S2D>5-Rz$V*-d&JuFY<m#Up+C
z(SOx<#bM?DEe_CjleDvP?8hN5%ldT;1<)QxxIq))6KpfHubq~2r?FnEY=gz=hu!?S
zU!rN8W;;jpm4B)_z9#jRSHmB8aB98C@9aDAt8<q5eu&RDn6dP_8RRc8;>7sAe5vMo
z;c>nwIm$=uy2(Wpmt01n1<q$ht%q@|gzkHZaekc75toy^kKTWS$c<epk}r~v-Y<`|
z1Ky3$`$;}}zctc-^INI<KIwH7x1{zb!^IQ~Nd3{z|4#iW%S6gc>WRd)jlaqMUqADs
zqLTH!+V1~(zbp!KFlYA-k>PXhgeu4k-v}G6Di?kmP&{Z83dQ6F_*`+qPkP-Kk;f29
zzUg&j<aN3l(sq>4Z;bVflsoJFw!dBvC_RYukF-}p{Y8^isJ_Auq5PzrbV#q8NqG0P
z;nhtXKPRP&X)zcG4~{SCV7{K?n|#^#^4!8e=)5hU=+-qHUujB3x09CY=OKDudvJQV
zJI60rw0Zxshu7h|!j{)hFXH(By$*jea`SaqC8@r($$XCYN_y0MiS2c@^AlUPY?ACd
zgX6=NSABW#<29JNep$st7mkmckrH(1*flsb|Cw8Drf__(t}}bY&Ag^|-eUKT`CgqK
z$MGE-zG&f4=o<WfWo6GVV>rH3V#wfC53a(7mG4?tAI0(Dli$dcEmzgfXY3Hr3(Oh9
z@f)g}EjwQCDlBc%+HChgj$h{SY*6a=EAZg0wx_Q2<M_$-Hn^;YSJcjP?AGMv{iN<3
zAMJiQ&b{9im>BLfvbZJ3FLd=P@%4XK;GX`8yDE3&_)GO|Uhj&!takom?{U$s8@J~8
z{l_boem?gy?A(3hmSz?l-~4q8$;fV(q1&mWeodQj{QKPlr%KAGpBGtt?~|gn>T!J4
zHx=LQJaP#hIzHZ}oEgU_kG~hvan>c+^~{=)M%6jKz1#M24J_5qmpnV~*xAPwIDYU8
zm)J+IFT!MQL(zb;9DlIHBcHW?7hzDRYt=WF;`r*LlU8h)dr|E?%9TM$``Z`hc>9kF
zz=R4H;f0gU%9JwV_~X`hp7HPHaBfudUjG<ye8QdMHHVy*tDRrDaaW7Bk7Lup)fcx$
z+S$lqcyfcO{XeDyiwODux|_*ijqZip-+q%0rh8v$6Bh&3&bu6KUiqlkgLKgD;3RSX
z#t^O#@6tc*MmlI_amcG*@CEp|ealV#E-Icc_4(uT=@-<_$82bH{9CwhI;c`Mx%sjS
z=V8G|@g`OFrvs7UxK}<S&%>gd=H9XXI~{~rHa|SX=)Br_nhU!(sJ?rp;yljOGlly0
zJqMc|etZ7>!gR1IcB=K3gtM?tUs2q>ndxAO+_UcZ9cR_f-<*{4w(WDrbO7_&r|peC
z10UU~(Cd+HI#{{*%h}={XW-cAB?(ssq$|$rJGWfKpMjwr4cq+x^Kjmlvn;==V>(!H
zAuwk4`qS`u&v&aww^BTRe&B!tZBN5zD+g|FY@QD8yt+QN?7IN9^K%OS<?9$wIUNKA
zuCWZ_1K`Gmw`xx;lMe1KsQP*Or&F;0(F!vv6jhuT>i1XDzAdNJ&Tk9-*x~fhq%?5(
za&-Q!aVKHT32pYq#H0ad@#4?6dr!iLF+QbdN2Y<N-NK3;wmzwL-dp9>S6yd3N&}N8
z4}JLb_zBo}<O2H-H`73eMl~m$wLbyJ6n%Zu`eGWWeydo6>Sa%;oi|!8tozrZN7BH$
z)1%6)m~|Y^IBV}SY)=}HtWNH@tif?uw_f{)d$y#3i0~3Si@iRkc0Orn6V9sK(loI3
z^&ZZ*+cDVib<E$L=cIuVE`v>PrXPj>+jYru{M0nixMjg56HXjeJI{3H{@Yu(*rtJd
z(dBk_sd^MjR+zTl+h0-MtG%bVJ@to6i#XNT-#rb)A8h!FU+b@S{^{AWhH_c+G+=hc
z|GSNWKfHZxgoA7SG~g80@1EaTKiGblVaI{h6#S4JI3jcLQ#&t}b52>BR4fhHNe(65
z5B7z#D$WbPkS`4cIz1nGxrZ-ow05iA%%oHRhA(?}^~w>o^Hn<+_4j%lkqU;qf3hL{
zln*SB7HKp6c`Eoe$j`ocV;|V=VUu#j@23LOV#|Hz`yN(1kM;K1Zkwx}O9kZyKCRpP
z&>>j)l4B1q|5Om*H`DxKokP&~ap&k#`%}R`F3YUOA3CUZervux{%`(Sn+hgYG~f5s
z_W(3lf26j>(o|4qN7ILXjSs*{ZF_WBIy)82X|cxO?*MPL^Iq%E_$#=eLn?57@?`0`
ztNURkue)vA3{3?eJ6}zV?ztb%|7L0H&?go2OscT2#glz%=fi%i>a=f1t5h)j`hkk&
zKkbFp4R4<Jk)#5V)zdouGxoxZ1A96ht&s{W+s{lMX|z}EJlQ*E4P15<PXz|4d-s&9
zxd*OnRs8(Ye5v4bjK}-Rr*^}9R`IjPB&LAjsju#I=(=0&{Mq|a@7f!^NdceQ<zEmf
z-36mNCpEwEBn4dj*y4j@fn9KN{jUq=-%bIEWjv}D-o6vk=hc?nx+wn9i4^eo+?cT^
z?s!4VvDU?<2U9?$-SRfj!3&zLnKnqYBLyTZuQ%?KftT9(wvC4v*cEV30WAuRGjo~o
zH@te|QNfwBQ$S3uu!5&fZijKBZpq(BQ-E#TdX8P=wyB+mTgK{j!_9+Iz=Xap-#g6S
z3OAPD{V2R=3h*C1sei(yEwFi)8@4sur+`7`4mVyE-lBGX?zc{^q7`*gfM-vG*MWOH
zVX=*&hDR!;fbl8gB)4Kc;3W4ESFV*z0o{ukzG&FRL+!lXcEzU+c$$<9R@?g3t@CXo
zoPMBMsfW?Ypk6W2t1CS=!kz8vjK1_b8FYSXQRcSy2DS5fMSe3q=if{QU1rEn)Nx#|
zIB&M*m_8ShL931Twk^H24o**s*;eLgGN`}GJ#|Q<b!zANZhq-qZvN(Euw$e7g=zWL
zLeD-Oy4GBk3_b)MD0F$U4AzqNKOuKZ2I9MItZQ6eqjvu9NVm?VgB_B=vIF;y$=0rd
zN#Rw#_8O85JX4J}jdfcIqnl3*J>M%CTyGS+bENYMwex}-M2_*B*)$m#T-%f{X#6tR
zw7b1qpjk4gG5=n%zLVVH))F%s7Ot2K=1#LcSlo4~+WEpaPwh;aXpju5cGzw`c=KXd
zdEV%>HQ$p!&%Fmq#`!OT?IS$>_C_TE+s2bGzQ4Cn?L1<SB_YvA?j`{X>((bzD=dH&
z>|#gkx|jsOwJ86(-R48bmAm^aKb{24eESr7Fvm^p{Nl@YMV*^(NdkX02p!eu!yMT5
zQTtUXE0e&;MHOWar_P2CkBMSV%u52DcgC(;W;9FfyyLlhdkoKSp9IQ}snq_7l`A~?
z>}v0=1Czjm@Bf@wS718a-DmfJ8r>A-7yT`}aot7jeB{6@hbNV8m;^#c88wZaFcn^H
z)Xrg6l_YR{b86XQeWt+mwuM1(sU&cwaHUP3`N?YMDO>oC88kW}5p=3>Akg1^A{;h0
z!fo^WM6lex%9N6R6JW|km%f)?Bm#$(|2XXVJWlQW<>1tU=cDC`z<BxYhWP^=V1*eg
zT7Es62zsq=U$#!WF)(_5`MU4+Bm!{k<&1i7Mys9I{Kzz@zQ21S@L3gAWX_b4@VvQ6
z^#wB%!8XYfkD7KPpvRT9{X0)c1TGcK`;T_8Q#;?e-NzE<-n|n+RM!Ch#FAm~gyf;T
zQ@cdqc1l`2N-`97Z#=Z`WAj9C*niHkYsH7Cod?}6UwFr$(up8!@IzU7r$KP7WZ1G6
zMvC{1ZI(FnqYZR>vfg83LIU`9w!`*_oi=LcM;`{EPij0*01abr_jmTUfnzJ0e;RZ<
z0d(JbdsDw6gP_s18|&Q8CxCkcMJ{!B3{pF9y6Vm%{rB!j0Ox|b&2pbO1lo*#XuVsO
z06vdf*Z<q0p>T&)=|k)1CjgtJb}g$`8?JUf^_1Gxs&pTf05;xB@wX1MgHDIztCh4)
z0CP|KG*9|A0`e~hb-2|r0Sq#Gv&y-~D7EvfT}_^BFI6)EjIC&U?aga@XbDo!?kbZ2
zmfpG<W;4kF4vPysSIsB^yh@1Rcb6ZlcK&tA>PzY+N5+E>37@%3cPGG2^|sIL6A}+<
zRF69QzmSP=z2Z>s&4KY?c-7FeN8F@p=VcdL_rFGI-thp!*YB;KO@U{pM+BAM91p%8
zI<b38WoKynynmS4vUuRoqfkS87Z<hjwU<WrEM0DFJa~4@Bh;jyD-7#6zf!VIJTU#z
zqGy}uGvJCPCx_p)j0bb4f3Gv&d6wFF+{e92_OLXI2R&z`4*2kFF6_8?KxlN?crfhX
z!1)quH`sDhh|?0Ic;ISY%IR8%`D*8P-{Uv-cYhlPDm1%wBe=~%ShU0ZwO@ndz&$7`
zIKy%gj5uhrspr)=u>WPzlx{;8tDW~9zQ*}-@a{MezCYj2`#YAxjss5IH&`DB@;BOH
zbM&4&+;-Y|s(3*h5Syh0l`FVh?R@akeP*<&Zx;v7J!xtY;<*xr8k+Q|*gFn95EYs;
z?DHyE)8f{d<Ti1j!-6`JYrWT~ohN>B!KEcDD#U@K=eK{fD!UdY-W+wSRiQZGZm|0L
zuE%TP3bXjOFA~0k(@z}Bm)*Kf?fmiTM$QXA1%C&Pi^uNR+;{`DbS~It_?7Qq;<Hfa
zi$)vam!2QvE*$+1qF$yK0IxQxomZY;`mL?i8pZPot*1!$c)%|G+n;iq^BqLrn%4%b
z_Jn&+G|#tl{C6<DOWi}8=WJFx-~8bU_|>n=cX0cSnWVPUR%p7%!p*DccMz~=hsR8(
zZLr&s+Gplf`wqIiKGI$`b-UVm=$*UA6fK$_3l`rw{k{F#9WccAQ0cQDW5KH#l|$!R
zdqJ+9?_Qf{v7q7n`jv*h@lrcKy~0MPHX~2Of;D9=QkoRn1x=p}XnJ>7EEr>9)Rf<~
z3%;~`zp=I~7PP;5VbQoYyVcHH4?O(Q`N+6f;JNRiPwP2*;H?_f-aZ)^3;M~6M#Way
z3;(k5+#Ay|7A&eaD*5A$y=v#PZ+~|s;d7-}@bJNnB^RpghkJH3S@5`UEO^&tV7&vu
z`=NKhpiYMqzJX59rWA-<?yYv7`-|;9jx`>917~~934WS>09G3jaWU}1H?YjE-t7hO
zApAS!edi&EzJZO4Mko6(IHY#|dmRhcTRj(j1E)f~8wZy1fq(yBd)FQpRh7p_T$6QG
zNY}*3VohCkeT^1nn|dZ=6<u8)*oKc3G67Q(UubPAif#BP$MU17t&O?$@so=Vii9q3
zKwbujZ2K|HM`5j|@*3cilG3^i=lryLE_dz#GcYs$<_`|%-rv3FoZtOD&b{aNJ-@Py
z-1x{EvT=OK!zV7q)8w>QuGP(|A={o`{zKT2gR*?_J<gun9qv~{Lft+bS*bPBQNFte
zz4e<KqKnzKG~`1g9X4<K6@8Z);^+JOd8r-=vi$M`Y4DfJFISU<Nq!f{+(@9tjLdG6
z{#{L08VZiBUReLWb7kw=o~hO3(};EdOZ@L)Sw8yHPwa4y+Ez{8+Y{Dp*!d&$c=472
zoujKsc*J#M%=AP`kFD9hbwM@hr(N)uYhNYG^4H(k@@~khqpL|KId$iOOGoMaYlFW?
z9aK#&>`IxG7nVeaJoTRgnqJjJlX`npvMEWH?|#$RlbO>htH=j4%SR9RWis{j`q3-&
zTot*Ud*bnwxMYgB8N|&XZU%8Ph?_y&Ok=ni^?0P_?mybsDfb^$=SNfHvj$qp2fesd
z<q&n%xT{Clx*M#2Iz*?j6RpMpxxTQ7i<ZYfJ<0XP$>X2AnSLle5YU^*-j?Z<u6Nbr
z_LHkWIHg>v+2PneXop=@7pxhV7`LP47p3#KI=F|e_rrxaH!+T{rEzZ1U&oDeD=A`;
zMr|ACwvz9c4*U{Jh;vIn&g=ukuU+MJ_f&S=hhKQrS4w9*!s<mwSqc8C=y54sEw6j^
zXUD~UI9gpPj!Ybfn%y=0mW>zUM;&dP*dP0HK(xC)yt(0m9@su!Qp6oF2?YIadLAtu
zpBG8>SB^K(+n1ruyt{~VHyoCZL%c26F@1t%I?!vq@I1x)xk_EhE{c8t<P74{;Vn(g
zUWNywy>t~>Y0fz0?{zJ;TeTdSA=P7Nv0F<1V<pe?0z9@o-n;R5```k`$6($t<mD;i
z`DNg^T3#MCJ1}_&YX|T-FitG*4-SC820t99mFlm^1@OZ~b1l<>JQ?1%OefyoS?lWd
zw7Kgy`wz1BTi`t7kle1uIuz@a@_q&V(>fec?F&q9PENHfKZqmB{Bi0WjV{0D?C%oi
zW>drwA&v-fM2I6o91-G(5J!YKBE%6PjtFr?h$BKA5#optM}#;c#1SEmC^xc-M2s3x
zQ5yR%8dzDnWW~%XvMRk>KxqGzdR+Gh{d$h6A~(hj8W@n1LPvVLefx>Fij4T<&MtA`
zsdW5HF)97MtH}9Yiz8=sNuxg5+w-<}sH)!=<l?p6iD`878~J^1R+>o6JN<lKnoa4(
z_fKCKd(K1(^}+K3JEYTF<xvlA$TE?MSB8B>XQk6CeNUV_e$YflXf_3{IGRouAGjQN
zG0sFLef!Y#RUR4iK*uRNP3ug={n@9dpLsQd_SPbf2ysM+BSIVz;)oi<5vi{C{+X?(
zz<QrUEgE5G?9KuCI;HhE5wqlU(6~43h3V$z))0E+Ucb+-`JeXfq7Si^J?VO#qpwr4
z_3wN8<mx6$`!7Ba(3`kPmg&UnT)l2Txmxd(a;0X6jnkO?fIPXneBk(`98iikyDG#x
z*gjs_99gR2yQ6C~fE!6)%jB^nhLuY38?{HJ@f#3#3ZrveoJj0W*7LT-i5U2P+Q2Wd
z^yLR&e{hKJ6`zj}(D8lAflfb*JwM}YseLID_T9ENzDDd{Th-O<sqb_)UTKt6TlE&(
zZEX(7=O0WGKrfn(!-%zV-l1>pp<2}*#C1X(X7o&^A4(U)JS1|4WjfF=uOQ2GO6@q?
zak*SE3vIuAgPX0CyT~?SrJ8^5#CrvL2PNb?P&?oJd<S~YH-m&P;`^tq^9`!`rJ&VJ
zPe2bPox_UDb3kdHU&QBUzyRChO4M+P&5cLb8reJ=JOLYlsArCrM=*#pdg<Koc@gzV
zSJZGbqW?^P*|NU_|7i@dOb1*G7+{$W<igp@*7u)vR^rgmf0||M-uGJ#$C(%;A}8v#
zt^KE5KA48UyrKE?0fUp`W*YxKBj!VrJfxph0{Hiv@V6}67sSmVZl>xnjV@!~?A5b-
zoS>M02lMY>{vFJ}gZXzb{|@Hg!TdXze+TpLVE!G<zk~UAF#itb-@*Jln15$S65?h|
zBxqXQPsL+1>B3Hr#=lrsNhVIeFgol&CN)Ndt}vdiBoAhPXA1GmqJua5peg^Tl4J#L
zD#(w@qNcqYc6W)YB#9Z__wD~Wi{|*;+}&+lC7IA;{rI}E*);y<Gc)SED@pLpPxj^=
z&ZY$$-X5A?Q9%xQbkz5JEQj7H4vt=vQ9;(uy58aY%{lbB+yB_#Yg+}m^5%rP5p_BA
z^Rj&&yFx2S;l<)*y(4m||DK_bwtJ?6^v}8Y`sQzP=?lTfbcmZl+zjGo5I55pZU)wS
zVI9}7#`5|JtmFF5u}lYXAI50QbU$0daq55Vx1KlLDdP^*`O(xoMmj6`pm#_)Qj3>@
zbxY7gha3cU%5EGG<E6lUYSv5b=fZiR%gmivGCwTb{13zn=*$P3@ssN*YuTFoN#>p@
z<KhCizDoNq##K04yp)KG3ljJ*vKaK$_&4b|MY{W~Za=x2;FNNu#t+!P9?JF{Z)%R$
zGapmdv&ozveI+91=z6B=4?knkmY`pf!C}SaIUxE&Eg!MxCaImZbv`M@{1eGc9*d5%
zQi(2|vhyVwtQ2vA%)k7iUc1?^=u<|g^gw<232h#<%MSSbiSfac?=$Yh>D&c6Z4A3E
z?mbS=>0*D9-p}(LcE03*<^J%?+q_V-A38l75AeCoJB%G@RSwAgEt3S8FE<?*m(JsA
znOwnO<bGBUryQz;?>z8c;j_-No`c_w-6|cgI4ICB?H0>);&skim&+Bi(DuW{WmzeA
zk!`|K(NFvEUL<pW1D&!p&PQw)<XS7_-LU=pG8>WT*X44dY;V{ZUBiu(=k-P_<M<uZ
z-l+NK+K-w20XqY7g5~|e0VVx4ZI}8H+v~50546?yu4F6ibgV-?U|;^e62RYAlsbpF
z8N|&XZsra-7`FOs`NL^+NJqrYAZ`Y6Gl-i(+zjGo5I2Lk8N|&XZU%8Ph?_y&4B}=G
zH-oqt#LX=JHR5K@d{nn(?-P0Sk&Ty4i~5z5_%-D_(hlX(sfixX54%}LzHC2q?%RX%
zX~fVezxw868OiOulW5}eY3`fGGf}(Bi2qsdm_dUJ=o=(<*sTR+WKEB>P2Cd;=;}*<
zPntinjGPGgD^2n{Mg1a%-_Gt@My4(qxY;}96#YE&={+7LrDXQ?g<)g-3u!^0fuqMJ
zl#=r=uPS`Du#kRzE8^XVb)}?}rr#LNjMLQpp{It#1eTH~qx3W4zCTTOc1=1H+ZS;&
zh?_y&Ok=ni`FgLt;{f~<*ovuieM}yI(l&7bj=Mi@0Oy-Qva}`g&DQ(-sjYtpea4<6
z#4SiVhZUFSfVlqc!$)j7D%EdW$IZy?mKdkf-2HK{ay?`Bg9j+;nPcKT%+Y}CXGoUv
zAI1SCj4mL_G9Bo(@t9>gxUaL9_LD1<Q_7V(KbnevthJI4ju{W86fb}EIMWaJR933_
zO??`p1ATRNy}|xBjRT_J1Uu`?<p=KONO5D9@mre1Y+et@FE#?_@SFGY2SyA0hvxt#
z=<8W}fA~Wfw^}Lp6V}HgT9&H^oKmio@~`dBWO5V48LH+d4BugVv6}3;@_EruX!&3#
zb2);s+ZsOzc?CT)n7nFvd!04I8=7pV=~<8cbJt9pI&MnfoB;nZesj!Sn4^}T1M(H9
zeP^jp@xj}XmXL=??R@j|4rq&+#<gyqK{Y$zo`>^AE)w@;A=&|$F9zjG=OM2wTJf5h
zLCYf<n*b-MW7Q`AWq7-GcX7VPV%8qq*DIUVhYS5^azieDfqZ@>{m#y8PJJG3&Iftf
z-0xXlr`Qt;kZ@Q&FAtZTw|7O{pQQGrRF1{=VctDC-!h$eovYW|POgF-Qm$a$;&q1G
zC2p@eYjh2wV9mUxYds~{x5=+B3t#lwg!xOCy|^TV+f;kVRLIA^-M#bg5T$!td)`_)
zwRVn`%v^nb(Rr#whsTiv$N}U4asWAi9B4TQ)ZYJ%{QU~nZCl<V-ocqTp!ojJ=eIBx
z?cE<c0`piM5HzuDr*UhT3kUP-?#{1;_x_g6ul?uG<@0O#`wWwK7)SHJ>+R?7D^d-T
YzX!)QkwdLZF6WP})>*!A%5|3i1F>xE!2kdN

literal 0
HcmV?d00001

diff --git a/test/python/srvo3_transp.py b/test/python/srvo3_transp.py
index dcbf613b..04893c25 100644
--- a/test/python/srvo3_transp.py
+++ b/test/python/srvo3_transp.py
@@ -20,11 +20,12 @@
 ################################################################################
 
 from numpy import *
+from h5 import HDFArchive
 from triqs_dft_tools.converters.wien2k import *
 from triqs_dft_tools.sumk_dft import *
 from triqs_dft_tools.sumk_dft_tools import *
 from triqs.utility.comparison_tests import *
-from triqs.utility.h5diff import h5diff 
+from triqs.utility.h5diff import h5diff
 
 beta = 40
 
@@ -32,7 +33,7 @@ Converter = Wien2kConverter(filename='SrVO3', repacking=True)
 Converter.convert_dft_input()
 Converter.convert_transport_input()
 
-with HDFArchive('SrVO3_Sigma.h5', 'a') as ar:
+with HDFArchive('SrVO3_Sigma_transport.h5', 'a') as ar:
     Sigma = ar['dmft_transp_input']['Sigma_w']
     SK = SumkDFTTools(hdf_file='SrVO3.ref.h5', mesh=Sigma.mesh, use_dft_blocks=True)
     SK.set_Sigma([Sigma])
@@ -47,4 +48,4 @@ SK.hdf_file = 'srvo3_transp.out.h5'
 SK.save(['seebeck','optic_cond','kappa'])
 
 if mpi.is_master_node():
-    h5diff("srvo3_transp.out.h5","srvo3_transp.ref.h5") 
+    h5diff("srvo3_transp.out.h5","srvo3_transp.ref.h5")

From 7d02483a54c1ede0f5ac50a90aa03c6b7d645583 Mon Sep 17 00:00:00 2001
From: Alexander Hampel <ahampel@flatironinstitute.org>
Date: Fri, 22 Apr 2022 18:29:38 -0400
Subject: [PATCH 11/11] change NiO tutorial scripts to reflect changes to sumk

---
 doc/tutorials/images_scripts/NiO_local_lattice_GF.py | 8 ++++----
 doc/tutorials/images_scripts/nio.py                  | 4 ++--
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/doc/tutorials/images_scripts/NiO_local_lattice_GF.py b/doc/tutorials/images_scripts/NiO_local_lattice_GF.py
index e8683ce3..584def07 100644
--- a/doc/tutorials/images_scripts/NiO_local_lattice_GF.py
+++ b/doc/tutorials/images_scripts/NiO_local_lattice_GF.py
@@ -12,9 +12,9 @@ import warnings
 warnings.filterwarnings("ignore", category=FutureWarning)
 
 filename = 'nio'
-SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False)
-
 beta = 5.0
+SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False, beta=beta)
+
 
 # We analyze the block structure of the Hamiltonian
 Sigma = SK.block_structure.create_gf(beta=beta)
@@ -79,12 +79,12 @@ for ik in mpi.slice_array(ikarray):
         add_g_ik.zero()
         add_g_ik << SK.downfold(ik, 0, bname, G_latt_KS[bname], gf, shells='csc', ir=None)
         gf << gf + add_g_ik
-        
+
 G_latt_orb << mpi.all_reduce(
                 mpi.world, G_latt_orb, lambda x, y: x + y)
 
 mpi.barrier()
 
-if mpi.is_master_node(): 
+if mpi.is_master_node():
     ar['DMFT_results']['Iterations']['G_latt_orb_it'+str(iteration_offset-1)] = G_latt_orb
 if mpi.is_master_node(): del ar
diff --git a/doc/tutorials/images_scripts/nio.py b/doc/tutorials/images_scripts/nio.py
index 88268fb3..da7f8440 100644
--- a/doc/tutorials/images_scripts/nio.py
+++ b/doc/tutorials/images_scripts/nio.py
@@ -17,9 +17,9 @@ warnings.filterwarnings("ignore", category=FutureWarning)
 
 filename = 'nio'
 
-SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False)
-
 beta = 5.0
+SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False, beta=beta)
+
 
 Sigma = SK.block_structure.create_gf(beta=beta)
 SK.put_Sigma([Sigma])