{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "matplotlib.use(\"Pdf\")\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import warnings \n",
    "warnings.filterwarnings(\"ignore\") #ignore some matplotlib warnings\n",
    "\n",
    "# numpy\n",
    "import numpy as np\n",
    "\n",
    "from h5 import HDFArchive"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this basic example we will perform a VASP calculation for SrVO$_3$, build PLOs for the Vanadium t$_{2g}$ orbitals, and load them as SumK object, which can be then used to perform a DMFT calculation. \n",
    "\n",
    "__Note: This example works with VASP version 6.5.0 or newer with hdf5 support enabled__\n",
    "\n",
    "## VASP setup\n",
    "\n",
    "First we setup the VASP [INCAR link](INCAR) file by specifing the LOCPROJ, EMIN, EMAX and LORBIT flags:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "SYSTEM = SrVO3\n",
    "NCORE = 4\n",
    "LMAXMIX=6\n",
    "EDIFF = 1.E-10\n",
    "\n",
    "# DOS energy window\n",
    "NEDOS = 2001\n",
    "\n",
    "# Smearing procedure\n",
    "ISMEAR = -5\n",
    "\n",
    "# the energy window to optimize projector channels\n",
    "EMIN = 3.9\n",
    "EMAX = 7.1\n",
    "\n",
    "# use the PAW channel optimization\n",
    "LORBIT=14\n",
    "\n",
    "# project to V d\n",
    "LOCPROJ = 2 : d : Pr\n",
    "```\n",
    "Moreover we prepare a [KPOINTS link](KPOINTS), [POSCAR link](POSCAR), and a POTCAR file. For the POTCAR file please use the VASP provided PBE pseudopotentials: `Sr_sv`, `V`, and `O`. \n",
    "\n",
    "Now VASP is executed, which should converge in roughly 27 iterations. Afterwards you should find the files LOCPROJ and PROJCAR in you directory. \n",
    "\n",
    "## PLOVASP\n",
    "\n",
    "First import the PLOVASP module of DFTTools:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: could not identify MPI environment!\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Starting serial run at: 2025-02-21 14:05:20.237950\n"
     ]
    }
   ],
   "source": [
    "# import plovasp converter\n",
    "import triqs_dft_tools.converters.plovasp.converter as plo_converter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, create a configuration file for plovasp [plo.cfg link](plo.cfg):\n",
    "\n",
    "```\n",
    "[General]\n",
    "DOSMESH = -3.0 3.0 2001\n",
    "\n",
    "[Shell 1]\n",
    "LSHELL = 2\n",
    "IONS = 2\n",
    "EWINDOW = -1.4 2.0\n",
    "\n",
    "TRANSFORM = 1.0  0.0  0.0  0.0  0.0\n",
    "            0.0  1.0  0.0  0.0  0.0\n",
    "            0.0  0.0  0.0  1.0  0.0\n",
    "\n",
    "```\n",
    "where the energy window of the t$_{2g}$ bands is specified by `EWINDOW` and the `TRANSFORM` flag picks the correct three orbitals out of the five Vanadium $d$ orbitals [see the guide for the ordering of orbitals](../../guide/conv_vasp.html). \n",
    "\n",
    "Now run PLOVASP:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[WARNING]: Reading from vaspout.h5\n",
      "Read parameters: LOCPROJ\n",
      "0  ->  {'label': 'dxy       ', 'isite': 2, 'coord': array([-0.5, -0.5, -0.5]), 'l': 2, 'm': 0}\n",
      "1  ->  {'label': 'dyz       ', 'isite': 2, 'coord': array([-0.5, -0.5, -0.5]), 'l': 2, 'm': 1}\n",
      "2  ->  {'label': 'dz2       ', 'isite': 2, 'coord': array([-0.5, -0.5, -0.5]), 'l': 2, 'm': 2}\n",
      "3  ->  {'label': 'dxz       ', 'isite': 2, 'coord': array([-0.5, -0.5, -0.5]), 'l': 2, 'm': 3}\n",
      "4  ->  {'label': 'dx2-y2    ', 'isite': 2, 'coord': array([-0.5, -0.5, -0.5]), 'l': 2, 'm': 4}\n",
      "  Total number of ions: 5\n",
      "  Number of types: 3\n",
      "  Number of ions for each type: [1 1 3]\n",
      "\n",
      "   Reduced number of k-points: 120\n",
      "    Total number of k-points: 3375\n",
      "   Total number of tetrahedra: 1667\n",
      "eigvals from EIGENVAL\n",
      "\n",
      "  Unorthonormalized density matrices and overlaps:\n",
      "  Spin: 1\n",
      "  Site: 2\n",
      "  Density matrix                                                  Overlap\n",
      "   0.6117996  -0.0000018  -0.0000000  -0.0000018  -0.0000007       0.9363450   0.0000000   0.0000000  -0.0000000  -0.0000035\n",
      "  -0.0000018   0.6117247   0.0000000  -0.0000018  -0.0000000       0.0000000   0.9363441   0.0000000  -0.0000000   0.0000000\n",
      "  -0.0000000   0.0000000   0.6323387  -0.0000000   0.0000000       0.0000000   0.0000000   0.6770489   0.0000000   0.0001153\n",
      "  -0.0000018  -0.0000018  -0.0000000   0.6117247  -0.0000000      -0.0000000  -0.0000000   0.0000000   0.9363441  -0.0000000\n",
      "  -0.0000007  -0.0000000   0.0000000  -0.0000000   0.6323382      -0.0000035   0.0000000   0.0001153  -0.0000000   0.6743177\n",
      "\n",
      "  Generating 1 shell...\n",
      "\n",
      "    Shell         : 1\n",
      "    Orbital l     : 2\n",
      "    Number of ions: 1\n",
      "    Dimension     : 3\n",
      "    Correlated    : True\n",
      "    Ion sort      : [2]\n",
      "Density matrix:\n",
      "  Shell 1\n",
      "Site diag : True\n",
      "    Site 1\n",
      "     0.3333529    -0.0000025    -0.0000025\n",
      "    -0.0000025     0.3332565    -0.0000025\n",
      "    -0.0000025    -0.0000025     0.3332565\n",
      "      trace:  0.9998660024816521\n",
      "\n",
      "  Impurity density: 0.9998660024816521\n",
      "\n",
      "Overlap:\n",
      "  Site 1\n",
      "[[ 1. -0.  0.]\n",
      " [-0.  1. -0.]\n",
      " [ 0. -0.  1.]]\n",
      "\n",
      "Local Hamiltonian:\n",
      "  Shell 1\n",
      "    Site 1 (real | complex part)\n",
      "     0.5835197    -0.0000000    -0.0000000 |     0.0000000    -0.0000000     0.0000000\n",
      "    -0.0000000     0.5835196    -0.0000000 |     0.0000000     0.0000000    -0.0000001\n",
      "    -0.0000000    -0.0000000     0.5835196 |    -0.0000000     0.0000001    -0.0000000\n",
      "\n",
      "Evaluating DOS...\n",
      "  Shell 1\n",
      "    Total number of states: [[[1.99988129 1.99954392 1.9994287 ]]]\n",
      "  Storing ctrl-file...\n",
      "  Storing PLO-group file 'vasp.pg1'...\n",
      "  Density within window: 0.9999999999316307\n"
     ]
    }
   ],
   "source": [
    "# Generate and store PLOs\n",
    "plo_converter.generate_and_output_as_text('plo.cfg', vasp_dir='./')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PLOVASP created one shell with three orbitals, which are equally filled by 1/3, one electron in total. Additionally we calculated the density of states. Both in VASP and PLOVASP. The later stores the data in the file pdos_x.dat, which can be simply plotted with [matplotlib](https://matplotlib.org/):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "plo_dos = np.loadtxt('pdos_0_0.dat')\n",
    "\n",
    "with HDFArchive('vaspout.h5','r') as vaspout:\n",
    "    vasp_dos = vaspout['results/electron_dos/dos']\n",
    "    vasp_dos_energies = vaspout['results/electron_dos/energies']\n",
    "    vasp_fermi = vaspout['results/electron_dos/efermi']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmUAAAF8CAYAAACOvszpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABurUlEQVR4nO3dd3gU1f4G8He2b9qmk0IaIZTQe+9ViqKi4AUF4WJDLCgqVqxcr41rQ/QKiIIoXik/UYpIk95Cb4FAQnpPNmXr/P5Y2GRJQgqb7G7yfp5nH2fPzsx+E2L2zZkz5wiiKIogIiIiIoeSOLoAIiIiImIoIyIiInIKDGVEREREToChjIiIiMgJMJQREREROQGGMiIiIiInwFBGRERE5AQYyoiIiIicAEMZERERkRNgKCMiIiJyAg4NZbt27cL48eMREhICQRCwbt0662sGgwEvvvgiOnToAHd3d4SEhOChhx5CSkqK4womIiIiqicODWVFRUXo1KkTvvjiiwqvFRcX4+jRo3jttddw9OhR/Prrrzh//jzuvPNOB1RKREREVL8EZ1mQXBAErF27FhMmTKhyn0OHDqFnz564evUqwsPDG644IiIionomc3QBtZGfnw9BEODt7V3lPjqdDjqdzvrcbDYjJycHfn5+EAShAaokIiKipkoURRQWFiIkJAQSSe0uSLpMKCstLcWLL76IBx54AF5eXlXut3DhQrz55psNWBkRERGRraSkJDRv3rxWx7jE5UuDwYB7770X165dw44dO24Zym7uKcvPz0d4eDiSkpJueRyRK7l48aKjSyCqVmZhKTJ+fxsTpHutbRcnbndgRfYVExPj6BLICRUUFCAsLAx5eXnQaDS1Otbpe8oMBgPuv/9+XL16FX/99Ve1wUqpVEKpVFZo9/LyYiijRsPDw8PRJRBVqwhKXJbHwEu+z9rWmH52+ZlCt1KXIVNOHcpuBLKLFy9i+/bt8PPzc3RJRERERPXCoaFMq9UiPj7e+jwhIQFxcXHw9fVFcHAwJk6ciKNHj+K3336DyWRCWloaAMDX1xcKhcJRZRMRUQ0JcIoRMkQuwaGh7PDhwxgyZIj1+dy5cwEA06ZNw4IFC7BhwwYAQOfOnW2O2759OwYPHtxQZRIRUR0IRh3ule52dBlELsOhoWzw4MG41X0GDXkPgslkgsFgaLD3o6ZFLpdDKpU6ugyiBiUzFKK95IqjyyByGU49pqwhiKKItLQ05OXlOboUauS8vb0RFBTE+fKo6RDNFZtEkf8PEFWhyYeyG4EsMDAQbm5u/GVBdieKIoqLi5GRkQEACA4OdnBFRA1DEE0V2kQA/C1LVLkmHcpMJpM1kPHOTqpParUaAJCRkYHAwEBeyqQmQTSXhbJ1pr54xvAkfmcqI6qSQxckd7QbY8jc3NwcXAk1BTd+zjh2kZoKc7lQZoLlDxHei0lUtSYdym7gJUtqCPw5o6ZGNBut2ybR8nHjHGvIEDknhjIiIqoXorlsoL8JN0IZUxlRVRjKyOUJgoB169Y5ugwiuoloKuspe0C2HVdU/+DlS6JbYChzMePHj8fo0aMrfW337t0QBAEnTpywtj366KOQSqVYs2ZNhf2Li4sxf/58REdHQ6VSISAgAIMGDcL69eut+wwePBiCIEAQBKhUKsTGxuLLL7+ssr4rV65AEATExcXV6uuq63E1MX36dOvXIJfL0axZM4wYMQJLly6F2Vzxlv29e/dizJgx8PHxgUqlQocOHfDxxx/DZLK9k2znzp0YOnQofH194ebmhpiYGEybNg16vd7uXwORKxIrnRLDAYUQuQiGMhczc+ZMbN26FdeuXavw2rJly9C9e3d07NgRgCV0rV69Gi+88AKWLl1aYf/HHnsMv/76Kz777DOcO3cOmzZtwsSJE5GdnW2z36xZs5CamoozZ87g/vvvx+zZs/Hjjz/WzxdYT0aPHo3U1FRcuXIFf/zxB4YMGYKnn34a48aNg9FY9tf82rVrMWjQIDRv3hzbt2/HuXPn8PTTT+Odd97B5MmTrZdezpw5g9GjR6N79+7YtWsXTp48ic8++wwKhaJCeCNqqkxixXGUzGREVWMoczHjxo1DQEAAli9fbtOu1WqxZs0azJw509q2Zs0axMbG4qWXXsKuXbuQlJRkc8yGDRvw8ssvY8yYMYiMjES3bt0wZ84czJgxw2Y/Nzc3BAUFoUWLFliwYAFiYmKsS2DdLCoqCgDQpUsXCIJgXQ7LbDbjrbfeQvPmzaFUKtG5c2ds2rSp2uMOHTqEESNGwN/fHxqNBoMGDcLRo0dr/X1TKpUICgpCaGgounbtipdffhnr16/HH3/8Yf1eFhUVYdasWbjzzjvx9ddfo3PnzoiMjMQ///lPfPfdd/jll1/w888/AwC2bNmCoKAg/Pvf/0b79u0RHR2N0aNH45tvvrFOf0HU1OV6tsLbhik2bWamMqIqMZS5GJlMhoceegjLly+3GTC7Zs0amEwmPPDAA9a2b7/9FlOnToVGo8Edd9xRIcgFBQXh999/R2FhYa1qUKvVVV6iO3jwIADgzz//RGpqKn799VcAwH/+8x989NFH+PDDD3HixAmMGjUKd955Jy5evHjL4woLCzFt2jT8/fff2L9/P2JiYjBmzJha11yZoUOHolOnTtb32rJlC7Kzs/H8889X2Hf8+PFo1aqVtYcwKCgIqamp2LVr123XQdRYmcxiJVOSMZURVaVJTx5bmfGf/Y3MQl2Dv2+ApxL/N6d/jfadMWMGPvjgA+zcudPao7Rs2TLce++90Gg0AICLFy9i//791sAxdepUzJ07F6+++qp1aoavv/4aU6ZMgZ+fHzp16oT+/ftj4sSJ6NevX6XvazKZ8OOPP+LEiRN45JFHKv86AgIAAH5+fggKCrK2f/jhh3jxxRcxefJkAMD777+P7du3Y9GiRfjiiy+qPG7o0KE25//666/h7e2NnTt3Yty4cTX6ft1KmzZtrGPwLly4AABo27Ztlfve2Oe+++7D5s2bMWjQIAQFBaF3794YNmwYHnroIXh5ed12XUSNQWW9YmZ2lRFViT1lN8ks1CGtoLTBH7UJgm3atEHfvn2t48Ti4+Oxe/dum0uXS5cuxahRo+Dv7w8AGDNmDPLz8/HXX39Z9xk4cCAuX76Mbdu2YeLEiTh9+jQGDBiAt99+2+b9vvzyS3h4eECtVmPWrFl49tln8fjjj9e43oKCAqSkpFQIe/369cPZs2dveWx6ejpmzZqFmJgYaDQaeHl5QavVIjExscbvfyuVrcNXk1v2pVIpli1bhmvXruHf//43QkND8d5776Fdu3ZITU21S21Ers4kihBu6hkT2VNGVCX2lN0kwFPpEu87c+ZMzJkzB1988QWWLVuG6OhoDBo0CIClR+u7775DWloaZLKyf2KTyYSlS5di2LBh1ja5XI4BAwZgwIABePHFF/HOO+/grbfewosvvgiFQgEAmDJlCl555RWo1WoEBwdDImm4LD9t2jRkZ2fjP//5DyIiIqBUKtGnTx+73eF49uxZ63i2Vq1aWdv69u1b6b6xsbE2baGhoXjwwQfx4IMP4u2330arVq3w1Vdf4c0337RLfUSuTFMQj1fkq2zaePclUdUYym5S00uIjnb//ffj6aefxqpVq7BixQo8/vjj1h6fG+PEjh07ZrPG4qlTp/Dwww8jLy8P3t7elZ43NjYWRqMRpaWl1lCm0WjQsmXLGtV145jydyB6eXkhJCQEe/bssQZHANizZw969uxZ5XE39vnyyy8xZswYAEBSUhKysrJqVEt1/vrrL5w8eRLPPvssAGDkyJHw9fXFRx99VCGUbdiwARcvXqzQi1iej48PgoODUVRUZJf6iFydSpdRoY09ZURVYyhzUR4eHpg0aRLmz5+PgoICTJ8+3frat99+i7Fjx6JTp042x8TGxuLZZ5/FypUrMXv2bAwePBgPPPAAunfvDj8/P5w5cwYvv/wyhgwZUudxUYGBgVCr1di0aROaN28OlUoFjUaDefPm4Y033kB0dDQ6d+6MZcuWIS4uDitXrrzlcTExMfj+++/RvXt3FBQUYN68eXW6u1Gn0yEtLQ0mkwnp6enYtGkTFi5ciHHjxuGhhx4CALi7u2PJkiWYPHkyHnnkETz55JPw8vLCtm3bMG/ePEycOBH3338/AGDJkiWIi4vD3XffjejoaJSWlmLFihU4ffo0Pvvsszp974gaHXPF6WHESuYGJCILjilzYTNnzkRubi5GjRqFkJAQAJYxWBs3bsS9995bYX+JRIK7774b3377LQBg1KhR+O677zBy5Ei0bdsWc+bMwahRo6zTPtSFTCbDp59+iiVLliAkJAR33XUXAOCpp57C3Llz8dxzz6FDhw7YtGkTNmzYgJiYmFse9+233yI3Nxddu3bFgw8+iKeeegqBgYG1rmvTpk0IDg5GZGQkRo8eje3bt+PTTz/F+vXrbXoTJ06ciO3btyMxMREDBgxA69at8cknn+CVV17B6tWrrb2RPXv2hFarxWOPPYZ27dph0KBB2L9/P9atW2fTG0jUpFUWythRRlQlQWzkC5EVFBRAo9EgPz+/Qu9PaWkpEhISEBUVBZVK5aAKqamw58/b+fPn7VQVUf1JOrQBw+PfAQC8bZiKb01jsGpSJPzcGsdFmtatWzu6BHJCt8od1WFPGRER1Y9ylyqNsPRIN+5uAKLbw1BGRET1QjSXLWFmuv5xY2YqI6oSQxkREdWPcguSm/lxQ1Qt/l9CRET1Qywb6P+e/FtcUf0DorHhV0whchUMZUREVD/Eindf3jzDPxGVYSgjIqJ6kS0PrdDGtS+JqsZQRkRE9SLBowveN0y+qZWTxxJVhaGMiIjqhclcyYLk7CgjqhJDGRER1QuzCPxl7mLT1sjnKye6LQxlRERUL8yiiHNiuE0bMxlR1RjKXNT06dMhCAIEQYBCoUDLli3x1ltvwWg0Wl+fMGHCLc+RlJSEGTNmICQkBAqFAhEREXj66aeRnZ1d7fs/++yzuOeee+zxpRBRI9U76xccVT5yUytTGVFVGMpc2OjRo5GamoqLFy/iueeew4IFC/DBBx/U6NjLly+je/fuuHjxIn788UfEx8fjq6++wrZt29CnTx/k5OTc8viDBw+ie/fu9vgyiKiRkptK4CtobdpEkQP9iarCUObClEolgoKCEBERgccffxzDhw/Hhg0banTs7NmzoVAosGXLFgwaNAjh4eG444478OeffyI5ORmvvPJKpcfp9XrI5XLs3bsXr7zyCgRBQO/eve35ZRFRY1FJAOPlS6KqMZQ1Imq1Gnq9vtr9cnJysHnzZjzxxBNQq9U2rwUFBWHKlCn46aefKh2QK5PJsGfPHgBAXFwcUlNTsWnTJvt8AUTUqAhi2dqXW01d8bZhKowShQMrInJuDGWNgCiK+PPPP7F582YMHTq02v0vXrwIURTRtm3bSl9v27YtcnNzkZmZWeE1iUSClJQU+Pn5oVOnTggKCkJhYSEGDx6M2NhYdOzYEWvWrLntr4mIXJ9Qrqfsa+M4fGsaA7NE6cCKiJybzNEFOKW9nwP7vqh+v+BOwD9W27atmgykHq/+2D6zgb5P1q2+63777Td4eHjAYDDAbDbjH//4BxYsWFDj4+t6a/qxY8fQqVMn63OZTIZFixahc+fOSEtLQ7du3TBmzBi4u7vX6fxE1DgI5ZZZMl3vAzDz+iVRlRjKKqMrBApTqt9PU3EJERRn1exYXWHt67rJkCFDsHjxYigUCoSEhEAmq9k/Z8uWLSEIAs6ePYu77767wutnz56Fj48PAgICKj0+Li7OJpQFBwcjODgYgOXyp7+/P3JycqyhLC4uDo8//jiKi4vxj3/8A3/99Rc2b95c2y+XiFxNuZ4yEy/MEFWLoawySk/AM6T6/dz8K2+rybFKz9rXdRN3d3e0bNmy1sf5+flhxIgR+PLLL/Hss8/ajCtLS0vDypUr8dBDD0EQhEqPP3nyJO69995KXzty5AhMJhPCwsIAAAaDAdOnT8fq1avRpk0b3HnnnejYsWOtayYi1yMp11MWKOQhRMyC2RzswIqInBtDWWX6Pln3S4s3X850oPz8fMTFxdm0+fn5ISwsDJ9//jn69u2LUaNG4Z133kFUVBROnz6NefPmITQ0FO+++26V5zWbzTh//jxSUlLg7u4OjUYDwHIDwUMPPYRvvvnGuu/atWvRp08ftGnTBoBlvFr79u3t/8USkdMpf/nyG8XHAIA/DOsBcGgDUWXYn9yI7dixA126dLF5vPnmmwCAmJgYHD58GC1atMD999+P6OhoPPLIIxgyZAj27dsHX1/fKs/7zjvvYPny5QgNDcU777wDANDpdJgwYQJeeukl9O3b17rviRMn0LlzZ+vz06dPs6eMqIkoH8qszJynjKgq7ClzUcuXL6/29er2iYiIqHafykydOhVTp061PhdFEdOnT8fQoUPx4IMP2uzr6+uL+Ph4AJaQuG3bNvz666+1fk8icj273Ebhl9wYLFJ8WdZYbpoMIrLFnjK6bXv27MFPP/2EdevWoXPnzujcuTNOnjwJwBLgtm/fjk6dOmHDhg3o0aMHFArOU0TUFMTLW2GduT/+z1Q2wbRgrqT3jIgAODiU7dq1C+PHj0dISAgEQcC6detsXhdFEa+//jqCg4OhVqsxfPhwXLx40THFUpX69+8Ps9mMuLg466NDhw4ALDcjHD58GMeOHYNUKq3Qk0ZEjZf5+uwXRkjLGrnMElGVHBrKioqK0KlTJ3zxReVzgv373//Gp59+iq+++goHDhyAu7s7Ro0ahdLS0gaulOrqgw8+QPv27dG1a1coFAr885//dHRJRNRATNfnJDPZhDL2lBFVxaFjyu644w7ccccdlb4miiIWLVqEV199FXfddRcAYMWKFWjWrBnWrVuHyZMnN2SpVEcLFiyo1YS2RNR4BOhT0V7IQishydpW6eB/IgLgxGPKEhISkJaWhuHDh1vbNBoNevXqhX379jmwMiIiqokHtEvxm/JVdJQklDWaOdCfqCpOe/dlWloaAKBZs2Y27c2aNbO+VhmdTgedTmd9XlBQUD8FEhHRLUkqGT/Ggf5EVXPanrK6WrhwITQajfVxY2Z5IiJqWAIqGdTPy5dEVXLaUBYUFAQASE9Pt2lPT0+3vlaZ+fPnIz8/3/pISkqqcl8iIqo/5ZdZulv3JobqPkSuW7QDKyJybk4byqKiohAUFIRt27ZZ2woKCnDgwAH06dOnyuOUSiW8vLxsHkRE1PAkKAtlZ8VwXBZDYJRynkKiqjh0TJlWq7XO9g5YBvfHxcXB19cX4eHheOaZZ/DOO+8gJiYGUVFReO211xASEoIJEybYtQ7x+m3bRPWJP2fU1AjlxpTdmBaD/xsQVc2hoezw4cMYMmSI9fncuXMBANOmTcPy5cvxwgsvoKioCI888gjy8vLQv39/bNq0CSqVyi7vL5fLAQDFxcVQq9V2OSdRVYqLiwGU/dwRNXblL1+arl+YYSgjqppDQ9ngwYNv2XsgCALeeustvPXWW/Xy/lKpFN7e3sjIyAAAuLm5QRCEenkvarpEUURxcTEyMjLg7e0NqVRa/UFEjYCk3ED/+6Q7IYUZbqXDAXBcGVFlnHZKjIZy46aBG8GMqL54e3vf8iYVosamfCh7X/4NAOAvbRQYyogq1+RDmSAICA4ORmBgIAwGg6PLoUZKLpezh4yaHGml61xySgyiqjT5UHaDVCrlhyYRkR09q1qA5NwSPCLbiMdk/weAk8cS3QpDGRER1Qst3JEDOVJFX2sb174kqprTzlNGRESu7cZ9XKbyHzUMZURVYigjIqJ6cWNEWflQJnBBcqIq8fIlERHVi7uNf8AgLcRY6YGyxkoH/xMRwFBGRET1ZLJxA4LlmbaNHOhPVCVeviQionpRfp6yG3j5kqhq7CkjIqJ6lyr6wgQJ9BL7LJNH1BgxlBERUb240VOWIvqir+5zAMDTfgEIc2RRRE6Mly+JiKheiShbU9jMBcmJqsRQRkRE9UKAJYHZhjKmMqKqMJQREVG9sIYykT1lRDXBMWVERFQvboSyMEkmPpf/B1KYYcgZAmCiYwsjclIMZUREVC8kKOsWG3d9Atk9JVGOKofI6TGUERFRvbiACLiZtQgRsuAraAEAEq59SVQlhjIiIqoXjwqvIVdvQmchHuuUrwMABIYyoipxoD8REdWLG4P6jeUXJGcoI6oSQxkREdUL8fr0F2aGMqIaYSgjIqJ6cWOYvxFSaxtDGVHVOKaMiIjqxRK8C6nCAGm5hckZyoiqxlBGRET1oivOQSUxoFhUWtsYyoiqxsuXRERUL27M428od/lSIpor35mI2FNGRET1Q7h+2bIQbths7AEjJJCq2qOzY8sicloMZUREVC9u9JTlih54wfgoAOAuTw1DGVEVePmSiIjqxY21L8tPicEFyYmqxlBGRET1ovzalzeYRKYyoqowlBERUb2QCJYAJrKnjKhGOKaMiIjsr1yPmEIw4IRyJqQwIyGjHYCvHFcXkRNjKCMionpQFspMkMFLKAEAKMy6Si5qEhHAUEZERPXkC+OdEADkSv3RAZcAABLRBE4fS1Q5hjIiIrI7EQI+ME4GAPjKpXjJvAxSQYQEDGVEVeFAfyIisrvyA/rlUsG6KDmXWSKqGkMZERHVK7VMAtP1UCZhKCOqEi9fEhGR3ZnNIpTQAwA8ZDKYrvcBsKeMqGoMZUREZHeCsRjnVdMBAMdLO1hDGXvKiKrGy5dERGR3Yrl5yizD+3n5kqg6DGVERGR3N6+mZBauD/SH2QHVELkGXr4kIiK7E8Wy8GUWJPhAOgsFJXqIcjXmOrAuImfGUEZERHYn3tRVdljZE5eL9JAbBYYyoio49eVLk8mE1157DVFRUVCr1YiOjsbbb79d4X92IiJyLqK5rKdMhARqueXjxmASoTfxdzhRZZy6p+z999/H4sWL8d1336Fdu3Y4fPgwHn74YWg0Gjz11FOOLo+IiKog2qxwKcBTUdYHoNWZ4Ovm1B8/RA7h1P9X7N27F3fddRfGjh0LAIiMjMSPP/6IgwcPOrgyIiK6pXJjykRBQHNZProIV6ARtCgt8gHc/BxYHJFzcurLl3379sW2bdtw4cIFAMDx48fx999/44477qjyGJ1Oh4KCApsHERE1rPLLLIkQMK54HdYq38ByxQcQci45rjAiJ+bUPWUvvfQSCgoK0KZNG0ilUphMJrz77ruYMmVKlccsXLgQb775ZgNWSUREFZSbj0yEAJPCy/rcVJLviIqInJ5Th7Kff/4ZK1euxKpVq9CuXTvExcXhmWeeQUhICKZNm1bpMfPnz8fcuWX39hQUFCAsLKyhSiYiIgAGmRcm6V6DIIiI8A3ACGW89TWRoYyoUk4dyubNm4eXXnoJkydPBgB06NABV69excKFC6sMZUqlEkqlsiHLJCKim5ilChwQ2wIiIJW5A6qMshd1DGVElXHqMWXFxcWQSGxLlEqlMJs5IzQRkTMrP3ORIAASlbf1uUTPsb5ElXHqnrLx48fj3XffRXh4ONq1a4djx47h448/xowZMxxdGhER3YK5XCoTBAFSN431uUxf6IiSiJyeU4eyzz77DK+99hqeeOIJZGRkICQkBI8++ihef/11R5dGRES3IBhLMVJyCCIEBBuaQ+EWan1NbmBPGVFlnDqUeXp6YtGiRVi0aJGjSyEiolqQl2bja8UnAIB9RQMhd3vZ+prKpHVUWUROzanHlBERkWvSZiVat01SFeRKtfW5XCx1RElETo+hjIiI7E6Xedm6bQjsCLlcCbMoAACkZqOjyiJyak59+ZKIiFyUodi6KXfTQCqVoJvxWxSY5Aj3UeMrB5ZG5KzYU0ZERHYnMZZYt2VKdwCAUeYGE6QoNYpVHUbUpDGUERGR3UmNZT1lMqUbAEAls3zk6Iyca5KoMgxlRERkd1JTWU+ZXG3pKVPKLGPK2FNGVDmOKSMiIruTlwtlCpUllN1r/hNqWTKUohHAWw6qjMh5MZQREZHdhRotU2IYRCkkbn4AgJHmnWgnOwcAOG18DTKZ3GH1ETkjhjIiIrKrwlI98oz+8JbkIlsaCMgtc5QZBaV1H31pCWQeDGVE5TGUERGRXZ3LMuBVwzwAIv7RWopp19uNEoV1H6OhFICXI8ojcloc6E9ERHZ1JuPGjP0CIoICrO1Gicq6bdCVgIhsMZQREZFdnc0sW0YpNrAsiJmkZZcvDXoutUR0M4YyIiKyG5PJjPgMy4Lj/m5SBLiXjZIxlwtlZj17yohuxjFlRERkN5mpV3BA+k8cF6Jx3HMEBCHK+lr5UGYysKeM6Ga16imbOHEiNm3aBFHkxH9ERFRRceJRuAk69JGeQVtVns1rZmnZpUwzL18SVVCrUJabm4uxY8ciPDwcr7/+Oi5fvlxfdRERkQtSZ520bktDOtu+KCt3+ZI9ZUQV1CqUbdu2DZcvX8bMmTPxww8/ICYmBkOHDsWqVaug0+nqq0YiInIRocVnAAAmUYBPZGeb14pUIdhnisV2UydoBU8HVEfk3Go90D8iIgILFizA5cuXsXXrVoSEhGDWrFkIDg7G7NmzceTIkfqok4iInFxhYQGizZaZ/BOkkZBdX17phquBQ/GA4VU8bHgRl9XtHFEikVO7rbsvhw4dih9++AFpaWlYuHAhVq9ejV69etmrNiIiciE5CccgESxjjtM8KoYu1fUFyQEuSk5Umdu++zIhIQHLly/H8uXLkZ+fj+HDh9ujLqew91IWnv/5OAa3CcR7d3dwdDlERE5NSI2zbhsCO1V4XVkulOlM5oYoicil1KmnrLS0FD/88AOGDh2KmJgYrFixAjNnzkRCQgI2bdpk7xod5h/fHEBKfilWHUhEUk6xo8shInJq/gWnrNteUV0rvK6Uln3klBrYU0Z0s1r1lB08eBBLly7FTz/9hNLSUtx9993YtGkThg0bBkEQqj+BCyvWmxxdAhGR0zIZjYgxXAAEIBM+8PBrXmGfwOKL2Kh4GyrocTp7BIAnG75QIidWq1DWu3dvdOrUCW+//TamTJkCHx+f+qrL6Ugad+YkIrotmdcuIFawzNJ/RdUWXpX8oS6XydBOchUAkGxIa9D6iFxBrULZ4cOH0bVrxS7ppkBn5PgHIqKqHCwJwUuln6Gb5AL6RYehsk8KiVpj3VYZCxuuOCIXUasxZeUD2e7duzF16lT06dMHycnJAIDvv/8ef//9t30rdBJ6DkolIqrSmYxSpMIPv5n7wKNF70r3UXmUXV1RmxjKiG5Wp4H+//vf/zBq1Cio1WocO3bMOnFsfn4+3nvvPbsW6Cz07CkjIqrSxWzL54BCKqCFr7LSfZRKFXSiHADgxlBGVEGdQtk777yDr776Ct988w3kcrm1vV+/fjh69KjdinMkk9n2ziCGMiKiyulNIlIKDACAcI0CsioG4QoSCQoEDwCAh6htsPqIXEWd5ik7f/48Bg4cWKFdo9EgLy/vdmtyCgUlBpvnDGVERJXLSz6PF6U/4ZQ5CnLPbrfct1DwQICYC09Ri7yGKY/IZdQplAUFBSE+Ph6RkZE27X///TdatGhhj7ocLv/mUMYxZURElUs+jEdlGwEAfwiPA6h6CaViiSdgAtwEHYz6UsgUqgYqksj51eny5axZs/D000/jwIEDEAQBKSkpWLlyJZ5//nk8/vjj9q7RIW4OZToj5ykjIqqMMu+iddvs3/aW+5ZIPazbpUV59VUSkUuqU0/ZSy+9BLPZjGHDhqG4uBgDBw6EUqnE888/jzlz5ti7Rrswm2s3e3SFnjJeviQiqpSm6Kp12yO41S331cm8AP31bW0ePHyC6rM0IpdSp1AmCAJeeeUVzJs3D/Hx8dBqtYiNjYWHh0f1BztIn39tw9hu0fhgYscarT7AUEZEVAOiiGCDJZSlib4I8L31pOJnvAZgR34z5MMdgwRv+DVEjUQuolaXL8PDw/Hkk09iy5YtMBqNUCgUiI2NRc+ePZ06kAFAkc6EX45cQ2IN17CsePmSoYyI6Gbm4mx4oQgAcE3aHNJqlj9J9umF/5rGYo1pMLLNng1RIpHLqFUo+/7776FUKjF79mz4+/tj0qRJWLlypUvdcXn4Sm6N9uNAfyKi6mnT4q3bOeqIavf3VEqt24U6jtUlKq9WoWzQoEH46KOPcPHiRezZswedO3fGZ599hqCgIAwdOhSLFi3C5cuX66tWu7iaXVSj/W6eEkNnYCgjIrqZPqMslJV6RlW7v6ey7GOnUM/fq0Tl1enuSwBo164d5s+fj/379yMhIQGTJ0/Gtm3b0L59e7Rv3x4bN260Z512c7WOly/ZU0ZEVJEsr+wPcYlfdLX7eylEBCAXMcI1SLUp9Vkakcup00D/mwUHB+ORRx7BI488gqKiImzZsgVKZeXLbDhallZXo/1yi/U2zznQn4ioIs/isjsvPYNbVrt/qD4Bh1SzAQC7M8cAWFBPlRG5njqFsqNHj0Iul6NDhw4AgPXr12PZsmWIjY3FggULcPfdd9u1SHsq1lc/hiEppxibT6fbtDGUERFVdMTUCjkmM3yEIgT4+Ve7v8K97O5MhaGgPksjcjl1unz56KOP4sKFCwCAy5cvY/LkyXBzc8OaNWvwwgsv2LVAeyupQSh7d+PZCm0MZUREtoxmEW+W3Iuphlcw2+3Dau+8BAClu7d1W8VFyYls1CmUXbhwAZ07dwYArFmzBgMHDsSqVauwfPly/O9//7NnfXbjrbZ0ClbXU6Y3mrHzQmbFdo4pIyKykVJgwI2/V8O9FTU6RuXmAaNo+ehxYygjslGnUCaKIsxmy/+Jf/75J8aMGQMACAsLQ1ZWlv2qA5CcnIypU6fCz88ParUaHTp0wOHDh2t9HrXCEspKDLcOZRfSC637dGqusbZzmSUiIluJeWVjbyNqGMokEgnyYZnX0t3MUEZUXp1CWffu3fHOO+/g+++/x86dOzF27FgAQEJCApo1a2a34nJzc9GvXz/I5XL88ccfOHPmDD766CP4+Nx6xujKqOSWL7W6y5d/nEq1bveJLhsfwcuXRES2knO01u0In5qFMgDQCpZQ5ilqq9mTqGmp00D/RYsWYcqUKVi3bh1eeeUVtGxpuePml19+Qd++fe1W3Pvvv4+wsDAsW7bM2hYVVf08OJVRy2UAjCjWGyGKYqVLLZnNIr7Yfsn6vHtEWfjjjP5ERLbGXnkPTypP4bIYgkK3TwDUbGWXIoknYAY8hRKYjAZIZfL6LZTIRdQplHXs2BEnT56s0P7BBx9AKpVWckTdbNiwAaNGjcJ9992HnTt3IjQ0FE888QRmzZpV5TE6nQ46Xdm0FwUFlrt71AopACPMoiVgqeQV67x5Goxu5UIZe8qIiGz565PgK2jhjgRcrMUVjBKpJZQBgK44H25e1d+1SdQU1Hny2Ly8PPz3v//F/PnzkZOTAwA4c+YMMjIy7Fbc5cuXsXjxYsTExGDz5s14/PHH8dRTT+G7776r8piFCxdCo9FYH2FhYQAAT1VZ/rx5tv4bsovKQlnbYC9o1GV/vXGgPxFRGZPRgGBzGgDgmiQYcnnN/8YvlZWteanT5tm7NCKXVaeeshMnTmDYsGHw9vbGlStXMGvWLPj6+uLXX39FYmIiVqxYYZfizGYzunfvjvfeew8A0KVLF5w6dQpfffUVpk2bVukx8+fPx9y5c63PCwoKEBYWBn8PJXB90dyMQh0CvVQVjs0qLOthG9QqABKJALlUgMEkcpklIqJy8tOvQiFYxuhmKcPgXYtjDXIvoMSyrS/Ks3dpRC6rTj1lc+fOxcMPP4yLFy9CpSoLN2PGjMGuXbvsVlxwcDBiY2Nt2tq2bYvExMQqj1EqlfDy8rJ5ALgeyiwyCyuf1T+z3Gz//h6WQasKqeVbxJ4yIqIyRRll42+L3cNrdez+wPsxUPcJOpZ+jURVa3uXRuSy6tRTdujQISxZsqRCe2hoKNLS0m67qBv69euH8+fP27RduHABERERtT5XgFfZnUGJVax/maUtu3x5I8Qp5VIU6U0cU0ZEVI6Ye8W6bfaOrN3BbgFIFC03WxXqq9mXqAmpU0+ZUqm0DqAv78KFCwgICLjtom549tlnsX//frz33nuIj4/HqlWr8PXXX2P27Nm1PldskJd1+2BCTqX7ZNv0lFlCmbWnjKGMiMhKrS27YqEIqH4h8vI8lWUfPYU6zgFJdEOdQtmdd96Jt956CwaDZcC8IAhITEzEiy++iHvvvdduxfXo0QNr167Fjz/+iPbt2+Ptt9+2TsdRW22DveCptHQM7r+cDVEUK+xTfrFyf8/rly9lvHxJRHQz39Ik67ZnYGStjvVUlt39Xqjn71aiG+oUyj766CNotVoEBgaipKQEgwYNQsuWLeHp6Yl3333XrgWOGzcOJ0+eRGlpKc6ePXvL6TBuRSaVoGeULwDLXZYXMypOWlj+8qWf+/WeMhl7yoiIbuZlygUA5IkeULt7VbO3LV+JFlOlWzFbug6RObvrozwil1SnMWUajQZbt27Fnj17cPz4cWi1WnTt2hXDhw+3d3121SfaD9vOWabs2HE+A62aedq8fuPypUQAfN1tB/pzmSUiIgu9SUSwWAAIQKGkdoEMAHyEIrwjt0wKvr9gEIC77FwhkWuqU0/ZihUroNPp0K9fPzzxxBN44YUXMHz4cOj1ertNh1EfhrUtWwLq/46nVnj9Rk+Zr7sCUollEKry+vJMBpMIs7niJU8ioqYmU2vAdP2L+Kf+OfziPaPWx5efLNZTb7+5LYlcXZ1C2cMPP4z8/PwK7YWFhXj44Ydvu6j6EuXvjg6hlkXGTybn41hirvU1URStU2KUnz7jRk8ZwHFlREQAUGQw46jYCn+au+Gypk+tj3fz8EI6LMNJmhsTIZr5u5UIqGMoq2rtyGvXrkGj0dx2UfVpSq+y+XQe+GY/Xl9/CtvOpiOnSG8dNxbgWS6UyRjKiIjKKyk3mbZaXvGzoCZS5JEAAI1QhPxs+02lROTKajWmrEuXLhAEAYIgYNiwYZDJyg43mUxISEjA6NGj7V6kPd3TtTk+3x6Pa7klKDWYsWLfVfx0KAmLp3a17hPl727dVpYPZRzsT0SEEmPZUA61rG6r9RW4RwF5RwEAhann4R0QYpfaiFxZrULZhAkTAABxcXEYNWoUPDw8rK8pFApERkbadUqM+qCQSfDu3R0wbelBa5vOaAlnN5QPZeV7ynQMZUREkBdewyjJUWSLnvAROgDXL0XWhtG3JZBn2TZnXQQwxJ4lErmkWoWyN954AwAQGRmJSZMm2Syx5EoGtQrAW3e1w+vrT1vbdpzPtG53i/CxbitlZfPpsKeMiAgIyj2EJYpPAQBbip4HUPtVVtTNYoDLlm23gku33pmoiahTv/O0adNcNpDd8FCfSFx89w60vmlaDE+lDLHBZbd42wz0ZygjIoJMl2fdNim963QOTXAMTNeXWgosvXL7RRE1AnUKZSaTCR9++CF69uyJoKAg+Pr62jxchVwqwWvjbBc8H9Q6ALJyQUzBMWVERDYUhrK7780qn1vsWTWZUo1kSTAAIMJ8DXqD0S61EbmyOoWyN998Ex9//DEmTZqE/Px8zJ07F/fccw8kEgkWLFhg5xLrV/8Yfzw+OBohGhWGtgnE6zeFNNu7LzmBLBGR3FBo3ZaoPG+x561dVnfAHlM7rDQNQ0puYfUHEDVydZrRf+XKlfjmm28wduxYLFiwAA888ACio6PRsWNH7N+/H0899ZS966xXL45ugxdHt6n0NQ70JyKyJTGXLUknVajrfJ5dUXOx9Eg2AOCFQhkiA2+7NCKXVqeesrS0NHTo0AEA4OHhYZ1Idty4cdi4caP9qnMC5ceUMZQREQGCqSyUKRR1H1/cwldh3U7I0d1WTUSNQZ1CWfPmzZGaalmmKDo6Glu2bAEAHDp0CEql8laHupwbyywBHFNGRAQAMnNZgFLcxu/8KJ+yYxNy9bfYk6hpqFMou/vuu7Ft2zYAwJw5c/Daa68hJiYGDz30EGbMqP06aM6Md18SEdmSlrt8qVTW/fKln5sUHgoJABF5OVwDk6hOY8r+9a9/WbcnTZqEiIgI7N27FzExMRg/frzdinMGnNGfiMiWVLSEMpMoQCmX1/k8giDgA9VS9Bb2QmMuQlzxX3B3c6/+QKJGqk49Zbt27YLRWHb7cu/evTF37lzccccd2LVrl92KcwZc+5KIyJZZBMyiAB0UkNdxmaUbfGQG+AhaSAQReelJdqqQyDXV6f+mIUOGICcnp0J7fn4+hgxpXEtllA9lpQZOiUFENEP6L7TQ/YB++Pa2z2V287duF+Zl3mJPosavTqFMFEUIglChPTs7G+7ujavrOdCz7M6ixJxiB1ZCROQcSo1mAALkMkW1+1ZH7l424XhpQdZtn4/IldVqTNk999wDwDIOYPr06TZ3WppMJpw4cQJ9+/a1b4UO1rbckktnUgocWAkRkXMovT6+ViW/vUuXAKDUNLNuSwuv3fb5iFxZrUKZRqMBYOkp8/T0hFpddteNQqFA7969MWvWLPtW6GC+7goEa1RIzS/FmdSCKnsJiYiaAlEUUWIQAQBq2e3/LlQFtwVOWLaDi8/f9vmIXFmtQtmyZcsAAAEBAViwYAHc3NwAAFeuXMG6devQtm1b+Pv73+oULqldiBdS80tRWGpEUk4Jwv3cHF0SEZFD6PU6LJR9jWKoUGSOBjDtts4n+kRBCzU8UIJWxgvINpshSG6/B47IFdXpJ//YsWNYsWIFACAvLw+9e/fGRx99hAkTJmDx4sV2LdAZxIZorNunU/JvsScRUeNmKCnEZNkOzJBtQl/jwds/oSBBvLw1ACBAyENuWsLtn5PIRdU5lA0YMAAA8Msvv6BZs2a4evUqVqxYgU8//dSuBTqD2HLjyk5zXBkRNWGG0qKybWndl1gqL9uni3Vbm2CHoEfkouoUyoqLi+Hp6QkA2LJlC+655x5IJBL07t0bV69etWuBzqBdSLnB/qkMZUTUdOmLy34HGqT2udteFt7Tui3J4rgyarrqFMpatmyJdevWISkpCZs3b8bIkSMBABkZGfDy8qrmaNfT3EcNjdoyazUvXxJRU2bQZpdtK33sck7/qI542fQo+usW4YWShyCKol3OS+Rq6hTKXn/9dTz//POIjIxEr1690KdPHwCWXrMuXbpUc7TrEQTBegkzvUCHLK2umiOIiBonc1FZKDOrfG+xZ83JZHKcD7wD18RAZJeYkVJgsMt5iVxNnULZxIkTkZiYiMOHD2PTpk3W9mHDhuGTTz6xW3HOJDaE48qIiKTFaWVP3O13t32noLIplo6nldjtvESupM73HQcFBaFLly6QlLt1uWfPnmjTpo1dCnM27WxCGS9hElHT5JZbNuZLGRBtt/N2LBfKTjCUURNVq3nKmrJ2NtNisKeMiJqevNxsdNEfAQQgD54ICLVfKIvxV2Ko/CR6iSfQPiUZovkrzldGTQ5DWQ1FB7hDIZNAbzTjLEMZETVB2iuHoRSMAICT3sMQKJXb7dwyiYA5qk3oYjgGANiXngDfYPuFPiJXwD9DakgmlaBNkGUakITsIhTpjA6uiIioYZ3WB+FDw31YZ+qLgqA+dj9/tk9n63bRZc5XRk0PQ1kt3BhXJorAWc5XRkRNzKHiQHxuuhvPGJ6EIqqf3c8vC+9h3XbPPGr38xM5O4ayWojluDIiaqJEUcSFLMt0QAqpgFAv+126vCEgqhOKRCUAILrkBESz2e7vQeTMGMpqwXa5Jd6BSURNR2qhEelay7CNtgEqSCWC3d9DJpPjgiIWABCIXOSmcx1MaloYymqhbbAnhOu/h7jcEhE1Jann9qOP5DRkMKJrqLr6A+qo/LgyLceVURPDUFYLbgoZWvhb1nq7kKaFwcSudSJq/PQlWgy79B5+VLyLI8rHMMBPW2/vVX5cmQfHlVETw1BWS+1DLePK9CYzx5URUZOQuf1zNINleaVLitYIDQmrt/fiuDJqyhjKaql7ZNlab4cSchxYCRFR/UtNOI2BeesAADpRDl3fF2Edx1EPbh5XlpPGcWXUdDCU1VLPcqHsAEMZETViJqMBQQffg0yw9FbtC5wM35AW9f6+VwOH4RvjGMzUP4cjee71/n5EzoIz+tdSTKAHvN3kyCs24NCVHJjNIiT1cBcSEZGjpf/1BYaYLwIArgohCBz0SIO8r9D2Trx7qSsAoCRVxMjYBnlbIodzqZ6yf/3rXxAEAc8884zDapBIBPS43luWX2LAhYxCh9VCRFRfUs/tx8CsHwEARlGCK91egVSubJD3jvFTwkctBQDsTypCutbQIO9L5GguE8oOHTqEJUuWoGPHjo4uBb2iOK6MiBqvosIctDn2JqSCCADYGfgggmJ6VHOU/UglAsa1sdxUZRaBP87zpipqGlwilGm1WkyZMgXffPMNfHx8HF0OepYLZfsZyoioERFFEWnbFiPo+t2WJ2TtETz40Qav445WXvAVCjFTuhH3XXgWJoO+wWsgamguEcpmz56NsWPHYvjw4Y4uBYBlZn8PpWU43t74LJjMooMrIiKyj/87l49nc+/G98bhyBU9oRvyLqSyhh9+7Ocmw+deP+A1+Up0xTmkn9jS4DUQNTSnD2WrV6/G0aNHsXDhwhrtr9PpUFBQYPOwN5lUgv4t/QEAucUGxCXl2f09iIga2sWsUnx9MAs6KPCacQY29/wOXv7BDqunOOYu67ZfwjqH1UHUUJw6lCUlJeHpp5/GypUroVKpanTMwoULodForI+wsPqZ5HBom0Dr9o7zGfXyHkREDaVIb8Z7O9JhuD5X6z3tNOjSsrlDawqJ7YdEBAEAOhpOoIBrYVIj59Sh7MiRI8jIyEDXrl0hk8kgk8mwc+dOfPrpp5DJZDCZTBWOmT9/PvLz862PpKSkeqltcOsA6/Z2hjIicmGi2YzEzYug1l4BALTyV2JGN3/HFgVAKpXibOBY6/PSuDUOrIao/jl1KBs2bBhOnjyJuLg466N79+6YMmUK4uLiIJVKKxyjVCrh5eVl86gPgV4qtA+1nPtUcgGS80rq5X2IiOrb1X0/Y4J2NTYoXsU9iv14eXAQ5FLnmH9R0+Ue6ETLmLb2OZthNugcXBFR/XHqUObp6Yn27dvbPNzd3eHn54f27ds7ujyMig2ybq87luzASoiI6ib58mkMuvoZAMBd0GFcaw2CPeUOrqqMj28ADqv6WLZRiMyTHPBPjZdThzJnd3fXUOv2L0euQRR5FyYRuY78wkJEHXgFKsEyOevfXmMQ2n28g6uqqDDmHuu23+W1DqyEqH65XCjbsWMHFi1a5OgyAADNfdzQN9oPAJCQVYQjV3MdXBERUc0YzSLyN7+LKKQAAC5LIuE78iUHV1W58Ng+uIIQAEA7wynkpsQ7uCKi+uFyoczZ3Nu17O6kH/ZfdWAlREQ1d3LnLxhh+AsAUAwlCoYshFRes7vcG5pUKsG5ZpYB/1pRhbjTp2t0HK9ekKthKLtNYzsGw8fNMv7itxOpSC8odXBFRES3djXhAsanfmp9fiZ2HtSB0Q6sqHpBPSbiU/N96Kf7FB+mtENmkfGW+x9JLsakHxPw4e70BqqQ6PYxlN0mlVyKKb0iAFguB3y/j71lROS8SkpLELn/VbgLlrsYj3mPgGenOx1cVfXUnj5IbjMD+fCA0Qz8fUV7y/1f3ZqCfJ0ZW+MLcSWXd2ySa2Aos4MH+0RYbx9feeAqSg0V508jInIGu3b/hZai5Y/HREkoVMNecXBFNdc/wt26fSr91tMQlV/9Lk176141ImfBUGYHzbxUGNfRMgg1t9iAtZweg4ic0K4ELT6+1hpTDC8jUQxEWr+FkCjcHF1WjUX5KOGhsHxsJSVeRtq1yzU6Lr+EfyiTa2Aos5MZ/aKs20v/TuAAUyJyKhlaA/6z17L6yD5zO/xf9++had7GwVXVjlwq4MF2Crwh+w5bFc/Dff9HNTouu5g9ZeQaGMrspENzDXpE+gAALmZose0sl14iIucgiiI+3pMBrd6ysOWgKA8Mi/F2bFF1NLatL8bIDkMhmNDdcAjXzuyv9piMam4KIHIWDGV29MjAsruXPv3rInvLiMgppJzZi7bpv0GAGYHuMjzVJwCC4BzLKNWWXOmGU5EPW58Hn/wcZrMlbOpNItIKDTCLItzkZR9vGVpDg9dJVBcMZXY0vG0g2gZb1sM8cS0fOy9kOrgiIiJATNiO9+Tf4mfFW3iigwAPZcV1g11JSM97kCCEAQBizRdw9cgmJObpcc/Ky5j2y1VMW3MVxQazdf9CvbmqUxE5FYYyOxIEAU8NbWl9/sHm8zCa+MuAiBwrpOgsAKCH5AK6h3s5uJrbJ5HKkNjucevzmAtL8NjaBBhMlqsTN1+u1Bl51YJcA0OZnY1qF4Q2QZ4AgNMpBfh6d83uDiIiqg8moxHRpgQAQKIQDLm7j4Mrso+g9sNwRhYLAIiWpOIB6V9V7ltq5B/H5BoYyuxMIhHwr3s7QnJ9uMaiPy8iPqPQsUURUZOVlZ4IhWDpOUpXRlWztwsRBOR3e8r69B35MuxXPomVI4F2gbbLRbGnjFwFQ1k96BzmjX8OaAEA0BvNeOGXEzCZ+UuBiBpecUZZb32Je5gDK7E/3xZdEefW1/rcHNwFQT5eWDAsGCO9ErFe8Sqel/0EHXvKyEUwlNWTuSNaIcrfMvv00cQ8fLf3imMLIqImScgtW/rNpIlwYCX1w33I8yhW+KHAtyMM7R+A0S0QXiopvpB+jE6Sy3hSth4hpmTeDU8ugaGsnqjkUvzrng7W5//efA5JOcUOrIiImiJ1UVkok/s1osuX15m9muPanf9DVr83URrQ0dpuUvtatwOFXOtNAETOjKGsHvVq4YeH+lj+Mi01mPHljksOroiImhrv0mvWbY9mLRxYSf0R5W4weITatBU162nd9kIxSjiujFwAQ1k9e25Ea+ti5T8eTMSqA4kOroiImgqTWYTGYJkvMUf0gtqzcdx5eSveF35C+Obp8Du7wtoWKOTiXGapA6siqhmGsnqmcZNj1oCyv05fXnsSn3O2fyKqZ6Io4r3tqXhB/wgAIF0eWs0RjYNb2iGoc87YtAUKediVoHVQRUQ1x1DWAOaNao1HB5YFsw+3XMDm0+kOrIiIGrt9iUX4O7EYB8S2+Mk4GLnRExxdUoMoDB9eoS0QedibWAQ9x5WRk2MoawCCIOClO9rg4X6R1rbt57hgORHVj1KDGV/sz7I+T+n+Ivy63u3AihpOUWh/mCUKmzYPoRjFBjOOpvBmK3JuDGUNRBAEPDeytfX5zguZyNbqsOF4CraeSedyTERkN+tOpCOr2DJhbLdQN9zRxreaIxoPs9wDRSF9bNqUsCxIzkuY5Oxkji6gKfFQytC6mSfOpxciraAU3d750/paqLcaYzoEYXT7YHSLaPyDcYmofuRlpmDO+QehkY3Af03j8ESvcAiC4OiyGlRh2DB4Xttpfa4WLKHsbAYH+5NzY09ZA3tlbFvr3ZjlJeeV4JvdCZi0ZB/OphY4oDIicnWG0mKodrwGX6EQz8h+xUeBm9Bco6j+wEZGGzrA5nmRIgAAkFpoQImBVyXIeTGUNbCBrQKwbe5g3NMlFDf+eJVKykKa0Sxi4R/nHFQdEbmqlNO74b12EjoZTwIAsqFBswEzHFyVY4hyd6T2fBWlvm2QPOAD/BrynKUdwNU8vWOLI7oFXr50gHA/N3w8qTPeGN8OaoUUJXoT/rXpHH48aJnDbNeFTOy+mIkBMQEOrpSInJ02LwulOz/EkOK/rG2lohynur6NIM+mM5bsZgXRd6Ig+k4AQIuCPGv7X5cK0SZAVcVRRI7FnjIH0rjJoZBJoHGTY+E9HbBoUmfra+/9fo6LmBNRlUSzGVf2r0XU75PQr1wgOyVti8MDliKoTW8HVudcuoe64caokfVn87E1nkNEyDkxlDmROzuFoH2oFwDgbGoB1h5LdnBFROSMcot1yF33HEYlLISvUAgAKBDdsDX8WUgnLkVAWOtqztC0hGoUeKJ32ZWHz/ZmIuf63alEzoShzIlIJAJeHtPW+vyjLedRajA5sCIicjZ/X9Hi0XXJuFjkZm07qOqHi3esRni/ByCRSB1YnXNpduhfCN35LEL+fgnj2mgwMsYTAKAzidhwNt/B1RFVxFDmZPpG+2Nom0AAQGp+Kb79O8HBFRGRM9DqTPj3rnS8vT0N+TozFhofwHHEYHur16G5+xN4+AQ5ukSno8y9AI+UPfBM+gsSfSGmdfGzXsbcdqmQy92R02Eoc0Lz72iDGzdkLt5xCRmFzjG3zrXcYmw6lQadkb13RA3p2pn9+N/a1dh2qdDa1jXcD6V3LUdIt3EOrMy5lfq1s26rsk/D312GzsFqAEBGkREJubwTk5wLQ5kTimnmiUk9wgAAWp0Rb/7fmWqOqH8Gkxl3f7kXj/1wBF9uv+TocoiaBJNBj+S/FmNw3NN4xfwVwoV0uMkleH5AIF4bGgRvN7mjS3RqJeVCmTr7FACgV5i7tW1/YlGD10R0KwxlTmreqDbwdbdM+rjxRCr+POPYBcwvZWqRWagDAPxn20WH1kLU2JnMIs4f2Qb3X+7H0PRlkAoi3AUd5nltw5IJYRjR0qvJzdJfF6X+HazbqixLKOsdXi6UJTGUkXNhKHNSvu4KvDaubND/a+tPobDU4LB6lDIOHiaqbyaziNPHDyD/58dw54X5iEAKAMAsCtjuOwkt73oZgR7sHaspg3sIjEpvANd7ykQRzTzkiPKx/MF7PkvHuzDJqTCUObEJnUMxsNX15UHyS/Hh5vMOq4VzphHVH5NZxNnj+6D9+RHcc2YO+ojHrK+dlbbB9h5LEDLqOUiknO+7VgQBpX7tAQBSfQHk2iQAQO8w9paRc2Ioc2KCIODdCe2hllt6qVbsv4rDV3IcUovBxPXiiOrDvsQifPPrRkw48zR6iset7RmCPw60fQWS+75D85jOjivQxdmMK7t+CbNvRFko23FZ2+A1EVWFoczJhfm64bmRrQAAogg89eMxZGt1DV7HzaHMzJ4zottSqDNhwbZULNiWivWFrXDVbJkKJ01ohoOt5yH3vvXw7nwXwLFjt8VmXFn2aQBAjJ8SoV6Wy8An0kqQVcRLmOQcGMpcwPS+kegR6QMASMkvxX1f7UN8RmE1R9mXwWQbwvTsOSOqswytAc/8dg37rt/9Z4YEq9wfxP7WLyLvvl+h6XofIOXYMXso9Y2FCEuwvRHKBEHAkBYeACyLlO9IaNjfp0RVYShzATKpBJ//oyuaeSkBAJezijDm07+xZOelBrusePP76IwMZUR1tW/nRuQWWNZf1CgleHlwM9x91/3w6XovBIYxuzIrPJDZ5Rmk9l6AjG7PWduHtPC0bn9zKJu9ZeQUOGrURTTzUuGXx/pixvJDuJihhd5oxsI/zmHtsWS8d08HdA33qXDM7ydTseV0GgwmEV3CvdG7hR/0JjO6hHlXezu93mhGYk4RgjRqeChlFUKZnqGMqE4yky/j+fyFeEypxjLJBHQeOxtBngxi9Sm3zQMV2pprFIjyUVgnkP3lVC4e6xVQYT+ihsRQ5kLCfN3wf3P64/1N57B87xWIInAurRD3Lt6Lf/QMxwuj2kBzfTLJ306k4MlVZXdwbTyZat1+ZGALmzU2b7bvUjaeXHUU2UV6uCmkWHhPB3ipbD80ePmSqG6Mx76HRBDhhWJ0CFQxkDnQ/R188P4uyxyQexOL8GhPf87/Rg7Fy5cuRiWX4o3x7bD2iX6IDfYCYLkBYOWBRAz7eAc2nbKEr6W3WDPz612XsfFEaqWv6Y1mPL/mOLKLLH89FutNeOanOPx4MNFmPwN7yohqxWQy4+r2bzGg4HcAQJGohF/PSQ6uqmkbGu1pXXYpXWvEpRwuu0SO5dShbOHChejRowc8PT0RGBiICRMm4Px5x83V5Uw6h3ljw5P98OrYtnBTWKbMyNLq8fjKo3jhl+M4mpgHAAjyUuG3Of3x6KAWNse/8MtxXMqseCv47ydTkZxXYtMmisCWm1YUYE8ZUc0ZjUbkrH8BI9OWQCpYbpo54H8PVB4Vhx1Q/ZDotfBK2AhN/Fqb9gGRHtbtPVc5PQY5llOHsp07d2L27NnYv38/tm7dCoPBgJEjR6KoiJP9AZYbAP45oAX+nDsI3SMsv9xFEfj58DXrPvNGtUb7UA3m39EWCQvHYELnEABAkd6Ex384gmK97eDWuKQ86/Z3M3pibIfgSt+bY8qIasZsNiNz41vor9tlbTuo7IOgoXMcWFXTIph0aLFhPIL3vwn/k18DZpP1tT7h7rhxwfLvq/xsIcdy6lC2adMmTJ8+He3atUOnTp2wfPlyJCYm4siRI44uzamEeKuxfEZPBHmpbNpn9o/CPV1Drc8FQcB793RAq2aWvwwvpGvx8q8nIYqWv9wzC3XYdTHTun/7EC+8MT4WUknFMRa8+5KoeiajEdm/vY7BxZsAAEZRgl2hj8BrwieQyjikt6GIUiWKA7sBAGSl2VBnxllf83OTITbQ8rszMU+PpHxewiTHcepQdrP8/HwAgK+vb5X76HQ6FBQU2DyaAg+lDD892huxwV7wc1fgtXGxeHVs2wqDVt0UMiye2g3u1y95rotLwX1f7cPkr/eh98JtuJxp+UuxZaAH/DyUCPRSYWCMf4X3Y08Z0a2VGsw4+cdi9C/aAgAwiQL+bjUfzQb+E4LEpX71NgqF4SOs216JW21eKz/D/64EXsIkx3GZ3wxmsxnPPPMM+vXrh/bt21e538KFC6HRaKyPsLCwBqzSsSL83LHxqf44/OpwzOwfVeVdRNEBHvjgvk7WicIPX83F/ss51vUtPZQyvDuh7Hs8sVvF7yHHlBFVzmQWsTW+AI+tT8SHWX2QKXrBIEqxPeYVBHe/y9HlNVna0AEwSy1zPXokbQfMZUM3BkZ6WC9hbo0v5Fq/5DAuE8pmz56NU6dOYfXq1bfcb/78+cjPz7c+kpKSGqhC5yAIQo1u6R7TIRhfP9gdzX3U1rbmPmrMHhKNzc8ORK8Wftb2YW0DKxzPnjIiW9q8DJz5ey1m/pqID3dnILXQiEx440XzbOxp9ybCetzp6BKbNFHuhqKQfgAAmS4Xbullw2ACPeToHGL5XZhaaMCaU7kOqZHIJQY1PPnkk/jtt9+wa9cuNG/e/Jb7KpVKKJXKBqrMtY2IbYZhbQJxLbcECpkEzbyUlQY6lVyK18fF4q3fzljbGMqILGPGUs7shlf8OnQpPYhuggkf6RYBsPwh0ylIjal9xyJYo3BonWRRED4Cnkl/AQA8E7eiOLiX9bWpnX1xPDUZZhH48XguRsd4wVvtEh+R1Ig49U+cKIqYM2cO1q5dix07diAqKsrRJTU6EomAcD+3aveb0d/yvb8RzHKKGn5RdCJnkZNxDUXHfkH7nM2IRbal8frfM/dJd2JH4EOY0skXHYLUVZ+EGlxRSD+YZWpIjCXwvLYD6d1ftK4x2r6ZGuPaaLDhbD5KjSI2nMvHQ138qjkjkX05dSibPXs2Vq1ahfXr18PT0xNpaWkAAI1GA7Wav+waWpugsrXiEnOKHVgJkWNkpSVCd3gF+hX8DqVgO51MFjQ46TMSPTreg9EhoVWcgRxJlKmgDR0Ar6tbINUXwD39oPWSJgBM6uCD387lwywCG88V4P72PlDJXWaUDzUCTh3KFi9eDAAYPHiwTfuyZcswffr0hi+oiSvfo8ZQRk1Foc6EnQlahJ5ejMn6Xy2N13vFTKKAOGV35Le4E806DEWIjEsmObvC8BHwumq5I9Yz8S+bUObvLkPP5u7Yn1SEvFITtsYXYnxbjaNKpSbIqUPZjfmzyDkEa9SQSQQYzSISc0qqP4DIRRkMBhxOLsa2hBLsTyqCwSRirKQ5Jl8fGlYiKnDIbwLcuk+Bp18wPG59OnIiRcG9UdSsB4qb9YA2tH+F1+9o5YX9SZapgS7lcJgGNSynDmXkXKQSAc191LiSXYzE7CKIosjFe6lR0Rbmo2D3YnTM+xN/G2Zht7mb9bU/zV1xVmiJ9MB+8Oo6Ec28mzmwUqorUarEtaFfVPl6bLOySbiTCziRLDUshjKqlWCNJZQV6U0o1pvgruSPELmeosIc5CfHQ5edAFn+VXiXJiHQkIzW5utLlAnABOnf+NPcDRqVFENaeGBESy8IvisRzD9EGjVPhQS+ailySkw4lV6K3BIjfHgXJjUQ/qRRrfi4l42ZySsxMJSRS9Hr9Sjd9Aq6aXdBIlQ9PEInyuHj4YZ3ugWhS6g7ZJUsNUaNkyAIGN7SEz+fzINZBI6mlGBYtGf1BxLZAW8roVrxdiubbym3iF375BpEUURcajGe/iMNe/L9Kg1kOlGOeIRhn9swnBrxI/wnvI8eYR4MZI2YIu8S/OM+h1vaIZv27qFlNzUdS+FNTdRw2M1BteKtLuspyy8xOLASopopNpjx3vY0HEq2fLguwr1oJhQgTCODTtMCCr9IKAOioPYJhUIuQ9Ur61Jj4pZ2CGHbZwMAFAVXUBzUw/pa20A1lFIBOpOIYynFHD9LDYahjGrF260slGWzp4ycXGF2Crbv2IxDBQOsbc01CvgOeQMBPlz5oykrDuwCk9wDUoMWnsm74Hl1KwojLIuWK6QC2gepcSS5GFnFJiTlGxDuzVUZqP7x8iXVSrivu3X7YnqhAyshurW0hFMI2zID8/SLcadkD9zkEszs7ofPx4chkoGMJDIURoy0PvU5v8rm5S7BZROU8xKmLb3RjA3HU5CYze+LvbGnjGqlXYiXdft0SoEDKyGqnMloQPLf32NAyrdQCZZL7M8q1uHOUfchzLf6JcWo6chu9zC84y0TAquzT0NanAmTWwAAoGuIG3B9Ca3fzheglb8KbQIqXx+4KfnrXDrmrTmB7CI9VHIJFk3qgtHtgxxdVqPBnjKqleY+anipLFn+dEq+g6shKmPQleDK3p/hvmYiRqR+ZQ1k56QxyLljCQMZVWB0a4as9v+0Pve89pd1O8pXAY1KCgBIzNPjmY3XcCCp6fYMnUkpwBMrj2DG8sPWoSulBjMeX3kE/919mZO92wlDGdWKIAiIvd5bll6gQ5aWM16TA4kispLOI+X399H817EYdfVDRCDV+vJOrzthmvAtPLwDHVgkObPCsGHWbc/EslAmEQSbS5gA8Pb2VDQ1yXklePanOIz5dDd+P5lW4XVRBN7ZeBavrT8Fo8nsgAobF16+pFprF6LB/ss5ACx/PQ1sFeDgiqgpMZpFxGfrcPBaEQ5cLcIPxU8hSMi12eeULBbZnWYjqFWPKs5CZKHXtIDOMwLKwqtQZ8ZBWpIFk9ofANA73B07ErTWfY1NKHMUlBrw5fZLWLonAfpyX7ifuwIv3dEG93Ztjv9su4j/bLsIAPhhfyKSckrw+T+6wFPFNWDriqGMau3mcWUMZVSftPnZyE04CknacYjaDDxR/Ch0prJLJdtlnfGAbDv0ogyH3PrD1O4+NGvZHUFNfOwP1ZAgoDB8GJSnl0KACM+kv5DX6n4AQM/m7hV2N4siJI34ZyslrwSrDyXhh/1XkVPuDntvNznmDI3BlF7hUMktl3WfHdEKEX5uePF/J2Awidh5IROTluzHsod7oJmXqqq3oFtgKKNaaxeisW5zXBnZk9lkRPa18yhJioN79kk0Lz6L1rC9ZKQy/QO6ckuAH/AYjkDfVvDtfBf8Nf4NXTI1AtqwYfCO/xWF4cNREtDZ2u5ebsmlG1ILDQj1alzTY5jMInZdyMTKA1fx17kMmMsND1NIJXi4XySeGNISGnXFHrB7ujZHiLcaj35/BPklBpxJLcA9X+7FdzN6oGUgV0KoLYYyqrXoAHcoZBLojWac4R2YVANpZ/dAfXYNmusuwggZrilb4rfIl+GhkMBdIYVKKiLyymr0zliNtii65bkGuiVCG9QTHYPU6NncHf7uLRvoq6DGSufdEpcm/A5IKn4kLhwVgkfXJVmfJ+ToG00oS8svxbq4ZHy/7yqS80psXpNKBIzrGIznR7au9iaZ3i388L/H+2L6soO4lluC5LwS3PPlXiyd3gPdIzkdc20wlFGtyaQStA3yxPFr+UjILkJesd5m+SWiG0wGPXL/eBODirbatGeXuGPVcdtxYBsVf0IjsQ1kelGGeGk0MjxjYWrWCV6RXfCkX3C9101NjCAAQuUfh5E+Srw2JAhvb7cMcn97exr+mB7tspcws7U6bDuXgfVxydh7KRs33zQZrFFhco9wTOoRhiBNzS9Btgz0wK9P9MXDyw7hdEoBCkqNmL7sEH74Zy90DvO27xfRiDGUUZ30jPLF8Wv5EEVgT3w2xnbkByXZMhkNKNj4CvqV7LS2FYpqSGDGBXPzCvu7oRQmUcBxRRdk+3WHonkn+EV0hFyhBH+6qKEJxlKIMksoifK1/aPzXGYpYgPVlR3mNFLzS7AnPht747Nw+Gou0gtKobvFnQqDWwdgaq8IDG4dAJm0bhMzBHqq8NOjffDY90fwd3wWtDojpi09iNWP9EbbYK/qT0AMZVQ3A2IC8M3uBADA7ouZDGVkoyA7FYq/XkNv4wkAlh6vHWFPonmf+yBIpAjILcQ7JTIU6c0oNphRpDfjTME0ZLTqDe+A5qg4vJqoYUj0BfA9+wO843/FldHfw+gejGBP27FUaYVGxDrZLCsms4ijibnYeiYd286m41LmrYcBAECknxvu6hyKu7uEItLfPv/XeShl+Oah7pi+7CAOJOQgv8SAqf89gJ8e7YOWgR7Vn6CJYyijOukZ5WsdV7b7YhYX7CUAlskkLx38DUOvLoK3YJlKQC/KsLPN64joOtq6X6ifBqEVjp7YYHUSVcXn3I/wO7McAOB1dQtyYqdBIgjoF+GOPVctQSet0ODACsuIoiWIrTl8DVvPpFe5HrFaLkW4rxvUCilUcgk6NvfGHe2D0DnMu15+b6sVUnw7vQce/PYAjiXmIbtIjyn/3Y81j/ZFuB8ncb4VhjKqE5Vcip6Rvvg7PgvJeSVIyCpCiwD+FdRUlRjM+L9z+fjlVB466E24R2EJZBnwwZmubyK8TW8HV0hUM9rmA+F/+lsAgLywbID/3bHe1lBWbHDshGU5RXr8evQafjqUhIsZ2gqvSwSgS7gP+rX0R/+W/ugc5g2FrGHnivdQyrD84Z544Ov9OJNagPQCHR74Zj9+fqwPQr2d+9KvIzGUUZ0NiPHH3/FZAIAV+65iUKsA9I/xh7yO4xHIcUxGIzLO78WR4mbQKgMwoqUX/Nyq//VQUlSIXWeu4r8X3VCgs3xQ/Y32OGKOgcE9GIqhLyPYi3dfkeswuIdYtxXlQplSVtajVGps+CWF9EYzjibm4pcj17DheIrNhK6ApTdsYCt/jIwNwtA2gfBxd/zNVxq1HN/P7InJX+/HxQwtkvNKMPW/B7D+yX7w4gSzlWIoozob0yEYC/84BwBYvvcKlu+9gjlDW+K5ka0dXBnVRsqZPWhx4kMMFpOxVv8E1pv7Y92ZfLw9PBgx/pXffVVcmIec/d+je+avMIih+Fj/BgABAoBBUZ7I7vAVIvw4RxG5HrPCCwb3YMiLUqHOPgmJXguzwgOqcj1Npbc5tb8oiigoNSK/2ICCUgMMJjPaBnvBaBaRWWhZvi4xuxhXsotwOasIV7KKcDmzCCUGU4Vz9Yj0weQe4RjTIRhqhfS26qoPfh5KrJzVC5OX7MflrCIkZBXhjfWn8cmkzo4uzSkxlFGdhfm6oXUzT5xPL7S2ffZXPEOZixDNZiTv+hZDUv4LiWD5yz9MyAQA5JaY8NzvyXgs6Dzae2gh94+CVOkGY2E2hKR96Jr7B7oIxYAA9BAuoK/kDNRRvTG5ow/CvR3/FzpRnQkCtKED4HPhZwhmI9xT96EwYoRNT5muFj1lBpMZJ67lY//lbJxNLcDV62GrsNRY5xK9VDLc1z0MD/QMd4nB84GeKnw3oyfG/Gc3CnVGrD2WjMGtA3BX54ojS5s6hjK6LQNi/G1CGbmO5O1fYVjGcuD6Z805SQyC2/ZFy1Ql4rN10JlEdMv4HwZmnwSu3nTw9WOMogQH3IfgycEd4RfcrCHLJ6o32tCB8LnwMwDAI3kXCiNG2PSU6W6x8LbBZMbZ1ALsu5SNfZezcSghB0X6ij1ctSERgOY+bugW4YMBMf4Y3T4IbgrX+vgO83XDO3e3x9Or4wAA8389id9PpiI6wAMtAjwQHeCOFgEela4a0JS41r8qOZ3+Mf74798Jji6Dailp3y8YnrHc+nyb3xSEDH8S3SVSdDCasexINnbHZ6GXcK7S43WiDIc9hkLa85/wD4psmKKJGkhxYFeY5B6QGrRwT9kLmI1Qlr98aRCv/9eEC+mFOJVcgJPJ+Tidko9zaYUVxnuVJxGAUB81wn3d4OuuhKdKhowCHa7lFsPHTQF/TyUCPJQI8VYhyt8dkf7uCPNxa/CB+vXhrs6h2H4uA+viUlCsN2Hz6XQA6Tb7+Hso0SLAHdHXg5oltLmjuY8bpJLGf4c/Qxndll5Rfo4ugWrBrCtC+s6vMDz7J2vbn4EzEDbsMetzpUyCx3oFYEZXb+yN/xCmtFNQFKcCZiNEmRp675bwbTMI/j5ONlETkb1IZCgK7guvxC2QGgrhlnEMYrPu1pcvZpdixv+uIl17CUbzrS9lBngq0aeFH/pE+6FHpA/Cfd0bRcCqq7cntIdJBDafSoO+kh7HLK1lTN3BhBybdoVMgkg/N2tIu9HD1iLAvVHdNMBQRrdFrZDCSyVDQbnxEcV6o8t1rTcVxb+9iMH6g9bnf7uPsglk5SnkcoS07QO07dNQ5RE5DW3zgfBK3AIAcE/dh+KgHlDKBOiMIooNIooNFecqEwQgyt8d7UM06BHliz4t/BAd4M45HMvxVMnx2QNdYDKLSM4twaVM7fVHES5f/2+WVlfhOL3RjAvpWlxIrzgFSICn0nr580Zoiwn0QKi32uW+94Io3rzyVeNSUFAAjUaD/Px8eHlxmYf6cCalAGM+3W19/sjAFnh5TFsHVtT4nT9/vtbHpOeXwOu3GegsuQQASBCaQzvhByhUnMyR6GYSfSH8Tn0DbeggqHLOQqFNwgt592LLNUsvl0omICrAE7HBXmgf6oX2oRq0DfaCh5J/kN6u/BIDLmdqcTmzCJfK/fdqdnGlvWtV8XNXoFOYNzo190bHMA06NfeGbwNMFXI7uYOhjOwi8qWNNs+v/GusgyppGuoSyv64kI9P96RjjnQthkiPQT/6E2i4uDfRLSlzLyBy01QAQG7EGOxrMx/eKil81VK0adPGwdU1LUaTGddyS3A5S4tLGUU2/83SVr6awc3Cfd2uBzUNOod5o12Ixu5TidxO7mCkJ7sQBKBxx3vXF5+tgxkS/Md0L5qPfAod/TirNlF1VFknrds+V39HdN8FjiumiZNJJYi8fvPD0JvycH6xAZeytLiUYbkEeja1AMev5SGv2PYyc2JOMRJzivF/x1MAAFKJgFbNPNGpuQadwrzRsbkGrZp5OmwSdIYysgsGMud3KbtsnEa0L+cSI6oJUap0dAlUAxo3ObqG+6BruI+1TRRFJOYUIy4pDyeu5eN4Uh5OJudDV+7uWJNZxNnUApxNLcDqQ5YVHJQyCdqHWsKZRAB83BTwdVfAz8PyX193BfzclfB1V9j9pg2GMqImQG8w4nKuJZQFecjg7oQzfxM5I0Vhos1zwVgKUVb5ShfkXARBQISfOyL83K0T1RpMZlxIL8TxpHycuJaHuKQ8XEgvRPmbaHVGM45czcWRq7nVvoenUgZfjxtBzfJfN6HuC9YzlJFdBHmpkFZQan2uM5qglPGD31mkHV6LLZLvsE3WFVd8xwGIdHRJRC7BqA6wea4oTITOp5WDqqHbJZdK0C5Eg3YhGvyjVzgAy4wBp1MKcDwpD8evWcLa1eziGp2vUGdEoc5os79ZV7NjK8NQRnbxxZSuuHfxXuvzEj1DmbMQzWZEXP0FEZIMzJBsws7w8Y4uichl5MXcC58Lq62LkysKrjKUNTJuChl6RPqiR6SvtS23SI/LWVrIpRLkFhuQU6RDtlaPnCLLI7uo3LZWZzMt1O1gKCO76Bbhgyh/dyRkFQEASgwmeDu2JLou9fgWDBEtqy5ckEQjqGX3ao4gIitBgowuzyAg7jPovSJhUmocXRE1AB93Bbq5+1a/43UGkxm518NaUlo2Ri2q2/sylJHddAn3toay4ttc643sw2TQI+r8EuvzKzEPI8LFJlMkcrSi0AEoCh3g6DLIicmlEgR6qRDopUKoe93P03TXeiC7U8vLLleWMJQ5hbQdXyFSTAYAnJG2Rnjn4Q6uiIiIqsJQRnZTPpSVGhjKHC3j8nEMyPwRAGAUJcjq/jwECf+XJyJyVvwNTXZTflbkEoYyh9LmZ6PlgfmQC5Z/h13+k9CsRScHV0VERLfCMWVkNypevnQK+pIiqDY/i2BkAQDOSFohaMhsB1dF5Np8z3wHzaV1kBhLcW3Qx9D5cn1fsj/2lJHd2IwpY0+ZQ5QYzHhnRxqy9XIAQBY0KBr+AaRyzuBPdDsk+kIotMmQlWZDYqj7PFREt+ISoeyLL75AZGQkVCoVevXqhYMHDzq6JKqEh7Ks47XQTnO2UM2lFRow9/drOJAm4p+G57HWPBBnen0ILy46TnTbRFnZWrESY4kDK6HGzOlD2U8//YS5c+fijTfewNGjR9GpUyeMGjUKGRkZji6NbuLjXtYbk1ukd2AlTUtpcSGO79mIR9cl4nKO5fsuVyghjHiH48iI7MRcbmkliYmhjOqH048p+/jjjzFr1iw8/PDDAICvvvoKGzduxNKlS/HSSy/V/EQluYC8mktqCk9AWu5bYtQBNeqmFgC1t22TTguYa7D+lVQBKG6a1KQkD0ANVviWuwOycpelTEZAX1j9cQCg1ADl78QzlADG0qr3v0GQAiov27bSAkA0IUBWAi9oAQC5ORnIz7GdZFGUKm3+2oQoQtDlWzareDvx+krnosIDkJT924gmPYRK/m0qWxhdVHlDLP8OhiIIpuv/NlW8sQgAEjnEm/9tSgsA0VzJ/rYnEmVqQKYsq8dsglDu3+ZWX6+o8AQk5VZDMJYCBtt/m6zsbJSWlkKy9yP0M+xHB5jwtel9xKM5gj3leGt4MMK9ecmSyF7M5X53yYozIS3OBAq8AHd/QCov21FfDOgKqj+hIAE8Am3bSvIsnzvVkasA1U2T2GozUaPPDaWX5fgbTIbrnzk14OZn+7mhL7J8dlRHIgXUPrZtpfmAuQbDXGQqQOFW9lwULZ/nNaH0tP23MeoBQ1HNjr25Xn0RYKpBZ4NEXv0+t+DUoUyv1+PIkSOYP3++tU0ikWD48OHYt29fpcfodDrodGU/1Pn5lg/9gvfbA8pqJs2csQVoFlv2/PhPwO/PVV+oewDw1DHbtnVPAGc3VH9s+/uA8Z/Ytn3cAbgeVm7pzs+BdhPKnqccA76r4RI6T58A3MrNVvz3f4DdH1R/XFBH4OHfbdtW3A0kH0ILAH/f+BYfvf4oZ4lxLD433mN9rkIpDqtqNgD9Qd2LOCaWLW0yTHIE/1F8We1xxaICPXWLbdpel32H+2W7qj12q6krnjXY1rdZ+QJChexqj33T8CDWmAZbn0cKqfhN+Wq1xwHASN37SBH9rc+nSP/EfPmPNvt0KLd9I67NMP2CvyKfw0NdfOAu00OrZW8lkb0IpQLcdJbQo9r/CVT7P0EBYOfPjSdv43Ojq4t9btwDJB+q/tg+c4DBL5Y91xcBH7WuWb1TfgXCe5Y9P/c7sPaR6o+TqYF5F23b/ngBiFtV/bExo1Ew8mMAZZ0KtSI6seTkZBGAuHfvXpv2efPmiT179qz0mDfeeEOE5c8FPvjggw8++OCDD4c8Ll26VOvc49Q9ZXUxf/58zJ071/o8Ly8PERERSExMhEajcWBljVtBQQHCwsKQlJQELy+v6g+gOuP3umHw+9xw+L1uGPw+N4z8/HyEh4fD17fma2fe4NShzN/fH1KpFOnp6Tbt6enpCAoKqvQYpVIJpVJZoV2j0fCHsAF4eXnx+9xA+L1uGPw+Nxx+rxsGv88NQ1KHFVSc+u5LhUKBbt26Ydu2bdY2s9mMbdu2oU+fPg6sjIiIiMi+nLqnDADmzp2LadOmoXv37ujZsycWLVqEoqIi692YRERERI2B04eySZMmITMzE6+//jrS0tLQuXNnbNq0Cc2aNavR8UqlEm+88UallzTJfvh9bjj8XjcMfp8bDr/XDYPf54ZxO99nQRTrcs8mEREREdmTU48pIyIiImoqGMqIiIiInABDGREREZETYCgjIiIicgJNKpTdeeedCA8Ph0qlQnBwMB588EGkpKQ4uqxG58qVK5g5cyaioqKgVqsRHR2NN954A3o912K0t3fffRd9+/aFm5sbvL29HV1Oo/LFF18gMjISKpUKvXr1wsGDBx1dUqOza9cujB8/HiEhIRAEAevWrXN0SY3SwoUL0aNHD3h6eiIwMBATJkzA+fPnHV1Wo7N48WJ07NjROjlvnz598Mcff9TqHE0qlA0ZMgQ///wzzp8/j//973+4dOkSJk6c6OiyGp1z587BbDZjyZIlOH36ND755BN89dVXePnllx1dWqOj1+tx33334fHHH3d0KY3KTz/9hLlz5+KNN97A0aNH0alTJ4waNQoZGRmOLq1RKSoqQqdOnfDFF184upRGbefOnZg9ezb279+PrVu3wmAwYOTIkSgqKnJ0aY1K8+bN8a9//QtHjhzB4cOHMXToUNx11104ffp0jc/RpKfE2LBhAyZMmACdTge5XO7ochq1Dz74AIsXL8bly5cdXUqjtHz5cjzzzDPIy8tzdCmNQq9evdCjRw98/vnnACwriYSFhWHOnDl46aWXHFxd4yQIAtauXYsJEyY4upRGLzMzE4GBgdi5cycGDhzo6HIaNV9fX3zwwQeYOXNmjfZvUj1l5eXk5GDlypXo27cvA1kDyM/Pr9PirEQNTa/X48iRIxg+fLi1TSKRYPjw4di3b58DKyOyj/z8fADg7+R6ZDKZsHr1ahQVFdVqWcgmF8pefPFFuLu7w8/PD4mJiVi/fr2jS2r04uPj8dlnn+HRRx91dClE1crKyoLJZKqwakizZs2QlpbmoKqI7MNsNuOZZ55Bv3790L59e0eX0+icPHkSHh4eUCqVeOyxx7B27VrExsbW+HiXD2UvvfQSBEG45ePcuXPW/efNm4djx45hy5YtkEqleOihh9CEr+DWSm2/1wCQnJyM0aNH47777sOsWbMcVLlrqcv3mYioJmbPno1Tp05h9erVji6lUWrdujXi4uJw4MABPP7445g2bRrOnDlT4+NdfkxZZmYmsrOzb7lPixYtoFAoKrRfu3YNYWFh2Lt3b626F5uq2n6vU1JSMHjwYPTu3RvLly+HROLyfwM0iLr8THNMmf3o9Xq4ubnhl19+sRnfNG3aNOTl5bF3vZ5wTFn9e/LJJ7F+/Xrs2rULUVFRji6nSRg+fDiio6OxZMmSGu3v9AuSVycgIAABAQF1OtZsNgMAdDqdPUtqtGrzvU5OTsaQIUPQrVs3LFu2jIGsFm7nZ5pun0KhQLdu3bBt2zZrQDCbzdi2bRuefPJJxxZHVAeiKGLOnDlYu3YtduzYwUDWgMxmc60yhsuHspo6cOAADh06hP79+8PHxweXLl3Ca6+9hujoaPaS2VlycjIGDx6MiIgIfPjhh8jMzLS+FhQU5MDKGp/ExETk5OQgMTERJpMJcXFxAICWLVvCw8PDscW5sLlz52LatGno3r07evbsiUWLFqGoqAgPP/ywo0trVLRaLeLj463PExISEBcXB19fX4SHhzuwssZl9uzZWLVqFdavXw9PT0/r2EiNRgO1Wu3g6hqP+fPn44477kB4eDgKCwuxatUq7NixA5s3b675ScQm4sSJE+KQIUNEX19fUalUipGRkeJjjz0mXrt2zdGlNTrLli0TAVT6IPuaNm1apd/n7du3O7o0l/fZZ5+J4eHhokKhEHv27Cnu37/f0SU1Otu3b6/053fatGmOLq1Rqer38bJlyxxdWqMyY8YMMSIiQlQoFGJAQIA4bNgwccuWLbU6h8uPKSMiIiJqDDjQh4iIiMgJMJQREREROQGGMiIiIiInwFBGRERE5AQYyoiIiIicAEMZERERkRNgKCMiIiJyAgxlRERERE6AoYyIiIjICTCUERHdhuzsbAQGBuLKlSt2P/fkyZPx0Ucf2f28ROScGMqIqEmZPn06BEGo8Bg9enSdzvfuu+/irrvuQmRkZI32Hz9+fJXvtXv3bgiCgBMnTgAAXn31Vbz77rvIz8+vU21E5FoYyoioyRk9ejRSU1NtHj/++GOtz1NcXIxvv/0WM2fOrPExM2fOxNatW3Ht2rUKry1btgzdu3dHx44dAQDt27dHdHQ0fvjhh1rXRkSuh6GMiJzOf//7X3Ts2BFqtRoajQZDhw616/mVSiWCgoJsHj4+PrU+z++//w6lUonevXvbtJvNZixcuBBRUVFQq9Xo1KkTfvnlFwDAuHHjEBAQgOXLl9sco9VqsWbNmgoBb/z48Vi9enWtayMi18NQRkRO5ddff8ULL7yA1157DefPn8fevXvx3HPPObqsSu3evRvdunWr0L5w4UKsWLECX331FU6fPo1nn30WU6dOxc6dOyGTyfDQQw9h+fLlEEXResyaNWtgMpnwwAMP2JyrZ8+eOHjwIHQ6Xb1/PUTkWAxlRORUzp8/j4iICIwYMQLh4eFo164dxo4da9f3+O233+Dh4WHzeO+996yvy2QydO7c2fooKSmp9DxXr15FSEiITZtOp8N7772HpUuXYtSoUWjRogWmT5+OqVOnYsmSJQCAGTNm4NKlS9i5c6f1uGXLluHee++FRqOxOV9ISAj0ej3S0tLs9eUTkZOSOboAIqLyZs2ahZ9++gm+vr5wc3PDyZMnERUVVWG/l156Ce+///4tz3X27Fm0adOmQvuQIUOwePFimzZfX1/rtre3N+Li4qqttaSkBCqVyqYtPj4excXFGDFihE27Xq9Hly5dAABt2rRB3759sXTpUgwePBjx8fHYvXs33nrrrQrvoVarAVjGrxFR48ZQRkROw2AwYPLkyejbty++/fZbaDSaKu9qfO655zB9+vRbnq9FixaVtru7u6Nly5a3WS3g7++P3NxcmzatVgsA2LhxI0JDQ21eUyqV1u2ZM2dizpw5+OKLL7Bs2TJER0dj0KBBFd4jJycHABAQEHDb9RKRc2MoIyKnsXbtWsTHx+PPP/+sdt+AgIB6Cyp5eXno3LkzAKB79+7473//W+l+Xbp0qXBnZGxsLJRKJRITEysNWTfcf//9ePrpp7Fq1SqsWLECjz/+OARBqLDfqVOn0Lx5c/j7+9f9CyIil8BQRkROQ6/XIzU1Fd9//z0GDBgArVaLPXv2YObMmZDJ7PfrSqfTVRijJZPJrMGnppcvR40ahfnz5yM3N9d696anpyeef/55PPvsszCbzejfvz/y8/OxZ88eeHl5Ydq0aQAADw8PTJo0CfPnz0dBQUGVvX67d+/GyJEj6/7FEpHLYCgjIqcxefJkHDt2DC+//DLS09Ph6+uLYcOG4dFHH7Xr+2zatAnBwcE2ba1bt8a5c+dqdZ4OHTqga9eu+Pnnn21qfPvttxEQEICFCxfi8uXL8Pb2RteuXfHyyy/bHD9z5kx8++23GDNmTIUbBgCgtLQU69atw6ZNm2pVFxG5JkEsf082ERHB398fWVlZNdp348aNmDdvHk6dOgWJxL43tC9evBhr167Fli1b7HpeInJO7CkjIroNY8eOxcWLF5GcnIywsDC7nlsul+Ozzz6z6zmJyHmxp4yIiIjICXDyWCIiIiInwFBGRERE5AQYyoiIiIicAEMZERERkRNgKCMiIiJyAgxlRERERE6AoYyIiIjICTCUERERETkBhjIiIiIiJ8BQRkREROQEGMqIiIiInABDGREREZET+H/fdT8hGX2YtwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 700x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "fig = plt.figure(figsize=(7,4))\n",
    "\n",
    "plt.plot(vasp_dos_energies-vasp_fermi, vasp_dos[0,:], lw=2, label=r'VASP total DOS')\n",
    "\n",
    "plt.plot(plo_dos[:,0],plo_dos[:,1]+plo_dos[:,2]+plo_dos[:,3], '--', lw=2, label=r'PLO $t_{2g}$')\n",
    "\n",
    "plt.axvspan(-1., 1.5, facecolor='gray', alpha=0.3)\n",
    "\n",
    "plt.xlim([-3,3])\n",
    "plt.ylim(0,12)\n",
    "plt.xlabel(r'$\\epsilon - {\\text{E}_\\text{F}}$ (eV)')\n",
    "plt.ylabel(r'states/eV')\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the gray area highlights the energy window for the PLOs. The total DOS of VASP (blue) coincides with the PLO DOS in the window, as we re-orthonormalized the projector functions in the given window, picking up also Oxygen weight. This setting is closest to the result of maximally localized Wannier functions created with [wannier90](http://www.wannier.org/) without running the actual minimization of the spread. Note, for a proper comparison one can use the hydrogen projector in VASP by using the the line `LOCPROJ=  2 : d : Hy`, instead of `Pr`. \n",
    "\n",
    "\n",
    "## Converting to hdf5 file\n",
    "\n",
    "Finally we can run the VASP converter to create a h5 file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading input from vasp.ctrl...\n",
      "{\n",
      "    \"ngroups\": 1,\n",
      "    \"nk\": 3375,\n",
      "    \"nkibz\": 120,\n",
      "    \"ns\": 1,\n",
      "    \"kvec1\": [\n",
      "        0.2602750065719439,\n",
      "        0.0,\n",
      "        0.0\n",
      "    ],\n",
      "    \"kvec2\": [\n",
      "        0.0,\n",
      "        0.2602750065719439,\n",
      "        0.0\n",
      "    ],\n",
      "    \"kvec3\": [\n",
      "        0.0,\n",
      "        0.0,\n",
      "        0.2602750065719439\n",
      "    ],\n",
      "    \"nc_flag\": 0\n",
      "}\n",
      "\n",
      "  No. of inequivalent shells: 1\n"
     ]
    }
   ],
   "source": [
    "# import VASPconverter\n",
    "from triqs_dft_tools.converters.vasp import *\n",
    "\n",
    "\n",
    "# create Converter\n",
    "Converter = VaspConverter(filename = 'vasp')\n",
    "# run the converter\n",
    "Converter.convert_dft_input()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting h5 file `vasp.h5` can now be loaded as sumk object via:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SumK\n",
    "from triqs_dft_tools.sumk_dft_tools import SumkDFTTools\n",
    "\n",
    "SK = SumkDFTTools(hdf_file='vasp.h5', use_dft_blocks = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here one should carefully determine the block structure manually. This is important to find degenerate orbitals and spin-channels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: No Sigma set but parameter with_Sigma=True, calculating Gloc without Sigma.\n",
      "found 1 blocks of degenerate orbitals in shell 0\n",
      "block 0 consists of orbitals:\n",
      "  up_0\n",
      "  up_1\n",
      "  up_2\n",
      "  down_0\n",
      "  down_1\n",
      "  down_2\n"
     ]
    }
   ],
   "source": [
    "G = SK.extract_G_loc(transform_to_solver_blocks=False)\n",
    "SK.analyse_block_structure_from_gf(G, threshold = 1e-3)\n",
    "for i_sh in range(len(SK.deg_shells)):\n",
    "    num_block_deg_orbs = len(SK.deg_shells[i_sh])\n",
    "    mpi.report('found {0:d} blocks of degenerate orbitals in shell {1:d}'.format(num_block_deg_orbs, i_sh))\n",
    "    for iblock in range(num_block_deg_orbs):\n",
    "        mpi.report('block {0:d} consists of orbitals:'.format(iblock))\n",
    "        for keys in list(SK.deg_shells[i_sh][iblock].keys()):\n",
    "            mpi.report('  '+keys)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This minimal example extracts the block structure by calculating once the local Green's functions and then finds degenerate orbitals with a certain threshold in `Gloc`. Afterwards the result is reported, where 1 block is found with size 6 (3x2 orbitals for spin), where a all 6 orbitals are marked as degenerate. This is indeed correct in cubic SrVO$_3$, as all 3 t$_{2g}$ orbitals are degenerate. Note: for a magnetic calculation one has to break the symmetry between up and down at this point manually. Moreover, one can reduce the threshold for example to `1e-5` and finds that then the degeneracy of the 3 t$_{2g}$ orbitals is not found anymore, resulting in only two degenerate blocks for spin up and down, each with size 3x3.\n",
    "\n",
    "Afterwards the exact same DMFT script as in the [Wien2k tutorial](../srvo3.html) can be used. For a more elaborate example including charge self-consistency take a look at the [VASP NiO example](../nio_csc.html)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}