Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
{
"cells": [
{
"cell_type": "markdown",
"id": "f2adfe85-18c8-4430-b370-1b4d6f1b05b3",
"metadata": {},
"source": [
"# Interatomic potentials\n",
"\n",
"In general interatomic potentials can be written as a sum of functional terms depending on the positions of the atoms in a structure. Then the energy $U$ of the system is\n",
"\n",
"$U = \\sum_{i=1}^N U_1(\\vec{r_i}) + \\sum_{i,j=1}^N U_2(\\vec{r_i}, \\vec{r_j}) + \\sum_{i,j,k=1}^N U_2(\\vec{r_i}, \\vec{r_j}, \\vec{r_k}) +...$\n",
"\n",
"The one body term only matter when the system is within an external field. Most classic interatomic potentials don't use 4th and higher order body terms, pair potentials only use the 2 body term. As a general rule potentials that include higher order body terms can be more accurate but are slower.\n",
"\n",
"There are many different forms for $U_i$, which cover a wide range of different run time and accuracy necessities.\n",
"\n",
"Simple pair potentials (f.e Lennard-Jones, Morse, Buckingham) only contain very few parameters, typically less than 10. In many body potentials (f.e. EAM, MEAM, Tersoff) the typical number of parameters is 10-50 and for machine learning potentials the number of parameters can reach several thousands.\n",
"\n",
"# Fitting\n",
"\n",
"In the fit process the parameters of the chosen functions for $U_i$ are optimized. For this purpose an objective or cost function is defined and minimized. In general the objective function is defined as\n",
"\n",
"$\\chi^2 = \\sum_i w_i r_i$\n",
"\n",
"where $w_i$ is a weight and $r_i$ is a residual that describes the difference to target values. This residual can be defined in different ways, so it is not possible to simply compare the residual for different fitting processes or codes. A more in depth explanation and some examples can be found on https://atomicrex.org/overview.html#objective-function.\n",
"\n",
"The minimization can be done with local or global optimization algorithms.\n",
"Generally local optimization algorithms should all be able to find the local minimum coming from some initial parameter set, so the \"best\" local algorithm is the one finding the minimum in the shortest time. Typically used local algorithms are f.e. (L)BFGS or Nelder-Mead.\n",
"Examples for global algorithms are evolutionary algorithms or simulated annealing. For most problems it is impossible to tell a priori which global algorithm will give the best results, so using global algorithms typically involves testing many of them.\n",
"\n",
"# EAM potentials\n",
"\n",
"EAM potentials are pair functionals. \n",
"In a generalised form they are equal to Finnis-Sinclair, effective medium theory or glue potentials. Their total energy can be written as\n",
"\n",
"$E = \\frac{1}{2}\\sum_{ij}V(r_{ij}) + \\sum_i F(\\rho_i)$\n",
"\n",
"with\n",
"\n",
"$\\rho_i = \\sum_j \\rho(r_{ij})$\n",
"\n",
"The original functions for V, $\\rho$ and F were derived from different theories, but they can be chosen freely.\n",
"\n",
"# Fitting code\n",
"\n",
"Fitting is done using the pyiron interface to the atomicrex code https://atomicrex.org. It can be used to fit different types of classic interatomic potentials:\n",
"- pair potentials\n",
"- EAM\n",
"- MEAM\n",
"- Tersoff\n",
"- ABOP\n",
"- ADP (in development)\n",
"\n",
"It allows to fit different properties (energies, forces, lattice constants, elastic properties, etc.) and implements the LBFGS minimizer. Additionally it offers an interface to the nlopt library which implements several global and local optimization algorithms and the ability to apply arbitrary constraints to parameters."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "f4075857-362a-400f-b043-2fc5ad54a407",
"metadata": {},
"outputs": [],
"source": [
"from pyiron import Project, ase_to_pyiron\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import pandas as pd"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "d70cf5ae-8e25-444c-b24f-2d676ef8d2f4",
"metadata": {},
"outputs": [],
"source": [
"pr = Project(\".\")"
]
},
{
"cell_type": "markdown",
"id": "6d83afb5-b99d-4081-9883-ce67c3557f5d",
"source": [
"### Get the training data\n",
"Load a job that contains structures with energies and forces from DFT"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "71b278fc-7886-4863-87f2-50c31e9d3787",
"tc = pr['../../introduction/training/Al_basic_atomicrex']"
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
]
},
{
"cell_type": "markdown",
"id": "57cc1f71-1255-4d15-9c31-e598098fe1e4",
"metadata": {},
"source": [
"### Have a look at the training data"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "25b99c5a-e8a2-42de-aaf0-2794d675385b",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>V</th>\n",
" <th>E</th>\n",
" <th>space_group</th>\n",
" <th>crystal_system</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>16.484415</td>\n",
" <td>225</td>\n",
" <td>cubic</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>14.835973</td>\n",
" <td>225</td>\n",
" <td>cubic</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>15.165661</td>\n",
" <td>225</td>\n",
" <td>cubic</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>15.495350</td>\n",
" <td>225</td>\n",
" <td>cubic</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>15.825038</td>\n",
" <td>225</td>\n",
" <td>cubic</td>\n",
" </tr>\n",
" <tr>\n",
" <th>...</th>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" </tr>\n",
" <tr>\n",
" <th>303</th>\n",
" <td>28.793489</td>\n",
" <td>1</td>\n",
" <td>triclinic</td>\n",
" </tr>\n",
" <tr>\n",
" <th>304</th>\n",
" <td>29.358067</td>\n",
" <td>1</td>\n",
" <td>triclinic</td>\n",
" </tr>\n",
" <tr>\n",
" <th>305</th>\n",
" <td>29.922645</td>\n",
" <td>1</td>\n",
" <td>triclinic</td>\n",
" </tr>\n",
" <tr>\n",
" <th>306</th>\n",
" <td>30.487223</td>\n",
" <td>1</td>\n",
" <td>triclinic</td>\n",
" </tr>\n",
" <tr>\n",
" <th>307</th>\n",
" <td>31.051801</td>\n",
" <td>1</td>\n",
" <td>triclinic</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"<p>308 rows × 4 columns</p>\n",
"</div>"
],
"text/plain": [
" V E space_group crystal_system\n",
"0 16.484415 -3.753374 225 cubic\n",
"1 14.835973 -3.704532 225 cubic\n",
"2 15.165661 -3.723358 225 cubic\n",
"3 15.495350 -3.737149 225 cubic\n",
"4 15.825038 -3.746438 225 cubic\n",
"303 28.793489 -3.016706 1 triclinic\n",
"304 29.358067 -2.979471 1 triclinic\n",
"305 29.922645 -2.942740 1 triclinic\n",
"306 30.487223 -2.906539 1 triclinic\n",
"307 31.051801 -2.870883 1 triclinic\n",
"\n",
"[308 rows x 4 columns]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEMCAYAAADal/HVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA3iElEQVR4nO3deXxU5dXA8d9JCGFfNCCyR18BycYuCIgIBUQUsOCuIK8gLnVrsSBVqVVLhVqxFpWioL60FBUQFQuCQcCNfUcUJAoBIQQDRBLIct4/ZiYmYWYySWYyk8n5fj7zycydZ+49T5Q5uc9z73lEVTHGGGM8iQh2AMYYY0KbJQpjjDFeWaIwxhjjlSUKY4wxXlmiMMYY45UlCmOMMV5VC3YAgRATE6OtW7cOdhjGGFNpbNy48ZiqNnL3XlgmitatW7Nhw4Zgh2GMMZWGiHzv6T0bejLGGOOVJQpjjDFeWaIwxhjjVVjOUbiTk5PDwYMHyc7ODnYoVU6NGjVo3rw5UVFRwQ7FGFMGVSZRHDx4kLp169K6dWtEJNjhVBmqSnp6OgcPHiQ2NjbY4RgTnrYtgJVPwYmDUL859HsCEm/w2+6rTKLIzs62JBEEIsL5559PWlpasEMxJqws3pzKtGV76HLyY6ZWf42anHG8ceIAvP+A47mfkkWVmqOwJBEc9ns3xr8Wb05l0sLtpGZkMaHagl+ShEtOluMMw0+qVKKoLFatWsWQIUPcvjd48GAyMjIqNiBjTEiZtmwPWTl5ADSVY+4bnTjot+NZoqhkli5dSoMGDYIdhjEmiA5lZP3yXGPcN6rf3G/Hs0ThweLNqfSc+gmxEz+k59RPWLw5tdz7fPPNN0lMTCQpKYnbb7+d0aNH88477xS8X6dOnYLnJ0+eZPjw4bRv357x48eTn58POO46P3bsmNv9GWOqhqYNahY8fy73Bk5r9aINomo6JrT9pMpMZpeGa/zPdWqXmpHFpIXbARjWsVmZ9rlz506eeeYZPvvsM2JiYjh+/DiPPPKIx/br1q1j165dtGrVikGDBrFw4UJGjBjhdX/GmKphwsC2Bd9RS/J7QQ78PmoBTSUdCcBVT3ZG4Ubh8T+XrJw8pi3bU+Z9fvLJJ4wYMYKYGMdp4nnnnee1fbdu3bjooouIjIzk5ptvZu3ateXanzEmfAzr2Iw/X59AswY1EWBjvV+xfthqZEoGPLzDr0kC7IzCrcLjf75s94WqnnP1T7Vq1QqGlFSVs2fPFrxXvG3x1+72Z4ypOoZ1bFbmEY7SsjMKNwqP//my3Rf9+vVjwYIFpKenA3D8+HFat27Nxo0bAXjvvffIyckpaL9u3Tr2799Pfn4+//nPf+jVq1eJ+zPGmECwROHGhIFtqRkVWWRbzahIJgxsW+Z9xsXFMXnyZPr06UNSUhKPPPIIY8eO5dNPP6Vbt2589dVX1K5du6B9jx49mDhxIvHx8cTGxjJ8+PAS92eMMYEgqhrsGPyuS5cuWnw9it27d3PppZf6vA/XXY+HMrJo2qAmEwa2rbDTvHBU2t+/MaZiichGVe3i7j2bo/CgIsf/jDEmlAVl6ElERorIThHJFxG3GczZbpCI7BGRvSIysSJjNMYY4xCsOYodwPXAak8NRCQS+AdwNdAeuFlE2ldMeMYYY1yCMvSkqruhxGJx3YC9qvqds+18YCiwK+ABGmMqHZtXDJxQvuqpGXCg0OuDzm1uicg4EdkgIhuspLUxVUvhaqrKL9UU/FF6xwQwUYjIChHZ4eYx1NdduNnm8RItVZ2lql1UtUujRo3KFrQxplIKRDUF84uAJQpV7a+q8W4e7/m4i4NAi0KvmwOH/B9pxUhJSSE+Pj7YYZTLlVdeSfHLjo0JBYGopmB+EcpDT+uBS0QkVkSqAzcBS4IckzEmBAWimoL5RbAujx0uIgeBHsCHIrLMub2piCwFUNVc4H5gGbAbWKCqOyssyG0L4G/xMKWB4+e2BeXeZV5eHmPHjiUuLo4BAwaQlZXFvn37GDRoEJ07d6Z37958/fXX5Obm0rVrV1atWgXApEmTmDx5MgBPPfUUXbt2JT4+nnHjxuG6YXL9+vUkJibSo0cPJkyYUHD2kp2dzZ133klCQgIdO3YkOTkZgLlz53L99dczaNAgLrnkEh599NGCOO+55x66dOlCXFwcTz75ZLn7bYxbfvw3FohqCqYQVQ27R+fOnbW4Xbt2nbPNo63/UX36AtUn6/3yePoCx/Yy2r9/v0ZGRurmzZtVVXXkyJH61ltv6VVXXaXffPONqqp++eWX2rdvX1VV3bFjh7Zr106XL1+uHTp00DNnzqiqanp6esE+b7vtNl2yZImqqsbFxelnn32mqqq///3vNS4uTlVVp0+frqNHj1ZV1d27d2uLFi00KytL58yZo7GxsZqRkaFZWVnasmVL/eGHH4ocIzc3V/v06aNbt25VVdU+ffro+vXry9T/Uv3+TVhbtOmgPvmnJ/TnJ2L8+m9s0aaDevmfV2rr33+gl/95pS7adNCPUYc/YIN6+E61O7PdWfmUY83Zwlxr0JajfG9sbCwdOnQAoHPnzqSkpPD5558zcuTIgjZnzjjWvo2Li+P222/n2muv5YsvvqB6dcfCJMnJyTz33HOcPn2a48ePExcXR+/evTl16hSXX345ALfccgsffPABAGvXruU3v/kNAO3ataNVq1Z88803gKOwYP369QFo374933//PS1atGDBggXMmjWL3NxcDh8+zK5du0hMTCxzv41xcV2d9LH8H7UizhZ9s5z/xqyaQuBYonDH01qz5VyDNjo6uuB5ZGQkR44coUGDBmzZssVt++3bt9OgQQOOHDkCOIaR7r33XjZs2ECLFi2YMmUK2dnZBcNP7nh7r3g8ubm57N+/n+nTp7N+/XoaNmzI6NGjyc7OLmVPjXHPdXVS0+jAr/Ns/CeUJ7ODx9Nas35cgxagXr16xMbG8vbbbwOOL/WtW7cCsHDhQtLT01m9ejUPPPAAGRkZBV/YMTExZGZmFiyj2rBhQ+rWrcuXX34JwPz58wuOccUVVzBv3jwAvvnmG3744QfatvU8bnvy5Elq165N/fr1OXLkCB999JFf+2yqNtdVSBWxzrPxH0sU7vR7wrHmbGF+XoPWZd68ebz22mskJSURFxfHe++9x7Fjx5g4cSKvvfYabdq04f777+fBBx+kQYMGjB07loSEBIYNG0bXrl0L9vPaa68xbtw4evTogaoWDCnde++95OXlkZCQwI033sjcuXOLnEkUl5SURMeOHYmLi2PMmDH07NnT7302VZfrKqSKWOfZ+I+VGfdk2wLHeOmJg46/cvy8Bq2/ZWZmUqdOHQCmTp3K4cOHmTFjRpCj+oWVGTdQdD366yLW8mg1xzrP2bWaUOvq8s0BmvKxMuNlkXhDpfqf9sMPP+TPf/4zubm5tGrVirlz5wY7JGPO4ZpsnrZsD+9n9GJjrV9ZTaZKwBJFmLjxxhu58cYbgx2GMSWyq5MqH5ujMMYY45UlCmOMMV5ZojDGGOOVJQpjjDFeWaKoQqZMmcL06dMBeOKJJ1ixYoXX9kuWLGHq1KkVEZoxJoTZVU9V1FNPPVVim+uuu47rrruuAqIxxoQyO6Pw4MPvPmTAOwNIfCORAe8M4MPvPizX/lJSUmjXrh133XUX8fHx3HrrraxYsYKePXtyySWXsG7dOo4fP86wYcNITEyke/fubNu2DXCcCYwZM4Yrr7ySiy66iBdffLFgv88//zzx8fHEx8fzwgsvFGx/8803SUxMJCkpidtvv/2ceEaPHl1QAqR169Y8+eSTdOrUiYSEBL7++mvAUYr8/vvvB+DIkSMMHz6cpKQkkpKS+Pzzz8v1+zDGVB52RuHGh999yJTPp5Cd56itdPjnw0z5fAoA11x0TZn3u3fvXt5++21mzZpF165d+de//sXatWtZsmQJzz77LC1atKBjx44sXryYTz75hDvuuKOgYODXX39NcnIyp06dom3bttxzzz1s27aNOXPm8NVXX6GqXHbZZfTp04fq1avzzDPP8NlnnxETE8Px48dLjC0mJoZNmzYxc+ZMpk+fzuzZs4u8/8ADD9CnTx8WLVpEXl4emZmZZf49GGMqF0sUbszYNKMgSbhk52UzY9OMciWK2NhYEhISAEcZ8X79+iEiJCQkkJKSwvfff8+7774LwFVXXUV6ejonTpwA4JprriE6Opro6GgaN27MkSNHWLt2LcOHD6d27doAXH/99axZswYRYcSIEcTEOAqvnXfeeSXGdv311wOO8ucLFy485/1PPvmEN998E3BUmnXVkjLGhD8benLjx59/LNV2XxUuxhcREVHwOiIigtzcXLclwUXknM+6SoJ7qtOlqgWfK21srn0bY4yLJQo3mtRuUqrt/lK4JPiqVauIiYmhXr16XtsvXryY06dP8/PPP7No0SJ69+5Nv379WLBgAenp6QA+DT2VpF+/frz88suAY0nXkydPlnufxpjKwRKFGw92epAakTWKbKsRWYMHOz0Y0ONOmTKFDRs2kJiYyMSJE3njjTe8tu/UqROjR4+mW7duXHbZZdx1110FJcInT55Mnz59SEpK4pFHHil3bDNmzCA5OZmEhAQ6d+7Mzp0Vt3y5MSa4rMy4Bx9+9yEzNs3gx59/pEntJjzY6cFyzU9UdVZmvOwWb05l2rI9HMrIommDmlZt1QSElRkvg2suusYSgwm6wus3AKRmZDFp4XYASxamwtjQkzEhzLXGdGFZOXlMW7YnSBGZqsgShTEhzLXGtK/bjQkESxTGhDDXGtO+bjcmEIKSKERkpIjsFJF8EXE7eeJslyIi20Vki4hs8NTOmHA1YWBbakZFFtlWMyqSCQPbBikiUxUFazJ7B3A98KoPbfuq6rEAx2NMSCq8xrRd9WSCJShnFKq6W1VtNs6NZ599tuB5SkoK8fHxQYtl1apVDBkyxO17gwcPJiMjo2IDqqKGdWzGZxOvYv/Ua/hs4lWWJEyFC/U5CgWWi8hGERkX7GACSVXJz88vkijKK5ClOJYuXUqDBg0Ctv9wsnhzKlOefpKDT1xM/pQGnP5LO9i2INhhGeOzgCUKEVkhIjvcPIaWYjc9VbUTcDVwn4hc4eV440Rkg4hsSEtLK3f8J95/n2+v6sfuS9vz7VX9OPH+++XeZ/GS4CkpKVx66aXce++9dOrUif/93/8lKyuLDh06cOuttwKOchljx44lLi6OAQMGkJXluNply5YtdO/encTERIYPH85PP/0EwJVXXsljjz1Gnz59mDFjBitXrqRjx44kJCQwZswYzpw5AzhKiz/22GP06NGDLl26sGnTJgYOHMjFF1/MK6+8UhDzyZMnGT58OO3bt2f8+PHk5+cXfP7YMceIYEklzauyxZtTWbtoJo/mzKR5xDEiUGplHSb3vd9YsjCVh6oG7QGsArr42HYK8Dtf2nbu3FmL27Vr1znbPMlYskR3J3XQXW3bFTx2J3XQjCVLfN5HcRs2bND4+HjNzMzUU6dOafv27XXTpk0qIvrFF18UtKtdu3bB8/3792tkZKRu3rxZVVVHjhypb731lqqqJiQk6KpVq1RV9fHHH9cHH3xQVVX79Omj99xzj6qqZmVlafPmzXXPnj2qqnr77bfr3/72N1VVbdWqlc6cOVNVVR966CFNSEjQkydP6tGjR7VRo0aqqpqcnKzR0dG6b98+zc3N1f79++vbb79d8Pm0tDTdsWOHtmnTRtPS0lRVNT093W3/S/P7DyeX/3mlHnj8ItUn6537eD4u2OEZUwDYoB6+U0N26ElEaotIXddzYACOSfCAO/q3F9DsomXGNTubo397ocz7LFwSvE6dOgUlwVu1akX37t09fi42NpYOHToAjhLgKSkpnDhxgoyMDPr06QPAqFGjWL16dcFnbrzxRgD27NlDbGwsbdq0cdvOtXpdQkICl112GXXr1qVRo0bUqFGjYP6hW7duXHTRRURGRnLzzTezdu3aIvF98sknpS5pXpUcysiiqXi4FuPEwYoNxpgyCtblscNF5CDQA/hQRJY5tzcVkaXOZhcAa0VkK7AO+FBV/1sR8eUePlyq7b5QDzW1XGtJeOKuvHhJXPv0dMzi+y5c8tz12nWc4uXKi7/WMpQ0r0qaNqjJIY1x/2b95hUbjDFlFKyrnhapanNVjVbVC1R1oHP7IVUd7Hz+naomOR9xqvpMRcVX7cILS7XdF55KghcXFRVFTk6O133Vr1+fhg0bsmbNGgDeeuutgrOLwtq1a0dKSgp79+712s6bdevWsX//fvLz8/nPf/5Dr169irwfiJLm4WTCwLa8wE2c1upFtudG1oB+TwQpKmNKJ2SHnoKp8cMPITWKlhmXGjVo/PBDZd6nu5LgDRs2PKfduHHjSExMLJjM9uSNN95gwoQJJCYmsmXLFp544twvnRo1ajBnzhxGjhxJQkICERERjB8/vlRx9+jRg4kTJxIfH09sbCzDhw8v8n4gSpqHk2Edm9Fr+L08F3UvB/NjyEc4XfNCqg39OyTeEOzwjPGJxzLjIlLSyjQCHFbVNn6Pqpz8UWb8xPvvc/RvL5B7+DDVLryQxg8/RP1rr/V3qFWGlRk3JrSVtcz4PlXtWMKON5crshBW/9prLTEYYwzeh55+7cPnfWljjDGmEvOWKB4RkZ7ePqyq3/k5HmOMMSHGW6L4FpjurOD6FxHpUEExGWOMCSEeE4WqzlDVHkAf4DgwR0R2i8gTIhJyE9jGGGMCo8TLY1X1e1X9i3Ni+xZgOLA74JEZY4wJCSUmChGJEpFrRWQe8BHwDTaJXWoZGRnMnDmz3G2CqU6dOsEOwRgTBB4ThYj8SkReBw4C44ClwMWqeqOqLq6g+MJGeRNFXl5eIMIyxpgSebuP4jHgXzgqtla5ugzffPUjX7y3j8zjZ6hzXjQ9hl5Mm8ualHl/EydOZN++fXTo0IFf/epXNG7cmAULFnDmzBmGDx/OH//4x3PaXHPNNfzxj3/kwgsvZMuWLezatYthw4Zx4MABsrOzefDBBxk3zrFMx2uvvcZf/vIXmjZtyiWXXEJ0dDQvvfQS33//PWPGjCEtLY1GjRoxZ84cWrZsyejRo6lXrx4bNmzgxx9/5LnnnmPEiBFkZmYydOhQfvrpJ3Jycnj66acZOrQ0leGNMWHHU1lZLVriuxdwp/N5IyDWl88F61HeMuN7vjysr/wmWV+6e2XB45XfJOueLw/7vI/i9u/fr3FxjrLSy5Yt07Fjx2p+fr7m5eXpNddco59++mmRNqqOMt+1atXS7777rmCbq4z36dOnNS4uTo8dO6apqanaqlUrTU9P17Nnz2qvXr30vvvuU1XVIUOG6Ny5c1VV9bXXXtOhQ4eqquqoUaN0xIgRmpeXpzt37tSLL75YVVVzcnL0xIkTqqqalpamF198sebn56tq0RLopVVVy4wbU1ngpcx4iWtmi8iTQBegLTAHiAL+D/B6j0Vl9sV7+8g9m19kW+7ZfL54b1+5zipcli9fzvLly+nY0XHje2ZmJt9++y0tW7Y8p223bt2IjY0teP3iiy+yaNEiAA4cOMC3337Ljz/+SJ8+fQpKfI8cOZJvvvnG0ZcvvmDhwoUA3H777Tz66KMF+xo2bBgRERG0b9+eI0eOAI4/HB577DFWr15NREQEqampHDlyhCZNyt9vY0zlVGKiwHGVU0dgEzgqvLrWiQhXmcfPlGp7aakqkyZN4u677y6yPSUl5Zy2hcuQr1q1ihUrVvDFF19Qq1YtrrzySrKzs0ssJ15Y4ZLghUuLu/Yxb9480tLS2LhxI1FRUbRu3ZrsYmtzGGOqFl+qx551npYoFCwiFNbqnBddqu2+qFu3LqdOnQJg4MCBvP7662RmZgKQmprK0aNHi7Rx58SJEzRs2JBatWrx9ddf8+WXXwKOs45PP/2Un376idzcXN59992Cz1x++eXMnz8fcCSB4mXC3R2jcePGREVFkZyczPfff1/mPhtjwoMviWKBiLwKNBCRscAK4J+BDSu4egy9mGrVi/5qqlWPoMfQi8u8z/PPP5+ePXsSHx/Pxx9/zC233EKPHj1ISEhgxIgRnDp1qkibCRMmnLOPQYMGkZubS2JiIo8//njBynjNmjXjscce47LLLqN///60b9+e+vXrA46hqjlz5pCYmMhbb73FjBkzvMZ56623smHDBrp06cK8efNo165dmftsjAkPHsuMF2kk8iscS5EKsExVPw50YOXhjzLj/r7qKdAyMzOpU6cOubm5DB8+nDFjxpyzdkQwWZlxY0JbWcuMF3AmhpBODv7W5rImIZ0YipsyZQorVqwgOzubAQMGMGzYsGCHZIwJEx4ThYh8oKpDvH3YlzamYkyfPj3YIRhjwpS3M4peIrLEy/sCtPdzPMYYY0KMt0Thy+24Z/0ViDHGmNDkMVGo6qcVGYgxxpjQ5MvlscYYY6owSxQVpKTqsZdffrnXz48ePZp33nkHgLvuuotdu3Z5bf/KK6/w5ptvlj5QY4wpxpdaT0OApaqaX1Jb45krUdx7771Ftufl5REZGcnnn3/u875mz55dYpvx48eXOsbKZPHmVKYt28OhjCyaNqjJhIFtGdaxWbDDMiYs+XJGcRPwrYg8JyJ+uWNKRKaJyNcisk1EFolIAw/tBonIHhHZKyIT/XFsX+1ek8ys++7krzddy6z77mT3muRy7a9wCfGuXbvSt29fbrnlFhISEoCiiwI999xzJCQkkJSUxMSJ53b7yiuvxHVDYZ06dZg8eTJJSUl07969oLjflClTCi6Z3bt3L/379ycpKYlOnTqxb9++cvUl2BZvTmXSwu2kZmShQGpGFpMWbmfx5tRgh2ZMWPJlKdTbcBQF3Idj3ewvRGRcOQsDfgzEq2oijhXzJhVvICKRwD+Aq3FchnuziFTI5bi71ySzfNZLnDqWBqqcOpbG8lkvlStZTJ06lYsvvpgtW7Ywbdo01q1bxzPPPHPOENJHH33E4sWL+eqrr9i6dWuRaq/u/Pzzz3Tv3p2tW7dyxRVX8M9/nltd5dZbb+W+++5j69atfP7551x44YVl7kcomLZsD1k5RRdyysrJY9qyPUGKyJjw5tMchaqeBN4F5gMX4qgou0lEflOWg6rqclXNdb78Emjuplk3YK+qfqeqZ53HrpAVdNbMf5Pcs0UrxeaePcOa+f4b8y9ePtxlxYoV3HnnndSqVQugoHS4J9WrV2fIEMc9j507dz6nAu2pU6dITU0tKOdRo0aNgn1XVocyskq13RhTPr6smX2tiCwCPsGxFkU3Vb0aSAJ+54cYxuBYi7u4ZsCBQq8POrd5inOciGwQkQ1paWnlCuhU+rFSbS+LwuXDC1PVIqXASxIVFVXQPjIyktzc3CLvl6YEeWXRtEHNUm03xpSPL2cUI4G/qWqiqk5T1aMAqnoax5e8WyKyQkR2uHkMLdRmMpALzHO3CzfbPH7rqeosVe2iql0aNWrkQ7c8q3t+TKm2+7TPEkqIuwwYMIDXX3+d06dPA3D8ePlWoa1Xrx7Nmzdn8eLFAJw5c6Zg35XVhIFtqRkVWWRbzahIJgxsG6SIjAlvJV71pKp3eHlvpZf3+nvbr4iMAoYA/dT9n70HgRaFXjcHDnmP1j9633QHy2e9VGT4qVr1aHrf5PFXUaLCJcRr1qzJBRdc4LbdoEGD2LJlC126dKF69eoMHjyYZ599tszHBXjrrbe4++67eeKJJ4iKiuLtt9/moosuKtc+g8l1dZNd9WRMxSixzLiInOLcv+RPABuA36rqd6U+qMgg4Hmgj6q6HScSkWo4Jrr7AanAeuAWVd1Z0v79UWZ895pk1sx/k1Ppx6h7fgy9b7qDS3v39fnzpigrM25MaCtvmfHncfwl/y8cw0E3AU2APcDrwJVliOklIBr42Dm+/qWqjheRpsBsVR2sqrkicj+wDIgEXvclSfjLpb37WmIwxhh8SxSDVPWyQq9niciXqvqUiDxWloOq6v942H4IGFzo9VJgaVmOYSqpbQtg5VNw4iDUbw79noDEG4IdlTFVmi+T2fkicoOIRDgfhf/Vht8lNSZo1i95layF98OJA4A6fr7/gCN5GGOCxpdEcStwO3AUOOJ8fpuI1ATuD2BspgpZvDmVphufoyZF718hJ8txhmGMCRqvQ0/Ou6PvUdVrPTRZ6/+QTFU0bdke1uDhPpUTBys2GGNMEV7PKFQ1D+hcQbGYKuxQRhaH1MN9KvXd3bhvjKkovgw9bRaRJSJyu4hc73oEPLIwU94y4/60atWqgrIfoaJpg5o8l3sDp7V6ke1ZRDsmtI0xQeNLojgPSAeuAq51PkLrW6YS8JQo8vIcxe1KU2Y8HE0Y2JaPI/swMecuDubHkK9Cqsawo9Of7KonY4LMlzuz76yIQELNz5uPcnJZCnkZZ4hsEE29ga2p3bFxmfdXuMx4VFQUderU4cILL2TLli3s2rWLOnXqkJmZSX5+Pvfffz+ffvopsbGx5OfnM2bMGEaMGMHKlSv53e9+R25uLl27duXll18mOjqa1q1bM2rUKN5//31ycnJ4++23adeuHevWreOhhx4iKyuLmjVrMmfOHNq2Dc0yF7/cbV2d3hm97G5rY0KILwsXtQFeBi5Q1XgRSQSuU9WnAx5dkPy8+SgZC79FcxxrNeVlnCFj4bcAZU4WU6dOZceOHWzZsoVVq1ZxzTXXsGPHjnMqyC5cuJCUlBS2b9/O0aNHufTSSxkzZgzZ2dmMHj2alStX0qZNG+644w5efvllHnroIQBiYmLYtGkTM2fOZPr06cyePZt27dqxevVqqlWrxooVK3jsscd49913y/6LCbBhHZtZYjAmBPky9PRPHOtF5ACo6jYcd2eHrZPLUgqShIvm5HNyWYrfjuGpzPjatWsZOXIkERERNGnShL59HXeH79mzh9jYWNq0aQPAqFGjWL16dcHnrr/eMW1UuNT4iRMnGDlyJPHx8Tz88MPs3FlhN7YbY8KIL4milqquK7Yt123LMJGXcaZU28vCW5nx0mx3iY6OBoqWGn/88cfp27cvO3bs4P333yc7O7scERtjqipfEsUxEbkY513YIjICOBzQqIIsskF0qbb7wtcy47169eLdd98lPz+fI0eOsGrVKgDatWtHSkoKe/fuBRwVYfv06eN1XydOnKBZM8dQzty5c8scuzGmavMlUdwHvAq0E5FU4CHgnkAGFWz1BrZGoor+aiQqgnoDW5d5n4XLjE+YMMFju1//+tc0b96c+Ph47r77bi677DLq169PjRo1mDNnDiNHjiQhIYGIiAjGjx/v9ZiPPvookyZNomfPngVXVxljTGmVWGa8oKFIbSBCVUv+szjI/FFm3N9XPZVGZmYmderUIT09nW7duvHZZ5/RpEmTCjl2oFiZcWNCW7nKjItINPBroDVQzbXspqqGdQGe2h0bV1hiKG7IkCFkZGRw9uxZHn/88UqfJIwxlZsvZcbfw7FQ0UYoXrHNBIJrXsIYY0KBL4miuaoOCngkxhhjQpIvk9mfi0hCwCOpAL7Oxxj/st+7MZWbL4miF7BRRPaIyDYR2S4i2wIdmL/VqFGD9PR0+9KqYKpKeno6NWrUCHYoxpgy8mXo6eqAR1EBmjdvzsGDB0lLSwt2KFVOjRo1aN7cSoUbU1l5TBQicpWqfqKq34tIrKruL/Te9cD3FRKhn0RFRbktmWGMMcY7b0NP0ws9L15J7g8BiMUYY0wI8pYoxMNzd6+NMcaEKW+JQj08d/faGGNMmPI2mX2RiCzBcfbgeo7ztQ32G2NMFeEtUQwt9Hx6sfeKvy4VEZmGY0nVs8A+4E5VzXDTLgU4BeQBuZ7qkBhjjAkcj4lCVT8N4HE/Biapaq6I/AXHwki/99C2r6oeC2AsxhhjvPDlhju/U9Xlqupa/OhLwC6yN8aYEBWURFHMGOAjD+8psFxENorIuAqMyRhjjJMvd2aXiYisANzVx56squ8520zGsazqPA+76amqh0SkMfCxiHytqqvdNXQmknEALVu2LHf8xhhjHHxZj+JjYKRrsllEGgLzVXWgt8+pav8S9jsKGAL0Uw8FmFT1kPPnURFZBHQD3CYKVZ0FzALHwkXejm2MMcZ3vgw9xRS+IklVfwLKtaKPiAzCMXl9naqe9tCmtojUdT0HBgA7ynNcY4wxpedLosgXkYKxHBFpRflvuHsJqItjOGmLiLzi3HdTEVnqbHMBsFZEtgLrgA9V9b/lPK4xxphS8mWOYjKOL2zX5bJX4JwLKCtV/R8P2w8Bg53PvwOSynMcY4wx5VdiolDV/4pIJ6A7jruyH7b7GowxpurwOPQkIu2cPzsBLYFDQCrQ0rnNGGNMFeDtjOIRHENMf3XzngJXBSQiY4wxIcVbCY9xzp99Ky4cY4wxocaX+ygigWuA1oXbq+rzgQvLGGNMqPDlqqf3gWxgO5Af2HCMMcaEGl8SRXNVTQx4JMYYY0KSLzfcfSQiAwIeiTHGmJDkyxnFl8AiEYkAcnDcS6GqWi+gkZkKsXhzKv/+92LaHVxD3bxMouqdx4A77uTS3nYNgzHGwZdE8VegB7DdU/E+Uzkt3pzKP+e+Q+8jyUQ5lwfJPXmcj175O4AlC2MM4NvQ07fADksS4eff/17MlT+uKEgSLpp7ljXz3wxSVMaYUOPLGcVhYJWIfASccW20y2Mrt91rkun0/X+J8FDf8VS6VWkxxjj4kij2Ox/VnQ8TBlbOnUU1L0WA654fU4HRGGNCmS9FAf8I4FwbQlU1M+BRmYDavSaZM5mnPL4v1arT+6Y7KjAiY0woK3GOQkTiRWQzjkWDdjrXr44LfGgmEHavSeajl1/w2ubq8b+xiWxjTAFfJrNnAY+oaitVbQX8FvhnYMMygbB7TTJLX/ormpfnsU2NunUtSRhjivBljqK2qia7XqjqKufSpGFl/ZJXabFpGo01jaPSiAOdJtD1uruDHZbf+HImAXDVqHKtSWWMCUO+JIrvRORx4C3n69twTG6HjfVLXiV+4x+oKWdBoAlp1N/4B9ZD2CSLlXNneT2TAIiuY2cTxphz+TL0NAZoBCx0PmKA0QGMqcK12DTNkSQKqSlnabFpWpAi8q/d/5rKmcyTXttEVKtGv9F2NmGMOZcvZxT9VfWBwhtEZCTwdmBCqniNNc1RmOSc7WFwL8G2Bfz3vdV4/ZtAhEHjH7SzCWOMW76cUUzycVuldVQaedhe+e8lmDl1NvnusqBTRLVqDL7vEUsSxhiPPJ5RiMjVwGCgmYi8WOitekCu+09VTgc6TaC+a47CKUurc6DzBJoEMa7yeuHeseTlReH2dMnJziSMMSXxNvR0CNgAXAdsLLT9FPBQAGOqcF2vu5v14Lzq6RhHJYYDnSv3VU+71ySTm34Y8ZIkwAr/GWNK5m3N7K3AVhH5l6rmuLaLSC/gGeC+CoivwnS97m5wJoYmzkdl9t9XZpSQIiDpV4MrJBZjTOXmSwmPHBHpANwC3IDj0tiFAY7LlMOK2TPJz/U8OqjA+c1b0P+ueysuKGNMpeVxMltE2ojIEyKyG3gJOACIqvZV1b+X56Ai8icR2SYiW0RkuYg09dBukIjsEZG9IjKxPMesSrZ+vNTjewpE1K7PnX99ueICMsZUat6uevoa6Adcq6q9nMnB+x1bvpumqomq2gH4AHiieAMRiQT+AVwNtAduFpH2fjp+2Foxe6bX9wV45PV5FROMMSYseEsUvwZ+BJJF5J8i0g9vl8+UgqoWvvurNritd90N2Kuq36nqWWA+MNQfxw9Xu9ckez2bAMfd18YYUxoeE4WqLlLVG4F2wCrgYeACEXlZRAaU98Ai8oyIHABuxc0ZBdAMx3CXy0HnNk/7GyciG0RkQ1paWnnDq5SW/qPktaTs7mtjTGmVeMOdqv6sqvNUdQjQHNgClDhfICIrRGSHm8dQ534nq2oLYB5wv7tduAvHS5yzVLWLqnZp1Mj9DXThbM5v7wEfVqu1y2GNMaXlSwmPAqp6HHjV+SipbX8fd/sv4EPgyWLbDwItCr1ujuPeDlPM7jXJHD94oMR2dWOqXgI1xpSfLyU8/E5ELin08jocE+fFrQcuEZFYEakO3AQsqYj4Kptlr75YciOwVeuMMWVSqjMKP5oqIm2BfOB7YDyA8zLZ2ao6WFVzReR+YBkQCbyuqjsrIrjFm1OZtmwPhzKyaNqgJhMGtmVYR4/TI0G1YvZM8nJySm6IDTsZY8omKIlCVX/tYfshHPWlXK+XAt4v4/GzxZtTmbRwO1k5jiuBUzOymLRwO0BIJouSrnJysWEnY0xZBWXoKZRNW7anIEm4ZOXkMW3ZniBF5Nmc397jU7tq1aNt2MkYU2bBGnoKWYcyskq1PVhWzJ7p8wR275vusGEnY0yZWaIopmmDmqS6SQpNG9QMQjSe+TLkNPj+31qCMMaUmw09FTNhYFtqRkUW2VYzKpIJA9sGKaJzlVSmAxyVYS1JGGP8wc4oinFNWIfyVU8lnU1IZKRVhjXG+I0lCjeGdWwWUomhsAV/mlxim6vveSjwgRhjqgwbeqpEVsyeyYEdW722sSEnY4y/WaKoREoacoqMirIhJ2OM31miqCR8mcAeePcDFRCJMaaqsURRCfiyzoQNORljAsUSRSXw8ex/eH3frnIyxgSSJYoQt2L2THKys722saucjDGBZIkihNmQkzEmFFiiCGFr5r/p9f3I6GgbcjLGBJwlihB2Kv2Y1/cHjnW3gqwxxviX3ZkdQnavSWbN/Dc5lX6MuufHUKNOHbJPnXLb1oacjDEVxRJFEOxek8zKubM4k+lIAjXq1qVt997s/HQluWfPAHDqWJrHzyf9arANORljKowligq2YvbMcyaos0+d8qlsuK0tYYwJBksUFciXq5g8qRvTiHH/mOPniIwxpmQ2mV1Bdq9J5qOZfyvz50ua2DbGmECxRFEBdq9JZvmsl9D8/DLvo+75MX6MyBhjfGeJogKsmf9mwSR1WfW+6Q4/RWOMMaVjiaIClHvYSMQmsI0xQWOJogJERUeX6/NJ/a/2UyTGGFN6ligCbPea5BKL+nlzXvMWds+EMSaognJ5rIj8CRgK5ANHgdGqeshNuxTgFJAH5Kpql4qM06X4HdOluZehpHpN3lSLjubOv75c5s8bY4w/BOuMYpqqJqpqB+AD4AkvbfuqaodgJonls15y3CmtyqljaSyf9RK71yT79PnyzE8MsFpOxpgQEJREoaonC72sDWgw4vCFuyuWcs+e8flMoayXtUZGR9sEtjEmJATtzmwReQa4AzgBePpGVGC5iCjwqqrO8rK/ccA4gJYtW/otTk9nBN5qMbnsXpPM2TLOT1hlWGNMqAjYGYWIrBCRHW4eQwFUdbKqtgDmAZ6+FXuqaifgauA+EbnC0/FUdZaqdlHVLo0aNfJbP7ydEXgbfloxeyZLX/prQeG/0rDKsMaYUBKwRKGq/VU13s3jvWJN/wX82sM+Djl/HgUWAd0CFa8n3m508zT8VJ6aTlYZ1hgTaoIyRyEilxR6eR3wtZs2tUWkrus5MADYUTER/sLbX/aehp8+nv2PMh1LIiIsSRhjQk6wrnqa6hyG2oYjATwIICJNRcT1p/gFwFoR2QqsAz5U1f8GI9i6MZ6HstwNP5X1vgnVkJ3TN8ZUYUGZzFZVb0NNg53PvwOSKjIuT3rfdAdLX/qr2/fmv/wqKz5TJgxsy7COzcp1HCv8Z4wJRXZntg+8DT/VzcskNSOLh/+zhT8s3s6K2TPLfBwr/GeMCUW2cJGPJCLCY5nwUT+8xecNL+Orj7+hwbGVSBn2b1c6GWNClSUKH3lKEgLUy8vk2cjzaFIzFup0Lnjv25Ob2Hx8RYn7jq5T1yaxjTEhy4aefORtQvuKC0bSpGYsIlLkcUm9TnQ8r3+J++43epw/QzXGGL+yROEjb/MHriRRnIjwP/U6et1vi/gkG3IyxoQ0SxQ+urR3X2rUrVvqz4mXGQuJjOSGx58pT1jGGBNwlihK4apR4yht/UL10F6Bq+95qNwxGWNMoFmiKIVLe/flfxr8RPFk8WPWfrc3y6kqe09uPnc7sO/8RBtyMsZUCpYoSmnl+X1IrH+IwsmiRmRtwJEYCj/cXfWkwK568Qwce18FRm2MMWVniaKUluT34lDjiwqSxYCmo2lQvfE5VzwpSvqZ1CKfdZ1JXHfvb8p9F7cxxlQUSxSl1LBWFE/mjilIFq4kUVyERJDYsI/zlZIPrG91NS/PfNaShDGmUrFEUUpPXhvn+Jk7hvdjBnltW6taPUBJi2zIzNh7uPnmYYEP0Bhj/MwSRSkN69iM27o7VtBbkt/La9vTuSf5b0x/5re8mVu7t7QzCWNMpWQlPMrg6WEJAPzflz9wPE85L5Jzhp/yVZkWFcWh2pfywvAESxLGmErLEkUZPT0sgS6tzuOmhduYn1WL8yJ/eS8f+JNkMeTGeGZbgjDGVHKWKMphWMdmHs8UZldwLMYYEyg2R2GMMcYrSxTGGGO8skRhjDHGK0sUxhhjvLJEYYwxxitxV/W0shORNOD7YMfhRQxwLNhBBFi49zHc+wfh38dw7x+Uro+tVNXtUp5hmShCnYhsUNUuwY4jkMK9j+HePwj/PoZ7/8B/fbShJ2OMMV5ZojDGGOOVJYrgmBXsACpAuPcx3PsH4d/HcO8f+KmPNkdhjDHGKzujMMYY45UlCmOMMV5ZogggEWkhIskisltEdorIg87t54nIxyLyrfNnw2DHWl4iEikim0XkA+frsOqjiDQQkXdE5Gvnf88e4dRHEXnY+f/oDhH5t4jUqOz9E5HXReSoiOwotM1jn0RkkojsFZE9IjIwOFGXjoc+TnP+f7pNRBaJSINC75Wpj5YoAisX+K2qXgp0B+4TkfbARGClql4CrHS+ruweBHYXeh1ufZwB/FdV2wFJOPoaFn0UkWbAA0AXVY0HIoGbqPz9mwsUX6/YbZ+c/y5vAuKcn5kpIpGEvrmc28ePgXhVTQS+ASZB+fpoiSKAVPWwqm5yPj+F48ulGTAUeMPZ7A1gWFAC9BMRaQ5cQ9FlOMKmjyJSD7gCeA1AVc+qagZh1Ecca9PUFJFqQC3gEJW8f6q6GjhebLOnPg0F5qvqGVXdD+wFulVEnOXhro+qulxVc50vvwSaO5+XuY+WKCqIiLQGOgJfAReo6mFwJBOgcRBD84cXgEdxLO7nEk59vAhIA+Y4h9dmi0htwqSPqpoKTAd+AA4DJ1R1OWHSv2I89akZcKBQu4PObZXdGOAj5/My99ESRQUQkTrAu8BDqnoy2PH4k4gMAY6q6sZgxxJA1YBOwMuq2hH4mco3DOORc5x+KBALNAVqi8htwY2qwombbZX63gERmYxj+Huea5ObZj710RJFgIlIFI4kMU9VFzo3HxGRC53vXwgcDVZ8ftATuE5EUoD5wFUi8n+EVx8PAgdV9Svn63dwJI5w6WN/YL+qpqlqDrAQuJzw6V9hnvp0EGhRqF1zHMNvlZKIjAKGALfqLzfLlbmPligCSEQEx7j2blV9vtBbS4BRzuejgPcqOjZ/UdVJqtpcVVvjmCj7RFVvI7z6+CNwQETaOjf1A3YRPn38AeguIrWc/8/2wzGfFi79K8xTn5YAN4lItIjEApcA64IQX7mJyCDg98B1qnq60Ftl76Oq2iNAD6AXjlO7bcAW52MwcD6OKy6+df48L9ix+qm/VwIfOJ+HVR+BDsAG53/LxUDDcOoj8Efga2AH8BYQXdn7B/wbx5xLDo6/pv/XW5+AycA+YA9wdbDjL0cf9+KYi3B957xS3j5aCQ9jjDFe2dCTMcYYryxRGGOM8coShTHGGK8sURhjjPHKEoUxxhivLFEYY4zxyhKFMVWciPxdRDaJSNdC2y4VkVecpdXvCWZ8JvgsURhThTmLGzYG7sZR8gEAVd2tquOBG4AuQQrPhAhLFKbSEZHhIqIi0s75uoGI3OvnY3xeirarii8CIyIPicjMEj6XWdb4SktEWotIlohsKbxdVX8GLgRWAS8W+8x1wFpgpYjUFJEtInJWRGIqKGwTIixRmMroZhxfYDc5XzcA/JooVPXyUjT/d6FYXG5ybg8l+1S1Q+ENInI+jvUnTgF5hd9T1SXO38Otqprl/GylLZRnys4ShalUnCXbe+KoaeP6cp4KXOz8i3eas90jzmU9d4jIQ85trZ1LRM52bp8nIv1F5DPn0pjdCh0ns9DzO5zLSm4VkbfchPUOMEREol3HwVGue627OIr1p3WxZSx/JyJTShOviNwmIuuc/X+1lCuz/QHHWhQ7gfaF4rhSRF4UkVeBpaXYnwlHwS5qZQ97lOYB3Aa85nz+OY5y362BHYXadAa2A7WBOji+BDs62+UCCTj+SNoIvI6jTv9QYHGhfWQ6f8bhKKAW43zttjAe8CEw1Pl8IjDNUxyFj+Em9t8BU5zPS4wXuBR4H4hyfmYmcIeb+Iocp9C2r5z7ewkY68PvP8X1u7BH1XnYGYWpbG7Gse4Fzp83u2nTC1ikqj+raiaO9RV6O9/br6rbVTUfxxf3SlVVHF/ord3s6yrgHVU9BqCqxZfWdCk8/OQadvIWh69KircfjoS03jn/0A/Hiny+eBp4yrm/3TiSojHnqBbsAIzxlXM8/SogXkQUiMRRxr34pLG7lbxczhR6nl/odT7u/z0Ivq0Cthh4XkQ6ATVVdZOIXOHD53IpOgRco5TxCvCGqk7y4VgFRKQDcD3QS0T+4Tzu9tLsw1QddkZhKpMRwJuq2kpVW6tqC2A/0BKoW6jdamCYcyGe2sBwYE0Zj7kSuMGZpBCR89w1cp4xrMIxNOSaxPYljiNAYxE53znHMYTSWQmMEJHGrvhEpJUPn/sLcK3z99gaSMLOKIwHlihMZXIzsKjYtndxDPV85pzwnaaqm4C5OFbv+gqYraqby3JAVd0JPAN8KiJbgee9NP83ji/c+c7PlhiHOpYefcr5/gc4Fg8qTXy7cExILxeRbcDHOC539UhErgJqq+rKQvs5gmOtbLeJ0FRttnCRMVWA80qsD1Q1vpz7SQG6uOZsTNVgZxTGVA15QP3iN9z5ynXDHRCFY37EVCF2RmGMMcYrO6MwxhjjlSUKY4wxXlmiMMYY45UlCmOMMV5ZojDGGOOVJQpjjDFeWaIwxhjjlSUKY4wxXv0/xEaCskssIIAAAAAASUVORK5CYII=\n",
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"tc.plot.energy_volume(crystal_systems=True)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "f039bd7e-e62e-43d6-9c22-6a4bb1c68c73",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>a</th>\n",
" <th>b</th>\n",
" <th>c</th>\n",
" <th>alpha</th>\n",
" <th>beta</th>\n",
" <th>gamma</th>\n",
" <th>V</th>\n",
" <th>N</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>4.039967</td>\n",
" <td>4.039967</td>\n",
" <td>4.039967</td>\n",
" <td>90.000000</td>\n",
" <td>90.000000</td>\n",
" <td>90.000000</td>\n",
" <td>65.937658</td>\n",
" <td>4</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>3.900545</td>\n",
" <td>3.900545</td>\n",
" <td>3.900545</td>\n",
" <td>90.000000</td>\n",
" <td>90.000000</td>\n",
" <td>90.000000</td>\n",
" <td>59.343892</td>\n",
" <td>4</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>3.929227</td>\n",
" <td>3.929227</td>\n",
" <td>3.929227</td>\n",
" <td>90.000000</td>\n",
" <td>90.000000</td>\n",
" <td>90.000000</td>\n",
" <td>60.662645</td>\n",
" <td>4</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>3.957496</td>\n",
" <td>3.957496</td>\n",
" <td>3.957496</td>\n",
" <td>90.000000</td>\n",
" <td>90.000000</td>\n",
" <td>90.000000</td>\n",
" <td>61.981399</td>\n",
" <td>4</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>3.985366</td>\n",
" <td>3.985366</td>\n",
" <td>3.985366</td>\n",
" <td>90.000000</td>\n",
" <td>90.000000</td>\n",
" <td>90.000000</td>\n",
" <td>63.300152</td>\n",
" <td>4</td>\n",
" </tr>\n",
" <tr>\n",
" <th>...</th>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" </tr>\n",
" <tr>\n",
" <th>303</th>\n",
" <td>7.008729</td>\n",
" <td>7.232998</td>\n",
" <td>10.690023</td>\n",
" <td>121.764941</td>\n",
" <td>89.204727</td>\n",
" <td>90.106498</td>\n",
" <td>460.695818</td>\n",
" <td>16</td>\n",
" </tr>\n",
" <tr>\n",
" <th>304</th>\n",
" <td>7.054241</td>\n",
" <td>7.279967</td>\n",
" <td>10.759441</td>\n",
" <td>121.764941</td>\n",
" <td>89.204727</td>\n",
" <td>90.106498</td>\n",
" <td>469.729069</td>\n",
" <td>16</td>\n",
" </tr>\n",
" <tr>\n",
" <th>305</th>\n",
" <td>7.099174</td>\n",
" <td>7.326338</td>\n",
" <td>10.827974</td>\n",
" <td>121.764941</td>\n",
" <td>89.204727</td>\n",
" <td>90.106498</td>\n",
" <td>478.762321</td>\n",
" <td>16</td>\n",
" </tr>\n",
" <tr>\n",
" <th>306</th>\n",
" <td>7.143545</td>\n",
" <td>7.372128</td>\n",
" <td>10.895651</td>\n",
" <td>121.764941</td>\n",
" <td>89.204727</td>\n",
" <td>90.106498</td>\n",
" <td>487.795572</td>\n",
" <td>16</td>\n",
" </tr>\n",
" <tr>\n",
" <th>307</th>\n",
" <td>7.187372</td>\n",
" <td>7.417357</td>\n",
" <td>10.962497</td>\n",
" <td>121.764941</td>\n",
" <td>89.204727</td>\n",
" <td>90.106498</td>\n",
" <td>496.828823</td>\n",
" <td>16</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"<p>308 rows × 8 columns</p>\n",
"</div>"
],
"text/plain": [
" a b c alpha beta gamma \\\n",
"0 4.039967 4.039967 4.039967 90.000000 90.000000 90.000000 \n",
"1 3.900545 3.900545 3.900545 90.000000 90.000000 90.000000 \n",
"2 3.929227 3.929227 3.929227 90.000000 90.000000 90.000000 \n",
"3 3.957496 3.957496 3.957496 90.000000 90.000000 90.000000 \n",
"4 3.985366 3.985366 3.985366 90.000000 90.000000 90.000000 \n",
".. ... ... ... ... ... ... \n",
"303 7.008729 7.232998 10.690023 121.764941 89.204727 90.106498 \n",
"304 7.054241 7.279967 10.759441 121.764941 89.204727 90.106498 \n",
"305 7.099174 7.326338 10.827974 121.764941 89.204727 90.106498 \n",
"306 7.143545 7.372128 10.895651 121.764941 89.204727 90.106498 \n",
"307 7.187372 7.417357 10.962497 121.764941 89.204727 90.106498 \n",
"\n",
" V N \n",
"0 65.937658 4 \n",
"1 59.343892 4 \n",
"2 60.662645 4 \n",
"3 61.981399 4 \n",
"4 63.300152 4 \n",
".. ... .. \n",
"303 460.695818 16 \n",
"304 469.729069 16 \n",
"305 478.762321 16 \n",
"306 487.795572 16 \n",
"307 496.828823 16 \n",
"\n",
"[308 rows x 8 columns]"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABIQAAAG/CAYAAAApVwz/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA3E0lEQVR4nO3dfbxlZ1kf/N9lQgigHgUC1bw4oROjkdSAY0RrLVXUYBxCRSVBIVqaMW1B7aPV8GgfU6s1PI9WG8XiWCFRMUjRQoaE4ksLSIWa8OYkxEDAgQxBEkCOCFSIXs8few0cTs6ZnJnzsvY56/v9fPZnzr73Wmtfa79cs9a17/te1d0BAAAAYDo+a+wAAAAAANhaCkIAAAAAE6MgBAAAADAxCkIAAAAAE6MgBAAAADAxCkIAAAAAE6MgxDGrqjOq6q+r6oQtfM5XV9U/36rnA6apqv5RVd0+dhzAsfP9HV9VPb6qDo8dB8wTuWnzVNWVVfWbY8exnSkIbXNDoeQvq+qBy9oPVdUTNuM5u/s93f3Z3f23a12nqr66qj5aVZ+zwmNvrqpnbWyUwHY35LGPV9VHqurDVfXHVXV5VW3a/13d/UfdffayGDYll8IUHe93aqUfhqqqq2r3kfvLv7/rVVUnD7nn61d47Oer6qXr2PaG55blr8dWGOM5YTPITZ9af1OOe2rmXVX1to3eNuujILSNVdWuJP8oSSd50rjRHF13vz7J4SRPWdpeVY9Ock6S68aIC5h7e7v7c5J8UZKrkvxokl8bNyRgCrr7/yT57STPWNo+9JC+JMm1Y8RVVSeO8bzAfNimuenrkjwiyaOq6iu3KCTWQEFoe3tGkjckuSbJpUcaq+o3kpyR5MAwtOtHhvYnVdWtQ0X51VX1pUvWOVRV/6aq/nToyfNrVfXIqnrl8Ov8H1TV5w/L7hoq3ycO9x9aVS+sqruG3kovWyXea7MscQ33b+juD1bV11TVTVW1OPz7NSttZHnXwBXieXVV/dTQm+Cvq+pAVT2sql5UVX81bHvXkvW/pKp+v6o+VFW3V9V3ruXFB7ZOdy929/VJnprk0qp6dFU9sKp+tqreU1Xvr6rnV9WDkk8PW6iqH6qqu6vqfVX1vUe2V1XfUlVvG/Lbe6vqh5euN/x9n1xaVTdU1bOXxjbkzSdv0UsBO1JVfX5VvaKq7hmOJV5RVacNj/10Zj+A/dLwXfylqnrtsOpbh7an1rLhSlV1elX97rDND1bVLy157J9V1W3Dc72qqr5oldCuTfKUqnrwkrZvzuwY+pVVtTAcM71vyCU/VUuG1FfVZcPzfGTIOY9dx3Haj1bVnyb5aB1DUWidufJhw3HUkeOnn6qq1w2P3ec9WLLeMeVemFdy04blpkuTvDzJjVly3jps49VV9e+r6n8N8fxeVT18yePPqKp3D6/Vv62j9GKqqsfV7Bzww1X11qp6/JLHvqdmvZQ+UlV/XlXftUqs09Ldbtv0luSOJP8yyVck+WSSRy557FCSJyy5/8VJPprkG5M8IMmPDOuftGT5NyR5ZJJTk9yd5E1JHpPkgUn+R5KfGJbdlVmvpBOH+zdkVqX+/GHb/3iVeE8f4jxjuP9ZmfUaenKShyb5yyRPT3JiZtXtv0zysGHZVyf558PfVyb5zSXbXR7Pq4d9+/tJFpK8Lcnbkzxh2PavJ3nhsOxDktyZ5HuHxx6b5ANJvmzs99fNbeq35XlsSft7kvyLJL+Q5Pohf3xOkgNJfmZY5vFJ7k3yk0Ne+pYkH0vy+cPj70vyj4a/Pz/JY5esd3i1GJJ8Z5L/veT+lyf54JFc6ubmdvTbUb7XD8usF/GDh+/zf03ysiWPf+o4YElbJ9m95P6nvr9JTkjy1iQ/P/xff3KSrx0ee/JwnPClw//9P57kj48S89uTfPeS+9cl+YXh75cl+ZXhOR6R5E+SfN/w2HckeW+Sr0xSSXYn+aKVXoes7TjtLZkdSz1olTg/4/VY0r6eXPni4fbgzHp035nkdffzHhxz7nVzG/smN21qbnpwkr8a8sFTMjvXOmnZa/jO4bkeNNy/anjsnCR/neRrk5yU5GczO598wvD4lRnOCzM7h/3g8DyfNcT8wSSnDK/DXyU5e1j2C+J8L92th9B2VVVfm9kQipd09xsz+xI97SirPDWznji/392fzOzL9KAkS3vh/GJ3v7+735vkjzI76Xlzd/9Nkv+WWXFoeRxfkOSJSS7v7r/s7k9292tWCqC770zymiTfPTR9Q2ZJ8IYkFyZ5R3f/Rnff293XJfmzJHvX9ILc1wu7+53dvZjklUne2d1/0N33ZpbIj+zLtyY51N0vHJ73TUl+J8m3H+fzApvvrsxObC5L8q+7+0Pd/ZEk/yHJxUuW+2SSnxzy0o2ZHVCcveSxc6rqc4fc9aY1PvfLk5xVVWcN95+e5Le7+xPr3CeYtO7+YHf/Tnd/bPg+/3SSf7yOTZ6f5AuT/Jvu/mh3/5/uft3w2PdlVhC5bTgu+A9JzjvKL/G/nqGHc1V9bpKLklxbVY/M7BjoB4fnuDuzk7wjeeifJ/l/u/umnrmju9+9ynOs5Tjt6u6+s7s/vtYXoaoqx5krh94ET8nsB8GPdffbsrahKJuRe2EUctOG5KZvS/I3SX4vySsyK3ZduGyZF3b324dtvCTJeUP7tyc50N2vG461/p/MCm4r+e4kN3b3jd39d939+0luzqxAlCR/l+TRVfWg7n5fd9+6ynYmRUFo+7o0ye919weG+7+VZd3vlvnCJJ/6onf332X2K8+pS5Z5/5K/P77C/c9eYbunJ/lQd//lGuNeOmzs6Ul+a0gunxHf4N3L4jsWa92XL0ryVUO3wg9X1YeTfFeSv3eczwtsvlMzO5h4cJI3Lvnu/vfMfgU64oPDAdURH8unv/tPyewA4d1V9Zqq+uq1PPFQIH9Jku+u2eTWlyT5jfXsDJBU1YOr6leGYQF/leS1ST6vjv+KpqcnefeyHHDEFyX5T0tyx4cy+5V8tWOOX0/yT6rq1MxOTu7o7jcP23lAkvct2davZPZr/JEY3rnGeNdynHbnGre11Ck5/lx5Sma5dunzriWGDc+9MBa5aUNy06WZdWK4dziO+t3c97z1L5b8vTRnfOHS7Xf3xzLr9bOSL0ryHcvO6742yRd090czK25dntnrckNVfcn9xD0JJqXbhmo27vs7k5xQVUe+PA/MLDl9eXe/NfetnN6V5Nwl26jMksF71xnOnUkeWlWf190fXsPyv5vkl6vqn2RWLX78kviWV7/PyOygZbmPZnZwc8R6ijd3JnlNd3/jOrYBbJGaTUR4amZdoX80s+6+x5zHuvumJBdV1QOSPCuzIs/pKy26Qtu1mRWBXpfkYz2bNB9Ynx/KrBfJV3X3X1TVeUnenNnJULL6L8KruTPJGVV14gonXncm+enuftFaNtTd76mqP8rsB6MnZnYSdmQ7f5Pk4auc3N2Z2fD1FTe77P5ajtOO9TVIZkMzPp7jy5X3ZDb867TMhqYkK+fJNTuG3AvzQm5aR26q2XxLX5/k/Ko6cnGhByc5uaoevqRzw2rel0/3MDxyHvywVZa9M8lvdPdlKz3Y3a9K8qphGz+V5FczmwNq0vQQ2p6enORvMxtTed5w+9LMhnkd6X3z/iSPWrLOS5JcWFXfMPwn/EOZJYo/Xk8g3f2+zIZk/XLNJl17QFV93VGW/2iSlyZ5YWbV8ZuHh25M8sVV9bSqOrFmExOek1m3wuXekuTrquqMqlpI8px17MIrhud9+hD7A6rqK2vJZGnA+Krqc6vqWzOby+I3h8L3ryb5+ap6xLDMqVX1zWvY1klV9V1VtTD0UPyrzHLqSpbn0gwFoL9L8nPROwiOxwNqdtnkI7cTM5ub4+NJPlxVD03yE8vWuc93cZW2I/4ksxOJq6rqIcPz/MPhsecneU5VfVmS1Gzy1e+4n5ivzayA8Q+TvCj51DHQ7yX5uSFHfVZV/f2qOjKc5L8k+eGq+oqa2b1k6MdmHaedtPS1zeyk9bhyZXf/bWY/5F059JL4ktz34iBHew8+wzHmXhiD3LTxuenpmRWUz86nz1u/OLN5ZC9Zw/ovTbK3ZhcfOinJv8uni3HL/eaw7DdX1QnDa/v4qjqtZhdLelJVPWSI/68j/yRRENquLs1snOV7uvsvjtyS/FKS7xqS188k+fGhu9wPd/ftmY2r/MXMfi3am9nlnDdi3ounZzYm/M8ym4z6B+9n+Wsz6w10pIqd7v5gZvP5/FBm3QB/JMm3rlQ1HsaD/naSP03yxqxcNFqTYSzwN2U2pvauzLorPjezHlfA+A5U1Ucy+9Xnx5L8x8wmgU9mPYTuSPKGmnXj/oMs+RXpfjw9yaFhvcvz6bnNlvuMXLqk/dcz+8XsN1deDTiKGzM7wTpyuzKziY8flNkxyhty3x7C/ynJt9fsyjtXD21XZjZfxodr2RVCh2LG3swmS31PZicfTx0e+2+Z/V//4iEH3JLZr+tH89LMJkH+w+Fk64hnZDbR6dsyuxjGSzObrDTd/V8zm2/kt5J8JLOejQ8d1tus47Rb85mv7fdmfbnyWZldoOMvMiuAX5fZydQRV2aV92AVa829MAa5aeNz06VJfnnpOetw3vr8HH26kwyx3prk2Zn9IPi+Id6785l56Miyd2Y2j9L/nVkPxzuT/JvMah6fldl55l2ZDcX7x5ldnGnyqvt4ep8CwHRV1TOS7Ovurx07FoCtUlXPTfL3uvt+T+QANlpVfXaSDyc5q7v/fORwdgQ9hADgGFTVgzP7VWn/2LEAbKaq+pKq+gfDsJLzkzwzsyvPAmyJqto7DFt9SGZXOTuY2eXu2QAKQgCwRsO8G/dkNsb+t0YOB2CzfU5m8wh9NLO5RH4uyctHjQiYmosyG+p1V5KzklzchjltGEPGAAAAACZGDyEAAACAiTlx7ACS5OEPf3jv2rVr7DCAdXrjG9/4ge4+Zew4jpdcBDuDXATMA7kImAdHy0VzURDatWtXbr755rHDANapqt49dgzrIRfBziAXAfNALgLmwdFykSFjAAAAABOjIAQAAAAwMQpCAAAAABOjIAQAAAAwMQpCAAAAABOjIATsWFX15Kr61ap6eVV909jxAAAAzAsFIWBbqaoXVNXdVXXLsvYLqur2qrqjqq5Iku5+WXdfluR7kjx1hHABAADmkoIQsN1ck+SCpQ1VdUKS5yV5YpJzklxSVecsWeTHh8cBAACIghCwzXT3a5N8aFnz+Unu6O53dfcnkrw4yUU189wkr+zuN620varaV1U3V9XN99xzz+YGDwAAMCcUhICd4NQkdy65f3hoe3aSJyT59qq6fKUVu3t/d+/p7j2nnHLK5kcKAAAwB04cOwCADVArtHV3X53k6q0OBgAAYN6N2kOoqvZW1f7FxcUxwwC2v8NJTl9y/7Qkd40UCwAAwNwbtSDU3Qe6e9/CwsKYYQDb301JzqqqM6vqpCQXJ7l+5JgAAADmljmEgG2lqq5L8vokZ1fV4ap6Znffm+RZSV6V5LYkL+nuW8eMEwAAYJ6ZQwjYVrr7klXab0xy4xaHAwAAsC3pIQQAAAAwMTu6h9CuK25Y87KHrrpwEyMB5llV7U2yd/fu3WOHwjbm/xyYY1euMl/llS5sAtuW7zWs244uCAGsRXcfSHJgz549l40dC7DzVNWTk1yY5BFJntfdvzduRDvXaoXZQydvcSAAsA0YMgYAcIyq6gVVdXdV3bKs/YKqur2q7qiqK5Kku1/W3Zcl+Z4kTx0hXACA+1AQAgA4dtckuWBpQ1WdkOR5SZ6Y5Jwkl1TVOUsW+fHhcQCA0SkIAQAco+5+bZIPLWs+P8kd3f2u7v5Ekhcnuahmnpvkld39ppW2V1X7qurmqrr5nnvu2dzgAQCiIAQAsFFOTXLnkvuHh7ZnJ3lCkm+vqstXWrG793f3nu7ec8opp2x+pADA5JlUGgBgY9QKbd3dVye5equDAQA4Gj2EAAA2xuEkpy+5f1qSu0aKBQDgqBSEAAA2xk1JzqqqM6vqpCQXJ7l+5JgAAFakIARMXlXtrar9i4uLY4cCbBNVdV2S1yc5u6oOV9Uzu/veJM9K8qoktyV5SXffOmacAACrMYcQMHndfSDJgT179lw2dizA9tDdl6zSfmOSG7c4HACAY6aHEAAAAMDEKAgBAMwBw1cBgK2kIAQAMAe6+0B371tYWBg7FABgAhSEAAAAACZGQQgAAIAkhq/ClCgIAQAAkMTwVZgSl50HAABgRzj32nNXbD946cEtjgTmnx5CAAAAABOjIAQAAAAwMQpCwOSZPBGYB3IRALCVFISAyTN5IjAP5CIAYCspCAEAAABMjIIQAAAAwMQoCAEAAABMjIIQAAAAwMQoCAEAAABMjIIQAAAAwMQoCAEAAABMjIIQAAAAwMQoCAEAzIGq2ltV+xcXF8cOBQCYAAUhAIA50N0HunvfwsLC2KEAABOgIARMnl/lAQCAqVEQAibPr/IAAMDUKAgBAAAATMyJYwcAAMDx2XXFDSu2Hzr5aSuvcKWhsQDAjB5CAAAAABOjhxAAAOt27rXnrth+8NKDWxwJALAWm9JDqKqeXFW/WlUvr6pv2oznAAAAAOD4rLkgVFUvqKq7q+qWZe0XVNXtVXVHVV2RJN39su6+LMn3JHnqhkYMAAAAwLocSw+ha5JcsLShqk5I8rwkT0xyTpJLquqcJYv8+PA4AAAAAHNizQWh7n5tkg8taz4/yR3d/a7u/kSSFye5qGaem+SV3f2mlbZXVfuq6uaquvmee+453vgBAAAAOEbrnVT61CR3Lrl/OMlXJXl2kickWaiq3d39/OUrdvf+JPuTZM+ePb3OOAAAtrWq2ptk7+7du8cOhWV2XXHDiu2HrrpwxXYTbAOwHax3Uulaoa27++ru/oruvnylYhAAAJ+puw90976FhYWxQwEAJmC9BaHDSU5fcv+0JHetc5sAAAAAbKL1FoRuSnJWVZ1ZVScluTjJ9esPC2DrVNXeqtq/uLg4digAAABb4lguO39dktcnObuqDlfVM7v73iTPSvKqJLcleUl337o5oQJsDsM0AACAqVnzpNLdfckq7TcmuXHDIgIAgO3gylV+SDjzjK2NA3awVSd1P3mLA4EdaL1DxtbFMA0AAACArTdqQcgwDQAAAICtN2pBCAAAAICtpyAEAAAAMDEKQgAAAAAToyAEAAAAMDEKQgAAAAATc+LYAQAAkFTV3iR7d+/ePXYok3Huteeu2H7w0oNbHAkAbL1RewhV1d6q2r+4uDhmGAAAo+vuA929b2FhYexQAIAJGLUg5MAHAAAAYOuZQwgAAABgYhSEAAAAACZGQQgAAABgYlxlDIDJ2HXFDWte9tBVF25iJAAAMC49hAAAAAAmRkEImLyq2ltV+xcXF8cOBQAAYEsoCAGT190HunvfwsLC2KEAAABsiVELQn6VBwAAANh6oxaE/CoPAAAAsPVcZQwAANg051577ortBy89uMWRALCUOYQAAAAAJkZBCAAAAGBiFIQAAAAAJsYcQgCwgl1X3DB2CAAAsGn0EAIAmANVtbeq9i8uLo4dCgAwAQpCAABzoLsPdPe+hYWFsUMBACZAQQgAAABgYkYtCOkaDQAAALD1Ri0I6RoNAAAAsPVcZQyAbc3VwAAA4NiZQwgAAABgYhSEAAAAACbGkDEAjtuxDNc6dNWFmxgJAABwLPQQAgAAAJgYBSEAAACAiVEQAgAAAJgYBSEAAACAiVEQAiavqvZW1f7FxcWxQwEAANgSCkLA5HX3ge7et7CwMHYoAAAAW2LUgpBf5QEAAAC23qgFIb/KAwAAAGw9Q8YAAAAAJkZBCAAAYIerqidX1a9W1cur6pvGjgcYn4IQAADANlRVL6iqu6vqlmXtF1TV7VV1R1VdkSTd/bLuvizJ9yR56gjhAnNGQQgAAGB7uibJBUsbquqEJM9L8sQk5yS5pKrOWbLIjw+PAxN34tgBAAAwu/pqkr27d+8eO5Qkya4rblix/dDJT1t5hTPP2MRogJV092uratey5vOT3NHd70qSqnpxkouq6rYkVyV5ZXe/aWsjBeaRHkIAAHPA1VeBDXJqkjuX3D88tD07yROSfHtVXb7SilW1r6purqqb77nnns2PFBiVHkIAAAA7R63Q1t19dZKrj7Zid+9Psj9J9uzZ05sQGzBHFIQAACbi3GvPXbH94KUHtzgSYBMdTnL6kvunJblrpFiAOWbIGAAAwM5xU5KzqurMqjopycVJrh85JmAOKQgBAABsQ1V1XZLXJzm7qg5X1TO7+94kz0ryqiS3JXlJd986ZpzAfDJkDAAAYBvq7ktWab8xyY1bHA6wzSgIHYfVLsO6kkNXXbiJkQAAAAAcO0PGAAAAACZm1IJQVe2tqv2Li4tjhgEAAAAwKaMWhLr7QHfvW1hYGDMMAAAAgEkxZAwAAABgYrbdpNLHMqEzAAAAAPe17QpCAAAAbI6q2ptk7+7du8cOBXac1Tq4jHV1cgUhAACYN1euMsfmlS7Gwubq7gNJDuzZs+eysWMBNpc5hAAAAAAmRkEIAAAAYGIUhAAAAAAmRkEIAAAAYGIUhAAAAAAmRkEIAAAAYGIUhAAAAAAmRkEIAAAAYGJOHDsAgLFV1d4ke3fv3j12KABwVOdee+6K7QcvPbjFkQCw3ekhBExedx/o7n0LCwtjhwIAALAlFIQAAAAAJkZBCAAAgCSzofRVtX9xcXHsUIBNpiAEAABAEkPpYUpGLQipPgMAAABsvVELQqrPAAAAAFvPkDEAgDmg5zQAsJUUhAAA5oCe0wDAVlIQAgAAAJgYBSEAAACAiVEQAgAAAJgYBSEAAACAiVEQAgAAAJiYE8cOYF7suuKGsUMAAAAA2BJ6CAEAAABMjIIQAAAAwMQYMgbAZzCEdvMdy2t86KoLNzESAACmSg8hAAAAkiRVtbeq9i8uLo4dCrDJ9BACAICJO/fac1dsP3jpwS2OhLF194EkB/bs2XPZ2LEAm0sPIQAAAICJURACAAAAmBhDxgAAAOAoDKtkJ9JDCAAAAGBi9BACAICR7LrihhXbD528xYEAMDl6CAEAAABMjIIQAAAAwMQoCAEAAABMjIIQAAAAwMQoCAEAAABMjIIQAAAAwMQoCAEAAABMzKgFoaraW1X7FxcXxwwDAAAAYFJGLQh194Hu3rewsDBmGAAAAACTcuLYAQAAMB27rrhhxfZDV124xZEAwLSZQwgAAABgYhSEAAAASGKeV5gSBSEAAACSmOcVpkRBCAAAAGBiTCoNwJZYbSLZlZhcFgAANpeCEAAAACTJlasMlTvzjK2NA7aAIWMAAAAAE6MgBAAAADAxCkIAAAAAE6MgBAAAADAxCkIAAAAAE6MgBAAAADAxCkIAAAAAE6MgBAAAADAxCkIAAAAAE6MgBOxYVfWoqvq1qnrp2LEAAADMEwUhYFupqhdU1d1Vdcuy9guq6vaquqOqrkiS7n5Xdz9znEgBAADml4IQsN1ck+SCpQ1VdUKS5yV5YpJzklxSVedsfWgAAADbg4IQsK1092uTfGhZ8/lJ7hh6BH0iyYuTXLSW7VXVvqq6uapuvueeezY4WgAAgPmkIATsBKcmuXPJ/cNJTq2qh1XV85M8pqqes9KK3b2/u/d0955TTjllK2IFAAAY3YljBwCwAWqFtu7uDya5fKuDAQAAmHd6CAE7weEkpy+5f1qSu0aKBeAzuOIhADCPFISAneCmJGdV1ZlVdVKSi5NcP3JMwA7miofATlVVe6tq/+Li4tihAJtMQQjYVqrquiSvT3J2VR2uqmd2971JnpXkVUluS/KS7r51zDiBHe+auOIhsAN194Hu3rewsDB2KMAmM4cQsK109yWrtN+Y5MYtDgeYqO5+bVXtWtb8qSseJklVHbni4dvub3tVtS/JviQ544wzNjbYbe7ca89dsf3gpQe3OBKO2HXFDSu2Hzr5aSuvcKbPNMA80kMIAGBjuOIhALBt6CEETF5V7U2yd/fu3WOHAmxvrngIAGwbeggBk2esPLBBXPEQANg2FIQAADaGKx4CANuGghAAwDFyxUMAYLszhxAAwDFyxUMAYLvTQwgAAABgYhSEAADmQFXtrar9i4uLY4cCAEyAghAAwBxwxUMAYCspCAEAAABMjEmlgcmrqr1J9u7evXvsUBjsuuKGsUMAAIAdTUEImLzuPpDkwJ49ey4bOxYAgJ3u3GvPXbH94KUHtzgSmDZDxgAAAAAmRkEIAAAAYGI2vCBUVY+qql+rqpdu9LYBAAAAWL81FYSq6gVVdXdV3bKs/YKqur2q7qiqK5Kku9/V3c/cjGABAHaqqtpbVfsXFxfHDgUAmIC19hC6JskFSxuq6oQkz0vyxCTnJLmkqs7Z0OgAACaiuw90976FhYWxQwEAJmBNBaHufm2SDy1rPj/JHUOPoE8keXGSi9b6xFW1r6purqqb77nnnjUHDAAAAMD6rGcOoVOT3Lnk/uEkp1bVw6rq+UkeU1XPWW3l7t7f3Xu6e88pp5yyjjAAAAAAOBYnrmPdWqGtu/uDSS5fx3YBAAAA2ETr6SF0OMnpS+6fluSu9YUDsPVM5AoAAEzNegpCNyU5q6rOrKqTklyc5PqNCQtg65jIFQAAmJq1Xnb+uiSvT3J2VR2uqmd2971JnpXkVUluS/KS7r5180IFAABgM+k5DdOxpjmEuvuSVdpvTHLjhkYEAADAKLr7QJIDe/bsuWzsWIDNtZ4hY+um+gwAAACw9dZzlbF1U30GAJipqr1J9u7evXvsUGBbOPfac1dsP3jpwS2OBGB7GrWHEAAAMya4BwC2koIQAAAAwMQoCAEAAABMjIIQAAAAwMQoCAEAAABMjIIQAAAAwMSMetl5l1cF5oFcxDzbdcUNa1720FUXbmIkAADsJKP2EHJ5VWAeyEUAAMDUGDIGAAAAMDEKQgAAAAAToyAEAAAAMDGjTioNAMDM5Ce4v3KVedzOPGNr44Ak51577ortBy89uMWRAGwePYQAAOaACe4BgK2kIAQAAAAwMQpCAAAAABMzakGoqvZW1f7FxcUxwwAAAACYlFELQsbKAwAAAGw9Q8YAAAAAJkZBCAAAAGBiFIQAAAAAJubEsQMAGFtV7U2yd/fu3WteZ9cVN6x52UNXXXgcUQEAAGwePYSAyTPBPQAAMDUKQgAAAAAToyAEAAAAMDHmEAIAmAPHM58Z7CSrzc936OSnrbzCmWdsYjQAO9+oPYSqam9V7V9cXBwzDACA0ZnPDADYSqMWhBz4AAAAAGw9cwgBAAAATIyCEAAAAMDEKAgBAACQxDyvMCUKQgAAACQxzytMiYIQAAAAwMQoCAEAAABMjIIQAAAAwMQoCAEAAABMjIIQAAAAwMQoCAEAAABMzKgFoaraW1X7FxcXxwwDmDi5CAAAmJpRC0LdfaC79y0sLIwZBjBxchEAADA1howBAAAATIyCEAAAAMDEKAgBAMwB85kBAFtJQQgAYA6YzwwA2EoKQgAAAAAToyAEAAAAMDEKQgAAAAAToyAEAAAAMDEKQgAAAAAToyAEAAAAMDEKQgAAAAAToyAEAAAAMDEKQgAAAAATc+KYT15Ve5Ps3b1795hhbKpdV9yw5mUPXXXhJkYCAAAAMDNqD6HuPtDd+xYWFsYMAwAAAGBSRu0hBAAAAFtttZEch07e4kBgROYQAgAAAJgYBSEAAACAiVEQAgAAAJgYBSEAAACAiTGpNDB5VbU3yd7du3dvyvZXm7RwJYeuunBTYmAafNYAAFgrPYSAyevuA929b2FhYexQAAAAtoSCEAAAAMDEKAgBAMyBqtpbVfsXFxfHDgUAmAAFIQCAOWD4KgCwlUwqDQAAwHFb7aIGh05+2sornHnGJkYDrJUeQgAAAAAToyAEAAAAMDEKQgAAAAAToyAEAAAAMDEKQgAAAAAToyAEAAAAMDGjFoSqam9V7V9cXBwzDAAAAIBJGbUg1N0HunvfwsLCmGEAAAAATIohYwAAAAAToyAEAAAAMDEKQgAAACQxzytMiYIQAAAASczzClOiIAQAAAAwMQpCAAAAABOjIAQAAAAwMQpCAAAAABOjIAQAAAAwMQpCAAAAABOjIAQAAAAwMQpCAAAAABOjIAQAAAAwMQpCAAAAABOjIAQAAAAwMQpCAAAAABOjIAQAAAAwMSeOHQDA2Kpqb5K9u3fvHjuUY7LrihvWvOyhqy7cxEjg+PgMAwCMRw8hYPK6+0B371tYWBg7FAAAgC2hhxAAwBzYrr0VmU+r9cA7dPLTVl7hzDM2MZo5duUqPwZN9fUAJkUPIQCAOaC3IgCwlRSEAAAAACZGQQgAAABgYhSEAAAAACZGQQgAAABgYka9ypiraRy/1a4csZJDV124iZGsjXjny07fPwAAAI5u1B5CrqYBAAAAsPUMGQMAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIlREAIAAACYGAUhAAAAgIk5caM3WFUPSfLLST6R5NXd/aKNfg6A+yMXAfNALgIA5tWaeghV1Quq6u6qumVZ+wVVdXtV3VFVVwzN35bkpd19WZInbXC8wITJRcA8kIsAgJ1grUPGrklywdKGqjohyfOSPDHJOUkuqapzkpyW5M5hsb/dmDABkshFwHy4JnIRALDNrakg1N2vTfKhZc3nJ7mju9/V3Z9I8uIkFyU5nNnBz5q3D7AWchEwD+QiAGAnWM8cQqfm0794JbMDnq9KcnWSX6qqC5McWG3lqtqXZF+SnHHGGesIg/uz64ob1rzsoasu3JTt7mSb9fqyZjsqF/lesd1N+DO8o3IR7DSr5aZDJ29xIEzCqp835wLbwpTev/UUhGqFtu7ujyb53vtbubv3J9mfJHv27Ol1xAFMm1wEzAO5CADYVtbTdflwktOX3D8tyV3rCwfgmMlFwDyQiwCAbWU9BaGbkpxVVWdW1UlJLk5y/caEBbBmchEwD+QiAGBbWetl569L8vokZ1fV4ap6Znffm+RZSV6V5LYkL+nuWzcvVGDq5CJgHshFAMBOsKY5hLr7klXab0xy44ZGBLAKuQiYB3IRALATjHr506raW1X7FxcXxwwDAAAAYFJGLQh194Hu3rewsDBmGAAAo/NDGbCZqupRVfVrVfXSsWMB5sOoBSGAeeAkDJgHfigDjlVVvaCq7q6qW5a1X1BVt1fVHVV1RZJ097u6+5njRArMIwUhYPKchAEA29Q1SS5Y2lBVJyR5XpInJjknySVVdc7WhwbMOwUhAACAbai7X5vkQ8uaz09yx9Aj6BNJXpzkorVsr6r2VdXNVXXzPffcs8HRAvNGQQgAAGDnODXJnUvuH05yalU9rKqen+QxVfWclVbs7v3dvae795xyyilbESswojVddh4AAIBtoVZo6+7+YJLLtzoYYH7pIQQAALBzHE5y+pL7pyW5a6RYgDk2akHIlX0AAAA21E1JzqqqM6vqpCQXJ7l+5JiAOTRqQciVfQAAAI5PVV2X5PVJzq6qw1X1zO6+N8mzkrwqyW1JXtLdt44ZJzCfzCEEADAHqmpvkr27d+8eOxRgm+juS1ZpvzHJjVscDrDNVHePHUOq6p4k797ip314kg9s8XMeKzFuDDFujLXE+EXdve0uSXHkJCzJU5O8Y+Rw1mI7fF6Oh/3aXuZ5v7ZlLjpipOOi+zNv7/e8xZOIaa2mFJNcNF/m8bO30aawj8k09nMj93HVXDQXBaExVNXN3b1n7DiORowbQ4wbYzvEOBU79b2wX9vLTt0vVjZv7/e8xZOIaa3ExFim8D5PYR+TaeznVu2jq4wBAAAATIyCEAAAAMDETLkgtH/sANZAjBtDjBtjO8Q4FTv1vbBf28tO3S9WNm/v97zFk4hprcTEWKbwPk9hH5Np7OeW7ONk5xACAAAAmKop9xACAAAAmCQFIQAAAICJmURBqKpOr6r/WVW3VdWtVfUDQ/tDq+r3q+odw7+fP3KcJ1TVm6vqFfMY3xDT51XVS6vqz4bX86vnLc6q+tfD+3xLVV1XVSePHWNVvaCq7q6qW5a0rRpTVT2nqu6oqtur6ptHjPH/G97rP62q/1ZVnzdmjFNQVRcMr+kdVXXFCo9XVV09PP6nVfXYoX3FPDcPjnefljz+GblxXqxnv1bKpVsb/erWuV/3yb9bGz2bYd6+g/P4/ZmHz/6xHmuMGNOqxxZjxLPksR+uqq6qh29VPGyeqvqB4ft4a1X94NA2V+csx+NYv+fb8Xh9lX38juG9/Luq2rNs+W23j8n8nHtNoiCU5N4kP9TdX5rkcUn+VVWdk+SKJH/Y3Wcl+cPh/ph+IMltS+7PW3xJ8p+S/Pfu/pIkX55ZvHMTZ1WdmuT7k+zp7kcnOSHJxXMQ4zVJLljWtmJMw2fz4iRfNqzzy1V1wkgx/n6SR3f3P0jy9iTPGTnGHW14DZ+X5IlJzklyyfBaL/XEJGcNt31J/vPQvlqeG9U69+mI5blxdBuwXyvl0tGtZ7+Okn/Z/ubtOzhX3585+uxfkzUea4wc04rHFiPGk6o6Pck3JnnPFsbCJqmqRye5LMn5meWIb62qszL+92EjXJP5P6dYr2ty3328Jcm3JXnt0sZtvI/JnJx7TaIg1N3v6+43DX9/JLMDh1OTXJTk2mGxa5M8eZQAk1TVaUkuTPJfljTPTXxJUlWfm+TrkvxaknT3J7r7w5mzOJOcmORBVXVikgcnuSsjx9jdr03yoWXNq8V0UZIXd/ffdPefJ7kjs//QtjzG7v697r53uPuGJKeNGeMEnJ/kju5+V3d/IsmLM3utl7ooya/3zBuSfF5VfcFR8tzYjnufklVz4zw47v06Si6dB+t6v7Jy/mUbm7fv4Bx/f0b/7B/jscZoMR3l2GKUeAY/n+RHkrjazs7wpUne0N0fGz5rr0nyTzN/5yzHbDucU6zXKnnjtu6+fYXFt+U+JvNz7jWJgtBSVbUryWOS/O8kj+zu9yWzolGSR4wY2i9k9h/R3y1pm6f4kuRRSe5J8sKh6/h/qaqHZI7i7O73JvnZzH7heV+Sxe7+vXmKcYnVYjo1yZ1Lljuc+Tix/2dJXjn8Pa8xbndreV3vd5lleW5s692nX8h9c+M8WM9+rZZL58Fx79dR8i/b2y9kvr6Dc/f9mfPP/jwe/yy19NhiFFX1pCTv7e63jhkHG+qWJF9XVQ+rqgcn+ZYkp2f+vw/Ha7udU2yknbyPW3LuNamCUFV9dpLfSfKD3f1XY8dzRFV9a5K7u/uNY8dyP05M8tgk/7m7H5Pko5mzrpbDmNmLkpyZ5AuTPKSqvnvcqI5ZrdA26i9WVfVjmQ1JetGRphUW86va+q3ldT3qMnOY5457n+Y8N67nvZrnXLqe92sn5F+WmNPv4Nx9f3z2j88KxxZjxPDgJD+W5P8ZKwY2XnffluS5mQ2/+e9J3prZZ21qpnC8viP3cSvPvSZTEKqqB2R2kvSi7v7dofn9S4YlfEGSu0cK7x8meVJVHcqse/7XV9VvzlF8RxxOcri7j/Q6eGlmB2XzFOcTkvx5d9/T3Z9M8rtJvmbOYjxitZgOZ/YrxhGnZcRhF1V1aZJvTfJd3X0k8cxVjDvIWl7XVZdZJc+NbT37tFpunAfr2a/Vcuk8WM9+rZZ/2b7m8Ts4j9+fef7sz+Pxz2rHFmP4+5kV8t46fM5PS/Kmqvp7I8bEBujuX+vux3b312U2LOcdmdPvwwbYFucUm2TH7eNWn3tNoiBUVZXZWPPbuvs/Lnno+iSXDn9fmuTlWx1bknT3c7r7tO7eldlkUf+ju797XuI7orv/IsmdVXX20PQNSd6W+YrzPUkeV1UPHt73b8hsLpV5ivGI1WK6PsnFVfXAqjozs4lb/2SE+FJVFyT50SRP6u6PLXlobmLcYW5KclZVnVlVJ2WWD65ftsz1SZ5RM4/LbGjC+46S58Z23Pt0lNw4D9azX6vl0nlw3PuV1fMv29Q8fgfn9Pszz5/9uTv+OcqxxZbr7oPd/Yju3jV8zg8neezwOWMbq6pHDP+ekdlkxNdlDr8PG2Tuzyk20Y7ax1HOvbp7x9+SfG1mXar+NMlbhtu3JHlYZjOxv2P496FzEOvjk7xi+Hse4zsvyc3Da/myJJ8/b3Em+XdJ/iyz8cO/keSBY8eY2X9C70vyycwONp55tJgy6778ziS3J3niiDHekdl41SPfm+ePGeMUbkNuevvw2v7Y0HZ5ksuHvyuzq0C9M8nBzK5qs2qeG3t/1rNPy7bxqdw4L7f17NdKuXTs/dmg/bpP/h17f9w27HMxN9/Befz+zMNnf5X/x+fx+GfVY4sx4ln2+KEkDx/78+S2Ie/1H2VWLH5rkm8Y2ubqnOU492vuzyk2aR//6fD33yR5f5JXbed9PMp+bvm5Vw0bBwAAAGAiJjFkDAAAAIBPUxACAAAAmBgFIQAAAICJURACAAAAmBgFIQAAAICJURACAAAAmBgFIQAAAICJURBiS1XVL1bVm6rqK5e0fWlVPb+qXlpV/2LM+ICda6X8c5Rl5SVgbh1LPgOmq6p+sqoOVtXbq2rfOre1q6o+XlVvWdJ28ZCLfnC4/6CqektVfaKqHr6+6NkKCkJsmap6SJJHJPm+JN96pL27b+vuy5N8Z5I9I4UH7GCr5Z/VyEvAvDrWfAZMU1V9c5LHJDkvyVOSPHkDNvvO7j5vyf2Lk3xlksdV1Wd398eHx+/agOdiCygIcUyqandVHVzW9sCq+vOqOmdJ230qyN390SRfkOTVSa5eto0nJXldkj8c7qsuA8dspdyTrJ5/qurcqnrFstsjhsc+lZfkJGAtqurxVfUbG7StNeczOQqmoaq+rKr+YOjx82+H3oKr9RR8UpJrkjwgybOS/M5RtntWVR2qqt3D/QdU1Vur6rT7C2n4t5f8zTZy4tgBsO28K8npVfVZ3f13Q9u+JK/p7rctW/YzKshV9bAkD07ykSR/u3TB7r4+yfVVdUOS3+rujyc5r6oObc5uADvY8l+vVs0/3X0wq/zCvjQvdfdvRU4C7t95Sd68gdtbUz5z3AQ7X1WdnOS/JvmOzM7J/izJG7v7plVW+YokNyX5YJJDSf71atvu7ndU1f4k35zkjswKSC/v7sP3E9bvJrk5yW9290fWvjfMCz2EOCZDEeg9SXYls1+kkvxQkivXsPqPJ/nZJLcmWdqb6PFVdXVV/UqSGzc4ZGCHqap/WVW3VNW7q+rZa1xtxfxzlOeQl4AkSVV9e1W9Yfi1/HVVdcpRFv/yJKdW1f+uqndV1eOPst0vrKrfqao3V9WfVdX5awzpmPIZsGM8Icmbu/vWoQh8UpKfW2nBqvqsJKd19zVJHp7kjUn+r/vZ/i1Jzq6qhyZ5ZpLn3l9A3X1tdz+mu1eMg/mnIMTxuC3Jlwx//6sk13f3oaOtUFW7knxNkt8e1v+yI49196u7+/u7+/u6+3mbEjGwI1TVU5J8Y2Zj4h+X5Ceq6qi9XY+Wf1YjLwFL/M/uflx3f3mS389sbrHVnJfkI939VUkuT/LvV1poyFuvTPLC7n5Mksdmlp+O6njyGbBjPCbJm5JZQTnJX3f3/1pl2bOTvCP5VA/C/5XkhPvZ/tuH9a5M8rPD8FR2OAUhjsdtmVWPPzuzgtBPr2Gdn0ryk93dcQADHL/vT/Kj3f3J7n5fkk/m/v8vk3+A9fieqvqTqnprkn+Z5P+stNBQ5HlYkv8wNL0ls1/mV/LkJLd19yuSpLs/tsbhFvIZTNffJDkyp8/PZNZDKElSVX9YVacuWfYxSR5YVSdU1QOTPC3Jy46yfJK8M7Pi9PlJfn3jw2cemUOI43Fbkq9P8gNJXtTd7z/awlV1XpJvS/K1VfW8JCcnOXi0dQCWq6oHJPkH3f324f4XJPlAd3/iKOucF/kHOE5V9YzMTo6+vrv/uqpem9lQrZWck+SOJTnpsUneusqy5yV5wzHGcl7kM5iy30ry8qq6PcmvZFbw+YXMhoLtTvKhJcuel+RBmRV5PpDkl7v7rcmnhpMtXz7d/cmq+qskVyyZK5YdTkGI43FbkisyG8f62DUs/9wke7v7yBXEHpmNnXARmIZzkixU1aMymxzxZ5L84v2sI/8A63Fukj8eikFPyWy41sFk9gt7kmd093uHZb88yZnDr/EPSPITGSZxXWHZvxiWz/D4Kd19z/3EIp/BhA0TPH/F8vaqenSS3xmGhh3xmCRP7+5bVtjUOSssf8QDkrxmI+JlezBkjONxe2YHSPu7e/FoC1bV1yd5yJGDlyQZehQ9ZJiwDGCtHpPkRUmuS/KnSd7T3ftXW1j+ATbAtUm+v6r+KMkXJ3lXd390lV/YvzyzHPXHSf4kydXd/YZVlr0mySOr6tbhsvJffbQg5DNgNd19S3cvnzD6SzK7Ctlalz8yR9m7hyGpK/nbzH6Ye8tqsVTVg4bHH5BEL6NtoFZ/v+H4DQnlFd396HVu51CSPd39gY2IC9i+hm7Rr+/u3z7KMruyAbnnKNs/FDkJJm/4Rf6frXRStZ5ll623K8eYz+QoAI6FHkJslvutIB+N6jKwgvMym6T1aNaVe1YjJwFLrfYL+3qXXWbN+UyOAuB46CEEAAAAMDF6CAEAAABMjIIQAAAAwMQoCAEAAABMjIIQAAAAwMQoCAEAAABMjIIQAAAAwMQoCAEAAABMzP8PUC52vwGAreUAAAAASUVORK5CYII=\n",
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
"text/plain": [
"<Figure size 1440x504 with 4 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(20,7))\n",
"tc.plot.cell()"
]
},
{
"cell_type": "markdown",
"id": "4cb440b6-aa79-4369-80d6-141dcfc37696",
"metadata": {},
"source": [
"### Create an atomicrex job"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "56951cee-038f-44e7-bb11-b20841036f81",
"metadata": {},
"outputs": [],
"source": [
"job = pr.create.job.Atomicrex(\"FitAl\", delete_existing_job=True)"
]
},
{
"cell_type": "markdown",
"id": "9807a140-dc20-4c4a-b22c-dcc89c789425",
"metadata": {},
"source": [
"### Training Data\n",
"Set the training data. This can also be done structure by structure to set other fit properties and weights, but here we simply load the TrainingContainer"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "adbcf0a6-1ef8-4697-93d5-bd9346aad1ab",
"metadata": {},
"outputs": [],
"source": [
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
]
},
{
"cell_type": "markdown",
"id": "1d275454-0b14-4862-b4f0-8b7166852f96",
"metadata": {},
"source": [
"Set the potential type. In this case an EAM potential"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "e90476e0-54cf-4a6e-9424-bf282517a4b5",
"metadata": {},
"outputs": [],
"source": [
"job.potential = job.factories.potentials.eam_potential()"
]
},
{
"cell_type": "markdown",
"id": "1cd689cb-0d35-4b5c-93e6-e0112a5382b2",
"metadata": {},
"source": [
"### Functions\n",
"Reminder: $E = \\sum_{ij}V(r_{ij}) + \\sum_i F(\\rho_i)$ with $\\rho_i = \\sum_j \\rho(r_{ij})$\n",
"\n",
"It is necessary to define a pair potential, an electronic density function and an embedding function.\n",
"For all of those it is possible to choose between different functional forms.\n",
"Classic pair potentials are physically motivated and have a very limited number of paramaters that can often be derived from an experimentally measurable quantity.\n",
"Splines or polynomials offer more flexibility, but can lead to unphysical oscillations or overfitting. Compared with the machine learning potentials shown later the number of parameters is very low no matter which functions you choose.\n",
"\n",
"In this case a generalized morse function is used for the pair interaction. It has the form\n",
"\n",
"$(\\frac{D_0}{S-1}exp(-\\beta \\sqrt{2S}(r-r_0))-\\frac{D_0S}{S-1}exp(-\\beta\\sqrt{2/S}(r-r_0)))+\\delta $\n",
"\n",
"The parameters in the morse potential can be derived from phyiscal quantities, 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.8 $\\mathring A$ so it is taken as initial value.\n",
"In the case of analytic functions the initial parameter choices should not matter too much, since the functional form is constrained.\n",
"\n",
"The electronic density will be a cubic spline. The embedding function will be $-A*sqrt(\\rho)+B*rho$, which can be defined as a user function.\n",
"\n",
"The pair function and the electron denity and their first derivatives are required to smoothly approach 0 at the cutoff distance $r_{c}$\n",
"For this purpose the pair function is screened by multiplying with the function:\n",
"\n",
"$\\Psi(\\frac{r-rc}{h})$ where $\\Psi(x) = \\frac{x^4}{1+x^4}$ if $x<0$ else $\\Psi(x)=0$\n",
"\n",
"For the spline it is necessary to set a node point with y value 0 at the cutoff "
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "96a1900f-f4d2-43b6-87dd-cf5570517836",
"metadata": {},
"outputs": [],
"source": [
"morseB = job.factories.functions.morse_B(\"V\", D0=0.15, r0=3.05, beta=1.1, S=4.1, delta=-0.01)\n",
"morseB.screening = job.factories.functions.x_pow_n_cutoff(identifier=\"V_screen\", cutoff=7.6)\n",
"morseB.parameters.D0.min_val = 0.05\n",
"morseB.parameters.D0.max_val = 1.55\n",
"\n",
"morseB.parameters.r0.min_val = 2.6\n",
"morseB.parameters.r0.max_val = 3.1\n",
"\n",
"morseB.parameters.S.min_val = 1.5\n",
"morseB.parameters.S.max_val = 4.5\n",
"morseB.parameters.delta.max_val = 0.005"
]
},
{
"cell_type": "markdown",
"id": "252c78d0-9230-4e23-83a1-20304f1630b5",
"metadata": {},
"source": [
"It is also possible to plot most of the functions. This can help to judge if the initial parameters are reasonable"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "1c380b21-78e8-404c-a9ae-a5b6d6ea85fd",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(<Figure size 720x504 with 1 Axes>,\n",
" <AxesSubplot:xlabel='r [$\\\\AA$]', ylabel='func(r)'>)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAGzCAYAAABjHhDPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkQklEQVR4nO3de3Scd33n8c93ZqTR1ZJv8kW+xUmwIbZJYgU2CVAwaTcEGrq0oVC23aVp3XZpF5bdbpfSLdvbKV0Kp7d0u1nu3VCg0Gxowr3hWpIQ2yRxHMdu4sSOYyeW7diWZVvSzHz3j3kkj21FHluX7zwz79c5Oppn9EjzZeCEd57Lb8zdBQAAgBiZ6AEAAAAaGTEGAAAQiBgDAAAIRIwBAAAEIsYAAAACEWMAAACBwmLMzFrM7Idm9rCZbTOz34uaBQAAIIpFrTNmZiap3d2Pm1mTpO9Lere73x8yEAAAQIBc1At7uQKPJ5tNyRcr0AIAgIYSFmOSZGZZSZslXSbpNnd/YJx9NkraKEnt7e3rV69ePbNDoqY9f+yUDgwMaW1vV/QoAACcYfPmzQfdff759gs7TXnGEGbdku6U9Bvu/uiL7dfX1+ebNm2asblQ+z7zwB799p1bdd/7NmhRV2v0OAAAjDGzze7ed779auJuSnc/Iunbkm6MnQRps7i7RZK078jJ4EkAALg4kXdTzk+OiMnMWiXdIOnxqHmQTr3d5aNhzx45FTwJAAAXJ/KasUWSPpVcN5aR9Hl3vztwHqTQoiTGODIGAEiryLspH5F0VdTroz505HPqam0ixgAAqVUT14wBk7G4u5UYAwCkFjGG1OvtbuGaMQBAahFjSD2OjAEA0owYQ+ot7m7V0ZMjOj5UiB4FAIALRowh9RYnd1Tu5+gYACCFiDGkXm+y8OuzxBgAIIWIMaTe4rGFX4kxAED6EGNIvZ7OFmUzxkX8AIBUIsaQetmMaeGsFu1jeQsAQAoRY6gLvd2tevYFjowBANKHGENd6J3dyjVjAIBUIsZQF5bNadO+oyc1XChFjwIAwAUhxlAXls1pkzt3VAIA0ocYQ11YNrdNkrT70GDwJAAAXBhiDHVh2ZxyjD1z+ETwJAAAXBhiDHVhfkde+VxGe4gxAEDKEGOoC5mMaemcNmIMAJA6xBjqxvI5bdp9iBgDAKQLMYa6sXROm545fELuHj0KAABVI8ZQN5bNadPgcFGHB4ejRwEAoGrEGOrG6B2VXDcGAEgTYgx1Y3StMWIMAJAmxBjqxtLZSYxxET8AIEWIMdSN1uasejrzHBkDAKQKMYa6soy1xgAAKUOMoa4sS5a3AAAgLYgx1JWlc9q0/9gpDRWK0aMAAFAVYgx1ZfncNrlLe184GT0KAABVIcZQV1hrDACQNsQY6spojHHdGAAgLYgx1JX5nXnlcxk+MBwAkBrEGOqKmbG8BQAgVYgx1B2WtwAApAkxhrqzbG75yJi7R48CAMB5EWOoOyvmtuvEcFEHBoaiRwEA4LyIMdSdlfPbJUlP9h8PngQAgPMjxlB3Vs7vkCTt6h8MngQAgPMjxlB3Fs1qUUtThhgDAKQCMYa6k8mYLpnXoV0HOU0JAKh9xBjq0sr57RwZAwCkAjGGunTpvHbtfeGEhgrF6FEAAJgQMYa6tHJ+h0ouPhYJAFDziDHUpdHlLXaxvAUAoMYRY6hLl8wbXWuM68YAALWNGENd6mxpUk9nXk8dJMYAALWNGEPdKt9RyWlKAEBtI8ZQt1bO79AujowBAGocMYa6tXJeu46cGNHhweHoUQAAeFHEGOrWpWOfUcmpSgBA7SLGULdOL2/BqUoAQO0ixlC3lsxuU3M2oyf5jEoAQA0jxlC3shnT8rltHBkDANQ0Ygx1jeUtAAC1jhhDXVs5v0N7Dp9QoViKHgUAgHERY6hrK+e1a6To2nOYDwwHANSmsBgzs6Vm9i0z225m28zs3VGzoH5dvqBTkvQvBzhVCQCoTZFHxgqS/rO7v1TSv5L0LjN7WeA8qEOX95TXGtv53EDwJAAAjC8sxtx9v7tvSR4PSNouqTdqHtSn9nxOy+a06fHniTEAQG2qiWvGzGyFpKskPTDOzzaa2SYz29Tf3z/jsyH9XrKgkyNjAICaFR5jZtYh6YuS3uPux87+ubvf7u597t43f/78mR8QqbdqYYeeOjiooUIxehQAAM4RGmNm1qRyiN3h7v8QOQvq16qFs1QoOYu/AgBqUuTdlCbpY5K2u/tHouZA/VuV3FG5k+vGAAA1KPLI2PWSfl7SBjN7KPm6KXAe1KlL5rWrKWt6nOvGAAA1KBf1wu7+fUkW9fpoHM25jFbO6+AifgBATQq/gB+YCasWdmoHpykBADWIGENDWLWwU3tfOKnjQ4XoUQAAOAMxhobwEi7iBwDUKGIMDWH1wiTGuG4MAFBjiDE0hN7uVrU1Z7mjEgBQc4gxNIRMxsofi8RpSgBAjSHG0DBWLejUDo6MAQBqDDGGhvGShZ06NDisg8eHokcBAGAMMYaGwUX8AIBaRIyhYYzG2GP7jwVPAgDAacQYGsbcjrwWdbVo67NHo0cBAGAMMYaGsqa3S48SYwCAGkKMoaGsWdylXQcH+VgkAEDNIMbQUNYumSV36bF9XDcGAKgNxBgayprFXZLEqUoAQM0gxtBQema1qKczT4wBAGoGMYaGs7a3izsqAQA1gxhDw7mit0tP9h/XiWEu4gcAxCPG0HDW9nap5NJ2Fn8FANQAYgwNZ21v+SL+rXs5VQkAiEeMoeEsmJXXvI5mbX2WI2MAgHjEGBqOmWlNb5e27ePIGAAgHjGGhrS2t0v/cuC4To0Uo0cBADQ4YgwN6YrFXSqWXI9xET8AIBgxhoa0dkn5Iv5trDcGAAhGjKEhLe5q0dz2Zj30DDEGAIhFjKEhmZmuWjZbP9rzQvQoAIAGR4yhYa1fPlu7Dg7q8OBw9CgAgAZGjKFhrV8+W5K0ZTdHxwAAcYgxNKx1S7qUy5g2c6oSABCIGEPDamnK6orFszgyBgAIRYyhoV29fLYe3ntEI8VS9CgAgAZFjKGhrV8+W6dGStrO4q8AgCDEGBra1cvKF/Fv5lQlACAIMYaGtri7VYu6WogxAEAYYgwN7+rls/WjPUeixwAANChiDA1v/bLZevbISe0/ejJ6FABAAyLG0PCuHlv89UjsIACAhkSMoeG9bNEs5XMZrhsDAIQgxtDwmnMZvXxptx58+nD0KACABkSMAZKuXTlX2/Yd1dETI9GjAAAaDDEGSLru0rkqufTAU4eiRwEANBhiDJB05bJutTRl9IMniTEAwMwixgBJ+VxW16yYo/uIMQDADCPGgMS1l87VjucHdPD4UPQoAIAGQowBiesunSdJun8XR8cAADOHGAMSaxbPUmc+x3VjAIAZRYwBiVw2o1eu5LoxAMDMIsaACtdeOk9PHRzUviN8TiUAYGYQY0CF6y6dK0kcHQMAzBhiDKiwakGn5rQ3c90YAGDGEGNAhUzGdO3KufrBkwfl7tHjAAAaADEGnOU1L5mn/UdPacfzA9GjAAAaADEGnOV1q3okSf+0/UDwJACARhAaY2b2cTM7YGaPRs4BVOqZ1aI1vbP0rceJMQDA9Is+MvZJSTcGzwCcY8PqBdqy5wW9MDgcPQoAoM6Fxpi7f1fS4cgZgPFsWN2jkkvf2dkfPQoAoM5FHxk7LzPbaGabzGxTfz//x4iZsa63S/M6mnUvpyoBANOs5mPM3W939z5375s/f370OGgQmYzptat69J2d/SoUS9HjAADqWM3HGBBlw+oeHT05oi17jkSPAgCoY8QY8CJeffk85TLGqUoAwLSKXtri7yTdJ2mVme01s1sj5wEqdbY06RWXzGGJCwDAtIq+m/Lt7r7I3ZvcfYm7fyxyHuBsG1b3aMfzA9pz6ET0KACAOsVpSmACN65ZKEm6Z+v+4EkAAPWKGAMmsGR2m65c2q17tu6LHgUAUKeIMeA83rRukR599piePjgYPQoAoA4RY8B53LR2kSROVQIApgcxBpzH4u5WrV8+W3c/QowBAKYeMQZU4Y1rF2n7/mN6sv949CgAgDpDjAFVGD1V+WWOjgEAphgxBlRhYVeLrlkxm+vGAABTjhgDqvTGtYv0+HMDeuLAQPQoAIA6QowBVbpp3SJlM6YvbH42ehQAQB0hxoAq9XS26HWrevTFLXtVKJaixwEA1AliDLgAb+1bov6BIX17R3/0KACAOkGMARfgdat7NK8jr89teiZ6FABAnSDGgAvQlM3op9f36t7HD+jAwKnocQAAdYAYAy7QLeuXqlhy3bmFC/kBAJNHjAEX6LKeDvUtn63PbXpG7h49DgAg5Ygx4CK8tW+pdvUPavPuF6JHAQCkHDEGXIQ3rlukjnxOf3v/7uhRAAApR4wBF6E9n9Nb+5bqnkf267mjXMgPALh4xBhwkd55/QqV3PWp+56OHgUAkGLEGHCRls5p00+8bKE+88AenRguRI8DAEgpYgyYhFtffYmOnhzRF1nmAgBwkYgxYBL6ls/WuiVd+sT3n1KpxDIXAIALR4wBk2BmuvVVl2jXwUF9e+eB6HEAAClEjAGTdNPaRVo4q0V/851dLAILALhgxBgwSU3ZjH7ttZfqh08d1n1PHooeBwCQMsQYMAV+9pqlWjirRR/5xk6OjgEALggxBkyBlqas3rXhMm3a/YK+/8TB6HEAAClCjAFT5K19S9Tb3crRMQDABSHGgCmSz2X1rtddph/tOaLv7OyPHgcAkBLEGDCFfmb9Ei2Z3aoPf30n644BAKpCjAFTqDmX0X/5iVXa+uxRfWHL3uhxAAApQIwBU+zNVy7W1cu69T+/ukMDp0aixwEA1DhiDJhiZqb/cfMVOjQ4pL+894nocQAANY4YA6bBuiXdumX9En3in5/Srv7j0eMAAGoYMQZMk9/816uVz2X1B3c/xlIXAIAXRYwB02R+Z17vueFyfWtHv/7xkf3R4wAAahQxBkyjf3/dCr18abc+cNej6h8Yih4HAFCDiDFgGuWyGX34lnUaHC7qd/7fVk5XAgDOQYwB0+yynk6998dfoq9te15fenhf9DgAgBpDjAEz4JdfvVJXLu3WB760Tc8dPRU9DgCghhBjwAzIZkx/esvLNVwo6dc/s0UjxVL0SACAGkGMATPksp4O/fFb1mrT7hf0J195PHocAECNIMaAGfTmK3v1C9cu10e//5S+spXlLgAAVcSYmV1rZreZ2SNm1m9me8zsy2b2LjPrmokhgXry/je+VC9f2q3f/MIjeuIAq/MDQKObMMbM7CuSfknS1yTdKGmRpJdJ+h1JLZLuMrObp3tIoJ7kc1n99TuuVj6X0Ts/+UMdGOCCfgBoZDbRukdmNs/dD074B6rYZ6r09fX5pk2bZuKlgGn38DNH9Lbb79elPe367MZr1ZHPRY8EAJhCZrbZ3fvOt9+ER8bc/aCZZc3smxPtczEDAo3u5Uu7dds7rtL2/QP6D3dwhyUANKrzXjPm7kVJJ7g+DJh6G1Yv0B/91Bp9d2e/3vv5h1UgyACg4VR7XuSUpK1m9g1Jg6NPuvt/nJapgAbytlcs05GTI/rgVx5XoVjSn7/tKjXnuNEZABpFtTF2T/IFYBr86o9dqlzG9If3bNfIHZt12zuuVj6XjR4LADADqooxd//UdA8CNLpfevVK5XMZ/fe7tunWT27Sbe+4Wl2tTdFjAQCm2fmWtvhHM/tJMzvn/xHMbKWZ/b6Z/eL0jQc0lp+/doU+9DPrdP+uQ3rLX/+zdh8aPP8vAQBS7XwXpvyypFdLetzMHkwWe73XzJ6S9L8lbXb3j0/7lEADuaVvqf721lfq0OCwfuq2f9YDuw5FjwQAmEYTrjN2xo5mK1Re9PWkpJ3ufmLSL252o6Q/l5SV9FF3/+BE+7POGBrJ0wcHdeunHtTTh07oP91wuX7ttZcpm7HosQAAVZqSdcYq/tglkp5z9/vc/SFJpSTOJjNgVtJtkt6g8qr+bzezl03mbwL1ZMW8dt35rut109pF+tOv79Q7Pnq/9h89GT0WAGCKVXv//N9LqlwAqZQ8NxmvkPSEu+9y92FJn5X05kn+TaCuzGpp0l+87Up96GfW6ZG9R3Xjn31Pn3twj0ql6o5oAwBqX7UxlkuCSZKUPG6e5Gv3SnqmYntv8hyACmamW/qW6u7feJVesqBDv/XFrfrZ2+/TzucHokcDAEyBatcZ6zezm939S5JkZm+WNNmPQRrv4pdz/nXfzDZK2ihJy5Ytm+RLAum1cn6HPrfxWn1hy1798Ze366Y//55+7pXL9OsbLlNPZ0v0eABSrFRyFd1VLLncNfZ49PlSyVVKni+Vyj8rusvdVSypvK9f3PMlH3189hyjv3fu8+5SyaVS8jdHH5dcciU/L/k5+7iS7+4qlSp+x12u09tjv5Ps49K5r5PM7ec8d3q7WtXG2K9KusPM/krliHpG0i9c4H/XZ9sraWnF9hJJ+87eyd1vl3S7VL6Af5KvCaRaJmN6a99S3fDSBfrw13foMw/s0d9v2qt3Xr9CG1+zUt1tkz1gDWA87q6RomukWEq+Kh+fuT1ccBVKpx+PFEvl7YJrpFRSseQqFMthUSi5iqXy71duF0qntwvFM7eLxeT5c/5WqeJ3ztyu/D5SPPP5YsovezCTMmbKmGSyM7YzlmxnTKbR7dM/y1j57MO4v2M29ni83xnveyaTOeN3qv7PUO3dlOX/wNaR/M6kz4+YWU7STkmvl/SspAcl/Zy7b3ux3+FuSuBMTx8c1Ee+sVNfenifWpuyuqVviW591SVaPrc9ejTgohVLruFCScOFkoYKRQ0VSslXMXmu/DX285GShoslDY0Uk+/Jz1/kuZFiSYWijz0+ZzuJqeHC6cgqzECwZEzKZTPKZUzZjCXfy9u5rFU8nyl/z5a3m87artzvzOdO/+2mbHk7m5GyZspk7PT3zOnwyI5tW7J9/ufH/paZMsnfH/f5itfM2OjjiZ8ffb3KaLqQ6Jlp1d5NWVWMmVle0k9LWqGKo2nu/vuTmFFmdpOkP1N5aYuPu/sfTbQ/MQaMb8dzA/o/39ulux56VoWSa8OqHr31mqXasLpHTVk+5xKT4+4aKpR0aqSokyNFnRwufz81kjw3tl0c2+fUSGls38rnT46UdGq4qFOF8nNDhdIZMTVUmJrwMZPyuYzyuazyuYyak698LqumrKkpm6n4fvpxc7YcMGc/35Qt/34uk2znMmrOloNn9HFTNqNc8jvNye/kznrclMRWLpNRtiKwRuMD9WWqY+yrko5K2iypOPq8u394MkNeKGIMmNiBY6f06ft26/ObntGBgSHNbW/WzVcu1hvWLNL65bNZp6xBuLtODBd1fKhQ/jpV0OBQQQND5e/jPX/8VEGDw+Xvx4cK5YgqlJLvRV3ASZQxTVlTS1NWrU1ZtTZn1ZLLqqU5q9amjFqbsuVQasqMRVM5lk5H0+nHGeWbsmrOZsr7j36f4HdyGavpIyZoDFMdY4+6+5opmWwSiDGgOoViSd/9l359/sG9unfHAQ0XSprX0awbXrpA1102T9eunKv5nfnoMVHB3XVypDgWQ2cE01gkFXV8aOSMx4NDxdORlcTV8eFCVfGUzZg68rnTXy05tedz6shn1dacU2tTVi1JOJUjKjsWVy0VP2ttrnxuNLzKR4mARlZtjFV7Af8PzGytu2+d5FwAZkAum9GG1Qu0YfUCHR8q6Ns7Duirjz6nex7Zr88+WF5R5vKeDl136Vxdc8kcXbG4S8vntHGa5AK5u06NlDSQRFFlSJ1xxGmcuBo46/nBoUJVd19lM6b25qw6W5rUns+qI5/TrJacertb1N5cDqrRuGrP59TZknvR5/O5DEePgBpQ7ZGxxyRdJukpSUMq31Hp7r5uesc7E0fGgMkpFEvatu+YfvDkId2365AefOqwTo6Urzxob87qpYtm6YrFs3TZgk4tn9OmFXPbtbi7pa6OcIxe/3R2DJ19RGpw6NxgOl5xBGp0u5qAypjKAZSE0NlhdPaRqc5xn8+qM9+kliYCCkiLqT5NuXy8591990XMdtGIMWBqDRdK2vn8gB7bd0zb9h3Vtn3HtH3/MQ0Oj10aqlzG1Du7VQtmtainM6/5yde8jrxmteTU2dI0FgydLTl15puUz2UmdZRtorvphpLrmAaHykeYBoeKOlH5fbioE0Pl03iV25XXUFVzK7+Z1NFceepunCNO+aw68k3qyGfL+511BGr0fWltyhJQQAOa6tOU6V6EBMC4mnMZrent0preLo0u+1cquQ4MDGn3oUHtPnRCuw+Xvx8YGNJj+47pwMCQjg8Vzvu3R2+nr7w7LZuxisURz1poMVnvaDJ307U2ZdWeXO/Uns+pvTmrrtYmLe5qSZ7LjgXSGWFVccRq9HFrU5bTtgBmRLUxdo/KQWaSWiRdImmHpCumaS4AQTIZ08KuFi3satErV84dd58TwwUdOj6sgVMFDZwaGTvidCzZHi6U122qXAyzUCo/l82cXhtodBHG0UUTsxl7kTvrKu6my2XU1pwtx1YSXG1JPHG3KIA0qirG3H1t5baZXS3pV6ZlIgA1r605p7Y51f67HABgIhd1Va67b5F0zRTPAgAA0HCq+ldbM3tvxWZG0npJ/dMyEQAAQAOZ8MiYmf1t8vB3JXUmX3lJd0t68/SOBgAAUP/Od2RsfbKsxR5Jf3nWz9oknZqWqQAAABrE+WLsbyR9VeW7JysX+DKV765cOU1zAQAANIQJT1O6+1+4+0slfcLdV1Z8XeLuhBgAAMAkVXU3pbv/2nQPAgAA0Ijq5wPnAAAAUogYAwAACESMAQAABCLGAAAAAhFjAAAAgYgxAACAQMQYAABAIGIMAAAgEDEGAAAQiBgDAAAIRIwBAAAEIsYAAAACEWMAAACBiDEAAIBAxBgAAEAgYgwAACAQMQYAABCIGAMAAAhEjAEAAAQixgAAAAIRYwAAAIGIMQAAgEDEGAAAQCBiDAAAIBAxBgAAEIgYAwAACESMAQAABCLGAAAAAhFjAAAAgYgxAACAQMQYAABAIGIMAAAgEDEGAAAQiBgDAAAIRIwBAAAEIsYAAAACEWMAAACBiDEAAIBAxBgAAEAgYgwAACAQMQYAABAoJMbM7BYz22ZmJTPri5gBAACgFkQdGXtU0lskfTfo9QEAAGpCLuJF3X27JJlZxMsDAADUjJq/ZszMNprZJjPb1N/fHz0OAADAlJq2I2Nm9k1JC8f50fvd/a5q/4673y7pdknq6+vzKRoPAACgJkxbjLn7DdP1twEAAOpFzZ+mBAAAqGdRS1v8GzPbK+laSfeY2dci5gAAAIgWdTflnZLujHhtAACAWsJpSgAAgEDEGAAAQCBiDAAAIBAxBgAAEIgYAwAACESMAQAABCLGAAAAAhFjAAAAgYgxAACAQMQYAABAIGIMAAAgEDEGAAAQiBgDAAAIRIwBAAAEIsYAAAACEWMAAACBiDEAAIBAxBgAAEAgYgwAACAQMQYAABCIGAMAAAhEjAEAAAQixgAAAAIRYwAAAIGIMQAAgEDEGAAAQCBiDAAAIBAxBgAAEIgYAwAACESMAQAABCLGAAAAAhFjAAAAgYgxAACAQMQYAABAIGIMAAAgEDEGAAAQiBgDAAAIRIwBAAAEIsYAAAACEWMAAACBiDEAAIBAxBgAAEAgYgwAACAQMQYAABCIGAMAAAhEjAEAAAQixgAAAAIRYwAAAIGIMQAAgEDEGAAAQCBiDAAAIBAxBgAAEIgYAwAACESMAQAABAqJMTP7kJk9bmaPmNmdZtYdMQcAAEC0qCNj35C0xt3XSdop6X1BcwAAAIQKiTF3/7q7F5LN+yUtiZgDAAAgWi1cM/aLkr4SPQQAAECE3HT9YTP7pqSF4/zo/e5+V7LP+yUVJN0xwd/ZKGmjJC1btmwaJgUAAIgzbTHm7jdM9HMz+3eS3iTp9e7uE/yd2yXdLkl9fX0vuh8AAEAaTVuMTcTMbpT0W5J+zN1PRMwAAABQC6KuGfsrSZ2SvmFmD5nZ3wTNAQAAECrkyJi7XxbxugAAALWmFu6mBAAAaFjEGAAAQCBiDAAAIBAxBgAAEIgYAwAACESMAQAABCLGAAAAAhFjAAAAgYgxAACAQMQYAABAIGIMAAAgEDEGAAAQiBgDAAAIRIwBAAAEIsYAAAACEWMAAACBiDEAAIBAxBgAAEAgYgwAACAQMQYAABCIGAMAAAhEjAEAAAQixgAAAAIRYwAAAIGIMQAAgEDEGAAAQCBiDAAAIBAxBgAAEIgYAwAACESMAQAABCLGAAAAAhFjAAAAgYgxAACAQMQYAABAIGIMAAAgEDEGAAAQiBgDAAAIRIwBAAAEIsYAAAACEWMAAACBiDEAAIBAxBgAAEAgYgwAACAQMQYAABCIGAMAAAhEjAEAAAQixgAAAAIRYwAAAIGIMQAAgEDEGAAAQCBiDAAAIBAxBgAAEIgYAwAACESMAQAABCLGAAAAAoXEmJn9gZk9YmYPmdnXzWxxxBwAAADRoo6Mfcjd17n7lZLulvS7QXMAAACECokxdz9WsdkuySPmAAAAiJaLemEz+yNJvyDpqKTXTbDfRkkbk80hM3t0BsarV/MkHYweIqV47yaH929yeP8mh/fv4vHeTc6qanYy9+k5KGVm35S0cJwfvd/d76rY732SWtz9A1X8zU3u3jeFYzYU3r+Lx3s3Obx/k8P7Nzm8fxeP925yqn3/pu3ImLvfUOWun5F0j6TzxhgAAEC9ibqb8vKKzZslPR4xBwAAQLSoa8Y+aGarJJUk7Zb0q1X+3u3TN1JD4P27eLx3k8P7Nzm8f5PD+3fxeO8mp6r3b9quGQMAAMD5sQI/AABAIGIMAAAgUCpizMxuNLMdZvaEmf236HnSxsw+bmYHWKPtwpnZUjP7lpltN7NtZvbu6JnSxMxazOyHZvZw8v79XvRMaWNmWTP7kZndHT1L2pjZ02a2NfnovU3R86SNmXWb2RfM7PHkn4HXRs+UFma2Kvnf3ejXMTN7z4vuX+vXjJlZVtJOST8uaa+kByW93d0fCx0sRczsNZKOS/q0u6+JnidNzGyRpEXuvsXMOiVtlvRT/O+vOmZmktrd/biZNUn6vqR3u/v9waOlhpm9V1KfpFnu/qboedLEzJ6W1OfuLFp6EczsU5K+5+4fNbNmSW3ufiR4rNRJOuZZSa90993j7ZOGI2OvkPSEu+9y92FJn5X05uCZUsXdvyvpcPQcaeTu+919S/J4QNJ2Sb2xU6WHlx1PNpuSr9r+N8AaYmZLJL1R0kejZ0FjMbNZkl4j6WOS5O7DhNhFe72kJ18sxKR0xFivpGcqtveK/zNEADNbIekqSQ8Ej5IqyWm2hyQdkPQNd+f9q96fSfqvKi8DhAvnkr5uZpuTj9ZD9VZK6pf0ieQ0+UfNrD16qJR6m6S/m2iHNMSYjfMc/2aNGWVmHZK+KOk9Z33QPc7D3YvufqWkJZJeYWacKq+Cmb1J0gF33xw9S4pd7+5XS3qDpHcll2ygOjlJV0v6X+5+laRBSVyzfYGS07s3S/r7ifZLQ4ztlbS0YnuJpH1Bs6ABJdc6fVHSHe7+D9HzpFVyiuPbkm6MnSQ1rpd0c3Ld02clbTCz/xs7Urq4+77k+wFJd6p82Quqs1fS3ooj2V9QOc5wYd4gaYu7Pz/RTmmIsQclXW5mlySF+TZJXwqeCQ0iuQD9Y5K2u/tHoudJGzObb2bdyeNWSTeIjz+riru/z92XuPsKlf+5d6+7/9vgsVLDzNqTm26UnF77CUncUV4ld39O0jPJp+VI5eueuHHpwr1d5zlFKcV9HFLV3L1gZr8u6WuSspI+7u7bgsdKFTP7O0mvlTTPzPZK+oC7fyx2qtS4XtLPS9qaXPckSb/t7l+OGylVFkn6VHI3UUbS592dJRowExZIurP871PKSfqMu381dqTU+Q1JdyQHQnZJemfwPKliZm0qrwTxK+fdt9aXtgAAAKhnaThNCQAAULeIMQAAgEDEGAAAQCBiDAAAIBAxBgAAEIgYAwAACESMAQAABCLGAKCCmf2lmW0xs2uiZwHQGIgxAEgkH5vTo/KK2W8KHgdAgyDGADQcM1thZicrPuJKkuTugyp/hNO3Jf1Fsm+rmT1kZsNmNm/GhwVQ94gxAHXPys7+592T7n7lWfvNldQmaUBSUZLc/WSy374ZGBVAAyLGANSl5OjXdjP7a0lbJC2t4td+R9KfStom6WXTOR8AjCLGANSzVZI+7e5XufvuiXY0sxWSrpP0OUnbJV0x/eMBADEGoL7tdvf7q9z3DyX9vru7iDEAMygXPQAATKPBanYysyslvUXSq8zsNkktkrZO41wAMIYYAwDpTyT9pLv/kySZ2QJJP4odCUCj4DQlgIZmZhsktY+GmCS5+/OS2s1sTtxkABoFR8YA1CV3f1rSmir2u1fSveM83zUNYwHAOTgyBqARFSV1nb3o63hGF32V1CSpNM1zAWhAVr5xCAAAABE4MgYAABCIGAMAAAhEjAEAAAQixgAAAAIRYwAAAIGIMQAAgEDEGAAAQKD/D+Ek1yBJ3W++AAAAAElFTkSuQmCC\n",
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
"text/plain": [
"<Figure size 720x504 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"morseB.plot()"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "b93e8548-cdcf-4d51-bbc0-ae7e5646cd02",
"metadata": {},
"outputs": [],
"source": [
"rho = job.factories.functions.spline(identifier=\"rho_AlAl\", cutoff=7.6, derivative_left=-0.2, derivative_right=0)\n",
"## the spline requires node points and initial values\n",
"init_func = lambda r: np.exp(-r)\n",
"nodes = np.array([1.0, 2.2, 2.7, 3.2, 4.0, 5.0])\n",
"init_vals = init_func(nodes)\n",
"rho.parameters.create_from_arrays(x=nodes, y=init_vals, min_vals=np.zeros(len(nodes)), max_vals=np.ones(len(nodes)))\n",
"# set node point at cutoff\n",
"rho.parameters.add_node(x=7.6, start_val=0, enabled=False)\n",
"rho.derivative_left.max_val = 0.0"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "56710703-5e6e-485e-a9d9-bc2c4bedac5b",
"metadata": {},
"outputs": [],
"source": [
"# User function for embedding term\n",
"F = job.factories.functions.user_function(identifier=\"F\", input_variable=\"r\")\n",
"F.expression = \"-A*sqrt(r)+B*r\"\n",
"F.derivative = \"-A/(2*sqrt(r))+B\"\n",
"F.parameters.add_parameter(\"A\", start_val=3.3, min_val=0.0)\n",
"F.parameters.add_parameter(\"B\", start_val=1.8, min_val=0.0)"
]
},
{
"cell_type": "markdown",
"id": "865bca67-41c4-4bb1-b357-ba2c4d37ce8b",
"metadata": {},
"source": [
"Assign functions"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "ce82c058-99b1-4044-89e4-d5b9fc6191d5",
"metadata": {},
"outputs": [],
"source": [
"job.potential.pair_interactions[morseB.identifier] = morseB\n",
"job.potential.electron_densities[rho.identifier] = rho\n",
"job.potential.embedding_energies[F.identifier] = F"
]
},
{
"cell_type": "markdown",
"id": "3fa3875d-c6ac-4a93-8917-0219be6bd90e",
"metadata": {},
"source": [
"Set a fit algorithm"
]
},
{
"cell_type": "code",
"id": "d0c70020-3999-482a-9fc3-037d6b4e5be9",
"metadata": {},
"outputs": [],
"source": [
"job.input.fit_algorithm = job.factories.algorithms.ar_lbfgs(max_iter=10000)"
]
},
{
"cell_type": "code",
"id": "6e161ce8-8d33-4eea-8d41-889287ccb7c5",
"metadata": {},
"source": [
"job.run()"
]
},
{
"cell_type": "markdown",
"id": "cd20a5b9-74f8-4609-96f8-d8c84f7a6976",
"metadata": {},
"source": [
"Have a look at some of the outputs"
]
},
{
"cell_type": "code",
"id": "8a895b83-d90e-489f-8382-393a520ae0a3",
"metadata": {},
"source": [
"plt.plot(job.output.iterations, job.output.residual)\n",
"job.output.residual[-1]"
]
},
{
"cell_type": "code",
"id": "74465658-bc69-486c-a327-4cd7a5790ec7",
"metadata": {},
"source": [
"job.potential"
]
},
{
"cell_type": "code",
"id": "2fe49eb2-bccb-4ba6-9a9f-8af61f6715c5",
"metadata": {},
"source": [
"job.plot_final_potential()"
]
},
{
"cell_type": "markdown",
"id": "fd58d03d-79bb-4a27-a7bd-371e6dba48a4",
"metadata": {},
"source": [
"Acces the plotting interface common to the different fitting codes usable via pyiron"
]
},
{
"cell_type": "code",
"id": "1230ba4e-4a4b-42fe-9355-5e1360fab63d",
"metadata": {},
"outputs": [],
"source": [
"plots = job.plot"
]
},
{
"cell_type": "code",
"id": "7f94c7c7-eda9-4d39-8e82-a2e498b069cb",
"metadata": {},
"source": [
"plt.figure(figsize=(14,7)) # Increase the size a bit\n",
"plots.energy_scatter_histogram()"
]
},
{
"cell_type": "code",
"id": "92a3bc17-de4a-459d-8ccf-4b0e59283f80",
"metadata": {},
"source": [
"plt.figure(figsize=(14,7))\n",
"plots.force_scatter_histogram()"
]
},
{
"cell_type": "code",
"id": "6bc97148-4858-4493-b527-e57037a9f7cb",
"metadata": {},
"source": [
"plt.figure(figsize=(14,7))\n",
"plots.force_scatter_histogram(axis=0)"
]
},
{
"cell_type": "markdown",
"id": "42b4825a-9f46-4cf2-9ab6-3ab2da665eba",
"metadata": {
"tags": []
},
"source": [
"### Short Test"
]
},
{
"cell_type": "code",
"id": "e233952f-ff1d-42d9-bb61-ea94c4461368",
"metadata": {},
"source": [
"lmp = pr.create.job.Lammps(\"Al\", delete_existing_job=True)\n",
"lmp.structure = pr.create.structure.ase.bulk(\"Al\", cubic=True)\n",
"lmp.potential = job.lammps_potential\n",
"lmp.calc_minimize(pressure=0)\n",
"lmp.run()"
]
},
{
"cell_type": "code",
"id": "822bae18-b13c-4de6-949b-beb34ffe60b8",
"metadata": {},
"source": [
"lmp.get_structure()"
]
},
{
"cell_type": "code",
"id": "a1eb8451-356b-4768-b903-fc4843db40d5",
"metadata": {},
"source": [
"ref = pr.create.job.Lammps(\"AlRef\")\n",
"ref.structure = lmp.get_structure()\n",
"ref.potential = job.lammps_potential\n",
"\n",
"murn = pr.create.job.Murnaghan(\"AlMurn\", delete_existing_job=True)\n",
"murn.ref_job = ref\n",
"murn.run()"
]
},
{
"cell_type": "code",
"id": "bce8a797-d906-4c0b-b674-eeb7f09169d3",
"metadata": {},
"source": [
"murn.plot()"
]
},
{
"cell_type": "code",
"id": "f82f8b1e-bc91-4c16-a508-784f6ba1a5be",
"metadata": {},
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
"source": [
"murn.fit_polynomial()"
]
},
{
"cell_type": "markdown",
"id": "af0bb904-c780-4aeb-a975-a4b800c10c2b",
"metadata": {},
"source": [
"### Improving the potential:\n",
"- try with different starting parameters / global optimization\n",
"- change weights of structures and fit properties\n",
"- change functions used"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1c654e1d-8c92-4b43-a78a-b0222b5f6128",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
}
},
"nbformat": 4,
"nbformat_minor": 5
}