From b3b199bf40a31d77cc93baef928fbbcf07a037f9 Mon Sep 17 00:00:00 2001 From: Manuel Zingl Date: Fri, 19 Dec 2014 11:53:06 +0100 Subject: [PATCH] Restore everything which was lost in rebase. --- python/converters/wien2k_converter.py | 197 +++++++--------- python/sumk_dft_tools.py | 319 +++++++++++--------------- test/SrVO3.h5 | Bin 863936 -> 838840 bytes test/SrVO3_Sigma.h5 | Bin 240448 -> 240448 bytes test/srvo3_transp.output.h5 | Bin 4792 -> 13272 bytes test/srvo3_transp.py | 8 +- 6 files changed, 224 insertions(+), 300 deletions(-) diff --git a/python/converters/wien2k_converter.py b/python/converters/wien2k_converter.py index a3c691b5..ed415d09 100644 --- a/python/converters/wien2k_converter.py +++ b/python/converters/wien2k_converter.py @@ -54,7 +54,7 @@ class Wien2kConverter(ConverterTools): self.bands_subgrp = bands_subgrp self.transp_subgrp = transp_subgrp self.fortran_to_replace = {'D':'E'} - self.vel_file = filename+'.pmat' + self.pmat_file = filename+'.pmat' self.outputs_file = filename+'.outputs' self.struct_file = filename+'.struct' self.oubwin_file = filename+'.oubwin' @@ -373,166 +373,137 @@ class Wien2kConverter(ConverterTools): Reads the input files necessary for transport calculations and stores the data in the HDFfile """ - - #Read and write files only on the master node if not (mpi.is_master_node()): return # Check if SP, SO and n_k are already in h5 ar = HDFArchive(self.hdf_file, 'a') - if not (self.dft_subgrp in ar): raise IOError, "No %s subgroup in hdf file found! Call convert_dmft_input first." %self.dft_subgrp + if not (self.dft_subgrp in ar): raise IOError, "convert_transport_input: No %s subgroup in hdf file found! Call convert_dmft_input first." %self.dft_subgrp SP = ar[self.dft_subgrp]['SP'] SO = ar[self.dft_subgrp]['SO'] n_k = ar[self.dft_subgrp]['n_k'] del ar - # Read relevant data from .pmat file - ############################################ - + # Read relevant data from .pmat/up/dn files + ########################################### + # band_window_optics: Contains the index of the lowest and highest band within the + # band window (used by optics) for each k-point. + # velocities_k: velocity (momentum) matrix elements between all bands in band_window_optics + # and each k-point. + if (SP == 0 or SO == 1): - files = [self.vel_file] + files = [self.pmat_file] elif SP == 1: - files = [self.vel_file+'up', self.vel_file+'dn'] + files = [self.pmat_file+'up', self.pmat_file+'dn'] else: # SO and SP can't both be 1 - assert 0, "convert_transport_input: reading velocity file error! Check SP and SO!" + assert 0, "convert_transport_input: Reading velocity file error! Check SP and SO!" - vk = [[] for f in files] - kp = [[] for f in files] - bandwin_opt = [] - + velocities_k = [[] for f in files] + band_window_optics = [] for isp, f in enumerate(files): - bandwin_opt_isp = [] - if not os.path.exists(f) : raise IOError, "File %s does not exist" %f - print "Reading input from %s..."%f - - with open(f) as R: - while 1: - try: - s = R.readline() - if (s == ''): - break - except: - break - try: - [k, nu1, nu2] = [int(x) for x in s.strip().split()] - bandwin_opt_isp.append((nu1,nu2)) - dim = nu2 - nu1 +1 - v_xyz = numpy.zeros((dim,dim,3), dtype = complex) - temp = R.readline().strip().split() - kp[isp].append(numpy.array([float(t) for t in temp[0:3]])) - for nu_i in xrange(dim): - for nu_j in xrange(nu_i, dim): - for i in xrange(3): - s = R.readline().strip("\n ()").split(',') - v_xyz[nu_i][nu_j][i] = float(s[0]) + float(s[1])*1j - if (nu_i != nu_j): - v_xyz[nu_j][nu_i][i] = v_xyz[nu_i][nu_j][i].conjugate() - - vk[isp].append(v_xyz) - - except IOError: - raise "Wien2k_converter : reading file %s failed" %self.vel_file - bandwin_opt.append(numpy.array(bandwin_opt_isp)) - - print "Read in %s file done!" %self.vel_file + if not os.path.exists(f) : raise IOError, "convert_transport_input: File %s does not exist" %f + print "Reading input from %s..."%f, + R = ConverterTools.read_fortran_file(self, f, {'D':'E','(':'',')':'',',':''}) + band_window_optics_isp = [] + for ik in xrange(n_k): + R.next() + nu1 = int(R.next()) + nu2 = int(R.next()) + band_window_optics_isp.append((nu1, nu2)) + n_bands = nu2 - nu1 + 1 + velocity_xyz = numpy.zeros((n_bands, n_bands, 3), dtype = complex) + for _ in range(4): R.next() + for nu_i in range(n_bands): + for nu_j in range(nu_i, n_bands): + for i in range(3): + velocity_xyz[nu_i][nu_j][i] = R.next() + R.next()*1j + if (nu_i != nu_j): velocity_xyz[nu_j][nu_i][i] = velocity_xyz[nu_i][nu_j][i].conjugate() + velocities_k[isp].append(velocity_xyz) + band_window_optics.append(numpy.array(band_window_optics_isp)) + print "DONE!" # Read relevant data from .struct file - ############################################ - if not (os.path.exists(self.struct_file)) : raise IOError, "File %s does not exist" %self.struct_file - print "Reading input from %s..."%self.struct_file + ###################################### + # lattice_type: bravais lattice type as defined by Wien2k + # lattice_constants: unit cell parameters in a. u. + # lattice_angles: unit cell angles in rad + + if not (os.path.exists(self.struct_file)) : raise IOError, "convert_transport_input: File %s does not exist" %self.struct_file + print "Reading input from %s..."%self.struct_file, - with open(self.struct_file) as f: + with open(self.struct_file) as R: try: - f.readline() #title - temp = f.readline() #lattice - latticetype = temp.split()[0] - f.readline() - temp = f.readline().strip().split() # lattice constants - latticeconstants = numpy.array([float(t) for t in temp[0:3]]) - latticeangles = numpy.array([float(t) for t in temp[3:6]]) - latticeangles *= numpy.pi/180.0 - print 'Lattice: ', latticetype - print 'Lattice constants: ', latticeconstants - print 'Lattice angles: ', latticeangles - + R.readline() + lattice_type = R.readline().split()[0] + R.readline() + temp = R.readline().strip().split() + lattice_constants = numpy.array([float(t) for t in temp[0:3]]) + lattice_angles = numpy.array([float(t) for t in temp[3:6]]) * numpy.pi / 180.0 except IOError: - raise "Wien2k_converter : reading file %s failed" %self.struct_file - - print "Read in %s file done!" %self.struct_file - + raise "convert_transport_input: reading file %s failed" %self.struct_file + print "DONE!" # Read relevant data from .outputs file - ############################################ - if not (os.path.exists(self.outputs_file)) : raise IOError, "File %s does not exist" %self.outputs_file - print "Reading input from %s..."%self.outputs_file + ####################################### + # rot_symmetries: matrix representation of all (space group) symmetry operations - symmcartesian = [] - taucartesian = [] - - with open(self.outputs_file) as f: + if not (os.path.exists(self.outputs_file)) : raise IOError, "convert_transport_input: File %s does not exist" %self.outputs_file + print "Reading input from %s..."%self.outputs_file, + + rot_symmetries = [] + with open(self.outputs_file) as R: try: while 1: - temp = f.readline().strip(' ').split() + temp = R.readline().strip(' ').split() if (temp[0] =='PGBSYM:'): - nsymm = int(temp[-1]) + n_symmetries = int(temp[-1]) break - for i in range(nsymm): + for i in range(n_symmetries): while 1: - temp = f.readline().strip().split() - if (temp[0] == 'Symmetry'): - break - - # read cartesian symmetries - symmt = numpy.zeros((3, 3), dtype = float) - taut = numpy.zeros(3, dtype = float) + if (R.readline().strip().split()[0] == 'Symmetry'): break + sym_i = numpy.zeros((3, 3), dtype = float) for ir in range(3): - temp = f.readline().strip().split() + temp = R.readline().strip().split() for ic in range(3): - symmt[ir, ic] = float(temp[ic]) - temp = f.readline().strip().split() - for ir in range(3): - taut[ir] = float(temp[ir]) - - symmcartesian.append(symmt) - taucartesian.append(taut) + sym_i[ir, ic] = float(temp[ic]) + R.readline() + rot_symmetries.append(sym_i) except IOError: - raise "Wien2k_converter: reading file %s failed" %self.outputs_file - - print "Read in %s file done!" %self.outputs_file + raise "convert_transport_input: reading file %s failed" %self.outputs_file + print "DONE!" - - # Read relevant data from .oubwin/up/down files - ############################################ + # Read relevant data from .oubwin/up/dn files + ############################################### + # band_window: Contains the index of the lowest and highest band within the + # projected subspace (used by dmftproj) for each k-point. if (SP == 0 or SO == 1): files = [self.oubwin_file] elif SP == 1: files = [self.oubwin_file+'up', self.oubwin_file+'dn'] else: # SO and SP can't both be 1 - assert 0, "Reding oubwin error! Check SP and SO!" - bandwin = [numpy.zeros((n_k, 2), dtype=int) for isp in range(SP + 1 - SO)] - + assert 0, "convert_transport_input: Reding oubwin error! Check SP and SO!" + + band_window = [numpy.zeros((n_k, 2), dtype=int) for isp in range(SP + 1 - SO)] for isp, f in enumerate(files): - if not os.path.exists(f): raise IOError, "File %s does not exist" %f - print "Reading input from %s..."%f + if not os.path.exists(f): raise IOError, "convert_transport_input: File %s does not exist" %f + print "Reading input from %s..."%f, + R = ConverterTools.read_fortran_file(self, f, self.fortran_to_replace) - assert int(R.next()) == n_k, "Wien2k_converter: Number of k-points is inconsistent in oubwin file!" - assert int(R.next()) == SO, "Wien2k_converter: SO is inconsistent in oubwin file!" - + assert int(R.next()) == n_k, "convert_transport_input: Number of k-points is inconsistent in oubwin file!" + assert int(R.next()) == SO, "convert_transport_input: SO is inconsistent in oubwin file!" for ik in xrange(n_k): R.next() - bandwin[isp][ik,0] = R.next() - bandwin[isp][ik,1] = R.next() + band_window[isp][ik,0] = R.next() + band_window[isp][ik,1] = R.next() R.next() - - print "Read in %s files done!" %self.oubwin_file + print "DONE!" # Put data to HDF5 file ar = HDFArchive(self.hdf_file, 'a') if not (self.transp_subgrp in ar): ar.create_group(self.transp_subgrp) # The subgroup containing the data. If it does not exist, it is created. If it exists, the data is overwritten!!! - things_to_save = ['bandwin', 'bandwin_opt', 'kp', 'vk', 'latticetype', 'latticeconstants', 'latticeangles', 'nsymm', 'symmcartesian', - 'taucartesian'] + things_to_save = ['band_window', 'band_window_optics', 'velocities_k', 'lattice_type', 'lattice_constants', 'lattice_angles', 'n_symmetries', 'rot_symmetries'] for it in things_to_save: ar[self.transp_subgrp][it] = locals()[it] del ar diff --git a/python/sumk_dft_tools.py b/python/sumk_dft_tools.py index 087a9bb9..82a4e12e 100644 --- a/python/sumk_dft_tools.py +++ b/python/sumk_dft_tools.py @@ -1,4 +1,3 @@ - ################################################################################ # # TRIQS: a Toolbox for Research in Interacting Quantum Systems @@ -501,72 +500,72 @@ class SumkDFTTools(SumkDFT): """ Reads the data for transport calculations from the HDF file """ - - thingstoread = ['bandwin','bandwin_opt','kp','latticeangles','latticeconstants','latticetype','nsymm','symmcartesian','vk'] + thingstoread = ['band_window','band_window_optics','velocities_k','lattice_angles','lattice_constants','lattice_type','n_symmetries','rot_symmetries'] retval = self.read_input_from_hdf(subgrp=self.transp_data,things_to_read = thingstoread) return retval - def cellvolume(self, latticetype, latticeconstants, latticeangle): + def cellvolume(self, lattice_type, lattice_constants, latticeangle): """ Calculate cell volume: volumecc conventional cell, volumepc, primitive cell. """ - a = latticeconstants[0] - b = latticeconstants[1] - c = latticeconstants[2] + a = lattice_constants[0] + b = lattice_constants[1] + c = lattice_constants[2] c_al = numpy.cos(latticeangle[0]) c_be = numpy.cos(latticeangle[1]) c_ga = numpy.cos(latticeangle[2]) volumecc = a * b * c * numpy.sqrt(1 + 2 * c_al * c_be * c_ga - c_al ** 2 - c_be * 82 - c_ga ** 2) det = {"P":1, "F":4, "B":2, "R":3, "H":1, "CXY":2, "CYZ":2, "CXZ":2} - volumepc = volumecc / det[latticetype] + volumepc = volumecc / det[lattice_type] return volumecc, volumepc - def transport_distribution(self, dir_list=[(0,0)], broadening=0.01, energywindow=None, Om_mesh=[0.0], beta=40, with_Sigma=False, n_om=None): - """calculate Tr A(k,w) v(k) A(k, w+q) v(k) and optics. - energywindow: regime for omega integral - Om_mesh: contains the frequencies of the optic conductivitity. Om_mesh is repinned to the self-energy mesh - (hence exact values might be different from those given in Om_mesh) - dir_list: list to defines the indices of directions. xx,yy,zz,xy,yz,zx. - ((0, 0) --> xx, (1, 1) --> yy, (0, 2) --> xz, default: (0, 0)) - with_Sigma: Use Sigma = 0 if False + def transport_distribution(self, directions=['xx'], energy_window=None, Om_mesh=[0.0], beta=40, with_Sigma=False, n_om=None, broadening=0.01): + """ + calculate Tr A(k,w) v(k) A(k, w+Om) v(k). + energy_window: regime for omega integral + Om_mesh: mesh for optic conductivitity. Om_mesh is repinned to the self-energy mesh! + directions: list of directions: xx,yy,zz,xy,yz,zx. + with_Sigma: Use Sigma_w = 0 if False (In this case it is necessary to specifiy the energywindow (energy_window), + the number of omega points (n_om) in the window and the broadening (broadening)). """ - # Check if wien converter was called + # Check if wien converter was called and read transport subgroup form hdf file if mpi.is_master_node(): ar = HDFArchive(self.hdf_file, 'a') - if not (self.transp_data in ar): raise IOError, "No %s subgroup in hdf file found! Call convert_transp_input first." %self.transp_data - - self.dir_list = dir_list - + if not (self.transp_data in ar): raise IOError, "transport_distribution: No %s subgroup in hdf file found! Call convert_transp_input first." %self.transp_data self.read_transport_input_from_hdf() - velocities = self.vk - n_inequiv_spin_blocks = self.SP + 1 - self.SO # up and down are equivalent if SP = 0 + n_inequiv_spin_blocks = self.SP + 1 - self.SO # up and down are equivalent if SP = 0 + self.directions = directions + dir_to_int = {'x':0, 'y':1, 'z':2} + # k-dependent-projections. + assert self.k_dep_projection == 1, "transport_distribution: k dependent projection is not implemented!" + # calculate A(k,w) ####################################### - # use k-dependent-projections. - assert self.k_dep_projection == 1, "Not implemented!" - - # Define mesh for Greens function and the used energy range + # Define mesh for Greens 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]) + mesh = None mu = self.chemical_potential n_om = len(self.omega) - print "Using omega mesh provided by Sigma." + print "Using omega mesh provided by Sigma!" - if energywindow is not None: - # Find according window in Sigma mesh - ioffset = numpy.sum(self.omega < energywindow[0]) - self.omega = self.omega[numpy.logical_and(self.omega >= energywindow[0], self.omega <= energywindow[1])] + if energy_window is not None: + # Find according window in Sigma mesh + ioffset = numpy.sum(self.omega < energy_window[0]) + self.omega = self.omega[numpy.logical_and(self.omega >= energy_window[0], self.omega <= energy_window[1])] n_om = len(self.omega) # Truncate Sigma to given omega window + # In the future there should be an option in gf to manipulate the mesh (e.g. truncate) directly. + # For we stick with this: for icrsh in range(self.n_corr_shells): Sigma_save = self.Sigma_imp_w[icrsh].copy() spn = self.spin_block_names[self.corr_shells[icrsh]['SO']] @@ -576,179 +575,133 @@ class SumkDFTTools(SumkDFT): for iL in g.indices: for iR in g.indices: for iom in xrange(n_om): - g.data[iom,iL,iR] = Sigma_save[i].data[ioffset+iom,iL,iR] # FIXME IS THIS CLEAN? MANIPULATING data DIRECTLY? - + g.data[iom,iL,iR] = Sigma_save[i].data[ioffset+iom,iL,iR] else: - assert n_om is not None, "Number of omega points (n_om) needed!" - assert energywindow is not None, "Energy window needed!" - self.omega = numpy.linspace(energywindow[0],energywindow[1],n_om) - mu = 0.0 + assert n_om is not None, "transport_distribution: Number of omega points (n_om) needed to calculate transport distribution!" + assert energy_window is not None, "transport_distribution: Energy window needed to calculate transport distribution!" + assert broadening != 0.0 and broadening is not None, "transport_distribution: Broadening necessary to calculate transport distribution!" + self.omega = numpy.linspace(energy_window[0],energy_window[1],n_om) + mesh = [energy_window[0], energy_window[1], n_om] + mu = 0.0 + # Check if energy_window is sufficiently large if (abs(self.fermi_dis(self.omega[0]*beta)*self.fermi_dis(-self.omega[0]*beta)) > 1e-5 or abs(self.fermi_dis(self.omega[-1]*beta)*self.fermi_dis(-self.omega[-1]*beta)) > 1e-5): - print "\n##########################################" - print "WARNING: Energywindow might be too narrow!" - print "##########################################\n" + print "\n####################################################################" + print "transport_distribution: WARNING - energy window might be too narrow!" + print "####################################################################\n" + # Define mesh for optic conductivity d_omega = round(numpy.abs(self.omega[0] - self.omega[1]), 12) - - # define exact mesh for optic conductivity - Om_mesh_ex = numpy.array([int(x / d_omega) for x in Om_mesh]) - self.Om_meshr= Om_mesh_ex*d_omega + iOm_mesh = numpy.array([int(Om / d_omega) for Om in Om_mesh]) + self.Om_mesh = iOm_mesh * d_omega if mpi.is_master_node(): print "Chemical potential: ", mu - print "Using n_om = %s points in the energywindow [%s,%s]"%(n_om, self.omega[0], self.omega[-1]) - print "omega vector is:" + print "Using n_om = %s points in the energy_window [%s,%s]"%(n_om, self.omega[0], self.omega[-1]), + print "where the omega vector is:" print self.omega - print "Omega mesh interval ", d_omega - print "Provided Om_mesh ", numpy.array(Om_mesh) - print "Pinnend Om_mesh to ", self.Om_meshr + print "Calculation requested for Omega mesh: ", numpy.array(Om_mesh) + print "Omega mesh automatically repinned to: ", self.Om_mesh - # output P(\omega)_xy should have the same dimension as defined in mshape. - self.Pw_optic = numpy.zeros((len(dir_list), len(Om_mesh_ex), n_om), dtype=numpy.float_) - - ik = 0 - - spn = self.spin_block_names[self.SO] - ntoi = self.spin_names_to_ind[self.SO] - - G_w = BlockGf(name_block_generator=[(spn[isp], GfReFreq(indices=range(self.n_orbitals[ik][isp]), window=(self.omega[0], self.omega[-1]), n_points = n_om)) - for isp in range(n_inequiv_spin_blocks) ], make_copies=False) - mupat = [numpy.identity(self.n_orbitals[ik][isp], numpy.complex_) * mu for isp in range(n_inequiv_spin_blocks)] # construct mupat - Annkw = [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)] + self.Gamma_w = {direction: numpy.zeros((len(self.Om_mesh), n_om), dtype=numpy.float_) for direction in self.directions} + # Sum over all k-points ikarray = numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): - unchangedsize = all([ self.n_orbitals[ik][isp] == mupat[isp].shape[0] for isp in range(n_inequiv_spin_blocks)]) - if not unchangedsize: - # recontruct green functions. - G_w = BlockGf(name_block_generator=[(spn[isp], GfReFreq(indices=range(self.n_orbitals[ik][isp]), window = (self.omega[0], self.omega[-1]), n_points = n_om)) - for isp in range(n_inequiv_spin_blocks) ], make_copies=False) - # construct mupat - mupat = [numpy.identity(self.n_orbitals[ik][isp], numpy.complex_) * mu for isp in range(n_inequiv_spin_blocks)] - #set a temporary array storing spectral functions with band index. Note, usually we should have spin index - Annkw = [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)] - # get lattice green function + # 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) + 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)] - G_w << 1*Omega + 1j*broadening - - MS = copy.deepcopy(mupat) - for ibl in range(n_inequiv_spin_blocks): - ind = ntoi[spn[ibl]] - n_orb = self.n_orbitals[ik][ibl] - MS[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb].real - mupat[ibl] - G_w -= MS - - if (with_Sigma == True): - tmp = G_w.copy() # init temporary storage - # form self energy from impurity self energy and double counting term. - sigma_minus_dc = self.add_dc(iw_or_w="w") - # substract self energy - for icrsh in range(self.n_corr_shells): - for sig, gf in tmp: tmp[sig] << self.upfold(ik, icrsh, sig, sigma_minus_dc[icrsh][sig], gf) - G_w -= tmp - - G_w.invert() - - for isp in range(n_inequiv_spin_blocks): - Annkw[isp].real = -copy.deepcopy(G_w[self.spin_block_names[self.SO][isp]].data.swapaxes(0,1).swapaxes(1,2)).imag / numpy.pi - - for isp in range(n_inequiv_spin_blocks): - kvel = velocities[isp][ik] - Pwtem = numpy.zeros((len(dir_list), len(Om_mesh_ex), n_om), dtype=numpy.float_) + for isp in range(n_inequiv_spin_blocks): + # Obtain A_kw from G_w (swapaxes is used to have omega in the 3rd dimension) + A_kw[isp].real = -copy.deepcopy(G_w[self.spin_block_names[self.SO][isp]].data.swapaxes(0,1).swapaxes(1,2)).imag / numpy.pi + b_min = max(self.band_window[isp][ik, 0], self.band_window_optics[isp][ik, 0]) + b_max = min(self.band_window[isp][ik, 1], self.band_window_optics[isp][ik, 1]) + A_i = slice(b_min - self.band_window[isp][ik, 0], b_max - self.band_window[isp][ik, 0] + 1) + v_i = slice(b_min - self.band_window_optics[isp][ik, 0], b_max - self.band_window_optics[isp][ik, 0] + 1) - bmin = max(self.bandwin[isp][ik, 0], self.bandwin_opt[isp][ik, 0]) - bmax = min(self.bandwin[isp][ik, 1], self.bandwin_opt[isp][ik, 1]) - Astart = bmin - self.bandwin[isp][ik, 0] - Aend = bmax - self.bandwin[isp][ik, 0] + 1 - vstart = bmin - self.bandwin_opt[isp][ik, 0] - vend = bmax - self.bandwin_opt[isp][ik, 0] + 1 - - #symmetry loop - for Rmat in self.symmcartesian: - # get new velocity. - Rkvel = copy.deepcopy(kvel) - for vnb1 in range(self.bandwin_opt[isp][ik, 1] - self.bandwin_opt[isp][ik, 0] + 1): - for vnb2 in range(self.bandwin_opt[isp][ik, 1] - self.bandwin_opt[isp][ik, 0] + 1): - Rkvel[vnb1][vnb2][:] = numpy.dot(Rmat, Rkvel[vnb1][vnb2][:]) - ipw = 0 - for (ir, ic) in dir_list: + # loop over all symmetries + for R in self.rot_symmetries: + # get transformed velocity under symmetry R + vel_R = copy.deepcopy(self.velocities_k[isp][ik]) + for nu1 in range(self.band_window_optics[isp][ik, 1] - self.band_window_optics[isp][ik, 0] + 1): + for nu2 in range(self.band_window_optics[isp][ik, 1] - self.band_window_optics[isp][ik, 0] + 1): + vel_R[nu1][nu2][:] = numpy.dot(R, vel_R[nu1][nu2][:]) + + # calculate Gamma_w for each direction from the velocities vel_R and the spectral function A_kw + for direction in self.directions: for iw in xrange(n_om): - for iq in range(len(Om_mesh_ex)): - if(iw + Om_mesh_ex[iq] >= n_om): - continue - - # construct matrix for A and velocity. - Annkwl = Annkw[isp][Astart:Aend, Astart:Aend, iw] - Annkwr = Annkw[isp][Astart:Aend, Astart:Aend, iw + Om_mesh_ex[iq]] - Rkveltr = Rkvel[vstart:vend, vstart:vend, ir] - Rkveltc = Rkvel[vstart:vend, vstart:vend, ic] - # print Annkwl.shape, Annkwr.shape, Rkveltr.shape, Rkveltc.shape - Pwtem[ipw, iq, iw] += numpy.dot(numpy.dot(numpy.dot(Rkveltr, Annkwl), Rkveltc), Annkwr).trace().real - ipw += 1 - - # k sum and spin sum. - self.Pw_optic += Pwtem * self.bz_weights[ik] / self.nsymm - - self.Pw_optic = mpi.all_reduce(mpi.world, self.Pw_optic, lambda x, y : x + y) - self.Pw_optic *= (2 - self.SP) + for iq in range(len(self.Om_mesh)): + if(iw + iOm_mesh[iq] >= n_om): continue + self.Gamma_w[direction][iq, iw] += (numpy.dot(numpy.dot(numpy.dot(vel_R[v_i, v_i, dir_to_int[direction[0]]], + A_kw[isp][A_i, A_i, iw]), vel_R[v_i, v_i, dir_to_int[direction[1]]]), + A_kw[isp][A_i, A_i, iw + iOm_mesh[iq]]).trace().real * self.bz_weights[ik]) + for direction in self.directions: + self.Gamma_w[direction] = (mpi.all_reduce(mpi.world, self.Gamma_w[direction], lambda x, y : x + y) + / self.cellvolume(self.lattice_type, self.lattice_constants, self.lattice_angles)[1]) / self.n_symmetries - def conductivity_and_seebeck(self, beta=40): - """ #return 1/T*A0, that is Conductivity in unit 1/V - calculate, save and return Conductivity + + def transport_coefficient(self, direction, iq=0, n=0, beta=40): """ + calculates the transport coefficients A_n in a given direction and for a given Omega. (see documentation) + A_1 is set to nan if requested for Omega != 0.0 + iq: index of Omega point in Om_mesh + direction: 'xx','yy','zz','xy','xz','yz' + """ + if not (mpi.is_master_node()): return + + assert hasattr(self,'Gamma_w'), "transport_coefficient: Run transport_distribution first or load data from h5!" + A = 0.0 + omegaT = self.omega * beta + d_omega = self.omega[1] - self.omega[0] + if (self.Om_mesh[iq] == 0.0): + for iw in xrange(self.Gamma_w[direction].shape[1]): + A += self.Gamma_w[direction][iq, iw] * self.fermi_dis(omegaT[iw]) * self.fermi_dis(-omegaT[iw]) * numpy.float(omegaT[iw])**n * d_omega + elif (n == 0.0): + for iw in xrange(self.Gamma_w[direction].shape[1]): + A += (self.Gamma_w[direction][iq, iw] * (self.fermi_dis(omegaT[iw]) - self.fermi_dis(omegaT[iw] + self.Om_mesh[iq] * beta)) + / (self.Om_mesh[iq] * beta) * d_omega) + else: + A = numpy.nan + return A * numpy.pi * (2.0-self.SP) - if mpi.is_master_node(): - assert hasattr(self,'Pw_optic'), "Run transport_distribution first or load data from h5!" - assert hasattr(self,'latticetype'), "Run convert_transp_input first or load data from h5!" - - volcc, volpc = self.cellvolume(self.latticetype, self.latticeconstants, self.latticeangles) - - n_direction, n_q, n_w= self.Pw_optic.shape - omegaT = self.omega * beta - A0 = numpy.zeros((n_direction,n_q), dtype=numpy.float_) - q_0 = False - self.seebeck = numpy.zeros((n_direction,), dtype=numpy.float_) - self.seebeck[:] = numpy.NAN - - d_omega = self.omega[1] - self.omega[0] - for iq in xrange(n_q): - # treat q = 0, caclulate conductivity and seebeck - if (self.Om_meshr[iq] == 0.0): - # if Om_meshr contains multiple entries with w=0, A1 is overwritten! - q_0 = True - A1 = numpy.zeros((n_direction,), dtype=numpy.float_) - for idir in range(n_direction): - for iw in xrange(n_w): - A0[idir, iq] += beta * self.Pw_optic[idir, iq, iw] * self.fermi_dis(omegaT[iw]) * self.fermi_dis(-omegaT[iw]) - A1[idir] += beta * self.Pw_optic[idir, iq, iw] * self.fermi_dis(omegaT[iw]) * self.fermi_dis(-omegaT[iw]) * numpy.float(omegaT[iw]) - self.seebeck[idir] = -A1[idir] / A0[idir, iq] - print "A0", A0[idir, iq] *d_omega/beta - print "A1", A1[idir] *d_omega/beta - # treat q ~= 0, calculate optical conductivity - else: - for idir in range(n_direction): - for iw in xrange(n_w): - A0[idir, iq] += self.Pw_optic[idir, iq, iw] * (self.fermi_dis(omegaT[iw]) - self.fermi_dis(omegaT[iw] + self.Om_meshr[iq] * beta)) / self.Om_meshr[iq] - - A0 *= d_omega - #cond = beta * self.tdintegral(beta, 0)[index] - print "V in bohr^3 ", volpc - # transform to standard unit as in resistivity - self.optic_cond = A0 * 10700.0 / volpc - self.seebeck *= 86.17 - - # print - for im in range(n_direction): - for iq in xrange(n_q): - print "Conductivity in direction %s for Om_mesh %d %.4f x 10^4 Ohm^-1 cm^-1" % (self.dir_list[im], iq, self.optic_cond[im, iq]) - print "Resistivity in direction %s for Om_mesh %d %.4f x 10^-4 Ohm cm" % (self.dir_list[im], iq, 1.0 / self.optic_cond[im, iq]) - if q_0: - print "Seebeck in direction %s for q = 0 %.4f x 10^(-6) V/K" % (self.dir_list[im], self.seebeck[im]) - + + def conductivity_and_seebeck(self, beta=40): + """ + Calculates the Seebeck coefficient and the conductivity for a given Gamma_w + """ + if not (mpi.is_master_node()): return + + assert hasattr(self,'Gamma_w'), "conductivity_and_seebeck: Run transport_distribution first or load data from h5!" + n_q = self.Gamma_w[self.directions[0]].shape[0] + + A0 = {direction: numpy.full((n_q,),numpy.nan) for direction in self.directions} + A1 = {direction: numpy.full((n_q,),numpy.nan) for direction in self.directions} + self.seebeck = {direction: numpy.nan for direction in self.directions} + self.optic_cond = {direction: numpy.full((n_q,),numpy.nan) for direction in self.directions} + + for direction in self.directions: + for iq in xrange(n_q): + A0[direction][iq] = self.transport_coefficient(direction, iq=iq, n=0, beta=beta) + A1[direction][iq] = self.transport_coefficient(direction, iq=iq, n=1, beta=beta) + print "A_0 in direction %s for Omega = %.2f %e a.u." % (direction, self.Om_mesh[iq], A0[direction][iq]) + print "A_1 in direction %s for Omega = %.2f %e a.u." % (direction, self.Om_mesh[iq], A1[direction][iq]) + if ~numpy.isnan(A1[direction][iq]): + # Seebeck is overwritten if there is more than one Omega = 0 in Om_mesh + self.seebeck[direction] = - A1[direction][iq] / A0[direction][iq] * 86.17 + self.optic_cond[direction] = A0[direction] * 10700.0 + for iq in xrange(n_q): + print "Conductivity in direction %s for Omega = %.2f %f x 10^4 Ohm^-1 cm^-1" % (direction, self.Om_mesh[iq], self.optic_cond[direction][iq]) + if not (numpy.isnan(A1[direction][iq])): + print "Seebeck in direction %s for Omega = 0.00 %f x 10^(-6) V/K" % (direction, self.seebeck[direction]) + def fermi_dis(self, x): + """ + fermi distribution at x = omega * beta + """ return 1.0/(numpy.exp(x)+1) diff --git a/test/SrVO3.h5 b/test/SrVO3.h5 index 64ebdab9c7f208d827be6bb958250fbb72e35dd5..1b226b0a8d2790a44f9573b6bf2b1d096f818bdf 100644 GIT binary patch delta 20635 zcmch9eLz%I`ag5;Fx+8A9ef2b8MTts615i;OO)E!tQbW%#5D)eB3}{F@RdocC9!32 zf`e&aN3E6AVifgKR*uqc1n!D2yUS*qv$_=t#+1EaR(_xJoOAA6+~)V&e||szsPmkA z&gXod=RD7I&b^#7{(s(Cw`GhwQLZ~--9P+tl;lpU?XUbsk!@Z}`Tmuc`fVeO@*}-GwH>`HB>C=M9?z0H<7&PA4o1;h(xXgzT~^+}u6Da3$G(Y0Z=vFG z*c0W>y$oqUq`7v7-C{tbnnB&LKr;R%?WlDQd03{^o{Ul3YabnY$te43H8}eYTC8=A zIKfC{jG8WC+L zQuo-WwJL%AmI(UKGchODu}4G&U)lWlSTD-x*%;>~#*JHMJ+5#|-9VAM4=s zl=n|iQV(d|?Vi9Ko9sV`quX*7Ik1fuiM^+y$le23?aWP*9j`*MVsos$Hhu21krWtB zw|M>&W+8Y6u8FR774Bp99=U?syXG0W-7~z1k;y2!J!JQRYYt;uPtBk<&+&DGQGp%} z(ELm>Huq!{y;bXT>8weIXx_kvAUz&|rVE~D$Ysy7@n|0@{|~j%E z#{D?2;D5#wo-`j0Os-#uJmDF$pAqipJ&`2cD=D?!{c9BCh68s^@z~#PN5z)RHqNNc zJCY@b8SXtX+gUR!m{x()}~US&K%mZcgPH9DjkOf`%! zA{EvCfZ-hj7D?b&hC)`#cSdEBy$>1&u{O;Q8mw|3aO!E9k7nF&h>#;(%xDoVG9}HL zKx~r@itJ9ef!LwSynRhJNXth?_t9lHoYCam#SF~gAoHrojuL)6Z#B+Qh|CZ z+pwBF$e3+d&YX#MNS4=m+xsMn6#EZ$&FGxo)7?CwV;qYZs{H3nkrb zDg~lbkEivWI@d@+qwF80?B9Qv0g{zD{T>$8iaMw>>U#qT<|%e!AEy*j1CaBhl$FG0 zl+4i^By{qmF$!6|A@WhrGjBdfI+sSpk(l2}f0yhJs+WEz4c09_5f4!R_t3@1_A6Sf zbW2HcFb$)^;g*;~kO>3ic9OYSV&RxkCGRA09%+G*h8J;uC)-K2M@p4lpIAtnN4iN$ z97v*Hk_O5)zeOGXl2l=ks!i(g3h8T>cQRC|o$T7e;Aqw)2VRy~+AOPvk?ObdN00+m z?>31gPBQW@((9^!yTsGx(My0bb0;T7lPFoRi%DwPE=jjcs;RF@<8eY6M&cfUaebB* zhTF$6&R@2)lQVDd!6~l`Z!%s$a;s<6oPO$CHBvkaPGsY%FW!p&%h$IEH2v%G=qKy0%N_a^7fG&w?rRVM_cNSqMnC&918>~6mjn#MNr z^0$0C=(O$mj%#%CJLyH1QfB`JCfEFvnB0*oz}v5L)6s22MXijHga$uukkR$&g0Ffd znawRJ8w?m@O)_&3fB(E6Q!MMn z&hTdF47fJ6UZ_>GlD^BJ?hh7msfVEY@D>UzJCn*DU5Qh^uf(n}GIGnzPbHGTy zcq{LMo^oZV4nVC-?7|^9LU71fMZIby=TI|}b09ED?PTvLJx2a$!2$w~u^_XrLzJhD zWzBF3rjFyLq%)P9!U14DcdECImowxA%4xRI6n>e+Wye`opU(Zw(a^PB5tpkPS$z?2 z=l`s0yH+kMU%-a$wMxI7NzSEm8Z1*x6pnCm1H^OM{d)EfaEqdkab+!@cKI#o(c%Jv#?%bvhGpirV9?3Cu<+m^*Nr!IiP=>nZY@r zFI8oCnaAeGIR`ZT{2V@5S&kq?{dr?9udd0JKa`98=WehUmz7UgeHFKB(hptRWpTOa zM^<0L+Z6@7w)1k?|23hOKhfL}N2DImtvu=gD7V z%BRM#gdSMQ))?b%vT8oJ9Ihwn3-rdhu!ujyrS8gNVO8q?a=wsgN^i)QlVwZMyNrI- zQ46sY%Fp$08QYas#5g9~$CEvUoE95*MVE1l*tj@DU$5jOFmB}4U4~t)pW4GN;a!mA zw7dWrasxw9!FTmJ-}Rinw_5jtp2O5lIu3JQEE=TvuY{$r!oxpy<0diYAS$&XG-F{GS)t3+J2{=_wK<#*2pKeFL$69%4 zYD+uHkO*CbNT#1ls6p(pE*d!g$T^@7T>5}Bg%qQk*yu)1p+|q?PW4nH&owHnf4Y6i zu2VXl51i&*a0#kEqo+`EwhIM+i=IOEIo=U7%z}?uN44x@VPDiD^%G%+lm8O=9(S9) zKGl&J{29l>HE7J|oWz80EIZY!pLfZaN$p}TFb12x6ohFEp6lQ&u6O8J5Z_m_R}QKc z^~f=~k3{mmk@35SI^kP6#$XP9z@pSiCeBud64!Tf>eRqvEHIP#J8IxoBbKVU+?mUs zXNg4Ai$CXbUM729$7S!Y)t=wWk4kdjF#TdcTvug9ZazfI9^b=?>bol68QOWCzQ*l= zA%5{kovFV5i9G^0*rNJ=mVe87!55)JZyL4F9$P&2{;_JTY+S>3Pq-3oh%|*NF>i>L)mw$ zvXi_WZ7iVa75%iol>w_0tVYdGw~|ZUjC(Kz!w0t=B=`=z_&5KCq7!Eiqavrz#B%EM z@Tn*{8_T{P#w0@)83n^s+Mr42WT;+5f@!caH`y3pwZmpy&h@3Z{udMl`b+_umPEcT zq)$OO+Yk!vz_JTI49bZaSaxoNE>1FP9%us}4A}=3bnv}>*dQzo9qG%5Jl{8@R~z{# zmfI0sye*D5f@pljfbVAI`ycnbx4%Gj4KTjS-uf`|d)>lXqRj^1%4dd}|C2dNGzs1c z=5L*Y43lz!dFnjIb|^4+JVDEq8)dUvKGbMo`a|m6p{CwuJb|LxV|a)qbq$61M;M=_ zU$MwZ_&}rP-EQ0_v(GfjC?o%rqc)5-mPGyHha%z_XJlWB$i{I-_E|~wk2C5%W05oW za)Z(DSfu_wF3}H}YWf6F=$Dc6Ckf;|lLT^Ll0e4KP*)z{SY_>ieDPO$>9e4}FgWh=xM3hb$c0=uM8V4qxJyv$yQ&%pP-(jv|kChdhVT$9%byn;0X zuX>HZ>sZI}=$B2h*TwNpx{NP`nWd~v0&nvsf!DA}-~kNhN1dw4&G9-mlvfxNLsAp2{&AfGN|5jox22api{secXovJ^w>M*SY706C4W-RPbl62K$cu=dh8L zeo8>EKP8~+P6?A-Kf@>QXc;)ek-xGb)boEe?hoT!?tclq#(#0V!xsb`O97v9YA8vc z>4~~N3rn=+3jy5o5&*Lla7jx6?VMRpm zr$r!-whH7tD@Rr#SzYEnH49J|#4=R1EH*M+_DdZg09OrQfHCU+0g*hNq6zC0A_p<= z^5;TMt+ng9CJhS9wd7_2ef(y@we#jkhxE&_o<3AGUq4hduNx|ihSA+Lj1S#0ERtms zaKCJI?1;$6!ujAm0<`fS4th9QU@{*}iNviLT{@bGG@eU15w7uJm1~(GfP1C_F!RAw zVYjk&F#pUn&K9i5*aB#_`r-Xy+1ecfG}|Fe;C2YQ;E0gmGC3&H-8_xc1!v4wW2Q%L z|K+Sfd@}^{r5OU*o(;&l#d6Fn0ZrKtd7MkkB_9uKiOw7WJ$jCS&YQzY(|o5wzT3y# za{D4~3Gm2l^^3)kkA+$2QVRs=ssaWIZdoAA%iMB#SGQcD=bQ9&SiU8z1oZJ$f^X+4 zVPy0(%Q}uVa$RWVW~fEQk=`&VBzBVkow`YYmTVHB7~@;ZI4C4g%^hZ_N_kkO&X)w} z%P$Gg=9dI0F7g)%2gQBgnJT>v99tr@!h}Wjs{;ALs{%RpH9%(W;JXBIKD?TPWmI#% z6H-js>Xz5TPFirEKppmmKy|($_`>c4{;28WpvXRVzn*Vke`J_NkGKvCmAQ&>aT zP72V6P6|-hNdb!Mz2`Is#W-GnMo0JX*|77LbzT5(J}-b9&I@3)>A4FW7`QVxgZ)~z z`s;s(%?+)s0(9^f0yOIj0g5jH3)%(ScD?VmTn@`L>8b#oe^r2buL@99`oND|X|{9w zNm!Zwx!eA8*cmG_nE2v)z+mEwYpcO@zxD#qM1C-uSWbtXn=aJyOr~(_s^2U?e=rNs zq-YL`NvzIl+Q;7k^8WzgF{ZDHE#R0ReSJ;geAf{#pf|<~Xn(xMm%de<8OX8FF2@r% zdqhsPhb3J#L?E9YB9N7vHDtP;^&iTyCJfbEq9oCDmK}A0moH`^9J+rBYpz15hubj( z?-QdAn;=GYPJmH$OUbfHLLs`8!|`ZJB5yxo zvXVP8OzYLExuyXuIbjmqut>-G@ZzwPvzCe>HZK)JG%OXAdlZ`Zt1hPNrsbS-$MUe8 zy(occ_$>0*?8{n~Q% z$={jYHP8~NnGrROtbPF|3#4hp{|8ehxe^jzhi;=bnpVr$t>#7O#u}#n8yeP^fie_j zi4#y#+)PquyG<3W!}lIB&h;v2VX5TpM^>7aGaTD(a-4dsQj9R+WznJLZ>DJ^bDJRM zD}nKvKtv)P+f2`~6F+LvYbJ5bCiNX?Q@R_v(6c}!bu*l`k*hV4GF~l2)mv>^$UpDc z_Gu8!Z*YrMZh_8Q-xQq#HV}K&aM2T#VpT_tX$=QxKA-{E>!IPrgId2o!=N({3DI>( zOo(@x6sN`<2EExQCihz!1e_MiB)g7?KK{jk82^qCK}j~Jp8PFw>Zx}$1czTkz*9!K z$3+*H2S}XLM9?vlWY+Vip43O<)be`MaDHsb|A7YM`cSZS^#cx*8wDULS+9w{SKxD* zxydw{>yq^s4Z=SOIxc7y9UZHwbu=-ll3M4Kc-p)UaCV&1aDrzw9N*sr6?BBve;0jx zGXU|5COX|FaM%`8J=do5V+{h`!E#=7L}%F$5~;IX6QXJlm|o%#b^q2N(2z}^35eil zoD@3Ckhp;gnv(H26WhRnUJwWQv1P7?#vTQa{CPg7S)gVZfwYB;T5X*g)e zYga@|`!~Wo?u~%h_KkSp+79-(?OV+r!K)es8dCmIv_wO$)^7GuB(trzLKMa!{)Llb>M3+g@Y@#1{OC~v@i4eS;qPHnBPF)%kCEO#quLglWJvmlD;8fIWA_Okz z7RR5e&bX+TxU+W-&>%2yEeQew19yWaf}y~}gy`|kiBregHD&6CXfQa1O}7XbrrWJT zWV#7a4cw|RDH^6h;Ov!-(E5xJ2EZU#eY<#o;qun)8jIvn8bp3C2@Yh9iMq)0Q|B1% z2}XhCZh?_>w*W_luZah${T_{4-8ii`x<%7?(K|Sv_pY1@<~edNuZIj9r(V5RBUUs~ z13;%JO%ni^L*)B~h{jtjMAdu0#vwUfgFvE_9}<0VX6iMO`SE1&6j!K~lQo}7bqE}w z%>63=D{zmdsiF(IaayJj+2WfilsWvE03bM-y=I7p7?PzSaR$f796d2x0I2?18Ujw} znb`sXm*vS$A;Lny7YpyMIMtA&A-LuW1k7EYCqxfi=qyhP5te{VvOyE`n<0pgJ*nBE zWTAi|7%*EGi;iGoVH@iArAzn&)we{Ok@RI+Z*T)x36_iAuH{;9YT3yv_=C+#och8F zO$%j}27qo5zeWI{8+bGkW5}{ra7bFKacFR95SRhatrrj&Lzx>ik*;J-Lc|?_@@&-m ze357_(O_U`t~ri6vs579WS`MQwgcFtDdR5FAaLki<)R}lRq>lOakH3I+yQLX5FC{n z0>)@=mFR*?{y8DSlyioGvyf8Dw`d5qS2P5Sks;d!0>;QzO~e@K2#BoetM% zovRLLozZN24~othDMJovBHa#d)kObgSQ-9tNV7!NTN(};Z^04X5#!~ekO;=h%y;-h z)%8x)$J`~H$20_3TZF;s3g1>QNZ_)P`#u-Zsl4w;iEX95K|{dV7}F^F;cV>H#O6kU zj|*H*Q`A2SVh`+W(Z;fWB+$@-+ddM#(SdW%39&gAOt#~k zHtBf*4FVl_^##!qr}2s=!Zc1J3or7?BesVisl0wsGkxOc8V-hJ$`=9$VQV!JCwTZJ z{#12b(!8McvIc>B;cFcN!X6aL;k32|`2)0l;#ZpOtH03@Fb5v_P9Oxn69|}d_I@wM zZ~i{YC7GK?v-6s8aN>QVU3s4Ujos-v1?O_Beu(NVvGvsVqlS-Wt^Zlz<06o9U5L)> zVmPXO#c+jhuV{oXScv7 zj4(ZvkwjuvD1FJuLZt_JbD82N0}GU-@}o9#-!|nol2D)|t0&EhS1$LDBB>RospO+> z%3vI|vbz|isJj@2u%+52n6RZfMu_MYtA(g~W3+zGUL1^|F^2aRUC`>J zwU3}s*GGd&sxZwVANAub3A#*Xtk#9@`kOQ{7%LF+;xq(IrK{sb7jHaogNb0l0A8=U z2ZVaiMRV6c(F0$e;uD038%>W8)w+aGKf1|m8YCLRc5*72G+3iP5V#&0EXE5A*7|vG z7Usgdu;3PffLZ#Y5LM+?tqZQV@reQ~HBtE{l6enKQt2H6BG@hsDjvf;q=cXIak-+g zRK;cciLjwnb+)DlU8fh%$yq#i z_R~)-p1b5}Y@m~8FUc=_Dt8&1jfFXdAk3Scv-rt{^Po=1xl0x=E6iD3$eU^Pg)8&t z(VmNEFI%~2(Y(Tfr)aN;h=L`B`kI~*E9Na+lKWI4Jf6LPHZYK5&rAJ?^Hs%O+`NFz zU@Lt9|6+f`k!R*?v+oOPjH|`w+y_bSyS^hcFBL7RnEh)G}qY z6s;R}uJVlBT8QnO^PpYla@PGeKy4_%vX9ocETd(JO|vuwN)@DPsUa3GD$0t5;nDMA z=kt_O*`^A05>sl1bZ`;qic|`zV(2fT-&|ETigA?{g9Jka&>-nIR=B%Xv`adP2w?Y~>$XMSXuZ6N(q%~ev;Rn~M_&)OtOFAd&c16=T#d+~QgJt(j z8awM~x$+gOSt`a_yciwvJPnJr#GJV(yD z%ADu-z>JQI8#y`nE!MhBvYAZW2#Lh;h9W!H!t)B^U18`+hNYTIWP6r{ymCC6p#lDs z^7j-BKGpq$b&o-^D`fnSRyKS`z*6oB?CUT4#ztc&_GjF{TBFnBv=;r&J;|)PCq=u9 zNzwb!tD>FHu_Wss?-`#Ei*`l@?8l=+3b2}PmM%3A{^QK6 zQgjxP@@u=8Xi6__YZcm_;?`>akz~X3jnRi2E$w|N9~Cr;#FNpp2qfYqD!j90Z8J}}FJTE@!3kz-)QJ8v(3-l)<;71MQ zme=u)fw#~s?Z%Sh4YO2@r6AP32BmjBuHva=Pdn?$ZXX`G%twXGIoiQ}IN;u)8|zY# z32Kk8^lhTAa3FRy@f$(^8w>oUlrGZv{^7nCa&msBWg)BTw0z47^u}-f53CQtD>`Sc zaf#jslGER&+==%AgIPL)#iS)xGx_RA%TPJ_Ev#^UYBEdjK(T`4{{Z&`hX2d}&_Ls^ zbJB3NLZ?`^!ODyXmp97C6IR7l%k&CmQLN`FNNiE2b8aMSHNQ+LjbQX)wf!)l~-?50sb$kTn z3iRZj1Gvg~&1}d`oEBa)8#7R+$fN|`5bgo(j1FBPWrX<^P);~tA$#rC+wn^MObH6~ zD!pZp^^ieMgj@FL8;L`NZ!rEBDLv{OhE&f|TLMA)=QaACoQwoQ$r~J=oDm7cIY4eA zy@qN8+sTkbO=u&dZ{uw7a!J}Su2%6{>(1!V^-FY&*r0`+d);t{bm{-V`U~zM1$S^< zB;DGM{DL_qGg;aOLFJ-XU~eM+GVq8tETX-Y4k$BwHzqLv6~eU}#s)^(rO_N0FKS&I zqvspCUQ79!o3Dx+Smpl!=lCW(>%iBX^{`|;YwvfI8JgPXnV1V^U2ipemTqFJk^9>a zZ+|D{O_Ii0*-tg#9c537&dxi>Gh|faz`cAJ3_#trRODU4Z<#thJ7R{)5kH6LUG*d9 zePN=WcUA;$5uIuzES^7E5gQx9McRN4n<*uYSz>nn7}y)JqRpzl@Y!-6==3t(y#r zftykGDb^%8zZ=`#!nzMRhT1gvfdqf>pUUSMubeH*)ERdDG)|jdIomv46Sk5)opL{g z|Fi(w((LN&n zv*_D1H8yRe@weP`cqQ`OEN)-OCf!)(M_JEmdI%mQaNREK(H)(GpkTW++#Xz~>hDUit(l3hOQ>ia!!@E|1lN zKev;%rJDYte@R|$PPDRRrvGAYGn@qX2qa50pO^e= zf}B9y3&2fo<0>e!YdGon36SQVxh$LOlhXHtq;^Ez>Dsafg0u zf=vFC=?Q79m87hZ?@-44w_22p=31P7TCYX=-AEF*8uuV|zuhd2#>#vb=yMkqal4TW z#R}s0(e1`6#!iKwMd-#H2CofIUc*^q%9*)VuY>M)Fi7ch@)JrLU66*=d1g|<@d4adZ3y*Rsm&zWr?j#?UTAPqE zyQcagACKODd~})=tvY^h-Ny`x&d}y&4b&dx`ivLtVbs)-9$|O~0AVj_To0c#y1tFK$+1^Lui5ZjXJzHT;g1Qy zS@ccLLw_nRVm5o#8vwi<#aOo81m6^4YTxIG1 E0&=_uBLDyZ delta 22980 zcmd6Pdw3K@)_1ylGU>Sh0|-pGWdMZ;0^#Bs0VRO0Sx_?)&;TM6ZblIpAW9@iCa`K& za5G^b6+o9E=n}bQAV^1%WhA;pq7uZ{%?IupVbw@b22l*6@|`+W)zg!W&%6J8&-3vI zsZ(9&SHC*v)Mcuxx_Z-y@IOcTla=t0cKZf>lOU%QRdrV%RHYE<__-}9;CHA>@SLq= zUoTSM*`!v59k-a2K;NROi#KePl`cs|fyE;etAgDQB+y*ap;&%RsklVT2V5>yss0X= z!CW=qb|oui5kyLPD8W*-&1I7isj^QyOpuJviw+n5 zf+nkcx4h3tWDJ`sW8um^Tgd4TOx*(0Mt0_`o}quME{$L4H?pP;$hckdlWb}SUxq3bk#Ws7QWAErce<3+A8b{1 z9lOYZipd|-3XpxkHKj@QU9rW{e}Qq?&VYNTvnp-sOI9pWJkvCR9KRy{DX`1ahvrnB z@VqT44d2_ypovodgs7q`&^X(nc&;Kd#~fAheNL0a)k#&QG*b9IQ1Z-4Q6f!{y%G3S zrO*9iJT{ctZPDZNEkX$NUTLlJpy-0<7;^D*tUjtEJDRJe7mZTT zL1gVdrSU7;D4Rz@d&SjW2617)r8uiv*3~8j8h5)$!_R8>fOSvTK;sQ|6Ig6Ym21xw zteAXOV0u3=A%yji>%}09*vI&RuSv=svRW0~w^B7h7)%Mc-n_utVKz^%dhF0cN}%y? z0{w^){YZ-xNPU|>n*DayK-j}7F4^a*Dtp_$3SqW3$G=aOrAe&vf^XudX^}MP786#c zx$c(Uk}ydce=p^-l-HAxNp{^M^=KrqnX}DUXL1EEMTTa}OSNr#j1bN$0g+TK9Tr5hwX`GB4U$f~Hf#heIwo|0 z(^O%b>_Fonr{<}bZV)+3w5C(sT7iouUN@ z)OGP^V-p)>nKhk!))gA{fGV$TkI*3M%~+y2dh=okUDS_DT5+;$@KN3nc`n+KgN>Bu zh;a<~+%12Kjzx&Pi2|Y^!&p)1L1RV94~tS@azjPVM|nO0YcJJGo{$aRRF?{+4$F;7 zEi2`SNr@Y!?%Q{(M0_nuqI{1^)|%uY&PghUG5bF*u6`m*m1)CT5-v)_TPJmtpOpin zyFY0uyN^BdkUpc;B*nGY9_Z1%my-S*d+hy8JpQ~??iO%#agmlK35lfR^YTCB@LKKs z^KxIqtP}hb?J79ygkGS@Qoo#{_=_;ty+_ru{IVV6hm7s1Tp*e2W!8H$%9L`F7?2mR ze(JfRxJX4nPE%a{?4&s$_mXR86YB=Kw-PehwLu%?jgp)`Njthx{*Dbq8Jc{7Y~RG- z=%`8dZ={Wj<^kIS(i^$~h@Ih>``{nWPnO zmko2ImiDSV2K~=4l6MM>S&4S2Za2pm`Sf0UaQ;pQ>fqPGg*0zg5|(z4orWi}g@i+4x#+1Uh$+ z%=ZjDW*_G~s72E14Jh#Wy!g2bWZ6lB#_EqaSkguS3!LIjq)P@Z>tlJYp$$0= zybZ`)KEoBq_66H?{*)_X2=x1>5vJ*%97Y4w9}M=5d@1akaRDs4_G<&L&F2}s2JSPu z$a_I{D+Q{@j=h0TaB(SIcF$({4Kf?^)dj7}aoN3)8+O19t~?+s-SzMJKMvFR~$M)`#F+%P+#YL!W|hzs^lZ z*AW}56zoc9@S~D~u}^1wEhsC@x1_i)VBBt2Sbz|WiBqaU@2rLEQ4-%-nSiS*t=gtA z2Wv?k6n(|lOlEiDr7$R#b!KP;d%>C;xP&2+-bGnw2oYbR!h0P?MUS4WtOik`NB0FX zy0-!1U=l~rpd|(lAKk<`U{^oaU#QSXxjaBfG)$|LV>m=^5gaOSg~(hxlyj&Y$~h2l z(gm_>m=VJ}T(E!+$5@afZ$ecQZfC{N1(Qc{Q&LaO8O;H}pXJ)kW0VZ#v5Ky-aa7Oc zh?C`~^LUHcW6L?)F*&sZ`H4uj@@JeyZmoGl2^n0{dn&#s)Kyr&d+ zT&--z<%+&6zl__q{?6EPEnIevWBEnAob%S$a%o)l+{W^&c{%rkD36X4gHj&Lh%3ut z1uJ;LqV!n65SQyGlHobZcR1`*W7wh(>|{$v(+;xYac(&rb}Yt+`N|77Gad7Xj7>=!u) zoa&eT+1T{d;+XLA1_TKh=y+{7$n_vx3!`^X8Iv?v(HMg~`DL!o7CtpwSuwB#EDwUZ z0Lz1IoD#-(MR^Q~#=i;@Fw$`nY(Cs+fY`i?+vUhEg{?B#8l`cc-~i*>1(JV&Q{8&N zXw}UJjRZWJ80;>-X&iEg|BaOj(BW;u5r#zQL_{*fT*4MO-;Qa4`(4ffW8nNT&J>0i zgJZ47IfV}0P379j1VhijxY;8 zWr|wyr@DRVypi@VVTKc*i!mRUoSnZgkm&m*$HIAN#8;fe)=FEscJ-^6piH?Se1Uzi z?pt0K`{0?2oW=EvMiwM=NeL?cpKaQqx0Nn3$@yNvZ4Yhi4@yT#@!Rbr^NO-?vhO9{ zUD?7WuCc!fQx#k;t66;<*6MAL!L zy{Z-IX}9_jev9Fxt|U8`R*6?2zwrzUjdpI~CJP zwsgWNXidB+y2M#;HnAlR9oEYen1K5pq)O%G0F<;pX0j5W#gwdgHq)OWHZiZdH&Ky9 zEGAW{eHqe+1L1on6RXXD+uhCtZLHm-$LAJuzMW|&4!y9|wv9yIf=vSN7DV%O0JNs9 zEWKP+m4?4yI@H0GqO=5Qy3_$h%KVUA81OukN$MPot#;92dXl%6!Zxk|^_EhqVoT=J zLDNP=bM&Lg8!+vI9YLiw!1BG2AGaIRuv5r(TQIJ7bzzmT(dSTCUgc+9qs?k2r=8r2 z7}?E?rql5qDxsJb`V>xj2%{+^}{>?a`h{?0csHSFJg`*UMb zPkhN~>-lsdxz^bPLc zkzg}GPMXE(AxP~k3X(QQfIL3O2;!T=L14FI9Er4#S2fps(<>_L1=pc@fyG8Nhu5T< zNNTP?o1JS!%gY5c#>8<7m~2_jb&o7JeM2E>Em|bcl&+Ahl>%e+N+U+uN*%-XXTUgA zz%gJ;GK%3SG^N`}-a3I%zRrkIy-t{h?qiZO#T){Bj}QbclU`y9nuuB|Fou^JG15x~ z2Ks!#CXN9J#Pd` zJuiSzb^15Ds*z0Xx2B+-6kQV-`>z=>YOlp$od1<$fSsZ-^m65BSx?%$ABIL9;!)-@58!#$}bRwRW2i&ls*)x;zowlgPgiCzKi~*L2J4(-UBP0 z0ivGu1B~@l3=s98@sAIT$29@DIMB#7^_F-V3AYK1)3+HhBDV<)773#{ITi_HI3gO( zH-=Je87qMLrWt`!(*zKUgb6}*f{{$^U6f3kTVOoyHe&ePF&O2U90REZrWoz!m>LhA zt$vyS`e~XG$T1y&41@iMnZ{=IeME4hZ%u-aU_8tgNaJT4k+NoUN^DS`tH;AaW3!qT zVLU7q7%wa}Vg#3pX0dp9GBzHb;@pNNo2P5Tmd8H?@vuT5eYC=e)U-k%p^>etcteJ+ zGP21lh==Z$x=vusUT4I}TPHBEb8IQ*)gLJ~VyGpwMZOIJW9J4VMtFn3z^UUL;TX7V zI90|aj2ibQ>haoF1kiV{7=hGR0f;#&M4)pFs^GdA6-G8?uTeIRJpyFZ9wUf*j}C&9 z;j$11L2`5U8R<0cqj3{FD3FdHG$J(~V34Fvc(DkcN8DjY0$JEAo zp9!4u&x|Ie^sa>`Tak0@+g)6VrYuy`jy(GOfvI<{-|U&4D=0>nXgx2e&;;s)P6qR zjANhH!ED~AkVcD8w_1(rE-P2Z-dt@r?`Hlb-oM7zkr7?ZY}UZQ4RtkxLA*T#(&ipU zq)-n*1nqFDH&?|L9Zh0v*a|AmWrh}2_Y**8`x${+`UxOd3GfAKK5!HVhf&cmC6$e)wK&F#T1JgE*5V!uwHOA(W#f%}XHPKlZJc0+;Unl4NXOkq zq(--ngwvTa#cZ72E`x|%(T3islvGoefa)~O2$eidG#wi@@nOzl@x#W3RnMdtDcJ&J zdbSb6o6Rw3k0o!;H`~dO4D)KOY_7Q{_brYIYZe-K9b7~q^Ogw6@+C&d>Lq#$F)YsH za@Fg(MqcSpQIPP{0_6PD29Sd*Vw%?X8QwH(Xs1!tl0 z3f2@0Jv<69+AqZLccAyQU+Dc95PGGXpgLPA7yzmxLrcw1vU(jiD~a0Mr9xrsW}!er zpFl;KFA3yO5r}mlVnZ&zWPX+%htTp~HH#ApX>UQ9wL64L+Ipx{snA8rJYR@fu)_QV z|CQ_5twYrB;if~vTlNa&Mh9SZt`wqY8np9)N}$btC)e-S5nMG;^4tMYvhf3`b;j#L z^u5mU0}YyIszQD$Ryhj36=070POLW5TTLGxdHE{MD65TItO=HM?kO4 zIVx0q0pM^(7om+{*&6;dkamBfR#Iag#7}&LkLfVJ_e31|x&ent#|0p^u|^j|FH3zi zOPzTVZ%fvPIz)IpC@!cMitc=B91ukqOd&H)@`sx5r0$W($2tPU0S{&jf&706ER5jG zx`?4StifEv+tKo=4uOGS`%EZen5>D4G)%7PqU$|y?d^@`e2!E7Zyg8iSofvCiG0bK zVXPebO0=W#E4>}==BP3?;hb-TGMey=E@GHue#@W2DJZS%Tl2l#o=ul@6g1$q%R0*U z!Z`jH0KWNq!NKu^`BUD2$W>h%EvWn~w9$h3y2!?WsEA|0FCxBP4<(3x=9Lq;OVav5 zt)rC$K8gfi0kKLK{j20e?UIrpdWkndhe)vq2#k2??%_E~j=l=eCT z404%dfQR2t5^=mMEUGku6bY z($Z7M!TxGX(o53*x<(hlSP))9)LgEF4cx8O{d5TQUfoRsg0-!`5Lw%V2qQ`NzFf#Ul{K5Qz796aTNFAnYdrdMNNg0uFj*V$8BlIU2__p5(+?3x5 zTpXLO=?`3=n4s_6rcezV;~hG)ym7iR2F2P535H=q zxl2Hy?eam4`umN$bPB2IItr4VG*RfGV`_Afh4Lgp#4{;Dbf7f1fB>>Q@V%EHI_f41 zEDYianL=dKZl)k{@F4-fX0^`KL_zG3Yop>cj*o-#xak5w3(wRM(3PiV2?U&vCwYVj zBR`)Vg4+`{DO*SI%@qhZS_S3{4V>9*j|vgSe289MjQHbWetuN%K1ELm7=ry}%Oar& z))ls)5x;gZf1rgH>n=@SrfWk0z`XdR(Dps4Yg5ZkdWt{5X|m|41TjsjD*z!ndGuH* z05A#yx`;i<_Ke_=@{G=*)~9p8@$bxP0fD_IbB!+2sjE&Gamu&+QMW{Wk&Xj(^QL39 zGu8?S?6{|d$fo>tdRzR(Is#U`y+r8Y4Ao=3F0L0&#VLQi?o@ZFj(}ZsPMJ`_IsS|g z!R4GF&Tb`}bOguCIs$f&ep>|sc8@K(h~1-ekUtFsXT$m@s2juis(``%JZOhNaqkc{ zV;AXB!5?U86}m?nb_r#IUToeYv>khecIeB5Ddf|=2`;v5^c;mrII2ULe#>q>nyU~&>_%x3l4D|>@VkpM6kcicuN4%K~vlGmTn2p+d2%) zFWlUH*Bg1c&DP-4=Lb?9O1fOh4U&0S|21>do zDUQ*qmj6@2FquuDp`Ue7G-A!K0u+b+(bt9Oxh@#eCT|?1^5qom(;K&Qsl?;l*r;w` z@sCkmCy9a)NoA3y<;bddVV7)DYn1}_vbuh3($Ssp@AUurclD1`rcRyu`dIN-=m}Zi ziN7NMH~tprKtkca{rg)Lx6T)``0E^e*OY60u8w7Y$*&1&|8~B97Oa-OiBBO*Eb1L4 zM;v6)52{6*Wl>ucvcsl2NyczfCOO?s?Ta=pZ7*z^*Iw9^u!V?2H7-Ot3NK*D=L=B_ zcGT58ojDdk&kwpms9?lz7b1?GJ-Ud8XMlA|+Hlisa=IHQNw8C7I&~HLTc%DI zBTj*klc*!$Op)J1s04fPGC0mH=*jcpg>+Pdt_imH78>|tut$;*aTy&DqE?+0)uT)3 zx;~;Htdl2`@qKmby@BhYzM{UyzPeuUMqw@-@fO@95O5MYCqzx{udCp4qDQg-OBrPL z#S`~%(`3?Vh}t*0U|2guV8N^IJ_RFK_w}=&jH}(jK+Bdp6(z2Md_6)dMm zBh&^VG(A=M{ttfp9LrBZ&w4?au$>BQ`Ppw^Q@e!H^8w;wHbH?Y~3=*(kgq7 zoTHN6V^kOJ5{^A!3y|rf)hyZ9W^P%VxuMwHLb$4=+Jy5Rd+=OEXryY9$Ea|uf}YY* z{CChKiM%UY0N^TBDy@d3 z5IDP0_InCp$>M(Xl5r!|B;xqOY?1F`;1?8&d^aZX7ESO@S`i6V*<2(tMt!1H<6W^0 z;wDLLGuIuPTNt2*jyhH*`S(jxEoBYNp~>>Ob-4OujMBC$R$6ynv^K70BTL$v#z^aax)s5$jYnPjPQi19t?I zm~|$rOO*PJ(8!s@|Cr=ZQoYz<2R#=Q*RdUqT0Ae5)fv-xUMS}Vnhu2Feq&BH77x+l zuDO`5Zx1H)&xiEgr14>Og_3TK)|Zi~E+=>149CJ!-?F6~!`8TSB$u=hbwttg~0 zwq()T)l|e@Yq7%pb{BCKBsY>SJ=OL7SE2ms-~lo3Vt{4~NbG z!zheP_jQqF85^s{Iq^`kB01X`NY~C&|3VhcRmYJI?I7k(*x|pibJd|_$R+g;$fj3IS#T8@$e}aq4H-iZ+Aqk=l-QgOe*<}hMUQ@r z#&6$}tV)v!dE*g6kA};g6EjSr9X=Y5Z%2YE_(@X%+gBLkKf!H9T_ILNJX-4x36@KDAC>7i35;Co5^i zkk6?&M#fNSRl5%S+h&)3^ht&`Xj&(>Rv>J%GUZiQ8n3WJpe4G_hTK;QaZw2rN_V&zIttb8OpRz`m(-xdX6z1fg~z1bi`?|r%?Jb?T- z4?4(+@1f^hzop(Mj%QQO0SDj0KrM_Z>J`59<+D=X90!~#rn+@ zfAF8QgobwiT6Mh4{8bh^wk=^(UPC1d=_WqII5j4VDK9LRw0z1-|2umwro5Eiu~S}{ zFPSSiUi!e;DbIPF=~u8xVPIT!V(gR`KFJE?aDX_zLT%B;4KWXp*&6F!STLT1P{gTD zwIz$Pzp`VoKP6b?pJ~$O=1&^0p%~<%<2oxSt7@M}3agOMQgP|gPL8GXh;%C~jXkfh zVnux5D&tVz!dP{E!Sozm8V`vsWaOro+FGy1z1&u=v(4O;HgnTsa|=o1+wuJ*mWh)k z-1j!n@U)Amu#lwH$POj=H`>W@(G)s@NxEpNI7E}gP1jP@Z$bL70jRpo_SDDb7M6s) z1|3gpuDDcj?Zj$usgim#j)J&U@ejf@&LhEFFrD&=%_8r`;{GWxARWXcu2k}>A$_-@ zvb@%nxxKJ1N`8l`*u;V*^<6~31&ewVABhEvvxbVeV1Zk=`ho?2X82>V`AWN@&D?OC zxeaaPPUCK4Q#P!3>G^)8>99&{g_5TvJX3GoBJTYAhJ9!rf+WPMf*D*j)I56|b_B{;pnFF;3HEeOQ`K(?u>! z<23C*hP-i_uINkiahi614A~2HOH7n2f-$BjB<=%}#TxCeY_3MPm9-*7T<=0iR5&TP zK3SJ_Lk>7u7j?xnPS)xFV0`Fg9U6@JI9Ufzv?3IWt%Vz|w$0qeHgj9r%yn#yX(-1{ zj?FC$cym?ao+86_iuaK)PU^l7Xj-4tQ+DAanABmhu^ZAx$EU?2@_Mq`%=NaJThwN5 zFgCXkzlkA7nETke-eHyd6#R<~sqFvgR)QFE|G(IlvE3)ll*I4d0!NJ%=Yt|XE7ae~ zb{w5T{(rlZG%)VJ*hxYvIE@bF(7zN?1^Mboyfpa zwa6NNiG8*vG$Al1$tfy(fD~O-Q#v~09Ly(2na4cjy`R;at?pG!WIF7GTT1`<9GJeI zkNl!)kzXwS0+L#)uCUgB#>~lO9kgk*ulYw&+ND3hs^nc{?L7;u|BgxA)d}AL z1)?ozPkhf$5g1{_(M1cwpb(a-`=y&ssIl zv$hLesWG&Smu&HP7l_@*&Mu>+D@ zx7pZTF0Hbh!{KG6@>f^^6!yQ$1-Q9Hzuu$deT(?dq`?JI{SsL-8MVt$-NfCT+~VEFGLP?N*#s|fc^`%T=0wrMyYxM7%}5j!C*YArDmB1}F`` zU`CL@Fsw?j+T4iu;$C>eKzH+-oG#q)Gc@VVw|E)&?0}*sz)j@Cw)^ml9nFVr?@Dr( zi5%}`pRYN?wgb$K$i`!#YuZ5LUt!$0;_}ur6*dp!F&EY^@Cz7ok&@2_X^6NQxQI+| z3eoF-o!;eU!b|M^{?3|;kA@%PkKn?0xixx^8ODHlmUMpQ-V@l?aexwbv?i=%b%3w^ zd>~hn*>xNS2|GI@A6Kt<3;Dc@Jt>g(M+?ce!zP0NJ?I~uqHI&UAlo!Eyr4Psfq~cY z4;i3#@k3!jdSSk~p1TR@MsNC4x=NhOPVx)lhsfRkgQz3XP8NJ@JBjW`j`ck$3p8qD zFHAy|N;za$SKv#H8txIL8R?BQF?g5#gHvFiaA@EX-fj;3A43z$Pjj2#1$4vb3mDl! zoa{|-XZ%y!V1<6k<7GLZ?3P0)xjA{HjrqLw9H%2os>HiL1^*o{z?%qu&(o6=jYe{4JDCCHu5h-piaFzy~K#9Z1fU0%$EFA`ehWsuxRe1q)K3GwVL+yVF@ zg3E(BG_C}cp`s$%jOa&>zz6f?RWLrxl!vH!|1JOI<-C6mAHl1{{@Jtw`QZX1A6e5V z@-4^Qa50R1Ph*m;GM<4nUuE3N9Rc4#XvdZa_^pavw{g-q-W*EhQo%)Axi$LzB{C>t zn=d~Pab02>qOSfg+aNdk$fz;)kMt-~f`3HH)cm1m8PD&*(1n;w&F^1>xe=QC93(fA z$eqmm>`OqIX^bHD(b?lzoA6T^%f@r+IKD)`n*x(1ytE``6YTHfvW9+EBlB)v7BAM% z`Mp7tB;UgdvX3~OOXu~`d-nHFk*!+j5BA;6ei#Ci?qvnE#rJZ%BDKfw=jJZHUpF^? z2+1vKrXLE~@&H#p@_63mFo3qdxr~t|5ubCsC%w`E`nmC z%I9-?Q>6Wm>cTGen!UEw+hf>H*StTeS<((~j(w7r^OO3fWN}#7FC|A{bW$P)2K@Bc z&H!3s0r(z6V{Gmz#RYHhUhKnuF+wx870k|FG7EPy0ISmYG~0gZ3Fc4vGo~HA zAno?!Rt&yvjZ2Nb5BY2ryf87y7Ykx#c$@b$@gIm+r3$iSiptie{54$YL$if!TMREa z=u%r=9Ic_jfQWBgT(NJx;3Z@0j8HU#zCm4JR?-p3=BC8v!q)sl5C`fhtQBXVeOQx!25h#c z?8j_6UOHZDRcU%G5UJ^}SYH8)Bi5>olp2%q+QjX>6_k7vw;y*8l(j diff --git a/test/SrVO3_Sigma.h5 b/test/SrVO3_Sigma.h5 index ecfecdc3087e5710dc85769bdcbfee162f8ecd87..0b9d0bc7d2417d848de4726ecdaabd57397a49d8 100644 GIT binary patch delta 62 zcmX?bjqkuUz6}M;lP3v?Z9MRZak_vbqr~I}Rtfgt%=FyE`0~mBS2%8dz$(DPDA4?o Sv;89{BM>uf|H#SgzYqW$I2dRE delta 56 zcmX?bjqkuUz6}M;8#iz;F$zqsWR>S+WPkvVi5u-VKVTJLVdQB3$l3mplM#rSwtwVg H_Fo79`$rMR diff --git a/test/srvo3_transp.output.h5 b/test/srvo3_transp.output.h5 index 8aa6d1af3cacd32f84d8e342265bd9e2d37cc3a3..950baa0336e42e7b6a92540b3fa6c89a9b76ecbf 100644 GIT binary patch literal 13272 zcmeI2&ubGw6vt;1qoMw4QbBD!#2<%%AcFO18X+kW!NxYAw2PB7MGJ)c+3xcK?l{SvB|%x3X*)j|_cX6-$*hJ*yWtR;h&Z*`KfT11@sq-crKy zVV*Z5@LRTH_hT1x6I-k)|ABRE8tYi|7$M1ljJ=h=Q{3@CYB-MdN# z(sa6JlL6up)bAKgkF8Vh+vt+ic7mUY-QZr^nPMMlqpY+p8p}3J)(<=`) z^Chx5$aS&t`P#>KiQV-VZ@!(rwH<3b&on;W-Hr`nOP&z|LO=)z0U;m+#)AOGpT6#? zLOn$F<~U>bGTpb3A3?LCh+kOG9@ju^VuVxtnF!9WeJ;)|}<@`8rcq&^-CH5`K+OPpmpVNJ*MS0PNOjMYW8q@!z%6dD*~x$I%W`l(YvM$N2~&%E z#~lu2h6Lo_@%VvpCp;`q^J~Axl*iwihjW-xcxLGxA4@YY$6#IZYW^qyjPrOQKEs_q zhBR%-YlVOi5CTHLV*;1(9qk#uqrJs@X3Q-v*ui%ri3>dKo{=@~<{uI_j%?RS#e{$m5CVr1_yhU! B3|#;K delta 169 zcmcbSzC(3_2Gb6qg<9;Bt0V*zH5fpEk%