tutorials/21-compare-radial-atomic-potentials.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute and compare different radial potentials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using JAC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We shall here describe how different atomic (radial) potentials can be calculated and applied within the JAC toolbox.\n",
    "\n",
    "Any (radial) potential is handled in JAC by the data struct `Radial.Potential` that (may) contain a name, the radial function $Z(r)$ as well as the radial grid on which the potential is defined. The potential function `Zr` typically contains the effective charge $Z(r)$ as seen by an electron, so that the full radial potential is given by $V = -\\:Z(r) / r$. This convention has been frequently applied in atomic physics and has also been adopted in JAC.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\texttt{struct  Radial.Potential}  ... defines a struct for a local radial potential.\n",
       "\n",
       "\\begin{verbatim}\n",
       "+ name           ::String            ... A name for the potential.\n",
       "+ Zr             ::Array{Float64,1}  ... radial potential function Z(r) = - r * V(r) as usual in atomic structure theory.\n",
       "+ grid           ::RadialGrid        ... radial grid on which the potential is defined.\n",
       "\\end{verbatim}\n",
       "\\rule{\\textwidth}{1pt}\n",
       "\\texttt{Radial.Potential()}  ... constructor to define an 'empty' instance of the radial potential.\n",
       "\n"
      ],
      "text/markdown": [
       "`struct  Radial.Potential`  ... defines a struct for a local radial potential.\n",
       "\n",
       "```\n",
       "+ name           ::String            ... A name for the potential.\n",
       "+ Zr             ::Array{Float64,1}  ... radial potential function Z(r) = - r * V(r) as usual in atomic structure theory.\n",
       "+ grid           ::RadialGrid        ... radial grid on which the potential is defined.\n",
       "```\n",
       "\n",
       "---\n",
       "\n",
       "`Radial.Potential()`  ... constructor to define an 'empty' instance of the radial potential.\n"
      ],
      "text/plain": [
       "  \u001b[36mstruct Radial.Potential\u001b[39m ... defines a struct for a local radial potential.\n",
       "\n",
       "\u001b[36m  + name           ::String            ... A name for the potential.\u001b[39m\n",
       "\u001b[36m  + Zr             ::Array{Float64,1}  ... radial potential function Z(r) = - r * V(r) as usual in atomic structure theory.\u001b[39m\n",
       "\u001b[36m  + grid           ::RadialGrid        ... radial grid on which the potential is defined.\u001b[39m\n",
       "\n",
       "  ────────────────────────────────────────────────────────────────────────────\n",
       "\n",
       "  \u001b[36mRadial.Potential()\u001b[39m ... constructor to define an 'empty' instance of the\n",
       "  radial potential."
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "? Radial.Potential"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Perhaps, the simplest choice of a radial potential is the **nuclear** potential which, of course, depends on the nuclear model. Let us first simply consider a $\\;^{12}$C $\\;$ carbon atom with a Fermi-distributed nuclear charge and a standard radial grid.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Fermi nuclear model for Z = 6.0 with mass = 12.18, radius R = 2.4934845703596404 fm and nuclear spin I = 0, dipole moment mu = 0.0 and quadrupole moment Q = 0.0."
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid     = Radial.Grid(true)\n",
    "nucModel = Nuclear.Model(6., \"Fermi\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Then, the corresponding nuclear potential is generated by simply calling:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "nuclear-potential: Fermi-distributed (radial) potential ... defined on 392 grid points ...\n",
       "Zr:    [0.0, 0.015619818207380716, 0.035905804451892086, 0.060430845202356614, 0.08495453726302814, 0.10523736995823374, 0.11777879390742439, 0.12582486891780031, 0.14610064533238984, 0.17888509371708336  …  0.31112937694872056, 0.32313641267043025, 0.35185818903226784, 0.39828287685804126, 0.4543674906380097, 0.5103989493317181, 0.5566993479831148, 0.5853078885893802, 0.6022668932973185, 0.6428176395098482]  ...  [6.0, 6.000000000000002, 6.000000000000001, 6.000000000000002, 6.000000000000001, 6.0, 6.000000000000001, 6.0, 6.000000000000002, 6.000000000000002  …  6.000000000000002, 6.000000000000001, 6.000000000000002, 6.000000000000002, 6.000000000000001, 6.000000000000001, 6.0, 6.000000000000001, 6.000000000000001, 6.000000000000001]\n",
       "Radial grid:  rnt = 2.0e-6,  h = 0.05,  hp = 0.0,  NoPoints = 392,  ntL = 69,  ntS = 71, orderL = 7,  orderS = 8,  nsL = 62,  nsS = 63,  mesh = JAC.Radial.MeshGL(), ...  \n",
       "r:    [1.780504510281598e-8, 9.042759120177923e-8, 2.0787030685460475e-7]  ...  [564.4302351890452, 595.097239730825, 614.0606632064549]\n",
       "wr:   [4.530145586285996e-8, 9.785739473901763e-8, 1.3358660616640612e-7]  ...  [34.882547080539766, 25.552825070730517, 11.829256033233948]\n",
       "tS:   [0.0, 0.0, 0.0]  ...  [618.7099715607404, 618.7099715607404, 618.7099715607404]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nucPotential = Nuclear.nuclearPotential(nucModel, grid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "As seen from the output, the effective charge Z(r) vanishes near to the origin but quickly increases to the full charge $Z = 6.\\;$ at larger $r$-values, i.e. near to the surface and outside of the nucleus\n",
    "\n",
    "To calculate **atomic** potentials, we also need to know the electron density around the nucleus. Most often, this density refers to a particular atomic level. Therefore, we first need to compute some of the low-lying levels of, for instance, neutral carbon. We shall here just compute the levels of the $1s^2 2s^2 2p^2$ ground-state configuration. We use the *simplest* approximation for these computation but will explain in some later tutorials how we could control and refine such simple structure calculations, if needed.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "... in performSCF ...\n",
      ">> include Configuration: 1s_1/2^2 2s_1/2^2 2p_1/2^0 2p_3/2^2 \n",
      ">> include Configuration: 1s_1/2^2 2s_1/2^2 2p_1/2^1 2p_3/2^1 \n",
      ">> include Configuration: 1s_1/2^2 2s_1/2^2 2p_1/2^2 2p_3/2^0 \n",
      "(Re-) Define a new standard subshell list.\n",
      "Start SCF process with hydrogenic orbitals.\n",
      ">> (Re-) Define a storage array for various B-spline matrices:\n",
      "Nuclear model = Fermi nuclear model for Z = 6.0 with mass = 12.18, radius R = 2.4934845703596404 fm and nuclear spin I = 0, dipole moment mu = 0.0 and quadrupole moment Q = 0.0. \n",
      "Generate hydrogenic orbital for subshell 1s_1/2 \n",
      "  -----------------------------------------------------------------------------\n",
      "   Index    Subshell     Energies [a.u.]    Dirac-E  [a.u.]     Delta-E / |E|    \n",
      "  -----------------------------------------------------------------------------\n",
      "      1      1s_1/2      -1.80086332e+01    -1.80086350e+01    +1.02396837e-07    \n",
      "      2      2s_1/2      -4.50269966e+00    -4.50269857e+00    -2.42835410e-07    \n",
      "      3      3s_1/2      -2.00098004e+00    -2.00095940e+00    -1.03150269e-05    \n",
      "      4      4s_1/2      -1.12562938e+00    -1.12543844e+00    -1.69628100e-04    \n",
      "      5      5s_1/2      -7.21117584e-01    -7.20234830e-01    -1.22414829e-03    \n",
      "      6      6s_1/2      -5.30019154e-01    -5.00139888e-01    -5.63739359e-02    \n",
      "      7      7s_1/2      -4.90072472e-01    -3.67436826e-01    -2.50239816e-01    \n",
      "      8      8s_1/2      -3.59948993e-01    -2.81311119e-01    -2.18469492e-01    \n",
      "      9      9s_1/2      -2.64671561e-01    -2.22265641e-01    -1.60220918e-01    \n",
      "     10     10s_1/2      -1.92978410e-01    -1.80031939e-01    -6.70876634e-02    \n",
      "     11     11s_1/2      -1.39236869e-01    -1.48784504e-01    +6.85711654e-02    \n",
      "     12     12s_1/2      -9.97730204e-02    -1.25018733e-01    +2.53031453e-01    \n",
      "     13     13s_1/2      -7.12821673e-02    -1.06523685e-01    +4.94394588e-01    \n",
      "     14     14s_1/2      -5.10014594e-02    -9.18486436e-02    +8.00902262e-01    \n",
      "     15     15s_1/2      -3.67807216e-02    -8.00097188e-02    +1.17531672e+00    \n",
      "     16     16s_1/2      -2.69852768e-02    -7.03205343e-02    +1.60588524e+00    \n",
      "     17     17s_1/2      -2.03591822e-02    -6.22904546e-02    +2.05957548e+00    \n",
      "     18     18s_1/2      -1.59386307e-02    -5.55612291e-02    +2.48594744e+00    \n",
      "     19     19s_1/2      -1.30118085e-02    -4.98663309e-02    +2.83239048e+00    \n",
      "     20     20s_1/2      -1.10621734e-02    -4.50041540e-02    +3.06829222e+00    \n",
      "     21     21s_1/2      -9.87491289e-03    -4.08199215e-02    +3.13369940e+00    \n",
      "     22     22s_1/2      +3.34227908e-01    -3.71932146e-02    +1.11128100e+00    \n",
      "     23     23s_1/2      +3.69962786e+00    -3.40292102e-02    +1.00919801e+00    \n",
      "     24     24s_1/2      +1.30389838e+01    -3.12524195e-02    +1.00239684e+00    \n",
      "     25     25s_1/2      +3.55194076e+01    -2.88021434e-02    +1.00081088e+00    \n",
      "     26     26s_1/2      +8.57637919e+01    -2.66291266e-02    +1.00031049e+00    \n",
      "     27     27s_1/2      +1.93321549e+02    -2.46930634e-02    +1.00012773e+00    \n",
      "     28     28s_1/2      +4.17176371e+02    -2.29607143e-02    +1.00005504e+00    \n",
      "     29     29s_1/2      +8.72261367e+02    -2.14044706e-02    +1.00002454e+00    \n",
      "     30     30s_1/2      +1.77341765e+03    -2.00012468e-02    +1.00001128e+00    \n",
      "     31     31s_1/2      +3.49574464e+03    -1.87316200e-02    +1.00000536e+00    \n",
      "     32     32s_1/2      +6.63118506e+03    -1.75791539e-02    +1.00000265e+00    \n",
      "     33     33s_1/2      +1.20091868e+04    -1.65298645e-02    +1.00000138e+00    \n",
      "     34     34s_1/2      +2.06924979e+04    -1.55717933e-02    +1.00000075e+00    \n",
      "     35     35s_1/2      +3.40310129e+04    -1.46946655e-02    +1.00000043e+00    \n",
      "     36     36s_1/2      +5.38302322e+04    -1.38896135e-02    +1.00000026e+00    \n",
      "     37     37s_1/2      +8.26243463e+04    -1.31489512e-02    +1.00000016e+00    \n",
      "     38     38s_1/2      +1.24029065e+05    -1.24659907e-02    +1.00000010e+00    \n",
      "      :       :    \n",
      "     55     55s_1/2      +6.32629348e+07    -5.95061791e-03    +1.00000000e+00    \n",
      "     56     56s_1/2      +9.50494830e+07    -5.73998987e-03    +1.00000000e+00    \n",
      "     57     57s_1/2      +1.48001126e+08    -5.54035017e-03    +1.00000000e+00    \n",
      "     58     58s_1/2      +2.43113574e+08    -5.35094755e-03    +1.00000000e+00    \n",
      "     59     59s_1/2      +4.54506293e+08    -5.17109386e-03    +1.00000000e+00    \n",
      "     60     60s_1/2      +1.07334753e+09    -5.00015784e-03    +1.00000000e+00    \n",
      "  -----------------------------------------------------------------------------\n",
      ">> Use hydrogenic orbital from this symmetry block also for 1s_1/2.\n",
      ">> Use hydrogenic orbital from this symmetry block also for 2s_1/2.\n",
      "Generate hydrogenic orbital for subshell 2p_1/2 \n",
      "  -----------------------------------------------------------------------------\n",
      "   Index    Subshell     Energies [a.u.]    Dirac-E  [a.u.]     Delta-E / |E|    \n",
      "  -----------------------------------------------------------------------------\n",
      "      1      2p_1/2      -4.50269888e+00    -4.50269857e+00    -6.91708938e-08    \n",
      "      2      3p_1/2      -2.00096714e+00    -2.00095940e+00    -3.86904886e-06    \n",
      "      3      4p_1/2      -1.12549933e+00    -1.12543844e+00    -5.41025770e-05    \n",
      "      4      5p_1/2      -7.20400486e-01    -7.20234830e-01    -2.29951072e-04    \n",
      "      5      6p_1/2      -5.00841865e-01    -5.00139888e-01    -1.40159360e-03    \n",
      "      6      7p_1/2      -3.91485428e-01    -3.67436826e-01    -6.14291102e-02    \n",
      "      7      8p_1/2      -3.49598890e-01    -2.81311119e-01    -1.95331771e-01    \n",
      "      8      9p_1/2      -2.61761475e-01    -2.22265641e-01    -1.50884823e-01    \n",
      "      9     10p_1/2      -1.91479999e-01    -1.80031939e-01    -5.97872353e-02    \n",
      "     10     11p_1/2      -1.38459760e-01    -1.48784504e-01    +7.45685495e-02    \n",
      "     11     12p_1/2      -9.93780913e-02    -1.25018733e-01    +2.58011008e-01    \n",
      "     12     13p_1/2      -7.10836783e-02    -1.06523685e-01    +4.98567429e-01    \n",
      "     13     14p_1/2      -5.09019393e-02    -9.18486436e-02    +8.04423265e-01    \n",
      "     14     15p_1/2      -3.67303975e-02    -8.00097188e-02    +1.17829711e+00    \n",
      "     15     16p_1/2      -2.69592306e-02    -7.03205343e-02    +1.60840286e+00    \n",
      "     16     17p_1/2      -2.03451836e-02    -6.22904546e-02    +2.06168063e+00    \n",
      "     17     18p_1/2      -1.59307447e-02    -5.55612291e-02    +2.48767306e+00    \n",
      "     18     19p_1/2      -1.30071225e-02    -4.98663309e-02    +2.83377113e+00    \n",
      "     19     20p_1/2      -1.10591358e-02    -4.50041540e-02    +3.06940965e+00    \n",
      "     20     21p_1/2      -9.87350741e-03    -4.08199215e-02    +3.13428783e+00    \n",
      "     21     22p_1/2      +3.90206612e-01    -3.71932146e-02    +1.09531672e+00    \n",
      "     22     23p_1/2      +3.31827545e+00    -3.40292102e-02    +1.01025509e+00    \n",
      "     23     24p_1/2      +1.08803000e+01    -3.12524195e-02    +1.00287239e+00    \n",
      "     24     25p_1/2      +2.82338195e+01    -2.88021434e-02    +1.00102013e+00    \n",
      "     25     26p_1/2      +6.59854243e+01    -2.66291266e-02    +1.00040356e+00    \n",
      "     26     27p_1/2      +1.45687867e+02    -2.46930634e-02    +1.00016949e+00    \n",
      "     27     28p_1/2      +3.10533661e+02    -2.29607143e-02    +1.00007394e+00    \n",
      "     28     29p_1/2      +6.45434359e+02    -2.14044706e-02    +1.00003316e+00    \n",
      "     29     30p_1/2      +1.31210512e+03    -2.00012468e-02    +1.00001524e+00    \n",
      "     30     31p_1/2      +2.60237480e+03    -1.87316200e-02    +1.00000720e+00    \n",
      "     31     32p_1/2      +4.99986575e+03    -1.75791539e-02    +1.00000352e+00    \n",
      "     32     33p_1/2      +9.21993649e+03    -1.65298645e-02    +1.00000179e+00    \n",
      "     33     34p_1/2      +1.62105279e+04    -1.55717933e-02    +1.00000096e+00    \n",
      "     34     35p_1/2      +2.71669430e+04    -1.46946655e-02    +1.00000054e+00    \n",
      "     35     36p_1/2      +4.36433180e+04    -1.38896135e-02    +1.00000032e+00    \n",
      "     36     37p_1/2      +6.77761587e+04    -1.31489512e-02    +1.00000019e+00    \n",
      "     37     38p_1/2      +1.02595088e+05    -1.24659907e-02    +1.00000012e+00    \n",
      "     38     39p_1/2      +1.52429653e+05    -1.18348904e-02    +1.00000008e+00    \n",
      "      :       :    \n",
      "     55     56p_1/2      +7.44463010e+07    -5.73998987e-03    +1.00000000e+00    \n",
      "     56     57p_1/2      +1.11991148e+08    -5.54035017e-03    +1.00000000e+00    \n",
      "     57     58p_1/2      +1.74326907e+08    -5.35094755e-03    +1.00000000e+00    \n",
      "     58     59p_1/2      +2.87668439e+08    -5.17109386e-03    +1.00000000e+00    \n",
      "     59     60p_1/2      +5.37617989e+08    -5.00015784e-03    +1.00000000e+00    \n",
      "     60     61p_1/2      +1.26790739e+09    -4.83755953e-03    +1.00000000e+00    \n",
      "  -----------------------------------------------------------------------------\n",
      ">> Use hydrogenic orbital from this symmetry block also for 2p_1/2.\n",
      "Generate hydrogenic orbital for subshell 2p_3/2 \n",
      "  -----------------------------------------------------------------------------\n",
      "   Index    Subshell     Energies [a.u.]    Dirac-E  [a.u.]     Delta-E / |E|    \n",
      "  -----------------------------------------------------------------------------\n",
      "      1      2p_3/2      -4.50053984e+00    -4.50053930e+00    -1.19800113e-07    \n",
      "      2      3p_3/2      -2.00033044e+00    -2.00031959e+00    -5.42524482e-06    \n",
      "      3      4p_3/2      -1.12535230e+00    -1.12516853e+00    -1.63298442e-04    \n",
      "      4      5p_3/2      -7.21812733e-01    -7.20096642e-01    -2.37747308e-03    \n",
      "      5      6p_3/2      -5.88640690e-01    -5.00059921e-01    -1.50483597e-01    \n",
      "      6      7p_3/2      -4.94091695e-01    -3.67386470e-01    -2.56440710e-01    \n",
      "      7      8p_3/2      -3.59119769e-01    -2.81277385e-01    -2.16758837e-01    \n",
      "      8      9p_3/2      -2.63766279e-01    -2.22241949e-01    -1.57428502e-01    \n",
      "      9     10p_3/2      -1.92354920e-01    -1.80014668e-01    -6.41535525e-02    \n",
      "     10     11p_3/2      -1.38872484e-01    -1.48771528e-01    +7.12815335e-02    \n",
      "     11     12p_3/2      -9.95719755e-02    -1.25008738e-01    +2.55461062e-01    \n",
      "     12     13p_3/2      -7.11733657e-02    -1.06515824e-01    +4.96568597e-01    \n",
      "     13     14p_3/2      -5.09427169e-02    -9.18423497e-02    +8.02855350e-01    \n",
      "     14     15p_3/2      -3.67485977e-02    -8.00046017e-02    +1.17707904e+00    \n",
      "     15     16p_3/2      -2.69671150e-02    -7.03163180e-02    +1.60748389e+00    \n",
      "     16     17p_3/2      -2.03482925e-02    -6.22869395e-02    +2.06104011e+00    \n",
      "     17     18p_3/2      -1.59315437e-02    -5.55582679e-02    +2.48731227e+00    \n",
      "     18     19p_3/2      -1.30067433e-02    -4.98638131e-02    +2.83368932e+00    \n",
      "     19     20p_3/2      -1.10582249e-02    -4.50019952e-02    +3.06954964e+00    \n",
      "     20     21p_3/2      -9.87114174e-03    -4.08180568e-02    +3.13508972e+00    \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     21     22p_3/2      +1.24996175e-01    -3.71915928e-02    +1.29754185e+00    \n",
      "     22     23p_3/2      +2.69273197e+00    -3.40277908e-02    +1.01263690e+00    \n",
      "     23     24p_3/2      +9.55804164e+00    -3.12511703e-02    +1.00326962e+00    \n",
      "     24     25p_3/2      +2.55081379e+01    -2.88010381e-02    +1.00112909e+00    \n",
      "     25     26p_3/2      +6.04337789e+01    -2.66281441e-02    +1.00044062e+00    \n",
      "     26     27p_3/2      +1.34480686e+02    -2.46921860e-02    +1.00018361e+00    \n",
      "     27     28p_3/2      +2.88122942e+02    -2.29599276e-02    +1.00007969e+00    \n",
      "     28     29p_3/2      +6.01200790e+02    -2.14037625e-02    +1.00003560e+00    \n",
      "     29     30p_3/2      +1.22657485e+03    -2.00006072e-02    +1.00001631e+00    \n",
      "     30     31p_3/2      +2.44237088e+03    -1.87310403e-02    +1.00000767e+00    \n",
      "     31     32p_3/2      +4.71496718e+03    -1.75786269e-02    +1.00000373e+00    \n",
      "     32     33p_3/2      +8.74391356e+03    -1.65293840e-02    +1.00000189e+00    \n",
      "     33     34p_3/2      +1.54668470e+04    -1.55713539e-02    +1.00000101e+00    \n",
      "     34     35p_3/2      +2.60705119e+04    -1.46942628e-02    +1.00000056e+00    \n",
      "     35     36p_3/2      +4.20943262e+04    -1.38892433e-02    +1.00000033e+00    \n",
      "     36     37p_3/2      +6.56501441e+04    -1.31486103e-02    +1.00000020e+00    \n",
      "     37     38p_3/2      +9.97331750e+04    -1.24656760e-02    +1.00000012e+00    \n",
      "     38     39p_3/2      +1.48629201e+05    -1.18345992e-02    +1.00000008e+00    \n",
      "      :       :    \n",
      "     55     56p_3/2      +7.64518397e+07    -5.73989154e-03    +1.00000000e+00    \n",
      "     56     57p_3/2      +1.16562068e+08    -5.54025693e-03    +1.00000000e+00    \n",
      "     57     58p_3/2      +1.85376103e+08    -5.35085904e-03    +1.00000000e+00    \n",
      "     58     59p_3/2      +3.19183797e+08    -5.17100978e-03    +1.00000000e+00    \n",
      "     59     60p_3/2      +6.47519682e+08    -5.00007789e-03    +1.00000000e+00    \n",
      "     60     61p_3/2      +2.01326894e+09    -4.83748345e-03    +1.00000000e+00    \n",
      "  -----------------------------------------------------------------------------\n",
      ">> Use hydrogenic orbital from this symmetry block also for 2p_3/2.\n",
      "(Re-) Define the standard grid with 392 grid points.\n",
      ">>  for various B-spline matrices:\n",
      "\n",
      "Iteration 1 for symmetries ... \n",
      "  1s_1/2::  en [a.u.] = -7.5200261e+00;   self-cons'cy = 4.1086e-01  [1.0000e+02 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -1.1797638e-02;   self-cons'cy = 9.9477e-01  [1.0000e+02 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -2.6580075e-01;   self-cons'cy = 8.8852e-01  [1.0000e+02 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -3.7294812e-01;   self-cons'cy = 8.4695e-01  [1.0000e+02 for sym-block kappa = -2]\n",
      "\n",
      "Iteration 2 for symmetries ... \n",
      "  1s_1/2::  en [a.u.] = -1.0847250e+01;   self-cons'cy = 1.8115e-01  [1.0149e+00 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -9.2554877e-01;   self-cons'cy = 9.7483e-01  [1.0149e+00 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -8.0889267e-02;   self-cons'cy = 5.3336e-01  [1.0028e+00 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -9.3222493e-02;   self-cons'cy = 6.0005e-01  [1.0022e+00 for sym-block kappa = -2]\n",
      "\n",
      "Iteration 3 for symmetries ... \n",
      "  1s_1/2::  en [a.u.] = -9.8287731e+00;   self-cons'cy = 4.9259e-02  [1.0036e+00 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -4.6270685e-01;   self-cons'cy = 3.3340e-01  [1.0036e+00 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -2.2971831e-01;   self-cons'cy = 4.7915e-01  [4.7915e-01 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -1.9405288e-01;   self-cons'cy = 3.5099e-01  [3.5099e-01 for sym-block kappa = -2]\n",
      "\n",
      "Iteration 4 for symmetries ... \n",
      "  1s_1/2::  en [a.u.] = -9.8305270e+00;   self-cons'cy = 8.9215e-05  [9.9572e-01 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -4.3112119e-01;   self-cons'cy = 3.5338e-02  [9.9572e-01 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -1.3712735e-01;   self-cons'cy = 2.5240e-01  [2.5240e-01 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -1.5264566e-01;   self-cons'cy = 1.1943e-01  [1.1943e-01 for sym-block kappa = -2]\n",
      "\n",
      "Iteration 5 for symmetries ... \n",
      "  1s_1/2::  en [a.u.] = -9.9219196e+00;   self-cons'cy = 4.6269e-03  [9.8274e-01 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -4.7105767e-01;   self-cons'cy = 4.4267e-02  [9.8274e-01 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -1.6105986e-01;   self-cons'cy = 8.0260e-02  [8.0260e-02 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -1.5588522e-01;   self-cons'cy = 1.0500e-02  [1.0500e-02 for sym-block kappa = -2]\n",
      "\n",
      "Iteration 6 for symmetries ... \n",
      "  1s_1/2::  en [a.u.] = -9.8785043e+00;   self-cons'cy = 2.1926e-03  [2.5578e-01 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -4.5450705e-01;   self-cons'cy = 1.7882e-02  [2.5578e-01 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -1.5900870e-01;   self-cons'cy = 6.4085e-03  [6.4085e-03 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -1.5898520e-01;   self-cons'cy = 9.8452e-03  [9.8452e-03 for sym-block kappa = -2]\n",
      "\n",
      "Iteration 7 for symmetries ... \n",
      "  1s_1/2::  en [a.u.] = -9.8876275e+00;   self-cons'cy = 4.6156e-04  [7.4534e-03 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -4.5725498e-01;   self-cons'cy = 3.0139e-03  [7.4534e-03 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -1.5681807e-01;   self-cons'cy = 6.9362e-03  [6.9362e-03 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -1.5695988e-01;   self-cons'cy = 6.4104e-03  [6.4104e-03 for sym-block kappa = -2]\n",
      "\n",
      "Iteration 8 for symmetries ... \n",
      "  1s_1/2::  en [a.u.] = -9.8887820e+00;   self-cons'cy = 5.8378e-05  [4.0564e-02 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -4.5802540e-01;   self-cons'cy = 8.4173e-04  [4.0564e-02 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -1.5814161e-01;   self-cons'cy = 4.2023e-03  [4.2023e-03 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -1.5756386e-01;   self-cons'cy = 1.9203e-03  [1.9203e-03 for sym-block kappa = -2]\n",
      "\n",
      "Iteration 9 for symmetries ... \n",
      "  1s_1/2::  en [a.u.] = -9.8870070e+00;   self-cons'cy = 8.9755e-05  [1.9165e-02 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -4.5727888e-01;   self-cons'cy = 8.1560e-04  [1.9165e-02 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -1.5777061e-01;   self-cons'cy = 1.1744e-03  [1.1744e-03 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -1.5752403e-01;   self-cons'cy = 1.2639e-04  [1.2639e-04 for sym-block kappa = -2]\n",
      "\n",
      "Iteration 10 for symmetries ... \n",
      "  1s_1/2::  en [a.u.] = -9.8877423e+00;   self-cons'cy = 3.7184e-05  [4.3381e-03 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -4.5755039e-01;   self-cons'cy = 2.9679e-04  [4.3381e-03 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -1.5778226e-01;   self-cons'cy = 3.6927e-05  [4.3854e-05 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -1.5746560e-01;   self-cons'cy = 1.8550e-04  [1.8550e-04 for sym-block kappa = -2]\n",
      "\n",
      "Iteration 11 for symmetries ... \n",
      "  1s_1/2::  en [a.u.] = -9.8876063e+00;   self-cons'cy = 6.8786e-06  [3.4489e-04 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -4.5751199e-01;   self-cons'cy = 4.1965e-05  [3.4489e-04 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -1.5782460e-01;   self-cons'cy = 1.3416e-04  [1.3416e-04 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -1.5749971e-01;   self-cons'cy = 1.0827e-04  [1.0827e-04 for sym-block kappa = -2]\n",
      "\n",
      "Iteration 12 for symmetries ... \n",
      "  1s_1/2::  en [a.u.] = -9.8875788e+00;   self-cons'cy = 1.3910e-06  [7.3183e-04 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -4.5749594e-01;   self-cons'cy = 1.7540e-05  [7.3183e-04 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -1.5780205e-01;   self-cons'cy = 7.1466e-05  [7.1466e-05 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -1.5749035e-01;   self-cons'cy = 2.9699e-05  [2.9699e-05 for sym-block kappa = -2]\n",
      "\n",
      "Iteration 13 for symmetries ... \n",
      "  1s_1/2::  en [a.u.] = -9.8876100e+00;   self-cons'cy = 1.5789e-06  [3.1856e-04 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -4.5750887e-01;   self-cons'cy = 1.4125e-05  [3.1856e-04 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -1.5780777e-01;   self-cons'cy = 1.8123e-05  [1.8123e-05 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -1.5749058e-01;   self-cons'cy = 7.2018e-07  [1.7775e-06 for sym-block kappa = -2]\n",
      "\n",
      "Iteration 14 for symmetries ... \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  1s_1/2::  en [a.u.] = -9.8875981e+00;   self-cons'cy = 6.0317e-07  [6.3336e-05 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -4.5750453e-01;   self-cons'cy = 4.7378e-06  [6.3336e-05 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -1.5780785e-01;   self-cons'cy = 2.5151e-07  [5.1385e-07 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -1.5749168e-01;   self-cons'cy = 3.4875e-06  [3.4875e-06 for sym-block kappa = -2]\n",
      "\n",
      "Iteration 15 for symmetries ... \n",
      "  1s_1/2::  en [a.u.] = -9.8875999e+00;   self-cons'cy = 9.4910e-08  [9.7914e-06 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -4.5750501e-01;   self-cons'cy = 5.2256e-07  [9.7914e-06 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -1.5780707e-01;   self-cons'cy = 2.4696e-06  [2.4696e-06 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -1.5749110e-01;   self-cons'cy = 1.8243e-06  [1.8243e-06 for sym-block kappa = -2]\n",
      "\n",
      "Iteration 16 for symmetries ... \n",
      "  1s_1/2::  en [a.u.] = -9.8876005e+00;   self-cons'cy = 3.0982e-08  [1.3049e-05 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -4.5750534e-01;   self-cons'cy = 3.5398e-07  [1.3049e-05 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -1.5780744e-01;   self-cons'cy = 1.1929e-06  [1.1929e-06 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -1.5749124e-01;   self-cons'cy = 4.5393e-07  [4.5393e-07 for sym-block kappa = -2]\n",
      "\n",
      "Iteration 17 for symmetries ... \n",
      "  1s_1/2::  en [a.u.] = -9.8876000e+00;   self-cons'cy = 2.7663e-08  [5.1990e-06 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -4.5750511e-01;   self-cons'cy = 2.4333e-07  [5.1990e-06 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -1.5780736e-01;   self-cons'cy = 2.7230e-07  [2.7230e-07 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -1.5749125e-01;   self-cons'cy = 1.0951e-08  [3.3630e-07 for sym-block kappa = -2]\n",
      "\n",
      "Iteration 18 for symmetries ... \n",
      "  1s_1/2::  en [a.u.] = -9.8876002e+00;   self-cons'cy = 9.7199e-09  [8.7749e-07 for sym-block kappa = -1]\n",
      "  2s_1/2::  en [a.u.] = -4.5750518e-01;   self-cons'cy = 7.4948e-08  [8.7749e-07 for sym-block kappa = -1]\n",
      "  2p_1/2::  en [a.u.] = -1.5780735e-01;   self-cons'cy = 1.9275e-08  [2.4890e-08 for sym-block kappa = 1]\n",
      "  2p_3/2::  en [a.u.] = -1.5749123e-01;   self-cons'cy = 6.4179e-08  [4.1326e-07 for sym-block kappa = -2]\n",
      ">>> Sign changed for orbital 2p_3/2\n",
      "Compute CI matrix of dimension 2 x 2 for the symmetry 0^+ ...   ... done.\n",
      "Compute CI matrix of dimension 1 x 1 for the symmetry 1^+ ...   ... done.\n",
      "Compute CI matrix of dimension 2 x 2 for the symmetry 2^+ ...   ... done.\n",
      "\n",
      "  Eigenenergies:\n",
      "\n",
      "  Level  J Parity          Hartrees                    eV                   [eV]\n",
      "\n",
      "     1    0 +     -3.769232142213917e+01    -1.025660308207744e+03    -1.025660308207744e+03 \n",
      "     2    1 +     -3.769222590220146e+01    -1.025657708977847e+03    -1.025657708977847e+03 \n",
      "     3    2 +     -3.769203577987614e+01    -1.025652535485861e+03    -1.025652535485861e+03 \n",
      "     4    2 +     -3.763845411653460e+01    -1.024194504161081e+03    -1.024194504161081e+03 \n",
      "     5    0 +     -3.755793951144176e+01    -1.022003590161652e+03    -1.022003590161652e+03 \n",
      "\n",
      "  Energy of each level relative to immediately lower level:\n",
      "\n",
      "  Level  J Parity          Hartrees                    eV                   [eV]\n",
      "\n",
      "     2    1 +     9.551993770173794e-05    2.599229897408343e-03    2.599229897408343e-03 \n",
      "     3    2 +     1.901223253213402e-04    5.173491985339009e-03    5.173491985339009e-03 \n",
      "     4    2 +     5.358166334153935e-02    1.458031324780311e+00    1.458031324780311e+00 \n",
      "     5    0 +     8.051460509284425e-02    2.190913999429243e+00    2.190913999429243e+00 \n",
      "\n",
      "  Energy of each level relative to lowest level:\n",
      "\n",
      "  Level  J Parity          Hartrees                    eV                   [eV]\n",
      "\n",
      "     2    1 +     9.551993770173794e-05    2.599229897408343e-03    2.599229897408343e-03 \n",
      "     3    2 +     2.856422630230782e-04    7.772721882747352e-03    7.772721882747352e-03 \n",
      "     4    2 +     5.386730560456243e-02    1.465804046663058e+00    1.465804046663058e+00 \n",
      "     5    0 +     1.343819106974067e-01    3.656718046092300e+00    3.656718046092300e+00 \n"
     ]
    }
   ],
   "source": [
    "basis     = Basics.performSCF([Configuration(\"1s^2 2s^2 2p^2\")], nucModel, grid, AsfSettings())\n",
    "multiplet = Basics.performCI(basis, nucModel, grid, AsfSettings());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "From this multiplet, we just take the *lowest* level (the ground level) in order to calculate the density for the (electronic) potential. In JAC, we have different choices for generating such an **electronic** potential that just arise from the electronic density, and as seen from ? compute at the REPL:\n",
    "\n",
    "   •    (\"radial potential: core-Hartree\", grid::Radial.Grid, level::Level) \n",
    "       ... to compute a (radial) core-Hartree potential for the given level; \n",
    "           cf. Basics.computePotentialCoreHartree. A potential::RadialPotential \n",
    "           is returned. \n",
    "\n",
    "   •    (\"radial potential: Kohn-Sham\", grid::Radial.Grid, level::Level) \n",
    "       ... to compute a (radial) Kohn-Sham potential for the given level; \n",
    "           cf. Basics.computePotentialKohnSham. A potential::RadialPotential is \n",
    "           returned.\n",
    "\n",
    "   •    (\"radial potential: Dirac-Fock-Slater\", grid::Radial.Grid, level::Level) \n",
    "       ... to compute a (radial) Dirac-Fock-Slater potential for the given level; this \n",
    "           potential is rather simple but includes some undesired self-interaction and \n",
    "           exhibits a wrong asymptotic behaviour. Cf. Basics.computePotentialDFS. \n",
    "           A potential::RadialPotential is returned.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "coreSubshells = Subshell[1s_1/2, 2s_1/2];    subshells = Subshell[1s_1/2, 2s_1/2, 2p_1/2, 2p_3/2]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "DFS (radial) potential ... defined on 392 grid points ...\n",
       "Zr:    [0.0, 2.6955607059564223e-6, -5.434596604675834e-7, -5.82282910371863e-7, -1.6059227517515662e-6, -4.559757581708367e-6, -6.406366344836484e-6, -6.917100009381071e-6, -6.793141197899899e-6, -7.566326471167558e-6  …  -1.6521180096761336e-5, -1.6997984793507585e-5, -1.840888391790231e-5, -2.1396074388405442e-5, -2.5194990715107997e-5, -2.831799102219741e-5, -3.063234020923076e-5, -3.21115206158104e-5, -3.3027169688644986e-5, -3.533447353281762e-5]  ...  [-6.0, -6.000000000000003, -5.999999999999999, -6.000000000000004, -5.999999999999999, -6.000000000000002, -5.999999999999999, -6.0, -6.000000000000001, -5.999999999999996  …  -6.000000000000004, -6.000000000000001, -6.000000000000003, -6.000000000000005, -5.999999999999999, -5.999999999999999, -6.000000000000002, -6.0000000000000036, -6.0, -5.999999999999998]\n",
       "Radial grid:  rnt = 2.0e-6,  h = 0.05,  hp = 0.0,  NoPoints = 392,  ntL = 69,  ntS = 71, orderL = 7,  orderS = 8,  nsL = 62,  nsS = 63,  mesh = JAC.Radial.MeshGL(), ...  \n",
       "r:    [1.780504510281598e-8, 9.042759120177923e-8, 2.0787030685460475e-7]  ...  [564.4302351890452, 595.097239730825, 614.0606632064549]\n",
       "wr:   [4.530145586285996e-8, 9.785739473901763e-8, 1.3358660616640612e-7]  ...  [34.882547080539766, 25.552825070730517, 11.829256033233948]\n",
       "tS:   [0.0, 0.0, 0.0]  ...  [618.7099715607404, 618.7099715607404, 618.7099715607404]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "level  = multiplet.levels[1]\n",
    "potCH  = compute(\"radial potential: core-Hartree\", grid, level)\n",
    "potDFS = compute(\"radial potential: Dirac-Fock-Slater\", grid, level)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "All these potentials are represented as *effective charge* $Z(r)$ and can be added to the nuclear potential, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "nuclear-potential: Fermi-distributed+DFS (radial) potential ... defined on 392 grid points ...\n",
       "Zr:    [0.0, 0.015622513768086672, 0.03590526099223162, 0.06043026291944624, 0.08495293134027639, 0.10523281020065203, 0.11777238754107956, 0.12581795181779093, 0.14609385219119195, 0.1788775273906122  …  0.3111128557686238, 0.32311941468563676, 0.35183978014834993, 0.39826148078365287, 0.45434229564729456, 0.5103706313406959, 0.5566687156429055, 0.5852757770687644, 0.6022338661276299, 0.6427823050363154]  ...  [0.0, -8.881784197001252e-16, 1.7763568394002505e-15, -2.6645352591003757e-15, 1.7763568394002505e-15, -1.7763568394002505e-15, 1.7763568394002505e-15, 0.0, 8.881784197001252e-16, 6.217248937900877e-15  …  -2.6645352591003757e-15, 0.0, -8.881784197001252e-16, -3.552713678800501e-15, 1.7763568394002505e-15, 1.7763568394002505e-15, -1.7763568394002505e-15, -2.6645352591003757e-15, 8.881784197001252e-16, 2.6645352591003757e-15]\n",
       "Radial grid:  rnt = 2.0e-6,  h = 0.05,  hp = 0.0,  NoPoints = 392,  ntL = 69,  ntS = 71, orderL = 7,  orderS = 8,  nsL = 62,  nsS = 63,  mesh = JAC.Radial.MeshGL(), ...  \n",
       "r:    [1.780504510281598e-8, 9.042759120177923e-8, 2.0787030685460475e-7]  ...  [564.4302351890452, 595.097239730825, 614.0606632064549]\n",
       "wr:   [4.530145586285996e-8, 9.785739473901763e-8, 1.3358660616640612e-7]  ...  [34.882547080539766, 25.552825070730517, 11.829256033233948]\n",
       "tS:   [0.0, 0.0, 0.0]  ...  [618.7099715607404, 618.7099715607404, 618.7099715607404]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "add(nucPotential, potDFS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "As seen from this output, the effective charge of the overall atomic potential goes to *zero* as required for a neutral atom. We can display these potentials also graphically by using `JAC.plot()`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **The following part of this nootebook is still under construction.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.6.2",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}