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.simpson(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) 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 $V_{H}(r) = U(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') xlim([0,20]) legend(loc='best') show()
plot(R, Z-U0/2,label='effective charge') xlim([1,20]) 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,'.-',label='odeint') plot(R,U2,'-',label='Numerov') xlim([0,20]) legend(loc='best') show()
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) show()
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+2*Z, label='-Uh+2Z') #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() ax1.set_xlim([0,20]) ax2.set_xlim([0,20]) show()
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.simpson(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) #poly = Polynomial.fit(R[:2], ur[:2], deg=1) 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=22 E0=-1.0*Zatom**2 Eshift=0.1 # sometimes energies can be positive!!! Esearch = -logspace(-4,log10(-E0+Eshift),200)[::-1] + Eshift #Esearch = -logspace(-4,log10(-E0+Eshift),500)[::-1] + Eshift plot(Esearch,'o-') ylim([-1.2*Zatom**2,1]) show() #ylim([-0.2,0.])
We will first demonstrate algorithm on oxygen $Z=8$ and than on Potassium ($Z=19$), which has first nontrivial "aufbau" principle in which 4ドルs$ orbital gets filled before 3ドルd$.
R = linspace(1e-8,10,2**14+1) Zatom = 8 # 30 #29 works 20 # 19 mixr = 0.5 E0=-1.2*Zatom**2 Eshift=0.1 # sometimes energies can be positive!!! Esearch = -logspace(-4,log10(-E0+Eshift),500)[::-1] + Eshift nmax = 4 #3 exc = ExchangeCorrelation() Uks = -2*Zatom*ones(len(R)) # here I changed from -2 to -2*Zatom Eold = 0 Etol = 1e-7 for itt in range(100): 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.simpson(pot, x=R) Etot = epot + Ebs print('Total density has weight', integrate.simpson(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()
Found bound state at E= -32.000000000H l= 0 Found bound state at E= -8.000000000H l= 0 Found bound state at E= -3.555555556H l= 0 Found bound state at E= -1.999999922H l= 0 Found bound state at E= -8.000000000H l= 1 Found bound state at E= -3.555555556H l= 1 Found bound state at E= -1.999999948H l= 1 Found bound state at E= -3.555555556H l= 2 Found bound state at E= -1.999999978H l= 2 adding state (0, -32.000000000317144) H with fermi= 1.0 adding state (0, -8.000000000006459) H with fermi= 1.0 adding state (1, -7.999999999996606) H with fermi= 0.6666666666666666 Total density has weight 8.0 Itteration 0 Etot[Ry]= -330.37370870782803 Etot[Hartre]= -165.18685435391401 Diff= 330.37370870782803 Found bound state at E= -14.313207892H l= 0 Found bound state at E= -0.018196371H l= 0 adding state (0, -14.313207892090661) H with fermi= 1.0 adding state (0, -0.01819637127074424) H with fermi= 1.0 Total density has weight 5.999999999999999 Itteration 1 Etot[Ry]= -112.66205148843432 Etot[Hartre]= -56.33102574421716 Diff= 217.71165721939371 Found bound state at E= -19.162955385H l= 0 Found bound state at E= -1.850079828H l= 0 Found bound state at E= -0.396156915H l= 0 Found bound state at E= -0.049235428H l= 0 Found bound state at E= -1.288042563H l= 1 Found bound state at E= -0.278105803H l= 1 Found bound state at E= 0.005596814H l= 1 Found bound state at E= -0.172494716H l= 2 adding state (0, -19.16295538482916) H with fermi= 1.0 adding state (0, -1.8500798283089868) H with fermi= 1.0 adding state (1, -1.2880425629058951) H with fermi= 0.6666666666666666 Total density has weight 7.0 Itteration 2 Etot[Ry]= -157.2415904747316 Etot[Hartre]= -78.6207952373658 Diff= 44.57953898629728 Found bound state at E= -18.754350439H l= 0 Found bound state at E= -1.264387318H l= 0 Found bound state at E= -0.150146103H l= 0 Found bound state at E= -0.713173010H l= 1 Found bound state at E= -0.066038436H l= 1 Found bound state at E= 0.017327547H l= 2 adding state (0, -18.75435043859003) H with fermi= 1.0 adding state (0, -1.2643873184645649) H with fermi= 1.0 adding state (1, -0.7131730095705805) H with fermi= 0.6666666666666666 Total density has weight 7.5 Itteration 3 Etot[Ry]= -151.30480523623677 Etot[Hartre]= -75.65240261811839 Diff= 5.936785238494821 Found bound state at E= -18.735348427H l= 0 Found bound state at E= -1.044169033H l= 0 Found bound state at E= -0.059543197H l= 0 Found bound state at E= -0.502432661H l= 1 Found bound state at E= 0.007681966H l= 1 adding state (0, -18.735348426756776) H with fermi= 1.0 adding state (0, -1.0441690325248625) H with fermi= 1.0 adding state (1, -0.5024326607512009) H with fermi= 0.6666666666666666 Total density has weight 7.75 Itteration 4 Etot[Ry]= -149.77979609137316 Etot[Hartre]= -74.88989804568658 Diff= 1.525009144863617 Found bound state at E= -18.751895759H l= 0 Found bound state at E= -0.951526061H l= 0 Found bound state at E= -0.020672905H l= 0 Found bound state at E= -0.414448169H l= 1 Found bound state at E= 0.039778296H l= 1 adding state (0, -18.751895759307637) H with fermi= 1.0 adding state (0, -0.9515260606492291) H with fermi= 1.0 adding state (1, -0.41444816886492225) H with fermi= 0.6666666666666666 Total density has weight 7.875000000000001 Itteration 5 Etot[Ry]= -149.32265739553657 Etot[Hartre]= -74.66132869776828 Diff= 0.4571386958365906 Found bound state at E= -18.758178791H l= 0 Found bound state at E= -0.908711810H l= 0 Found bound state at E= -0.002441622H l= 0 Found bound state at E= -0.373839481H l= 1 adding state (0, -18.758178791143) H with fermi= 1.0 adding state (0, -0.9087118101631386) H with fermi= 1.0 adding state (1, -0.37383948062795047) H with fermi= 0.6666666666666666 Total density has weight 7.9375 Itteration 6 Etot[Ry]= -149.12252511525537 Etot[Hartre]= -74.56126255762769 Diff= 0.20013228028119556 Found bound state at E= -18.759659685H l= 0 Found bound state at E= -0.888769863H l= 0 Found bound state at E= 0.006342592H l= 0 Found bound state at E= -0.354922279H l= 1 adding state (0, -18.759659684930984) H with fermi= 1.0 adding state (0, -0.8887698631062732) H with fermi= 1.0 adding state (1, -0.3549222790517335) H with fermi= 0.6666666666666666 Total density has weight 7.968750000000001 Itteration 7 Etot[Ry]= -149.02941391088433 Etot[Hartre]= -74.51470695544216 Diff= 0.09311120437104137 Found bound state at E= -18.759595553H l= 0 Found bound state at E= -0.879476669H l= 0 Found bound state at E= 0.010599573H l= 0 Found bound state at E= -0.346100228H l= 1 adding state (0, -18.75959555344527) H with fermi= 1.0 adding state (0, -0.8794766686520796) H with fermi= 1.0 adding state (1, -0.3461002278509865) H with fermi= 0.6666666666666666 Total density has weight 7.984375 Itteration 8 Etot[Ry]= -148.98561034622242 Etot[Hartre]= -74.49280517311121 Diff= 0.04380356466191415 Found bound state at E= -18.759193488H l= 0 Found bound state at E= -0.875141677H l= 0 Found bound state at E= 0.012656318H l= 0 Found bound state at E= -0.341980732H l= 1 adding state (0, -18.759193488179804) H with fermi= 1.0 adding state (0, -0.8751416771908085) H with fermi= 1.0 adding state (1, -0.3419807322528575) H with fermi= 0.6666666666666666 Total density has weight 7.992187499999999 Itteration 9 Etot[Ry]= -148.96484989322687 Etot[Hartre]= -74.48242494661343 Diff= 0.020760452995546075 Found bound state at E= -18.758830598H l= 0 Found bound state at E= -0.873119049H l= 0 Found bound state at E= 0.013642362H l= 0 Found bound state at E= -0.340056328H l= 1 adding state (0, -18.75883059792217) H with fermi= 1.0 adding state (0, -0.873119049257727) H with fermi= 1.0 adding state (1, -0.3400563282194268) H with fermi= 0.6666666666666666 Total density has weight 7.996093749999999 Itteration 10 Etot[Ry]= -148.9549960008792 Etot[Hartre]= -74.4774980004396 Diff= 0.009853892347678084 Found bound state at E= -18.758580120H l= 0 Found bound state at E= -0.872175469H l= 0 Found bound state at E= 0.014110267H l= 0 Found bound state at E= -0.339157384H l= 1 adding state (0, -18.758580120030583) H with fermi= 1.0 adding state (0, -0.8721754689833804) H with fermi= 1.0 adding state (1, -0.33915738416210267) H with fermi= 0.6666666666666666 Total density has weight 7.998046875 Itteration 11 Etot[Ry]= -148.9503233116414 Etot[Hartre]= -74.4751616558207 Diff= 0.004672689237793293 Found bound state at E= -18.758425980H l= 0 Found bound state at E= -0.871735485H l= 0 Found bound state at E= 0.014329859H l= 0 Found bound state at E= -0.338737620H l= 1 adding state (0, -18.758425980451836) H with fermi= 1.0 adding state (0, -0.8717354854980826) H with fermi= 1.0 adding state (1, -0.33873761986087486) H with fermi= 0.6666666666666666 Total density has weight 7.999023437499999 Itteration 12 Etot[Ry]= -148.94811253490417 Etot[Hartre]= -74.47405626745208 Diff= 0.0022107767372290255 Found bound state at E= -18.758337005H l= 0 Found bound state at E= -0.871530478H l= 0 Found bound state at E= 0.014431841H l= 0 Found bound state at E= -0.338541737H l= 1 adding state (0, -18.758337005010496) H with fermi= 1.0 adding state (0, -0.8715304775010679) H with fermi= 1.0 adding state (1, -0.3385417365720025) H with fermi= 0.6666666666666666 Total density has weight 7.99951171875 Itteration 13 Etot[Ry]= -148.94706964465044 Etot[Hartre]= -74.47353482232522 Diff= 0.0010428902537285012 Found bound state at E= -18.758287697H l= 0 Found bound state at E= -0.871435050H l= 0 Found bound state at E= 0.014478762H l= 0 Found bound state at E= -0.338450408H l= 1 adding state (0, -18.75828769696511) H with fermi= 1.0 adding state (0, -0.8714350495188786) H with fermi= 1.0 adding state (1, -0.33845040813826976) H with fermi= 0.6666666666666666 Total density has weight 7.999755859375 Itteration 14 Etot[Ry]= -148.94657933222322 Etot[Hartre]= -74.47328966611161 Diff= 0.0004903124272175319 Found bound state at E= -18.758261134H l= 0 Found bound state at E= -0.871390685H l= 0 Found bound state at E= 0.014500165H l= 0 Found bound state at E= -0.338407876H l= 1 adding state (0, -18.75826113419569) H with fermi= 1.0 adding state (0, -0.8713906850538338) H with fermi= 1.0 adding state (1, -0.33840787558163754) H with fermi= 0.6666666666666666 Total density has weight 7.999877929687502 Itteration 15 Etot[Ry]= -148.94634966594393 Etot[Hartre]= -74.47317483297196 Diff= 0.00022966627929577044 Found bound state at E= -18.758247121H l= 0 Found bound state at E= -0.871370093H l= 0 Found bound state at E= 0.014509844H l= 0
Found bound state at E= -0.338388097H l= 1 adding state (0, -18.758247120902563) H with fermi= 1.0 adding state (0, -0.871370093218127) H with fermi= 1.0 adding state (1, -0.3383880970006329) H with fermi= 0.6666666666666666 Total density has weight 7.99993896484375 Itteration 16 Etot[Ry]= -148.9462425486771 Etot[Hartre]= -74.47312127433855 Diff= 0.00010711726682188782 Found bound state at E= -18.758239847H l= 0 Found bound state at E= -0.871360554H l= 0 Found bound state at E= 0.014514179H l= 0 Found bound state at E= -0.338378916H l= 1 adding state (0, -18.758239846615176) H with fermi= 1.0 adding state (0, -0.8713605539634439) H with fermi= 1.0 adding state (1, -0.3383789158812842) H with fermi= 0.6666666666666666 Total density has weight 7.999969482421876 Itteration 17 Etot[Ry]= -148.946192826413 Etot[Hartre]= -74.4730964132065 Diff= 4.9722264094498314e-05 Found bound state at E= -18.758236119H l= 0 Found bound state at E= -0.871356145H l= 0 Found bound state at E= 0.014516100H l= 0 Found bound state at E= -0.338374663H l= 1 adding state (0, -18.758236118977212) H with fermi= 1.0 adding state (0, -0.8713561452115111) H with fermi= 1.0 adding state (1, -0.33837466327112675) H with fermi= 0.6666666666666666 Total density has weight 7.9999847412109375 Itteration 18 Etot[Ry]= -148.9461698660551 Etot[Hartre]= -74.47308493302755 Diff= 2.2960357910051243e-05 Found bound state at E= -18.758234230H l= 0 Found bound state at E= -0.871354114H l= 0 Found bound state at E= 0.014516938H l= 0 Found bound state at E= -0.338372699H l= 1 adding state (0, -18.758234229685463) H with fermi= 1.0 adding state (0, -0.8713541138000558) H with fermi= 1.0 adding state (1, -0.338372699074085) H with fermi= 0.6666666666666666 Total density has weight 7.99999237060547 Itteration 19 Etot[Ry]= -148.94615933875014 Etot[Hartre]= -74.47307966937507 Diff= 1.0527304965535222e-05 Found bound state at E= -18.758233281H l= 0 Found bound state at E= -0.871353181H l= 0 Found bound state at E= 0.014517297H l= 0 Found bound state at E= -0.338371795H l= 1 adding state (0, -18.758233280553725) H with fermi= 1.0 adding state (0, -0.8713531807964497) H with fermi= 1.0 adding state (1, -0.3383717945421041) H with fermi= 0.6666666666666666 Total density has weight 7.9999961853027335 Itteration 20 Etot[Ry]= -148.94615454594503 Etot[Hartre]= -74.47307727297252 Diff= 4.792805100350961e-06 Found bound state at E= -18.758232807H l= 0 Found bound state at E= -0.871352754H l= 0 Found bound state at E= 0.014517448H l= 0 Found bound state at E= -0.338371379H l= 1 adding state (0, -18.758232807338246) H with fermi= 1.0 adding state (0, -0.8713527537778318) H with fermi= 1.0 adding state (1, -0.338371379317186) H with fermi= 0.6666666666666666 Total density has weight 7.999998092651366 Itteration 21 Etot[Ry]= -148.94615237534052 Etot[Hartre]= -74.47307618767026 Diff= 2.170604517459651e-06 Found bound state at E= -18.758232574H l= 0 Found bound state at E= -0.871352560H l= 0 Found bound state at E= 0.014517509H l= 0 Found bound state at E= -0.338371190H l= 1 adding state (0, -18.75823257359043) H with fermi= 1.0 adding state (0, -0.871352559800071) H with fermi= 1.0 adding state (1, -0.33837119006426036) H with fermi= 0.6666666666666666 Total density has weight 7.999999046325683 Itteration 22 Etot[Ry]= -148.94615141284305 Etot[Hartre]= -74.47307570642153 Diff= 9.624974666166963e-07 Found bound state at E= -18.758232459H l= 0 Found bound state at E= -0.871352472H l= 0 Found bound state at E= 0.014517532H l= 0 Found bound state at E= -0.338371104H l= 1 adding state (0, -18.758232458546455) H with fermi= 1.0 adding state (0, -0.871352471986121) H with fermi= 1.0 adding state (1, -0.3383711040735523) H with fermi= 0.6666666666666666 Total density has weight 7.999999523162842 Itteration 23 Etot[Ry]= -148.94615099142393 Etot[Hartre]= -74.47307549571197 Diff= 4.214191164919612e-07 Found bound state at E= -18.758232402H l= 0 Found bound state at E= -0.871352432H l= 0 Found bound state at E= 0.014517540H l= 0 Found bound state at E= -0.338371065H l= 1 adding state (0, -18.75823240186256) H with fermi= 1.0 adding state (0, -0.8713524319928878) H with fermi= 1.0 adding state (1, -0.33837106472823075) H with fermi= 0.6666666666666666 Total density has weight 7.999999761581421 Itteration 24 Etot[Ry]= -148.9461507928682 Etot[Hartre]= -74.4730753964341 Diff= 1.985557389616588e-07 Found bound state at E= -18.758232375H l= 0 Found bound state at E= -0.871352415H l= 0 Found bound state at E= 0.014517542H l= 0 Found bound state at E= -0.338371048H l= 1 adding state (0, -18.75823237527206) H with fermi= 1.0 adding state (0, -0.871352415131519) H with fermi= 1.0 adding state (1, -0.33837104805581364) H with fermi= 0.6666666666666666 Total density has weight 7.999999880790711 Itteration 25 Etot[Ry]= -148.94615073141514 Etot[Hartre]= -74.47307536570757 Diff= 6.145305064819695e-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. Check Z=29 first and than check what you need to change to get Z=22-28 working.