{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analyzing a lammps trajectory\n",
    "\n",
    "In this example, a lammps trajectory in dump-text format will be read in, and Steinhardt's parameters will be calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyscal as pc\n",
    "import os\n",
    "import pyscal.traj_process as ptp\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we will use the `split_trajectory` method from `pyscal.traj_process` module to help split the trajectory into individual snapshots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "trajfile = \"traj.light\"\n",
    "files = ptp.split_trajectory(trajfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`files` contain the individual time slices from the trajectory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'traj.light.snap.0.dat'"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "files[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can make a small function which reads a single configuration and calculates $q_6$ values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_q6(file, format=\"lammps-dump\"):\n",
    "    sys = pc.System()\n",
    "    sys.read_inputfile(file, format=format)\n",
    "    sys.find_neighbors(method=\"cutoff\", cutoff=0)\n",
    "    sys.calculate_q(6)\n",
    "    q6 = sys.get_qvals(6)\n",
    "    return q6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a couple of things of interest in the above function. The `find_neighbors` method finds the neighbors of the individual atoms. Here, an adaptive method is used, but, one can also use a fixed cutoff or Voronoi tessellation. Also only the unaveraged $q_6$ values are calculated above. The averaged ones can be calculate using the `averaged=True` keyword in both `calculate_q` and `get_qvals` method. Now we can simply call the function for each file.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "q6s = [calculate_q6(file) for file in files]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now visualise the calculated values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f5e424c0f60>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAy0klEQVR4nO2dd5gUVdbG3zMzDDkzCMwAQ5ScQXIQFBAU4xrWLIt5ddfPFYyYUXfVNSKGNcMaUFgBURFBJA5IFJABBhhAGHJm0v3+6O6hprvCreqqrjDn9zw8TFdXOLf61ls3nHMuCSHAMAzD+J8ktw1gGIZh7IEFnWEYJiCwoDMMwwQEFnSGYZiAwILOMAwTEFLcunCdOnVEZmamW5dnGIbxJcuXL98nhEhT+841Qc/MzERWVpZbl2cYhvElRLRN6zsecmEYhgkILOgMwzABgQWdYRgmILCgMwzDBAQWdIZhmIDAgs4wDBMQWNAZhmECAgs6w5QBhBBYvu0g/JAue/7vedi+/4TbZvgSFnSGsZG5G/Zi5prdbpsRwzerd+OyNxfiyxU73TbFkOvfW4r+L8x12wxf4lqkKMMEkZveXwYAyJkwwmVLSrNt/3EAwNZ9x1y2hHESbqEzDMMEBBZ0hmGYgMCCHjAWZu/D9FW73DaDYRgX4DH0gHHNO0sAABd1bOCyJQzDJBpuoTNMGcIHXotMHLCglwH2HzuNcVPX4HRhkdumMIwnOFVQhHd+3oKi4mC94VjQywDPztqAyUu343+rvOcfzSQWIrct8AavzNmEp2asx9QVuW6bYitSgk5Ew4hoIxFlE9FYjX0GEtFKIlpHRPPsNZOJh2LuZ/uK46cLcex0oSPn5qoQ4uip0P09WRCsXquhoBNRMoDXAQwH0AbA1UTUJmqfGgDeAHCREKItgCvsN5Vh/IcQApPmb8beI6ekj2k3fjbaPTbbQauYoCLTQu8BIFsIsUUIkQ9gCoBRUftcA2CqEGI7AAgh9tprJsP4k017j+GZmRtw56crpI9RtqLPfngWhr083zZ7tIZcCouKceRUgW3XYdxBRtDTAexQfM4Nb1PSEkBNIvqJiJYT0fVqJyKiMUSURURZeXl51ixmGB9RWBRS50gX3yynC4ux4Y+jttmjNeTyjy9Wo8P470o+n8wv8kROmj8On3I0UVfQhqBkBF3tnR59G1IAdAUwAsBQAI8QUcuYg4SYJIToJoTolpaWZtpYJj78kGlPjeveXYJz//WT22b4GjKYDZ36a+mkXeOnr8Mdn6zAr9sPOmmWIT2fnRN3oq4Dx/OxfFvpctg1Odz96R/wwBerUVQssCzngD0njQMZQc8F0FDxOQNAdChiLoBvhRDHhRD7AMwH0NEeE5myzs+b9mFL3nG3zfA1si/zwycLsPPQSew8dBIAHJucTSSXT1yIy95cWGqbXW2bvKOn8d+sHXjzp2xcMXERFm3eb8+JLSIj6MsAtCCiJkSUCuAqANOj9pkGoB8RpRBRJQDnAFhvr6nBZnPeMWSOneHoW96oleYE2/efwIHj+Qm/LmONwf+ahz4TfnTbDFvRawzY9Uhs2hvKYrnHxOS3ExgKuhCiEMBdAGYjJNKfCSHWEdFtRHRbeJ/1AL4FsBrAUgDvCCHWOmd28Pglex8AYPpK+/OwTHUxB3b/F+ai//NlN7e1V/y+ZV/m+46ddtgSxkmkcrkIIWYCmBm1bWLU5xcAvGCfaUwieXTaWtSoWA5/P/9s288dhG67VbwybWF1/sQr9sty/HQh9hw5haZpVaT291v5jOBI0ThYtHk/doXHGv2A3kP94aJteOXH7ARawzD2c9N/luHcfxnHNdrdc/JKCgEW9Di4+u3FGPKi/4Ni846e6Wb3fGaOr15SWvy6/SAem7bWt549dmN1/sQrQ0ayLHXJ0+Sb1SEXT7fvFwu6BXYcOIEb/7MUAHAi3/+hw92f/qHk7z+OnMK0OMfxez87B89/uyFes+LiiomL8MGibSh0qeVUWFQMwP0HPF6U78MT+YU4XoaHz2Rwu/3Agm6BZ2etx08bEx8Y9cjXa9H3Oec9EERMmIE5dh0+hTd+2myTNf7kb5+tctuEuFB7EbV9bDba6qQkOHg8H/dO+dUzcyYb/ziKtTsPq37ntvA6BQu6j/ho8TbkHjyJlTsOWTreDbdFt4l+cAuKivHHYeddy/7nwKpR+YXFuO+zVa4NiRmJ4Otzs/H1yl2YvGS79DnfXbAVbzr08h/68nyMfHWBI+f2KizoFiDV4NnEcfHrv1g6bsKsDfg8a4fxjgFm3NQ16PnsHJzIT2wr0o4W4U8b9+LLFbl4dNo63f0Wbt6XcH/oLk9+j3cWbDV93JPf/IbnHBieM/p9g9q28bWgr915GOt3H7HlXEIIfLfuDxTbMOZ6qqCoZAw1wsHj+ThVUIRfsvdh9AdZpSbrxn652vAhjcZKqPG+Y6dx/xerTR3jd6If3B/W7wEAnC4oVtnbPRZu3ocfN+yx5VzXvL3Etpbp9e8txSqJHqETwWPZ4WAdK5zzzBwbLZHH7ReFrwV95KsLMPzfP1s6NnvvUfR/fm5JRfzq150Y89FyfLgoB0Cogv6+x1pSpFaPfIsbwpOmETo/+T2umLgIoz/Iwg/r95TKwzxlmflW8ytzNuGKiYssRZZu228cRj/4Xz/h2Vnmgn3zjp7Gl8u9tWCAVss40UOoRg/6NW8vwc3vZ9l2PaXnUrws3lI6nD0eoTWD0oOs9SPfmnINtJoMze/4WtDjYeK8Ldh+4ERJi23PkdADsDvcVR328nyc/5L1tKW/ZMfmdFiz83CJkGsLjcB9n63Cki36OSE27Q29bPYeMf/gDnjhJ8N9Nucdx1vztpg6718+zMJ9n2tPBgoh8O6Crapi89a8zaZfIHpoCaibDahbP8rCN6vtH1tPNG6U4WRBkS+WUHR7sjWwgr7r0EmpccR5G/OQOXYGNoVb4weO5ePQiXzs1WvhOKgK+YXF+HJFLq57d6nufm5XHDX0FnG4//NVaDJuJp785jfcPTk2N/izszaYfoGokb33GFbuOGT6/mzdd1zTI8IOhABmr9uDuz791bFrWCXSK/U6XqzzMgz650947cdNCbmWVOi/H+kdTjCUM2GE6veRyvHtuj8AACvCaUI/X56Lz+MYNohO06kFUUh84l0ezokxOyceHOU9PXLSue5wpJuekmTuxgz6508AtOuLVZwZU7X3BzKav/GpjkoxY81u3NA707bzqf3eW/cdxz+/+x2XdslAgxoVbbuWGr5soX+xPBddnvxeat/x0/Ura+T+21Vpo9N06jHkxXkxwzolYmogBHaI7vHThej4+HfGO0oia5JV0/86+Vdc8oach0/kwSoo8tbkpx8RIjTnFEHWy8vuRF9OvFiWbk1cZGnvBGSx9KWgP/TVGulZ9fcX5qhujzd4RoaComKccngR2js+WaHrmZOzT3sCNHvvMRw+Gbvs2Oxwr8WItTsPY6vi/E53iaev2oVftx8ydUzbx2Zjc96ZSbyDJ0Ll3apzX9xEaatW3dkWXsFnTe5hjHptAU5qRCvPWL07JvWB7E/086Z9JX9/sDAHQ1480/CQeXb+t2oXuj31A7IcCMUXQmDy0u1o/9hsU15p/Z+fa9jA8zu+FHQ7ibTknOgZX/jqArR65FvV7+wUv9OFsa3QvUdD49kDw0MJZlidqz6WXFQscPP7y0paNSNfXVAyVGGGyL3+eVNo/mLHAeeWGAOADbtjvZXmrLfHRdBuBisSS2mlHd609xiOnCrAE9+sw6rcw1ijMfZ/56crMHdj/Mv7/qExNzJt5U5kjp2B3IOxv1+kjvym4lace/BEXBOcd3yyAuOmrsHR04UoKJbvgW0/cEKzgRfNweP5+G1Xadt3Hz6JZ2aut8W12SkCJ+iTl25H5tgZtp3v503WQ/ztXAtSD7UW04/rjR9ks9Vy/7HT+HHDXs0Fj2V7PZG9vgiPq2dtc7bbKyDwWdYOTFtpnBc+nt/bbvQm5juM/w4bJerXvmP52HHgBDaHXQ3taLhEhly+Di9bZ8a9N7+wGH2fm4u//Xel6etGehuz1v6h2BZyq7R7Qvui1xfggldKu0T/7b8rMWn+Fvy6w91l+fTwpaDrTTSNm7rG4jnVT3rnJ7HipXV5GcGI8OGibarbJYfQS4mnWms/kW2IyINmttdhZ69o056juikR/vHFatwzZWXJZy1TjbyLvMQRCV/rDxfloN/zc0vWDLWjXuQXFemeS+/5LAy3qOduMP/ivOk/y1S3XzVpselAKqO6t+NAbHqFgiJr9TyR+FLQbSHqR7Ejzep/TQQImRlbLC4WeOgray+qeNm2/zh2HjqJqyctxpFTsePtAPC1iReZEjufi/Nemm8pJcL+Y6cdybviFdbutCeSWsnrczdj7oa9JcIWPUkqhHZvLZ7HLEvSg0wGNTMmzNqgGv8xeel23bkoNdbkHkbm2BkJz7sTWLdFWUKV0R5pscNFreQhUZwrO+8YPolKeGT0YMg8OC9+/7vhPsogpBmrQ13dvKOnMXfDmSGdfUdDE9QebrjEELk/Yz5aLu1qGtf1bLg7Wr1IMxGUdvWK5v2uaGFLnnTv0VPo8bQ7IfkyTJy3GRPnbca6x4eWbBNCYNzUNahRqRyaSa6CBACfLAn1wEvdpwRQZlvoJY+Ajf1+M0m74n28NynCr9XOJSMg801WNuU535VIxPT2fPVAoUhvKFFRm2rBPJGy2BXGfqqgCN+peAfZmchNqxc54dsNKCwqxoHj+VijMaFdcg7brDnD4ROxPTdluef/noesnAPYvt/ZyW+7uGdKbH05pCijzD10a1imzLfQreBkGtrIQ1tYpF8jSrsLxu7rdIVSe2FEX/PpmfaF8itZvu0AujauVfJZJnmUFmpum2p8sDAHs9buxpQxvVS/f/Kb32J6UVoUFhVj16FTEBB4f2EOruzeUNpeNVbtOITmD81Co1qVsN1hj6EIykfgXoMJzuvfC81LfH6b+r2LMPzfPyfMg0TvCV6n8G7x8ni5Gr4U9HhaPatzD6Fdg+qx50xwmjStq0UmryIr7cz7PQ/rdllrdW3S8T6wEnCj9ZI5c+vkan+094/WQ/PzpjykJMV2Ii97cxEA4NIu6XjxT50wymI6YS3yC4uRmlL6uo8Z+C9rCanai++F2RvxlqL38p9fcswbacIGp9CbFLUyxGRX5lSZ60xdcWbe51RBUanPStRKsW7nYXTMqBFTR7yALwXdLEXFAp9n7UCTOpVx5aTFuH/omZXt88M+3PvjjGo7mV9kKsf2nA1y/sE3vGfsdRGK5DuGDX+ceSDemr8ZD3+9VvOYl38wHj+P5rW59i0ivXbnYcPWsZHHydQVO/HinzpZun6+iu9+hH9+txEPXtDa0nm1+H3PmaGdBdn7dPb0D35drzU6Q+vLP2zCxHn6i2wo23vj//cbNucdx5MXtzO8VqKTwZUJQf90yTY8Mm0dOmaEWuYvzN6ISzqnl9pHywVMxjUMAHo8/QOOaiy99et25yfdoherVnO7UpJ70P7Z933H5HNi25Wv+8q3FqluLzAYsvrPLzl47MK2qt+Z9WgAtHt4ar1Jq1GqXlpxavbaP9CsrvokoRDq5VYLQLITqy8YZVoDrXMKUTpP06rcQ5au5TTe6zM4QGRC46ByYsPEj19ULPDBwpySVp3aY6Ul5gBwyRvy+V2UqE02qeLPhlIJCzfrpwrWY0kCc3HYhZWFxXcfdmfZOS12WVjG72//9eY6qyctpOcwEyeSSMqEoCeFM+9ZzWw4eel2PDZ9Hd7+Wd9rw26ufnux1H5WKs+0lfH5Xqvle7fKFx5bFMOLHD9d6KshjkidtHvB6OhYg48XqwfoxYNSrDXvuId6S0oCI+inC4s0u06Re698HnaYGHKIBNRoBdY0GTdT+lxmUMuDoYaPnnNDJi+VX2DYSRZt3o93NF7gQCiTYObYGfj+t8TkhLn1o+UJuY4dKLXu+W83Sh0j+7KKjgZ+akb8nlQFheYTmJ2w+UVlF74UdLWX46NfryuVEU5JEsW20M0Ek0QOm732DyzZst9zL+cA6TnGTV2DTXuOxj1JHS9Xv71YVywiHhkfGCR7+nKFPb2PzXneyw4ZlIbEUp2oba0ybpKMX0i0VgRiUvTwiQLM0Vlg95U5odVC4q2AOftP4MpJi1G3avn4ThRA9FqzZjnvpfmoUamcbedzksjQgtpzu/vwSUzSCK4KAnreOmZdi7UyfJrBrheM0vKWD8/S3O+4B1vpgRD0jk/oL9JgZRJKD93l6VzA7bHVrfuOSwfVyHJIdkLYAbSG1pTICJZRcJhZ7PZyOXyyAEclymoWvVwuWtgdSxAPspZ3e+oHw33yba4DRgRC0GWxOinqtmAa4XZ65uhV4f2C1u+6eMuZLviOAyfQsFYlw3M5IYxOc96L8zzXOPETMt4xj+jEgjiBL8fQrWJV94x8mt3GLRepCF4c35VBZkJNy4ti6dbQSywyVGC3N4dTKBs1Tom51+aYzBKP+R8v3obDJwo0Yw1e+v53bNvv3PPiS0G3esOtHldoYlUUN/C4eZ5FJsGY1qvylR9DUbNHw4FniejE2dFTLCoSyN57FIdOyAeBWcHOpGSy/K4TIJQo1u06gpGv/YxlOepOF/+eswk3va+e190OytSQS1Ap8viQkJ+REdHdh09iWwLyqIyxwXVxzc7DeOfFrahXrYINFmlj1Gs0k/JXlmEv/2y8k8MIYRylrZd2Il5Y0CXwul6+/4txS5OxRvRvf/RUQamVjwCg17Pqq7l7ceghElmrtU6oXRw1SJmRbyE5XKKwEgUbwe3hT18OuSQaj+s5pgd4xR2v8fXKXfhRMrFaWUWI+CORg4wTeZQiSAk6EQ0joo1ElE1EY1W+H0hEh4loZfjfo/abGj9WJ4G83kJ3Gy+2RO2Cf3rGCdQWQ7EDwyEXIkoG8DqA8wDkAlhGRNOFEL9F7fqzEGKkAzaWorhY4LjNfuVGuN2NMsLJMTkZgvzCiy5bot3Q/IhssiutFa3KAtscWr1JpoXeA0C2EGKLECIfwBQAoxyxRoKPl9ifjMcQjwvWQReDcIKOgLA8geeldLeJRGatWsC5Fa3cxO3GjYygpwNQLmefG94WTS8iWkVEs4hINdE0EY0hoiwiysrLs7Z4qpPjTwwTjRDAR4ty3DaDYaSQEXS1Zkb0e2gFgMZCiI4AXgXwtdqJhBCThBDdhBDd0tLSTBmqZ4zTWI0wZfxPQVExDnAPiJFEVimc6rzJCHouAOUqthkASk1hCyGOCCGOhf+eCaAcEdWxzUqXcTu0nnEPK2uvRrjvs5X2GcL4Atn5LKfaiDKCvgxACyJqQkSpAK4CMF25AxHVo/CAIRH1CJ/Xnwk+VJCJKGSCSTwvc2VOGIZJBIZeLkKIQiK6C8BsAMkA3hNCrCOi28LfTwRwOYDbiagQwEkAVwmvZ7RiGAm4FjN+QipSNDyMMjNq20TF368BeM1e0zQom44DjEsICHyRtcN4R4YxgVNj6Bz6zzA6TF2x020TGEYaDv1nGIYJCCzoDMMwAYEFnWEYJiD4TtDdSJzPMAzjB3wn6AzDMIw6LOgMwzABwXeCXkYT2DEMwxjiO0FnGIZh1PGdoJ84rb9WIcMwTFnFd4JewKkPGYbxOTsPObOug+8EnZMlMQzjdz5e7MzKa74TdM+vB8cwDGOAU8sT+k7Qi91dD5lhGCZunHLW85+g85gLwzCMKr4T9CIWdIZhfI6ba4p6CtZzhmH8jlM5qXwn6AzDMH6HW+hhOPSfYRhGHd8JekoSKzrDMP6GvVzC9Gxa220TGIZh4oL90MNUKc/rWjMMw6jhO0FnGIbxOzzkwjAMExTYyyUEu6EzDON3uIXOMAwTEHhSlGEYhtHFd4LOof8Mw/gdjhRlGIZhdGFBZxiGSTA8KcowDMPo4kNB50F0hmH8DXu5MAzDBAQecmEYhmF0kRJ0IhpGRBuJKJuIxurs152IiojocvtMLA27LTIM43ecWkrTUNCJKBnA6wCGA2gD4GoiaqOx33MAZtttJMMwTJAoLnZJ0AH0AJAthNgihMgHMAXAKJX97gbwJYC9NtrHMAwTOJwaaZAR9HQAOxSfc8PbSiCidACXAJhon2nq8IgLw5Rd2qVXc9sEWzh6utCR88oIutqEbLSuvgzgASFEke6JiMYQURYRZeXl5UmayDAMEyIlif049JC5O7kAGio+ZwDYFbVPNwBTiCgHwOUA3iCii6NPJISYJIToJoTolpaWZs1ihmHKLL2b8RKUesgI+jIALYioCRGlArgKwHTlDkKIJkKITCFEJoAvANwhhPjabmND13LirAzD+IE7BzV32wRbqFMl1ZHzGi7QKYQoJKK7EPJeSQbwnhBiHRHdFv7e8XFzhmEYAEhyKk1hgnHIycVY0AFACDETwMyobapCLoS4MX6zGIZhYgmInqPIRbdFT1GrsjNdFYZhmETh1FyA7wS9F0+KMAzjc2pUKufIeX0n6AzDlF2CMuTiFCzoDMMwCYfT5zKMJ+nYsIbbJjAMAJ8Kep0q5d02gWFK+O+Ynm6bUGYox5Giuvjy7giOLmI8RIVyyW6bwDAAfCroDOMVrjmnkdsmlCl4UlQfXwq6Xvt87eNDDY/nSsHYQY1K5fDMJe3dNoPxIU5pkD8FXWfIpUp5qeBXVSZe29XysU5wQ6/GcZ/jgWGtbLCEMaJhrYpum1AmcGpx5aDgS0Ef2rZeXMc3qK7+8J3bqm5c57UbO6KD3aj/ZXGScHi7+m6bwNjIiPbO/p68SLSCAS3jS71bWyPTmfDY8hl22OPH9ky/FnV0v0+v4b3WMPeEgoVf50Z8KejxtjqTk/whc3bl76kaxzBUxXLJSEnw/fJLt1pppV/qFBNsfCro8T08iRYoq9jlnfnCFR0tH/vuDd3sMcIERm6p9atXSJAl+nirP8cAQPfMmpaPbdvgzPJ29TxSx8ziT0G3cMy/r+pU8rdaTmWviIQSoyT4PTJrSZ0nnoni3s31hz/UsFPo1N69FVPZ77ss0qVRDcN9alayno11imLup1laFcz9v4G6+zeqVcnytdjLRYGVJPfKRp9fBOGWvk30d/BHRyMualXmqGC/cVY1936zeISyaoXSGRCb1KkcpzWJx5+CHofVlVITPyZsFbvGZWUq+b1DWsR1vFNMvb23exdPEBMuDZYvu1OB3FpDrco5orIeRO5LQSeNpun5bc4CADx2YZuY74w8RhKhWWa7aJEKXEmjRyFjM5FcJb93SEsTlskTLVaPjGyDs8+qqrrvbQOaxWxrVDv2nvll0lSWFmdVcSw/thOMV3m+lDilqVpzK2skggmdQK0aVqtgfXjTDnwp6FpKFrnBqSmxxfLCm7ulhpBZ5cruDVW3u90DaZoW6qr+64qO6NGk9Dj/sHb1UKGcerXT2r7+iWG4d0gLdJLMavjEqLZ4+cpO0vZaRatO3dg70/Fru8mQcMNJDS+6lCaK7//WHwvGniu1r1ajNF58KehaY+heEG09ejaVm8SMRq20E6/tgku7ZKj2RmKOT7C+161aATkTRuCyrhm2nK9iajLuHdJSd1gomlGdGuCL23rZcn0ttO5rVZdbaWq8cHkHLHtoCF66MuTx1CGjuiPXEUK4+hwmsq5HX6p6pXKoksotdNNk1CzdCri0c7rhMV4Xe6vc1Ed/4tSplkA0aVXlJsKEELiok/7vpfVbmfkJiQjdVLyAruiagV/GnotHRxq/CBOJ07+SQOg3yqwd6j3FM2xV2WXRskJdyfoZL0ku9459KejRLobnS6QCiIiB1u3WquALHhhkwjJvoBS+iqnJ6NpY3Tc3U2V8GgAm/8V86L6ZF+bNfTJVt8u+fIz2urBDA83vXriiI9JrVEQbhc+xHnPuG6D5nX2NBPLNvEBqShJqGi7U7kzrSe0eDWktl65jcGvtYSI9Pri5Bzpq9Gai7UlU40kPXwp67I2LrwINaJmGlxV+6hE+vuUcZNS07msazcUSPQkl+gE28pVHK1+3lo+5nQtxq5VAS7wMUx1ofB098aolOMoEWuc0qRUzvq9Gs7QqhvvEjz+6j09f0g4z/9rXcD+nesNqz0NlyRgLq3maBrRMw7S7jMsMmBvuYT90E/RrHpvrpVJqMga3qou3rouNfPzg5h7ortI9tyuD3sanhiFnwgjLKy2pCWD/luYDfmLOK71f7J6rHjs/7utbJfp2NKtrXnSJSGqozgopXlxVR0Jkf3tC31vkz+c0RvO6+hP7aVXLa15q7HBr+W70Jlpl67Adz4sR7rfPAyPopW9lo9qVkDNhRKltSQS8e2N39DVI/GQ3zdIqo3yKvYFMT4xqi0o2jGPKthKiW84Tr+2K6hXtd7OLvDi0WupeS56mxV/6N8Foo6CwUnhBCmBLndIL2lO6pfZpLt8LjNRTIzde3XM4cI/HX9S29DXC9snEj3C2xTgZUmoM7czt9JuL2ax7+uG6nvHnSQesV/Jh7eJLX6yFkWBHBMfqBFd0eZ3q9lZKTcHDJiddnZ8UdfdluPzhIaU+16goH6If/TvN+Gs/9NUYLtSq00781gNappV6SUUuoeflFHl2qzjkCRVoQd/09PCSv1OSlUU9U7mVvs2pyd66HWqPYOv61UxNoNlRkb0w2QOExr2fv7wDHruwrfHOKsjci3bpsZOlWjlzZNe2vbKberyAVS7pnI7alVNL+e3fPrAZPrqlh6nzfHtvP818JZ+MPseSbX2b11F94daOGm6MJyajSZ3KuKyr+nCZ2ovLDvfVOfcNwH9u6q57PbX61TQqfcAlXdLx8IjWuPtceRdcM3hLwSwSCaSJDigqZyDQtw5oilGdznhEEFGpVYKcGFaQRVkRzMqpUmiiNefhEa3jsCp+9DTQ6MVBRPhTt4bSE2Gx5ze+3jd394sVJMkfYIiGJ8U/hp1teGyr+iGBu3VAU13PGgB46cpOWP7IeVjwwJkglgeGtdLMPTKsbT2MVPH8aVWvmuYxvS1OjN8xsDn+7/xQefs2r6M5PHLrgKaq26MFUIlevdbj7HpV426SNEurgkFn60+sRuqT8lrRc3PJRBjdr6ljC4v7z6EUsW/CQa3q4vaBzfCXfuqVROUMAICujWrGuh4pPteII3ObDE3rVMaWfcdVv5twWQfL5z23VV38sH6v6nej+zXFUzPWA7C7G6ox7m2hp2/2GMteFQbl/+xW/ZZddN1567quKCwuNm8GARXC8yzdG9eS9qyRnWSfeJ380ooTr+2C+horesmQlESlGlYLx56LE/lFJZ9/GXsuDhzL16x7D49sjZvfz5K+nl5vtVJqMk7kF4EosW6hetXR6YGvQLTQk5MIDwxrhVqG/rHWeOXqzo6cFwiNB/dqWrskgs8Oejatjat7hFZcMVOPG2v4pZdVjO5H9JBLchKpToBH9qpaIQXt09V9miOT9W7/BsPa1UfHqBQLTepUNmVXScwHhRpFDRReKuk1KqJ9RnXtsW6dt6yaKEf/BsrjX7qyE9rUr4ZKEq1hrd6VGUhFTaOHgJx+rQRC0J3gZkUEZsRtqn16dcuuV1osfWgIJo/piUs6x4bJOx3dGl25Pr+tl+lxWO2zybFQkfsiuvIb5YMHgE//co70S0u2lWb3bY/8jqnJSXhVo3FwY+9MLH1wMFpojC2r5SeKxolWaNsG1TD3/wZi3v2JC7CbdmcfW84ztG09zLynH5KSnJsFUp6Zov5X3d9hRS+Tgp5WNSQUam5a54SDTdQ8OVKSSTUjoJKrNBJmxYXJSlC/ekXc0jcTdauWx3k6iZSiqVu1Avq1iG+9VrM0qFExJpVDBJm8972b1dF98X15e298Hp4Us/wsaZy/uQX/dy2ICHWraS+yYiaPjR6RYZqeBkFVRIScCSMw46/9VL+XWWxC//za37WqXxUNqlfARR1D4/5mWu11NXKxG1Ulq0J7xyB9PbjHoSymWpRJQX9kZBtMuLS9qi/s8Pb1kfXwENUoQpkW84MuTzreOagZLmhfD83rVsXSh4agblVtkdBr9U24tD2GSaRUMEa+vVvih25jE7lr45qaq9jIPsNaLn93DzYvsm5H+DesVQlz/28g7h9qPFGrx7OXWp/jAUrf+4nXlh7jL5+SjIXjBuOfV3TEtT0b4aHwMyXjVfTgBdaeP6t1rppiUQy1HpJML9NOyqSgV0pNwVU9Gml2UaMnm2Qfwl5NzXkG3D5Q/+1upe/fr0WaYdf7h7/3BwBc2kU7G+JVPRqdmUyTKr91FW4angQsH5U+1z7xC9tms5jKrpwVeSG4IebRAXZAaEw8JU4XXb3YGVl3zgil0iYrzpuakoSnLm6P2ipzY32b10ESATdExZFoeY8kYlI0+go9m9ZKuCt0mRR0p5g8pieqpKagf8s0w7S2ORNG4ApJ/2S7q2LzulWRM2EEWteXS1Dl9NT8a9d0xvs3dS9xF4z0HBrXklsCrFPDGlK9iRi3RcmH3Mj9VR5SbQna9ft6ZSGuM5OickMlA1qeGearpeNZpjymbrUK2PLsCOkc+UbYoffRLppv/LkriAhPXdwu/pNLIlVTiWgYEW0komwiGqvy/SgiWk1EK4koi4jkstn4DBldS0oifHhzD/TTSDHwyehz8OKf5DxaIrPm1Vz0h08E1SqUw0CFj2+tyql4+/pumHS9nLtdakqSrmtekzpVcEXXDLzx59L7yD7DDRxetMEun+R6OmPwRpiZa7EbpVBHe9goMdvyTzRaL7BrezZWDVhzAkM/dCJKBvA6gPMA5AJYRkTThRC/KXabA2C6EEIQUQcAnwGw1x1Eg5ZnJSIbnn300QhZVqNahXJ4/KK2OLdVXfR7fq6DVhkQR+slWqxkF+i2U2CSkwgvXGGfW6hpdHQoOYnke0oqPHdZe7y/cBsA68MKq8efj4omXyrxSqtHOhOO4sYLSKaF3gNAthBiixAiH8AUAKOUOwghjokz1ldGAvOBTrszcZ0BM5XQym+5evz5eDyc8Cfi93tD70w0NLkWqTuo352MmpXw1nVdsfTBwZh+Vx/LGSedwNDzQWN7ZE1U2QdW6Zcdzb0mJlbV1qS9snsjzLqntCeK2WXgqlUoZ+OwEqSmLKy8e5wYB9dyI9VisMU0vIlCJlI0HcAOxedcADGJHojoEgDPAqgLIHYmJrTPGABjAKBRo0ZmbVVFtsVnB1qP76Cz0zB3Y17c569WoRyu79UY1/ZsLJWxzVtoi9vQ8Pi2nlueXbSJo7WrRbRun1W9AjbuOWp43LQ7+2BZzoGSz8pftHHtSnjmkvboaWIiffa9/ZFfqB+JOvve/qhbtTw6P/m99HmtIFM73fLoefnKTlix/aDUvmqeXnatyevGoiUygq5mVczTK4T4CsBXRNQfwJMAhqjsMwnAJADo1q2btwfEFBj9LJOu7xbzoFn9LYkIyX7Tco+w6enh0p4niaBjwxro2LAG9hw5pfq9meE3IUKNF6MGzNn1Qr2HJ0a1dTQXUdxDLg7+Thd3Tje9mEyED2/uoZnfxixuDLnICHouAKU7RgaAXVo7CyHmE1EzIqojhNgXr4FqeOiZBRDygrC1y+oThrTWzhnjBnb/BmqP43s3disZs5Z9XN2Yy7u+V2ZCr3dTn0x0blTT0Ws4IZDRp+zfMv7AOl3vHodnD2SegGUAWhBREyJKBXAVgOnKHYioOYVLQURdAKQC2G+3sW4R+YF8NwriMJNUVn+Kh8bhBYydWpE+GvmUAWf+PrfVWZYfSSLrLVuvNWKiaZ9evSSy043c64+ObIP/SS4VFw+3GkSKA+564xi20IUQhUR0F4DZAJIBvCeEWEdEt4W/nwjgMgDXE1EBgJMArhRe9zEyQYf06ri5TxPcpLG4sRrBKb028axw3jGjekw2y66Na+KHvw9AszR7urxGxN1aSuBv7LX6pDfHE8kcqRWhaxW9lu/NplaIsk6PJrWQVrW84VyGW0ilzxVCzAQwM2rbRMXfzwF4zl7TvENSEuFRg0ChIONE41Br4V0786N4hVKLILhoh53o5S3v1aw2xl/YBpd11Y5E9hpmfpcl4wbrn0s1K6RJgyxS9gZ+E4RbXWSz6VfbNkhMwIMXseq2aPl6IJ+simpMtGgpW+xEhBv7NEHVCv4JiDPzuyQlUVy9UydhQbfAL4qUr15i41PDTC9O8OXtvbFm/Pm6+3w8+hxc2sWa10AQiG5dmX1ZOxnu7wXG9G+KC9rXd9sMU9x3XkvN3PR+xpeC7vYal2YDN5xm8l96Yv79g1QXVzCiQrlkw5ZU98xaePFPnSxa5z8GnR3ydDByEZSd/IukbrDqSud1Hrygte+8vO4e3AL/u/vMsJ+divLq1Z3Rq2ltx1c8U8OXS9D5iUQEF/SSWP9x2p19MHejd1wMvcBQjYReT13cHn8d3ELzRWc28KRK+RSse3woKpZLxo6DJwAY52+pnJqM44ql24LOkgcHo7A4GANSfZrXMRVjYCf+eq36iKZ1quDmPk0wycR6jk7SsWEN3Otgsv1Ir8Xq4sJmefv6brigfXz52iuUS8bicYNj0hinpiQho+aZuYjod/Izl7THDb0ao7+JxUAql09BUhKhUa1KuH/o2XjnBn2Xzy/v6C197iBwVrUKnuv5OoHT7TtuoTtEWfOMieRoualPYtzHzmtzli0JvOpVr4AHhrXCmz9tlj6mbrUKeHyUtZSoRIQ7BzW3dCzDGMGCzjCSjE6QrzOTOAa1SsOlXdLjXsHJiES5LbKg+4Sf/zEIBUXeDGYAvB/JaMSno8+B0RBuy3rqCziXZZomKAjMKcqnJAdqwt+Xgu4l8XBkUWgVvJ5C12uRjGbprTOJ5UbZImP46TUqYuehk4k3QIJv7u6rucA34w48KRonEy6Lb7FchlGjSvkU5EwYgYs6NXDbFE3apcemb3CCyPqnlcv7sv2ZUPgOWeSZS9rjZEHZcSsry7jZI/RQZ9Q1OmZUx9jhrXB5nKkE/n5eS7z4/e82WeVNWNAtcs059izQwTCMPkSE2ySyHBrx18EtWNAZJsLUO3rjLI1Vh7w0r8EwXiOzTiX8tvsIKjm8whoLOiNNF50FDPw+KcowTvL85R1xaecMNE1zNpuoLydFuTHIMIyfqFI+BUNsCIQzglvoTFykJJHmMAxjHwFaL4ZxEF+20BnvsOHJYZj/j0Eln4PYe7qlbxNUq5CCgWfHv94kUxo70jcwZ+AWOhMXKT5Lm2qF1vWrYfX4oa7akIisnYlm3eNDUT4l+PUnkbCgM4wPCOKQi1uBQpUd9jRxExZ0hvEwAWyYu8pr13RGh/QabpvhGL4U9CB2P/3Ooxe2ARHQt4U7if2DSgAb5q4ysoN3UynYgS8FnfEejWtXxjs3dHfbDIYp0/CMBMN4GO6MMmZgQWcYDxNZfLkseBMx8cNDLgzjYW7t3wwn84twY+9Mt01hfAALOsN4mIqpyRh3QWu3zWB8gi/7cTysyDAME4svBZ1hGIaJhQWdYRgmILCgMwzDBAQWdIZhmIDAgs4wDBMQfCnoHD3HMAwTiy8FnWEYhomFA4sYhmHioEXdKmifUd1tMwBICjoRDQPwbwDJAN4RQkyI+v7PAB4IfzwG4HYhxCo7DWUYhvEi3/99gNsmlGA45EJEyQBeBzAcQBsAVxNRm6jdtgIYIIToAOBJAJPsNpRhGIbRR2YMvQeAbCHEFiFEPoApAEYpdxBCLBRCHAx/XAwgw14zGYZhGCNkBD0dwA7F59zwNi1uATBL7QsiGkNEWUSUlZeXJ29l7HksH8swDBNUZARdTT1VF8YiokEICfoDat8LISYJIboJIbqlpaXJW8kwDMMYIjMpmgugoeJzBoBd0TsRUQcA7wAYLoTYb4958TPnvgFI4hY9wzBlABlBXwagBRE1AbATwFUArlHuQESNAEwFcJ0Q4nfbrYyDZmlV3DaBYRgmIRgKuhCikIjuAjAbIbfF94QQ64jotvD3EwE8CqA2gDfC49uFQohuzpnNMAzDRCPlhy6EmAlgZtS2iYq/RwMYba9pDMMwjBk49J9hGCYgsKAzDMMEBBZ0hmGYgMCCzjAMExB8LeipKb42n2EYxlZ8mz734RGt0b8lR5syDMNE8K2gj+7X1G0TGIZhPAWPWTAMwwQEFnSGYZiAwILOMAwTEFjQGYZhAgILOsMwTEBgQWcYhgkILOgMwzABgQWdYRgmIJAQqsuDOn9hojwA2yweXgfAPhvN8QNc5rIBl7lsEE+ZGwshVMPkXRP0eCCirLK2IhKXuWzAZS4bOFVmHnJhGIYJCCzoDMMwAcGvgj7JbQNcgMtcNuAylw0cKbMvx9AZhmGYWPzaQmcYhmGiYEFnGIYJCL4TdCIaRkQbiSibiMa6bU88ENF7RLSXiNYqttUiou+JaFP4/5qK78aFy72RiIYqtnclojXh714hIkp0WWQgooZENJeI1hPROiK6J7w9yGWuQERLiWhVuMyPh7cHtswRiCiZiH4lom/CnwNdZiLKCdu6koiywtsSW2YhhG/+AUgGsBlAUwCpAFYBaOO2XXGUpz+ALgDWKrY9D2Bs+O+xAJ4L/90mXN7yAJqE70Ny+LulAHoBIACzAAx3u2wa5a0PoEv476oAfg+XK8hlJgBVwn+XA7AEQM8gl1lR9r8D+BTAN0Gv22FbcwDUidqW0DL7rYXeA0C2EGKLECIfwBQAo1y2yTJCiPkADkRtHgXgg/DfHwC4WLF9ihDitBBiK4BsAD2IqD6AakKIRSJUGz5UHOMphBC7hRArwn8fBbAeQDqCXWYhhDgW/lgu/E8gwGUGACLKADACwDuKzYEuswYJLbPfBD0dwA7F59zwtiBxlhBiNxASQAB1w9u1yp4e/jt6u6chokwAnRFqsQa6zOGhh5UA9gL4XggR+DIDeBnAPwAUK7YFvcwCwHdEtJyIxoS3JbTMflskWm0sqaz4XWqV3Xf3hIiqAPgSwL1CiCM6Q4SBKLMQoghAJyKqAeArImqns7vvy0xEIwHsFUIsJ6KBMoeobPNVmcP0EULsIqK6AL4nog06+zpSZr+10HMBNFR8zgCwyyVbnGJPuNuF8P97w9u1yp4b/jt6uychonIIifknQoip4c2BLnMEIcQhAD8BGIZgl7kPgIuIKAehYdFziehjBLvMEELsCv+/F8BXCA0RJ7TMfhP0ZQBaEFETIkoFcBWA6S7bZDfTAdwQ/vsGANMU268iovJE1ARACwBLw924o0TUMzwbfr3iGE8Rtu9dAOuFEC8qvgpymdPCLXMQUUUAQwBsQIDLLIQYJ4TIEEJkIvSM/iiEuBYBLjMRVSaiqpG/AZwPYC0SXWa3Z4YtzCRfgJB3xGYAD7ltT5xlmQxgN4AChN7MtwCoDWAOgE3h/2sp9n8oXO6NUMx8A+gWrjybAbyGcASw1/4B6ItQ93E1gJXhfxcEvMwdAPwaLvNaAI+Gtwe2zFHlH4gzXi6BLTNCnnerwv/WRbQp0WXm0H+GYZiA4LchF4ZhGEYDFnSGYZiAwILOMAwTEFjQGYZhAgILOsMwTEBgQWcYhgkILOgMwzAB4f8BIR7ERoBELAQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.hstack(q6s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding a clustering condition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now modify the above function to also find clusters which satisfy particular $q_6$ value. But first, for a single file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys = pc.System()\n",
    "sys.read_inputfile(files[0])\n",
    "sys.find_neighbors(method=\"cutoff\", cutoff=0)\n",
    "sys.calculate_q(6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now a clustering algorithm can be applied on top using the `cluster_atoms` method. `cluster_atoms` uses a `condition as argument` which should give a True/False value for each atom. Lets define a condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def condition(atom):\n",
    "    return atom.get_q(6) > 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above function returns `True` for any atom which has a $q_6$ value greater than 0.5 and `False` otherwise. Now we can call the `cluster_atoms` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sys.cluster_atoms(condition)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The method returns 16, which here is the size of the largest cluster of atoms which have $q_6$ value of 0.5 or higher. If information about all clusters are required, that can also be accessed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "atoms = sys.atoms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`atom.cluster` gives the number of the cluster that each atom belongs to. If the value is -1, the atom does not belong to any cluster, that is, the clustering condition was not met."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusters = [atom.cluster for atom in atoms if atom.cluster != -1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can see how many unique clusters are there, and what their sizes are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_clusters, counts = np.unique(clusters, return_counts=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`counts` contain all the necessary information. `len(counts)` will give the number of unique clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Cluster ID')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWeElEQVR4nO3dfZQldX3n8fdn8AFBFFgagkg7wFE8iorYZjVmFTEkJCDEJwy7cvAhGTdHDWp8wKjBrHFFo0bUHM3sESEnBB9RUVRAIpDNSZBhAgICojghKDooGhEVxPnuH7dmaZrunuqHqkt3vV/n9OlbdevW71vUzHd+/O6vvr9UFZKk4Vgz7gAkSf0y8UvSwJj4JWlgTPySNDAmfkkamPuMO4A2dtttt1q7du24w5CkFeXSSy/9QVVNzNy/IhL/2rVr2bBhw7jDkKQVJcm/z7bfoR5JGhgTvyQNjIlfkgbGxC9JA2Pil6SBMfFL0sB0lviTnJJkc5IrZ+x/RZJrk1yV5J1dtS9Jml2XPf5TgcOm70jydOAo4LFV9WjgXR22L0maRWeJv6ouAm6ZsfuPgZOq6vbmmM1dtS9Jml3fT+4+AvhvSd4G/AJ4TVVdMtuBSdYB6wAmJyf7i1AamLUnnL3NYzaddHgPkagvfX+5ex9gF+BJwGuBjyfJbAdW1fqqmqqqqYmJe5SakCQtUt+J/0bgzBr5KrAF2K3nGCRp0PpO/J8BDgFI8gjgfsAPeo5BkgatszH+JGcABwO7JbkROBE4BTilmeJ5B3Bcudq7JPWqs8RfVcfM8dYLumpTkrRtPrkrSQNj4pekgTHxS9LAmPglaWBM/JI0MCZ+SRoYE78kDYyJX5IGxsQvSQNj4pekgTHxS9LAmPglaWBM/JI0MCZ+SRoYE78kDYyJX5IGprPEn+SUJJub1bZmvveaJJXE9XYlqWdd9vhPBQ6buTPJ3sChwA0dti1JmkNnib+qLgJumeWtvwZeB7jWriSNQa9j/EmOBL5TVZf32a4k6S6dLbY+U5IdgDcCv93y+HXAOoDJyckOI5OkYemzx78fsA9weZJNwEOBjUl+bbaDq2p9VU1V1dTExESPYUrS6tZbj7+qrgB237rdJP+pqvpBXzFIkrqdznkG8C/A/kluTPKSrtqSJLXXWY+/qo7Zxvtru2pbkjS3eXv8SdYkObqvYCRJ3Zs38VfVFuDlPcUiSepBmzH+85oSC3sn2XXrT+eRSZI60WaM/8XN75dN21fAvssfjiSpa9tM/FW1Tx+BSJL6sc2hniQ7JHlTkvXN9sOTHNF9aJKkLrQZ4/8IcAfwG832jcBfdhaRJKlTbRL/flX1TuCXAFX1cyCdRiVJ6kybxH9HkgfQlFFOsh9we6dRSZI602ZWz1uALwF7JzkdeArwoi6DkiR1p82snnOTXAo8idEQz/EWVpOklavNrJ7zq+qHVXV2VX2+qn6Q5Pw+gpMkLb85e/xJtgd2AHZLsgt3faH7IOAhPcQmSerAfEM9LwVeySjJX8pdif8nwN90G5YkqStzJv6qOhk4Ockrqur9PcYkSepQm+mc30uyE0DzBO+ZSQ7qOC5JUkfaJP43V9WtSX4T+B3gNOCD3YYlSepKm8T/q+b34cAHq+qzwP229aEkpyTZnOTKafv+Ksk1Sb6W5NNJdl5U1JKkRWuT+L+T5G+Bo4EvJLl/y8+dChw2Y995wAFV9VjgG8AbFhCrJGkZtEngRwPnAIdV1Y+BXYHXbutDVXURcMuMfedW1Z3N5r8CD11QtJKkJWuT+HcDNgC3J5kE7gtcswxtvxj44lxvJlmXZEOSDTfffPMyNCdJgna1es5mVKAtwPbAPsC1wKMX22iSNwJ3AqfPdUxVrQfWA0xNTdVi25Ik3V2bWj2Pmb7dTOV86WIbTHIccATwjKoyoUtSz9r0+O+mqjYmeeJiGktyGPB64GlV9bPFnEOStDTbTPxJXj1tcw1wELDNQfckZwAHM6r1cyNwIqNZPPcHzksC8K9V9T8XHrYkabHa9Ph3mvb6TkZj/p/a1oeq6phZdn+4ZVySpI60GeP/iz4CkST1Y76yzJ+jWW5xNlV1ZCcRSZI6NV+P/129RSFJ6s18ZZkvBEiyI/DzqtrSbG/H6AtaSdIK1ObJ3fMZrcS11QOAL3cTjiSpa20S//ZV9dOtG83rHeY5XpJ0L9Ym8d82feGVJE8Aft5dSJKkLrWZx/9K4BNJvtts7wk8v7OIJEmdajOP/5IkjwT2Z1So7Zqq+mXnkUmSOtGqVk+T6K/c5oGSpHu9NmP8kqRVxMQvSQPTaqgnyV7Aw6Yf3yytKElaYdqUZX4Ho1k8Xwd+1ewuwMQvSStQmx7/7wP7V9XtHcciSepBmzH+6xktsC5JWgXa9Ph/BlyW5Hzg//f6q+pPOotKktSZNon/rOZnQZKcwmhR9c1VdUCzb1fgY8BaYBNwdFX9aKHnliQtXpsnd09b5LlPBT4A/N20fScA51fVSUlOaLZfv8jzS5IWYb4VuD5eVUcnuYJZVuKqqsfOd+KquijJ2hm7j2K0ADvAacAFmPglqVfz9fiPb34fsYzt7VFVNwFU1U1Jdp/rwCTrgHUAk5OTyxiCJA3bfCtwbU3Q/95fOHdrfz2wHmBqamrOtX8lSQvTd8mG7yfZE6D5vbnn9iVp8PpO/GcBxzWvjwM+23P7kjR4C0r8SXZJMu+XutOOPQP4F2D/JDcmeQlwEnBokuuAQ5ttSVKP2tTquQA4sjn2MuDmJBdW1avn+1xVHTPHW89YYIySpGXUpsf/4Kr6CfBs4CNV9QTgt7oNS5LUlTaJ/z7NF7FHA5/vOB5JUsfaJP7/BZwDfLNZf3df4Lpuw5IkdaVNyYZPAJ+Ytn098Jwug5IkdafNl7v7AK9gVFht+gpcR3YXliSpK22qc34G+DDwOWBLp9FIkjrXJvH/oqre13kkkqRetEn8Jyc5ETiXuy/EsrGzqCRJnWmT+B8DHAscwl1DPdVsS5JWmDaJ/1nAvlV1R9fBSJK612Ye/+XAzh3HIUnqSZse/x7ANUku4e5j/E7nlKQVqE3iP7HzKCRJvWnz5O6FSfYAntjs+mpVuYCKJK1Q2xzjT3I08FXgeYwKtV2c5LldByZJ6kaboZ43Ak/c2stPMgF8Gfhkl4FJkrrRZlbPmhlDOz9s+TlJ0r1Qmx7/l5KcA5zRbD8f+OJSGk3yKuAPGT0IdgXwoqr6xVLOKUlqZ5s996p6LfC3wGOBxwHrq+p1i20wyV7AnwBTVXUAsB3wB4s9nyRpYdqUZX5HVb0eOHOWfUtp9wFJfgnsAHx3CeeSJC1Am6GeQ4GZSf53Z9nXSlV9J8m7gBuAnwPnVtW5M49Lsg5YBzA5ObmYpiR1YO0JZ2/zmE0nHd5DJFqsOYd6kvxxkiuA/ZN8bdrPt4GvLbbBJLsARwH7AA8BdkzygpnHVdX6qpqqqqmJiYnFNidJmmG+Hv8/MPoS9+3ACdP231pVtyyhzd8Cvl1VNwMkORP4DeDvl3BOSVJLcyb+qvpP4D+BYwCS7A5sDzwwyQOr6oZFtnkD8KQkOzAa6nkGsGGR55IkLVCbJ3efmeQ64NvAhcAmljCds6ouZvTw10ZGUznXAOsXez5J0sK0eRDrL4EnAd+oqn0Y9dD/eSmNVtWJVfXIqjqgqo6tqtu3/SlJ0nJok/h/WVU/BNYkWVNVXwEO7DYsSVJX2kzn/HGSBwIXAacn2Qzc2W1YkqSutOnxHwX8DHgV8CXgW8AzuwxKktSdNvX4b2tebgFO6zYcSVLXrLIpSQNj4pekgZmvZMP5ze939BeOJKlr843x75nkacCRST4KZPqbVbWx08gkSZ2YL/H/OaMaPQ8F3jPjvQIO6SooSVJ35qvV80ngk0neXFVv7TEmSVKH2kznfGuSI4GnNrsuqKrPdxuWJKkrbYq0vR04Hvh683N8s0+StAK1KdlwOHBgVW0BSHIa8G/AG7oMTJLUjbbz+Hee9vrBHcQhSepJmx7/24F/S/IVRlM6n4q9fUlasdp8uXtGkguAJzJK/K+vqu91HZgkqRttevxU1U3AWR3HIknqwVhq9STZOcknk1yT5OokTx5HHJI0RK16/B04GfhSVT03yf2AHcYUhyQNzrw9/iRrkly5nA0meRCjL4g/DFBVd1TVj5ezDUnS3Obt8VfVliSXJ5msqhuWqc19gZuBjyR5HHApcPy0BV8ASLIOWAcwOTm56MbWnnB2q+M2nXT4otuQpJWkzRj/nsBVSc5PctbWnyW0eR/gIOCDVfV44DZGxeDupqrWV9VUVU1NTEwsoTlJ0nRtxvj/YpnbvBG4saoubrY/ySyJX5LUjTbz+C9M8jDg4VX15SQ7ANsttsGq+l6S/0iyf1VdCzyDUQ0gSVIPtpn4k/wRo7H2XYH9gL2ADzFK2Iv1CuD0ZkbP9cCLlnAuSdICtBnqeRnw68DFAFV1XZLdl9JoVV0GTC3lHJKkxWnz5e7tVXXH1o0k92G0ApckaQVqk/gvTPJnwAOSHAp8Avhct2FJkrrSJvGfwGje/RXAS4EvAG/qMihJUnfazOrZ0iy+cjGjIZ5rq8qhHklaodrM6jmc0SyebzEqy7xPkpdW1Re7Dk6StPzazOp5N/D0qvomQJL9gLMBE78krUBtxvg3b036jeuBzR3FI0nq2Jw9/iTPbl5eleQLwMcZjfE/D7ikh9gkSR2Yb6jnmdNefx94WvP6ZmCXziKSJHVqzsRfVZZRkKRVqM2snn0Y1dZZO/34qjqyu7AkSV1pM6vnM4xWy/ocsKXTaCRJnWuT+H9RVe/rPBJJUi/aJP6Tk5wInAvcvnVnVW3sLCpJUmfaJP7HAMcCh3DXUE8125KkFaZN4n8WsO/00sySpJWrzZO7lwM7dxyHJKknbXr8ewDXJLmEu4/xL2k6Z5LtgA3Ad6rqiKWcS5LUXpvEf2JHbR8PXA08qKPzS5Jm0aYe/4XL3WiShwKHA28DXr3c55ckza3Nk7u3ctcau/cD7gvcVlVL6am/F3gdsNM87a4D1gFMTk4uoSlJK8XaE85uddymkw7vOJLVbZtf7lbVTlX1oOZne+A5wAcW22CSIxiVer50G+2ur6qpqpqamJhYbHOSpBnazOq5m6r6DEubw/8U4Mgkm4CPAock+fslnE+StABthnqePW1zDTDFXUM/C1ZVbwDe0Jz7YOA1VfWCxZ5PkrQwbWb1TK/LfyewCTiqk2gkSZ1rM6uns7r8VXUBcEFX55ck3dN8Sy/++Tyfq6p6awfxSJI6Nl+P/7ZZ9u0IvAT4L4CJX5JWoPmWXnz31tdJdmL0pO2LGM3Eefdcn5Mk3bvNO8afZFdGT9b+D+A04KCq+lEfgUmSujHfGP9fAc8G1gOPqaqf9haVJKkz8z3A9afAQ4A3Ad9N8pPm59YkP+knPEnScptvjH/BT/VKku792jzANShtikRZIEoz+edGK4m9ekkaGBO/JA2MiV+SBsbEL0kDY+KXpIEx8UvSwJj4JWlgTPySNDAmfkkamN4Tf5K9k3wlydVJrkpyfN8xSNKQjaNkw53An1bVxqbO/6VJzquqr48hFkkanN57/FV1U1VtbF7fClwN7NV3HJI0VGMt0pZkLfB44OJZ3lsHrAOYnJzsN7CB6LqwWJfnb3PupZy/Sys5di2fcf45GNuXu0keCHwKeGVV3aO+f1Wtr6qpqpqamJjoP0BJWqXGkviT3JdR0j+9qs4cRwySNFTjmNUT4MPA1VX1nr7bl6ShG0eP/ynAscAhSS5rfn5vDHFI0iD1/uVuVf1fIH23K0ka8cldSRoYE78kDYyJX5IGxsQvSQNj4pekgTHxS9LAmPglaWDGWqRtNVhoIbKFHL/QIk5dF31ayUXXurxPXeviz8G9NfYuzr+U+9rl39dxsscvSQNj4pekgTHxS9LAmPglaWBM/JI0MCZ+SRoYE78kDYyJX5IGxsQvSQMzrsXWD0tybZJvJjlhHDFI0lCNY7H17YC/AX4XeBRwTJJH9R2HJA3VOHr8vw58s6qur6o7gI8CR40hDkkapFRVvw0mzwUOq6o/bLaPBf5rVb18xnHrgHXN5v7AtcsYxm7AD5bxfPdmXuvq5LWuTst9rQ+rqomZO8dRnTOz7LvHvz5VtR5Y30kAyYaqmuri3Pc2Xuvq5LWuTn1d6ziGem4E9p62/VDgu2OIQ5IGaRyJ/xLg4Un2SXI/4A+As8YQhyQNUu9DPVV1Z5KXA+cA2wGnVNVVPYfRyRDSvZTXujp5ratTL9fa+5e7kqTx8sldSRoYE78kDcygEv+QSkUk2ZTkiiSXJdkw7niWW5JTkmxOcuW0fbsmOS/Jdc3vXcYZ43KZ41rfkuQ7zf29LMnvjTPG5ZBk7yRfSXJ1kquSHN/sX3X3dZ5r7eW+DmaMvykV8Q3gUEZTSi8Bjqmqr481sI4k2QRMVdWqfPAlyVOBnwJ/V1UHNPveCdxSVSc1/7DvUlWvH2ecy2GOa30L8NOqetc4Y1tOSfYE9qyqjUl2Ai4Ffh94Iavsvs5zrUfTw30dUo/fUhGrSFVdBNwyY/dRwGnN69MY/UVa8ea41lWnqm6qqo3N61uBq4G9WIX3dZ5r7cWQEv9ewH9M276RHv9Dj0EB5ya5tCl/MQR7VNVNMPqLBew+5ni69vIkX2uGglb88Md0SdYCjwcuZpXf1xnXCj3c1yEl/lalIlaRp1TVQYyqoL6sGS7Q6vFBYD/gQOAm4N1jjWYZJXkg8CnglVX1k3HH06VZrrWX+zqkxD+oUhFV9d3m92bg04yGula77zdjp1vHUDePOZ7OVNX3q+pXVbUF+D+skvub5L6MEuHpVXVms3tV3tfZrrWv+zqkxD+YUhFJdmy+MCLJjsBvA1fO/6lV4SzguOb1ccBnxxhLp7YmwsazWAX3N0mADwNXV9V7pr216u7rXNfa130dzKwegGZq1Hu5q1TE28YbUTeS7Muolw+jshz/sNquNckZwMGMyth+HzgR+AzwcWASuAF4XlWt+C9F57jWgxkNBxSwCXjp1nHwlSrJbwL/BFwBbGl2/xmjse9VdV/nudZj6OG+DirxS5KGNdQjScLEL0mDY+KXpIEx8UvSwJj4JWlgTPxa9ZL8WpKPJvlWkq8n+UKSRyRZO73i5QLP+cIkD1liXC9M8oHm9fSqjNclOTPJo5ZyfmkuJn6tas2DMp8GLqiq/arqUYzmS++xxFO/EFhQ4k+yraVO/7qqDqyqhwMfA/4xycQi45PmZOLXavd04JdV9aGtO6rqsqr6p+kHTe99N9ufT3Jwku2SnJrkymZ9g1cleS4wBZze9NAfkOQJSS5siuKdM63EwAVJ/neSC4Hj2wZdVR8DzgX++9IuX7qn3hdbl3p2AKNa54t1ILDXtDr4O1fVj5O8HHhNVW1oaq68Hziqqm5O8nzgbcCLm3PsXFVPW0TbG4FHLiF2aVYmfml+1wP7Jnk/cDajXvhM+zP6B+a80cgS2zGqrLjVxxbZ9mwVZaUlM/FrtbsKeG6L4+7k7kOf2wNU1Y+SPA74HeBljFZIevGMzwa4qqqePMe5b1tQxHd5PLDqls3U+DnGr9XuH4H7J/mjrTuSPDHJzKGXTcCBSdYk2ZumHG6S3YA1VfUp4M3AQc3xtwI7Na+vBSaSPLn5zH2TPHopQSd5DqOqqmcs5TzSbOzxa1WrqkryLOC9zXqtv2CU5F8549B/Br7NqFrilYzG12G0SttHkmztJL2h+X0q8KEkPweezOj/Kt6X5MGM/l69l9H/bSzEq5K8ANixieGQqrp5geeQtsnqnJI0MA71SNLAmPglaWBM/JI0MCZ+SRoYE78kDYyJX5IGxsQvSQPz/wDLaJymrQJulgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(range(len(counts)), counts)\n",
    "plt.ylabel(\"Number of atoms in cluster\")\n",
    "plt.xlabel(\"Cluster ID\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can finally put all of these together into a single function and run it over our individual time slices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_q6_cluster(file, cutoff_q6 = 0.5, format=\"lammps-dump\"):\n",
    "    sys = pc.System()\n",
    "    sys.read_inputfile(file, format=format)\n",
    "    sys.find_neighbors(method=\"cutoff\", cutoff=0)\n",
    "    sys.calculate_q(6)\n",
    "    def _condition(atom):\n",
    "        return atom.get_q(6) > cutoff_q6\n",
    "    sys.cluster_atoms(condition)\n",
    "    atoms = sys.atoms\n",
    "    clusters = [atom.cluster for atom in atoms if atom.cluster != -1]\n",
    "    unique_clusters, counts = np.unique(clusters, return_counts=True)\n",
    "    return counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "q6clusters = [calculate_q6_cluster(file) for file in files]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the number of clusters for each slice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'number of unique clusters')"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA23klEQVR4nO3deXyU9bX48c/JBkmAhISwhSUJIrJvkYRNEResdatLW6toBfSqvb3aWlptf7f7rba23t7aFUGl1bpU0bZWDYoba9i3sLgkbIEsEAIhZM/5/TETGpCQCcwzzyzn/Xo9L2aezMxzHJMzz3yf7zlfUVWMMcZEjii3AzDGGBNYlviNMSbCWOI3xpgIY4nfGGMijCV+Y4yJMDFuB+CLHj16aEZGhtthGGNMSFm3bt1BVU07dX9IJP6MjAzWrl3rdhjGGBNSRGT36fbbUI8xxkQYS/zGGBNhLPEbY0yEscRvjDERxhK/McZEGMdm9YhIZ+BDoJP3OC+r6g9EJAV4EcgAdgFfVNXDTsVhjDEd8dqGYh7L28n+yhr6Jsczd8YQrh+b7nZYfuXkGX8dMF1VRwNjgCtFJBd4CFiiqoOBJd77xhjjutc2FPPwoi0UV9agQHFlDQ8v2sJrG4rdDs2vHEv86nHMezfWuylwHbDQu38hcL1TMRhjTEc8lreTmoamk/bVNDTxWN5OlyJyhqNj/CISLSIbgTLgbVXNB3qp6gEA778923ju3SKyVkTWlpeXOxmmMcYAsL+ypkP7Q5WjiV9Vm1R1DNAPmCAiIzrw3Hmqmq2q2Wlpn6k4NsYYv+ubHH/a/VFRwtPLi6iuawxwRM4IyKweVa0E3geuBEpFpA+A99+yQMRgjDHtmTtjCCIn74uLFvp3j+dH/9zGxEeW8OibOyg5UutOgH7iWOIXkTQRSfbejgcuA3YA/wDu8D7sDuDvTsVgjDEdcWFmCqrQrXMMAqQnx/OLm0bz/txLWHTfJKYM7sG8Dz9l6i/e5ZsvbmTb/qNuh3xWnGzS1gdYKCLReD5gXlLV10VkJfCSiMwG9gA3OxiDMcb4bHFBCQCvfW0yWWldTvrZuAHd+f2t49lz6DhPLS/ipbV7WbShmCnn9WDO1EwuPj8NOfXrQpCSUFhsPTs7W607pzHGaV+et5KK6noWf+Pidh975HgDz63ezcIVuyg9Wsf5vbowZ0oW143tS6eY6ABE2z4RWaeq2afut8pdY4wBKqrrWV1UwYzhvX16fFJCLPdNO4+l357Or24eTZQI335lM5MffY8nlnzM4ep6hyM+eyHRj98YY5z2zrZSmhWfE3+LuJgobhzfjxvGpbP8k0M8ubSQX739Eb97/xNuHt+f2VMyyeiR6FDUZ8cSvzHGAHkFJaQnxzO8b7ezer6IMGVwD6YM7sFHpVXMX1rIi2v28mz+bi4f2ou7Lsoie2D3oLgOYInfGBPxjtU1svSTg9yWM9Avifn8Xl35xU2j+daMIfx5xW6ezd/N4m2ljO6fzF1TM7lyeG9iot0babcxfmNMxHt/Zxn1jc1cOaJjwzzt6dm1M9+aMYQVD03nJ9cN58jxev7zrxuY9sv3WbCsiGMuFYRZ4jfGRLy8glJSE+MYP7C7I6+fEBfDzIkZLHlwGn+aOZ4+SZ35yeuegrBH3tzOgSOBbQlhQz3GmIhW19jEezvKuHpUH6KjnB1/j44SZgzvzYzhvdm4t5Inlxby5IeFLFhaxDWj+zJnaibD+yYBzraHtsRvjIloKz45xLG6Rmb4eZinPWP6J/O7r4xjb8Vxnl6+ixfX7OHVDcVMGpTKiL7d+Muq3dQ0NAP/bg8N+CX521CPMSai5RWU0KVTDJMGpbpy/P4pCXz/mmGsePhSHvrcBRSWVzNvadGJpN/Cn+2hLfEbYyJWU7Py9rZSLrmgp+vVtknxsdxz8SA+/PYlbT7GX+2hLfEbYyLW2l0VHKqu58oOFm05KS4mivQ22kO31Ta6oyzxG2MiVl5BKXExUUwbElxrfsydMYT42JO/gcTHRjN3xhC/vL5d3DXGRCRVJa+ghKnn9SCxU3ClwpYLuDarxxhj/Khg/1GKK2u4/7LBbodyWtePTfdboj+VDfUYYyJSXkEJUQKXDe3ldigBZ4nfGBOR8gpKmJCZQkpinNuhBJwlfmNMxCksP8ZHpceCajZPIFniN8ZEnLyCUgCusMRvjDGRIa+ghFH9kvw2Lz7UWOI3xkSUkiO1bNxb2eGVtsKJJX5jTERZvK0E6PgSi+HEscQvIv1F5D0R2S4iBSJyv3f/GBFZJSIbRWStiExwKgZjjDlVXkEJg9ISOa9nF7dDcY2TBVyNwIOqul5EugLrRORt4BfAj1T1TRG5ynt/moNxGGPOwMm+78HmcHU9qworuOfiLLdDcZVjiV9VDwAHvLerRGQ7kA4o0LKacRKw36kYjDFn9tqGYh5etIWahibA/33fg82SHWU0NWtED/NAgMb4RSQDGAvkAw8Aj4nIXuCXwMOBiMEY81mP5e08kfRb+LPve7DJKyihb1JnRqYnuR2KqxxP/CLSBXgFeEBVjwL3At9Q1f7AN4AFbTzvbu81gLXl5eVOh2lMRGqrv7u/+r4Hk+P1jXz4UTlXDO+NiLNLLAY7RxO/iMTiSfrPqeoi7+47gJbbfwNOe3FXVeeparaqZqelBVfLVGPCRVvz2MNxfvsHO8upa2yO+GEecHZWj+A5m9+uqo+3+tF+4GLv7enAx07FYIw5s69NH/SZfZ1iovzW9z2YvFVQQveEWC7M6O52KK5zclbPZGAmsEVENnr3fRe4C/g/EYkBaoG7HYzBGHMGVTWNAKR17cTBqjoUGNs/Oewu7NY3NvPujjI+N6I3MdFWvtShxC8i3YH+qrq5vceq6jKgrYG08R05rjHG/xqamnlmxS4mZqXy/N25AHzv1S38bd0+Dh6ro0eXTi5H6D8rCw9RVdtowzxe7X70icj7ItJNRFKATcDTIvJ4e88zxgS3N7Yc4MCRWuZMzTyxb9aUTOobm3l21W4XI/O/t7aWkBgXzeTzergdSlDw5TtPknc2zg3A06o6HrjM2bCMMU5SVRYsKyKrRyKXDOl5Yv+gtC5cekFP/rJyN7WnTPMMVU3NytvbSpl2QU86n7KObaTyJfHHiEgf4IvA6w7HY4wJgLW7D7N53xFmTckkKurkEdnZUzM5VF3PPzaGR23lhj2HOXiszoZ5WvEl8f8IyAM+UdU1IpKFzcQxJqTNX1pIckIsN47r95mfTcxKZVifbsxfVoiquhCdf721tYS46CguGWLTwlucMfGLSDSei7mjVPU+AFUtVNUbAxKdMcbvdh+qZvG2Um7NGUB83GeHPkSE2VMy+aj0GEs/PuhChP6jquRtK2Hyeal07RzrdjhB44yJX1WbgGsDFIsxJgCeXr6LmCjh9okZbT7mmtF96dm1E/OXFQUuMAdsP1DF3ooaG+Y5hS9DPStE5LciMlVExrVsjkdmjPG7I8cbeGntXq4Z3Zde3Tq3+bi4mCjumJTBhx+Vs7OkKoAR+tdbBSVECVw2rJfboQQVXxL/JGA48GPgV97tl04GZYxxxvNr9nC8vonZUzLbfexXJgygc2wUT4XwWf/ighKyM1LCqibBH9pN/Kp6yWm26YEIzhjjPw1NzTyzfBeTBqUyvG/73Sm7J8Zx0/h+vLqxmPKqugBE6F+7Dlazo6TKhnlOw5cCrl4iskBE3vTeHyYis50PzRjjT29sOUDJ0ZMLttoza3LoFnTlFXiWWLzChnk+w5ehnmfwTOfs673/EZ6e+saYEHGiYCstkWnn92z/CV5ZaV24bGhPnl0VegVdeQUljEjvRv+UBLdDCTq+JP4eqvoS0Aygqo1AaP0GGBPh1uzyFmxN/mzBVntmTfEUdP19Y7FD0flf2dFa1u+pZMYwG+Y5HV8Sf7WIpOJZMhERyQWOOBqVMcavzlSw1Z4TBV1Li0KmoCtvWykAM0ZY4j8dXxL/N4F/AINEZDnwZ+C/HI3KGOM3uw5W8/b2Um7LGXjagq32iAhzpmbycdkxPgyRgq7FBSVk9UhkcM8ubocSlHxJ/AV4Fk6ZBPwHnqmdO5wMyhjjP08vL/IWbA0869e4epS3oGtpoR8jc8aR4w2s/PSQLbF4Br4k/pWq2qiqBaq6VVUbgJVOB2aMOXeegq19XDs6nZ5nKNhqT0tB19KPDwZ9QdeSHaU0Niszhttsnra0mfhFpLeIjAfiRWRsq6rdaYBdJjcmBPx19R5qGnwr2GrPrTkDiI+NZsGy4D7rzysooXe3zozul+x2KEHrTCtwzQC+CvTDU63b8p2pCs8SisaYINbQ1MzCFbuYfF4qw/p2O+fXS07wFHS9uGYvc2dcQFrX4KuGralv4oOPyvlidv8Oz16KJG2e8avqQlW9BPiqqk5vVbV7raouCmCMxpizcKJga0qW317zzskZNDQHb0HXBx+VU9vQbNW67fBljL+fd+lFEZH5IrJeRK5wPDJjzFlTVZ5cWkhWWiIXn++/PvRZ3hW6grWga3FBCckJsUzITHE7lKDmS+Kf5V168QqgJ3An8KijURljzsnqogq2Fh9l9mlW2DpXs6dkcai6ntc2BFdBV0NTM+9sL+XSC3oRG+1Laotcvrw7Lb81V+FZc3dTq33GmCA0f1kR3RNiuWFsxwu22pOblcLwvt2Yvyy4CrpWFR7iaG2jzebxgS+Jf52ILMaT+PNEpCve9g1nIiL9ReQ9EdkuIgUicn+rn31dRHZ69//i7MNv22sbipn86LtkPvQvJj/6btCdnRjjlKKD1byzvZTbcs+uYKs9LQVdn5Qd44OPyv3++mcrr6CE+NhoLvLj0Fa4OtOsnhazgTFAoaoe97ZvuNOH5zUCD6rqeu+HxToReRvoBVwHjFLVOhHxvWOUj17bUMzDi7ZQ4x2DLK6s4eFFWwC4fmy6vw9nTFB5enkRsVFRzDyHgq32fH5kXx59cwcLlhUxbYjf/4Q7rLlZWVxQyrQhaXSO9f+HXbjx5Yx/CtAFGCUiF+Gp3E1u70mqekBV13tvVwHbgXTgXuBRVa3z/qzs7EJv22N5O08k/RY1DU08lrfT34cyJqhUHq/nb2v3ce2YvvTsevYFW+1pXdC1o+SoY8fx1Ya9lZRV1dlsHh/5kvjnttr+G/gn8MOOHEREMoCxQD5wPjBVRPJF5AMRubCN59wtImtFZG15ece+Tu6vrOnQfmPChT8LttrzlQnegq6l7q/QtbighNho4ZIL3P/2EQp8WYHrmlbb5cAIoNTXA4hIF+AV4AHv7KAYoDuQi+fD5CU5TUMNVZ2nqtmqmp2W1rExu77J8R3ab0w4qG/0FGxNOa8HQ/uce8FWe5IT4rg5ux9/37jf1RW6VJW8ghImDupBUnysa3GEkrOZ87QPT/Jvl4jE4kn6z7Uq+toHLFKP1XguFPc4izjaNHfGEOJPGeeLj41m7owh/jyMMUHljS0HKD1aF5Cz/RZ3Ts6kobmZv7hY0LWztIpdh47bbJ4OaPfirog8gbcXP54PijHAJh+eJ8ACYLuqPt7qR68B04H3ReR8IA7wa6/Xlgu4j+XtpLiyhtho4ZEbRtqFXRO2VJX5ywoZ5OeCrfZk9kjk0gt68eyq3dw3bZArF1bztpYiApfbEos+8+WMfy2wzrutBL6jqrf58LzJwExguohs9G5XAU8BWSKyFXgBuEMdmAx8/dh0lj80na9PP49mhcvsl8KEsfwTBVtZAe9RM2dqJhXV9bzq0pTpvIISxg/o7ujF7HDT7hm/qi48mxdW1WW0XejlyweHX+RkpvLEu5+wdldFUEw7M8YJ85d6C7bGBf5bbU5mCiPSu7FgWRFfCnBztL0Vx9l24Cjfu2powI4ZDs7UlnmLiGw+zbZFRDYHMshzMW5gMrHRQn5RhduhGOOIwvJjLNlRyszcga4MtYgIc6ZkeQq6Pg5sQVdeQQmATePsoDOd8V8dsCgclBAXw6h+yawqPOR2KMY44unlu4iNiuI2Bwu22nPVyD6egq6lRVwSwG/WeQUlDO3TjQGptkRIR5ypLfNuVd3tfUxpq/tlhFivnpzMFLbsO8Lx+ka3QzHGryqP1/Pyun1c53DBVntaCrqWfXKQ7QcCU9BVXlXH2t2HbTbPWfDl4u7fOLk3T5N3X8jIyUqlsVlZt/uw26EY41cnCramBm4KZ1taCrqeWhaYgq63t5WiCleOsGGejvIl8ceoan3LHe/tOOdC8r/sgd2JjhIb7jFhpXXB1gW9nS/Yak9SQuyJgq6yqlrHj/dWQQkDUxMY0qur48cKN74k/nIRubbljohch5/n3TstsVMMI9OTyC+0C7wmfPxry35PwVYQnO23aCnoenalswVdR2sbWPnpQWYM781pCv9NO3xJ/PcA3xWRPSKyB/gOcLezYflfTlYKm/ZVUlMffKsGmcAKh5bdqsr8pUWc17MLFw8OnjbEmT0SuWxoL/7i8Apd7+0oo6FJbTbPWfKlV8+nqpoLDAOGq+okVf3U+dD8KzcrlYYmZcMeG+ePZC0tu4sra1D+3bI71JL/qsIKCvY7s8LWuZozJZPDxxtYtN659/StrSX07NqJsf2THTtGOPO5V4+qHvO2Vw5J2QO7EyXYOH+EC5eW3QuWFZKSGMcXgrANyYTMFEamJ7FgWSHNzf5foau2oYn3d5Zz+bBeQfehFyoiZmHKrp1jGZGexCor5IpYqhoWLbsLy4/xzvYybnOpYKs9LSt0fVpe7cgKXUs/PkhNQ5PN5jkHEZP4wTOff+PeSkfHHk3waWxq5vXN+7n+9yto6/yzV7fQ6fPy1PIi4qKjmJnrXsFWe64a2Yfe3Tozf1mh31/7ra0ldOscQ25Wqt9fO1K0m/hFJEFE/ltEnvTeHywiIVnVm5uVSn1jMxv3VrodigmAY3WNnqUBf/k+//nXDVQer+fGcel0jv3sr31tQyOflAX/SGbrgq20rp3cDqdNsdGegq7lnxzya0FXY1MzS3aUcunQXsRGR9R5q1/58s49DdQBE7339wE/dSwiB2VnpCA2zh/2Dhyp4ZE3tzPxkSX85PVt9O7WmT/eNp53H5zGr744hkdvGEV6cjwCpCfH8+AV5xMTHc2Nf1jJut3BPRT4XP4eahuag2oKZ1tOrNDlx4Ku1UUVVB5vsNk858iXxdYHqeqXROQWAFWtOd2KWaEgKT6WYX262Xz+MFWw/wjzlxbxz037aVblcyP6MGdqJmMHdD/pcdePTf/M2gzXj0nn9qdWc+v8fH57y7igbOPdUrA1dXBwFGy1Jykhli9m9+Ovq/fw7RlD6OmH4bS3CkroHBsV0DUHwpEvZ/z1IhKPdzEWERmE5xtASMrJTGX9nsPUNdo4fzhoblbe21HGV55cxed/s4y8ghJmThzIB3Mv4Xe3jvtM0m9L/5QEXr5nIkN6deXuv6zlhdV7HI68417fvJ+yqsCusHWu7pycSWOz+mWFruZmZXFBKRcNTiM+LvguaocSXxL/D4C3gP4i8hywBPi2o1E5KDcrhbrGZjbvO+J2KOYc1DY08eKaPVzx6w+585k1FJZX89DnLmDlw5fyg2uG0z+l490aU7t04q935XLR+Wk8tGgLv1nyMQ6sEXRWWgq2BvfsElJnuxk9Erl8qGeFrnMtntxcfISSo7U2m8cPfFmI5W0RWY9ncXQB7lfVkGrZ0NqETO84/6eHuDAjxe1wTAdVVNfz7Krd/HnlLg4eq2don248/sXRXD2qL3Ex536xL7FTDE/ens1Dr2zh8bc/ovRoLT++bgTRLs8XX1l4iG0HjvLoDSNDrkXBnKlZLN5WyqIN+7g15+xnIr21tYSYKOHSC4JvGC7U+LLm7kXemy1THoaJCKr6oXNhOSc5IY4hvbqSX1TB190OxvissPwYC5YV8cr6fdQ2NDNtSBp3T81i4qBUvyfC2OgofnnzKHp268Qf3v+Ug8fq+L8vj3V1zvyCpUWkJsaF5LrRF2Z0Z1S/JBYsK+KWCwecVdGVqrK4oITcrFSSEmIdiDKy+HJxd26r252BCXjW353uSEQBkJuVyotr9tLQ1GxTwoKYqrK6qIInlxaxZEcpsVFRfGFsOnOmZjLY4Y6MIsJ3rryAnl078ePXtzFzQT7zb7/QlaTzafkxluwo4/5LBwdlwVZ7RITZUzK5/4WNvP9RGdPP4oz9k7JjFB6s5s4Qur4RzHzp1XNNq+1yYARQ6nxozsnNSqGmocnG+YNUY1Mz/9y0n+t/t5wvzVvFut0VfP2S81j+0HR+ftMox5N+a3dOzuSJW8ayae8Rbv7TCg4cCXyF79Pegq3bgrhgqz1XjexDn6TOZz21862tniUWrwjC2VahyJcz/lPtw5P8Q9aETE/F36rCQ4wf6NusD+O8qtoGXlyzl6eX76K4sobMHon89PoR3Diun6uzOK4e1ZeUxDju/vM6bvj9Cv48a0LAPnwOV3sKtq4fG9wFW+1pKeh69M0dbNt/lGF9OzYdNW9bCWMHJIdUhXUw86Vy9wkR+Y13+y2wFNjkw/P6i8h7IrJdRApE5P5Tfv4tEVER6XH24Z+dlMQ4zu/VxRZgd8HpWiLvr6zhZ29sZ9Ij7/LTf20nPTmeeTPHs+SbF3Nb7sCgmLo3aVAPXvyPXBqblZv+uJK1uwLzu/PX1d6CrSlZATmek265cAAJcR0v6Np3+Dhbi49ypRVt+Y0vZ/xrW91uBJ5X1eU+PK8ReFBV14tIV2CdiLytqttEpD9wOeDaZOmczFQWrd9n4/wB1NISuaU7ZnFlDQ++tIlmVaKihM+N6M1dU7MYHaStdof3TWLRvZO4w1vo9ZtbxjpaQVrX2MQz3oKtIb1Df5UpT0FXf57L3813rvS9oCuvwDOybNW6/uPLGP/CVttzPiZ9VPWAqq733q4CtgMtUxL+F08tgGuTpHOzUqmub2JrsY3zB8rpWiI3qZLQKZoP5k7jt18ZF7RJv0X/lAT+ds9ELujTjXufXcdf8507d3l90wHKq+qYMzX0z/Zb3Dk5g8Zm5c8dWKErr6CEIb26ktEj0cHIIosvQz1bRGTzabYtIrLZl4OISAYwFsj3LuNYrKpnHC4SkbtFZK2IrC0v939r1wmZnjn8NtwTOG21Pj5e10S/7h0vuHJLapdOPH9XDhefn8Z3X93Cr9/5yO+FXqrK/GWegq2LBgd8NNQxA1MTuWJYL57N962g6+CxOtbuqmCGFW35lS9jHG/iqdy91bu9AbwMXA1c096TRaQL8ArwAJ7hn+8B32/veao6T1WzVTU7Lc3/lYppXTsxKC2RfGvYFjB9k+M7tD+YJcTFMO/2bG4a349fv/Mx3311K41NzX57/ZWferpazpmaGXIFW+2ZMzWLyuMNvLJ+X7uPfWdbKc0KM4bbbB5/8iXxT1bVb6vqFu/2EDBDVXer6hm/r4lILJ6k/5yqLgIGAZnAJhHZBfQD1ouIKx/nOVmprN112K9/sKZts6dkfGZffGw0c2cMCXwwfhAbHcVjN43ia5cM4vnVe7j3ufV+W+th/jJPwdZ1Y0KvYKs92QO7M7pfEk8tK2p3ha68ghL6dY9nWJ/gb0oXSnxJ/IkiMqXljohMAtodbPN28FwAbFfVxwG8Hxw9VTVDVTPwTA0dp6olZxX9OcrNSqWqrpFtfuwXbtp24EgtAvTq1ulES+RHbhgZktWoLUSEuTMu4EfXDued7aXcNj+fyuP15/San5Yf490dwbvC1rkSEWZNyaTwYDXvf1TW5uOqahtY/skhrhzeO+y+9bjNl1k9s4GnRCTJe78SmOXD8yYDM4EtIrLRu++7qvpGR4N0Sm7LOH9hBaP6JbsbTJg7VtfIC6v3cvXovjxxy1i3w/G7OyZl0KNLJ77x4kZu/uNKFs6acNZDWE8tKyIuJrQLttpz1cg+PPrmDuYvLWqzkve9neXUNzXb+L4DfJnVs05VRwOjgNGqOqZltk47z1umqqKqo7zPGXNq0vee+bvW8K1nt85k9kgkv8jG+Z320pq9VNU1hlRL4Y76/Kg+LJw1gZIjtdzw+xV8VNrxFb0qqut5Zf0+vjAmPaQLttoTGx3FVydlsOLTQxTsP/3MuryCEnp0iWOcj621je/aTPwicpv332+KyDeBOcDsVvfDQk5mCquLKmhqZ6zRnL2mZuWp5UVcmNGdMUE+XfNcTRyUykv3TKRZlZv+sILVHZw19tf83SGzwta5+vKEtgu6ahuaeH9HGZcP6+16Z9RwdKYz/pZx/K5tbGEhNyuVo7WNfl0X1JxscUEJ+w7XhEX1qS+G9unGovsm0aNrJ25bkH+iz0x76hqbWLhyNxedn8b5AexH5JakeE9B1z837af0aO1JP1v+yUGq65tsNo9D2kz8qvon778/Ot0WuBCdlZNl8/mdNn9ZEQNSErg8ghps9euewMv3TGJYn27c99w6nvVhBap/thRshfFw2KlmeVfo+vPKXSftzysooWunGCYNCp8ahmDiSwFXmoh8V0TmichTLVsggguEPknxDEhJsPn8Dlm/5zDrdh9m1uSMiPvKnpIYx1/vymHakJ78v9e28vjinW0WenlW2Crk/F5dmBpGBVvtGZCawIxhvXkuf8+Jgq7Gpmbe2V7G9KE9/bK4jvksX97VvwNJwDvAv1ptYSMnM4XVuyranVNsOm7BsiK6do7h5uz+bofiioS4GObNHM8Xs/vxm3c/4eFFW05bN7Li00PsKKli9pTwK9hqz+ypmScVdK3ZdZiK6nrrzeMgX6ZzJqjqdxyPxEW5Wan8bd0+dpZWMdQKRfxm3+HjvLnlAHdNzSKx09l0AA8PMdFR/PzGUfTs2pnfvvcJB4/V8cQt407qOrogjAu22tO6oOsrEwaQV1BCXExUSK0tHGp8OeN/XUSucjwSF50Y57fhHr9auGIXIsIdkzLcDsV1IsK3ZgzhJ9cNZ8mOMm6dv4rD1Z5Cr0/KPAVbMyeGZ8FWe0SE2VOzKDxYzbs7ylhcUMJFg9Mi+mTBab68s/cD3xWROqABz4Lrqqphc2rcr3sC6cnx5BdV8NXJkXNhzUlVtQ28sHovnx/ZJyR78Thl5kRPodf9L2xkxq8/RARKj9YBkNolzuXo3PO5Eb1Jjo/hP/6yjiZVjjc08dqG4pCu6g5m7SZ+VQ3/eWV4hnve21mGqkbcGKsTXlq7j6q6RuZEwHz0jvrcyD4U7D/Cb9/79KT9P/vXDrp2io3IZPevzQc4VtdEk/fid+XxBh5etAUgIt8Pp/kyq+ei022BCC6QcrJSqKiu5+OyY26HEvIam5p5enkREzJSrBVGG17dsP8z+2oamngsb6cL0bjvsbydNJ4yuSKS3w+n+TLUM7fV7c7ABGAdMN2RiFyS612HN7/wUEQUzzhp8bZS9h2u4b+vHuZ2KEGrrbUJ2tof7uz9CCxfevVc02q7HM9C66XOhxZY/VPi6ZPUmVVWyHXO5i8tZGBqApcNjZyCrY4Kp7UJ/MHej8A6m+qIfXiSf1gREXKzUskvPOT31ZQiybrdh1m/p5I7J0VewVZHzJ0xhPhTZvCE8toE58rej8Bqd6hHRJ7g32vjRgFjgDMumxiqcjJTeHVDMZ+WV3Nezy5uhxOSnorwgi1ftVywfCxvJ/sra+ibHM/cGUMi9kKmvR+B5csY/9pWtxuB531dcD3U5GR5x/mLDlniPwt7K47z5tYD3HVRZBds+er6semW2Fqx9yNwfJnOuTAQgQSDjNQEenbtRH5hBbfmhO8iGE5ZuGIXUSJ81Qq2jAlq1gGplZZx/lU2zt9hVbUNvLBmL58f1Yc+SXZBzphgZon/FDlZKZRV1bHr0HG3QwkpL67Zy7EwX2HLmHBxphW4/uL99/7AheO+nFbz+Y1vPAVbu5iQaQVbxoSCM53xjxeRgcAsEekuIimtt0AFGGiD0hLp0aWTLczSAXkFpRRX1kTUAiLGhLIzXdz9I/AWkIWnUrf1pGz17g87IkJOVsqJcX7r29O++cs8BVuXWsGWMSHhTEsv/kZVhwJPqWqWqma22sIy6bfIzUzhwJFa9lZYuXh71u0+zIY9lcyanGkFW8aECF9aNtwrIqNF5D+92yhfXlhE+ovIeyKyXUQKWq4ViMhjIrJDRDaLyKsiknyO/w1+1zKff1WRjfO356llRXTrHMNN4/u5HYoxxke+dOf8L+A5oKd3e05Evu7DazcCD3q/NeQCXxORYcDbwAhVHQV8BDx8tsE7ZXDPLqQkxpFfaOP8Z9JSsPWVnIFWsGVMCPHlr3UOkKOq1QAi8nNgJfDEmZ6kqgeAA97bVSKyHUhX1cWtHrYKuOlsAneSiJCT6RnnN217xluwdcckK3YzJpT4Mo9fgKZW95s4+UJv+y8gkgGMBfJP+dEs4M02nnO3iKwVkbXl5eUdOZxf5GSmUFxZw77DNp//dI7WNvDimr1cbQVbxoQcXxL/00C+iPxQRH6I5yx9ga8HEJEuwCvAA6p6tNX+7+EZDnrudM9T1Xmqmq2q2WlpgV90+UTfHhvuOa2XThRshfV1fmPCki8Xdx8H7gQqgMPAnar6a19eXERi8ST951R1Uav9dwBXA7dqkPZGGNKrK8kJseTbBd7PaCnYyslMYWS/JLfDMcZ0kE9X5FR1PbC+Iy8sngnwC4Dt3g+Plv1XAt8BLlbVoB1HiYoSJmSksMrO+D/jrYISiitr+ME1tsKWMaHIyV49k4GZwHQR2ejdrgJ+C3QF3vbu+6ODMZyTnKxU9lQc58ARm8/f2oJlRWRYwZYxIcuxOXiquozTXwR+w6lj+ltOpqczRX5hhfUJ92op2PrxdcOtYMuYEHXGM34RiRaRdwIVTLAZ2qcbXTvH2LTOVhYsKyQpPtYKtowJYWdM/KraBBwXkYi8ghcd5ZnPbw3bPPZWHOetrSV8JWcACXFWsGVMqPLlr7cW2CIibwPVLTtV9b8ciyqI5GSm8s72MsqO1tKzW2e3w3HV08u9BVsTM9wOxRhzDnxJ/P/ybhEpJ8szzr+qqIJrR/d1ORr3eAq29nDN6L70TorsD0BjQp1Pa+6KSDwwQFV3BiCmoDKsTze6dvKM80dy4n9x9V6q65tshS1jwoAvTdquATbi6c2PiIwRkX84HFfQiImOIjuje0SvyNXY1MwzKzwFWyPSI/JyjzFhxZd5/D8EJgCVAKq6EYio076crFQ+La+mvKrO7VBc0VKwNWeqtWcwJhz4kvgbVfXIKfuCss2CU1rm86+OwNk9qsqTS70FWxf0dDscY4wf+JL4t4rIV4BoERksIk8AKxyOK6iMSE8iMS46Iufzr99zmE17K5k9JZMoK9gyJiz4kvi/DgwH6oDngaPAAw7GFHRio6MYn5ESkQ3b5i8tIik+lhutYMuYsOFLd87jqvo94FLgElX9nqrWOh9acMnJTOGj0mNUVNe7HUrA7Dl0nLyCEm61gi1jwoovs3ouFJEtwGY8hVybRGS886EFl9yslnH+yDnrf3pFEdFRwh2TMtwOxRjjR74M9SwA7lPVDFXNAL6GZ3GWiDIyPZn42OiIadN8pKaBl9bs5ZpRfekV4RXLxoQbXxJ/laoubbnj7bpZ5VxIwSkuJorxA7tHzAXeF9fsobq+iVlWsGVM2Gkz8YvIOBEZB6wWkT+JyDQRuVhEfg+8H7AIg0hOZgo7S6uoPB7e4/yNTc08s3wXuVlWsGVMODrTFbtfnXL/B61uR9Q8/hY5WamoeubzXzG8t9vhOObNrSXsP1LLj68b4XYoxhgHtJn4VfWSQAYSCkb3T6JTTBSrCsM38asq85cWktkjkelWsGVMWGp3jp6IJAO3AxmtHx8pbZlb6xQTzbgB3cN6Pv+63YfZtO8IP7l+hBVsGROmfLm4+waepL8FWNdqi0g5WSlsO3CUIzUNbofiiPlLi0hOiOXGcbbUpDHhypeqnM6q+k3HIwkROZmpqH7M2l0VYbfY+O5D1eRtK+G+aYOsYMuYMObLGf9fROQuEekjIiktm+ORBamxA5KJi4kKy2mdTy/fRUyUcLutsGVMWPMl8dcDjwEr+fcwz9r2niQi/UXkPRHZLiIFInK/d3+KiLwtIh97/+1+Lv8BgdY5Npox/ZPDbh3eIzUNvLR2L9eMtoItY8KdL4n/m8B53srdTO/mS2P2RuBBVR0K5AJfE5FhwEPAElUdDCzx3g8puZkpbC0+QlVt+Izzv7hmD8dthS1jIoIvib8AON7RF1bVA6q63nu7CtgOpAPXAQu9D1sIXN/R13ZbTlYqzQprdx12OxS/aPAWbE3MSmV4XyvYMibc+XIFrwnYKCLv4WnNDHRsOqeIZABjgXygl6oe8L7GARE57WRxEbkbuBtgwIABvh4qIMYN6E5stLCq6BCXhMFc95aCrZ9cbwVbxkQCXxL/a97trIhIF+AV4AFVPSri29xwVZ0HzAPIzs4Oqkrh+LhoRvdLJj8MGra1FGxl9UjkkiGh/yFmjGlfu4lfVRe295i2iEgsnqT/nKou8u4uFZE+3rP9PkDZ2b6+m3KyUvjjB4VU1zWS2Cl0pz6u3X2YzfuO8FMr2DImYvjSj79IRApP3Xx4nuBp6bxdVR9v9aN/AHd4b98B/P1sAndbTmYqTc3K2t2hPc4/f2mht2DLVtgyJlL4cqqa3ep2Z+BmwJd5/JOBmXgWb9no3fdd4FHgJRGZDezxvl7IGT+wOzFRQn7hIS4+P83tcM7K7kPVLN5WytemnUd8XLTb4RhjAsSXoZ5TK5V+LSLLgO+387xlQFtjB5f6Fl7wSuwUw8h+SSE9n//fBVsD3Q7FGBNAvjRpG9fqbhSebwBdHYsohORkprJgWSHH6xtDrsVB64KtnlawZUxE8SVbte7L3wjsAr7oSDQhxnOB91PW765kyuAebofTIS+stoItYyKVL0M91pe/DdkDuxMdJeQXHQqpxN/Q1MwzK3YxaZAVbBkTiXwZ6ukE3Mhn+/H/2LmwQkPXzrGM6Nst5Obzv7HlAAeO1PI/X7CCLWMikS8tG/6Op81CI1DdajN42jds3FtJbUOT26H4RFVZsKyIrLREpp1vBVvGRCJfxvj7qeqVjkcSonKzUpj3YSHr9xxm0qDgH+5Zs8tTsPU/X7CCLWMilS9n/CtEZKTjkYSo7IwUooSQGe6Zv7SQ7gmx3DDWCraMiVS+JP4pwDoR2Skim0Vki4hsdjqwUNGtcyzD+nYL+nV4X9tQTM7P3mHxtlIam5S8ghK3QzLGuMSXoZ7POR5FiMvJTOXZVbupbWiic2zwVcC+tqGYhxdtocZ7HaKqrpGHF20B4PqxtrauMZGm3TN+Vd19ui0QwYWK3KxU6hqb2bS30u1QTuuxvJ0nkn6LmoYmHsvb6VJExhg3+TLUY9oxISMFEYK2fUNxZc1p9+9vY78xJrxZ4veDpIRYLugdfOP8zc3KI29ub/PnfZPjAxiNMSZYWOL3k5zMFNbtPkx9Y7PboQCe6txv/W0Tf/qgkMmDUomPPfl/dXxsNHNnDHEpOmOMmyzx+0luViq1Dc1s3lfpdihU1zUye+FaFm0o5ltXnM+zc3J45IZRpCfHI0B6cjyP3DDSLuwaE6FCq6VkEJuQ6VmiIL+oguwMX5YrcMbBY3XMemYNBfuP8vMbR/KlCz3rFV8/Nt0SvTEGsDN+v0lJjGNIr66sKnRvnH/PoePc9IcVfFRaxbyZ408kfWOMac0Svx/lZHnG+RuaAj/Ov7X4CDf8YTmVNQ08NyeXS4f2CngMxpjQYInfj3KzUjle38SW4iMBPe6yjw/ypT+tpFNMNC/fM4nxA7sH9PjGmNBiid+PTozzB7Bvz983FnPnM6vpn5LAovsmcV7PLgE7tjEmNFni96MeXTpxXs8uAZvPP39pIfe/sJFxA7rz4n9MpJctoWiM8YElfj/LyUxhTVEFjQ6O8zc3Kz97Yzs//dd2rhrZm4WzJpAUH+vY8Ywx4cUSv5/lZqVSXd9Ewf6jjrx+fWMz33xpI/M+LOT2iQN54pZxQdkYzhgTvBxL/CLylIiUicjWVvvGiMgqEdkoImtFZIJTx3dLTlbLfH7/D/ccq2tk9sI1vLZxP3NnDOFH1w4n2hZTMcZ0kJNn/M8Ap67c9QvgR6o6Bvi+935Y6dm1M1k9Ev1+gbe8qo5b5q1ixaeH+MVNo/jaJechYknfGNNxjiV+Vf0QODX7KdDNezsJ2O/U8d2Uk5XC6qIKmprVL6+3+1A1N/1xBR+XVfHk7eP5YnZ/v7yuMSYyBXqM/wHgMRHZC/wSeLitB4rI3d7hoLXl5eWBis8vcrNSqaprZPuBcx/n37LvCDf+YQVHaxp4/q5cpl9ghVnGmHMT6MR/L/ANVe0PfANY0NYDVXWeqmaranZaWlrAAvSHnMxUgHNu37D043K+PM9bmHXvJMYOsMIsY8y5C3TivwNY5L39NyDsLu4C9E7qzMDUhHNamOW1DcXc+fQaBqQmsui+SQxKs8IsY4x/BDrx7wcu9t6eDnwc4OMHTE6mZ5y/+SzG+Z/8sJAHXtzIhRkpvPgfuVaYZYzxKyencz4PrASGiMg+EZkN3AX8SkQ2AT8D7nbq+G7LzUrlSE0DO0qqfH5Oc7Py09e38T9vbOfzo/rwzKwL6dbZCrOMMf7lWD9+Vb2ljR+Nd+qYwSQnyzPOn190iGF9u7XzaE9h1rf+tol/bNrPVydl8P2rhxFlc/SNMQ6wyl2HpCfH0697vE8XeI/VNTLrmTX8Y9N+vnPlBfzgGkv6xhjn2ApcDsrNSmXJ9lKam7XNRF5eVcedz6xm+4EqfnnzaG4a3y/AURpjIo2d8TsoJzOFw8cb+Ljs2Gl/XnSwmhv/sIJPy6qZf0e2JX1jTEBY4ndQbqtx/lNt3lfJTX9YwbG6Rp6/O5dLhvQMdHjGmAhlid9B/brH0zep82fG+T/4qJwvz1tFfFw0L98zkTH9k90J0BgTkSzxO0hEyM1KZXVRBaqe+fyL1u9j9jNryEhNZNG9k8iywixjTIDZxV2HxcYIB4/Vk/XwG3TtHMPR2kYmDUrlTzPH09Xm6BtjXGBn/A56bUMxr23wNCBV4GhtI9ECN4xNt6RvjHGNJX4HPZa3k7rGk5dgbFL433fCtlOFMSYEWOJ30P7Kmg7tN8aYQLDE76C+yfEd2m+MMYFgid9Bc2cMIf6UhdDjY6OZO2OISxEZY4zN6nHU9WPTAc9Y//7KGvomxzN3xpAT+40xxg2W+B12/dh0S/TGmKBiQz3GGBNhLPEbY0yEscRvjDERxhK/McZEGEv8xhgTYaSla2QwE5FyYPdZPr0HcNCP4YQ6ez/+zd6Lk9n7cbJweD8GqmraqTtDIvGfCxFZq6rZbscRLOz9+Dd7L05m78fJwvn9sKEeY4yJMJb4jTEmwkRC4p/ndgBBxt6Pf7P34mT2fpwsbN+PsB/jN8YYc7JIOOM3xhjTiiV+Y4yJMGGd+EXkShHZKSKfiMhDbsfjFhHpLyLvich2ESkQkfvdjikYiEi0iGwQkdfdjsVtIpIsIi+LyA7v78lEt2Nyi4h8w/t3slVEnheRzm7H5G9hm/hFJBr4HfA5YBhwi4gMczcq1zQCD6rqUCAX+FoEvxet3Q9sdzuIIPF/wFuqegEwmgh9X0QkHfgvIFtVRwDRwJfdjcr/wjbxAxOAT1S1UFXrgReA61yOyRWqekBV13tvV+H5o47oRQJEpB/weWC+27G4TUS6ARcBCwBUtV5VK10Nyl0xQLyIxAAJwH6X4/G7cE786cDeVvf3EeHJDkBEMoCxQL7Lobjt18C3gWaX4wgGWUA58LR36Gu+iCS6HZQbVLUY+CWwBzgAHFHVxe5G5X/hnPjlNPsieu6qiHQBXgEeUNWjbsfjFhG5GihT1XVuxxIkYoBxwB9UdSxQDUTkNTER6Y5nZCAT6Askisht7kblf+Gc+PcB/Vvd70cYfmXzlYjE4kn6z6nqIrfjcdlk4FoR2YVnCHC6iDzrbkiu2gfsU9WWb4Ev4/kgiESXAUWqWq6qDcAiYJLLMfldOCf+NcBgEckUkTg8F2j+4XJMrhARwTN+u11VH3c7Hrep6sOq2k9VM/D8XryrqmF3VucrVS0B9orIEO+uS4FtLobkpj1ArogkeP9uLiUML3SH7WLrqtooIv8J5OG5Mv+Uqha4HJZbJgMzgS0istG777uq+oZ7IZkg83XgOe9JUiFwp8vxuEJV80XkZWA9ntlwGwjD1g3WssEYYyJMOA/1GGOMOQ1L/MYYE2Es8RtjTISxxG+MMRHGEr8xxkQYS/wmbIlIqohs9G4lIlLsvX1MRH4foBimtXT/FJFrI7lLrAkeYTuP3xhVPQSMARCRHwLHVPWXLsbzDyK0iNAEFzvjNxHnlLPwH4rIQhFZLCK7ROQGEfmFiGwRkbe8rS4QkfEi8oGIrBORPBHpc5rXvdnbw32TiHx4mp9/VUR+673dS0Re9T52k4hM8u6/TURWe7+Z/MnbXtwYv7LEbwwMwtOi+TrgWeA9VR0J1ACf9yb/J4CbVHU88BTwP6d5ne8DM1R1NHBtO8f8DfCB97HjgAIRGQp8CZisqmOAJuDWc/2PM+ZUNtRjDLypqg0isgVPe4+3vPu3ABnAEGAE8LanfQvReFr2nmo58IyIvISnudeZTAduB1DVJuCIiMwExgNrvMeJB8rO/j/LmNOzxG8M1AGoarOINOi/+5g04/kbEaBAVc+4HKGq3iMiOXi+PWwUkTEdjEOAhar6cAefZ0yH2FCPMe3bCaS1rEMrIrEiMvzUB4nIIFXNV9XvAwc5uS34qZYA93qfF+1dBWsJcJOI9PTuTxGRgX7+bzHGEr8x7fEu3XkT8HMR2QRs5PQ92h/zXhTeCnwIbDrDy94PXOIdXloHDFfVbcD/AxaLyGbgbeAzF5GNOVfWndMYYyKMnfEbY0yEscRvjDERxhK/McZEGEv8xhgTYSzxG2NMhLHEb4wxEcYSvzHGRJj/DwpAsxR02iHSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(len(q6clusters)), [len(x) for x in q6clusters], 'o-')\n",
    "plt.xlabel(\"Time slice\")\n",
    "plt.ylabel(\"number of unique clusters\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also plot the biggest cluster size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Largest cluster size')"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4R0lEQVR4nO3dd3hU95Xw8e/RqEsDAo2QkOio2GBMMXFBxImdgo3rOs3e9ORN1m96WdYmySa2N8UJqU6ctk7i7Juss05is7FNjEuKY3AcgwSmCgRGAgmEBJJQb3PeP2YEQh5JI2lm7pTzeZ55NHPn3rlHA5oz91fOT1QVY4wxZrgkpwMwxhgTnSxBGGOMCcgShDHGmIAsQRhjjAnIEoQxxpiAkp0OIJQ8Ho/OmzfP6TCMMSZmbN++vUlV8wI9F1cJYt68eWzbts3pMIwxJmaISM1Iz1kTkzHGmIAsQRhjjAnIEoQxxpiALEEYY4wJyBKEMcaYgOJqFJMxobSxso4Nm6uob+miMCeDdWvKuHl5kdNhGRMxliCMCWBjZR3rH9lFV98AAHUtXax/ZBeAJQmTMKyJyZgANmyuOpscBnX1DbBhc5VDERkTeZYgjAmgvqVrXNuNiUeWIIwJoDAnY1zbjYlHliCMCWDdmjLSks//88hIcbFuTZlDERkTeWFLECIyW0T+LCL7RGSPiHzSv32DiOwXkZdF5FERyRnh+CMisktEdoiIFVgyEXXz8iLWXlRw9nFORgpfu2WJdVCbhBLOK4h+4LOqeiFwOfBREVkEPA1cpKoXAweA9aO8xlWqukxVV4YxTmMC6vMqBVPSSU9J4i2XzLLkYBJO2BKEqh5X1Qr//TZgH1Ckqk+par9/t78Ds8IVgzGTUVnbwiXzplEyw82BhjanwzEm4iLSByEi84DlwIvDnvoA8McRDlPgKRHZLiIfHuW1Pywi20RkW2NjY0jiNeZEazd1LV1cMmcapfmWIExiCnuCEJFs4PfAp1T1zJDtn8fXDPXrEQ4tV9UVwLX4mqeuDLSTqv5UVVeq6sq8vIBrXhgzbhW1zQCsmDuNsoJsGs700NLZ63BUxkRWWBOEiKTgSw6/VtVHhmx/L3A98E5V1UDHqmq9/+dJ4FHg0nDGasxQFTXNpCUnsWjmFEry3QAcaGh3OCpjIiuco5gE+BmwT1W/PWT7NcAdwI2q2jnCsVki4h68D7wZ2B2uWI0ZrqK2mSVFU0lNTqLMnyCqrJnJJJhwXkGUA+8GrvYPVd0hImuBHwBu4Gn/th8DiEihiGzyH5sPPC8iO4F/AE+o6pNhjNWYs3r6B9hdd4ZL5k4DYObUdNxpyRy0BGESTNiK9anq84AEeGpTgG2DTUpr/fcPA0vDFZsxo9ldd4beAS/L5/gShIhQWuCm6oQlCJNYbCa1McNUnu2gzjm7bXAk0whdZsbEJUsQxgyzvaaZ2dMzmOFOP7utND+b5s4+Gtt7HIzMmMiyBGHMEKpKRW0zK/zNS4MGO6oPnLCRTCZxWIIwZoj61m4azvS8KkGUFgwOdbV+CJM4LEEYM0RFjb//YViC8GSnkZuVagnCJBRLEMYMsb2mmYwUFxfMdL/quZL8bJsLYRKKJQhjhqisbebiWVNJcb36T6Ms383BhnYbyWQShiUIY/y6+wbYU3+GFXOnBXy+tMBNe08/9a3dEY7MGGeEbaKciT0bK+vYsLmK+pYuCnMyWLemLKHWQHj5WCv9Xn1V/8OgcyOZ2iiypUdNArArCAP4ksP6R3ZR19KFAnUtXax/ZBcbK+ucDi1izlZwnZMT8PkSq8lkEowlCAPAhs1VdPUNnLetq2+ADZurHIoo8ipqmpmXm0ludlrA56dmpFAwJZ0DVnLDJAhLEAaA+paucW2PN74Jci0jNi8NKi1w2xWESRiWIAwAhSO0qY+0Pd4cPd1FU3vPiB3Ug8rys6k+2c6A10YymfhnCcIAsG5NGa6k84vvZqS4WLemzKGIIutc/8MYVxD5bnr6vdSeDriUiTFxxRKEAWBhXjYDXiU7zXV22z03LU6YUUwVtc1kpbooK3j1BLmhBp+30t8mEViCMKgqdz+2B092KlvXv4FfvO81QOI0L4EvQSydnfOqq6jhimdkA1aTySQGSxCGP+ysZ1tNM+vWlDElPYVL508nOUl4vrrJ6dAiorO3n33H286uIDeazNRk5kzPtI5qkxAsQSS4zt5+vrZpPxcVTeGtl8wGICstmeVzctiSIAli59FWBkaZIDdcab7blh81CcESRIL74Z8PceJMN3fdsPi85pXyYg+76lpp7exzMLrIGOygXj7CBLnhygqyOdzYQW+/N4xRGeM8SxAJrPZUJz/922FuWlbIynnTz3uuvNiDKrxwOP6vIipqmlmQl0VOZmpQ+5fmu+n3Kq80dYQ5MmOcFbYEISKzReTPIrJPRPaIyCf926eLyNMictD/M+B1vYhcIyJVIlItIneGK85E9pVNe3GJcOe1F7zquWWzc8hKdcV9P4SqUnm0hUuCbF4CX4IAK7lh4l84ryD6gc+q6oXA5cBHRWQRcCfwrKqWAM/6H59HRFzA/cC1wCLgNv+xJkS2VDexeU8DH71qITOnvnq0UooricsW5LK1+pQD0UXOkVOdnO7oHXOC3FAL8rJwJYmV3DBxL2wJQlWPq2qF/34bsA8oAm4Cfunf7ZfAzQEOvxSoVtXDqtoL/MZ/nAmB/gEvdz+2h9nTM/g/r10w4n6rFuZyuKmDujgutzHSCnKjSUt2Md+TZUNdTdyLSB+EiMwDlgMvAvmqehx8SQSYEeCQIuDokMfH/NsCvfaHRWSbiGxrbGwMadzx6tcv1nKgoZ3Pr11EeoprxP1Wl3gA4no00/baZtxpyZT45zcEqyzfbQnCxL2wJwgRyQZ+D3xKVc8Ee1iAbQGL36jqT1V1paquzMvLm2iYCaO5o5dvP32A8uJc1izOH3Xfsnw3nuzUuE4QFTXNLJuTQ9IYE+SGK813U3O6k67egbF3NiZGhTVBiEgKvuTwa1V9xL+5QURm+p+fCZwMcOgxYPaQx7OA+nDGmii+9XQV7T39fOmGxYiM/qEoIpQXe9hSfSoul9ls7+nnQEPbuJqXBpXmZ6MK1SfbwxCZMdEhnKOYBPgZsE9Vvz3kqT8A7/Xffy/wvwEOfwkoEZH5IpIK3Oo/zkzC3voz/PeLtbz78rlnR+KMpXyhh6b2Hg40xN8H4c6jLXiVcXVQDyotsJFMJv6F8wqiHHg3cLWI7PDf1gL3Am8SkYPAm/yPEZFCEdkEoKr9wMeAzfg6tx9W1T1hjDXuDdZbmpqRwqffWBr0ceX+foh4HO66vaYZEd+Q3vGaOz2T1OQkm1Ft4lrY1qRW1ecJ3JcA8IYA+9cDa4c83gRsCk90iWfTrhO8+MppvnzzRUzNTAn6uKKcDOZ7stha3cQHV88PY4SRV1HbTMmMbKZmBP9+DEp2JVGcl21XECau2UzqBNDVO8BXN+3jggI3t106Z9zHr1qYy98Pn6JvIH5KS3i9SmUQK8iNpqzAbXMhTFyzBJEAfvLcIepaurjrxsVjlrMOZHWxh47eAXYebQl9cA453NROa1ffpBJESX429a3dnOmO/3pVJjFZgohzdS1d/Pivh7huyUwuX5A7ode4YmEuIrAljmZVV9S0ABProB5U5u/ot34IE68sQcS5r27ahyqsX/vqekvByslM5aLCqXE1H6KitpmpGSks8GRN+DUGR4LF4wgvY8ASRFx78fApnnj5OLe/biGzpmVO6rXKiz1U1DbT0dMfouicVVHbzPIJTJAbqigng6xUly0/auKWJYg4NeBV7npsL4VT07n9dQsn/Xqriz30e5V/HDkdguic1drVx4GG9nFVcA0kKUkosZIbJo5ZgohTv3mpln3Hz/C56y4kI3XkekvBWjlvGqnJSWw5GPvNTDv8ne2T6X8YVJqfbQnCxC1LEHGotbOPb26u4rL507luycyQvGZ6iouVc6fFxYS5ippmkgSWTmCC3HCl+W6a2ns51d4z+cCMiTKWIOLQd545QGtXX1D1lsajvNjD/hNtNMX4h2FFbTOl+W6y0yY/T7SswDqqTfyyBBFnDjS08f/+XsNtl85hUeGUkL52ebGv7MbWQ7E73NXrVXbUtnBJCJqX4NxQV2tmMvHIEkQcUVXueWwvWakuPvvmspC//pKiqbjTk2O6H+LgyXbaevonNUFuqDx3GjmZKVZyw8SloBKEiKwWkff77+eJSHwV5YkTT+1t4PnqJj7zplKmZ6WG/PVdScKqhbk8X90Us+W/K2r9K8iF6ApCRCidYSU3THwaM0GIyJeAO4D1/k0pwK/CGZQZv+6+Ab78xF5K87N51+Vzw3ae1cUe6lq6qD3dGbZzhNP2mmamZ6UyL3dy80KGKi3wFe2L1aRpzEiCuYL4J+BGoAPOVl0NbjEBEzE/e/4Vjp7u4ovXLybZFb6Ww1XFsV3+u6K2mRVzckLaeV+W76atu5+GM7HdeW/McMF8kvSq76uRAojIxGsTmLA40drN/X+u5s2L8s+uIx0uCzxZzJyaztYYrMvU0tnL4cYOloeo/2HQYMkN64cw8SaYBPGwiPwEyBGRDwHPAA+ENywzHl9/cj/9XuUL1y0K+7nOLkN6qAmvN7aaVCprWwBC1kE96GxNJuuHMHFmzAShqt8Efodvbeky4Iuqel+4AzPB2V7TzKOVdXzotfOZE8J29dGUF+fS0tnH3uNnInK+UNle04wrSVg6e2pIX3daVip57jS7gjBxJ5hO6n8H9qvqOlX9V1V9WkQ+HIHYzBi8Xt8yovlT0vjI64sjdt7yhbHZD1FR28yFM91kpoZ+IcUyq8lk4lAwTUwfBzaLyFVDtt0epnjMOPxu+zFePtbK+msvJCsEs4KDNWNKOqX52TFV/nvAq+w8OrkV5EZTmu/mYEN7zDW7GTOaYBJEHXANcK+IrPNvC90QEDMhZ7r7+Mbm/Vwydxo3LSuM+PlXLfTw0pHTdPcNRPzcE1F1oo2O3oGwJYiygmy6+gY41twVltc3xglBjYdU1VrgdcAiEfktkDHWMSLycxE5KSK7h2z7HxHZ4b8dEZEdIxx7RER2+ffbFtyvkli+/+xBTnX0cleI6y0Fa3Wxh+4+79mJZ9Fuuz/OUJXYGK7ERjKZOBRMgtgGoKrdqvp+4C9AMNN0H8R35XGWqr5DVZep6jJ8nd6PjHL8Vf59VwZxroRyqLGdX2w5wtsvmc2SWaHtcA3WZQum40qSmBnuWlnTjCc7jVnTxvxuMyElM7IBq8lk4kswo5g+NOzx/aq6IIjjngMCri4jvq+8bwceCjJOM8R/PL6XjBQX664Jfb2lYLnTU1g6a2rMdFSHY4LcUO70FIpyMixBmLgyYoIQkYf9P3eJyMvDb5M872uBBlU9OMLzCjwlItvHGjElIh8WkW0isq2xsXGSYUW/P+1v4C9VjXzyjSV4stMcjWV1sYeXj7XQ2tXnaBxjaWrv4cipzpDVXxpJWYHblh81cWW0K4hP+n9eD9wQ4DYZtzH61UO5qq4ArgU+KiJXjrSjqv5UVVeq6sq8vLxJhhXdevu9/Mfj+1iQl8V7rpjndDiUF3vwqm/t62g2OEEuXP0Pg0rz3Rxu7KBvwBvW8xgTKSMmCFU97r/bBBxV1RogDVgK1E/0hCKSDNwC/M8o5673/zwJPApcOtHzxZNfbHmFV5o6+PfrF5Ga7Hyl9uVzppGR4or64a4Vtc0kJwlLisLbX1Oan03vgJeaUx1hPY8xkRLMp8xzQLqIFAHPAu/H1wE9UW/EN/HuWKAnRSRLRNyD94E3A7sD7ZtITrZ18/0/VXP1BTO4qmyG0+EAkJqcxKXzp0d9P0RFTTOLC6eQnjL5tblHc7Ym0wlbXc7Eh2AShKhqJ75v/d9X1X8Cxiz6IyIPAS8AZSJyTEQ+6H/qVoY1L4lIoYhs8j/MB54XkZ3AP4AnVPXJ4H6d+LXhySp6+gf49+vDX29pPFYXezjU2MGJ1m6nQwmob8DLzmMtIS/QF0jxjGySxEYymfgRzPRbEZErgHcCgx/yYx6nqreNsP19AbbVA2v99w/ja8YyfjuPtvDb7cf4lysXMN8TXcV0VxXnArCluom3XDLL4Whebf/xNrr7vGHvfwBIT3ExLzfLEoSJG8FcQXwS32JBj6rqHhFZAPw5vGGZQV6vctdje/Bkp/GxqyNXbylYFxZMYXpWatT2Q4R6BbmxlOa7bbKciRvBXAk8h68fYvDxYeAT4QzKnLNxRx2VtS1seOvFuNNTnA7nVZL8y5BuOeRbhtSJWd2j2V7TTP6UNAqnpkfkfKX52Ty19wTdfQNh7/MwJtycHwpjRtTe08+9f9zP0tk5vGVF9DXfDFpd7KHhTA+HGqOvc7aitplL5k6LWOIqLXDjVaLyvTBmvCxBRLH7/1zNybYe7rphEUlJ0fXNfKjywWVID0ZXM9PJtm6ONXeFrUBfIGX+kUwHGyxBmNg3aoIQEZeIfDpSwZhzjjR18LO/vcItK4oiMgJnMmZPz2TO9Ey2HIquCXMVNS0AEX3/5nmySHGJ9UOYuDBqglDVAeCmCMVihvjyE/tIcQl3XnOB06EEpbzYw98PnaI/imYRV9Q2k+pK4qKiKRE7Z4oriYV52bb8qIkLwTQxbRGRH4jIa0VkxeAt7JElsOcONPLMvgY+dnUJM6ZEpnN1ssqLc2nr6eflulanQzmroqaZi4qmkJYc2c7iEhvJZOJEMPMgVvl/3jNkmwJXhz4c0zfg5Z7H9zI3N5MPrJ7ndDhBW+VfhnTLwaaItvmPpLffy8t1rbzn8rkRP3dZfjaP7ayno6c/oiv9GRNqwQxzvWqsfUzo/NcLNVSfbOc/37My4t98J2N6ViqLC6ew5VATH39DidPhsPf4GXr7vRGb/zDUYMmNgyfbWTY7J+LnNyZUxkwQIpIPfBUoVNVrRWQRcIWq/izs0UXAxso6Nmyuor6li8KcDNatKePm5UWOxHHvk/s50dpNWnIS7d3RXUI7kPJiDw9uOUJnbz+Zqc5+c95e458g58DVTFmBL0EcONFmCcLEtGD+ih8EfgF83v/4AL5KrDGfIDZW1rH+kV10+ddVrmvp4o7fv0xzZy9rFhdELI7Ne05w7x/309Pv6+Dt6ffyuUd3IyKOJKuJKi/28NPnDvPSkWZeV+ps6fWK2maKcjIoiNAEuaFmT8skPSXJ+iFMzAsmQXhU9WERWQ+gqv0iEhsr1Y9hw+aqs8lhUE+/l7sf28vdj+11KCqfrr4BNmyuiqkE8Zp500h1JbG1usnxBFFZ0+xI8xL4ZpeXzHBbTSYT84JJEB0ikouvYxoRuRyInqEqk1Df0jXic19/y5KIxXHH73cF3D5afNEoMzWZ5XNyHC//fby1i/rWbv6Pg53lpflu/nYw/lc4NPEtmATxGeAPwEIR2QLkAW8La1QRUpiTQV2AD+GinAze8Zo5EYvjvmerA8ZRmJMRsRhCZXWxh289fYDTHb1Mz0p1JIbBCXKRqOA6krKCbH5fcYyWzl5yMp15H4yZrGDmQewBXodvuOu/AIuB/eEMKlLWrSkjY1hBtYwUF+vWlCVkHKFQXuIb7vqCg7OqK2qbSUtO4sKZkZsgN9zgSKYDVnLDxLBgEsQLqtqvqntUdbeq9uFbCCjm3by8iK/dsoSinAwE35XD125ZEvF2/2iJIxQuLpqKOy3Z0WamitpmLp411dFlWQdHMllHtYllIzYxiUgBUARkiMhyYLBa3BQgMwKxRcTNy4ui4oM4WuKYrGRXEpctyHVsfYjuvgF217XygfL5jpx/UMGUdNxpyVZyw8S00fog1gDvA2YB3+JcgmgDPhfesEwsW12cyzP7Gjh6upPZ0yP7XWJPfSt9A+rYCKZBIkJpgZXcMLFtxAShqr8Efikib1HV30cwJhPjVvv7IbZUN3HrpZHr7IdzHdTRUO6jNN/Nk7uPR+VCSsYEI5hG2lkiMkV8HhCRChF5c9gjMzFrYV42M9xpjvRDVNQ2M3t6BnnutIife7iy/GyaO/tobO9xOhRjJiSYBPEBVT0DvBmYAbwfuHesg0Tk5yJyUkR2D9l2l4jUicgO/23tCMdeIyJVIlItIncG+buYKCEirC728MKhU3i9GrHzqirba5qj4uoBfKvLARw4YSOZTGwKJkEMXhuvBX6hqjuHbBvNg8A1AbZ/R1WX+W+bXnUyERdwP3AtsAi4zV//ycSQ8mIPpzp62R/BTtq6li5OtvU4Ov9hqMGhrtYPYWJVMAliu4g8hS9BbBYRNzDmqjCq+hxwegIxXQpUq+phVe0FfoMtWhRzBpchjeRoporaFiA6+h8APNlp5Gal2kgmE7OCSRAfBO4EXqOqnUAqvmamifqYiLzsb4IK9JdcBBwd8viYf1tAIvJhEdkmItsaG620QbQomJrOwrwsthyKYIKoaSYjxcUF/qadaFCa7+bASUsQJjYFkyBWA9nAxSJyJb6Z1DkTPN+PgIXAMuA4vuGzwwVqvhqxIVtVf6qqK1V1ZV6eswXizPlWF3t48fBpevsjswxpRW0zS2dPJdnl3AS54coK3Bw40YZq5PpijAmVYP6S1g25/TvwGHDXRE6mqg2qOqCqXuA/8TUnDXcMmD3k8SygfiLnM85aVeyhq2+AytrmsJ+ru2+AvfVnoqZ5aVBpvpuO3oGAtbaMiXZjJghVvWHI7U3ARUDDRE4mIjOHPPwnYHeA3V4CSkRkvoikArfiKxZoYszlC3JJksj0Q7x8rJV+r0ZhgsgGsNLfJiZN5Fr8GL4kMSoReQhfzaYyETkmIh8EviEiu0TkZeAq4NP+fQtFZBP41psAPgZsBvYBD6vqngnEaRw2NSOFi2flsCUChfsGV5BbPicn7Ocaj5LBkUw21NXEoGCWHP0+5/oAkvD1H+wc6zhVvS3A5oCr0KlqPb5RUoOPNwGvGgJrYk95cS4//uth2rr7cKenhO08FbXNzPdkkZvt/AS5oaZmpDBzajoH7QrCxKBgriC2Adv9txeAO1T1XWGNysSN8mIPA17lxcMTGfEcHFWlsrY56q4eBpXmW00mE5vGvILw12QyZkJWzJlGekoSWw418cZF+WE5x9HTXTS190Zd/8Og0vxsXjh8igGv4kqymkwmdoxW7nsXgYeXCqCqenHYojJxIz3FxWvmTQ9rR/X2Wt/VSfQmCDe9/V5qTnWwIC/b6XCMCdpoVxDXRywKE9fKiz3c+8f9nDzTzYwp6SF//YqaFrJSXWcX6Yk2g3EdaGizBGFiyoh9EKpao6o1/n0ahjw+SXC1mIwBfBPmALaGaTRTRW0zy+bkRG3zTfGMbERs+VETe4LppP4t59deGvBvMyYoi2ZOISczJSzlvzt6+tl3PPomyA2VmZrMnOmZ1lFtYk4wCSLZXzQPAP/91PCFZOJNUpKwaqFvGdJQl5zYeawFr+L4CnJjKZnhtqJ9JuYEkyAaReTGwQcichPg3Ir0JiaVF3s43trNK00dIX3dysEKrrOjO0GUFWTzSlNHxOpSGRMKwSSI24HPiUitiNQCdwAfDm9YJt6sDlP574qaZhbmZTE1M3yT8EKhNN9Nv1dDniCNCadgajEdUtXL8S3es1hVV6nqofCHZuLJnOmZFOVkhLQfQlWpqI2eFeRGMziSyfohTCwJuhaTqrarqv3vNhMydBnSgRAtQ/pKUwfNnX1Rs4LcaBZ4sklOEuuHMDElegrnm7hXXuLhTHc/u+taQ/J6Z1eQi4EEkZqcxDxPll1BmJgyZoIQkVdVPwu0zZixrFqYCxCyZqaK2mbc6ckUx8jks7J8t5X9NjElmCuIF4LcZsyoPNlpXFDgZmuIliGtqGlm2ewckqJ0gtxwpfluak930tU74HQoxgRlxAQhIgUicgmQISLLRWSF//Z6IDNSAZr4srrYw0tHmunum9yHZFt3H1UNbTHR/zCorCAbVag+aTOqTWwYrRbTGuB9+Jb8/Bbnymu0AZ8Lb1gmXpUXe3jg+VfYdqSZ1SWeCb/OzqOtqEZvgb5ASvPPjWRaMmuqw9HEro2VdWzYXEV9SxeFORmsW1PGzcuLnA4rLo2YIPxlvn8pIm9R1d9HMCYTxy6dP53kJOH56qZJJYjtNc2IwLIoXQMikLm5WaQmJ1k/xCRsrKxj/SMv09Xnm3BY19LF+kd2AViSCINg+iBmicgU8XlARCpE5M1hj8zEpay0ZFbMmTbpfoiK2mZKZ7iZEsZV6kLNlSQU52VTZUNdJ+zeP+4/mxwGdfUNsGFzlUMRxbdgEsQHVPUM8GZgBvB+4N6wRmXi2qriXHbVtdLS2Tv2zgF4vb4V5FbMzQltYBFQVuC25Ucn6MndJzhxpjvgc/UtXRGOJjEEkyAG+x7WAr9Q1Z1YuW8zCauLPajCCxMs/324qZ0z3f0sj6H+h0Gl+W7qW7s5093ndCgxo7Wrj888vIPbf7WdFFfgj57CnIwIR5UYgkkQ20XkKXwJYrOIuDm//HdAIvJzETkpIruHbNsgIvtF5GUReVREckY49oiI7BKRHSKyLcjfxcSIpbNzyEp1sWWCzUzba5qB2OqgHlRW4JuzYVcRwdlS3cS1332O/91RzyeuLubeW5aQkeI6bx9XkrBuTZlDEca3MdekBj4ILAMOq2qniOTia2Yay4PAD4D/GrLtaWC9qvaLyNeB9fiK/wVylapa1dg4lOJK4rIFuWypntgVREVNCzmZKSzwZIU4svArmeEfyXSinUvmTnc4mujV1TvA15/cz4Nbj7DAk8Xvbr/i7BWjKynp7Cim7LRk2nr6Y2YuTKwJJkEovkJ91wP3AFnAmOtGqupzIjJv2Lanhjz8O/DWoCM1caW82MOf9p+krqWLonE2D1TUNrM8hibIDVWUk0FWqstGMo1i59EWPv3wDg43dvC+VfO445oLyEg9d9Vw8/KisyOW+ga8vOMnL/D5R3axfHYOs6fbFK1QCqaJ6YfAFcBt/sdtwP0hOPcHgD+O8JwCT4nIdhEZtbS4iHxYRLaJyLbGxsYQhGUiYaLlv1u7+jh4sj0mm5fAt3hSiZXcCKhvwMu3nz7ALT/aSlfvAL/64GXcdePi85LDcCmuJL5363IAPvGbSvoGbL2NUAomQVymqh8FugFUtZlJrignIp8H+oFfj7BLuaquAK4FPioiV470Wqr6U1Vdqaor8/LyJhOWiaDS/Gw82WnjThCVtf7+hxiaQT2c1WR6tYMNbdzyw63c9+xBblpayJOfujLoeTKzp2fy1VuWUFnbwnefORDmSBNLMAmiT0Rc+L7VIyJ5BNFJPRIReS++5qp36gjrT6pqvf/nSeBR4NKJns9EJxGhvNjXDzGeZUgraltIEl9Hd6wqLXDT1N5LU3uP06E4zutVHvjbYa77/vMca+7kR+9cwbffsYypGeOb33LD0kLesXI2P/zLIbaGYe3zRBVMgrgP34f0DBH5CvA88NWJnExErsHXKX2jqnaOsE+Wf6QUIpKFb/7F7kD7mthWXuyhqb2HAw3B1yaqrG2mrGAK2WnBdJ9Fp9J830imRL+KONbcyT8/8He+/MQ+rizxsPnTV3LtkpkTfr0v3biI+Z4sPvU/OzhlyTckgllR7tfAvwFfA44DN6vqb8c6TkQewlf1tUxEjonIB/GNanIDT/uHsP7Yv2+hiGzyH5oPPC8iO4F/AE+o6pMT+N1MlCv390MEW/57wKtU1rawIobKawRS5q/JlKiLB6kqv912lGu++zd2HWvlG2+5mP98z0pmuMcc+zKqzNRkvn/bclo6+1j3u5fHdWVqAhvza5iITAdOAg8N2ZaiqqPO9FHV2wJs/tkI+9bjm2eBqh4Glo4Vl4l9RTkZzPdksaW6iQ+unj/m/gdPttHe0x9TFVwDyXOnkZOZwoEErOra1N7D+kd28fTeBi6dP51vvW1pSEceLS6cyufWXsBdj+3lwa1HeH/52P+vzMiCuU6vAGYDzfhmUOcAx0XkJPAhVd0evvBMvCsvzuXRijr6BrykuEa/oK2oaQFic4LcUCJCab474a4gntx9gs8/uou27n4+v/ZCPrh6fliGKr931Tz+drCJr23az6Xzp7O40CrnTlQwfRBPAmtV1aOqufhGFj0MfATfEFhjJmx1sYeO3gF2Hm0Zc9+K2mamZ6UyNzf2x7qX5bupamhLiGaQM919fPbhndz+q+0UTE3n8U+s5kNXLgjbPBYRYcPblpKTmcLHH6qks7c/LOdJBMEkiJWqunnwgX+y25Wq+nfAlh41k3L5glxEguuHqKhpZsWcHERib4LccKX52bR1949YfC5ebK1u4prvPMfGHXV84upiHv1I+dl1McJpelYq333HMl5p6uDuP+wN+/niVTAJ4rSI3CEic/23fwOa/UNfbVaKmZSczFSWFE1l6xhlN5o7ejnc1BHT8x+GOrt4UJw2M3X3DXD3Y3v45wdeJC3Fxe9uv4LPvLmM1ORgPnJCY1Wxh4+8fiH/s+0oj+2sj9h540kw/1r/jG9VuY3+22z/Nhfw9nAFZhJHebGHitpmOnpGbgqoPBq7BfoCGUwQB8cxxDdW7DzawnX3/Y1fbDnCe6+Yy6ZPvNaxyrufemMpy+fk8LlHdnH0dMCR9WYUoyYI/1XCd1X146q63H/7uKo2qmqvqlZHKE4Tx8oXeuj3Kv945fSI+2yvacaVJFwcJ0t1TstKZYY7jao4mgvRN+DlO/5SGZ3+Uhl333TRqKUywi3FlcR9VopjwkZNEKo6AOSJyKRKaxgzmpXzppGanDRqP0RFTQsXznSTmRq7E+SGKyuIn5Ib1Sd9pTK+9+xBbhxnqYxws1IcExfMX9sRYIuI/AHoGNyoqt8OV1AmsaSnuHjNvGkj1mXqH/Cy81gLb7tkVoQjC6+SGW7++x81eL0ak5VpwVcq48GtR/j6k/vJTHXxw3euYO0kZkOHyw1LC3n+YBM//Mshyhd6WFUcHckr2gXTB1EPPO7f1z3kZkzIrFroYf+JNhrbXl0ioaqhjc7egbjpoB5UVpBNd5+Xo82x2TZe19LFOx94kXse38vqYl+pjGhMDoOGluI43TGx5W4TzZhXEKp6dyQCMYltdbGHDZur2HqoiZuWFZ33XEUMryA3msGO6gMN7czNjf7FjzZW1p1dqGdqZgpdPf0ku5L4+luW8PaVs6N++PFgKY5/un8r6367kwfeuzLqY3bamFcQIpLnXyp0k4j8afAWieBM4rioaCpT0pMDDnetqG0hz53GrGnxte5wydkEEf39EBsr61j/yC7qWrpQoKWzjz6v8pk3l/KO18yJmQ/axYVTWb/2Ap7df5Jfbj3idDhRL5gmpl8D+4H5wN34+iReCmNMJgG5koRVCz08X930qtnFFbXxM0FuqOy0ZGZNy4iJuRAbNlfR1Tdw3javws+fP+JMQJPwvlXzuPqCGXx103721p9xOpyoFkyCyFXVnwF9qvpXVf0AcHmY4zIJqLw4l7qWLmpOnWuTb2rvoeZUZ9w1Lw0qjZHFg+pbusa1PZqJCBveerG/FEeFleIYRVALBvl/HheR60RkOb6Jc8aE1GD57y2Hzo1mOtv/EGcd1INK890camyP+vH5ee7AVXUKx7meeLTIzU7ju+9YxmErxTGqYBLEl0VkKvBZ4F+BB4BPhTMok5jme7IonJp+3nDXitoWUlzCkqL4mCA3XFlBNn0DSs2pjrF3doiqBlygKSPFxbo1ZQ5EFBpWimNswSwY9LiqtqrqblW9SlUvARZGIDaTYESEVcUeth46hdfr64eoqG1mUeFU0lOcm40bTudqMkVvyY0ndh3ncFMHb71kFkU5GQi+tTy+dssSbl5eNObx0cxKcYxuopWzPhPSKIzxW13soaWzj73Hz9A34OXlY7G/gtxoFuZlkyREbcmNrt4BvvrEPi6cOYWvv+Vittx5Na/cex1b7rw65pMDWCmOsUw0QcTXcBITNVYV5wK+8t/7jp+hu88btx3U4JtFPi83K2oXD/rxXw9R39rNXTcswhWjs73HMrQUx/eeOeh0OFFlooVt4n+VE+OIGe50SvOz2VLdRLq/NHSsLzE6lmgdyVTX0sWP/3qI6y6eyWULcp0OJ6xuWFrI3w42cv9fqllVnMuqhVaKA0a5ghCRNhE5E+DWBhRGMEaTYMqLPbx05DQvHD5FwZT0mB0pE6zSAjdHTnXQPWyegdO+umkfIvC5tRc6HUpE3HXjYuZ7svi0leI4a8QEoapuVZ0S4OZW1TGvPETk5yJyUkR2D9k2XUSeFpGD/p8BvxqKyDUiUiUi1SJy58R+NROrXElCd5+XzXsaaOnqZWNlndMhhVVZvhuvwqHG6Omo/vvhUzzx8nFuf91CiuI8QQ8aLMXR3NHHv/1uZ0IsBzuWcC7v9CBwzbBtdwLPqmoJ8Kz/8Xn8a1Dcj2/t60XAbSKyKIxxmiiysbKOX/295uzj7j4v6x/ZFddJojQ/G4iekhsDXuXux/ZSlJPBv1yZWAMWB0txPLPPSnFAGBOEqj4HDF8B5ibgl/77vwRuDnDopUC1qh5W1V7gN/7jTALYsLmK7r7zR5J09Q2wYXOVQxGF3zxPFikuiZqhrg/9o5Z9x8/wubUXOrrYj1OsFMc5kVsg1idfVY8D+H/OCLBPEXB0yONj/m0BiciHRWSbiGxrbGwMabAm8uKppEOwUlxJLMzL5mAUXEG0dvbxraequGz+dNYuKXA6HEdYKY5zIp0gghFoLN2IjYGq+lNVXamqK/Py8sIYlomEkTqk476jOt8dFXMhvvPMAVq7+rjrxsVxVxxxPHKz0/iOvxTHPY8lbimOSCeIBhGZCeD/eTLAPseA2UMez8K3aJFJAOvWlJExbNZ0rJd0CEZZgZtjzV209zj3bbXqRBv/7+81/PNlc7hw5hTH4ogW5cUe/u/rFvKbl47y+MuJ+REU6QTxB+C9/vvvBf43wD4vASUiMt+/Fvat/uNMArh5eRFfu2VJ3JV0GEvJDF9HtVPNTKrKPY/vITstmc++Kb6T8Xh8+k2lLJudw/oELcURthXgReQh4PWAR0SOAV8C7gUeFpEPArXA2/z7FgIPqOpaVe0XkY8BmwEX8HNV3ROuOE30uXl5UdwnhOHKCs4tHrTcgZnjm/c0sKX6FHfdsIhpWakRP3+0SnEl8f3blrP2e3/jk7+p5OF/uYJkVzS2zIdH2BKEqt42wlNvCLBvPbB2yONNwKYwhWZM1Jk9LZP0lCQONER+JFN33wBf2bSX0vxs3nX53IifP9oNluL4+EOVfPeZg/xrnDd3DpU4qdCYKJaUJI6V3PjZ869w9HQXX7phcUJ9Ox6PG5YW8vaVs7j/L9VsHbJeSbyz/w3GRInSfHfElx890drN/X+uZs3i/LMLNpnAErEUhyUIY6JEaX42J9t6aI7gh8+9f9xHv1f5wnVWrGAsmanJ3HdrYpXiCFsfhDFmfAYXDzrQ0BaR6qnba06zcUc9H7uqmNnTM8N+vnhwUdFU7rz2Au55fC/L7nmaM119FOZksG5NmSMDKzZW1rFhcxX1LV1hicOuIIyJEmdHMp0Mf0e116vc9Ye9FExJ5yNXJVa9pcmalplCkkBrVx+Kryy6E/XCNlbWsf6RXdS1dIUtDruCMCZKFExJx52eHJHFg367/Si76lr53q3LyEy1j4Hx+OZTB/AOa13q6hvgjt+/zCMRTBIvHj5FT3/gumWhuoqw/xnGRAkRoSwCJTfOdPexYXMVK+dO48altrTLeI1UF6yn38uZrr6IxTE8OQwKZd0ySxDGRJGSfDd/3H0cVQ1bLaT7njnIqY5eHnz/pQldb2miCnMyqAvwIVyUk8HGj5ZHLI7ye/8UMI5Q1i2zPghjokhZfjYtnX00tvWE5fWrT7bz4NYjvGPlbC4qmhqWc8S7aKkXFok47ArCmChSerbkRjszpqSH9LVVlf94fC8ZKa6Emg0caoPt++EcPRQtcViCMCaKlPmHulY1tLG6JLQT1/60/yR/PdDIF667EE92WkhfO9FES72wcMdhTUzGRJHc7DQ82akhH8nU0z/Afzy+lwV5WbzninkhfW0TvyxBGBNlSmaEfiTTg1uOcORUJ1+8fhGpyfZnb4Jj/1OMiTJlBW4ONrSFrJTDybZuvv+nat5wwQxeXxZolV9jArMEYUyUKc1309E7EHAI40R848kqevoH+ML1Vm/JjI8lCGOiTFmBb3W5UJT+3nG0hd9tP8YHVs9nvidr0q9nEoslCGOiTMngSKYTk6vJ5Ku3tIc8dxofv7okFKGZBGMJwpgoMyU9hZlT0yd9BfFoZR07jrZwxzUXkJ1mI9rN+FmCMCYKTXbxoPaefu59cj9LZ+dwSxSM1zexyRKEMVGorMBNdWM7A8PLhgbpB3+qprGth7tuWERSktVbMhMT8QQhImUismPI7YyIfGrYPq8XkdYh+3wx0nEa46TSfDe9/V5qTnWM+9gjTR38/PlXeMuKWSyfMy0M0ZlEEfGGSVWtApYBiIgLqAMeDbDr31T1+giGZkzUKBuyutyCvOxxHfvlJ/aS4hLuuMbqLZnJcbqJ6Q3AIVWtcTgOY6JK8YxsRMY/kumvBxp5Zt9JPnZ1SciL/ZnE43SCuBV4aITnrhCRnSLyRxFZPNILiMiHRWSbiGxrbGwMT5TGRFhGqos50zPHNZKpb8DLPY/tYV5uJh9YPS98wZmE4ViCEJFU4EbgtwGergDmqupS4PvAxpFeR1V/qqorVXVlXl5eWGI1xgml+e5xJYj/eqGGQ40dfOG6RaQlu8Y+wJgxOHkFcS1QoaoNw59Q1TOq2u6/vwlIEZHQ1j42JsqV5bt5pamDnv6BMfc91d7Dd585wJWlebzhQqu3ZELDyQRxGyM0L4lIgfjXQhSRS/HFeSqCsRnjuJL8bPq9yitNY49k+uZTVXT1DvDF6xfZMqImZBxJECKSCbwJeGTItttF5Hb/w7cCu0VkJ3AfcKuGqrSlMTGirGCw5MbozUy761r5zUtHee+qeRTPGN+IJ2NG48j8e1XtBHKHbfvxkPs/AH4Q6biMiSYLPNkkJ8mo/RCqvnpL0zNT+cQbrN6SCS2nRzEZY0aQmpzEfE8WBxpGHur6h531bKtpZt2aMqZmpEQwOpMILEEYE8VKC0YeydTZ28/XNu3noqIpvG3l7AhHZhKBJQhjoljpDDe1pzvp7O1/1XM/+sshTpzp5q4bFuOyeksmDCxBGBPFygqyUYXqk+c3Mx093clPnjvMTcsKWTlvukPRmXhnCcKYKFaaH3gk01ee2IdLhDuvvcCJsEyCsARhTBSbm5tFanISB4dcQWytbuLJPSf4yOsXMnNqhoPRmXhnCcKYKOZKEkpmZJ+9gugf8HL3Y3uZNS2DD125wOHoTLyzBGFMlBtak+m//1FLVUMbX7juQtJTrN6SCS9LEMZEudJ8N8dbu6k91cm3njrAqoW5rFlc4HRYJgHYSubGRLnTHT0AXLnhzwC8tsRj9ZZMRNgVhDFRbGNlHf/1wvnrad33bDUbK+scisgkEksQxkSxDZur6On3nretq2+ADZurHIrIJBJLEMZEsfqWrnFtNyaULEEYE8UKcwLPcxhpuzGhZAnCmCi2bk0ZGcOGs2akuFi3psyhiEwisVFMxkSxm5cXAb6+iPqWLgpzMli3puzsdmPCyRKEMVHu5uVFlhCMI6yJyRhjTECWIIwxxgRkCcIYY0xAliCMMcYEZAnCGGNMQKKqTscQMiLSCNSMuWNgHqAphOHEMnsvzmfvx/ns/TgnHt6LuaqaF+iJuEoQkyEi21R1pdNxRAN7L85n78f57P04J97fC2tiMsYYE5AlCGOMMQFZgjjnp04HEEXsvTifvR/ns/fjnLh+L6wPwhhjTEB2BWGMMSYgSxDGGGMCSvgEISLXiEiViFSLyJ1Ox+MkEZktIn8WkX0iskdEPul0TE4TEZeIVIrI407H4jQRyRGR34nIfv//kSucjslJIvJp/9/JbhF5SETSnY4p1BI6QYiIC7gfuBZYBNwmIoucjcpR/cBnVfVC4HLgown+fgB8EtjndBBR4nvAk6p6AbCUBH5fRKQI+ASwUlUvAlzArc5GFXoJnSCAS4FqVT2sqr3Ab4CbHI7JMap6XFUr/Pfb8H0AJOxCBCIyC7gOeMDpWJwmIlOAK4GfAahqr6q2OBqU85KBDBFJBjKBeofjCblETxBFwNEhj4+RwB+IQ4nIPGA58KLDoTjpu8C/AV6H44gGC4BG4Bf+JrcHRCTL6aCcoqp1wDeBWuA40KqqTzkbVegleoKQANsSftyviGQDvwc+papnnI7HCSJyPXBSVbc7HUuUSAZWAD9S1eVAB5CwfXYiMg1fa8N8oBDIEpF3ORtV6CV6gjgGzB7yeBZxeJk4HiKSgi85/FpVH3E6HgeVAzeKyBF8TY9Xi8ivnA3JUceAY6o6eEX5O3wJI1G9EXhFVRtVtQ94BFjlcEwhl+gJ4iWgRETmi0gqvk6mPzgck2NERPC1Me9T1W87HY+TVHW9qs5S1Xn4/l/8SVXj7htisFT1BHBURMr8m94A7HUwJKfVApeLSKb/7+YNxGGnfbLTAThJVftF5GPAZnyjEH6uqnscDstJ5cC7gV0issO/7XOqusm5kEwU+Tjwa/+XqcPA+x2OxzGq+qKI/A6owDf6r5I4LLthpTaMMcYElOhNTMYYY0ZgCcIYY0xAliCMMcYEZAnCGGNMQJYgjDHGBGQJwiQ8EckVkR3+2wkRqfPfbxeRH0YohtcPVowVkRsTvbKwiQ4JPQ/CGABVPQUsAxCRu4B2Vf2mg/H8gQSesGmih11BGDOCYd/q7xKRX4rIUyJyRERuEZFviMguEXnSX6IEEblERP4qIttFZLOIzAzwum/zryGwU0SeC/D8+0TkB/77+SLyqH/fnSKyyr/9XSLyD/+Vzk/8peuNCSlLEMYEbyG+8t83Ab8C/qyqS4Au4Dp/kvg+8FZVvQT4OfCVAK/zRWCNqi4FbhzjnPcBf/XvuwLYIyIXAu8AylV1GTAAvHOyv5wxw1kTkzHB+6Oq9onILnylWZ70b98FzAPKgIuAp33leXDhKwU93BbgQRF5GF+Rt9FcDbwHQFUHgFYReTdwCfCS/zwZwMmJ/1rGBGYJwpjg9QCoqldE+vRcnRovvr8lAfao6qhLcarq7SJyGb6rkR0ismyccQjwS1VdP87jjBkXa2IyJnSqgLzBtZpFJEVEFg/fSUQWquqLqvpFoInzS84P9yzwf/3Hufwruz0LvFVEZvi3TxeRuSH+XYyxBGFMqPiXrX0r8HUR2QnsIPAaARv8ndu7geeAnaO87CeBq/zNWtuBxaq6F/gC8JSIvAw8DbyqM9yYybJqrsYYYwKyKwhjjDEBWYIwxhgTkCUIY4wxAVmCMMYYE5AlCGOMMQFZgjDGGBOQJQhjjDEB/X/9lbHHhZJwEwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(len(q6clusters)), [max(x) for x in q6clusters], 'o-')\n",
    "plt.xlabel(\"Time slice\")\n",
    "plt.ylabel(\"Largest cluster size\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final thing to do is to remove the split files after use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "for file in files:\n",
    "    os.remove(file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using [ASE](https://wiki.fysik.dtu.dk/ase/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above example can also done using ASE. The ASE read method needs to be imported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.io import read"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = read(\"traj.light\", format=\"lammps-dump-text\", index=\":\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above function, `index=\":\"` tells ase to read the complete trajectory. The individual slices can now be accessed by indexing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Atoms(symbols='H500', pbc=True, cell=[18.21922, 18.22509, 18.36899], momenta=...)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "traj[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the same functions as above, but by specifying a different file format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "q6clusters_ase = [calculate_q6_cluster(x, format=\"ase\") for x in traj]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will plot and compare with the results from before,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Largest cluster size')"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4R0lEQVR4nO3dd3hU95Xw8e/RqEsDAo2QkOio2GBMMXFBxImdgo3rOs3e9ORN1m96WdYmySa2N8UJqU6ctk7i7Juss05is7FNjEuKY3AcgwSmCgRGAgmEBJJQb3PeP2YEQh5JI2lm7pTzeZ55NHPn3rlHA5oz91fOT1QVY4wxZrgkpwMwxhgTnSxBGGOMCcgShDHGmIAsQRhjjAnIEoQxxpiAkp0OIJQ8Ho/OmzfP6TCMMSZmbN++vUlV8wI9F1cJYt68eWzbts3pMIwxJmaISM1Iz1kTkzHGmIAsQRhjjAnIEoQxxpiALEEYY4wJyBKEMcaYgOJqFJMxobSxso4Nm6uob+miMCeDdWvKuHl5kdNhGRMxliCMCWBjZR3rH9lFV98AAHUtXax/ZBeAJQmTMKyJyZgANmyuOpscBnX1DbBhc5VDERkTeZYgjAmgvqVrXNuNiUeWIIwJoDAnY1zbjYlHliCMCWDdmjLSks//88hIcbFuTZlDERkTeWFLECIyW0T+LCL7RGSPiHzSv32DiOwXkZdF5FERyRnh+CMisktEdoiIFVgyEXXz8iLWXlRw9nFORgpfu2WJdVCbhBLOK4h+4LOqeiFwOfBREVkEPA1cpKoXAweA9aO8xlWqukxVV4YxTmMC6vMqBVPSSU9J4i2XzLLkYBJO2BKEqh5X1Qr//TZgH1Ckqk+par9/t78Ds8IVgzGTUVnbwiXzplEyw82BhjanwzEm4iLSByEi84DlwIvDnvoA8McRDlPgKRHZLiIfHuW1Pywi20RkW2NjY0jiNeZEazd1LV1cMmcapfmWIExiCnuCEJFs4PfAp1T1zJDtn8fXDPXrEQ4tV9UVwLX4mqeuDLSTqv5UVVeq6sq8vIBrXhgzbhW1zQCsmDuNsoJsGs700NLZ63BUxkRWWBOEiKTgSw6/VtVHhmx/L3A98E5V1UDHqmq9/+dJ4FHg0nDGasxQFTXNpCUnsWjmFEry3QAcaGh3OCpjIiuco5gE+BmwT1W/PWT7NcAdwI2q2jnCsVki4h68D7wZ2B2uWI0ZrqK2mSVFU0lNTqLMnyCqrJnJJJhwXkGUA+8GrvYPVd0hImuBHwBu4Gn/th8DiEihiGzyH5sPPC8iO4F/AE+o6pNhjNWYs3r6B9hdd4ZL5k4DYObUdNxpyRy0BGESTNiK9anq84AEeGpTgG2DTUpr/fcPA0vDFZsxo9ldd4beAS/L5/gShIhQWuCm6oQlCJNYbCa1McNUnu2gzjm7bXAk0whdZsbEJUsQxgyzvaaZ2dMzmOFOP7utND+b5s4+Gtt7HIzMmMiyBGHMEKpKRW0zK/zNS4MGO6oPnLCRTCZxWIIwZoj61m4azvS8KkGUFgwOdbV+CJM4LEEYM0RFjb//YViC8GSnkZuVagnCJBRLEMYMsb2mmYwUFxfMdL/quZL8bJsLYRKKJQhjhqisbebiWVNJcb36T6Ms383BhnYbyWQShiUIY/y6+wbYU3+GFXOnBXy+tMBNe08/9a3dEY7MGGeEbaKciT0bK+vYsLmK+pYuCnMyWLemLKHWQHj5WCv9Xn1V/8OgcyOZ2iiypUdNArArCAP4ksP6R3ZR19KFAnUtXax/ZBcbK+ucDi1izlZwnZMT8PkSq8lkEowlCAPAhs1VdPUNnLetq2+ADZurHIoo8ipqmpmXm0ludlrA56dmpFAwJZ0DVnLDJAhLEAaA+paucW2PN74Jci0jNi8NKi1w2xWESRiWIAwAhSO0qY+0Pd4cPd1FU3vPiB3Ug8rys6k+2c6A10YymfhnCcIAsG5NGa6k84vvZqS4WLemzKGIIutc/8MYVxD5bnr6vdSeDriUiTFxxRKEAWBhXjYDXiU7zXV22z03LU6YUUwVtc1kpbooK3j1BLmhBp+30t8mEViCMKgqdz+2B092KlvXv4FfvO81QOI0L4EvQSydnfOqq6jhimdkA1aTySQGSxCGP+ysZ1tNM+vWlDElPYVL508nOUl4vrrJ6dAiorO3n33H286uIDeazNRk5kzPtI5qkxAsQSS4zt5+vrZpPxcVTeGtl8wGICstmeVzctiSIAli59FWBkaZIDdcab7blh81CcESRIL74Z8PceJMN3fdsPi85pXyYg+76lpp7exzMLrIGOygXj7CBLnhygqyOdzYQW+/N4xRGeM8SxAJrPZUJz/922FuWlbIynnTz3uuvNiDKrxwOP6vIipqmlmQl0VOZmpQ+5fmu+n3Kq80dYQ5MmOcFbYEISKzReTPIrJPRPaIyCf926eLyNMictD/M+B1vYhcIyJVIlItIneGK85E9pVNe3GJcOe1F7zquWWzc8hKdcV9P4SqUnm0hUuCbF4CX4IAK7lh4l84ryD6gc+q6oXA5cBHRWQRcCfwrKqWAM/6H59HRFzA/cC1wCLgNv+xJkS2VDexeU8DH71qITOnvnq0UooricsW5LK1+pQD0UXOkVOdnO7oHXOC3FAL8rJwJYmV3DBxL2wJQlWPq2qF/34bsA8oAm4Cfunf7ZfAzQEOvxSoVtXDqtoL/MZ/nAmB/gEvdz+2h9nTM/g/r10w4n6rFuZyuKmDujgutzHSCnKjSUt2Md+TZUNdTdyLSB+EiMwDlgMvAvmqehx8SQSYEeCQIuDokMfH/NsCvfaHRWSbiGxrbGwMadzx6tcv1nKgoZ3Pr11EeoprxP1Wl3gA4no00/baZtxpyZT45zcEqyzfbQnCxL2wJwgRyQZ+D3xKVc8Ee1iAbQGL36jqT1V1paquzMvLm2iYCaO5o5dvP32A8uJc1izOH3Xfsnw3nuzUuE4QFTXNLJuTQ9IYE+SGK813U3O6k67egbF3NiZGhTVBiEgKvuTwa1V9xL+5QURm+p+fCZwMcOgxYPaQx7OA+nDGmii+9XQV7T39fOmGxYiM/qEoIpQXe9hSfSoul9ls7+nnQEPbuJqXBpXmZ6MK1SfbwxCZMdEhnKOYBPgZsE9Vvz3kqT8A7/Xffy/wvwEOfwkoEZH5IpIK3Oo/zkzC3voz/PeLtbz78rlnR+KMpXyhh6b2Hg40xN8H4c6jLXiVcXVQDyotsJFMJv6F8wqiHHg3cLWI7PDf1gL3Am8SkYPAm/yPEZFCEdkEoKr9wMeAzfg6tx9W1T1hjDXuDdZbmpqRwqffWBr0ceX+foh4HO66vaYZEd+Q3vGaOz2T1OQkm1Ft4lrY1qRW1ecJ3JcA8IYA+9cDa4c83gRsCk90iWfTrhO8+MppvnzzRUzNTAn6uKKcDOZ7stha3cQHV88PY4SRV1HbTMmMbKZmBP9+DEp2JVGcl21XECau2UzqBNDVO8BXN+3jggI3t106Z9zHr1qYy98Pn6JvIH5KS3i9SmUQK8iNpqzAbXMhTFyzBJEAfvLcIepaurjrxsVjlrMOZHWxh47eAXYebQl9cA453NROa1ffpBJESX429a3dnOmO/3pVJjFZgohzdS1d/Pivh7huyUwuX5A7ode4YmEuIrAljmZVV9S0ABProB5U5u/ot34IE68sQcS5r27ahyqsX/vqekvByslM5aLCqXE1H6KitpmpGSks8GRN+DUGR4LF4wgvY8ASRFx78fApnnj5OLe/biGzpmVO6rXKiz1U1DbT0dMfouicVVHbzPIJTJAbqigng6xUly0/auKWJYg4NeBV7npsL4VT07n9dQsn/Xqriz30e5V/HDkdguic1drVx4GG9nFVcA0kKUkosZIbJo5ZgohTv3mpln3Hz/C56y4kI3XkekvBWjlvGqnJSWw5GPvNTDv8ne2T6X8YVJqfbQnCxC1LEHGotbOPb26u4rL507luycyQvGZ6iouVc6fFxYS5ippmkgSWTmCC3HCl+W6a2ns51d4z+cCMiTKWIOLQd545QGtXX1D1lsajvNjD/hNtNMX4h2FFbTOl+W6y0yY/T7SswDqqTfyyBBFnDjS08f/+XsNtl85hUeGUkL52ebGv7MbWQ7E73NXrVXbUtnBJCJqX4NxQV2tmMvHIEkQcUVXueWwvWakuPvvmspC//pKiqbjTk2O6H+LgyXbaevonNUFuqDx3GjmZKVZyw8SloBKEiKwWkff77+eJSHwV5YkTT+1t4PnqJj7zplKmZ6WG/PVdScKqhbk8X90Us+W/K2r9K8iF6ApCRCidYSU3THwaM0GIyJeAO4D1/k0pwK/CGZQZv+6+Ab78xF5K87N51+Vzw3ae1cUe6lq6qD3dGbZzhNP2mmamZ6UyL3dy80KGKi3wFe2L1aRpzEiCuYL4J+BGoAPOVl0NbjEBEzE/e/4Vjp7u4ovXLybZFb6Ww1XFsV3+u6K2mRVzckLaeV+W76atu5+GM7HdeW/McMF8kvSq76uRAojIxGsTmLA40drN/X+u5s2L8s+uIx0uCzxZzJyaztYYrMvU0tnL4cYOloeo/2HQYMkN64cw8SaYBPGwiPwEyBGRDwHPAA+ENywzHl9/cj/9XuUL1y0K+7nOLkN6qAmvN7aaVCprWwBC1kE96GxNJuuHMHFmzAShqt8Efodvbeky4Iuqel+4AzPB2V7TzKOVdXzotfOZE8J29dGUF+fS0tnH3uNnInK+UNle04wrSVg6e2pIX3daVip57jS7gjBxJ5hO6n8H9qvqOlX9V1V9WkQ+HIHYzBi8Xt8yovlT0vjI64sjdt7yhbHZD1FR28yFM91kpoZ+IcUyq8lk4lAwTUwfBzaLyFVDtt0epnjMOPxu+zFePtbK+msvJCsEs4KDNWNKOqX52TFV/nvAq+w8OrkV5EZTmu/mYEN7zDW7GTOaYBJEHXANcK+IrPNvC90QEDMhZ7r7+Mbm/Vwydxo3LSuM+PlXLfTw0pHTdPcNRPzcE1F1oo2O3oGwJYiygmy6+gY41twVltc3xglBjYdU1VrgdcAiEfktkDHWMSLycxE5KSK7h2z7HxHZ4b8dEZEdIxx7RER2+ffbFtyvkli+/+xBTnX0cleI6y0Fa3Wxh+4+79mJZ9Fuuz/OUJXYGK7ERjKZOBRMgtgGoKrdqvp+4C9AMNN0H8R35XGWqr5DVZep6jJ8nd6PjHL8Vf59VwZxroRyqLGdX2w5wtsvmc2SWaHtcA3WZQum40qSmBnuWlnTjCc7jVnTxvxuMyElM7IBq8lk4kswo5g+NOzx/aq6IIjjngMCri4jvq+8bwceCjJOM8R/PL6XjBQX664Jfb2lYLnTU1g6a2rMdFSHY4LcUO70FIpyMixBmLgyYoIQkYf9P3eJyMvDb5M872uBBlU9OMLzCjwlItvHGjElIh8WkW0isq2xsXGSYUW/P+1v4C9VjXzyjSV4stMcjWV1sYeXj7XQ2tXnaBxjaWrv4cipzpDVXxpJWYHblh81cWW0K4hP+n9eD9wQ4DYZtzH61UO5qq4ArgU+KiJXjrSjqv5UVVeq6sq8vLxJhhXdevu9/Mfj+1iQl8V7rpjndDiUF3vwqm/t62g2OEEuXP0Pg0rz3Rxu7KBvwBvW8xgTKSMmCFU97r/bBBxV1RogDVgK1E/0hCKSDNwC/M8o5673/zwJPApcOtHzxZNfbHmFV5o6+PfrF5Ga7Hyl9uVzppGR4or64a4Vtc0kJwlLisLbX1Oan03vgJeaUx1hPY8xkRLMp8xzQLqIFAHPAu/H1wE9UW/EN/HuWKAnRSRLRNyD94E3A7sD7ZtITrZ18/0/VXP1BTO4qmyG0+EAkJqcxKXzp0d9P0RFTTOLC6eQnjL5tblHc7Ym0wlbXc7Eh2AShKhqJ75v/d9X1X8Cxiz6IyIPAS8AZSJyTEQ+6H/qVoY1L4lIoYhs8j/MB54XkZ3AP4AnVPXJ4H6d+LXhySp6+gf49+vDX29pPFYXezjU2MGJ1m6nQwmob8DLzmMtIS/QF0jxjGySxEYymfgRzPRbEZErgHcCgx/yYx6nqreNsP19AbbVA2v99w/ja8YyfjuPtvDb7cf4lysXMN8TXcV0VxXnArCluom3XDLL4Whebf/xNrr7vGHvfwBIT3ExLzfLEoSJG8FcQXwS32JBj6rqHhFZAPw5vGGZQV6vctdje/Bkp/GxqyNXbylYFxZMYXpWatT2Q4R6BbmxlOa7bbKciRvBXAk8h68fYvDxYeAT4QzKnLNxRx2VtS1seOvFuNNTnA7nVZL8y5BuOeRbhtSJWd2j2V7TTP6UNAqnpkfkfKX52Ty19wTdfQNh7/MwJtycHwpjRtTe08+9f9zP0tk5vGVF9DXfDFpd7KHhTA+HGqOvc7aitplL5k6LWOIqLXDjVaLyvTBmvCxBRLH7/1zNybYe7rphEUlJ0fXNfKjywWVID0ZXM9PJtm6ONXeFrUBfIGX+kUwHGyxBmNg3aoIQEZeIfDpSwZhzjjR18LO/vcItK4oiMgJnMmZPz2TO9Ey2HIquCXMVNS0AEX3/5nmySHGJ9UOYuDBqglDVAeCmCMVihvjyE/tIcQl3XnOB06EEpbzYw98PnaI/imYRV9Q2k+pK4qKiKRE7Z4oriYV52bb8qIkLwTQxbRGRH4jIa0VkxeAt7JElsOcONPLMvgY+dnUJM6ZEpnN1ssqLc2nr6eflulanQzmroqaZi4qmkJYc2c7iEhvJZOJEMPMgVvl/3jNkmwJXhz4c0zfg5Z7H9zI3N5MPrJ7ndDhBW+VfhnTLwaaItvmPpLffy8t1rbzn8rkRP3dZfjaP7ayno6c/oiv9GRNqwQxzvWqsfUzo/NcLNVSfbOc/37My4t98J2N6ViqLC6ew5VATH39DidPhsPf4GXr7vRGb/zDUYMmNgyfbWTY7J+LnNyZUxkwQIpIPfBUoVNVrRWQRcIWq/izs0UXAxso6Nmyuor6li8KcDNatKePm5UWOxHHvk/s50dpNWnIS7d3RXUI7kPJiDw9uOUJnbz+Zqc5+c95e458g58DVTFmBL0EcONFmCcLEtGD+ih8EfgF83v/4AL5KrDGfIDZW1rH+kV10+ddVrmvp4o7fv0xzZy9rFhdELI7Ne05w7x/309Pv6+Dt6ffyuUd3IyKOJKuJKi/28NPnDvPSkWZeV+ps6fWK2maKcjIoiNAEuaFmT8skPSXJ+iFMzAsmQXhU9WERWQ+gqv0iEhsr1Y9hw+aqs8lhUE+/l7sf28vdj+11KCqfrr4BNmyuiqkE8Zp500h1JbG1usnxBFFZ0+xI8xL4ZpeXzHBbTSYT84JJEB0ikouvYxoRuRyInqEqk1Df0jXic19/y5KIxXHH73cF3D5afNEoMzWZ5XNyHC//fby1i/rWbv6Pg53lpflu/nYw/lc4NPEtmATxGeAPwEIR2QLkAW8La1QRUpiTQV2AD+GinAze8Zo5EYvjvmerA8ZRmJMRsRhCZXWxh289fYDTHb1Mz0p1JIbBCXKRqOA6krKCbH5fcYyWzl5yMp15H4yZrGDmQewBXodvuOu/AIuB/eEMKlLWrSkjY1hBtYwUF+vWlCVkHKFQXuIb7vqCg7OqK2qbSUtO4sKZkZsgN9zgSKYDVnLDxLBgEsQLqtqvqntUdbeq9uFbCCjm3by8iK/dsoSinAwE35XD125ZEvF2/2iJIxQuLpqKOy3Z0WamitpmLp411dFlWQdHMllHtYllIzYxiUgBUARkiMhyYLBa3BQgMwKxRcTNy4ui4oM4WuKYrGRXEpctyHVsfYjuvgF217XygfL5jpx/UMGUdNxpyVZyw8S00fog1gDvA2YB3+JcgmgDPhfesEwsW12cyzP7Gjh6upPZ0yP7XWJPfSt9A+rYCKZBIkJpgZXcMLFtxAShqr8Efikib1HV30cwJhPjVvv7IbZUN3HrpZHr7IdzHdTRUO6jNN/Nk7uPR+VCSsYEI5hG2lkiMkV8HhCRChF5c9gjMzFrYV42M9xpjvRDVNQ2M3t6BnnutIife7iy/GyaO/tobO9xOhRjJiSYBPEBVT0DvBmYAbwfuHesg0Tk5yJyUkR2D9l2l4jUicgO/23tCMdeIyJVIlItIncG+buYKCEirC728MKhU3i9GrHzqirba5qj4uoBfKvLARw4YSOZTGwKJkEMXhuvBX6hqjuHbBvNg8A1AbZ/R1WX+W+bXnUyERdwP3AtsAi4zV//ycSQ8mIPpzp62R/BTtq6li5OtvU4Ov9hqMGhrtYPYWJVMAliu4g8hS9BbBYRNzDmqjCq+hxwegIxXQpUq+phVe0FfoMtWhRzBpchjeRoporaFiA6+h8APNlp5Gal2kgmE7OCSRAfBO4EXqOqnUAqvmamifqYiLzsb4IK9JdcBBwd8viYf1tAIvJhEdkmItsaG620QbQomJrOwrwsthyKYIKoaSYjxcUF/qadaFCa7+bASUsQJjYFkyBWA9nAxSJyJb6Z1DkTPN+PgIXAMuA4vuGzwwVqvhqxIVtVf6qqK1V1ZV6eswXizPlWF3t48fBpevsjswxpRW0zS2dPJdnl3AS54coK3Bw40YZq5PpijAmVYP6S1g25/TvwGHDXRE6mqg2qOqCqXuA/8TUnDXcMmD3k8SygfiLnM85aVeyhq2+AytrmsJ+ru2+AvfVnoqZ5aVBpvpuO3oGAtbaMiXZjJghVvWHI7U3ARUDDRE4mIjOHPPwnYHeA3V4CSkRkvoikArfiKxZoYszlC3JJksj0Q7x8rJV+r0ZhgsgGsNLfJiZN5Fr8GL4kMSoReQhfzaYyETkmIh8EviEiu0TkZeAq4NP+fQtFZBP41psAPgZsBvYBD6vqngnEaRw2NSOFi2flsCUChfsGV5BbPicn7Ocaj5LBkUw21NXEoGCWHP0+5/oAkvD1H+wc6zhVvS3A5oCr0KlqPb5RUoOPNwGvGgJrYk95cS4//uth2rr7cKenhO08FbXNzPdkkZvt/AS5oaZmpDBzajoH7QrCxKBgriC2Adv9txeAO1T1XWGNysSN8mIPA17lxcMTGfEcHFWlsrY56q4eBpXmW00mE5vGvILw12QyZkJWzJlGekoSWw418cZF+WE5x9HTXTS190Zd/8Og0vxsXjh8igGv4kqymkwmdoxW7nsXgYeXCqCqenHYojJxIz3FxWvmTQ9rR/X2Wt/VSfQmCDe9/V5qTnWwIC/b6XCMCdpoVxDXRywKE9fKiz3c+8f9nDzTzYwp6SF//YqaFrJSXWcX6Yk2g3EdaGizBGFiyoh9EKpao6o1/n0ahjw+SXC1mIwBfBPmALaGaTRTRW0zy+bkRG3zTfGMbERs+VETe4LppP4t59deGvBvMyYoi2ZOISczJSzlvzt6+tl3PPomyA2VmZrMnOmZ1lFtYk4wCSLZXzQPAP/91PCFZOJNUpKwaqFvGdJQl5zYeawFr+L4CnJjKZnhtqJ9JuYEkyAaReTGwQcichPg3Ir0JiaVF3s43trNK00dIX3dysEKrrOjO0GUFWTzSlNHxOpSGRMKwSSI24HPiUitiNQCdwAfDm9YJt6sDlP574qaZhbmZTE1M3yT8EKhNN9Nv1dDniCNCadgajEdUtXL8S3es1hVV6nqofCHZuLJnOmZFOVkhLQfQlWpqI2eFeRGMziSyfohTCwJuhaTqrarqv3vNhMydBnSgRAtQ/pKUwfNnX1Rs4LcaBZ4sklOEuuHMDElegrnm7hXXuLhTHc/u+taQ/J6Z1eQi4EEkZqcxDxPll1BmJgyZoIQkVdVPwu0zZixrFqYCxCyZqaK2mbc6ckUx8jks7J8t5X9NjElmCuIF4LcZsyoPNlpXFDgZmuIliGtqGlm2ewckqJ0gtxwpfluak930tU74HQoxgRlxAQhIgUicgmQISLLRWSF//Z6IDNSAZr4srrYw0tHmunum9yHZFt3H1UNbTHR/zCorCAbVag+aTOqTWwYrRbTGuB9+Jb8/Bbnymu0AZ8Lb1gmXpUXe3jg+VfYdqSZ1SWeCb/OzqOtqEZvgb5ASvPPjWRaMmuqw9HEro2VdWzYXEV9SxeFORmsW1PGzcuLnA4rLo2YIPxlvn8pIm9R1d9HMCYTxy6dP53kJOH56qZJJYjtNc2IwLIoXQMikLm5WaQmJ1k/xCRsrKxj/SMv09Xnm3BY19LF+kd2AViSCINg+iBmicgU8XlARCpE5M1hj8zEpay0ZFbMmTbpfoiK2mZKZ7iZEsZV6kLNlSQU52VTZUNdJ+zeP+4/mxwGdfUNsGFzlUMRxbdgEsQHVPUM8GZgBvB+4N6wRmXi2qriXHbVtdLS2Tv2zgF4vb4V5FbMzQltYBFQVuC25Ucn6MndJzhxpjvgc/UtXRGOJjEEkyAG+x7WAr9Q1Z1YuW8zCauLPajCCxMs/324qZ0z3f0sj6H+h0Gl+W7qW7s5093ndCgxo7Wrj888vIPbf7WdFFfgj57CnIwIR5UYgkkQ20XkKXwJYrOIuDm//HdAIvJzETkpIruHbNsgIvtF5GUReVREckY49oiI7BKRHSKyLcjfxcSIpbNzyEp1sWWCzUzba5qB2OqgHlRW4JuzYVcRwdlS3cS1332O/91RzyeuLubeW5aQkeI6bx9XkrBuTZlDEca3MdekBj4ILAMOq2qniOTia2Yay4PAD4D/GrLtaWC9qvaLyNeB9fiK/wVylapa1dg4lOJK4rIFuWypntgVREVNCzmZKSzwZIU4svArmeEfyXSinUvmTnc4mujV1TvA15/cz4Nbj7DAk8Xvbr/i7BWjKynp7Cim7LRk2nr6Y2YuTKwJJkEovkJ91wP3AFnAmOtGqupzIjJv2Lanhjz8O/DWoCM1caW82MOf9p+krqWLonE2D1TUNrM8hibIDVWUk0FWqstGMo1i59EWPv3wDg43dvC+VfO445oLyEg9d9Vw8/KisyOW+ga8vOMnL/D5R3axfHYOs6fbFK1QCqaJ6YfAFcBt/sdtwP0hOPcHgD+O8JwCT4nIdhEZtbS4iHxYRLaJyLbGxsYQhGUiYaLlv1u7+jh4sj0mm5fAt3hSiZXcCKhvwMu3nz7ALT/aSlfvAL/64GXcdePi85LDcCmuJL5363IAPvGbSvoGbL2NUAomQVymqh8FugFUtZlJrignIp8H+oFfj7BLuaquAK4FPioiV470Wqr6U1Vdqaor8/LyJhOWiaDS/Gw82WnjThCVtf7+hxiaQT2c1WR6tYMNbdzyw63c9+xBblpayJOfujLoeTKzp2fy1VuWUFnbwnefORDmSBNLMAmiT0Rc+L7VIyJ5BNFJPRIReS++5qp36gjrT6pqvf/nSeBR4NKJns9EJxGhvNjXDzGeZUgraltIEl9Hd6wqLXDT1N5LU3uP06E4zutVHvjbYa77/vMca+7kR+9cwbffsYypGeOb33LD0kLesXI2P/zLIbaGYe3zRBVMgrgP34f0DBH5CvA88NWJnExErsHXKX2jqnaOsE+Wf6QUIpKFb/7F7kD7mthWXuyhqb2HAw3B1yaqrG2mrGAK2WnBdJ9Fp9J830imRL+KONbcyT8/8He+/MQ+rizxsPnTV3LtkpkTfr0v3biI+Z4sPvU/OzhlyTckgllR7tfAvwFfA44DN6vqb8c6TkQewlf1tUxEjonIB/GNanIDT/uHsP7Yv2+hiGzyH5oPPC8iO4F/AE+o6pMT+N1MlCv390MEW/57wKtU1rawIobKawRS5q/JlKiLB6kqv912lGu++zd2HWvlG2+5mP98z0pmuMcc+zKqzNRkvn/bclo6+1j3u5fHdWVqAhvza5iITAdOAg8N2ZaiqqPO9FHV2wJs/tkI+9bjm2eBqh4Glo4Vl4l9RTkZzPdksaW6iQ+unj/m/gdPttHe0x9TFVwDyXOnkZOZwoEErOra1N7D+kd28fTeBi6dP51vvW1pSEceLS6cyufWXsBdj+3lwa1HeH/52P+vzMiCuU6vAGYDzfhmUOcAx0XkJPAhVd0evvBMvCsvzuXRijr6BrykuEa/oK2oaQFic4LcUCJCab474a4gntx9gs8/uou27n4+v/ZCPrh6fliGKr931Tz+drCJr23az6Xzp7O40CrnTlQwfRBPAmtV1aOqufhGFj0MfATfEFhjJmx1sYeO3gF2Hm0Zc9+K2mamZ6UyNzf2x7qX5bupamhLiGaQM919fPbhndz+q+0UTE3n8U+s5kNXLgjbPBYRYcPblpKTmcLHH6qks7c/LOdJBMEkiJWqunnwgX+y25Wq+nfAlh41k3L5glxEguuHqKhpZsWcHERib4LccKX52bR1949YfC5ebK1u4prvPMfGHXV84upiHv1I+dl1McJpelYq333HMl5p6uDuP+wN+/niVTAJ4rSI3CEic/23fwOa/UNfbVaKmZSczFSWFE1l6xhlN5o7ejnc1BHT8x+GOrt4UJw2M3X3DXD3Y3v45wdeJC3Fxe9uv4LPvLmM1ORgPnJCY1Wxh4+8fiH/s+0oj+2sj9h540kw/1r/jG9VuY3+22z/Nhfw9nAFZhJHebGHitpmOnpGbgqoPBq7BfoCGUwQB8cxxDdW7DzawnX3/Y1fbDnCe6+Yy6ZPvNaxyrufemMpy+fk8LlHdnH0dMCR9WYUoyYI/1XCd1X146q63H/7uKo2qmqvqlZHKE4Tx8oXeuj3Kv945fSI+2yvacaVJFwcJ0t1TstKZYY7jao4mgvRN+DlO/5SGZ3+Uhl333TRqKUywi3FlcR9VopjwkZNEKo6AOSJyKRKaxgzmpXzppGanDRqP0RFTQsXznSTmRq7E+SGKyuIn5Ib1Sd9pTK+9+xBbhxnqYxws1IcExfMX9sRYIuI/AHoGNyoqt8OV1AmsaSnuHjNvGkj1mXqH/Cy81gLb7tkVoQjC6+SGW7++x81eL0ak5VpwVcq48GtR/j6k/vJTHXxw3euYO0kZkOHyw1LC3n+YBM//Mshyhd6WFUcHckr2gXTB1EPPO7f1z3kZkzIrFroYf+JNhrbXl0ioaqhjc7egbjpoB5UVpBNd5+Xo82x2TZe19LFOx94kXse38vqYl+pjGhMDoOGluI43TGx5W4TzZhXEKp6dyQCMYltdbGHDZur2HqoiZuWFZ33XEUMryA3msGO6gMN7czNjf7FjzZW1p1dqGdqZgpdPf0ku5L4+luW8PaVs6N++PFgKY5/un8r6367kwfeuzLqY3bamFcQIpLnXyp0k4j8afAWieBM4rioaCpT0pMDDnetqG0hz53GrGnxte5wydkEEf39EBsr61j/yC7qWrpQoKWzjz6v8pk3l/KO18yJmQ/axYVTWb/2Ap7df5Jfbj3idDhRL5gmpl8D+4H5wN34+iReCmNMJgG5koRVCz08X930qtnFFbXxM0FuqOy0ZGZNy4iJuRAbNlfR1Tdw3javws+fP+JMQJPwvlXzuPqCGXx103721p9xOpyoFkyCyFXVnwF9qvpXVf0AcHmY4zIJqLw4l7qWLmpOnWuTb2rvoeZUZ9w1Lw0qjZHFg+pbusa1PZqJCBveerG/FEeFleIYRVALBvl/HheR60RkOb6Jc8aE1GD57y2Hzo1mOtv/EGcd1INK890camyP+vH5ee7AVXUKx7meeLTIzU7ju+9YxmErxTGqYBLEl0VkKvBZ4F+BB4BPhTMok5jme7IonJp+3nDXitoWUlzCkqL4mCA3XFlBNn0DSs2pjrF3doiqBlygKSPFxbo1ZQ5EFBpWimNswSwY9LiqtqrqblW9SlUvARZGIDaTYESEVcUeth46hdfr64eoqG1mUeFU0lOcm40bTudqMkVvyY0ndh3ncFMHb71kFkU5GQi+tTy+dssSbl5eNObx0cxKcYxuopWzPhPSKIzxW13soaWzj73Hz9A34OXlY7G/gtxoFuZlkyREbcmNrt4BvvrEPi6cOYWvv+Vittx5Na/cex1b7rw65pMDWCmOsUw0QcTXcBITNVYV5wK+8t/7jp+hu88btx3U4JtFPi83K2oXD/rxXw9R39rNXTcswhWjs73HMrQUx/eeOeh0OFFlooVt4n+VE+OIGe50SvOz2VLdRLq/NHSsLzE6lmgdyVTX0sWP/3qI6y6eyWULcp0OJ6xuWFrI3w42cv9fqllVnMuqhVaKA0a5ghCRNhE5E+DWBhRGMEaTYMqLPbx05DQvHD5FwZT0mB0pE6zSAjdHTnXQPWyegdO+umkfIvC5tRc6HUpE3HXjYuZ7svi0leI4a8QEoapuVZ0S4OZW1TGvPETk5yJyUkR2D9k2XUSeFpGD/p8BvxqKyDUiUiUi1SJy58R+NROrXElCd5+XzXsaaOnqZWNlndMhhVVZvhuvwqHG6Omo/vvhUzzx8nFuf91CiuI8QQ8aLMXR3NHHv/1uZ0IsBzuWcC7v9CBwzbBtdwLPqmoJ8Kz/8Xn8a1Dcj2/t60XAbSKyKIxxmiiysbKOX/295uzj7j4v6x/ZFddJojQ/G4iekhsDXuXux/ZSlJPBv1yZWAMWB0txPLPPSnFAGBOEqj4HDF8B5ibgl/77vwRuDnDopUC1qh5W1V7gN/7jTALYsLmK7r7zR5J09Q2wYXOVQxGF3zxPFikuiZqhrg/9o5Z9x8/wubUXOrrYj1OsFMc5kVsg1idfVY8D+H/OCLBPEXB0yONj/m0BiciHRWSbiGxrbGwMabAm8uKppEOwUlxJLMzL5mAUXEG0dvbxraequGz+dNYuKXA6HEdYKY5zIp0gghFoLN2IjYGq+lNVXamqK/Py8sIYlomEkTqk476jOt8dFXMhvvPMAVq7+rjrxsVxVxxxPHKz0/iOvxTHPY8lbimOSCeIBhGZCeD/eTLAPseA2UMez8K3aJFJAOvWlJExbNZ0rJd0CEZZgZtjzV209zj3bbXqRBv/7+81/PNlc7hw5hTH4ogW5cUe/u/rFvKbl47y+MuJ+REU6QTxB+C9/vvvBf43wD4vASUiMt+/Fvat/uNMArh5eRFfu2VJ3JV0GEvJDF9HtVPNTKrKPY/vITstmc++Kb6T8Xh8+k2lLJudw/oELcURthXgReQh4PWAR0SOAV8C7gUeFpEPArXA2/z7FgIPqOpaVe0XkY8BmwEX8HNV3ROuOE30uXl5UdwnhOHKCs4tHrTcgZnjm/c0sKX6FHfdsIhpWakRP3+0SnEl8f3blrP2e3/jk7+p5OF/uYJkVzS2zIdH2BKEqt42wlNvCLBvPbB2yONNwKYwhWZM1Jk9LZP0lCQONER+JFN33wBf2bSX0vxs3nX53IifP9oNluL4+EOVfPeZg/xrnDd3DpU4qdCYKJaUJI6V3PjZ869w9HQXX7phcUJ9Ox6PG5YW8vaVs7j/L9VsHbJeSbyz/w3GRInSfHfElx890drN/X+uZs3i/LMLNpnAErEUhyUIY6JEaX42J9t6aI7gh8+9f9xHv1f5wnVWrGAsmanJ3HdrYpXiCFsfhDFmfAYXDzrQ0BaR6qnba06zcUc9H7uqmNnTM8N+vnhwUdFU7rz2Au55fC/L7nmaM119FOZksG5NmSMDKzZW1rFhcxX1LV1hicOuIIyJEmdHMp0Mf0e116vc9Ye9FExJ5yNXJVa9pcmalplCkkBrVx+Kryy6E/XCNlbWsf6RXdS1dIUtDruCMCZKFExJx52eHJHFg367/Si76lr53q3LyEy1j4Hx+OZTB/AOa13q6hvgjt+/zCMRTBIvHj5FT3/gumWhuoqw/xnGRAkRoSwCJTfOdPexYXMVK+dO48altrTLeI1UF6yn38uZrr6IxTE8OQwKZd0ySxDGRJGSfDd/3H0cVQ1bLaT7njnIqY5eHnz/pQldb2miCnMyqAvwIVyUk8HGj5ZHLI7ye/8UMI5Q1i2zPghjokhZfjYtnX00tvWE5fWrT7bz4NYjvGPlbC4qmhqWc8S7aKkXFok47ArCmChSerbkRjszpqSH9LVVlf94fC8ZKa6Emg0caoPt++EcPRQtcViCMCaKlPmHulY1tLG6JLQT1/60/yR/PdDIF667EE92WkhfO9FES72wcMdhTUzGRJHc7DQ82akhH8nU0z/Afzy+lwV5WbzninkhfW0TvyxBGBNlSmaEfiTTg1uOcORUJ1+8fhGpyfZnb4Jj/1OMiTJlBW4ONrSFrJTDybZuvv+nat5wwQxeXxZolV9jArMEYUyUKc1309E7EHAI40R848kqevoH+ML1Vm/JjI8lCGOiTFmBb3W5UJT+3nG0hd9tP8YHVs9nvidr0q9nEoslCGOiTMngSKYTk6vJ5Ku3tIc8dxofv7okFKGZBGMJwpgoMyU9hZlT0yd9BfFoZR07jrZwxzUXkJ1mI9rN+FmCMCYKTXbxoPaefu59cj9LZ+dwSxSM1zexyRKEMVGorMBNdWM7A8PLhgbpB3+qprGth7tuWERSktVbMhMT8QQhImUismPI7YyIfGrYPq8XkdYh+3wx0nEa46TSfDe9/V5qTnWM+9gjTR38/PlXeMuKWSyfMy0M0ZlEEfGGSVWtApYBiIgLqAMeDbDr31T1+giGZkzUKBuyutyCvOxxHfvlJ/aS4hLuuMbqLZnJcbqJ6Q3AIVWtcTgOY6JK8YxsRMY/kumvBxp5Zt9JPnZ1SciL/ZnE43SCuBV4aITnrhCRnSLyRxFZPNILiMiHRWSbiGxrbGwMT5TGRFhGqos50zPHNZKpb8DLPY/tYV5uJh9YPS98wZmE4ViCEJFU4EbgtwGergDmqupS4PvAxpFeR1V/qqorVXVlXl5eWGI1xgml+e5xJYj/eqGGQ40dfOG6RaQlu8Y+wJgxOHkFcS1QoaoNw59Q1TOq2u6/vwlIEZHQ1j42JsqV5bt5pamDnv6BMfc91d7Dd585wJWlebzhQqu3ZELDyQRxGyM0L4lIgfjXQhSRS/HFeSqCsRnjuJL8bPq9yitNY49k+uZTVXT1DvDF6xfZMqImZBxJECKSCbwJeGTItttF5Hb/w7cCu0VkJ3AfcKuGqrSlMTGirGCw5MbozUy761r5zUtHee+qeRTPGN+IJ2NG48j8e1XtBHKHbfvxkPs/AH4Q6biMiSYLPNkkJ8mo/RCqvnpL0zNT+cQbrN6SCS2nRzEZY0aQmpzEfE8WBxpGHur6h531bKtpZt2aMqZmpEQwOpMILEEYE8VKC0YeydTZ28/XNu3noqIpvG3l7AhHZhKBJQhjoljpDDe1pzvp7O1/1XM/+sshTpzp5q4bFuOyeksmDCxBGBPFygqyUYXqk+c3Mx093clPnjvMTcsKWTlvukPRmXhnCcKYKFaaH3gk01ee2IdLhDuvvcCJsEyCsARhTBSbm5tFanISB4dcQWytbuLJPSf4yOsXMnNqhoPRmXhnCcKYKOZKEkpmZJ+9gugf8HL3Y3uZNS2DD125wOHoTLyzBGFMlBtak+m//1FLVUMbX7juQtJTrN6SCS9LEMZEudJ8N8dbu6k91cm3njrAqoW5rFlc4HRYJgHYSubGRLnTHT0AXLnhzwC8tsRj9ZZMRNgVhDFRbGNlHf/1wvnrad33bDUbK+scisgkEksQxkSxDZur6On3nretq2+ADZurHIrIJBJLEMZEsfqWrnFtNyaULEEYE8UKcwLPcxhpuzGhZAnCmCi2bk0ZGcOGs2akuFi3psyhiEwisVFMxkSxm5cXAb6+iPqWLgpzMli3puzsdmPCyRKEMVHu5uVFlhCMI6yJyRhjTECWIIwxxgRkCcIYY0xAliCMMcYEZAnCGGNMQKKqTscQMiLSCNSMuWNgHqAphOHEMnsvzmfvx/ns/TgnHt6LuaqaF+iJuEoQkyEi21R1pdNxRAN7L85n78f57P04J97fC2tiMsYYE5AlCGOMMQFZgjjnp04HEEXsvTifvR/ns/fjnLh+L6wPwhhjTEB2BWGMMSYgSxDGGGMCSvgEISLXiEiViFSLyJ1Ox+MkEZktIn8WkX0iskdEPul0TE4TEZeIVIrI407H4jQRyRGR34nIfv//kSucjslJIvJp/9/JbhF5SETSnY4p1BI6QYiIC7gfuBZYBNwmIoucjcpR/cBnVfVC4HLgown+fgB8EtjndBBR4nvAk6p6AbCUBH5fRKQI+ASwUlUvAlzArc5GFXoJnSCAS4FqVT2sqr3Ab4CbHI7JMap6XFUr/Pfb8H0AJOxCBCIyC7gOeMDpWJwmIlOAK4GfAahqr6q2OBqU85KBDBFJBjKBeofjCblETxBFwNEhj4+RwB+IQ4nIPGA58KLDoTjpu8C/AV6H44gGC4BG4Bf+JrcHRCTL6aCcoqp1wDeBWuA40KqqTzkbVegleoKQANsSftyviGQDvwc+papnnI7HCSJyPXBSVbc7HUuUSAZWAD9S1eVAB5CwfXYiMg1fa8N8oBDIEpF3ORtV6CV6gjgGzB7yeBZxeJk4HiKSgi85/FpVH3E6HgeVAzeKyBF8TY9Xi8ivnA3JUceAY6o6eEX5O3wJI1G9EXhFVRtVtQ94BFjlcEwhl+gJ4iWgRETmi0gqvk6mPzgck2NERPC1Me9T1W87HY+TVHW9qs5S1Xn4/l/8SVXj7htisFT1BHBURMr8m94A7HUwJKfVApeLSKb/7+YNxGGnfbLTAThJVftF5GPAZnyjEH6uqnscDstJ5cC7gV0issO/7XOqusm5kEwU+Tjwa/+XqcPA+x2OxzGq+qKI/A6owDf6r5I4LLthpTaMMcYElOhNTMYYY0ZgCcIYY0xAliCMMcYEZAnCGGNMQJYgjDHGBGQJwiQ8EckVkR3+2wkRqfPfbxeRH0YohtcPVowVkRsTvbKwiQ4JPQ/CGABVPQUsAxCRu4B2Vf2mg/H8gQSesGmih11BGDOCYd/q7xKRX4rIUyJyRERuEZFviMguEXnSX6IEEblERP4qIttFZLOIzAzwum/zryGwU0SeC/D8+0TkB/77+SLyqH/fnSKyyr/9XSLyD/+Vzk/8peuNCSlLEMYEbyG+8t83Ab8C/qyqS4Au4Dp/kvg+8FZVvQT4OfCVAK/zRWCNqi4FbhzjnPcBf/XvuwLYIyIXAu8AylV1GTAAvHOyv5wxw1kTkzHB+6Oq9onILnylWZ70b98FzAPKgIuAp33leXDhKwU93BbgQRF5GF+Rt9FcDbwHQFUHgFYReTdwCfCS/zwZwMmJ/1rGBGYJwpjg9QCoqldE+vRcnRovvr8lAfao6qhLcarq7SJyGb6rkR0ismyccQjwS1VdP87jjBkXa2IyJnSqgLzBtZpFJEVEFg/fSUQWquqLqvpFoInzS84P9yzwf/3Hufwruz0LvFVEZvi3TxeRuSH+XYyxBGFMqPiXrX0r8HUR2QnsIPAaARv8ndu7geeAnaO87CeBq/zNWtuBxaq6F/gC8JSIvAw8DbyqM9yYybJqrsYYYwKyKwhjjDEBWYIwxhgTkCUIY4wxAVmCMMYYE5AlCGOMMQFZgjDGGBOQJQhjjDEB/X/9lbHHhZJwEwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(len(q6clusters_ase)), [max(x) for x in q6clusters_ase], 'o-')\n",
    "plt.xlabel(\"Time slice\")\n",
    "plt.ylabel(\"Largest cluster size\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the results are identical for both calculations!"
   ]
  }
 ],
 "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.6.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}