LDA Neon

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  = 10
Zb = 0

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

Nmo = [[3],[2]]
N   = [[6],[4]] # Why do we need molecular orbitals (?)


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

#Grid Options
NP = 7 #Number of points per block
NM =  [10,10] #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      -173.70297      +4.28212e-02       +1.00000e+00
    2      -129.49817      -2.47347e+00       +3.41355e-01
    3      -131.96403      -1.23654e+00       +1.98279e-02
    4      -130.22035      -8.18264e-01       +1.33902e-02
    5      -129.02979      -6.47305e-01       +9.22699e-03
    6      -128.88827      -5.67080e-01       +1.19854e-03
    7      -128.33113      -5.28687e-01       +4.34144e-03
    8      -128.27004      -5.11861e-01       +4.76242e-04
    9      -128.23873      -5.04640e-01       +2.44115e-04
   10      -128.23216      -5.01002e-01       +8.11805e-05
   11      -128.22964      -4.99340e-01       +4.27014e-05
   12      -128.22927      -4.98549e-01       +2.21693e-05
   13      -128.22934      -4.98177e-01       +1.14056e-05
   14      -128.22949      -4.98001e-01       +5.81501e-06
 Total Energy: -128.2294916671537

Compare Againts NIST Total Energy: Etot = -128.233481