laplacian-eigenstates.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy.special import spherical_jn as j_l\n",
    "from scipy.special import spherical_in as i_l\n",
    "from spherical_bessel_zeros import Jn_zeros\n",
    "from scipy.integrate import quadrature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "l_big = 26\n",
    "n_big = 26\n",
    "a = 6.0  # Radius of the sphere"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6 13\n"
     ]
    }
   ],
   "source": [
    "z_ln = Jn_zeros(l_big, n_big)  # Spherical Bessel zeros\n",
    "z_nl = z_ln.T\n",
    "\n",
    "E_nl = z_nl**2\n",
    "E_max = E_nl[6-1, 0]\n",
    "n_max = np.where(E_nl[:, 0] <= E_max)[0][-1] + 1\n",
    "l_max = np.where(E_nl[0, :] <= E_max)[0][-1]\n",
    "print(n_max, l_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def R_nl(n, l, r):\n",
    "    return j_l(l, z_nl[n, l]*r/a)\n",
    "\n",
    "def N_nl(n, l):\n",
    "    # Normalization factor for LE basis functions\n",
    "    def function_to_integrate_to_get_normalization_factor(x):\n",
    "        return j_l(l, x)**2 * x**2\n",
    "    integral, _ = sp.integrate.quadrature(function_to_integrate_to_get_normalization_factor, 0.0, z_nl[n, l])\n",
    "    return (1.0/z_nl[n, l]**3 * integral)**(-0.5)\n",
    "\n",
    "def get_LE_function_python(n, l, r):\n",
    "    R = np.zeros_like(r)\n",
    "    for i in range(r.shape[0]):\n",
    "        R[i] = R_nl(n, l, r[i])\n",
    "    return N_nl(n, l)*R*a**(-1.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Feed LE (delta) radial spline points to Rust calculator:\n",
    "\n",
    "n_spline_points = 101\n",
    "spline_x = np.linspace(0.0, a, n_spline_points)  # x values\n",
    "\n",
    "spline_f = []\n",
    "for l in range(l_max+1):\n",
    "    for n in range(n_max):\n",
    "        spline_f_single = get_LE_function_python(n, l, spline_x)\n",
    "        spline_f.append(spline_f_single)\n",
    "spline_f = np.array(spline_f).T\n",
    "spline_f = spline_f.reshape(n_spline_points, l_max+1, n_max)  # f(x) values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_LE_function_derivative(n, l, r):\n",
    "    delta = 1e-6\n",
    "    all_derivatives_except_at_zero = (get_LE_function_python(n, l, r[1:]+delta) - get_LE_function_python(n, l, r[1:]-delta)) / (2.0*delta)\n",
    "    derivative_at_zero = (get_LE_function_python(n, l, np.array([delta/10.0])) - get_LE_function_python(n, l, np.array([0.0]))) / (delta/10.0)\n",
    "    return np.concatenate([derivative_at_zero, all_derivatives_except_at_zero])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "spline_df = []\n",
    "for l in range(l_max+1):\n",
    "    for n in range(n_max):\n",
    "        spline_df_single = get_LE_function_derivative(n, l, spline_x)\n",
    "        spline_df.append(spline_df_single)\n",
    "spline_df = np.array(spline_df).T\n",
    "spline_df = spline_df.reshape(n_spline_points, l_max+1, n_max)  # df/dx values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"splines.txt\", \"w\") as file:\n",
    "    np.savetxt(file, spline_x.flatten(), newline=\" \")\n",
    "    file.write(\"\\n\")\n",
    "\n",
    "with open(\"splines.txt\", \"a\") as file:\n",
    "    np.savetxt(file, (1.0/(4.0*np.pi))*spline_f.flatten(), newline=\" \")\n",
    "    file.write(\"\\n\")\n",
    "    np.savetxt(file, (1.0/(4.0*np.pi))*spline_df.flatten(), newline=\" \")\n",
    "    file.write(\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rascaline import SphericalExpansion\n",
    "\n",
    "hypers_spherical_expansion = {\n",
    "        \"cutoff\": a,\n",
    "        \"max_radial\": int(n_max),\n",
    "        \"max_angular\": int(l_max),\n",
    "        \"center_atom_weight\": 0.0,\n",
    "        \"radial_basis\": {\"Tabulated\": {\"file\": \"splines.txt\"}},\n",
    "        \"atomic_gaussian_width\": 100.0,\n",
    "        \"cutoff_function\": {\"Step\": {}},\n",
    "    }\n",
    "\n",
    "calculator = SphericalExpansion(**hypers_spherical_expansion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(100, 1, 6)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGdCAYAAAAMm0nCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABY3ElEQVR4nO3dd1zV9eLH8dc5h42CA2UoigNHLlQUcGSDQrNhWam30syGM41yNTRvmWZqZpqroU21oXXNSCO1VHDj3rkVUFNQkHm+vz/6RZerlkfRLwfez8fjPO7ly+f7Oe/vudzOu++0GIZhICIiIuJErGYHEBEREXGUCoyIiIg4HRUYERERcToqMCIiIuJ0VGBERETE6ajAiIiIiNNRgRERERGnowIjIiIiTsfF7ABFwW63c/z4ccqWLYvFYjE7joiIiFwBwzA4d+4cQUFBWK2O7VMpEQXm+PHjBAcHmx1DRERErsKRI0eoWrWqQ+uUiAJTtmxZ4I8PwMfHx+Q0IiIiciXS09MJDg4u+B53RIkoMH8eNvLx8VGBERERcTJXc/qHTuIVERERp6MCIyIiIk5HBUZEREScTok4B0ZERC7NMAzy8vLIz883O4qUUjabDRcXlyK/zYkKjIhICZWTk8OJEyfIzMw0O4qUcl5eXgQGBuLm5lZkc6rAiIiUQHa7nQMHDmCz2QgKCsLNzU03+pQbzjAMcnJyOHnyJAcOHCA0NNThG9ZdjgqMiEgJlJOTg91uJzg4GC8vL7PjSCnm6emJq6srhw4dIicnBw8PjyKZVyfxioiUYEX1b7si1+J6/B3qL1tEREScjgqMiIjIFbBYLCxcuBCAgwcPYrFYSEpKMjVTaaYCIyIi4qDg4GBOnDhBw4YNzY7yj2bPnk25cuWuy9zLly/HYrFw9uzZ6zL/31GBERGRYi0nJ8fsCBex2WwEBATg4qJrYcxyVQVm6tSphISE4OHhQUREBGvXrr3s2G+++Ybw8HDKlSuHt7c3YWFhfPLJJ4XGGIbBiBEjCAwMxNPTk+joaPbu3Xs10S5pb8o5Js5fQsKMASR8NJTEz19n7deT2LD4A7au+Ia9m37h6NEjpF3IxW43iux9RUTEcbfccgv9+/dn0KBB+Pn5ERMTA8DEiRNp1KgR3t7eBAcH07dvX86fP1+w3qFDh7jnnnsoX7483t7eNGjQgMWLFxf8fvv27dx99934+PhQtmxZ2rZty/79+wFYt24dd9xxB35+fvj6+tKuXTs2btx42Yz/ewjpzz0R8fHxhIeH4+XlRatWrdi9e3eh9b799luaNWuGh4cHNWvWZNSoUeTl5V32fex2O//+97+pWrUq7u7uhIWFERcXV/D7S+0BSUpKwmKxcPDgQZYvX07Pnj1JS0vDYrFgsVh49dVXAQgJCeG1116jW7dueHt7U6VKFaZOnXrZbQQ4e/YsFouF5cuXc/DgQW699VYAypcvj8Vi4fHHH7/sthQ1h6vjvHnziI2NZfr06URERDBp0iRiYmLYvXs3lStXvmh8hQoVeOmll6hXrx5ubm4sWrSInj17Urly5YI/ynHjxjF58mTmzJlDjRo1eOWVV4iJiWHHjh1FcrnVgVMZbEjaSKzbx5cdMya3GzPy78FmtdDcK4UxxmQyXCuQ7e5HnldlrL5BuFUIpqx/DSoG16Vc+Yq6p4KIOBXDMLiQa84deT1dbQ79M3POnDn06dOHVatWFSyzWq1MnjyZGjVq8Ntvv9G3b1+GDBnCe++9B0C/fv3Iycnhl19+wdvbmx07dlCmTBkAjh07xs0338wtt9zCzz//jI+PD6tWrSooD+fOnaNHjx68++67GIbBhAkTuOuuu9i7dy9ly5a94twvvfQSEyZMoFKlSvTu3ZsnnniiYBt+/fVXunfvzuTJkwvK09NPPw3AyJEjLznfO++8w4QJE5gxYwZNmzblww8/5N5772X79u2Ehob+Y55WrVoxadIkRowYUVCm/vxMAN566y1efPFFRo0axY8//sjAgQOpU6cOd9xxxz/OHRwczNdff03nzp3ZvXs3Pj4+eHp6/uN6RcXhAjNx4kSeeuopevbsCcD06dP5/vvv+fDDDxk2bNhF42+55ZZCPw8cOJA5c+awcuVKYmJiMAyDSZMm8fLLL3PfffcB8PHHH+Pv78/ChQvp2rXrVWxWYTX8vLmjZRiJx7tgzc3AlpeJLS8T1/wMPPPO4W0/R5qtPORDvt3AI/MEtdx+g/zfIAtIA078Nd+43If50PoAweW9CPM9T+fsb7FUqIGHf20qhTQioFooVpvtmnOLiBSlC7n53DTiR1Pee8e/Y/Byu/KvnNDQUMaNG1do2aBBgwr+e0hICK+//jq9e/cuKDCHDx+mc+fONGrUCICaNWsWjJ86dSq+vr7MnTsXV1dXAOrUqVPw+9tuu63Qe82cOZNy5cqxYsUK7r777ivOPXr0aNq1awfAsGHD6NixI1lZWXh4eDBq1CiGDRtGjx49CvK99tprDBky5LIFZvz48QwdOrTgu/DNN99k2bJlTJo0qdDekstxc3PD19cXi8VCQEDARb9v3bp1wXd3nTp1WLVqFW+//fYVFRibzUaFChUAqFy58nU7z+ZyHCowOTk5bNiwgeHDhxcss1qtREdHk5CQ8I/rG4bBzz//zO7du3nzzTcBOHDgAMnJyURHRxeM8/X1JSIigoSEhEsWmOzsbLKzswt+Tk9P/9v3DfUvS2in9kD7y44ZC7yam8/ZzFzOnGrA5sN1yD57gvxzqVjOJ+OWcYIy2SlUyEvluOFHVq6dvannqXJqE5Fu8yAV2AWsgCzDlWMuVTnrVYND1R/Eu/7t1A/woWp5T6xW7bUREfknzZs3v2jZTz/9xJgxY9i1axfp6enk5eWRlZVFZmYmXl5ePPvss/Tp04clS5YQHR1N586dady4MfDHYZW2bdsWlJf/lZKSwssvv8zy5ctJTU0lPz+fzMxMDh8+7FDuP98PIDAwEIDU1FSqVavG5s2bWbVqFaNHjy4Yk5+fX2gb/lt6ejrHjx+ndevWhZa3bt2azZs3O5TrcqKioi76edKkSUUy9/XmUIE5deoU+fn5+Pv7F1ru7+/Prl27LrteWloaVapUITs7G5vNxnvvvVfQ7pKTkwvm+N85//zd/xozZgyjRo1yJPoV8XC1EeBrI8A3BGqFXHbcm3n5DDqbxeHfM0k75Ebivm64nztM+awjBOUfx8OSS638A3DuAB9vrMe368sB0NZtDy+5zSPNty6WoDAqhkZQrV5zXN3ci3xbRET+l6erjR3/jjHtvR3h7e1d6OeDBw9y991306dPH0aPHk2FChVYuXIlvXr1IicnBy8vL5588kliYmL4/vvvWbJkCWPGjGHChAkMGDDgHw9t9OjRg9OnT/POO+9QvXp13N3diYqKcvgE4v8uSH8eMrPb7QCcP3+eUaNG8cADD1y03tWeLvHnDeIM46/zN3Nzc69qrhs5d1G4IadPly1blqSkJM6fP098fDyxsbHUrFnzosNLV2r48OHExsYW/Jyenk5wcHARpf1n7i42Qvy8CfHzhjp3wH/tasvLzeHIod2cPriNrOM78LO2psHv3uxNOU9o/n7q5e2E0zvh9ELYCtmGK3tca3CmXEPONepBaMMWVKvgpfNrRKTIWSwWhw7jFCcbNmzAbrczYcKEgi/W+fPnXzQuODiY3r1707t3b4YPH86sWbMYMGAAjRs3Zs6cOeTm5l5yL8yqVat47733uOuuuwA4cuQIp06dKtJtaNasGbt376Z27dpXNN7Hx4egoCBWrVpVcFjqz6wtW7YEoFKlSgCcOHGC8uXLA1x0bxo3N7fLPo08MTHxop/r169/0dxNmza97NyAKU87d+gv2c/PD5vNRkpKSqHlKSkplzy29ier1VrwP1hYWBg7d+5kzJgx3HLLLQXrpaSkFOxu+/PnsLCwS87n7u6Ou3vx3Gvh4upGcO1GBNf+4xhs5P8vz823c/RgTdbvakLekU2UObOdatl78LFkUidvD5zawwNxjdn4QwZ+Zdzo4neItl6HKHfTLdRq3EZ7aUSkVKtduza5ubm8++673HPPPaxatYrp06cXGjNo0CA6dOhAnTp1OHPmDMuWLSv4Mu7fvz/vvvsuXbt2Zfjw4fj6+pKYmEjLli2pW7cuoaGhfPLJJ4SHh5Oens7gwYOL/ITUESNGcPfdd1OtWjUefPBBrFYrmzdvZtu2bbz++uuXXGfw4MGMHDmSWrVqERYWxkcffURSUhKfffZZwecSHBzMq6++yujRo9mzZw8TJkwoNEdISEjBDoQmTZrg5eVVcLhq1apVjBs3jk6dOrF06VK+/PJLvv/+e+CPZxhFRkYyduxYatSoQWpqKi+//HKhuatXr47FYmHRokXcddddeHp6FjpJ+LoyHNSyZUujf//+BT/n5+cbVapUMcaMGXPFc/Ts2dNo166dYRiGYbfbjYCAAGP8+PEFv09LSzPc3d2NL7744ormS0tLMwAjLS3tijMUB/l5ecaRvVuNdYtmGglTnzY6v7vMqP3i90b1oYuMj196wDBG+hjGSB8jY0QlY8sb7YzVHw4xdiTGGTnZWWZHF5Fi7sKFC8aOHTuMCxcumB3FYe3atTMGDhx40fKJEycagYGBhqenpxETE2N8/PHHBmCcOXPGMAzD6N+/v1GrVi3D3d3dqFSpkvHYY48Zp06dKlh/8+bNxp133ml4eXkZZcuWNdq2bWvs37/fMAzD2LhxoxEeHm54eHgYoaGhxpdffmlUr17dePvttwvWB4wFCxYYhmEYBw4cMABj06ZNhmEYxrJlywplMQzD2LRpkwEYBw4cKFgWFxdntGrVyvD09DR8fHyMli1bGjNnzrzsZ5Gfn2+8+uqrRpUqVQxXV1ejSZMmxg8//FBozMqVK41GjRoZHh4eRtu2bY0vv/zyovft3bu3UbFiRQMwRo4caRiGYVSvXt0YNWqU8dBDDxleXl5GQECA8c477xSae8eOHUZUVJTh6elphIWFGUuWLDEAY9myZQVj/v3vfxsBAQGGxWIxevToccntuNzf47V8f1sMw3Doxifz5s2jR48ezJgxg5YtWzJp0iTmz5/Prl278Pf3p3v37lSpUoUxY8YAf5yvEh4eTq1atcjOzmbx4sUMGzaMadOm8eSTTwJ/nFU9duzYQpdRb9my5Yovo05PT8fX15e0tDR8fHwc2ZxiJys3n+3H00hb8znlD35PjcwtlON8oTEZhgfDqn1GWN1atKntRx3/MjrkJCKFZGVlceDAAWrUqFFkT/+VkiUkJIRBgwYVurrrernc3+O1fH87fDC0S5cunDx5khEjRpCcnFxwU50/T8I9fPhwoadOZmRk0LdvX44ePYqnpyf16tXj008/pUuXLgVjhgwZQkZGBk8//TRnz56lTZs2xMXFlcr/03m42mhevQJU7w/0x56fz4Fd60ndtgzXI6upcX4jvxtl+M/ebP6zdwcAE70+oqqvGy5176R25D34lKto7kaIiIhcZw7vgSmOStIemH9iz89nz/59rDjhwsp9p9h0IIU1tifxtvxxWXmuYWOP+02cC76NKpGdCQ5tYnJiETGD9sDIP3H2PTAqME4uKzubPYmLydweR9CpX6lmP1bo9z+53My65uO48yZ/woLLY9N9aERKBRUYKU6KxSEkKV483N1p3O5+aHc/AMd+28nRtQvxPLiUeheSWH8hiBkrfmPGit+oUSaPMeUWUqZpZ+q1jMHlMjd0EhERKe5UYEqYKjXrU6VmfWA46WdP03jvSe7dd4Flu1MJy/yVyLxvYOk3nF7qy76Kt+Dd9EHqR96FTU9UFRERJ6JvrRLMp1xF7mpRkbtaQE6enS1rXVi79hB1zv5CRdKoePpb+OlbTv5Unv2V76TsbbHcVLeurmgSEZFiTwWmlHBzsRLe6nZodTu5OdlsTVzMhU1fUffMMipxhoop82k1uw3uFU/QKawKnZtUplrlcmbHFhERuSQVmFLI1c2dRjffDzffT052Fkm/LuD4zkTSUitz4XQm78TvJeyXJ8nwyONcvYeod3t3XZotIiLFigpMKefm7kFYdDfCortxS04eS7an8MP6nbQ+ug233HzYup2sLaNZ79sOz4gnuCmqA5b/us+PiIiIGfRNJAW83Fzo1LQKM56K5uzTG0is+SwHrcF4WHIJT/+JBkv/xbHXbuKnuZNJTc8yO66ISCGvvvrqZZ+hdzm33HJLkd8H5UpyrFq1ikaNGuHq6kqnTp2K9P0d9fjjj5ue4WpoD4xcUuUqNajc/TUM+yj2JP3KmVUf0ODUEqpygqlbD/LMlp+Jrl+Zx1oG0aq2P1abzezIIlLKvfDCCwwYMMChdb755ptLPp36eouNjSUsLIwffvjhhj388ODBg9SoUYNNmzYVKljvvPMOznhLOBUY+VsWq5U6zdpBs3Zknk9j3dI5HD1en/wjOfy4PQW/XZ8R4vYDR2t1oV77PpTzu/xTyUVErgfDMMjPz6dMmTIOl4EKFSpcp1R/b//+/fTu3ZuqVaua8v7/zdfX1+wIV0WHkOSKeZXxpcX9z/JJvztY8tzN9IiqzgMuq6lqnCBy3yQ8323I2knd2LdltdlRRcSJZWdn8+yzz1K5cmU8PDxo06YN69atK/j98uXLsVgs/PDDDzRv3hx3d3dWrlx50aGbvLw8nn32WcqVK0fFihUZOnQoPXr0KHS45H8PIYWEhPDGG2/wxBNPULZsWapVq8bMmTML5Rs6dCh16tTBy8uLmjVr8sorr5Cbm3tF23bw4EEsFgunT5/miSeewGKxMHv2bGbPnk25cuUKjV24cGGh21r8uX2ffPIJISEh+Pr60rVrV86dO1cwxm63M27cOGrXro27uzvVqlVj9OjRANSoUQOApk2bYrFYuOWWW4CLDyFd6ecfHx9PeHg4Xl5etGrVit27d1/RZ1BUVGDkqtTxL8uo+xpS74UfWdtwJPtstXC35NLy7GJqf9OBHaNbsyFuDnn5drOjisj/ysm4/Cs3y4GxF65srIOGDBnC119/zZw5c9i4cSO1a9cmJiaG33//vdC4YcOGMXbsWHbu3Enjxo0vmufNN9/ks88+46OPPmLVqlWkp6ezcOHCf3z/CRMmEB4ezqZNm+jbty99+vQp9OVctmxZZs+ezY4dO3jnnXeYNWsWb7/99hVtW3BwMCdOnMDHx4dJkyZx4sSJQg83/if79+9n4cKFLFq0iEWLFrFixQrGjh1b8Pvhw4czduxYXnnlFXbs2MHnn39e8LDltWvXAvDTTz9x4sQJvvnmm0u+x5V+/i+99BITJkxg/fr1uLi48MQTT1zxdhQFHUKSa+JdthwtH4zFsA9i14afyfhlKo3TV3BT7jZ+WvURbTdV5fFWIXRtWQ1fTz26QKRYeCPo8r8LvRMe+fKvn9+qDbmZlx5bvQ30/P6vnyc1gszTF497Ne2Ko2VkZDBt2jRmz55Nhw4dAJg1axZLly7lgw8+YPDgwQVj//3vf3PHHXdcdq53332X4cOHc//9fzxqZcqUKSxevPgfM9x111307dsX+GNvy9tvv82yZcuoW7cuAC+//HLB2JCQEF544QXmzp3LkCFD/nFum81GQEAAFosFX19fAgIcO+xut9uZPXs2ZcuWBeCxxx4jPj6e0aNHc+7cOd555x2mTJlCjx49AKhVqxZt2rQBoFKlSgBUrFjxsu/ryOc/evRo2rVrB/xRJjt27EhWVtYNe/aWCowUCYvVSr0W0dAimtRjB9n/w2TmHa/GibQsxvywi3nxibwR8AvV7oolKKSu2XFFpJjav38/ubm5tG7dumCZq6srLVu2ZOfOnYXGhoeHX3aetLQ0UlJSaNmyZcEym81G8+bNsdv/fs/wf+/NsVgsBAQEkJqaWrBs3rx5TJ48mf3793P+/Hny8vJu2IOEQ0JCCsoLQGBgYEG2nTt3kp2dze23337V8zvy+f/35xQYGAhAamoq1apVu+r3d4QKjBS5ylVCqPzkRJrm5vNd0nHeX/kbD5z+nMjUReR/NI8NZdvhE/0CoWFtzY4qUjq9ePzyv7P8zxWFg/f9zdj/OQth0Narz3QVvL29r8u8/3tVksViKSg9CQkJPPLII4waNYqYmBh8fX2ZO3cuEyZMuKb3tFqtF10JdKnzav4um6en5zVlcNR/Z/nzXJ1/KodFSefAyHXj4Wrj4RbB/DjoZiLveJit7s2wWQyan19O6MK72TrmVrb9+i3GDfyDFxHAzfvyL1cPB8Z6XtlYB9SqVQs3NzdWrVpVsCw3N5d169Zx0003XfE8vr6++Pv7Fzr5ND8/n40bNzqU53+tXr2a6tWr89JLLxEeHk5oaCiHDh26pjnhj8M7586dIyPjr3OGkpKSHJojNDQUT09P4uPjL/l7Nzc34I/P4XKK6vO/EbQHRq47i8VCWLv7oN19/LZtDb8veYuwtHgaZW+E+O5s+aURx+6ZT0zDQKxWPUhSpDTz9vamT58+DB48mAoVKlCtWjXGjRtHZmYmvXr1cmiuAQMGMGbMGGrXrk29evV49913OXPmzDU9sDY0NJTDhw8zd+5cWrRowffff8+CBQuuer4/RURE4OXlxYsvvsizzz7LmjVrmD17tkNzeHh4MHToUIYMGYKbmxutW7fm5MmTbN++nV69elG5cmU8PT2Ji4ujatWqeHh4XHQJdVF+/teb9sDIDVWzYQThsV+R2jORNX6duWC4sf5CFfp8vomYSb/wbdIx8vLyzI4pIiYaO3YsnTt35rHHHqNZs2bs27ePH3/8kfLlyzs0z9ChQ+nWrRvdu3cnKiqKMmXKEBMTc00nmd57770899xz9O/fn7CwMFavXs0rr7xy1fP9qUKFCnz66acsXryYRo0a8cUXX/Dqq686PM8rr7zC888/z4gRI6hfvz5dunQpOEfGxcWFyZMnM2PGDIKCgrjvvvsuOUdRff7Xm8Vwxtvv/Y/09HR8fX1JS0u7YSdSSdH4PfUYc9cdYdq6NM5l5dHMsofJHjM43rgvTTs+g6ubu9kRRZxSVlYWBw4coEaNGjfsqpDizm63U79+fR5++GFee+01s+OUKpf7e7yW728dQhJTVahchb4dq/Do7bl8knCIWiumUtU4QdXNr3B8y7scbdiPpvf0UZEREYcdOnSIJUuW0K5dO7Kzs5kyZQoHDhzgX//6l9nRpAjoEJIUCz4ervS7tTZtn/+cxNqDOEU5goxUWm4dSeqYRqxbMJm83ByzY4qIE7FarcyePZsWLVrQunVrtm7dyk8//UT9+vXNjiZFQIeQpFjKyjxP0oKJhO59n4r8cROsjdZGHLx7LveFVcGmk31F/pYOIUlxcj0OIWkPjBRLHl5liHxkBJ4vbCWx9iDO4MNX2S2Jnb+Z9pN+IW7rcV1+LSJSiqnASLHmVcaXyEdH4fb8Vqrd/hS+nq7sTT1P3Bfvsu+Nlmz9ZYGKjIhIKaQCI07Bu2w5et9Wn1+G3MqAW2vR23URoXl7afTz4+wYewt7Nq4wO6JIsVQCzhKQEuB6/B2qwIhT8fV05fmYelTq+wOJlbuQY7jQIGczdb67l43j7+Xovm1mRxQpFv68zXtm5mUexChyA/35d/i/j0K4FjqJV5xa8uG9HPn6ZZqf/RGrxSDXsLGkan9adn2JSmV16bWUbidOnODs2bNUrlwZLy+va7oDrcjVMAyDzMxMUlNTKVeuXMFDH/90Ld/fKjBSIvy2bQ3nFr1Ek6x1dMl+hW2uDel7a216tamBh6vtnycQKYEMwyA5OZmzZ8+aHUVKuXLlyhEQEHBRiVaBUYGR/5e07hdGrLGy5egfl173L7OMWxrXplnHp7DaVGSkdMrPz7/kk41FbgRXV1dsl/nnrwqMCoz8F7vd4LvNx/nwh9XMze6HlyWbPS51yL/jdepHxJgdT0RE/p/uAyPyX6xWC52aVmH+oA5srvkk5w1P6uTtof4PD7NhQieSD+81O6KIiFwjFRgpsTy8yhDV4w2y+qxjTYV7sRsWmp9bRrkPokj44Hkyz6eZHVFERK6SCoyUeH4BwUQ8+wkHOi9mu1sjPCy5NDs8m+6TF/Hd5uO6T4aIiBPSOTBSqhh2O5t+nMPKDZuZeP4OACJqVGB0dCVq16ptcjoRkdJF58CIXCGL1UqzDj15eugEYu+og4erlayDa6n2cUsS33uKtDOnzI4oIiJXQAVGSiUPVxvP3h7KT7HteDpgD26WfCJT55P3TjPWLZyi5yuJiBRzKjBSqlUt70XHgVPZeutHHLZWoSJptEh6iZ1jb+bgzvVmxxMRkctQgREBGrV7gIChG0mo0Z9Mw52bcrZSZe6dLH9/KJk5eWbHExGR/6ECI/L/3Nw9iOoxmrQnfmWTVytcLfn8cCCf6AkrWLojxex4IiLyX3QVkshlbFj+HYMSvThyNhuA3rVO07NDa/yr1jI5mYhIyaCrkESug+a33MuS2Fvpe0stfKzZPHb0VcrMasWaL0aTn6fDSiIiZlKBEfkbnm42hrSvxzc963POtTLeliwido/jt7GR7N+aaHY8EZFSSwVG5ArUDq1PneGrWHPTy6TjRWjeXqp/1YGEmc+SdSHD7HgiIqWOCozIFbLabEQ8PJicpxPZWOZmXCx2oo7P4eS4cNbv/M3seCIipcpVFZipU6cSEhKCh4cHERERrF279rJjZ82aRdu2bSlfvjzly5cnOjr6ovGPP/44Foul0Kt9+/ZXE03kuvMLqk6zF/7DplZTOUl5NucF8+Ccnby4YCvpWblmxxMRKRUcLjDz5s0jNjaWkSNHsnHjRpo0aUJMTAypqamXHL98+XK6devGsmXLSEhIIDg4mDvvvJNjx44VGte+fXtOnDhR8Priiy+ubotEbpCmdz6K28D1bGr0MgCfrzlM14nfsXnZVyYnExEp+Ry+jDoiIoIWLVowZcoUAOx2O8HBwQwYMIBhw4b94/r5+fmUL1+eKVOm0L17d+CPPTBnz55l4cKFjm8BuoxazJew/zTDvtnC4PQx3G1bwzrf9tTpMQXfCpXMjiYiUmzdsMuoc3Jy2LBhA9HR0X9NYLUSHR1NQkLCFc2RmZlJbm4uFSpUKLR8+fLlVK5cmbp169KnTx9Onz592Tmys7NJT08v9BIxU1StisQNaI1fUA3shoUWaXHkTg5n05JPzY4mIlIiOVRgTp06RX5+Pv7+/oWW+/v7k5ycfEVzDB06lKCgoEIlqH379nz88cfEx8fz5ptvsmLFCjp06EB+fv4l5xgzZgy+vr4Fr+DgYEc2Q+S68PRwI7LPDPZ0/JLD1ir4cZamq/uxfmJn0k7rTr4iIkXphl6FNHbsWObOncuCBQvw8PAoWN61a1fuvfdeGjVqRKdOnVi0aBHr1q1j+fLll5xn+PDhpKWlFbyOHDlyg7ZA5J/Va3kHlQevIyGwO/mGhfD0n8h5N4K1CcvMjiYiUmI4VGD8/Pyw2WykpBT+t8mUlBQCAgL+dt3x48czduxYlixZQuPGjf92bM2aNfHz82Pfvn2X/L27uzs+Pj6FXiLFiYenN1HPvMvee77msLUK2YaNJ749xeAvN+tKJRGRIuBQgXFzc6N58+bEx8cXLLPb7cTHxxMVFXXZ9caNG8drr71GXFwc4eHh//g+R48e5fTp0wQGBjoST6TYqRd+O5VfWMt/Gr5LhsWLLzccJWbiCpISfzY7moiIU3P4EFJsbCyzZs1izpw57Ny5kz59+pCRkUHPnj0B6N69O8OHDy8Y/+abb/LKK6/w4YcfEhISQnJyMsnJyZw/fx6A8+fPM3jwYBITEzl48CDx8fHcd9991K5dm5iYmCLaTBHzeHiVoc9DdzH/mSiqV/TilozFhMXdz5qpvbiQcc7seCIiTsnhAtOlSxfGjx/PiBEjCAsLIykpibi4uIITew8fPsyJEycKxk+bNo2cnBwefPBBAgMDC17jx48HwGazsWXLFu69917q1KlDr169aN68Ob/++ivu7u5FtJki5msRUoEfBrblrqBMACJOfsXJCRHs2bjC5GQiIs7H4fvAFEe6D4w4my3LvyZg+QtU5nfyDCvrqj1B+GNv4Oqm0i4ipccNuw+MiBSNxrd0xn1AIuvL3v7HM5WOvM+BcW04sm+r2dFERJyCCoyISXwr+hP+/DdsaDGBdLyplbuXIXN+5ou1hykBO0ZFRK4rFRgRkzXv+CQXnvyVmRWeJyG3NsO/2cpTH2/gdHqG2dFERIotFRiRYsC/ai2eGfAyL91VHzeblQO7NpI1MYwtejCkiMglqcCIFBNWq4Wnbq7Jwn6teanMIqqQSuMVvUh872myszLNjiciUqyowIgUMzcF+dDquc9Z49cZgMjUeRx9qxWHdm4wOZmISPGhAiNSDHl4lSGi/4cktZ3BGXyolX8A/7kxrJn/FobdbnY8ERHTqcCIFGNht3clv/cqtng0x8OSS8SO1/lo+lukZep5SiJSuqnAiBRzfgHVaDh4KYmhsSTYb2L04frcNflX1h/83exoIiKmUYERcQJWm43IR0ZS5qnFBFcsy7GzF3hs5q8s++xN8vPyzI4nInLDqcCIOJFGweVZ9Gxb7m9ahVjrXG7d+wY737qdU8cPmR1NROSGUoERcTJl3F14u0sYYeFtyDTcaZidhGVmW7b+ssDsaCIiN4wKjIiTatGpH6ceWcJv1hAqkkaD+J4kzBpEXm6O2dFERK47FRgRJ1atThhBL6xiTcVOWC0GUcc+Yu+4W0g5ut/saCIi15UKjIiT8/AqQ8SAOWxoMYHzhidVcn7jyQ9WsXx3qtnRRESuGxUYkRKieccnOfvYT7xVdhhbL1Sk5+x1jP9xN3l5+WZHExEpciowIiVI1doNeWngAB6NrIZhwKYVC9g57jZOnThsdjQRkSKlAiNSwni42ni9UyPe7dKIMa4f0ignCWa0Zduq/5gdTUSkyKjAiJRQ9zStBo/M54C1On6cpf6Sx0iY8yL2fB1SEhHnpwIjUoJVqxNGwPOrWFeuAzaLQdSBqWwZfxdpv580O5qIyDVRgREp4Ty9y9Ji0FzWNhpFtuFK2IVEMt5tzc79B8yOJiJy1VRgREqJlp0HceSBhRy3+JOQV4f7PtrJvHU6uVdEnJMKjEgpUrtJG7wHrCK+5hBy8gyGfr2VkfMTyLqQYXY0ERGHqMCIlDK+FSoxpUdbBsfUxWax03brcI6Mv5kTh3abHU1E5IqpwIiUQlarhX631mbeQ0E0t+0jNH8fnh/dxtYV35gdTUTkiqjAiJRi4c2ak/3EMva6hFKO8zT4+QkSZg/TpdYiUuypwIiUcgHVQgl+fgVrK9zzxwMhD05j84R7OJf2u9nRREQuSwVGRPDw9Kbls5+yttGr5BguNM1cxa7J97Mv9bzZ0URELkkFRkQKtOz8HAfu+ZIjBPLvzAfpNHUVP25PNjuWiMhFVGBEpJC64bfh8dwGPEPCOZ+dxzOfbOCTb77TeTEiUqyowIjIRSr5evPZkxE80boGjSy/8fDmJ9gy/i7Sz542O5qICKACIyKX4WqzMuKem3iptRcAYRcSOTu5LYd2J5kbTEQEFRgR+QeRd/ficKdvSKEi1ezHqPh5e5J++sLsWCJSyqnAiMg/Cm16M7beK9jh2pAylguErexNwkdDMex2s6OJSCmlAiMiV8QvIJjaL8Szxu8BAKIOTeejmePJzMkzOZmIlEYqMCJyxdzcPYjo/xFrG44kzt6S1w7W54H3VnPk90yzo4lIKaMCIyIOa/lgLBV7zqViGQ92JZ/joSnL2LYm3uxYIlKKqMCIyFVpUaMi3/VvQ8OgsgzOnUbdxQ+xZv5bZscSkVJCBUZErlpQOU++fKol1cu54WrJJ2LH66yZ0pPcnGyzo4lICacCIyLXxNPTg+bPfUVCjf7YDQsRp75h94Q7STudYnY0ESnBVGBE5JpZrFaieoxmS5v3yDA8aJidxLkpN3No10azo4lICaUCIyJFJuyOf5H68Hcct1SmqpGM7YuHWb7zmNmxRKQEUoERkSJVo0EEHn2Ws8UtjMG5T/PEx0l8sPIAhmGYHU1EShAVGBEpchUqV6He4J+p3rw9dgNeW7SDdz5fSE52ltnRRKSEUIERkevCzdXGmAca8XLH+tS2HuOJPX3ZO+EOzp5KNjuaiJQAV1Vgpk6dSkhICB4eHkRERLB27drLjp01axZt27alfPnylC9fnujo6IvGG4bBiBEjCAwMxNPTk+joaPbu3Xs10USkGLFYLDzZtiZv3emHFYMGOVs4P7WdnmgtItfM4QIzb948YmNjGTlyJBs3bqRJkybExMSQmpp6yfHLly+nW7duLFu2jISEBIKDg7nzzjs5duyvE/vGjRvH5MmTmT59OmvWrMHb25uYmBiysrS7WaQkaHrLA5zquojjFn+qGsmU/6IDW3/51uxYIuLELIaDZ9ZFRETQokULpkyZAoDdbic4OJgBAwYwbNiwf1w/Pz+f8uXLM2XKFLp3745hGAQFBfH888/zwgsvAJCWloa/vz+zZ8+ma9eu/zhneno6vr6+pKWl4ePj48jmiMgN9HvqMVJmPUj93B3kGVY2NHiRiIcHmx1LRExyLd/fDu2BycnJYcOGDURHR/81gdVKdHQ0CQkJVzRHZmYmubm5VKhQAYADBw6QnJxcaE5fX18iIiIuO2d2djbp6emFXiJS/FWoXIWaz8ezzvdOXCx2Ina8ztdz3iHfriuURMQxDhWYU6dOkZ+fj7+/f6Hl/v7+JCdf2Yl5Q4cOJSgoqKCw/LmeI3OOGTMGX1/fgldwcLAjmyEiJnL38CJ84DwSQvqSaK/PsJ0hPPXxes5n55kdTUScyA29Cmns2LHMnTuXBQsW4OHhcdXzDB8+nLS0tILXkSNHijCliFxvFquVqMfH8PsDX2J1cePnXak8/N6vJB87ZHY0EXESDhUYPz8/bDYbKSmFn3GSkpJCQEDA3647fvx4xo4dy5IlS2jcuHHB8j/Xc2ROd3d3fHx8Cr1ExPncFRbM/Gei8CvjzgOnZ2Cb1Y69Sb+aHUtEnIBDBcbNzY3mzZsTHx9fsMxutxMfH09UVNRl1xs3bhyvvfYacXFxhIeHF/pdjRo1CAgIKDRneno6a9as+ds5RaRkaBJcjm+fasJtbjupxBmqLniAjT9+YnYsESnmHD6EFBsby6xZs5gzZw47d+6kT58+ZGRk0LNnTwC6d+/O8OHDC8a/+eabvPLKK3z44YeEhISQnJxMcnIy58+fB/64T8SgQYN4/fXX+e6779i6dSvdu3cnKCiITp06Fc1WikixVsW/EpUGLmOLRws8LTmErR5A4icjMOx2s6OJSDHl4ugKXbp04eTJk4wYMYLk5GTCwsKIi4srOAn38OHDWK1/9aJp06aRk5PDgw8+WGiekSNH8uqrrwIwZMgQMjIyePrppzl79ixt2rQhLi7ums6TERHnUta3Ajc9v5g1M54h4tQ3RO5/hzVTDtC8zwe4uLqZHU9EihmH7wNTHOk+MCIlh2G3s2buG7TcPR6rxWC9Z2vqPruQsp4qMSIlzQ27D4yIyPVmsVqJ/NfLbG49hQzDg/fTWvDQjESOn71gdjQRKUZUYESkWGp656McfHQ1G7zbsiv5HJ2mrmLbkdNmxxKRYkIFRkSKrQahtVjQtxV1/MtgO3ecMu+3ZvPPc82OJSLFgAqMiBRrVct78VWfVozy+4kQywkarujNmvnjzI4lIiZTgRGRYs/Hw5Vbn53F2nJ3YbMYROwYTeL0vtjz882OJiImUYEREafg6uZOi2c/IyGkDwCRyZ+R9Pb9ZGWeNzmZiJhBBUZEnMYfz1Aay/pmb5JjuNDs/AoOvn0HZ05d2cNkRaTkUIEREacTfm9v9sZ8Qhre5Gdn8MhHGzl0OsPsWCJyA6nAiIhTatDqLs52XcSLnq+w47TBA++tZtPhM2bHEpEbRAVGRJxW9XrNmNX/HhpW8eF0Rg5x77/CpiWfmh1LRG4AFRgRcWqVy3ow7+ko+lY7ynDrJzRZ1Z81c8eYHUtErjMVGBFxet7uLsQ+1ZM1Fe7FajGI2DWWhBn9dJm1SAmmAiMiJYKLqxst+88hoUY/AKJOfMrGSQ+RnZVpcjIRuR5UYESkxLBYrUT1eIN1TceQa9gIPxfPvokxpJ05ZXY0ESliKjAiUuK0uK8vu27/iPOGJw1ytjB95hQ9zVqkhFGBEZESqdHN95Hy4EKm2B5l2pkW3P/eKnYlp5sdS0SKiAqMiJRYtRpFcv+A8YRWLkNKeja9pi1la0Kc2bFEpAiowIhIiValnCdf9W5F6+revG28Sd24R1j//SyzY4nINVKBEZESz9fLlQ8eb4GLjz9uljzC171A4mejzI4lItdABUZESgUPT2+aDFrAmkoPAhC5dyKJ057RvWJEnJQKjIiUGjYXF1r2mUVirYEARKbMZdOkB3WvGBEnpAIjIqWKxWol8rF/s77Zm+QaNpqf+5nEyd05l5VrdjQRcYAKjIiUSuH39mbnbR9w1KjEqDMxPDwjkdT0LLNjicgVUoERkVKrcbv7SeuVQLp3CDtPpHP/e6v57cRJs2OJyBVQgRGRUq1BtUp806c1IRW9CElfS9kZLdi9/mezY4nIP1CBEZFSr1pFL77uHcUQ78VU4gzV/tOFzT/PNTuWiPwNFRgREaBiWQ9qP/sdWzxa4GnJocGKPqxdMNnsWCJyGSowIiL/z7tsOerHfs863xhcLHZabn6FhDkvYtjtZkcTkf+hAiMi8l9c3dwJHziXhMBHAYg6MJW1057SDe9EihkVGBGR/2GxWol6ZiqJdV4A4GhyMgPnbiI7TyVGpLhwMTuAiEhxFfmvV1i1tD4vL3fjwtYUzlxYz/THmlPGXf/oFDGb9sCIiPyN1nc8yIweUXi52Vi9L5XFk/pwOuWo2bFESj0VGBGRf3BznUrMfTqSEZ5f8fCF+VyYHs2x33aaHUukVFOBERG5Ao2rluP2fz3PcUtlqhoncP+4Pfu3JpodS6TUUoEREblCwaGNcXvqJ/bbauDHWSp/1YntCT+YHUukVFKBERFxgF9Qdfz6/8QOt0aUtVygdtxjbPzxU7NjiZQ6KjAiIg7yLe9HzUFxbPJqjbsllzqrn+eblZvNjiVSqqjAiIhcBQ+vMjR6biGJ5e9jYG4/YhcdZcrPezEMw+xoIqWCCoyIyFVycXUj4tk51Gv3MADjl+xh0jfLdNdekRtABUZE5BpYLBYGx9TjlbtvItiSwiNberLpnYfIyc4yO5pIiaYCIyJSBHq1qcH4NjbKc47m6fHsersjmefTzI4lUmKpwIiIFJGIjj3YcctMMg13Gmet58ikOzh7KtnsWCIlkgqMiEgRanLrgxy++3PS8KZu3m7S3osm5eh+s2OJlDgqMCIiRaxei2jOdPmOVCpQ3X4E4/0YDu3bbnYskRJFBUZE5DoIqR9Ofs84jliCSLWX4ZHP97P1qM6JESkqV1Vgpk6dSkhICB4eHkRERLB27drLjt2+fTudO3cmJCQEi8XCpEmTLhrz6quvYrFYCr3q1at3NdFERIqNwOp18e7zE29WHM3RTBtdZyawet8ps2OJlAgOF5h58+YRGxvLyJEj2bhxI02aNCEmJobU1NRLjs/MzKRmzZqMHTuWgICAy87boEEDTpw4UfBauXKlo9FERIqdCpWrMKN3DK1qVSQjJ59f5oxg049zzI4l4vQcLjATJ07kqaeeomfPntx0001Mnz4dLy8vPvzww0uOb9GiBW+99RZdu3bF3d39svO6uLgQEBBQ8PLz83M0mohIsVTG3YWPerbguZpHGWb7jMarB7L2q4lmxxJxag4VmJycHDZs2EB0dPRfE1itREdHk5CQcE1B9u7dS1BQEDVr1uSRRx7h8OHDlx2bnZ1Nenp6oZeISHHm7mKj/xO9WFvhHmwWg5bbRpEw5yUMu93saCJOyaECc+rUKfLz8/H39y+03N/fn+Tkq7/XQUREBLNnzyYuLo5p06Zx4MAB2rZty7lz5y45fsyYMfj6+ha8goODr/q9RURuFJuLCy36f0xClccBiDowhTUz+urRAyJXoVhchdShQwceeughGjduTExMDIsXL+bs2bPMnz//kuOHDx9OWlpawevIkSM3OLGIyNWxWK1EPfUOiaGxAESmfMGGyf8iLzfH5GQizsWhAuPn54fNZiMlJaXQ8pSUlL89QddR5cqVo06dOuzbt++Sv3d3d8fHx6fQS0TEmUQ+MpJ1YW+QZ1hpkRbH1PdnkpWrPTEiV8qhAuPm5kbz5s2Jj48vWGa324mPjycqKqrIQp0/f579+/cTGBhYZHOKiBQ3LTr1Y2vb9xif3423D9Wgx4drSc/KNTuWiFNw+BBSbGwss2bNYs6cOezcuZM+ffqQkZFBz549AejevTvDhw8vGJ+Tk0NSUhJJSUnk5ORw7NgxkpKSCu1deeGFF1ixYgUHDx5k9erV3H///dhsNrp161YEmygiUnw1je5Gm56jKevuwpoDv/PU9KWcTjlqdiyRYs/F0RW6dOnCyZMnGTFiBMnJyYSFhREXF1dwYu/hw4exWv/qRcePH6dp06YFP48fP57x48fTrl07li9fDsDRo0fp1q0bp0+fplKlSrRp04bExEQqVap0jZsnIlL8RdasyBdPR9L7gxUM+30UmdMzyXl8IYHV65odTaTYshiGYZgd4lqlp6fj6+tLWlqazocREad1eP8OXD+5l0BOkkoFLnT5iur1m5sdS+S6uZbv72JxFZKIiEC1WjdhfWoJB63BVOZ3fOfdy56Ny82OJVIsqcCIiBQj/lVqUq7vT+xxqUM5zlP124fZ+su3ZscSKXZUYEREiplyfgFUGbiUbe5N8bJkUzf+CTYs+dzsWCLFigqMiEgx5F22HKHPLWaj982k4cXg5ZnMW3f5R6yIlDYqMCIixZS7hxdNnlvAh/Vm8Zs9kKFfb2XGiv1mxxIpFlRgRESKMZuLC0O6xtC7XS0AVv84j9WzBukhkFLqOXwfGBERubEsFgvDOtQjyHaWB1dNwutYNmumnCa870fYXPSPcSmdtAdGRMRJdL8zkm2Nh5NvWIj4/Ts2T3qAnOwss2OJmEIFRkTEibTs/BxbIt8mx3Ch2fkV7Hr7LjLPp5kdS+SGU4EREXEyTTv0ZPftH5JpuNM4awOHJ8WQdjrF7FgiN5QKjIiIE2p0830cuWcuaXhTL28n380cSWq6DidJ6aECIyLipOqG38bvD3/LN9Y7eTXtLjpPX82h0xlmxxK5IVRgREScWI2bWtCi/xyqVizLkd8v0GXaSvbv3mJ2LJHrTgVGRMTJBVfw4sveUdTzL0Ns1lQqfdGeXWuXmh1L5LpSgRERKQEql/VgXs8mNPY8iQ8ZVPv+EbYs/9rsWCLXjQqMiEgJ4VuuHNUGxrHFIxwvSzb1lj3FhsUfmR1L5LpQgRERKUG8yvhS77nv2VD2Vtws+TRd8xxrv5podiyRIqcCIyJSwri5exA28CvWVLwPq8Wg5bZRrP7kVbNjiRQpFRgRkRLI5uJCy36zSQjqQY5hY8ZON8b8sBPDMMyOJlIkVGBEREooi9VK1NOTWRD5JSvsTZix4jeGf7OVfLtKjDg/FRgRkRKuS4fbebNzI6wWSFy/ltWTHiM7K9PsWCLXRM9hFxEpBbq0qIavm4U63zxPzfQTbH27IzX7L8C7bDmzo4lcFe2BEREpJdo3CSbj9jfJNNxplL2Ro+/cqYdAitNSgRERKUX++yGQdfN28/vUOzh1/KDZsUQcpgIjIlLK/PkQyJOUp4b9EDmz7uTYb9vNjiXiEBUYEZFSqMZNLcjtEcdRSwBBRgp7PxnIruR0s2OJXDEVGBGRUiqoRj08nlnKcte2DLrwJA9PT2DDoTNmxxK5IiowIiKlmF9ANZo+t4Ba1auRnpXHo++vYc3GDWbHEvlHKjAiIqWcr5crn/RqSbs6lehsj6Ppt3ewYfH7ZscS+VsqMCIigpebC7Mea859FQ7//0MgX2Dtl+PNjiVyWSowIiICgJurjWYD57OmYqc/HgK5/TUS5ryEYbebHU3kIiowIiJS4I+HQH5EQpWeAEQdmMKamf1VYqTYUYEREZFCLFYrUU9NIjE0FoDI5M9YN/kR8vLyTU4m8hcVGBERuaTIR0aytsnr5BsWlp30od8Xm8jKVYmR4kEFRkRELqvl/QNIjPmODyz38eP2FJ6YvY7z2XlmxxJRgRERkb/XutXNzO7ZkjLuLmzef5TVE7tx5uQJs2NJKacCIyIi/6hVLT++eCqStz3f587sJaRPiyb5yD6zY0kppgIjIiJXpFFVX+p2HUsyflS3H4UPYji8J8nsWFJKqcCIiMgVq143DJ6I47C1CgGcouzn97A36VezY0kppAIjIiIOCagWSpneS9nrEkp50glc8BDbVv3H7FhSyqjAiIiIwypUrkLgs0vZ7taEMpYL+C55jiVbD5sdS0oRFRgREbkqZXzKU+u5H1hVJoancp6n9+dbmb/+iNmxpJRQgRERkavm4elNxHNzady8FXYDhny1hc/jVpgdS0oBFRgREbkmLjYrb3ZuzDM316SlZScPJHQmYYaenyTX11UVmKlTpxISEoKHhwcRERGsXbv2smO3b99O586dCQkJwWKxMGnSpGueU0REiheLxcLwu+ozsMEFPCy5RJ34hHXvPkpebo7Z0aSEcrjAzJs3j9jYWEaOHMnGjRtp0qQJMTExpKamXnJ8ZmYmNWvWZOzYsQQEBBTJnCIiUjy1fvQV1jYaRb5hoeWZ79k66X6yLmSYHUtKIIthGIYjK0RERNCiRQumTJkCgN1uJzg4mAEDBjBs2LC/XTckJIRBgwYxaNCgIpsTID09HV9fX9LS0vDx8XFkc0RE5DrY+OMnNFw9CDdLHtvdmlCt30LK+lYwO5YUM9fy/e3QHpicnBw2bNhAdHT0XxNYrURHR5OQkODQG1/LnNnZ2aSnpxd6iYhI8dEs5jH23Dmb84YnDXI2kzI5mtOnUsyOJSWIQwXm1KlT5Ofn4+/vX2i5v78/ycnJVxXgauYcM2YMvr6+Ba/g4OCrem8REbl+Gra+h+QHvuJ3fNiV48fDH23j6JlMs2NJCeGUVyENHz6ctLS0gteRI7rvgIhIcVS7SRvOP7aE8V7Ps/90Fp2nrWZPyjmzY0kJ4FCB8fPzw2azkZJSeDdgSkrKZU/QvR5zuru74+PjU+glIiLFU7Va9Znbtx11/MuQmn6BjdN6sWvtUrNjiZNzqMC4ubnRvHlz4uPjC5bZ7Xbi4+OJioq6qgDXY04RESleAnw9mP9MFC/5raQrP1L9+3+x+ef5ZscSJ+bwIaTY2FhmzZrFnDlz2LlzJ3369CEjI4OePXsC0L17d4YPH14wPicnh6SkJJKSksjJyeHYsWMkJSWxb9++K55TREScXzkvNx7p8zKbPVviacmhwYpnWPfte2bHEifl4ugKXbp04eTJk4wYMYLk5GTCwsKIi4srOAn38OHDWK1/9aLjx4/TtGnTgp/Hjx/P+PHjadeuHcuXL7+iOUVEpGTw9C7LTc8tYv2URwhPX0qLTcNJPH+KyEdGmB1NnIzD94EpjnQfGBER52LPz2ftjD5Eps4DICGoO5FPvoPF6pTXlshVumH3gRERESkKVpuNiN7TSajRD4Bmxz5n4uffkZev5yfJlVGBERERU1isVqJ6vMHaRqOIzevLu9tc6f3pRrJy882OJk5ABUZEREzVsvMg7v1XP9xcrPy0M4UhMxaQduaU2bGkmFOBERER08U0CODjJ1pS2/0sw04O5vd3b+fk8YNmx5JiTAVGRESKhciaFZnxUG3cLHZq2A+SO+sOjuzdbHYsKaZUYEREpNio1TCC7O5xHLUEEmSkUuazu9mzcYXZsaQYUoEREZFipUrN+nj2/om9ttqUJ52q3z7E1l8WmB1LihkVGBERKXYq+lclcOBPbHVvipclm7rxvVj94zyzY0kxogIjIiLFUhmf8tSJ/YH1ZW9jvxFE72UWPlh5wOxYUkyowIiISLHl7u5Js0Ff8V3YLNLx5rVFOxj7wy4Mu254V9qpwIiISLFmtdkYcn8kg2PqApCzcgrr3ulGbk62ycnETA4/zFFERORGs1gs9Lu1NtVtp7gz/nPc0vLZPLEjof2/xquMr9nxxATaAyMiIk7j7psj2XHze1ww3GiStY4jk+7gzMkTZscSE6jAiIiIUwm7vSuH7v6Cs5Shbt5uzr93G8cP7DI7ltxgKjAiIuJ06rWIJq3b9yRTiWDjOG5zYti/ZZXZseQGUoERERGnVL1uGLanf+I3awh+nGXu11+xap8eAllaqMCIiIjTqhQUgt+zPzPTdyCzsqN5/KO1fLf5uNmx5AZQgREREafmU64iPQaMpGOjQHLzDV7+YiUr5r1tdiy5znQZtYiIOD13FxvvdmuKfxkbMRueJmLnLhKn7aLl0+9htdnMjifXgfbAiIhIiWC1WnjlnkYYoTEARKbMZdOkzmRnZZqcTK4HFRgRESkxLFYrkY/9m/XN3iTHsNH83DL2TYwh/exps6NJEVOBERGREif83t7svv0jzhueNMjZwunJt5F6TA+CLElUYEREpERqdPN9JHf+hlOUo4b9IMff78ae5HSzY0kRUYEREZESq3bjVuQ8voRt1nq8kNWTztMTSPxNh5NKAhUYEREp0YJC6lL1+V8oV60R57Ly6P7BWn5K3Gh2LLlGKjAiIlLilfN259MnI2jfIIDmxlba/nAniZ+NMjuWXAMVGBERKRU8XG1MfaQZfYOP4G7JJXLvRBLfe4r8vDyzo8lVUIEREZFSw2a10OaZd0is/RwAkanz2fL2fWRlnjc5mThKBUZEREoVi9VK5KOvsqHFBHIMF5pmrOTgxGjOnDxhdjRxgAqMiIiUSs07Psm+9p+Sjjf18nZy/r3bOHL0sNmx5AqpwIiISKl1U1QHznRdxAkqsT6vBp0+3MWmw2fMjiVXQAVGRERKter1muHyzM/MqfQCpzNz6TYrkbitOpxU3KnAiIhIqVcpsBqfPXMzt9atRE5uHvnze5D4+etmx5K/oQIjIiICeLu7MKt7OK/XPUBH2xoi97yly6yLMRUYERGR/+dis9KtR38Saz4L/P9l1hPv5ULGOZOTyf9SgREREfkvFquVyO6vsb7F+D8us85cxZG3b+NU8hGzo8l/UYERERG5hPCOT7Gvw2ecpQx18vaQM+M2Du3SM5SKCxUYERGRy7gpsj3nHlnMUUsAFe1neGnuKhL262nWxYEKjIiIyN8IDm2CV59ljC0/ipVZNen+4Rq+2XjU7FilngqMiIjIP6hQOYih/XrTsXEgufkGH3z5Las+HIJht5sdrdRyMTuAiIiIM/BwtfFu16bU9jXourY/gYd/Z/2k/TTq+wnuHl5mxyt1tAdGRETkClmtFp7r2JwjjQeSa9gIT/+J3yboQZBmUIERERFxUMvOg9gd/RHpeFE/dzsZ793Ckb2bzY5VqqjAiIiIXIWGbe/jTJdFHLdUpqqRjM9nHdi+6nuzY5UaKjAiIiJXqXr95rg98zO7XerhSwbHf5zI/PW64d2NcFUFZurUqYSEhODh4UFERARr16792/Fffvkl9erVw8PDg0aNGrF48eJCv3/88cexWCyFXu3bt7+aaCIiIjeUX0Aw1WPjWVL+XzyX05shX23hzbhd2O2G2dFKNIcLzLx584iNjWXkyJFs3LiRJk2aEBMTQ2pq6iXHr169mm7dutGrVy82bdpEp06d6NSpE9u2bSs0rn379pw4caLg9cUXX1zdFomIiNxgHl5liB7wHk/c1hiAacv38cW0f+sZSteRxTAMhypiREQELVq0YMqUKQDY7XaCg4MZMGAAw4YNu2h8ly5dyMjIYNGiRQXLIiMjCQsLY/r06cAfe2DOnj3LwoULr2oj0tPT8fX1JS0tDR8fn6uaQ0REpCh8s/EoexaMYZjtU/a6hFLuia+oFBRidqxi6Vq+vx3aA5OTk8OGDRuIjo7+awKrlejoaBISEi65TkJCQqHxADExMReNX758OZUrV6Zu3br06dOH06cvf6vm7Oxs0tPTC71ERESKgweaVeXejvdwhrKE5u3FmHkr+zavNDtWieNQgTl16hT5+fn4+/sXWu7v709ycvIl10lOTv7H8e3bt+fjjz8mPj6eN998kxUrVtChQwfy8/MvOeeYMWPw9fUteAUHBzuyGSIiItfVTZHtyey+lEPWYCrzO0HfPMDGHz8xO1aJUiyuQuratSv33nsvjRo1olOnTixatIh169axfPnyS44fPnw4aWlpBa8jR3TGt4iIFC9Vatan/LMr2OLRHC9LNs0S+pMw5yU9fqCIOFRg/Pz8sNlspKSkFFqekpJCQEDAJdcJCAhwaDxAzZo18fPzY9++fZf8vbu7Oz4+PoVeIiIixY1PuYrc9Hwca/w6A9Dit/cY/8kCsvMufYRBrpxDBcbNzY3mzZsTHx9fsMxutxMfH09UVNQl14mKiio0HmDp0qWXHQ9w9OhRTp8+TWBgoCPxREREih0XVzci+n/Imvov8lp+D6bu9KDbzEROnss2O5pTc/gQUmxsLLNmzWLOnDns3LmTPn36kJGRQc+ePQHo3r07w4cPLxg/cOBA4uLimDBhArt27eLVV19l/fr19O/fH4Dz588zePBgEhMTOXjwIPHx8dx3333Url2bmJiYItpMERERc0V0GUp0j5fw8XBh4+GzDHh3Pvu3Jpody2k5XGC6dOnC+PHjGTFiBGFhYSQlJREXF1dwou7hw4c5ceKvh1q1atWKzz//nJkzZ9KkSRO++uorFi5cSMOGDQGw2Wxs2bKFe++9lzp16tCrVy+aN2/Or7/+iru7exFtpoiIiPnahlZiYb/WNK5oMDprNIFf3cumJZ+aHcspOXwfmOJI94ERERFnkvb7SQ7PeIhG2ZsASAzpS0T30VisxeLamhvmht0HRkRERK6db4VK1H/+R9b4PQBA5MH32Ph2Z9251wEqMCIiIiZwcXMnov9HrGnwCrmGjebnfuboxFtIObrf7GhOQQVGRETERBEPvcCemE//uHNv/j62fdCXjYfPmB2r2FOBERERMVmDVndxocdPrHEJZ8iFHnSdmcjXG46aHatYU4EREREpBoJq1KPB4B9pdlMdcvLsPP/lZuZ9OoO83ByzoxVLKjAiIiLFRBl3F2Y82pxnb6vNA9Zf6LJvCDvHx5B2OuWfVy5lVGBERESKEavVQuyddXn05vpkGu40yt7IuSk3c3DnerOjFSsqMCIiIsVQs/aPc+LB/3DcUpmqRjKV5nbUTe/+iwqMiIhIMVWrUQSefX9hu1sTvC1ZNF3dj4QPnseer4dBqsCIiIgUY+UrBVLnhaUkVn4YgKgj7zPu/Tmcy8o1OZm5VGBERESKOVc3dyL7zmJd2Ggm5T/E9AP+dJq6iv0nz5sdzTQqMCIiIk6iRaf+3Pr0eAJ9Pdh/MoNnpnzLpp/mmh3LFCowIiIiTqRJcDm+69+G1tXLMN4YT9OVz5D4fmypOy9GBUZERMTJVCrrzuxekWT5Nwcg8ugHbB3fnrTfT5qc7MZRgREREXFCrm7uRPZ7n3VNx5BluNLkwlrOvduGA9vXmB3thlCBERERcWIt7uvL0Qe+LbhfjP/8e1i/aKbZsa47FRgREREnV7tJa7z6/cpW92Z4WbLxXDuF177bTG6+3exo140KjIiISAlQzi+AmwYvZVXVp+mdO4gPVh+l28xEUtOzzI52XajAiIiIlBA2FxdaP/kWIx7rSFl3F9YfOsPHk4azIzHO7GhFTgVGRESkhLnjJn++G9CGhyseIDb/I+r80I3Ez0Zh2EvOISUVGBERkRKohp83r/Z5jI2+t+NisRO5dyKbJtzHubTfzY5WJFRgRERESiivMr40H/Qla+q/SI5ho1nGL5yd1LpEXGqtAiMiIlKCWaxWIroM5cA9X5GMH8HGcQLm383a76abHe2aqMCIiIiUAnXDb8O930q2eDTH05LDl2v2M/ybLWTlOucjCFRgRERESonylQJp8MIS/lP/Lb6y38IXa4/wwHurOXjynNnRHKYCIyIiUorYXFy4p8vTfPxESyp6u3H8xDFypkSxMW622dEcogIjIiJSCrUNrcT3z7ZlRMWfqWM5QrPEgSS+9xQ52c5x4zsVGBERkVIqwNeDewa+S0LgowBEps7n0FttOPbbTpOT/TMVGBERkVLM1c2dqGemsrnNdNLwJjRvLz4f38qmYn5ISQVGREREaBLdjQtPrGCXS33KcoGmiQNZ+OFYsvOK51VKKjAiIiICQEC1UGoNWUFC4KMctldixJ6aPPDeag6cyjA72kVUYERERKTAn4eUDj60FJtXObYfT+fuyb/w69KFZkcrRAVGRERELnJzoxosHtiWljUqcF/+Utqu6sG6SV3JPJ9mdjRABUZEREQuI9DXky+eiuSu2h7kGxZanP2BUxOj2L810exoKjAiIiJyeTarhTY932BXzOekUoFq9mNU/epuEr94A8NuNy2XCoyIiIj8owat7sK132qSPCNxt+QSuftNNr/Vgd9Tj5mSRwVGRERErkj5SoE0GfwDa+oNI9twpVHmGobO+IaVe0/d8CwqMCIiInLFLFYrEV2Hc+yh75ns2YelGTV57MM1jFm8k5wb+GRrFRgRERFxWM2GEfSOfZ1/RVTDMOCnX39l75ttOLQ76Ya8v8UwDOOGvNN1lJ6ejq+vL2lpafj4+JgdR0REpFT5cXsyPl8+RBRbuGC4saXhUFp2jsVi/fv9JNfy/a09MCIiInJNYhoEUPvpT9jq3gxPSw4R218jaXzH63qCrwqMiIiIXLNKQSE0GPITiaHPk2O40DRzNfb3WrH557nX5f1UYERERKRIWG02Ih8ZwZEHv+egtRp+nKXJL8/w8ZwZZObkFe17FelsIiIiUurVahRJwOBEEvy7sd5eh1E7A+k4eSWbDp8psve4qgIzdepUQkJC8PDwICIigrVr1/7t+C+//JJ69erh4eFBo0aNWLx4caHfG4bBiBEjCAwMxNPTk+joaPbu3Xs10URERKQY8PD0JqrPdHIe+ZbKvt4cOJXBv6avYNmcf5Obk33N8ztcYObNm0dsbCwjR45k48aNNGnShJiYGFJTUy85fvXq1XTr1o1evXqxadMmOnXqRKdOndi2bVvBmHHjxjF58mSmT5/OmjVr8Pb2JiYmhqysrKvfMhERETFdq7pBxA28mXubBDHQ+iW3HpjAoTejOLhz/TXN6/Bl1BEREbRo0YIpU6YAYLfbCQ4OZsCAAQwbNuyi8V26dCEjI4NFixYVLIuMjCQsLIzp06djGAZBQUE8//zzvPDCCwCkpaXh7+/P7Nmz6dq16z9m+vMyrOPHj+syahERkWJqy5I51E8aTTnLebINF5b796R9vwlXdRm1iyODc3Jy2LBhA8OHDy9YZrVaiY6OJiEh4ZLrJCQkEBsbW2hZTEwMCxcuBODAgQMkJycTHR1d8HtfX18iIiJISEi4ZIHJzs4mO/uv3U/p6ekABAUFObI5IiIicoMFlLHw/j0edKwDUUdmXvU8Dh1COnXqFPn5+fj7+xda7u/vT3Jy8iXXSU5O/tvxf/6nI3OOGTMGX1/fgldwcLAjmyEiIiImST5vcPcXF3jqPxc4l33199J1aA9McTF8+PBCe3XS09MJDg7WISQREREnsmf3Tni7xVWt61CB8fPzw2azkZKSUmh5SkoKAQEBl1wnICDgb8f/+Z8pKSkEBgYWGhMWFnbJOd3d3XF3d79oube3N97e3le8PSIiImKeWrXrXPW6Dh1CcnNzo3nz5sTHxxcss9vtxMfHExUVdcl1oqKiCo0HWLp0acH4GjVqEBAQUGhMeno6a9asueycIiIiUro5fAgpNjaWHj16EB4eTsuWLZk0aRIZGRn07NkTgO7du1OlShXGjBkDwMCBA2nXrh0TJkygY8eOzJ07l/Xr1zNz5h8n7lgsFgYNGsTrr79OaGgoNWrU4JVXXiEoKIhOnToV3ZaKiIhIieFwgenSpQsnT55kxIgRJCcnExYWRlxcXMFJuIcPH8b6X0+fbNWqFZ9//jkvv/wyL774IqGhoSxcuJCGDRsWjBkyZAgZGRk8/fTTnD17ljZt2hAXF4eHh0cRbKKIiIiUNA7fB6Y4upbHcYuIiIg5ruX7W89CEhEREaejAiMiIiJORwVGREREnI4KjIiIiDgdFRgRERFxOiowIiIi4nRUYERERMTpqMCIiIiI01GBEREREafj8KMEiqM/byacnp5uchIRERG5Un9+b1/NQwFKRIE5ffo0AMHBwSYnEREREUedPn0aX19fh9YpEQWmQoUKwB8PknT0A5DC0tPTCQ4O5siRI3qu1DXSZ1k09DkWHX2WRUefZdFIS0ujWrVqBd/jjigRBebPp1/7+vrqD6mI+Pj46LMsIvosi4Y+x6Kjz7Lo6LMsGn9+jzu0znXIISIiInJdqcCIiIiI0ykRBcbd3Z2RI0fi7u5udhSnp8+y6OizLBr6HIuOPsuio8+yaFzL52gxrubaJRERERETlYg9MCIiIlK6qMCIiIiI01GBEREREaejAiMiIiJOp0QUmKlTpxISEoKHhwcRERGsXbvW7EhO55dffuGee+4hKCgIi8XCwoULzY7klMaMGUOLFi0oW7YslStXplOnTuzevdvsWE5p2rRpNG7cuOBGYVFRUfzwww9mx3J6Y8eOxWKxMGjQILOjOJ1XX30Vi8VS6FWvXj2zYzmtY8eO8eijj1KxYkU8PT1p1KgR69evv+L1nb7AzJs3j9jYWEaOHMnGjRtp0qQJMTExpKammh3NqWRkZNCkSROmTp1qdhSntmLFCvr160diYiJLly4lNzeXO++8k4yMDLOjOZ2qVasyduxYNmzYwPr167ntttu477772L59u9nRnNa6deuYMWMGjRs3NjuK02rQoAEnTpwoeK1cudLsSE7pzJkztG7dGldXV3744Qd27NjBhAkTKF++/JVPYji5li1bGv369Sv4OT8/3wgKCjLGjBljYirnBhgLFiwwO0aJkJqaagDGihUrzI5SIpQvX954//33zY7hlM6dO2eEhoYaS5cuNdq1a2cMHDjQ7EhOZ+TIkUaTJk3MjlEiDB061GjTps01zeHUe2BycnLYsGED0dHRBcusVivR0dEkJCSYmEzkD2lpaQBX9aAy+Ut+fj5z584lIyODqKgos+M4pX79+tGxY8dC/7wUx+3du5egoCBq1qzJI488wuHDh82O5JS+++47wsPDeeihh6hcuTJNmzZl1qxZDs3h1AXm1KlT5Ofn4+/vX2i5v78/ycnJJqUS+YPdbmfQoEG0bt2ahg0bmh3HKW3dupUyZcrg7u5O7969WbBgATfddJPZsZzO3Llz2bhxI2PGjDE7ilOLiIhg9uzZxMXFMW3aNA4cOEDbtm05d+6c2dGczm+//ca0adMIDQ3lxx9/pE+fPjz77LPMmTPniucoEU+jFimO+vXrx7Zt23SM/BrUrVuXpKQk0tLS+Oqrr+jRowcrVqxQiXHAkSNHGDhwIEuXLsXDw8PsOE6tQ4cOBf+9cePGREREUL16debPn0+vXr1MTOZ87HY74eHhvPHGGwA0bdqUbdu2MX36dHr06HFFczj1Hhg/Pz9sNhspKSmFlqekpBAQEGBSKhHo378/ixYtYtmyZVStWtXsOE7Lzc2N2rVr07x5c8aMGUOTJk145513zI7lVDZs2EBqairNmjXDxcUFFxcXVqxYweTJk3FxcSE/P9/siE6rXLly1KlTh3379pkdxekEBgZe9C8i9evXd+iQnFMXGDc3N5o3b058fHzBMrvdTnx8vI6TiykMw6B///4sWLCAn3/+mRo1apgdqUSx2+1kZ2ebHcOp3H777WzdupWkpKSCV3h4OI888ghJSUnYbDazIzqt8+fPs3//fgIDA82O4nRat2590S0m9uzZQ/Xq1a94Dqc/hBQbG0uPHj0IDw+nZcuWTJo0iYyMDHr27Gl2NKdy/vz5Qv8WceDAAZKSkqhQoQLVqlUzMZlz6devH59//jnffvstZcuWLTgXy9fXF09PT5PTOZfhw4fToUMHqlWrxrlz5/j8889Zvnw5P/74o9nRnErZsmUvOgfL29ubihUr6twsB73wwgvcc889VK9enePHjzNy5EhsNhvdunUzO5rTee6552jVqhVvvPEGDz/8MGvXrmXmzJnMnDnzyicpmguizPXuu+8a1apVM9zc3IyWLVsaiYmJZkdyOsuWLTOAi149evQwO5pTudRnCBgfffSR2dGczhNPPGFUr17dcHNzMypVqmTcfvvtxpIlS8yOVSLoMuqr06VLFyMwMNBwc3MzqlSpYnTp0sXYt2+f2bGc1n/+8x+jYcOGhru7u1GvXj1j5syZDq1vMQzDKOJiJSIiInJdOfU5MCIiIlI6qcCIiIiI01GBEREREaejAiMiIiJORwVGREREnI4KjIiIiDgdFRgRERFxOiowIiIi4nRUYERERMTpqMCIiIiI01GBEREREaejAiMiIiJO5/8AFrXO0RigItcAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import ase\n",
    "from equistore import TensorMap, TensorBlock, Labels\n",
    "\n",
    "def get_dummy_structures(r_array):\n",
    "    dummy_structures = []\n",
    "    for r in r_array:\n",
    "        dummy_structures.append(\n",
    "            ase.Atoms('CH', positions=[(0, 0, 0), (0, 0, r)])\n",
    "        )\n",
    "    return dummy_structures \n",
    "\n",
    "# Create a fake list of dummy structures to test the radial functions generated from rascaline.\n",
    "\n",
    "r = np.linspace(0.1, a-0.001, 100)\n",
    "structures = get_dummy_structures(r)\n",
    "\n",
    "spherical_expansion_coefficients = calculator.compute(structures)\n",
    "\n",
    "block_C_0 = spherical_expansion_coefficients.block(species_center = 6, spherical_harmonics_l = 0, species_neighbor = 1)\n",
    "print(block_C_0.values.shape)\n",
    "\n",
    "block_C_0_0 = block_C_0.values[:, :, 0].flatten()\n",
    "spherical_harmonics_0 = 1.0/np.sqrt(4.0*np.pi)\n",
    "\n",
    "\n",
    "all_species = np.unique(spherical_expansion_coefficients.keys[\"species_center\"])\n",
    "all_neighbor_species = Labels(\n",
    "        names=[\"species_neighbor\"],\n",
    "        values=np.array(all_species, dtype=np.int32).reshape(-1, 1),\n",
    "    )\n",
    "spherical_expansion_coefficients.keys_to_properties(all_neighbor_species)\n",
    "\n",
    "plt.plot(r, block_C_0_0/spherical_harmonics_0, label=\"rascaline output\")  # rascaline BUG?\n",
    "plt.plot([0.0, a], [0.0, 0.0], \"black\")\n",
    "plt.plot(r, get_LE_function_python(n=0, l=0, r=r), \"--\", label=\"original function\")\n",
    "plt.xlim(0.0, a)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "??\n",
      "TensorMap with 28 blocks\n",
      "keys: ['spherical_harmonics_l' 'species_center']\n",
      "                  0                   1\n",
      "                  1                   1\n",
      "                  2                   1\n",
      "                  3                   1\n",
      "                  4                   1\n",
      "                  5                   1\n",
      "                  6                   1\n",
      "                  7                   1\n",
      "                  8                   1\n",
      "                  9                   1\n",
      "                 10                   1\n",
      "                 11                   1\n",
      "                 12                   1\n",
      "                 13                   1\n",
      "                  0                   6\n",
      "                  1                   6\n",
      "                  2                   6\n",
      "                  3                   6\n",
      "                  4                   6\n",
      "                  5                   6\n",
      "                  6                   6\n",
      "                  7                   6\n",
      "                  8                   6\n",
      "                  9                   6\n",
      "                 10                   6\n",
      "                 11                   6\n",
      "                 12                   6\n",
      "                 13                   6\n"
     ]
    }
   ],
   "source": [
    "import torch\n",
    "\n",
    "n_max_l = []\n",
    "for l in range(l_max+1):\n",
    "    n_max_l.append(np.where(E_nl[:, l] <= E_max)[0][-1]+1) \n",
    "\n",
    "def cut_to_LE(map: TensorMap, E_nl, E_max) -> TensorMap:\n",
    "    LE_blocks = []\n",
    "    for idx, block in map:\n",
    "        l = idx[0]\n",
    "        counter = 0\n",
    "        for n in block.properties[\"n\"]:\n",
    "            if E_nl[n, l] <= E_max: counter += 1\n",
    "        LE_values = torch.zeros((block.values.shape[0], block.values.shape[1], counter))\n",
    "        counter_LE = 0\n",
    "        counter_total = 0\n",
    "        labels_LE = [] \n",
    "        for n in block.properties[\"n\"]:\n",
    "            if E_nl[n, l] <= E_max: \n",
    "                LE_values[:, :, counter_LE] = torch.tensor(block.values[:, :, counter_total])\n",
    "                labels_LE.append([block.properties[\"species_neighbor\"][counter_total], n, l])\n",
    "                counter_LE += 1\n",
    "            counter_total += 1\n",
    "        LE_block = TensorBlock(\n",
    "            values=LE_values,\n",
    "            samples=block.samples,\n",
    "            components=block.components,\n",
    "            properties=Labels(\n",
    "                names = (\"a1\", \"n1\", \"l1\"),\n",
    "                values = np.array(labels_LE),\n",
    "            ),\n",
    "        )\n",
    "        LE_blocks.append(LE_block)\n",
    "    return TensorMap(\n",
    "            keys = Labels(\n",
    "                names = (\"lambda\", \"a_i\"),\n",
    "                values = map.keys.asarray(),\n",
    "            ),\n",
    "            blocks = LE_blocks\n",
    "        )\n",
    "\n",
    "spherical_LE = cut_to_LE(spherical_expansion_coefficients, E_nl, E_max)\n",
    "print(\"??\")\n",
    "\n",
    "print(spherical_LE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.7.7 ('virtualenv': venv)",
   "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.7.7"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "144ec3b5cd37a648e305ae2c09b9516066242b1bc5e9f237ba5a8f13b411ab3c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}