From 5f9001a59de01dc1fcc0645f401fb4533e1e8ad4 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Mon, 4 Nov 2013 16:59:03 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10935 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- doc/Section_start.html | 20 +- doc/Section_start.txt | 20 +- doc/dump.html | 8 +- doc/dump.txt | 8 +- doc/fix_phonon.html | 18 +- doc/fix_phonon.txt | 17 +- src/USER-PHONON/fix_phonon.cpp | 1854 +++++++++++------------ src/USER-PHONON/fix_phonon.h | 370 ++--- src/dump_cfg.cpp | 130 +- src/dump_cfg.h | 4 +- src/dump_custom.cpp | 2 +- tools/amber2lmp/README | 8 + tools/amber2lmp/amber2lammps.py | 16 + tools/msi2lmp/README | 7 + tools/msi2lmp/TriclinicModification.pdf | Bin 78328 -> 0 bytes 15 files changed, 1265 insertions(+), 1217 deletions(-) delete mode 100644 tools/msi2lmp/TriclinicModification.pdf diff --git a/doc/Section_start.html b/doc/Section_start.html index cc9106ebf6..230c8644ed 100644 --- a/doc/Section_start.html +++ b/doc/Section_start.html @@ -189,6 +189,7 @@ within the LAMMPS code. The options that are currently recogized are:

The read_data and dump commands will read/write gzipped files if you -compile with -DLAMMPS_GZIP. It requires that your Unix support the -"popen" command. +compile with -DLAMMPS_GZIP. It requires that your machine supports +the "popen" function in the standard runtime library and that a gzip +executable can be found by LAMMPS during a run.

-

If you use -DLAMMPS_JPEG, the dump image command will be -able to write out JPEG image files. For JPEG files, you must also link -LAMMPS with a JPEG library, as described below. If you use +

If you use -DLAMMPS_JPEG, the dump image command +will be able to write out JPEG image files. For JPEG files, you must +also link LAMMPS with a JPEG library, as described below. If you use -DLAMMPS_PNG, the dump image command will be able to write out PNG image files. For PNG files, you must also link LAMMPS with a PNG library, as described below. If neither of those two defines are -used, LAMMPS will only be able to write out text-based PPM image +used, LAMMPS will only be able to write out uncompressed PPM image files.

+

If you use -DLAMMPS_FFMPEG, the dump movie command +will be available to support on-the-fly generation of rendered movies +the need to store intermediate image files. It requires that your +machines supports the "popen" function in the standard runtime library +and that an FFmpeg executable can be found by LAMMPS during the run. +

Using -DLAMMPS_MEMALIGN= enables the use of the posix_memalign() call instead of malloc() when large chunks or memory are allocated by LAMMPS. This can help to make more efficient use of diff --git a/doc/Section_start.txt b/doc/Section_start.txt index f905bd2e2a..f785d4d103 100644 --- a/doc/Section_start.txt +++ b/doc/Section_start.txt @@ -183,6 +183,7 @@ within the LAMMPS code. The options that are currently recogized are: -DLAMMPS_GZIP -DLAMMPS_JPEG -DLAMMPS_PNG +-DLAMMPS_FFMPEG -DLAMMPS_MEMALIGN -DLAMMPS_XDR -DLAMMPS_SMALLBIG @@ -194,18 +195,25 @@ within the LAMMPS code. The options that are currently recogized are: -DPACK_MEMCPY :ul The read_data and dump commands will read/write gzipped files if you -compile with -DLAMMPS_GZIP. It requires that your Unix support the -"popen" command. +compile with -DLAMMPS_GZIP. It requires that your machine supports +the "popen" function in the standard runtime library and that a gzip +executable can be found by LAMMPS during a run. -If you use -DLAMMPS_JPEG, the "dump image"_dump.html command will be -able to write out JPEG image files. For JPEG files, you must also link -LAMMPS with a JPEG library, as described below. If you use +If you use -DLAMMPS_JPEG, the "dump image"_dump_image.html command +will be able to write out JPEG image files. For JPEG files, you must +also link LAMMPS with a JPEG library, as described below. If you use -DLAMMPS_PNG, the "dump image"_dump.html command will be able to write out PNG image files. For PNG files, you must also link LAMMPS with a PNG library, as described below. If neither of those two defines are -used, LAMMPS will only be able to write out text-based PPM image +used, LAMMPS will only be able to write out uncompressed PPM image files. +If you use -DLAMMPS_FFMPEG, the "dump movie"_dump_image.html command +will be available to support on-the-fly generation of rendered movies +the need to store intermediate image files. It requires that your +machines supports the "popen" function in the standard runtime library +and that an FFmpeg executable can be found by LAMMPS during the run. + Using -DLAMMPS_MEMALIGN= enables the use of the posix_memalign() call instead of malloc() when large chunks or memory are allocated by LAMMPS. This can help to make more efficient use of diff --git a/doc/dump.html b/doc/dump.html index eb83e900ad..3d2e4f4a24 100644 --- a/doc/dump.html +++ b/doc/dump.html @@ -97,8 +97,8 @@ dump 2 subgroup atom 50 dump.run.bin dump 4a all custom 100 dump.myforce.* id type x y vx fx dump 4b flow custom 100 dump.%.myforce id type c_myF[3] v_ke -dump 2 inner cfg 10 dump.snap.*.cfg id type xs ys zs vx vy vz -dump snap all cfg 100 dump.config.*.cfg id type xs ys zs id type c_Stress2 +dump 2 inner cfg 10 dump.snap.*.cfg mass type xs ys zs vx vy vz +dump snap all cfg 100 dump.config.*.cfg mass type xs ys zs id type c_Stress2 dump 1 all xtc 1000 file.xtc dump e_data all custom 100 dump.eff id type x y z spin eradius fx fy fz eforce @@ -237,9 +237,9 @@ extended CFG format files, as used by the package. Since the extended CFG format uses a single snapshot of the system per file, a wildcard "*" must be included in the filename, as discussed below. The list of atom attributes for style cfg must -begin with either "id type xs ys zs" or "id type xsu ysu zsu" or +begin with either "mass type xs ys zs" or "mass type xsu ysu zsu" since these quantities are needed to -write the CFG files in the appropriate format (though the "id" and +write the CFG files in the appropriate format (though the "mass" and "type" fields do not appear explicitly in the file). Any remaining attributes will be stored as "auxiliary properties" in the CFG files. Note that you will typically want to use the dump_modify diff --git a/doc/dump.txt b/doc/dump.txt index 743acf148a..f68415f422 100644 --- a/doc/dump.txt +++ b/doc/dump.txt @@ -85,8 +85,8 @@ dump myDump all atom 100 dump.atom dump 2 subgroup atom 50 dump.run.bin dump 4a all custom 100 dump.myforce.* id type x y vx fx dump 4b flow custom 100 dump.%.myforce id type c_myF\[3\] v_ke -dump 2 inner cfg 10 dump.snap.*.cfg id type xs ys zs vx vy vz -dump snap all cfg 100 dump.config.*.cfg id type xs ys zs id type c_Stress[2] +dump 2 inner cfg 10 dump.snap.*.cfg mass type xs ys zs vx vy vz +dump snap all cfg 100 dump.config.*.cfg mass type xs ys zs id type c_Stress[2] dump 1 all xtc 1000 file.xtc dump e_data all custom 100 dump.eff id type x y z spin eradius fx fy fz eforce :pre @@ -225,9 +225,9 @@ extended CFG format files, as used by the package. Since the extended CFG format uses a single snapshot of the system per file, a wildcard "*" must be included in the filename, as discussed below. The list of atom attributes for style {cfg} must -begin with either "id type xs ys zs" or "id type xsu ysu zsu" or +begin with either "mass type xs ys zs" or "mass type xsu ysu zsu" since these quantities are needed to -write the CFG files in the appropriate format (though the "id" and +write the CFG files in the appropriate format (though the "mass" and "type" fields do not appear explicitly in the file). Any remaining attributes will be stored as "auxiliary properties" in the CFG files. Note that you will typically want to use the "dump_modify diff --git a/doc/fix_phonon.html b/doc/fix_phonon.html index 0c46909990..e013e09928 100644 --- a/doc/fix_phonon.html +++ b/doc/fix_phonon.html @@ -25,8 +25,14 @@

  • Nwait = wait this many timesteps before measuring -
  • map_file = file containing the mapping info between atom ID and the lattice indices +
  • map_file = file or GAMMA +
      file is the file that contains the mapping info between atom ID and the lattice indices. 
+
    +
      GAMMA flags to treate the whole simulation box as a unit cell, so that the mapping
+  info can be generated internally. In this case, dynamical matrix at only the gamma-point
+  will/can be evaluated. 
+
  • prefix = prefix for output files
  • one or none keyword/value pairs may be appended @@ -43,7 +49,8 @@

    Examples: 

    fix 1 all phonon 20 5000 200000 map.in LJ1D sysdim 1
-fix 1 all phonon 20 5000 200000 map.in EAM3D 
+fix 1 all phonon 20 5000 200000 map.in EAM3D
+fix 1 all phonon 10 5000 500000 GAMMA EAM0D nasr 100

    Description:

    @@ -124,6 +131,13 @@ which lattice point; the lattice indices start from 0. An auxiliary code,     latgen, can be employed to generate the compatible map file for various crystals.

    +

    In case one simulates an aperiodic system, where the whole simulation box +is treated as a unit cell, one can set map_file as GAMMA, so that the mapping +info will be generated internally and a file is not needed. In this case, the +dynamical matrix at only the gamma-point will/can be evaluated. Please keep in +mind that fix-phonon is designed for cyrstals, it will be inefficient and +even degrade the performance of lammps in case the unit cell is too large. +

    The calculated dynamical matrix elements are written out in energy/distance^2/mass units. The coordinates for q points in the log file is in the units of the basis vectors of the diff --git a/doc/fix_phonon.txt b/doc/fix_phonon.txt index 9590ded071..50285d06be 100644 --- a/doc/fix_phonon.txt +++ b/doc/fix_phonon.txt @@ -17,7 +17,12 @@ phonon = style name of this fix command :l N = measure the Green's function every this many timesteps :l Noutput = output the dynamical matrix every this many measurements :l Nwait = wait this many timesteps before measuring :l -map_file = file containing the mapping info between atom ID and the lattice indices :l +map_file = {file} or {GAMMA} :l + {file} is the file that contains the mapping info between atom ID and the lattice indices. :pre + + {GAMMA} flags to treate the whole simulation box as a unit cell, so that the mapping + info can be generated internally. In this case, dynamical matrix at only the gamma-point + will/can be evaluated. :pre prefix = prefix for output files :l one or none keyword/value pairs may be appended :l keyword = {sysdim} or {nasr} :l @@ -30,7 +35,8 @@ keyword = {sysdim} or {nasr} :l [Examples:] fix 1 all phonon 20 5000 200000 map.in LJ1D sysdim 1 -fix 1 all phonon 20 5000 200000 map.in EAM3D :pre +fix 1 all phonon 20 5000 200000 map.in EAM3D +fix 1 all phonon 10 5000 500000 GAMMA EAM0D nasr 100 :pre [Description:] @@ -111,6 +117,13 @@ which lattice point; the lattice indices start from 0. An auxiliary code, "latgen"_http://code.google.com/p/latgen, can be employed to generate the compatible map file for various crystals. +In case one simulates an aperiodic system, where the whole simulation box +is treated as a unit cell, one can set {map_file} as {GAMMA}, so that the mapping +info will be generated internally and a file is not needed. In this case, the +dynamical matrix at only the gamma-point will/can be evaluated. Please keep in +mind that fix-phonon is designed for cyrstals, it will be inefficient and +even degrade the performance of lammps in case the unit cell is too large. + The calculated dynamical matrix elements are written out in "energy/distance^2/mass"_units.html units. The coordinates for {q} points in the log file is in the units of the basis vectors of the diff --git a/src/USER-PHONON/fix_phonon.cpp b/src/USER-PHONON/fix_phonon.cpp index 83f6679076..5e3113300a 100644 --- a/src/USER-PHONON/fix_phonon.cpp +++ b/src/USER-PHONON/fix_phonon.cpp @@ -1,927 +1,927 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ -/* ---------------------------------------------------------------------- - Contributing authors: - Ling-Ti Kong - - Contact: - School of Materials Science and Engineering, - Shanghai Jiao Tong University, - 800 Dongchuan Road, Minhang, - Shanghai 200240, CHINA - - konglt@sjtu.edu.cn; konglt@gmail.com -------------------------------------------------------------------------- */ - -#include "string.h" -#include "fix_phonon.h" -#include "atom.h" -#include "compute.h" -#include "domain.h" -#include "fft3d_wrap.h" -#include "force.h" -#include "group.h" -#include "lattice.h" -#include "modify.h" -#include "update.h" -#include "citeme.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; -using namespace FixConst; - -#define INVOKED_SCALAR 1 -#define INVOKED_VECTOR 2 -#define MAXLINE 512 - -static const char cite_fix_phonon[] = - "fix phonon command:\n\n" - "@Article{Kong11,\n" - " author = {L. T. Kong},\n" - " title = {Phonon dispersion measured directly from molecular dynamics simulations},\n" - " journal = {Comp.~Phys.~Comm.},\n" - " year = 2011,\n" - " volume = 182,\n" - " pages = {2201--2207}\n" - "}\n\n"; - -/* ---------------------------------------------------------------------- */ - -FixPhonon::FixPhonon(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) -{ - if (lmp->citeme) lmp->citeme->add(cite_fix_phonon); - - MPI_Comm_rank(world,&me); - MPI_Comm_size(world,&nprocs); - - if (narg < 8) error->all(FLERR,"Illegal fix phonon command: number of arguments < 8"); - - nevery = force->inumeric(FLERR, arg[3]); // Calculate this fix every n steps! - if (nevery < 1) error->all(FLERR,"Illegal fix phonon command"); - - nfreq = force->inumeric(FLERR, arg[4]); // frequency to output result - if (nfreq < 1) error->all(FLERR,"Illegal fix phonon command"); - - waitsteps = ATOBIGINT(arg[5]); // Wait this many timesteps before actually measuring - if (waitsteps < 0) error->all(FLERR,"Illegal fix phonon command: waitsteps < 0 !"); - - int n = strlen(arg[6]) + 1; // map file - mapfile = new char[n]; - strcpy(mapfile, arg[6]); - - n = strlen(arg[7]) + 1; // prefix of output - prefix = new char[n]; - strcpy(prefix, arg[7]); - logfile = new char[n+4]; - sprintf(logfile,"%s.log",prefix); - - int sdim = sysdim = domain->dimension; - int iarg = 8; - nasr = 20; - - // other command line options - while (iarg < narg){ - if (strcmp(arg[iarg],"sysdim") == 0){ - if (++iarg >= narg) error->all(FLERR,"Illegal fix phonon command: incomplete command line options."); - sdim = force->inumeric(FLERR, arg[iarg]); - if (sdim < 1) error->all(FLERR,"Illegal fix phonon command: sysdim should not be less than 1."); - - } else if (strcmp(arg[iarg],"nasr") == 0){ - if (++iarg >= narg) error->all(FLERR,"Illegal fix phonon command: incomplete command line options."); - nasr = force->inumeric(FLERR, arg[iarg]); - - } else { - error->all(FLERR,"Illegal fix phonon command: unknown option read!"); - } - - ++iarg; - } - - // get the dimension of the simulation; 1D is possible by specifying the option of "sysdim 1" - if (sdim < sysdim) sysdim = sdim; - nasr = MAX(0, nasr); - - // get the total number of atoms in group - ngroup = static_cast(group->count(igroup)); - if (ngroup < 1) error->all(FLERR,"No atom found for fix phonon!"); - - // MPI gatherv related variables - recvcnts = new int[nprocs]; - displs = new int[nprocs]; - - // mapping index - tag2surf.clear(); // clear map info - surf2tag.clear(); - - // get the mapping between lattice indices and atom IDs - readmap(); delete []mapfile; - if (nucell == 1) nasr = MIN(1,nasr); - - // get the mass matrix for dynamic matrix - getmass(); - - // create FFT and allocate memory for FFT - // here the parallization is done on the x direction only - nxlo = 0; - int *nx_loc = new int [nprocs]; - for (int i = 0; i < nprocs; ++i){ - nx_loc[i] = nx/nprocs; - if (i < nx%nprocs) ++nx_loc[i]; - } - for (int i = 0; i < me; ++i) nxlo += nx_loc[i]; - nxhi = nxlo + nx_loc[me] - 1; - mynpt = nx_loc[me]*ny*nz; - mynq = mynpt; - - fft_dim = nucell*sysdim; - fft_dim2 = fft_dim*fft_dim; - fft_nsend = mynpt*fft_dim; - - fft_cnts = new int[nprocs]; - fft_disp = new int[nprocs]; - fft_disp[0] = 0; - for (int i = 0; i < nprocs; ++i) fft_cnts[i] = nx_loc[i]*ny*nz*fft_dim; - for (int i = 1; i < nprocs; ++i) fft_disp[i] = fft_disp[i-1]+fft_cnts[i-1]; - delete []nx_loc; - - fft = new FFT3d(lmp,world,nz,ny,nx,0,nz-1,0,ny-1,nxlo,nxhi,0,nz-1,0,ny-1,nxlo,nxhi,0,0,&mysize); - memory->create(fft_data, MAX(1,mynq)*2, "fix_phonon:fft_data"); - - // allocate variables; MAX(1,... is used because NULL buffer will result in error for MPI - memory->create(RIloc,ngroup,(sysdim+1),"fix_phonon:RIloc"); - memory->create(RIall,ngroup,(sysdim+1),"fix_phonon:RIall"); - memory->create(Rsort,ngroup, sysdim, "fix_phonon:Rsort"); - - memory->create(Rnow, MAX(1,mynpt),fft_dim,"fix_phonon:Rnow"); - memory->create(Rsum, MAX(1,mynpt),fft_dim,"fix_phonon:Rsum"); - - memory->create(basis,nucell, sysdim, "fix_phonon:basis"); - - // because of hermit, only nearly half of q points are stored - memory->create(Rqnow,MAX(1,mynq),fft_dim, "fix_phonon:Rqnow"); - memory->create(Rqsum,MAX(1,mynq),fft_dim2,"fix_phonon:Rqsum"); - memory->create(Phi_q,MAX(1,mynq),fft_dim2,"fix_phonon:Phi_q"); - - // variable to collect all local Phi to root - if (me == 0) memory->create(Phi_all,ntotal,fft_dim2,"fix_phonon:Phi_all"); - else memory->create(Phi_all,1,1,"fix_phonon:Phi_all"); - - // output some information on the system to log file - if (me == 0){ - flog = fopen(logfile, "w"); - if (flog == NULL) { - char str[MAXLINE]; - sprintf(str,"Can not open output file %s",logfile); - error->one(FLERR,str); - } - for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); - fprintf(flog,"# group name of the atoms under study : %s\n", group->names[igroup]); - fprintf(flog,"# total number of atoms in the group : %d\n", ngroup); - fprintf(flog,"# dimension of the system : %d D\n", sysdim); - fprintf(flog,"# number of atoms per unit cell : %d\n", nucell); - fprintf(flog,"# dimension of the FFT mesh : %d x %d x %d\n", nx, ny, nz); - fprintf(flog,"# number of wait steps before measurement : " BIGINT_FORMAT "\n", waitsteps); - fprintf(flog,"# frequency of the measurement : %d\n", nevery); - fprintf(flog,"# output result after this many measurement: %d\n", nfreq); - fprintf(flog,"# number of processors used by this run : %d\n", nprocs); - for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); - fprintf(flog,"# mapping information between lattice index and atom id\n"); - fprintf(flog,"# nx ny nz nucell\n"); - fprintf(flog,"%d %d %d %d\n", nx, ny, nz, nucell); - fprintf(flog,"# l1 l2 l3 k atom_id\n"); - int ix, iy, iz, iu; - for (idx = 0; idx < ngroup; ++idx){ - itag = surf2tag[idx]; - iu = idx%nucell; - iz = (idx/nucell)%nz; - iy = (idx/(nucell*nz))%ny; - ix = (idx/(nucell*nz*ny))%nx; - fprintf(flog,"%d %d %d %d %d\n", ix, iy, iz, iu, itag); - } - for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); - fflush(flog); - } - surf2tag.clear(); - - // default temperature is from thermo - TempSum = new double[sysdim]; - id_temp = new char[12]; - strcpy(id_temp,"thermo_temp"); - int icompute = modify->find_compute(id_temp); - temperature = modify->compute[icompute]; - inv_nTemp = 1.0/group->count(temperature->igroup); - -return; -} // end of constructor - -/* ---------------------------------------------------------------------- */ - -void FixPhonon::post_run() -{ - // compute and output final results - if (ifreq > 0 && ifreq != nfreq) postprocess(); - if (me == 0) fclose(flog); -} - -/* ---------------------------------------------------------------------- */ - -FixPhonon::~FixPhonon() -{ - // delete locally stored array - memory->destroy(RIloc); - memory->destroy(RIall); - memory->destroy(Rsort); - memory->destroy(Rnow); - memory->destroy(Rsum); - - memory->destroy(basis); - - memory->destroy(Rqnow); - memory->destroy(Rqsum); - memory->destroy(Phi_q); - memory->destroy(Phi_all); - - delete []recvcnts; - delete []displs; - delete []prefix; - delete []logfile; - delete []fft_cnts; - delete []fft_disp; - delete []id_temp; - delete []TempSum; - delete []M_inv_sqrt; - delete []basetype; - - // destroy FFT - delete fft; - memory->sfree(fft_data); - - // clear map info - tag2surf.clear(); - surf2tag.clear(); - -} - -/* ---------------------------------------------------------------------- */ - -int FixPhonon::setmask() -{ - int mask = 0; - mask |= END_OF_STEP; - -return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixPhonon::init() -{ - // warn if more than one fix-phonon - int count = 0; - for (int i = 0; i < modify->nfix; ++i) if (strcmp(modify->fix[i]->style,"gfc") == 0) ++count; - if (count > 1 && me == 0) error->warning(FLERR,"More than one fix phonon defined"); // just warn, but allowed. - -return; -} - -/* ---------------------------------------------------------------------- */ - -void FixPhonon::setup(int flag) -{ - // initialize accumulating variables - for (int i = 0; i < sysdim; ++i) TempSum[i] = 0.; - - for (int i = 0; i < mynpt; ++i) - for (int j = 0; j < fft_dim; ++j) Rsum[i][j] = 0.; - - for (int i =0; i < mynq; ++i) - for (int j =0; j < fft_dim2; ++j) Rqsum[i][j] = std::complex (0.,0.); - - for (int i = 0; i < 6; ++i) hsum[i] = 0.; - - for (int i = 0; i < nucell; ++i) - for (int j = 0; j < sysdim; ++j) basis[i][j] = 0.; - - prev_nstep = update->ntimestep; - neval = ifreq = 0; - -return; -} - -/* ---------------------------------------------------------------------- */ - -void FixPhonon::end_of_step() -{ - if ( (update->ntimestep-prev_nstep) <= waitsteps) return; - - double **x = atom->x; - int *mask = atom->mask; - int *tag = atom->tag; - int *image = atom->image; - int nlocal = atom->nlocal; - - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - double *h = domain->h; - double xbox, ybox, zbox; - - int i,idim,jdim,ndim; - double xcur[3]; - - // to get the current temperature - if (!(temperature->invoked_flag & INVOKED_VECTOR)) temperature->compute_vector(); - for (idim = 0; idim < sysdim; ++idim) TempSum[idim] += temperature->vector[idim]; - - // evaluate R(r) on local proc - nfind = 0; - for (i = 0; i < nlocal; ++i){ - if (mask[i] & groupbit){ - itag = tag[i]; - idx = tag2surf[itag]; - - domain->unmap(x[i], image[i], xcur); - - for (idim = 0; idim < sysdim; ++idim) RIloc[nfind][idim] = xcur[idim]; - RIloc[nfind++][sysdim] = double(idx); - } - } - - // gather R(r) on local proc, then sort and redistribute to all procs for FFT - nfind *= (sysdim+1); - displs[0] = 0; - for (i = 0; i < nprocs; ++i) recvcnts[i] = 0; - MPI_Gather(&nfind,1,MPI_INT,recvcnts,1,MPI_INT,0,world); - for (i = 1; i < nprocs; ++i) displs[i] = displs[i-1] + recvcnts[i-1]; - - MPI_Gatherv(RIloc[0],nfind,MPI_DOUBLE,RIall[0],recvcnts,displs,MPI_DOUBLE,0,world); - if (me == 0){ - for (i = 0; i < ngroup; ++i){ - idx = static_cast(RIall[i][sysdim]); - for (idim = 0; idim < sysdim; ++idim) Rsort[idx][idim] = RIall[i][idim]; - } - } - MPI_Scatterv(Rsort[0],fft_cnts,fft_disp, MPI_DOUBLE, Rnow[0], fft_nsend, MPI_DOUBLE,0,world); - - // get Rsum - for (idx = 0; idx < mynpt; ++idx) - for (idim = 0; idim < fft_dim; ++idim) Rsum[idx][idim] += Rnow[idx][idim]; - - // FFT R(r) to get R(q) - for (idim = 0; idim < fft_dim; ++idim){ - int m=0; - for (idx = 0; idx < mynpt; ++idx){ - fft_data[m++] = Rnow[idx][idim]; - fft_data[m++] = 0.; - } - - fft->compute(fft_data, fft_data, -1); - - m = 0; - for (idq = 0; idq < mynq; ++idq){ - Rqnow[idq][idim] = std::complex(fft_data[m], fft_data[m+1]); - m += 2; - } - } - - // to get sum(R(q).R(q)*) - for (idq = 0; idq < mynq; ++idq){ - ndim = 0; - for (idim = 0; idim < fft_dim; ++idim) - for (jdim = 0; jdim < fft_dim; ++jdim) Rqsum[idq][ndim++] += Rqnow[idq][idim]*conj(Rqnow[idq][jdim]); - } - - // get basis info - if (fft_dim > sysdim){ - double dist2orig[3]; - for (idx = 0; idx < mynpt; ++idx){ - ndim = sysdim; - for (i = 1; i < nucell; ++i){ - for (idim = 0; idim < sysdim; ++idim) dist2orig[idim] = Rnow[idx][ndim++] - Rnow[idx][idim]; - domain->minimum_image(dist2orig); - for (idim = 0; idim < sysdim; ++idim) basis[i][idim] += dist2orig[idim]; - } - } - } - // get lattice vector info - for (int i = 0; i < 6; ++i) hsum[i] += h[i]; - - // increment counter - ++neval; - - // compute and output Phi_q after every nfreq evaluations - if (++ifreq == nfreq) postprocess(); - -return; -} // end of end_of_step() - -/* ---------------------------------------------------------------------- */ - -double FixPhonon::memory_usage() -{ - double bytes = sizeof(double)*2*mynq - + sizeof(std::map)*2*ngroup - + sizeof(double)*(ngroup*(3*sysdim+2)+mynpt*fft_dim*2) - + sizeof(std::complex)*MAX(1,mynq)*fft_dim *(1+2*fft_dim) - + sizeof(std::complex)*ntotal*fft_dim2 - + sizeof(int) * nprocs * 4; - return bytes; -} - -/* ---------------------------------------------------------------------- */ - -int FixPhonon::modify_param(int narg, char **arg) -{ - if (strcmp(arg[0],"temp") == 0) { - if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); - delete [] id_temp; - int n = strlen(arg[1]) + 1; - id_temp = new char[n]; - strcpy(id_temp,arg[1]); - - int icompute = modify->find_compute(id_temp); - if (icompute < 0) error->all(FLERR,"Could not find fix_modify temp ID"); - temperature = modify->compute[icompute]; - - if (temperature->tempflag == 0) - error->all(FLERR,"Fix_modify temp ID does not compute temperature"); - inv_nTemp = 1.0/group->count(temperature->igroup); - - return 2; - } - -return 0; -} - -/* ---------------------------------------------------------------------- - * private method, to get the mass matrix for dynamic matrix - * --------------------------------------------------------------------*/ -void FixPhonon::getmass() -{ - int nlocal = atom->nlocal; - int *mask = atom->mask; - int *tag = atom->tag; - int *type = atom->type; - double *rmass = atom->rmass; - double *mass = atom->mass; - double *mass_one, *mass_all; - double *type_one, *type_all; - - mass_one = new double[nucell]; - mass_all = new double[nucell]; - type_one = new double[nucell]; - type_all = new double[nucell]; - for (int i = 0; i < nucell; ++i) mass_one[i] = type_one[i] = 0.; - - if (rmass){ - for (int i = 0; i < nlocal; ++i){ - if (mask[i] & groupbit){ - itag = tag[i]; - idx = tag2surf[itag]; - int iu = idx%nucell; - mass_one[iu] += rmass[i]; - type_one[iu] += double(type[i]); - } - } - } else { - for (int i = 0; i < nlocal; ++i){ - if (mask[i] & groupbit){ - itag = tag[i]; - idx = tag2surf[itag]; - int iu = idx%nucell; - mass_one[iu] += mass[type[i]]; - type_one[iu] += double(type[i]); - } - } - } - - MPI_Allreduce(mass_one,mass_all,nucell,MPI_DOUBLE,MPI_SUM,world); - MPI_Allreduce(type_one,type_all,nucell,MPI_DOUBLE,MPI_SUM,world); - - M_inv_sqrt = new double[nucell]; - basetype = new int[nucell]; - - double inv_total = 1./double(ntotal); - for (int i = 0; i < nucell; ++i){ - mass_all[i] *= inv_total; - M_inv_sqrt[i] = sqrt(1./mass_all[i]); - - basetype[i] = int(type_all[i]*inv_total); - } - delete []mass_one; - delete []mass_all; - delete []type_one; - delete []type_all; - -return; -} - - -/* ---------------------------------------------------------------------- - * private method, to read the mapping info from file - * --------------------------------------------------------------------*/ - -void FixPhonon::readmap() -{ - int info = 0; - - // auto-generate mapfile for "cluster" (gamma only system) - if (strcmp(mapfile, "GAMMA") == 0){ - nx = ny = nz = ntotal = 1; - nucell = ngroup; - - int *tag_loc, *tag_all; - memory->create(tag_loc,ngroup,"fix_phonon:tag_loc"); - memory->create(tag_all,ngroup,"fix_phonon:tag_all"); - - // get atom IDs on local proc - int nfind = 0; - for (int i = 0; i < atom->nlocal; ++i){ - if (atom->mask[i] & groupbit) tag_loc[nfind++] = atom->tag[i]; - } - - // gather IDs on local proc - displs[0] = 0; - for (int i = 0; i < nprocs; ++i) recvcnts[i] = 0; - MPI_Allgather(&nfind,1,MPI_INT,recvcnts,1,MPI_INT,world); - for (int i = 1; i < nprocs; ++i) displs[i] = displs[i-1] + recvcnts[i-1]; - - MPI_Allgatherv(tag_loc,nfind,MPI_INT,tag_all,recvcnts,displs,MPI_INT,world); - for (int i = 0; i < ngroup; ++i){ - itag = tag_all[i]; - tag2surf[itag] = i; - surf2tag[i] = itag; - } - - memory->destroy(tag_loc); - memory->destroy(tag_all); - return; - } - - // read from map file for others - char line[MAXLINE]; - FILE *fp = fopen(mapfile, "r"); - if (fp == NULL){ - sprintf(line,"Cannot open input map file %s", mapfile); - error->all(FLERR,line); - } - - if (fgets(line,MAXLINE,fp) == NULL) error->all(FLERR,"Error while reading header of mapping file!"); - nx = force->inumeric(FLERR, strtok(line, " \n\t\r\f")); - ny = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); - nz = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); - nucell = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); - ntotal = nx*ny*nz; - if (ntotal*nucell != ngroup) error->all(FLERR,"FFT mesh and number of atoms in group mismatch!"); - - // second line of mapfile is comment - if (fgets(line,MAXLINE,fp) == NULL) error->all(FLERR,"Error while reading comment of mapping file!"); - - int ix, iy, iz, iu; - // the remaining lines carry the mapping info - for (int i = 0; i < ngroup; ++i){ - if (fgets(line,MAXLINE,fp) == NULL) {info = 1; break;} - ix = force->inumeric(FLERR, strtok(line, " \n\t\r\f")); - iy = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); - iz = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); - iu = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); - itag = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); - - // check if index is in correct range - if (ix < 0 || ix >= nx || iy < 0 || iy >= ny || iz < 0 || iz >= nz|| iu < 0 || iu >= nucell) {info = 2; break;} - // 1 <= itag <= natoms - if (itag < 1 || itag > static_cast(atom->natoms)) {info = 3; break;} - idx = ((ix*ny+iy)*nz+iz)*nucell + iu; - tag2surf[itag] = idx; - surf2tag[idx] = itag; - } - fclose(fp); - - if (tag2surf.size() != surf2tag.size() || tag2surf.size() != static_cast(ngroup) ) - error->all(FLERR,"The mapping is incomplete!"); - if (info) error->all(FLERR,"Error while reading mapping file!"); - - // check the correctness of mapping - int *mask = atom->mask; - int *tag = atom->tag; - int nlocal = atom->nlocal; - - for (int i = 0; i < nlocal; ++i) { - if (mask[i] & groupbit){ - itag = tag[i]; - idx = tag2surf[itag]; - if (itag != surf2tag[idx]) error->one(FLERR,"The mapping info read is incorrect!"); - } - } - -return; -} - -/* ---------------------------------------------------------------------- - * private method, to output the force constant matrix - * --------------------------------------------------------------------*/ -void FixPhonon::postprocess( ) -{ - if (neval < 1) return; - - ifreq = 0; - int idim, jdim, ndim; - double inv_neval = 1.0 /double(neval); - - // to get - for (idq = 0; idq < mynq; ++idq) - for (idim = 0; idim < fft_dim2; ++idim) Phi_q[idq][idim] = Rqsum[idq][idim]*inv_neval; - - // to get - for (idx = 0; idx < mynpt; ++idx) - for (idim = 0; idim < fft_dim; ++idim) Rnow[idx][idim] = Rsum[idx][idim] * inv_neval; - - // to get q - for (idim = 0; idim < fft_dim; ++idim){ - int m = 0; - for (idx = 0; idx < mynpt; ++idx){ - fft_data[m++] = Rnow[idx][idim]; - fft_data[m++] = 0.; - } - - fft->compute(fft_data,fft_data,-1); - - m = 0; - for (idq = 0; idq < mynq; ++idq){ - Rqnow[idq][idim] = std::complex(fft_data[m], fft_data[m+1]); - m += 2; - } - } - - // to get G(q) = - q.q - for (idq = 0; idq < mynq; ++idq){ - ndim = 0; - for (idim = 0; idim < fft_dim; ++idim) - for (jdim = 0; jdim < fft_dim; ++jdim) Phi_q[idq][ndim++] -= Rqnow[idq][idim]*conj(Rqnow[idq][jdim]); - } - - // to get Phi = KT.G^-1; normalization of FFTW data is done here - double boltz = force->boltz, kbtsqrt[sysdim], TempAve = 0.; - double TempFac = inv_neval*inv_nTemp; - double NormFac = TempFac*double(ntotal); - - for (idim = 0; idim < sysdim; ++idim){ - kbtsqrt[idim] = sqrt(TempSum[idim]*NormFac); - TempAve += TempSum[idim]*TempFac; - } - TempAve /= sysdim*boltz; - - for (idq = 0; idq < mynq; ++idq){ - GaussJordan(fft_dim, Phi_q[idq]); - ndim =0; - for (idim = 0; idim < fft_dim; ++idim) - for (jdim = 0; jdim < fft_dim; ++jdim) Phi_q[idq][ndim++] *= kbtsqrt[idim%sysdim]*kbtsqrt[jdim%sysdim]; - } - - // to collect all local Phi_q to root - displs[0]=0; - for (int i = 0; i < nprocs; ++i) recvcnts[i] = fft_cnts[i]*fft_dim*2; - for (int i = 1; i < nprocs; ++i) displs[i] = displs[i-1] + recvcnts[i-1]; - MPI_Gatherv(Phi_q[0],mynq*fft_dim2*2,MPI_DOUBLE,Phi_all[0],recvcnts,displs,MPI_DOUBLE,0,world); - - // to collect all basis info and averaged it on root - double basis_root[fft_dim]; - if (fft_dim > sysdim) MPI_Reduce(&basis[1][0], &basis_root[sysdim], fft_dim-sysdim, MPI_DOUBLE, MPI_SUM, 0, world); - - if (me == 0){ // output dynamic matrix by root - - // get basis info - for (idim = 0; idim < sysdim; ++idim) basis_root[idim] = 0.; - for (idim = sysdim; idim < fft_dim; ++idim) basis_root[idim] /= double(ntotal)*double(neval); - // get unit cell base vector info; might be incorrect if MD pbc and FixPhonon pbc mismatch. - double basevec[9]; - basevec[1] = basevec[2] = basevec[5] = 0.; - basevec[0] = hsum[0] * inv_neval / double(nx); - basevec[4] = hsum[1] * inv_neval / double(ny); - basevec[8] = hsum[2] * inv_neval / double(nz); - basevec[7] = hsum[3] * inv_neval / double(nz); - basevec[6] = hsum[4] * inv_neval / double(nz); - basevec[3] = hsum[5] * inv_neval / double(ny); - - // write binary file, in fact, it is the force constants matrix that is written - // Enforcement of ASR and the conversion of dynamical matrix is done in the postprocessing code - char fname[MAXLINE]; - sprintf(fname,"%s.bin." BIGINT_FORMAT,prefix,update->ntimestep); - FILE *fp_bin = fopen(fname,"wb"); - - fwrite(&sysdim, sizeof(int), 1, fp_bin); - fwrite(&nx, sizeof(int), 1, fp_bin); - fwrite(&ny, sizeof(int), 1, fp_bin); - fwrite(&nz, sizeof(int), 1, fp_bin); - fwrite(&nucell, sizeof(int), 1, fp_bin); - fwrite(&boltz, sizeof(double), 1, fp_bin); - - fwrite(Phi_all[0],sizeof(double),ntotal*fft_dim2*2,fp_bin); - - fwrite(&TempAve, sizeof(double),1, fp_bin); - fwrite(&basevec[0], sizeof(double),9, fp_bin); - fwrite(&basis_root[0],sizeof(double),fft_dim,fp_bin); - fwrite(basetype, sizeof(int), nucell, fp_bin); - fwrite(M_inv_sqrt, sizeof(double),nucell, fp_bin); - - fclose(fp_bin); - - // write log file, here however, it is the dynamical matrix that is written - for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); - fprintf(flog, "# Current time step : " BIGINT_FORMAT "\n", update->ntimestep); - fprintf(flog, "# Total number of measurements : %d\n", neval); - fprintf(flog, "# Average temperature of the measurement : %lg\n", TempAve); - fprintf(flog, "# Boltzmann constant under current units : %lg\n", boltz); - fprintf(flog, "# basis vector A1 = [%lg %lg %lg]\n", basevec[0], basevec[1], basevec[2]); - fprintf(flog, "# basis vector A2 = [%lg %lg %lg]\n", basevec[3], basevec[4], basevec[5]); - fprintf(flog, "# basis vector A3 = [%lg %lg %lg]\n", basevec[6], basevec[7], basevec[8]); - for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); - fprintf(flog, "# qx\t qy \t qz \t\t Phi(q)\n"); - - EnforceASR(); - - // to get D = 1/M x Phi - for (idq = 0; idq < ntotal; ++idq){ - ndim =0; - for (idim = 0; idim < fft_dim; ++idim) - for (jdim = 0; jdim < fft_dim; ++jdim) Phi_all[idq][ndim++] *= M_inv_sqrt[idim/sysdim]*M_inv_sqrt[jdim/sysdim]; - } - - idq =0; - for (int ix = 0; ix < nx; ++ix){ - double qx = double(ix)/double(nx); - for (int iy = 0; iy < ny; ++iy){ - double qy = double(iy)/double(ny); - for (int iz = 0; iz < nz; ++iz){ - double qz = double(iz)/double(nz); - fprintf(flog,"%lg %lg %lg", qx, qy, qz); - for (idim = 0; idim < fft_dim2; ++idim) fprintf(flog, " %lg %lg", real(Phi_all[idq][idim]), imag(Phi_all[idq][idim])); - fprintf(flog, "\n"); - ++idq; - } - } - } - fflush(flog); - } - -return; -} // end of postprocess - -/* ---------------------------------------------------------------------- - * private method, to get the inverse of a complex matrix by means of - * Gaussian-Jordan Elimination with full pivoting; square matrix required. - * - * Adapted from the Numerical Recipes in Fortran. - * --------------------------------------------------------------------*/ -void FixPhonon::GaussJordan(int n, std::complex *Mat) -{ - int i,icol,irow,j,k,l,ll,idr,idc; - int *indxc,*indxr,*ipiv; - double big, nmjk; - std::complex dum, pivinv; - - indxc = new int[n]; - indxr = new int[n]; - ipiv = new int[n]; - - for (i = 0; i < n; ++i) ipiv[i] = 0; - for (i = 0; i < n; ++i){ - big = 0.; - for (j = 0; j < n; ++j){ - if (ipiv[j] != 1){ - for (k = 0; k < n; ++k){ - if (ipiv[k] == 0){ - idr = j*n+k; - nmjk = norm(Mat[idr]); - if (nmjk >= big){ - big = nmjk; - irow = j; - icol = k; - } - } else if (ipiv[k] > 1) error->one(FLERR,"Singular matrix in complex GaussJordan!"); - } - } - } - ipiv[icol] += 1; - if (irow != icol){ - for (l = 0; l < n; ++l){ - idr = irow*n+l; - idc = icol*n+l; - dum = Mat[idr]; - Mat[idr] = Mat[idc]; - Mat[idc] = dum; - } - } - indxr[i] = irow; - indxc[i] = icol; - idr = icol*n+icol; - if (Mat[idr] == std::complex(0.,0.)) error->one(FLERR,"Singular matrix in complex GaussJordan!"); - - pivinv = 1./ Mat[idr]; - Mat[idr] = std::complex(1.,0.); - idr = icol*n; - for (l = 0; l < n; ++l) Mat[idr+l] *= pivinv; - for (ll = 0; ll < n; ++ll){ - if (ll != icol){ - idc = ll*n + icol; - dum = Mat[idc]; - Mat[idc] = 0.; - idc -= icol; - for (l = 0; l < n; ++l) Mat[idc+l] -= Mat[idr+l]*dum; - } - } - } - - for (l = n-1; l >= 0; --l){ - int rl = indxr[l]; - int cl = indxc[l]; - if (rl != cl){ - for (k = 0; k < n; ++k){ - idr = k*n + rl; - idc = k*n + cl; - dum = Mat[idr]; - Mat[idr] = Mat[idc]; - Mat[idc] = dum; - } - } - } - delete []indxr; - delete []indxc; - delete []ipiv; - -return; -} - -/* ---------------------------------------------------------------------- - * private method, to apply the acoustic sum rule on force constant matrix - * at gamma point. Should be executed on root only. - * --------------------------------------------------------------------*/ -void FixPhonon::EnforceASR() -{ - if (nasr < 1) return; - - for (int iit = 0; iit < nasr; ++iit){ - // simple ASR; the resultant matrix might not be symmetric - for (int a = 0; a < sysdim; ++a) - for (int b = 0; b < sysdim; ++b){ - for (int k = 0; k < nucell; ++k){ - double sum = 0.; - for (int kp = 0; kp < nucell; ++kp){ - int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; - sum += real(Phi_all[0][idx]); - } - sum /= double(nucell); - for (int kp = 0; kp < nucell; ++kp){ - int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; - real(Phi_all[0][idx]) -= sum; - } - } - } - - // symmetrize - for (int k = 0; k < nucell; ++k) - for (int kp = k; kp < nucell; ++kp){ - double csum = 0.; - for (int a = 0; a < sysdim; ++a) - for (int b = 0; b < sysdim; ++b){ - int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; - int jdx = (kp*sysdim+b)*fft_dim + k*sysdim + a; - csum = (real(Phi_all[0][idx])+real(Phi_all[0][jdx]))*0.5; - real(Phi_all[0][idx]) = real(Phi_all[0][jdx]) = csum; - } - } - } - - // symmetric ASR - for (int a = 0; a < sysdim; ++a) - for (int b = 0; b < sysdim; ++b){ - for (int k = 0; k < nucell; ++k){ - double sum = 0.; - for (int kp = 0; kp < nucell; ++kp){ - int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; - sum += real(Phi_all[0][idx]); - } - sum /= double(nucell-k); - for (int kp = k; kp < nucell; ++kp){ - int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; - int jdx = (kp*sysdim+b)*fft_dim + k*sysdim + a; - real(Phi_all[0][idx]) -= sum; - real(Phi_all[0][jdx]) = real(Phi_all[0][idx]); - } - } - } - -return; -} -/* --------------------------------------------------------------------*/ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: + Ling-Ti Kong + + Contact: + School of Materials Science and Engineering, + Shanghai Jiao Tong University, + 800 Dongchuan Road, Minhang, + Shanghai 200240, CHINA + + konglt@sjtu.edu.cn; konglt@gmail.com +------------------------------------------------------------------------- */ + +#include "string.h" +#include "fix_phonon.h" +#include "atom.h" +#include "compute.h" +#include "domain.h" +#include "fft3d_wrap.h" +#include "force.h" +#include "group.h" +#include "lattice.h" +#include "modify.h" +#include "update.h" +#include "citeme.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +#define INVOKED_SCALAR 1 +#define INVOKED_VECTOR 2 +#define MAXLINE 512 + +static const char cite_fix_phonon[] = + "fix phonon command:\n\n" + "@Article{Kong11,\n" + " author = {L. T. Kong},\n" + " title = {Phonon dispersion measured directly from molecular dynamics simulations},\n" + " journal = {Comp.~Phys.~Comm.},\n" + " year = 2011,\n" + " volume = 182,\n" + " pages = {2201--2207}\n" + "}\n\n"; + +/* ---------------------------------------------------------------------- */ + +FixPhonon::FixPhonon(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) +{ + if (lmp->citeme) lmp->citeme->add(cite_fix_phonon); + + MPI_Comm_rank(world,&me); + MPI_Comm_size(world,&nprocs); + + if (narg < 8) error->all(FLERR,"Illegal fix phonon command: number of arguments < 8"); + + nevery = force->inumeric(FLERR, arg[3]); // Calculate this fix every n steps! + if (nevery < 1) error->all(FLERR,"Illegal fix phonon command"); + + nfreq = force->inumeric(FLERR, arg[4]); // frequency to output result + if (nfreq < 1) error->all(FLERR,"Illegal fix phonon command"); + + waitsteps = ATOBIGINT(arg[5]); // Wait this many timesteps before actually measuring + if (waitsteps < 0) error->all(FLERR,"Illegal fix phonon command: waitsteps < 0 !"); + + int n = strlen(arg[6]) + 1; // map file + mapfile = new char[n]; + strcpy(mapfile, arg[6]); + + n = strlen(arg[7]) + 1; // prefix of output + prefix = new char[n]; + strcpy(prefix, arg[7]); + logfile = new char[n+4]; + sprintf(logfile,"%s.log",prefix); + + int sdim = sysdim = domain->dimension; + int iarg = 8; + nasr = 20; + + // other command line options + while (iarg < narg){ + if (strcmp(arg[iarg],"sysdim") == 0){ + if (++iarg >= narg) error->all(FLERR,"Illegal fix phonon command: incomplete command line options."); + sdim = force->inumeric(FLERR, arg[iarg]); + if (sdim < 1) error->all(FLERR,"Illegal fix phonon command: sysdim should not be less than 1."); + + } else if (strcmp(arg[iarg],"nasr") == 0){ + if (++iarg >= narg) error->all(FLERR,"Illegal fix phonon command: incomplete command line options."); + nasr = force->inumeric(FLERR, arg[iarg]); + + } else { + error->all(FLERR,"Illegal fix phonon command: unknown option read!"); + } + + ++iarg; + } + + // get the dimension of the simulation; 1D is possible by specifying the option of "sysdim 1" + if (sdim < sysdim) sysdim = sdim; + nasr = MAX(0, nasr); + + // get the total number of atoms in group + ngroup = static_cast(group->count(igroup)); + if (ngroup < 1) error->all(FLERR,"No atom found for fix phonon!"); + + // MPI gatherv related variables + recvcnts = new int[nprocs]; + displs = new int[nprocs]; + + // mapping index + tag2surf.clear(); // clear map info + surf2tag.clear(); + + // get the mapping between lattice indices and atom IDs + readmap(); delete []mapfile; + if (nucell == 1) nasr = MIN(1,nasr); + + // get the mass matrix for dynamic matrix + getmass(); + + // create FFT and allocate memory for FFT + // here the parallization is done on the x direction only + nxlo = 0; + int *nx_loc = new int [nprocs]; + for (int i = 0; i < nprocs; ++i){ + nx_loc[i] = nx/nprocs; + if (i < nx%nprocs) ++nx_loc[i]; + } + for (int i = 0; i < me; ++i) nxlo += nx_loc[i]; + nxhi = nxlo + nx_loc[me] - 1; + mynpt = nx_loc[me]*ny*nz; + mynq = mynpt; + + fft_dim = nucell*sysdim; + fft_dim2 = fft_dim*fft_dim; + fft_nsend = mynpt*fft_dim; + + fft_cnts = new int[nprocs]; + fft_disp = new int[nprocs]; + fft_disp[0] = 0; + for (int i = 0; i < nprocs; ++i) fft_cnts[i] = nx_loc[i]*ny*nz*fft_dim; + for (int i = 1; i < nprocs; ++i) fft_disp[i] = fft_disp[i-1]+fft_cnts[i-1]; + delete []nx_loc; + + fft = new FFT3d(lmp,world,nz,ny,nx,0,nz-1,0,ny-1,nxlo,nxhi,0,nz-1,0,ny-1,nxlo,nxhi,0,0,&mysize); + memory->create(fft_data, MAX(1,mynq)*2, "fix_phonon:fft_data"); + + // allocate variables; MAX(1,... is used because NULL buffer will result in error for MPI + memory->create(RIloc,ngroup,(sysdim+1),"fix_phonon:RIloc"); + memory->create(RIall,ngroup,(sysdim+1),"fix_phonon:RIall"); + memory->create(Rsort,ngroup, sysdim, "fix_phonon:Rsort"); + + memory->create(Rnow, MAX(1,mynpt),fft_dim,"fix_phonon:Rnow"); + memory->create(Rsum, MAX(1,mynpt),fft_dim,"fix_phonon:Rsum"); + + memory->create(basis,nucell, sysdim, "fix_phonon:basis"); + + // because of hermit, only nearly half of q points are stored + memory->create(Rqnow,MAX(1,mynq),fft_dim, "fix_phonon:Rqnow"); + memory->create(Rqsum,MAX(1,mynq),fft_dim2,"fix_phonon:Rqsum"); + memory->create(Phi_q,MAX(1,mynq),fft_dim2,"fix_phonon:Phi_q"); + + // variable to collect all local Phi to root + if (me == 0) memory->create(Phi_all,ntotal,fft_dim2,"fix_phonon:Phi_all"); + else memory->create(Phi_all,1,1,"fix_phonon:Phi_all"); + + // output some information on the system to log file + if (me == 0){ + flog = fopen(logfile, "w"); + if (flog == NULL) { + char str[MAXLINE]; + sprintf(str,"Can not open output file %s",logfile); + error->one(FLERR,str); + } + for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); + fprintf(flog,"# group name of the atoms under study : %s\n", group->names[igroup]); + fprintf(flog,"# total number of atoms in the group : %d\n", ngroup); + fprintf(flog,"# dimension of the system : %d D\n", sysdim); + fprintf(flog,"# number of atoms per unit cell : %d\n", nucell); + fprintf(flog,"# dimension of the FFT mesh : %d x %d x %d\n", nx, ny, nz); + fprintf(flog,"# number of wait steps before measurement : " BIGINT_FORMAT "\n", waitsteps); + fprintf(flog,"# frequency of the measurement : %d\n", nevery); + fprintf(flog,"# output result after this many measurement: %d\n", nfreq); + fprintf(flog,"# number of processors used by this run : %d\n", nprocs); + for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); + fprintf(flog,"# mapping information between lattice index and atom id\n"); + fprintf(flog,"# nx ny nz nucell\n"); + fprintf(flog,"%d %d %d %d\n", nx, ny, nz, nucell); + fprintf(flog,"# l1 l2 l3 k atom_id\n"); + int ix, iy, iz, iu; + for (idx = 0; idx < ngroup; ++idx){ + itag = surf2tag[idx]; + iu = idx%nucell; + iz = (idx/nucell)%nz; + iy = (idx/(nucell*nz))%ny; + ix = (idx/(nucell*nz*ny))%nx; + fprintf(flog,"%d %d %d %d %d\n", ix, iy, iz, iu, itag); + } + for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); + fflush(flog); + } + surf2tag.clear(); + + // default temperature is from thermo + TempSum = new double[sysdim]; + id_temp = new char[12]; + strcpy(id_temp,"thermo_temp"); + int icompute = modify->find_compute(id_temp); + temperature = modify->compute[icompute]; + inv_nTemp = 1.0/group->count(temperature->igroup); + +return; +} // end of constructor + +/* ---------------------------------------------------------------------- */ + +void FixPhonon::post_run() +{ + // compute and output final results + if (ifreq > 0 && ifreq != nfreq) postprocess(); + if (me == 0) fclose(flog); +} + +/* ---------------------------------------------------------------------- */ + +FixPhonon::~FixPhonon() +{ + // delete locally stored array + memory->destroy(RIloc); + memory->destroy(RIall); + memory->destroy(Rsort); + memory->destroy(Rnow); + memory->destroy(Rsum); + + memory->destroy(basis); + + memory->destroy(Rqnow); + memory->destroy(Rqsum); + memory->destroy(Phi_q); + memory->destroy(Phi_all); + + delete []recvcnts; + delete []displs; + delete []prefix; + delete []logfile; + delete []fft_cnts; + delete []fft_disp; + delete []id_temp; + delete []TempSum; + delete []M_inv_sqrt; + delete []basetype; + + // destroy FFT + delete fft; + memory->sfree(fft_data); + + // clear map info + tag2surf.clear(); + surf2tag.clear(); + +} + +/* ---------------------------------------------------------------------- */ + +int FixPhonon::setmask() +{ + int mask = 0; + mask |= END_OF_STEP; + +return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixPhonon::init() +{ + // warn if more than one fix-phonon + int count = 0; + for (int i = 0; i < modify->nfix; ++i) if (strcmp(modify->fix[i]->style,"gfc") == 0) ++count; + if (count > 1 && me == 0) error->warning(FLERR,"More than one fix phonon defined"); // just warn, but allowed. + +return; +} + +/* ---------------------------------------------------------------------- */ + +void FixPhonon::setup(int flag) +{ + // initialize accumulating variables + for (int i = 0; i < sysdim; ++i) TempSum[i] = 0.; + + for (int i = 0; i < mynpt; ++i) + for (int j = 0; j < fft_dim; ++j) Rsum[i][j] = 0.; + + for (int i =0; i < mynq; ++i) + for (int j =0; j < fft_dim2; ++j) Rqsum[i][j] = std::complex (0.,0.); + + for (int i = 0; i < 6; ++i) hsum[i] = 0.; + + for (int i = 0; i < nucell; ++i) + for (int j = 0; j < sysdim; ++j) basis[i][j] = 0.; + + prev_nstep = update->ntimestep; + neval = ifreq = 0; + +return; +} + +/* ---------------------------------------------------------------------- */ + +void FixPhonon::end_of_step() +{ + if ( (update->ntimestep-prev_nstep) <= waitsteps) return; + + double **x = atom->x; + int *mask = atom->mask; + int *tag = atom->tag; + int *image = atom->image; + int nlocal = atom->nlocal; + + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + double *h = domain->h; + double xbox, ybox, zbox; + + int i,idim,jdim,ndim; + double xcur[3]; + + // to get the current temperature + if (!(temperature->invoked_flag & INVOKED_VECTOR)) temperature->compute_vector(); + for (idim = 0; idim < sysdim; ++idim) TempSum[idim] += temperature->vector[idim]; + + // evaluate R(r) on local proc + nfind = 0; + for (i = 0; i < nlocal; ++i){ + if (mask[i] & groupbit){ + itag = tag[i]; + idx = tag2surf[itag]; + + domain->unmap(x[i], image[i], xcur); + + for (idim = 0; idim < sysdim; ++idim) RIloc[nfind][idim] = xcur[idim]; + RIloc[nfind++][sysdim] = double(idx); + } + } + + // gather R(r) on local proc, then sort and redistribute to all procs for FFT + nfind *= (sysdim+1); + displs[0] = 0; + for (i = 0; i < nprocs; ++i) recvcnts[i] = 0; + MPI_Gather(&nfind,1,MPI_INT,recvcnts,1,MPI_INT,0,world); + for (i = 1; i < nprocs; ++i) displs[i] = displs[i-1] + recvcnts[i-1]; + + MPI_Gatherv(RIloc[0],nfind,MPI_DOUBLE,RIall[0],recvcnts,displs,MPI_DOUBLE,0,world); + if (me == 0){ + for (i = 0; i < ngroup; ++i){ + idx = static_cast(RIall[i][sysdim]); + for (idim = 0; idim < sysdim; ++idim) Rsort[idx][idim] = RIall[i][idim]; + } + } + MPI_Scatterv(Rsort[0],fft_cnts,fft_disp, MPI_DOUBLE, Rnow[0], fft_nsend, MPI_DOUBLE,0,world); + + // get Rsum + for (idx = 0; idx < mynpt; ++idx) + for (idim = 0; idim < fft_dim; ++idim) Rsum[idx][idim] += Rnow[idx][idim]; + + // FFT R(r) to get R(q) + for (idim = 0; idim < fft_dim; ++idim){ + int m=0; + for (idx = 0; idx < mynpt; ++idx){ + fft_data[m++] = Rnow[idx][idim]; + fft_data[m++] = 0.; + } + + fft->compute(fft_data, fft_data, -1); + + m = 0; + for (idq = 0; idq < mynq; ++idq){ + Rqnow[idq][idim] = std::complex(fft_data[m], fft_data[m+1]); + m += 2; + } + } + + // to get sum(R(q).R(q)*) + for (idq = 0; idq < mynq; ++idq){ + ndim = 0; + for (idim = 0; idim < fft_dim; ++idim) + for (jdim = 0; jdim < fft_dim; ++jdim) Rqsum[idq][ndim++] += Rqnow[idq][idim]*conj(Rqnow[idq][jdim]); + } + + // get basis info + if (fft_dim > sysdim){ + double dist2orig[3]; + for (idx = 0; idx < mynpt; ++idx){ + ndim = sysdim; + for (i = 1; i < nucell; ++i){ + for (idim = 0; idim < sysdim; ++idim) dist2orig[idim] = Rnow[idx][ndim++] - Rnow[idx][idim]; + domain->minimum_image(dist2orig); + for (idim = 0; idim < sysdim; ++idim) basis[i][idim] += dist2orig[idim]; + } + } + } + // get lattice vector info + for (int i = 0; i < 6; ++i) hsum[i] += h[i]; + + // increment counter + ++neval; + + // compute and output Phi_q after every nfreq evaluations + if (++ifreq == nfreq) postprocess(); + +return; +} // end of end_of_step() + +/* ---------------------------------------------------------------------- */ + +double FixPhonon::memory_usage() +{ + double bytes = sizeof(double)*2*mynq + + sizeof(std::map)*2*ngroup + + sizeof(double)*(ngroup*(3*sysdim+2)+mynpt*fft_dim*2) + + sizeof(std::complex)*MAX(1,mynq)*fft_dim *(1+2*fft_dim) + + sizeof(std::complex)*ntotal*fft_dim2 + + sizeof(int) * nprocs * 4; + return bytes; +} + +/* ---------------------------------------------------------------------- */ + +int FixPhonon::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"temp") == 0) { + if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); + delete [] id_temp; + int n = strlen(arg[1]) + 1; + id_temp = new char[n]; + strcpy(id_temp,arg[1]); + + int icompute = modify->find_compute(id_temp); + if (icompute < 0) error->all(FLERR,"Could not find fix_modify temp ID"); + temperature = modify->compute[icompute]; + + if (temperature->tempflag == 0) + error->all(FLERR,"Fix_modify temp ID does not compute temperature"); + inv_nTemp = 1.0/group->count(temperature->igroup); + + return 2; + } + +return 0; +} + +/* ---------------------------------------------------------------------- + * private method, to get the mass matrix for dynamic matrix + * --------------------------------------------------------------------*/ +void FixPhonon::getmass() +{ + int nlocal = atom->nlocal; + int *mask = atom->mask; + int *tag = atom->tag; + int *type = atom->type; + double *rmass = atom->rmass; + double *mass = atom->mass; + double *mass_one, *mass_all; + double *type_one, *type_all; + + mass_one = new double[nucell]; + mass_all = new double[nucell]; + type_one = new double[nucell]; + type_all = new double[nucell]; + for (int i = 0; i < nucell; ++i) mass_one[i] = type_one[i] = 0.; + + if (rmass){ + for (int i = 0; i < nlocal; ++i){ + if (mask[i] & groupbit){ + itag = tag[i]; + idx = tag2surf[itag]; + int iu = idx%nucell; + mass_one[iu] += rmass[i]; + type_one[iu] += double(type[i]); + } + } + } else { + for (int i = 0; i < nlocal; ++i){ + if (mask[i] & groupbit){ + itag = tag[i]; + idx = tag2surf[itag]; + int iu = idx%nucell; + mass_one[iu] += mass[type[i]]; + type_one[iu] += double(type[i]); + } + } + } + + MPI_Allreduce(mass_one,mass_all,nucell,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(type_one,type_all,nucell,MPI_DOUBLE,MPI_SUM,world); + + M_inv_sqrt = new double[nucell]; + basetype = new int[nucell]; + + double inv_total = 1./double(ntotal); + for (int i = 0; i < nucell; ++i){ + mass_all[i] *= inv_total; + M_inv_sqrt[i] = sqrt(1./mass_all[i]); + + basetype[i] = int(type_all[i]*inv_total); + } + delete []mass_one; + delete []mass_all; + delete []type_one; + delete []type_all; + +return; +} + + +/* ---------------------------------------------------------------------- + * private method, to read the mapping info from file + * --------------------------------------------------------------------*/ + +void FixPhonon::readmap() +{ + int info = 0; + + // auto-generate mapfile for "cluster" (gamma only system) + if (strcmp(mapfile, "GAMMA") == 0){ + nx = ny = nz = ntotal = 1; + nucell = ngroup; + + int *tag_loc, *tag_all; + memory->create(tag_loc,ngroup,"fix_phonon:tag_loc"); + memory->create(tag_all,ngroup,"fix_phonon:tag_all"); + + // get atom IDs on local proc + int nfind = 0; + for (int i = 0; i < atom->nlocal; ++i){ + if (atom->mask[i] & groupbit) tag_loc[nfind++] = atom->tag[i]; + } + + // gather IDs on local proc + displs[0] = 0; + for (int i = 0; i < nprocs; ++i) recvcnts[i] = 0; + MPI_Allgather(&nfind,1,MPI_INT,recvcnts,1,MPI_INT,world); + for (int i = 1; i < nprocs; ++i) displs[i] = displs[i-1] + recvcnts[i-1]; + + MPI_Allgatherv(tag_loc,nfind,MPI_INT,tag_all,recvcnts,displs,MPI_INT,world); + for (int i = 0; i < ngroup; ++i){ + itag = tag_all[i]; + tag2surf[itag] = i; + surf2tag[i] = itag; + } + + memory->destroy(tag_loc); + memory->destroy(tag_all); + return; + } + + // read from map file for others + char line[MAXLINE]; + FILE *fp = fopen(mapfile, "r"); + if (fp == NULL){ + sprintf(line,"Cannot open input map file %s", mapfile); + error->all(FLERR,line); + } + + if (fgets(line,MAXLINE,fp) == NULL) error->all(FLERR,"Error while reading header of mapping file!"); + nx = force->inumeric(FLERR, strtok(line, " \n\t\r\f")); + ny = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); + nz = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); + nucell = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); + ntotal = nx*ny*nz; + if (ntotal*nucell != ngroup) error->all(FLERR,"FFT mesh and number of atoms in group mismatch!"); + + // second line of mapfile is comment + if (fgets(line,MAXLINE,fp) == NULL) error->all(FLERR,"Error while reading comment of mapping file!"); + + int ix, iy, iz, iu; + // the remaining lines carry the mapping info + for (int i = 0; i < ngroup; ++i){ + if (fgets(line,MAXLINE,fp) == NULL) {info = 1; break;} + ix = force->inumeric(FLERR, strtok(line, " \n\t\r\f")); + iy = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); + iz = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); + iu = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); + itag = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); + + // check if index is in correct range + if (ix < 0 || ix >= nx || iy < 0 || iy >= ny || iz < 0 || iz >= nz|| iu < 0 || iu >= nucell) {info = 2; break;} + // 1 <= itag <= natoms + if (itag < 1 || itag > static_cast(atom->natoms)) {info = 3; break;} + idx = ((ix*ny+iy)*nz+iz)*nucell + iu; + tag2surf[itag] = idx; + surf2tag[idx] = itag; + } + fclose(fp); + + if (tag2surf.size() != surf2tag.size() || tag2surf.size() != static_cast(ngroup) ) + error->all(FLERR,"The mapping is incomplete!"); + if (info) error->all(FLERR,"Error while reading mapping file!"); + + // check the correctness of mapping + int *mask = atom->mask; + int *tag = atom->tag; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; ++i) { + if (mask[i] & groupbit){ + itag = tag[i]; + idx = tag2surf[itag]; + if (itag != surf2tag[idx]) error->one(FLERR,"The mapping info read is incorrect!"); + } + } + +return; +} + +/* ---------------------------------------------------------------------- + * private method, to output the force constant matrix + * --------------------------------------------------------------------*/ +void FixPhonon::postprocess( ) +{ + if (neval < 1) return; + + ifreq = 0; + int idim, jdim, ndim; + double inv_neval = 1.0 /double(neval); + + // to get + for (idq = 0; idq < mynq; ++idq) + for (idim = 0; idim < fft_dim2; ++idim) Phi_q[idq][idim] = Rqsum[idq][idim]*inv_neval; + + // to get + for (idx = 0; idx < mynpt; ++idx) + for (idim = 0; idim < fft_dim; ++idim) Rnow[idx][idim] = Rsum[idx][idim] * inv_neval; + + // to get q + for (idim = 0; idim < fft_dim; ++idim){ + int m = 0; + for (idx = 0; idx < mynpt; ++idx){ + fft_data[m++] = Rnow[idx][idim]; + fft_data[m++] = 0.; + } + + fft->compute(fft_data,fft_data,-1); + + m = 0; + for (idq = 0; idq < mynq; ++idq){ + Rqnow[idq][idim] = std::complex(fft_data[m], fft_data[m+1]); + m += 2; + } + } + + // to get G(q) = - q.q + for (idq = 0; idq < mynq; ++idq){ + ndim = 0; + for (idim = 0; idim < fft_dim; ++idim) + for (jdim = 0; jdim < fft_dim; ++jdim) Phi_q[idq][ndim++] -= Rqnow[idq][idim]*conj(Rqnow[idq][jdim]); + } + + // to get Phi = KT.G^-1; normalization of FFTW data is done here + double boltz = force->boltz, kbtsqrt[sysdim], TempAve = 0.; + double TempFac = inv_neval*inv_nTemp; + double NormFac = TempFac*double(ntotal); + + for (idim = 0; idim < sysdim; ++idim){ + kbtsqrt[idim] = sqrt(TempSum[idim]*NormFac); + TempAve += TempSum[idim]*TempFac; + } + TempAve /= sysdim*boltz; + + for (idq = 0; idq < mynq; ++idq){ + GaussJordan(fft_dim, Phi_q[idq]); + ndim =0; + for (idim = 0; idim < fft_dim; ++idim) + for (jdim = 0; jdim < fft_dim; ++jdim) Phi_q[idq][ndim++] *= kbtsqrt[idim%sysdim]*kbtsqrt[jdim%sysdim]; + } + + // to collect all local Phi_q to root + displs[0]=0; + for (int i = 0; i < nprocs; ++i) recvcnts[i] = fft_cnts[i]*fft_dim*2; + for (int i = 1; i < nprocs; ++i) displs[i] = displs[i-1] + recvcnts[i-1]; + MPI_Gatherv(Phi_q[0],mynq*fft_dim2*2,MPI_DOUBLE,Phi_all[0],recvcnts,displs,MPI_DOUBLE,0,world); + + // to collect all basis info and averaged it on root + double basis_root[fft_dim]; + if (fft_dim > sysdim) MPI_Reduce(&basis[1][0], &basis_root[sysdim], fft_dim-sysdim, MPI_DOUBLE, MPI_SUM, 0, world); + + if (me == 0){ // output dynamic matrix by root + + // get basis info + for (idim = 0; idim < sysdim; ++idim) basis_root[idim] = 0.; + for (idim = sysdim; idim < fft_dim; ++idim) basis_root[idim] /= double(ntotal)*double(neval); + // get unit cell base vector info; might be incorrect if MD pbc and FixPhonon pbc mismatch. + double basevec[9]; + basevec[1] = basevec[2] = basevec[5] = 0.; + basevec[0] = hsum[0] * inv_neval / double(nx); + basevec[4] = hsum[1] * inv_neval / double(ny); + basevec[8] = hsum[2] * inv_neval / double(nz); + basevec[7] = hsum[3] * inv_neval / double(nz); + basevec[6] = hsum[4] * inv_neval / double(nz); + basevec[3] = hsum[5] * inv_neval / double(ny); + + // write binary file, in fact, it is the force constants matrix that is written + // Enforcement of ASR and the conversion of dynamical matrix is done in the postprocessing code + char fname[MAXLINE]; + sprintf(fname,"%s.bin." BIGINT_FORMAT,prefix,update->ntimestep); + FILE *fp_bin = fopen(fname,"wb"); + + fwrite(&sysdim, sizeof(int), 1, fp_bin); + fwrite(&nx, sizeof(int), 1, fp_bin); + fwrite(&ny, sizeof(int), 1, fp_bin); + fwrite(&nz, sizeof(int), 1, fp_bin); + fwrite(&nucell, sizeof(int), 1, fp_bin); + fwrite(&boltz, sizeof(double), 1, fp_bin); + + fwrite(Phi_all[0],sizeof(double),ntotal*fft_dim2*2,fp_bin); + + fwrite(&TempAve, sizeof(double),1, fp_bin); + fwrite(&basevec[0], sizeof(double),9, fp_bin); + fwrite(&basis_root[0],sizeof(double),fft_dim,fp_bin); + fwrite(basetype, sizeof(int), nucell, fp_bin); + fwrite(M_inv_sqrt, sizeof(double),nucell, fp_bin); + + fclose(fp_bin); + + // write log file, here however, it is the dynamical matrix that is written + for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); + fprintf(flog, "# Current time step : " BIGINT_FORMAT "\n", update->ntimestep); + fprintf(flog, "# Total number of measurements : %d\n", neval); + fprintf(flog, "# Average temperature of the measurement : %lg\n", TempAve); + fprintf(flog, "# Boltzmann constant under current units : %lg\n", boltz); + fprintf(flog, "# basis vector A1 = [%lg %lg %lg]\n", basevec[0], basevec[1], basevec[2]); + fprintf(flog, "# basis vector A2 = [%lg %lg %lg]\n", basevec[3], basevec[4], basevec[5]); + fprintf(flog, "# basis vector A3 = [%lg %lg %lg]\n", basevec[6], basevec[7], basevec[8]); + for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); + fprintf(flog, "# qx\t qy \t qz \t\t Phi(q)\n"); + + EnforceASR(); + + // to get D = 1/M x Phi + for (idq = 0; idq < ntotal; ++idq){ + ndim =0; + for (idim = 0; idim < fft_dim; ++idim) + for (jdim = 0; jdim < fft_dim; ++jdim) Phi_all[idq][ndim++] *= M_inv_sqrt[idim/sysdim]*M_inv_sqrt[jdim/sysdim]; + } + + idq =0; + for (int ix = 0; ix < nx; ++ix){ + double qx = double(ix)/double(nx); + for (int iy = 0; iy < ny; ++iy){ + double qy = double(iy)/double(ny); + for (int iz = 0; iz < nz; ++iz){ + double qz = double(iz)/double(nz); + fprintf(flog,"%lg %lg %lg", qx, qy, qz); + for (idim = 0; idim < fft_dim2; ++idim) fprintf(flog, " %lg %lg", real(Phi_all[idq][idim]), imag(Phi_all[idq][idim])); + fprintf(flog, "\n"); + ++idq; + } + } + } + fflush(flog); + } + +return; +} // end of postprocess + +/* ---------------------------------------------------------------------- + * private method, to get the inverse of a complex matrix by means of + * Gaussian-Jordan Elimination with full pivoting; square matrix required. + * + * Adapted from the Numerical Recipes in Fortran. + * --------------------------------------------------------------------*/ +void FixPhonon::GaussJordan(int n, std::complex *Mat) +{ + int i,icol,irow,j,k,l,ll,idr,idc; + int *indxc,*indxr,*ipiv; + double big, nmjk; + std::complex dum, pivinv; + + indxc = new int[n]; + indxr = new int[n]; + ipiv = new int[n]; + + for (i = 0; i < n; ++i) ipiv[i] = 0; + for (i = 0; i < n; ++i){ + big = 0.; + for (j = 0; j < n; ++j){ + if (ipiv[j] != 1){ + for (k = 0; k < n; ++k){ + if (ipiv[k] == 0){ + idr = j*n+k; + nmjk = norm(Mat[idr]); + if (nmjk >= big){ + big = nmjk; + irow = j; + icol = k; + } + } else if (ipiv[k] > 1) error->one(FLERR,"Singular matrix in complex GaussJordan!"); + } + } + } + ipiv[icol] += 1; + if (irow != icol){ + for (l = 0; l < n; ++l){ + idr = irow*n+l; + idc = icol*n+l; + dum = Mat[idr]; + Mat[idr] = Mat[idc]; + Mat[idc] = dum; + } + } + indxr[i] = irow; + indxc[i] = icol; + idr = icol*n+icol; + if (Mat[idr] == std::complex(0.,0.)) error->one(FLERR,"Singular matrix in complex GaussJordan!"); + + pivinv = 1./ Mat[idr]; + Mat[idr] = std::complex(1.,0.); + idr = icol*n; + for (l = 0; l < n; ++l) Mat[idr+l] *= pivinv; + for (ll = 0; ll < n; ++ll){ + if (ll != icol){ + idc = ll*n + icol; + dum = Mat[idc]; + Mat[idc] = 0.; + idc -= icol; + for (l = 0; l < n; ++l) Mat[idc+l] -= Mat[idr+l]*dum; + } + } + } + + for (l = n-1; l >= 0; --l){ + int rl = indxr[l]; + int cl = indxc[l]; + if (rl != cl){ + for (k = 0; k < n; ++k){ + idr = k*n + rl; + idc = k*n + cl; + dum = Mat[idr]; + Mat[idr] = Mat[idc]; + Mat[idc] = dum; + } + } + } + delete []indxr; + delete []indxc; + delete []ipiv; + +return; +} + +/* ---------------------------------------------------------------------- + * private method, to apply the acoustic sum rule on force constant matrix + * at gamma point. Should be executed on root only. + * --------------------------------------------------------------------*/ +void FixPhonon::EnforceASR() +{ + if (nasr < 1) return; + + for (int iit = 0; iit < nasr; ++iit){ + // simple ASR; the resultant matrix might not be symmetric + for (int a = 0; a < sysdim; ++a) + for (int b = 0; b < sysdim; ++b){ + for (int k = 0; k < nucell; ++k){ + double sum = 0.; + for (int kp = 0; kp < nucell; ++kp){ + int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; + sum += real(Phi_all[0][idx]); + } + sum /= double(nucell); + for (int kp = 0; kp < nucell; ++kp){ + int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; + real(Phi_all[0][idx]) -= sum; + } + } + } + + // symmetrize + for (int k = 0; k < nucell; ++k) + for (int kp = k; kp < nucell; ++kp){ + double csum = 0.; + for (int a = 0; a < sysdim; ++a) + for (int b = 0; b < sysdim; ++b){ + int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; + int jdx = (kp*sysdim+b)*fft_dim + k*sysdim + a; + csum = (real(Phi_all[0][idx])+real(Phi_all[0][jdx]))*0.5; + real(Phi_all[0][idx]) = real(Phi_all[0][jdx]) = csum; + } + } + } + + // symmetric ASR + for (int a = 0; a < sysdim; ++a) + for (int b = 0; b < sysdim; ++b){ + for (int k = 0; k < nucell; ++k){ + double sum = 0.; + for (int kp = 0; kp < nucell; ++kp){ + int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; + sum += real(Phi_all[0][idx]); + } + sum /= double(nucell-k); + for (int kp = k; kp < nucell; ++kp){ + int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; + int jdx = (kp*sysdim+b)*fft_dim + k*sysdim + a; + real(Phi_all[0][idx]) -= sum; + real(Phi_all[0][jdx]) = real(Phi_all[0][idx]); + } + } + } + +return; +} +/* --------------------------------------------------------------------*/ diff --git a/src/USER-PHONON/fix_phonon.h b/src/USER-PHONON/fix_phonon.h index 96bc04b91f..e9e0aa79ad 100644 --- a/src/USER-PHONON/fix_phonon.h +++ b/src/USER-PHONON/fix_phonon.h @@ -1,185 +1,185 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - www.cs.sandia.gov/~sjplimp/lammps.html - Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ -/* ---------------------------------------------------------------------- - Contributing authors: - Ling-Ti Kong - - Contact: - School of Materials Science and Engineering, - Shanghai Jiao Tong University, - 800 Dongchuan Road, Minhang, - Shanghai 200240, CHINA - - konglt@sjtu.edu.cn; konglt@gmail.com -------------------------------------------------------------------------- */ -#ifdef FIX_CLASS - -FixStyle(phonon,FixPhonon) - -#else - -#ifndef FIX_PHONON_H -#define FIX_PHONON_H - -#include -#include "fix.h" -#include -#include "stdio.h" -#include "stdlib.h" -#include "string.h" - -namespace LAMMPS_NS { - -class FixPhonon : public Fix { - public: - FixPhonon(class LAMMPS *, int, char **); - ~FixPhonon(); - - int setmask(); - void init(); - void setup(int); - void end_of_step(); - void post_run(); - double memory_usage(); - int modify_param(int, char **); - - private: - int me,nprocs; - bigint waitsteps; // wait these number of timesteps before recording atom positions - bigint prev_nstep; // number of steps from previous run(s); to judge if waitsteps is reached. - int nfreq, ifreq; // after this number of measurement (nfreq), the result will be output once - int nx,ny,nz,nucell,ntotal; // surface dimensions in x- and y-direction, number of atom per unit surface cell - int neval; // # of evaluations - int sysdim; // system dimension - int ngroup, nfind; // total number of atoms in group; total number of atoms on this proc - char *prefix, *logfile; // prefix of output file names - FILE *flog; - - double *M_inv_sqrt; - - class FFT3d *fft; // to do fft via the fft3d wraper - int nxlo,nxhi,mysize; // size info for local MPI_FFTW - int mynpt,mynq,fft_nsend; - int *fft_cnts, *fft_disp; - int fft_dim, fft_dim2; - double *fft_data; - - int itag, idx, idq; // index variables - std::map tag2surf, surf2tag; // Mapping info - - double **RIloc; // R(r) and index on local proc - double **RIall; // gathered R(r) and index - double **Rsort; // sorted R(r) - double **Rnow; // Current R(r) on local proc - double **Rsum; // Accumulated R(r) on local proc - - int *recvcnts, *displs; // MPI related variables - - std::complex **Rqnow; // Current R(q) on local proc - std::complex **Rqsum; // Accumulator for conj(R(q)_alpha)*R(q)_beta - std::complex **Phi_q; // Phi's on local proc - std::complex **Phi_all; // Phi for all - - void readmap(); // to read the mapping of gf atoms - char *mapfile; // file name of the map file - - void getmass(); // to get the mass of each atom in a unit cell - - int nasr; - void postprocess(); // to post process the data - void EnforceASR(); // to apply acoustic sum rule to gamma point force constant matrix - - char *id_temp; // compute id for temperature - double *TempSum; // to get the average temperature vector - double inv_nTemp; // inverse of number of atoms in temperature group - class Compute *temperature; // compute that computes the temperature - - double hsum[6], **basis; - int *basetype; - - // private methods to do matrix inversion - void GaussJordan(int, std::complex*); - -}; -} -#endif -#endif - -/* ERROR/WARNING messages: - -E: Illegal fix phonon command... - -Self-explanatory. Check the input script syntax and compare to the -documentation for the command. You can use -echo screen as a -command-line option when running LAMMPS to see the offending line. - -E: No atom found for fix phonon! - -Self-explanatory. Number of atoms in the group that was passed to -fix-phonon is less than 1. - -E: Can not open output file %s" - -Self-explanatory. - -E: Illegal fix_modify command - -Self-explanatory. - -E: Could not find fix_modify temp ID - -Self-explanatory. - -E: Fix_modify temp ID does not compute temperature - -Self-explanatory. - -E: Cannot open input map file %s - -Self-explanatory. - -E: Error while reading header of mapping file! - -Self-explanatory. The first line of the map file is expected to -contain 4 positive integer numbers. - -E: FFT mesh and number of atoms in group mismatch! - -Self-explanatory. The product of the 4 numbers should be exactly the -total number of atoms in the group that was passed to fix-phonon. - -E: Error while reading comment of mapping file! - -Self-explanatory. The second line of the map file should be a comment line. - -E: The mapping is incomplete! - -Self-explanatory. - -E: Error while reading mapping file! - -Self-explanatory. - -E: The mapping info read is incorrect! - -Self-explanatory. - -E: Singular matrix in complex GaussJordan! - -Self-explanatory. - -W: More than one fix phonon defined - -Self-explanatory. Just to warn that more than one fix-phonon is defined, but allowed. - -*/ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: + Ling-Ti Kong + + Contact: + School of Materials Science and Engineering, + Shanghai Jiao Tong University, + 800 Dongchuan Road, Minhang, + Shanghai 200240, CHINA + + konglt@sjtu.edu.cn; konglt@gmail.com +------------------------------------------------------------------------- */ +#ifdef FIX_CLASS + +FixStyle(phonon,FixPhonon) + +#else + +#ifndef FIX_PHONON_H +#define FIX_PHONON_H + +#include +#include "fix.h" +#include +#include "stdio.h" +#include "stdlib.h" +#include "string.h" + +namespace LAMMPS_NS { + +class FixPhonon : public Fix { + public: + FixPhonon(class LAMMPS *, int, char **); + ~FixPhonon(); + + int setmask(); + void init(); + void setup(int); + void end_of_step(); + void post_run(); + double memory_usage(); + int modify_param(int, char **); + + private: + int me,nprocs; + bigint waitsteps; // wait these number of timesteps before recording atom positions + bigint prev_nstep; // number of steps from previous run(s); to judge if waitsteps is reached. + int nfreq, ifreq; // after this number of measurement (nfreq), the result will be output once + int nx,ny,nz,nucell,ntotal; // surface dimensions in x- and y-direction, number of atom per unit surface cell + int neval; // # of evaluations + int sysdim; // system dimension + int ngroup, nfind; // total number of atoms in group; total number of atoms on this proc + char *prefix, *logfile; // prefix of output file names + FILE *flog; + + double *M_inv_sqrt; + + class FFT3d *fft; // to do fft via the fft3d wraper + int nxlo,nxhi,mysize; // size info for local MPI_FFTW + int mynpt,mynq,fft_nsend; + int *fft_cnts, *fft_disp; + int fft_dim, fft_dim2; + double *fft_data; + + int itag, idx, idq; // index variables + std::map tag2surf, surf2tag; // Mapping info + + double **RIloc; // R(r) and index on local proc + double **RIall; // gathered R(r) and index + double **Rsort; // sorted R(r) + double **Rnow; // Current R(r) on local proc + double **Rsum; // Accumulated R(r) on local proc + + int *recvcnts, *displs; // MPI related variables + + std::complex **Rqnow; // Current R(q) on local proc + std::complex **Rqsum; // Accumulator for conj(R(q)_alpha)*R(q)_beta + std::complex **Phi_q; // Phi's on local proc + std::complex **Phi_all; // Phi for all + + void readmap(); // to read the mapping of gf atoms + char *mapfile; // file name of the map file + + void getmass(); // to get the mass of each atom in a unit cell + + int nasr; + void postprocess(); // to post process the data + void EnforceASR(); // to apply acoustic sum rule to gamma point force constant matrix + + char *id_temp; // compute id for temperature + double *TempSum; // to get the average temperature vector + double inv_nTemp; // inverse of number of atoms in temperature group + class Compute *temperature; // compute that computes the temperature + + double hsum[6], **basis; + int *basetype; + + // private methods to do matrix inversion + void GaussJordan(int, std::complex*); + +}; +} +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal fix phonon command... + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: No atom found for fix phonon! + +Self-explanatory. Number of atoms in the group that was passed to +fix-phonon is less than 1. + +E: Can not open output file %s" + +Self-explanatory. + +E: Illegal fix_modify command + +Self-explanatory. + +E: Could not find fix_modify temp ID + +Self-explanatory. + +E: Fix_modify temp ID does not compute temperature + +Self-explanatory. + +E: Cannot open input map file %s + +Self-explanatory. + +E: Error while reading header of mapping file! + +Self-explanatory. The first line of the map file is expected to +contain 4 positive integer numbers. + +E: FFT mesh and number of atoms in group mismatch! + +Self-explanatory. The product of the 4 numbers should be exactly the +total number of atoms in the group that was passed to fix-phonon. + +E: Error while reading comment of mapping file! + +Self-explanatory. The second line of the map file should be a comment line. + +E: The mapping is incomplete! + +Self-explanatory. + +E: Error while reading mapping file! + +Self-explanatory. + +E: The mapping info read is incorrect! + +Self-explanatory. + +E: Singular matrix in complex GaussJordan! + +Self-explanatory. + +W: More than one fix phonon defined + +Self-explanatory. Just to warn that more than one fix-phonon is defined, but allowed. + +*/ diff --git a/src/dump_cfg.cpp b/src/dump_cfg.cpp index e5adf3e252..05f33997c9 100755 --- a/src/dump_cfg.cpp +++ b/src/dump_cfg.cpp @@ -13,6 +13,7 @@ /* ---------------------------------------------------------------------- Contributing author: Liang Wan (Chinese Academy of Sciences) + Memory efficiency improved by Ray Shan (Sandia) ------------------------------------------------------------------------- */ #include "math.h" @@ -34,7 +35,7 @@ using namespace LAMMPS_NS; -enum{INT,DOUBLE}; // same as in dump_custom.cpp +enum{INT,DOUBLE,STRING}; // same as in DumpCustom /* ---------------------------------------------------------------------- */ @@ -42,26 +43,23 @@ DumpCFG::DumpCFG(LAMMPS *lmp, int narg, char **arg) : DumpCustom(lmp, narg, arg) { if (narg < 10 || - strcmp(arg[5],"id") != 0 || strcmp(arg[6],"type") != 0 || + strcmp(arg[5],"mass") != 0 || strcmp(arg[6],"type") != 0 || (strcmp(arg[7],"xs") != 0 && strcmp(arg[7],"xsu") != 0) || (strcmp(arg[8],"ys") != 0 && strcmp(arg[8],"ysu") != 0) || (strcmp(arg[9],"zs") != 0 && strcmp(arg[9],"zsu") != 0)) error->all(FLERR,"Dump cfg arguments must start with " - "'id type xs ys zs' or 'id type xsu ysu zsu'"); + "'mass type xs ys zs' or 'mass type xsu ysu zsu'"); if (strcmp(arg[7],"xs") == 0) if (strcmp(arg[8],"ysu") == 0 || strcmp(arg[9],"zsu") == 0) - error->all(FLERR,"Dump cfg arguments can not mix xs|ys|zs with xsu|ysu|zsu"); + error->all(FLERR, + "Dump cfg arguments can not mix xs|ys|zs with xsu|ysu|zsu"); else unwrapflag = 0; else if (strcmp(arg[8],"ys") == 0 || strcmp(arg[9],"zs") == 0) - error->all(FLERR,"Dump cfg arguments can not mix xs|ys|zs with xsu|ysu|zsu"); + error->all(FLERR, + "Dump cfg arguments can not mix xs|ys|zs with xsu|ysu|zsu"); else unwrapflag = 1; - // arrays for data rearrangement - - rbuf = NULL; - nchosen = nlines = 0; - // setup auxiliary property name strings // convert 'X_ID[m]' (X=c,f,v) to 'ID_m' @@ -105,8 +103,6 @@ DumpCFG::DumpCFG(LAMMPS *lmp, int narg, char **arg) : DumpCFG::~DumpCFG() { - if (rbuf) memory->destroy(rbuf); - if (auxname) { for (int i = 0; i < nfield-5; i++) delete [] auxname[i]; delete [] auxname; @@ -117,7 +113,8 @@ DumpCFG::~DumpCFG() void DumpCFG::init_style() { - if (multifile == 0) error->all(FLERR,"Dump cfg requires one snapshot per file"); + if (multifile == 0) + error->all(FLERR,"Dump cfg requires one snapshot per file"); DumpCustom::init_style(); } @@ -160,78 +157,57 @@ void DumpCFG::write_header(bigint n) fprintf(fp,"entry_count = %d\n",nfield-2); for (int i = 0; i < nfield-5; i++) fprintf(fp,"auxiliary[%d] = %s\n",i,auxname[i]); - - // allocate memory needed for data rearrangement - - nchosen = static_cast (n); - if (rbuf) memory->destroy(rbuf); - memory->create(rbuf,nchosen,size_one,"dump:rbuf"); } -/* ---------------------------------------------------------------------- - write data lines to file in a block-by-block style - write head of block (mass & element name) only if has atoms of the type -------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- */ void DumpCFG::write_data(int n, double *mybuf) { - int i,j,m,itype; + int i,j,m; - double *rmass = atom->rmass; - double *mass = atom->mass; - - // transfer data from buf to rbuf - // if write by proc 0, transfer chunk by chunk - - for (i = 0, m = 0; i < n; i++) { - for (j = 0; j < size_one; j++) - rbuf[nlines][j] = mybuf[m++]; - nlines++; - } - - // write data lines in rbuf to file after transfer is done - - double unwrap_coord; - - if (nlines == nchosen) { - for (itype = 1; itype <= ntypes; itype++) { - for (i = 0; i < nchosen; i++) - if (rbuf[i][1] == itype) break; - if (i < nchosen) { - if (rmass) fprintf(fp,"%g\n",rmass[i]); - else fprintf(fp,"%g\n",mass[itype]); - fprintf(fp,"%s\n",typenames[itype]); - for (; i < nchosen; i++) { - if (rbuf[i][1] == itype) { - if (unwrapflag == 0) - for (j = 2; j < size_one; j++) { - if (vtype[j] == INT) - fprintf(fp,vformat[j],static_cast (rbuf[i][j])); - else fprintf(fp,vformat[j],rbuf[i][j]); - } - else { - - // Unwrapped scaled coordinates are shifted to - // center of expanded box, to prevent - // rewrapping by AtomEye. Dividing by - // expansion factor restores correct - // interatomic distances. - - for (j = 2; j < 5; j++) { - unwrap_coord = (rbuf[i][j] - 0.5)/UNWRAPEXPAND + 0.5; - fprintf(fp,vformat[j],unwrap_coord); - } - for (j = 5; j < size_one; j++) { - if (vtype[j] == INT) - fprintf(fp,vformat[j],static_cast (rbuf[i][j])); - else fprintf(fp,vformat[j],rbuf[i][j]); - } - } - fprintf(fp,"\n"); - } + if (unwrapflag == 0) { + m = 0; + for (i = 0; i < n; i++) { + for (j = 0; j < size_one; j++) { + if (j == 0) { + fprintf(fp,"%f \n",mybuf[m]); + } else if (j == 1) { + fprintf(fp,"%s \n",typenames[(int) mybuf[m]]); + } else if (j >= 2) { + if (vtype[j] == INT) + fprintf(fp,vformat[j],static_cast (mybuf[m])); + else if (vtype[j] == DOUBLE) + fprintf(fp,vformat[j],mybuf[m]); + else if (vtype[j] == STRING) + fprintf(fp,vformat[j],typenames[(int) mybuf[m]]); } + m++; } + fprintf(fp,"\n"); + } + } else if (unwrapflag == 1) { + m = 0; + double unwrap_coord; + for (i = 0; i < n; i++) { + for (j = 0; j < size_one; j++) { + if (j == 0) { + fprintf(fp,"%f \n",mybuf[m]); + } else if (j == 1) { + fprintf(fp,"%s \n",typenames[(int) mybuf[m]]); + } else if (j >= 2 && j <= 4) { + unwrap_coord = (mybuf[m] - 0.5)/UNWRAPEXPAND + 0.5; + fprintf(fp,vformat[j],unwrap_coord); + } else if (j >= 5 ) { + if (vtype[j] == INT) + fprintf(fp,vformat[j],static_cast (mybuf[m])); + else if (vtype[j] == DOUBLE) + fprintf(fp,vformat[j],mybuf[m]); + else if (vtype[j] == STRING) + fprintf(fp,vformat[j],typenames[(int) mybuf[m]]); + } + m++; + } + fprintf(fp,"\n"); } - nlines = 0; } } diff --git a/src/dump_cfg.h b/src/dump_cfg.h index 654582b1c6..d686b23287 100755 --- a/src/dump_cfg.h +++ b/src/dump_cfg.h @@ -31,14 +31,12 @@ class DumpCFG : public DumpCustom { private: char **auxname; // name strings of auxiliary properties - int nchosen; // # of lines to be written on a writing proc - int nlines; // # of lines transferred from buf to rbuf - double **rbuf; // buf of data lines for data lines rearrangement int unwrapflag; // 1 if unwrapped coordinates are requested void init_style(); void write_header(bigint); void write_data(int, double *); + }; } diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 4e4b4224e4..e3682bbaca 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -44,7 +44,7 @@ enum{ID,MOL,TYPE,ELEMENT,MASS, TQX,TQY,TQZ,SPIN,ERADIUS,ERVEL,ERFORCE, COMPUTE,FIX,VARIABLE}; enum{LT,LE,GT,GE,EQ,NEQ}; -enum{INT,DOUBLE,STRING}; +enum{INT,DOUBLE,STRING}; // same as in DumpCFG #define INVOKED_PERATOM 8 diff --git a/tools/amber2lmp/README b/tools/amber2lmp/README index cf43a171a6..0c712f54ba 100644 --- a/tools/amber2lmp/README +++ b/tools/amber2lmp/README @@ -8,6 +8,14 @@ dump2try99.py same as dump2trj for LAMMPS 99 format ------------------------- +Modifications in file amber2lammps.py, by Vikas Varshney +Dated Nov 4, 2013 +Email address: vv0210@gmail.com + +added support for flags used in current version of AMBER + +------------------------- + Modifications in file amber2lammps.py, by Vikas Varshney Dated July 5, 2005 Email address : vv5@uakron.edu diff --git a/tools/amber2lmp/amber2lammps.py b/tools/amber2lmp/amber2lammps.py index 08bdf38d9e..cc74cf2bb2 100644 --- a/tools/amber2lmp/amber2lammps.py +++ b/tools/amber2lmp/amber2lammps.py @@ -8,6 +8,7 @@ # # Modified by Vikas Varshney, U Akron, 5 July 2005, as described in README # Bug Fixed :Third argument in Dihedral Coeffs section is an integer - Ketan S Khare September 26, 2011 +# Modified by Vikas Varshney, Oct 8, 2013 to include additional flags (Atomic_Number, Coulombic and van der Waals 1-4 factors which are included in newer vesions of .top and .crd files in amber12. #============================================================ @@ -554,6 +555,11 @@ class Amber: for i in range(self.NATOM): self.CHRG.append(eval(Pop(Item_list,0))) + print 'Reading Atomic Number...' + self.ANUMBER = [] + for i in range(self.NATOM): + self.ANUMBER.append(eval(Pop(Item_list,0))) + print 'Reading Atomic Masses...' self.AMASS = [] for i in range(self.NATOM): @@ -619,6 +625,16 @@ class Amber: for i in range(self.NPTRA): self.PHASE.append(eval(Pop(Item_list,0))) + print 'Reading 1-4 Electrostatic Scaling Factor...' + self.SCEEFAC = [] + for i in range(self.NPTRA): + self.SCEEFAC.append(eval(Pop(Item_list,0))) + + print 'Reading 1-4 Van der Waals Scaling Factor...' + self.SCNBFAC = [] + for i in range(self.NPTRA): + self.SCNBFAC.append(eval(Pop(Item_list,0))) + print 'Reading Solty...' #I think this is currently not used in AMBER. Check it out, though self.SOLTY = [] for i in range(self.NATYP): diff --git a/tools/msi2lmp/README b/tools/msi2lmp/README index ebb52f6bdd..a70fc3471b 100644 --- a/tools/msi2lmp/README +++ b/tools/msi2lmp/README @@ -1,6 +1,13 @@ Stephanie Teich-McGoldrick (Sandai) is the current maintainer of the msi2lmp tool. She can be contacted at steichm at sandia.gov +08 Oct 2013 Axel Kohlmeyer + +Fixed a memory access violation with Class 2 force fields. +Free all allocated memory to better detection of memory errors. +Print out version number and data with all print levels > 0. +Added valgrind checks to the regression tests + 02 Aug 2013 Axel Kohlmeyer Added rudimentary support for OPLS-AA based on diff --git a/tools/msi2lmp/TriclinicModification.pdf b/tools/msi2lmp/TriclinicModification.pdf deleted file mode 100644 index a88bb902f0da05d2ea9e776986c6732e0e37b480..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 78328 zcmbrl1#lfb(l&a`%oOu6Gcz+Y#>~u&am>ujcFfGov17)V*)cn2X1dP1`|j@d?^pHT zx_3_1NF`~e)jDd;^Qc>7ilX9lO!O?UWc_Ld>BegMfF%>m6wl^^a^7F$wyEvH|+QNE(Ok`@ut~Fx}{ko-k zD1~6xQd(5nuqhE&<}3}=ZA5ZZT0-M-J-IbCOPqT?^d`N0=|_>knOr~$2B{Gzp4#)R zCM}xUr?_5TBDaoOH*sBb4K7=iYi6?59eM9q+$Y=*Z1bvbG3(%Xys+;{=JHfm-8p~D zA8Y75X&-HDldVoTt@8Y0eB{{q?Cfk^owGvAv9Q+3B#?RYc-HrY#6|Y$GQY%p=GWS< z#U}&r_$ZRI?k>i6Rt{e_-!^PfeiQpKYU?VDGN`+K&gb1D%dI+&nB`co2Cq0Hle@@z z-wG`IwIzm~CDOeSDEzAmJ<4@+uN@S;R*$b$1!wfA!2sI9$a`-65c&!|R$cY^690R) z9(>BkLGL-rBy1gfsNbePOW5v;s+hH#<#?~kN}aD)9eK9JWJ&;>eJbo=t2(3?uxF8j>OO7KW>(}mn1*5C@|{50 zN<5QGR<>Nav^4Rpnn}>@W-5!#WfbX*L+1wMr3o_d3iGW4-zB|=6z#1EjQ=-~9GM2A z@Vh)?yGx3yKT(HXODvS&jJqaDJEI`=;-c) zna3qGCp3AGai=`2`f>{#3XDx3B4@0o)U*5nFO51V=ExjLcj;5^mVCA>jA!md=k>Dz z>dzfdy$mO#Q+jg~(AIQbu3DUt!-PUf#4~PiB!74;6Ds_U69mATaktH2eSDF4woD*v znqwR9lWUkA{AYh$&0Q06a4XtoZ;moo#k2%;3ZzrKZw#U-w`7_N>*N(C%jhMSq{iQ4 zBzmmc_YqV22LZ(yS3z}?gCZo{E||QRhxFDK9n2LOPR5fz8DGsjbqzYyF}Y^qa)cDT zFoKH83CG0g3qM6wZ=wS40@@o%m~4L0+^q$mx;Hb>qu>6F)M@^`%0ZEsOP4ZDUzD68 zucdxAnVYgGC#tA=t!-Z@ac(G3pF@y}X+C4^yo?w2jl9`gG-nGGeg3^|Yk86gYXy-+ z`IsKdskDIIv&Sul9Kl&uQXN6o$+Kc#we3c z;600sdgt}(Goj7r?JDS4G9}q2!CJ0jWlcr|#sVre3RsAPR}fe;t?LI$yffZu$Y_y- zuLlRBJl#Y2r2vezBdlBB;DEW%lqVCoE@Y)e4zCzYgUT~T5ZI+re=b5bYH)E{Qi<;G z$4c&n-ZmnN{vpuO4DN0})kfmWTfCzZ)%P0VT#ULP3V#m;=YxU!<#&H0QKf+oX66=F z>)mpvbsZv-YCjtrvs$ely`s05tFn?{eHpGJz}O)W(HWq7dfiud`xQsBboN!zhqk>^ zZ~2yVr(expdu%~Xv?;nEMSHnD;v4z$MHx*cIE)lL7O6snjKR6d*GdDp<-Br;1h4q# zR#z>r=chiks|*vh2|f4iij2dM9@H?X95(`yl7>aenX$3x`Mi%gb?`kDN*oX@hA5?ep-h`l;5&2wX5v!Alk*U@ zM~je&;FTrc*cfPR$rk})Nn=&}yiXS6JH;RUd`|b0r>sA(-gnMriOCL7!@XQB$m3|I z@uy932YUFxlEHpy_VCvN)vQAiX{`DIrHDL5G$9fKjX>>{z)UzCr)Z0eI!EDc)-w|z z!zCy%A5rH8_GpHes9FI}E)159+a5T|>iLi=M7~Tn7*6iRl0xnipHFUZ&IhA=U8M6Z zi0B4}Qy6j#um8p=CQqtr#cJ&l78W@t<+mpdEYWh~R5;*w;9LYnd;td|;Z0gW2EFFt ziCPrum!Qa>h({!DK~4zz=U*n->MIV2g7M~-n=L1xt@^8FRre}0cq@!z;z+S2y_^!Z zNhqUCb1Z@WJR21nqSKFGclU0qxwEsO3;om>gBTM{&+30*iqY{1I|&QqJ66ViHw?uB z-$oRwC>PM?6QaY`=LqZDP2Y~~W3fr(#O1P)I-i~7UPV%s)SzmHzpT}@=KSs^xKr_H zt6AASx`kCIX`>q=YDl0T2nBg))%hSwpldZorv<(kGwHQ4N4+Poe=O0hfZrLs*3=9dZ(YD^rv-(Y2Y zP(7x=z92X_E~jB!ZG2dgOneFbf{TSuKJ-!(6|o6Ksd`3rGcQfSop7$wN-#jH>%AY( zzMj$MvUqyLCCIt!@0s22?vB#-Wf+8uC770@u&`Zl-%PYdf|L}%{9EB(eS4D6O){R| z=e6ids5P^(l^_@+KJQdbcVY%7!TD|N0 zBXZ#aegPz zN%F}D)Oueqg}R0hl*5AcQ^hVOnE^drU{(iDi+K?|my2q<%_gwesgL*Qs$0bQGIe1RzL0`x|nrUtEdGe(+e@lMO0jET20E!wI2T zUt`Qac=Sw&r}=kp;P>?{4u3goXdh(gJv9<}7tkEu07m!ype}NCIt;^cs>;0kA>`ot z^)TaVf^`}e@XWU0#dt<=z2(My4Y$~On)*rhZ0vQC{A9Fy5BQDQGm5q!Fx$4~6`1}I z03nG@J4Oe~0iK96Y`WvXjLeA9IEE485H)RYYF zX@=;Wo()@l-}dL%J?{o@n#8rzH<^6LT|seD6s4qu&K#H9A60C>e3I+~za|X4jUZ)2 z$gMBzD>N_(+!1nM?2j}(w;n!?#EF;^0_{0$lgB<&@|`8STjwU#egx&aQ zthzPP$*L@7s>jepR;LWLQ5~%;TzyXJJQI6iNDVM+g*CM^`6ufA$o&D!e-QGYl$nEz z?N9Ar`TqdYs-6y}Kn6ua^S|y+rgknsmVbb1Wm9K+S0`grXCV7O3PkMfTt4cZfq#Jb z2kMqHHL)}lw)X&PGkz4XGcf}>xR`ZeKY;f?<@|Z>Z@{kTWN)lu>H^gM5GE=PWKcEr zZ~-z%*?tHX{%0!k&r}Mi3uF+nx3PCpaWFJC1^&VNBF;=ejz3~yKeEg~rhoL2pZ^2% zfAs3FCgRNWj6kk`bWNP;BhUHIJj+L(=}#AxKS%+PK|{gF%GCJJGq4}>nSe}8e^>px z=YN&`-SvNRjSv0)Qu`kx_V-&DR9uZ*{u;x-bOADawlw)5Cd@2MKn6)uOLGesAQKnY zhh`sR=j34jhoHc^IJtg|_J3%|^v@9$wsdh;G<6cOw{@_$`}0?D{+ru~im3h_7`zk;L009OL4h{zKaf5(>fP{jFhWfaW;NW23QIJqkQIJqj z&@l0^(9m%(P*6S*eZs-VCnO|9!y+anCLqNlAS8JI4uFFO3;_;6Ku`eh{Qx`w2mtIu z#D|c-FK`G*Fep$E=#N@HWB>>l7#PSOSr8C^8i0U;fde2=AW@0XpwNk#gfK9fl?+Kf z`zOt@2%9-3?vN|18O5s@=S|PPuK=L`yj1}3HdaI>t(^2|k9@Rm>UIUcC9 z_Pj|{2FNggtQ^(&q46(Ax>VBY%C}{tT0EjHD<#s+L-o@2@;S#fq|MDbF*l;1Iy9@b z;S>RJLXrnrYCBjRELR0$M(6f!qRn2H8KOy4_H1u10H?wWPdrDJ!~c}u8K}NYfBk8( zZg-V@EVY!PzlFWIAepE_egV0OO#kF)RPU#B#At1nQI{F$nj0J<58j}oEt!l-oR>+A zN__87vRmbO-8&%I<^0ySLA}{yeStai?h7TCCLd?{C=yuWmV!X&xBcP;q|cnjZUts- zUDj=~{#f$s8ga&}>po<5T_BS7eF-qin(nc8*HKaPZc>P5S?$j_wmF+=P6o}d>9>sA z!R=I}#Vbc)L}Zl@QoMnUJ4Wta(f{ix{!25n4%zOB*$la{w=cGJD<- zKTc&yY(41)KM)E(Yo?Xg16k_1kd2;xE44|QOCaw~p5qP;kot3I+>RB1!V-3U2(OGk z7TPwX&z-lEJ-+Qc>T253Lnp;xs z|3Dc^Kp7Ktrktsgp7ag?xZ(|iMXzd)wp|p=P-8^w4(tzso?{>qg->TD%;*4U$t+<3 zev*?&oaEhlU7AZ7H(1p_ag%;Pij30>gGD#;t|?0mOTGh`w%dB=*kNW+2r#aesSdD@ z^O{ZbtW-P)FF;1*C32Guo;AAHDU)^h8aRhBzVmfsSdnN5Q<1CD*(pH z6f=1UwM+5`sUb5Qe3v=A^S6(STK}BLmF^aILuXjF!0Y%_;zj$;^ajOMSXYC5mRAg0 zNm9Q53!KSrh-MC^?-RgQf{AG=CvkyYM+)Lut|E+?sv&i!?TaAY2 z#CJ=V=T-k?zAA^qHLCbE7Zy=nL7uuk@cZtA%GB0xuCk`6ZZVylC0ax-P>~CzlCwG! zPGc%rh|U`fB3O0yhlC%UXwR+u#!TMijq)u&?p&t>REpUnD66VvEeZsUx0!axFtYJEdIl#bTYLqy&g- zDvW^W_;OF)HG%4`T#YT@>dd+~3VYrcV&aRS<h_~%fI7l{rP z9nzT=S2-DtN~%>oB#_+T+p#7U0vc)F5{FGq;J0Ah?S!~~X%mSzuNrG_HRHx&7dFlS z0N{lNhyuxy|6`>H_GuFA)-utmwt?D~ltcuUMr{Ug;RQfLT~uG>4+Y~%L@n9Y7SUjn zCQ~L;iRMU=CP?ij$hj3!so5Sri@4v{l!ov;Kf}dfpuz~F2UXZPpraO1a8AMy6j%Mr zo@zqZEko3=Cy9zPrIZ(%07-@L2~t;}b@=&fai@9m%eN#_snj%;S}SH6wvhS|3{>A8 z>$+hB{H35iDJe=f*G42G1XFE3MOu>m%+50>#N=Z+h&Yy3u5XE!AfI2M>z-wnYl|eU zAU!;5Ka$2X9X=NomkIWELm%OS00zWQahAESc4sFWx3k~AhNKx$xoO9n6t<7)JhYBr z_kzWEf|hv`r`osGd3UT7v$m?3hhm2fo#X*#&}Mgj%#mVYq7CZY?ny7`*~x0 zvG|O}(E4Ms!rHF+w`o7UGH60W6fnR;X)OzhUt+R;L5ewR@pLe{)ceF^e8ldc?Ig4i zqezZ@Rg;L1)D<)-3qCF}F)n&uxK1anb7@m|JFlOKo5e{>#t$N>ibxO{2D|qFd~vR# zsG0WRC+V6I6SKN^r#pePjyHifs=QikkvLAGeAyl|!$SU%Y9ke>v3Ee6pau!QB`GK$ zGmEC2vD5KHPz8CJdFQwj%_X;!plR5Mi|e{NM$^wmJ85MhGJQlw)$G^5{!tTeu@;%u ziCHQYc>&W{H&nZ*=9Mg!aN82}L1kuDzXwwZKvogd7A&#lEVtfSNwLDLqAAtE+}YzL z>h_$pa>6$}@NdXLQH~ceG$PWI&x8-3<=+9HpV_R zB`iyW{8PMMCb{&^S{_R4hSZ1?@sx?Zu9UdGDg7ny0Jgmd+uzyWYSe{h7fm%x*JcS^ zz(6*)QgBu>^Kn+Q^EfSa;Dds%3JbA0#0grhlcn4S95EwV8{6zc*Rod~+PK9{TdP_T z*g9aSuXyIa3$y&jJQXl(o^so=tJgLh_n_n<3M;i80s+Pq0boV|Sx*u1%rxiinDe7L zszE!(QWk-?m%ywy{^abFa+=dSpJII7QSIS7Uj2G4Ef{Q50a6IZm*ZG#WK(bAU!|vx zYUGC{B`y3sa&ef?=5*boQb5#Z9S6Xd#fCZ9f`w*JI8g)^9K6tQ{h!ynDSj}9T zjte-H@*=UEP4Iru;(G+bEamkYE;ql`J1U-}S5al)Th>u;(4q31W0y%yxr8KFw{|Rd zW*YxbqgGmNv;|byr|(gWma(z7 zMw`sV*-ZAP>?B~;HPO<^^MZ*P3`r|{EsdAok;6d@l8u+iV7@ODhu(V2oS$Se-Y_gE@r;|d3N`Rz#LOhVI_^3 zG8=WK1m)km`%=+}8w7Kx{B1Itk`y$5tW1w(FN;tI=TUBEsUo7PthJvd@lqU#qMUPm z)H#58k&rHD%Psd%3Kl7N?;p6;XPQduhnt{7BtLMp&B5D*by76DR`N$%`non`n6 zwTXA|)quNfGjK2x4HS%8h$>2lhEsuzE`c!zDwqmj>$2W9rWS(H^=Ky-HwF|S>kZ~V z&sZb*Stcg^ebYRmy!*EgpI2v#OEyTVF4sAcEt$`*%~DcnDJm#G=_v{92wq|)z8vHs z|L6%W-;$O|E4zFCkfx+U^>iEO=Mf3cM|E+1vnz=9d!t7Efg1VsY~gu3E8KRCav5iBQh!KiLwsg zSv-g9P{M-|gUIJ#J*Ydw;T>z>woZYw%FSQdtU1nh?l`EWz0Md8RndocY-qrcc-PT- ziNeri%{!sSJiJ#vO{=+KbLF8-@A+KtrSH;S{DhfTu4|Yq0TNegbb+G*i ze2s$oEav!=_F3Y$dkH-U*A~_t3m7SB7~z&8m?<@pdq@_^#&x1Q%ObkkZ)?F<7^CQ8 z{RJ>${li*AU%q4+&ZF+MEw>^UJh!nOJ7idf*)aYvZO*i%6h?GP^7rF<`mFW;>PxC( z<;^yDzW$n9!Jy?^)ZyEZ5!T)87`?s~>KKJ*yLyRt5T*Q_S^9Rzmb+&wkXO>y9(!)A z`KCVmgya3c7r*Gm&PfWVcxx;Kk(z9i73Mx5J6R#-*6Y+Pq}VlAFWy#MY)~UqDbm>{ zHi=R6wStHw#SfRg?H>trMVmirEBv-RvnzOBTJTM7d<)6>|5|Kin<%?cbuVrva|)-T zyKL|F+nHoFVK!X~%+ZV*s}~qW=>+M-jF|E7Pft=g!L)shWfjENp89~D3%YauF+-k~ z)+E%Rouf*m!3Ms}{H@}1x3j5aD(4-ooAJ*kom!Fz0-eoS$jn}u_@Bc{Uvr$Cq@wN` zyGk14twl&$iKJsno*3Nq^o4~RW~q<-+Zu+=`R(kUSr2a~&Ro0YBE~GC|)_0;qld zg{F(Bl;}BROxql%9a{TC%IG1}WX{RzkERk-61N3}lwt~J%X|hAp5_QFw|_*ZN1uuK zp(=_yVFV?{@*|BF3ueM@_Qr$ic-jP)K@94^$?s;N#19H>M#M{RIfP(Qk!}gy;Uz9;e1637W=IxGr1jqfHvF~BL)AczOhPDTV_S2jamwV!YU~-a9D{z zbWtNdJz`NRDUKyUHNS}RWi*^Lx2(GHpBzPLU)G3NS&tsIU`xifY=AHaHBO$twI)++ z5MyO!VsPW(E6d0zR%79XgAh0Pm8Q&3VGOV=_-Z$!c^UW(!R%d`f$FBhk)Pwi(mucBBUpT|H@! z>P9PoAE4^~nkyP?KMSv2>^ebRy8afcN{L&>uM-S<22u<`I&lh=RgbgU;8jR zMV%y~UczxZ{>5|kA%m)w@S8{~deV0-^KNGYr|ZS($OkaZL2HCTcgS9icdRq9I8 z8f_$mnv*5*r({?rY>!A^Wd(EQE4H14^ zq+7}~DH5FIEUkG?W4Eh#UK!5Ss&W*(qDuJkVBIh!l&Ky3zi&PyPwgCWmNm(hopazy zxZizE?6N|vK7E^XBj)BD&RHs!h%;sN>zTdzfj7MOaxj?q>guD)g4ncd{}P{8D7VV4mH!Si(+~tcX#an)_Qf zi!NX5Ij8lA@bFdvcrk;rOfUTYuw9pXC~X_A4lW|wq=WH=V|jX=GWU<|e(Pj=FYzVX z6C0vJA|h%cDQtBeI4;u5C?C6jD?puBC1c{Yy$e~K+&-6c)D)twC9?!otUEaI;@{%t zPoJf$?@nW1&DI|aeb0)1q+%M2Ejek%Eve5G zXTbqp>iIk|Rw`=PWo);W)(Vn1aF%O&S&4Ov<`o;S zt{9zNpREVYcm7t1Apz+JmKIFmCiGV>mbShVyvt{ARVorDnmFhtpD=8aoHMNN7G zfGC&!Yr4`*oOm%i<8U7q*eM3ZosNDfO*O?_td;gz3|yk{(v#4X8pbOFo&hd3eWu$} z_LditXHHj%IcOA!TJW?Zh4V*MJtmn$;^L}S@s$&x6(tdp1DyTrTp$x#57#C8^6B>n z!ZDj_TR&sJ-Y~x1je?tF7ICHoC~qpO3JRF|Z*AYJ4@5+&l@DPP3wUMBX&ee}(FK($ zQ1GY6geD-3$7;>YRJ3zU-Qv3Metrim6#bAbY!n?z_ncUxDvmTJq{X1g<#*OC&r~_A zaLYe?`2oc7xOYJKRE8jr0?^Q0YO%k~^Dj&LA2g~Rl>|L8x0u9*Jm50d}z6fp?> ziM+D=@M-`VL`~f+jZKv$KKwcVw}7nd6n9mXWz5j2E{?9BsKh|B8yqhf;0%u>2DnnV zl%FDKND#T3D1UR){P-0^-bhy9(qeYl2WW~W=B~h`x^(bqY^D4IMjf~eJ z&0ot0EsM+M=RS+=0L7O8@IDs}5I@#D4fM0)UL91lvHp9|PskvM8X(`d%*?8}82N{Y*MmX?qmOofqDH|wM$%Hrck^&GC$T$ObD}yFH!{e!3ZlUCbkfK+%?rLj znRcm4WdQc>DQ9ONvSa~_=;`|72jqD+;ml&_*k|EapaAHkey&a@#K=c3iyVMQvH5NU zc({~NG0I|^<=*dMv!}}jl}95cnxr3cTVWouQ#Xz9F^pJ#HTEfDAY_`h=@#*)Rd}bT znhX0Dk3Gx32)?Sk@(p`Et6h?xOyvkg#q?< zZvq?SJHO@oiJbxp`)jr}x`nGgz$chXu~iCc*q#p{LIOr=C*TF@bF7Wb1UuZ_X;Bql z;~dNE^>*+?2&WTawD6Ihko*#~cH|m z5keOOz#7XaE|{k}wK!C`2yX>ko`+)Gj6I1N$|6e8<*;zH0ev5YvxLfp1!$#}V^wO(v0(I#%b!ZSo`c|$SnSk0Yyghb;TMKl>%-s>?oPh zbp7fEu?E_;sdZpWh9^9Dh(q78xA7dlwtmuKwqvUXZ$Vm&>hL@Mfv`n>>Cs8Ofqou< z)a!g<^o6J!))&v0;8!f+U=Kwcj25_c01{COkkp1$6&VNm0QxlurH@s~zf{46I0V_V zU&n~5E@nxnQwon650QoNEg zDM3IniIN;$K^jj?cv+rHx=qki&{L#H8KKl}PNO`nJzGnmOVmf&he0BdCG|RWXQV2z zBHlXrIMtg>A!#wyJyoARmqMS&SNg5^r4>}ESFl&bzo1qweqR&eSSp<&y2RQvQ}r9Z=3*Vtz$#A?U5VX}CAzp1K>?w4E7sSi(v zZtXf|hkVD>9pWAFGx;+U6kn)EXcyxBl8YN{;l3&@H zlT?!klLc__SQl7xbf2829k*KXoZ{(C=%e&K9P;MWR+D$Mj-STm7NZZ``fdBym-v=n zTc2r=WRT(|B8QSn5)W{u;Pk>$Yd*(+u4zfQH#)aEFX89r59g2ON9%U$uJB=aFnk4l zb$eRBUc5MXT6|Uj=LD~W`3+kLE)K2^#Sg~}g9_ma{-YP9w6YAFn<9Z&_KrcKv+{VV_7vZT1A0xFiW>9PsSpvH&BSz$86h*}z57CwqDh99g2 zj3Vu$mh)$7X!(NE@v@Z8nr0noF>kWBl2N3YWPzG+4OaCV<`Sm4qe|l_`!~B8yLMyv zBiAF(U)jvq++EGP-QGf8dfCpb5vzh1s&!V|rW*qpALHB?!xrV{H--x`f=M{FyQttA}) z()9WKxITUtblO#PSM=<@=?3k-wtm=tt&wh|ski*IdQ-P*Z~Q9%s-@%eNw>P8wxa=( zWvoeTBXenemvoJl#C*w7;gLU1Y{q`qk$1~!_R7%eV$UyI-&v4oumLE30(87i-{E^j z>nb7J0mHSe${$n5`A6+N=DzyAqNkFV@lz7P5=BFOv99wxbBi3M97*hb(FucfJB~y4 zuba0b-e?tviYx_f58#|cT$kHJ4}|AA@Vv{O5kBRI6tB5=%|EpwwDK0#7Q&kg*)6?g zy?-70&1U0I;?oIazJIxix=YV!7&Vufil4;D%H+y*|9T}n7onvwdV!!LCq@6{tL0(Sxoy{>L&Iyz zPE|}dgwMf~?lZ&2UAxa-%zHL2SR~}r8GX0+yWSe#q0h$aG*lUqp1`;l=c~o5@ySj; z1*`&fPPstU>-@u_Sy|W3$+ty=sfIA~*pOBchWF&h8igMaQ)N>$;*H`*5>*mYQIApU zUUfGE7h_SWvpc~(RBz4q$kAO1rt97Vc~hOn&dqx?7$X zmqct-bPL@1zvev~-weeco+vyhBlHQWDnA9Vfvzu?us!;6?dLXZB& ztAApI{$SXTFe51u5n)4TQxo7{7*-jm`*$0rzZ?B0ru}EI(Ekri8>~8IjV6Tt#%@>k zlNegd8nNCZj#x^w!a{70{T0Q?8h(|qSuD@*MR-H*WQ z!h}IOSWAp@qk_^X4~y<{!cMli=;-_8*Xs?Ab02%RUjPKB8P;1r@Lhx<<$M%Nn=isv zFJXX~d>A2CGPTtcc~VHBK<7vAgC#E{e{J<2-p}s!?o}Bx%^6TZuY#Kr8OzkQ=U9%g zny8>*yub3OhFqY3?JQlL#B#}ddh|9S(eU+a5`Ph7qnq@VYgL&bn-vU{X?dHJ%qXnnzfvpc+A6bqfU)6D$9^o>r`c?f(I=~^0bJWDp&sj@SVNurZ zp*$u13PNhM>0XKT1*_ja&3Hh__pb%=$jUttHL8W)H74fU#%pc&&oR!_thR1^lmR- z%qRMdx2c%O-69T%!czz6FB6tzD&V#aY=@w;lS8OxY~!Cu=M4-pI&T?t1++c7#C_p! zm6EW>XH>)Qefy|QShL#@?Vn;-hhDID?;{R!DM1-*Hiy&OZB-o4k@HN8?3}AXv3jO` zE}HGE3%+u`fNStpU{)!!%20MufaCBpAXt^+^%v2XBZ7yat9z%YuRnLmMh0Qxt0WN8}R+ph+tA z2h77t;v*~JQAWUWVilHBHfapNg!<~&3O~PVz1Q68FO03;XF1<<+5bwc1oOii#0&?g z1OV02sO}Jq2Wg|ii+cAKBPT|XnM8y^9RnHVhM$ifVD;|7!1dQ^``WrQ-Z=zi!=bN( z@=gJgN~z`%%v%k()oPr?OfgVI?3-82oB&0}Fy~h6EY+|@P`^Fr`u^pgkPK@OVS+fr zU>M&(T>L;7A=JNPastdZ5j*__pdnF!khCCpap21SpJ<^x4e%+!2Mk!25UG3h6yc79 zNH}4&zDqe^>Ht1}XZs0t79hk3n&aIQ;wb^5g?|BYZxXg#A=I%Cn^}u9rjN0Q&z%I#bpye_y(U8 z4*3n`w-GCfoP_fX|o$zFNjkSsA==x;2_8*&{p7l0T}*57o^~5Mldko zOonWEVT$slM69SQ(7HjmK|+0>w(!OKl9+JJg1L;pYLE{0&eY9oz^I4g_h;0}pFpw) zMH}aa)r->QLrjMc6VsEDqo$(4h41xqZ#!;#8WB1MJ!FUym&MCUB9aOy!+(&E1XbY` zIV<9(L@%iyg)q$Nl|s$Hj<~=P8z=I%xX#2cXg038V17KJHxpyEnWmZmF;c1mUzPES%#Vvl!=kBO5{(WB5MTBr2U zd&;?y!N$h+!Inr@#O@x~N{3A+Va75IG1*W`F5OpEu)<@D{J{HuEAm}(K2$U`*0jg; zxHOqm9a7DDcUdQ3%Pmf$P76-7T-9qhZPaaMZO2YidsJ8ASF(>}S8Z2&S3|osqp?(q z=usHF7)!7f80F~ll<<_Yl+-jUl#n#V)G3q%%AJJ|qT!V$mAh5_YSfA=H7@d&8tLC` zzMYhfE9uqfRD99%6b&vk%I{N_Qm0aX6e^N1)lrjj&9O>4f}N+HS1Av)0%tR3OUV$d zQavrz>REZ#Txm#a@MvJ4w90bzS-`zrJ~}(%y(PXSzvaKZ$6-aVk4BHC!=b~uXJe9c zo=}|F%kXBMvW}ThPe0|&Wa+TDn_hNrRHSpRVysHmdTwfOB5!nU^l1`X;c1$(j6Ezp zY&|ADwwdLsVkxhvJT4zwTrc((zo(K*9xQ~hz-W|TY0XK=$<<5IYtoC|DC`dZTKKj5 zYvP*^m?c<ZxS`BMWaR21~)jsg_i1&l5&bfK~FylJhP6F$v^`X0MlFf_lmFL;Qo)Il+c{DROOoT0>+^Q@C>=Ij{m*Dag9tnDwuk0W0aA3GnLu0wTibeHJo z=>w}-7d;kzssru3>}*;M(P`-G&v1DjnA!4 z4>K1cTQwxS#L^^Z;u+#*Vs4`UU^-plcMOBUwWn@@Wp;NN7sC$ySZZvqA__^N{U!dF=J3{#CqMB?L>m6e# z%%Kjj7qdH}c}AVLk}nA;T1k`0c%@vV)-pIVc2iVySO@QSKVRvSK9ilxeN9Boasuu( z0JVN5Q77Z%dI^`3%A&Xx9~M?g9%E>R!HNYQ${()I4h?{8x5JYfJBnKu!@fA5PclzD zky}%)vUN~DpI(*j4XIOBBUk^PWSUw|yXIV^xPMMwlAaum;qB^l^%SW{t0@0WQjuD9 zreR!LY%npb^{WYLQL08p^InFZuZQ-1{^|FvEyN0hHC#x~Jk?d=Bqq_O3_o&*C=Xn~ZCTYx7-8wk9Wd1XRh@?$kV0*-mFx znj9B1*Von;?_19;*DAVNJkFneUVK{0eJV1oqpU~o6YtE=eN}PQR&V)so}aiU+PB+d zp@N_wkwWlKD&0$K1PWbC8$oP`tv>PLxj(FYP=!WzEjD$IGBylC5Q+L2-dCo#p|=7M z37#0H7`8W)o7w7}(V@tU`2DQ?%ap>D-i5S<4YW53?v=yN#`ZjvoQ2IZ@Y zXWy{HLl)AxcPG3(-j6morGxq|P+-V>{=QC5>OIzTAb2+toDreLuTY@Qa&sD6h^|{WpWK>~CEvNrxzqbpri^*vng-B|!q|r-Yz&OpU(MQgGjm zZ??Xere&g|Ak-oD^S4Rm2xh-6)Hg{V`S#E1j@Ss8`izh2L>t&)zGW)2_FuUH7b1B@J#K@6TXzD{plxW!2i-W;3RyEKQ8 zcyhM^6eeiWOC{ch#W)^>a;^*dd{3XvVGOD`r^YwFbCrBdap|PZcQ>j#rNBUbZS!<6FX=jGvpvR{y|ZDSsu~8DtZU zeS=v_nWSnrRZ>%--AX2d-spi29e(Z=ec4#%Zg_E5YiJA6Z#KX?Xo9QK*Z)TTe}VNs z0G^qP@n5Hk{yN0+U(o*_2Xg+!y_x?UqWM3e|342L{okBM^EzWO_$?^w_n#d+q2jbt zUpupPZAAURa3_u9Xl;l->B=;tdYViPho9DrlfNgI z^o4(R(nk_8e&W-!i`a4B-_f3ZI{1;hdOm)Q)x3slZTBlJyN#*koz=J1l~3}UtBCf> znUw`pWH;}nZ1=4}YayQ`!k!K$hKWhQeMkP*xY6{Fyz_{7BI6<(WpOdtMA5q0#rw1?dvk z#PivUk5^C&t40wE@xN8(t)5qu1J-|4w%Q`8Ew_Xc*n1kLk#E(~hQ_H8VD>e=e4gFy2HU$=YNVx8UJzJ~Bf@fcH!G$Ww7e z!zx8SN*(3NClrzx_#I~e6YpGh=c9?lH+Of_8sPkRPFZ9$+op=Cd}LcB@uT&dym>Zh zm4jNE6zVmMvOa+iJhT58v`#Y5)q7IhS}n1Bn_&2(TjX?@?ht{mTWqu0J*9=J|Gi;r za43b0hIE_VQT<^E{VV`Wx?z@<` z7vEz{LItL~2sRO2eh4_d93v{K5rVpV7{1QzG(eC9(U|onhx34=K%R(!JNci&NMLc0 z)3yn1p=7Om*#*2)ocfttW>(Q~PUH!E^3VulIb|Lxj?bpbH1(KPTmRjYh{%b^fR4DN zZc68)!o4ZI=V!c`HFguYkuyUfZZ&NLB^${56}wd|pfhVY#SO}hGwHcHWY)Ec=EU$I z>+DWL$3ywvp7s?jM((!Ck+#${Gj=7t9g9hA$M^$ko!=B-XjR1q zS?jmK?q7lilJOnJ!;@q)>WSIIl#-dw1INEH_kwre8Y_`-Hcbez#vvT_vN2@yaEOyM zj-I$_1zC0!IRf7{Gz(KC%v_U0+nJeY**s+PH055NHK^?+qTx|Yg%jwLafMez#Y(M` zCX&0Z1~a71r?-JleHloRA#Cw(RUFs*Y@5-?V1%ZS=GU^+l7O6R0J?6-*cV90N4r~` ziSWcn%>fb(8OKqjZVCX^u}Ka3eOSLv`ih@Hj5m_l!=L!VoFeK8XABeGZiN{(jof3S z4p49J^uB*>ao<0C%XD|>D0ch%~uRjc>jRr|M}_TED?y_*a7!V1rc z`h5zJn-P!Ye(*#U0e$B5?a%!MnOrRHMJ8fQy_72+SyTWSY|_<=EhSVU7s6*O%}G=R z@3sGcHW#hx%MoCryQM0$A~|D1A~Q!KnM+-32$XF!YCR4=NBt>65$30G#;+O5f!e); z_pxzrZw(N%Z*LxSBdB4e6MWP+>+=0Ly!co{Fy}ns@KRpzY#G-8*`yA?FuV%<`t${U z3_5E%9qM&adsJ0wxf~6^Rddd8?%Os6zw~*SH!QN|IT>wITHy>h43C#+jSR}DZCkgT zhNi|4*wPE^9SGQ;kruio_bx` z8_0$&1@m@XEkmTCYlWirC|^J&s1|Aa3~?e?N4kMYiV|ncC6p}A16f9-0YgTjA=`zP zhUdIKW${8v7(jl2TL!F4+_PkmCe)r{VYW!eoQ6{MI-Jb9t?aj>=`Pq}s5M@a2NYcX1HJOD`yX^&npuqAPfB7s&F4cgA z>~)JL$6+vs0=JgSrec(g+A?jYnoyU!WBzqalm#RqEKOFFJ}*`$VjRyf#-;YM6!VL0 ziHwQSorI=S>0H+|#j}Q@Cku_I92yw3%bX}-fiF6#=|!)Zpr^JXgV<~g zMfGdmyc_8?yA5Rj{opG0)8ZsklMjPCX%^ z!c=Ostub96IOY|5{+FpTD*xa?Fy^(izf?;G@JKycP5G*JmRVVx1TM%_ zn+Mv|x9!LhA&rlBPqNM=;tMjfWJaAbWTu(&_4s2(5nu2y!~n5aX0Gan9RcT&Rys=J z0_)8&@!cs+#3U}>{GXjnzdKDBGLugDE>2SqW}ZYFsZL8SKKk?&H#h+5R}9Z5>p_OVBalYJH0U z-T3e!=8K98zUyCV!H1Y?aWcHr2vohSzNDzUpy?R?i}lEc)IJmfh+`_lEk1R~mkR2O zsHc%_oy$s_z-6xi%w>geTnDSCq8ugcpW2t(9!?vQemu+(7?`4EhpF#Sl90)qsX)0Z1S+jSlJ(ZwlV#8rz zaxRohkm!J{7_kQ*(;r2rrpiGydN2M8c@7zi@Yip<4MECYW-lCDy4hr<5U9y$5hd}g zMDCoNs@p_O6MqQnzBoaeVLENe9P>eEICrl>g7xwc#@(ldK$lcJEMc*m@YjogND|w; zap>=PTNU8yt`H^3_!+IZ1?#0L+_6`2qWs&SW8v7Gh&;2X~|~&{0+{VzsGyA9g@xR zKrkNtr_?U7F8GJJO*9=ifK#sub7Jm7Tw;&xQpcpbO~NY{OOJ^2Jj4~gv(hC9onxN& zzKI=V;3t*9WZ_D85DRJMhQTsWmWvZqGVzi7tS@Y**~WOTf5u!gIo=ncH{sSXR#&oC z2FsQ4pl1^#OA`=BJ+cQk#jat=H8Hed((h=i*Scb@o;@>v#Ee?7rdH6+(_iF$`h;`M z-*qajpCa}CD=qs7>Ex(VRaVKralTq!CWY%AilI zqC2K7(}zt-oUhrh`s5r)t<7g1Hp(;LMux(zVKd@|o6?=H8Pe*~X&nT|ECS4S3$7|E zz^uJDE8_Yk0N}`>L*zd=oUX+<9n#{!CJBc8pmt3I5wcL*$C${n&o_dXMFe$Ps@--| z7j^%rMMXr0g$Xvq8d+d9XpdS*{-L12^yqjzBzI^fK9ud2XdOXuu8aiDi&nRx$V}g- zhBrDW(}&E~0BcwN0V(o6lXfx(lQ1TUmB7^phgoh0sc1nabw7JMJO>;=XFqv1{3|rI zl|F`tfQ#Udq0TI*bM~1%?W&PX@uyfb>aQrS2k=pT|3jS2;`6i;V`@OGT41no zjQrzwjzN#Gk*u45?rzM&T2ozSN^%n=dXQ)e8jA-BAM8$tqj>p)V$*H|%75l{6fxU8 z6G1liCO$yQc4!SE2XZtnplB)|k%N^NrAl}!GZv|`U=;jVl|mZHumvTz_GH18_=ZS0 z-_M}CbV#DjJ5AfLSu~Y1Sw!X5;Aey`TB)7G>jc(tjZ^Pnv;1B!H=pX`(HBSq6kOT= zur2*p@`;I+p5b2^DaL=`9safB{C}ZcF#bD?;XiN}|F&oSFT4xr|Hixc2WakJT#f%q zsQVv*7yk#a#lNvB{s))^G~IuLSulUYEExZZz59k){FitC4@qste-SYLFL=C4jlZ_* ztcaawRPOV@a5RbqBzWmf96wt@)d@n1AiRu<;Vs)&YygjRZxXl?8Ddop#uRK6@>Iwi z-WeWuro$gZn_+ZnV z(Fmg$D-DxjUAZGgoZO*HURZO?0JOJ&nmg$)voC_GP4g8<(nq(vPYO4~FOBzx>>g@Z zl*qyRkJgeozDrP>@T_Zl>`mo6ay;)Tz03D=sxNdU-yyi9N@PWfc_7thas;3}Qd_-a*SwdgmsW^)S7G&!aI8 zeNECz8?GxI2J)o?c@-teSht-L5S%}1%gCEJkn-UCW(A9jw$N8hUbn(0nY9~S7F=sd zYj|Pt2{5wPlNzO*)k4WSEaeET&_B2x_kH(bcIS#S$%!_5!XRlmI6LKJURSkvLTQFm z$k?d3{p#K9Fqgw1h$9enCQ(0ObispPH#d7tB@2rL9uf)M8v9C!$1B26Y6o_-=9G!JC(~|IJ2mM z48JRa9Bf0Tfkr$+494?tb?;kFW<>tmtABn-RMC({5mR$k%E{gg2xgNjASiX7DNC-j z`c*z(F#r^YdR&WD=)V!>m0_PesL}lLnpjd_FPbevW1W3{lvpe>C z53^`f(S_$@nsN=>RbE;5UHgtskS&R(!zTTCAG-d~T>`oT*R4VFoY9=yEQ>0!%M1Riq`Whzq4{<1GE{wa|JZ6lfq`Z35M zge5~rw@AKLt^h(^=)FM@h|3Zt(yTqrrMQFjq5B`O82s|T^OP$bz6sKIl-)fTElRs1$p66GZ#yg^YWZ>A!;A=_ z47#@ax}yUoDR{*I2;!~6KdDJeqL_ns##u&G&4GXogtxPoH{do2q5v=>%DXwE)>lrn zba{Ao;Yz7ex;RTMA==0ZhgF&du@KLN+h!k5YFn&0*8NQKS1mH2!9A zl_dJ>sz6Gj`+$B}Ej-Dx9ka9Xgxx(ef=rLH?gaT`Y9`7KmvuQmw9Dj)2yNIctp`Cr!PX?Jf!Uh&v9@L8Jtj?Q6C*`zrIpwqu|P+ zvtmYf8a@smJvyRvAUemeLJ(?S{7udIFwj27$f65n-WN3utu>s7+1yyt#sH0M<&D@3*)H zO{G6u90Tm9NHA?F+{c*eS(Krq~ zLJ@`NQ2UHDEFLBiKV!#+xMZqgUKq8FiefD*kAyO4MAF*F^2<|i?2Q8#_T&P<-X=Aa zycb`|7J*q!PI#S(o7}#BoI9t0{tpiRUy2h3R(6j6@Vfuff&D{G^8e)E{~Z$bpAP#!N>yyU%(8LDhuwy1PSv~1~J$*((=ZnaLvAk3(s;ULudh{uc zlvY<865*GRjB>TfN$n$9(*1C{!^t#WCMTIm!V4~Ruv9BeK!g2duL%|7l5Yd~QI#$; zu@?S2OqCY8tX`}S-_FvERS{iAqJ`fEly9fI!dTt;1AW078eL|hjob{M(<^+P%*g!R zaugv!sgSamDPj@_Tra<}l-@_j0O%gZuG=c`}FAx7$H! zPD^bvTvv&rmEY3!n+vJCSwf5!4jyDKs47jnJM9xqqmzGiN2H>N=<6^TOhQ5f6p;zn zNDw55=%I8_N!#{Ubug;E?{i2{C>g@yQc2=26U@^mG&q*WXV0~EiI5)QxKc|E$lca5 z6Qkj~5_|+*iExBtmdT%e#(BxgR~)2sf?MEs6u!qPhS(({zvA#7wey7Z7RtUHQUII_ zZ-CXC4+(+PB|=|=fMj(}=glEJe&+%i`2kit$lFuH%YUwaMf}9;#kVoCX?v!o)GB9; zjJhISK#_iPV9G>mX9Dpiw1!X^3;%d>X--sUV^8-8|2_1Mb6oCsTx1c9i0Q~03;~f# zo&8PHN=PllCPES<<@Gz-OAqv?4*zSk-u43YKT-5(e}|)*c35!F#kIV}ze|HZ*t`JG z^rayySHiDoUMdB4bz9;0E(Hej-G6m$V8n*?VNRhgXJRbHspP9`w$5ptMP>ss&k4wW zeXr5?r2e@|%C0e}4HJXSY_rGpiVvp~LJ-sRPsYmmle@Da!-J9Rf zu+QyjT@8^&@4Br6aH4mu*ocSnpC7DZ_7vVv4~B0jc%671xx?jBgTEgdkr9QI1n}G% z8aNberGPPWF8{6Sqt*%0XoIhJ#*JNZc^C7I@xaHorlRVL65=%_8gl|d{Rrd&fr$4{44N7zmIw?Z z0K^~*?VPrJO#2B(zSaDC;2d@R89HeP)e%-;5V#YV-V+SL11#E!ph_qucslHwMi5It zDY*9h;+(xv;T-xbpa@vi8a!stqJe{A}b_j45oiHp;x_?Nl#fms~31##B|G zXC5!&&gmdwGOggrq~LxjBZHfsU|(G4DqY)U%M@!2qhqTQ<^=3pI3%Gg!liB2w#c*@ z+Xk!&(R>A#QBhM8Q8%&ImaOfJfl5nmg&)%i!$j5oz4EfAh$y0fNNveDj#+V5$6iEE z#ubi$$PUUhZYIJ-fr-5$lEh`iG-m(_F*&@%qg>=UeZdg<8>uK6oKba3N5m`_mO`2_ zHWx)Iduu?E$4uLLfz|Q8PoUya+nSSyX5HoEpy=)sZS~z!D>#Q$Fv1&xLOq2wUefc) z)Rz+Lj@&+Jb5sUB$BVVDQC*zbtOtNNzB>gmlhXtLBPy#5zbX1NlN=ELBoT1z_OR#C z%51J~v^NaTDbYZz@aon8(oZ13Wt1Z9!G+8c%yfY_+(6a)VXiie(!i7|9&4^<%oyg% zcyQs`_;V3qO(hi4##>>QXtm)v5537FAJQfqF&EIB*g?43v14;!6Bp2LBuK~%MJ$Md zMM5}3A{>wxcKkv+Q$Vz~zj!BOi)t6g*opL=ryW|3a%--|6-$R9yY-uz#EKK2-T=Wi zldjB0Ew_tWm=ZLX-E8u%;#=`lVQ#Y$C9V!7D|Ek8VLH5ZeyuPzvR7C z+fBdgt7!cLQM?I7Dv!Ye04K4E$gcf_v2H|P4J>;@$rX;fLir^Fa_96kcd(a!E4nhQ};DQb!VOp7{&}?x6-!GJO!4bysB`=ygOj zu?mgd_l3a~=^dbp3{De9-Uw+O$Yy3ke>6>l+MM`e;fwt-u7V6+pT$kmNuMw|GSFZ} zm3g6|*_fO8C;`nT%{wRm5twt!3^nYqOHT5k7D(rwmX`eaXX!4%2dr)}7B0}5oa>mo zE&zvwZceGiihMZ(uhqN`cq%G;*jgtJQ6qb0iP@Ny_`3-|3h7BD# zTkJ=_)W-?kh-ICuC$?K9F{>g!czmrfQfsDeQdg2R>CnEGIG-bCYcO?D<;lni*Q{G{ zHc2&X2tWP*czkG|G z3b+Lrt`K>xVpAv&5e=zmM=jeovVwd3$;bUw3@EL0d^Ad2Zg^_MY@uWRrU91xf>^mABon+x42&*iI-CWn0t_q-u&R~$WeagVtg)r?b^Jx%xb_ou`Z zjDDW5oRDZiw+Pl0Bsf@YwuO-$2vv)HRT{HWOe zcLkv0N*Q;VwOW>vuBvRQYTv~6v0bKD-nWi5|EGkFw2e$-;p!Yd3c?~{2Qa4!jVB*& z_lx>f!h)uYM>G6zfx6|~Ar<+!I?gnh7#kt~r+%$S$V=PVsY z*5J-vf{!uiKm3orjvXuj%5)cmZZ(t~KtTuHKc?n}r7eDf6x3Zf6Sv1U)eD*gamK(E zvOjFTm})3_Q3QE7gsc5J+*+Ar`1GJ3P)_6V{pWUIAPp|S`N zRt5<>ZdRn@AncUbOpa0rc*l6*XH~-I8Ocf-1cvsdK~N-%_dX%K{nHS*DgIdFsSZ7f z#Rd{3HfJV)PBts~+6n3fgbhb=O70_RO9-=^JtDDvBkD+(`D=>MErzl!{5_Eytb#(N zN?}YH!6$Pfo~;&VRW{`vR#Z(0rnI1CM01sdC&zbelXxkPxGvS3geQ6v-+)d$e_v_T zof93yo>)UJsbDE4an$$)!5KdGm_(5TAHqe8RAux{uz5aiF{k0Ck;)Hl1?Dp_M(<0!-89L^f&j#MnUqJnIqr{bd+t68^r zyc`JZ+4AfDi+v_X@)g;oLR*4}T=p)A^)6}5$ta4xIi6**N4)f&rcZAlhBc9Sxvt(T zFvbj6MO~uweZE+z0@_ZF4>5i>D2nrseag^ca{5ub>)?n&S;a)G3D?39sR)KmBV+B? z6X(FRk@fopV^|J;UlaLZXCXoSReFaSUAE)AReuo6*QYZj9hH@KHIXU2gRr{UxLtPz zg9-vZEU2R($ABN*Q)doGEg5Njw)L9hWZwyN<1ht9U3~0K9GSuNct>v6n*p9@7i&)F zB0c%(jL1Cm=czT4{P(Vtl2;8WTpF`4Eh(NkMZT>#%Nb3@wR~<({0K4xga{F?;{=EJ z*lr*73WhZ9U$^jpCbAbV zM$xwkX|!#A+<$6--7<#d4V$ z2}OID4-#g;SGY;{I)1%hbEJ6iCzgFdjjgjqo?jglb{!UT!RC%@EhFgb>FHD1)61b| z+q$aE9UgP`9iB|8mL$n#@zy$f`y@IcYJ$Sl$oZp@bj00gEK|{BGPQ)Yt?7qsMizh4 zdS4QW^!_FhekrbYc9Kj3^M>7?l#&eOM>Ok)-itC9P zVW_Fe%y2fF1#sI37(_CZBAV1~C)uEgh+HoZQ&L4{Qu9O?*@0?~>>LP|5}ch}+0rC& zCMGz-utK#($9n`08^vNTxA15YG>bVzGF_`arcJ)BOD;KdcMz7%<= zl=vJ4tw{pc%|D6`spq3Ip6kn{d9vllu}|w-jVGr^6lsMpr)sq&=1q+^F3AwR$+hAg zkv>a-OUzs`B7aLN2j0*YTVL8Qf7QwQh z_>1JVFtB;O)Cae&q+{I-<<^V6*iVnkqz4b>^|7(orESLAc2c4(7^XkgCa;X8qb&uj zngA6OcYZak?_W;VfQ=))YqM36pcuf@Bq*8^OD63;vAftz~<1WOrg*OBm zR?2v#E`>gbDZF4$l;@uWEg9Ekfz6$;1#8#iKELA>5p@I~JL=gTNMeQL0B6Wv@K>yr z&^Ca$-CriQYsNqWZjr5}+$#7B7zZDiTY-+~bW0q_7f4d~A&zw~mS8zBMa1)}5@(L8 z{Ov;=8QOQFX|U&PSPorMaPf6A!{9@Mpe93d_;Ypku(Mulad`Jmh%{jTEFl6UP+^;x z{cz?$_ULr~Skc@F<4@?Aene=baDsC50Hx$O#<<_!Q|M*nBmmu*J1Pm(sxk@uS#(bY zKXX)&VgBK=udkQ(V{bKxFm1g^V-62oECDE7_z?$6FCk8{?!pR1Y`8vH1Pr+p#Y0D! zMx3U}<`tVeYrVDV2EN>fC0)l)9y9mJzRNM5F`g+8-uw5lrE`1~3OLEXa*k{vi-$6! zTv$*dc%e-PTS_(DtQ?O5Gn3EThLt3WtPBk8@0WOIUfi2|%+@Z)fq(E-SNMY)m^=kp zoNL}5aPo%9Neo$s$tj>PR1QGn;rXvw1mI33-k~m}ZPYHzFPtw*sx%ctI$Yruu5@1} zz0Ht_)WJtuOXfXf)rK{Y6S+7VDkT9%6QgUtQzM+6YC*UtX#*7IIL#MqZ%Z~pmEExL z8ucVdWWpum?mO;wCf-$Kbj)Ne&183ZT-?EI@*EB;E|;415-oeC@)bPKv%c_H4tmlO z2hazQnPaise*4N@Y;|n8JeuznRnmQMqDqMPJ$`DmH3lY{j~g#iIHbCNQ*k7JT5vUW zR5b+HHd$d=@gKgqh6owUHmq3_-8Ej)r)vcC;ngy|-BnrRx4P_Xy{4iNPoD0w`>N=Q z2V$LR3ZaDki3qRObJ8CzAue^eXa9Ka>8eaGdf%~LdfO`KlER!BGgHE!(aOn9A35eR zr5|^+Y6Pt|Xi7J*BTI$@O@?W7UB_ z1AebdB+{0BGr4Hm))GIwCf_ZC*2z6LANbHu>|Xqc&+DDC1v3Uy{^R?IrK{rZ=Wz%sgev)Q$y$I}CbU#077L6UzWoaB`7 z9AK3Yo+&IaL{;H+JZXOa_w7TpJ1{ZaN%yF1AjRSlU>HA)8H31>XbQLFSH5g#EICE&>w>8tB1Crpyc+cw zCuBe};)LY^10C1baCX*8l_-cNr2fiaUeebeiW3ogLU#S6Zb?M~y{D#!k@(1R>JPuurKD?S+qz??QGjQ*LHT3bm(OF>Yt;dW z*!6plmb9mDMO$(ha=fWw5eR|vd2FNgx-}B5ww+k=R zg)+0OR*Qu~nR6S;xrnUnCay$ZzC_8+y1g#+lfF`~ZXPhL@prq!%A*=bSQvLSD>E{# zhmwM|j0Kqeg>zZgZUtZ6BC=0Me_vWI;MZX`pPIP^fgC5~d3hHRv{!?5WFMoOK}x$6-k#;jgJn64f?7Hp;9H z=JfY(2xBI_HIh|fx&MKtT3de+u$AS0MUFbQj}2ZK-a)?$u<*9tHmmz#1T(i(S7YpG z$?7K?wW@z6o^9*KDRPLtQD89~v>}-iP!i1s8RB8XZY(B)0TkX6&{?3Bx^fgD%1mxF z16d+4z!l=*ouZtvpINzgLrBrks=>TX*U+Xd-67A$kuhwt>7?D<;;1gPJTdlMzjkdN zMTrLU{aIV-jPLq+VhxTTcg2Lc8Z;w!}Svy>l=at zV?KQH+|u@@6R>hvBdcdJr4+_J4SUpEkun>;ki5R-EaO*kT{t#p(7v_z^;cegL-y>ANWMf z#Zg94)t;VB`5roFuc=chIZj+(Az|S>TD``{kD6`@aRjSStX*$Zk2`lTa`O~O)DF!jWe6{ z_2$FRQB&1Cl}AF^@;d)7uapm6f4Ovz@tU8CQqFopBM7%wRR^)r}DH+ zdbx?FvqFN%2TojmWpGkUj*Nrku4YE=P_f=V%MA_BJ8pxY*0xPs_9F!G$vh38LisJc zr<*Gr3i$w~W10&mp*rFs(pczU3qu~JF(_p)D3GREQUJ$CM}}!w5Ys(Z_JyucS1350 zAu?5{oMrTSz5c`GkX=RF_(ffVb450vU41*xS2?5u##B zU9XukSt5DG_m>oH4|{D{oB)Fs2}mvQ%2teVGs#FCfWBOYRh;fGM=vCyUAiqbiA)9w-qQyw}+o6DlWVf(>;te6*3dkQCri{dp{Pl1Zn>9nn(nt>mxks6Jy^Sh!&CA$z{zvVW{Ck@-G%75L0?}9!lVjxv(dPYzIe^Qi(}qZ> zEUsSx9Lq|B&(_|^=EJt|UfhJ}#TjTJ)FDFOG*WQB=b#u!VNrd|qcsejf!C^s0PXx=Lw0m9l0~|bbYs5L3{cCt2Pr&N85wNn z@Xb6!Hg3zmi9F%2LT$Ie1?DcNd6?Hq(Qz4l((@5YzPA<|PO@`o7UhjxEE*tj&tP%S zzizC2gC|X6`nyR_Ir-muL{9O5RK8f>vc3YT?X6l6HGBn=$!HaT#rpZvGgwe%2w=!r zF}%orRuKIShpzasN$ub8z|`=gYJcDBmUI}ne~*{9>A zrIX$CLML;+<@|B3`z)_=wP+2V=aYXeg_G{>?XL7dNy*IjYgNZ)sCAwdn(FAs#2B%O zvvks_2oCcfIH<70q(B$!#+4pqr^G)sYo_o@iUkJHjK~>ij}{p_Z=e%*1|Bg}lJ+56 z92Fl;UZm-=mLX}+l8@?$IeoBgX6pLt&f!*u&i>4E^mFW6EKF_V?l$I`4nD^^dtAE3 zgrf{O>;Wj2Lehxbx(RfYupN2qSP_SJ(BxMK>exbjPTyss&w_uBQ`_EPJ({{NiM+x! znQLRF6`0M5F$PM}TN6rbez7?*q8G8WHBYilG!HkglD`ool)>@R`IS}&pkB%v2QusL zN=&wbPK=3LVoIvb^oy`g7P7zETP}&iO{=z;7h$~`)@o-63_yn}orO@*xqJyjshphK z`n{0Ckiprugv0UD7k8euseaJoK_I0h|*%##{R9<~=JfJhyC&Q8uWn z*{W(Q8mpQ&F1^NHEThVSQfAS@!99=xs@DZ#KXJZV@Mr;=63cQG$Ae8-_r+O(TJUKl zr2!K&jlvN^(So}niGOoBQt|^f8=Ggh;sQ@2Uz>W-sg$4n_i_%ZaR&=f7e}_C$s7udPud%DptDo4j1I4S__{hwjMugZR)RRX!sg# z_|xQ9J&x>-`&+H-Qdy1qR-H~1v1!RB9b80`D0OFw|L*87^Gdfuw_hABotj`Hb+fx!a`ILsd&y6 zRwCuEI593v!KmPd5C$mCAoz{>%TmmvJMWTjmL@|Tp8S-5_j5m++Iq07?oIdgjMe2K za8Ff>*4gRgX4Pt2ZKr1^clEg*yguBv*@5Vo0!0ArzV*f9L78YBehYOAQNxnwR5BUJ zxEchK{(LXNqJfI8mD`hM?GF*#!b>L{x#&Q#zb+?4lQ%w^S~*Q}lr3K4NOjE$AuYZ% zo}|$jT}y~%<)CnA^U88_jn|n4eS=CcCe(*F^eKzrXcUnIl3wQ7r*z7I4*E;1Z*c>E zRd#}6i?f#e=wq?7=`??Rnn`Zg@_7}__r>dVr#7bAz*B!ZICKnp@oGh@^>H+y`c(in z|9oVm>_zbK%Ua4$=E8R;NBLuWQwYHZd<2Yq1LDDq$6K%EqMC0XQr#VntY~pmkw(Qt z-Iibp5p6)5+jM+2Darm7%7&l1l|KXr!=5D2hkjQa5Q32Hg*Au+e)ThtSNq(Wk+GUScNnIh*Kt+*j*>FQRa*@tO6mw0EhMzs+vQBh_I)K@iTKiYj$Shjw z=ykz7jZtN0u#u`oM55+P;~_l`L(M|4DhShi@6Ba)AgQeckTdAldrGZ7KgiGLGs^(- zGua;M5WjWa#NX!)x#Sth^`pcXN)zB0(}ZuiLf8q?&~Z$R7XH~P7+UfdbwZ3-k3u)3 zCs!P(mIYF+9Rpv6HTDAq-oME8!qLmY!&E-9WInH3t03^bqCY=AcffKfO}?W0PW5|qT1Io?7eSup=z$k%amv;Jb`S~S$OQal2yRzXYGq;95Y zR-ZSlmbgVv4(jwEQ;MR}R#nxM6!b6{*8&BHG`t^Dp=!=9o2I&_oa|FZ7xPd#V5-J? znTs6=FSngkoDUOmunOHWV7YTug-6HR9L1uZf)NW^AsaobEy_qSVd|if^bLG;=LKnp zBe!#hca#b9S86S!WD=0+HJ_jp75Er?|B<6G3KFV8GoO$KZ3aid6QmlQdW4ovX~t`< zzA?WsX~^8om_IsXB+p6t=CBSR?oJ$N?p|T^nr!2uv-_AhoTob$v3k|xbbOwOF5u09 zsg6X*WYpi~Nmt@Z_yD|_>bxkc-1skYw5NGV?PYGRd%C9oUG@8WeM7aQFF103 zua|gG6MPuk{4}nAcZUJVfR5Jx+%;W%qn)6tJK(sxEnsZy*FXyw6qKnQG_BVTbn!oTpo zA}eWa%z+($FR3gLMf}1eWl0N==U^kZPIT<6lZfD=CB30=7oy*d5ChCZmI~X@gZDDg zuf?992?_v`DNF(=ab%ou6qHy< zRC*5-nB*76$lPbPz?I%_fm-G*w*cMeEde^1VC}NuKiXHiLTaZzh1wwR6{zl-W0pPm z_e=|h;CLyd3wthJP91S+R4=!sV^PC)um?)jE2#*2;?`sdJ9cn3yNIWKcp#vz;Lp^f z)#Iya;oscKY5c5oN2R^Gq?x5H7k3>?^vc#VwazoK4X=k|+|uTRkF1WdDoI?9duPgr z@%vF+AFn09m-(YFN*EAeLiOP@`e-u8_AL3tni=~5 zy&=}Xi{%9w`dLUo;O8h&mCH&|HA(A@MMqjYF%d|=U_?gY)#N#0pMa0AO(El1Eb|p8 z{3n_QSUvRf@n1FAQi0QecpHjf7j})%wDe*G6=2gq{{;Tx4<+=yHF?|{qjM`A{j4fp zPQVYT-_XgYhk+0fT+~M+xn0ETQjNr}a(n6)$Ea~rQ!$&x>twW$upk4Vn#^|ClRZF) zKzqoPM^S1ee3zSw{8hi|k5Y;oMqP%(fKht?mDiF-u`WnrfRD}uwI>v;>$f3O-D>hn zsKte0jD)geW$(5f!}n>sHL@G-a&cW8VB%5#BR=hd2Q0Rg`XrqddYWaJNR1y=5qpEFu#_ zX53ogr9-2JDiEW?hR)7p21q#D?Ml=5kwn3fmMLHR`tTHOi>~bKd_D^AVqf$GR9Ybi z32L^J6qO+qQ%!_AOdYM3$dGvG4~8BJ*|aingy<%1Xzuo+DN6fx`aOPs~q>J1J{qZ*XQ-iRM^X9z?occgAQ^ZGdz}70@rl z-+Jkou;UEg%)4}ZP+8^GWa9?T_1gIz>usB_T+)RP$Y-SRrOLE*iEHFiYQtkTROi^% zqA*I4EKOaH4AuppPEigCHlZ8d`0aYDmdlpwv^+NQI^zb-Pl_p83MojM3UmXZ6*HWh z3|Xg$kO0y672x5fHI^>fi7Hierwg++1F5r(BH>z}4J|e;nZpuW9RS?h2S!1rrhFeK zfoNJzs3>j(OE_+r1f6K2X`#xGx|)mjOq>@G@2#fNMT7?QIHKFXcj)A41Z;~7w;Sts z-2?>^w0yQL3Zc>BHMno_DbH4hGT8&4!|jLzTd0b=!v&~B#!$lIvZ&o#CTvbC3|6i_ zW_G@SPHw#v-KWR$TD3=%6&@b>G?r>IEJCw5(&Z_0YtV;WUA-+|W_%LJ%k{~(Z7oj` zIy;~w)v}I+v<4;?jTOJg_I6slw;boop`F`wv$x@)DN&{ljZ5I7f3I$xEVxO*w&Vne0;c|WOuXN5_pT@{**-~|A{g|CDoPkOz%?-*eEJA%ujB+^yM0wdv zv2)_+Eq|R(^CRTRhOC`9v*W79!(^QZRBHrvS=Lkf#*duAfcGv#o{wAGzj>2%^H&sa zUD+=Kc3M6NYnH|rIT_7Bp-^C(e@9EJyc{nSneXiOE((^`^n87_4VG<_@$EGN;a7q7 zgWG=L=ZPeBg3CU)vJWCrwJ$<#6C?oD)+etl?>*W~XWdmCbAussA)Dr3mN=M$WK{hq zF)xq;m0j6sMueSbhXH1x1#QPGn}J@3L;{sF9RjZ8;5Jc#8MfnKC?yAB9{@EEqbYtL zS}qc@)A<@)%a_EWmCE*Q-ANns+K+xu;cc=#JNC*Je^7jke%Tle*_@!F#PPAQgu3W- zb@GeIJ<}{L+Zsi38jyWaR;T&`=gC9AJ04&vz;ujgkWPvSe-mfP&7CC;9=kmy?t>0N zmy5xIQ@W$2Yh@Y2+Y9GB5(MgyE-K+dxti>m&ny)*sxT;WAPX2wg~c*sH1`(8pTtgF z-VLgdr=&eS>+vYhdkz0JBT*ZiS59+h19M)#QU6Xwv$R&_wqBy6wXE2nyK3CJ&@|%E z-c)s+l-F*fg{#H0bDrcL{oX>whQ~Z0Wswmxntp*IDlDkK3;uydwK$DyT}5A3W7igz3zG%Xa!?^$#KD{1kFGOcic(x%PX5u<_O@ zX(*}rF#~tj@dUVIS=OOAoso|W1B*!5+eGu7skIHz=V+9J02HN1_uJL`)uPdKECUk1 zXb!w2T)Do_Z5wR&k#W=pQr~;*zKBq!z<=W^gsy!2Z54nRJ6e@f~}9bmsTdztiJxH zj%%X3C;b|$U~|U_4;ibZ%>X#k*G8!`by;>{$zuL5*3LOPlW*(uv2EM7I<{@ww%tj` zwrx94Y}+=vV|LV;^zXj+ox5h%%$h%D)jC!CoPDZlpR7~QlT^0*535uvAuW{U9sO<{a=$2zlbXrrv<+e*3hD5~bKmGJ7vQNtLh6A5)vDU#-)EO<{mIL&s|k^Cty5{a5N|XYT>_PKozC z(Q$D~VCh5>dQ{rZpvMZq6hB5GzR4#2u=XNz#uN{PhQ!z$HySWI*+kS`!poo&Hdc4NR5F0Zin~0+etVoSXq#MR5!vCrv1|eY4RpU>J~{|N0Zr#JGSU zuQ##tq+}^;7$Wc6V`t;1?soU{dyV32Ic@suD`8uU45Q$a?FL>B`@scjVD6^!eY(M1 z+P=4NDoV0dWF`?r+EBWXpOAW_V|9*-@m>Nn0dxDY*WsWGcMQ;V_?2Jgxw90~=FBok&CmCi-E)1_uBj=6sctKx6 zqJSqq9W%ncGFM^FN<&={M|&SE0RGjlox3E+K*AR}ps=gDQ@obCSG1?yH?e`4d6SE= zVXG6kmR?}ZUy|!hZL^|U4Nd)x=6?H+RY{bzmR0$HzD4iipK?As4T6^I=76Ye>ir04 z$nCn&;$S$skuAa*2KVL*n%#k+u-MS<8;j@H>mUoPy!%COMN6Q^;ggAYQ|g>rIm&ph z1{krIc=w(w44;`hyXZhG1JvG?!ycr68rT~Mpo2r6F-?B_0rd7c#_ zM|_4oMpAGXx?R)}ozAwqNfxux&Rdx{sBCy_)>?72K_gz8`O^jFHI5j|sd&o?l(a-^ zR5D;lKs&Gog)dG&lnU4_j3p-=+eGL%ZkBdF0Mk@Vl>JYcbAmk`rj9#X@K`dwPuOYx z@_mSmhp%MCz|96AN2k9Thlm)cZeaGOg;LR+PG+CEjWmx0Y}K7YY6?7&LLtXs^buRj zeH*x1=PwoG+L1A1iZ?RePOf>G3Xc!?VR!OGs(y0T>Sv@q7`QyF#LUJn%k~6=67Q;o zW60I$ayeWb-gfEWgbB3SKYXnw2heM3n)Z~Aip9vh%h=;qDD*a;Mygn>>MiB@YoaVj zi>v8w`3|Q?)q;xFL5fNQT5W~5sRhK-s7n+qNK|MfcS^DhkRV@phlYh>NVJMB0du$P zXxpOeW7x8AtW$i;?9zZZ#jR_dcGErVrdY%FntkoL{^+CRHhAkn@Xhq#O^gLAl>S~c zuq{*JayDTW+c05!{LvLKAz-aOZIY%Y?kqwiW3Dh}e~U@bqBvcslSyBL!EK_~3We1I z1zLxK$Wi%HnE5$UfR=G6%)3wS4Q`+?$?~(B1myP(I~7B{s+3r?oAOxlre$Jo zp>ycUJ{?>3ijmCm03yrj122NUY9I_c2glEYd=P37s3?rxqKdKX#qRfIw42biX#*k9 zedG^@rKI{-x%&Q?g-*kWATH`XQ8Zy3t_Z(Ao(7cxV(@BTYAuk}6&gxY392!it^G*y zAF-!xeO7U0n`u{_)~6_ZkK1*KbcW0D9`pf!G*%e`2vR2lcDPSwIJK2`OLzRz?fKRm zH}^T?dc3h*U2R`s&&kfz;^oIDG?$;25q{Py+I#B4ZBPP)F+v2$4tZNf4UdXAf6U_& zU6LYB0vp@QOXK;Usu4d7)k{bwc3LTdwEhyj?@k4c{_?_vJ>4(?|SpsRRV66aA zBd*iOn>U?wLc6Msxv6cZp!J;ol=-Y#qo1D1Z=zR#6^c>4S2&6UWL6B=U;yc){U;4OZ*{70W`WY403Ml#(-gif@b3^HKyV@ zBuE;@V@tQ;MH-U9$pUp}eU_u+M@2gb62v;T+sx*(A=Tde=)l|Y;CgVR%HD+Fr}pb{&9^@&>)4unW`|7*wK)Nh zdZ#|OP}4;xr_?%4_YGU%8ohT<07h5>Ry-SbIIM^d(T#ioEAz&hF?Bm!bCuWjw@?NP z)ak=*RAHAjS+i8DFaY-wwwJ_<$xL|AtTVq$z^^x*J=P$v9Z53eqca|sF{jnJB=ik$ zZ^(N7)`aREjx9f_Gv|A!;R3brbt_(NKhJa-9U?(snI-M1OTpu z2ov0vwclJ3tc*|s``B&$$ItuCSVLS*hJDZ0Mq+*N`jJwQm2p%{5P=Lww>S8RoA_4KPJIRg~SH?SP@rjgRi3E;~{a6 zuppW9WAqA zDG`Mi^#7RynVkfi>RAa6`OAi`Bn# zhKhHoX%`4F(p6pQQ)6{E@MNS62*%G6VIhD5^$+pG8&tqr8LJqpu0mr$p+dJO?Z<6G z?XLmfr8`AcUurAUx`OY+k8{RxLmmTiK!({o-I2;R;~DYJEqiWSRZ3`|(*l~VCzslm zKGFmZs*p~5nA&Qx8ajUYam;K4EYwux79{y9G!G6c4Nzi9D-)LZ&lP@qo-KKm}_J1ru^s+eKjAgdxZ zEhR{GJRSS(nERr{=~m|Ax~@TXr~~Q|LU${5vE7?l=1qcDu_J6Lcxs67%0THFsOm0d}SX&@(=T+`1?yZ5!G zlS|{VUy8Y&ES~)bEdk;)g1yhfeD)swBx6?QX0P*mJ*h-7BnP9GBHz37>fp)#w1yo& zCoUE<8N2)|3|9X~1;+U_M$`dA8(Ifr(-J^UMW3i{3&kUC)SAw3J54(bK#5JEbE56f z_PJ~u9vekxCsXO@XHK-|j3%jV^6*x)Dt z^{sX7{Jt_5p|z@;EZGDpTQq6_ZQwt7GjmU~8?G;TXxy{$N1LN`U z3!RgEv0@@+&pD?dr>O+@F=UAx6F97#J}gK}EM3ff>>pf18eE4bQ%Hdvy+T{pGr-;o@A!Rx_uXmn{LQhw(xq)w1KyzE@SQDc zLy_#c-$UeS(^QGA53~FUJ4d}tCwO=iSsBjWb`3M*p|WS$3iXmriP1cr-BYMVSVB}NmNA`W(%e@D`A9K=z~`{ci;r<%j5R5swKqD_`YiOKgP#4Wr?ji|#EhoxZ z$et0TG!a8MAjLDb1}0)g%0o{Hjt(k(Gw_$3ZVUlNrf)5@lMcX%z$hVCGW#QoGXPibdBc*9lXi(jQh070VBnt!M4j zb}SN(sqnk{{Si!RozOfzAlDqrip5ATN?uvn$={9-8#k!)60sAC?v!-LoyFD7;m1r* z;d_iVS+KJLe}J2b=!G>BCpFNLxWw8#I_Bcoy+5xHy0_NTu9$6|wpR|v7|NuPg`(pB z{yp;a?wsfn+uYkpHMmvfE@uQbk)lxTY z$+H5RBC$riH@$aTXgPTNtwwhmVYoY7x`$3dst7q4cjxF1s{y%wf8wPwnm)xgBU=`sPsS*t|{ zcB)zq%5a0!IWkng5Kz6blU6$!ZT_ zJPd_K+}Dhhb*qPhunSyudH9ap#F($D`+zfClVf;$lYXjiEO3b+1UKl9RGs&}1f$%` zb!#>k-7=yrHMCv>y(!s_nXw`bM<6z1-S>`Szj{!nF;#lkg_DGApg9E(lr@rdFRH60);O8Ct!#aybM&rOi)d!; z>T}D{22dAm|0SE1xBLEByCv&ldt%<~+y2RjY)}3kZySvR6!EamsB(YaAEz8w6J~S? zBv+d1F>RXgOGU%6y$f7PASf#DN)3zzB$&tJ`B$D3Ud2_Cl+WW#F*x=pN;j10a)Ia# zbEQcnFT^~xSR_?c19qk%?MBu3~Sl zUu=r(5VMN+^CX=mmO99^!FsFu&ESw2X8gRAt+Lf&!*~!!tbanGzZhpnJ(Lndsf@i zgtA~AWWlqEf<{V^dJW(57$?16fMJO)kn|{Le8LNkD3dT4Izx?cIiMKGgq6`G;ara! zJvl{!3;3^Ed0Z+cPtMIH2YOYi>%wIs>Ul^f>dA>2^~drFs(h{g-1(4lcHfRQ0ZddN z>HnC$9c84|U=xjAWAj~Bj=y0LX0vA2g$12Gx9 zm1dR5*u(8sd72eUdQY+wXO7DEu*)8_fu^KmcFOP{37^Y8oR`j80o@E^7^;sOQyhS0 z{WHwFP=hfH3}Q9SP*vb2CH2aBC=4^|y(sEvDd8`>0%TRd6k~~`7)CcQ$b)sPN(_nt zBQtDs!ZQ8WKFcUM+$Xhw&c*JD!t1w>M#7yTg1WwGR<6M8{;oX58{vMEkYE953?!k3 z3G99nlgWcR#}7b0%Le64BCPKrayofvKs>qs!23Y&)F^!Y1Z6-XM_1*dDMu#!O~%r} z+<{BQ;Kc+3i$^>dx1ILDECfsk=!$eheFpi?&F>aM4gn5z*0wbI_V!~dZvX3<*Wn@> zhYM??Uv*kszpvUlLp9Es&c=o{#c-kLzy&=#<(nbiJP==^s)DCjvCa(xW+uJyp_?Y$Q`y&Xe4=b1r#*Sv8m{sX`khs_r`k5 zRvH6@2Ov^jV16wx(3z#aSXZrMCc?|d!+99h#^1VtV_^{mnaW~97&VXKba82%o`wBX zDGVSPTop3?T}ZlEaf*ha_8eT+CB!Of)OqFy#}=GsFra3BboqVA8WNmjzy}vnQpgXA zCrGn=LaQ46Av^#L_*eB|LN*@L!nZLsiw2KGy`Ch+?5vyf@H}-DEEI*{C1(JHg-$8- z1lXe5gc?Q0J;e@_1RnojhfsmA0o6OA&)v$dzhF1#`(>+1Mbg^_N43Ldx4j51zn6yZ z=-tw&UihawE4<34jziDm+iuSb1Wz;YH;NZcAx!xhy^>v{Wp7JND^m6jhRot@T z$*Kp{6Zml*yloVAgd6y=?j;W^`=VQzU#h!HTQ4KvD_$GIyTMyCXPZ}uAA}IQ2zicq zP=s-%MLX%)b11p|cR+USz7-r|Z&P<|UCHK78%9Ta6voT2sLTzJFmag`<}1_FC+baj z=rUO7lh7TPF4^qEx>>`*2Rw(%u~VGJ^xsWK(~vS4AF%Rl?^!axINrKqq3;%xj%`Ak ziig!gn(3mmxS_L9)^xU`*dg^di9la6J&S^15JuiXD)y;JgGBGwCx()rc!b*VxVd&!!8PpHtQ*wW$)gL?M_*X z&wL3yKAUEsoiox31pE&tguA|n3fZZ(AvK3;)dfjiwTRVY1(~uY&MmxXg!|d_j7e%q zikyPbk}qUc6}4Fdu_BGwN! zzX!fVb8&C57j9zE0^ylfVbE2bKPX@IPg#T$ULmaz7k2mbGw42LraE6=PKlb?qWvz6 zaoIi$Ia$_oI?g8+fzUrkn7kcL)#cqgKR|YymI3VSHz`8m{$iLCW*jfuCBkfng-B;$ zho{X}WvVhvIY(IB7$|ejr47*fN{$rjty>+t^YD{q@KwxEmzw@lF;Y=laZ-_5u~O0c zGV)UTa`KY;^1DmXr>lES36GJEf$ma=Yk12t13OQvk**caUbB};z_c51v~%a<&1dC5 z@+x|pL=HIYjqz!qPEYWH>11uVk~Q)WZ{BYr+qj9g9AqtopFI-FUDn&rq#ijSsnS9a zKO|URU}YSe{s|uB%>R-gsAni^sBLJmCn@L#MiL|zWCF~M($3lk(mS%@34}S|T`uxu z!h|zmrZVWEid}%ygAQ8Qx_Mc1Qc=5}vz+{to$295FbV4gx7X{dqEPO=HH^`TD=~r> zYvWvAMh=3QHxQV6SXX=0rU3<;9$(B=6V477is+QQi~Zvc03s{~4Fziqgd1!ZHo&sg zNvZ1698L1feWbV&C_4fxZ~B{;=zdv0*2_%jw$gWvHH-Eet(Kmd4SCg3>YOmiPR56& zEF0_OTGa{m4shJO;+=okeVQ z#)dcSQOBjTM+3t0eke?GpR;Bz4xRQcHvSBa#MK@ zfZe!JfFK-`Jup#bdd=pVj&a+13_|tdVza7q0jEx&y1RXedfsnUw-8WlOjpTvGt<6l z+^DK$FaWqe=kjq!B;J_BS%;YO&6wMkhE=JwIb6LCmpl6&fn<_dt~>n!x92e@gf_O9 zokP}Q8=cfsRybE<@Ngpu$G~yvQqMeW=~KQOVoRWLNZ1_(U(Jt?^b;m{>fo9WBS(cn z$OF!c5C2ZADmlQ0jy_ zte5?scfh~}mG|7a0e6a2+LB+rZij&TEOIa!Mszo^qHD~`a~r{D|H+^OiUH?4c7=*= z;$-;kF@i6=hnI1==JhZ|A$p$f662KpaAcdg{-fY23N!Y6ap3A|!_Vg$`H*6(BVoA@ zNDMkv@)e}Yq}0I7!_Y4>Aw7mhZT4FSy7P?+4GQLKr;`xT^DsXQIrt0?YuPo2MGQ2I zTuXw4BmA2+r!tahLg8?pn+jKNW-`&hGG-?<2k3!VXw~y$mbCXl7d%OLOLB7h!3oa7 zs4p|NQEn(}LCnZ_>AEXNmPNLOwoq-+YD3Bhdg-qz+%t$){m$k3`S$S#rjA+Hv=^_t z_LtQou|bCWUq@{x&*OyC21yFj7UNbf++ik9sC6BoL(?BqCQ#Arkym%>v_K_4tIZ${ zl1Y5e!2Eb>zHwD>+|LFIK?LeWr9t;c0RBWMmr|+2LKzN+L;nUQ!IQ_6kB_{s4g8p@ zcSp~dvQHY1mcv|Y42{ErGe*)19Za)k%GygVRo%~92hIv)U-K(t;T*uFkEKl7yt9Yflo#7tLc`t zX!N=9VYGBJ~7cQvzX3%=Q)_03kbVB!j z%>JAKCmNfBKKebLa73GHf|VU)3ByKxA7fWvLp*Hfpdfo~inPr2k{Q};k3C#w!wBKb zle`9GS%IW}6NU3qN&gm%;7?x1x8y&@#r4|iHKlm$t7wkYzvFOJVP}3}>EVngY=y-$ zuO>q3Rg5W`%dXveNe==Y3cfK}g;wHU26wAgikrusMc}Pv0lA> z;s7?*sA_3>uBW|wEOw?YUQN1fGy&;)9!hr0SSBDw2ONP7%7A_;O%zTcs|w3(5)e7j z`ao5wN(v=%WoqWg%2sPq!4Dyz0N?~+$xtXO<-&!+h2V0K$i+%E9q#7jEBjeLfc#4W zp`MSq1%`Xq+ey#cmCu)AV9TK~%eyGR>g0l`>}abx4RyVquR--3A&)7SePk&u=i>rM zz!20sU)FlGn^jZmBZnI`3z_hH(?;)fGyVKs?ADv}^>}yuTfKEF^(1P8A~YPnOkLq} z5p+md6@!a*0_TU_lAQj;_eZ~6g#3iOA2zA@bj^4DiWm03f(*&9;B#TC`y%@>&fQh* z;l8)DfskzqZBz)x1_=%-m_3a0qZ6wM_eph+G}D33Nliiq7ju;0gh)lcOHJt8gqG8~q(Teby6a1Mt=UlyCY-NWrdO!?58T3tW899Sh}OOEu*D@#bp&#(qKEK?3dl z_4c8^cC^k!5saDHYxUyBQN+oU*1sYNp3|2jvGFVMx93W+Jm1*|q>97GZy?vPAF8XQvY@OWm@hSZ0i zQDk{-nNmzvi+$1VaO%231L2yukEcO(KIX;a&ZO-WPE*U)a77&M*3R0+s+>0Fmr2=t^DpXh|LCiA{$!3%i_ZE;G zftHQ9Wgyi$kd}Xa*6{R{2(=8|DbSY5UA57egFJqyUge6f6fG+?I97o`A#mftyZAu) zXn_*LzHRBre5Fd<>p;w<`t_G*mm|!oL%>F(|6xDd zkP0*T4lUv%?FOd-37a@i7Ewk(++!h0TXV6?u|u&a1_5RQUBWfe9ac&}N3S1yL2>4? zEHYfBpGz76T#5XnsfDk)q+=9}SX>*%9lV1#{0W=Eu>UWm>xiULz4R z_I{YDGk%ochLAZ)wF~Y*b&gx(6^cc~=J6%$g>HIOT?K9Ech3TMmv-)+4)^ey=AC>^|)6HRVV(=+m$G8@}A%NsmegsBLGvlF5trCv=qz_GK!}w`R9&%2(qG z0{&*mz|kp{_->#+2<@If#)vb1jOvKc2BD6vkW(#{V8zK^i5@{=M_ zt-^8+W%@zPB(|Q@Za06&{@6c&19&*Akoo$P;{|b4hSY=Z46ED2={yVaV=H)xw2(TL z_Udg$r6}%R!akRKNk>zHH5t!$)@vfZjLA8w3mBxY9 z%1Yq$iPngLl`Y$0IbLNrajb(9o{`8K6wPF0rH`evradG=brese^d$wJVz()(GGxn( zk>$=j<%1UiVVoTP<Yi+&7$k214W|sO7SZ3beaSw+eK$ zZ2XF{0AMFTZ8+)fIPR}MrD}-%-yZVL!FCR#>%j# z1JvQVOE!akzc7Z^^)-y<RB| z5JnWT@T*8rLzCKd^lugk!$T4h3`hh48j=<@?<6^?HoU9hZAC|Hb)yY}B&wwRo+GiF zo!1+H9JtyL;CJ)0sm-s&&TTuS-pEw`TxkZy;twkU?QLe{?6qNb4OQoCMf++W1A~QEA4XdfO2jv{fnW<-eQN&n5gS>WDbGes~Le z6roJ=l{@(4HR@C6i?4i?6R}3aFsl|p@)WO~6$okQuh&;PxD^NljGAO|#%{4c7zpd2 z*Z%Ju+kfb(e;KSS%$)51O{!=8-`TH!dGMlUE+)=aj;;>QU&6>PCT3r>*DuwzFYv38 zqokRYg{3PI=NIUm@o%@kK6(~rmM;*it&znS$(7}c9{+U-VPOYPA{}}T4t64X)-QTI z3)9yb%uGxiUp)H1et9E1Ga^Q9aVZHc8CoG{D{Ob8PzQ1i>1?KPZ z|I79@zW-0#|K&=5$M)~L{2lLq3j+QpQTbo+=6^I07**ViUA-K?W~HjLo7uk{!bUD; z|8oAXGxcxl|9nhf{+TUbf(`%7#jmH}zXTn`t(;w4MJ$b+zn+V)d294v6&B{NnJZ>* z;_%Nj{+h!k$_{o$_Wy4A&w2V+cj4dYTvIHWMk`&NElJZ0m!jYk6*lZ~JFxKSANX_UZ{C$S(OV33p z6(Ha`6sUcn5-}Urb~Yu}^c$}kIt}G1w7btr&C=#q0eR5;8-K2z0ncPUbsySq4Nn_3 zbuBX)Eu7UlLbAnKP@-B-+vGW-e!@1u{(y^-mIB+adLPeaavxEr*2uWR~g3i)q%q}gKO z=DlOT%stBV06>rwP)BKh>UM|C!YU-8&fm@8szoxzaA*z^0qockL}UB^x6>!}B4`ge;le=xXcuF$z~AYs^9ILm4D z38MV%ih<;B>deX57`ak%wW^ac*L7$7G#8{%Qsj}(B^WmcexTB%4m}R}&e9mEa;iqs zX*DZ6w88NH`6?W*5_t!B3aVS(j|DC(6|jvBfEJ?%B@{D7!I3h# z0`M(2p+)XfRnJ9JLg_&vWWGl~Z8DwPrH^sIt%s*qdVy2#_-_UdPE!WxOjxVj_N39o zS)hV8wiTkQt_hac3%BQrN-C~FaUnQj8v;32)MuwJd(nWZcNF3tE*E)q?~x3b+aqli zgVZbCB;mEJcP5j6o9R?lk^)!f)w`*bs}9Z>npdT?R6*~bLl#ns=yelxIJ)1=Ow>ak z4gCg2?^LM9sqGcKViFAOW2>ma5+LYMduGQ3ecv$kLL`IaQk>5chra)c#$#O%Eho(2 zRCJ|HET|heJ-~&2g^^1(J@#oyl!=SSPqI~J{1NXFxnE|kW%Dy|U72|+AG5TOgL*3j zzFh*p%Y%rgBA89|JBe~^pWvaKq@A!dhIjnW{@O>wGVgj7zI0W@_Bbeyf)j>jYE*lQ zoB{8JsCrPV8TDE;s0qXn+7ts-X0qR4j8LS6z|_>zAE4f_c^RB$9_7or0O_C4gIp}I)Hs1!ueG9AI=U$4uu z6GB284Vn(KB`@HgAbYki*pIm1r?Gt9;QlbcvJPH&J}K6(a^qQznGxnU4PWRrf>Mvf z)?-R>^LU2C-)&nsG!ZK!dV;(h2Z{uy1cs~Z8-PQu?*Rdw{L!4aIEgZ2EHySWi7@Qf z97tkKtq0KJz^uVucqK+L{DDeG{LPC?;T1Rw>usVvf*Ep{b=p zaBj!YP+@_?8zbg*3>qOK=Tdh0M&4URLZqK~WzrknAx1ernFbJg7ksyHrzc!aq>57s zeF?)pc0m9Yj}Ld6W+$R6A`G_`l27!|2v9bbYcxW$IrwE7fDH2vwMKphi^Gmv!J7sVB}cO1lxG!VT&du^dm(Jan|Yp=;&W(wQ?b zaY`m0H8_6fzlUf9Lsy_8P%L~JaY!bfzmV3^PhsrnRE;lTf8(u9G0nURr72+_`l$)H8L!&H#s+)NWeO0~vcRXnU? z8gVErV#CltWMyXIJ}6*Pu;}fb7mNodVAW-RAB3e{+kGal8bTJoVFH^SdtZ6O}ygi}5qp$-ov zm_@qL9s;NbR{J2$6_aCl&^Nd(BXkce`Y)^cfFdXXOdJJE3Y`=>nv5y2PqGqA%8)E5 zL6!M6a2_(6-|(RV%dhVYRGg6MRrPNbxJY663|yR%$vbGi%J>1K2otNj?@W#Ps|p*f zkAC>+WwTb(EEf()bz3w>wO57|@aIh3Qpruek9sZGv&2V5iYlD5M!uuX3b^ z{%J|0X=3Ln0yZVC8C{R_^-1X;dWzI+{gb_TM`Wy?hFqc9+QP=e^4TVDQTnZgRNH)c zYUN_?VwJZn;mJmvjczG^aHwsohTG(D3bM$O?dWs`$cMz zDOnHr_XXJp-|PXBJ$xzjR8Y857zae2O@-J{=ScjJt@^sHrQM# z@HnD2O7v*Z<}w&cAl+OVU-mV_QdmjFqw@RFJziU`*;zSo%)<3-)p{(VRGR3!6>_X{ zpqAmeas?2Q8sk$7iOAw$l{p?p^y=dpx8isQ-sN~@H*b^$UPd@W;XH(Sjs%*Kz4JZ4 zeg&r8-k^HJxHgdxl0F;o5Z3}@C<1yvkQm(H3J6e#xVM4;I0|lFtSvQ1jT`0%F<1m} zLTGuY9TlE(ky}KjYrqn*a{qLJP^*Fy7$D@hzompS#3qx*heSS2gha4P> z4_#N0Rnks6;widS!>7YN(>;=-aC$Lz5u#Z>N00Is9Z%9s{^B_skJ6~C=qcV*G!?Tw zB^P^TM?SKbxauOyp_S+l=pr6*^Tgt>U)AXEOXpLMZdeu?*^+rV8CfLZT`99c$NlYN7h#3p0maU@A^}{?8du;L>q%CS%O4m?pC}C zXg9*rQ!=t}!qdNph;GKOjaSRLcy+9;rwcYES+hDeflo^kLt2(*i3@jOAC0r^!1*NG zVS@8)eCHPY7Lp;l?5-4|Wk;#JQ~N1eUeO>B3S|ca#a+b?T#Nbd z95W=%FPZV50Tt+r#;%43oiU6YH9y$gL^IWiT*=zD+ZVDf^l@Mw9dt)&g{6o$8eG?< z5z?w9Yon?yW$mNc!d@vWJwH`HNa}L2tS%@SIi?+c1t@9Fr)QSS+D6h@$-*tIwy^u! zR0Dc&s%oxwqt!;g*&LcIS9yN2se{(lo|>@VG~oQ!e3b6yQD3B)d9Jbh;9;R>#($gG z2V%wP8^q_{A7X3U;!a_Gm`kqrFqtt*nByD5k|C2zu##d0siU!@vfESRvgn!!HvXlV z;&GVlIAhdlQ7tXrK+h7Ry9PCItCOhaqMf+!6}7=tMTQ=~w|lvc^QU}Sf1$Y&`E#sJ zZE%C#hVfyHc%T0BQS0qG*X`FzjaHPUo-I5v{dl4m{rFGr^a%B+0v;Sxsw~F>S)`># z|Lz_uCf4W)EDE93yM52Egrxj@V-OC%=G*_NfldxJdT74L&gznIF$>(gx~54v|}$;Lx*veEn+tK z%G#ujhIpY1@3k)^+#lR{NAqljhrx5y$dU`yILa_?&wfF2w?`r*>9Su?bQq&?V2o#X zAuUv~P{>d4|x?O?t~?)dQ@e=)@?!)ftU_(%T^$vxG40{xI=X?(|$TxxvV&9leUA-}W03-^criA(fPvaERnd@CaOj8oQADd-NJ zixl}Y#Kg-ZsAbbsGj$=62~s|~R4>|lf`~-07iKRAm&kdeTND(E)x>XOgG-cMp-Z4g z)QP1OD-5z}xr36dW5UjSS_Ird{=sCl=7;+gE z86IQDra#_cejNu49vmOopWnIqok(P08u|&yz&}v==~u>IVRy6No1R$QXo!!CDjq3E zHzjY0=Tg?if0BI;bC_`?$c*tuzsum{6XSpn?;3xt+qXm8AwsR+aG!F=Ul&Rw6GKQM zJ;JiQr(_Z4&zH{>u~;Z^{r16VZaleeGhi|Ryf2p!DFEC}S`5UdEt5+TKISl9b;M3n z@aeoET~|)FJ=5Z5(jKm>n?piPrf;KSF6` zuI9o2RpnLY)#sI29XbXyz{?iKW~jVs9X@~2+?^0WSkybhhwrG|HTh@c4}+v|t|GlF zv;vHjks=qvyIS&aIhRm2)eHt^1onmOch?3sg$>fnB>ngq%b#qM>^Yf;EmmI??%H-wjs>KN4o%^FRuN561y{r&;e6AgO7@T%PGgR>GzLz9bTe8`meDM zdheIXxk+IqyE3OG_c*1uHqPv3nKo=bYytDu&C46=mwFqVEuS^NiaVOYSP_Gsd%e>U zEgbPSC$kSIym5wsHSyJ}dBgnD!&Po~i}fW<4CuNDTwIv*&i2l5#AJn$b!Fn>X!HSu-E2r!{OfZ7fq z_nzU8{JrV*M4sjsU0J%nugu9E$(>-F%x_BRnNy-xo1JQ(>|N5m0DOMBP?Z^x|8#l8?D$~NHw{s3 zi2cI%lI|<+wuVcHt#FWjO*YjXqdST+z|4~KDn&VCi+uFExgn)mh%rm3d(lClk! z%PZ>(ZF!E;DLa94K>~xQL0-b-m~*tEoVJyK`yNuocT+nwe^~ERl?>S;uLtTqx4ls9 zz4zLbowDd?1%u3tA$90c$5B(oo$yWCa9vRV6D;itv$U=&M}_q4Fy2mnaC&iFRg-6E z&6q9nN?9FgFR?&>F4aNIhPZHY*o*YXQq<2<+1eac74*95bWMGZ0$1Pa;`VYT?Rdg* z*zn*(kPdGDywyW@ZhZGF+m>`W@euv8Wy!)mD(qb}t1Bq-B77-koQ|VJ<=i+rq5eb% zG*>baD(gS)S*AY@>4lctv}r88mIYc=8|hQ(w3WG4C4NeehG&u#996|6%n|9aL9U42 zLB=GAQh901%_=W;TW>_Wa>0cvCV7~dP)U6gmK#68)p%laV@r)K`w$;|7<)huTevOW zH56U=p}K&2eo>ZNp(@9uuz2Jz?rKiQd>+TfLGFi~U-$(OwHm)O{gvHXI+E?krNqN8 zRmK=aCx=q|F|tHq+Q}oU%*t1x77kq#?4CuoDg1F5?sG+uJiZW{ZsleAl1RIeIDqmn zpWHmmWN)(IaxZngs^yBn=QYV`;}V)nab1bzu(jAko{X2f_yoU|&2&wiAV}C}20o3? zH;wgbf1wJ?S&?`VE2@*Y6GUu8hgM!SQOjsEj^?V1D*cmmk%ZQ#_PlT9w@;$!mE3;v zZXKMN_>W%dCDvYg6#11?Q4|>*m$#f^iz_)$e)6L?ce(T5{_Cmaw!OnVsZ1HY8JfU_ z7QdE~H?qsLPyEHSCQh+3fc=oo5=K>g$)5+T7ey+n;Ga(br#7y%I1q zISWd2W0`C1GCpf`TXRv-Tt&7uY0lGkD2|wtZ_y@^m$aFV(W82-gJrX*sWJAxQnO+f zi36#ZL2#_)C-Bou4?~IZ2)?m=MzWSK`+-_9xrCWi?AbIosF6{Qj+|uVxH(U>-R)-P zR7hB6AVv&opD_MBXZEeC#cr%FE8|A7zB)hCVo5*O$AKPC$R%#gy1Y8K4Mv-=KQlbJ zY0Wmh--f0$BhP_rH8a%?2(`~ytedN0G6>%@hOx`pO*KxIyG~D)OcsU?Z>zw-5w|l} z)fK@YxezTEZO^_Wuajgd>iWCGun!<+Q_F`xDU82#Sx73_ z=kP}ezE<>N9WmZK8l+&4D2ofC?reT2#DR0B$=VJf(0&(?*CYhn=!N8IMbzf)hfj;b zjo(W>R#1$ol45VpdA8#JfaB-`47e+i&Lpo$=h^TyYuE_PZI`|+2*S7F=yS@m*Rtzb zd(ZO@^O{Cz1yCATEw1Z}6Sj95>fH@PJH{z|>$B7@!Z&@443RX}CW5F*FC8jdl(eG2 zD^V&8@82WeN%yU_Z@yA~vEBZn*eJ&}`x49#-PP0E(h7nhl5!KGgw4{Tyd;JW87cC! zW^xagN;();OGyLG^kK_7vWz$=MvH!qyeHmp9-v+WW^hGTo8%2>^0M4;x>r~25-5i3 z69SkUbY*BO1nc}rySFL*?!fs+3{8zrk}VU$8UV!3c?*sG^~?2HSIbw+*R8x&pDFrn zbE%G$@30{no9dk!uXPANt+i@OwGHfY`^_ygAEr;XUw7*@*m~za2i`f}l1}74vdptC z;Jq?1ufR(Bu3~*W&ZmTIGE%hAkE#Jb(72d!pE3kyy#~%;dv6C;?=KHf`SE`p7#VZh zQC?m)4;;3RYs-RKr zFb`sTzX^Ica2+jxETFnBRKWATDGGbpl{@qz@Sb#=n^?Pt-`mYu>!Y&Fr-ViSXu0t5 zG@@p5ipM~tOD$&`tMyecOQT{wH!^7pGEZ)yHIMiq#7ssZru2Y~u(&;^zWlv;b@~8M zh4RPt?R+ERxR?*G)v8%-`Gq}2n9c9!{H}q5){7ML!e;V#z_Gl4C5;hfI#UQxjRpqc z1Hm2Q@w>L~y3dx+md93IRxT@_hqqzDHS;a-wZ4-^t^6Tk0*PGACyG<2N2f!lLE827 z^YI(FhpH}>{Ud+wz30_oPf5!{81NXT^CMhMf3 zZ?OB@tsoo+e~2%Tu}!~nE3F%yH3G;iX;_(2cfTz%gcl3!yNd`F4-q!!#;>yOw!`9j zrn=~A=xx(sXK~^D;(;ECZSh!Mzi3%f$F@0*B}Ey7E3<9!MWC*NL5tF*{duGWqSfO1B!$(2fPkryV`zQ;w(pkNW zQsvjlw~9XsACF(T1kS_kyOSR`tAwAn>wODhO*KkAGHHj^s#1V@nMUR6rchb-CyFhzehO|`W@>hO}<#Wsa zLfEok=96#&urN^o=c!VmR!L?^crys|Jyp(&0i7R4!)LCeA~Viq2WZ(Lu6YN|>%q<$ zQxiwtUafRvT@_QK|AVu)fNo=1(5(~4jvX^I#>~vj%pfz%%*>8CW@ctPW@bBPW@ct) zwtwfIxpU9Vo0&KNTavn@R&}*h)oV#=?fsRVY8mPzMWl%=_Pt!MAop7qFe(tG#g&`e zElm#0EADl{;$#KB%_Yp#_G&Dn=dJD<)L2_>l@qyVrff&v9c^G77!20szhu_W0=hgo zdqyS~nJ3|3nf--9N<&MHUsJ}*?wYnxdpf*k$*SFSP#TBSeX)4*o%?MlU$m(4Xzj+D z!G8@|we@Ljjas7|riOO<{XtaiX=EAv+;tNttroy{*YR6o@5<3PhEOU8bChT@Yn_&^ zAb2WDiH6>^<80Pbv!LFdU)>oNG4JiB#}3?(JgM^z@-|<#a}&6C@#1}*G<*zB(3IG9 zz5nSwl9V*=>6X|bA1FAe(`%no1lMwtl5@EU(!hO?a;m>u{z|iF$4=-9<`hJ;L-kuY zO>v?4Fwj0bzsRN|IOSn-5^H;*wI$x>5fmtJ>>u$h9jTf0M&=+j3rd^qM10U-0@E@Qy6+&Lq98zO+q}o26XJnL_S(tG{oY-?_st zGY%9}-U4#{n}wj{6vTaXd`lV(K^pjgtZ1V?_V&v8p7Ez3KFDo4l^Lb$YkBl=i24gv z5=_^q(U-+`>QEKh_L_#*K`i3}Ua?fuNhn_G$1xx^WBtl8(f?(v*&q?KVgo6Gz-y7Hgdj$7u+nSH}t*2kqx{;7W69+|c*p z0h=o0fW|JOco#G7#73mEc&p+$#q=8DeS|~)T;y8YEshf~6}5fJm5Tf&`k^iTHVs>c zAX73;TqC_yoZ7fy;Ma_G7Wk)z+Lw?h&9Bb#GSHFBmitnJ@88uFq=5Kfg}tv-7Z&|x zemMEefYD~lL&aR{i5%w+9s9wp=H}x)X>Kx{vGDdNp4%3`hFfZL@^2xuX%*EAXm57s zYl48qS!>kuuNAnzCBUlPev0V^g}7*L#YQ3*$rU@hZ3!vlsT6BqM(1VC5I9O06IE$! zu}1SbNV%jUKCp@Yq=W9(LK*IlMjE``>cLSnI&9XHGA#;f=UZ;E`h=zAJ9=IqJhe%i zLXcnIjkMd%B{fol9xkAS#tfOm2WAtVGoFsk>Q5pYe0-^IuM15fyDP&-{ykjV*hS|( z_%rEmU7Jcx`8jYH1Qu=R!(nJZP5c(*D|6tjh2KM*SH))5Tjl1*Ol^MdK?^18dM@Jc zA)ApW{tpRnmfhrc*$+N%27az>UVfV_P7N-pZ&re|)AbmDO7otQTS^$#>Qo_1LURPH zTk61Jhs48WX!VNh5k%wI-bSyG?6n+}B=MPp+plQXp}1{MCn69G@^lBpu__wbRY?>! zaX_C#b;rHNXwf4IY7u+QLRH9Z5!~zpTj9NLxiUYDgp=n!IU9N6j3tepVYJvTnJ$_> zj!w7Y8`CigllPJl$QRB7AaWGpg@*8K^22DoSH zT1I*6=v!!nXTn?*6&tG(08aEkg#i<+wEFtf8wf5*bvGtYj#sOv*BfuKpmFyAxDFFk zR_~^0@U$(Vy*LNR=0cQ|8yR4y;FMr@8T~+D6hGa?H4X(F znNI|taTezU+JrxdF$fht+8u&fzyn$Tq(H9EN;qvqdVpG?a?>Df8GUhNz7zjZFTsY{ zFG|4)3~@H0qG?!9`Vs`90>I17QOnm2p%%h6dQU3uxHU1-x4Mq^SLaMmC;v9~D;r}9 z7Ish1r^%bhsglDshh5&OWU0NEeH}Z~yH7W@n5F~lJ18^>2eVTeh&9WQ zBdAb=(NH!SWkNx_udy$7;9qtPj$=WpnT4A6dMkfferZ}q)mvQtgf&VF;%~$sDGEiMgV)Fj+7OcK!wbF8vEm`5Z`;e+Ah#C6<>Px=)D)(<7c& z!Bep_6w>zQ3Uf{d7cLvLKpzjFI8yp)hLY3kQY(jCXuHx=RIaLHC<7|ZDFwYx;g_4Q zn$_`S+;7_Yh2bj0PQnyu>rxfezMoyg9m-w8?b{D%`nBmB*DJl^*!-%{^os4UfYs>K z}X%*^s8_BzK&ICVQHfhPNuf;^oDfD6=d0afYGxt)qngAG5h&HzoEi7zHU7Z@`s zk;3bk+0pjto&`n)x&`h9tiRtgmK8{(>1S*e4VfZpj8I1uq&?R~9>9=5mD6wSSELb_ z#hwRro75uuk(ORNN91nr=#vQPDIimlEQvSMXVA)J<~svcOQ&-kz1x(|Fd5adJ1?KZ z9qOj<4(Yy>PiX=hhdkq*mCI)`@K2Az!KXS?rT8iZ61xvS2lUB2pdKx!xe=t};xxDF zRP`IL%y=!I0xB96Ik?>R-foI-c&Ch~&NVY!(&D+eTkQMW1!dnL&{)-e)Xi@bJq)}n z?A59vIgydfDc>odX|%FF;_AaCYsctZy~%qGKS%Eou8KxPYG#^a$7087GM+N7j;{He z#*1?EsR#hOxK9(jg63H^^K2o zYYZ5Cj4h*L$V|kLOVLBx7l(QtX~a1cJCh3965V+tBS?D}7AxFu*p?qb5e-1`;-%Nq z8kzpFjEGKX&Sft|Y)&^8)Gs5D@DoZQc5@0%znVq1r24^#Rprx?-co^!C&~>L&Gokw zXpVhj?N9Pg)-#?@0dNbXDA#+|>+UC)!+^xs7qC_&L|bFzlRX7n9+Q4YF3MJ_Hi=^H zY;knNJE{gQ=(z{YigWCF?g(krt%AYwM(ctJlaxcV8}6KWdmeFdkum0w-l96yDwMCX zRlV2&wo0X{Wj)@s`{|WqXeK3!aHKkV^G#I!MOXAC=dn)M6c-_`@ka4_@!Z4IDNYY3 zC4HzrIe+_wo~M?6rxEvFoai`Myh%JnQ%Ehgp3$CjxsEAo&}s7+-cBUF5Ks(4V`zU) z_Gq_1cP{sSDmKns@wsKuZ41~yAsx);#r@7OCSFr@hwz%a2foWaWz~)Fig2B!Q)oXW zxzL+p)O}*>J`%(f^g@f}1#}y>>AGpUW2|AHPdQe=C(sP(MQv!#uy5N#Z1m;T#Q{_Z zgn~!C^3khN->FCj1*o+;g?tV2V7FB_1%%)7LG*FhT!qxosiXv~X`)c5{8U`Y06k(Q z3_W#YOU9_pjUZE#zOu-#%aS1n8`vxWDN!t7(~T^{!-9=kF} z{TjS=Lm@ooR>3e90`Yxe7|1N7NhQzw3zN!83X=N7!_Fj6pmv+iG5b($K<2Fo0u+VL z(4{TNZ|BDvo$@G)|M4uG03vA<=&luLucT5o_lYIU! zd_i@&v90Vb6O?`1@BYe)vz_BVRJQeW4JVVNa8}-Et?-f^X*^ED^^ymUtkx+!o-m16+Sorj1Q_jN`SWws^N97uU>6bf z8vs+mKx-@Yu<^*7$O&0RNob1e1dDIbcF>TBv>K=2FRxbnWfwEY`u$i#F`T- zlx-_Ix}Bdc#6`<3b6(@2eMhr9O-pbd=cOFkFn;Hey2D8Gw92l)`|jpKaby|jCcOU3 zc=ajIOVwNRy;vLai1Q}(#*-`xG#vFz{xhi-!AdiN80Q&GIOutpk0}A55?2^+EHacw z%kyrazrETu3r6poufI9+`w@uMPc6FJ^Wl`IEdq*yac_3TFt2ZT0!GGeu&eoF6TxKrY*xUCy;8+mlYBx>*DjEGt(y$tc^lN!a0LxwIOulIIZK5irN3`~%j z%ctxk<-POa^EU=5D~C&Gsg6^jHB(?Z{zpsPtK}3nNers>mbVH9YD!EsN zj`649Gl1KP<{08w5)Q&S)$k9PVU4!qW%P`)Y>qKU34WbwqXz#8W@udTR60Lt8fT99 zI3ecYX|wUhbZ%gU=r)^gh$F-nT+ypk39?)3g%ww)z(aK-?k&Z~4Ae})d z#Ah|4%I238P53f0+E#NTE3T=bbTf&yK@T$GP{{rvXGlw^Gcs@fV%aK?qln!^0u_%I zm%Kh|t9naCi3{+eT~(lZx!-pjeXc@DY^gf)^yLLT4_dwal6+gJ%%TQmZ)Z;gN{`Z{ z)NaN~5mni7;lS*{`Kwxu{_VcdmnJfebV*J-Q<@=I}#EN)@#jF8i7GrQm5WwLjj8ij4d5Zh*rhEike!X1mnEMjAw0<;(O0O)pS4$m3``<3IAe3ai?L6!lC__#H?jw&%HGUJf_&CxSDWK@F3=qtHwnr?E*P{h zK1nQdC{vNsz!Lh6({vf@G;$2K+IGA&X~xZ+0tY!a5-+L!>h9bYg9+{aGFP=2V`*7% zh!}55v3VRjv#BR>F^F!+BSv-{@k*dyz`5%g-d**&@{C|Q>kaNib}IAMo9^R8org(+ z^RZXASBb0qeT=wzVT+`Fi_9l$AEcA{t#g8b>NBrpoCNzq*#d zsD(Nt+95dxC@RF3vn#vqksA!o3Cy!9_e7x@1$Ya}$4)UkUwBQqUWZ6|rM%|^1>v;4 zU%$X4n6z_OB^e)h8<*wNe=5~L=U9SnMaMpO@_(*jH1h(H6G3do?BZAlT z2dLYZ5oS!2jnIwTFzvk%ijF(pLAfokpYqyWFt_ccJ~b_y^vt89-PmhX_3DQYRsV9q zcvln5bi}05d%%4Y7Twnt8M>t@&-E5O-)a{{+jl5ew?pVxuBZUCg_}`$JB`8v53dC+ z@5lO3Th5!z4)-_{wc1cLmYca=@eU_%nr@ayp+tBQ1HazAn&{q-;ph%fUp&kdnXl&} za}&w81&qa($IOWfB)5%LPfeKyxi6SOrYy7-gayy08g1|P;%ihOo7Zhv`lufhrk*^X_Mhfxk;t(eN}@Rba=WNO#I4qD zQ?Xs@T&;V!3W&{t!h3Yv2q;r%t3-*deAGRsBCbX&?Gt%s6P4|#uOhTv>_wAK)*unV zYp$<4I!D61D|kmTrvm34RyFwM&W|y-z=WWc-M<&iW#PT>D7X2A*sZK)aqvX+g5_%W z0z7(rbpMK}cMoh2We+?IZqa!kCb<>{fw+U-&2jDO0QDJlKi#!q_t-o~I>l&P>e|BG z-bFowg|>rE-QX7_U{_iA*Hs;+pW$>++n#`;wzjrVYGz@Wj=|DK>!72-Wphkf`LZOg zBe$F#{yAUIiuC#I%Bg4G4k_{AC8q2`nLnXzLrmzpYU821oF0Ya=A5rPl&a4bOb^t{ z3LAhcs<~|2^wC9~B_je4ji$d%UA+4RcPmCMz$B(E@n(kz+>&83b@^3?5GI-ti=_5EGH)z*TBq5i43 zV=+4b*5CW=qy>Jd#V$}xdyn;tFtH7@p&qVtOvLd1W3-R1m*(hmBPB=&QMN>P^26-8{5{f!=CIEh4?%T#0;n8uH3Wd{z zm0RG}H7^EUq0{fitt8lwxW`{3zVmifE?>`sD?uB4%U9!Z!Cu{vO$KrEf03MRW}rqY z$b~_?s&n(t78AIkKA?zMmt@vIF?zAoKhT`D%Vow^!Y;Rr*^z>e+v6-Nauy_&j0tIK zX}K)boS$DYR*<>nwymlAlW`?I?nE`s?J&(=84f8yL3*xOmntjUdfcV|;2LSi*gD_d zSpPK5KfmLgW#(%)9@tNPMcdtoswDGzcbw?iQ1JUOCE9E6D*W8zVk}koaNHtvStT*Z zF8Ah9Xn9d~9XB6X)dUq{H!MI?SOz@m5_OAMtFhdgG39qQOryO2l}zpYc5zE2BTsit zuuSs`;oMasVf)kHq({yo-6GjkT(yn|?Gn-o?Q#J}=BSF1^cK>sjNkmq&Q2q^a6I+~ z0_Rir1r-0TOecDuFdd)f3ke!BEB)tObgO&|FGx0e-87ZGM^zUOM7kNS|2E5DKq?yXhPK|u{wtB?r zO7Z>l;ce*x+|G86RejyTUmN7(sTGvflywbQ;jp0Wp^V)Uz1}%l=P_>^)+Ao_9f-KS z?WtJz(_@C%rIGNmD#BYYm7ZuOz;yr{lpAF9Q`o)4D38Ni>bRyZdD{rs`=nn#DF<-v zWMX?O`XeS9K-l*Z<}O!84djgO89G7l_-@CcC4t!FuQ?K|TvzkWa-w6OsCEdhnpD&vn215JT zl}?(rw`eAWN>Pw2XvfwWpr-_c$7irLC%KFZ>x2f&B8+jyct{I_C1TWadEf)5jg*Kb zsJ=QEdfX9Fz1~%yxPP4>ey&HmshTz0?DpSc5y`wu$b)9dfM$`|G`0<5tkh&Qegg9# z7+)VOsc(yWl&eazyoFud6g0e_+iqJ({%pu5MUa{A*oENqV?;d^evp!Cr&%BzjOBix zDW=`rt@CD>rEaWghPBM>FrEi~t?TWW3(HK$O+Mpve2@8?Z_p)#P8VtnuaV+Aa!?f7 zOojggnt-+kf)XF(mwWX7DTld#_`mSd{>UtU1Q&V+ zc1ET@&;NPA_(y*Fr}pwMN9SL{&425q(X;&TytEG+jhc~_<-<$+kecY}=;{B1m!>SB z^6`V|FE5Q+z}mv_Kgnr-oV0()X&*21Z*m$d`yVUk|1PJ|f7JUU5dA|=gJk^U_$c82 zO-}nGr2TpK|COBfPo?Rv^6z@rKZU1%mVaveo5uDxP3kZC?62j&+W70mG04_)rB`u}{zf69N? zwEo!%TwEVU)<3+gzf7$EQPcV_2kt++T7M3%zqwj}cJ%S^zj3wb89t1|f0|i;e8d0V z)%v?w_HV8h)1S`!X#C$@EqXS_k0bPB;2nT=AEkw!15gO~A%X$_jK95u9Z=5_(&cln zo4M|Xt9796p4s2JKSD&rD6Tt3Cr%=Th=_RAB_@*r+C>QJ4;m~EAEyM52oF}gffR;1 zwG>HF7qAxcT`8F;2t zKc?2*n%bJ@G>FY_^^Gj8mPP||bR7xQ?8I`dDoy6;F|r@XZB~c3p`#A^oIe^bpM)tZ z7^jPaKdl*~(HtaM7i@xTZUw}B%%4HEc@q$txUOIDx(|^^qi%O9=Q@-ahg;RSosj*u z$t_XB>O`fjI~Az{u}`H{ZPAoG#x0?$vDBvf6F2oUb`9%ql}W03L9|bp2T)`|i})a4 zz^*nZK5bH&fAN`eW1k&xD|8S@ZF!*#UC#mdugk3Kly2|mD_X9v*>C2&B0(W~mws*2 z=DXdwtSXkXqf*1`IEvd#QXM11k5XT-zsHi1eOGejm|07+=|paeJ)Yh zSViLfy>>iL0W~6~)`Z~doz+^s5iY^&o=LG=bI(U;LS#b9CSl{df&ICRD4-2g<6F*5 zoZOTKdrX_ut=poWr*wFuQO|X_OFIG}ypcnw1}4hV4@GGo!R!!0_XL4D_tf7(!00*u zm0_9qG-eUiG<2jueA;j?VMWmOY)D0082Z=uuga5DC1Lrc-y;D0Se-P!UqEd^zFn!( zmlEi16wvWy?JL1cee`2xfG)2On3Md^udBo6Y2Q9RUHamTe{LWD_xg~-ydHJo(O#hQ{Te&)`2P+?GFh&3b^*kot&@?PzopSr>#8u;wV zakcXXlXq1_A!>_luj$`r?x$bX=BJ@A$OjPNLD+s+evJ?)^l=tIXC(6m-&FH|7xmtm zSF)7!&igS@@kH2zuQ`#)oHp_DMioC%_p0{j&+$}4@l9(L@69o5Dx`)$YLAF7T*!uT zv*d~JUg1&-Y4jjcNs&Ugo-T*?*oQA4CrFTVS-iFwo~07+yM?j$M(Rb@k5?b>An?Vp zRi6sNX^>2zuhfm47q*7L4b?e(!BHF+gPRhcdEbh!8RJ7xFIBN*%|*1f!>Z%Mubn=aQ#Ru`14@M5p`Y&mcQC5b4UTqht4sDmAmc@agjJ-0`_5 zO|^Wnj26?qJK!HVfMb$pWhx(xOR3|W;7gaqK5`Kt;*=9i2OMvoz$qDFP@Bzwm=*g_ z`W1ozCeS5B&Y_Tn%OS!rJSpHqX1|nn4YwuRSqH9g3Y_84kgB($_HI?PQL%*fepds9s`)`|OXnppJbCJzBZLd)OLUN} z+k;ZjPo4#ju8uARV$I~kS{aC;JUMj|$j`opOg(7zQysIy>P&EDB4=L&P7?~|dgkB} zBb&cgIui*+P#OJ}g|B;rLMN^aPbg_oN}(Lk$jN+EX0`U^=Ua<+42Qpl`mfacmol+b=CzMVD4}TtbEaU4=fO#xKb9{DeZkMU6FObdu zVIuG;{j0Pw1WoRupJkMc(joGJl17zuVORZxN;*Y~B@;3%(x3_CPta+63nt_~CcnG` z_m~ZM@nvVCk10-H*JFjG5M?x0EdtK*V}>7J-Z6FLl$Y6aA+a!1pQ+?#6{MjjY-6e& zS_v>@-XAQB-8Y%v$QNofa_Z%jXK=rF>ue6I4AED_Q`}VWTUgTT8{a1F{FD$PK{}{* z6rOO_IMV)X(=zGnSDh{>+!Wsw?}jkgy4Dz4y(%c&zFZDbaMpV&C=3a7FIs7b%ru$5Gm`!#VOrk-2=|LYsT>~&NY{E+$ zL<%$rTgsv*=LYnx%W5&+R#dt(m-TxZD!bJQ%%UM?v4xfV2UE-uX!U%>P(N-_x@@i! zdV~)mt$4g$0c6|{W|H!Q<`o4KIu3nCJAgC)RHgzP{( z>5`oMNR$giMgQpyy?h9Ii?P_5+4XJ2N9_W#c_2dvti!jMjjA1qM;MitpE}G>M1d`M z3-Z#k^^4HX$(ai-Tk&3c6vXzXo(m2@BQpZe5AmUmj{Qy61>TcNAgC_?0**~|dl;Bs z%G{35p$qf!nOK{e2Y%cKYZIh{^i2aZ7@DzYouv`Ri@z6>-ucTK_-TB9C{qVZ@sh;yIOZL(opfo;2hzv-^&8Re}cZ%tVr$?Yoa zYN+ya_?{_f(2VkRLzqqEc(0c!VUTB=-;4E~+B?8om{A=5jA#joDbxdV+wf^L>yg$^ zxe}J;x3wur9&G3#sPZ!({J${y$_!{Dq4%N6*p!c@82S?^$ z%L7S!a>s8SG%kVo3H%p{1?pDHU5Ucoh=uxv^=-xN(aYs+-q7>VcltNC7kdsAPE$~V z0UwzyN3>CM7~3mxmk6(*O==hLwfyOqAR@;fuu{F>)BG7PC^Y?78$YvdjbHM=AiX6o zCI!J5Osf*MWUO|>))2Oa>2!hzP-XKKq3l7hLJSAE_kcV}Bb=i(gmU(zZ?Rwc?!-Gq z0ljziceb9Z0;*7|I{GYa0UknJylWWmAu{oVCbt|LZd-%2!K9H<4ll5F zpK(7>w!+~}K)i#Y(ILc6h;}HA5Gn zi@#ub`DSjj_a$e7ss04ku+&GDWv_Fee7=KU9?|~J&fXPK&AMpn6iDeoM2eT}wc8p4 zE643Y%1Z%_6cfBW?Vswe+;b@!@TCjSzXQS>jbEt;lwNBD0juY^=gSs>!5GaTDNuq{ z!db#ud^wL_Eh=_L{}J|4`!Vs#dWV%NI3WsH0>m5r4)vPEQ2P$!Kl%XTx)WOkHoqz5 zpbuBi_kB2U9YV_@=7Mf{V4!YBJfDH!m%$sG#m2yj0PQ>Z8+HbX>=R@oRDbvA1EcyO znhm))0M=sf*)?qT5RLc&BB9IG%l~GIerrw@+Yd<>k(260jf=?{N78tHgR_e{M%m>v zuI{tHQ??l2)~Uo*%iMJ)AHBy|;-z9OB`Q=~Eg>NH6QG)ML$b>{2KocM6A=WH+3~X2 znb(6HRjhCis}&+WHbAj~_a3v_P;;R~{SAQk`+Vu@B4&XY(OU>9+pj#NmOl>dB*)1l z+M3w_mAh4AptH?!X5hV%d?N7?&7hlJMQV!~0mQqaZw#aEw?0&g5=Oj1Ozone4;Z_V zF?(l>$%eX+JwU_V;0^8=ckU~wF(%MTgQ|!YiPhy;4t@8X7tiDoRnQ|yfA*Jymc7Q> z#I`o}JNfQ4+8(IW9h(>6Tt$a#jF}d_+3%$;y zm)6N9_QhyWoLtHf!&1~uwxBFMC_v5Hu~%mKN{ zr}6}hN;;V(g@QP7leTIi?n?Qk&3IEn%5B= z@$n}Js=a(7-D`SCHCe2l|A=_I#(ob0Wi;Xq(r@&sh{83 z)6HB}!TF`FJ3Q*clbZhV#a=j8Q1Y41l?xNWE}{Q+l6+R&KCHms%q;f}UOC^$MZ~iI z#bo3eybV=v52KPSYlRw@j&F1NQVsazSxr;-0cicNmtvw z-oVSmtAAyW>^zIbEH}gAsDUJ<^(gNA`D03t@}1^kNkv!LvI@o7V0X8<;>j=5sk9}j zm@^tCrr9h-G{Z&a=f%0@=Y{zKI=C=0c%u5+E#qvA5i zn988<3iC|d2NK1KmEDb4f_x`KX5=JArPX!5-!PS7Uzq=T4GS16kjjS{4Y72{Ho1!R~@1|Q$l;m?MQ6#@^! z)JvA#<80k{st|Ma0}&P6+dJ${wMHU`IzDGM{s@jYl)a;Egb;`)$sk%VT?sb zP@7MOo|xsC^o6c)S{6me#w_=-;OeW|vr{pWBK~J=LkBR17$Yceu(1&95~_qr!jICR z;cjWV%y_WM7T^}B3B5X`!T8`20&XVRA{oS`Q;WYx=&Rp(!NRC8aWukZd7F~X>7E25 zP%K2s(*=SMzEsce0U!c)0=zm^uJZ-5FMb`Mb;O@VaVJF2ng9G@5%ESRBrBM@jAIoQPr;T1~6)9o_T{O-w-kxAR3%A@>P?J?AASCrLg=i^KN$S*q) z%l~ea+40MHT6KgB!Zg1Q+d4Kml0J8#Fiy5B)2J!sm5qcuJkn+m*I$uRjqre!tx2*Yk z)wbaP*AOM8Opcgd(wQ(v<}`GM$Q(u3axf{glo(I^4PQ;{4o#_TXPD%y!3ikKehr*Y zR^eWH-A#tlS~bqv`VO3S{#9496^2Woc}8wN$~nDznLyyf+GHWcQ(Lc3>ENuF=%taB zWkmS}WTw>r+f54fL1Uy*zgiD3K$L$Icz^^Cfj_;z)(H_yvX?TJZpH;oa(LqC4yL zz_MklF`b*peUkYRbk_Fr&EQ+OGntErp=r^`il5w+lH_2A0Mr>3&1LZ`#^9ZOvh)Y; z+OO88T1Y4({tnM#qH^Ed{Jxp#yV>PCv#iS$>hw?I&q2Vcp#ugwuIPYCp@rg_dfF91 z^zk$w|5}h%$%cwuRaOq4ea33Lf0!=>Vg zD}O2_2C3xE-KyaBv6pe$q>7~|itDHQQ^CY;o?ZLxOs}HE{Nr4C`Pm7hL!|I>?iTdx z(0p)ijdv-wkd9hLVJWl;!GU#pK(!|L?d|V!l>>?%e7r5TX?C0hB&l}H(u90*)gMQ< z(`V`Jule=VElaYiQMSuN1y7vD>r4idGgv)sRe4ukvDU9lPEPw_-7S^pf$pdo>T8bA zFYcys4T$$=Ijf0~Ep<~fN?~m3s!$z1GO!KR8jlpVhZncSEgQD=e(yB9mor`-Dcsw^ zu(AX7@zrj9lW;T@SYH>#M$h1&?tWn1w?q>;<*97KECa>*AvFc|Nt$JKZF{&RWOA?1 zP*-G~zno#(kkiceYLR!xPVuW4D)Wf$H^NPej84UiAAbZPO0g@7=&z&~4qp_{SwTwS zh7XD&&RHQ7cfxa9bBG&OClwE`MJ-%LCuw3m=e+NZyW8#{jc3^OjxDJl8Q;hxYiP8f zy4&BsU52t|{I2o@Z-BKK{)@8^`M9(Ejap4ZSdp=8`W zP>evp$WG}EeJIm;J|@IM@F;iGq`wy1Q2!7-xqR*4x0PCxfV$4Cua&@Nb>lzNjBv|7 z%0r))&fLtDJ)@IAy+;y8q01ArdRC}G5;PutyzS*KVfM@Ho>al;muDPhFYT^q;9zbS zO{xG%;Mulgj~=Zocyt`%A#!wF?B4hI?;^fd?(zxi<#{Itgr)?q7BK|{D!$y~Yv%h_ zEoQvB`zslI(IU?(&tA{$Ef#n^yr#vyNw8Q9)1A34O5^O28;_L(85XqFtC| zaqoE&u4*kC3>u!{vAeoYgpUdKLo?9YtUH$@8q)gXow@tlGRZPsu0WfUHXIi0T~A;O zM@5Zh50$hsgoZXgQ7}x2JCuP0Bd9qwYeKL@hGuVB964qmacprMZ&?8{O1e;-zx_Jy zr`f|a1WWmnQV?Cd5OpLWBKR#DNLM?7!SJm%sWv7MGUtK;r`D>Lb{R-RM19ohcdR+o z>Bm4#oNTj7F=2HK8B@@b7h;GjMHk;zeh6$OcS3}sw=w9yyIW?>Uq#?_y|4UL{q#$% z=}mrXEp(q0@cx2&nh^?D89xt4I|cb`DHQgor7YNYeTpI~(WiPYSN+(v!m;F*j2sYr zmd$3*$ZCA!Li_}p*{HP!Sd&JKTl;kmNnt}`%p3UBIm z7iD+o@ud8Grmfx9Au70q1G?Bv-yY32*#$f9+XW3Ya=%Q(@88c=v3dQpH(t->XMQY* z<6Z_ACo0sRg_KWF8bU;SksRnrM1`3gy5nbsUAHXSRaS^`eY;(Lc7$7BPCqa=FpdqT zX1%KjFgRgT6XH2r;IwMwB!SoZL9A+xt_eFo;V z$TN@q)y3~FPP%{-DU@g!wPB-9?BlCO({e^NRPUhURBpX!x`2pXa9NL98kr$qGxhY~dVL$vyT8}D@zVKf zd#+#0Bj2OMJ-_SN0x7M=GL87zL+qXuwqe$&0a`G#*9w4hCN#vCDSIoJvkvDJ)V<%} zN|+-hn`H(5fZ`r22>4}+6yk>CGKb(yH*4F|xe~!88sJwF!3?`c=ag7snu7O=rGP!? zWzgn`WKAaN=?sJMAW71mcDE1}{vP=(xo+6^#-Nms5jA^0hg116cU|T&U}(kfvfT1k zOPMh=>#tNJS=`S})g0)Ei>vWje3hR9<=Ge*kzicLgE0@WFW?MvSF@%fJpOeS9WP@w?Gc2vIsiG4ljxXr2?Qrcnc3%b#9j}P)Bp=ae zHYr?VN6@!r&?H7E=t!G2e_KSs-%sqmwT*u4>_JC;+1ViGG64Ji{ra$AWp(ztchb2MUlQ7!8IJ2iwr+;n(looV*UVr zY8HDwz_1$C79-b??yjn7r&7&Cugu8<2sKN2b3Dv7lChSGEX@(sLm)1$Kz;>`Nvchh zCeZ;VS!tUPI;6AT0$E*wC2{Ws%y+whRmy2ugpJ{8oV{ubZky=_ z8P?0@xk=X1`%JGZ!rg#XhbKuga;m|N0K1=XmUxRda}tOw`zv@|!HQ11LK?5pG5&*1 zu?0K}v6Eeec*~aprT40(Zuz811bFuAwKh_CXGAI!>A@FpZnIH!gBo5fx5c^Ro8P`oa6&v(9=&j~44^vRYSr@lAo0TIQ{P#$T=`|w8jV3k zzEG;r9Q?%-C%X}gImIb0mfBOGy-6Sv6V{?KTkknT#G!qLzO&i{!@!abMI8VbWSxPUuBX$SXlow_g9ui7}C{Nzw|rHRkTZ7)BJ>=PF~tt48cV+&I2PAra` z2YyG#`gDR48UlLWIO|b(!##Zryu`8AA4NLIyz8( za(1E`r%`7(`4T+LiR!4$2_$)uHJn;P{y!O{e?S_4cpElmCWe0^NFPDd|E~zrM^g5G z6^+f#@*kqHnOW%mLXha0KJb{2Xzc%pAn{9y$tx)S4MF-V8~Z zqUbNy=$~=e|Kf~32&j+F{dW%P-`eWkD#7fw4zOT_54oJmGVeHrs4hv;}5J%1j*Z(DIcvpaJ-W8crP*Wp+KH1I)MfWUS^ zC9qb2o>t}KyV!Q7z*;j@2;wtQt2ep0GFV%5y}a}e&rwV4S-yX&9M_zhYNbmK%7jm# z>hv|tKqeefuWmK4WU%)9rgvJ60`>LS(wGiUb*?dPejvyzN-j|?$9jN2XIuFyo_5G-{e$Q10j#E^oGqLbsP~6IXtd5+9VR+!v zM@hsrzaQD`Kb-!FQ2_2X6xv1amIbgyW>_GsI#FyRpX&>t2B5p3M_eL9 zU(;E zN_^c!81q|uwyPJ8i+z1|Vi5m*)_%111LE3bq&MyMBmgOH5w7BZVouMd%kTGHJbZ$=RK z{T!4TE;~d|(5fsYvNKU5KLHUR)KGvQoCfVIXEq0_f1uOglDf}_S|rMVXeQP@oI{2r zH!IV)J6_iX{>ysoa(VlYvNmPHmO%_0(7-TCDp#gt@@KuN8` z7%k&ja*kLU;llYthn$Fcd7Oijsvglz{8;3cB0Bsls=c~e0X0m*FdsymfMEDKTgO_f zJ{JWv?xs&id+xcsBNzQYz0?JjE+R?g&2eT%vCktCz>f75^xbNiRMFxKw6P9Jb~4603_&z{_D%D*)U&Q&E`rTlsqoykLCQ3mt& zTHz&C`(<_1%(HSMR*@GGO~?7AMvH2>xBBenz|)cO`()wtac%%rv{*$&=vm5UI@TfK zc#PkyvgZx}dPelAk?YFh)bN+ z4BR43@cxajPWLru8(;$%rK9H^HT>yr{|SPDsTF!&G;!1M{<);2`4yc2_HMRf^M>>} zM5U6TMHxlj1DyaeB;;3ik$<&S4l{4=4qOg&i$8-DR`9e@#!299z$$thnr)cH902io zeOJ8?+;6rl3{W7_>`2$%j!4>{rKrZPvh18|e1c^9dZBS-PRo)X9r0K;SG_ zj1Q(Sq15yJ+n4Q_NuGOZ(^nU4?1@Un$^|yncwQQB>=J=9AWm}S0+;G;JFu3vEE5~O zRH1Ug*FclyFd}01;YUkz-O-fw>_&|nyc3@mle35x6H^(%6t(#5q>Z97cLGd@su7)1 zKO$*tKzC)OFM-_(agQXhn(OD!88!~$NGW-7)5_djo<_*RS~KN-AQF^QWwkCiWYGj1 zF7-)}63-{AA3cTWJO*@YDGxETy#S)O&^NS@I#t9LyuB<7bdd*qOP%8mAlJnU^r#sp0~{C=Cb z*i?XSq2bXUzG~Z5Co}Eaw|Y@GjQt0w-&v%DyLgg)Jp>>5&tFXpt@e-z zHp7joCt3A>vyM9^Zt%^nF4{{uCD8-upJ%Ena8-=>`Rx(T+UpE0Zr;+{Eix*x(Fphn z`e=M#l>89AAel?^UL0Ebez^V$>MhnouDA zN)&4?KJ0cYgTGT<@;HMZ$4_^yhja}G{7V`%kC#wl`k6t~)fBTBw zS7>^^xnieB24)b4)=5vAX<1Xa&_!8BT%SJD#=z=z>*t9-Rg^F>5v-IJ*Y#vZo-hbV zM8&#)K-n@m7DNvYW@OtpRqo}gX>Dj8Xdj{GuJMMWqTSo3RZ_@=Lv_hbm~ms$I0uAo za-@-7w`q=DdziU?A!?t9E`Yt^cpdNg*3W>1_zmrl%vw3BS7GzO_lO3zlvkOu+KwSx2yX^I z!ILX~A%4Lp87GCN*H*9J$3J28=Jw_eC<;NgJe?)mNH$kG5}U=?4%}X(vdwc%ccoAj zZCnLx=eWAgDe!uIGxQ+xl{J1y+V0zw7z{kugWaPN~jof^6(_&b((>c1vsuEUK{W9mZ5Kn;*)h zcEp&xc~`$?FVsCI^Oc*8+`rXMKCkp7+xK#QD7;p0qjlbFD$w}>xdKX1iIn#A?VZ>*Be&c>Vcyv$0R5>|J_^*Yp^ zyPs@GTe$e9Bbuokl=l5uag^uLD8=iUs*-eSnL?FtE#3SB@127O>Kd<_L`4af#YpVf z)t(ZQQ*!0S!6Is)VV(YtQ>*pynd3-;E)+ol_cpIX*3D4o?onB{A+nZxB~QwCzD|f* zSGqhWcxh;y{@zXJwZR*rSJf`z@H%V#`Kp(;A|sMM!n^*Cw)d4at%*tA(AR|35~7_Q zemHQnOTXAJ(UTwCnJQYowPrH#z=zyP&d!>m;Mvv=rXc8$?3TFQiy9ZJmbLJpfuo&n zZOOWoK7^0aGq&X&L%sOyk7ak`hL&Nec;p=KJIQ(E7A-Q92Kum1>>n)7WH6?aVH+Q* zH-z>4@gcXUHQql-0)I193({&U%rY(cCNU8-Qe-$b^fbbJK~8~;vWl{?X3ML|aWDI1 zwj=IKGa1O>fZ;9+GD?w%ilv)+u-_DzmfpE9QQQ*Qf`}U4WZK#qX|W0Q%MXZB#g8K` zkg3o+F1mzy)!TO8&48O|3q(H2#}n7(W#+3-`iQX4DNnid%X}+s8@EpNqIq}z7iC|I z>33aR7iG1mYqud>Dh4}o>guycCt5YkOG{S!OMOyXI_@yyD!a)pF6VYn88yt!5#oX; zFFUR#%iOR>h85h?=W(8$2vTHE3B^bqI2?u2vhha93s^{(J@T-BqM+#rF8_Gw9T>1F z3JqgvuqxwQcd2f+l#<24V@u1p@48KSz*Ueu&_3KPB`vvR1#`Xs{aI0oIOWKg-rU}Z zE6fenMwb5N@J(1_b+m|h_BTVl$PTtwd6iNpFQ40V?|K?`kbwyflHTq+iH>1wt17+h zu?gAE#ao;dfDB0DmpRaym}EzPrC+7tk=^Ix3t2#CoS@{^OLF}2?mkA}2wJLaxIQCC zq&=@q6PWubAH?0_mBv|W@_0OI?*h%`sgO_)h8Z;v5oX2bvnLB7ubL+EMefy|d@-xh zzdZ$9>>C-LL$o@5y%0g8WI8e+NBP`0zL%Xw`o0F9QSZXN*9W;r__;b=WYfv3Lu0PO z&IJ?CdZnd=?v5@{-|sf-C4bFK#EO*yV3*9?En@NI>0mhk{>YjLYQk zl7NP&7Lk^cPD4owpV~!6B`#4jKJ(h%*=h^qqbZT2`304q6jYxtc^t$?)DhRl-0Ckx z56aK%ea4Dq0Bo6mXO1wZ*eOmOfVhs{w4o7U){az|NIoM3IpzJLNY9uNZEeNXO(&svWPbCrB1SI1m^BdKU!#+g#*8GC2T#s>~K8@2c={D z3EWk7kzT@oftMAX!xOfcpZwGDilwC`rnJkELc!N6}%XXSY%_#}l|Rn7_w4zkJAhX7uv1fU+Guij1Q|Jf>nC$V@SWQR z?QI*uvuh3R%jLR?cBC;!!Zi(zzKi&?pEmK{wxKLE(#4c3N%M6ACTHUK6zV@GD)2`+ zS4D-;bzf)UxO`Q=Nh(L=xgN{BnR)_??d7Sb#gB3hx{_8pxvj8j)y#Rd7lQH#Lppe6 zj75qi_ec!|NAt;$?+5~jc6qAE3@oLY%#L=6=OpT}>T>zE>EgHq;2R@X^q9h<_Z58N zqGP77xft$j)f~+I+=mCOw6pmMMGKJ9l4WnO9NpdiEc^E9`g8j&`IjQ&pHv{Zb#azI;-##)7fy7jgZg2>3R+fENmE$+$7+sf!D%t|s3nx-b?2=LMwgVLcjX}5Lc_h~0Hff-$uXgzuC#A7RQewq z2Hi`SQ9iuS!aQV1?=e84@w@c4zd6nsq~F|yPl}YS>Kc4sOEXgj<*U+N@8EQ>)nX7v z8+CW+am}67S}M%Awa$}Q-%bzsm}mm=n;MKT&|!?sl$+xlWf)0zDcfze-3zZzZ-39- zs$p0@ZI$aXsGpcV@78fyK>WGjif(%%vjmese}-k689?t{E`A-2`?}6B;DOfm=51qe zhT7(f_Qk|Yl3mw<9IST5io7LAi=*6F`F8jCp1daJ{6~uwE*(vXNFLh2J~_j{p(-m% zN&$L}hy8Xj>=|ov!;gglfk>{Vj!FwIc~B8%6TYY0Hbc z^1f|X2I=B$+-Hb`RkOdQl$l_tWnv|K9_K|xd{lJ}ajHhNnpSjj#1wSIjbUm^hD^OX zDdZn0f_fp(fNP+wFEs4}-!jAUZvq zL}mre`*9|v5uy%7oS_BWG+zNAr36g|&+k7!y8FG@O-yp7adlhAY9OmG zA>)_*J2c~gkwRZ@xB$Z~b1GLd3WU@Oj7h;dwkV_m-h;=o?GpV8;BceOTxQS)V*_)g zl#?CbwUh4$P+nf&h$dc>ebk`p^rkzHjSyY$6G;0=K0TXUiQ5)+XYa%YcWT99cnzzn zu5hEbg|MC|v*l8cFgr`6jfeK^*z6I16< zGE@i)#N|L3H13$iJ%M9D&P&RIl#q5Va>4F72f6E6Er6X#zU1w+Py{vX1%!+6o6Dd;Z72i<%O}Qk&ENu-#_S=*WbZc zkbeBgU%ia@g?enaiN2OV(rtVBG}^c@X3>ux-#50~wB7rQafZ`x4go)o-9UKE!Ke-jLoj>d01kDg821zlN-zTfBk2ayG86mS zon&hq{Y&XG1Y^xJOtvqfAkR$Mm9qW zZnc#-!8VpJMKW)S3eqQ?#^yqE{O2X?+RG#FJgFlvOgCJs-G6ZBQqfl=-(f`r17ViV z+M{5#H4-6fW1UNk!WV`6slV1zyiqhMCGC=|dXh*kCHd6?dGp}P&Ouv;82Y#-GSYD+ zt;{tp6S0gv^;^yJt zd8711F`S<~>mb6>@qV%uZVqLR{B9AF58~A_@ACo&G4k47bhjB|GI5mtvh;ZEV0^kB zsX6^3TL7%1H7uLFW+?EpAU|I@sS?T4qqo3jwZ?xU&fD(wInkPzIm^?xeLr$g*!BJ9 zfc&5n>uhBDt+a|Q;)1mjHRWh4m~dqq0|t?u$*OF+3U2SZ(q;`t|$Xc>wxE^ zXQyV;T-_X1EPS2S4nC^>wRR$*Zt^GQ7*L4Wj|`i%(Q+te@29&YGjF5!eltXUD)*Xn zEMWxWoXnh{j4?XA0!*NZQ0B6cDD#Hb|GW0FAhC=JAw z8Yj^|Iyy=c>1dPYyeu}-gXB|rLEn5iwdD3=S)Cdo!M$A*5BROIIAa57Vz{IJIWO8w%z9v>zgX3bezX&Py&=?FYUxp+0)7e$U|(4G#H2A#Hg$%{K!=!0 zCE@pLLq*?diZ}qNw@OkY-;-LR2+ZfbCi!{41uc7GjX##ok0NAPHWebj;Hqe1YE zkyP6ZQpFVnjoxv4@%D^7KL&tr6j)T1nOi_Qtl%ZoaBaLuQT4cSxj*NL`IEermYhDh zx0atcMr9Pjna#Oo<=ocMK;xa6Ski;7LRFKB4S&K$Yazpc=gJ>#Zj$y&e6-GI~xaesp(O;*5Ee33ufK?^}5`VW9~^?n)u=!*|Oz;XKyvD zF|lMjgWN{=W&Pa+E*ez*Zgnq`3VGPnilztNIDzh|P%reep5bM@!8@Hj5Ql=yn&Lrf{RlBq6Nkj zPd|41SwOh+;Y>rW;$b>!GO=#q!(>!RNsW)%^Ayr;9+O(1oo|ts;e%~os3Yi3K26W+ zSOh$_p;VUrp4vocf7ldYJxV0>ZlIoqI8k2nONSCCfh~y_t}Y6N?&?V?D@XG^HP06l zNPzAVKV5>5xV*iVnA~>_mh$FW4=;Pr!`?X@`BqYrkJ{0-XENI@7DWtqHX$62kjW= zcd1+tzq*tAN(x*y$%Cn`3gT5$DmP3mGf^hLii%7daFXpF3ZAZ1-38m4*~!6PsevPU zFtzhYTtzQ1E{inzwS>OUt;hocp6sEm2>JERx0uSS-*?5SP$A4A8=OJ6 znW`M^r#;hz#4Ta$jl1pHSIy_=D2`ZE7brbuCV{1m*hhch0Huu}mw`gEgKC8N{t=xJ ziv`AhiU+}6?V!3mcv)$((I6DJQE-pGPJ&y-6VIkkuOsz(8;x-K?5?^&tay%o@9owZ5v%9W9wwLKL2NgUjM%S`V&P;kjb8PGhNw@=4T&^ zVC&a64sR-Dh}f%iZ-1`c9r`nT?e+Ic@#BP|ImG z>Qnn)pQqpcvj0s}IJNy3MFCqo{r`cY@Lz7rd3wSx<(c1AtNH#eEdAe%+(g<#PP4#{jKrBj+IpbsA>Z4VtFukV0XX(bvH9- zS5JVNl#C=+C&$vm^`Ab!RaV%wENziCJjzH}tUeER8hdU4mI#BkF>|uD{FgQK)X%2` z5B7hc-^3nvWm}{($^~WR3eZJ4Tl{12oAkmiB>)CuwSllz>w!Qp90bw_0@<+de=L8~ zF#c`y-)GbSz|Qf{UHm=gUt%`4uGsSSe1EIZh`YJkpqv4Gzo<8-#1lJ9q${@GI(Fe$ zs)(91%EAqa)r#T!Tc$=5D-q+0awdSD4v97RtBrN46jn(F06kqg7Q}k{Ok2_Nfjf4K zE&%B1im?Xb;%Z6)zfl6{FY3rY3Pc{xmR8tf0EH4@<&pk=0bo!lR;tJf@C$>(_z~Ff zIeh_~{>Fe{#A%`WGZ;3=PfOXK!LW_fBXAZ2^CQmp1ww$J^K*c}kkdFl(-+DQhhzP9 z2IGSQPlc7vU@(3l|M|W^2=*SIvjZZar>vnfeL)cLDShWG27>_4+W`@PQS4uc*2LWT#?zz4&Fh4dN&e&mvtoX3E@eBrn{pl;L9T*z~XY8;ZI1Q&W z7#zlT%3eByfx$r79~kyzoZlPP?vzt^t}h7o$C|+q*!eYs!F-_e7#|-PnVGigJNmX=Y5aGPLIKvHADHK!1HrpBjfbwp0R_$q5S9f3xmV{h$|SJ zA9j9T7(e)r*o8w8fA|lJ{o#8!42*p^pIr+ah6Vq`{zG3l=zJW&u^TxbYgp_LKV#4S zALks4ohFwve!*h!f7`FCvze`fr8B{464bQy#HK4OQ&Izkas@!J-aAc8@=jJL0Cx3% zT}uHB076z^ONfP~g}H?#(h7-0AP^t~0*W+)SeYRa=0LEO2*Lks0x