{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Calculating elastic constants\n",
      "=============================\n",
      "\n",
      "*Read previous tutorials to understand all material presented here.*\n",
      "\n",
      "The initial setup is identical to the one used in \"Remote calculation\". \n",
      "Scroll down below first cell for the material of the tutorial."
     ]
    },
    {
     "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",
      "\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 qe-util package\n",
      "from qeutil import RemoteQE\n",
      "\n",
      "# Access info\n",
      "import host"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Elastic library\n",
      "---------------\n",
      "Calculation of elastic constans requires inclusion of additional library: `elastic`. This library extends standard crystal object, and provides new features related to elastic properties of the crystal. This library must be first installed in the system. See the `Install` document for the information how to do it. Here we just import the library. It **must** be imported before the creation of the first crystal.\n",
      "\n",
      "Calculator definition\n",
      "---------------------\n",
      "Since the calculation of elastic constants requires deformation of the crystal the calculator *must not* try to use symmetry for the calculations. Otherwise each calculation will be done with different internal setup and the results will be inconsistent. Thus we pass the `use_symmetry=False` parameter to the calculator constructor. The calculation of elastic constants requires a large number of individual calculations to be performed, thus only `RemoteQE` calculator is suitable for this type of calculations."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Import additional library for elastic constants calculations\n",
      "import elastic"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Stup the SiC crystal\n",
      "# Create a cubic crystal with a spacegroup F-43m (216)\n",
      "a=4.3366\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])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "qe=RemoteQE(label='SiC-elast',\n",
      "               kpts=[2,2,2],\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,\n",
      "               pseudo_dir='../pspot',\n",
      "               use_symmetry=False,\n",
      "               procs=8)         # Use 8 cores for the calculation\n",
      "\n",
      "print \"Local working directory:\", qe.directory"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Local working directory: calc/SiC-elast.NTcCnS\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Assign the calculator to our system\n",
      "cryst.set_calculator(qe)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Run the calculation to get stress tensor (in Voigt notation, GPa) and pressure (in GPa)\n",
      "print \"Stress tensor     (GPa):\", cryst.get_stress()/ase.units.GPa\n",
      "print \"External pressure (GPa):\", cryst.get_pressure()/ase.units.GPa"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Stress tensor     (GPa): "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.149  0.149  0.149 -0.    -0.    -0.   ]\n",
        "External pressure (GPa): -0.149\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Equation of state\n",
      "-----------------\n",
      "\n",
      "The `elastic` library provides an equation of state calculation procedure, similar to the one we have used in the first tutorial. It automatically runs a number of single point calculations for given range of volumes (from `lo` to `hi` - the values are scaling factors). The default range [0.98,1.02] is equivalent to +/-2% compression/expansion. The data is collected and the Birch-Murnaghan equation of state:\n",
      "\n",
      "$$\n",
      "P(V)=\\frac{B_0}{B'_0} \\left[\\left(\\frac{V_0}{V}\\right)^{B'_0}-1\\right]\n",
      "$$\n",
      "\n",
      "is fitted to the data points. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fit=cryst.get_BM_EOS(lo=0.96,    # lower bound of the volumes range\n",
      "                     hi=1.04,    # higher bound of the volumes range\n",
      "                     n=5)        # number of volume points used in the fit\n",
      "\n",
      "print \"\\nA0=%.4f A  ;  B0=%.1f GPa  ;  B0'=%.2f \" % (fit[0]**(1.0/3), fit[1]/ase.units.GPa, fit[2])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Launching:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 2"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 3"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 4"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 5"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "     Done:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 2"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 3"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 4"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 5"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "\n",
        "A0=4.3358 A  ;  B0=224.8 GPa  ;  B0'=3.74 \n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Quality test\n",
      "------------\n",
      "\n",
      "We can easily check the quality of the fit by plotting the fitted function together with the calculation data. The calculated data points are saved into the `cryst.pv` attribute of the crystal object. The array contains following data: volume, pressure, lattice constants: A,B,C."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get the data from the crystal object\n",
      "pv=cryst.pv.T\n",
      "\n",
      "# Definition of the B-M EOS function\n",
      "func = lambda v, v0, b0, b0p: (b0/b0p)*(pow(v/v0,-b0p) - 1)\n",
      "\n",
      "figsize(12,7)\n",
      "plot(pv[0],pv[1],'+',label='Calculated points',markersize=10,markeredgewidth=2)\n",
      "\n",
      "v=linspace(min(pv[0])*0.99,max(pv[0])*1.01, 50)\n",
      "plot(v, func(v,fit[0],fit[1],fit[2]),label=\"B-M fit\");\n",
      "xlabel('Unit cell volume ($A^3$)')\n",
      "ylabel('Pressure (GPa)')\n",
      "legend();"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAuEAAAG+CAYAAAA9cNotAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdclvXi//H3dTMciLgQEVBUHDjTBFNDcWLuUYk5E8cx\n9WSZp6xMMq302LDMMnOmkDYUc+AgsfJoaFpaWGmmIe6Jioub+/eHP/lCDpRxX4zX8/Hgce5xXff9\nvj338bz98Pl8LsNms9kEAAAAwG4sZgcAAAAAChtKOAAAAGBnlHAAAADAzijhAAAAgJ1RwgEAAAA7\no4QDAAAAdmZqCY+OjlatWrVUvXp1TZ069ZbnlyxZogYNGqh+/fpq3ry5du/enfacr6+v6tevr4YN\nGyowMNCesQEAAIBsMczaJ9xqtapmzZrauHGjvLy8FBAQoMjISPn7+6cds3XrVtWuXVtubm6Kjo5W\neHi4tm3bJkmqUqWKfvzxR5UpU8aM+AAAAECWmTYSHhcXJz8/P/n6+srJyUmhoaGKiorKcEzTpk3l\n5uYmSWrSpIkOHz6c4XmuMwQAAID8yLQSnpiYKB8fn7T73t7eSkxMvOPxc+fOVceOHdPuG4ahtm3b\nqnHjxpozZ06uZgUAAABykqNZb2wYxj0fu2nTJs2bN09btmxJe2zLli3y9PTUyZMn1a5dO9WqVUtB\nQUG5ERUAAADIUaaVcC8vLyUkJKTdT0hIkLe39y3H7d69W0OHDlV0dLRKly6d9rinp6ckyd3dXT16\n9FBcXNwtJdzPz09//vlnLn0CAAAA4IZq1app//7993y8aSW8cePG2rdvnw4ePKiKFStq6dKlioyM\nzHDM33//rZ49e2rx4sXy8/NLezw5OVlWq1Wurq66dOmS1q9fr4kTJ97yHn/++SfzxnGL8PBwhYeH\nmx0DeQzfC9wO3wvcDt8L3M79zPKQTCzhjo6OmjlzpkJCQmS1WhUWFiZ/f3/Nnj1bkjR8+HBNmjRJ\nZ8+e1YgRIyRJTk5OiouL07Fjx9SzZ09JUkpKivr27av27dub9VEAAACA+2JaCZekRx55RI888kiG\nx4YPH552+5NPPtEnn3xyy3lVq1bVTz/9lOv5AAAAgNzAFTNR6AQHB5sdAXkQ3wvcDt8L3A7fC+QE\n0y7WYw+GYTAnHAAAALnufnunqdNRAAAAcluZMmV09uxZs2OggChdurTOnDmT7ddhJBwAABRo9AHk\npDt9n+73e8accAAAAMDOKOEmOHflnNkRAAAAYCJKuJ2du3JOdWbV0Y4jO8yOAgAAAJNQwu2sVNFS\n+qDjB+r2WTcdOnfI7DgAACCfs1gsOnDgQLZeIzw8XP3798+hRHcXGxsrHx+fHH/dv//+W66urvlm\n/j8l3ATda3XXuGbj1Cmik85fOW92HAAACjXDuPFj1vmSFBERocaNG8vV1VUVK1ZUx44dtWXLluy9\n6H24n0uuDxo0SBMmTMjFNFlTqVIlXbhw4Z4+y8GDB2WxWJSammqHZLdHCTfJ002eVivfVuq1rJeu\nWa+ZHQcAAJjk7bff1jPPPKOXX35ZJ06cUEJCgkaOHKmVK1faLUN+GT3OaWZ+bkq4SQzD0Lsd3lUx\np2L616p/FdovPwAAhdn58+c1ceJEzZo1S927d1exYsXk4OCgTp06aerUqZKkuLg4NW3aVKVLl1bF\nihU1evRoXb9+/bavd/nyZY0dO1a+vr4qVaqUgoKCdOXKldtOAfH19dU333xz29d57LHH5OnpqVKl\nSqlly5aKj4+XJH388ceKiIjQtGnT5Orqqm7dukmSjhw5ol69eql8+fKqWrWq3n///QyZBg0apDJl\nyqhOnTravn37Xf9MLBaL3n//fVWrVk3u7u76z3/+k9aTbDabJk+eLF9fX3l4eGjgwIFKSkqSdOvo\ndnBwsF555RU9/PDDKlmypEJCQnT69GlJUosWLSRJpUqVkqurq3744Qft379fLVu2VKlSpeTu7q7Q\n0NC75swuSriJHCwOiuwVqZ+P/6zXv3vd7DgAAMDOtm7dqitXrqhHjx53PMbR0VEzZszQ6dOntXXr\nVsXExGjWrFm3Pfa5557Trl27tHXrVp05c0b//e9/ZbHcvu7dbdpGp06dtH//fp08eVKNGjVS3759\nJUnDhg1T37599fzzz+vChQuKiopSamqqunTpooYNG+rIkSOKiYnRu+++q/Xr10uSXn31Vf311186\ncOCA1q1bp4ULF2Y6ZWTFihX68ccftXPnTkVFRWnevHmSpPnz52vhwoWKjY3VgQMHdPHiRY0aNeqO\nrxMZGakFCxboxIkTunbtmqZPny5J+u677yTd+EfQhQsX1KRJE02YMEEdOnTQuXPnlJiYqH//+993\nzZhdlHCTlXAuoVV9VunjnR8rYk+E2XEAACjwbs7h/udc7js9ntPnp3f69GmVK1fujkVZkho1aqTA\nwEBZLBZVrlxZw4YN0+bNm285LjU1VfPnz9eMGTPk6ekpi8Wihx56SM7OzvcWJp1BgwbJxcVFTk5O\nmjhxon7++WdduHAh7fn0v8Hfvn27Tp06pZdfflmOjo6qUqWKhgwZos8++0yS9Pnnn+ull15SqVKl\n5O3traeffjrTGQDPP/+8SpUqJR8fH40ZM0aRkZGSpCVLlqSN9Lu4uOiNN97QZ599dtu53YZh6Mkn\nn5Sfn5+KFi2qxx9/XD/99NMt+W9ydnbWwYMHlZiYKGdnZzVr1uy+/9zuByU8D/B09dSqPqs0JnqM\nvjv0ndlxAACAnZQtW1anTp266wLBP/74Q507d5anp6fc3Nz00ksvpU2rSO/UqVO6cuWKqlWrlq1M\nVqtVL7zwgvz8/OTm5qYqVaqkvf7tHDp0SEeOHFHp0qXTft544w2dOHFC0o2pKumnwlSqVCnTDP88\n/siRI5Kko0ePqnLlyhmeS0lJ0fHjx2/7OhUqVEi7XaxYMV28ePGO7zlt2jTZbDYFBgaqbt26mj9/\nfqY5s4MSnkfU86inJT2X6LHPH9Pvp343Ow4AAAWWzZbxJ7PHc/r89Jo2baoiRYpo+fLldzxmxIgR\nql27tvbv36/z589rypQpty3t5cqVU9GiRbV///5bnnNxcVFycnLafavVqpMnT972/SIiIrRy5UrF\nxMTo/Pnz+uuvv/7/57vxof45laRSpUqqUqWKzp49m/aTlJSkVatWSZI8PT31999/px2f/vad/PN4\nLy8vSVLFihV18ODBDM85OjrKw8Mj09dM73bTYTw8PPTxxx8rMTFRs2fP1lNPPZXtrR/vhhKeh7Sr\n1k5TWk9Rp4hOOnnp9v/DAAAABYebm5smTZqkkSNHKioqSsnJybp+/brWrl2r559/XpJ08eJFubq6\nqnjx4vrtt9/04Ycf3va1LBaLBg8erGeffVZHjx6V1WrV1q1bde3aNdWoUUNXrlzRmjVrdP36dU2e\nPFlXr1697etcvHhRRYoUUZkyZXTp0iW9+OKLGZ738PDIUE4DAwPl6uqqadOm6fLly7Jarfrll1+0\nY8eNCxM+/vjjeuONN3Tu3DkdPnw4w6LNO5k+fbrOnTunhIQEvffee+rdu7ckqU+fPnrnnXd08OBB\nXbx4US+++KJCQ0PvOJ3nTtNe3N3dZbFY9Oeff6Y99vnnn+vw4cOSbizYNAzjrtOEsosSnseENQpT\n7zq91e2zbrp8/bLZcQAAQC579tln9fbbb2vy5MkqX768KlWqpFmzZqUt1pw+fboiIiJUsmRJDRs2\nTKGhoRlGctPfnj59uurVq6eAgACVLVtW48ePl81mk5ubm2bNmqUhQ4bI29tbJUqUyDDlwzCMtNcZ\nMGCAKleuLC8vL9WtW1dNmzbN8B5hYWGKj49X6dKl1bNnT1ksFq1atUo//fSTqlatKnd3dw0bNixt\n15KJEyeqcuXKqlKlijp06KABAwZkujCzW7duevDBB9WwYUN17txZgwcPliQNHjxY/fv3V4sWLVS1\nalUVL148Q6n/5+v+88/p5v3ixYvrpZdeUvPmzVWmTBn98MMP2rFjhx566KG0XV/ee+89+fr6Zv5f\nYBYZtgK8N55hGPly679UW6r6ftVXKakpWvroUlkM/q0EAEBWZdYHbva0rFaG7J6PjCwWi/bv36+q\nVauaHeW27vR9ut/eSbvLgyyGRfO7zdfRC0c1fuN4s+MAAFCg3c8c7tw4H4UTJTyPKupYVCtCV2j5\nb8v10Y6PzI4DAABgF/dy2fmCwNHsALizcsXLaU3fNQqaH6SKrhXVtWZXsyMBAADkKqvVanYEu2Ak\nPI/zK+OnqNAoha0M07bD28yOAwAAgBxACc8HAr0CtaDbAnX/rLv+OP2H2XEAAACQTZTwfKJTjU6a\n0nqKOizuoGMXj5kdBwAAANlACc9HwhqFaWCDgeoU0UkXrl6QdGNbpEKyfgEAAKDAYJ/wfMZms2nY\n18OUkJSgr/t8LWdHp///uMnBAADIowpiH4B52Ce8kDIMQx92/lCOFkcN/XqoJP5SAQAA92b58uXy\n8fFRyZIl9dNPP6lu3br69ttvzY5VKFHC8yFHi6OWPrpU8SfjpdYTzI4DAACywdfXV8WLF5erq6vK\nlCmjzp076/Dhw3c8PjY2VhaLRT179szw+M8//yyLxaJWrVrd8dznnntOs2bNUlJSkh544AH98ssv\natGihSQpPDxc/fv3z5kPhUxRwvOZm3PASxRx0fZnVkl1lkqNP0p7nDniAADkL4ZhaNWqVbpw4YKO\nHj0qDw8PjR49+q7nuLu7a9u2bTpz5kzaYwsXLlSNGjXueLEbm82mv//+W7Vr187R/MgaSnh+dqm8\ntDhaajlJqhlldhoAAJBNRYoUUa9evRQfH3/X45ydndW9e3d99tlnkm5c4GbZsmXq27fvbeclX716\nVa6urrJarWrQoIGqV68u6cYofExMjKKjo/XGG29o6dKlcnV1VcOGDXP+wyEDSng+Y7Nl/NHZalJk\nlMqFDdH//t76f48DAIB842ZxTk5O1tKlS9W0adNMz+nfv78WLVokSVq3bp3q1q2rihUr3vbYIkWK\n6OLFi5Kk3bt3a9++fZJujMIbhqEOHTroxRdfVGhoqC5cuKBdu3blxMfCXXDZ+oLgSIAWdl+oHkt7\naPOgzapZrqbZiQAAwD2y2Wzq3r27HB0ddenSJZUvX17R0dGZnte0aVOdOXNGf/zxhxYtWqSBAwcq\nOTk5WznYRcZ+KOEFRMfqHfV6m9f1yJJHtGXwFnm6epodCQCAfMN4NfsLqmwTs1ZgDcNQVFSUWrdu\nLZvNphUrVqhly5aKj4/X1atXVadOnbTjkpKSMpzbv39/vf/++4qNjdWCBQu0ePHibH8O2AclvAAZ\n3HCwjlw4og5LOmjzoM0qVbSU2ZEAAMgXslqgc5phGOrRo4eGDx+uLVu2qGfPnrpw4cIdj+/Xr5+q\nV6+ugQMHqmjRotl+b9gPJTyf++dvjV4KekknLp1Q18iuWtdvnYo5FTMnGAAAuGc3p4HYbDatXLlS\nZ8+elb+/f6bnValSRd9++62qVq2a7QwVKlTQxo0bZbPZKOR2wMLMAsYwDL3b4V15l/RW6JehSklN\nMTsSAADIRJcuXeTq6io3NzdNmDBBixYtumsJT1+SmzVrpgoVKqQ9frcCfbfnHnvsMUlS2bJl1bhx\n4/v9CLhPXLa+gLpmvaaukV3l6eqpeV3n8S9aAEChVZj7AHIel63HXTk7OOvLx7/U3pN79cLGF8yO\nAwAAgHQo4QWYi7OLVj+xWl//8bWm/2+62XEAAADw/7Ews4ArW7ys1vVbp4fnP6xyxctp0AODzI4E\nAABQ6Jk6Eh4dHa1atWqpevXqmjp16i3PL1myRA0aNFD9+vXVvHlz7d69+57Pxf/xcfPRun7r9MLG\nF/T171+bHQcAAKDQM21hptVqVc2aNbVx40Z5eXkpICBAkZGRGVYCb926VbVr15abm5uio6MVHh6u\nbdu23dO5Egsx/ikuMU6dIjrpq8e/UlDlILPjAABgF/QB5KR8vzAzLi5Ofn5+8vX1lZOTk0JDQxUV\nFZXhmKZNm8rNzU2S1KRJEx0+fPiez8WtAr0CFdEzQr2W9dLu47szPwEAAAC5wrQSnpiYKB8fn7T7\n3t7eSkxMvOPxc+fOVceOHbN0Lv5Pu2rtNLPjTD2y5BEdOHvA7DgAAACFkmkLM+9n3+pNmzZp3rx5\n2rJly32fi1s9XudxnUo+pfaftteWwVvkUcLD7EgAAOSa0qVL0x2QY0qXLp0jr2NaCffy8lJCQkLa\n/YSEBHl7e99y3O7duzV06FBFR0enfeh7PVeSwsPD024HBwcrODg4Zz5APvdUwFM6eemkOizpoE0D\nN6lU0VJmRwIAIFecOXPG7AgogGJjYxUbG5vl801bmJmSkqKaNWsqJiZGFStWVGBg4C2LK//++2+1\nbt1aixcv1kMPPXRf50osxMiMzWbT09FP68ejP2p9v/VycXYxOxIAAEC+lG8WZjo6OmrmzJkKCQlR\n7dq11bt3b/n7+2v27NmaPXu2JGnSpEk6e/asRowYoYYNGyowMPCu5+L+GIahdzu8qxpla6jH0h66\nmnLV7EgAAACFgmkj4fbASPi9SUlNUegXoUq1pWrZY8vkaOEaTgAAAPcj34yEI+9wtDhqSc8lunT9\nksJWhinVlmp2JAAAgAKNEg5JUhHHIvrq8a/055k/9fTap/kNAgAAQC6ihCONi7OLVj+xWlsStmjC\npglmxwEAACiwKOHIwK2om9b1W6cv936p/275r9lxAAAACiRW4OEW7i7u2tB/g4LmB6lkkZIa3ni4\n2ZEAAAAKFEo4bsu7pLc29N+glgtaqmSRkupTr4/ZkQAAAAoMSjjuyK+Mn9b1W6e2i9qqhHMJdanZ\nxexIAAAABQJzwnFXdcvX1dd9vlbYyjB989c3ZscBAAAoECjhyFSAV4CWPbZMvb/orW2Ht5kdBwAA\nIN+jhOOeBPsGa2H3heoa2VU7j+40Ow4AAEC+RgnHPetYvaNmd56tjks6as/xPWbHAQAAyLdYmIn7\n0sO/h65ZrylkcYhiBsTI393f7EgAAAD5DiUc96133d66ar2qdp+2U+ygWPmV8TM7EgAAQL5CCUeW\nDGgwQFdSrqjNojbaPGizfEv5mh0JAAAg36CEI8uGPThMV1OuphVx75LeZkcCAADIFyjhyJbRTUbr\nqvWqWi9src2DNsvT1dPsSAAAAHkeu6Mg255r9pwGNBigtp+21clLJ82OAwAAkOdRwpEjXm7xsnrU\n6qF2n7bTmctnzI4DAACQpxk2m81mdojcYhiGCvDHy3NsNpvGbRinzYc2a2P/jXIr6mZ2JAAAALu4\n395JCUeOstlsGr12tHYe3al1/dbJtYir2ZEAAAByHSU8HUq4OVJtqRr+9XDtO7NPq59YLRdnF7Mj\nAQAA5Kr77Z3MCUeOsxgWfdT5I1UuVVldIrso+Xqy2ZEAAADyFEo4coWDxUHzus6Td0lvdY7oTBEH\nAABIhxKOXONgcdD8bvPlXdKbEXEAAIB0KOHIVTeLeEXXiur2WTddvn7Z7EgAAACmo4Qj1zlYHLSg\n2wJ5uHio62ddKeIAAKDQo4TDLhwsDlrYfaHKu5RnRBwAABR6lHDYzc0iXq54OXVf2p0iDgAACi1K\nOOzK0eKoRT0WqUyxMuqxtIeupFwxOxIAAIDdUcJhd44WR33a41OVKlqKIg4AAAolSjhM4Whx1OKe\ni1WySEn1XNqTIg4AAAoVSjhM42hx1JKeS1TCuYR6LeulqylXzY4EAABgF5RwmOpmEXdxcmFqCgAA\nKDQo4TCdk4OTInpFyK2om7pGduXKmgAAoMCjhCNPuLlYs0KJCuoU0UkXr100OxIAAECuoYQjz3C0\nOGp+t/mqWqqqHlnyiC5cvWB2JAAAgFxBCUee4mBx0Jyuc1THvY7aL26v81fOmx0JAAAgx1HCkedY\nDIs+7PShAioGqO2nbXX28lmzIwEAAOQoSjjyJMMwNKPDDLWo1EKtF7XWqeRTZkcCAADIMZRw5FmG\nYWh6++nqUK2DWi1spROXTpgdCQAAIEeYWsKjo6NVq1YtVa9eXVOnTr3l+d9++01NmzZV0aJF9dZb\nb2V4ztfXV/Xr11fDhg0VGBhor8iwM8Mw9Hqb19WzVk8FLwjW0QtHzY4EAACQbY5mvbHVatWoUaO0\nceNGeXl5KSAgQF27dpW/v3/aMWXLltX777+vFStW3HK+YRiKjY1VmTJl7BkbJjAMQ6+2elVODk5q\nuaClvhn4jbxLepsdCwAAIMtMGwmPi4uTn5+ffH195eTkpNDQUEVFRWU4xt3dXY0bN5aTk9NtX8Nm\ns9kjKvKIl1u8rCGNhqjlgpY6dO6Q2XEAAACyzLQSnpiYKB8fn7T73t7eSkxMvOfzDcNQ27Zt1bhx\nY82ZMyc3IiIP+k/z/2h04GgFLwzWgbMHzI4DAACQJaZNRzEMI1vnb9myRZ6enjp58qTatWunWrVq\nKSgoKIfSIS8b89AYFXEoopYLWmp9v/Xyd/fP/CQAAIA8xLQS7uXlpYSEhLT7CQkJ8va+93m+np6e\nkm5MWenRo4fi4uJuW8LDw8PTbgcHBys4ODjLmZF3jAgYIRdnF7Ve1Fprnlijhp4NzY4EAAAKkdjY\nWMXGxmb5fMNm0sTqlJQU1axZUzExMapYsaICAwMVGRmZYWHmTeHh4XJ1ddXYsWMlScnJybJarXJ1\nddWlS5fUvn17TZw4Ue3bt89wnmEYzBsv4L6M/1IjVo/QitAVaubTzOw4AACgkLrf3mnaSLijo6Nm\nzpypkJAQWa1WhYWFyd/fX7Nnz5YkDR8+XMeOHVNAQICSkpJksVg0Y8YMxcfH68SJE+rZs6ekG2W+\nb9++txRwFA69aveSi7OLun3WTUsfXarWVVqbHQkAACBTpo2E2wMj4YXH5oOb9ejnj2p+t/nqXKOz\n2XEAAEAhc7+9kytmokBo6dtSq/qsUtjKMC39ZanZcQAAAO7KtOkoQE5r4t1EG/pvUIfFHZR8PVlP\nNnzS7EgAAAC3RQlHgVLfo742Ddykdp+208VrFzW6yWizIwEAANyCOeEokA6eO6i2i9pqcMPBejHo\nRbPjAACAAu5+eyclHAXWkQtH1O7TdupWs5umtJ6S7QtEAQAA3AklPB1KOE4ln1LI4hA1826mGY/M\nkMVgLTIAAMh5lPB0KOGQpPNXzqtLZBdVcquk+d3my8nByexIAACggKGEp0MJx03J15P1+OePS5KW\nPbZMxZ2Km5wIAAAUJOwTDtxGcafiWt57uUoXK62QxSE6d+Wc2ZEAAEAhRglHoeHk4KSF3ReqUYVG\nCl4QrGMXj5kdCQAAFFKUcBQqFsOidzu8q17+vRQ0P0h/nf3L7EgAAKAQooSj0DEMQxNaTtCYJmMU\nND9Iv5z4xexIAACgkGFhJgq1yD2RGrNujFb0XqGmPk3NjgMAAPIpdkdJhxKOe7F231oNXDFQi3su\nVvtq7c2OAwAA8iF2RwHu0yPVH9Hy3svVf3l/Lft1mdlxAABAIeBodgAgL2heqbk29N+gjks66uzl\nsxreeLjZkQAAQAHGdBQgnT/P/Kn2i9vryQee1EtBL8kwDLMjAQCAfIA54elQwpEVxy4e0yNLHlEz\n72Z675H35GBxMDsSAADI4yjh6VDCkVVJV5PUY2kPlSpaSkt6LlFRx6JmRwIAAHkYCzOBHFCySEmt\neWKNnCxOXOYeAADkOEo4cAdFHIsooleEGlZoqKD5QUpMSjQ7EgAAKCAo4cBdWAyL3gl5RwPqD1Dz\nec219+ResyMBAIACgDnhwD369OdPNW7DOH3V+ys182lmdhwAAJCHsDAzHUo4clr0/mgNWD5Ac7vO\nVZeaXcyOAwAA8ggWZgK5qINfB61+YrWGrRqmT3Z+YnYcAACQTzESDmTBvtP71GFJBy7qAwAAJDEd\nJQNKOHLTzYv6POT1kGZ2nMlFfQAAKMQo4elQwpHbkq4mqdeyXirmWEyRvSLl4uxidiQAAGAC5oQD\ndlSySEmtfmK1yhQro1YLW+n4xeNmRwIAAPkAJRzIJmcHZ83vNl8dq3dU07lN9fup382OBAAA8jim\nowA5aP6u+RofM15fPP6FHq70sNlxAACAnTAnPB1KOMyw/s/16vdVP33Q8QM9Vucxs+MAAAA7oISn\nQwmHWX469pO6RHbR002e1timY9nCEACAAo4Sng4lHGZKOJ+gjhEdFVw5WO92eJctDAEAKMAo4elQ\nwmG281fOq9eyXirhXEIRvSJU3Km42ZEAAEAuYItCIA9xK+qmNX3XyK2om1otbKUTl06YHQkAAOQB\nlHAglzk7OGtBtwXqUK0DWxgCAABJTEcB7Gruzrl68ZsXtezRZWrp29LsOAAAIIcwJzwdSjjyopgD\nMXriqyc0te1UDXpgkNlxAABADqCEp0MJR1619+RedY7srNA6oXqt9WuyGMwMAwAgP8tXCzOjo6NV\nq1YtVa9eXVOnTr3l+d9++01NmzZV0aJF9dZbb93XuUBe5u/ur21h2xR7KFZ9vuyjy9cvmx0JAADY\nkWkj4VarVTVr1tTGjRvl5eWlgIAARUZGyt/fP+2YkydP6tChQ1qxYoVKly6tsWPH3vO5EiPhyPuu\npFzR4KjBOnD2gKJCo+RRwsPsSAAAIAvyzUh4XFyc/Pz85OvrKycnJ4WGhioqKirDMe7u7mrcuLGc\nnJzu+1wgPyjqWFRLei5RSLUQPTT3If164lezIwEAADswrYQnJibKx8cn7b63t7cSExNz/VwgrzEM\nQ6+2elWvtXpNrRa20rr968yOBAAAcplpJdwwDFPOBfKqfvX76cvHv9TAFQP10Y6PbnneMG78AACA\n/M/RrDf28vJSQkJC2v2EhAR5e3vn+Lnh4eFpt4ODgxUcHJylvIA9BFUO0veDv1fniM764/Qf+m+7\n/8rB4mB2LAAA8A+xsbGKjY3N8vmmLcxMSUlRzZo1FRMTo4oVKyowMPC2iyulG0Xa1dU1bWHmvZ7L\nwkzkV2cvn1WvZb3kWsRVS3ouUQnnEmmj4HylAQDIe/LVPuFr167VmDFjZLVaFRYWpvHjx2v27NmS\npOHDh+tOfVu7AAAgAElEQVTYsWMKCAhQUlKSLBaLXF1dFR8frxIlStz23H+ihCM/u2a9pqdWP6Xt\nR7ZrZehK+ZauLIkSDgBAXpSvSnhuo4QjP7sx8m2THnpXaj5NWvaFlND8luP4igMAYD5KeDqUcORn\nGRZh+q2VegyUNkyTfhqU4Ti+4gAAmC/f7BMO4O5stnQ/+x6R5m+WWkzWc+vGKcVqTXsOAADkP4yE\nA/mEYUgqdlqtPnxMxZ2KK6JXhEoWKWl2LAAAIEbCgYLtclmt67dOPiV91GxuMx04e8DsRAAAIAso\n4UA+4+TgpFmdZmlE4xFqNreZNh/cbHYkAABwn5iOAuRjG/7coH7L+2lyq8ka+uBQs+MAAFBosTtK\nOpRwFAZ/nP5DXSK76BG/RzS9/XQ5Wky7EC4AAIUWJTwdSjgKi7OXzyr0y1BJ0me9PlPpYqVNTgQA\nQOHCwkygECpdrLRWP7FaddzrKGBOgH498avZkQAAwF3c00j4pUuXlJCQIMMw5O3tLRcXF3tkyzZG\nwlEYLfp5kZ5b/5w+7vKxutfqbnYcAAAKhRybjnLhwgXNmTNHn332mU6dOiUPDw/ZbDYdP35cZcuW\nVd++fTV06FCVKFEix8LnNEo4CqvtidvVa1kvDW44WK+0fEUWg196AQCQm3KshLdp00ahoaHq2rWr\nPDw8Mjx37NgxrVy5UkuXLlVMTEz2EuciSjgKs2MXj+nRZY/K3cVdi7ovkmsRV7MjAQBQYLEwMx1K\nOAq7a9ZrGr1mtL5P+F5RoVHyK+NndiQAAAqkXCnhZ86c0b59+3T16tW0x1q0aJG1hHZECQdu+GjH\nR3pl0yta1GOROvh1MDsOAAAFTo6X8Dlz5ui9997T4cOH9cADD2jbtm1q2rSpvvnmm2yHzW2UcOD/\nfHfoO/X+oreeeegZPdfsORmGYXYkAAAKjBzfonDGjBmKi4tT5cqVtWnTJu3atUtubm7ZCgnA/oIq\nB+mHIT9o6a9L1fervkq+nmx2JAAACq1MS3jRokVVrFgxSdKVK1dUq1Yt/f7777keDEDO83Hz0XdP\nficHi4MenvewDp07ZHYkAAAKpUxLuI+Pj86ePavu3burXbt26tq1q3x9fe0QDUBuKOZUTIu6L1K/\n+v3U5JMm2vDnBrMjAQBQ6GQ6J/zkyZNyd3eXJMXGxiopKUkdOnSQs7OzXQJmB3PCgbuLPRirPl/2\n0dNNntbzzZ9nnjgAAFmUYwszv/76aw0ePFiOjo5ycHDQ0qVL1bx58xwLag+UcCBzh5MOq9eyXvJy\n9dKC7gtUskhJsyMBAJDv5NjCzBdffFHfffedjh49qi+//FLjx4/PkYAA8hbvkt76dtC3Ku9SXoFz\nArX35F6zIwEAUODdsYQ7OjqqVq1akqQmTZrowoULdgsFwL6KOBbRR50/0n+a/0ctFrTQl/Ffmh0J\nAIACzfFOT5w8eVJvv/122rB6+vuGYejZZ5+1W0gA9jG44WDV96ivXst6KS4xTlPaTJGj5Y5/TQAA\ngCy645zw8PDwDIu0bpbvm/85ceJEu4XMKuaEA1lzKvmU+nzZRzabTZG9IuXu4m52JAAA8rRcuWx9\nfkUJB7LOmmrVy9+8rIhfIvTFY18owCvA7EgAAORZOVbCf/nlF/3555/q1q2bJGnMmDE6f/68DMPQ\nqFGj1KhRo5xJnIso4UD2fbX3Kw1fNVxvtHlDQxoNMTsOAAB5Uo7tjvLCCy+oXLlyaffXr1+vzp07\nKzg4WJMmTcpeSgD5Rk//nvp20Ld6a+tbGhw1mMvdAwCQA+5Ywo8ePZphX3BXV1f16tVLAwYM0MmT\nJ+0SDkDe4O/ur+1Dt+tyymU1ndtU+07vMzsSAAD52h1L+D+3JPzhhx/Sbp84cSL3EgHIk0o4l1BE\nzwgNf3C4ms1rpq/2fmV2JAAA8q07lvCKFStq27Zttzy+detWeXl55WooAHmTYRh6KuAprX5itZ5d\n96zGrhur69brZscCACDfuePCzLi4OPXu3VuDBg1So0aNZLPZtHPnTi1YsEBLly5VkyZN7J31vrEw\nE8g9p5NPq//y/kq6mqSljy6VV0n+cQ4AKLxydIvC48ePa+bMmYqPj5ck1alTRyNHjpSHh0f2k9oB\nJRzIXam2VL3x3RuauX2mFvdYrDZV25gdCQAAU7BPeDqUcMA+Yg7EqN/yfhoVMErjg8bLYtxxphsA\nAAVSjm1R2KlTJ33++edKTr51O7JLly5p6dKl6tixY9ZSAihQ2lRtox1Dd2jt/rXqHNFZp5NPmx0J\nAIA87Y4j4SdOnNDMmTP1xRdfyMHBQZ6enrLZbDp27JhSUlLUu3dvjRw5Uu7uefdy1oyEA/Z13Xpd\n42PG64v4L7TssWUK9Ao0OxIAAHaRK9NRjh07pkOHDkmSKleurAoVKmQ9oR1RwgFzLN+7XMNXDdf4\nh8drzENjZBiG2ZEAAMhVzAlPhxIOmOfA2QMK/SJUnq6emt9tvsoUK2N2JAAAck2OzQkHgOyoWrqq\nvh/8vaqWqqpGsxtp2+FbrzsAAEBhxUg4gFwX9VuUhq0apnHNxunZps+yewoAoMDJlZHw5ORk/f77\n71kOdSfR0dGqVauWqlevrqlTp972mH//+9+qXr26GjRooF27dqU97uvrq/r166thw4YKDGTxF5CX\ndavVTXFD4vRF/BfqGtmV3VMAAIVepiV85cqVatiwoUJCQiRJu3btUteuXbP9xlarVaNGjVJ0dLTi\n4+MVGRmpvXv3ZjhmzZo12r9/v/bt26ePP/5YI0aMSHvOMAzFxsZq165diouLy3YeALmrcqnK+vbJ\nb1WrXC01nN1QW/7eYnYkAABMk2kJDw8P1w8//KDSpUtLkho2bKgDBw5k+43j4uLk5+cnX19fOTk5\nKTQ0VFFRURmOWblypQYOHChJatKkic6dO6fjx4+nPc9UEyB/cXZw1vT20zWr0yz1XNZTb37/plJt\nqWbHAgDA7jIt4U5OTipVqlTGkyzZn8+ZmJgoHx+ftPve3t5KTEy852MMw1Dbtm3VuHFjzZkzJ9t5\nANhP5xqdtWPoDn39x9fqFNFJJy+dNDsSAAB2lWmbrlOnjpYsWaKUlBTt27dPo0ePVrNmzbL9xve6\nb/CdRru///577dq1S2vXrtUHH3yg7777LtuZANiPj5uPYgfGqoFHAzX6uJFiD8aaHQkAALtxzOyA\nmTNnavLkySpSpIj69OmjkJAQTZgwIdtv7OXlpYSEhLT7CQkJ8vb2vusxhw8flpeXlySpYsWKkiR3\nd3f16NFDcXFxCgoKuuV9wsPD024HBwcrODg429kB5AwnBye92fZNBfsG64kvn9CQRkP0SstX5GjJ\n9K8mAABMFRsbq9jY2Cyff9ctClNSUtSuXTtt2rQpy29wt9euWbOmYmJiVLFiRQUGBioyMlL+/v5p\nx6xZs0YzZ87UmjVrtG3bNo0ZM0bbtm1TcnKyrFarXF1ddenSJbVv314TJ05U+/btM344tigE8o1j\nF49pwPIBunT9kiJ6RqhyqcpmRwIA4J7db++863CTo6OjLBaLzp07d8u88OxydHTUzJkzFRISIqvV\nqrCwMPn7+2v27NmSpOHDh6tjx45as2aN/Pz85OLiovnz50uSjh07pp49e0q6Ueb79u17SwEHkL9U\nKFFB0f2i9db/3lLAnADN6jRLj9Z+1OxYAADkikwv1tO1a1ft2rVL7dq1k4uLy42TDEPvvfeeXQJm\nByPhQP60PXG7+nzZR22qtNE7Hd5RcafiZkcCAOCu7rd3ZlrCFyxYcNs3ubl1YF5GCQfyr6SrSXpq\n9VPaeXSnPnv0M9X3qG92JAAA7ijHS3h+RgkH8r9Pf/5Uz65/VuEtw/VUwFP3vLMSAAD2lOMlvEqV\nKrd9k5y4YE9uo4QDBcO+0/sU+mWofEr6aG7XuSpbvKzZkQAAyCBHF2ZK0vbt29NuX7lyRV988YVO\nnz6dtXQAkAXVy1bX1rCtejHmRT0w+wEt7rFYLX1bmh0LAIAsy9J0lEaNGmnnzp25kSdHMRIOFDzR\n+6M1OGqwBjYYqFdbvSpnB2ezIwEAkPPTUX788ce0OZipqanasWOHPvzwQ/3888/ZS2oHlHCgYDpx\n6YQGRw3WsYvHtKTnEtUsV9PsSACAQi7HS3hwcHBaCXd0dJSvr6+ee+451ayZ9/9PjxIOFFw2m00f\n7fhIr8S+osmtJmvYg8NYtAkAMA27o6RDCQcKvr0n96rvV33l4+ajT7p8IncXd7MjAQAKofvtnZbM\nDpgxY4aSkpJks9kUFhamRo0aad26ddkKCQA5xd/dX9uGbFOtsrX0wOwHFL0/2uxIAABkKtMSPnfu\nXJUsWVLr16/XmTNntGjRIr3wwgv2yAYA98TZwVlT203V4h6LNezrYXp67dO6fP2y2bEAALijTEv4\nzWH11atXq3///qpbt26uhwKArGhVpZV+/tfPOnrxqALmBGj38d1mRwIA4LYyLeEPPvig2rdvrzVr\n1igkJERJSUmyWDI9DQBMUbpYaS19dKnGNRunNova6J2t7yjVlmp2LAAAMsh0YWZqaqp27dqlatWq\nqVSpUjp9+rQSExNVv359e2XMMhZmAoXbgbMH1O+rfiruVFzzu82Xj5uP2ZEAAAVUji/M3Lp1q2rW\nrKlSpUrp008/1eTJk+Xm5patkABgD1VLV9W3T36rVr6t9ODHD2rx7sX8wxwAkCdkOhJer1497d69\nW7t379agQYM0ZMgQLVu2TJs3b7ZXxixjJBzATTuP7lT/5f1V2722Pur0kcoWL2t2JABAAZLjI+GO\njo4yDEMrVqzQyJEjNXLkSF24cCFbIQHA3hp5NtKPw35UpZKVVP+j+lq7b63ZkQAAhVimJdzV1VWv\nv/66Fi9erM6dO8tqter69ev2yAYAOaqoY1G9FfKWlvRcohGrR+hfq/6li9cumh0LAFAIZVrCly5d\nqqJFi2revHmqUKGCEhMTNW7cOHtkA4BcEewbrJ//9bOupFzRAx89oK0JW82OBAAoZO7psvUHDx7U\n/v371bZtWyUnJyslJUUlS5a0R75sYU44gMws37tcI1aPUFjDME0MnihnB2ezIwEA8qEcnxP+8ccf\n67HHHtPw4cMlSYcPH1aPHj2ynhAA8pAe/j30879+1p4Te9Tkkyb65cQvZkcCABQCmZbwDz74QN9/\n/33ayHeNGjV04sSJXA8GAPbiUcJDUaFRGhUwSq0WttK0LdNkTbWaHQsAUIBlWsKLFCmiIkWKpN1P\nSUmRYRi5GgoA7M0wDIU1CtP2odsVvT9aD89/WL+f+t3sWACAAirTEt6yZUtNmTJFycnJ2rBhgx57\n7DF16dLFHtkAwO58S/lq44CN6levnx6e/7De2foOo+IAgBx3T5et/+STT7R+/XpJUkhIiIYMGZIv\nRsNZmAkgO/4886eejHpSNtk0v9t8+ZXxMzsSACCPut/eedcSnpKSorp16+q3337LkXD2RgkHkF2p\ntlS998N7mvLdFE1sOVFPBTwli5HpLxEBAIVMju6O4ujoqJo1a+rQoUPZDgYA+ZHFsGjMQ2O0ZfAW\nReyJUJtFbfTX2b/MjgUAyOcynY4SFBSkXbt2KTAwUC4uLjdOMgytXLnSLgGzg5FwADnJmmrV21vf\n1rT/TdPkVpM17MFh+WJqHgAg9+XodBRJ2rx5syRleFHDMNSyZcssRrQfSjiA3BB/Ml6DVgySW1E3\nze06V5XcKpkdCQBgshwr4ZcvX9ZHH32k/fv3q379+ho8eLCcnJxyLKg9UMIB5JaU1BRN2zJN72x7\nR5NbTdbQB4cyVxwACrEcK+GPP/64nJ2dFRQUpDVr1sjX11czZszIsaD2QAkHkNt+PfGrBq8cLBcn\nF83pMkfVylQzOxIAwAQ5VsLr1aunPXv2SLqxS0pAQIB27dqVMynthBIOwB6sqVa9u+1dvfH9G5rQ\nYoJGBY6Sg8XB7FgAADvKsd1RHB0db3sbAJCRg8VBY5uN1dawrfpy75dqsaCFfjuVP7d2BQDYxx1H\nwh0cHFS8ePG0+5cvX1axYsVunGQYSkpKsk/CbGAkHIC9pdpS9eH2DzUxdqLGNRunsc3GytHCQAYA\nFHQ5vjtKfkYJB2CWg+cOaujXQ3X28lnN7zZf9TzqmR0JAJCLcvRiPQCArPEt5av1/dZrROMRar2o\ntV6NfVXXrNfMjgUAyCMYCQeAXHY46bD+tepf+vv835rbda4CvALMjgQAyGFMR0mHEg4gr7DZbIrY\nE6Gx68eqb72+mtRqklycXcyOBQDIIUxHAYA8yDAM9a3fV3tG7NGJ5BOq+2Fdrdu/zuxYAACTMBIO\nACZYt3+d/rX6XwqqFKS3Q95WueLlzI4EAMgGRsIBIB8I8QvRnhF7VK54OdWdVVcReyIYNACAQsTU\nEh4dHa1atWqpevXqmjp16m2P+fe//63q1aurQYMGGa7YeS/nAkBeVsK5hN4OeVtf9/laU7dMVaeI\nTjp07pDZsQAAdmBaCbdarRo1apSio6MVHx+vyMhI7d27N8Mxa9as0f79+7Vv3z59/PHHGjFixD2f\nCwD5RYBXgHYM3aGgSkF68OMHNWPbDFlTrbccZxg3fgAA+Z9pJTwuLk5+fn7y9fWVk5OTQkNDFRUV\nleGYlStXauDAgZKkJk2a6Ny5czp27Ng9nQsA+YmTg5PGB43X/8L+p+W/LVezec205/ges2MBAHKJ\naSU8MTFRPj4+afe9vb2VmJh4T8ccOXIk03MBID+qUbaGvhn4jYY0HKI2i9po/MbxSr6ebHYsAEAO\nM62EG/f4O1UWKgEobCyGRUMfHKrdI3br0PlDqjurrqL3R5sdCwCQgxzNemMvLy8lJCSk3U9ISJC3\nt/ddjzl8+LC8vb11/fr1TM+9KTw8PO12cHCwgoODc+YDAEAu83StIClCqrZOjxx4Sno0QIp+R4bh\nmeE4xioAwP5iY2MVGxub5fNN2yc8JSVFNWvWVExMjCpWrKjAwEBFRkbK398/7Zg1a9Zo5syZWrNm\njbZt26YxY8Zo27Zt93SuxD7hAPK3DL8wdEqWWkyWGs2RNk2Sfhwu2W78MpO/5gDAfPfbO00bCXd0\ndNTMmTMVEhIiq9WqsLAw+fv7a/bs2ZKk4cOHq2PHjlqzZo38/Pzk4uKi+fPn3/VcAChIMv5dXlyG\n8bq0u6+avzlcVtsize48W/U96psVDwCQDVwxEwDyiZsj49bUVM3dOVcvffOSnnzgSb3S8hW5OLuY\nGw4ACjmumAkABdzNhZt7RuxR4oVE1f2wrtbsW2N2LADAfWAkHADyiZsj4f/8a23Dnxs0YvUIPVDh\nAb0T8o583HxuPRkAkKsYCQeAAspmu/0izHbV2mnPiD2qW76uGs5uqOn/m67r1uv2DwgAuGeMhANA\nAbL/zH6NWjNKh5MOa1anWWpRuYXZkQCgULjf3kkJB4ACxmaz6au9X2nMujFqXaW1prWdJo8SHmbH\nAoACjekoAFDIGYahXrV7ae/IvfJw8VDdD+tq1vZZsqZazY4GAPj/GAkHgALulxO/6KnVTyn5erI+\n7PShArwCzI4EAAUO01HSoYQDwA02m02f7v5Uz298Xj1q9dCU1lNUulhps2MBQIHBdBQAwC0Mw9CA\nBgMU/1S8DBny/8Bf83bNU6ot1exoAFAoMRIOAIXQj0d+1Ki1o5RqS9XMR2YyRQUAsonpKOlQwgHg\nzlJtqVr08yKNjxmvLjW66PU2r6tc8XJmxwKAfInpKACAe2IxLBr0wCD9NvI3uTi5qPYHtfVB3AdK\nSU0xOxoAFHiMhAMAJN3YRWX02tE6e/msZnacqYcrPWx2JADIN5iOkg4lHADuj81m07Jfl+m5Dc8p\n2DdY09pOk6erp9mxACDPYzoKACDLDMNQ77q9tXfkXnm7eqveh/X01v/e0nXrdbOjAUCBwkg4AOCO\n/jj9h56OfloHzx3UOyHvqINfB7MjAUCexHSUdCjhAJB9NptNq/et1jPrnlHNsjX1dsjbqlG2htmx\nACBPYToKACBHGYahzjU665cRv6hl5ZZqNreZxq0fp/NXzpsdDQDyLUo4AOCeFHEsonHNx+mXp37R\nmctnVOuDWpq7cy5X3QSALGA6CgAgS3Yc2aGno5/W1ZSrmtFhhppXam52JAAwDXPC06GEA0Dustls\nivwlUs9vfF5BlYI0rd00eZf0NjsWANgdc8IBAHZjGIaeqPeEfhv5m6qVrqYGHzXQpM2TlHw92exo\nAJCnMRIOAMgxB88d1LgN4/TD4R/0Rps31KdeH1kMxnsAFHxMR0mHEg4A5vj+7+/1zLpnZDEserv9\n28wXB1DgUcLToYQDgHlSbalasnuJXvzmRTXzaaY327ypKqWrmB0LAHIFc8IBAHmCxbCof4P++n3U\n76rjXkeN5zTWCxtfUNLVJLOjAYDpKOEAgFxV3Km4Xmn5ivaM2KPjl46rxvs1NHvHbKWkppgdDQBM\nw3QUAIBd7Ty6U8+ue1anL5/WW+3fUvtq7c2OBADZxpzwdCjhAJA32Ww2rfhthcZtGKcaZWtoWrtp\nqlu+rtmxACDLmBMOAMjzDMNQD/8eih8Zr/bV2qv1wtYaunKojl44anY0ALALSjgAwDTODs4a89AY\n/T7qd5UuVlp1P6yriZsm6uK1i2ZHA4BcRQkHAJiudLHSmtZumn4c9qP2n93P4k0ABR5zwgEAec6O\nIzs0bsM4Hb94XFPbTlXnGp1lGIbZsQDgjliYmQ4lHADyL5vNptX7Vus/G/6j8i7lNb39dDWu2Njs\nWABwW5TwdCjhAJD/paSmaN6ueQqPDVewb7CmtJ7ClTcB5DnsjgIAKFAcLY4a9uAw/TH6D9UoW0ON\n5zTWmOgxOnnppNnRACDLKOEAgHyhhHMJhQeHK/6peKWkpqjWB7X02ubX2EkFQL5ECQcA5CseJTw0\ns+NM/TDkB8Wfilf196tr1vZZum69bnY0ALhnzAkHAORrPx75UeNjxuuvc39pSusperT2o7IYjDEB\nsC8WZqZDCQeAwmPDnxv0QswLshgWvdnmTbWp2sbsSAAKkXyxMPPMmTNq166datSoofbt2+vcuXO3\nPS46Olq1atVS9erVNXXq1LTHw8PD5e3trYYNG6phw4aKjo62V3QAQB7Vrlo7bR+6Xc81fU7DVw1X\nyOIQ7Tq6y+xYAHBbppTwN998U+3atdMff/yhNm3a6M0337zlGKvVqlGjRik6Olrx8fGKjIzU3r17\nJd34l8azzz6rXbt2adeuXerQoYO9PwIAIA+yGBb1rttb8SPj1bVGV3WM6KjeX/TWH6f/MDsaAGRg\nSglfuXKlBg4cKEkaOHCgVqxYccsxcXFx8vPzk6+vr5ycnBQaGqqoqKi055lmAgC4E2cHZ40MHKl9\no/epgUcDNZvbTENWDlHC+QSzowGAJJNK+PHjx+Xh4SFJ8vDw0PHjx285JjExUT4+Pmn3vb29lZiY\nmHb//fffV4MGDRQWFnbH6SwAgMKthHMJvRj0ovaN3qfyLuX1wOwH9Ez0Mzpx6YTZ0QAUcrlWwtu1\na6d69erd8rNy5coMxxmGIcMwbjn/do/dNGLECP3111/66aef5OnpqbFjx+Z4fgBAwVG6WGm93uZ1\n/frUr0pJTZH/B/56ZdMrOn/lvNnRABRSjrn1whs2bLjjcx4eHjp27JgqVKigo0ePqnz58rcc4+Xl\npYSE//u1YUJCgry9vSUpw/FDhgxRly5d7vhe4eHhabeDg4MVHBx8H58CAFCQVChRQe93fF9jm41V\neGy4qr9fXc81e06jAkepuFNxs+MByEdiY2MVGxub5fNN2aLwP//5j8qWLavnn39eb775ps6dO3fL\n4syUlBTVrFlTMTExqlixogIDAxUZGSl/f38dPXpUnp6ekqR33nlH27dvV0RExC3vwxaFAIC7iT8Z\nrwmbJmjb4W16OehlhTUKk7ODs9mxAORD+WKf8DNnzujxxx/X33//LV9fXy1btkylSpXSkSNHNHTo\nUK1evVqStHbtWo0ZM0ZWq1VhYWEaP368JGnAgAH66aefZBiGqlSpotmzZ6fNMU+PEg4AuBfbE7fr\npW9e0r4z+zSx5UT1q99PjpZc+2UxgAIoX5Rwe6GEAwDux7eHvtWETRN07OIxTWw5Ub3r9JaDxcHs\nWADyAUp4OpRwAMD9stlsivkrRhM2TVDS1SRNCp6kHv49ZDFM2VAMQD5BCU+HEg4AyCqbzaa1+9dq\nwqYJSrWlalLwJHWu0fmuu3cBKLwo4elQwgEA2WWz2RT1e5QmbJqg4k7FNSl4ktpXa08ZB5ABJTwd\nSjgAIKek2lL1+a+fK3xzuMoVL6dJwZPUqkors2MByCMo4elQwgEAOc2aalXEngi9uvlVVXKrpIkt\nJ6qlb0uzYwEwGSU8HUo4ACC3pKSmaPHuxZr87WTKOABKeHqUcABAbqOMA5Ao4RlQwgEA9kIZBwo3\nSng6lHAAgL1RxoHCiRKeDiUcAGAWyjhQuFDC06GEAwDMlpKaoiW7l2jyd5PlWcJTE1pMUNuqbdln\nHChgKOHpUMIBAHlFSmqKlv6yVFO+m6KSRUpqQosJ6li9I2UcKCAo4elQwgEAeY011aqv9n6l1759\nTU4OTno56GV1q9VNFsNidjQA2UAJT4cSDgDIq1JtqVr5+0q99u1ruma9ppeDXtajtR+Vg8XB7GgA\nsoASng4lHACQ19lsNq3dv1avffuazl4+q5eCXlKfen3kaHE0OxqA+0AJT4cSDgDIL2w2m2L+itFr\n376mw0mH9ULzFzSgwQAVcSxidjQA94ASng4l/P+1d+9RVZf5Hsc/W0S8FoiC11EzQS47QMVLZjoq\n5LTU1DMikCOh4S2PzpQrq1lnltN00cZKcqY0U9Fg8K5d7GqT2tFMRFAxkLwgNgIqqIOKqLDPHx2Y\njYCCwv7B5v1a67eE335+P77b9fj44bee/TwAgPpo16ldeu2715RyNkVzH56rqF5RatGkRbl2JZ/p\n5GhDVagAABdUSURBVL86wHiEcCuEcABAfZZ4JlGv/+/r2nVql2b3m61nAp+RSzOX0tcJ4UDdQQi3\nQggHANiD1HOpWrh7oT5J/0RRvaL0h/5/kHtLd0I4UIdUN3eyHhIAAHWcV1svxYyJUeLURF2+flle\nf/fSrM9mSfefMro0AHeJJ+EAANQTpfv6tMyW+i+Wei2X0kdJ//uCdL5naTv+6wNsj+koVgjhAAB7\nUm5zzaYXpL5/l/q9I2UOlHbPk37uTwgHDEAIt0IIBwDYs5JQfrnwilYmrdSb37+prs5dNW/gPI14\ncIRM5VI7gNpCCLdCCAcA2LNbP5h5o+iG1h1Zpzd2vyGTyaR5A+cpxCeEjX8AGyCEWyGEAwDsWWWr\no5Tswrlw90JlXsrUcwOe0+SAyWru2Nz2RQINBCHcCiEcAGDPqrJE4d6f92rh7oXac3qPngl8Rs8E\nPiPX5q62KRBoQAjhVgjhAAD8IvVcqv6656/amrZVk/wm6ff9f6+uzl2NLguwG6wTDgAAyvFq66WV\nT6zUoRmH1MShiXq/31thm8KUeCbR6NKABokn4QAANED/Lvy3licu1+IfFqtH6x6a+/BcjXhwhBqZ\neD4H3A2mo1ghhAMAcHslK6os2rNIN4pvaO6AuQo3h8upsZPRpQH1CiHcCiEcAICqsVgs2n5iuxZ9\nv0iHcw5rdr/ZmtZ7mlyauRhdGlAvEMKtEMIBAKi+g9kH9eb3b+rT9E81yW+S5vSbo24u3YwuC6jT\n+GAmAAC4J37t/LRm7BodmnFIjo0c1Wd5H43fMF7fn/7e6NIAu8GTcAAAcFv5hflalbxKi/culntL\ndz3b/1mN9RrLTpyAFaajWCGEAwBQc4qKi/Tx0Y/11t63dPrSac3uN1tTAqbo/qb3G10aYDhCuBVC\nOAAAtSPhXwl6e+/b+vL4l4rwi9DsfrPZ/AcNGiHcCiEcAIDadfrSaS3Zt0QrklZoaLeh+kP/P2hA\npwEymUxGlwbYFCHcCiEcAADbyC/MV0xyjKJ/iJZLMxfN6TdHIT4hauLQxOjSAJsghFshhAMAYFtF\nxUX67KfPFP1DtH4896Nm9JmhaX2mya2Fm9GlAbWKEG6FEA4AgHFSzqbonR/e0YYfN2hMzzGa02+O\n/Nv5G10WUCvqxTrheXl5CgoKkoeHh4KDg3Xx4sUK202ePFnu7u4ym813dT0AADCOr5uv3h/1vn76\n75/k0dpDo+JHaXDMYG1O3ayi4iKjywMMZUgIX7BggYKCgpSenq5hw4ZpwYIFFbaLjIzUF198cdfX\nAwAA47Vp3kYvDnpRJ2af0Mw+M7VozyJ1f6e7Fu1ZpAsFF4wuDzCEIdNRevbsqZ07d8rd3V3Z2dka\nMmSI0tLSKmybkZGhUaNG6fDhw9W+nukoAADUTQn/SlD0D9Ha9tM2hXiHaFbfWTK7m+98IVBH1Yvp\nKDk5OXJ3d5ckubu7Kycnx6bXAwAAYwV2DFTsuFilPpOqTvd10oi4ERoSM0Sbftykm8U3jS4PqHW1\ntt9sUFCQsrOzy51/9dVXy3xvMpnuaS3Re70eAAAYp13Ldvqfwf+jFx55QVvStmjxD4v1+y9/rxl9\nZujpXk+zqgrsVq2F8K+//rrS10qmkbRr105ZWVlyc6veP7DqXD9//vzSr4cMGaIhQ4ZU62cBAIDa\n5+jgqBCfEIX4hCg5O1l/2/c3ef7NU6M9R2tW4CwFdgw0ukSgjB07dmjHjh13fb0hc8Kff/55ubq6\nat68eVqwYIEuXrxY6YcrK5oTXtXrmRMOAED9lXs1VyuTVurd/e/KvYW7ZvWdpfHe4+XU2Mno0oBy\n6sU64Xl5eQoJCVFmZqa6du2q9evXy9nZWWfOnFFUVJS2bdsmSQoLC9POnTuVm5srNzc3vfzyy4qM\njKz0+nJvjhAOAEC9V1RcpG0/bdOSfUt0KOeQJvtP1rQ+09TVuavRpQGl6kUItxVCOAAA9iU9N11L\n9y/VmoNrNKDzAM3sM1OPPfiYGpkMWWsCKEUIt0IIBwDAPl29cVVrU9bq7wl/18VrFzW993RFBkSq\nTfM2RpeGBooQboUQDgCAfbNYLEo4k6B3E97VR0c/0mjP0ZrZZ6b6duzL6mmwKUK4FUI4AAANR+7V\nXK1KXqX39r8n56bOmtFnhsJ8w9SiSQujS0MDQAi3QggHAKDhKbYU68tjX2pZ4jLtOrVL4eZwTes9\njR05UasI4VYI4QAANGynL53WiqQV+uDAB+ri3EXTe0/Xb71/q2aOzYwuDXaGEG6FEA4AACTpZvFN\nbUvfpqWJS5XwrwRN8pukqb2nqmebnkaXBjtBCLdCCAcAALc6ceGElicu16rkVfJq66XpvadrrNdY\nNXFoYnRpqMcI4VYI4QAAoDLXi65ra9pWLd2/VEfOHVGEX4SiekWph2sPo0tDPUQIt0IIBwAAVXH0\n/FF9cOADrTm0Rt5tvTW111SN9Rqrpo2bGl0a6glCuBVCOAAAqI7rRdf1UdpHWn5guZKykzTRPFFR\nvaPk3dbb6NJQxxHCrRDCAQDA3Tpx4YRWHFihVcmr9IDLA4rqFaXxPuPV3LG50aWhDiKEWyGEAwCA\ne1Wyssr7B97X3p/3KtQnVFN7T5VfOz+jS0MdQgi3QggHAAA16fSl01qZtFIrklbIrYWbpgRMUZg5\nTM5NnY0uDQYjhFshhAMAgNpQVFyk7Se2a0XSCn11/CuN8hylKQFTNLjLYJlMJqPLgwEI4VYI4QAA\noLadv3pesYditSJphQpuFGhywGRF+EWo430djS4NNkQIt0IIBwAAtmKxWJRwJkErDqzQhh836OHO\nD2tKwBSN9BgpRwdHo8tDLSOEWyGEAwAAI1y5fkUbf9yoFUkrdDT3qCaaJyoyIFK+br5Gl4ZaQgi3\nQggHAABGS89NV0xyjNYcXKP2rdor0j9Sob6hat2stdGloQYRwq0QwgEAQF1RVFykr098rVXJq/TF\nsS804sERivSPVNADQXJo5GB0ebhHhHArhHAAAFAX5RXkaW3KWq1KXqWs/CxN8pukp/yfkoerh9Gl\n4S4Rwq0QwgEAQF2XcjZFMckxij0Uq+6tuyvSP1Ljvcfr/qb3G10aqoEQboUQDgAA6osbRTf0+bHP\nFZMco29OfqPHezyuCL8IDX9guBo3amx0ebgDQrgVQjgAAKiPzl89r3Up67T64Gr9/O+f9aT5SUX4\nR7C6Sh1GCLdCCAcAAPVd6rlUrTm4Rh8e+lBuLdwU4RehcHO42rZoa3RpsEIIt0IIBwAA9qKouEjf\nZnyr1QdX65Ojn+jRLo8qwi9CIz1Gyqmx022vNZl++ZNYVHsI4VYI4QAAwB7lF+ZrU+omrT64Wody\nDum3Xr/V7/x+p4GdB8pUkritEMJrHyHcCiEcAADYu8xLmYo7FKcPD32ogpsFmmieqIkPTZRnG8/S\nNoTw2kcIt0IIBwAADYXFYlFSdpJiD8UqPiVene/rrN899DtN8J0g95Zu/9/G4CLtGCHcCiEcAAA0\nRCaHm1K3byS/DyWPT6XMR6RDE6W0J6SbzUrbEZNqDiHcCiEcAAA0RGWmhTe5LPXcIj0UK3XcJx19\nQjocLp0cKksR64/XFEK4FUI4AADAf0L5mX9nad2RdfrH4X8o81KmJvhM0JMPPanADoEVfqATVUcI\nt0IIBwAAqPiDmem56Yo/HK+4w3GyyKJw33CFm8PLfKATVUcIt0IIBwAAuP3qKBaLRYlZiYo7FKe1\nR9aqY6uOCjeHK9Q3VB1adbBtofUYIdwKIRwAAKDqSxSWbAj0j8P/0Na0rQpoH6BQn1D9l/d/qXWz\n1rVfaD1GCLdCCAcAALg7BTcK9NlPn2ntkbX66vhXGvSrQQrzDdNoz9Fq5dTK6PLqHEK4FUI4AADA\nvcsvzNdHRz/S2pS1+i7zOz3W/TGF+obq8R6Pq2njpkaXVycQwq0QwgEAAGpW7tVcbU7drPiUeCVl\nJ2m052iF+oRq+APD5ejgaHR5hiGEWyGEAwAA1J6s/CytP7Jea4+s1bG8Y1r1xCqN9BhpdFmGIIRb\nIYQDAADYxskLJ9XcsbncW7obXYohCOFWCOEAAACwhermzka1WEul8vLyFBQUJA8PDwUHB+vixYsV\ntps8ebLc3d1lNpvLnJ8/f746deqkgIAABQQE6IsvvrBF2QAAAECNMCSEL1iwQEFBQUpPT9ewYcO0\nYMGCCttFRkZWGLBNJpOeffZZJSUlKSkpSSNGjKjtkmFHduzYYXQJqIPoF6gI/QIVoV+gJhgSwj/+\n+GNFRERIkiIiIrR169YK2w0aNEguLi4VvsY0E9wtBk9UhH6BitAvUBH6BWqCISE8JydH7u6/TNp3\nd3dXTk5Ote+xZMkS+fn5acqUKZVOZwEAAADqoloL4UFBQTKbzeWOjz/+uEw7k8kkU8leqlU0Y8YM\nnTx5UsnJyWrfvr2ee+65miwdAAAAqF0WA3h6elqysrIsFovFcubMGYunp2elbU+ePGnx9fW9q9e7\nd+9ukcTBwcHBwcHBwcFRq0f37t2rlYcbywCjR4/W6tWrNW/ePK1evVpjxoyp1vVZWVlq3769JGnL\nli3lVk8pcezYsXuuFQAAAKhphqwTnpeXp5CQEGVmZqpr165av369nJ2ddebMGUVFRWnbtm2SpLCw\nMO3cuVO5ublyc3PTyy+/rMjISE2aNEnJyckymUzq1q2bli1bVjrHHAAAAKjr7HqzHgAAAKAuMmR1\nlJp29OjR0o17AgICdP/99+udd96p8qZAsE8V9Yvo6Gg2e4Jef/11+fj4yGw2Kzw8XIWFhYwXqLBf\nMF4gOjpaZrNZvr6+io6OllT1TQdhvyrqF9UdL+zuSXhxcbE6duyoffv2acmSJWrTpo2ef/55LVy4\nUBcuXKh0YyDYN+t+sXLlSrVq1UrPPvus0WXBABkZGRo6dKhSU1Pl5OSkCRMm6PHHH9eRI0cYLxqw\nyvpFRkYG40UDlpKSorCwMCUkJMjR0VEjRozQ0qVLtWzZMsaLBqyyfhEbG1ut8cIunoRb2759ux58\n8EF17ty5ypsCwf5Z9wuLxcJmTw3YfffdJ0dHR129elU3b97U1atX1aFDB8aLBq6iftGxY0dJYrxo\nwNLS0tSvXz81bdpUDg4OGjx4sDZt2sR40cBV1C82b94sqXrjhd2F8LVr1yosLExSzWwKBPtg3S9M\nJhObPTVgrVu31nPPPadf/epX6tChg5ydnRUUFMR40cBV1C+GDx8uic3hGjJfX1999913ysvL09Wr\nV/XZZ5/p559/Zrxo4CrqF6dPn5ZUvfHCrkL49evX9cknn2j8+PHlXrubTYFgH27tF2z21LAdP35c\nixcvVkZGhs6cOaPLly8rNja2TBvGi4anon4RFxfHeNHA9ezZU/PmzVNwcLB+85vfyN/fXw4ODmXa\nMF40PJX1i5kzZ1ZrvLCrEP7555+rd+/eatu2raRffjvNzs6W9Mva4m5ubkaWB4Pc2i/c3NxKB82n\nn35a+/btM7hC2NL+/fv18MMPy9XVVY0bN9a4ceP0/fffq127dowXDVhF/WLPnj2MF9DkyZO1f/9+\n7dy5Uy4uLvLw8CBfoEy/cHZ2lqenp9q2bVut8cKuQnh8fHzplAPpP5sCSbqrTYFgH27tF1lZWaVf\n326zJ9innj17au/evSooKJDFYtH27dvl7e2tUaNGMV40YJX1i5KgJTFeNFRnz56VJGVmZmrz5s0K\nDw8nX6BMv9iyZYvCw8OrnS/sZnWUK1euqEuXLjp58qRatWolqfJNgdBwVNQv2OwJb7zxhlavXq1G\njRqpV69e+uCDD5Sfn8940cDd2i+WL1+up59+mvGigXv00UeVm5srR0dHvf322/r1r39NvkCF/aK6\n+cJuQjgAAABQX9jVdBQAAACgPiCEAwAAADZGCAcAAABsjBAOAAAA2BghHAAAALAxQjgAAABgY4Rw\nAAAAwMYI4QCAeiklJUU7duzQSy+9ZHQpAFBthHAAQL107Ngx9ejRo3T7aACoTwjhAIB6acyYMcrP\nz1dgYKDRpQBAtRHCAeAeZGRkyGw2lzk3f/58vfnmm7e9buDAgZKkS5cu6b333quxelq1aiVJatmy\nZY3cr6buUx2FhYUaPHiwLBZLudeuX7+uiRMnSpLeeOMNde7cWceOHVN6enq5ezz66KMqLi62Sc0A\nUF2EcACoYSaT6Y5tdu/eLUm6cOGC3n33XUNqsOV9qiMuLk4jR46s8GevW7dOe/fulfTLLzIHDhyQ\nk5OTunTpUqadk5OTBg0apK1bt9qkZgCoLkI4ANSSU6dOycvLS1OnTpWvr68ee+wxXbt2TdJ/njC/\n8MILOn78uAICAjRv3rxy91izZo38/Pzk7++vSZMmlZ6PjY1Vv379FBAQoOnTp1fpie+LL75YJvBb\nP7F/6623ZDabZTabFR0dXe7aW5/4L1q0SH/+859LX+vZs6ciIyPl6empJ598Ul999ZUGDhwoDw8P\nJSQkVKvm+Ph4PfHEE+XOX7lyRcXFxbp8+bKkX0L4oEGD9Morr8jJyalc+9GjRys+Pv6Ofy8AYARC\nOADUomPHjmnWrFlKSUmRs7OzNm3aJOk/T5gXLlyo7t27KykpSQsXLixz7ZEjR/Tqq6/q22+/VXJy\ncmk4Tk1N1fr167Vnzx4lJSWpUaNGiouLu2MtEyZM0Pr160u/37Bhg0JDQ5WYmKiYmBjt27dPe/fu\n1fLly3Xw4MHb3uvWp9THjx/X3LlzlZaWpqNHj2rdunXavXu3Fi1apNdee01paWlVqrmoqEgpKSny\n8PAo91pcXJxCQkLUrFmz0l9mbsff31979uy5YzsAMEJjowsAgPqssukaJee7deumhx56SJLUu3dv\nnTp1qky7iuY9l/jnP/+pkJAQtW7dWpLk4uIiSfrmm2+UmJioPn36SJIKCgrUrl27O9bq7++vs2fP\nKisrS2fPnpWLi4s6duyojRs3aty4cWrWrJkkady4cdq1a5f8/PzueM8S3bp1k4+PjyTJx8dHw4cP\nlyT5+voqIyOjyjWfP3++dF77reebNGmiZs2aqU2bNsrJySk3BeVWTk5OKi4u1rVr19S0adMqvxcA\nsAVCOADcA1dXV124cKHMudzcXD3wwAOSVGaahIODQ5We4JYwmUyVhvSIiAi99tpr1a53/Pjx2rhx\no7KzsxUaGlrhz7FYLOV+uWjcuHGZ6SMFBQVlXrd+n40aNVKTJk1Kv75586YsFkuVa67oPb///vtq\n0aKFli1bpmvXruncuXN3DOGVvRcAqAuYjgIA96Bly5Zq3769vv32W0lSXl6evvzySz3yyCPlwqTF\nYil3rlWrVsrPz6/w3kOHDtWGDRuUl5dXem9JGjZsmDZu3Khz586Vns/MzKxSvRMmTFB8fLw2btyo\n8ePHS1LpBxgLCgp05coVbd26VYMGDSpznbu7u86ePau8vDwVFhbq008/rdLPK1HVmtu0aVM657tE\nZmamPDw8NGfOHE2bNk2BgYHKycm5488sLCyUg4NDhfPFAcBohHAAuEdr1qzRX/7yFwUEBGjYsGGa\nP3++unXrJqnsdBWTyVT6fcmfrq6uGjhwoMxmc7kPZnp7e+uPf/yjBg8eLH9/f82dO1eS5OXlpVde\neUXBwcHy8/NTcHCwsrOzy9y3sqe/3t7eunz5sjp16iR3d3dJUkBAgJ566in17dtX/fv3V1RUVOlU\nlJL7ODo66k9/+pP69u2r4OBgeXt7l3tv1m597XY1W3NwcJCvr6+OHj0q6ZcpOSNHjpSrq6skKTk5\nWampqVq3bl1poK9MUlKSBgwYcNs2AGAUk+V2ExIBALCxmJgY5eTkVLhaTHW89NJLCgwM1NixY2uo\nMgCoOYRwAECdcv36dQ0fPlw7d+686/nchYWFCgoKuqd7AEBtIoQDAAAANsaccAAAAMDGCOEAAACA\njRHCAQAAABsjhAMAAAA2RggHAAAAbIwQDgAAANgYIRwAAACwMUI4AAAAYGOEcAAAAMDG/g9tNvpM\ny40r9wAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x40cd890>"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Elastic constants\n",
      "-----------------\n",
      "\n",
      "`Elastic` module provides mainly one central function: the procedure for calculation of elastic constants of a crystal. Since the calculation of elastic constants involves many individual calculations (from 10 to even hundreds) it is beneficial to have some control over the execution of individual calculations. Thus the `get_elastic_tensor` provides two modes of operation. \n",
      "\n",
      "* Fully automatic (default) - does not provide any access to the component calculations. In case of any problems with the execution of individual jobs we have no chance to correct the situation. This may be *very* wastefull since the only option is to repeat the whole calculation. But the usage is extremely simple - just call `cryst.get_elastic_tensor()` to get the elastic tensor and Birch coefficients for the given structure.\n",
      "\n",
      "* Staged mode - recommended, particularly for large systems with many individual jobs. This mode provides a three-step procedure:\n",
      "    1. Run `cryst.get_elastic_tensor(mode='deformations')` to get a list of deformed crystals to calculate separately.\n",
      "    2. Run `elastic.ParCalculate(system_list,calculator)` to calculate the stresses of individual structures\n",
      "    3. Pass the list of systems to the `cryst.get_elastic_tensor(mode='restart',systems=system_list)` to get the elastic tensor.\n",
      "\n",
      "We will run the calculation in the staged mode which is safer and gives us a chance to intervene if something goes wrong."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "deformations=cryst.get_elastic_tensor(mode='deformations')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "elastic.ParCalculate(deformations,qe);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Launching:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 2"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 3"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 4"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 5"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 6"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 7"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 8"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 9"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 10"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "     Done:"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 2"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 3"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 4"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 5"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 6"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 7"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 8"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 9"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 10"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Cij,Bij=cryst.get_elastic_tensor(mode='restart',systems=deformations)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Finally we got the elastic tensor of the crystal. Depending on the symmetry of the crystal the number of constants changes.\n",
      "The general ordering of $C_{ij}$ components is (except for triclinic symmetry and taking into account customary names of constants - e.g. $C_{16} \\rightarrow C_{14}$:\n",
      "\n",
      "$$\n",
      "   C_{11}, C_{22}, C_{33}, C_{12}, C_{13}, C_{23}, \n",
      "   C_{44}, C_{55}, C_{66}, C_{16}, C_{26}, C_{36}, C_{45}\n",
      "$$\n",
      "\n",
      "To get the ordering for your crystal drop from the above list the positions irrelevant for its symmetry. Thus for our cubic crystal the order of constants will be: $ C_{11}, C_{12}, C_{44}$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"Elastic tensor (GPa) [C_11, C_12, C_44]:\", Cij/ase.units.GPa"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Elastic tensor (GPa) [C_11, C_12, C_44]: [ 392.951       145.439       274.90982353]\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    }
   ],
   "metadata": {}
  }
 ]
}