From aeaebb04aebf694f18f92007480a482d76c43bf6 Mon Sep 17 00:00:00 2001 From: sobeck Date: Tue, 15 Sep 2020 16:36:57 +0200 Subject: [PATCH] Revision of Wannier90Converter: 1) write the Hamiltonian in Bloch space for charge self-consistent calculations, and 2) spin-orbit coupling if bloch_basis = True: * if "seedname_u.dat" (and in case of disentanglement "seedname_u_dis.dat") present, write hopping in Bloch basis * "proj_mat" transforming from Bloch to orbital space * diagonal hoppings are directly read from "seedname.eig" * fermi weights and band_window of Wannier Hamiltonian are read from DFT output and "seedname.nnkp", written into new subgroup "dft_misc_input" * automatic calculation of "density_required" * implemented for Quantum Espresso (read from "seedname.nscf.out" if verbosity = 'high') and VASP (read from "OUTCAR"/"LOCPROJ") * spin-orbit coupling SO = 1 implemented * substitute k_mesh and bz_weights with kpts and kpt_weights, respectively (previous names kept for compatibility) * updated tests --- .../triqs_dft_tools/converters/wannier90.py | 535 ++++++++++++++---- test/python/w90_convert_hloc_diag.ref.h5 | Bin 128048 -> 240862 bytes test/python/w90_convert_wannier.ref.h5 | Bin 128048 -> 106377 bytes 3 files changed, 427 insertions(+), 108 deletions(-) diff --git a/python/triqs_dft_tools/converters/wannier90.py b/python/triqs_dft_tools/converters/wannier90.py index 7f4cf8ab..0cadcf91 100644 --- a/python/triqs_dft_tools/converters/wannier90.py +++ b/python/triqs_dft_tools/converters/wannier90.py @@ -24,14 +24,15 @@ # Wannier90 to HDF5 converter for the SumkDFT class of dfttools/TRIQS; # # written by Gabriele Sclauzero (Materials Theory, ETH Zurich), Dec 2015 -- Jan 2016, -# and updated by Maximilian Merkel (Materials Theory, ETH Zurich), Aug 2020, +# updated by Maximilian Merkel (Materials Theory, ETH Zurich), Aug 2020 -- Jan 2021, +# and by Sophie Beck (Materials Theory, ETH Zurich), Sep 2020 -- Jan 2021, # under the supervision of Claude Ederer (Materials Theory). # Partially based on previous work by K. Dymkovski and the DFT_tools/TRIQS team. # # Limitations of the current implementation: -# - the case with SO=1 is not considered at the moment # - the T rotation matrices are not used in this implementation -# - projectors for uncorrelated shells (proj_mat_all) cannot be set +# - in bloch_basis mode, the number of wannier functions per shell has to be equal +# to the dim of the shell # # Things to be improved/checked: # - the case with SP=1 might work, but was never tested (do we need to define @@ -46,11 +47,11 @@ import numpy -import math -from h5 import * -from .converter_tools import * -from itertools import product import os.path +from itertools import product + +from h5 import HDFArchive +from .converter_tools import ConverterTools import triqs.utility.mpi as mpi class Wannier90Converter(ConverterTools): @@ -59,8 +60,8 @@ class Wannier90Converter(ConverterTools): """ def __init__(self, seedname, hdf_filename=None, dft_subgrp='dft_input', - symmcorr_subgrp='dft_symmcorr_input', repacking=False, - rot_mat_type='hloc_diag'): + symmcorr_subgrp='dft_symmcorr_input', misc_subgrp='dft_misc_input', + repacking=False, rot_mat_type='hloc_diag', bloch_basis=False): """ Initialise the class. @@ -74,11 +75,15 @@ class Wannier90Converter(ConverterTools): Name of subgroup storing necessary DFT data symmcorr_subgrp : string, optional Name of subgroup storing correlated-shell symmetry data + misc_subgrp : string, optional + Name of subgroup storing miscellaneous DFT data. repacking : boolean, optional Does the hdf5 archive need to be repacked to save space? rot_mat_type : string, optional Type of rot_mat used Can be 'hloc_diag', 'wannier', 'none' + bloch_basis : boolean, optional + Should the Hamiltonian be written in Bloch rather than Wannier basis? """ self._name = "Wannier90Converter" @@ -93,11 +98,13 @@ class Wannier90Converter(ConverterTools): self.w90_seed = seedname self.dft_subgrp = dft_subgrp self.symmcorr_subgrp = symmcorr_subgrp + self.misc_subgrp = misc_subgrp self.fortran_to_replace = {'D': 'E'} # threshold below which matrix elements from wannier90 should be # considered equal self._w90zero = 2.e-6 self.rot_mat_type = rot_mat_type + self.bloch_basis = bloch_basis if self.rot_mat_type not in ('hloc_diag', 'wannier', 'none'): raise ValueError('Parameter rot_mat_type invalid, should be one of' + '"hloc_diag", "wannier", "none"') @@ -113,14 +120,14 @@ class Wannier90Converter(ConverterTools): - dft_subgrp - symmcorr_subgrp - in the hdf5 archive. + in the hdf5 archive. """ # Read and write only on the master node if not (mpi.is_master_node()): return - mpi.report("Reading input from %s..." % self.inp_file) + mpi.report("\nReading input from %s..." % self.inp_file) # R is a generator : each R.Next() will return the next number in the # file @@ -139,8 +146,10 @@ class Wannier90Converter(ConverterTools): else: # some default grid, if everything else fails... nki = [8, 8, 8] - # read the total number of electrons per cell + # read the total number of electrons per cell if not in bloch basis + # in bloch basis, this is later calculated from the partial occupations density_required = float(next(R)) + # we do not read shells, because we have no additional shells beyond correlated ones, # and the data will be copied from corr_shells into shells (see below) # number of corr. shells (e.g. Fe d, Ce f) in the unit cell, @@ -149,6 +158,7 @@ class Wannier90Converter(ConverterTools): # l, dim, SO flag, irep): corr_shells = [{name: int(val) for name, val in zip( corr_shell_entries, R)} for icrsh in range(n_corr_shells)] + try: self.fermi_energy = float(next(R)) except: @@ -169,17 +179,40 @@ class Wannier90Converter(ConverterTools): for ish in range(n_shells): shells.append({key: corr_shells[ish].get( key, None) for key in shell_entries}) - ### - SP = 0 # NO spin-polarised calculations for now - SO = 0 # NO spin-orbit calculation for now + + # Determine if any shell requires SO + if any([corr_shell['SO']==1 for corr_shell in corr_shells]): + SO = 1 + SP = 1 + print('Spin-orbit interaction turned on') + else: + SO = 0 + SP = 0 + # Only one block supported - either non-spin-polarized or spin-orbit coupled + assert SP == SO, 'Spin-polarized calculations not implemented' + charge_below = 0 # total charge below energy window NOT used for now energy_unit = 1.0 # should be understood as eV units - ### + # this is more general - n_spin = SP + 1 - SO + n_spin_blocs = SP + 1 - SO + assert n_spin_blocs > 0, 'Input error, if SO=1, SP must be 1.' + + if SO == 1: + for shell_list in [shells, corr_shells]: + for entry in shell_list: + entry['dim'] *= 2 + dim_corr_shells = sum([sh['dim'] for sh in corr_shells]) - mpi.report( - "Total number of WFs expected in the correlated shells: %d" % dim_corr_shells) + mpi.report('Total number of WFs expected in the correlated shells: {0:d}'.format(dim_corr_shells)) + + # build the k-point mesh, if its size was given on input (kmesh_mode >= 0), + # otherwise it is built according to the data in the hr file (see below) + # If output is in bloch_basis, we use k mesh from seedname_u.mat for consistency + if kmesh_mode >= 0 and not self.bloch_basis: + n_k, kpts, kpt_weights = self.kmesh_build(nki, kmesh_mode) + self.n_k = n_k + self.kpts = kpts # determine the number of inequivalent correlated shells and maps, # needed for further processing @@ -190,15 +223,6 @@ class Wannier90Converter(ConverterTools): shells_map = [inequiv_to_corr[corr_to_inequiv[ish]] for ish in range(n_corr_shells)] mpi.report("Mapping: " + format(shells_map)) - mpi.report("Subtracting %f eV from the Fermi level." % self.fermi_energy) - - # build the k-point mesh, if its size was given on input (kmesh_mode >= 0), - # otherwise it is built according to the data in the hr file (see - # below) - if kmesh_mode >= 0: - n_k, k_mesh, bz_weights = self.kmesh_build(nki, kmesh_mode) - self.n_k = n_k - self.k_mesh = k_mesh # not used in this version: reset to dummy values? n_reps = [1 for i in range(n_inequiv_shells)] @@ -207,33 +231,37 @@ class Wannier90Converter(ConverterTools): for ish in range(n_inequiv_shells): ll = 2 * corr_shells[inequiv_to_corr[ish]]['l'] + 1 lmax = ll * (corr_shells[inequiv_to_corr[ish]]['SO'] + 1) - T.append(numpy.zeros([lmax, lmax], numpy.complex_)) + T.append(numpy.zeros([lmax, lmax], dtype=complex)) spin_w90name = ['_up', '_down'] hamr_full = [] + umat_full = [] + udismat_full = [] + bandmat_full = [] # TODO: generalise to SP=1 (only partially done) rot_mat_time_inv = [0 for i in range(n_corr_shells)] # Second, let's read the file containing the Hamiltonian in WF basis # produced by Wannier90 - for isp in range(n_spin): + for isp in range(n_spin_blocs): # begin loop on isp # build filename according to wannier90 conventions - if SP == 1: + if SP == 1 and SO == 0: mpi.report( "Reading information for spin component n. %d" % isp) - hr_file = self.w90_seed + spin_w90name[isp] + '_hr.dat' + file_seed = self.w90_seed + spin_w90name[isp] else: - hr_file = self.w90_seed + '_hr.dat' + file_seed = self.w90_seed # now grab the data from the H(R) file mpi.report( - "The Hamiltonian in MLWF basis is extracted from %s ..." % hr_file) - nr, rvec, rdeg, nw, hamr = self.read_wannier90hr(hr_file) + "\nThe Hamiltonian in MLWF basis is extracted from %s files..." % file_seed) + (nr, rvec, rdeg, nw, hamr, u_mat, udis_mat, + band_mat, k_mesh_from_umat) = self.read_wannier90data(file_seed, kmesh_mode) # number of R vectors, their indices, their degeneracy, number of # WFs, H(R) - mpi.report("... done: %d R vectors, %d WFs found" % (nr, nw)) + mpi.report("\n... done: %d R vectors, %d WFs found" % (nr, nw)) if isp == 0: # set or check some quantities that must be the same for both @@ -241,7 +269,14 @@ class Wannier90Converter(ConverterTools): self.nrpt = nr # k-point grid: (if not defined before) - if kmesh_mode == -1: + if self.bloch_basis: + mpi.report('Reading k mesh from seedname_u.mat file') + kpts = k_mesh_from_umat + n_k = len(kpts) + kpt_weights = numpy.full(n_k, 1/n_k) + self.n_k = n_k + self.kpts = kpts + elif kmesh_mode == -1: # the size of the k-point mesh is determined from the # largest R vector nki = [2 * rvec[:, idir].max() + 1 for idir in range(3)] @@ -249,9 +284,10 @@ class Wannier90Converter(ConverterTools): # wannier90 convention: if we have nki k-points along the i-th direction, # then we should get 2*(nki/2)+nki%2 R points along that # direction - n_k, k_mesh, bz_weights = self.kmesh_build(nki) - self.n_k = n_k - self.k_mesh = k_mesh + n_k, kpts, kpt_weights = self.kmesh_build(nki) + self.n_k = n_k + self.kpts = kpts + # set the R vectors and their degeneracy self.rvec = rvec @@ -273,9 +309,10 @@ class Wannier90Converter(ConverterTools): # we assume spin up and spin down always have same total number # of WFs - n_orbitals = numpy.ones( - [self.n_k, n_spin], numpy.int) * self.nwfs - + # get second dimension of udis_mat which corresponds to number of bands in window + # n_bands_max corresponds to numpy.max(n_orbitals) + n_bands_max = udis_mat.shape[1] + n_orbitals = numpy.full([self.n_k, n_spin_blocs], n_bands_max) else: # consistency check between the _up and _down file contents if nr != self.nrpt: @@ -286,8 +323,9 @@ class Wannier90Converter(ConverterTools): "Different number of WFs for spin-up/spin-down!") hamr_full.append(hamr) - # FIXME: when do we actually need deepcopy()? - # hamr_full.append(deepcopy(hamr)) + umat_full.append(u_mat) + udismat_full.append(udis_mat) + bandmat_full.append(band_mat) for ir in range(nr): # checks if the Hamiltonian is real (it should, if @@ -322,39 +360,92 @@ class Wannier90Converter(ConverterTools): "Rotations for spin component n. %d do not match!" % isp) # end loop on isp - mpi.report("The k-point grid has dimensions: %d, %d, %d" % tuple(nki)) + # Reads misc input needed for CSC calculations + if self.bloch_basis: + if os.path.isfile(self.w90_seed + '.nscf.out'): + fermi_weight_file = self.w90_seed + '.nscf.out' + print('Reading DFT band occupations from Quantum Espresso output {}'.format(fermi_weight_file)) + elif os.path.isfile('LOCPROJ'): + fermi_weight_file = 'LOCPROJ' + print('Reading DFT band occupations from Vasp output {}'.format(fermi_weight_file)) + else: + raise IOError('seedname.nscf.out or LOCPROJ required in bloch_basis mode') + + f_weights, band_window, self.fermi_energy = self.convert_misc_input(fermi_weight_file, + self.w90_seed + '.nnkp', n_spin_blocs) + # Get density from k-point weighted average and sum over all spins and bands + density_required = numpy.sum(f_weights.T * kpt_weights) * (2 - SP) + print('Overwriting required density with DFT result {:.5f}'.format(density_required)) + print('and using the DFT Fermi energy {:.5f} eV\n'.format(self.fermi_energy)) + + if not self.bloch_basis: + mpi.report("The k-point grid has dimensions: %d, %d, %d" % tuple(nki)) + # if calculations are spin-polarized, then renormalize k-point weights - if SP == 1: - bz_weights = 0.5 * bz_weights + if SP == 1 and SO == 0: + kpt_weights *= 0.5 - # Third, compute the hoppings in reciprocal space - hopping = numpy.zeros([self.n_k, n_spin, numpy.max( - n_orbitals), numpy.max(n_orbitals)], numpy.complex_) - for isp in range(n_spin): - # make Fourier transform H(R) -> H(k) : it can be done one spin at - # a time - hamk = self.fourier_ham(self.nwfs, hamr_full[isp]) - # copy the H(k) in the right place of hoppings... is there a better - # way to do this?? - for ik in range(self.n_k): - #hopping[ik,isp,:,:] = deepcopy(hamk[ik][:,:])*energy_unit - hopping[ik, isp, :, :] = hamk[ik][:, :] * energy_unit - - # Then, initialise the projectors - k_dep_projection = 0 # we always have the same number of WFs at each k-point - proj_mat = numpy.zeros([self.n_k, n_spin, n_corr_shells, max( - [crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], numpy.complex_) + # Third, initialise the projectors + k_dep_projection = 0 # at the moment not really used, but might get important + proj_mat = numpy.zeros([self.n_k, n_spin_blocs, n_corr_shells, max( + [crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], dtype=complex) iorb = 0 - # Projectors simply consist in identity matrix blocks selecting those MLWFs that - # correspond to the specific correlated shell indexed by icrsh. + # Projectors are either identity matrix blocks to use with Wannier basis + # OR correspond to the overlap between Kohn-Sham and Wannier orbitals as + # P_{nu,alpha](k) = # NOTE: we assume that the correlated orbitals appear at the beginning of the H(R) # file and that the ordering of MLWFs matches the corr_shell info from # the input. - for icrsh in range(n_corr_shells): - norb = corr_shells[icrsh]['dim'] - proj_mat[:, :, icrsh, 0:norb, iorb:iorb + - norb] = numpy.identity(norb, numpy.complex_) - iorb += norb + for isp in range(n_spin_blocs): + # now combine udismat and umat + u_total = numpy.einsum('abc,acd->abd',udismat_full[isp],umat_full[isp]) + # transpose and write into proj_mat + u_temp = numpy.transpose(u_total.conj(),(0,2,1)) + for icrsh in range(n_corr_shells): + dim = corr_shells[icrsh]['dim'] + proj_mat[:, isp, icrsh, 0:dim, :] = u_temp[:,iorb:iorb+dim,:] + iorb += dim + + # Then, compute the hoppings in reciprocal space + hopping = numpy.zeros([self.n_k, n_spin_blocs, numpy.max(n_orbitals), numpy.max(n_orbitals)], dtype=complex) + for isp in range(n_spin_blocs): + # if bloch_basis is True, use Kohn-Sham eigenvalues as hamk + # this ensures that the calculation of the band-correlation energy + # is consistent with SumkDFT's calc_density_correction + if self.bloch_basis: + # diagonal Kohn-Sham bands + # TODO: test for system with self.nwfs > dim_corr_shells + hamk = [numpy.diag(bandmat_full[isp][ik]) for ik in range(self.n_k)] + + # Sanity check if the local Hamiltonian with the projector method + # corresponds to W90 result + wannier_ham = self.fourier_ham(hamr_full[isp]) + for ik in range(self.n_k): + proj_mat_flattened = proj_mat[ik, isp].reshape(self.nwfs, numpy.max(n_orbitals)) + downfolded_ham = proj_mat_flattened.dot(hamk[ik].dot(proj_mat_flattened.conj().T)) + + if not numpy.allclose(downfolded_ham, wannier_ham[ik], atol=self._w90zero, rtol=0): + mpi.report('WARNING: mismatch between downfolded Hamiltonian and ' + + f'Fourier transformed H(R). First occurred at kpt {ik}:') + + with numpy.printoptions(formatter={'complexfloat': '{:+.3f}'.format}): + mpi.report('Downfolded Hamiltonian, P H_eig P') + mpi.report(downfolded_ham) + mpi.report('\nWannier Hamiltonian, Fourier(H(r))') + mpi.report(wannier_ham[ik]) + break + # else for an isolated set of bands use fourier transform of H(R) + else: + # make Fourier transform H(R) -> H(k) : it can be done one spin at a time + hamk = self.fourier_ham(hamr_full[isp]) + # finally write hamk into hoppings + for ik in range(self.n_k): + hopping[ik, isp] = hamk[ik] - numpy.identity(numpy.max(n_orbitals)) * self.fermi_energy + hopping *= energy_unit + mpi.report("Subtracting {:.5f} eV from the Fermi level.".format(self.fermi_energy)) + + # bz_weights required by triqs h5 standard but soon to be replaced by kpt_weights + bz_weights = kpt_weights # Finally, save all required data into the HDF archive: # use_rotations is supposed to be an int = 0, 1, no bool @@ -367,18 +458,25 @@ class Wannier90Converter(ConverterTools): things_to_save = ['energy_unit', 'n_k', 'k_dep_projection', 'SP', 'SO', 'charge_below', 'density_required', 'symm_op', 'n_shells', 'shells', 'n_corr_shells', 'corr_shells', 'use_rotations', 'rot_mat', 'rot_mat_time_inv', 'n_reps', 'dim_reps', 'T', 'n_orbitals', 'proj_mat', 'bz_weights', 'hopping', - 'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr'] + 'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr', 'kpt_weights', 'kpts'] for it in things_to_save: ar[self.dft_subgrp][it] = locals()[it] - def read_wannier90hr(self, hr_filename="wannier_hr.dat"): + if self.bloch_basis: + # Store Fermi weights to 'dft_misc_input' + if not (self.misc_subgrp in ar): + ar.create_group(self.misc_subgrp) + ar[self.misc_subgrp]['dft_fermi_weights'] = f_weights + ar[self.misc_subgrp]['band_window'] = band_window+1 # Change to 1-based index + + def read_wannier90data(self, wannier_seed="wannier", kmesh_mode=0): """ Method for reading the seedname_hr.dat file produced by Wannier90 (http://wannier.org) Parameters ---------- - hr_filename : string - full name of the H(R) file produced by Wannier90 (usually seedname_hr.dat) + wannier_seed : string + seedname to read H(R) file produced by Wannier90 (usually seedname_hr.dat) Returns ------- @@ -392,21 +490,30 @@ class Wannier90Converter(ConverterTools): number of Wannier functions found h_of_r : list of numpy.array = Hamilonian matrix elements in the Wannier basis - + u_mat : numpy.array + U_mn^k = unitary matrix elements which mix the Kohn-Sham states + udis_mat : numpy.array + U^dis(k) = rectangular matrix for entangled bands + band_mat : numpy.array + \epsilon_nk = Kohn-Sham eigenvalues (in eV) needed for entangled bands + If not self.bloch_basis unused and therefore None + k_mesh : numpy.array + The k mesh read from the seedname_u.mat file to ensure consistency + If not self.bloch_basis unused and therefore None """ # Read only from the master node if not (mpi.is_master_node()): return + hr_filename = wannier_seed + '_hr.dat' try: - with open(hr_filename, "r") as hr_filedesc: + with open(hr_filename, 'r') as hr_filedesc: hr_data = hr_filedesc.readlines() - hr_filedesc.close() except IOError: mpi.report("The file %s could not be read!" % hr_filename) - mpi.report("Reading %s..." % hr_filename + hr_data[0]) + mpi.report('reading {:20}...{}'.format(hr_filename,hr_data[0].strip('\n'))) try: # reads number of Wannier functions per spin @@ -415,11 +522,79 @@ class Wannier90Converter(ConverterTools): except ValueError: mpi.report("Could not read number of WFs or R vectors") + k_mesh = None + + if not self.bloch_basis: + # For kmesh_mode == -1, size of automatic k mesh known when rvec_idx + # have been read. For kmesh_mode >= 0, kpts have been determined already + if kmesh_mode >= 0: + n_k = self.n_k + else: + # first, read u matrices from 'seedname_u.mat' + u_filename = wannier_seed + '_u.mat' + with open(u_filename,'r') as u_file: + u_data = u_file.readlines() + # reads number of kpoints and number of wannier functions + n_k, num_wf_u, _ = map(int, u_data[1].split()) + assert num_wf_u == num_wf, '#WFs must be identical for *_u.mat and *_hr.dat' + mpi.report('reading {:20}...{}'.format(u_filename,u_data[0].strip('\n'))) + + # Reads k mesh from all lines with 3 floats + k_mesh = numpy.loadtxt((line for line in u_data if line.count('.') == 3)) + assert k_mesh.shape == (n_k, 3) + # Reads u matrices from all lines with 2 floats + u_mat = numpy.loadtxt((line for line in u_data if line.count('.') == 2)) + assert u_mat.shape == (n_k*num_wf*num_wf, 2) + + mpi.report('Writing h5 archive in projector formalism: H(k) defined in KS Bloch basis') + + try: + # read 'seedname_u_dis.mat' + udis_filename = wannier_seed + '_u_dis.mat' + # if it exists the Kohn-Sham eigenvalues and the window are needed + band_filename = wannier_seed + '.eig' + wout_filename = wannier_seed + '.wout' + + with open(udis_filename,'r') as udis_file: + udis_data = udis_file.readlines() + disentangle = True + except IOError: + disentangle = False + mpi.report('WARNING: File {} missing.'.format(udis_filename)) + mpi.report('Assuming an isolated set of bands. Check if this is what you want!') + + # read Kohn-Sham eigenvalues from 'seedname.eig' + mpi.report('Reading {}'.format(band_filename)) + band_data = numpy.loadtxt(band_filename, usecols=2) + + if disentangle: + # reads number of kpoints, number of wannier functions and bands + num_k_udis, num_wf_udis, num_ks_bands = map(int, udis_data[1].split()) + + assert num_k_udis == n_k, '#k points must be identical for *.inp and *_u_dis.mat' + assert num_wf_udis == num_wf, '#WFs must be identical for *_u.mat and *_hr.dat' + + mpi.report('Found {:22}...{}, '.format(udis_filename,udis_data[0].strip('\n'))) + udis_data = numpy.loadtxt(udis_data, usecols=(0, 1), skiprows=2) + + # read disentanglement window from 'seedname.wout' + with open(wout_filename) as wout_file: + for line in wout_file: + if 'Outer:' in line: + content = line.split() + index = content.index('Outer:') + 1 + dis_window_min = float(content[index]) + dis_window_max = float(content[index+2]) + break + mpi.report('and {} for disentanglement energy window.'.format(wout_filename)) + else: + num_ks_bands = num_wf + # allocate arrays to save the R vector indexes and degeneracies and the # Hamiltonian rvec_idx = numpy.zeros((nrpt, 3), dtype=int) rvec_deg = numpy.zeros(nrpt, dtype=int) - h_of_r = [numpy.zeros((num_wf, num_wf), dtype=numpy.complex_) + h_of_r = [numpy.zeros((num_wf, num_wf), dtype=complex) for n in range(nrpt)] # variable currpos points to the current line in the file @@ -437,7 +612,7 @@ class Wannier90Converter(ConverterTools): ir += 1 # for each direct lattice vector R read the block of the # Hamiltonian H(R) - for ir, jj, ii in product(list(range(nrpt)), list(range(num_wf)), list(range(num_wf))): + for ir, jj, ii in product(range(nrpt), range(num_wf), range(num_wf)): # advance one line, split the line into tokens currpos += 1 cline = hr_data[currpos].split() @@ -445,8 +620,7 @@ class Wannier90Converter(ConverterTools): if int(cline[3]) != ii + 1 or int(cline[4]) != jj + 1: mpi.report( "Inconsistent indices at %s%s of R n. %s" % (ii, jj, ir)) - rcurr = numpy.array( - [int(cline[0]), int(cline[1]), int(cline[2])]) + rcurr = numpy.array([int(cline[0]), int(cline[1]), int(cline[2])]) if ii == 0 and jj == 0: rvec_idx[ir] = rcurr rprec = rcurr @@ -457,16 +631,57 @@ class Wannier90Converter(ConverterTools): "Inconsistent indices for R vector n. %s" % ir) # fill h_of_r with the matrix elements of the Hamiltonian - if not numpy.any(rcurr) and ii == jj: - h_of_r[ir][ii, jj] = complex(float(cline[5]) - self.fermi_energy, float(cline[6])) - else: - h_of_r[ir][ii, jj] = complex(float(cline[5]), float(cline[6])) + h_of_r[ir][ii, jj] = complex(float(cline[5]), float(cline[6])) except ValueError: mpi.report("Wrong data or structure in file %s" % hr_filename) + # first, get the input for u_mat + if self.bloch_basis: + u_mat = u_mat[:, 0] + 1j * u_mat[:, 1] + u_mat = u_mat.reshape((n_k, num_wf, num_wf)).transpose((0, 2, 1)) + else: + if kmesh_mode == -1: + n_k = numpy.prod([2 * rvec_idx[:, idir].max() + 1 for idir in range(3)]) + # Wannier basis; fill u_mat with identity + u_mat = numpy.zeros([n_k, num_wf, num_wf], dtype=complex) + for ik in range(n_k): + u_mat[ik,:,:] = numpy.identity(num_wf,dtype=complex) + + # now, check what is needed in the case of disentanglement: + # The file seedname_u_dis.mat contains only the bands inside the window + # and then fills the rest up with zeros. Therefore, we need to put the + # entries from udis_data in the correct position in udis_mat, i.e. + # shifting by the number of bands below dis_window_min + if self.bloch_basis: + # reshape band_data + band_mat = band_data.reshape(n_k, num_ks_bands) + else: + # Not used in wannier basis + band_mat = None + + if self.bloch_basis and disentangle: + # Determine which bands are inside the band window + inside_window = numpy.logical_and(band_mat >= dis_window_min, + band_mat <= dis_window_max) + n_inside_per_k = numpy.sum(inside_window, axis=1) + + # Reformats udis_data as complex, without header + udis_data = udis_data[:, 0] + 1j * udis_data[:, 1] + udis_data = udis_data.reshape((n_k, num_wf*num_ks_bands+1))[:, 1:] + udis_data = udis_data.reshape((n_k, num_wf, num_ks_bands)) + + #initiate U disentanglement matrices and fill from file "seedname_u_dis.mat" + udis_mat = numpy.zeros([n_k, num_ks_bands, num_wf], dtype=complex) + for ik in range(n_k): + udis_mat[ik, inside_window[ik]] = udis_data[ik, :, :n_inside_per_k[ik]].T + assert numpy.allclose(udis_data[ik, :, n_inside_per_k[ik]:], 0) + else: + # no disentanglement; fill udis_mat with identity + udis_mat = numpy.array([numpy.identity(num_wf,dtype=complex)] * n_k) + # return the data into variables - return nrpt, rvec_idx, rvec_deg, num_wf, h_of_r + return nrpt, rvec_idx, rvec_deg, num_wf, h_of_r, u_mat, udis_mat, band_mat, k_mesh def find_rot_mat(self, n_sh, sh_lst, sh_map, ham0): """ @@ -588,7 +803,7 @@ class Wannier90Converter(ConverterTools): return succeeded, rot_mat - def kmesh_build(self, msize=None, mmode=0): + def kmesh_build(self, msize, mmode=0): """ Method for the generation of the k-point mesh. Right now it only supports the option for generating a full grid containing k=0,0,0. @@ -604,7 +819,7 @@ class Wannier90Converter(ConverterTools): ------- nkpt : integer total number of k-points in the mesh - k_mesh : numpy.array[nkpt,3] of floats + kpts : numpy.array[nkpt,3] of floats the coordinates of all k-points wk : numpy.array[nkpt] of floats the weight of each k-point @@ -614,30 +829,26 @@ class Wannier90Converter(ConverterTools): if mmode != 0: raise ValueError("Mesh generation mode not supported: %s" % mmode) + assert len(msize) == 3 + # a regular mesh including Gamma point # total number of k-points - nkpt = msize[0] * msize[1] * msize[2] - kmesh = numpy.zeros((nkpt, 3), dtype=float) - ii = 0 - for ix, iy, iz in product(list(range(msize[0])), list(range(msize[1])), list(range(msize[2]))): - kmesh[ii, :] = [float(ix) / msize[0], float(iy) / - msize[1], float(iz) / msize[2]] - ii += 1 + nkpt = numpy.prod(msize) + kpts = numpy.array(list(product(range(msize[0]), range(msize[1]), range(msize[2])))) + kpts = kpts / numpy.array(msize) # weight is equal for all k-points because wannier90 uses uniform grid on whole BZ # (normalization is always 1 and takes into account spin degeneracy) - wk = numpy.ones([nkpt], dtype=float) / float(nkpt) + wk = numpy.full(nkpt, 1/nkpt) - return nkpt, kmesh, wk + return nkpt, kpts, wk - def fourier_ham(self, norb, h_of_r): + def fourier_ham(self, h_of_r): """ Method for obtaining H(k) from H(R) via Fourier transform The R vectors and k-point mesh are read from global module variables Parameters ---------- - norb : integer - number of orbitals h_of_r : list of numpy.array[norb,norb] Hamiltonian H(R) in Wannier basis @@ -648,14 +859,122 @@ class Wannier90Converter(ConverterTools): """ - twopi = 2 * numpy.pi - h_of_k = [numpy.zeros((norb, norb), dtype=numpy.complex_) + h_of_k = [numpy.zeros((self.nwfs, self.nwfs), dtype=complex) for ik in range(self.n_k)] ridx = numpy.array(list(range(self.nrpt))) for ik, ir in product(list(range(self.n_k)), ridx): - rdotk = twopi * numpy.dot(self.k_mesh[ik], self.rvec[ir]) - factor = (math.cos(rdotk) + 1j * math.sin(rdotk)) / \ - float(self.rdeg[ir]) + rdotk = numpy.dot(self.kpts[ik], self.rvec[ir]) + factor = numpy.exp(2j * numpy.pi * rdotk) / self.rdeg[ir] h_of_k[ik][:, :] += factor * h_of_r[ir][:, :] return h_of_k + + def convert_misc_input(self, out_filename, nnkp_filename, n_spin_blocs): + """ + Reads input from DFT code calculations to get occupations + + Parameters + ---------- + out_filename : string + filename of DFT output file containing occupation data + nnkp_file : string + filename of Wannier postproc_setup run + n_spin_blocs : int + SP + 1 - SO + + Returns + ------- + fermi_weights : numpy.array[self.n_k, n_spin_blocs ,n_orbitals] + occupations from DFT calculation + band_window : numpy.array[self.n_k, n_spin_blocs ,n_orbitals] + band indices of correlated subspace + fermi_energy : float + the DFT Kohn-Sham Fermi energy + """ + + # Read only from the master node + if not (mpi.is_master_node()): + return + + assert n_spin_blocs == 1, 'spin-polarized not implemented' + assert 'nscf.out' in out_filename or out_filename == 'LOCPROJ' + + occupations = [] + if 'nscf.out' in out_filename: + occupations = [] + with open(out_filename,'r') as out_file: + out_data = out_file.readlines() + # Reads number of Kohn-Sham states and total number of electrons + for line in out_data: + if 'number of electrons' in line: + n_total_electrons = float(line.split()[-1]) + + if 'number of Kohn-Sham states' in line: + n_ks = int(line.split()[-1]) + + if 'Fermi energy' in line: + fermi_energy = float(line.split()[-2]) + + # get occupations + for ct, line in enumerate(out_data): + if line.strip() == 'End of band structure calculation': + break + + assert 'k = ' in out_data[ct + 2], 'Cannot read occupations. Set verbosity = "high" in {}'.format(out_filename) + out_data = out_data[ct+2:] + + # block size of eigenvalues + occupations per k-point + n_block = int(2*numpy.ceil(n_ks/8)+5) + + for ik in range(self.n_k): + # get data + k_block = [line.split() for line in out_data[ik*n_block+2:ik*n_block+n_block-1]] + # second half corresponds to occupations + occs = k_block[int(len(k_block)/2)+1:] + flattened_occs = [float(item) for sublist in occs for item in sublist] + occupations.append(flattened_occs) + else: + # Reads LOCPROJ + with open(out_filename, 'r') as file: + header = file.readline() + n_ks = int(header.split()[2]) + fermi_energy = float(header.split()[4]) + + occupations = numpy.loadtxt((line for line in file if 'orbital' in line), usecols=5) + occupations = occupations.reshape((self.n_k, n_ks)) + + # assume that all bands contribute, then remove from exclude_bands; python indexing + corr_bands = list(range(n_ks)) + # read exclude_bands from "seedname.nnkp" file + with open(nnkp_filename, 'r') as nnkp_file: + read = False + skip = False + for line in nnkp_file: + if line.strip() == 'begin exclude_bands': + read = True + # skip one more line that contains total number of excluded bands + skip = True + continue + elif line.strip() == 'end exclude_bands': + read = False + continue + elif skip: + skip = False + continue + elif read: + # wannier index -1 + corr_bands.remove(int(line)-1) + + # For now, it is only supported to exclude the lowest and highest bands + # If bands in the middle are supposed to be excluded, this doesn't work with the band_window + # We'd need to manually add rows of zeros to the projectors + if numpy.any(numpy.diff(corr_bands) != 1): + raise NotImplementedError('Can only exclude the lowest or highest bands') + band_window = numpy.array([[(min(corr_bands), max(corr_bands))]*self.n_k]*n_spin_blocs) + + included_occupations = numpy.array(occupations)[:, corr_bands] + # Adds spin dimension without changing the array + f_weights = included_occupations.reshape(included_occupations.shape[0], 1, + included_occupations.shape[1]) + + return f_weights, band_window, fermi_energy diff --git a/test/python/w90_convert_hloc_diag.ref.h5 b/test/python/w90_convert_hloc_diag.ref.h5 index 6beacdd6a398097c605138f0f087dc41a5ae1205..02faef6cd4996f222fe6e6613f97911db5226cf2 100644 GIT binary patch delta 46548 zcmb?E30zFu|1%9*ETK|FlQo2pHZ3Dd*)oxYmMKvZrUmVDB`GD9p(HhBtF#MAHIqou zGDVxvwC|awX=?WQ-`mXKdE4*(-@87|J?GwY&pF@qe9!ls^Oe&l+4^Z#)JjBa1|nB| zDFQw%I&%JzW;xj8WseZvgIwT~i6mhb<{czVsEYElVuXn(QNHGbFySrAuW%Apl#`pA z0>|g5N;zDdDtwhHYT-_-umv~_0fSOXqAzE|W>2Or?h-z*7d`*8a(#P5V=RX!UEv$dQobF=)O^O*TfrN(gA6$q`qE+RM88V z!#V63DD8>$I0xNNLqFhk+~o{@kCGPNK^YlcA}gV~1)rp00iXXJ85sAWpK=6s!VC%6 zNfH0Nm=T?;A|vM4)akjov*7qB6;)|n5y@3^b0;RBNF?pX&xl``ThamF`Da}8s9H!L z`^kPmIQbK<)D;%MZ`iE_pMs|;t4S~yifhNaEwO;_{@#0nda{g|{rZ%_-r~1dNlDnN zLaAD1DG75Vl*(BuC9y>%I+r-f2tU4t7~w|_VuT+#)be=KKkgqcLx`TggJbnA5xn@t z(nx%Xv=rXlUkWhYhTc2*TtsnU85nGe2rJuN94F}~Aa$1pkqPG-cb***HiwNJwq}E~+`}r7VBd*)w2UCVLBqNy6k|QZN&&#~BBZ zALqm2lYb} z1&9ZsswFN78<^xt{2XaTv7F*V?JeM2zk?0|>-O8M_s>}&?0p-4z^=bugc>fEWPY#P9j{-l;1>raU1&>1fI!9z? zsfrZ9?MV#;gZ)xJ@c$A5NU#6d2zs6TTc-H4a^ET;fcL;)_Rt-?#0o~tB`K7rurOi-ZY#eBh6vZmz{X^maI?#Ka>!P z!=ScBPU^p%7zmpqMixII1;XFRsqYvAaRk&J74~ekNL&9GID({c227|X;b00Hg)taM z>>0PyC%iG3$p%HEF)*pXgiI$CVqsc1c8ZZM983RLIe&=x-!#$ABSJv!g*TXxRCuSA z3#3~ol@tu-hId$roVBk?SVB}FDQ?LmnN&gR8kA`GNqn6H5|n6oGd>VRI>gq8okWC? z@IDEpK%s%}#4kRE#H(sa!;eFB*s{+8?)e?@5FJ#1n_cPntYZqQP8*q`#XbfU?HXt#rt7`(B+ zBEC&Z3KGWI#uo4&3=ZKL=K347OZ}V`+EgRQz6;)LdnJR5Vj&jc;4jAa$r8Vj?Q517 zBxe607O1A7sq`u-4^d^J+j`MmY^MKa@*Q2GyDI4)e4qyom16xnIm5x%N&T=gMTp6f zw}1zI$3(XE&`&uCO(%c6ggOn`**OQk<0WdH#c$MboZ0^<-~W?%&%y`Hv6NU>_Pu+c zDW>wr(r%P1r3Um`^=|Z&|QT=C){d z|ItxskE~agS}nEvx1WUDh(3mUj7CUImYNQ}h7#cH zYVv}E${r5>JAT6Hw8Oa*lkcDxQeZjgLP$)h4W3Ma&KWQZ{(&+Klfp+y$}3b^+hfG$ zhW9aHDRX>#A`)cCDGE@vf?qDX=v+1Fwj5F8c$yC~H;4I$lT&<-)nkdLS(8i#wu7HE zp>!R*;gxOlL2Lf07TmOz_MC^@XNH|K?_xIhU!b<$Ep)@}BgcO!)xONVzNh0fE}Pu^ zdJHWWefWt@_X-=s5I64TI;2N z!oNUPqneOnzeg76ApMT43WuMO_46G3IUX94pc&Eo8$yyp^CUL@yE%%8*^7wIEym|8 zm^K3ea}!}w|Kwx=nIvKBu*y@y4=#cqw{3~e?ZKD$%*D6dH8uUmd?1_xAS5LCLBh%* zFUTbQB_WVtF!Nu14ACW5Ce^4@k`PMq@877=ogi`36PYXWbB3!y;@;l0>uxobG2c%K zvwh34cOJLm-B{D)2*mE}N3U1~q?@GY&Z=XkW-q@M>h!L&!A;vHByP{51G@u?9UIJa zNxg4TcAxLWx-nqa6VzS>{Cnb_yC($Y_r(2K62uDQuO-12l7z!}5pMoU5~54&58KP^ zThokizBFI`v^h)4$0=vZOFm-C*wU`J+A#!d#pdZV{~zrB?-8;%=YK(ncry4uAY@wd z|Adgkw%%~{Fx_Q!OT4}&!eyQd5}&Nb@IT>{NKK59eLm0C-1t)|&;16~RbH|19CHTs zXlCJOJ%OH$-AT380{p7SJ_h*e%R@8QwY(le!2FLd7kltvldLl7(hd`OYlSXt$i~;M z{%Ob?lu-G31RRuD@pD!H#Kys#=v+-`q`_A# zJ0RAC&?<{YbehD7Pj0R@e9|9h^rKhMo9pG1Z?;^UEper8GI?8a1^`Lwe`y0ij!LmN z(7XP zSyiaL?S;J)|KhjomY%fe5-T9~nN0S)n@pDYNWoR{i$Me0KIX{37^QlS2vs@cGlKrU z59#LrHKX9&XD^%_tER$qMPn7hN;p;>7WM%aDvgg$S}a^;xREgVZbPW(-4>BmA^H^d ztB=52gtQ9ba+WtZa!y$+EG(Rdh(Cl#FQ$n+et6!}|CGv*KX0zctp&9R4F`Y2Mf?BO z!wpp=?lL5Af5t|1p=AdH4OGZYAvOQtfgL;k$oK*cPc`Yr4zMcd!n?=(mGN~HzrcSs zenA}i*#FyNK=4K8KV&gL&`gr#S3JRB;G(PTlsg2N!$ z|8)@ubUDgwMZmV?z&?Ra8@%$J06=VtCyLN6Mv1T*NYWu&AME!38xxX`+=YTE#snl* z&=Uw7z$7}r#PC2w%(-xwD0) zDG1P{H%wp1WKfY`-DIl5?_l>tZZmO(-#s&sb11qb6q>ip`yZ+yN?|E77-VHsa58yT}I@oj+&pY>~DS@jR$Ue$^kK z%Su`}@dwM`+Y_RI|F=f?YuzaBf><+U)eE24i(Yuc@p^VEJ&RVjMHogY$!~(KSK4Cr zXCNZAqjR;TyTFlR^1^S3Q^c=NRQfwP#JjC{k3eAH__c~i1msD8{OJFK&H6VcrLYq} z3nXN+dc$wo`uLlbbBxYBJ0`o;QFdm;!zsB=r_L@}vUu^%tA@#zrzrUy2Y~v0!$+5o z2UzWyiNCiy{FJNfp0&@O?KrbzKksXzPtOaw2~mf{-%oI-8cxJHCr;=&w{nR0TZ-Qn zYoEyOdNM-!Mz~Xj*?2;;sZ&w8K4GC=K#?Xr+q1{V;haEJ}8t(_E0dHWurJ=@~-A}$iNHDrOIQWb+Qaml*6X-0mOWGZ7wXf*O z*ZJryi;e`IM#^Stig$d`xvqgYVD-x}Mh0o$t?ZopH|YGV!Y9N1qkXQ<{&7{#rG01U z%Ht}XMSS=x&ew}C^72~v%NPty>*tdn{rG-&X}O(%{K5j4Y&z?f6&vu+scm*U%J!d!$N`M^jX^y*J6u)9g3N|TDE3FiV# zq8YEb39M^N*~MGOwg9Y$-s#cI+?_@8%h(UuD(JDjv=2;H=$pzpPxuFs2W7zy#SYA5 zE+pcXvQxK?UG#Pv+ko$j9EhmVKmcw;!oXti){1Tl$wGFJeNYXWjQmTtWUd=BKQDI1 z#f?e!^6HT-b4i?_kv+8MJ1akQX5^C=x!f}BabGSeKm`HQL*8;_fx=B%o1702Hq>z1 zWoH>JK%OgZT}qngJr?u7uu=!DESp-ADriy0<-R&=&CrI#5GV*77DW6|WF!w&^3}+h^3N+^XJ~^8H zMKXa^?~0_+9@UWFvh`kZHLt~WZZ1_KkVutE*6uxnhoUqedEdWcMfKDk`sAO*Ykaf| zLHaa{cxQM|Z1b(oR84~EaU(O*0CRT%CS@~_t-f$^%|zUbu?j*_V#%2LIN{rsp618x z*H)3_v4KQB@42s`z}&y?GWFuxVGq0OzT@2cZcKfOb0pc*p>2;&lJ|TcAji}2Jw~e( zD_5#!OIUd?zw&so{H7j%!h-xIoy#so?n|BF1Egx&?H=Pq^#zXZ0VzCxZ%3KQjYD{# z%`g1gwuJ^VEX{zAam|hsn%qsTD=Dn2Dl19y8V&Jm+=6T}5mUU6LD<=26?cuGubxeM zK(Q+CD$hshQd0)E@g!)y%peUQ{-xS@NPUf@fW&WW+mGM7zPz7jd~)MTy9FpMjkRl@ zq6SY(lAqLHJnzg8A!YLTrv(GCk2V8zOUkxJs(}v8Wql1{<|6-r?34aRbjkM`WqTOU z6WiiVm^dP-j+(-c>T@*e=e~+VGV+&SrktOmVT&A4Y^d=4#AS_3N=9~kP67Kf(T#-| z_N^@UQlDBC%P5$<43)b(tdr z)yv8XU5|bgaB3!e9@zremo0@OEF$T8P>eSN~or7J&{_;KeYb<{z zv+ZT`aGrPW!|dWweqEV=oQ_@k3_o2@96vb$4a3fBTr4O`aEbMD90HCh(&ruKy?YZ} z#L;y??>epDcEVv3Y+^ZBzwkS?XKRY zt)F(0fQ}0jJEOy=KR5w=I#>{L_W`rVI1GzhS&1*rz5nIQ`AoaAwgNo!-Xl9X>g z=UXN!KRnKN`eJ%zEz-@f@|1DMDg@9rB6&6oeT9Ezu*sgtJ+rHz}7sF^L%O1;P{cz`VTtOtMp$U<^)|vt!D-;klYc?kU24R4%VU#Ga@MsMAakP zu(yGkCbT=gu4VjX%$`7MQCmFQ`p{f-JrcLfxB3AJh6SRM_~UZSMfaGu?vOu*p}{`; z0-8P8BiGp5meBN>L5+D;4^S()-`v_ODMnLCEO}BFZh>kQ)1529l|KPB)W)jf#8;Id zAAGQ08H7sB{k-Fw>EgZFB-Drt5q$|e07@h?KB%eQdXalpiK0N;3n!2{wiW@TdpUXo zKoMo`dETQ#gd0z%G4F&PT}!^^T{y0yr~mnyVoWjZnN7jUmtBj=)FaQRUMS#U?S=WQ z)V*BEaK(2Aq|u5`6JFQ_)X=|NRMR(5yb`vJ@-2h+?UD40vzMY!XSk0jwFf^JBQ zpGW)SR_Z2aEVh+KeT7qIjX15%`*wbMR3VN<2>O=5<3 zgYM!5^mkWT^|O~F=iU93dZk7^{|zNyfiSFY@sL_Ef7GEWWE(bOex>7u>^{0h-hPj^ zdw}{E8_$LL9T0|&f8p-C>a zK>^y$hqFzFJj;+IIKO@(A5J=>#OwmLU9>9BYAQ$4W#lmG#SVsYx9zgKAK?lTWyh)q zFHt^(VeaVwPH)hsmW7XAS4_WG>XVUGO3QiH_i$|O^evQH%_iSNOy2yFW4thXmHg7- zCj9+3DsT1}wy96|KU?ahXCZqR9B}UVrk84`GJ0vDI@S4&*`>--{#rEW>HFUU)|f2a zm}y16xE_KL&v-Bv8=}HD_ciA}n2@Ya+02~eu816}cnRX*svOnOA<}gVHbModz6pX3~)DFr5Q6rKe*k@T} zE}Z`;s}r;(v>1&CxYGWb7qK#szDpy`Ea5FtW?E1LRmcUmaVr`IK>#=fosun-u-n}C zsEb~YN-L`wyQtU`5sXvIP8-0{G z@LQbr_OaLe7c<8>X1Rb}<|xWQg>#}9sWp)uvz~wn+Oder%TRc4S4#Jn9*X?zAowE9 z>2+<^07%#UFBhO>oSblhT|kEYer49xunP2yQva4?M!o_Yd6tZ(KyKwVUQpqZJsU<* z(H+Z1n#BS%tFbFbT`z-*C-OwxpkhVqXibyPH!H>~ z3cQOfMJl@dPEDt=++zdSi(}ugnM&jH#tLOvm9HMD9M4E|eu9!&M)SnR79_$m){b#z zjmDxotGPY4eFZZ{)mOBt8h3;gT=SNmv%jQCfrYB}n$x8r2(ax&BunlB4POw0fE;jV zkc@Uj70+c|&5jlXs_1sWadPD^l2K@4Q-#j7XfCM4FAfHWU8#L-*M-m=k(>oY#f|CM zpf0mU^YxsJwmN-4y76|JjHy3U+5L*TFWNtNKJt}7U*I@34;xKTy2!nJ9X*Woio)mx z@-(@|k9#kfYZ}=sK3@Dd=U6C^pH`7yiaf_Je^im&i-a z@@iJO7PLaNo?n1qh;(Z>s>w$hOn6)6q;`MIAX_xGUpD3)RCf%-FVpQrHq! z?43vo>t$VOj+|l6l7Z+~-P>?wFO?cixbf_j*LvoEx~^84=v%gZB1{>=LQFY z+KKXj#y}9OD?p%z$%X2vtipm+`%yJ%(%rW19)NwgQU6Nonhg~df`hdzO@7J*jxsQ- z_4E26>WdzWT^VVm&iw(~-Q7#{DJL6_FVC@?o??%QXUzzVgi$<0*fK;kP6EX+`yyLZqmlTf&@$J4GGsL4O!Yxp`)mbEI4kxmcV zE(_ccYR78r#e(}Z7owywQYO1GGV67yb6HX6D^gG`EwdFp2=`|P@r%A(;*!lCN|u2` zbX%=DC>EZiWLzNAqA7RvbPQj6E~jo^b7TvhHY~fj^#Kp30VvHfEP61eItM61jmF7q zRBI@4Y`t8gk*uq$!NtEC$~ivZyn0kY2Jp*N3yxfV3T4qb#!u>uaOIsM4iEaU5Z#k% zFUrOL<3;Qy<~JGQ)kM8fX9G269tl;9+Q~#Teo8?}F@qlP=h=@|VpRo*q7OENG5)+` zUBK$~%Gt4~onJJvj43U^9H4~y5&4v|2g4tSGl%anyFnppO9k%nj8FNqClGk*IlgL} z>~0Ky(HIXCk?G7||JDZJE5F0cHQf1LgIoQ|bs$q)SwE?6p4A*4XZOm=PtQ6Tb)M${ z%csl@{_kE+dl<9Q=G)4c?uL9)!ABtfnk#ppzYe*uBR+P_&gX2^=D2{e>irR@mUmX? z?KtEVoq6vVaugt~L8ob&%6TUHYe^-_(a!SnSCx+*8{Y;W>4iI61Dm|=;1bTPG4SpE zggx|*yL|ULjmBfjBsw~a(6~2YY3#zxn9cOB54lI=>T*VrmW`u3G6!2Jl=C5gzck$?l=UTT)4-N59V!M* zfe$u>8!cSYDR2cmN7%WhuMlX{^P3|E4Zmou_T#zC)E+L?y%}+yL`}PBTcxva=my>C z^UTqAZM(c?`Ss*wGOId#4Cs0!_DB`Jca%1OZZ%t!D8CGvOdhBtARZYt^1uPn>=IV1 zy>P@ik^*=Le0Dxs!&*+h&sXzN1Wk6Lo4Wsl()71k=Iy;+-7Uk(3rHstvNYX~R3-X* z9Qf*0K0ZUX!%Chq*PlKyZRU|f$!C@4vp^))P)Fu-AFal*5N!-Ae-F=q zt-%F2tq+;8Jygvy!a_-$&MKVG0%k@Uv#W=F`6@7l%J9^YyNzm84?Pov%Ui{JpHUl> zNZveBxW7NHFLT;e!s?EvUPXr=Hrw`jrG3c$G9#eZf*oUJW=KDbEjom7dP`OwvZZ>v zYvxKbwc&(D)3YO)jFIaff*!E%dCx2^kefg94mlgkx9Pth`ypg#n&Cln->6W@zO@Hc zi5s$js=dXm>7KaxY=jYJ$85vnU9#8hH2bjw>~AZR8F!;eUwNlo8SixF(fs5V#JDE9 zZgb+jY7b%~>JX?%rt(FmX6MqQdwC}-O|SqRB}aEmTc_+65I(%9@LpMy3c$&<{FJVq z*oifJO*rUU@F2v1W>)c$)_LV5XkUcbn)m}Rw!ZWA;lgS2Bn6jHy12bOIKKyUl`Qbh z+Z<*xg zeBdZ2pcs;w)XWk6MF1NY3l^kQrH~9^Bw2aVmwM|f{Asw|(!<~k0XLn!^CA*(Zj(LG zzxL)AfjXWOP4b@MFB_N*XBtTppQ$_rhk$nMEW2%36(F3|)=L;CJ|HMwIvi~&J6CyT z1Ft`D=~ymY*YeKEa;(ny$mAbz3}1A3 za|rN)xTTd&pDXfKUlyvZ?e)? zA*@SZNz`NZ6yI0L7=HOonf&#BC5`nFU(kb%??6Ag%^j7;i|wbkRV+^UP+N9-2U=w+ z>61L^CfYfnp=eC7^pnXpZ1NhDYJ#&4L!E9+AvhSNC-(w}Km$L2y(aFoNEix&j zbhoT<7z?-e<%O%Zh8lT$^(~J4fMCkbxr6bWYBT33;GHm@%>>>io~NM zJ$Ee$a3EFJ<|KuCFp^|OmDz7)#KcyEKIs4$7bKdsyebY%YtP86~8telLuV@jh z4)uC&ZT4!^&~S>+u`O;*m!|8C96il#>!G0gd=9klx^R_*of7N2(f#29>xj!h;!QON zo!B&T;nDE$&a=-A8dyqpk&2yUuG_ulH|P!P%0Lez6{zZ-L&-g!fj%0h#I5X#E#Kbe zbJkRb<6#<|n#CZNGVX0pOh(x@yPZV6y}ITAtx|zhHxyVJ?ygd|#Km%;gJRXjdv&^} zHhQgMS=Fl-GJ8(i99sr#l%_^HzgmmK_6r!>uPNv0xBBiJEZaPvt8G8c&I_IE|1O(a zoM1V)))=TtH5iGwv3$&T%9DJ*;*DEt8n~NI&=8N;tx9It2j}>=a(!1ki-zp_R?a^=g9DOaF)@8eUpkhtu;1VLwi2Aad$5^SBt7p z)!B_DYcK$(QewoTNY{o51PWmXcd$JH04|ZFpu=2y=LmpbO?)rS40>C>LTMsD=^7VKz?OkME0-FOtt{#Us!6iy zLz$TLsrN;anMUG`I?}!>dJJ;*dye-frzFZH_H^)KQa!!Q)>@?!k%yeDJ334o4$B0* z#2>+H$cC=cYI;0@_KP^f&YCB9lHaC{aYx_Blr% zeF|Bqc#&%laP>ui3QUs z%E@kpd-BX=ve}&D+~Q!qhLxnY+vqzJOyxlrvVK8G1ritzSxAZ~19p45U6Vj>DXI9{5CKi;F$Z)N3F43ax zi#J*0kvCs2($ne7hK$FGq{bvAhqh-e_5mY>t}UNC-t!uOr+SThY@5Pn7h+iSA+LaV zdw-Xv&y<54T(R)3VNM zt$I1UG>gAFIAWlG!jZ-C^j@mLdujS69Ny|XQg4bG?P)p%RIF@8;MDz2xR1;Fe5zqL zSTF3?m!kx+3G`=AakPo*H)BIAR<6Avb=v9lSJ($iNh6P0dcJM$PlY|4KK-<%ib2;vM}Z8lxl&h z_Ix_cdDU)N;MCoiH71uZ9GF3yF+V$ zPfM+Qr@N4sW=?0)`oT$F{RAi}Ns?_VHO{WaEobB(dqa6?c_e|}q%yKd-F*A`X07^6 z-rAK$38WnvDB|*yD{$0(5wph)zHE@-ClLYjYvX`sTQWZ?qjP4q>3SyOE#&RGLdssV z(ry5|3R~}^cc!*?fR|jYBv41EZ4H>!MrVBTS-k8r0aa>E$h!w5xRpX5Knl8u zYNNX5ogUy{ao-F9%nrlE zr~ALTU{r9FMZSljkQsc26L7uA`pu=LGJ_#3owxB@JTP~47$t{bTh}{Ipzv$rT97HU zOp11EZ1HCYqtaGXWN(U-mZ0A}_LEKb($-#6?M8y#y%Fzyv;l4@!xtG_6Y${i&aFm) z1tE!_J{E5pwC`1OdVjRi2G1!jI&73fSF}FC=y=WFmB0*Mc(&In2B0-BdFOb=tF`Jp zzpef0Na<4a24_K(=U8n~0B|lht3quIS%fv2>-ydxMb>8&mDlx>sza@zoy*;uGm&oH zI_0#?d473zHR{uBpUV*|Yn7d8y%A&=Oh{jpT4RfV4NK7}tI9Mb(i zwhFRuc#*Cf78Z2@+Hs#-sn~A#kqIrz+2QV{9D&_gy7Ic0tBQ(_wd*URs%`PN%5nbh zHxx_^Z_LhJJI?ofvhc~o;M3r|o_C5t{K|Z@w$t4~6pe%NuY-m_2Mo(OS&_)Cq~9HV zoS>ms(WGuzA&3agaPOu!qvn()5-TL4NPgKzhGg$k zW;y|jF%x4>)R^h!xxTN5lq!h9?-1hNl&zy)sjW8Tw9rsn73B8ccCo?UXBLb zjc?dNd}`CXZaTQ`b%{%m!+NBqVdAFyR%DN@&%Fo{E;g2eG>~=|ev+%CnM^^9x zUnAf6^tD?(1np$wQsnU!JhM_tl*YSe%I;=P_pnEkb3XDdiHHFmkGfoAbAYMj+m|w? zih`5Q%ai)KX1$6}Nx@49FhFi%nEPh0o`*|Z`&w`g${O-7+o+f6Hg@r5P5g4GY9%fJ z{!U?pTy-jPXF{&$0of6q5O4*RZN1+QAYE!|t+gNcFw2^?+D||V%jP}mWCSeFcOSsA zak)pnYJ6b5p7oLJJ=x_cjC@i_wl%FXJl^|Z4f#YDwkbNh@Iu!|NBUcoR0_XVej=O- zY6A{vK3D8zgTGQ3{&PTuo3)?&elJ@&gd`n+x`L|M^o4Rgm7`29K`!G4zJh-C15#Uy zq!q9RvSa#|{jt#cZ7F48Q9WghueH{cbC2LuAWu}{(DuWFagher@zO)iR*g0Z&&)dt zMvTY#-Y9@^te{}qb6mM4-`NPGbG$CJyZF#(`Z&%<8PF&g*czzg*Y33jZFj)1sa^0b z|1d!1y+~2%d(uK4z2|PWdf!=8Ecy9FIsQ0>7TA37ZeN>(4 z^+J+MeL?On?#Tkd=n?lPZ>wYVBUUveOgyJ6yb}ri6Ovt}wTt!4@#a=acJ#xurq|A# zb1(ejKovB-rxRx0w46tc;*aA1^=#eyR17djk_h9k*h)_yc@!Ku%}(oVFl*uGTsm3f zWldVEN$Zmdf^-os-C($K@6M=O!#u#&D5V;Zy}4_2tseMfNAJrmbEV$=OE z2gMF@zwl}pvz~p$*!o}c?@U{oKrH5;ciEuGdbuiJKQimZy(c?q8t&}prn1|Hc(uy_ zyQ;2ORL08^>Jr!2p3R(yDOEZP=*9x8qn^)e^eKJs22}-Razk9QKhdP$3I{jrI^f%8 zT|%8&rjg=(9+!Nk@}-)Ut^X}cQ{tuNh@sTV+>_q-y7SRb~5) zh7GN5?K=9#I{LMzo!PoNTmGqDNqo0KUvj|*#Q2eQ8+LYY%T-ys!>7%XUfWB#YX$jU zkCHad_f^Uoi;9LEq~r&rUChtxF#LHA)Y$xYZ%h2g_WR%!H(`#%o%jHX zA)IN=8U+XZC~#}d6ejO5O&xTTE#yC+4y`lTP$<5SCxupt3#YrhYm@bzAuu9sE4RKB zrqlYVUPHWSl}3u@OGSU}6UfLrU&F<$2gbdO8w3jBID%y#fwoOhwRP-%h6yU{6J^et zCD-^d_H|zi31wmO>b1!GgDe0nuw;yEvxQVO@faCwQRD#lQ>{}j5z>1a%%%owm?g9LAvYaQX7G7s}I5%{p6L7+8-K?A_;B25`q z*cE4Q^DuWu6J&Xq;X3Na>OR4-mx8P59CMWeSZwZ^zSy1HNpE(H0RbU=0IXtJU#v0SCdhy~Z!=uFXfWm0xLu zBrew((O%BmpWFd@Xzt{8m*mgl#l)Hg?92?g-Y2<bKO;M^HaQ5SOdwIn>*jF#@1VHcK#lan$eQS9a=QD%q{qGEdt6$NicDBP-jm1-vQ8w#e*3@l37)(`d?4RI@PG~0bx8UWwzGt#AogwuPJV1@nvX9|(UQJ{ zhdz1faQ9dAi3X<)ea^eK0afx1vWxZa0Z#gC=X1qAyT&TbAqFIcYO+ZZN0zYhL@OS~ zK&^kmWgGj}=`Ye|WeceM#ZB?I(es8Xl|xsK$@TMReI?wS%G*<1OFAdZlWRB#2x1?A zQtMRgn^Vg(Os}KdXEpM60rx1qmSdb|rXklAH%xoM-rH$8c*)?ZvpI~Xrfw-vY~rC_ z^$GgX?0`)`d>yLkOP7_l4k!9SfzdjrER3pVFUI(|CaaTjga5LyO>+kC+$p0m3f^_? zA!Xlix#*XIne4Tf`DH`3K|msTO<%!s_TWmpr8om#X^*c-Xnrz)@rx(SY)LM9SpOn2!1EwLY^u%my9x*a%J(Zc7<(xiR(* z+WbRyg46rlIO7qT|5UppNVMlUohCICn&c;l~dXxZv_JyfFFXlGQ3X1v2$Wr zb{4W8Tq80Iq9|BYw-%_;p$d5GIwIsFIqj`W~s=!^=&$a$ru&GxXPNH|heg26?!A@I&=8guzP9BFy zi3+S#4BnVRT5e6|#fF5qeMK&zWFGo_E2dB5c{6Xei~n>S*^)u)EXX39JTXyOM3znW zowAEIz$S1WUhm@6TA=D`3i-}nE{!UTRnIJWvcSnk!H7dDA#^_j1$bx=+`?t~<$1Jh=~T|8Nc?Cc5Jj?Wf?#~(J{ z>1PuPQd1n9c(P;%u;cYswDWaGcs)*$^&~1iRZpv#sj(;Y3CKJ-pV{fwDY0)3aqWx_ zHwfG;Db#uau(N^EUN3MP<>)psEVcGy&X;BgqF?e;=$**0NVZN?KP4W0eTE^W?ZT=@ z9Lf-zAUSaoY-%d5=^pXa#(_Cm>1<^01TRAWi}JLB_xe1&&xh6yKgmCJ=7J&v!klp6Jp>hL{jsqdPXfCg4{sU(8dnz8P&q)B!Je)A zOnJ(QyPz#&CSa@6_^Oi&_bb}5?L=-w}hoe%P5EyJF-UluF^%0QN&j(_OuBw5`R z8LRG_6SndqKs?$4(qOs3x%`q;)ZUGQz%s6JlfdnVRlnf?I`0G6mwL)h9(ps-3^+9k zbQyt0TyW9>;uWMpsvz2UgKVa=;qY!qUKO=*MjU@eGy9Ck@YKpl;L)c0yahY{6yj1Hea8 zB!}d4ro10W`Bb*KUSMG8^MD)W`|!=#dclm7nvZ%Dho%Pe18kc7qS(l^L11?p8)-n^ z(a(+aEz0XCGpF+Bkqd^c>c|Hcl6q?nGEn}@2Gp|d;f@g>o3?Qd8cKl1ArJj6BmCZh zD?JMDQUzBrPjVUg9Ns7CdXT*GfxeMnQYV-(^#FQh$*yCeUd17FMn@y>wcejTU>K^R+YVUWm;~@c zSoZW!dbF0T<*5Sc+lo3~{_XEiU-9CR)0`P7o;fv|+-3E$^dJwtxPED>PoHOKH1^yx zMH4|{gUQ%7QsRKXV*t=i6GRgNDdDNv_aNx;PNXRr?`)0xeD8GRNm5kiVz$f(!U5ZO zM{#}0fgggn#nJ4Zm=NjE1m*Y3VH$?NReoi zuiv<<!M_Ck)1Fr5W}zz>-1gS6LU!;|sEVZ3&X;g8;Z+;Ua2e zRw}rmX5%Qd8xdUXmj)cSlW(xCmkj`Sl2nYV%TrM60cOxxx2)MKDl#>Dv>VIsQ@cdE zy7_@oFk3lI&@Kbq{seT!E^Ct=iFD?VtmvplGA=+naN3~FjkCZV$_$+Yf~6QukkJ`B zpw67hxL1#S^5p>I8gM^xe{~bs0u!O;+qMpmfz+y-tpFcQI?u=w=-I}d=#u?t20;6~ zGGZau)D?dAHD2a2CT|ME6720P;3H(&KGFDwMeL{C2?Vz}AbKqJbLN9ivcP≶y0D znbK;>QSP>C5}cfzjav4XJ~7%@BcLe?h{Nakk@`0bpVOCf6l!XtulW%Kovew|4Y>Ke zf(Pfl#o^hwYfA3+PA*}~BlhBRivpB1SqyGr!CRkYb)^QBrI)L&iuT2^e##>};2}EL2QiZze{4DlekkjW_QOMr-`8nnMRetyjRj ziYg)t+)cw@UckxFVp^)}v(xXVNqX(xJUM9tXJ2v}soo!%Z!Y+uGRYdWu89abaJ0{V zr4<1!Pg3p`+;lu~jSj4-?8oMzqLNW)eLdx%k(MK1 z&I@ZHFVJQAm0>_Dow1vW^_qdZA4k%$LoW)cA*T(dOytXe$O#1Zz=Ig0Nd}iEtN~X$ zyn#*hZY=mh1(gh%A|!|6wbOuF=Iw#TYBvp-C3Fya%a}wOpNh$Oy=N>xcSAnN9b;$2WZR^m-3{=xKXK=mcxi zA&J4xfW)BBp1Akg94->t@5P`$zJtEZtZAiy)!_o--1(6T>K;CXiMi`#n1!M5TX{Nr z;`4~DJiWvh8ynR!b@ZL7D8E4J_6v=z0-gPx^*ifYTdDc&l0^*r^Y*4itn(y5{h*PA z_s5sf^bK3z`N%n}DFxe7w`j<3;3JwM`H3`)!6 z(4&y*i1$g%Pa`$ZjRd`Su5;GSfL04!JmKKR*sDIJq>sLX-meF*f@=tG?f5XDQwFpM*-!hi;mX`^OM0Z z5c&@A+;kp%-9HKF+P4^=;R}Y_17pvEL8avm3y=ByU(3L_~z*ZgmfGOoM!+r!PjIp!@q?*{`TWC>3vuMG6xAOTr*}A>@98$@s}b`2Q|!k0w@#Md--( zZ^HIO=SD!`e>X-6p@u012|*E^uY~;?Do%Xx@=z=kqNr`Qltjf$Fts>IgaNn3zy1EO zD6~}@9JnsXgRwTYnm_9PIg%yiW_alCHIu|-X9$m5% z0MC*}ypo2g&Oj6%NpU+T;JpY*x?7f%#JhQU&#`gO_SzFRw(KEw$~pra(^|vzo$IV> zf^_#Mt)O@$Ql)xb5&PXo(@U*it&BJqR8U2To3jpv{*zealSkE|^R3V&4jwKh-LKEZ zdGTFb9tti3UE=5BNZS8Hta0H{d2vOPH2tX}@$-KcZ2X_Z8W$ba|4)OB{{(;-6Mq_O zT-5&WP$v#mBxWaZP(slO@ZXjHIoLQv3UPTT!~}HxT8}u^_{LvpCXTl7>vLf>u)hyB z4hClV;}{<&9fWn*sHs4543bJhYC8EgEfH*W+pZjFLz&IvbCKLZA8okgy z5GpLOa6vNDs3=6aIJyb=2)`&^8lU5-8%>nTFQLLu!U0#ea~Q7^F3u; zI1b=m%lW66Zc#|E=@!k52BVsj{i{JQ_Gaanb3<}Wt6A2}-n46KU;9H0BH?nO#xPrA zt|5^kff6-^UQ}DHa1q=3W3$Z~WveAyl(X~Pa~{n*&+gZ*Uwm^fi0D1%cYo(O&vTyV zIeI?N>enNK?ZZ1o?XftK-{!W&a_XNB0Zt=x2De*H&N)L7iB^A4o8B(c0o zVtEyM$P+!XiF?Eo5&xrgtJ&ydtoXn7^t-mzyOn-MgC&QyciQN^kv;!Uk5}Q5#~Dt9 zV4$w-$A%NX^xs`K@3yzUX??nV`Os(2X20>t;roY*KQAjpY5#&l`MaMSJh$cak>ctP z)?PFGczwwVgkWrM+v=+TxzB8aKQ%0tltvu*nP* zAN82&iFe#-Y-EUPHa2zBYJeu9pOTd6Og7W+7%{szFM}ja2PB#K(6;NYp7YO@eRb@@#j*M`ztw&?im}T-A6=Jx>CJeBMZ_B*&i$b@BBuuqrx=xx zB9*ZdfL6(zJ^@QSI`@6SAc8_e5DX$HGz7sQfYg5Cnq=3JpOph@guvbtnu$=lQJIh@j9A z1cL|)4M8x7pwJKmg9r)@K`@A*&=3TJ2nr2BFo>Yg5Cnq=Dh)vf5fmDNU=Tr}AqWN$ z6dHnH5J90K2nG=p8iHUDL7^cC1`!k*f?yCqp&JJLPOB;?P4Q>LPHP?A}BNj!61S{Ll6uiC^Q7YAc8_e z5DX$HGz7sQfJJLPHP?A}BNj!61S{Ll6uiC^Q7YAc8_e5DX$HGz7sQ zfJJLPHP?BB(S38AMQM2!cTbg@zy)L{Ml5fLOjv~CK&|`GcTxCSYFUa|(2nHuk$5XFf^>p-=q-A7z zg}!FQ;kwZ1DHHpaxH)@oCPhnPAG>cSWminyE*kg>dhzyp8ikFQ>oIxzvfTY^Y(e%c zaMH-=^unKS2z6Ea7FopJ7@0U&e|y73xGCV+9^+B3-N~V@_34Yh^wJ8K?_cFteN7d% zU!lIezDY$;%Q2gjcN`0^P$xrug}zzk6ok48(qA61Svdz5UeS^rE3;5vFwsdQI*B># WBpga(yG2F@3*EXu8Bx|WKs2bHa6qhLp*8qe-2!!FAq;na`alfHE(JTU zOv5R4^2CM{kJm_2dpL%(Qf2w!*dvqdZe}Y<*?9N*88~%V!HGE<4z5>>eyUUn#}yE@ zj+c2?QDF~h^T0M*oiSY>?JdJX!GSeO8=#SL_kQ@`8%KwG`Obble_!#RZ}ge+)|>X_ zdyB7eIq2?r`u8i6r@1V>^s5DZNBzTPeEwn2LB}VJpPS>J)elvl`sTVRU+w+;lLf!H z)@Dx7d!G30$^AFD9Q@sp4Rir6OFw1{a0kJQ72*Qr0Rw<(|LO!#*7ISLw)A_S(^UNg*1kuWzU?r@U_LzbN`z6 z&)57G!>2nwJN|xK8b#Tt=DmUuYeT2zHVpOr{F=GID+70?tlEErZ%}&CULhgeC!rF8Mx26dQq+3(P4Sk#^Zm8-{Y(wFg1Cs!Z#Y3pmn}f7YUk*~f z&A<9A6hbH#TZcNfPnfJoJb$9o!m&b|3t}!Zl}IHcAO$h&@Bn`o(wPT zC+G*>rpM?5DSj?v(%;yC`zCAcY*-4B%35d(fh2v8N-pprou#yUW-D?qkM`^ZcGWS5 znd;rOt8(Q&9mv6X-Oej_=-r?{+;10Wx9n=(Lp3k#Xfg&GVYqhBt|ok4 zTTY^QQMLqElH_-=VZPfC$kFc@dR(T=+!})YMeWQM8fl=E$$YN#IE=4kEk9p70OEG# zAp7y|k_$#e2F0VRnRU*_7a3e&u(m;w6UCg#f%&wTtW@N{+;HqNaZ#0V_;=OJaLApM zsk!JqXS1MTU*FoPYj@P|+0tAap!-)C2pp~}Qlk))0sO(tg|h46H1-)!%$EF*%l?_N zqjF?ysF9@{YVZ*bHME37bu6Y(lkU+_hvcZ}lxB2~s4~EsXFD!fKO8Fqn_6yxZRaTo zaPJE?E^Z?IT)i#|EMBjo6&!2%V}KPUyLXr5AvqyVKlyCt5uV+Wsi8Q{_2Twg$6=>o^x)mY^2Xl^j_uxb;t zr1h)wUVw-4~0$Ibc?M2~uJ9aizA(3fvCE=YFHreN)k>IQ%^=mVa4(CSN>tC-R4_>{g z2)j#ERACH~6j>jR?KBFdCOV{1p^K`EFAQg_+S`z+x%z`l4UHv>k#W<;3piQY#GkcJ z)(ey~UdFD4W1Z{*oV@IjT_PUNnd>>uS^4vfGZlFmWh$zs9eT&=TWOEn3@D0gtH?2U z*b|TF(|F<;t@4akc}J^!qKYzS#oV~eIU$r7b6vF-b50-UjJ5x7PB7;D;?96Mq~Xo= z3y~e8`XZwaCF>Y%3XN#)iEPOpUbs#f*=`)tN1R)Y1nmE^djBL*w&YVEm zG84r&E?Y`^Ig$g77F){8-v(P<;$C^24F)I^G&e_f(;ltdJx#!vJg|%QC}WA8v_~V& zhihrVvq+1GUYAbEWPrwAR5YTh9#OS4KXV>s&Zm`9$uxg(WJ=4YXT%2|R zbmsJ_Gk~3I``|$2#vZ9A$2l0t7^Kk;dh>Lpj)P9wXgmCzync4Jn2KJV&LBwu0eI2vz z0ewH5ha7)}_DCg_T?c88o|fYPy=l`O`)p%67VJ|0qg`;)*#3ep#b!Y?PNu6)S0!N1 zs0V8MtIPqs6xa_jk^hf>_{?u@vhydg z*b0OxaaMep&Na+%G3SY}>D(6*oV)JryeDKjPkRld`9e>nt1I2Pr_*#kc*=B61F7p% zBSLNN$;_U7f9Nkl`luK+@VHN#OIqIoAVu?W z?W|`=rZ=9;NH;{%u+SD|N@=i#hLN|_OAa=_bE6-vGeG%r*tn6w0ID?{J3z8bL3-!) zdK6-q*1~cO6X(&3Jx^Al93iyts_cn1I?Sf9?+M-rYy_J7XUxbw-OBOXovEv6dN{u1 z(l6)BP9@wMztpG7w#G|8n=h|0#o`Np+$*0a@8ffnLADQGJEvG>DR={sQe`~tMnU36jSwB#!$WkE-TKafYx}nS7tKBk=l}o! diff --git a/test/python/w90_convert_wannier.ref.h5 b/test/python/w90_convert_wannier.ref.h5 index 3829a31c0a07154b2cd76fa9f64469ff951f5fcb..6f1a22684e48bf51f0ea8076bb0ad64f66d6fd84 100644 GIT binary patch literal 106377 zcmeEv2|QF?AOA?Dk|wEC%4nsPk|trMQqd|&D48fBF_vT-GnI-KQ7T(#MOj9vD9j{l zwxKMQW$a^#!5C)D_Mb8L_B?r0@BjIG|L@~{=K7d(&pE%deeds_^EODiCO zn9Y#&l92U00PhL>P(l2-apQ)POaf*^aj|ogg=K+HNa&mKr~H?c97FvvvkK$?+rL7- zng|I707#ko>ly_BbTN{@rh5m_c@+xj052q-B=kQ4!WDoB@&u6L=lSXa!Ff_Yj|W5V z12;-rg*RZX#-2jIn-@UvA91>ozvr*%L;N$H&n7}n**sK}k@?D3T3FTJ@fA4Z{6`x0 zmopuGJOV^`@SnyX<4&*$k00>3FT&#oJYxRG@%XFrb>hE?i|=0^`~DLN!FBiQbYbr2 zoP_RF$mXFNF^B|YCVzC;fp)Sp=P#b-iir(Lpws!W!nsg5kAU;ycz*cu3qqfe=I8S3 z(hoz~Lyid@ATME_tqX;L)A-JDw38Fc+1m2-Y3HFpP+kVtdrCn16hUEsx>SJswSc=0 za8&^ptS?9h1dr&9RIUo})WL-2XV7+nduO!M&{Fd>5OKn0DENCj|J`zDfEoW=uEOTg zpC41^3X2xW3j!R(#omo+JcRQ31jv4%11u;KK3%XsAW_1*UH66IU_bjT=#0O(T7kIo zEa459-zA`*KPeD~gIuH}@r%s(QK0^V*6?e(|5ohp^pf{?32UbC->RUf-2W;9zvB-E zFd+?_w*ISi78U;g4}p>Nb;dh|bvOo8WhAVymH~fDsD<*;ctDp1fwao~x|#|=^fB;e1afd+YCmiQ0w5UdhnwYp zP4{t(Y4BHO`$iM3w zKtHg}$lp~IUOL)(dxq9(z`z5p+d%nYo6YedJqDO=csQ8P@Nm!%x?n!Qwm9FOM0oA7 zaURNlBH#zZL3(hRI6NGD%EQAakC2bhkOlHD<$~KzV7ozm5x8wNJRICE8XgWVkB5gt zht&scFF4MF^#{X2b<*%~P!BXb955I7{+z-3&m5uLSwQ&kdd?mpTy}(f!1-)=eE@Pu zA_x>1AaX?rh!7AVAVNTdfCvE*0wM%N2#636As|BFzly-<>!t1=5H<*0p$mOuR{skw zs7H883PAX=4I>Z$!Qgs%zuK?qK)t}`^}nY3yT^_`Us-Pzws~yw)#&5lAi{&s5P#@@ z-NrN?zuTvsv`xtU;P~)&@BgCujV=OHA-Y2uib19h*<+$m$I*5Ur!C!3{8vEu*YE5y z1^QmlE)&#w3(|pleZjrp0{x&~|95|%8_e?(%m?P5oNK0OCm z?2S$NVz}MH{bBh3_%t981t4I~mfMhRjEg92%HTrmb7U99~Ab-p&bz>S2 z@V!I_xUkos1IC@vcrQ^f*Nim3KmfaXaI6{a`2}oO9`F`fop4#tU)66OD)Rs3eM`Zb z9n60cU=67RtOG_Pe}S#nNcM*e{h{hZjtQL+>59^Rq#N!Bl7F{P1pR_%Cb&So#PIb_ zB4Ai3CM1Jw1`OAx0|s^AJfgctNC46e=dbVc67YOLy5Y}D0Vv0J7~OE(0kXfX1M5B7TqS>zbNDInhoWDf8=2=E;NfMm!4 z?axR5Iks8|%Nd)`arE)%J|)8A&+99O2oF9({ITAI*oye$&v;B4(|G(&{)C(mK5bC0 z5c)*F{}cp<>rXZU4c|OW?)#no#6VA&(8kHqL2w7o`wGI~fBP3q2PPPv(v=c|-4%~zw3$2$=oe182y4$m9Yc#Nib1?2*bhg z-`D2Xbl=DSlfS0>KK@(&nr^gXN1w0ejA_0aeLNIIc<`C_kMUnggvSqfsEF|R0gsJi z8jsPG5B_{a4JZ$Q_65V|fzcngN`RRXoH>OhjurV}^zrySrt$cle(AG=@cDvr$Z-A= z`Tw^O7)f910PLrOtsiZD>FD>QyO7ba6A)lNyuT;`O}X#%W1B!c z!PoTH0nw^}3+@YmbU&v4SD%0$v^V&({$RLV@H=^7IG?2qeRb_$%L79|K5N|%hWy-B z0D@!A{yZKG1>2kF28R6H1^F!4<>&EW=yc%5cj&Z6F!Vp~j{f+iJErY(Mjwy;B0TtP z{Kxv)*jvceV2>2~{<(7bF%cd=@Q3x7#$))rx>3m95C&lX33wr+Er){g!f<=j(U2De zX9N##7bq{3cz`Vuy5PJxTJpk=xsU7v<%Qp^6Ym4|5DFt_C5rhUMqng94j%TM(b4AP z;pKwk!gTD9V2}2=F#7V0!~fLy49eZXxF4Sb`24~3jqZ=pf5A7hyi-XMd}lb|FXZ{3 zPvYQ&K7|_kGNG_63lEg5rHz&K89>WzjdpObvA2TQ^S}OX>F8o}_Rk*^cQMmZKnaM8E$q1V(y)p%p3w;@EHS5Fa5P?8idA!34-gAVARwjQM#&Um)5S z{P4lg$-Wy_hS!foAPi3d^3w$d9}+_`W@F*af^M z_1kJ02mVL`J~z+><+&e6%^%ZjBUYx_^Yv8Ol#u0mm?!Uc>v)AeyAYlapL2K z3Xp^-4Ee=`#qsY2@nVn@mM2e}pRpW@6XY_9Pbav<$Bo|^D5MqSDEJo?{e5{deE(`R z?dSyGz8if#UK(&x(Dx?rOMyU!x8v`sml~?Q`nIZMw1~9hgU=a>?mh{Zoo1AgT*iyb^;<7IMdT&iUJOuKfccQvVj`P6y7s2~Za% z+W#G8KxOgWFT4?IRQ~v)@7rN3`(fnYPt-|>A0G3ge-7X^P5m)^xK9BLxC!2jO?M=X z3ZVEJu~8abKSqHHBY-{s7@~e86Zu34h!7AVAVNTdfCvE*0wM&)C<4PROO$|ioVx#G zgH!GOf7IeU4Z zd#?6I9wd~z&pWPh(O1{{Uix%IwWh7Fy93tlRD~1tsgDB!veH};_*~APin%VeLS-b` zdd;8cesxPo_%RkxDttdJ7t})hxR(AgjmRfLK!kt@0TBWs1Vjjk5D+0CLO_JT*g#+m z?)l8e3!(kL+VgR<5GMNFo{yW@aSOF}de?>tDeTps8mYDN zo)h3u$j2Ki)J9sZ{ZXr5bIoaBXA`6nt=2#!M6L({5dtCvLi~EA zk}8f0v$k9SGLCcJ@WENvyY7%Yd%ka*SCm6N*tE_ zOBspa;S-E=!nj1>P>9hud_rLUCRh*(Okjt-Ysdp@vX>~lQvd@*?mr5F;U=itvwz)u zu$!pd`ZXPBf_mrOujv$k=qDM!ru)0cj=qt3)|fUjAALMb088v24F$%E(f#kn;hJ0_ zhyU2RzdP59>Ni>ljC6Rzfl$AE$aAC&|Tfs_Wq{)#jD{d0SfupU2h@L1`e!%e}! zdCmnmR070U1-JoU5zqxA{-Qhj&tdPFehxx;MDR+P@G2UH|6hdpAD0YVR}4QKWi-|m zf^p|}^O~-8@nA@n>A}@8vC`8v6+RPT@QOXd>u>uM2>F&;`TAA^-8i42}W=o1kdH zzo6*vgWT}(4SXehDKJ0m2igPq{xbLXfT5e=<4#fOkMSbkzlXqR9wH&bFA|h1;z`2t zp8)gtXvh^a`C0ysZ%4B|vwi`gV7&g>;e^wEk><})%?*F{2|I`&WOzUQ-eGrX?- z4dUY_szU#$y5LSwRl)x+wJD>y3ZBg`Y&af5j^js&E)fDE1Vjjk5D+0CLO_Io2!a0z z1l+_L&Cfw1wrmxLY)*uH0XH8feHW{(hngui@ekQl5CB9&1jd(w3w$g7cjEQ$=)ujw za$s}NL74sb-~V?pTjE)_FgI|~BdC&~`ag);!xurqO~{Fm=|IDGb$zeWPj3eqU|0jjfjE7yI2a*vMF@xx5FsE!K!kt@f&Uo<+{6Oa z&6J$e9Q?Eb=FD6VQJJ&ZXe1e8{f58RY?=bNB`-`5${1iMXcz&y;C{DA2Kj4vi&BXY z5FsE!K!kt@0TBWs1Vjjk5D+0CLO_Io2!a0w0;7Lgr1=A3JATv<{OcY|z<;B;zbzBo zkVT?|8jxN0g-J$ZLspibm7Ca|a920O?J_z32e$2>G)Y%xk=mk*)5ncB{tS7rQ&Ly) zRO^F*hYF@1RKrtBA9ruwb?%bXUa^}A>z__pbYa`J@!L5Gu4MJ<1x@u>ZXm_2Dm979 zKB>3P`3?QPdulWOwy#q+d2lrmzhYi}=^3v_?~Wi}Uxh7idwUAgUTxL4a9}X!o&~a? z`&ySv=kk!_2~kxY1O&I`;mY{faA*7{#JBzq%B$I_B=jpd;X}oqs+x{gWLRgLncDl)sC3XTn1y6Q0cPqI)L#~FRr+HhWYsV$;u11U&Dq3G>hFqZ& zjMCHWk-Ahfl<2(MdS)AC&@ByyntLYpOEkP@KnGLPk-h`(c{HkrDuIrx?J~n1@7p<; z*0HR%t8qKmqhuhTeWDr;S!QC_1xx;HyP|g`q>lAPnVk|D{~t;iZpS5<;nuiM6WJxCb5f4vBYl1 z)KIkXo(n#(b1HKyNhlKSO0F!%0qIYLl95f3EKIR@aq^}d=&fGwj%IFS&!hlo7t1+a z6PFm-MRVsmA~RuK9&aMuu-%+-uYs@)0;wWjy}gKsB)%s!1?aCc(Pj`#~wPc=SDJnpy84sYANb zXOM#BqIgAcSK=G9uH;9R)4oZq*h#rGNiCk7yj%k}znOr0ZF>Ccoq1Z~RN067E~X9S zezT0VEY?8^`TTE6>LJbVwPdPSO}ajB%Q zUT}(M`hhonTaG9p8j{JcFWyafu@$vu`dhR^M$rp+h4bu=r`R~R$gRpZ%9yTHxpl;* z^&c7w=<<_`q^nXV#Vc+tyGVFj0{6*I-8#TDK_AqN>(GPsc7Nl&1JO=nn|k&3`bUd-cR9qGLw&IHY{P*yT~AS_N3c3 z#M*xDx!up*7PI^=dC_O=teKqYMAc>3TMeEcvWzOXX-aOUD%j%YQ*JcFbsihDnm0Dr zlg)DYkNQ>v>rPKni$^DWvVln_M!JyrLgEMs!FaPFyYi#lqHpY-PZK}n2hDGg`NorO zedJcD$_*%le(g8fl0jcbT<*ugNicVCH>X)nVRuLMc<;c}E9O4uP&VeSB2HND={2oU zds>5dEw#QF!n~g$VM8UUGA|E6rL)8tCPBwB-moQ%DD)B5tTm^369UQI+~pb{h@^V& z@vw2tavAc1%$41gy*HchEQsdHm0PGPqc(_Umf-i-YWTkeU@LF z$=ybkJ~%5m-0ogDEoo{1^cMRRt;+F09LITnaqg$HY*WGmZ0pyvnl9%snS+wghN{H5 z_5{u_!MtNhEN_r&Ptv?XpNv*^cHxnZVjwqTMtv*X_B*z_vz zaWSVqP%cJh7BY1%M#ZE(1Z}CE0r=QT_%t`5lhx^65+y43qfRWPH2 zUvwBu*WhiO)>q2(Y1m7FaRSOUdbKy)gPC9~C$w(YI*FNdC@_wu7m!#<_1)K!~0l~ZItHMZ@E2g5pVoPMc(JnRdS)`s{sq)-@d=&j&jz=Ep zOhhs*(^0@8$0)nkZXcW$#tFLV>uPt;fAMLgt~`pxmU$dZo*9m7)6_3wRK}BV@1d6* zv)?N#fCXH{F3HOQ2kf_A>b2S6eZaC`TbDl|t_kPBs)FQb1 zq@LHz{d+!hSHA0Ay={D+1e~|gKQaT}s`qp=b3m2RKdaGYX(N;~jeVZ8mVMQ;!<_@Hv&PalSW0i&RC4ydPbR?C0!AQQ*xW$JAz!=h3wn5jJfmY{3QG zrz!VHVjQScbwNSCiEPTn$t-r`#+ST2`09GpCSoA0y$%b*I5*lx;x3G*h2?uxJ|!=@ zj3qiOeY&Sz`Ky=brY#?-u2i#AB~ke|LsItb+Z1Te)5AS}K*ZVYy86XgZEuTV-vB;o zXYA$LR*o}Ki^0xXLS|LVD`1dXPYXf@;vW@~8gI5(-FC0uxt27LRnpkqL}MpBtl9xEDYybHCnOI*pC^5u}q$!?kw32huD?hWVKtC5Wu z+w7o@y}nH>gS*TkGu%2vB>UJ|vg5=1$O-=NcWPQh)jq_A4#%yIUK2<-iOSr-b`!44 zJ;~ErwyTjE(<*H#die&domV>}FFx(Y=NjF!^GR+?&%n)%ez zmB{yz%iC_y4mxMVcCxq`aE%w7VSQA zCVE52yHaI0~yKSI|@UG&U-pN)E) z&FNUwGTc*{aycV7D&%R$4e8h;)kiGM;KvlM`A?VGKyW=nda$-+R$SNW553)$YgPNI zcBA8%guOUUy_s7*{W=mSMP3CpH zxj}DgMQUP?vd=IbzLF&8xz1nej(JnNUJD_YxHh=L0CoM?U`3J}?1Os3`zKUro6JJb zyHskgRNr04M%($uMwaKyoEL2C&@YL`!3Wnj#7m>oS$)sl8|P(rkY1w?hpOttMNC1y zT#pSnB2AWDM(|kD_j;XR7T8OH6KM2}(OxO7_tCb3CSM_QArifo z_-eXI(I8!+is-cI2x&)F@(tk}?l|8QnQ6G2k^Hm=KP}V4-a-qPNwU2jS)ZK$xbho$ z7aW%uaiKvrQ2ADeoS&QPwlv~J*uJ?)Ykq&O>UZ-*n-knb-xw!=;%=gh%bcNcv# ziaeDNSvavw=5bJ`Y|oZOuQem=8zwMDaUOk}47rZqgXB{G(3#q}=K{mi0l z>>t!VC|v~|Jd)))tvF~^ME57ptNlgnsfM4Cnr1tvlRZ9j%NN&49EM)QH+i^oKeIw! zB&BzvLSP*`*J4YWtghlqkL>X488}j{;83L5zVZ2O<<0a+_m?^D`I=jKj7PqL4eY}Fwtk$qNk#~!c&F&aH&w$`=BJ$C1vCrf)hzW#C zF}NnG*LAd5@A-x+q-J&tJ&(mQ?sj~66t?~vep|W&$#zH_DUn_@u1UN{7BBgN(|^#* z>%HeQlD8ukX2m+9{F)^DwLh$}oBe!1^~tzSBtszt*@UA#xbBht47u9Z>LcTNyMCha z2`aHml9%%k5eU1{s5zDJhQAh7CvPx|n?v_>9OP_udG4YU=I3(l?1gE$ZE_VJH&yMG zYVpB&>CXwyb#GJdaP^zo9;MU^k;nHE1`Fp78h!mYC6FlnVJ{|W= zPu*2fb10@V;vTKJvC%U~^Hkaa&lSZVLeZtWOinP%l#SG?AL6*l>Bg^1m7k#3b{45= zEveM#)b}eWm{;O4->USXq19Dl%FK@CPds3IJbT|c#V4WC2lTe(DyyONqNhl`IMKGz zY~Nm&7C(toufKihLri(&^4KEdkrTU#5yWcVgx4;ANlm7jJ(4i}eCTD}?wv2(D_?l# zR~&tylWSMgf1vfGNn|{w{z!Vxk~_9>#}Dj-R+T?{)x?J|d_tzY>y;?_LH(W*xYXd3 zBsaYTLr!b$B?*`37uxV@XvV8z>xPxmsH4@bcQ}4#9Bhc}j-$xeTc=HBIiGpDJ&f-4 zNN>ZsL8aPt@mG?$O$A0ecTgQS7bI7tc(&$`50}KYTC?vXUpB5>9+7n2tY{#&mUzq( zk2JjrTxf+c%p#}{lsvK`X41@0xybbjCDAEpco!<3rBvgNoXszVZV70&W@ zXY$=~gD2)SRjsaqIVE&TT)_m~hT`AjyaULE=v#RI3fy$o1ca`K7PDR*6Hw>)a@N#GONhzIFAG%ciD;$KOj_~MO8lG^Y;eBhCipX*d;WJ}_ja7cghMYqskFJz;14F0 zUa5|O{sM^{iC&rk-M7-H*t3dw^EE83<9VZ8Q9F+>d2?0`t??$p@?e8=F%BI|@x1HM zuf@sM8nqv&*CyYEy-rUquiMG< zw6|ThG{z?r;R_2wZguE|b$G8XPn($i+OGIjrm5%Q*v`G^xT>3ui_OSerH+;icGNAF zXd-7rj?hT5t8j*gQ%6Oy-Qu9h`O zrJ-&(r|cprw`kd0^tqk7DkZfpJSAl9+`;!*8Q)%=jbJoPNXys0OdjkdQOueM$LhVe{ zjoR0H&a%@QgCF`Q=Y=p2_~%DrBlqCmy6pNKe)k3Wl9TG;JHA1DeGe5LzvFe>`l+nI z>K%7DcC7=F&mI!P1rMBPO>|{f1r2=7XiJW-4vxVZ+&M?Y3Orx<+WmTo- za-(Vn;c)#*j%BFhC(6?Rv8_|*>KV)zl`wM6^!|G z8~xM|=YnEXiX-o>#Rh!Yi+K`QD~ozl(OUl0&n%r{8d;`=2wUuY(tqdE>VAVZF&K6& zCZOOnqV@b`s8bCuSK?!=c#~ZJ%5Cg#G#WO&&VJjMFt}dcTm&wzB3$n@s`_m0p*NYb z)NBtbo>*uD$A+K4l*GRCQI>)&DM;Qd(8nxLAU{?ol^u-rVg@PUb~5P)YZ^7^(lGLE zjEXz+kzF#hT_F;A>qO^s++|mW2T^_;*IuQk-0)qJ+I6TDKOcrdhHzYPMw~<2V=~}I z%9JDSF!F|W>1OmDC+z2?+Fs1comkW5*44qUIA;?5(ML|SztLJ@b^R;9_WZg2#3IJ@ z6O8*Ngq_`6xKfDf%dpw=kgaOI)T1(}<@a4vO8I1ACPaf(`}x|K`sL|K89vpQp&n1d z(N5FYRiz9g*d-WN6&qkFO*UUfSSMJK&*?ee{+usG2(&}KL|h0XzkCDf%sNK*+j{Q} zlq=aD{vY6W?MWRdlj&Hu4~#7ju_-+h$foJ{&{@<_+Y@&!@Ifz=rDGS)EK9EeQ6PxswH9}m{knHSmIh(a>mb(SRa2!u+JDw3 z4E^mB?DK751uu?dL^l-BKVnzAA8hu|j%{CIK)$t}v0d|S-Sc7+qkC@e!o}^DhsqT5 zsYknyoMo(z#}i!cYBzLDII3=Z_#nNU_~r%e`8V^@?88^M1P}O@*&dG6tnA$U6}{q& z70U@(5>ZBEep)~+#jK&L)TDtNZZu*m(@H|5T%@->Hy&n!i@zEsXC zoV%)!X>|!(-CkR`!T+$5bHul}Yj>;H57g#4SuBl3?%uEG97c;x7*D`M$m+1JbFP!? z!q?IFEqhe!n0^qIJK%A<-6NiTYJHdN1>3kyF74%Ibw?iFQ{`2owL5MfX!23EYTjK= zo*VrvgwbvSe{qMg)XQj#2g=866YOyvrerex)U5V9D(sJ$`mj9WLF5#bAzk=Ee?lQG z`GkHoAg{%)-b{-wRAL;SH?7A4)7b^9@ zb%6%XfS@&z%4hTOu$hz|RmQRmc+Yai-3x=5Yvoj@n99Za_~diMw$ljfGR1TmOx}6g z3I@MdP7qAAj)?o;*>bJ}BKCN(l~a;eZJB!=0Bk22>z`C$?13IK;e(7ZvEbiuyRN1bE0ZFZ!=#?2hYrE01_sB(5 z<$l|Hlxjv7i^ph5CwBMjJrA9M*r8#LIGUi`id8mN`(ol|()Lou(X6*WmF)MygZ4`ACDY4y@j7T& z*@=P|YU@zp*NH{Dzbw->dt1|1d;Ns_n=|+YoLe{Yo~j-4Q*y#xN%;Dt`~0`?8ATlv z`!?XnWKDex&Edo+8;h64NzXssl8HusyqWun_>tYP@eRS|HS$^X+(m)O{`n!X+^KI; zwuMqP`sD4-xP7gslESy&p+zbm{$hF6^li0)3YnU@q>^Q_U_tzfj~?%{uG&Aqp;y0{ z(edHIA&Pyu0`2P2!HVo@k(Q=7T`P-@E2gFBvS$s8oBdq8UzFQy+AFW>VR`aVWtv{E zQc^^dO+(-H*XTWmkw!ckskr0LjlJtS=xJzq|7%96MIL?QlH|RbQWMqR=5OU3esG}Z zq*U-Tf&O6jrO;;eU93iqm)_FW)M;J@#%Xxo0=d#d^UIn)R*+~FsA!j^S*~&oo-;gP zCLUGu()UUB%<8lBCRoo#H+tWm(qi>E&pKBRca_0C(Xg^$D|PiFqbion**cEw8Ne(QTCY zz_Da0?=rzUKU?vQ7z^bFnAAgM(CtRo_s>M->Juuz75P3bC09 zl_G|a+!6d?ov&uQF0bfg+*Cls7gCl$-A-vSF-%|g+PaO~*+;GVb#PLMv}66GH2+T@ z<8nXkDCI6DFJXE&*!2%ghSe#j`9G>d8sVGa%H#_i`u6OY9&XQBi?Z#Jkv#+W$Nt;+ zdV|N6(IkIz)_}q{ScP(ivkf!D=J5=6W5>(a^i67p=P=Hd=mxr}V0$3)GTDhm%fJ&X zkZPUL%_zmg-E`v~zh^x&{X5x?d^z^|jDs*;SUI#?jO6*=#~=sQd%hk&vo7;2gCq+j z1$vxl4jZ7l+{GCjzJlg(>drhGtk0?JN8L_jw=`G1Y)S2}aioy74;GspU2(|oF;(%# zLCef8yA}7mrd2FFSk$#MkQnY+ZK`u;y5z2?HcetePUFp8?4FFlQ*#f6lKii#VUGCG zNcRyB_q4hed@2r)Lt93a^j+JPc_v|H=G*n;bW^ipHpw+4w)l?un}_MHS<1U8dA?b+ z$}Z}f_0M+RoVGCVL_8vgn%4BdBeMKcT) z{WKn~Zfco_L^K5LAHY1@b@49CeMb(~*EaF12g0Hw#N7KiVP1Zq9CqV3n)Va6FLj9w ziW5}2@`yIRUas9=b%#{(*MS|itCISP>-rksOZeS+@M`7lh;Yx4RmA23qh9w!s41KW?W-BdQ9bx-F=U|{w&dwv z3|%^F;P@jW`?VfLtgMB&hE1{kH?t+~S` z-S|V*p}^c?YovU~&Km~!5X`62JUGJ$?t{yH2Ti&TJK>m>JfopX%7e9;U;l8_k$H$E z{p-2c+T#`lpphyF)3w}k{xszc<;!H|#EUA}eHCq3|M*l(sVX+65rz$8iUY ziuY-Kq)>`p=(pvfW#^~pQ2To9Qon87X!!Jy>of;$dv;0cM1FfM%<^So-wCBuZ|E(W z7=I&1lJJJFe0j@a^oMXMdG(hKojPygPcVbfe$7!V@yp4(%6E|aneIvV&g`BR+wJB* zEdV+Xb8cFpj2lyHJoD~&#O-sqmGgMJdNX(DK;J%IfcjD-7}Mk`P}=;>smTc4hLB#0 zH_EHilD4JzGL9t2YuE42OwEDz8fLjZx`Lk?fG&Zp5L7CKE`3UmQijT3Qngy z8enZj(90$I+SJ9fQwf7>4Z|?k^;17`j`$bsO5h8st`^O0%6mU3?|T?Q*uSzS^y#&* z6)h;cVrM-!t76&6olTMQ7ulb0uvH?FyH*V*6%v?CiEFo(q1N}*4y3r&x)AUP|E<-w z&z_JuM~DV;90QMQYX}}S`)GXIgqNr{W-G(gE4vSz3i^WE%R;x4=kMeyAg9AtEGB46 zqskPIPNaW`iw@V-(&A<&n4FqRUP@G)$80|$`-uEWKbz5xnqLv9UD#8e3hlFQl3^?C zRz;jk4t3j2`Sf0^f!HqO zRl~IoSK_v%D*)2z!&|2X=*ccZlZIP(+zbP z{Jo|8dpJZuWa!{>qecv8)29wkLi8PLSC6gNGA4vALU6KP_(yM&TG{`Q_H3KhN5gvZ zmf5h#0tqohD_=a#fM*v=I?bo)Xq8=CT_@KbJ&{z1MaVB-aj9u}@QNap?aV{d(%>Dq zHL#h2d36b6*_4Gm`%s#HTvsY`XGIgGIHt1Lh(NoPaEY;DD`uy1IWCDvy<+mD9Uf}* z;N};6)4)u|v0lZ0PQ5K)N>UElHdtZ&$%pvS7WWCHyBdHgX z*cX@~nBXAwYwjUo!6^01%sbfNT=jFz8)?D3dFR~kwg>w^aE~-=^;f>g40(AZo~c|! z@8z4fQUzaGOO?NJT`I+K1zTT7hLA5@3`b_XuO;bj=RRR>Cqx9* zS%pTn*|GL=&#~37GC#m?Ds$HAWz!9{(-arhHrffdi~6fl7riP~{AA9hKJJ8S!LZi; z?whkJtH{(apR$Z#-b1#%6>BezPB-r$`BWZH-?SN#LR&YaO|Lr>kl%>p-H8rNu`I%+ z$8Pn)-gkf8a+R-`^w0HXv*sgo^N||JhzQ;|<^Awo&(3pcsC_HV<0-p>TN`K+T`M11 zOz0WVW?XUGfvMx`KCK&MdxBfE)%mg=oM4LYOjfXLbH4hu=f1GZkE@aXkv!YJURVR^VL2{ZpyUbRF;W3H|nx!=37GjcsW(+}TC z)ZM_NC)9DyyO!|oSq?sio~3vhyrC`R)$mpov4b>V4{RU7K6zscvs?8@=b3kwwW8y^ zu^C&=@{+!7;Uw_-k;Dzfow%;Uys+I;Rq#Uv@m$Klgv#BV^Vw$dZ|$@PXFL<{z(sMr z3UzebU7v8V?Xc~$5!C#rlTDeMsJ`6jMvR^OccT zl?C}6w?Dyf?CVW-=Oyjiid$3kXg_gf%i^Z#eg088Jf*%}?YWv)xH)mHXf~ycNV6&9 znI2%XDx+7KJ2%56ujTsLR2;zgVk@6V$lKX-zn&XF+J+*x@RFjqR_#`gxfDN>tkOcO z-Oos+yy{pMwP)|$le5K3c$%Gl-5(birA5UM9cv^q!vi)w+P-M8;3@CRTF%#2T!Iop zo2N>LCu|`q(BT!Ed8&hQPeOQyzKE`*w{Nj398zXAmy%5MU~z1}9QHLVxyH3ppH0(uEg`2&;#{yTRW_DkZkWkQ znY9NOTSeqW)$ro@>oT);=sgqr16K~So=s|!8;m@jt-{AF5x==8qLHr!$t<5`7yYu3 zyg$vGxwTNc{ZebL=?!e5N(EULlXrzVDL?t^GIEm&W5BSFGBHo04@u-5 z5G}aZ7&~aEZ?@v(FK|oi>y&#I!rU%1b=|KQHJ*dE;?_KU$5%n5Uy*-R@m--X|4u7-FUIU=EbJFC_Xoe zr+(->Un+9wn9%cJHsczjeSG)z!B?3XJ|4mFJ(n=H5`+*Y#qX;(my=xQ6I2R?{vp4h zX`as}IZmPHnMzzUuV8y7uA`j`#vwoM+C$H-49CeH@#@K9d-h74((I=wej`l@qk7D% zAdbhJKIcxM7O_*~ZZPYg2g4F5vAK;W#8i0+Z_IjOZNv*6SIL=I3=xPDw&xb}7O}NF zd83`d%>AJ_eEgvo?9bN+2PgDKZ>l222b1LencoP(t`U|kWO{5c4sF}wqsF(%+|9D= zbDbGQFz42=Jjrg6NEKd0iAgh!DMq8J`CM(w^U010is8+m>Ae`wJ|v3`qsF3Nxm*uX z42mtpyf(FLNpex%KxIEd6UJk^Yhm$y{`lzU@s%u&_!)#+$7CnUK(wRSPcx%+W32z{-hOa>yGubtMgWM;H2{yLaj$#UHin!7Ol4DVRN z0!81JrCz+n<1xESWV$(0g-|EID|ihgdrcQqmpY>w+4P=yY(~4`p?g+EOj+jHC$Bom zF((({i7?HYS%L=A^Wqdu%izo?+G6rLb8fm8yGj{{b)_opT?u=M@NspQ!mBuDuuk8Nbi^U0Xa zUgfvZkrolYf^nR$$C{W2ulCnCwV)UY)sKXJroFvYEbFosb&1h_X#=d9Ie!!W0-uZc z`&JWSX(Tw=KZ6PR&L~g5XsQiHZes zil9`15ut)Q%mI17ZD2y3_Y%atM5aRW7mDhK7IBw__i$zpXo*i@fMz?ahP)uI{RnO2 z4d4B*A+GI0a3pg!P4)O%WrMl2bv&oCK)bihRL8f|x-|EB^i>(VTGUczw#2>&X3iCunttW9 z`yD*nnmWZxvo;FGy}#y~J1mrHe4b&JNGZRr{BUqidM%f7yd-Gu;IqZGlG11cS??aU z(#_e}Ja6}R?!mmbHVv^TpIYi#{7I^w@73q**>M%dtXYqz@fUs4NE*`@FQa%vF8x+R zo2RDzjZ#X;F50$1LxG+IWE1o6 zmswjru`3s44QR6ON;2y`z243}sWmx#+uXw4aqqECUNfn-z60uoxI26+1OBcQz7w9i zg>MJo+CXgR#3CKsLtSqX6$ok%f_kfxBUUygx12mWkm1v`uVGw8q)bPy$I=#)69T&r z?-MNAD#9B~VGEMylJrEW9m#7hGdM}O(?VR4s<{})tbkAF%E1lw2^r|G-zd(^&Dd$T z(Jfx-UhtEIS$uKBvL}@1(UmOogiVP}^U_p%MG~X70Ig*% z)*$HTUyU9hQw3{5>6!d3(_X!gY@>-t1f$X@KiTiI9)C;py_sHIC#9}MuP3+}+Eer` zI{NKHk3EC$Bxa)rx!t|(i&}`*+~A_vkb$K{_l}#dc>A2IU>VQz`DSmj&-h%sn&(dc zigljf4xPHZw?D2k=)(3Lp}EV+r`{KtebC@N;9H(|$gVa$Y^qOON%zND3>dYYjzaA; zAQF_x&^?CbhwkAVl>Bo$aZ%WYRjou;6_MDp3%WM54|#mS(SuSvu0F|gg{KFK$G0GU z%2TflIpSS$`m|El10J&_pj(Zo6Eu72H?CST5ym+dN=Fu$pXI6c6&_T*N@blg^Lsju zw(V;>`g+IxsKcB$YGpB6SWZH&XUx7R+`(S{@~{jRoLuqlQFyMAY-%UhYA@$aEA$Mj zzLdA>DUp373qGIo^?iHpw|EcuVdtphAv>bvx2w%6--~>1f_Bm=&F`aFaHK7+wBhci z!Dd-b>1A%J`)k}%dS-Mk8=5Uk zLxgrb;aX7|o_T-DTwutTP%^8EW78gN*6)t@Sx388MubJ5or6?}2rpy#;@3cc`nQaLLK9@~&c3 zd#G#PW)NF8-D2N=id!$%Q{=~`>@Km5ARj6BVC(zOw(X(Dco3tHF)vQY(5aF-cpJSQ zUkSl%Jju3gpkw%}FTuDx*OqbFy`*T}MB0Ltb^fuzyu;|hqUPw<+{fBn3-|Mdg}sIlGWI=NmQFJu+A8t|rAy3bdwi+dSwVx51&3+U z>E)m0q`W;=KeJ)aFF$x5x(pVmv7$S}RD<7Dz7OvjC3<^kGwzt+?57MKjlNGVXSBV+ zOy}EdY1O+fZK|^DKF?15%rFY$*YQ@1Mt)V#oiaQJ*Sq%u`h4$}bI@BY>Z=fpa-)2B zsK4<9o*BP7H18bXBSIU8)K-_?+ON+2YfTm5F( z#xn=3r5k9yjkupzM)=Fe?gq$9(rGoB#mk|ZB}1_aVeRdI^WlwGB?~=^XOTYX8uv(bDW#Y_!!Ns$(ePZ2C!g@*Ju%)7r$+O> zXWbC%uRzdttgYbEE_>D|pOL>)P`1)tobqJ~mpQxoLqB?}nqv_ex==owKb;-WTGw%9 zmew9Zn|fZdSq|~-A~tI+TP1_B<(>EwYo`Y%zX zfDpkiMI9Eyox`cmIxg4(*SSY(Wn(AO(UaTBdNHn_^r=I(kA;-M9?vDcc5ucyTI6G< zacJZt)5%t=~EVMG3h5o8sEZGFBpl%%6ffw$Yn<~w@Aj~a~!zA6iWc*DTcihlljAMlMK2+ri! zbiz|2ArSw=zot_Fiite(Yr3D;=I;|a`h2B4rul00@qml);IsQA;M`YnF{G8S$BfN! zM{xF<;7DcOunX+bDD$7S4vZDK|LX|sKs(u)pMgN;^3OLG%sGPpKXU{)TL?bt1|G@{ z&N;ur223}c=6C!EHx{N;7@kn%|3@Gol#|51%Sn2`==+aoyr_}CML>7^mi+=EV}KHr zlOvD+IO3=P@nEPf;QxPhEbT3wtlUv9_BQ-MlYiNxY=>N1)Co%ml!Fucl;!a=HfZ~y zK>m-xP;L3Ho-pLNwYig(CCb9`G}=`Va>CNy+2)Kpf7R>gV&i0aLP!OHIJ?`~q0kQe z{GF{WPoH-F9uD~IQOD6vPJcxH=z~C9oGnpK=riVmx(E{U?@@w*9&k}-Z0syiHuh)1 zeC_!)aNrZ3u(1Q~gfRyJ0(%tN$-?H0Iam`xONTN9;w?N-u9h}d)@Oj4TB98tZ0xQ0 z^|!GXH0|skRUJw(RNpgbFctz*@q+|W5L<^cKV-=F59JKvJsejMFAy!TypYtk>n zEHa@f+nR!Ctq2TiFeJ2y8iPQB5E>Hg71hu$!f2^zkfx?UXlkgSr9Ys#r52ZPXedG? z?tAxJc=aZMA;tF%UhcWabMO7|)N{_ec=qMnqk+j)CrERmN0SW@PBtQZrKmsH{0wh2J zBtQZrKmt`IfU1l>85kbPz_b!o8AYrQ`#1{ko3^%F3Hu_}hyAI9WjAJ2&c4VQ=uZ5| rgYwp7-@H0c3kQ)eL5sW!@)?EN?2!NokN^pg011!)36KB@{3QY(F;)9) literal 128048 zcmeHQ2Ygh;_P!Jm0U;n_M8s8E?tJGaGZ3=zZSqK5-}`ZQ_O>(M`ObIF-JP>} zG9)Yw`ffJ!|_7yHH;Y;uP-WJ6frXXXBhr? zzdWw5AJ?N>>(++*NV-04ahh@J2KlNaz2V#bs}Klj-6Aw&!hJAn;?sWcjKCN)I58|H zHhy#>-XLB&>402rv*T|_F?peA@KW-(c((ivM_7^l$12P7X&&Ro$#I(87y--VixiA~ zn}Ur1@iz^8jR!GPFUS`~+{&9(C;gz4jsQ9DcuYUL<5B|Ezcd*?2}85}0v zbbz-M7dnq|8{bAlpomez@{jx{BU4F43CCjY443h`yn|zp(`Z=QX!tkqzQO%&#{Qy3 zi(03Y#PP+7gu^xDA$j8pS*~yH-P$WNvucL@kA~UtFsY?9i1;c;>AXU&>~uKk@Z$6W zTfk|M>=XWUn7*Ihr>(x9#@X$s`{ZN?olo1`+C}Hlmh7VY)BCig{rWh&tx`Z zn0Y6ga5}&9?5KVZCx2{@_b*yiY(=;8g7(kZG$H2izc1MP8zk;l&X0*J6JK{mu)Y7{ z6IGP)8`|{-e(`YK0m&cR*O&IC+vi*N4H;$9{jF`%BkDNI{=ck(G1ETfJR0cd)BE-| z^Zu`Q4?iB!@r-RxYUJIAPCT>S9tWM=6#eTPXKed_d!uiKA5%};<23F3Cj*DRS@Hi) z*!v%sWBYnd<@N0GV^Krbg;YLmk9GVTKC?!_|LDS76ukBQ-Dj!(y{!2I=iBySzkUCx z&n~<}XFUVi(G2Z{QeBXf``LMt?;{U4E>*&dY&JR9pMF+->E|vZWhLpfu9SV<8I6Jy8{GHGx-%&o zi~hH7dmoQ@A0BR(W1Ra_odVF2CAoJJCb=Pa}`Ajz8|LRBn8M*Q)eoJw_+r$& z_Hl&Tl6$Uw-1WZK`@Bji(7f_K9;rS&+(D4}Ia*qQ#)H=D>G`OP`xE_)@5GB3x4Gw9 z(&*8xUDxhm?!|-;gbj*Fj0hVQ867nu$}nzmC!n&{aZ(vV$0tU|#RkVjCiYH*nl>bjsh@;ExnSJz2A{w@DzSJxduJw-TpyRP%TPVUc97hbj2{Up(puB)wHk0jT> zws};KxawZ}II2gI>RtOdivNKByY_L{`(E$!>PUg+mG|*D=EK7s1ex{7sRE4$)g$_I zbY9gX*?lh(jpuNygmL4WQ{9r@P%07WmalHP-j&(A_GMSM%+NlMx&G^xvjti|^tz=_ z(6#+cbxY8L*FKKwmXJo*K91u5l;-zqpLPB1^?sgCFVMX5J|35Rc(`ewS+^KB$ta)* zlB2Aj-(JtpOH{Y$pTp;O-Qs8|ar3lpS>H@W-B-6rJpL{JW>>d3@qH8FaQ)XUB@48E z=yi)<3qt;?G_w%%0f##L>@$mQI;f{jL zy2VkT@u0dzUuTe4bqn3kQAJ9L6YC)UyZd#F8I{fVP2T5!K!NsmcAs~Rlgn-mPwU2{ zcG7>8Uwn1rb)Oo&YG8JC<2dw3(apTQE-*X(y1L0bddlCk9pv#;_xSir9{YD}b^g-F zgh9|XLLs~7o|uFvi)!7w$RqXjgbuoZwxv4XuHV61TF2+#U7irsQ*Mc$M-W$Z+$D_t zVXp)CNsg8;=(znDpB*aSC%xW3Ks_!UW500u{`LN^_w#ySfzIoaS0#)9I4(T&+8zwV z?)a|PK8|=A*y-BGIj}nv@sQoo#6M@}?2eG{6Yr0;zxHw7-s^o{1r=ysc^?m_4-YqK zGoK$re0V5$g!=GM@aR*Z@u25@y;K{{~}k$e|8?i z@A=XBadY%^+T^~gB~K}ufCe8I*+9eV4c*;Y*&%Grc|CWEVt7`+pRRpg8 zx;Cjm>xW*~`VCg`*GE*>I)>mvuf;qKwQ)_+vjrbk}usLbc&8R6GH z&f9yvU&k#e(7f_K9?3pD+@#H{YuEemQ1D3c;i2G>TA=Zux>o-_IKS)KlDS-u6ptVF z@zu57&aegJdfC;rC5Nf@U;lM&T7lLNy{=D(qmqyxeUG)XkwW;!Fysmmo z#+wf8ITnW(*=(q;-5;CH0eGAABwXXswa@zR_jRoBw}^uHVN-`y`i`yA2>w7JBJrU*i@s-|G5tNv6nU8M-$I6fd=7B`=q(namwe{h6Fi08@iXKLTED&oN#7|EGWV;u7*fk^+IQe`VbE$={7EW0Y}MmtkR%aU$v(EE# z%coROr!^+LX`A_%q?Oy)1Ieq*-(U8hVDGx>3wbZYNE?6E-S(Rfyb79Vzjf6YM2CRK zuDZLxym|fbYDOsu+W(eU*AE{4rG0p0eJ)*pUEo3Wyq*{GuAaa8--|Yo596UUIj#4Y zF-0Coy5+0G0NlT1qkk?QfO$J6s{nYvPdBH(r=8osXBdP(Wl9-c+OO9ihN0#Wb(a6T zkzThPnR;zoXqMwNKj~*V@J9M{eUCb(_&+sGeZQcM^5D6?K=aD`c!c=ya0f-^=i;FS z8V{;lPT}wT8PJpbu3KizkxJuf-Qt`rk0%}U)hz{0ssB5B^}59mmpg{ff#~XictM{Z zgwT$Exb?oC_9@W%A?ubR_{_#I0%!e0VE^T(kq+yhXXdvK`z?@;@KlHQ&y&adbok#% z;OulbJW0jl`q$y`LZ!p08vd#`e533TARX5KK5%~P@RB9cLQi!#bg?|%r^DHKws60n zoerOUTAjD|b@_mea}c^+S2>QjDPz5V>I+TiAsIzzWV2X`Sbl_WU@ydL%Fe(*T&bXfm9I>55GNKr{A`t$5WE`40g$nEFVmu{ATsEJzVt6~4* zYOdh>uMaL_+@!55h?aPixH=vg-|1)kN9zd?j0VNEj#o~{|8pJiUo9;$Xlo61T>aes z%J?sSxwSF+$c_77`My2G3;KM1KgOx}yo~sue||Ot#|4#?FX%Y^ca7uedh`K2#`-1& zZ=A(?EPa1L8UOeCdv*U@`MHb)v`7(^uf)TD&aW96f6mnXUy|J5U!7luF;3Z^-ve5& z;1A*z0mF;jpyENiqkuahz{H*SDX1{_*)IKN0Zp zEkDuD{k?qt$vtNIc}8LPPhp>5esU!*?&b4OUZCo|Yx?}-JuLI$N`uvlhXy^W3KL6w%v-~{6=b!vUJNNhU`6u_7<>wiN-9LqW ze)-9jyttRoKY4+w_pa&lkN2?5i!1s3lNYFZ@0vdUcn{0GxKd&FPhp>5UT`8muj%to zexjZGd-?p6d(86l44;4U6YbpJ%jcimW0s$16n6g<_W9)}SMuUsKL6wes@}V%&p+P7 zGB2*=^G{x&>b-0F{Np_=^WsW{-9LqWetE%({Jf^mKlzDv?(gOEPwp|x&og}f$xpO% ze=nbZa*tVlo>AESQ`qO1pIpg{d-?p67pQvgnm+${56ir`lFvVRfvWee>GO~Gu*{1q z6?Xp=_W9)nC-U=}KL6w=+PS}%&p)}xEI-fi`6oZo&i%c7{>eRN`FTcR_fKJ;Uw(2W zFYe{@PhOzvy=(gX<2@|%;z~aM~UOxZi1*+b=rq4g#!!j?f{>e|YbAK1D;0MC6!!V$1t;?Jnm+&JC)&Bcm(M@B z$1Fe3@cAb{(a!z7eE!KjX8CzWVfRmApI?4*B`@yf^G{x&>b-0F{Np_=^WsWA|KtU# z-n*vHKi-L^|;#m0U>Zk1U_a2|(msEd?ef*pQcZ8p5y5HWf&*y_1 z%!4Ui7uoxDe06(clYd*jYg;`h*{|ED<9+$mgm3%&=d68w;A@ZfNPF$^gFS!Ue`e!O zd;FoakJJ6Wl6c$uYyJ7l=u^uh?ENh@{B(Su(eTv$qx(;ve=Cg-vo(B=YIuI9`A7HP z>zcpLYW&st-CpCfZm%9cUup58^H-nG8I8}ny)jz+L~8M(+o$9GhsFmTU)^4vZ?Mny zKlp3=zpBOu;AeY_rjI-kK0d%ilK^!U*I1%0s3hd$W%*Y!ciSGQN^n{J;8kD+w;@;5TWG{-Co^4y8ra~>-^Q{qsO0auby9Zz1R7x^FhZ~ zw^xrJ=!K1su3vh7)Ad2mmwNu!^O0^J@}GVGKd-KDw!d|Ib-o??eBU0w`!ef=`C~p! z-}TLOTtATAT6sp){ZA~!_1dv_0|xxQc+EzB`+om;^L%={1N%RA=jlCPwwbqW`efd&e6J~;OD754 z-o%+brv!GP`8dqdeG+b+^;rt;7ht}Z_~L@!U)n8xSan0&TNb^=-`D2T?!|H6&MUE) zcRBvccgND76RN$2@mNHIo`JQG+xT4gs(Is5e%o+;!}$5Rb6J zZ9Z)d*G;TZbegEPd++1R@0-tce06&d*7NT_susCWE~g8^$|ptGSl#+hvboABTNnx5QfjZ~p!E_lCW^rJJbn*V<2? z3=HSQ2Y<;2uHh@%s(5nMKLKWv!Z*Xv=Fi3ezYirJIN^Iz;wwT_JjDdXKLO@>g>U|5 zy5_GYDt~o8aKiT;$!DS4tH)2ptnuRNZ)*JI@tDss%x5q7++OmTf3Mi<$F%riFDdcT zkUazcAB25`H~7IdJ_sFO-CmV%9ZUiL!#<($gQu%}5Wv^wlght9^RUJ@ZGM_QJV5$@ z>oqhz)CZfdx_@-NK>XqQ?M%`G@L4nlKSn{{zy}U}6X#TanF{}0_cI+2;%^T6V8@?s zuNptLzQO(g_CECOBJ}M8=v!szbw?SGg7^XdaI2n+XII3B?k`=xpbxel>-wPM3;x;m zs(fq3puYn8%NI+2bK*DgITZTlp*|3vT+?qeUGbm8)b;VGil41-z*lro@f50mveh@5 zzo+^jDE@v>>@^H6en=m1{fefpZ|`9~w?VHpK64NCO6gY2q3eU5FV*~i zWj@mU&!KOkGxDFt2hukiU+9}{ugbR$rmheA{MG#I#|~)t)rCHckov%%gdXpQ9;^N_ zn9esP{)|9#IP$|s$Pc4Y|4f$k4=21SAM5t2{0lV8D)mqLt~a0$=cPV~w=kbd-~&Ux z?j-T$YW?HK6n%5stNJs*{0{Zc3FK#*59RMx$k%tv`bX&a>h`Mss>7PVKGHA3`%U28 z0QFBNsSn~srT(cy`UZaZnR@*L`#9_)yeXfEOUQS3%lbzUA8v)dYWVV|u$TC%`p02D zrPM!uOr5{MhJp|8qaOYi_#Q`nso}{VRqCIHrkbA`P(I@F`8|LIg1@oIZ*Qx75VO_# zg6sCG@e^n&^@QPPCaU}uG@qlG5Ait!^}HG%y1u11RpX_ON%p-Z{m+AdcL(T0LzNF) z$5*#ks1GXK7y~+8O2i$`e5@_ z^^albdZFfLKl32`UlDpx5B!Kz^^G3@-l{&-VJiQsY4H}oJk$p}{&aiQ__6g3_@pCE_Kv34%_ToXd+-DPk?|0IuNEJMq1QiJyb0nj>9MX4I=zW)OL#)IY?X2+EFpg{U*{tx{Fm!v>_({LQpJcJR)#t-tKON;z?F}4{ z@ePQdV2qEY>h`JpFpOUnKBNWd`0D)A`6ey0;|~UQ!~N5MufrtE}`wv7v$K;Q6F%AvZwk6`%_a-D)DC+A1i!*DHZ&{@xU_`d~$sZ`-p$p z_{&N0mW<>5F+akiI_x(L@XN**{Il)V`KH?UI^`2C@k|N&N#Qf;5%Ggq=!8A#LGm5a zGvceR-}ZcQe5O9Cd@u|>f2j5vhKKq_cmO}Lo3{FVn2o1Aq9Mjdq`g<<*VdT-S3zoi zJ`$wz!7z|NfNyfDhx$hL1Ft~MzHYcb=`qDW&5!h#@}(2!lU|TM{G|EcFjTygQ`LND z>l@iedQ9UWoo~e}e)~Vy|v3_P4>?#|^fNM_z5%+V#HRQ+oH0_+jqbBE`J(^wv{Pi%v7x#(o!9%k#!Os^P=f>@5Qh zy>?&*>H>wO!V`cw{*;X&+}f5i}gM9NizSl z?~h;9ue+3QZ08qzqVZE)xA$^nzh!S!rTG}fq<&?pu6=9=?02xOyC)r(bI+^dlVUr{ zb!>Wo>+_lNLU8i$H_qeX8y=mu_3;(3uMV3d?ZtI^VLs&7cC756UB?e44&Z0DKQ{OG z@22tTC0DJj6*`_b!+gHOeEi2wiGDtM0sn3O_bYnFOy@ejy1h5VzdUY>gXYuDoC^B} z!oC+YybsPcChKE=?a?!S?wc>I2y;RR3hFZ@@jk)ch~T ztNaizNq)3t@);PHY)%R5B0&0Kb3DzQ;#3`6LvVvob+vmiOy*YFi&?;z;0?jO~Ew!YbTc&HBo`T+aes`*IIuW9W~ z=!30qlVyJ9^EG|b@YV8ZfT_lht#7cutx5Vt`GNEp`e5e^RUhh@lrJT}0*QaPZe~7u zNsTv-d?LbtcQ*QF@3Df7dN2z}T%_yb*F^ngHUvK2QSm?v!z&9EBq9gK! zHNO`$eRDq#&djP8Mn`r5`dAD3?kD7jA*g@4$@)hW?I-J>Rmc~*y(<44Y^qZK80IGA z=MRye=U_g4FrP8VAHY}g&DD}UDeVQWqisEq{_AO0fxaC_ex8f`Tn+O%hxt@SzWzn# z5253$+pGA?)zuuQ*++O&e%^z8cNqE3vd_`iBz�t1mN@{6zJ-VZc5C`#A6x&!OJj zgZgJP@*VMk>fz16_aoGo<6-ZBtLlTFR&N`II{yGOK;?rVe3ziUyj{gp`pXyy{~*5_ z^-cH>`x(0h`WCP9S4>s)PH6aQ`d}D({HXp+-!)j}FLz)*sDEre6TXs9P0Sft`4&Ba z$@l5shbUfb`w(yZ5#SvOyk7?&-bFpE1km$<j{z{cK)#JE83lX zulT1W)oX~C2gyH7*Dumz@K@;ipyLbv+4d@YbL$)UYuooM@*nYgucmJv>I314`|A3w z=c6Vj>Ejb>eh?bI9Qneo*Hr&xt8cI)fO)77q`x??qRKZtzh=>U%Fh$=`NT|>-(r>vXy*L2E`7dW7~5Zo+|sG`P}t#MqXP%7G-gi1^&ISuEnSR#`#18K zq!pL$X|RKT{ldeo+uXRD|9s)+zr%{W#cQ=KvO4^YHR6-Rgl9hfbBj2+t-~EM#5{iFXxw! z>{K1(t^9f+1_Lq2X-&!I3Phb4IVa(ft<`a$i3~k!ALWP)FJgL}6cm857<|8$G ze+&5bV6*Cn`0TmSH;t+MiKu7J|8nQhy$Lf?8}}4l(~nMmHDVyw@zw1$R!qJ;JecOw!Hm}ItEu5Vw_TeQyiG zdEagyh_Vu2K0)Fu>PmY#>=h3`SY-NFhu855fqy)4|H`f6CutvteS|kJCGi%|ynXZZ zCYxP6W_#NU2TsivmBELh5_*Rnmat%+R`p3^y{TIOW`D=XOy@Bs^iLbZ_JiAFe z1^H)=^be1R|K5;%6ZPS*@zP)78Q$#JByh_XOhST^T=mya?Ss!uu)7U$GE;s4Mxvb$oSu6~4K;G1$jTBi_CR zKY%x1s`7y&9{3}OhnvB_iHL8FZ(N_BK7Z%~&!s+a=!5W3-$?JwiWZ+=k$M;2mYuZt zaR~a?j>&l2ig6pJ^GoL+_-5nbp+4C0r^bhasqtg$8|(=%-qqG;eJ=ztMJX}pxNg%Wj@Kodyd8DWAb{= z_NIsWz=^NqAJT8)8|5Qi9~C|rJ(!AbxMr`zRQ;2!zS;irP#-w-;fj2u=T}|tRsK3v zKDgK0DEOweW-8yZ)i?OFl}Y*qeGt$GL3*5xKG^wE&;J_VgdT4Se@_iyv$OKSHG=ZJ zjjzhTR;I!?x4zZPINxDD**E&@hNahFzJct!PV>sn{W*2=G7f26sy_nGHkDZ2&d_cJH082EwX;6DC+rS?N+oO(|@x4ix@4^MbT z+#XXRtaXL8e8R6icZU8pR=iZPVx0ulx&yEw?f;adN|$w~wqB3!UZ* zVJRI%;%xs1w*TH+JQsQB{$FY=a33$~H-0*qKUSs0kG=O#5Xa}2eX!2?8RBeKJ~#%k z6a`;@)n31F_RFiC9{;5k&Brj#)eG$wUStjIb(rD(M~>_fJ0IgZEEs$^EBU}*0=|nR zzG9ZNm&d~17pi?Qc4F6qyxi2SmzqyM#yc+y+WFb4B+(_f>55MdEf=@8E;pe}uQlS& zn+BJjzkG{Ws_-{P<%3bryj#JyGt=!Y8M|iqygD?WK=YOrQI%F5TMhdISP1yLMDmvp zR{0>-NPNX6;2R2iLt*cd{<9~IuCkBcd3eZegG1iq>&|@Hy5^~|yy(mY$J#tUmG4sG ztqm)emA}pqrtrZP%5;2ndpCS_;|mvOV?M4SCX?|d9Ei8#lE3c!DEaFGA3m0R;FTo4 z+$d$?%Qs7Vg(2+~uYH_YJpSb+ymF<`^cI8u!(GsaAEZ9;HP8pdn|M^}15e4FJ_zW8 z@KE32uR!yAi_fo%XD>8tboT7Ke4)jUjbhllJ68=_ww=FX;k`=i`SFJ(QL~=snm&jP z(1(LkA2`L^rrha+fIbKh^^MDVrnb4!ioX>8L;8>JcRI359JKV!Oj}Z<_-?s7MOXFLVb?y-!>Vtqj@D%7nd8rRVi#HB^ z;D19O&Z&G5n}F{Hi7#iccZ#%Ec&Klfe;w09ec;fCE9_kieYivFgIJGvD*}C30{+HG z{)$83!!ME#ybKp72AbsG_2LXNH&CgV-W{G7@=O8)Zd;6n+?2M&F+@rAzG_UigZ`u%ygB0IV~cZYca`KZEzr^0qu zz1u7gJ#UtH`y1^mx%-7w{KRriAAqNr20Tf>tH3{$kD%{71@kX~`NvuLVV`)Y>Y>`j zM$I(aTKOY|A78L-$HLo3n(M55@f@!m@!ra@L%K8FUgF=-M|t)2j|8&8$frfil>Tz~ zM{U?DF^3+Ddo+D`8SzjW@lX-@ ziQ@ZJ@a;D6t%Q}|hVi2HBR}J@)y-B`{#(k29j+7@+%blYwesC;UZd!qIjatHn-9=e zPWX~OkiDe0RG)S^zxTfHYc*&4QEz{GyxGzZq61mSRpt>uJu~}S&e*phpb;kbLczOf4V-Lhu$U2 z{3lGA?|2&Woo+AX*FFp1zAK>1#a&O}c_D4&^L_e!`RgvJZz+7l!Qj{87Dd>6cze-Z z#l{+Q&DIiMvFOaL^?#T-f|0%cu$StCbkqk!k*_wQo@G+M#1!NQ;_p=Ck7}Ae==kdP zQa*aS%A1SV+}(q%m+>a}yTe{v@<%VWNb;B8&OY7o)aXDq4t$uha(;{0CmGZ8{czMT zYf-;cM!nrr*4tdu2i^tsUuV>RJHX#Exzh&$eK0-LHx7Mcq;HpyAD)7~jX?fLK)#p> zeItG74tt;1>TilS(t8i}fkPjz(Bo|MjYHqe|ET)LzXhLHg3pcB`a*1j{+3krL6AP| zQu7gqJ_zUogFXo8gQ@ES^v%W>`exg!>sx8mXQXdi)p}D@0iSO|zVlEYIP{?@=>t8# zQoPagE%bq}NB#E}^4;&KAJ+okA5mZ4tlBF))He=&L;aJZK5*!Rxk~%|fIeTifX^33 zBHlJYA7-fhDQqYYKhf&<7LxAfOLsHu~VMm*wXLvCy;i(1$9} zhi-Ttr+7nsCWyb&p>MaU`X-=nHonj|+g@GYzBn{(-ouQZUma$zT{YJ{P+|`1r9gJ# zLTJUaw>%@{^M%36PkyS`wl$l$qhr~hm!t5QmAfqDE9&3y<&aTNz>ZN^&oP*_{N#<5>xU=tHEUhf@A`Ez*X@0|qX;&e(R>2g zy4ZSuZLhT+*YP(u-VxEY{JK}gcPECtT=sT$sfARc;TRZli&ksv` z`R&qPd7ZS$^80YtH1@&5L3cj&!cy_a-@CqQe)uUlUn5NJ4-YpnAF)!-M;!bjaEWO7 zthlwt{e2dVnl3zEy_*-k0;|zhnNjRQ?K@5BM$S0N-HA zXRh0;#*b5r7wT&mZc5z3%``b5@hte=TH-5K1K)SSw_6Y&)_KXJ<$R9KG2u^<3V-VM z5#FBz@6W-9Ajt==^GoHQlPP?2ug8G@ZTr9v6>lQ|{xFQA8orGr|HKy9D{-Lp!PKWS zMxK9aj!E;ci1~Byce^}ac&HBo`XE&QWUFt`2g59F@fp8|mHd!d&Em%q$vfA>Y^KGB zm%u~!?*)NhRkG?sN$_v8n}3Km@D2CLQ6C6jPxTG<__LK({Jn$w3}DK9!h5j0tawB} zuOq^Qe>nUjG(D5oaXke8wlp>WbK(c|!L~#C*+KCDdkho&vyW5wW^`gzv)X6eW6o0W z4lz$ye8%@*;~rri>I3Pmv?GgtlfL`Ge;<%4U+75sv}+4d@YbL-puzOU{*Z}-3Z zvC0Xe>?1k2) z*9*(9&ib=I#+yKes=FSn@<*%PfZ@)clTO;s{G$yly9|^hg1r7 zKM@z@#-~$LKKcBo-8H87sx>6wDc=6|Me*S=lSGvw_kZ?Um;Subc1|>;elrAQzQ57>q`Co zjaWlzFaKKFE9dE&1U-mi$#bfG_<8pv@t^ds+f4mrvZ!{U%3GngP7zD$9bA0NO=I}3 z*Y9e*V&`J9vhA~*8b;0HgO9dq_{x-M9$GQ+r>8sj7dpPWz2o}#m^H)x{%Z&u@l}n# z7l^&mUZa`$YWdEkAO0yy+U4wHp8Y)T(aEo4K52=jzlQI#uy+IOmGRdR^9eOCkBu0% z{o+6|Ah7h&s;7JM@=weh`B%+`+`MVwd$+Z&C1xBhS8ZzBhxppa{V$$;yo+d;l@CTk z^Zu;(8V_mqmM}}-Gppeb_3Zl8tY5G3_Y2F8k@g!7w(XOK+drB-5%a5K8kNC^Pb@w( zCwv!4e8r#AUfvJ(TJt#t|G47Ki4{kmsaj!?*ySiSv%!QXdEE|8A6R_H8qT)8{rkB_ zi^a1Fe*>5+D<6!O=G$5EHTp8$-h{Y`dxJaJ^I?D0|7i6M(a*vDdS-3#H$?K6zXCq2 zmV6K;CBC8~@O@0$%ePB=<#nRbKM)>g4hpI}?%m$ueB8~KR+s&(o8VdF&5gIh-_Y*t z+pK&zHR=j{b$cV`KDlu4ESitM=`Z6=FvJ`9EBb=JPk;~4N&#B^Ezk6d#YSyz7y4$~ ztLq!-ckREN2OfQ;I;)I)^iuhrK|cgEHvfR0_pev;-uS!y8R6F*@w8Oahp(-C(NpXI zo>th5ak2XCrrvq-rg-K;K1zfCI?DNrAzAg&@Hd~c@&nfCx%t}7${+KDjKgl`J0&yn zaU-tVOZ?k2r^SLBlkQ>fA)nq*{qD_erZ+b={t8*orHN9|J6T`6%n^U$QN+hi#D~29 zWOj3#q=WZ=*v!07#*3hQYVCJ_?z7y=&(kIDuI?=IbVmLzF6~bl%?4ZfI!fsH>h_Yq zD&NwmeB~+ko3Bd$3-Z6re@i&=gZM!D&;s~Q0>0p%0RMnz6q~rcL#qR0x|^2$&H=2O zmEVR-`;8H1dn^BiW9(o{tb7+HommnvK(V_93^`ylG=vv@wR{E_k)>-}}o>DW&b&D&)?E9iOQB;pPFC6GUezx*5I3(^N2 zU)^5HN6zj&M^+z0zdO={ok6@=_P9owk>Ia2c7Oh}Nb-TNmHCseMSXo6>^+U=8+rUd zHtdCiPnU?jANWjUYaW|tbZeGqQvHW~Bslmx6?!hu&ru&3^g(#2Zyfq&YV`u234WYF z{^*Q+u?_l0diH8oeJDPerMEckxILh&Sx40e>F=EN!Gu0|s&5?n#%dt`zLE70Z(}{5 zq=}o5zlR}T8}C^4h3Ei$`vYIpucT+-X9VjKIekElYWN*3=KvG>Af!J=Fz5runfd^I zv+;$#+4k!CwiEen81h>~6L{Yh=mAM?BH|KbVgpQ(ZVnG)M&|JRM^ zm(lH2*WbJ2Lh1jyj}1Y8(kVp2H zY3NT1*5-2*{YstDucYIv+pDZscmJMn8_hl)@9F6Gsf~Ug;LCw8e*pboqtTxejsCCk zvj2-8L4ObIyD&^!p5w{a?HS_%KA`$43C)3h1wE2YY9upHBCW>OViWTJe{A z{j-XnAbh8yzwWS#CoiG;$M)Y|^cyx({l!nH{I%x;e)AeCpM`F(8b1!E`ZLHR{?TK@Ob`de*%gFk~v9|ZJ)Lm!AAT-OKU7wH4>kMu#$M~Zz$ zZ^`D?^?Iqkf$fz zyUi7w0a>_T1NILuz}YOQ~+GlM;X_0N+Q zZ@F++qXt*ttJ_QK181Bm{mm0=>+9>E-RnRy`t4TH`e#n-UC|%GAIJL9QdmCvy* z0qdWs--iqIAL#bd`sY~*zZLtc-96@1tbcCYe#Nt~$)V`~qVjgJsyerzBN67^{zi+z3X#W zKROcqUw&BsJOcVg_O^q)snCZq=#Qw1euGQ7(+3WHxI&M!(Kil#V=rNS?$21COZv6~ z{cJ&4|NJ=Ck4}O9l0Mu8J^Mr2%MWY)5zq$#eK4U90{UR{f%L)dzX86`H``ua-^ytH zY1D5=>vQXY&!y1+<)J=s=z~f6upj+-8xe2RZx90huEqM?H^7IVupVCax2Eu4(fp}ukGo2UA~p%12}4`L?t;T-xSRw3RtK_6CvzwubVukk@>_zDm8O+er5d<=as zp$`K3VCM_chj!3k;46~V_0OUK^r18QBZ?v3qM&cG|9OjWq2KOal@A>HX5$Nev+dRO zjr3c7?$nUIk9;J@4mMfoe|MO&9}xZT2P4=#3t#lVy9ToDRz61mds>(|3;y{I{#mE$ zyO@vpuf}@y?&x>^68+A>=pTL({llHDd_0p^ll{0tx0m=w_Mb&Qm1Bc#u=IiY->Wm( z4~YKv@a9bRhob-8{d*MawfNYI{Q4&1r;X;%?dZ>*h5qa@=yzTu``<-5^uG_1{lh}X zSGSk^CGC%4GQXk!-TikUhFJVS|9kQ<^Dc`I=zkBH#1oj*dq^EXz13HacCK3bm@Usq4nURQ|L|L$k^@5%Z3 zu~?=5y}H>Ad^m&OL!jRmn2FyPSS{`4-C^$zWKa`&+pH|WIJt-@HjVM#$gOSJv zWyLpbz!mn=^8@&0_+``oJ{!N^Fk0mU;cNE~cgOPq*}EFg50z0L%tpSdh?86`QVB(b$oSuDIdxCG+GuY%(3KjVy4z9#pSF9|ZJ)Lmve6!F{~^eGFY6pl>$5&^Oy&UEdzi`pE;4|Duuqq`qzAsDHSJ z`XHbW;uYw_&3Jz8iFoUX=W+e{Q2xG$%?JKG@cjk#CFSGEs3)mDB|W45WEvy>F#qbP zf09KG`e6V52-!=YA85}He3lw-Hh)q7aH{9P2T>jP#-X0q?afBtV1IxKeXxIz2KvCE z4_wy=JP!!^ya4uc(ua3dec%*tRR5OJ`t53<{wafcSjQLoX4|Xlo4UT+VRljaaSY=N zrT;5{l|;W!RrLE*NB_*Hvj0mwiGG~)R7_qd_#FzhQ2``Q5SXMwkMy#4psoU8=+ zPzC*f%Yg4$iLW4g(f`F&{{)zw6}~ylO^Uz#%s0S?#Y#WmXxJOhj!1hs;c3|yzg+Hj z&7*!IaWnk)Df&Au!C!ZSzbDYYSRC`Y6Z4r4KDR_aqi(MjKTPqbs|}NQxv&j2UsL

*MpJ|LdRFPSLK^wAl?Mo z_wLp3c7N|8=h=k%o%j#H7x679zOCcZh;M!6dH$GR+E8{J^QV5x2cQqn%lUE-^+7-% zxaObi^v%}CTeJ3moqdSa&f5QFbT(ywEBe1&Ew12~HNSb_-^_U?_$I(NVeQXZABfMM z>YKc-;m>4$EA@Z1V7FwQ&reb8{;d5Dfi;KY86#h&NY%6Z*iR4+8oi zpbx+?5PWb_JmCER6Z&S`tMaXf*&wTZ#yw2O`%Kn;pU_7<)dvB6koV1^-=C=QCa6F3 z0_Lyq!TtM(D!xur!PEVFf*$G{=1;#*pQAo-#2=@A!_vr~8GOo`U)!Ufiux^8{su91 zKK)fbJNjOs$6EZT{&cl5DIb%4&<7RoEc$Q~`uiL3rT$QgZ!12sX%9ai9&u&@OEQ|HeiF^Ft{mKz&{!Qzj*SxUu+2`t- zwB8l{1pFCo{phW7y^^5y`qq8-EihZJ`{b7cRb2M^XIf8O66^2EmbLm>q@8K;Y>{05 zjCl+<+sgINn8(9rrBc@VXQA6m>z^xs^zDRATe`7}SpR&Vv226im>On%ZGA4Scct}# z+WKccSzZ4Odls7u>+XGLS=)Jx)(bYqdUk8Rj%ahPvOagP*=BE+^{(mhX0%*C%5{8o zdue@O)pm~_8aAs9+pVsDw%6xkeW1M_9_ycZC*b=j);Bh|+ImdaXmeDLO&`TP*hA`8 zT5r}yt|ylBb;hs=xjq;3e~`(3pD<4A;k#jdpzN>9>_4IP&$M2f)(1vUY+7#p6J70o z!_dBWx#qU0WUj&b=XVA*F4g(*VCd*@;)C?RG2EP=b$z+}bFg-by^GCqUoD7vb>v)o z{qyHoPrMoH8`og{-7=NGTwf1Q>z~o@WA_`9y|n&0{Pl0AKmXeUtRB`s@9xwfEZ}?k zI}eGb?Ej+m`_bkrd$O#DcmHnb!@$>x_0OlX_7nA3XfAvJ&dT+~q6*^e zWvs`g_0NCF{u!aIe-=8vy1leMu-<@$EiSI2zuP;Al|j5&{5BHUQ1Dldn-CuatshMW zzIUkWwZW&cX1mLe#P{3WjalbA7n*BluW%fSdeoePc%y#3Zs6}>tlzhe&siTB^nrV* zZ@j8ed9y1{`pVzSCt07mw+#=*KU;6Ugucz$ur&j5>|#t zf2H*{p${@1oiQf#0r8os570LoU+9}{udZ*eU_CMQOSI6|=ca4xT_yi=)CUHAfPIMQ)|1$LLNAO_**26>JY3$XZ${&#o3wyfvSFP_|h$$U-y?={U}#K#wi z59_+&3*~%_g$(^Jg7m&7`rkKN`FV_Z8~x5p(C-|N{`X)jU-#!azPi2SFFBt;bEKtT z=zn+iHXT;}L;t%m#5`^B0rNSPz^MQIOW^yWjBl|X`N_JjyPpTwXQ}9apJe5?o+4$s z?0-LVqt*Y;C!pVXA^M$VJmL2c(jraSA3^x`Vv4yzU(pQwod`b6mV6M%pPZft zzH21=^^d?_dEEr?4Sz>`IsD~%%xsBxI|F^Mgx_!Y0{ZBoJ_zW8@KE2N(+wE)$I$Z- z^|ziv{*d~%L@Xa`^}pNwI485}!$kP2&qDJB#2fX)d8iKp`e4`F61Qyi4S6Tn?r%l? z!zq7Lf97_zUx{n_%h$o))yV(WeZxn>U%@8yfkPiS^g%!$gsu;if292lNRMrMb$yFL zeMa^Bapb?%$ba8y`OZUq5YPvC-z@q)4fXbJ^dr;nH#~^?uQmGTs9%cuhrdUCIT-d* zzdk*GlfQJmK>T6;hT)+;aOi`WB<&T@GkKjX`IX|$L_JH-cT~SWj()bTsP9h$-)B+J z>-J`&Z?L}s=>vy82(+)*#Q0nLkNg*Owutr8A8#D1EML%g+&DQ-!~TF}@bK zpnK;o!9|SXh99Qk`bP30yIGaQj+Er=j*#y=usiv?YhRCemh|nlkGuBW{`Y%}7$w|M z=Kc@UC}K#YXjaAC2xq+ZbH6S3I5UnmZg7t>-Y;6TsC8&b!`~VgD-ux7Isiw@*M?Cc z%XMGfyS0~LB+;(xYU_PG0)2S6J)FsppaP9Yc_@zm85to#$&$Wh45wR-GQ#LK_gqUF zJ-W5)+C3~JxXlA$gCY_m!bU|#M~#Sb*A?#5%38FjEkZNq+6NP+d>DQ~@`7*(#0z>igm(Ow+o;&6gdyX@ zM#sh^T33n<8)j{Xg$;^|4~tKT8yXdv7!wyux9M(KYZ%?T$af;6BNB#0g$;}v9yd1Q zok3Btqhb=rhb2Uf937JoHHgj|HGag1u()_5Hf&UM)bQb>?DsK>4U3FRNXXo&jtpb; zsHm`nxWtGITQW{|zYZG_kx1LH#F!CLVKK2|=(@3P3*y}e4vHCpS8`7e91t58moP9U zF@h|~n5lI^>>v0@*x0C;A<>DjG&(LmJ|=dE+y0o?jH!*uG}Su8vNth~_8N4Q`>l*! uR&LEWFq9mK`}j8h3;|g`Qnh0kmn41v*~R|FXHh