Solution of the atom within LDA¶
In this excercise we will solve the multielectron atom in LDA approximation.
We will test it on He and oxygen by computing the total energy and charge density.
We will plot charge density and compute the total energy, which will be compared to the reference data at NIST database: https://www.nist.gov/pml/atomic-reference-data-electronic-structure-calculations/atomic-reference-data-electronic-7
We want to solve the Schroedinger equation for an atom with nucleous charge Z. We will approximate electron-electron interaction with an effective potential, which is computed by so-called ”local density approximation” (LDA). In this theory, the classical (Hartree) potential is treated exactly, while the rest of the interaction is ”approximated” by so called local exchange-correlation functional. We will not go into details of this functional, we will just use it here.
The Schroedinger equation we are solving is \begin{eqnarray} (-\frac{\hbar^2}{2m}\nabla^2-\frac{Z e^2}{4\pi\varepsilon_0 r} + V_H(r) + V_{xc}(r))\psi(\vec{r})=\varepsilon \psi(\vec{r}) \end{eqnarray}
The first two terms are appearing in Hydrogen problem, and we already coded them. The Hartree is the electrostatic potential, and the exchange-correlation potential has an approximate form, which depends only the charge density of the system. We will use the module excor.py, where the function $V_{xc}(\rho)$ is tabulated. We will use it as $V_{xc}(r)== V_{xc}(\rho(r))$, corresponding to the local density approximation.
Using atomic units, the equation is: \begin{eqnarray} \left(-\nabla^2-\frac{2Z}{r}+V_H(r)+V_{xc}(r)-\varepsilon\right)\psi(r)=0 \end{eqnarray} and when spherical symmetry of an atom is taken into account: \begin{eqnarray} u''(r)- \left(\frac{l(l+1)}{r^2}-\frac{2Z}{r}+V_H(r)+V_{xc}(r)-\varepsilon\right)u(r)=0 \end{eqnarray}
First we take the code from the Hydrogen project and repeat.
from scipy import integrate
from scipy import optimize
from numpy import *
from numpy.polynomial import Polynomial
from pylab import *
from numba import jit
@jit(nopython=True)
def Numerov(f, x0, dx, dh):
""" Given precumputed function f(x) solved the differential equation
x''(t) = f(t) x(t)
input: x0 = x(t=0), and dx = dx/dt(t=0)
"""
x = zeros(len(f))
x[0] = x0
x[1] = x0 + dh*dx
h2 = dh**2
h12 = h2/12.
w0 = x0*(1-h12*f[0])
w1 = x[1]*(1-h12*f[1])
xi = x[1]
fi = f[1]
for i in range(2,len(f)):
w2 = 2*w1-w0 + h2*fi*xi
fi = f[i]
xi = w2/(1-h12*fi)
x[i]=xi
(w0,w1) = (w1,w2)
return x
@jit(nopython=True)
def fShrod(En,l,R):
return l*(l+1.)/R**2 - 2./R - En
def ComputeSchrod(En, R, l):
f = fShrod(En,l,R[::-1])
ur = Numerov(f, 0.0, -1e-7, -R[1]+R[0])[::-1]
norm = integrate.simps(ur**2, x=R)
return ur/sqrt(abs(norm))
def Shoot(En, R, l):
ur = ComputeSchrod(En, R, l)
ur = ur/R**l # expecting urn \propto R
#f0,f1 = ur[0],ur[1]
#f_at_0 = f0 + (f1-f0)*(0-R[0])/(R[1]-R[0]) # extrapolation to zero
#return f_at_0
#poly = polyfit(R[:4], ur[:4], deg=3)
#return polyval(poly, 0.0)
poly = Polynomial.fit(R[:4], ur[:4], deg=3)
return poly(0.0)
def FindBoundStates(R, l, nmax, Esearch):
n=0
Ebnd=[]
u0 = Shoot(Esearch[0],R,l)
for i in range(1,len(Esearch)):
u1 = Shoot(Esearch[i],R,l)
if u0*u1 < 0:
#print 'Sign change at', Esearch[i-1], Esearch[i]
Ebound = optimize.brentq(Shoot,Esearch[i-1],Esearch[i],xtol=1e-15,args=(R,l))
Ebnd.append( (l,Ebound) )
if len(Ebnd)>nmax: break
n += 1
print('Found bound state at E=%14.9f' % Ebound)
u0 = u1
return Ebnd
def cmpKey(x):
return x[1]*1000 + x[0] # energy has large wait, but degenerate energy states are sorted by l
def ChargeDensity(Bnd,R,Z):
rho = zeros(len(R))
N=0.
for (l,En) in Bnd:
ur = ComputeSchrod(En, R, l)
dN = 2*(2*l+1)
if N+dN <= Z:
ferm = 1.
else:
ferm = (Z-N)/float(dN)
drho = ur**2 * ferm * dN/(4*pi*R**2)
rho += drho
N += dN
print('adding state', (l,En), 'with fermi=', ferm)
if N>=Z: break
return rho
%matplotlib inline
Esearch = -1.2/arange(1,20,0.2)**2
R = linspace(1e-6,100,2000)
Z=28
nmax = 5
Bnd=[]
for l in range(nmax-1):
Bnd += FindBoundStates(R,l,nmax-l,Esearch)
Bnd = sorted(Bnd, key=cmpKey)
z = 28. # like Ni atom
rho = ChargeDensity(Bnd,R,Z)
plot(R, rho*(4*pi*R**2), label='charge density')
xlim([0,25])
legend(loc='best')
show()
Found bound state at E= -0.999999943 Found bound state at E= -0.249999990 Found bound state at E= -0.111111108 Found bound state at E= -0.062499999 Found bound state at E= -0.039999942 Found bound state at E= -0.249999998 Found bound state at E= -0.111111111 Found bound state at E= -0.062500000 Found bound state at E= -0.039999957 Found bound state at E= -0.111111111 Found bound state at E= -0.062500000 Found bound state at E= -0.039999977 Found bound state at E= -0.062500000 Found bound state at E= -0.039999992 adding state (0, -0.9999999428189236) with fermi= 1.0 adding state (0, -0.2499999899816187) with fermi= 1.0 adding state (1, -0.2499999979714578) with fermi= 1.0 adding state (0, -0.11111110835536761) with fermi= 1.0 adding state (1, -0.11111111057681711) with fermi= 1.0 adding state (2, -0.1111111111469023) with fermi= 1.0
Hartree term¶
The Hartree term is treated exactly in this approximation.
It describes the electrostatic interaction of one electron with the cloud of all electrons (including the electron itself).
This is classical electrostatic term. Namely, any electron feels not only the nucleous with Coulomb attraction $-2Z/r$, but also the charge of all (other) electrons $\rho(\vec{r}')$ through the same Coulomb force, but it is repulsive, hence we should have $$-2Z/r \rightarrow -2Z/r + \int d\vec{r}' \frac{2\rho(\vec{r}')}{|\vec{r}-\vec{r}'|}$$
Mathematically, this term comes from the mean-field approximation (lowest order) of the Coulomb repulsion between electrons: \begin{eqnarray} && \frac{1}{2}\int d\vec{r} d\vec{r}' \psi^\dagger (\vec{r})\psi^\dagger (\vec{r}') v_c(\vec{r}-\vec{r}') \psi(\vec{r}')\psi(\vec{r}) \rightarrow\\ && \int d\vec{r} \psi^\dagger(\vec{r}) \psi(\vec{r}) \int d\vec{r}' \langle\psi^\dagger(\vec{r}') \psi(\vec{r}')\rangle v_c(\vec{r}-\vec{r}') \equiv \int d\vec{r} \psi^\dagger(\vec{r}) V_{H}(\vec{r}) \psi(\vec{r})\nonumber \end{eqnarray} with \begin{equation} V_H(\vec{r}) = 2 \int d\vec{r}' \frac{\rho(\vec{r}')}{|\vec{r}-\vec{r}'|} \end{equation} where $2$ is due to Rydberg units sinc $v_c = 2/r$.
For any atom, the electron density is spherically symetric and hence $V_{H}$ depends only on radial distance. (In solids, the hartree potential should be expanded in spheric harmonics to sufficiently high $l$, maybe $l=6$).
Step 1: Using $\rho(r)$ computed in previous homework, compute the Hartree potential.¶
This is usually achieved by solving the Poisson equation. From clasical electrostatic we know \begin{eqnarray} \nabla^2 V_{H}(\vec{r}) = -8\pi \rho(\vec{r}) \end{eqnarray} If we express $\nabla$ in spherical coordinates, we get: \begin{equation} \frac{1}{r^2}\frac{d}{dr}(r^2 \frac{d V_H}{dr})= -8\pi\rho(r) \end{equation} which simplifies to \begin{equation} U^{''}(r) = -8\pi r \rho(r) \end{equation} where $U(r) = V_{H}(r) r$.
This second order differential equation has the following boundary conditions $U(0)=0$ and $U(\infty)=2 Z$.
The two point boundary problem does not require shooting because we know solution to the homogenous differential equation $U^{''}(r)=0$. The Hartree potential can be obtained from any particular solution by \begin{equation} U(r) = U_p(r) + \alpha r \end{equation} where $\alpha = \lim_{r\rightarrow\infty}(2 Z-U_{p}(r))/r$.
def FuncForHartree(y,r,rhoSpline):
""" y = [U,U']
dy/dr = [U', -8*pi*r*rho(r)]
"""
return [y[1], -8*pi*r*rhoSpline(r)]
from scipy import interpolate
rhoSpline = interpolate.UnivariateSpline(R, rho,s=0)
U0 = integrate.odeint(FuncForHartree, [0.0,0.5], R, args=(rhoSpline,))[:,0]
alpha = (2*Z - U0[-1])/R[-1]
U0 += alpha * R
plot(R, U0,label='U(r)')
plot(R, ones(len(R))*2*Z, label='2Z')
legend(loc='best')
show()
Numerov again¶
Poisson equation does not have the first order derivative, hence it can also be more efficiently solved by the Numerov algorithm.
We have Poisson equation, which has the form \begin{equation} x^{''}(t)= u(t) \end{equation} and the Numerov algorithm, as appropriate for the Poisson equation, is \begin{eqnarray} x(h)+x(-h) = 2x(0)+h^2 u(0)+\frac{2}{4!}h^4 x^{(4)}(0)+O(h^6) \end{eqnarray} and the approximation for the forth order derivative is \begin{equation} x^{(4)}\sim \frac{u_{i+1}-2 u_i+u_{i-1}}{h^2} \end{equation}
Inserting the fourth order derivative into the above recursive equation (forth equation in his chapter), we get
\begin{equation} x_{i+1}-2 x_i+x_{i-1}=h^2 u_i +\frac{h^2}{12}(u_{i+1}-2 u_i+u_{i-1}) \end{equation}If we switch to a new variable $w_i=x_i-\frac{h^2}{12}u_i$ we are left with the following equation
\begin{equation} w_{i+1} -2 w_i + w_{i-1} = h^2 u_i+O(h^6) \end{equation}The variable $x$ needs to be recomputed at each step with $x_i=(w_i+\frac{h^2}{12}u_i)$.
@jit(nopython=True)
def NumerovUP(U, x0, dx, dh):
x = zeros(len(U))
x[0] = x0
x[1] = dx*dh + x0
h2 = dh*dh
h12 = h2/12
w0 = x[0]-h12*U[0]
w1 = x[1]-h12*U[1]
Ui = U[1]
for i in range(2,len(U)):
w2 = 2*w1 - w0 + h2*Ui
Ui = U[i]
xi = w2 + h12*Ui
x[i] = xi
w0, w1 = w1, w2
return x
U2 = NumerovUP(-8*pi*R*rho, 0.0, 0.5, R[1]-R[0])
alpha2 = (2*Z-U2[-1])/R[-1]
U2 += alpha2 * R
plot(R,U0)
plot(R,U2);
Step1 routine HartreeU¶
For generic density the following routine will work:
def HartreeU(R, rho, Zatom):
"""Given input charge density it returns Hartree potential in the form VH(r)*r
"""
U2 = NumerovUP(-8*pi*R*rho, 0.0, 0.5, R[1]-R[0])
alpha2 = (2*Zatom-U2[-1])/R[-1]
U2 += alpha2 * R
return U2
U2 = HartreeU(R,rho,Z)
plot(R,U2);
Step 2 : Compute the exchange correlation potential.¶
Note that $V_{xc}(r)=V_{xc}(\rho(r))$ is unquely determined by the electron charge density $\rho(r)$. If we know $\rho$, we can instantly compute $V_{xc}$ by the module provided parametrized function.
Download the module excor.py from
http://www.physics.rutgers.edu/~haule/509/src_prog/python/homework5/
and import it in your code.
Instantiate the ExchangeCorrelation object by
exc = ExchangeCorrelation()
and used it, for example, by
exc.Vx(rs(rho[i]))+exc.Vc(rs(rho[i]))
where $r_s = ({4\pi\rho/3})^{-1/3}$.
Be careful: The energy unit in "excor.py" is Hartree and not Rydergs. Hence, you need to multiply energy or potential by 2.
from excor import ExchangeCorrelation
exc = ExchangeCorrelation()
@jit(nopython=True)
def rs(rho):
"1/rho = 4*pi*rs^3/3 => rs = (3/(4*pi*rho))**(1/3.)"
if rho < 1e-100: return 1e100
return pow(3/(4*pi*rho),1/3.)
Step 3: Find bound states using Hartree and XC potential.¶
Add the Hartree potential and the exchange-correlation potential to the Schroedinger equation and find bound states of the atom.
The Schroedinger equation is \begin{equation} u^{''}(r) = \left(\frac{l(l+1)}{r^2}-\frac{2 Z}{r} + V_{H}(r)+V_{XC}(r)-\varepsilon\right)u(r). \end{equation} or \begin{equation} u^{''}(r) = \left(\frac{l(l+1)}{r^2}+\frac{U_H - 2 Z +r V_{XC}}{r}-\varepsilon\right)u(r). \end{equation}
We will use notation: $$U_K=U_H - 2 Z +r V_{XC}$$ so that \begin{equation} u^{''}(r) = \left(\frac{l(l+1)}{r^2}+\frac{U_K}{r}-\varepsilon\right)u(r). \end{equation}
Vxc = array([2*exc.Vc(rs(rh)) + 2*exc.Vx(rs(rh)) for rh in rho]) # factor of 2 due to change of units from Hartree to Rydbers
Vx = array([2*exc.Vx(rs(rh)) for rh in rho])
Vc = array([2*exc.Vc(rs(rh)) for rh in rho])
Uks = U2 - 2*Z + Vxc*R
fig,(ax1, ax2) = plt.subplots(2, 1)
ax1.plot(R, -Uks, label='-U_total');
ax1.plot(R, -U2, label='-Uh')
ax1.plot(R, -U2+2*Z, label='-Uh+2Z')
ax2.plot(R,-Vx*R,label='-Ux')
ax2.plot(R,-Vc*R,label='-Uc')
#ax2.plot(R,U2,label='Uh')
ax2.legend(loc='best')
ax1.legend(loc='best')
ax2.grid()
The effective potential $U_{KS}== r(V_H+V_{xc}-2 Z/r)$ has the following features:
- at small $r$ it is just bare nucleous potential $U_{KS}=-2 Z$ because $V_{KS}=-2 Z/r$
- at large $r$ the screening of nuclear charge is perfect, and $U_{KS}=0$ exponentially fast beyond the charge clound, hence $U_{KS}\approx 2(\rho_{inside}-Z)$ vanishes exponentially.
@jit(nopython=True)
def fShrod(En,l,R, Uks):
return (l*(l+1.)/R +Uks)/R - En
def ComputeSchrod(En, R, l, Uks):
#f = fShrod(En,l,R[::-1],Uks[::-1])
#ur = Numerov(f, 0.0, -1e-10, R[0]-R[1])[::-1]
f = fShrod(En,l,R, Uks)
ur = Numerov(f[::-1], 0.0, -1e-10, R[0]-R[1])[::-1]
norm = integrate.simps(ur**2, x=R)
return ur/sqrt(abs(norm))
def Shoot(En, R, l, Uks):
ur = ComputeSchrod(En, R, l,Uks)
ur *= 1/R**l # expecting urn \propto R
#f0,f1 = ur[0],ur[1]
#f_at_0 = f0 + (f1-f0)*(0-R[0])/(R[1]-R[0]) # extrapolation to zero
#return f_at_0
poly = Polynomial.fit(R[:4], ur[:4], deg=3)
return poly(0.0)
def FindBoundStates(R, l, nmax, Esearch, Uks):
n=0
Ebnd=[]
u0 = Shoot(Esearch[0],R,l,Uks)
for i in range(1,len(Esearch)):
u1 = Shoot(Esearch[i],R,l,Uks)
if u0*u1 < 0:
Ebound = optimize.brentq(Shoot,Esearch[i-1],Esearch[i],xtol=1e-15,args=(R,l,Uks))
Ebnd.append( (l,Ebound) )
if len(Ebnd)>nmax: break
n += 1
print(f"Found bound state at E={Ebound/2:14.9f}H l={l:2d}")
u0 = u1
return Ebnd
## From before we are using
# R = linspace(1e-6,100,2000)
# Esearch = -1.2/arange(1,20,0.2)**2
nmax=5
Bnd=[]
for l in range(nmax-1):
Bnd += FindBoundStates(R,l,nmax-l,Esearch,Uks)
Bnd = sorted(Bnd, key=cmpKey)
Found bound state at E= -0.301181700H l= 0 Found bound state at E= -0.114301424H l= 0 Found bound state at E= -0.034082204H l= 0 Found bound state at E= -0.003612354H l= 0 Found bound state at E= -0.313798778H l= 1 Found bound state at E= -0.118616033H l= 1 Found bound state at E= -0.034746947H l= 1 Found bound state at E= -0.003345161H l= 1 Found bound state at E= -0.300009838H l= 2 Found bound state at E= -0.108653489H l= 2 Found bound state at E= -0.028845107H l= 2 Found bound state at E= -0.279116709H l= 3 Found bound state at E= -0.093722034H l= 3
Step 4: Compute the new electron density¶
by filling the lowest $Z$ eigenstatates.
# This is modified from Hydrogen
def ChargeDensity(bst,R,Zatom,Uks):
rho = zeros(len(R))
N=0.
Ebs=0. # sum of all eigenvalues of KS equation.
for (l,En) in bst:
ur = ComputeSchrod(En, R, l, Uks)
dN = 2*(2*l+1)
if N+dN <= Zatom:
ferm = 1.
else:
ferm = (Zatom-N)/float(dN)
drho = ur**2 * ferm * dN/(4*pi*R**2)
rho += drho
N += dN
Ebs += En * dN * ferm
print('adding state', (l,En/2), 'H with fermi=', ferm)
if N>=Zatom: break
return (rho,Ebs)
Z=28
rho_new, Ebs = ChargeDensity(Bnd,R,Z,Uks)
adding state (1, -0.3137987779463955) H with fermi= 1.0 adding state (0, -0.30118170042796616) H with fermi= 1.0 adding state (2, -0.30000983841183604) H with fermi= 1.0 adding state (3, -0.27911670884400286) H with fermi= 0.7142857142857143
Evaluate the total energy¶
Once we see that the code converges, we can insert a new step betwen Step 4 and Step 5 to compute the total energy of the system. The total energy can be obtained by \begin{eqnarray} E^{LDA}_{total} &=& \sum_{i\in occupied}\int d\vec{r} \psi_i^*(\vec{r})[-\nabla^2]\psi_i(\vec{r}) +\nonumber\\ &+& \int d\vec{r} \rho(\vec{r}) [V_{nucleous}(\vec{r})+\epsilon_H(\vec{r}) + \epsilon_{XC}(\vec{r})]\nonumber\\ &=& \sum_{i\in occupied}\int d\vec{r} \psi_i^*(\vec{r})[-\nabla^2+V_{nucleous}+V_H+V_{XC}]\psi_i(\vec{r}) \nonumber\\ &+& \int d\vec{r} \rho(\vec{r}) [\epsilon_H(\vec{r})-V_H(\vec{r}) + \epsilon_{XC}(\vec{r})-V_{XC}(\vec{r})]\nonumber\\ &=& \sum_{i\in occupied}\epsilon_i + \int d\vec{r} \rho(\vec{r}) [\epsilon_H(\vec{r})-V_H(\vec{r}) + \epsilon_{XC}(\vec{r})-V_{XC}(\vec{r})]\nonumber\\ &=& \sum_{i\in occupied}\epsilon_i + \int d\vec{r} \rho(\vec{r}) [-\epsilon_H(\vec{r}) + \epsilon_{XC}(\vec{r})-V_{XC}(\vec{r})] \end{eqnarray} Here we used \begin{eqnarray} && E_H[\rho] \equiv \int d\vec{r}\; \rho(\vec{r})\; \epsilon_H[\rho(\vec{r})]\\ && E_{xc}[\rho] \equiv \int d\vec{r}\; \rho(\vec{r})\; \epsilon_{xc}[\rho(\vec{r})]\\ && V_H[\rho]\equiv \frac{\delta E_H[\rho]}{\delta \rho(\vec{r})}\\ && V_{xc}[\rho]\equiv \frac{\delta E_{xc}[\rho]}{\delta \rho(\vec{r})}. \end{eqnarray}
Here the Hartree energy is the classical electrostatic potential energy, i.e., \begin{equation} E_H = \int d\vec{r} d\vec{r}' \frac{\rho(\vec{r})\rho(\vec{r}')}{|\vec{r}-\vec{r}'|}, \end{equation} hence \begin{eqnarray} \epsilon_H(\vec{r}) = \frac{1}{2} V_H(\vec{r}) = \int d\vec{r}' \frac{\rho(\vec{r}')}{|\vec{r}-\vec{r}'|}\\ \end{eqnarray} and we also know \begin{eqnarray} V_H(\vec{r})= \frac{U_H(\vec{r})}{r} \end{eqnarray} Hence $$\epsilon_H(\vec{r})=\frac{1}{2}\frac{U_H}{r}$$
The exchange-correlation energy can be obtained by a call to the method of ${ExchangeCorrelation}$ object.
For oxygen, the total energy in this implementation is : -74.4730769 Hartree, while NIST shows -74.473077 Hartree. This is in excellent agreement in all digits provided by NIST.
Zatom=8
E0=-1.2*Zatom**2
Eshift=0.5 # sometimes energies can be positive!!!
Esearch = -logspace(-4,log10(-E0+Eshift),200)[::-1] + Eshift
plot(Esearch,'o-')
ylim([-1.2*Zatom**2,1])
(-76.8, 1.0)
We will first demonstrate algorithm on oxygen $Z=8$ and than on Potassium ($Z=19$), which has first nontrivial "aufbau" principle in which $4s$ orbital gets filled before $3d$.
R = linspace(1e-8,20,2**14+1)
Zatom = 20 # 20 # 19
mixr = 0.5
E0=-1.2*Zatom**2
Eshift=0.5 # sometimes energies can be positive!!!
Esearch = -logspace(-4,log10(-E0+Eshift),200)[::-1] + Eshift
nmax = 4 #3
exc = ExchangeCorrelation()
Uks = -2*ones(len(R))
Eold = 0
Etol = 1e-7
for itt in range(1000):
Bnd=[]
for l in range(nmax-1):
Bnd += FindBoundStates(R,l,nmax-l,Esearch,Uks)
Bnd = sorted(Bnd, key=cmpKey)
rho_new, Ebs = ChargeDensity(Bnd,R,Zatom,Uks)
if itt==0:
rho = rho_new
else:
rho = rho_new * mixr + (1-mixr)*rho_old
rho_old = copy(rho)
U2 = HartreeU(R, rho, Zatom)
Vxc = [2*exc.Vc(rs(rh)) + 2*exc.Vx(rs(rh)) for rh in rho]
Uks = U2 - 2*Zatom + Vxc*R
# Total energy
ExcVxc = [2*exc.EcVc(rs(rh)) + 2*exc.ExVx(rs(rh)) for rh in rho] # eps_xc(rho)-V_xc(rho)
pot = (ExcVxc*R - 0.5*U2)*R*rho*4*pi # (eps_xc-V_xc-0.5 U_H/R)*rho * d^3r
epot = integrate.simps(pot, x=R)
Etot = epot + Ebs
print('Total density has weight', integrate.simps(rho*(4*pi*R**2),x=R))
#print('Total Energy=', Etot/2.)
print('Itteration', itt, 'Etot[Ry]=', Etot, 'Etot[Hartre]=', Etot/2, 'Diff=', abs(Etot-Eold))
if itt>0 and abs(Etot-Eold)< Etol: break
Eold = Etot
#plot(R, U2, label='U-hartree')
#plot(R, Vxc, label='Vxc')
#plot(R, Uks, label='Uks')
#show()
plot(R,rho*(4*pi*R**2))
xlim([0,5])
show()
/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/ipykernel_19754/3467380318.py:10: RuntimeWarning: overflow encountered in square norm = integrate.simps(ur**2, x=R) /Users/haule/anaconda3/lib/python3.11/site-packages/scipy/integrate/_quadrature.py:514: RuntimeWarning: overflow encountered in multiply y[slice1] * (hsum * /Users/haule/anaconda3/lib/python3.11/site-packages/scipy/integrate/_quadrature.py:510: RuntimeWarning: overflow encountered in add tmp = hsum/6.0 * (y[slice0] *
Found bound state at E= -0.500000000H l= 0 Found bound state at E= -0.124987114H l= 0 Found bound state at E= -0.049918048H l= 0 Found bound state at E= 0.016713154H l= 0 Found bound state at E= -0.124994607H l= 1 Found bound state at E= -0.051611420H l= 1 Found bound state at E= 0.008086683H l= 1 Found bound state at E= -0.053967564H l= 2 Found bound state at E= -0.005535033H l= 2 adding state (0, -0.4999999999999946) H with fermi= 1.0 adding state (0, -0.12498711431297779) H with fermi= 1.0 adding state (1, -0.12499460664692982) H with fermi= 1.0 adding state (2, -0.05396756442457235) H with fermi= 1.0 Total density has weight 20.0 Itteration 0 Etot[Ry]= -56.46605274547102 Etot[Hartre]= -28.23302637273551 Diff= 56.46605274547102 Found bound state at E= -45.710175753H l= 0 Found bound state at E= -18.036813308H l= 0 Found bound state at E= -8.533495295H l= 0 Found bound state at E= -4.311115427H l= 0 Found bound state at E= -45.707060807H l= 1 Found bound state at E= -18.019757254H l= 1 Found bound state at E= -8.513487979H l= 1 Found bound state at E= -17.981446754H l= 2 Found bound state at E= -8.470663757H l= 2 adding state (0, -45.71017575302781) H with fermi= 1.0 adding state (1, -45.707060807478015) H with fermi= 1.0 adding state (0, -18.036813308357747) H with fermi= 1.0 adding state (1, -18.01975725379621) H with fermi= 1.0 adding state (2, -17.981446754405038) H with fermi= 0.4 Total density has weight 20.0 Itteration 1 Etot[Ry]= -1415.9009163025169 Etot[Hartre]= -707.9504581512584 Diff= 1359.434863557046 Found bound state at E=-170.041614028H l= 0 Found bound state at E= -26.493156166H l= 0 Found bound state at E= -7.052823895H l= 0 Found bound state at E= -2.567760463H l= 0 Found bound state at E= -25.016519856H l= 1 Found bound state at E= -6.410905650H l= 1 Found bound state at E= -2.320880181H l= 1 Found bound state at E= -5.097650320H l= 2 Found bound state at E= -1.863907428H l= 2 adding state (0, -170.04161402792897) H with fermi= 1.0 adding state (0, -26.49315616633762) H with fermi= 1.0 adding state (1, -25.01651985620556) H with fermi= 1.0 adding state (0, -7.052823895041433) H with fermi= 1.0 adding state (1, -6.410905650217471) H with fermi= 1.0 adding state (2, -5.0976503197204455) H with fermi= 0.2 Total density has weight 20.000000000000004 Itteration 2 Etot[Ry]= -1674.0272758080575 Etot[Hartre]= -837.0136379040288 Diff= 258.1263595055407 Found bound state at E=-152.142793735H l= 0 Found bound state at E= -17.426151876H l= 0 Found bound state at E= -2.916304329H l= 0 Found bound state at E= -0.749800998H l= 0 Found bound state at E= -14.956211356H l= 1 Found bound state at E= -2.219502402H l= 1 Found bound state at E= -0.563955035H l= 1 Found bound state at E= -1.065212262H l= 2 Found bound state at E= -0.279039611H l= 2 adding state (0, -152.1427937348216) H with fermi= 1.0 adding state (0, -17.426151875846312) H with fermi= 1.0 adding state (1, -14.956211356444634) H with fermi= 1.0 adding state (0, -2.916304328724621) H with fermi= 1.0 adding state (1, -2.219502401693711) H with fermi= 1.0 adding state (2, -1.0652122615217565) H with fermi= 0.2 Total density has weight 20.000000000000004 Itteration 3 Etot[Ry]= -1435.425852319622 Etot[Hartre]= -717.712926159811 Diff= 238.60142348843556 Found bound state at E=-146.419037947H l= 0 Found bound state at E= -15.166166955H l= 0 Found bound state at E= -1.783564720H l= 0 Found bound state at E= -0.282486268H l= 0 Found bound state at E= -12.477382815H l= 1 Found bound state at E= -1.130102235H l= 1 Found bound state at E= -0.168846203H l= 1 Found bound state at E= -0.214548236H l= 2 Found bound state at E= -0.039867014H l= 2 adding state (0, -146.41903794733756) H with fermi= 1.0 adding state (0, -15.166166954651871) H with fermi= 1.0 adding state (1, -12.477382814556874) H with fermi= 1.0 adding state (0, -1.7835647204889968) H with fermi= 1.0 adding state (1, -1.1301022351447578) H with fermi= 1.0 adding state (0, -0.28248626836652746) H with fermi= 1.0 Total density has weight 20.0 Itteration 4 Etot[Ry]= -1358.7724996174811 Etot[Hartre]= -679.3862498087406 Diff= 76.65335270214086 Found bound state at E=-145.089884322H l= 0 Found bound state at E= -15.140545494H l= 0 Found bound state at E= -1.768949616H l= 0 Found bound state at E= -0.209342490H l= 0 Found bound state at E= -12.404203037H l= 1 Found bound state at E= -1.100086618H l= 1 Found bound state at E= -0.107204348H l= 1 Found bound state at E= -0.155906703H l= 2 Found bound state at E= -0.013627680H l= 2 adding state (0, -145.08988432238033) H with fermi= 1.0 adding state (0, -15.140545493702785) H with fermi= 1.0 adding state (1, -12.404203037263269) H with fermi= 1.0 adding state (0, -1.7689496162816907) H with fermi= 1.0 adding state (1, -1.100086618357816) H with fermi= 1.0 adding state (0, -0.20934249001468003) H with fermi= 1.0 Total density has weight 20.0 Itteration 5 Etot[Ry]= -1355.930837575852 Etot[Hartre]= -677.965418787926 Diff= 2.841662041629206 Found bound state at E=-144.447668773H l= 0 Found bound state at E= -15.086599507H l= 0 Found bound state at E= -1.735001622H l= 0 Found bound state at E= -0.169950130H l= 0 Found bound state at E= -12.332827823H l= 1 Found bound state at E= -1.061714708H l= 1 Found bound state at E= -0.076208196H l= 1 Found bound state at E= -0.114761338H l= 2 Found bound state at E= 0.000093512H l= 2 adding state (0, -144.44766877287543) H with fermi= 1.0 adding state (0, -15.086599507340976) H with fermi= 1.0 adding state (1, -12.332827823155622) H with fermi= 1.0 adding state (0, -1.735001621525317) H with fermi= 1.0 adding state (1, -1.061714708028648) H with fermi= 1.0 adding state (0, -0.16995012996904738) H with fermi= 1.0 Total density has weight 20.0 Itteration 6 Etot[Ry]= -1353.4539129621244 Etot[Hartre]= -676.7269564810622 Diff= 2.476924613727533 Found bound state at E=-144.162123936H l= 0 Found bound state at E= -15.061995444H l= 0 Found bound state at E= -1.718218090H l= 0 Found bound state at E= -0.153068519H l= 0 Found bound state at E= -12.302450721H l= 1 Found bound state at E= -1.043368218H l= 1 Found bound state at E= -0.063033548H l= 1 Found bound state at E= -0.096118259H l= 2 Found bound state at E= 0.008005757H l= 2 adding state (0, -144.162123936473) H with fermi= 1.0 adding state (0, -15.061995443504403) H with fermi= 1.0 adding state (1, -12.302450720662986) H with fermi= 1.0 adding state (0, -1.718218089516831) H with fermi= 1.0 adding state (1, -1.0433682176538581) H with fermi= 1.0 adding state (0, -0.15306851868308669) H with fermi= 1.0 Total density has weight 19.999999999999996 Itteration 7 Etot[Ry]= -1352.309105992004 Etot[Hartre]= -676.154552996002 Diff= 1.1448069701202712 Found bound state at E=-144.036444793H l= 0 Found bound state at E= -15.052483847H l= 0 Found bound state at E= -1.711176758H l= 0 Found bound state at E= -0.146168310H l= 0 Found bound state at E= -12.291214371H l= 1 Found bound state at E= -1.035745940H l= 1 Found bound state at E= -0.057528874H l= 1 Found bound state at E= -0.088404182H l= 2 Found bound state at E= 0.012986556H l= 2 adding state (0, -144.03644479323992) H with fermi= 1.0 adding state (0, -15.052483846511437) H with fermi= 1.0 adding state (1, -12.291214370912371) H with fermi= 1.0 adding state (0, -1.711176757828648) H with fermi= 1.0 adding state (1, -1.0357459400474252) H with fermi= 1.0 adding state (0, -0.1461683095483543) H with fermi= 1.0 Total density has weight 20.0 Itteration 8 Etot[Ry]= -1351.8288468345079 Etot[Hartre]= -675.9144234172539 Diff= 0.4802591574962207 Found bound state at E=-143.980772259H l= 0 Found bound state at E= -15.048940338H l= 0 Found bound state at E= -1.708311167H l= 0 Found bound state at E= -0.143367233H l= 0 Found bound state at E= -12.287257061H l= 1 Found bound state at E= -1.032666096H l= 1 Found bound state at E= -0.055209816H l= 1 Found bound state at E= -0.085264899H l= 2 Found bound state at E= 0.016256265H l= 2 adding state (0, -143.980772258525) H with fermi= 1.0 adding state (0, -15.048940338478968) H with fermi= 1.0 adding state (1, -12.287257061334225) H with fermi= 1.0 adding state (0, -1.7083111668731947) H with fermi= 1.0 adding state (1, -1.0326660958019334) H with fermi= 1.0 adding state (0, -0.14336723271393864) H with fermi= 1.0 Total density has weight 20.0 Itteration 9 Etot[Ry]= -1351.629568738381 Etot[Hartre]= -675.8147843691905 Diff= 0.19927809612681813 Found bound state at E=-143.955883561H l= 0 Found bound state at E= -15.047636754H l= 0 Found bound state at E= -1.707145905H l= 0 Found bound state at E= -0.142224887H l= 0 Found bound state at E= -12.285920339H l= 1 Found bound state at E= -1.031424104H l= 1 Found bound state at E= -0.054216255H l= 1 Found bound state at E= -0.083984480H l= 2 Found bound state at E= 0.018435478H l= 2 adding state (0, -143.95588356126763) H with fermi= 1.0 adding state (0, -15.0476367537906) H with fermi= 1.0 adding state (1, -12.285920339470987) H with fermi= 1.0 adding state (0, -1.7071459045565336) H with fermi= 1.0 adding state (1, -1.031424103865588) H with fermi= 1.0 adding state (0, -0.14222488709815914) H with fermi= 1.0 Total density has weight 20.0 Itteration 10 Etot[Ry]= -1351.5462801680246 Etot[Hartre]= -675.7731400840123 Diff= 0.08328857035644432 Found bound state at E=-143.944656717H l= 0 Found bound state at E= -15.047162989H l= 0 Found bound state at E= -1.706669430H l= 0 Found bound state at E= -0.141754849H l= 0 Found bound state at E= -12.285497565H l= 1 Found bound state at E= -1.030921424H l= 1 Found bound state at E= -0.053780746H l= 1 Found bound state at E= -0.083458401H l= 2 Found bound state at E= 0.019884740H l= 2 adding state (0, -143.94465671651386) H with fermi= 1.0 adding state (0, -15.047162988513845) H with fermi= 1.0 adding state (1, -12.285497565148965) H with fermi= 1.0 adding state (0, -1.7066694298935265) H with fermi= 1.0 adding state (1, -1.0309214243276112) H with fermi= 1.0 adding state (0, -0.14175484941232647) H with fermi= 1.0 Total density has weight 19.999999999999996 Itteration 11 Etot[Ry]= -1351.5111204412322 Etot[Hartre]= -675.7555602206161 Diff= 0.0351597267924717 Found bound state at E=-143.939550077H l= 0 Found bound state at E= -15.046993699H l= 0 Found bound state at E= -1.706473072H l= 0 Found bound state at E= -0.141559025H l= 0 Found bound state at E= -12.285380031H l= 1 Found bound state at E= -1.030716803H l= 1 Found bound state at E= -0.053584560H l= 1 Found bound state at E= -0.083240030H l= 2 Found bound state at E= 0.020835015H l= 2 adding state (0, -143.93955007661577) H with fermi= 1.0 adding state (0, -15.046993698637914) H with fermi= 1.0 adding state (1, -12.28538003095776) H with fermi= 1.0 adding state (0, -1.7064730720257413) H with fermi= 1.0 adding state (1, -1.030716803292129) H with fermi= 1.0 adding state (0, -0.141559025215455) H with fermi= 1.0 Total density has weight 19.999999999999996 Itteration 12 Etot[Ry]= -1351.4961252860926 Etot[Hartre]= -675.7480626430463 Diff= 0.014995155139558847 Found bound state at E=-143.937209305H l= 0 Found bound state at E= -15.046934585H l= 0 Found bound state at E= -1.706391320H l= 0 Found bound state at E= -0.141476064H l= 0 Found bound state at E= -12.285356961H l= 1 Found bound state at E= -1.030632832H l= 1 Found bound state at E= -0.053493540H l= 1 Found bound state at E= -0.083148140H l= 2 Found bound state at E= 0.021444129H l= 2 adding state (0, -143.9372093050684) H with fermi= 1.0 adding state (0, -15.046934585371533) H with fermi= 1.0 adding state (1, -12.285356961276797) H with fermi= 1.0 adding state (0, -1.7063913197773042) H with fermi= 1.0 adding state (1, -1.0306328320104254) H with fermi= 1.0 adding state (0, -0.14147606359097104) H with fermi= 1.0 Total density has weight 19.999999999999996 Itteration 13 Etot[Ry]= -1351.4896615303428 Etot[Hartre]= -675.7448307651714 Diff= 0.006463755749791744 Found bound state at E=-143.936128738H l= 0 Found bound state at E= -15.046914582H l= 0 Found bound state at E= -1.706356848H l= 0 Found bound state at E= -0.141440163H l= 0 Found bound state at E= -12.285358885H l= 1 Found bound state at E= -1.030598000H l= 1 Found bound state at E= -0.053450096H l= 1 Found bound state at E= -0.083108804H l= 2 Found bound state at E= 0.021823476H l= 2 adding state (0, -143.93612873771556) H with fermi= 1.0 adding state (0, -15.046914581816887) H with fermi= 1.0 adding state (1, -12.285358885080857) H with fermi= 1.0 adding state (0, -1.706356847560514) H with fermi= 1.0 adding state (1, -1.030597999754475) H with fermi= 1.0 adding state (0, -0.14144016312558877) H with fermi= 1.0 Total density has weight 20.0 Itteration 14 Etot[Ry]= -1351.4868442543043 Etot[Hartre]= -675.7434221271521 Diff= 0.002817276038513228 Found bound state at E=-143.935626695H l= 0 Found bound state at E= -15.046908096H l= 0 Found bound state at E= -1.706342085H l= 0 Found bound state at E= -0.141424224H l= 0 Found bound state at E= -12.285364846H l= 1 Found bound state at E= -1.030583348H l= 1 Found bound state at E= -0.053428832H l= 1 Found bound state at E= -0.083091610H l= 2 Found bound state at E= 0.022052098H l= 2 adding state (0, -143.93562669515623) H with fermi= 1.0 adding state (0, -15.046908096120836) H with fermi= 1.0 adding state (1, -12.285364846415135) H with fermi= 1.0 adding state (0, -1.706342084797461) H with fermi= 1.0 adding state (1, -1.030583347565921) H with fermi= 1.0 adding state (0, -0.14142422442489025) H with fermi= 1.0 Total density has weight 20.0 Itteration 15 Etot[Ry]= -1351.4856018962191 Etot[Hartre]= -675.7428009481096 Diff= 0.0012423580851645966 Found bound state at E=-143.935392086H l= 0 Found bound state at E= -15.046906120H l= 0 Found bound state at E= -1.706335652H l= 0 Found bound state at E= -0.141416945H l= 0 Found bound state at E= -12.285369648H l= 1 Found bound state at E= -1.030577082H l= 1 Found bound state at E= -0.053418209H l= 1 Found bound state at E= -0.083083919H l= 2 Found bound state at E= 0.022185195H l= 2 adding state (0, -143.93539208589306) H with fermi= 1.0 adding state (0, -15.046906120300644) H with fermi= 1.0 adding state (1, -12.285369648242195) H with fermi= 1.0 adding state (0, -1.7063356523907418) H with fermi= 1.0 adding state (1, -1.030577081589071) H with fermi= 1.0 adding state (0, -0.141416944835085) H with fermi= 1.0 Total density has weight 20.0 Itteration 16 Etot[Ry]= -1351.485047645418 Etot[Hartre]= -675.742523822709 Diff= 0.000554250801087619 Found bound state at E=-143.935281879H l= 0 Found bound state at E= -15.046905571H l= 0 Found bound state at E= -1.706332794H l= 0 Found bound state at E= -0.141413520H l= 0 Found bound state at E= -12.285372689H l= 1 Found bound state at E= -1.030574348H l= 1 Found bound state at E= -0.053412814H l= 1 Found bound state at E= -0.083080389H l= 2 Found bound state at E= 0.022260123H l= 2 adding state (0, -143.93528187868924) H with fermi= 1.0 adding state (0, -15.046905570543528) H with fermi= 1.0 adding state (1, -12.285372689393382) H with fermi= 1.0 adding state (0, -1.7063327936724986) H with fermi= 1.0 adding state (1, -1.0305743479030367) H with fermi= 1.0 adding state (0, -0.14141351957391055) H with fermi= 1.0 Total density has weight 20.0 Itteration 17 Etot[Ry]= -1351.4847974405957 Etot[Hartre]= -675.7423987202978 Diff= 0.0002502048223504971 Found bound state at E=-143.935229868H l= 0 Found bound state at E= -15.046905436H l= 0 Found bound state at E= -1.706331495H l= 0 Found bound state at E= -0.141411859H l= 0 Found bound state at E= -12.285374422H l= 1 Found bound state at E= -1.030573127H l= 1 Found bound state at E= -0.053410037H l= 1 Found bound state at E= -0.083078727H l= 2 Found bound state at E= 0.022301084H l= 2 adding state (0, -143.9352298679875) H with fermi= 1.0 adding state (0, -15.046905436450055) H with fermi= 1.0 adding state (1, -12.285374421851651) H with fermi= 1.0 adding state (0, -1.7063314946863342) H with fermi= 1.0 adding state (1, -1.0305731267904832) H with fermi= 1.0 adding state (0, -0.14141185919570198) H with fermi= 1.0 Total density has weight 20.0 Itteration 18 Etot[Ry]= -1351.4846830785493 Etot[Hartre]= -675.7423415392747 Diff= 0.00011436204636083858 Found bound state at E=-143.935205224H l= 0 Found bound state at E= -15.046905411H l= 0 Found bound state at E= -1.706330893H l= 0 Found bound state at E= -0.141411035H l= 0 Found bound state at E= -12.285375354H l= 1 Found bound state at E= -1.030572570H l= 1 Found bound state at E= -0.053408595H l= 1 Found bound state at E= -0.083077926H l= 2 Found bound state at E= 0.022322974H l= 2 adding state (0, -143.93520522387615) H with fermi= 1.0 adding state (0, -15.046905411356422) H with fermi= 1.0 adding state (1, -12.285375353635676) H with fermi= 1.0 adding state (0, -1.7063308933670132) H with fermi= 1.0 adding state (1, -1.0305725698723494) H with fermi= 1.0 adding state (0, -0.14141103521402812) H with fermi= 1.0 Total density has weight 20.0 Itteration 19 Etot[Ry]= -1351.4846303856943 Etot[Hartre]= -675.7423151928472 Diff= 5.2692854978886317e-05 Found bound state at E=-143.935193501H l= 0 Found bound state at E= -15.046905403H l= 0 Found bound state at E= -1.706330605H l= 0 Found bound state at E= -0.141410613H l= 0 Found bound state at E= -12.285375831H l= 1 Found bound state at E= -1.030572305H l= 1 Found bound state at E= -0.053407836H l= 1 Found bound state at E= -0.083077527H l= 2 Found bound state at E= 0.022334499H l= 2 adding state (0, -143.93519350074973) H with fermi= 1.0 adding state (0, -15.046905403494092) H with fermi= 1.0 adding state (1, -12.285375831264126) H with fermi= 1.0 adding state (0, -1.7063306047397166) H with fermi= 1.0 adding state (1, -1.0305723052840483) H with fermi= 1.0 adding state (0, -0.1414106130994221) H with fermi= 1.0 Total density has weight 20.000000000000004 Itteration 20 Etot[Ry]= -1351.484605658407 Etot[Hartre]= -675.7423028292035 Diff= 2.4727287382120267e-05 Found bound state at E=-143.935187908H l= 0 Found bound state at E= -15.046905399H l= 0 Found bound state at E= -1.706330465H l= 0 Found bound state at E= -0.141410395H l= 0 Found bound state at E= -12.285376071H l= 1 Found bound state at E= -1.030572178H l= 1 Found bound state at E= -0.053407435H l= 1 Found bound state at E= -0.083077326H l= 2 Found bound state at E= 0.022340511H l= 2 adding state (0, -143.93518790834992) H with fermi= 1.0 adding state (0, -15.046905399002416) H with fermi= 1.0 adding state (1, -12.285376071157) H with fermi= 1.0 adding state (0, -1.7063304650893316) H with fermi= 1.0 adding state (1, -1.03057217822073) H with fermi= 1.0 adding state (0, -0.14141039494462587) H with fermi= 1.0 Total density has weight 20.0 Itteration 21 Etot[Ry]= -1351.4845940562832 Etot[Hartre]= -675.7422970281416 Diff= 1.160212377726566e-05 Found bound state at E=-143.935185231H l= 0 Found bound state at E= -15.046905393H l= 0 Found bound state at E= -1.706330394H l= 0 Found bound state at E= -0.141410279H l= 0 Found bound state at E= -12.285376187H l= 1 Found bound state at E= -1.030572114H l= 1 Found bound state at E= -0.053407222H l= 1 Found bound state at E= -0.083077222H l= 2 Found bound state at E= 0.022343632H l= 2 adding state (0, -143.93518523119468) H with fermi= 1.0 adding state (0, -15.046905392504758) H with fermi= 1.0 adding state (1, -12.285376186926461) H with fermi= 1.0 adding state (0, -1.7063303941503092) H with fermi= 1.0 adding state (1, -1.030572113718685) H with fermi= 1.0 adding state (0, -0.1414102789922587) H with fermi= 1.0 Total density has weight 20.000000000000004 Itteration 22 Etot[Ry]= -1351.4845884625638 Etot[Hartre]= -675.7422942312819 Diff= 5.593719379248796e-06 Found bound state at E=-143.935183948H l= 0 Found bound state at E= -15.046905388H l= 0 Found bound state at E= -1.706330359H l= 0 Found bound state at E= -0.141410218H l= 0 Found bound state at E= -12.285376243H l= 1 Found bound state at E= -1.030572081H l= 1 Found bound state at E= -0.053407108H l= 1 Found bound state at E= -0.083077168H l= 2 Found bound state at E= 0.022345249H l= 2 adding state (0, -143.9351839476714) H with fermi= 1.0 adding state (0, -15.046905387534638) H with fermi= 1.0 adding state (1, -12.285376242836591) H with fermi= 1.0 adding state (0, -1.706330358682587) H with fermi= 1.0 adding state (1, -1.0305720814813182) H with fermi= 1.0 adding state (0, -0.14141021780661645) H with fermi= 1.0 Total density has weight 19.999999999999996 Itteration 23 Etot[Ry]= -1351.4845858674255 Etot[Hartre]= -675.7422929337127 Diff= 2.5951383122446714e-06 Found bound state at E=-143.935183327H l= 0 Found bound state at E= -15.046905380H l= 0 Found bound state at E= -1.706330337H l= 0 Found bound state at E= -0.141410182H l= 0 Found bound state at E= -12.285376265H l= 1 Found bound state at E= -1.030572061H l= 1 Found bound state at E= -0.053407044H l= 1 Found bound state at E= -0.083077136H l= 2 Found bound state at E= 0.022346088H l= 2 adding state (0, -143.9351833271321) H with fermi= 1.0 adding state (0, -15.046905379991674) H with fermi= 1.0 adding state (1, -12.2853762654924) H with fermi= 1.0 adding state (0, -1.706330336794817) H with fermi= 1.0 adding state (1, -1.0305720611904292) H with fermi= 1.0 adding state (0, -0.14141018169658603) H with fermi= 1.0 Total density has weight 20.000000000000004 Itteration 24 Etot[Ry]= -1351.4845843857597 Etot[Hartre]= -675.7422921928799 Diff= 1.481665776736918e-06 Found bound state at E=-143.935183032H l= 0 Found bound state at E= -15.046905379H l= 0 Found bound state at E= -1.706330329H l= 0 Found bound state at E= -0.141410166H l= 0 Found bound state at E= -12.285376280H l= 1 Found bound state at E= -1.030572055H l= 1 Found bound state at E= -0.053407014H l= 1 Found bound state at E= -0.083077123H l= 2 Found bound state at E= 0.022346520H l= 2 adding state (0, -143.93518303227134) H with fermi= 1.0 adding state (0, -15.046905379367823) H with fermi= 1.0 adding state (1, -12.285376279807599) H with fermi= 1.0 adding state (0, -1.7063303294180896) H with fermi= 1.0 adding state (1, -1.030572054585891) H with fermi= 1.0 adding state (0, -0.14141016621793162) H with fermi= 1.0 Total density has weight 20.0 Itteration 25 Etot[Ry]= -1351.4845838385568 Etot[Hartre]= -675.7422919192784 Diff= 5.472029442898929e-07 Found bound state at E=-143.935182889H l= 0 Found bound state at E= -15.046905378H l= 0 Found bound state at E= -1.706330325H l= 0 Found bound state at E= -0.141410157H l= 0 Found bound state at E= -12.285376286H l= 1 Found bound state at E= -1.030572050H l= 1 Found bound state at E= -0.053406996H l= 1 Found bound state at E= -0.083077115H l= 2 Found bound state at E= 0.022346744H l= 2 adding state (0, -143.9351828892915) H with fermi= 1.0 adding state (0, -15.046905377799803) H with fermi= 1.0 adding state (1, -12.285376285604356) H with fermi= 1.0 adding state (0, -1.7063303247022823) H with fermi= 1.0 adding state (1, -1.0305720502393798) H with fermi= 1.0 adding state (0, -0.14141015675062057) H with fermi= 1.0 Total density has weight 19.999999999999996 Itteration 26 Etot[Ry]= -1351.4845836045583 Etot[Hartre]= -675.7422918022792 Diff= 2.3399843485094607e-07 Found bound state at E=-143.935182817H l= 0 Found bound state at E= -15.046905374H l= 0 Found bound state at E= -1.706330319H l= 0 Found bound state at E= -0.141410149H l= 0 Found bound state at E= -12.285376285H l= 1 Found bound state at E= -1.030572045H l= 1 Found bound state at E= -0.053406985H l= 1 Found bound state at E= -0.083077108H l= 2 Found bound state at E= 0.022346861H l= 2 adding state (0, -143.93518281704658) H with fermi= 1.0 adding state (0, -15.046905373595223) H with fermi= 1.0 adding state (1, -12.285376285021568) H with fermi= 1.0 adding state (0, -1.7063303191524133) H with fermi= 1.0 adding state (1, -1.0305720448874949) H with fermi= 1.0 adding state (0, -0.14141014937882082) H with fermi= 1.0 Total density has weight 20.0 Itteration 27 Etot[Ry]= -1351.484583297219 Etot[Hartre]= -675.7422916486095 Diff= 3.0733940548088867e-07 Found bound state at E=-143.935182785H l= 0 Found bound state at E= -15.046905375H l= 0 Found bound state at E= -1.706330319H l= 0 Found bound state at E= -0.141410148H l= 0 Found bound state at E= -12.285376288H l= 1 Found bound state at E= -1.030572045H l= 1 Found bound state at E= -0.053406981H l= 1 Found bound state at E= -0.083077107H l= 2 Found bound state at E= 0.022346921H l= 2 adding state (0, -143.93518278543843) H with fermi= 1.0 adding state (0, -15.046905374652523) H with fermi= 1.0 adding state (1, -12.285376287868589) H with fermi= 1.0 adding state (0, -1.706330319433535) H with fermi= 1.0 adding state (1, -1.0305720452334444) H with fermi= 1.0 adding state (0, -0.14141014782359984) H with fermi= 1.0 Total density has weight 20.0 Itteration 28 Etot[Ry]= -1351.4845833280006 Etot[Hartre]= -675.7422916640003 Diff= 3.078162080782931e-08
Homework¶
Compute charge density and total energy for oxygen atom, as we have done in class. Check that your code agrees with NIST values.
Get the code working for potassium (Z=19) and check that LDA respects aufmau principle of the periodic table.
How far (in periodic table) can you push the algorithm? What part needs improvements to extend the algorithm to higher Z? Currenly Z=21 is breaking down due to some instability.