examples/greedy_atomic_density.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "pycharm": {
     "is_executing": false
    }
   },
   "source": [
    "# Greedy Method For Fitting Beryllium's Atomic Density\n",
    "\n",
    "\n",
    "This file will showcase how to use the greedy method for optimizing the Kullback-Leibler divergence functions\n",
    "that iteratively selects the next best $(n + 1)$ basis functions from the previous $n$ basis functions.\n",
    "\n",
    "There are three different choices of picking the next set of basis-functions, the one we're going to use\n",
    "in this example is the `pick_one`. This generetes the next set of initial guesses by taking the average between\n",
    "each consequent exponent, and multiplies the smallest and largest exponent by a constant $C$.  For example,\n",
    "for three exponents $\\{\\alpha_1, \\alpha_2, \\alpha_3\\}$, the three next initial guesses are\n",
    "\n",
    "$$\\mathcal{C} = \\bigg\\{ \\{C\\alpha_1, \\alpha_2, \\alpha_3 \\}, \\big\\{\\alpha_1, \\frac{\\alpha_1 + \\alpha_2}{2}, \\alpha_3 \\big\\}, \\big\\{ \\alpha_1, \\alpha_2, \\frac{\\alpha_2 + \\alpha_3}{2}, \\alpha_3 \\big\\}, \\big\\{\\alpha_1, \\alpha_2, C \\alpha_3\\big\\} \\bigg\\}. $$\n",
    "\n",
    "The greedy algorithm would optimize each initial guess in $\\mathcal{C}$ and picks the best \"local\" answer found and further\n",
    "optimizes it.  Then, this optimal answer is then used for the next `5` basis-functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Specify a radial grid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is to specify a one-dimensional (radial) grid.\n",
    "\n",
    "The ClenshawCurtis grid is chosen here and is recommended for radial grids since it centers the points where the electron density is thought to be high.\n",
    "However, just to make sure let's add additional points [25, 50, 75]. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from bfit.grid import ClenshawRadialGrid\n",
    "\n",
    "atomic_numb = 4\n",
    "grid = ClenshawRadialGrid(atomic_numb, num_core_pts=400, num_diffuse_pts=500, \n",
    "                          extra_pts=[25, 50, 75])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Get Slater density."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to get the atomic density of Beryllium.\n",
    "\n",
    "Additionally, let's integrate it spherically to make sure the integral of the atomic density is indeed close to the atomic number $4$ and let's plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The integral of potassium is:  3.9999910287814453894\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAi0ElEQVR4nO3debRcVZn38e/PEIQQIEACIpJEaRABBfVK49A0QxREbfV19oJR7DdNUNF2xNZ+hWXT2rRTqwts1ECEOKCCosYhoBEHRG8kBEJAFJMgQXIJIASUIXneP86uULmp4dzh1Kmq8/usVauqzlTPvkk9Z599du2tiMDMzKrjMWUHYGZmneXEb2ZWMU78ZmYV48RvZlYxTvxmZhXjxG9mVjFO/NYxkmZK2ihpUtmx5CXp+5Lmduiz/kPSnZL+3InPaxPLGZIuSq9nSwpJ26X3HfubWDGc+CtK0lJJd0t67IjlqyXNKeIzI2JtREyNiE2j2U/SGyVtSieNjZL+KOl8SQcUEWe9iHhhRCysi+PnRXyOpH2BdwEHRcTjGqw/StLmur/BbZLOLCKWdur/JtabnPgrSNJs4B+AAP6p3GhyuyoipgK7AnOAvwLLJB1SblgTZhawISLWt9hmXTpxTgWeB7xZ0svG8mG9dNVlE8+Jv5reAPwKuADYcsku6UJgJvCdVKt8b1r+T5JWSronXSk8pW6f1ZLeI2mFpPslfVHSXqk54D5Jl0vaLW07sslg91RzX5euPr7VLvCI2BQRf4iIU4GfAmfUxXKEpF+mOK+VdFTduqWSPizpFymuH0mantbtIOkiSRvSvr+RtFfdfv+cyvw54Nnpb3OPpGdJuqNWnrT9KyQtbxS7pF0lfUnSsKQ1kj4o6THpCmsJ8Ph07Aty/B3+CPwSOKju+AdKWiLpLkk3SXp13boLJJ0rabGk+4F3jib2EeVYKumf0+stTULp/ch/46WpCeuXqWzfkbSHpEWS7k1/69ntPtMmlhN/Nb0BWJQex9WSXEScBKwFXpJqlmen5pSvAO8AZgCLyU4M29cd7xXA84EDgJcA3wf+DZhO9n/stCZxXAhMAQ4G9gQ+OcpyXEJ25YKkfYDvAf8B7A68G/impBl1278eeFP6rO3TNpCd/HYF9gX2AE4hu6LYIiJWpeVXpb/NtIj4DbAhlb3mxFSuRj6TPudJwD+S/Tu8KSIuB17IozX6N7YruKT9geeSncCRtBPZyePLqXyvA86RdPCI8p8F7JxiGU3s4/Fa4CRgH2A/4CrgfLJ/p1XAhwr4TGvBib9iJD2PrFnh4ohYBvyBLCE08xrgexGxJCIeBj4G7Ag8p26bz0TEHRFxG/Az4OqIuCYiHgQuBZ7eII69yZLdKRFxd0Q8HBE/HWVx1pElD8iS1uKIWBwRmyNiCTAEnFC3/fkR8buI+CtwMXBYWv4wWcL/u3RFsSwi7s0Zw8L02UjaHTiOLPluJTWtvAZ4f0TcFxGrgY+TJcS8Hp+uNO4FfgdcDdTuObwYWB0R50fEIxHxW+CbwCvr9v92RPwi/X3+ljf2CXB+ukr7C1ml4A8RcXlEPAJ8nQb/P6xYTvzVMxf4UUTcmd5/mbrmngYeD6ypvYmIzcCtZLW3mjvqXv+1wfupDY67L3BXRNydP/Rt7APclV7PAl6VEuM9ku4hawffu277+t4yD9TFdSHwQ+CrqdnpbEmTc8ZwEfASSVOBVwM/i4jbG2w3newqY03dsjVs/XdsZ1260tgFmEb2t63dZJ0F/P2I8g8C9TeKbx1j7OM1lv8fVqDt2m9i/ULSjmRf8El6tMvgY4Fpkg6NiGvJbvjWWwc8te4YIkvat40znFuB3SVNi4h7xniMl5NdYdSOd2FE/N/RHiRdyZwJnJnamxcDNwFfHLlpg31vk3RViuUk4NwmH3Mn2ZXFLOCGtGwmY/w7RsRfJH0Z+FpadCvw04h4fqvdxhh7K/eTNdfVbNMjybqPa/zV8jJgE9kNwcPS4ylkyfMNaZs7yNqgay4GXiTp2FQLfhfwINmNxTFLNcvvk7VD7yZpsqQj2+0naZKkJ0r6DHAUWcKGR2uvx6VtdlDWBfIJOY55tKSnpuaYe8kSdKMup3cATxhxfwPgS8B7yU6QlzYp7yayv+VZknaWNAt4Z4p71FIt/bXAyrTou8ABkk5Kf8vJ6ebzU5ofJV/sbSwHjlT2G41dgfeP4RjWYU781TKXrL11bUT8ufYAPgsMpp4YHwE+mJoL3h0RN5G1A3+GrNb6ErKbvw9NQDwnkSXZG4H1ZDeQm3m2pI1kiXkpsAvwrIi4DiAibgVeSnZTeZisBvwe8v0ffxzwjXTsVWS9hRol5B+TJdo/S7qzbvmlZDX5SyPi/haf8zayGvItZG3zXwYW5IivptbrZyNZM9HuZM05RMR9wAvITgbryJq1/ovsiq6VvLE3lO6lfA1YASwjOwGNmbLeY4PjOYa1J0/EYjZ+kv4A/EvqodNTejl2GxvX+M3GSdIryNrPf1x2LKPVy7Hb2Pnmrtk4SFpKds/kpNTjqWf0cuw2Pm7qMTOrGDf1mJlVTE809UyfPj1mz55ddhhmZj1l2bJld0bEjJHLeyLxz549m6GhobLDMDPrKZLWNFruph4zs4px4jczqxgnfjOzinHiNzOrGCd+M7OK6YlePWMibbvMP1YzM+vTGn+jpN9quZlZhRSW+NN46L9WNun1SklnpuVnSLpN0vL0OKHdsczMbOIUWeN/EDgmIg4lm/DjeElHpHWfjIjD0mNxgTFsa9Gijn6cmVm3KSzxR2Zjejs5PcpvZD/xxLIjMDMrVaFt/GkKvOVksystiYir06q3SlohaYGk3ZrsO0/SkKSh4eHhIsM0M6uUQhN/RGyKiMOAJwCHSzqEbELn/ciaf24HPt5k3/MiYiAiBmbM2GaMITMzG6OO9OqJiHvI5kk9PiLuSCeEzcDngcM7EYOZmWWK7NUzQ9K09HpHYA5wo6S96zZ7OXD9hH+4++ubmTVV5A+49gYWSppEdoK5OCK+K+lCSYeR3ehdDfxLgTGYmdkIhSX+iFgBPL3B8pOK+szc9tkHbrut7CjMzErRn7/cbWfdurIjMDMrTTUTv5lZhfVv4n/848uOwMysK/Vv4ncbvplZQ/2b+M3MrKHqJv6DDy47AjOzUlQ38d9wQ9kRmJmVorqJ38ysovo78R90UNkRmJl1nf5O/CtXlh2BmVnX6e/E345n4zKzCqp24j/55LIjMDPruGon/oceKjsCM7OOq3biNzOroP5P/PPnlx2BmVlX6f/Ef845ZUdgZtZV+j/xt3PqqWVHYGbWUU78555bdgRmZh3lxG9mVjGFJX5JO0j6taRrJa2UdGZavrukJZJuTs+7FRWDmZltq8ga/4PAMRFxKHAYcLykI4DTgSsiYn/givS+WO7ZY2a2RWGJPzIb09vJ6RHAS4GFaflC4GVFxbCFe/aYmW1RaBu/pEmSlgPrgSURcTWwV0TcDpCe92yy7zxJQ5KGhoeHiwwT5swp9vhmZl2k0MQfEZsi4jDgCcDhkg4Zxb7nRcRARAzMmDGjsBgBuOKKYo9vZtZFOtKrJyLuAZYCxwN3SNobID2v70QMZmaWKbJXzwxJ09LrHYE5wI3AZcDctNlc4NtFxbAVT8piZgbAdgUee29goaRJZCeYiyPiu5KuAi6W9GZgLfCqAmN41MqVIHXko8zMullhiT8iVgBPb7B8A3BsUZ87Zqee6t4/ZlYJ/uVujYduMLOKcOI3M6uYaiX+Y7uvhcnMrNOqlfgvv7zsCMzMSletxN+Ox+Y3swpw4q/nG7xmVgFO/GZmFVO9xO8bvGZWcdVL/L7Ba2YVV73E346HaDazPufEP5KHaDazPufEb2ZWMdVM/J6D18wqrJqJ36NwmlmFVTPxt+MbvGbWx5z4G/ENXjPrY9VN/J6Ny8wqqrqJ/5RTyo7AzKwU1U387W7weqROM+tThSV+SftK+omkVZJWSnp7Wn6GpNskLU+PE4qKYVw8UqeZ9anCJlsHHgHeFRG/lbQzsEzSkrTukxHxsQI/28zMmigs8UfE7cDt6fV9klYB+xT1eWMydSps3Fh2FGZmHdWRNn5Js4GnA1enRW+VtELSAkm7NdlnnqQhSUPDw8PFBPa5zxVzXDOzLlZ44pc0Ffgm8I6IuBc4F9gPOIzsiuDjjfaLiPMiYiAiBmbMmFFMcIODrdf7h1xm1ocKTfySJpMl/UURcQlARNwREZsiYjPweeDwImMYF/+Qy8z6UJG9egR8EVgVEZ+oW7533WYvB64vKgYzM9tWkTX+5wInAceM6Lp5tqTrJK0Ajgb+tcAY2vNInWZWMYqI1htIbyVrqrm7MyFta2BgIIaGhor7gFbDN1x0Uft7AWZmXUjSsogYGLk8T43/ccBvJF0s6fjUhFMdJ59cdgRmZhOqbeKPiA8C+5O1178RuFnSf0rar+DYusNDD5UdgZnZhMrVxh9Ze9Cf0+MRYDfgG5LOLjC2ztlhh7IjMDPrmLaJX9JpkpYBZwO/AJ4aEfOBZwKvKDi+zvjCF1qvX7SoM3GYmXVAnhr/dOD/RMRxEfH1iHgYIPXDf3Gh0XVKu5u3c+d2Jg4zsw7Ik/ifGBFr6hdIuhAgIlYVElW32bSp7AjMzCZMnsR/cP0bSZPImnn6i9v5zawimiZ+Se+XdB/wNEn3psd9wHrg2x2LsFPczm9mFZHnB1wfiYj3dyiehgr/AVdNq58o7LSTh3A2s57S7AdcTcfjl3RgRNwIfF3SM0auj4jfTnCM3e3++8uOwMxsQrSaiOWdwDwaD5scwDGFRFQmCdpcAZmZ9bqmbfwRMS89H93g0X9JH+CUU1qv9wTsZtYH8vyA61VpzlwkfVDSJZKeXnxoJTjnnNbrPQG7mfWBPN05/z3Nmfs84DhgIeA5C83MelSexF/79dKLgHMj4tvA9sWFVLKpU8uOwMysUHkS/22S/hd4NbBY0mNz7teb2k3A7nZ+M+txefrxTwGOB66LiJvT1IlPjYgfdSJA6GA//pp2Uw6454+Z9YAxT8QSEQ+Q/VL3fkkzgcnAjRMfopmZdUKrfvwASHob8CHgDmBzWhzA0wqMq1w77AB/+1vZUZiZFSJPW/3bgSdHxMER8dT0aJv0Je0r6SeSVklaKentafnukpZIujk97zbeQky4duP2HHxw6/VmZl0sT+K/FfjLGI79CPCuiHgKcATwFkkHAacDV0TE/sAV6X13aTc+/w03dCYOM7MCtG3qAW4Blkr6HvBgbWFEfKLVThFxO3B7en2fpFXAPsBLgaPSZguBpcD7Rhu4mZmNTZ4a/1pgCVnf/Z3rHrlJmg08Hbga2CudFGonhz2b7DNP0pCkoeHh4dF83MSYP7/1enfrNLMe1bY755YNpZ0iYtRDVEqaCvwUOCsiLpF0T0RMq1t/d0S0bOfveHfOGnfrNLMeNubunJKeLekGYFV6f6ikNoPabNl3MvBNYFFEXJIW35F+C0B6Xp+zDGZmNgHyNPV8imyMng0AEXEtcGS7nSQJ+CKwasT9gMuA2uzlc+nm2bw8HaOZ9aFcQy9ExK0jFuWZffy5wEnAMZKWp8cJwEeB50u6GXh+et+d2nXrnDOnM3GYmU2gPL16bpX0HCAkbQ+cRmr2aSUifg40ayQ/Nn+IJRochBNPbL7+iis6F4uZ2QTJU+M/BXgLWVfMPwGHpfdmZtaD8ozVc2dEDEbEXhGxZ0ScGBEbOhFcVzi2zcWJu3WaWY9pmfglHZ1m3FqZHt+QdFRnQusSl1/eer1n5TKzHtM08Ut6EbAA+A7wemAQWAwsSDdpzcysB7Wq8b8HeFlEnB8R10bE8ohYALyMqg2x0G5WLjf3mFkPaZX4H5f67G8lIlYAexUXUhdqNyuXm3vMrIe0SvythmcY9dANPa3daJ1mZj2kVT/+/SRd1mC5gCcVFE/3mjoVNm5svv7UU+GcXCNZmJmVqukgbZL+sdWOEfHTQiJqoLRB2uotWtT6x1zgQdvMrKs0G6Qt9+icZeqKxA8erdPMesqYR+e0Ou0GbXPvHjPrAU78o9Fu0Db37jGzHuDEPxru3WNmfSDPRCxLJE2re7+bpB8WGlU3c3OPmfW4PDX+6RFxT+1NRNxNk3lyK8HNPWbW4/Ik/s2SZtbeSJoFVLf7ipt7zKzH5ZmI5QPAzyXV+u0fCcwrLqQesN128MgjzdfPmdN+VE8zs5LkGY//B8AzgK8BFwPPjIjqtvEDXHBB6/WemcvMulirYZkPTM/PAGYC64DbgJlpWXXlae5ZtKj4OMzMxqBVjf+d6fnjDR4fa3dgSQskrZd0fd2yMyTdNmLy9d7Ubmauk0/uTBxmZqPUtI0/Iual56PHeOwLgM8CXxqx/JMR0fbE0fUuv7z1EA4PPdS5WMzMRqHtzV1Jk4AXAbPrt4+IT7TaLyKulDR7nPH1tkWL3AvIzLpOnu6c3wHeCOwB7Fz3GKu3SlqRmoJ2a7aRpHmShiQNDQ8Pj+PjCjR/fuv1c+d2Jg4zs1FoOzqnpBUR8bQxHTyr8X83Ig5J7/cC7iT7HcCHgb0jom1jeNeMztmIR+w0sy41ntE5vy/pBRMRRETcERGbImIz8Hng8Ik4blebM6fsCMzMtpIn8f8KuFTSXyXdK+k+SfeO5cMk7V339uXA9c227Rntmnvcp9/MukyexP9x4NnAlIjYJSJ2johd2u0k6SvAVcCTJf1J0puBsyVdJ2kFcDTwr+MJvivkmW7RtX4z6yJ5hmy4Gbg+RjlVV0S8rsHiL47mGD3joIPghhuar3et38y6SJ7EfzuwVNL3gQdrC9t156yUlSvb3+Q1M+sSeZp6/ghcAWzPxHTnrKaDDy47AjMzIEeNPyLOBJC0c/Y2NhYeVS+aP7/1WPytmoLMzDoozwxch0i6hqwHzkpJyyS5+jpSnpu8HrjNzLpAnqae84B3RsSsiJgFvIusD76N1G7gtg98oDNxmJm1kCfx7xQRP6m9iYilwE6FRdTL2k2+snZtZ+IwM2shT+K/RdK/S5qdHh8ku+FrjbT6QddjHuPmHjMrXZ7EfzIwA7gkPaaTDdpmjZxzTvMmn02bYN48J38zK1WexD8nIk6LiGekxzuA5xccV2+7/HK46CKYNGnbdQ884LZ+MytVnsT//pzLrN7gIGze3Hid2/rNrERN+/FLeiFwArCPpE/XrdoFeKTowPrCzJmwZk3j5WZmJWlV418HDAF/A5bVPS4DJmSY5r531lkwZcq2yzdudDu/mZWm1Zy71wLXSvpyRDxcWy7pecBZwFs6EF9vq027+Pa3w4YNjy7fsCG7yVu/jZlZh7Rt44+IhyUdJulsSavJZs66sfDI+sXgIEyduu1y3+Q1s5K0auM/AHgt8DpgA/A1sqkaj+5QbP2j2c1c3+Q1sxK0qvHfCBwLvCQinhcRnwE2dSasPtPsZq5/0GVmJWiV+F8B/Bn4iaTPSzoW8KDzY9HsJq9/0GVmJWia+CPi0oh4DXAgsJRsmsS9JJ07UZOvV8bgIJx3nn/QZWZdIc/N3fsjYlFEvBh4ArAcOL3owPqOf9BlZl0izy93t4iIuyLifyPimHbbSlogab2k6+uW7S5piaSb0/NuYwm6Z7mt38y6wKgS/yhdABw/YtnpwBURsT/ZdI7VunJwW7+ZdYHCEn9EXAncNWLxS4GF6fVC4GVFfX5Xclu/mXWBImv8jewVEbcDpOc9m20oaZ6kIUlDw8PDHQuwcK3a+tesca3fzArX6cSfW0ScFxEDETEwY8aMssOZWK0GaXOTj5kVrNOJ/w5JewOk5/Ud/vzu0KytH9zkY2aF63TivwyYm17PBb7d4c/vDrW2/mbcvdPMClRY4pf0FeAq4MmS/iTpzcBHgedLuplsFq+PFvX5XW9wEGbNarzO3TvNrEBNB2kbr4h4XZNVTSakraCzzsra9B94YOvlte6d4GGbzWzCde3N3Upw904zK4ETf9nade+cPdvNPmY2oZz4u0Gr7p1r1riLp5lNKCf+btCqeye42cfMJpQTfzeotfU36+UD7uJpZhPGib9bDA7C6tXu4mlmhXPi7zYewdPMCubE323cxdPMCubE343cxdPMCuTE363cxdPMCuLE363cxdPMCuLE363ydPH0xC1mNgZO/N2sXRdPcJOPmY2aE38v8MQtZjaBnPh7QbuJW9zTx8xGwYm/V7SauAXc08fMcnPi7yV5evrMnevkb2YtOfH3kjw9fTy0g5m14cTfa/L09PENXzNroZTEL2m1pOskLZc0VEYMPa9ds49v+JpZE4VNtp7D0RFxZ4mf39tqk7DPnZs17zRSu+Fbv72ZVZ6benrZ4CAsXNj+hu+JJ7r2b2ZblJX4A/iRpGWS5jXaQNI8SUOShoaHhzscXg/Jc8MX3N3TzLZQRHT+Q6XHR8Q6SXsCS4C3RcSVzbYfGBiIoSHfCmhr9uwswbcza1Z2j8DNP2Z9TdKyiBgYubyUGn9ErEvP64FLgcPLiKPvtLvhW+Pav1mldTzxS9pJ0s6118ALgOs7HUdfytvsA277N6uwMmr8ewE/l3Qt8GvgexHxgxLi6E+1fv4XXeTav5k11PHEHxG3RMSh6XFwRJzV6RgqwbV/M2vC3Tn7mWv/ZtaAE38VuPZvZnWc+KtiLLX/k04CyScBsz7jxF81o6n9137j4SYgs77ixF9Fo639g5uAzPqIE3+Vjab2X+MmILOe58RfdWOp/bsJyKynOfFbZmTtX8q3n5uAzHqOE789qlb7j4ALLxx9E9Cb3gTTp8NjHuMTgVkXc+K3xsbSBPTww7BhQ3bi8L0As67lxG+tjbUJCLa+F+CrAbOu4cRv7Y2nCajGVwNmXcOJ30ZnLE1AjTS7Gpg+3VcGZgVz4rexqW8CkmCPPWD77cd2rPqrgQ0btr4ycBOR2YRz4rexq9X+N2+GO++EBQvGdi+glZFNRCNPBKeemj37xGCWmxO/TZxG9wLGezUw0sgTwbnnZs+N7h34pGDWkBO/FaMTVwON1N87GHlSyHO1sGiRTxbW9xS1L0oXGxgYiKGhobLDsImyaBF84AOwdi3svjvcdx889FDZUcHkydlJqT6WyZNhl13grrtg5kw44QRYvDiLvdH7s9KEcrXy1ZYNDpZTJqs0ScsiYmCbFRHR9Y9nPvOZYX3soosiZs2KkCL22CN71F5vv31EVmfvjcfkydvGPHnyo2WaNSti/vxHyztrVlb++r9Bs2Vj3a/ZMut7wFA0yKmlJHLgeOAm4PfA6e22d+KvsJEnhV47EYznZNGurHn3y3MymoiTyHj2H+8xO71tmctHoWsSPzAJ+APwJGB74FrgoFb7OPHbFs1qwZAtKzuJ9+pju+2yk0G7ZZMnb7tsxx0jTj45e56oE9OUKdm/7ZQp3bdtmctHmfybJf6Ot/FLejZwRkQcl96/PzU5faTZPm7jt1zq7x2MbH/Pcy+hURu/WTeZNSvrNJFTszb+Mnr17APcWvf+T2nZViTNkzQkaWh4eLhjwVkPq+9JtHo1nHNO455FUvY8f/7W788/f+ttxtINdfLkieu6ajbS2rUTcpgyEn+jvnzbXHZExHkRMRARAzNmzOhAWNb3Wp0YVq/O1jfrhtrsZDERJ49GJ4s8J5C8+3XiZDRpUrnH7PS2ZS2fObPx8lEqI/H/Cdi37v0TgHUlxGHWXruTxUScPEaeLJotG+t+eU5G4zmJTJmSzcQ2cuymTh2z09uWubzWXXi8GjX8F/kAtgNuAZ7Iozd3D261j2/umk2wMnvgdGtPnW7rvVNgr55SfsAl6QTgU2Q9fBZERMvTmG/umpmNXrObu9uVEUxELAYWl/HZZmZV57F6zMwqxonfzKxinPjNzCrGid/MrGJ6YlhmScPAmjHuPh24cwLDKYvL0T36oQzgcnSbIsoxKyK2+QVsTyT+8ZA01Kg7U69xObpHP5QBXI5u08lyuKnHzKxinPjNzCqmCon/vLIDmCAuR/fohzKAy9FtOlaOvm/jNzOzrVWhxm9mZnWc+M3MKqZvEr+k4yXdJOn3kk5vsF6SPp3Wr5D0jDLibCdHOQZT/Csk/VLSoWXE2Uq7MtRt9yxJmyS9spPx5ZWnHJKOkrRc0kpJP+10jHnk+D+1q6TvSLo2leNNZcTZiqQFktZLur7J+l75frcrR2e+343Gau61BzkmcAdOAL5PNgPYEcDVZcc9xnI8B9gtvX5ht5UjTxnqtvsx2Sitryw77jH+W0wDbgBmpvd7lh33GMvxb8B/pdczgLuA7cuOfUSMRwLPAK5vsr7rv985y9GR73e/1PgPB34fEbdExEPAV4GXjtjmpcCXIvMrYJqkvTsdaBttyxERv4yIu9PbX5HNYNZN8vxbALwN+CawvpPBjUKecrweuCQi1gJERDeWJU85AthZkoCpZIn/kc6G2VpEXEkWVzO98P1uW45Ofb/7JfHnmcA91yTvJRttjG8mq+V0k7ZlkLQP8HLgcx2Ma7Ty/FscAOwmaamkZZLe0LHo8stTjs8CTyGbAvU64O0Rsbkz4U2YXvh+j1Zh3+9SJmIpQJ4J3HNN8l6y3DFKOprsP8bzCo1o9PKU4VPA+yJiU1bJ7Ep5yrEd8EzgWGBH4CpJv4qI3xUd3CjkKcdxwHLgGGA/YImkn0XEvQXHNpF64fudW9Hf735J/HkmcO+FSd5zxSjpacAXgBdGxIYOxZZXnjIMAF9NSX86cIKkRyLiWx2JMJ+8/6fujIj7gfslXQkcCnRT4s9TjjcBH42sYfn3kv4IHAj8ujMhTohe+H7n0onvd7809fwG2F/SEyVtD7wWuGzENpcBb0h3/48A/hIRt3c60DbalkPSTOAS4KQuq1nWtC1DRDwxImZHxGzgG8CpXZb0Id//qW8D/yBpO0lTgL8HVnU4znbylGMt2VULkvYCngzc0tEox68Xvt9tder73Rc1/oh4RNJbgR/y6ATuKyWdktZ/jqz3yAnA74EHyGo5XSVnOf4fsAdwTqoxPxJdNDJhzjJ0vTzliIhVkn4ArAA2A1+IiIbd9MqS89/jw8AFkq4jazJ5X0R01TDHkr4CHAVMl/Qn4EPAZOid7zfkKkdHvt8essHMrGL6panHzMxycuI3M6sYJ34zs4px4jczqxgnfjOzinHit56XRvhcLun6NMrktFHuv1TSQHq9uN3+klZLmt5k+XVplMsfSXpci2MMSPp0m8+ZJunUnMUwy82J3/rBXyPisIg4hGwArLeM9UARcUJE3DOOWI6OiEOBIbJRL5t9zlBEnNbmWNMAJ36bcE781m+uIg3OJenwNKb5Nen5yWn5jpK+msY8/xrZODukdVtq85K+lQZfWylp3ijjuBL4O0k7SDo/XQlck8ZgqY3j/930+ow0TvtSSbdIqp0QPgrsl65m/lvS3pKurLu6+Yfx/KGsuvril7tmAJImkQ098MW06EbgyPTr1TnAfwKvAOYDD0TE09K4KL9tcsiTI+IuSTsCv5H0zVGMnfJispEu3wIQEU+VdCDwI0kHNNj+QOBoYGfgJknnAqcDh0TEYal87wJ+GBFnpbJOyRmL2Vac+K0f7ChpOTAbWAYsSct3BRZK2p9spMbJafmRwKcBImKFpBVNjnuapJen1/sC+wPtEv9PJG0iG8bhg8D5wGfSZ90oaQ3ZcM4jfS8iHgQelLQe2KvBNr8BFkiaDHwrIpa3icWsITf1WD/4a6oVzyKbZarWxv9h4Cep7f8lwA51+7Qcq0TSUcAc4Nmpzf6aEfs3c3S63/CGdK8g77jTD9a93kSDSlmaxONI4DbgQnXn+P/WA5z4rW9ExF+A04B3p1rxrmRJEuCNdZteCQwCSDoEeFqDw+0K3B0RD6QmmiPGGFb9Zx0AzARuyrnvfWRNP6T9ZwHrI+LzZM1ZXTmvrHU/J37rKxFxDdm8sq8FzgY+IukXZCNT1pwLTE1NPO+l8bjzPwC2S9t8mGwavLE4B5iURr78GvDG1KSTpywbgF+kG7n/TTaq43JJ15Ddq/ifMcZkFefROc3MKsY1fjOzinHiNzOrGCd+M7OKceI3M6sYJ34zs4px4jczqxgnfjOzivn/+01Fh4oWUUgAAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from bfit.density import SlaterAtoms\n",
    "\n",
    "# Grab information about potassium.\n",
    "be = SlaterAtoms(\"be\")\n",
    "\n",
    "# Get the total electronic density.\n",
    "density = be.atomic_density(grid.points, mode=\"total\")\n",
    "\n",
    "print(\"The integral of potassium is: \", grid.integrate(density * 4.0 * np.pi * grid.points**2.0))\n",
    "\n",
    "# Plot only the first 500 points.\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(grid.points[:500], density[:500], \"ro-\")\n",
    "plt.xlabel(\"Radial Points\")\n",
    "plt.ylabel(\"Atomic Density\")\n",
    "plt.title(\"Atomic Density of Beryllium.\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Construct Greedy-KL Object."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All the pieces are in place to provide to the \"KLDivergenceFPI\" class. Here, want to make sure that `spherical=True` since we're dealing with radial atomic densities. The `pick-one` attribute specifies how the next set of basis-functions are chosen and finally the `l_maxiter` and `g_maxiter` specifies the number of iterations needed for the local choice and global choice.  We specify a large number for the global choice so that it optimizes and refines the best found local choice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bfit.greedy import GreedyKLFPI\n",
    "\n",
    "greedy = GreedyKLFPI(grid, density, \"pick-one\", l_maxiter=2500, g_maxiter=5000, spherical=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "------------------------------------------------------------------------------------------------------------------------------------------------\n",
      "|# of Functions | # of S-type | # of P-type |  Integration  |      L1       |   L Infinity   | Least-squares |Kullback-Leibler|Change Objective|\n",
      "------------------------------------------------------------------------------------------------------------------------------------------------\n",
      "|       1       |      1      |      0      |   3.999983    | 3.758501e+00  |  3.524177e+01  | 7.943230e+00  |  4.733178e+00  |      nan       |\n",
      "|       2       |      2      |      0      |   3.999907    | 7.253813e-01  |  2.631367e+01  | 1.031142e+00  |  1.394058e-01  |  4.575385e+00  |\n",
      "|       3       |      3      |      0      |   3.999924    | 3.799149e-01  |  1.722084e+01  | 1.043090e-01  |  4.636776e-02  |  9.294167e-02  |\n",
      "|       4       |      4      |      0      |   3.999953    | 1.323686e-01  |  1.620842e+01  | 7.630855e-02  |  4.721414e-03  |  4.135387e-02  |\n",
      "------------------------------------------------------------------------------------------------------------------------------------------------\n",
      "Successful?: True\n",
      "Number of s-type: 4\n",
      "Number of p-type: 0\n",
      "Termination Information: Max number of functions reached,  4\n",
      "Objective Function: 0.004717692886119254\n",
      "Coefficients [0.56955324 1.56881298 1.2997459  0.56187837]\n",
      "Exponents [ 0.10919557  0.256764    6.49534329 28.40436813]\n",
      "[0.56955324 1.56881298 1.2997459  0.56187837]\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAtg0lEQVR4nO3deZhcVbnv8e8vSYeQBASSMJiQbg4iEELSQEQEhDCIyCB40XtPbGa8LUEGRT2geBQH1MPRg8eBcEBBIBFUREEBFbgJgwPSwSQQgqKYhACSpglTEiDpfu8ftTup7tSwu9M1df0+z1NPVa3aw7uru969au1VaykiMDOz+jGk0gGYmVl5OfGbmdUZJ34zszrjxG9mVmec+M3M6owTv5lZnXHir2OSXpP0LyXc/qWSZpdq+/VK0lWS/r3ScaQl6d2S/lLpOGwjJ/46IGmppLVJou++vTUiRkfEU8kyP5T0lRzrHVmZqIvrdTxdvY6xpQz7ny4pJP1br/LTJT1Yqv1GxNkR8eW+rpf1f/CqpJck/V7S2ZJKmgci4oGI2L1XHFX7f1UPnPjrx/FJou++PVvpgDZX9vEAy+l5jHO6l5M0rEQhnAa8mNzXiuMjYiugEfg6cBHwg8qGZOXmxF/Hktrq2yS1Ai3AvyW15V9KuhGYCPwyKfu3ZJ0DkpriS5IWSpqetb1dJN2X1CjvBsYW2PcSScdlPR8m6QVJ+0oaIWm2pI5kPw9L2qEPxzVd0gpJF0n6J3Bdrlp49/Enj7eQ9A1JyyU9nzSnbFlgHyOBDwIfA3aTNC0p3xO4CnhX8r69lJS/RdINktolLZP0ue6adhLb7yRdkRzvU5IOTMqflrRS0mlZ++7x7UzSCZIWSHpF0t8lHV3sPYqIlyPiduD/AKdJmlzsfch6Xz+ZxPScpDOy4jhG0uPJ3/8ZSZ/KXi95vMn/laQ7JJ3X6/1dJOnEYsdh/ePEb0TE1cAc4PKktnx8RJxCz1r05ZLGA3cAXwG2Az4F/EzSuGRTPwLmk0n4X6ZwTfgmYEbW8/cCL0TEI8l6bwF2BsYAZwNr+3hYOyYxNgKtKZb/D+DtQDPwNmA88PkCy58EvAb8FPgNcCpARCxJ4v1D8r5tkyz/HTLH9C/AocnyZ2Rt753AIjLH+yPgZuAdSSwnA9+VNLp3EJL2B24APg1sAxwCLE1xvCTx/glYAbw7KSr2PuyYHMd44Czge5K2TV77AfDR5BvFZOD/5djfJv9XwPXJMXYf09Rk+3emPQ7rGyf++vGLpDb5kqRf9HMbJwN3RsSdEdEVEXcDbcAxkiaSSVT/HhFvRMT9wC8LbOtHwPuTmjPAh5MygHVkEuDbIqIzIuZHxCt9jLUL+EISS8GThiQB/xf4RES8GBGvAl8F/rXAaqcBP46IziTuGZIa8mx/KJma9Wci4tWIWAp8Ezgla7F/RMR1yfZ+TOak96Uk/t8Cb5JJxL2dBVwbEXcnf5NnIuKJQsebw7PAdinfh3VJXOsi4k4yJ7/ds16bJGnriFiVnMTTuI3Mt6bdkuenkHlv3+zjcVhKTvz148SI2Ca5ndjPbTQCH8o6gbwEHAzsBLwVWBURq7OWX5ZvQxHxN2AJcHyS/N/PxsR/I5la9M2SnpV0eb6kWkB7RLyectlxwEhgftZx/Top34SknYHDyHxLgkziGgEcm2f7Y4Hh9Hw/lpGp1XZ7PuvxWoCI6F22SY2fzAni73n2m9Z4Mtcq0rwPHRGxPuv5mqy4TgKOAZYlTX7vSrPziHgD+AlwctL8NYPM/4CVSKkuelntyTVMa++yp4EbI+L/9l5QUiOwraRRWcl/Yp7tdutu7hkCPJ6cDIiIdcAXgS9KaiLzlf8v9O0iZO/9riaT1Lrj3THrtRfIJNa9IuKZFNs+JYn5l5lKMpBJ/KcCv8ix7xfI1IYbgceTsolAmn0V8zSwa39XlvQOMon/Qfr+PvQQEQ8DJyQn6XPJJPOdcy2ao+x6Msn+QWBNRPyhr/u39Fzjt27Pk2l/LlQ2m0wN/b2ShiYXYadLmhARy8g0+3xR0nBJBwPHF9nnzcBRwEw21vaRdJikvZMmklfIJM3OzTo6WAjsJalZ0gjg0u4XIqILuAa4QtL2SQzjJb03z7ZOJXNias66nQQcK2kMmfdtgqThyfY7ySTByyRtlZwkLyTzfm6uHwBnSDpC0pAk7j2KrSRpa2Uurt8MzI6IR/vxPmRvb7ikFklvSU7cr5D/b7bJ/1qS6LvINIG5tl9iTvzW7Qdk2mezrwF8DfhcUvapiHgaOAH4LNBOprb5aTb+H32YzEXKF4EvkLnomFdEPAf8ATiQTLt2tx2BW8gkjyXAfWxmkoyIvwJfAu4BniRTs8x2EfA34I+SXkmW273XMkg6AGgCvhcR/8y63Z6sP4PMRc3FwD8lvZCseh6Zbx1PJfv+EXDt5hxTclx/InOR+ArgZTLvVWOBVX4p6VUyf7tLgP+i50XmVO9DHqcAS5P1zibrgm0vPf6vsspvAPZmYE6IVoA8EYuZVQNJpwKtEXFwpWMZ7FzjN7OKSy7wnwNcXelY6oETv5lVVHINoZ1M2/+PiixuA8BNPWZmdcY1fjOzOlMT/fjHjh0bTU1NlQ7DzKymzJ8//4WI2OSHiDWR+Juammhra6t0GGZmNUVSzl/Pu6nHzKzOOPGbmdUZJ34zszpTE238Zlad1q1bx4oVK3j99bQDoVopjBgxggkTJtDQkG4QWyd+M+u3FStWsNVWW9HU1ETWSKVWRhFBR0cHK1asYJdddkm1zuBt6pE2vZnZgHr99dcZM2aMk34FSWLMmDF9+tY1OBN/vn9C/3OaDTgn/crr69+gZIk/Gav9T8pMyL1Y0heT8kuTiZgXJLdjShWDmZltqpQ1/jeAwyNiKpmJKo5OxjIHuCIimpNbeSdUnjOn+DJmVjNGj841I2VPH/nIR3j88czkZ1/96ld7vHbggQf2ex+SOOWUjVMnr1+/nnHjxnHccccV3Wa2pqYmXnjhhc1eJq2SJf7IeC152pDcKj8i3Mn55oYws1KbMweammDIkMx9ueph3//+95k0aRKwaeL//e9/3+/tjho1iscee4y1a9cCcPfddzN+/Pgia1VeSdv4k+n5FgArgbsj4qHkpXMlLZJ0raRt86zbKqlNUlt7e3spwzSzMpgzB1pbYdkyiMjct7YOXPKfN28e06dP54Mf/CB77LEHLS0tdI8+PH36dNra2rj44otZu3Ytzc3NtLS0ABtr86+99hpHHHEE++67L3vvvTe33XZbqv2+733v44477gDgpptuYsaMGRtee/HFFznxxBOZMmUKBxxwAIsWLQKgo6ODo446in322YePfvSjZI+SPHv2bPbff3+am5v56Ec/Smfn5s46mkNElPwGbAPMBSYDOwBDyZx0LgOuLbb+fvvtF32S+b/KfzOzAfH4449veHzBBRGHHpr/tsUWuT+OW2yRf50LLigew6hRoyIiYu7cubH11lvH008/HZ2dnXHAAQfEAw88EBERhx56aDz88MM9lu+9/rp16+Lll1+OiIj29vbYddddo6urK+c62esuXLgwTjrppFi7dm1MnTo15s6dG8cee2xERJx77rlx6aWXRkTEvffeG1OnTo2IiPPOOy+++MUvRkTEr371qwCivb09Hn/88TjuuOPizTffjIiImTNnxvXXXx8REY2NjdHe3p73fcj+W3QD2iJHTi1LP/6IeEnSPODoiPhGd7mka4BflSMGM6usN97oW3l/7L///kyYMAGA5uZmli5dysEHp5vJMSL47Gc/y/3338+QIUN45plneP7559lxxx0LrjdlyhSWLl3KTTfdxDHH9Oyr8uCDD/Kzn/0MgMMPP5yOjg5efvll7r//fm699VYAjj32WLbdNtPwce+99zJ//nze8Y53ALB27Vq233779G9ASiVL/JLGAeuSpL8lcCTwH5J2iswk2wAfAB4b8J1HuOumWZl961uFX29qyjTv9NbYCPPmDUwMW2yxxYbHQ4cOZf369anXnTNnDu3t7cyfP5+GhgaamppS941///vfz6c+9SnmzZtHR0fHhvLIMdFVd9fLXF0wI4LTTjuNr33ta6nj7o9StvHvBMyVtAh4mEwb/6+AyyU9mpQfBnyihDGYWZW47DIYObJn2ciRmfJyamhoYN26dZuUv/zyy2y//fY0NDQwd+5cluU6S+Vx5pln8vnPf5699967R/khhxzCnOQixrx58xg7dixbb711j/K77rqLVatWAXDEEUdwyy23sHLlSiBzjaAvcaRVshp/RCwC9slRfkqOxctr/Hh45plKR2FWV5JrqVxyCSxfDhMnZpJ+d3m5tLa2MmXKFPbdd98NyTcTXwvHH38806ZNo7m5mT322CP1NidMmMAFF1ywSfmll17KGWecwZQpUxg5ciTXX389AF/4wheYMWMG++67L4ceeigTJ04EYNKkSXzlK1/hqKOOoquri4aGBr73ve/R2Ni4mUfdU03MuTtt2rTo80QsxZp6auC4zardkiVL2HPPPSsdhpH7byFpfkRM673s4ByywczM8hq8if+tb610BGZmVWnwJn634ZuZ5TR4E7+ZmeVUv4l/r70qHYGZWUXUb+JPRuozM6s39Zv4zazmdXR00NzcTHNzMzvuuCPjx4/f8PzNN98ckH1Mnz6d3XffnSlTprDHHntw7rnn8tJLLw3Itru1tbVx/vnnA5kfem3OiKFpDO7EnwzDamZVYoDHZR4zZgwLFixgwYIFnH322XziE5/Y8Hz48OF9GrKhcNhzWLRoEYsWLWKLLbbghBNOGJDtdps2bRrf/va3ASf+zbd4caUjMLNupR6XOXH66adz4YUXcthhh3HRRRdx6aWX8o1vbBgbksmTJ7N06VKg70MgDx8+nMsvv5zly5ezcOHCgtsYPXo0l1xyCVOnTuWAAw7g+eefB+CnP/0pkydPZurUqRxyyCFAJtkfd9xxLF26lKuuuoorrriC5uZmHnjgAXbZZZcNQ0y88sorNDU15Rxyoi8Gd+IvxrNxmQ2cj38cpk/PfzvrLFizpuc6a9ZkyvOt8/GP9yuUv/71r9xzzz1885vfzLvMkiVL+PGPf8zvfvc7FixYwNChQ3sM4ZDP0KFDmTp1Kk888UTBbaxevZoDDjiAhQsXcsghh3DNNdcA8KUvfYnf/OY3LFy4kNtvv73Htpuamnp8c3n3u9/N9OnTN4z3f/PNN3PSSSfR0NDQr/elW30n/jPPrHQEZvWjHOMyJz70oQ8xdOjQgstkD4Hc3NzMvffey1NPPZVq+91D3RTaxvDhwzdMwbjffvtt+JZx0EEHcfrpp3PNNdekmmTlIx/5CNdddx0A1113HWeccUaqGAspy3j8VWuALv6YGdUxLnNi1KhRGx4PGzaMrq6uDc+7h1ru7xDInZ2dPProo+y5556sXLky7zYaGho2DL2cPUT0VVddxUMPPcQdd9xBc3MzCxYsKLi/gw46iKVLl3LffffR2dnJ5MmT+xRvLvVd4zez8qnQuMxNTU088sgjADzyyCP84x//APo3BPK6dev4zGc+w84778yUKVP6tY2///3vvPOd7+RLX/oSY8eO5emnn+7x+lZbbcWrr77ao+zUU09lxowZA1Lbh3pI/DNnVjoCM4PM+MtXX52p4UuZ+6uvLvm4zCeddBIvvvgizc3NzJo1i7e//e1AzyGQp0yZwnve8x6ee+65nNtoaWlhypQpTJ48mdWrV2+Yj7cv2+j26U9/mr333pvJkydzyCGHMHXq1B6vH3/88fz85z/fcHG3e/+rVq3qMZ/v5hi8wzJnKzREcw0cv1m18rDM5XHLLbdw2223ceONN+Zdpi/DMtd3Gz/AOefAlVdWOgozs5zOO+887rrrLu68884B26YT/6xZTvxmVrW+853vDPg2B38bv5mVVC00Fw92ff0blCzxSxoh6U+SFkpaLOmLSfl2ku6W9GRyv22pYjCz0hoxYgQdHR1O/hUUEXR0dDBixIjU65SyqecN4PCIeE1SA/CgpLuA/wXcGxFfl3QxcDFwUQnjyPTsmTWrpLswq0cTJkxgxYoVtLe3VzqUujZixAgmTJiQevmy9OqRNBJ4EJgJ3ABMj4jnJO0EzIuI3Qutv9m9ejJB5H/NtRUzG4QqMtm6pKGSFgArgbsj4iFgh4h4DiC53z7Puq2S2iS1lbw2ceSRpd2+mVkVKWnij4jOiGgGJgD7S0r9W+OIuDoipkXEtHHjxpUsRgDuvbe02zczqyJl6dUTES8B84CjgeeTJh6S+5XliMHMzDJK2atnnKRtksdbAkcCTwC3A6cli50G3FaqGHrwpCxmZkBpe/XsBFwvaSiZE8xPIuJXkv4A/ETSWcBy4EMljGGjxYsLX+A1M6sTJUv8EbEI2CdHeQdwRKn2228eusHM6oR/udvN/fzNrE448ZuZ1Zn6SvxHVF8Lk5lZudVX4r/nnkpHYGZWcfWV+Is555xKR2BmVnJO/Nl8gdfM6oATv5lZnam/xO8LvGZW5+ov8fsCr5nVufpL/MV4iGYzG+Sc+HvzEM1mNsg58ZuZ1Zn6TPwzZ1Y6AjOziqnPxO9ROM2sjtVn4i/GF3jNbBBz4s/FF3jNbBCr38Tv2bjMrE7Vb+I/++xKR2BmVhH1m/iLXeD1SJ1mNkiVLPFL2lnSXElLJC2WdEFSfqmkZyQtSG7HlCqGzeKROs1skCrZZOvAeuCTEfGIpK2A+ZLuTl67IiK+UcJ9m5lZHiVL/BHxHPBc8vhVSUuA8aXaX7+MHg2vvVbpKMzMyqosbfySmoB9gIeSonMlLZJ0raRt86zTKqlNUlt7e3tpArvqqtJs18ysipU88UsaDfwM+HhEvALMAnYFmsl8I/hmrvUi4uqImBYR08aNG1ea4FpaCr/uH3KZ2SBU0sQvqYFM0p8TEbcCRMTzEdEZEV3ANcD+pYxhs/iHXGY2CJWyV4+AHwBLIuK/ssp3ylrsA8BjpYrBzMw2Vcoa/0HAKcDhvbpuXi7pUUmLgMOAT5QwhuI8UqeZ1RlFROEFpHPJNNWsKk9Im5o2bVq0tbWVbgeFhm+YPbv4tQAzsyokaX5ETOtdnqbGvyPwsKSfSDo6acKpH2eeWekIzMwGVNHEHxGfA3Yj015/OvCkpK9K2rXEsVWHN9+sdARmZgMqVRt/ZNqD/pnc1gPbArdIuryEsZXPiBGVjsDMrGyKJn5J50uaD1wO/A7YOyJmAvsBJ5U4vvL4/vcLvz5nTnniMDMrgzQ1/rHA/4qI90bETyNiHUDSD/+4kkZXLsUu3p52WnniMDMrgzSJf5eIWJZdIOlGgIhYUpKoqk1nZ6UjMDMbMGkS/17ZTyQNJdPMM7i4nd/M6kTexC/pM5JeBaZIeiW5vQqsBG4rW4Tl4nZ+M6sTaX7A9bWI+EyZ4smp5D/g6lboJwqjRnkIZzOrKfl+wJV3PH5Je0TEE8BPJe3b+/WIeGSAY6xuq1dXOgIzswFRaCKWC4FWcg+bHMDhJYmokiQo8g3IzKzW5W3jj4jW5P6wHLfBl/QBzj678OuegN3MBoE0P+D6UDJnLpI+J+lWSfuUPrQKuPLKwq97AnYzGwTSdOf892TO3IOB9wLXA56z0MysRqVJ/N2/XjoWmBURtwHDSxdShY0eXekIzMxKKk3if0bS/wD/G7hT0hYp16tNxSZgdzu/mdW4NP34RwJHA49GxJPJ1Il7R8RvyxEglLEff7diUw6454+Z1YB+T8QSEWvI/FJ3taSJQAPwxMCHaGZm5VCoHz8Aks4DvgA8D3QlxQFMKWFclTViBLz+eqWjMDMriTRt9RcAu0fEXhGxd3IrmvQl7SxprqQlkhZLuiAp307S3ZKeTO633dyDGHDFxu3Za6/Cr5uZVbE0if9p4OV+bHs98MmI2BM4APiYpEnAxcC9EbEbcG/yvLoUG5//8cfLE4eZWQkUbeoBngLmSboDeKO7MCL+q9BKEfEc8Fzy+FVJS4DxwAnA9GSx64F5wEV9DdzMzPonTY1/OXA3mb77W2XdUpPUBOwDPATskJwUuk8O2+dZp1VSm6S29vb2vuxuYMycWfh1d+s0sxpVtDvnhgWlURHR5yEqJY0G7gMui4hbJb0UEdtkvb4qIgq285e9O2c3d+s0sxrW7+6ckt4l6XFgSfJ8qqQig9psWLcB+BkwJyJuTYqfT34LQHK/MuUx9IlU/GZmVo/SNPV8i8wYPR0AEbEQOKTYSpIE/ABY0ut6wO1A9+zlp1GC2bzSJvWiJwBPx2hmg1CqoRci4uleRWlmHz8IOAU4XNKC5HYM8HXgPZKeBN6TPK+ovCeAYt06jzyyJPGYmZVSml49T0s6EAhJw4HzSZp9ComIB4F89ekj0odYPlLmmu6G0ZlbWuDkk/OvcO+9ZYnLzGwgpanxnw18jExXzBVAc/J8UJo1C0aOrHQUZmalk2asnhcioiUidoiI7SPi5IjoKEdwlbJ2bVbTzxFFvpy4W6eZ1ZiCiV/SYcmMW4uT2y2SppcntP4bqF6WEnDPPYUX8qxcZlZj8iZ+SccC1wK/BD4MtAB3AtcmF2mrWkTP2zbb9G877vZpZoNNoRr/p4ETI+K6iFgYEQsi4lrgRGpwiIVVqzaeBPrqFYrMyuXmHjOrIYUS/45Jn/0eImIRsEPpQiq9iOJN99nO5ioKni/c3GNmNaRQ4i80PEOfh26oNvfck772fxNFRus0M6shhfrx7yrp9hzlAv6lRPGUXUS6dvxXGc3WvJZ/gXPOyfoBgJlZ9co7SJukQwutGBH3lSSiHMoxSFux5D+DOczh5Ly/SAM8aJuZVZV8g7SlHp2zkso1Omex5N+FnPjNrGb0e3TOelIsb79OkUHb3LvHzGqAE38vhZL/WXzfvXvMrOY58eeQL/m7d4+ZDQZpJmK5W9I2Wc+3lfSbkkZVBfIlfzf3mFmtS1PjHxsRL3U/iYhV5Jknd7DJNe2um3vMrNalSfxdkiZ2P5HUCIVz32Bx5ZXQ0NCzzM09Zlbr0iT+S4AHJd0o6UbgfuAzpQ2rerz55qZl64rNX+OZucysiqUZj//XwL7Aj4GfAPtFxKBv4882e3bP56fzw8JfeTwzl5lVsULDMu+R3O8LTASeBZ4BJiZldaOlBSZN2vg8VXPPnDmlC8jMbDMUqvFfmNx/M8ftG8U2LOlaSSslPZZVdqmkZ3pNvl4TFi+GIVnv1m85onCt/8wzSx2SmVm/5G2sjojW5P6wfm77h8B3gRt6lV8REUVPHNWos3PjsA5Hcw9dhQZwyHVxwMysCqTpxz9U0vslnS/pwu5bsfUi4n7gxQGJsork6uKZl5t7zKwKpenV80vgdGAMsFXWrb/OlbQoaQraNt9CkloltUlqa29v34zdDawrr9zY5PM9ZhZu7jnttHKEZGbWJ0VH55S0KCKm9GvjUhPwq4iYnDzfAXiBzO8AvgzsFBFFG8PLNTpnX3Q3+RQasTMAecROM6uQzRmd8y5JRw1EEBHxfER0RkQXcA2w/0BstxK6p24smtbdp9/MqkyaxP9H4OeS1kp6RdKrkl7pz84k7ZT19APAY/mWrXb33JO5v7JAc4+AcJ9+M6syaRL/N4F3ASMjYuuI2Coiti62kqSbgD8Au0taIeks4HJJj0paBBwGfGJzgq+0CDiPFNMtutZvZlWkyNgDADwJPBZ9nKorImbkKP5BX7ZRCyZNgkWPT2IKj+ds6++u9aeY1tfMrCzSJP7ngHmS7gLe6C6MiP8qWVQ1ZPFikBYX7tNPpmdni8d3M7MqkKap5x/AvcBwBqY756Azc2bxi7yTT96rLLGYmRWTerJ1SVsBERGvlTakTVVjd87erhp6Dh/tmlWwa+ePZodr/WZWNv3uzilpsqQ/k+mBs1jSfEmuvvZydmfxi7wPXeBf8ppZ5aVp6rkauDAiGiOiEfgkmT741kvbNvkHbhNwYcclHsXBzCouzcXdURExt/tJRMyTNKqEMdWsd6y6h1D+i7wTWc5erZnHbvIxs0pJU+N/StK/S2pKbp8jc8HXctDM/D/o6mIIJ6yZwyWXlDUkM7Me0iT+M4FxwK3JbSyZQdsslyuvREfkbvIZRifX0MpBy9zeY2aVkybxHxkR50fEvsnt48B7ShxXbbvnHjR7NusZuslLo1jDV3Fbv5lVTprEn2ti9bqZbL3fWloYSlfOl3ZmOa2tHq7fzCoj78VdSe8DjgHGS/p21ktbA+tLHdhgoMaJsGzZJuXLmciaNXDJJb7Ia2blV6jG/yzQBrwOzM+63Q4MyDDNg95ll8HIkT2KAhjFa8xgDsuXVyYsM6tvhebcXQgslPSjiFjXXS7pYOAy4GNliK+2JdX5jlMuYLvoQGT684+jg2toZQgwZ06La/1mVlZF2/gjYp2kZkmXS1pKZuasJ0oe2WDR0sIW243eZCiHUazhK3GJ2/rNrOzyJn5Jb5f0eUlLgO8CT5MZ2+ewiPhO2SIcBEa/mLtNZyLLN7T1m5mVS6Ea/xPAEcDxEXFwkuw7yxPWIDNxYs7iLoYwgzksW+Zav5mVT6HEfxLwT2CupGskHQGeT6RfclzkhY0/6JrBHDf5mFnZ5E38EfHziPg/wB7APDLTJO4gadZATb5eN1pa4Oqr6RqS/wddbvIxs3JJc3F3dUTMiYjjgAnAAuDiUgc26LS0MCRy/6BrIplrAO7eaWblkOaXuxtExIsR8T8RcXixZSVdK2mlpMeyyraTdLekJ5P7bfsTdM0q0tY/ZIibe8ys9PqU+Pvoh8DRvcouBu6NiN3ITOdYX98cirT1/+9Ot/WbWemVLPFHxP3Ai72KTwCuTx5fD5xYqv1XpaStn6Fu6zezyilljT+XHSLiOYDkfvt8C0pqldQmqa29vb1sAZZcSwt05W7rb2SZh3Iws5Ird+JPLSKujohpETFt3LhxlQ5nYOVp6xdwDa18mDlu7jGzkil34n9e0k4Ayf3KMu+/OuRp6wcP5WBmpVfuxH87cFry+DTgtjLvvzp0t/Xn4aEczKyUSpb4Jd0E/AHYXdIKSWcBXwfeI+lJMrN4fb1U+696LS3Q2JjzpeyhHJqaXPM3s4GVd1jmzRURM/K8dESp9llzLrsMWlthzZoexd3dOwFuWtZCa+ahh282swFRtRd360KKoRwAN/uY2YBy4q+0AkM5NLKMf9DkLp5mNqCc+KtBge6dTSxzF08zG1BO/NWgQPdOcBdPMxtYTvzVoLt7Z55ePuAunmY2cJz4q0VLCyxd6i6eZlZyTvzVJsVsXcuW4WYfM+s3J/5q4y6eZlZiTvzVyF08zayEnPirlbt4mlmJOPFXK3fxNLMSceKvVim6eDayjBPWzHFbv5n1iRN/NSvSxbN74pYDl81xF08zS82JvxYUmbjlq1ziLp5mlpoTfy1Imn0iz8vdPX3c7GNmaTjx14qWFlSgyae7p8+By9zTx8wKc+KvJSl6+tzAadxzhpO/meXnxF9Lsnr65Gv2GUYn313XykMXOPObWW5O/LUm6emTr9kHMjX/CzsucU8fM8upIolf0lJJj0paIKmtEjHUvCLNPo0sY96yJjf7mNkmKlnjPywimiNiWgVjqF0FBnODjRd83exjZr25qaeWtbQw5IbrWT+88AXf/+44mRXDmnjwHJ8AzKxyiT+A30qaL6k11wKSWiW1SWprb28vc3g1pKWFYdcWvuArYELnMvaZ1erkb2YoIl+6KOFOpbdGxLOStgfuBs6LiPvzLT9t2rRoa/OlgGJeG9vE6I5lBZcJ4JmhjSxtvYyDr2wpT2BmVhGS5udqTq9IjT8ink3uVwI/B/avRByDzej/vqxgsw+49m9mFUj8kkZJ2qr7MXAU8Fi54xiUUjT7dBvFGg6a5bZ/s3pUiRr/DsCDkhYCfwLuiIhfVyCOwSnp5/+7mbNZjWv/Zrapsif+iHgqIqYmt70i4rJyx1APDr6yhT/PvJoVQ137N7Oe3J1zEDv4yhYmrHft38x6cuKvA679m1k2J/460Z/a/4GzTiEkPOiP2eDixF9n+lL7H0IggGXLWH+mp/cyGyyc+OtQX2r/3Ya9uYY4+WTX/s0GASf+OtaX2j+wofYfJ58CbgIyq1lO/HWuP7V/dZ8m3ARkVpOc+A3oWfvvAroy9fuiupuA3AvIrHY48dsG3bX/m2YHZzXcyFLSNwFN6FzGtFln0DFkLF0a4hOBWRVz4rdNtLTAkde1ML1xKS2kbwIawTrGRAdDiA3dQbsknwTMqowTv+WUDPnDsbNbOLfhapbStyYgyHQHHYK/DZhVGyd+Kyi79j+U4BTSNwFl87cBs+rhxG9Fddf+I+CY2X1vAsol37eBDo31NwOzEnPitz7pPgn8KLJ7AYkOxvAGw/u1zexvA2Po6PHNwE1EZgPPid/6rbsX0JDoYky8wMMzr+1zd9BiejcR9T4RzNvrHFYMa/KJwawPnPhtwGw8EQS/n3njgHwb6K33ieDQx2cxoXNZzmsHvU8K1x45h6YmGDLEPzq2+laRydb7ypOt174Hz5lD09WX8NbOZYAY0ufLw30X0ON7x+s08CpbM4YXWc5E7uQYjh9yJ+O7lvPs0Iksbc3MCZSJc2OZJ6W3WpVvsnUnfiu7jSeB5axiO0bzKlvwZtnjyHViEOoRy+s0sFpbs228yLNDJ/K33Y/hbX+5c8OJofdznzysmjjxW9XqfSJAsG28WNGTQj69TxYDcfLId7LoXdbf9Xwyql/5Ej8RUfYbcDTwF+BvwMXFlt9vv/3C6tMDM2fH00MboxPFC4yJ1xkekelZGgHRlfW4Vm69Y15LwybHlausv+vlK3tBY6ITxdNDG+OBmbN7vNd9Lev9tyr3Nsu9bCXL+wJoi1w5OFdhKW/AUODvwL8Aw4GFwKRC6zjxW7feH4a5k2Ymz4lOVPMnhUrd3mBYvEFD0bLXadikbDVbxn27nRmr2bJH+eacmF5jZMydNDNeY2TVLVvJ8r4m/3yJv+xNPZLeBVwaEe9Nnn8m+ebxtXzruKnH0shuMurdNLJK27G1XqWha2MTTFC8mcYqp/ffp5qWrVT5iqGNTFi/NGesueRr6qlEd87xwNNZz1ckZT1IapXUJqmtvb29bMFZ7cr+XcGE9UuZvvjKjb8z6HqBhhuuhcbGzCQyjY385YiZG7qcrhjaSNvM67J+i5C7G2rvalLv56/TMGBdV2tJeauP9eutncsHZDuVqPF/CHhvRHwkeX4KsH9EnJdvHdf4rVIKfYtIc8E11wXqNN80cpX1d71yfJNZz1CG0VmxbZZ72UqV13KNfwWwc9bzCcCzFYjDrKhC3yJyPT/4ypYCv2jOfLO4b1Lhbxr5yvq7XppvMrm+qaQtW81IHpzUusnYTeXaZrmXrWR5d8Vis+Vq+C/lDRgGPAXswsaLu3sVWscXd80GViV74FRrT51q671Tyl49FenHL+kY4FtkevhcGxEFT2Nu6jEz67t8TT3DKhFMRNwJ3FmJfZuZ1TsP0mZmVmec+M3M6owTv5lZnXHiNzOrMzUxOqekdmBZP1cfC7wwgOFUio+jegyGYwAfR7UpxXE0RsS43oU1kfg3h6S2XN2Zao2Po3oMhmMAH0e1KedxuKnHzKzOOPGbmdWZekj8V1c6gAHi46geg+EYwMdRbcp2HIO+jd/MzHqqhxq/mZllceI3M6szgybxSzpa0l8k/U3SxTlel6RvJ68vkrRvJeIsJsVxtCTxL5L0e0lTKxFnIcWOIWu5d0jqlPTBcsaXVprjkDRd0gJJiyXdV+4Y00jxP/UWSb+UtDA5jjMqEWchkq6VtFLSY3ler5XPd7HjKM/nO9dYzbV2I8UE7sAxwF1kJjE6AHio0nH38zgOBLZNHr+v2o4jzTFkLff/yIzS+sFKx93Pv8U2wOPAxOT59pWOu5/H8VngP5LH44AXgeGVjr1XjIcA+wKP5Xm96j/fKY+jLJ/vwVLj3x/4W0Q8FRFvAjcDJ/Ra5gTghsj4I7CNpJ3KHWgRRY8jIn4fEauSp38kM4NZNUnztwA4D/gZsLKcwfVBmuP4MHBrRCwHiIhqPJY0xxHAVpIEjCaT+NeXN8zCIuJ+MnHlUwuf76LHUa7P92BJ/GkmcE81yXuF9TXGs8jUcqpJ0WOQNB74AHBVGePqqzR/i7cD20qaJ2m+pFPLFl16aY7ju8CeZKZAfRS4ICK6yhPegKmFz3dflezzXZGJWEpAOcp691NNs0ylpY5R0mFk/jEOLmlEfZfmGL4FXBQRnZlKZlVKcxzDgP2AI4AtgT9I+mNE/LXUwfVBmuN4L7AAOBzYFbhb0gMR8UqJYxtItfD5Tq3Un+/BkvjTTOBeC5O8p4pR0hTg+8D7IqKjTLGlleYYpgE3J0l/LHCMpPUR8YuyRJhO2v+pFyJiNbBa0v3AVKCaEn+a4zgD+HpkGpb/JukfwB7An8oT4oCohc93KuX4fA+Wpp6Hgd0k7SJpOPCvwO29lrkdODW5+n8A8HJEPFfuQIsoehySJgK3AqdUWc2yW9FjiIhdIqIpIpqAW4BzqizpQ7r/qduAd0saJmkk8E5gSZnjLCbNcSwn860FSTsAuwNPlTXKzVcLn++iyvX5HhQ1/ohYL+lc4DdsnMB9saSzk9evItN75Bjgb8AaMrWcqpLyOD4PjAGuTGrM66OKRiZMeQxVL81xRMQSSb8GFgFdwPcjImc3vUpJ+ff4MvBDSY+SaTK5KCKqaphjSTcB04GxklYAXwAaoHY+35DqOMry+faQDWZmdWawNPWYmVlKTvxmZnXGid/MrM448ZuZ1RknfjOzOuPEbzUvGeFzgaTHklEmt+nj+vMkTUse31lsfUlLJY3NU/5oMsrlbyXtWGAb0yR9u8h+tpF0TsrDMEvNid8Gg7UR0RwRk8kMgPWx/m4oIo6JiJc2I5bDImIq0EZm1Mt8+2mLiPOLbGsbwInfBpwTvw02fyAZnEvS/smY5n9O7ndPyreUdHMy5vmPyYyzQ/Lahtq8pF8kg68tltTaxzjuB94maYSk65JvAn9OxmDpHsf/V8njS5Nx2udJekpS9wnh68CuybeZ/5S0k6T7s77dvHtz3iirX4Pil7tmAJKGkhl64AdJ0RPAIcmvV48EvgqcBMwE1kTElGRclEfybPLMiHhR0pbAw5J+1oexU44jM9LlxwAiYm9JewC/lfT2HMvvARwGbAX8RdIs4GJgckQ0J8f3SeA3EXFZcqwjU8Zi1oMTvw0GW0paADQB84G7k/K3ANdL2o3MSI0NSfkhwLcBImKRpEV5tnu+pA8kj3cGdgOKJf65kjrJDOPwOeA64DvJvp6QtIzMcM693RERbwBvSFoJ7JBjmYeBayU1AL+IiAVFYjHLyU09NhisTWrFjWRmmepu4/8yMDdp+z8eGJG1TsGxSiRNB44E3pW02f+51/r5HJZcbzg1uVaQdtzpN7Ied5KjUpZM4nEI8Axwo6pz/H+rAU78NmhExMvA+cCnklrxW8gkSYDTsxa9H2gBkDQZmJJjc28BVkXEmqSJ5oB+hpW9r7cDE4G/pFz3VTJNPyTrNwIrI+IaMs1ZVTmvrFU/J34bVCLiz2Tmlf1X4HLga5J+R2Zkym6zgNFJE8+/kXvc+V8Dw5JlvkxmGrz+uBIYmox8+WPg9KRJJ82xdAC/Sy7k/ieZUR0XSPozmWsV/93PmKzOeXROM7M64xq/mVmdceI3M6szTvxmZnXGid/MrM448ZuZ1RknfjOzOuPEb2ZWZ/4/t4xQbxLL7vMAAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "max_numb_functions = 4\n",
    "factor = 50.0\n",
    "result = greedy.run(factor, max_numb_funcs=max_numb_functions, disp=True)\n",
    "\n",
    "from bfit.model import AtomicGaussianDensity\n",
    "model = AtomicGaussianDensity(grid.points, num_s = result[\"num_s\"], normalize=True)\n",
    "print(result[\"coeffs\"])\n",
    "model0 = model.evaluate(result[\"coeffs\"], result[\"exps\"])\n",
    "\n",
    "numb_pts = 500\n",
    "plt.plot(grid.points[:numb_pts], model0[:numb_pts], \"bo-\", label=\"Initial Model\")\n",
    "plt.plot(grid.points[:numb_pts], density[:numb_pts], \"ro-\", label=\"True Density\")\n",
    "plt.xlabel(\"Radial Points\")\n",
    "plt.ylabel(\"Atomic Density\")\n",
    "plt.title(\"Fitted vs True Atomic Density\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "It becomes apparent that more basis-functions are needed in order to capture the cusp better. Free-free to re-run and modify the calculation as needed."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "language": "python",
   "display_name": "Python 3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}