Helium 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  = 2
Zb = 0

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

Nmo = [[1]]
N   = [[2]]

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       -3.62587      -3.79963e-01       +1.00000e+00
    2       -3.05643      -5.22432e-01       +1.86311e-01
    3       -2.89905      -5.58398e-01       +5.42842e-02
    4       -2.85419      -5.67581e-01       +1.57196e-02
    5       -2.84080      -5.69832e-01       +4.71360e-03
    6       -2.83662      -5.70305e-01       +1.47100e-03
    7       -2.83525      -5.70355e-01       +4.83381e-04
    8       -2.83477      -5.70328e-01       +1.75088e-04
    9       -2.83459      -5.70298e-01       +8.00280e-05
   10       -2.83452      -5.70279e-01       +3.77715e-05
   11       -2.83448      -5.70268e-01       +1.81939e-05
   12       -2.83447      -5.70262e-01       +8.87283e-06
 Total Energy: -2.834468599168821