Structural_optimization.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Structure optimization\n",
      "======================\n",
      "\n",
      "This tutorial is concerned with structural optimization. Except for high-symmetry structures - like the cubic crystal we are using in this tutorial the determination of the correct structure for the investigated system is a crucial part in the workflow of the computation. The Quantum-Espresso package is able to do a full structural optimization of the crystal and the QE-util provides an interface to this mode of operation. To continue skip over the initial setup to the second cell below."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Import the basic libraries\n",
      "\n",
      "# ASE system\n",
      "import ase\n",
      "from ase import Atom, Atoms\n",
      "from ase import io\n",
      "from ase.lattice.spacegroup import crystal\n",
      "from ase.units import GPa, Bohr, Rydberg\n",
      "\n",
      "# Spacegroup/symmetry library\n",
      "from pyspglib import spglib\n",
      "\n",
      "# iPython utility function\n",
      "from IPython.core.display import Image\n",
      "\n",
      "# Import the remote execution tools from the qe-util package\n",
      "from qeutil import RemoteQE\n",
      "\n",
      "# Access info\n",
      "import host\n",
      "\n",
      "qe=RemoteQE(label='SiC-structure',   # A name for the project\n",
      "               kpts=[3,3,3],    # k-space grid\n",
      "               xc='pz',         # Exchange functional type in the name of the pseudopotentials\n",
      "               pp_type='vbc',   # Variant of the pseudopotential\n",
      "               pp_format='UPF', # Format of the pseudopotential files\n",
      "               ecutwfc=70,      # Energy cut-off\n",
      "               pseudo_dir='../pspot',\n",
      "               use_symmetry=True,\n",
      "               procs=4)         # Use 8 cores for the calculation\n",
      "\n",
      "print qe.directory"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "calc/SiC-structure.GHl7tH\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Preparations\n",
      "--------------------\n",
      "\n",
      "To start any calculation we need to create a crystal. We will again use a cubic, zinc-blende (F-43m) SiC crystal again with a lattice constant determined in our first tutorial. But this time we will distort it randomly to have a low-symmetry structure with many potential degrees of freedom and configuration displaced from the equilibrium position. Our task is to find the equilibrium structure - which hepefully will be a high-symmetry cubic crystal with a F-43m space group symmetry."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a=4.344\n",
      "cryst = crystal(['Si', 'C'],\n",
      "                [(0, 0, 0), (0.25, 0.25, 0.25)],\n",
      "                spacegroup=216,\n",
      "                cellpar=[a, a, a, 90, 90, 90])\n",
      "\n",
      "# Assign the calculator to our system\n",
      "cryst.set_calculator(qe)\n",
      "\n",
      "# Verify the symmetry\n",
      "print \"Symmetry group:\", spglib.get_spacegroup(cryst)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Symmetry group: F-43m      (216)\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"Stress tensor  (Voigt notation, GPa):\\n\", cryst.get_stress()/GPa"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Stress tensor  (Voigt notation, GPa):\n",
        "[ 0.054  0.054  0.054  0.    -0.     0.   ]"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let us distort randomly the positions of atoms in the unit cell as well as the unit cell itself. We will use standard deviation of 0.05 (5% of the unit cel size). This should generate a completely non-symmetric crystal (group P1)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cryst.rattle(stdev=0.05)\n",
      "cryst.set_cell(diag(1+0.01*randn(3))*cryst.get_cell(), scale_atoms=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Verify that indeed we have a low symmetry structure\n",
      "print \"Symmetry group:\", spglib.get_spacegroup(cryst)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Symmetry group: P1         (1)\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Display the structure\n",
      "ase.io.write('crystal.png', cryst, format='png', show_unit_cell=2, rotation='115y,15x', scale=30)\n",
      "Image(filename='crystal.png')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAL4AAADVCAYAAADpR0oEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXd4FNXXx7+zu7MzW9IIEAJp9FCTSAeBgAiCgCIo0mxg\nRcGOov6Ir1jAgoqIiiCC2GmKUhQB6SDSIUAgIQkEQiCQZHenz/vHZCVs6ia7O7PJfp6Hh2Qye+/J\n5uydc889BQgQIECAAAECBAgQIEBtglBbgACq0RpACwCxxf/sALYD+FNNoXyFQW0BAniFegCa47pS\nCwA+crnnFQANAJwt/mcH8COA5QAe9pmkKhFY8f0PA4AmuK7UEQDec7nnfQDJuK7UhwB8VYWxZQDt\nABwDEAzADOCCJ4QOEKAyzADiAQwC8AiAIJefPw0gC8A2AMsAvAlA56G55RJfDwFwGcBLHho7QB2G\ngGKGJAG4E8rKXZKBABgAp6DY2gsBNPKhfLLL980A3F/8NQFA70NZAvgROijK3BOKjV2SKACFAK5B\nMT1+BdDZ5R49PLd6V4eUCn42AMBBALf6RpQAWoKC4glp4XKdgGIXs1Bs4t0A7nW5Rw8g1NsCehEC\nwEgA3yKwN6x1BEPZMLryNYDzUBQ7HcDrZdzTFADtPdE0RwsAnwOIVFuQABVDAAgv4/rzAA4AyAdg\nA7C6jHvaAohGwMYtSSiA2QBOIOAWVxUDStvHwwFsgPLHcUBZtV0f0+0A3ASgfhk/C1A5phL/T0Bg\ncfA6HQB8A2ArFF81C8XcKEk8gNugrNhWn0rnn6TU4LXRUP4WB3D9wxDATcIBfADlFHEvgFwAQ13u\niYbiakuGovCkD+Wrrbi6M92FANClxPeuC1GdhgAwDcCnAH4DcATAky73BAF4AcA9ALpB8WWr6ear\nK9RU8UsSCeAigAUofQinOt6waUdCMUOcR+q7AEx3uScFQB6uH6mfAVDkBVkCuIcMz+pEKIBHoYRQ\nCFAWL8mD41cbd3/JZCj2slOpr0H5xUryMhQb7yyATCgby8waSRnAV3ha8UsSCmAPgFkAFgMQvTSP\nVzgHJU7kQwDjcKM9F8D/SfHy+F0BrIVvwzA8ggzFHh+mtiAB/B4DlAOwjmoLUhVcNz9vQjnFDLgK\nA7iLAcoiehFKUoymcVX8WCi+9P+pIEuA2oGlxNd3Q6OLaHnuLueGaCiAwQicgPorKSrOrQewFMrp\neoKKcpRJSiU/HwjFi/O290UJ4AU86cevLp1xPdgvWk1B3IWEkg4HKIkMrskWAbSLFhTfCQFgB5RY\nq5Yqy+I290FJX3tabUECVAktKT6gLKKTcT2px69M6BgA/Yq/phEIYdUyWlN8V9ZCJS9iSg1fPwnA\nUSgb4ADaI0VtASohFsoGeLivJ/ZE9N4wlI7dCRCgOkyDj7yInn4UDoUStxHl4XED1A2GAUgFMNHb\nE3la8YOhnP6u8/C4AeoOJK4nvvSBlxZRb21+nLH2TaFEe9ZkA6yDsqFuC6V+TQco7jBzTQQM4Bc8\nB8WLON7TA6d4ekAXWgPYBKVKWFXttggAY0mS/NhsNu/X6/UOk8lkCw0NLQgPD78WFhZ2LSgoqFCn\n0/E0TZ81mUw/A5gCoJMbc9QVUtQWwAPE4Lr7syHKWUS1+IcnoOzeM6DUsWkDJY/T9Z5eNE2/IAjC\noBYtWnDNmze3RkdHE1FRUbBaS3u7BEFATk4OsrOzkZ6ezhw/flzkef4CwzCzoNSKsXn1t/IPvBmP\nrwbvABgBxZmyvOQPtP5LJkCx/9dD2cBIBEGMpijqLYqiGvbv39/cpUsXwmx234qRJAknT57E5s2b\ni9LS0nQ6nW4xx3HToSTX1FX8TfEpKD79yy7X10AxdZtAWdBOwyV3xB9+yWAoBUw3UxS1xGKx9Bw9\nerSlVatWIAjPiJ+fn49169Yx+/fvt3EcNw7KB60uojXFp6HUEi3JWwD6Q7EKwqFUgXbNAuwFpUBB\nVvHrS/1e/pDAXUgQhJ4kyZNJSUnJPM9b8vLyIEmeS90MCwvDmDFj6EmTJoVbrdYVFEV9AyDEYxP4\nD2VViPMWZX3A7gPwC5QanVehxOq4sg7KJrYLFG+Oq9IDSoOLUyj9oalw8opIgW83QDqj0fiZ1Wod\n+8ADD1hiYmKQlZWFrVu3YvTo0dDrPV+7iGEYrFy5ktm/f/9ljuN6IpAv7Cl6QykH48zXjoSSh1sy\n97YXlA2pM1/7MjzjSSy14rur+L58FOopivquYcOGQ5544gmLyVS6RlFhYSG+++47DB48GNHRno1g\n3bx5s/j777/nFyv/KY8O7j+YoOyzbgIQYTAYggDoRFG0ybJ8FUrV531QzNGJuK7UsVCqK58sMdZN\nUFZpp1JnwneVNfxG8QmKohZHRkaOeuKJJ8xGo7HMm0RRxK5du7Bp0yY899xzKOvDURN27NghrVq1\n6jLHcZ2g2Iu1HQJAT5qmHyMI4maWZaPCwsLsZrPZpNPpSJZl4XA40KpVK5hMJj49Pd2ek5NjIgjC\nJsvyFUEQ/oDiPUmD8n5ppZJCClwsFU0qPkmSr9WvX3/a1KlTLTRdeQFiSZKg0+kgSRK2bNmCHj16\noCqvqwp//fWXsG7dumyO49pCqb9Z2xgN5fzkFp1O15miKP2tt95KtmjRQhcZGQlJkrBmzRqEhYWh\nXr16CAsLQ0RExH/vryRJuHz5Ms6ePYvt27cXZmdnA8AXPM9/AsUlrUm0aON3NBqNu6ZPn24KDXWv\nnDzDMFixYgWOHz+O559/HiEhntmffvnll/aTJ09+yXHcVI8M6Du648Y6SMeh1LVxogPwN0EQXRs2\nbCjddNNNVOfOnREeXlZR6aqRm5uLrVu3srt375YJgljOsuxkaNBFrCXXFQCQFEUdGTFiRMvu3btX\nW7ZLly6hQYMGAIDs7Gw0adKkRq7PoqIizJw508EwzC0AdlZ7IM/SFEqlZ6dS81A6GZbkfShVoJ0V\n6w5AsckBII6iqO/q1avX4f7777c0auTZUjclnASadBFrSvFJknw9Njb2+cmTJ5s94aPnOA5z5syB\nxWLB+PHj4e4TpCQHDhzAt99+m8VxXGt43+QJgWJ+lPSAPOdyz3QonhKnUh8HsKoqg+t0ukl6vf7D\ngQMHUv379zd4wzvm5MSJE1i6dKmd5/mVLMs+AqWtqOpoqZZ5Q51O9/PkyZNNntqk6vV69OjRAwRB\nICYmBnq9HoIgQKdz//iiUaNGOHPmDHnlypVCWZZrsuoboMSTJADoCyWp4ghu/DCNBfB/UJKtJSge\nkL240bW3FUrXwzVQYptSqzA3QZLkO0FBQTOmTJliSUhI0FXnvXCH+vXro2fPnmROTk7Lq1evDhVF\n8QcoJd19SQqAzSUvaGbF1+v1ryYlJU0fP36812qri6KId955Bx07dsSAAQPc9gKlp6fjs88+y2FZ\nNgrlFz+1QFFs52r9J5Qjcyf9ACzB9ZX6LJSSjLluCeM+hNFo/DA0NHTiU089ZQkK8m0BY0mS8PPP\nPzP79u07w7JsDwAFPpy+xie3KR4T5Ub0er1+at++fb3aUECv12Py5MkoKipCYWGh26+PjY2FxWIJ\ngnJaeCdKl79rDOXQZTWuny66RsxtgrKS3wyl/uh0eF/pQZLk/4KDgydNnTrV50oPADqdDnfffTed\nlJTUnKKoP6Fy8wituDOHNWrUaNlLL73k07+ILMv49ttv0bFjR7Rv3x6yLKOgoAD5+fkwGAylDsX+\n97//weFwQBCEIlmWN0HpQ1uyZxZR/E8TpbBLMMhisax48cUXzZ7ydFUXSZKwaNEix8mTJ7/hOO4R\nH02rzQMsk8n0+/Dhwwf36NHD00NXiCzLSE1NxerVq0HTNLKysmCxWBAWFoaEhAT079//hvtZVjFN\nX3nlFVYQhEgoDeO0TojRaDw9ceLE8NattVGi0m63Y+bMmXa73T4EwBYfTJkCF2tFE6U+JEnq0ry5\naz9k70MQBNq0aYPWrVvjl19+weTJk0GS5XcUoigKABAREcGcO3fuJgAbfSRqtaEo6pOEhASLVpQe\nAMxmM8aOHWtesmTJ9xzHtYD3cyFSXC9oITozXBTF4Pr166smgE6nw+bNmytU+pI0bdrUBCWDS+sM\nIEnyrrvuuktz/Xfbt2+Pdu3ahVAU9b4a87ur+N4IW+0UGRnJeNut5kliY2ONZrO5r9pyVIbJZJp1\n1113mT0VvuFpRo4caRJF8QEADXw9t+peHYIgOhevoH5DdHQ0RFHU+orfDkCbjh2123fBarWiY8eO\nkl6vn+TruVVfZkmSbBwWFqZ6q85BgwZV+d7g4GAIgqDpRBWKoqbefPPNpMGgiW1cuSQnJ5v0ev0z\n8O5haorrBdUVX6fTWatqW3uTwYOrXtXQaDRCkiT1hS4fqyiK43v16qVtrQcQExODsLAwGt4tKznD\n9YLqpg40dHrsJlqWu0/jxo35msQm+ZIePXoEURQ10pdzuqv4pT45NUWSJLsgCJ4e1qsUx/vwastR\nHgRBdG7evLnfFNCKiYmBTqfz6SGO6qaOIAiXioqKtHbSWSE2mw16vV6zdXhMJlNybGys5s0cJ02a\nNAHLss0AlJ1qV3NKeSNVV3xJkvZnZGSo3tV87dq1Vb43OzsbBoPhqBfFqRGCICRGRflPHV6KohAS\nEsJA8UR5gxTXC6orPoB92dnZpCyr25Ng/fqq50lkZWWJDMP44qi9OoRIkmStSRaVGhR/UNv7aj4t\nHGBliaIoXLumuey0cjl9+rRNkqQ9astRDkEkSfKeKrblK8xmswE+7HSiBa+OTJLk4YyMDC8M7XlE\nUcS5c+coXE/h0xq0wWDwqz0TABiNRj18GKqsiQ2Q3W7/atOmTYn16tUz8zwPnueh0+lgMBhAkiSC\ngoIQEhLisZKBNeHo0aPQ6/UnoPRj1SKCJ6vM+QpRFCUoecPeIAUaiM40QMn872w2UT2NBkPPIgfT\nwl6YT6795WeYKCNoioQoymA5HgzL4UJePgACMTHRiGwShaioaERHR3v0w1DVk9tNmzYVOhyOWZXf\nqRoOQRC0lFJaJTiOk1BByb8aMgMqKT4BoH+I1fyCjWH7N6wXzHRq24zo1rGlNalNU3RsFYMgS/lP\nOVmWcS73CvYfz8D+4+n459BurPj5RxCEDp06d0b3Hj1rVBIDqNrJbW5uLorrxiyv7F4VucLzvIFl\n2f/CqP2BS5cu8QDO+Wo+dxU/Be7Z+aE6HfGAmaaerx8WFDJl3GDLXbd2I0KDLG4d9xMEgaiIcERF\nhGNYshIbJssyTp3NwcKVm/Hxhx8gLi4O3Xr0Qnx8fLWSyavC5s2bWQCfw/fJ0u7A0zR9Jjs7u7Ua\nOQ7VQZZlnD9/noYP903eysBqbjXTKbwgjhrYs6P4xL2DLN0TWnrNRrczLJb/sRvzf/gTl6/Z0L1H\nL/Tu0weeDNDKzMzE3LlzC3iebwXgoscG9gJGo3HhkCFDHkpOTlZblCqRl5eH2bNnX+E4zls+WK+X\nCdcZDPqnTZTx0OQxg+49vPI9euk7T1l6JHquln1ZmGkKE4b1wfalr+PbdybDfjkTH37wHjIzq17o\nuKIDLJ7nsXjxYhvP849C40oPABzH7UhPT9fsybIrWVlZMBgMB704hVdPblsEmem9bZtFzdy69P/M\nLz88wtAw3LeRuwRBoFO7ZvjxvalIeexOfLVwAX5b8yuqEgtU0QHW2rVrOZvNthXADx4U15vsOHXq\nFOEv3p3U1FTG4XBs8OIUKa4XPHGA5VzlD06bdGfCpq9mWFrEeLYcnbsQBIFRA7tj97czoWfz3V79\nS3L8+HFs3bqVYVn2AWi/5b2T47Ispx87dkxtOSqFYRj8+++/hCzLS3w5b3U2tyUJCTLTa5tGRXRc\nNPNxs9oK70rD8BB8N/spLP9jN154fwH6JvdDcr/+lb+wmNOnT+Orr76y8zx/G/zAxCmGABDlcDj+\nWLF8ecs9u7Ybc3MvgeN4cDwHnhdgMBhgJEmQRhL169VDoyZRaNIkCjExMQgPD/eac6As9u7dK+v1\n+o08z/v0XKQmu78GVjO9ddTA7nHvvXAfpddrIeynNM7Vv2diawx7cjZsNhuG3D600j3HiRMnsHDh\nQjvHcXdCO4Viy6OBQa97KMhiGu5guQ6kQW/o0DKG757QytipXTPExzWG2USBNpIwGknwggCGVc5I\nTmddxP7jGfjn2FH8se43FNociI2OQnRsU3Tt1q3GbuKKkGUZmzZtsjEM867XJlFIQQ3r4ztpbDFR\nOx+959bI1x4bSWrhRLUqXL5aiDunvI+wiCYYcdeoG5R/7dq1GDx4MGRZxs6dO+WVK1faeJ6/HcDf\n6klcIQSAbkEW03McLwwbltxJuvOWLqakNk3RuEFYtZ0JefkFOHjiLDbsOIzv1273qpv40KFDWLZs\nWUZxSLI3zcgaF5QCgHCLifpn6oQhUS8+dIcmQh7cocDmwPAn30W9RjEYNvyOGxQkPz8f33zzjS07\nOzuLZdlRALQYemwhCIwJspimmWkqcvKYQaZxQ3vr6oV4Pr7rBjdxgR3duvdEt27dYLFYajy2zWbD\nzJkzHQ6HYxCUArjepMaK/5bVTN/x4J3JLd6Ycq+3kga8Tv61Igx67G20bJOAWwcOhCiK2L17t7xq\n1SpGkqR3BEF4C4AW08IGmWnjsq4dWlJPjbvN2q9rO5/Y47Is499j6fj8541Yv/0ghtw+DF27dq2R\ni3rx4sWOY8eOfc1x3OMeFLU8aqb4Jsoo33lLF/bT1yZR/mLelMfFy1fR/6E30LBxjHjixAlWkqRU\nhmHuh1KyW2uEWM30PBNlHLHg9UfNyV29la9ROUdOZeGR1xdAT1kx8u57qtVz4PDhw1i6dOl5juNa\nwjf18kspvjvBTKMahgffs+LD5w0k6XcWTimsZhpJbeLw6TerBYeDuVVQnP0/qi1XGQwy08YtIwZ0\nTfjpg2dNrZs2VlWYhuEhuG94b+RdvoK5X34Ds9niVseZrKwsLFiwwM5x3HAAZ7wr7X8QqGZ9/AYm\n2nhq9dwXQ7p2aOFxqdTk+XeXsD+s3bGm0M6MhLYqJ9BWM73ARBnvUnuVL4+Sq//oMWNhtVa8z8jJ\nycFHH33kYBhmHICVvpGybKpkIAZZ6IUP3plsqm1KDwCvPzmasphpb9Z0qQ5BVjO9uU/nNiP3/TxL\nk0oPAO1bRmPL4v8hOTEWn86bi6tXr5Z7b0ZGBj766CM7y7IPQWWlB6pm6oxqWC/4mW9mTTGRBr8L\n864UI2nATW2bksvWbAOAuVC/pWe41Uxvu6N/lzYLUh41mSht+xD0Oh36d2sPQeDx6aLvEd+m7Q1e\nH0mSsHnzZvHbb79lOI67BxpQeqByxQ8x0cZN37/3dFBcE5/X9fQZ0Y3q450vVyHIYorleEHNWPtg\nq5necd/wPi3ef+E+yp8K6Xbv2BJWM4U5n3+Ddu07wGw2Izc3F/Pnz7cdOnToCMdxtwDYrbacTipU\nfJ2OePTWHh0HThk/RMvl8jwCw3LYczitpSBKX8J3reZLYrKa6c0jb+0W/+7zE/zSa5YUHwfSoMe8\nRd8h7/IV8YcffmAKCgpeEQThYQBXVBQtBS6b24oUn7Ca6J9nPzchPCZSvdr1viK5azukn8sVTqSf\nZ0RJ8nXpECLIQq8e0KND1/mvTaL9aaV3pUv75igotGH1+i1FDMt1lWV5DdQP7tsElwDLit7h5Hoh\n1tCeia28K5KGePSeW2mDXj8FPs5FJghiXP3Q4D6fz3jE5M9K72T6IyPQI7GVwUQbH1RblvIo910O\ntppeeGrcYIs/PnKrS8dWsWgeHUECGOrDaSNpIzn/qzefsFDG2mFREgSBea9ONBt0uicBdFZbnrIo\nT/Eb87zY757BPeuO1hfz1LjbgkKCzNN8NB0RZKaXPn7vQCoxPs5HU/qGiPBQfDDtftpqpn8CoHbW\ne9UysIyk4eFRg7ojuILKB7WNtxcoXrY7+neBJEkJALyeqU0QxLj6YcHdp028o3Ys9S6MGtid6J7Q\nsqGJNr6hsigprhfKVHyLiRo6LLmTNhsneYlZC5V2tZSRRHKXdiIAb/e4qnUmjitaNnnKUnzC5mDb\nJbaO87UsmqF7QiurxUT19OYcJtr44vhhvY21zcRxJSI8FK8/NZoOsZq9nWziFmUpfjOLiZJ9nSiu\nJRLj40AaDN5UfJMsy5OeuHeQto9lPcS9g3sSgih2BxCntixOylL8TonxcVqMRfcq0ybe8d/XHVvF\noMjBtID33Jp3d2rbDE2jGnppeG1hpimMH9qboCnySZVESHG9UErxKSPZrXtCK5+Va9YKLz884r+v\ngywmNKwXzECp8elxgq2ml54ce1udeo8fvnsARYB4BOp4eCpv/mamjX2S2sT5/ylKDenUthkB73Qv\n70waDLEDeyZ4YWjt0iKmERJaxwLAKLVlAcpQfF4Q41rFRqohi6bo2CrWatDrWnt6XKuZfubxewdq\ntiqFN3ly3G1BwVbTS2rLAZSh+IIo0mZa7fMG9THRRpCkIcjDwxKCKA4bfVvP2hffXQVu65UIjhNa\nAfB18FflB1iSJJMUVTv9yhXhPMByQhtJkHp9zcsJ3Ei00WAwREXU8/Cw/oHBoEeb5k0c8I4JWREp\nrhdKKb4oSnqqFuTUuovzAMuJ0WgAoSM8fXTdqUOrGL/rT+VJuie0MusIQvXDrFKKr9frRE4Q1ZBF\nU/C8CFmSPZqNZSQNXXsktvL0U8Sv6NS2GRkSZPb2qXillFJ8nY7gGYZTQxZN4WA58KLo0VLbFhPV\nJ6lN0zpp3ztJjI8DxwtJastResXX6Vg7G1B8B8OB5wVPZmIRDpbrUNtDFCqjWVRDiKIUDMCXjXhT\nXC+UUnwjaTibdjbHJ9JoiZIntwBwJC3LJojSSQ9O0VCv05GNG4R5cEj/Q6fToWVspANABx9OW/kB\nloPhtu5PzVA7VcznlDy5BYB9R09L8GxPpmCLiRbq8sbWSWiwGQA87Sp2i1KKz3D8rl0HT6mRbK0Z\niuwMLl6+RsOzRWNpmiL9o0WJlzHTFAEfNnMui7KOD/cdSE2ve8eKJTh8MhMWE30Gnm04TBtJQ517\nkpaFiTLqoEHFP11oY3R5+QU+F0ZNSh5gHTiRAV4Qdnh4CkmSA3oPAKLSm8uXPvMqpR5KFhN1/EBq\nhvfF0RAlD7B2HTxlszlYTys+w3J8wMAH4GA4Cb6tWJfieqFMk8bOsGt/+3u/lpsYew2OF7BpzxEd\nPN8JxcFyfJ02IZ3YGU4GwKgpQ5l/CJYTPv9+7XYU2VWVTRXWbN4HAsRRAJ50ZQLApWuFdkoInIrj\n/KUrBIBcNWUobwXKMhoMW39av7POGaVzl60tvFZkn+WFoQtpyph3IsOnzf00h4PhkH3hignAYR9O\nm+J6odxH77Ui++yPl621yXVkQzZt4h04djobJzLOiwBWV/qCaqDX6f6pa3snV46mZcFiprLgW1On\n8gOsEmy8dKWgcM/hNC/Kox1efngEPv/xD1aUpHmo3I1pBBAGoDGACADBqEKvgauFti37jp6pk3sn\nJwdSMyBL8i615ajojyU5GO69ed+t90WPItUptDnw47qdMssJ811+ZATQG8DTJpNpBU3TWTqdzk6S\nZA5N06coiko3GAyX9Hq93Ww2HyJJ8lMAEwC0LGOafbsOnap7G6cS7DmSZi+wObapLUeFgfeiJH31\nx45DKfuPpyOpTVNfyaQK73y5iiMN+t8dLM4VX4o2GAxPEATxeGhoKNG8eXNjXFwcHRUVhcjISOj1\n+huiLO12O7KzsztkZWW1P3PmzPgzZ84YAKQ5HI5ZAH4GwALYn5Z5wSwIIgy1sMlGVdh7+LQAz4aC\nVItK/coEQYyNjaz/xe7v36q1Fb/2HknD8MmzCxws1wJANE3Tb4ui2KdLly7o3bs3HRnpfg6yKIo4\nevQoNm3aVJidnQ0An/M8/0aQhT668qMXojq393qFQs1x5VoRWt8+leMFMQS+tfFT4LLBrcqyc0SU\npGSbg41J7tqu1i1TDobD7Y+/bb9yrWg9SZJtSZL8bOjQoa0nTJhAJiQkGIKCqhdLpdPpEBERge7d\nu1OJiYmUzWa7KS8v73GO59exHN96aHKn2rmKVMCXyzdKuw+nreJ44TsfT73Z9UJVTxIbmSjjyd8/\nezmotpk8r3z0Hbdo5V97HAx/c3x8vH3s2LHm4OBgr8yVmpqKJUuW2CWRp1J//VAfGlx3krEkSUKb\nYc/YLl6+NgCApje3JbnAcPxjD77yqY3lPBm3pS57j6RhwU8bRUnWdQaARx991GtKDwDx8fF47bXX\nzGGhYfpv1ni7i7222LTnKBwMlwON9MGq8hG6LMvfXb5WuOOlOcu42uDbz8svwKhn5nBGmtY/++yz\nNIAataivKiaTCXePvhef/fgnJKnuRCnPXbauqMDmmAX12wIBcEPxAciFNuaen9bvynx7wUq/Xvav\nFdnRa/xrnASd7rnnnjNGRERg0KBBPpu/adOm0JMUtuw95rM51eTs+UvYdfCkDoCvbXsnKa4X3N2s\nMhwv/HggNWM8ZSQtXTu08LugK5uDRd8HZjBXChj9M888YwgPV1I/W7Ysy+3uHQiCAAgdfvvzb4y7\n/WafPGnUZMYnP7LHz5z7ihfEVSqJ4Fbzt7JIAXDR5mB7vbVgZe67i37h/cnsuVpow62T3rBlX8zX\nT5o0yVC/vnrdHLt27Yrcqw58tWqzajL4gm3/puKn9TvtdoZ9VW1ZSuKu4jtjHs7aHWynD5f+lv3q\nx9/7hc2fe/kaBkx8w3bm3KXcrl27Sc2aNVNVHr1ej7tHj8Hr85cjMydPVVm8hc3BYuJr8+12hrsf\nQL7a8pSkJqZKjs3Bdvl69Zbjgx550372/CWPCeVpfvv7X3S9d7r97PlLv1MUHTF8+HBNFAeNjIxE\n7z598cTMRfCHxcNdXpv7PWtzsL8B+FVlUarW/M0NLhfZmc4HT2S82WPcq46Fy/+StPQHzL9WhAem\nz3M8MuPznKuFttsJnb7fhAkTzEZj6UYka9euVUFCoF+//jh/uajWmTzb/k3F979vtxfZmUfVlgXu\nhCW7gcBE4VAqAAAW1klEQVRywlt2B9t5xrwfU7Wy+v/2979IGjXN/sfOQ1/bHGwLAKENGjSgWrRo\nUeb969ev962AxZQ0edKzVc3N8BjXiuyaNXGcuKv4pR4ZJThWZGcSnKv/R0t/l64WeLQCX5U4mpaF\nCS/Nda7yt9kc7OMA7CaTaVq/fv1UreVSHpGRkRh022AMf+pdXMi7qrY4NcLOsLjjydm2QjuzBOqb\nOOVSHa9ORThX/07vL/51dfywp5lHU75weDv5guMF/LxhF/reP6NwwKQ38jfsOPhe8SrvPB5tJcty\nQmJiolflqAk9e/ZCwk1dcfsTs5B/zT/LGnG8gNHPzrGfOpvzu93BTlZbnorwVtBZHssLPwiitOBk\nxvmiH9buSPh5wy7RTFNU85gIkAbPlCHPzMnDnCW/8Q++8im3YcehA2dz8l4QBPFhUZQ2okQyicFg\neLFXr1692rZtW+7vK8uyT335ZdG0WTPk5l3G3K9XYMQtXWExaWIPXiUcDId7X/jQ/u+x9L9tDnYU\nfFs+pDJS4BKo5quTEz2AISFW84s2B9s9JjLc3qVDC7JbhxamxPg4tG0eBZqquPPl5auFOJCagX+P\np8u7Dp4qOpCariuyM4Rep19qZ9iPABwv77Vms/mf8ePHd2rb1iu93DyKLMuY8/670Mk8fp33EqIb\n+bK2avUosDkwYsq79hPp59cV2Zl74dlCXJ5Ahouuu6v4Kajc3KkME4COADoHW0w3EwTR3cawUY0b\nhNlDrGbZRBtBUyREUQLD8nCwHC5evmootDF6i4k65mC4rQzH74KSzHAaQGUBL4TBYLDNmDHDVN0Q\nY1/z119/IS/zBM6ev4R5r05Ectd2aotULkdOZeG+6Z/YLuZd/dHmYCeh8r+HGtRY8UsN4CFMAOKh\nFBI1AaChPCoZKIWHclE1JS+L5maz+dBbb71l9pCsXufkyZPYvWUDXp40HE+/sxi39uyI/3tyNIIs\nqlbduwFeEPD+4jX8R0t/51ienypJ8iJoJACtDErprVZibRwA9kMp4rQeSpWDNQD+BLAdwClUfyXp\nFB0d7VcNq6OionDsdBb6dW2P7ctmghdE9Bz3Kjbv8WQN2+pz5FQWeo17zTbvu3U7HSzXRpLkhdCu\n0peJVhTfm0TWr1+/0l2iWgdYZWE2m0EQBIocDEKsZnzyykTMmXY/nnxzIZ6ZtRhquIkBZQP7zper\n+FsnvWFLy7owtdDGJAPIUkUY9/D4ya0/QJMkWan3Sq0DrPIwkiRY9voecUCPjti+bCYAIGHkC3jq\nzUXwVY2e9OxcvDznW77F4KfkT75dd8APV/kU1wvu+hUrOsDSKnq9Xu93cb+EjnBWFf6PEKsZc6Y9\ngOkP34Vv1vyNCS/NRUR4CCaOvAUjbulSqWfMHURRwoYdB/HJt+sK9x07QwD4kmH5swD6AMj22EQq\n4a7ip3hDCC/DcBwnwntnFl6B5wXQ5VS1aFAvGM/cNxRTxg3Bhh0HsXD5X3j14+/Rv2s7JLZpisT4\nOHRsFePWZphhORxNy8aB1AzsOZJm/3PnYVkQxIzicoo/4XpVhA+L/28K4GkAbwDwu/DSutDQtqCo\nqIiHUhjKLxBFESzHwVzJAZZer8Pg3kkY3DsJGecvYdu/qTiYmoFVG/fgaFoWmkTUQ2LrOLRu2hhm\nEwXaSMJoJCEIIhwsB4blkHrmvH3P4TQh68Jlk8VMZcmSvKu44NM2VFzfsgCKp2QHgLYAtOxASIHL\nol2bFN9pzpS0O4cAuC8zM7PSpc+XqYeVceHCBURF1Ic7dYziGjdAXOMGGD+0NwBAEESkpp/HgRMZ\nOJ15AXlXC8GwPFiOh5E0wEga8P3v2/lCOzMXwAoAh64VulUe+zKAKQCsUJSeAjCseCyt+fJnoBYp\nfjsAzwCIBRBT/K8dgDMl7skEMC8/P78Hz/M6kixfkQYPHuxFUd0jKysLNW0LajDo0b5lNNq3jC73\nnqiIcP17i39pWmhj9tRgKmdgUSMALwJ4AUAyfNv4wW08HaTmKeoBeAvAMiiBZpkAXDXTBuAfAO8D\nuBNAfdyo9ABwBMAPFEVl5eSUbmHKMAyuXbuGvLw85OTk4OLFi7hy5QoKCwshCOo9uc+fy0aX9t6v\nXzR+WG8dz4vDobzfNeUsgO4AnsN1pW/igXG9grsrfqlHRjUgoDwiY0v8+xLAZyXuEaC8eeuhvKFn\ngf9qWjrJcHlNuciyvOvkyZNNbTYbsrKycPF8NjKzslBks8NiokEV276yLP1nDtgdDBpHNkKTqChE\nNo5CdHQ0GjduDIOHAuwq4lxWJhLH9fb6POGhQejUrhm3ff+JPgA8kQguQdkbAEAkgAMAvgfwMq4/\nGTSBN/6Kw6CYHE6l3osb65PLUN6U8wB2QlHqUy5jFEDxFtSUdhYT9TQgD934xwZ0jG+Gzm3jMLJH\nTyTFN0VckwbQ6ZSH3tsLVt7Q69bBcDiSlokDqRnYeyQdv67YgczzuYhq0hiJN3VGly5dQNO0B0S8\nkQsXLuDqtavo3M43tTV7JLay7Dmc1tULFRByALQBMBlKwVxAWfTU8P2XcsNXJ1ZnIq4rdQGU1bsk\nKQDMUBQ6E8AxKHE2vsII4M6QIPM0AG0euXsAOfq2noZmUQ3/U/KyCO3+AK7uWlzhwA6Gw57DaViw\nYhO27DmKxKQk9OjZC40bN/aY8KtWrkCnpiF47bGRHhuzIn77+188OXPhjvwCWy8vTxUK5WnwNpT6\nOqpugN1d8XOgKPYvUGzvE2Xck1IzkapNA8pIPqPTEY+3bR6lnzJucNCQPkkei/0HABNtRN8ubdG3\nS1tcyLuKr1dvwcKFXyCsXj1063EzEhMTK/xwVQbLstj3zz/45FlPPOyqRlJ8UzhYrgO8vxpfBfA4\ngP8D8BcUXVINd1d8CsBTUDaRKzwvTrUgANxjooxfjBrUnZo8ZhAV39T9PVVVVvyyEAQRa7cdwOyv\nfoUAEqPuuRfOIlXusm3bNhReOI0f3ptarddXB1mWETvgCUeBzdEavo270QP4AMB8AKk+nBeA+14d\nFsB7uK70KVA2LmrFyza0muk10Y3CF/46b1rw3OkPVUvpa4LBoMew5E7YvOh/GD0gCR99+AG2bdvm\ndl3MwsJCbPxjA16aONxLkpYNQRDo0CqGB9DZpxMrupcNxXJo7eW5Ulwv1PQY/xyUtjetUEYNci9C\nABhtoox/PHRXv1bfvPOUKaZxgxoNKMsyendqU+3X63QEuie0xJDeSfj6p9+w9599aNqsOczmytMA\nZFnG998tw/A+HTBmiLdN7dKcyDhv2H0o7Siu5yj7AgnKqe8XuO6xGwblqeNpX3KNSwi6cgrAKCh2\nGwAMAtCvhmNWBh1koVc6V/mZU+41eiI4q6RHpya0btoYfy18FffckoiPPvwA+/ZV3vVm//79KLhy\nCdMfvtMjMrhLsMWkN+h1aqWnXSv+Xw9lEU0F0MHbk3oqLNn5XDcAWATPuCLLIshqpjf36dR24N4f\n3rZotZ2OXq/DM/fdjt8/fQnrf/8V27aV3+vs0qVL+GX1SixIeditEAVPQhlJkAaD2l0qRAD3ABgL\nwNlqM8Jbk3k6Hv83KCmEzoOlGHhO+DCrmd45vF/nhCVvP2nyZAiut2jfMhobvpiOXdu2YOOff5T6\neX5+Pr74fD5ef2IUbmqrXi1PmiKh1+usqglwI9uhHF4SUOry/ALA42+ONxJRWFy32QYAOArlEKMm\nWK1mesuYIb1aznt1Iq3X+0/+TFyThvhjwXQc2r8XW7Zs/u96Xl4ePvt0HqaOHYj77+irnoAABFGC\nJMla678rQ2mzugXX3e7VzasodYDl7fP3RVDyaJ2fWCOUjYs7Lg/aaqY3DO3bqeXs58Yb/bGWfGSD\nMKyZNw0DH34TNG1CgwYNsGzpErw8aRgmjbxFbfHAsBwEUdRUSEExLJRYLCeroPTP+hDuBcGluF7w\nRXRmGq7bbA8AeAzA81AOMSrFaqa/6JXUOnHeqxNpf1R6J9GNwvHLJy9gwCQlfXD+a5MwpE+SylIp\nMCwPjhe0qPiuPAfgHQC3AVhZk4F8bTMsgHJkXdVn+yDaSI5c8H+PmfzJvCmPlrGR+OSVh1Av2ILk\nrtopbpWWecEGJXZK66RB8SI6lf5pVNOL6GttkqGksTmD1m6DshEuawMcYqaNyxa8/qg5WEP1ZGrK\nsOROuKltU8z8TCsH38DeI6dFKCHe/kY2gIUAHnT3hWovo3sA2KHE3d+A1Ux/MmJAN0u/bu19L5WX\nmf3ceKz4czd2HXQNSvU9doZFzqV8E5TcBX/jZygRoD8Wf98dZS+iKa4X1Fb8KwCehbLyA4r7834A\nt9FG8q63nxnr+bhfDRAeGoR3n5+AyTO/hJ1R15ly+GQmLCYqAwCnqiDVh4WSlAQomV9HAYxxuWeG\ny/eqK74T5xG1BcDjJsr4S20zcVwZltwJifFxeHtBjfZoNeZAagYEUdqpqhCe4x0A3QD8W/x9PZSj\n41pRfCfHASzqntCSr40mjitvTh2DJb/8jauF6lRGA4A9h9NsRXam/KNl/+M0rofLvwwlC6wUWlN8\nIshimjZl3GC/KfBaExrVD8Ut3drj+993qDK/KErOJtO7VBGg5uihFBp2ZTmUBKjHoZwh7Xe9QWuK\n38NqpiP6dtGOq8/bTBx5Cxau2KhK18ONuw+DE4RsVFw/R03Kih5OgRJteQaKY+SjMu6ZB2A0gMZQ\nyp/c5HqDpsqLBFtMzz85ZpCpJllM/kbPxFYw6PX4e99x9O3s2w/8J8vWFRUUOWb5dNKKGQOlYoYz\ntfUiANf+TdughE+fhRLCXJZ3oNLDUS0pfkNOEAaPG9q77mg9lESQiXf1x8Llf/lU8TPOX8Kew2kE\ngB98NGUPKJGXTqVuDKUWT8mWQTlQwhKc+dplpSf+6QlhNKP4pEE/8Y5+XeTQYLWjY33PPYN74o3P\nlyPnUj4iG4T5ZM6FyzdyOh3xFRRzoabEArgPN5aMGYLroSqA4rlLA7ARilKfRemYrc0ekKVKaEbx\nrWZ6xF23dq29/ssKCLaY0PumNthx4ARG3trd6/MxLIevVm6W7Az3cRVutwIYhxuVOgU3rrwklC42\ne6EcKjlrIZVkb/E/TaAVxSfsDNuupmXz/JnE+FjsP57hE8X/4qeNoo4g9kDJoLsDSuVjp1KvArCk\nxO0EgC5QFNlZ4Ougy5BpAF7xstgeRSuK39REU4gID1VbDtVIatMUc5b85tExD6Rm4NTZHGRduIys\nC3loEROJQb0S8PaClayD5R4qvm0klPS/s1ByYF1jdgoBTPKoYBpAK4rfKbF1rJbLTHudxPg4HDpx\nFpIkVak2z7ncKziVkYPMC3nIunAZBr0e0ybeccM9y9ZsRV5+IaIbhaNt8yh0aBmLh1791MaL4iu4\nXuTrPs//NtpHE4pPGQ3duie00krqmyqEhwYhOMiM9OxcNGoQhvTsXGQVK/XFvKt47fFRN9y/9Je/\nse3fVEQ3Ckd0ZH20io0sNea7z0+44fu5y9aK6dm5qYIgzvXqL+MHaELxzTTVJ6lNXJ1yY5ZFYus4\nHDiRgaJ/UzH/hw2IblS/WLHDIYoSSuYkvDTJvYoMaZkX8NYXK1kHy42G9urX+xxNKD7HC819XQhK\ni8Q3a4yTGTl4+eERuP/OZI+Ny7AcHnhlnquJU6fRxCoriBJtNdfKCGS3sJpoOBjPRgfzgoBxL35s\nzzh36a+AiXMdTSi+KEmkWjVltARNGcFwfOU3VhFJkvDIjC8cuw+n7S2yM6MQMHH+QxuKL0p6yqgJ\nq0tVKKMBrIcUnxcEPPTafMcfOw8dKbIzQ+C/iSZeQRPaptMRkiBIOlIT0qgHL4gweuBNcDAcxk37\n2L7ncNruIjtzOzTej0oNNLHiG3Q6nmEDCxLDcjUuI3jsdDaSH0yx7TmctqHIzgxCQOnLRBNrrF6v\n4xwsb/JNeJZ2cbB8uU2dK0MQRMxZ8pvwwddrWF4QnhVEaQHUabvjF2hC8Y1G8nz6udyQxg3rtuqn\nn8tFz4RWbr/u2OlsPPjqp7bzuVcOOFhuLJToxwAVoAlTh+eF7QdSM9QWQ3UOpmYgwY1AvQKbA7MW\nrhZueej/bKczLzxbaGN6I6D0VUITK77Nwe7YdfDkvZPHDKqzYQs2B4uMc5fQplnlB3nHTmfjsx/+\nYH5avxMkqf/TwXKTEVB4t9CE4gPYt+/YmTptjx4+mYn4Zk3K9epwvIBfN/+DucvWFZ7IOC9IkjSP\n5YT5DtYvSv9pDq0o/rHcKwV0oc2BoFpcS6ciDqRmIDE+9r/vL+RdxYHUDOw/ni7tOniq6J+jp40G\ng/7wtUL7LCg14z130lUH0YriC1YTnXboZGabXkne7gOmTfYcTsO53MvcsMmzmEMnzhoYjpfNtPGw\nzcH+zfHCHijZSwFzxkNoRfHB8fzGP3ceatkrqbVmZPIVgiBi054jXH6BbTaUAkj7AJxlOb5Om3/e\nRBNeHQCwM9yni1Zs4jm+7uWjrNt+AJIknwTwGpRiSBkI+OC9imYUH0r5wKNrNlfeJbC28cmydYXX\niuxaqm9T69GS4uNakX3W3GVrC9WWw5ekZV7AgRMZgFKdIICP0JTiA1h9IuO8eOx0ttpy+IwFP//J\nQekUw6gtS11Ca4rPi6I07/Mf/9BaBz6vYGdYfPPrVolh+U/UlqWuoTXFB8sL839ct1M8nXlBbVG8\nzsffrBX0et3fANLVlqWuUVY1WrUpJAgwOw+euvm+4X38sr1nVThyKguTZy502BzsQAAFastT19Ci\n4kOS5N02BzPSRBsbdmnfQnNPpZrCCwKGPznbdim/YKosY4va8tRFtKpUUqGNGf3GZ8u52mjyvL94\njXAhL3+fJMmL1JalrqLJFb+YKwTgqG0mTwkTpz8CJo5qaFnxFZPHzoyUZLlBr6TWWn06VZlCmwMj\nprxry7saMHHURtOKD0DmeOG3fUdP318/LMiUGB/nt8s+w3IYMeU9+5nsi8s5XijVfjKAb9G64gNA\nAS+Iq7b8c2xCs6iGdJtmUX6n/LwgYOwLH9kPpGb8aXOwYxGob6M6/qD4AHBFEMS1G3YcGhvVsB7V\nvmWM3yi/g+Ew+rk59r1HTu8ssjN34npP3wAq4i+KDwAXBUFc/dfuI2MsJorq3L655m3+ApsDdzw1\n2374ZOa6IjszAoHkEc3gT4oPAHm8IP648+DJu4+dzib7dm5L0pRRbZnKZPehUxg+eZY9++Llb20O\n9j7c2OQsgMr4jcnggtViouYYScO4z2Y8bBrUy7UjpHo4GA4pn/7ELvllC+NguIcArFBbpgCl8VfF\nd9LPbKK+G9I7Kfi9FyaYQoPU7Zi4+9ApPPTqfFtBkX1joZ2ZCCBPVYEClIu/Kz6grP4fUkZyzLvP\njzcP79cZpMG32YsX8q5iztdr2CW//s04GG4ilCyqABqmNii+k+QQq/l9gkD8w3cPMD54Zz+DNyuz\nybKM7ftP4NPv1tv+2n1EZzDovy+yMy8isMr7BbVJ8Z20s5iopwVRGte3S1tp8phBlj6d2sBTIQ8F\nNge+/327PHfZWtvVQtsVm4N9V5LkJQiEH/gVtVHxnQQRBMZZzaZpBIGIhNZxfI+EVtakNnG6xPi4\nKnUQFwQRJzLO40BqBv45cprddegUczrroslEGf+8VmR/F8AWBJLC/ZLarPhOCAAxADqRBn03q5nu\nY2e4DrSRJNo2j+JDgy06E2XUmWijTpIk2c6wkp3hpOwLl+XTWRfNlJG8pNcR/1wttG+BUvZjPwKr\nu99TFxS/LJwfhg5QWtaboLSkF6HUk2cA5CKg5AECBAgQIECAAAECBAjgF/w/p8P6SHvkreIAAAAA\nSUVORK5CYII=\n",
       "prompt_number": 6,
       "text": [
        "<IPython.core.display.Image at 0x45c2510>"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Structure check\n",
      "---------------\n",
      "\n",
      "Let us check the effect the distorsion has on the crystal by calculating the stress tensor and forces acting on atoms in the distorted structure. Note that both forces as well as stresses are far from zero."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"Stress tensor  (Voigt notation, GPa):\\n\", cryst.get_stress()/GPa\n",
      "\n",
      "# Print also the forces (eV/A)\n",
      "\n",
      "print \"\\nForces on atoms (eV/A)\"\n",
      "print \"======================\"\n",
      "print cryst.get_forces()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Stress tensor  (Voigt notation, GPa):\n",
        "[-1.445 -1.858  2.36  -4.608 -3.146 -1.636]"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "\n",
        "Forces on atoms (eV/A)\n",
        "======================\n",
        "[[-0.00182192 -0.00303699 -0.00235983]\n",
        " [-0.00233322 -0.00147804 -0.00041353]\n",
        " [-0.00242831 -0.0023165   0.00020909]\n",
        " [-0.00177738  0.00091272 -0.00011439]\n",
        " [ 0.00153439  0.00353868  0.00321501]\n",
        " [ 0.00208721  0.00029438 -0.0003175 ]\n",
        " [ 0.00309936  0.00178758 -0.00212825]\n",
        " [ 0.00163988  0.00029818  0.00190939]]\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Optimization procedure\n",
      "----------------------\n",
      "\n",
      "Our goal is to find the the structure as close as possible to minimum of energy (i.e. with as small as possible stresses on the unit cell and forces on atoms). The Quantum Espresso has a special mode of working where it minimizes forces and stresses.\n",
      "The QE-util package provides full access to this mode by the way of `calc` parameter of the calculator object. There are two variants of the optimization procedure:\n",
      "\n",
      "* `'relax'` changes just the atomic positions keeping the cell fixed\n",
      "* `'vc-relax'` which changes atomic positions as well as unit cell shape and size\n",
      "\n",
      "We will first try to optimize the internal degrees of freedom (atomic positions) and in the second step optimize both atomic positions and the unit cell parameters. This is actually a recommended procedure in the case of complicated structures. Running just `vc-relax` optimization may not converge at all in such cases or, even worse, give you a *false minimum*.\n",
      "\n",
      "**Note:** The optimization runs take much longer then the simple energy calculations. Be patient."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Switch to the atomic position relaxation mode\n",
      "qe.set(calc='relax')\n",
      "\n",
      "# Switch off the use of symmetries. \n",
      "qe.set(use_symmetry=False)\n",
      "\n",
      "# Force recalculation by clearing the results from the previous calculation.\n",
      "qe.reset()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Run the calculation and get the stresses and forces at the end.\n",
      "# The structure in cryst is *not* modified\n",
      "print \"Stress:\\n\", cryst.get_stress()/GPa\n",
      "print \"\\nForces:\\n\", cryst.get_forces()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Stress:\n",
        "[ 0.276 -0.392  4.008 -0.027 -0.018 -0.012]"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "\n",
        "Forces:\n",
        "[[ -2.02115509e-05  -2.47100075e-05  -1.84702755e-05]\n",
        " [ -3.78708907e-06   2.69222867e-06  -1.12200827e-05]\n",
        " [ -1.30597594e-05  -6.89003624e-06   4.93017778e-06]\n",
        " [ -1.27279953e-05  -4.46656372e-06   2.32351547e-06]\n",
        " [  3.25844458e-05   9.91480688e-06   4.82866497e-06]\n",
        " [  2.21889115e-06   1.71758900e-05  -4.63847430e-06]\n",
        " [ -5.96358599e-06  -1.83975441e-05   1.44568246e-05]\n",
        " [  2.09470326e-05   2.46804482e-05   7.79042766e-06]]\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "After running the first cycle, we note that the stresses have been reduced mostly to the diagonal part of the tensor (first three numbers in the Voigt notation) and forces have been diminished to the level of $10^{5}$ eV/A. We can improve this result by lowering the convergence condition to the $10^{-8}$ eV/A (note the translation of units to the units required by Quantum Espresso). We change the condition and re-run the calculation on the same structure."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "qe.set(forc_conv_thr=1e-8*Rydberg/Bohr)\n",
      "qe.reset()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"Stress:\\n\", cryst.get_stress()/GPa\n",
      "print \"\\nForces:\\n\", cryst.get_forces()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Stress:\n",
        "[ 0.276 -0.392  4.008 -0.    -0.     0.   ]"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "\n",
        "Forces:\n",
        "[[  2.11582259e-07  -5.12620252e-07  -3.57045062e-07]\n",
        " [ -2.83535784e-07  -2.86258350e-07  -9.60676801e-08]\n",
        " [ -3.31764094e-07   4.87728222e-07   8.94557345e-08]\n",
        " [  2.74979149e-07   2.48142429e-07   3.50433116e-07]\n",
        " [  5.40623787e-07   5.09508749e-08  -1.34572540e-07]\n",
        " [ -9.15171058e-07   5.61237499e-07   1.48963245e-07]\n",
        " [ -4.84616718e-07  -4.15385758e-07   6.10632622e-08]\n",
        " [  9.87513521e-07  -1.33794664e-07  -6.18411382e-08]]\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Updating the structure\n",
      "----------------------\n",
      "\n",
      "We can see that the forces dropped by another two orders of magnitude and the stress got concentrated even more on the diagonal elements. Thus, let us update the positions in the `cryst` object with a our new positions. The found positions are present in the `results` dictionary of the calculator object. We can check the symmetry of the resulting crystal - but it will probably be different from the desired F-43m. Which should be expected - the unit cell is still the same deformed unit cell we have generated at the start."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Update the positions using calculated values\n",
      "cryst.set_scaled_positions(qe.results['atomic_positions'])\n",
      "\n",
      "# Check the symmetry. Probably not the F-43m !\n",
      "print \"Symmetry group:\", spglib.get_spacegroup(cryst,symprec=1e-4)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Symmetry group: P1         (1)\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Full cell optimization\n",
      "----------------------\n",
      "\n",
      "We have approched optimal atomic positions in the unit cell. It is time to search for the optimal shape and size of the unit cell itself. We need to switch to the `vc-relax` mode and run the optimization again. Then update the structure again - this time both unit cell and atomic positions and check the symmetry."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "qe.set(calc='vc-relax')\n",
      "qe.reset()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"Stress:\\n\", cryst.get_stress()/GPa\n",
      "print \"\\nForces:\\n\", cryst.get_forces()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Stress:\n",
        "[-0.001 -0.001 -0.001 -0.    -0.     0.   ]"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "\n",
        "Forces:\n",
        "[[  2.37252165e-08  -1.78911469e-08  -1.78911469e-08]\n",
        " [ -2.41141545e-08  -7.77875952e-09   2.06137127e-08]\n",
        " [ -3.50044178e-08   7.77875952e-10   2.95592862e-08]\n",
        " [  1.16681393e-08   1.94468988e-08   3.22818520e-08]\n",
        " [ -1.88634918e-07   1.07735819e-07  -3.01037993e-07]\n",
        " [  1.92913236e-07  -9.17893623e-08  -3.16984450e-07]\n",
        " [  1.73466337e-07   9.02336104e-08   2.69922955e-07]\n",
        " [ -1.54019438e-07  -1.01123874e-07   2.83535784e-07]]\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Update the crystal\n",
      "cryst.set_cell(qe.results['cell'])\n",
      "cryst.set_scaled_positions(qe.results['atomic_positions'])\n",
      "\n",
      "# Check the symmetry\n",
      "print \"Symmetry group:\", spglib.get_spacegroup(cryst,symprec=1e-4)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Symmetry group: F-43m      (216)\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Finishing\n",
      "---------\n",
      "\n",
      "After a single run the forces/stresses as well as the symmetry will not be converged well enough. You will need to repeat the cycle of updating the crystal and optimizing it again (two cells above) one or two more times - until you reach the appropriate level of forces/stresses and symmetry. \n",
      "\n",
      "We can conclude the run by \"snapping\" the structure to the high-symmetry values by rounding the positions and sizes to appropriate number of decimal places. \n",
      "\n",
      "Finally we check the symmetry of the structure and verify the stresses and forces by running a single point (`scf` - Self Consistent Field) calculation on the final structure - which should result in stress components below 0.01 GPa and forces below $10^{-6}$ eV/A as well as recovery of the expected F-43m symmetry group of the crystal."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Round the sizes and positions to get to the high-symmetry structure\n",
      "cryst.set_cell(np.round(qe.results['cell'],4))\n",
      "cryst.set_scaled_positions(np.round(qe.results['atomic_positions'],3))\n",
      "\n",
      "# See the structure\n",
      "print \"Unit cell:\\n\", cryst.get_cell()\n",
      "print \"\\nAtomic positions:\\n\", cryst.get_scaled_positions()\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Unit cell:\n",
        "[[ 4.3342 -0.     -0.    ]\n",
        " [-0.      4.3342  0.    ]\n",
        " [-0.      0.      4.3342]]\n",
        "\n",
        "Atomic positions:\n",
        "[[ 0.004  0.994  0.997]\n",
        " [ 0.004  0.494  0.497]\n",
        " [ 0.504  0.994  0.497]\n",
        " [ 0.504  0.494  0.997]\n",
        " [ 0.254  0.244  0.247]\n",
        " [ 0.754  0.744  0.247]\n",
        " [ 0.754  0.244  0.747]\n",
        " [ 0.254  0.744  0.747]]\n"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Vierify the symmetry\n",
      "print \"Symmetry group:\", spglib.get_spacegroup(cryst)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Symmetry group: F-43m      (216)\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Switch to the single point energy calculation mode\n",
      "qe.set(calc='scf')\n",
      "qe.reset()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Verify the final stresses and forces\n",
      "\n",
      "print \"Stress tensor  (Voigt notation, GPa):\\n\", cryst.get_stress()/GPa\n",
      "\n",
      "# Print also the forces (eV/A)\n",
      "\n",
      "print \"\\nForces on atoms (eV/A)\"\n",
      "print \"======================\"\n",
      "print cryst.get_forces()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Stress tensor  (Voigt notation, GPa):\n",
        "[-0.024 -0.024 -0.024  0.     0.     0.   ]"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "\n",
        "Forces on atoms (eV/A)\n",
        "======================\n",
        "[[  3.50044178e-09  -8.80944516e-07   6.17633506e-07]\n",
        " [  2.48920305e-08   8.67720624e-07  -6.30857397e-07]\n",
        " [  5.79517584e-08  -7.93433471e-07  -5.72127763e-07]\n",
        " [  7.38982154e-09   8.83278143e-07   6.89975969e-07]\n",
        " [ -5.62404313e-07  -1.03457502e-06   1.72299523e-07]\n",
        " [  2.12749073e-07   1.23021082e-06   8.94557345e-08]\n",
        " [  5.25844144e-07  -1.72688461e-06  -2.76534901e-07]\n",
        " [ -2.69534017e-07   1.45501697e-06  -8.98446725e-08]]\n"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Display the structure\n",
      "ase.io.write('crystal.png', cryst, format='png', show_unit_cell=2, rotation='115y,15x', scale=30)\n",
      "Image(filename='crystal.png')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAANcAAAC9CAYAAAAgEuiJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXdYFFf3x7+zZWYbgmAUERDBgiUKgr13jYklRmM05jXF\nmGh68sZoTEL6q+aniSZRY9cYE1NNTGyx94ZgxYJUURFE2Dr998eAIoK03Z1Z2M/z8OjuzM49LHPm\nnnvuKYAXL168ePHixYsXLyDkFqCG4g8gBkCUWqWqS2o1Jo1aZQIBkeV4K8NyZkEQcwGcABAPIF9W\nab24BK9yVR8KQA8VQXTw9TH0YliuPcfzvpFNGtk7PthUX9/fl9RRWugoLQDAQbOw0wyu3bhFHz51\nyXEx7aqBIjU5apX6eL7ZulsEjgHYB4CT85fyUn28ylV1GutI7RSCIF4ID66P7jEtdTGtmpBRkWFo\nGhoIlUpVoYtwHI8LaVeRkJSK42cu03uOn6OvXL/JcTy/gGG5xQCuuvbX8OIqvMpVOVQABvqaDP9l\nOa7r+Id7qCY91o9sHhbk1EFOXUzH4vXbHL9sPUSQWs2/BRb7HAB7AIhOHciLS/EqV8V51Kinvgms\n52d69cmHTKMGdoZRT7l0wHyLDT9tOiAuWLvJmldgzbXYHM8D2OrSQb04Da9ylU89H6NuqY9RP2DR\nB88berSPBEG492sTRRFbDyRi6sfLbA6G/dVic7wEoMCtQnipNF7luj+P6nXk8qeG9dJ9MOUxyqBz\n7UxVHvkWG96Zu9axYcdRi83BjId3FlM0XuUqnXo+Rt3SOkZD/+WfvGjs1LaZ3PLcxY7Dp/H8B4u9\ns5jC8SrXvTQz6Mh9Ex7p6Rs3dQyl15Fyy1Mq+RYbpv3fWsdfu45dt9rpbgCuyC2Tm/AFEA0gFICu\n8EcLgAbgAGADkATgNABGJhkBeJWrJG31FLn789fH1Zk4onfFfOkyM3fVRu6LFX/l2hx0NwDJcsvj\nZOoAaA8gxs/H0JsXhBia4eo1axxoiwgJVBn1lFqvI9WkVkM4aFawOWi+wGIXzlzKEK/l3NIb9VQK\nxwsHLTbHfgDHAZwCwLpLeK9y3SFaT5G7v5n5rOnRAZ086ntZ/tsO4d2vfsy300wnABdLHDYAaAwg\nCNLNle1u+SoJAaCTj1H/JsOwjzRtHOjo3La5LqZ1OBUVGYbmjRtCo1GXexGrncbpi+lISErF4ZMX\nrUdPJws5eWZRhLjQQbPfAkh3xy/iBWih15GHFr0/yXd43w4e+Z2s+H2n8O5X63JtDqY9gEwAjwFY\nCMAEIKPwvU4AlgB4D4BZLlnLwABgXB2jfppeRzZ8adxg/fiHe6j8fU1OG+B8ShaW/PIv/cPf+0St\nRrMv32KbDWA7AMFpgxTDI28kJxNs0JEnZr/5pP+Tj/T0CFOwLL5a8w83e/mGLKudbg8pfEoHaaYq\n2nwWAXwO4CNI6xMC8m9Mmww66hNBFJ7r0q65+NK4waY+HVtXOMKlKljtNH7eclCc//0ma/bNfLPN\nQb8nCOJyOPm7qO3KRfgYdfunjB3UYfqkkRq5hXEGb85ezazfcmCT2eoYUcphEXf+5iYAuwB8AeAn\nyKNkfQ06at1DPaPrvP/iY7rQhvXcOrgoijh88hJe/d8Ka9b1mwlmm2McnGgu1mrl0qhVLzRr3PCL\nvas/MlbEjvcE7A4GsWPesV7JvvkUgN9KHC6uXADQG8BnAMYBSHWLgBImo576iiK1Yxd9MMkwsGs7\nNw59LxzHY97qv7m5qzbSNMu+IQjiEjjhYVOblauxniLP7Fz5gTGySSO5ZXEqhxIvYuQrc/LtNNMU\nQE6xQyWVqzgEgE8BrMC9ThFn0tego9YN7RXtM+etCXo/H6MLh6ocZ5Mz8fTMb502i9WMx3XlIXyM\nun/e+M/DocP6xHr0Oqs0ggMDcCOvgLiQmtWKYbkfix0iIJmCpaGBlIO2HFLkxzUni0XoSO17Pkb9\nt0s/muz3+lMPa3WUsvYQH/Cvg4kjepGiiKBjp5Mnczx/AEBaVa9XK2eummgOlqQc8/B+BAC4CWmW\n6wdgPyTnR3UgDDpybn1/3+c3LZ5haPhA3WpezvXsPnoWT/z3K5vNQY8B8HdVrlFblKsngBGQ9nvC\ndaS27c6VcaqW4TXLHCxJoXmYZ6eZQFQ+WoEAsB5ABwCjIG3CVgXCZNAtCQkMGPv3wulGZ7rWXc3R\n05cw8pUv7FY7PUEUxV8r+/maYBIFAZgOYBGATQDOAmhS4hwVgCwAPxIE8W+fTq0dNV2xAKBzu2Zo\n3TRYA2BkFT4uAhgN4D8ALhe+V+kpx6Aj54YEBozdsmSmRykWAHRo0xT/LJquN+qpNQAGVfbzSreJ\nCADPQLo5/gPgVUgxZKeLneMPaWa6AGAbgGUAUgDwxc5JBXAAwDkfg27prDefrOdut69c+PoYqJ1H\nzkTQDLsYQBzKXnOVRRrumIWbAQyHVIogr7wP6kjtzPoBfq9tXTLTqCTHRWVoEOCHrlEttL/9e/hR\nluN3Q9qQrxBym4VdIQVhNi78uQhgZolzvgaQC+mPnAbgJIAbVRyvd0hgwF8nf//C5O6cLLlgOQ7N\nhrxqu2W2doIU/lSdX1wH4BUAOwEcLefc3n4+xr8P/vCJR6yxymPrgUT8Z8Y3N+0Opikq8GABXDtz\nhQAYDOAhSPsoTwH4pcQ5YwFEQvJMHYb0VL1e4px/Ct8/AWlGslVVoDom/Tf/fXpYq9jWEbVDswCo\nVSpYHTSRmJRqYDk+GsCH1bgcB8nBkVX4ejmAFpDWY8UL6pgMOmrPko8m141uWdJC90wiQgKReS2X\nuJR+rRnDciXv41Kp6k2mB9ARd2acxgBmQzLNingEwERIs0164b9/QL5wmyAdpb184Z/5VB2jXiYR\n5CErOw/Rj71tpxlWD+daKy0A/A+SN21p0ZtGPbVkcPeo8cs+frFGfdEWmwPtR0+zZefmV8iDWFbI\nTwtIToEixdlf4mL+AD7BHVPtCO6tvfdX4Y8iIICxI/t1FGubYgFAUP266Bbdgt9x+HT5J1eO85DW\nw0UKOwGAP6nVjPu//z5V475ok0GHZR+9YHj8zXmrbRUwD0vzFmoAJAL4ElIUtR137/IDUmJeDwBP\nAngXwHe415xTFH51jP16d2itk1sOuejX+UGTjtKecNHli6wRQa8j/2/h+5MMfnU804FRHj1iWmLM\n4K4Gk0G3qLxzS1MuDlKhyxQAXwH4GNJ6yKNhOT4mKjJMbjFkIyoyDDqKdOlsYtRTMUO6RzGDu0e5\nchjZ+eSVsTqDjhoKKTazTO63zzUE0gwGSMGdnZwjmizUZVjOv2looNxyyEbb5qGw2BwRkFLiXYGB\n54Xn358yusaZgyUxGXR4Z9IIg6/JMP1+51VkE5kAcAlSCM0AZwgnA+1bhAXZ1OqasGdeNXyMetT3\nr+MA0MpFQzzesW1TMSzoARddXlmMGdSFYDiuJ4Dgss6pyN0mQnK5Noe0vwFISuYxmxcEENO5XbMa\n/0Qtj5hW4QSk4FynU8ekn/bSuMGeFYJRDUwGHcYO6UZQpHZKWeeUpVyl7YVYcWcv4yFInqL+1RPR\nPfjWMfaMaR2urBBsGejUtpnJqKe6uuDSHXSkNrhfpwddcGnlMnlMf0pFEFMAlHpvlaVcceVc93VI\nEdNnCl/Xh/zRHmWiIohGjer7yy2G7AQ38AdFakOdfV0fo+7NKWMH6Wub2R3ZpBFaRQSrIAWF30N1\nvo1TuNOBYz6ktjftq3E9lyGKor6ohU9thiK1gBQA4ExMNMONmDDMs+uPVJWXxg328fMxvF7aMWd9\nIeMh7dDXKXytqFlMEERd4Y1Vq9FLyYnOVq7oJsH1HQF+Pk6+rGfQq0MrWO10NErRJWcpFw8pPXxX\n4ev5kAqfKMLpIUJUqV1YTchTKDTbnF2IJ6bTg83kLaIvI/6+Jvj6GDgA99Q8r+qaqzw+BeAD4Llq\nXscpqAiCZlhvo0YHzQJSxI3T8PUx9IxtE15rI18AILplEwGleGHLUq4PqjneNQCTAcwpfP04pCKV\nspiLBEE47LSsZcMVQeF3UOWsgtIQBbFDbY58AYDObZuZdJS2c8n33WUrXYdU5fV/bhrvLgiCMBdY\nnPrA9kgKrHYIguDM5uYmO800qA1Z3fcjumUTQk+RPUq+765CmLsgeRJ9C19HQiqIX9HGAQSkdP5Y\nlUrVQafTdSIIwlcURR1BEBwAB8/zKQ6HYx+k3KKTKFZUxWqnD5+6mN5pcPeoWr3wSkxKZfLNNmfG\niUaFhzSwaTUa3/JPrblERYbBaqcjIU1Wt0tjl6Vc1UmoKwseUlUhQFK0BZDMxvvNZlEURb0hCMJw\nlUpFNmrUiAkPDzcGBwerdTodNBoNBEEAy7LIy8vrkpqaOiIlJYXLy8vTUxSVbLfbvwKwlmbYw4cS\nL1hwx5tZKzmYcEEUq15opjRCw4Mb1OoHFiA5NVQEoYI0edxOQylLueJcLM8PAHYACCt8bYQU/UFD\nisgfrdfrp6lUqvCePXtSsbGxan9/fxAEcd+Fc/fu3Q0AwLIskpOTI3fv3v3FxYsX5xIE8deJc6m1\n2hfP8wIupF2l4Fzl0hn1VK1XLgDQajW8g2Hvuj/lrI9+DXcKTz4OKS9sKUmSLwYFBdXt16+fqVWr\nVlCrK1+JQKvVIjIyEpGRkcZbt25h//79j+7etVOTk1eAenVr5+R1Kf0atBo1GJarUP2HCqI36Eiv\ncgGgSA1vtt69h6iU5gPr1Gr1wwRBfDhu3DhtVJTz8oH8/PwwdOhQTUbaZSQkpaJ/l7ZOu7YncSIp\nxRWX1ZJajVe5AGjUahEl0nmU8MW0IknyfMuWLQfFxcVpo6KicPr0aaxYsQI5OSUToKtOo+BQHD/r\nkhvMIzh2+jJjtdPOvqzD5mD48k+r+TAsp0KJysSu2kSuKLFarfbgqFGjgp977jmDySRlLDRv3hzB\nwcH47bfKVGG+P61at8FPmw9BFOVuR+V+WI7Dz1sP8gC+dfKl7XYHc9s7lpZ1A9sOnMTy33bgw29/\nxldr/nHycMqlNOUqyyz8AK5XsPZarXbnhAkTTG3b3m2qkSSJAQMG3FaEnJwcnDp1Cj169IBGUzVL\nNjw8HKwA7ItPQo+YltUW3pP4e3c8IOIcgKnVuIwJ0hZKKKSiRSEADhVYbLeVa+3f+3DsdDJCAgMQ\nEhiAyPCgasntKYiiCJphtSixQS/XmquZVqvd8eSTT96jWMUpKtxJEAQuXbqE+Ph4vPHGG6hKQU+C\nINCxc1cs+nl7rVOuBWs3m/Mttln3OaUBpHLNPdVqdQSAuhqN5gpBENksy+bwPH8Ckpt5Mu5U/EoD\ncOFMcuZt62fGpKpUzfZ8Uq/cgFarKWA5/q5WuHIol4aiqD+GDh3q065dxZqeBQQEYNKkScjPzwdB\nEOA4DpmZmQgLC6vUwB06dMAnH/+DqzfyUBOqwFaEpJQrOHMpQ4CUezcA0qyTCkCv1+un8jzfQRAE\noyAIpMlkQp06dQiTyYSIiIh2er0eVqtVTE1NtWRkZKjsdruGoijBbrfvFEVxJQBL7i0zmW+xwddk\nkO+XlJmEpFRQWk2ircSa1u3KpdFopgUFBTXu0aNHpZ0pvr5SIEBOTg5WrlyJxo0b46mnnqqwu16n\n0yE6uj1WbtiN6c+Vmt/mkdwyW5FxLRcZV3PBchyG9+1w+9ji9f/SDoa1AtgAyZuVo9FomtarV4/o\n3bu3T9OmTREQEHA/a4CAFIQNi8WC9PT06IMHDzY/d+7c52q1+hc9qc48eT4tvLZZA8WJP3uZy7fY\ndpV8v6y78n5N0qpDK41G8/3UqVMNBkPVn3Qmkwldu3aFVqtFw4YNAQA0TVdoPebnVxcr1v6K5x/r\nC0/ozSUIAq7n5uPc5UwcTryI8JAG0BR7mHy08Bc8+94i7Dl2FimZ10EQBHrGSjVobpmtePHDJRzH\n89EajUZUqVQj27dv7zdmzBjT0KFDqeDgYBgMhgqb2SRJ4oEHHkD79u3Jzp07a7RabevLl1N8m4UG\nqjq1vSfjotbw2Xe/ma9k5y1AiY6c7oxSJ3Q6XcIjjzzSplu3bk7dAsjNzcXcuXPRr18/9OzZs1wl\nW7NqJbq2CsLHL49xphhVgmE5ZGXfRPq1XIQH10dwg4Dbx/bFJ2HUa1/Ax6hHSGA9hAQGYN60/6B4\nYiLLcdCo1aUqyOS47xx/7jy2WQDRPCgoqPGECROM/v7OLXdw8OBBWK5dwtpZLzv1up6CKIoI6j2Z\nttNMY5QojOtOs7AzSZIRXbp0cfreWkBAAF555RXs3r27QuePeHQU5n4xB8P7xCC2TYSzxbkLi82B\njGu5UBEEWjS54z3jOB5Ro/6L67n5aFDPDyGBAZj2zPC7lKvjg02RuvVb6HVl19bRlvEg2bo/EX/s\nOMqJIAaNGDGc6tKli8oVnV0aN26MH/fuLP/EGkrqlRsgCMKGUipOu025dDrdW3369NGrXJQR3KBB\nA4wZI81EVqsVq1evxqBBgxAeHn7PuT4+Pnhk+Ag8/+ESHPj+I1S1N68oiriZb8GNvAKUbFr+ymfL\n8deu43DQLIID/fHYwC6Y9uzw28c1GjU2f/cuAgP8yjRPSa2mSiU88wosmDhzIVfH10/z4osv6pw9\nWxUnMDAQNgeD0xcz0KZZiMvGUSp/7DjKq9WqLaUdK+tRFgfn7nM9oNVq0z/88ENdddZaFUUQBMTH\nx2PLli149dVXUbQ5XRxRFLFm9Up0a9WoTPOQ43hczbmF+v51ULwGx89bDmLO8j+ReT0XJKlFy/BG\n2LRoxl2fTc26AaOOQr26PlXaOqgqoigi9vHpPC2oxSlTpmrc8X1v3boF9bR2zJ8+0eVjKQmeF9Dy\nkdes2TcLekNqCHgXZf3VxfscqzRqtfqd6Ojo95988km3FuYUBAEqlQqCIGDnzp3o3LkzjMY7DQJu\n3ryJ+V/Owy9zX7ttHl7OuI6XPluOjKs5uJ6bjwA/E/6Y//ZdJl1Wdh5uma0IDgyA0rqmPDX9Gxw5\nlyG+9tprhE7nnuz7/Px8zJk9C2c2fFGrXPJb9idg0geLzxdY7JGlHXeLWUhR1LCoqCi334VFJijH\nccjLy8Nnn32GiIgI3Lp1C3l5ebDb7YiNjcWE6d9g65J3ERIYgPoBvpj27HCEBtZDowb+kmlWgqD6\ndRFUX3n7ZN/9/C92HT+PadOmuU2xAGmLJLJFc/z4z35MHuOpFc8rz9c/bLEUWOxlbs67Y+YitFpt\nwcyZM01F+1RykZubi40bN6Jnz57w9/eHj48PVCoVdu/ahRNHD2LbknfxgL9npqQknk/FgEmfYfLk\nyYiIcK2TpjQuXryITX/+huPrP3OrGSwXqVey0fmJd60Ohq2PMuqSuCMqPlStVqvlVixA8iqeOHEC\nTZo0ga+v7+2ZrVfv3mj5YBQeeWkOcvIKZJay8mRcy8XDU2bzHTp04ORQLABo2rQpGF7EnuPnZBnf\n3Sz5ZTurUhErcJ+CP5WpFV9VYoKDg1knXs8lDBw0GGHNWmLg85/hSvbN8j+gEC6lX0PX8TNpQq3h\nhw8fLlt+HkEQ6N23P96e+wNqehm78ylZWP7bTtbmYL6433nuSDlpHhwcrPhVLkEQGDzkITzYviMG\nTPoUZy5lyC1SuRw9fQn9nv3IzvLijVGjRpEkKW+viZiYGOiMvpi17E9Z5XAlPC/gmfcWWjmenwYp\neLlM3GEW6nU6nVIynsulT5++6DtgCPo9+xFmL9vAcZzycgEdNIOZ839khk2dXZBvtk1XqVT+bdq0\nkVssEASBRx8bg6W/7UBCUqrc4riEBT9s4tOv5pxhOb7c3DiXK5darTZUpQ6Gqxg0aFC558TGxoLl\nBHb+2k1Huz35nlVJs9ix08noOHaGddWGXdvtNNOMoqj2vXr1IpXyHfv6+mLow8Mw+aOlNc48PJ+S\nhVlLN9Bmq30sipVQKwuXKxfP83aeV87Tf8iQIRU6TxRFlcXm6Hsx/err/Z/92Dp7+QaO5eS7WYpm\nq0emzipIv5rzrNnqGAqA5nl+dJcuXRRlGXTo0KHGmYclzMEK1Ytwx5rLTtO0Rz3CeJ6HIAgqALQg\niEvsNNNqwdrNh2NGv2NZ/vtO0WJzlHsNZ3HLbMXCH7eKUaPevj1bAfgJ0nZJj+DgYKZOHWVtHxSZ\nh8t+34WtBxLL/4DCEUUR7361jknPunGqIuZgEa6qFV+c5KysLKfWJ3c12dnZoCgqG9INDADpZqu9\nR/rVnEc//Obnbc0eesXx+qxVdFLKFZfJkHg+DS9+tMQe+fBrjs+X/v7ntZxbQwpnq+yic1QqVYfw\n8HDjfS4jG76+vpj49DN4Pm4JDiScl1ucajF7+QZ2zV97Ms02x1BUwBwswh3mxLGMjAzFmC2bNm0q\n1zTMyMiAWq0uGSsmAtiWb7FtAxDyw9/7Xvzxn/1TWkUEq54d1denQ5sIRIQ0QFUDkwuLduLoqUtY\nvH6bOeXKDYbj+a8YlvvOQbP3RFwDgE6n6xUaGqqY77YkYWFheGL8BIx7+2v8Ou91xLS+N4ha6Xyz\nbgv/1ZpNOTYH3R13KkZXCHf8YVIYhoHZbIaPj/wN0rZs2VKucqWnpzM2m+1++SsZNMPOABB3/Ozl\nEckZ157heKE9y/G+kU2C7F2imhtiWoVr2zZvjAA/E3QUCV1h4K+dZuCgWdzIK0Di+TQcP5PMHEy8\naL+YdtVAaTW5KrXqaL7ZthTAP7jTg7pUOI5rFxKi7Ej0Fi1aYPTjYzHq9XlY/fkU9PSQjGVRFDF7\n+Qa2ULG64k4X1QrjjlrxIkmSZzMyMmJbtWrlxMu6jsuXL9tRsbLPDID1t8y29YWvAxLPp8WcupAe\n6+tj6MVxfDuW4028IGg5ntcCgEatZtVqFaPVaMxqtSo+v8C6W5Qiqk8wLHerEmL6CoLg9ORHV9C6\ndWuQE57CU9O/wYIZT+OR3ve0slIUHMfj3fnrmO//2pthc9A9UAXFAtxUK56m6Y2JiYltWrVqpfgm\naWazGdnZ2RSAo1X4eC6ArYIobs0rsH5W4hgBQGQ5HizHFzWiqw4mrVbLEgThETXwmzVrhmeeex5T\nP1mMTXtP4PPXxykygv58ShaenvmtNeNazimrnR6KSpqCxXFLxV2e55fEx8cTDof7vGxV5cCBA7xa\nrf4ZgLnckyuHs6uRUmq1usKLayUQGhqKqOhoJGdcQ7fxM7H90Cm5RboNx/GYu2oj13tinO18atbb\nZqujG6qhWID7MpGzNBrN9qNHjw7p0aOHrCHT99tE5nkee/bsoWmanutGkaoKJwiCx4Wfq1QqjOjX\nEZFhQXj58xXo27E1Pnn1CVnz4orNVqftNDMWUum5auO2WvF2u33Ozp07rXKXk76fM+PMmTMQBOEy\ngAT3SVRl7BzHKSMsoxJwLAsdqUWfTm1wYO0nUKlU6DruXSxevw35Fvfu2KRkZmP6vB+Y3k/HWQtn\nq65wkmIB7q0Vv9tqtV49duyYIou1cxyHP//802q32z+SW5YKcpPnecJm86gtRNy8mYOQQKkITx2j\nHl++MxHLPn4Rh09eRNuRb+G1/63EqYvpLhuf5wVs2nsCQ1/83NJl/LuWlX/s+sbuYFryvPAtKrGH\nVRHcWStepGn6iV9++WVv8+bN9UrI7yrOli1bWLPZfBjAL3LLUkF4iqIuZmZmtmnevLncslQIURSR\nlp6Jkg3KO7Vthk5tm+F67i2s3rAHY9/8EsEN/PHMo33Rp2Nr1A+o3r0iCAIuZ2Zjw46j/MKfttIs\ny6cWlvdejxLNE5yJW2poFEer1f4vPDz85RdeeMEgR8ZqaZvIGRkZmD9/voVl2RYAstwuVBXRarUL\nBw8ePLlfv34esfa6ceMGli7+Fuc3zrvveRzHY/P+BKz5cw8On7oEA0WiXWQYoiLDEN0yDFEtwspU\nuCJFSkhKxfGzl9lDiRds55Kv6NVqVYGKILaYbY55cG53zTJx++4+y7Lvp6amjj5y5EiTTp06uf2m\nKLmJzDAMVq1aZWVZdio8SLEAgGXZA5cvX36yX79+95a3UiAZGRmIimxS7nkajRoP94rBw71ipNku\n6wYSklKRkJSKb9dtRcL5VDAMe3tzXiN1zITZaoeDYQUdqb1JajUn8s22XYIoHgMQD8B5zd4qiByh\nMwxN0yN++eWXA76+vqbIyFIL57gFjuOwZMkSm9ls3gxgjWyCVJ39ycnJao7jqtxayZ1cunAeg2Kb\nVuozBEEgrFF9hDWqjxH9OgKQzEurnQbNsHAwLFiWA0VqcfjkRbz86fJUs83R1GqnZV/buyPNvzRO\nsSw7ePny5dYLFy64eKjS4TgOK1eutKenpx+kafoJOH8fyh1cJgji9MmTJ+WWo1xsNhsST57EuIe6\nV/taBEHAZNAhwM8Hjer7I6xRfTR8oC6G9YmF0aBrAKBz9SWuPmW5cne5YewMnuf3JSQkjGnQoAEZ\nGBjohiGBzZs3o0+fPvjuu+9sKSkpe2maHgYpjMkj4Tju5vXr1x/y8fEhz549iwsXLiApKQlJSedw\n6VIyUlNTkZ6eDqvVCoqi4M6Sa8XZv38/mjYwYdzQbi4bgyAIiKKoOXbmcl2G5daX/wnXooSFcAxJ\nkpujo6NNI0eO1Ln6j7927VokJSXZGIb5mabp51BOcKwCCQLQldRqOhr1VC+bg2mjI7WG9q2bokVY\nIIx6EnqKBKnVgGY42GkGFhuNs5ezcOpCKlRqNUJDQhAYFIyQkBA0adLkrkKprkAQBMyZ9TlWfPQ8\nOrdzTTeUz5f8jumTRuJmvgUtH3mdphk2BMANlwxWQZSgXADgS1HU11qt9tEJEyYYWrRo4fQBaJrG\nH3/8QR8/ftzKMMwESFHnnoIKwIA6Jv1/OY7v1qFNU6Zzu2am6JZNVFGRYQis51ehi4iiiIxruUhI\nSkX8uRQcO5OKhKQUtG7VEp27dkdYWJhLag4eP34cxw/uxv41H7qspqFf54m4dWglAGDMG3MtWw+c\nnApgtUsGqyD3C9yNc58YyKdpegJN02uWLVu2tmnTpoa+ffsaIiIiqv3HsFqtOHLkiPDvv/86OI77\ni2GYFwBkmhfGAAAX/0lEQVRUJvpcTgLUKtXTOkr7ZmA9P+MrTz5kemxgZ8Kop6o0vRMEgdCG9RDa\nsB6G9YkFAOTlW7D273347ud1INRadOrSDTExMaAoyim/QEFBAf7a8Ad++/J1txUL7dY+0rg3PqmL\n3cHIqlxu3+eqAD4EQUwkSfK/RqOxbp8+fYwxMTFEZRoKCIKAzMxM7Nmzx56YmEio1epNDodjNoBD\nrhPbqfiZDLp5HMePHdIjWpjyxEBDbOvqP2juhyAI2HX0LBat/xcHEy+iS9euGDBgILTaqgfdi6KI\n1StXoGfbEHw4dbQTpb2X4jPX7mNnMXHGN6fzCqwPunTQclCichVBAOil0+neZlm2n9FoZEJDQ9Gk\nSRNjo0aNCJ1OB61WC0EQwLIs8vLykJaWxqakpNiuXr2qV6vVuRzHfc3z/FIUS433AB4y6Mg1owd1\nMb73wiiqXl3318fIuJaLd+atw4kLGRjz+BOV7j1dxPHjx3Fwz3bsX/PhXV1iXEFx5bpltqLZkFcY\nluMNAGSrjqRk5SqOGkALALEkSXbRaDQdAZhEUdQB4AiCcADItNvtu8Q7m4a5ZVwrDu41eSuKn8mg\nW2jQUcOWfvyCQQkZu39sP4I35nyP9jGxGDR4SKVmseTkZKxetQIb5r91T7iTKyhyaBTRYuirluu5\n+Z0hNVqXBU9RLmeixN/t9mz16atPUCaDcnJKb9wswGuzVldqFktLS8OKZUuw4uMX0KeTPMVKn3jr\nS8umfQmyOjXc3XBcCcTB9ZvkFYXQU+RHdesY562Z9XKdF8YM0JTWskhOjHoKowZ0RKMHfDHrm5XQ\nG4wIDg4u8/wLFy5gzaoVWPzBc+jfpa0bJb2bjKs52sMnL2ZwvLBZLhnckubvpVQIg56a37Ce39P/\nLJpuaBBQMXe6XIzo1xGtm4Zg2EtzQDsc6Nmr113HOY7Dv9u24vChg/j+fy+he3v5wtoAoK6viSBJ\nTT0HI18PELclSyoIJcxahMmgWxEeXP/p7cveNypdsYpo1rghtnw3A8eP7Mf2f7fdfj8jIwNffTkX\nnPk6Dq79WBbF+nzJ73e9pkgtVCqVrEU6lGWDuIc4mccnDDrqy9CG9R7btGiG0UdhbV/LI7RhPWxZ\nPAP9n/sULMvC7rDjZEICPn9tLB4f3FW2xnezlm24y6Ghp7RQEYSsBVNr48wlKzpSO7NBgO+zf3/7\njscpVhGB9fyw8du3ceTQAei4Ahxc+zHGDummqI6SgiBCFEVZmxS4M83fC9CDorTTNy2eYazr6xEp\nWGUSFvQA/l44HcfOXIbNTsstzj3YaQbifbo+ugN31Ir3ImEw6qkfv5n5rL6isYBKp2V4I7w58WFM\n/XQZBEFZVd5ohoXAC1Y5ZaiNZmGcHIMa9dQX/Ts/WPfhXsquNltZJo8eAFEQ8d3P22WVY9qzw+96\nfTXnlmCnGVkzy2ujcskxK/fQaNQT570z0TMXWfdBrVbhm5nPYs6KP3E5o9R+EW6huDMDAA4lXLCw\nHF+VqslOQ65M5NrEbXPQ38PXWWURERqoOPMw8UKaFm4qRFMWXoeGi9GR2hl9O7WpceZgSSaPHgCB\nF/DjpgNyi4KrN/LAMJyAchqCu5raaBa6E5IgiJdmTh5V48zBkqjVKrz19DAsXr8NclRVLr6JnJCU\nCr2OPAWZ66LURuVyp8k7snXTYFWLJkFuHFI++nVug1tmG+LPVqhlsFOZtWzD7f/Hn0sRLDbHHrcL\nUYLaqFxx7hrI12SY9vL4IfJ3/HMTKpUKzz7aB0t/lddzuOvIGQvL8bInxnrXXK6jDUGgxUM9o+WW\nw62Mf7gH/tl7AjfzLbKMn5KZjVMX01UAZIuGL8K7iewijHrq1Umj+5NaDyjW6UwC/HwwpEc01m7c\nK8v4S375lyEIYhkAuywCFMNrFroGDcvx458e0ad2aVYhT4/o7XblmvbscNgdDFZt2C3YHcwCtw5e\nBrVRudwxK7cM8PPhgurXdcNQyiOmVTjSr+agwOq+yWP6pJH4ffsRaNTqYwCS3TbwffBuIruG2JhW\nTZQTIu5mNBo1WjcNQWJSqlvHXbB2kznfYpvt1kHvg9eh4QKMeqpr53bNa2Y4RgWJigxDghuVa+/x\nc0jLynFAQcVea6NZ6HK0Gk0Pd1Q8UjJRkWFIPO+eAAmrncakDxbbbA76WchYSq0ktVG5XG3yaix2\nR9O2LRq7eBhl486Z670FP9I5t8xpAP5yy4AVpDYqV5yLr9+ynp+PSs7u9EqgRVgQsrJvutypsS8+\nCev+2W/jOF7+Qo8l8K65nM9Djer711pnRhEajRoN6vkhOzffZWNY7TSee3+Rze5g/uOyQaqBdxPZ\nyWg0mhEmo6xFhxSDniLhoF1X2mzmV+toi82xCQozB4uolZucLqal0eCcDiGejo7SgnZR3cAFazdx\n67cczLba6ecL31Lc9pF3zeVc6gqCYNKovc8sANCo1eB45zvvVm/YLXy+5Pc8q53uBuBm4dtxTh+o\nmtRG5XKlydve19eXdtXT2tOw04zTu5us3rBbmDZvbb7NwXQDkOHUizuZsh6xiptiPYSYhg0bkg7G\nY1ssOxWaZqGjnKdcC9Zu4j5f8nue3cF0B3DRaRd2Ed5a8U6EJMnG/v7+moKbV+UWRRGYbXYY9dXv\n2GK103hv/o/0T5sPZHvCjFVEbVwcVGRWLuoHFqPVaruQJBkOwFD4vl0QhHy73X4IUgGUeBS2gVWp\nVCZ/f3+ciD8GQRCgUtVGq1si95YZVhuNRtUMXj6QcB7PzFxos9gcm612ehLurLFKEgeFTQq1Ubni\nynhfB+Axg8HwKk3TDxoMBjY0NBTh4eHGunXrEiRJgiAIcBwHu92OzMzMhy9fvmzPzs42qNXqXFEU\nfxUEoa5Op4PRYMDlzGw0DQ1046+lLBKSUtG2ReMqP2CsdhoffP0TvXbjPrudZiYC2FDORz6AhyhX\nHBQmqAsJ02q1LwF4PiQkBL169fJp1qwZDIZy/ekkAFIQBFy/fj3w6NGjz+/fv1+1d+9e1PXzQ0JS\naq1XrqrGV+6LT8Jz7y+yWWyOzXaaud9spWhqY2fJIkitVvsRQRCvdOrUSd2jRw+yfv361bogx3FI\nSEjAxo0b8fjADpj1xngniep5THhnAYb37YDHBnau0PkOmsEfO45i/vebzKlXsmmbg3kO5c9WxVHc\nPVsbzUIAiKEoan2TJk0Cn3jiCb2vr69TLqrRaBAbGwuj0Yhj++Ut0iI3CUmp+GDK6HLPS72SjaW/\n7mBW/rGLV6tU8YX5WP8A4FwupIupbcpFEASxkyTJjqNHj9bFxMQQrmh7ExISgu8vZ4LnBajVtc+p\nkZNXgHyzDeHB91oC13JuISEpFSfOpQg7j5yxnLyQpiIIYnlhav6lagyruO2j2qRcKpIklzIM02vG\njBlw1mxVGiaTCYEN6uPfQycxqFuUy8ZRKr9uO4xWEcH4afNB0AyLqzfyhEOJFy2J51M1dpoVDTry\njNVO72ZY7hCATXBOMZk4J1zDqdSWTWSCJMlVgYGBI9PT012qWEV07NwNi9Zvr3XKJYoiFqzdzN8y\nW+Onf/nDFUEQLHYHk8Vy/BFIWxdpNMPKWgnXXdSKlBOSJOfVq1dv5JQpU9zWxjM6OhrxZ1OQeiXb\nXUMqgn3xScg3W7MsNkenWwXWkQUW+wSW46cB+BVAKmQuMe1OasOCYChFUZOmTp1q1OmqHy1QUUiS\nREyHDlj22063jakEvl23xWp10HPgfiWKc/N45VLTlasuSZKrn3rqKYPRKE1agwYNctvgXbp0xfcb\n98FB145Yw6s38rDjyBmVIIirZBhecTmINVq5KIpaFBsba2jWrNnt94YMGeK28R944AE0atQIv2w7\n7LYx5WT5bzt5rUb9E4ACuWVRAjV5zTWQoqihw4cPd58tWAp9+w9E3Le/yFY73V0kp1/D1+s20xab\n4xO5ZVEKNTbNX6/Xfzhs2DAjRcmbFRweHo4H20bhzTnfyyqHK+F5Ac+8t9DKcfxMKKTarRKoqWZh\nc1EU20VF3esG37Rpk9uFGfLQUBw8lYyNu2XtIuoyFq3fxl/OvH6e5fivZBRDcdtHNVK5SJJ8pWvX\nrhpNKR1GtmzZIoc8GD1mLF6btbrGmYfJ6dfwyeJfabPV8TgAORsix8k4dqnUxFrxelEUJ3bv3t25\n+eXVJCIiosaZhyXMweqELtVIaqJDI9bf35/39/eXW457GPLQUBw+nYKFP22TW5RqI4oiZnz5A305\n8/pZmc1BxVITzcKYJk2akHILURokSeK55yfj/1b9gx/+3ie3ONVi1rIN7Nq/92WarY5BkNccVCxq\nuQVwNnq9/u3OnTs/GBwcXOpxURRRfN/L3ej1erSIbInZ36xEaAN/tIwoXU4ls2DtJv7/Vm3Mttrp\nzgBy5JankDgAu2SW4S5qnHJpNJov+vfvX6dOnTqlHpdTsYowmUxo1qw5Zn29En4+hipn7LobURTx\n6Xe/sfNW/33dZqe7Argit0zF2AmF+Qpq2pqLoGk6KDBQ+en1jRo1wgtTpuKTJRvwxYo/IQjKtqwY\nlsObs1fTi37almqz0+0BpMstk9KpaWn+lEqlss6dO9djZuS8vDws+vZrBNWrgxWfTkFY0ANyi3QP\nJy+k4el3v7Vm3yw4bLbaR6Gw2pXCUNw9W9McGpRKpVL2FFCCunXrIqZDJ5y6mL6vy7h32UXrt4lK\nmcUYlsMni39lB0361HI5M3uK2WrvD2UqFqAwkxCoecrFCoJw399JjgiN8hAEATwv7LQ7mLafLPr1\n9MBJn1pTs27IKtPJC2noMu5d6+L12/bZaaaFKIqroexcrDi5BShJTVMuhyAIqvs9+eWI0CgPhmF4\nURQdAJIsNkf0yQtpH3cZ9679nblrmeT0a26VJfF8Gl74cInjzmzl6Acgy61C1BBqWpq/SFHUjRs3\nbtRv0KCB3LJUmCtXrthwJ+CVZ1huFlj8tGrD7pdXbtg1KTqyCV4aP9hnUNd20Gicv5x00Aw27DiG\nr77/x5J6JVvLcvwiluM/B1C70qidTI2rFa9Wq+MzMjIGe5JyZWZmqiDVlyhOqp1m3gTw7sHEC6PP\nJmdO06hVTSaPGUAN7NZO3TK8UbU6iNgcNE5dSMfG3ceZFb/vEoqVNfsbUlkzAsD3AH4G8CeUbRIq\nkppU/UkDIMhms6WlpaUxsbGxpUZpuDMTuSJYLBbQNK1C2akaDgBr8i22NQDaz1+76eWv123uY3Mw\njcKCHrB1fLCptkObCH10yyYID2kAPaWFtjBgWRRFsBwPm53G+dQsJCSl4sipS7ajp5O5rOw8vdFA\npbEst8kmlTUrrWvIWgBfALAAUHohxjgobFJQlOuyHPQAQiB5q4qbK1EA/gDQEFK0wK3g4OCQt956\ny8f9Ilaec+fOYc2aNfE2my2mkh81AGgHIMbXZOgBoJPNQTfkeF4LENCoVTzH82oChKhWq1ijnkrn\neeGg2ebYD2mWPA2ArsA4GgA8pJnrSQB7AaRVUlZ3oDhXvFJqxRMA/CDdMFdKvL8XQLPC4xkAZgBY\nX+yc8wD6FR5jAPheu3btmsVigclkcoPo1ePUqVM0wzBVcWHaABwEcDDfYvv67kOiluV4EgAtQuQE\njscts62qIhavfBsG4CtISqY8t6vCcFcmsgpAaYXYPwRwCkA+pKfhzBLHRQCvAIiGNHM1xd2KBUgF\nJZMhKRYA5Gs0mg2HDx9WxmbRfaBpGkePHgXHcYudfGkWgBXOLwn9CYC2AA4Uvm4BQFGpPUrCWa54\nEvfOgoMh2emXID1lj+Heafs3ABMANAbgC+DFUq4dD8kVXGFlcTgc/7dr1y67UjZjy+L48eOiRqPZ\nCw9p5lbIFUgPQwCYBuAMgJ7yiaNcKqpcJc8LBLAO0hPsCqQFb2yJcy4AmAVgKIC6AEJxr8cpEUAC\ngLxSjlWHoyzLZiYlJd1zQCmbyKIoYseOHVa73T5bblmqwTMAXsYdZZNzFlPc9tH9lOsPACcg9UZ6\nv8QxMySX7TQAXSCtlQ6VOOcygK2Q1kTOqAVeKRwOx8d//vmnlS/RTV4pm8iJiYkwm805UL4Xrjy2\nQHpIAsASAD9AWpu5mzgZxrwvpSmXBpLXLQqSjd0C9z4VrJD2QPZCio5WYruXH/Ly8hK3b9+uONnM\nZjN+/PFHO03TT6BmJRq+BMli6Su3IEqgNOXiIDkfXoU089wA4LYa605EpGl63LZt2+irV5XVAPyn\nn36yCYLwHe6d7T0dC6QZZHnh61cgmY2KzAx3NWWZhSKkrn5mSE6IPQCWQdpL8iTSeJ5/c9WqVVaO\nU8YEFh8fj4sXL+YyDPOO3LK4gV2Q1tzuKBQa54YxKkVFA9V+AhADIAWAeyNJq4koivEsy3bPyMgI\njoqK0gLyZSMnJydj5cqVNoZhBsOzPIRV5Tqk5cMeSBvRnQEEwTUZzIrLRK7qjvZKAPshTf/8/U9V\nBDqKona0bt06evz48Tq12v25lCkpKVi4cKGNYZhh8HwnRlUZBuBbSPdPyT3N6qK4CI2qCtMBUszZ\nQgA/Ok8cl2KkKGpzRERE+6efftqg1brPa5yUlITly5fbGIYZBWCz2wZWJkYAwZC8yPUgrfGdkYBZ\nY5Sr6LMEJG/Xk5AiLRLv+wn5oSiKWmcwGAZOnDjR2LhxY5cOxjAMNm7cSB88eNDBsuzDADy7nprz\neQrAHADvAfiumtdSnHJV1z4q2vjtCGmqvwFlKxjP83xLh8Ox9Pjx44PsdjsRERGhdoWZmJKSgq+/\n/tqampq6lWGY/gDu3dH2kghpJveDFEhcnYghAgorreZMTfeFtEeWC6AlgExI3kalUfSEq09R1Aq9\nXt9r4MCBhpiYGMIZHVEyMjKwa9cu+8mTJxmWZZ+BFOLlpWKMhbRX9iYAj29q5qpp9GMAzwGYAuB3\nF41RVUqaD/31ev00nue7d+zYkejevTtV2dJsDMMgISEBO3bsMN+8eZPmeX4ez/OLIEW3eKk4akim\nYncAz8osS7VxpY0aC2maPwJAByl3SAnZrGXZ5qFarXYKgMlqtZoKDg5mIyIiTCEhISp/f39oNBqo\nVCpwHAe73Y6srCykpqbaU1JS2Js3bxooijpkt9tnQUrF8AQPqicQDWA8gE8hxZ96FO5aAM6B9EW9\nBSlQV07KW/gSkJIyY9VqdUeKonqKohgkiqJOFEUVQRCMSqWyCoIQ73A49kJKPDwJGeInawENIO1d\ndYD0sFbCw1lxaCGZiFPlFgQK3Mn3Ui6Gwn/rAHgUpT8c49wmjcJ5DtK6zCNS8b0ohhaQLJ/duDd/\n0DurFRIKYA2AFXIL4sXjUENKcwKkGayoTYxXuUpQ9PSJhhTgqahNQC+KpymkrZ8V8CpXmXQHcA7e\nmczL3agANAJQcgNyDKTyAmZIWdDpUFjQrtLQQKryBEg79iEuGifORdf1Unko3JsrGAZgB6SiQzSk\nLIy2Jc5pVPier4vlq5EMhjTdO7sKFeA1H9xJaUmSSyGVhMuCpDyvljhuBNAf0oNW51LpajEhAMYV\n/p+A86oDe5XLOZS2Pn4GUpJtUdGhraWc8zCkZUAIamBnU0+kH4CzkP4w1XV6eJWrasRAKkBTVHTI\ninvXQj0BjATQHkAAvA4qj4CA5E1cher/weKqLU3NwwdSiFFR0aE0AANLnBMJ4AUAQwC0AqD8csZe\nqkRTSJ5FVzk9ahoTAcwF8Cuk4qwlS+XpCt/7D4DeAJrAW0m32nhql5PrkFJatgNoDal8c22lE4A+\nkKoWN4ZU3fixEufUBXAVUrWpNEhVkIvjAPCRa8WsfXi6TUxB8jjpIDk/VkOZNRSrSiAkM6xxsZ/R\nkDypRTwGScHSCn+SIa1PvciMpytXEUGQFCsIUla0RV5xKoQadxQntPDfJbi7eM2DkDIJ0or97IM0\n03hRODVFuQDpd4mCVIIbACJQekO5OLjHqdEGUkZ2kfKcB/BNseMEpI6NWbijOHtQO0qu1QpqknIV\n5wFI9Rm2Q6r4Wry6kDMKmdSD5DErMtVMAEoW+ZxXeKxIcY5CKkfnpZbgqQ6N8rgBKT1hMqT9GEAy\nwyqSIUxCmvWKTLXGABbg7mKowyBtmKZBimsr6SAAgNerIriXmkNNnblKYoTU5+srSKZZe9xRnCOQ\nwnGKaA2pw0vxdc538LBKw168uJN2ADZCKhqTCGm9swBADzmF8uLFixcvXrx48eLFixcvXrxUnv8H\nGvA9psumzwAAAAAASUVORK5CYII=\n",
       "prompt_number": 23,
       "text": [
        "<IPython.core.display.Image at 0x4a742d0>"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}