WorkshopPotentialEAM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ranking-inside",
   "metadata": {},
   "source": [
    "# Fitting an EAM potential\n",
    "EAM potentials are pair functionals. \n",
    "In a generalised form they are equal to Finnis-Sinclair, effective medium theory or glue potentials. Their total energy can be written as\n",
    "\n",
    "$E = \\sum_{ij}V(r_{ij}) + \\sum_i F(\\rho_i)$\n",
    "\n",
    "with\n",
    "\n",
    "$\\rho_i = \\sum_j \\rho(r_{ij})$\n",
    "\n",
    "The original functions for V, $\\rho$ and F were derived from different theories, but they can be chosen freely.\n",
    "\n",
    "Fitting is done using atomicrex https://atomicrex.org. In the fit process an objective or cost function is minimized. The objective function is defined as\n",
    "\n",
    "$\\chi^2 = \\sum_i w_i r_i$\n",
    "\n",
    "where $w_i$ is a weight and $r_i$ is a residual that describes the difference to target values. This residual can be defined in different ways, so it is not possible to simply compare the residual for different fitting processes or codes. A more in depth explanation and some examples can be found on https://atomicrex.org/overview.html#objective-function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "honey-element",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from pyiron import Project, ase_to_pyiron"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "governing-madagascar",
   "metadata": {},
   "source": [
    "### Import the training data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "constant-respect",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>id</th>\n",
       "      <th>status</th>\n",
       "      <th>chemicalformula</th>\n",
       "      <th>job</th>\n",
       "      <th>subjob</th>\n",
       "      <th>projectpath</th>\n",
       "      <th>project</th>\n",
       "      <th>timestart</th>\n",
       "      <th>timestop</th>\n",
       "      <th>totalcputime</th>\n",
       "      <th>computer</th>\n",
       "      <th>hamilton</th>\n",
       "      <th>hamversion</th>\n",
       "      <th>parentid</th>\n",
       "      <th>masterid</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>759</td>\n",
       "      <td>finished</td>\n",
       "      <td>None</td>\n",
       "      <td>df1_A1_A2_A3_EV_elast_phon</td>\n",
       "      <td>/df1_A1_A2_A3_EV_elast_phon</td>\n",
       "      <td>/home/niklas/pyiron/projects/</td>\n",
       "      <td>pyiron_potentialfit/datasets/imported_datasets/Cu_database/</td>\n",
       "      <td>2021-02-08 10:33:52.341472</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>zora@cmti001#1</td>\n",
       "      <td>TrainingContainer</td>\n",
       "      <td>0.4</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>760</td>\n",
       "      <td>finished</td>\n",
       "      <td>None</td>\n",
       "      <td>df3_10k</td>\n",
       "      <td>/df3_10k</td>\n",
       "      <td>/home/niklas/pyiron/projects/</td>\n",
       "      <td>pyiron_potentialfit/datasets/imported_datasets/Cu_database/</td>\n",
       "      <td>2021-02-08 10:33:53.993230</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>zora@cmti001#1</td>\n",
       "      <td>TrainingContainer</td>\n",
       "      <td>0.4</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>761</td>\n",
       "      <td>finished</td>\n",
       "      <td>None</td>\n",
       "      <td>df2_1k</td>\n",
       "      <td>/df2_1k</td>\n",
       "      <td>/home/niklas/pyiron/projects/</td>\n",
       "      <td>pyiron_potentialfit/datasets/imported_datasets/Cu_database/</td>\n",
       "      <td>2021-02-08 10:33:54.435308</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>zora@cmti001#1</td>\n",
       "      <td>TrainingContainer</td>\n",
       "      <td>0.4</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    id    status chemicalformula                         job  \\\n",
       "0  759  finished            None  df1_A1_A2_A3_EV_elast_phon   \n",
       "1  760  finished            None                     df3_10k   \n",
       "2  761  finished            None                      df2_1k   \n",
       "\n",
       "                        subjob                    projectpath  \\\n",
       "0  /df1_A1_A2_A3_EV_elast_phon  /home/niklas/pyiron/projects/   \n",
       "1                     /df3_10k  /home/niklas/pyiron/projects/   \n",
       "2                      /df2_1k  /home/niklas/pyiron/projects/   \n",
       "\n",
       "                                                       project  \\\n",
       "0  pyiron_potentialfit/datasets/imported_datasets/Cu_database/   \n",
       "1  pyiron_potentialfit/datasets/imported_datasets/Cu_database/   \n",
       "2  pyiron_potentialfit/datasets/imported_datasets/Cu_database/   \n",
       "\n",
       "                   timestart timestop totalcputime        computer  \\\n",
       "0 2021-02-08 10:33:52.341472     None         None  zora@cmti001#1   \n",
       "1 2021-02-08 10:33:53.993230     None         None  zora@cmti001#1   \n",
       "2 2021-02-08 10:33:54.435308     None         None  zora@cmti001#1   \n",
       "\n",
       "            hamilton hamversion parentid masterid  \n",
       "0  TrainingContainer        0.4     None     None  \n",
       "1  TrainingContainer        0.4     None     None  \n",
       "2  TrainingContainer        0.4     None     None  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_pr = Project(\"../../datasets/imported_datasets/\")\n",
    "if len(data_pr.job_table()) == 0:\n",
    "    data_pr.unpack(\"../../datasets/Cu_training_archive\")\n",
    "data_pr.job_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "dirty-measurement",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_job = data_pr.load(\"df1_A1_A2_A3_EV_elast_phon\")\n",
    "df = data_job.to_pandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "referenced-julian",
   "metadata": {},
   "outputs": [
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "Are you sure you want to delete all jobs from 'WorkshopPotential'? y/(n) y\n"
     ]
    }
   ],
   "source": [
    "pr = Project(\"WorkshopPotential\")\n",
    "pr.remove_jobs()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "voluntary-limit",
   "metadata": {},
   "source": [
    "### Create an atomicrex job"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "entertaining-jacksonville",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = pr.create_job(pr.job_type.Atomicrex, \"PotentialDF1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "raising-clear",
   "metadata": {},
   "source": [
    "### Add the structures that should be fitted.\n",
    "It is possible to assign different weights to certain structures or properties, depending on what should be investigated using the potential. Here every structure has the same weight, but the force vector with N*3 values is normalized to have the same total weight as the single value energy. Therefore it is divided by the number of atoms. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "located-individual",
   "metadata": {},
   "outputs": [],
   "source": [
    "for id, row in df.iterrows():\n",
    "    struct = ase_to_pyiron(row.atoms)\n",
    "    s = job.structures.add_structure(struct, f\"id{id}\", relative_weight=1)\n",
    "    s.fit_properties.add_FitProperty(\"atomic-energy\", target_value=row.energy/row.number_of_atoms, relative_weight=1)\n",
    "    s.fit_properties.add_FitProperty(\"atomic-forces\", target_value=row.forces, relative_weight=1/row.number_of_atoms)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "angry-leader",
   "metadata": {},
   "source": [
    "### Define the type of potential and necessary functions.\n",
    "In this case an eam potential is fitted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "functional-formation",
   "metadata": {},
   "outputs": [],
   "source": [
    "job.potential = job.factories.potentials.eam_potential()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "realistic-karaoke",
   "metadata": {},
   "source": [
    "It is necessary to define a pair potential, an electronic density function and an embedding function.\n",
    "For all of those it is possible to choose between different functional forms.\n",
    "Classic pair potentials are physically motivated and  have a very limited number of paramaters that are derived from a experimentally measured quantity.\n",
    "Splines or polynomials offer more flexibility, but can lead to unphysical oscillations or overfitting. Compared with the machine learning potentials shown later the number of parameters is very low no matter which functions you choose and the problem is highly non linear.\n",
    "\n",
    "In this case a generalized morse function is used for the pair interaction. It has the form\n",
    "\n",
    "$(\\frac{D_0}{S-1}exp(-\\beta \\sqrt{2S}(r-r_0))-\\frac{D_0S}{S-1}exp(-\\beta\\sqrt{2/S}(r-r_0)))+\\delta $\n",
    "\n",
    "The parameters in the morse potential can be derived from phyiscal quantities, but in this case they are just educated guesses. For example $r_0$ is the equilibrium distance of a dimer. The nearest neighbor distance in fcc Cu is about 2.5 $\\mathring A$ so it is taken as initial value.\n",
    "In the case of analytic functions the initial parameter choices should not matter too much, since the functional form is constrained.\n",
    "\n",
    "The electronic density and embedding function will be splines. Depending on the properties that are calculated other functional forms could give better results. The inital parameters require more testing and hand tuning than the parameters of analytic functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "interpreted-orange",
   "metadata": {},
   "outputs": [],
   "source": [
    "V = job.factories.functions.morse_B(identifier=\"V_CuCu\", D0=0.35, r0=2.5, beta=2, S=2, delta=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "mathematical-gasoline",
   "metadata": {},
   "outputs": [],
   "source": [
    "V.parameters.D0.min_val = 0\n",
    "V.parameters.D0.max_val = 2\n",
    "V.parameters.r0.min_val = 1.5\n",
    "V.parameters.r0.max_val = 3.0\n",
    "V.parameters.S.min_val = 1.1\n",
    "V.parameters.S.max_val = 10.0\n",
    "V.parameters.delta.min_val = -1\n",
    "V.parameters.delta.max_val = 1\n",
    "V.parameters.beta.min_val = 0.1\n",
    "V.parameters.beta.max_val = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "written-commission",
   "metadata": {},
   "source": [
    "Additionally a screening function needs to be defined for the morse potential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "discrete-terminology",
   "metadata": {},
   "outputs": [],
   "source": [
    "V.screening  = job.factories.functions.exp_A_screening(identifier=\"V_cutoff\", cutoff=7)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "wireless-parts",
   "metadata": {},
   "source": [
    "The electron density is chosen to be a spline function. The cutoff has to be defined. Derivatives left and right are optional, they default to 0. For the right cutoff this is fine, since the forces should smoothly go to 0. For the left this is not necessarily the best choice, since the function value should increase at very close distances. Very large absolute values will lead to osciallations and should be avoided."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "authentic-expression",
   "metadata": {},
   "outputs": [],
   "source": [
    "rho = job.factories.functions.spline(identifier=\"rho_CuCu\", cutoff=7, derivative_left=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bored-afternoon",
   "metadata": {},
   "source": [
    "For a spline function it is necessary to define node points. They can be equally spaced or sampled with higher density around turning points, f.e. the first neighbor distance.\n",
    "Too few nodepoints mean low flexibilty, too many lead to overfitting. This requires some testing to find an optimal choice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "hidden-wildlife",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.5 , 1.58, 2.67, 3.75, 4.83, 5.92, 7.  ])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rho_nodes = np.linspace(0.5, 7.0, 7).round(2)\n",
    "rho_nodes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "binary-devil",
   "metadata": {},
   "source": [
    "The nodes need initial values. The electron density should be proportional to $e^{-r}$, so this function is chosen to calculate them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "comparative-brush",
   "metadata": {},
   "outputs": [],
   "source": [
    "decaying_exp = lambda r: np.exp(-r)\n",
    "rho_initial = decaying_exp(rho_nodes)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "gentle-infrastructure",
   "metadata": {},
   "source": [
    "Additionally it is a good idea to define limits for the node points. This is optional for local minimizers, but the fit can quickly run away without limits. Global optimizers typically require them to constrain the sampled space.\n",
    "\n",
    "A density can't be negative so the lower limit is set to 0. The upper limit is chosen to be 3 times the initial values. These choices aswell as the choice for $e^{-r}$ as initial values are somewhat arbitrary, but don't matter much. The electron density from single atoms does not directly influence the calculated energies and forces, instead the summed up density at some place is used in the embedding function, so the final numerical values are an interplay between electron density and embedding function. Since the latter will also be a spline function it can only be defined for a certain range of rho values as node points. Therefore it is better to limit the range of electron density values and define larger limits for the embedding function instead. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "funny-trinidad",
   "metadata": {},
   "outputs": [],
   "source": [
    "rho_mins = np.zeros((len(rho_nodes)))\n",
    "rho_maxs = 3*rho_initial.round(6)\n",
    "rho.parameters.create_from_arrays(rho_nodes, rho_initial, min_vals=rho_mins, max_vals=rho_maxs)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "promising-draft",
   "metadata": {},
   "source": [
    "Finally the last node point at the cutoff range is set to 0 and fitting is disabled to prevent a discontinuous change of energy at the cutoff."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "mexican-absence",
   "metadata": {},
   "outputs": [],
   "source": [
    "rho.parameters[\"node_7.0\"].start_val = 0\n",
    "rho.parameters[\"node_7.0\"].enabled = False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "standard-relative",
   "metadata": {},
   "source": [
    "$-\\sqrt(\\rho)$ can be used as initial guess for the embedding energy, which is taken from second moment approximation tight binding. \n",
    "The node points have to be chosen in a range compatible to the electron density. This can be estimated by calculating it for a densely packed structure.\n",
    "Alternatively atomicrex writes the maximum electron density of all structures to the output. This can be used as a hint for the node points for consequent fits.\n",
    "Everything else is similar to the electron density."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "large-rating",
   "metadata": {},
   "outputs": [],
   "source": [
    "F = job.factories.functions.spline(identifier=\"F_CuCu\", cutoff=5)\n",
    "F_nodes = np.linspace(0.0, 5.0, 7).round(2) #9 is worse 11 is worse 7 is best\n",
    "F_init = -np.sqrt(F_nodes)\n",
    "F_maxs = np.zeros(len(F_nodes))\n",
    "F_mins = -np.ones(len(F_nodes))*5\n",
    "F.parameters.create_from_arrays(F_nodes, F_init, F_mins, F_maxs)\n",
    "F.parameters[\"node_0.0\"].enabled=False\n",
    "F.parameters[\"node_0.0\"].start_val = 0\n",
    "F.parameters[\"node_5.0\"].enabled=False\n",
    "F.parameters[\"node_5.0\"].start_val = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "several-mercy",
   "metadata": {},
   "source": [
    "The functions have to be assigned to the potential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "heavy-acoustic",
   "metadata": {},
   "outputs": [],
   "source": [
    "job.potential.pair_interactions[V.identifier] = V\n",
    "job.potential.electron_densities[rho.identifier] = rho\n",
    "job.potential.embedding_energies[F.identifier] = F"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "alien-chancellor",
   "metadata": {},
   "source": [
    "### Define fitting procedure\n",
    "Finally a few parameters need to be set that influence the fitting process.\n",
    "The minimization can be done with different algorithms. Atomicrex itself implements the BFGS algorithm. Additionally the algorithms from the nlopt library can be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "enormous-segment",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The job PotentialDF1 was saved and received the ID: 819\n"
     ]
    },
    {
     "data": {
      "application/json": {
       "error": "None",
       "iterations": "array([   1,    2,    3, ..., 1998, 1999, 2000], dtype=uint32)",
       "residual": "array([1.39371e+03, 1.39371e+03, 1.39371e+03, ..., 1.52231e-01,\n       1.52231e-01, 1.52231e-01])"
      },
      "text/plain": [
       "Output({'error': None, 'residual': array([1.39371e+03, 1.39371e+03, 1.39371e+03, ..., 1.52231e-01,\n",
       "       1.52231e-01, 1.52231e-01]), 'iterations': array([   1,    2,    3, ..., 1998, 1999, 2000], dtype=uint32)})"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## Define the atom types of the potential\n",
    "job.input.atom_types.Cu = None\n",
    "## Limited number of steps for the workshop\n",
    "job.input.fit_algorithm = job.factories.algorithms.ar_lbfgs(max_iter=2000)\n",
    "job.run()\n",
    "job.output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "vanilla-chocolate",
   "metadata": {},
   "source": [
    "Plot the resiudal over steps to see how the calculation converges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "aging-backing",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fdb336c5ee0>]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAWUUlEQVR4nO3dcZBdZ33e8e8jrS2wwWDjtaNIMpI9KiCTEsyOa0pC6TjENiGW2wwZeUjRJJ5RwzgNtM0kVj1T6Ew1Q5omTWhrZwS4iIZgVAJjtQkpHiUpacfGWRsbWxaOZQS2kLAWMwEDRrbkX/+4R9qrvbsr7d7du/a538945577nvfc89PZ9bNn33Pfe1JVSJKGw7KlLkCSNDiGviQNEUNfkoaIoS9JQ8TQl6QhMrLUBZzK+eefX2vXrl3qMiTpReXee+/9dlWNTm1/wYf+2rVrGR8fX+oyJOlFJck3pmt3eEeShoihL0lDxNCXpCFi6EvSEDll6Ce5LcnhJA9Ns+43klSS87vatibZl+SRJFd1tb8pyYPNug8nycL9MyRJp+N0zvQ/Dlw9tTHJGuDtwONdbRuATcClzTa3JFnerL4V2AKsb756XlOStLhOGfpV9UXgO9Os+k/AbwLdH9O5Ebi9qo5U1X5gH3B5kpXAOVV1V3U+1vMTwHX9Fi9Jmpt5jeknuRb4ZlU9MGXVKuCJrucHmrZVzfLU9plef0uS8STjExMT8ymRj/+//ex64OC8tpWktprz5KwkZwE3Az873epp2mqW9mlV1XZgO8DY2Ni8PvD/g//zYZYFrn3Dj89nc0lqpfmc6V8CrAMeSPJ1YDVwX5Ifo3MGv6ar72rgYNO+epr2RfPet13CyDLfnCRJ3eacilX1YFVdUFVrq2otnUC/rKq+BewCNiVZkWQdnQu291TVIeDpJFc079p5D3DHwv0zJEmn43Tesvkp4C7gNUkOJLlhpr5VtQfYCTwM/DlwY1Uda1a/F/gonYu7jwGf77N2SdIcnXJMv6quP8X6tVOebwO2TdNvHHj9HOvrS8182UCShlJrB72d+SVJvVob+pKkXoa+JA0RQ1+ShkirQ7+8jitJJ2lt6PsZnpLUq7WhL0nqZehL0hBpdeg7pC9JJ2tt6MfpWZLUo7WhL0nqZehL0hAx9CVpiLQ69MvZWZJ0ktaGvpOzJKlXa0NfktTL0JekIWLoS9IQaXXoexlXkk7W2tD3Oq4k9Tpl6Ce5LcnhJA91tf1Okq8m+UqSzyV5Zde6rUn2JXkkyVVd7W9K8mCz7sOJ76+RpEE7nTP9jwNXT2m7E3h9Vf194G+BrQBJNgCbgEubbW5JsrzZ5lZgC7C++Zr6mpKkRXbK0K+qLwLfmdL2hao62jy9G1jdLG8Ebq+qI1W1H9gHXJ5kJXBOVd1VnRlTnwCuW6B/wyy1L/YeJOnFZSHG9H8F+HyzvAp4omvdgaZtVbM8tX1aSbYkGU8yPjExMb+qHD2SpB59hX6Sm4GjwCePN03TrWZpn1ZVba+qsaoaGx0d7adESVKXkflumGQz8E7gypr8kJsDwJqubquBg0376mnaJUkDNK8z/SRXA78FXFtVP+xatQvYlGRFknV0LtjeU1WHgKeTXNG8a+c9wB191i5JmqNTnukn+RTwNuD8JAeAD9B5t84K4M7mnZd3V9WvVtWeJDuBh+kM+9xYVceal3ovnXcCvZTONYDPI0kaqFOGflVdP03zx2bpvw3YNk37OPD6OVXXBy/jSlKv1s7IlST1MvQlaYi0PvS9e5YkTWpt6Ds3S5J6tTb0JUm9DH1JGiKGviQNkdaHvtdxJWlSa0M/Ts+SpB6tDX1JUi9DX5KGSOtD3yF9SZrU2tB3cpYk9Wpt6EuSehn6kjREDH1JGiKtD30/ZVOSJrU29L2OK0m9Whv6kqRepwz9JLclOZzkoa6285LcmeTR5vHcrnVbk+xL8kiSq7ra35TkwWbdhxPfVClJg3Y6Z/ofB66e0nYTsLuq1gO7m+ck2QBsAi5ttrklyfJmm1uBLcD65mvqay4KR/QladIpQ7+qvgh8Z0rzRmBHs7wDuK6r/faqOlJV+4F9wOVJVgLnVNVd1bmy+omubRaFf0dIUq/5julfWFWHAJrHC5r2VcATXf0ONG2rmuWp7dNKsiXJeJLxiYmJeZYoSZpqoS/kTnd+XbO0T6uqtlfVWFWNjY6OLlhxkjTs5hv6TzZDNjSPh5v2A8Carn6rgYNN++pp2iVJAzTf0N8FbG6WNwN3dLVvSrIiyTo6F2zvaYaAnk5yRfOunfd0bbOonJslSZNGTtUhyaeAtwHnJzkAfAD4ELAzyQ3A48C7AKpqT5KdwMPAUeDGqjrWvNR76bwT6KXA55uvReM7QiWp1ylDv6qun2HVlTP03wZsm6Z9HHj9nKqTJC0oZ+RK0hAx9CVpiLQ+9Ms5uZJ0QutDX5I0ydCXpCFi6EvSEGl96Ds5S5ImtTb0nZslSb1aG/qSpF6GviQNEUNfkoaIoS9JQ6S1oZ9p79siScOttaEvSepl6EvSEGl96Ds5S5ImtTb0nZwlSb1aG/qSpF6GviQNkb5CP8m/TLInyUNJPpXkJUnOS3Jnkkebx3O7+m9Nsi/JI0mu6r98SdJczDv0k6wCfh0Yq6rXA8uBTcBNwO6qWg/sbp6TZEOz/lLgauCWJMv7K//UvHOWJE3qd3hnBHhpkhHgLOAgsBHY0azfAVzXLG8Ebq+qI1W1H9gHXN7n/mfkdVxJ6jXv0K+qbwL/EXgcOAR8t6q+AFxYVYeaPoeAC5pNVgFPdL3EgaatR5ItScaTjE9MTMy3REnSFP0M75xL5+x9HfDjwNlJfmm2TaZpm3bspaq2V9VYVY2Njo7Ot0RJ0hT9DO/8DLC/qiaq6jngs8A/BJ5MshKgeTzc9D8ArOnafjWd4aBF5eQsSZrUT+g/DlyR5KwkAa4E9gK7gM1Nn83AHc3yLmBTkhVJ1gHrgXv62P+snJwlSb1G5rthVX0pyWeA+4CjwJeB7cDLgJ1JbqDzi+FdTf89SXYCDzf9b6yqY33WL0mag3mHPkBVfQD4wJTmI3TO+qfrvw3Y1s8+JUnz54xcSRoirQ99r+NK0qTWhr53zpKkXq0NfUlSL0NfkoaIoS9JQ6T1oV9OyZWkE1ob+s7IlaRerQ19SVIvQ1+ShkjrQ98RfUma1PrQlyRNMvQlaYgY+pI0RAx9SRoirQ9952ZJ0qTWhn6cnSVJPVob+pKkXoa+JA2R9oe+Y/qSdEJfoZ/klUk+k+SrSfYmeXOS85LcmeTR5vHcrv5bk+xL8kiSq/ovf5baFvPFJelFqt8z/T8A/ryqXgu8AdgL3ATsrqr1wO7mOUk2AJuAS4GrgVuSLO9z/5KkOZh36Cc5B3gr8DGAqnq2qv4O2AjsaLrtAK5rljcCt1fVkaraD+wDLp/v/iVJc9fPmf7FwATw35J8OclHk5wNXFhVhwCaxwua/quAJ7q2P9C09UiyJcl4kvGJiYk+SpQkdesn9EeAy4Bbq+qNwA9ohnJmMN0w+7SXWatqe1WNVdXY6OhoHyVCeSVXkk7oJ/QPAAeq6kvN88/Q+SXwZJKVAM3j4a7+a7q2Xw0c7GP/s3JuliT1mnfoV9W3gCeSvKZpuhJ4GNgFbG7aNgN3NMu7gE1JViRZB6wH7pnv/iVJczfS5/b/AvhkkjOBrwG/TOcXyc4kNwCPA+8CqKo9SXbS+cVwFLixqo71uX9J0hz0FfpVdT8wNs2qK2fovw3Y1s8+58oPXJOkSa2dkeuQviT1am3oS5J6GfqSNEQMfUkaIq0Pfa/jStKk1oa+d86SpF6tDX1JUi9DX5KGSOtDv5ydJUkntDb0HdKXpF6tDX1JUi9DX5KGiKEvSUOk9aHvZVxJmtTa0Pc6riT1am3oS5J6GfqSNEQMfUkaIq0PfSfkStKkvkM/yfIkX07yv5rn5yW5M8mjzeO5XX23JtmX5JEkV/W771MUtqgvL0kvRgtxpv8+YG/X85uA3VW1HtjdPCfJBmATcClwNXBLkuULsH9J0mnqK/STrAZ+DvhoV/NGYEezvAO4rqv99qo6UlX7gX3A5f3sX5I0N/2e6f8+8JvA811tF1bVIYDm8YKmfRXwRFe/A01bjyRbkownGZ+YmOirwHJ6liSdMO/QT/JO4HBV3Xu6m0zTNm0iV9X2qhqrqrHR0dH51TevrSSp3Ub62PYtwLVJ3gG8BDgnyR8BTyZZWVWHkqwEDjf9DwBrurZfDRzsY/+SpDma95l+VW2tqtVVtZbOBdq/qKpfAnYBm5tum4E7muVdwKYkK5KsA9YD98y7cknSnPVzpj+TDwE7k9wAPA68C6Cq9iTZCTwMHAVurKpji7B/SdIMFiT0q+qvgL9qlp8Crpyh3zZg20Ls87R5HVeSTmjtjFznZklSr9aGviSpl6EvSUOk9aHvkL4kTWpt6MfpWZLUo7WhL0nqZehL0hAx9CVpiLQ+9L1zliRNam3oOzlLknq1NvQlSb0MfUkaIq0Pfe+cJUmTWhv6DulLUq/Whr4kqZehL0lDxNCXpCHS+tB3cpYkTWpt6Ds5S5J6zTv0k6xJ8pdJ9ibZk+R9Tft5Se5M8mjzeG7XNluT7EvySJKrFuIfIEk6ff2c6R8F/nVVvQ64ArgxyQbgJmB3Va0HdjfPadZtAi4FrgZuSbK8n+IlSXMz79CvqkNVdV+z/DSwF1gFbAR2NN12ANc1yxuB26vqSFXtB/YBl893/5KkuVuQMf0ka4E3Al8CLqyqQ9D5xQBc0HRbBTzRtdmBpm1ReR1Xkib1HfpJXgb8CfD+qvrebF2naZs2k5NsSTKeZHxiYmJ+dTknV5J69BX6Sc6gE/ifrKrPNs1PJlnZrF8JHG7aDwBrujZfDRyc7nWrantVjVXV2OjoaD8lSpK69PPunQAfA/ZW1e91rdoFbG6WNwN3dLVvSrIiyTpgPXDPfPcvSZq7kT62fQvwz4AHk9zftP0b4EPAziQ3AI8D7wKoqj1JdgIP03nnz41VdayP/Z+WcnaWJJ0w79Cvqv/LzB9meeUM22wDts13n3PikL4k9WjtjFxJUi9DX5KGSOtD/+kfHV3qEiTpBaO1oX/J6NkAfOSLX1viSiTphaO1of+mV5/Hq84+k89++Zsc+u4zS12OJL0gtDb0AW7+udcB8NePfnuJK5GkF4ZWh/7bN1wIwN2PPbXElUjSC0OrQ/9lKzrTEJ55btHngEnSi0KrQz8Jl130Sr5/xHfwSBK0PPQBzjpzxDF9SWq0PvRf/pLOEI+fwSNJQxD6P7H6FQAcOfr8ElciSUuv9aF//GLuN5764RJXIklLr/Whf8noywB490fvZu+h2W7sJUnt1/rQf/PFr+Kf/6OL+fb3n+X6j9zNj3z7pqQh1vrQX7YsbL3mdfzyW9bydz98jmv+4K/Zd/j7S12WJC2J1of+cR/4+Uu59d2Xsf/bP+A//8WjS12OJC2JoQl9gGt+YiUAj3/Hi7qShtNQhT7AP37NKM8867i+pOE0dKF/zkvP4LvPPLfUZUjSkhh46Ce5OskjSfYluWnQ+3/JyHIOffdHfh6PpKE0MsidJVkO/Ffg7cAB4G+S7KqqhwdVw09e9Eo+Pf4Eb/h3X2DrNa/l7BUjrBhZxpkjyzhj+TLOOnM5r/2xcxh9+YpBlSRJAzPQ0AcuB/ZV1dcAktwObAQGFvq/cNlq7nrsKXY9cJB//6d7Z+2bQOh8WueyQOg0pFm3LDmxPgDH27q2O973pGVyYvsT+2nal+Xk1zv+OpKGz5/++k+xYmT5gr7moEN/FfBE1/MDwD+Y2inJFmALwEUXXbSgBZw5sowPX/9GfvcX38APjxzjB88e5cjR5zl67HmePfY833vmKHsOfpfvPfMcBVRBUVTB881y8x9V1Wnr6gPwfNVJbcf71pS+1fTtfr1O2+Qyfk6cNLSa078FNejQn+5f0BNrVbUd2A4wNja2KLF3xvJlvOKsZbzirDN61r35klctxi4lackN+kLuAWBN1/PVwMEB1yBJQ2vQof83wPok65KcCWwCdg24BkkaWgMd3qmqo0l+DfjfwHLgtqraM8gaJGmYDXpMn6r6M+DPBr1fSdIQzsiVpGFm6EvSEDH0JWmIGPqSNERS9cKe8plkAvjGPDc/H/j2ApazUKxrbqxrbqxrbtpa16uranRq4ws+9PuRZLyqxpa6jqmsa26sa26sa26GrS6HdyRpiBj6kjRE2h7625e6gBlY19xY19xY19wMVV2tHtOXJJ2s7Wf6kqQuhr4kDZFWhv5S3nw9yZokf5lkb5I9Sd7XtH8wyTeT3N98vaNrm61NrY8kuWoRa/t6kgeb/Y83becluTPJo83juYOsK8lruo7J/Um+l+T9S3G8ktyW5HCSh7ra5nx8krypOc77knw4fd7vcoa6fifJV5N8JcnnkryyaV+b5Jmu4/aHi1XXLLXN+Xs3oGP26a6avp7k/qZ9IMdslmwY7M9Y5zZ+7fmi85HNjwEXA2cCDwAbBrj/lcBlzfLLgb8FNgAfBH5jmv4bmhpXAOua2pcvUm1fB86f0vYfgJua5ZuA3x50XVO+d98CXr0Uxwt4K3AZ8FA/xwe4B3gznTvFfR64ZhHq+llgpFn+7a661nb3m/I6C1rXLLXN+Xs3iGM2Zf3vAv92kMeMmbNhoD9jbTzTP3Hz9ap6Fjh+8/WBqKpDVXVfs/w0sJfOvYFnshG4vaqOVNV+YB+df8OgbAR2NMs7gOuWsK4rgceqarYZ2ItWV1V9EfjONPs77eOTZCVwTlXdVZ3/Oz/Rtc2C1VVVX6iqo83Tu+nchW5Gi1HXTLXNYkmP2XHNWfEvAp+a7TUWuq5ZsmGgP2NtDP3pbr4+W+gumiRrgTcCX2qafq35c/y2rj/hBllvAV9Icm86N58HuLCqDkHnhxK4YAnqOm4TJ/+PuNTHC+Z+fFY1y4OqD+BX6JztHbcuyZeT/J8kP920DbquuXzvBl3bTwNPVtWjXW0DPWZTsmGgP2NtDP3Tuvn6oheRvAz4E+D9VfU94FbgEuAngUN0/ryEwdb7lqq6DLgGuDHJW2fpO9DjmM7tM68F/kfT9EI4XrOZqY5BH7ebgaPAJ5umQ8BFVfVG4F8Bf5zknAHXNdfv3aC/p9dz8snFQI/ZNNkwY9cZ9t9XXW0M/SW/+XqSM+h8Uz9ZVZ8FqKonq+pYVT0PfITJIYmB1VtVB5vHw8DnmhqebP5cPP7n7OFB19W4Brivqp5salzy49WY6/E5wMlDLYtWX5LNwDuBdzd/5tMMBTzVLN9LZxz47w2yrnl87wZ5zEaAfwp8uqvegR2z6bKBAf+MtTH0l/Tm68144ceAvVX1e13tK7u6/RPg+LsKdgGbkqxIsg5YT+cizULXdXaSlx9fpnMh8KFm/5ubbpuBOwZZV5eTzr6W+nh1mdPxaf48fzrJFc3Pwnu6tlkwSa4Gfgu4tqp+2NU+mmR5s3xxU9fXBlVXs985fe8GWRvwM8BXq+rE8MigjtlM2cCgf8bmeyX6hfwFvIPOlfHHgJsHvO+fovOn1leA+5uvdwD/HXiwad8FrOza5uam1kdYgHdUzFDXxXTeCfAAsOf4cQFeBewGHm0ezxtkXc1+zgKeAl7R1Tbw40Xnl84h4Dk6Z1M3zOf4AGN0gu4x4L/QzHxf4Lr20RnvPf4z9odN319ovr8PAPcBP79Ydc1S25y/d4M4Zk37x4FfndJ3IMeMmbNhoD9jfgyDJA2RNg7vSJJmYOhL0hAx9CVpiBj6kjREDH1JGiKGviQNEUNfkobI/wdEaXp5kEKZcAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(job.output.iterations, job.output.residual)\n",
    "#plt.ylim(0,5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "informed-formula",
   "metadata": {},
   "source": [
    "Finally it is a good idea to have a look at the final potential. This can reveal unphysical behavior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "resident-gnome",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<Figure size 576x1296 with 3 Axes>,\n",
       " array([[<AxesSubplot:title={'center':'Cu F'}, xlabel='$\\\\rho $ [a.u.]'>],\n",
       "        [<AxesSubplot:title={'center':'Cu rho_CuCu'}, xlabel='r [$\\\\AA$]'>],\n",
       "        [<AxesSubplot:title={'center':'Cu V_CuCu'}, xlabel='r [$\\\\AA$]'>]],\n",
       "       dtype=object))"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeIAAAQVCAYAAABuRsmlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABVvklEQVR4nO3deZhc5Xnn/d9dVb0vavWitdVqCQFiR7jNvuMFDAFv48EOtmPHQzKxE/t1Mh479puZZDKZN8mMx+MtDmM7xoEYxwY7DoYYY8DYBgEtIRBCaN9aW3dr632rut8/qiRaTbfUUp/qp6vq+7muvmo5p855DoX62+dU1SlzdwEAgDBioQcAAEAhI8QAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxkOPM7ANm1mpmPWa218weNbOrT3NZbma9mWX1mNnhiIcLYAxCDOQwM/u0pC9J+itJcyU1Sfq6pDumsNiL3L0y81Mz5UECOCHjzFpAbjKzWZJ2S/qIu/9ggnm+I6nN3b+QuX29pPvcvXGC+V3Sme6+ORtjBvBG7BEDuesKSaWSfhR6IABOHyEGcledpE53H4l4uavN7HDm58sRLxvAGInQAwBw2g5IqjezRMQxvoRD08D0YY8YyF3PShqQ9M4TzNMrqXzU7XnZHBCAU0eIgRzl7kck/Zmkr5nZO82s3MyKzOwWM/ubzGxrJL3DzGrNbJ6kTwUaLoAJEGIgh7n7FyV9WtIXJHVI2iXpE5J+nJnlHyW9JGm7pMckfX/aBwnghPj4EgAAAbFHDABAQIQYAICACDEAAAERYgAAAiLEAAAEFOTMWvX19d7c3Bxi1QAATLtVq1Z1unvDeNOChLi5uVmtra0hVg0AwLQzsx0TTePQNAAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAkYXYzOJm9qKZPRzVMgEAyHdR7hF/UtL6CJcHAEDeiyTEZtYo6VZJ34xieQAAFIqo9oi/JOkzklITzWBmd5tZq5m1dnR0RLRaAABy25RDbGa3SWp391Unms/d73H3FndvaWhomOpqAQDIC1HsEV8l6XYz2y7pAUk3mtl9ESwXAIC8N+UQu/vn3L3R3Zsl3SnpCXe/a8ojAwCgAPA5YgAAAkpEuTB3f0rSU1EuEwCAfMYeMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAAU05xGa2yMyeNLP1ZrbOzD4ZxcAAACgEiQiWMSLpj919tZlVSVplZj9391cjWDYAAHltynvE7r7X3VdnrndLWi9p4VSXCwBAIYj0NWIza5a0QtJz40y728xazay1o6MjytUCAJCzIguxmVVKelDSp9y9a+x0d7/H3VvcvaWhoSGq1QIAkNMiCbGZFSkd4fvd/aEolgkAQCGI4l3TJulbkta7+xenPiQAAApHFHvEV0n6oKQbzWxN5ucdESwXAIC8N+WPL7n7ryVZBGMBAKDgcGYtAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAAUUSYjO72cw2mNlmM/tsFMsEAKAQTDnEZhaX9DVJt0g6V9L7zezcqS4XAIBCkIhgGZdK2uzuWyXJzB6QdIekVyNY9km9uqdLm9q7VZKIqyQRU3EipoqShOoqilVXWazy4ig2EQCA7IiiUgsl7Rp1u03SZWNnMrO7Jd0tSU1NTRGsNu2RtXv11Sc3Tzi9rCiuebNK1VRbrua6ci2uq9CShgqdN79ac6pLIxsHAACnI4oQ2zj3+RvucL9H0j2S1NLS8obpp+t3r16id65YqMGRpIZGUhocSalnYEQHe4d0oHdInT2D2ndkQDsO9mr1jkPqHhw59tj6yhKdt6Ba5y+s1psWz9abFtdqVllRVEMDAOCkoghxm6RFo243StoTwXInZXZFsWZXFE9qXnfXwd4hbW7v0bo9XXp1b5fW7enS3/9yq0ZSLjPpnHnVunRJra44o05XLatXZQmHtgEA2WPuU9s5NbOEpI2SbpK0W9ILkj7g7usmekxLS4u3trZOab1R6h9K6sVdh/TCtkN6fvsBrdpxSAPDKRXFTZcuqdUNZ8/RDcvnaGl9hczGOwAAAMDEzGyVu7eMO22qIc6s4B2SviQpLunb7v7fTzT/TAvxWEMjKa3acUhPbWjXkxvatXF/jyRpaUOFbrtgvm69cIHOmltJlAEAk5L1EJ+qmR7isdoO9enJ19r16Cv7tHLrAaVcWjanUrdeMF+3X7xAZzRUhh4iAGAGI8QR6uge1L+t26efvrxHz207KHdpRVON3teySLddOF9VpbzZCwBwPEKcJe1dA/rRi7v1g1Vt2tzeo9KimG45f77+3ZsadcUZdRy6BgBIIsRZ5+5as+uwfrCqTf/60h51D4xoaX2F7rp8sd7b0qhq9pIBoKAR4mk0MJzUI2v36h9X7tCLOw+rvDiud65YqA9dsVjL51WHHh4AIABCHMjatiP67rPb9ZOX9mhwJKVLm2v10auX6K3nzlU8xmFrACgUhDiwQ71D+sGqXfruszvUdqhfzXXl+t1rluq9lzSqrDgeengAgCwjxDPESDKln63br3ue3qKX2o6otqJYH7x8sT50xWLVVZaEHh4AIEsI8Qzj7np+20H9319t1ePr21WSiOm9b2rU3dcu1eK6itDDAwBE7EQh5kTKAZiZLltap8uW1mlze4+++aut+kFrm773/E7dcfFCffyGM7RsTlXoYQIApgF7xDNEe9eA/u+vtuq+lTs1MJLULefP08dvWKbzFswKPTQAwBRxaDqHHOwd0rd/vU33PrNd3YMjumn5HH38xmW6pGl26KEBAE4TIc5BR/qH9d1ntuvbv9mmQ33DumpZnT5xw5m6fGktZ+wCgBxDiHNY7+CI/um5nfr7p7eqs2dQb26erT+66UxdvayeIANAjiDEeWBgOKnvv7BL3/jlFu09MqAVTTX6oxvP1PVnNxBkAJjhCHEeGRxJ6oer2vT1J7do9+F+XbBwlv7wxmV667lzCTIAzFCEOA8NJ1P60erd+uqTm7XzYJ/OmV+tP7pxmd5+3jzFOH0mAMwohDiPjSRT+slLe/TVJzZra2evzppbqU/ceKZuvWA+57MGgBmCEBeAZMr18MvpIG9q79HShgp94oZluv2iBUrEY6GHBwAFjRAXkFTK9W/r9unLv9ik1/Z1a3FduT5+wzK9a8VCFRFkAAiCEBegVMr1+Pr9+vITm/TK7i41zi7TH1y/TO9500KVJPjGJwCYToS4gLm7ntzQri//YrPW7Dqs+bNK9R+vP0Pva1mk0iKCDADTgRBD7q5fberUl3+xSa07DmlOVYl+77oz9IFLm/hOZADIMkKMY9xdz249oC//YpNWbj2o+spi3X3tUv32ZYtVUcKXcQFANhBijOv5bQf1lSc26VebOjW7vEgfu2apPnTFYlWVFoUeGgDkFUKME1q985C+8otNenJDh2aVFemjVy3R71zVrFllBBkAokCIMSkvtx3WV57YrJ+/ul9VJQn9zlXN+uhVSzS7ojj00FAghkZSOtI/rCP9w+oaGNbQSEqplCvprmTKlYjFVFYcU1lRQmXFcVWXJjS7vJizyWHGI8Q4Jev2HNFXn9isR1/Zp9KimN59SaM+elWzls2pCj005DB3V2fPkDa392hbZ6/2HO7XniP92nO4X3uPDKi9a1D9w8lTXm4iZqqvLNGc6hLNqy7VkoYKLa2v0JL6Si1tqFB9ZUkWtgY4NYQYp2Xj/m59+9fb9NCLuzU0ktJ1ZzXod69eomvO5CsYcWIDw0m9urdLa9uO6NU9Xdrc0aPN7T060j98bJ54zDS3qkQLaso0v6ZMc6tKVFNepOqyIs0qK1J1aZFKEjHFYqZ4zBQzU8pdfUNJ9Q8l1T88oiN9w2rvHlR796A6ugfVdqhPOw/2aTj5+u+1udUlumDhLJ2/cJYuWDhLFy+qUR1xxjQjxJiSAz2Duv+5nfruszvU2TOoM+dU6qNXL9G7Vizks8hQKuXa3NGjF7Yf1Nq2I3q57Yg27u/WSCr9u6W2oljL5lSmfxrSl0sbKjSvujQrp18dSaa05/CAtnam4//K7iNau/uItnb26uivuzPnVOrypXW6bGmtLltSp4YqwozsIsSIxOBIUg+/tFff+vU2vbq3S7PLi/TvWhbp/Zc2aUl9RejhYZqMJFNav7dbz207oOe3HdQL2w/qUF96T7emvEgXLJylCxtn6YKFNbpo0SzNqy6dEUdQegZH9OqeLrXuOKiVWw9q1faD6h1KHwpfPq9KNy6foxuXz9HFi2o4PzsiR4gRKXfXc9sO6t5ntuvnr+7XSMp11bI6/fZli/XWc+dyTus84+7a2tmrX27o0NObOtS6/ZB6BkckSYvrynVpc60uXVKrNzfXanFd+YyI7mSMJFN6ZU+Xnt1yQE9taFfrjkNKplw15UW69swGve28ubpx+RyVF/P5ekwdIUbWtHcN6J9bd+l7z+/S7sP9aqgq0b9vWaQ7L12kxtnloYeH09Q7OKJnthzQLze265cbO7TrYL8kaWlDha48o06XLqnTpc21mjerNPBIo3Okf1i/3tSpJ15r1y83tquzZ0ilRTHduHyObr1ggW5Y3kCUcdoIMbIumXL9cmO77l+5U09uaJdLuvKMOr17RaNuPn8eZ+2a4dxdG/f36KkN6fC+sP2ghpOuiuK4rlxWr+vOatB1ZzVoUW1h/HGVTLle2H5QP315rx59ZZ86ewZVVhTXjefM0btXLNS1ZzVw5AenhBBjWu0+3K8ftO7SQ6t3a+fBPpUXx3XL+fP1nksW6vKldXzmc4boGhjWbzZ16pcbO/TLjR3ae2RAUvr10uvOatB1ZzeoZXGtihOFHZxkyvX8toP66do9enTtPh3oHVJ9ZbHeefFCvedNjTpnfnXoISIHEGIE4e5q3XFID65q009f3qvuwREtrCnT7Rcv0K0XzNd5C6pz5vXEfJBKuV7d25UO74YOrdqZfk20qiShq8+s1/VnN+jasxo0f1ZZ6KHOWMPJlJ7a0KEfrtqlJ15r13DSdd6Car3nkka9a8VCTn6DCRFiBDcwnNRjr+7Xg6va9OvNnUqmXM115XrHBfP1DqKcNYd6h/SrzZ16akO7nt7Yqc6eQUnS+QurM4eb52hFUw2HWU/Dwd4h/WTNbj24erfW7j6i4kRMt14wX3dd3qRLmmbz/zOOQ4gxoxzsHdJj6/bpp2v36pktB45F+ZYL5ust58zVxYtqFOfw9WkZSab0UtuRY4ebX247LPf0x4quObNB15/VoGvOqtecqvx5k9VMsH5vl773/E49tHq3egZHtHxelX77sia9c8VCvkQFkggxZrDxolxbUazrz2rQDcvn6NqzGvjyiZPYe6RfT2/s0NMbO/WrTR3qGhhRzKSLFtXo2jPTr/Ve1MgfN9Ohd3BE//rSHt333A69srtL5cVx3XHxAv32ZYt1/sJZoYeHgAgxcsKRvmE9valDT7zWrqc2tOtQ37DiMVPL4tm6elm9rlxWpwsbOYx6pG9Yz28/qGe3HNCvN3do4/4eSelTOV53Vvp13quX1aumnNcrQ3q57bDuX7lT//LSbg0Mp/Tm5tn6nSuX6O3nzeWEIQWIECPnJFOuNbsO6YnX2vXkax16dW+XJKm8OK43N9fqyjPqdMUZdTpnfnXeh7lrYFgvbEuHd+W2A1q3p0vuUnEipjc3zz72Wu9Zcyt5XXIGOtI/rB+07tK9z27XroP9mj+rVB+8YrHe/+Ym3txVQAgxct7B3iE9t/WAntlyQM9uPaDN7em9wJJETBcsnKUVTTW6eNFsXdxUowWzZsYpFU9HKuXa0tGjF3ce1ou7DmvNrsPasK9LKZeK4zGtaKrRFWfU6YqldbpoUQ3n+s4hyZTridfa9Z1ntuk3mw+oJBHTu1Ys1O9c1azl8/gIVL4jxMg77V0DWrntoNbsPKw1uw7plT1dGhpJSZLqK0u0fF6Vzs78LJ9XpTPnVKmseGZFa2A4qc3tPXptX7c27OvSq3u79PKuI+rOnD6yqjShixfVaEXTbF2+tFaXNM0mvHliw75ufeeZ7frRi20aGE7piqV1+shVzbrpnLm8lp+nCDHy3tBISuv3dmnNrsNauzv97T8b93drYDgdZzNpwawyNdWWp3/qyrWotlyNs8vUUFmi+sqSyEPt7uoeHFF716B2HerTzgPpr+jbebBPWzt6tP1An5KZbygqTsR01txKXdSYDu/Fi2q0tL6Ck5/kuUO9Q/p+6y5995nt2nNkQItqy/Shy5v1vjcv4k2KeYYQoyAlU66dB/u0YV+XXtvXrZ0H+rQjE8KO7sE3zF9RHFdDVYlqK4pVWVqkiuK4yosTqihJX8YsHXST6eiR74HhpPqHk+ofSql/eES9g0kd7B3SgZ5BdfYOHdtLP6q0KKam2nItrqvQ8nlVWj6vWmfPq1JzXTlv4ClgI8mUHnt1v/7hN9v0wvZDKi+O6z2XNOrDVzZr2ZzK0MNDBAgxMEbf0IjaDvVr9+F+dXYPqqNnUJ3dQ+rsGdSB3kH1DibVN5QOa+/QiPoGk0q5y5Xe0z36r6Y0EVdZcVxlRenL8uK4aiuKVVdRovqqYtVnLptq03vgDZUlOfv6NabHK7uP6DvPbNdP1uzRUDKla89q0EeuatZ1ZzZwhCSHEWIAyDGdPYP63nM79Y8rd6i9e1BL6yv04Sub9Z43NaqSL1HJOYQYAHLU0EhKj76yV//wm+1as+uwqkoS+ncti/ThKxdrcV1F6OFhkggxAOSBF3ce0nee2a6fvrxXSXfdtHyOPnLVEl15Rh0vecxwhBgA8sj+rgHdv3KH7n9upw70DumsuZX6nSuX6F0rFs64j+khLWshNrO/lfRbkoYkbZH0EXc/fLLHEWIAmLqB4aQefnmv/uE327RuT5dmlRXpzjcv0p2XNmlJPYetZ5Jshvhtkp5w9xEz+2tJcvf/fLLHEWIAiM7R7/7+h99s08/W7Vcy5bpiaZ3ef1mT3n7eXJUk2EsO7UQhntJb79z9sVE3V0p671SWBwA4dWamNzfX6s3NtWrvGtAPVrXpgRd26o++96JmlxfpPZc06s5Lm/hM8gwV2WvEZvavkr7v7vdNMP1uSXdLUlNT05t27NgRyXoBAG+USrl+s6VT33t+px5bt18jKdelzbV6b0ujbjl/Ht+TPM2mdGjazB6XNG+cSZ9393/JzPN5SS2S3u2TKDuHpgFg+nR0D+rB1W164Pmd2n6gTyWJmN567ly9a8VCXXtWQ95/g9lMkNV3TZvZhyX9vqSb3L1vMo8hxAAw/dxdL+46rH95cbf+9eW9Otg7pNqKYt124XzdftECXdI0m7N3ZUk236x1s6QvSrrO3Tsm+zhCDABhDSdT+tWmDv3oxT16bN0+DY6k1FBVoredO1dvP2+eLl9ap+JE4e4pD44kI32TWzZDvFlSiaQDmbtWuvvvn+xxhBgAZo7ugWE9uaFDP3tln57c0K6+oaSqSxO66Zy5umH5HF2zrF6zK4pDDzPrDvUO6adr9+rHL+5WSVFM93/s8siWnc13TS+byuMBAOFVlRbp9osW6PaLFmhgOKlfb+rUv63bp8fX79ePXtwtM+nChbN0zZkNuvasBq1oqsmb15X3HunXY+v262fr9um5bQeVTLnOmlup95zbKHefljOWcWYtAMC4kinXy22H9fTGTj29qUNrdh1WMuUqL47r4kU1ammuVcvi2VrRVJMz78IeGE6qdfsh/WZLp36zuVMvtx2RJJ3RUKG3nzdP77hgvs5bUB15gDnFJQBgyo70D+vZLZ1aufWgXth+UOv3dinlUsyks+dV6/wF1Tp3QbXOnV+tcxZUqzpwnN1d+7oG9NKuw3qp7Yhe3HlIq3ce1tBISomY6eJFNbph+Ry9/bx5Wf+MNSEGAESuZ3BEL+48pBe2H9KLOw/p1T1dOtA7dGx64+wyLamvUHNdhZrrK9RcV66m2nLNqSpVdVkisr3OgeGk9ncNaOfBPm3t6NWWjh5t6ejRxv096ugelCQlYqbl86t0+ZI6XbWsXm9eUjutXyeZtdeIAQCFq7IkoWvObNA1ZzZISu+BdnQPat3eLr26p0uv7evW9s5e/XjnbnUPjhz32OJETA2VJWqoKlFtRbHKi+OqKE6ovCR9GYuZ5K6ju4ojKVfv4Ih6BkeOXXZ0D2p/16CO9A8ft+yqkoSWzqnUtWc26MLGWbqwcZbOmV+t0qKZeapPQgwAiISZaU51qeZUl+qGs+ccu9/ddahvWNsP9GrXwT51dA++/tMzqP1dA+ofSqp3aER9g+nLlB9dpmSS4jFTRUlCFcUJVZYkVFESV3NdhS5fWqe51aWaW12qBTWlWjanUg2VJTn1tZCEGACQVWam2opi1VYU65Km2aGHM+Pkx/vPAQDIUYQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgCIJsZn9iZm5mdVHsTwAAArFlENsZoskvVXSzqkPBwCAwhLFHvH/lvQZSR7BsgAAKChTCrGZ3S5pt7u/FNF4AAAoKImTzWBmj0uaN86kz0v6U0lvm8yKzOxuSXdLUlNT0ykMEQCA/GXup3dE2cwukPQLSX2Zuxol7ZF0qbvvO9FjW1pavLW19bTWCwBArjGzVe7eMt60k+4RT8Td10qaM2ol2yW1uHvn6S4TAIBCw+eIAQAI6LT3iMdy9+aolgUAQKFgjxgAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIHP36V+pWYekHREusl5SZ4TLC4ltmXnyZTsktmWmypdtyZftkKLflsXu3jDehCAhjpqZtbp7S+hxRIFtmXnyZTsktmWmypdtyZftkKZ3Wzg0DQBAQIQYAICA8iXE94QeQITYlpknX7ZDYltmqnzZlnzZDmkatyUvXiMGACBX5cseMQAAOYkQAwAQUE6F2MxuNrMNZrbZzD47znQzsy9npr9sZpeEGOdkTGJbrjezI2a2JvPzZyHGeTJm9m0zazezVyaYnkvPycm2JVeek0Vm9qSZrTezdWb2yXHmyYnnZZLbMuOfFzMrNbPnzeylzHb8+Tjz5MpzMpltmfHPyVFmFjezF83s4XGmTc9z4u458SMpLmmLpKWSiiW9JOncMfO8Q9KjkkzS5ZKeCz3uKWzL9ZIeDj3WSWzLtZIukfTKBNNz4jmZ5LbkynMyX9IlmetVkjbm8L+VyWzLjH9eMv+dKzPXiyQ9J+nyHH1OJrMtM/45GTXWT0v6p/HGO13PSS7tEV8qabO7b3X3IUkPSLpjzDx3SPqup62UVGNm86d7oJMwmW3JCe7+tKSDJ5glV56TyWxLTnD3ve6+OnO9W9J6SQvHzJYTz8skt2XGy/x37sncLMr8jH2nbK48J5PZlpxgZo2SbpX0zQlmmZbnJJdCvFDSrlG32/TGf5CTmWcmmOw4r8gc/nnUzM6bnqFFLleek8nKqefEzJolrVB6r2W0nHteTrAtUg48L5lDoGsktUv6ubvn7HMyiW2RcuA5kfQlSZ+RlJpg+rQ8J7kUYhvnvrF/hU1mnplgMuNcrfS5SS+S9BVJP872oLIkV56Tycip58TMKiU9KOlT7t41dvI4D5mxz8tJtiUnnhd3T7r7xZIaJV1qZuePmSVnnpNJbMuMf07M7DZJ7e6+6kSzjXNf5M9JLoW4TdKiUbcbJe05jXlmgpOO0927jh7+cfdHJBWZWf30DTEyufKcnFQuPSdmVqR0uO5394fGmSVnnpeTbUsuPS+S5O6HJT0l6eYxk3LmOTlqom3JkefkKkm3m9l2pV8evNHM7hszz7Q8J7kU4hcknWlmS8ysWNKdkn4yZp6fSPpQ5p1ul0s64u57p3ugk3DSbTGzeWZmmeuXKv1cHZj2kU5drjwnJ5Urz0lmjN+StN7dvzjBbDnxvExmW3LheTGzBjOryVwvk/QWSa+NmS1XnpOTbksuPCfu/jl3b3T3ZqV/Bz/h7neNmW1anpNE1AvMFncfMbNPSPqZ0u86/ra7rzOz389M/4akR5R+l9tmSX2SPhJqvCcyyW15r6T/aGYjkvol3emZt/HNJGb2PaXfIVlvZm2S/ovSb97IqedEmtS25MRzovRf+h+UtDbzOp4k/amkJinnnpfJbEsuPC/zJd1rZnGlo/TP7v5wLv7+0uS2JReek3GFeE44xSUAAAHl0qFpAADyDiEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIgTxhZs1m1j/qxBdRLrvM0t8rOzQDT1UI5DRCDOSXLZmT8UfK3fszy53R5z4GchEhBnKImf2Bmb1iZjvM7A8nMf+PzWyVma0zs7vHmd5sZq+Muv0nZvZfIx42gBPImXNNA4XOzN4j6a1KfydvvdLnX/47dx85wcM+6u4HMyfnf8HMHnT3GXXyfaDQEWIgd/yRpP/g7sOS9prZsE5+VOuPzOxdmeuLJJ2pGfYtOEChI8RADsh8J++F7r4xc3u+pE53HzrBY65X+ivqrnD3PjN7SlLpmNlGdHzMx04HkGW8RgzkhnMlzTKzpWYWk/Q/JH3lJI+ZJelQJsLLJV0+zjz7Jc0xszozK5F0W6SjBnBShBjIDSsk3S/pe5JelrTT3e85yWP+TVLCzF6W9N8krTw6wcweMbMFmcPcfyHpOUkP641f8P6ImS2IbjMAjMX3EQM5wMy+JOlZd//+CeZplvSwu5+fxXFsl9Ti7p3ZWgdQaNgjBnLDxZLWnGSepNKHr0823yk7ekIPSUWSUlEvHyhk7BEDABAQe8QAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIgRxkZtebWVvocQCYOkIMZImZfcDMWs2sx8z2mtmjZnZ16HGdiJnNN7NvZcbbbWavmdmfm1nFdDweKESEGMgCM/u0pC9J+itJcyU1Sfq6pDsiWHZiqsuYYLm1kp6VVCbpCnevkvRWSTWSzsj244FCRYiBiJnZLKW/4/fj7v6Qu/e6+7C7/6u7/6fMPN8xs78c9ZgTHmo2Mzezj5vZJkmbRt3/x2bWntkD/cjoMZjZd82sw8x2mNkXzOxk/94/Lalb0l3uvl2S3H2Xu3/S3V82s+bMOI79IWBmT5nZxyJ6PFCQCDEQvSsklUr6UcTLfaekyySdm7k9T9IsSQsl/a6kr5nZ7My0r2SmLZV0naQPSfqITuwtkh5y99P9msOpPh4oSIQYiF6dpE53H4l4uf/D3Q+6e3/m9rCkv8jsbT8iqUfS2WYWl/TvJX3O3bsze6f/S9IHJzHuvVMY31QfDxQkQgxE74Ck+iy8lrtr7HrGxL5PUqWkeknFknaMmrZD6T3nEzkgaf4UxjfVxwMFiRAD0XtW0oDSh5In0iupfNTteZNYrk9y/Z1K7y0vHnVfk6TdJ3nc45LedYLXknszlxONe6qPBwoSIQYi5u5HJP2Z0q/ZvtPMys2syMxuMbO/ycy2RtI7zKzWzOZJ+lSE609K+mdJ/93MqsxssdJvpLrvJA/9oqRqSfdmHiMzW2hmXzSzC929Q+mY32VmcTP7qI5/N/RUHw8UJEIMZIG7f1Hp+H1BUofSh5U/IenHmVn+UdJLkrZLekzS9yMewh8qvQe6VdKvJf2TpG+fZMwHJV2p9N70c2bWLekXko5I2pyZ7T9I+k9KH4Y+T9IzUT0eKFTmPtmjXQAAIGrsEQMAEBAhBgqImX0jc8rNsT/fCD02oFBxaBoAgICycs7ak6mvr/fm5uYQqwYAYNqtWrWq090bxpsWJMTNzc1qbW0NsWoAAKadme2YaBqvEQMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAEFFmIzSxuZi+a2cNRLRMAgHwX5R7xJyWtj3B5AADkvUhCbGaNkm6V9M0olgcAQKGIao/4S5I+IykV0fIAACgIUw6xmd0mqd3dV51kvrvNrNXMWjs6Oqa6WgAA8kIUe8RXSbrdzLZLekDSjWZ239iZ3P0ed29x95aGhoYIVgsAQO6bcojd/XPu3ujuzZLulPSEu9815ZEBAFAA+BwxAAABJaJcmLs/JempKJcJAEA+Y48YAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICAphxiM1tkZk+a2XozW2dmn4xiYAAAFIJEBMsYkfTH7r7azKokrTKzn7v7qxEsGwCAvDblPWJ33+vuqzPXuyWtl7RwqssFAKAQRPoasZk1S1oh6blxpt1tZq1m1trR0RHlagEAyFmRhdjMKiU9KOlT7t41drq73+PuLe7e0tDQENVqAQDIaZGE2MyKlI7w/e7+UBTLBACgEETxrmmT9C1J6939i1MfEgAAhSOKPeKrJH1Q0o1mtibz844IlgsAQN6b8seX3P3XkiyCsQAAUHA4sxYAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICAcj7E6/Yc0fee36mhkVTooQAAcMpyPsQ/WbNHn3torW74n0/pvpU7NDiSDD0kAAAmLedD/Nlbluvej16qudUl+sKPX9F1f/OU7n1muwaGCTIAYOYzd5/2lba0tHhra2uky3R3PbPlgP7P45v0/PaDmlNVot+77gx94NImlRXHI10XAACnwsxWuXvLuNPyJcSjPbvlgL78i016dusB1VcW6+5rl+quyxervDiRtXUCADCRggvxUc9vO6ivPLFJv9rUqdqKYv2Ha5bqg1csVmUJQQYATJ+CDfFRq3Yc0pd/sUm/3NihmvIifezqJfrQlc2qLi2atjEAAApXwYf4qDW7Dusrv9ikX7zWrurShD569RJ95MolmlVOkAEA2UOIx3hl9xF9+Reb9Nir+1VVktDvXNWsj161RLMrioONCQCQvwjxBF7d06WvPrlJj6zdp4riuD58ZbM+ds1S1RJkAECECPFJbNjXra88sUk/XbtXZUVxffDyxfrYNUvVUFUSemgAgDxAiCdpc3u3vvrEZv3kpT0qTsT025ct1u9du1RzqktDDw0AkMMI8Sna2tGjrz25RT9es1vxmOkDlzbp965bqvmzykIPDQCQgwjxadre2auvP7VZD63erZiZ3vfmRv3H65dpYQ1BBgBMHiGeol0H+/T1p7boh6t2SZLe+6ZG/cH1y7SotjzwyAAAuYAQR2T34X5946kt+v4Lu5Ry17svWaiP37BMi+sqQg8NADCDEeKI7TsyoG/8cov+6fmdSqZcd1y8QJ+4YZmWNlSGHhoAYAYixFnS3jWgv396q+5/bocGR1K6/qwG3XX5Yl1/9hzFYxZ6eACAGYIQZ1lH96C+++x2PfDCLnV0D2phTZk+cFmT3vumRs3lo08AUPAI8TQZTqb02Lr9um/lDj279YDMpCvPqNMdFy/UzefP40smAKBAEeIAtnb06Mcv7taP1+zRzoN9Kk7EdOPZc/TWc+fqhuVzOI0mABQQQhyQu+vFXYf1Ly/u1qOv7FN796BiJr1p8Wy95Zy5uubMBi2fV6UYrykDQN4ixDNEKuV6Zc8RPb6+XY+/ul+v7u2SJM0uL9JlS+p0xRl1unxpnc6cU0mYASCPEOIZas/hfj275YCe2XJAK7ce0O7D/ZKkiuK4LmicpYsaa3RhY40uWDhLjbPLiDMA5Kish9jMbpb0fyTFJX3T3f+/E81PiN/I3bXrYL+e23ZAL7cd0ctth7V+b7eGkilJUllRXGfMqdCyhkqdObdKZzRUqnF2mRpnl2lWWZHMiPR4UinXUDKlwZGUhpMpDY26HDp224/dP3a+kZTL3ZVyvX6po9dd7srcl74eM1MiZorFTHGT4seupy8TMUvfZ+nLeMxUnIipJB5TcSKmosxlcSKm4viYy0RMiZjxXAM5KKshNrO4pI2S3iqpTdILkt7v7q9O9BhCPDlDIylt2NetV/Yc0eb2Hm1q79GW9p5je85HVRTHtXB2mRbWlGnerFLVVhSrtqJEdRXFmevFmlVWpIqShMqL4ypJxIL8Mj8axaOhGxxJaWA4qYHhlAZGkpnrmdujL0fS1weHk+ofO09mGYOZ+16fnp42NJKa9u3MJjOpKH58uEuKXo91ydGIJ+LHrpe8YVpMJYn4scgf//hRj0sc/0dBSVH8+PXEYxylASbpRCFORLD8SyVtdvetmZU9IOkOSROGGJNTnIjpgsZZuqBx1nH39wyOaFtHr9oO9Wn34X61HerX7sP92n2oX2t3d+lQ35CSqYn/wIrHTOXFcVWWJFRWHFdRLKZ4zFQUT++hJeLpPa9EPKaYadTe4NE9wPSeoY5ddyVTrqGkazj5xj3P4WQ6wCca02SUJGIqLYqrtCimsqK4SoviKimKqzQRU015sUqLMtMTcZUVx1VSFFNpIn5caIrjr+91FsVfj83r95lKRt0+uidrSu/tmklmptjoS6XvlyR3KemuZNLTl6nX//scvT6ScqVS6ekjmf82QyNv3FsfvXc+evrgmPtGP34wc3mkf1iDw8lxpw2OJDXFp+KY9H+viaL++h8F6aMDMcVjUiLz/9uxHzPF45nLzH2J0dNPMI9GPQeZmzI7/vka/bwdf78kjXkuJ3q8MvePWl760ZnLY3+PjDfN3jCvZW6N/Xt4vGmvL89GzfPG9R4/7/jrPbqc8dY7dqzjrXf08jTeWEdt83jrPX4B4yxz1BjGG/txyxlz/3j3jTfvydZ13LzTtMMSRYgXSto16nabpMsiWC4mUFmSGDfQR6VSrq6BYR3oHdLB3iEd6BlS98Cw+oaS6h0aUd9g+rJ3cER9Q0mNJNNxGEmlMtfTv7D7hpJKuSs2NjxHf1HFTEUxk+n1Q6zp0NlxsTsauKJRMUzELBPVdFhLMgE9GtOyUdNKi8LtxeerkeTrAR8cFerBkeRxgR8cHj1fctR8Y+dJviH4R6cd6R9O/+GR+RlJpZRypS9T6ctkSkqmUsfmOfpHTPqlgdD/tVCIzp1frUc+ec20rCuKEI/32/EN/3TM7G5Jd0tSU1NTBKvFRGIxU015sWrKi3VGQ+jRYCZKxGNKxGMqz4GPs6dGhflYpDN/PLo8c2Qm/Tr90aM37jp29MalY0dzfMzt1Kh5Rz8+fcTg6O3jH3/0aIJnfs35sdvp6Uevj74yet7R82vMY3zMA9847/jrHb2c46eNmdfHX+/YsY633jdsx5iVHTfGce47bl1jxjt23OON/Y2PP/G8xy/31JdVX1ky/sKyIIoQt0laNOp2o6Q9Y2dy93sk3SOlXyOOYL0ACkAsZorJVBQPPRIgO2IRLOMFSWea2RIzK5Z0p6SfRLBcAADy3pT3iN19xMw+IelnSn986dvuvm7KIwMAoABEcWha7v6IpEeiWBYAAIUkikPTAADgNBFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAAU0pxGb2t2b2mpm9bGY/MrOaiMYFAEBBmOoe8c8lne/uF0raKOlzUx8SAACFY0ohdvfH3H0kc3OlpMapDwkAgMIR5WvEH5X06EQTzexuM2s1s9aOjo4IVwsAQO5KnGwGM3tc0rxxJn3e3f8lM8/nJY1Iun+i5bj7PZLukaSWlhY/rdECAJBnThpid3/Liaab2Ycl3SbpJncnsAAAnIKThvhEzOxmSf9Z0nXu3hfNkAAAKBxTfY34q5KqJP3czNaY2TciGBMAAAVjSnvE7r4sqoEAAFCIOLMWAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABRRJiM/sTM3Mzq49ieQAAFIoph9jMFkl6q6SdUx8OAACFJYo94v8t6TOSPIJlAQBQUKYUYjO7XdJud39pEvPebWatZtba0dExldUCAJA3EiebwcwelzRvnEmfl/Snkt42mRW5+z2S7pGklpYW9p4BANAkQuzubxnvfjO7QNISSS+ZmSQ1SlptZpe6+75IRwkAQJ46aYgn4u5rJc05etvMtktqcffOCMYFAEBB4HPEAAAEdNp7xGO5e3NUywIAoFCwRwwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQkLn79K/UrEPSjggXWS+pM8LlzVRsZ35hO/ML25lfot7Oxe7eMN6EICGOmpm1untL6HFkG9uZX9jO/MJ25pfp3E4OTQMAEBAhBgAgoHwJ8T2hBzBN2M78wnbmF7Yzv0zbdubFa8QAAOSqfNkjBgAgJ+V8iM3sZjPbYGabzeyzoceTDWb2bTNrN7NXQo8lm8xskZk9aWbrzWydmX0y9JiywcxKzex5M3sps51/HnpM2WJmcTN70cweDj2WbDKz7Wa21szWmFlr6PFki5nVmNkPzey1zL/TK0KPKWpmdnbmeTz602Vmn8rqOnP50LSZxSVtlPRWSW2SXpD0fnd/NejAImZm10rqkfRddz8/9HiyxczmS5rv7qvNrErSKknvzMPn0yRVuHuPmRVJ+rWkT7r7ysBDi5yZfVpSi6Rqd78t9Hiyxcy2S2px97z+fK2Z3SvpV+7+TTMrllTu7ocDDytrMo3ZLekyd4/y3BfHyfU94kslbXb3re4+JOkBSXcEHlPk3P1pSQdDjyPb3H2vu6/OXO+WtF7SwrCjip6n9WRuFmV+cvcv4gmYWaOkWyV9M/RYMHVmVi3pWknfkiR3H8rnCGfcJGlLNiMs5X6IF0raNep2m/LwF3chMrNmSSskPRd4KFmROWS7RlK7pJ+7ez5u55ckfUZSKvA4poNLeszMVpnZ3aEHkyVLJXVI+ofMyw3fNLOK0IPKsjslfS/bK8n1ENs49+XdnkWhMbNKSQ9K+pS7d4UeTza4e9LdL5bUKOlSM8urlxzM7DZJ7e6+KvRYpslV7n6JpFskfTzzclK+SUi6RNLfufsKSb2S8vJ9OZKUOfR+u6QfZHtduR7iNkmLRt1ulLQn0FgQgcxrpg9Kut/dHwo9nmzLHNp7StLNYUcSuask3Z557fQBSTea2X1hh5Q97r4nc9ku6UdKv2yWb9oktY06evNDpcOcr26RtNrd92d7Rbke4hcknWlmSzJ/vdwp6SeBx4TTlHkT07ckrXf3L4YeT7aYWYOZ1WSul0l6i6TXgg4qYu7+OXdvdPdmpf9dPuHudwUeVlaYWUXmzYXKHKp9m6S8+4SDu++TtMvMzs7cdZOkvHoj5Rjv1zQclpbShxpylruPmNknJP1MUlzSt919XeBhRc7Mvifpekn1ZtYm6b+4+7fCjiorrpL0QUlrM6+fStKfuvsj4YaUFfMl3Zt5R2ZM0j+7e15/vCfPzZX0o/TfkUpI+id3/7ewQ8qaP5R0f2bHZ6ukjwQeT1aYWbnSn8b5vWlZXy5/fAkAgFyX64emAQDIaYQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgqYmX3FzFab2ZtDjwUoVIQYKFCZ0zHOUfrsQXn7XcHATEeIgTxnZs1m1j/qtKGSJHfvVfp0m09J+nJm3jIzW2NmQ2ZWP+2DBQoQIQbyjKWN/be9JfO1i6Pnq5NULqlbUlKS3L0/Mx/fYgZME0IM5IHMXu96M/u6pNU6/utBJ/IFSf9T0jpJ52ZzfAAmRoiB/HG2pO+6+wp333GiGc2sWdKVkr4vab2k87I/PADjIcRA/tjh7isnOe9fSvoLT3/9GiEGAsrp7yMGcJzeycxkZhdLerekq83sa5JKJa3N4rgAnAAhBgrPX0v6LXf/hSSZ2VxJL4YdElC4ODQNFBAzu1FSxdEIS5K775dUYWa14UYGFC72iIE84O7bJZ0/ifmekPTEOPfPysKwAEwCe8RA/ktKmjX2hB7jOXpCD0lFklJZHhcASZZ+0yQAAAiBPWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYmCamdkHzKzVzHrMbK+ZPWpmV5/Gcl4zs4+Oc/8nzax1Eo9/u5k9bWbdZtZhZr80s9tPYf1TejyANEIMTCMz+7SkL0n6K0lzJTVJ+rqkO05jcfdK+tA4938wM+1E43ivpB9I+q6kxsxY/kzSb01mxVN9PIDXcYpLYJqY2SxJuyV9xN1/MME835HU5u5fyNy+XtJ97t44zryNkrZLOsPdd2TuO0fSS5IWuHvnBOswSTskfcXd/3aCef6rpGXuflfmdrOkbUqfgzo5lce7+8h4jwEKFXvEwPS5QlKppB9FsTB3b5P0pNJ7wEd9SNIjE0U442xJiyT98DRXPdXHAxiFEAPTp05SZ8R7hPcqE2Izi0n6bZ3ksHRmHJK09zTXOdXHAxiFEAPT54CkejOL8nvAH5I038wul3S9pHJJP53EOCRp/mmuc6qPBzAKIQamz7OSBiS98wTz9Cod06PmnWiB7t6n9CHiDym9Z/yAuw+dZBwbJO2S9J7THMdUHw9gFEIMTBN3P6L0O4u/ZmbvNLNyMysys1vM7G8ys62R9A4zqzWzeZI+NYlF3yvp3ysdxpMdlpan36H5aUn/r5l9xMyqzSxmZleb2T2jxnGtmTVl3mT2uageD+B4hBiYRu7+RaUj9gVJHUrvWX5C0o8zs/yj0u963i7pMUnfn8Rin5Z0RNJud39hkuP4odLx/qikPZL2S/pLSf+Smf7zzLpflrRK0sNRPh7A6/j4EgAAAbFHDABAQFG+exPADGJmPRNMusXdfzWtgwEwIQ5NAwAQUJA94vr6em9ubg6xakRow75ulRfHtai2/OQzA0ABW7VqVae7N4w3LUiIm5ub1dp60i+HwQz37q//RmXFcd3/sctDDwUAZjQz2zHRNN6shdPWUFWizu6TnTsCAHAihBinrb6yRB09g6GHAQA5jRDjtNVXluhQ35CGk6nQQwGAnEWIcdoaqkrkLh3s5fA0AJwuQozTVl9ZIknq6ObwNACcLkKM09ZQlQ5xJ68TA8BpI8Q4bQ3sEQPAlBFinLb6qmJJUmcPrxEDwOkixDht5cUJVRTH2SMGgCkgxJiS+qoSXiMGgCkgxJiShsoS9ogBYAoIMaakvpI9YgCYCkKMKWmoKlE7e8QAcNoIMaZkbnWJjvQPa2A4GXooAJCTCDGmZG51qSSpvYu9YgA4HYQYU3I0xPu6BgKPBAByEyHGlMyblQ7xfkIMAKeFEGNK5lYRYgCYCkKMKakuS6i0KEaIAeA0EWJMiZlpbnWp9vFmLQA4LYQYUza3upQ9YgA4TYQYU0aIAeD0EWJM2bzqEu3vGpC7hx4KAOQcQowpm1tdqoHhlLr6R0IPBQByTmQhNrO4mb1oZg9HtUzkhqMn9djfzeFpADhVUe4Rf1LS+giXhxxx7OxaRwgxAJyqSEJsZo2SbpX0zSiWh9wyr5qTegDA6Ypqj/hLkj4jKTXRDGZ2t5m1mllrR0dHRKvFTDCnukQSIQaA0zHlEJvZbZLa3X3VieZz93vcvcXdWxoaGqa6WswgpUVx1ZQX8cUPAHAaotgjvkrS7Wa2XdIDkm40s/siWC5yyNyqUu07wtm1AOBUTTnE7v45d29092ZJd0p6wt3vmvLIkFPm15RqX1d/6GEAQM7hc8SIxIKaMu05zKFpADhVkYbY3Z9y99uiXCZyw8KaMh3sHVL/UDL0UAAgp7BHjEgsqEl/hGnPEQ5PA8CpIMSIxIJZZZKkPYcJMQCcCkKMSCyoIcQAcDoIMSIxb1apzKTdvGELAE4JIUYkiuIxza0qZY8YAE4RIUZkFs4uI8QAcIoIMSKzoKZMuwkxAJwSQozILKgp1d7DA0qlPPRQACBnEGJEZmFNmYaSKXX2cs5pAJgsQozIvP5ZYt45DQCTRYgRGT5LDACnjhAjMgsJMQCcMkKMyFSXJVRZklDbIUIMAJNFiBEZM1Pj7DLtOtgXeigAkDMIMSK1uK5cOwkxAEwaIUakmmrTIeazxAAwOYQYkWqqLdfgSEodPXyWGAAmgxAjUotqyyWJw9MAMEmEGJFqyoR4xwFCDACTQYgRqcbZ5TJjjxgAJosQI1LFiZgWzOIjTAAwWYQYkVtUW8YeMQBMEiFG5Jpqy3mNGAAmiRAjcovrKtTZM6i+oZHQQwGAGY8QI3JHP8K06yDnnAaAkyHEiNzrH2HqDTwSAJj5CDEit5jPEgPApBFiRG52RbFmlxdpa2dP6KEAwIxHiJEVSxsqtaWDQ9MAcDKEGFmxtL5CWwkxAJwUIUZWLG2oVGfPoLoGhkMPBQBmNEKMrFjaUCFJ7BUDwEkQYmTFGcdCzBu2AOBECDGyoqm2QvGYsUcMACdBiJEVxYmYFs0u4yNMAHAShBhZs7Shkj1iADgJQoysWVpfoW2dvUqlPPRQAGDGIsTImqUNlRocSWnPEb78AQAmQoiRNcvmVEqSNrXzOjEATIQQI2vOmpsO8cZ93YFHAgAzFyFG1tSUF2tudYk27CfEADARQoysOmtulTYSYgCYECFGVi2fV6VN+3uU5J3TADAuQoysOmtulQZHUtpxgM8TA8B4CDGy6ux5VZLE4WkAmAAhRladOadKZtKGfXyECQDGQ4iRVWXFcS2uLWePGAAmQIiRdWfNreIjTAAwgSmH2MwWmdmTZrbezNaZ2SejGBjyx9nzqrSts1cDw8nQQwGAGSeKPeIRSX/s7udIulzSx83s3AiWizxx3oJqJVOuDZxhCwDeYMohdve97r46c71b0npJC6e6XOSP8xbMkiS9sudI4JEAwMwT6WvEZtYsaYWk56JcLnJb4+wyzSor0iu7CTEAjBVZiM2sUtKDkj7l7l3jTL/bzFrNrLWjoyOq1SIHmJkuWDhLawkxALxBJCE2syKlI3y/uz803jzufo+7t7h7S0NDQxSrRQ45b2G1Nuzr1tBIKvRQAGBGieJd0ybpW5LWu/sXpz4k5KMLFs7ScNL5PDEAjBHFHvFVkj4o6UYzW5P5eUcEy0UeOf/oG7Y4PA0Ax0lMdQHu/mtJFsFYkMcW15WrqjShtbuP6M7QgwGAGYQza2FamJnOW1DNHjEAjEGIMW0uaqzRq3u7OMMWAIxCiDFtLlk8W8NJ1zpO7AEAxxBiTJtLmmZLklbtOBR4JAAwcxBiTJuGqhI11ZYTYgAYhRBjWr1p8Wyt3nlY7h56KAAwIxBiTKtLmmrU0T2otkP9oYcCADMCIca0umQxrxMDwGiEGNPq7LlVqiiOE2IAyCDEmFaJeEwrmmbr+W0HQw8FAGYEQoxpd8UZddqwv1udPYOhhwIAwRFiTLsrz6iTJK3ceiDwSAAgPEKMaXfBwlmqLEnomS2EGAAIMaZdIh7TpUtq9SwhBgBCjDCuPKNO2zp7tfcInycGUNgIMYK4IvM68TOb2SsGUNgIMYI4Z161aiuK9atNHaGHAgBBEWIEEYuZrj+rQU9t7FAyxXmnARQuQoxgblg+R4f7hvXiTs6yBaBwEWIEc+1ZDYrHTE+81h56KAAQDCFGMLPKivSmxbMJMYCCRogR1I3L5+i1fd3ac5iPMQEoTIQYQb3lnDmSpMfW7Qs8EgAIgxAjqGVzqnTW3Eo9spYQAyhMhBjB3XrBAr2w46D2HRkIPRQAmHaEGMHdeuE8uUuPvrI39FAAYNoRYgS3bE6Vls+r0k9fJsQACg8hxoxw6wXz1brjkNoO9YUeCgBMK0KMGeGdKxZKkh5ctTvwSABgehFizAiLast11bI6/WDVLqU49zSAAkKIMWO8r2WR2g71a+VWvhoRQOEgxJgx3n7ePFWVJvT91l2hhwIA04YQY8YoLYrrXSsW6tG1+9TRPRh6OAAwLQgxZpQPX9msoWRK963cEXooADAtCDFmlDMaKnXj8jm6b+UODQwnQw8HALKOEGPG+d2rl+hA75B+8tKe0EMBgKwjxJhxrjyjTsvnVekbv9yiJB9lApDnCDFmHDPTJ286U1s7evWTlzjBB4D8RogxI739vHk6Z361/s/jmzSSTIUeDgBkDSHGjBSLmf6ft5yp7Qf69NBq9ooB5C9CjBnrrefO1cWLavS3j21Q98Bw6OEAQFYQYsxYZqY/v/08dXQP6itPbA49HADICkKMGe2iRTV6X0ujvv3rbdrc3h16OAAQOUKMGe8zNy9XZWlCf/KDl3njFoC8Q4gx49VXlui/3XG+1uw6rL9/emvo4QBApAgxcsJvXbRAt144X196fKNe2nU49HAAIDKEGDnjL+84X3OqSvX7961SZw/fzgQgPxBi5IzZFcX6+w++SQd7h/Tx+1drcIQvhQCQ+wgxcsr5C2fpr99zoZ7bdlCf/v5LnIsaQM6LJMRmdrOZbTCzzWb22SiWCUzknSsW6gu3nqOfrt2rP31orVLEGEAOS0x1AWYWl/Q1SW+V1CbpBTP7ibu/OtVlAxP52DVL1dU/rC8/sVk9QyP64vsuUkkiHnpYAHDKphxiSZdK2uzuWyXJzB6QdIckQoys+vTbzlZlaUJ/9chrOtAzqK9+4BLVV5aEHhZmiFTKNZJyJVOukVRKyZRrOHn87WPTk+n7Ui6l3OXu6eup9OWx2+5KHp2eSt8eOz39eCmZev16atT04+c9+tj0fJLkktwl1+v3Sen50ve/Pl3HrqevHJ2WXs748+vYfa9P1zjrdPdxp41eh47eN+ag1OuPfP3xx0/XCaePt4yT3JSPWcjJ1nGyxzfOLtdnb1n+xoFlQRQhXihp16jbbZIuGzuTmd0t6W5JampqimC1gHT3tWeooapEn31wrW798q/05TtX6LKldaGHVXDcXYMjKQ0MJ9U3lFT/cFL9Yy4HMteHkykNjqQ0lExpaGTUT3LM5UTXk6nj4nksqMnXwzqcSo37yz3fmUmm9OlhpaPXJVN6gh03n407v0bdl7n79duZx0pjp2XWMWYsx91+w1jH3jP+9px4mVNb5xtGMOqOoZHpO3lQFCEe77/mG/4JuPs9ku6RpJaWlgL8J4JsedeKRp09t1p/cP8q3fl/V+quyxbrP918tqpLi0IPbUYbSabUPTCinsGRY5c9g8PqHhh1O3PZNTB87Hrv4Mix2I4O7+mGLxEzFSdiKk7EVBSPqTgeU0nmdnEifbs4EVNlaSI9PRFTUcwUj8WUiJkScVPi6O24KR5L306Muf36ZWzUY9LzxUdNN5NiZpmf9C/vmKW/Eezo7Xhm+rF5Yxozf+YxZsced9z8o+Y9et/rkXs9ktLxkRsdzuPmn0TUMHNFEeI2SYtG3W6UtCeC5QKTdu6Caj38R9fofz22Qfc+s12PvrJPf3D9GfrAZU0qLcr/1457B0d0oGdIh/rSP0f6h3Wod0iH+obT1/uGdLhvWIf70vcd7htS18DISZcbM6myJKGq0iJVlSZUWZJQTXmxFtTEVVYUV2lxXOVFcZUVx1ValL6vvPj12+XFmfky10uK4sfCWpIJbzxGRFDYbOxx8VNegFlC0kZJN0naLekFSR9w93UTPaalpcVbW1untF5gIi+3HdZfPbJeK7ce1JyqEn3gsia9r2WRFtSUhR7aKekfSqqzZ1AdPYPq6B5MXz/ucujY7b6hiT9TXV2a0OyKYtWUFammvFg15UWanbmszgQ2HdlMbEsTqipJX5YVxdnbAiJgZqvcvWXcaVMNcWYF75D0JUlxSd929/9+ovkJMabDs1sO6O9+uUW/2tQhk3TZkjrddM4c3bh8jpbUVwQJzOBI8lhAx4/r65HtGRx/j7W2olj1lcVqqCpRfWWJGipLVF9VorqKYtVWFB8X21llRexxAjNA1kN8qggxptOug33659Zd+tm6fdq4v0dSOmYXL6rR+QuqtbiuQs315VpYU66a8qJJH8p2dw0Mp3S4P30oOH3od1hd/cPq7B291zpwLLwTHQ6eVVZ0fFxHXTYcvawqUW1FsYrinIcHyDWEGMjYeaBPT2/q0Jpdh7Vm12Ft7ejR2POBlBbFVFNWrOJE7NgbeuIx01AypcHh9Dt+B0eSGhxOv4N3IhXF8WMBHb33elxsq0pUX1nMZ6CBPHeiEEfxZi0gZzTVleuuusW66/LFktIfUWg71KcdB/q050j/sTc0He4b1sixz5mmNJJ0FSViKk3EVVKUfqNRSSKuWWVFqikvSl+WFam6LH29rrJY5cX88wJwcvymQEErTsS0tKFSSxsqQw8FQIHixSYAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACCgKYXYzP7WzF4zs5fN7EdmVhPRuAAAKAhT3SP+uaTz3f1CSRslfW7qQwIAoHBMKcTu/pi7j2RurpTUOPUhAQBQOKJ8jfijkh6NcHkAAOS9xMlmMLPHJc0bZ9Ln3f1fMvN8XtKIpPtPsJy7Jd0tSU1NTac1WAAA8s1JQ+zubznRdDP7sKTbJN3k7n6C5dwj6R5JamlpmXA+AAAKyUlDfCJmdrOk/yzpOnfvi2ZIAAAUjqm+RvxVSVWSfm5ma8zsGxGMCQCAgjGlPWJ3XxbVQAAAKEScWQsAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAiIEAMAEBAhBgAgIEIMAEBAhBgAgIAiCbGZ/YmZuZnVR7E8AAAKxZRDbGaLJL1V0s6pDwcAgMISxR7x/5b0GUkewbIAACgoUwqxmd0uabe7vzSJee82s1Yza+3o6JjKagEAyBuJk81gZo9LmjfOpM9L+lNJb5vMitz9Hkn3SFJLSwt7zwAAaBIhdve3jHe/mV0gaYmkl8xMkholrTazS919X6SjBAAgT500xBNx97WS5hy9bWbbJbW4e2cE4wIAoCDwOWIAAAI67T3isdy9OaplAQBQKNgjBgAgIEIMAEBAhBgAgIAIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGACAgQgwAQECEGACAgAgxAAABEWIAAAIixAAABESIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAERIgBAAjI3H36V2rWIWlHhIusl9QZ4fJmKrYzv7Cd+YXtzC9Rb+did28Yb0KQEEfNzFrdvSX0OLKN7cwvbGd+YTvzy3RuJ4emAQAIiBADABBQvoT4ntADmCZsZ35hO/ML25lfpm078+I1YgAAclW+7BEDAJCTCDEAAAHlfIjN7GYz22Bmm83ss6HHkw1m9m0zazezV0KPJZvMbJGZPWlm681snZl9MvSYssHMSs3seTN7KbOdfx56TNliZnEze9HMHg49lmwys+1mttbM1phZa+jxZIuZ1ZjZD83stcy/0ytCjylqZnZ25nk8+tNlZp/K6jpz+TViM4tL2ijprZLaJL0g6f3u/mrQgUXMzK6V1CPpu+5+fujxZIuZzZc0391Xm1mVpFWS3pmHz6dJqnD3HjMrkvRrSZ9095WBhxY5M/u0pBZJ1e5+W+jxZIuZbZfU4u55faILM7tX0q/c/ZtmViyp3N0PBx5W1mQas1vSZe4e5UmojpPre8SXStrs7lvdfUjSA5LuCDymyLn705IOhh5Htrn7XndfnbneLWm9pIVhRxU9T+vJ3CzK/OTuX8QTMLNGSbdK+mbosWDqzKxa0rWSviVJ7j6UzxHOuEnSlmxGWMr9EC+UtGvU7Tbl4S/uQmRmzZJWSHou8FCyInPIdo2kdkk/d/d83M4vSfqMpFTgcUwHl/SYma0ys7tDDyZLlkrqkPQPmZcbvmlmFaEHlWV3SvpetleS6yG2ce7Luz2LQmNmlZIelPQpd+8KPZ5scPeku18sqVHSpWaWVy85mNltktrdfVXosUyTq9z9Ekm3SPp45uWkfJOQdImkv3P3FZJ6JeXl+3IkKXPo/XZJP8j2unI9xG2SFo263ShpT6CxIAKZ10wflHS/uz8UejzZljm095Skm8OOJHJXSbo989rpA5JuNLP7wg4pe9x9T+ayXdKPlH7ZLN+0SWobdfTmh0qHOV/dImm1u+/P9opyPcQvSDrTzJZk/nq5U9JPAo8JpynzJqZvSVrv7l8MPZ5sMbMGM6vJXC+T9BZJrwUdVMTc/XPu3ujuzUr/u3zC3e8KPKysMLOKzJsLlTlU+zZJefcJB3ffJ2mXmZ2duesmSXn1Rsox3q9pOCwtpQ815Cx3HzGzT0j6maS4pG+7+7rAw4qcmX1P0vWS6s2sTdJ/cfdvhR1VVlwl6YOS1mZeP5WkP3X3R8INKSvmS7o3847MmKR/dve8/nhPnpsr6UfpvyOVkPRP7v5vYYeUNX8o6f7Mjs9WSR8JPJ6sMLNypT+N83vTsr5c/vgSAAC5LtcPTQMAkNMIMQAAARFiAAACIsQAAAREiAEACIgQAwAQECEGCpiZfcXMVpvZm0OPBShUhBgoUJmzQM1R+qQFefsVhcBMR4iBPGdmzWbWP+psZZIkd+9V+ixfT0n6cmbessyXoQ+ZWf20DxYoQIQYyDOWNvbf9pbMtz2Nnq9OUrmkbklJSXL3/sx8fHkKME0IMZAHMnu9683s65JW6/hvJZvIFyT9T0nrJJ2bzfEBmBghBvLH2ZK+6+4r3H3HiWY0s2ZJV0r6vqT1ks7L/vAAjIcQA/ljh7uvnOS8fynpLzz9rS+EGAgop78GEcBxeiczk5ldLOndkq42s69JKpW0NovjAnAChBgoPH8t6bfc/ReSZGZzJb0YdkhA4eLQNFBAzOxGSRVHIyxJ7r5fUoWZ1YYbGVC42CMG8oC7b5d0/iTme0LSE+PcPysLwwIwCewRA/kvKWnW2BN6jOfoCT0kFUlKZXlcACRZ+k2TAAAgBPaIAQAIiBADABAQIQYAICBCDABAQIQYAICACDEAAAERYgAAAiLEAAAE9P8DwBoSJ+DLvO4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x1296 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "job.plot_final_potential()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "adult-democracy",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The job PotentialTest was saved and received the ID: 820\n",
      "The job strain_0_9 was saved and received the ID: 821\n",
      "The job strain_0_92 was saved and received the ID: 822\n",
      "The job strain_0_94 was saved and received the ID: 823\n",
      "The job strain_0_96 was saved and received the ID: 824\n",
      "The job strain_0_98 was saved and received the ID: 825\n",
      "The job strain_1_0 was saved and received the ID: 826\n",
      "The job strain_1_02 was saved and received the ID: 827\n",
      "The job strain_1_04 was saved and received the ID: 828\n",
      "The job strain_1_06 was saved and received the ID: 829\n",
      "The job strain_1_08 was saved and received the ID: 830\n",
      "The job strain_1_1 was saved and received the ID: 831\n",
      "job_id:  821 finished\n",
      "job_id:  822 finished\n",
      "job_id:  823 finished\n",
      "job_id:  824 finished\n",
      "job_id:  825 finished\n",
      "job_id:  826 finished\n",
      "job_id:  827 finished\n",
      "job_id:  828 finished\n",
      "job_id:  829 finished\n",
      "job_id:  830 finished\n",
      "job_id:  831 finished\n"
     ]
    }
   ],
   "source": [
    "lmp = pr.create_job(\"Lammps\", \"template\", delete_existing_job=True)\n",
    "lmp.structure = pr.create_ase_bulk(\"Cu\", cubic=True)\n",
    "lmp.potential = job.lammps_potential\n",
    "murn = lmp.create_job(\"Murnaghan\", \"PotentialTest\")\n",
    "murn.run(delete_existing_job=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "vocal-heather",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEbCAYAAAAWFMmuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABBXklEQVR4nO3dd5xU1fnH8c+XpQmCKEVR1AUDERFEXSxRVJCosSEoihpjV4z6U1HsJpJEE4kKsUuwYOwaWzTGKIK9gIpKERsoCCpC6EXK8/vj3HFmh5nZ2d2Zndnd5/163Re333Mvs/e5555zz5GZ4ZxzzuVSg0InwDnnXN3jwcU551zOeXBxzjmXcx5cnHPO5ZwHF+eccznnwcU551zOeXBxVSLpakn3Fzodzrni5MGlCEiaJelHSW2S5k+WZJJKC5Q0lwVJbSS9IWmBpEWS3pK0V4b1t5L0tKSFkuZIGpK0vK+k9yUtkfSlpDMSlg2WNEPSYknfSxorqWWKY3SWtCrxAUDSDpImSfpfNLwkaYeE5ZJ0XXQeCySNkKRoWTtJD0maGx37DUm7pzm/e6Lf7c8S5m0m6RFJP0TDA4npljQ6Oq/1kk5K2p8k/UnSN9GxJ0jqlrB8QnSuy6JhRsKy0igtyxKGq7I554Ttx0taIekTSf0Slu0XpTdx3yemuibZqOBYlycdZ2V07DaZ9llQZuZDgQdgFjADODdhXvdongGlVdhnwzyn+Wrg/kJfu8peB0BAg0psX+H6QFPg54SHNQFHAAvT/R8A44FRQCNgp2jdPtGyRsBi4MxoX72AZcBO0fKtgTbR+MbAA8BNKY7xX+C1xP8joBVQGu23BPg/4KOE5WdGv7kOwFbANGBItKwTMBRoH217BvADsHHScfcGXo1+tz9LmH9blKaWwCbAS8CNCcvPBvYHJgEnJe3zaGBulIYS4M/A+wnLJwCnpbnWpVFa0v1fpD3naPlbwI3ARsCRwCKgbbRsP2BODn+raY+V5u/v5UL/fWUaPOdSPP4B/CZh+kTgvsQVoie00xKmT5L0esK0STpb0mfAZ9GT1RxJF0ZPufMknZyw/iGSPoiekGdLujrpeL+R9FX0RHeVQg6rX8IqjSXdJ2mppKmSyhK2vVTSF9GyaZIGJKdb0vXRE/RMSb/K9kJJ2lLSPyXNj7b9v4RlV0t6XNL9kpYAJ0XX7RpJbwArgE6SfiFpYvQkPFHSL5Kuc7n1M6XHzFaZ2QwzW0+4ca8DNgU2S5H2jQk3pWvMbI2ZfQg8DpwSrbIZ4Qb8DwsmAtOBHaJjzTazHxJ2uQ74WcI0kgYTbkzjktK5yMxmWbg7KcW2JwI3mNkcM/sGuAE4Kdr2SzO70czmmdk6MxsNNCYE1dhxGwI3A+ekuEwdgafMbImZLQaeBH7KfZjZrWY2DliVZtvXozSsA+6PXY8cSHvOkroAuwC/N7OVZvZP4GPCjb9CkjaRdFf0d/dNlPsqSbNu1seKclYnAGMre7I1yYNL8XgbaCmpa/QDPIbwR1RZRwC7E//j24LwpLgVcCpwq6RNo2XLCQGtFXAIcJakIyC8QiE8bR5PeFqN7SPR4cDD0fbPALckLPsC6B1tNxy4X1L7hOW7E54Y2wAjgLtiryOiwPRsqpOT1AD4F/BhlJ79gfMlHZiwWn/CDbsV4ckewh/jGUALYCnwHHAT0JrwtPicpNYJ+0hc/ytJz0q6NFWaEtL2EeHm+Awwxsy+T7Va0r+x8R0BzOw74CHgZEklkvYEtgUSHyL2lrQ4Oo8jCbmg2LKWwB+ACzOkc1GUzpuBaxMWdSNc15gPSQgASfvoSQgunyfMvgB41cw+SrHJrcChkjaNfn9HAs+nS2OSh4GfSeoiqREhIPwnaZ0/K7xue0PSfin28VX0oHVP0qukTOfcDfjSzJamWQ7QTtJ30UPOSEnNE5aNBdYSAvjOwAHAaaSWzbFiegObA/9Ms6/iUOiskw8/vRbrB1xJyPIfBLwINCThtRhJ2X/CE9brCdMG9E2Y3g9YScIrAeB7YI806RgFjIzGfwc8lLCsGfAj0C+avhp4KWH5DsDKDOc4GeifkO7Pk/ZtwBZZXKvdga+T5l0G3JOQrleTlk8A/pAwfQLwbtI6bxG9jklev5L/l02BY4ETM6zzOuHG3pTwtLoQmJGw/DDgO8KNaS1wepr9bBWdb5eEeX8DLkm4FilfXQLNgd8ChyTMWwdsnzDdOfp/UdK2LQlP1ZclzNuaEGg2SfgtJr4W25LwKmx9NLwINE5zbU5Kmtc4Oi+LrsdMoGPSb6IF0IQQeJYC20XLNgbKCH9LmxMeOl7I5pyj38nbSWm5Brg3Gt+C8LtvQMhdvQrcGS3bHFgNbJSw7bHA+DT/HxmPlTT/rlTzi23wnEtx+QdwHOHme1/mVdOanTS9wMzWJkyvIPzBIWn3qABxfvQkPISQk4BwM/hpX2a2AliQtO9vk/bbNHo1EnulNlmhgHsR4cm8Tapto30TS1cFtgW2jO032vflhD/mmORrkDxvS+CrpOVfUT5nlmofFbLwiuwh4FJJO6VZ7XjCzWg2cDshdzUHQNL2wCOEHGVjwpPrxZIOSXGsbwhP8A9H2/YkPKSMzCKdy4E7gPsktYtmLyMEjpiWwDKL7mjRMTYi5BzfNrM/J6w7ihCQF6c55GPAp4Qg0JKQs802Z/57QtnT1oSAPBx4WVKz6FzeMbOlZrbazMYCbwAHR8uWmdkkM1trIVd4DnCA4pUJMp1z8rLY8qXRvr81s2lmtt7MZgIXA0dF621LKD+bl/A7vRNoBxC9Ro4Vzveu6Fgx0fUfRJG/EgN/LVZUzOwrwlPZwcATKVZZTnjKj9ki1W4qccgHCa9wtjazTQg3m9jrmnmEQk7gpx916w32kIKkbYG/E/6QW5tZK2AK5V8FVdVsYKaZtUoYWpjZwQnrpLoGifPmEv74E20DfFPBPiqjEWnKaszsKzM71MzamtnuhOv6brR4R0Iu5oXopjWD8AovXZlUQ2C7aHw/QgH215K+BS4CjpT0fpptGxB+T7GgOpVQwSBmp2geAJKaAE8RrtOZSfvaH/irpG+jYwO8Jem4hH3daWbLzWwZ4bd2MNnZCXjEQrnIWjO7l1Cmla7cJZbzSLeMhOWZznkqoXyuRZrlmY47m5BzaZPwO21pZt0AzKybmW0cDa9V4lgDCTndCWnSUDQ8uBSfUwmvtpanWDYZGCipmUI1z1OreawWwEIzWyVpN0KuKeZx4LCo4Lsx4Wkx2+DQnPCHNh9AoRLBjtVMa8y7wBJJl0jaKCqX2FFSr0rs499AF0nHSWoo6RjCjSplOU9FJO0RlYM0jtJ0CSEn9U6a9btKahGt/2vCu/gbo8UfAJ0VqiNL0nbAoUTlApKOl7RNtGxbwquTWMH9aEKg6RkNdxAC04HRtr+UtHN0zVpGx/wfocIAhNzyUIWq0lsSym3ujbZtRPhNrAR+Y6HyQqIuhJth7NgQXu89GY1PBE6Lrs9GhPKsn8o6omvRlPAbaySpaVS+Ftt2kKTNJTWQdAIheH8uqZWkA6P1G0o6HtgHeCHa7+6Sfh5t15pQzjYhIYeV9pzN7FPC39zvo/0PAHoQlXUoVJiJ/V9sDfwFeDradh6hdtwNklpGx99O0r6kUNGxEpwI3JeYmyxahX4v50O8zCXF/OQylzaEH+xSQtb/ajYsc0l8z70fSVUlE49FyMJ/Fe3vWUKBfGLV1ZOArwmvw64iPLH2jpZdnbRuKQlVPgk3vYWE6qo3Aq8QlReRVFaUnHbCa67nM1yvLQmF3t8Sbo5vU74s6P6k9SeQVFWVUGX2PUK13/eAvStY/3ng8jTp2Zdwo1wanfMrwD4Jy48HpiZMn08IvMsJZQxlSfs7mpDTW0p4XXYdUXXo6LrOibadQwgordOkK/n/aBDwCeEVzHxCkO2RsFyEyhULo2EEUXlLdI5GeP25LGHonebYyb/FjoTXaQuiff8H6Jx0zS1p2C9a1pRQIWAesAR4HzgoWtaWEHyWEmrIvQ38MmG/xxLeBiyPtr+PhLK9TOec8LueQAiqM0j4OyVUzf4muiazCeVoLRKWb0J47TmH8Dv7ABic4Xed9ljR8q2IKggU+p6VzRD74TiXkUIV2kWEG8LMAifHOVfk/LWYS0vSYdEruObA9YQaQrMKmyrnXG3gwcVl0p9Q+D2XUEVzsHlW1zmXBX8t5pxzLuc85+Kccy7nGhY6AcWgTZs2VlpausH8b7+F5s2hRYsNt6mspUth+XLYItWXKc45Vwu99957P5hZ21TLPLgApaWlTJo0aYP548fD0UfDHXdAnz5V339sP48+Wr39OOdcMZGU3NLFT/y1WAZ9+oSAcPTRIUBUhQcW51x95MGlAtUJMB5YnHP1lQeXLFQlwHhgcc7VZ17mkqXEAFNRwPDA4lzNWLNmDXPmzGHVqlR9jLlcadq0KR06dKBRo0ZZb+PBpRKyCTAeWJyrOXPmzKFFixaUlpYi5aLRbZfMzFiwYAFz5syhY8eOWW/nr8UqKdMrMg8sztWsVatW0bp1aw8seSSJ1q1bVzp36MGlChIDzL//DV9/7YHFuULxwJJ/VbnG/lqsirp0gd694dBDoXNnWLjQA4tzzsV4zqWKli+HJ58EM/j0Uzj2WA8szrmK7bfffik/2q4Jv/jFLypcZ+ONs+ltvGIeXKqoSxcoK4tPjxlT9Q8tnXOuJrz55ps1diwPLlU0fnzIscSUlFTvS37nXO00a9Ystt9+e0488UR69OjBUUcdxYoVKxg3bhw777wz3bt355RTTmH16tXltrvrrru44IILfpr++9//ztChQ5k1axZdu3bl9NNPp1u3bhxwwAGsXLkSgMmTJ7PHHnvQo0cPBgwYwP/+9z8g5IYuuOAC9tlnH7p27crEiRMZOHAgnTt35sorr/zpGLFcybJly9h///3ZZZdd6N69O08//XTuL0whur8kdLc6FVhPUhev0fJtCF2oXlTBfi4idIfaJqGb0JWEvqgnA3dkk55dd93VKuPll83atDF78UWzTp3Mwssxs6FDw/yXX67U7pxzVTRt2rSfxmN/h/kYMpk5c6YB9vrrr5uZ2cknn2x//OMfrUOHDjZjxgwzMzvhhBNs5MiRZma277772sSJE23ZsmXWqVMn+/HHH83MbM8997SPPvrIZs6caSUlJfbBBx+YmdmgQYPsH//4h5mZde/e3SZMmGBmZldddZWdd955P+3z4osvNjOzUaNGWfv27W3u3Lm2atUq22qrreyHH34wM7PmzZubmdmaNWts8eLFZmY2f/5822677Wz9+vXl1sl0rePXnEmW5r5aqJzLFGAg8Gqa5SMJfZanJWlr4JeEPt4TfWFmPaNhSLVTmiSxVli/fnD22fFlL70EjzziORjn6putt96avfbaC4Bf//rXjBs3jo4dO9KlSxcATjzxRF59tfztrnnz5vTt25dnn32WTz75hDVr1tC9e3cAOnbsSM+ePQHYddddmTVrFosXL2bRokXsu+++Kfd5+OGHA9C9e3e6detG+/btadKkCZ06dWL27Nnljm1mXH755fTo0YN+/frxzTff8N133+X0mhQkuJjZdDObkWqZpCOALwk5m0xGAhcTci41IlV141NOgWbNwvhHH0HDhtVv7NI5V7tUtTr0aaedxr333ss999zDySef/NP8Jk2a/DReUlLC2rVrK9xXbJsGDRqU275BgwYbbP/AAw8wf/583nvvPSZPnszmm2+e81YOiqrMJeqr/RJgeAXrHQ58Y2YfpljcUdIHkl6R1DvDPs6QNEnSpPnz51eYtnTfsbRqBSecEJ+++ebctKbsnKucfL4Yq8jXX3/NW2+9BcBDDz1Ev379mDVrFp9//jkA//jHP37KcSTafffdmT17Ng8++CDHHntsxmNssskmbLrpprz22msZ95mNxYsX065dOxo1asT48eP56qu0LedXWd6Ci6SXJE1JMfTPsNlwYKSZLcuw32bAFcDvUiyeB2xjZjsDQ4EHJbVMtR8zG21mZWZW1rZtyr5uflLRB5KJr8aefBJmz/YA41x90rVrV8aOHUuPHj1YuHAhF1xwAffccw+DBg2ie/fuNGjQgCFDUr+lP/roo9lrr73YdNNNKzzO2LFjGTZsGD169GDy5Mn87nepboMVO/7445k0aRJlZWU88MADbL/99lXaT0bpCmNqYgAmkFCgD7wGzIqGRcBC4JykbboD3yest5ZQ7rJFRftPN2Qq0I8V3ldUSN+nT/w555JLKr+9c67yUhUy17SZM2dat27dqrz9IYccYi+99FIOU5QftaVAPyUz621mpWZWCowCrjWzW5LW+djM2iWsNwfYxcy+ldRWUgmApE5AZ0L5TZVUpkmX886Lj48eDStWhHHPwTjnUlm0aBFdunRho402Yv/99y90cnKuIMFF0gBJc4A9geckvZDFNmMklVWw2j7AR5I+BB4HhpjZwqqmc+LE7Jt0OfRQiDUY+r//wf33x5fFAszEiVVNiXOuWJWWljJlypRKb9eqVSs+/fRTHnvssTykqvBk2ZRW1XFlZWWWi+YYRo6EoUPDeLdu8PHH4G3qOZc/06dPZ/vtt/fGK/PMzPjkk0/o2rVrufmS3jOzlA/9RfVarLY75RSINcszdSqMG1fY9DhX1zVt2pQFCxbgD8n5Yxb6c2natGmltvNWkXNok03gpJPglqiUaNSo8KGlcy4/OnTowJw5c8jmcwJXdbGeKCvDX4uRu9diAJ99Fhq1jJkxo/y0c87VFf5arAZ17gyHHBKf/tvfCpcW55wrFA8ueRAr1Ae4997QkZhzztUnHlzyoE8f6NEjjK9YEb57cc65+sSDSx5I5XMvN98MP/5YuPQ451xN8+CSJ4MHw+abh/G5c6GOfiflnHMpeXDJkyZN4Jxz4tM33phd66rOOVcXeHDJoyFDIPbd0fvvwyuvFDY9zjlXUzy45FGbNvCb38Snb7ihcGlxzrma5MElzxIL9p99FqZPL1xanHOupnhwybOf/xyirq2BUPbinHN1nQeXGnDRRfHx++6Db78tXFqcc64meHCpAXvvDbvvHsZ//BFuvbWw6XHOuXzz4FIDpPK5l9tug+XLC5ce55zLNw8uNWTAgHhPlQsXwt13FzY9zjmXTx5cakhJSfmaYzfcAGvXFi49zjmXTx5catApp4RvXwC++goefbSw6XHOuXzx4FKDmjWDc8+NT48Y4U3COOfqpoIEF0mDJE2VtF7SBr2YSdpG0jJJF6XZ/mpJ30iaHA0HJyy7TNLnkmZIOjCf51EVZ58dggzAhx/CCy8UNj3OOZcPhcq5TAEGAq+mWT4SeL6CfYw0s57R8G8ASTsAg4FuwEHAbZJKcpTmnGjdGk4/PT593XWFS4tzzuVLQYKLmU03sxmplkk6AvgSmFqFXfcHHjaz1WY2E/gc2K3KCc2TCy4IBfwAEybAO+8UNDnOOZdzRVXmIqk5cAkwPIvVz5H0kaS7JW0azdsKmJ2wzpxoXqpjnSFpkqRJ8+fPr1a6K2vbbeHYY+PTf/lLjR7eOefyLm/BRdJLkqakGPpn2Gw44XXXsgp2fzuwHdATmAfE2htWinVTFpmb2WgzKzOzsrZt21ZwuNy75JL4+FNPwbRpNZ4E55zLm4b52rGZ9avCZrsDR0kaAbQC1ktaZWa3JO37u9i4pL8Dz0aTc4CtE1btAMytQjrybscdQ4OWzzwTpq+7DsaOLWyanHMuV4rqtZiZ9TazUjMrBUYB1yYHFgBJ7RMmBxAqCAA8AwyW1ERSR6Az8G5+U111l10WH3/gAZg1q2BJcc65nCpUVeQBkuYAewLPSaqwQq6kMQnVlkdI+ljSR0Af4AIAM5sKPApMA/4DnG1m6/JyEjmwxx7Qp08YX7cOrr++sOlxzrlckflXfJSVldmkSZMKcuwXX4QDDgjjTZuG3MvmmxckKc45VymS3jOzDb5VhCJ7LVYf9esHZdF/zapV3pmYc65u8OBSYBJcfnl8+rbbYMGCwqXHOedywYNLEejfH7p1C+PLlsFNNxU2Pc45V10eXIpAgwZwxRXx6b/9DRYvLlx6nHOuujy4FImjj4bOncP44sXeFbJzrnbz4FIkSkrKl73ceKN3heycq708uBSR44+H0tIwvmAB3H57QZPjnHNV5sGliDRqBJdeGp/+619hxYrCpcc556rKg0uROekk2DpqHe377+GOOwqaHOecqxIPLkWmSZPybY6NGOG5F+dc7ePBpQidcgp06BDGv/sORo8ubHqcc66yPLgUoSZNype9XHcdrFxZuPQ451xleXApUqeeCltFfWh++y3ceWdh0+Occ5XhwaVINW1aPvfyl7942Ytzrvbw4FLETjutfNnLbbcVNj3OOZctDy5FaMQIGD8+5F4S2xy77rrQsGVljB8f9uecczXJg0sR6tUrtDU2fnyoObbNNmH+Dz/ALRt0+pze+PFhP7165SedzjmXjgeXItSnDzz6aAgMb7wBV10VX/bXv8KSJRXvIxZYHn003pWyc87VFA8uRSoxwGy7LXTqFOYvXAijRmXe1gOLc67QPLgUsViAOe44GDQoPv+GG0KQScUDi3OuGBQkuEgaJGmqpPWSylIs30bSMkkXpdn+aknfSJocDQdH80slrUyYX+tb5ooFmLvuirc5tmRJ6kJ6DyzOuWJRqJzLFGAg8Gqa5SOB5yvYx0gz6xkN/06Y/0XC/CG5SGyhxQLMokXxeTfdBPPmxac9sDjniklBgouZTTezGamWSToC+BKYWqOJKnJ9+sCTT0LDhmF65Uq45pow7oHFOVdsiqrMRVJz4BJgeBarnyPpI0l3S9o0YX5HSR9IekVS7wzHOkPSJEmT5s+fX92k14j994c//CE+PXo0PPSQBxbnXPHJW3CR9JKkKSmG/hk2G0543VXRp4K3A9sBPYF5wA3R/HnANma2MzAUeFBSy1Q7MLPRZlZmZmVt27atzKkV1KWXQrduYXzNGjj5ZA8szrni0zBfOzazflXYbHfgKEkjgFbAekmrzKzcp4Nm9l1sXNLfgWej+auB1dH4e5K+ALoAk6p0EkVICs3A7LtvmF69Gtq1K2yanHMuWVG9FjOz3mZWamalwCjg2uTAAiCpfcLkAEIFASS1lVQSjXcCOhPKb+qUdetCl8gxZ5xRuLQ451wqhaqKPEDSHGBP4DlJL2SxzZiEassjJH0s6SOgD3BBNH8f4CNJHwKPA0PMLM0XIbVTrPD+1lvj8958E26+uXBpcs65ZDKzQqeh4MrKymzSpOJ/c5ZcK+z44+HBB8Oyhg3hhRegb9/CptE5V39Ies/MNvhWEYrstZhLL1V14z/+Mf56bO1aGDAgrOecc4XmwaUWSPcdS6dOcOaZ8enNNou3puycc4WUMbhI2lPSrdH3JPMlfS3p35LOlrRJTSWyPqvoA8mrroKNNw7js2aFqskeYJxzhZY2uEh6HjgNeAE4CGgP7ABcCTQFnpZ0eE0ksr7K5sv7du1g2LD49AMPwNixHmCcc4WVtkBfUhsz+yHjxlmsUxsUY4F+ZZp0Wb4cfvYz+PbbMP2nP8EvfuFf7jvn8quqBfpXS/pFph3XhcBSjCrbVljz5uWbhbnuuvAVf6w/GM/BOOdqWqbg8hlwg6RZkq6T1LOG0lTvTZxY+RzHySfDDjuE8aVLYfjweGvKEyfmJ53OOZdOhd+5SNoWGBwNTYGHgIfN7NP8J69mFONrsap49lk47LAwXlICU6bA9tsXNk3OubqrWt+5mNlXZnZd1BjkcYTmVqbnOI0uBw45JJ7bWbcOLkrZ1ZpzzuVfhcFFUiNJh0l6gNCB16fAkXlPmas0KXSBLIXp556DF18sbJqcc/VTpqrIv5R0NzAHOAP4N7CdmR1jZk/VUPpcJe28M5x0Unz6wgtDLsY552pSppzL5cBbQFczO8zMHjCz5TWULlcNf/oTNGsWxj/+GO65p7Dpcc7VP2mDi5n1MbO/m9lCSXtLOhl+ata+Y80l0VXWllvCJZfEp6+8EpYsKVx6nHP1TzZlLr8ndD18WTSrEXB/PhPlqu+ii2CrrcL4d9/BNdcUNj3Oufolm4YrBwCHA8sBzGwu0CKfiXLV16wZjBgRnx45Ej77rHDpcc7VL9kElx8tfAxjAJKa5zdJLleOPTY0AwOwZo1XTXbO1Zxsgsujku4EWkk6HXgJ+Ht+k+VyQYK//S0+/cwz8N//Fi49zrn6I5uPKK8ndBn8T+DnwO/MzDvVrSXKykLTMDHnnx9yMc45l0+ZvnNRbNzMXjSzYWZ2kZm9mGodV7yuvTbe58v06XDLLYVNj3Ou7suUcxkv6VxJ2yTOlNRYUl9JY4ETq3JQSYMkTZW0XtIG7dJI2kbSMklpSwmitM2I9jMiYf5lkj6Plh1YlfTVNVtsAb/7XXz697+PN8/vnHP5kCm4HASsAx6SNFfSNEkzCa0lHwuMNLN7q3jcKcBA4NU0y0cSmppJSVIfoD/Qw8y6AddH83cgNLDZLUr/bZJKqpjGOuW88+KNWC5dWv47GOecy7VMH1GuMrPbzGwvYFtgf2BnM9vWzE43s8lVPaiZTTezGamWSToC+BKYmmEXZwF/MbPV0f6+j+b3J7TYvNrMZgKfA7tVNZ11SePGcNNN8en77oM33ihcepxzdVs2tcUwszVmNs/MFuUzMVE150uA4RWs2gXoLekdSa9I6hXN3wqYnbDenGheqmOdIWmSpEnz58+vbtJrhV/+Eo5MaHL0nHNg7drCpcc5V3dlFVyqQtJLkqakGPpn2Gw44XXbsgp23xDYFNgDGEaoLi0gVQWDlB3WmNloMyszs7K2bdtmcUZ1ww03wEYbhfHJk+H22wuaHOdcHdUwXzs2s35V2Gx34KiogL4VsF7SKjNLrt80B3gi+rjzXUnrgTbR/K0T1usAzK1COuqsbbeFK64I7Y1B+Peoo6B9+8KmyzlXt2TTttg5kjaticSYWW8zKzWzUmAUcG2KwALwFNA3Sl8XoDHwA/AMMFhSk6hxzc7AuzWQ9FrloougS5cwvmSJf7nvnMu9bF6LbQFMlPSopINy8W2LpAGS5gB7As9JeiGLbcYkVFu+G+gkaQrwMHCiBVOBR4FpwH+As83MezNJ0qQJ3HZbfPrBB+HllwuXHudc3aPwZqmClUJAOQA4GSgj3MDvMrMv8pu8mlFWVmaTJk0qdDJq3HHHwUMPhfGf/xw+/DAEHuecy4ak98xsg28VIfvaYgZ8Gw1rCYXpjyd+vOhqnxtugBZR+9YzZpRvRdk556ojmzKX/5P0HjACeAPobmZnAbsCR2bc2BW19u3L9/NyzTXw6aeFS49zru7IJufSBhhoZgea2WNmtgbAzNYDh+Y1dS7vfvtb6BV9JbR6NQwZAlm8KXXOuYyyCS6jgKWSNksYGkH40j6vqXN5V1ICo0eHfwHGjw9f7zvnXHVkE1zeB+YDnxLaFZsPzJT0vqRd85k4VzN69gxN8cdceCHUk0YLnHN5kk1w+Q9wsJm1MbPWwK8ItcV+C9yWcUtXawwfDttE7V8vWFA+2DjnXGVlE1zKzOyn71DM7L/APmb2NuAVV+uI5s3hjjvi0w8+CM89V7j0OOdqt2yCy0JJl0jaNhouBv4XNWW/Ps/pczXoV7+C44+PT591Vmie3znnKiub4HIcoY2up6Jh62heCXB0vhLmCmPUKGjTJozPng2XXVbQ5DjnaqmMwSXKnYwys3PNbOdoONfM5pvZj2b2eQ2l09WAESPg44/L9/ty663w2muV39f48f5RpnP1WcbgErXL1VZS4xpKjyugXr3g6KNh883hkEPi8085BVasyH4/48eH/cS+n3HO1T/ZvBabBbwh6SpJQ2NDntPlCqBPH3j0UTjmGPjNb6BlyzD/88/jTfRXJBZYHn007M85Vz9lE1zmAs9G67ZIGFwdFAswZ58NZ54Znz9qVMXdIntgcc7FZNUqMoQuiM1seZ7TUxD1tVXkTGKBorQUYpemc+fQe2WzZunX98DiXP1RrVaRJe0paRowPZreSZJ/PFnHxXIwX34Z7xb5s89CL5bJPLA455Jl27bYgcACADP7ENgnj2lyRaJPH3j88Xi7YxBej40fH5/2wOKcSyXb/lxmJ83y3h3riT594OmnoVGj+LyTTw7dI3tgcc6lk01wmS3pF4BJaizpIqJXZK5+6NsXHngAYh1cf/UVDB7sgcU5l142wWUIcDawFTAH6BlNu3pk0KDy1ZGffz40bumBxTmXSta1xeoyry2WvW7dYNq0MN62bfiif/PNC5sm51xhVLe2WFtJl0saLenu2FDNBA2SNFXSekkbJEzSNpKWRa/g0u3jXEkzov2MiOaVSlopaXI03JFue1d548fDd9/BxhuH6fnzw9f7/nzinEvWMIt1ngZeA14idwX5U4CBwJ1plo8Enk+3saQ+QH+gh5mtltQuYfEXZtYzR+l0kVjh/WOPwbp18Mtfhvn//jfcfnvoLtk552KyCS7NzOySXB401j2yYiXECSQdAXwJZPpg8yzgL2a2Otrf97lMnysvVa2woUPhxhvDeKzspWvXgiXROVdksinQf1bSwXlPCaEVAOASYHgFq3YBekt6R9IrkhKbSOwo6YNofu8MxzpD0iRJk+Z7n75ppatufM010L17GF+zBg47DFatKkwanXPFJ5vgch4hwKyStETSUklLKtpI0kuSpqQY+mfYbDgw0syWVbD7hsCmwB7AMOBRhWzQPGAbM9sZGAo8KKllqh2Y2WgzKzOzsrZt21Z0OvVSpu9YmjYNvVU2ifoi/eILOO64mk+jcy61ESPKf/BcHVXpQqPC4GJmLcysgZk1NbOW0XTKG3bSdv3MbMcUw9MZNtsdGCFpFnA+cLmkc1KsNwd4woJ3CT1itjGz1WYWa0ngPeALQi7HVVI2H0juuGP81RjAk0+GHI1zrvBiXWhUN8BUtQuNbGqLSdKvJV0VTW8tabeqJTMzM+ttZqVmVkpoduZaM7slxapPAX2j9HQBGgM/RDXbSqL5nYDOhPIbVwmV+fL+rLPgiCPi01ddFbZzzhVWrH3A6gSY6rTCkc1rsduAPQldGwMsA26t3GHKkzRA0pxov89JeiGLbcYkVFu+G+gkaQrwMHCihQ929gE+kvQh8DgwxMwWViet9U1lf0wS3HUXdOgQps3g17+Gl17KbzqdcxWrToCpdvNOZpZxAN6P/v0gYd6HFW1Xm4Zdd93VnNnLL5u1aRP+raxXXjFr0MAshBezZs2qth/nXO5V9m872/WBSZbmvppNzmVN9KrJIHxUSSjjcHXMxIlVf0rZZx+4+ur49IoVoT0y51zhVSYHk6sGabMJLjcBTwLtJF0DvA5cW/VDumJ18cXV+zFdfnn840oIrSl/80310+Wcq75sAkwuWzrPprbYA8DFwJ8JVX2PMLPHqndYVxeVlMD990P79mH6hx/gmGPCdzDOucLLFGBy3YWGN1yJN1yZa6+8EprpXx+9PD3/fBg5sqBJcs4liAWSfv1CIOncuWqBpVoNVzpXWfvuC9cmvDgdNcqrJztXTPbdF3r3hocfhjPPhP79c983kwcXlxcXX1z++5dTToHp3sWccwW3bl0IKE8+GZ+35Zaw3365PY4HF5cXEtx7L/zsZ2F6+XIYMAAWLy5ospyr19auhZNOgjFj4vO6doUFC2DChNwey4OLy5tNNoEnnoCNNgrTM2aEDyzXe0V252rc6tWhgs3998fnHXAAfPRR9b/kT8WDi8ur7t3h7oSu5Z59FoZX1Oa1cy6nVqwI5SpPPBGfd+ihobvyhg1z01RMMg8uLu8GD4Zhw+LTf/hD+fe9zrn8WbwYDjwQXkhoZOuoo+CZZ6BBQgTIdYDx4OJqxJ//XP4DyxNOCNlx51z+fP99CBqvvx6fd9JJIYik6KsxpwHGg4urESUlodpjp05hevny0MHY996HqHN58dVXsPfe8MEH8Xm//S3cc0/qwBKTqwDjwcXVmM02C1nxFi3C9Ndfw8CBoaDROZc706fDXnvBZ5/F5w0bBrdm2Z59LgKMBxdXo7p1g4ceij85vfFGqHPvDUU4lxtvvRVyLInt+l19deV7kqxugPHg4mrcIYfAX/8anx47tvwX/c65qnnuOdh/f1iY0IvV9dfD739ftf1VJ8B4cHEFMXRo+Go/5sorQ47GOVc1Y8eG6sYrV4bpZs3gttvgwgurt99YgJk4sXLbecOVeMOVhfLjj/CrX8HLL4fpxo1h3LiQpXfOZccs1Ma84or4vG23hf/+F7p0ye+xveFKV5QaN4Z//jM0PwEh2PTvH77kr47x4yv/ftm52mjdOjj77PKBpXt3ePPN/AeWinhwcQXVqlV4T9yqVZheuDB88DVvXtX2F2tKvFevXKXQueK0YgUceSTcfnt8Xp8+8NproSHKQvPg4gquY8eQhW/SJEx/9RUcfDAsWVK5/eS6syPnitV334VWjJ9+Oj7v2GNDcy6bbFKwZJVTkOAiaZCkqZLWS9rgfZ2kbSQtk3RRmu0fkTQ5GmZJmpyw7DJJn0uaIenAPJ6Gy6FevUK7R7HmKCZPrtw3MB5YXH3xySewxx7lC9iHDQsNUsYe0IpBoXIuU4CBwKtplo8Enk+3sZkdY2Y9zawn8E/gCQBJOwCDgW7AQcBtkkpymG6XRwcfDH//e3x63LjQivK6dZm388Di6ouXX4Y994RZs8J0gwbhw8gRI8q3E1YMCpIcM5tuZimLbSUdAXwJTK1oP5IEHA3EKrH2Bx42s9VmNhP4HNgtJ4l2NeKUU+CPf4xPP/44DBmS/iNLDyyumI0YkbtWhocNC+3zLVoUpps1C6/Ffvvb3Ow/14oq1klqDlwCZNsoe2/gOzOLNXKwFTA7YfmcaF6qY50haZKkSfPnz69qkl0eXHEFnHdefHrMGLj00g3X88Diil2vXtVvo2vdutCy+PXXx/tCat8eXn01NJtfrPIWXCS9JGlKiqF/hs2GAyPNbFmWhzmWeK4FIFVzbCmfec1stJmVmVlZ27ZtszycqwkS3HhjaDk5ZsQIuOaa+LQHFlcbVLcJlaVLQ1/3jzwSn9ezJ7z7Luy6a86SmRcN87VjM+tXhc12B46SNAJoBayXtMrMbkleUVJDQrlN4iWeA2ydMN0BmFuFdLgCa9AA7ror9EXxzDNh3pVXhl4td97ZA4urPRIDTGV+s7NmQd++MHNmfN5hh8GDD8LGG+clqTlVVK/FzKy3mZWaWSkwCrg2VWCJ9AM+MbM5CfOeAQZLaiKpI9AZeDefaXb506hReGLrl/CYcuGFcPjhHlhc7VLZHMwrr8BOO5UPLMOGhU72akNggcJVRR4gaQ6wJ/CcpBey2GZMUrXlwZR/JYaZTQUeBaYB/wHONrMK6hq5Yta0KTz1VHg1ELNsGXz+ecGS5FyVZBNgzEJ7YH37xr/zatwY7r03vBouqUV1X71tMbxtsdrg2WdhwABYuzY+b/RoOP30wqXJuapIV164ejWce2756vibbx6aSNprr5pPZza8bTFXq40fDyefHP7IyhJ+xmecUf4P0bnaIFUOZs4c2Gef8r/nsjKYNKl4A0tFPLi4opb4lHf44fDiixsGmGx713OuWCQGmFGjQs2vdxNKh3/961DVuEOHgiWx2jy4uKKV6vVBq1YbBphzzinf+ZhztcF++4WGJy+4AL7/Psxr0CAEm/vuCzUjazMPLq4oZfqOJRZg9tgjPu/ii0NXrl6E6GqDxYtD23l33hmf16xZaN7lvPPi3YDXZh5cXNHJ5gPJVq1CS8r77RefN3x46OEy9hWzc8Xogw9Czvupp+Lzttoq1IysS79dDy6uqFTmy/sWLUJfMAcmtH09ahSceCKsWZPXZDpXaWZwxx2h4cnEqvQDBsCXX4Z29KrbVEwx8eDiikZVmnSJNd531FHxefffH/5gV6zITzqdq6wlS0J/K2edVb4biauuCl1NNG5c/aZiio0HF1cUqtNWWJMm8PDDcOaZ8XnPPRf2Eysoda5Q3nkntAeW2D5YSQmMHQt/+EP5detSgPHg4orCxInVa9KlpCR095rYl/i774ZXEJ99ln475/Jl/Xr4y19g773LN+PStGn4KPg3v0m9XV0JMP6FPv6Ffl1z223hS+dY4Wjr1uHVWW39GM3VPrNnh+AxYUJ8XvPmoarx009n9xBVG1r+9i/0Xb3y29+G99ix7wQWLAhtNd1/f2HT5eqHhx+G7t3LB5Yddgivb7MNLFD7czAeXFyd1L9/+IOMddXz44+hf5grrqhb1T1d8Vi4EI47LhTcL14c5jVoEH53338faoNVNgdSmwOMBxdXZ+2+eyh36dYtPu/aa8NX0bEWZ53Lheefhx13hIcS2mkvLQ1V459/vnqvtmprgPHg4uq00lJ480341a/i8556Knzd/+mnue3jfPz4sD9XfyxeHFrmPvhgmDcvPv/EE+HDD2HlytyUmcQCzMSJ1dtPTfLg4uq8li1Db5ZDh8bnTZ8Ou+0WPrbMxRNhrPC1V6/q7cflRz4eImK5lTFj4svatQsPL/feG353F1+cu8L4Pn3C/moLDy6uXmjYEG64IRTqN20a5i1eHLpO3n//6gWY2lCrp77r1St3DxGDBsG4cSG3MiehH9yBA2HKlFDe5zy4uHrm+OPhjTdgm23i8x55BNq3D1/5V/bm44GldshFucXLL4fAsWZNaNcupk2b8Bt6/PF4BRLnwcXVQ7vsAu+/H548Yz7+OPRyecQR2d98PLDULtUJMPffDwcdBEuXlq8McswxMG1a2GddaMk4lzy4uHqpdWv417/gmmtCdVEIN40lS0LQeeGFzNt7YKmdKhtgVq+GU04J1YkTG0Pt0CH8fh5+2HMr6RQkuEgaJGmqpPWSNvi6U9I2kpZJuijN9o9ImhwNsyRNjuaXSlqZsOyOPJ+Kq8UaNIDLLw/vz7fcMj5/1So45JDyBbWJPLDUbtkGmOefh+22g3vuic+TQud006bBoYfmP621WaFyLlOAgcCraZaPBJ5Pt7GZHWNmPc2sJ/BP4ImExV/ElpnZkFwl2NVd++0Xqo0efnh83rp1oYrp6aeH8RgPLHVDpgDz5ZehbOXgg+Gbb+LzY10R33xz6O7BZVaQ4GJm081sRqplko4AvgSmVrQfSQKOBh6qaF3nMmnTJlQhvfXW8t3LjhkTqptOneqBpa5JDjBLlsAll0DXrqHqeswmm4SA8s475bvXdpk1LHQCEklqDlwC/BJI+UosSW/gOzNLbPe2o6QPgCXAlWb2Wu5T6uoiKbRL1q9faHTwnXfC/E8+gR49QhXmJ5/0wFKX9OkDDz4Ycq0NG8KiReWXn3pqaNWhXbuCJK9Wy1vORdJLkqakGDLVAh8OjDSzZVke5ljK51rmAduY2c7AUOBBSS3TpO8MSZMkTZo/f36Wh3P1QZcu8Prr8Kc/QaNGYd769aHzsaFD4TV/XKkTzMLDwrnnwrJl5QNL167h4WLMGA8sVVXQJvclTQAuMrNJ0fRrwNbR4lbAeuB3ZnZLim0bAt8Au5rZnOTlqfafjje579K55x4444xQTTnRySeHr7TbtClMulz1jB8fGjF9660Nlw0cGL5Z8arFFas1Te6bWW8zKzWzUmAUcG2qwBLpB3ySGFgktZVUEo13AjoTym+cq7Tx40NzG88/DzfeGP+yH0LQ6dwZbrqpfBVVV9zefju89uzbt3xgadYs9Ldy6aXw6qvlm8t3VVOoqsgDJM0B9gSek1TBVwUgaUxSteXBbFiQvw/wkaQPgceBIWa2MFfpdvVHYuF9v35wwQWhocu9946vs2gRnHce7LQT/Oc/4TWLK05vvgkHHhh6Jh03Lj6/ceOQU9loo/Ddyp//XDtbIC5KZlbvh1133dWci3n5ZbM2bcK/qVxzjVmDBmYhnMSHvn3NJk7M/jjXXZf+GFVJ83XX5WZf+VST57x+vdm4cWb9+m34f1VSYnbqqWYPPZT6/7qi34ALgEmW5r5a8Bt7MQweXFxMtjeV//zHrHlzs6ZNN7xxHX202bRpuTtWTe2nJtTEOa9da/bPf5rtttuG/zcNGpgdf7zZjBkVp6U2XddC8eDiwcVlobI3k5dfNttsM7NDDw1Pwok3Mcns2GPNpk/P7TFzvX0h5Oucly41u+kms06dUgeVE04IQaUyaaiN17cmeXDx4OIqUNWbSGy7e+81Gzhww5uaZHbkkZlfl1X32LXxxpfLc/70U7PzzzfbZJMNr3+TJmZDhph98UXVj12br3O+eXDx4OIyyOWT9Lvvmh188IY3uViZzL/+ZbZuXfXTUBdueNU559WrzR5/3OyAA1Jf6003Nbv0UrO5c6t3zOpuV9d5cPHg4tLIVxnA22+nDzLbbWc2cqTZwoVVS0t101xMFQkqe8533202bJhZu3apr23nzma33mq2bFnVj1XdtNYnHlw8uLg08n2j/fBDs+OOS127rGlTs9/8xuy110LNptg+8l3IXGwVCSraz2OPhcoTXbqkDiiS2WGHhUoWqXKFNZnW+saDiwcXV2CzZoUn7latUt8gt9vO7Pe/N/vss/Q3sFze2IrtKT55f99+a3bHHWY9e6a+XmC21VZmV11lNnNmxfsvptxaXeLBxYOLKxLLloWb5s47p79p7rKL2WmnhXKD2A0xH0/MxVT+sH692T33hBxKt24hN5Lq2jRubDZoUCi7WrMmd8d3VZMpuBS0bbFi4W2LuZpmBu+9B6NHhy/CFy9OvV5JCfTqBdOnh/UOOCC36ahsNwK57HZgwYLQzMpLL4Umdr76KvV6Euy7Lxx7LAwaBJtuWr3jutzJ1LaYBxc8uLjCWrUqND1y332he+V0bZU1awb77AP77x+aodlll9B8SXVlGzCqG1i++QbeeCMMr70GkydnbjJnm23gwgvhqKPK9xTqiocHlwp4cHHFYtGiEGgeeyy0V5apUcymTUOupqws9JK4yy6hMc2GVeilqaLAUZnAsn59yIVMmRJ6+Jw0KQyJvTqm0rJlaKftgw/gzDNh7FjvmK3YeXCpgAcXV2xiN/OhQ+Hrr+Hpp2HevIq3a9w49EfTtWsINB07QqdO0KEDtG+fuXvedAEkeb4ZLF8Oc+eGtM2eDTNnwmefhWHGjNA/SkVKSmC33cI+DzgAfvwRjjsufhzv+bP4eXCpgAcXV0zS3VQffDA80e+xR+jn/csqdCax8cahD5pNNw3DxhuHFoGbNQs5nnnzQhnI/vuH9b7+OrzG2n77EAwWLoTvv4eVKyt/7GbNQjDZa6/40LJl5nP2AFPcPLhUwIOLKxaVeT21/faht8QPPoD33w//VvTqqaa0bg3du8OOO4ZXdmVlITdVUrLhurl8JedqlgeXCnhwccUgFwXrixfDJ5+E2mVffhleV82cGV5hzZ0Lq1fnJq1NmsAWW4RC9623Dv927hyGLl1C18DZ9ORYU5UJXH54cKmABxdXaDVRJdgsVBhYuBD+978wLF8OK1aE11xr14Z1pk+Hu+8OvTWOHw/DhoXaaS1bQqtWIXBsvHH1uwEuZDVolxuZgkvBP2AshsE/onSFVEwfMybvM1/NnRTTObuqw7/Q9+DiilMxNcNSE83O5GJ/HmCKhwcXDy6uCBVTY4o10WBmMe7HVU+m4NKgJt/POefiJk7MTflBnz5hPxMnVm37bMoyYsc4+uiwflUVyzm7/CtIgb6kQcDVQFdgNzOblLR8G2AacLWZXZ9i+57AHUBTYC3wWzN7N1p2GXAqsA74PzN7oaL0eIG+q6+8UN1VR6YC/ULlXKYAA4FX0ywfCTyfYfsRwHAz6wn8LppG0g7AYKAbcBBwm6QUNeudc1UJFLnKwbi6ryDBxcymm9mMVMskHQF8CUzNtAsg+raXTYC50Xh/4GEzW21mM4HPgd1ykmjn6pDq5EA8wLhsFFWZi6TmwCXA8ApWPR/4q6TZwPXAZdH8rYDZCevNiealOtYZkiZJmjR//vxqpdu52iQXr7Y8wLiKVKH91OxIegnYIsWiK8zs6TSbDQdGmtkyZf5C6yzgAjP7p6SjgbuAfkCqjVIWKpnZaGB0lNb5kpJ7k2gD/JApEXVcfT9/qLPXYMvNYfmKvn0XL81i5QquwSYt+vZt3gzmfper1BWZOvobqJRM12DbdBvlLbiYWb8qbLY7cJSkEUArYL2kVWZ2S9J6JwLnReOPAWOi8TnA1gnrdSD+yixTWtsmz5M0KV1BVX1Q388f/BqAX4P6fv5Q9WtQVK/FzKy3mZWaWSkwCrg2RWCBEDD2jcb7Ap9F488AgyU1kdQR6Ay8m99UO+ecS5a3nEsmkgYANwNtgeckTTazAyvYZgxwR1Rt+XTgb5IaAquAMwDMbKqkRwnVmNcCZ5vZujyeinPOuRS84co0JJ0RlcvUS/X9/MGvAfg1qO/nD1W/Bh5cnHPO5VxRlbk455yrGzy4OOecyzkPLhFJJZI+kPRsNP1XSZ9I+kjSk5JaFTiJeZV8/gnzL5JkktoUKm01JdU1kHSupBmSpkZV5Ou0FH8HPSW9LWly9NFxnW7xQtIsSR/Hzjeat5mkFyV9Fv27aaHTmS9pzr9K90IPLnHnAdMTpl8EdjSzHsCnxFsBqKuSzx9JWwO/BL4uSIpqXrlrIKkPoUmhHmbWjdAaRF2X/DtI2Y5fHdfHzHomfNtxKTDOzDoD46Lpuiz5/Kt0L/TgAkjqABxC/GNMzOy/ZrY2mnyb8EFmnZTq/CMjgYtJ08pBXZLmGpwF/MXMVgOY2feFSFtNSXMN0rXjV5/0B8ZG42OBIwqXlJpX1XuhB5dgFOEmuj7N8lPI3EpzbTeKpPOXdDjwjZl9WKhE1bBRbPgb6AL0lvSOpFck9SpIymrOKDa8BueTuh2/usqA/0p6T9IZ0bzNzWweQPRvu4KlLv9SnX+irO+F9T64SDoU+N7M3kuz/ArCB5kP1GjCakiq85fUDLiC8BqkzsvwG2gIbArsAQwDHlUFjd7VVhmuQawdv62BCwjt+NVle5nZLsCvgLMl7VPoBNWwtOdf2XthQb7QLzJ7AYdLOpjQ+VhLSfeb2a8lnQgcCuxvdfeDoA3OH/gH0BH4MLqXdgDel7SbmX1bsJTmT8rfAKGtuiei//t3Ja0nNOJXF5vRTncNDiN1O351kpnNjf79XtKThC47vpPU3szmSWoP1NnXo2nO/9Uq3QvT9X9cHwdgP+DZaPwgQjMybQudrkKcf9L8WUCbQqevAL+BIcAfovEuhO4cVOg01vA1mA7sF43vD7xX6PTl8bybAy0Sxt+M7gN/BS6N5l8KjCh0Wmv4/Kt0L/ScS3q3AE2AF6On97fNbEhhk+Rq2N3A3ZKmAD8CJ1r0l1ePpGzHr47aHHgy+ntvCDxoZv+RNJHwSvRUQs3JQQVMYz6lO//PqcK90Jt/cc45l3P1vkDfOedc7nlwcc45l3MeXJxzzuWcBxfnnHM558HFOedcznlwcc45l3MeXJxzzuWcBxfn6gFJN0t6P9b4pqSuku6Q9LikswqdPlf3eHBxro6T1JzQku+ZhPahMLPp0VfWRwNlGTZ3rko8uDiXBUkTJB2YNO98Sbdl2GZZ/lO2wTE3iroHKInNM7PlQHtgAnBTwrqHA68TOsBCUmNJr0ZNvThXLR5cnMvOQ8DgpHmDo/nF5BRCS87rYjMktQaaAUuBn+ab2TNm9gvg+Gj6R0KgOaZGU+zqJA8uzmXnceBQSU0AJJUCWwKvSxoqaUo0nJ+8oaTSqPHL2PRFkq5OWPaJpDHR9g9I6ifpjajP9t2i9X4t6d2ob/M7E3MmSY4Hnk6adyWho6+pwA7R/vaTdJOkO4F/J6z7VLQP56rFg4tzWTCzBcC7hObHIeRaHgF2AU4Gdid0Kna6pJ0rufufAX8DegDbA8cBewMXAZdL6krITexloS/7daQIAJIaA53MbFbCvFLgF1FapwPdovOZYGb/Z2ZnmtmtCbuZAtT1HjddDfDg4lz2El+NxV6J7Q08aWbLzWwZ8ATQu5L7nWlmH5vZekLuYlzUtP/HQCmhH5VdgYmSJkfTnVLspw2wKGnenwh90hgJwSWd6HXaj5JaVPIcnCvHC+6cy95TwI2SdgE2MrP3s+wGdy3lH+SaJi1fnTC+PmF6PeFvVMBYM6uo//qVifuW1BMYCOwt6dZo2cdZpLcJoe8W56rMcy7OZSnKmUwgdCIWK8h/FThCUrOoyu8A4LWkTb8D2klqHZXZHFrJQ48DjpLUDkDSZpK2TZG+/wElkmIB5jrgMDMrNbNSYCcqyLlEhf/zzWxNJdPoXDmec3Guch4ivPoaDBDlXu4llMcAjDGzDxI3MLM1kv4AvAPMBD6pzAHNbJqkK4H/SmoArAHOBr5Ksfp/CTmV9UBzMxuXsJ/vJDWXtJmZLUxzuD6UL+B3rkq8J0rn6pCoMsFQMzuhits/AVxmZjNymzJX3/hrMefqkCjXND5DVeW0otpmT3lgcbngORfnnHM55zkX55xzOefBxTnnXM55cHHOOZdzHlycc87lnAcX55xzOefBxTnnXM79PwwynoLhVcSYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "murn.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "attached-palestinian",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "140.0186998321462"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "murn[\"output/equilibrium_bulk_modulus\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "injured-rainbow",
   "metadata": {},
   "source": [
    "### Same cane be done for the 1000 structures dataset\n",
    "The final parameters of the 100 structure fit can be used for the 1000 Structure fit. This can speed up the fitting process and often leads to better results, especially if the initially guessed values are far from the optimum. \n",
    "\n",
    "In general it is a good idea to start with few structures and try around with different functions and initial parameters. This is much faster than using all structures from the beginning and gives good guesses for the initial values of the parameters. It also allows to use global optimization with millions of steps in short time spans.\n",
    "\n",
    "### This can take long and writes 1000 seperate POSCAR files, TAKE CARE\n",
    "Run these jobs with more cores or more time after the workshop. Also increase the number of iterations for better results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "focal-rehabilitation",
   "metadata": {},
   "outputs": [
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "Are you sure you want to delete all jobs from 'PotentialDF2'? y/(n) n\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No jobs removed from 'PotentialDF2'.\n"
     ]
    }
   ],
   "source": [
    "pr = Project(\"PotentialDF2\")\n",
    "pr.remove_jobs()\n",
    "j = pr.create_job(pr.job_type.Atomicrex, \"PotentialDF2\")\n",
    "j.potential = job.potential.copy()\n",
    "## Use the final parameters as starting values for the new fit\n",
    "j.potential.copy_final_to_initial_params()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "bibliographic-wonder",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = data_pr.load(\"df2_1k\").to_pandas()\n",
    "for id, row in df.iterrows():\n",
    "    struct = ase_to_pyiron(row.atoms)\n",
    "    s = j.structures.add_structure(struct, f\"id{id}\", relative_weight=1)\n",
    "    s.fit_properties.add_FitProperty(\"atomic-energy\", target_value=row.energy/row.number_of_atoms, relative_weight=1)\n",
    "    s.fit_properties.add_FitProperty(\"atomic-forces\", target_value=row.forces, relative_weight=1/row.number_of_atoms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "taken-remove",
   "metadata": {},
   "outputs": [],
   "source": [
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "loved-princess",
   "metadata": {},
   "outputs": [],
   "source": [
    "j.input.atom_types.Cu = None\n",
    "j.input.fit_algorithm = j.factories.algorithms.ar_lbfgs(max_iter=100000)\n",
    "\n",
    "## if possible increase number of cores\n",
    "#j.server.cores = 16\n",
    "t1 = time.time()\n",
    "## Uncomment if you want to run the job\n",
    "#j.run()\n",
    "t2 = time.time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "sunrise-brother",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.1457672119140625e-05"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t2-t1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "therapeutic-treasure",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/json": {
       "error": "None",
       "iterations": "array([   1,    2,    3, ..., 5594, 5595, 5596], dtype=uint32)",
       "residual": "array([758.612 , 758.612 , 758.612 , ...,  58.7461,  58.7461,  58.7461])"
      },
      "text/plain": [
       "Output({'error': None, 'iterations': array([   1,    2,    3, ..., 5594, 5595, 5596], dtype=uint32), 'residual': array([758.612 , 758.612 , 758.612 , ...,  58.7461,  58.7461,  58.7461])})"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "j.output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "composite-porter",
   "metadata": {},
   "source": [
    "This is the result if the initilly guessed values are taken instead of the fitted ones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "regulated-document",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The job PotentialDF2_BadStartParams was saved and received the ID: 832\n"
     ]
    }
   ],
   "source": [
    "j2 = pr.create_job(pr.job_type.Atomicrex, \"PotentialDF2_BadStartParams\", delete_existing_job=True)\n",
    "j2.potential = job.potential.copy()\n",
    "j2.input.atom_types.Cu = None\n",
    "j2.input.fit_algorithm = j.factories.algorithms.ar_lbfgs(max_iter=100000)\n",
    "j2.structures = j.structures\n",
    "## if possible increase number of cores\n",
    "j2.server.cores = 16\n",
    "## Uncomment if you want to run the job\n",
    "j2.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "greek-infrastructure",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/json": {
       "error": "None",
       "iterations": "array([   1,    2,    3, ..., 5289, 5290, 5291], dtype=uint32)",
       "residual": "array([5717.28  , 5717.28  , 5717.28  , ...,   58.7461,   58.7461,\n         58.7461])"
      },
      "text/plain": [
       "Output({'error': None, 'residual': array([5717.28  , 5717.28  , 5717.28  , ...,   58.7461,   58.7461,\n",
       "         58.7461]), 'iterations': array([   1,    2,    3, ..., 5289, 5290, 5291], dtype=uint32)})"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "j2.output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "twenty-collins",
   "metadata": {},
   "source": [
    "With this choice of functions and initial parameters starting directly from all structures gives the same residual. In a previous iteration of the potential it was about 7 times worse, so it is a good idea to test this."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}