notebooks/hcp_workfunction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Workfunction of hcp (0001) surfaces\n",
    "\n",
    "In this notebook, we will show how to calculate the workfunction of selected hcp(0001) surfaces using [VASP](https://www.vasp.at/). Please keep in mind that the parameters used here give no converged results. They have been chosen to demonstrate the workflow using inexpensive calculations. For converged results, parameters such as lattice parameters, plane-wave energy cutoffs, reciprocal space sampling or the need to perform spin polarized calculations have to be carefully chosen "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "%matplotlib inline \n",
    "import matplotlib.pylab as plt\n",
    "import pandas as pd\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyiron.project import Project"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "pr = Project(\"hcp_workfunction\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating the Workfunction of Mg(0001) \n",
    "\n",
    "### Structure creation\n",
    "We use the `create_surface()` function which uses the ASE surface generator to build our surface slab structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "99239e67eb69452d8bda5e262b3fb7fb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>NGLWidget</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "NGLWidget()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    " # Now we set-up the Mg (0001) surface\n",
    "a = 3.1919\n",
    "c = 5.1852\n",
    "\n",
    "# Vacuum region to break the periodicity along the z-axis\n",
    "vac = 10\n",
    "size = (2, 2, 4)\n",
    "Mg_0001 = pr.create_surface(\"Mg\", \n",
    "                            surface_type=\"hcp0001\", \n",
    "                            size=size, \n",
    "                            a=a, \n",
    "                            c=c, \n",
    "                            orthogonal=True, \n",
    "                            vacuum=vac)\n",
    "Mg_0001.plot3d()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using selective dynamics\n",
    "\n",
    "We use selective dynamics to restrict relaxation to the surface atoms (first and last Mg layers). We use the advanced array indexing options available in the NumPy package (see [here](https://docs.scipy.org/doc/numpy-1.13.0/reference/arrays.indexing.html)) to detect which atoms are at the surface and then freeze the rest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initially freeze all the atoms\n",
    "Mg_0001.add_tag(selective_dynamics=[False, False, False])\n",
    "\n",
    "# Find which atoms are at the surface \n",
    "# (based on the z-coordinate)\n",
    "pos_z = Mg_0001.positions[:, 2]\n",
    "z_min, z_max = np.min(pos_z), np.max(pos_z)\n",
    "eps = 1e-4\n",
    "relax_indices = np.argwhere(((pos_z - eps) > z_min) \n",
    "                            & ((pos_z + eps) < z_max ))\n",
    "relax_indices  = relax_indices.flatten()\n",
    "\n",
    "# Now allow these atoms to relax\n",
    "\n",
    "Mg_0001.selective_dynamics[relax_indices] = [True, True, True]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup and execution\n",
    "\n",
    "To automate the calculation we define a function that has as input the project object, structure, job_name, Fermi smearing width, the type of k-point sampling and the plane-wave energy cutoff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_ham(proj, basis, name, sigma=0.1, mesh=\"GP\", encut=350):\n",
    "    ham = proj.create_job(pr.job_type.Vasp, name)\n",
    "    ham.set_convergence_precision(electronic_energy=1e-7, \n",
    "                                  ionic_energy=1e-2)\n",
    "    # Setting fermi-smearing\n",
    "    ham.set_occupancy_smearing(smearing=\"fermi\", width=sigma)\n",
    "    # Ionic_minimization\n",
    "    ham.calc_minimize(ionic_steps=100, \n",
    "                      electronic_steps=60, \n",
    "                      retain_electrostatic_potential=True, \n",
    "                      pressure=None)\n",
    "    ham.structure = basis\n",
    "    ham.set_encut(encut=encut)\n",
    "    if mesh == \"GP\":\n",
    "        # Only the Gamma point\n",
    "        ham.set_kpoints(scheme=\"GP\")\n",
    "    elif len(mesh) == 3:\n",
    "        ham.set_kpoints(mesh=mesh)\n",
    "    return ham"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "ham_vasp = get_ham(proj=pr, \n",
    "                   basis=Mg_0001, \n",
    "                   name=\"Mg_0001\", \n",
    "                   sigma=0.1, \n",
    "                   mesh=\"GP\", \n",
    "                   encut=350)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Submitting to the queue (optional)\n",
    "\n",
    "If you use a cluster installation of pyiron, you can send the created jobs to the cluster by specifying the name of the queue and the number of cores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# queue = ham_vasp.server.list_queues()[-1]\n",
    "# ham_vasp.server.queue = queue\n",
    "# ham_vasp.server.cores = 20\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Choosing an appropriate executable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['5.3',\n",
       " '5.3_col',\n",
       " '5.3_col_mpi',\n",
       " '5.3_mpi',\n",
       " '5.4',\n",
       " '5.4.4',\n",
       " '5.4.4_gam',\n",
       " '5.4.4_gam_mpi',\n",
       " '5.4.4_mpi',\n",
       " '5.4.4_ncl',\n",
       " '5.4.4_ncl_mpi',\n",
       " '5.4.4_std',\n",
       " '5.4.4_std_mpi',\n",
       " '5.4_gamma',\n",
       " '5.4_gamma_mpi',\n",
       " '5.4_mpi']"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ham_vasp.executable.available_versions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since this example uses the $\\Gamma$ point only, we can use the VASP Gamma-only version. If you use more k-points choose an appropriate executable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "ham_vasp.executable.version = \"5.4_gamma\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Execution\n",
    "\n",
    "The job is ready for execution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "ham_vasp.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Post processing\n",
    "\n",
    "To analyze the results we ensure that the job is finished (the `if` statement in the first line). We then compute the work function by subtracting the Fermi-level from the vacuum level\n",
    "\n",
    "$\\Phi = V_{vac} - \\epsilon_F$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "wf: 3.37343565133\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEPCAYAAABV6CMBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJztnXecXGd573/P9LIzs322a9VlyZIley13gwtgjLFpFxtiPjjANTihOMEXMCa5NwklBEgMcTAIAg4thmDHwQVsbIyNccHqlqxetu/O7s6WmZ0+894/zpzVStoy5Zzznt19vp+PPtJOOfNoz5nze5/6khACDMMwDFMoFtkGMAzDMAsLFg6GYRimKFg4GIZhmKJg4WAYhmGKgoWDYRiGKQoWDoZhGKYoWDgYhmGYomDhYBiGYYqChYNhGIYpCptsA/SgtrZWtLe3yzaDYRhmQbFjx45hIUTdfK9blMLR3t6O7du3yzaDYRhmQUFEnYW8jkNVDMMwTFGwcDAMwzBFwcLBMAzDFAULB8MwDFMULBwMwzBMUSwI4SCi64joEBEdJaLPybaHYRhmKWN64SAiK4B/A/BWAOsBvI+I1su1imEYZumyEPo4tgI4KoQ4DgBE9CCAmwC8LtUqhlkCpLM5RBIZTCYzSGZySGVySGVzSKazSGXzP+cfS2VyyAmBbA7ICQEhBHIC+ccERP7f6mO53Kl/qxABBJr2b+Vv5Wea93UEgs1KsFktsFkINgvBbrUoj1kssOefs1vyr7ES3HYrvA4bPE7lb5fdctpnMWezEISjGUD3tJ97AFx05ouI6HYAtwNAW1ubMZYxzAIkmxPoHY3jxMgkBsbjGJxIYmAigdBEAsPRFCKJNCKJDCKJDOLprGxzDYcI8Dps8LlsqPY6TvvTXOlGa7UHrVUetNV4UOFcCLdQ7VkI/+uZpF+c9YAQ2wBsA4COjo6znmeYpUg2J3CgfwK7usewq2sU+3sncGJkEqlM7rTXVXsdCPpdqK1woKnSBZ/TDr/bBp/LDp/LBq/TBqfNAqfNAofNAofVqvxts8BhtcBpV/62WggWIlgsUP4mgoUUb0F5TnmcCLDmn6f880IoXgmgfMGFEFNfdCEA9adpDsrU4+pjOSGQyQqkczlksmf8e9rfafW5bA6JdBaTqSxiqQxiqSxiyQwmU1mMx9MYnUxhZDKFzpEYRqJJTKZOF9LaCgfOafRjY3MAm1oCuGh5Daq8Dh3OpLlYCMLRA6B12s8tAPok2cIwpieWyuC5Q0P47YFBPHswhNFYGoByk9vUUok3rq3DijovltdWoDHgQr3fCafNKtlqRTxOjxCZL1w0HkujKxyb+nN8KIp9fRPY9vxxZHICRMCmlkq8YU0dbjyvEavqfbJN1gUSwtyLcyKyATgM4BoAvQBeBfB+IcT+2d7T0dEheFYVs9Q4OTyJH73Uif/a0Y1IIoOA246r19XjjWvrcH5bFVqq3By714lEOov9feN4/vAwnj8yhD3dY8gJoGNZFW6+sBVvP68JLrt8cZ4PItohhOiY93VmFw4AIKLrAdwLwArgB0KIL831ehYOZilxZDCCrz15CE+9PgibhXD9xkbcsrUVW9urYbOavnByUTIUSeLhnT34+avdOD48iTqfE3e8YSVuvXgZHDbznpNFJRzFwsLBLAViqQz++anD+MEfT8DrsOFDly/Hn13Uhnq/S7ZpTB4hBF46NoJv/e4IXj4exoo6L770jo24ZGWNbNNmpFDhWAg5DoZhzmBf7zg+/rOdODkSw/svasNdb16L6iWQlF1oEBEuXVWLS1fV4ncHB/F3j76O93//ZXzk8uX4zHXrYF+gHiELB8MsMB7Z1YvPPLQX1R4HHrz9Yly8wpyrV+Z0rl4XxMUravDlJw7ge384gb0947j/1gsWpOAvTLljmCXK954/jjt/vhtbWivx+CcvZ9FYYHgcNnzxHRtx782bsat7DO+5/0UMjCdkm1U0LBwMs0D4/h+O40tPHMDbNjXiRx/eipoKp2yTmBJ5x5Zm/OwjFyEUSeLmbS8tOPFg4WCYBcBDO3rwxccP4PqNDfjWLVtM0XfBlEdHezV+/OGtGI4k8ecPvIpoMiPbpIJh4WAYk7OraxR3P/waLltVg3tv3gKrhXsxFgtb2qrw7VsvwOHBCD7xs53I5RZGlSsLB8OYmPBkCh/7yQ4EA07c977zTd0DwJTGG9bU4f++fT2ePTSEbX84LtucguCrkGFMihACX3jkNYQnU/jOrRcsiRlIS5UPXLwM129swNeePIRdXaOyzZkXFg6GMSmP7u3HE68N4M5r12BDU0C2OYyOEBG+8q5NqPc58bmHXkM6m5v/TRJh4WAYExJJpPH3j76O81oC+OiVK2SbwxhAwG3H3924AYcGI/ieyUNWLBwMY0Lue/YohqNJ/P1N5/K8qSXEmzc04C0bgvjWM0cwOGHeEl2+IhnGZHSNxPCDF07gPRe04LzWStnmMAZzz/Xrkc0JfPOZI7JNmRUWDoYxGfc9ewREhP/zlrWyTWEk0Fbjwfu3timTdYeiss2ZERYOhjER3eEYHt7Zi/dvbUOQp9wuWT5xzWo4bRbc97ujsk2ZERYOhjER3/79UViI8LE3rJRtCiOR2gonbr6wFb/a04e+sbhsc86ChYNhTEJ4MoWHdvbiPR0taAiwt7HU+fDlyyEA/PCPJ2SbchYsHAxjEn7+ajdSmRxuu7RdtimMCWip8uBtGxvxs1e6EEmkZZtzGiwcDGMCsjmBn7zciUtW1GBN0CfbHMYk/Pll7ZhMZfHonn7ZppwGCwfDmIBnD4bQOxbHBy9dJtsUxkRsbq3E2qAPP3+1S7Ypp8HCwTAm4KGdPaitcOCac4KyTWFMBBHhlq2t2NMzjtf7JmSbMwULB8NIZjyWxjMHQnj7eU0Ldg9qRj/euaUZDpsFv9jeLduUKfgqZRjJPP5aP1LZHN65pVm2KYwJqfQ4cM26ejz+Wj+yJtmvg4WDYSTzyK5erKzzYmMzT8BlZuZtmxoxFEni1ZNh2aYAYOFgGKmEIgm82hnGjec1g4h39mNm5up19XDZLXh8rzmqq1g4GEYizxwIQQjgLedyUpyZHY/DhmvWBfHrfeYIV7FwMIxEnto/gLZqD9Zy7wYzD9dvbMRwNIXtJghXsXAwjCSiyQz+eHQEb14f5DAVMy9XrqmFzUJ49tCQbFNYOBhGFs8dGkIqm8ObNzTINoVZAPhcdlzYXo3fHwrJNoWFg2Fk8ftDIQTcdlywrEq2KcwC4ap1dTg4EJE+MdfUwkFEXyOig0S0l4j+m4h4OzRmUSCEwAtHh3HZqhpYLRymYgrjqrX1AIDfSw5XmVo4APwWwLlCiE0ADgO4W7I9DKMJx4Ym0T+ewGWramWbwiwgVtVXoLnSjWclh6tMLRxCiKeEEJn8jy8DaJFpD8NoxQtHlBXjFavqJFvCLCSICFeuqcXLx0ekluWaWjjO4EMAfi3bCIbRgheODqOt2oO2Go9sU5gFxsUrahBJZHCgX97QQ+nCQURPE9G+Gf7cNO019wDIAPjpHMe5nYi2E9H2oSH55WoMMxuZbA4vHw/j8tUcpmKK56LlNQCAl4+PSLPBJu2T8wghrp3reSL6IIAbAFwjhJjVNxNCbAOwDQA6Ojrkt1YyzCwc6I8gmszg4hU1sk1hFiANAReW13rx8vERfOSKFVJskO5xzAURXQfgswBuFELEZNvDMFqgDqq7sJ3LcJnSuHhFNV45EZaW5zC1cAC4D4APwG+JaDcRfUe2QQxTLjs6R9Fc6UZjwC3bFGaBouY5ZG3uJD1UNRdCiFWybWAYLRFC4NWTYVyyksNUTOlc2F4NANjRGcbGFuPH8Zvd42CYRUV3OI5QJImO/BefYUqhMeBC0O/E7u4xKZ/PwsEwBrK9k/MbTPkQEba0VmEXCwfDLH62d47C57JhTT2PUWfKY3NbJTpHYhiJJg3/bBYOhjGQ13rGsaklAAvPp2LKZEurMrpvT4/xXgcLB8MYRDKTxcGBCZzLe4szGrCxJQCrhbCri4WDMYBUJodQJCHbjCXH4YEo0lmBTc085JkpH4/DhrVBn5QEOQvHEiORzuK2H/4JW7/0DN773ZdwcEDevJulxt5e5Qu+SUL5JLM4Oa+1Enu6xzDHUA1dYOFYQuRyAnc+uBsvHhvB+7a24Wgoirv+aw9yEqdsLiVe6xlHpceOlipu/GO0YUOTHxOJDHoN3tiJhWMJ8cdjw/jN/gF87q3r8JV3bcTf3HAO9vVO4Fd7+mSbtiTY2zOOjc0B3l+c0YwNTX4AwH6DO8hZOJYQD/6pG1UeO/78snYAwE3nNePcZj++9uQhZLI5ucYtchLpLA4PRrCRE+OMhqxr8MNCMHz0CAvHEmEkmsRTrw/gXee3wGmzAgAsFsJfvnEVesfiePXkqGQLFzeHByPI5ARXVDGa4nZYsaKugj0ORh8e2tmDdFbglgtbT3v8yjV1cNgsePrAoCTLlgYHByIAgHUN3PjHaMv6Rr/hmzqxcBTJcDSJm+57Ae+5/0Xs7Fo4q/Qn9w9iU0sAq4On37i8ThsuW1mD374+aHhlRqns7RnDDf/6B/zlT3fi+cMLY9OuQwMRuOwWLKvxyjaFWWSsb/KjdyyOsVjKsM9k4SiC4WgSN3/3JRwajKAzHMO7739xam8FMxNLZbCnewyXrZp5x7lr1wfRFY7hSChqsGXFMxZL4Y6f7ET/WAJ/OhnGx36yA4MT5u9JOTQQwep6H6zcMc5ojJogNzLPwcJRBN997hg6R2L4jz/fit99+g2o8TrwrWeOyDZrXnZ0jiKTE7PuOHftOUEAwG9fN3+46p5H9iEUSeD7H+zAQx+7FJmswFd/c1C2WfNycCCCtRymYnRgfWNeOAwMV7FwFEgyk8Uvd/TgTeuDuGhFDXwuOz58+Qr84cgw9kiaUFkoLx8fgdVC6Fg280TWoN+FdQ0+vHLC3N5TaCKBX7/Wjw9fvgJb2qrQVuPBhy5fjod39mJf77hs82ZlOJrEcDTJ+Q1GF2oqnKjxOnDUwIgBC0eBPLl/EKOxNN5/UdvUY7de3Aa/y4bvPn9MomXz88rxMDa1BOB1zr5v12ZJHajF8NjefuQE8J4Lmqceu+ONK2G3Eh41cS/KoanEuF+yJcxiZVV9haGhZhaOAvnPV7rQWu3GZStP5Ql8LjvesaUZzx4cQjKTlWjd7MRSGezpGZs1TKVyXmslxuNpdI6Yd2v3/9nThw1NfqyaNpI84Lbj4hU1pq4KUyuqOFTF6MXqYAWODEYMW/ixcBTAeCyNl0+M4J1bWs4ah33l6jrE01nsMGkfxJ7ucaSzAluXz73j3Hkt8kY0F8KJ4Uns6R7DTZubznru2nOCODY0ieND5kzuHxqYQI3XgTqfU7YpzCJlTdCHiUQGoYgxe3OwcBTAqyfDEAK4dIZ9oi9ZWQO7lfDcEXOWhaoJs3Ob5m48WxOsgMtukbYV5Xw8nU/cv23T2cJxzTn1AIBnDoQMtalQjoaiWB2skG0Gs4hZVa9cX0cGjVk8kZlj2qXi8/nEBRdccNpj733ve/EXf/EXiMViuP766896z2233YbbbrsNw8PDeM973nPac13hGOKrrsaJn38JQwN9+MAHPnDa86/3T6D1yvdix/c+h0OHDuGjH/3oWcf/whe+gGuvvRa7d+/GnXfeedbzX/7yl3HppZfixRdfxOc///mznr/33nuxefNmPP300/jiF7941vPf/e53sXbtWjz66KP4xje+MfX4saEoxmJpvP78Y2htbcXPf/5z3H///We9/5e//CXu+OVhHP/j4/B1//Gs55944gl4PB58+9vfxi9+8Yuznv/9738PAPj617+Oxx577LTn3G43fv3rXwMA/uEf/gHPPPPMac/X1NTgoYceAgDcfffdeOmll057vqWlBd633IkD/RGc3/c/2L1792nPr1mzBl3rb4XfbUfljh/i8OHDpz2/efNm3HvvvQCAW2+9FT09Pac9f8kll+ArX/kKAODd7343RkZGTnv+mmuuwd/8zd8AAN761rciHj99oNwNN9yAu+66CwDwxje+8azfzRHvRtz2kdtxz1tWFn3tAcAdd9yBm2++Gd3d3WddewDw6U9/Gm9/+9tNd+2p/PjHP5732qutrcUDDzyABx544KznZV97P/nJTwAAd95554zX3rZt2wAAt99+u7RrbyiSROuGDrTXeNEQcE09X+x977nnntshhOg464VnwB5HAUzE02iv8cJlt874fKXbju7RmCn3uIilsvA4Zrb7TM5rqUTvWBxmXEvs6BybtSoMUDrgd3eNIWuySb/pbA6JdBYr6tjjYPSjtsIBm9WCeNqgXKsQYtH9ueCCC4RWRBJpseLux8XXnzw462te6xkTyz77mHhkV49mn6sF6UxWrL7nCfGlx18v6PX/s7tXLPvsY2J/77jOlhXHyeGoWPbZx8SPXzo562se3aPY/lrPmIGWzc8rx0fEss8+Jn53cFC2Kcwi5z33/1H8r/tfLOsYALaLAu6x7HHMw/aTYWRzAhctn70qaW2DDw6bxfBBY/NxYngSqUwO5zQWVs2zJh+HPxKK6GlW0ezoVAoPOtpn9zjMmtxXE/ar2ONgdGZVvQ+HDfrusnDMw47OUVgthPOXzb7dp91qwTkNPtM1oamJ8UL7B5bXemEh4JjJRo9s7xyFz2nD6vrZBbClyo0qjx17u811Do4PT8Jhs6CpkjdvYvRlZZ0XY7E0Rif1n1nFwjEPB/onsKLWC49j9uY5ANjQHMC+3nFTNdAdHIjAbiWsLHC167RZ0VbtwVGTlbXu7BzFlmVVc855IiJsbKnEXpOJ97FQFMtrvDyjitGd9vwAzRMjk7p/FgvHPBQ6Y0jdwrFn1NgtHOfiQP8EVtZVwGEr/DSvqq8wdHTBfKgbIG0uYJ/uTc0BHB6MIJ4yTzPm8eFJrKznibiM/rTXKtdZJwuHXKJJRQgKmTGk9kns7zPPivdoKIo1weK6lVfWV+DE8KRpdgQ8GooiJ4C1BYTbNrUEkM0JQ4e9zUUqk0NXOIYVtZzfYPSnrdoDCwEnhvWf/sDCMQfFzBha26CMzN7Xa46bVjKTRd9YfGoVUiir6iqQzgp0hc0xeuTwoHIO1hTQQLcpnyA3S66pKzyJbE6wx8EYgsNmQXOVGyeH2eOQysEBRQQKCVW57Fasrq/APpN4HN3hOHICWF7rKep9ageqWcJVhwejsFupIAEM+p3wuWymsV1d+bXz5k2MQbTXeHGSQ1UKRHQXEQkimnknIp04NBBBhdOGlqrCKmI2NAUM3zR+NtRVR7E7zq1UhcMkCfIjgxGsqK2A3Tr/pUqkFAIcM4ntqtfGu/4xRrG81osTw5O6F+mYXjiIqBXAmwB0Gf3ZamKcqLCKmJX1XoQiSUQSaZ0tmx911bG8yJuW32VHg99lmlX7ocFIUXOeTCUcI5PwOW2o8thlm8IsEZbVeBFJZBDWuSTX9MIB4F8AfAaAoXWuQggcKnLXNjUJetKA5NR8nBieRMBtR5XXUfR7V9QpqxbZTOaLE4pJ8K+s92Jwwhzi3RWOobXaU/DCg2HKRQ1N6x2uMrVwENGNAHqFEHsKeO3tRLSdiLYPDZU/qXYomsR4PI019YWvdlfUKav748PyV7ydI7GiE+MqrVUedIfllxWrXk8xwqGK9/Eh+cLXGY5hWU1xOSaGKQc1n6b34lW6cBDR00S0b4Y/NwG4B8DfFnIcIcQ2IUSHEKKjrq6ubLvUDY2Kufkuq/GAyBw3rRPDk2gv8abVVuPBcDQpvR/iUBEVVSqr6s0h3tmcQE84jjYWDsZAWqs9sFpId49j7nZoAxBCXDvT40S0EcByAHvyrn4LgJ1EtFUIMaC3XapwFJPYdNqsaKly47jkME8inUXfeBztNS0lvV8tBugejRXdB6Ilx0JROKyWos5BW7XSpX0sJPccDE4kkMrm0FbNwsEYh91qQWPAhW6dy+mlexyzIYR4TQhRL4RoF0K0A+gBcL4RogEo3ZcWApqLnDG0orYCJySvdrvDMQihVFiUQmv+Zqf3xTcfXeEYWqrdRY3rcNgsWFbtkZ4gn1p4VHNFFWMsLVVudOs8wcK0wiGbzpEYmqvcRY3rAPLlcEP6l8PNxckpb6m01W5rlTmEo3MkhmUlrNhX1HmlC4f6u2OPgzGa1ioPekZN6HEQ0cen/Xv2sbEakvc8ho34LCCf2CxhtbiizovJVNawvX9nQr1oSr1p1VY44LZb0SUxQS6E0r1eyv9hZV0FTg7HpG7q1BmehM1CaKp0zf9ihtGQlioPBieSSOi4qVOpHseyaf++WwtDzEbXyGRJiU21qkfmirdvLA6X3YLqEkpxAaWRrrXajW6dVy1zMRpLI5rMoK2E5rnWag9S2ZzUHRm7wnE0V7lhK6BxkWG0pLVaCa/3jem38Cv1qrYQ0RVEZAEw+w5HC5TxeBqjsXRJVUnL1ZJciZVVvWNxNFW6y+ofUEpy5QmHOuGzFI9DTe7LnFTcFY5NhfwYxkimcpQ6Xv+lCsdnAJwH4HsAfqWdOeaga0QN9RS/2m30u2C3ktTVeu9Youik/pm0VnvQMxqXlqs5Na6jdOHolSgc/WNxDlMxUji1cNLvHlRSOa4QIgvgPo1tMQ2dYXXOU/E3LYuF0FTplnrT6huLY93a+rKO0VrtQTSZwWgsXXLIqxxU8S5l1d5cqbxH7wThbKSzOQxFk2gI8K5/jPEEffnFq445ypKEg4j+H4CLAPQC2CWE+DctjZJN50h5yeXmSjd6dYwvzkUyk8VQJInmAgczzkar2ssRjskRjnAM9T4n3A5r0e91O6yorXBIOwdDkSSEABoD7HEwxmOxEJor3bounEoNVVUCeBnAlwCs1c4cc9A1EkNthRNeZ2n9kc0SPY7+MSUhXO4e1y35lb6sm2+54zqUL44c2wcmlHPQ4GfhYOTQWu0xZY4jDMAKIJT/96Kibzxe1oq9pcqDUETfcrjZUCspyo2vq6vlgXE5lUnd+QGBpdJS5ZEnHPnfWQN7HIwkWqrc6NGxuKVU4fg6gO8A+BYAc+xcpCF9Y3E0lfGlV0WnX8JNV/UQWirLq+ip9NjhtFmmVs9GkkhnMTCRKKvruqVKCRfmJPRyTAkHexyMJFqqPBiZTCGWyuhy/FKF418B3ADgfwsh/kVDe6QjhED/eAKNZSQ21YomGeGq3rE4iIBgwFnWcYgIDQGXFPEbGE9ACJTl9TVXuZHK5DAcNb4Rc2AiAafNgkreh4ORhHoP6hvT5/tbknAIIT4M4ACA7+VHny8aJuIZxFLZskI9U+WgY8ZX9fSNxVFX4YTTVnxS+Uwa/C4MShCOvvF8uK0Mr+/UoEbjxbt/PIGGgIv34WCkoYaa+8f1uf5LHTlyJYCNAOIAvqipRZJRb1rleBwNARcsJKcBrW8sUXZFlUpjwIX+CeP/D2qop7GMBL/M5P7geILDVIxU1OIYvSIGpYaqNgDYC+AeIcQmDe2RjqrQjWV4HHarBQ1+l5RQVV++a1wLGgJuDI4nDc8T9GuQI1BddRm9HP0TcU6MM1Kp9yuh6n5ZoSoiuo+ILp3+mBDifiHEC0KIRZcYV29a5dbgN1e50WPwancqP6PRarfB70Qqm0M4pu/+xWfSPx5HlcdeUg+HitdpQ6XHruu8npkQQmBwPMnCwUjFabOitsIpNVR1BMA3iOgkEX2ViDbrYolJ6B9LwGoh1PvKFA4JvRyRZAbxdBZBrYQjH64zuiS3fyyhSdd1g9+FgXFjk+PhyRRS2RyHqhjpNOpY3DKvcAghvimEuATAG6D0bPyQiA4Q0d8S0RpdrJJI33gcQZ+zqM2DZqK5yo2BiQQy2ZxGls1PaEK5SapuarnI6uXoG0+UlRhXCfpdGDS4nFgtX+aucUY2inBITo4LITqFEF8VQmwB8H4A74RSWbWo6B9LlJWUVWkIuJHNCYxMGhfmCeVvWuV6SypquKXf6JvveLysHJNKg99leB+KKrJaeX0MUypNlW75yXEishPR24nopwB+DeAwgHfrYpVE+sfjmqwW1VCFkat1dfMorTyO2grF8xrQadUyE/FUFqOxdFlVbSrBgAvD0STSBnp9pzwOHnDIyKUh4EIkkUE0qX0TYCHJ8TcR0Q+g7Pl9O4AnAKwUQtwshHhEc4skoiaXtahKCuZv3kaueNWwjFarXauFEPQ5Dc0TaBnqafC7IIQydNAoBsYTsJCyiyLDyGSql0OHApFCPI7PA3gJwDlCiLcLIX4qhJC3S5GOhCdTSGZymnocIQOFIxRJwuOwoqLE4Ywz0RBwYcDAXg71IteiKqkhYLx4D4wnUO9z8c5/jHRUr1ePcNW8dxghxFWaf6pJOVWKW77HUaOGeQz2OLSOrTcEXDg4ENH0mHOhnoMmLUJV+d+Fkd3vAxMJBDkxzpgAPbvHeVk0DXWPai1Wu0pJr7FhnlAkiTqfNvkNlXqfC0MTxv0f1ItcE49DzTMZ7HFo1UfDMOUQ9LtApM+8qlJHjjRobYgZGMzfIIMaJZfr/a4pMTKCkA4eR53PqfSHpIwZEd8/nkC11wGXvfxZW1UeB+xWY72+gfycKoaRjcNm0a0JsFSP4wlNrTAJA+MJECnVRFrQ4HcaVlUlhEAokkS95h6HcjyjBFDLcJsl38hpVKgqmswgksywcDCmoTHgmloQa0mpwrEox36GIgnUeJ2wa5TYNLKPIJpUpvpq5S2p1KtJfoMqk7QWPyW5b8w5GNBoXA3DaEW9T58m2FLvkN/T1AqTMDiRnKrE0YJ6v1JHrddmKtOZ6uHQqPlPZcrjMCjPoXgcGgqHX58V10xoXQ7NMOUS9Dt1WfSVuh/Ht7U2xAwMjCcQ1PDGqyZnjbhxqTctrZr/VIwMVWVzAsPRlKbiF/S78htD6T/hV6sBmQyjFUG/K99moG2OkquqphGKJKZCM1rQYOCsJ7XJTevVbpXHAZuFDAlVjUwmkc0JbT2OgBPxdBYTCf29PvY4GLOhfpe0boJl4ciTzuYwHE1petNSj2XEan3K49A4OW6xEOp8TkNCVepn1GnscQAwZNihOg5ei4owhtGCep2iHiXtx7EY0WPFHjRwXlVoQvuucZVnJ703AAAfO0lEQVR6n9MQ8VM/Q+scB2DMORgY174cmmHKQQ29az3BwvT7cRDRJ4joEBHtJ6J/0utz1MobLfdR8Lns8DqshlT1DOarkfTY57rO5zJk3tOpsfA6hAsNOAcDEwnObzCmQl2Eae1xm3o/DiK6CsBNADYJITYA+LpenxXSKbls1J4QoQlt8zPTqdepMuNMVHe6TqM+GsDYsSPc/MeYjWqv2gQrKcchaT+OOwD8oxAimbchpNcHneoa1/aLHzSoHFSP5j+Vep9T2dkuo+948lBE6Rp32LRLvbnsVlR67Lp7HKmMkiNr8PM4dcY8EClNsDJCVaoBMvbjWAPgCiJ6hYieI6IL9fqggYkE7FZCtUfbcdgNAZdBOQ794utqeexwVF8BHJzQR/waDPD61ONr2QfEMFoQ9DsxqHGOct5MKhG9CcD7ALwNwJ8APAjgdq1GqxPR0wBmmn11T96+KgAXA7gQwC+IaIWYoSifiG6Hsl8I2trairZjcEIZh20pc8vYM1HCPAnkckLzY6tEkxlMprK6ehyA4tVosVfJbAxpXA6tEjSgg/+UcLDHwZiLoN+FI6GopscspATn8wB+BuAuIURY008HIIS4drbniOgOAA/nheJPRJQDUAtgaIbjbAOwDQA6OjqK7vYanEhoPlkWUFa76azAaCyFGg1j99MJ6dw/oOZ99N5bZHAiidVBn+bHbfC7sL9vQvPjTkdt/tOyuIJhtCDod+GFo8OaHrOQ5PhVQojv6SEaBfAIgKsBIJ+IdwDQ9jeQZ2Bcn4oYI0Z7qzkUvTwOVZD0TJDncgJD0aTms7YAZQvZkUl9t5A95XGwcDDmot7v1Hz0kdkbAH8AYAUR7YMSIvvgTGEqLVDmVGn/pa83oAFN7X/Qq6qqxusAkb7CMTKZQjYnNJ+1BRizhWz/eAIehxV+l/Z9NAxTDqd6ObS7/k19lQshUgBu1ftzIok0osmMLmEGVYz0rKw61f+gj8dhs1pQ43ViSMcmwCnx0yNcOG0LWb1yNAMTCTT4Xbr00TBMOUyfntBe69XkmGb3OAxBrXrSxePI3wj1rKwKRRJw263w6dA1rlKv89iRqTlPOpwDI3o5uIeDMSvqwmlQQ4+bhQOn8g9a7DV+JnarBbUVDl1DVYMTSdT79ekaV9G7CVD1yHTx+gzIMw2MJzgxzpiSqT11NLz+WTigf0WM3t3joYi24+BnQu95Verui3pUtlV7HXBYLboJRy4nMDjBHgdjTnxOG9x2q6b3IBYOnAoj6ZUjUHYC1DfHUaeT7Sr1PheGo0oCWw+03n1xOkSEer9Tt1DVyGQKmZxg4WBMCREh6Hdqeg9i4YASwqjxOnQbh12vu8eR1N/j8DuRzQmEJ1O6HF+ZLKuf+Om5je8A93AwJkfrexALB/Qfh92g0y5cADCZzCCazOjmLanovRPg4ERS1xtvMKDfzLAB7uFgTE7Qr+28KhYO6Nf8p6JWNehRlRSa2kdEX+Go06EWfDqDOk73BfIeh05byA6Mx5XPYOFgTErQ58TgRFKz65+FA/kafB2/9Ho2AZ7a+U//5Digj8eRyuQwMpnS1eNo8Lt020J2YCIBm4VQ6+UBh4w5Ceav/0hSm+t/yQtHIp1F2ICbFqBPE6BxHoeeXpP2O/+dSb1OG9oASlVe0K/9gEyG0Qq1P0qrcNWSFw71Rqinx6FnH4F6IWi5T/dMuOxWBNx2XXo5pvZCMeIc6FBZpXdin2HKJehTF07afH+XvHD05+PTejT/qVR67HDYLLqsdkORJFx2iyEzkvTq5ZjqGtdR/PTcQlbZMpbHqTPmJahxuHzJC8epihj9VoxqHbVeOY56nzEzkvTqHjdisqxeY0eEELpX5TFMuZwK1bLHoQmn5lTpu2JUq3q0JjShzyjymVC2oNReOAYmEnBYLajy2DU/topeW8hGkhnEUlldq/IYplw8Dht8Lht7HFrRP55AhdOGCh0HBAL6jR0ZjCR0r6hSqfc5MRTRrqRPJWTArC1Any1k1cWAnvkZhtECLe9BS144jJoxpJw07W+6Q/mbrhEE/S6ksjmMxtKaHteoAYF6bCGrCgd7HIzZUcaOsHBoQr/OzX8qevQRxFIZRJIZwzyOqQSzxiG3wYgxOQIlXKhtqI3HjTALhcaAW7Pv7pIXDqMSm8GA9k2Aar7BqByH1pUZKoMGngOtt5BVV3CcHGfMTmNACVVlNLj+l7RwZLI5DEWThngcp+qotbvpGtU1rqJHSWs0mcFkKmuI+KlbyGpZGdY/nkBthQMO25L+KjELgMaAGzkBDEXLv/6X9NWujgk3IsehR5jHqK5xlXqfE0Ta/h/03H3xTKa2kNXQ/sEJLsVlFgbqArlvrPzrf0kLh9r8Z1RiFljYHoc9v/e4tuE24/4P6n7j6nnXAqNyZAxTLo2V2i1el7RwGNF4pqKO7NByXtVQJAmnzQK/W/+ucZWGgHaVGYCxI8nV7u6+Me2Egz0OZqGgXv9aLJyWtHD0T5VSGjMuQuvNhPrHlVJiI7rGVbRuZBw0MMHvdyn9Olq46oAxAzIZRiv8Lhs8DiuHqsplYCIBh03fjuXpBAPaNqAZ1f8wHa0bGQcnEvC5bPA49PeaiAhNlS7NPA514dFcxXOqGPNDRGgMuDAwwR5HWag3XqNW7EGfU9PVev9E3PD4eoPfhdFYGom0NrsZGh3qaQy40adRjkMVIDV3wjBmpzHgZo+jXNRQj1E0BFwYjiY1qaPO5QQGx5O6z9g6E/Umr9XMqoEJY72mpko3+jUKVfXmhaOZhYNZIDQGtAk1L2nhMDrUE/S7kBNKGXC5hGMppLI5wz2OoMa9HCEDR6YAQHOlCyOTKU08pr6xOIi4+Y9ZODRWuhGKlN8EuGSFI5cT6B+PGxqfVm8wWtx0jex/mI6Wm1LlckKZFWawxwGcyk+UQ+9oHPU+Jzf/MQuGxoCyeB0sswl2yV7xoUgS6awwNMzQoGEvhxpfl5HjALTZ1yIcSyGTE4bnOABtSnL7xuOc32AWFI1TjcjlXf9LVjh6x2IAjK2ICQa0GzuirviN3nnO77bBbbdqsmKfGkluoHCoC4VeLYRjLMH5DWZB0TR1/Zf3/V2ywtEzqtw4Wgz84td4nbBaSJPkVP94AnYrocbr0MCywlFLWrVoIpraMtbAHEcwoIxNKTdBLoRA71ichYNZUEwtnEYXscdBRJuJ6GUi2k1E24loq1bHVoXDSI/DaiHU+5yadI+rU30tFuOa/1SaKt0ahXryfRAG3nydNitqK5xl2z8ymUIqk+NQFbOg8DptqPY60DMaK+s4phYOAP8E4O+EEJsB/G3+Z03oHYujymM3pPFsOlo10PWPG9/DodIUcE/d9MuhdzQOu5VQW2GcxwHkha9Mj0ldsbFwMAuNlir31MK5VMwuHAKAP//vAIA+rQ7cOxpHS5VHq8MVjFa7cA2MJwzv4VBpqnRjKJJEMlNeSWvfWByNAbfhXlNLZflfnFPNf1yKyywsFOFY3B7HnQC+RkTdAL4O4O7ZXkhEt+fDWduHhobmPbCs+LQW+14LIaROZVWnbA6WuZten6Rz0FLtRu9oHLlc6dv4cvMfs1Bpzi+cytnGWrpwENHTRLRvhj83AbgDwF8JIVoB/BWAf5/tOEKIbUKIDiFER11d3ZyfKYRA76ixPRwqwYALkUQGsVTpW8iOxtJIZnLShutpVZnUNyannLW1yoNUNofBSOkC3jeWgMehTDxmmIVES5UHyUyurEZkYwP8MyCEuHa254joRwA+lf/xvwB8X4vPDE+mEE9npawWg75TM/FX1FWUdAzZYRLV0ymnsiqTzWFgIoFmCf+HlvyCoTscL7mcuXs0htYqj6GTiRlGC9Trv2c0hjpfaflF6R7HPPQBeEP+31cDOKLFQafCDBI8jkYNdgKcqgirND5HA5xKCJdTmTQwkUBOyEkut1Yrv7dy4rxdI7Gp4zDMQkLN7ZaT55PucczD/wbwTSKyAUgAuF2Lg/aOyotPq2LVU8ZNV6bwAcqmVNVeR1mVVeqEThnCoZ737nBp50AIga5wDJetqtXSLIYxhKl70GIVDiHECwAu0Pq46i+sVUJVVWPADaLyTlrPaAweh9WwfURmotx9LWSOJHfZrQj6negu0eMYjiqhzrZqTowzC48Kpw1VHntZHrfZQ1W6cHJkEgG3HQEJN16HzYIGv6usk9Y7qlQjyYyvNwbKG08uuyqptcqD7nBp56Ar/762Gg5VMQuTlipPWYvXJSkcXeEYlkn80pfbgNM7JqcibDrNZXaP943FUe11wO2wamhV4bRWl/7FUQWnrdqrpUkMYxgtVe6SPW5giQpH50gMbRITmy1VnrJmxfSMxqcqI2TRVOlCJJnBRCJd0vt7x+JSm+daq9zoH48jXcK+BKrHIfscMEyptFV70BOOI1tiL9OSE450NofesTjaa+StFlvKuGlFkxmMx9PSKqpU1MqMUsM9ate4LFqqPciJ0irDusIxNPhdcNnleEsMUy7ttV6ksrmSS+qXnHD0jSkqKzM+3VLlRk6UVpLbK2E440yoHlsplUlCCHSH41KKE1Raq0q3vyss12NlmHJRF84nh0tb+C054Tg5ovyilkkOVQEoKcaoJtVlh0lay/A41KokmXkmdeFwcmSy6Pd2h7mHg1nYLK9VhONECdc/sASFoyv/i1omOVQFlFaSq1YjGbmPyEwEPHb4XbapeH8xdIWVcyBz1d7od8Flt+DkcHFfnEQ6i4GJBHsczIIm6HeWdP2rLDnh6ByJwWmzoL7EVnstKKeXo3c0DofVYvgo8ploq/GU5DWpYiNz1W6xENprvDhe5BenZzQGIYC2Gk6MMwsXIuX6Z+EokM58Ka6MDZBUyunl6BlVqpFk2q/SWuUpyePoHImBSH64bUWdFyeK/OIcDSmvX1ninDGGMQvtNV4OVRVK10jMFPX3LVVu9JSQmO0MT6JNYphtOm35Xohix5ObpSppea0XXeFYUdVtx4aiAFDygEqGMQvttV50h2MlleQuKeHI5QQ6w5NSk7IqbdVedIaLU3shBE4Ox7DcBPYDSklrKpNDKFLcvhxmSS6vqK1ANieKSvAfG4qiwe9ChdPU03oYZl6W13qQzoqSStKXlHD0jsWRSOewql7+anF5rQeDE0lMJgvfl2NkMoVoMoP2WvN4HACKDld1jsSkVrWpLK/LV5YUEa46Foqa4vphmHJRS3KLDdcCS0w4joQiAIA1Qflf/OW1ig3FlIOqiSyZzYvTOdXLUbhwxFNZhCJJU1QlrcgL8PGhws6BEALHhiaxss4cv3+GKYepklwWjrk5PKjEp1fV+SRbArTX5vsIimjAUU+wWTyOpkoXiIrzONSCADMMCKz0OFDtdRRcWRWKJBFNZrCSPQ5mEVDnc8LvsuHwYKTo9y4p4TgyGEW9zyllKu6ZnHITowW/p3MkBquFpFcjqThtVjT6XUUJR+eIOiBQvnAAyqqr0HNwLKS8jiuqmMUAEWFdgx+HBlg45uRIKII1QfneBgB4nTYE/U6cKMbjGJlES5Ubdqt5TtvK+oqpSqNCUENzZhGOFbXeqRLb+Tia/39yjoNZLKxt8OHQYARCFFdZZZ47kM7kcgJHTZbYLGa1CwCdI5OmyW+orKyrwLFQtOAL79hQFFUeO6q9Dp0tK4x1jX4MR5MYKqAy7FgoigqnTWrzKMNoydoGHyKJTNG7eS4Z4egbjyOWyprG4wAU4VBnZ82HWorbboLcwHRW1ldgMqWM4SiEo6EoVtf7pG5CNZ31jX4AwIH+iXlfe2gwgpX1FaaxnWHKZV2Dcj88NDD/9T+dJSMcR/KJ8dUmqKhSWV7rRXgyhfHY/HtaDEfNVYqrolYYHQ3N7zkJIXAkFDVVclkVjtfnEY5cTmB/7wTObfIbYRbDGMKavHAcLDLPsXSEI1+Ku9pEN62pBHkBJblqbsBswqGG/goRjpHJFMZiaVOdg4DHjuZKN/b3zS0cneEYIskMNrUEDLKMYfTH71Ku/2IT5EtGOPb2jKMp4EKlxxyxdUCZlQScqtaZC/XEminUBgB1FUpJXyEJcjN6fQCwvsmP1/vG53zN3p4xAMC5zSwczOJiTbCChWM2dnWNYcuyKtlmnEZ7jRcOmwUHC4gvHuifgM9lQ1NA3narM0FEWFlfUZDHcTTv9ZmpQAFQwlXHhycRS83exb+vdxwOm8V0ws0w5bKu0Y9jQ1Ek0tmC37MkhCM0kUDvWBxbWitlm3IaNqsFa4M+HOifX+0PDkRwToPflInZVXUVBZW0Hs1XJTX4zSV+65v8EAJzrrr29ozjnEa/qUqhGUYLzm+rQjorsKd7rOD3LIlvwa78L2RLm7k8DgA4p9GH1/sn5ixnzeUEDg1EsK7RnKvdlfUVGI4m503yH8mXQ5tN/NQE+b5Z8hy5nMD+vglsbObEOLP46MhHYl49GS74PUtCOHZ2jcJhteBcE37x1zf6EZ5MzTlhtncsjmgyg3UN5rMfODX7a76Qm9n6aFRaqtwI+p14+fjIjM+fGJlENJnBpmZzeawMowVVXgfWBCvwp5OjBb9nSQjHrq4xrG/yw2mTu//DTJxTQDmo2mNwjkk9jo35G+rentkTzIMTCYQiyam6cTNBRLhsVS1ePDo8494iLx1TBOX8ZSwczOLkwvZq7Oxk4Zginc3htZ5xbGkz55d+nSocc5SDHhyIgMh8FVUqdT4nmivd2N0ze4x0V5fy3PkmK1BQuWJ1LUZj6RkF/NmDIbRWu3lGFbNo2bq8GtEitnhY9MLx4rERxNNZXLyiRrYpMxJw29FS5Z6zc/lA/wSWVXvgNfHmQZtbK+dMru3uHoPdSlP5BLNx2apaAMAfjgyf9ngincUfjw3j6rX1psvNMIxWdLRXF/X6RS8cv9rdB5/ThjesqZNtyqyc0+if1eMQQmBvzzjWm7xjeXNrJXpG4xiOzpyr2dU1ivVNAenbxc5Gvc+FtUEfXjg6dNrjLx0bQSKdw9XnBCVZxjD601zpxoYi7jGmEA4i+l9EtJ+IckTUccZzdxPRUSI6RERvKea4iXQWT+4fwHXnNpj2hgUAW9oqcXx4EqHI2fOeusIx9I7FTesxqZzXquY5zvY6Mtkc9vaMm64c+kyuWF2LV0+OYmSa+P3uYAgehxUXLS9uRcYwC43HP3lFwa81hXAA2AfgXQCen/4gEa0HcAuADQCuA/BtIipYAZ49GEI0mcGNm5u0tFVzrlileEN/PDp81nMv5hOzl66sNdSmYjm32Q8LAbu7z06QHxqMIJ7OmjbPpHLL1lakszls+8NxAEB4MoVH9/bhytV1pl54MIzRmEI4hBAHhBCHZnjqJgAPCiGSQogTAI4C2Drf8QYmEghPpvDNZ46gtsKJS0y+Wt/Q5EeVx35WfB1QhKPe5zT9dqUehw1rgj7s6Dy7FnwqMW7CPprprKr34cbzmvCjFzsxHE3in35zEJFEBn/95jWyTWMYU2EK4ZiDZgDd037uyT82J0ORJC7/6u9wbCiKf7n5PNhM3u1rsSjloC8cGT6tEVAIgZeODePSlTULIjF71bp6vHw8jPBk6rTHf3cwhMaAyzQ7F87FJ69ZjWQmi2v/+Tk8+Go3PnRZu2mr2RhGFobdUYnoaSLaN8Ofm+Z62wyPzdhiTUS3E9F2Itpeacug2uvAt27ZgitWmzcpPp0rVtciFEniyLSZT4cHoxiOpnDpKnOHqVTetrER2ZzAk/sHph4biSbx/OEh3Li5aUGI38q6Cmz7QAfeem4DbtrchE9dy94Gw5yJYfWdQohrS3hbD4DWaT+3AOib5fjbAGwDgI6ODvHCZ68u4ePkcXle4J7cNzC1wv3vXb0gAi5fIMKxocmP5bVePLa3D+/b2gYAePy1fmRyAu/YPK+jaBquXR/Eteu5iophZsPcMRzgVwBuISInES0HsBrAnyTbpAvNlW5ctbYO3/vDcYzH0ghPpvCjl07ihk1NaKo0f4gHUDqwb9jUiJeOjUxViD2yqxfrGnxTHfIMwyx8TCEcRPROIuoBcAmAx4noSQAQQuwH8AsArwP4DYC/FEIUPvt3gfGZ69Yhkszga08dxDeeOoR4OotPXr1KtllFcdPmZhARPv6zXfj2749iZ9cY3rll4XgbDMPMD801lXWh0tHRIbZv3y7bjJL461/sxsM7ewEAN21uwjdv2SLZouL51Z4+fOrBXRBCyXv8883nmXJOGMMwp0NEO4QQHfO9zrwzLJYo/+/GDXjz+gb43TZcWOQYALNw43lNsFsIJ0Ym8bErV8JiMX9SnGGYwmHhMBl+lx3Xndsg24yyeevGRtkmMAyjE6bIcTAMwzALBxYOhmEYpihYOBiGYZiiYOFgGIZhioKFg2EYhikKFg6GYRimKFg4GIZhmKJg4WAYhmGKYlGOHCGiCICZNoaSTS2As3drkosZbQLYrmIwo00A21UMZrFpmRBi3r0oFmvn+KFC5q0YDRFtN5tdZrQJYLuKwYw2AWxXMZjRprngUBXDMAxTFCwcDMMwTFEsVuHYJtuAWTCjXWa0CWC7isGMNgFsVzGY0aZZWZTJcYZhGEY/FqvHwTAMw+jEohIOIrqOiA4R0VEi+pxsewCAiFqJ6FkiOkBE+4noU7Jtmg4RWYloFxE9JtsWFSKqJKJfEtHB/O/tEhPY9Ff587ePiP6TiFyS7PgBEYWIaN+0x6qJ6LdEdCT/d5VJ7Ppa/hzuJaL/JqJK2TZNe+4uIhJEVGukTXPZRUSfyN+/9hPRPxltVzEsGuEgIiuAfwPwVgDrAbyPiNbLtQoAkAHwaSHEOQAuBvCXJrFL5VMADsg24gy+CeA3Qoh1AM6DZPuIqBnAJwF0CCHOBWAFcIskcx4AcN0Zj30OwDNCiNUAnsn/bDQP4Gy7fgvgXCHEJgCHAdxtAptARK0A3gSgy2B7VB7AGXYR0VUAbgKwSQixAcDXJdhVMItGOABsBXBUCHFcCJEC8CCUEyEVIUS/EGJn/t8RKDfBZrlWKRBRC4C3Afi+bFtUiMgP4EoA/w4AQoiUEGJMrlUAlJ4nNxHZAHgA9MkwQgjxPIDwGQ/fBOA/8v/+DwDvMNQozGyXEOIpIUQm/+PLAFpk25TnXwB8BoCUBO8sdt0B4B+FEMn8a0KGG1YEi0k4mgF0T/u5Bya5QasQUTuALQBekWvJFPdC+QLlZBsyjRUAhgD8MB9C+z4ReWUaJITohbIC7ALQD2BcCPGUTJvOICiE6AeUhQqAesn2zMSHAPxathFEdCOAXiHEHtm2nMEaAFcQ0StE9BwRXSjboLlYTMJBMzxmmpIxIqoA8BCAO4UQEyaw5wYAISHEDtm2nIENwPkA7hdCbAEwCTmhlynyOYObACwH0ATAS0S3yrRpIUFE90AJ2f5Ush0eAPcA+FuZdsyCDUAVlHD2/wHwCyKa6Z5mChaTcPQAaJ32cwskhRPOhIjsUETjp0KIh2Xbk+cyADcS0UkoYb2riegnck0CoJzHHiGE6pX9EoqQyORaACeEEENCiDSAhwFcKtmm6QwSUSMA5P82TZiDiD4I4AYAfybk1/6vhCL+e/LXfQuAnUTUINUqhR4ADwuFP0GJAhieuC+UxSQcrwJYTUTLicgBJXn5K8k2Ib9q+HcAB4QQ/yzbHhUhxN1CiBYhRDuU39XvhBDSV9FCiAEA3US0Nv/QNQBel2gSoISoLiYiT/58XgNzFRT8CsAH8//+IID/kWjLFER0HYDPArhRCBGTbY8Q4jUhRL0Qoj1/3fcAOD9/zcnmEQBXAwARrQHggDmGHs7IohGOfBLu4wCehPKl/oUQYr9cqwAoK/sPQFnR787/uV62USbnEwB+SkR7AWwG8GWZxuS9n18C2AngNSjfGymdvkT0nwBeArCWiHqI6MMA/hHAm4joCJRqoX80iV33AfAB+G3+uv+OCWySzix2/QDAinyJ7oMAPmgCD21WuHOcYRiGKYpF43EwDMMwxsDCwTAMwxQFCwfDMAxTFCwcDMMwTFGwcDAMwzBFwcLBMAzDFAULB8MwDFMULBwMIxEi+lci2mn2oXYMMx0WDoaRRH7qbz2Aj0KZ58QwCwIWDoYxACJy58dlW9XHhBCTABoB/B7At/KvcxDR8/l9PxjGlLBwMIwxfAjK9NOs+gAR1UDZFCoCIAsoG1dB2cXvZhlGMkwhsHAwTJkQ0cemDbA8QUTPzvCyP8PZU2u/AGWDqP1QtjtWeST/eoYxJSwcDFMmQojvCCE2A7gQyqju08bn58f8rxBCnJz2WDuUPT1+DmWa84Zpb9mXPxbDmBIWDobRjm9C2dfk0TMerwVw5r7pXwTw9/nR2acJRz6clSIin57GMkypcAKOYTSAiG4DsAzKnjBnEgfgmvbazQDeBeByIvq3/HOvnfEeJ4CELsYyTJmwcDBMmRDRBQDuAnCFECJ35vNCiFEishKRSwiRAPBVAG8XQjyTf38QwK5px6sBoG5TyzCmg4WDYcrn4wCqATyr7CyL7UKIj5zxmqegeBg5AF5VNABACDFIRF4iqhZChAFcBeAJg2xnmKLhHQAZxgCIaAuAvxZCfKCA1z4M4G4hxCH9LWOY4uHkOMMYgBBiFxSPxDrX6/IVWI+waDBmhj0OhmEYpijY42AYhmGKgoWDYRiGKQoWDoZhGKYoWDgYhmGYomDhYBiGYYqChYNhGIYpChYOhmEYpij+P5LI2ZuPfDgCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2ab9e8daaf28>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "if ham_vasp.status.finished:\n",
    "    # Get the electrostatic potential\n",
    "    epot = ham_vasp.get_electrostatic_potential()\n",
    "    \n",
    "    # Compute the lateral average along the z-axis (ind=2)\n",
    "    epot_z = epot.get_average_along_axis(ind=2)\n",
    "    \n",
    "    # Get the final relaxed structure from the simulation\n",
    "    struct = ham_vasp.get_structure(iteration_step=-1)\n",
    "    r = np.linalg.norm(struct.cell[2])\n",
    "    z = np.linspace(0, r, len(epot_z))\n",
    "    \n",
    "    # Computing the vacuum-level\n",
    "    vac_level = np.max(epot_z)\n",
    "    \n",
    "    # Get the electronic structure\n",
    "    es = ham_vasp.get_electronic_structure()\n",
    "    print(\"wf:\", vac_level - es.efermi)\n",
    "    plt.plot(z, epot_z - vac_level)\n",
    "    plt.xlim(0, r)\n",
    "    plt.axhline(es.efermi - vac_level,\n",
    "                color=\"black\", \n",
    "                linestyle=\"dashed\")\n",
    "    plt.xlabel(\"z ($\\AA$)\")\n",
    "    plt.ylabel(\"V - V$_{vac}$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Looping over a series of hcp(0001) surfaces\n",
    "\n",
    "We now repeat the workflow for a set of hcp metals (the chosen lattice parameters are approximate). Note that if you use the same naming convention, pyiron detects that a job with the same name exists (\"Mg_0001\") and loads the output from this calculation rather than launch a new job with the same name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "hcp_dict = {\"Zn\": {\"a\":2.6649, \"c\": 4.9468}, \n",
    "            \"Mg\": {\"a\": 3.1919, \"c\": 5.1852}, \n",
    "            \"Co\": {\"a\": 2.5071 , \"c\": 4.0695},\n",
    "            \"Ru\": {\"a\": 2.7059 , \"c\": 4.2815}}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "vac = 10\n",
    "size = (2, 2, 4)\n",
    "for element, lattice_parameters in hcp_dict.items():\n",
    "    surf = pr.create_surface(element, \n",
    "                             surface_type=\"hcp0001\", \n",
    "                             size=size, \n",
    "                             a=lattice_parameters[\"a\"], \n",
    "                             c=lattice_parameters[\"c\"], \n",
    "                             orthogonal=True, vacuum=vac)\n",
    "    surf.add_tag(selective_dynamics=[False, False, False])\n",
    "    pos_z = surf.positions[:, 2]\n",
    "    z_min, z_max = np.min(pos_z), np.max(pos_z)\n",
    "    eps = 1e-4\n",
    "    relax_indices = np.argwhere(((pos_z - eps) > z_min) \n",
    "                                & ((pos_z + eps) < z_max ))\n",
    "    relax_indices  = relax_indices.flatten()\n",
    "    surf.selective_dynamics[relax_indices] = [True, True, True]\n",
    "    job_name = \"{}_0001\".format(element)\n",
    "    ham = get_ham(pr, surf, \n",
    "                  name=job_name,\n",
    "                  sigma=0.1, \n",
    "                  mesh=\"GP\", \n",
    "                  encut=350)\n",
    "    #ham.server.cores = 20\n",
    "    #ham.server.queue = queue\n",
    "    ham.executable.version = '5.4_gamma'\n",
    "    ham.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading and analyzing\n",
    "Now we iterate over all jobs in this project and calculate the workfunction. We also time how long the cell takes to execute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time: 9.250723838806152s\n"
     ]
    }
   ],
   "source": [
    "t1 = time.time()\n",
    "for ham in pr.iter_jobs():\n",
    "    if ham.status.finished:\n",
    "            final_struct = ham.get_structure(iteration_step=-1)\n",
    "            elec_structure = ham.get_electronic_structure()\n",
    "            e_Fermi = elec_structure.efermi\n",
    "            epot = ham.get_electrostatic_potential()\n",
    "            epot_z = epot.get_average_along_axis(ind=2)\n",
    "            vacuum_level = np.max(epot_z)\n",
    "            wf = vacuum_level - e_Fermi\n",
    "            element = final_struct.get_majority_species()[-1]\n",
    "            hcp_dict[element][\"work_func\"] = wf\n",
    "t2 = time.time()\n",
    "print(\"time: {}s\".format(t2-t1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compiling data in a table using [pandas](https://pandas.pydata.org/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    a [A]  c [A]  wf [eV]\n",
      "Co  2.507  4.069    5.569\n",
      "Mg  3.192  5.185    3.373\n",
      "Ru  2.706  4.282    5.305\n",
      "Zn  2.665  4.947    3.603\n"
     ]
    }
   ],
   "source": [
    "df = pd.DataFrame(hcp_dict).T\n",
    "df = df.rename(columns={'a': 'a [A]', \n",
    "                        'c': 'c [A]', \n",
    "                        'work_func': 'wf [eV]'})\n",
    "print(df.round(3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.5.4"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}