Beryllium 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  = 4
Zb = 0

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

Nmo = [[2]]
N   = [[4]]

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

#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      -17.52412      -2.50572e-02       +1.00000e+00
    2      -15.40680      -1.79335e-01       +1.37428e-01
    3      -14.74245      -2.07728e-01       +4.50636e-02
    4      -14.53904      -2.11670e-01       +1.39912e-02
    5      -14.47292      -2.10250e-01       +4.56819e-03
    6      -14.45718      -2.10086e-01       +1.37112e-03
    7      -14.44591      -2.07539e-01       +8.23717e-04
    8      -14.44550      -2.06687e-01       +4.47611e-04
    9      -14.44574      -2.06240e-01       +2.38841e-04
   10      -14.44601      -2.06010e-01       +1.25535e-04
   11      -14.44620      -2.05892e-01       +6.52941e-05
   12      -14.44632      -2.05832e-01       +3.36878e-05
   13      -14.44639      -2.05801e-01       +1.72731e-05
   14      -14.44643      -2.05785e-01       +8.81403e-06
 Total Energy: -14.446430904094559