{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "357d9be4",
   "metadata": {
    "tags": []
   },
   "source": [
    "# 原子間ポテンシャル\n",
    "\n",
    "　原子間ポテンシャルは、２つの原子を近づけていくと、ある距離までエネルギーが減少する。これは２つの原子に電子が共有されることが原因である。一方で近づきすぎると、電子同士の反発する力が強くなるためエネルギーは急上昇する。\n",
    " \n",
    "　このような原子間ポテンシャルは、第一原理計算と呼ばれる計算手法によって求めることが可能である。ここでは、[Materials Project](https://next-gen.materialsproject.org)に収録されている計算データを利用する。\n",
    " \n",
    "　Jupyter notebookを使用して、FCC構造のCuについて原子間ポテンシャル（全エネルギーと体積との関係）を図示してみましょう。\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c53d3a0e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from pylab import * #this includes numpy as np!\n",
    "from scipy.optimize import leastsq\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (5,4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fa6344a",
   "metadata": {},
   "source": [
    "　FCC構造のCuの計算データはmp-30として登録されている。体積をv、エネルギーをeとします。またvを原子間距離のrにするため、``np.cbrt(v)``としてvの立方根をとる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e8fc9386",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Energy (eV)')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdYAAAF3CAYAAAAYWmoRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAyzklEQVR4nO3de1xUdf7H8feICEhCrqPrEKMkqaT2SzfNdLVkvWYbKj82L6WZq4/tV+6mlReqVXy0RZaVZfe2X2kpriWm2UXaEjOtH9pKrVZYXgABLUQHJXcsnN8fs8xK3GbgzA1ez8djHjSH75n5zAl5c77ne75fk8PhcAgAABiilb8LAACgOSFYAQAwEMEKAICBCFYAAAxEsAIAYCCCFQAAAxGsAAAYiGAFAMBArf1dQKA7d+6ciouL1a5dO5lMJn+XAwDwE4fDoVOnTikmJkatWtV9XkqwNqC4uFhWq9XfZQAAAkRhYaFiY2Pr/H7QBGtSUpJyc3P13XffqX379hoxYoSWLl2qmJiYOvdJS0vT2rVrVVhYqDZt2uiKK67QAw88oIEDB7r9vu3atZPkPJBRUVFN/hwAgOBUXl4uq9XqyoW6mIJlruDHH39cgwYNksViUVFRke6++25J0s6dO+vcZ82aNerUqZO6deumM2fO6PHHH9frr7+ub7/9Vh07dnTrfcvLyxUdHS2bzUawAkAL5m4eBE2w/tymTZs0fvx42e12hYaGurVP1UH5+9//ruHDh3u0D8EKAC2bu3kQNF3B5ysrK9Pq1as1ePBgt0P17NmzeuGFFxQdHa3LL7+8znZ2u112u931vLy8vMn1AgBajqC63WbBggWKjIxUhw4dVFBQoI0bNza4z+bNm3XBBRcoPDxcjz/+uN5//32ZzeY626enpys6Otr1YOASAMATfg3WtLQ0mUymeh+7d+92tZ83b5727NmjrKwshYSEaNq0aWqoJzsxMVG5ubnauXOnxowZoxtuuEHfffddne1TU1Nls9lcj8LCQsM+LwCg+fPrNdbS0lKVlpbW2yYuLk7h4eE1th85ckRWq1U7d+7UoEGD3H7P7t27a8aMGUpNTXWrPddYAQBSkFxjNZvN9XbL1qfq74Hzr4e6u5+n+wAA4K6guMaak5Ojp556Srm5ucrPz9fWrVs1ZcoUxcfHVztbTUhI0IYNGyRJFRUVuueee/Tpp58qPz9f//jHPzRz5kwdOXJEv/vd73z7ASorpexsKSPD+bWy0rfvDwDwmaAYFRwREaHMzEwtXrxYFRUVslgsGjNmjNauXauwsDBXu7y8PNlsNklSSEiIvv76a61cuVKlpaXq0KGDBgwYoO3bt6t3796+Kz4zU7rjDunIkf9si42VnnhCSk72XR0AAJ8I2vtYfaVJ11gzM6WUFOnnh7hqzuE33iBcASBIuJsHQdEVHJQqK51nqrX93VK1bc4cuoUBoJkhWL1l+/bq3b8/53BIhYXOdgCAZoNg9ZaSEmPbAQCCAsHqLRaLse0AAEGBYPWWoUOdo3/rWhzdZJKsVmc7AECzQbB6S0iI85YaqWa4Vj1fvtzZDgDQbBCs3pSc7Lyl5qKLqm+PjeVWGwBopoJigoiglpwsjRvnHP1bUuK8pjp0KGeqANBMEay+EBIiDRvm7yoAAD5AVzAAAAYiWAEAMBDBCgCAgQhWAAAMRLACAGAgghUAAAMRrAAAGIhgBQDAQAQrAAAGIlgBADAQwQoAgIEIVgAADESwAgBgIIIVAAADEawAABiIYAUAwEBBE6xJSUnq0qWLwsPDZbFYNHXqVBUXF7u9/x/+8AeZTCYtX77ce0UCAFq8oAnWxMRErVu3Tnl5eVq/fr0OHDiglJQUt/Z988039X//93+KiYnxcpUAgJautb8LcNfcuXNd/921a1ctXLhQ48eP148//qjQ0NA69ysqKtLs2bO1ZcsWXXfddb4oFQDQggVNsJ6vrKxMq1ev1uDBg+sN1XPnzmnq1KmaN2+eevfu7dZr2+122e121/Py8vIm1wsAaDmCpitYkhYsWKDIyEh16NBBBQUF2rhxY73tly5dqtatW+tPf/qT2++Rnp6u6Oho18NqtTa1bABAC+LXYE1LS5PJZKr3sXv3blf7efPmac+ePcrKylJISIimTZsmh8NR62t/9tlneuKJJ/TKK6/IZDK5XVNqaqpsNpvrUVhY2OTPCQBoOUyOupLJB0pLS1VaWlpvm7i4OIWHh9fYfuTIEVmtVu3cuVODBg2q8f3ly5frzjvvVKtW//nbobKyUq1atZLVatXhw4fdqrG8vFzR0dGy2WyKiopyax8AQPPjbh749Rqr2WyW2Wxu1L5Vfw+cfz30fFOnTtWIESOqbRs9erSmTp2qW265pVHvCQBAQ4Ji8FJOTo5ycnI0ZMgQtW/fXgcPHtSiRYsUHx9f7Ww1ISFB6enpmjBhgjp06KAOHTpUe53Q0FB17txZPXv29PVHAAC0EEExeCkiIkKZmZkaPny4evbsqRkzZqhPnz7atm2bwsLCXO3y8vJks9n8WCkAoKXz6zXWYMA1VgCA5H4eBMUZKwAAwYJgBQDAQAQrAAAGIlgBADAQwQoAgIEIVgAADESwAgBgIIIVAAADEawAABiIYAUAwEAEKwAABiJYAQAwEMEKAICBCFYAAAxEsAIAYCCCFQAAAxGsAAAYiGAFAMBABCsAAAYiWAEAMBDBCgCAgQhWAAAM1NrfBcBNlZXS9u1SSYlksUhDh0ohIf6uCgDwMwRrMMjMlO64Qzpy5D/bYmOlJ56QkpP9VxcAoAa6ggNdZqaUklI9VCWpqMi5PTPTP3UBAGpFsAayykrnmarDUfN7VdvmzHG2AwAEhKAJ1qSkJHXp0kXh4eGyWCyaOnWqiouL691n+vTpMplM1R5XXXWVjyo2wPbtNc9Uz+dwSIWFznYAgIAQNMGamJiodevWKS8vT+vXr9eBAweUkpLS4H5jxoxRSUmJ6/HOO+/4oFqDlJQY2w4A4HVBM3hp7ty5rv/u2rWrFi5cqPHjx+vHH39UaGhonfuFhYWpc+fOvijReBaLse0AAF4XNGes5ysrK9Pq1as1ePDgekNVkrKzs9WpUyf16NFDs2bN0nfffVdve7vdrvLy8moPvxk61Dn612Sq/fsmk2S1OtsBAAJCUAXrggULFBkZqQ4dOqigoEAbN26st/21116r1atX68MPP9Sjjz6qXbt26Te/+Y3sdnud+6Snpys6Otr1sFqtRn8M94WEOG+pkWqGa9Xz5cu5nxUAAojJ4ahtyKlvpKWlacmSJfW22bVrl/r37y9JKi0tVVlZmfLz87VkyRJFR0dr8+bNMtV1RvczJSUl6tq1q9auXavkOu7/tNvt1YK3vLxcVqtVNptNUVFRbn4yg9V2H6vV6gxV7mMFAJ8oLy9XdHR0g3ng12AtLS1VaWlpvW3i4uIUHh5eY/uRI0dktVq1c+dODRo0yO337N69u2bOnKkFCxa41d7dA+l1zLwEAH7lbh74dfCS2WyW2Wxu1L5Vfw/U1637c8ePH1dhYaEswTjYJyREGjbM31UAABoQFNdYc3Jy9NRTTyk3N1f5+fnaunWrpkyZovj4+GpnqwkJCdqwYYMk6fTp07r77rv1ySef6PDhw8rOztb1118vs9msCRMm+OujAACauaC43SYiIkKZmZlavHixKioqZLFYNGbMGK1du1ZhYWGudnl5ebLZbJKkkJAQ/fOf/9SqVat08uRJWSwWJSYm6m9/+5vatWvnr48CAGjm/HqNNRgEzDVWAIBfuZsHQdEVDABAsCBYAQAwEMEKAICBCFYAAAxEsAIAYCCCFQAAAxGsAAAYiGAFAMBABCsAAAYiWAEAMBDBCgCAgQhWAAAMRLACAGAgghUAAAMRrAAAGIhgBQDAQAQrAAAGIlgBADAQwQoAgIEIVgAADESwAgBgIIIVAAADEawAABiIYAUAwEAEKwAABgqaYE1KSlKXLl0UHh4ui8WiqVOnqri4uMH9vvrqKyUlJSk6Olrt2rXTVVddpYKCAh9UDABoiYImWBMTE7Vu3Trl5eVp/fr1OnDggFJSUurd58CBAxoyZIgSEhKUnZ2tzz//XH/+858VHh7uo6oBAC2NyeFwOPxdRGNs2rRJ48ePl91uV2hoaK1tJk2apNDQUL366quNfp/y8nJFR0fLZrMpKiqq0a8DAAhu7uZB0Jyxnq+srEyrV6/W4MGD6wzVc+fO6e2331aPHj00evRoderUSQMHDtSbb75Z72vb7XaVl5dXewAA4K6gCtYFCxYoMjJSHTp0UEFBgTZu3Fhn2++++06nT5/WQw89pDFjxigrK0sTJkxQcnKytm3bVud+6enpio6Odj2sVqs3PgoAoJnyqCvY4XBo27Zt2r59uw4fPqwffvhBHTt2VL9+/TRixAiPQygtLU1Lliypt82uXbvUv39/SVJpaanKysqUn5+vJUuWKDo6Wps3b5bJZKqxX3FxsS666CJNnjxZa9ascW1PSkpSZGSkMjIyan0/u90uu93uel5eXi6r1UpXMAC0cO52Bbd258XOnDmjxx9/XM8884yOHz+uyy+/XBdddJEiIiL07bff6s0339SsWbM0atQoLVq0SFdddZVbRc6ePVuTJk2qt01cXJzrv81ms8xms3r06KFLL71UVqtVn376qQYNGlRjP7PZrNatW6tXr17Vtl966aX6+OOP63y/sLAwhYWFuVU/AAA/51aw9ujRQwMHDtRzzz2n0aNH13pdMz8/X2vWrNHEiRN13333adasWQ2+blVQNkbVifb5Z5fna9OmjQYMGKC8vLxq2/fv36+uXbs26j0BAGiIW13Be/fuVZ8+fdx6wbNnzyo/P1/du3dvcnFVcnJylJOToyFDhqh9+/Y6ePCgFi1apJKSEu3bt891hpmQkKD09HRNmDBBkrRhwwZNnDhRTz/9tBITE/Xee+9pzpw5ys7O1pAhQ9x6b0YFAwAkg0cF9+nTR7m5uW69cZs2bQwNVUmKiIhQZmamhg8frp49e2rGjBnq06ePtm3bVq3bNi8vTzabzfV8woQJeu655/Twww/rsssu01//+letX7/e7VAFAMBTbg9eatWqlfr166eZM2dqypQpio6O9nZtAYEzVgCA5IX7WHfs2KFf/epXWrhwoSwWi2666SZt3brVkGIBAGgu3A7WQYMG6cUXX9TRo0f17LPP6siRIxoxYoTi4+P1wAMP6MiRI96sE0aprJSys6WMDOfXykp/VwQAzYrHE0RERETo5ptvVnZ2tvbv36/Jkyfr+eef18UXX6yxY8d6o0YYJTNTiouTEhOlKVOcX+PinNsBAIZo8lzBp0+f1urVq3XPPffo5MmTqmxmZ0DN5hprZqaUkiL9/H931eQab7whJSf7vi4ACBJenyt427Ztuvnmm9W5c2fNnz9fycnJ2rFjR2NfDt5UWSndcUfNUJX+s23OHLqFAcAAbk0QUaWwsFCvvPKKXnnlFR06dEiDBw/WihUrdMMNNygyMtJbNaKptm+X6rsG7nBIhYXOdsOG+awsAGiO3A7WkSNHauvWrerYsaOmTZumGTNmqGfPnt6sDUYpKTG2HQCgTm4Ha0REhNavX6/f/va3CgkJ8WZNMJrFYmw7AECdGj146dtvv9WBAwd09dVXKyIiQg6Ho9ZVZoJdsxi8VFnpHP1bVFT7dVaTSYqNlQ4dkvijCQBq5bXBS8ePH9fw4cPVo0cPjR07ViX/7j6cOXOm7rrrrsZXDO8JCZGeeML53z//46fq+fLlhCoAGMDjYJ07d65CQ0NVUFCgtm3burZPnDhR7733nqHFwUDJyc5bai66qPr22FhutQEAA3k0KliSsrKytGXLFsXGxlbb3r17d+Xn5xtWGLwgOVkaN845+rekxHlNdehQzlQBwEAeB2tFRUW1M9UqpaWlLBAeDEJCuKUGALzI467gq6++WqtWrXI9N5lMOnfunB555BElJiYaWhwAAMHG4zPWRx55RMOGDdPu3bt19uxZzZ8/X/v27VNZWRkzLwEAWjyPz1h79eqlL774QldeeaVGjhypiooKJScna8+ePYqPj/dGjQAABI0mT8Lf3DWL+1gBAE1m6H2sBQUFHr15UVGRR+0BAGgu3ArWAQMGaNasWcrJyamzjc1m04svvqg+ffook/U9AQAtlFuDl7766is9+OCDGjNmjEJDQ9W/f3/FxMQoPDxcJ06c0Jdffql9+/apf//+euSRR3Tttdd6u24AAAKSR9dY//Wvf+mdd97R9u3bdfjwYZ05c0Zms1n9+vXT6NGj1adPH2/W6hdcYwUASO7nAYOXGkCwAgAkL07CDwAA6kawAgBgIIIVAAADBU2wJiUlqUuXLgoPD5fFYtHUqVNVXFxc7z4mk6nWxyOPPOKjqgEALY3HwVpRUeGNOhqUmJiodevWKS8vT+vXr9eBAweUkpJS7z4lJSXVHv/7v/8rk8mk//7v//ZR1QCAlsbjUcEXXHCBbrjhBs2YMUNDhgzxVl0N2rRpk8aPHy+73a7Q0FC39hk/frxOnTqlDz74wO33YVQwAEDy4qjgjIwM2Ww2DR8+XD169NBDDz3UYJes0crKyrR69WoNHjzY7VA9duyY3n77bf3+97+vt53dbld5eXm1BwAA7vI4WK+//nqtX79excXF+p//+R9lZGSoa9eu+u1vf6vMzEz99NNP3qhTkrRgwQJFRkaqQ4cOKigo0MaNG93ed+XKlWrXrp2Sk5PrbZeenq7o6GjXw2q1NrVsAEALYsgEEStWrNC8efN09uxZmc1m3XrrrVq4cKHatm1b735paWlasmRJvW127dql/v37S5JKS0tVVlam/Px8LVmyRNHR0dq8ebNMJlODNSYkJGjkyJFasWJFve3sdrvsdrvreXl5uaxWK13BANDCeX3mpaNHj2rVqlV6+eWXVVBQoAkTJuj3v/+9iouL9dBDD8lisSgrK6ve1ygtLVVpaWm9beLi4hQeHl5j+5EjR2S1WrVz504NGjSo3tfYvn27rr76auXm5uryyy9v+MOdh2usAADJ/TxwaxL+82VmZurll1/Wli1b1KtXL91+++266aabdOGFF7ra9O3bV/369Wvwtcxms8xms6clSJKq/h44/+yyLi+99JKuuOIKj0MVAABPeXyN9ZZbblFMTIx27Nih3NxczZ49u1qoSlK3bt107733GlWjcnJy9NRTTyk3N1f5+fnaunWrpkyZovj4+GpnqwkJCdqwYUO1fcvLy/X6669r5syZhtUDAEBdPD5jLSkpafDaaUREhBYvXtzoomp7vczMTC1evFgVFRWyWCwaM2aM1q5dq7CwMFe7vLw82Wy2avuuXbtWDodDkydPNqweAADq4vE11rpuPzGZTAoLC1ObNm0MKSxQcI0VACB58RrrhRdeWO8o3NjYWE2fPl2LFy9Wq1ZBM2MiAACG8DhYX3nlFd17772aPn26rrzySjkcDu3atUsrV67Ufffdp++//17Lli1TWFiY7rnnHm/UDABAwPI4WFeuXKlHH31UN9xwg2tbUlKSLrvsMj3//PP64IMP1KVLFz3wwAMEKwCgxfG4r/aTTz6p9Vaafv366ZNPPpEkDRkyRAUFBU2vDgCAIONxsMbGxuqll16qsf2ll15yTf93/PhxtW/fvunVAQAQZDzuCl62bJl+97vf6d1339WAAQNkMpm0a9cuff3113rjjTckOachnDhxouHFAgAQ6Bo1pWF+fr6ee+455eXlyeFwKCEhQX/4wx8UFxfnhRL9i9ttAACSl263+fHHHzVq1Cg9//zzSk9Pb3KRAAA0Nx5dYw0NDdXevXvdWk0GAICWyOPBS9OmTat18BJakMpKKTtbyshwfq2s9HdFABAwPB68dPbsWf31r3/V+++/r/79+ysyMrLa9x977DHDikMAysyU7rhDOnLkP9tiY6UnnpAaWEQeAFoCj4N17969+tWvfiVJ2r9/f7Xv0UXczGVmSikp0s/HuxUVObe/8QbhCqDFa/RC5y0Fo4L/rbJSiourfqZ6PpPJeeZ66JAUEuLT0gDAF9zNg0bPkv/tt99qy5YtOnPmjKT/LDyOZmr79rpDVXKexRYWOtsBQAvmcbAeP35cw4cPV48ePTR27FiVlJRIkmbOnKm77rrL8AIRIP79/9mwdgDQTHkcrHPnzlVoaKgKCgqqLXg+ceJEvffee4YWhwBisRjbDgCaKY8HL2VlZWnLli2KjY2ttr179+7Kz883rDAEmKFDnddQi4pqDl6S/nONdehQ39cGAAHE4zPWioqKameqVUpLSxUWFmZIUQhAISHOW2okZ4ier+r58uUMXALQ4nkcrFdffbVWrVrlem4ymXTu3Dk98sgjSkxMNLQ4BJjkZOctNRddVH17bCy32gDAv3l8u82XX36pYcOG6YorrtCHH36opKQk7du3T2VlZdqxY4fi4+O9VatfcLtNLSornaN/S0qc11SHDuVMFUCz55VJ+CWpV69e+uKLL/Tss88qJCREFRUVSk5O1u233y4LA1dahpAQadgwf1cBAAGJCSIawBkrAEDy4hmrJJ08eVI5OTn67rvvdO7cuWrfmzZtWmNeEgCAZsHjYH3rrbd04403qqKiQu3atas2P7DJZCJYAQCBww9jQjweFXzXXXdpxowZOnXqlE6ePKkTJ064HmVlZd6oEQAAz2VmOuc4T0yUpkxxfo2Lc273Io+DtaioSH/6059qvZcVAICAULUa18/nOK9ajcuL4epxsI4ePVq7d+/2Ri31SkpKUpcuXRQeHi6LxaKpU6equLi43n1Onz6t2bNnKzY2VhEREbr00kv17LPP+qhiAIBfVFY6142ubWxu1bY5c5ztvMDja6zXXXed5s2bpy+//FKXXXaZQkNDq30/KSnJsOLOl5iYqHvuuUcWi0VFRUW6++67lZKSop07d9a5z9y5c7V161a99tpriouLU1ZWlm677TbFxMRo3LhxXqkTAOBnnqzG5YVbBz2+3aZVq7pPck0mkyq99BfAz23atEnjx4+X3W6vEe5V+vTpo4kTJ+rPf/6za9sVV1yhsWPH6v777691H7vdLrvd7npeXl4uq9XK7TYAECwyMpzXVBuyZo00ebLbL+u19VjPnTtX58NXoVpWVqbVq1dr8ODBdYaqJA0ZMkSbNm1SUVGRHA6Htm7dqv3792v06NF17pOenq7o6GjXw2q1euMjAAC8xc+rcTV6oXN/WLBggSIjI9WhQwcVFBRo48aN9bZ/8skn1atXL8XGxqpNmzYaM2aMnnnmGQ0ZMqTOfVJTU2Wz2VyPwsJCoz8GAMCbqlbj+vmCIVVMJslq9dpqXG4H69ixY2Wz2VzPH3jgAZ08edL1/Pjx4+rVq5dHb56WliaTyVTv4/yBUvPmzdOePXuUlZWlkJAQTZs2TfX1ZD/55JP69NNPtWnTJn322Wd69NFHddttt+nvf/97nfuEhYUpKiqq2gMAEET8vBqX29dYQ0JCVFJSok6dOkmSoqKilJubq27dukmSjh07ppiYGI+6g0tLS1VaWlpvm7i4OIWHh9fYfuTIEVmtVu3cuVODBg2q8f0zZ84oOjpaGzZs0HXXXefaPnPmTB05csTtRdmZ0hAAglRmpnN08PkDmaxWZ6g2YjUuw6c0/Hn+GjHFsNlsltlsbtS+Ve9//kCj8/3444/68ccfawy2CgkJqTENIwCgGUpOlsaN8/nMS42aK9jXcnJylJOToyFDhqh9+/Y6ePCgFi1apPj4+GpnqwkJCUpPT9eECRMUFRWla665RvPmzVNERIS6du2qbdu2adWqVXrsscf8+GkAAD7jh9W43A7WqmueP9/mCxEREcrMzNTixYtVUVEhi8WiMWPGaO3atQoLC3O1y8vLq3YdeO3atUpNTdWNN96osrIyde3aVQ888IBuvfVWn9QNAGh53L7G2qpVK1177bWuIHvrrbf0m9/8RpGRkZKcXbLvvfeez2658RWusQIAJC9cY7355purPb/ppptqtGFlGwBAS+d2sL788sverAMAgGYhqCaIAAAg0BGsAAAYKChutwEAtDCVlT6//9QoBCsAILDUNmNSbKxzmsJGzJjka3QFAwACR2amlJJScz3VoiLn9sxM/9TlAYIV/lNZKWVnO9dOzM52PgfQclVWOs9Ua5teoWrbnDkB/7uCYIV/ZGZKcXFSYqJzQeLEROfzIPhrFICXbN9e80z1fA6HVFjobBfACFb4XjPo6gHgBSUlxrbzE4IVvtVMunoAeIHFYmw7PyFY4VvNpKsHgBcMHeoc/VvXAi8mk3M91aFDfVuXhwhW+FYz6eoB4AUhIc5baqSa4Vr1fPnygL+flWCFbzWTrh4AXpKcLL3xhnTRRdW3x8Y6twfBfaxuLxvXUrFsnMEqK52jf4uKar/OajI5/wEdOhTwf5UC8KIAnHnJ8GXjAENUdfWkpDhD9PxwDaKuHgBeFhIiDRvm7yoaha5g+F4z6OoBgLpwxgr/SE6Wxo0LuK4eAGgqghX+E8RdPQBQF7qCAQAwEMEKAICB6AoGABgrAG+V8SWCFQBgnCBfpNwIdAUDAIzBylWSCFYAgBFYucolaII1KSlJXbp0UXh4uCwWi6ZOnari4uJ69zl27JimT5+umJgYtW3bVmPGjNE333zjo4oBoAVh5SqXoAnWxMRErVu3Tnl5eVq/fr0OHDiglJSUOts7HA6NHz9eBw8e1MaNG7Vnzx517dpVI0aMUEVFhQ8rB4AWgJWrXIJ2Ev5NmzZp/PjxstvtCg0NrfH9/fv3q2fPntq7d6969+4tSaqsrFSnTp20dOlSzZw50633YRJ+AHBDdraUmNhwu61bg3ZiGHfzIGjOWM9XVlam1atXa/DgwbWGqiTZ7XZJUnh4uGtbSEiI2rRpo48//rjO17bb7SovL6/2AAA0oJksUm6EoArWBQsWKDIyUh06dFBBQYE2btxYZ9uEhAR17dpVqampOnHihM6ePauHHnpIR48eVUk9XRHp6emKjo52PaxWqzc+CgA0L81kkXIj+DVY09LSZDKZ6n3s3r3b1X7evHnas2ePsrKyFBISomnTpqmunuzQ0FCtX79e+/fv1y9+8Qu1bdtW2dnZuvbaaxVSz//Y1NRU2Ww216OwsNDwzw0AzRIrV0ny8zXW0tJSlZaW1tsmLi6uWndulSNHjshqtWrnzp0aNGhQva9hs9l09uxZdezYUQMHDlT//v319NNPu1Uj11gBwEPNdOaloFjo3Gw2y2w2N2rfqr8Hqq6l1ic6OlqS9M0332j37t26//77G/WeAAA3tPCVq4LiGmtOTo6eeuop5ebmKj8/X1u3btWUKVMUHx9f7Ww1ISFBGzZscD1//fXXlZ2d7brlZuTIkRo/frxGjRrlj48BAGgBgmKu4IiICGVmZmrx4sWqqKiQxWLRmDFjtHbtWoWFhbna5eXlyWazuZ6XlJTozjvv1LFjx2SxWDRt2jT9+c9/9sdHAAC0EEF7H6uvcI01CDTT6zkAAktQXGMFmoyVNAAEmKC4xgrUipU0AGNVVjpnUMrIcH5tARPmewPBiuDEShqAsTIzpbg457SEU6Y4v8bF8QdqIxCsCE6spAEYh94fQxGsCE6spAEYg94fwxGsCE4Wi7HtgJaK3h/DEawITqykARiD3h/DEawITqykARiD3h/DEawIXqykATQdvT+GY4IIBLfkZGncOGZeAhqrqvcnJcUZoucPYqL3p1EIVgS/Fr6SBtBkVb0/tc1itnw5vT8eIlgBAPT+GIhgBQA40ftjCAYvAQBgIIIVAAAD0RUMAMGONYkDCsEKAMGMNYkDDl3BABCsWJUmIBGsABCMWJUmYBGsABCMWJUmYBGsABCMWJUmYDF4CZAYVYngw6o0AYszViAzU4qLkxITpSlTnF/j4hj4gcDGqjQBi2BFy8aoSgQr1iQOWAQrWi5GVSLYsSZxQAq6YLXb7erbt69MJpNyc3PrbetwOJSWlqaYmBhFRERo2LBh2rdvn28KReBjVCWag+Rk6fBhaetWac0a59dDhwhVPwq6YJ0/f75iYmLcavvwww/rscce01NPPaVdu3apc+fOGjlypE6dOuXlKhEUGFWJ5qJqVZrJk51f6f71q6AK1nfffVdZWVlatmxZg20dDoeWL1+ue++9V8nJyerTp49WrlypH374QWvWrPFBtQh4jKoE4AVBE6zHjh3TrFmz9Oqrr6pt27YNtj906JCOHj2qUaNGubaFhYXpmmuu0c6dO+vcz263q7y8vNoDzRSjKhEIKiul7GwpI8P5lWv6QS8ogtXhcGj69Om69dZb1b9/f7f2OXr0qCTpl7/8ZbXtv/zlL13fq016erqio6NdD6vV2vjCEdgYVQl/41avZsmvwZqWliaTyVTvY/fu3VqxYoXKy8uVmprq8XuYfvYL0+Fw1Nh2vtTUVNlsNtejsLDQ4/dEEGFUJfyFW72aLZPDUdu9Br5RWlqq0tLSetvExcVp0qRJeuutt6oFYmVlpUJCQnTjjTdq5cqVNfY7ePCg4uPj9Y9//EP9+vVzbR83bpwuvPDCWvepTXl5uaKjo2Wz2RQVFeXmJ0PQYeYl+FJlpfPMtK5R6SaT84+7Q4f4OQwg7uaBX6c0NJvNMpvNDbZ78skn9Ze//MX1vLi4WKNHj9bf/vY3DRw4sNZ9Lr74YnXu3Fnvv/++K1jPnj2rbdu2aenSpcZ8ADQfVaMqAV/w5FYvfi6DTlDMFdylS5dqzy+44AJJUnx8vGJjY13bExISlJ6ergkTJshkMmnOnDl68MEH1b17d3Xv3l0PPvig2rZtqylTpvi0fgCohlu9mrWgCFZ35eXlyWazuZ7Pnz9fZ86c0W233aYTJ05o4MCBysrKUrt27fxYJYAWj1u9mjW/XmMNBlxjBWC4qmusRUW1T6nJNdaA5G4eBMXtNgDQrHCrV7NGsAJG4mZ/uItbvZqtZnWNFfCrzEznajnnj/aMjXWemfBLErVJTpbGjeNWr2aGa6wN4Bor3FJ1s//P/zlVdetxBtL8cO9zi8M1VsBXWNe15WEqQtSDYAWainVdWxamIkQDCFagqbjZv+WgdwJuIFiBpuJm/5aD3gm4gWAFmop1XVsOeifgBoIVaCpu9m856J2AGwhWwAjc7N8y0DsBNzBBBGAUbvZv/qp6J1JSnCF6/iAmeifwbwQrYCTWdQ0OTZncoap3orZZtpYvp3cCBCsQkJjVx3uMmHqS3gnUgykNG8CUhvA55hz2HqaeRBMwpSEQjJjVx3uY3AE+QrACgYJf/N7F5A7wEYIVCBT84vcuJneAjzB4CQgU/OJvWFMGdTG5A3yEM1YgUPCLv35NXaqNyR3gIwQrECiM/sVfWSllZ0sZGc6vwXxt1ohBXUw9CR8hWIFAYeQv/ua0ELeRg7qYehI+wH2sDeA+VvhcbfexWq3uz+oTiPdqNuXaaHa28w+Dhmzd6v6sV0zAgUZwNw8YvAQEmqbM6tPQ2Z3J5Dy7GzfO/SBpagg1dcILbwzqYupJeBHBCgSixv7i9+SWHXdev6mhWNfZc9W1UXfOnhnUhSATdNdY7Xa7+vbtK5PJpNzc3HrbZmZmavTo0TKbzW61B4KekWd3TR0wZNS1UUbzIsgEXbDOnz9fMTExbrWtqKjQr3/9az300ENergoIEEad3RkRikZNeMFoXgSZoArWd999V1lZWVq2bJlb7adOnapFixZpxIgRXq4MCBBGnd0ZEYpGnj0zmhdBJGiusR47dkyzZs3Sm2++qbZt23rtfex2u+x2u+t5eXm5194LMJxRC3EbEYpGXxtlqTYEiaA4Y3U4HJo+fbpuvfVW9e/f36vvlZ6erujoaNfDarV69f0AwxlxdmdEKHrj2mjVoK7Jk51fCVUEIL8Ga1pamkwmU72P3bt3a8WKFSovL1dqaqrXa0pNTZXNZnM9CgsLvf6egOGSk6XDh533dq5Z4/x66JD7XaZGhCLXRtFC+XWCiNLSUpWWltbbJi4uTpMmTdJbb70l03n/OCsrKxUSEqIbb7xRK1eurPc1Dh8+rIsvvlh79uxR3759PaqRCSLQYlWNCpZq71J29+y3qRNeAAHC3TwIipmXCgoKql3rLC4u1ujRo/XGG29o4MCBio2NrXd/ghVoJKNCkZmO0Aw0q5mXunTpUu35BRdcIEmKj4+vFqoJCQlKT0/XhAkTJEllZWUqKChQcXGxJCkvL0+S1LlzZ3Xu3NkXpQPBzagBQ8x0hBYkKILVXXl5ebLZbK7nmzZt0i233OJ6PmnSJEnS4sWLlZaW5uvygOBEKAIeCYquYH+iKxgAILmfB0Fxuw0AAMGCYAUAwEAEKwAABiJYAQAwEMEKAICBmtXtNt5QNWiayfgBoGWryoGGbqYhWBtw6tQpSWIyfgCAJGcuREdH1/l97mNtwLlz51RcXKx27dpVm6s4mJWXl8tqtaqwsJB7c93EMfMcx6xxOG6e89UxczgcOnXqlGJiYtSqVd1XUjljbUCrVq0anIs4WEVFRfEP10McM89xzBqH4+Y5Xxyz+s5UqzB4CQAAAxGsAAAYiGBtgcLCwrR48WKFhYX5u5SgwTHzHMescThungu0Y8bgJQAADMQZKwAABiJYAQAwEMEKAICBCFYAAAxEsDYz6enpGjBggNq1a6dOnTpp/PjxysvLq3efkpISTZkyRT179lSrVq00Z84c3xQbIBpzzDIzMzVy5Eh17NhRUVFRGjRokLZs2eKjiv2vMcfs448/1q9//Wt16NBBERERSkhI0OOPP+6jigNDY47b+Xbs2KHWrVurb9++3isywDTmmGVnZ8tkMtV4fP311z6pmWBtZrZt26bbb79dn376qd5//3399NNPGjVqlCoqKurcx263q2PHjrr33nt1+eWX+7DawNCYY/bRRx9p5MiReuedd/TZZ58pMTFR119/vfbs2ePDyv2nMccsMjJSs2fP1kcffaSvvvpK9913n+677z698MILPqzcvxpz3KrYbDZNmzZNw4cP90GlgaMpxywvL08lJSWuR/fu3X1QMbfbNHvff/+9OnXqpG3btunqq69usP2wYcPUt29fLV++3PvFBShPj1mV3r17a+LEiVq0aJEXqwtMjT1mycnJioyM1KuvvurF6gKXJ8dt0qRJ6t69u0JCQvTmm28qNzfXN0UGGHeOWXZ2thITE3XixAldeOGFvi1QnLE2ezabTZL0i1/8ws+VBI/GHLNz587p1KlTLfY4N+aY7dmzRzt37tQ111zjrbICnrvH7eWXX9aBAwe0ePFiX5QV0Dz5WevXr58sFouGDx+urVu3ers0Fybhb8YcDofuvPNODRkyRH369PF3OUGhscfs0UcfVUVFhW644QYvVheYPD1msbGx+v777/XTTz8pLS1NM2fO9EGVgcfd4/bNN99o4cKF2r59u1q3btm/st09ZhaLRS+88IKuuOIK2e12vfrqqxo+fLiys7M96lFprJb9f6mZmz17tr744gt9/PHH/i4laDTmmGVkZCgtLU0bN25Up06dvFhdYPL0mG3fvl2nT5/Wp59+qoULF+qSSy7R5MmTvVxl4HHnuFVWVmrKlClasmSJevTo4cPqApO7P2s9e/ZUz549Xc8HDRqkwsJCLVu2zCfBKgeapdmzZztiY2MdBw8e9Gi/a665xnHHHXd4p6gA15hjtnbtWkdERIRj8+bNXqwscDX256zK/fff7+jRo4fBVQU+d4/biRMnHJIcISEhrofJZHJt++CDD3xUsf819WftL3/5iyMhIcHgqmrHGWsz43A49Mc//lEbNmxQdna2Lr74Yn+XFPAae8wyMjI0Y8YMZWRk6LrrrvNylYHFqJ8zh8Mhu91ucHWBy9PjFhUVpX/+85/Vtj3zzDP68MMP9cYbb7SIf99G/azt2bNHFovF4OpqR7A2M7fffrvWrFmjjRs3ql27djp69Kgk5+K8ERERkqTU1FQVFRVp1apVrv2qRhiePn1a33//vXJzc9WmTRv16tXL55/B1xpzzDIyMjRt2jQ98cQTuuqqq1z7REREuLUQcrBrzDF7+umn1aVLFyUkJEhy3te6bNky/fGPf/TPh/ADT49bq1atalxL7NSpk8LDw1vMuInG/KwtX75ccXFx6t27t86ePavXXntN69ev1/r1631TtE/Oi+Ezkmp9vPzyy642N998s+Oaa65pcL+uXbv6tHZ/acwxu+aaa2rd5+abb/Z5/f7QmGP25JNPOnr37u1o27atIyoqytGvXz/HM88846isrPT9B/CTxv77PN/ixYsdl19+uddrDRSNOWZLly51xMfHO8LDwx3t27d3DBkyxPH222/7rGbuYwUAwEDcxwoAgIEIVgAADESwAgBgIIIVAAADEawAABiIYAUAwEAEKwAABiJYAQAwEMEKAICBCFYAAAxEsAJolM2bN6tbt24aMGCA9u/f7+9ygIDBXMEAGqVHjx565plntG/fPn3yySdau3atv0sCAgJnrADqdPz4cXXq1EmHDx+u8T2z2axLLrlE3bp1q7FUXkpKih577DEfVQkEFs5YAdTp7rvv1okTJ/TSSy/V+N6LL76oW2+9VR07dtS+ffvUoUMH1/e++OILJSYm6tChQ4qKivJlyYDfccYKtHBnz56tdfuZM2f00ksvaebMmTW+99NPP+mJJ57Q/PnzderUKbVv377a9//rv/5LcXFxWr16tVdqBgIZwQq0MMOGDdPs2bN15513ymw2a+TIkbW2e/fdd9W6dWsNGjSoxveee+45devWTbfffrt++OEHffPNNzXaJCUlKSMjw/D6gUBHsAIt0MqVK9W6dWvt2LFDzz//fK1tPvroI/Xv37/G9hMnTuj+++/X0qVLFRsbq+joaOXm5tZod+WVVyonJ0d2u93o8oGARrACLdAll1yihx9+WD179lRCQkKtbQ4fPqyYmJga2xctWqQJEybo0ksvlST16tVLn3/+eY12F110kex2u44ePWps8UCAa+3vAgD4Xm1noj935swZhYeHV9v25Zdf6rXXXtNXX33l2tanT59az1gjIiIkST/88EPTigWCDMEKtECRkZENtjGbzTpx4kS1bXPnztXJkycVGxvr2nbu3DlZLJYa+5eVlUmSOnbs2MRqgeBCsAKoVb9+/fTaa6+5nm/evFmfffaZ9uzZo9at//OrY9euXZoxY4a+//77aiG6d+9excbGymw2+7RuwN8IVgC1Gj16tFJTU3XixAldcMEFuuuuuzRv3jz17du3Wruq+1Q///xzjRgxwrV9+/btGjVqlC9LBgICwQqgVpdddpn69++vdevWqaKiQidPntTs2bNrtLNarWrbtq1yc3Ndwfqvf/1LGzZs0JYtW3xdNuB3zLwEoE7vvPOO7r77bu3du1etWrl/E8HTTz+tjRs3Kisry4vVAYGJM1YAdRo7dqy++eYbFRUVyWq1ur1faGioVqxY4cXKgMDFGSsAAAZigggAAAxEsAIAYCCCFQAAAxGsAAAYiGAFAMBABCsAAAYiWAEAMBDBCgCAgQhWAAAM9P+yyexwkads2wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 500x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Cu(mp-30)\n",
    "v = np.array([8.6664, 8.9586, 9.5623, 9.8741, 10.5179, 9.2571, 10.8500, 11.5351, 12.9905, 12.2484, 13.7620, 12.6158, 14.5635, 11.1890, 14.9757, 15.8231, 10.1926, 15.3955, 11.8882, 14.1590, 13.3726])\n",
    "r = np.cbrt(v)\n",
    "e = np.array([-3.3243, -3.4941, -3.7562, -3.8536, -3.9921, -3.6372, -4.0371, -4.0873, -4.0660, -4.0937, -4.0115, -4.0836, -3.9364, -4.0681, -3.8927, -3.7961, -3.9314, -3.8457, -4.0956, -3.9763, -4.0416])\n",
    "\n",
    "\n",
    "plot(r,e,'ro')\n",
    "xlabel('r ($\\AA$)')\n",
    "ylabel('Energy (eV)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7bdf1cb",
   "metadata": {
    "tags": []
   },
   "source": [
    " 原子間ポテンシャルは$r_0 \\sim 2.3$ の周りで極小値を取ることがわかる。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edad656e",
   "metadata": {},
   "source": [
    "## 二次関数でフィッティング\n",
    "\n",
    "　この原子間ポテンシャルはどのような関数で表現されるのであろうか。原子間に働く相互作用を２つの球体を繋ぐバネのようなもの（調和振動子）で考えた場合、原子間ポテンシャルは原子間距離$r$の２乗に比例する。\n",
    " \n",
    "　そこで二次関数 $ar^2 + br + c$でフィッティングしてみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "eea51bed",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Energy (eV)')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdYAAAF3CAYAAAAYWmoRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABPYUlEQVR4nO3dd1hT9/4H8HcII4AMZQhIAEWEKtaFA4uDuttbB6V11Vm9XXbY1tX2V+3tbam1raO2tcPrqqNWcGu1A1xoFQW3qCgbB4IEGQGS/P6IUJCVQJKTwPv1POfRnHxP8uGI+eS7RSqVSgUiIiLSCTOhAyAiImpKmFiJiIh0iImViIhIh5hYiYiIdIiJlYiISIeYWImIiHSIiZWIiEiHmFiJiIh0yFzoAIydUqlEZmYm7OzsIBKJhA6HiIgEolKpkJ+fDw8PD5iZ1V4vZWKtR2ZmJqRSqdBhEBGRkUhLS4Onp2etz5tMYh05ciQSEhJw584dtGzZEoMHD8bixYvh4eFR6zWLFi3Cli1bkJaWBktLS/To0QOffPIJevfurfH72tnZAVDfSHt7+0b/HEREZJpkMhmkUmlFXqiNyFTWCl66dCmCg4Ph7u6OjIwMvPvuuwCA2NjYWq/ZtGkTXF1d0a5dOxQVFWHp0qX49ddfcf36dbi4uGj0vjKZDA4ODsjLy2NiJSJqxjTNByaTWB+1a9cujB49GnK5HBYWFhpdU35T/vjjDwwaNEira5hYiYiaN03zgck0BVeWk5ODjRs3om/fvhon1ZKSEvzwww9wcHBAly5dai0nl8shl8srHstkskbHS0REzYdJTbeZN28ebG1t4eTkhNTUVOzcubPea/bs2YMWLVpAIpFg6dKl+P333+Hs7Fxr+YiICDg4OFQcHLhERETaEDSxLlq0CCKRqM4jLi6uovycOXMQHx+PgwcPQiwWY/LkyaivJTs0NBQJCQmIjY3F8OHD8fzzz+POnTu1ll+wYAHy8vIqjrS0NJ39vERE1PQJ2seanZ2N7OzsOsv4+PhAIpFUO5+eng6pVIrY2FgEBwdr/J5+fn6YPn06FixYoFF59rESERFgIn2szs7OdTbL1qX8+0Dl/lBNr9P2GiIiIk2ZRB/ryZMnsXLlSiQkJCAlJQXR0dGYMGECfH19q9RWAwICsH37dgBAQUEB3nvvPZw4cQIpKSk4c+YMZsyYgfT0dDz33HOG/QEUCiAmBti8Wf2nQmHY9yciIoMxiVHB1tbWiIqKwsKFC1FQUAB3d3cMHz4cW7ZsgZWVVUW5xMRE5OXlAQDEYjGuXLmCdevWITs7G05OTujZsyeOHDmCTp06GS74qCjgzTeB9PR/znl6AsuXA2FhhouDiIgMwmTnsRpKo/pYo6KA8HDg0Vtcvubwtm1MrkREJkLTfGASTcEmSaFQ11Rr+t5Sfu6tt9gsTETUxDCx6suRI1Wbfx+lUgFpaepyRETUZDCx6ktWlm7LERGRSWBi1Rd3d92WIyIik8DEqi/9+qlH/9a2ObpIBEil6nJERNRkMLHqi1isnlIDVE+u5Y+XLVOXIyKiJoOJVZ/CwtRTatq0qXre05NTbYiImiiTWCDCpIWFAaNGqUf/ZmWp+1T79WNNlYioiWJiNQSxGBg4UOgoiIjIANgUTEREpENMrERERDrExEpERKRDTKxEREQ6xMRKRESkQ0ysREREOsTESkREpENMrERERDrExEpERKRDTKxEREQ6xMRKRESkQ0ysREREOsTESkREpENMrERERDrExEpERKRDJpNYR44cCS8vL0gkEri7u2PSpEnIzMzU+PqXXnoJIpEIy5Yt01+QRETU7JlMYg0NDcXWrVuRmJiIyMhIJCUlITw8XKNrd+zYgb///hseHh56jpKIiJo7c6ED0NTs2bMr/u7t7Y358+dj9OjRKC0thYWFRa3XZWRkYNasWThw4ACefvppQ4RKRETNmMnUWCvLycnBxo0b0bdv3zqTqlKpxKRJkzBnzhx06tTJgBESEVFzZVKJdd68ebC1tYWTkxNSU1Oxc+fOOssvXrwY5ubmeOONNzR+D7lcDplMVuUgIiLSlKCJddGiRRCJRHUecXFxFeXnzJmD+Ph4HDx4EGKxGJMnT4ZKparxtU+fPo3ly5dj7dq1EIlEGscUEREBBweHikMqlTb65yxXXKrQ2WsREZFxEqlqy0wGkJ2djezs7DrL+Pj4QCKRVDufnp4OqVSK2NhYBAcHV3t+2bJlePvtt2Fm9s93B4VCATMzM0ilUiQnJ9f4fnK5HHK5vOKxTCaDVCpFXl4e7O3tNfzJqsq4X4T3t59HcnYB/npnIMzMNE/0RERkHGQyGRwcHOrNB4IOXnJ2doazs3ODri3/PlA5CVY2adIkDB48uMq5YcOGYdKkSZg2bVqtr2tlZQUrK6sGxVSbVjaWOJOSC1lxGQ5fu4uB/q46fX0iIjIeJtHHevLkSaxcuRIJCQlISUlBdHQ0JkyYAF9f3yq11YCAAGzfvh0A4OTkhMDAwCqHhYUF3Nzc4O/vb9D4rS3FeLaHJwBg49+pBn1vIiIyLJNIrNbW1oiKisKgQYPg7++P6dOnIzAwEIcOHapSu0xMTEReXp6AkdZuYm8vAMCfl28jK69I4GiIiEhfTGIea+fOnfHXX3/VW66+7uLa+lUNob2rHXq3bYW/b+Zg88k0vD2kg2CxEBGR/phEjbWpeKGPNwBgy8lUlCqUAkdDRET6wMRqQMM6ucG5hSXu5Mvx5+XbQodDRER6YBJNwU2FpbkZZg/pALFIhH5+LkKHQ0REesDEamATe3sLHQIREekRm4KJiIh0iIlVACVlSmw4noznVsWiqITLHBIRNSVMrAIQm4nww5EbOJWci93nNN+snYiIjB8TqwDEZqKKvtafT6QIHA0REekSE6tAng+SwtLcDOfS85CQdl/ocIiISEeYWAXSytYS/3rcHQCw4ThrrURETQUTq4AmPVyJafe5TOQWlAgcDRER6QITq4C6Sh0R2MYeJWVKbI1LEzocIiLSASZWAYlEIkwJ9sHgx1qju3dLocMhIiId4MpLAnsuSIrngqRCh0FERDrCGisREZEOMbEaibScQkTsu4ykuw+EDoWIiBqBidVIfLT7Ir4/fINTb4iITBwTq5GYFOwDAIg8nY4H8jJhgyEiogZjYjUS/do7o52zLfLlZdgenyF0OERE1EBMrEbCzEyEScHqBSPWxyZDpVIJHBERETUEE6sRCe/hCVtLMa7deYDjSfeEDoeIiBqAidWI2EksENbdEwCwNjZZ2GCIiKhBmFiNzJS+PrCTmMOrlU3V5mCFAoiJATZvVv+p4AbpRETGiCsvGZn2ri1w6v3BkFiI/zkZFQW8+SaQnv7POU9PYPlyICzM8EESEVGtWGM1QtWSanh41aQKABkZ6vNRUYYNjoiI6mQyiXXkyJHw8vKCRCKBu7s7Jk2ahMzMzDqvmTp1KkQiUZWjT58+Boq4cVQqFf6+fheHPvseqGmEcPm5t95iszARkRExmcQaGhqKrVu3IjExEZGRkUhKSkJ4eHi91w0fPhxZWVkVx759+wwQbeNFncnA2J9O4j9dxqDWiTcqFZCWBhw5YsjQiIioDibTxzp79uyKv3t7e2P+/PkYPXo0SktLYWFhUet1VlZWcHNzM0SIOjW0U2vYRqqQ5CTFUZ+u6JecUHvhrCyDxUVERHUzmRprZTk5Odi4cSP69u1bZ1IFgJiYGLi6uqJDhw6YOXMm7ty5U2d5uVwOmUxW5RCCncQCz3mqf7Y1PUbWXdjd3QARERGRJkwqsc6bNw+2trZwcnJCamoqdu7cWWf5ESNGYOPGjfjrr7/w5Zdf4tSpU3jyySchl8trvSYiIgIODg4Vh1Qq3F6pU54NBgBE+wbhZkuP6gVEIkAqBfr1M3BkRERUG0ET66JFi6oNLnr0iIuLqyg/Z84cxMfH4+DBgxCLxZg8eXKdS/+NHTsWTz/9NAIDA/HMM89g//79uHr1Kvbu3VvrNQsWLEBeXl7FkZaWptOfWRttW9tjkKMCKpEZ1vZ4puqTIpH6z2XLALG42rVERCQMkUrARWmzs7ORnZ1dZxkfHx9IJJJq59PT0yGVShEbG4vg4GCN39PPzw8zZszAvHnzNCovk8ng4OCAvLw82Nvba/w+unLsejYm/vQ3bEqLcfybKXCQF6ifkErVSZXzWImIDELTfCDo4CVnZ2c4Ozs36Nry7wN1Nes+6t69e0hLS4O7CfVJ9vV1QoCbHYpLbZD2y044FN5R96n268eaKhGRETKJUcEnT57EyZMnERISgpYtW+LGjRv48MMP4evrW6W2GhAQgIiICIwZMwYPHjzAokWL8Oyzz8Ld3R3Jycl477334OzsjDFjxgj402hHJBJh7bRecLWzgpmZSOhwiIioHiaRWK2trREVFYWFCxeioKAA7u7uGD58OLZs2QIrK6uKcomJicjLywMAiMVinD9/HuvXr8f9+/fh7u6O0NBQ/PLLL7CzsxPqR2kQN4fqTeFERGScBO1jNQVC97FWVlyqwOGrdzG0k+nNyyUiMnWa5gOTmm7TnMnLFAj9Igb/3nAaCWn3hQ6HiIhqwcRqIqzMxQj2dQIA/HTkhsDREBFRbZhYTciMkHYAgP0XbiE9t1DgaIiIqCZMrCako4c9nmjvBIVShbXHkoUOh4iIasDEamLKa61bTqUhv7hU4GiIiOhRTKwmZkAHF7R3bYEH8jL8ckq45RaJiKhmTKwmxsxMhBdD2gIArt7OFzgaIiJ6lEksEEFVjenWBt28HBHgJuy8WiIiqo41VhMksRAzqRIRGSkmVhN3K68YN7MLhA6DiIgeYmI1YdtOp6Pf53/hk72XhA6FiIgeYmI1Yd28HFGmVOGPy3dw/Q4HMhERGQMmVhPm69ICgx9rDQD48fBNgaMhIiKAidXkvdRfvWDE9vgM3JEVCxwNERExsZq4IJ9W6O7liBKFEuuOJwsdDhFRs8fE2gT8+2Gt9ecTqSiQlwkcDRFR88bE2gQM6egGHycblJQpcS49T+hwiIiaNa681ASIzUT4enx3eLa0RktbS6HDISJq1phYm4jOng5Ch0BERGBTcJN0Nu0+VCqV0GEQETVLTKxNiEqlwtQ1JzHqm2OITrwjdDhERM0SE2sTIhKJ0KG1HQBgVcwNgaMhImqemFibmOlPtIWFWISTyTk4nZIrdDhERM0OE2sT4+YgwZhubQAAqw4lCRwNEVHzYzKJdeTIkfDy8oJEIoG7uzsmTZqEzMzMeq+7fPkyRo4cCQcHB9jZ2aFPnz5ITU01QMTC+Xd/X4hEwO+XbuPabS7OT0RkSCaTWENDQ7F161YkJiYiMjISSUlJCA8Pr/OapKQkhISEICAgADExMTh79iz+7//+DxKJxEBRC6O9awsM7ahenP871lqJiAxKpDLReRm7du3C6NGjIZfLYWFhUWOZcePGwcLCAhs2bGjw+8hkMjg4OCAvLw/29vYNfh1DO5t2H6O+OYZ2LrbY/2Y/WJmLhQ6JiMikaZoPTKbGWllOTg42btyIvn371ppUlUol9u7diw4dOmDYsGFwdXVF7969sWPHDsMGK5AuUkf8/GJvHHyrP5MqEZEBmVRinTdvHmxtbeHk5ITU1FTs3Lmz1rJ37tzBgwcP8Nlnn2H48OE4ePAgxowZg7CwMBw6dKjW6+RyOWQyWZXDVIX4OcNcbFL/xEREJk/QT91FixZBJBLVecTFxVWUnzNnDuLj43Hw4EGIxWJMnjy51hWGlEolAGDUqFGYPXs2unbtivnz5+Nf//oXVq1aVWtMERERcHBwqDikUqluf2gBqBfnvy90GEREzYKgfazZ2dnIzs6us4yPj0+Ng43S09MhlUoRGxuL4ODgas+XlJTA1tYWCxcuxAcffFBxft68eTh69CiOHTtW4/vJ5XLI5fKKxzKZDFKp1OT6WMul5RTi+e+PI7+4DMfmPwkH65qbzomIqG6a9rEKugi/s7MznJ2dG3Rt+feBykmwMktLS/Ts2ROJiYlVzl+9ehXe3t61vq6VlRWsrKwaFJMxauNoDXuJBbLyirHheDJmPekndEhERE2aVolVpVLh0KFDOHLkCJKTk1FYWAgXFxd069YNgwcP1luz6cmTJ3Hy5EmEhISgZcuWuHHjBj788EP4+vpWqa0GBAQgIiICY8aMAaBuOh47diz69++P0NBQ/Pbbb9i9ezdiYmL0EqcxMjMT4ZWBvnjrlwSsPnoT00PawsaSmxoREemLRn2sRUVF+PTTTyGVSjFixAjs3bsX9+/fh1gsxvXr17Fw4UK0bdsWTz31FE6cOKHzIK2trREVFYVBgwbB398f06dPR2BgIA4dOlSldpmYmIi8vH82+h4zZgxWrVqFzz//HJ07d8ZPP/2EyMhIhISE6DxGY/avx93h7WSD3MJSbPq7aS+OQUQkNI36WKVSKXr37o2pU6di2LBhNU5xSUlJwaZNm7Bq1Sp88MEHmDlzpl4CNjRTncf6qF9OpWJe5Hm42FnhyNxQSCw4BYeISBua5gONEuuFCxcQGBio0RuXlJQgJSUFfn5Noy+vqSTWkjIlQr+IQcb9Inw8qhMmBfsIHRIRkUnR6QIRgYGBSEhI0OiNLS0tm0xSbUoszc3w8oB2AIALGaY7N5eIyNhpPI+1e/fu6NGjB7777rsq/ZhkOp4LkmLnK8FY7JwDbN4MxMQACoXQYRERNSkaJ9Zjx46he/fumD9/Ptzd3fHCCy8gOjpan7GRjkl270SXkC5AaCgwYYL6Tx8fICpK6NCIiJoMjRNrcHAwfvzxR9y6dQvfffcd0tPTMXjwYPj6+uKTTz5Benq6PuOkxoqKAsLDgYf/TndsW+Ksmx+QkaE+z+RKRKQTjVp5KSkpCWvWrMH69euRlZWFIUOGYN++fbqMT3BNYvCSQqGumT5Mqkd8uuLFZz+EZ95t/L76NYihAjw9gZs3ATFHCxMR1cQgu9v4+vpi/vz5eP/992Fvb48DBw405uVIX44cqUiqANAtMxE2pcW44STFnoB+gEoFpKWpyxERUaM0OLEeOnQIU6ZMgZubG+bOnYuwsLBa198lgWVlVXnYoqQIM07tAACseGIcFCKzGssREZH2tEqsaWlp+Pjjj+Hr64vQ0FAkJSXh66+/RmZmJn788Uf06dNHX3FSY7i7Vzs15fRuOBTlI8lJin3+T9RajoiItKNxYh0yZAjatm2Lb7/9FuHh4bh8+TKOHj2KadOmwdbWVp8xUmP166fuQxWJKk7ZlRThxTj1frZf9x0HpdRLXY6IiBpF48RqbW2NyMhIpKenY/HixfD399dnXKRLYjGwfLn675WS65TTu2FX/ABXXbyx/8PlHLhERKQDGifWXbt2YdSoURA//PC9fv06Dhw4gKKiIgCodcNxMhJhYcC2bUCbNhWnHOQFmH4tBlYiFW4FdBEwOCKipkPr6Tb37t3D888/j+joaIhEIly7dg3t2rXDiy++CEdHR3z55Zf6ilUQTWK6TWUKhXr0b1YW4O4OWc8+KFYArvbVN5MnIqJ/6G26zezZs2FhYYHU1FTY2NhUnB87dix+++23hkVLhiMWAwMHAuPHAwMHwt5WwqRKRKRDWu94ffDgQRw4cACenp5Vzvv5+SElJUVngZHhnU7JhbxUgb7tnYUOhYjIZGldYy0oKKhSUy2XnZ1dZdNxMi074jPw7Hex+GDnBSiU7C8nImoorRNr//79sX79+orHIpEISqUSS5YsQWhoqE6DI8MZ9JgrHG0scONuAXadzRA6HCIik6V1U/CSJUswcOBAxMXFoaSkBHPnzsXFixeRk5PDlZdMmJ3EAjP7tcOSA4lY/sc1PPO4B8zFjVrxkoioWdL6k7Njx444d+4cevXqhSFDhqCgoABhYWGIj4+Hr6+vPmIkA5nS1wetbC2RfK8QUfGstRIRNUSjdrdpDprcdJt6/HA4CZ/uuwLPltb4652BsDRnrZWICNDxdJvU1FSt3jwjg7UdUzWpjw9c7KyQnluEX+LShA6HiMjkaJRYe/bsiZkzZ+LkyZO1lsnLy8OPP/6IwMBARHHTbJNlbSnGrND2aONojZY2FkKHQ0RkcjQavHT58mV8+umnGD58OCwsLBAUFAQPDw9IJBLk5ubi0qVLuHjxIoKCgrBkyRKMGDFC33GTHo3v5YVxvaSwMufawURE2tKqj7W4uBj79u3DkSNHkJycjKKiIjg7O6Nbt24YNmwYAgMD9RmrIJpbHysREdVM03zAwUv1aM6JtUyhRFR8BkoVSkzs7S10OEREgtLbWsFCGTlyJLy8vCCRSODu7o5JkyYhMzOzzmtEIlGNx5IlSwwUtWn7/dJtzN12Dp/tv4L7hSVCh0NEZBJMJrGGhoZi69atSExMRGRkJJKSkhAeHl7nNVlZWVWO//3vfxCJRHj22WcNFLVpG9rJDf6t7ZBfXIZVh24IHQ4RkUkw2abgXbt2YfTo0ZDL5bCw0Gz06ujRo5Gfn48///xT4/dpzk3BAPDHpduYsT4OEgszHJoTitbcCYeImqkm1xRcWU5ODjZu3Ii+fftqnFRv376NvXv34sUXX6yznFwuh0wmq3I0Z4Mec0V3L0cUlyqx4s9rQodDRGT0GrS7jVDmzZsHW1tbODk5ITU1FTt37tT42nXr1sHOzg5hYWF1louIiICDg0PFIZVKGxu2SROJRJg3PAAA8MupNCRnC/fvT0RkCrROrK1bt8b06dNx9OjRRr/5okWLah1gVH7ExcVVlJ8zZw7i4+Nx8OBBiMViTJ48GZq2ZP/vf//DxIkTIZHU3ZS5YMEC5OXlVRxpaVx9qHc7Jwzo4IIypQpf/X5V6HCIiIya1n2su3fvxtq1a7Fnzx54e3tj+vTpmDx5Mjw8PLR+8+zsbGRnZ9dZxsfHp8ZkmJ6eDqlUitjYWAQHB9f5GkeOHEH//v2RkJCALl26aBVjc+9jLXchIw//2XMJc4f5I8inldDhEBEZnN7nsd67dw/r16/H2rVrcenSJQwbNgzTp0/HyJEjYW6u9W50WktLS4OXlxeio6MxcODAOstOnToVFy5cqFL71RQTKxERAQYYvOTk5ITZs2fj7Nmz+Oqrr/DHH38gPDwcHh4e+PDDD1FYWNjQl67m5MmTWLlyJRISEpCSkoLo6GhMmDABvr6+VWqrAQEB2L59e5VrZTIZfv31V8yYMUNn8RCgVJrkYHIiIr1rcGK9desWPv/8czz22GOYP38+wsPD8eeff2Lp0qXYvn07Ro8erbMgra2tERUVhUGDBsHf3x/Tp09HYGAgDh06BCsrq4pyiYmJyMvLq3Ltli1boFKpMH78eJ3F05zJiksRsf8yxv1wgsmViKgGWjcFR0VFYc2aNThw4AA6duyIGTNm4IUXXoCjo2NFmYsXL6Jbt24oKTH91XrYFFxV9gM5BnwejYISBb4e3w3PdNG+b52IyBTprSl42rRp8PDwwLFjx5CQkIBZs2ZVSaoA0K5dO7z//vtaB03Gz7mFFf7d3xcAsORAIkrKlAJHRERkXLSusRYWFsLGxkZf8Rgd1lirK5CXYcCSGGQ/kOOjkZ0wpa+P0CEREemd3mqsZWVl1VYmkslkyM/PbxJNv1Q/WytzvDnYDwCw4s9ryC8uFTgiIiLjoXVidXR0RMuWLasdjo6OsLa2hre3NxYuXAilkk2ETdm4nlK0c7bFvYISrDqUJHQ4RERGQ+vEunbtWnh4eOC9997Djh07sH37drz33nto06YNvvvuO/z73//GihUr8Nlnn+kjXjISFmIzzBuhXupw49+pKCpRCBwREVF1ydkFWH30pkHHg2i9ksO6devw5Zdf4vnnn684N3LkSHTu3Bnff/89/vzzT3h5eeGTTz7Be++9p9NgybgM7dgac4b5I6x7G1hbioUOh4ioms/2X8FvF2/h2u18fPbs4wZ5T61rrMePH0e3bt2qne/WrRuOHz8OAAgJCUFqamrjoyOjJhKJ8Fpoe7g7WAsdChFRNSdv5uC3i7dgJgKmh7Q12PtqnVg9PT2xevXqaudXr15dsRPMvXv30LJly8ZHRyblUqZM400RiIj0SalU4ZO9lwAAY3t6oUNrO4O9t9ZNwV988QWee+457N+/Hz179oRIJMKpU6dw5coVbNu2DQBw6tQpjB07VufBkpFQKIAjR4CsLMDdHaqQELy17Tx2JmRi9ZQgDHqstdARElEzt/tcJs6m58HWUoy3h3Qw6HtrnVhHjhyJq1evYtWqVUhMTIRKpcKIESOwY8cO+Pj4AABeeeUVXcdJxiIqCnjzTSA9veKUyNMTbu+uBGCOT/ZdRv8OLrAQN3i1TCKiRikuVeDz3xIBAK8M9IWLnVU9V+iWVom1tLQUQ4cOxffff4+IiAh9xUTGKioKCA8HHm3uzcjAa/MmYtu7W3DjbgE2nkjB1CcM159BRFTZ3nNZyLhfBHcHCV4MaWfw99cqsVpYWODChQsQiUT6ioeMlUKhrqnW1IeqUsG+pBBvx27C+70mYukf1zC6Wxs42lgaPk4iavbCureBrZUYZiKRIDMWtG6vmzx5co2Dl6iJO3KkSvNvNSoVxsb8Av8WIuQVlWL5n9cMFxsRUSUikQjDA90xtJObIO+vdR9rSUkJfvrpJ/z+++8ICgqCra1tlee/+uornQVHRiQrq94i5iolPnB9gEkPbLHheApe6OMNX5cWBgiOiAhIyymEvcQCDjYWgsahdWK9cOECunfvDgC4evVqlefYRNyEubtrVKxfB1cMsrTFpSwZbsuKmViJyCBUKhXe+fUsrt3Ox4rx3dDPz0WwWLROrNHR0fqIg4xdv36ApyeQkVFzP6tIpH6+Xz9E9CiFnZUFV2MiIoPZf+EWTt7MgcTCDO0E/kLf4DkR169fx4EDB1BUVAQAXBigqROLgeXL1X9/tGWi/PGyZYBYDFc7CZMqERlMcakCn+67DAD4d39ftHEUdjU4rRPrvXv3MGjQIHTo0AFPPfUUsh72vc2YMQPvvPOOzgMkIxIWBmzbBrRpU/W8p6f6fFhYldNKpQpbT6Xhl1Nc3pKI9OenIzeQnlsEN3sJXh5g+Ok1j9I6sc6ePRsWFhZITU2tsuH52LFj8dtvv+k0ODJCYWFAcjIQHQ1s2qT+8+bNakkVUDfNzI08h//uuYzsB3LDx0pETV7m/SJ8E63eunLBUwGwsdS6h1PntI7g4MGDOHDgADw9Pauc9/PzQ0pKis4CIyMmFgMDB9ZbbHigGzp52ONipgxfHEg02M4SRNR8ROy/gqJSBXr6tMTILh5ChwOgATXWgoKCKjXVctnZ2bCyMuyyUWTcxGYifDSyEwDgl7g0nE/PEzgiImpKFEoVWliJYW4mwqKRnYxmZorWibV///5Yv359xWORSASlUoklS5YgNDRUp8GR6QvyaYVRXT2gUgELd12AUslBbkSkG2IzESLCHseReaHo5OEgdDgVtG4KXrJkCQYOHIi4uDiUlJRg7ty5uHjxInJycnDs2DF9xEgmbsGIx/D7pds4k3ofUfEZCO/hWf9FREQaqnNP6Ed240K/furuLD3SusbasWNHnDt3Dr169cKQIUNQUFCAsLAwxMfHw9fXVx8xkolzc5DgjUF+AIDFv11BcalC4IiIyJTdeyDHa5vO4GZ2Qd0Fo6IAHx8gNBSYMEH9p4+P+rweiVQmMgF15MiRSEhIwJ07d9CyZUsMHjwYixcvhodH7Z3VDx48wPz587Fjxw7cu3cPPj4+eOONN7Ta1k4mk8HBwQF5eXmwt7fXxY/SLJWUKfHOr2cx/QkfdPNqKXQ4RGTC5m47i61x6ejm5YioV/rW3Lda225c5WVrmCJYH03zQYMS6/3793Hy5EncuXMHSqWyynOTJ0/W9uU0snTpUgQHB8Pd3R0ZGRl49913AQCxsbG1XjNz5kxER0fjp59+go+PDw4ePIhXX30VkZGRGDVqlEbvy8RKRGQ8Tqfk4tnv1J/7ka/0RQ/vGr6oKxTqmmltG4eUrxR386ZWzcJ6S6y7d+/GxIkTUVBQADs7uyrfFEQiEXJycrR5uQbbtWsXRo8eDblcDguLmhdcDgwMxNixY/F///d/Fed69OiBp556Ch9//LFG78PEqh8Z94vg4SAxmlF8RGT8FEoVnvn6KC5lyfB8kCc+D+9Sc8GYGHWzb32iozWaOlhO03ygdR/rO++8g+nTpyM/Px/3799Hbm5uxWGopJqTk4ONGzeib9++tSZVAAgJCcGuXbuQkZEBlUqF6OhoXL16FcOGDTNInFSzb6KvI3RJDKLOZAgdChGZkI1/p+BSlgz2EnPMGx5Qe0ENduPSqpyWtE6sGRkZeOONN2qcy6pv8+bNg62tLZycnJCamoqdO3fWWX7FihXo2LEjPD09YWlpieHDh+Pbb79FSEhIrdfI5XLIZLIqB+mWmUiEEoUSn+67jPuFJUKHQ0QmIPuBHEsOJAIA5gwPgFOLOtZN0HA3Lo3LaUnrxDps2DDExcXp5M0XLVoEkUhU51H5vebMmYP4+HgcPHgQYrEYkydPrnPx/xUrVuDEiRPYtWsXTp8+jS+//BKvvvoq/vjjj1qviYiIgIODQ8UhlUp18rPSP14MaQs/1xa4V1CCzx/+RyEiqsuPh28gv7gMgW3sMaGXV92Fy3fjqq2rSSQCpFJ1OT3Quo919erV+M9//oNp06ahc+fO1ZpiR44cqfFrZWdnIzs7u84yPj4+kEgk1c6np6dDKpUiNjYWwcHB1Z4vKiqCg4MDtm/fjqeffrri/IwZM5Cenl7rusZyuRxy+T/r2spkMkilUvax6tiJG/cw7ocTEImAqFf6cqQwEdVJXqbAD4duoH8HF3SROtZ/QfmoYKDqyGADjArWeoGImTNnAgD+85//VHtOJBJBodB8jqKzszOcnZ21DQHAP9vUVU6ClZWWlqK0tBRmZlUr5WKxuNpI5sqsrKy4NKMB9GnnhLDubRB1JgPvb7+AXbOegLm4wbsYElETZ2UuxusP58NrpHw3rjffrDo62NNTvcWllklVG1p/kimVyloPbZKqNk6ePImVK1ciISEBKSkpiI6OxoQJE+Dr61ulthoQEIDt27cDAOzt7TFgwADMmTMHMTExuHnzJtauXYv169djzJgxeomTtPPeU4/BwdoCl7JkWH+cGzgQUXWnU3JRqqi9MlQnLXbj0iXh99fRgLW1NaKiorBw4UIUFBTA3d0dw4cPx5YtW6rULhMTE5GX989C71u2bMGCBQswceJE5OTkwNvbG5988glefvllIX4MeoRzCyvMHe6P/+y+VGtXCBE1XzezCzD+xxNo52yLzTP7oKWtpfYvouFuXLqkcR/rU089hc2bN8PBQb3Q8SeffILXXnsNjo6OANQboPfr1w+XLl3SW7BC4DxW/VIqVcjMK4JnS8OPMici46VSqTBp9UkcvZ6Nfn7OWD+9l+Dz3nU+j/XAgQNV+jMXL15cZd5qWVkZEhM5wpO0Y2YmYlIlomq2x2fg6PVsWJmb4b+jAwVPqtrQOLE+WrE1kSWGyYTEJefghZ/+hqy4VOhQiEhA9x7I8fEedevnm4P94O1kK3BE2uEwTDIKCqUKcyPP4ej1bCz5jS0fRM3Zf/deRm5hKQLc7DCzXzuhw9Gaxom1fMGGR88R6YLYTIT/jgoEAPz8dwpOp+QKHBERCeHw1bvYHp8BkQj47NnHYWGC0/A0HhWsUqkwderUilG4xcXFePnll2Frq66i1zaflEhTfds7I7yHJ7adTsf8yHPY80YIrMz1uyExERmXNi2t0cunFTq1sUdXTRaCMEIajwqeNm2aRi+4Zs2aRgVkbDgq2LDuF5Zg8FeHkP2gBG8M8sPbQzoIHRIRGZhSqUJpaRmsjh9TL5Tv7q5eflCLLd70Qa/7sTYnTKyGt+98Fl7deAbmZiLsfj0Ej7nzvhM1dcWlCkgsHibOqKiaV0xavlzvizvURW/bxhHp24hANwzr1BplShW2na5lo2IiajLkZQqM/uYYPtp9EYW/Plzj99FNyjMy1OejooQJUgsmsfISNVEKBXDkSLWmHpFIhI9HBWJIRzc8272N0FESkZ5989d1XLmVj+x8Od5Y+x5sampIVanUC+i/9RYwapTgzcJ1YY2VhBEVBfj4AKGhwIQJ6j99fCq+jbraSxDew5Mjz4mauEuZMnwbkwQA+MhPhJY36phup1IBaWnqL+RGjImVDC9Ku6YeWXEpvom+DoWSwwGImpJShRJzI8+iTKnCsE6t8ZTyrmYXZmXpN7BGYmIlw1Io1IMSamvqAdRNPQ93SlIoVXj221gsOZCINcduGi5OItK7VTFJuJAhg4O1BT4eFQiRh7tmF7prWE4gTKxkWEeOVK+pVvZIU4/YTIRpT7QFACw5kIib2QWGiJKI9OxSpgwr/roGAPhoZCe42kvU4yw8PVHrdlciESCVqssZMSZWMixNm3AqlRvfS4qQ9s6Qlykx59ezbBImagIy7xdBYiHG0I6tMaqrh/qkWKyeUgNUT67lj5ctM+qBSwATKxmapk04lcqJRCJ89mxntLAyR1xKLpuEiZqAwR1b4/fZA/BpWOeqgxTDwoBt24A2j8wI8PRUnxdwHqumuEBEPbhAhI4pFOrRvxkZNfezikTq/0A3b1b7Vrr5ZCoWRJ2HpbkZ9r4eAr/WdoaJmYgMr5bpeELiAhFknBrR1DOupxQD/V1QUqbEp/su6zdOItK54lIFJv50An9duV1/YbEYGDgQGD9e/aeRN/9WxsRKhtfAph6RSITPn30czwd54qvnu+o/TiLSqS8OJOLY9XuYF3kehSVlQoejN2wKrgebgvXICJt6iEg/YpOyMeHHvwEA/5sahCcDWgsckfY0zQdc0pCEU97U00AqlQq/XbiFgf6usLZkQiYyVrLiUry79SwAYHwvL5NMqtpgUzCZrI92X8IrG8/gk32XhA6FiOqwaNdFZOYVw9vJBh88/ZjQ4egdEyuZrEGPuQIAfj6Rij8vazAYgogMQ6EAYmKAzZuxZ/MfiDqTATMR8NXzXWBr1fQbSplYyWT183PBiyHqVZnmbjuHu/lygSMiokc32Di9eisA4FX3MvTwbiVsbAbCwUv14OAl41Zcqt7H8cqtfIT6u+B/U3tyRxwioZRvsPFIWoluF4SQlARYbP3FJBZ4qA3nsVKzILEQY/m4brA0N0N04l1sOJEidEhEzVMdG2yE3oiDhVJRZYONpsxkEuvIkSPh5eUFiUQCd3d3TJo0CZmZmXVec/v2bUydOhUeHh6wsbHB8OHDce3aNQNFTIbi72aHBSMCAAAR+67g3gM2CRMZXKUNNs66+WFG2Ae4a+P4z/MmspeqLphMYg0NDcXWrVuRmJiIyMhIJCUlITw8vNbyKpUKo0ePxo0bN7Bz507Ex8fD29sbgwcPRkEBd0hpaqb29cHYIClWTwmCUwsrocMhan4ebpyRb2mN10fOxR9+ffBVv4m1lmvKTLaPddeuXRg9ejTkcjksLCyqPX/16lX4+/vjwoUL6NSpEwBAoVDA1dUVixcvxowZMzR6H/axEhFpICYGqtBQvPHMHOzuOABt8m5j35o34CB/pCITHd2o+etCatJ9rDk5Odi4cSP69u1bY1IFALlc3RwokUgqzonFYlhaWuLo0aO1vrZcLodMJqtykOm5mV2AQ1fvCh0GUfPRrx9+GTAWuzsOgFipwIpdS6omVRPZS1UXTCqxzps3D7a2tnByckJqaip27txZa9mAgAB4e3tjwYIFyM3NRUlJCT777DPcunULWXU0RURERMDBwaHikEql+vhRSI8uZOThXyuOYNbGM0i9Vyh0OETNwtXsQizq+wIA4N3DG9Aj88o/T5rQXqq6IGhiXbRoEUQiUZ1HXFxcRfk5c+YgPj4eBw8ehFgsxuTJk1FbS7aFhQUiIyNx9epVtGrVCjY2NoiJicGIESMgruMfdsGCBcjLy6s40tLSdP5zk375u9khwN0e+fIyvLbpDORlTX8UIpGQikoUmLXpDIqVIvSzV+KljL+rFjChvVR1QdA+1uzsbGRnZ9dZxsfHp0pzbrn09HRIpVLExsYiODi4ztfIy8tDSUkJXFxc0Lt3bwQFBeGbb77RKEb2sZqmzPtFeGrFEdwvLMWUYG98NCpQ6JCImqzUe4WYsuYkHsjLsO+NfnCxMW+SG2yYxCL8zs7OcHZ2btC15d8HyvtS6+Lg4AAAuHbtGuLi4vDxxx836D3JdHg4WmPp810xbe0prDuegl5tnfD04+5Ch0XUJHk52WDXrCeQmlMIF7uHo/JNdICSLphEH+vJkyexcuVKJCQkICUlBdHR0ZgwYQJ8fX2r1FYDAgKwffv2ise//vorYmJiKqbcDBkyBKNHj8bQoUOF+DHIwEIDXPHyAF8AwLzIc0jO5jQrIl0qVSgr/m4nsUAnDwcBozEeJpFYra2tERUVhUGDBsHf3x/Tp09HYGAgDh06BCurf+YsJiYmIi8vr+JxVlYWJk2ahICAALzxxhuYNGkSNm/eLMSPQAJ5d2gH9PRpiQfyMnx/OEnocIiajPziUjy94gj+d/RmrWNdmiuTncdqKOxjNQH1bJh+K68YG04k463BHWAhNonvkkRGTaVSYdameOw9nwV3BwkOzO4Pe0nNUx+bEpPoYyVqtKgo9fqkD5dSA6Aegbh8ecUIRDcHCeYMCxAoQKKm56cjN7H3fBbMzURYOaF7s0iq2uDXdzJd5TtpVE6qAJCRoT4fFVXtklKFEp/uu4yEtPuGiZHIlFTaRxUxMTUumB+blI3PflPPUf3wmY7o4d3SsDGaACZWMk117KRRca6GnTRW/nUdPxy+gVd/Ps3F+okqe2QfVYSGqh9X+oKaeb8Ir2+Kh0KpQlj3NpjUx1uwcI0ZEyuZpko7adSolp00ZvRri3bOtsjMK8arG89UGdVI1Gxp0PpTUqbEKxvP4F5BCTq62+PTMZ2593EtmFjJNGm6Q8Yj5ewkFvh+Ug/YWorx980cfLL3sh6CIzIhGrb+WECJZx53h3MLS3w/qQckFqa/4IO+MLGSaXLXcLGHGsr5tbbD0rFdAQBrY5OxNY7LVlIzpmHrj+joUczo1w7R7w6EtJWN4eIzQUysZJr69VOP/q2tKaqenTSGdnLDW4P9AAAfbL+A+NRcfUVKZNzqaf250NoX+ZbWFeXsOAK4XkysZJrEYvWUGqB6ctVwJ403nvTD0I6tYSEW4X5RqX7iJDJ2dbT+pDq0xgtjP8aYSV8i07G1AYMybUysZLrCwtQ7ZrRpU/W8hjtpmJmJ8NXYroh69QmE+rvqMVAiI1ZL60++pTVmPPsh7lvbw1akRKvQEIECND1cIIJMW1gYMGpUg3fSaGFlDn83u4rHd2TFcG5hBTMzjnakZqK89Sc8XJ1cVSqUiczw+sh5uOriDdcH9/D9IHdIJJZCR2oyWGMl0ycWq3fSGD9e/WcDt6eKS87BiOVHsOyPqzoNj8joPdL6898nZyDGNwiSshL81MsObuOaxz6qusIaK9FDyfcKca+gBCv+uo62LrYY081T6JCIDOdh68+GjX9h7aUSAMDSSb3weJc29VxIj2KNleih8B6eeGlAOwDA3G3ncOLGPYEjIjKsYiXwwy11i8+cYf4YwaTaIEysRJXMGxaApzq7oVShwr/Xx+H6nXyhQyIyGImFGJGv9MWcYf54daCv0OGYLCZWokrMzET46vmu6O7lCFlxGaauOYW7+VxTmJo2hfKfVZdc7SR4LbQ9lytsBCZWokdILMT4cXIQvJ1skJ5bhBV/XhM6JKK6abArTW3yikox+ptj2BGfobfwmhsmVqIaOLWwwpqpPRHewxPvP/2Y0OEQ1U6DXWlqIy9T4JWfT+N8Rh4i9l9GgbxM7+E2B0ysRLVo59ICXzzXpcpi46qaFionEkoD9iQup1Cq8PYvZxGbdA+2lmL8b2pP2FpxooguMLESaUClUmHJgSv4/ECi0KEQqTVwT2L10yos2nURe89nwUIswveTgtDJw0G/8TYj/HpCpIFTybn4JjoJAODcwgovhrQVOCJq9rTZk3jgwCpPrfjzOjacSIFIBCwd2xUhfs76jbWZYY2VSAO92rbCnGH+AICP91ziVnMkvAbuSXzixj0sfbi62EcjO+Ffj3voOrJmjzVWIkDdXFbPesOvDvTF/cIS/HjkJuZHnoOdlTlGdNZwX1giXWvgnsS927bCS/3bwcrcDJODfXQfF0Gk4miMOslkMjg4OCAvLw/29vZCh0P6EBWl7quq3Kzm6alemPyRHXJUKhXmR57HL3FpsBCLsHpKT/Tv4GLggImg/jLo46MeqFTTx7hIpP49vnmz2pfE8o99zlXVjqb5gE3B1LxpOapSJBLh07DOeLqzO0oVKrzy82lkP+ACEiQALfYkjkm8g1d+Po3iUsXDp0VMqnpkcolVLpeja9euEIlESEhIqLOsSqXCokWL4OHhAWtrawwcOBAXL140TKBk/Bo4qlJsJsLSsV0xKMAV/xkVCOcWVvqPlagmGuxJHJuUjZd/Po39F25h9dGbwsTZzJhcYp07dy48PDTrbP/888/x1VdfYeXKlTh16hTc3NwwZMgQ5Odz/VeCdqMqH2FpboafpgTh2R6elYqzV4UEEBYGJCcD0dHApk3qP2/eBMLCcOLGPUxfewrFpUo8GeCKmf3aCR1ts2BSiXX//v04ePAgvvjii3rLqlQqLFu2DO+//z7CwsIQGBiIdevWobCwEJs2bTJAtGT0GjiqslzlprQ7smKErzqO8+l5uoiMSDs17En89417mLZGnVQH+rvg24ndYWluUh/5Jstk7vLt27cxc+ZMbNiwATY2NvWWv3nzJm7duoWhQ4dWnLOyssKAAQMQGxtb63VyuRwymazKQU1UA0dV1mTxb4k4nZKLiT+dwNm0+42Li6iRTiXnYNraUygqVaB/BxeseqFHlRXESL9MIrGqVCpMnToVL7/8MoKCgjS65tatWwCA1q1bVznfunXriudqEhERAQcHh4pDKpU2PHAybv36qfuiahvEIRIBUqm6XD0+GtUJQd4tISsuwws//Y3TKbk6DpaarEYsoF+T4lIFXt8Uj8ISBfr5OeOHSUyqhiZoYl20aFHF6LTajri4OHz99deQyWRYsGCB1u/x6Mg3lUpV52i4BQsWIC8vr+JIS+NCAE2WFqMq69PCyhzrpvdC77atkC8vw+TVf+NUco5u46WmpxEL6NdGYiHGNxO7Y3gnN/wwKYhJVQCCzmPNzs5GdnZ2nWV8fHwwbtw47N69u0pCVCgUEIvFmDhxItatW1ftuhs3bsDX1xdnzpxBt27dKs6PGjUKjo6ONV5TE85jbQZqmscqlaqT6iPzWOtTWFKGGeviEJt0D9YPt5/jcnFUo/KpXo9+BJd/zj0c1aupAnkZF9HXM03zgUksEJGamlqlrzMzMxPDhg3Dtm3b0Lt3b3h6ela7RqVSwcPDA7Nnz8bcuXMBACUlJXB1dcXixYvx0ksvafTeTKzNhAYrL2mquFSBf284jcNX7yLAzQ573+gHsRnnDFIl5Ys71DYqvY7FHWqy62wmPtp1Eeum90JgGy6mry+a5gOT+Hrj5eVV5XGLFi0AAL6+vlWSakBAACIiIjBmzBiIRCK89dZb+PTTT+Hn5wc/Pz98+umnsLGxwYQJEwwaP5mA8lGVOqDeKL0HPtl7GS8N8GVSpeoasYD+ozacSMGHOy9ApQJ+OZXGxGoETCKxaioxMRF5ef9Md5g7dy6Kiorw6quvIjc3F71798bBgwdhZ2cnYJTUHFiZi/GfUYFVzl25JUOAG1s9CI2e6gWoW+W+/us6vvpdvaD+pD7e+GhkJ11ER41kEk3BQmJTMOnCbxdu4dWNpzGzfzvMGxYAM9Zim7eYGPVApfpER9dYY1UoVfh4zyWsjU0GALwxyA+zB/txmUI941rBREYk5V4BlCrg+0M38M6vZ1FSphQ6JBJSI6Z6FZcqMGvTmYqkuvCZjnh7SAcmVSPCxEpkAC8N8MWS8MchNhNhe3wGXlx3CvnFpUKHRUJpxFQvsZkIBSUKWIrNsGJ8N0x7oq1+YyWtsSm4HmwKJq3UM7o4OvEOXtt4BoUlCgS42WH11J5o42gtYMAkqAZO9XogL0PirXz08G6p/xipQpOabiMkJlbSmIb7up5Lv48X18Xhbr4czi2ssP/NfnCx4w45JkdXU7Q0eJ2/b9zDX1fuYP6IADb5CqhJTbchMnq1TfYv39e10mT/xz0dsfO1JzB97Sn0aefEpGqKNPwSpZF6pnptjUvD+9vPo1Shgr+bHcK6V5+3T8aFNdZ6sMZK9WrgZP8CeRmszM1gLlYPdcgvLoWtpTlHDBs7Ha+YVJsyhRKfH0jED4dvAACe7uyOL57rAmtLLlEoFI4KJjKUBu7ramtlXpFUSxVKzFwfhxnr4yDjoCbjpVCoa6o11UfKz731VqMX0s8pKMHUNacqkuqbg/zw9fhuTKomgomVqLF0MNn/YqYM8an38deVOxj9zTFcv5Ovo+BIpxr4JUobFzLy8MzXR3H0ejZsLMVYOaEbZg/pwJYME8LEStRYOtjXtavUEdte7gt3Bwlu3C3AyJXHsDMhQ0cBks7o4EtUfXILS5CVVwRvJxtsf/UJ/Otxjwa/FgmDiZWosXS0r2tnTwfsfj0EfX2dUFiiwJtbEvDBjvOQlzWuWZF0SAdfompSeahLPz8XfDuxO3a9FgJ/Ny6/aoqYWIkaS4f7ujq3sMKGF3vj9SfbAwB+PpGKBZHndRgsNYqOvkRVdu12PsZ8G4ukuw8qzg0PdIeDjUVjoyWBMLES6UJYmHo0aJs2Vc97emo9SlRsJsI7Q/2xZlpPuDtI8MpAXx0HSw2mwy9RKpUKP59IwTMrjyIh7T4+2n1Jt7GSYDjdph6cbkNa0eG+rgBQUqaEpfk/339/u5CFYF9nOFizNtMojf13auCKSeWyH8gxP/Ic/rh8BwDQz88ZS8d2hXMLzmk2Zlx5SUeYWEkQNXzwx6Xl4fnvj8PVToLF4Y9jQAcXoaM0Tbpa3KGByfnPy7cxL/I8sh/IYSk2w7wRAZjW14ejfk0AE6uOMLGSwdXywX/ukxV4805L3MwuAACM7+WF959+DC2suICaxgy0uENt/rh0GzPWxwEA/FvbYdm4rnjMnZ8rpoKJVUeYWMmg6vngL/plGxZb+VdsGdbG0RqfhnVm7VUTDVwhS5fKFEqErzqOXm1b4e0hHSCx4IIPpoSJVUeYWMlgtPjgP558H3O2nUV6bhEA4IU+Xvjv6M6Gi9UUNXJz8Ya4lVeMVYeSMH9EQEUSLVUoYSHmuFFTxCUNiUyNFqv6BPs64cBb/fFiSFuYicDmRE0YYHGHcgqlCmuO3cTgrw5hbWwyvo1JqniOSbXpY+cMkbHQ8oPf1soc//evjng+SAo/1xYVTx+7ng1rSzG6ezXBvTobM5pXT4s7POrvG/fw0e5LuJQlAwB083LE8E5ujXpNMi1MrETGooEf/JVX5ykqUWDutnPIuF+EZ7t5YI5DLtxyb+tk6o/gGjuat3xxh4yMmhfRL29q12Jxh8rScwsRsf8K9p5Tf/Gxk5hj3vAATOjlxRG/zQzbJIiMhQ5W9SkpUyLY1wkAEBmfiYG/5+KL7/Yhf9hT6v7bqCg9BG4A5YO6Hm0qL9/vVpOfS4eLO9TkiwOJ2HsuC2YiYGJvL8S8OxAv9PFmUm2GmFiJjIUOPvgdbCzwhTgJURveRVD6RRRbSLCy7zgMeOkn/M+tO4rHjje95KrLrdp0uEJWUYkCOQUlFY/fHNwBIe2dsfv1EHwypjOcuNhDs8VRwfXgqGAyuMas6lNpZLEKwEG/Plg8YCpuOHkCACJ/fhc9zAq0m1Kii9WkGvMa+hjN24h4CkvKsPFEKr4/fAP9/Zzx1diumr0nmTxN8wH7WImMTVgYMGpUwz74K40sFgEYdu0Enkw6hV87D8ZZ9w7okXGlotzRNoHo6uVY9wITulilqLGvoY/RvGKx1lNqZMWl2HgiFT8euVFRU41LyUVhSRlsLPlRSv/gbwORMWrABz+AGpOLhVKBCWcPYMLZAxXn7qbewvQ/imBlboYX+nhjWl8fuNpLql5Y22IV5f2amjSd6uI1DDSatzbpuYVYcywZv5xKwwN5GQDAq5UNZoW2x5jubTh9hqoxud8IuVyOrl27QiQSISEhoc6yUVFRGDZsGJydnTUqT2TyNEwuWQ4u8GxpjfziMnwXk4S+n/2FWZvO4O8b99R7g+qiX1NXfaN62KpNG3vPZWH10Zt4IC+Dn2sLfPFcF/z1zgA831PKpEo1Mrnfirlz58LDw0OjsgUFBXjiiSfw2Wef6TkqIiOhYRJ6/F8D8cfsAfhhUg8EebdEmVKFPeeyMPaHExix/Agu7Y3ReLGKWmmx4EWd9Dyat7LsB3L8ePgGDl68VXFufG8vDH7MFWum9cTB2f0R3sMT5kyoVAeTagrev38/Dh48iMjISOzfv7/e8pMmTQIAJCcn6zkyIiNRnoTCw9VJp3Jt8ZEkZAZgaCc3DO3khouZedhwPAU7EjKQdPcB3J1kFZelOLqh9YMcSMpKUE1d/Zq67BstH81bU1+thlu11aaoRIHfL9/G9jPpOHwtGwqlCt28HDH04aIO9hIL/DSlZ4Nfn5ofk0mst2/fxsyZM7Fjxw7Y2Njo7X3kcjnkcnnFY5lMVkdpIiPUgCTUycMBnz37OBY89RjOpOSi5a1/Nt1+61/vItHFG08mncLTV46i/80zsC0tVj9ZV9OzrvtGGzOoqwb7zmdhz7lMxCTeRWHJP83RXTwd8FwPKVQqFUS11fyJ6mASiVWlUmHq1Kl4+eWXERQUpNcaaEREBD766CO9vT6RQTQwCTlYWyA0wBXwcwI8PVF4Jxt3bR1RaGmNPY/1x57H+sNCUYpeaRcxMOc6ngzoBt/aXkwfKx01cFCXUqnC1Tv5CHD7Z4pE1Jn0io3Gpa2sMaZrG4zq1ga+Li1qexkijQg6j3XRokX1JrFTp04hNjYWv/zyCw4fPgyxWIzk5GS0bdsW8fHx6Nq1a73vo035mmqsUqmU81ip+Xk4olcF4KybH/b7P4H9HfoiteU/NcxRXT2wfFw3AOot0S5kyhDgZvfPdmjlo4KBmpul9bT/aXGpAhczZTidkoO45FycTsnFvYISHJozEN5OtgCA3y/dxtm0+xjWyQ2BbexZO6V6mcS2cdnZ2cjOzq6zjI+PD8aNG4fdu3dX+cVXKBQQi8WYOHEi1q1bV+draJuIK+MCEdSsPTIHVQXgZsceiH7t/xBj5YZRXdsgvId68YlLmTI8teIIzM1E8Gtth0APe7R3bYH218/D94v/oM3Vc7BQPmxy1XTBi3rIy9SvZ2WuTuS/XbiFpb9fxfW7D6BQVv1os7EU4+vx3TDosdaNek9qvkwisWoqNTW1Sl9nZmYmhg0bhm3btqF3797w9PSs83omVqJG0HCVokNX7+LtXxJwr6CGQU4AFgZYYppVNuDujsQOXbHu7zQ4WlvA0cYCNpbmsDQ3g5W5GSzEZujcxgHSVuqxFDezC/DbhVvILy6FrLgUd2Ry3M6X446sGLdlxfhmQneM6KyuRf9+6TZmro8DADjZWqK7d0sEebdEkE9LdG7jCEtzjualhmtSKy95eXlVedyihboPxNfXt0pSDQgIQEREBMaMGQMAyMnJQWpqKjIzMwEAiYmJAAA3Nze4uXEbJyKNaNivOaCDC+I+GIysvGKcS8/D5SwZku4+QNLdAty4+wDuQZ2BQPX/u6TzWdj0d2qtr/V5+OMVifXG3QdY/NuVWsum5BRW/L2nT0usnhKEjh72cLOXsHmXBGESiVVTiYmJyMvLq3i8a9cuTJs2reLxuHHjAAALFy7EokWLDB0eUZMnEong4WgND0drDA/858urUqmCslLjmK9LC7w12A/3C0txv7AEhSUKlCiUKClTolShhJOtZUVZr1Y2eLa7J+wk5rC3toCLnRVa21mhtb0E7o4SuFRa7N7RxpJNvSQ4k2gKFhKbgomICNA8H7DDgYiISIeYWImIiHSIiZWIiEiHmFiJiIh0iImViIhIh5hYiYiIdIiJlYiISIeYWImIiHSIiZWIiEiHmFiJiIh0iImViIhIh5rUIvz6UL6UcuVt64iIqPkpzwP1LbHPxFqP/Px8AIBUKhU4EiIiMgb5+flwcHCo9XnublMPpVKJzMxM2NnZNZm9HWUyGaRSKdLS0rhjj4Z4z7THe9YwvG/aM9Q9U6lUyM/Ph4eHB8zMau9JZY21HmZmZlU2U29K7O3t+R9XS7xn2uM9axjeN+0Z4p7VVVMtx8FLREREOsTESkREpENMrM2QlZUVFi5cCCsrK6FDMRm8Z9rjPWsY3jftGds94+AlIiIiHWKNlYiISIeYWImIiHSIiZWIiEiHmFiJiIh0iIm1iYmIiEDPnj1hZ2cHV1dXjB49GomJiXVek5WVhQkTJsDf3x9mZmZ46623DBOskWjIPYuKisKQIUPg4uICe3t7BAcH48CBAwaKWHgNuWdHjx7FE088AScnJ1hbWyMgIABLly41UMTGoSH3rbJjx47B3NwcXbt21V+QRqYh9ywmJgYikajaceXKFYPEzMTaxBw6dAivvfYaTpw4gd9//x1lZWUYOnQoCgoKar1GLpfDxcUF77//Prp06WLAaI1DQ+7Z4cOHMWTIEOzbtw+nT59GaGgonnnmGcTHxxswcuE05J7Z2tpi1qxZOHz4MC5fvowPPvgAH3zwAX744QcDRi6shty3cnl5eZg8eTIGDRpkgEiNR2PuWWJiIrKysioOPz8/A0TM6TZN3t27d+Hq6opDhw6hf//+9ZYfOHAgunbtimXLluk/OCOl7T0r16lTJ4wdOxYffvihHqMzTg29Z2FhYbC1tcWGDRv0GJ3x0ua+jRs3Dn5+fhCLxdixYwcSEhIME6SR0eSexcTEIDQ0FLm5uXB0dDRsgGCNtcnLy8sDALRq1UrgSExHQ+6ZUqlEfn5+s73PDbln8fHxiI2NxYABA/QVltHT9L6tWbMGSUlJWLhwoSHCMmra/K5169YN7u7uGDRoEKKjo/UdWgUuwt+EqVQqvP322wgJCUFgYKDQ4ZiEht6zL7/8EgUFBXj++ef1GJ1x0vaeeXp64u7duygrK8OiRYswY8YMA0RpfDS9b9euXcP8+fNx5MgRmJs3749sTe+Zu7s7fvjhB/To0QNyuRwbNmzAoEGDEBMTo1WLSkM173+lJm7WrFk4d+4cjh49KnQoJqMh92zz5s1YtGgRdu7cCVdXVz1GZ5y0vWdHjhzBgwcPcOLECcyfPx/t27fH+PHj9Ryl8dHkvikUCkyYMAEfffQROnToYMDojJOmv2v+/v7w9/eveBwcHIy0tDR88cUXBkmsUFGTNGvWLJWnp6fqxo0bWl03YMAA1ZtvvqmfoIxcQ+7Zli1bVNbW1qo9e/boMTLj1dDfs3Iff/yxqkOHDjqOyvhpet9yc3NVAFRisbjiEIlEFef+/PNPA0UsvMb+rv33v/9VBQQE6DiqmrHG2sSoVCq8/vrr2L59O2JiYtC2bVuhQzJ6Db1nmzdvxvTp07F582Y8/fTTeo7SuOjq90ylUkEul+s4OuOl7X2zt7fH+fPnq5z79ttv8ddff2Hbtm3N4v+3rn7X4uPj4e7uruPoasbE2sS89tpr2LRpE3bu3Ak7OzvcunULgHpzXmtrawDAggULkJGRgfXr11dcVz7C8MGDB7h79y4SEhJgaWmJjh07GvxnMLSG3LPNmzdj8uTJWL58Ofr06VNxjbW1tUYbIZu6htyzb775Bl5eXggICACgntf6xRdf4PXXXxfmhxCAtvfNzMysWl+iq6srJBJJsxk30ZDftWXLlsHHxwedOnVCSUkJfv75Z0RGRiIyMtIwQRukXkwGA6DGY82aNRVlpkyZohowYEC913l7exs0dqE05J4NGDCgxmumTJli8PiF0JB7tmLFClWnTp1UNjY2Knt7e1W3bt1U3377rUqhUBj+BxBIQ/9/VrZw4UJVly5d9B6rsWjIPVu8eLHK19dXJZFIVC1btlSFhISo9u7da7CYOY+ViIhIhziPlYiISIeYWImIiHSIiZWIiEiHmFiJiIh0iImViIhIh5hYiYiIdIiJlYiISIeYWImIiHSIiZWIiEiHmFiJqEH27NmDdu3aoWfPnrh69arQ4RAZDS5pSEQN0qFDB3z77be4ePEijh8/ji1btggdEpFRYI2ViGp17949uLq6Ijk5udpzzs7OaN++Pdq1a1dtR5/w8HB89dVXBoqSyLiwxkpEtXr33XeRm5uL1atXV3vuxx9/xMsvvwwXFxdcvHgRTk5OFc+dO3cOoaGhuHnzJuzt7Q0ZMpHgWGMlauZKSkpqPF9UVITVq1djxowZ1Z4rKyvD8uXLMXfuXOTn56Nly5ZVnn/88cfh4+ODjRs36iVmImPGxErUzAwcOBCzZs3C22+/DWdnZwwZMqTGcvv374e5uTmCg4OrPbdq1Sq0a9cOr732GgoLC3Ht2rVqZUaOHInNmzfrPH4iY8fEStQMrVu3Dubm5jh27Bi+//77GsscPnwYQUFB1c7n5ubi448/xuLFi+Hp6QkHBwckJCRUK9erVy+cPHkScrlc1+ETGTUmVqJmqH379vj888/h7++PgICAGsskJyfDw8Oj2vkPP/wQY8aMwWOPPQYA6NixI86ePVutXJs2bSCXy3Hr1i3dBk9k5MyFDoCIDK+mmuijioqKIJFIqpy7dOkSfv75Z1y+fLniXGBgYI01VmtrawBAYWFh44IlMjFMrETNkK2tbb1lnJ2dkZubW+Xc7Nmzcf/+fXh6elacUyqVcHd3r3Z9Tk4OAMDFxaWR0RKZFiZWIqpRt27d8PPPP1c83rNnD06fPo34+HiYm//z0XHq1ClMnz4dd+/erZJEL1y4AE9PTzg7Oxs0biKhMbESUY2GDRuGBQsWIDc3Fy1atMA777yDOXPmoGvXrlXKlc9TPXv2LAYPHlxx/siRIxg6dKghQyYyCkysRFSjzp07IygoCFu3bkVBQQHu37+PWbNmVSsnlUphY2ODhISEisRaXFyM7du348CBA4YOm0hwXHmJiGq1b98+vPvuu7hw4QLMzDSfRPDNN99g586dOHjwoB6jIzJOrLESUa2eeuopXLt2DRkZGZBKpRpfZ2Fhga+//lqPkREZL9ZYiYiIdIgLRBAREekQEysREZEOMbESERHpEBMrERGRDjGxEhER6RATKxERkQ4xsRIREekQEysREZEOMbESERHpEBMrERGRDv0/ARZl33KyrpoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 500x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#make a vector to evaluate fits on with a lot of points so it looks smooth\n",
    "rfit = np.linspace(min(r),max(r),100)\n",
    "a,b,c = polyfit(r,e,2) #this is from pylab\n",
    "\n",
    "plot(r,e,'ro')\n",
    "plot(rfit, a*rfit**2 + b*rfit + c,'--',label='parabolic fit')\n",
    "xlabel('r ($\\AA$)')\n",
    "ylabel('Energy (eV)')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "862ab803",
   "metadata": {},
   "source": [
    "　２次関数を用いたフィッティング結果をみると、平衡原子間距離 $r_0 \\sim 2.3$ を中心にして、赤丸のデータは対象な形をとっていないことがわかる。\n",
    " \n",
    "　赤丸の計算値は、右側にやや平らな形状をしている。これは原子間ポテンシャルが調和振動子では記述しきれない**非調和性**を持っていることを示している。熱膨張はこのような原子間ポテンシャルの非対称性（非調和性）によって現れる。\n",
    " \n",
    " \n",
    "　温度が上昇すると平衡原子間距離 $r_0 \\sim 2.3$ から振動するのであるが、\n",
    " \n",
    "* rが小さい領域のポテンシャルは急峻になっている\n",
    "\n",
    "* rが大きい領域のポテンシャルはなだらかになっている\n",
    "\n",
    "　これによってrが大きい領域へ平均原子間距離が動き、熱膨張を生じるということになる。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25a76c3c",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Morse ポテンシャルでフィッティング\n",
    "\n",
    "　この様な固体のミクロの振舞いを記述するモデルを考える。固体中の原子間に働く力を数学的に記述する2体間ポテンシャル (pair potential) にはさまざまなものがあるが，簡単なものとしてモース (Morse) 型がある。\n",
    " \n",
    "　Morseポテンシャルは，2原子分子の原子間距離$r$に依存するポテンシャルである.\n",
    "\n",
    "$$\n",
    "E_{\\text{total}}(r) = A + D\\exp\\lbrack - 2\\lambda(r - r_{0})\\rbrack - 2D\\exp\\lbrack - \\lambda(r - r_{0})\\rbrack\n",
    "$$\n",
    "\n",
    "　ここで$r$は原子間距離。$r_0$は平衡原子間距離。$A$、$D$、$\\lambda$はフィッティングパラメータである。Cuのエネルギーと体積の関係に、モース型ポテンシャルをフィッティングしてみよう。\n",
    "\n",
    "```\n",
    "*****here*****\n",
    "```\n",
    "に正しい式を記入して実行してください。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "289089ac",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#difinition of the functions\n",
    "def Morse(parameters,r):\n",
    "\n",
    "    A = parameters[0]\n",
    "    D = parameters[1]\n",
    "    r0 = parameters[2]\n",
    "    la = parameters[3]\n",
    "    \n",
    "    *****here*****\n",
    "\n",
    "    return E\n",
    "\n",
    "def objective(pars,y,x):\n",
    "    #we will minimize this function\n",
    "    err =  y - Morse(pars,x)\n",
    "    return err"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f624b3f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#now here are our initial guesses.\n",
    "A = -4.1\n",
    "D = 1\n",
    "r0 = 2.3\n",
    "la = 0.1\n",
    "\n",
    "x0 = [A, D, r0, la] \n",
    "\n",
    "morpars, ier = leastsq(objective, x0, args=(e,r)) #this is from scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9230e41",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(r,e,'ro')\n",
    "plot(rfit, Morse(morpars,rfit), label='Morse fit')\n",
    "xlabel('r ($\\AA$)')\n",
    "ylabel('Energy (eV)')\n",
    "legend(loc='best')\n",
    "\n",
    "#add some text to the figure in figure coordinates\n",
    "ax = gca()\n",
    "text(0.4,0.5,'A = %1.2f eV' % morpars[0],transform = ax.transAxes)\n",
    "text(0.4,0.4,'D = %1.2f eV' % morpars[1],transform = ax.transAxes)\n",
    "text(0.4,0.3,'r0 = %1.2f $\\AA$' % morpars[2],transform = ax.transAxes)\n",
    "text(0.4,0.2,'lambda = %1.2f $\\AA^{-1}$' % morpars[3],transform = ax.transAxes)\n",
    "savefig('a-eos.png')\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d4e64e8",
   "metadata": {},
   "source": [
    "　正しい結果は以下のようになります。\n",
    "\n",
    "<img src='fig/Morsefit.png' width=\"500\">\n",
    "\n",
    "二次関数の時よりもフィッティング結果が良好であることがわかる。ここで、\n",
    "\n",
    "* $r_0$の平衡原子間距離は2.29 Å\n",
    "\n",
    "* $A$=-1.31 eV、$D$=2.79 eV、$\\lambda$=1.80 Å$^{-1}$\n",
    "\n",
    "である。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aea026f0",
   "metadata": {},
   "source": [
    "## 課題\n",
    "\n",
    "　The Materials Projectからダウンロードした原子間ポテンシャルのデータをMoodleからダウンロードしてください（EOS.xlsx）。Jupyter notebookを使用して、これらのデータを読み込み、各元素の原子間ポテンシャルのグラフを作成しましょう。その結果から、それぞれの原子間の安定な原子間距離と最安定なエネルギーを求めてください。\n",
    " \n",
    "- 提出ファイルはJupyter notebook。\n"
   ]
  }
 ],
 "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",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}