Oxygen LDAΒΆ
[1]:
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 = 8
Zb = 0
#Set polaization. 1 Unpolarized, 2 Polarized
pol = 1
Nmo = [[2],[2]]
N = [[4],[4]]
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 -97.38518 +4.09431e-02 +1.00000e+00
2 -76.76286 -1.19860e+00 +2.68650e-01
3 -76.25430 -6.74205e-01 +1.30179e-02
4 -75.29033 -4.76172e-01 +1.28034e-02
5 -74.76910 -3.92386e-01 +6.97111e-03
6 -74.63917 -3.53979e-01 +1.74087e-03
7 -74.46442 -3.36203e-01 +2.34675e-03
8 -74.43862 -3.27945e-01 +4.49173e-04
9 -74.42824 -3.24290e-01 +2.42419e-04
10 -74.42553 -3.22561e-01 +1.28093e-04
11 -74.42502 -3.21757e-01 +6.66725e-05
12 -74.42513 -3.21380e-01 +3.42829e-05
13 -74.42533 -3.21204e-01 +1.74472e-05
14 -74.42549 -3.21121e-01 +8.80796e-06
Total Energy: -74.42549267778712
Compare againts NIST Energy: -74.473077