diff --git a/day_2/atomicrex/WorkshopPotentialEAM.ipynb b/day_2/atomicrex/WorkshopPotentialEAM.ipynb index d9e9888c544af43db8728bb2db642be598d8ff1b..678c7c901ff55d6219d4e27e1ab7ab83f95c429f 100644 --- a/day_2/atomicrex/WorkshopPotentialEAM.ipynb +++ b/day_2/atomicrex/WorkshopPotentialEAM.ipynb @@ -1,22 +1,55 @@ { "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": "political-robinson", + "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": "tracked-postcard", + "id": "constant-respect", "metadata": {}, "outputs": [ { @@ -60,54 +93,54 @@ " <tbody>\n", " <tr>\n", " <th>0</th>\n", - " <td>619</td>\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>import_database/Cu_database/</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>GenericJob</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>620</td>\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>import_database/Cu_database/</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>GenericJob</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>621</td>\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>import_database/Cu_database/</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>GenericJob</td>\n", + " <td>TrainingContainer</td>\n", " <td>0.4</td>\n", " <td>None</td>\n", " <td>None</td>\n", @@ -118,24 +151,29 @@ ], "text/plain": [ " id status chemicalformula job \\\n", - "0 619 finished None df1_A1_A2_A3_EV_elast_phon \n", - "1 620 finished None df3_10k \n", - "2 621 finished None df2_1k \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 timestart timestop \\\n", - "0 import_database/Cu_database/ 2021-02-08 10:33:52.341472 None \n", - "1 import_database/Cu_database/ 2021-02-08 10:33:53.993230 None \n", - "2 import_database/Cu_database/ 2021-02-08 10:33:54.435308 None \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", - " totalcputime computer hamilton hamversion parentid masterid \n", - "0 None zora@cmti001#1 GenericJob 0.4 None None \n", - "1 None zora@cmti001#1 GenericJob 0.4 None None \n", - "2 None zora@cmti001#1 GenericJob 0.4 None None " + " 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, @@ -144,16 +182,16 @@ } ], "source": [ - "data_pr = Project(\"../../datasets\")\n", + "data_pr = Project(\"../../datasets/imported_datasets/\")\n", "if len(data_pr.job_table()) == 0:\n", - " data_pr.unpack(\"Cu_training_archive\")\n", + " data_pr.unpack(\"../../datasets/Cu_training_archive\")\n", "data_pr.job_table()" ] }, { "cell_type": "code", "execution_count": 3, - "id": "authentic-substance", + "id": "dirty-measurement", "metadata": {}, "outputs": [], "source": [ @@ -164,7 +202,7 @@ { "cell_type": "code", "execution_count": 4, - "id": "vertical-simon", + "id": "referenced-julian", "metadata": {}, "outputs": [ { @@ -176,13 +214,22 @@ } ], "source": [ - "pr = Project(\"WorkshopPotential\")" + "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": "finite-planner", + "id": "entertaining-jacksonville", "metadata": {}, "outputs": [], "source": [ @@ -191,7 +238,7 @@ }, { "cell_type": "markdown", - "id": "fiscal-advocate", + "id": "raising-clear", "metadata": {}, "source": [ "### Add the structures that should be fitted.\n", @@ -201,7 +248,7 @@ { "cell_type": "code", "execution_count": 6, - "id": "funded-offense", + "id": "located-individual", "metadata": {}, "outputs": [], "source": [ @@ -214,7 +261,7 @@ }, { "cell_type": "markdown", - "id": "simplified-withdrawal", + "id": "angry-leader", "metadata": {}, "source": [ "### Define the type of potential and necessary functions.\n", @@ -224,7 +271,7 @@ { "cell_type": "code", "execution_count": 7, - "id": "worth-electricity", + "id": "functional-formation", "metadata": {}, "outputs": [], "source": [ @@ -233,40 +280,47 @@ }, { "cell_type": "markdown", - "id": "earlier-arrow", + "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 are not directly physically motivated and can lead to unphysical oscillations or overfitting.\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, while the electronic density and embedding function will be splines. Depending on the properties that are calculated other functional forms could give better results.\n", + "In this case a generalized morse function is used for the pair interaction. It has the form\n", "\n", - "The parameters in D0=3.5 and r0=1.8 are the approximate cohesive energy and the equilibrium lattice constant. Beta and S can also be derived from physical quantities but are chosen randomly in a typical range in this case. Delta is a parameter that shifts the whole function up or down. The initial parameter choices should not matter too much as long as they are somewhat reasonable since they will be optimized in the fitting process anyway." + "$(\\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": "grave-settlement", + "id": "interpreted-orange", "metadata": {}, "outputs": [], "source": [ - "V = job.factories.functions.morse_B(identifier=\"V_CuCu\", D0=3.5, r0=1.8, beta=2, S=2, delta=0)" + "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": "sharp-photographer", + "id": "mathematical-gasoline", "metadata": {}, "outputs": [], "source": [ "V.parameters.D0.min_val = 0\n", - "V.parameters.D0.max_val = 5\n", - "V.parameters.r0.min_val = 1\n", - "V.parameters.r0.max_val = 2.5\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", @@ -275,7 +329,7 @@ }, { "cell_type": "markdown", - "id": "rough-purchase", + "id": "written-commission", "metadata": {}, "source": [ "Additionally a screening function needs to be defined for the morse potential" @@ -284,7 +338,7 @@ { "cell_type": "code", "execution_count": 10, - "id": "practical-details", + "id": "discrete-terminology", "metadata": {}, "outputs": [], "source": [ @@ -293,16 +347,16 @@ }, { "cell_type": "markdown", - "id": "external-ready", + "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." + "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": "departmental-dynamics", + "id": "authentic-expression", "metadata": {}, "outputs": [], "source": [ @@ -311,7 +365,7 @@ }, { "cell_type": "markdown", - "id": "latter-wright", + "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", @@ -321,7 +375,7 @@ { "cell_type": "code", "execution_count": 12, - "id": "double-engineering", + "id": "hidden-wildlife", "metadata": {}, "outputs": [ { @@ -342,7 +396,7 @@ }, { "cell_type": "markdown", - "id": "handled-housing", + "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." @@ -351,7 +405,7 @@ { "cell_type": "code", "execution_count": 13, - "id": "developing-apache", + "id": "comparative-brush", "metadata": {}, "outputs": [], "source": [ @@ -361,7 +415,7 @@ }, { "cell_type": "markdown", - "id": "separated-journal", + "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", @@ -372,7 +426,7 @@ { "cell_type": "code", "execution_count": 14, - "id": "colored-discount", + "id": "funny-trinidad", "metadata": {}, "outputs": [], "source": [ @@ -383,7 +437,7 @@ }, { "cell_type": "raw", - "id": "preliminary-peace", + "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." @@ -392,7 +446,7 @@ { "cell_type": "code", "execution_count": 15, - "id": "lasting-reply", + "id": "mexican-absence", "metadata": {}, "outputs": [], "source": [ @@ -402,7 +456,7 @@ }, { "cell_type": "markdown", - "id": "saving-vegetable", + "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", @@ -414,7 +468,7 @@ { "cell_type": "code", "execution_count": 16, - "id": "overall-flower", + "id": "large-rating", "metadata": {}, "outputs": [], "source": [ @@ -432,7 +486,7 @@ }, { "cell_type": "markdown", - "id": "neutral-gateway", + "id": "several-mercy", "metadata": {}, "source": [ "The functions have to be assigned to the potential" @@ -441,7 +495,7 @@ { "cell_type": "code", "execution_count": 17, - "id": "major-capacity", + "id": "heavy-acoustic", "metadata": {}, "outputs": [], "source": [ @@ -452,7 +506,7 @@ }, { "cell_type": "markdown", - "id": "unusual-retirement", + "id": "alien-chancellor", "metadata": {}, "source": [ "### Define fitting procedure\n", @@ -463,25 +517,26 @@ { "cell_type": "code", "execution_count": 18, - "id": "productive-spare", + "id": "enormous-segment", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "The job PotentialDF1 was saved and received the ID: 622\n" + "The job PotentialDF1 was saved and received the ID: 819\n" ] }, { "data": { "application/json": { "error": "None", - "iterations": "19588", - "residual": "0.162222" + "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': 0.162222, 'iterations': 19588})" + "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, @@ -492,43 +547,235 @@ "source": [ "## Define the atom types of the potential\n", "job.input.atom_types.Cu = None\n", - "##\n", - "job.input.fit_algorithm = job.factories.algorithms.ar_lbfgs(max_iter=500)\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": "abandoned-reform", + "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": [ - "### Same for the 1000 structures dataset\n", - "The final parameters of the 100 structure fit can be used for the 1000 Structure fit. This speeds up the fitting process and often leads to better results.\n", + "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" + "### 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": 19, - "id": "preceding-march", + "execution_count": 24, + "id": "focal-rehabilitation", "metadata": {}, - "outputs": [], + "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()\n", - "j.input = job.input.copy()" + "j.potential.copy_final_to_initial_params()" ] }, { "cell_type": "code", - "execution_count": 20, - "id": "ancient-september", + "execution_count": 25, + "id": "bibliographic-wonder", "metadata": {}, "outputs": [], "source": [ @@ -542,8 +789,8 @@ }, { "cell_type": "code", - "execution_count": 21, - "id": "correct-israeli", + "execution_count": 26, + "id": "taken-remove", "metadata": {}, "outputs": [], "source": [ @@ -552,37 +799,35 @@ }, { "cell_type": "code", - "execution_count": 22, - "id": "multiple-harbor", + "execution_count": 27, + "id": "loved-princess", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "The job PotentialDF2 was saved and received the ID: 623\n" - ] - } - ], + "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", - "j.run()\n", + "## Uncomment if you want to run the job\n", + "#j.run()\n", "t2 = time.time()" ] }, { "cell_type": "code", - "execution_count": 23, - "id": "eligible-soviet", + "execution_count": 28, + "id": "sunrise-brother", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "538.0122790336609" + "2.1457672119140625e-05" ] }, - "execution_count": 23, + "execution_count": 28, "metadata": {}, "output_type": "execute_result" } @@ -593,22 +838,22 @@ }, { "cell_type": "code", - "execution_count": 24, - "id": "subsequent-diploma", + "execution_count": 29, + "id": "therapeutic-treasure", "metadata": {}, "outputs": [ { "data": { "application/json": { "error": "None", - "iterations": "3619", - "residual": "117.878" + "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, 'residual': 117.878, 'iterations': 3619})" + "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": 24, + "execution_count": 29, "metadata": {}, "output_type": "execute_result" } @@ -619,7 +864,7 @@ }, { "cell_type": "markdown", - "id": "collaborative-charleston", + "id": "composite-porter", "metadata": {}, "source": [ "This is the result if the initilly guessed values are taken instead of the fitted ones." @@ -627,44 +872,49 @@ }, { "cell_type": "code", - "execution_count": 25, - "id": "dutch-module", + "execution_count": 33, + "id": "regulated-document", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "The job PotentialDF2_BadStartParams was saved and received the ID: 624\n" + "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 = j.input.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": 26, - "id": "innocent-pasta", + "execution_count": 34, + "id": "greek-infrastructure", "metadata": {}, "outputs": [ { "data": { "application/json": { "error": "None", - "iterations": "607", - "residual": "8269.14" + "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': 8269.14, 'iterations': 607})" + "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": 26, + "execution_count": 34, "metadata": {}, "output_type": "execute_result" } @@ -674,12 +924,12 @@ ] }, { - "cell_type": "code", - "execution_count": null, - "id": "handed-tribute", + "cell_type": "markdown", + "id": "twenty-collins", "metadata": {}, - "outputs": [], - "source": [] + "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": {