{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 02. Preparing Coarse-grained Mapped Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prequisites\n",
    "\n",
    "This code requires the `gsd` package, which is available on conda forge."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "Coarse-graining a box of 1000 methanol molecules from the all atom simulation given in `CG_tutorial/meth.gsd` file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# disable GPU. Remove this if you've compiled HOOMD for GPU\n",
    "import os\n",
    "os.environ['CUDA_VISIBLE_DEVICES'] = '-1'\n",
    "\n",
    "\n",
    "import tensorflow as tf\n",
    "import hoomd\n",
    "import hoomd.md\n",
    "import hoomd.htf as htf\n",
    "import numpy as np\n",
    "import gsd, gsd.hoomd\n",
    "import matplotlib.pyplot as plt,matplotlib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build the Model\n",
    "\n",
    "Here we prepare the computations that will be executed at each step during the simulation. We have access to the neighbor list, positions, types, box dimensions of the simulation.\n",
    "\n",
    "We work in four steps:\n",
    "\n",
    "1. Create 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 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": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HOOMD-blue 2.5.2 DOUBLE HPMC_MIXED TBB SSE SSE2 SSE3 \n",
      "Compiled: 04/30/2019\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"
     ]
    }
   ],
   "source": [
    "# set-up the system\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",
    "# disabled particle sorting!\n",
    "c.sorter.disable()\n",
    "set_rcut = 10.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finding molecules...99.90%\n"
     ]
    }
   ],
   "source": [
    "#get mapping from molecule index to particle index\n",
    "molecule_mapping_index = htf.find_molecules(system)\n",
    "# get number of atoms\n",
    "N = sum([len(m) for m in molecule_mapping_index])\n",
    "# get number of molecules\n",
    "M = len(molecule_mapping_index)\n",
    "# get number of atoms in a molecule=MN\n",
    "MN = len(molecule_mapping_index[0])\n",
    "\n",
    "## TRY CHANGING DIFFERENT BEADS DISTRIBUTION FOR TESTING\n",
    "\n",
    "# create one bead mapping -> \n",
    "# 1 x 6: [1, 1, 1, 1, 1, 1] that means \n",
    "# all atoms contribute to CG bead equally\n",
    "# massess are accounted for in sparse_mapping\n",
    "# molecule_mapping = np.ones([1, MN], dtype=np.int)\n",
    "\n",
    "# create 2 bead mapping -> \n",
    "# 2 x 6: \n",
    "# [0, 0, 1, 1, 1, 1]\n",
    "# [1, 1, 0, 0, 0, 0]\n",
    "# 4 atoms in first bead, 2 in second\n",
    "molecule_mapping = np.array([[1, 1, 1, 0, 0, 0], [0, 0, 0, 1, 1, 1]])\n",
    "\n",
    "bead_number = molecule_mapping.shape[0]\n",
    "\n",
    "#create a mass-weighted M x N mapping operator \n",
    "cg_mapping = htf.sparse_mapping([molecule_mapping for _ in molecule_mapping_index], \n",
    "                                molecule_mapping_index, system=system)\n",
    "assert cg_mapping.shape == (M * bead_number, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define model\n",
    "class MappingModel(htf.SimModel):\n",
    "    def setup(self, CG_NN, cg_mapping, rcut):\n",
    "        self.CG_NN = CG_NN\n",
    "        self.rcut = rcut\n",
    "        self.cg_mapping = cg_mapping\n",
    "        self.avg_cg_rdf = tf.keras.metrics.MeanTensor()\n",
    "        self.avg_aa_rdf = tf.keras.metrics.MeanTensor()\n",
    "    def compute(self, nlist, positions, box):\n",
    "        # calculate the center of mass of a CG bead\n",
    "        box_size = htf.box_size(box)\n",
    "        mapped_pos = htf.center_of_mass(positions[:,:3], self.cg_mapping, box_size)\n",
    "        # create the mapped neighbot list\n",
    "        mapped_nlist = htf.compute_nlist(mapped_pos, self.rcut, self.CG_NN, box_size, True)\n",
    "        # compute RDF for mapped and C-C in all-atom\n",
    "        cg_rdf = htf.compute_rdf(mapped_nlist, [0.1,self.rcut])\n",
    "        aa_rdf = htf.compute_rdf(nlist, [0.1,self.rcut], positions[:,3], type_i=3, type_j=3)\n",
    "        self.avg_cg_rdf.update_state(cg_rdf)\n",
    "        self.avg_aa_rdf.update_state(aa_rdf)\n",
    "        return\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# build model. output_forces = False because \n",
    "# this model doesn't comptue forces\n",
    "max_neighbor_est = 256\n",
    "model = MappingModel(\n",
    "    max_neighbor_est, \n",
    "    CG_NN=max_neighbor_est, \n",
    "    cg_mapping=cg_mapping, \n",
    "    output_forces=False,\n",
    "    rcut=set_rcut,\n",
    "    check_nlist=True)"
   ]
  },
  {
   "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": 8,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----\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",
      "notice(2): Force mode is FORCE_MODE.hoomd2tf \n",
      "notice(2): Starting TensorflowCompute \n",
      "notice(2): completed reallocate\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:02:24 | Step 8 / 250 | TPS 0.770384 | ETA 00:05:14\n",
      "Time 00:02:35 | Step 19 / 250 | TPS 1.03705 | ETA 00:03:42\n",
      "Time 00:02:45 | Step 30 / 250 | TPS 1.07384 | ETA 00:03:24\n",
      "Time 00:02:56 | Step 41 / 250 | TPS 1.0386 | ETA 00:03:21\n",
      "Time 00:03:07 | Step 52 / 250 | TPS 1.03282 | ETA 00:03:11\n",
      "Time 00:03:17 | Step 63 / 250 | TPS 1.04259 | ETA 00:02:59\n",
      "Time 00:03:28 | Step 74 / 250 | TPS 1.04234 | ETA 00:02:48\n",
      "Time 00:03:39 | Step 85 / 250 | TPS 1.00901 | ETA 00:02:43\n",
      "Time 00:03:49 | Step 96 / 250 | TPS 1.04988 | ETA 00:02:26\n",
      "Time 00:03:59 | Step 107 / 250 | TPS 1.06528 | ETA 00:02:14\n",
      "Time 00:04:10 | Step 118 / 250 | TPS 1.03596 | ETA 00:02:07\n",
      "Time 00:04:21 | Step 129 / 250 | TPS 1.03033 | ETA 00:01:57\n",
      "Time 00:04:31 | Step 140 / 250 | TPS 1.0352 | ETA 00:01:46\n",
      "Time 00:04:42 | Step 151 / 250 | TPS 1.03773 | ETA 00:01:35\n",
      "Time 00:04:52 | Step 162 / 250 | TPS 1.0362 | ETA 00:01:24\n",
      "Time 00:05:03 | Step 173 / 250 | TPS 1.02768 | ETA 00:01:14\n",
      "Time 00:05:14 | Step 184 / 250 | TPS 1.0412 | ETA 00:01:03\n",
      "Time 00:05:24 | Step 195 / 250 | TPS 1.03206 | ETA 00:00:53\n",
      "Time 00:05:35 | Step 206 / 250 | TPS 1.05436 | ETA 00:00:41\n",
      "Time 00:05:46 | Step 217 / 250 | TPS 1.02297 | ETA 00:00:32\n",
      "Time 00:05:56 | Step 228 / 250 | TPS 1.03952 | ETA 00:00:21\n",
      "Time 00:06:07 | Step 239 / 250 | TPS 1.04558 | ETA 00:00:10\n",
      "Time 00:06:17 | Step 250 / 250 | TPS 1.0911 | ETA 00:00:00\n",
      "Average TPS: 1.03\n",
      "---------\n",
      "-- Neighborlist stats:\n",
      "64 normal updates / 0 forced updates / 0 dangerous updates\n",
      "n_neigh_min: 235 / n_neigh_max: 561 / n_neigh_avg: 420.825\n",
      "shortest rebuild period: 3\n",
      "-- Cell list stats:\n",
      "Dimension: 3, 3, 3\n",
      "n_min    : 175 / n_max: 292 / n_avg: 222.222\n",
      "** run complete **\n"
     ]
    }
   ],
   "source": [
    "### Hoomd-Sim code ###\n",
    "\n",
    "tfcompute = htf.tfcompute(model)\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",
    "tfcompute.attach(nlist, r_cut=set_rcut)\n",
    "\n",
    "#Hoomd production run\n",
    "hoomd.run(250)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "Now we load the RDF we computed and plot it. Note this is not actually coarse-grained, \n",
    "we're just looking at the mapped rdf which can then be used for coarse-graining."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEQCAYAAABFtIg2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABBNklEQVR4nO3dd3hb5dn48e+jZclT3jt7OcNxJgFSEqBAgECYJbRQxsss0Akd9C39taV93xbeFspoy4YCYY9AAjTssDLJcvaw472XbMlaz+8PKSbDO7bkcX+uS5etc55zdB850a3nPEtprRFCCCHaYwh3AEIIIQYuSRJCCCE6JElCCCFEhyRJCCGE6JAkCSGEEB2SJCGEEKJDIUsSSimrUmqtUmqzUipfKfW7dspEKKVeVErtVUqtUUqNClV8QgghjhXKmkQrcJrWejqQByxSSs07qsx/AXVa63HA34A/hzA+IYQQRwlZktABjuBTc/Bx9Ei+JcDTwd9fAU5XSqkQhSiEEOIoIW2TUEoZlVKbgEpgldZ6zVFFMoEiAK21F2gAEkMZoxBCiG+YQvliWmsfkKeUsgOvK6Wmaq239fQ8SqkbgBsAoqKiZk2aNKlvAxVCiCFuw4YN1Vrr5K7KhTRJHKK1rldKfQQsAg5PEiVANlCslDIBcUBNO8c/AjwCMHv2bL1+/fr+D1oIIYYQpVRhd8qFsndTcrAGgVLKBpwB7Dyq2HLgquDvlwAfapmBUAghwiaUNYl04GmllJFAcnpJa/22Uur3wHqt9XLgceDfSqm9QC2wNITxCSGEOErIkoTWegswo53tdx32uwu4NFQxCSGE6FxY2iSEEEOXx+OhuLgYl8sV7lAEYLVaycrKwmw29+p4SRJCiD5VXFxMTEwMo0aNQoY5hZfWmpqaGoqLixk9enSvziFzNwkh+pTL5SIxMVESxACglCIxMfG4anWSJIQQfU4SxMBxvH8LSRJCiCHpjTfeQCnFzp3f9LQvKChg6tSpAHz88ccsXry4w+N//OMfk5mZid/vP+Kc27dv75d4HQ4HN954I2PHjmXWrFksXLiQNWuOnpQi4J133mH27NlMnjyZGTNm8LOf/axfYgJJEkKIIWrZsmXMnz+fZcuW9fhYv9/P66+/TnZ2Np988knb9v5MEtdddx0JCQns2bOHDRs28OSTT1JdXX1MuW3btnHrrbfy7LPPsn37dtavX8+4ceP6JSaQJCGEGIIcDgefffYZjz/+OC+88EKPj//444+ZMmUKN998c1uS+eKLL1i+fDl33HEHeXl57Nu3j02bNjFv3jxyc3O58MILqaurA2DhwoX85Cc/Yfbs2eTk5LBu3Touuugixo8fz3//938f83r79u1jzZo13H333RgMgY/l0aNHc+655x5T9i9/+Qu//vWvOTQdkdFo5Oabb+7xNXaX9G4SQvSbf75VzP4yZ5+ec0y6jZvOy+q0zJtvvsmiRYuYMGECiYmJbNiwgVmzZnX7NZYtW8bll1/OkiVLuPPOO/F4PJx00kmcf/75LF68mEsuuQSA3NxcHnjgARYsWMBdd93F7373O+677z4ALBYL69ev5/7772fJkiVs2LCBhIQExo4dy09+8hMSE7+ZuzQ/P5+8vDyMRmOXsW3btq1fby8dTWoSQoghZ9myZSxdGpiwYenSpT265eR2u1m5ciUXXHABsbGxnHDCCbz33nvHlGtoaKC+vp4FCxYAcNVVV/Hpp5+27T///PMBmDZtGlOmTCE9PZ2IiAjGjBlDUVHR8VxeSElNQgjRb7r6xt8famtr+fDDD9m6dStKKXw+H0op7rnnnm4d/95771FfX8+0adMAaGlpwWazddrI3Z6IiAgADAZD2++Hnnu93iPKTpkyhc2bN+Pz+Y6pTTz00EM8+uijAKxcuZIpU6awYcMGpk+f3qN4ektqEkKIIeWVV17hyiuvpLCwkIKCAoqKihg9ejSrV6/u1vHLli3jscceo6CggIKCAg4cOMCqVatoaWkhJiaGpqYmAOLi4oiPj28777///e+2WkVPjR07ltmzZ/Pb3/6WQ3OaFhQUsGLFCm655RY2bdrEpk2byMjI4I477uBPf/oTu3fvBgKN7P/85z979brdIUlCCDGkLFu2jAsvvPCIbRdffHG3bjm1tLTw7rvvHtFgHBUVxfz583nrrbdYunQp99xzDzNmzGDfvn08/fTT3HHHHeTm5rJp0ybuuuuuTs7euccee4yKigrGjRvH1KlTufrqq0lJSTmmXG5uLvfddx+XX345OTk5TJ06lf379/f6dbuiBvtM3LKehBADy44dO8jJyQl3GOIw7f1NlFIbtNazuzpWahJCCCE6JElCCCFEhyRJCCGE6JAkCSGEEB2SJCGEEKJDkiSEEEJ0SJKEEGLIUUpxxRVXtD33er0kJyf3eNR0Xzl8ivKjlZWVHRHX2rVrOeWUU5g4cSIzZszguuuuo6WlBYB3332XuXPnMmnSJPLy8rjssss4ePAgALfffjsffvhhn8cu03IIIYacqKgotm3bhtPpxGazsWrVKjIzM8MdVrv++te/cv311wNQUVHBpZdeygsvvMCJJ54IBEaQNzU1sX//fm677TaWL1/eNuZh+fLlFBQUMGLECG677Tauv/56TjvttD6NT2oSQogh6ZxzzmHFihXAN7O6HrJ27VpOPPFEZsyYwUknncSuXbsAeOqpp1iyZAkLFy5k/Pjx/O53vwMCNYFJkybxve99j5ycHC655JK2b/cbNmxgwYIFzJo1i7POOouysrK27dOnT2f69Ok89NBDHcb56quvsmjRIiAwT9NVV13VliAALrnkElJTU/nzn//MnXfeecSguPPPP59TTjkFgJEjR1JTU0N5eflxv3eHk5qEEKLfvFL0IsUtfTvjaVZkNpdkX9ZluaVLl/L73/+exYsXs2XLFq699tq2eZYmTZrE6tWrMZlMvP/++9x55528+uqrQCCBbNu2jcjISObMmcO5555LUlISu3bt4vHHH+fkk0/m2muv5eGHH+ZHP/oRt912G2+++SbJycm8+OKL/PrXv+aJJ57gmmuu4cEHH+SUU07hjjvuaDfGAwcOEB8f3zYB4LZt27jqqqvaLZufn8/tt9/e6TXPnDmTzz//nIsvvrjL96e7pCYhhBiScnNzKSgoYNmyZZxzzjlH7GtoaODSSy9l6tSp/OQnPyE/P79t3xlnnEFiYiI2m42LLrqIzz77DIDs7GxOPvlkAK644go+++wzdu3axbZt2zjjjDPIy8vj7rvvpri4mPr6eurr69u+5V955ZXtxlhWVkZycnKPr62mpoa8vDwmTJjAvffe27Y9JSWF0tLSHp+vM1KTEEL0m+584+9P559/Prfffjsff/wxNTU1bdt/85vfcOqpp/L6669TUFDAwoUL2/YppY44x6Hn7W3XWjNlyhS+/PLLI/bV19d3Kz6bzYbL5Wp7fmga8CVLlhxTdsqUKWzcuJHp06eTmJjIpk2buPfee3E4HG1lXC4XNputW6/dXVKTEEIMWddeey2//e1v29aGOKShoaGtIfupp546Yt+qVauora3F6XTyxhtvtNUeDh482JYMnn/+eebPn8/EiROpqqpq2+7xeMjPz8dut2O329tqIc8991y78U2YMIGCgoK257feeitPP/00a9asadv22muvUVFRwc9//nP++Mc/smPHjrZ9h9pFDtm9e3eHvah6S5KEEGLIysrK4oc//OEx23/+85/zq1/9ihkzZhyzANDcuXO5+OKLyc3N5eKLL2b27MBEqRMnTuShhx4iJyeHuro6br75ZiwWC6+88gq/+MUvmD59Onl5eXzxxRcAPPnkk9xyyy3k5eXR0WzbUVFRjB07lr179wKQmprKCy+8wO23387EiRPJycnhvffeIyYmhmnTpnH//ffz/e9/n4kTJ3LyySezY8cOvvvd7wKBBLV37962ePtKyKYKV0plA88AqYAGHtFa339UmYXAm8CB4KbXtNa/7+y8MlW4EAPLYJ4q/KmnnmL9+vU8+OCDR2wvKChg8eLFbNu2rc9f8/XXX2fDhg3cfffdx32ejRs38oc//OGYfcczVXgo2yS8wM+01huVUjHABqXUKq319qPKrdZah2fEixBChNiFF154RHtJb3m9Xn72s5/1QURHCtuiQ0qpN4EHtdarDtu2ELi9J0lCahJCDCyDuSYxVA26RYeUUqOAGcCadnafqJTarJR6Ryk1JbSRCSGEOFzIu8AqpaKBV4Efa60bj9q9ERiptXYopc4B3gDGt3OOG4AbAEaMGNG/AQshekxrfUyXUREex3u3KKQ1CaWUmUCCeE5r/drR+7XWjVprR/D3lYBZKZXUTrlHtNaztdazezMQRQjRf6xWKzU1Ncf94SSOn9aampoarFZrr88RspqECnyteBzYobX+awdl0oAKrbVWSs0lkMSOv0VHCBEyWVlZFBcXU1VVFe5QBIGknZWV1evjQ3m76WTgSmCrUmpTcNudwAgArfU/gUuAm5VSXsAJLNXydUSIQcVsNjN69OhwhyH6SMiShNb6M6DTm5Ra6weBBzsrI4QQInRkxLUQQogOSZIQQgjRIUkSQgghOiRJQgghRIckSQghhOiQJAkhhBAdkiQhhBCiQ5IkhBBCdEiShBBCiA5JkhBCCNEhSRJCCCE6JElCCCFEhyRJCCGE6JAkCSGEEB2SJCGEEKJDkiSEEEJ0SJKEEEKIDkmSEEII0SFJEkIIITokSUIIIUSHJEkMID6/Zu3ORuodnnCHIoQQAJjCHYD4xuufVfH4O6UYDJA3NoaF0+2cNMVOlNUY7tCEEMOU1CQGiLLaVp59v4xZ42O45FsplFS38tdXirj8j9u4+9kDVDW4wx2iEGIYkprEAKC15oHXizEYFD+6OJvkOAtXn5XOzoMtfLKljvfW1/I/zxfwlxvGYzKqcIcrhBhGpCYxAHy0qY6v9zZxzVnpJMdZAFBKkTMyipvOy+JHF2az42ALz75fHuZIhRDDjSSJMGto9vKvFSVMyo7knBOS2i2zMC+eM2cn8NInFWza1xTiCIUQw5kkiTB7dEUJzU4fP7ooG6Oh41tJN5+XSVZSBPe8WEi9wxvCCIUQw5kkiTD6em8TH3xdx6ULUhmVZuu0rNVi5JeXj6LJ6eP/Xi7E79chilIIMZyFLEkopbKVUh8ppbYrpfKVUj9qp4xSSv1dKbVXKbVFKTUzVPGFmsvt5++vF5GZGMHlp6Z265gx6TauPyeD9bubeOPzqn6OUAghQluT8AI/01pPBuYBtyilJh9V5mxgfPBxA/CPEMYXUp9vq6e81s0PlmRhMXf/z7B4XhInTo7jyffK2HrA0Y8RCiFECJOE1rpMa70x+HsTsAPIPKrYEuAZHfAVYFdKpYcqxlDaX+bEbFJMHxPdo+OUUvzk4mxS7RbufHwfK9dW91OEQggRpjYJpdQoYAaw5qhdmUDRYc+LOTaRoJS6QSm1Xim1vqpqcN52OVDuYmSKFWMvxj3ERJr42y3jyRsbzQOvF/P314twe/39EKUQYrgLeZJQSkUDrwI/1lo39uYcWutHtNaztdazk5OT+zbAECkodzIqzdrr42NsJv7fVWP4zoIU3llbwy8f3Utto8z5JIToWyFNEkopM4EE8ZzW+rV2ipQA2Yc9zwpuG1LqHV7qHF5Gd9GjqStGg+KaRRnc+d1R7C9zcduDu9hb2tJHUQohRGh7NyngcWCH1vqvHRRbDnw/2MtpHtCgtS4LVYyhUlDhBOiy22t3fWuanb/9YDxKKf7yYiEeufUkhOgjoaxJnAxcCZymlNoUfJyjlLpJKXVTsMxKYD+wF3gU+EEI4wuZA2WBJDH6OG43HW10mo3bLsiiqLKV16V7rBCij4Rsgj+t9WdAp620WmsN3BKaiMKnoMJFXJSJ+Bhzn573hJw4Tpwcy/MfVLBwejwpdkufnl8IMfzIiOswOFDm7NNaxOFuXJwFaP719pBryhFChIEkiRDz+TUHK1191h5xtNR4C5eflsYX+Q2s29WrzmNCCNFGkkSIldW00urR/VaTALhofjJZyRE8vLyYVo80Ygshek+SRIgVVLiAvuvZ1B6zycAPzs+ivNbNy59U9NvrCCGGPkkSIXagzIlBwYiU/qtJAMwYF8OCXDsvfVJJaXVrv76WEGLokiQRYgUVLtITI7Ba+v+tv/7cTExGxdP/GXJDTYQQISJJIsT6s2fT0RJjzZw9J5HP8+upa5IpO4QQPSdJIoRcbh/lde5+bY842llzEvH54f2NtSF7TSHE0CFJIoQKK1xo3bcjrbsyIsXK1FFRvLuuhsBYRSGE6D5JEiF0oLz/eza1Z9GcREpr3GzZL4sUCSF6RpJECBWUO7FaDKTFh3a6jPnT7ERbjbyztiakryuEGPwkSYTQgXIXI1OtGAw9X2joeESYDZw2I57P8xtoaPaG9LWFEIObJIkQ0VpTUO487jUkemvR3ES8Ps0H0oAthOgBSRIhUtvkpbHFd1yr0R2P0Wk2JmVHSgO2EKJHJEmESEH5oTUkwlOTADh7biJFVa3kFzSHLQYhxOAiSSJEDvVsCmX316OdkmsnMsIgDdhCiG6TJBEiBeVOEmPNxESGbJ2nY1gtRk7Ni2f1tnqaWqQBWwjRtfB9Yg0zB8pDNx1HZ86em8iKNTX8+/1yJmVH4nT7cQUfY9JtnDg5LtwhCiEGEEkSIeD1aQ5WtjJzfGy4Q2FsRiSTsiN568tq3vryyH0mo+IfP5pIVnL4k5kQYmCQJBECJdWteH39u9BQT/zhmjGU17mxWYxYLQasFgPOVh83/G0nj64s5XdXjQl3iEKIAULaJELgwADo2XS4aJuJcRmRZCZFkBhrJspqJCnOwndPS2PtzkbWy7KnQoggSRIhUFDuwmCAzOSIcIfSqfNPSiIj0cIjK0rw+mQshRBCkkRIFFQ4yUqyYjEN7LfbYjJw/bmZFFW18vZX1eEORwgxAPT4U0spFaWUMvZHMENVYbkrbCOte+qESbHMHB/Ds++XUe+QbrJCDHddJgmllEEp9V2l1AqlVCWwEyhTSm1XSt2jlBrX/2EOXi2thxYaGhxJQinFDedm4nT7+ff7suypEMNdd2oSHwFjgV8BaVrrbK11CjAf+Ar4s1Lqin6McVArrAiOtE4dGI3W3TEy1cp585J4d20N+8uc4Q5HCBFG3UkS3wb+CCzWWvsPbdRa12qtX9VaXwy82F8BDnYFwek4Rg6SmsQh3/t2GtE2I/94qxifXxqxhRiuukwSWmtPMDks7qxMV+dRSj2hlKpUSm3rYP9CpVSDUmpT8HFXV+ccDAorAgsNpdpDu9DQ8Yqxmbj27Ay2HWjm+Q/Kwx2OECJMejKYbotS6rfAHw6vUfTAU8CDwDOdlFmtte4wGQ1GBce50FB+w1ZWlr1NmjWdOQlzmRAzCYMKTS+pM2clsO2Ag+c/rGBCViQn5MiUHUIMNz1JEgnAAuBmpdQaYAuwRWv9cncO1lp/qpQa1fMQBy+tNQcqnJzYiw/XVl8rr5e8wuqqT0i0JFHuLOOrmi+INcUyK2EOsxLmMDJyVL8mDKUUt16QzYFyF/e8VMjfb5lIRtLAHushhOhb3U4SWuvvACilIoApwDTgBKBbSaKbTlRKbQZKgdu11vntFVJK3QDcADBixIg+fPm+Ve/w0tjsY1QPR1oXNB/g6QNPUNVayempZ3BexgVoNPkNW1lfu5bVVZ/wUeUHmJWZdFtG4GHNIMOWgd0ST5QxikhTJGZlQanjWyo1wmzgN1eM4rYHdvOHZw/wtx+Mx2qRHtBCDBeqq1XKlFJKd1GoO2WC5UYBb2utp7azLxbwa60dSqlzgPu11uO7Oufs2bP1+vXruyoWFl/vbeLOx/fxP9eNJW9sTLtl/NpPs9dBk9eBw9vE7qadvFf2DnHmOL4/+homxEw65hinr4Wt9VsoaimizFVCqbOUBk/9MeVMykSkMZJUaxrjYyYwPmYio6PGYDaYe3wtG3Y38pun9rMg187PLxt53MlHCBFeSqkNWuvZXZXrTk3iI6XUq8CbWuuDh72AhUA32KsIdJN9qpexAqC1bjzs95VKqYeVUkla60E79PfQanSjUo/t2bS3aQ9PHniUBk8DmiPz65yEE/hO9uVEmiLbPa/NGMncxHnMTZzXtq3Z20yZs5RGbwMt3hZafM20eFto9jZT5DzIO2UrWFn2NiZlYlTUaOYlnsS8xJO6/WE/a0IsV56RxjP/KWdidhQXnJzc3bdBCDGIdSdJLAKuBZYppcYAdYCNQM+o/wD3aa2/Pt5AlFJpQIXWWiul5gbPP6iXUDtQ7sIebcIefeQ3d6/fy/MH/41SBs5OP5coUzQxphiiTTHYLXZSrWk9fq0oUxTjYjqueLV4W9jn2Msex252NObzbOHT7HPs5bIR3+12zeKyBansKW7hsZUlzBofQ3bK4OrWK4TouS6ThNbaBTwMPKyUMgNJgFNrXd+TF1JKLQMWAklKqWLgt4A5+Br/BC4h0CjuBZzA0u7cwhrICitc7dYiPq76kApXOTePu5WpcbkhiSXSFMk0ey7T7Ln49UWsLHuLd8pWUOYq5foxN2O32Ls8h8Gg+OGFI/ive7fz2MpSfne1TCkuxFDX7YZrpdQeYCuwGdiklNqktS7s7vFa68u72P8ggS6yQ4LfrymscHH23MQjtjd46llZ+hZT46aFLEEczaAMLM5YQpYtm6cLnuTPO+7m+rE3MyZ6bJfH2qNNXH5aGo+/U8rGPU3MHN9+W4sQYmjoSf/JfwHlBG4BnQ3kK6W2KqV+H6xhiMOU17pp9fiPqUm8UfwaPu3j4qzLwhTZN/LiZ3LHpF9iMUZw3+57+bL6824dd/5JSaQlWHh0RQk+mVJciCGtJ0niCq31D7TWD2qtbyLQaP0R0Aj8tV+iG8QKKoKN1odNx7HPsZe1tV9xeuoZpFhTwhXaETJsmfxi0p2Mj57Ac4XPsLm+6+Yli8nAf52dQUGFi/fWD+pmIyFEF3qSJBqUUm33R7TWm4AFWut7gZP7OrDBrm3OpmBNwq/9vHRwGXZzPGelnRPO0I4RaYrixnG3MDJqFE/uf4zC5oIujzl5ShxTR0XxzKpyml2+/g9SCBEWPUkSNwJPKqUeV0rdppR6EGgJ7htcExOFQEGFi7QES9vAs8+qP6XYWcRFWZcSYRx4o5YtBgs3jr2FWHMs/9z7ILXuzmsISimuPzeThmYvL3xUEaIohRCh1u0kobXeCcwF3gVSgL3AYqVUFPBC/4Q3eBWUOxkdvNXk8Dp4u+RNJsRMZGb8rDBH1rFYcyw3j7sNj/bwjz0P4PR1Pk34hKxITp8RzxufV1FW2xqiKIUQodSjiX+01j6t9cta699ore/TWtdorZu11nf3V4CDkdvjp6SmlVHBNSTeLVuB0+fk0uzLB/xI5XRbBteNuZFyVzmP738En+78VtLVZ6VjNCieeKc0RBEKIUJpYC+6PEgVVbnw+wNrSNS561hd9QknJJ5Ihi0j3KF1y6TYySwd+T12NObz0sFldDZcJSnOwqULUvhsWwP3vlRIU4sseSrEUNKTWWBFNx1qtB6VauO98lfRaM5OPzfMUfXMyUnfospVxaqKd0mKSOKMtEUdlr1sYSo+n+bFTyrYuKeJWy7I4uQp9tAFK4ToN1KT6AcFFS5MRoU1pokvqj/jxMT5JEYkhTusHjs/8wJmxc/mjZLX2FC7rsNyJqPi+2emc/8tE4iPMXP3swX8z/MF1Du6XItKCDHASZLoBwXlTrKTI1hVuRKFYlH6wOry2l0GZeDKUdcwNnoczxQ8yV7Hnk7Lj8uI5P5bJnDVmel8sb2BG/+2k/fW1+CX5U+FGLQkSfSDggoXmdlO1tR8ybeSFxBviQ93SL1mNpi5YewPSLAk8sjeh6lwdb6UqcmoWHpqKg/eNpHsFCv3vVrET/+5h93FLZ0eJ4QYmCRJ9LEmp5fqBg/ejC8wKmOn9/IHi2hTNLeM/yFKGXh4z99p8jR2eczIVCv33DCOO74zgsp6Nz9+eDcPvF5EY7M0bAsxmEiS6GOFFS7MMbVUmrewMOU04sxDY13opIhkbhp3Cw2eBh7e+wCVrq4H0CmlOG1GAo/+NIcLT07m3fU1XPd/O3j4zWI27G7E7e3NUulCiFDqcmW6gW6grUz39lfVvFz+OAkjivlD7p+INg2tWVK31m/mqQOP49VeTk89g7PSzun2CPLCCif/fr+c9buaaPX4sVkMzJoQwwk5ccyfGifLogoRQt1dmU6SRB/7v5Xr2Z/6CIvSzuG8zAvCHU6/aPA08Gbxa6yp/RK7OZ6Lsi5lZvysbg8UbPX42byviTU7Glmzs5GaRg/2aBOXLUzlnLmJWMxSwRWiv0mSCJMfrboPb+xe7pn1ZyJNUeEOp1/tc+zlpYPLKHYWMSFmEleM/H6Pu/r6/Zr8wmaee7+czfsdJMWZ+e5paZwxKwGTcWCPThdiMOtukpCvbH3I4XXgidtFbHPukE8QAGOjx/GLnF9z2YjvcrC5gP/Z8Qc21vUsYRsMimmjo/nf68fxp/8aS2Ksmb+/XsSNf9vB1gOOfopcCNFdMuK6D60u/wJl9DHRPC/coYSMQRk4JXkhk2On8OSBx3h8/yPsSNrOJVmX9Xi22xnjYsgbG82anY08tqKU3z69n3tvHM+YdFs/RS/6g9aaZpcPh9NHq0fT6vHj9vhp9fhJjDUzKk3+noOJJIk+orXm8+rVuGpSyckeGe5wQi4pIpmfTryDt0vfYlX5u+xz7OXa0deTFZndo/MopZiXE8fYDBs//cce7npqP3+9eTwpdpmNvi9orSmtcWM2KWIjTVgt39xM8Pk0xdWt7CttYX+Zk4OVrWQkWsgbF8O00dFEWY/sWNDk9LKjsIWdB5sprWmlpslDTYOH2iYPrZ6Ob2OPTrPy7ZkJnJoXT3zMsYtaNjm9eLyahHb2idCTNok+stexh7/tuofKdady/6UXkZE48NaMCJWdjTt4puAJHF4HU+NymRU/m6lxuT2uWRSUO7n9X3tIjLVw703jiLHJd5rjUdfk4YE3ivlye0PbtgizgbgoI7YII2U1rbi9gc8Dk1GRmRTRts1ggIlZkeSOiaGxxcv2wmYKKwJzlBkMkGq3kBhnJjHG3PYzJtKIxWwgIviwmBQHyly8v7GWXcUtGAwwa3wsk0dGUVbbSnFVKyXVrTQ0e1EKFuTaufy0NEakWNu9HnF8pOE6xJ468DgbqjdRtOIqXr1rFkbD8G50bfI08V75SjbWrafB04BZmZlqDySMXHseRtW97q5b9jfx6yf2M2lEJH+8Zqz0fOql1VvrefDNIpytfpYuTCUx1kxDs5eGFi8NDi/NLh8ZiRGMSbcxJsNGdrIVk1Hh9vjZcbCZTfscfL23iT3FLdgiDOSMiGLyyMBjYnZkj7svH6x08cHGWj74uo6aRg/x0SYykyPISrKSlRxBvcPL219V0+rxsyA3nu+dnkpWshWtNSXVrWwraGbbAQcFFS7SEyyMzYhkbIaNsek2EmKlBtIdkiRCqNnbzJ1b7sBYNY3WXafy0A8nhTWegcSv/exz7GVj3Xq+rttAk7eJDGsGl2QvZWJs996njzfX8ecXCjllmp1fLB2JYZgn4J5oavHy8PJiPt5cz/hMGz+7dGTbkrq94XL7sZhUn/0NfP5Am0VkxLFJpt7h5dXVlbz1ZTUer5+po6MpqnRR5wiM2rdHmxidZqOstpXyWnfbcfExJhbmxnPuvCQyk4Zvjb4rkiRC6KOKD3il+EWca77HhPiR/PLyUWGNZ6DyaR9b6jfxevEr1LhryLPP5KKsS7rVbfaVTyt5/J1SzpuXxE3nZQ7rRLGnpIWSqlZmTYzp8BZcU4uXD76u4+VPKmho9vLd09L4zsLUQdmtuN7h4dXVVazd2cjYDBtTR0eROzqazKSItrE5zS4f+8uc7Ct1sq3AwVfbG/D5Yeb4GBbPS2LuxFiM3bj2Vo8fj9eP2WTAbOy7ZDgQSZIIEa01d2//f1hUBJ88ew5XnJ7Gd09PC1s8g4Hb7+aDilX8p/wdtNZ8O+0szkw7C4uh4299WmueeLeMVz6tZNGcBG69IHvY3dLbVdTMcx9UsG5XYO4sk1Exe0IMC3LjmTc5lgizgfzCZt5dW8PqrfW4vZqcEZH8YEkW4zIiwxx9aNU2eXhvXQ0r19ZQ3eAhKc7MmHQbcVEmYiONxEaZiLGZaGzxUlrTSmlNK2U1bmoaj5ze3mAAs9FAbKSRZLuFFLs58DPOwqh0K5NHRA3aRNLdJCEtgcdpf/M+yl1lnBG7lI81ZEsjW5csBgtnp5/LvMQTeb34Vd4pe5u1NV+xdMT3mBw3pd1jlFJcuygds0mx7MMKXG7N7ZeO6Na3w8FuR2Ezz31QzoY9TcTYjFx1ZjrTRkfxRX4Dn2yp56sdjUSYDSTGmiitcRMZYeDM2YksmpPA2BAmB7/20+x1UO+pp95dj8PrwOVz0uJraftpUiZSrKmkRKSSYk0lKSIRo/rmY0hrjVd7UShMht5/PCXEmLn8tDS+syCVNTsbeH9jLZX1Hg6UOWlo9rY10EPg9lRGYgQzx8eQnmDBajHg8Wo8Ph386aex2UtVvYddRS18tq0Bry9wfGKsmflT41iQG8+kEZEDfnni3pCaxHF65sCTbK7/mrP0nfzt5XL++eOJjEyVfuA9sbtpFy8UPktFawWz4+dycfZ3iDXHdlj+xY8reOq9Mk6eEscvlo7EbBqajdml1a38a0UJa3c2Ehtl5OJvpbB4XtIR9+8PjVj/ZHMdpTVuFuTaWTDd3mfzYHn8Hho8DdS5a6n31FHvrqPZ23zYB78Tp6+FJk8jDZ4GvLr9WX4thggijTbcfjctvm+mjTdgIMoUjVd78fo9eHTgm7xCkWBJJN2WTqo1nTRrGhm2TLIjR3S700NnXG4/TS1eom2Bnl094fdr6h1etux38OnWetbvbsTj1aTYzcyeEEtagoXUeAspdgsp8Rbio00DMnnI7aYQaPE2c+eWnzMv6SRad57GS59U8Mbvcofsh1Z/8vg9/Kf8Xf5T/g5mg4ULMy/mxKSTMaj238s3Pq/iX2+XMHtCDP99xWgihlCvJ5fbz0sfV/Dyp5WYg+tznHdiUo8/zNpT6apgZ+MO4i0JZNgySbAkHPEBVu+uY1fTTnY27mBP0y7qPHXHnMOkTNiMNmzGyOBPGzHmWOxmO3EWO3azHbs5nhhzDDZjJFaj9YgPdofXQZWrksrWCipdFTR5mzAbzJiVGVPwp0d7qHSVU+4qp8JV3pZ8bEYbE2NymBw7hZy4ySRYEvFrP42eBqpaK6lqraLOXYvNaCPWHEecOS74047V2D+1/GaXj6+2N/Dplnq2FzbjcPmO2B9hVmQnWxmRYmVkauDnqDQrqfGWTpOH36+pc3hJiOmfJDPgkoRS6glgMVCptZ7azn4F3A+cA7QAV2utN3Z13nAmiY8qP+CVohf5Zc5veOo1LwcrXTz605ywxDJUVLjKWVb4LHscuxkXPZ4rRl1FckRKu2XfWVvDA28UMSrVypVnpDMvJ3ZAfmPrLq01X2xv4JG3S6is97BwejzXnZNB4nF26XT7W/m6biNfVn/OHsfuI/bZjDbSrRkkRSRzsKWQclcZEFhDZELMJNKtGdgtduItCcSb47Fb4vvtw7Yjfu2nurWaopaD7GzczvbGfOqDyctujqfZ62irgXQmKSKZsdHj2h6pEWn98u+l2eWjst5NZZ2byno3pTVuDla6OFjporrhmziT4sxMHxNN7phopo+NJsVuobTGzeZ9TWze52DzfgcNzV4mj4zi2kXpTBkV3adxDsQkcQrgAJ7pIEmcA9xGIEmcANyvtT6hq/OGK0kcarCOMETw85w7ueGvO8hKtnLXlaNDHstQo7Xmy5rPea34Zbx+L0syL2JByqnt1iq+yK/nsZWllNW6GZNu43unp3Li5Lg++8+vtcavwevT1Ds81DR6qW3yUNPoobHZS3aKlSkjo0ju5Yhwn1+zp7iFDbubWLOzgT0lTkalWfnB+VlMG927DwWtNTXuag62HGRX40421K3F6XOSFJHMSYnzmRE/C4e3iRJnMaXOEkqcJVS3VpJhy2JSTA4TY3PItGV2WIsLN6015a4ydjTmU9hSSJw5juSIFJIikkmOSCHeEo/L56LR00CDp4FGTwN17joKWw6wz7EXhzcwJ1i0KZosWzbptozAwxr4aTPaAn93fPi0H5/2YcCA2WA+7vek2eXjYKWLfSVOthxwsCWYCAAiIwy0tAbWWEmMDSSQjKQIVq6pprbJywk5sVx9ZnrbtCY+v2ZvSQtf73Ww9YADl9uPyagwGhRmo8JkUoxNt3HqjHjSE47tFDLgkgSAUmoU8HYHSeJfwMda62XB57uAhVrrss7OGa4ksbdpD3/bfQ/fG/l95safzAV3beaSU1K4+qyMkMcyVNW761h28Fm2NWxlbPQ4rhh5FSnW1GPK+XyaDzfV8cJH5ZTWuBmTbuWqMzOYO6njdo32eLx+nvugnDe/qMbt9aM1dPe/R3KcmZyRUUwZGcXsCbFkdNI/3+H0smZHI+t2NbJxTxNNTh9KwYSsSE7Li+fcE5K63SDv9Dkpd5VR7iynzFVCUctBilqKcAbv+5uVmbz4mZyUNJ9x0eMH7Ad/qGitqWytYJ9jL/sceylxFlPuLDuiJqJQaNr/wxswYDFYMBnMRATbWWymyLZbb9GmaBIjkki0JJEUkUSCJRGzoeOaoN+vOVjpYvN+BwXlLsZl2pg+5sjuvS63nze/qOLlTypoafVzyjQ7bq9my/4mml2BpDI6zYo92oTHp/H5NF6fptWjKapyoTVMHhnFqXnxnDLNTrTNSGW9m/RE66BLEm8D/6u1/iz4/APgF1rrTjNAuJLEk/sfI79xK3/K/Qvl1Zob/7aTO74zgtNmJIQ8lqFMa83a2q94uehFvH4P85NPYWpcLuOixx/T+8Xn03y0uY5lHwaSxVmzE7hxcWa37uUfKHNy78uF7C9zcco0O+mJERgMYFCBnlVGgyI+2kRCrJnEWDMJMWairAYKKlxsL2hme2Ez+YXNbV0oR6ZaOWlyHCdOiWNchi1w33pHI6u31rNxTxNenyY+xsSs8THMmhDLzHExxEYd25vH6/dS1VpJvaeOOncd9e566j11VLdWUe4qp8FT31bWpExk2rLIjhzR9siwZXb6ISUCt7Nq3NWUOUspd5XR6nNjUAaMyoBBGTEoA1prPMGGdY/fjcfvpdXvwhlsuHf6nDi9LTR5m45ovFco7JZ4Mm1ZZNmyyYwM/EyKSMKrPcHjnTi9TjzaTaw5jnhzQrtT2DQ2e3npkwre+rKa+BhzYELMcdFMHxODPbr9nmBV9W4+2lzHh1/XUVjhwmRUmE0KZ6ufd/93xtBNEkqpG4AbAEaMGDGrsLCwX+M+WpOnkV9v/QXfSl7ApdlL+WxbPX98roC/3zqB8ZnDqz96qNS763mt+CU212/Cq71EGCKYFJvDlNhpTLVPI85sbyt7qEbw0ieVpMVbuOM7I8kZ2f7U7T6/5tVPK/n3++XE2Iz88KJs5uX0bslZrTXldW7W7Gjki/x68gua8WvapsDw+gI9YOZPtTN/mp2JWZEd9rFv8TazuupTPq78gEbvkWuKx5hiSbAkkGZNJ82WFuz9k05SRFKf9PwRvXeoEb3aXU1NazXVrdVUuioodhZR6arAT+Cbf2e1FYAoYxQJlgTslgQijZFYjBYiDBFEGKyYlBmlwKe9+LQPr/bh0z6sBivR5miiTdFEm2KIMkXR7G2m3FlGmauUgoYSypxleA1OFIoH5/xj0I2TKAEOnzI0K7jtGFrrR4BHIFCT6P/QjvRVzRf4tI/5SacAUFQZmOgsO1mmAOgvdouda8fcQKuvld1Nu8hv3Mq2hq1srt+EOqgYFz2emQmzmWGfSYw5lqvPymD2hFjueamQ2/+1h6WnpnL5aWn4taa63kNFsGHxPxtq2V7YzPypcdx6QTZx7Xyb7y6lFOkJESyaF8O8mV6KGp1sKi5hX3UV08zJnDlmJnkjUzptL6lureajyvf5ovpz3P5WcmInc2HCPBIiEok3xxNnth/X+AHRvwzKgN0SaOAfFz3+iH1uv5syZyklzmJqWquJMEZgDfYOsxkjMStzW3fjWnctde4aatzVlPpctPpbafW1tttAb8CAURk7bbyPMESQZk1nZso0YkwxgSTBP7p1TQPpX9ty4Fal1AsEGq4bumqPCAe/9vNZ1aeMix5Pui3Q/nCwspUUu1nWaA6BCGME0+y5TLPnorWmzFXK13Ub2Vi3nhcPPs/LB19gQsxEptmnMzJlFPffNobH3q7k+Q8reOPzKpxu/xHtDDE2Iz+/bCQLp9t73Njt9Xspc5VR4iyipKWY4mBDsMPbdGTBJGgCHq95g1Gu0UyJm8qU2GlEGCOodFVQ2VpJpauCClc5ex17UCjmJJzA6alnkBmZdfxvmhgQLAYLI6NGMTJqVK/P4dd+Wv2tbYnBqIxt/2592ovD24zD24TD48DhdRBpiiTNmo7d3PN/34eELEkopZYBC4EkpVQx8FvADKC1/iewkkDPpr0EusBeE6rYemJn4w6q3dVHrF99sNIl0xmHgVKKDFsmGbZMzklfTKmrhI2169lQt56Xi14AAt+y0qdm8O3x6Tiqksk2TSI7LoXUeAup8WaSYi3dbiT2a3+wG+YOdjbtYL9jb9v9Z7Myk2HLJDduOsnW5ECXUUsCCZYEYk1xlDiL2dawlfyGrawofYu3S5cfce4oYxQp1lS+nXoWC1JOJd4S37dvlhgSDMqAzdj+YF2jMhEXHBtCH47nDVmS0Fpf3sV+DdwSonB67bPqT4g2xTDdPgMI9E4oqXYxfUzP1nYWfUspRaYti8zMLBZnLKHeU0dhcyFFLYUUthRS5NmBI3EdlaykxJrOFMtULMap2BmD0+umxXtYA6SvhRZfC06vkxZfM06fk3p3HXsce9p6DWXasjgl+VRGRo0kKzKblIjUTnsOHfoGeW7GeTR5GtnRuB2NJiUilWRrCtGmvu0DL0RfGUi3mwa8OncdW+u3cHrqmW09Rirr3bR6NCOOY/pl0beUUm3f5PPiA8n8UN/6/MZtbG/YxseVH/JBxaouz2XA0Na1Mc8+g0mxOUyMmURMJ9OGdCXGHMvcxOGzxK0Y3CRJ9MAX1Z+h0cxP/lbbtoOVrYA0Wg90Sqm2QVPfTj0Tl8/F7qadlDpL2hoPI43f9HePDPZ9jzBEDOpR3EIcL0kS3eTTPr6oXk1O7GSSIpLbth8M9mySNonBxWq0kmvPI9eeF+5QhBjQhvfwyx7Y1rCFek8985MXHLG9qMpFfLSJmEjJt0KIoUeSRDd9XrUau9nO1LhpR2wvqnTJGhJCiCFLkkQ3tHib2dG4ndkJc48Y0ap1YN4VaY8QQgxVkiS6YWvDFvz4yYufecT2uiYvzS6/1CSEEEOWJIlu+LpuI3ZzPCMjRx2x/WCVNFoLIYY2SRJdcPlc7GjMJy9+xjGDpYqkZ5MQYoiTJNGF/IateLWXPPvMY/YdrGwlMsJAQoz0bBJCDE2SJLrwdf1GYkwxjI0ed8y+A+VORqXZZLCVEGLIkiTRCbe/lfyGbUy3zzzmVpPfr9lf6mRsRh/OpCWEEAOMJIlObG/Yjtvfyozg/D+HK69143T7GZsuSUIIMXRJkujEpvqNRBmjGB8z4Zh9+8qcAFKTEEIMaZIkOuDxe9hav5lp9ukY1bEN0/tKnRgNyOyvQoghTZJEB3Y17cTldzEj/theTQD7SlsYkWLFYpK3UAgxdMknXAc21W3EarAyMSan3f37yqTRWggx9EmSaIdP+9hSv4mp9ty2xYUOV9vkoa7JyxhptBZCDHGSJNqxp2k3zb5mZrQzgA5gf+mhRuvIUIYlhBAhJ0miHZvqN2IxWJgcN6Xd/dKzSQgxXEiSOIpP+9hU9zVTYqdiMbQ/Bfj+Uidp8RairMZ29wshxFAhSeIou5t20eRtZHbC3A7L7JOR1kKIYUKSxFHW167FarAy5agV6A5pafVRWtsqjdZCiGFBksRhPH4Pm+o2khc/s91eTQAHypxoLe0RQojhQZLEYfIbtuLyu5idMKfDMvtKpdFaCDF8SJI4zPratcSYYpgQM6nDMvvLnMRGGUmMbb+mIYQQQ4kkiSCnz8nWhi3MjJ+NUXXca2lfqZOx6ZGyhoQQYlgIaZJQSi1SSu1SSu1VSv2ynf1XK6WqlFKbgo/rQhXblvpNeLW3015NXp+moMIlt5qEEMNGyNbdVEoZgYeAM4BiYJ1SarnWevtRRV/UWt8aqrgOWVe7lkRLIqOjxnRYpqjShdenZQ0JIcSwEcqaxFxgr9Z6v9baDbwALAnh63eoydPIrsYdzEqY2+ltpL3SaC2EGGZCmSQygaLDnhcHtx3tYqXUFqXUK0qp7FAEtrFuA378zOnkVhMEGq0jzAYyktofiS2EEEPNQGu4fgsYpbXOBVYBT7dXSCl1g1JqvVJqfVVV1XG/6PratWRYM8iwtZezvrGvtIXRaVaMBmm0FkIMD6FMEiXA4TWDrOC2NlrrGq11a/DpY8Cs9k6ktX5Eaz1baz07OTn5uIKqaa1hf/M+Ziec0Gk5rTX7ZQ0JIcQwE8oksQ4Yr5QarZSyAEuB5YcXUEqlH/b0fGBHfwe1oW4dALM6GUAHUFHnptnlZ4wkCSHEMBKy3k1aa69S6lbgPcAIPKG1zldK/R5Yr7VeDvxQKXU+4AVqgav7OSbW1XzF6KgxJEUkdVq2baR1uqwhIYQYPkKWJAC01iuBlUdtu+uw338F/CpU8eQ3bqPUVcr3Rn6/y7L7ypwYDDAqzRqCyIQQYmAYaA3XIaO1ZkXpchItSZyQOK/L8ntLnGQnW4kwD9u3TAgxDA3bT7xtDVs52FLIovRzMKrOK1Rur5+tBxxMGRUVouiEEGJgGJZJQmvNyrK3SOpmLSK/oBmX28/cibEhiE4IIQaOYZkktjVsCdYizu2yFgGwbmcjZpNi+tiYEEQnhBADx7BLElprVpS9RVJEMnO7UYsAWLurkeljorFaht3bJYQY5obdp97Whs0UtRxkUdq5nU4JfkhpdSsl1a3MlltNQohhaFglCa01K0vfIjkihbmJnY+wPmTd7kYA5kiSEEIMQ8MqSWxp2EyRs4iz07tXiwBYv6uRrOQIMhJlUj8hxPAzbJKEX/tZWfoWKREpnS4sdDiX28fm/Q7mTJBahBBieBo2SWJT/dcUO4s4O31xt2sRm/c58Hg1cyZJkhBCDE/DIkn4tI+3S98k3ZrR7VoEwLpdjVgtBhlEJ4QYtoZFklhbs4YKVznnZS7BoLp3yVpr1u1qZMa4aCymYfE2CSHEMYb8p5/H72FF2XJGRo4iNy6v28cdrHRRWe9hzsS4/gtOCCEGuCGfJD6vXk2du5bzMi/odP3qo63b1QTA7IkyyloIMXwN6STR6mvl3bIVjI+ewKSYnB4du25XI6PTrCTHWfopOiGEGPiGdJL4uPJDmrxNnJ95YY9qEc0uH/kFDhlAJ4QY9oZskmjxNrOq4j2mxuUyJnpsj47duKcJnx/p+iqEGPaGbJJ4v+I/OH0tnJexpMfHrtvVSLTVSE62dH0VQgxvQzJJVLoq+KjyA2bFzyErMrtHx9Y7vHy+rZ7ZE2MxGrt/i0oIIYaiIZckmjxNPLT371gMFpZkXtjj45//sByXx8/lp6X2Q3RCCDG4DKkk4fa38s99D9LgruemsbeSGJHUo+OLq1ysXFPNojmJjEix9lOUQggxeAyZJOHXfp468ASFzQVcPfo6RkeP6fE5nni3DIvZwJXfTuuHCIUQYvAZEklCa82rxS+xuf5rLs7+DnnxM3p8jq0HHHy5vYHvLEjFHm3uhyiFEGLwGRJJ4qPK9/m48kNOS/k2p6ac3uPj/X7NoytKSIozc8HJyf0QoRBCDE6DPknUumt5rfgV8uwzuTDrkl6d45Mt9ewpcXL1memyjrUQQhzGFO4AjpfD28SClFNZknlRt2d4PZzb4+ep90oZm2Hj1Lz4fohQCCEGr0H/tTndmsGl2UuxGHo3x9KbX1RRWe/hunMyMBhkXIQQQhxu0NckzIbeNTI7nF6e+6CCt76sYu6kWPLGymyvQghxtJAmCaXUIuB+wAg8prX+36P2RwDPALOAGuAyrXVBX8bg82veW1fD06vKaGrxsWhOItecld6XLyGEEENGyJKEUsoIPAScARQD65RSy7XW2w8r9l9AndZ6nFJqKfBn4LLjfW2X20dtk5eiShfPrCpjf5mLqaOjuGlxJmMzIo/39EIIMWSFsiYxF9irtd4PoJR6AVgCHJ4klgD/L/j7K8CDSimltdYdnbS4qpXb/7XnmO1+v6ax2Uetw4Oz1d+2PTnOzK8uH8m3ptl7NH24EEIMR6FMEplA0WHPi4ETOiqjtfYqpRqARKD68EJKqRuAG4JPHf9304RdPQnk37/qSek2SUfHMYgNpWuBoXU9Q+laYGhdz1C6FoCJ3Sk0KBuutdaPAI+E8jWVUuu11rND+Zr9ZShdCwyt6xlK1wJD63qG0rVA4Hq6Uy6UXWBLgMPn7c4Kbmu3jFLKBMQRaMAWQggRBqFMEuuA8Uqp0UopC7AUWH5UmeXAVcHfLwE+7Kw9QgghRP8K2e2mYBvDrcB7BLrAPqG1zldK/R5Yr7VeDjwO/FsptReoJZBIBoqQ3t7qZ0PpWmBoXc9QuhYYWtczlK4Funk9Sr6oCyGE6Mign5ZDCCFE/5EkIYQQokOSJLqglLIqpdYqpTYrpfKVUr8Ld0zHSyllVEp9rZR6O9yxHC+lVIFSaqtSalN3u/QNVEopu1LqFaXUTqXUDqXUieGOqbeUUhODf5NDj0al1I/DHVdvKaV+Evz/v00ptUwpNWjXN1ZK/Sh4Hfnd+ZtIm0QXVGBYdpTW2qGUMgOfAT/SWn8V5tB6TSn1U2A2EKu1XhzueI6HUqoAmK21HvSDnJRSTwOrtdaPBXsARmqt68Mc1nELTslTApygtS4Mdzw9pZTKJPD/frLW2qmUeglYqbV+KryR9ZxSairwAoEZMNzAu8BNWuu9HR0jNYku6ABH8Kk5+Bi0mVUplQWcCzwW7ljEN5RSccApBHr4obV2D4UEEXQ6sG8wJojDmABbcPxWJFAa5nh6KwdYo7Vu0Vp7gU+Aizo7QJJENwRvz2wCKoFVWus1YQ7peNwH/Bzwd1FusNDAf5RSG4LTtQxWo4Eq4MngrcDHlFJR4Q6qjywFloU7iN7SWpcA9wIHgTKgQWv9n/BG1WvbgG8ppRKVUpHAORw5yPkYkiS6QWvt01rnERglPjdYZRt0lFKLgUqt9YZwx9KH5mutZwJnA7copU4Jd0C9ZAJmAv/QWs8AmoFfhjek4xe8bXY+8HK4Y+ktpVQ8gclHRwMZQJRS6orwRtU7WusdBGbX/g+BW02bAF9nx0iS6IFg9f8jYFGYQ+mtk4Hzg/fxXwBOU0o9G96Qjk/wWx5a60rgdQL3WgejYqD4sFrqKwSSxmB3NrBRa10R7kCOw7eBA1rrKq21B3gNOCnMMfWa1vpxrfUsrfUpQB2wu7PykiS6oJRKVkrZg7/bCKyHsTOsQfWS1vpXWussrfUoArcAPtRaD8pvRABKqSilVMyh34EzCVSnBx2tdTlQpJQ6NDPn6Rw5jf5gdTmD+FZT0EFgnlIqMtiR5XRgR5hj6jWlVErw5wgC7RHPd1Z+UM4CG2LpwNPBHhoG4CWt9aDvOjpEpAKvB9cFMQHPa63fDW9Ix+U24LngLZr9wDVhjue4BBP3GcCN4Y7leGit1yilXgE2Al7gawb3FB2vKqUSAQ9wS1cdJKQLrBBCiA7J7SYhhBAdkiQhhBCiQ5IkhBBCdEiShBBCiA5JkhBCCNEhSRJCCCE6JElCiBBRSj2glNqolJoT7liE6C5JEkKEQHBgWQqBgWWDenp2MbxIkhCiDymlRimlnMFZg9torZsJjN7/GPh7sKwtuCCPWymVFPJghegGSRJCHAcVcPT/o33BWYMPL5dIYB2CJgJTO6C1dgbLDda1CcQwIElCiB4K1hZ2KaWeITChYKfz8Qf9N4E1CfKBKf0ZnxB9SZKEEL0zHnhYaz2lqxXXlFKjCEwt/SKB2UMlSYhBQ5KEEL1T2IN1zu8Gfq8Ds2lKkhCDikwVLkTvNHenkFIqj8Cc/fOVUg8BVmBrP8YlRJ+SJCFE//ozcL7W+n0ApVQqgfUIhBgU5HaTEP1EKXUaEHkoQQAEl/GMVkolhC8yIbpPahJC9JDWugCY2o1yHwIftrM9th/CEqJfSE1CiL7lA+KOHkzXnkOD6QAz4O/nuIToFVm+VAghRIekJiGEEKJDkiSEEEJ0SJKEEEKIDkmSEEII0SFJEkIIITokSUIIIUSHJEkIIYTokCQJIYQQHfr/u5kwqDuvtZEAAAAASUVORK5CYII=\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",
    "matplotlib.style.use('seaborn-muted')\n",
    "aa_rdf = model.avg_aa_rdf.result().numpy()\n",
    "cg_rdf = model.avg_cg_rdf.result().numpy()\n",
    "r = aa_rdf[1,:]\n",
    "plt.plot(aa_rdf[1,:], aa_rdf[0,:] / aa_rdf[0,-1],label ='All Atom C-C')\n",
    "plt.plot(cg_rdf[1,:], cg_rdf[0,:] / cg_rdf[0,-1], label='Mapped (CG)')\n",
    "plt.xlim(2.5,9)\n",
    "plt.ylim(0,3)\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.8.5"

  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}