yambo/1_YamboCalculation_G0W0.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "02206d5f",
   "metadata": {},
   "source": [
    "# YamboCalculation: setting up the first G0W0 calculation\n",
    "\n",
    "Once completed the [preliminary DFT step](../prerequisites/0_2_QE_starting_point.ipynb), we are ready to run a G0W0 calculation, for simplicity in the common Godby-Needs plasmon-pole approximation (PPA). In particular, we will compute the $\\Gamma$-$\\Gamma$ band gap of Silicon. \n",
    "\n",
    "The *p2y*+*yambo init*+*yambo* run is submitted by means of the *YamboCalculation* *Calcjob* class, implemented in the *aiida-yambo* plugin. This does not involve any post processing decisional step (e.g. error handling, automatic determination of quasiparticle band gaps and/or excitonic eigenvalues).\n",
    "\n",
    "Input file generation and output parsing is done using the [*Yambopy*](https://github.com/yambo-code/yambopy/blob/master/README.md) package as python dependence.\n",
    "\n",
    "As in the [*02_QE_starting_point*](../prerequisites/0_2_QE_starting_point.ipynb) module, we create a builder instance and then we submit it to the AiiDA scheduler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "960cd7e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from aiida import orm, load_profile\n",
    "load_profile()\n",
    "\n",
    "from aiida.plugins import CalculationFactory\n",
    "from aiida.engine import submit\n",
    "\n",
    "import yaml"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "208aa0b3",
   "metadata": {},
   "source": [
    "### Loading the YamboCalculation class and the corresponding builder:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "65c9d7ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "YamboCalculation = CalculationFactory('yambo.yambo')\n",
    "builder = YamboCalculation.get_builder()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76398a89",
   "metadata": {},
   "source": [
    "It is possible to inspect all the possible inputs that we can provide to the builder, using the tab completion feature:\n",
    "\n",
    "```python\n",
    "builder. # + tab completion\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb19f345",
   "metadata": {},
   "source": [
    "### Setting *yambo*, *p2y* codes and resources\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b4a41a89",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"../configuration/codes_localhost.yaml\", 'r') as stream:\n",
    "    codes = yaml.safe_load(stream)\n",
    "    \n",
    "with open(\"../configuration/resources_localhost.yaml\", 'r') as stream:\n",
    "    resources = yaml.safe_load(stream)\n",
    "    \n",
    "builder.code = orm.load_code(codes[\"yambocode_id\"])                    # yambo executable\n",
    "builder.preprocessing_code = orm.load_code(codes[\"yamboprecode_id\"])   # p2y   executable\n",
    "builder.metadata.options = resources                                   # computational settings"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0a460e4",
   "metadata": {},
   "source": [
    "## The yambo parameters\n",
    "\n",
    "The input parameters used to run a yambo simulation (what usually is written in the *yambo.in* file for an ordinary yambo calculation), must be provided as an AiiDA Dict data type. This can be done first creating a python dictionary (e.g. `yambo_parameters`), and then converting it into the right AiiDA type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "68fa7ff4",
   "metadata": {},
   "outputs": [],
   "source": [
    "yambo_parameters = {\n",
    "    'arguments': [                                      # <=== runlevels\n",
    "        'dipoles',\n",
    "        'HF_and_locXC',\n",
    "        'dipoles',\n",
    "        'gw0',\n",
    "        'ppa',],\n",
    "    'variables': {                                      # <=== parameters\n",
    "        'Chimod': 'hartree',\n",
    "        'DysSolver': 'n',\n",
    "        'GTermKind': 'BG',\n",
    "        'NGsBlkXp': [2, 'Ry'],\n",
    "        'BndsRnXp': [[1, 50], ''],\n",
    "        'GbndRnge': [[1, 50], ''],\n",
    "        'QPkrange': [[[1, 10, 4,5]], ''],\n",
    "        },\n",
    "    }\n",
    "\n",
    "builder.parameters = orm.Dict(dict=yambo_parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10b01120",
   "metadata": {},
   "source": [
    "Runlevels are included under the 'arguments' key, as a list of strings. Parameters are instead in the 'variables' dictionary, as you can see. The format of the `yambo_parameters` is the same as the one that you should use in a *yambopy* calculation (indeed, to write the input file, the code uses the `YamboIn()` class of *yambopy*)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8134d4c",
   "metadata": {},
   "source": [
    "**Advanced settings**: You can also provide a settings dictionary, where to specify some instructions (False is default for any of them)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5602fb58",
   "metadata": {},
   "outputs": [],
   "source": [
    "settings = {\n",
    "    'INITIALISE':    False,      # if True, run only p2y and yambo initialization\n",
    "    'COPY_DBS':      False,      # if True, copy the out folder of a previous yambo calculation\n",
    "    'COPY_SAVE':     False,      # if True, copy the SAVE folder of a previous yambo calculation\n",
    "    'RESTART_YAMBO': False,      # if True, soft-link the out folder of a previous yambo calculation\n",
    "}\n",
    "\n",
    "builder.settings = orm.Dict(dict=settings)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acb0514c",
   "metadata": {},
   "source": [
    "Finally, to retrieve additional files, you can provide their names in the input settings Dict:\n",
    "\n",
    "```python\n",
    "settings_retrieve = {\n",
    "\t\"ADDITIONAL_RETRIEVE_LIST\":[\n",
    "\t\t'LOG/l-*',\n",
    "\t\t],\n",
    "\t}\n",
    "\n",
    "settings.update(settings_retrieve)\n",
    "builder.settings = orm.Dict(dict=settings)\n",
    "```\n",
    "\n",
    "This selects the additional files that will  be retrieved and parsed after a calculation is finished. Default retrieved files include the report files ``r-*``, text outputs ``o-*``, logs, the quasiparticle\n",
    "database for GW calculations ``aiida.out/ndb.QP``, and the Hartree-Fock and local exchange\n",
    "db ``aiida.out/ndb.HF_and_locXC``. Quantities\n",
    "that you may want collect to further analyse are for example the dipoles or the dielectric function databases, output of a standard GW\n",
    "calculation. In the above example we just selected all the log files."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc748413",
   "metadata": {},
   "source": [
    "### Setting the QuantumEspresso parent nscf calculation\n",
    "\n",
    "Yambo simulations require a DFT starting point, i.e. the ground state wavefunctions as obtained in a non self-consistent calculation. This was computed in the previous module of the tutorial ([*02_QE_starting_point*](../prerequisites/0_2_QE_starting_point.ipynb))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c1c6955d",
   "metadata": {},
   "outputs": [],
   "source": [
    "nscf_pk = 397 # the pk of the previous nscf calculation or workchain.\n",
    "\n",
    "builder.parent_folder = orm.load_node(nscf_pk).outputs.remote_folder"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08f91edc",
   "metadata": {},
   "source": [
    "## RUN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b7774365",
   "metadata": {},
   "outputs": [],
   "source": [
    "run = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "05394759",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "uuid: 23b71131-a31d-4219-9b87-ce9eb39cbbd5 (pk: 2053) (aiida.calculations:yambo.yambo)\n"
     ]
    }
   ],
   "source": [
    "if run:\n",
    "    print('A calculation is already running -> {}'.format(run.pk))\n",
    "    print('sure that you want to run again?, if so, copy the else instruction in the cell below and run!')\n",
    "else:\n",
    "    run = submit(builder)\n",
    "\n",
    "print(run)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b99e811",
   "metadata": {},
   "source": [
    "## Inspecting the outputs\n",
    "\n",
    "suppose that your calculation finished successfully:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "46173eeb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run.is_finished_ok"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d243637",
   "metadata": {},
   "source": [
    "Then you can access the outputs via tab completion: \n",
    "\n",
    "```python\n",
    "run.outputs. #+TAB\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "447459e6",
   "metadata": {},
   "source": [
    "For example, we have the `array_ndb` ArrayData, which essentially contains the information on the run.output.QP_DB (the parsed ndb.QP - which can also be accessed directly):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f0d8160f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['E',\n",
       " 'Z',\n",
       " 'Eo',\n",
       " 'Sc',\n",
       " 'Sx',\n",
       " 'Vxc',\n",
       " 'Band',\n",
       " 'Kpoint',\n",
       " 'nbands',\n",
       " 'nkpoints',\n",
       " 'qp_table',\n",
       " 'E_minus_Eo',\n",
       " 'Kpoint_index']"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run.outputs.array_ndb.get_arraynames()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e29d857e",
   "metadata": {},
   "source": [
    "The quasiparticle corrections and the renormalization factors can be accessed:\n",
    "\n",
    "```python\n",
    "energies_DFT = run.outputs.array_ndb.get_array('E_0')\n",
    "QP_corrections =  run.outputs.array_ndb.get_array('E_minus_Eo')\n",
    "energies_QP = run.outputs.array_ndb.get_array('E')\n",
    "Z_factors =  run.outputs.array_ndb.get_array('Z')\n",
    "kpoint_band_array = run.outputs.array_ndb.get_array('qp_table')\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "12946a15",
   "metadata": {},
   "outputs": [],
   "source": [
    "energies_DFT = run.outputs.array_ndb.get_array('Eo')\n",
    "QP_corrections =  run.outputs.array_ndb.get_array('E_minus_Eo')\n",
    "energies_QP = run.outputs.array_ndb.get_array('E')\n",
    "Z_factors =  run.outputs.array_ndb.get_array('Z')\n",
    "kpoint_band_array = run.outputs.array_ndb.get_array('qp_table')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd01adcc",
   "metadata": {},
   "source": [
    "The `energies_DFT` and `energies_QP` are in Hartree units (as well as the other related quantities). We can use the [*ASE*](https://wiki.fysik.dtu.dk/ase/) module `units` to convert easily to eV and plot the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "893490c8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGxCAYAAACa3EfLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAyIElEQVR4nO3df3SU5Z3//9ckhBmoyUhIQ0YTJMIHNKaiAaNRW4g0GNRUelbPclj8seDuMYLHH3v6tYt7TNItpXukaqs11VXBPZhV+7GocWs0bflh69KAkWpIwQ8YxcoExOhMhCYMyf39gyYS5x6SCTNzz33P83FOjt537uR++06Ul9d139flMgzDEAAAgM2lWV0AAABALBBqAACAIxBqAACAIxBqAACAIxBqAACAIxBqAACAIxBqAACAIxBqAACAI4yxuoBE6u/v1/79+5WZmSmXy2V1OQAAYAQMw1B3d7fOOOMMpaVFHo9JqVCzf/9+FRQUWF0GAAAYhY8++kj5+fkRP2/bULN69WqtXLlSd9xxhx566KERfU1mZqak403JysqKY3VSKBTS66+/rvnz5ysjIyOu97IT+mKOvoSjJ+boizn6Ys4pfQkGgyooKBj8czwSW4aabdu26fHHH9f5558f1dcNTDllZWUlJNSMHz9eWVlZtv5FijX6Yo6+hKMn5uiLOfpizml9Ge7REds9KPzFF1/oH/7hH/Sf//mfmjBhgtXlAACAJGG7kZrly5fr6quv1re//W398Ic/POm1vb296u3tHTwOBoOSjifXUCgU1zoHvn+872M39MUcfQlHT8zRF3P0xZxT+jLS+l2GYRhxriVmnn32Wa1atUrbtm2Tx+PR3LlzdcEFF0R8pqa2tlZ1dXVh5xsaGjR+/Pg4VwsAAGLhyJEjWrx4sQKBwEkfH7FNqPnoo480e/Zsvf7665o5c6YkDRtqzEZqCgoKdOjQoYQ8U9Pc3KyKigpHzGPGCn0xR1/C0RNz9MUcfTHnlL4Eg0Hl5OQMG2psM/301ltv6eDBg5o1a9bgub6+Pm3ZskWPPPKIent7lZ6ePuRr3G633G532PfKyMhI2A83kfeyE/pijr6Eoyfm6Is5+mLO7n0Zae22CTXz5s3Tu+++O+TcP/7jP+qcc87RPffcExZoAABAarFNqMnMzFRxcfGQc1/72tc0ceLEsPMAACD12O6VbgAAADO2Gakxs2nTJqtLAAAAScLWoQYAAFivr99QS0eXDnb3KDfTo9LCbKWnJX7jaEINAAAYtaY2v+oa2+UP9Aye83k9qqkqUmWxL6G18EwNAAAYlaY2v6rXtw4JNJLUGehR9fpWNbX5E1oPoQYAAEStr99QXWO7zFbwHThX19iuvv7ErfFLqAEAAFFr6egKG6E5kSHJH+hRS0dXwmoi1AAAgKgd7I4caEZzXSwQagAAQNRyMz0xvS4WCDUAACBqpYXZ8nk9ivTitkvH34IqLcxOWE2EGgAAELX0NJdqqookKSzYDBzXVBUldL0aQg0AABiVymKf6peUKM87dIopz+tR/ZKShK9Tw+J7AABg1CqLfaooymNFYQAAYH/paS6VTZ1odRlMPwEAAGcg1AAAAEcg1AAAAEcg1AAAAEcg1AAAAEcg1AAAAEcg1AAAAEdgnRoAAGyqr99IikXvkgWhBgAAG2pq86uusV3+QM/gOZ/Xo5qqooRvT5AsmH4CAMBmmtr8ql7fOiTQSFJnoEfV61vV1Oa3qDJrEWoAALCRvn5DdY3tMkw+N3CurrFdff1mVzgboQYAABtp6egKG6E5kSHJH+hRS0dX4opKEoQaAABs5GB35EAzmuuchFADAICN5GZ6YnqdkxBqAACwkdLCbPm8HkV6cdul429BlRZmJ7KspECoAQDARtLTXKqpKpKksGAzcFxTVZSS69UQagAAsJnKYp/ql5Qozzt0iinP61H9kpKUXaeGxfcAALChymKfKoryWFH4BIQaAABsKj3NpbKpE60uI2kw/QQAAByBUAMAAByBUAMAAByBUAMAAByBUAMAAByBUAMAABzBNqGmvr5e559/vrKyspSVlaWysjK9+uqrVpcFAACShG3WqcnPz9ePf/xjTZs2TZL09NNP69prr9Xbb7+t8847z+LqAAAYqq/fYGG8BLNNqKmqqhpyvGrVKtXX12vr1q2EGgBAUmlq86uusV3+QM/gOZ/Xo5qqopTdwiARbBNqTtTX16df/vKXOnz4sMrKyiJe19vbq97e3sHjYDAoSQqFQgqFQnGtceD7x/s+dkNfzNGXcPTEHH0xl0x9eW3nAd3+7J9kfOV8Z6BH1etb9fCimbryvEkJqSWZ+nIqRlq/yzCMr/Y9ab377rsqKytTT0+PTjvtNDU0NOiqq66KeH1tba3q6urCzjc0NGj8+PHxLBUAkIL6DamuNV2fH5XC99CWJEOnj5VqSvrETNTIHTlyRIsXL1YgEFBWVlbE62wVao4ePap9+/bp888/1wsvvKAnnnhCmzdvVlFRken1ZiM1BQUFOnTo0EmbEguhUEjNzc2qqKhQRkZGXO9lJ/TFHH0JR0/M0RdzydKXP3Z0aclT24e9bv3S2bq4MDvu9SRLX05VMBhUTk7OsKHGVtNPY8eOHXxQePbs2dq2bZt++tOf6rHHHjO93u12y+12h53PyMhI2A83kfeyE/pijr6Eoyfm6Is5q/vy6ZFjI74ukXVa3ZdTNdLabfNKtxnDMIaMxAAAYKXcTE9Mr0N0bDNSs3LlSi1YsEAFBQXq7u7Ws88+q02bNqmpqcnq0gAAkCSVFmbL5/WoM9AT9qCwdPwpmzzv8de7EXu2Gak5cOCAbrjhBs2YMUPz5s3TH//4RzU1NamiosLq0gAAkCSlp7lUU3X8Oc+vPgc8cFxTVcR6NXFim5GaJ5980uoSAAAYVmWxT/VLSsLWqcljnZq4s02oAQDALiqLfaooymNF4QQj1AAAEAfpaS6VTZ1odRkpxTbP1AAAAJwMIzUAAIgNKJ2AUAMASHlsQOkMTD8BAFJaU5tf1etbhwQa6csNKJva/BZVhmgRagAAKauv31BdY7vpQnkD5+oa29XXb5ttElMaoQYAkLJaOrrCRmhOZEjyB3rU0tGVuKIwaoQaAEDKOtgdOdCM5jpYi1ADAEhZbEDpLIQaAEDKGtiAMtKL2y4dfwuKDSjtgVADAEhZbEDpLIQaAEBKG9iAMs87dIopz+tR/ZIS1qmxERbfAwCkPDagdAZCDQAAYgNKJ2D6CQAAOAKhBgAAOAKhBgAAOALP1AAAHKOv39D2vZ/ysG+KItQAABzhT5+6tPonW9QZ7B085/N6VFNVxGvZKYLpJwCA7b2284Ceei9tSKCRpM5Aj6rXt6qpzW9RZUgkQg0AwNb6+g398Ne7TD9n/O2vdY3t6us3TK+BcxBqAAC21tLR9bcRGvNnZwxJ/kCPWjq6EloXEo9QAwCwtYPdPTG9DvZFqAEA2Fpupmf4i6K4DvZFqAEA2FppYbbystz68gmaoVw6/hZUaWF2QutC4hFqAAC2lp7m0r9ddY6k8KdqBo5rqopYryYFEGoAALZ35XmTtHR6vyZluYecz/N6VL+khHVqUgSL7wEAHGHmREP/3z98S2//pZsVhVMUoQYA4BjpaS6VTZ1odRmwCNNPAADAERipAQDEXV+/oZaOLqaFEFeEGgBAXDW1+VXX2C5/4MvF79hoEvFAqAEAxE1Tm1/V61vDVpDpDPTo1vWtuuvb/0dTcr7G6A1iglADAIiLvn5DdY3tpkviDZx78Df/b/Acozc4VTwoDACIi5aOriFTTsPpDPSoen2rmtr8cawKTkaoAQDERbQbSA6M3tQ1tquv33zLA+BkbBNqVq9erYsuukiZmZnKzc3VwoULtXv3bqvLAgBEMJoNJA1J/kCPWjq6Yl8QHM82oWbz5s1avny5tm7dqubmZh07dkzz58/X4cOHrS4NAGCitDBbPq8nbD+mkYh2lAeQbPSgcFNT05DjtWvXKjc3V2+99Za+9a1vWVQVACCS9DSXaqqKVL2+VS5F2kPb3GhGeQDbhJqvCgQCkqTs7Mhbyff29qq3t3fwOBgMSpJCoZBCoVBc6xv4/vG+j93QF3P0JRw9MWe3vsybkaOHF83UD3+9S53B3mGvd0nK87p1YX5mVP+MdutLojilLyOt32UYhu2exjIMQ9dee60+++wzvfHGGxGvq62tVV1dXdj5hoYGjR8/Pp4lAgBO0G9Ie4MuBUPSJ3+VXv3LwNMPJ05OHf/jaOn0fs2caLs/mhBHR44c0eLFixUIBJSVlRXxOluGmuXLl+t//ud/9Pvf/175+fkRrzMbqSkoKNChQ4dO2pRYCIVCam5uVkVFhTIyMuJ6LzuhL+boSzh6Ys4pfXlt54Gw0Ruf1617F5yjK8+bFPX3c0pfYs0pfQkGg8rJyRk21Nhu+un222/Xyy+/rC1btpw00EiS2+2W2+0OO5+RkZGwH24i72Un9MUcfQlHT8zZvS/XXJCvBeefGfP9oOzel3ixe19GWrttQo1hGLr99tu1YcMGbdq0SYWFhVaXBAA4BelpLpVNnWh1GXAQ24Sa5cuXq6GhQS+99JIyMzPV2dkpSfJ6vRo3bpzF1QEAAKvZZp2a+vp6BQIBzZ07Vz6fb/Djueees7o0AACQBGwzUmPD55kBAEAC2WakBgAA4GRsM1IDALBGX78R87eUgHgg1AAAImpq86uusV3+wJd7Mfm8HtVUFamy2GdhZUA4pp8AAKaa2vyqXt86JNBIUmegR9XrW9XU5reoMsAcoQYAEKav31BdY7vpJpQD5+oa29XXz0scSB6EGgBAmJaOrrARmhMZkvyBHrV0dCWuKGAYhBoAQJiD3ZEDzWiuAxKBUAMACJOb6YnpdUAiEGoAAGFKC7Pl83oU6cVtl46/BVVamJ3IsoCTItQAAMKkp7lUU1UkSWHBZuC4pqqI9WqQVAg1AABTlcU+1S8pUZ536BRTntej+iUlrFODpMPiewCAiCqLfaooymNFYdgCoQYAcFLpaS6VTZ1odRnAsJh+AgAAjkCoAQAAjkCoAQAAjkCoAQAAjkCoAQAAjkCoAQAAjkCoAQAAjkCoAQAAjkCoAQAAjkCoAQAAjkCoAQAAjsDeTwBgA339hlo6uuT//LDeD7jU128ow+qigCRDqAGAJNfU5lddY7v8gZ6/nUnX//3JFtV+5zxVFvssrQ1IJkw/AUASa2rzq3p96wmB5rgDwV5Vr29VU5vfosqA5EOoAYAk1ddvqK6xXYbJ5wbO1TW2q6/f7Aog9RBqACBJtXR0hY3QnMiQ5A/0qKWjK3FFAUmMUAMASepgd+RAM5rrAKcj1ABAksrN9MT0OsDpCDUAkKRKC7Pl83rkivB5lySf16PSwuxElgUkLUINACSp9DSXaqqKJCks2Awc11QVKT0tUuwBUguhBgCSWGWxT/VLSpTnHTrFlOd1q35JCevUACdg8T0ASHKVxT5VFOV9uaLwzh1a8fffksc91urSgKRCqAEAG0hPc6ls6kSFQln69V/eZsoJMMH0EwAAcARbhZotW7aoqqpKZ5xxhlwul1588UWrSwIAAEnCVqHm8OHDmjlzph555BGrSwEAAEnGVs/ULFiwQAsWLLC6DAAAkIRsFWqi1dvbq97e3sHjYDAoSQqFQgqFQnG998D3j/d97Ia+mKMv4eiJOfpijr6Yc0pfRlq/yzAMW27v6nK5tGHDBi1cuDDiNbW1taqrqws739DQoPHjx8exOgAAECtHjhzR4sWLFQgElJWVFfE6R4cas5GagoICHTp06KRNiYVQKKTm5mZVVFQoIyMjrveyE/pijr6Eoyfm6Is5+mLOKX0JBoPKyckZNtQ4evrJ7XbL7XaHnc/IyEjYDzeR97IT+mKOvoSjJ+boizn6Ys7ufRlp7bZ6+wkAACASW43UfPHFF9qzZ8/gcUdHh3bs2KHs7GxNnjzZwsoAAIDVbBVqtm/frvLy8sHju+++W5J00003ad26dRZVBQAAkoGtQs3cuXNl0+eaAdhcX7+hlo4uHezuUW6mR6WF2ey/BCQZW4UaALBCU5tfdY3t8gd6Bs/5vB7VVBWpsthnYWUATsSDwgBwEk1tflWvbx0SaCSpM9Cj6vWtamrzW1QZgK8i1ABABH39huoa22U26T1wrq6xXX39TIsDyYBQAwARtHR0hY3QnMiQ5A/0qKWjK3FFAYiIUAMAERzsjhxoRnMdgPgi1ABABLmZnpheByC+CDUAEEFpYbZ8Xo8ivbjt0vG3oEoLsxNZFoAICDUAEEF6mks1VUWSFBZsBo5rqopYrwZIEoQaADiJymKf6peUKM87dIopz+tR/ZIS1qkBkgiL7wHAMCqLfaooymNFYSDJEWoAYATS01wqmzrR6jIAnMSoQ82HH36ozs5OuVwuTZo0SWeddVYs6wIAAIhK1M/UPPjggyooKNDZZ5+tsrIyXXLJJTr77LNVUFCghx56KA4lAgAADC+qkZp///d/15o1a7Ry5UpdeeWVmjRpkgzD0MGDB/Xaa6+ptrZWX3zxhf7t3/4tXvUCAACYiirUPP7443r66ae1cOHCIefPOOMMXXDBBZo+fbpWrFhBqAEAAAkX1fTTp59+qhkzZkT8/PTp0/XZZ5+dclEAAADRiirUlJaWatWqVTp27FjY544dO6Yf/ehHKi0tjVlxAAAAIxXV9NPDDz+s+fPnKzc3V3PmzNGkSZPkcrnU2dmpLVu2yO12q7m5OV61AgAARBTVSM03vvENvffee1q1apWysrLU0dGh999/X1lZWVq1apV27dql8847L161AgAARBT1OjWZmZmqrq5WdXV1POoBAAAYlahCzeWXX64rrrhC5eXlKisrk8fjGf6LAAAAEiCq6acZM2aooaFB8+bN04QJEzR37lz94Ac/0BtvvKFQKBSvGgEAAIYVVah58skntWfPHu3bt0+PPfaYpk6dqqefflpz5szR6aefroqKCq1evTpetQIAAEQU9TYJkpSfn68bb7xRTz75pPbu3asPP/xQd911l1paWlh4DwAAWGLUG1ru3btXmzZtGvz4/PPPVVZWpjlz5sSyPgAAgBGJKtSsXbtWGzdu1KZNmxQIBHTZZZdpzpw5Wr58uWbPnq0xY0adkQAAAE5JVClk2bJlmjx5su69914tXbpUGRkZ8aoLAAAgKlE9U/Pzn/9cl1xyiWpra5Wbm6uqqir95Cc/0fbt22UYRrxqBAAAGFZUoaa6ulrPPvus/H6//vCHP+iqq65SS0uLrrnmGk2YMEFXX3211qxZE69aAQAAIhrV20+SVFRUpOrqaj333HNqbW3VihUr9Pvf/1733HNPLOsDAAAYkVE92XvgwIEhbz699957Gjt2rC6++GKVl5fHukYAAIBhRRVqli9fro0bN2r37t0aM2aMLrroIl133XUqLy/XpZdeyrYJAADAMlGFmtbWVi1cuFDl5eW67LLLNH78+HjVBQAAEJWoQs3//u//xqsOAACAUxL1MzWGYeg3v/mN3nzzTXV2dsrlcmnSpEm67LLLNG/ePLlcrnjUCQAAcFJRvf308ccfq6SkRAsWLNCGDRv0/vvva8+ePdqwYYMqKys1e/Zsffzxx/GqFQAAIKKoRmpuu+02ZWdn66OPPpLP5xvyOb/fryVLlmj58uV68cUXY1kjAADAsKIaqfntb3+rBx54ICzQSJLP59OaNWv0m9/8JmbFmXn00UdVWFgoj8ejWbNm6Y033ojr/QAkXl+/oT92dOmtQy79saNLff2sWA5geFGN1IwbN05dXV0RP//ZZ59p3Lhxp1xUJM8995zuvPNOPfroo7rsssv02GOPacGCBWpvb9fkyZPjdl8AidPU5lddY7v8gR5J6fqv/7ddPq9HNVVFqiwO/x8qABgQVahZtGiRbrrpJj3wwAOqqKiQ1+uVJAUCATU3N+tf/uVftHjx4rgUKkkPPPCAli1bpltuuUWS9NBDD+m1115TfX29Vq9eHXZ9b2+vent7B4+DwaAkKRQKKRQKxa3OgXuc+FccR1/M0ZfjXtt5QLc/+yd9dVymM9Cj6vWtenjRTF153iRLaksW/K6Yoy/mnNKXkdbvMqLYifLo0aO644479NRTT+nYsWMaO3bs4PkxY8Zo2bJlevDBB+V2u0dX9TD3Hj9+vH75y1/qu9/97uD5O+64Qzt27NDmzZvDvqa2tlZ1dXVh5xsaGlhjB0gy/YZU15quz49KktlblIZOHyvVlPQpjZcsgZRy5MgRLV68WIFAQFlZWRGviyrUDAgGg3rrrbfU2dkpScrLy9OsWbNOeqNTtX//fp155pn6wx/+oEsvvXTw/I9+9CM9/fTT2r17d9jXmI3UFBQU6NChQ3GtVTqeKpubm1VRUaGMjIy43stO6Is5+iL9saNLS57aPux165fO1sWF2QmoKDnxu2KOvphzSl+CwaBycnKGDTVRTT/99a9/1W9/+1tdc801Ki8v17/+678OhobGxkaNGTNGP/jBD+K6XcJX18ExDCPi2jhut9t01CgjIyNhP9xE3stO6Iu5VO7Lp0eOjfi6VO3RiVL5d+Vk6Is5u/dlpLVHFWr+67/+S6+88oquueYaSdIjjzyi8847b/Dh4F27dsnn8+muu+6Kstzh5eTkKD09fXB0aMDBgwc1aVJqz7EDTpCbObL/GRrpdQBST1SvdD/zzDNaunTpkHMNDQ3auHGjNm7cqPvvv1/PP/98TAscMHbsWM2aNUvNzc1Dzjc3Nw+ZjgJgT6WF2fJ5PaZP00jHn7LxeT0qTeGpJwAnF1Woee+99zR9+vTBY4/Ho7S0L79FaWmp2tvbY1fdV9x999164okn9NRTT+nPf/6z7rrrLu3bt0+33npr3O4JIDHS01yqqSqSFP6Y8MBxTVWR0nlKGEAEUU0/BQIBjRnz5Zd88sknQz7f398/5MHcWPv7v/97ffrpp/rBD34gv9+v4uJi/frXv9ZZZ50Vt3sCSJzKYp/ql5ScsE7NcXmsUwNgBKIKNfn5+Wpra9OMGTNMP//OO+8oPz8/JoVFctttt+m2226L6z0AWKey2KeKojz9756Dev2NP2r+Ny9W2bRcRmgADCuq6aerrrpK9913n3p6esI+99e//lV1dXW6+uqrY1YcgNSUnubSxYXZmpVj6OLCbAINgBGJaqRm5cqVev755zVjxgytWLFC06dPl8vl0q5du/TII4/o2LFjWrlyZbxqBQAAiCiqUDNp0iS9+eabqq6u1ve//30NrNvncrlUUVGhRx99lNerAQCAJaIKNZJUWFiopqYmdXV1ac+ePZKkadOmKTub1ywBAIB1og41A7Kzs1VaWhrLWgAAAEYtqgeFAQAAkhWhBgAAOAKhBgAAOAKhBgAAOAKhBgAAOAKhBgAAOAKhBgAAOAKhBgAAOAKhBgAAOAKhBgAAOAKhBgAAOAKhBgAAOAKhBgAAOAKhBgAAOMIYqwsAkFh9/YZaOrp0sLtHuZkelRZmKz3NZXVZAHDKCDVACmlq86uusV3+QM/gOZ/Xo5qqIlUW+yysDABOHdNPQIpoavOren3rkEAjSZ2BHlWvb1VTm9+iygAgNgg1QAro6zdU19guw+RzA+fqGtvV1292BQDYA6EGSAEtHV1hIzQnMiT5Az1q6ehKXFEAEGOEGiAFHOyOHGhGcx0AJCNCDZACcjM9Mb0OAJIRoQZIAaWF2fJ5PYr04rZLx9+CKi3MTmRZABBThBogBaSnuVRTVSRJYcFm4Limqoj1agDYGqEGSBGVxT7VLylRnnfoFFOe16P6JSWsUwPA9lh8D0ghlcU+VRTlsaIwAEci1AApJj3NpbKpE60uAwBijuknAADgCIzUADbAJpQAMDxCDZDk2IQSAEaG6ScgibEJJQCMHKEGSFJsQgkA0bFNqFm1apUuvfRSjR8/XqeffrrV5QBxxyaUABAd24Sao0eP6vrrr1d1dbXVpQAJwSaUABAd2zwoXFdXJ0lat26dtYUACcImlAAQHduEmtHo7e1Vb2/v4HEwGJQkhUIhhUKhuN574PvH+z52Q1/MmfXlwvxM5WW5dSDYa/pcjUtSntetC/MzHdlPflfM0Rdz9MWcU/oy0vpdhmHY6inDdevW6c4779Tnn38+7LW1tbWDIzwnamho0Pjx4+NQHRBbf/rUpafeG5glPnFdmuP/2i6d3q+ZE231rzAARO3IkSNavHixAoGAsrKyIl5naaiJFDpOtG3bNs2ePXvwOJpQYzZSU1BQoEOHDp20KbEQCoXU3NysiooKZWRkxPVedkJfzJ2sL6/tPKAf/nqXOoNf/i77vG7du+AcXXnepESXmjD8rpijL+boizmn9CUYDConJ2fYUGPp9NOKFSu0aNGik14zZcqUUX9/t9stt9sddj4jIyNhP9xE3stO6Is5s75cc0G+Fpx/ZsquKMzvijn6Yo6+mLN7X0Zau6WhJicnRzk5OVaWANgCm1ACwPBs86Dwvn371NXVpX379qmvr087duyQJE2bNk2nnXaatcUBAADL2SbU3HfffXr66acHjy+88EJJ0saNGzV37lyLqgIAAMnCNovvrVu3ToZhhH0QaAAAgGSjkRrATvr6jZR9sBcArEKoAWKsqc2vusb2Ifs2+bwe1VQVqbLYZ2FlAOBstpl+Auygqc2v6vWtYRtRdgZ6VL2+VU1tfosqAwDnI9QAMdLXb6iusd10S4OBc3WN7errZwVgAIgHQg0QIy0dXWEjNCcyJPkDPWrp6EpcUQCQQgg1QIwc7I4caEZzHQAgOoQaIEZyMz0xvQ4AEB1CDRAjpYXZ8nk9ivTitkvH34IqLcxOZFkAkDIINUCMpKe5VFNVJElhwWbguKaqiPVqACBOCDVADFUW+1S/pER53qFTTHlej+qXlLBODQDEEYvvATFWWexTRVEeKwoDQIIRaoA4SE9zqWzqRKvLAICUwvQTAABwBEINAABwBEINAABwBEINAABwBEINAABwBEINAABwBEINAABwBNapQUrp6zdYFA8AHIpQg5TR1OZXXWO7/IGewXM+r0c1VUVsXwAADsD0E1JCU5tf1etbhwQaSeoM9Kh6faua2vwWVQYAiBVCDRyvr99QXWO7DJPPDZyra2xXX7/ZFQAAuyDUwPFaOrrCRmhOZEjyB3rU0tGVuKIAADFHqIHjHeyOHGhGcx0AIDkRauB4uZmemF4HAEhOhBo4XmlhtnxejyK9uO3S8begSguzE1kWACDGCDVwvPQ0l2qqiiQpLNgMHNdUFbFeDQDYHKEGKaGy2Kf6JSXK8w6dYsrzelS/pIR1agDAAVh8DymjstiniqI8VhQGAIci1MBWTnWbg/Q0l8qmToxjhQAAqxBqYBtscwAAOBmeqYEtsM0BAGA4hBokPbY5AACMBKEGSY9tDgAAI0GoQdJjmwMAwEgQapD02OYAADAStgg1H3zwgZYtW6bCwkKNGzdOU6dOVU1NjY4ePWp1aUgAtjkAAIyELULNrl271N/fr8cee0w7d+7Ugw8+qF/84hdauXKl1aUhAdjmAAAwErZYp6ayslKVlZWDx2effbZ2796t+vp6rVmzxsLKkCgD2xx8dZ2aPNapAQD8jS1CjZlAIKDs7JNPN/T29qq3t3fwOBgMSpJCoZBCoVBc6xv4/vG+j92cSl/mzcjR3P/zTW3/8DMd7O5VbqZbs8+aoPQ0l+37zO9LOHpijr6Yoy/mnNKXkdbvMgzDdot77N27VyUlJfrJT36iW265JeJ1tbW1qqurCzvf0NCg8ePHx7NEAAAQI0eOHNHixYsVCASUlZUV8TpLQ02k0HGibdu2afbs2YPH+/fv15w5czRnzhw98cQTJ/1as5GagoICHTp06KRNiYVQKKTm5mZVVFQoIyMjrveyE/pijr6Eoyfm6Is5+mLOKX0JBoPKyckZNtRYOv20YsUKLVq06KTXTJkyZfDv9+/fr/LycpWVlenxxx8f9vu73W653e6w8xkZGQn74SbyXnZCX8zRl3D0xBx9MUdfzNm9LyOt3dJQk5OTo5ycnBFd+/HHH6u8vFyzZs3S2rVrlZZmixe3AABAgtjiQeH9+/dr7ty5mjx5stasWaNPPvlk8HN5eXkWVobh9PUbauno0sHuHuVmenRhfqbVJQEAHMoWoeb111/Xnj17tGfPHuXn5w/5nA2fc04ZTW3+8Fews9y6Ks+lqyysCwDgTLaYw7n55ptlGIbpB5JTU5tf1etbwzaiPBDs1VPvpem1nQcsqgwA4FS2CDWwl75+Q3WN7TKLnAPnVr26S339hFIAQOwQahBzLR1dYSM0Q7nkD/SqpaMrYTUBAJyPUIOYO9h9skAT/XUAAIwEoQYxl5vpiel1AACMBKEGMVdamC2f1xO2o/aXDPm8bpUWnnzvLgAAokGoQcylp7lUU1UkSWHBZuD43gXnKD0tcuwBACBahBrERWWxT/VLSpTnHTrFlOd1a+n0fl153iSLKgMAOJUtFt+DPVUW+1RRlBe2ovBrTa9aXRoAwIEINYir9DSXyqZOHDwOhUIWVgMAcDKmnwAAgCMQagAAgCMQagAAgCMQagAAgCMQagAAgCMQagAAgCMQagAAgCMQagAAgCMQagAAgCMQagAAgCMQagAAgCMQagAAgCMQagAAgCMQagAAgCMQagAAgCMQagAAgCMQagAAgCMQagAAgCOMsboAjF5fv6GWji4d7O5RbqZHpYXZSk9zWV0WAACWINTYVFObX3WN7fIHegbP+bwe1VQVqbLYZ2FlAABYg+knG2pq86t6feuQQCNJnYEeVa9vVVOb36LKAACwDqHGZvr6DdU1tssw+dzAubrGdvX1m10BAIBzEWpspqWjK2yE5kSGJH+gRy0dXYkrCgCAJECosZmD3ZEDzWiuAwDAKXhQOElFerMpN9Mzoq8f6XUAADgFoSYJnezNpoqiPPm8HnUGekyfq3FJyvMeD0EAAKQSpp+SzHBvNjW3d6qmqkjS8QBzooHjmqoi1qsBAKQcQk0SGembTRVFeapfUqI879AppjyvR/VLSlinBgCQkmwz/fSd73xHO3bs0MGDBzVhwgR9+9vf1n/8x3/ojDPOsLq0mInmzabKYp8qivJYURgAgL+xzUhNeXm5nn/+ee3evVsvvPCC9u7dq+uuu87qsmIq2jeb0tNcKps6UddecKbKpk4k0AAAUpptRmruuuuuwb8/66yz9P3vf18LFy5UKBRSRkaGhZXFDm82AQAwerYJNSfq6urSM888o0svvfSkgaa3t1e9vb2Dx8FgUJIUCoUUCoXiWuPA94/mPhfmZyovy60Dwd6TvNnk1oX5mXGvP15G05dUQF/C0RNz9MUcfTHnlL6MtH6XYRi2WU//nnvu0SOPPKIjR47okksu0SuvvKKJEydGvL62tlZ1dXVh5xsaGjR+/Ph4ljpqf/rUpafeG5gVPHE66fiPaen0fs2caJsfGQAAp+zIkSNavHixAoGAsrKyIl5naaiJFDpOtG3bNs2ePVuSdOjQIXV1denDDz9UXV2dvF6vXnnlFblc5s+SmI3UFBQU6NChQydtSiyEQiE1NzeroqIi6umx13Ye0A9/vUudwS9r93ndunfBObryvEmxLjWhTqUvTkZfwtETc/TFHH0x55S+BINB5eTkDBtqLJ1+WrFihRYtWnTSa6ZMmTL49zk5OcrJydH06dN17rnnqqCgQFu3blVZWZnp17rdbrnd7rDzGRkZCfvhjuZe11yQrwXnn+noN5sS+TOwE/oSjp6Yoy/m6Is5u/dlpLVbGmoGQspoDAwwnTgS4yQDbzYBAICRscWDwi0tLWppadHll1+uCRMm6P3339d9992nqVOnRhylAQAAqcUW69SMGzdOv/rVrzRv3jzNmDFDS5cuVXFxsTZv3mw6vQQAAFKPLUZqvvGNb+h3v/ud1WUAAIAkZouRGgAAgOHYYqQmmfX1G45+SwkAALsg1JyCpja/6hrbh2xC6fN6VFNVpHkzRvdWFwAAGB2mn0apqc2v6vWtYbtqdwZ6VL2+Va/tPGBRZQAApCZCzSj09Ruqa2w33Z9p4NyqV3epn90MAABIGELNKLR0dIWN0JzIkOQP9GpvkGdrAABIFELNKBzsjhxoThS096aoAADYCqFmFHIzPSO6Lsu+22wAAGA7hJpRKC3Mls/rUaTJJZeO76o9NYuHagAASBRCzSikp7lUU1UkSWHBZuD43gXniOVqAABIHELNKFUW+1S/pER53qFTUXlej+qXlOjK8yZZVBkAAKmJxfdOQWWxTxVFeaYrCodCPCUMAEAiEWpOUXqaS2VTJ1pdBgAAKY/pJwAA4AiEGgAA4AiEGgAA4AiEGgAA4AiEGgAA4AiEGgAA4AiEGgAA4AiEGgAA4AiEGgAA4AgptaKwYRzfNTsYDMb9XqFQSEeOHFEwGFRGRkbc72cX9MUcfQlHT8zRF3P0xZxT+jLw5/bAn+ORpFSo6e7uliQVFBRYXAkAAIhWd3e3vF5vxM+7jOFij4P09/dr//79yszMlMvliuu9gsGgCgoK9NFHHykrKyuu97IT+mKOvoSjJ+boizn6Ys4pfTEMQ93d3TrjjDOUlhb5yZmUGqlJS0tTfn5+Qu+ZlZVl61+keKEv5uhLOHpijr6Yoy/mnNCXk43QDOBBYQAA4AiEGgAA4AiEmjhxu92qqamR2+22upSkQl/M0Zdw9MQcfTFHX8ylWl9S6kFhAADgXIzUAAAARyDUAAAARyDUAAAARyDUAAAARyDUJMB3vvMdTZ48WR6PRz6fTzfccIP2799vdVmW+uCDD7Rs2TIVFhZq3Lhxmjp1qmpqanT06FGrS7PcqlWrdOmll2r8+PE6/fTTrS7HMo8++qgKCwvl8Xg0a9YsvfHGG1aXZKktW7aoqqpKZ5xxhlwul1588UWrS0oKq1ev1kUXXaTMzEzl5uZq4cKF2r17t9VlWaq+vl7nn3/+4IJ7ZWVlevXVV60uKyEINQlQXl6u559/Xrt379YLL7ygvXv36rrrrrO6LEvt2rVL/f39euyxx7Rz5049+OCD+sUvfqGVK1daXZrljh49quuvv17V1dVWl2KZ5557Tnfeeafuvfdevf322/rmN7+pBQsWaN++fVaXZpnDhw9r5syZeuSRR6wuJals3rxZy5cv19atW9Xc3Kxjx45p/vz5Onz4sNWlWSY/P18//vGPtX37dm3fvl1XXHGFrr32Wu3cudPq0uKOV7ot8PLLL2vhwoXq7e219a6psXb//fervr5e77//vtWlJIV169bpzjvv1Oeff251KQl38cUXq6SkRPX19YPnzj33XC1cuFCrV6+2sLLk4HK5tGHDBi1cuNDqUpLOJ598otzcXG3evFnf+ta3rC4naWRnZ+v+++/XsmXLrC4lrhipSbCuri4988wzuvTSSwk0XxEIBJSdnW11GbDY0aNH9dZbb2n+/PlDzs+fP19vvvmmRVXBLgKBgCTx35K/6evr07PPPqvDhw+rrKzM6nLijlCTIPfcc4++9rWvaeLEidq3b59eeuklq0tKKnv37tXDDz+sW2+91epSYLFDhw6pr69PkyZNGnJ+0qRJ6uzstKgq2IFhGLr77rt1+eWXq7i42OpyLPXuu+/qtNNOk9vt1q233qoNGzaoqKjI6rLijlAzSrW1tXK5XCf92L59++D13/ve9/T222/r9ddfV3p6um688UY5ceYv2r5I0v79+1VZWanrr79et9xyi0WVx9do+pLqXC7XkGPDMMLOASdasWKF3nnnHf33f/+31aVYbsaMGdqxY4e2bt2q6upq3XTTTWpvb7e6rLgbY3UBdrVixQotWrTopNdMmTJl8O9zcnKUk5Oj6dOn69xzz1VBQYG2bt3quOHAaPuyf/9+lZeXq6ysTI8//nicq7NOtH1JZTk5OUpPTw8blTl48GDY6A0w4Pbbb9fLL7+sLVu2KD8/3+pyLDd27FhNmzZNkjR79mxt27ZNP/3pT/XYY49ZXFl8EWpGaSCkjMbACE1vb28sS0oK0fTl448/Vnl5uWbNmqW1a9cqLc25A4en8vuSasaOHatZs2apublZ3/3udwfPNzc369prr7WwMiQjwzB0++23a8OGDdq0aZMKCwutLikpGYbhyD9zvopQE2ctLS1qaWnR5ZdfrgkTJuj999/Xfffdp6lTpzpulCYa+/fv19y5czV58mStWbNGn3zyyeDn8vLyLKzMevv27VNXV5f27dunvr4+7dixQ5I0bdo0nXbaadYWlyB33323brjhBs2ePXtwFG/fvn0p/czVF198oT179gwed3R0aMeOHcrOztbkyZMtrMxay5cvV0NDg1566SVlZmYOjvB5vV6NGzfO4uqssXLlSi1YsEAFBQXq7u7Ws88+q02bNqmpqcnq0uLPQFy98847Rnl5uZGdnW243W5jypQpxq233mr85S9/sbo0S61du9aQZPqR6m666SbTvmzcuNHq0hLq5z//uXHWWWcZY8eONUpKSozNmzdbXZKlNm7caPp7cdNNN1ldmqUi/Xdk7dq1VpdmmaVLlw7+u/P1r3/dmDdvnvH6669bXVZCsE4NAABwBOc+xAAAAFIKoQYAADgCoQYAADgCoQYAADgCoQYAADgCoQYAADgCoQYAADgCoQYAADgCoQYAADgCoQZA0rr55pvlcrnkcrmUkZGhSZMmqaKiQk899ZT6+/sHr5syZcrgdQMf+fn5qq2tDTv/1Y8PPvjAun9AADHFNgkAktbNN9+sAwcOaO3aterr69OBAwfU1NSk1atX65vf/KZefvlljRkzRlOmTNGyZcv0T//0T4Nfm56ernHjxumLL74YPHfRRRfpn//5n4dc9/Wvf13p6ekJ/ecCEB/s0g0gqbnd7sGd288880yVlJTokksu0bx587Ru3TrdcsstkqTMzEzTHd5P3Nk8PT094nUA7I/pJwC2c8UVV2jmzJn61a9+ZXUpAJIIoQaALZ1zzjlDnoe55557dNpppw1+/OxnP7OuOACWYPoJgC0ZhiGXyzV4/L3vfU8333zz4HFOTo4FVQGwEqEGgC39+c9/VmFh4eBxTk6Opk2bZmFFAKzG9BMA2/nd736nd999V3/3d39ndSkAkggjNQCSWm9vrzo7O8Ne6b7mmmt04403Wl0egCRCqAGQ1JqamuTz+TRmzBhNmDBBM2fO1M9+9jPddNNNSktjsBnAl1h8DwAAOAL/mwMAAByBUAMAAByBUAMAAByBUAMAAByBUAMAAByBUAMAAByBUAMAAByBUAMAAByBUAMAAByBUAMAAByBUAMAABzh/wfZZukGvVlEQAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ase import units\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "plt.plot(energies_DFT.real*units.Ha,energies_QP.real*units.Ha,'o')\n",
    "plt.ylabel('G0W0')\n",
    "plt.xlabel('DFT')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "93308bf9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.001109  +0.00029073j,  0.10830379-0.00055046j,\n",
       "       -0.02391809+0.00031545j,  0.09867941-0.00054466j,\n",
       "       -0.04554427+0.00036696j,  0.07676777-0.00051919j,\n",
       "       -0.06044236+0.00037155j,  0.06342715-0.00049692j,\n",
       "       -0.08035547+0.00042069j,  0.07538369-0.00051789j,\n",
       "       -0.06892121+0.00040126j,  0.08728357-0.0005428j ,\n",
       "       -0.03420922+0.00033798j,  0.12298796-0.00059331j,\n",
       "       -0.10801804+0.00049053j,  0.02807045-0.00046628j,\n",
       "       -0.10317245+0.00047007j,  0.04541743-0.00049563j,\n",
       "       -0.12011653+0.00046856j,  0.13777733-0.00060036j], dtype=complex64)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "energies_QP"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eeb8aa19",
   "metadata": {},
   "source": [
    "The `output_parameters` Dict contains some information on the simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "91a49104",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'errors': [],\n",
       " 'timing': [' <---> [01] MPI/OPENMP structure, Files & I/O Directories\\n',\n",
       "  ' <---> [02] CORE Variables Setup\\n',\n",
       "  ' <---> [03] Transferred momenta grid and indexing\\n',\n",
       "  ' <---> [04] Dipoles\\n',\n",
       "  ' <---> [05] Dynamic Dielectric Matrix (PPA)\\n',\n",
       "  ' <38s> [06] Local Exchange-Correlation + Non-Local Fock\\n',\n",
       "  ' <52s> [07] Dyson equation: Newton solver\\n',\n",
       "  ' <01m-18s> [08] Timing Overview\\n',\n",
       "  ' <01m-18s> [09] Game Over & Game summary\\n'],\n",
       " 'has_gpu': False,\n",
       " 'memstats': [],\n",
       " 'warnings': [],\n",
       " 'Fermi(eV)': 6.467203,\n",
       " 'game_over': True,\n",
       " 'last_time': 78,\n",
       " 'para_error': False,\n",
       " 'time_error': False,\n",
       " 'ns_db1_path': '/tmp/tmp9w66h637',\n",
       " 'memory_error': False,\n",
       " 'p2y_completed': True,\n",
       " 'yambo_version': '5.x',\n",
       " 'requested_time': 3600,\n",
       " 'last_time_units': 'seconds',\n",
       " 'yambo_wrote_dbs': True,\n",
       " 'corrupted_fragment': False,\n",
       " 'X_par_allocation_error': False}"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run.outputs.output_parameters.get_dict()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f5a3fab",
   "metadata": {},
   "source": [
    "As output, also the `ndb.QP` database is stored:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "ed63f08f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<SinglefileData: uuid: 2b8fc5a9-74dc-4ac3-8f71-8cc339fe252f (pk: 2056)>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run.outputs.QP_db"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6e9ec2a",
   "metadata": {},
   "source": [
    "This can be used again: for example, we may be interested on submitting BSE simulations on top of these quasiparticle corrections. \n",
    "We can do the following (assuming the `builder.parameters` are now the one for a BSE calculation):\n",
    "\n",
    "```python\n",
    "    builder.QP_corrections = orm.load_node(2032)   # ---> QP_db node\n",
    "    params['variables']['KfnQPdb'] =  \"E < ./ndb.QP\"\n",
    "    builder.parameters = orm.Dict(dict=params)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a87f1f1",
   "metadata": {},
   "source": [
    "### How to access remote files.\n",
    "\n",
    "You may need to go in the folder where all inputs and outputs are collected. \n",
    "This can be easily done via the command line:\n",
    "\n",
    "```bash\n",
    "verdi calcjob gotocomputer <pk of the calcjob>\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cc046a4",
   "metadata": {},
   "source": [
    "## BSE simulations\n",
    "\n",
    "It is possible also to run BSE calculations, - to compute optical properties - \n",
    "as shown in the following sections of this tutorial.\n",
    "\n",
    "If you want to run a single G0W0/BSE simulation within the aiida-yambo plugin, we suggest to always use the YamboWorkflow workchain provided in\n",
    "the plugin and as shown in the following sections."
   ]
  }
 ],
 "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.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}