{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Enhanced Sampling Simulations\n",
"## with Ensembler"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this notebook, we give examples for enhanced sampling simulations with Ensembler. We show how to execute, visualize, and analyze those simulations for both 1D and 2D systems.\n",
"\n",
"The enhanced sampling technologies are briefly explained in order to directly use this notebook for teaching purpose.\n",
"\n",
"Maintainers: [@SchroederB](https://https://github.com/SchroederB), [@linkerst](https://https://github.com/linkerst), [@dfhahn](https://https://github.com/dfhahn)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Loading Ensembler and necessary external packages "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import os, sys\n",
"my_path = os.getcwd()+\"/..\"\n",
"sys.path.append(my_path)\n",
"\n",
"import matplotlib\n",
"matplotlib.use('TkAgg')\n",
"\n",
"import numpy as np\n",
"%matplotlib inline\n",
"from matplotlib import pyplot as plt\n",
"\n",
"from ensembler.potentials.OneD import fourWellPotential, harmonicOscillatorPotential, addedPotentials, metadynamicsPotential\n",
"from ensembler.potentials.TwoD import harmonicOscillatorPotential as harmonicOscillator2D, wavePotential as wavePotential2D, gaussPotential\n",
"from ensembler.potentials.TwoD import addedPotentials as addedPotentials2D, metadynamicsPotential as metadynamicsPotential2D\n",
"from ensembler.samplers.stochastic import langevinIntegrator, langevinVelocityIntegrator\n",
"from ensembler.system import system\n",
"from ensembler.ensemble import replica_exchange\n",
"\n",
"##Visualisation\n",
"from ensembler.visualisation.plotSimulations import simulation_analysis_plot"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Enhanced Sampling Simulations in 1D"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Unbiased System"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We start our walkthrough with an unbiased reference simulation of a four-well-potential and give a small introduction in how to set up a simulation with Ensembler. The four-well-potential is defined by the x-positions ($a$-$d$) and the y-position $ah$ - $dh$ of the four wells. If wished, the potential can be scaled in the y-direction using $V_{max}$. Note, that the energy is given in units of $k_BT$. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Simulation Setup\n",
"pot = fourWellPotential(Vmax=4, a=1.5, b=4.0, c=7.0, d=9.0, ah=2., bh=0., ch=0.5, dh=1.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here, we choose the stochastic Langevin integrator for the numeric integration of our system. The user has to set the step size $dt$ and the friction coefficient $gamma$. Note, that this integrator already contains a thermostat. The temperature of the simulation will be set during the system setup (see below)."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Simple Langevin integration simulation\n",
"integrator = langevinIntegrator(dt=0.1, gamma=15)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The system class wraps the integrator and the potential. Additionally, the initial position of the particle $position$ as well as the temperature parameter $temperature$ are set."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# Put Potential and Integrator together to generate the simulation system \n",
"system1 = system(potential=pot, sampler=integrator, start_position=4.2, temperature=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To start the simulation we define the number of steps and run `sys.simulate()`. The progress of the simulation is displayed by a progress bar."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\r",
"Simulation: Simulation: 0%| | 0/2000 [00:00\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
position
\n",
"
temperature
\n",
"
total_system_energy
\n",
"
total_potential_energy
\n",
"
total_kinetic_energy
\n",
"
dhdpos
\n",
"
velocity
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
4.200000
\n",
"
1
\n",
"
0.158629
\n",
"
0.158622
\n",
"
0.000008
\n",
"
1.595959
\n",
"
-0.003895
\n",
"
\n",
"
\n",
"
1
\n",
"
4.165484
\n",
"
1
\n",
"
0.108275
\n",
"
0.108275
\n",
"
NaN
\n",
"
-1.595959
\n",
"
NaN
\n",
"
\n",
"
\n",
"
2
\n",
"
4.125017
\n",
"
1
\n",
"
0.061324
\n",
"
0.061324
\n",
"
NaN
\n",
"
-1.321306
\n",
"
NaN
\n",
"
\n",
"
\n",
"
3
\n",
"
4.230917
\n",
"
1
\n",
"
0.211764
\n",
"
0.211764
\n",
"
NaN
\n",
"
-0.999128
\n",
"
NaN
\n",
"
\n",
"
\n",
"
4
\n",
"
4.118459
\n",
"
1
\n",
"
0.054943
\n",
"
0.054943
\n",
"
NaN
\n",
"
-1.841802
\n",
"
NaN
\n",
"
\n",
" \n",
"
\n",
""
],
"text/plain": [
" position temperature total_system_energy total_potential_energy \\\n",
"0 4.200000 1 0.158629 0.158622 \n",
"1 4.165484 1 0.108275 0.108275 \n",
"2 4.125017 1 0.061324 0.061324 \n",
"3 4.230917 1 0.211764 0.211764 \n",
"4 4.118459 1 0.054943 0.054943 \n",
"\n",
" total_kinetic_energy dhdpos velocity \n",
"0 0.000008 1.595959 -0.003895 \n",
"1 NaN -1.595959 NaN \n",
"2 NaN -1.321306 NaN \n",
"3 NaN -0.999128 NaN \n",
"4 NaN -1.841802 NaN "
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"system1.trajectory.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The simulation results can be visualized using the build-in visualizing tool. The left panel displays the energy surface as well as all visited states. The start and end-state are colored in blue and red, respectively. The middle panel shows the probability density distribution of the simulation as well as a boxplot of the distribution. This plot can be used to check if the system was simulated successfully. The rightmost panel shows the development of the force over time. \n",
"\n",
"Note that without enhanced sampling only the minimum around x=4 is sampled. It would require a very long simulation time to overcome the energetic barrier and hence a lot of computing time. Enhanced sampling methods were developed to speed up those slow processes. We will explore some of the most commonly used enhanced sampling methods in the subsequent notebook. "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"#plot\n",
"simulation_analysis_plot(system1, title=\"Position Langevin\", limits_coordinate_space=list(range(0,10)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Enhanced sampling"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Enhanced sampling methods can be divided into time-independent and time-dependent methods. Time-independent biases stay the same throughout the whole simulation whereas for time-dependent enhanced sampling the bias is updated during the simulation time."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Time-independent bias"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Umbrella sampling"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Umbrella sampling is a time-independent enhanced sampling method. In this method energy barriers are overcome by adding a bias potential. Note that umbrella sampling requires the choice of a reaction coordinate along which the bias is added. The choice of a suitable reaction coordinate is non-trivial for high dimensional systems. For our low dimensional 1D we can simply chose the x-axis.\n",
"\n",
"One of the most frequent umbrella sampling method uses hormonic potentials to restrain the sampling to a certain region of the potential. This is especially useful for sampling transition regions."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We start as in the unbiased case by defining the potential. However, we have to define two potentials; the original potential and the bias potential. The original potential is the same potential as defined above, the bias potential is a hormonic oscillator centered at $x_{shift}$ and force constant $k$. In this case we want to sample the transition region around $x$ = 5.5. Therefore, we set the $x_{shift}$ parameter to 5.5. The force constant $k$ defines how much we constrain the system. The higher the energy barrier, the more constrain is needed. \n",
"\n",
"To sample the full potential energy landscape, we can set up multiple simulations with different $x_{shift}$ parameters. For subsequent analysis of umbrella sampling simulations (e.g. using the weighted histogram analysis method WHAM) it is important that the sampling space of the different simulations overlap. The higher the force constant $k$, the more simulations are needed achieve the overlap."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"#Simulation Setup\n",
"origpot = fourWellPotential(Vmax=4, a=1.5, b=4.0, c=7.0, d=9.0, ah=2., bh=0., ch=0.5, dh=1.)\n",
"biaspot = harmonicOscillatorPotential(k=10, x_shift=5.5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `addedPotentials` function of the biasOneD class wraps any two 1D potentials together. Therefore, it is straightforward to generate the enhanced sampling system from the original and biased potential."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"#Add the bias and the original system\n",
"totpot = addedPotentials(origpot, biaspot)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"All subsequent steps are identical to the unbiased system. Note, that the starting position of the simulation has to match to the $x_{shift}$ parameter (i.e. should be reasonable close to $x_{shift}$) in order to avoid starting the simulation at very high energy states."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\r",
"Simulation: Simulation: 0%| | 0/2000 [00:00\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
position
\n",
"
temperature
\n",
"
total_system_energy
\n",
"
total_potential_energy
\n",
"
total_kinetic_energy
\n",
"
dhdpos
\n",
"
velocity
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
4.200000
\n",
"
1
\n",
"
8.609497
\n",
"
8.608622
\n",
"
0.000875
\n",
"
-11.404041
\n",
"
-0.04184
\n",
"
\n",
"
\n",
"
1
\n",
"
4.285465
\n",
"
1
\n",
"
7.699532
\n",
"
7.699532
\n",
"
NaN
\n",
"
11.404041
\n",
"
NaN
\n",
"
\n",
"
\n",
"
2
\n",
"
4.495405
\n",
"
1
\n",
"
6.021825
\n",
"
6.021825
\n",
"
NaN
\n",
"
9.870341
\n",
"
NaN
\n",
"
\n",
"
\n",
"
3
\n",
"
4.505998
\n",
"
1
\n",
"
5.958023
\n",
"
5.958023
\n",
"
NaN
\n",
"
6.117316
\n",
"
NaN
\n",
"
\n",
"
\n",
"
4
\n",
"
4.614714
\n",
"
1
\n",
"
5.418150
\n",
"
5.418150
\n",
"
NaN
\n",
"
5.928966
\n",
"
NaN
\n",
"
\n",
" \n",
"
\n",
""
],
"text/plain": [
" position temperature total_system_energy total_potential_energy \\\n",
"0 4.200000 1 8.609497 8.608622 \n",
"1 4.285465 1 7.699532 7.699532 \n",
"2 4.495405 1 6.021825 6.021825 \n",
"3 4.505998 1 5.958023 5.958023 \n",
"4 4.614714 1 5.418150 5.418150 \n",
"\n",
" total_kinetic_energy dhdpos velocity \n",
"0 0.000875 -11.404041 -0.04184 \n",
"1 NaN 11.404041 NaN \n",
"2 NaN 9.870341 NaN \n",
"3 NaN 6.117316 NaN \n",
"4 NaN 5.928966 NaN "
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"integrator = langevinIntegrator(dt=0.1, gamma=15)\n",
"system2 = system(potential=totpot, sampler=integrator, start_position=4.2, temperature=1)\n",
"\n",
"#simulate\n",
"sim_steps = 2000\n",
"cur_state = system2.simulate(sim_steps, withdraw_traj=True, init_system=False)\n",
"\n",
"print(\"Trajectory length: \",len(system2.trajectory))\n",
"print()\n",
"print(\"last_state: \", cur_state)\n",
"print(len(system2.trajectory))\n",
"system2.trajectory.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In order to visualize the resulting potential we again use the plotting function `simulation_analysis_plot`. The visualization shows that umbrella sampling samples the high energy transition region around x=5.5 very well. In contrast, the unbiased simulation (see above) was stuck at the minimum around x=4. "
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"#plot\n",
"simulation_analysis_plot(system2, title=\"Position Langevin\", limits_coordinate_space=list(range(0,10)), oneD_limits_potential_system_energy=[0,30])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##### Scaled potential"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`addedPotentials` can add any potential classes with a given dimensionality. In the special case that the added potential is identical to the original potential but a prefactor, one can scale the potential. In the example below we chose a four well potential that we will scale down in order to cross the energy barriers more easily. The procedure to define the original and bias potential are the same as described above."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\r",
"Simulation: Simulation: 0%| | 0/2000 [00:00\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
position
\n",
"
temperature
\n",
"
total_system_energy
\n",
"
total_potential_energy
\n",
"
total_kinetic_energy
\n",
"
dhdpos
\n",
"
velocity
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
4.200000
\n",
"
1
\n",
"
0.021572
\n",
"
0.019828
\n",
"
0.001744
\n",
"
0.199495
\n",
"
-0.059064
\n",
"
\n",
"
\n",
"
1
\n",
"
3.993458
\n",
"
1
\n",
"
-0.000150
\n",
"
-0.000150
\n",
"
NaN
\n",
"
-0.199495
\n",
"
NaN
\n",
"
\n",
"
\n",
"
2
\n",
"
4.085934
\n",
"
1
\n",
"
0.003545
\n",
"
0.003545
\n",
"
NaN
\n",
"
0.006084
\n",
"
NaN
\n",
"
\n",
"
\n",
"
3
\n",
"
3.935761
\n",
"
1
\n",
"
0.001858
\n",
"
0.001858
\n",
"
NaN
\n",
"
-0.085983
\n",
"
NaN
\n",
"
\n",
"
\n",
"
4
\n",
"
3.871948
\n",
"
1
\n",
"
0.007934
\n",
"
0.007934
\n",
"
NaN
\n",
"
0.063492
\n",
"
NaN
\n",
"
\n",
" \n",
"
\n",
""
],
"text/plain": [
" position temperature total_system_energy total_potential_energy \\\n",
"0 4.200000 1 0.021572 0.019828 \n",
"1 3.993458 1 -0.000150 -0.000150 \n",
"2 4.085934 1 0.003545 0.003545 \n",
"3 3.935761 1 0.001858 0.001858 \n",
"4 3.871948 1 0.007934 0.007934 \n",
"\n",
" total_kinetic_energy dhdpos velocity \n",
"0 0.001744 0.199495 -0.059064 \n",
"1 NaN -0.199495 NaN \n",
"2 NaN 0.006084 NaN \n",
"3 NaN -0.085983 NaN \n",
"4 NaN 0.063492 NaN "
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sim_steps = 2000\n",
"\n",
"#Simulation Setup\n",
"origpot=fourWellPotential(Vmax=4, a=1.5, b=4.0, c=7.0, d=9.0, ah=2., bh=0., ch=0.5, dh=1.)\n",
"# same potential but Vmax is different \n",
"biaspot = fourWellPotential(Vmax=-3.5, a=1.5, b=4.0, c=7.0, d=9.0, ah=2., bh=0., ch=0.5, dh=1.)\n",
"#Add the bias and the original system\n",
"totpot = addedPotentials(origpot, biaspot)\n",
"\n",
"integrator = langevinIntegrator(dt=0.1, gamma=15)\n",
"\n",
"system3=system(potential=totpot, sampler=integrator, start_position=4.2, temperature=1)\n",
"\n",
"#simulate\n",
"cur_state = system3.simulate(sim_steps, withdraw_traj=True, init_system=False)\n",
"\n",
"print(\"Trajectory length: \",len(system3.trajectory))\n",
"print()\n",
"print(\"last_state: \", cur_state)\n",
"print(len(system3.trajectory))\n",
"system3.trajectory.head()"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"#plot\n",
"simulation_analysis_plot(system3, title=\"position Langevin\", limits_coordinate_space=list(range(0,12)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Through the visualization we get the confirmation that we scaled the original potential down to 12.5% of its original height (first panel, compare the y-axis). Accordingly, the system has now enough thermal energy to cross all energy barriers and all four minima can be sampled."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Temperature Replica Exchange / Parallel Tempering"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Temperature replica exchange, also called parallel tempering, is an enhanced sampling method that runs multiple copies of the simulated system at different temperatures. After a specific number of simulation steps, the current coordinate is exchanged with a simulation at a different temperature. However, this exchange is only triggered if a certain condition, e.g. the Metropolis criterion, is fulfilled. This approach has the advantage that one can couple high temperature simulations, that cross energy barriers quickly, with lower temperature systems that sample the minima. Therefore, thermodynamic properties can be calculated with higher precision as more local minima can be explored.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In our example we will perform a temperature replica exchange simulation with two systems. The temperatures of these two systems are defined by the $T\\_values$ parameter. We first define how many simulation steps we run per replica and how often the system should try to exchange the coordinates ($simulation\\_steps\\_total\\_per\\_approach$ and $trials$). We wrap our system with the replica exchange conditions using the $replica\\_exchange.temperatureReplicaExchange$ class. Other replica methods covered in the Ensembler package are Hamiltonian replica exchange and Replica Exchange Enveloping Distribution Sampling (REEDS). REEDS is explained in the Free energy example notebook.\n",
"\n",
"For our example simulation we use a common metropolis monte carlo criterium for RE approaches exchanging the coordinates of replica i and replica j,\n",
"\n",
" $p_{ij} = min(1, e^{(H_{i}(r_j)+H_{j}(r_i))-(H_{i}(r_i)+H_{j}(r_j))})$\n",
"\n",
"Note that every second trial is not accepted, as the algorithm alternates the partner of the pairwise exchange. Therefore every second exchange is a border exchange, which is in this case of two replicas not exchanging at al.=l.\n",
"Subsequently, the simulation is performed as in our previous examples. As the langevine sampler is not calculating any kinetic energy, every exchange is accepted with this criterium."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"DO trials: 10 steps: 200\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"\r",
"Running trials: 0%| | 0/10 [00:00"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"stats= ensemble.exchange_information\n",
"replicas = np.unique(ensemble.exchange_information.replicaID)\n",
"trials = np.unique(ensemble.exchange_information.nExchange)\n",
"\n",
"import itertools as it\n",
"\n",
"fig, ax = plt.subplots(ncols=1, figsize=[16,9])\n",
"\n",
"replica_positions = {}\n",
"for replica in replicas:\n",
" replica_positions.update({replica: stats.loc[stats.replicaID==replica].replicaPositionI})\n",
"\n",
" x = trials\n",
" y = replica_positions[replica]\n",
"\n",
" ax.plot(x,y , label=\"replica_\"+str(replica))\n",
"\n",
"#plt.yticks(replicas+1, reversed(replicas+1))\n",
"ax.set_yticks(ticks=replicas)\n",
"ax.set_yticklabels(labels=replicas)\n",
"\n",
"ax.set_ylabel(\"replica Positions\")\n",
"ax.set_xlabel(\"trials\")\n",
"ax.set_title(\"Replica exchange transitions\")\n",
"if(len(replicas)<10): plt.legend()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now visualize which part of the potential energy surface our simulations sampled. For that reason, we color-code by the temperature of the replica.\n",
"We start both replicas in the global minimum at $r=4$. One replica is at the high temperature (T=4, red), one at the low temperature (T=0.1, blue). The replica at low temperature intensely samples this minimum, whereas the replica at high temperature is able to overcome the energy barriers surrounding the minimum. Therefore, only states from the high temperature replica are observed at the energy barriers.\n",
"\n",
"After 200 steps the coordinates between our two simulations are exchanged. If the high temperature replica crossed the energy barrier before, the low temperature replica inherits these coordinates and starts in a different minimum. This minimum is then intensely sampled. Using multiple exchange trials, the low temperature replicate can sample all four minima without the need to cross the energy barriers in the low temperature replica.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"C:\\Users\\benja\\anaconda3\\envs\\EnsemblerDev\\lib\\site-packages\\ipykernel_launcher.py:24: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbYAAAEyCAYAAACWKPW+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydeXhV1dW435WZkIFAAiIYwiAK6DUqWnAoiCBVHNOo/Tkhah1bRUFr/Sqg0pZWsGq1SKlT1X62pKhUbNUiggP6FS1EBlEEEkEZhBCGQMb1++OcJDfJvTc3yU3uTbLe5znPvefsaZ1z9jnr7L3XXltUFcMwDMPoKESFWwDDMAzDCCWm2AzDMIwOhSk2wzAMo0Nhis0wDMPoUJhiMwzDMDoUptgMwzCMDkVYFJuIaEu2cMhsGB0FEdnrPkuTfYQ96oatCodsnQkRudjrvdYt3PJ0JKzFZoSUQC9NwzCMtiBcii3Zz3aeV5ybA8QzDMMwDJ/EhKNQVT3g67iIHPLaLfUXzzCM1kFVJwPW2jbaNdYVaRiGYXQoOoRiE5EkEblHRD4Ukd0iUioiW0Xkf0XkewHSveqOB73q7p8lIq+LyA4ROSgiq0TkRhERrzQDROQpEdkkIodFpFBEficiqX7KyPYaIM4WkTQRmSUiX4jIIRHZJSKviMiINjzPc0RkkYh8KyKVIvKcV9xuInKFiLwsIhtdGQ+JyGYRedGfnNVlANXX4Xc+DH+yveI3OhYnIpPdOHt9hNUZeBeRPq7hw5ciUuLLyEhE4kTkFhFZIiI7RaRMRLa7sv/AnxyN0Zay1DfuEJERIvJ3914edsv8bXONEernHyDeBSLyNxEpcMv9zn1eHheRkT7i9xKRG9zzq05z0H0O/igiQ5sjr49y+onIIyKyRkT2udf/cxF5TET6+ogvIvJv95x3ikgvP/ke78qsIjLPT5x4EblVRN5272WpiGxzn9dfiMiARmRPF5HZ7nN3WJx3w9+9nxsfaRLd+vesiKx3r2n1e2FhY/W6/nMoIhNdefe6ef1H6r0D/eRznIj8xT3vw+K8Hx8XkSPdZ6L6+bg4QB5NuncBUdWI2YDRgLrbtUGmORnY5pXO1/Y/ftK+6oa/CtwBVPlJ/5gb/0ygyE+cT4AEH2Vke8U5D/jKT/pK4MY2OM9f+Ej3nFfcpY2UUQXcE6CMQFu2V/y97rHJAc55shtnr4+wi73yHQ3sql9evfgDgPWNyPcUENWMettmsgCPuuGrgKuAcj/pC4D+fuT1e+298/eTthvwr0Zk93W/NjeSpgy4soXvj2uBwwHK2A+M95GuD7DbjfOGj/B4IN8N/wLo6iPOMPw/29XbqwHqzXCg0E+6Q8D3/Zzzs42UqcCTAa5ZdV24E3ipmXlcApT6SbcdOM1r/+JQ3ju/MrWkIoV6o4mKDejvdWM2AFcD/YA0HEXwZ6/8rvCRvvplXICjWP4KnAJ0B06g7gN8PrADR4GdC2S45c/2ivNzH2V4K7bNwAHgLlfOnm6l2OiGVwIjW/E8t+EopsU4SjodONr7oQFexnmpXgB43PPMAs4BXqFWuY2uV0YCkOQl5z3uvvcW5RU/lIptq7tdC/QFjgAu8YqbBmxx434D3AIMco8fBzxC7UfNfc2ot20mC7WKZwdwEPgvMN69T4OA6ThKQnFexrEBXmZNUmw4Y/LLvc71BeD7btk9gdOBh4C1PtK+C8xx69FQoAcwELjIDat+gQ9u5rsjx0uut4AfuNc+3b0+H7lhB3yVAVzqlf4n9cJ+5x4vB07xkbYPsNONUwLMxHl/dHfrwHnAPOClAPVmk1svrgSOdK/nNUAxtc+9+Cj7dzjPfy5wItDLLfNM4BmvunSNn+u216v8CmCWe3+6A98D3veS8XQf6QdTq5C2ufIf4Z7DtTj11FvhN1BsLb13Ps+rOZWotTaartjecONuBFL9xKl+UAuBmHph3q2MZ32kTcT54qiu1PlAoo94r7txPvcR5q3YFN9fjL3dCqDAB618nq/4ekCacI/+4ObzeiMPil+FFWw8gldse4GsIGTeDWQ2UtZBoHsTr0mbyeJ1nxX4HEj2kf46rzg/bsq1J7Bi+4lXvncGOMcYf2F+4ovXM/REM+pkIrWt5L/5iRMHfOzG+V8/cZ6jVsEOdY+NpVY5/MJPujw3/DA+Pkz9XZd69WYncKSPNFd5xTmtGdfmHjftGj/he73yv8FHeCq1rdk/+ghfQG2LaqCP8GHUbYldXC88JPeuQZqmXqjW3GiCYsP52quucD8IEK8btd01Z9ULq37hlwE9/aSf7++m+Kl8qfXCvBVbg24Or3h3esU7upXOU4FBLbxH1d0KB4BoH+HhUGwPBcgjCecrWoGbA8SLpvYjZlITr0mbyUJdxXZZgDxW4f9DqbmKbYO/PFu6AVe4ea9vRtrr3bSl9Z+/evHOp1YBxfsIT8Zpuah7/XpT2/3/gZ/63genp0WBX7ag3kzxEyeeWsVwRzOuzZFeZfQKUBc+C5DHs26c/9Y7nkJt78BvAqR/wkuG+ootJPeu/taejUfG4HzplQMrxTGsaLDhNK+3uGlO9pPXalXd6Sdsk9f/f/uJ85XXf5+Dzy6vBgh7xev/6V7/Q3meX6rqxgAyACAix4pjEPOpO4hcKbVeXz5wo3XFeWgigX8GCBsJdHH/fxDg+nUB1rrx/F2/SJKlEqeV44/qujZcROKCFd4fItIHp9sJnC7I5uQx3DUUWesaB1R51auX3GhHNyPrMe7vSqAywHX93I0XDxxfPxNV3Y/TzV+J05W4GqeO7weuUtVKH2WPptYIr1nXxeVNXwdVtRSnFwb8vFtEpK+IzBSRj0Rkj4hUeF3XbV5RB/tKH6h8ly/9lD8ciHX//yNA+kUBwkJy7+oTlnlsIeIY9zcWpykbDBl+jm8PkKZ6bt0B9T+vznv+XRc/caD25jRAVbeIyGGcsap+XkGhPM/NjSUUkRuAJ3Ga/42RCnwdpEytSaDzOsbrf36Q+fm7fpEkS6GqlgRIV13X4nBaHgVBluePgV7/Vzc1sYhMxxn/C2hdB0SLSFKAZ80X1df1NBwlFAw+r6uqfiAiv8YxsqqO81NV9Xdfq69LKU6Ltrl8GyCs+j43eLeIyASccfGkIMrwabndgvK931NfBEgf6LqE7N55055bbIFukj/i/Rz39SXWnDgQ+ME92Eja6nDvShrK8zzk5zgAInICziB3HLAOuAHHgKQnTrdDMnCGV5JI+TAKdF6hvH7B0FayNFaXvBVDMC+9xvD2+BPsCwgAETkfmIHzbKzAMTAYivOCqvYmdIVXkqbWq1Df4/94/T8AvBYgbvV1OaBun1kzCeb9UufdIiJHUKvUvsax7D4Zx/Ci+nk9yitJoOva5PJxem2qCVQfA32ktMrzGSkvpuZQfSELVDUrnII0ga5BhntXhLY8zx/jfOxsxxkE31c/gojENkjVPIJ5CYSiflZfvwqgi6pWhCDPSJClsbrkrcxC4cHHW5k11a3dze5vPo4FboPzFpGWfExUX9fnVfXaFuSDiGQAf/Q6lATMBf6fnyTV1yVJRKSFyq2pXIUj32HgDFUtrB9BRNJbsXxvZdYV/8ot0IdVyO6dN+25xVY99nWUiPQIqyTBc6y/ABHJwumGhLrdRm15nh7391++lJrLcSEq67D7G6jrtncIyqm+fjEE0TffyoRSlkwRSQwQXl3XygjczRQs3mOzHr+xfFMdf2EAZd6SelV9XU9sQR7V/AlnLGkv8KB77EcicqWf+NXXJZ7AY1itQfV1/ciXUnMJ1fPqC+/3VKBzDxQWyntXQ3tWbNWGHFE48z3aA35n3dcL+8Drf1ueZ/VXc3SAOFcECAPHyKWxPKB2XDOQscC4RvIIhuU4L3dw5tWEk1DKEo1jKeaP6vq0UlXLAsQLClX9htpxu6uamDxgvXJ7AS5tpmgAb7u/HhFp9gtSRG4ELnR3b1XV6dQaPjwpIpk+kr2LY7UMjuFJWxKK57UlrKT2eQ9UFy8MEBaSe1efdqvYVHUdtZY8MwK5nQEQkaMbcwvTBpwrIuPrHxSR3sDP3d0PVbXaCqmtz7N6gPws1xKpft434Vj2BWK3+9tYa6t6HOMCEWnQtSUikwhBC0tV9+KYKwPcKiLnBIrvugAKxZhUW8jyoJ9rdx2OVR84c7NCxRPu7xki8lN/kUSkfhdydb06T0R8vYQfBHwpjWB5gdp696yIpAWKLCINWhAicjTO5HiAv6jq/7r/b8CZY5oK/FlE6rwzXYW/0N2dIiKnBig31EM/1df1JNdqtX55PwB+FOIya3B7daoV/63iw2WYiAzBuYb+aPG98ydcxGw0z/NI9eTBAzgz/k/C8WqQgfNw34AzwbkS/xOXXw1Qht+5VF5xvOeqZQcI2+zKeSfOg5yB43nkSzc8kOeRVj1PN94PvWR9F8dQJB1nkuVsN+91/s7VzeMvbtjXON4PknC63urLdLpXPsuBU3G8bwwBfoUzDlXtsaCxeWzdGjmvNGq9u5QDj+Mo6Az3Gg7FaYUswLFua9Jcv7aUhYaeRz7F8c6QjmOhN43W8zwSTV3PI89T68EmAxiBY/n4Wb10U7zS/B3HwKGHW4er50h516uA19DPdb2I2vmeW4BbcbrAuuF8ZJ2GM1n5U+D9emljqJ0AXEDDuajnecnmy51cH2onGR8AHsDpJkxzw8bjTMx/sTn1hto5iY/WOz7cK/1qt5yeOB5ofoFjyOR9XX15/WjpfNJjqZ1n9zXOWGQv95pPxOkGb8zzSLPvnV+Zm1qBWnOjeb4ij6f2RRFoK6PeBEvaXrGdR+0E0PpbY74iW/U83XhCrScBX9t/gbP8naubx+nUTlitv9W/Nk8GKOsFaiett0ixufGPwmklNnb9KgngOcRP3m0mC3V9RV6N8wHgK20hrecr8t+NyL23Xpp4nO51f/HfAi5ryjX0I9ulOMYcjV3XJfXSPeh1vUf5ybvaY0ypn3p/PLWu0vxtgXxFNlmxuWGzA5S3mbrvn5ArNjc8l9qPqfrbDpwP5Or980N57/xt7bYrshpV/QznK/dGHN+O23G+hA/hKJGFOF8OPdX3BMu25Bucr6w5OErqMPAdjjnx6ar6R38J2+I81alhP8KpyPk4D/E+nAfrPpyWRVEjeXyAMzb2Bo6boECy/MQ9n5U4c2X24ZiDT1TVq3EqckhQ1a9xfN/9CGcy/Dac8yvF+dJ8A+dL8UhV3RKqcltTFlV9ARiFU392uum/wnnZnaD+5161RPa9OPf3UpxuqG9xXmq7cD58HsPxB+mdphTHNdUDOHOaynDq0UfAbTi+AUMxDrgAx8H0dDfvPTj1bz/OhPengQk4vl4BEGclgvvc3dmqusxP9lNd2eOAl0QkwTvQfT6H4PiBfQ+nh6UMx2/oB24Zd7T0HOujqlNxxt4/xnmGDuI42P41Tot4S6jL9CFDHk4r/K849bDMLXeue9zbyMTnVJHm3LtAiKstjVbCHRP7r7t7oqoGXA7EMAIhIo/ivCBXq2rA8VbDiAREZDTOqiHgdK1/FSB6SGj3LTbDMAwjornA/d1LXReFrUZEKzYR+bk4C93tE2fRvX+IyHH14oiIzBCRb8RZEPNdERkWLpkNwzA6EyLSPUDY8ThLMwHkaRt1EUa0YsMxJvkDjlXMGJyB8n/Xu5D34Fhd/RRnLbWdwNu+zKANwzCMkDNfRBaIs5L3UeKsmD1ERO7GsaLtgmML8Ju2EiiiXWqpap05XyJyNc7Ce6cD/3Dna00GZqnq3904E3GU2xU4fg8NwzCM1iMaZ7HQXD/hh4AfaRAri4SKdmU84k5k/gY4U1XfdycEfgWcqqr/8Yq3GPhOVSf6yONGHEs8unbtevKxx/r1chUSSkpKWL9+PQBDhgwhMTGQFyTDCMzXX3/Nzp076dKlC0OHDg23OEYn5ZNPPvlOVTPAWY4Ix0p2FM7q3Rk4FroFOJ5FHlX/Lr9ahYhusfngMRzT8xXu/hHu74568XbgTIxsgGtS/0eA4cOH68qVK1tBTMMwjI6LiNSY8KvqSpwpOxFDu1FsIvIIzkS/M3zM06rf7BQfxwzDMIxOQKQbjwAgIr/DcdUyRlW9zUWrHekeUS9JTxq24gzDMIxOQMQrNhF5DMcQZIyq1l+BejOOchvnFT8Bx3/dh20mpGEYhhExRHRXpIg8ieMP72KgyF0xFpzVag+oqrqeGP5HRD7HWZ78FziOSP8SFqENwzCMsBLRig3HXx7AknrHH8BZah7gtzjzJJ7E8ab9MXCOqjZp+XrDMAyjYxDRik1VG11XzJ3JPoNaRWcYhmF0YiJasRlGpFFaWsqePXvYv38/lZXhXizCMEJHdHQ0ycnJdO/enfj4+MYTRDCm2AwjSEpLSyksLCQtLY2srCxiY2MJ/6LshtFyVJXy8nL27dtHYWEhmZmZ7Vq5RbxVpGFECnv27CEtLY309HTi4uJMqRkdBhEhLi6O9PR00tLS2LNnT7hFahGm2AwjSPbv309KSkq4xTCMViUlJYX9+9u37Z0pNsMIksrKSmJjY8MthmG0KrGxse1+/NgUm2E0Aet+NDo6HaGOm2IzDMMwOhSm2AzDMIwOhSk2wzDahGuvvZasrKya/S1btiAiPPfcc61e9nPPPYeINLrNmDGj2WVs3bqVn/70p4wcOZLExEREhC1btviMW1RUxA033EB6ejpdu3Zl7NixfPbZZw3iHT58mLvvvpvevXvTpUsXRo4cyfLly4OWaf78+Rx77LHEx8dzzDHH8NRTTzX39NoVNo/NMIyw0Lt3b1asWMHAgQNbvawJEyawYsWKmv1PP/2U2267jccff5xTTjml5njfvn2bXcbGjRv529/+xsknn8yZZ57JW2+95TOeqnLhhReyefNmfv/735OWlsavf/1rzjrrLFatWlVHhuuvv57Fixfz8MMPM2DAAJ588knGjx/PihUryM7ODijP/Pnzuemmm/j5z3/O2LFjWbJkCbfeeiuqyi233NLs82wXqGqn3U4++WQ1jGBZt25duEVoEw4fPtwq+U6cOFH79evXKnk3laVLlyqgb7/9dsjyrKysrPk/f/58BXTz5s0N4r366qsK6DvvvFNzbO/evZqWlqY//elPa46tWrVKAX3mmWdqjpWXl+vgwYP1ggsuCChLeXm5ZmRk6DXXXFPn+KRJk7RHjx5aVlYWMH1jdR1YqRHwDve3de6uyKKicEtgGGFlxowZiAhr1qxh/PjxJCUlcdlll9WEL1y4kBEjRpCYmEi3bt249NJLKSwsrJNHVlYWV111FfPnz2fQoEEkJCRw0kknsXTp0oBl++uKXLZsGePGjSM1NZWuXbtywgkn8PTTT9eEv/zyy4wZM4aMjAySkpI48cQTef7551t+MVpIVFRwr9NFixZx5JFHctZZZ9UcS01N5YILLuC1116rEy82NpbLL7+85lhMTAw/+tGPePPNNyktLfVbxooVK9i1axdXXXVVneNXX301u3fv5v333w/2tNolnVuxffNNuCUwOjv5+TBjBlx3nfObnx8WMS666CJGjRrFokWLuPPOOwF46qmn+OEPf8jQoUPJy8tj3rx5rFmzhlGjRjWYwLts2TIeeeQRfvnLX/Lyyy8THx/Pueeey4YNG5okx2uvvcbZZ59NWVkZ8+bN47XXXuO6666joKCgJs6mTZvIzc3lpZde4tVXX+WCCy7ghhtuCOn4UWVlJRUVFY1uTuOlaaxdu5bjjjuuwfFhw4ZRWFjIgQMHauL179+fxMTEBvHKysrYuHFjwDKABuUMGzYMgHXr1jVZ7vZE5x5jC/DFYxitTn4+zJ4NaWnQt6/TgzB7NkydCh5Pm4py++23c8cdd9TsHzhwgJ/97GdMmjSJZ555pub49773PQYPHszTTz/N5MmTa47v2LGDDz74gMzMTADOPvts+vXrx8yZM3nhhReCkkFVueOOO8jOzmbp0qU1LaCxY8fWiXfffffV/K+qqmL06NF8++23zJ07l5tvvrnpJ++Ds88+m2XLljUab+nSpYwePbpJee/Zs6eOEU013bt3BxzDkqSkpBoXbv7iBXJ7VR1WP30waTsCnVuxNeNryzBCxsKFjlKrfvlU/y5c2OaK7ZJLLqmzv2LFCvbt28eVV15JRUVFzfG+ffty7LHHsnz58jqKbcSIETVKDSA5ObmBwUZjbNiwgYKCAu69996A3Xpffvkl06ZNY/ny5Wzfvp2qqiqAkDrtnTdvXlBupY455pgm562qPidB12/9BRvPXxnQMSZbN4fOrdhEIC8PcnPDLYnRGSksdFpq3qSmOsfbmN69e9fZ37lzJ9CwtVRN/ZZAr169GsTp1asX27ZtC1qG3bt3A4EtEw8cOMC4ceNITExk1qxZDBw4kLi4OObOnVunZdlSBg0aFJQCiY6ObnLe3bt399liKnLH/Kuvbffu3RuMZ3rHq259+SsDnJaZ972tLjdQ2o5A51Zs0dHwxBOm2IzwkJnpdD96K4niYud4G1P/y75Hjx6AM/+relzGm+Tk5Dr7O3bsaBBnx44d9OnTJ2gZ0tPTAQIqwxUrVlBQUMB7773HGWecUXPcu1UZClqzK3LYsGE+pwKsW7eOzMxMkpKSauK98sorlJSU1BlnW7duHXFxcQwaNChgGeCMtXkrtuqxtaFDhzZJ5vZG51ZsAJs2hVsCo7OSk+OMqYHTUisudhTd9deHVy7gtNNOIzk5mY0bNzJx4sRG43/00Ud8/fXXHHXUUYCzEsLixYuZMGFC0GUOHjyYrKws/vSnP3HjjTf67EYrKSkBqOOMuqioqI41YShoza7ICy+8kGeffZZly5YxatQoAPbt28c//vEPrrjiijrxpk+fzoIFC2ruQUVFBX/9618555xzAna9jhw5kvT0dF566aU6re4XX3yR7t27c/rppzdZ7vZE51ZsFRWwb1+4pTA6Kx6PYyiycKHT/ZiZ6Si1Nh5f80VKSgoPP/wwt912G7t27eLcc88lNTWVbdu2sWzZMkaPHl3nJdyrVy/OOeccZsyYQXx8PL/5zW84ePAg999/f9BligiPPvooOTk5jBkzhptvvpmMjAzWr1/Pzp07eeCBBzjttNNISUnhtttu44EHHuDgwYPMnDmT9PR0iouLQ3b+zVFYAHl5eQB88sknAPzzn/8kIyODjIyMGiV24YUXMnLkSK666ioefvjhmgnaqso999xTk1d2djaXX345kydPpry8nP79+zN37lw2b97MSy+9VKfcQYMG0a9fP5YsWQI4iv+hhx7i1ltvpU+fPowdO5Z33nmHZ555ht///vfExcU16/zaDeGeSBfOrTe9dTXHq86e7WMKomHUpSNO0J4+fboCWl5e7jN88eLFOnr0aE1OTtaEhAQdOHCgTpo0SdeuXVsTp1+/fnrllVfq/PnzdcCAARoXF6fZ2dm6ZMmSOnnVn6C9efNmBfTZZ5+tE2/JkiU6evRo7dq1q3bt2lU9Hk+dScpLlizR7OxsTUhI0AEDBuhjjz1Wcx7B0hoTtFVVAZ/bqFGj6sTbvXu3Tpo0SdPS0rRLly46ZswYXbVqVYP8SkpK9M4779RevXppfHy8nnrqqbp06dIG8fr169egDFXVp556So8++miNi4vTQYMG6ZNPPhnUebT3CdqindgysJdk6vioXzG15wt4vn0z3OIYEc769esZMmRIuMWIOLKysjjjjDN48cUXwy2KESIaq+si8omqDm9DkZpEp56gHU0lafGHWLhnVLhFMQzDMEJEp1ZsAKmHtlMo/cIthmEYhhEiOrfxCEoxKWTK1zafzTCaib+lWQwjXHTqFlslMRTF9iSn5wfOfDbDMAyj3dOpW2zRVDI15lE8+7+ETUnhFscwDMMIAZ1asR3JNxxfvhNKqxwvJIZhGEa7p1N3RQLOJO2qKti7N2xLhhiGYRihwxRbNVVV8NBD4ZbCMAzDaCGduisSQGNiEHC6It97L9ziGIZhGC3EWmyVlWhVlbOEzaFD4ZbGMAzDaCGdWrEdogszmM518gwzyu4jn+NtnM0wWolrr722zsrRW7ZsQUR47rnnWr3s5557DhFpdJsxY0aLyvn666/Jzc0lNTWVlJQUcnJyfK6p5ovDhw9z991307t3b7p06cLIkSNZvnx5g3iPPPIIF1xwAb179w6JzB2RTt0VuZ1eFGkqfSu3UEQ3ZldMZupDr+JZEH7v6obR0enduzcrVqxg4MCBrV5W/dW8P/30U2677TYef/xxTjnllJrjgRY5bYySkhLGjBlDfHw8zz//PCLCL37xC8466yzy8/Pp2rVrwPTXX389ixcv5uGHH2bAgAE8+eSTjB8/nhUrVpCdnV0Tb/78+aSkpHDxxRfz1FNPNVvejkynVmwxVNKdYlSEbnGHoHI3C9/siic/PyKWDjGMcFBaWhpwra9QER8fz4gRI1q9HKBm6ZhqDh8+DMCQIUNCJsP8+fPZtGkTGzZsqFkE1OPxcPTRRzNv3jzuuusuv2lXr17NX/7yF5555hkmTZoEwKhRoxg2bBjTpk1j0aJFNXHXrl1LVFQUFRUVptj80Km7IqOjFK1el6isjNTSnRSW9nLWxzKMTsCMGTMQEdasWcP48eNJSkrisssuqwlfuHAhI0aMIDExkW7dunHppZc26FrLysriqquuYv78+QwaNIiEhAROOukkli5dGrBsf12Ry5YtY9y4caSmptK1a1dOOOEEnn766Zrwl19+mTFjxpCRkUFSUhInnngizz//fMsvRgtZtGgRI0aMqLOydf/+/Tn99NMbXQh10aJFxMbGcvnll9cci4mJ4Uc/+hFvvvkmpaWlNcejojr1azsoOvUVqqwSpKwM3KV7ikkhs/wrCGJJeMMIBfn5MGMGXHed8xuuId6LLrqIUaNGsWjRIu68804AnnrqKX74wx8ydOhQ8vLymDdvHmvWrGHUqFENVpdetmwZjzzyCL/85S95+eWXiY+P59xzz2XDhg1NkuO1117j7LPPpqysjHnz5vHaa69x3XXXUVBQUBNn06ZN5Obm8tJLL/Hqq69ywQUXcMMNN4S09VJZWUlFRUWjm/eyX2vXruW4445rkNewYcNYt25dwPLWrl1L//79SUxMbNWK4hgAACAASURBVJC2rKyMjRs3hubEOgmduiuyghj2kkYKe9mrqewljevlOdiwI9yiGZ2A/HyYPRvS0qBvXygqcvanTm37nvDbb7+dO+64o2b/wIED/OxnP2PSpEk888wzNce/973vMXjwYJ5++mkmT55cc3zHjh188MEHZGZmAnD22WfTr18/Zs6cyQsvvBCUDKrKHXfcQXZ2NkuXLq1pmYwdO7ZOvPvuu6/mf1VVFaNHj+bbb79l7ty53HzzzU0/eR+cffbZLAviA3fp0qWMHj0agD179pCWltYgTvfu3SkqKgqYT6C01eFG8HRqxXaE7CBNiimo6ksmBVzPn/DoWtgV7bx1bJzNaEUWLnSUWvX7rPp34cK2r3qXXHJJnf0VK1awb98+rrzySioqKmqO9+3bl2OPPZbly5fXUWwjRoyoUWoAycnJDQw2GmPDhg0UFBRw7733Buxu+/LLL5k2bRrLly9n+/btVFVVAYR0XHDevHkNWqW+OOaYY+rsi0iDOMEs5qyqzU5rNKRTK7YuXWBG6QMolbUHFaisdLyQLFgQNtmMjk9hodNS8yY11Tne1vTu3bvO/s6dO4GGraVq6rcuevXq1SBOr1692LZtW9Ay7N69GwhsmXjgwAHGjRtHYmIis2bNYuDAgcTFxTF37tw6LcuWMmjQoKCUSrSXj9m0tDSfLauioiKfrTFvunfv7nNaQHVLr7rlZgRHp1ZsHHEEbN4MOPqs5ntJBN55J1xSGZ2EzEyn+9H7nVdc7Bxva+q3Fnr06AE487+GDRvWIH5ycnKd/R07Gnbf79ixgz59+gQtQ3p6OkBAZbhixQoKCgp47733OOOMM2qOe7cqQ0FzuiKHDRvG2rVrG8RZt24dQ4cODZjPsGHDeOWVVygpKakzzrZu3Tri4uLqGKQYjdO5FVtaGhQUIJVeLbboaMeYJIhuCMNoCTk5zpgaOC214mJH0V1/fXjlAjjttNNITk5m48aNTJw4sdH4H330EV9//TVHHXUUAPv372fx4sVMmDAh6DIHDx5MVlYWf/rTn7jxxht9ds2VlJQAEBsbW3OsqKioUavDptKcrsgLL7yQqVOnsmnTJgYMGAA4lp8ffPABs2bNCpjPhRdeyPTp01mwYEHN9a6oqOCvf/0r55xzTptMv+hIdG7FBhAbC6rkVx3HQi6msLIf/SggR/6BjbAZrYnH4xiKLFzodD9mZjpKLRKGdlNSUnj44Ye57bbb2LVrF+eeey6pqals27aNZcuWMXr0aK644oqa+L169eKcc85hxowZxMfH85vf/IaDBw9y//33B12miPDoo4+Sk5PDmDFjuPnmm8nIyGD9+vXs3LmTBx54gNNOO42UlBRuu+02HnjgAQ4ePMjMmTNJT0+nuLg4ZOdff+wsGH784x/zxBNPcNFFFzFz5kxEhPvvv5+jjjqKm266qSZeQUEBAwcOZNq0aUybNg2A7OxsLr/8ciZPnkx5eTn9+/dn7ty5bN68mZdeeqlOOStXrmTLli01Y4vr1q0jLy8PgPPOO6+BZWWnRFU77XbyySer9uunq/Ho1TynP+V3+gtm6O38Tq/meV09+001jGrWrVsXbhFCzvTp0xXQ8vJyn+GLFy/W0aNHa3JysiYkJOjAgQN10qRJunbt2po4/fr10yuvvFLnz5+vAwYM0Li4OM3OztYlS5bUyWvixInar1+/mv3NmzcroM8++2ydeEuWLNHRo0dr165dtWvXrurxePSZZ56pE56dna0JCQk6YMAAfeyxx2rOI1iWLl2qgL799ttBpwmGgoICzcnJ0eTkZE1KStKLLrpIN2/eXCdO9XlPnz69zvGSkhK98847tVevXhofH6+nnnqqLl26tEEZEydOVJzRkwZb/bKaS2N1HVipEfAO97eJRrjVjYh8H5gKnAwcCUxS1ee8wp8D6veVfKyqjboTGD58uK485RRm/LE3e6pSSGOvkydCUWwGaUd2YcaWSSE6E6O9s379eoYMGRJuMSKOrKwszjjjDF588cVwi2KEiMbquoh8oqrD21CkJtEeJmgnAWuAOwB/7vf/DfT22s4LOvdbbqGQfnSjtj9dgdTKIgp3WL+2YRhGeyPiFZuqvqGq96lqHlDlJ1qpqm732oKfzejxkNllF8VRadTaRSrFVV3JLNto3v4NwzDaGRGv2ILkDBHZKSJfiMh8EenZlMQ5x66jqCqFIlKpQiiiG0WkkVO1wBndNwzDL1u2bLFuSCOi6AiK7V/ANcDZwBTgVOAdEfHZjygiN4rIShFZuWvXLgA852cyVR6hO0Vs5SjSKGKq/A5P9Hr48MO2Og/DMAwjBLR7c39Vfdlr9zMR+QQoACYADdz0q+ofgT+CYzwCQE4OngcfxBOzlioRtKqKKBGIjYPy8tY/CaPdoH5cHxlGRyHSDQqDoSO02Oqgqt8AW4Gjg07k8UD37lBVhbhzQzQmBqqqnOOGAcTFxXHokD/7JcPoGBw6dKjdTwjvcIpNRNKBPsC3TUr4859DXBy4Ck0OH4bSUoiKAnfyo9G5SU9PZ+vWrezZs4fy8vIO8WVrGOC00srLy9mzZw9bt26tcanWXon4rkgRSQKqHaVFAZkikg3scbcZwN9xFFkW8GtgJ/BKkwqaMgW++QZ54gk+02G8Qg6FMQPpt/sbcm59yvFCkpsbgjMy2iupqanEx8eza9cudu/eHXL/hIYRTmJiYkhISCAzM5OEhIRwi9MiIl6xAcMB76V4H3C354FbgONxjEe64Si3pcBlqtp0Z4/JyeSnnMGc/T8mmWL6xO6mqKoHs/fdyNRZi/GYYuv0JCQk1PhDNAwjMol4xaaq7+LleN8H40NWWGEhCw+MIy32IEmVxaBRpMUcgCpl4aYTzXekYRhGO6DDjbG1iMxMCiWLbpVFREVFOX7HqqpIlX0UxtmyEYZhGO0BU2ze5OSQmVFCcUUiUaWlRJeUwIEDFB+OJ7N/dOPpDcMwjLBjis0bj4ecOWdQlNKPvRVJKDheSKK6k7PmQZgzJ9wSGoZhGI1giq0entzBTM3Ko3v3KL6OzqKb7GWKzsZz+P/gwQfNd6RhGEaEE/HGI+HAU/w+HhaxqrIPr3AJj+pP6VdVQM6+V/A89BAsWBBuEQ3DMAw/WIvNF336kL+nL48whSLS6MvXFJHGbKaQ/9b2cEtnGIZhBMAUmy9+8hMWcjFpFJHGXqJQ0thLGkUs3D823NIZhmEYATDF5ovcXArpRyrFNYcUSKWYQj3KxtkMwzAiGFNsfsiM+ZZiUuscKyaVTApg7twwSWUYhmE0hik2P+QcnU8R3dlLN76lF/9iPIuZwA6OIH+hraxtGIYRqZhi84PnwUuZ2uUPlBLPO4xFgTH8mzipYHbRdeQ/9Gq4RTQMwzB8YOb+/sjNxVNQwBFTd3Iei0ljLwBCFGgUC9/pZr4jDcMwIhBrsQViyhQKYwfQjX0oAgioklqxh8L9aeGWzjAMw/CBKbZGyIz5hmK6ucsLKAoUk0xm5WYbZzMMw4hAGu2KFJHvt1LZZar6USvlHTJyen7A7MJLUa0ilWKKSWUvaVxf9TRMeQ/efjvcIhqGYRheBDPG9i7ONK5Qsx3o0wr5hhTPD45k6rzZLORi/ks2e+lGGntZSA78+xU8eXm2srZhGEYEEazxSBWwNYTl9iPw4qGRw6234pk/H6qq2MQA+lFAKsU1LramTltgK2sbhmFEEMEqtl2q2j9UhYpIVajyanU8HujShYUHc+jmutgCan4XbjjOrCMNwzAiCDMeCYZjj6WQfnTz5WKrqo8ZkRiGYUQQwbTYFgF7Qlxua+TZetx7L5mXfkYRaRwmjs8ZQjHdiOcwJ/IpLFzotOwMwzCMsNNoi01VLwamhbJQVb1YVa8LZZ6tSm4uOYlv8RX9Wc5oDtGFOErZRyrborLIX1UZbgkNwzAMl2C7IgtFZG2rShLheB7M5Si+IYV9lBFPFw4xincZEF3Awi+OD7d4hmEYhktTXGp1blcbU6ZQ+uz/MX7dW4hWAKAxsSAxFBZ0BzP7NwzDiAjMeKQJZOaeSnH6QEjsSmVsnOOFRLqRydcwa1a4xTMMwzAwxdYkcnKg6HACexP7QFQMe6J6sCe6Bzld/glffGHWkYZhGBGAKbYm4PHA1MGLSCvZxmeVw1hTNoR9JTEs3DOK/NJjbAFSwzCMCMAUWxPx3DuBnMoFpFbu4QRWcTz57NFuzK6cbAuQGoZhRABNUWwpInKviJwrIke2mkSRTm4uC5OvIU2LSJNiogS6xR8iLfYgC0t+4MxpMwzDMMJGU6wiuwC/rN4Rkd3AKndb7f6uV9X24y6rmRR2OZa+cZ9AXFcqKysRIJV9FEZlQeEb4RbPMAyjU9MUxVYOfAdUt9bSgbHA2V5xSkVkHbUKbxWQr6r7QiBrxJA5IJqig0eSVraTXeU9+LxyEN/Rk54VReTvzzLfkYZhGGGkKV2Re1S1L5ABjAPuBl4C1gKVON76E4CTgEnAY8AyoEhEvhKRvFAKHk5yftKHoi5H8kXlAD6q+h57NZVoLadPVSGz3zye/DlvhVtEwzCMTktTWmwAqOpuYIm7ASAiccBxwAlAtrt5gFQ3Sn8gq4WyRgye3MFMBW6fVEyFRpPBdxzLeo7Q79hbksbCB8vxjDvC/EcahmGEgSYrNl+oahnwqbvVICJZ1Cq6E0JRVqTgyR3MgJ+/x/f3/B9StJvtmsG7VWdSTDfYp+RMfR7PW3PCLaZhGEanIySKzR+qugXYArzamuWEi8w+FRRtUg5rBis4jQQOE0cpgjJ72SlMzbdGm2EYRltj89haQM5P+lBUlcIqTiSBwwhKKfGcxKekle00y3/DMIwwEKxiew3Y2pqCtEc8uYOZmjiXMuIoJY4EDjGSD+nFTlLZS+GyzeEW0TAMo9MRVFekql7ivS8iJ6jq6tYRqX3hOa8vF+UtoohuHCaO9QzhI0YSTxknfrYSx27GMAzDaCua2xX5qoh0D6kk7ZX77ycn9nW+YoDXIqRl7IvqxraiRPLzvgi3hIZhGJ2K5iq2fkCeiEQHm0BEjmpmWZGNx4PHA0dFbSeFA5SR4CxCWvUOA3QTC6etCreEhmEYnYrmKrbXgFHA48FEFpHTgf80s6zI5957KY3qwvjYd7iI1xjNu/RiB6m6h8L1B2HKlHBLaBiG0WlormK7GvgCuFlEbggUUUQmAf/G8VjSMcnNJfPoeIoruiLAdnryLqPII5dNZJH/6DvOCtuGYRhGq9MsxaaqB4CLgQPAEyIysn4ccXgE+BMQT5Ctu/ZKzoPZFEX34AsGsoLTKaYbMZTTh63MrppM/t0vhFtEwzCMTkGz57Gp6gbgGiAW+Lv3UjYikgL8E7gDx3ny9ap6ZwtljWg8uYOZevRrfEMfKoghmnLiKGctx7GBwfyh4Lxwi2gYhtEpaNEEbVV9DZgJHIFjKRkvIoOBj4FzgJ3AWar6bHPLEJHvi8giEdkmIioi19YLFxGZISLfiMghEXlXRIY1/6yaj+fBSxnAZs5kORXEE0UlKRSjRPFvHWMWkoZhGG1Ao4pNRO4UkbNEJM1XuKpOB94ATgZeBz4CjgE+AYar6ooWypgErMFp/R3yEX4PMAX4KXAKjjJ9W0SSW1hu08nNJTO5iFVkk8AhunAYAQSlR9ReFj60xlbYNgzDaGWCabHNwTH++E5EtojIKyIyXUQuFJFMN84VwEZgDNAN+F/gTFXd1lIBVfUNVb1PVfOAOouYiogAk4FZqvp3VV0DTASSXZnanJzpx7M7qheKoAiHSOAQiWT320Nh+ZG2wrZhGEYrE4xiewP4Bme9tUzgImAa8Aqw2V1J+xVgkxvnRWCiqh5uFYnr0h+nG7RmATRVPQQsB05rg/Ib4JlyDuPOrkJiotkflUZ5VDySEMV72wezqTCa/P9da602wzCMVqRRxaaq56vqUTgrZo8DplJ3gdE04CycMTWAK4EDIrJKRJ4XkbtE5GwRaQ1z/yPc3x31ju/wCquDiNwoIitFZOWuXbtaQSS4ZfYgjhkaw5CBpVSk9EDLqog+tJ8+B79k9leXkD9xjik3wzCMVkJUtfmJG19gFKC6AFXVFi2TIyIHgJ+o6nPu/mnAB0Cmqn7tFe9ZoLeq/iBQfsOHD9eVK1e2RCS/5Od9we23VbCzKJaM8m30ZDs76cUuetKTXTx+/J/w5L/YKmUbhmG0JiLyiaoOD7cc/mipVWSZqn6qqs+q6h2qOkpV04CBQA7wIPAPnJUBWmOJnO3ub/3WWU8atuLaFE/uYAaMPILcqgUcy3q+4FgO0YV0drGTDGZ/No78OW81npFhGIbRJFplPTZV3ayqr6rqDFW9WFX74XRZhprNOMptXPUBEUkAzgQ+bIXymkRmdneKNYUNHEs8h6kkmq8YyB7SnLlt08Oqew3DMDokbbbQqKoWNyediCSJSLaIZOPIm+nuZ6rTj/oocK+I5IjIccBzOB5R/hIq2ZtLTg4UxaSzk56UE00BWZSSQB+2oQj/PjjCWm2GYRghpj2soD0c+K+7dQEecP8/6Ib/FngEeBJYCfQGzlHV/W0val08Hpj6/f/Qk51spR/xHKYnO9hNOpsZwCES+MMMa7UZhmGEkmAmaFeJSIvnozU3T1V9V1XFx3atG65ul2dvVU1wx/nWhFLeluCZM5HHU+4nlWK6sYfdpHOIBIQq0tjDvw9Yq80wDCOUBNtik1YouzXyjDw8HjxPT2Ycb7OX7lQRRQKHSec7iuhOMSncPi3N3G0ZhmGEiGDN76PdhUI7hzIKNbm53NL7B7z97Tj6sI0KoikkC4BMtrCzpCuzb/qCqTjWlIZhGEbzCVaxpQNbQlhu8yfPtVM8j9/AuMuX8mnV8XxDX+I5TCrF7OAIKolmw54e/OH6lTw1+LAzOGcYhmE0i6Z0RYZ661zk5nLL5HiOid5Ed/bUTNg+QFcEZRMD+eu+8eRNeCbckhqGYbRrgnGpFdVK25GNld3R8MyZyNTfZNAzeg9b6YdQCQgHSOIgiZQRy31bbyL//PvCLaphGIZfejjO7iOWoFpsrhd9IwR4ppzD45csJRVnWp+zHGkMh4inhES+ZDBnL76dvDMfC7OkhmEYPsjLoztdMhuPGD6C7YrcKSJ/FpHLRSS18ehGIDz3X8y4qHc4SBIVCGXE4dwK5/vhOzK49f3LTLkZhhF5zJrFLnq3yO9vaxOsYksDrsLx5rFTRN4RkSkickzridaB8Xi45cp9JFBKLJWuJY330KOwi55Mev9a8k78ZevLk5/vGKyI+N9GjLAVCQzDIG/VQPaREtG9eMEqtgwcxfY3HHdVo3E8fqwTkS9F5BF3aZqI1uKRhOfPdzM2+T9UEgNE4xiKKt7K7QDJXL7qbq7pmhd6AfLz4eab4cgj4YQT4LPPAsf/+GMn3siRpuAMo7OSn88TlTcTVXfN54gjKMWmqkWq+hdV/X84nvNH47ix2oDjyX8yzmKf34nIX0XkGhFJbyWZOwy/eGYAmbIVqakk3h9BTguuihheKMnhfFkUGoWSlwc9ezpKat48+PZboFat1t8a8NFHcOqpMGdOy2UxDKN98dBDfEX/cEvRKE32Famqlaq6XFXvVtWh1Cq2d4AE4FLgWeBbEflQRO4TEZuY5QNP7mB+decuktmP1FElDRXcYs6n7wmpzRt3y8+HGTOgf3+49FJwF1gNqMBcfCq60lKYOhWuuabpshiG0W7Jy4MielBJdLhFCUiLFhptkJlIEjAeuAD4AU7rrrqArcDrwGJgiaqWhqzgZtKaC402hbxrXmPKCx6+5ijUy4ikIdWjcUqWFPDbqz4j988X+c/4zDPh/ff95NJ86kiXnQ3//W8LczQMI9LJz/uCqy8toZgUCsmlSj+N2HG2YM39f+kqrYCo6gFV/bvroLg3MBL4NfAZcBRwC87Coz9rtsQdkNw/X8Q/7lrOD1lA4DaU03pTotisWVz6wgWIVPnf3l/K+fwdaLx11pSZ9HXyWLUKzj+/iWdsGEZ7Y+H0VZQTSz8KiSfs7ZKABNsV+XPgCxH5cbBz2lyv+x+r6i9UNRtHsd0KvAGUNE/cjotnzkQWLIjlap4PMkUwzl2iWczFTHCVW6AcGiupPnWU2+LFTh+FYRgdlmXr0vmGI/mY71FOXLjFCUiwiu0gcATwFLBKRMY1Er8BqrpNVZ9S1QtUdXZT03cKcnP58+rhTOB1QudOU/gXF3jt+VFmMTFw8smwYAGoOtvq1XDccTXpAjJlillLGkYHJe/8Z/mEkzhEPOXEUBnhS3kGK91AYB5QCRwP/EtE3hCRIa0mWWfF4+F1vZC7Mv9OLIcIhYKrIsp/yyw6Gu66C8rLYeVKyM2tIwuffeYou4SEBunrSLZ1K8yc2WJZDcOIPGYtdj5wo1HiKY14Z7/BmvvvVNVbgOOARTjvyPHAahH5g4hktKKMnZI5BZeycvb7HEc+jdsuBqbBnBMR6N0bZs+GiorGTfdzc515bFFR/pVbVRW89ppNAzCMjkZeHl8yiCr36a+smXcbuTSpPamqX6jqxcD3gf/gLHtzE/CliNwjIpHd8drO8Ew5h880m9kT3iWRfTSvMinn8g/nb1SUo6RWrYJvvnG6D4MWxgO//S3QsOVXI1VZGUybZl2ShtGByJu2mgN05RCJVBJFNJUdo8VWH1V9X1VHAD8CNgMpONaPn4vI5SGUzwCmvD6Gg9qN1VfPAUppfBp19fFKJvAqr3/vt854WWWl063Y3PXepkxxui1jYvwrt5ISuPLK5uVvGEZEkZ8P931+DTFUUEUU5cRymPgIb6+FYB6b60brNuB/cBYkVeBj4C5V/ajFErYikTKPrSWIlEODyZJRhHB6YkPy8uDaa+HgwToVvI6yW7Cg7nidYRjtjpvP2ciLb/ekCuEQXah9yk9F9ZOIbbi12LRFVStU9TFgEI7/yFJgBPCBiLwsIlktLcPwj2osqlH1tlYuNDcXnnuuoSzeO/fc08pCGI1y/vn+nVqnp9sUDaNR3n1XqCQaAbpSQizlRFHlekqKXEJms6mq+1T1XmAw8ALOe+4yYL2IzBKR5FCVZUQAubkQF+e/r33LFhtrCxd5eRAX58wv9Mfu3Y57NXOLZvgjP59d5SkAlBJPGXHEUkYcZWEWrHFCPhlBVbeq6kTgJOBtIB64G/gy1GUZYWbUqAaHar7jVB1/kkbbcuaZjsIqLw8u/gsvmOcYwydzfvgBe+nGYbpQSRSVRLn/hRjKI7rJFnLFJiI9RORMnO7Iz4FdOB2zNiWgozF7Noj4b7UtXWrdXW3J+efX+Ab1t1qDT5OjxYut5WbUIS8PZm68nBjKa6wgq4ghhnISOUwC+yK62dbs9dNE5AhgCDC03lZ/uZqIHWA0WojHA2PHwttv1zlcs6pcVRXMmmVGJG1Bfn5N12Own9Leq//xwguQkWHzEA0Anpi+EyGKNIooJo1KoqgimiqiSGY/3/Hd1nDLGIigFJuInENd5TUE6FY/mtf/Upy12tYCa7w2o6MxezaMH49s397whVpVBZ9/Hg6pOh+hmGLxyCPOQrL2IdK5yc9n07pUKklhH6lEU+FOyRYqiOF83uApSvaGW8xABGXuLyJV+F4srBLYiKO0vJXYl6paGVpRQ09HMPePCPLy4NJLGyg2qbbA++tf7WXZmrjXHxq21gIvgOQjTlYWbN4cMtGM9kf+0T/kjI3zOUgS6qq0aMrpSgkJlPHWkddxwjf/+kRVh4dbVn80RbFtoW7ray2wXlUjuq81EKbYQsi4cbBkCd71qeaF2aeP40vSaB0GDKhRRn4V1tVXO78vvNAgXoO4rT5fxIhkxskbLGUslcRQOyIbRTRl/IafM2XB6cill0a0YgvWeCRZVQeo6oWqep+q/kVVV7dnpWaEmDlzHBNzr1WNal6P27aZcUJrkZfXuFK76y7485+dbcKEhuH10tq96rzkT3meZZyF00qr7nSLIppykihhytg17aL3JVgnyAdbWxCjnePxwHHHIf6+9l94wea1hZr8fLjuOiCAwcjxx9c1CHn9db/KrYYXXzRr1k7KH+YqlcQCEE0lsVQQQwWxlBNNVbsxLorsRXWM9sW99zY4VOeFaz4kQ0tuLuzf3+BwjcKKi3OUVH1efx0GDmxwuM4cxNtvD5WURjshP+8LFhyaQBVQSSzlxLh1Qiklnv4pu5vvZ7aNMcVmhI7cXJ9L29SwxgxjQ0Z+Pnzp+Dzw2wX5q1/5fxEtXBjYmfW337abr3Oj5eTnw23XH2Iv3aBGnUVRQYw7IVu59+ljwitkEzDFZoSWU08NHG4vy9Bw001AgC7IsWMDL0vk8dS0yvx+iPzqV82VzmhnPHRTIZ/sG+guSCxArVF7NDD66G3tYWitBlNsRmiZNw+o+7Ks8/I1xdZy8vLgo4YLZ9RRUMFc5zlznDG4etTcrz177H51AvLz4d8fJ1NFFPEcQlAEwXF/XMkRUTuZndc/3GI2CVNsRmjxeCA5gL/rXbvMiKSlzJoFBGitHX108GMhL74YuPv4f/7HDEk6OH+46gMOagKlxFNCElWI62dEiaOM84/5qr0MrdVgis0IPa6RiM9WW2ysM75jNI/8fFi/vsHhOoqpKYrI46npPvZ5v0pL4frrmyik0V7Iz4dXPxtABbHU1oAoKl3Dkb7yLbc82DuMEjYPU2xG6LnlFhgyxHdYRYVjlWc0j7lzoaLCv8FIamrTLdfmzYOEBP/h+/YFHq8z2i0PXfU5u0hHa2pR3X6AX121Dk/u4LYXrIWYYjNCj8cDL78M3Rx3otWPioKznMqqVda91VyWLnVavb6IiYH77296nh5PzdI1fsdGH3206fkaEU3enC288dlRVBFN7Z13DEeiqSCd78j980VhlLD5mGIzWgePB0aNQmJjG47fVFbC5MnhkKr9s3cv2rUrZSkpqNS7srNmNb9ldf/9EB3tP7yqyj5GOhjTZkRxmFjqmvRzSAAAGe9JREFUt9IEIZW9HNunJDyChQBTbEbrkZ3tf8HLbdvMiKSp5OVBURG6cyex+/ahsbFIRgb06OEYjLSku9DjgbPOAgK02m65pfn5GxFFft4XfH6gt9tai6L2rlehCAmU85NHB4VRwpZhis1oPXJygAAvSjMiCZ68PLjxRigrc+caQVRZGRQVOR5GxoxpeRmNmfZ/952NtXUQbrquzPXc760CBMcvZCWP3VXQruat1ccUm9F6NGbEUFjYNnJ0BO65x2mteR0ScLp1k5Lg1ltbXobr77Mmb5c6HyNPPtnycoywMufMV/lovy/jLudOXxL7BrlzTmtboUKMKTajdenvTOxs8KJMSYHMzHBI1P7Iz4ctW/yHDx0aOh9+L70UeKyttNS6kNsxeVM+ZNr7Z1O3+9Eb5f6LPmtjqUKPKTajdfntbx0l5otAE7mNWmbOBFXfJv4izlhmqPB44JJL6paBObPuEOTnc+sj/Smhq49AZ9219OhiPPe3T0tIb9q9YhORGSKi9bbt4ZbLcMnNhaefhv79URGqRKiMiXE8yM+caS6bguGf//QflpRUM5YZMhqbMrBmjd23dsiccf9iF73w3VIToqng3jsOtRsP/oFo94rNZQPQ22tr6ADPCB+5uZCZSdTRR1PVpQtUVaElJXDgADz4oHVtBSIvz7lO/pg2LfQvIo/HmehNgFbb7NmhLdNoXebM4aGd1+Pf5bVyR+KzTJnTpy2lajU6imKrUNXtXtuucAtk1GPbNvjuO2IOH0agdg7W/v3w0ENhFS2imTYN8LM0TWJi61kpNtZq27HDPkjaC3PmcP7UgRST5ieCksVXzHk+vU3Fak06imIbICLbRGSziLwsIgP8RRSRG0VkpYis3LXL9F+b0acP7N2LREU5k30rKtCKCqdL8q23wi1dZOK15ppPzjuv9cqeMgXSnRedz1abKkyd2nrlG6EhP58pU6tYzIX4a60JlTx89TratX1/PTqCYvsYuBY4F/gxcATwoYj08BVZVf+oqsNVdXhGRkbbSdnZ+clPnN+Kiob2WPv22ZiNL/7wB6fbVqRGodRct+a6z2oKc+c6rUJ/LF1qrbYIZ8oJb/IIdxGoC/KqzA/bressf4iq38Uv2iUikgRsAmap6iOB4g4fPlxXrlzZNoIZMGAAbN7su1stIQE+/rhDDFyHjN69Yfv2utdLxGktjR0Lb7/d+jLk5cFll/m3yszNhQULWl8Oo8kcI5/zBYMJsJQsqVLM3qpuTc5bRD5R1eEtEK9V6Qgttjqo6gFgLXB0uGUx6vHb30J0tO+urcOHa1aFNnBasPWVGjhKLTGx7Vq4ublw0kmuAFLXoTXAO++0jRxGkzhCvm1EqSlQxf13HmxDqdqODqfYRCQBOBb4NtyyGPXIzYVjjvEf/n//13ayRDqBDGomTGjblu2994IIotrwNdnZV9meMoUp8jCJcoAoqSRKKjgy+lvy5mwJizj5+dBDdrLDr1l/LROyv+0wVpANUNV2vQGzgVFAf+B7wOvAPqBfY2lPPvlkNdqYBQtUnXaHVnlt1cd0wYJwSxh+Vq8OfI1Wr257mbKyVEXqyFMjU1JSp7lvqxds0HGJ72k0pQqVXltVva1Ss/mw9p5FRalefXXryjb7TY1jnw9Z6m8VelevF1tUFrBSI+D9728LuwAtPgF4GfgGKAO2AX8HhgaT1hRbmBDx/9IeMiTc0oWf6dP9Xx+R8Mi0YIFqYqJ/ubKywiNXG3HX1ds1irIAiszXVqkT+HvtNarexo4N+cfJ1YM/ClKuSr0r6pEWlxfpiq3DGY80BTMeCRMjR8JHH1G/5gk41n7+lrrpLPz/9u49Oqry3OP498mde7gIWhW0alBLh7BOXae1UJUKXTb0qBjUHgjYZW1FrJ5mqKW0qK3VZdsM1iqe4lm9mPRim6lHa2hPgaVtQWspLTBiLSmKUFBAwIABcpt5zx97MpnMLSFO9t6z5/msNUuyJ848mUz2b9693/28o0db667FbYodVBo2LPMF24Np3rzYmmxbmcL/Mpc9TGISu5nLU/ga7/HUlHGA4MJnuLlhBscop69De8kMBYTpoiT9/3nuuda55wG+bqFgMzPmjetnfYYKmtlhLhzQc8XTySNKJVq9GkjzZ9jVld8LWgYC6UMN4KqrbC4ozooVUFpKiCmsZCnNnM/rvJ9fcj01NBD8zBrnahsE/uKHuaHhKo4xmlMPNUskuos1cbdedu2yPjAMG3ZK5yqDC59hpBxm6rzz+1mfYQJ7sxJquUBHbDpic8bw4XD8eOod+DnnWH/w+WjcODh8OH2wbdvm7CURs2dz77qP0Mz5vMIHKaWNMto4xigM0FAbwhdY5Fx9WRIoXsaXur4ZXbNsoHpGbJmkjaSzzoI1a2K/b5FO6FVPf8PWUMHL7NgmWXvv6IhNqVSiI4+UU//feCM/R22hUOZQKy52/jq/ujr2MIl9nE0pbYQpZA+T2M/pvMmZPLbyZG5ftB0MEiz+NF/rWoEZ8O6xZ2x2Fc/GFobt67sNEGIKQidCGNm7G5k6BZEIIhGsUJO4W/9qmcDerIZaLtBgU85YscLaUadz11321eIWfbWouvFGe+rIxOdjYsl+3uY0whSyl7PppJgCwpTQzno+Tsj/hNNVDozfT3Dez1jStZI2yujP4b3kWyR6X5gqnqaJ62Lf3VckhZjCVLaSHGCnGmY99U0YeZz9ZmJehRposCmn+Hxw+eVAmlHbrl25/cl/IDZuTH9fQQHU19tXSwZzbz+TYjrZx1kU0QVAmCLGcISxHOKp9SNz73e3cCHBlbtYwipa6O7Eke40jQG6qJq4HWMKEm6F0f8W02TmQk2NtWZeglRR9RRzyd4uOUJN1TvsP5qfax5qsCnn9LX0Sb51ImlrS38Ycto0m4tJzxdYxIqiB+mgmHZKKKKTERzlIONpoZyn+RShS2/NmXALzVnO9Q2fpIafcIhxRBB6Rl6JuqgtesQKrt39GAXV11tNvxsbrXPHKXQH3L+YNLAfoBfDCI6ybVsh9U0p2+XmBQ025RyfL3OT3Zdesq8WpwUC1lVO6SxbZl8t/VD94CXcwC85j9copY13Gcl4DjKKo5TSTt3xzxOa/y2ny+xT0P8i1WsW8muupoMihAiR2KHACPGHGCewj8aa3xDovPPUn6i62joKYYz1gW5scuhMZPd7+EmsUWRN5XaOmfJ8O/KYRINNOauqCsiwoGU+tGsKBmH58vSjtSlT3Hd9mN/P4qq9TKaZ0bQwid0M5SQtlBOmkM18iDu2f5ZQsNnpStMKzHmORSs/yOucRwfF0de/AAMUEqGILgoIU0ob1SVNrG08kZ0u+H4/HDpkhVxtLQwZAsBcniL9SDGVntCtqfgLxhRTvyXPEy1Kp/vrdH9nhUJQWZm+e/wZZ8CbbzpQmI1GjIDW1tQ//5gx1vIwLv0IHiqoZJH5IQahkE7eZSTltFBGG4c4jQ+d/iZLfzfLdeUHpv2Er2y9js6kqfgRiqLhUkCEclpYVfV/VDd9xpa6QoG1TF06k/6NOcIYk2EC1iDS6f5KZeLzwYIFQIpRW2EhHD7sRFX2CQSSQq2X6693bagB+Bb4uJpnuIw/MJzjlNNCmEJe4zyOMJod+0fw2Oe3Ol1mL8Fp9/O1rdfGhVryzMMiuqxQq9lkW6gB+PyzMaYoxaSUVDdnQi0XaLAp59XXWyGWKByGjg5vH45Mce6s12HIxYttK2VA6uuZW/Ib3mE0bzOeTgrZzTm0U8aZ7MMA618a5ppLAEJzlrN869zodH5InkJfQAERzmaPNVLz2AKc+UKDTbnDFVdAcTFSVBTbFBvFLFvmzXDz+60WYpm4eLTWzffTL7OUOsZzkL1MopQ2xnOAw4xjF+/nJGU8trLV+VmSwSD+NdN5jfNJnmzfs9Lcf/A0wemP2jpSU9mlwabcIRCAD3wg9Y6+qwvuucf+mgbbo48CpJ808uEP21nNwFVX46udzfe4g1EcpZwjHGYcJylDiDCaI6xnNqHrvu5YiaFgM7NvHM16PhHr39jzahtAKCDMLNbSWPVTfBtWOVSpygYNNuUOPh88YR2uSjlD8vhx5z/xZ1MwaB1mzSTaLDonBAL46m5iFutoYQwRChDCFBHmAGdYo7adM2HOHNtLCwbeoPoG4bnwDJLfXdbITTBUsoW6ih9CU5PtNars0mBT7uHzpezSEDN/vn21DLYvfQnIMFqrrc2Jw5C9+P0sZjVDaGM8b2Eooo0SjjOEEwzlF9xIcE0JzJhhW0mhYDPLl0XYF5lAmEJ6wqx3K6wLeYUfld6Ob0ejbbWpwaPBptxl5MikTbGd//bt3mmO/MYb6e8bOjRnzyn6aqbGRm3tFNNOGaW000URbZSxhFUEN54GCxcOei2haYu4cV4br3VN5CTJjQAEQyFdjOEIT078Kr62TYNek7KHBptylxUrgAztXqMjnZwWPaSadrT2hDtmEA5IfT2LK55nCG0M4wRDOU47Q2inGIhwiHHU8BOub/gkoYXfGbQy/IUPccnW/+ZVphChABN7hbtHa1awFdHF8pp9+HY/O2i1KPtpsCl38futVbQT9FrSJtfPtfXVA9NtXUZOkW9HI7MKnqeDEk4wDKLdF9spIwK0UcpTzOVTDdcRXPhMVp87FGxmduF6Hop8gY5Yh/74W8/HiXJauH/6Wvz1lVmtQTlPg025T4qu/71E23DlpFAINm1KP1pLcSg2Fy3+r1Lex1uU0k4xXXRRGF3fzNrlhCliD5O4ueFjBOf8KCvPGZjxNFfMG826yBVxa6l1v9I9I7YCDBfxCs9XPYR/wzVZeW7lLhpsyn3SnF+KhcHevbk7aluwwOr2ns7dd9tXyyDyBRaxomoLxXTSRll0JeqeWYgW4RijuHnNNQQn3Paefqf+skdZunEORxiLtVuLf56ejxGFdHEeO3my5nf4mh4Y8PMpd9NgU+7j88G55wIZRm1+v23lZE0gAC+/nH60Nnx4bv5caVQ3fYZVNZso52i0tXDib9P6+hijmH8wwPXT/kkosPaUniMUWMts+S0r2xfTs0Bn8nN0h9soWnigZge+eg+cq1VpabApd/r2t1NujoXCc8/l3szB++7LfP+999pShp2q669mVe3rjOJowj29z3d1UEpj5BouWTqD2fJbQpPnZXxc/6RGSuUkU5deyTo+QfpdWc9zjOEgq2t3apusPKDd/bW7v3tNmwZbrQa6KUc5w4fDu+/aXdXABIMwb1760ZrHVzEI+l/k5pUXcYxRWXzUDNc8xhgK6eQK2UDglxPxVVdk8fnzl3b3V2qgtmyB6dOTNsfCobU1N0ZtwSDc2cfilN/7nj21OKQ6cCk/mF7PSI6SsOJenMQZjH3d+mIYxrv8rerrrItcqaGWR3TEpiM29xs6FE6eTD3aKSmB9nYHijoFF18MO3Zg4iaN9NotDxkCJ07YXpYTggufobbBxz7Oil4EkNi3MVusvicPTl+rMx8HgY7YlHqvMo3aOjps6WIxYMEgvPpq+pmQIn2fe/OQ6vqradp2LnPP3sxwWkk/ehuI7hZZXZzDTn5eu1VDLU9psCn3q6sDMnymb2hwb6utW28FMnQZufJKT82E7A+fDxr3fIQX6jZxZcEfyE64GSDCv/MC2yo/yy5TQXXg0iw8rspFGmzK/Xw+a6IIaTr/A9xyi50V9U8wCIcPp99tT5gQC+185PPPZl14JrWVfwC6GHjAWaO0upHf5KVto/BtyeGWZCorNNhUbkgzFT62K2xpcd9EkkWLkjb1Gq1de23udfAfBIEtM2lsLAY6Sey6n1rv76ks2o4xJfiP3quvpwJ08ohOHsklF1wAO3cCaQ7tlZfDO+/YXVVqkydDc3P6Q5ClpbBpk+6IMxDpxLroOlEYY4rtLkfF0ckjSmXLr34V+2fKQ5ItLbau9ZVWMJgUakmWLNFQ64MxxRhTkOKmoaYy02BTuSPuXFuiWIhs3OhsuIVC1iKhCXqN1mpq3HfYVCkP0WBTuSXuXFvaWZIbNzo3S3LBAvjXv9IfggSor7exIKXyjwabyi1+f8oREbhglqTfn9TkOEmlrv2l1GDTYFO5JxCAWbOA5NGQY+fbAgFYuTJpc6/6KiqsNmFKqUGlwaZyU12d1bWDDOG2caM9Fz8Hg7F11NIegrzoItixY/BrUUppsKkc5fPBF78Y+zJtuH33u4N/vu2mm+DEicyHIL/xjcGtQSkVo8GmclcgkPZ8W0wkAvPnD96K26efDsePJ4Var6AdOxaqqwfn+ZVSSTTYVG4LBKCqCsgwauvogJkzsx9uY8fCgQOZQ62oCL7//ew+r1IqIw02lfuammIrAKQNt8OH4ZJLsnNYMhiEsjI4cqTvUHvwQR2tKWWzIqcLUCorNmyAESOgtRWh9yQOQzRwOjpg3jyYOBF27x7Y88yYYU1KIbmTYdL1aj//uYaaUg7QEZvyjgwXb/cKoT17rBmVp3o5wOTJ/Q+1ykoNNaUcosGmvMPvj51vg9Th1iuQNm6EUaP6ftxAwArC5ubY48RLCrUJE/R6NaUc5JlgE5HbRGSXiLSJyF9FxAXdcJXt4s63Qeq2W70C7tgxK7RGjLAWBe2eYDJnjrVdBJYuTbuQStLjjxkD+/e/959DKTVgngg2EbkBeBh4AJgGvAj8VkQmOlqYcsaGDUkjtz4DrrUVVq+GqVOtMFuzJvl7EiQ9ZkWFNUlFKeUoTwQbUAv82BjzP8aYV40xXwDeAhY7XJdySlNTUl/GvgIu1S2VlI9TU6OdRZRyiZwPNhEpAf4NWJtw11rgUvsrUq6xZUuvkVu3tKsC9CFloIlYF4lrx36lXCPnV9AWkfcB+4DLjDF/jNt+NzDfGDM54fs/B3wu+uUUYLtdtbrcOOCQ00UMlovhwiEwrD/f+zZwWj++rwPaX/b++8fT74tTpK9Fj8nGmBFOF5GOl65jS3VePym1jTGPA48DiMhmNy9vbid9LXqIyObd+loA+r6Ip69FDxHZ7HQNmeT8oUisT1Bh4PSE7eOBA/aXo5RSykk5H2zGmA7gr8CshLtmYc2OVEoplUe8cihyJdAgIpuAF4BbgfcBfXWffXywC8sh+lr00Neih74WPfS16OHq1yLnJ490E5HbgLuAM7BO6H8xfjKJUkqp/OCZYFNKKaXAA+fYlFJKqXgabEoppTwlL4NNGyaDiHxFRP4iIsdE5G0ReVZEpjhdlxuIyHIRMSLyqNO1OEFEzhCRJ6LvizYR+buIXOZ0XXYTkUIRuS9uX7FLRL4pIl6ZdJeWiHxMRH4tIvuifws3JdwvInKviLwpIidF5Pci8gGHyk2Sd8GmDZNjLgcew2o7NhPoAtaLyBgni3KaiHwYuAUIOV2LE0SkHGtmsQBVwEXAF4CDTtblkC8DS4A7gAuBO6Nff8XJomwyHGsS3p3AyRT33wX4sd4bl2C9P9aJiCu6keTd5BER+TMQMsbcErftn0DQGJMPb9iURGQ4cBS4xhjzrNP1OEFERgF/wwq2u4Htxpjbna3KXiLyAFZ7uo86XYvTRKQJOGyMWRS37QlgrDFmjnOV2UtEWoHbjTE/jn4twJvAo8aY+6PbhmCF21JjzGqnau2WVyM2bZic0Qis98M7ThfioMexPuA853QhDroG+LOI/EJEDorIVhG5PbozyzcbgStE5EIAEbkY6+jGbxytynnnYnV6iu1HjTEngT/ikv2o548VJxgHFJLcausAcKX95bjKw8BW4E9OF+IEEbkFOB+ocboWh70fuA14CHgQqAQeid6Xb+ccv4X1ge/vIhLG2l/eb4x5zNmyHNfdvjDVfvRMm2tJKd+CrVu/GibnCxFZCUwHphtjwk7XYzcRmYx1znVGtEVbPisANscdlt8iIhdgnVvKt2C7AVgI/CfwClbIPywiu4wxP3C0Mndw7X40rw5Fog2Tk4jIQ8CngZnGmNedrschH8EazW8XkS4R6QIuA26Lfl3qbHm2egv4e8K2V4F8m1wF8B2gzhjzpDHmZWNMA1b7vrw9Fx+1P/pf1+5H8yrYtGFybyLyMNan0ZnGmH84XY+DngY+iPWJvPu2GXgy+u98GsW9AExO2FYB7HagFqcNxfogHC9Mnu03U9iFFW6x/aiIlAEzcMl+NB8PRQ60YbKniMgqrPNJ1wDviEj3p69WY0yrc5XZzxjTArTEbxOR48ARY4zXFxJN9BDwooh8FfgF1iUxdwDLHa3KGc8Cy0RkF9ahyGlALeD55dKjs6TPj35ZAEwUkUqsv4k9IvJd4Ksi8g+gGfga0Ar8zJGCE+TddH/QhskAIpLuF/91Y8y9dtbiRiLye/Jwuj+AiFRhnXOcDOzBOrf2iMmznUX0mqz7gGuxDrO9hTWK/4Yxps3J2gabiFwOPJ/irieMMTdFZ8neA3weGA38GVjilg+CeRlsSimlvCvfjxUrpZTyGA02pZRSnqLBppRSylM02JRSSnmKBptSSilP0WBTSinlKRpsSimlPEWDTSmllKdosCmllPIUDTallFKeosGmlFLKUzTYlHIJEXlDRIyIXC4iZ4rIYyLyuoi0i8hWp+tTKlfk47I1SrldBdCItfjpCaDT2XKUyi06YlPKfQJYS6R81BgzzBgzHKh2uCalcoaO2JRyny5gljHmQPcGY8xOB+tRKqfoiE0p96mPDzWl1KnRYFPKff7kdAFK5TINNqXc522nC1Aql2mwKeU+YacLUCqXabAppZTyFA02pZRSnqLBppRSylM02JRSSnmKBptSSilP0WBTSinlKWKMcboGpZRSKmt0xKaUUspTNNiUUkp5igabUkopT9FgU0op5SkabEoppTxFg00ppZSnaLAppZTyFA02pZRSnqLBppRSylM02JRSSnnK/wNS91aKveJ5oAAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"#plot\n",
"keys = sorted(list(replica_trajs.keys()), reverse=False)\n",
"positions = np.linspace(0,10)\n",
"fig, ax = plt.subplots(ncols=1)\n",
"\n",
"\n",
"# plot potential energy surface\n",
"ax.plot(positions, ensemble.replicas[0].potential.ene(positions), c=\"grey\", alpha=0.7, zorder=-60)\n",
"\n",
"\n",
"colormap = { 0:'red', 1:'blue'}\n",
"\n",
"for traj in keys:\n",
" T = round(ensemble.replicas[traj].temperature,2)\n",
" #min_e = np.min(replica_trajs[traj].total_potential_energy[eqil:])\n",
" ax.scatter(replica_trajs[traj].position, replica_trajs[traj].total_potential_energy, label=\"replica T=\"+str(T), c=colormap[traj], alpha=0.5)\n",
"\n",
" ax.set_ylim([0,20])\n",
" ax.set_xlim([0,10])\n",
" ax.set_xlabel(\"r\")\n",
" ax.set_ylabel(\"$V/[kT]$\")\n",
" ax.legend()\n",
" ax.set_title(\"Temperature replica exchange\")\n",
" fig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Time-dependent bias"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Time-dependent biasing methods update the bias during the simulation time. A frequently used time-dependent method is metadynamics/local elevation. There, a gaussian potential is added to the positions that were already visited during the simulation. Therefore, visiting these positions again, is energetically less favorable then in the previous visit (energetic penal).\n",
"\n",
"Note, that in case of a Gaussian bias potential, the mean of the Gaussian is set to the current position of the particle and its width should be chosen small enough to avoid a big penal for neighboring states. Step by step the energetic minima are \"flooded\" and the particle can cross barriers more easily. In most applications the bias is only added every $n^{th}$ step to iterate between free diffusion and biasing steps."
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"#### Metadynamics / Local elevation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We first define the original four-well potential. To apply metadynamics/local elevation we use the `metadynamicsPotential` function. In the initialization we have to define the height ($amplitude$) and width ($sigma$) of the gaussian bias function added. This bias potential is added every $n\\_trigger$ steps to the current position.\n",
"\n",
"Adding more and more potentials every step leads to an energy function that demands more and more computation time every step. To avoid slowdown of the simulation the metadynamic bias is usually stored and calculated grid based. This allow a much faster simulation but comes at the cost of additional input parameter and small numerical errors in the bias force calculations. To initialize the grid, the user has to define the minimum x-position ($bias\\_grid\\_min$) and the maximum x-position ($bias\\_grid\\_max$) the grid should cover as well as the number of grid bins. Note, that no bias will be added to values below $bias\\_grid\\_min$ or above $bias\\_grid\\_max$.\n",
"\n",
"In our example system we want to sample all four energy minima. Therefore, it is sufficient to set the grid between 0 and 10, which covers all four minima.\n",
"\n",
"For the metadynamics/local elevation simulation we will reduce the simulation length by the factor of 10. This makes it easier to visually distinguish, and thus understand, the effect of metadynamics/local elevation."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\r",
"Simulation: Simulation: 0%| | 0/200 [00:00\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
position
\n",
"
temperature
\n",
"
total_system_energy
\n",
"
total_potential_energy
\n",
"
total_kinetic_energy
\n",
"
dhdpos
\n",
"
velocity
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
4.000000
\n",
"
1
\n",
"
0.003797
\n",
"
-0.001344
\n",
"
0.005141
\n",
"
0.0034275749807509046
\n",
"
0.101399
\n",
"
\n",
"
\n",
"
1
\n",
"
3.999182
\n",
"
1
\n",
"
0.298541
\n",
"
0.298541
\n",
"
NaN
\n",
"
-0.00342757
\n",
"
NaN
\n",
"
\n",
"
\n",
"
2
\n",
"
3.946611
\n",
"
1
\n",
"
0.303253
\n",
"
0.303253
\n",
"
NaN
\n",
"
0.0426495
\n",
"
NaN
\n",
"
\n",
"
\n",
"
3
\n",
"
3.974171
\n",
"
1
\n",
"
0.299241
\n",
"
0.299241
\n",
"
NaN
\n",
"
0.127614
\n",
"
NaN
\n",
"
\n",
"
\n",
"
4
\n",
"
4.072425
\n",
"
1
\n",
"
0.300863
\n",
"
0.300863
\n",
"
NaN
\n",
"
0.0388076
\n",
"
NaN
\n",
"
\n",
" \n",
"
\n",
""
],
"text/plain": [
" position temperature total_system_energy total_potential_energy \\\n",
"0 4.000000 1 0.003797 -0.001344 \n",
"1 3.999182 1 0.298541 0.298541 \n",
"2 3.946611 1 0.303253 0.303253 \n",
"3 3.974171 1 0.299241 0.299241 \n",
"4 4.072425 1 0.300863 0.300863 \n",
"\n",
" total_kinetic_energy dhdpos velocity \n",
"0 0.005141 0.0034275749807509046 0.101399 \n",
"1 NaN -0.00342757 NaN \n",
"2 NaN 0.0426495 NaN \n",
"3 NaN 0.127614 NaN \n",
"4 NaN 0.0388076 NaN "
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sim_steps = 200 # reduced simulation length\n",
"\n",
"#Simulation Setup\n",
"origpot = fourWellPotential(Vmax=4, a=1.5, b=4.0, c=7.0, d=9.0, ah=2., bh=0., ch=0.5, dh=1.)\n",
"\n",
"#Performe metadynamics\n",
"totpot = metadynamicsPotential(origpot, amplitude=.3, sigma=.21, n_trigger=10,\n",
" bias_grid_min=0, bias_grid_max=10, numbins=1000)\n",
"\n",
"integrator = langevinIntegrator(dt=0.1, gamma=15)\n",
"\n",
"system4=system(potential=totpot, sampler=integrator, start_position=4, temperature=1)\n",
"\n",
"#simulate\n",
"cur_state = system4.simulate(sim_steps, withdraw_traj=True, init_system=False)\n",
"\n",
"print(\"Trajectory length: \",len(system4.trajectory))\n",
"print()\n",
"print(\"last_state: \", cur_state)\n",
"print(len(system4.trajectory))\n",
"system4.trajectory.head()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Metadynamic systems are visualized with the default `simulation_analysis_plot` function. `simulation_analysis_plot` will display the resulting potential after the last simulation step."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(, None)"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABHcAAAFACAYAAAA26JCpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd5xcdb3/8ddnZvtuNnXTaIFQQihSQpVIBKTZpahYiCLYfteKDSzovSp2773oVfFCEC4oRQQUERACIk2KQqRLAiEhm02yKdt3Zz6/P75nds/uzsz2nS3v5+NxHlNO+8yZ3XNmPvP9fr7m7oiIiIiIiIiIyPiUKHQAIiIiIiIiIiIyeEruiIiIiIiIiIiMY0ruiIiIiIiIiIiMY0ruiIiIiIiIiIiMY0ruiIiIiIiIiIiMY0ruiIiIiIiIiIiMY0ruiIiIyKgws5+Z2VfyzL/AzH45mjGNRToOIiIiMlDm7oWOQURERCYZM1sGXOXuO4/CvpYDH3L3Y0Z6XyIiIiKFoJY7IiIiIiIiIiLjmJI7IiIi0ouZrTGzL5nZU2ZWb2aXm1lZbP65ZvaCmW0xs5vNbH70vJnZj8xso5ltM7MnzGz/aN4KM/sPM6sE/gjMN7OGaJpvZheZ2VWxfbzFzP5pZlvNbKWZ7dsjvvOj7W8zs9/E4xvA6/yAmT1tZjvM7EUz+3Bs3jIze8XMPhu9nlfN7AOx+TPN7BYz225mf4te232x+YvM7I7oGD1rZmdGzx9pZhvMLBlb9u1m9kR0v/M4mNkCM3MzO9vMXjazTWZ24UBfp4iIiExsSu6IiIhILu8BTgIWAnsDXwYws+OAbwNnAvOAl4BfR+ucCLwuWn4a8E5gc3yj7t4InAKsd/eqaFofX8bM9gauAT4F1AC3AreYWUlssTOBk4HdgQOB5YN4jRuBNwHVwAeAH5nZIbH5c4GpwE7AOcBPzGx6NO8nQGO0zNnRlIm/ErgDuBqYDbwb+KmZ7efuD0brHRfbz1nRsrkcA+wDHA98NZ7oEhEREVFyR0RERHK5xN3XuvsW4JuEBAWEpM9l7v6Yu7cCXwKOMrMFQDswBVhEqO33tLu/Ooh9vxP4g7vf4e7twPeBcuDo2DL/5e7ro/huAQ4a6E7c/Q/u/i8P7gFuB5bGFmkHvuHu7e5+K9AA7BO1ujkN+Jq7N7n7U8AVsfXeBKxx98vdvcPdHwNuAE6P5l9DdDzNbApwavRcLl9392Z3/wfwD+A1A32tIiIiMnEpuSMiIiK5rI3dfwmYH92fHz0GwN0bCK1zdnL3u4BLCK1aas3sF2ZWPYh999xHOopnp9gyG2L3m4Cqge7EzE4xswejrlNbCUmWWbFFNrt7R5b91ABFdD9G8fu7AUdEXcq2Rtt+D6GVD4RWOu8ws1LgHcBj7v4SuQ35tYqIiMjEpeSOiIiI5LJL7P6uQKbr1HpC8gLo7II0E1gH4O7/5e6HAvsRumd9Lsu2+xqus+c+LIpn3cBeQm5RYuUGQqugOe4+jdD9y/qxeh3QAcRH+4ofr7XAPe4+LTZVuftHAaKWPi8Ruqf11SVLREREJC8ld0RERCSXj5vZzmY2A7gA+E30/NXAB8zsoChB8i3gIXdfY2aHmdkRZlZMqCvTAqSybLsWmGlmU3Ps+1rgjWZ2fLStzwKtwP2DfC1mZmXxCSgBSokSNWZ2CqFmUJ/cPQX8FrjIzCrMbBHw/tgivwf2NrP3mVlxNB3Wo1bO1cAnCDWKrhvk6xIRERFRckdERERyuppQg+bFaPoPAHf/M/AVQquXVwkFl98VrVMNXArUE1qmbCa0jOnG3Z8h1Jh5Meq2NL/H/GeB9wL/DWwC3gy82d3bBvlajgaas0yfICSS6gktaG4ewDb/H6HY8gbgyuj1tEbx7yAkit5FaIW0AfgOIZmUcQ2wDLjL3TcN7mWJiIiIhEKHhY5BRERExhgzWwN8yN3vLHQs44WZfQeY6+5n97mwiIiIyDBSyx0RERGRQTCzRWZ2oAWHE4ZKv7HQcYmIiMjkU1ToAERERETGqSmErlXzgY3AD4CbChqRiIiITErqliUiIiIiIiIiMo6pW5aIiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5I6IiIiIiIiIyDim5M4YYmYXmZlnmdrNbKOZ3W1mnzKzqkLHOl7EjuHyEdj2ymjbK4awjUoz+5iZ3W5mG8ys1czqzewZM7vTzL5hZsebWdkwhj6umNmK6DivzDJvWew9XjDqwYmMIXmuIbmmiwodc9xk+38ejuuTme1tZj80s8fNbKuZtZnZ+ujxb6Lry37DGLaIyJgy2PPgUM/BZrY8s408y8w3s5+Z2b/MrCW2z4MGs0+RvhQVOgDplyKgBlgWTf/PzE5x9+eHcyexk9MH3H3FcG57OJnZMuDu6OHu7r6mcNEMnpkdCNwI7NFjVgkwDdgHOB74CvABYMVoxiciImOXmX0a+A5Q3GPWvGg6CDgzs/gohiYiMirG8nnQzKqBB4Bd8yyzEjgWuMLdl49OZDKRKbkzdu0HvBzdLwb2Aj4PnAYsBG4yswPdvaNA8ckQmNkM4A5gNtAGXAFcBawBGoH5wOHAW4BTChOliIxj8WtILm2jEYgMPzN7L/DD6OFLwI+Ae4BXCD8Q7A28DjgD2L8QMYqIjKRxcB48i5DYSQPnALcBDdG8pgLEI5OAkjtjV5O7N8QePwycbmZ3ACcA+wLvAK4tRHDjhbuP1V8rP0VI7ACc6e439Zi/GXgS+F8zmw9UjmZw44W7r0S/SItk0/MaIhPLN6Pb1cCh7l7fY/56YCXwjai1q4jIRFPQ82DUy2FFnkUOiG6fGMs9ImRiUc2d8ec/YvePL1gUMlQnRLfPZ0nsdOPu64e7C56IiIxPZrY3Xc38f5nlC003URJcRGTCGCfnwYrodlsB9i2TlJI7488/Y/d3zraAmZ1mZrfECvRuNLM/mdl7zaxXK4dMYeDYU5dnKby5IMt6883sO2b2DzPbZmbNZvaCmf2PmfWsIxNfb028mKeZnWlm90SFhJuiAmifMLNklnWdrno7AKt7xLmm5/K5iqWZWZmZvdHMfmFmq8ysIVaE7WYze1uu1zAMZkW32we7ATMrMrPXm9mPzeyx6D3IFN++PSr01usYxtbvVhDazE6J/k7qzGyHmf3NzN7XY51FZnZZ9B62mNlqM7vYzLK2LLIeBVLNbLaZ/ShWWK42KnY3qMJyPbefZX7n+2/Bh83sYTPbHr3G+3u+xhz7WWJm10fxNpvZ82b2fTObFb2uzH6WDeZ1iBSame0anYPdzHImnM3sC7G/91N7zOs8t0f/b+eY2QMWClzuMLMHzexDZr2vQwOMdUDXuGidbv+nForZf9XC9Wu75ShwGe3rd9F1odXMNpvZXdE5Je9nKDPbw8wuN7N10fnuJTO71MwWDuX103X9gKFdQ+Lnx2Iz+6yF6++O6D2728zekWd9M7MjzOxbZvZQ9PfTHh2jey1cx8v7EUeFmX062l9tdJzXRtv4vJntkmfdQb8/IjKuDct5MGMwnxEtR0HlzLUQWB49dax1/65yUTQ5od4OwNk9lun13cXMjjWzX1vXZ/AGC5/D7zWzr5nZoqEeB5kA3F3TGJmAiwCPpgU5lpkeW+bWHvPKCAV6Pc90JzClx3or+1inVzzA6YTaMLmWbwZOy/Ea1kTLXAT8T55t/F+WdfuKc02O5Zdn2daP+rG9KwHL8Toyx23FIN7rB6N1m4B5g/x7+WQ/4r8TKO8rfuCrebbx9Wj5kwl9hbMtcy9QlGUfy2LLHE/oB51t/XbgXTniXBEts7KP7ff6n4nN+xBwc57X+M08x/lsoCPHei8Dr489Xlbo84imyTvRj2tIH+u/O7b+h7PMPwhojeb/JMv8NXSd26/J8//2236cL3rFzyCvcdG6C2LLnA48nWXdg2LLTwVu72Nff862r2j948h9jdwKHBF7vHyA79Oi2Lo3DuHvJbON8wh1KnK9zv/Osf5b+zg+DjwO1OSJYQm5rwuZaUWW9Yb0/mjSpGl8T8NxHoytP6jPiITkjQPe4/k1fZybLqL79TrXtDy2zQv6sfyPC/2+aCr8VPAANMXejP4ld46OLfPLHvMui827OvrQNJPwgfzS2Lzf9VivHKiKzf9w9Dg+WWz5E4BUtOy9wJsJFelnEj7Q3h3NawUOzvIaMie9FwlFxv4TOBCYARxCKDScieWNPdatIhQYzsxf3CPOih7L5/zwDPw78GvgndF+5xEKGR9NSDq1R+t+Isd7sZIcHzz78V5/Mxbb44TESa8vO31s4yOEi9HZhOLLOwNzovf9O3R9sfhhH/Gvjm4vj47DjGgb90bPp4CTCF9I/kpI0swC9gR+FnsdH82yj2Wx+aujbXwM2CWK9V3A2mh+G7B/lm2siOav7GP7vf5nYvNejP4ev0YYhWwGsBR4NPYaD8iy/iF0JXZeJHwpnB0d648AW4B/xfazrNDnEU2Td2KIyZ1oG/8Xrd8I7BN7vpzQctQJiZFeSWO6zu3xc8pBhGvDEcDvYvF9O8v6ff0/D+oaF627IDb/FUKi+vPReWwWIUk7O1o2Qdd1rJGQ/N6f8OPKQuB8us6v12bZ1/zoXOeE+mnnATsRrjHvA16NHaOs16c+3iOje0Lkp8Beg3iv4+fmFPB9Qj2/WdF7cV9smWzJvjcSkmkfBY4CdovWfQ1wIbApWve3Ofa/F6G7ghN+ef9y7DjvShhQ4HLgZz3WG9L7o0mTpvE/Dcd5MLbuYD8jLs9so8fzFYTvJFdF8/9C9+8qJdFUFc3zaNme372Kou0tout71+3AGwifo2sI178zCDVYe11XNU2+qeABaIq9Gf1L7sQzy++NPX9o7Plf5Fj3O7FlTs4yv88PmoQi3C9Fy90GJLIskyT8YubAH7PMXxPb1/lZ5pfT9YU/2wfnZX0dp4G8pjzrnhetu5YsrXcYWnKnJvYaM9N2wgfliwkJsyH94gicSNeH3+o88Tvwoyzzp9H1wbsduB8oybJc5sL0YB/vVQdwVJZldo/t5w9Z5q9g6MkdB07PMn8OXV8CvptlfibRuAnYKcv8I+jeqmfZUN4zTZqGMtH9GtIz8d1ryrGNqbFz9CNAcfT8JXQlYQ/Jse6a2P5/mmW+ATfFzinze8zP+f/MEK9xdE/upIHj8xzHzLm/Azg2xzKvj7bjwBE95mWS3u2EIp89192b7q0glw/ivX5Xj/ObE5I0vwE+CxxGjlansW3E1/18lvmlhMEcnJCkKhtgjIujY5Amy5cuwvXOgR3AgXm2U9Tj8ZDeH02aNE2MaajnwR7rDeYz4vLM+jm2v4Icn19jy6ykj+8SwCeiZTYQXZM1aco1qT/yOGChtsq+ZnYZ4Us/hATLdbHFPhjdNgGfy7GprwIbo/vnDDKctxJ+UXPgg+6e7rmAu6cI2W+AkywM+53NS3QNYRhfvxm4Pnp42CDjHA5XRbc7Ez6MDxt3rwOOIXy4zZhCaBXzBUISrzaqzzB3kPu4Hagj/IJwVJ5FM7989lx/K+EXAghJvQvcPdvQyb+Jbg82s3wj8F3r7g9k2c9q4MfRw5PNbE6ebQzWfe5+fc8n3b2WrtfY7W/NwihlmaLlP3D3dVnWf4jQ+ktkrPkn4UtzzsnMpvVcyd23Ae8nfDE+FPi6mZ0MfDxa5Kvu/lgf+24GvpRl2w58mnD9KALeM4DXM5zXuD+4+5/zzP+36PaX7n5PtgXc/W7grujhWZnno3Ng5vGV7v5olnWfA36SZ/99cvdfA2cSRoPJWBA9931CUuZFM/toP2rPrAV+kGUfrYTrEYRfst/cc5k+YnyK0DLVCK16O5nZvnSdX7/h7k/k2U5Hj6cG/f6IyMQxjOfBAX9GHGWZ+pmb3L29gHHIOKDkztjVWSiY8MvXU8AHonkbgLdFH7wyjolu74w+nPcSLX9L9HDpIOPKfEB7GthuZlXZJuDZaDkjdG3J5s5syaFIZnSokfii38nM5kRFyP4aFWJsjx33xtiiw5rcAXD3l9z9DYQvUN8BHiI0C80oJ/QDftzMFueIvzoqgpkpRNkWL8ZGaCHUV/wPuPuOHPNejG5bCE30s/lXdFtC+AKQy+/yzLsxuk2QPxE1WH/KMy/X39qRdA2zfgu53TzYoETGIne/F/hu9PALhK5aELpqfjfrSt2tzHMdehHIfJF/7QDCGs5r3B9zzTCz2YQuPgB/yXWNi65zT0bLHRrbxAGERD3075w3aO5+HbAHoWvxlXSdrzMWELoqXG95iusDv49+lMlmJZAZhabX+2VmJWZ2npndFitqHL8GZb4Q9bwGxZM9V+aJref+hvr+iMgEMkznwcF8RhxNf49u9zOzb5rZ9ALGImNcvl/ZZWxpJSR4biYUN9zcY/5u0e1TfWwnM9rWHDMrydESI599otvFhF9/+6Mmx/Ov5lmnKbrtc6SNwTKzYwkfsPtzkpw6UnFEv4I/FsVUTOg/+yZCHYMaYC7wazN7TfTLN9GyiwkXpKyjpvWQL/4NeeY1R7ebsvx62nMZyP9+PZNn3tOx+7vlXGrwBvO3Fo/juTzrP5tnnkih7O7ua4aw/lcJXTszdbi2Ae/Lk5CPy/e/npn/Ggb2vz6c17jVedaPJyGuyrlUd/Fr3ILY/XzHoa9j1C9RQuvaaCL60L+U0FrlDELC/O2EFlPfH2gs7u5m9hyhC2q398vM5hG6ru7Xj1B7XoMyI4ZtdPd816Cehvr+iMgEMwznwYJ+H+mLu99tZrcQWk9eAHzezB4m/OByD3B3jx/8ZRJTy52xaz/Cr39TCIUry9z9EHe/KEtiB0INBQj9+POJJ2Sm5Fwqt8EkOUpzPJ/rl8IRF3VHuIGQ2NlAKMB4OKHg5VTCsamOrTIqiVB3b3f3v7n71wgJtMwvBgcQ+yU6avp/AyGxs4NQa+O1hKKd0+j621nbj/j78z70973KN8RxY64Z0UUps4+qXMsNQX/i7xl7Znj3jj6SoH39z4mMO1HT73/EnrrP3V/u5+o5/9cjmf+ZgfyvD+c1rjnH8zD0a1xl7H6+4zAi5w13r3f3m939XcA7CF3gIBSyz2Ww79eVhM8q7YRuXa8ndNueTtc16K/Rsj2vQZn3pr8/EmUM52cQEZmABnEeHMxnxNF2OqG780uE8+nRwBcJLVE3mNnXzaykgPHJGKGWO2NXk7sP5MNfA+FDT18fluPzB/qhCro+BN7k7m8bxPpjxemEUVZSwOvdvdcvl9nqUYwmd99kZhfQVVvpYEKWHkLh0UXR/dPc/Y5s2zCz6mzPF0hlrhlmVkpXn+KxkizJ/K0X9dHKbSSSUSIFZWZvoasrMMAbzeyd7v6bXOvE5Pxfj/Q3URM32tc4CKM9/j3nkn2vn+84jPh5w91vMrNbCaNa7W5mU3N0aRvw+2Vme9JVM+fj7n5pthWj7lHZZN6bgf7INNT3R0QmkQGcB8e06DPoxcDFZraIUMJgGaGl/wxCa9t9CfWGZBJTy52JY010u28fy2WaT28YRJcs6OrLevAg1h1LDoxun8yW2Insn+P50fTP2P2K2P1M/PV5Ejs7M4LdyQZhUZ558b/bl0Y6kH6Kx5GvZtGw12MSKaSoqPkvo4d3ArdG9/8nOq/0Jd//enz+QP7X10S3o3WNg8Fd59bE7uc7Dn0do+GS6xoSlzMWMzO6znHx9+vA2P1rc6xbTO7z4wvR7ewBDhow1PdHRCaf/pwHxw13f8bdL3f3swkt+DPn4DOiYvUyiSm5M3Fkit2ekKu1RtRc7809lo/L1FTJV3gxk0jY1cxeP+Aoh0e8Uny+WPPJNNPOt/5YGGEj/kUqPhrAeIk/Ll9Lr8y8NNBrRK0CeZCuprxvyrPcW0YhFpHRdBmhTskWwlCvHySMvDcduCL6wp/PMjPLmlg2s93pSgz8NdsyOQzHNa5PUdezTI2t5YPYxJN0tUrpzzlvpGWuIe3AphzLvClPodFj6apLF3+/4l2dcq37dnLXqbgrdv+9OZbpZRjeHxGZfPpzHiyUzHeaQX2fiUYY/nbsqdH64UDGKCV3Jo7LottKQrO9bL5OV7X3/80yP1PLZ16e/dxAVx2Xn/c1bLWZ7ZNv/iDFaw7lizWfTEHNRVHz8m7M7LXAuYPcdp/M7FtmlndYxair0pejh2m6fxjOxF8dFYbuue5ehKJrY8mZZtZrJKzoy96nooe3RUNPFpy7r6frmH8mKh7ajZktAd41qoGJjCAz+xhwavTwI+6+Lvqf/FD03HGEopT5lNP9w2Zm2wb8iFC7oIOuUbj6Yziucf31o+j2dWb2mXwLmtmU+LkhKjx/dfTwvWbWa6QmM9sb+H+DDc7MFprZf5hZvtEJMbODCPUmIIxglmsI3V2Az2ZZv5Su0dG20H3UwHhR6l5DpEejWuUcWS1qMXtn9PArZpazpWxUYy5u0O+PiEwMI3AeLJQ+v3uZ2V6Wfyj3hbH72eqyyiSi5M4EEY24lPnw+1Ezu8rMDjWzGWZ2oJn9nFB4C+Bmd78ty2Yei27PNrMjzKzSzIriH6yiZu7LCbVq9iIM0/0pM9vXzKaZ2dxo3U+b2YPA9SPwcl8Atkf3v2Bme5pZaRRrfzPfvyUkTIqBP5jZm6LYdzez84Hb6CpmPBJOBB42s7+Z2flmdrSZzY+O4Z5m9n7C0Oivi5b/ubuvja3/J7p+Hb7GzM4ys53MbGczO5fwq3UL4QP5WLEWuNXMPhrFOcfM3kmo9D+V8OvFFwoaYW9fJPyt1xCG3T3NzGZHx/o8wvvQ3yKzIqOpwvIMEx1NZfEVomT896KHv4qGmAXA3W8GMnVVvmVmB+TZ9xrCdehyM3tNdB06nDA64VujZb4fJVD7ZZiucf11KV2Jhx+Y2XVmdmJ0jZgenaPfbma/IJzXeg4R/g3C6GLFwO1mdm50fp9rZu8lDC9eN4T4yoELgXVmdo2Zvc/MFpvZTDObZWZLzOwbhBptZYRz2EV5trcG+LaZfc/MFkXbOTY6BpkfIS5w95bYOo/Q1U3rv6Lz+m7Ra3wnoQVmDfm73n2UcC2vBu4zswui1zEtuka80cx+Cfx3j/WG+v6IyPg33OfBQsl89zrGzN4RncOKoinTSvZC4PkomXV8dH6cbmZ7m9nH6bo2v8zYaf0uheLumsbIRDjpeDQtGMT6ZYQPz55nuhOYkmP9U/Kst6DHsqcC9X3sy4FHs+xnTTTvojyvZXlmGznmX5xjf2t6LJd5fnmWbXw5T9zrCaNV5Vt/ZTRvxSDeq7v6cewy0wqgOMs23k+4WGVbZxuhOX3OY92f+GN/k2vyLLMsz99JfN7xwLoc8bYD78qx/RXRMisHsu++3v/+vkZCl5Rcx3lt9Loyj48ZqfODJk19TXS/hvRn+l1s3WLCF3YntMqozrL9SkKXGAeeAEp7zO883wC/zrPfG4GiLNvv6/950Nc4whDlmWWW9eNYVhF+nOjPcXxLlvWPIxT/zXV+PrI/56ccse0FtPYztq3A6Tm2k1nmw4QvQLm28d851j+B8CNCtnVagXfSx3WGMErlq328hl7rDvX90aRJ0/iehuM82J9zMHk+I9L3d5UV5Pj8GltmHrm/Ty3vsZ18Ux1weKHfF02Fn9RyZwJx9xZ3fzthJKg/ALWEL811hFo57wPe4O5ZRxBx9z8SaovcTuiTmnNoQHe/FdidMCzffYRmgCnCaBpPA78i1BQYqV/MLiB05fkboQWLD3QD7v4fwBmE+BuBJuBZwpCuB7v7U8MWbW8nAMcQuhH8ifDrZguhq0I9IZP/E+AId1/uWZqRuvuvou3cTvj1s5XwpeznwCHufs8Ixj8Y/wIOJfwKu5oQ70bCaGCHufuvCxhbTu5+GeGL2I2E/6VWQlHPHxNez8bY4oMZnUdkLPg64e85DbzP3bf3XMDdGwn1UTqAA8jS9Srm3cBHgIcJ56fG6P65wDs8dF8akKFe4wa4rwZ3P52QpLmS8D/fFO2vltDi8EvA3h5aNfVc/y5CbaEVhB8L2gjJ4MuBJe7+4BBie57QKuZdwM8Ix3UT4X1pifZ3J/B5YC9376sFbSshSf1FQtKukfCerSR8Ifq3HHHcSRix5UZCK9HMa7wKONL7MbKauz9MKLr8RcIvzvWx7dwLfI7wq3XP9Yb0/ojI+DYC58GCcPdXCd+VrgZeoXtd0YwvEH7QvRL4B+FzZwchafUQ8DVgUXQ+lUnO3Af8nVhEpF/MbBlwd/Rwd3dfU7hoRoaZvY3w5Qagxt3HWrE+kVFhZmuA3YCvu/tFhY1G+mJmmQ+AH3D3FYWMRURERIZOLXdERIYmU0z0JSV2RERERESkEJTcERHJw/KMxGBmRxK6gkDoXiYiIiIiIjLqeg4vKSIi3V1rZhsIBWL/TqhHsROhptQXCIVot9I1PK+IiIiIiMioUnJHRCS/IuA90ZTNDuA0H8CwziIiIiIiIsNpUhdUnjVrli9YsKDQYYhMWDt27OC5554DYP/996e0tLTAEQ1cQ0MD9fX1NDQ00N7eTkdHB2ZGaWkp1dXVzJkzh+Li4kKHmdWjjz66yd1rCh3HWKDz/ch78sknaWtrY968ecyfP7/Q4UgfHn30UQB22203Zs2aVeBoZKh0vu+i872ITGT5zveTuuXOggULeOSRRwodhojIiDCzlwodw1ih872ITGQ633fR+V5EJrJ853sVVBYRERERERERGceU3BERERERERERGceU3BERERERkRFhZmvMzLNMf8ix/IIcy5882rGLiIwnSu6IiLGgi2YAACAASURBVMiYYGavM7ObzWxd9EF+eY/5ZmYXmdl6M2s2s5Vmtl+BwhURkf45DJgXmw4BHLi2j/VO7rHeXSMYo4jIuKfkjoiIjBVVwCrgk0BzlvmfBz4L/Bvhy8JG4A4zmzJqEYqIyIC4e527b8hMwKnAduC6PlbdHF/P3dtGPloRkfFLyR0RERkT3P1Wd7/A3a8H0vF5ZmbAp4CL3f0Gd18FnA1MAc4a/WhFRGSgonP5OcBV7t7Ux+K/NbONZvZXMzt9FMITERnXlNwZoPPOO4/ly5cXOgwRkclmd2AucHvmCXdvBu4Fji5UUCIiMiBvIJzPf5lnmQbgfOBMQiufPwO/MbP35lrBzM4zs0fM7JG6urrhjFdi0o2NbP3tjWz41rdINzbStnYtdZf8BHcvdGgiAhQVOoDxZtWqVTz99NOk02kSCeXGRERGydzotrbH87XATtlWMLPzgPMAdt1115GLTERE+utc4G/u/vdcC7j7JuAHsaceMbNZhK65V+VY5xfALwCWLFmiTMMIWX/Bhez4058AKD/gQFpWPcmWK37F1De/iZLdditwdCKi7MQAvfa1r2Xr1q088cQThQ5FRGQy6vmh3bI8FxZ0/4W7L3H3JTU1NSMfmYiI5GRms4G3ApcOYvWHgL2GNyIZiHRrKw333MO0M04nUV1N40MP0vjgQwC0rVlT2OBEBFByZ8De9KY3AXDzzTcXOBIRkUllQ3Q7t8fzs+ndmkdERMaeDwCtwK8Hse5BwKvDG44MRPPjj+MtLVQddxwVhx9Gw90raX32WQBaV68ucHQiAkruDNgxxxxDRUUF9957b6FDERGZTFYTEjxvyDxhZmXAUuD+QgUlIiJ9iwopfwj4tbvv6DHv22b259jjs83sLDPb18z2MbPzgY8D/z26UUtc41//CsXFVB5+OJVHHkVq8+bOeW2r1xQuMBHppJo7A5RMJtlvv/34+99zdhUWEZFBMLMqYM/oYQLY1cwOAra4+8tm9mPgQjN7BngO+DKh8ObVBQlYRET6axnh/P6eLPPmAQt7PPdlYDcgRTjff9Dds9bbkdHR8Ne/UnHQQSQqK6k88ggAElVVlOyxh7pliYwRarkzCEcffTSbN2/m6aefLnQoIiITyRLg8WgqB74e3f9GNP+7wA+BnwCPEL4QnNjzV2ARERlb3P1udzd3fzjLvOXuviD2+Ap3X+zule5eHdVOU2KngDrq62l96mkqXxsGpyxZuJCiuXOpPOpIShcupE3dskTGBCV3BuGUU04B4KabbipwJCIiE4e7r4w+/Peclkfz3d0vcvd57l7m7se6+6oChy0iIjKhZWrrlB94IABmxm4rLmfu175GyYIFdGzcSKqhsZAhighK7gzKcccdR2lpKStXrix0KCIiIiIiIiOm9fkXACjZc8/O50oWLKBo1ixKdl8AaMQskbFAyZ1BKC4uZt9991XdHRERERERmdBa//UCiepqimpqes0r3X13QMkdkbFAyZ1BOuqoo6itreXFF18sdCgiIiIiIiIjovWFFyjdc0/CoGfdFe+2GwBtL63pfM7b2ki3to5WeCISUXJnkE4++WRAdXdERERERGRicnfang/JnWwSJSUkp07tNjT6+gsuZN2nPzNaIYpIRMmdQTrxxBMpKSnhrrvuKnQoIiIiIiIiwy61eTOpbdso3bPnaPVdkjNn0rGlvvNx20sv0frMM6MRnojEKLkzSGVlZey999489thjhQ5FRERERERk2LW+EIop52q5A1A0Y0a3ljvpHTtor6vD0+kRj09Euii5MwRHHXUU69evZ/Xq1YUORUREREREZFhlGymrp+SMGXRs2dL5OLVjB7S3k6qvz7mOiAw/JXeG4C1veQsA119/fYEjERERERERGT6eTrP91lspqqnJOlJWRnLmDFKx5E56+3YAOmprRzxGEemi5M4QnHTSSZSVlXHHHXcUOhQREREREZFhs+3G39H8+OPUfOpTWUfKyiiaMZPU1q14Rwfp1la8vR2AdiV3REbVmE/umNnrzOxmM1tnZm5my3vMXxE9H58eHI3YiouLOeCAA3jkkUdIq0+piIiIiIhMAN7ezsYf/pDygw9m6tvflnfZ5Izp4E5q69bOVjsAHbUbRzrMCSe1bRu13/0ejQ8/XOhQZBwa88kdoApYBXwSaM6xzJ3AvNh06uiEBsceeyz19fU8/vjjo7VLERERERGREdP48MOkNm9m5ofOwRL5vzIWzZwJQMfmLaHeTqS9dsOIxjjRtL/6Ki++5a1suewytl57XaHDkXFozCd33P1Wd7/A3a8HcjWPaXX3DbFpS47lht1pp50GwA033DBauxQRERERERkxDX/+M1ZeTuVrX9vnsskZMwBI1W8hHUvuqOXOwDQ++BAdtbUU1dTQvnZtocORcWjMJ3f66Rgz22hmz5nZpWY2e7R2fPjhhzNt2jTuvvvu0dqliIiIiIjIiPB0mh1/vouqY15Loqysz+WLouROx+bNpLZHyZ1EQgWVByjd3ARA2QEH0PbKKwWORsajiZDcuQ14P3A88FngcOAuMyvNtrCZnWdmj5jZI3V1dUPeeSKR4JBDDuGJJ54glUoNeXsiIiIiIiKF0F5by5Zf/YqO2lqmnHBCv9ZJRt2yUpu3kN4Rau6U7LorHRuV3BkIbwrJndJ99ia1eTPpxsYCRyTjzbhP7rj7r939Znd/0t1vAU4B9gHemGP5X7j7EndfUpNnSL+BOP7442lqauLOO+8clu2JiIiIiIiMpqbHHuPFN72ZjRd/h+KddqJq2bJ+rZecOjW01Knf0tlyp3SvPWkfYres9tqNPLd0KS3PPTek7YwX6UxyZ889AWh7ZV0hw5FxaNwnd3py9/XAK8Beo7XPM844A4CbbrpptHYpIiIiIiIyLNpeWcfLHzyHopkzWXDD9Sz8020hadMPlkiQnD49tNxpCMmdkoULSW/fTro513g4/Yhp9WpSdZtoe+GFQW9jPEk3NWMVFZTsuhsA7a+o7o4MzIRL7pjZLGAn4NXR2udee+3FvHnz+Mtf/jJauxQRERERERkWrc8+g7e0MP/ib1O+335YUdGA1i+aMYOOLVHNnaIiShYsABhS3Z1U1MUrNUm6J6WbmkhUVFCyy84AtKmosgzQmE/umFmVmR1kZgcR4t01erxrNO/7ZnaUmS0ws2XALcBG4MbRjPPwww/nmWeeoaGhYTR3KyIiIiIiMiSprdsASM6aNaj1kzNnktq8hdSO7SSnTOksspzaunXQMaV3hO9Vk6X2TLqpiUR5OYmpU0lMmUL72lBUuf6666i9+DsFjk7GgzGf3AGWAI9HUznw9ej+N4AUcABwE/AccAXwLHCUu+/IurUR8sY3vpGOjg6uu+660dytiIiIiIjIkGSSMMlp0wa1ftGM6aS2bCG9o4HElCkkKiqArjoyg5EpzjxpkjvNzSQqKjAzinfZmbZX1uLubP75L9hy1VVD6uImk8OYT+64+0p3tyzTcndvdveT3H22u5e4+27R86Pehu3MM8+kqKhIdXdERERERGRcSW3bBkVFJCorB7V+Uc1s2mtrSW3dSnKYkjuZ4sxD2cZ4km5q7DxuJTvvQvvaV2h97jnaX3kFOjpo/scTBY5Qxroxn9wZL6ZOncrixYu5//77Cx2KiIiIiIhIv6W2bSM5dSpmNqj1S/fZB29poeWf/yRRPUzJncnWcifqlgVQvMvOtK9dS/1V/wdmYEbTI48UOEIZ65TcGUYnnHACdXV1PKJ/PBERERERGSdSW7f2e3SsbMoW7xu2U19PsmoK1pncGXxXonSm5U7j5Gi5403NJCrDcZv29rdjJSVsve46yg85hNJ99qHp0fAds+mxx3npAx9QNy3pRcmdYfTud78bgGuuuabAkYiIiIiIiPRPatu2QdfbAShduBArLgaIWu6E7l1Da7mTSe5MopY7UVKsdM892ek//xOKi5n65jdTceihNP/9H6Sbmnj1wgtpeuBB2tevL3DEMtYouTOMlixZQk1NDXfeeWehQxEREREREemXobbcseJiSvfaC4DklGoS5WXAEAsqb5983bIs6pYFUHXMa9n7gfuZ9s4zqVhyKN7czEvvfR9tq1eH5ZtbAPC2Ntad/zlann2uIHHL2KHkzjA7+uijeeqpp9i2bVuhQxERmVDMLGlm/25mq82sJbr9DzMrKnRsIiLSm5ldZGbeY9rQxzoHmNk9ZtZsZuvM7Ks22EIw0m9DbbkDULrvIgASU6qwZBIrLx+eljuTpaByc3Nni6eMZFUVZkbl0tdR/eY301FXR8keewDgzeG4tDzzDNt//3uaHnxg1GOWsUXJnWH21re+lY6ODq699tpChyIiMtF8Afg48AlgEfDJ6PGXChmUiIjk9SwwLzYdkGtBM6sG7gBqgcMI5/vPAZ8Z+TAnt6G23AEo23cxEFruACTKy0k3Db7VzXhruePurL/wQpoee3zg66ZSeEtLZ7esnpJVlez0ve+y11/uZf63vgnQWXOn5amnAUiNk+MkI0fJnWF2xhlnUFxczM0331zoUEREJpqjgVvc/RZ3X+PuNwM3A0cUOC4REcmtw903xKa6PMu+B6gAznb3Ve5+A/Ad4DNqvTNy0q2teHMzyWlDTO5ERZUT1VPCbUXFpKq5k25sZNsNv6Xh7rsGvm6UqEnEumXlYuVRseqoW1bL00937l8mNyV3hllVVRX7778/DzzwAOl0utDhiIhMJPcBrzezRQBmthg4Dri1oFGJiEg+e0Tdq1ab2a/NbI88yx4F/MXd48MA/QmYDywYySAng3RbG40PPdzr+VRUTmKo3bLKDzqIORdcwJTjjweGltzxdJp0QwMwfpIWmTg7Nm0e+LrRccqMlpVPZz2jzm5ZSu5IoOTOCDjllFPYvHkz9957b6FDERGZSL4DXAk8ZWbtwD+BK9z9p9kWNrPzzOwRM3ukri7fD8UiIjJCHgKWA6cA5wJzgfvNbGaO5ecSumTF1cbmZaXzff9sv/VWXj77bNprux/i1NatAEPulmWJBDPe/z6SU7pa7vggkzvpxkZwx8rKSDc14e5Dim00ZJJRHVsGntzJHKdc3bLiMq17vLkZT6VojQopT5Yh4yU3JXdGwAc/+EEArrzyygJHIiIyobwTeD9wFnBIdP9jZnZOtoXd/RfuvsTdl9TU1IximCIiAuDuf3T3a939CXe/E3gT4fvH2flW6/HYcjwf34/O9/2Q2hySDplkTkZ6mFru9BRa7jT3vWAWmXo7xXPnQjqNt7QMZ2gjItNyJjWYljtRtywbSLespmbaVq/uPDZquSNK7oyAhQsXsnDhQg2JLiIyvL4HfN/df+3uT7r7lcAPUUFlEZFxwd0bCK0u98qxyAZ6t9CZHd32bNEjA5Srhk3HMLXc6SlROfhuWZlYi+aGP4fxkLjobLmzeQjdsvrVcifqltXSTMvTz4TnqqvHxTGSkaXkzgg58cQTefnll1m1alWhQxERmSgqgFSP51LoWiYiMi6YWRlhtMNXcyzyALA0Wi7jDcB6YM3IRjfxpbdHyZ0oCdH5/Ii23BlkcifecofxkdzJjFaV2rx5wN3IBpLcsWQSKynBm5tp/dcLkExSvv/+4+IYycjSB+IRcvbZobXp5ZdfXuBIREQmjFuAL5rZG81sgZm9nTA87o0FjktERLIws++b2bFmtruZHQFcD1QCV0Tzv21mf46tcjXQBKwws/3N7B3AF4Ef+ngoujLG5Wq501lQeZhb7tgQkjvpTMudeeMnuZOJ0dvbO+Pv97pR97X+JHcgM8x8M+nt20lWVZGcNnVcHCMZWUrujJAjjjiCuXPncttttxU6FBGRieLfCF8Mfgo8DfwAuBS4sJBBiYhITjsD1wDPAr8FWoEj3f2laP48YGFmYXffRmipMx94BPgJ4Vz/w1GMecLK1LFJ9Wi5k9q6FSsuxvqZWOivRPlQWu6E5Ejx3HkAQxpSfaDaX32Vl85e3qt7Vcuzz7Hpf/6H5n/8I+t66Yau5MpAR8waSMsdCLV50i0tpHY0kJgyhURlpZI7ouTOSDruuON45plnWLduXaFDEREZ99x9h7t/yt13c/dyd9/D3S9w97FfZVFEZBJy93e5+3x3L3H3ndz9NHd/KjZ/ubsv6LHOk+7+Oncvc/d57v51tdoZHjlb7mzdRmLaVMws22qDlqiowFta8FTPHtV9S+/IdMuaEx6PYuKi6eGHaXroIZr//vduz2+57H+p+8//Ys0738X2P91OurGRLb+6svP1xWNMbd40oH2mm8K6A2q509xEeseOkNypqOzV3U4mHyV3RtBZZ51FOp1mxYoVhQ5FREREREQmsUzCJN7CBELLneHukgVdiYrMSFADkWm5U4iCym1rXwGgvccP9K2r11B+8MEUzZvHtptuov6aa6j91rc6k0Dx5MpAiyp7dIwS/RgtK7OcNzWTathBsqoqtNxpasLT6QHtVyYWJXdG0Mknn8zUqVO56aabCh2KiIiIiIhMYpmESa+WOzt2kKweweRO48C7VKV3bCdRUUFyyhSgq1jxaGhfuzbcxpI77k7biy9Stu8iqk98A4333Uf9b64Ny23YEGJubIRE+Ho90OROuqkJzLCysr4XBqwidMtKx7plhe0Mbuh5mRiU3BlByWSSY489lscff5y6urpChyMiIiIiIpNUZ7esnqNl7djRmUQZTonKkNzx5oEnd1Lbd3RLWnhT06i13ml7JbTcaYsld1KbNpFuaKBkwe5MOfFEvK2tMwnUsTF8z0s3NlA0ezaYkRpocqexiUR5OZbo39fzRFk56ebm6L2r6kruTOK6O+5Ouq2t0GEUlJI7I+z9738/HR0dXHrppYUORUREREREJiFvb8ejor3ZWu4kRiK5k2m5M4hiyKmtW0lOn96ZtNh2y+957sijaH/11WGNMZuuljvrO59rXb0agJI99ghds2pqsIoKrLSUjtraEHNjI8nqapLTpw+8oHJz84AKWifKy/HmJlINDSSqpii5A9RfeSX/OvEkfBIneJTcGWFve9vbmDZtGtdff32hQxERERERkUkoPkJWqrFHy53t20lWj7HkzpYtJKdPw4qKsNJSWlatwtvbaX7yyeEOs5t0SwsdGzeCWbduWW2r1wBQuvsCLJFgzpe+yNwLL6Bo7hw6NobkTrqxkURlJUUzZ9KxZeDdsvpbTBkgURENhd7QQKJKLXcgdKPr2LCBxoceLnQoBaPkzghLJpMcf/zxPPHEE6xfv77vFURERERERIZRZhh06J4AcPeo5U71sO8zUxx4MMmdjq31FE2fEbYTJS4AWp99bniCyyGT0CnddxHp7dtJRcetbfVqrKyMonlhaPbqU09l2mmnUTx7Du21G4FQqDpRVUVy1kxSgxgKfSDJHSsrp2PLFkin1S0r4u3tAOy4444CR1I4Su6MgnPOOYdUKsXPf/7zQociIiIiIiKTTKaYMslkt9GyvKkJUqkRabljQ2q5U09y+nQgltxJJGh97tlhiy+btqhLVuURRwLQHv0437r6RUoWLOhVE6dozpzQ0odYy50ZMwdeULl5gC13ysu7RtiKd8tqmrzJnUy9nR133dU5PP1ko+TOKDjppJOYOXMmv/3tbwsdioiIiIiITDKZYdCL58zp1rojU2R5ZGruZFqTDCy54+3toVDwjCi5U1GBVVRQ9brX0TLSLXeiYdArjzwiPI5a8rStXkPJ7gt6LV80ezYdtbWhmG9DA4nKCopmDTy501FXN6Dh6BMVXUOmJ6ZUdRavnswtd4ha7qQ2beocnn6yUXJnFCQSCU466SSeeuopVkfFuEREREREREZDpuVO0bx53UbLSkfJnWT1CHTLqhxcy52O+voQU9Ryp+Kww5h+xhmUHXgA7S+/PKgERt0lP2HNe9/bKznk7t0et7+yFquooOyAA8LjdetIt7bS/sorlO6+R6/tFs+Zjbe1kdq6lXRjI8mqKpIzZ4XRvfr5utONjbS9uJqyxYv7/XqsrCu5k5yigsoQWu5k/mZGujbTWKXkzig599xzSafT/PSnPy10KCIiIiIiMomkMi135s3L3nKnauwUVE7VbwWgaEaouTP3yxcy50tfpGyffQBofeGFAW3P02nqr76a5kceZfXpp7P+wgtp+Mt9bL3ht7yw7PW88ulPdyaUWlevpmSXXUhOn45VVNC+bh2tz78A6TSl0f7jiubMAaCjtjZWUDnE3bFlS7/ia3n6aUinKTtg/36/pkw9IwjvXTKT3OkxzP1k4m3tFM2ejZWX0/HqhkKHUxBK7oySZcuWscsuu3DdddcVOhQREREREZlE0lHLneJ5c/G2ts7hojMFg0ek5k5JSajx0zzQ5E7Ucmfa9G7PZ5IrLc/2XXen8f77Wfuxj/PyBz9I86OPktqyhTkXXMD0M05n+x9uZe255/LqhReSrK5mx51/5qV3nxVG43r875S/5jWYGSW77Ubr8y/Q+uwzAJTts3ev/RTNDsmdtqh3RqKyiuTMmeF1bNrUr9fbvGoVAOX7DyC5UxFvuVMV6huZkZrELXe8vQ0rKaF4zhzaN0zO5E5RoQOYTE477TR+/OMfc88993DssccWOhwREREREZnA3J2WJ58MSZxEgqKa2QCkGhspKinp7JY1EjV3zIxERcUgWu6EFi+ZLjYZxfPnk6iqouXJJ+HMMwHY9oc/kNq8hRnvf1/X+tu28fKHP0KyqopUfT1ta1+BZJKpb30LyalTqfnUp2h99lnSrW1UHn0U2353E69ecAFbb7iB9I4dVBy2BICKgw9i2+9uomTBblhFBcW77tor1qLZ4Xi2vvgiQNRyZxYwgJY7T66iaN48imbN6vcx6tZyZ8qUrmM9mZM7be1YSQlF8+bSMUmTO2q5M4o++clPkkgkuOSSSwodioiIiIiITHBNDz3EmjPfybbf/S7UZqmqArpqs3S13Bn+mjsQumbFR+fqj0wXqaIZ3ZM7lkhQufQYdtx1N55KseOuu1n/uc+z+dJLuy23484/Q3s7u/zi51QsWUL72rVUHHJIZ8HiZHU1FYcdRtUxr8USCaYcfxwkk9Rd8hMAKpaE5E75oYeSbmpi+61/pGyvvXqNlAVQPLsGCAWXARJVVV3dsvrZcqdl1SrK99+vX8t2HouyeLesqs7bSZ3caW/Hiospnjtv0rbcUXJnFC1YsIBDDjmE2267jbaoKaSIiIiIiMhIaLz/ASDUhElUV5Oo6l54dyRb7gAU77JzZ5el/kptyXTLmtZrXvVJJ5HavJltN97Ius9+Ftzp2Ly529DX22+7jeKdd6Zs//2p+cxnAJjyhhNy7i85dSoVhx5KatMminfemeJ584CuJE9q61ZKFy3Kuq6VlJCcMYPW558HQhHpzm5Z/RgxK7V9O20vvUTZ/gf0uWxcZ7esqMVO2HflgEcmm0i8rQ0rKaZ43lw66urwjo5ChzTqlNwZZWeffTYNDQ1cccUVhQ5FREREREQmsKaHHuq8321UpajwbmrHDqy0lERJyYjsv2zRvrQ8+2y35EtfUvX1JKqrseLiXvOqli7Fysp49atfw5JJZp53HqTTnUOPp7ZupfGBB6g++STMjIpDDmb3m29i+rvfnXefVce9HuhK6EAYNr54552j19G7mHJG+WteQ2tUByhRWUmitJTElCl0bOo7udP6r3/1uf1sMt2yElVVnS2KQnJnErfcaWvDiksomjs3/E1s3FjokEZdn8kdM3vdCE1HjsYLHGvOOeccKioquOyyywodioiIiIiITFDpxkaaV62i+tRTAEhUV5Ps0S0rvX0HiREoppxRtu++eFMTbS+/3O91UvVbSE7v3WoHQgKjaulSSKeZff75lEcjTHVsrANgx113Q0cHU046uSuGvffOmiiKm3LCCVhxMZVLj+n2fMWhhwJkHSkro/qUrn1ljm/RjBl0bOk7udNRWxuWnzuvz2XjMt2yElOqOp+b9Mmd9vbQcmfuXIBJ2TWrPwWVVwI+AvveAOw0Atsd08rLy3nDG97A73//e1avXs3uu+9e6JBERERERGSCaXrscUilmHraaZBIUrLrrllb7iSnjEy9HYCyfUN3ptann6Y0y/ee9tqN1F99NTPe997OgsId9fUUTZ+Rc5uzPvoRSvfai2lnnkFLNNJUR91GYD8a7/sLRTU1lA2whk3Jzjuz58q7Sc7ovt8pJ59E0+OPU5ajWxZA1XHHYSUleFtb5/FNzppFatNm0i0tACTKyrKum0lAFM+dM6B4M92ykrEh7BOVlbT3s4jzRORtbSRKopY7QPurrxY4otHX325ZaeDlYZwsmialL3zhC6RSKS6++OJChyIiIiIiIhNQ08MPQXExFQcfzE7f/x41n/i3zuK7qc6WO9tJjlC9HYDSPfeE4mJann6m1zx359UvfZHNP/85q99xGs2r/hliq9/aa6SsuLLFi6n5xL9hiQRFNaGgccfGOjyVouGv91N5zDGYDfyrZtHMmb3Wm/L617Pn7X/qrGuTTbKqisqlS4Gu4sZFM2fSsXkz6z75Kdadf37OdTs21GLl5SQGWNC6s1tW7L1LTp9GR/0kTu60t0NxcWfNpMk4YlZ/h0Kvc/dha2JiZunh2tZ4dNRRR7F48WKuu+46LrnkEor7aCYoIiIiIiIyEE2PPU754sXdEhOZliXbb7qZlifCEOn5EilDZSUllC5cSPOTT1J/zTVULl1KSVTHpv6q/6Px/geY8cEPsv2Pf2Tdpz/NHjffRGrLFsr2W9yv7Xe29qmro+XJJ0lv20ZVj65Vo2HmB5ZjiURnEejkzBl03LeB9rVrKZqXu8tVx8ZaimfPHnAyKtMtK9MNDKCopobU5i14KoUlk4N4FeNbpuVOprZU+4baQoc06lRQuUDOPfdc6uvr+eUvf1noUEREREREZALx9nZa/vlPyg96TbfnM4mepkceYet119H63HMj2nIHQt2dpgcfskxkxAAAIABJREFUZMPXv8HaD3+Eji1bqL34O9R+85tULl3K7M+dz/xvf5v2tWvZ+L3vkaqvp6ifCScrLiY5YwYdGzfS8Jf7IJGg4qijRvT1ZFOxZAk7//d/dSZVimbOIt3UhLe307FhA+7Zq5y0b6jt7EY0EIny0M0r3nKnqKamW3HpySYzFDpA0dy5dGxQt6xsbgZuG+b9jsQ2x5WPfOQjVFdX87Of/azQoYiIiIiIyATS8txzeEsLZQce2O15SyYpmjuX8iWhULC3tY3YMOgZVUuPITFlCjPPPZe2F1/k+WOXsWXFCqafdRY7/+QSzIzKI49g6mnvoP7qa/D29gEVGC6aPZuOujoa//pXyg7Yv9+JoZFUNGtm531vayO1dWvW5To2bKBozuwBb9+SyTDKWaygcvHssJ2OuroBb28iCEOhh1HfiufPp+nxv9Nw7700/OUvtNdOjpGz+uyW5e5vM7Odh3On7v624dzeeFRWVsbb3/52rrjiCh599FEOjSqxi4iIiIiIDEXLE08AUP6ag3rN2+P3t5CoqGDNme+kZdUqkiM4WhZA9amnMuWUUzAzimbNpOGee5n1sY92G3YcYN43vsG0t70NT6UoP+SQfm+/qKaGtjVraHv5ZWae+6HhDn9QOgszJxKhNc2GDb2STp5O015XR/WcgbfcAZj1/z7eOZoXEKs/tBH2G1hB6fEi0xoqOXVqt+fdPRoKPbTcqfn4x1j3mc+y9rwPA+FvcKcf/mDU4x1t/e2W9bKZ/XNEI5mEvvjFL5JIJPj3f//3QociIiIiIiITRPM/niA5cybFO83vNS9ZVYUlElQtWwZAYgRHy8rI1JSZcfbZ7HrZ//ZK7EBojVJx2GFUHnkkiagFRn8Uza6hbfVqSKWoPPzwYYt5KDK1gCqPPhrIPix3assWaG+naIAjZWXMOvdcKmJJsM7kzgRuuVP73e+y9iMf7T2jowOgs+VO+UEHscfvb2GnH/2Q0n337Tz+2/90O+3r149avKNtIDV3Ct++bYJZtGgRS5cu5Y9//CPrJ/AfmYiIiIiIjJ7mJ56g/MAD8xbqnXLc64EwytJ4VhR1R6K4mPKDDy5sMJHSvfem8uijmfmh0JKoo7Z3cd/26LniOYNL7vQULy49UXVsrMt6LL2tDQAr7koKJioqqD7lFEp3353U5s14WxvrPv1pNl92+ajFO9pUULnAvvKVr9DW1sZFF11U6FBERMY8M5tnZleYWZ2ZtZjZU2Z2bKHjEhERGStSDY20vfgi5QcekHe5ssWL2eXSXzD11FNHKbKRkWmxUn7ggZ1DhBdasqoqtFA6bAkkk1lb7mSSFEWD7JbVk5WUkJw2bUInd7y1lXRLS+/n29sBOrtlxSWjYek7Nm2CdJrW5/8/e/cdHmWVPXD8e6dPekIaHYKASkdFUJSigKKrIBZWLNhQd1m7axdRsa7YsGJHBRFEUX+sFRdEpSiIIE16SSNl0qbP/f3xzoQkpGcmM0nu53nmIZl53/eeKExmzpx7zo6QxxkuEZ/cEUKcLoRYKoQ4KISQQoipVR4XQoiHhBCHhBB2IcQPQogWs8nwjDPOoF+/fsyfP5+SkpJwh6MoihKxhBAJwCpAAOcAxwH/AtpGlzxFURRFqQdPtpZIMHbpUuexMaedVj4evaUKNBKOGnJSmCM5mtDrtYbP1YzlDiR8GtNQuSaGlBQ8uYeDdr1II10uZDXJHV+gcqea7XyGdu3wlZTgPngQAOfOnaENsgJPQQElP65qtvUiPrkDxACbgJsBezWP/xu4He0F/kloL/K/EUKEtjNYEN11112UlJTwxBNPhDsURVGUSPZvIFNKeYWUco2UcreU8jsp5ZZwB6YoiqIcTQhxjxBirRCiyF9x+bkQom8d53Tzf6Bb9XZWc8Xd3KSUFH76adAm+gTe3BuSU4JyvUhn7tkTYTQSO2pUuEOpljEtDXd2dZU7OWAwYGjXrpqzGseQmqo1VG6lfC6tcqfqaHnpqq1yR2tu7diyFQDv4cN4CgpCHKkm/8032X/99ficzmZZL+KTO1LK/5NS3iulXAT4Kj4mtE2ktwBPSCkXSyk3AVcCscClzR9t4/z973+nS5cuvPrqqziqyUQqiqIoAEwAVgshPhJC5AghNgghpovaGgooiqIo4TQSeBk4BRgNeIBvhRBJ9Tj3LKB9hdv3IYox7GyffELm3fdQ8OGHQbme57A/uZOSHJTrRTpT1670Xv8b1ipj3yOFIT0dT1Y2ZevX487MLL/fnXkIQ2oKQq8P3lopKa17W5bLDT5f+Tas8vvdtVfuADi2bi2/z/XXXyGM8gj77xvB68VXWtos6zUkuRMnhLhbCHG2EOLotuvh0R1IB74O3CGltAMr0H6JtAg6nY67776bvLw8Zs2aFe5wFEVRIlUG8A9gFzAOeB54AvhndQcLIaYJIdYJIdbltuIXOoqiKJFKSjlOSvm2lHKTlPIP4HIgBTi1HqfnSSmzKtxcoY02PFz79pE16zGAoPUCCby5DzTYbQuEwRDuEGpkTEvDfeAAe6+4ktwXXiy/371vP6ZOnYO6liElBc/hw0ifr+6DWyDpr4CpujWrvHKnluSOs0Jyx9kMyR3p9WLfrA0cj8TkjhWYBXwB7Pd/avq1EOIpIcQUIUQfIURzVwIFuk9V3cSYXeGxSiL1xf71119P9+7deemllyhtpv/5iqIoLYwO+E1KeY+Ucr2U8m3gBWpI7kgpX5dSniilPDElpW2UpiuKokS4WLTn8vrsifjE/35jlRDiwhDH1awcW7awY8RIXAcOkP/+++DxYD3hhOAldw7nIkwmdLEtpktFq2ZIT9cqTdxuXHv3lt/v2r8fU9e6+yI1aK2UFPB48BYWBvW6kSIwFctnr5LcqaOhMmjJU31yMrqYGJw7Qp/cce7ciSwrAyIzueMGMtEaWQogGTgTrd/Ne8BGoMSfOHnDXyo/XAgRF+ygqyGrfC+quU87MEJf7Ot0Oh588EEKCgp4+OGHwx2OoihKJMoE/qxy3xYguK+MFEVRlFB5HtgA/FzLMSXAHcDFwHjgO+AjIcRlNZ0QqR/e1sT+++94srMp+vL/KPnf/4gaejLRpwzDfeAAPv+bwabwHj6MITm51jHoSvMxde4EgD4lGff+/YA20cybl4exc5CTO6na+9vWujWrvHLHWbVyp5ZtWUnaLlDpdmNITcHco0ezVO44/thU/nUkJnfypZSd0EopxwB3Ah8AmwEvWkLFAgwGrkJ78v4fUCCE2CmEWBTMwP0CnamqVumkcnQ1T8S74oor6NmzJ6+++io2my3c4SiKokSaVUDvKvf1AvZWc6yiKIoSQYQQs4HhwCQppbem46SUh6WUz0gpf5FSrpNSPgi8htZUv6ZzIvLD25q4M7W3MAUffoh77z5iTh+BuWdPkBLnzl1Nvr4n9zD6NtJvpyWIGTmSrh9+QOLf/44nNxef3Y57/z4ATPWYaNYQgbHwnpzWmdzxNaJyRxcVhYiKAsCYkoqp5zHNMjHL/sfG8q8jMbkDgJQyzz+d5Bkp5eVSyv5oE61OBK4BXgR+BIo4UuXTHZgYvLDL7UZL8IwJ3CGEsACnAT+FYL2Q0ul0PP744xQVFXHzzTeHOxxFUZRI8ywwVAhxnxDiGCHERcBNwEthjktRFEWphRDiWeDvwGgpZWOyF6uBnsGNKnzcmYcA8GRrn0XHjByB+RjtxwtGRYHn8OE2MymrJRAGA1GDB2PyV+m4Dx7EtU+r4DF1CXLPHf8WJG9+XlCvGykCFTrSYa/2fmE6OrkDR/67GFJTMHXqjDcvLyhVcrVxbPyjPNnmKykJ6VoBQemRI6V0SSl/8zdMu1lKOUJKmYjW/PIC4GFgaWOuLYSIEUIMFEIM9Mfbxf99F6nNQHsOuFsIcYF/tOI7aOWcwWk338wmTZrEqaeeygcffMCmTZvqPkFRFKWNkFKuRZuYdTGwCa0P3ANok1gURVHaPCHE0XsSwkwI8TzaFNvRUsqtdR1fg4FoW3NbBU9mFob27QEw9eiBqVMnTF06I0ymoPTd8fi3ZSmRJbA9y7V/P659WtGxMciVO4H+Mp68/KBeN1IEtmX5HJVHi5dPyzJW/xQY2JplSEnB2FGbDVVxclkoOPfswTpwAADeZqrcCWlbcSnlHmAP8GkTLnMisLzC9zP9t3eBqcBTaM2eXwIS0TL7Y6WUxU1YM6xeffVVBg8ezLRp0/jppxZXgKQoihIyUsovgS/DHYeiKEqEKhVCbAPOllLuD3cwQoiX0CZkTUBr1RBopVAipSzxH/M4MERKeYb/+yvRen2uB3zA39Aa59/VzOGHjDszk6gTTgAhsA4aCGjVHaYePZqc3JFuN96CApXciUDGzlqVjnv/Adz79qNPSkIfExPUNXQxMQijEU/e4aBeNxJIjwf8U8COrtwJTMuqvnJHX165k4qxY0cA3IcOYe7RIySxektKkWVlmLp1B5pvW1bkzozzk1L+gLa1q6bHJfCQ/9Yq9O3blylTpvDOO+/wwQcfMGXKlHCHpCiKoiiKokQ+N9A1EhI7fv/w//ldlftncuS1e3ug6jus+4GuaH09twNXSynfD1GMzUr6fLizs4lrn07q7bdXeszc8xjKfv4FKWWjmyF78gtASgyq507E0SclIaKicB3Yj2vfPkydg7slC0AIgT45GW8rrNwJVO1AzZU7umoaKkOFbVkpKRg7+Ct3Dh4KRZgAeHJzADB16wpEXs+dz4ADoQxEqey5554jJSWFm2++mYKC+kyLVBRFURRFUdq4XdTyoWhzk1KKGm4PVThmqpSyW4Xv35VSHi+ljJZSxvkbJbeKxA5oW6Zwu8u3ZVUUNfgEPLm5uPcePSegZOVK9k2bhs9uP+qxytfXGumqyp3II4TA1KkT7v0HcO3fhzHIY9ADDElJeFphz51AM2WopXKnmobKAPp2R7ZlGVJSwGDAfSiUyR3/v8P0dITFgq80tP19AuqV3JFSTpRSDgl8L4QYELqQFID4+HjmzJlDXl4eV111VbjDURRFURRFUSLfQsAqhDg93IG0ZdlPPkXhJ0uqfczj7/NhrC65c7L2dqv0l9VHPVb0xZeUrljJ4Vdfq3Vt72FtO45K7kQmY+fO2Df9gSczC1OXriFZQ98uqXVW7lRI7hw9LavmUegAxnTt35shPR2h12NMT2+e5E5KCrro6Iir3KnqUyFEUlAjUY5y8cUXc+655/LZZ5+xYMGCcIejKIqiKIqiRLZngH3Ai0KI4DbzUOrNtnQpxd9X3YmmCYxBry65Y+rWDUNaGqWrfznqMbt/0EreW2/h3LUb0N7sZj3yKDnPzEb6e5F4AsmdFjASvi0ydeqEN/cwxvbtSbggFMOkwdAuGU9e66vcqZjckU5HtY/VVLkTP+F8us57D2NqKgDGDh1CmtzxVvh32BKSO12BRUIIfX1PEEIEf1NhG/Duu++SmprK9ddfz65djZkcqSiKoiiKorQR04EX0CbWrhNCnBLmeNocKSVemw1foa3ax921VO4IIYgeejJlq9egtRXVeEtKce3aReKll4LHQ9GXXyJdLvbfcCMFH3xA3ty5HLrrbrzFxZStWQuAXlXuRKSY0aOJGTGCrvM/LO/9EmyGdkl48/Iq/R1qDSr13DmqcifQULn6yh2dxULUSSeVfx/M5E7Zb+spXLSo0n2e3FyE0Yg+IaFFJHc+A0ag/fKokxDiVGBtI9dq05KSkliwYAFlZWX87W9/w+3/i6soitJcInG0rqIoilKtx4H/AFFAT2CFEGKVEOLfQoihQoio8IbX+vlKS8HjwVtUVO3jnqxMRFQUuri4ah+POnko3vx8Sn74ofw+x+bNICUxo0Zi6tIF544dlK5eQ+lPP5F2332k3HILRZ9/zo7hp2H77DMSL7sMndkcih9PaaLok4fQ+bVXMaalhWwNfVI7pNuNr6QkZGuEQ6VtWTVV7tSQ3KnK2KEDnpyc8qRQU+S/+y7ZT/+n0n2e3Fz0KckIIdBFR0V8cudytM71Nwghrq3tQCHEVcC3gKoNbKRRo0Zxzz338Oeff3LhhRfi85ddKoqiNJNSIcQmVYGpKIoS8ZYCe/xfC7TX+sPQkj6rAJsQYosQYr4/4TNWCJEanlBbJ6+/Yqem5I77UCbG9u1rnIYVM+J0DGlpHLjxH2TOeAgAxx8bAbD07Yu5V0+c27drCR+07SbJN1xPt48XEn3acNLuuZu0++4N8k+ltCSGZG0yVGCLXnMr+/VXDtx8S/lWwWDxVajckdVV7ggB+vptLDJ27AD+yXVN5dq3D5/Npo1q9/Pk5pZvjYz4yh0pZQkwASgB5gghhlU9RmhmA28AZupZ5aNU76GHHuL8889n6dKl3HzzzeEOR1GUtiXSRusqiqIo1ZBSTpBS9gASgNPQtmnNRaugtwN6oDdwCVrCZxkQusYTbZC3sFD7s6bkzsGDtW7HMbRrR4//LiNxyhQKP/qI4uXLsf+xCWOnThgSEzH37IVr3z7K1v+GsWsX9LGxAFj79aPznDkkXXllo8eoK62DPklL7njzw9NUuWTlSoq/+gpfcXFQrxuYiAXgO2palgthMtX7736wxqFLKcun23ltR7ZiVkzu6JsxuWNo7IlSym1CiCuAT4DFQogTpZSHAIQQcWjd+segvSm4UUr5djACbqt0Oh2LFy9m+PDhzJkzh/j4eB599NGQrllQUMBPP/3Enj17MJlMDBo0iMGDB6PTNbbgS1GUFmoX0C3cQSiKoij1I6UsRqvUWRW4T2jveo4BBlS4DQQ6hSPG1spr05I7sqys/A1ngJQS1759xA8aVOs1dFYraXf9m7I1qzn077uQTicxZ4wGwNyrJ/h8lK76idgxZ4buB1FaLIN/7He4mioHEpy+sjL08fFBu650VajccTgrP+Z219hMuTqBnleerMwmxeTNy8NXpo059+bnY2jnr5rKycV64omAVrnjLW2eLXKNTu4ASCk/E0I8CjyANkHrNLRmy5+hfSqQDVwgpfy5yZEq6PV6vvnmG0499VRmzZpFUVERzz33XFCTLUVFRcyePZvFixezZcsWvF5vpcfbt2/PtGnTeOCBB9DXs+xNUZQWbyEwQwhxupRyRbiDURRFUWDlypUMGDCAuBp6t1Qlte6qO/y38u6fQoiE0ETYNgXe2IJWvWNITkZ6PFoPFLsdX0kJpq5d6ryOMJlo/9jjZM+ahbl3b5IuvwwAc69e2gEeD9Y+fULyMygtm75deCt3AlsTA0mPYKnUc8dReVuWr0oitS46f9LJW9S06iLXvn3lX3sKCjD7Y/HabEe2ZUVF4ysN7n+LmtSZ3BFC3ApsADZIKQuqPi6lnCGEOAE4G/gCOAGtFPRXYIKU8mBwQ27bYmJi+OWXXxg1ahQvvvgif/75J59++ikxMU2bdrl582ZmzZrFZ599RllZGSkpKUycOJHTTz+dnj174nA4+Omnn/joo4+YOXMmCxcuZNmyZXTt2jVIP5miKBHsGeAqtNG6p/q35iqKoihhNGLECIQQgddiPYQQDwK/o71m31vf60gpC+s+SqmvilszAsmd3Jdeovj/ltH+iccBMHapO7kDYO3Xl24L5le6z9SlC8JkQrpcWI4/PniBK62GITERAM/hcFfu2Os4smEC07KE2Yx0HN1QuSGVO3r/e2dvcfXbJ+urYnLHW+DfkllhDDpolTuyrAzp8yFCvAOmPld/Bq0h8mEhxB4hxBIhxAwhxHlCiMAz06XAX8BotMTOfOA0ldgJDavVysqVK7n88sv57rvv6N27N59++mmjrvXpp59yyimn0L9/fxYsWECfPn2YP38+WVlZfPzxx/zrX//irLPOYsKECTz11FPs2rWLBx54gO3btzNkyBC2bdsW5J9OUZQIpEbrKoqiRJikpCSklOzZswe0198z0Nol7BJC5AshlgshnhVCTBVCDBRC1P+dj9JolSp3/Ikex+bNuPbupWz1akBL0DSWMBgw9egBoJI7SrWE0Yg+Ph5vfpiSOwVaPUiwK3d8/sodfVzcUZU70u1uUOWOMBoRUVH4ipv2eaW7UnJHq5Ty5OYCYEhOBrTkDgT/v0d16pPc+T+0RmsC6AKcDzwILAF2CyHy/F/v8h/zPnCllNJR/eWUYDAajbz33nvMnTuXsrIyJk6cyLBhw1i2bFmd07T27NnDXXfdRUZGBhMnTmT9+vVMnDiRDRs2sGbNGiZPnlzjVi+9Xs/DDz/MwoULKSgo4MwzzyQnJycUP6KiKIDDERFPpWq0rqIoSoQ5fPgw+/btY+nSpaC9Vl8C7PY/nACcDtwEvIlWUV8ihNgghHhHCHGrEGKUECIxHLG3Zr6KlTv+r937DwBg+/wLEAJjp6a1OYoaNAhzr17oE9SOOqV6+nbt8OSFa1tWoHInuE2EpVNL7uji45D2qg2VG5bcAa16p8mVO3v3YfCPtQ8ktdxZ2gQuQ4o2iLA8udMMTZXr3JYlpTwXQAiRBAziSPO1gcCxQCIwCpD+U6YAFwshtqKVhgZuG6WUucH+Adq6a6+9lgsuuICbbrqJRYsWMX78eFJSUjj55JPp27cvqampGAwGMjMz2bZtG+vXr2fv3r34fD7S09O54447uOeee0hKSmrQuhdccAFz587lqquu4qKLLmL58uWq0bKihMA555xDbm4uGzZsCOe/saVAP7SmysJ/GwYM9T/uE0L8hbaFdz1HtvKqzK+iKEoIderUiU5aoiBTSnkhgBDiVLSeOmmAF8hHS/YYgf5oz+eXB64hhNgvpezWvJG3Xt7CQm0ks5T4ioqQPh/ug9pmBtfOnRg7dEDXwDehVaXefVelyUGKUpUuJqbZJjRVJKWs1FA5qNf2b8vSx8Ufde2GNlQG0MXGNrlyx7VvH+YePfCVlODxJ3fsG39HGI2Yj9Eq7HT+LWARkdwJkFLmA9/5bwAIIUxAXyonfPoD8f4/+wOXBS7RkPWU+ktKSuL999/nmWee4fnnn+eLL77gm2++4Ysvvqh0nMFgICMjg7///e9ce+21nH766U16s3jllVeyYsUK3nrrLR566CEefvjhpv4oiqJUYLPZWLVqFcOHDw9r8lRKOQFACBGL9rxecdJKX7SKnt7+28WB01DP+YqiKM1KCNEP+BIoAM4FvpJSeoUQOrTn7AuB64DkCqd1bvZAWzFvoQ1DejqezEy8tiI8ubmVGsHWt99ObXQmEzQxQaS0bjqrFZ89uD1v6kPa7eVJmKrVNU2+ttu/LSs+Hm+VSWBVJ9PVhz42Fl8Qeu7EjT8b1/795T13ytatw9K/PzqLBQBdtFbgHlHJnepIKV3Ab/5bOSFEdyonfAahfnGEXFpaGo899hiPPfYYPp+PLVu2kJubi8vlokuXLmRkZGAK8i+C1157jRUrVvDss88yZcoUevfuHdTrK0pb9vrrr+N0OrnqqqvCHQqgRusqiqK0AI8BscB4KeVPgTullD60ysr1QoingFeAycCrQOMaNyrV8tpsmDp18id3bLj37wfAlJGBa9euJvXbUZT60lmtuHOym33dij2ngt5zp7xyJ67864CGNlQGrXKnYrwN5S0uxuf/9+5ITMSbn4+vtBTH5j9pd801R9aJpG1ZjSGl3I2257f8l4UQInhD7pU66XQ6+jTDeESDwcA777zDaaedxg033MB3332ntmcpSpDMmzePxMREJk+eHO5QaqRG6yqtkdPrpMhZhM1po9BZSJmnjBGdR4R0zbWZa0kwJxBvjifeHI9Jrz6VVxplOFBWMbFTlZTSBlwqhLAA1wOfNVdwbYG3sBDzcceii47GW2TDdUDrtxM3fjyH58yp1xh0RWkqXZQVGeRpVfURyuSOdLlAr9emT1XtueN2I8wNr9yp2BC5oTz+vrOGtHT0iQl4cw9TtmEDeDxEnXhi+XEtPrlTHf8vEqUVOvXUU7n44otZuHAhS5YsYdKkSeEOSVFavM2bN/PHH39wxRVXoNfrwx1Og6nRukpLUuYuK0/k2Fw2HJ6jG5kv37e83tcb1WVUg44HKHGXUOIu4UCJ9kbQarBqiR6TluyJMqq+5Uq9mYQQeimlt47j7gYmADcCX4U+rLbBa7Ohj49HFx+Hz1akNVMWgvgJ55P/3ntYBw0Kd4hKGyDCtC2rUnKnNNg9d7StV8JiOXpalsuFLjamQdfTxcbiLWl8z53y5E5qCobEJJw7dlC2bh3o9ZX+neurSe6UrPwRQ0oylmOPbfT61VH9EJSgeO6551iyZAmPPvoo48ePx2q1hjskRWnRnn/+eQBuueWWMEeiKK1PqbuUQmchhY5CbE4bLp+r7pOamd1jx+6xk1WaBYBJZyLeEk+COYEEcwLRxugwR6hEqA1oU7ImUqGasjpSyu1CiGJgSHME1hZIn09L7iQkoI9PwFtUBAf2Y0hPx9SpE73XrA53iEobobNGhSW5E2gqDAR9felyoTOZ0FksSKcT6fMh/DtGGtNQWR8Xi6+o8T133NnatjdjWhr6xES8+QWU/bIay/HHo4858js6ULnjtWlrSa+XrBkzMHbsSNd57zV6/erUuX9GCOETQhwM5qKhuKYSXunp6UybNo0NGzbw3nvB/UuqKG2Nz+fj008/pXfv3gxSn/ApSpNIKSlxlXCg+ACbD29m1cFVrM1ay46CHeTacyMysVMdl89FblkuOwp2sDZrLasOrmJz3mYOlhykxFWCtkNSUXgNbaLhs0KIbrUdKISIAWKAho1MVWrkKykBnw99fAL6uDi8NhuuAwcxdewY7tCUNiZcDZUDlTvCZArBtiynVrlj1RoVywp9dwKJn4bQxcQi3e6j+vfUxJ2dzYFbbsVbXAyAJ0cbBG5ITUWflIh0OLCvX0/MyMrbuPWJiRg7dODwa6/h2rOHkuXLcR86ROJllx21RlPVtzmKCPrKobmmEkYzZ84kKiqKF198kaImZEGVo/l8Pr744gtuuukm7ruj9jcmAAAgAElEQVTvPrZs2RLukJQQWrJkCbm5uVx66aXhDkVRWhyX10WePY89tj1szN3IqkOrWJe9jr8K/yLXnovb1zrGB7t97vJkz7rsdaw6tIqNuRvZY9tDnj0Pt7d1/JxKw0gpFwBLgY7AWiHEdUKImj7OnoH2ejyvhseVBgq8sdUn+JM7RVpDZWNnNVdGaV66KCu43Uh38/4uCPwbMLRPD0nPHWE2ozNryZ2KW7MaW7kD1Lt6p3TVTxT/97+UrVkDgCc7G11cHDqrFUNiYvlxsWeOqXSeMBjo/MZc8HrZe/kV5L44B0P79sSeMbpB8dZHfbdl6YUQnVEJGaUWSUlJXHfddbzwwgt8/PHHXFOhS7jSeJmZmZx77rn89tuRoXRPPvkkV155Ja+//nqL7Mei1O7VV1/FZDIxffr0cIeiKBHLJ32UucsodZdS6iml1FVKsbsYl7dlVOIEm8fnId+RT74jv/w+k95ErDGWaFM00YZooo3RRBmj0Ak1+KCVmwy8AVyKNg3rMSHEl8CvaCPSOwDnA0MBCXwSpjhbHa9NazEa6LnjPnAQ6XBg6t49zJEpbY3wt8jw2e3oG5j0aApvoQ1dbCz62Dh8ZcFtIOwL9NwJVO5UTO40YhS6LkZL7niLSzCkpNR5vDsrEwDHn1uIPeMMPDk5GFK18/T+5I6xaxfMvXoeda45I4Ou773Lgen/wrltGym33oowBL9DTn2vmAzsCeK6qna4lbrzzjuZM2cOb7zxBpMmTSIhQQ3MaYrs7GyGDBlCZmYmN954I3feeSfZ2dnccccdvPXWWxw4cIBly5apCWWtSHZ2NitWrOC0004jKUlVyittm5QSh9dR3n/G7rZT5imjzFNWbdNjpTKX10WeN488x5HCDIHAYrBgNViJMkRhNVqxGrSbRW9BCPU5XksnpXQAlwkhvgJmAt2AK4DLKxwW+B/9K3BfswbYijm3bQPAkJSIPj4e6XCgi40l4aILwxyZ0tborFoTfp/dro0Odzjw5ORg6hLcaW0+/3YoX2kphZ8swZOTgz4hAV1UVNCndZVX7lj8lTv2KsmdBo9C1xowu/fvI//tt0i7++7y/jjV8WRpPXYcf/6pfZ+TgzE1DQB9ovaaPW7MmBp/j5p79qTb4kUUffl/xJ/3twbFWl/1Te6o3/RKvXTs2JGLLrqITz75hC+//JIpU6aEO6QWy+v1Mn78eA4dOsTcuXO5+uqrAejevTsrVqzgsssuY/78+Vx33XW8+eabYY5WCZb//Oc/uFwubr/99nCHEtGEEPcCs4CXpJSqxKkF80kfDo+DMk8Zdo8dh8dR/qfD68AnfeEOsVWRyPJkWT75lR7TCR0WvaU8+VMxCWQxWFTFTwsjpZwnhPgQOBs4CxgApPkf/gttBPpbUkq1h68JfGVlHLr7HqwDB5L32mtY+vfH0rcvpau1rRvJ//hHpS0bitIcdFH+yh3/1qismQ9T8sMP9Pr5p6Ct4c7JYeeZY+jy5hu4c3LInjUL9Hosffqgs1rx5OYGbS3QeuwIkxFhCfTcqbotq4Gj0OPiALB98SVFn39O7NhxxJw2vMbjyyt3/O0x3Dk5RJ98MgDmXr2IGT2ahIsvrn3NmBgSL6n9mKaoM7kjpdQJIR4DHpNSNn5WmNJm3H333SxYsIB33nmH8ePHk6h+oTXKvffey2+//cZtt91WntgJ0Ol0vP/++xw6dIi3336bMWPGMHny5DBFqgSLz+fj/fffp1u3bpx99tnhDidiCSGGAtcBG8Mdi1I/Xp/3SPVNhVuZp6zNbqOKRD7pK6+Mqo5Zby6v8ql60+vUFuFI5B+H/oX/poRA2bp1FH/9NcVff40wmejw+GMIg4GYUSNxZ2WSOEX1z1Oan86/LUva7bizsrB9/jl4PEivFxGklg7u/fuRLheOrduQHo92p9eLPiEeXXQUvr3B77mjMx1duSM9HnxlZbVW3VRHF6NV7th//x04krypiSdTm17pycrCk5eHJzcXQ2oqAPqYaDq//FKD1g+F+lbu3A1MFULMAN6QaiSDUosBAwZw2mmnsW7dOlauXMl5550X7pBanJ07d/LCCy9w/PHH8/TTT1d7jE6n45NPPuHYY49l+vTpjB8/njh/BlppmRYuXEhWVhYzZ84MdygRSwgRD3wAXAM8GOZwlCrcPveRPjjuUi1R4C7D6a3fJIqW7tslqbzxVAYAk4cN5dp/7+LMiTlhjip4nF4nTq+TQmfhUY+Z9WaijFFEGaKINh7p72PUNV+vB0UJh7L160Gvp+N/nkYXF4e5Rw8ALL160X7GjDBHp7RVFXvu2ObNA3/yxWe3o/cnNZrK6x977snJ0ZI7ej1IqVWqGY1Bb6jscznRx8ahi44pX7c8DikxpCQ36HqByh33vn3a9fzJm5q4s7IwHdMD1187KV21CjweDGmpDf0xQqq+yZ1SIB2tKdt0IcQdUspvQheW0tJNnz6dSy65hEWLFjF69GhigvQk0lZMmzYNt9vNm2++WWs/naSkJJ577jmmTJnCP//5T+bNm9eMUSrB9vzzzxMVFcUtt9wS7lAi2evAIinl90IIldwJI4/PQ5GriCJnESXuEopdxS0qiXNev/MothWH7PrZB63Muhlm3Rz8a8fGx7L0j6XBv3ATBBI/BRRUut+sNxNrii2/xZniMOiC30RSUcLFvn4Dlt69iVMVt0oEKe+5U2bHtmgxwmjUxn6XlQUvueOfjOXJyUF6vRg7dCD1jjswdetK4eLFQR/FLl1uhMmEtV9fDGlpFC78iLizxuE5fBgAfXLDkjuBhsoB7qyakzveklJ8xcXETr6EvL92UvLDDwDllTuRor6/XXugNWS7BugH/NffpO12KaWayawcZcKECaSkpLB27Vp+/vlnxowZU/dJCgCrV69m+fLlTJw4kaFDh9Z5/KWXXsprr73G/Pnzue222xg0aFAzRKkE25YtW1izZg0XXHCBqsCqgRDiOuAYKjcFVZpRmbuMXHsuefY8il3FyBY8H6HYVszyfcuDdr3Jw4Yy+uBiHuM+urCPfXThXmbxfcdJLPj5l6CtAzCqy6igXi+UnF4nTruTw3btxbdAEGuKpZ21HSnWFKKMUWGOUGkOQoh/AHcC7YHNwC1SypW1HN8PmAMMAfKB14BHIm33gPR4sG/cSMLEieEORVEqCfTc8RYW4LXZyitOZBCraTyByp3cHKTXhyElhbhxY/3rR+ErK0NKGbRG/dLpRJhNCKORxClTyJ09G8e27eXJHUO7BiZ3oqNApwOf1tvPU8u2rMBj5l69MfXoQdE33wJgTEur8ZxwqFdnPClljpTyRqAvsBStwfI44HchxMtCiLpnhyltislk4uqrr2b79u1899132IOcuW3N7rrrLgwGQ43bsaozd+5cdDodt956awgjU0LpwQe1IpQHHnggzJFEJiFEb+AxYIqUsl5NWoQQ04QQ64QQ63KD3NSvrSlyFbEhZwNrstaw27abIldRi07shMIZBxczl2l0Yy86JN3Yy1ymccbBxeEOLaJIJEWuInbbdrMmaw0bcjZQ5CoKd1hKCAkhLgGeR3sOHwT8BCwTQlQ7tkcIEQd8A2QDJwE3oSWGbmuWgNE+pXfu2l3ncY5t25BlZVgHqw/WlMgS6LnjPqQlJYwdOgAEdatUoHLHnZOjjQWvME5cZ40CrxfpCl5fPel0ovOPO0+8+CKE1UrBhx/iOaxNhGzotiwhBLpYf/WOXo+7lm1Zbv+kLGP7dDo8Nqs8YRVplTsNGnsgpdwupZwAnA6sRav8uR7YIYT4txCiYS2qlVZt2rRp+Hw+Vq9ezZo1a8IdTouwZs0aVqxYwXnnnUdGRka9z+vVqxcTJ07kf//7Hz/4ywSVliM7O5vPP/+cYcOG0b9//3CHE6mGAcnAJiGERwjhAUYA//B/b656gpTydSnliVLKE1NS1GcQjZVZksn67PXV9llRjnhCfw/RVH7RHE0ZT+jvCVNELUOhs5D12evJKq2914HSot0GvCOlnCul3CKl/BeQCdxYw/FTgCjgSinlJinlYuBJ4DYRrBKAWrgOHGDPxRez69xzKVxce3LWvn4DAFGqalqJMIGeO+7MQwAY2/uTO0H8wN1bENiWlVupuTBolTsQ3GSSdLkQJu3lnj4hgahBA3Fs2YLnsPYBnqFduwZfM7BFzdq/P+7sbGoqDgxU7hjS07EOGECHp54k+vTTKiW0IkGjZlpKKX+UUg4FJgO7gTjgcWCrPzuvKGRkZDB27Fg2btzIzz//jCuImdvWasaMGej1ep588skGnzt79mwsFgt33HFHCCJTQunhhx/G6XSqqp3afYq2LXhghds6YIH/a/UEEwI2p43tBdtVlU49dPIeaND9yhESybb8baqCpxXyf/B7AvB1lYe+Bk6p4bRhwEopZcV3oV8BHYBuwY4RoHDRIvLfe4+cZ2aze+IFePLysA4eROZ992P77LMaz7Nv/B1Damp5VYSiRIpAcsWdGcLKHf+2LF9REb6SksqVO/71g7kNzOdyIUxHakmMnbvg3rcP7+HD6KKiytdsCJ2/FUL0Kacgy8rwFVX/e8idmQVCYPQnsOLOOosur78etMljwdKkjnZSyoVCiE+AfwL3oT3hfiiEuBm4TUoZ3E3mSotz7bXX8vXXX/Pnn3/y66+/MmzYsHCHFLGysrL47rvvGDFiBD38kxYaomPHjlx66aW89dZbLF++nFGjWk5PhrbMbrfzwQcfcPzxxzNu3LhwhxOxpJSFQKXSESFEKZAvpdwUnqhav9223a02sSNnxMGb5wfvgjNiq71bACODuQ5a7D8E9YrhJ5HsKtzFwNSB4Q5FCa5kQI+2xaqibODMGs5JB6pmRbMrPHbUfikhxDRgGkCXLtXu9qpV/rvv4tzxFwhB7JgxpN52K8YOHdh39TVkznwY68CBmLp2xbl7N/q4uPIKAcfmP7H06dPg9RQl1ALbsjxVt2WVBn9bVkCl5E50iCp3zEcKtU2dO+EtLMS5Zw/6Bm7JCtDHxGBITcXc8xhAa6qsj48/6jh3Vib65HaVkkuRqMnjCqSUHuB5IcTbwL1o+2KHAquEEB8Dd0sp9zR1HaVl+tvf/kZiYiI7duzg559/ZsiQIegjLMMZKZ588kncbjd33nlno6/xyCOP8P777/Pggw+ycmWNfQqVCPLkk09is9mYPXt2uENRlKPoRet9vhYzi4LaUDl1ybf0vvs/6O1HJoZ5rWa2PXEHORNreg/bOKO6jGL5NUG9ZERQU7RatapZYlHNfXUdX9392p1Svo42TZETTzyxwRnp7p9+iq+4GBkY4+zX4emn2HX+BPZedRUxp55K4eJPiB0zhk7PP4evtBTXrl1qSpYSkYRejzCZjlTudAzFtqwC9AkJ5UkeQ+rRlTtBTe44nQiTsfx7Y6fOANg3/I75mGMadc3YcePw2goxpqcDWqWTpXfvSsd4CgooW7MWU4eOjYy8+TRqW1Z1pJRFUsq7gV7APLQn34uBLUKIJ4QQ1X+kpbRqFouFyZMns379enJzc9m4cWO4Q4pIXq+XefPm0bNnT8aOHdvo63To0IEJEyawatUq1q1bF8QIlVBwOBzMmTOHLl26MHXq1HCH0+JIKUdKKaeHO47WrE9yH45NOpZESyI6EbSXDK1SzsQz2fbEHTg6piGFwNExLSSJndZGJ3QkWZI4NulY+rRTFRCt0GHAi1ZxU1EqR1fzBGTVcDy1nNMkQq9Hn5BQKbEDYGzfns6vvoKhXTKFHy9CFxODY5NWLOrYuhWkxNLn+FCEpChNprNa8ZWUANrfZQBfWWnQru8tLMTcq1f595UbKlv96wUnuSM9HvD50FWo3DF27qStYbM1qt8OQNJlU0j55z8x+JM7xcuWseeSyWwdNJidZ53NgZtuZu+lU/BkZ5PSAgbXBP2VmpTygJTySmAwWqd7M1qH+x3BXktpGaZOnYrT6WT//v389NNPNTaqasvmzZtHXl4eV199dZPHBT722GPo9Xruu+++IEWnhMqjjz5KXl4e9913HzqdeuOsRB6d0JEenc6AlAEM7zicwWmD6ZnYk/bR7Yk3x2PUGeu+SBuSM/FMfvl5Af/b+z2//LxAJXaqMOqMxJvjaR/dnp6JPRmcNpjhHYfTP6U/6dHpQRuXq0QO/3TDX4ExVR4agzY1qzo/A6cJISxVjj8E7Al2jHWJGjyYbgs/otfqX2h3zTW4Dx7EW1yMY/OfAFiOV0lJJTIJf/WMiIpCn5AABDHZ4vXitdkwV6hyqZjcCawdrEoh6dSqYgMNlQFMFbZgGpIbty2r/PyUFNDrsX22FM/hwyRMmoSpe3ecO3cijEY6vTSH6KEnN2mN5hD0+lchRDvgeOA4YCtao8sU/01pg0466SSOO+44/vzzT3r27Mm2bds49thjwx1WRHnllVeIjY3ln//8Z5Ov1aNHD8aOHctXX33F1q1b1X/rCFVaWsrLL79M9+7dufbaa8MdjqLUSSd0xJniiDPFVbrf7XNj99ixu+04vA4cHgd2jx2n14nT68QnfWGKuG6jurTM3mSx8ZFZDK0TOsx6M2a9GavBisVgwaK3YDVasRqsKhnYds0G5gkh1gCrgBvQmiO/CiCEeBwYIqU8w3/8h8AM4B0hxKNouwLuBmbKMH1CKIRAHx+PubdWpeDcsQPH5s3oU5IxpkXWKGRFCQhUz+gT4hEWCwiBDFKyxVtUBFJi6txZu7bXW55Aggrbsir0+AkkegJxNYTPP5inYs8bfWws+vh4vDZbg8egVyX0eozt2+NzOuj63rsYO0b+FqzqNDq5I4RIR0vgHF/lVvW/rPoYpo0TQjB16lTuuusuxo8fz48//kjv3r3VJ3R++/btY926dZx//vnExgbnBfusWbP473//y/3338+iRYuCck0luB5++GEKCgqYPXu2qtpRWjSjzojRZDwq6QMgpcTlc+HwOHB6nTg8DhxeB06PU/vT68Tj84QhaoLab6eqUV1GhfT64WLQGcoTNya9CaveitmgJXMsBgsmnUn9bleOIqX8yP/h7/1Ae2ATMF5Kudd/SHugR4XjbUKIMcBLaFMRC4Bn0JJEYRXoxeHctg3Hn5uxqqodJYIdSe4kIITQtmkFqaFyYAy6PjFRG4Hu8VR6/tdFRQOVK4UO3nobGPR0njOn8eslVG52bOzcGa/Nhr6R27Iq6vjsbPQJCS02sQP1TO4IIcZSOYFzHJBQ9bAKXzuBbcBmtCfwwE1poy677DLuueceDh06hNFoZO/evXTr1i3cYUWE2bNn4/P5mD49eK1DBg4cyPDhw/n88885dOgQHdSIzoiSm5vLyy+/TM+ePbniiivCHY6ihIwQorySoyYen6c86ePwOMq/tnvsODwOvNLbjBG3bXqhx2KwEGWIwmKwlCdtLHoLFoNFNTxWGk1K+TLwcg2PTa3mvj+A00McVoMZ0tPRxcVR/M03OP/aSaxqpqxEsEByx+CvqBHRUUHbluUt1Mag6xMTtWbN7sof1FQ3LcuxfRvS7kBK2eAPAjw5OQAYUipXyhk7d8KxaROG5KZvErL269fka4RbfX9L/5cj3ekr/p/wAn+hJW4qJnJ2SKlejSlHdOjQgbFjx/Ltt98yffp0Vq1apZI7aJ9qL1y4kIyMjKCPLn/ooYcYPXo0M2bMYO7cuUG9ttI0N910EyUlJTz77LOqakdp8ww6AzGmGGKIqfZxl9elbfvy2CnzlGF3+//02CN6y1ek0gkdVoOVKEMUVqP/T//3Rr3aNqUotRFCYOnVi9KffkaYzSRedFG4Q1KUGomoI5U7oG2VClYPHG9BQfm12z/yCPgq/z4OND6WDm096fHgyc4Brxfv4cOV+vPUhyfXn9xJrZzcMXXW+u40dVtWa9GQj2D2UrkKZzOwxd8oTVHqNHXqVCZPnozJZOKvv/4iKyuL9PSqwxDalmXLlpGZmck999wT9FL2UaNGMWDAAObPn88zzzxDXNzRWyaU5rdhwwY+/vhjRo4cyTnnnBPucBQl4pn0Jkx6E/HmyqXYPunD7rFT6i6lxF1CiauEYlcxbp87TJFGHqPOSKwpVkueGbWbxWBRk88UpQnMvXtTtm4dCRdd1OA3qIrSnHRWrXqmPLljDWbljn/8eWL125iE0QhGIz67AwBPbi54tdoPx7btxDQ0uZNTfXLHctyxCKMRo9qlANQ/uRMrpQze3LQgEkI8hNZ0raJsKWXbzhpEoPPPP5/4+HhWr15N//79WbVqFZMmTQp3WGH14osvYjKZuOmmm0Jy/XvvvZdLLrmEWbNm8eSTT4ZkDaVhpk2bhk6n45VXXgl3KIrSoumEjmhjNNHGaFI58mLP7rFT5Cyi0FlIobMQuyc4n1K2BFGGKOLN8SSYE4gzx2E1NLxppaIotYs6eQi2zz+n3TVXhzsURalVxZ474K/cCVZyp0LlTo3rm81Ip5bccWdmlt/v3LaNmOGnNmg9T04Ouqgo9DHRle6PPessegw+odGj0FubeiV3IjWxU8E2YGSF79WWsAhksViYPHky7737DpecYmTjJgejzG+SNPhf0KHt7Vm22Wx8//33nHLKKSGrYLrwwgvJyMjgzTff5JFHHsFUocO80vxeeOEF1q5dyw033KCmmClKiFgN2mSmtOg0ABweB3mOPHLLcrE5bUjCMmwnJASCeHM8KVEptLO0w2Kw1H2SoihNEjd2LLFnnIHQ68MdiqLUSlfNtixvcVGTr2vfuJHSNWsQJlP5yPPqCKulvHLHfcif3NHpcG7f3uA13Tk5R1XtgLZVUk2sO6K11OV6pJRZFW654Q5Iqd7Uc3pjdzg5uPlrdMLHz38ZYMUF8MfD4Q6t2b300ku4XC6uv/76kK2h0+m49dZbycvL49lnnw3ZOkrdsrKyuP/+++ncuTPPP/98uMNRlDbDYrDQMaYjA1MHMrTDUDLiM1p8RYvVYCUjPoOhHYYyMHUgHWM6qsSOojQjldhRWgJRtXLHakU2sXLH53KxZ8pllK5YiXXQoFrbSujMlgqVO4cAsA4ciKMRyR1PTm61yR2lstaS3MkQQhwUQuwWQiwQQmSEOyCleidbl9K7o4EFP2QzIPZ31hcNoMRlgE2z4NCycIfXrObNm0dKSgoXhbgZ34033khaWhpz5szB51PNR8Plsssuo7S0lLlz56oKKkUJE7PeTJe4LgxJH0L/lP6kWFNaTP8ZndCREpVC/5T+DEkfQpe4LrVOIVMURVHatqN67kRFNXkUuregANxu0u69l67vvlPH+kcqdzyZmeji44kaPAjXX38h3Q3rj+fJVcmd+mgZr2hqtxqYCpwNXAekAz8JIardeCeEmCaEWCeEWJebqwp8mpso3cPU0+DHbZDuWolX6lldeDJIF2y4O9zhNZtff/2VrVu3MmHCBPQh/vRHr9dzww03cODAAd59992QrqVU7+mnn+a7775j8uTJjBs3LtzhKEqbJ4QgyZJEn+Q+DOswjN6JvUm0JEZcokcndCRZkuid2JthHYbRp10fkixJQW/AryiKorQ+R/XciW76tCxvfj4AhrS0Oo8VFis+/7Qs96FMjO3bY+7dG+l249qzp95rSinx1LAtS6kssl7FNIKUcpmUcqGUcqOU8lvgXLSf68oajn9dSnmilPLEFNXhvvnFdOfyUz3oBCz9pZDjY7aw1nYSTq8ZihpeotdSzZ49GyEEt99+e7Osd9dddxEfH89TTz3VLOspR6xfv54HHniAHj168Oabb4Y7HEVRqjDqjLSPac+AlAGc0uEU+ib3pX10+7Btc7IarLSPbk/f5L6c0uEU+qf0p31Me4w6NaZcURRFqT9dtNZ8WJ+YCGjbtJraUNkTSO4kJda9vtmMdDgBcGdlacmdXr0AGrQ1y1dcjHQ41HS6emjIKPQWQUpZIoTYDPQMdyxKNY67k47ZPzCmn+S9lbDqnJX8WXI864qGcGq7teGOrlm43W6+/PJLBgwYQO/evZtlTavVyuWXX86cOXP44osvOPfcc5tl3bYuPz+f888/HyEEn3zyCRaL6omhKJHMoDOQbE0m2ZoMaJO3bE4b+Y58SlwlDbpWlKHmJpMVCSGIMcaQaEkkwZygeucoiqIoQRE7bizS68HYqROgbcuSTifS62103yhvvn9KVlJSnccKqxWvzQZo07KiBg/G3L07GAw4t22Hc86p15pHxqCr5E5dWl1yRwhhAY4Floc7FqUaHc6GqK5MPX0Pf58D23dnkpGwi18KTubkLsWt7y9kNd555x1sNhtXXXVVs647Y8YM3njjDR555BGV3GkGbrebsWPHcuDAAV577TX69+8f7pAURWmgwOSt9Oj0Bm+FOrnDyQ06XsrWM8VLURRFCT9DYiJJU6aUf6+L0ip5fHY7+piYRl3TW6BV7tQnuaOzmPFkO/CWlOKz2TC0T0eYTJgzMnBu21bvNQPJHaPallWnFr8tSwjxHyHECCFEdyHEycAiIBpQzUUi1ZCXOX9YKvFR8M4KODXxR0q8MfwuJ4Q7smbx5ptvEh0dzbRp05p13eTkZCZNmsSaNWv45ZdfmnXttsbn83HhhRfy66+/cvvtt3PdddeFOyRFUZpIShnSm6IoiqKEks4/trwpTZU9BQWg06GPj6/zWK3njgNPljYG3ZjeHgBzr1713paV9+ZbFCz4CED13KmHFp/cAToB84FtwCeAExgqpdwb1qiUmnU4G2u/fzJ5mJ5Fa6CdPpsO1lx+2lyE78D/hTu6kNq9ezfr1q1j3LhxYdmi88gjj2AwGHjggQeafe3m5HK52LRpE+vXr2fnzp3NOiUskNhZunQpEydO5Omnn262tRVFURRFURSlOroorcGyr6y00dfw5hegT0hA6OpOI+gsFqTdjreoCDjS+8fcuxeezMzyLVs1kS4XOc88Q/HXXwOonjv10OKTO1LKyVLKDlJKk5Syo5RykpTyz3DHpdQh55/duLMAACAASURBVAeuPKsTdhcs3pjMqR33ku+KY+uqd8IdWUjNnj0br9fLrbfeGpb1u3fvzrhx41i+fDlbt24NSwyhsmnTJq644gq6du2K1WqlX79+DB48mGOOOQaLxcIxxxzDZZddxpIlS0KW7LHb7YwbN44lS5Zw3nnnsWjRopCsoyiKoiiKoigNEajckU2YmOXNz0dfj2bKAMJiwed04ivVkkm6aG19i7/nqHPHjlrPd2dmgs9H3PizSbnllvL4lZq1+OSO0kKV7Gbo8Un06mjmna+zOdawgiRjHj/uTkQuTofvRsOhZeGOMqi8Xi+LFy8mIyOD4cOHhy2ORx55BJ/Pxz333BO2GIJp69atjBgxgv79+/P+++9jsVi45JJLuP/++5k5cyZ33HEH48ePR0rJ/PnzueCCC0hJSWHKlClB3Z62Y8cOTjjhBL799lumTJnCkiVL0NXjUw1FURRFURRFCbXybVlNmJjlKcjHkFC/5E6gcsdXog0kCPT5MfuTO446+u649u0HIPHSS0m+4frGhtymtIX+tUokiumOsGcydXQ0987LZ9dBN6cmruLznPPYXRBHhm4HrJ0OJ83RmjC3AkuWLCEzM5MHH3wwrHEMGjSI0aNH88UXX7B582b69OkT1ngay+fzcf/99/PMM8/g8/mYNGkSM2bMoG/fvjWek52dzdy5c1m4cCHz58/nww8/pHv37lx00UXcdNNNdOzYscFxeL1ennrqKR599FHcbjePPvoo9913X1N+NEVRFEVRFEUJKp01sC2r8ckdb34B5mOOqdexwmJGut14i4q19f2j2Q2pqejj43Fs2VLr+e4DWnLH2Llzo+Nta9THykp4HHcneF1cPrQMIeC9H730j/2dWH0R3+eNRLpsoDfBltbTr+TVV1/FYDBw0003hTsUZs+ejZSS2267LdyhNEp+fj4jR47k8ccfp3fv3qxZs4aPP/641sQOQFpaGvfffz8bN25kx44dTJ8+HafTyVNPPUXXrl0ZNmwYr732Gg6Ho84YHA4HL7zwAhkZGdx7772kpqayYsUKldhRFEVRFEVRIo4or9xpwrasgoJ6b8vSWbRkkifvsPa9P7kjhMDSrx+OjX/Uer5r336E2ax67TSASu4o4dHhbDhpDp0SvYzpC++uAB1eRid/x0FHR34v6AFeN5TsCXekQXHo0CF+/PFHRo4cSbt27cIdDv379+ecc87h22+/Zc2aNeEOp0F2797NoEGD+PHHH7nuuutYv349gwYNavB1evTowYsvvsj+/fv56quvOPfcc/njjz+44YYbSEhIYMCAAVxzzTU88cQTLFiwgMWLFzN37lzuvPNORowYQUpKCjfffDPFxcXMnDmT7du3M3To0BD8xIqiKIqiKIrSNOWj0BtZuSO9XryFhRjqMQYdQFi14THew3na+v7kDoB1wACcf/2Ft6Tm5s7uA/sxdupUr+bNikZty1LCp8PZoDcx9XQ3l74EP2yBUcdv5FfbiXxzeAw9o14lOrl7uKMMihdffBGn08nNN98c7lDKPfvssyxbtozbb7+dlStXhjucetm4cSNnnnkmhYWFvPLKK1x/fdP33+p0OsaOHcvYsWNxOBy8/fbbLF26lA0bNvDWW29Ve05MTAwDBgzg8ssv56qrrsJkMjU5DkVRFEVRFEUJFZ0/2eJzNK5yx2uzgZToE+uX3NGZtfU8eXkIiwVhOJJ6sA7oDz4fjk2biB56crXnu/btx6S2ZDWISu4o4eVxMeFESIyG17+H0X0k56Z+ztz90/g0+zwuTfoDEe4Ym8jj8fDxxx+TmprK+PHjwx1OuYyMDC666CI+/PBDli1bxtlnR3Zvo82bNzNy5Ejsdjvz589n0qRJQV/DYrFw4403cuONNwJQUFDA5s2b2b17N6AldY477jh69eqlmiUriqIoiqIoLYbOoiVbpL3u9gPV8ebnA9R/W5Y1kNw5jM7fTDnA0q8fAPaNG6tN7kgpce3fT9TJQxoVa1ulkjtKeOn0WE1urjwNXvoGsm2QFp/D2OSvWJZ7Dv87UMrIcMfYRF999RU7d+7klltuibiEwHPPPcfSpUv517/+xfbt2yMuvoA9e/ZwxhlnUFZWxmeffca4ceOaZd3ExESGDx8e1ulmiqIoiqIoitJUwp/c8Tkbl9zx+JM79d6W5e+54z2cVz4GPcCQmIipWzfsv/9e7bne/HxkWRmmTqpypyEi852c0nbE9gLg+jPA7YW3/6fdfVL8WgbEbuB/2QP5/vvv8fl8YQyy8aSU/9/encdJUV77H/+cGWZjGXaEAQHFXUD4SRS3qASNoChq3ILeaBY0iolL0EvilhijRiR6o9ELyXXf16sEQSUacbmiUQRUUAzKjqDszMbM+f1R1dA03UPPMD01PfN9v1716ulaTxVDT/Wp8zwPd911F2bG5ZdfHnU4O+jcuTNXXnklX3zxBbfcckvU4SS1atUqjj76aL755hsefPDBBkvsSONjZuPM7D0zW29mq8zsRTOruRdtEREREcFyc7G8vF2o3FkDQG77dDtULgCCZlm5rVrvsLzooP6UfvQR7r7DsopFiwDI66nkTm0ouSPRGhgkFPYrgWP2h//+B1RXgxmcvNsLDOj0FTNmzOBvf/sbH330Ed9++y1lZWVJPwQao3nz5jFjxgyOPPJIevbsGXU4SV177bX06tWLW2+9lW/DjHxjsX79eo466iiWLFnC3XffzZlnnhl1SBKtY4C/AIcDQ4AtwKtmlt4jJBEREZFmzAoL61y5U7U2ltypXeVO9YYN23WmHFN4YF+qVq+mavXq7eavnjSJxRdeBEBB7951irW5UrMsiVbJMCAfqOCi78HZd8HLc+CEgyDHqjm53f3s0WIAr687heeff37rZmZGYWEh7du3p0+fPgwaNIji4uLITiOVu+++m02bNnHVVVdFHUpKubm53HHHHZx66qn89Kc/5dlnn406JCAYavzYY49l/vz53HLLLYwePTrqkCRi7r5d2ZaZnQesA44AXowkKBEREZEskVNYWPfKnbXrAMht1zbNYxVs+7n1jpU7BXvvBUD5F19sN9z5+hdeoEXnznS55WbyldypFSV3JHr9fgNzfsup36mmczHcOz1I7gCYOf3bfEjf1rNYkX8UKzv9hNKCvShb9TGlK99j1bpc3pyxlLfffJ1ju37I4Xtuwg4YGyaNorVixQomT55Mt27dGlVHysmMHDmS4cOH89xzz/HCCy9w8sknRxpPVVUVJ5xwAh988AFXXXUVV199daTxSKPVhqACdU3UgYiIiIg0dlZURHVZ3ZI71Rs3YPn55BQU7HxltlXuAEkrd/L79AGgfMEXtBo8eOv8LWvX0uaYY2kzZEid4mzOlNyR6PW7DoD8eX/iJ0ev5Y+TYck30KPjtlVyzCmpfp+S9Utgzx/Bqgeg+GtovZG1le2Ytur7vLp8ECs2zuWUr8+gRb+rtu43Kg899BBfffUVN910U6PtqDjeAw88wD777MOFF17IkCFDaJ0kw94QqqurGTlyJP/85z/52c9+xq233hpJHJIV7gRmAe8kW2hmo4HRQKNtFikiIiLSUHIKCvA6Jneq1m8gp02b9I8VjpYFkNN6x+ROi86dyWnblvIFn2+d5+5UrVmbdr8+sr3G/41Tmod+18EZaxg9NB8HJr6WZJ2qzbDp3zD3t1C6DKo3AtAuby1ndnuCoR1fYe6Gvjy17BSqPr4Vlr3UoKcQb+nSpTz11FPk5+dz4YUXRhZHbXTq1Inbb7+dFStW8JOf/CSyOH70ox8xefJkTj/9dO69997I4pDGzcwmAEcCp7t7VbJ13H2iuw9y90Gd48p9RURERJqjXa3cya1Fcic2Ohckr9wxMwr69KFiwRfbjrFpE2zZQm67dnWKsblTckcalT367MdJA+GeV6G0IsVKXg3V238omcERHd7ixM6T+WzTPrywbBj+yW2ZDziFKVOm8NFHH3HmmWfSsWPHnW/QSFxwwQUcf/zxPPnkk9x3330NfvxLLrmEhx9+mKFDh/Lkk09mRcWTNDwz+xNwDjDE3f8ddTwiIrIjM+tgZn82s3lmVmpmi83sHjOr8cbIzM43M08yFda0nYjsXNDnTmmdtq3asLF2lTs7Se4AFOy1F+VfbEvuVK2p3Yhcsj19c5LGZeAtXDEMVm+Ah9+s/eaD2r3PsR3/wewN/Zg2v3Mko2p99tlnPPfcc1RUVDBmzJgGP/6ueuqppygpKeHSSy9l3rx5DXbcX//61/zlL3/h8MMPZ8qUKUrsSFJmdifwQ4LETsP9goqISG2VAN2Bq4B+wLnAd4HH0th2M9AtfnL3upUbiMhWVlhAdXl5nbatXr+e3Dbpd9tgcX3z5Kbo7qFgrz5UrVnDlm++AaBq7dpg/TQ7bZbt6duTNC4lwzj6gBYM7A1/eikYFr22jmr/Boe2+z/eXX0Ab9z/Y/jfPWD6kAZpprVlyxamTJnCzJkzGTx4MIceemjGj1nfiouLeeaZZ6isrGT48OGsW7cu48f83e9+x80338zAgQOZPn06eXl5GT+mZB8zuxu4gKBqZ42ZdQ2naDqIEhGRlNx9rruf5u4vuPsCd/8nMBYYamY7G+LU3X1F/NQAIYs0eTmFRXhZHSt3Nm4kp036oxNbTs7WBE9NlTsQdKoM2yp3Wqhyp06U3JFGx9oewBXDjE+XwbTZddje4PudpjGg+CNeX9Sbd2d0hbNfg+7DoUsLuCtzw5K/8cYbvP3223zzzTeMHTs2Y8fJtMGDBzNhwgQWLlzIcccdR2VlZcaONXbsWK6//nr69u3L66+/TmGhqq4lpYsJRsiaDiyPm34VZVAiIpK2YqCcoDKnJkVm9pWZLTGzyWY2sAFiE2nyrLCA6rI6Vu5s2EBOLSp3YFvTrJxWybfL7xMbDn0BEF+5oz536kLJHWl8Bt7CmYflUtI+qN6pCzNnRJcX2a/sU6a2OYE39zsSB1hVBVfeBr/qWe+VPIsXL2bGjBnMmjWLPn36cMopp9Tr/hvaJZdcwuWXX857773HySefTFVV0j5r66y6uprRo0czfvx4Bg0axDvvvENxcfpPA6T5cXdLMd0QdWwiIlIzM2sH3AhMcvctNaw6H/gxcApBpWYZ8JaZ7V3Dvkeb2ftm9v6qVavqM2yRJiWnsIjqulbubNhAbuv0+9yBbZ0qp6rcadGlM1ZYSOXiJcEx1OfOLlFyRxqfkmHkt+zApSfk88pcmL2obrvJsS2c/j9P02/2bKYPHcrTZ5zBplatoAJ4YDG8PhweNXh+j11O9GzcuJGnn36ab7/9ls8//5wrrriC3NzcXdpnYzBhwgTOPvtspk6dyrBhw+qtgmfTpk0cf/zxTJo0iaOOOooZM2ZENvS6iIiIpM/Mfp+iw+P46ZiEbVoBLwJLCfrgScnd33H3B9x9lrvPAM4CvgAurWEbjY4okoacokK8tPbdV3llJV5aSk5x7ZI7Wyt3kgyFDsGIWXndu1O5JEjubFm7FnJyatVxs2yj5I40Tu0OZPSIXrQuyuEP/1v33bT4uopTn32W7736KvP33Zc/X3oprwwdysrc3dja1/LmL+HNc2qf4JnzO3iqPRWPFPHY3WMp3bSe+fPn07FjR84///y6B93IPPLII4waNYpXXnmFwYMHs6tPxGbPns2AAQOYPn065513Hq+99pqaYomIiGSPO4D9dzLNjK0c9osWu8k6qbYdI7t7FfA+kLJyR0TSYwWFte5QecuaNVRt3AhQ75U7AHk9ulOxdFvlTm67dpgGVqkTXTVpnPYfS4eWVVx6UnuefBc+WVLH/XQCA458800uvPde9lqwgHcOP5x7f/5zJiy8kkeWjuKV1cfx0be9+fr9u6lOtwfnOb+DuTeyqbSKBxefy/KyLvStvJ9p06Zx8cUX07JlyzoG3Pjk5OTw8MMPc8UVVzBr1iz69evH3//+91rvp6qqiuuvv55DDjmExYsXM378eB588MEmUeEkIiLSXLj7aneft5NpM4CZtQGmArnAcHffWNvjmZkB/Qn6WBORXZBTVAiVlXia1fibP/iQz484ktKPPgq2r2VFTaxyJ9VoWQD53XtQuWQpAFVr16m/nV3QIuoARJIqGQbfuYsrqm7iv154i98/D4/WZVTxM4G/AhXQefVqfvD002xo35oFF+zFl933YGX5biws3YMqbwErocWs69gtfyUlRSvp0aUlPQaeRft9RxLcV2zjn/6JuesPYNqq4ymrLuTMbk9y7d++pE0hXHZW3/q4Ao3O7bffzsEHH8zFF1/MiBEjOPHEE7n99tvZZ599atyuurqahx56iOuvv56vvvqKfffdlyeffJL+/fs3UOQiIiLS0MLEzssEnSiPBFqFzbMAvnX3inC96cBMdx8Xvr8e+D/g83DbXxAkd37esGcg0vRYYREA1eXl5KYxOm3Zxx9DdTWlH84CILeWzbKsKDhezZU7Pahev56q9eu3Vu5I3Si5I41XyTA6nTaMS44zbvs7XHca7FdSy30cEb4+CawGOkGbMzcy8OBZDCT4kKr2HFZXdGRFeTeWl3dlRXk3Zq3rx3tr8+Gz2bQs/IQePfekc+fO5OTksH7FJyxceAHrt7Rlt/wVnNf9Qb5e+TVPz4RrRkKHFffCgWfW44VoPH74wx/yve99jwsuuIApU6YwdepUDjnkEE477TSGDBlC//79yc3NZenSpcycOZMXXniBqVOnsmLFCtq2bcuNN97IuHHjVK0jIiLS9B0MDA5//ixh2bHA6+HPfYDFccvaAROBrsA64EPgu+4+ExHZJTmFwdDkXlYGafR3Wbkk+K9ZPn9+sH0tm2Xl7GQodAiaZQXHWkLV2rXk9ehRq2PINkruSKN35XC46xW46Xl46OI67OAItiV5ksixaroUrKJLwSr6E4y9Xu3G1xVdWMIglpbtzuJv2rJgwQKqq6tpk7eJkoLlHN9pGge0/hQzZ8xz0KYQLh9msPHLupxm1thtt92YMmUKH374Iddddx2vvfYab7/99tblZoZv7dAI+vTpw4033shll12mTpNFRESaCXd/naB1/M7W653w/nLg8sxEJdK8ba3cKUuv66uKRWFy57MgP5tby6HQragIKyzEWqROO+SHyZyKJUuoWrOGwr4H1uoYso2SO9LodWmXx8VDK7l9Cow9Cfr3zPwxc8zpWrCSrrzEoNbFcMSj8Olt+NfvYOZQva0jsllfsq1qp3UOtO6d+QAbgYEDB/Liiy9SUVHBtGnTmDlzJosXL6ayspIOHTqw3377MWLECHr2bIB/MBERERERqdHWyp3S9IZDj1XuVC5bFmxfXFzL4xXWWLUDbK3UqVyylKq1a2mhYdDrTMkdafyK92fcyXP562vVXPUYTL26IQ9eDZVrg2HTCR8/bStKwR1+9Sh0aA1XDDcwg/3HNmSAkcvPz2fEiBGMGDEi6lBERERERCSF2OhV1WU7HzHL3alYvP2oNjV1jJxM+3POpuXgQ2tcJ7e4mJziYso//xyvqFCfO7tAo2VJ4zfwFjp07Mw1I41ps+GVObEF0ffb8tJHMP1juO5UaF9cBMUHBJ1Bi4iIiIiINCI5YQfHXrbzyp0tX68K+uaJ63g5p5bJnaKDDqLdyJE7XS+ve3fK5gZf8nJVuVNnSu5I41cyDA67jzFnHkjvzvCrR42qwj0hJ9rCsy1VMPZR2Gs3+PlQg8KuMPCWSGMSERERERFJxsIOjtOp3Ik1ySoKR7jNadmyxr5zdkV+jx6Uf74AgNwOHTJyjOZAyR3JDiXDKDhlDrf89mpmL3ImTl0F1RWRhjTpNfhkKfzxHMgvagvfuUtVOyIiIiIi0ijFKneq06jciXWm3OrQQ4Jt29RupKza6DTmEjqNGcNu11xDq8MPz9hxmjr1uSNZ5cyLbmbSM+/zn4+8zsiBeXRrG02CZ/kaGPcEDDkQRn4nN+hwWYkdERERERFppHLCPnc8jdGyKhcvgpwcig4+GIDc4swldwr33ZfCfffN2P6bC1XuSFYxM+655x7KK6v45QNbIovjsoegrBLu/TFYuwOV2BERERERkUZtW4fKO0/ulC/4gryuXcnv1QuAnNaZS+5I/VByR7LO3nvvzTWj9uCpd6t58QNr8ONP/gCefBeuOQX27l4IA9TPjoiIiIiING5bK3dKa07urJ44iQ0vv0yro44ir0sXMCOnTe06U5aGp+SOZKWrrr+DfrsbP53krFzXcMf9eh389K/QtwdcNQI4YJyqdkREREREpNFLp3Kn/N8LWTVhAsXDh9H1mt9g+fm06NaVFhrFqtFTckeyUn7vk3n0ym6sK4UfTwT3zB/THS6YCGs3w2NjIL9FDvS7LvMHFhERERER2UWx0bLi+9zZsmoV3z76KB5+oSqdNQuATpdcgoXDoPe48046XfqLBo5WakvJHclafXvnM/6HMGUWTJiS+ePdMTU41m3nQN/dgV7nZP6gIiIiIiIi9cDMsKKi7Sp3Vv35Llb+7kYqly4FoGzuHHJatSJ/jz22rlPUrx/5Pbo3eLxSO0ruSPZqvQeXHAenfQeuegxempW5Q02bDb96BE4dBGOOB3qNgiMeztwBRURERERE6llOQQEeDoVetXEj6yZPBqBy2TIASmfPobBvXyxHqYJso38xyV77j8Ushwcvgv494ey74OMl6WxoQIu0D/PpUjjrz0G1zoM/N+zYKUrsiIiIiIhI1gkqd8oBWD95Mr55MxAkd6orKiibP5+i/v2iDFHqSMkdyV4lw6DXObQqhBeuhJYFMPRm+Gx5TRsZ4EB6w6gvWBHss6AFvHCF0brQ1IGyiIiIiIhkpZzCwq2VO2uffY78vfoAsGX5csrnzYPKSgr7KrmTjZTckex2xMPQ77fs3rmQ6eOgqhqG3ASfpKzgSb/n5XnLYMgfoLwSpl9bSK8ueZBXXC9hi4iIiIiINDQrLKS6rJzqsjLKPv6YNkOHktupE5XLllE6Zw6AKneyVJNJ7pjZxWa20MzKzOxfZnZU1DFJA+l3HZxdygHnTuHV35ZQ5XDYDfDKnLrv8rWP4bDroawSXh2XQ9/u1UA17Hd5fUUtIiIiIiLSoGKVO+Xz50NVFYUHHkhet25ULl1G2ew55HbqRIuuXaMOU+qgSSR3zOws4E7gD8BA4G3gJTPrGWlg0rBKhtH/sqW8+/L99Oycx/dvhasfCypvIG/rao+8dQ69f7mQnFFV9P7lQh55a9uoVxVb4DdPBk2xurWDd2/MZ8AeOdCiJfS9VkOfi4iIiIhI1rLCAqpLyyidOxeAor59ySspoXL5ckrnzKGoXz/MLOIopS7S71W2cbsCuN/dJ4XvLzWzE4CfA+OiC0ui0PM7P+Kdj0/nirP24Y+Tl/P0+/nccG4JZ/X7kqfePYfRf53E5opWAHy1ujej/zqJLVXVtMh9ghuegQUr4cfH5PKnyw+n+OQ3Ij4bERERERGR+pHTqhXlixZTNmcuuR070qJrV/JKStj4j3/gW7ZQfNKJUYcodZT1yR0zywcOBsYnLHoZOLzhI5LGoHXr1kyc9DfOmPRjrn5wDf8x/ksua22UV1axueIRoISgc+VlbK6YyY8nvkq1B6NuTbkql2GHdIJByguKNEZmdjEwFugGfAxc5u4zoo1KREREpPFrM3QoG1+dzvrVq2k5+FDMjLxu3fDKSgCK+qm/nWzVFJpldQJygZUJ81cCOzQWNLPRZva+mb2/atWqhohPolIyjON+9j+8f+9hTLl2N044pAubyl8GLgRGACcBo4GnqfajeeFK48ObCxj23QPhsPs0KpZII6RmuCIiIiJ113b4cFp06YKXl1N0YF8A8rqXbF1e2LdvVKHJLsr6yp04icMgWZJ5uPtEYCLAoEGD0h86SbJTyTBySoYx7DgYBrzZZTmLVpUDXwPVwG5AL3p1+ZoR45+JNFQRSYua4YqIiIjUkeXn0/68c1l1+wQK+x4IQF63bsHr7rvTon37KMOTXdAUkjurgSp2rNLpwo7VPNLM/eG6lYz+1T5sLu+9dV7Lgs3cdO0KkhR6iUgjoma4IiIiIruuw7nnklPUktZHBQNM55UElTtqkpXdsr5ZlrtXAP8CjktYdBxBub7IVqPGDGDi+M/o1WUFZtX06rKCieM/Y9SYAVGHJiI7p2a4IiIiIrsop6iIDueOwvKCEYVz2ral+KSTaDvylIgjk13RFCp3ACYAD5nZTOAt4CKCHnPvjTQqaZRGjRnAqDGxd11RxY5I1lEzXBEREZF6YmZ0H39b1GHILmoSyR13f8LMOgLXEIyeMhcY7u5fRRuZiIjUIzXDFRERERFJIuubZcW4+1/cvbe7F7j7we7+RtQxiYhI/VEzXBGR7GJmr5uZJ0yPp7Hd6Wb2iZmVh6+nNkS8IiLZrMkkd0REpFmYAJxvZj81s/3N7E7UDFdEpDG7j6CyPjZdWNPKZnYY8ATwCDAgfH3KzA7NcJwiIlmtSTTLEhGR5kHNcEVEss5md19Ri/UvA15z95vC9zeZ2bHh/HPqPToRkSZClTsiIpJV1AxXRCSrnG1mq83sYzMbb2ZtdrL+YcDLCfOmAYdnJjwRkaZBlTsiIiIiIpIJjwJfAcuAA4GbgYPYse+0eF3ZsZP8ldQwvKmZjQZGA/Ts2XMXwhURyV6q3BERERERkbSY2e+TdJKcOB0D4O4T3X2au89x98eBs4ChZvb/dnIYTzxsknnbVg6OM8jdB3Xu3HlXTk9EJGuZe8rPySbPzFYRPE2orU4EQ/I2ZTrHpkHn2DTU9Rx7ubvuctmlz3sRkWzQYJ/3ZtaJ4O9STRa5++Yk2+YAFcAod38ixf4XAX9299vi5o0Fxrh7rzTi0/19/dO1SU3Xpma6PqnV+/19s26WVdc/gmb2vrsPqu94GhOdY9Ogc2wamsM5ZpqSXCIi9cPdV1P3L2v9gFxgeQ3rvEPQbOu2uHnHAW+nGZ/u7+uZShe2pwAAErpJREFUrk1qujY10/VJLRPXplknd0REREREpP6ZWR9gFDCFIBl0AHA78CHwVtx604GZ7j4unHUn8IaZjQOeA04FjgWObLjoRUSyj/rcERERERGR+lYBfI9gpKv5wH8RjII11N2r4tbrA3SLvXH3t4GzgR8Bs4H/AM5y93cbKG4Rkayk5E7dTIw6gAagc2wadI5NQ3M4R5GkzGxy2EHrr6KOJV1xMV8Z4bGz5npJ0+Tui939aHfv6O4F7r6Xu//S3b9NWK+3u5+fMO9pd9/P3fPdfX93f7YBQtbf2tR0bVLTtamZrk9q9X5tmnWHyiIiItK4hZ2r7g58391fjjqedJjZEqA7cJy7v1qP+/0B0Bd4zd3/mWKdrLteIiIisuvU546IiIg0SmbWCthE0KTjo4jDSYuZdSRI7ADMqufdXwv0B+alOHbWXS8RERGpH6rcEREREaknZjYUeAVY6u496nG/RcB6ggdze7v7gvrat4iIiGQ/9blTC2Z2sZktNLMyM/uXmR0VdUz1yczGmdl7ZrbezFaZ2Ytm1jfquDLFzH4d9ktwV9Sx1Dcz62ZmD4T/jmVm9omZHR11XPXFzHLN7Ma4/48Lzez3Zpa11Yhm9l0ze8HMloa/l+cnLDczu8HMlplZqZm9bmYHRhSuiKR2UPha31U7AwkSO+uAL+p53yLNVlO/v09HeH/hCdOKuOXN6h6kPu7JzKzAzP5sZqvNbFO4v3pL+EcljWtzf5Lfpf9LWKepXpudfpfO9O+OkjtpMrOzCIZm/APBDdbbwEtm1jPSwOrXMcBfgMOBIcAW4FUz6xBlUJlgZoOBnxGMwtCkmFk7giFGDTgR2B+4FPg6yrjq2dXAJcAvgP2AX4bvx9W0USPXGphLcC6lSZZfBVxJ8G/5HYJ/z1fMrE2DRShSD8xsqJk9YWaLzKw8vAn63MyeNrOfJqz7fHhj+J9J9hNbNs7MWpnZdWb2sZltNrPFZjbBgmZKsfWPMLNnzGx5+CVuppkNSbJfC2NyMzs0xTnsHi6vTvJ/cED4ukOzKDMrMrNTzOye8EvkCjOrMLOV4fkckWSbwWbmbBs6ui1QnXDj3H5n1ytuf8eZ2ZNm9mV4/Vea2WOpvqhZXAfNZpZnZpea2bvhNVphZs+a2Z6pjifSmDWT+/t0zScYtSw29Ytb1tzuQerjnuwO4HTgHOAooBiYbGa5GYy7Iezs2gC8yva/S8MTljfVa3MMO/8undnfHXfXlMYEvAtMSpj3OXBz1LFl8JxbA1XAiKhjqefzakvw1HMI8DpwV9Qx1fP5/QF4K+o4MnyOk4EHEuY9AEyOOrZ6Or+NwPlx7w1YDvwmbl4RsAG4MOp4NWlKZwp/j/8CeNy0jmC45Nj7+Qnb/DucPyzJ/mLLfgEsCH9eH95Mxfb3eLjub4FqoDJcJ7a8FOiZsN+9wmVVQMsU53JyuM4XSZbNCZedkWTZBQnnvx4oi3u/w99c4IfAijBWDz8fVsRNH6V5vYqABxOOFX8tNgPfS7LdonD5TwgeiDhBvz6b47ZdDLSO+ndMk6baTjTD+/sU1+EGYG6KZc36HqQu92QE3zUqgFFx6+we/h36ftTnlKlrE867nxrux5vLtQnPa7vv0g3xu6PKnTSYWT5wMJA46sTLBJm5pqoNQXXXmqgDqWcTgafd/R9RB5IhI4F3wyfjX5vZLDMbY2YWdWD16E3gWDPbD8DMDiBI1k2JNKrM2QPoStxnkLuXAm/QtD+DpGn5GfBz4FvgR0Abd28LFAI9CZ5iPRZb2cyKgd7h2+2qYBKW/Y4gATHA3YsJnnA9HC47w8xuBcaG+28XrnMEQWKjMIwpXqxZ1WfuvjnFuQwMX7er/jSzAoJqQkjeLKstcCswGGjr7sXuXhge8x8Ef3PHx2/g7o+6e1eCL6EAV7p717jpoCTXJPF65QCPAucBKwmeBrYMr8WewHSCG8xHE6qd2hPcVBLGvYHgM6c10Ar4AUEyrQdwVtIrJdJINeP7+1T2tKCpzUIzezyuIk/3INtL53ocDOQlrLMY+JTmcc2ODL+DfGZmk8ysS9yy5nRtEr9LZ/x3J2v7p2hgnYBcghuieCuBoQ0fToO5k+Dm9J2oA6kvZvYzgqey50UdSwbtCVwM/Am4haCJwJ/DZU2lf6FbCT4wPzGzKoLPspvc/S/RhpUxXcPXZJ9B3RHJDj8JXy929ydiM929mqDyI/HzqT/BU67V7r4sxTIIKmWOd/ct4f42m9kvgHMJbqrGAsPdfWrcMd82s/sIysr3Sdh3rFlVTX3m/L/wNbHpVV+Cz6ONBNVE23H3O5LtzN1nm9mZwCpgHzPr4O7fJqzWP8Ux45enul5XEiT+vwEGu/uXccdeaGZnE1T9dCGoSool2frH7WMBcKy7V8TNe8bMJof7TryOIo1dc72/T+Zd4HyCkfi6ANcAb4fNNXUPsr10rkdXgoqN1UnW6UrTNhV4FlhI8MDh98A/zOxgdy+neV2bxO/SGf/dUXKndhKHFrMk85oEM5sAHAkc6e5VUcdTH8xsX4ImS0cl3Jw2NTnA++4e63/mQzPbm6BPmqaS3DkL+A+C5gofE3wZu9PMFrr73yKNLLOazWeQNEmxPiwq01w/VkGTLJkRW1YJnBtL7MS4+xozKyWoRpkYn9iJE0uA5CXMTye5k7RyJy6u2R7WUtfCJoKy61wgP36BmfUC2ofL56TYPun1CitxrgnfXhmf2Ilx99Vm9gFwNEFCJ5bcie1zC/DDFH87Y08k9Vkk2arZ/21195fi31vQAe6/CaosY53hNvvrlKAu16PJXzN3fzzu7Rwz+xfwFUE/oM/WsGmTujY7+S6dsd8dNctKz2qCDFpitqwLO2besp6Z/YmgZHuIu/876njq0WEET2nmmtkWM9tCcCN7cfi+INrw6s1y4JOEeZ+y7YtVU3AbMN7dH3f3Oe7+EDCB7O5QuSaxESuaxWeQNFlzw9f7zOxqM9vZE990kjuvuvtXiQvNrCNBYgcgVcK3V/i6KMW+kyZ3wn3HmiolxpayM+Vw287hub8WdkZcHusUmaBPnVyChFXiE7tYTAvcfVOK80l1vX5A0FRtFfA4qcU+S+KTXbHKnddquB+INd34soZ9izRGzer+vjbcfSPBw7O90T1IonSuxwqCz/NONazTLISVpEsIfpegGVybGr5LZ/x3R8mdNIRPqv4FHJew6DiCXvWbDDO7k6AaYoi7z4s6nnr2PEHP/wPipvcJbnYHEHRe1RS8BeybMG8fgqx5U9GS4IYsXhVN9zNtIcGH/dbPIDMrJOhBv0l9BkmTNoZgNJZigiajS8zsAzP7jZklKzVOJ7kzOcWxYkmWr939vRTrxEaCiSWdCEe0iCVuUlXuxJpkbSR4sp3suDtsa2anAZ8RnPsxwG4EnSmvDKcN4aqfJ1YiUfO12Nk63wtfnw5L4lNpG77GJ5Zi+3yxhu12uI4i2aA53d/XVniPsR/BA0Pdg2wvnevxL4JEffw6PQhGsG1W18zMOhE0OVoezmrS12Yn36Uz/rujZlnpmwA8ZGYzCb48XwSUAPdGGlU9MrO7CfqiGQmsibvZ3hhm8LOau68F1sbPM7NNwLfu3pRuSv9E0E76N8ATBM0HfgH8OtKo6teLwH+a2UKCJ0sDgSsIRoLJSmbWmqA/KAiSVD3NbADB7+ciM7sD+I2ZzSP4gngNwZfLRyMJWKSW3H2+mR1E0KfL6QRJh4Hh9Gsz+0GsaUDYAXDfcNPETovjl32Q4nCxpMSHNYQUq0qJT8TEtlvh7qmekA0KX+ckaXqVbJ+Y2WHAkwRP4x4H7iFoPrs5bp27CJrPJksq7awiKOX1ipu/s79zsaHQZ4f7zI2b968Ux+0JdCAoFa8p8STSWDX5+/t0mNl4gnurRQQVAtcSdJr+gLt7c7sH2dV7MndfZ2Z/A24zs68J+jubQPD5+mrDnk39qunahNMNwDMEyZzewM0Ew30/B03+2tT4XTqd/0u7fH0yPQRYU5oIOqn9EignuNH5btQx1fP5eYrphqhjy+A5v04TGwo9PK8TCW60y8IPjl8AFnVc9Xh+bYA7CKqRSgmenv8BKIw6tl04p2NS/P+7P1xuBH8wl4f/rv8E+kYdtyZNdZ0IEh2nhH9XnWB0qtiyfcJ5FUB+wnaxZTUNVR4b8vuWFMt7x+2jVdz8X4bzX6oh7nfCde5JmL9nOH8LUJSw7KVw2aQU+zRgabjO2CTLY0O9n5hi+5qu18Jw2Xk1nFP/cJ1yghHFIHhqH7tGSYc5p4Yh4TVpypapqd/fp3kNHifoh6wi/Cx6Bjggbnmzugepj3sygtEY/0zw5XwzQfJs96jPLZPXhqA59DSCZE4FwX36/Ynn3YSvzU6/S2f6d8fCHYiIiIg0ODMbDvwdKHX3luG8MwgqXea4e/+E9WPL5rv7fon7C9eZTdBc6BzfvnPH2PKRBE8RP3P3fePmx6pn7nL3S5Ns15+gssYIRv26J27ZqQSdRX7q7gckbLeeICk91N2nJ9nvD4CnwrfHufurcctaETTZMqCnB0OiJm5f0/WaQ1C9c5W735a4bbjOIwRl5E+7+xnhvLMIvvBtd40StruWYCj6Z9399GTriIiISMNoqv1TiIiISCOQRmf1sWZD8R3Bp9PfTqoOj/MJqk5SrlPD/mPDq3dJst9WBKNIxdapTWfKsc6dOybZb0/gv+NmJcbcPe6YiUOcx9R0vWL7Oz7Zhmb2Q4LEzhbg+rhFSZuYJaixuZiIiIg0HCV3REREJJMuMrN3zOynZhbrrBgz62BmY4Hfh7MmxG1T5+QOQT8xeQSlzJ+lWCfVPuaHr6ea2RlmlhdOJxA0x+oQLnd2HJK8piHUY+dxg5ntD0HSK6yOeSduvaXunjhSVnzfP6mqY2q6XveFr0PN7AYzaxMev4uZ/R54IFx+kbsnS7DV1G9ROsPGi4iISANQckdEREQy6SBgMDAJWGRmm8xsHUFb8j8S9L1zrbs/mrAN1C25E0s4zHH36hpiSrb/hwiGDM8jaOZUStAm/qXw9ZpwvYXuviFh25oSHdcTJIT2Bz4Jz38DQbOnWWzrvHWHbd19HfBK+PYJM1sfDqO+MOxIuabzwd3/AfxXXBzrzGwtQdLoNwQVOz9398Qh42us3AmTRHukOq6IiIg0LCV3REREJJMmAFcDUwk6P3egIPz5PuA77h6r3sHM2rNtOPLEkbLil+2syVU6SYnt1nH3NQRDkk4mGBJ8E/AeQT88hxMMX54srnZAz1THdfe/AycBMwk6bd0CvAmMcvcTCTpjrumcRgF3E3SOXBDG8Y27V9d0veKO/0uCETzeJqhoagHMA+4EDnL37UYGqsV1NmCtu3+VYh0RERFpIOpQWUREREREREQki6lyR0REREREREQkiym5IyIiIiIiIiKSxVpEHYCIiIiIiIhINjGzHwB9gdfc/Z9RxyOi5I6IiIiIiIhI7VxLMLLgvKgDEQE1yxIRERERERFJm5kVAQeEb9+PMhaRGCV3RERERERERNI3kKAVzDrgi4hjEQGU3BHJODN73szczMaZWWszu87MPjazUjPbYmaFUccoIiIiIiI1M7PBZubAW+GstkB1eK8fm9pHGKI0Y+pzRyTz+oev5cAcoDdQClQDn7t7WURxiYiIiIhI+vYEVhIkdQqBTcDGuOUr3X1NFIGJmLtHHYNIk2VmxcBawAg+/L8ELgVed3c3s17u/lWEIYqIiIiISC2Y2evA0cBF7v7fEYcjAqhyRyTT+hMkdgCWA0e4+7rYQiV2RERERESyTqwy/6NIoxCJoz53RDLroPDVgR/HJ3ZERERERCS7mFkvoD1BFwtzIg5HZCsld0QyK5bcedfdZ0QaiYiIiIiI7KrY/f0Cd98UaSQicZTcEcms2If//0YahYiIiIiI1IfY/b2aZEmjouSOSIaYWQ7QN3z7RpSxiIiIiIhIvRgQviq5I42KkjsimbMX0JKgPe7siGMREREREZFdF6vcmRVpFCIJlNwRyZz49rgbI41ERERERER2iZm1AvYM3+rhrTQqSu6IZI6y+iIiIiIiTUd3wMKfl0UZiEgiJXdEMieW3Pkw0ihERERERKQ+rIz7+fTIohBJQskdkcxR5Y6IiIiISBPh7uuAV8K3T5jZejNbYWYLw8FURCKjX0CRDDCz9sDu4Vsld0REREREmoZRwN3AQqAA2A34xt2rI41Kmj1z96hjEBERERERERGROlLljoiIiIiIiIhIFlNyR0REREREREQkiym5IyIiIiIiIiKSxZTcERERERERERHJYkruiIiIiIiIiIhkMSV3RERERERERESymJI7IiIiIiIiIiJZTMkdEREREREREZEspuSOiIiIiIiIiEgWU3JHRERERERERCSLKbkjIiIiIiIiIpLF/j9CDAxPdZnCvgAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"#plot\n",
"simulation_analysis_plot(system4, title=\"position Langevin\", limits_coordinate_space=list(range(0,10)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the first panel of the visualization we can see how the system's energy minimum around x=4 was flooded and the particle can cross neighboring energy barriers more easily. The longer one simulates, the flatter the whole energy surface become. Note however, that artifacts can arise once the system leaves the grid defined by $bias\\_grid\\_min$ and $bias\\_grid\\_max$. Therefore, these parameters have to be selected very carefully.\n",
"\n",
"\n",
"The Ensembler package also contains a build-in functionality to make short movies of the simulation using `animation_trajectory`. An example code is showing the beginning of the trajectory is given below."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"