### Fluorine LDA

import numpy as np
from CADMium import Psgrid
from CADMium import Kohnsham

#Distance of the nucley from grid center
a =  1.0

#Nuclear charges on centers AB
Za  = 9
Zb = 0

#Set polaization. 1 Unpolarized, 2 Polarized
pol = 1

Nmo = [[3],[2]]
N   = [[6],[3]]  # Why are we considering molecular orbitals (?)

optKS = {
        "interaction_type" : "dft",
        "SYM" : False,
        "FRACTIONAL" : True,
        }

#Grid Options
NP = 7 #Number of points per block
NM =  [14,14] #Number of blocks [angular, radial]
L = np.arccosh(15./a)#Maximum radial coordinate value
loc = np.array(range(-4,5)) #Non inclusive on upper bound

#Create and initialize grid object
grid = Psgrid(NP, NM, a, L, loc)
grid.initialize()

#Kohn Sham object
KS = Kohnsham(grid, Za, Zb, pol, Nmo, N, optKS)
KS.scf()

print(f" Total Energy: {KS.E.E}")

 iter    Total Energy     HOMO Eigenvalue         Res

-----------------------------------------------------------

    1      -119.69484      -3.88326e-02       +1.00000e+00
    2      -100.08111      -1.05584e+00       +1.95978e-01
    3      -101.61887      -5.61417e-01       +3.48165e-02
    4      -100.30326      -4.66725e-01       +1.62265e-02
    5      -99.59278      -4.32251e-01       +7.40656e-03
    6      -99.37333      -4.31746e-01       +3.45818e-03
    7      -99.18043      -4.07035e-01       +1.94490e-03
    8      -99.12508      -4.04224e-01       +7.16056e-04
    9      -99.10622      -4.02651e-01       +3.38915e-04
   10      -99.09828      -4.01929e-01       +1.61795e-04
   11      -99.09512      -4.01587e-01       +7.80343e-05
   12      -99.09384      -4.01426e-01       +3.80132e-05
   13      -99.09334      -4.01350e-01       +1.87328e-05
   14      -99.09314      -4.01314e-01       +9.41549e-06
 Total Energy: -99.09314377654036

Compare againts Nist Total Energy: Etot = -99.099648