{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 03. Coarse Grained Simulations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Coarse-graining a box of 500 methanol molecules from the all atom simulation given in `CG_tutorial/meth.gsd` file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "# reduce logging\n",
    "tf.get_logger().setLevel('ERROR')\n",
    "from tensorflow.keras import layers\n",
    "import hoomd\n",
    "import hoomd.md\n",
    "import hoomd.htf as htf\n",
    "import numpy as np\n",
    "import gsd, gsd.hoomd, pickle\n",
    "import matplotlib.pyplot as plt,matplotlib\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build the computational graph\n",
    "\n",
    "Here we prepare the computations that will be executed at each step during the simulation. We have access to placeholders for the neighbor list, positions, types, box dimensions of the simulation.\n",
    "\n",
    "This graph works in three steps:\n",
    "\n",
    "1. Creaet a mapping matrix that describes how to group atoms together in a coarse-grained system given the all-atom system\n",
    "2. Create a coarse-grained mapped trajectory, where atoms/forces are grouped together using the mapping matrix\n",
    "3. Compute the radial distribution function of this new mapped trajectory\n",
    "4. Compute the radial distribution function of the C-C from the all-atom sysetm for comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Graph building code ###\n",
    "\n",
    "def make_train_graph(NN, rcut, NN_aa, N_hidden, system):\n",
    "\n",
    "    #mamximum number of neighbors=0, graph will not compute forces to use in TF\n",
    "    graph = htf.graph_builder(NN_aa, output_forces=False)\n",
    "    assign_ops = []\n",
    "    prints = []\n",
    "    #get mapping from molecule index to particle index\n",
    "    molecule_mapping = htf.find_molecules(system)\n",
    "\n",
    "    for m in molecule_mapping:\n",
    "        #ensure all molecules are of the same length\n",
    "        assert len(m) == len(molecule_mapping[0])\n",
    "    N = sum([len(m) for m in molecule_mapping])\n",
    "    M = len(molecule_mapping)\n",
    "    #number of atoms in a molecule=MN\n",
    "    MN = len(molecule_mapping[0])\n",
    "    # create one bead mapping. mapping matrices must be created for other types of mappings\n",
    "    molecule_mapping_matrix = np.ones([1, MN], dtype=np.int)\n",
    "    #create a mass-weighted M x N mapping operator \n",
    "    cg_mapping = htf.sparse_mapping([molecule_mapping_matrix for _ in molecule_mapping], \n",
    "                                    molecule_mapping, system=system)\n",
    "    assert cg_mapping.shape == (M, N)\n",
    "    #calculate the center of mass of a CG bead\n",
    "    mapped_pos = htf.center_of_mass(graph.positions[:,:3], cg_mapping, graph.box_size, name='com-mapped-positions')\n",
    "    #create the mapped neighbot list\n",
    "    mapped_nlist = htf.compute_nlist(mapped_pos, rcut, NN, graph.box_size)\n",
    "    #create a non mass-weighted M x N mapping operator\n",
    "    force_cg_mapping = htf.sparse_mapping([molecule_mapping_matrix for _ in molecule_mapping], molecule_mapping)\n",
    "    #calculate the coarse grained forces\n",
    "    mapped_force = tf.sparse.matmul(force_cg_mapping, graph.forces, name='mapped-forces')\n",
    "    #compute RDF for mapped and C-C in all-atom\n",
    "    rdf = graph.compute_rdf([0.1,rcut], 'rdf', nbins=200, nlist=mapped_nlist)\n",
    "    aa_rdf = graph.compute_rdf([0.1,rcut], 'aa-rdf', nbins=200, type_i=3, type_j=3)\n",
    "    graph.running_mean(rdf,'avg-rdf') \n",
    "    graph.running_mean(aa_rdf,'avg-aa-rdf') \n",
    "    #save the graph\n",
    "    graph.save('CG_tutorial/cg_model')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running the simulation\n",
    "\n",
    "Here we begin a simulation of methanol. This code is a little complex, but not really about hoomd-tf. This is the details of setting-up the force field to simulate methanol (e.g., treating electrostatics, dispersion, thermostat). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HOOMD-blue v2.5.1 DOUBLE HPMC_MIXED SSE SSE2 \n",
      "Compiled: 03/04/2020\n",
      "Copyright (c) 2009-2019 The Regents of the University of Michigan.\n",
      "-----\n",
      "You are using HOOMD-blue. Please cite the following:\n",
      "* J A Anderson, C D Lorenz, and A Travesset. \"General purpose molecular dynamics\n",
      "  simulations fully implemented on graphics processing units\", Journal of\n",
      "  Computational Physics 227 (2008) 5342--5359\n",
      "* J Glaser, T D Nguyen, J A Anderson, P Liu, F Spiga, J A Millan, D C Morse, and\n",
      "  S C Glotzer. \"Strong scaling of general-purpose molecular dynamics simulations\n",
      "  on GPUs\", Computer Physics Communications 192 (2015) 97--107\n",
      "-----\n",
      "HOOMD-blue is running on the CPU\n",
      "notice(2): Group \"all\" created containing 6000 particles\n",
      "Finding molecules...99.90%\n",
      "Note: Backed-up CG_tutorial/cg_model previous model to CG_tutorial/cg_model/previous_model_7\n",
      "notice(2): Started TF Session Manager.\n",
      "-----\n",
      "You are using PPPM. Please cite the following:\n",
      "* D N LeBard, B G Levine, S A Barr, A Jusufi, S Sanders, M L Klein, and A Z\n",
      "  Panagiotopoulos. \"Self-assembly of coarse-grained ionic surfactants\n",
      "  accelerated by graphics processing units\", Journal of Computational Physics 8\n",
      "  (2012) 2385-2397\n",
      "-----\n",
      "notice(2): -- Neighborlist exclusion statistics -- :\n",
      "notice(2): Particles with 5 exclusions             : 6000\n",
      "notice(2): Neighbors included by diameter          : no\n",
      "notice(2): Neighbors excluded when in the same body: no\n",
      "** starting run **\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "*Warning*: charge.pppm: system is not neutral and unscreened interactions are calculated, the net charge is -0.000357628\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "notice(2): charge.pppm: RMS error: 0.000406376\n",
      "Time 00:00:13 | Step 28 / 200 | TPS 2.77151 | ETA 00:01:02\n",
      "Time 00:00:24 | Step 57 / 200 | TPS 2.80709 | ETA 00:00:50\n",
      "Time 00:00:34 | Step 85 / 200 | TPS 2.78746 | ETA 00:00:41\n",
      "Time 00:00:44 | Step 115 / 200 | TPS 2.97403 | ETA 00:00:28\n",
      "Time 00:00:54 | Step 144 / 200 | TPS 2.87635 | ETA 00:00:19\n",
      "Time 00:01:04 | Step 171 / 200 | TPS 2.69084 | ETA 00:00:10\n",
      "Time 00:01:14 | Step 199 / 200 | TPS 2.68526 | ETA 00:00:00\n",
      "Time 00:01:15 | Step 200 / 200 | TPS 3.30671 | ETA 00:00:00\n",
      "Average TPS: 2.80062\n",
      "---------\n",
      "-- Neighborlist stats:\n",
      "52 normal updates / 0 forced updates / 0 dangerous updates\n",
      "n_neigh_min: 0 / n_neigh_max: 533 / n_neigh_avg: 211.555\n",
      "shortest rebuild period: 3\n",
      "-- Cell list stats:\n",
      "Dimension: 3, 3, 3\n",
      "n_min    : 168 / n_max: 294 / n_avg: 222.222\n",
      "** run complete **\n",
      "notice(2): Force mode is FORCE_MODE.hoomd2tf \n",
      "notice(2): Starting TensorflowCompute \n",
      "notice(2): completed reallocate\n",
      "notice(2): TF Session Manager has released control. Starting HOOMD updates\n",
      "** starting run **\n",
      "Time 00:01:26 | Step 214 / 700 | TPS 1.39365 | ETA 00:05:48\n",
      "Time 00:01:36 | Step 230 / 700 | TPS 1.53088 | ETA 00:05:07\n",
      "Time 00:01:47 | Step 247 / 700 | TPS 1.63033 | ETA 00:04:37\n",
      "Time 00:01:57 | Step 264 / 700 | TPS 1.62376 | ETA 00:04:28\n",
      "Time 00:02:08 | Step 281 / 700 | TPS 1.58996 | ETA 00:04:23\n",
      "Time 00:02:19 | Step 298 / 700 | TPS 1.61569 | ETA 00:04:08\n",
      "Time 00:02:29 | Step 315 / 700 | TPS 1.6277 | ETA 00:03:56\n",
      "Time 00:02:39 | Step 332 / 700 | TPS 1.63728 | ETA 00:03:44\n",
      "Time 00:02:50 | Step 349 / 700 | TPS 1.62727 | ETA 00:03:35\n",
      "Time 00:03:00 | Step 366 / 700 | TPS 1.68034 | ETA 00:03:18\n",
      "Time 00:03:10 | Step 383 / 700 | TPS 1.6327 | ETA 00:03:14\n",
      "Time 00:03:20 | Step 400 / 700 | TPS 1.69372 | ETA 00:02:57\n",
      "Time 00:03:31 | Step 417 / 700 | TPS 1.63694 | ETA 00:02:52\n",
      "Time 00:03:41 | Step 434 / 700 | TPS 1.64186 | ETA 00:02:42\n",
      "Time 00:03:51 | Step 451 / 700 | TPS 1.64749 | ETA 00:02:31\n",
      "Time 00:04:02 | Step 468 / 700 | TPS 1.64073 | ETA 00:02:21\n",
      "Time 00:04:12 | Step 485 / 700 | TPS 1.64216 | ETA 00:02:10\n",
      "Time 00:04:22 | Step 502 / 700 | TPS 1.68894 | ETA 00:01:57\n",
      "Time 00:04:32 | Step 519 / 700 | TPS 1.68217 | ETA 00:01:47\n",
      "Time 00:04:43 | Step 536 / 700 | TPS 1.65193 | ETA 00:01:39\n",
      "Time 00:04:53 | Step 553 / 700 | TPS 1.67843 | ETA 00:01:27\n",
      "Time 00:05:03 | Step 570 / 700 | TPS 1.64553 | ETA 00:01:19\n",
      "Time 00:05:13 | Step 587 / 700 | TPS 1.64451 | ETA 00:01:08\n",
      "Time 00:05:23 | Step 604 / 700 | TPS 1.68972 | ETA 00:00:56\n",
      "Time 00:05:34 | Step 621 / 700 | TPS 1.67757 | ETA 00:00:47\n",
      "Time 00:05:44 | Step 638 / 700 | TPS 1.69154 | ETA 00:00:36\n",
      "Time 00:05:54 | Step 655 / 700 | TPS 1.68949 | ETA 00:00:26\n",
      "Time 00:06:04 | Step 673 / 700 | TPS 1.71408 | ETA 00:00:15\n",
      "Time 00:06:15 | Step 690 / 700 | TPS 1.64928 | ETA 00:00:06\n",
      "Time 00:06:21 | Step 700 / 700 | TPS 1.66766 | ETA 00:00:00\n",
      "Average TPS: 1.64154\n",
      "---------\n",
      "-- Neighborlist stats:\n",
      "108 normal updates / 1 forced updates / 0 dangerous updates\n",
      "n_neigh_min: 307 / n_neigh_max: 483 / n_neigh_avg: 409.472\n",
      "shortest rebuild period: 4\n",
      "-- Cell list stats:\n",
      "Dimension: 3, 3, 3\n",
      "n_min    : 189 / n_max: 259 / n_avg: 222.222\n",
      "** run complete **\n",
      "notice(2): Sending exit signal.\n",
      "notice(2): Shutting down TF Manually.\n",
      "notice(2): TF Queue is waiting, sending None\n"
     ]
    }
   ],
   "source": [
    "### Hoomd-Sim code ###\n",
    "\n",
    "g = gsd.hoomd.open('CG_tutorial/meth.gsd')\n",
    "c = hoomd.context.initialize('--mode=cpu')\n",
    "system = hoomd.init.read_gsd(filename='CG_tutorial/meth.gsd')\n",
    "c.sorter.disable()\n",
    "set_rcut=10.0\n",
    "\n",
    "#call the graph building code\n",
    "make_train_graph(128,10.0,256,48, system)\n",
    "\n",
    "with htf.tfcompute('CG_tutorial/cg_model', device='/CPU:0') as tfcompute:\n",
    "    nlist = hoomd.md.nlist.cell()\n",
    "\n",
    "    #set-up pppm\n",
    "    charged = hoomd.group.all()\n",
    "    pppm = hoomd.md.charge.pppm(nlist=nlist, group=charged)\n",
    "    pppm.set_params(Nx=32, Ny=32, Nz=32, order=6, rcut=set_rcut)\n",
    "\n",
    "    #set-up pair coefficients\n",
    "    nlist.reset_exclusions(['1-2', '1-3', '1-4','body'])\n",
    "    lj = hoomd.md.pair.force_shifted_lj(r_cut=set_rcut, nlist=nlist)\n",
    "    forces = [lj]\n",
    "    lj.pair_coeff.set(\"opls_156\", \"opls_156\", sigma=2.5, epsilon=0.03)\n",
    "    lj.pair_coeff.set(\"opls_156\", \"opls_157\", sigma=2.96, epsilon=0.05)\n",
    "    lj.pair_coeff.set(\"opls_156\", \"opls_154\", sigma=2.79, epsilon=0.07)\n",
    "    lj.pair_coeff.set(\"opls_156\", \"opls_155\", sigma=5.0, epsilon=0.0)\n",
    "    lj.pair_coeff.set(\"opls_157\", \"opls_157\", sigma=3.5, epsilon=0.07)\n",
    "    lj.pair_coeff.set(\"opls_157\", \"opls_154\", sigma=3.31, epsilon=0.11)\n",
    "    lj.pair_coeff.set(\"opls_157\", \"opls_155\", sigma=5.92, epsilon=0.0)\n",
    "    lj.pair_coeff.set(\"opls_154\", \"opls_154\", sigma=3.12, epsilon=0.17)\n",
    "    lj.pair_coeff.set(\"opls_154\", \"opls_155\", sigma=5.59, epsilon=0.0)\n",
    "    lj.pair_coeff.set(\"opls_155\", \"opls_155\", sigma=10.0, epsilon=0.0)\n",
    "\n",
    "    #set-up bonds\n",
    "    harmonic = hoomd.md.bond.harmonic()\n",
    "    harmonic.bond_coeff.set(\"opls_156-opls_157\", k=340.00, r0=1.09)\n",
    "    harmonic.bond_coeff.set(\"opls_154-opls_157\", k=320.00, r0=1.41)\n",
    "    harmonic.bond_coeff.set(\"opls_154-opls_155\", k=553.00, r0=0.95)\n",
    "\n",
    "    #set-up angles\n",
    "    harm_angle = hoomd.md.angle.harmonic()\n",
    "    harm_angle.angle_coeff.set(\"opls_154-opls_157-opls_156\", k=70.0, t0=1.90)\n",
    "    harm_angle.angle_coeff.set(\"opls_155-opls_154-opls_157\", k=110.0, t0=1.89)\n",
    "    harm_angle.angle_coeff.set(\"opls_156-opls_157-opls_156\", k=66.0, t0=1.88)\n",
    "\n",
    "    #set-up dihedrals\n",
    "    dihedral = hoomd.md.dihedral.opls()\n",
    "    dihedral.dihedral_coeff.set(\"opls_155-opls_154-opls_157-opls_156\", k1=0.0, k2=0.0, k3=0.45, k4=0.0)\n",
    "\n",
    "    group_all = hoomd.group.all()\n",
    "    kT = 1.9872/1000\n",
    "\n",
    "    #NVT Simulation in Hoomd\n",
    "    im = hoomd.md.integrate.mode_standard(dt=5.0/489.0)\n",
    "    nvt = hoomd.md.integrate.nvt(group=group_all, kT=298.15 * kT, tau=350 / 48.9)\n",
    "    nvt.randomize_velocities(1234)\n",
    "\n",
    "    #equilibrate\n",
    "    hoomd.run(200)\n",
    "\n",
    "    #communicate positions, neighbor list and forces to TensorFlow model\n",
    "    tfcompute.attach(nlist, r_cut=set_rcut)\n",
    "\n",
    "    #Hoomd production run\n",
    "    hoomd.run(500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "Now we load the RDF we computed and plot it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Load saved rdf tensor from the graph\n",
    "train_vars = htf.load_variables('CG_tutorial/cg_model', ['avg-rdf', 'rdf-r', 'avg-aa-rdf'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEQCAYAAABIqvhxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nOzdd3iUVfrw8e+ZyaT33gMptNAJva8iRQQUWMGGoouyyq7uqrvqqmvdtfwsu+rrWrEgWEGUjgjSIaGGEgjpvfcy7bx/TIiJCZBAyiScz3XlgjnPeZ65Q8k9pwspJYqiKIrSGprODkBRFEXpelTyUBRFUVpNJQ9FURSl1VTyUBRFUVpNJQ9FURSl1VTyUBRFUVqtw5KHECJECPGzEOKkEOKEEOLPzdQRQoj/CCEShRDHhBBDG1xbJIQ4W/e1qKPiVhRFUZoSHbXOQwgRAARIKQ8JIVyAOGCOlPJkgzozgGXADGAk8KaUcqQQwhOIBWIAWXfvMCllcYcEryiKojTSYS0PKWW2lPJQ3e/LgVNA0G+qzQY+lRb7APe6pDMV2CKlLKpLGFuAaR0Vu6IoitJYp4x5CCF6AEOA/b+5FASkN3idUVd2oXJFURSlE9h09BsKIZyBb4EHpZRl7fD8JcASACcnp2F9+vRp67dQFEXptuLi4gqklD6XqtehyUMIocOSOFZIKb9rpkomENLgdXBdWSYw6Tfl25t7Dynle8B7ADExMTI2NvaK4+5uissNSMDTRdfZoSiKYmWEEKktqdeRs60E8CFwSkr52gWqrQXuqJt1NQoolVJmA5uA64QQHkIID+C6ujKllRKzqrjvzdP88c0EUnOrOzscRVG6qI4c8xgL3A78TghxpO5rhhDiPiHEfXV11gNJQCLwPvBHACllEfAccLDu69m6MuUScov1nMuqQm80czqtksfeP4edjQatBh774BzpeTWdHaKiKF1Qh03V7QxXe7eVlJK7XjlFbrEejQY0QuDjpuNf90SiN5j52/uJCAFv3N8LHzfbzg5XURQrIISIk1LGXKpehw+YKx0nJaeG3GI914/0wtnBhvJqI7f8zh8vV8tYx4t3R7DsrTN8vSOPP84K7uRole7OYDCQkZFBTY1q7VoDe3t7goOD0ekub+xTJY9u7GCCZTLbwgYJo6Ee/g5MHuzB5thCbr3GHzcn9c9BaT8ZGRm4uLjQo0cPLEOgSmeRUlJYWEhGRgY9e/a8rGeova26sYMJZUQEOjSbOM6bO96XWoPkx70FHRiZcjWqqanBy8tLJQ4rIITAy8vrilqBKnl0U+XVRk6mVTK8t+tF64X52TOyrytr9+VTozd3UHTK1UolDutxpX8XKnl0U4fOlmM2w4hLJA+AeeN9Kas0sTVOTWBTur81a9YghOD06dP1ZSkpKfTv3x+A7du3M3PmzAve/+CDDxIUFITZ/OuHrTVr1nDy5MkL3nMlKioquPfee4mIiGDYsGFMmjSJ/ft/uzmHxYYNG4iJiaFfv34MGTKEv/71r+0SE6jk0W0dPF2Gi4OWXiGOl6wb3cOJPiGOrNiWo6buKt3eypUrGTduHCtXrmz1vWazmdWrVxMSEsKOHTvqy9szedxzzz14enpy9uxZ4uLi+PjjjykoaNrNHB8fzwMPPMDnn3/OyZMniY2NJTIysl1iApU8uiWzWRJ7ppxhvVzRai7dNBVC8NC8UATwtw8SSVMJROmmKioq2LVrFx9++CGrVq1q9f3bt28nOjqapUuX1iefPXv2sHbtWh555BEGDx7MuXPnOHLkCKNGjWLgwIHceOONFBdbNgCfNGkSDz30EDExMfTt25eDBw9y0003ERUVxT/+8Y8m73fu3Dn279/P888/j0Zj+XHds2dPrr/++iZ1X375ZZ544gnOb8mk1WpZunRpq7/HllLTa7qhs5lVlFYaGdHn0l1W54X62vPvP0Ty9/cT+fv7ibxybxRB3nbtGKVyNXv3hwySstt2h4PwAAfuu+HiU86///57pk2bRq9evfDy8iIuLo5hw4a1+D1WrlzJwoULmT17No8//jgGg4ExY8Ywa9YsZs6cybx58wAYOHAg//3vf5k4cSJPPfUUzzzzDG+88QYAtra2xMbG8uabbzJ79mzi4uLw9PQkIiKChx56CC8vr/r3O3HiBIMHD0ar1V4ytvj4+Hbtpvot1fLohnYcLUGrgWFRLq26L9TXnpf+EInBJHnj2zTM5u67gFS5Oq1cuZIFCxYAsGDBglZ1Xen1etavX8+cOXNwdXVl5MiRbNrUdJek0tJSSkpKmDhxIgCLFi3il19+qb8+a9YsAAYMGEB0dDQBAQHY2dkRHh5Oenp6k+dZK9Xy6Gb0BjNbDxUxJtod18tYtxHia8890wN547t0thwqYmqM16VvUpRWulQLoT0UFRWxbds2jh8/jhACk8mEEIJXXnmlRfdv2rSJkpISBgwYAEBVVRUODg4XHVxvjp2dpUWv0Wjqf3/+tdFobFQ3Ojqao0ePYjKZmrQ+3n77bd5//30A1q9fT3R0NHFxcQwaNKhV8Vwu1fLoZnbFl1BebWL6iMv/oT9lmCfRPZz4cH0WJRXGS9+gKF3AN998w+23305qaiopKSmkp6fTs2dPdu7c2aL7V65cyQcffEBKSgopKSkkJyezZcsWqqqqcHFxoby8HAA3Nzc8PDzqn/vZZ5/Vt0JaKyIigpiYGJ5++mnObyWVkpLCunXruP/++zly5AhHjhwhMDCQRx55hBdffJEzZ84AlsH9d99997LetyVU8uhm1h8oJNDLlkHhzpf9DI1GsGxOCNV6M//7MUN1XyndwsqVK7nxxhsblc2dO7dFXVdVVVVs3Lix0UC1k5MT48aN44cffmDBggW88sorDBkyhHPnzvHJJ5/wyCOPMHDgQI4cOcJTTz112XF/8MEH5ObmEhkZSf/+/bnzzjvx9fVtUm/gwIG88cYbLFy4kL59+9K/f3+SkpIu+30vRW2M2I2k5lZz3xsJ3D09kHkTmv7jaq0VP+Xw+dYcxvV346/zQ7G3vfSgnaJcyKlTp+jbt29nh6E00NzfidoY8Sq04UAhNlrBtUM92+R5t/zODwdbDR9uyCKr8Cz/XBSudt9VFAVQ3Vbdhsks2Xa4mDHRbrg7t81nAiEEN4335ZlF4WQX6XlxRQomU/dtqSqK0nIqeXQTyTnVlFebGNmKtR0tFdPblT/dGMLp9CpWbc9t8+critL1qOTRTRxPqgBgwBUMlF/MpEEeTBrkwRfbckhIr2yX91AUpetQyaObOJ5cQYCnbbuOSdw/OwgvFx2vfJmmduBVlKtchyUPIcRHQog8IUT8Ba4/0uBs83ghhEkI4Vl3LUUIcbzu2tUzfaqFzGZJfHJlu7U6znN2sOEv80PJLKzl083Z7fpeiqJYt45seSwHpl3oopTyFSnlYCnlYOAxYIeUsuEe4ZPrrl9yCtnVJjW3hvJqEwN6tm/yABgc4cL1I71Ysyefk6mq+0rpWoQQ3HbbbfWvjUYjPj4+rV4l3lYabgX/W9nZ2Y3iOnDgABMmTKB3794MGTKEe+65h6qqKgA2btzIiBEj6NOnD4MHD+bmm28mLS0NgIcffpht27a1eewdljyklL8ALT0wYiHQ+v2Sr1LHkuvGOzogeQAsnh6Ij5uO179No9aguq+UrsPJyYn4+Hiqqy2bMm7ZsoWgoKBOjqp5r732Gn/4wx8AyM3NZf78+bz00kskJCRw+PBhpk2bRnl5OfHx8SxbtoxPPvmE06dPc+TIEW699VZSUlIAWLZsGf/+97/bPD6rG/MQQjhiaaF826BYApuFEHFCiCWdE5n1Op5cga+7Dj+PjlmD4Win5cG5oWTk1/Llz2r2ldK1zJgxg3Xr1gG/7pJ73oEDBxg9ejRDhgxhzJgxJCQkALB8+XJmz57NpEmTiIqK4plnngEsLYc+ffpw66230rdvX+bNm1ffGoiLi2PixIkMGzaMqVOnkp2dXV8+aNAgBg0axNtvv33BOL/99lumTbN01rz99tssWrSI0aNH11+fN28efn5+vPTSSzz++OONFvvNmjWLCRMmABAWFkZhYSE5OTlX/GfXkDUuErwB2P2bLqtxUspMIYQvsEUIcbquJdNEXXJZAhAaGtr+0XYyKS3jHTG9WreD7pUaEunChIHurNmdz43jfHBxtMZ/Soq1+ib9SzKq2nYH2WDHEOaF3HzJegsWLODZZ59l5syZHDt2jMWLF9fvQ9WnTx927tyJjY0NW7du5fHHH+fbby2fYw8cOEB8fDyOjo4MHz6c66+/Hm9vbxISEvjwww8ZO3Ysixcv5p133uHPf/4zy5Yt4/vvv8fHx4cvv/ySJ554go8++oi77rqLt956iwkTJvDII480G2NycjIeHh71GyfGx8ezaNGiZuueOHGChx9++KLf89ChQ9m9ezdz58695J9PS1ldywNYwG+6rKSUmXW/5gGrgREXullK+Z6UMkZKGePj49OugVqDtLxaSiuN7T5Y3pwFk/2o1ptZu6fpqWaKYq0GDhxISkoKK1euZMaMGY2ulZaWMn/+fPr3789DDz3EiRMn6q9NmTIFLy8vHBwcuOmmm9i1axcAISEhjB07FoDbbruNXbt2kZCQQHx8PFOmTGHw4ME8//zzZGRkUFJSQklJSX2r4Pbbb282xuzsbC7n51dhYSGDBw+mV69evPrqq/Xlvr6+ZGVltfp5F2NVHxeFEG7AROC2BmVOgEZKWV73++uAZzspRKtzvIPHOxrq6e/A6H5urNmTz5xxPjjZq72vlJZpSQuhPc2aNYuHH36Y7du3U1hYWF/+5JNPMnnyZFavXk1KSgqTJk2qvyZE41M5z79urlxKSXR0NHv37m10raSkpEXxOTg4UFPz64me57dbnz17dpO60dHRHDp0iEGDBuHl5cWRI0d49dVXqaioqK9TU1ODg4NDi967pTpyqu5KYC/QWwiRIYS4WwhxnxDivgbVbgQ2SykbTuPxA3YJIY4CB4B1UsqNHRW3tTuTUYW7sw0Bnp2z59SCyX5UVJtYt7+AwjIDP+zNJz654tI3KkonWrx4MU8//XT92RznlZaW1g+gL1++vNG1LVu2UFRURHV1NWvWrKlvbaSlpdUniS+++IJx48bRu3dv8vPz68sNBgMnTpzA3d0dd3f3+lbLihUrmo2vV69e9QPeAA888ACffPIJ+/fvry/77rvvyM3N5dFHH+WFF17g1KlT9dfOj7ucd+bMmQvO6rpcHdbykFIubEGd5Vim9DYsSwI65nSTLuhcVhURAQ5NPv10lF7BjsT0cuHzrTks35SNlODhbMNHj/TD3tYae0UVBYKDg/nTn/7UpPzRRx9l0aJFPP/8803OCR8xYgRz584lIyOD2267jZiYGFJSUujduzdvv/02ixcvpl+/fixduhRbW1u++eYb/vSnP1FaWorRaOTBBx8kOjqajz/+mMWLFyOE4Lrrrms2PicnJyIiIkhMTCQyMhI/Pz9WrVrFww8/TF5eHhqNhgkTJjBt2jT8/Px48803ueOOOygrK8Pb25vQ0ND6QX2DwUBiYiIxMW27ykFtyd6F6Y1mbnr6GHPH+3LXtMBOi+NcVhXv/pjJ4HAXArxseeWrNBZPC2D+RL9Oi0mxPl15S/bly5cTGxvLW2+91ag8JSWFmTNnEh/f7NrnK7J69Wri4uJ4/vnnr/g5hw4d4rnnnmtyTW3JfpVKy63BZIbIIMdOjSMi0JFXlkTVv/75SDFf78hjxkhvNQ6iKJfpxhtvbDQec7mMRiN//etf2yCixlS/Qhd2Lsuy0CkioG0Hwq7UHdcFUF5tYs3u/M4ORVHaxJ133tmk1QHQo0ePdml1nHfPPfdc8TPmz5+Pu7t7G0TTmEoeXVhiVjUOdhr8O2mw/EKighwZE+3GdzvzKCjVd3Y4iqK0A5U8urDzg+UaTdsOlpfoizGajZesd7zkKDk1zW+QeMeUAKSEh/7fWZKzq9s0PqXr6s5jrF3Nlf5dqOTRRZnMkqTsGiIC27bLqlhfzFPxj/P3Yw+zIuVTkiuSmq2XUZXO/869w7uJb2MwG5pcD/Oz59X7opAS/vruWeLOlLVpnErXY29vT2FhoUogVkBKSWFhIfb29pf9DDVg3kVlFdRSazATEdi2g+Wnyk5gkiainKOIKz7I3sLd3NHjLkZ4jaqvI6Xk24yv0Gl05NfmsS13K1MDpjd5VniAA2/8MYonlyfx8pepfPK3ftjbqgH0q1VwcDAZGRnk56uxMGtgb29PcHDwZd+vkkcXlVg3WB7Zxi2PU2UncdO5syTij9Saa/nfuXf4NOVjTNLEaG/LoqhjpUc5U57A70MWklB+mo056xjhNQoPW48mz/N2s2XZnBD++u5ZNhwo5MZxvm0ar9J16HQ6evbs2dlhKG1EdVt1UeeyqtDZCEJ8L7/Z+VtmaSah7BR9XfshhMBea8/SyPvp49qXz1M/YXnyhxwqjmV1xjf42wcwzmcCc4PnI6VkdcbXF3xuvzAnBoY78+3OfPRqC3dF6RZU8uiizmVV08PPHhtt2w2Wp1elUWmqpK9rv/oyW40d90bczwSfSZwsjefDpPfIr81jbvDv0QotXnbeTPGfSlxx7EV3SV0w2Y/CMgNbDrX0SBdFUayZSh5dkJSSc1nV7TDecRKA3i59GpXrNDpuDr2Ffw16lQd7Pczd4Uvo5xZdf32S7zXYCBv2Fu6+4LMHRzjTO8SRr3fkYTSpAVNF6epU8uiCcov1lFebiAxq2/GO02UnCXEIwUXn2ux1rdAS5dKLoR6Ndy5wsnFikPtgDhbub3bmFVh2Gl0w2Y/cYj274lu2s6iiKNZLJY8u6FSaZcfMPiGX1/LYmb+D7zNXozf/uoCvxlRDUuU5+rhGX+TOCxvtPZZKUyXHS49esM6I3q74uuvYGqe6rhSlq1PJows6nVaJva2GHn6tb3mYpZkfMtewOWcDr5x6kaxqywExiRVnMElTo/GO1ujt0hcPnQd7Cy7cdaXRCH43xJPDieUUljXfQlEUpWtQyaMLOp1eRa9gR7SXMVieWplCpamS8T4TKTeW8+9Tz/H4sUf5JPkjdEJHuHPEZcWkERpGeY/hVNlJivUXbllcM8QDs7RsnqgoStelkkcXU2swcy6r6rK7rE6UHUcguCFwDo/3e5rJvtcQ7dqf/m4DuSlkPjqN7rJjG+U1BolkX+HeC9YJ9rGnb6gjWw8VqZXGitKFqUWCXcy5rGpMZugT6nRZ958ojaenUzhONpb7bwye12axedv50MulN3sLdjPVfzoa0fxnk2uHevLfNRkkZlUT1cnbySuKcnlUy6OLOZVmOaH3cloepYZS0qpSiXYbcOnKl2ms93gK9QWcLjt1wTrjB7qjsxH8pNZ8KEqXpZJHF3M6rRJ/D1s8XFrfvXSy9ARAuyaPQe5DcLZxZlfBjgvWcXGwYXRfN7YdLqZGb2q3WBRFaT8dljyEEB8JIfKEEM2enCKEmCSEKBVCHKn7eqrBtWlCiAQhRKIQ4u8dFbM1Op1eRZ/Qyx/vcNO5Eexw+ZuhXYpOo2O011iOlxyjRH/h9Ryzx/pQXm1i00HV+lCUrqgjWx7LgWmXqLNTSjm47utZACGEFngbmA70AxYKIS5vPmkXl1+qp6DUQJ+Q1o93mKSR02UniXYdgBBte/7Hb43xHo8ZM3sLd12wTr8wJ6J7OPHdLrXiXFG6og5LHlLKX4DL+Zg5AkiUUiZJKfXAKmB2mwbXRZw+vzjwMgbLE8oSqDZVE+3Wv63DasLX3pc+Ln3Znb8Ts7zwRojzJ/iSV2Lgl2Nq2q6idDXWNuYxWghxVAixQQhxfqlzENBwx72MurJmCSGWCCFihRCx3e3cgJOplehsBOEBrdtJ1yzNrM1ajYfOg34dkDwAxvlMoNhQzPeZ311wSu7w3q6E+dnz9S95atquonQx1pQ8DgFhUspBwH+BNZfzECnle1LKGClljI+PT5sG2JlMZsnO4yUMiXBBZ9O6v7b9hftIr0pjTvBcbDUdc975IPchjPeZyNbczXyT8WWzyUGjEcyb4EtKTg2xZ8o7JC5FUdqG1SQPKWWZlLKi7vfrAZ0QwhvIBEIaVA2uK7uqHDpr2dJjSoxnq+6rMdXwQ9ZqejqFM8xjeDtF15RGaLg55BZ+53st2/O28V3GN83WmzTIA2d7LbvVZomK0qVYTfIQQviLupFcIcQILLEVAgeBKCFETyGELbAAWNt5kXaOrXFFuDpqGdmn+R1vL2RLzkZKDaXMDf59uw+U/5YQgpuC5zPGexzb836ioLZpN6KNVjAowpnDieWq60pRupCOnKq7EtgL9BZCZAgh7hZC3CeEuK+uyjwgXghxFPgPsEBaGIEHgE3AKeArKeWJjorbGpRXG9lzspTJgz1a1WX1S952NuasZ4TnKHo6h7djhBcmhOD6gBvQCA1bcjY3W2dolAt5JQYyC2o7ODpFUS5Xh21PIqVceInrbwFvXeDaemB9e8TVFew4UoLRJLl2WMu7rDbnbOD7zNUMdBvELWG3t2N0l+Zu68FIr9HsK9zNjMDrcdO5N7o+NMoFgLiz5QT7tN2xuoqitB+r6bZSLmzLoSLCA+yJbOHJgQcK9/N95mpiPEdwT8S9V7TZYVu51m8qJmliW+7WJtf8Pe0I9LLl0Fk1aK4oXYVKHlYut7iWMxlVXDOkZa0OszSzMftHgh1CWNRjMVphHXtf+tr7MtQjhp35O6gyVja5PjTKlWNJFRiMF14XoiiK9VDJw8ql5lrGAVq6MPBYyRFya3O57iK72naW6/ynUWuuZXdB05XnQyKdqdGb609JVBTFulnXTxelicyCGgCCvO0uWVdKyaacDXjb+TDEY2h7h9ZqwY4hhDtFsLdgV5OZVYMiXNBoUF1XitJFqORh5TILanG21+LqqK0vM5qNvHzqRfYV7mlUN6H8NGlVqUzxm2p1rY7zxniPJ7c2l3MViY3Kney19Alx4tDZsk6KTFGU1rDOnzBKvcyCWoK87Rqt0cipySa1KoWVqZ+TWZUBWBLK+uwfcdO5MdJrdGeFe0lDPYZhr7Fnd8HOJtdG9nXlbGY1m2ILOyEyRVFaQyUPK5dVWNukyyqz2pIwbIQNHya/R7G+iLfOvsG5irNcHzDLKmZXXYid1o4Yz5EcLo6jyth4fOOmcb4Mi3LhP6vTOXBatUAUxZqp5GHFag1m8koMTZJHRlUGOqHjDxH3kVeTyz/j/0FyZRKLeixmrM/4Toq25cZ6j8MgDfySv509Bbv4IOl/pFamYKMVPH5rDyICHHjxixTOZKjBc0WxVip5WLGsQstMqyCf37Y80glwCKSPaz9mBc3BVefGn3v9lRFeozojzFYLcQwl2CGEH7LWsCL1Uw4Xx7E5ZwMAjnZanlkUjpuTlle/SkWvpu4qilVSycOKnd+uI8jr1+QhpSSzOoOgutMAr/OfzrP9XyTcOaJTYrwcQgjmhy5gmv8MHu3zOJN9r+F46TEq69Z/eLjo+NONIaTn1/Llz7mdHK2iKM1RycOKnU8egQ26rcqMpVQYKwhy/PUo2Y7e8LAtRDpHcUPQHMKcejDSazQmaSKu+GD99WG9XJk82IOvduSRmlvdiZEqitIclTysWGZBLZ4uNjja/TpNN6NudlV7nkPe0YIdQgi0D+RA4b5G5UuuD8LRTsOb36VjNqsddxXFmqjkYcXOT9NtVFY30yqoGyUPIQQjvEaRXJlEXs2v3VTuzjYsui6AU2lVavBcUayMSh5WrNnkUZWBh60njjatP8fcmg33HIlAcKBof6Py8QPd0Whg36nSTopMUZTmqORhpSqqjZRWGgnybrxFeUZ1erfqsjrP3daD3i592FewmzLDr2s8XBxsGNDTmb0n1boPRbEmKnlYqcwCPdB4TyuD2UBeTW636rJqaGrADCqMFbx86gXSqlLry0f1dSMtr4YsdViUolgNlTysVHMbImZXZ2HG3GimVXfSy6U3f+nzKACvnX6ZU2UnARjdz3L07l7VdaUoVkMlDyuVWVCLRoC/p219WUZ1OgBBDiGdFVa7C3UM49G+T+Bp68k36V9ilmb8POzo6W/PvpMqeSiKtejIM8w/EkLkCSHiL3D9ViHEMSHEcSHEHiHEoAbXUurKjwghYjsq5s6UWVCLj7sttg3OLE8sP4utxhYfO59OjKz9uepcuT5wFjk12RwpOQzAqH5unEytpLTS2MnRKYoCHdvyWA5Mu8j1ZGCilHIA8Bzw3m+uT5ZSDpZSxrRTfFYlr0TfqNWRWH6W/UV7Ges93mq3W29LQzyG4Wvnx8bsdUgpGdXXDbOE3fElnR2aoihAh51RKqX8RQjR4yLXGx5OsQ/onh37LZRfamBQuDMAerOeFamf4GXrzQ2Bczo5so6hERqmBkzns5TlxJcep3/QAPw9bfnvmgy+3J5LvzBnQn3tCPCyY1gvF1wcrOO4XUW5WljrR9i7gQ0NXktgsxAiTgix5GI3CiGWCCFihRCx+fn57RpkezGZJEXlBnzcLS2PdVlryavN45aw27HTXvpEwe5iuOcIvGy92JD9IwCvLIlk6Q1B9A5xIj6lgk+35PDSqlSe/DipkyNVlKuP1X1cE0JMxpI8xjUoHielzBRC+AJbhBCnpZS/NHe/lPI96rq8YmJiuuSeFkUVBsxm8HHTUVBbwE+5WxjrPZ4+rn07O7QOpRU2TA2YwRepn7Ej/2cm+f6OWWN8mDXGMuZTozexenc+n27OITm7mp4BDp0csaJcPayq5SGEGAh8AMyWUtYfJyelzKz7NQ9YDYzonAg7RkGJAQBvN1sOF8cikUz1n97JUXWO0V5j6e82kO8yvia1MqXRNXtbLdeP8MZGK9gcV9Q5ASrKVcpqkocQIhT4DrhdSnmmQbmTEMLl/O+B64BmZ2x1FwVlluTh46bjcPEhQh3D8LLz7uSoOodGaLi9x5242rjyYdJ7TU4fdHWyYVRfV7YdKcKgzv5QlA7TkVN1VwJ7gd5CiAwhxN1CiPuEEPfVVXkK8ALe+c2UXD9glxDiKHAAWCel3NhRcXeG/BLL6nIbp3JSq1IY4jG0kyPqXM42ziwOX0Kxvog1md82uX5djBdllSZ1dK2idISzEqgAACAASURBVKCOnG218BLX7wHuaaY8CRjU9I7uK7/UgL2thrPVxwAY7H51Jw+AcOcIRnuP5WDRfm4MnouD1rH+2tAoF7xcdWyOK2Jsf/dOjFKxBkaTpFpvokZvxs3RBlud1XSwdCtWN2CuQEGpHm83HUdKDhPoEISvvV9nh2QVxnqPZ3fBTmKLDjLeZ2J9uVYjuGaIB9/8kkdWYS2BXlfPjDSlsWNJ5fzz02Sqay1dmIFetvx3We9GZ+IobUOlZCuUX2rAy6uWpIpE1epoINQxjCCHYPYU7GpybdoIL+x0Gh58+4zaxuQqlZ5Xw3OfpeDtquPemUEsnhZAdpGeTzdnd3Zo3ZJqeVihglIDERHJlCOv+vGOhoQQjPEex9fpq0ivSsPfPoBvM77CXefOtIDr+c8Dvfj3qlSe+SyZIZEu2NtqsNNpuO1a/ybnoihdX2mlkYMJZbg52eDhbMMLK1KwsRE8d1c4fh6Wv++8EgNr9xYwebAnvUMcL/FEpTVU8rAyxroFgqGuifja+RJgH9jZIVmVEZ4jWZPxLdtyt1JiKOFM+WlsNbZM9r2WYB97XlsaxWdbcog7Yxk8zynSk11Uy2v3RaHRdL2z3pXGpJRkFtTy474CNh4sotbw6ww7O53gpT9E1icOgDunBrDnRAn/WZ3Of+7vhVar/g20FZU8rExRuQEpJZW6dPq7DEcI9Y+9IUcbJwZ7DOVA0T40aJjoM5kd+T9zvPQoMZ4jsLXRcPf0QO6ebkm62w4X8cpXaWw8WMiMkVfndOeuymiS/POTJHJL9Lg726DTChIzqymvNqHVwOTBHswc5YPRZCa7SE94gAPhv1ko6mSv5Y+zgnl+RQoLXognMtCB3iFOjOzjSq8QR7TqA8VlU8nDyuSX6NG5FmEUNYQ7R3R2OFbpGr8ppFelcVPwfPq69uNIyWHiig4S49l07ejkwR5sji3i443ZjO7nhoeLrhMiVi7Hmt35xJ0tZ3hvV6r1JiqqTYzt70avYCeG9XLB1/3XjUOje1z4OWOi3XhsYRiHEys4l1XF17/k8uX2XNycbJgwwJ2pwz2JCPy1S0tKyepd+Xy/J5/eIU6M7ufGiD6uONmrQfeGVPKwMgVlBuy9cgAId4rs5GisU4hjKE9GP1P/ephHDL/kb6fKWIWjTeN+bSEED8wJZumbCby3LpNHbw5TrbkuILdYz+dbcxjV15Wn7wi/omcJIZgw0IMJAz0AyxHPsWfK2XuylI2xhfywr4CIQAcmDnRnRB9XvtiWyy/HSugT4kh8SgU7jxfj6Ghi2jBfZo72odiQR2JJOqFu/gz2u3r/j6rkYWXySyzJw0nr3O3P7WgrMZ4j2Ja3lSMlhxjjPa7J9WAfexZM9uPzrTmE+tmzcLJ/J0SpXEpKTjU6Gw0Bnra8+0MGAPfd0Pabazs72DBpkAeTBnlQXmXk5yPFbD1UxEcbs/loYzYaAXdNC2D+BF/OlJ1hedIKysyFHAGONNiDU5ZqGJC4kCWjJ1yV3V8qeViZglI9Dt45RLpEqk/ILRTqGIa3nQ9xRQebTR4ACyf7kV1Yy6ebc3BxsGHmKDX+YS2qa028vz6LDQcs29k52GmorjVz9/RA/DxsL3F388zSzI78nynWFzMz8AZsNY1n22VWZ7Ahex06oSOyVxQPDQoio7iCYxm5eLlrCPQqY1XaFnYV/IK3rTezfW6ktMpASk4l7jY+BDj6s7XkG47breThr2p4cPJEwvzsr/jPoitRycPK5FaUYBNYSk8nNd7RUkIIYjyGsylnA2WGMlx1rk3qaDSCB+eGUlFj4p21Gfh72BLTu2k9pWMkZ1eTUVBLbrGedfsLyC3WM3e8DyE+9pzJrMJslswZ27KWd7G+mLjig9SaaujpFIGLzoWv0laSVHkOgJOl8dwdvgRfez9ya3LZlf8LO/O3Y6+1Ryu0HCja9+vDdECl5UsgmOx7DTcEzvn1KIQGPWiTDT158dirFIev5rnjO3A9pyPaJ4w7Im5Hp+n+Y2sqeViZApkKQITz1duXejmGe41kU84G1mR+yx097mq2jo1W8NjCHix57RRr9xao5NEB8kr0bDxYSJCXHQMjnMkpsoxlHEuqqK8T6GXLy3+IpH9Py+FnU4d7tejZ6VVprM74hjPlCUgkAoHEcgqDg9aRO3rchavOjeXJH/LSqRcAMEgDAsF4n4nMDJyNo9aRvNpccmtycNG54q5zx0boMEkjWmGDi87lgu/vpHPi8UEP81XK1yQYCsgrqeaQ3T6O/VzEFOfb6RviQoCXXbdd3a6Sh5Wpsk3HQWoJcQzt7FC6FH/7AKYFXM+G7B/p7zaAoR7Nn1Zsp9MwYaAHq3flUV5lxMVR/RdoDyaz5Md9BSzflE2NvvFuxx4uNiy5PpDBES74uOtwstfWd9FmV2exPvtHEsvP4O8QQLBDCKO8xhDk+OvYh1ma2Zq7iR+z1uKkdWJ6wEyGe47AVedGSmUSOTU5DPEYipvOss/Z4/2eZH3Wj9hqbAl2DCHCORLvBuOJfvb++Nlf3jiYk40Td0XeCZGQVVDLpyfWke6+kS+Tv6BgzSRA4O2mY0BPZwaFOzOmv1u3OfWye3wX3YTBaAaXLJxNgVdFs7etTQ+YwcnSeFamfk64UwTuth7N1psw0J1vfsljz8lSpsa07FOu0nKpudW88W06p9OriOnlwv2zg6msMXEsqQKdVsO1wzyxt7XsjFRuKCe2+CR5NbmkV6UTX3oMW40tA9wGkV+bx878HfySv525Ib9nvPdEzlac4fvM70ipTGaI+1AWht2Ok41T/Xv3ce1HH9d+jeJx07mzMOy2dv++A73t+PvEm/gxU8sG1tE3XBBcPp3MLDsOJ5bz85FiPt2SzV/mhTKsV/OtXrNZklVYi85Gc9njPR1FSNklD9trkZiYGBkbG3vpilYis7CCF849TKRmLH+Jub2zw+mS8mpy+dep5/C182NuyHyinHs3mXggpWTxq6cI8rLj+cVqbKmtmEySlT/n8OX2PBztNNw7M4jJgz2anfhRaaxka+5mtudtQ2+uRSDwtPVkqMdwrvWfgrONpbuo3FDOpykfcbLsBL52vuTV5uGuc2d20E0M9xxplZNKpJT8lLeFdVlrkVIyNWAGv/O9lqRMI69/m0Z6Xi0zRnoxNtqdqGAHSsqNxJ4p43BiBafSKqmoNgHQO8SRyYM86vdtA8uEmvfWZTE1xvOCCehKCSHipJTNN90b1mtt8qg7kKlGSmm63OA6SldLHlvOHGNN+VtcY38nN0WP6exwuqyjJUf4Mm0FpYZSIpyjuC1sEb72vo3qfLwxi2925rHy8f64OqkGeFv4fk8+7/6QyeTBHiy5Pgh358Z/rmZp5kx5AgeL9nO4OA69Wc9Qjxiu8ZtCoEPQBVvb57updhfsYoLPJMb7TMRWY92fygGK9UV8k/4VR0oO4aZzZ0bATIa6jeaTjbn8sK+A3/7oDfGxI7qHE31CnSirNLL9aAlJ2dX09Lfn8Vt6APDER+fIKzFgoxX849YejOzr1uZxt1nyEEJogAXArcBwoBawAwqAdcD/pJSJVxxxO+hqyeM/cSs5bd7OX0L/RaSfZ2eH06UZzAb2FOxiXdZanGyceLjPY426NxKzqlj23zP8+aYQprVwgFa5MKNJsviVk/h62PLqvVGNrmVVZ7K/cC+xRQcoMZRgr7FniMcwJvtdQ5BD26/jsDaJFWf5PuM7kirPEekcxR8iliIN9iRmVnM2swoney0xvVzru6kKaguoNFbgqnMlIQle/zoTg1GisxEIIfjbgjA+3phFck4ND80NYVgvV1wdtW3WCmvL5LED2Ap8D8RLKc115Z7AZOAWYLWU8vMrjrqNdbXk8ej+pykogXd/97Q6wKaNnKtI5D9nXiPcOYIHov6MVlg+DUspuef/TuHnYceLd6uuqyu19VAR/9l0mBnX57Moelb9YV2HiuP4KOk9BIJot/4M9xzFAPeBXaLl0JaklBwo2scXqZ/haevJPeH3UWwo5lTZCSqNldhp7Cwts4oECmrz6++zETbM9r2NDet9KCw18Oyd4QT72FNRbeSJj5I4k2E5ltnZXkugtx2BXrYEetkR6GVHgJcdqbk17DlRwrmsamJ6u3LNEA8G9HS+6CahbZk8dIAJ+LuU8sUL1ZFSGi71Zh2tKyWPUkMpjx97hJqzY/hwwZ2dHU63sr9wL5+mfMw47wmNBk4/2ZzNV9tzeetPvenp73CRJygXYzAZWbZ6BZqwvQiNmV4uvflj5J8o1hfz0qnn8XcI4L6IBy467fVqkVhxlvcS36HSVAmATuhw1bmhN9dilmZ6OkfQ17UfnraelBvK2Vu4m/SqNJZFPURPp8hGK9lrDWYOny0ju0hPVqGezIJasgtrySvRY27wYz3A05aIQAfizpZTXWvG113H74Z4cs0QD4J9mi5sbPMxDyHEHinlFXXECyE+AmYCeVLK/s1cF8CbwAygCrhTSnmo7toi4B91VZ+XUn5yqffrSsljb8EePk9djm38Hby+qPlV0srlW5PxLVtyN3FvxB8Z6D4YsJwHce/rp/F11/H6UrVd9+UwmA386+hr5JrPEcxARof04+v0VQz1iCGvJpdifRF/7/cPPG1V1+B5+bV5xBXFEuYURqRzr4vOrKwwVvB/p/9NlamKR/o8Vj/F2CzNrEz7nAOF+/Cx8yXIIZhR3qPp6xqN3mgmt+4oAh83W3r42yOEoEZvZt+pUn46VMShs+WYJdw/K5iZo3/dbaG61oSjvU2LkkdrRgqPCSGeBp4733V1GZYDbwGfXuD6dCCq7msk8P+AkXVdZE8DMYAE4oQQa6WUxZcZh9U5WRaPudYJf/vu3wfcGWYGzuZk2QlWpa0g0jkKRxsn3JxsWDoriH+vTGX17nzmTfC99IOuQjV6E2l5tUQGOqDRCCprTHyxLYcdR0tw6PcTtsHnqD05hUdvmYdWKzCY9azJ/A6ApZEPqMTxGz52vkwLmNGius42ziyNXMYrp//F6wmvMi/k9wx0H8ynKR8TW3SAYR7DqTXXklB+mtjiAwzxGMZNwfMI9HEn2Meu0TiIva2mfk+vojID/12TzttrMxBaM9cM8eDHvcV8sS2nxd9Ha5KHJzARWCqE2A8cA45JKb9u6QOklL8IIXpcpMps4FNpaQ7tE0K4CyECgEnAFillEYAQYgswDVjZivitlkmaOFV2ksrsUAI81Il37cFGY8NtPRbxyql/8U3GV/Wr0CcMcGfH0WI+25LNqL6uzTbjr2Zms+SFFSnEninH08WGmN6u7D9VRlmVkUFjkigLiMe9ZCzzRk+rb7ld6zcVENhp7ejvNrBzv4FuwNfejweiHuTz1E/4IOl/uOs8KDEUMyvwRqYGTAcsLcCtuZvYmL2ew8Vxv95r50e0W3/6uw2gl0sfNMIyllqry6fn+F8oDD/HjxSw7pCWysJQeg2KbnFcLU4eUsrfAwgh7IBoYACW1kGLk0cLBAHpDV5n1JVdqLwJIcQSYAlAaGjXWKWdXJlEtamKquww/EKvroHEjhTqGMZ1/tPYmLMeb1sfJvlOxtHGiftnh3Df66d57vMUXl4SiZuaultv7d4CYs+Uc8MobwrLDWw/UkxUsCPXXlvJ+tIt9HPpz9Kht9f/UALLXmNT/Kd2YtTdT5hTD/7e9x/sLtjJ5pyN3Bg8j2v9rqu/rtPomB4wkxjPERwpPoRRmjBJI6mVKezM38HPeT/hZ+fHZL9ryavJY3veT9hqbIn0CSctJYxKYyXeYSmUt2Li7CX/lwghhGwwMCKlrAUO1X01W6czSSnfA94Dy5hHJ4fTIidL4xFoqMoNxs9TJY/2NC3gejKq01mXvZYtuRsZ5z2BG4Pn8eTtPXjy4yT+8dE5/v2HSHXwD5CUXc2HG7IY1deVpbOCEEJQZahmddbX/FCwixCHEO7qeXejxKG0H63QMsFnEhN8Jl2wjo+dL1P8pzUq05trOVZylK25m1mVtgKBYIz3OG4InIOLzgXZSyKEwCzN7Cvcwzu836J4WvIR62chxLfA91LKtPOFQghbYBywCPgZy3jGlcoEQhq8Dq4ry8TSddWwfHsbvJ9VOF12Cg8ZijTaWf2WBF2dTqNjaeQyMqrS+Sl3C9vytuJo48j08Jk8cWtPnv0siaeWJ/HkbT1wd776toj5MWst+wv34qRxITVDi1d4FH+6aSZCCE6WWsaMivSFXOc/jRkBN6htdLoAW40dMZ4jGOYxnOTKJOy19gQ6/Npxc35cRCM0FzzSoDktSR7TgMXASiFEOFAMOAAaYDPwhpTycIvf8eLWAg8IIVZh6RIrlVJmCyE2AS8KIc5vVnQd8FgbvWenMkkTmdUZeNeMQCPAx00lj44Q7BhSP+6xLusHIp2jGNGnN39b0IOXv0zlvjcSeGBOMJGBDvx0uJgz6VXMn+RL/x7OnRx5+9met40N2T8S4dSL5Cw9Zoc8nAad4/+lxeNr58vhkkP42fnxYO+HiXSOuvQDFasihGjTo60vmTyklDXAO8A7dWs+vIFqKWVJa99MCLESSwvCWwiRgWUGla7ufd4F1mOZppuIZaruXXXXioQQzwEH6x717PnB864utyYXozRiKPPG202HjZou2mGEENwcegsplcl8nPwhj/V9kvED3AnxseP/vk7jhRUpdfUsi7Ae++AcS2cFMWNE9zhIKqeotv5c7xT9SarDvyNUF03erukkJVXxt4Wh2AeeZV32Wo6XHmNGwEyu85+uWhsK0Lp1HmeB48BRsJzIKGXd4RNWqius8zhYtJ/lyR9iG78IXa0vLy9Rn+g6WkZVOq+c/he9XHqzNHIZGqHBaJKs219Ajd7M5MEeONppeGlVKrFnyonu4YSHsw2uTjbMHe9LoFfnzJCrqjWRXVhLqK89OpuLjzuUVxkpKjeiEVBUbuCHvQXsOVmKlOAenI7n8I3oy9zJ3D4HadLxl3mhTBlm2SLHVDf4+tvT+JTuqaWLBFszreR/WM7RKsSyHmOFECIZWI1l7YfVrTDvCjKrMrARNuRnuzA0Uv3n7AzBjiHMC7mZVWkr+Cl3M1P8p2GjFcwe0/gku38uCmfVz7kcTCgjLa+W3OJydseX8syicHqHOLZ7nHklemITyjiUWM7ptCoKyyz/5SYNcufRm8Oa7G1kNEmOnitnU2wRe0+WYjT9+kHR2UHL/Ik+uEYe5eeSHwl0CGJxv/tJCbTBRisaHZSlFVq0Qk0gUBprTfK4TUo5+PwLIcS7WMZCyoDXgGVtHNtVIaM6Az+7AM6UmvFzV+MdnWWc9wQSyk+zNnMNEc5RzfYNazWCW6/x59ZrLAcHZeTX8I+Pk/jb+4k8cWsPhl/iZMKKaiOn0ixnXLR2E7u4M2U8/UkSJjP4uOkYFO5MmJ89eSV61u0vZGiUK1OGeSKlZMfREn45XsLRc+VU1ZpxddRywyhveoc6IiXotILQsCo25K1hW8kRhrgP5fYed2GntcO/36VjURRoXfIoFUIMlFIeA5BSHhFCTJRSDhJCHLrUzUrzMqvSCbPrA6Cm6XYiIQS3ht1OWmUqHyW/z9/7PlF/psSFBPvY89rSKJ5ansQLK5J58/7ehPk1v8hQbzTz9CfJnEyt5Lk7w1t1BG5KTjUvfJFCqK89jy3sgZu7npzabCKdQ5BSkJZXwztrMwjwtOWrHXkcTCjD113HxEEexPRyIaa3K7Y2GgxmA8mVScQWHWDVmd3YaHTMCbqJa/2mWuW5GIp1a03yuBdLV9URLGMevbEMagOon3qXocxQRpmxDAet5ZOsmqbbuRy0jtwdvoTXE17hvXPvsizqwUsODnu66Hj2znD++GYCL3+ZwutLezXZEVlKydvfZ3AytRJnBy3LN2czNMrlojubSimpqjWTX6Ln6U+ScLDV8MyicDzdtPzf6ddIrUrBy9aLsd4TeGDeaB5+K51H3kvE1kZw3w1B3DDKmypzBcdLjrM2O53M6gxSKpPRm/VohZbxPpOYHjADF506x125PK1ZYX5aCDECuAkYiGVG1NN1h0Otaqf4urXMasuieVFl6Vv3V8mj04U59eC2Hov4OPkDVqWt4LawRZf8VO7pouMv80J4+pNkPt6UzaLr/DmRUklGfi22Og2ZBTVsji1i4WQ/grztePXrNHafKGX8APcmzyqvMvLVjjx+3FdQf/a3nU7Dq/dG4uNuy8bsdaRWpXCt33WkVaWxNms1221+4ua5N3LycCC3XutFtUMqy1NWc7TkCEZpRCd0BDoEMdprLH1c+xLl0hsHrdpFWLkyrdqHoe70wK9puiXJ820W0VUkoyoDgJpiL2y0lXi6qimQ1iDGcwS5NTmsz/6RUkMpXrZeuOncmeA7sdmurBpTDb3Czcwa7c2a3fn8sDcfU4OtQ7V2VYzp781t1/ojga925PHZlmwGRzizenc+O44W42SvxdNVR3xyBVW1ZiYOdCcqyBEXRy19QpwI8bUnvSqd9dk/MtQjhhuD5wGQVpnKF6mfsq36EwKGBPJmVh5GacRR68g4nwmM8hpLkEOQWgWutDm1iU8nyqzOwF3nTmGhDT5uukZ79Suda0bADdSYajlRdpz0qjQqjRXszN/ObT3uJNrt19ME0qpS+eDcu5QaSpkzYj61xgjcHG0YGO5MRKADJ8ri+SpnBZVae7bmTWG8z0TuuM6f5z9P4bZ/nUBvlAyNckEAuUV6BkW4cNu1/k3OFykzlPFp8oc4aZ24OfSW+vJQpzAe6fsYP+Vu5UTpcaLdfkdvl75EuVx8q29FuVKtPsO8K7H2dR4vnPgnnnZeJGyaioOdhn/dE9nZISkXkFGVzvLkD8muyWKA20B6OIWjERrWZa3F2cYFP3s/EspPM8R9KNMCrifAIZCjJYdZnvwhgfaBOOtcOFV2EgetIxN9JhO3IxJpcOD2a/2JDGo6zbfCWI5W2GCvsedQcRxfpq2g1lzLkog/NkpeitLW2mOdh9KGDGYDOTU5DHAfxM5S/SWneSqdK9gxhL/1fYJ1WT9wuDiW46XHAOjj0pc7e96Dk40TP+VuZm3mGg6XHMJWY4fBrKenUzh/jFqGg9aR1MoUNuVsYGPOOmz72jEv+PdE+oQ3eh8pJZtyNvBD1hrActKcQRoIc+zB7T3uJMAhsMO/d0VpjkoenSSnJhszZvzsgiguN+Kr9rSyejqNjjnBNzEn+CaqTVUU6YsIsA+sH0+Y4j+NGM8RJFacJbkiGY0Q3BA4BzutZfFnmFMPlkQsJbs6i6/TV7EqbQWBjkH0dLIkEJM08VXaKnYV7GCIxzB6OPakzFiKp60X430mqoV6ilVRyaOTpFQmA+BsCgAK8XZX/dNdiYPWkSCHpt1NHraeDPccyXDPkRe8N8AhkHvC7+Nfp57j46T3+XvfJ6kyVbEy9TNOl5/iOv9p3BA4Rw1yK1ZNJY9OcqL0OF62Xhgr3IBCfNXq8quKo40jd/W8h9cTXuG/Z18nuzoLjdBwS9jtjPUe39nhKcolqeTRCQxmAwnlpxnlNYaCAiNg2XJCubqEO0dwQ9Bsvs9czSD3wcwPWYiHrcelb1QUK6CSRyc4U56A3qynv9sAjp3TA+CtxjyuSlP8pjHMYzhedt1jm3fl6qE6VTvBidJ4dEJHlEtvCkoMuDppsbdVfxVXIyGEShxKl6R+YnUwKSUnSo/R27UPthpb8kr0aqaVoihdjkoeHSy3NocCfQHRbgMAyC814KNmWimK0sWo5NHBTpTGAxDtWpc8SvTq3HJFUbqcDk0eQohpQogEIUSiEOLvzVx/XQhxpO7rjBCipME1U4Nrazsy7rYUX3qcAPtAvOy8qKwxUVVrxkdN01UUpYvpsNlWQggt8DYwBcgADgoh1kopT56vI6V8qEH9ZcCQBo+obniSYVdUpC/kbHkCU/2nA5ZWB6C6rRRF6XI6suUxAkiUUiZJKfVYzgCZfZH6C4GVHRJZB9mdvwuAMd7jAMgrsZxBrQbMFUXpajoyeQQB6Q1eZ9SVNSGECAN6AtsaFNsLIWKFEPuEEHMu9CZCiCV19WLz8/PbIu42YZJG9hTsJNqtf/3UzPxS1fJQFKVrstYB8wXAN3WHT50XVrdN8C3AG0KIiOZulFK+J6WMkVLG+Pj4dESsLXK05AhlxjLG+0ysL8svMaDVgIeLSh6KonQtHZk8MoGQBq+D68qas4DfdFlJKTPrfk0CttN4PMTq7czfgaetF/1cfz2LIb9Ej5erOgRKUZSupyOTx0EgSgjRUwhhiyVBNJk1JYToA3gAexuUeQgh7Op+7w2MBU7+9l5rlVOTzZnyBMZ6j2+0U2peqV5tiKgoSpfUYclDSmkEHgA2AaeAr6SUJ4QQzwohZjWougBYJRsfcdgXiBVCHAV+Bv7dcJaWtdtTsAsNmvqB8vMKSgxqTytFUbqkDt0YUUq5Hlj/m7KnfvP6n83ctwcY0K7BtROTNHKgcB8D3Afhqvv1tECzWVJQZsBXDZYritIFWeuAebcRXxpPubGc0V5jGpWXVBgxmqRaIKgoSpekkkc721ewGxcbV/q59W9Unnd+gaA6x0NRlC5IJY92VGYoI770OCO9RjU5fzqn2JI81IC5oihdkUoe7ehA0T7MmBntNbbJtbMZVdjaCIJ97DohMkVRlCujkkc72lewh55O4fg7BDS5djK1kl7Bjuhs1F+Boihdj/rJ1U6K9UVk12Qx1COmybVag5nErGr6hjl1QmSKoihXTiWPdpJUkQRAhHNkk2v/v717D47rLO84/n1WlmRdfJNt2cSSbPmSEMcJuajh2hQSAgmlTqdlwFBmAqVNmRKg5Y8OmXYopO0MtPzRC2mHTEgboMEkpsy4EEigDiktuSkmJfFdF+saW7ZkYu2uLt7dp3/scZClla3Vrs56j3+fGY13z3uO9zn2SD+d877nfY8MJEmlna0tCg8RKU8KjwXSneik0ippqm2a0XagJwHAEztkWQAAEH1JREFU6xUeIlKmFB4LpCveyYa6Vips5nOYB3oTrFtZzfL6UJ/RFBEpGoXHApjMTNKX7KW1fubEv+7O/p4kV66vLUFlIiLFofBYAL2JHjJk2Fg3MzxeGZ7k1URKneUiUtYUHgugK9EJQGt964y2/b3Z/g51lotIOVN4LICueCdrqtdQv2jJjLYDPQlqq2O0NC4uQWUiIsWh8Cgyd6c70ZmzvwOyVx5XttQR0wJQIlLGFB5FdmJiiHgqzsYc4XEmlaF3aJzN69RZLiLlTeFRZGf7O3J1lg+cnCCTgfVrdMtKRMqbwqPIuuPd1FTUsGbx2hltvUPjALQ0ajJEESlvCo8i60v20Fzbcs5a5Wf1HB8nZtC0WlceIlLeQg0PM7vNzA6ZWYeZfTZH+0fM7ISZvRh8/cGUtjvN7EjwdWeYdc9V2lMMjPXTXNuSs71naJy1DVVUVyqzRaS8hTY/hplVAPcBtwL9wPNmttvd90/b9dvufve0YxuAvwTaAAdeCI49FULpc/bK2CukPEVL7fqc7b1D4+rvEJFICPNX4BuBDnfvcvdJYCdwxxyPfTfwI3cfCQLjR8BtC1TnvPUmewBozhEeZ1IZBk9O0NJYE3ZZIiJFF2Z4rAP6przvD7ZN97tm9gsz22VmzXkeW1J9yV4Wxxazunr1jLbB4QnSGXWWi0g0XGw33/8T2ODu15C9ungo37/AzO4ys3Yzaz9x4kTRCzyfvmQvTbXNs3aWg4bpikg0hBkeA0DzlPdNwbbXuPuwu08Ebx8AbpjrsVP+jvvdvc3d21avnnkFsFDSnqY/2TdrZ3nvkEZaiUh0hBkezwNbzKzVzKqAHcDuqTuY2dTFvrcDB4LXjwPvMrMVZrYCeFew7aJxfPwYZ/zMrJ3lPcc10kpEoiO00VbunjKzu8n+0K8AHnT3fWZ2L9Du7ruBT5nZdiAFjAAfCY4dMbO/IhtAAPe6+0hYtc/F2c7yWcNjaFyTIYpIZIS6lJ27PwY8Nm3b56a8vge4Z5ZjHwQeXNACC9Cb6KUqVk3j4jUz2s6OtHrL1mUlqExEpPh0D6VI+pI9NNU05ews/9VIK115iEg0KDyKIOMZ+sf6aKmb7eHA7BiAFo20EpGIUHgUwStjg0xmJmft7+gcHCMWg2aNtBKRiFB4FEFnogOATfWbc7Yf7k+wYU2NRlqJSGTop1kRdMY7WFa5jJVVq2a0ZTLOob4kr2/WAlAiEh0KjyLojHewsX4zZjOXlu0/OUFyIsMVCg8RiRCFR4FGJoc5NTky6y2rQ31JAIWHiESKwqNAXfHssrOzh0eCmuqYpiURkUhReBSoM95BdayadTVNOdsP9SW5vKmWitjMW1oiIuVK4VGgzngHrXUbqbCKGW0TZzJ0HxvjiibdshKRaFF4FGAsnWRwbIBN9VtytncOjpHOqL9DRKJH4VGArngXjp+3vwPg9c11YZYlIrLgFB4F6Ip3ECPGhrrWnO2H+pKsWlZJw9LKkCsTEVlYCo8CdMY7aKptproi99Kyh/qSumUlIpGk8JintKfpSR6ltW5TzvaR0TMcOzWpW1YiEkkKj3kaHBtgMjPJxvqNOdv3Hc32d1zdqvAQkehReMxTd7wLgNa63OHxUnec6soYmy7TbSsRiR6Fxzx1JTpZumgpDVUrc7a/3B1n6/paFlXo4UARiZ5Qw8PMbjOzQ2bWYWafzdH+GTPbb2a/MLP/MrP1U9rSZvZi8LU7zLpz6U500Vq/KedkiKNjKY4eH2fbhvoSVCYisvBCCw8zqwDuA24HtgIfNLOt03b7OdDm7tcAu4C/ndI25u7XBl/bQyl6FqNnTnNy4sSst6z2H03gDttaFR4iEk1hXnncCHS4e5e7TwI7gTum7uDuT7p7Mnj7DJB7wqgS604E/R2zdJa/fDTBogrTMF0Riawww2Md0DflfX+wbTYfA34w5f1iM2s3s2fM7LcXosC56k50ESM267KzL3XHuaK5VisHikhkLSp1AbmY2YeBNuA3pmxe7+4DZrYR2GNmL7l7Z45j7wLuAmhpaVmQ+rrjXTTXtlAVq5rRNj6ZpmMgyftualyQzxYRuRiE+avxANA85X1TsO0cZvZO4M+B7e4+cXa7uw8Ef3YBPwGuy/Uh7n6/u7e5e9vq1auLV33g7MOBG+tzPxx4oDdJOoM6y0Uk0sIMj+eBLWbWamZVwA7gnFFTZnYd8FWywTE0ZfsKM6sOXq8C3grsD63yKQaS/UxmJs/7fEfM4Mr1ejhQRKIrtNtW7p4ys7uBx4EK4EF332dm9wLt7r4b+DugHng0GALbG4ysuhL4qpllyAbeF929JOFxJH4YYNZp2PceGeWK5lrqFs9c30NEJCpC7fNw98eAx6Zt+9yU1++c5bifAVcvbHVzc2T0EI3Va1hetXxG2+lEisP9SX7vlrUlqExEJDwaDpSHjGfoiB9hy5LLc7bvPTKKO7RdviTkykREwqXwyMPAWD9j6bFZw6P98GmW1laweZ2e7xCRaFN45OHI6CEAttTPDI9Mxtl7ZJTrtyyhIqb5rEQk2hQeeTg8epjV1Y0sr1oxo63r2Bin4iluuHxpCSoTEQmXwmOOLtTf0X5oFIDrt6i/Q0SiT+ExRwNjA4ylkzlvWQG8cPg0my6roWGJ1isXkehTeMzRa/0dOa48fhk/w/7ehEZZicglQ+ExR4dHD7KqahUrqhpmtD3630PgcMt1M9tERKJI4TEH8VSc/af3cfXyN8xoGzl9hu89fZJ3XLuC5sbFJahORCR8Co85aB95jrSnedPKt8xoe+Sp46Qyzodu1lPlInLpUHjMwbPDT9NU00xTbfM520+8Osn3nx3m1usbuGxVdYmqExEJn8LjAgbHBuhN9uS86tj55HEAPqirDhG5xCg8LuCZ4aeJEaOt4cZzth8/NcET7SO8u62BNStmLgolIhJlCo/zSHua54afYduyq1lSee4w3If3HMcMdrxjTYmqExEpHYXHefz0xFOMpk7zxmm3rAZPTvDjvSO858aVrFqmqw4RufQoPGZx8PR+vtP3CFct3cY104boPrznGJUVxvvfrqsOEbk0KTxyOD5+jAe67mft4rV8dOMfErNf/TPtOxrnyRdP8VtvXq2pSETkkhXqSoIXu7Sn+NnJ/+X7g7upsAo+vvluaipqXmsfHJ7g3m92s7ahive/vbGElYqIlJbCI3Dg9D4e7d3J8YnjbKrfzAdaPsTK6lWvtY+Opfj8Q114Br5w50aW1OifTkQuXaH+BDSz24B/ACqAB9z9i9Paq4GvAzcAw8AH3P1o0HYP8DEgDXzK3R8vRk2jZ06zq/8R2keeo7G6kbs2/THXLHsDZsariRQ/eG6Yg30JDvYmSYyn+Zvf30TTak1DIiKXttDCw8wqgPuAW4F+4Hkz2+3u+6fs9jHglLtvNrMdwJeAD5jZVmAHcBVwGfBjM7vc3dPzrefY+Cs8NfQkzw4/TcpT3P669/LutbdTGaskk3Eebx/mwR8MEh9P07SqmrbLl3DzdQ1cs7F+vh8pIhIZYV553Ah0uHsXgJntBO4ApobHHcDng9e7gK+YmQXbd7r7BNBtZh3B3/f0+T6wd3SQT/7kC3gshVcmYFESLAMYmOOZCiYGtxA/cgMPJRv4hh0kFjPSaWd0LM221jruvqOJ9WtqzvcxIiKXnDDDYx3QN+V9P/DG2fZx95SZvQqsDLY/M+3Ydbk+xMzuAu4K3sa/8o7PH5pvwT8EvvxHF9xtFXByvp9xEYrS+UTpXCBa5xOlc4Fonc8Vc9kpcr2+7n4/cH9Yn2dm7e7eFtbnLbQonU+UzgWidT5ROheI1vmYWftc9gvzOY8BYOq0tE3Btpz7mNkiYBnZjvO5HCsiIiEJMzyeB7aYWauZVZHtAN89bZ/dwJ3B6/cBe9zdg+07zKzazFqBLcBzIdUtIiLThHbbKujDuBt4nOxQ3QfdfZ+Z3Qu0u/tu4GvAN4IO8RGyAUOw3yNkO9dTwCcKGWlVZKHdIgtJlM4nSucC0TqfKJ0LROt85nQulv3FXkREZO40t5WIiORN4SEiInlTeMyTmS02s+fM7P/MbJ+ZfaHUNRXKzCrM7Odm9r1S11IoMztqZi+Z2YtzHXp4sTKz5Wa2y8wOmtkBM3tzqWuaLzO7Ivg/Oft12sz+pNR1zZeZ/Wnw/f+ymX3LzMp67iIz+3RwLvsu9P+iPo95Cp58r3P3uJlVAv8DfNrdn7nAoRctM/sM0AYsdff3lrqeQpjZUaDN3cv+wS0zewj4qbs/EIxUrHX3X5a6rkIFUxYNAG90955S15MvM1tH9vt+q7uPBYN6HnP3fyttZfNjZtuAnWRn75gk+5z0x929I9f+uvKYJ8+KB28rg6+yTWIzawJ+E3ig1LXIr5jZMuAmsiMRcffJKARH4BagsxyDY4pFQE3wXFotMFjiegpxJfCsuyfdPQU8BfzObDsrPAoQ3OZ5ERgCfuTuz5a6pgL8PfBnQKbUhRSJA0+Y2QvBlDXlqhU4AfxrcEvxATOrK3VRRbID+Fapi5gvdx8Avgz0Aq8Ar7r7E6WtqiAvA79uZivNrBZ4D+c+nH0OhUcB3D3t7teSfeL9xuCyr+yY2XuBIXd/odS1FNHb3P164HbgE2Z2U6kLmqdFwPXAv7j7dUAC+GxpSypccPttO/BoqWuZLzNbQXbS1lays33XmdmHS1vV/Ln7AbIzmT9B9pbVi2SXwMhJ4VEEwW2EJ4HbSl3LPL0V2B70E+wEbjazb5a2pMIEvxXi7kPAd8nexy1H/UD/lKvaXWTDpNzdDux19+OlLqQA7wS63f2Eu58B/gN4S4lrKoi7f83db3D3m4BTwOHZ9lV4zJOZrTaz5cHrGrLrlBwsbVXz4+73uHuTu28geythj7uX7W9QZlZnZkvOvgbeRfaSvOy4+zGgz8zOznR6C+cuY1CuPkgZ37IK9AJvMrPaYADNLcCBEtdUEDNrDP5sIdvf8fBs+0ZuVt0QvQ54KBgxEgMecfeyH+IaEWuA72a/n1kEPOzuPyxtSQX5JPDvwa2eLuCjJa6nIEGg3wpceMGDi5i7P2tmu4C9ZKdN+jnlP03Jd8xsJXCG7DRQsw7O0FBdERHJm25biYhI3hQeIiKSN4WHiIjkTeEhIiJ5U3iIiEjeFB4iIpI3hYdIiZnZP5nZXjP7tVLXIjJXCg+REgoemGsk+8BcWU+DL5cWhYdICMxsg5mNBbMwv8bdE2RnK/gJ8I/BvjXBQkmTZrYq/GpFLkzhIbIALGv691dnMAvz1P1Wkl0HYpTsFBe4+1iwXzmvDSERp/AQKZLg6uKQmX2d7ESMs66FMMVfkF0TYh9w1ULWJ1JMCg+R4toC/LO7X3WhFfLMbAPZKby/TXY2VoWHlA2Fh0hx9eSxjv1fA/d6dnZShYeUFU3JLlJcibnsZGbXkl0v4W1mdh+wGHhpIQsTKSaFh0hpfAnY7u4/BjCzNWTXgxApC7ptJRIyM7sZqD0bHADBcqz1ZtZQuspE5k5XHiJF4u5HgW1z2G8PsCfH9qULUJbIgtCVh0g40sCy6Q8J5nL2IUGgEsgseGUi86BlaEVEJG+68hARkbwpPEREJG8KDxERyZvCQ0RE8qbwEBGRvCk8REQkbwoPERHJm8JDRETy9v9fRc5SsKEGRgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Plot RDFs between Mapped CG beads (center of masses) and the reference (C-O in all atom system )\n",
    "\n",
    "matplotlib.style.use('seaborn-muted')\n",
    "plt.plot(train_vars['rdf-r'],train_vars['avg-aa-rdf']/ train_vars['avg-aa-rdf'][-1],label ='All Atom C-C')\n",
    "plt.plot(train_vars['rdf-r'], train_vars['avg-rdf'] / train_vars['avg-rdf'][-1], label='Mapped (CG)')\n",
    "plt.xlim(2.5,9)\n",
    "plt.ylim(0,2)\n",
    "plt.xlabel(r'r [$\\AA$]')\n",
    "plt.ylabel('$g(r)$')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}