Skip to content

Instantly share code, notes, and snippets.

@dwhswenson
Created January 27, 2019 23:35
Show Gist options
  • Save dwhswenson/bb79a137a1d65629c22e7b00aa569d76 to your computer and use it in GitHub Desktop.
Save dwhswenson/bb79a137a1d65629c22e7b00aa569d76 to your computer and use it in GitHub Desktop.
Examples for using the Double Well Dimer testsystems in OpenMMTools
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Double-Well Examples\n",
"\n",
"The symmetric double-well \"bond\", where there are two metastable \"bond distances,\" is a widely-used example system in the study of rare events. This notebook will explain the `openmmtools.TestSystem` subclasses that have been implemented to facilate use of variants of the system.\n",
"\n",
"Note that this requires a few packages in addition to the requirements for OpenMMTools. All requirements are conda-installable (in channel `conda-forge`) and included in the import statements in the next cell."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from __future__ import print_function\n",
"import sys\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from simtk import unit as u\n",
"from simtk import openmm as mm\n",
"\n",
"import nglview as nv\n",
"import mdtraj as md\n",
"\n",
"from openmmtools import integrators, testsystems"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The system we'll describe here is a symmetric quartic double well dimer in a bath of WCA particles. The potential for the double well is:\n",
"$$\n",
"V_{dw}(r) = h \\left(1 - \\left(\\frac{r - r_0 - w}{w}\\right)^2\\right)^2\n",
"$$\n",
"where $r$ is the distance between the particles, and the parameters $r_0$, $w$, and $h$ define the minima and maximum of the quartic potential such that the two wells have minima at $(r_0, 0)$ and $(r_0 + 2w, 0)$, and the barrier has a maximum at $(r_0 + w, h)$. The simple bistable nature of this system, combined with the fact that the order parameter is obvious, makes it an attractive system when developing and validating new approaches for rare events.\n",
"\n",
"The potential for the WCA interaction is:\n",
"$$\n",
" V_\\text{WCA}(r) =\n",
" \\begin{cases}\n",
" 4 \\epsilon \\left( \\left( \\frac{\\sigma}{r} \\right)^{12} - \\left( \\frac{\\sigma}{r}\n",
" \\right)^6 \\right) + \\epsilon & \\text{if $r\\le 2^{1/6} \\sigma$} \\\\\n",
" 0 & \\text{if $r> 2^{1/6} \\sigma$}\n",
" \\end{cases}\n",
"$$\n",
"where $\\epsilon$ provides an energy scale and $\\sigma$ provides a length scale. This potential is repulsive at short range, and zero at longer range. Note that the WCA interaction applies to all pairs of atoms, including the dimer.\n",
"\n",
"This model is usually described in units of $\\epsilon$ and $\\sigma$, however, OpenMM still requires real units. In practice, we'll talk about energies in terms of $k_B T$, and distances in terms of $\\sigma$. In addition, there is a mass unit $m$, and a time unit $\\tau = \\sqrt{\\sigma^2 m / \\epsilon}$. \n",
"\n",
"In this first section, we'll set all the parameters for the systems we'll build. The actual parameters we use for $\\sigma$, $\\epsilon$, and $m$ are based on the Lennard-Jones parameters for argon, but after the initial definition, we prefer to define other quantities in terms of $\\sigma$, $\\tau$, and $k_B T$."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# energy in OpenMM is energy/mol, so k in kT is k_B*N_A\n",
"kB = u.BOLTZMANN_CONSTANT_kB * u.AVOGADRO_CONSTANT_NA\n",
"\n",
"# these are the parameters from the default WCA Fluid\n",
"epsilon = 120. * u.kelvin * kB\n",
"sigma = 3.4 * u.angstroms\n",
"mass = 39.9 * u.dalton\n",
"tau = np.sqrt(sigma**2 * mass / epsilon)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We'll define the temperature such that $k_B T = 0.824 \\epsilon$ (this is a commonly-used temperature for this system.)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"temperature = 0.824 / kB * epsilon"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we'll define the parameters for the double-well. We typically like to think of the height of the barrier in terms of units of $k_B T$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"h = 6.0 * kB * temperature\n",
"r0 = 2.**(1./6.) * sigma\n",
"w = 0.3 * sigma"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we set a few last parameters we'll re-use:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"density = 0.96\n",
"timestep = 0.001 * tau\n",
"collision = 2.5 / tau"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# convenience for converting distances to units of sigma\n",
"sigma_nm = sigma.value_in_unit(u.nanometer)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exploring the potentials\n",
"\n",
"To explore the potentials, we'll create a system of just 2 particles. We can manually set the positions to see the interactions.\n",
"\n",
"First, let's see how we can set this system up. There are a lot of arguments to the initialization, although we defined most of them above:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"two_particles = testsystems.DoubleWellDimer_WCAFluid(\n",
" ndimers=1,\n",
" nparticles=2,\n",
" density=density,\n",
" mass=mass,\n",
" epsilon=epsilon,\n",
" sigma=sigma,\n",
" h=h,\n",
" r0=r0,\n",
" w=w\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's make that easier. Most of these, we will re-use over and over again, so we'll put them in a dictionary for convenience:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"# this creates a dict we can re-use to simplify parameters to the potential\n",
"potential_kwargs = {'density': density,\n",
" 'mass': mass,\n",
" 'epsilon': epsilon,\n",
" 'sigma': sigma,\n",
" 'h': h,\n",
" 'r0': r0,\n",
" 'w': w}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can use that dictionary, and only specify the things that are specific to this setup:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"two_particles = testsystems.DoubleWellDimer_WCAFluid(ndimers=1, nparticles=2, **potential_kwargs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is an `openmmtools.TestSystem`. It contains the `openmm.app.Topology` and the `openmm.System` objects necessary to create an `openmm.app.Simulation`, as well as the `positions` to set after creation. The `Simulation` also needs an integrator, so we'll use one that is available in OpenMMTools."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"integrator = integrators.BAOABIntegrator(temperature, collision, timestep)\n",
"sim_2_particles = mm.app.Simulation(two_particles.topology,\n",
" two_particles.system,\n",
" integrator)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can directly set the positions now. Usually, we'll use the positions given by the `TestSystem`, and then energy minimize and run our dynamics. However, right now we want to look at the potentials, so we'll directly set the distance between the two particles in order to calculate the potential energy surface."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"# trick is that nonbonded WCA is a different force group than the bonded DW interaction\n",
"results = []\n",
"for x in np.linspace(0.9, 1.9, 300):\n",
" position = np.array([[0, 0, 0], [x, 0, 0]]) * sigma\n",
" sim_2_particles.context.setPositions(position)\n",
" V_WCA = sim_2_particles.context.getState(getEnergy=True, groups={0}).getPotentialEnergy() / (kB * temperature)\n",
" V_dw = sim_2_particles.context.getState(getEnergy=True, groups={1}).getPotentialEnergy() / (kB * temperature)\n",
" V_total = sim_2_particles.context.getState(getEnergy=True).getPotentialEnergy() / (kB * temperature)\n",
" results.append((x, V_WCA, V_dw, V_total))\n",
" \n",
"pos, wca, dw, total = zip(*results)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAENCAYAAAAG6bK5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXd4m9XZh++jZVm25D1jO3b2cPbCZBAS9ihhEwqFUEoTRqFAKaWDUdpS2lL6QYFSRspoCCthhZBAgJBFBmQvEmfZieO9LVnjfH/ITgNkeLzSq3Hu6/Jly9b7PD8nlp73nPMMIaVEoVAoFNGHQW8BCoVCodAHFQAUCoUiSlEBQKFQKKIUFQAUCoUiSlEBQKFQKKIUFQAUCoUiSlEBQKFQKKIUFQAUCoUiSlEBQKFQKKIUk94CTkRqaqrMz8/XxpjPA2WbICEH4tK0sdlGbbObAzXN9Muw01Rfy549exg8eDBWq1VTPwqFIjzYX9lAvctHYY+EoPtet25dpZSyQ29yIpRbQYwePVquXbtWG2M+LzyUAqfdA6ffp43NNpburOBHL6zmjZlFVO9YwznnnMOyZcsYP368pn4UCkV4cMrt/+SQz8G+J64Num8hxDop5eiOPDekVwCaYjBCbBI0V2luOjnOAkBVYyuDBw3ikUceITc3V3M/CoUiPGh0C0w49ZZxUqInAADYUqCpQnOzKfH+AFDd1EpuYR6//OUvNfehUCjCBycmrLj1lnFSousQOC4NmgK3AqhuciGlpLi4mMOHD2vuR6FQhAceQwxxptDdXm8nulYAcalQsUNzszEmI/ExJqqaWgHo378/v/jFL/jjH/+ouS+FIpRwu92UlJTgdIb+dkcw+ffVhZjxsm3btoD5sFqt5OTkYDabu2wjygJAGuxdFhDTyXEWqptaEUKQnJxMVZX2Kw2FItQoKSnBbreTn5+PEEJvOSGB1ydxH6wjLc5MVlJcQHxIKamqqqKkpISCgoIu24m+LaCWavB6NDfdHgAAUlJSVABQRAVOp5OUlBT15n8UXp8PgBhL4O6vhRCkpKR0e+UVXQEgvi01NgCZQClxFqoa/xcAqqurNfehUIQi6s3/2zhdbYe/bYEgUGjx7x5dAaC9AKypXHPTyXEWqppc/q/VFpBCEbU0t92V+wKw06A10XcGAAFJBU2O928BSSm59dZbaWho0NyHQqH4NlVVVUydOhWAsrIyjEYjaWn+1/nq1auxWCzfen51dTWvv/46M2fOPKFdj8dDamoqtbW1ndbk9vgAAzHm0H97DX2FWnIkAFRqbjo1Lga3V9Lg8nDmmWdqbl+hUHyflJQU1q9fD8ADDzxAfHw8d99993GfX11dzTPPPHPSANAd3N62AGDpenZOsAjqFpAQIlEI8aYQYrsQYpsQoiiY/olL9X8OYDFYVWMrlZWVLF++HLc79AtBFIpI5dFHH6WwsJDCwkKeeOIJAO6991527NjB8OHDuffee6mvr2fKlCmMHDmSoUOH8v7773fbr8fnAykxm4zdthVogr0C+AewUEp5mRDCAtiC6t2aCAZzQAJAanwMAJWNLj75eD4/+clP2L9/v2oJoYgqJk+e/L3vXXHFFdx88800Nzdz3nnnfe/n119/Pddffz2VlZVcdtll3/rZZ5991iUdq1ev5tVXX2X16tV4vV7Gjh3LaaedxiOPPMKuXbuOrBrcbjfvvPMOdrud8vJyxo8fzwUXXNAln+14fQC+sDgcD9oKQAjhACYBzwNIKVullJ3fYOueiLZq4MCtACobXKSm+lcaFRXa+1EoFCfniy++4NJLL8Vms2G325k2bRrLln2/BkhKyS9/+UuGDh3KWWedxYEDB6is7N4WsTXWhjUMtn8guCuAXkAF8KIQYhiwDrhdStkURA3+baAAnAGkHbUC6NV2CNXdPyRF6FNcepin5n3O8m8qqHEbcIkYjNLDyMH9Gdsnk8J0C5MH9cAaYzm5sQjgRHfsNpvthD9PTU3t8h3/d+lol+OXXnqJuro6vvrqK0wmEzk5Od3OrfdKMBnDI8EymCpNwEjgaSnlCKAJuPe7TxJC3CSEWCuEWBuQO+i4NGgMTBooQGVjq1oBRAE7yur5+dz1nPHEGt4sieMgyfikIEG0YBUe6lw+nv58NzPnbqX/3XO58oHnOVihakOCxaRJk5g3bx4tLS00NjbyzjvvMHHiROx2+7cy9Orq6khPT8dkMrF48WJKS0u77dvV6kaGQQooBHcFUAKUSCm/bHv8JscIAFLKZ4FnwT8PQHMVcWlQ+Y3mZk1GA0k2M5WNLtLSegAqAEQiZZU1XPnHOeyz9MRmMTJ9dDaFdieXTz0H03cO/RpdHp6Z9ynPL23kS2ceRb//kB8OsvLwTRdjMITHHWK4MnbsWKZPn86YMWMAmDVrFkOGDAFg9OjRDBkyhPPPP58777yTCy+8kNGjRzNy5Ej69u3bbd9eCR53a7ftBIOgDoQRQnwB3Cil3CGEeACIk1L+4njP13QgTDsf/RrWPA+/PuQ/E9CQMx/7nF5pcTz9w5G88cYbjB49mt69e2vqQ6EfLy/4gt98uAdsSRRaqnjl3qtJiuvY1s5/Fy7nt+9uxevIpqeoYsEDVxMXE/5Z2Nu2bWPgwIF6ywgZ3B4P28oaiTe46ZWt7eTBY3Gsf//ODIQJ9m3IbcCrQoiNwHAg+O0y49LA0wKt2h89pMbHUNnYisFg4Morr1Rv/hHE3U++zm+WVCKk5MGJCXzw++s7/OYPcPU549n2jxkUxZVzgBQufmo5+6qCe/ylCDztbSDMxtBPAYUgBwAp5Xop5Wgp5VAp5TQpZU0w/QMQn+7/HIhUUHsMVY3+dhCrV69m+fLlmvtQBJ/f/Xcpb5bEEdtcxif3ns31F0zqkh2L2cSc387g5R+Po6LBxVmPLmThyg0aq1Xoicvt3/sPhxoAiLZeQBDQauCUOAuVbQ3h7rnnHn71q19p7kMRXP69tJiXNjYwKg3W/nUGvXMyu21zfJ9U/vGDfJwtTn762lbeW7pOA6WKUKDV4wUIizYQEJUBIHDVwGn2GBpdHpxuL2lpaeoQOMy55bE5/GHBNs4fmsXcO84l3mbVzPak4f148dph4PNy65vbWbNll2a2FfoRY40FIE7Dv5VAEoUBIHAN4VLbisEqGlwqAIQ5f5+zkPfL4khsLuEfVw4PSF73lDGFPHn5QDBauPKfn7O/TNWNhDsenz+pRtUBhCq29hWA9rUAKXH+YrCqplbS0tKorq7G6/Vq7kcRWD5auYHHVzdgbCxn4e+uCOiL+cJJo7i7KAFvXCq3vroOny/058gqjk9DUzMGJIYwaAMB0RgAzFaIcQSmI6i9rRq4rR2ElFINhgkzquoamPnKWvB5mHvLZDJTkwLu87YrzuL+CwexscLDU5+praDOYjQaGT58OIMHD2bYsGE89thj+LoxjCU+Pv6Y37/++ut58803T3itq9WD9Hqpra0lJSXlSEXyypUrEUJQUlIC+AvQkpOTT6jzgQce4K9//WuHfXeF6AsAELBq4PYtoMpGFxdffDHLli3D4XBo7kcROH70+Hv44tO5syiZMYP7BM3vjIl9uGhYNn9btIN/vrE4aH4jgdjYWNavX8+WLVtYvHgxCxYs4MEHH9RFiw+BwEdiYiKZmZlHhsKvWLGCESNGsGLFCgBWrVrFuHHjdC8IjM4AEJ8e0I6gVU2t5OTkMH78eGJiYjT3owgM7204yJaWBH7QJ4bbrzo7qL6FEPz2vD7QUM6jSw+r84Aukp6ezrPPPsuTTz6JlBKn08mMGTMYMmQII0aM4NNPPwVg9uzZ3HrrrUeuu+CCC77Vh+iuu+5i5MiRTJ069ZhneevWreO0005j1KhRnH322Rw6dAgAKQSGtt2f8ePHH3nDX7FiBT//+c+/9fjUU08FYPfu3ZxzzjmMGjWKiRMnsn37ds3/XY5HeOQqaU18BhzeorlZq9lIfIyJigYXDQ0NzJ8/n1NOOUWT8nJFYNm2p5Rfz9vE8NxE/n5DcMdUtJOaYOdPFw3klx9XcPkf5/Ll/92ii46u8uB7W9h6sF5Tm4OyHdx/4eBOXdOrVy98Ph/l5eW88sorAGzatInt27dz1llnsXPnzhNe39TUxMiRI/nb3/7GQw89xIMPPsiTTz555Odut5vbbruNd955h7S0NObOncuvf/1rXnjhBRBGjPi3dU499VSWLl3KjTfeSHFxMZdffjn/+te/AH8AaE8Tv+mmm3jmmWfo27cvX375JTfffDNLlizp1O/cVaJ0BZABjYcDYjo13kJloz8A/OhHP+KTTz4JiB+Fdvh8Pq7663zqmpw8fGF/XTM4rjqriOHmMg7b8nnohXd00xHutO+9L1u2jGuvvRaAAQMG0LNnz5MGgPZKfoBrrrnme22kd+zYwebNmznzzDMZPnw4Dz/8MCUlJXi8XhACk9G/BGhfAezZs4f8/HysVitSShobG1m3bh1jx46lsbGRFStWcPnllzN8+HB++tOfHllNBIPoXAHYM8BVD63NYNF2Jo2/HYSLlJQUQLWEDgd+99x86uz5nBpbRmFeqt5yePXXP2LY3S/z/PpYbq5tJDXx2IeSoUZn79QDRXFxMUajkfT09OO2hTaZTN86gD1RC+jvDnaRUjJ48GBWrlz5re+72orAkhL85359+/alpqaG9957j6Ii/6py1KhRvPjiixQUFBAfH099fT2JiYlHBtQEmyhdAbRVcwZgFZAaH0NVYysxMTE4HA5VCxDilFXW8PLmZoz1B3nx3mv1lgNAvM3Kw9MKMdiSeHrZfr3lhBUVFRXMnDmTW2+9FSEEkyZN4tVXXwVg586d7N+/n/79+5Ofn8/69evx+XwcOHCA1atXH7Hh8/mOZNz897//ZcKECd/y0b9/fyoqKo4EALfbzZYtW/B422oAjjrYLSoq4h//+MeRAFBUVMTjjz9+ZP/f4XBQUFDAG2+8AfiDy4YNwWsPEr0rAPAHgOQCTU2nxFv4co+/H5AqBgt9rvvLXIQtlwcm9QypoS3TzypiU9MmZq/Yy7RhmQzJTdZbUsjS0tLC8OHDcbvdmEwmrr32Wu68804Abr75ZmbOnMmQIUMwmUzMnj2bmJgYxo8fT0FBAUOGDKGwsJCRI0cesRcXF8eWLVsYNWoUCQkJzJ0791v+LBYLb775Jj/72c+oq6vD4/Fwxx13kJrVNv7V5wX8E8HGjx/PggULGD3a35yzqKiI4uLiIwEA4NVXX2XWrFk8/PDDuN1urrrqKoYNGxbAf7H/EdR20J0lIO2gAco2wzPj4fL/wOBpmpr+++Kd/OOTb/jmD+cycfypOBwOFi1apKkPhTZ8c7iBsx77jCzXflY8HnoHrnXNbsb9/gNk3WG2PfkT3VMGj4VqB/0/9pVVUecx0ic1Fps1ONl/4dYOOjSIP2oFoDHtxWDVTa28/PLL/swARUjypw+3Ex9rYf5DN+gt5Zgk2MycntKIy5HDQy+8q7ccxUlw+6fBExMm84AhWgOALQUMJmgo09x0+2zgigYXffv2JScnR3Mfiu7zyqIvWbK9nJsn9yHdEau3nOPy99uuwFBfxuz1tTQ2d29WrSKweHwSfD6MIbhSOx7ho1RLDAaISw9INXC6wx8AyhucrFmzhj//+c8dHlCtCA4+n48H569HNtdw7bjQDtDWGAu3jM+G+DRu/vtrestRnACvBGTXW1DoQXQGAPAfBDdqvwJIb9sCKq938dlnn3Hvvfd+awi1Qn9+/+K7uB05XJhvID429Cu17/rhudjq9/F5hZXaxha95SiOhzBgFOF1sxe9ASA+AxoCkwYK/2sJDWo4fCjhdLUy+6saRMNh/nrLZXrL6TB/uno8wmrnv2tL9ZaiOA4mswV7nLZ1RYEmugNAAFYAVrORhFgz5Q0u0tP94yfLy7XfalJ0jbuefANpT2fGqJSQSvs8GRdNGMbp/dN4dmkx9S2testRHAO3V4bNHIB2wkutltgz/S2hvR7NTafbYyhvcJKR4c82UgEgNHB7fSytshHTWMZvrr9Qbzmd5pbTCqhtdnPDn17SW0rI8POf/5zHH3/8yOOzzz6bG2+88cjju+66i8cee4ydO3dy3nnn0adPHwYOHMgVV1zB4cP/2wG4/fbb6dGjR5fbSLs9HnxS4m51df2X0YHoDQDxGYAM2GjIiqNWAEf/oSn04931B2nwWXjy5vNDMqf+ZIzulYa9YS9rGhyUHK7SW05IcOqppx7psOnz+aisrGTLlv81elyxYgWjRo3i/PPPZ9asWezatYtt27Yxa9asI1uzPp+PefPmkZuby9KlS7uko8XlBsCgDoHDBHt7O4jAHASXN7jIzs6mpKSEGTNmaO5D0Tncbg9/XbCRAZl2zhiYobecLvPApWMQMXHc+sRbeksJCY5uubxlyxYKCwux2+3U1NTgcrnYtm0bO3bsoKioiAsv/N+q7/TTT6ewsBCATz/9lMLCQmbNmsWcOXO6pMPV6g8AFrOxm79RcAlqKwghxF6gAfACno5WqwWE9mKwABwEpzuslDe4MBgM9OjRQ3P7is7z0IvvcqgphvPyGr/X3CucuHTKWB586598bU6lvLqO9OQEvSX9jw/vhbJN2trMHALnPnLcH2dnZ2Mymdi/fz8rVqygqKiI0tJSVq5cSUJCAkOHDmX79u2MGjXquDbmzJnD9OnTueiii7jvvvtwu92YzZ0r5nK5vYAgxhxe3XX0WAGcLqUcruubPxxVDRyYFUCrx0e908MzzzyjqoF1xufz8d/1VYiGcn559Vl6y+k2d503BGGN555n5ustJSRoXwW0B4CioqIjj4/uuXMsWltbWbBgAdOmTcPhcDBu3LgutW5xt83+DqfEAojWZnAQ0BVAmr09FdR5ZEl5ww2h2W4gGvjrfz/E68jmB5n1WMLsDu1YXHf+JJ5bM4+dogcery90Mk9OcKceSNrPATZt2kRhYSG5ubn87W9/w+FwcMMNN1BeXs7nn39+zGsXLlxIXV0dQ4YMAaC5uRmbzcb555/fKQ2y7V5arQBOjAQWCSHWCSFuCrLvb2OyQGxyQPoBpR1VDJaenq6ygHTm+eX7kU1VPDLzEr2laMZvLy/iYJ2LBZu1X8GGG+PHj+f9998nOTkZo9FIcnIytbW1rFy5kqKiIq6++mpWrFjBBx98cOSahQsXsmnTJubMmcNzzz3H3r172bt3L3v27GHRokU0Nzd3SoM5xorFaAi77cVgB4DxUsqRwLnALUKISd99ghDiJiHEWiHE2oAXUNkzAxIA0u1WACoaXWRkZKgsIB1ZtmUfroQ8JqR7g9ahMRicMTCDDBvc99KSLqcuRgpDhgyhsrKSU0455VvfS0hIIDU1ldjYWN5//32eeOIJ+vbty6BBg5g9ezYOh4OPPvroW3f7cXFxTJgwgffee69TGtyhtBLrBEFdr0gpD7Z9LhdCzAPGAku/85xngWfB3w46oILsmVB/UHOzR68AMjIyqKmpobW1FYslvPYHI4G3t9Rgsxh57NZL9ZaiKQaDYFRcLQuak/jnmx9z2xXhf7bRVYxGI/X1355FPHv27G89HjBgAAsXLvzetdXV1d/73ttvv91pDU0tTswGgPCY3tZO0EKWECJOCGFv/xo4C9gcLP/HxJENDdrP33RYTcSYDJQ3OElPT8dgMKjRkDqw93A17204yBWjc8lIcugtR3P+dNM0ZEsdz3y2S28pUY/EEHaN4CC4W0AZwDIhxAZgNfCBlPL7ITmYOHr4W0J73ZqaFUKQ7vAXg82YMYPW1lays7M19aE4OT/7vzdxe7xcMTxdbykBIcEexyhHE02OnixZo++9VDTj8XrBYMBkCK/9fwhiAJBSFksph7V9DJZS/iFYvo+LIxuQgTkIjvcXg1ksFozG8CoOiQTqGprY0OjA0XiAQXlpessJGL+//hyk183Dc7tWwaoF0d7u3NlWBWwJ8hmAFv/u4XdqoSX2trvyAJwDpNv9xWC1tbXMnDmTJUuWaO5DcXx+89w7iFg7M0/vp7eUgDK4dy4FpjpKLbk0ubTva3UyrFYrVVVVUR0EnDpUAUspqaqqwmq1dstOeCWtao2jPQBo32I33RHDyuIqzGYz//rXvygoKGDKlCma+1F8H5/Px4JvmjDQzMxLIr/+4m8zL+TSp1cyf30pPxzXM6i+c3JyKCkpieqW5w3NLupcPjw2IzXl2p8pHg+r1drtiYMqAECAVgAx1LW4MVqs2Gw2VQsQROZ+/CVeRzZnJ1eHZdO3zjIyL4mBGfH834INTB+TG9Tf2Ww2U1BQEDR/och/Vuzl/ne3sObXZxzJAAwXIv/VcSJik8AUG9BU0IoGVQsQbNbVxhJjhAdmdK6aM1wRQtCHUg67jDz37md6y4k6yhucGA2ClLjwS/OO7gAghH8VEIAAkOHw7821zwVQASA41LW4eX/TIS4ZlUdWWrLecoLGAzdciHQ28szH2/SWEnXMW7gEX3MdBpUFFIYEOACU1bnIzc3V3L7i2Nz79Fs43T6uHBVdabcpCXb6Waqpis1ha3GJ3nKiijoXGFvDc+63CgABCgCZ7QGg3snrr7/O4sWLNfeh+DY+n4+Fu5sw1ZcyvGeK3nKCzj2XjEcYTfz+lc53s1R0nRYsxApta4mChQoAjmxoOAga91NJtJmxmAwcrndqaldxfF54bynSnsl5fcOrHF8rzjxlKJa6A3xdFxvVaZnBxmu2kRB+2/+ACgD+amCfB5q1bdUghCDTYaWszsknn3zCBRdcoNpBBJh/fbwF6WridzMu0FuKbtx/7Zk4zQ7W7K3RW0pU0NjsRMQmkGLr3ACZUEEFAHuW/3MAagEyHVbK6p1UV1fzwQcfcPCg9ltNCj8Hyiopj8mip6ggNTHy+v50lItH98QeY+LVVXv0lhIVlFb79/6H9c3TWUnXUAEggLUAGQlWDtc7ycz0zx8uK1O92wPFextKEaYYbj5npN5SdMVmMdHPWs/8dfvYX6ZWnIGm0euv/j29KDz/7lQAcLTN7A1EALDHfCsAHDoUvCrBaOPj3Y30y4jnyjNPOfmTI5wrx+QiTDH8/j8L9JYS8eyv8K8AshK615JBL1QAiEsDgykwW0AJVpxuH3FJ/mZkKgAEhhWbd/PV/louG5UTdhOZAsFlU8ZirD/Ikn0uvaVEPG8v9Pf4svrCM9lDBQCDwd8ULoC1AA0eI4WFhWogTIB4+JXFSJ+XiXnheRemNQaDgdPyLHgd2by3dJ3eciKasjon0u0iLzNVbyldQgUACFwtQML/agE2bdrEnXfeqbmPaKfV7WFLUzyOxgMMzO+ht5yQ4b6rz0J6PTz5wRq9pUQ01U4fwlUftj2nwlO11iTkQO1+zc22F4MdrgvP5WE48OQbixFxSVw8IktvKSFFn7wsBib4qLDl4/WpmoBA0eAxEuNr0VtGl1EBACAx138G4PNqajbd4W8IV1bv5M9//jMXXBC9+emB4tWVxUhnI7+4+ly9pYQcP7voFKpbvCzfpbKBAkWrMRaHMfxGQbajAgBAQq6/GKxB2zTNGJOR5DgLZfVOysvLWbJkiarQ1JCK2kYqY7IpMFZhj4vVW07IMWVAOjYTPPKaGkYUCKSUGOKTGRqmNQCgAoCfxLb/wLoDmptOt8dwuM5JVlYWLS0t1NfXa+4jWvloexXCZOHB69Xd/7Gwmo0kN+1jS52J8uo6veVEHNVNrfgwMGFUod5SuowKAOBfAUBgzgESrBxucB4ZCq9SQbXjra9K6J9hZ1Jhvt5SQpYbpgxGmK088sqHekuJOPaU+dttpMeHb3afCgDgPwOAgB0El9W5yMryH1KqAKANy9Zv5+v9tYxO86nc/xNw/fkTobGSD7er3kBas3jFWgBqDoZv2w0VAAAscWBLCcgWUIbDSlWTix65eUyYMAGTKbqncGrF428vA+DikWrWwokwGAwMS3DRbM/hq+3FesuJKPYergVgUEH4zp4IegAQQhiFEF8LId4Ptu8TkpALtdoHgMwEK1JCXGo2X3zxBRMnTtTcR7Th8/lYV2Ugpu4Aowf11ltOyHPnJRMQwsBrK3brLSWiOFjThPR5GVTQvcHseqLHCuB2IPTm1iXmBmQFcKQYrC58c4VDjXeXrkPaM5hcEKe3lLDgtFGDGNUzia9rLSoLTUMqmzzgrMcao84AOoQQIgc4H3gumH47RGJP/wpA4xdIdoI/PfFgrZPzzz+fWbNmaWo/Gnnmw3VIr4dfXHWG3lLChmkjevBNeSPrdquOtFpR5zZgdjfpLaNbBHsF8DhwDxB6lRMJueBpgSZti2ayEv0rgIO1LdTV1bF9+3ZN7UcbPp+kLCaHHFMDffJU9W9Hmdo3Cenz8rvn3tFbSsTgyMxjYH54/w0GLQAIIS4AyqWUJ+xOJYS4SQixVgixtqKiIkjq+F8mUJ22mUAOqxl7jIlDdf5UUJUF1D3W7quhtlVwz5Wn6y0lrMhOcRDfWMLWxlh8Go8/jVaavCZGDuilt4xu0ekAIISIE0IYu+BrPPADIcRe4DVgihDile8+SUr5rJRytJRydFpaWhfcdJEjtQDanwNkJVoprW0hKytLBYBu8tQHq7GaDJwxMENvKWHHGf2SID6VOR+t0FtK2FNV10SDy0NijN5KusdJA4AQwiCEuFoI8YEQohzYDhwSQmwRQvxFCNG3I46klL+SUuZIKfOBq4AlUspruqVeSwJYDZyVEMuhOn8AqK+vp6kpvPcN9aKx2cmnu+qIrysmLkal03aWu646E+lx88InG/WWEvas3vINAOV7wntLtyMrgE+B3sCvgEwpZa6UMh2YCKwCHhFChM4beVeJTYQYR0CKwbITYzlU62TEiBFcddVVOJ2qO2hX+OdbnyCs8Vw6pqfeUsKSvMxUEltK2dXqwO326C0nrNmyx7+S75+brrOS7tGRALBWSvl7KeVGKeWRzUMpZbWU8i0p5aXA3M44lVJ+JqUMvdaYAaoFyE6wUtXUymlTzmDOnDmkpKRo7iMaeGP1XqSzkdsuU9k/XWXG1KEIWxLr9tfqLSWs2XWoCoChfcK3BgA6FgBOetompXRroEV/EvMCsgLISvSngh5qmwugDuE6T019IxXmTHJFJfFauW93AAAgAElEQVQ2Nfmrq9x0QRGxZiPvb1JnUd2hpMpfBFbYO3w7gYJqBfFtkvKhZq/2tQBtqaDFZTUkJSXx2GOPaWo/Gnjhw1UIi5XLx/XRW0pYY7OYGJsTyxurdtPiVDODu0p5kwfRUovFHN5nUR0JAMOEEHuEEO8KIf4ohJguhBgihDAHXF2wSe4F7iZoLNfUbHsxWLVT4vP5KCkp0dR+NFBKKvYYIz+ZNllvKWFPb3MtLsw8M0/NCegqaT370Sc7/LdyOxIANuJP4XwSqALOAl4EKoUQmwOoLfgkF/g/12jb3a+9HcShOic5OTkqAHSSVo+PxdsOc3ZhFjZrmOfdhQC3X34G0tXM3FWqN1BXqfMYKQzjJnDtdGj9IqU8CBwEFrV/T/h78EbWejypLQBU74G8UzQzazUbSY23cLC2hR49eqgA0EmeevNjGpwexmRG3qJTDxLscWT6yjloSKehqUVNU+skrW4PB2uaicsP/15UHVkB/PNY35R+vtFYj74k5oEwaL4CAH8twEG1AugSb6z6BtnazLkj8vWWEjFcNqYAQ0wcT771id5Swo6txQfwIajYt1NvKd3mpAFASvmtxm1CiIFCiDPbvo6sWweTBRw5/hWAxmQnWjlU28J5553HNdeEf9lEsHC6WimRyaS5y3HE2/SWEzHMumQKtDazYn+z3lLCjo27/TdwvbOSdVbSfbqSBfQ0UCiEeBd4SQjxkMaa9CU5P3ArgNoWLr30Uh555BHN7UcqL7y3FBHr4Lyh4b/fGkrE26xcekpf9nkSaPWotOTOsH2/P0lkcEF4N4KDrgWATVLKv+Nv7HY5EP5h8GiSCgK2Amhq9VLv9NDY2KiqgTvI3BU7kW4Xt102VW8pEcd5Q7Kod3r4fJuqCegMe8vrABjZP19fIRrQlQBwihDiSSBPCDEEiKyBrMkF0FwJznpNzWa3FYN9tnoDdrud998PrYFooYjPJ6mK7UGuuYG0JIfeciKO8X1SEB4X9z83X28pYcWhOifS1UR2Wvjf+3a6ikFKOaZtsMso4HIgshqzJB2VCpo1TDOzPdoCgM+aCKAOgjvA1wdqaPSaePhK1fohEFjNJlJdhyg1ptLsdKkU2w6S038o3qrIaOjY4RWAEOLh9q+llCVSyneAB0Oyp093SD4qFVRDcpL8B5i1HiNWq1UFgA7w8qebMRkEUwaGd8OtUObCETkIazz/fuczvaWEDVUtkj4RcAAMndsC6iGEmN7+QAiRBnysvSSdSQpMMVhqvAWr2UBpTYtKBe0APp+P+ev2Yqvfh8Oq8v8DxW2XTUW2Onlj1S69pYQNxWU1mFsb9JahCZ3ZAvop8JEQYjcg8VcD/zIgqvTE6gBbiuYrACEEOUk2DtQ0qwDQAd5cshriUjgtI7x7rYQ6SY54Ut1lHDAk43S1hvWA82BQXl2HUxqpPah9oogenPTVJYR4CfgK+Bq4Bfgv4AGmSSkj87YhqSAgqaA5SbGU1LRw00030draqrn9SOI/n2xA+rK44/IpekuJeK45bTD/WNvEmr01TOyvJq2diNVb/e0zemUk6KxEGzqyBfSftufdgP/NPx+oAa4RQlwWOGk6ktwLqoo1N5ubZKOkpoXp06dz3XXXaW4/UvD5fGytt2BrLKV3TqbeciKemT+YSKzZyKJtQZzBHaZs3F0KwOD8yPi77Egl8CdSyseklNdJKYcDqcBdwG5Au4Y5oURqX6gvgVZtT/pzkmKpa3FTUdvItm3bVC3Acfho1SakPZ2J+fF6S4kKYi1GTukZz1urd9OqJoWdkJ2l1QCMHligsxJt6MhM4G/l+UspPW3TwV6WUt59rOeEPaltY46rtO2W2J4J9MaCTxg0aBCbN0dWM1Wt2NEUiwDuvupMvaVEDVmewzT7TLy04Au9pYQ0B6qbkW4n/fLCvwoYOjgTWAhxmxDiW6NvhBAWIcQUIcR/gMjaz0hpDwDa9rrLTfbXAhjsaQAcOKD9+MlIYOGWMkbnJ9EvLzKW2eHAzy6bgvS4+O8X2/SWEtLkF46iV0YCBkNkzNLqyG9xDuAF5gghDgohtgohioFvgOnA36WUswOoMfik9AYEVGobANpXAO4Yf1Xrvn37NLUfCXy2bivbyxoYk6WyUYJJZmoSiS2H2O1y4PF49ZYTspTWOumVHhkHwNCxMwCnlPIpKeV4/FW/U4GRUsqeUsqfSCnXB1xlsDHHQmKu5gEgyWbGZjFS4xLExcWpAHAMnn5vJQCTe0fOiyxcOGtgGiIuiVc/Wq63lJDE5/Ox61ANNFfpLUUzOrWOkVK6pZSHpJS1gRIUMqT01XwLSAhBbpKNAzUt9OzZk71792pqPxJYV+7FVFfK2MK+ekuJOm6/fCrS6+H9DaV6SwlJ9hwsxyNMtFZHTvO8oFXZCCGswFIgps3vm1LK+4Plv9Ok9oWvVvkHxGt4xt1eC/CHP/yBxMREzexGAqs3f4PH0YNx1sh5gYUTORkpTOqfzp6qFqSURFpuR3dZs9VfG9Q3OzLaQEDXuoF2FRcwRUo5DBgOnCOECN000tS+/gHx9Qc1NesPAM1MmzaNyZMna2o73Hly/jIAZl1YpLOS6OXCYTmU1LSwubRObykhx+a9/huTob0iZzZFZ5rB3SqESOqqo7YRko1tD81tH7Kr9gJOgDKBcpJsNDg97CkpY9GiRaoW4Ch2t9qJaalk8qjBekuJWqYMSAXp476nXtdbSsix+1ANAGMG9dJZiXZ0ZgWQCawRQrwuhDinK7n/QgijEGI9UA4sllJ+2VkbQaO9FkDjg+D2VNC3Fy3l7LPPZteuyOym0VnKG5wcdMequ3+dSbXHEttQwuY6Mz6fmhR2NBUtEulspEd6it5SNKPDAUBK+RugL/A8cD3wjRDij0KI3p2w4W2rJs4BxgohCr/7HCHETUKItUKItRUVOpam27PAEh+AAOBPBTUm+HuuqEwgP2+t2oWUcG5hZBTYhDMT8uOR9nQWrtygt5SQIrvfUIZE0PYPdD4LSAJlbR8eIAl4UwjxaCft1AKf4a8x+O7PnpVSjpZSjk5LS+uMWW0RAlL6aL4F1DMlDgCXWdUCHM3/zfsCU0s1/TJU+we9uf3SyUjp498L1+otJaQoqW4mL8WmtwxN6cwZwM+EEOuAR4HlwBAp5Sz8k8Eu7cD1aUKIxLavY4EzgO1dUh0sUvtBxU5NTcbHmEiNt1DjNmKxWFQAAHaXlNEc34P+tmaVeRICFPbOJaahlA1V6v+iHY/Hy57yelxV2iaF6E1nVgCpwCVSyrOllG9IKd0AUkof0JGpYFn420psBNbgPwMI7cG46QP8TeE0ng+cl2xjX1UzeXl5KgAAj7+xBGEwcv3U4XpLUbRxzqAMfI4siisiY/BJd9m4ax/SYMTojKwSqA7XAUgpf3eCn520gYiUciMwoqP+QoL0Qf7P5dsgb5xmZvNT4lhVXMULL7yArttcIcKn39SCwc2lU87VW4qijXt+eDbv/vlTPtpSzqzJdr3l6M6XW/zt4Qt7Rta8hA4HACHEncf4dh2wLiLbQcBRAWCrpgGgZ0oc89aXMuaUyVjNRs3shiMlh6toiMumH4cipsFWJJCTZGNQZjyvLdvOrMkdzvOIWDbtPQwkMmZgvt5SNKUzr7jRwEygR9vHTcBk4N9CiHu0lxYCJOT6M4HKte2Q2DPFhpSwcuNOnn/+eVwul6b2w4lVB5oQRrNK/wxB7LW72NcIa7dq2xY9HCkub0B63IwaGDk1ANC5AJCCvwncXVLKu/AHhDRgEv600MjDYIC0Af4VgIb0bMskWLxyPTfeeCP79+/X1H44sXhbORmOGKZNGKa3FMV3mHmBv1D/iflqRkCjsGFurcNijqwZ1Z0JAHnA0YNs3UBPKWUL/jYPkUnGIH8AkNoVLbenghKfCkRvKmhFTT0fbzlIUa4Ng0FlnIQap48ejKH+EKtKVLV6fGY+k0cN0luG5nQmAPwXWCWEuF8IcT/+VNA5Qog4QNtb5FAifZC//WuTdkVpSTYzdquJFpM/5724WPv5w+HAE29+ghcj/WKb9ZaiOA5Dk3047T3Ysjt6hxdJKdlb1fS/G7cIokMBoK3tw2zgJ0At/sPfmVLKh6SUTVLKHwZOos6kD/R/1nAbSAhBfkoclS4DFouF3bujc491wcaDyJZ6brhwkt5SFMfhxnNGI4SBf32wSm8purHxm3043T5ayiNvpd6hDS0ppRRCzJdSjgLWBVhTaNGeCXR4K/SarJnZvBQbW0rrKCgoiMoAUN/YTIU5gxzfYawxavpXqHLeqcPJ+aKSalveyZ8coaza7F+hZ9giL0utM7/RKiHEmIApCVXi08GWqvlBcH6KjZKaFt6aN59nn31WU9vhwD/fWoKwxHLJmMjKqog0DAYDPxiRx8riKmqaWk9+QQSycY+/DfTYQfn6CgkAnQkAp+MPAruFEBuFEJvaqnojn/SBAUgFjcPjk8Sn55KcHDkDJjrK8v1NyNZmfnrx6XpLUZyESQV2vD7JH2a/q7cUXSiuaE8Bjbx6iM4EgHOBXsAU4EL87R8uDISokCNjsH8F4NNuWHZ+24HSF19v5/7776empkYz26GO2+vjgC+Ji8f2Id5m1VuO4iSM7ZsFTVV8tOWw3lJ04VCDB0NLdcSlgELnAsB+YCJwnZRyH/5hLpFVF308soaBuxmqtNurL0j1B4D1xYd46KGH2L49tPviacmyneXUtbg5b4hq/RwOGAwG+libqLf14GBFtd5ygo6MT4vI/X/oXAB4CigCprc9bgD+qbmiUCSrrUjpkHb90VPjLTisJlrMCQBRdRD8u2ffQnhbmdg3VW8pig5y9aRBCJOZx1//RG8pQUVKSaslgfNPG6u3lIDQmQAwTkp5C+AEkFLWANGRvpHaD4wxcEi7lkdCCHqlxVPVakQIETW1AK1uD/t9SSS1lBJribwldaTyo3MnIJtrWbStXG8pQeVwvZMWt5f8CJsD0E5nAoBbCGGkbY6vECINiI6ZcUaz/xxAwxUAQO+0ePZUNdOjR4+oWQHM/mApIjaBc9Tkr7DCZDIyOMFNQ3wuza0eveUEjTc+8rfB8NZG1hyAdjoTAP4PmAdkCCH+ACwD/hgQVaFI1jAo26hpS4je6XEcrndR0G9g1PQDem3ZdqSnlZ9dNkVvKYpO8pvrLsCLkc936DiqNcisLy4DYGSfHjorCQydmQn8KnAP/jf9g8A0KeUbgRIWcmQNA2cd1GpXDdgr1d8K4s9PvciSJUs0sxuqeDxedrvsJLQcJDM1SW85ik4ytiCZJJuZ15ZHT8LCrvJGpKuZwt65eksJCJ0ZCRkDjAQS8HcGvVwIcdwhMRFHAA6C+6T7M4HKmmVUjEL8en81Ii6Zy4v66i1F0QVMRgPxtcV8trOShqYWveUEhXKnwOKsjthZFZ35rd4BLsI/DL7pqI/oIH0QGEyaBoC85DiMBsGXW/dxww03RPw5wMfbKzEZBD+7RBV/hSs/GJmLsNh48q3oyAZqNtlJtkTumUdnAkCOlPJKKeWjUsq/tX8ETFmoYbb6ZwNoGAAsJgN5yTb2Vrfw4osvsnFj5BZW+3w+Xl+5k1G5dhJsZr3lKLrIzIun4HM18c66yGuM9l1qm10IWxKjInT/HzoXAFYIIYYETEk4kDUMDq7X9iA4LY7KVv9YyJ07d2pmN9SY//laatwmUpqj47A7UrHHxZLhKeeQSKXZGbljQAD2Vvm3uaZNOUVnJYGjMwFgArBOCLEj6noBtdNjJDRXQq12b2K90+LZX+MkMysroquBX1z8NdLn5Y7L1PZPuHPh8ByENZ7n3v1cbykBZeNef81Dn/R4nZUEjs72AuoLnEW09QJqJ6etGWrJGs1M9kqLo9Xjo/fQsezYsUMzu6HG5jozsY2l9OuZrbcURTe57fIziDFCmSlTbykBZc4HS5BeN9mOyK13PWkAaB/43tb/Z6yUcl/7B/DTjjoSQuQKIT4VQmwTQmwRQtzeddk6kT4YzDZNA0DvNP/dRfaAkfh8kVlX99HKDUh7BqfmRWY1ZbSR5IjnjEFZLNpajten3XZoqFHa4MXQXBXR8yo6sgK46qivf/Wdn53TCV8e4C4p5UDgFOAWIUR4Ddk0miB7pKYBoG+6HYDTp13NqlWROXVpzrLtSOnj9ksm6y1FoRGn5tmobHTxysLleksJGA3YSBCRPQ+5IwFAHOfrYz0+LlLKQ1LKr9q+bgC2AeF3vJ4zGg5tBLc2fxgJNjOZDivfHG7UxF4oUmbOYkROAsP69dRbikIjpg7MQHpa+c+SyDwGbGx24otLIcce2f2qOhIA5HG+PtbjDiGEyAdGAF925XpdyRkDPre/LYRG9Mu0s6W0lrPOOou5c+dqZjcU2F3RyPayBi4aGZmVlNFKZmoSCS2H2O2yR+TW5dKvtyEMRgb0iOyK9Y4EgGFCiHohRAMwtO3r9sedTgsVQsQDbwF3SCnrj/Hzm4QQa4UQaysqQrDnSM5o/+cDqzUzOSDTzp7qFr5Ytpw1a7TbXgoF7nvK3y1kan/V+jnSmNo/BRGXzGuLVuotRXMqPTEAnHfqMJ2VBJaTBgAppVFK6ZBS2qWUprav2x93qqJHCGHG/+b/qpTy7eP4e1ZKOVpKOTotLa0z5oODPRMS8jQ9B+iXYafV46PXsHERlwq65rAXS90B8lLtektRaMwdl09Bej289OkmvaVoTlmLAZNBMGFoZLctCVqDC+FvdvM8sE1K+Viw/AaE3DGaBoABmf43x6wBoyIqFfTTtVvwObI4pUeM3lIUAaBnVhqJreWUiFSkhsWRocCqbfvJS4zBYorMHkDtBPO3Gw9cC0wRQqxv+zgviP61I3cc1JdCjTbl8H3S4xECbFl9KC4uxuWKjArLp97zbw387JJJOitRBIpfXn02jcSy7VCD3lI0ZV1xGfUHIms1fiyCFgCklMuklEJKOVRKObztY0Gw/GtK/gT/533apMBZzUb/kPjEbE4//XRqa2s1sas3X1VIzPUljB7UW28pigBx1uAMDALmrYmciXYHyioRcSn0TrHqLSXgRPb6JlCkDYTYJNirXQ50/ww7jQY7H3/8MRkZGZrZ1Ys9FQ14Hdmc2T9FbymKAJIaH4PdWc7zi77SW4pmLFq9GYCRvSO70hlUAOgaBgP0HA97v9DMZL9MO3urmnC6vRGxn/rhlsMA/Pq683VWogg0I9MN+OwZLP4yMg6DV23z9/qaOmqAzkoCjwoAXSV/gn86WO0BTcz1z7Djk3D+9B9z9dVXa2JTT15btoMh2XZ6JMbqLUURYG5vO+P514LwK+s5FtvLGpCuJkYOKNBbSsBRAaCraHwO0L8tE8gTn8FXX4X3cnrZ+u3sbwR77S69pSiCwPD+BZjrS/i6IvxXrgAJPQfSP8MesVPAjibyf8NAkT4YrImwd5km5gpS44g1G7H16MeuXbtoaQnfkXv/fMcfFG+bNl5nJYpgMSbTjNeRzfIN4Z3GLKVkX42bcQOio3JdBYCucuQcQJsAYDQIBmU7cMam4fP52Lp1qyZ29WDNYS+mulKKhvbXW4oiSNx+6WkAYb8K2Li7hAaXh8zYyGtvcSxUAOgOBROhZo9m9QCF2Q4OOU2AYNOm8DxQW7lxJx5HD8ZkGvWWoggi4wb3ZlCWgyU7q/SW0i3mf74OgETf97rURCQqAHSH3lP9n3drMyB7cI8EnB7JFT++hZ49w7Nz5jML/D2Sbps2QWclimAzscDOV/tr2bBzr95SuszqXYeRXg/nTxiut5SgoAJAd0jt6+8LtEubAFCYnQDA1bfey+mnh+foxMPWPAamxXDqMLX9E22MTPN3h3/8rfAdFbm3zouxqYJEe+SOgTwaFQC6gxDQZyoUfw5ed7fN9c2Ix2IysLm0jtLSUg0EBpcdZQ1sL2vgqqI+ektR6MDZRcMw1B9i2YHwHKLi8/lotCSTaYmMViwdQQWA7tJnKrQ2aNIe2mw0MDDTzsIvt5CTk0NlZaUGAoPHvU+/iZA+zi2M/ApKxbEZmyFwO3L4Yn349dH5avtehNVOYbZDbylBQwWA7lIwCQwm2PWxJuYG90igyuefnRtOB8E+n4+vq4zYGg6Q7oj8HiqKY3PXZf5soMff1q5KPlhUeP1Fi9f/YIrOSoKHCgDdxZrg7w6qUQAozE6g2QOmhAzWr1+vic1g8MrC5RCfwtkDkvWWotCRMYP7YK0vYXuLI+xammwqrcNkEFHRA6gdFQC0oM9U/4jI+oPdNlXYw7/8zBw0lnXr1nXbXrB48ZNNSI+Lu6efpbcUhc78+tqzaDLGs+VgeKVSzvt8HSkmF1Zz9KQwqwCgBf3bGp5t/6D7pjLtWIwGMgefwtq1a7ttLxg4Xa0UexJJdh4iO02tAKKdC4f1wGwUvLFGm/qYYODz+TjYYsKkwU1cOKECgBak9YeUPrD9/W6bijEZGdzDQUKv4Tz44IMaiAs8n249iIhN4MpTeuktRRECJNosZPiqmb1kE61uj95yOsSqTd8grPFRdQAMKgBogxAw8EJ/W4jm6m6bG5mXxP4mwcWXXq6BuMDz8Td12K0m7rhSbf8o/JyWH4uwJfLv+Z/qLaVDzF+2EYAzRvXTWUlwUQFAKwZcCD4P7Pyo26ZG5iXR6vHx5scr+frrrzUQFzhq6hv5YEMJ5wzOjKq9U8WJuXv62cjWZl5ZtlNvKR3iy90VyFYnF04YobeUoKICgFZkjwB7tibbQCN7JgJw32P/5pFHHum2vUDy1zkf4fRCHuV6S1GEEEmOeDK95Rw0pFPX0KS3nJNSbUjE7q7EGmPRW0pQUQFAKwwGGHC+vy1Ea/f+4LMSYslKsJLUd1TIHwS/s+Ewsqman04Lz9YVisAxvag3IsbGY68t0lvKCWl0eWgyJ3DDD6Lvb1gFAC0ZfDF4WjTJBhqZl4TbkUNxcTHV1d0/VwgEG3buoyE+h8L4Zixmk95yFCHGzIunEGfwsId0vaWckK/31+CTMLpnkt5Sgo4KAFqSVwQJubBxbrdNjchLpFFaMMYnh2w9wJ/mfIwwGLn7UtX5U/F9rDEWrp3Yn+V76ihvCN3+QH9/+V2QkiHZdr2lBB0VALTEYIChV8DuJdDYvT3xkW13IzE9BvLll6E3a1VKydoqM5a6A5w+erDechQhymWjeuD1SR6Y/aHeUo7L9koXhoYykuKjr4VJ0AKAEOIFIUS5EGJzsHzqwtArQfpg05vdMlOYnYDVbOCHP7+fO+64QyNx2vH1gVo8thR+Pq1IbymKEKZPuh1L3QE+3F6Dzxd6U7Za3R6arGnkWFv1lqILwVwBzAbOCaI/fUjrD1nDu70NZDEZGJOfzN5mC/Hxodeb/I21JcSajVw7Wd39K07M2f0cSHsGL3+ozfhULXl/2VcIi41xvVL1lqILQQsAUsqlQGieZmrN0Cvh0Ho4vKVbZop6p7DzcCO/+O1D7NgROsO2a+obmbN8J4MdrcTHqMNfxYn59bXnIludPLs49Lrbvr3c/xqdfsZonZXoQ8idAQghbhJCrBVCrK2oqNBbTtcYdhUYY2DtC90yc2pv/13JU29+zMKFC7VQpgl/fGkBmK1MKbDpLUURBmSmJpHtLaPEkEFFTWg1iKs2pxHnbWTkgOhsYxJyAUBK+ayUcrSUcnRaWprecrqGLdmfErphLrgau2ymMNuBPcZE2uDxfPFF6PRXf2dLNTRWcJPK/Vd0kBunDkZYYpmzLHQGxbR6fOxpNHHZhOjdxgy5ABAxjPmxf1LYpte7bMJkNDCuVzLW/GF88cUXIdFf/b2l62hNyGV8hsRkUq0fFB1jxgWT6JUWx2cHuj86VSs+Wb+bFreX8X2ic/8fVAAIHDljIGMIrHkBuvHGXdQ7FafJTlWLj2+++UZDgV3jL/O/RHpaeWjGuXpLUYQRBoOBa8b15Kv9tXy2sVhvOQA8/fYSpM9LYXqM3lJ0I5hpoHOAlUB/IUSJEOLHwfKtC0L4VwGHN8G+5V02U9QrBYDEfmPZvXu3Vuq6RKPLw+HYngy2u+idEz1TkxTacNHQDKTbxa+eD42agG01YGk8RHZa9FUAtxPMLKDpUsosKaVZSpkjpXw+WL51Y+iVYEuFZY932cSATDtp9hiuuOMhzj1X37vu+V+X4vLCwzPO1lWHIjxJcdjI9R3ioDmLPaWHddVScriKVnsWAxJ1laE7agsokFhscMpM2LUYyrqWAmcwCKb0T+eLbypp9ehXSOPz+fj7e2vpm2plRG6Uv2oUXeYX08YhTDH85oUFuup4/v0vEAYjF4zuq6sOvVEBINCM+QlY7N1aBUwdmE6Dy8PwMy/VrS3Efz5YRpXXSl9RhhBCFw2K8Oei00ZjqTvA8sMCj8erm44FG0uRzkauPW+8bhpCARUAAk1sIoyeAVvehqqu7eFP6JtKjElw2JjOhx/qs3/6xEebkK4mHphxgS7+FZHDpUNTIT6Nf727VBf/Hq+P5sTenFrgwGaN3gNgUAEgOBTdCiYrLHm4S5fbLCbG90kjcfBEPtShIGzFhh1U2XIZaKkmPTkh6P4VkcWvr7uAJKuBL2vjdPG/em81DS4v100drov/UEIFgGBgz4CiW/yrgINdG/E4dWA63tgk1u8uo6qqSmOBJ+Y3L38CUvLHGWrmr6L7xNus3HhaX5burGDrwbqg+//7659gEpJJ/cK00FRDVAAIFqf+DGKT4eMHunT51AEZAFj7jGXRouBNWKprcbPP2INMd2nUlssrtOeacT0xSg8zHn01qH59Ph+rS13E1u/DZlF9rFQACBZWB0z6BRR/1qXB8ZkJVkb1TKTHqdPIzc3VXt9xePXLfXiFiRd+8cOg+VREPgk2MwW+g5RZcvhyc/AKHOd9thYRn8LkPtGb+380KgAEkzE3Qmp/WHA3tDZ3+vKLR+TQbE4gqdeQAIj7PhU19TyzZCeT+qUxuIfa+1doy5+uPxOk5FezF0J/9j8AAA4dSURBVAfN5zML1yG9bm6/VPWxAhUAgovJAhf8HWr3w9JHO335+UOyMBsFL32+ja1btwZA4Le548m3qG+VXNRbLZUV2jNmcB+yWkvYTSabdu0PuL9mp4udrgQSm0vpk5cVcH/hgAoAwSZ/PAy/BlY80enisKQ4C5P7p/Payt3c+6v7AiTQT1VdA8sqY4it38elp40MqC9F9PKX66eAwcTtz3wQcF/vfLkTEZvAFWOCt4Ua6qgAoAdnPgS2FHjrRnC3dOrSi0f0QNgS+XRLKTU1NQESCHc88SYiNoGfnzUgYD4UiokjBjI6xUNJTE8O1XXutdBZlh5oJSXOwi+uOT+gfsIJFQD0IC4Fpj0NFdth0W87demUAenEmQXWwqnMmzcvIPL2l1WytNJGbP1+bpo2JSA+FIp2Hv/JOfgkPLlkV8B8HCivYfHWw0wb0QOLamN+BBUA9KLPVDjlFljzb9j6Tocvs5qN/LCogLj+43n6pa7PGjgRj7yzDmJsPHzZqIDYVyiOJjfZxnkDEnh11V4WrlgfEB/3PDMPj08yMcccEPvhigoAenLG/ZAzFubNhEMbOnzZ9afmI4RgjymXPXv2aCqptLaFjw/4uGhoFpdOGaupbYXiePxsci9kawt3vbJCc9sNTS2sqLQQW7+fycOju/nbd1EBQE9MMXDlK/4CsTnTof5Qhy7LTozlnMEZZJx6CWnZ2h5o/Xz25wDcc94gTe0qFCeiT14Wk1KaaHL05C+vaNsp9L5n5yFsSdw0sUBTu5GACgB6Y8+A6XPAWQcvXQSNFR267KeT+9LU6mPumgP4fNq0iX70lQ9YXeZlfGIDOUlq4LsiuDx953RoKOefq8qpqmvQxKbH4+X9XS6M9Qe5/UrVyuS7qAAQCmQNhatf99cHvHQRNFWe9JLhuYmMy0/ij2+v4b7fPdhtCRU19Tz1ZRWioYwnbv5Bt+0pFJ0l3mblF5NzID6NGX97UxObT7+3AmlP54ohSRgM6u3uu6h/kVAhfzxc/RpU74bnz+xQ6+h7zh2A1xLHC8v3Ull58qBxIq54+BWIS+E3Z/cm3mbtli2FoqvccvmZFKV52exO56v93Utzdnt9zNvtJcdh4v4b1E3NsVABIJToNRmue8+/HfTcGVD8+QmfPqpnMqfm2bCO+AH3PfzXLrt98Pn57DH3pLdnHz/+weQu21EotODZW84l02Hl9jlfceBw1zvfPrloM8WVTTx0yXCsMRYNFUYOKgCEGrlj4ceLIS7Vvx30ye/B6znu0x+dPg6T2cw7pVbWr+98Ct3+qmZe22MmtqWcdx66oTvKFQpNsFvN/O3yIRyoauK8h17r0uSwFRt28PjH39DH5uT0/ukBUBkZqAAQiqT0hps+gxE/hC/+Cs+eBvtXHfOpOUk27jqjN7G9x3D7k291ys3Bqnpu+M8azCYTi+6/Um39KEKGoj7pTE2upcGez7TfPtepaxubnVz3r8/A5+ORy4arEaYnQAWAUMUSBxf9E658FVpq4YWz4fUfQdnm7z111tRBjMy0UJpexMrdHVsy7y+rZOJv57KnvJ5nfzSa3GSV9aMILf59zzVkO/eyWeYw7TfPdijbzefzcfZ9z+N25PCjgUZGD+odBKXhS1ADgBDiHCHEDiHELiHEvcH0HbYMvABuXQ2T7oHdn8Iz4+Hli2HLfPC4ADAYBLNnTiY/NY4f/2cNDz87FynlcU1+tb2Y0x+ahycujVlDTJzSKyVYv41C0WEMBgOf/eWnZLbsZb2nB//f3r3HyFWWcRz//na3SwvFtktLEUq7sHIxINW2BFCRgtz/gIAoAtKIl1UBNSHKxaBIGsNNjTFYsRCyXBKIIlZQxCKRS8SmN3qltoK1ZamkLS0tpaXd3Xn845wlw8JOh2XnTGfO75NssufMe+Y878zkPOe8c+Z5L771d/QU+v9cFwrB6dfO5JWhrRytTqa3n59htLVJpQ4Ug7ojqRFYBZwGdALzgIsiot+6xlOmTIn58+dnEl9N2LEZ5t4JC+6BrZ3QPDwpKXHIZ+DASazf+yOcftsTbO5qom37Ch68sZ39R7e8vflbO3fxvV89xKNrG6GhgSsm7c3VLoxle7ju7h7O/mEHq+IAjm0dxbWntTG5bew72ry8aTs/+MNSnv33Rtq61/LErd/I7W2fkhZExJSy2maYAE4AfhwRZ6TL1wFExE39beME0I9CT3I18K9HYeXjsO3VZH1jMz0jW1m4oYHVhbG83tVEc3QxqqWFoR8azUsbtvHmrgKNO7dw7qQJtI3zl2NWO5a9spWnVq5nR1fyGR7VHDRIbNkZ7NprBE2NDUw9YgwTx42k5of9h+wDJ1w+oE331ARwAXBmRHwtXb4UOC4iruzTrh1oBxg/fvzkNWvWZBJfzYpIfkC2biG8shA2r4at69i+/r+w6w2G0M0QDc4vhc0sI/vsD98f2FSZ7ycBZDnV03vl5Hdln4iYCcyE5Aqg0kHVPAlGTUj+jjrv7dXv+Eq3UIB4/7fSme3JCoUChULQ5PLOA5ZlAugEiiuXjQPWZbj//GpowDd8Wb1paPSn+oPK8vWbBxwm6RBJzcAXgUcy3L+ZmRXJ7AogIrolXQn8FWgE7o6I5Vnt38zM3inLISAi4jFgcIt9m5nZgHgIzcwsp5wAzMxyygnAzCynnADMzHLKCcDMLKecAMzMcsoJwMwsp5wAzMxyygnAzCynnADMzHIqs/kABkLSBqDWJgQYDWysdhAZc5/zwX2uDRMiYkw5DffoBFCLJM0vdzKGeuE+54P7XH88BGRmllNOAGZmOeUEMPhmVjuAKnCf88F9rjP+DsDMLKd8BWBmllNOAAMg6W5J6yUt6+dxSfqlpBclLZE0KesYB1sZfb4k7esSSc9Jmph1jINtd30uanespB5JF2QVW6WU02dJUyUtkrRc0tNZxlcJZXy2R0h6VNLitM+XZR1jpTgBDEwHcGaJx88CDkv/2oFfZxBTpXVQus+rgZMi4hhgOvUxdtpB6T4jqRG4hWSu63rQQYk+SxoJzADOiYijgM9nFFcldVD6fb4CeCEiJgJTgZ9Jas4gropzAhiAiHgG2FSiybnAvZGYA4yU9OFsoquM3fU5Ip6LiM3p4hxgXCaBVVAZ7zPAt4HfA+srH1HlldHni4GHI2Jt2r7m+11GnwPYV5KA4Wnb7ixiqzQngMo4CHi5aLkzXZcXXwX+Uu0gKk3SQcB5wB3VjiVDhwOjJD0laYGkadUOKAO3Ax8F1gFLge9GRKG6IQ2OpmoHUKf0HutycbuVpJNJEsCnqx1LBn4BXBMRPcnJYS40AZOBzwLDgH9KmhMRq6obVkWdASwCTgHagCckPRsRW6sb1gfnBFAZncDBRcvjSM4e6pqkY4C7gLMi4rVqx5OBKcCD6cF/NHC2pO6ImFXdsCqqE9gYEW8Cb0p6BpgI1HMCuAy4OZJ75l+UtBo4Ephb3bA+OA8BVcYjwLT0bqDjgS0R8b9qB1VJksYDDwOX1vnZ4Nsi4pCIaI2IVuAh4PI6P/gD/BE4UVKTpL2B44AVVY6p0taSXPEgaSxwBPCfqkY0SHwFMACSHiC5G2C0pE7gBmAIQETcATwGnA28CGwnOYOoaWX0+UfAfsCM9Iy4u9aLaJXR57qzuz5HxApJjwNLgAJwV0SUvE12T1fG+zwd6JC0lGR495qIqLUKoe/JvwQ2M8spDwGZmeWUE4CZWU45AZiZ5ZQTgJlZTjkBmJnllBOAmVlOOQGYmeWUE4DVvLQWf299+sWSrpLUkD72XIntRkq6PLtI37X/YZKeTktKl2rXLOkZSf7hpg0qJwCrBzsi4uNpffrTSH6FfQNARHyyxHYjgaolAOArJKWVe0o1iohdwJPAhZlEZbnhBGB1Ja1P3w5cmdZi2gYgaR9Jf06vEJZJuhC4GWhLrx5uS9vNSsscL5fUnq5rlbRC0p3p+tmShvXuU9K0dCa0xZLuS9d9SdLc9Ll/089Z/iUktXV6n2dieqb/gqSCpJB0Y/rwrLS92aBxKQireZK2RcTwPus2k1RsfCkihkv6HHBmRHw9fXwEMAr4U0QcXbRdS0RsSg/w84CTgH1J6jpNiYhFkn4LPBIR90s6iqQI3qciYqOkFmAscCtwfkR0SZoBzImIe4v20wysjYgD0uWhJCWHp0XEXEnTgaHA1RERaQJ5NSLGDPoLaLnlKwCrV30L9C8FTpV0i6QTI2JLP9t9R9JiklnNDiaZ1hNgdUQsSv9fALSm/58CPNRbHCwiNpFUjpwMzJO0KF0+tM9+RgOvFy2fCiyMiN4Sw0uAlrQEMekw0S5J++6+62bl8ZdKVnckHQr0UDRNY0SskjSZ5PuBmyTNBu7ts91UkgPxCRGxXdJTJGfhADuLmvaQTIYCSaLpexkt4J6IuK5EmDuKnhvgaJIk1WsSsLDPNnsBb5V4TrP3xVcAVlckjSGZovH2KBrflHQgsD0i7gd+SnKAfYNkeKfXCGBzevA/Eji+jF0+CXxB0n7pflrSdRdI2r93naQJxRul8yc3pkM/AK8Bx6TtDwfOBx4sin8/YENEdJX3Spjtnq8ArB4MS4dahpBM1n0f8PM+bT4G3CapAHQB34qI1yT9Q9IykjmMrwe+KWkJsJJkGKikiFgu6SfA05J6gOcj4suSrgdmp7ejdgFXAGv6bD6bZOrMvwEPAOeksWwELuozq9rJJPNMmA0afwlsViWSPgFcFRGXltH2YeC6iFhZ+cgsLzwEZFYlEfE88PdyfggGzPLB3wabrwDMzHLKVwBmZjnlBGBmllNOAGZmOeUEYGaWU04AZmY55QRgZpZTTgBmZjn1f4/W7UaY2orwAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# plotting\n",
"plt.plot(pos, total, 'k--', label='Total')\n",
"plt.plot(pos, dw, label='Double Well')\n",
"plt.plot(pos, wca, label='WCA')\n",
"plt.ylim(-0.5, 6.5)\n",
"plt.ylabel(\"Energy ($k_B T$)\")\n",
"plt.xlabel(\"Distance ($\\sigma$)\")\n",
"plt.legend(loc=1);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Of course, this matches exactly with our expectations. We set $r_0 = 2^{1/6}\\sigma \\approx 1.12\\sigma$, and we see the condensed-state minimum of the double well at $(r_0, 0)$. We set the height of the barrier $h=6k_B T$, and we see the barrier maximum at $(r_0 + w, h) \\approx (1.42, 6.0)$. The extended-state minimum of the double well is at $(r_0 + 2w, 0) \\approx (1.72, 0)$.\n",
"\n",
"For the WCA potential, $\\epsilon$, in units of $k_B T$, is approximately $1.21$. By construction, the inner double well minimum is at the location where the WCA potential switches to $0$. However, note that at smaller values of the distance, the dimer *does* feel both the WCA interaction and the dimer interaction. The dimer interaction dominates at any practical energy, but this does make a slight difference from some implementations, such as Swenson and Bolhuis."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Double Well Dimer with 2 Dimers\n",
"\n",
"Now let's do some actual dynamics with this system. We'll use have 2 separate dimers (similar to Swenson and Bolhuis) with a total of 216 particles. In this example, we'll use the positions from the `TestSystem`, then energy minimize, and then run a little MD."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"dw_dimers = testsystems.DoubleWellDimer_WCAFluid(ndimers=2, nparticles=216, **potential_kwargs)\n",
"integrator = integrators.BAOABIntegrator(temperature, collision, timestep)\n",
"sim_dw_dimers = mm.app.Simulation(dw_dimers.topology, dw_dimers.system, integrator)\n",
"sim_dw_dimers.context.setPositions(dw_dimers.positions)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"sim_dw_dimers.minimizeEnergy()"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"n_steps = 200000\n",
"interval = 1000\n",
"reporters = [\n",
" mm.app.DCDReporter(\"dw_2_dimers.dcd\", interval),\n",
" mm.app.StateDataReporter(sys.stdout, 10*interval, step=True, potentialEnergy=True, \n",
" temperature=True, elapsedTime=True, remainingTime=True,\n",
" totalSteps=n_steps)\n",
"]\n",
"sim_dw_dimers.reporters.extend(reporters)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"#\"Step\",\"Potential Energy (kJ/mole)\",\"Temperature (K)\",\"Elapsed Time (s)\",\"Time Remaining\"\n",
"10000,5.435455322265625,96.46224241638895,0.0005002021789550781,--\n",
"20000,3.506977081298828,96.19163981663785,7.222485065460205,2:10\n",
"30000,1.1279618740081787,97.83939633447592,14.777429103851318,2:05\n",
"40000,2.1191678047180176,105.92029436584274,22.226277351379395,1:58\n",
"50000,0.9582363963127136,92.17251093712606,29.623764991760254,1:51\n",
"60000,1.97712242603302,107.47657367879991,36.93014121055603,1:43\n",
"70000,9.438933372497559,99.1488061429421,44.14107918739319,1:35\n",
"80000,2.5515003204345703,96.42551559307948,51.70799016952515,1:28\n",
"90000,4.100320339202881,90.12446376486484,59.216212034225464,1:21\n",
"100000,2.4359078407287598,97.72176174679997,67.03333830833435,1:14\n",
"110000,3.1212961673736572,100.43919548986091,75.06573009490967,1:07\n",
"120000,2.851691246032715,94.15675063489913,82.89549231529236,1:00\n",
"130000,1.6578476428985596,95.38218833424088,90.46070718765259,0:52\n",
"140000,4.061891555786133,101.59969287007205,98.28142023086548,0:45\n",
"150000,2.02421498298645,91.19585724972109,107.32349109649658,0:38\n",
"160000,6.398697853088379,87.73605908399132,115.76265120506287,0:30\n",
"170000,0.8174424767494202,93.56348555985139,124.14089322090149,0:23\n",
"180000,5.09033203125,93.42113230207633,134.3850440979004,0:15\n",
"190000,3.7300450801849365,100.73669182142467,142.2352740764618,0:07\n",
"200000,0.9542685151100159,97.01255568681678,149.87825226783752,0:00\n"
]
}
],
"source": [
"sim_dw_dimers.step(n_steps)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"traj = md.load(\"./dw_2_dimers.dcd\", top=md.Topology.from_openmm(dw_dimers.topology))"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "09a60ee72b2242ea9c3c132d9ad92c9d",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"NGLWidget(count=200)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# movie of trajectory in NGLView, if available\n",
"# note selection to make the dimer atoms blue\n",
"view = nv.show_mdtraj(traj)\n",
"r0_nm = r0.value_in_unit(u.nanometer)\n",
"view.add_ball_and_stick(\"Ar\", radius=r0_nm)\n",
"view.add_ball_and_stick(\"0-3\", color=\"blue\", radius=r0_nm)\n",
"view"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x152159d3c8>,\n",
" <matplotlib.lines.Line2D at 0x152159d518>]"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsvXecXFd99/8+d3rb2V7Ui23JlizZuIHBYEyxqYbgBBwCocUhyUPyJISUhwSS8EvhR0mjhQeMIQkldDDVGIyNsQHZlossW32llVbbd6f38/xx7pl7Z3dni3ZWml2d9+ul12p37szcuXPv53zO53zPuUJKicFgMBhWF9a53gGDwWAwNB4j7gaDwbAKMeJuMBgMqxAj7gaDwbAKMeJuMBgMqxAj7gaDwbAKMeJuMBgMqxAj7gaDwbAKMeJuMBgMqxDvuXrjzs5OuWnTpnP19gaDwbAieeihh0allF3zbXfOxH3Tpk3s2bPnXL29wWAwrEiEEP0L2c7EMgaDwbAKMeJuMBgMqxAj7gaDwbAKMeJuMBgMqxAj7gaDwbAKMeJuMBgMqxAj7gaDwbAKMeJuWH5OPgwDD53rvTAYzivO2SQmw3nED/8KKmV46w/O9Z4YDOcN56VzPzGe4Zt7T57r3Th/yIxBduJc74XBcF5xXor7F391nD/64l5yxfK53hWHXAIK6XO9F8tDdgJyU+d6LwyG84rzUtxTuRIAI8n8Od4TF//zBvjOO8/1XjQeKZW45xPnek8MhvOK81Lc0wXl2E8ncud4T1wkTkFiFUZFhTSUC1DMQLl4rvfGYDhvmFfchRC3CyGGhRBP1Hk8LoT4thDiUSHEPiHEmxu/m40lU1DOfaiZxL2Yg0LmXO9F43Fn7Tnj3g2Gs8VCnPsdwE1zPP4HwJNSyt3A9cCHhBD+pe/a8pHO2859qonEvZSFYvZc70XjqRH3yWV9q2/uPclXHxpY1vcwGM6YU4/A3e87a283r7hLKe8FxufaBIgJIQQQtbctNWb3lgft3IebKXMv5qC4CgdUa8S9zqDqI/8NmblOsYVxx8+P8Z8PLmipa4Mmefpc78H5w57PwH0fhNLZ0Z1GZO4fAS4GTgGPA38kpazMtqEQ4jYhxB4hxJ6RkZEGvPWZYZz7WcQt7rMNqiZOwTd/Hx794pLfaipTbK4KqGZn4CH40DZ4/CszHxvap+YmNAupYdj7+XO9F0vj9GPqZz55Vt6uEeJ+I7AXWANcBnxECNEy24ZSyk9KKa+UUl7Z1TXvXaKWjabL3MslqJRWv7jP5ty1Y2/AYPJktki+NKuvOP94+vtwx8vVuVWPiaPq53ffBSmX2Ro9BB+/Fp66c3n3cTE8+HH4xu8pkV+JlEsw9KT6/1mqHGuEuL8Z+JpUHAKOAtsb8LrLhq6WaRpxL9miXkir0sHVxHwDqlrwk4NLehspJVNZ49wBdR7d+cdw7D41gaweumHNTcLdf+v8fXCv+jl+ZPn2cbEM/Er9XKkVZaMHoGzHMSvIuR8HXgAghOgBtgFNdFbMJJPXzj2PXEYx/cG+07zh07+gUpnnPYp2IyPLq69ccD7nrgdZE0sT91S+RLkiz764//Cv4Au3nt33nI+f/zskT6n/zyUkmTFAwObnwvCTzt+H96ufS/xOGkalrNYnAhXjrUR0JAPNI+5CiC8ADwDbhBADQoi3CiHeLoR4u73J+4BrhRCPA3cDfy6lHF2+XbYZfEx1H2ejXIKf/OOsg3SViiRTLBP2e8gWyyRyyzf2+8DhMe47OMqxsXkGSkuuOOZsDqoW58/5x1J5nji5hNml2QmI9gCi2h0dTua4/O9+yM8OjjqCP48je/fXH+fOx+pf2JMZ1Sjmimc5lhk5ACebaFE0KeHnH4Fwp/p9rgggMwahNgjEas8DLe7JJhHS4f3OdTFd3MslOHrf3D1eKeGhz0I+tXz76KaQgfS0HtPpx53/N4u4SylvlVL2SSl9Usp1UspPSyk/IaX8hP34KSnli6WUl0opd0op/2v5dxv45h/AXe+Z/bHTj8FP/wme+OqMh3KlMlLC5s4IAMPLGM2MpQsAPDowTwlg0bUPZzN3//rvwtd+Z85N/uPeI7zpM7888/fITkC4AwItVSF/6NgEE5kin/9lP2TtY5M8PecF+rWHT/Lj/fXz1qmsEvdssdzw3lg6P4cBKOdVDjxXtu1m8oQSo+UiMw6FJGx8lvp9TnEfVd+NL1K79IV28c3i3HUkAzPFff+34LMvh4M/rP/8kafh238IT31nefZvOj/6G7jjpbV/G3xUNaTQPOLetKSG69dN6xNgbKaz15UyWtyXc5bqeFplbI+emMf5up37PBOZsoUyr/zIz/jVsaWXDjLytBKbORhLFZjIFM9cMLMTEGqHoCPuj9s9gR/tHyafsj9HOV+3HLJYrpAtlpnIFOq+jXbuQEMHVe87OMJlf/dDTk3WaXRLBUBCegHVXz/9APz7M+Bzr3QatUajxy46L1I/54tlwh3gCzmmopCGiWO1r3WuGfiV2s/4+pniPvK0+vngx+o/X383Z2t9o6EnYOywY1akVM5947PV7ytoQPXsIyVkx+ufuPoEGD044yFdKbOlKwqo3H3mRuMNWcRrLKXEaO+JSSoVyXCyTkNS49znFvcjoykeG5hi3wKikscGJnn9px7kxHid10wNz/t+mYLKss9YMLMTEGqFYLw6oPr4ySmiAS+FUoUTp1wCUicGSNrR2USm/njEZNYR/nwDo5mfPDVCsSw5NlrnfCjb7zufEJYK8JP/T0UgsjL3QOdS0HXrCxL3cVvcw464jzwFSOjcpl6rGcohB34F666CljUz4ztt4I7c48RJ08napqFwdhwzE/1QKTrjTf0/V0b0wher341znwO9Xkk9AdYnwCzirp37Ftu5z1ox87mb4a73Lnx/fvEf8NW3zfizjmWeHEzwD9/dz3P//59U44MaajL3ucX21KTa34WI7Y/2D3P/oTHe8OlfzGxYykV10s/TU0jZkUSm4LrI73oPHLt/3vcHbHFvs8V9CiklT5yc4iU7e1nfHuL08JCzbZ3BsmROHbNZj52N27nnSo0TJN1DGk3X6TXoCojU0OyPa7T4a/e22ElbxdzMzLj/5zO/v+nOfa4lHzJjEG63nbtdqaUF8oIXqAH+hfRIlhMplQvuvliJ+/RGdOwQrLkcvCH45f+d/TV0Q9oIUU2Pzr18dSnv6I8+do/8p4olL/11sLxG3OfE/rKK2TonrhaJqRMzTn7t3NsjfsJ+D+OzXbTjR2Dy+ML3p//+GZmflJKJdIFNHWEKpQqf+tlRcsUKjw/M4rgX4dxPTqjHCwsQ94NDSdrCPgancnz4hwdqH7RPvEx67i6izpuruXOpAPf/K9z3oXnfH3DEPdAC+SlOTmaZyBTZtS7Os7d2UslMVLPI7z3w8Kwvkciq996e/hX84wbVvS7lYd83ql1ft/A3qmImlS+x75T6vkbrzWYuLdC568d7L1U/F+vcv/dn8PnfcH4fPwqfeQk8Nm3yV9W5X6h+1hMSKZ1Yxh9WvYlyQYm7Nwgb7Mz+XFSnVMrw8efA/juVCZFl8EehZa3aH3fcMXZYufq+XTA208wB9cW9XFQNwmI+4xdfr8pM6zE1gJq0jx0dT6nzdOdr1HEOxIy4z4ndzSrn6ox+V78sCeOHax7SNe4Rv+CfvR9h/fBPap9bzEEhtbh1UPIpW3CchiKRLVGqSJ6/vRuAWFDd9GrWwdVFZO6nphbu3J8eSnL15nZ2rYvPqNg5NaCm6Vul7OwNnI127FktmPaxrxz5KS/+h2/NLaTFLJRyNc5dV97sXBunLeInVE4hO7dRQXD48OwXp3bum4qHID+lBvoOfB++/NvVKoRacW9MLLP3+CS6inUsPVPcDwwlkdVYZp5p/Pqc7NmhfmYX6dwn++G0a+0+Pcg4NS2mSA46FTDekDpes1FIKTHXsQwoYzF+FNq3QHyd83pnm9wUDD2uvlttdnxhiPWp3/W1mRpWUUvHBUr86/XkM7bTdvd8CmlVwvrdP1XVRQuhUlHFGnM1BnpiGEB6GJ74mrq+L3+D+psR93mwu7RBmVUHfDqJkyozhBnRjK5x7xjfy42V+7hx8BO1r5GxqzgXM+BVsE8alxsbtcVg97pWXrKzl79/9aVs6Yzw6IlZXncR1TInJzL8oedrBDL2RTf46KwDRflSmf6xDBf1xOhpCc5YauHHD+0DICiKfO+x+ottpaY7d/szWrLExakH5l7Cwe6+7h0THE5atrgn8FiCi/taaA35iJGmHOog6WmjszJGqez6Lu78YzhyT7VctQv72OWTzrGeUvs+6RpszTbIuf/y2DiWUA3zaLK2AXz6dJIX//O95HL29zVd3AuZWhHQIqnFfbGxTD6phFp/11rcp8dBydNKBGFuIUnb57keUAV17hWSqpelX2MxrjafnH0pg8WiBxyLGed68IVULANOFY/O2zu2QiBav9RxunOXUs12PXy3GuzXk7ZA7f+/XzH7+i/JU3bjMkdvd8K1tlFqGE78AqK9sPYZ6m+BFiPuc1Fx15BOjzGkVCfk5usAQWXkAB++6wAn7WoH7dw7jnwDgL7icTjyY+f5ZzKyrk+qjFPer91wR9TPx3/rCl65ew271sXnd+7z1LmXxo7yJ76vsGvoG+oC/eTz4fabZkzLPjKSplyRXNgToy8eZHAqV614SeSK7D/oVBL9YO9R6qFFPZvLqWObdj7jTZ5fzT3L1xb3bz2d5VenK5BPcuD0FFs6IwR9HtrCflpEmrw3xqjooE+MV8dEKJdgz+1w8C4StnPvElrcp5wLxB6Ercncp4v741+B7/9l/f2sw55j41zc18K6tvAM566rrEoF+/NPF/f7PgifvN6JEBKnVNwR3wDCs3jnrj+vrm4a2KN+Tp+OnxyEWK/6/1zirhsXt3MvZNS5HIhCtFvt53w9kuRpVUMO6jh/9a1Otc2Zoq+nYqbWubesVf/XDY7ulXdcOLdznz6guu9r8OQ34Ya/hktvUQapUlbn3N1/pxqN2arItFGcq9pl4hh4/CpbTw2r3zu2ghDq8UDMVMvMRT7hmiNVmNZaZ8bVIFfHBRBfT/rUU/zb3Qer90zNFEr4KRI59G1+Fb6OcdEGD37Ceb5uOOaLZR74mJqdCM5J4xI+XSnTHnFWP969vpWhRJ7xx++Cf7nUcQD1nPvUSTsbHYcP74CBhyhOqYu5L/UEnPilyiNHnoLPvaqmB3JgSO3TNtu550uVanSx59g4sbIzKLS//3Q1W55OOl8mTI4rv/IsePQLVRd0KHwZ11uPMjI5RyNoC8jBhJexYgBkhUxqiu6WAADxsI84abKeGEO00SMmSNljItULoJCuVst0Yv8tn6weu/5jh3mof5zJbJFOT4Y4qZni/sh/qVK5hQ4C2zw5mGDXulY6o342jN0Px35WfSxhH8tqLJOaJoLjR5Sr1oNvyUHlhi1LDWIu1rnrc2VqQJ0velLMXM492FLfZWo3G+l0OfeMup78UbA8avLZfLHMI/+pashTI05Dk17iHEbdILkn2flC0GJ/rqSr1NkTUBGSP1q/Gsbt3ItZ+M6fwtor4do/VIOxhZQS7n1fV/EXqPG66eiegj6mJ36lBrXdTPZD60aIdKlYZvwotG12HjexzNzkEq4R/Omttf7iW9aoQaUxNZDYP6ocQDpf5nprL1Zukj1tL+dO7wvh0F1ODKPddylXK7qgFv4ZUnEGT39XDfiAy7k7PQrt9Doigerfdq1rBWDiyR/B5HEKQ0+pGKImc7c/z/5vwz9fAv0/pzS0HxIDJA8/iJVV+7cusx9OPKgcwo3/CMP7YMQpBTswlMRrCTZ3RuiLq4t30I5QxtNFuoQjyj2hCq/66P38n68/zpf3nKhGMYWhpymVS1xlPY2/MKEExf6M94ZfTEgU8ByfdnK7sS+U45UOhgrqOBTTk3QHK/ClN9BbOE5E5MmICCcr7fSJsWps5u6aJ6c791yi+vgj+57k7779JFOZIh8NfJR/831kZuauK0B++v76+zp91zMFJjNFtnRG6IwGeEPyU/AjZ/0V3VD6pd1jmO5wtdDpgfnEoBMrhNrVcdQN90LQgjB1QuW+laJy3W7nXikrsV+Qc7fP1XC7K3PPqvPPr8qEaembZUZoUa3OqHNqOxYjNeT0epe6fHNV3N2xTFjFGwjXPJbDanzA8oA/Un9tJr0/+aR6bnYcrnobeLywxo5LTu6Bn31YiTKoaDeXgMe+7Lymdu6FpDrWP/obNRHQ/Z4Tx6DNFveJftXot21yHjfiPjfFlCuWmX6g9Bffshbia/GmlbM5ag8oZgolrrAOIj0BTnVcwxPlDWp7fRG6XYfLvQ9NZeF/3gh3/omzXdVd2uLuKhsbn8W571jTgtcSFIbUSfKhL9/N//n649VGRApLncz5JHz3z9STxg7y872qQXnw8f10CPWe4UpKdYN7d8F2ezacy1keGEqxuTOC32vRGw8CzhLHE+kCnS5x/8837uQVu9bwjUdO8q6vPMb1H7iHn+x5DN8nnsVbPN/jWZaesXiqeqH81LqGvPTRPngvdRk9QNnyMSC7GCmpBqacnWSbdRL2f4t1/V9XLyuiHCu2ERcZMimXgAMU0tVqGSeWcZx7R3mUfacSDCcy7JJPsUUMkneXQmbG7QtsMxz9qertgD2lvf74xlG7rn1TZ4TOqJ+2ygRy5KnqhayiIolf37ogPVI7S3W6uCdPOY463K4c/f5vw4cvnj/6qJSduG7qhJO3X3STcoflInz2lcpFy7Irc58j362K+7QBVR3LgHodt3Mv5eE/nqvy6h++W0WXekA3Peyc/zoGkRIe/hwc+tHM9y+XagoQaqheV+5YJghev2q4dOwzelBFHqD2uVKaPSuvinuqNo4CZQB9Ebjnn9TM3Be9DxCq0Xr0i/C1t8Epu4rLXY2TT6hjOHm8doG1iX4l5tFudXMOgHbj3BdMJeWOZdL8/NAo9x20TyxdY9qyBgIteOyLQk9CSefL9HiSiGg3kVCIQ4V2tb3uhrlyc+3mH+qf4A3/9FkYO4hMnuLeAyPI9LBTIWN3ze+4aw/3HlD7MZYuEAt68XudQxz0edi9vpXglDoZKhPH+fojJxmfmiIr/eQIKMG594PVHkhp8iQPPaHEdXLoBB24YpDESVh/NbRuUF3Bo47QHhhKclFPDMARdzsnnsgUapx7m7fEh197GU/8zY185e3PojXs4+s/uAshy/yG5x6eadm9leRpdXyCcU7n/fyisp0NEw/W/6JGDzERXE8ZD0mUgFj5BL0edXLHhlVuPFIMcrysLrbyuC2G05x7q79CXGScx+zH+6wJShVJe+44IZmlR4yTK7jq4bVrv8GO0I4/qETjP56nlpGtg64u2twZoTts0SpSiHyien4lsiVCXrCEZFS0q1JCd014Wot7vz0ONOjECiE7ljn1iOoh6vy8Hm4xmDyhBuniG1TDXimp1zn6U2c5jqpzb6mf72bGVK8v0DItlkk6zr1tk4oVdKN15B4lgNtfrn6f6HfMVGrEMUaZMfWcb71D/XvgozPf//t/Dp//9bk/7/QBVYDuS9QM0EJGia0eoNb7PL0nXy46FUP5pNPwaHG3PNC3W13/666C3a9TcdTUgD2hCzjyU/Vz9BAI+3rOJZzXOmyP2aXHlCFs3aheQ5s+E8ssguw4RelR/y+k+JcfHVQOGNTJJuy80B/FX8liUWE4mSdTKJEplOiyEhDpJBb0cqxkf8lV5+66QO1B1ScHE9xkKbckk0O8+fYH1MVZKTkXMeDPj1en1o+lC3S4XDtH74Uvvp7nbY3TW1ICsVaMUixLfvz4cXL4SVX8VAopOPpTKpueS8bfwQN7HyeUV/vUwSRdYoo0QbLCPtnXXaV+brpO1dtXKoyl8vSPZdixVi2r3x0LIIQTy0xkCvRYCSUyUHWFliW4clM7L7y4h66cilQusk5ymaUao8LESZ46cgzCHSRyRe6t7KKv0A8/+2f49I0zyzhHD3BCqEGwKakmjcVJqeMP+IZUlUJ/xs+AVN1hOWF/DznHvSVyRS5pcTmyfJJMQl1Y672qAd4p1KCwX5RrF23S66RseKYqDUwNKSdfzpM7/hB3PTn75KOjI2ne5f0Sm07eSZ/PEYzkcbW6XyJXpMNO3E4LVe5azd2LOWdAfvK4cunlPMTsWCbcpoRBl80NPjrrPrg/b5XJfrU2zebrlDsEp8em37OmWmYOcQ93qIE+7dwz46qR8qvvip4dar+1M933DVXS+uz/7exLwo5l3M49Mw4Hf6B6Et7g7FUsY4fqzyitEXfXgCpA70615MDgXrWvvbvU3/U+T8/d9ZhHuFOd53ofw23ONrqS5cZ/VMcjvlaJ+6g9N+TIPaqRmTqhGhewnfu483hyCP77NYCA9dc48Q5Mc+4tZ+1m8StS3D25CQalEqZKPkUqX+LEeFZNs0+cUs7F8qiTG4igRO3YaIZ0oUynSECki5aQj3FiSG+I0YFDjCTztjDYI9t2LHNiPMNLPKo7b5XzbBDDCD1RwVVn3C6SyLHD8PT3GE/nnUimXILvvBOeupOX+x8mJJTT3+Ib55rN7ZQLWQrCT1b6SaWSkBxi3NvN4VwMEoNc06m27xZTdIgEU552jvjsSSrrr1Y/N1+nTuThfezpVyf01ZvUMfJ5LLqiAYaqmXuBDjHlZIHTRLkt7GOzHKDkjZCXqj7/VOQSrPRpRodOUg51MJVV4g6o7PHEg7VzCspFmDjK/lIvsaCXYdlqf4ZJ2qQSIWH3eI4mvQxItYqhlbB7UFXnrgZUN/hdF20uQT6tvptAOc0lHRaXWk7Fjyc9qG5pdtd7lPsKxFVMF+1SF7cdmST793Lbf+6pZvpuTo6M8bveO/Hu+zK9HkcgP/Klb/H06SRT2SIdIXUOnMQWWR2vuA3C5HFXVKhjmQ4lDOOLFHd/TLn07Lhapjfao/6uxd2remhV5x60Y5lZc+gxx736beHURsW+bqpCNrxP9VCf/g5sexl0XqD+PrTPaVBS02IZXRK49oqZRQ/6M6VHZy9lrjegCkrMywWn5LJPi3sd567jp7aN6qc2cfqzg2qsfutrsN42SvF1qoemM/bjD9o32pBOQ5AYVOMelleJ+6deoBqd130eNlzjNLyBuLNgGDjH9iy49xUp7r7CZNXpFTOJal3zA4fHlMtotXN0Ozvc0ak+5rGxNJl8iTYSEO6kJegFBMXYOh594nH+5UcHVOwQX6+erwdZT+3lYus4wzHVBdwuXLNXE06NeIeY4poTn4Yvv5mxZJ7NwTQ8+iX4xcerLmDTUTWjMCHDbPVP8NbnbCZiFYhEomQIkkxMQXqYSaud07KdqzpyXN6mRHmtd4o13hQpTyv3+58NG6519nXTc9TPb/w+4fvfj98ruHRdvLpvl0UneO3hd0FuikQ6S6tMOOI+rZy0LeLnAusUE9EL+HHlcnLSxxMtz8Mri2y1Bsn5WskUyhyQ63hUXoBcc7l6YvI0J8YzvO32n5EZOgSVEg+nu7h6UzujxKlg0S0maKnUTt8+mPAwSpy89OFP2Y2l7dxlIUMyV6LP44qj8gk8hSQFu/f2vL4Sl1pHKXmUSPnTgwzd+yk1k3b/ndC9XTmyaI9y7rYIdxROEpQ5DgzNFJ/Q6YfwUoapATpxxl4uEgMcHU2RyBZps537EewSPS0GOm/3BpWY6Nw65hpQLeedRa/cddazoYWg+2LlVkH11LS4n/iFcoq7X6eqRyJaWOx1bGYrEXSLu3bF+m5MWii7tqkYYuhJJWC5KbjkZnuSVLy2UiR52nHJmXH1mT0BdX7O5tzzSTU+MFtV2lzOvWen+vnEVyDY6pz/dcXddtf6XJ/od+IoTbRLLbegaVmnGt7Uadj4HPVdffdP1WPrn6l+6qqaLc+3J4QV4c3fdca/9HfTttEpgwQj7vMRKk054p5LVmuxf35oRHXDbcdR8asD+ex1ykEfG0uTzhdprUxCpJOWoA+AdKiP7soQ/WMZ5SbsQZrHDh1DHn+Qd5x8J0OylQe6XwvAdssl7rZzT8gwPZ4UvZmnoZSlmB7n5dlvwNdvUyWTa6+Anp1Y/cpl3V/ZQVd5mBfv6OWm7XEikSh54UckBqBSYky0clq24884YtQqp9gZz5LytvFN/8vgLd8jV6qom4HE16nSrnyS607dznVrLAJeT3U3Xyx+yTOyD8LhnyD1RVxP3MN+toqTDPo28J7im/jNwrs55VFusE+MM4W6MNa2hrk5/7ekbr5dPTFximMP3cVH+l9F6qEvAXCg3MfVm9sp42HKaqWbSaLFWnE/lPQisTgpOwik1fGs2OKeTiVI5Ip0W+piGLM6IZ/AV0zSj3LCV7Vn2CGOMbnmOgAC2SGiSdsVp04rUQQleq6SPQvJNjHA06drLzQpJesT9iDa1ABtUgnQ0UoPF4kTjKQKJHKlqriflq1KEPQNGbQD7ttdx7nbcVgpq7JzV4MzK1oIemwn3XGBig60OyykVK33i/8e3vpDNfAI9YWkUlHXiY4LtCvW+60jDl9IVaMMP6lqwwNx2Pp89VjbRtdSvMKu1LJ7CJkxu2qnRxmses4dZl+7Rj9WmCVz77hANRq5KeXaq/Xj0dk/q3burdq596vGyS2404mvVa4c4Io3qcbg1MPw3D9T8R44g7qX3aoGYX/nblVWqdGxjDuSASPuc1LKE6hkOWl348vZRHWK/MHDB9WXbl/MU2V19V3QCl2xAMdG08h8Cj9FiHRVlwQY9vSwVoxyYqJW3O96+ACZb/wxU5Uwryn8LQfYBMBOj6sG1r5whzy9dDPBmqISfn92mN7KsCrfesF74OaPqa40UPKGeVxchL+UhNwU3nIe4Qth+SPEs6onMCRbGaYdKzehBtE8foQsE0z0k/a2kS9VKJQqXPtPP+brj9hu98XvI3/dnwPwzL7ar3Zbxc5N++/Hk9Xibp/w02KZLm+aLpHgCGsZoY19nu0M2TEYwHBFXUgX9kQBwVDF7iEkTxMde4ygKNL5iCqVOyL7qj2IIalq2UPFcWdhK2AKJSYnZSeRrDqeuZRqAKyScu66UmbA6oNcAn85zXFLubZrrX1ERY7YpS+lhEXXovkuAAAgAElEQVRXcj8RmXYctY4Xot1KwFy14bt8Azx9ujaXHksXuFzag8iFFNGUaih+XtnJheIkY4kMyWyRtoASs2zZHpTT8Yp+/XVXKWHrv1854Kgdl4ScY8mOm9XPU3O4dx1R6c9hn0fVJQZAVX0EorDmMud52p1Oz92Hn1TXiV7ETMc5usehhVK/56lHVGXPJa8Er92itW10TEHXNnXTEgCEcvC6rt9fR9z1mMps90StRnKuWEZ/To/XaeT6djvPqWbu05y7HvTU5/pEf20kMxt66QVQgv28v4Cb3g83vNs5plrco73w7D+sfQ44DW9bE4u7EOJ2IcSwEOKJOo+/Swix1/73hBCiLIRon23bhmB3s3Q3vpxLkS6U6IoF6MzYma99EQznlTPvCxbZ3BHh2GgGf15P3ugiZjv3/lIH7SJFZtJeqyLWS84K0SZS+CcO893y1QzILgaK6qS/2HLEXdqxzJivj6DM4kF1m7vkGG1lu6G47p0qGrAvSk/Xhbzj125QLzB5Qg3A+cL4Q1EiUp2cJ0txUgG79S9lnaqASpGMr41CqUI6X2I8XeDwiHPxHE6pz/SM7tqcdV1eRQblY/fTmz+m/thh5/bTZsXqwdQni0qMulsCDFac3PBkXnWRL+xWx+N0SiqnkjxFMK2Oh1Upkgt2kSJMS9BHNODlRClOj5jAlx1VDqxlLUV85FH7fFp0E8upCCOXVOLur+RIZIu0ywnSnhZGKzFIDeOhzOnAJgCCD98OlpfAhdczSjubE8pR/n3hdVSe80645FVqx6PdyPQogwOHGaOVrAhzTWSQp6Y59/7TY+wWh8lG1AXrOfUwaYLslVsJiiKV8SNMZYvEfeoYZ8pe5SJHD6oIQveMdD772Jdgx6sdRx12XR4XvxIQc+fuWgjWXaXE8uJXqt+FUJEC1DSWVbQQpUfgqe86VS86Ttl4rfM6vrDjou0eL6DOu8RJJdC7Xuv83V273bfbWR2zbaMdy5xW2X8gqjJyd9mju7RzNueuG4NiWv3z+JWoV/fJjmZ63eIerX2uZrpzTw7WNq6z0WILteVTn+d574Jn2jeeC2pxt2MZd57uJr5eNSLa6WuqDW4TiDtwB3BTvQellB+QUl4mpbwM+Evgp1LKBtxJog72l5WyWkgRpJRNIiU8e2sHFwlbdG3nftoW9y5/kU2dYQ4OJ/HlHHFvCakT5umsGuzbLm13G+4kISNcJE7gk3mOyV4sASeyAUp4WSMd5ycnlWueDNa23D1ignhh2JkyDepiEhai40JCnXaLPnVCibc3SCjiXFTH81EKoV7nuX2OI8v52ymUKmTssYaEa0Bw/6T6zNtbXDXXuQRt2X4mZQRr+Enean2HqdB6e+U/MaPeuzWtnOqetBKOrmiAwbKT3x/NKqd3oV1qOZTIqQs5MUgoPUB/pZuSL0YypnpAQZ9FW8THsGyjR0xiZUZVY9C1nZw3BgiiAS/jvh6ipXEoZinYA6ZeyohKkXh5koyvnbFSsJphF0JdqnqjkIRr3g5tmxi1Omgrqu9nX2UTiWv/QsUDANFuBJKRw3sZpo1i58VcbB3n6aFkzc1IEoceICBK5C+5Rf1h8FEmRSsTMXXf9/j44yRyReJ+1ZDnpIdKzy5AqjK99LDaL914IuB5f+4cYLdz7NmhDMDQNO/knv6uhaBjK/zlgBONgJPt6lUg3Wghuu9D8MVb1VLWydOqJxFf74xNgRL3es4d1HmsnT44Yhnprj3HO7cpt5wcsp27fU67RdctbLPNZq2uAVNRPQwdyWh078TdS5krc/eFXdUrsrZxnQ3twtu3gMdX+5g3oGIh7dzrvVYgCn92BLa9ZNrf6/SmloGF3GbvXmChYn0r8IUl7dF82N0sT6SDjAxSzNrT7Htb2GYNkPJ3VQ/4yYzKnNt9BV5yaR8TmSJee4anKoVUX9zDSXXAdwkl7mlfG2PlELvtEsCjspcLu2NMZEuMW6qlLnjVSSvt+vhUWEUEKamEr49xgtkhld9pgnF46Qfgmb8HrfZA0ORx27kHiUadQZ4juQhlXdIGNV3QnL+dfKlM1o6jprKOkB9MqAYrUkmqEf17P1idTPH58gsQSC6x+unf/Fo1Fd4fmRHLhKYOk5c+9ibjWALaIwGmCjAq1f49nVRdc11HP5TMqcHC5CCx7Emekht46Dn/wSMXvwuAgNdDe9jPkGxTk6e0uF/zu/yy7/UAtIZ9TPjsxmxqgHLWGUANkSNWGicf7GSqElQDcYAViiuRivbC9X8BwJhHxXVpGeA07Uxmitz+s6Pc9rk91YHG7dYAF23dSsvGy1ibP8JkpkBiz//Ah7ZDMUvplCqrjV3+GrUD+QTBtj5+4+U3kRYR1kw+TLEsabGdewEvhW57Kd/Bx1QsE+2xxVOo9Uu6trkOsC0I0R51/Nu31K7HcuAH8C871WuBq1omOjMrnkvcdQRw9F7VoJx6WE12OvYzx7VrfGEnwtARB6jSQ1BL1louudDO3Z39A3RdpGr381PO53N/hun/nytzB2Xm9GCq5rLfgjd+s/Yz183c7RuSuBus+cQ90qVc+2zHFFSjqUsu6zn3ejRTLLNQhBBhlMOfeePSBlJOqxMw0NJFiiAV+yB1Rv1c4jnJKf+m6rbH00rcQ5UM11/UxeUbWmkX9kGNdBHxe7AE7EsrV7rbUrHOgWSABBFiQjnao5Vedq2LM54uMCrVtrnIWvLShyejTs5CixLrfXITaU+cHdYxhCzVuhpQ057XXalOIF1NUcqCN0TEFveCJ8yJlIXlfq4u+QIKgU4KpUp1DZWEa7nbp6bs7mt2Qt1f8sfvU2uBo8S9KHzkpY/JbfYEEl94RixjjR2gX6yhgkUk4CUS8DCUyHHazt0HcspJ9bYEaQl6VQlqrBcSp4jlBzkhuxhqvZzhsLo4Aj6L1rCfIewLQVaUIFx0I09vfROgBnGnAnZjNnkc6VoTJUyecHGUUqiLpHQudG84Dq/4N/itr1Yvmklb3A/LNYBgKlvkgSNj/PDJIQbL6vj6KeKJ9ULPDgKlJGsYI/PU3apHMHmc0sRxcgTw9O5QFznQ0b2WF+9cy9HwpWzLK/F3xN1HIdyjvtPBR1UsE+lWd6B6w9dUg+5GC0L7FvWzbZMSd917+NWn1E89IzKfdNZ7mU6sT51H2km70UJSKcHOW+A3v6TKVTOjznrtGrc79ruEsH0L3PIZeO6f1m6vxb1lreOKhcf5THrfArPEJW7XWlfc7UYsPTbTufuCsOX62r95g2pcY7pzTw2p9XPc1THzxTKWpXJ0vUzvdPRr+SLOGMRCWYniDrwCuH+uSEYIcZsQYo8QYs/IyJnd4SXZvpN3FW/D27GRDEGkXWYV8Qu2coKDcn1126Gc3aXKJxGlPH9+XQcdevGpSCdCCGJBHyPEKeDleZZySo+Ne6qTbrLSTyXay9q2EFPZYjWeqIQ7SaBOurz0IexKiH2VTYzQxjMse5Gh6eLuHAzV/Zs6UXXugZB6z4SnndFUgVhrm3OhtW+tdnELwXbypUp1INkdyzw+BhWEPe3ejo9GnoJoL6XYer7F9dxevolYm+34fKGZk49GnmbAu94+rl7Cfi8TmSKnpRKlcdR+tIS8XLu1k7v3D1OJ9UFmFH8lxwnZTb5Yrq45H/B6aI/4GbJr3fXxByXqoJx7KuiIu6fonPxhkSeYH8fX0k0K50IPRFph3RWOuwSmfMpFKnFXa8CMpVQe/NnHXJ8z2l3Nbrdbx7FG7MlOkyfwpwdJ+LuVmOqel+1OT7dewRZxii4mifnU8S9IL8WyVL2rE79QjYR2s1tvmOnuPF71N7e45xOqQZ4acKbr6yqbfMIRhek8+4/g1i/OLvxuQdtyvRrzeekH1eCkO9qB+uIOsPPXVK/TTXw9IOz5A/ZnDXeoyUKaWK8Ty+TrxTK2DlTKau2e8aPqcR1dzebcZ0MI9V7TM/epAWdhMc18zh1UEcS2Omm0jrsW8jrT8UcAseLE/XXME8lIKT8ppbxSSnllV1fXXJvWZdTby5fL19Pe3klaBrHslrqjcIoABfbm11S3TZYsiti3tbrvg1zz/VfwaxdaVAIt1RY3FlRleF+K/TYPWM/gwdaXc89ITLku4LjoY31HtDohaaiivlgR6SJhNwApgnjbNlKIruWeym6OF+POOigtzv7MIL5eZau2c9dd2KFKnEK5Qlc0YDuzkLq47IuoFOqiVJHVElC9iNVkpsBYpkTB26KEIjVkOz4f9O3mgu4o78y+mfeXbq2KKv5IbSlkMQuTxxkJKCcYCXgI+5VwDGtxl2qNnJDPw0t39TGczNNfcC7+E7KLQrlSXeMl4LVoDfuqz1cvrL7/1pBqgNvCfvLhHkp4YLIffzFJUtiNCBm8pTSR1m6SLnEPxWZeXAm/et1DlbXVY6Nvd/jf+1wLwcV6q1UXzwydpDWpXHJu9CjtpWHy4TXOdwTV+GOyR010ucp6iojHFne8FMoVdRu1sYNq5qk7qpiNW2533HC1BvuYWsFSVlT5nV63PJ+sL+6t62cKtUY/R1iwyc7Lr3yzyu3dA6LgxCcevzPwOxe+ILzq43D1bU5dfaSrdjyhrnO3hc0XccR9wF646/GvqMd13JQZnenc6+GPqPc5+ZBasVFKZZ7iG9Rn8tgue75qmfnQjWaode7tZkMINR/BHdMtEw0RdyFEHHge8M1GvN5cTNg3ZVjTGiRNEKukxL09pSKVX2Z6Kdo3fMgUymRFWH3hY4cQ6WG2TtyP5ZoarGvdf7nmDXyk+294d/l3uP/oJLFW5UAC3RfyxmdtrIrhCOoL9bR0V9dLScsgsZY4I297iHsru2sqS2aUSLlpXV/j3PVJfCyvLsquWEA1DrEeZxKOsKgE1T5oUdcLa+nFrirBVlvch9Ug3Ov/B170d1zQ7biXNj171heuFfexQ4BkMqwGfKMBLxFb3B+RFzDoWUOSEPGQDyEEN2zvxu+1+PmwM/CknHulujpjwGtVM/cqtiDEw+p5rWEfkWCAQdEN40cJVtIkfeo76BWqMxhr6yKNkwdHWmaKezKgRPmgVOI+mS1Wb5GXrATIYpf9RbuV+LVt4kbfIwSkEv7JwaP0iTGsVvt709+ffc7I3t2kZYBrrP1EvOrzFfCp2x7ueq1zE+T5xH3rDbXOHZS4P/FV5bBbNzornM4l7nNheVTjvvaKWuftrjzRaAGd7trn4rJb1WxV/VkjnbVuNtYzexWLjmXatzjiftRev2XocUA6g+DZiYWLu75hx/f/Ut3oJTuhzm39HepjOF8sMx/6WJ7p67z6E2ocZplZSCnkF4AHgG1CiAEhxFuFEG8XQrzdtdmrgR9KKee+00QD0DfBWNMaIk0IX1kJU0vyABLBgcoaTtk35sgUyuSssLo4dCXA1ImadR90rXtPLMD69jCH7ZtcXLxZVRJsumgXN1+2turcR+xoIRDvIWHnv2lCtIb91bVkTuts2ROY2yXEN9irCeaVO7e7nzq+6IoGVPb3/HfbO9sL4U78PrXP+u5DOpbR4m7pVQdTQ0pEt94A3dur4u6xhD07FzX13B3L2LMm03FV6RL2ewkH1LZfLl/Pv+34H0DQYjvuaMDL9Rd18b1+Z6DvpOwkX1LO3e+1EELQFvEzToyysOODqnPXsYyfSMDDcdmDHD9MhCy5oBKNtUINgnsiHQRjjltqaZ05mDUY3cFbC+/kR5UrABhO5EgXylzUoz572m9/H9oZ9uxkY+6p6vMLI4foZpJIt51ha2GwBayzJcoTcjOXWP1ELOXci3iVoRACXvGv6j6p66+ZsW910Xn5wB41k3nrDapRrzr3OWKZ+djxKrjyrfNvpwU0sAhx1wRbVU8j0ukInjeo/l4d6JzFubdvdqpl9OJcuiRUfz+ysrBYBpxlfyeO2T2oY+rvunhB78uZxCludCyz2MHUs8wsTXgtUspbF7DNHaiSyWVnbWuINz97E5s7IxyRAQIVJeSRqYPkYxvI5oIcH8+wsSNCrlgm7wnbtceuyRIRJxfUFTM9LcFq/PDqy9fS2m6XotkTmqrO3R5Q9cW6SQsnlmkL+wn6VIQxVLZPnpY1c8+EqylFC1ZPYt2AdMYC0PNCZ5vn/DFceguBSdUmT9rOXQ+uHh1N47EE/liHqnxIDUPPpdWna3FvCyvXrd43rBqCR7+o8slyARCUWrcAJ9WAqt/Jc3Vte7VxAJ63rYsPPRmBICSsVjIEyZfK5IsVAvaqmG1hPxKLXKCLSH6kemHoRrMj4idfKnO00s219n1tS5EeSMAmnz2jNdRGLJ6FHFSkIN468yIN+j18yxZ2UHekAnj15ev42E8O4Yl1w9hJl7jvgKfUuvyTrTtoHX0YS0haeuxS1apzt8U9GmBQRlgvRgh7nMy9ek/bljXwdmfp5QURiKrG7vEvq983XKsW1ep/QP2eTzrrxSyWm2dZkXE2tIAuxrlrLAsuvFEtgaEFL9Zr5+CzOHc9WN6xVQ36Z8ZhwF6KuTo5yNXzWXAsE1Uxjh5rOm4fPx2t6QZyybGM7dyX2kgsMytuhurOtXHe+4odrImHyBAkJJW4B8efRnap+vbj48qJZgplCp6Icj6pYeUuoMa561r37pYAO9bGCfk8/MHzL1CuA9RAJo4I6WUPiK8jb5dDpmWIVjte6Ii6Bg7nimTAcRRgZ+7qAht2O3c3fbtg+8uqywi7bwqdyBY5MppmfVtIOffMmPrMrovkwqq4uzJVn+3c93wGfvL3cPgn0LaRlpi96FrAQ8ivjpEQsFWLe8iJYVpDypVXLD9DHiWayrlXqksgtEXU9oVQj2pc7bK63niQf33dZbzqsrVE/F6OlHsQ9tRvT4sStNdvtxuXUBvtbeqCShGkIxqccUiDPqch6owGqhO8LuqJ8sh7XkRbd60T15PDTshujlrraSmqmKAay1zwQlVp0qsayc6YnzRBQuQJWjpzt2OZpdC2SdXHe4NqVqReS71SsWOZlnlfYkmcSSzj5tbPw5VvUdm2+x6s+vVmlEIKp8ey/9vKVLjr6KOuxmzBzj3qrNcDzlK81fVnGhXLrAznvuLEXWNZgpwVxkeJCFk8E4cJrN2JzyM4Me7EMiVvRAldfsrJQ13OvcXl3G/c0csj73kRmzsjcOGL4Lo/VWWLUBXvfXIzmTf+ALZcT9GnTpYUwerj7ZFAtWSwbqWMJu4Sd7dzJ47XEsRDvlmfpt3wlOu+oVPZIkdH0mrfQ20qfpJlx6ECHdEAbWGfk7eDalCKGTUIKCvKQXVuo83+PG7nHg/56LFvkecW93DAAwjyLRvpt5cEyBdVLKP39dK1cV62q4/g2h0zZlPefNla4mEfkYCXfuk0Rv42dfyspL28Qridzk7VuOqZr9PR4t4S9NIR8Vejqo5oAK/HUt9JoMURHbtiZji0hYcnXcKmv5v4Orjl09WGtyMSICsDREQOr30XpgJ2tcxS0Ln72iuVQLasUeubZMbOPHNfDD57LONMYpnpxPqcXqk3oEzV9AHVQMw5N/d8Wg36X/Z6Z5szce6BaO340bH71TWlHXZATZg7o4HQmvfR4t7czn3eWKaZKXpCIGGXdQQhy4ieS1jXFlZ110C2UKIUi8K4neNddKMSdtessWrm3qJcYNX5hdvhBX9d3S7o8xDxeyhWJKHN14AQlPwtkIe8Fa461M6In8f0wOFclTKgLgLhUSLsDcGGZyKv+h32/eISOoMBLGv2SEcL5mS2VtyPjaW5Zku72ne9euC0gb2bL1urBmo1vohz31kEIKHroqq7jwaczL0t7K/2JtzCGrUff/z6T/Pxu1WclS+pUsiAT+1rLOjjo7/5DCj+c3US0nSiAQ/90mmMop22wOpbuYXa6Omyp/yLyKzHR39/rWE/8ZCvGpdU19Z/zh/Djl9z4rK2zdCyjnTHszhwIAX6Y9VpmP1ei5I3TIR89a4/BbyNce7gTC7Szjdx8iyJ+xKdu5vX/pfjbnU0Mz1zD8ScHvTgo3DtO2pKWt2mZFHVMhrLp6rQOi+qXVws1Dp72ehiWEop5FlkhYt7BEpwlceuKe++hHVtCU5MZJBSkimWKfsizj1KY33wyn+veY1NHRFagl56W2Z28afTFvFTKksnrw60QBLKXuekao/4GRNxys/6X3h2vmbuF/R4lYhMHVfOPRhHvOyDdB+6D88cfSrdkOgBVYAjo2kyhTKbOiLgcXUX3RcJ8Dev3FH7Yr6Qsy7IZa+Hvf8F3ZdU3X3Y76k699awj7awn7DfQ7ergdBjFePebobLY0CWQqliZ+7TLiRf/eMcCXg54XLusS5bYFNDqhEMtLC+W1CUHnKuY+4maDcmrWFfTe+iI2qLe6y3Nr+2LPijvVxRhPiD34Z7/q9yZP45ogB/hHAhVxX3Il4K5dkbrAVTFXd7cpFuXOzqpRUl7l3T1rkJxGonF+Wn1N86tqqZzc98e3VF0yo14r6IWAZUIUPfbtULdfeOd75m9jV4FktgZcQyK1rcS54wlOBKzyHVUrdvpT2yj+PjGfKlClJCxec6WSMza+tffflabtzZS8g/f2vebou7RtglUWXXe7zk0l5iQR+eG1++sA/Rul6Ju9dxJ//r+Rcgqd/N98/i3PfZd4Ba3x6CQn1xn/liLpG8+m2qumLzc2kbV68dDXirx6Y15MOyBF/7/WtZ0+rsr3bu6Xy5Wv6oq2W02C6ESMBLHj+Dsp0+MY433K4y6FKuukzrho4IKUIUvbOLXdDrREg61gr7PYT9c5zqHh8RD+zeuQvuoXbJiFkIhFugAOQmkQhKeCiUlhjLXPxK1YPaZK/4qJcH1rd6W/bM3RbQRsQy0/FHa++QpMcQwu3wTtfdmIIt6hpNj9h3ibLsaplFDKiCuqa6ttni7hr32v4y9W+pdF6o9Ka6dlBzsqLFPe9TschzeQg6LwGvn2jASypXqq67It0n6yxCZ1miKk7z8crda6i4FpjyhO3szvUeN2zv4Ybt8wiqG+0sXI72Zbv66mys8Lsy92jASypf4olTqgJhfVsYEm5xn6fe2u2K2rdU16Re0+phY0eYi/taiPidWAZge2+t0GjhzBRK1SURqrGMd+Hirr+HftlDnxhXF7s/4og7EPJ7OG21UgjO3iV2xzJ6HKRz+sB0PbSou93eLLz0igvgbiA7gfT4AaEmMS2FYIsqe9VEupW46QHC5Xbu/iVUy8z72pGZsUywTu6ta98DMRUZFpILd+76Omzd4Dj01rm/yzOi+2J49+DMRcWajBUt7odCu/iLibexOzrJrS9WFZuxoI9krlRdMVG4L4pZnPtieNt1W2p+90XaZr7HYtEnn3eB7gTXgGq2yIaOMKmREvsHlbivawtD0VVrPN++6Ys63FEz0SXs9/LTd6mZj8NJNcGnNTz7zEUtyql82SXuqlqmtc6g8GzoeKe/0sMzrf32zZsjwFhNvpn/tTtY19o562tUYxmXc69GMvPhC0HXxbXrhM9CrMUWpuyEmtEJFJeauU/H41Vm5NDd6vf5xm+WSrUUcva4a0lMv2FHPlm/AW3fota29/rtpTGSi8/cWzc64h7fUH/7pdDkwg4rXNx9gSBfLN/AUy2t3HqBKqOKBdVU8GoeXZ2V1rawadWLef+YLTBLGX3XVQVzZNHT0c69VJHEQz5CPg+ZQpnOaEBFKHp/ot1z19mDc1G3b6m7STTgxRLUDsS6CPoshIBkrlitGskXK+SLZQJ1nlPvfQB+Li7j1zcUsHxOeag739x+6VV1XyPkc8YHquIeWcTiTr97r1MyWw8tIplx27mzdOc+G7oc8rp3zlzoq9FUJzEtQw/BH1VLAGvmGiC+6m3OTa/1Pi04c7dfs3WDGpjeecvMBcbOI1a0uOt8NRJw8nItEMP2lHNP0P7C58uezwBP9zZ+t/C/ubjvBfNvXI/tL4exw8oxLhD3IGXI56El5CVbLKu8HRwhXMhn1hfO9DvGuAj7vXzuLddw6dr4rI8LIYj4vdXZw6BimUKpQsC38MqEiP3djWx4CdZb/rZ2/xZYdqZjmXjI54plFtGoL2hdFXufsuPVNYqWXC0zG9f9ieodPOONjX/t6TRyQHU6gdjMSUzTFyLTrLuyWn5cc7u/haC3b9uoYq5bPn1m+7tKWNnibnfj3YNlWtxHtLiH7Hx4iZHMbLSGffygcjVXR5dwQYTb4UV/u6in+F05dsjnIR7yMZTIq7wd7DxTLEzc/fM7d4DnXDh7DKKJBDyMptzirmbNLiZzjwVVTf3zt7nGCfQFu8DKhEC1WsZfrZZZcCyzULQAZicQtlssLodzv/gVjX/NeixnLOOPKrf+wMdU1UwxvbAeQtW5L1Dcu7apsYo1zzjzfV1FrGhx113wsKvSRdetzxD3ZXDu2hEuyhk2ALdghvyeas35hnb7ArUsldG2LiBv1JNX5hH3+Yj4vYyl89XfnRmqCxf3gNfDj//0+toBUC064YWJuz4WnVGXuC8mllkIWgCzE9Wp7PnlcO5nk67t6jZ+7ps8Nwqduf/sw85aMgsS93Dtz/no2ArvOnhm+7gKWdHirgfPapy7Le5DCTUI6A/b3b/5qkbOgK1dUT72+mdww/bGv/ZcTHfuWsSqsQzAG78FkQWsodG3W3X7L3jh/NvOQSTgxDJeS6g699Isde7z0DN9vsEsmftc7FjTwmfedBXXXdjFWCqPzyPsG3k3EC3usoJYzljmbBJuh7f9aHle2x9RNwxx35hjUeK+8GIDg8OKFnft3N0LW8UCSuiGE8pFLqe4CyF46aVzly0uB243HPZ7KFVscW9zOZzOCxb2Yv7wjIldZ0LY76lO9VczQ8tq+YFF1LnPiu5ZLDBzF0LwfLux7W4JsuevXlSzyFlDcOXSwuPH5xHLE8usFtw33L7wRjj4g4XV7S92QNVQw4pdWwZmj2W0c9fle/62tbDl+WqN7FWC3zV9Nej3VMVrffu5uwh0vT2odWcyhTLFslxULLLTaIcAACAASURBVDMri3Tu09HrzjcUdy7tCeDzWCvfuS8nuv688yJ1y8G1V9TcNrIuix1QNdSwop17dUDVNQkpFqytlgmFw/DGb5z9nVtGhBD4PRaFcoWQz0NnJEBL0EtffOHllI3G/R3EQz6OjSkXH1xEtcys+JYm7suCL0R1HR6vH7/XWp5SyNWC7ulseb6qZPmdHy/seca5L4mVLe7embHM9FLI0FLFpUkJ2IIS9nt447UbufmyNWrVw3OE+zuIh3zVez0v3bnb7q2ZFmkSwrmlm8eP32OZWGYu9PrpW29Y3PMWWy1jqGFlxzK2oIRcA6oBr4XPowb0Al4LT52VFVc6elA15PMQ8HroXsDCZ8tJxOXcdX05sOgB1Rm0bValnctQyrokXPcc9XmslV8ts5xsfDa8/qtqVdbFEN+gbrjtPbfn9kplIbfZu10IMSyEeGKOba4XQuwVQuwTQvy0sbtYn9kGVIUQ1bsrhRewGNhKRTvi0FwLYp1Fpjt3zZKd+6W3wDufaj73psXdG1C9KCPu9bEsuPCF88+Wns6Vb4F3PLT45xmAhTn3O4Cb6j0ohGgFPga8Ukq5A/j1xuza/FRLIact/KWjmTlXAlzhuJ17MxCZlrlrllwtI0TzCTvMcO4mllkGPN6l31jjPGbeK09KeS8wPscmvwl8TUp53N5+eI5tG4peyKp92oJWWtwXs9zsSkPHHSF/c3zG6QOqmiXHMs1Kde1we0DVOHdDk9EIZbgIaBNC3COEeEgIUXchDCHEbUKIPUKIPSMjI/U2WzDXbG7ni7c9k51ra2tmdTnk+eHcm+MzumOZlkbGMs2KK5Yx1TKGZqQRV54XuAJ4GXAj8NdCiFlvdyKl/KSU8kop5ZVdXUsfIBNC8MwtHTPqmHXd90JuwLFSqYp7k3xGHct4pq2Pv+rFXU9iWurNOgyGBtMI2zcAjEop00BaCHEvsBs40IDXPiOczL05hG850KLZLJ9R39BDVe84gr6YVSFXFD5H3P1eD1Ouu2IZDM1AI2zVN4HrhBBeIUQYuAbYP89zlhUnllmlwkIzDqiq/Qj6rJqc/Xxw7n6P1fibdRgMS2Re5y6E+AJwPdAphBgA3ot9j3gp5SeklPuFEN8HHgMqwKeklHXLJs8GUXt9mSXPjmxitGg2y2eMVAexPTUVMs2yfw2nmrn78XsbcJs9g6HBzCvuUspbF7DNB4APNGSPGkDsvHDuM9fVOZfUiLs7llm1zl1XywTUUhDGuRuajFV55cXOg2qZpnPu/vM4lvGaOndD87Eqrzw9oNosefRy4PdaTbW8Qtg1oOo/HwZUXbGMWRXS0IysanFvlshiOdjWE2PHmgWsiX2W8Hst/B7rPIpltHMPmElMhqZkVeYW58PaMr997SZ++9pN53o3aggH1CJmWtAtoe7KtCqpirvPTGIyNCWr0lbpzL1Z8ujzhWjAS9jvqUYxAa+n8TfKaBbcM1TttfWlNBOZDM3DqnTuGzrC7FzbwqXr4ud6V84r3vPyS+iNB6vOfcmLhjUz8fUgPBBfh/+UhZRQqkh8nlXamBlWHKtS3FuCPu58x3XnejfOO168oxcAKSWWWMV5O0D7ZviLfgjE8B88DECxXMF3Dm+YYjC4MWeioeEIIezsfZXHYgF142ct6GZQ1dBMGHE3LAsBn7Wql1x2o0s/jbgbmonz4+oznHX8Hmv1O3cbv3bupmLG0EQYcTcsCwGftbozdxfGuRuakfPj6jOcdQJez+qulnFRFXfj3A1NxPlx9RnOOrGgl1jAN/+GqwA9oGpu2GFoJlZlKaTh3POBW3adP5l71bmXz/GeGAwORtwNy8IF3bFzvQtnDT2gmjeZu6GJMLGMwbBE/F41K7VYNrGMoXmYV9yFELcLIYaFELPeXUkIcb0QYkoIsdf+957G76bB0Lz4PSp+MtUyhmZiIbHMHcBHgM/Nsc19UsqXN2SPDIYVhs7czQ07DM3EvM5dSnkvMH4W9sVgWJHoxcKMczc0E43K3J8lhHhUCPE9IcSOBr2mwbAi8FrqMipXTOZuaB4aUS3zMLBRSpkSQrwU+AZw4WwbCiFuA24D2LBhQwPe2mA493hs527E3dBMLNm5SykTUsqU/f/vAj4hRGedbT8ppbxSSnllV1fXUt/aYGgK9N2mSkbcDU3EksVdCNEr7NvtCCGutl9zbKmvazCsFPRNyssVk7kbmod5YxkhxBeA64FOIcQA8F7AByCl/ARwC/B7QogSkAVeJ839xgznEca5G5qRecVdSnnrPI9/BFUqaTCclzjO3Yi7oXkwM1QNhiXiMc7d0IQYcTcYlohx7oZmxIi7wbBETJ27oRkx4m4wLBHbuJtYxtBUGHE3GJaIEAKvJUwppKGpMOJuMDQAjyWMczc0FUbcDYYG4LUEZbOeu6GJMOJuMDQA49wNzYYRd4OhAXg9lqmWMTQVRtwNhgZgnLuh2TDibjA0AFMtY2g2jLgbDA3AEsa5G5oLI+4GQwPweoTJ3A1NhRF3g6EBmMzd0GwYcTcYGoDXElSMuBuaCCPuBkMD8FiWce6GpsKIu8HQAFS1jBF3Q/Mwr7gLIW4XQgwLIZ6YZ7urhBBlIcQtjds9g2FlYDJ3Q7OxEOd+B3DTXBsIITzA+4EfNGCfDIYVh6lzNzQb84q7lPJeYHyezd4BfBUYbsROGQwrDY8lKJmFwwxNxJIzdyHEWuDVwCcWsO1tQog9Qog9IyMjS31rg6FpMHXuhmajEQOq/wL8uZSyPN+GUspPSimvlFJe2dXV1YC3NhiaA1MtY2g2vA14jSuBLwohADqBlwohSlLKbzTgtQ2GFYFHmHuoGpqLJYu7lHKz/r8Q4g7gTiPshvMN49wNzca84i6E+AJwPdAphBgA3gv4AKSU8+bsBsP5gKmWMTQb84q7lPLWhb6YlPJNS9obg2GF4jEDqoYmw8xQNRgagJmhamg2jLgbDA3AzFA1NBtG3A2GBmCcu6HZMOJuMDQAUy1jaDaMuBsMDcA49/ObXLFMplA617tRgxF3g6EBqLVlTCnk+cr77nySt312z7nejRqMuBsMDcA49/ObExNZTk/lzvVu1GDE3WBoAKZa5vwmky+RLzVXz82Iu8HQADzGuZ/XZAplCk0WyxlxNxgagNc49/OaTKFE0Yi7wbD68FjqUqoYgT8vyRTKFEwsYzCsPrweAWDc+3mKEXeDYZXisZS4m9x9cUykC0i5so+ZlJJ0oUSpIpuq52bE3WBoAF5LO/fmcm/NzEgyz9X/8CPuPzR2rndlSeRLFXT71EyDqkbcDYYGYJz74hlN5SmWJaemsud6V5ZEOu/MTDXibjCsMhznbsR9oeSK6rbLzVYfvlgyBef20c2Uu88r7kKI24UQw0KIJ+o8frMQ4jEhxF4hxB4hxHMav5sGQ3Ojq2WMc184uaISwnyxPM+Wzc2KFXfgDuCmOR6/G9gtpbwMeAvwqQbsl8GwojDOffHkS6vDuaddC4Y1U637vOIupbwXGJ/j8ZR0hrsjgDm7Decdls7cy+b0XyjauTeT2z0TsivYuc+LEOLVQoingO+g3LvBcF5hqmUWz6px7q4B1Wb6LA0Rdynl16WU24FXAe+rt50Q4jY7l98zMjLSiLc2GJoCUy2zePI6cy+tosx9JcUyi8GOcLYKITrrPP5JKeWVUsoru7q6GvnWBsM5RTv38gqfkHM2ya0S5+4W92ITfZYli7sQ4gIhhLD//wzAD6zsWQkGwyLRzr1kMvcFUy2FLDaPIJ4J7jswNZNz9863gRDiC8D1QKcQYgB4L+ADkFJ+AngN8EYhRBHIAq+VK30+scGwSPTaMiaWWTha1JtJEM+EZi2FnFfcpZS3zvP4+4H3N2yPDIYViK5zN6WQC6cay6zwOnd3KWQzibuZoWowNACvGVBdNNVJTE0kiGdCJn8eDKgaDOcrHlMKuWicUsiV7dwzhTJq1NE4d4Nh1WGc++JZNc69UCIWUAm3ce4GwyrDY5YfWDSrp1qmTGvYD6yyUkiDweCaxGRKIReMduzN5HbPhEyhRFvYBzTXZzHibjA0AOPcF4+z5O/KztzTece5m8zdYFhleM2Sv4umuvzACo9lssUyLSHbuRtxNxhWFx6z/MCiWU0Lh0X8Hvwei0ITxXJG3A2GBuBUy6xsoTqb5Jp44bCbP3o/H7vn0IK2zRTKhP1e/F7LOHeDYbVh1pZZPHqGaqFUodlWLDkynOLoSHre7aSUZAolIgGPEvdy8zRURtwNhgZg1pZZPHpAtSKbbyA6X6osKC7KlypUJIR0LGOcu8GwujDVMovHLZ7NlLtXKpJCubKguEgvGhbxe/l/7Z15kGRHfec/WfXq7qquvrtnuqdnRjOj0aB7RhdCWMAEIJZF2MIEGIzYtSG8tlljll1jEwtsrIxXOCBiN4wly16MOSTEAlq0BgwjkNCCzhlppJE099Vz9H3Vfef+8V6+elVd1ceouvrY/ER0dPerqle/yqr65i+/+ct8HkOQX0UjNy3uGk0D0NUySyeTL+L3mO22mjYPU7XqmUVU8airMOnMXaNZp+jMfWlIKcnkS0T8Zgnhasrcl3KFKGfm7jXcq+p1aHHXaBqArpZZGio7bg2sQnFfQolmIpsHcEyomo9J5QrsuWcfjx8ZW75AF0CLu0bTAJYjc/9vPznCf/7frzTsfKsJZXm0rsLFP0rUF2PLjMdzAHS2+PC6hb23zPBsholEjuNj8eULdAG0uGs0DWA59pZ5+tQkD+8/Z/u66wnlsUfszL3SAimtoL21lK2IxxNZALrDvorMfTppin4yu3JzCQuKuxDia0KIMSFEzRRCCPEhIcTL1s9TQohrGh+mRrO6cYvGZ+6xdJ5cocSTx8Ybds7VgsqOa9kyuUKJm/7q5zzy4vkViS2zhG0RJuJZhID2kLdiQnU6Zdo1zuurNpvFZO5fB945z+2ngd+QUl4N/FfggQbEpdGsKVwugUtAqYGLcWJpUyD2HR5t2DlXC6rGPeI390F3CunwbJrxeJaTYwsvIlLEMvmGjXDmy9xfGJrmtYsx+/+JRJa2oBfD7cJruMhXZ+65VZy5SymfBKbmuf0pKeW09e8zQH+DYtNo1hSGy9WwzF1Kyawl7r84MkZhFW0l2whUdlzLlrkwkwYgsQSx/sNvvcBnHznUkNjm29Ds8z98lS/99Ij9/3g8S1eLDwBPReZuintqBS21Rnvuvwf8pN6NQoiPCyH2CyH2j4+vv6Gm5v9v3C7RsDr3dL5IoSS5ZiDKTCrPi+dmGnLeakZjmWU570IoMa81oXpxxowpnlm8MA5NpexO4fXHVv8KUVPJnD2iAjNz7wyb2/16DZf9mKnUGsjcF4sQ4i2Y4v5n9e4jpXxASrlHSrmnq6urUU+t0awKDJdo2N4yKmt/8/ZOAI6N1q66ePrkJF/Zd+ySnuOFoWlu+uLPef5M3YH5sjE3c3eKuynS8Ux+7gPrEMvkl9QZzIfqeHLF0pzOOpbJV4woJhI5Oq3M3VdjQnW1e+4LIoS4GvgH4E4p5WQjzqnRrDXcbjFvnfsDT57klQuzizqXEvcdPWG8houhyVTN+/3o0EW++viJSxoxHBk2O4xHD16c937LUblS9txr2DLTS7NlpJTE0o0U98rJXUWxJIlnCiQcz1Nty9ieuzWh6qyWkVI2tQrodYu7EGIT8APgd6WUl5ZCaDTrAMMl6nruxZLkiz8+wvcOLK4CJJY2BaQt6GWwPciZydqTi6lckWJJMpnMLjneoSmzw/jpqyMVolMolrj3X44wMpthLJ7h+nv28bNXR5Z8/vmYr1rm4uzSxD2RLVCSS8v0543N4bU7Ox0l6iquZLZAOl+kM2yKe0W1TI3M/eHnz3Hrvb9omsAbC91BCPEQcDvQKYQ4D3we8ABIKe8HPgd0AH8rzHKwgpRyz3IFrNGsVubz3JPWl3zS+tIvhMrcIwGDwY4gZ+tk7ikrMxyLZekO+5cU77lp85xj8Swvnpth92AbAEdH49z3xEli6TwD7UFmUnkODE3z9jf0Lun882Fn7oG51TL2hOoiM/GYQ3SllFg6dMk4Bd3Z6cSszkM9z4RV465sGed+7mpC1Zm5Hx9LMDybIZEr2COW5WRBcZdSfnCB238f+P2GRaTRrFHmq5ZRlsFEfHEZtpq0aw14GOwI8asTEzWFS3UaY/EM0LqkeM9Npbh2IMqrF2f56asjtrhPJkxh+sELF+iystIzE4svS1wMai/3altGSln23BeZuc9aFkhJmiOZkG9BWZsXp6BnHBuaqQ63JM0J73HrvVRt5DXM979UkjXr3NV7OpvKN0Xc9QpVjaZBuFz1d4VUlsFi7RM7c/d72NwRJJMvMVajY0jnypn7fNSyAs5NpXjDhgjXb2qrmFRVMabzRYamUvgMF2cmao8cLhWVqUf8HoSgYvFPJl/CZ7iWkLmX7ZjFWjlfffwE77vvKWZS5ZHUQ88NcXw0XncrYmeVTCJbcGTuZrWMx+2yHzNTo1pGxTmbbox9tBBa3DWaBjFf5q6ESmXFC6GEIOw3GOwIAdS0ZpR41BJ+xZPHxtl9zz5OO7LveCbPdMq0XbrCPjv7dca4qy9CyOvmrt39nJlMNtQrVpm7z+PC5yghVJOp23tazHJQa4IyVyhxos4+LU6xXMyk6k8ODfPXPz3K/rPT/MG3DpArlCgUS/zFI4d4+PlzFdm60y6q6EQyBcatdupyVMuAWR5ZktAW9JArlOxJVjWPEmvQ3MBCaHHXaBqE6bnXrpZRojOVyi1qQdJsOk+Lz8Bwu9hsiXutSdW0NewfjWX40cvD/Pb9T825ZN3RkTjTqTxf/PFh+9i5KVNEN7UHaQ14KgRyPJHF63bxwEd28+DHbmZXX4RsocRIA2viVSmkz3DhdTvE3bJkdvSEgbJn/a1nzvKu//6rmpOmsQpxn184k9kC//F7L3Pdpihfuutqnjk1xXeeH2IqlUNK0+aqsGUKc20ZMDP3ccfWA2DaMlBeO7CxLQCUtwVWoh7TmbtGs7YwakyoKiFXX2wpy2Vy8xFLF+xKkg1RP4ZLcLaGuDsz932vjfD8mWlmqs4/kzYzzH2vjfLrExNAuVJmoC1INOhhJp23O4XJRI6OFi/9bUGuGYiypdPqXBrou2cLRXyGCyEEPo/b9tyV376z1xR31W4vnZ8hV6xtTVWL7nwMz6ZJZAt89I2bef8NA4S8bs5MpOzRSiJbrLhwSDZf4lMPH+TvfnnSzrzBzNwnElnara0HwKyWAexOsD8aBMq+u7ZlNJo1SnW1zFgsw41f/DnfP3C+auHLwr77bDpP2Np3xXC76G8LcKaGLVP23DMcGTFti9F4ZYY9k8oT8Rt0tvh46LkhAM5blTID7QFaAx6KJWl3FJOJLB2Wjwyw2RL303XKMQH+8NsH+O7+cwu+LkU2X8LvcQNm9q7sj/PTaQIeN/1tpjCqdlM1+bVsrZjDilnIllGP7wiZVkpX2Md4IlsW90y+ynMv8qsTEzx5fHyOtz8Rz9qVMlD23Eet+Q+VuavRh+octLhrNGuM6jr3L//sGFPJHMfG4hWisxjfPZbJ25k7wGBHaE7mLqW0q2UuzGQ4OZ4A5k6uzqTydIZ9bO0K2bcNTaUI+w1aAx77eZToTCZztvgB9EX81qRqbXE/NZ7gx4dG+OXRyi1Fpucp+8zki7ZH7fTcT08k2NwZsju2RLZAtlC0X9tkjY6xYqIzM38nqkpRVefVFfYxFsvYk8jJbLGqWqZEPFNgeCYz15ZJZO1KGSjbMmMqc7dtmQKlkrQtIy3uGs0aw5m5Hx6O8d0DZiY7lchViE51xUzSEjAnsXSluG+IBhiZrczIs4USUoLHLZhIZO2LM1fvFzOTztEW9NId9lklk2alzEBbECEErQFT6FSFx0S8MnN3uQSDHUFO16mY+YV1taHh2fLeLo+9Nsqev3yMc1O1H2NeP1Vl7uXL052aSLK1K0SLVc6YyBQ4MZawO81agh1L5+3dJVV2fWIswQ1/+Ri/rNouWXUOHZZP3h32M57I2mWNqjNRFotaqHRhJs1sOk/I67bvN5HI2pUyUBZ325axRh/JbJFkzlxoZcbbnC0JtLhrNA3CcLnsvWUefHaIgMfNQHuA6VSOeCZvX4pvvMo3/p1/eJZ7/vlwxbHZdN7edwWgN+JnIpGr6ATUFrcDlogoqn3p6WSeaMBDd9hv33Z+Om1nls7MXUrJRDJnV4AoNneE6q6SfczaktjZ+Tz60kWKJVn3MZl8yb44ts/jIlsoki0UOTeV4jJH5h7PFmxLBsy9XKqZTefZEDVfi7JxTozFkRJ+/PJwxX1V5t4WKmfu4/GsfTyZK5DNl+zFVepiHNlCibOTKft54pkCE/FchS3jnFD1uIWd1adyhQrrSGfuGs0aozpzv3JDK4PtISaTOeKZAj0RPx63qFilWipJDg/H7DrzF4amOTwcm5O597Waq0+dlouqwlCeuMctaPEZjMYyFfuYzKbztAY9dIV9pHJFEtkCF2fSticcDVrinjI3xcoVShWZO5jVK2cmksQzeZ47PcXnfviKuS1xKs/zZ6bxGi5G41mKJUm+WOLxo2Y2P+qIt1As2RPMY/FMefGPVS0zNJmiJOGy7hZafGZMiUyBIyMxfIaLiN+ouU4glsnTFvQS9Lpt+0vtLPnzI2MVJZyTiRzRoMf2x7vCPuKZgl2CaY6iSnbH6hwpHBuN0xX2YbgEY7EM6Xyx0paxzjk0maIt6KXFZ2b5yVyxoopHi7tGs8Yw3IJCqYSUkqOjcXb0ttAe8jKVzBHPFgj7DTpCvopVqhPJLLlCiVPjSfLFEn/68EH+6MEXSOaKFasYeyxxd5Yj2uJulUpe1tVCX6uf0ViGR1+6yO579pHJF5lJlW0ZMC2LZK7IxujczL16wlHxpu2dFEqSX5+Y5G+fOME3nj5LLFPgyePjFEuSd1/dR7FkLsl/9tSULbJOi+hT332JP37wRcAsedzQaj6/WS1Tsn31rZ0ttNiee54jI3F29ITpjvhrzlfMWh1h2G/Y9peyiCYSWQ45NmubSubs0kUo16gfHo5Zz2faMqrtnaOsVK5Ia8BDyGfYk9u1MveLsxlu295F0Gu+hlS2YFsxHrfQ4q7RrDVU5j4SyxDPFLi8J1wW94y55Lwz7K3I3FXGmCuWOHB2mrOTKU6Nm1ZGa6C8jF5l7k7rQ02mbuk0bZmdvWF6Iqb18usTE0yn8pyeSJLMFU1bJmIK0UvW3vBK3FXmPpPO25lxdea+e7CNsM/g/7x00S6nHJnNcGQkhtsleIe178zwbIZ9r43g97gIed325CLAwXMzHDw3Q65gljQqi8Oslily0nrdW7pChLxuhDDtj8PDcXb2hukIeWtXy6QLRAIGLT7DtmWGZzN0hLy4BHzzmbN88+kzZZ/c0XF1WW1yyposzuRLJLNFu8Or9vhbAx5afIZtN3XWyNwB3re7n5Al7slc0Z703RANNK3O/fVtwqDRaGzcwqyWOWqVJO7oCTOdMreinU7m6W8LEPC6KwTDeYGJ7z5fWUro9Nx7InPFPe2wZaJBDzdu6WD/2SmePZW0vX91SbhoyGtbCActcVfiGvC47YxyPG6KZ2eV5+5xu7h1Wyc/OlT2sC/OprkwnaY34rf9+5HZNE8cG+dN2zo5O5mybZl8scSFmTTFkuTsZBIpy52L2gf91HiSnojPnkxt8Rmcnkgykciysy9CKlfk8EiMalTm3uL32BOqw7MZdvSEKZRKfO/Aeb534DyG28VkMse2rhb7sSpzd5awTiVzDHYErYnqys4kYo0Q1P76XTUy9/62ADdtabcngVPZgh3XQFuQIzVew3KgM3eNpkGozF198S/vDdsTd0NTKVr8Bh0tldmnytyFgB8dGsbtEtxmXaDD6blH/AZBr5thZ+ZuZanRgJenPvNWPnDDgJW5ZzhqxfCaZTeoCVUoZ+5K3M2KGXOVqsrcq8Ud4C07zQvsKPEdnslwcSbDxmiAPstiefVijLOTKW7Y3E5PxG/X3A/PZGwBfebUZMXzB71uxuNZXj4/w9bOsvCGfQb7z5hX8LyiL0xHi3fOxmu5Qol03rRRIn5H5j6Tpi/q5967rua+D11PwOPm+GiCqWSuYlTS7ci8VcXNVDKHz3DhM9y2LaMmwyN+g5DPsCtf1FWYwOwkwczaXS6B11p968zc+9sC9sT1cqPFXaNpEKbnLjkyEqcn4iMa9Nold+l8kbDfoKvFx0Qia3+5L8ykCfsNLutqIVsocXlPmLtv2QxgCyaYAtwb8Vd42GlrJWXQ5yboNXC5BN1hH/mitDfisjP3oIdowIPHLTg1kcRruOzYwOxIZlNlz93pSyt+Y0c3QljiJcws/YI1MdsW9OA1XPzLK+a+79cMROmO+OwJ4LNT5aqZp06a4q4mdD988yDZQonjYwku6w7Z92vxG/Ycw87eCB0hH7FMoeICGiojjlh2STxToFiSjMazbGgNsLWrhTuu6uOy7hDHRuNMp3IVr7ujxYel2/bEdK5Ywme48Rku2x/f2hWy20l1bkJAe7B8rsGOIF/+7Wv42G1b7WMhn7uiWmagPUi+KO33bjnR4q7RNAi3y0XJytzV3ihtji9/2G9WrGQLJdt3vzCdZmM0wOXWcvvrNkXZu6uHJz59O7s2RCrO39vqr6glVysfg1btNZTtGzDFR9kYbUEvLpewM/INrX5crvL2wXbmnsgS8Ru2xVD9/N/52M186u076A77OTedZiRmZu5CCPpa/RwfS+AScNXGVnsUIaWs2PRMZe5qHuHq/ij3vPdKoLynDJRHCD0RH+0hr50lTznmLJxbI6sJ1bG4OUroi5bbYltXCy8OTSOlKegKt0vQbnnwm9rLJaU+j8uuwwfYbsUVCXjsyd6OUHnrAbO9BXft7q/YcjjoNUhmzcw96HXbnWYzJlW1uGs0DcJwCaZSOY6PJrjcEgOnBdDiWVDziQAADVxJREFUM7je2jP9udNm6eOFGbPeXN3/2oEoUM4inZiZu7N6w8wGVVUGmEII5uTelRta7X1mlMWjbAhliSiiQS8z6RxjVUvqq7lpawcRv4feVj8Hz81QLEk7A1divb07TMhn0GONIqZTeYamUngNF2G/wXQqT2eLt0I8379ngO//u1t4/56BcntZFSs7e81OTlXwOOcsnBc1afF5iGfydhnkBsfIZ1t3i729QvVksZqLGOxwiLvhslfQ+j0uW/gjfg9hS7znaydFOXM3J9SrVwMvJwuKuxDia0KIMSHEK3Vu3ymEeFoIkRVCfLrxIWo0a4ON0QAzqTyFkuTWbaZv7szcI36Dqze20uIzeOqkWXGiMvebtrTjM1zcvLWj7vl7rTJHVbetSiGdmbvy1bd1t7DJIVblRTvm7dXirjL314ZjFdlzPTZE/fYWwupcykZSHZQaRYzGMgxNphhoC9ibkFU/P8DuwfYKwVciurPPjEetBnVWGym7Q2XuyVzRnqSuyNy7y15+teXUbYt7uUP1GW579BL2e+x4I1YppBnPwuIe9BqW525W9KgSy2asUl1MtczXgb8BvlHn9ing3wPvbVBMGs2a5D+8fQcfu20rIZ/bHq63BcuTomG/B8Pt4sYt7Tx1cpLZdJ54tsDGtgA3be3g0BfeUdMOUfS1+imUJBNJ85J6qZy5TN7jsAZUueMVfRHaQ+ZzGy5hL5vvqpO5twY8jMxmyBclv3PjpgVfq3M+QFW99FqZ+zWWuHc7xP3sVIpN7UFCPoOXz89WZNX1ULbMrj4rc7fE1Lm/jJoc7msN2KtaT1iTyX1VmbuiWpTtzL29MnNXHU3Yb3DD5ja2dobY0lneGsG5gKkeIZ+bVLZAoVhafZm7lPJJTAGvd/uYlPJ5oDnFmxrNKkUIQWvQU+HDGm6X/YVW4nPL1g5OjSd5YcisBNlobQ07n7DD3HLIVK5A0OeuuI/PcPOne3fw4Zs32QIeDXrty/OpLHVjtPJ6q60Bj703zXWb2hZ8rcqCMc9lPs8GW9xbrXjN5zIz9ySDHSHb+lBWznwob9u2ZazMXdkyhWKJh54b4rbtnWyIlsX96GickNdtV78AbGoP4bbmGKoz996IH5egYqSjLiICZqe8szfCLz59O+0hr/08nS1zJ52rCXoNUrmiacvU2KRtOdF17hrNMtMR8toX3wC45TLTern/iZPA4oQOypnoxZkMV/ebE6pBj3vO/f5k73agvPQ/6hg9zJe5gznBeNXGha/FqmLpCHkJWKOCO6/bSMBr2Jm2eq7Dw3GSuSKb2oO2YNeyZarZ3BGks8VrV6qEfQZet4t//PUZ9r02yo1b2hmezfBf3vMGAHvLgoPnZuizJnkVXsNlbX6WrLDKAO5+42b2bG6rEH2/x21n7s5OAliSLRPymp57ScK2LmN1Ze6NRAjxcSHEfiHE/vHx8YUfoNGsA5TfHbb81l19EbrDPp49PUVHqCxeC3FZdwiv4eLZ02a1STpfIDjPxaDtFaiOevlrB6L0RHx2NqxQHcDO3rAt1vOh/GxnxxTxe3jf7n5bVH2Gm7agh32vmRuLbWoPcpn1WgcW0aF96KZB/u9/eqttOwkhePc1fXRH/FycyfDVx0+yodXPW3d2A87MPscdV/bOOd+2rhbag147g1d0hX3cfnk3PsNczGXGXs7cqy9mvRRbJugziGUK9kZwqnNbd5m7lPIB4AGAPXv2LH8Vv0azCmi3xd38urlcgn/+xJvIFkpsiAbmiE09gl6DN23r5LHDo3zu3bvMzH0eIVYCHHVkqldubOXZv9g7574qo7xuU3RRsSjPfOMCGfiGaIBXL8a4dVsHb9zWQcDj5u9+d7ctyPPhcok5Hc1X3n8tYC7guu+Jk1wzELVtsBs3t/OPH72Bq/pba2bVn3jrdvsKVPUI+QxmUnmzzt3huTtp8S8+c7+iN8yDz5oXSIn4Pbhdgjuu7K0ou1wutC2j0SwzHVXiDuXJxqWy94oefnFkjONjCdK5+cW9I+TFa7gqbJl6qNHFdQML++1gZq0+w1XhU9fi3ruuJpktcOOWdjujV/vQvB5CPoNPv+PyimMul+At83QaV/W3clX//JZTyKvE3YXf9twrZfKWrR18cu92btravmCcH755kM4WH/f/8iQ3bjHvf9+Hdy/4uEawoLgLIR4Cbgc6hRDngc8DHgAp5f1CiF5gPxABSkKITwK7pJTN2UBBo1nldLb4cLuEPZx/Pbztim54xLweajJXqFi0VI0Qgj+/YydXLsJDv7Y/yl/91lW8+5q+RcXhdgke/NjNdmljPRbz3KsJ9R75PC58nvKEqhO/x80n9+5Y1PmEENxxVR93XLW4dm0kC37apJQfXOD2EaC/YRFpNOuMj9wyyPWD0YoqmkulJ+Ln6v5WHjs8umDmDvBvbt2yqPO6XIIPLqIE0snuwcVl+WuJkK98dSifUduWWSvoFaoazTLTHfHz1p09DTvf3it6OHhuhpFYZkFx1ywNVQlj1rnXztzXClrcNZo1xt4repDSXKHq3HpA8/pRWbrO3DUaTdO5oi9sV6nozL2xqAts+CsWMWlx12g0TUAIwd4rzKqQUAMmaTVlyraMcxGTtmU0Gk2T2LvL9PADNVaoai4dZ7WMqrF3XjRlLaG7fY1mDXLz1g4++sbN3H5510qHsq5QmbvX7eKOK3vxuAUDTVhwtBxocddo1iAet4svWHuqaBrHv7qqj0KxRDToQQjBb163dqu8tbhrNBqNxaaOIJ942/aVDqMhaM9do9Fo1iFa3DUajWYdosVdo9Fo1iFa3DUajWYdosVdo9Fo1iFa3DUajWYdosVdo9Fo1iFa3DUajWYdIqRcmUuZCiHGgbOX+PBOYKKB4TSS1RqbjmtprNa4YPXGpuNaGpca16CUcsF9J1ZM3F8PQoj9Uso9Kx1HLVZrbDqupbFa44LVG5uOa2ksd1zaltFoNJp1iBZ3jUajWYesVXF/YKUDmIfVGpuOa2ms1rhg9cam41oayxrXmvTcNRqNRjM/azVz12g0Gs08rDlxF0K8UwhxVAhxQgjxmRWMY0AI8bgQ4rAQ4lUhxJ9Yx78ghLgghDho/bxrBWI7I4Q4ZD3/futYuxBinxDiuPW7bQXiutzRLgeFEDEhxCdXos2EEF8TQowJIV5xHKvZRsLkf1ifuZeFENc3Oa6/FkIcsZ77ESFE1Dq+WQiRdrTb/U2Oq+77JoT4c6u9jgoh3rFccc0T28OOuM4IIQ5ax5vZZvU0ojmfMynlmvkB3MBJYCvgBV4Cdq1QLH3A9dbfYeAYsAv4AvDpFW6nM0Bn1bEvAZ+x/v4McO8qeC9HgMGVaDPgzcD1wCsLtRHwLuAngABuBp5tclxvBwzr73sdcW123m8F2qvm+2Z9D14CfMAW6zvrbmZsVbd/GfjcCrRZPY1oyudsrWXuNwInpJSnpJQ54DvAnSsRiJRyWEr5gvV3HDgMbFyJWBbJncA/WX//E/DeFYwF4G3ASSnlpS5ke11IKZ8EpqoO12ujO4FvSJNngKgQoq9ZcUkpfyalLFj/PgM0/dpvddqrHncC35FSZqWUp4ETmN/dpscmhBDA+4GHluv56zGPRjTlc7bWxH0jcM7x/3lWgaAKITYD1wHPWof+2BpWfW0l7A9AAj8TQhwQQnzcOtYjpRwG80MHdK9AXE4+QOUXbqXbDOq30Wr63P1bzOxOsUUI8aIQ4pdCiNtWIJ5a79tqaq/bgFEp5XHHsaa3WZVGNOVzttbEXdQ4tqLlPkKIFuD7wCellDHgPuAy4FpgGHNI2GxulVJeD9wB/JEQ4s0rEENdhBBe4D3A/7IOrYY2m49V8bkTQnwWKADftg4NA5uklNcBnwIeFEJEmhhSvfdtVbSXxQepTCKa3mY1NKLuXWscu+R2W2vifh4YcPzfD1xcoVgQQngw37RvSyl/ACClHJVSFqWUJeDvWcbhaD2klBet32PAI1YMo2qIZ/0ea3ZcDu4AXpBSjsLqaDOLem204p87IcTdwLuBD0nLoLVsj0nr7wOY3vaOZsU0z/u24u0FIIQwgN8CHlbHmt1mtTSCJn3O1pq4Pw9sF0JssbK/DwCPrkQglpf3P4HDUsqvOI47PbLfBF6pfuwyxxUSQoTV35iTca9gttPd1t3uBn7YzLiqqMimVrrNHNRro0eBj1jVDDcDs2pY3QyEEO8E/gx4j5Qy5TjeJYRwW39vBbYDp5oYV7337VHgA0IInxBiixXXc82Ky8Fe4IiU8rw60Mw2q6cRNOtz1oxZ40b+YM4oH8PscT+7gnG8CXPI9DJw0Pp5F/BN4JB1/FGgr8lxbcWsVHgJeFW1EdAB/Bw4bv1uX6F2CwKTQKvjWNPbDLNzGQbymBnT79VrI8zh8letz9whYE+T4zqB6cWqz9n91n3vst7jl4AXgH/d5Ljqvm/AZ632Ogrc0ez30jr+deAPqu7bzDarpxFN+ZzpFaoajUazDllrtoxGo9FoFoEWd41Go1mHaHHXaDSadYgWd41Go1mHaHHXaDSadYgWd41Go1mHaHHXaDSadYgWd41Go1mH/D+vZjhDK1CxgAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"dists = md.compute_distances(traj, atom_pairs=[[0,1], [2,3]]) / sigma_nm\n",
"# TODO: use \\tau as time\n",
"plt.plot(dists)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Double Well Chain with 3 Atoms\n",
"\n",
"This is similar to the system studied by Rogal and Bolhuis. However, this test system allows one to create a polymer chain of arbitrary length (up to all the atoms in the system, although longer chains may have strange knots in them in the initial conditions)."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"dw_chain = testsystems.DoubleWellChain_WCAFluid(nchained=3, nparticles=216, **potential_kwargs)\n",
"integrator = integrators.BAOABIntegrator(temperature, collision, timestep)\n",
"sim_dw_chain = mm.app.Simulation(dw_chain.topology, dw_chain.system, integrator)\n",
"sim_dw_chain.context.setPositions(dw_chain.positions)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"sim_dw_chain.minimizeEnergy()"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"n_steps = 200000\n",
"interval = 1000\n",
"reporters = [\n",
" mm.app.DCDReporter(\"dw_2_chain.dcd\", interval),\n",
" mm.app.StateDataReporter(sys.stdout, 10*interval, step=True, potentialEnergy=True, \n",
" temperature=True, elapsedTime=True, remainingTime=True,\n",
" totalSteps=n_steps)\n",
"]\n",
"sim_dw_chain.reporters.extend(reporters)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"#\"Step\",\"Potential Energy (kJ/mole)\",\"Temperature (K)\",\"Elapsed Time (s)\",\"Time Remaining\"\n",
"10000,2.5537092685699463,100.29120066576628,0.00019502639770507812,--\n",
"20000,2.1571922302246094,91.33579188418442,9.913586139678955,2:58\n",
"30000,2.2482588291168213,93.78675791043933,23.54122304916382,3:20\n",
"40000,2.972226619720459,92.15136639187402,33.33867406845093,2:57\n",
"50000,3.5334603786468506,95.35766224773101,42.58765697479248,2:39\n",
"60000,5.468897342681885,106.27403643390603,50.56366300582886,2:21\n",
"70000,1.6885178089141846,98.51765307882634,57.91592788696289,2:05\n",
"80000,2.315932512283325,89.27350485476714,65.74040508270264,1:52\n",
"90000,0.6174218058586121,94.8512353787532,78.91926217079163,1:48\n",
"100000,4.524473190307617,91.12570584523502,90.78647994995117,1:40\n",
"110000,3.710754156112671,102.85923485737591,100.28000497817993,1:30\n",
"120000,2.2611424922943115,101.24817117750801,109.54292488098145,1:19\n",
"130000,2.3109912872314453,104.83908478215933,120.58458399772644,1:10\n",
"140000,2.853100299835205,106.97028275117917,128.04100584983826,0:59\n",
"150000,1.930417776107788,98.67910609727373,135.7251842021942,0:48\n",
"160000,1.8435741662979126,110.50581158125924,143.44128799438477,0:38\n",
"170000,5.065273761749268,104.22540018243073,151.14959692955017,0:28\n",
"180000,1.3028106689453125,103.35716539153508,158.36259698867798,0:18\n",
"190000,1.7675387859344482,99.71697304493554,166.1741020679474,0:09\n",
"200000,6.894049644470215,105.01972818885046,174.20511603355408,0:00\n"
]
}
],
"source": [
"sim_dw_chain.step(n_steps)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"traj = md.load(\"./dw_2_chain.dcd\", top=md.Topology.from_openmm(dw_chain.topology))"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "1e6f051d4c294cbca352fdf743efa8f3",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"NGLWidget(count=200)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# movie of trajectory in NGLView, if available\n",
"# note that the atom selection different from the dimers!\n",
"view = nv.show_mdtraj(traj.center_coordinates())\n",
"r0_nm = r0.value_in_unit(u.nanometer)\n",
"view.add_ball_and_stick(\"Ar\", radius=r0_nm)\n",
"view.add_ball_and_stick(\"0-2\", color=\"blue\", radius=r0_nm)\n",
"view"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x152214cbe0>,\n",
" <matplotlib.lines.Line2D at 0x152214cd30>]"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsvXeYJFd97v851Tn35Lg5Sbta7UpaAQIFsgCTBdgYePiZbBywjW1sX/Nzgnt9bYwDYJLBJJORMQIMQgKhHHYVdiVtThN2Qk9P59xV5/5xTlV3T17trFaS+32efWZ2urq6usJ73vN+wxFSStpoo4022nhmwbjQB9BGG2200cbqo03ubbTRRhvPQLTJvY022mjjGYg2ubfRRhttPAPRJvc22mijjWcg2uTeRhtttPEMRJvc22ijjTaegWiTextttNHGMxBtcm+jjTbaeAbCfaE+uLu7W65fv/5CfXwbbbTRxtMS+/btm5FS9iy33QUj9/Xr17N3794L9fFttNFGG09LCCFOr2S7ti3TRhtttPEMRJvc22ijjTaegWiTextttNHGMxBtcm+jjTbaeAaiTe5ttNFGG89AtMm9jTbaaOMZiDa5t9FGG208A7EsuQshviiEmBZCPLrI6zEhxE1CiEeEEI8JIX5j9Q/z3DE6W+S2w9MX+jDaaKONNp4UrES5fwl42RKv/xbwuJRyF/B84B+EEN5zP7TVxRfuPMnvfuOhC30YTwylFPzHG2Hm6IU+kqcEfnFomsfOZC70YbTRxlMay5K7lPJ2YHapTYCIEEIAYb1tfXUOb/WQKdXIV+qc1wXBH/s+lM8D6ez/Nhy9GU7ctvr7fhriz/7zAJ/55YkLfRhttPGUxmp47p8ELgbOAAeAD0gprYU2FEK8RwixVwixN5FIrMJHrxz5Sh1LQqW+4KGdO7Jn4Dtvh0dvXP19P/S1xmf8T8PBm+Dj26GSB8C0JNO5CoXKU04/rDrKNZNyzbzQh7EiPDqeIZGrLPialJJDk9kn+YjaWA1yvx54GBgEdgOfFEJEF9pQSvk5KeUeKeWenp5l+96sKmwyKFXP08NS0INVZfmb2LQkhydzK9vv5AGY3A+AzI4/0aN7+iJxCLLjMK76ECXzFUxL/o8g9z/67n5+75sPX9BjqJkW7/7KXh4dX3xGKqXkbV+4j3++9ciCrz84kuJl/3RH20p7krEa5P4bwI1S4RhwErhoFfa7OqgW4OPb2Zx7AIDi+VJCRe1cVYvLbvq9fWO8/J9vZya/sNJpwcPfQBoejlmDJCdW1C/omYVaSf0cuReAqaw6Z4XqU4TcJ/bDqTvPy65HZouMzC5/P51PTGXL/OzxKe46NrPoNplSjVSxxkS6vODr4/rvM/nqeTnGNhbGapD7CPAiACFEH7ANeOoYorlJyI4zXFbByNL5IoWSTe75ZTd9ZCyNJSG53M0uJRz8AYm+qzkk1yJyF9aWuevYDH/4nUee3A+tacJwyF39v1hZhUE6PQqPfu/c9nHrX8EPf//cj2UBFCv1ZQexHzxyht/++oPn5fNB2ZkA6VJt0W1GZ9UAPL2ILZMuqvv8vM2a21gQK0mF/AZwD7BNCDEmhHinEOJ9Qoj36U3+BniuEOIAcCvwISnl4sP8k42Ksj/8dWWXlKrnyXO3lXttYaV1+5EEN3z6bso107FkcuWmB+b+z8OB77a+aepRyIzyePQaJmUH4cq0IvwLhO8/NM53941RbCKc9351L996YOT8fah9PsceALPOVE6R+6oo95/+KXz3nWCew75mT0Lh/Nzuxaq5rP10+5EEP9w/QaV+fogzX9bkXlxciNizi8U891RB3efnM35QMy0s69yfjQdHUrznK3upm+eJJ55ELNvPXUr55mVePwO8dNWOaLWhyT1oKr+veN6Ue0r9XMSWue9kkn2nUzwymm6Qu/3g1qtwy19C5wbY+YbGmw7/BBDcZezBksfxyTKU0xDoOD/fYRkcmVLHncxXCXa6qZsWP3t8ipDPza9eufb8fKhty1TzMPUoU9kgsArKPTcJh34MSLXvQPzs92FZkBkFswaWCYbr3I5pDvKV+rKkndTW3mSmzLqu0Kp+PjTu0XRxCeWeUvf8TL6CZUkMQ7S8ni5p5X4eyf03/v0BBuN+/u4Nu85pP/edmOXmx6eYLVTpjfpX6eguDJ75Faqa3MOW+nnePfdaYcGXZ3LqBv/+w+POA5PTqoiRuxXBTB9s2BAAh38Mw3t4NONjUnYBUJkdOz/HvwwsS3JkSllOdqxgJl/FkpBtnrLnJuEHv9P6Pc4FtSL4Yur3kXuZzjaU+1mntTYr9Ie+ClLfC5UVBrfnIj8JZhWQUEo/sX0sgWK1TrlmLakikwV1X42nSqv++dBQ7qkllPuoVu51Sy64nT0wnC9bZjpb5s5jMxybXt4SXQ62bZstLz6YPV3wP4bcYyjSbb7B7jiaYDq3SiRUWjqgmtCE+P2HGr65Y8sc/Zn6adVh+jH94iSceRC2voxTyQIZj8oumhg9vjrHe5YYS5Uc5TWrCWVSE2221ESaJ26DB7/S+B7ninpZzWiCXZA46HjuloRy7Symztkz8Hcb4dCPlLW17yvg0rV2K4iTcO9n4LPXwieugIweYFNNAe7S/FKQumnxwo/dxo8PTKz8ODUqdZOaqQavwhKkaMdtxtPnidxXpNwbn51YIEnAJvzFlPtUtrxkNs5y+IWuPLfvy3OBfYzZ8lMkYH8OeOaTe1WRe1xo5a4flJpp8Y4vPcAnbj22Op/jZMssrNxtP7L5BneU+9GbmfRqW+OMTn2bVN0eKkNXMZEps2nTVgBSk0tnzNx9bIZX/PMd/NsdJ5b1ONPFKqkVPhCHpxrq1iYUm2gzzcrdVrC2TXWuqJXAE1RWVCXnZMvAWfruD34FKhmOP76PY6dGIDPCTPezAJDlFeRg7/t3SJ2C5DEYvU/9Ld10LYrJeW/JlGqcmCk8oRTAZttpMd9dSkmyoM7HmUUyVc4VtnLPLBlQLdIT8QEwnZ1P7kspdykl7/+PB3n/fzzxoPAtB88DuS/xfZ8ueOaTe8Umd1u5q5s1katQMyX7Tp87CRUqdWp5/XAvYsskchWifhXiGIz5cRlCPTizJ2HmCN80X0RORGBCk7smjjHZDcDu7duwpKCYWDp4+eBIiscnsnzkRwf51C+WHrg++O1H+INvryyP+kgTuc9oQlmQ3O0K3dWyKWpF8ATAF9HkXsbvUbftinz34iyYNcx9Xwbg5w8e5I++8nMAjpu9AMymlyrA1ihnYPOL1e9JnQyWaiX3YrXueODQUL25J6ACmweuYrXOv91xgt/82r6WwbhYNZ3Zy3j6/KRM2hbiYraMZUnGUyWuWKviQAsFVdNLKPd7jqtY1FKDx1Io10zuPDqDyxBky3VqpkWqUG25DmcDW/ydi3I/Np3nl0ee3CLNhfC0Jfd/ve0YX7nn1PIb2uSOmnoXqybUqyQmRwE4NJmlUFEP5RON5v/vHx9kekrbLQvYMpalFNZLd/QDcNFAlLDPrWyZ8X0A/Ky8jaOuzQ3lnhkFw8PRUhiATf0dpIw4VmbpQqZC1cTjEmzuDS/rQZ6cKTCRWZniOzKVYygeIOh1MauV+6R+b4s/6ZD7uQ+atx6colYuOORulTIkC1U2dKtzsqxyz07Ax7bCJ6/EpdNIN4cruMqKzEelIvfpxMLZLlLKhq9fzkJkACKDMKutsfRpEDqIWkzy0R8d5M2fv9d5v03qT0QFFpoGrnzF5JdHEvz3o5O85lN3OR53cyrtqin3B78Cv/y7puNQ36FcsxZ8PqZyZaqmxeXrVEB6IVvGTqNciNw/8fNji762Euw9laJUM3nBNmVbpopVPvidR3jlJ+58QgRvf8fcOXjun/j5Ud7zlb3O4H6h8LQkd8uSfPq243zkRweX9xo1uQdEFR9VRe73f5aLb3wxBhaWVJks7/j4t/nHH+17QsdzcCJL2NJT+wVSITOlGjVTsn0gyusvH+I1uweJ+N3q4dckOGnGeExuUEHVegXSIxAb5tSsemjXd4fI+3oRuTN8b9/YosHEUtUk6HXTHfYuWyQ1nas4xPPRHz3Oe7+6d9FtD0/m2NYfoTPkdYJ4tuderJrU7KDfKin3umnxnq/uI5vLaXKPYmr7ZGO3ygpZtko1fRqsGmTPUPAPcNBaw4ZAmU5t0R2vqSD17Ox8SwXgzZ+/l7/970MqG6ZWAH8MujZB0ib3Eejbrn4vJnlwJM2pmaJzbWxyP1flXqjUyZRqbOoJkcxX+OsfPg40ZlBhn5szK/XcE0fg317SsBHn4sB34dEbneuZbzr2hXz3kaS63y/qjxL0uubZMqYlHVVenmPLjKWK3HMiSWfIS7VuYT6BVEY7Znb5OjVzSBVqnNKi5QPffPis92lbRy1xpLPEeKpEpW5x68GpJ7yP1cDTktyPJ/LkynWqdYt/+tnCJc8OmjIh4uSVQkidxlvN0Il67f/+cD9fNj/Eyw/8vkpvWwSFH/4ZxR/8sfrPqTvh8R8AMDKTIyY0qS/gudtqpifi4+Nv2s1rdg8R8XvU1E+Te4YQD5nrFRlNPaaII76G08kCXSEvUb+H7sENrHWn+OB3HuGH+xcO0hUqdYJeF+8tfo6L079c/LtU6uQrdYd4Hh3Psu/0woQ8ni5xIlFgS1+YrrDPGTRsWwaarJlyuvXnUsiMwcTCRVEz+ap6MGslh9wtTe4bbHJfJvuimlVe7NRrv8XndnyFGTroMvJ0CDWjOVhW5J5OLzzLODqVV7ad7ck75K7trtRp6LkY3H7MQpLj03mqpuVM6W3ltmzmRSHp9M4BIHmcbd97Cf0knf2kizUuHY7z/hdsdipG7RnUjsEo4+lSY8BPHFk8A+jkL2Hsfjh998KvF5PUzTo7//Kn3HM82aI+7ZTGZtjB1DWdQXojvnnKPVuqOaUZxTnXyx4sNvWo6/lE1Lt9fGs7VYpsslBhIlNmKB7gzmMz3H707OyRhi3zxJW7PRu+6ZGzD6SvJp6W5P7QiCKOF17Uy/ceHHOmqQuimdxFXuW5a3W5xpNmY0+IztmHiIsCu81Hyd3xr4vuKvHIT5l6+KfqIbrj43DLX5Ip1TCLihzq7qAi9zmq2vYh7aAT9SqXG0fIV2pQnKXuDlHHzQMVO6j6kKqejK9lIlNmMB4AINS7iTVimqjPYO+phZVXsWYS9Lp4bvYnXFFa5AGmUU2Yq9QddeXLj2F9/sXq853tyrzl8/fi8xi84fJhukLeRrZMpozQKc2O9XA2tszPPwpffd2ChVm2IvNYZeouP/giCE2AGzUZFJdR7tOTyor5yZibw2kXFW8cfy3tDOr7C8pKKGTnH6uUkmy5xqlksTFQ+WPQuUllxhRmIDsGHesg2EU+NU1Vq91ErgLHf8H2O36L97u+v7BKPngTJA6r37/6Grj5zxuvHb2ZUOYouww1QyhU6qSLVWIBD++8egPDHQH+7qeHnWDqpcMxKnXLmVHxhRfD3Z9Uv2fG1Gcd+pH6f1rHbHS/onkozFCrqRTMQ5PZRi0GjWKkZtx3IknY52YoHqAn4mM6W26ZVTZXts4lb5tIu0LquXgiqZI2ua/pUOR+OlmkVDN59e5BAE4kFo6BLYZztWVMSzKZLeNxCW4/krigKZVPS3J/cCRFLODhfddtwpIgb/kLuPnDC2/c1MiryyioClX9sG4NFrhibQcvNB6iLjzcae4gcPvfNJTaHATrKcJmiodG01CYhvw0I8miowTTnn6VO222Khxb6XaHNbnv/xYfTf4B3uIklFJUPIpkRmUPlj+usjHykxBfR65cJxrQtWY92xD1Mld3Fzm4SOOxYqVOyOvCa5WIWNlF4wjTTao7X66TL5b5J++nMMYfgPs+67z2jz87ykSmzJd+41ls6YvQFfI2ZctUWKcV0zzlvgJyLyROqiyT2fndKuzpvZ8KyYobfBHctTwg2doXUe/XZHAmXeL1/3pXy0wCoJpTyv2BKcGpZAECXbgrKTpEjqrhZ9b0U5EeyoWMIj2bEFEec82UzOQrlHKanG3lDkoBSwvi6yDQSSndWAgmkavAI99kYOJW/tjzbV5Q+O/WL1crwXffAXf+oyp+mj7UmA2AE4cZFioWkCvXyZbrxAIe/B4XN1w+zP6xtFP2v3M47pwHTC1esjpd84svh2+9Fb7561Smj6pYDqieOHMhJRRnMHU9QDJfJV+u0RVSKaOZOcq9Ujf5yWOTvHRHH163QW/Ez2S2zA2fvpuP36wGruZA7HxyV5/TFfbqc/4EyL1cx20IBuKq4OjghHp2L9ZxrSWF3wIoLmDLfPWeUy0L/ZxOFnjhP9zWsm/TkliWZDpXxrQkr941RNW0uO3whQusPi3J/aGRNLvXxIkFPAB0jN7SUCZzUclT8KqMk7WBCqVaQ7lv9Od55a5BXuHbD+uv5ofGC3GbZcg3rdhUmIFyBmlZxKwsHeT49n2n1d+rOUamk06wdkroTpdzrJl5yl2TWaA0DaUUBVdEbymo9lwKR36i/htbQ65cI+JT35Me1Y/tOdEEhyaySiEdvKml/L1QNYl6TASSTpFbtCS8uQ9Itlzj9eUbudI4Qim6UdlNlTxSSu4/PMprNsIV2tPsCvtIFirkta2zRROtk11QWsBzz02pvH1QKlb/np3WRDM23+ufypVxYeIVJhNFwB/FwOSSHo8zk7E99wPjGR4cSc/LUDDzMxSlj/vHy5xOFnFFuhGVLMOuNBnduLRkBDGqOSoPfgNu/l/OcTdnb0xP6/vBF1XKHZj6ycfU3wYuhWAnZr5xDWbyFUifZiK2m4r0OK0vHJx5SAmA2ZMq/96qNc4POOQ+pMn9TEaReDyo7oMdg1GkhLuPzxDyuhxbYzxVgrratpxNIM06MjPKA9Y2AB64/56llXs5DVYdU1uT9nUe1gN4ao7nfvuRGXLlOq/apVRyT8TH6WSRB0fS7BtRg7udKdMR9Mwj75Kj3BW5z7VtVoJCpU7I56YjqPbx+Bl1rgdiftZ2BhmZVTGQb+8dXVEn0Uaee00fU52/+dFBvnz3KWebvadSnEgUuKXJU/+zGw/w1i/c5wS2r9qkLb8lir/ON5525J4t1zgynePytR2a3CWBwri6aa0Fbo5KjoxP3XwD3pK6gfQDvMab47ruHEPmGO5tL2OwT5OzrfYnH4VPXgk/+F1SmTQ+UcMtLG4/cASpW/wmJkbp0AG6U3XdFqBagDv/yQm8JXIVvG7DSYVEt+71V5NQSqkUSI1c546GtRFfS65cJ2K/r0fluu/wTJAt15keOaxU2d3/ol6//e8ZLB6m06Nu4g5yiwZVm8k9nS/yFn7Iz83dPLj7r1Xw8OBNjMwWeX3hG/zF1O8623aFvNRMyYmRMV5gPMTWPpW5klnKlvmv34Lvv1/9/t8fgm+8GcuSROuaEMcemH982Qp+1IMxlpNYXnWOrhryEPSqDBU76GhbQvvH5vj8hSSzRJjOVSjVTIIxlR2z2TVJwlSEKL1hQqJMdkb7o3V1Xpqn07OzetDwx6BjPRaCvvzjsOE6GNgFwS5c5Vl69eCdyKmA+KxngCwB/GYe05Icm84pQtNN0EidbJBt8+CnB/9BkSRKgV959A/4pvdvuHJMpXNuH1QD08OjabrCPob0YDeeLmFWFbkfPXmaz/x0HwJJsu+5ABx6dB8yPaoyfLLjyuufc74ALNNOF66SL9cZ7lD7nxtQvemRM8SDHq7erMSTI15o+M72ewZigXm2iz3z6tIz2rnKfjJT5qePTbb8bWZOVluuUifsc+NxqefLVu790Qa5Pzya5o+/u5+bH2/dF6jZQrON1LBl1Dm453iSat1irKlQy+6lc//JWec9N+0/wwOnZjmdVMLOviZPJEi8Wnjakfv+0QxSwuXr4kQDbrrJ4rYqYNV41yd/wH0nktRMix/tn1DTvkqOWU8fAF2uIsWqidQE1CfScOxWteMtL2HDoNounU4pz/vLr1L+auokM1ONFMSB2hjC0sUdiXHWBdSNfLiky+QzY3DLX8DDXwfUw94T9iFsg1qnMwZrKWRplrRs9ASZje1ofFmH3LVyD3RAuJ91plK8qQNa4Y/cq0ji5x/hucVfEHNrtSRyzOSrfP2+Ee482pru11yZWzv4Y3pElq+aL+Ggezt0rIdHvsGdx2ZYLyYJVaaddgL2FFre/zm+4PkY2zvVLZQp1VQw2h4Ym8k9fbqx0Eh2HBKHOHlmihC62nQB5T6dqzAYVA/GqaxksqzOwZ5+Nz63gcsQTp67PWvYP9ZaLOQqJ5mVjYEz0qnIfa0cZ8ZUg5I7ECVMiXJGq7C6XXXbILJMSpOgPwYef2OGdu0fqp/BLoL1DM/a0InHJZjN5iB7hoSrl6wMERFFcuUar/3U3fzpjfsbRVD5KUgcVL9XcyqoOq6KeSquEEMiwXW+I1xauJvt4jSbRr4DKOKIBTxYEjpDXmIBD163wXSuQiarzn+wnuF7d6i02pdeczUlXxc9hSOIwjSsU2TP5JxgdlHdI9Yc5d4T9uF1Gy0q1LIktx6c4mU7+vG41D2wqSeM12XwnI2dTGlyTznk7p+nzEtzbJm55P+tB0Z539f2tfTXueHTd/Oxnx52/l/Q5A7qXNgDRm/Ux9ouRe72fTE3A6ZSN7nq/9zKN+4fbZwCx5ZRx23bKs0B62Zyl1Jy17EZnTEmndmjPSBeQG5/+pF7wOvi+h197FoTJ+Bxsd7dIK3cxHHe9eW9vPlz9/JbX39QVa5VcqSNTqq46RBaOWlfuIdZlZnij0PnRtZrch+dnFbZBKVZ6L8UcpNkZhqR721Go79LOTXBhqC66UdMNRVzKhd1mXoiX6G7SdXYvmeHTEMpRcIMOap+KqRb4QsXZriffKVJuQP0bKOjeBIAz6nb1N/GH4QjPwXAbRaJudSNGRUlZjI5Pvqjx/nS3eo93PKXcM+nSDSlrHUf+RYTspPbrUuZyFZUsc7EI9x9LMmwW5O1nqnYKsuYfhxDSLZ0KdLNlmqa2CV14UE2Z8vkp50BlXIGakVOPa4Iblx2IaYONBqEaSRyZQZD6sk4nZXcOaqOd2ePgRCCoNc1T7kfnMi2EIG3kiJrxPC41KDa3atmcH5ZJolSv4FwjKhRZsauU9DxkmblXsg2kTvwCFu5V+6E9dcAUPHGCcsC2/uDdId91GdHAcmE0UeOABFKnE4WyVfq/OCRMczT96p7DuDk7U3naUpbMoIjsecxLJJc4h7FQvAD8yp85QRIiRCC7QPq+LvDXgTwp/7vYcwcJpXRytVT5A0XKR/aCHfj69vGta4D6nMueqX6Odd3txec0TPgRK7i3H8dQQ+pYpWv3nua0dkik9kyharJJUMx5+3X7+hj74dfzAu29VKomuTKNdLFKoaA3qh/ni1jE2mntmVKtVbyLddNpGyof9OSLWQNKqAa9rtb9tMZ8uJzu1jTGaRatxz7ZG7e+ViqRKpY45dHlO0mpWyxZaSUTmuDYtV0jsMm92ShyvFEnp8+NolbN0v7+aFpwj63Y6Gd12U9l8HTjtyvWNfBZ9+2h6jfgxCCzZ5GJsJaYwqJmq4CZAolqBXIygBZESFGnlqliNDqLFqfhZkjyssWwiH3yZlEQ4EOXQGFBIXZRk+YraIx0pu5KYZ8JaRwMY1+YO2ptk3uWrkDSt1qFdsjMlBKMV0LsrlXKckpo089+LEh7BqVqI4tqDdtw5U8wtqYh8HZ+yE6rDzbuz8BgNcqEnU1FNaRkyMUqmZjWrn/O3DsVqZzFbpCXnpJMZy8m++Y12LiUrnr0UEop9l3bJwhh9zVTW77o6GM6o+/Nq5UY7ZUcwbN02a3Ose1kup4WU5TyiZ52xfuczKLKidUJs/N5h41C5pDNNO5CkN6QiM9Ab59QD3QA371gIW8bsdDtYm4ZkoOTjQCzf5ampKng4sHonhcgs6efue1lIwQ8blxB6Jc3CkI1tVxVSvqwbVVXjzooZpPgTDAq67R79V/i7dU/pi6lmWT9RCGkOzslHSHfbiy6vqPyV5yMkhEFDme0Jk+nMFVScMlN6gDOXlH40vnJhS591zEuGc9HSLHbo5y2urllOzHMCvOfblDWzNdIR9MHuA36t9h88ytpG3lbmZ53xV61hLsxujZSofOEmJgF8TWNKqhNYopRYJSk/tEpowlVR59PODlrmNJPvz9R/nafac5rfPb1zd1ohRCEPV76I+pQWUyUyZdrBELeAh5XfNsl8Jccp/TjttumGZnZ6WLVaSEY4lG2mi+YhLytZJ7v+7maAf77z6uBue5mSt2jv5DI2mklFTqlpO4lSvXOZ4oMJYqcc0WZTvZNTUjs0WetaETgLuOJbnl4DQv3zmgCxPrDMb9GHqW3rZlzgEbPerCSQTDIsH3fvO5/OT3rgWglNfTMctP3ogQlXnc1UZwK1iZUcu49aiAUySmLlgyOdvwjnu2gbRwzTSmgttEQ7n7K0n6PEWkv4OC1C1CdVl6MXGSl/7jLzk5U2j4kYWEImNgnZhCSIszVT9benVgsmKS7H8eydgO52acq9yp5nlr9CECVgFpWwN6tuC3SkSMBrkfOnEK0JkUlqkIpF5hKltmU2+YHcYpDCxuM3fj9xiq6jQ6pPZVniJu6cEzbyt3Lx7qrJFqJuMXJrGARx2rPmcjuvKTUsqZ6gdFhXuPTmLpeEd8RtkPN1t71LZN6Zegcuj7tC3za1dtJY96UIVObQ35XA45ZEt1fG51Kzf77mEzTdUb51WXDvKS7X24w42lHWdlhN6oD3wRoqLMWp960GfTijjsGMLOoRhmMa2CqYZBzbQo1yQmLscOGimrKfjWaJWeiA9fXt0fp6xu8gSJUnSqhd/Yo4RBYtPr9TmabbRwzk2qQOfALia19XOZuZ/Dci0JqRWyDvbbvntX2AvHbgGgXimooi9AIBsZOKFu6N7aOLnxNbDhWmVJ6hjDgbEMn/7xffq9Fp0hr0NMYb9Soja5HZvKO97yuq4gc2GT62S2TKpYJR70EtTk3qxkS1VVkxHyqvt7LvnbjdNscm/+aVefKltGxWDsoKo9uNi57/b3yM8pJrMV+HRO5caXmgabYtXkTp0j/+vPUinKY6kipapJIlfh2i3d9EV9/OVNjzFbqPINCpvBAAAgAElEQVTKSwec2dRALOCQe9uWOQesMWbIGVGyvn7WiWk29YTY3BvG7zGoFhTZpC0/RVeUsMzhqSlyH5ddeApnVCqeJndbmWUys0ohubwQVxc2nG4id0MpMylc9Ig0PfUJRHyYuks95KROqd0VJ5lIFTCEYKc9fc02BobNQnnvM2aQjT0hhFAWw7vy7+Pdhfc7QZ1oC7kr2+Yt2c9jSsGPree0PLgBWSbSpNxl0VYtdfLJcZWqWS8znauwqSfMGqHIYlT2sqU3osg9MqCP74zKHgJnyt4Z8rJBTOAR+kE0lTLLlBrkPuqQexor18goGPDk8dQVKWyrKa/5UWsDNW8MmgZP05LM5Kv0BpRy27NliA+84nL1olauIZ/byXPPlmts6A7RHfbynb1jfOoXxygWcvhlmbq/i3dfu5F/fcsVEOx0PiNFhN6IX13zSg5fVReT5dX+s03k7qpmkdqSaa42tQeAY3lFKr1Gnp6wj2j5DBhuRqsxap4wEdEg9xuihxiXXdw41a8GDIA1z1Y/E4fU4Nu3nTMocvdS45BcQ8KeFWpy3zGojqcr7HPI3aoUyReaiqHsPPpgl3OP1HGp67vj9VDJcvjO/8SyJN97cIwO1Hc3kI37FVpsBoCj03lOJYt4XMLJXGrGXOUeD3rwe13IOQvUF6uqJsPvUeQ8d5U0m5RtUk829dWxz2e+3OS5a+++Tw8ug/EAza3l59oyzUsYPjSSdgYXOzD+4EiasM/Nszcqu3UsVXLes7YrxLuv2chLt/fx2bddwUu39zkD7mDcj6GZ1Xoq2zJCiC8KIaaFEI8u8vofCSEe1v8eFUKYQojOhbY9HxgkwZTRS9IzwAZXArcO7kT9HqpFdbOmTD9ld4ywmcVTU8rmmBxW6gaayD2ERFApZKkVM+rhi6ipfG/xKDXcyEAnnSKvtouuo1tkiOWOIfp2EInqB0LbMm5M3nt5iMf/+np+/dm6QElbNZXwsJPHnJJh+mN+wj432XKN06kK04W6U0jhBFQB+naAN0IwHOcz0d/lT388QrH/SgDM7osJiRJB0XgInKk4kDyjsjCsWolMqcZgzM9G9wxl6SFBjG39EaayZaywIvedRlP+ubZlfG4X79zalE9u1Yj63YrotCpvVu7NsYpndTQepi6Ro+YKkCNALrKxQUSoh9m0JN0+TQSeIC+/Yov63V58xdtQ7plSjY3eFF8OfYKPJH6X2q0f5Z5HVIm+bCJ03D7QWTfNyp1CAkPqgSKXh5N38K67X8BXff+XZ7sOE6ZIxa0G/uZAqx1gfDylyEUUpumOeOmqTSBjw2SqFviiRChyLJHHS42uqbvZH3g2//XIhApcA/TvBLcfjv8CgLuyvYxZXc7nHLTWUvTo75FXg+XWvjAffuV2Xn1RyAnQinqZTLap/mHmsIoTuDzQtRmAM7ILKQzYeB01XweHbv0yX77nFD/cP0GXUM+LC4tLh+eQe0AR57VbexhNFTk0mWVNZxDXnIU5oEGuk5kyU9kyXSEvAYfAG+q8WDUJeF0EdPbTXOVe14FdO1e+ueujbc3YqZAAncFWW8brNpzBpzPkndcG4nSyyMbuED63wUMjKefz7eN/cCTFpp4QHUGVoTWebiL3ziDvumYjn33bHq7f0Y8QwrHKBpuU+1Pdc/8S8LLFXpRS/r2UcreUcjfwp8AvpZQraLO3Oui1phmXPUwa/QyLRp5zNODB1DnXqbqPiidGwMwR1It2THrXNXbSYwcxBaY7SJgSucws+KOOiu2vj1FwxRBhRVxFd5yyv4+tYkwFunq3E4vOV+fP7680smTAyZSp9O7CEOrCp2WYnoiPqN/DZKasp51V52ZssWUCHfChU4jffoBfefuHyFfqfMf3enjlP1Ht2ESICiHRCJZ2ihyDWknlppR1Y1YVOfdGfawzkozJHkCwrS9C3ZIkXYpYdokmcs83zu2vrmvK4zdV5WS2VJ9ny5RyM6QTjSyj3ZHWbBYZ7gcEqcCGFnK3M3m6ffph1O0H1InTtkyT5x7Kn+YfEu9jR3EvO9d08HvuG4kf+CIAItjd8pm2ek8RUQ+xLwI0HsBcoQCJQwTMHHvEIa4c+wpRUSRtqSn+XOUupeS22TgFVxz2fYmesI8hkaAWWUO+XMcVjBMSFcaSOZ7vP4qoFXBf9DIen8iSCw6rHcXXQbgPeUZZVV88GmCsHlUqGzgk11L1a0tJK3chBO+8egM90/eqdQCAgKgwNdsUyJ45Bvb3j6+lbngZs7pVINPlYbz/xbzY2Mc//OhhZvIVNgWV7SIWUO4vvLiX118+xK/uWYOUKkVw/SIrP/k9LjqCHo5M5zmWyHPpcLxB7rVmcq8T8rqbiL/Vc59ry9jKXQjVGkJKSb5aJ6LJvUN77gOxxgpKazuDdAQ9bOkNz7NlRmeLbOwJc8lQjIdH087A0xdVyn0sVWJTbxghBEPxAOOpVnKfi91r1OxqXXeoyXNf8BQ9KViW3KWUtwMrJes3A984pyM6G0hJV22S01Y3Y/TQJVNOV8ao341V1svC1byUfV2EarNO6XkuoisNvWHHYwYQvgghypRyKUUooV5A4MKi7O2EkHrIcq44RW8XmwytTHsvpiOuvVNpURfqhrsoMCf3OjsO7gB0b3H+lCbMmo4g0YCHx3QRRqlmOo25WpQ7gMutAsDdIdZ2Brk/2wl7foOaK0hIlAg0kXsHOV50sQoUl2eV32vptMbeiJ9hMc2o7MHnNhz/dLLkpuaONJS7MJRyN+sqF9pO3wMw6/NsmTEUuU9OTlKYbcoy8imLqIx6eFwxlb0y7V+nvHldpm/n4Mc9+slwB9R39gSdzwj63E62xZ7SXQRkEd7zC4x3/owcIdZNqxRXT7Ths6s3qoGrq3eQK9d3anJvIF8oOJk7+907CWSP0eUqMVlRx9xclp4p1TiTKZMouzi45d1w4jYuKj/MGpGgGBwmX6njCSqSDFhFrvc+Am4/u659FYaAQxVNvB3rIDKAkBZZGeT+WT+5KmQ8PVQNPyOyFxHoAMPtKHcHJ34BviiF6CaCVJhJNQ2g9ZLy2wEMF2f6X8Tt1qUOWR7sfQUhUeEfXJ8k7oMtYXXeXVhcrP1jUJ779Tv6+fibdjt1DZW6taDfbqMv6ufnB6eQEvas71hQndvK3WUIvG6D4pxsGTugarc5tnvpbOuLcDyRV2nNEke528H+viZyf+91m/jzX9lOxO9paaUgpcq8WdcVZPeaOAfGM45Y6GtaXs+OhQ13BBhLlRidLRLxqeyhudjSF+HG9z+XV1zS79hBT2lbZqUQQgRRCv8cl5M/CxRm8MoKJ2pdnNC9uW1LJBrwOP5souYlH1yDIetcZCiCq3eqaSrdW6FJWRv+CHF3hVohrZS7y42lCb3u73TIIeuKU/A0Tfl7t9PTEcOSal9HpbJh3Lk5LXozYxAbwhNtZG709PSxpjNI1O9u8QFPzSiF3KLc52Bzb8TptV4xgoQoE5DqITWFiw6R59qtPXhcApnWMwqdLTQQ9zMgpxmVvcQCjSyHiUyJgr+XLl2cRdcW5bnf+yn42GY4eotjb2BWiToB1TQWgs5BNXAmpiepZhpkZM+sHrPUuTGiA7gMwRm3nkVp9W6nacbcmkg92tfVPd0Bwj6X8zBeWjtAwr9O2WuGwbHATrpN9Vm+WF/rCdPK/VPvegkv2d7nxFlsFItFJ8g44tuCSJ2m30gzUvRgWbKlz3emVOOQLpoRe94JkUGuuP8P6BEZkp5+aqbEHVRqLiKKPNd6ENZfQ29nJ1v7Ijxc7lcDZ9dmcrqK+qSxllzZZDxdYsa3jsnQRUgMYiGfEhqF6Zbj5fTdyrP3R1XRV31Od8immcvRa/6Zz5ivdsj9sO8S/qL2dq537eW7w9924g4uIVsKkpwKaWBdV8hJ+1tMuYNSz4WqidsQ7F4TX9SWsQvSAh7XvK6RNdtz1ymIs4UKEb+b7QNRjk3nHQ/dJvfnburmAy/awrM3NJ7L67b2cMMVw0T9btXLSSORV4VtazuDDMUDVOqWI6aa1061s9iGOgKqgd5MgTWdwdbZeBMuX9uB26XSdYV46tsyK8WrgLuWsmSEEO8RQuwVQuxNJFah54Im8hGrm0cLrX531O9xMitmqj5KYUUol2o16u67WG1vWzL2MfrCdLoreOp5xwqoBvTAEexylHtaxMi6lFK3/HGI9DMQD1LUqvSk2U3FE2ssyWYjMwbRIbzxBrk/7xIV7GpJeQTVD4WlyX1rX5iTMwVqpkXFCBCijF8XB1UDfXQbOXatiTEQC+AuqBRMOxV0bbBGWOYZlT3EAh6n981soUpWL+snXT41y8gnFJF4Qioou0HleDdsmRqylCZHiDUDA9QxyKQSkJ+mjM4/rioV/5i1Xh1HpJ+I381pQ9sTOqhqL7gccS1O7kGvUu5mvcZlHORMfI9zTpJdVzi/h+JzyV172XaGiq+V3EvFAtRLmBjMBDcBkrCZIVHzc2Q615JOlynWOKR7/GwZ7oY3fIHK8HO43dzJw4HnAGAE1H3ZS5qB+jisVcHTLX0RvlZ4Fvzm3RAd5HRFDZbhNTvVtatb/HDjh7l5+98CuvVAuHdOa4ykCsKuuwq3L0RANKp6Gyeg4d3btsVs0U4trHGj51fgqt9m85kfOLMCj7Dwe1xOoDLcdP953QbrdVfOpZS7LRR2DMUIet2Ocp/JV/jANx8ikatoclf7Dnpd84qcTMeWUYNtslClK+Rlc1+YiUzZWU/Afj4CXhe//5KtToC2GWG7xbbGaJO9YhdR2enCfU0Dm0Pu8SCZUo3bjyScNMjlYAiB+Qwh919jGUtGSvk5KeUeKeWenp6epTZdGZzVinoYregbTWeHRANuXDUdUSdAPboegEvESYrSR1/fAOz69Ua+sQ1fhLAo4zMb5J7U8WFvtMch91nipA2lymTPdhCqeVFJk3tSRrGiw6rS1Ua9okrLY8O4tHefkwFetkuRW3SO/XJypoDPbeBzz79ZbWzVPvmpmQJlEcAtLAK1DHiC+OP9vGKTl96In8G4n2BZPbyGVaE34iNYGHfOXzzocVLJUsUaGbdWfOFe9a8wDZMHYNvL4U9G4XkfUK9bNaJ+VS1ZK6bJyCAdIS8lV4TM7DSu8gwTbmV7eXPqXDwm1bUgMkDU72HU7FSWS+IIddPixgfHVbWnpQO3Hn1tfdEmz10VMRVP7SUsyiR7nuWcE3ONqsCsS4NoxxzPvXsrdG5UQUZ7n00oV0pQK1PBSza00fl7Vga5+1iyhSDSpRqHJ9UiJlG/B9Y9F9evfo3/r/6n/DKjBhWXtmU2GbpOQluAW3rDjGSqFOPKnjtRVuTeueEyZ/8i3IvQAf1YwCb3Jltm5B71c93z8PhDBKjiRw8+9iDWpNxt28K2N1LFqrrmV/22akkgLXAHMHQMoluTXsjXev9t0YS3lHK3rY0rdU8iW7nfdWyG/3r4DPeeSDqpkPbriwVUZwu2cq/SGfI6n2u3GrBTKZdC2OcmX24sqm7n6a/tCjpdKW1yt5W712WwRleabuhW9+DrLxvif/3Kxct+HoAhngGpkEKIGHAd8F+rsb8VQ6v0cdlNSmoFphcqjvo9uosgFPAjYgOYLh9hUSZLkDWdQXjdp2HLi1v36VWee8AqgD/KdLbMPQlFBN19Q46HOSOjpDTpG3rBhsFYgKJUN0paxPB1rVPHOLFfecV3/Ys6vktucAaJghFxmm/Z3R9tJTIyW5zvt8+BrSyOTucpCXUj+qtJ8AQRwS48FTXVHowH6Kgp1eeWdTZ0+p3B0VbuAa8Lny4zTxqKHESkX9kBxaSKF/TvBI9fZZ4AmDVHEZazs2RlkM6QFxHowFfL4CklyfsGwOVD6FnMUe92pOGGnouIBtxkKpbK5pg5zM8en2I8XeKdV2/QLQ9E47N8EcdqC/rcSAnlo6pnfWHgKuecdGzaQ1H6SBGmM9yYYgPwvN+D993VdL31fWOoc14pqeZbZTyUYxuUbQK4gnHuPznreO69EZ9eIzXPpt6G+g94XWzsCXPfSSUyvEFFblvt2oioijNs6Q0jZaMl7YGsunbRdTvx6pz9oM/t5HDHAt75yn3kHnD5YPAyDG+IsFHFb2dKRdTnOJ47DeVuZ5+kijXlHUcHGiIn3KtmZqg0S+8C4sJu2jfUMT8N0oYd1NyzvtM5L9BIYUwVqxSabRmva14Fqx1QTTXlt3eGGr10jkypfYWXmNnaCPvd1C3ppGKenCkghGrl0FDuivB7Ij6EUOsG2Nl3L7q4j2++5zl87I27nHYLy8EQ4qntuQshvgHcA2wTQowJId4phHifEOJ9TZu9DrhZSnl2zZPPFekRar4OCgTIEcTCcIJy0YCHAGVMVwATF/2xIBVtzWRkyBmR58EXJiQLBGQJfFH+z38fYtJSCt0INWyZGRllBqXKhLZ4BuJ+x5ZxR3sx4mtU8PGz18A/XQp3fAy2vxY2vwh0QU0w1nj4bOW+S7dwrZmyNcd9AWzuDSOEWgaviHqgPKUZ8AaVv6xnMmtiHrpkCulW22zs9DjFVqOy17GEOoJeUsUq08JW7n3OsQKK3AEMPeiYVWd6XsolycogHUEv/t5NXGqcpFtkVMzCHwNpYgk3b3/1SxF/cAg2v4ioX1k69GyDxBG+eNdJ1nYGefHFfY31U21/s8mWCWlSMEbu4og11OKtb+zv4AFrG5Oyk47QAsFob5OdYAdUNRnWqyWsaomy9BIKhpx0xVC0k5HZIrmyaqncGfKSLtY4PVNkwxx74pLBqLOQtzesruUWXdNgf05jUM6RKdb4VvYS7l3/flzrnst6vb+Q1+X4ycqW6VOxj9wUnLoLTt0Bw3vU4OcJEtKrjZmGt2HHNCn3iM+NxyWcrJO0LjAC1Ews2KW6XIJKVgh5nUyUZrzj6g38/IPXLUly12zp4dW7BrlaV3fayt1OYUzmq5SqJgGP23l9ri3jKPdiVS8GrmyZQYfc7fjL8uRuf49cuU4iV+Er95zmWes78XtcDrmPa+Ue9roJ+9zONQLwuAyes7ELY4HUz8VgCIH1VK5QlVK+WUo5IKX0SCmHpZRfkFJ+Rkr5maZtviSl/LXze6gaJ+9QDb0KM4rcI8rSkBjUffEW5R6hSMWlpnB9UT/1mArcFUTIKVWeB2+YmDmrpqb+KI+Mpuns1znqwW4n7/2MGeWEsZ6Pu94JO9/kfGZFKPKMdQ/Apb8Ku98Cr/oX1cbAF4Hr/7falz8OhodoR4M4bYLd3Bt2btil/HZQaWdrO4Mcnco7FbLuUkJ548EuZ7Db5C/gEpJqVJ2DjXEXpEeoukJkCDntk+NBD6lijUmpPelwnzOgAarXDqgCLwCz5qyMZBXTZFDn1r3j1awTU/SKNK5ILwT0ABmI8ZrLhtWAocvVs+UaRAaQhWkeOJXijVcMq/xpexUm58vGGuTuc+OmTnT6Ae6xtrfEK7pCXj7iej8fEr+/pKUFNDz3qEp59VCjVCpQkl41k+pWNRCBSCfj6RLZUo2I30Ms4OHkTJ5cpc66OfaEXWAE4A+r87jF7kekP8cOTB6dyrN/PE2WEPXnfRBcbud8hnzuBrkHNLlbdfja6+FLr1CrWK3VMxZPAL/23KXb37Blmjx3IYQavAsNz93J+ujbDn90vHF9LZMXXNSrgs5z4HEZTo+hxTAYD/Avb77MuY9t5W5bH7OFqkqF9DWUe6lm8rV7T/Omzyq7yVbu1bpFoWqSKlTpDHvpCnnxug0OTzXuheVgz4Bz5Rp/ddNjlKomH33dJUCjsnVMV+D6vQYffuV23nPtxoV3tkI8I2yZJxXSVM2WJg9AegQZW+u8ZPo7Wjz3bcYo025Fxv0xP1bHBgCqnuii0W58YdxS+5a+KOlSjVpIT3HDfTB8Jd9Y99fcVt9Fvmbx49BrVFaNhqX94d7+YRi+Al77r3DF2+FtN8IfHoWYTrsUQk2Bg42Hz1bpwx2NqeJytgyodK2j0znymtyNQqKh3Kt5qFfYEdZtiaU6H+tiitwLwSFAtJB7plhjzNQVkbYtA8ovtsnC9qxNtZhDn69Gd3WM07JPTf8vehWWXjw6EO93mm45DbPs7xxwqxx5XwRRVz3cd3gn4N5Pa3Kfo7K1LdMV9rFTnMRjlrjX2t4SrxBCEO5dSz60jmVhe+6hHkzDh486xUKBCqrbot1mORTrIlOqMZktEw24iQU8HNeWik3GNnYMNe6HYFSR+5BIqrbFeqZgByaPTuedRlh2bvnGHjXgBLWChCbPHWDqUbjy3XDNB2HPO9TfvEH8UpG78AQW9NyBljVwU83KXZ24xixJWrz5WWv52xsuXf4crgC2crddCrtvTaDZc6+a7D01ywF9PupNSeKnkwXqlppNGIZgMOZ31ipYiXK3tzmeKPDD/RO8+9oNbNZpjh6XoXoI1S0Mobz2N+1Zw6418aV2uSwM4yluyzzl0Kdtgcn9kB7B1dF4gEWw01GqcXeNneIk+9iO120ohdKpyF36YvN268DbyHuuukNkSjUSfc+DN3xRpZwJwYnel5CtyZamRTakR2cSrFmAWOYOKK/8R7j6953/2upzTWfQCX45qzAtgS19YU4kCqRN9R5RzSlSDOuMnEdvZOOpbwLwi5ReuCRqQGlWpXdCY+ETbcscr/eSNyKqyZRNKrYlA03kXkUIwasjR/BQ5zZrtzrXoS6k7pq4ft36Bqn7W8+9o9y19x2ixEVTP4Kf/Iny+JuVu23LSMnFAxGuMlQV6n3WxfPO01ufvY63PHsty8L23INdSLcXH1XKpQJlPGrA6FGWW7RLKdjDkzlHuduYmzWyY6DxHSOhEDWhVa72221s6Q1zbDrPI6Np1ncFiWkV3VDuLnYOxXjrc9aqxR+aB9mXfgRe9P83xIIniEdWCYpKK7mHWhMXOkNKuddNi1y53tJWAFCBVXB899XC3AwW298OelqV+3Su4izO3dx0y/bq7Rl3s9+/InLXwslezOPK9a0ZL/Z+g1734sLvLGEIsdAKkk8ann7kHupSvuWJ26BewtPdIFFXqMvpI96f3Y9HmNxc2ExfVPVSd3epaZahg1wLoqmoZbYewLQksWBABZx0w4ig1025ZpEr1xzv1zkGnyb3tStQjVuvb3icqEDV1Zu72bOuw5n2NucYL4ZNPWHqluRw88p23hBc8npY9zz4/vsQB77Nz3rfwf6qIoOhsAHVAkKTW0O5e0kVa0xXPfz51pvUMYZ7VWBxYFfTF9WKTzdBe4HrIbIyyF5rqxO4c13yOvUz2qTcA3OVu4di1aSuFXqICmGhc7XPPDSf3KUF1QK9ET/XeA9xyFrDLNF5aaQ3XDHMe6/btOy5w637B/XtQLj9eKlTqxSpSK/a547Xwes/T3StGtimcyrX2iZFlyFUcL4JsaCHNZ163VufC6FndsYC5H5ypsDNj0/xnI2NGdxVG7vYPhBlW18Ev8fFR167Uynszg2KfK/7kApqN0Ofp+cOGhiegFpIZOvL1GyzCR16DVx7bdOO4Bx7UgeQkatbWulzGy19Xmx/O+hreO6lqslUtkzdkkgpqZnSETlzyX0wpr6vyxD4PcvTmD0AHJ5S5D7c0XrNunXGzEJplE8UhriwXSGXH/Keiui/BI7/HABXx3rCPkGxWscd7oIp1Tq2K3EfdWlwZ2UTOwbUjRDoV1Psvt7exffdlPc8UfEANUdR2bAj/Ml8lW39rRWOm4d6kRkDb7iLs0Vf1M/X3qXyoLsdW2b5S2QvGH1guqnCzxNUBP+W78JNvwvd25Bdb6Xy9c+r/brrUC3giajva/fj6Ah6SBer1C13wxLyReBt/wkDuxv7b7JlkJJLi/dzm7UTDE8jCLfr11UWytqr4LH/VH+bp9zVtmURIAyqN47UhVyVrKpOteG0IFBN3S7nMF+3rlMO1wrS4RbF7zwEwsC461/wiRpWtUiZEN1+jyLRS9/EYKZRHNSs3Ic7AgsGFncMxJjKVpTnH4xBKdFSCQ3w/It6+cljk7z2siHeftV65+9rOoP8+APXzD/O6CB88HBrgNuGHhy7jTy4/LDuKvVvDrpCXmaLVacvzjzlbmhyW2hVs3OAEIKAx+UUNtnVos3ZMqWa6aRD1kxJ3bLoifhIFqoOudtpi86i8V7XipS2/RzZLaGH5jQ8s23QgHf19K7rAtsyT1Ny3wlHb1a/x9cSC0zicxuIYJfy3KUkPHkfj8oNFAg45ciujnWw4VrWX/7ixffdVLE4VtLkHliY3KdzZWdtURvhrdeClXNU/hOFfROvxHPf1K2OeazoAlvQebUP7A3CDf8GwPMqdb7u0hZBvQK1IrFojB/+ztVO06N40EPdkqSLtdaBZePzWz+0KVuGyf2EazPcZt1AR8jbeNjcXrjsrep3x5aZr9wBCihy7/HVcdeaFjVuVu524VFJrZTlk2XutS4m4nOfVRbDPLh04y+Pn7CrrrJl6GixenojftyGoG6pDKaYVrxzg6k23vKctWzq1a/Zg5IOptq4fG0HN//+dWd3rAsROzRiE8UUxIYXfXtHUGX5zOSrzv9bcJ6UOygCL1RNtvZFeFznqDfnuecrdcfGqJkWdVMy3Bng0GSOR/QaDXbnR5ucV/J8NG93KlmgO+x1vH4bji3jWT1KFEJc0IDq05fcbcTXEPHPqIsV7ASzAsVZ3BMPca91vdpcNwLC5Ya337T0vpuU++l8U6ZCE+yqunLNmlfgwa5fU//OEV1nodxjQQ9dIS/5QtNU3TufdEI+N9dfuh4eQwUrqwXwhlpW02kOsC354DRly3BGLfpwr3Wx05lvHpyA6nzPHSBv+ekDhoNmY6EUaA2o2j5yYcap/D0o182zZJ4wXD72DIdwz4CMdrc0h3IZgv6Yn7FUqUW5z02DtHHNlh6u2aKJ2P7Oc2yZVYU9CJZmoXvzopvZ95Xd2mJxz331yd22PHYOxZrIvWHLNIvcmmlRsyw6dUfGM5kyF/VHnIvzYAoAACAASURBVHa8tuc+7/lbBPZ2UsJQx/xrZtugfu/q2jIXMhXy6UnudlA1oBo/DcUDqsw3oIMkp25HWDUOCJXG1twIaFk0BVRP5JSKic8hrGDTDbCSNKwnAsdzXwG5g7JmHihUsDAwsFpJsQlvfu4WRe71skPuzWhWcksWhzTbMrqdQV4GGFqgoRLQ8NrneO7298vp4q+BQN1JdwRalbtdkFNMOouAzMgY61eo3paF20e3X4LbJD7cMy8APhQPaHJ3O+S+mHJvgZ1NFTmf5K6vdzXfamXNgT0jtNsmPJnK3X5uLhmO8a29oy1/m6ukq6aFaUo8LoNvvecq3C7BRf0RZ1bo2DIrfP58bhdet0G1bjnrmzbDtkEDK/DvV4oLXcT09CT3zg0qj1svpPH3b9ylyopHdWn2aZUnO+lbD7VGn4sVQSv3ujQ4nlYXZq66ab4Rz8nrXQLdoZWnQgJs7A7zwKkUZREgKAuthTrN0EVMqruinDcINHe7W7KASgjlp5tVh9wDgeDi9QOLKXdNkpm67ujnq6uFohH6+Jquna3ci0koJJDuAIYvtKKMohXB7Vd2Vb20IEHaVkA04GFtZxBDwK41S2Re2XBsmSdBucP8YGsT7OP96WOTwEKeu73KxOp67qDUeUfQw3CT370YuddNSc2SuF2CncPzz7FdAbuSTBkbEZ+bZL3a8vk2mrNlVgtG25Z5AjBcKotDPywOodgP/8g94PKSCwxDvtTS33lZ6GyZPAGn4GKu596sFoLnSblftraDt1+1TqXArQB2ULXiChCsF9TgtxDsUn5dDzBXubfaMst8N5dXZcvoLop//prdDHRGFt52Gc99RpN7t7cG2Rz0XgzTj7cOPrbnXpyFwgwi1M1r1g+pFZVWA26vahldKzfOUxNstRj1q0Kjhz780nnB9gXxZNgyzddxCeU+3KFW/TqRKOA2xHxyPK/K3c1QR6BFADTbMs1QnruFe5HYld/jojvsPTty97tJFqoLKnd7RjP3OM4FhnFhW/4+Pckd4I3/Pv9vti0z9Sh0byMsfEDpLG0ZpdxzMkiyUMXnNualRzXfAOEVen5ni4DXxV+95pIVb28XvlQNTYbLKXddDzDflmmQ1bKzBsPTsGWEi1/ZvUReee/FqgNnczolKksnFvBwy7ECNwCdnqpS7ltfNp/cXR5FlFq5E+rmo6/byarB7Vf7NiutSljD9nntQW9FxA4weJmq/AyefQbVitF8vAsMTM24dksPJxIF4kHv/EyT85TnDvCH12/DknIOubtaftqwA6pu1+KB8j95+cUMxlf+bIedIsGFPHd1TKubCtm2ZVYP9pJq0oKerUQL6uE7K2WnlXtBqBtg3rSVJ8dzP1vYyr3uDkKVxZW7TQK2cp9jyzTPUpZX7h5ty7gbg8ZiCPfCb903fxeG4OrN3fzksQksj6DDKEGtoFTu6z6neqc0w86IKsw0iqtWCy4vlHUwd4HvY3dDHIgtrowXxM43qH/nEy22zNLHd82Wbr5096kF7+3zqdztzLJC06IZth0zl1SrdUnNspbsX/OGKxbPCloItspfWLnbtszqkbvrAtsyT78ipqUQaEpL7N5GXPcot7vsrQguD7h8lAxFjvbakc1oJvSnCrmv7QwS9LqwPDrbZ1HlvrQt43YZDqkvO+V1ebVyryyrFpfCtVu7MS3VvTNu6RmFLwK7fhW65hQiNZP73CX0zhVuv7PS00IEuWd9J7f94fNbVil6yqB5MF9moH32xi7chlhwNaHzlefejKBXBTddhsCryXuuHVK3LExLLrhG6xNFWBcELtTNMh704jIEwVWciQvRtmVWDy4P+GJQyUDPNn5z+2Zee9nQ8u+bC1+YSlU9LAtNvVsCqk8Rcve4DH7w21czcPO/Q4ZFs2Ww89x1Je9C23UE1WLCy9oyLrcid1heuS+Ba7fq9sf4iVX0Ii5zVkhyEOxWa9RqW2ZV4fY20jAX+T7ru1eQHXMhsMKAKqh79tW7B53CtRacR+XufIQQKnW3UndsoWbvvVQzqdYtaqbEs4rkHg966A57FwyaugzBJ998WUta8LniQneFfGow02oi2KHIvXsL2/oj8ypIV4ToIOmMmvLPDaZCox8GrGyhgCcLm3vD/L/2zjxIkqu+859fVVd3Tx/VPTPdc1+6RhpJSIPUCLQSIIEsCa1ACLNYQgZxrVZr2F1wEAYbB8KB8a5W2GF7MSjGWCHwsgKDxBGAbBa8MFjoYNA5OtExmmnN0a05+pzume5++8fLV5lVXWdXdmVmze8T0ZE1ma8yf53V881ffd/vvccS7/ctUucOWEFOtZT03MH67rsPVZu5H/MWeShRJVMFq3uWsHllFxOH21kx5c1X3lbic+tYbleEmp1eBHFv90WtjodVJOR57pVto796z9biBxaxzj3I0o7WvKzWjQxd09vOC8MTTB2312+pcu70avjYpaeWtXLe9prVJY8thKhHqDaXLQNep6rYdT8Xyg13872+/wjMH8AE9g/OWT3VDqJoGE6sS2Xu4HccBtsH6O2wVQgVvxLnqmWm6hbDt5+zhpmWTlJjtkSvtLgvsw9vmDcpVt2kA9ZShew3dqQz/qjhemIPzAq5mCzrzM+gnefuBhi5aQjKdajWyqa+zrw5fBabqEeoNp+4d/ZB7/rSnnM1dK8k4y1sXLTTCb/jJS6ee47WCp475It7kYdAf3fb/EUuipHybJnZY3V57gAfe8upbF6/Kjcff9nM3RG2uAd/hyqy39jhPvN6Ym+A5w5wzdY1eVl0X1cbp6/s5t95pb+Tx2yna6bOaTyiREeohs2bP+1ndnXgOhULR6c6OltbODJ5PDaeew4n7qWqZcCK+9yrXvv57T7xO5s5OD5d+VquQ9XM1p25i0j+YtXViHvYpYUtCc7cwT6op0bqzNwX33MH+A8D6/P+3Z5J8y+feBNP7R3lf9z7TG7JvTA7VBtN7EshReQO4GpgyBhTtPBaRC4B/hrIAK8aY2qcDSlE1p0fymmcuJeat2RJa5p0SmirpRKnEWTX2P/k7WUqOoIiVkTc1/YumTdrXlGc5z43U3fmbmPpKv46iGbupXG+ez2xL2KdezVkPBvGLbmXCdGWaTR2sY7orl9N2nkn8CXg68UOikgv8GXgSmPMbhEJufg4GtyEVsU8d7C2TLXTjTaUre+FU95SvtbZZdmS9icAWwjpFpg5Zjs3S4lxLQQfNKUy92Anatgdqkn23MH/tpaAzL0Urq7d99xjljzVQCriUshq1lDdDhwq0+S9wD3GmN1e+6EybRODb8uUFvfYWTJgO9Z615dv4zLU1q75q0PVdC3XoVq/5w7kC3olW6a1q+JgnZrJy9yTKO4hZO45zz0icfe+CR/1MvcWtWUWfv0QzrEZWCoiPxeR34jI+0M4Z+TkhpiXyNy7AosXJw4nAvV0OoNvy4RQLWPj8bL/lnZ/1slC3CjksLN2d91ir5OC+1wTnbkX2jLJzdzTIlE9I4FwOlRbgPOBtwJLgPtF5AFjzHOFDUXkJuAmgA0bqljfMkLedFo/H7n4pJKjEf/g0lMZ8ZYqSxwuQy1XLlkNrlpmZjokcXcLjJSxeNp6rJ0U9uhUKOhQTaLnHkK1TE7cI/LcU/m2TJI7VJthhOogthN1ApgQke3AucA8cTfGbAO2AQwMDETY1VCZ5V1t/OnVZ5Y8ft6GMuuwxh0nxKFk7q7OPQxbxhP1UpYM2Kn2OpaF35kKybdlcqWQdXwWUWfuni0z1QwdqiLMRJi6h/Gd5/vAG0WkRUQ6gNcDT4dwXmWxCHru9ZDOhDK3TA63UEo5cQfYegOceU391ysknfTMfUn+diE0qM69FIW2TKkpf5NAOu7VMiJyF3AJ0Ccig8At2JJHjDG3G2OeFpF/Bh4H5oCvGmN2Ll7ISt24r+312jJuVsjZsMTdrTlaQdx/58/qv1Yx3O8g6dKef5zJ2TIJ9tw9MZ9chBGqjSb2towx5voq2twG3BZKRMrik8vcw7Blpr0RqiHYGNXYMouJuy9JtGQgnMw94jr3VEpIp4SjboRqgjtUo544LLl3Tlk4Oc+9TlsmlfGWxCPcQUxh1MwvhFxHc0LFfc1rYd0F9Y1diDhzB2vNNEOHauxtGaUJCataJp2xde4QbilkVJm789yTODoV4Kxr7U89RFznDjZbP9oUHaoxH8SkNCFhVssUe71Q4mLLJDVzD4MGzQpZjtaAuCe5Q1VEmFVbRmkoTrzKTS5WDcFOx2bI3FsSnrmHQcSeO9jMvRk6VFMCESbuasuckOQy9zDFPQTPvT0Ll30Otryj/nMtBHdfTujMPXrPvSUtHJ1sghGqES/WoeJ+IhJmtUzunCEJ4sWfCOc8CyHp1TJhEHGdO1hbZnrGPlyS3KEqTTC3jJI0WkKyZVIh2zJRk1Zxj0PmHszWk71Yh67EpDSaXOYepi0TQodq1KQzgJzgtkxj1lAtR6bFz9aT7LmntVpGaTgtYc0K2WSZu4h98J3QHarRZ+7BCpkki3szTPmrJA03XW5nneuq5HnuIXSoxoGWthM7c3fCGrHn7kh6KWTSp/xVksa618HN98GqoqsmVk+zZe4Ar30fbLgw6iiiIwaZe9PYMqmYzy2jNCEi9Qs7FHSoNknmfsUXoo4gWmJS5557neDMXW0ZJbnkjVBtEnE/0YlB5t4snrsdoRrd9VXclYUT9iAmJXriUOcetGUSXOduR6hq5q4kkWb03E90YpC5O1smnRKkngXcIybqEaoq7srCacZqmROdONS5e+Ke5KwdrOeuE4cpycRl7qmM/3VeSTYxmBXSiXuS55UBeyujnDis4t0TkTtEZEhEii6dJyKXiMiIiDzq/Xw2/DCVWOKqZTRrbx5i4Lm7OdyT3JkKkI64WqaaUsg7gS8BXy/T5pfGmKtDiUhJDs6WUXFvHmLkuSd5ABPYJQNjPbeMMWY7cKgBsShJw9ky2pnaPMSozj3pnrsIzDZBh+qFIvKYiNwrImeVaiQiN4nIDhHZMTw8HNKllchIqy3TdMQgc29tIlsm6aWQDwMbjTHnAv8L+F6phsaYbcaYAWPMQH9/fwiXViLF2TI6gKl5iMEaqi1N0qGa+Cl/jTGjxphx7/WPgYyI9NUdmRJ/tEO1+YhB5t4stkxKSHYppIisEm+kgYhc4J3zYL3nVRKAeu7NR07c41Atk/DM3Xs4RWXNVKyWEZG7gEuAPhEZBG4BMgDGmNuBdwP/WURmgKPAdSZKo0lpHFot03yIABKt597SLJm7jX/O2IU7Gk1FcTfGXF/h+JewpZLKiYZm7s1JKh1xnbsn7gnvUHXPptk5E8lasMn+3qNESyptv8Y3wxJ7io+kIp4V0gphkqf7BXLz4kQ1kCnZd0+JnnSrZu7NhqQj9dxztkzCM/d0znOP5voq7kp9pDLquTcbkop0UhTflkm2PDknRjN3JZl0LIMOrXxtKuLiuTdJh2pUo1R1mT2lPj7wQ2jvjToKJUwk2moZZ8c0i7hHdStV3JX66N0QdQRK2ETtuTfNCFW7VVtGUZR4EHG1TNOUQqaitWVU3BVFySdyz92KYhS14WGS0lJIRVFiRUwy96TXuec8dy2FVBQlFkg6FuKeeFsmMEI1kutHclVFUeJL5Jm7N0I18R2qassoihInUimtcw+BlI5QVRQlVkScubvpB9JNYsto5q4oSjyIuM69WSYOy41QVc9dUZRYELXn3iQThzlbJqrFmFTcFUXJJ+I699Z0ChFoz6QjiyEMnC0T25WYFEU5wYh4Vsj2TJrbf/98zt+4NLIYwiDqicMqZu4icoeIDInIzgrtXicisyLy7vDCUxSl4UgqUs8d4IqzVtHXleyppHOlkBE5XNXYMncCV5ZrICJp4FbgX0KISVGUKInYc28WYl8tY4zZDhyq0Oy/AHcDQ2EEpShKhETsuTcLiR/EJCJrgWuB26toe5OI7BCRHcPDw/VeWlGUxUAz91BIN0G1zF8DnzKmsklnjNlmjBkwxgz09/eHcGlFUUIn4jr3ZkEitmXCqJYZAL7prfTdB1wlIjPGmO+FcG5FURpNxNUyzYI/K2RCxd0Yc5J7LSJ3Aj9UYVeUBKOeeyj4I1SjuX5FcReRu4BLgD4RGQRuATIAxpiKPruiKAkj4jVUmwU3e0JsbRljzPXVnswY84G6olEUJXrUcw+FxFfLKIrSZGi1TCgkYRCToignEuq5h0I6YltGxV1RlHw0cw8FUVtGUZRYEfEaqs2Ceu6KosQLrZYJhbR67oqixAr13EMh6hGqKu6KouSjnnso+LZMRNeP5rKKosQWrXMPhagHMam4K4qSj2buoZDWDlVFUWJFKh1dL2ATIWrLKIoSKzRzD4XcSkwRqbuKu6Io+cRgDdVmwF+sQ8VdUZQ4oJl7KGi1jKIo8ULr3ENB1JZRFCVWaOYeCmrLKIoSL7TOPRRib8uIyB0iMiQiO0scv0ZEHheRR0Vkh4hcHH6YiqI0DF1DNRSSMP3AncCVZY7/DDjXGLMV+BDw1RDiUhQlKtRzD4XYD2IyxmwHDpU5Pm785b07AX3kK0qS0VkhQ8FfiSmm4l4NInKtiDwD/AibvSuKklTUcw+F2Hvu1WCM+a4x5gzgncDnS7UTkZs8X37H8PBwGJdWFCVstFomFKSZJg7zLJxTRKSvxPFtxpgBY8xAf39/mJdWFCUs1HMPhdh77pUQkVPFmyFHRM4DWoGD9Z5XUZSIkBRgtGKmTqK2ZVoqNRCRu4BLgD4RGQRuATIAxpjbgd8F3i8ix4GjwO8FOlgVRUkakrZbM+e/VmrGlULORqTuFcXdGHN9heO3AreGFpGiKNHizGIzB6i4LxQ3QjWqXFdHqCqKkk9uCSH13eshaltGxV1RlHzyMndloaQSMEJVUZQTiZznrpl7PYgIIgkfxKQoShOhmXtopETUllEUJSakvMxdPfe6SYnaMoqixIVc5q4VzfUiIsyquCuKEgty4q6Ze72kRSJ7Rqq4K4qSj3ruoZHSDlVFUWKDeu6hkVJbRlGU2KCZe2ikUmrLKIoSF7TOPTS0WkZRlPigmXto2Dp3FXdFUeJAznNXca+XVEqYjeg2qrgripJPlJn73CzMTDf+uotESnRWSEVR4kKUde73/Q185aLGX3eRUFtGUZT4EGXmvvdhOPRC04yOTYnaMoqixIUo69yP7LYPleOTjb/2IpBKxdiWEZE7RGRIRHaWOH6DiDzu/fxKRM4NP0xFURpGlJn7kd12Oz3W+GsvAnG3Ze4Erixz/CXgzcaYc4DPA9tCiEtRlKiIqs59ahSOHravm0jcZ+M6iMkYsx04VOb4r4wx3ifCA8C6kGJTFCUKopoV0mXtANOjjb22Y/A38O0PwMyx0m0O74K/2QojgxVP10yDmD4M3BvyORVFaSRRraGaJ+41Zu4v3w97H60/hpd+Dk9+F3ZtL93mwFNw+CUYeqbi6VIi8fXcq0VELsWK+6fKtLlJRHaIyI7h4eGwLq0oSphU67nfeTU8/k8Lu8b22+D7H83fV4+43/tH8JM/XVgsQdx1n/5h5TbBbxc/+zx89+Z5TVMikY0FC0XcReQc4KvANcaYg6XaGWO2GWMGjDED/f39YVxaUZSwqcZzP34Udv0SXvrFwq6x6z546Zf5++oR96NH4NBLC4slyPS43T7749IjdJ2oB8X9iX+Cl381r+nlx/+VP9j3x/XHtQDqFncR2QDcA7zPGPNc/SEpihIp1WTuE6/a7cgrC7vG1Mh8X/3Iy9C5wr6uVdynRmD0lfpHt7rrjh+w3y7+9QvzRd61mfLiH3nFPpiKxPyamSc4e/LXfOaex3hk9+F5xxeTlkoNROQu4BKgT0QGgVuADIAx5nbgs8By4MsiAjBjjBlYrIAVRVlkqqlzn/S+oI/WI+5jttPW6oYV95VnwotDtYn73Jz3oDBWZPtOW1hMYK/bswHG9sHP/8LuO/tdsGJLfhvwH0677/f/Hfx9gB4zSoo5fvDQc2RaWnjthqULj61GKoq7Meb6Csc/AnwktIgURYmWajL3yUDmXiBoVTE1AnMz1t5p7bD7juyG9W+w9kYt4u6EHaw1U4+4HxuDnrVw5X+HA09agR95JV/cj3nWjcvcdz9gt3MzMDMFmSW5plljf48sE7ww7L3v+BRk2hceY5XoCFVFUfLJee7lxN2rjj4+AVNHaju/MVbcIWBxjNif3g3Q1r0Acfc4XKfvPj0GrV2w5WrY+l67b3RwfpvgdZ24B495ZI1t0yMTPD80bn/3L26Gn36uvjirQMVdUZR8avHcoXbf/fhRmDtuXzuBHN1nt9k1tYu7e1CArUGvh+kxe32A7tX2XhT+fkHP/egROLATlp/q7wuQnfPFfd/IFBOH9sL0CHStqi/OKlBxVxQln2rq3CcDRXG1+u5BMXZiOL7fbrtX1Sfu9VbMTI/74p5usSJc+PsFq2WGnwUMnPSm/GMAszN0GmvFZJkAYN8L3iwufafWF2cVqLgripJPtZ67s2+qGKmZR1CMnRiOHbDbrlXQll2YuGfXzbdlhp+zo06rJZi5g/XfC3+/oJU0MWRf9232jgXEfeoIKa8voEesuI8MPmWPLa+jX6BKVNwVRcmnGs994lVYfoptW2vmHhTAnLh7tkz3Si9zr2H6ASfua7ZaWyY4IvRHf2inEyjkie/At96Xv29u1vYhBMU9W0bcp0dhwhuMuexk/9irz8OLv/D7JYAeL3OfHX4OWtqhZ331v98CUXFXFCWfahbrmDwEXSutL12r517UljlgOzLbuhdgy3jnWH2urVYZ8yyemWnY8xCM7M4TWsAOwHr6B9b/d7gqmLzMfZ19eAUfGNOBahnX97D0JO/YmK2P/84H86yrHplgTU87bSMvwrJTfOtrEVFxVxQln1RB5r7zHttxGGTyVehYZjtAq83c7/tbuP/LBbaMJ+Jj++3DAhbuua/2Zhs/+Fu7feU3MOsNatr3mBVjl4W78wezcrevtcvf17POPjCCD4dg5j4+BO090NnnxTJq+w8mD8KhF3Nv6ZUJzl7bw7Kju+03ngag4q4oSj7pjN3OTMORPTYLfeLb+W0mD0JHn/WkqxH3/Tvhp7fAw1/LL53M2TL7bWcqzBd3Y2D3g9ZKAeujf+dD/syNUyNWkNd6YycHd9jty/f559j3mLVovv5O77re+fOmPCiSuWfX2q0rh3TWTcsSW9c+Mgid/f57psf8bP4V6/VPmQzL0kc5vb+d1XP7mVu2+J2poOKuKEohrkxvdC+M7LGvjwaGzs/N2ky2Y7nnSb9iM9iZ6eIZtzHwz5+23wRGBv1vAalMYLh/gbjPTlvxNga+8W6443K4+8O24/W5e2Hn3b4wT4142fNy21G55yG7f9d9sOIsO+J0z4PwzI98b99ZOcXms2nL+vt6PHF31pNrk11jt4desOKezljBD/rwnri/bFayNDXJ5taDtMgcI52bSt/7EFFxVxQln9YOm5WPDPq2RdBKOXoYMNaKWPUaK8RfPA3+fAXcdmqeHQHA4K+tx9232fraI3tsp2LHcnteY6xou4eKE9dj43ZKgud/Cutfb/cdeTmwWpMrSRzx37P+9VbIZ45Zkd90Eaw+B5691y7dd2wcZmcCtsweP053vraALZP1lqcYLRB3J/qHd/mWTFu3V0HjZe4HnmRa2jhgltIrk2xkLwD7Wxe/MxVU3BVFKUbPOk/cPfELWimuo7BjOZx7Hdz8b/DWz8KbP2Wz8/u/nH+uoaftdusNdnvgKZtpt2etoE6PWaujO+C5gz2259f29QU32e3hl+cvxecyd4D1F8DRQ/DQNnvOjRd5XnywQ3Q0YMsExL1Yh2pnv/2GUejVO9GfPWbbuPcd2e13RM8dZzzVwyid9DDOquP2WrtZTSOoOLeMoignID3r4OALvh0RzNxdZtqx3G5Xvcb+gM1wH/nfcOmf2A5XsNm2pP3se+ip/JGo44Ead8j3r/c8CJlO2HyFf65i4u7e667x01tseeLpV/lzvbR227ljpkZstg8lbJmAuKdS+Z3G7gHgMnfwxb09a22aABPpLCOmky4zQfv4iwyZXvYcbaMRaOauKMp8etbbrD2XuQfEPZi5F3Lhx2DmKOz4B3/f4V3Qux6WbrT/PjZuM+22rPW+xwKjU8GvVpkeg8GHYO15VnA7V5QWd5e5922G9l7b2XnZn0FLq83mV58Lr/ey/6kjJWyZItUyYMs9XYzOuskWEfe2bj/D96YunmjpYYROOufGyBx6lhdYx76Rqfn3bRFQcVcUZT4966wIH/BGVOaJu5e5O685yIotsPJs304Ba6Us3WRLHVOeWZCzZQKZe3eB5z6231bZrL/A/rt3A+x9xHrnUFzcUymbrZ/yVtjydrtvyVL4T9vhlLd45z1g7aNMh+1gnXXz3BSxZcDaRS7GnOceWCo657ln/fJRL+bxdA8jppM0s8iBnexr3ci+kUBt/SKi4q4oynyceI0W6VCd8DL3JcuKvze7Fsb2+v8+vAt6N9r6+W6vyqS9xx+J6ipYugo895e2W/96nSfuSzdasXe4+dOnRu2DwnHtV+D3754/DbF7aLhsfcUWK8a5ztJRW/HiSkEdXav86REKq2UgkLkHYlj3OgAm0zZzB2D2GEc6T9HMXVGUCCkcHh+c7fDQC9aqKDUnefeqgI0xbjP9pZu883p2RnsPtPX4tkxLu599O3F/6nt26wklvRvI7xgdg2MT9gHg3usoNr+8a5MT9zPtNmjztHXNf1/3SuvRHz/qi3v3asC7RtCWsReHtecDMNnSw6jpyJ1qsncz+46ouCuKEhVB22HpJr9kEWD4Geg/vfR7s2tsrffMMeuRu3MEz+tsmeMT8Opz9mHiBNll4dPjcMVf2Pp1sNm/I91mhdZ9oygU92LkxN37NrLybLt1FTPHxudbMuB31o7t98XdffOA/A5VsH0R/WeApDmSWeln7gArdRuiegAAC0hJREFUzmBobIqZ2cVfNbuiuIvIHSIyJCI7Sxw/Q0TuF5FpEflk+CEqitJwOvsh3WpfrzzbZsfHJuySdsPPQv+W0u/t9kr9xvf786u7ztRsMHMP2C8b3uC/v7UT3vOPcPMv4cKP+vt7N9htWw9kV9cu7m1ZQHxxX3GG3QbLHIuJuyvRHD9g22Q6rcXUlrVVQO293vkDYt/VDzf/G4/0Xs6I8cS9axXLlq9kzsDQWJ1rvVZBNZn7ncCVZY4fAv4r8MUwAlIUJQakUr4QO/tiasRaGscny2fuTtxH99nOVPAn1nKZe1vW96hnpmDTxfnnOPMd+UvbgZ/9B1dryg08ylKRVMq2c2Le0Wc7W920vdNjtlyykMLM3Yl4W7ftTHWTgOXE3etgXXkmJpXxM/f+01ndY62sRvjuFcXdGLMdK+Cljg8ZY34NHA8zMEVRIqZnnZ0h0gn51Ii3OAXWdihF1hP3sX02c2/LWhF15wTflnFsvKi6eBBP3LMFmXtvdb9Te9bvwHXlleMBcS+auXvi7jJ316Y961sy4D9gAvvSKfEz9xVbWN1rxX1/HMQ9TETkJhHZISI7hoeHG3lpRVFqZcUWWzfuBiNNjVi/HSpk7l4lydg+67n3bvT99FXn2Ex45Vm+SPZssHXwlWhpsxn+pov8SpvcgKql1f1O7T1+uWJ7FrqqEPcly2wJZy5z9zpdt94A593otysi7iLCKJ28etYH4dzrWJ21A6oaUQ7Z0BGqxphtwDaAgYEBU6G5oihRctnnbIWI6xR14t610hf8YnQss3796F5bJ7/2PP9Yz1r4pJf9H/Pq1TdVkbU7PvBDu937SPHRrZUIevOt3Vbc9z7ixTNevFomlbK/c2Hmfv6N+e3cN5EuX9xTAiCMX/oF+vo6yRpDR2u6IbaMTj+gKEpxWjvtj7M+nLiXy9rBZundq+zUuyO78ztFg3TZUZyc9ObaYwtOXdCWtZOdVYMT90yHt0bqSpu5G2N/v2KZO9h2I3vs1AlnXVs6JphnywB0tVupFRG2/9GlLO1orS7eOtBSSEVRyuMEceqIVylTxm93dK+B3ffb1xsvLN6mdz3cfB+c83u1x+TEfWy//5CohsJa+s5+m7EffslOAlZq+bvuVfDyr6wVdNrlxdv0bbYTnAWOpzw7qrvdz6P7utpyor+YVMzcReQu4BKgT0QGgVuADIAx5nYRWQXsALLAnIh8HDjTGFPDIoiKosQW5yXve9wKoaueKUd2NWDse109eTFWlTlWNiZvzveRPdVbMuD/Lm7rRsW6OeBdVU8hXSvtfDXpVjj5kuJt0hm46ra8XSLQmk7R1pKuPsaQqCjuxpjrKxzfD6wr10ZRlATT0mptjJe223+75ezK4coh17/eX7YvTJw4H3weTr2s+vcVZu4u69/zoN26cstCXMXMpouL+/IluPzMVXS3Zyo3XATUc1cUpTLtPdY/T7XYSpdKOHEvZcnUixPnqRE/+66GkuL+ELkyy2K4dqddUVOYF56ynAtPKTJ7ZgNQz11RlMq4THnFFluSWAknkhsvLt9uwfEEOj4XIu6ussWbmpcDT9o6+pYSHZ1rB2zJ5para481IjRzVxSlMk4Uq7FkAM7493DD3f50vWETFPfuGjz39gLPvbMPOwGYKW3JgF2q7xNP1BhktGjmrihKZXLivrW69ukMnHZZ8dkZwyAvc6+jWiad8Wv2y4l7AlFxVxSlMrWK+2ITnEumlmqZnLgH3+/ZOiruiqKccLT32HlmqulMbQQLtmUKMnfwBx01mbir564oSmUGPgRrtlY/EnSxceKcyviTklVD92roOz2/78Bl7stK1LgnFBV3RVEqs+rshQ84WgwyHfabRNfK2nz9zBL42EP5+5xnX2oAU0JRcVcUJXmI2Oy9ls7UUpz1LjvQqpZvAAlAxV1RlGTSlq3Nby/FuvPtT5Oh4q4oSjK59E/y13pV8lBxVxQlmWx9b9QRxBothVQURWlCVNwVRVGaEBV3RVGUJkTFXVEUpQlRcVcURWlCKoq7iNwhIkMisrPEcRGRvxWR50XkcRE5r1g7RVEUpXFUk7nfCVxZ5vjbgNO8n5uAr9QflqIoilIPFcXdGLMdOFSmyTXA143lAaBXRFaHFaCiKIpSO2EMYloL7An8e9Dbt6+woYjchM3uAcZF5NkFXrMPeHWB711s4hqbxlUbcY0L4hubxlUbC41rYzWNwhD3YlOymWINjTHbgG11X1BkhzFmoN7zLAZxjU3jqo24xgXxjU3jqo3FjiuMaplBYH3g3+uAvSGcV1EURVkgYYj7D4D3e1UzbwBGjDHzLBlFURSlcVS0ZUTkLuASoE9EBoFbgAyAMeZ24MfAVcDzwCTwwcUKNkDd1s4iEtfYNK7aiGtcEN/YNK7aWNS4xJii9riiKIqSYHSEqqIoShOSOHEXkStF5FlvROynI4xjvYj8PxF5WkSeFJH/5u3/nIi8IiKPej9XRRDbLhF5wrv+Dm/fMhH5vyLyW2/b8DXFROT0wH15VERGReTjUdyzYiOvS92jRo7CLhHXbSLyjHft74pIr7d/k4gcDdy32xscV8nPTUT+2Ltfz4rIFYsVV5nYvhWIa5eIPOrtb+Q9K6URjfk7M8Yk5gdIAy8AJwOtwGPAmRHFsho4z3vdDTwHnAl8DvhkxPdpF9BXsO9/Ap/2Xn8auDUGn+V+bM1uw+8Z8CbgPGBnpXuE7VO6F1v2+wbgwQbHdTnQ4r2+NRDXpmC7CO5X0c/N+3/wGNAGnOT9n003MraC438JfDaCe1ZKIxryd5a0zP0C4HljzIvGmGPAN7EjZBuOMWafMeZh7/UY8DR28FZcuQb4mvf6a8A7I4wF4K3AC8aYl6O4uCk+8rrUPWrYKOxicRljfmKMmfH++QC23LihlLhfpbgG+KYxZtoY8xK22OKCKGITEQHeA9y1WNcvRRmNaMjfWdLEvdRo2EgRkU3Aa4EHvV0f875W3RGF/YEdRPYTEfmN2FHBACuNV6LqbUNYNr4uriP/P1zU9wxK36M4/d19CJvdOU4SkUdE5Bci8sYI4in2ucXpfr0ROGCM+W1gX8PvWYFGNOTvLGniXvVo2EYhIl3A3cDHjTGj2InTTgG2Yqdg+MsIwrrIGHMedlK3j4rImyKIoSQi0gq8A/i2tysO96wcsfi7E5HPADPAN7xd+4ANxpjXAn8I/B8RyTYwpFKfWyzul8f15CcRDb9nRTSiZNMi+xZ835Im7rEaDSsiGeyH9g1jzD0AxpgDxphZY8wc8Pcs4tfRUhhj9nrbIeC7XgwH3Fc8bzvU6LgCvA142BhzAOJxzzxK3aPI/+5E5EbgauAG4xm0nu1x0Hv9G6y3vblRMZX53CK/XwAi0gK8C/iW29foe1ZMI2jQ31nSxP3XwGkicpKX/V2HHSHbcDwv7x+Ap40xfxXYH/TIrgWKzoO/iHF1iki3e43tjNuJvU83es1uBL7fyLgKyMumor5nAUrdo0hHYYvIlcCngHcYYyYD+/tFJO29Phk77faLDYyr1Of2A+A6EWkTkZO8uB5qVFwBLgOeMcYMuh2NvGelNIJG/Z01otc4zB9sj/Jz2CfuZyKM42LsV6bHgUe9n6uAfwSe8Pb/AFjd4LhOxlYqPAY86e4RsBz4GfBbb7ssovvWARwEegL7Gn7PsA+XfcBxbMb04VL3CPt1+e+8v7kngIEGx/U81ot1f2e3e21/1/uMHwMeBt7e4LhKfm7AZ7z79SzwtkZ/lt7+O4GbC9o28p6V0oiG/J3pCFVFUZQmJGm2jKIoilIFKu6KoihNiIq7oihKE6LiriiK0oSouCuKojQhKu6KoihNiIq7oihKE6LiriiK0oT8f0lc2SkhPU5uAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# here the atom pairs are different than the dimers\n",
"dists = md.compute_distances(traj, atom_pairs=[[0,1], [1,2]]) / sigma_nm\n",
"plt.plot(dists)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.1"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment