Commit c27f7075 authored by Marco Govoni's avatar Marco Govoni
Browse files

Improving tutorials

parent e8d01ac9
%% Cell type:markdown id: tags:
This tutorial can be downloaded [link](http://greatfire.uchicago.edu/west-public/West/raw/master/doc/tutorials/west_100.ipynb).
%% Cell type:markdown id: tags:
# 1.0 Getting Started: GW calculation
%% Cell type:markdown id: tags:
In order to compute the GW electronic structure of the silane molecule you need to run `pw.x`, `wstat.x` and `wfreq.x` in sequence.
%% Cell type:markdown id: tags:
## Step 1: Ground State
%% Cell type:markdown id: tags:
The ground state electronic structure of silane molecule with [QuantumEspresso](https://www.quantum-espresso.org/) is obtained by running `pw.x`. The pseudopotential files for **Si** and **H** in UPF format can be downloaded from: [QE-PP](https://www.quantum-espresso.org/pseudopotentials) database, or from [SG15](http://www.quantum-simulation.org/potentials/sg15_oncv/upf/) database. Check out the `pw.x` [input description](https://www.quantum-espresso.org/Doc/INPUT_PW.html) in order to generate an input file for QuantumEspresso called `pw.in`.
The ground state electronic structure of silane molecule with [QuantumEspresso 6.1.0](https://www.quantum-espresso.org/) is obtained by running `pw.x`. The pseudopotential files for **Si** and **H** in UPF format can be downloaded from: [QE-PP](https://www.quantum-espresso.org/pseudopotentials) database, or from [SG15](http://www.quantum-simulation.org/potentials/sg15_oncv/upf/) database. Check out the `pw.x` [input description](https://www.quantum-espresso.org/Doc/INPUT_PW.html) in order to generate an input file for QuantumEspresso called `pw.in`.
%% Cell type:markdown id: tags:
Download these files in your current working directory:
%% Cell type:code id: tags:
``` python
%%bash
wget -N -q http://www.west-code.org/doc/training/silane/pw.in
wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/H_ONCV_PBE-1.0.upf
wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/Si_ONCV_PBE-1.1.upf
wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/H_ONCV_PBE-1.2.upf
wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/Si_ONCV_PBE-1.2.upf
```
%% Cell type:markdown id: tags:
Let's inspect the ``pw.in`` file, input for ``pw.x``.
%% Cell type:code id: tags:
``` python
%%bash
cat pw.in
```
%%%% Output: stream
&control
calculation = 'scf'
restart_mode = 'from_scratch'
pseudo_dir = './'
outdir = './'
prefix = 'silane'
wf_collect = .TRUE.
/
&system
ibrav = 1
celldm(1) = 20
nat = 5
ntyp = 2
ecutwfc = 25.0
nbnd = 10
assume_isolated ='mp'
/
&electrons
diago_full_acc = .TRUE.
/
ATOMIC_SPECIES
Si 28.0855 Si_ONCV_PBE-1.1.upf
H 1.00794 H_ONCV_PBE-1.0.upf
Si 28.0855 Si_ONCV_PBE-1.2.upf
H 1.00794 H_ONCV_PBE-1.2.upf
ATOMIC_POSITIONS bohr
Si 10.000000 10.000000 10.000000
H 11.614581 11.614581 11.614581
H 8.385418 8.385418 11.614581
H 8.385418 11.614581 8.385418
H 11.614581 8.385418 8.385418
K_POINTS {gamma}
%% Cell type:markdown id: tags:
Run ``pw.in`` on 2 cores.
%% Cell type:code id: tags:
``` python
%%bash
mpirun -n 2 pw.x -i pw.in > pw.out
mpirun -n 2 pw.x < pw.in > pw.out
```
%% Cell type:markdown id: tags:
The output file ``pw.out`` contains information about the ground state calculation.
%% Cell type:markdown id: tags:
## Step 2: Screening
%% Cell type:markdown id: tags:
The static dielectric screening is computed using the projective dielectric eigendecomposition (PDEP) technique. Check out the ``wstat.x`` [input description](http://www.west-code.org/doc/West/latest/) and generate an input file for WEST called ``wstat.in``.
%% Cell type:markdown id: tags:
Download this file in your current working directory:
%% Cell type:code id: tags:
``` python
%%bash
wget -N -q http://www.west-code.org/doc/training/silane/wstat.in
```
%% Cell type:markdown id: tags:
Let's inspect the ``wstat.in`` file, input for ``wstat.x``.
%% Cell type:code id: tags:
``` python
%%bash
cat wstat.in
```
%%%% Output: stream
input_west:
qe_prefix: silane
west_prefix: silane
outdir: ./
wstat_control:
wstat_calculation: S
n_pdep_eigen: 50
%% Cell type:markdown id: tags:
Run ``wstat.in`` on 2 cores.
%% Cell type:code id: tags:
``` python
%%bash
mpirun -n 2 wstat.x -i wstat.in > wstat.out
```
%% Cell type:markdown id: tags:
The output file ``wstat.out`` contains information about the PDEP iterations, and the dielectric eigenvalues can be found in the file ``west_prefix.wstat.save/wstat.json``.
The output file ``wstat.out`` contains information about the PDEP iterations, and the dielectric eigenvalues can be found in the file ``<west_prefix>.wstat.save/wstat.json``.
%% Cell type:code id: tags:
``` python
import json
import numpy as np
import matplotlib.pyplot as plt
with open('silane.wstat.save/wstat.json') as json_file:
data = json.load(json_file)
ev = np.array(data["exec"]["davitr"][-1]["ev"],dtype="f8")
iv = np.linspace(1,ev.size,ev.size,endpoint=True)
plt.plot(iv,ev,'o-',label="XXX")
plt.xlabel("i")
plt.ylabel("PDEP eigenvalue")
plt.show()
```
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAgbklEQVR4nO3de3RddZ338fc3yWl76IXQ0oYmUAoIlVYuHSIX6wUKWmQQKl6e4RkVZ2SqS30WOGO11VnqjDOC9Bkfx2dcs1YfcaxrmFGWlIKKVCggeAFJKTaENpRboUmbpJQmLU3SXL7PH2enJOk+JyfN2Wefy+e1Vlb2/u19zvnu1XR/z++6zd0RERFJpyLuAEREpLApUYiISEZKFCIikpEShYiIZKREISIiGVXFHUAUTjzxRJ8/f37cYYiIFI3NmzfvdffZYcdKMlHMnz+fhoaGuMMQESkaZrYz3TE1PYmISEZKFCIikpEShYiIZKREISIiGSlRiIhIRiU56klEpJxs2NLCmo3NtO7vprY6ycplC1i+uC5n769EISISk3Q3+PGUA6xe30h33wAALfu7Wb2+ESBnycJKcZnx+vp61zwKEYlCLm7uQ+XDb/AAUxIVfODcWn6+tZWevsER5VefW8svRpVPqqwgUWW80TvAaHXVSX63amnW12Vmm929PvSYEoWIlLrx3sTTvQY46uaeTFTyoQvquGtzS1blkyqNSxfM5jfP7aW3/82bfq4Z8NKtf579+UoUIlKuwr65Z7q533LdOcDRCWFyVQWJSuNgyLf3QpTLGoX6KESkYOWimWfNxuYRN3yA7r4B/vPxV476vO6+Ab70s60AHB4Y+W2/t3+Q3v7cXJcBtdVJWvZ3H3Ws0oyBkC/w6cqrkwl6+wePSnhDNaCcxKsahYgUovHWBMLKE5VG30D097jx3tzrgiQ20esbXgOa6Kgn1ShEpOikqwnc8fgrjL71pivPlCQy3cSB0G/76b69j/fmPvxGHnaDrz915rjKIXcjnMKoRiEiBWdg0DnjK/fl7P2SicoJ91GM9e39WDrMC4k6s0XkmOVyOGg25VefN5dHtnfQ3HYgNJ4Kg8GQ21a68rphfRUTHfVUiDf4XFGiEClR+biJ56IdfTzlADOnJrjmvFp++uSrdA+bN3CsbfilfIPPFSUKkRIUfhOv4APn1XLP060jxuhPrqrgsrfO5uHtHVmVJyqNd581m98//9pRN/F8qD1+Cr9ffXnOEp6MTYlCpAQtufWh0A7XUjDeyWIycZkShVaPFSlCnd19kScJ480RQKNVmkVaXpvmcyUeGh4rkmcTXU6i+rgEfQPpl34Y75j+dOW1ORzrfyzDR6VwKFGI5NHofoWhlT4bdu4bccMcvgIojByq+fqhPszgqredxMPNHZHerHM51v9Y5gZIYVAfhUge5bJf4ViGfaoTWNJRZ7ZIgTht1S+Pmj18rNThK7mkzmyRAvDocx1pj6Xr1K2rTqbtUFaHr+SL+ihEsnAsM3iHl89IVtHZ3c/cGZPZd6hvxJyFbDp1wzqU1eEr+aJEISUt6kdNpjs2unO6s7ufCoObrziLyYnKY+rUVR+CxEV9FFKycrH8xOSqCqoqjDcOHz07OVGZai4azzLW432YjEi+FFwfhZnNNLMHzGxH8PuEkHPON7M/mFmTmW01s/8RR6xSvDI9sCbb8t7+wdAkAakEMd5nHbSW6ExqKW1xdWavAja5+5nApmB/tEPAJ9x9EXAl8F0zq85fiFLsor4pZ+po1oxjKSVxJYprgXXB9jpg+egT3P05d98RbLcC7cDsfAUoxW9GMrwLbrzLSVQnEyQTlSPKhjqTVy5bEHrs+otOSfsakWITV6KocffdwfYeoCbTyWZ2ITAJeCHqwKQ0bN75Ol1BB/JwmW7i6cq/cc0ibrnuHOqqk0fWPxpaunr54rrQY/+0/Jy0rxEpNpF1ZpvZg8BJIYe+Cqxz9+ph577u7kf1UwTH5gKPADe4++MZPm8FsAJg3rx5F+zcufPYg5eitvdgL1d/77dMTlTw6feczvcfekEzlEXGUHAzs82sGbjU3XcPJQJ3P6pObmYzSCWJb7n7z7J9f416Kl/9A4N8/PY/8tQrr3P3Z5ewsHZG3CGJFIVMiSKueRT3AjcAtwa/7xl9gplNAu4GfjyeJCHFbyKrqw6to3T9hacoSYjkSFw1ilnAncA8YCfwUXffZ2b1wGfc/UYz+xjwH0DTsJd+0t2fHuv9VaMoXuOd+zD0sPuw16hPQCR7Bdf0FDUliuI13tVVhzqrB0P+jDW5TSR7hdj0JBJqvHMfwhLEsb6XiITT6rFSUNJNSNPqqiLxUaKQgrJy2QKmVI38sxxrAlu6SW+a3CaSG2p6koKyfHEdTa2d/L/HXsJAq6uKFAAlCik4M6YkANj6jfcxPdgGjsyEDpPpmIhMjJqepOA0tXZx2olTRyQJEYmPEoUUnGdaOzVZTqSAKFFIQek81Meu17tZpEQhUjCUKKSgNO3uBOBttcfHHImIDFGikILS1NIFoBqFSAFRopCC0tTayUkzpjBr2uS4QxGRgBKFFJSm1i7eVqfahEghUaKQgtF9eIAXOg6yUP0TIgVFiUIKxrY9XQy6+idECo0ShRSMptZUR/bb6lSjECkkShRSMJpaOqk+LkHt8VPiDkVEhlGikILR1NrFotoZWJolxUUkHkoUUhD6BgZp3nNAE+1ECpAShRSEHW0HOTwwqDWeRAqQEoUUhKbW1NIdi1SjECk4ShRSEJpauzhuUiWnnTg17lBEZBQlCikITa2dnD13BpUV6sgWKTRKFBK7wUHn2WDEk4gUHiUKid3OfYd44/CARjyJFKjYEoWZzTSzB8xsR/D7hAznzjCzXWb2b/mMUfLjmZZUR7ZGPIkUpjhrFKuATe5+JrAp2E/nm8CjeYlK8q6ptYtEpXFWzfS4QxGREHEmimuBdcH2OmB52ElmdgFQA/w6P2FJvjW1dnJWzXQmVaklVKQQxfk/s8bddwfbe0glgxHMrAL4F+CLY72Zma0wswYza+jo6MhtpBIZdz+ydIeIFKaqKN/czB4ETgo59NXhO+7uZuYh530WuM/dd421/o+7rwXWAtTX14e9lxSgPV097HvjsCbaiRSwSBOFu1+R7piZtZnZXHffbWZzgfaQ0y4B3mVmnwWmAZPM7KC7Z+rPkJhs2NLCmo3NtO7vprY6ycplC1i+uC5j+T/+4lkA/u9DOzg+mWD54rqYr0JERjP3eL58m9ka4DV3v9XMVgEz3f1LGc7/JFDv7p8f673r6+u9oaEhd8HKmDZsaWH1+ka6+waOlCUTlXzogjru2tySdfkt152jZCESAzPb7O71YcfG7KMwsxozu93MfhXsLzSzT+UgrluB95rZDuCKYB8zqzezH+Tg/SWP1mxsHnHTB+juG+A/H39lXOVrNjZHHquIjE82ndk/AjYCtcH+c8DNE/1gd3/N3S939zPd/Qp33xeUN7j7jSHn/yib2oTEo3V/d0G9j4jkTjaJ4kR3vxMYBHD3fmAg80uk3MyZMTm0vDLNIIR05bXVyZzFJCK5kU2ieMPMZgEOYGYXA52RRiVFZ/6s444qSyYquf6iU0gmKrMuX7lsQaRxisj4ZTPq6W+Be4EzzOx3wGzgw5FGJUXlubYDPPny61x61onsaH/jqNFN9afODB31lK5cRApLVqOezKwKWAAY0OzufVEHNhEa9ZRfN657kide3MejX7qME6ZOijscETkGmUY9jVmjMLNPjCr6MzPD3X+ck+ikqP3xpX08uK2dlcsWKEmIlKhsmp7ePmx7CnA58BSgRFHm3J1bf7WNOdMn89dLTos7HBGJyJiJwt3/1/B9M6sGfhJVQFI8fv1sG0+9sp9brjuH5KTKsV8gIkXpWBYFfAPQ18cy1z8wyG33b+eM2VP5yAUnxx2OiEQomz6KnxMMjSWVWBYCd0YZlOTPsazPtGZjMy3BxLi/XjKfqkotDy5SysYc9WRm7xm22w/sdPddkUY1QRr1lJ2w9ZkmV1WwbFENG5va6O0fHFm+sIaNz44sTyYquOW6czWsVaTIZRr1FNuigFFSosjOklsfOlIzmIi66iS/W7U0BxGJSFyOaXismR3gzSanEYdIPUJCT5opclqfSUSykTZRuLseYFziaquToTWKSjMGQmqa6cq1PpNIacu6F9LM5pjZvKGfKIOS/Fi5bAEVo9bm0/pMIjJaNqOeriH13OpaUk+hOxXYBiyKNjSJ2uVnz8Edpk2u4o3efq3PJCKhspmZ/U3gYuBBd19sZpcBH4s2LMmH3+7YiwO331DPRafPGnFs+eK60ASQrlxESlc2TU997v4aUGFmFe7+MBDaMy7FZdP2dmZMqeKCU0+IOxQRKWDZ1Cj2m9k04FHgDjNrJzU7W4rY4KDzSHM7ly6YowlzIpJRNneIa4FDwBeA+4EXgA9EGZREb2tLJ3sPHmbpW+fEHYqIFLhsahSfBn7q7i3AuojjkTx5aFsbFQbvOWt23KGISIHLpkYxHfi1mT1mZp83s5qog5LoPdTczgWnnqBnSIjImMZMFO7+D+6+CPgcMBf4jZk9GHlkEpk9nT0809LFZWp2EpEsjKcXsx3YA7wG6A5TxB5ubgfg8reqcigiYxszUZjZZ83sEWATMAv4G3c/N+rAJDoPbW+nrjrJWTXT4g5FRIpANp3ZpwA3u/vTufpQM5sJ/BSYD7wMfNTdXw85bx7wgyAGB65y95dzFUc56ukb4Lc79vLhC07GzMZ+gYiUvWz6KFYDjWZWm8O1nlYBm9z9TFI1lVVpzvsxsMbdzwYuJNX8JRPwxEv76O4bYOnZaj0Ukexks9bT54FvAG3A0BNrHJhI89O1wKXB9jrgEeDLoz53IVDl7g8AuPvBCXyeBB7a1kYyUcklo5bsEBFJJ5ump5uBBcEyHrlS4+67g+09QFiv6lmkZoWvJ/WM7geBVe4+EHIuZrYCWAEwb54Wtw3j7mza3s6St8xiyqhVYEVE0slm1NOrQOd439jMHjSzZ0J+rh1+nqcesRf2gKQq4F3AF4G3A6cDn0z3ee6+1t3r3b1+9mxNIgvzfPtBdr3ezVKNdhKRccimRvEi8IiZ/RLoHSp09+9kepG7X5HumJm1mdlcd99tZnMJ73vYBTzt7i8Gr9lAahXb27OIuSxt2NKSdgnwDVta+No9zwDwvU07OG5SpVaBFZGsZJMoXgl+JgU/uXAvcANwa/D7npBzngSqzWy2u3cASwE9CJvwhACwen0j3X2plrmW/d18+a6tvPZGKrffdn8zvf2pLqY9XT2sXt8IoGQhImMyD3m0ZeiJZse5+6GcfKjZLOBOYB6wk9Tw2H1mVg98xt1vDM57L6mHJhmwGVjh7ofHev/6+npvaCjNnLJhS8uIhABQVWFUGBweyO7fckhddZLfrVqa6xBFpAiZ2WZ3D32ERDajni4h1dwzDZhnZucBn3b3zx5rQEHH+OUh5Q3AjcP2H2Bio6tKzpqNzSOSBED/4PgSxJDWkOdli4iMlk1n9neBZaSW7sDd/wS8O8KYJIPx3tzrqpPUVSdDj9WmKRcRGS6rtZ7c/dVRRaFDVCV66W7u1ckEyVFDXpOJSlYuW8DKZQvSHhMRGUs2ndmvmtk7ADezBHATsC3asCSdlcsWsGr9Vnr6Bo+UJROVfOOaRQBpRz2NdUxEJJ0xO7PN7ETgX4ErSHUq/xq4KccT8HKqlDuzAb79q+38+29eAFJNS7rpi8hETagz2933An+Z86jkmE1Ppv7Z/vS193H8cYmYoxGRUpfNqKfvhRR3Ag3uHjb/QSLWuKuTU2cdpyQhInmRTWf2FOB8YEfwcy5wMvApM/tuZJFJWlt3dXJO3fFxhyEiZSKbzuxzgSVDi/GZ2b8DjwHvBBojjE1C7HvjMC37u7nhHafGHYqIlIlsahQnkJpsN2QqMDNIHL3hL5GoNLak1mc8p6463kBEpGxkU6O4DXg6eByqkZps9y0zm0pq6W/Jo8Zd+wFYVDcj3kBEpGxkM+rpdjO7j9QT5gC+4u6twfbKyCKTUFt3dXL6iVOZMUUd2SKSH2mbnszsrcHvPwPmknouxavASUGZxKCxpZNzTlZHtojkT6Yaxd+SemLcv4Qcc1LLfksedRzoZXdnj0Y8iUhepU0U7r4i+H1Z/sKRTJ4JOrLPPbk63kBEpKyMOerJzI4zs783s7XB/plmdnX0ocloW3d1YgaLatWRLSL5k83w2P8ADgPvCPZbgH+KLCJJq7FlP2fMnsbUydkMVhMRyY1sEsUZ7n4b0AcQPOXOIo1KQm3d1cm56p8QkTzLJlEcNrMkqQ5szOwMNNEu79q6emg/0KsRTyKSd9m0YXwduB84xczuAJYAn4wyKDla466hjmwlChHJr2wm3D1gZk8BF5NqcropWHpc8mhrSycVBgvnKlGISH5l1SsaPKTolxHHIhk07trPmXOmk5xUOfbJIiI5lNUzsyVe7q4Z2SISGyWKIrC7s4e9Bw+rf0JEYpFpracpZnazmf2bmX3azHI6eN/MZprZA2a2I/h9QprzbjOzJjPbZmbfM7OyG5r75tLiShQikn+ZahTrgHpSDyd6P+FrPk3EKmCTu58JbAr2RzCzd5AaZXUu8Dbg7cB7chxHwWvc1UlVhXH2XM3IFpH8y1RLWOju5wCY2e3AH3P82dcClwbb64BHgC+POsdJPYp1EqkRVwmgLcdxFLytLZ2cWTOdKQl1ZItI/mWqUfQNbbh7fwSfXePuu4PtPUDN6BPc/Q/Aw8Du4Geju28LezMzW2FmDWbW0NHREUG48XB3Gnft14xsEYlNphrFeWbWFWwbkAz2DXB3H7MdxMweBE4KOfTV4Tvu7mbmIa9/C3A2cHJQ9ICZvcvdHxt9rruvBdYC1NfXH/VehWDDlhbWbGymdX83tdVJVi5bwPLFdRnLb/nVNl4/1Mf9TXu45IxZLF9cF/dliEiZybTM+ITbOdz9inTHzKzNzOa6+24zmwu0h5z2QeBxdz8YvOZXwCXAUYmi0G3Y0sLq9Y109w0A0LK/m9XrG2nYuY+7NreMWd7Z3cfq9Y0AShYikleZRj0tHbZ92qhj1+Xgs+8Fbgi2bwDuCTnnFeA9ZlZlZglSHdmhTU+Fbs3G5iM3/SHdfQPc8cQr4ypfs7E58lhFRIbL1Efxv4dt3zXq2N/n4LNvBd5rZjuAK4J9zKzezH4QnPMz4AVSI6/+BPzJ3X+eg8/Ou9b93aHlnqaRLF15uvcREYlKpj4KS7Mdtj9uwbIgl4eUNwA3BtsDwKcn+lmFoLY6SUvITb7CYDAkKaQrr61ORhCdiEh6mWoUnmY7bF/GsHLZAhKVI/NrMlHJ/7xoHslRw14zla9ctiDyWEVEhstUozjdzO4lVXsY2ibYPy39yyTM8sV1PNLczoanWzEYMbqp/tSZoaOe0pWLiOSTeZrGcDPLOAPa3X8TSUQ5UF9f7w0NDXGHcZRv3beNdb9/me3fvJIyXIlERAqYmW129/qwY5mGxx5JBGY2OygrnZlsMdjT2UPNjClKEiJSVDINjzUz+7qZ7QWagefMrMPMvpa/8EpLW1cPJ82YEncYIiLjkqkz+wvAO4G3u/tMdz8BuAhYYmZfyEt0Jaatq4c5MybHHYaIyLhkShQfB65395eGCtz9ReBjwCeiDqzUuDttXb2qUYhI0cmUKBJhz8YO+ikS0YVUmrp6+unuG6BGiUJEikymRHH4GI9JiLauHgBqjleiEJHiks3qsUNDdIbG0RqpZ0TIOAwlCjU9iUixiXT1WHnTns6gRqHObBEpMmkThZlNAT4DvAXYCvwwogcYlYX2A70A6qMQkaKT7TOzryL3z8wuK3s6ezg+mdDjTEWk6MT5zOyyosl2IlKs4nxmdlnRZDsRKVaRPjNb3rSnq4ezaqbHHYaIyLhp1FMeDAw6HQd6OUlzKESkCGVqepIc2Xuwl0GHOeqjEJEipESRB5psJyLFTIkiDzTZTkSKmRJFHrQFk+1UoxCRYqREkQdtnT1UVhizpqlGISLFR4kiD9q6epg9bTKVFXoEqogUHyWKPNjT1aP+CREpWrEkCjP7iJk1mdmgmdVnOO9KM2s2s+fNbFU+Y8yltq4eLQYoIkUrrhrFM8B1wKPpTjCzSuD7wPuBhcD1ZrYwP+HlVluXJtuJSPHKtIRHZNx9G4BZxjb7C4Hng+d0Y2Y/Aa4Fno08wBzq6Rugs7tPNQoRKVqF3EdRB7w6bH9XUBbKzFaYWYOZNXR0dEQeXLaOPAJViUJEilRkNQozexA4KeTQV939nlx/nruvBdYC1NfX+xin540m24lIsYssUbj7FRN8ixbglGH7JwdlRUWT7USk2BVy09OTwJlmdpqZTQL+Arg35pjGrS2oUWhBQBEpVnENj/2gme0CLgF+aWYbg/JaM7sPjjws6fPARmAbcKe7N8UR70S0dfWQTFQyY0os4wZERCYsrlFPdwN3h5S3kno+99D+fcB9eQwt54Ym240xwktEpGAVctNTSdBkOxEpdkoUEdNkOxEpdkoUEXL3oOlJiUJEipcSRYQ6u/s43D+oRCEiRU2JIkJ7ujTZTkSKnxJFhNq6NNlORIqfEkWE2jq1zpOIFD8liggNNT3NUdOTiBQxJYoItXX1MHPqJCZXVcYdiojIMVOiiFBbVw9zpqs2ISLFTYkiQppsJyKlQIkiQnu6eqiZrkQhIsVNiSIi/QOD7D3YS41qFCJS5JQoItJxsBd3TbYTkeKnRBERTbYTkVKhRBGRPZpsJyIlQokiIm1dShQiUhqUKCLS1tVDVYUxa+qkuEMREZkQJYqI7Akm21VU6BGoIlLclCgi0t6lobEiUhqUKCKiyXYiUiqUKCLS1tWj5TtEpCRUxR1AoduwpYU1G5tp3d9NbXWSlcsWsHxxXcbyb9+/nQM9/ax/ahfnn1LN8sV1cV+GiMgxM3ePO4acq6+v94aGhgm/z4YtLaxe30h338CRsmSikg9dUMddm1uyLr/lunOULESkoJnZZnevDzsWS9OTmX3EzJrMbNDMQgMzs1PM7GEzezY496Z8x7lmY/OImz5Ad98Adzz+yrjK12xsjjxWEZGoxNVH8QxwHfBohnP6gb9z94XAxcDnzGxhPoIb0rq/O7Q8XR0sXXm69xERKQaxJAp33+buGb9mu/tud38q2D4AbAPy2n5TW50MLa9MMzUiXXm69xERKQZFMerJzOYDi4EnMpyzwswazKyho6MjJ5/7xfedxeh7fzJRyfUXzSOZqMy6fOWyBTmJR0QkDpGNejKzB4GTQg591d3vGcf7TAPuAm52965057n7WmAtpDqzxxluqOqpk3DghOMS7D/UN2J0U/2pM0NHPaUrFxEpVpElCne/YqLvYWYJUkniDndfP/Goxmftb17kpBlTePRLlzGpamTla/niutAEkK5cRKRYFWzTk5kZcDuwzd2/k+/Pb9zVyR9efI2/WjL/qCQhIlJO4hoe+0Ez2wVcAvzSzDYG5bVmdl9w2hLg48BSM3s6+LkqXzGufexFpk2u4vqL5uXrI0VEClIsM7Pd/W7g7pDyVuCqYPu3cFRfcl68uu8Q9zXu5lPvPI0ZUxJxhCAiUjDUphLi9t++hAF/tWR+3KGIiMROiWKU/YcOc2fDq1xzXi1zj9f8BxERJYpR7njiFQ4dHuBv3n163KGIiBQErR4b2LClhdvu305rZw+Tqypo3nOAs+fOiDssEZHYKVFw9Cqxvf2DrF7fCKA5ESJS9tT0RPpVYrXqq4iIEgWQfnVXrfoqIqJEAaRf3VWrvoqIKFEAsHLZAq36KiKShjqzebPDWqu+iogcTYkioFVfRUTCqelJREQyUqIQEZGMlChERCQjJQoREclIiUJERDIyd487hpwzsw5g5xinnQjszUM4hUbXXV503eVlItd9qrvPDjtQkokiG2bW4O71cceRb7ru8qLrLi9RXbeankREJCMlChERyaicE8XauAOIia67vOi6y0sk1122fRQiIpKdcq5RiIhIFpQoREQko7JLFGZ2pZk1m9nzZrYq7niiZGY/NLN2M3tmWNlMM3vAzHYEv0+IM8ZcM7NTzOxhM3vWzJrM7KagvNSve4qZ/dHM/hRc9z8E5aeZ2RPB3/tPzWxS3LFGwcwqzWyLmf0i2C+X637ZzBrN7GkzawjKcv63XlaJwswqge8D7wcWAteb2cJ4o4rUj4ArR5WtAja5+5nApmC/lPQDf+fuC4GLgc8F/8alft29wFJ3Pw84H7jSzC4Gvg38H3d/C/A68Kn4QozUTcC2Yfvlct0Al7n7+cPmT+T8b72sEgVwIfC8u7/o7oeBnwDXxhxTZNz9UWDfqOJrgXXB9jpgeT5jipq773b3p4LtA6RuHnWU/nW7ux8MdhPBjwNLgZ8F5SV33QBmdjLw58APgn2jDK47g5z/rZdboqgDXh22vysoKyc17r472N4D1MQZTJTMbD6wGHiCMrjuoPnlaaAdeAB4Adjv7v3BKaX69/5d4EvAYLA/i/K4bkh9Gfi1mW02sxVBWc7/1vWEuzLm7m5mJTk+2symAXcBN7t7V+pLZkqpXre7DwDnm1k1cDfw1ngjip6ZXQ20u/tmM7s05nDi8E53bzGzOcADZrZ9+MFc/a2XW42iBThl2P7JQVk5aTOzuQDB7/aY48k5M0uQShJ3uPv6oLjkr3uIu+8HHgYuAarNbOgLYSn+vS8BrjGzl0k1JS8F/pXSv24A3L0l+N1O6svBhUTwt15uieJJ4MxgRMQk4C+Ae2OOKd/uBW4Itm8A7okxlpwL2qdvB7a5+3eGHSr1654d1CQwsyTwXlL9Mw8DHw5OK7nrdvfV7n6yu88n9f/5IXf/S0r8ugHMbKqZTR/aBt4HPEMEf+tlNzPbzK4i1aZZCfzQ3f853oiiY2b/DVxKaunhNuDrwAbgTmAeqaXYP+ruozu8i5aZvRN4DGjkzTbrr5Dqpyjl6z6XVMdlJakvgHe6+z+a2emkvmnPBLYAH3P33vgijU7Q9PRFd7+6HK47uMa7g90q4L/c/Z/NbBY5/lsvu0QhIiLjU25NTyIiMk5KFCIikpEShYiIZKREISIiGSlRiIhIRkoUInlkZr+POwaR8dLwWBERyUg1CpE8MrODY58lUliUKEREJCMlChERyUiJQkREMlKiEBGRjJQoREQkIw2PFRGRjFSjEBGRjJQoREQkIyUKERHJSIlCREQyUqIQEZGMlChERCQjJQoREcno/wOZ4imkU5kxXAAAAABJRU5ErkJggg==)
%% Cell type:markdown id: tags:
## Step 3: GW
%% Cell type:markdown id: tags:
The GW electronic structure is computed treating the frequency integration of the correlation part of the self energy with the Contour Deformation techinique and by computing the dielectric screening at multipole frequencies with Lanczos iterations. Check out the ``wfreq.x`` [input description](http://www.west-code.org/doc/West/latest/) and generate an input file for WEST called ``wfreq.in``.
%% Cell type:markdown id: tags:
Download this file in your current working directory:
%% Cell type:code id: tags:
``` python
%%bash
wget -N -q http://www.west-code.org/doc/training/silane/wfreq.in
```
%% Cell type:markdown id: tags:
Let's inspect the ``wfreq.in`` file, input for ``wfreq.x``.
%% Cell type:code id: tags:
``` python
%%bash
cat wfreq.in
```
%%%% Output: stream
input_west:
qe_prefix: silane
west_prefix: silane
outdir: ./
wstat_control:
wstat_calculation: S
n_pdep_eigen: 50
wfreq_control:
wfreq_calculation: XWGQ
n_pdep_eigen_to_use: 50
qp_bandrange: [1,5]
n_refreq: 300
ecut_refreq: 2.0
%% Cell type:markdown id: tags:
Run ``wfreq.in`` on 2 cores.
%% Cell type:code id: tags:
``` python
%%bash
mpirun -n 2 wfreq.x -i wfreq.in > wfreq.out
```
%% Cell type:markdown id: tags:
The output file ``wfreq.out`` contains information about the calculation of the GW self-energy, and the corrected electronic structure can be found in the file ``west_prefix.wstat.save/wfreq.json``.
The output file ``wfreq.out`` contains information about the calculation of the GW self-energy, and the corrected electronic structure can be found in the file ``<west_prefix>.wstat.save/wfreq.json``.
%% Cell type:code id: tags:
``` python
import json
import numpy as np
with open('silane.wfreq.save/wfreq.json') as json_file:
data = json.load(json_file)
k=1
kindex = f"K{k:06d}"
print(kindex)
print()
bandmap = data["output"]["Q"]["bandmap"]
eqp = data["output"]["Q"][kindex]
for i, b in enumerate(bandmap) :
print(f"band = {b}")
for key in eqp:
if("re" in eqp[key]):
print(f" {key} = {eqp[key]['re'][i]:.3f}, {eqp[key]['im'][i]:.3f} eV")
else:
print(f" {key} = {eqp[key][i]:.3f} eV")
print()
```
%%%% Output: stream
K000001
band = 1
sigmax = -17.606 eV
vxcl = -11.250 eV
vxcnl = 0.000 eV
hf = -6.357 eV
z = 0.661 eV
eks = -13.236 eV
eqpLin = -16.273 eV
eqpSec = -16.135 eV
sigmac_eks = 1.765, 0.015 eV
sigmac_eqpLin = 3.577, 0.076 eV
sigmac_eqpSec = 3.439, 0.116 eV
sigma_diff = -0.023 eV
band = 2
sigmax = -15.766 eV
vxcl = -11.243 eV
vxcnl = 0.000 eV
hf = -4.523 eV
z = 0.869 eV
eks = -8.231 eV
eqpLin = -12.151 eV
eqpSec = -12.046 eV
sigmac_eks = 0.011, 0.003 eV
sigmac_eqpLin = 0.730, 0.009 eV
sigmac_eqpSec = 0.707, 0.009 eV
sigma_diff = -0.003 eV
band = 3
sigmax = -15.765 eV
vxcl = -11.243 eV
vxcnl = 0.000 eV
hf = -4.523 eV
z = 0.869 eV
eks = -8.231 eV
eqpLin = -12.148 eV
eqpSec = -12.044 eV
sigmac_eks = 0.013, 0.003 eV
sigmac_eqpLin = 0.732, 0.009 eV
sigmac_eqpSec = 0.709, 0.009 eV
sigma_diff = -0.003 eV
band = 4
sigmax = -15.765 eV
vxcl = -11.242 eV
vxcnl = 0.000 eV
hf = -4.523 eV
z = 0.869 eV
eks = -8.230 eV
eqpLin = -12.148 eV
eqpSec = -12.043 eV
sigmac_eks = 0.013, 0.003 eV
sigmac_eqpLin = 0.732, 0.009 eV
sigmac_eqpSec = 0.709, 0.009 eV
sigma_diff = -0.003 eV
band = 5
sigmax = -0.587 eV
vxcl = -2.090 eV
vxcnl = 0.000 eV
hf = 1.503 eV
z = 0.980 eV
eks = -0.466 eV
eqpLin = 0.666 eV
eqpSec = 0.665 eV
sigmac_eks = -0.349, -0.000 eV
sigmac_eqpLin = -0.373, -0.000 eV
sigmac_eqpSec = -0.373, -0.000 eV
sigma_diff = 0.000 eV
%% Cell type:markdown id: tags:
Explanation of keys:
- `k` : k-point
- `s` : spin
- `b` : band
- `sigmax` : $\langle k,s,b|\Sigma_x|k,s,b\rangle$, Exchange self-energy
- `vxcl` : $\langle k,s,b|V_{xc}^{local}|k,s,b\rangle$, Local exchange and correlation
- `vxcnl` : $\langle k,s,b|V_{xc}^{non-local}|k,s,b\rangle$, Non-local exchange and correlation
- `hf` : $\varepsilon_{k,s,b}+\langle k,s,b|\Sigma_x-V_{xc}^{local}-V_{xc}^{non-local}|k,s,b\rangle$, perturbative Hartree-Fock
- `z` :
- `eks` : $\varepsilon_{k,s,b} = \langle k,s,b|H_{KS}|k,s,b\rangle$, Kohn-Sham eigenvalue
- `eqpLin` : $\varepsilon_{k,s,b} + \langle k,s,b|\Sigma_c(\varepsilon_{k,s,b})|k,s,b\rangle$, Kohn-Sham eigenvalue
- `eqpSec` : $\varepsilon_{k,s,b} + \langle k,s,b|\Sigma_c(\varepsilon_{k,s,b})|k,s,b\rangle$, Kohn-Sham eigenvalue
- `sigmac_eks` : $\langle k,s,b|\Sigma_c(\varepsilon_{k,s,b})|k,s,b\rangle$, Correlation self-energy evaluated at $\varepsilon_{k,s,b}$
- `sigmac_eqpLin` : $\langle k,s,b|\Sigma_c(\varepsilon_{k,s,b})|k,s,b\rangle$, Correlation self-energy evaluated at $\varepsilon_{k,s,b}$
- `sigmac_eqpSec` : $\langle k,s,b|\Sigma_c(\varepsilon_{k,s,b})|k,s,b\rangle$, Correlation self-energy evaluated at $\varepsilon_{k,s,b}$
%% Cell type:code id: tags:
``` python
from westpy import ElectronicStructure
```
%%%% Output: stream
_ _ _____ _____ _____
| | | | ___/ ___|_ _|
| | | | |__ \ `--. | |_ __ _ _
| |/\| | __| `--. \ | | '_ \| | | |
\ /\ / |___/\__/ / | | |_) | |_| |
\/ \/\____/\____/ \_/ .__/ \__, |
| | __/ |
|_| |___/
WEST version : 3.1.1
Today : 2021-05-19 21:55:09.144532
%% Cell type:code id: tags:
``` python
es = ElectronicStructure()
es.addKey("eks","Kohn-Sham energy (eV)")
es.addKey("eqp","Quasiparticle energy (eV)")
```
%% Cell type:code id: tags:
``` python
for i, b in enumerate(bandmap) :
es.addDataPoint([1,1,b],"eks",eqp["eks"][i])
es.addDataPoint([1,1,b],"eqp",eqp["eqpSec"][i])
```
%% Cell type:code id: tags:
``` python
es.plotDOS(k=1,s=1,energyKeys=["eks","eqp"],energyRange=[-20,5,0.01])
```
%%%% Output: stream
Requested (emin,emax) : -20 5
Detected (emin,emax) : -16.134751796731177 0.6648562238096314
output written in : dos.png
waiting for user to close image preview...
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0d0lEQVR4nO2deZhU9Znvv28tvbI12CDYEnQe4aoQ2WQJA8IAihKXMFE0MZHLxOUJGHVmYswTb1ySiUkmZq5xHa+iSSTG4IbBUVxxRWSxjQtgAEVAlO4GBHqpruW9f5ylTp2q6q6qU3XW9/M8/VTVOadO/ar6d8739y6/90fMDEEQBCG4hJxugCAIguAsIgSCIAgBR4RAEAQh4IgQCIIgBBwRAkEQhIATcboBhXDUUUfxiBEjnG6G0BPtLcCXu4EhJwPhKqdbY4mWwzF8fqgLIwbVo2+NJy4RQcjJxo0bW5m5sbfjPNHLR4wYgQ0bNjjdDKEn7j8D2HUQWHQ3MGKa062xxL+veBePbtyNm/95DBaeOtzp5ghCyRDRzkKOE9eQUB7CUeUxdsjZdpSBSIgAADLFRggKIgRCeUl0Od2CspESIRACggiBUB5SCeUxGXe2HWWASLEI4smUwy0RBHvwRIxA8ACppPKYiDnbjjKg6gBiiaSzDSkD8Xgcu3fvRleXfyw1IZuamho0NTUhGo2W9H4RAqE86BaBD4RAfYzFvW8R7N69G3379sWIESN0S0fwF8yMtrY27N69G8cdd1xJ5xDXkFAeWL1pJrqdbUcZSKlR4ljC+0LQ1dWFQYMGiQj4GCLCoEGDLFl9IgRCeQiFlcek94WgO6EJgfddQwBEBAKA1f+xCIFQHkKql9EPQqAGif1gEQhCIYgQCOXFB8HibtUS6BYhEAKCCIFQHrSsIR9YBPGk4hpKykSCitKnT5+S3xuLxfCb3/wGkyZNwtixY3HOOefgjTfe6PE9V199NV599VUASrWC1tbWoj93xYoVOPnkkxEKhcpe7eDVV1/F+PHjEYlE8Oijj+rbW1paMG/evLJ+lhkRAqE8sH+EQLMERAjcSSwWw1lnnYVYLIbnn38ezc3NuPXWW3HTTTfh8ccfz/metrY2vPXWW5gxY4alzx49ejQef/xxy+fJxfDhw/Hggw/iW9/6Vsb2xsZGDB06tFehs4KkjwrlQbcIvD+hTBOChM+E4Ka/foAPPytvCZCThvXDDWef3OMxDz30EH73u9+hu7sbkydPxl133YVwOKzvb21txdlnn43rr78e48ePx8KFC3Ho0CEkEgncfffdmD59esb5brnlFpx//vm44oor9G0nnHACVq5ciTlz5uDMM89EbW1txnsee+yxnKPqzs5OLFiwAAsWLMCll17a6/c98cQTez1m7969vX6HXGiFNUOh7PH5eeedh+XLl2PatMrU8RKLQCgPWvooez/TJpYUi6BcbN68GY888gjeeOMNNDc3IxwOY/ny5fr+L774AvPnz8fNN9+M+fPn409/+hPOOOMMNDc3491338XYsWOzzvk///M/uPzyy7Ft2zZMnz4dp512Gn7wgx/gnXfewfnnn49nnnkm6z1vvPEGJkyYkLHtyJEjOPvss3HRRRfpIjB9+nSMHTs26++FF14o+Dvn+w4LFy7Mee4//OEPvZ5z4sSJeO211wpuQ7GIRSCUB80i0CaWeZi4bhH4K1jc28i9Erz44ovYuHEjTj31VADKCHzw4MEAlFnPs2fPxp133onTTjsNAHDqqadi8eLFiMfjOO+887KEoKWlBcceeyyICNdddx1uu+02nHjiiZg5cyYWLFiAUaNG4f33389qx969e9HYmFmN+dxzz8W1116Lb3/72/q2ctxs832HRx55pORzDh48GJ999pnltuVDLAKhPGiWQMr7FkG3WARlg5lxySWXoLm5Gc3Nzdi6dStuvPFGAEAkEsGECROwevVq/fgZM2bg1VdfxTHHHINFixblHC1rbqW2tjaMHz8etbW1mDlzJgBg3759utAYqa2tzZpwNW3aNDz77LNgQ5nZclgE+b6DFYugq6sry91VTkQIhPKQ8o8QpFQB8FuMwAlmz56NRx99FPv27QMA7N+/Hzt3KiXyiQjLli3Dli1b8Ktf/QoAsHPnTgwZMgSXXnopvve972HTpk0Z52tsbMSuXbvAzGhoaEBzczO6urrwyiuv4ODBg/j973+Pr3/961ntOPHEE7Ft27aMbTfffDMaGhqwZMkSfdtrr72mi5bxb86cOT1+zz179mD27Nk9fodHHnkk57m/+93v9vo7fvTRRxg9enSvx5WKCIFQHtg/rqEkS/pouTjppJPw85//HKeffjq++tWvYu7cudi7d6++PxwO4+GHH8ZLL72Eu+66C2vWrMEpp5yCcePG4ZFHHsFVV12Vdc5Zs2bhgQcewC233IKlS5di3rx5mDp1Ku655x78+te/xqBBg7LeM3/+fKxZsyZr+2233YbOzk5ce+21BX2fJ554Ak1NTVi7di3mz5+PM844A4DieopEFE97Id8hF+vXr0dTUxNWrFiByy+/HCefnHblvfzyy5g/f35B5ykJZnb934QJE1hwOb8dzXxDP+Y/X+x0Sywz7Zcv8ld+tIovunet002xzIcffuh0E8pOe3s7T58+ne+66y7u7OxkZuadO3fyfffd1+P7pk2bxgcOHKhIm26//XZeuXJlRc7NzDx9+nTev39/j8fk+l8D2MAF3GPFIhDKg2YRsPcDrOIacjd1dXVYvXo12traMGPGDIwZMwZLlizByJEje3zfrbfeik8//bQibVq6dCnOOeecipy7paUF//qv/4qGhoaKnB+QrCGhXPgoa0hcQ+6ntrYW119/Pa6//vqC3zN58uQKtqhyNDY24rzzzqvoZ4hFIJQHP8UItIraIgRCQBAhEMqDn7KGdIvA+24uQSgEEQKhPPjKIlBjBEmxCIRgIEIglAdt9OwHiyAlMQKvsH37dixevBijR4/GhAkTcM011+DAgQN5j+/s7MRpp52GZDKJNWvW5JxzUAiLFy/G4MGDK5Lb/5Of/ATHHntsVnXWO+64A8uWLSv75wEiBEK50LOGfCAEEiz2BOvWrcMFF1yAhQsX4t1338X69esxbdo0zJs3D21tbTnfs2zZMixYsCCj6F0pLFq0CM8++6ylc+Tj7LPPxttvv521ffHixbj99tsr8pkiBEJ50FxCfnANsaSPlpOHHnpIXzfg8ssvRzKpDBYeeOABjBw5EpMmTcKll16KpUuXAlBusldccQUmTpyIkSNHYtWqVVnnTCaTuPLKK/HXv/4VZ5xxBsLhMEKhEL75zW/iF7/4BX7605/mbMvy5ctx7rnnZm1fv349xo0bh+3btxf0nWbMmIGBAwf2eMyKFSswevRonHLKKUWVrZ4yZQqGDh2atb2urg4jRozIKRJWkfRRoTz4KViserl8ZxE8cx3w+XvlPefRY4Azf5l3t7H6aDQaxfe//30sX74cc+fOxQ033ICNGzeif//+mDVrFsaNG6e/75NPPsHbb7+N7du3Y9asWdi2bRtqamr0/S+++CLmzp2LYcOG4b777sOdd96J8ePHIxaL4aGHHsJNN92U1Zbu7m7s2LFDL/es8eabb+LKK6/EypUrMXz4cLz88su45pprst5fV1eHN998s+Cf5uabb8bq1atxzDHH4ODBgwCArVu3YuHChTmPX7NmDQYMGNDjObUqpJMmTSq4HYUgQiCUBx8VnZN5BOUjX/XRdevWYebMmXpF0IULF+Kjjz7S33fBBRcgFArhhBNOwPHHH48tW7ZkVCJ99913MWXKFLS0tOCPf/wj1q5di/feew8XXnghAGDo0KFoaWnJqDja2tqadaPdvHkzLrvsMjz33HMYNmwYAKWERXNzs+XvPm3aNCxatAgXXHABFixYAAAYNWqUpXMPHjwYW7Zssdw2MyIEgnWMaZZ+cA35dWZxDyP3SsFq9dFbbrklY/uTTz7Z4/uIqMfXgFKnaMeOHZg6dSpqampw6qmn4qijjgIAHDhwIGsmbq4KpEOHDkVXVxfeeecdXQjKZRHcc889WLduHZ5++mlMmDABGzduRGtrqyWLoFJVSEUIBOsYA8QeDxanDDd/mUdgndmzZ+Pcc8/FNddcg8GDB2P//v04fPgwJk+ejKuuugptbW3o168fVqxYgVNOOUV/34oVK3DJJZfg448/xo4dOzBq1KiM844ePRpvvfUWli5dirVr1yIWi+GDDz5Aa2srXnrpJQwbNkwvAqfR0NCAZDKJrq4u3c00YMAA3H///Zg7dy7q6+sxc+ZMSxbBHXfcAUApObF9+3ZMnjwZkydPxjPPPINdu3Zh7NixliyCjz76qCKrlEmwWLCO0R3kcYtAcwsBPrQIHCBf9dGhQ4fixhtvxNSpUzFt2rSsJSCHDx+OSZMm4cwzz8Q999yTER8AgDlz5uDpp59GKpXCt771LUyZMgV33nknxowZg8ceeyxvds3pp5+O119/PWPbkCFDsGrVKixZsgTr1q0r6HtddNFFmDp1KrZu3Yqmpibcf//9AIAtW7bo1U9/+MMfYsyYMRg9ejS+9rWvZQhdT1x77bVoampCR0cHmpqa9PUbAGWltblz5xZ0nqIopDJdKX8AlgHYB+B9w7aBAJ4H8Hf1saGQc0n1UZcTO6JUHr2hH/N/jXa6NZbo7E7wV360ir/yo1V84v95xunmWMYr1UcfeOABXrJkCTMzX3LJJbxixYpe3/PKK6/wpEmT+K233mJm5kQiwWvWrOE1a9bkfc/GjRv54osrVyF3/vz5HIvFKnLuTZs29dh2t1YffRCAebXo6wC8yMwnAHhRfS14nQyLwNvuFG0OAZFYBG5nxowZePDBB/G73/0OY8eOxfjx4/HEE09k1PE3M378eMyaNUtPYS03q1atQlVVVUXO3draip/97GcVOTcxV66zE9EIAKuYebT6eiuAmcy8l4iGAljDzKN6OgcATJw4kTds2FCxdgoW6TwA/GqE8rzP0cC/b3W0OVY4Ektg9A2rURUJIZlibP/FWU43yRKbN2/OcrsI/iTX/5qINjLzxN7ea3eMYAgza8sTfQ5gSL4DiegyItpARBtaWlrsaZ1QGj7KGtIyhqrCihBUcqBkF374DkLPWP0fOxYsVv1XeVvPzPcy80RmnmjMBRZciJYpRCHfZA1VRZRLw+veoZqaGrS1tYkY+BhmRltbW1ZAvRjsTh/9goiGGlxD+2z+fKESaDGCcLXnJ5RpWUPRsJK3nkwxwqHsHHav0NTUhN27d0Osan9TU1ODpqamkt9vtxA8BeASAL9UH1fa/PlCJdCsgHAVkIo72xaLaBZBNKxYBF6fXRyNRnHcccc53QzB5VTMNUREDwNYC2AUEe0mon+BIgBziejvAOaorwWvo1kBkSrfWASaaygpLhUhAFTMImDmi/Lsml2pzxQcgo2uIf8Ei42vBcHPyMxiwTpa1lA46oNgsfKouYZSIgRCABAhEKxjjBEAnp5Upk0oi2jBYnENCQFAhECwTsokBB62CtJZQ2IRCMFBhECwjm4RRJVHDweMU6YYgZSZEIKACIFgHR9aBHrWkAiBEABECATrpPxjESRTmRPKUhIjEAKACIFgHXOw2MMWgTlrSCwCIQiIEAjW0SeUVauvvZs1pLuGtGCxWARCABAhEKxjDhZ72CJIZpWYcLI1gmAPIgSCdczBYg/HCDQLIBpJF50TBL8jQiBYx1cxAn8VnROEQhAhEKxjLDEBeNoiMMcIZGaxEARECATr+MoiUB7FIhCChAiBYJ2sGIF3I6xZJSbEIhACgAiBYB0fZQ3pMQIJFgsBQoRAsI6PsobM6xFI0TkhCIgQCNbxkUVgrjUkReeEICBCIFhHzxryvkWQlT4qMQIhAIgQCNbxUdaQduOPhNSic2IRCAFAhECwTlb1Ue9mDWn3fSlDLQQJEQLBOj6yCMyuIUkfFYKACIFgHR9mDUnROSFIiBAI1vFh1pC2ME3Cw24uQSgUEQLBOj7MGtJiBOIaEoKACIFgHR/FCLKKzolBIAQAEQLBOn7KGjIHiyVrSAgAIgSCdVIJ5dEPFoFMKBMCiCNCQETXENEHRPQ+ET1MRDVOtEMoE7prSFuz2MNCoN73tWCxzCMQgoDtQkBExwD4AYCJzDwaQBjAhXa3Qygj5oVpPGwRMMsKZULwcMo1FAFQS0QRAHUAPnOoHUI5MAeLvWwRmLKGRAiEIGC7EDDzHgC/AfApgL0AvmTm58zHEdFlRLSBiDa0tLTY3UyhGLQbfyiiPLJ3g8WyMI0QRJxwDTUAOBfAcQCGAagnoovNxzHzvcw8kZknNjY22t1MoRg4CVAYCIWV1x62CFIpc/qoCIHgf5xwDc0B8DEztzBzHMDjAL7mQDuEcpFKKiJAanfycIxAmzegr1AmFoEQAJwQgk8BTCGiOiIiALMBbHagHUK58JFFkC5DLfMIhODgRIxgHYBHAWwC8J7ahnvtbodQRlIp1SJQhcDDFkEqxQhRej0CmVksBIGIEx/KzDcAuMGJzxYqgM8sgnCIENKFQJRA8D8ys1iwTioJhEIGi8C7N88UM0KkiEA4RBIjEAKBCIFgHR9ZBIprSBUCInENCYFAhECwjs+yhsKqWygUknkEQjAQIRCs4yeLgJVgMaBZBCIEgv8RIRCs46OsoWSKdYsgHBIhEIKBCIFgHU4qbiEfWARa1hAgQiAEBxECwTp6jMAHWUMpyRoSgocIgWAdH8UIjK6hEJHMLBYCgQiBYB3dIiAA5O0YgXkegQiBEABECATrcCrtFgqFPW0RMCPDIhDXkBAERAgE66QSabcQhb1tEaQM6aMhcQ0JwUCEQLCO5hoCPG8RJJn1OkORECEhQiAEABECwTpasBhQLQJvZw2FSZtZTDKzWAgEIgSCdTIsgpC3LQLjhDKZWSwEBBECwTrGYLHHYwTG6qOhkBSdE4KBCIFgHT/FCDJKTEjROSEYiBAI1tFKTACetwiSDD1YLK4hISiIEAjWybIIvOtPUYLFyvOwBIuFgCBCIFgnK2vIuxaBeYWyRFKEQPA/IgSCdXyWNRSSmcVCwBAhEKzjs6yhsMEikJnFQhAQIRCsoy1eD/gsa0gsAiEYiBAI1vFRjMCYNSRlqIWgUJQQEFGUiMYR0eBKNUjwID7NGoqIRSAEhB6FgIjuIaKT1ef9AbwL4A8A3iGii2xon+AFMiyCkLctAuPCNJI1JASE3iyC6cz8gfr8fwP4iJnHAJgA4NqKtkzwDtri9YDnYwQZ6aMk8wiEYNCbEHQbns8F8CQAMPPnlWqQ4EH8FCMwB4slRiAEgN6E4CARfZ2IxgGYBuBZACCiCIDaSjdO8Ag+yhoyF50THRCCQG9CcDmApQAeBHC1wRKYDeDpUj+UiAYQ0aNEtIWINhPR1FLPJbgAP61HkFFrCGIRCIEg0tNOZv4IwLwc21cDWG3hc28D8Cwzf5OIqgDUWTiX4DR+qz6qZg2FxDUkBIRe00eJ6EwieoWIWtW/V4jorFI/UM0+mgHgfgBg5m5mPljq+QQXkPJX1pBxqUoRAiEI9GgRENGlUNxD1wLYoG6eCOCXRNTEzPeW8JnHAWgB8AARnQJgI4CrmLnd9NmXAbgMAIYPH17Cxwi2wUkgpHYlj1sE5hITMo9ACAK9WQTXADidmV9i5kPq30sAzlT3lUIEwHgAdzPzOADtAK4zH8TM9zLzRGae2NjYWOJHCbZgDBb7KGtIZhYLQaE3ISBm3m/eyMxtFj5zN4DdzLxOff0oFGEQvIoxWOwDiyAktYaEgNGbEBxS3TcZqNsOl/KBaubRLiIapW6aDeDDUs4luARjsNjjWUNKsNhQhlosAiEA9BgjAPBvAJ4iogeg+PIBJUZwCYCLLXzulQCWqxlDO6DMWha8SoZF4O31CFKMjAll4hoSgkBv6aOvE9FkAN8HsEjd/CGAKVZmFzNzMxRBEbyOVmAu5I+ZxakUgwxF5xIiBEIA6M0i0Fw5PyWiRvV1S8VbJXgH7abvkxhB0pA1FJI1i4WA0Fv1USKiG4moBcBWAFuJqIWIfmpP8wTXo930fZg1FJYYgRAQCkkfnQZgEjMPZOaBACYDmEZEpaaPCn4ip0Xg3WCxMWtIqzXEYhUIPqc3IfgOgIuY+WNtAzPvgBIo/m4lGyZ4BN0i8EeMwJg1pD2KUSD4nd6EIMrMreaNapwgWpkmCZ4iyyLwbtYQM2cWnVOvDnEPCX6nmPUIitknBAUfZQ1p9/t0iYmQul2EQPA3vWUNnUJEh3JsJwA1FWiP4DV0i8D76xFoI3/NEtAeJYVU8Du9zSMI29UQwaP4KEagjfxDhlpDgLiGBP/TaxlqQegRH2UN6UJgqD4KQGYXC75HhECwRpZF4N31CLSRfySUKQRSeE7wOyIEgjV8NLNYE4IQZbqGxCIQ/I4IgWANH2UNpYPF6RXKALEIBP8jQiBYw09ZQ+ZgsfqYSIoQCP5GhECwRq6sITDgwVG0ZtxEQuaZxd77LoJQDCIEgjVSCeXRuGYx4EmrIKEqQdiUNSTpo4LfESEQrGEOFmsuIg/GCfRwh8k1JBaB4HdECARrmIPFHrYItBiBPrNYn1DmVIsEwR5ECARrZFkE4cztHsKcPipF54SgIEIgWMMcLPayRaBPKFMuC63onAiB4HdECARr6MFis0XgPX9KvqJzMo9A8DsiBII1NBeQD7KGzLWGpOicEBRECARraBaBD7KGzDOLw5I1JAQEEQLBGr7MGsqcUCYWgeB3RAgEa3CumcXwhUUQkjLUQkAQIRCsYXYNedki0ISAMovOyQplgt8RIRCskTIFiz2cNaSN/M0ziyVrSPA7IgSCNczpo162CNi0MI2sRyAEBMeEgIjCRPQOEa1yqg1CGdBG/j7IGkqYLAIpOicEBSctgqsAbHbw84Vy4COLIGWKEYSkDLUQEBwRAiJqAjAfwH1OfL5jpFLA9peAeKfTLcmiO5HCG9tawcXe9HKuRwBPWgT55hFI0Tl3EUsk8drfW0pz2aVSwLYXgXhX+RvmYZyyCP4vgGsB5L3EiOgyItpARBtaWlpsa1hF+fgV4I/fAN683emWZHHf6zvw7fvWYd3H+4t7ow9nFmcJgVgEruLBNz7Bd+5/G698VMJ94eM1wEMLXHkNOontQkBEXwewj5k39nQcM9/LzBOZeWJjY6NNraswX+5WHvdscrYdOdi08yAA4JPW9uLemDWz2LtZQ4m8FoH3vouf+XDvIQDA7gMdxb/5wE7l8Yv3ytgi7+OERTANwDlE9AmAPwP4JyJ6yIF22I/mEkq4zzWkZcoc6IgX90Z9ZrFmEahdyoMWQVYZalmPwJVoHqGi+yoAHPpMeQxXl69BPsB2IWDmHzNzEzOPAHAhgJeY+WK72+EImgC4MEYQSyg37oMd3cW9UQ8Wq13JwzECs2tI1zTJGnIVR7oUAdjfXmRfBYC4akV0fVnGFnkfmUdgJ3H3CsGXncrFdaBYITAvTOPhGEHSvHi9xAhciWYJFD1oAdLXXkdbGVvkfSJOfjgzrwGwxsk22Io2GnGxELTHiryB551Z7EUhUJQgJEXnXE1XXOlb7d0l9LGEmi0ULyG+4GPEIrATLWXNlUKguHg6uhPFvdFH8wg0iyBMmRaBzCNwF5oQdJYiBJoAdBeZFOFzRAjsRLcI3NcJNQHojBd5cbEpWOzhrCHNBRTSVyhTi84lRQjcRFdc6VtFD1oAw2BMLAIjIgR2opml3e7qhMxc+ihLTx9Vu5KHs4ayZhaLReBKutTEho6SXEOqNe6ya9BpRAjsRHMJJWOuulHGk6yn5BV9caWSigioN08vxwgS5sXrJUbgSmKqRdBVrPUKZFoEIvA6IgR2YjRHXeSj1EZYQCmuoWTaLQR4OkaQLkOtvJasIffBzOWxCMBpC10QIbAVY30TF/kotREWUKJrSLMCAE9bBOalKkNShtp1dCdT+kC+tGCxIVFD3EM6IgR24laLQLUC+tdGi7cIUinfWARZM4ul6Jzr0ALF1ZEQOuLJ4oskZgzG3HMNOo0IgZ0kugCovnQ3WQSqqd1QpwhBURdXKpH2pQCezhpK6TECzSJQtkutIfcQUwcqA+urkEwx4sVmdCU6gZr+ynOxCHRECOwk3gHUDVKfu2cugTbKaqivAnP6dUFwMtM15OGsIXPROSJCOEQSI3ARsYTSNwfUVQEowT0U7wJqB6rPRQg0RAjsJN4J1A1MP3cJaYtAubiKysZImYLFHo4RpJhBpAiARphIXEMuQuubA2qjyutEkf0sYbgGJVisI0JgJ/EuoLZBee6iTqhZACVdXKlEOi4AeD5GEDaIAKAYODKPwD3ofbVO7avFDFqScaW/ategiwZjTiNCYCfxDoNZ6p5OqF1M/TQhKMo1lPKNRZBk1t1CGopFIELgFrRBygDdei2ir2rXXK1YBGZECOwiGVdujm62CEoZZaUS6VnFgLctgmS2EIRCIgRuQncNldJXtWtOLIIsRAjsQgtMuThGoLuGio4R5JpH4D3HepKzXUMREQJXoQ1a+pfSV3WLwH2DMacRIbCLuGk04qJOmLYISjC3s2YWezdrKJVivb6QhmQNuYusQUuiiL4qFkFeRAjsQrMIXNgJjRPKgBKCxT6ZWZxIsT6HQCNEJDOLXYQlN6ZYBHkRIbALvRMOUB5d1Am1G39/7eIqJjfbRzOLE0lGJJzDIhAhcA3pQUsJqc5mIYi75xp0GhECu9CKXUXrgUiNqyyCmNnvWnT6aK6ZxR4UghTrlUc1QiSuITdhDhbHinFjatdgdV+lnybccw06jQiBXWg3/mitIgQuswiqIiHURpWbuLWZxZpF4L1gcSKVymkRiGvIPaRnFpcwaNEsgGiNch2KRaAjQmAXeiesVTuhe0YjsXgKNZEQanQhKNYi8Mc8gkQyO0YQCRFkgTL30BVPggjoW1NK+qh6zUW0wZh7rkGnESGwCy1Y7EaLIJ5EdTSMmmhIfV3EaN5cYsLDWUOJVArRsMk1FCIpOuciYokUaiJh1ERK6KtiEeRFhMAu4obRiNssgkQKNdEQaiIlWATJOBCOZG6jsGctAplZ7G6UQUsIkXAIkRCJRVAmRAjsIuHiGEE8iZpIGKEQoSocKjJYHAdC0cxtobAnLYJ4ihHJaRE41CAhC62vAkBNNFxaiYlorWIViEWgI0JgFxmd0F1maVc8qccHqqOh4jIxknEgbBICj1oEyVQK0awJZVJ0zk10xVO6C7MmWuSgxRini9RKGWoDIgR2kZU15B6zNJZIoTqiXVxha8FiQLUIvDeMjotryPVkDFoiRfbVRKcySAlHFYvARVa504gQ2IUeI3BfoMp4cdVEQyXECPxhESSSuYPFYhG4h1gihWpDXy3Keo13KdceoFoE7hmMOY0IgV3EO5TOR+RKi0Azt2uL9bvmjBGEPBkjSKSyZxZHQoSE5I+6BiVGUKL1muhUrj1ALAITtgsBER1LRC8T0YdE9AERXWV3GxwhYRiNuCxQ1RVPotoYgCvG75pMZFsEoagiEB4j1zwCmVnsLroSKYP1GkZnUSUmzBaBe65Bp3HCIkgA+DdmPgnAFABLiOgkB9phL/FOIFqnPI/Uus4iqNYCcMX6XVPx7BhBuEoRCI+RSKWySkzIzGJ3EYsnDfGsIt2Y8Q6TReCea9BpbBcCZt7LzJvU54cBbAZwjN3tsJ14p9L5ABdaBCndIqiOhopzDeWKEYSjQLK7jC20h7xF58QicA0Z8axIkW7MDKu8zlXXoNM4GiMgohEAxgFYl2PfZUS0gYg2tLS02N62shPvzDRLE52AS24wsUTSkJJXStaQWQiqvCkEUoba9WSmjxbpxsy4BsUiMOKYEBBRHwCPAbiamQ+Z9zPzvcw8kZknNjY22t/AcqMFi4G0ZZCIOdceAzGDRVATDeuFvQoi2Z09s9irQpBMZU0oE4vAXSiDFkNfLdYi0F1DtcogxoMuzErgiBAQURSKCCxn5sedaIPtJEyBKsAVI5JUitGdNIyyIqWkj1ZlbvOoayieYkRzuIYka8g9KBZBianOZosAcMU16AacyBoiAPcD2MzMv7X78x0j3pEOFmsWgQt8lNroPyNrqNCLi1ldqtLkGopUe1IIkqncE8pkHoE7YGZ0JSykj8Y7My0CwBXXoBtwwiKYBuA7AP6JiJrVv7McaIe9xLvSAuAii0C7kIzT9gtOyUuqKaJZrqFoep+HiCdzZw3JzGJ30J1MgRkZE8qKXrNYLIKcRHo/pLww8+sAqNcD/YYxfdSFFoHR79oVT4GZoRhvPaDNFcgVLO5uL3dTK04ime0aUmYWO9QgIYOubqWv1hqyhpIpRjzHjPCcGF1DYhFkIDOL7cI4q9GFFoGx1hCAwgLGukWQQwgS3nQNZQWLCWIRuAQtQ8g4aAGKKJue6Epfe2IRZCBCYBcZoxH3WgSaIBSUjZFSMy6yLALvBYuZlaB5dvXRkAiBS8jlxlS2F+geMs/lAVxxDboBEQI7YM4MFmujEhcUvcpnERSUn503RuC99NFuddEBzf+sURUhfZ/gLFrsqlYvmV6ERZCMK4kNEfdl7rkBEQI70OYLmEcjLuiEuWIEQIEXV94YQbXngsXd6u9QZXINVYVD+j7BWbSRv7mvxgoZtOhLxYpFkAsRAjvQO6HZInC+E2ZbBEWY23ljBN5zDelCEMm8JKqj4cJuNELF0fuqYc6Lsr2Avqpday6M07kBEQI7MK5FALjbIihm3WLN0tG+l4aHXUNmIRCLwD2YXUNFWa/6UrHuy9xzAyIEdqBZBFX1yqOrLYJiLi7TKEvDwxZBtdkiiISQYqX8hOAssXi+rKEC/jfd2jVosspdMBhzAyIEdqDl1GtC4GaLQHMNFTIK1i2C6sztHrQIYnlcQ9rrouovCRUhO0aguYYKGLRo12DUdA26YDDmBkQI7MALMYJoGS2CSLWSWuqhdYvzBYs1C0HcQ86T1zVUULBYG4yJRZALEQI76Da5hsIRZTEXF3TCjm5lLkB9lZICWtQoK69FoAaPPbRKmV5zKSt9VHktKaTOkzWPIFKCa0gfjFUDIFcMxtyACIEddB9RHjUhAFyzVF57LPcoq6AJZQlTEFxDq0bqIfdQ3vTRYibYCRWlMytGUMSgRY/T9VEeXbh2uJOIENiB2TUEuGapvI7uBOqqwgipM2qLMrd7ihEAniozoaWIZqWPaq6hpKSQOo028tf+J0VNKOs2uYYA160U6CQiBHZgDhYDrrEIjsSSqK9OzwwuW9YQ4EmLwJw1VFVMrrpQUTq7E6iNhvViiJpFUFAgXw8WG4TAZWuHO4kIgR243CKor0r7xYuapJN3HoFqISScF7pC0UtM5LUIRAic5kgsiT416UFLVTgEokJdQzkGY2IR6IgQ2EHO0Yg7OmF7LIm6qvTFFQmHEAlRkRaByTWkXWyaAHqAfDOLqyRryDUciSXQx2C9EpG6gH0hrqEOgMKZq+lFaj01WKkkIgR20N2uiIBx0ZOoO8zS9lgC9dWZmTLamgS9ks8i0ISg2ztCYA5EalQXM9NaqCi5+mptVYF9Nd6h9EvjGhvRGlcUfnQDIgR2YKw8quESi6CjO5ERIwCKWKUs0aWOskzVR3UhOFKmVlaeDjV7yvxbaDeejm4RAqc5Ekvoac4aBa+x3d2e6RYCxCIwIEJgB90dmdkKgCIMLhiNtHcnsy6u6khYn87fI8ZV14xo2zzkGjoSU+ZT1JksAu230fYLztFucg0BigVX0KBFs8qNRGs91UcriQiBHXR9CdT0z9xW01/Z7jDtMSV91Eif6ggOF3Lj6zoEVPfN3q7lantouUrtdwiZFqbRbjztIgSOcySWyAgWA0DfmggOdxXSV917DboBEQI76DwA1DZkbqsbCHTud6Y9Bo7Esl1DA+qiONhRQOpn7BBQ0y97u2b9eEkIupNZvwOQdhWJEDhPe86+WlVYX+06CNQOyNxW2wB0OH8NugERAjvoOgjUDMjcVjtQ8aE7OOkqkUzhcFcC/Wsz1xMYWF+FAx0FlIeIHc5tEUQ9KASxzDRajapICFXhEI7EJEbgJMyMLzvjWX11QF20sL6abzDW9SWQkv+tCIEd5OqE2uik84DtzdH4slO5gAbWV2VsL3iUlU8IPJg+mmu0qVFXHRaLwGGOxBKIJxkD6zL7akNdFQ4U0lc7D+YejIGVfQFHhKDSMKtCMCBze91A5dFB95B2ATXUmy+uKA52xMHcy6LtsUNAdQ7XUDiqZGR4yP/aniN7SqO+KiJC4DD723P31QF1URzuSvS8XgSz0hezBmPqawcHY25BhKDSxDuUUgs5RyNwtBPub1ctghyjrESKew8Y57MIAKD+KKCjrRzNtIVDnQn0zSMEfWsiOFRIQFKoGJoQDKzPdA01qH1Xs25zEjukLFyfyzUEuCJW5zQiBJXm0F7lsd+wzO1ap3QwWLW/XZkQNqAu2+8KAAfbe7i4WDWpzZkYGvVHAe0t5WimLexv78agPlU59zX2rUbrkZjNLRKMaNbrwPrMWexaX+0xTnDoM+Wx79GZ211wDboFEYJKc2i38mgWgr5D1f2f2dseA58dVCbTDBtQm7G9sa9ysX1xuIfJNrHDyszoPkNy769v9IwQMDPa2mMY1Kc65/7GvtVoOSxC4CRaXx3cN/N/pPfVQz301S/3KI/9mzK3a333sHPXoFsQIag0Wifsd0zm9j6DleyaAx/b3yaVT/d3oE91BA0mi+Arg5Rg7862HoK9hz9XHs2jLI36RqC9tRzNrDgHO+KIJxmD6nNbBIP71qDlcKz3mIlQMXbt70BVJISj+2WWMymor+YbjPUbptQe2u/cNegWRAgqTds2ZTUy82iECGgY4Wgn/HR/B5oaavWyvhrHDKhFOET4pLWH9M+DO5VH8/fS6H8scHivJ+oN7dyvtFG7qZgZ0q8a3ckUWo94p6y239jZ1oFjG2qzJvwd3a8GVeEQdrb10FfbtikVcTUrXCMUBgZ8xdHBmFtwRAiIaB4RbSWibUR0nRNtsI3PNgGDT8qu0AkAQ0YDezYq/nabYWY07zqIk4ZlZ/1URUL4h8Z6vLv7YP4T7PtQeTxqVO79Q04GOAW0bLHe2Arz0eeHAQDHN+YWgv91tPIbvf+Zd7Kg/AQz451dB3Di0Oy+Gg4RRh7dB827DuY/wZ53gKPHpNfJMDLkZGDPJkeuQTeRO02ighBRGMCdAOYC2A1gPRE9xcwf2t2WDLSOoHeInl73fCyzmsp26DNg51pgwiJwio1HKoyYgfB7f0Hy7y+A/2F2xj5zv2TD3p76rHEfZ36a3gaA8fbH+3GwvQtfO74BSCb07doJZhzfD39evwu79u3HsQ11md+RU8DWZ4BBJwD1g3I35JjxyuPmv6aflxmjqybze+c5JmO7to2x+oPPMbC+CsflsQjGNPVHJER47oMvMHNkY5YF5VXMrq7sPtfL8T2+N8e5Mw5K96f0Zu0aYRivreZPD2L/oXb844gR6Yq3hv2nHdcHv1+7E5/sbcEI/X+o7v9yN7BrHTDlCuTk+NOAD58EPn5VeR5QyG6/JxFNBXAjM5+hvv4xADDzLfneM35YlF+/VInwk/oPTj8i4zXybDcfHyJ7vnc7V2N+9y/wCQ/N2leDGF6o/iGaqBUxjiCJMBi9tb3n72/cXunv+PPUIjzE8/TX5pvxHaFbcXpoPbo5jDgiYLj3BqrNIM5HLJFEPJn5m5eKnVccGT6NcmzPt9/cz4p9v13XVyG0cw2+wf+JPRicta8OnXiK/g1DqQ0xjiKJkKv7abH0uXnfRmae2NtxTgjBNwHMY+bvqa+/A2AyMy81HXcZgMvUl6MBvG9rQ93LUQC8EYWtPPJbpJHfIo38FmlGMXOeyT5pbHcNFQoz3wvgXgAgog2FqFoQkN8ijfwWaeS3SCO/RRoi2lDIcU4Ei/cAONbwukndJgiCIDiAE0KwHsAJRHQcEVUBuBDAUw60QxAEQYADriFmThDRUgCrAYQBLGPmD3p5272Vb5lnkN8ijfwWaeS3SCO/RZqCfgvbg8WCIAiCu5CZxYIgCAFHhEAQBCHguFYIiOg/iWgLEf2NiJ4gogGGfT9Wy1NsJaIzHGymLRDR+UT0ARGliGiiYfsIIuokomb17x4n22kH+X4LdV+g+oURIrqRiPYY+sJZTrfJbgJVuqYXiOgTInpP7Qu9ppC6dh4BgOcB/FgNLv8KwI8B/IiIToKSaXQygGEAXiCikczs54VH3wewAMB/59i3nZnH2tscR8n5WwS0X5j5L2b+jdONcALXlq5xllnMXNDEOtdaBMz8HDNrRXDegjLfAADOBfBnZo4x88cAtgGY5EQb7YKZNzPzVqfb4QZ6+C0C1y+EDCYB2MbMO5i5G8CfofQJoQBcKwQmFgN4Rn1+DIBdhn271W1B5TgieoeIXiGi6U43xkGkXwBLVVfqMiJq6P1wXyH//0wYwHNEtFEt19MjjrqGiOgFALlWNvkJM69Uj/kJgASA5Xa2zW4K+S1ysBfAcGZuI6IJAJ4kopOZ+VDFGmoDJf4Wvqen3wXA3QB+BuUG8DMAt0IZQAnB5B+ZeQ8RDQbwPBFtYeZX8x3sqBAw85ye9hPRIgBfBzCb0xMefFmiorffIs97YgBi6vONRLQdwEgABdUXcSul/Bbwab8wUujvQkT/D8CqCjfHbfj+/18MzLxHfdxHRE9AcZ3lFQLXuoaIaB6AawGcw8zGZa6eAnAhEVUT0XEATgDwthNtdBoialSDZCCi46H8FjucbZVjBLpfEJGxzvk3ELxqvVK6RoWI6omor/YcwOnopT+4OWvoDgDVUMwaAHiLma9g5g+I6C8APoTiMlri98wQIvoGgNsBNAJ4moia1fUcZgC4mYjiAFIArmDm/Q42teLk+y2C2C9M/JqIxkJxDX0C4HJHW2MzJZau8StDADyh3jcjAP7EzM/29AYpMSEIghBwXOsaEgRBEOxBhEAQBCHgiBAIgiAEHBECQRCEgCNCIAiCEHBECATBBoioVi0DEu7hmJfNVVOJ6GoiuludM9JjCqAglIoIgSD0AhGVY77NYgCP9zK34WEoE6GMXAjgYWZuAbCXiKaVoS2CkIEIgeAriOhiInpbrcP+34aZ10eI6D+I6F0ieouIhqjbG4noMSJar/5NU7ffSER/JKI3APxRPe55dS2E+4hoJxEdRUQ3E9HVhs//DyK6KkfTvg1gpeG4H6qf9zciuknd/CiA+erMWBDRCCgltV9T9z+pnkcQyooIgeAbiOhEAAsBTFPXaEgifeOshzI7/RQoNVcuVbffBqWO/6kA/hnAfYZTngRgDjNfBOAGAC8x88lQbtjD1WOWAfiu+vkhKCP4h0ztqgJwPDN/or4+HUoJjEkAxgKYQEQz1FnhbwM4U33rhQD+YqiztQFAkCvMChXCzSUmBKFYZgOYAGVREgCoBbBP3deNdCG2jVAWMAGAOQBOUo8HgH5E1Ed9/hQzd6rP/xFKDR8w87NEdEB9/gkRtRHROChT+99h5jZTu44CcNDw+nT17x31dR8owvAq0u6hlerjvxjetw+KhSAIZUWEQPATBOD3zPzjHPvihpF1Eum+HwIwhZm7Mk6kCEN7gZ97H4BFUEpEL8uxvxNAjamdtzBzrhXnVgL4LyIaD6COmTca9tWo5xKEsiKuIcFPvAjgm2oNdhDRQCL6Si/veQ7AldoLtXBbLt4AcIF6zOkAjAu/PAFgHoBToRQ9y4CZDwAIE5EmBqsBLNYsDyI6RmszMx8B8DIUQXnYdKqRCF5VUcEGRAgE36CuT3s9lJWZ/gZl3euhPb8LPwAwUQ3afgjgijzH3QTgdCJ6H8D5AD4HcFj93G4oN++/9JAV9BwU9xKY+TkAfwKwlojegxJz6Gs49mEApyBbCGYBeLqX7yMIRSPVRwWhAIioGkBSLXc8FcDdakBaCxJvAnA+M/89z/vHA7iGmb9joQ2vAjhXtTAEoWxIjEAQCmM4gL+oN/1uqFlHRHQSlCD0E/lEAACYeZM6YSxcyjoJRNQI4LciAkIlEItAEAQh4EiMQBAEIeCIEAiCIAQcEQJBEISAI0IgCIIQcEQIBEEQAs7/Bxiiduq76xi9AAAAAElFTkSuQmCC)
%% Cell type:code id: tags:
``` python
%%bash
wget -N -q http://www.west-code.org/doc/training/silane/wfreq_spec.in
```
%% Cell type:code id: tags:
``` python
%%bash
cat wfreq_spec.in
```
%%%% Output: stream
input_west:
qe_prefix: silane
west_prefix: silane
outdir: ./
wstat_control:
wstat_calculation: S
n_pdep_eigen: 50
wfreq_control:
wfreq_calculation: XwgQP
n_pdep_eigen_to_use: 50
qp_bandrange: [1,5]
n_refreq: 300
ecut_refreq: 2.0
ecut_spectralf: [-1.5,0.5]
n_spectralf: 1000
%% Cell type:code id: tags:
``` python
%%bash
mpirun -n 2 wfreq.x -i wfreq_spec.in > wfreq_spec.out
```
%% Cell type:code id: tags:
``` python
import json
import numpy as np
import matplotlib.pyplot as plt
with open('silane.wfreq.save/wfreq.json') as json_file:
data = json.load(json_file)
k = 1
kindex = f"K{k:06d}"
print(kindex)
print(data["output"]["P"].keys())
bandmap = data["output"]["Q"]["bandmap"]
eqp = data["output"]["Q"][index]
freqlist = np.array(data["output"]["P"]["freqlist"], dtype="f8")
spf = data["output"]["P"][index]
for i, b in enumerate(bandmap) :
print(f"band = {b}")
eks, eqpLin, eqpSec = eqp['eks'][i], eqp['eqpLin'][i], eqp['eqpSec'][i]
sigmax, vxcl, vxcnl = eqp['sigmax'][i], eqp['vxcl'][i], eqp['vxcnl'][i]
sigmac_eks = eqp['sigmac_eks']['re'][i]
sigmac_eqpLin = eqp['sigmac_eqpLin']['re'][i]
sigmac_eqpSec = eqp['sigmac_eqpSec']['re'][i]
z = eqp['z'][i]
print(eks,eqp['sigmac_eks']['re'][i])
bindex = f"B{b:06d}"
sigmac = np.array(spf[bindex]['sigmac']['re'], dtype="f8")
plt.plot(freqlist,sigmac+sigmax-vxcl-vxcnl,'b-',label=r"$\Sigma_{xc}(E)-V_{xc}$")
plt.plot(freqlist,freqlist-eks,'r-',label=r"$E-\varepsilon_{KS}$")
plt.plot(freqlist,sigmac_eks+sigmax-vxcl-vxcnl+(1-1/z)*(freqlist-eks),'g-',label=r"$\Sigma^{lin}_{xc}(E)-V_{xc}$")
plt.legend()
plt.title(kindex+" "+bindex)
plt.xlabel("frequency (eV)")
plt.ylabel(r"$function (eV)$")
xmin, xmax = min(eks, eqpLin, eqpSec), max(eks, eqpLin, eqpSec)
ymin, ymax = min(sigmac_eks, sigmac_eqpLin, sigmac_eqpSec), max(sigmac_eks, sigmac_eqpLin, sigmac_eqpSec)
ymin += sigmax - vxcl -vxcnl
ymax += sigmax - vxcl -vxcnl
plt.vlines(x=eks,ymin=ymin-0.2*(ymax-ymin),ymax=ymax+0.2*(ymax-ymin),ls="--")
plt.vlines(x=eqpLin,ymin=ymin-0.2*(ymax-ymin),ymax=ymax+0.2*(ymax-ymin),ls=":")
plt.vlines(x=eqpSec,ymin=ymin-0.2*(ymax-ymin),ymax=ymax+0.2*(ymax-ymin),ls=":")
plt.xlim([xmin-0.2*(xmax-xmin),xmax+0.2*(xmax-xmin)])
plt.ylim([ymin-0.2*(ymax-ymin),ymax+0.2*(ymax-ymin)])
plt.show()
#for key in eqp:
# if("re" in eqp[key]):
# print(f" {key} = {eqp[key]['re'][i]:.3f}, {eqp[key]['im'][i]:.3f} eV")
# else:
# print(f" {key} = {eqp[key][i]:.3f} eV")
print()
```
%%%% Output: stream
K000001
dict_keys(['freqlist', 'K000001'])
band = 4
-8.230102882308456 0.012646784475642808
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAY0AAAEWCAYAAACaBstRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABTvElEQVR4nO3dd3hURRfA4d9JgVClC0gNzUCAAKEjvVkDCqhIl6b0qoAgHQUVKYogRQUEBBXlw0S69BZ6h0SagPRO+nx/3A0G2Cy7ybYk8z7PPiR3bjkLZE/unZkzopRC0zRN06zh4eoANE3TtJRDJw1N0zTNajppaJqmaVbTSUPTNE2zmk4amqZpmtV00tA0TdOsppOGpmmaZjWdNDS3JyKnRaRhgu/fEpEbIlLH9H2AiISKyH3TnwEJ9hUR+VRErplen4qIJGh31LH1RGS9iNwSkdNPeX9FRESJyF3T618R+VpEvBPsk0NEfhWReyJyRkRaP3aO1qbt90RkuYjkcMaxCfaba3oPxS29Vy3l00lDS1FEpD3wFfCyUuovEUkH/AYsALID3wO/mbYDdAWaAeWBcsCrQDfTuRx57D1gLjDIhreXTSmVGSgLVAd6JGj7CogCngXeAWaISBlTLGWAmUBbU/t94GtHHxtPRGoBxWx4n1pKppTSL/1y6xdwGmiI8YF9FQhM0NYY+AeQBNvOAk1NX28FuiZoexfY7uhjE2xrCJx+yvsrAijAK8G2icAs09eZMD64SyZonw98Yvp6PPBjgrZipv2zOPJY0/dewF6MpKqA4q7+/6Jfjn3pOw0tpXgPGA00UErtTrC9DHBAmT7BTA6Ytse370/Qtv+xNkcdm2Qikh9oAmw3bSoJxCilTliI5WGcSqkwTB/2Dj4WoB+wUSl1wPZ3qqVEXq4OQNOs1AhYDxx8bHtm4NZj225h/KZsrv0WkNnUN+HIY5PiqqnL5BlgG7AsQRy3bYgzYXuso44VkYIYd3+VLL8tLTXRdxpaSvEexm++sxN2RgN3gayP7ZsVuJNIe1bgrukOwZHHJkUupVQ2ICOwBfgzkTieFmfCdkce+yUwWin1eNLRUjGdNLSU4l+gAfACj3bUHgbKPZZIypm2x7eXT9BW/rE2Rx2bZEqpB8B3QDURyQWcALxEpISFWB7GKSK+QHrTcY48tgEwSUQuicgl07ZtiY2w0lIJV3eq6Jd+Pe2FqSPc9HUh4G9gsun7dMAZoA/Gh11P0/fpTO3dgaPAc0B+jA+87k441gPwAV40bfeJbzPz/oqQoCPcdL5PgIuYOtqBxcAijM7pmhiPicqY2spgPEZ6wdS+AFic4PyOOjYPkDfBSwHVgAyu/j+jXw78eXR1APqlX097JUwapu+LAueACabvKwChwANgD1Ahwb6CMRLpuuk1kUdHPDnq2LqmD9GErw2JvL/4pHHX9LoJ/AVUTrBPDmA5xlDes0Drx87R2rT9HsZQ4BzOOPax8+jRU2ngFf9bjKZpmqY9le7T0DRN06ymk4amaZpmNZ00NE3TNKvppKFpmqZZLVXPCM+VK5cqUqSIVfveuwcnT4IIlCgBGTM6NjanOHIEPD2hVClXR6JpaUb4lXsA+ObO5OJIki40NPSqUiq3ubZUnTSKFCnC7t27n76jybFj0KQJnD8Pv/0G9eo5MLinuHDzAQD5s2Wwqe2hixchf36YMAE+/NAhMWqa9qQ3Z24DYEm36i6OJOlE5ExibfrxVALPPw9btkChQtC0Kfz8s+ti6bdkH/2W7LO57aE/TRUomja1a1yapqVtqfpOIykKFICNG+GVV6BlS5gxA7p1c34cveqXSFLbQyEhkDcvlC//9H01TbMbq34+UzCdNMzIkQPWrDGSRvfucPkyfPSR0d/hLLVK5EpSGwAxMbBqFQQFOTdoTdOe/vOZwumkkYiMGWH5cnj3XRgxwkgcU6aAh5Me6J29dh+AQjmf7JG31AbArl1w4wa8+KLD4tM0e4uOjub8+fNERES4OpRkiY6NA8Db0/2f/vv4+FCgQAG8vb2fvrOJThoWeHvDd99B7tzwxRdw5Qp8/z2kT+/4aw9aZqyLY64zzVIbAMHBRnZr2NB8u6a5ofPnz5MlSxaKFCmCpOA75LArdwEoljuziyOxTCnFtWvXOH/+PEWLFrX6OJ00nsLDAz77DJ59Fj74wLjj+OUXyJbNsdft16hkktoAoz+jWjXjOZumpRAREREpPmGkJCJCzpw5uXLlik3H6aRhBREYPBjy5YNOneCFF+CPP6BgQcdds5pvziS1ceUK7N4No0Y5ICpNcyydMJwrKX/f7v/QzY20bWs8+TlzBqpXh4OPLzxqR2FX7j68zbWljVWrQCk91FbTNIfQScNGDRvCpk3G53KtWrBunWOuM/SXgwz9xXxWstRGSAjkygWV9LLNmqbZn348lQTly8O2bcbgpKZNjc7y1nZe4HJw08RLfyTaFhdnTOpr0sR5w7w0TXtE3qw+rg7BofQnSxIVKgSbN0ONGvDOOzBxonH3YS+VCuegUmHzHdmJtu3ZY/Rp6KG2mpZkXl5elCxZkoCAAAICAhg9erRVxz148IA6derg4yVkSu/FzJkzyZcv38PzBAQEsH//fmrXrk1MTEyS4wsLC6Ns2bKPbIuMjKRo0aIcPpzsJeqfSt9pJEP27MYv9u3bGyOrzp415nJ4eib/3Mcv3QGgVN4s1rcFBxu99o0bJz8ATUujPv30U7Zs2cLPP/9sU0fx3Llzef3114mIURATw8GDBxk7dizvvvvuI/s1aNCAJUuW8M477yQpvqJFi3L+/Hni4uLwMD1RmDVrFrVr16ZMmTJJOqct9J1GMqVPDz/+CAMGwFdfGbPIHzxI/nlH/HaIEb8dsq0tJAQCA42JJZqmJcmAAQPIkCEDY8eOtem4hQsXEhQUxKXbEVy6HcGBAwcICAh4Yr9mzZqxcOHCJMfn4eFBoUKFOH36NGDc4Xz++eeMctKISX2nYQfxczkKFoR+/aBBA/j9d6M/OqmGvuRnW9uNG7B9OwwblvSLapqb6NsX9u2z7zkDAuDLL63bd/bs2ZQvX566devywgsvPHX/qKgowsPDKVKkyMORjYcPH6Zjx44P7wbef/99unbtir+/P7t27UriuzD4+flx7NgxfH19+eqrr3j11VexdhmI5NJJw4769DEKHr7zjtHX8ccfULx40s5VvmA229pWrzY6wvVQW01LtqNHjxIbG0t5Kwt+Xr16lWwJZvxe+Oc8uXPn5sCBA0/s6+npSbp06bhz5w5Zsvz3iLlhw4ZcunTpif3HjRtHUFDQI9v8/Pw4fvw4tWvXZvr06ezYscPKd5Z8OmnY2RtvGMVlg4KMSdnLlxtDc211+MItAMrkf8a6tpAQo5OlSpWkhK1pbsXaOwJHiIyMpEOHDsyYMYOsWbMCcOXKFQYPHsyYMWMYMWIEM2bMYMiQIYgIhQsXpm3bto/UzDpx9DB+fok/LYiMjMTH59FRVmvWrLE6Rj8/P9auXcuUKVN45513ePbZZ5+IcebMmXh4eDBo0KCHcfbu3dvGv40n6aThADVrGk+KXn7ZeFQ1d65x92GL0SuOAObrSz3RppSRNBo1Ai/9T6ppyTFixAiqVatG4wQDSnLnzk2hQoUYMGAAc+bM4dtvvyUoKIg6deo83Cc2NvZh4jh25BDPP/+82fNfu3aNXLly2VQk8HF+fn5MmDCBNWvWEBoaajZGb29vpk+f/kScyaU/YRykeHFjLsfrr0ObNnDqlFEt19rBGCNeLW1924EDxkp9eqitpiXbZ599RvHixR92YteuXZvx48cTHh6Ol5cXmTNnJjQ0lO7duz9yXOPGjdm8eTM1a9fl3KljbN60ieDgYMAo17Fp0yYyZ87M+vXrefnll5MVY8mSJTl48CDjxo17+Fjs7t27j8QImI0z2ZRSqfZVqVIl5WqRkUq1a6cUKNWmjVIREQ64yIQJxgUuXHDAyTXNOY4cOeLqEMyKjo5WHTt2VKdPn1YTJ05U69evV8uXL1edOnVSAwYMUNeuXVNKKRUaGqratGnz1PM1b95cHT9+3OExKqXMxvk4c3/vwG6VyOeqKHvOSHMzgYGBypY1wh1FKRg3DoYPN4od/vor5LRQcxBg/7mbgPlO7yfa6taFW7dg7157haxpTnf06FGL/QApwdy5c3n9zdZ4enqSxefJx09RUVEsXryYdu3auSA688z9vYtIqFIq0Nz+ep6GE4gYK//9+CPs2GEUOzx50vIx4/84yvg/jj697fZtY2FzPWpK01yuU6dOXLsfw+U7kWbb06VL51YJIyl0n4YTvf22UX4k4ciqxIaAjw7yT/Q8j7StXWss76qThqZpTqDvNJwsfmRVrlxGxdwFC8zvVypvFrMlRJ5oCwmBLFmMiSGapmkOppOGC8SPrKpRw1ij4+OPjXl5CYWeuU7ometmj3/YppRRb6phQ2NtWk3TNAfTScNFcuQwih126ACjR8Nbb8H9+/+1Tww5zsSQ42aPfdh29CicO6eH2mqa5jS6T8OF0qUzJv6VKWMsJxsWBr/9ZpQiGf962USPe9j2w0zjT92foWlu47lsGVwdgkPppOFiIjBwIDz/vNFRXrmy0UFetWrmRI8pltvUFhJiZBxHLlauaZpNfLztsDaCG9OPp9zEK68Y/RwZMkCdOjDqyztsD79mdt/t4dfYeegsbNyo7zI0zc3cfhDN7QfRrg7DYVyeNERkjIgcEJF9IrJKRPInst9EETksIkdFZKrYsjpKCuHvDzt3QtWqMLJfFjr1fPBEBznA5NUnWD/jJ4iK0klD09zMlbuRXLlrfp5GauDypAFMUkqVU0oFAP8DRjy+g4jUAGoC5QB/oDJgvwpcbiRXLqPK+VttYjgaXIAWLeDu3Uf3mdSiPO9HnoKMGROf6KFpWpKYW6b14MGDrg7Lbbi8T0MpdTvBt5kAc3VNFOADpAME8Ab+dXx0rpEuHfz4gxfVAqF/f6O0+u+/GxMDAQrlyADrV0P9+sbSgZqm2U1iy7RqBpcnDQARGQe0A24B9R5vV0ptE5H1wEWMpDFdKWW+xkYCx68dZ+u5rdQomPImvm05dZVKL8HKUrl4883/OsirV4fda3YSGB5urDGraamRC5fuO3DgAB07dkz25UJCQvjwww8BSJ8+Pdu2bXu4il9K5pR3ICJrROSQmVcQgFJqmFKqILAQ6Gnm+OKAH1AAeA6oLyJmn8uISFcR2S0iu+9H3afm3JoELQ7i0GXz6227q2nrTjJt3UmaNjVmkGfNatQlnD8fDs5ebOyk+zM0ze7il2mNfzQ1a9asJJ2nV69eBAcHs2/fPnbs2JEqEgbgXqXRgULAITPbBwHDE3w/Ahj8tPNVqFhBjds4TmWdkFXJSFHtfm2n/r7xt4UCw+7jnxv31T837j/8/to1perXNyqgHyzcVEUVK+HC6DTN/tyhNPrZs2dVqVKlLO7ToEEDVaZMmSdey5cvV0opFRkdqyKjY1Xv3r1Vrly5VJ8+fR4e26xZMxUXF6fGjRunvvzySxUZGal69uyp+vTpo9q3b6/Onz/vyLdnlq2l0V3+eEpESiil4mu+BgHHzOx2FugiIhMwHk/VAb582rk9xIOhLwylW6VufLL5E6btnMbiQ4t5L/A9hr4wlDyZ8tjrbdhd/scmCOXIYUzL+KD3A3y/+Yv/FexC7WtPL7GuaZr1Dh48+NTy7E9bljWdlwdbt25FKcXFixfxMq2mee7cOTJnzsy7775L+/btqVOnDqtWrSJ37tyMGPHE+B+35Q73S5+YHlUdABoDfQBEJFBEZpv2WQaEAQeB/cB+pdQKay+QM2NOJjWexKnep2hXrh3Tdk6j2NRijNwwktuRt59+AhfYcPwyG45ffmSbtzd80WwjGXnA3AtNqFzZWLRP0zT7OHDgQKLLtFrr5v0oFvy4mJIlS+Ll5YVSitu3bxMaGsr69etp1arVw+VXa9SoQWRkJG3btmVBYtVL3U1ityCp4ZXYyn1HrxxVbyx5QzESlWtiLvXlti9VRLQjltRLulbfbFWtvtn6ZEPfvirSO516sf9OlT+/UhkzKrV0qfPj0zR7c4fHU61bt1bPPfecKl++vCpfvrwKCAhQd+7csekcpy7fUT+HrFf+/v6qbNmyqkqVKmr37t3qo48+Utu2bVNBQUHq0qVLjxwTGxurqlWrZs+3YjVbH0+5/IPdka+nLfe68/xO1eD7BoqRqEKTC6nv9n6nYmJjLB7jLP/efqD+vf3gyYZSpVREg0bq39sP1IULSlWvbvwrDhmiVIx7hK5pSeIOScMeTl2+o05dfjLRNG/eXEVHR6tz586poKAgdf/+ffXee++pXr16qY4dO6rff//dBdHq5V4fYe1yr2vC1/Dhmg8JvRhK6dylGV9/PK+Veg23m3T+99/g62sMG+zTB4DISOjVC7791ih2++OPYFpnXtNSlNSw3CtA2BVjNu7DGnFuTi/3mgQNfRuyq8sulrZcSkxcDM2WNKPm3JpsPLPRZTGtOfIva448Nn8xJASArcUDH7alTw+zZsE33xgzyatUgSNHnB2tpmlphU4aJiJCi9ItOPz+YWa9Moszt85Q57s6vLTwJfZd2uf0eL7dFM63m8If3RgSAkWLMuWcPNHWrRusX28sGV61qjERUNM0zd500niMl4cXXSp14VSvU0xsOJHt57dTYWYFWv/cmrDrYU6LY0abSsxoU+m/DVFRxnrgTZsyo23go20mtWrB7t1GmfXmzWHkyCdXBNQ0zbEK58hI4RwZXR2Gw+ikkYgM3hkYVHMQ4X3CGVprKMuPLef5r56nx8oeXLxz0eHXz5EpHTkypftvw+bNcO8eNG36ZFsCBQrApk3Qvj2MGmUkj1u3HB6upmkmXp4eeHmm3o/W1PvO7CSbTzbGNRhHWO8wulTswqw9syg+rTjD1g7jZsRNh1035NBFQg4lSE4hIcZEjfr1n2x7jI8PzJsHU6bAH39AYCDoIp2a5hzX70Vx/V6Uq8NwGJ00rJQvSz6+fvlrjvY4SlCpIMZvHo/vFF8mbZnEg+gHdr/evC2nmbfl9H8bgoONMuiZMz/ZZoYI9O4N69YZpdWrVTNGVmma5lg37kdx475OGppJ8RzF+fGNH9nTdQ9VC1Rl8JrBlJhWgtl7ZhMTF2O363zbPpBv25tGvJ0/D4cOGWNqH297ihdegD17oGJFeOcdI5FEpd7/z5qmOZhOGklUIV8Fgt8JZn379RR8piBdVnTB/2t/lh1Zhj3mvmT18Sarj7fxjWmobXxV20farJAvn3HH0a8fTJtmLMNx4UKyQ9Q0LQ3SSSOZ6hapy9ZOW1n+5nI8PTxpubQlVWZXYW342mSdd8X+C6zYb/pkDwkxerjLlHmyzUre3vDFF7B4sbFMQcWKxhLjmqZZp0OHDgAMHz7ctYG4mE4adiAiBD0fxIHuB/gu6Dsu37tMw/kNaTS/EbsvPH1GujkLtp9hwfYzEB1tzNpr2tToqEjYlgRvvgk7dsAzzxh3HF98Aam4KICm2czLy4uSJUs+XE9j9OjR3L9/n0yZMnHp0iWio6P5559/qFevHpMnT+bNN998eOyDBw94O6gpsbGxiS4bGxUVRe3atYmJSdrj7LCwMMqWLfvItsjISIoWLcrhw4eT9d6tklh9kdTwelrtKUd5EP1ATd42WeWamEsxEtXipxbq2JVjNp3jfmSMuh8Zo9TGjUZxqWXLnmxLhlu3lGre3Dh1q1ZK2ViTTdPszl1qT3322WeqefPmKi4u7uG2TZs2qdmzZ6vg4GC1aNEitXLlSjVr1iyllFHkMN706dPV5MmTVWxsnOrRo4eaPXu22WuMHDlSLViwIEnxxcbGqmzZsqnY2NiH26ZOnaratWuXpPPZWntK32k4gI+XD32r9SWsdxgf1/mYkFMhlPm6DF1XdOX87fNWnSNDOk8ypPM0Hk15ekLDhk+2JUPWrPDzz/Dpp7BsmVF+5Ji5lUw0LY0ZMGAAGTJkYOzYsQ+37dq1i8qVK7Nv3z4CAgLYv38/L7xgLB6asEbdwoULadasGR4ewoEDBwgICDB7jWbNmrFw4cIkxefh4UGhQoU4ffo0YNzdfP7554waNSpJ57P5+k65ShqVNX1WRtYdSVjvMHpU7sF3+76jxLQSDF49mOsPrls89te95/l173ljqG2NGsbzpMfbkkkEBg82nn5dvWqsQ75sWbJPq2kp3uzZs5k/fz6bNm0CjCVgy5Qpw8mTJylZsiSnTp2iZMmSXL16lbx58wIQFRVFeHg4WXLl49rdSIvLxvr7+7Nr164kx+fn58cx0295X331Fa+++ipFihRJ+hu2gctX7ksL8mTKw5QXp9Cvej8+3vAxn239jFmhsxhcczB9qvYhU7pMTxyzeOc5nrl1jeZ798L48U+0ATSvUMAu8dWvD6Gh0LKl8erdGyZNgnTmJ51rmsP1Delr95pvAXkD+LLpl1bte/ToUWJjYylfvjxgJBGAOXPmPPJnrly5+OyzzwC4evUq2bJl4+aDaC78c57cuXNzIJFV0jw9PUmXLh137twhS5YsD7c3bNiQS5cuPbH/uHHjCAoKevi9n58fx48fp3bt2kyfPp0dO3ZY9b7sQScNJyqSrQjfN/ueQTUGMWzdMIatG8a0ndMYUXsEnSt2xtvzv2G0CzpXRX743vjGNNQ2YZu9FSxojKb64AOj8vr27fDTT1C4sN0vpWluLTIykg4dOjBjxgyyZs1q9XEZMmQgIiICgBNHDz+1zHtkZCQ+Pj6PbHvaUrLx/Pz8WLt2LVOmTOGdd97h2WeftTrOZEussyM1vFzVEW6tzWc2q1pzaylGoopNKaZ+PPCjio37r3NLvfWWUnnzKpWgQ84Zli1TKmtWpbJnV2rFCqdeWkvD3KUjfPDgwapLly5JOrZAgQLq8NkrauBHI9WHH36Y6H5Xr15VpUqVSmqI6sCBAyogIED5+vqqGzduPNx++fJl1aFDB3Xu3DnVsWNHFRERofr166f69++vpkyZYvZceuW+FJQ0lFIqLi5OrTyxUpWbUU4xEhXwTYAKPhmslmwLVxHPZFOqffsnjvlp11n1066zDo3r1CmlKlQw/ocMHqxUdLRDL6dpbpM0PDw8VMmSJR8u+dqrVy+rP4w7deqkvl/6u3rt9VYWl41dunSp6t+/f5JjjIiIUJ6enuqTTz55om3EiBGqVatW6s6dO2ratGlqw4YNFs9la9LQj6dcTER4qcRLNC3elMWHFjN8/XBeXPgihWJLUjDLTaqbSocktCzU6ARvGVjQYXEVKwZbt0LfvjBxImzbBosWwXPPOeySmuYWYmNjzW4vVKgQAwYMYM6cOXz77bcEBQVRp06dR/bp0aMHoydM4otv5lhcue/HH3/kk08+SXKM6dOnNzvP4+7du4SHh+Pl5UXmzJkJDQ2le/fuSb6OOXr0lJvwEA9al23N0R5Hmf7idCI9/qFGZ2gW+T2HLz86YWdJt+os6Vbd4TH5+BgrAi5YYNSvqlDBGGmlaWmNuQ/jmjVrPrFfxYoVee3FRhTJkSHRc0VFRdGsWTNKlixp1xhjYmLo3bs3Y8eOJSAggA0bNtCsWTO6devGwIEDuX7d8ohNa+k1wt3U3RqBTCl6mYllb3En8g7tyrdjVN1RFM7mmp7pY8egRQtjKdkRI2D4cGP6iKbZi7uuER4TE0PXrl35+OOP+emnn6hcuTK3bt3i999/J3v27AwdOpQcOXK4Oswks3WNcJ003NHVq6g8eTjYuR/PfTmUCZsnMH3ndBSK9wPf5/lM7ciaPidvVynk1LDu3YMePeD776FBA1i4EJw5aENL3dw1adjqyh1jBFXuLD5P2dM92Jo09OMpd7RqFaIUS3OXJmfGnHzW+DNO9jpJ23JtmbpzKj3X12LC5jHcibzj1LAyZYLvvoM5c2DLFuNxlS56qGmPuh0Rw+0I+y2T4G500nBHISGQKxdjxnR8uKngMwWZ/dpsDr13iKDnX+TgndkUm1qMqTumEhkT6dTwOnUyih5mzgz16sHYsZBI36GmaamMThruJi7OSBqNG4PHk/88frn9WNZqGTs676Dss2XpE9KHUtNL8cP+H4iNc94nd7lyxizyt982+jcaNdJrdGhaWqCThrvZuxeuXGFz8UDmbzttdpf5205z/Gwe1rRdw6o2q8iZMSftl7en/Dfl+f347zirnypLFpg/H+bONe48AgL+Wy9K07TUSScNd2P61F2c3Y81Ry+b3WXN0cusOXoZEaFRsUbs6rKLn1r8RFRsFEGLg6g1rxabzmxySrgi0LEj7N4NefMaK9IOHqyXlNWSJjUMzPEQwSNB5Vt3lpS/bz16yt3UqgUREcansI2iY6OZt28eo/4axYU7F3ipxEuMrz+e8nnLOyDQJz14AAMGwIwZRqn1RYvA19cpl9ZSgb///pssWbKQM2fOR8qNa46hlOLatWvcuXOHokWLPtKmh9ymFDduQK5cMHQojBmT5NPcj77P9J3TmbB5ArcibvF22bcZU28Mvtmd8wm+bBl07mysCDh7tlE5V9OeJjo6mvPnzz8s+qc5no+PDwUKFMDb2/uR7TpppBRLl0KrVrB5M3NVfgA61Sr6xG5zN/+daFtCNx7cYNLWSXy5/Uui46LpVqkbH9X+iLyZ89o/9secPg1vvWX0dXTtalTOzZD4JFlNSzWmrj0JQO8GJVwcSdLpeRopRUgIZMsGVauyNewqW8Oumt3NUltC2TNkZ3yD8ZzqfYrOFTrzze5vKDa1GB+t+4hbEbfsHPyjihSBTZuMUuuzZhmPq44cceglNc0tbDl1lS2nnv7zmVLpOw13oRQUKAA1axoLWTjAyWsnGbFhBIsPLSZHhhwMqTWEHpV7kMHbsbcAf/4JbdvC3bswbZoxz0M/stZSqzdnbgNwSn04R9F3GinBwYPGRAczVW3tpUTOEix6YxGhXUOpnL8yg1YPouT0kszZM4eYOMfNYG3SBPbvN1at7dzZeGx144bDLqdpmgPppOEu4ic4NGkCwKyNYczaGGZ2V0tt1qiYryIhbUJY124dz2V5js4rOuP/tT8/H/nZYUMe8+Uz7jjGj4dffoHy5XUJEk1LidwmaYjIABFRIpIrkfb2InLS9Grv7PgcLjjYmGad3+gA33PmJnvO3DS7q6U2W9QrWo9t727j1zd/xUM8aLG0BVVnV2Vt+Npkn9scT08YMsSoW5UunVGC5KOPIDraIZfTNJfInjEd2TOmc3UYDuMWfRoiUhCYDTwPVFJKXX2sPQewGwgEFBBq2s/iQ44U06dx5w7kyGFMckjGwizJERsXy/wD8xmxfgTnbp+jkW8jJjSYQKX8lRxyvTt3oHdvowBi1apGxdxixRxyKU3TbJQS+jQmA4MxEoI5TYDVSqnrpkSxGmjqrOAcbu1aiImBpq57S54ennQI6MCJXif4ovEX7Lm4h8BvA2m1tBUnrp2w+/WyZIF582DJEmOtjoAAo+S6G/wOo2maBS5PGiISBPyjlNpvYbfngHMJvj9v2mbufF1FZLeI7L5y5YodI3WgkBDjU7RGjYebvt5wiq83nDK7u6W25PLx8qFf9X6E9wlnRO0R/HHyD0p/VZpuK7rxz+1/7H69Vq3gwAGjzHqHDkYBxJs37X4ZTXOaT0OO8WnIMVeH4TBOSRoiskZEDpl5BQFDgRH2upZSapZSKlApFZg7d257ndZxlDL6Mxo0MB70mxy5cJsjF26bPcRSm71kTZ+VUfVGEd4nnB6VezBv3zyKTyvOB6s/4PoD+ywbGa9QIVi/3iixvmyZ0Um+yTmlszTN7vacucGeM6l3eKBTkoZSqqFSyv/xFxAOFAX2i8hpoACwR0Qen7L8D1AwwfcFTNtSvmPH4OzZJ4baTm9dkemtK5o9xFKbveXJlIcpL07heM/jtCzdkklbJ+E7xZcJmyZwL+qe3a7j6QnDhhmd5F5eULeusaxsTOpdy0bTUiSbk4aIZBIRu6wOrZQ6qJTKo5QqopQqgvHYqaJS6tJju/4JNBaR7CKSHWhs2pbyPTbU1l0VzV6UH5r/wP7u+6lduDZD1w2l+LTifLP7G6Jj7Tf8qWpV2LfPmAw4Zgy88AKEh9vt9JqmJdNTk4aIeIhIaxFZKSKXgWPARRE5IiKTRKS4IwITkUARmQ2glLoOjAF2mV6jTdtSvuBg8PODwoUf2Tx17cmHNWweZ6nN0co+W5bf3/6dTR03USx7Md5b+R5+X/mx+NBi4lScXa6RJYsxqmrxYjh61Ogk/+EH3Umuae7AmjuN9UAxYAiQVylVUCmVB6gFbAc+FZE29gjGdMdx1fT1bqVU5wRtc5VSxU2vefa4nsvduwd//WV2Fnj4lbuEX7lr9jBLbc5Sq1AtNnXcxP/e/h8ZvTPy9s9vEzgrkD9P/Wm3CYJvvmnMJA8IgPbtjU7za9fscmpNc5h8z/iQ7xkfV4fhME+dpyEiZZVSB5+yj7dSyu2maLn9PI2VK+GVV2DVKmO91BQqTsWx6OAihq8fzt83/6ZukbpMaDCBagWq2eX8sbHw2WfGsrK5chkrBbpwdLKmpXrJnafxm4jMM03AM8sdE0aKEBICGTMaD+5TMA/x4J1y73Cs5zGmvTiNI1eOUH1OdZovac6RK8kvbevpaVTL3bnTmAP54ovQo4dxo6ZpmnNZkzSeB/YAG0VkioikgHGsKURwsFFLw+fJW9kvVh3ni1XHzR5mqc2V0nmmo2eVnoT1DmNMvTGsDV9L2Rll6fhbR87eOpvs8wcEGAsa9u8PX38NFSsaiUTT3MmoFYcZteKwq8NwmKcmDaVUlFJqGuCHMcFup4iMEZGsDo8uNTt1CsLCEq1qe+FWBBdumV/BzFKbO8icLjMf1f6I8D7h9K3al0UHF1FiWgn6/9mfq/eTt86Ajw98/jmsW2csL1ujBowcqetXae7DGfOoXMnm2lMikgnoA3QEZiqlPnNEYPbg1n0a06dDr15G8kjlRZfO3jrLqA2j+G7/d2TyzsTAGgPpX70/mdNlTtZ5b9406lfNnw+VKxt/lipln5g1Lan0ehr/naSIiDQFOgOFgDvAePuEmAYFB0Px4qk+YQAUeqYQc4LmcOi9QzQq1oiPN3yM7xRfpu2YRmRMZJLPmy2bMRT3p5+Mm7YKFYzHVnporqY5jjXzNA6IyDXgV6A98AywFmgHJO9XxbQqIsKom2FhwSVL9WtSam0bv9x+/NzqZ7a/ux3/PP70DunN8189z/z984mNi03yeVu2NNawql3b6CB/6SVjPStN0+zPmjuNZkAupVQFpdTbSqnRSqmlSqlDSqkoB8eXOm3caDyQtzBu9Ob9KG7eN//Xa6ktJahaoCpr263lzzZ/kt0nO+2WtyNgZgArjq9I8hyP/PmNm7evvjKmvpQtC0uX2jlwTbOCb+5M+ObO5OowHMbqPg0RKQF8CDxQSvV0aFR24rZ9GvHDf65fN4bcpmFxKo5lR5YxbN0wTl0/Rc2CNfmk4SfUKlQryec8ftwoQ7Jrl1E1d9o0yJnTjkFrWipnr/U05gPLgNqmk/qLyA92iC/tCQ6GOnXSfMIAY45HqzKtOPL+Eb55+RvCb4TzwrwXeOXHVzjw74EknbNUKaPw4ejRxt2Gvz+sWGHnwDUtjbIlaXgopYKBWACl1CHA3yFRpWanTxuVbS30ZwCMW3mEcSvNT4yz1JZSeXt60y2wG6d6n+KTBp+w5dwWAr4JoM0vbQi/YXvFQm9vYwb5rl2QJw+89hp07Ai3bjkgeE1LYMgvBxjyS9J+4UkJbEkaF0SkKKbV9UREgAwOiSo1+9NUnPcpdTAiouOIiDZfANBSW0qX0TsjH9T6gPDe4XxQ8wN+OfoLz09/np5/9OTfu//afL6AACNxDBtmDMn19zeqtmiao4RfuUf4ldRbrsCWPo0iGOt4lwaGYSy36qGUaumw6JLJLfs0mjUzan///TeIuDoat3fhzgVG/zWa2XtmG6sKVuvHwBoDecbnGZvPtWuXUfjw6FHo1g0mTTIq6mqaPel5GiZKqdMYiaI34Av8BbS1R4BpRlSUsR74iy/qhGGl/Fny880r33C0x1FeKfkKYzeNxXeqL59v/ZyIGNtmxVeuDHv2wMCBMGsWlCsHGzY4Jm5NS62smafx8NNNKRWjlFqmlBqulPpaKRXx+D6aBVu2wN27VpVotVS/JrXXtjGnRM4SLG6xmNCuoVTOX5mBqwdSYloJ5u6dS0yc9cv7+fgYdxibNhkrBNarZ8wqv3/fgcFrWipi1XoaItJLRAol3Cgi6USkvoh8jzHpT3uakBCjh7Z+fVdHkmJVzFeRkDYhrGu3jvxZ8vPu7+9SdkZZfjn6i01zPGrWNJ4S9uplDMkNCICtWx0WtpaGlM6fldL5U29pPmvW0/ABOgHvYDyWuoHRAe4BrAK+VkrtdXCcSeJ2fRrlyhkLQqxb5+pIUgWlFMuPLWfouqEcu3qMKs9VYUKDCdQvaltSXr8eOnUylmofMMAYqmum8LCmpRnJ6tNQSkWYHkXVxKg51QCooJQqrJTq4q4Jw+38849R6+IpQ20164kIzf2ac/C9g8x5bQ4X7lygwQ8NaLKgCaEXQq0+T716cOAAdO5sPLqqWBF27HBg4JqWgtlSsLA+8A0wEGguIpVEJL3DIkttrBxqG2/48kMMX37I5ra0yMvDi04VOnGy10k+b/w5uy/sJvDbQN5c9iYnrp2w6hxZssDMmcYTxDt3jJLrgwYZ1V40zRZ9F++l7+LU+7u0LfM05gIrMNYF9wVGAGmrNzY5goPhueeMiQJW8PH2wMfb/D+Ppba0zMfLh/7V+xPeO5zhtYez8sRKSn9Vmu7/686FO9ZVMGzSBA4fhi5djCVmy5c3Os01zVoXb0Vw0Y3Xu0kuW+Zp/KWUquPgeOzKbfo0YmKMvowWLWD2bFdHk2b8e/dfxm4cy8zQmXh5eNG7am8+qPkB2TNkt+r4deuMR1anT0PPnjB+PGTWdZ21p9DzNP6zUUT66eG1SbB9u1G/wspHU5p9PJv5Waa9NI1jPY/xRuk3mLhlIr5Tfflk8yfcj376GNv69Y2+jl69jDWzypY1ptloWlpmS9IoDbwHXBSRlSIyTkTcdja4WwkJAU9PaNjQ6kMs1a9J7bVt7M03uy/zm89nX/d91CxYkyFrh1B8anFm7p5JdKzldWIzZ4YpU4xq9unSGf+E3brpGlZa2mXLjPA3lFIlgaIY/RkngWqOCixVCQ6G6tWNpeaslC1jOrJlTGdzm5a4cs+W43+t/8fGDhvxze5L95XdKf11aZYcWkKcslzLq1YtY17HoEHGE0Z/f+OfVdMeV7FwdioWtu4RaEpkS59GDqAfkAc4AvyglLrhwNiSzS36NP79F/LmhXHjYOhQ18aiPaSUYuXJlQxZO4RDlw9RMV9FJjSYQCPfRjztCezOnca8jsOHoV07mDwZcuRwUuCa5gT26tNYjLEu+AogI7BZRKrYIb7ULb6kqu7PcCsiwislX2Fft33Mbz6f6w+u02RBExr80IAd5y1P0qhSBUJD4aOPYOFCKFMGli93Ttya5mq2JI3cSqmJSqn/KaUmAK8CUx0UV+oRHGws6BAQYNNhA5fuZ+DS/Ta3abbx9PCkTbk2HOtxjKlNp3Lo8iGqzanG60te5+iVo4kelz49jBljVM7NmxeaN4e33oLLl50YvOaWus8Ppft86yeXpjS2JI3rIlI2/hulVDjGHYeWmNhY406jaVPwsG1eRf5nfMj/jPlaFpbatKRJ75WeXlV7EdY7jNF1R7MmfA3+M/zp9Fsnzt46m+hxFSoYj6vGjoVffwU/P/j+e0jiUudaKnDjfhQ37ke5OgyHsaVP43mM5V43AQcxRlM9p5Rq7rjwksflfRo7dkC1avDjj8Zi1VqKcfX+VSZsmsD0XdMRhB6VezDkhSHkypgr0WOOHDEmBW7daoyymjkTfH2dGLTmFvQ8DROl1DGgIrAeozN8P6A/CS0JCTHWzWjc2NWRaDbKlTEXnzf5nJO9TtK6bGu+3PElvlN8GfPXGO5G3TV7TOnSxuzxr74yfl/w94fPPzfmdmpaamHNehojTX/WBNIrpX5SSo1USn0bv56GlojgYKPXNGdOmw+1VL8mtde2cSeFninE3KC5HHzvIA19GzJiwwiKTS3GtB3TiIp98hGEhwe8/74xsqphQ2PBp2rVjOG6mpYaWHOnYaq0Rx9gh4icFJHfRGSMntxnwbVrxsPuJFa19c2dGd/c5mtWWGrTHKN07tL88uYvbHt3G6Vzl6Z3SG9KTS/FggMLiI2LfWL/ggXht9/gp5/g3DkIDIQPP9QFENOCmsVzUbN44o8xUzqb1ghXSp02VbYtA5QF/JVSgxwZYHK4tE9j0SJo3dooIVK1qmti0BxCKcWqsFUMWTuEvZf2UjZPWcY3GM/LJV42O8fj+nVjUuDcuVC8uLHUbL16Lghc06xkr3kavwAopSKVUnuUUt8DP9sjwFQpJMR4LBVo9u9dS8FEhCbFm7C7624Wv7GYBzEPeHXRq9T+rjabz25+Yv8cOWDOHFizxhhVVb++UQjxhltPjdU086zp02glIp8AWUTET0QSHjPLXoGIyAARUSLyxH2diASIyDYROSwiB0TkTXtd1yHi4oz1Mxo3NmpOJUHPH/fQ88c9NrdpzuMhHrzp/yZH3j/CNy9/Q9j1MF6Y9wKvLnqVA/8+WRusQQOjAOLgwfDdd8bw3GXL9PDc1Kb93J20n7vT1WE4jDV3GlswyoZkBz4HTonIHhH5H2CXJ7QiUhBoDCQ2IP4+0E4pVQZoCnwpItnscW2H2LfPKB+SjFngltYZTu1rEKc03p7edAvsxqnep5jQYAKbzmwi4JsA2v7alr9v/P3IvhkzwqefGpMCn3sOWraEZs3g/HnXxK7ZX0R0LBHRT/ZzpRpKKatewAtAZtPXOYFKQCZrj3/KuZcB5YHTQC4r9t8PlHjafpUqVVIuMW6cUqDUpUuuub7mUtfuX1MfrP5A+Yz1Ud6jvVWvP3qpS3ee/L8QHa3UpElKZcigVJYsSk2bplRMjAsC1uyq1TdbVatvtro6jGQBdqtEPldt6dOYCcSYEs014AyQ7O48EQkC/lFKWVUXw1TvKh0Qlkh7VxHZLSK7r1y5ktzwkiY42Fho+tlnXXN9zaVyZMjBJw0/4VSvU3QM6MjXu76m2NRijFg/gtuRtx/u5+VlDMk9eNAYlturl1EMea8eTa25MVuSRoRKMC9DKXUVGG3NgSKyRkQOmXkFAUMxSq1bc558wHygo1Lma1krpWYppQKVUoG5c+e25rT2dfMmbNuW5KG28SzVr0nttW1Si+eyPsfMV2dypMcRXi75MmM2jsF3ii+Tt00mIua/KU7FihldYD/+CGfPGmMn+veHu+bnEGqaS9mSNMJF5PFPQqsWdVBKNVRK+T/+AsIx1ufYLyKngQLAHhHJ+/g5RCQrsBIYppTabkPczrVmjVFzKplVbSsWzkbFwtlsbtPcT8mcJVnSYgm7u+ymUv5K9F/Vn5LTSjJv7zxi4ozp4iJGpZmjR41SJJMnGx3lv/3m4uA1mzXwy0MDvzyuDsNhbJmnkQ8IxugU344xVyOvUirIbsEYiSPQdBeTcHs607VXKKW+tPZ8Lpmn0bmzMSTm6lXj+YOmPWbd3+v4cM2H7LqwC79cfoyrP45mzzd7ZI7Htm3GCoEHD0JQEEybZkwY1DRnsFftqYsYnd8/A7mBA0Bru0RohogEishs07etgNpABxHZZ3oFOOraSaaUMT+jUSOdMLRE1S9anx2dd/Bzq5+JU3G8/tPrVJ9TnQ2nNzzcp3p1Y82OiRONQsl+fsbdh65jpbma1XcaKZHT7zQOHoRy5Yz1QN99N1mn6vz9LgBmt69sU5uWssTExfD9vu8Z+ddIzt8+T5NiTRjfYDwV81V8uM/p09CjB/zxh7Esy6xZUFn/07stXeX2v5PUF5E5IvK5iHQUkUqmkiJavJAQ4087rNJXo1guahQzX7/GUpuWsnh5ePFuxXc50fMEnzX6jF0XdlFpViXeWvYWp66fAqBIEfjf/2DpUmORp6pVjZFWt265NnYtbbKlT+M00BfwBsqZXmWUUsUdFVxyOf1Oo359oy/jwJOzgTXNGrcibvHZ1s/4YvsXRMVG0blCZ4bXGU7+LPkBuH3bWGZ2+nRjxcCpU+GNN4yOdM096DuN/5xRSi1XSi1VSg1XSgW5c8Jwujt3YPPmZA+11dK2Z3yeYUz9MYT1DqNbpW7M3jub4lOLM2TNEG48uEHWrEai2LHDSBotW8IrrxiPsDTNGWxJGhtFpJ+YK+Opwbp1EB1tl0dTYLl+TWqvbaNB3sx5mf7SdI73PM7rfq/z6ZZP8Z3qy6ebP+V+9H0qVzYq70+eDH/9ZSwANX48REa6OnIttbMlaZQG3gMuishKERmn19NIICQEMmeGmjXtcrqGfnlomMhYb0ttWurim92XBa8vYG+3vdQoWIMP135IiWklmBU6CyXR9O1rzO146SUYNswYh7FmjaujTtteKZePV8rlc3UYDmPz6CkRyYCRQMoCZZVSAxwRmD04rU9DKWMx6PLlYflyx19PS7M2ntnIkLVD2HpuKyVylGBs/bG0KN0CD/EgJAR69oSwMGjVCr74wiiKqGm2stfoqd0iMgfoBmQDVrpzwnCq48eNh8p2ejSlaYmpXbg2mztu5ve3fie9V3reXPYmlb+tzKqwVTRpojh0CEaNMmaSP/+8sUZ5dLSro05bHkTF8iAq9Va5teXx1GvAUozSId2A06YRVZodh9rGe2f2dt6Zbb5aiqU2LfUTEV4t9Sr7uu3jh2Y/cO3+NZosaEKDHxpw4NpORoyAI0egTh2jIGLFirBpk6ujTjs6zNtJh3mpt8/RlhnhF5RSIUqpiUqpVkAgMMdxoaUgwcHGr3VFitjtlK+Uy88r5fLb3KalHZ4enrQt35bjPY8zpekUDl0+RNXZVXnjpzeIzHKUFSuMp6W3b0Pt2tC+vbHMi6Ylhy2Ppwon/F4pdRQoafeIUpr7943hK3Yeavt2lUK8XaWQzW1a2pPeKz29q/YmrHcYo+qOYlXYKvxn+NP593epWPccR47AkCHGsvWlSsFXXxk1NTUtKaxZ7nWK6ctfReS8iGwSka9FZDzg79jwUoANG4xxjro/Q3OxLOmzMKLOCMJ7h9Onah8WHFxAiWkl+HjLQAZ8dI0DB6BSJaOzvEoVY66HptnKmjuN2vFfKKUKAB2B1RiLIL3qoLhSjpAQyJDBuP+3ozdnbns4s9SWNk3LnSk3XzT5ghM9T/B22beZvH0yvlN9WXZ5LMv/uMuiRXDxolEUsVs3uHbN1RFrKYk1SWOtiGwDnhWRTsAzwB9KqTlKKb2ycUgI1KsHPj52PW2LSgVoUamAzW2aFq9wtsLMC5rHge4HqF+0PsPXD6f4tGJc9Z3OgcNR9O0Lc+YYj6xmz4Y4s8uaabZK7T+fVs3TEJFiwHrge4z5GWWAKOCQUupNh0aYDA6fpxEWBsWLG3UdevVy3HU0zQ62ndvGh2s/ZOOZjRTNVpQx9cZQRr1Nr54ebN5sVM6dNs0oiKilbcmep6GUCgMammpONVNKlQCqApPtGGfKEz/U1gH1pqJj44iONf+rn6U2TUtM9YLV2dB+A8HvBPOMzzO0+bUN7bdV4IPZK/nhB8W5c8Za5Z066VFWyXH9XhTX70W5OgyHsWXI7YnHvr/r1suuOkNwsLHAc3H7121sM3sHbWab76m01KZplogITYs3JbRrKIveWMS9qHu8uvgVvo2pw/frtjBwIMyfDyVLwpdf6omBSfHeglDeWxDq6jAcxpbJfVpCERGwfr3Dqtq+VaUgb1Uxv76npTZNs4aHePCW/1sc7XGUGS/P4OT1kzT5qRbHK77Gz5sOUq0a9OsHFSoYtTg1LZ5OGkm1aZMxR8NBQ22bVyhA8wrmO9MstWmaLbw9veke2J1TvU4xvv54Np7ZSLM/y/Ns93Z8s/hv7t2DBg2MWlZnz7o6Ws0d6KSRVCEhkD491K3rkNNbql+T2mvbaM6XKV0mhrwwhPA+4QyqMYilR5bS60QpXprWh0EjL7NihVH0YOxY4yZbS7t00kiqkBBjbkamTA45vaX6Nam9to3mOjky5ODTRp9ystdJOgR0YOaer/ja25euiz6m0Su3GT4cypSBFSuM4s5a2qOTRlKcPWtUhHPgLPA21QrTplphm9s0zR4KZC3ArFdncfj9w7xU4iWm7h/N1srFeO/7L/HOEMFrr8HLL8OJE08/V1qT2n8+bV5PIyVx2DyNWbOMqbRHjoCfn/3Pr2luZveF3QxdO5TV4asplLUQ1aNGsXJ8WyIfeNK/v7FueebMro5Ssxd7rRGuxQsOhkKFjIe8DnI7IprbEebHO1pq0zRHCMwfyKq2q1jTdg15MudhSURH8o8pR63Oy/n0U0WpUrBwoX5kBXDh5gMu3Hzg6jAcRicNW0VFwdq1xlBbBy6X3uX73XT53vxdkqU2TXOkBr4N2Nl5J8taLkM8Yln/bHPKTK5BZv+/aNMGatQw1i5Py/ot2Ue/JftcHYbD6KRhq61b4c4dh1e17VizCB1rFrG5TdMcTUR4o/QbHHr/ELNfnc3NuHOcqFEX/09e5MSdvVStaqzdceGCqyPVHEEnDVuFhICXlzF43YGa+uejqb/5xekttWmas3h5ePFuxXc52eskkxpN4h92cL1lRfw+eptFf56iZEkYP14P0U1tdNKwVUgI1KoFWbI49DKW6tek9to2WsqSwTsDA2sMJLxPOMNeGMaZDL+j3vcjd8f3GDbhIn5+sGyZ7u9ILXTSsMWFC7B/v1MWXLJUvya117bRUqZsPtkYW38sYb3D6FqpK+fzzCb94GLcqzaUlm1vUq8e7Nvn6ii15NJJwxZ//mn86aB6Uwl1ecGXLi/42tymaa6WN3Nevnr5K471OMYbZZpz5fkJZBziy650E6lQ5QHdusHly66O0nFS+8+nnqdhi1atYMsWOH/eoSOnNC012XdpH8PWDeOPk3+QKS4/D4JHkulERz4e7kWvXpAunasj1B6n52nYQ0wMrF5tPJpyQsK4fCeCy3fM9yBaatM0dxOQN4CVrVfyV4e/KFeoMHEvdyWuexkGzl1KGf84/ve/1NXfEXblLmFX7ro6DIfRScNaO3bAzZtO6c8A6PXjXnr9uNfmNk1zV7UL12ZLpy389tZvFCnkDa1aca5pFV7ts5oXXzQKLKQGQ385yNBfDro6DIfRScNaISHg6QmNGjnlcu/VLcZ7dYvZ3KZp7kxEeK3Ua+zvvp/vm31PXt+r0K4xa55rSNkmu+jdG65dc3WUmiU6aVgrJMRYCzNbNqdcrm6pPNQtlcfmNk1LCTw9PGlXvh3Hex7nyyZfkq3kfuI6V2Havy0oWvkYn38OkZGujlIzx22ShogMEBElIrks7JNVRM6LyHRnxsbly7B7t9MeTYHl+jWpvbaNlnak90pPn2p9+LtvOCPrjCRjuT+5064MAzd2pmTgeX7+OXX1d6QGbpE0RKQg0Bh42tpgY4CNjo/oMatWGX86YahtPEv1a1J7bRst7cmSPgsf1/2Y033D6VutN96V5nOueXFazBxE9QbX2LXL1RFq8dwiaQCTgcFAor9TiEgl4FlglbOCeig4GPLkMRZMdpJe9UvQq34Jm9s0LSXLnSk3k5tO5kTv47Sp8BZS83N2VvOlyoBxvNX2XopYcja1/3y6fJ6GiAQB9ZVSfUTkNBColLr62D4ewDqgDdDQtE/PRM7XFegKUKhQoUpnzpxJXoCxsfDss/DSS/DDD8k7l6ZpNjl0+RAfrvqIlWG/wd1n8do6nAF1uzDsw3SOruSTprl8noaIrBGRQ2ZeQcBQYMRTTvE+8IdS6vzTrqWUmqWUClRKBebOnTv5wYeGGsM5nNifAXD22n3OXrtvc5umpSb+efz5X5vlbOm0hSrFShLTuCef3vWjwIs/8s3MOGJiXB3hkw5fuMXhC7dcHYbDOCVpKKUaKqX8H38B4UBRYL/pLqMAsEdE8j52iupAT9M+nwHtROQTZ8ROSIgxma9xY6dcLt6gZfsZtGy/zW2alhrVKFiD7d3+4o/Wf1CicBZuN3qH9/ZVpPiLf/Dnn+7VUz56xRFGr0glk07McPnjqYQSezz12D4dsPB4KiG7lBGpXh3i4ozJfU60PdwYrF7NN6dNbZqW2sWpOBYfWkL/FR/xb3Q4nK5N1bsTmDOyBmXKuDo6eHPmNgCWdKvu4kiSzuWPp5JCRAJFZLZLg7h2zViGzImjpuJV882ZaFKw1KZpqZ2HeNC67Nuc/eAoUxp/RZaix9nhXxP/cUG06nmIf/91dYSpm1slDaVUkfi7DKXUbqVUZzP7fGfNXYZdrFlj3GU4uT8DLNevSe21bTTNGuk809G7+vtc/DCMoVXHka7EBpbmKkeBXu0ZPP40D/RUJodwq6ThdoKDIUcOqFzZ6Ze2VL8mtde20TRbZEqXiXFNh3Lhg3A6+Q0gzm8Jkx6UIk+7vkybe5nYWFdHmLq4VZ+GvSWrTyMuDvLnh3r1YNEi+wZmhdAz1wGoVDiHTW2altadv32e7otHsfLCXIjKSJ7wAcxo25/XX87qlOunhp9PS30aOmkkZu9eqFgRvvsO2re3a1yapjne0SvH6Dh/ODvuLIN7uSj57zB+6PUeVQPTuzo0t5ciO8JdLiTE+LNJE5dc/vilOxy/dMfmNk3TDH65n2d7/6Vsbr+Tks+U54RvP6otLEnNHt8T/rfjnlmFnrn+8G4jNdJJIzHBwUbZkLyPTxlxjhG/HWLEb4dsbtM07VE1i1Tm+LA1/Np8NXmz5mZrng4U/7wczYf8xo0b9n/SMjHkOBNDjtv9vO5CJw1zbt2CrVtdMtQ23tCX/Bj6kp/NbZqmmdesXEMujNzF1/WWkuWZGJb7NCPPkJr0nLRRl2G3gU4a5qxda9SccsFQ23jlC2ajfMFsNrdpmpY4EeG92i24NuYwwwO+xTPHWb66X4fsPV9iwrx9xMW5OkL3p5OGOcHBkDWrseiSi1iqX5Paa9tomqN5eXgxOqgzN0adpHOhiUTl2c7QsxXI1a01C/445erw3JpOGo9TyugEb9QIvL1dFoal+jWpvbaNpjlLBu8MfNtxEJeHhfNqtqHczLucttv9KNzjfTbsvujq8NySHnL7uEOHoGxZ+PZb6PzEhHSnib+TKJP/GZvaNE1LuvArF3ln5hi2R30LsekoH9GXH3sMorRvNqvPkRp+PvU8DVt89hkMGgTnzkGBAo4JTNM0t7Yr7BRt547geLpF8CA79byHsLBPT/LlyuDq0JxCz9OwRXAw+Pu7PGHsP3eT/edu2tymaVryVS5WnGPjfmTFy3vJF1uN9V6Dee6TErScMJu79y0v4rH55FU2n0y0UHeKp5NGQnfvwqZNLh1qG2/8H0cZ/8dRm9s0TbOfVwIDuDDpD76tsYGsFGRZVBeyD/Pn/enLiI42/5Rm2rqTTFt30smROo9OGgmtXw/R0S4dahtvdJA/o4P8bW7TNM3+Ojeqw41JWxlTejnenp7MuNaSrAOrMGr+GlLxE36zdNJIKDgYMmWCmjVdHQml8mahVF7ziyBbatM0zTFEhI9aBnH70wP0fO47YtJfZmR4I7L3bsSs/yVzsbcURCeNeEoZSaNBA0jv+oJmlurXpPbaNprmzrw8PZnWuT03Rh+nVdbJ3Mm0j26hlcnXpyW/b0m95UPi6aQR78QJOH3aLR5NgeX6Nam9to2mpQSZfXxY0q8v/3wQRgPPj7mUOYSgVWX4Y9fXXL1y09XhOYwechtvyhTo2xfCw6FoUYfGZY34lfmK5c5sU5umaa5x4p/LtJ4xnlCPryHOk0qxvfjx/Q8pWTDlrauh52lYo2lT407j2DGHxqRpWuq27ehp2s/7mJMZ50NkVhr4DGZh7z48mz2Tq0Ozmp6n8TQPHsBff7nFUNt428OvsT38ms1tmqa51j2Vga87TOS3pgfIF1mHtQwj/4TivPnZDO49iHZ1eMmmkwYYCSMiwm36MwAmrz7B5NUnbG7TNM21vt0Uzrebwnmtmj8XvviNWVW3kCW6BD/de59sH/nRY8YiomNSbjldnTTAGDXl4wO1a7s6kocmtSjPpBblbW7TNM29dGlag+uf/8WoUivxUpn4+nJrsg6qyPAfgomLS3ndAzppgFHVtl49yOA+dWUK5cxIoZwZbW7TNM39eHgII956iTuT9tIr/0JiPe8w9u+XyNa3LlN/2ZaiJgjqpBEebgy3daNHU2C5fk1qr22jaamVl6cHU7u05ua4o7yTfTr3Mhynz8Ea5OndjIWrD7s6PKvopBESYvzpRp3gYLl+TWqvbaNpqV3G9OlY0LsHlz86xcsZxnIty3rabClLwV4dCNl+xtXhWaSH3L76Khw5AqdOgYhzArPChZsPAMif7clHZpbaNE1zraT8fJ69eo02Mz5hU9Q0QFHy9vt8/+5QqpXL7aAoLdNDbhMTGQnr1hl3GW6UMMD4D5fYfzpLbZqmuVZSfj4L5crJxuGTONj1JOU92nIi21SqL/alUv9RHA2/46BIkyZtJ43Nm+H+fbfrzwDYcPwyG45ftrlN0zTXWrH/Aiv2X0jSsf4FC7Jv1Gw2vX2Y4tKEPc+MpPRMX+p+MIVzFyPtHGnSpO2kERwM6dIZI6fczIwNYczYEGZzm6ZprrVg+xkWJLNfotbzz3Ny3DJ+e2UH+TzK8VfGvhSeVIrXhv/Ateuxdoo0adJ2n4a/P+TLB6tXOy8oK12+EwFAniw+NrVpmuZab87cBsCSbtXtcj6lFPM2rmFg8IfcyLAHj6tleCfveL7u8yqZMzvmsbru0zDn3Dk4fNgtH02BkRASSwqW2jRNS11EhE51GnF1wi4mBv5EhsxRzI8JIsfAWvT5fBMREc6NJ+0mDTcdahtvzZF/WXPkX5vbNE1LnTzEg0Evt+TG2MN84DcLyXGaqXdrk63Hywybup+oKCfF4ZzLuKHgYChYEPz8XB2JWfH1a2xt0zQtdfP29OaTVl24MeokXYp8Sky+rYy/EUD2zu8wYWYY0Q6uiZg2+zSioyFnTnj7bZg50/mBWeH6PePXhhyZ0tnUpmmaazn75/P6/Ru8v3ASS899SRzRZD3VlbGNhvN+u7x4eibtnCmiT0NEBoiIEpFcibQXEpFVInJURI6ISJEkX2zbNrhzx237M8D4D5fYfzpLbZqmuZazfz5zZMzO4i7jOTcwjBef7cLtErPofaoYed76iLkLbxFn54K6bpE0RKQg0Bg4a2G3H4BJSik/oAqQ9IkKwcHg5WWsB+6mQg5dJOTQRZvbNE1zraW7z7F09zmnXzd/1nz80eNrjvc6Ss3cr3HdfxzvHvDluZafsWjZA7sVRXSLpAFMBgYDZt+WiJQGvJRSqwGUUneVUveTfLWQEKhZE7JmTfIpHG3eltPM23La5jZN01xrWeh5loWed9n1S+Yqzua+i9jdeQ/lc1fhUrlBtN5akiLN57D895hkJw+XJw0RCQL+UUrtt7BbSeCmiPwiIntFZJKImH1aJyJdRWS3iOy+cuXKkztcvAj79rn1oymAb9sH8m17s48ULbZpmqYBVHquAvsGBrOmzXqK5XmOsxU603y1PyWb/UxIiEpy8nBK0hCRNSJyyMwrCBgKjHjKKbyAF4CBQGXAF+hgbkel1CylVKBSKjB3bjPFvv780/jTTYfaxsvq401WH2+b2zRN0xJqUKwuJz/YxtI3fiVfXg9OVWzBi79UpWzQWtavt/18TkkaSqmGSin/x19AOFAU2C8ip4ECwB4RyfvYKc4D+5RS4UqpGGA5UDFJwQQHG7PAy5VL6ttxCkv1a5JT20bTtLRHRGjh34xzQw7y7cvzyF7wEocrNaT+vMZUfi2ULVusP5dLH08ppQ4qpfIopYoopYpgJIeKSqlLj+26C8gmIvG3DvWBIzZfMCbGKBnStKnbVbV9nKX6NfaobaNpWtrj6eFJ58AOXBhygon1J5Op+B52Vwqk1pRWvNDsBDt3Pv0cbjVPw3S3EaiUuioigUB3pVRnU1sj4HNAgFCgq1LK4hzIJ+ZpbNsGNWrAkiXQqpWD3oV9PIgyipJlSPdk142lNk3TXCsl/XzejrzNJ399zufbPicqLgL2vEsDrxGsXV4g0XkabpU07O2JpDFiBIwbB1euQI4crgtM0zTNjVy+d5mRa8cxa+8MYmM8YWyE+0/uc4qQEKhWLUUkjF/3nufXveaH7Vlq0zTNteZvO838baddHYZN8mTKw9evTeFk7+O8VbalxX3TTtK4cgV273b7obbxFu88x+Kd5icIWWrTNM21/nfgIv87kDIn3xbNXpRFrX6wuI+Xk2JxvVWrQCm3H2obb0Hnqklq0zRNc6S0kzSCgyF3bqiYtJG6zubtmfhNoKU2TdM0R0obnz5xccakviZNwCNlvGVL9WtcVdtG0zQtZXyCJteePXD1aorpzwDL9WtcXdtG07S0K1UPuRWRK0BqnwWXC7jq6iCcTL/ntEG/Z9cprJQyU4cplSeNtEBEdic2njq10u85bdDv2T2ljcdTmqZpml3opKFpmqZZTSeNlG+WqwNwAf2e0wb9nt2Q7tPQNE3TrKbvNDRN0zSr6aShaZqmWU0njRRKRFqKyGERiTOtPRK/vZGIhIrIQdOf9V0Zpz1ZeM85RWS9iNwVkemujNHeEnvPprYhInJKRI6LSBNXxehIIlJeRLaZ/j+vEJGsro7J0UQkQES2i8g+EdktIlVcHVNCOmmkXIeA14GNj22/CryqlCoLtAfmOzswB0rsPUcAwzHWkE9tzL5nESkNvAWUAZoCX4uI+6/6Y7vZwIem/8+/AoNcHI8zTARGKaUCgBGm792GThoplFLqqFLquJnte5VS8QuIHwYyiEh650bnGBbe8z2l1GaM5JGqJPaegSBgsVIqUin1N3AKcKvfSO2kJP8lzNXAGy6MxVkUEH9H9QxwwcK+Tpd2qtymTW8Ae5RSka4ORLO754DtCb4/b9qW2hzGSJDLgZZAQZdG4xx9gT9F5DOMX+xruDacR+mk4cZEZA2Q10zTMKXUb085tgzwKdDYEbE5SnLec0qVFt9zQpbeP9AJmCoiw4HfgShnxuYoT3nPDYB+SqmfRaQVMAdo6Mz4LNFJw40ppZL0H0VECmA8/22nlAqzb1SOldT3nJIl8T3/w6O/dRcwbUtxrHj/jQFEpCTwsuMjcjxL71lEfgD6mL5ditGv4zZ0n0YqIyLZgJUYnYdbXByO5ji/A2+JSHoRKQqUAHa6OCa7E5E8pj89gI+Ab1wbkVNcAOqYvq4PnHRhLE/QM8JTKBFpDkwDcgM3gX1KqSYi8hEwhEf/ozVWSl12fpT2ldh7NrWdxug8TGdqa6yUOuKSQO3oKe85/vFNDNBXKRXsqjgdRUT6AD1M3/4CDFGp/ENLRGoBUzCeBEUA7yulQl0b1X900tA0TdOsph9PaZqmaVbTSUPTNE2zmk4amqZpmtV00tA0TdOsppOGpmmaZjWdNLQ0Q0R6i8hREVno6ljsSUS+FJHaFtrbi8iix7blEpErpnkei0WkhOMj1VIDPeRWSzNE5BjQUCl1/rHtXkqpGBeFlSwikhNYqZSqZmGfrEA4UEgpdd+0rTtQRSnVSUTqAG2UUl2cErSWouk7DS1NEJFvAF8gWET6ichIEZkvIluA+SKSW0R+FpFdpldN03E5RWSVaU2L2SJyxvRbehEROZTg/ANFZKTp62IiEmJaz2STiDxv2v6diEwVka0iEi4iLRIc/4FpzYj9IvKJ6Rx7ErSXSPh9Am8AIQn2qyQif5mu/aeI5FNK3Qb+Al5NcNxbQPzdxyagoYjoskLaU+mkoaUJSqnuGOUZ6imlJps2l8a483gbYwbuZKVUZYwP4vh6Px8Dm5VSZTDqeRWy4nKzgF5KqUoYa3x8naAtH1ALeAX4BEBEXsSo5FpVKVUemGiqGXZLRAJMx3UE5pm5Vk0g1HQeb4zZ4y1M154LjDPttwgjUSAi+TFKjq8z/d3EYZRWL2/Fe9PSOP2bhZaW/a6UemD6uiFQWkTi27KKSGagNsYiSCilVorIDUsnNB1TA1ia4FwJ1zNZbvqQPiIizya49rz4R0dKqeum7bOBjiLSH3gT8+tl5AOumL4uBfgDq03X9gQumtpWYizUlBVoBfyslIpNcJ7LQH5MCUjTEqOThpaW3UvwtQdQTSn1yEJOCT74HxfDo3fqPgnOc9O06po5Cdc2SfTkJj9j3OmsA0KVUtfM7PMgwbUFOKyUqv74TkqpByISAjTHuOPo/9guPqZzaZpF+vGUphlWAb3iv0nwWGgj0Nq07UUgu2n7v0AeU59HeozHTZj6D/4WkZamY0REnvbYZzXGHUVG0zE5TOeKAP4EZmD+0RTAUaC46evjQG4RqW46j7dpXZV4izCSxbPAtsfOUxJjaVlNs0gnDU0z9AYCReSAiBwBupu2jwJqi8hhjMdUZwGUUtHAaIxy5KuBYwnO9Q7wrojs57+V5xKllArBKHW+W0T28eha5wuBOIykZs5KoK7pPFFAC+BT07X38eiqb6sxHkEtSVgp1vSY7IFS6pKlODUN9JBbTbOJqQR7oFLqqpOuNxB4Rik13MI+m4FXlFI3k3iNfsBtpdScpEWppSW6T0PT3JSI/AoUw1iIx5IBGKO6bibxUjeB+Uk8Vktj9J2GpmmaZjXdp6FpmqZZTScNTdM0zWo6aWiapmlW00lD0zRNs5pOGpqmaZrV/g+QXjrrCxHgGgAAAABJRU5ErkJggg==)
%%%% Output: stream
band = 5
-0.46551944347590346 -0.34895815035534355
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZEAAAEWCAYAAACnlKo3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABRD0lEQVR4nO3dd3yN5//H8dcny1YratUsgiBGUWqrakuNDqVq/uxSq7qt0pYvVXtUjarVUkrVqt2aQYwQtWvvLYmI6/fHfZBEIkNOTs7J5/l4nEfOfd3jvE/Euc59Xfd9XWKMQSmllEoIN0cHUEop5by0ElFKKZVgWokopZRKMK1ElFJKJZhWIkoppRJMKxGllFIJppWIUkqpBNNKRDkdETkuInUiLL8rIldFpLpt2U9EdojIHdtPvwjbiogMFZHLtsdQEZEI6+21b00RWSsi10XkeCzvL7+IGBG5ZXucF5HxIuIZYZssIrJQRG6LyAkRaR7lGM1t5bdFZJGIZLH3viJSQ0TuR8h9S0RaPem9KuenlYhyarYPqXHA68aY9SLiBfwO/AxkBmYAv9vKAToAjYDSQCmgAdDRdix77nsbmAp8FI+3l8kYkx4oCbwIdI2wbhxwF3gWeA+YICIlbFlKAJOA923r7wDj7b2vzRljTPoIjxnxeL/KGRlj9KEPp3oAx4E6WB/gl4DyEdbVBU4DEqHsP6Ce7fkmoEOEde2ALfbeN0JZHeB4LO8vP2AAjwhlw4DJtufpsD7Ii0RYPxP41vb8a2B2hHWFbNtnsPO+NYBTjv770EfSPvRMRDmrzsAgoLYxxj9CeQlgj7F9qtnssZU/WL87wrrdUdbZa98EE5FcwCvAFltREeCeMebfJ2R5mNMYcwTbh7+d9wXIbmt+OyYiI0UkXQLesnIiWokoZ/Uy1ofq3ijl6YHrUcquY32Tjm79dSC9rW/DnvsmxCURuYZ1hnMbmB8hx4145Iy43p77BgF+QE6gFlAO+C6G96ZchFYiyll1xvpmPCVi5zZwC8gYZduMwM0Y1mcEbtnOIOy5b0JkM8ZkAtIC/wArYsgRW86I6+22rzHmnDFmvzHmvjHmGNAXePPJb1E5O61ElLM6D9QGqhK54zcQKBWlYillK3+wvnSEdaWjrLPXvglmjAkGpgOVRCQb8C/gISKFn5DlYU4RKQiksu1nz30fi45+xrg+R3fK6EMf8X1g61i3Pc8LHANG2pa9gBPAh1gffh/Ylr1s6zsBB4DcQC6sD8BOSbCvG5AaeNVWnvrBumjeX34idKzbjvctcBZbxz0wF5iD1dldBatZqYRtXQmsZqeqtvU/A3MjHN9e+9YE8gECPAesBaY5+u9FH3b+/+joAPrQR3wfESsR23IB4CTwjW25DLADCAZ2AmUibCtYVzpdsT2GEfmKKnvtW8NWMUR8rIvh/T2oRG7ZHteA9cALEbbJAizC6iv5D2ge5RjNbeW3sS49zmLvfYFeWP03d2z/HqOBDI7+e9GHfR8PvtUopZRS8abtlUoppRJMKxGllFIJppWIUkqpBNNKRCmlVIJ5ODpAUsuWLZvJnz+/o2MkyNGLtwEo6K0jSSjl9O7fh127IHduyJHD0WlitWPHjkvGGO+o5Snu6qw0uTOZfXt2UyhrPkdHibemkzYDMK/jiw5OopSKj13/XQWgTN7MjwqvXIGsWWHUKOje3UHJ4k5EdhhjykctT3HNWSH3b1L4++J0nfU/wsLDHB1HKZUCtJy6jZZTt0UuDAmxfqZOnfSBElGKa87Km7YEF07nZ7xHX+b0m8nMphN5vVRlR8eKk261Cse+kVIq2ela4/nHC4ODrZ9OXomkuDMR78xeXBqzmCahi7gafI36C6tQdVgHLt2+4uhosXqpcDZeKpzN0TGUUvHUqUYhOtUoFLnwwZlImjRJHygRpbhKBCBdOljwdUO2td5Pnv968/ftqeQc4sOg338iOfcRBZ65TuCZqKN0K6WSuy1HLrPlyOXIhS7SnJUiK5EHXiiVnv+mDGdIvh3I1UL0D2jFc1/UYtvRIEdHi9agJfsZtGS/o2MopeKp/Ux/2s/0j1yozVmuQQQ+a1uac1//Q42bkzgdHkDF6aVoPLofd+4GOzqeUsoF9KxThJ51ikQu1OYs15Ilsxtrh3dgxRtBZD3blEVXvyJrv5JMWbvS0dGUUk6u7UsFaPtSgciF2pzlmupWfpbzE2fSI/Nf3A1xo/2GVyje712OXDjr6GhKKSe17uAF1h28ELlQm7Ncl7s7jOxem2N991D66gAOmEUUHuVDxx/HEX4/3NHxlFJOptucXXSbsytyoZ6JuL68uVIT8H1/Zr20l7RXKzD51Adk++RF/tix0yF5+tYrSt96RR3y2kqphPvkVR8+edUncqH2iaQczV8pzJXvV9I81Wyu8x8NFr9Ala96cOnmjSTNUS5fFsrly5Kkr6mUenrvVczHexWjDLWkzVkpi5eXMOuTZuzpEETBK53YFD6anIOLMeCX+Ul2b8mOE1fYcSL53xSplIpsxb5zrNh3LnKhNmelTL7PZ+LImHGMKL4Zt+DsDDzwNnk+rs+2f4/Z/bWHLT/IsOUH7f46SqnE1Wf+bvrM3x25UCuRlK3XOxW59M12Xg7/jjOe66n4UwkaDv+G4Lt3HR1NKZXM9K9fnP71i0cuDA4GT0/rSh4nppXIU8iQzoOVg3qy9s0DeF97lcW3PyPLZ2WYsnKjo6MppZKRt8o/x1vln4tcGBLi9GchoJVIoqhR9jnOj1nAR7n+4K65TfvN1Sj+aVuOnrvk6GhKqWRg8e7TLN59OnKhViIqIhEY1v51Tn4aSNk7H3PAcybPf+9Dx4nTuH8/+Q7qqJSyv88X7uPzhfsiF4aEOP3lvaCVSKLLlS0dO4Z+y9yau0gX7MPk823J2qc6S7cHPvWx+zUoTr8GxWPfUCmVrAxp7MuQxr6RC4OD9UxExaxpTV+ujNhA8/RTuO4VSP0lflTp/ymXb9xJ8DFL5HqGErmeScSUSqmk8Ebp3LxROnfkQm3OUrHx9HBjVu927OkYRMFbLdjk9i05virBV3P/TNDx/j50ib8PaT+LUs5mvv9J5vufjFyozVkqrnwLeHPku2mMKLkOt/DU9Dv4Os/1ehv/f0/HvnMEY9YcYsyaQ3ZKqZSyl4F/7GfgH1HmAtLmrLgTkakickFE9sWw3kdENotIqIj0ibLuuIjsFZEAEfGPUD5ARE7bygNE5DV7v4+n1atJdS4N3k0dGcKptH/wwgwfmgwbRcjde46OppSyo+FvlWb4W6UjF2pzVrxMB+o9Yf0VoDswPIb1NY0xfsaY8lHKR9rK/YwxCWsjSmIZ0nqxqt9nrHkrkGy3X2JhcA+yfFyR6Su3OzqaUspOXvHNwSu+OSIXaiUSd8aYDVgVRUzrLxhjtgNhSZEnOajpV5Dz3/1Jzzy/EOp5ljabKlLy4w84cV7nUFfK1czaeoJZW09ELtQ+kSRjgJUiskNEOkRZ94GI7LE1l2WO6QAi0kFE/EXE/+LFi/ZNGw9ubsJ37d7mxEdBlA79gH2pJ1BwhA/dJs3Ve0uUciHfLgvi22VBkQu1TyTJvGSMKQu8CnQVkWq28glAIcAPOAuMiOkAxpjJxpjyxpjy3t7e9s4bb3m8MxLwzWhmVt9K6rDcjD3XDO9e9Vi143Ck7b5uUpKvm5R0UEqlVEKNaVaGMc3KRC7U5qykYYw5bft5AVgIVLAtnzfGhBtj7gM/PCh3Zi1qlefqsK28lW40V9Jupu4iX2oN/Irrt0IBKOSdnkLe6R2cUikVXzWKZqdG0eyRC7U5y/5EJJ2IZHjwHKgL7LMt54ywaeMH5c7Oy9OdX/t0Y1e7IPIGN2Qt/fDuV5rh89fy1/7z/LX/vKMjKqXiaerfx5j6d5TpIlykOcsjKV5EROYANYBsInIK6A94AhhjJopIDsAfyAjcF5EeQHEgG7BQRB5knW2MWW477DAR8cPqMzkOdEyK95JU/Arl4sTweXzzaxv6b+3KR4G1SLvibV56/v+oU7yuo+MppeJh5F//AtD2pQJWwb171sMFKhFJqln5kovy5csbf3//2DdMRq7eDKbRd1+zIXwo3E1P8+xDmd69HZ4eyfpEUills+XIZQAqFcpqFdy+DenTw9Ch0LevA5PFnYjsiOY2i+TdnKUsmTOkYX3/r6jqPgfPG8WZfbMDWXpX5Zf1exwdTSkVB5UKZX1UgcCjWQ21T0QlpZw5ctGw5DA6Zp/B7TT/0nRNWcp/9hHnr952dDSl1BNMXHeEieuOPCoIDrZ+ukBzllYiTsbNTZjYuSWHugfhE9KGHamGk3tIcT6bsdjR0ZRSMRi37jDj1kW4ZN9F5lcHrUScysimfoxs6gdAoVxZOTD0Bya+8Dee9zPyzfGG5OzZiE2B/zk2pFLqMT+1rcBPbSPchaDNWcoRcmVKQ65Mkf/oOr5Whctf7+Q1z6GcS7eSKrOLU//rEdwJSTEjyCiV7JXJm5kyeSMMqqHNWcoRluw+w5LdZx4rT5vak6Wf9WXju/t59k5Nlob1Ieun5Zm8bLMDUiqlohq9+hCjV0eYxkGbs5Qj/LzlBD9vORHj+pd883NmxGL65vuNMI8rdNxahRIfd+L4uatJmFIpFdUPG4/yw8ajjwq0OUslV25uwtDWjfnvk/2UDe3J/tRTKDTSh66TftZBHZVykHkdKzGvY6VHBdqcpZK7XFkzsOObEcyu4U+a0AKMP/c+3r3rsHLnQUdHUyrFKZ7zGYrnfOZRgTZnKWfRrKYfV4dvomm6CVxJvYNXFpai5qD+3LgT4uhoSqUYw1ccZPiKCF/gtBJRzsTTw425fTpZgzreeot1ZhDZ+pVkxKJVjo6mVIowY/NxZmw+/qhA+0SUI0xoUY4JLcoleH+/53NwYuQsBhdZhbkv9Nldl4J9m7H/v3OJmFIpFdWirlVY1LXKowLtE1GOkCWdF1nSeT31cT5vVodz/fdQOaw/x1L9hu9EH1qNncC98PBESKmUiuqxuYC0OUs5wq/+J/nV/2SiHCvrM6n5Z/AAfq+3lww3y/PT5S5k6VuZRVt2JcrxlVKPfPPnfr75c/+jAq1ElCPM33GK+TtOJeox36hShCvfr6LtM7O45X6cxsvK8+LAnly6cTNRX0eplGz2tpPM3hbhC2BwMHh5gZvzfwQ7/ztQT83dXfixR3P2dw3i+esd2WJGkXNIMb6a/xspbb4ZpexheY+qLO9R9VGBi8yvDlqJqAh88mXm0Pfj+a7EJtxCstEv8E3yftKAXceOOzqaUk4td6a05M6U9lGBViLKlfV8uxIXh/hTK2wEpzzWUfbH4jQdPZS793RQR6USYsDiQAYsDnxUEBLiEpf3glYiKgYZ03uwenAvVr5xgCxXXuWXq5+Q5bMyzP7nb0dHU8rpLNh5igU7I/RnBge7zJmIzrHuRILvWpfgpvFyT9LXvX8fuo35gwknPsA8c4KytGPJh0PJlSlr7Dsrpbhy6y4AWdLbLtFv2BBOnICAAMeFiiedY90FpPFyT/IKBKwLSMZ9WJ9DPQLxudSXneEzyDusKJ/9Ml073pWKgyzpvR5VIKDNWcoxZm4+zsyIQycksUJ503FgzFDGld6Jx7WifHOgDbk+q8Gmw/tj31mpFOzzhXv5fOHeRwUu1JyllYgT+WPPWf7Yc9bRMejyZkkuDttIvdApnLu/jyo/+dFw9OfcvnvH0dGUSpYW7z7D4ogTyunVWfEjIlNF5IKI7IthvY+IbBaRUBHpE2XdcRHZKyIBIuIfoTyLiKwSkUO2n5kfP7Kylwzp3Vj2dTvWvBlEtrPNWXz1a7z7+zJ143JHR1Mq2dn6WW22flb7UYFWIvE2Haj3hPVXgO7A8BjW1zTG+EXp1PkEWG2MKQysti2rJFazgjfnJk6ne8a1hN7xot2aVyn11Tscv3La0dGUSjbSenmQ1svjUYH2icSPMWYDVkUR0/oLxpjtQHxuRGgIzLA9nwE0SnBA9VTc3WFUzxoc6r2b4ue+Ym/oEgp9V4yec0cTfl8HdVSq7/zd9J2/+1GB9okkKQOsFJEdItIhQvmzxpgHHQTngGeTPpqKqGDeVARO+IKJpffhdb4y3x/8kJz9KrLuX+e8pFqpxLJs3zmW7Ysw5YILNWd5xL6Jw71kjDktItmBVSISZDuzecgYY0QkxmtNbZVPB4C8efPaN60dzev4oqMjxEnHtwvR7JVlNBv8K3/e/5CasytQL2tX5vzfYDKleSb2AyjlYvYOeCVygTZnJR1jzGnbzwvAQqCCbdV5EckJYPt54QnHmGyMKW+MKe/t7W3vyArImFFYOuwd1r4ZRLZjXVl+eRw5Bhdj/Ppf9N4SpbQ5K2mISDoRyfDgOVAXeHCF12Kgle15K+D3pE+YtCZvOMLkDUccHSNearz4DGd/HMOHabdx90pOuq5rSvGvX+Pfi871PpR6Gj3nBdBzXoC1cO8ehIdrJRIfIjIH2AwUFZFTItJORDqJSCfb+hwicgroBXxh2yYjVj/H3yKyG9gGLDXGPLiG9FvgZRE5BNSxLbu01QcusPpAjCdcyZaHB3zftzyHPtpGseOjCLr9Dz5jfOk6bwih90IdHU8pu1sTdIE1Qbb/uy40IRUkUZ+IMaZZLOvPAXmiWXUDKB3DPpeB2tGtU8lToYLuBE7tzsRZb9JzRQ/Gu3/B3K9+ZlazidTzqe7oeErZze7+dR8tPKhEtE9EqfgTgc4tcnN29K+8cmkpV66H8uq8Grw8rjUXb190dDyl7C842PrpImciWokoh8icGZaPeY0VDfeRJfBT/jo/mzxDi/LduincN/cdHU+pRPXB7J18MHunteBizVlaiTiR1J7upPZM+lF87aluzbScmfk1XdwDCDvtS+/17Sn6bTX2no92hBylnNI/hy/xz+FL1oKLNWfpfCIq2ThwwNBwwHQO5f8ISX2ddiV68X2jfqTzSufoaEolnm3boGJF+OMPeP11R6eJM51PRCV7xYoJQXPaMLLgQTz3t2LKgWHk/qY4vwUucXQ0pRKPNmcpRxm9+hCjVx9ydAy7cnODHh2zcmLMFGoe28j1Cxl4c/4bVJ/QhJPXTzo6nlIJ0vEnfzr+ZGsB0UpEOUqkdlUXlyMHrJn+Egvr7STT9m/ZcHo5Bb8rxpC133Hv/j1Hx1MqXnb8d5Ud/121FlysT0QrEZWsNWrgxcm5H9P6TiD3Dtfgiw29Kfy/8mw5tcXR0ZSKM/8vXsb/i5etBb3EV6mklT49TPuuAFs/XEK+rQs4fuESL06pTMt5nbkafNXR8ZSKH23OUsoxKlQQDv3ehIHeB3Df3oOZ+yeT738+zNw9Swd1VMla2+nbaDt9m7WgzVnKUTKn9SJzWi9Hx3AoT0/o90kGDo7+jgq7/bl5Kh8tF7Wg8sSX+ffyv46Op1S0Ak/fIPD0DWvBxZqz9D4R5bSMgRk/hdN12mTuvPgp7qmC+eSlT/mixiek9nCN/6DKBX37LXz6Kdy541RnI3qfiHI5ItC6lTvHf+1MkzNBhO97kyH/DKTwd6X46+hfjo6nVPS0T0Q5ytDlQQxdHuToGMmOtzcsmJGD5e1n8+yKlZw6ZXh55su8Pfc9zt867+h4StFy6lZaTt1qLQQHQ6pU1rcgF6CViBPZeeIqO0/o1UgxeeUVOLLqZT702gvr+zF//3wKjizKRP+JOqijcqjD529x+Pwta8GF5lcHrUSUi0mXDr4fnhr//w3EZ+0e7hwuR+elnSk/oTK7z+12dDyVQm36tDabPrVNf5TSKxHblLWuNZSscjnlysGeNUX5tvhfeC7+mYDjxyg7qRy9VvTm1t1bjo6nUrKQEKfqUI9NrJWIiLiJSHMRWSoiF4Ag4KyI7BeR/4nI8/aPqVT8eXrCxx8LB355j6p7grjv/3+M3PIdhb8vxsIDC/XeEpVkmv+wheY/2EZZCA5OcWcia4FCwKdADmPMc8aY7MBLwBZgqIi0sGNGZZPzmdTkfMZ1/viSSqFCsG5ZZqa/NZEMv2zi/PEsNPmlCfVnv8Hxa8cdHU+lAP9ducN/V+5YCy7WnBXrfSIiUtIYszeWbTyNMWGJmsxO9D6RlO3CBeje4x7zjo9CavUnVSrDwFr96VmpJ57uno6Op1KCunXh1i3YtMnRSeLlae4T+V1EponIczFt4CwViFLZs8Pc2R4s/aI3ORbsJySwLh//9TGlJ5Thn//+cXQ8lRKkwOYsH2AnsEFERomIt50zqRgMXBLIwCWBjo7hEl57DQ5uy0u37Athzu/8e/wmL017ifaL23P5zmVHx1Mu5p2Jm3hnou3Mw8Was2KtRIwxd40xY4BiwElgm4h8JSIZ7Z5ORbL/zA32n7nh6BguI0MGGD0aNk9/g8KrAuGfj/hx5zSKjPFhRsAM7XhXieb8zVDO3wy1FlJaJfKAMSbEGDMc8AWCgR0i0sduyZRKIpUqwe7t6fmq2jDcp+zk+tHCtP69NTVn1OTAxQOOjqdcwPqParL+o5rWQkq7xPcBEckvIvWA/wPyAjeBr+O471QRuSAi+2JY7yMim0UkNLqKSUTcRWSXiPwRoWy6iBwTkQDbwy+u70WpqLy84IsvYO9fpai0/29YPJl/Du+h9MTSfLHmC4LDgh0dUbmKlNYnIiJ7ROQysBBoBTwDrAZaAunj+DrTgXpPWH8F6A4Mj2H9h0B0Xwk/Msb42R4BccyiVIx8fGDDejcmtG9Pqh+CuL/nXYZsHILveF9WHF7h6HjKSTUe/w+Nx9su3EiBzVmNgGzGmDLGmGbGmEHGmF+NMfuMMXfj8iLGmA1YFUVM6y8YY7YDj13lJSJ5gNeBKXF5LVdW0DsdBb3TOTqGy3Nzg06dIGhHdl4L/Qmmr+HMKU/qzapH0/lNOXPzjKMjKidzIziMG8G2j7eU1pxljDlqjDEiUlhEfhSRsUkRLILvgb5AdCPoDbGdKY0UkVQxHUBEOoiIv4j4X7x40V457e6bJqX4pkkpR8dIMfLkgd9/h3lDa5Lh593IukEs2Pc7PmN9GLN1DOH3wx0dUTmJ1b1rsLp3DWsSnJTWnBXBTGA+UA1ARHxF5Ce7pLIRkfrABWPMjmhWf4p1+fELQBbg45iOY4yZbIwpb4wp7+2tVyiruBOBd96BoMBUtMr/JeFj9hF+4kW6L+9OxSkV2XEmuj9NpWJw7x7cv59iKxE3Y8wyIBzAGLMP60ote6oCvCEix4G5QC0R+dn2+meNJRSYBlSwcxaH+/S3PXz62x5Hx0iRsmSBadNg1bznyb5yOfw6lwOnTlNhSgU+XPYhN0L10msVswZjNtJgzEaXm5AK4leJnBGRAoABEBEB7NqwZ4z51BiTxxiTH3gXWGOMaWF7/ZwRcjQCor3yy5UcvXiboxdvOzpGilanDuzbK/R6pSnBw4NIs7cLY7aNwWesD78G/qr3lqhohd67T+i9+48qkZTUJxJBD+AHIIeItME6M4jTB7eIzAE2A0VF5JSItBORTiLSybY+h4icAnoBX9i2ie1mxlkishfYC2QDBsfjvSiVYOnSwYgRsHX9MxT6dwxm8laCL+bgnfnv8Nrs1zh69aijI6pkZmXP6qzsWd3qDwGXOhPxiOuGxpjjtvtEGgGlgfXA1Dju2yyW9eeAPLFssw5YF2G5VlxeWyl7eeEF8PeHYcNeYOBX20jz4njWyReUGF+CL6t9SZ/KffBy93J0TJWcpMTmLFtzEQDGmHvGmPnGmC+NMeONMSFRt1EqJfH0hM8/hz0BHpS7152Q4QfIeP51Pl/zOX4T/Vh/fL2jI6pk4LVRG3ht1IYU25y1VkS6iUjeiIUi4iUitURkBtZNiMrOiufKSPFcOmRZcuTjA+vXw7hvc3Nn2nxSzV/K+cvB1JhRg9aLWnPxtvNeWq4SkQs2Z8VlPpHUQFvgPaAgcBWrQ90NWAmMN8bssnPORKPziSh7O3kSOneGpSvvkPPdr7j4/HAyps7IsDrDaFOmDW4S71mplatYvx5q1IDVq6GWc7XIJ3g+EdvAi+ONMVWwxsyqDZQxxuQzxrR3pgpEqaTw3HOwZAnMmp6WsGXfwMQA0t8pwf8t+T+qT69O4AUdzj/FSol9Ig+ISC1gItAHaCwi5Z50l7hKfD3m7qLHXK2znYEING8OBw5A01ol+G/genJtm8reswfwm+THJ399wp2wO46OqZJI3ZHrqTtyfYrtE3lgKrAEa171gkA/QL9SJaGz10M4ez3E0TFUPGTLBj//DEuXCm572nB9SBBFglsy9J+hlBhfgqX/LnV0RJUEUnm4kcrDzSX7ROJTiZwwxiyyDb74pTGmoTHmebslU8qFvPYaBAZCl9bZ2P/Nj+T4cwMmNC3159TnzV/e5NSNU46OqOxoSbeqLOlWNWU3Z2FNj9tTL+dVKmEyZoRx42DDBsh4rSonPttFuavf8Oe/f1JsXDFGbh7Jvfv3HB1T2VMKb84qDnQGzorIUhEZIiJv2ymXUi6ralXYvRs+7etFwNhPSDcjkMJeVem1shcv/PACW09tdXRElchqj1hH7RHrHjVnpXKd7uT4TI/7pjGmCFAAqz/kEFDJXsHU48rmy0zZfJkdHUMlgtSp4euvYft2eC59QXb1WUqlE/M5d+MCL/74Il2WduFayDVHx1SJJGMaTzKm8XTJM5FY7xN5uKFIFqAnkB3YD/xkjLlqx2x2ofeJqOQmLMwai2vAAEiT6QYv9O3H6ltj8E7rzXevfEcz32ZoK7KL6N8fBg2yhoN3sn/TBN8nEsFcrHnVlwBpgb9FxOWHX1fK3jw94ZNPrCYu38IZWdX7eyoEbOPZ1M/x3m/vUffnuhy6fMjRMVVieDAhlZNVIE8Sn0rE2xgzzBjzhzHmG6ABMNpOuVQ0Os3cQaeZOgmSqypa1DZ0yjjYt6ochz/dwltpx7Lt9DZKTijJoPWDCL0X6uiYKgGq/28t1f+31uXmV4f4VSJXRKTkgwVjzFGsMxKVRK7eucvVO3Ga1l45KTc36NIF9u2DqlXcmd+3K0VXBlE7d2P6r+tPqYmlWHNsjaNjqnh6NkMqns2QKsVXIl2BOSIyQUS62OZaP2KnXEqlaPnywbJlMGMGHA7IyerOc2jtsZzw++HU/qk27y98n/O3zjs6poqjXzpV5pdOla1KxIU61SF+V2cFAWWBtVid67uBJ84TopRKOBFo2RL274cGDWD6F6+Q/qe9tC/yJfP2zcNnnA+T/Cdx39x3dFQVVw/6RFxIXOYTGWD7WQVIZYz5xRgzwBjzw4P5RJRS9pMjB/z6KyxYAOdOpWFay0G0Dd1D6exl6LS0E1WmVmH3ud2Ojqme4KWha3hp6JoU25y1wvbzQ2CriBwSkd9F5Cu92TBpVXk+G1Wez+boGMpBmjSxBnRs0QImDfHh3LDV9POdyZErRyg3uRx9Vvbh1t1bjo6popE3S1ryZknrks1Z8blPJL9titxUQAmgJOBrjPnIngETm94nolzBihXQoYM1d0nHHlcIqfoJ0/f8wHMZn2PMq2No6NPQ0RFVdF56ybpbffVqRyeJt8S4T+Q3AGNMqDFmpzFmBrAgsQIqpeLulVesK7i6dIGJI7Owoc9kRpf+h0ypM9FoXiPemPMGJ66dcHRMFVVKbM4SkXdE5Fsgg4gUE4k0Ldtk+0VTUbWauo1WU7c5OoZKJjJkgLFjrQEd3d2he+PKVNi1g6+q/o/Vx1ZTfHxxhm8aTlh4mKOjpniVv1lN5W9Wu2RzVlzORP7BGuYkMzACOCwiO0XkDyDYnuFUZCFh4YSEhTs6hkpmHgzo+NFHMO1HTya27MMYnwPUKViHj1Z9RLnJ5dh0cpOjY6Zozz+bnuefTZ8yr84yxpw2xvwENAbeMcYUBF4G+gPONUmwUi4qTRoYNgy2bIHMmaHdm3nJ8Mfv/PTqIq6FXKPK1Cp0WNKBK8FXHB01RfqpbUV+alsxZTZnRTAJuAdgjLkMnABq2iOUUiphXngBduywBnOcNw9612/IkFz76fNiH6bumorPWB9m7p5JXC+oUYkshVciIRHvCzHGXAIGxXVnEZkqIhdEZF8M631EZLOIhIpIn2jWu4vILlsz2oOyAiKyVUQOi8g8EfGKx/tRyiV5eVmDxe7cCXnzQst303N00v9Y8eYOCmUpRMtFLan9U22CLgU5OmqKUXHIX1Qc8leK7RN54KiIvBqlLD4f2tOBek9YfwXoDgyPYf2HwIEoZUOBkbZpeq8C7eKRx+nULpad2sWyOzqGchIlS1rNW0OHwtKl8HbV0nTy/IcJr09k17ldlJpQii/XfElwmHZt2luJ3BkpkStDyuwTiaAb8I2IzBaR7iIyiXiMnWWM2YBVUcS0/oIxZjvw2KUkIpIHeB2YEqFMsPpk5tuKZgCN4prHGXWoVogO1Qo5OoZyIh4e0Lev1fFerBi0buXGkv4dWd04iKa+TRm8cTAlJ5Rk5ZGVjo7q0qa2rsDU98qAMSm3EjHGnAXKYd0b4g3sAZrbKVdU3wN9gYiDBGUFrhljHkxKfQrIHd3OItJBRPxFxP/ixYt2DapUclS0qHUp8KhRsG4d1Cj/LFUvzmRVi79wd3PnlZ9f4d3573L25llHR3VdLjirIcTvTARjTLgxZoEx5ktjzDhjzG17BXtAROoDF4wxCZ5Iwxgz2RhT3hhT3tvbOxHTJa2mkzbTdNJmR8dQTsrdHbp3h717oXx56NgRvu1Ym9/r7WFgjYEsClqEzzgfxm0bR/h9vZQ8MZUfvIq6Xy+3FlLqmYiI1BKRH0VkhIi0EZFytiFQ7K0K8IaIHMeaXbGWiPwMXAYyiYiHbbs8wOkkyKOUUytY0Bp1Y9Ik2LYNypVOReY9/djdcS8Vclfgg2Uf8OKPL7Lz7E5HR3UZ5fJmplx22xlISq1EgKlYU+NuAQoC/YBAe4SKyBjzqTEmjzEmP/AusMYY08JY1yiuBd6ybdoK+N3eeZRyBSLW2FuBgVC9unWG8n9NCjO2wkpmN5nNf9f/44UfXqDH8h7cCL3h6LhOb1LL8nzzehFrIQU3Z50wxiwyxvxqa85qaLsqKk5EZA6wGSgqIqdEpJ2IdBKRTrb1OUTkFNAL+MK2TcZYDvsx0EtEDmP1kfwYj/ejVIr33HPWlVvTp1tjcfn5CadXNCOwcxCdynVi9NbRFBtXjPn75+u9JU8r2HYVXAo+E9kgIj1tV0XFmzGmmTEmpzHG03Zm8aMxZqIxZqJt/TlbeUZjTCbb8xtRjrHOGFM/wvJRY0wFY8zzxpi3jTE6AbVS8SQCrVpZk1+98oo1fMprtTLRJf84tvzfFrKny87bv75N/Tn1OXb1mKPjOqUyg1bSctx6ayEFVyLFgc7AWRFZKiJDdD6RpFW/VE7ql8rp6BjKReXMCQsXwpw5cOQIlC0Lq6ZXYFPr7Yx8ZSQbTmyg+PjifLPxG+6G33V0XKdS5flslHvWNftE4jyfyMMdRNJgVSglgZLGmN72CGYvOp+IUrG7cAG6dYNffoEyZWDaNMha4BQfLv+Q3w78RrFsxZhYfyLV8lVzdFTnsXw5vPoqbN4MlSo5Ok28PfV8Irb7LH4EOgKZgKXOVoE4u+C74QTf1Usvlf1lz26NvbVgAZw5Y10SPOW7PMxptIA/mv1B8L1gqk+vTtvf23LpziVHx3UO2ifCG8CvWEOddASO2y67VUmk9bRttJ6m84mopNOkiXUF17vvwsCBVmWS4+brBHYJ5OMqHzNzz0x8xvowbdc07Xh/gtIDV/LxLNv/3ZRaiRhjzhhjlhtjhhlj3gHKo1dDKeXysmaFmTNh8WK4fBkqVoQhA9IysOq37Oq4C59sPrRd3Jbq06sTeMHuV/07pVo+2fHztlUeKfUSXxHJF3HZGHMAKJLoiZRSyVKDBtZlwO+/D19/bXW83znuy4Y2G5jSYAqBFwPxm+THp399yp2wO46Om6yMbOpHs5K20TJS2pmIiIyyPV1ou3djo4iMF5GvAV/7xlNKJSeZM1ud7H/+CTduwIsvwqefuPFe8XYEdQ2iRakWfPvPt5QYX4I/D/3p6LjJy4Oxs1JaJQI8vPzCGJMHaAOswhrBt4GdcimlkrFXX7XOStq1s2ZU9PODw3u8mdZwGutarSONRxpen/06b/3yFqdv6GhEJQesYNTSvdZCCmzOWi0im4FnRaQt8Azwp+1mwVP2jacieqtcHt4ql8fRMZQC4JlnYPJkWLnSuvCoShXo3RteyF6dgE4BDKk1hKWHluIzzodRW0Zx7/692A/qol71zUGJzJ7WnZ2eno6Ok6jidJ+IiBTCGqdqBtb9ISWAu8A+Y0xTuyZMZHqfiFKJ7+ZNa96SiROhcGGYOhVeegmOXj1K1z+7svzwcsrmLMvE1yfyQu4XHB3XMT76CMaNgzvO2V/0VPeJGGOOAHVsY2Y1MsYUBioCIxM5p3qCK7fvcuW23imskp8MGWDCBGt04LAwqFYNevSAZ70K8mfzP/nlrV84d+scFadU5IM/P+B6yHVHR05Sd+7eI+z2HZfrD4H4XeL7b5TlW8aYLYkfScWk88876PxzgqdVUcruatWy5ivp2tWaAKt0adi4UXi7xNsc6HqAbhW6McF/Aj7jfJi7b26Kubek4ter+X3LUZfrD4F4TkqllFKxSZ8exoyxZlA05tFQ8+73MjLq1VFs+79t5MmYh2YLmvHKz69w+MphR0e2uzdK5+L5DO4p+0xEKaXio3p12LPHqkDGjIFSpayKpVyucmxpt4Uxr45hy6kt+I735av1XxF6z3UH4R7SuKR1s6FWIkopFXfp0lnNWuvXWxcm1awJH3wAwXfc+aDCBwR9EERDn4b0W9eP0hNLs/bYWkdHtosrt+5y99Ydbc5SSqmEqFbNOivp0QPGj4eSJWHtWsiVIRfz3prHsveWEXY/jFo/1aLlwpZcuH3B0ZETVfXha9n571k9E1GO1aJSPlpUyhf7hkolQ2nTwsiRsGEDeHhYnfBdu8KtW1Dv+Xrs67yPz6t+ztx9c/EZ68PkHZO5b+47OnaieLNsHnKlwiUrkXjPJ+Ls9D4RpRzvzh34/HOrqSt/fvjxR6upC+DAxQN0XtqZ9SfW82KeF5lYfyKlni3l0LyJolw5yJULlixxdJIEeer5RJTjnbkWzJlrwY6OodRTi3hW4u4e+aykmHcx1rZay4xGMzh05RBlJ5Xlo5UfcfvubUfHTrDT1+7ofSLK8XrOC6DnvABHx1Aq0bz0EuzebfWVTJhgXcG1di2ICC1LtySoaxBt/NowfPNwio8vzuKDix0dOUHqfb+R8+evaSWilFKJ7UlnJVnTZuWHN37g7zZ/kzFVRhrObUijuY347/p/jo4dL80rPMczbuEuWYlonwgQFhbGqVOnCHkwVHMydfGmdR29d4ZUDk7y9FKnTk2ePHnwdLHB6NTTidpXMnUq1KhhrQsLD+P7Ld8zYP0AAAbWGMiHFT/E091J/oayZIEWLWD0aEcnSZCY+kS0EgGOHTtGhgwZyJo1KyLioGSxO3LxFgCFvNM7OMnTMcZw+fJlbt68SYECBRwdRyVDGzdC27Zw+LB1X8m331r3nACcuHaCbsu6seTfJZR6thQTX5/Ii8+96NjAsThy8RYF8nrj1q2bNXa+E3Jox7qITBWRCyKyL4b1PiKyWURCRaRPhPLUIrJNRHaLSKCIDIywbrqIHBORANvDL6H5QkJCkn0F4kpEhKxZsyb7Mz/lOFWrWn0lH35oDXxbqpTV3AWQL1M+fn/3dxY2XciV4CtUnlqZjks6cjX4qmNDP0GjsX/jFhLiks1ZSdUnMh2o94T1V4DuwPAo5aFALWNMacAPqCcilSKs/8gY42d7BDxNQGeoQLzTp8I7vfM3ZYFz/L6VY6VNC99/bw2VAtYwKh9+CLdvW38/jXwasb/LfnpW6smPu36k6Nii/Lzn52Q5qGPb8rmtJ3rHesIYYzZgVRQxrb9gjNkOhEUpN8aYW7ZFT9sj+f2FJJGMaTzJmMZJ2n+VSiQP7nbv1s3qTrBGBrbWZUiVge9e+Q7/Dv4UyFyA9xe+T+2fanPw0kHHho6iZ9XnrCd6JpL0RMRdRAKAC8AqY8zWCKuHiMgeERkpIjF+RReRDiLiLyL+Fy9etHdkuwkJCyckLNzRMZRKcunSWRXI2rVw/751VtKz56P5nfxy+LGp7SbGvzaenWd3UmpiKfqv7U/IveTRZHrwhG0YF61Ekp4xJtwY4wfkASqIiK9t1aeAD/ACkAX4+AnHmGyMKW+MKe/t7W3vyAnm4eFBkSJF8PPzw8/Pj0GDBkVaf/paMKejudkwODiY6tWrEx5uVTCTJk0iZ86cD4/j5+fH7t27qVatGvfuJXyK0iNHjlCyZMlIZaGhoRQoUIDAwMAEH1epuKpRwzor6dzZaury84NNm6x17m7udH6hM0EfBPFW8bcYtGEQJSeUZNWRVQ5MbOk+bbP1RJuzHMcYcw1rit56tuWztuauUGAaUMGB8RLF0KFD8fX1ZdeuXQQEBNCvX7847Td16lSaNGmCu7s7AHv37mXw4MEEBAQ8fJQuXZratWszb968BOcrUKAAp06d4v79R+MZTZ48mWrVqlGiRIkEH1ep+Eif3upsX70a7t61bljs08ea5x0gR/oczGoyi1Xvr0IQ6v5cl+YLmnPu1jmHZW5RKrv1RM9EkpaIeItIJtvzNMDLQJBtOaftpwCNgGiv/HImvXv3Jk2aNAwePDhe+82aNYuGDRs+XN6zZw9+fn6PbdeoUSNmzZqV4Hxubm7kzZuX48ePA9YZ0IgRIxg4cOCTd1TKDh7MotihA4wYAWXKwNYIjd11CtZhT+c99K/enwUHFuAz1ofx28cTfj/pm4TfL5PDeuKClYhHUryIiMwBagDZROQU0B+rkxxjzEQRyQH4AxmB+yLSAygO5ARmiIg7VoX3izHmD9thZ4mINyBAANApMbL26AEBAYlxpEf8/KxT77iYMmUKpUuXpkaNGlStWjXW7e/evcvRo0fJnz//w7LAwEDatGmDm5v1HaFLly506NABX19ftm/fHv83EEGxYsUICgqiYMGCjBs3jgYNGkR6baWSUoYMMHEivPkmtGsHlSvDRx/BgAHW53Vqj9QMqDGA5iWb02VpF7r+2ZXpAdOZVH8SZXKWSbKcB49foChoJZJQxphmsaw/h9XnEdUeINp/aWNMrUSIluwcOHCA8PBwSpcuHaftL126RKZMmR4unzx5Em9vb/bs2fPYtu7u7nh5eXHz5k0yZMjwsLxOnTqcO/f4qf6QIUMineGAVYkcPHiQatWqMXbsWLZG/OqnlIO8/LJ1VtK7Nwwdag2UO2MGlLfdGlckaxFWvb+KOfvm0HNFT8r/UJ7uFbozqOYgMqTK8OSDJ4Khi3YzFVyyTwRjTIp6lCtXzkS1f//+x8ocISQkxJQsWdKsWLHiYdmFCxdM69atzcmTJ02Llq3MxWs3Tc+ePU2vXr3MqFGjzJUrV0y+fPkebr906VLTqFGjGF8ja9as5u7duwnOOG/ePNOhQwczePBg89lnnz2WsU2bNubu3bvm3r17kXJGJ7n83pVr+fNPY3LnNsbd3ZjPPzcmJCTy+it3rphOSzoZGSAm94jcZsH+Beb+/ft2zbT0mynGgDFbttj1dewJ8DfRfKY6/EM9qR/JuRLp27evad++/WPl/fr1M++88465efOmGTNmjFm3bl2k9Xny5DHBwcHGGGO++eYb88knn0R7/EuXLpmiRYs+VcY9e/YYPz8/U7BgQXP16tVoMxpjos0ZVXL5vSvXc/WqMa1aWZ9wJUsas3Pn49tsPrnZlJpQyjAA8/qs182xq8fsF2jBAitMQID9XsPOYqpEknXHekozfPhw1q9f//Cy3O7du3Pr1i2OHj2Kh4cH7l6p2bp9O1WqVIm0X926dfn7778B68qsmTNnPjxGmTJluHXLul9z7dq1vP7660+VsUiRIuzdu5cOHTo8bEaLmDF9emtcrx07djyWU6mkkikTTJ8OixfDxYtQoQIMHAhhEW5nrpSnEjs67GBE3RGsO76O4uOKM/TvoYSFh8V02AQ7dMJ2f5o2Zzn/IzmfiUQVFhZm2rRpY44fP26GDRtmfl74p5kwY45p27at6d27t7l8+bIxxpgdO3aYFi1axHq8xo0bm4MHD9o149q1a40xxixatOixnFEl19+7ci2XLxvTvLl1IlC2rDF79z6+zX/X/jON5zY2DMCUGFfCbDyxMVEzfNmghxXgxIlEPW5SQpuznK8SierwhZvm8IWb0a778ccfzb1792LcNzQ01MyYMcNe0RLEWX7vyjUsWGCMt7cxXl7GfP21MWFhj2+zOGixyTsyr2EApu2itubi7YuJ8tqbeg6wPm7Pn0+U4zlCTJWINme5iLZt2z682TA6Xl5etGzZMgkTKZW8NGkCgYHQsCF89pl1k2JQUORtGhRtwP4u++lbuS8/7fkJn7E+TNs1zfrG/RRezG2bvsEFm7O0ElFKpRje3vDLLzB3Lhw6ZN2gOGIEhEe4/zCdVzqGvjyUnR12UjRbUdoubkuNGTXYf3F/gl/32Elbn4gL3ieilYhSKsVp2tQ6K6lb1xoypUYNawKsiEo+W5KNbTbyQ4Mf2Ht+L34T/fh89efcCbsT79dbtfME4eIGHklya16S0krEieTImJocGV3vm4xSjpAjByxaZN2UuHevNcT8uHHWKMEPuIkb/1f2/zj4wUGalWzG139/je94X5YdWhav13rpuQyYVKnABefR0UrEiaRL5UG6VK73TUYpRxGBli1h3z5rNsUPPrDOTv77L/J23um8mdFoBmtbrSWVRypem/0ab//6NqdvnI7T6xTP7IVHurR2eAeOp5WIE7kdeo/boQkfyl0pFb08eWDZMpg82RrE0dcXpk6FqP3pNfLXYHen3QyuOZg//v2DYuOKMXrr6FgHdTx15goh7l52fAeOo5WIEzl3I4RzN5LHJDtKuRoRaN/eatoqV84a0LF+fThzJvJ2Xu5efF7tc/Z13kfl5yrz4fIPqTilIv5n/GM89u7D57iQ+PcwJgtaiSilVAT581tzlYwaZc2k6OsLs2c/flZSKEshlr23jLlvzuXMzTNU+KEC3f7sxvWQ648ds3S21DyTJWPSvIEkppWIUkpF4eYG3btb00IULQrvvQdvv20NoRKRiNDUtykHuh7ggwofMG77OIqNK8Yvgb9EurckT2p4JrP9Rwt2BK1EkpHoprXdu3evo2MplWIVKQJ///1oePkSJawruqJ6JvUzjH51NNvabyNXhlw0nd+UerPqceTKEQAuXrrO5XDX/Lh1zXflpKKb1jbqnOZKqaTl7g59+8KOHfDcc9C4Mbz/Ply9+vi25XOVZ+v/bWV0vdFsPrkZ3wm+DN4wmGPnL3PwumteFKOVSDIS07S2D+R6JjW5non9PpHly5c/PJOpWLFipDnRlVIJ4+sLW7ZA//4wZ461vHz549u5u7nTrWI3DnQ9QIMiDfhy7Ze8/+YxzpV0vXtEAORpx4RxNuXLlzf+/pGvojhw4ADFihWzFhw4P27WrFnJnTv3Y9PaxlfhwoXZsGEDOXPmjPe+SSnS710pJ+LvD61awf791hzvw4dbU/VGZ9mhZXSd3JBj6cNoVboV/3v5f3in807awIlARHYYY8pHLdczkWQi4rS2D5qyolYgNWvVpniJEvj6+kZ6/P7775G2e+211yhVqhQ9evR4WNa4cWOMMXz99deMGjWKu3fv0q1bN3r06EHr1q05fTpuN00ppaxpd3fssIZM+eEH6273DRui3/bVwq+ydV52up4pyuy9syk6tihTdk7hvnGNFgK9/TmqOJwx2MPevXtj/VY+ZZ5VWRTyTh/jNps2bcIYw9mzZ/GwjdNz8uRJ0qdPT7t27WjVqhXVq1dn5cqVeHt7069fv8R7E0qlIKlTw//+Z40K3Lq1Nf5Wjx4wZEg0g/VevUN1/1J0GfwbnZd2pv2S9kwLmMbE1ydS8lnn7vfUM5FkYs+ePfj4+Dz1cX799VeKFCmCh4cHxhhu3LjBjh07WLt2Le+88w7Vq1cHoHLlyoSGhvL+++/z888/P/XrKpVSvfQS7N4NnTvDyJHWjYrbt0feJpOEU6FYLop7F2ddq3VMbzidfy//S9nJZfl41cfcvnvbMeETgVYiycSTprWNj2bNmjFp0iRKlSpFpUqVOHToEDt27GD+/PlMnDiR8+fPA5A+fXqGDBnCjBkzGDduXGK/HaVSlHTprMEbV66EGzfgxRehXz+4e9da73k3lHy5swLWvSWt/FoR1DWI1qVbM2zTMIqPL86Sg0sc+A6eQnQzVbnyw1VnNnySxo0bm7CwMHPy5EnTsGFDc+fOHdO5c2fTrVs306ZNG7N48WI7pI2ds/zelYqPq1eNef99ayLDMmWM2bvnvjFgdrf7MNrtN57YaEqMK2EYgGk8t7H579p/SZo3rohhZkO9OgvnuUroyEXrzORJfSLOxFl+70olxMKF0LEjhF4L5npYWkbWakPP1VOj3TYsPIyRW0YyYN0A3MSNgTUG8mGlD/FwSz7d1g6/OktEporIBRHZF8N6HxHZLCKhItInQnlqEdkmIrtFJFBEBkZYV0BEtorIYRGZJyKuOUymTe5MacidyfWm11TKFTVubA0x/0Zda9DUsH+Lc+RI9Nt6unvSt0pf9nfdT438Neizqg/lJ5dny6ktSZg4YZKyT2Q6UO8J668A3YHhUcpDgVrGmNKAH1BPRCrZ1g0FRhpjngeuAu0SM3Byk9rTndSeMc+jrpRKXrJnh58mW5XImUsZKFUKJkx4fDDHB/Jnys+SZkv47Z3fuHTnEpV/rEynPzpxNTia2+OTiSSrRIwxG7AqipjWXzDGbAfCopQbY8yDHmZP28OIiAC1gPm2dTOARomdOzm5ERzGjWAXHU9aKRclIcEAvPJ/t6lSBbp0gVdfhZhuzRIRGhdrzIGuB/iw4of8sPMHfMb5MHvvbJJj94NTXJ0lIu4iEgBcAFYZY7YCWYFrxpgHA9KcAnLHsH8HEfEXEf+LUYfhdCIXb4Vy8Vaoo2MopeIjxDoTWX32EitWWFdxbdwY8xDzD2RIlYGR9Ubi396ffM/k473f3uPlmS/z7+V/kzB87JyiEjHGhBtj/IA8QAUR8Y3n/pONMeWNMeW9vZ1vuAGllBOzVSLt6xRDxDoTCQiAYsWsIebfeQcuXYp59zI5y7C53WbGvTYO/zP+lJxQkgHrBhByL3lMUOcUlcgDxphrwFqsvpXLQCYReXD5Qh5Ax+5QSiUvtkqkaMFnHxYVLmydjXzzDfz+u3VW8scfMR/C3c2dLi90IeiDIN4q/hYD1w+k1IRS/HX0L3unj1Wyr0RExFtEMtmepwFeBoJs1y2vBd6ybdoK+D3agyillKMEW30ivx+4HKnY3R0++cS6u/3ZZ6FBA2t63ps3Yz5UjvQ5mNVkFitbrMRgeHnmy7z323ucu3XOnu/giZLyEt85wGagqIicEpF2ItJJRDrZ1ucQkVNAL+AL2zYZgZzAWhHZA2zH6hN5UGd/DPQSkcNYfSQ/JtX7UUqpOLGdiczdF31/bOnSsG2bVaFMnQqlSsU8mOMDLxd6mb2d9/JltS+Zv38+PmN9mLB9gmMGdYzuDkRXfjjzHeuhYeEmNCzc0TESjbP83pV6KvPnGwNm//KNsW76zz/GFCpkjIgxvXsbExwc++GDLgaZWjNqGQZgKv5Q0ew6u+vpM0eDGO5YT/bNWQpat24NwFcD++Plof9kSjkVW3NWsQh9IjGpXNnqdO/UCUaMsAZz3LnzyfsUzVaUv97/i5mNZ3Ls2jHKTS5HrxW9uBn6hHaxRKSfSMmIh4cHRYoUeTgI46BBg7hz5w7p0qXj3Llz3LwTwv7Dx6hZsyYjR46kadOmD/cNDg6mevXqhIeHxzhX+927d6lWrRr37iVsms4jR448Nl1vaGgoBQoUIDAw8Kneu1Iuy9acNdX/bJw2T58exo+HZcvg2jWoWBEGD4Yn/bcVEVqUakFQ1yDal23PyC0jKT6+OAsPLLT7vSVaiSQjQ4cOxdfXl127dhEQEEC/fv3YuXMnZcuWJSAggAJFS7B5+06aN29Oz549H84XAjB16lSaNGmCu7t7jHO1e3l5Ubt2bebNm5egfAUKFODUqVORptudPHky1apVo0SJEk/9/pVySbZKZNqu+HV+16sHe/fCW2/Bl19aQ84fOvTkfTKnyczE+hPZ1HYTWdJkockvTXhj7hscv3Y8geFjp5VIMtK7d2/SpEnD4MGDH5Zt376dF154gYCAAIr5luRA4F6qVq0KWN8+Hpg1axYNGzYEnjxXe6NGjZg1a1aC8rm5uZE3b16OHz8OWGc/I0aMYODAgU/eUamUzFaJTO5QJd67Zslizec+Zw78+6810/b48THfoPjAi8+9iH97f4a/PJw1x9ZQYnwJhv0zjLDwxB/xIvkMEZlM9Fjeg4BzAYl6TL8cfnxf7/s4bTtlyhRKly5NjRo1qFq1KoGBgXTv3p1Ro0bxZpsunDh2lCJFinDp0iVy5MgBwN27dzl69Cj58+cHIDAwkDZt2kQ7V7uvry/bo86YEw/FihUjKCiIggULMm7cOBo0aPDwdZVS0XjQJ5I/9j6RmLz7LlStCu3aQdeu1r0lU6dC7mjH6LB4unvSu3Jv3i7xNt2Xdefjvz5m5p6ZTHx9IlXyxr9Ci4meiSQzBw4cIDw8nNKlSwNWpeLu7s6PP/6Im5sb334/Hjc3N7Jly8bw4dZYlZcuXSJTpkxA7HO1u7u74+Xlxc0oF6PXqVPnsbnbo5u/vVixYhw8eJBbt24xduxYvvjiCzv+NpRyASEh3HdzZ/jqGIbwjaPcua1+kgkT4O+/oWRJmDs39v3yPpOXRe8uYlHTRdwIvcFL016i/eL2XL5zOfad40DPRKKI6xmDPYSGhtK6dWsmTJhAxowZ47xfmjRpCLGdMsdlrvbQ0FBSp04dqeyvv+J252uxYsVYvXo1o0aN4r333uPZZxP+7UqpFCEkhGAPL2ZsPk6fV4o+1aFErCu3ateGli2hWTNYtMhq4sqS5cn7NvRpSO2CtRm0fhDfbf6ORQcXMfzl4bQs3TJS03h86ZlIMtKvXz8qVapE3bp1o12fL0ta8mVJ+1h55syZCQ8PJyQkJNa52i9fvky2bNnw9PRMUMZixYqxbds2pk6dykcfffSw/OLFi7Rp04ZTp07Rtm1bQkND6dWrF71792b06NEJei2lXEJICKkzpGVR18RrQnowbMqQIbBggXVWsmJF7Pul90rPsJeHsbPjTgpnKUzr31tTc0ZNDlw8kOAsWokkI8OHD2f9+vUPL8vt3r17pA/nDu3/j/B7YdF+ONetW5e///471rna165dy+uvv57gjEWKFGHv3r106NDhYRMagLe3N3nz5n2Y64cffqBhw4aMGDGC7t27J/j1lHJ6wcG4p0mT6DOSenjAZ5/B1q2QKZN1NVfXrnDnTuz7lnq2FH+3/ZvJ9Sez5/weSk8szRdrviA4LDj+QaK7A9GVH854x3q/fv3MO++8Y06cu2y+HT7SrFu37rFtduzYYVq0aBHrsRo3bmwOHjyY6Blv3rxpWrRoYZo3b26MMaZ169YmLCzsifsk99+7UomiWTNzOVc+8/XSQLu9RHCwMT17GgPGFClizNatcd/3/K3z5v3f3jcMwBQcVdAsP7Q82u3QO9ad061btzh69CgeHh6EuXmx3d+fKlUePy0uW7YsNWvWJDw8PMZj3b17l0aNGlGkSJFEzXjv3j26d+/O4MGD8fPzY926dTRq1IiOHTvSp08frlyJcS4ypVxfSAiX7gmzt52020ukTg3ffQerV1sXg1WuDP37Q1gcrujNni47PzX+iTUt1+Dp5km9WfVoOr8pZ26eidNri7Hz3YzJTfny5Y2/v3+ksgMHDsTaGe0I9+7do0OHDvTv359ffvmFXIV9uXnjOtvXryJz5sx89tlnZImtNy0ZS66/d6US1auvcvfCRS6u3kDuTI/3aSa2a9ege3eYORPKl4eff4aicezPD70XyrB/hjFk4xBSeaRicM3BdHmhC+5u7ojIDmNM+aj7aCWC83yYHblo9W0kdtuqozjL712pp1KzJoSHxz40byKbPx86drTOTIYNs/pL4noR1uErh+mytAurjq6iXM5yTKo/ifK5y0dbiWhzllJK2VNICIdvhjNgcdKOL/fWW7BvH9SoAd26WR3vZ+LWQsXzWZ5nRYsVzH1zLqdvnqbClAoxbquViFJK2VNICCeD77Ng56kkf+mcOWHp0kc3KPr6wi+/xG1fEaGpb1MOdD1A5/KdY9xOKxEnUiBrOgpkTefoGEqp+AgOpnKJPKzvU9MhL//gBsWAAChSBJo2hRYtrL6TuMiUOhNjXxsb43qtRGycoW/IzU1wc0v4naXJiTP8vpVKFCEhpMqQjizpvRwao3Bh62xk0CBruJSSJWHNmqc/rlYiQOrUqbl8+XKy/2C7fCuUy7dCHR3jqRljuHz58mNDryjlkkJC2Hb2Dp8v3OvoJHh4WMPKb94M6dJZw6f06vVwoOGEHTPx4jmvPHnycOrUKS5ejH4O5OTi4k2rAvHOkMrBSZ5e6tSpyZMnj6NjKGV/ISEcvB7G4t1nGNK4ZOzbJ4EXXrBmTOzbF0aOhJUrrUuBY5hB4om0EgE8PT0pUKCAo2PEasCkzQDM6+jn2CBKqbgLDubdqoV587Pajk4SSdq0MHYs1K8PbdtChQrw1VfQpw+4u8f9ONqcpZRS9nL/Pty9i2f6dKT1Sp7f2R/MoPjGG/DJJ9ZtLbZ55+JEKxGllLKXUKsJ+s9DV+k7f7eDw8Qsa1b49VeYMcO6iqtUKet5XLqJtRJRSil7sfVY77kUyrJ98ZtjPamJWHOU7Nlj9Y20bg1vvw2XY5m7KsUNeyIiF4ETiXS4bMClRDqWvTlLVs2ZuJwlJzhP1pSaM58xxjtqYYqrRBKTiPhHN5ZMcuQsWTVn4nKWnOA8WTVnZNqcpZRSKsG0ElFKKZVgWok8ncmODhAPzpJVcyYuZ8kJzpNVc0agfSJKKaUSTM9ElFJKJZhWIkoppRJMK5F4EJEsIrJKRA7ZfmaOYbu8IrJSRA6IyH4RyZ/EUeOc1bZtRhE5JSIxTxpgJ3HJKSJ+IrJZRAJFZI+INE3CfPVE5KCIHBaRT6JZn0pE5tnWb3XEv7UtR2w5e9n+FveIyGoRyZccc0bY7k0RMSLisEtp45JVRN6x/V4DRWR2Ume0ZYjt3z6viKwVkV22f//XEjWAMUYfcXwAw4BPbM8/AYbGsN064GXb8/RA2uSa1bZ+FDAbGJsccwJFgMK257mAs0CmJMjmDhwBCgJewG6geJRtugATbc/fBeY54HcYl5w1H/wdAp2Ta07bdhmADcAWoHxS54zH77QwsAvIbFvOnkxzTgY6254XB44nZgY9E4mfhsAM2/MZQKOoG4hIccDDGLMKwBhzyxhzJ8kSPhJrVgARKQc8C6xMmliPiTWnMeZfY8wh2/MzwAXgsTtn7aACcNgYc9QYcxeYa8sbUcT884HaIpLUM4fFmtMYszbC3+EWwBHj8Mfl9wnwFTAUeIpZLp5aXLK2B8YZY64CGGMuJHFGiFtOA2S0PX8GiONM63GjlUj8PGuMOWt7fg7rwzeqIsA1EfnNdvr4PxGJx8DKiSbWrCLiBowA+iRlsCji8jt9SEQqYH3jOmLvYEBu4GSE5VO2smi3McbcA64DWZMgW7QZbKLLGVE7YJldE0Uv1pwiUhZ4zhizNCmDRSMuv9MiQBER+UdEtohIvSRL90hccg4AWojIKeBPoFtiBkieYxM7kIj8BeSIZtXnEReMMUZEors+2gOoCpQB/gPmAa2BHxM3aaJk7QL8aYw5Zc8vz4mQ88FxcgIzgVbGmPuJmzJlEJEWQHmguqOzRGX7UvMd1v8XZ+CB1aRVA+vMboOIlDTGXHNkqGg0A6YbY0aIyIvATBHxTaz/Q1qJRGGMqRPTOhE5LyI5jTFnbR9o0Z2+ngICjDFHbfssAiphh0okEbK+CFQVkS5YfTdeInLLGBNjh6eDciIiGYGlwOfGmC2Jme8JTgPPRVjOYyuLbptTIuKB1VwQy7iniS4uORGROlgVd3VjjCPmWY4tZwbAF1hn+1KTA1gsIm8YY/yTLKUlLr/TU8BWY0wYcExE/sWqVLYnTUQgbjnbAfUAjDGbRSQ11uCMidL8ps1Z8bMYaGV73gr4PZpttgOZRORBm30tYH8SZIsq1qzGmPeMMXmNMfmxmrR+SuwKJA5izSkiXsBCrHzzkzDbdqCwiBSwZXgXK29EEfO/Bawxth7MJBRrThEpA0wC3nBQ2z3EktMYc90Yk80Yk9/2N7kFK29SVyCxZrVZhHUWgohkw2reOpqEGSFuOf8DagOISDEgNZB4c4En9dUEzvzAauteDRwC/gKy2MrLA1MibPcysAfYC0wHvJJr1gjbt8YxV2fFmhNoAYQBAREefkmU7zXgX6w+mM9tZYOwPtzA+g/5K3AY2AYUdNDfZmw5/wLOR/j9LU6OOaNsuw4HXZ0Vx9+pYDW/7bf9X383meYsDvyDdeVWAFA3MV9fhz1RSimVYNqcpZRSKsG0ElFKKZVgWokopZRKMK1ElFJKJZhWIkoppRJMKxGVYolId7FGWp7l6CyJSUS+F5FqT1jfSkTmRCnLJiIXbaMSzxWRwvZPqlyBXuKrUiwRCQLqGGNORSn3MNY4WE5HRLICS40xlZ6wTUasm+LyGtugjCLSCahgjGkrItWBFsaY9kkSWjk1PRNRKZKITMQaPnuZiPQUkQEiMlNE/sEaW8hbRBaIyHbbo4ptv6xizRUTKCJTROSE7Vt8fhHZF+H4fURkgO15IRFZLiI7RGSjiPjYyqeLyGgR2SQiR0XkrQj7fywie0Vkt4h8azvGzgjrC0dcjuBNYHmE7cqJyHrba6+wDTFzA1gPNIiw37vAg7OTjUAd2zAuSj2RViIqRTLGdMIaErumMWakrbg41plJM6w5VkYaY17A+mCeYtumP/C3MaYE1lAseePwcpOBbsaYcljDy4yPsC4n8BJQH/gWQERexRrOu6IxpjQwzBhzBLguIn62/doA06J5rSrADttxPIExwFu2154KDLFtNwer4kBEcmEN2bHG9ru5j3UHfuk4vDeVwuk3DaUeWWyMCbY9rwMUjzC6cUYRSQ9UA5oAGGOWisjVJx3Qtk9l4NcIx0oVYZNFtg/t/SLyYBj8OsC0B01NxpgrtvIpQBsR6QU0xZpLIqqcPBoXqSjWgIarbK/tjjWhF1iDWY63NW29AywwxoRHOM4FrAnAdjzp/SmllYhSj9yO8NwNqGSMiTQxksQ8ZP49Ip/Zp45wnGvGGL8Y9os4mm5s4/EvwDoTWgPsMMZEN1pwcITXFiDQGPNi1I2MMcEishxojHVG0ivKJqltx1LqibQ5S6norSTC5D0RmpE2AM1tZa8CD+aEPw9kt/WZpMJqnsLW/3BMRN627SMiElsz0SqsM460tn2y2I4VAqwAJhB9UxbAAeB52/ODgLdYc0ggIp4iUiLCtnOwKo9ngc1RjlME2IdSsdBKRKnodQfKi8geEdkPdLKVDwSqiUggVrPWfwDGmlNiENZIvquAoAjHeg9oJyK7gUCinxL2IWPMcqzhvP1FJIDIM0/OAu4T83TGS7ENT26s6VLfAobaXjsAq2ntgVVYTVbzTITLNG3NasHGmHNPyqkU6CW+Sj0VETmONVz5pSR6vT7AM8aYL5+wzd9AfZPAGfZEpCdwwxiT6BOpKdejfSJKOQkRWQgUwpro7El6Y101di2BL3UNaxpipWKlZyJKKaUSTPtElFJKJZhWIkoppRJMKxGllFIJppWIUkqpBNNKRCmlVIL9P7jlWeWQZx0+AAAAAElFTkSuQmCC)
%%%% Output: stream
%% Cell type:code id: tags:
``` python
```
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment