Unverified Commit 6dd414bc authored by Jan Janssen's avatar Jan Janssen Committed by GitHub
Browse files

Merge pull request #48 from pyiron/updates

Updates
parents a41b5fea ae1a5ae5
%% Cell type:code id:political-robinson tags:
%% Cell type:markdown id:entitled-break tags:
# Fitting an EAM potential
EAM potentials are pair functionals.
In a generalised form they are equal to Finnis-Sinclair, effective medium theory or glue potentials. Their total energy can be written as
$E = \sum_{ij}V(r_{ij}) + \sum_i F(\rho_i)$
with
$\rho_i = \sum_j \rho(r_{ij})$
The original functions for V, $\rho$ and F were derived from different theories, but they can be chosen freely.
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
$\chi^2 = \sum_i w_i r_i$
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 id:boring-simon tags:
``` python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pyiron import Project, ase_to_pyiron
```
%% Cell type:code id:tracked-postcard tags:
%% Cell type:markdown id:southwest-southwest tags:
### Import the training data
%% Cell type:code id:adaptive-yacht tags:
``` python
data_pr = Project("../../datasets")
if len(data_pr.job_table()) == 0:
data_pr.unpack("Cu_training_archive")
data_pr.job_table()
```
%%%% Output: execute_result
id status chemicalformula job \
0 619 finished None df1_A1_A2_A3_EV_elast_phon
1 620 finished None df3_10k
2 621 finished None df2_1k
0 759 finished None df1_A1_A2_A3_EV_elast_phon
1 760 finished None df3_10k
2 761 finished None df2_1k
subjob projectpath \
0 /df1_A1_A2_A3_EV_elast_phon /home/niklas/pyiron/projects/
1 /df3_10k /home/niklas/pyiron/projects/
2 /df2_1k /home/niklas/pyiron/projects/
project timestart timestop \
0 import_database/Cu_database/ 2021-02-08 10:33:52.341472 None
1 import_database/Cu_database/ 2021-02-08 10:33:53.993230 None
2 import_database/Cu_database/ 2021-02-08 10:33:54.435308 None
project \
0 pyiron_potentialfit/datasets/imported_datasets/Cu_database/
1 pyiron_potentialfit/datasets/imported_datasets/Cu_database/
2 pyiron_potentialfit/datasets/imported_datasets/Cu_database/
timestart timestop totalcputime computer \
0 2021-02-08 10:33:52.341472 None None zora@cmti001#1
1 2021-02-08 10:33:53.993230 None None zora@cmti001#1
2 2021-02-08 10:33:54.435308 None None zora@cmti001#1
totalcputime computer hamilton hamversion parentid masterid
0 None zora@cmti001#1 GenericJob 0.4 None None
1 None zora@cmti001#1 GenericJob 0.4 None None
2 None zora@cmti001#1 GenericJob 0.4 None None
hamilton hamversion parentid masterid
0 TrainingContainer 0.4 None None
1 TrainingContainer 0.4 None None
2 TrainingContainer 0.4 None None
%% Cell type:code id:authentic-substance tags:
%% Cell type:code id:whole-dragon tags:
``` python
data_job = data_pr.load("df1_A1_A2_A3_EV_elast_phon")
df = data_job.to_pandas()
```
%% Cell type:code id:vertical-simon tags:
%% Cell type:code id:portuguese-trance tags:
``` python
pr = Project("WorkshopPotential")
#pr.remove_jobs()
```
%%%% Output: stream
Are you sure you want to delete all jobs from 'WorkshopPotential'? y/(n) y
%% Cell type:code id:finite-planner tags:
%% Cell type:markdown id:vanilla-asset tags:
### Create an atomicrex job
%% Cell type:code id:characteristic-corruption tags:
``` python
job = pr.create_job(pr.job_type.Atomicrex, "PotentialDF1")
```
%% Cell type:markdown id:fiscal-advocate tags:
%% Cell type:markdown id:competent-produce tags:
### Add the structures that should be fitted.
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 id:funded-offense tags:
%% Cell type:code id:streaming-distance tags:
``` python
for id, row in df.iterrows():
struct = ase_to_pyiron(row.atoms)
s = job.structures.add_structure(struct, f"id{id}", relative_weight=1)
s.fit_properties.add_FitProperty("atomic-energy", target_value=row.energy/row.number_of_atoms, relative_weight=1)
s.fit_properties.add_FitProperty("atomic-forces", target_value=row.forces, relative_weight=1/row.number_of_atoms)
```
%% Cell type:markdown id:simplified-withdrawal tags:
%% Cell type:markdown id:acoustic-developer tags:
### Define the type of potential and necessary functions.
In this case an eam potential is fitted.
%% Cell type:code id:worth-electricity tags:
%% Cell type:code id:compliant-pharmacy tags:
``` python
job.potential = job.factories.potentials.eam_potential()
```
%% Cell type:markdown id:earlier-arrow tags:
%% Cell type:markdown id:respective-writer tags:
It is necessary to define a pair potential, an electronic density function and an embedding function.
For all of those it is possible to choose between different functional forms.
Classic pair potentials are physically motivated and have a very limited number of paramaters that are derived from a experimentally measured quantity.
Splines or polynomials offer more flexibility, but are not directly physically motivated and can lead to unphysical oscillations or overfitting.
Classic pair potentials are physically motivated and have a very limited number of paramaters that can often be derived from an experimentally measurable quantity.
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.
In this case a generalized morse function is used for the pair interaction. It has the form
$(\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 $
The parameters in the morse potential can be derived from phyiscal quantities, here 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.
In the case of analytic functions the initial parameter choices should not matter too much, since the functional form is constrained.
The electronic density and embedding function will be splines. Their initial parameters require more testing and hand tuning than the parameters of analytic functions. Depending on the properties that are calculated other functional forms could give better results.
%% Cell type:code id:electric-corps tags:
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.
``` python
V = job.factories.functions.morse_B(identifier="V_CuCu", D0=0.35, r0=2.5, beta=2, S=2, delta=0)
```
%% Cell type:markdown id:retired-motor tags:
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.
Pre defined functions like the morse function can be plotted to see the influence of the initial parameter values
%% Cell type:code id:grave-settlement tags:
%% Cell type:code id:affected-tamil tags:
``` python
V = job.factories.functions.morse_B(identifier="V_CuCu", D0=3.5, r0=1.8, beta=2, S=2, delta=0)
V.plot()
```
%% Cell type:code id:sharp-photographer tags:
%%%% Output: execute_result
(<Figure size 720x504 with 1 Axes>,
<AxesSubplot:xlabel='r [$\\AA$]', ylabel='func(r)'>)
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAmMAAAGzCAYAAABjHhDPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAlCUlEQVR4nO3deXzc913n8fdnJI2O0X34lB3FTpw096G2SdOkVwoplLTbFpbwoCzHYuDB0S48OMsuUJbHwsIWdqELhLa7LdsDaJstvWjTbUt6JGnsHE4cO47t2LF8SpZlybo189k/Zmwrjm2NdX1m5vd6Ph56SDPz08yn08R+5Xd8x9xdAAAAiJGKHgAAACDJiDEAAIBAxBgAAEAgYgwAACAQMQYAABCIGAMAAAgUFmNmVmdm3zOzp8xsu5n9QdQsAAAAUSxqnTEzM0kZdz9lZjWSvi3pPe7+SMhAAAAAAaqjXtjzFXiqcLOm8MUKtAAAIFHCYkySzKxK0lZJV0j6oLs/ep5tNkvaLEmZTObWq6++enmHRNnbf3xMkzNZbVrZFD0KACBBtm7dOuDuXXNtF3aY8iVDmLVKekDSL7v7Mxfarre317ds2bJsc6Ey/NZntunrO4/pe++7O3oUAECCmNlWd++da7uSuJrS3YckfVPSPbGToBK1NNRoaGxapfAfHgAAnCvyasquwh4xmVm9pLsl7YyaB5WrrSGtqWxO49PZ6FEAAHiZyHPGVkv6aOG8sZSkf3T3LwTOgwrVWl8jSRoam1ZDOvQ0SQAAXibyasptkm6Oen0kR2tDPsZOjE1pTWt98DQAALxUSZwzBiyl1oa0JOnk2HTwJAAAvBwxhop3ds8YMQYAKD3EGCpeW2HP2ND4VPAkAAC8HDGGitcy6wR+AABKDTGGildXU6W6mpSGxtgzBgAoPcQYEqGtIc2eMQBASSLGkAgt9TWcwA8AKEnEGBKhtaGGw5QAgJJEjCER2jNpDY2zZwwAUHqIMSRCW0Nag6PsGQMAlB5iDInQnklraGxK2ZxHjwIAwEsQY0iEtoa0ci4Nc6gSAFBiiDEkQnsmvwr/ICfxAwBKDDGGRDgdYyc4bwwAUGKIMSTCmT1jxBgAoMQQY0iEttN7xjhMCQAoMcQYEqG94fSeMU7gBwCUFmIMiVCfzn9YOHvGAAClhhhDYrSz8CsAoAQRY0iMtgwxBgAoPcQYEqOdGAMAlCBiDInR1pDmnDEAQMkhxpAY7BkDAJQiYgyJ0daQ1sjEjKazuehRAAA4gxhDYrRnaiSx8CsAoLQQY0iM9kytJOkEC78CAEoIMYbEaCvsGeO8MQBAKSHGkBh8WDgAoBQRY0iMMzHGOWMAgBJCjCExTn9Y+MDIZPAkAACcRYwhMaqrUmprqNHxUWIMAFA6iDEkSkdjrY6f4jAlAKB0EGNIlI5MmhgDAJQUYgyJ0tlYqwEOUwIASggxhkTpaGTPGACgtBBjSJSOTK1Ojk9raobPpwQAlAZiDInS2ZRf3oLPpwQAlApiDInSUfh8yoFTnDcGACgNxBgSpbMxv2eM88YAAKWCGEOidDTm94yx8CsAoFQQY0iUDvaMAQBKDDGGRGmqrVa6KqUBYgwAUCKIMSSKmRXWGuMwJQCgNBBjSJyOxjRXUwIASgYxhsTpyNTq+CiHKQEApYEYQ+J0NtZyAj8AoGQQY0iczsa0+k9Nyt2jRwEAgBhD8nQ11WpqJqeRyZnoUQAAIMaQPF1N+YVfjw1zEj8AIB4xhsTpKqzC3z9CjAEA4hFjSJwVzYUYY3kLAEAJIMaQOF2NdZKkY8MTwZMAAECMIYGa66uVrk6xZwwAUBKIMSSOmamrsZZzxgAAJYEYQyJ1NRFjAIDSEBZjZrbOzL5hZjvMbLuZvSdqFiQPMQYAKBWRe8ZmJP2au79C0m2SftHMrgmcBwmyghgDAJSIsBhz98Pu/njh5xFJOyStjZoHydLVlP+w8OlsLnoUAEDClcQ5Y2bWI+lmSY+e57HNZrbFzLb09/cv+2yoTKdX4ecDwwEA0cJjzMwaJX1G0nvdffjcx939fnfvdfferq6u5R8QFWlFU36tMQ5VAgCihcaYmdUoH2Ifd/fPRs6CZDnz+ZQjLPwKAIgVeTWlSfqwpB3u/oGoOZBMp2OMPWMAgGiRe8bukPRuSW80sycLXz8QOA8SpLMxLYkYAwDEq456YXf/tiSLen0kW211lVobanSUw5QAgGDhJ/ADUVY11+noMHvGAACxiDEk1ormOh0dZs8YACAWMYbEWtVcqyMniTEAQCxiDIm1qrlOA6cmNcMq/ACAQMQYEmtFc51yLg2wCj8AIBAxhsRa1Zxfhf8I540BAAIRY0isVS2FGOO8MQBAIGIMibWimY9EAgDEI8aQWJ2ZWlWnjD1jAIBQxBgSK5UyrWiq5ZwxAEAoYgyJtqK5TsdYhR8AEIgYQ6Ktaq5jzxgAIBQxhkRb1VKno5wzBgAIRIwh0VY212lkckajkzPRowAAEooYQ6KtLCxvwaFKAEAUYgyJtrqlXpJ0eIgYAwDEIMaQaGtb8zF26OR48CQAgKQixpBoK1vyhykPDRFjAIAYxBgSrba6Sl1NtRymBACEIcaQeGta6jhMCQAIQ4wh8da01nOYEgAQhhhD4q1uqdehoQm5e/QoAIAEIsaQeGta6zQ+ndXJ8enoUQAACUSMIfHWFJa3OMihSgBAAGIMiXc6xriiEgAQgRhD4q1pqZPEwq8AgBjEGBKvs7FWNVWmQ+wZAwAEIMaQeKmUaVVLHctbAABCEGOApDUt9ZzADwAIQYwBkrrbGnTwBDEGAFh+xBggaV17vY6OTGhyJhs9CgAgYYgxQPk9Y+7iJH4AwLIjxgBJ69rya40dGBwLngQAkDTEGCCpu71BktTHeWMAgGVGjAGSVjXXqTplOnCCPWMAgOVFjAGSqlKmNa317BkDACw7Ygwo6G6rVx97xgAAy4wYAwrWtTXowCB7xgAAy4sYAwq62+o1cGpSE9OsNQYAWD7EGFCw7swVlRyqBAAsH2IMKOg+s9YYhyoBAMuHGAMK1nfk94y9yMKvAIBlRIwBBV2NtWpIV2nf8dHoUQAACUKMAQVmpvXtDXrxOHvGAADLhxgDZunpyLBnDACwrIgxYJbLOvNrjWVzHj0KACAhiDFglp6OjKayOR0+yRWVAIDlQYwBs1xWuKJyP+eNAQCWCTEGzHJZR0YSMQYAWD7EGDDL6uY6patT2s9J/ACAZUKMAbOkUvnlLbiiEgCwXIgx4Bw9HQ0cpgQALBtiDDjHZYW1xnIsbwEAWAbEGHCODV0ZTUzndIjlLQAAyyA0xszsI2Z2zMyeiZwDmG1jV6MkaW8/540BAJZe9J6x/y3pnuAZgJc4HWN7+k8FTwIASILQGHP3hyQNRs4AnKuzMa2mumr2jAEAlkX0nrE5mdlmM9tiZlv6+/ujx0ECmJk2djWyZwwAsCxKPsbc/X5373X33q6uruhxkBAbujLsGQMALIuSjzEgwsauRh0ZntCpyZnoUQAAFY4YA85jY1f+MypfYO8YAGCJRS9t8UlJD0u6ysz6zOxnIucBTuOKSgDAcqmOfHF3vy/y9YELuawjo6qUEWMAgCXHYUrgPNLVKfV0NGjX0ZHoUQAAFY4YAy5g08om7TrKnjEAwNIixoAL2LSySfuOj2piOhs9CgCgghFjwAVctapJ7tLuY+wdAwAsHWIMuIBNK5skifPGAABLihgDLqCno0HpqpSeI8YAAEuIGAMuoLoqpQ1dGe06QowBAJYOMQZcxFWruKISALC0iDHgIjatbNLBoXGNTExHjwIAqFDEGHARr1idP4l/J4cqAQBLhBgDLuKa1S2SpGcPDQdPAgCoVMQYcBErm2vVnklr+6GT0aMAACoUMQZchJnp2jXNevYwe8YAAEuDGAPmcM3qZu06ckrT2Vz0KACACkSMAXO4Zk2zprI5PhYJALAkiDFgDteuaZYkbeckfgDAEiDGgDlc3tmoupoUV1QCAJYEMQbMoSplesXqZj1zkCsqAQCLjxgDinBjd6uePnhSM5zEDwBYZMQYUIQb17VofDqr3f2cxA8AWFzEGFCEG7tbJUlPHRgKnQMAUHmIMaAIPR0ZNddV68kDnDcGAFhcxBhQhFTKdOO6Vm3rG4oeBQBQYYgxoEg3drdq55ERTUxno0cBAFQQYgwo0g3dLcrmnCUuAACLihgDinTz+jZJ0tb9J4InAQBUEmIMKFJXU616Ohq0hRgDACwiYgy4BL097dqyb1DuHj0KAKBCEGPAJXhlT5tOjE1rT/9o9CgAgApBjAGX4NbL2iVJW/cPBk8CAKgUxBhwCTZ2ZdTWUKPH9nHeGABgcRBjwCUwM/X2tOuxfewZAwAsDmIMuES3bejQ/uNjOjg0Hj0KAKACEGPAJXrNxg5J0nd3DwRPAgCoBMQYcImuWtmkjkxaD+85Hj0KAKACEGPAJUqlTLdt7NB39gyw3hgAYMGIMWAe7tjYqaPDk6w3BgBYMGIMmIcz543t4bwxAMDCEGPAPFzW0aD17Q361+f6o0cBAJQ5YgyYBzPTG67q0nf2DGhiOhs9DgCgjBFjwDy9/uoVmpjO6dEXWAAWADB/xBgwT7dv6FBtdUrf2HksehQAQBkjxoB5qqup0ms2duibzxFjAID5I8aABXjj1Su07/iYdh8biR4FAFCm5owxM7vdzD5oZtvMrN/MXjSzL5nZL5pZy3IMCZSqN1+zSpL0L88cCZ4EAFCuLhpjZvZlSf9e0lck3SNptaRrJP2upDpJnzOze5d6SKBUrWqp0y3rW/VlYgwAME/Vczz+bnc/d1XLU5IeL3z9NzPrXJLJgDLxlutW64++tEMvHh/T+o6G6HEAAGXmonvG3H3AzKrM7GsX22bxxwLKxz3XFQ5Vbj8cPAkAoBzNec6Yu2cljXF+GHB+69obdP3aFn3+KWIMAHDp5jpMedqEpKfN7EFJZz4Z2d1/ZUmmAsrM229eqz/8wrPafWxEV6xoih4HAFBGil3a4ouS/qOkhyRtnfUFQNIP3bhaKZMeeOJg9CgAgDJT1J4xd//oUg8ClLMVTXW688ou/d8nDunX3nyVUimLHgkAUCbmWtri82b2Q2ZWc57HNpjZ+83sp5duPKB8vOOWtTo4NK5HXjgePQoAoIzMdZjyZyXdKWmnmT1WWOz162b2gqS/lbTV3T+y5FMCZeD7r12llvoafeLRF6NHAQCUkYsepnT3I5J+Q9JvmFmP8ou+jkva5e5jC31xM7tH0n+XVCXpQ+7+xwt9TiBKXU2V3nVrtz728D71j0yqq6k2eiQAQBko6gR+M7tc0hF3f9jdn5SUK8TZvJlZlaQPSnqL8qv632dm1yzkOYFoP/bq9ZrOuv5xy4HoUQAAZaLYqyn/SVJu1u1c4b6FeJWk3e6+192nJH1K0tsW+JxAqI1djbp9Q4c+/sh+TWdzc/8CAEDufuYrl1v4V3YRvmayuQV/FavYdcaqC8F0+k2bMrP0pb7Z51grafbugz5Jr17gcwLhfvq1l+tnP7ZFX3r6sN5209rocYB5y+Vc07mcZrKumaxrKpvTTOH2dDan6cL3mcJfXGdv538+u11+m2zOlfOzf9nlf9Z57nPl/EL3n75v1uN+9i/gnJ/9XS98lyR3yeXK5fLf3c/el/9+9ndcp+Pg7OM5z9+nC26X/znn5zy/e+GxwusX7j/r7Hwvvafw86wH/Mx9F3989o25nmv28/msey80jy5h29Ovd77XwksVG2P9Znavu/+zJJnZ2yQt9GOQznft/8v+bzKzzZI2S9L69esX+JLA0nvT1Su0sSujv/3Xvbr3xjUyY5kLLA531+RMTmNTWY1NzWhsKqvRyZnC7fx9o5P57xPTWU3O5DQ5k8v/PJ3T5Mw5980U7pvOaWLm9Db5+6ZmcsoF/MVpJlWZKWWmVKrwc8pUlbKzP1v+9kseP31f4fdSZjIzWeE5Tafvk0z576mUZErlfy48ppdsJ5mZ8ivVFH5n1u/nv+zMdvnXOOe+M88z+zXOvlb+1tn/7WfvO3vj/Nue/8+V8z3H+X7/ZfcXbrzkWeeY54LPNce2szcu5jnmazH+6F3oU/zKnxS3XbEx9vOSPm5mf6X8bAck/cS8JjurT9K6Wbe7JR06dyN3v1/S/ZLU29tLU6PkpVKmzXdt0G9+5ml96/kB3bWpK3oklAB31/h0VsPjMxqemNbw+LROjk8Xfp7RcOHnk+PTZ7eZmNbYZFajUzMam8xqbDqr7CUUUsryF5bUVqdUW12l2prUmZ/ralKqr6lSa31N4f78dnU1VUpXp5SuSqm6ylRTlVJ1Kv+9pspU/ZLbp7cxVadSF9gm/70qZapOnQ2nlwaWzsQU//GCSlLsxxQVu+jrHkm3mVmjJHP3kXlPdtZjkq4sXBxwUNKPSvqxRXheINzbb16rP3/wef3F13bpzis7+QumQp2anNHAyKSOj07p+KlJDY5OFX6e0uDo5Kyf819Tc5xDUl9Tpeb6ajXX1ailvkZdjbXKdFQrk65WQ22VGtJVakhXK5OuUkNt4f50/v5MbfWZ7/XpKtXXVKmauAHKQlExZma1kt4pqUdS9el/ud39/fN9YXefMbNfkvQV5Ze2+Ii7b5/v8wGlpLa6Sr/8piv0vgee0TeeO6Y3Xr0yeiRcgmzO1T8yqSPDEzpyckJHhyd0ZHhCR0/mv5/+eXQqe97fb0hXqaMxrfZMrVa11OnaNc1qb0yrrSGt5rqaM8HVXJ+Prua6ajXV1ShdXew1VQAqSbGHKT8n6aTyn0c5uVgv7u5fkvSlxXo+oJT8SO86/e2/7tWffWWXXr9pBR+RVELcXYOjUzpwYlwvDo7pwOCY+k6M6cDguA6cGNOhoXFNZ196OLA6ZVrZXKeVzbW6elWTXrepSyub69TZWKuOxrQ6Mml1NNaqI5NWXU1V0P8yAOWo2Bjrdvd7lnQSoMLUVKX0q2/epPf+w5P69ON9+pHedXP/EhbVqckZ7e0/pT39p7S3f/TM9xcHxzR2zl6tjkxa3e0Nun5ti37w+tVa21avVc11hQCrU0cmTVADWBLFxth3zex6d396SacBKszbblqjjz28T//1X3bqnutWqbnuZR/zikVwanJGzx0Z0c4jw9p5eER7CgF2dPjsjvyqlOmy9gZd3pnR7Rs7tK6tQevbG7SuvUHdbfXK1Bb7xyEALK5i//R5raSfLHwm5aTyV1S6u9+wZJMBFcDM9Af3Xqd7P/htfeCru/T7914bPVJZc3e9ODimHYeHteNwIb6OjGj/8bOfztZUW60rVjbqtVd0aeOKjDZ2NWpjV0br2zOckwWgJBUbY29Z0imACnZ9d4vefdtl+ujD+/QD16/Wqy5vjx6pbBwbntBTfSf11IEhPdU3pKcPntTQ2LSk/BpCl3dmdN2aFr3rlm69YnWzrl7dpLWt9VxBCKCsFBtjrO8FLMBv3nO1vvHcMf36p5/Sl99zpxrSHBI718R0Vtv6TmrL/kE9+eKQtvWd1JHhCUn5Q4ybVjbpnmtX6YbuVl27plmbVjapPs2J8gDKX7F/I3xR+SAzSXWSLpf0nCSOuQBFyNRW60/fdaPu+7tH9DuffVp//m9vSvzem5Nj09qyf1CP7TuhLfsGta3v5Jl1uC7vzOjVG9p1Q3erbuxu0bVrWggvABWr2EVfr59928xukfRzSzIRUKFu29Ch/3D3Jn3gwV3q7WnXj992WfRIy+rE6JQe3ntc390zoMdeOKHnjubXjq5Oma7vbtFP3tGjV/a069bL2tSeWehH3wJA+ZjXsRJ3f9zMXrnYwwCV7pfecIW27j+h3//n7bqso0F3Xlm5H5U0NjWj770wqO/uOa7v7B7Qs4eH5S5l0lW6taddb71htXp72nXTulb2egFINPMiPkLdzH511s2UpFsltbv79y/VYOfT29vrW7ZsWc6XBBbd8MS0fuRvHlbfiXF98mdv0/XdLdEjLYrpbE5PHhjSd3YP6Lu7j+uJAyc0nXXVVJluXt+mOzZ26o4rOnTjulbVVHFVI4DKZ2Zb3b13zu0uFmNm9vfu/m4zG5L054W7ZyTtk/QZd59YhFmLRoyhUhw+Oa53/fXDGpmY1t//zKt147rW6JHmZd/AqL71fL8een5AD+85rlOTMzKTrlvTotds7NBrrujUK3vauGABQCItVow9q/yyFp+X9PpzH3f3wQXMeMmIMVSSvhNjuu/vHtGJ0Wn95X036w1Xr4geaU7DE9N6eM9xPbSrX996fkAvDubX9+puq9ddm7p05xWdun1jh1obOOcLABYrxn5F0i8of/XkodkPKb/o64aFDnopiDFUmiMnJ/QzH31MOw4P69e+7yr9/Os2qqqEPnInm3Nt6xvSt54f0EO7+vXEgSFlc65Mukq3b+zIB9iVXerpaEj81aEAcK5FibFZT/bX7v4LizLZAhBjqERjUzP69U9v0xe3HdZtG9r1X95xgy7vzITN03diTN9+fkDfen5A3949oJPj0zKTrl/bojuv7NSdV3bplvVtrGYPAHNY1BgrFcQYKpW765+29un9n39WUzM5/dQdPdp81wZ1NNYu+WsfGBzTI3uP69EXBvXI3uPqOzEuSVrZXKs7r+zSXZu6dMfGjmWZBQAqCTEGlKFjIxP64y/v1ANPHFRddZXefvMa/XDvOt28rnVRDgNOTGf17OFhPV34iKFHXxjUwaF8fLU11OhVl7frtg0des3GTm1a2cihRwBYAGIMKGO7j53S/Q/t0eefOqzx6axWt9Tp9g0dur67RTd0t2pjV0bNdTVKnef8slzOdWJsSsdGJrVvYFR7B0a1t39U2w+d1PPHTimby/8739mY1it78vH16g3t2rSi6bzPBwCYH2IMqADDE9N6cPtRPfjsUW198YT6RybPPFaVMrVn0mqsrVY258rmXJMzOQ2OTip3zr/WK5pq9YrVzbqhu0XXrW3R9WtbtLqljj1fALCEio0xFv8BSlhzXY3eeWu33nlrt6T81Zfb+oZ04MS4jp+a1ODolEansqoyKZUy1Van1JGpVWdjWp1NterpyKinM6PGWv5VB4BSxZ/QQBlZ1VKnVS2roscAACwirk0HAAAIRIwBAAAEIsYAAAACEWMAAACBiDEAAIBAxBgAAEAgYgwAACAQMQYAABCIGAMAAAhEjAEAAAQixgAAAAIRYwAAAIGIMQAAgEDEGAAAQCBiDAAAIBAxBgAAEIgYAwAACESMAQAABCLGAAAAAhFjAAAAgYgxAACAQMQYAABAIGIMAAAgEDEGAAAQiBgDAAAIRIwBAAAEIsYAAAACEWMAAACBiDEAAIBAxBgAAEAgYgwAACAQMQYAABCIGAMAAAhEjAEAAAQixgAAAAIRYwAAAIFCYszMftjMtptZzsx6I2YAAAAoBVF7xp6R9A5JDwW9PgAAQEmojnhRd98hSWYW8fIAAAAlo+TPGTOzzWa2xcy29Pf3R48DAACwqJZsz5iZfU3SqvM89D53/1yxz+Pu90u6X5J6e3t9kcYDAAAoCUsWY+5+91I9NwAAQKUo+cOUAAAAlSxqaYt/Y2Z9km6X9EUz+0rEHAAAANGirqZ8QNIDEa8NAABQSjhMCQAAEIgYAwAACESMAQAABCLGAAAAAhFjAAAAgYgxAACAQMQYAABAIGIMAAAgEDEGAAAQiBgDAAAIRIwBAAAEIsYAAAACEWMAAACBiDEAAIBAxBgAAEAgYgwAACAQMQYAABCIGAMAAAhEjAEAAAQixgAAAAIRYwAAAIGIMQAAgEDEGAAAQCBiDAAAIBAxBgAAEIgYAwAACESMAQAABCLGAAAAAhFjAAAAgYgxAACAQMQYAABAIGIMAAAgEDEGAAAQiBgDAAAIRIwBAAAEIsYAAAACEWMAAACBiDEAAIBAxBgAAEAgYgwAACAQMQYAABCIGAMAAAhEjAEAAAQixgAAAAIRYwAAAIGIMQAAgEDEGAAAQCBiDAAAIBAxBgAAEIgYAwAACESMAQAABCLGAAAAAhFjAAAAgYgxAACAQCExZmZ/amY7zWybmT1gZq0RcwAAAESL2jP2oKTr3P0GSbsk/XbQHAAAAKFCYszdv+ruM4Wbj0jqjpgDAAAgWimcM/bTkr4cPQQAAECE6qV6YjP7mqRV53nofe7+ucI275M0I+njF3mezZI2S9L69euXYFIAAIA4SxZj7n73xR43s38n6a2S3uTufpHnuV/S/ZLU29t7we0AAADK0ZLF2MWY2T2SflPS69x9LGIGAACAUhB1zthfSWqS9KCZPWlmfxM0BwAAQKiQPWPufkXE6wIAAJSaUriaEgAAILGIMQAAgEDEGAAAQCBiDAAAIBAxBgAAEIgYAwAACESMAQAABCLGAAAAAhFjAAAAgYgxAACAQMQYAABAIGIMAAAgEDEGAAAQiBgDAAAIRIwBAAAEIsYAAAACEWMAAACBiDEAAIBAxBgAAEAgYgwAACAQMQYAABCIGAMAAAhEjAEAAAQixgAAAAIRYwAAAIGIMQAAgEDEGAAAQCBiDAAAIBAxBgAAEIgYAwAACESMAQAABCLGAAAAAhFjAAAAgYgxAACAQMQYAABAIGIMAAAgEDEGAAAQiBgDAAAIRIwBAAAEIsYAAAACEWMAAACBiDEAAIBAxBgAAEAgYgwAACAQMQYAABCIGAMAAAhEjAEAAAQixgAAAAIRYwAAAIGIMQAAgEDEGAAAQCBiDAAAIBAxBgAAEIgYAwAACESMAQAABAqJMTP7QzPbZmZPmtlXzWxNxBwAAADRovaM/am73+DuN0n6gqT/FDQHAABAqJAYc/fhWTczkjxiDgAAgGjVUS9sZn8k6ScknZT0hotst1nS5sLNSTN7ZhnGq1SdkgaihyhTvHcLw/u3MLx/C8P7N3+8dwtzVTEbmfvS7JQys69JWnWeh97n7p+btd1vS6pz998r4jm3uHvvIo6ZKLx/88d7tzC8fwvD+7cwvH/zx3u3MMW+f0u2Z8zd7y5y009I+qKkOWMMAACg0kRdTXnlrJv3StoZMQcAAEC0qHPG/tjMrpKUk7Rf0s8X+Xv3L91IicD7N3+8dwvD+7cwvH8Lw/s3f7x3C1PU+7dk54wBAABgbqzADwAAEIgYAwAACFQWMWZm95jZc2a228x+K3qecmNmHzGzY6zRdunMbJ2ZfcPMdpjZdjN7T/RM5cTM6szse2b2VOH9+4PomcqNmVWZ2RNm9oXoWcqNme0zs6cLH723JXqecmNmrWb2aTPbWfgz8PbomcqFmV1V+Ofu9Newmb33gtuX+jljZlYlaZekN0vqk/SYpPvc/dnQwcqImd0l6ZSkj7n7ddHzlBMzWy1ptbs/bmZNkrZKejv//BXHzExSxt1PmVmNpG9Leo+7PxI8Wtkws1+V1Cup2d3fGj1POTGzfZJ63Z1FS+fBzD4q6Vvu/iEzS0tqcPeh4LHKTqFjDkp6tbvvP9825bBn7FWSdrv7XnefkvQpSW8LnqmsuPtDkgaj5yhH7n7Y3R8v/DwiaYektbFTlQ/PO1W4WVP4Ku3/AiwhZtYt6QclfSh6FiSLmTVLukvShyXJ3acIsXl7k6Q9FwoxqTxibK2kA7Nu94m/DBHAzHok3Szp0eBRykrhMNuTko5JetDdef+K9xeSfkP5ZYBw6VzSV81sa+Gj9VC8DZL6Jf2vwmHyD5lZJnqoMvWjkj55sQ3KIcbsPPfxX9ZYVmbWKOkzkt57zgfdYw7unnX3myR1S3qVmXGovAhm9lZJx9x9a/QsZewOd79F0lsk/WLhlA0Up1rSLZL+2t1vljQqiXO2L1Hh8O69kv7pYtuVQ4z1SVo363a3pENBsyCBCuc6fUbSx939s9HzlKvCIY5vSrondpKycYekewvnPX1K0hvN7P/EjlRe3P1Q4fsxSQ8of9oLitMnqW/WnuxPKx9nuDRvkfS4ux+92EblEGOPSbrSzC4vFOaPSvrn4JmQEIUT0D8saYe7fyB6nnJjZl1m1lr4uV7S3eLjz4ri7r/t7t3u3qP8n3tfd/cfDx6rbJhZpnDRjQqH175PEleUF8ndj0g6UPi0HCl/3hMXLl26+zTHIUop7uOQiubuM2b2S5K+IqlK0kfcfXvwWGXFzD4p6fWSOs2sT9LvufuHY6cqG3dIerekpwvnPUnS77j7l+JGKiurJX20cDVRStI/ujtLNGA5rJT0QP6/p1Qt6RPu/i+xI5WdX5b08cKOkL2Sfip4nrJiZg3KrwTxc3NuW+pLWwAAAFSycjhMCQAAULGIMQAAgEDEGAAAQCBiDAAAIBAxBgAAEIgYAwAACESMAQAABCLGAGAWM/tLM3vczF4ZPQuAZCDGAKCg8LE5K5RfMfutweMASAhiDEDimFmPmY3P+ogrSZK7jyr/EU7flPQ/CtvWm9mTZjZlZp3LPiyAikeMAah4lnfun3d73P2mc7brkNQgaURSVpLcfbyw3aFlGBVAAhFjACpSYe/XDjP7n5Iel7SuiF/7XUl/Jmm7pGuWcj4AOI0YA1DJrpL0MXe/2d33X2xDM+uR9BpJ/yBph6Rrl348ACDGAFS2/e7+SJHb/mdJ73d3FzEGYBlVRw8AAEtotJiNzOwmSe+Q9Foz+6CkOklPL+FcAHAGMQYA0p9I+iF3/3+SZGYrJT0ROxKApOAwJYBEM7M3SsqcDjFJcvejkjJm1h43GYCkYM8YgIrk7vskXVfEdl+X9PXz3N+yBGMBwMuwZwxAEmUltZy76Ov5nF70VVKNpNwSzwUggSx/4RAAAAAisGcMAAAgEDEGAAAQiBgDAAAIRIwBAAAEIsYAAAACEWMAAACBiDEAAIBA/x8XHUTgZ5J39QAAAABJRU5ErkJggg==)
%% Cell type:markdown id:forced-listing tags:
Additionally it is a good idea to define limits for the parameters. 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.
%% Cell type:code id:sunset-beginning tags:
``` python
V.parameters.D0.min_val = 0
V.parameters.D0.max_val = 5
V.parameters.r0.min_val = 1
V.parameters.r0.max_val = 2.5
V.parameters.D0.max_val = 2
V.parameters.r0.min_val = 1.5
V.parameters.r0.max_val = 3.0
V.parameters.S.min_val = 1.1
V.parameters.S.max_val = 10.0
V.parameters.delta.min_val = -1
V.parameters.delta.max_val = 1
V.parameters.beta.min_val = 0.1
V.parameters.beta.max_val = 10
```
%% Cell type:markdown id:rough-purchase tags:
%% Cell type:markdown id:raised-kansas tags:
Additionally a screening function needs to be defined for the morse potential
A screening function needs to be defined for the morse potential
%% Cell type:code id:practical-details tags:
%% Cell type:code id:republican-rabbit tags:
``` python
V.screening = job.factories.functions.exp_A_screening(identifier="V_cutoff", cutoff=7)
```
%% Cell type:markdown id:external-ready tags:
%% Cell type:markdown id:diverse-smith tags:
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 id:departmental-dynamics tags:
%% Cell type:code id:described-blond tags:
``` python
rho = job.factories.functions.spline(identifier="rho_CuCu", cutoff=7, derivative_left=-1)
```
%% Cell type:markdown id:latter-wright tags:
%% Cell type:markdown id:solved-debut tags:
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.
Too few nodepoints mean low flexibilty, too many lead to overfitting. This requires some testing to find an optimal choice.
%% Cell type:code id:double-engineering tags:
%% Cell type:code id:hazardous-lying tags:
``` python
rho_nodes = np.linspace(0.5, 7.0, 7).round(2)
rho_nodes
```
%%%% Output: execute_result
array([0.5 , 1.58, 2.67, 3.75, 4.83, 5.92, 7. ])
%% Cell type:markdown id:handled-housing tags:
%% Cell type:markdown id:blocked-volleyball tags:
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 id:developing-apache tags:
%% Cell type:code id:geographic-legislation tags:
``` python
decaying_exp = lambda r: np.exp(-r)
rho_initial = decaying_exp(rho_nodes)
```
%% Cell type:markdown id:separated-journal tags:
%% Cell type:markdown id:outdoor-seven tags:
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.
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. 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.
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 id:colored-discount tags:
%% Cell type:code id:collectible-viewer tags:
``` python
rho_mins = np.zeros((len(rho_nodes)))
rho_maxs = 3*rho_initial.round(6)
rho.parameters.create_from_arrays(rho_nodes, rho_initial, min_vals=rho_mins, max_vals=rho_maxs)
```
%% Cell type:raw id:preliminary-peace tags:
%% Cell type:raw id:bottom-potato tags:
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 id:lasting-reply tags:
%% Cell type:code id:narrow-hearts tags:
``` python
rho.parameters["node_7.0"].start_val = 0
rho.parameters["node_7.0"].enabled = False
```
%% Cell type:markdown id:saving-vegetable tags:
%% Cell type:markdown id:miniature-department tags:
$-\sqrt(\rho)$ can be used as initial guess for the embedding energy, which is taken from second moment approximation tight binding.
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.
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.
Everything else is similar to the electron density.
%% Cell type:code id:overall-flower tags:
%% Cell type:code id:younger-freight tags:
``` python
F = job.factories.functions.spline(identifier="F_CuCu", cutoff=5)
F_nodes = np.linspace(0.0, 5.0, 7).round(2) #9 is worse 11 is worse 7 is best
F_nodes = np.linspace(0.0, 5.0, 7).round(2)
F_init = -np.sqrt(F_nodes)
F_maxs = np.zeros(len(F_nodes))
F_mins = -np.ones(len(F_nodes))*5
F.parameters.create_from_arrays(F_nodes, F_init, F_mins, F_maxs)
F.parameters["node_0.0"].enabled=False
F.parameters["node_0.0"].start_val = 0
F.parameters["node_5.0"].enabled=False
F.parameters["node_5.0"].start_val = 0
```
%% Cell type:markdown id:neutral-gateway tags:
%% Cell type:markdown id:alert-stake tags:
The functions have to be assigned to the potential
%% Cell type:code id:major-capacity tags:
%% Cell type:code id:banned-allocation tags:
``` python
job.potential.pair_interactions[V.identifier] = V
job.potential.electron_densities[rho.identifier] = rho
job.potential.embedding_energies[F.identifier] = F
```
%% Cell type:markdown id:unusual-retirement tags:
%% Cell type:markdown id:positive-frost tags:
### Define fitting procedure
Finally a few parameters need to be set that influence the fitting process.
The minimization can be done with different algorithms. Atomicrex itself implements the BFGS algorithm. Additionally the algorithms from the nlopt library can be used.
Finally atom types and fit algorithm have to be assigned.
Atomicrex itself implements the BFGS algorithm. Additionally the algorithms from the nlopt library can be used.
%% Cell type:code id:productive-spare tags:
%% Cell type:code id:existing-average tags:
``` python
## Define the atom types of the potential
job.input.atom_types.Cu = None
##
## Limited number of iterations for the workshop
job.input.fit_algorithm = job.factories.algorithms.ar_lbfgs(max_iter=500)
job.run()
job.output
```
%%%% Output: stream
The job PotentialDF1 was saved and received the ID: 622
The job PotentialDF1 was saved and received the ID: 848
%%%% Output: execute_result
Output({'error': None, 'residual': array([1.39371e+03, 1.39371e+03, 1.39371e+03, ..., 1.52231e-01,
1.52231e-01, 1.52231e-01]), 'iterations': array([ 1, 2, 3, ..., 1998, 1999, 2000], dtype=uint32)})
%% Cell type:markdown id:regulation-foster tags:
### Plot the results
%% Cell type:markdown id:black-assets tags:
Plot the resiudal over steps to see how the calculation converges
%% Cell type:code id:partial-burner tags:
``` python
plt.plot(job.output.iterations, job.output.residual)
plt.yscale("log")
```
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAZWElEQVR4nO3de5AcZ3nv8e/TM3vR/brYji5IsoVBJ2DZLDZXkxSEyOQYEThwrLjqmIMiYxIfkspJKqYIqZyqk6QcV6gUhYJLCY5J4mNFmASkOgKTEIODsYPWxsIWQmglhL2Wba2k+CJLq92defLH9OzOjnaWnZ3py3b/PlVbu9PT8/ajntFv3nm7521zd0REJPuCpAsQEZF4KPBFRHJCgS8ikhMKfBGRnFDgi4jkRDHpAqayfPlyX7NmTdJliIjMGo8++uhJd++Z7L5UB/6aNWvo6+tLugwRkVnDzH7a6L5UDumY2fVmtuPFF19MuhQRkcxIZeC7+x53v3nRokVJlyIikhmpDHwREWk/Bb6ISE4o8EVEckKBLyKSEwp8EZGcyGTg/81DP2HP/uNJlyEikiqp/uLVTP2fPT+ko2Bcf8XPJV2KiEhqZLKH/7F3rsPMki5DRCRVMhn4hoEu5CUiMkE2A1+dexGRC8QW+Gb2OjO708zuM7OPR709VxdfRGSClgLfzO4ysxNm9mTd8k1mdsjM+s3sNgB3P+jutwAfBnpb2e7PrAvQtdlFRCZqtYd/N7CpdoGZFYDtwHXABmCLmW0I73sf8B3gmy1ud0qmIXwRkQu0FPju/iBwum7x1UC/ux9192FgJ7A5XH+3u78VuLFRm2Z2s5n1mVnf4ODgjOoyDFcXX0RkgijOw18BPF1zewC4xsx+AfgA0AXsbfRgd98B7ADo7e2dUWrroK2IyIWiCPzJ4tbd/VvAtyLY3qTUvxcRmSiKs3QGgFU1t1cCTc1z0OoVr3TQVkTkQlEE/j5gvZmtNbNO4AZgdzMNtHzFK43piIhcoNXTMu8FHgYuN7MBM9vq7qPArcD9wEFgl7sfaLLdlnv4IiIyUUtj+O6+pcHyvUxxYHYa7e4B9vT29m6baRthO5pTR0QklOmpFTSOLyIyLpWB3/qQTiXxlfciIuNSGfitHrQd7+Er8kVEqlIZ+K3SqL2IyIUyGfhV6t+LiIxLZeC3PIavg7YiIhdIZeC3PoZfPWirxBcRqUpl4LeLevgiIuNSGfitDumIiMiFUhn47TotU0RExqUy8Fs19sUrDemIiIzJZuBXz9LRQVsRkTHZDPykCxARSaFUBn67DtpqSEdEZFwqA79tc+m0sSYRkdkulYHfqvGDtop8EZGqbAa+evgiIhfIZOCLiMiFMh34GtERERmXycA3jemIiFwglYHf+iUOK/TFKxGRcakMfM2lIyLSfqkM/HbRGL6IyLhMBv74kI6IiFRlM/BNX7wSEamX0cCv/Fbci4iMy2bgJ12AiEgKZTLwqzSiIyIyLpWB3/L0yNUxfA3qiIiMSWXgt3we/lhDbStJRGTWS2Xgt0pfvBIRuVAmA79KHXwRkXGZDPzxC6AkXIiISIpkM/DHzsNX4ouIVGUz8MPf6uGLiIzLZuCHif/S0EiyhYiIpEgmA//SnvkAfPG7P024EhGR9Mhk4PeuWUpnMeDe7z3FueFS0uWIiKRCJgMf4OPvvBSABw6dSLgSEZF0yGzg/9o1qwH4j7PDCVciIpIOsQW+mb3fzP7KzL5qZu+JensLuzsA+M7hk1FvSkRkVmgp8M3sLjM7YWZP1i3fZGaHzKzfzG4DcPevuPs24CPAf29lu9Mxp7MAQKB5FkREgNZ7+HcDm2oXmFkB2A5cB2wAtpjZhppV/iC8P3I/v2IhQyM6aCsiAi0Gvrs/CJyuW3w10O/uR919GNgJbLaK24Gvuftjjdo0s5vNrM/M+gYHB1spj+5igXMKfBERIJox/BXA0zW3B8Jl/wt4N/DfzOyWRg929x3u3uvuvT09PS0V0t1R4LtHTrXUhohIVhQjaHOyQXN3988Cn41gew2Vyk4x0Bi+iAhE08MfAFbV3F4JHG+mgZaveBW66tWLNX2aiEgoisDfB6w3s7Vm1gncAOxupoFWr3hV1VkoUCo7pbJiX0Sk1dMy7wUeBi43swEz2+ruo8CtwP3AQWCXux9ost229PA7ipXhnJFSuaV2RESyoKUxfHff0mD5XmBvC+3uAfb09vZum2kbAJ2FyvvZ+dEy3R2FVpoSEZn1Mju1AkBnsfLPe+7FoYQrERFJXioDv11DOj3zuwB45oWz7ShLRGRWS2Xgt+ug7esuWQhoXnwREUhp4LfLqqVzedtly/j2jwf5yvefSbocEZFEpTLw2zWkUwiM7b92FQDf/JHmxReRfEtl4LdrSAdg8dxOrli5iIPPvtSGykREZq9UBn67LejuoP/EGZ4+rYO3IpJfuQj8bdeuA+A37mk4SaeISOalMvDbNYZf9c7X9HBpzzyeeOZFDhxvT5siIrNNKgO/nWP4Vbd/8A2ALnkoIvmVysCPwpWrlwDogigiklu5CfxCYHQUjKERTaQmIvmUm8CHymRqPzl5JukyREQSkcrAb/dB26pXhkuaSE1EciuVgR/FQVuAjasW01XUNMkikk+pDPyozOsqUHZd/UpE8ilXgR+YKfBFJLdyFfhmhi5vKyJ5lavADwxcPXwRyalUBn5UZ+kE6uGLSI6lMvCjOksnMDSGLyK5lcrAj4rG8EUkz3IV+BrDF5E8y1ng67RMEcmvHAZ+0lWIiCQjV4FvOmgrIjmWq8APzFDei0he5Szw1cMXkfxKZeBH+8UrBb6I5FMqAz+qL16ZGWVd8EpEciqVgR8VnYcvInmWs8A3RnRepojkVK4Cf6RcZvDl8wyPalxHRPInV4G/ZG4nAOdGSglXIiISv1wF/qolcwAoaVhHRHIoV4FfCAxQ4ItIPuUq8IMw8HUuvojkUa4Cv2Dq4YtIfuUq8AMN6YhIjuUq8Ks9fA3piEgexRb4ZrbOzL5gZvfFtc16OmgrInnWUuCb2V1mdsLMnqxbvsnMDplZv5ndBuDuR919ayvba5UO2opInrXaw78b2FS7wMwKwHbgOmADsMXMNrS4nbYohoE/qh6+iORQS4Hv7g8Cp+sWXw30hz36YWAnsHm6bZrZzWbWZ2Z9g4ODrZR3gUBn6YhIjkUxhr8CeLrm9gCwwsyWmdmdwJVm9slGD3b3He7e6+69PT09bS2sOoZ/7OTZtrYrIjIbFCNo0yZZ5u5+Crglgu1N25plcwHYd+w0v/KGS5IsRUQkdlH08AeAVTW3VwLHm2kgqiterb9oARct7OIbB55ra7siIrNBFIG/D1hvZmvNrBO4AdjdTANRXfEKYH5XkfndUXywERFJt1ZPy7wXeBi43MwGzGyru48CtwL3AweBXe5+oMl2I+nhA7z24oU6aCsiudRSV9fdtzRYvhfY20K7e4A9vb2922baRiNBYCjvRSSPcjW1AkDBdFqmiORTKgM/yiGdIDAFvojkUioDP8qDtsXANLWCiORSKgM/SoXANLWCiORSKgM/0iEdM8oKfBHJoVQGfpRDOoXAKGlIR0RyKJWBH6XAdNBWRPIpd4FfDDSkIyL5lMrAj3IMvxAYrwyXGBoptb1tEZE0S2XgRzmGv3x+FwCf+9f+trctIpJmqQz8KL3/yhV0FgM+90A/h59/OelyRERik7vA71nQxd9vvQaAP/xqU3O6iYjMarkLfIA3rVnCgq4iDx89Rf+JM7hO0xSRHEhl4Ed50DZsn+03XgXAuz/zba694wG++vgzkWxLRCQtUhn4UR60rbr2NT184aZePvXe19FRCPitnY/zwKETkW1PRCRpqQz8uLzrdRex7dp13P2RqwH4vS/9IOGKRESik+vAr1q9bC43XrOak2fOc/LM+aTLERGJhAI/dMWqxQD8WKdqikhGKfBDG8PAf+TIqWQLERGJSCoDP+qzdCZzac98AM2VLyKZlcrAj+MsnXqFwOgsBpo6WUQyK5WBnxTNpCkiWabAr1EwXf5QRLJLgV8jUA9fRDJMgV9Dlz8UkSxT4NcoBEapnHQVIiLRUODXKJhRKivxRSSbFPg11MMXkSxLZeAn8cUrgCCAssbwRSSjUhn4SXzxCqpDOgp8EcmmVAZ+UoLA2L3/eNJliIhEQoFfY9GcDgBGNJAvIhmkwK/x/o0rAHj69NmEKxERaT8Ffo3Fcys9/GOnXkm4EhGR9lPg17j84gUADI1oSEdEskeBX6O7WABgaKSUcCUiIu2nwK/R3VEJ/GMnNaQjItmjwK+xeG4H3R0Bn/3Xfs4OjyZdjohIWynwa3R3FNj69rUAfOvQYMLViIi0lwK/zofeuArQOL6IZE8xrg2Z2TzgL4Fh4Fvufk9c225GsWCALmYuItnTUg/fzO4ysxNm9mTd8k1mdsjM+s3stnDxB4D73H0b8L5WthulYlDZJaMlBb6IZEurQzp3A5tqF5hZAdgOXAdsALaY2QZgJfB0uFpqx0sKQaWHr3nxRSRrWgp8d38QOF23+Gqg392PuvswsBPYDAxQCf0pt2tmN5tZn5n1DQ7Gf+C0Q0M6IpJRURy0XcF4Tx4qQb8C+Efgg2b2eWBPowe7+w5373X33p6engjKm1q1h68hHRHJmigO2toky9zdXwH+57QaMLseuP6yyy5ra2HT0VEIx/DVwxeRjImihz8ArKq5vRJoapL5pC6AAuM9/EPPvRT7tkVEohRF4O8D1pvZWjPrBG4AdkewnUgUw8Af0ZCOiGRMq6dl3gs8DFxuZgNmttXdR4FbgfuBg8Audz/QZLuJXNM23DbrXzVf17YVkcxpaQzf3bc0WL4X2NtCu3uAPb29vdtm2kYrioVAY/gikjmaWmESxcAY1WUORSRjUhn4SQ7pQGV6BfXwRSRrUhn4SZ6lA9UevgJfRLIllYGftGIQMKqpFUQkY1IZ+BrSERFpv1QGftJDOmbG9596QVe9EpFMSWXgJ23d8nkAnDoznHAlIiLto8CfRO+aJQCc01WvRCRDUhn4SY/hz+koALrMoYhkSyoDP+kx/Grgf6f/ZCLbFxGJQioDP2mvX1l5o3lIgS8iGRLbRcxnkwXdHWxctZiH+k+x/YF+Vi+dy6I5HbxqYRevvXhh0uWJiMyIAr+BW955Kbd//Ufccf+hCctfv2IRuz72FuZ0FhKqTERkZsxTOA1wzRWvth0+fDjRWl4eGuGZF85xZmiUO799lH85+Dyrl85l18fewsWLuhOtTUSknpk96u69k96XxsCv6u3t9b6+vqTLmODGv36Eh/pPAfCZD1/BB65a+TMeISISn6kCXwdtm/T3W6/hK7/5NgB+Z9d+tj/Qr9M3RWRWUOA3yczYuGoxO29+M1euXswd9x/itZ/+Oju/91TSpYmITElDOi362hPP8vF7HgNg46rFvHrZXN566TI2b1xBd4cO7IpIvDSkE6HrXn8Jj336l/jQG1dSDIw9+4/z+19+gv/7/3+YdGkiIhOk8rTMmrN0ki5lWpbO6+SOD10BVKZjeO2nv07fsf+g79hpetcsTbg6EZEKDelE4He/tJ/7Hh0AKlfPmtdVpKsY0FEImN9V5OJF3XzkbWv4xctflXClIpI1Oi0zZuWyc/zFc3zjwPMMnjnP2fOjDJfKnB8p88rwKPcfeB6ARz75Lp3LLyJtNVXgp3JIZ7YLAmPlkrl89O1rJ73/kaOnuGHHI3xi5/fZ9bG3xFydiOSVDtom4M3rlrF0Xif7n34h6VJEJEcU+Am5dv1yzo+WOf2KrqolIvFQ4CfkPf/lYgC+e0RTMItIPBT4CblqdeUyig8fOZVwJSKSFwr8hFy0sAuAe/79KQ4993LC1YhIHqQy8JO+pm0czIw/+dXX090R8Mt/8SC//sV9jJbKSZclIhmm8/ATduLlIf73rv382+GTdHcE3HXTm3jrZcuTLktEZil98SrlymXn9+77AV9+rPLt3J9fsZBl87roKAR0Fm3sG7obVy1mxZI5vGHlYuZ1FjCzhCsXkbRR4M8SZ4dH+et/+wn7jp3mpaFRRkbLjJQqPyfPDHPm/OiE9YuBUSwYxSBgybwO3vTqpfzc4jl0FgOqbwXV94T6N4ex5eGa47envn/88dZg/br765bT8HHTq6P2V2CG2fhvM6tbXlkzCO+r/J5svcptq2tvrORJ9ttk9zZ+jDVYXrv+5O1M2NI01p9um43/Te3995hBVzGgq1igI+y8FANTZyVC+qbtLDG3s8gn3rV+0vtGS2WeOn2Wn54+yw+Pv8T50TKjpTKlsjNScvoHz/DQkZOcePk8KX4PFwFgQVeRngVdY2+wtW/aQd2b+NhtxpcXAqO7oxA+lp/xxh++/UzSTvVNvhi2110M6O4s0F0sVG53BMzpqPzd1RHQ3VEYu93dEdBdLDCns0BXMZgVb2IK/FmiWAhY1zOfdT3zp5x0zd0plSuJ72PLwt/hkvo3hEb3jz9+8vaY4eMcr3t8k3X4+Pplr6xddse9sk7t7fHllXbL4TqVXVS9PfHx5Zr9M1Zr3X6rf0+t/aQ84b7ptNWg3UZteoMNNG5n8hqmu+3G6//sf0/ZnaHRMsNhB2WkVGa45Jx+5TwvnB3Bw22Vy+PPZfU5qH8Oa38Pj5Z5aWiEUjl8zuue38rtyV4H1W1M3N5o2RkaKTE0MvMTJ7qK428IncXgwk+3XPjJmEk+RVXX2fuJd9BZbO95NQr8jDGrDPOISPPcnfOj5bHwHxopcW6kNH57tMTQcKnye8L9Zc6H650bKTE8WnnjmOyNt74DM2G9mgcEEfw3VuCLiITMwqGdjF6tLpXn4YuISPsp8EVEckKBLyKSEwp8EZGciC3wzWydmX3BzO6La5siIjJuWoFvZneZ2Qkze7Ju+SYzO2Rm/WZ221RtuPtRd9/aSrEiIjJz0z0t827gc8DfVheYWQHYDvwSMADsM7PdQAH407rHf9TdT7RcrYiIzNi0At/dHzSzNXWLrwb63f0ogJntBDa7+58C/3WmBZnZzcDNAKtXr55pMyIiUqeVL16tAJ6uuT0AXNNoZTNbBvwxcKWZfTJ8Y7iAu+8AdoSPGTSzn86wvuVAGq8fqLqao7qao7qak8W6Xt3ojlYCf7Iv/jactsvdTwG3NLMBd+9ptqgqM+trNGNcklRXc1RXc1RXc/JWVytn6QwAq2purwSOt1aOiIhEpZXA3wesN7O1ZtYJ3ADsbk9ZIiLSbtM9LfNe4GHgcjMbMLOt7j4K3ArcDxwEdrn7gehKbdqOpAtoQHU1R3U1R3U1J1d1pfqKVyIi0j6aWkFEJCcU+CIiOZG5wG9muocItr3KzB4ws4NmdsDMfitc/kdm9oyZPR7+vLfmMZ8Maz1kZr8cYW3HzOyJcPt94bKlZvbPZnY4/L0kzrrM7PKaffK4mb1kZr+dxP6abPqQmewfM3tjuJ/7zeyz1uKFThvUdYeZ/cjMfmBm/2Rmi8Pla8zsXM1+uzOquqaorennLqZ99g81NR0zs8fD5bHssymyId7XWOX6j9n4oTKtwxFgHdAJ7Ac2xLj9S4Crwr8XAD8GNgB/BPzuJOtvCGvsAtaGtRciqu0YsLxu2Z8Bt4V/3wbcHndddc/dc1S+NBL7/gKuBa4Cnmxl/wDfA95C5XsqXwOui6Cu9wDF8O/ba+paU7teXTttrWuK2pp+7uLYZ3X3/znwh3HuMxpnQ6yvsaz18Meme3D3YWAnsDmujbv7s+7+WPj3y1TOXloxxUM2Azvd/by7/wTop/JviMtm4Ivh318E3p9gXe8Cjrj7VN+sjqwud38QOD3J9qa9f8zsEmChuz/slf+Zf1vzmLbV5e7f8MpZcgCPUPkOTENR1NWotikkus+qwt7wh4F7p2qj3XVNkQ2xvsayFviTTfcwVeBGxipzD10J/Hu46NbwI/hdNR/b4qzXgW+Y2aNWma8I4CJ3fxYqL0jgVQnUVXUDE/8TJr2/oPn9syL8O676AD5KpZdXtdbMvm9m3zazd4TL4q6rmecu7treATzv7odrlsW6z+qyIdbXWNYCv6npHiIrwmw+8GXgt939JeDzwKXARuBZKh8pId563+buVwHXAb9pZtdOsW6s+9EqX9x7H/ClcFEa9tdUGtUR9377FDAK3BMuehZY7e5XAr8D/D8zWxhzXc0+d3E/p1uY2LGIdZ9Nkg0NV22w/ZbqylrgJz7dg5l1UHlC73H3fwRw9+fdveTuZeCvGB+GiK1edz8e/j4B/FNYw/PhR8TqR9jqFNZx78frgMfc/fmwxsT3V6jZ/TPAxOGVyOozs5uozEp7Y/jRnvDj/6nw70epjPu+Js66ZvDcxbnPisAHgH+oqTe2fTZZNhDzayxrgZ/odA/h+OAXgIPu/pma5ZfUrParQPXsgd3ADWbWZWZrgfVUDsi0u655Zrag+jeVg35Phtu/KVztJuCrcdZVY0KvK+n9VaOp/RN+JH/ZzN4cvhb+R81j2sbMNgG/D7zP3c/WLO+xynUqMLN1YV1H46or3G5Tz12ctQHvBn7k7mNDInHts0bZQNyvsZkedU7rD/BeKkfAjwCfinnbb6fy8eoHwOPhz3uBvwOeCJfvBi6pecynwloP0YYzJxrUtY7KEf/9wIHqfgGWAd8EDoe/l8ZZV7iducApYFHNstj3F5U3nGeBESq9qK0z2T9AL5WQO0LlokEWQV39VMZ3q6+xO8N1Pxg+v/uBx4Dro6pritqafu7i2Gfh8ruBW+rWjWWf0TgbYn2NaWoFEZGcyNqQjoiINKDAFxHJCQW+iEhOKPBFRHJCgS8ikhMKfBGRnFDgi4jkxH8CQW0wXvV+inEAAAAASUVORK5CYII=)
%% Cell type:markdown id:appreciated-bracelet tags:
Finally it is a good idea to have a look at the final potential. This can reveal unphysical behavior
%% Cell type:code id:demographic-marijuana tags:
``` python
job.plot_final_potential()
```
%%%% Output: execute_result
(<Figure size 576x1296 with 3 Axes>,
array([[<AxesSubplot:title={'center':'Cu F'}, xlabel='$\\rho $ [a.u.]'>],
[<AxesSubplot:title={'center':'Cu rho_CuCu'}, xlabel='r [$\\AA$]'>],
[<AxesSubplot:title={'center':'Cu V_CuCu'}, xlabel='r [$\\AA$]'>]],
dtype=object))
%%%% Output: display_data
![](data:image/png;base64,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=)
%% Cell type:markdown id:permanent-convention tags:
### Run a lammps simulations with the fitted potential
%% Cell type:code id:inclusive-brick tags:
``` python
lmp = pr.create_job("Lammps", "template")
lmp.structure = pr.create_ase_bulk("Cu", cubic=True)
lmp.potential = job.lammps_potential
murn = lmp.create_job("Murnaghan", "PotentialTest")
murn.run()
```
%%%% Output: stream
The job PotentialTest was saved and received the ID: 849
The job strain_0_9 was saved and received the ID: 850
The job strain_0_92 was saved and received the ID: 851
The job strain_0_94 was saved and received the ID: 852
The job strain_0_96 was saved and received the ID: 853
The job strain_0_98 was saved and received the ID: 854
The job strain_1_0 was saved and received the ID: 855
The job strain_1_02 was saved and received the ID: 856
The job strain_1_04 was saved and received the ID: 857
The job strain_1_06 was saved and received the ID: 858
The job strain_1_08 was saved and received the ID: 859
The job strain_1_1 was saved and received the ID: 860
job_id: 850 finished
job_id: 851 finished
job_id: 852 finished
job_id: 853 finished
job_id: 854 finished
job_id: 855 finished
job_id: 856 finished
job_id: 857 finished
job_id: 858 finished
job_id: 859 finished
job_id: 860 finished
%% Cell type:code id:straight-honduras tags:
``` python
murn.plot()
```
%%%% Output: display_data
![](data:image/png;base64,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)
%% Cell type:code id:thousand-caribbean tags:
``` python
murn["output/equilibrium_volume"]**(1/3)
```
%%%% Output: execute_result
Output({'error': None, 'residual': 0.162222, 'iterations': 19588})
3.632164318463751
%% Cell type:markdown id:abandoned-reform tags:
%% Cell type:markdown id:impressed-exhaust tags:
### Same for the 1000 structures dataset
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.
### Same cane be done for the 1000 structures dataset
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.
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.
### This can take long and writes 1000 seperate POSCAR files, TAKE CARE
Run these jobs with more cores or more time after the workshop. Also increase the number of iterations for better results
%% Cell type:code id:preceding-march tags:
%% Cell type:code id:pretty-marina tags:
``` python
pr = Project("PotentialDF2")
#pr.remove_jobs()
j = pr.create_job(pr.job_type.Atomicrex, "PotentialDF2")
j.potential = job.potential.copy()
## Use the final parameters as starting values for the new fit
j.potential.copy_final_to_initial_params()
j.input = job.input.copy()
```
%% Cell type:code id:ancient-september tags:
%%%% Output: stream
Are you sure you want to delete all jobs from 'PotentialDF2'? y/(n) y
%% Cell type:code id:proprietary-breed tags:
``` python
df = data_pr.load("df2_1k").to_pandas()
for id, row in df.iterrows():
struct = ase_to_pyiron(row.atoms)
s = j.structures.add_structure(struct, f"id{id}", relative_weight=1)
s.fit_properties.add_FitProperty("atomic-energy", target_value=row.energy/row.number_of_atoms, relative_weight=1)
s.fit_properties.add_FitProperty("atomic-forces", target_value=row.forces, relative_weight=1/row.number_of_atoms)
```
%% Cell type:code id:correct-israeli tags:
%% Cell type:code id:convinced-gambling tags:
``` python
import time
```
%% Cell type:code id:multiple-harbor tags:
%% Cell type:code id:still-counter tags:
``` python
j.input.atom_types.Cu = None
j.input.fit_algorithm = j.factories.algorithms.ar_lbfgs(max_iter=100000)
## if possible increase number of cores
#j.server.cores = 16
t1 = time.time()
j.run()
## Uncomment if you want to run the job
#j.run()
t2 = time.time()
```
%%%% Output: stream
The job PotentialDF2 was saved and received the ID: 623
The job PotentialDF2 was saved and received the ID: 861
%% Cell type:code id:eligible-soviet tags:
%% Cell type:code id:photographic-division tags:
``` python
t2-t1
```
%%%% Output: execute_result
538.0122790336609
850.7202744483948
%% Cell type:code id:subsequent-diploma tags:
%% Cell type:code id:normal-newspaper tags:
``` python
j.output
```