{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 05. Using a Trajectory instead of HOOMD simulation\n", "\n", "In this notebook we show how to use an existing trajectory, possible from another simulation engine, to do computations. Reading is done via `MDAnalysis`, which is required to run these examples." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import hoomd\n", "import hoomd.htf as htf\n", "import tensorflow as tf\n", "import MDAnalysis as mda\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "# disable GPU\n", "import os\n", "os.environ['CUDA_VISIBLE_DEVICES'] = '-1'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Building our Model\n", "\n", "We demonstrate two types of computation. The first is to compute a potential energy function from the trajectory. This could then be trained by force matching (not shown). The second is that we compute RDFs from the trajectory." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "class TrajModel(htf.SimModel):\n", " def setup(self):\n", " self.avg_chrdf = tf.keras.metrics.MeanTensor()\n", " self.avg_ohrdf = tf.keras.metrics.MeanTensor()\n", " def compute(self, nlist, positions): \n", " # pairwise energy. Double count -> divide by 2\n", " inv_r6 = htf.nlist_rinv(nlist)**6\n", " p_energy = 4.0 / 2.0 * (inv_r6 * inv_r6 - inv_r6)\n", " # sum over pairwise energy\n", " energy = tf.reduce_sum(p_energy, axis=1)\n", " # get forces\n", " forces = htf.compute_nlist_forces(nlist, energy)\n", " \n", " # now get RDF\n", " # For reference, type indices in this case are: {C:0, H:1, N:2, O:3} \n", " # compute C-H RDF\n", " chrdf = htf.compute_rdf(\n", " nlist, [0, 15], positions[:, 3], \n", " nbins=20, type_i=0, type_j=1)\n", " # compute O-H RDF\n", " ohrdf = htf.compute_rdf(\n", " nlist,[0, 15], positions[:, 3], \n", " nbins=20, type_i=3, type_j=1)\n", " # average the RDFs\n", " self.avg_chrdf.update_state(chrdf)\n", " self.avg_ohrdf.update_state(ohrdf)\n", " return forces\n", "model = TrajModel(256)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run from trajectory\n", "\n", "We load the trajectory with `MDAnalysis` and then call the `iter_from_trajectory` command. It runs over the trajectory computing the graph and constructs neighborlists according to the r_cut. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "universe = mda.Universe('test_topol.pdb', 'test_traj.trr')\n", "for inputs, ts in htf.iter_from_trajectory(256, universe, r_cut=25):\n", " result = model(inputs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Analysis\n", "\n", "Now we plot our RDFs. Note that the trajectory is very short to save on space in the repo, so the RDFs are under-sampled. Also, HOOMD-TF makes no attempt at normalizing the RDF" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAApBUlEQVR4nO3deZhU1Z3/8feXpqHZGtlkFcE9ijsxRrNp1KCTiFsixkScMeOYn46JJpkxcUaNmSdjMtkTJ5GoMYlxiLskokaN2TQqrQKCuKAiNoLS3YDQxdLL9/fHudUURfVe996q5vN6nnq66ta9db8o1KfPOfeca+6OiIhIvn5pFyAiIqVJASEiIgUpIEREpCAFhIiIFKSAEBGRghQQIiJSkAJCREQKUkCIRMzs02ZWY2abzGy1mT1gZh9oZ98pZuZm1j9v+y1m9l/tHPMRM2uNPn+jmb1kZv+Yt4+bWWO0T72ZPWpmZ+ft8ycz2xLtk328v7d/fpF8CggRwMwuB34AfBMYC0wG/heYWeRTveXuQ4Fq4DLg52a2f94+h0b77A/cAvzEzK7O2+cSdx+a8/h7kesUoX/nu4j0bWY2HLgW+Ed3vzvnrd9Fj6LzsITBfDNrAA4BXiqwTx3wazPbDNxqZj9x9/o46hEpRC0IEXg/UAXck9QJzayfmZ0KjAaWd7L7fYRf5o6KvTCRHAoIERgF1Ll7cw+OrTOz9dkH8OlO9p8Q7beZEEiXu/tzHR3g7k1AHTAyZ/OPcs77bA/qFumUAkIE6oHR+QPOufIGhCfnvDXa3XfLPoDbOjnXW9F+1cCPgOM7K87MKoExQEPO5ktzzntEZ58h0hMKCBH4O7AVOK29HfIGhFf29oTuvhX4d+BgM2v3vJGZQDPwdG/PK9IdCgjZ5bn7BuAq4HozO83MBptZpZmdbGbfjvG824DvRufeiZmNNLNzgeuBb2mAWpKmgBAB3P27wOXAfwBrgTeBS4B7Yz71zcBkM/tEzrZFZraJMHj9OeAydy8YIiJxMt0wSEREClELQkREClJAiIhIQQoIEREpSAEhIiIF9Zm1mEaPHu1TpkxJuwwRkbLyzDPP1Ln7mELv9ZmAmDJlCjU1NWmXISJSVszsjfbeUxeTiIgUpIAQEZGCFBAiIlJQnxmDEBEplqamJmpra9myZUvapRRNVVUVkyZNorKyssvHKCBERPLU1tYybNgwpkyZgpmlXU6vuTv19fXU1tYyderULh8XaxeTmc2Ibsy+3MyuKPD+5Wb2gpktjm7OvmfOey1mtjB6zIuzThGRXFu2bGHUqFF9IhwAzIxRo0Z1u0UUWwvCzCoIyxSfCNQCC8xsnru/kLPbc8B0d8+Y2eeBbwNnR+9tdvfD4qpPRKQjfSUcsnry54mzBXEUsNzdX4vWvZ9LuPFJG3d/zN0z0csngUkx1lO6lt4Dm9amXYWIyA7iDIiJhDX1s2qjbe25AHgg53WVmdWY2ZPt3XHLzC6M9qlZu7ZMv2A3r4M7zodnbkm7EhEpMWvWrGHWrFnsvffeHHnkkZxyyim8/PLLO+xzzTXX8J3vfGeHbVOmTKGurq7X5y+JQWoz+wwwHfhwzuY93X2Vme0F/NHMnnf3V3OPc/c5wByA6dOnl+eNLRqj/4nv1qZbh4iUFHfn9NNPZ/bs2cydOxeARYsW8fbbb7PffvslUkOcAbEK2CPn9aRo2w7M7ATgSuDD0X16AXD3VdHP18zsT8DhwKv5x5e9toB4K906RKSkPPbYY1RWVnLRRRe1bTv00EMTrSHOgFgA7GtmUwnBMAv4dO4OZnY4cAMww93fydk+Asi4+1YzGw0cSxjA7nsyCgiRkvbAFbDm+eJ+5riD4eTrOtxlyZIlHHnkkV36uO9///vceuutba/feqs43yexBYS7N5vZJcBDQAVws7svNbNrgRp3nwf8DzAUuCMaYV/p7qcC7wFuMLNWwjjJdXlXP/Udmeg+9O/u1LgSEemSyy67jC9/+cttr4u1snWsYxDuPh+Yn7ftqpznJ7Rz3BPAwXHWVjKyXUyb10HTZqgclG49IrKjTn7Tj8tBBx3EnXfeudP2K6+8kvvvvx+AhQsXxlqD1mJKW7YFAepmEpE2xx9/PFu3bmXOnDlt2xYvXsyMGTNYuHBh7OEACoj0NeZciqaAEJGImXHPPffwyCOPsPfee3PQQQfx1a9+lXHjxiVWQ0lc5rpLy9TDwOGwdYMCQkR2MGHCBG6//fYO97nmmmt22rZixYqinF8tiLRl6sIVDaCBahEpKQqItDXWw257QNVwtSBEpKQoINLkHloQg0dB9UQFhEgJcS/PxRna05M/jwIiTU0ZaN4CQ0ZD9QR1MYmUiKqqKurr6/tMSGTvB1FVVdWt4zRInabsFUyDo4BYvTjdekQEgEmTJlFbW0vZLgJaQPaOct2hgEhTdpmNbBdT4zvQvA36D0i3LpFdXGVlZbfuvNZXqYspTY3RJLlsFxPAxtXp1SMikkMBkabsLOrBo7YHhAaqRaREqIspTdkupiGjoWVbeL5RASEipUEBkabGOuhXCQOroTraphaEiJQIdTGlKTsHwiyExIChCggRKRkKiDRlGkL3EoSQ0FwIESkhCog0NUYtiKxh49WCEJGSoYBIUyYvILTchoiUEAVEmhrrt3cxQehi2rgGWprTq0lEJKKASEtLU7gHxOC8gPCWMKNaRCRlCoi0ZCfJDcnrYgJ1M4lISVBApKUxZx2mrLbZ1LqSSUTSp4BIS9syG7ldTGpBiEjpUECkJXeZjazBI6FioFoQIlISFBBpaSzQgmibLKcVXUUkfQqItGRbEING7LhdcyFEpEQoINKSqQ/hUJG3XqKW2xCREqGASEtj3Y7dS1nVE8JNg1pbk69JRCSHAiItmbxZ1FnVE8K9IbJXOYmIpEQBkZb8hfqyNBdCREqEAiItmfpOAkID1SKSrlgDwsxmmNlLZrbczK4o8P7lZvaCmS02s0fNbM+c92ab2SvRY3acdSautbWDLqbsZDm1IEQkXbEFhJlVANcDJwMHAueY2YF5uz0HTHf3Q4A7gW9Hx44ErgbeBxwFXG1medeDlrEt68OifIUGqYeMgX791YIQkdTF2YI4Clju7q+5+zZgLjAzdwd3f8zdM9HLJ4FJ0fOPAQ+7e4O7rwMeBmbEWGuy2pbZKNDF1K9CNw4SkZIQZ0BMBN7MeV0bbWvPBcAD3TnWzC40sxozq1m7dm0vy01QoZVcc2kuhIiUgJIYpDazzwDTgf/pznHuPsfdp7v79DFjxsRTXBzaVnIt0MUE2+dCiIikKM6AWAXskfN6UrRtB2Z2AnAlcKq7b+3OsWWr0EJ9ubLLbbgnV5OISJ44A2IBsK+ZTTWzAcAsYF7uDmZ2OHADIRxyb6P2EHCSmY2IBqdPirb1DYXuBZGregI0ZcJgtohISvp3vkvPuHuzmV1C+GKvAG5296Vmdi1Q4+7zCF1KQ4E7zAxgpbuf6u4NZvYNQsgAXOvuDXHVmrhMA1QOgcpBhd/PnQuRv5ifiEhCYgsIAHefD8zP23ZVzvMTOjj2ZuDm+KpLUaau/QFq2PHGQWMPSqYmEZE8JTFIvctpb6G+rGHjw09dySQiKVJApCHTzjpMWcPGAaa5ECKSKgVEGjIN7V/BBFBRCUPHqgUhIqlSQKShvZVcc1VPUAtCRFKlgEjatkZo3txxCwIUECKSOgVE0jqbA5Gle1OLSMoUEElrW6ivCy2Ire/C1o3x1yQiUoACImltC/V1FhDZuRBak0lE0qGASFqXu5h061ERSZcCImmZ7gaExiFEJB0KiKRl6qFfJVQN73i/ttnUCggRSYcCImnZORBhccL2VVaF/dTFJCIpUUAkLVPf+QB1luZCiEiKFBBJa6yDwSO7tq/mQohIihQQScvUdz4HIkv3phaRFCkgkpap614X0+YGaNocb00iIgUoIJLU0gRbNnSjBZFz4yARkYQpIJKUie6a2uUxCM2FEJH0KCCSlJ0k1+UuJrUgRCQ9CogktS2z0cWAyE6W26iAEJHkKSCS1N0WxMChYca1WhAikgIFRJLaxiA6WYcpl+ZCiEhKFBBJynYxDeriIDVoLoSIpEYBkaRMHQwaARX9u37MsPFqQYhIKhQQSWqs6/oAdVb1RNj0DjRvi6cmEZF2KCCSlKnv3vgDRHMhHDatiaUkEZH2KCCS1J2VXLM0F0JEUqKASFL2XhDdoVuPikhKFBBJaW3tYQtCy22ISDoUEEnZugG8pfstiKrhUDlEASEiiVNAJKWxPvzs7lVMZrqznIikItaAMLMZZvaSmS03sysKvP8hM3vWzJrN7Ky891rMbGH0mBdnnYloW2ajmy0IUECISCq6MWOre8ysArgeOBGoBRaY2Tx3fyFnt5XA+cCXC3zEZnc/LK76EtfdhfpyVU+E1/9S3HpERDoRW0AARwHL3f01ADObC8wE2gLC3VdE77XGWEdpyGS7mHrYgti4GlpboF9FcesSEWlHnF1ME4E3c17XRtu6qsrMaszsSTM7rdAOZnZhtE/N2rVre1FqArq7kmuu6glhgHvTO8WtSUSkA6U8SL2nu08HPg38wMz2zt/B3ee4+3R3nz5mzJjkK+yOxvpwNVLloO4fq0tdRSQFcQbEKmCPnNeTom1d4u6rop+vAX8CDi9mcYnL1PVsgBo0WU5EUhFnQCwA9jWzqWY2AJgFdOlqJDMbYWYDo+ejgWPJGbsoSz1ZhylLy22ISApiCwh3bwYuAR4ClgG3u/tSM7vWzE4FMLP3mlkt8EngBjNbGh3+HqDGzBYBjwHX5V39VH56spJr1uBRUDFALQgRSVScVzHh7vOB+Xnbrsp5voDQ9ZR/3BPAwXHWlrhMPez+np4dq8lyIpKCUh6k7lt6slBfLt16VEQSpoBIwrYMNG/uZUDo1qMikiwFRBJ6MwciKztZzr04NYmIdEIBkYTeLLORVT0RWrZtn5EtIhIzBUQSerPMRpbmQohIwhQQScgGRG+7mEAD1SKSGAVEEtq6mHrRghimFoSIJEsBkYRMHfTrH+4O11NDdwerUAtCRBLTpYlyZrY7YbmLCcBmYAlQ4+59f5nuYsjOgTDr+Wf0q4Bh4xUQIpKYDgPCzI4DrgBGAs8B7wBVwGnA3mZ2J/Bdd3835jrLW6ahd1cwZWkuhIgkqLMWxCnAP7v7yvw3zKw/8HHCHePuiqG2vqM3K7nmqp4Aby/p/eeIiHRBh2MQ7v4VoNbMPlXgvWZ3v9fdFQ6d6c1Cfbmyy21ospyIJKDTQeponOHfEqil7+rNUt+5qidAUwa2rO/9Z4mIdKKrVzE9YmZfNrM9zGxk9hFrZX1FS1P4Qu/NHIgszYUQkQR1dbnvs6OfF+dsc2Cv4pbTB2Uaws+itCCyNw5aDWMP6v3niYh0oEsB4e5T4y6kzyrGQn1ZWm5DRBLUYReTmX2gk/erzWxacUvqY4qxDlPWsHGAqYtJRBLRWQviTDP7NvAg8AywljAPYh/gOGBP4EuxVljuirGSa1ZFJQwdqxaEiCSiw4Bw98uiwegzCfeNHkeYSb0M+Jm7Px5/iWWuGAv15arWbGoRSUanYxDu3mBm1cBi4PnsZuAAM2t094Ux1lf+si2IQUW66Kt6ItS/WpzPEhHpQFcvcz0SuAgYT1iP6V+AGcDPzUxzJDqSqYeq3aCiqxeMdaJ6gloQIpKIrgbEJOAId/+yu3+JEBi7Ax8Czo+ptr4hU1e87iUIAbF1A2zdWLzPFBEpoKsBsTuwNed1EzDW3TfnbZd8xVpmIyt3LoSISIy62u/xG+ApM7svev0J4DYzGwK8EEtlfUWmHkYWcT5h7lyIMfsV73NFRPJ0daLcN8zsAcI9IQAucvea6Pm5sVTWV2TqYdL04n2eltsQkYR0eeQ0CoSaTneU7dyjhfqK2MWUvfXoRgWEiMRLtxyN05b10Npc3EHqyqowK1stCBGJmQIiTo3ZZTaKGBCgS11FJBEKiDgVcx2mXNUTtdyGiMROARGntpVcixwQw7TchojETwERp2Iu1JeremJonTRtKe7niojkiDUgzGyGmb1kZsvN7IoC73/IzJ41s2YzOyvvvdlm9kr0mB1nnbHJtiCK3sWkK5lEJH6xBYSZVQDXAycDBwLnmNmBebutJCzVcVvesSOBq4H3AUcBV5vZiLhqjU2mASoHw4DBxf1czYUQkQTE2YI4Clju7q+5+zZgLjAzdwd3X+Hui4HWvGM/Bjzs7g3uvg54mLA4YHkp9jIbWW3LbSggRCQ+cQbERODNnNe10baiHWtmF5pZjZnVrF27tseFxiZTV/wBagj3hABdySQisSrrQWp3n+Pu0919+pgxY9IuZ2dxtSAGDoOBw9WCEJFYxRkQq4A9cl5PirbFfWzpyDQUf4A6S5PlRCRmcQbEAmBfM5tqZgOAWcC8Lh77EHCSmY2IBqdPiraVl2LfCyKXAkJEYhZbQLh7M3AJ4Yt9GXC7uy81s2vN7FQAM3uvmdUS7nd9g5ktjY5tAL5BCJkFwLXRtvKxLQNNGbUgRKRsFek+mIW5+3xgft62q3KeLyB0HxU69mbg5jjri1XbLOq4WhATYdPb0NIEFZXxnENEdmllPUhd0uJahymregLgsHFNPJ8vIrs8BURc4lrJNUuT5UQkZgqIuMTexZRz61ERkRgoIOLSGNM6TFlqQYhIzBQQccnUQ7/+UDU8ns+v2i2s86SAEJGYKCDikqkLrQezeD7fLLrUVV1MIhIPBURcGuvjG6DO0lwIEYmRAiIumfp4FurLVT1RASEisVFAxCXbxRSn6gmwcTW0tsR7HhHZJSkg4hLXSq65qieAt0BjCS51LiJlTwERh5Ym2LI+vjkQWW03DtJAtYgUnwIiDpvXhZ9JdDGBxiFEJBYKiDjEPUkuS7ceFZEYKSDiEPcyG1mDR0HFAHUxiUgsFBBxaGtBxBwQZjBsvFoQIhILBUQcskt9x92CAM2FEJHYKCDikA2IQSPiP5eW2xCRmCgg4tBYFxbTS+JOb9nlNtzjP5eI7FIUEHHI1CXTvQShi6ll2/ZWi4hIkSgg4pBJYKG+LN04SERiooCIQ2N9/HMgstrmQqxO5nwisstQQMQhUxf/Sq5Z2RbEhjeTOZ+I7DIUEMXmnmwX09DdYbfJ8PgPt8+/EBEpAgVEsW3ZAK3NyQ1S96uAT/4SNr0Dd5wfFgoUESkCBUSxZa8mSmoMAmDiEfCJH8KKv8If/jO584pIn9Y/7QL6nKSW2ch32DmwehE89VMYf2h4LSLSC2pBFFvbQn0JtiCyTvoGTPkg/O4LsOrZ5M8vIn2KAqLY2rqYEm5BQJi5/clbwsD1bz8TxiVERHpIAVFsSd0Loj1DRsPZt4ag0qC1iPSCAqLYMvVQORgGDE6vhgmHwak/hjceh4e+ll4dIlLWYg0IM5thZi+Z2XIzu6LA+wPN7LfR+0+Z2ZRo+xQz22xmC6PHz+Kss6ga69LpXsp3yKfg/ZfA03PguVvjPVdLkxYLFOmDYgsIM6sArgdOBg4EzjGzA/N2uwBY5+77AN8HvpXz3qvuflj0uCiuOosuU5/OAHUhJ3wdpn4Yfn8Z1D5T/M9vbYUnfwbXTYa/fb/4ny8iqYqzBXEUsNzdX3P3bcBcYGbePjOBX0bP7wQ+amYWY03xy9SlN/6Qr6I/nPULGDYuDFpvfLt4n732ZfjFDHjw3wGDp38OrS3F+3wRSV2cATERyF0gqDbaVnAfd28GNgDZb9epZvacmf3ZzD5Y6ARmdqGZ1ZhZzdq1a4tbfU81JrjMRlcMGQWzboPN6+CO2dC8rXef19IMf/0e/OwDsPYlOP0GOP2nsPEtWP5ocWoWkZJQqoPUq4HJ7n44cDlwm5lV5+/k7nPcfbq7Tx8zZkziRRaU5L0gumrcwTDzJ7Dy7/DgTkNBXbdmCdz4UXj067DfSXDx03DoLNjv5BCKz/6y888QkbIRZ0CsAvbIeT0p2lZwHzPrDwwH6t19q7vXA7j7M8CrwH4x1loc2zLQlCmdLqZcB58Fx1wKNTfBs7/q3rHN2+Cxb8KcD4f7Tnzyl+FS2mFjw/v9B4SgePlBzb0Q6UPiDIgFwL5mNtXMBgCzgHl5+8wDZkfPzwL+6O5uZmOiQW7MbC9gX+C1GGstjjTWYeqOE66BvY6D+78Eby7o2jGrngnB8OdvwbQzQ6vhoNN23u+I88IihYvmFrNiEUlRbAERjSlcAjwELANud/elZnatmZ0a7XYTMMrMlhO6krL9Hx8CFpvZQsLg9UXu3hBXrUXTtsxGiXUxZfWrgLNuhmHjo0HrNe3v27Q5LPx34wmweT2c81s4Yw4MHll4/zH7w6SjQutEl7yK9AmxLtbn7vOB+Xnbrsp5vgX4ZIHj7gLuirO2WDSmuMxGVw0eGQatbzoRbj8PZv8+dBHleuMJuO8SaHgVjpgd1niqGt75Zx9xHsy7BN58CiYfHU/9IpKYUh2kLk/ZLqZSbUFkjZsGM68PX+QP/Nv27Vs3wfyvwC9OhtYmOO8+OPVHXQsHgINOhwFD4dlfx1O3iCRKy30XU7aLqb1umFIy7QxYszhMcBt/KIzYE+Z9Idy69H0XwfH/CQOHdu8zBw4NIbHkLpjx31C104VnIlJGFBDF1FgH/fpD1W5pV9I1x/8nrHk+DFp7C4zaB/7pwd51Dx0xG577NSy9G448v2ilikjy1MVUTJn6cAVTuUwG71cBZ94IUz4AH7gMLvpb78cOJk2HMQeom0mkD1ALopgyJTaLuisGjYDZ+Vcf94JZGKx+6Gvw9gswNn/5LREpF2pBFFNjXXmMP8TtkFnQrzJ0NYlI2VJAFFMpLrORhiGj4IBTwqS55q1pVyMiPaSAKKZy7GKKyxHnweYGePH+tCsRkR5SQBRLS3NYMVUtiGCv46B6krqZRMqYAqJYNkcrgZTqOkxJ61cBh58Lrz4G61emXY2I9IAColgas5PkFBBtDjs3/HzuN+nWISI9ooAolnJZZiNJI/aEvT4CC3+ju82JlCEFRLG0LbOhgNjBEZ8Ny3e89qe0KxGRblJAFIu6mAo74ONhMl53b1IkIqlTQBRL282CNFFuB/0HholzL96/fTl0ESkLCohiydSHRfoqKtOupPQc8dmwfPji36ZdiYh0gwKiWBo1i7pdYw+CiUfqbnMiZUYBUSyZOo0/dOTwz8LaZeEe1yJSFhQQxdKoZTY6NO1MqByswWqRMqKAKJZMfVikTgqrqt5+t7mtm9KuRkS6QAFRDO5aqK8rDv8sbNsEL9ybdiUi0gUKiGLYsiFcpaMxiI5NPhpG7au7zeXbugm2ZZI/rzu88jD85pPwf+fAwtsg05B8HVKydEe5YtAyG11jFi55ffgqWPsSjNm/95/55tPw+8tg/GEw45tQNbz3n5kU9zAm8+AV4QZLh50D0y+AMfvFe97mrfD8HfDET8KFA8MmgPWDl+aHe6pP+SC85xNhkuOwsfHWIiVNLYhiaJskp4Do1KHnhC+h3i4D3rwNHv0G3PwxaFwLi26Dnx4LK/5WnDrj1lgPv/0M/O7ScAnwfifBgpvg+vfCLR+HpfdCS1Nxz7l5Hfz1e/CDQ+C+i8OKu6fPgS8sgsuWwD//EY7517D67v2Xw3f3h5s+Bn+/Hta9UdxapCyY95Hr0qdPn+41NTXpnPzF+TD3HLjwTzDh8HRqKCdzz4WVT8Lly6D/gO4f/86LcM+FsHpRWDF2xnWw9kW451+g4XU45hI47j+gsqr4tRfD8kfh3v8XfrE44Wo4+mLo1w82rYXnfgU1t8CGlTB0HBx5Phw5G6on9Px8696AJ38aWitNjeFeHcdeGn6a7by/e/jv+cI8WPY7ePv5sH38oaFl8Z6Z8bdyJDFm9oy7Ty/4ngKiCJ79Fcz7V/ji87Db5HRqKCcvPwS3fQo+9Ws48NSuH9faCk/9DB65BgYOhU/8MHxhZW1rhD/8B9TcDLsfCGfMgXEHF738HmvaAo9+HZ78Xxi9P5x5I4w/ZOf9WlvC2MCCG2H5I6H754BT4L2fg6kfLvylXshbz8ETPw6tETOYdlYIz+7+N2l4DZb9HpbNg9oFYdvo/aOw+EQIjq7WVEwtzfDmk/DKH6D/IJh2RnG6LXcxCoi4/fV74R/+11bDgMHp1FBOWprhB9PCF9W5d3TtmA21cO/n4fW/wH4z4NQfw9DdC+/78h9g3iVhwPX4K+GYS0N3SprefgHu+hy8sxSOuhBOvBYqB3V+XMPr8MwvwsD+5oYwyD/9n8J4xaARO+/f2hpC5YkfwYq/woBhMP18eN/nYfjE3v853n0rrKu1bF7ozvPW8EvRvifBnsfAnsfCsHG9P097tmwIf76XHgzBsGV9GL9pbQYcxk4LQTHtTBgxJb46+hAFRNweujL81nrl6nTOX44e/Qb87XvwxSUdf3G5hwHV+78cvgRm/He433Vnv7E21sPvvxi+yCa/H07/WTpfGK2t8PQN8PDVYS7IzP8N4w3d1bQlXB684CaofTr8xnzwWfDeC0K3ZtvA849D91D1RDj68+G/VVwD9431YWB72e/gjcfDJcwAI/feHhZ7HhMCpDctjIbX4eUHw7neeCL8PRg0MvyisP8M2Pv40Hp84b4wz+bNp8JxE6eHoDjotN510fVxCoi43XMRrHgcLns+nfOXo4bX4UeHhbGCD3+l8D6ZhjBYuvQe2ON94Ut+5F5dP4d7WCBw/lfCb7ozroPDP5Ncd8jGNWGs4dVHo1bPT2DomN5/7upFISievwOaMiEg3l0Nm9bA2IPDQPO0M5JdOLKlGdYsCl/g2ceW9eG96kkhKKYcG0Jj1D4d/z9obYHaGnj5gdBSWLssbB+9P+x/cnhMem/7rcL1K2HJ3SEs1iwGLJx32hlw4Gma0JpHARG3W88KV9L8y5/TOX+5uuXj4R/zpQvDIG2u5Y/AvReHNa6O+xoc+8WedxOtXxm+qFf8Ffb/hzB2UYwv6o68eH8Yl9qWgY/9V7h8tdjBtGUDLJoLz90KQ8aE8YX2Bp6T1toavtjfeCK0LlY8Do3vhPeGjMlpYRwbxouaMvDqH0NL4eWHwv93qwj77X9KaCl055eDrLpXorC4E+peDp+593GhZXHAP5TXZdExUUDEbc5HwiS5z9yVzvnL1eI74O7PwXn3hVuTQvhCffgqWPBzGHNAGGgef2jvz9XaCk/9FB75eujq+cSPwsBvsW1rhIe+Bs/cAuMOgTNv0hU/EFpz9a+GsMi2MDasDO9VDYemzdCyLTzf58TQStjnBBi0W/HO//bSEBRL7gq/NFQMCGMn084Ilxr3qwyXYPfrH34ZaXsevS6F4I1BagFhZjOAHwIVwI3ufl3e+wOBXwFHAvXA2e6+Inrvq8AFQAtwqbs/1NG5Ug2IHxwMk4+BM25I5/zlqmlzuNZ+nxPhrJug9plw+Wr98nDp50evKv6lqm+/AHdfGC7dPOI8+Ng3YeCw4nz2qmfh7n8OX4THfgGOu7Jnl/HuKtavhDf+DiufgMohoZUw+f3xd425h1WFl9wVWheb1nTtOCsQGtnng0aEsZYRe4afu+25/XmxWynNW2HT26FbcePq0JU5cGjoPu2BjgIitpnUZlYBXA+cCNQCC8xsnru/kLPbBcA6d9/HzGYB3wLONrMDgVnAQcAE4BEz28/dW+Kqt1ca67XMRk9UDoJDzoZnfhkuXX38RzBsPJw3D/b6cDznHHtgmBD2p/+Gx38Ar/05tFImH93zz2xtCZ/12Ddh6FiY/TuY+sFiVdx37TY5PA49O9nzmsGk6eFx0n+FOTnrVoTB79bm8P+z7Xmh13mPluZwhdm6N0I3ZnawPqtqt7zgmLL9+W6Tt1/52NoSuqqzX/obV+8YAhuj59mJubkmHN7jgOhInEttHAUsd/fXAMxsLjATyA2ImcA10fM7gZ+YmUXb57r7VuB1M1sefd7fi15lpgF+cXLPj3cPk4808NUzh38Wnp4Df/t+CIuTv128boX29B8QJqjt97Ewue4XJ4fLR3vahbB1E7xbCwedAR//XuHLT6U09asIg+dTji3O57mHGevrVoQW0vo3QnCsfyMsL/PKw9C8ZcdjhowJ3Vub3oadfge2cDn3sPEwfI8wOF89IVxKPGx89HNCbLc6jjMgJgJv5ryuBd7X3j7u3mxmG4BR0fYn847d6VpIM7sQuBBg8uQeTlDrV9H7yTXjDob3dGPCl2w3/hD46NXhypbuTJorhslHw0V/g798B9a93osPMjjgmnDZaR/tp5YuMgtf1oNHwsQjdn6/tTUM1q9fGQXHivCztQWqx0df+tGjejwM2R0q0lsyr6wX63P3OcAcCGMQPfqQquHwKd3EJlUfvDy9cw8cBid+Pb3zy66lX7/ot/5xsMdRaVfTqTgX61sF7JHzelK0reA+ZtYfGE4YrO7KsSIiEqM4A2IBsK+ZTTWzAYRB53l5+8wDZkfPzwL+6OGyqnnALDMbaGZTgX2Bp2OsVURE8sTWxRSNKVwCPES4zPVmd19qZtcCNe4+D7gJ+HU0CN1ACBGi/W4nDGg3AxeX7BVMIiJ9lCbKiYjswjqaB6EbBomISEEKCBERKUgBISIiBSkgRESkoD4zSG1ma4GO7qw+GqhLqJzeKJc6QbXGpVxqLZc6QbV2ZE93L7j+fZ8JiM6YWU17I/WlpFzqBNUal3KptVzqBNXaU+piEhGRghQQIiJS0K4UEHPSLqCLyqVOUK1xKZday6VOUK09ssuMQYiISPfsSi0IERHpBgWEiIgU1OcDwsxmmNlLZrbczK5Iu572mNkeZvaYmb1gZkvN7Atp19QRM6sws+fM7Pdp19IRM9vNzO40sxfNbJmZvT/tmtpjZpdF/++XmNn/mVlV2jVlmdnNZvaOmS3J2TbSzB42s1einyVxr9V2av2f6O/AYjO7x8x2S7HEbE071Znz3pfMzM1sdBq1ZfXpgDCzCuB64GTgQOAcMzsw3ara1Qx8yd0PBI4GLi7hWgG+ACxLu4gu+CHwoLsfABxKidZsZhOBS4Hp7j6NsET+rHSr2sEtwIy8bVcAj7r7vsCj0etScAs71/owMM3dDwFeBr6adFEF3MLOdWJmewAnASuTLihfnw4I4Chgubu/5u7bgLnAzJRrKsjdV7v7s9HzjYQvsp3uw10KzGwS8A/AjWnX0hEzGw58iHDfEdx9m7uvT7WojvUHBkV3VxwMvJVyPW3c/S+Ee7bkmgn8Mnr+S+C0JGtqT6Fa3f0P7t4cvXyScJfKVLXz3xTg+8C/AalfQdTXA2Ii8GbO61pK9Es3l5lNAQ4Hnkq5lPb8gPAXuDXlOjozFVgL/CLqDrvRzIakXVQh7r4K+A7ht8bVwAZ3/0O6VXVqrLuvjp6vAcamWUw3/BPwQNpFFGJmM4FV7r4o7Vqg7wdE2TGzocBdwBfd/d2068lnZh8H3nH3Z9KupQv6A0cAP3X3w4FGSqcbZAdR//1MQqhNAIaY2WfSrarrolsFp/4bb2fM7EpCd+5v0q4ln5kNBr4GXJV2LVl9PSBWAXvkvJ4UbStJZlZJCIffuPvdadfTjmOBU81sBaHL7ngzuzXdktpVC9S6e7YldichMErRCcDr7r7W3ZuAu4FjUq6pM2+b2XiA6Oc7KdfTITM7H/g4cK6X5gSwvQm/ICyK/n1NAp41s3FpFdTXA2IBsK+ZTTWzAYRBv3kp11SQmRmhr3yZu38v7Xra4+5fdfdJ7j6F8N/zj+5ekr/puvsa4E0z2z/a9FHCfc5L0UrgaDMbHP1d+CglOqCeYx4wO3o+G7gvxVo6ZGYzCN2ip7p7Ju16CnH35919d3efEv37qgWOiP4ep6JPB0Q0KHUJ8BDhH9vt7r403aradSzwWcJv5AujxylpF9UH/CvwGzNbDBwGfDPdcgqLWjl3As8CzxP+bZbOkgtm/wf8HdjfzGrN7ALgOuBEM3uF0AK6Ls0as9qp9SfAMODh6N/Wz1ItknbrLClaakNERArq0y0IERHpOQWEiIgUpIAQEZGCFBAiIlKQAkJERApSQIiISEEKCJEEWaB/d1IW9BdVJGZmNiW6J8mvgCXsuPyLSMnSRDmRmEWr874GHOPuT6ZcjkiXqQUhkow3FA5SbhQQIsloTLsAke5SQIiISEEKCBERKUiD1CIiUpBaECIiUpACQkREClJAiIhIQQoIEREpSAEhIiIFKSBERKQgBYSIiBT0/wGf/kXLPs3BYgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzSklEQVR4nO3deXxU9bn48c9DdsgGCVsWBVlUSKzViAjqpbXtRWvBumtbVECvrd7qrVzF9neptcu9tr/W3vurta8Wd23F4oYVRKtdLmtBi0JAICBqAkhIJoGQzGR7fn/MTDoMk33OzJnkeb9eeTFzzpnzfaKEJ+e7PF9RVYwxxphwQ+IdgDHGGHeyBGGMMSYiSxDGGGMisgRhjDEmIksQxhhjIrIEYYwxJiJLEMYYYyKyBGEMICI3ishWEWkUkYMi8rCI5HbzmT+LyMKwY7NEpLKLz+wTkSYRaQi087iIZIacf1xEmkXkaOBrm4j8p4jkhMXaFrhH8OsX/fj2jYnIEoQZ9ETkLuAB4N+BHGA6cDLwhoikOtDkl1Q1EzgT+DRwb9j5H6tqFjASuCkQz1oRGRZyzXpVzQz5ut2BOM0gZwnCDGoikg18D/hXVX1NVVtUdR9wNTAO+KpTbavqQWA1/kQR6bxXVTcBc4A8/MnCmJixBGEGuxlAOvBC6EFVbQBWAp93qmERKQIuBiq6uk5VjwJvABc4FYsxkViCMINdPnBYVVsjnDsQON+V/xGRuuAX8IcetPmSiBwFPgYOAd/twWf2AyNC3k8PbVdEpvfgHsb0iiUIM9gdBvJFJDnCubGB84jIr0IGhL8dcs03VTU3+AVc2oM2LwuMMcwCTqP7JARQCNSGvN8Q2q6qbujBPYzpFUsQZrBbD/iAy0MPBmYWXQy8CaCqt4YMCP8oGg2r6l+Ax4H/29V1gVg+B/xvNNo1pqcsQZhBTVXr8Q9S/z8RmS0iKSIyDngOqASecjiEnwOfF5FPhZ8QkTQRORt4CfAAjzkcizHHsQRhBj1V/THwbfy/yR8BNuIfH7hIVX0Ot10NPAksCTl8d2CMoiZw7m1ghqoeczIWY8KJbRhkjDEmEnuCMMYYE5ElCGOMMRFZgjDGGBORJQhjjDERRVoclJDy8/N13Lhx8Q7DGGMSyttvv31YVUdGOjdgEsS4cePYvHlzvMMwxpiEIiIfdnbOupiMMcZEZAnCGGNMRJYgjDHGRDRgxiCMMSYaWlpaqKysxOv1xjuUqEpPT6eoqIiUlJQef8YShDHGhKisrCQrK4tx48YhIvEOJypUlZqaGiorKxk/fnyPP2ddTMYYE8Lr9ZKXlzdgkgOAiJCXl9frpyJLEMYYE2YgJYegvnxPliASnKry4u4XaWxpjHcoxpgBxhJEgttTt4cl65bw8p6X4x2KMSaKKisrmTt3LpMmTWLChAnccccdNDc3n3Dd448/zu23337csVmzZkVl4bAliARX460BYNvhbXGOxBgTLarK5ZdfzmWXXcbu3bvZtWsXDQ0NfOc734lpHJYgEpzH6wGg/HB5nCMxxkTLW2+9RXp6OjfddBMASUlJPPjggzz66KM0NsauO9mmuSa4Wm8tAHvr99LY0sjQlKFxjsiYgeOBvz3A+7XvR/Wep404jXum3dPlNeXl5Zx99tnHHcvOzuakk06ioqKCM84447hzy5YtY82aNR3vKyoqohKrPUEkuDpfHQCKsr1me3yDMcbExTXXXMOWLVs6vsrKyqJyX3uCSHC13lpShqTQ0t5CeU05ZWOi8xfDGEO3v+k7ZcqUKSxfvvy4Y0eOHOGjjz7if//3f5k3bx4AK1eudDQOe4JIcB6vh8LMQsYOG2sD1cYMEBdddBGNjY08+eSTALS1tXHXXXdx4403ctttt3U8KRQUFDgahyWIBOfxeRiePpyS/BLKa2yg2piBQER48cUX+f3vf8+kSZOYPHky6enp/OhHP4ppHNbFlOA8Xg8nZZ3ElLwpvPHhG9T76slJy4l3WMaYfiouLuaVV17p9robb7yRG2+88bhjf/7zn6MSgz1BJDiP9x9PEGDTXY0x0eNoghCR2SKyU0QqRGRxhPNpIrIscH6jiIwLOXeGiKwXkXIR2Soi6U7GmojatZ06Xx3D04czJW8KANtqbBzCGBMdjiUIEUkCHgIuBqYA14nIlLDLFgAeVZ0IPAg8EPhsMvA0cKuqTgVmAS1OxZqojjYfpU3bGJ42nOzUbMZlj7MnCGOiQFXjHULU9eV7cvIJYhpQoap7VbUZeBaYG3bNXOCJwOvlwEXiLzn4BeA9VX0XQFVrVLXNwVgTUnAV9fD04QBMyZtiTxDG9FN6ejo1NTUDKkkE94NIT+9dR4yTg9SFwMch7yuBczu7RlVbRaQeyAMmAyoiq4GRwLOq+uPwBkTkFuAWgJNOOinq34DbeXzHJ4iS/BJWfrCS6sZqRg4dGc/QjElYRUVFVFZWUl1dHe9Qoiq4o1xvuHUWUzJwPnAO0Ai8KSJvq+qboRep6q+BXwOUlZUNnHTfQ8EyG6EJAqC8ppxZQ2fFKyxjElpKSkqvdl0byJzsYqoCikPeFwWORbwmMO6QA9Tgf9r4q6oeVtVGYCVwloOxJqRgF9OItBGAv8bLEBliC+aMMVHhZILYBEwSkfEikgpcC6wIu2YFcEPg9ZXAW+rv+FsNlIrI0EDi+CfACg2FCdZhyk3PBSAjOYMJuRNsHMIYExWOJQhVbQVux/+P/Q7gOVUtF5H7RWRO4LJHgDwRqQC+BSwOfNYD/Ax/ktkCvKOqrzoVa6Kq9daSkZxBRnJGx7GSvBLKD5cPqAE2Y0x8ODoGoaor8XcPhR5bEvLaC1zVyWefxj/V1XTC4/UwPG34ccdK8kt4seJF9h/bT2FmYZwiM8YMBLaSOoEFV1GHmpo3FbAd5owx/WcJIoEFC/WFmjx8MilDUmzBnDGm3yxBJLBIXUwpSSmcOvxUG6g2xvSbJYgEFqmLCWBq/lS212ynXdvjEJUxZqCwBJGgGlsa8bZ5IyeIvKkcaznGviP7Yh+YMWbAsASRoIJrIEakjzjhnJX+NsZEgyWIBBVcRZ2blnvCuVNyTiEjOcN2mDPG9IsliAQVrMMU6QkiaUgSp4843aa6GmP6xRJEggp2MUUagwD/QPX7te/T0m7baBhj+sYSRIIKr+QariSvBF+bjz11e2IZljFmALEEkaA8Xg/JkkxWSlbE8zZQbYzpL0sQCcrj85Cbnot/A74TFWcVk5WaZQvmjDF9ZgkiQXW2SC5IRJiaN9WeIIwxfWYJIkF5vJ6OjYI6U5Jfwm7PbnxtvhhFZYwZSCxBJKhgF1NXSvJKaNVWdtbujE1QxpgBxRJEgqr11p5QqC/c1Hwr/W2M6TtLEAmopb2Fo81HIy6SCzV66Gjy0vNsRbUxpk8sQSSgel890PkaiCARoSS/xAaqjTF9YgkiAQUXyXU3BgH+yq576/dyrOWYw1EZYwYaSxAJKFior7tZTOAfh1CU7TXbnQ7LGDPAWIJIQB6fP0F018UEtqLaGNN3jiYIEZktIjtFpEJEFkc4nyYiywLnN4rIuMDxcSLSJCJbAl+/cjLORBN8guhJghiRPoKCYQU2UG2M6bVkp24sIknAQ8DngUpgk4isUNXQvo4FgEdVJ4rItcADwDWBc3tU9Uyn4ktkwQSRk5bTo+un5k+1qa7GmF5z8gliGlChqntVtRl4Fpgbds1c4InA6+XARdJZcSHTodZbS3ZqNilDUnp0/dS8qVQ2VFLnrXM2MGPMgOJkgigEPg55Xxk4FvEaVW0F6oG8wLnxIvJ3EfmLiFzgYJwJp85X1+0aiFAd4xDWzWSM6QW3DlIfAE5S1U8D3wJ+KyLZ4ReJyC0isllENldXV8c8yHjprlBfuNPzTgcsQRhjesfJBFEFFIe8Lwoci3iNiCQDOUCNqvpUtQZAVd8G9gCTwxtQ1V+rapmqlo0cOdKBb8Gdar21Efei7kx2ajbjssfZOIQxplecTBCbgEkiMl5EUoFrgRVh16wAbgi8vhJ4S1VVREYGBrkRkVOAScBeB2NNKB6vp1ddTOAfqLaprsaY3nAsQQTGFG4HVgM7gOdUtVxE7heROYHLHgHyRKQCf1dScCrshcB7IrIF/+D1rapa61SsiURVqfPV9aqLCfyVXQ81HeJQ4yGHIjPGDDSOTXMFUNWVwMqwY0tCXnuBqyJ87nngeSdjS1RHmo/Qpm3dVnINF6zsWn64nFEnjXIiNGPMAOPWQWrTid4skgt12ojTSJIk24LUGNNjliASTG/KbITKSM5gQu4EG4cwxvSYJYgE09cnCPCvhyivKUdVox2WMWYAsgSRYHpTyTXc1Lyp1PnqqGoIn21sjDEnsgSRYIJdTD3ZCyJcxxakNg5hjOkBSxAJptZbS0ZyBhnJGb3+7OTcyaQMSbFxCGNMj1iCSDB13rpeT3ENSklK4dThp1rJDWNMj1iCSDC1vto+DVAHTc2fyvaa7bRrexSjMsYMRJYgEozH6+nT+ENQSX4Jx1qOsa9+X9RiMsYMTJYgEozH6+nTDKagkjwr/W2M6RlLEAmmL3WYQo3PGU9GcoZVdjXGdMsSRAJpam2iqbWpXwkiaUgSp4843aa6GmO6ZQkigXSsou7jLKagkvwSdtbupKW9JRphGWMGKEsQCaSvdZjCleSX4GvzsaduTzTCMsYMUJYgEkhHmY1ebhYUbmpeYEW1jUMYY7pgCSKB9KdQX6jirGKyU7MtQRhjumQJIoHUev2b6vVmP+pIRISpeVNtqqsxpkuWIBJIna+OZEkmOzW73/eamj+VCk8F3lZvFCIzxgxEliASSHAVtYj0+14leSW0ais7PTujEJkxZiCyBJFAar39q8MUqqP0t41DGGM6YQkigXi8nn6vgQgaPXQ0+Rn5bK/ZHpX7GWMGHksQCaS/ZTZCBQeq7QnCGNMZRxOEiMwWkZ0iUiEiiyOcTxORZYHzG0VkXNj5k0SkQUQWORlnoqj11kbtCQL83Uwf1H/AsZZjUbunMWbgcCxBiEgS8BBwMTAFuE5EpoRdtgDwqOpE4EHggbDzPwNWORVjImlpb+FI85F+L5ILVZJXgqLWzWSMicjJJ4hpQIWq7lXVZuBZYG7YNXOBJwKvlwMXSWCKjohcBnwA2GR9oN5XD/RtL+rOBAeqbQtSY0wkTiaIQuDjkPeVgWMRr1HVVqAeyBORTOAe4HtdNSAit4jIZhHZXF1dHbXA3Shaq6hDjUgfQcGwAqvsaoyJyK2D1PcBD6pqQ1cXqeqvVbVMVctGjhwZm8jipKMOUz82C4pkar4NVBtjInMyQVQBxSHviwLHIl4jIslADlADnAv8WET2AXcC3xaR2x2M1fVqff4yG9F8ggB/ZdeqhqqOBGSMMUFOJohNwCQRGS8iqcC1wIqwa1YANwReXwm8pX4XqOo4VR0H/Bz4kar+wsFYXc+JLib4R2VXq8tkjAnnWIIIjCncDqwGdgDPqWq5iNwvInMClz2Cf8yhAvgWcMJUWONX560DICctJ6r3nTx8MgAf1H8Q1fsaYxJfspM3V9WVwMqwY0tCXnuBq7q5x32OBJdgar21ZKdmkzIkJar3zU3LJTMlk4+Pftz9xcaYQcWtg9QmjMfnieoaiCARoSiriMqjlVG/tzEmsVmCSBAer6ff+0B0piiziMoGSxDGmONZgkgQHp8n6gPUQUVZRVQdraJd2x25vzEmMVmCSBAerzNdTODfgrS5vZnqxoG92NAY0zuWIBKAqlLnjV4l13BFmUUANlBtjDmOJYgEcKT5CK3a6twYRJY/Qdg4hDEmlCWIBFDnqwNwrItp7LCxDJEhNpPJGHMcSxAJwKlV1EEpSSmMHTbWniCMMcexBJEAar3O1GEKVZRZZGMQxpjjWIJIAB1PEFHcTS6cLZYzxoSzBJEAPD5nu5jAnyBqvbU0tjQ61oYxJrFYgkgAHq+HjOQMMpIzHGsjONXVxiGMMUGWIBKAx+txtHsJ/IvlAOtmMsZ06FE1VxEZBcwECoAmYBuwWdVqM8RCra82qntRRxJcC2ED1caYoC4ThIh8Bv8eDSOAvwOHgHTgMmCCiCwHfqqqRxyOc1BzchV1UHZqNlkpWfYEYYzp0N0TxCXAzar6UfiJwBahlwKfB553IDYT4PF6OCXnFEfb6Cj7bWMQxpiALhOEqv67iAwRkatV9bmwc63AS04GZ/ycrOQaqiiriN2e3Y63Y4xJDN0OUgfGGe6OQSwmgqbWJppam2KWIKoarOy3Mcavp7OY/igii0SkWERGBL8cjcwA/9iL2ulZTOCf6trS3sKhxkOOt2WMcb+e7kl9TeDP20KOKeBsx7ih1ud8mY2g0JlMY4aNcbw9Y4y79ShBqOp4pwMxkQXLbDhVyTVUceY/1kKcM+Ycx9szxrhbl11MInJ+N+ezRaQkuiGZUMEE4dReEKHGZI4hSZJsJpMxBuh+DOIKEVknIktE5IsiMk1ELhSR+SLyFPAHoNP6DyIyW0R2ikiFiCyOcD5NRJYFzm8UkXGB49NEZEvg610R+XJ/vslE5nSp71ApQ1IYM2yMrYUwxgDdT3P9t8Bg9BXAVcAY/CupdwC/UtW1nX1WRJKAh/Cvk6gENonIClXdHnLZAsCjqhNF5FrgAfzjHduAMlVtFZGxwLsi8kpgau2g4vF5SJZkslOzY9KeVXU1xgT1ZJprLZANvAe8AawBDgOniciZXXx0GlChqntVtRl4Fpgbds1c4InA6+XARSIiqtoYkgzS8Q+ID0oer4fc9FxEJCbtFWXaYjljjF9Pp7meDdwKjMVfj+lfgNnAb0SkszUShUBoYZ/KwLGI1wQSQj2QByAi54pIObAVuDXS04OI3CIim0Vkc3V1dQ+/lcTi8XpiMv4QFCz7fazlWMzaNMa4U08TRBFwlqouUtW78CeMUcCFwI1OBKaqG1V1KnAOcK+IpEe45teqWqaqZSNHjnQijLjz+DwxmcEUFJzqat1MxpieJohRgC/kfQswWlWbwo6HqgKKQ94XBY5FvCZQ2ykHqAm9QFV3AA3AoJwt5fHGpsxGUEfZb+tmMmbQ6+lCuWeAjSLycuD9l4DfisgwYHsnn9kETBKR8fgTwbXA9WHXrABuANYDVwJvqaoGPvNxYJD6ZOA0YF8PYx1Qar21MVlFHdSxcZA9QRgz6PV0odz3RWQV/j0hwD8msDnw+iudfKZVRG4HVgNJwKOqWi4i9+PfS2IF8AjwlIhUALX4kwjA+cBiEWkB2oFvqOrhPnx/Ca21vZUjzUdi+gSRk5ZDVmqW7QthjOnxEwSBhLC52wuP/8xKYGXYsSUhr734p8+Gf+4p4KnetDUQ1fnqgNisgQhlM5mMMWBbjrpaLBfJhSrKKqLqaPhwkTFmsLEE4WIddZjSYls4tzirmKqGKtra22LarjHGXSxBuJjHF6jD5PB+1OGKsqzstzHGEoSrxbKSa6iOmUw2DmHMoGYJwsWCCSInLSem7dpiOWMMWIJwtVpvLdmp2aQMSYlpu2OG+ct+21RXYwY3SxAuVueri/kMJvCX/R47bKx1MRkzyFmCcDGP1xPTVdShrOy3McYShIvV+mrj8gQBliCMMZYgXM3jjW0l11BFmUV4fB4amhvi0r4xJv4sQbiUqlLnrYvpXhChgjOZqhpsRbUxg5UlCJc62nKUVm2NWxdTR9lv62YyZtCyBOFS8VokFxR8grCprsYMXpYgXCpehfqCslOzyU7NtqmuxgxiliBcqiNBxGmaK9hMJmMGO0sQLhUs1BevJwiwfSGMGewsQbhUrbcWiG+CsLLfxgxuliBcyuP1kJGcQUZyRtxiKMoqorW9lU8aP4lbDMaY+LEE4VJ1vvitgQiyqq7GDG6WIFyq1hu/MhtBti+EMYObJQiX8ng9cU8QY4aNIVmS7QnCmEHKEoRL1fnqYr4XdbjkIcmMzRxri+WMGaQcTRAiMltEdopIhYgsjnA+TUSWBc5vFJFxgeOfF5G3RWRr4M/POhmnG9V6a2O+F3UkRZm2FsKYwcqxBCEiScBDwMXAFOA6EZkSdtkCwKOqE4EHgQcCxw8DX1LVUuAG4Cmn4nQjb6uXptamuJXZCFWUZWshjBmsnHyCmAZUqOpeVW0GngXmhl0zF3gi8Ho5cJGIiKr+XVX3B46XAxkikuZgrK7ihlXUQUVZRdT56jjafDTeoRhjYszJBFEIhHZeVwaORbxGVVuBeiAv7JorgHdU1RfegIjcIiKbRWRzdXV11AKPNzesog4KzmSyst/GDD6uHqQWkan4u53+JdJ5Vf21qpapatnIkSNjG5yD4l2oL1Sw7LcNVBsz+DiZIKqA4pD3RYFjEa8RkWQgB6gJvC8CXgTmqeoeB+N0nY4yGy7pYgJbLGfMYORkgtgETBKR8SKSClwLrAi7ZgX+QWiAK4G3VFVFJBd4FVisqmsdjNGV3PQEkZWaRU5ajqMJ4q+Vf6WxpdGx+xtj+saxBBEYU7gdWA3sAJ5T1XIRuV9E5gQuewTIE5EK4FtAcCrs7cBEYImIbAl8jXIqVrep89WRJElkp2bHOxTA2aquuz27ue3N23h6x9OO3N8Y03fJTt5cVVcCK8OOLQl57QWuivC5HwA/cDI2N6v11pKblouIxDsUwN/NtKNmhyP3Xlu1tuPPW864xZE2jDF94+pB6sHKDWU2QhVnFbO/YT+t7a1Rv/fa/f4E8W71uzaV1hiXsQThQh6fxxWL5IKKMoto1eiX/W5qbeKdT96hJK+ENm3jbwf+FtX7G2P6xxKEC7ntCcKpmUybD26mub2Zr5/5dYalDOt4mjDGuIMlCBfy+Dxx3wsilFMJYu3+taQlpTFtzDSmjZnG2qq1qGpU2zDG9J0lCJdpbW+l3lfvqi6m0UNH+8t+R3km09qqtZSNLiM9OZ2ZBTPZf2w/+47si2obxpi+swThMnW+OsAdayCCkockU5BZENXV1Psb/MlgZuFMAGYUzgBg3f51UWvDGNM/liBcps5bB7grQUCgqmsUu5iC4w0zC/wJojirmJOzT+6Y9mqMiT9LEC7TUajPBWU2QkV7sdy6qnWMGTaG8TnjO47NKJjBpoOb8LWdUJfRGBMHliBcpqMOkwufIOp99RxpPtLve7W0t7DhwAZmFsw8bjHgzIKZeNu8vPPJO/1uwxjTf5YgXCZYh8lNg9Twj5lMVUf7X/Z7a/VWGloamFEw47jj54w5h+QhyTYOYYxLWIJwmWAXU05aTpwjOV40y36v3b+WJEliesH0444PTRnK2aPOtvUQxriEJQiX8Xg9ZKVmkTIkJd6hHKcw07/XUzTGIdZVraM0vzRiMcIZhTPY7dnNocZD/W7HGNM/liBcxuN1V5mNoKzULHLTcvs9k8nj9VBeU94xrTVccFaTzWYyJv4sQbiMx+tx3QymoKLM/k91Xb9/PYp2JIJwk4dPJj8j38YhjHEBSxAu4/G5qw5TqKKs/k91Xbt/LTlpOUzNmxrxvIgwo2AG6w+sp629rV9tGWP6xxKEy7itUF+o4qxiDjQc6HPZb1Vl/f71TB87naQhSZ1eN7NgJvW+erbXbO9rqMaYKLAE4SKq6n+CcGsXU5a/7PfBYwf79Pldnl1UN1V32r0UdF7BeQjCmv1r+tSOMSY6LEG4yNGWo7S2t7r2CaIoM1DVtY/dTMFxhfD1D+GGpw9nSt4U1lXZOIQx8WQJwkWCdZjcOIsJ+l/2e23VWibmTmT0sNHdXjujYAZbD2+NysptY0zfWIJwkWCZDTftBRFq9NDRJA9J7tNiucaWRt459E633UtB5xeeT5u2sfHAxl63ZYyJDksQLuLWMhtBSUOSKMws7NMTxOZPNtPS3tLp+odwpSNLyUzJtPUQxsSRJQgX6ajk6tIxCOh7Vde1VWtJT0rn7NFn9+j6lCEpnDv2XNbud3aXOVXluZ3P9Xng3ZiBzNEEISKzRWSniFSIyOII59NEZFng/EYRGRc4nicifxKRBhH5hZMxuknwCcLVCaKP+0Ks27+OsjFlpCWl9fgzMwpmcPDYQT6o/6DX7fXU+v3r+f6G7/OTTT9xrA1jEpVjCUJEkoCHgIuBKcB1IjIl7LIFgEdVJwIPAg8EjnuB/wAWORWfG3m8HtKT0slIzoh3KJ0qyiziSPMR6n31Pf5M5dFK/+5xPRx/CAruNudk8b6l25YC8MaHb7Cvfp9j7RiTiJx8gpgGVKjqXlVtBp4F5oZdMxd4IvB6OXCRiIiqHlPVNfgTxaDh5lXUQR0zmXrRzdQxvbWH4w9BhZmFjMse51iC2HJoC5sObmJh6UJSk1J5dNujjrRjTKJyMkEUAqHTXSoDxyJeo6qtQD2Q19MGROQWEdksIpurq6v7GW78uXkVdVCw7HdvupnWVq1l7LCxjM8e3/3FYWYWzmTzwc14W6P/u8LSrUvJTcvl5tKbuXzS5byy5xUbizAmREIPUqvqr1W1TFXLRo4cGe9w+i0REkRH2e8eJoiW9hY2HtzIzMLjd4/rqRkFM/C1+aK+y9zO2p38pfIvfPX0rzI0ZSg3Tr0RgCfKn+j6g8YMIk4miCqgOOR9UeBYxGtEJBnIAWocjMnV3FxmIygzNZPhacN73MX0XvV7HGs51uvxh6Cy0WWkDEmJejfTI9seYWjyUK497VoACjILuOSUS1i+a3nHehRjBjsnE8QmYJKIjBeRVOBaYEXYNSuAGwKvrwTeUifnNLpcrbfW9U8Q0LuZTGur/LvHnTv23D61NTRlKGeNPiuq5b8/OvIRq/et5prTrjlu574FJQvwtfl4evvTUWtroGtua+ap7U+xy7Mr3qEYBziWIAJjCrcDq4EdwHOqWi4i94vInMBljwB5IlIBfAvomAorIvuAnwE3ikhlhBlQA4q31UtTa5NrF8mFKsoq6vFq6rX713LGyDPISs3qc3vnF5xPRV1F1MYHHit/jGRJZt6UeccdPyX3FD538ud49v1naWhuiEpbA9nO2p1c++q1/HjTj1n0l0W0tLXEOyQTZY6OQajqSlWdrKoTVPWHgWNLVHVF4LVXVa9S1YmqOk1V94Z8dpyqjlDVTFUtUtUBXfu5zlcH4PouJvBPdT147CAt7V3/g1DrrWVHzY5ui/N1Jzj7KRpPEZ8c+4SXK17my5O+TH5G/gnnF5Qu4GjLUZbtXNbvtgaqtvY2Ht/2ONe9eh01TTXcXHozH9R/wNM77MlroEnoQeqBpKMOU3pufAPpgeKsYtq0rdvf6IO7x51feH6/2puUO4lRGaOiUnbjye1P0q7tHYPS4abmTWVGwQye3P6kIzOnEt3+hv0sfH0hP337p1xQeAEvzn2Rb571TWYVz+Lhdx+2WWADjCUIl3B7HaZQPa3qum7/OnLTcjl9xOn9ak9EOK/gPDYc2NCvXebqvHX8ftfvuWT8JR3fQyQLSxdS663lpYqX+tzWQKOqvLLnFa5YcQXba7Zz/4z7+flnft7x9/Wec+6hXdv56eafxjlSE02WIFyiow5TgnQxQdeL5dq1nbVVazlv7Hld7h7XU+cXns+R5iNsq9nW53v89v3f0tTaxPyS+V1eVza6jDNHnslj2x7rthttMKjz1nHXX+7i22u+zeThk3l+zvN8edKXj5u2XJRVxILSBby27zWrwDuAWIJwiUSowxQ0augoUoakdDlQvcuzixpvTa9XT3dm+tjpCNLnbqZjLcd4ZsczfLb4s0wcPrHLa0WEhaUL2X9sP6s+WNWn9gaKNVVruHzF5fzp4z9x51l38ug/P9rp09f8kvkUZRbxo40/sgHrAcIShEt4vB6SJKlfs31ipSdlv4P/kPd3gDooNz2XkvySPq+HWL5rOUeaj7CwdGGPrr+w6EImD5/MI1sfoV3b+9RmImtqbeIHG37A1//4dXLScvjdF3/HgtIFXT4NpiWlsXjaYvbW7+WZHc/EMFrjFEsQLlHrrSU3LZchkhj/Swqzuk4Q6/avY/LwyYwaOipqbc4omMG2w9t6VSgQwNfm44nyJzh37LmUjizt0WeCTxF76/fyp4/+1JdwE9bW6q1c/crVLNu5jHlT5vHspc9y2ojTevTZfyr+J2YV+QesPzn2icORGqclxr9Gg0Cdry4hupeCutoXore7x/XUzMKZtGs7Gw5s6NXnXq54meqmam4uvblXn/v8yZ+nOKuYpVuXOronhVu0trfy8LsP87VVX6OptYmlX1jKv5/z770q0Q5w97S7aW1vtQHrAcAShEt4vJ6EmMEUVJxVzNHmoxF/m990cBOt7a1RG38IKs0vJSslq1fjEK3trTy27TFK80uZNmZar9pLHpLM/JL5bKvZ1uuklGg+PPIh81bN45dbfsns8bN5Ye4LfV79XpxVzILSBazat4q/HfhblCM1sWQJwiWCXUyJoquprmuq1pCRnMFZo86KapvJQ5KZXjC9V7vMrd63msqGShaWLuxTscA5E+YwKmMUS7cu7fVnE8Xzu57nqleu4sMjH/KTC3/Cf13wX2SnZvfrnvNL5lOYWegfsLaZYAnLEoRLJMJeEKGCU10/bjhxJtO6/esoG11GalJq1NudUTCDQ42H2FO3p9tr27WdpVuXMjF3IrOKZ/WpvdSkVOZNncffDv6Nd6vf7dM93OxX7/6K+9bfx5kjz+SFOS8we/zsqNw3PTmdxdMWs6d+D7/d8duo3NPEniUIF2htb+WI70hCdTF19gTx8dGP+ejoRx27wUVbcFyjJ7OZ/lr5VyrqKphfMr9fg/9XTb6KnLScAfUUoar8zzv/w0NbHmLOhDk8/LmHGT1sdFTbmFU8iwuLLuSXW37JocZDUb23iQ1LEC5Q76tH0YTqYhqWMowR6SNOSBDrqvz1kqI9QB00NnMs43PGd1uXSVX5zXu/oTCzkIvHX9yvNoemDOUrp3+FP3/8Z3Z7dvfrXm6gqvzs7Z/xm62/4YpJV/D9md+PymLGSBafs9gGrBOYJQgXSKQyG6EizWRau38thZmFnJx9smPtzizw7zLX1NrU6TWbDm7ivcPvMb9kPslDkvvd5vWnXc/Q5KE8su2Rft8rnlSVBzY9wOPlj3Ptqdey5Lwljk6tLs4uZn7pfFZ+sJJNBzc51o5xhiUIF+gos5FAYxBw4r4QLW0tbDywkRkFM/o0INxTMwtn0tzezNufvN3pNUu3LiU/I5+5E8O3Qe+bnLQcrj71alZ9sKrHpc7dpl3b+f6G7/PMjmeYN2Ue3z732zFZd7OgZIENWCcoSxAukEhlNkIVZRVx4NiBjh/6LdVbaGxtdKx7Kejs0WeTOiS10+mu2w5vY/2B9cybMq/Xc/i7Mm/KPJIkice2PRa1e8ZKW3sb3133XX6/6/csLF3IorJFjibxUOnJ6dxzzj1U1FXwux2/i0mbJjosQbhAR4JIgEJ9oYoyi2jXdg42+Es8r9u/jiRJYtrY3q036K2M5AzOHn12p+MQS7cuJSs1i6tPvTqq7Y4cOpLLJl7GSxUvUd1YHdV7O6m1vZXvrP0OL1W8xDc+9Q2++elvxiw5BM0qnsUFhRfwy3d/mVD/7QY7SxAuUOtLnL0gQgVnMgWnuq6tWsunRn4qJvWkZhbOZG/9Xg40HDju+J66Pbz50Ztcf9r1DEsZFvV2byq5iTZt48ntT0b93k5oaW/hnr/ew6t7X+WOs+7g62d+PebJAfylS+6ddi8tbS389O2BM2C97fA27v7r3TzwtwcGZGkRSxAu4PF6yErNImVISrxD6ZXirGLAP9W1pqmGHbU7HJveGq6z6a6PbnuUjOQMvnL6VxxptzirmIvHX8yynct6XRMq1prbmln050W8/uHrLCpb1ONChU4pzi7mppKbeHXvq2w+uDmusfTX1uqtfOOP3+C6V69jTeUafvf+77jkhUv44YYfDqhNkyxBuECdty7hZjDBP8p+Vx6t7OjucXr8IWhC7gRGDR11XDdTVUMVr+59lSsnX+noeM6CkgU0tTbx2/fduwDM1+bjzj/dyVsfv8W90+7lhqk3xDskwL+la8GwAn648YcJOWD9bvW73PrHW7l+5fW8d/g9vvnpb/L6la/zypdf4dIJl7J813IueeESfrDhBwMiUViCcIFaX2KV2QgaIkP8Zb8b/AlieNpwTs/r3+5xPSUizCyYyYb9G2htbwXg8W2PIyLMmzLP0bYnDZ/ErOJZPLPjGRpbGh1tqy+aWpv41zf/lTVVa1hy3hKuP/36eIfUISM5g3um+Qesn33/2XiH02NbDm3h1jdu5asrv0r54XLuOOsOVl+xmpvPuJnM1EyKs4r53ozv8YfL/8CcCXN4ftfzHYkivBs0kViCcAGPN7HKbIQqyirioyMfsW7/OqYXTI9pufIZhTM42nKUbYe3cbjpMC/sfoG5E+YyZtgYx9teWLqQel89y3ctd7yt3mhsaeS2N29jw4EN3D/zfq6afFW8QzrBZ4o/w/mF5/PLLe4fsP77ob9zy+u38LVVX2N7zXb+7ex/Y/UVq1lYujDiGFdhZiH3zbiPVy9/lcsmXsbzu5/nkhcv4f7197O/YX8cvoP+sQThAolWyTVUUWYROz07qfXWcn7h+TFt+7yx5zFEhrCmag1PbX+KVm3lppKbYtL2p0Z+imljpvFE+RM0tzXHpM3uNDQ38PU/fp13PnmH/7zgP7ls4mXxDimi4IC1r83Hz97+WbzDiejtT95m4esLmbdqHjs9O/nW2d/itSteY37JfIamDO328wWZBSw5bwkrv7ySyydezosVL/LFF7/I99Z/j6qGqhh8B9HhaIIQkdkislNEKkRkcYTzaSKyLHB+o4iMCzl3b+D4ThH5ZyfjjCdV9RfqS7AprkGh209Ga/e4nspJy6Ekv4Q3P3qTZTuX8YWTv+DoCu5wC0sXcqjpECv2rIhZm5050nyEf3njX3iv+j0euPABvnjKF+MdUpdOyj6Jm0pu4g97/9DlgsdY23xwMwtXL+TG125kt2c3i8oWseryVdxUclOPEkO4sZlj+Y/z/oNVl6/iiklX8HLFy1z6wqXct+6+LjfccgvHEoSIJAEPARcDU4DrRGRK2GULAI+qTgQeBB4IfHYKcC0wFZgN/DJwvwGnoaWB1vbWhO1iCs5kOnX4qeRn5Me8/ZkFM6moq+BYy7GYz9KZPnY6JXklPLrt0Y5xkHio99Vz8+s3s712Oz+d9VP+eVxi/D61sHRhx4B1PP/7gb80y/zV87lp9U1U1FWwqGwRr13xGjdMvaFPiSHcmGFj+D/T/w8rL1/JlZOvZMWeFXzpxS/x3XXfdfXK/P4XqencNKBCVfcCiMizwFxge8g1c4H7Aq+XA78Q/yTtucCzquoDPhCRisD91kc7yF2eXdz9l7ujfdsea1X/D0aiJojgE0S0NwfqqRkFM3j43Ye5sOhCTh1xakzbFhEWnrGQO/90J3NemkPqkOiXN+8Jj89DQ3MD//2Z/+bCogvjEkNfZCRncPc5d3Pnn+/kSy9+Kaqr3nujub2Zj49+TH5GPnefczdXTr6SjOQMR9oaM2wM35n+HRaWLuTRbY+yfNdyXq54mZOzT0bo+/qU8wvPZ9E5i6IYqZ+TCaIQCE2NlUD4FlUd16hqq4jUA3mB4xvCPlsY3oCI3ALcAnDSSSf1Kcj0pHROyT2lT5+NltL8Us4be15cY+irCTkTmF8yn2tOvSYu7ZfmlzK/ZH7Uai711meKP8O8KfM4cCx+M1WSJImrT72ac8acE7cY+uqzJ32WO8+6k/Ka8rjG8ZXTv8IVk64gPTk9Ju2NHjaae8+9lwWlC3h6x9P97m6K5t7vocSpvXZF5EpgtqouDLz/GnCuqt4ecs22wDWVgfd78CeR+4ANqvp04PgjwCpV7XTKSFlZmW7enNiLb4wxJtZE5G1VLYt0zslB6iqgOOR9UeBYxGtEJBnIAWp6+FljjDEOcjJBbAImich4EUnFP+gcPt1jBRBc4nkl8Jb6H2lWANcGZjmNByYBtvu5McbEkGNjEIExhduB1UAS8KiqlovI/cBmVV0BPAI8FRiErsWfRAhc9xz+Ae1W4DZVbXMqVmOMMSdybAwi1mwMwhhjei9eYxDGGGMSmCUIY4wxEVmCMMYYE5ElCGOMMRENmEFqEakGPuziknzgcIzC6Y9EiRMsVickSpxgsTol1rGerKojI50YMAmiOyKyubORejdJlDjBYnVCosQJFqtT3BSrdTEZY4yJyBKEMcaYiAZTgvh1vAPooUSJEyxWJyRKnGCxOsU1sQ6aMQhjjDG9M5ieIIwxxvSCJQhjjDERDfgEISKzRWSniFSIyOJ4x9MZESkWkT+JyHYRKReRO+IdU1dEJElE/i4if4h3LF0RkVwRWS4i74vIDhFx7dZ9IvJvgf/320TkdyISm+3NekBEHhWRQ4FNvoLHRojIGyKyO/CnK/bN7STWnwT+DrwnIi+KSG4cQwzGdEKcIefuEhEVkdhv9B5iQCcIEUkCHgIuBqYA14nIlPhG1alW4C5VnQJMB25zcawAdwA74h1ED/w38JqqngZ8CpfGLCKFwDeBMlUtwV8i/9r4RnWcx4HZYccWA2+q6iTgzcB7N3icE2N9AyhR1TOAXcC9sQ4qgsc5MU5EpBj4AvBRrAMKN6ATBDANqFDVvaraDDwLxGfz4m6o6gFVfSfw+ij+f8hO2IfbDUSkCPgisDTesXRFRHKAC/HvO4KqNqtqXVyD6loykBHYXXEosD/O8XRQ1b/i37Ml1FzgicDrJ4DLYhlTZyLFqqqvq2pr4O0G/LtUxlUn/00BHgTuBuI+g2igJ4hC4OOQ95W49B/dUCIyDvg0sDHOoXTm5/j/ArfHOY7ujAeqgccC3WFLRWRYvIOKRFWrgP+L/7fGA0C9qr4e36i6NVpVDwReHwRGxzOYXpgPrIp3EJGIyFygSlXfjXcsMPATRMIRkUzgeeBOVT0S73jCicilwCFVfTvesfRAMnAW8LCqfho4hnu6QY4T6L+fiz+pFQDDROSr8Y2q5wJbBcf9N97uiMh38HfnPhPvWMKJyFDg28CSeMcSNNATRBVQHPK+KHDMlUQkBX9yeEZVX4h3PJ2YCcwRkX34u+w+KyJPxzekTlUClaoafBJbjj9huNHngA9UtVpVW4AXgBlxjqk7n4jIWIDAn4fiHE+XRORG4FLgK+rOBWAT8P+C8G7g56sIeEdExsQroIGeIDYBk0RkvIik4h/0WxHnmCISEcHfV75DVX8W73g6o6r3qmqRqo7D/9/zLVV15W+6qnoQ+FhETg0cugj/Pudu9BEwXUSGBv4uXIRLB9RDrABuCLy+AXg5jrF0SURm4+8WnaOqjfGOJxJV3aqqo1R1XODnqxI4K/D3OC4GdIIIDErdDqzG/8P2nKqWxzeqTs0Evob/N/Itga9L4h3UAPCvwDMi8h5wJvCj+IYTWeApZznwDrAV/8+me0ouiPwOWA+cKiKVIrIA+C/g8yKyG/8T0H/FM8agTmL9BZAFvBH42fpVXIOk0zhdxUptGGOMiWhAP0EYY4zpO0sQxhhjIrIEYYwxJiJLEMYYYyKyBGGMMSYiSxDGGGMisgRhTAyJn/3cmYRgf1GNcZiIjAvsSfIksI3jy78Y41q2UM4YhwWq8+4FZqjqhjiHY0yP2ROEMbHxoSUHk2gsQRgTG8fiHYAxvWUJwhhjTESWIIwxxkRkg9TGGGMisicIY4wxEVmCMMYYE5ElCGOMMRFZgjDGGBORJQhjjDERWYIwxhgTkSUIY4wxEf1/tg5NMJRacMUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot C-H rdf\n", "chrdf = model.avg_chrdf.result().numpy()\n", "plt.figure()\n", "plt.plot(chrdf[1, :], chrdf[0, :], label='C-H', color='C1')\n", "plt.title(r'C-H RDF')\n", "plt.xlabel(r'r')\n", "plt.ylabel(r'g(r)')\n", "plt.legend()\n", "plt.show()\n", "\n", "# Plot O-H rdf\n", "ohrdf = model.avg_ohrdf.result().numpy()\n", "\n", "plt.figure()\n", "plt.plot(ohrdf[1, :], ohrdf[0, :], label='O-H', color='C2')\n", "plt.title(r'O-H RDF')\n", "plt.xlabel(r'r')\n", "plt.ylabel(r'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.7" } }, "nbformat": 4, "nbformat_minor": 2 }