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 ”hidden” into approximate 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})=E \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.
First we take the code from the Hydrogen project and repeat.
from scipy import *
from scipy import integrate
from scipy import optimize
from numpy import *
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
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-8,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.999922109 Found bound state at E= -0.249990190 Found bound state at E= -0.111108201 Found bound state at E= -0.062498772 Found bound state at E= -0.039999314 Found bound state at E= -0.250000016 Found bound state at E= -0.111111117 Found bound state at E= -0.062500003 Found bound state at E= -0.039999959 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.9999221089559623) with fermi= 1.0 adding state (0, -0.24999019020653063) with fermi= 1.0 adding state (1, -0.25000001561170354) with fermi= 1.0 adding state (0, -0.11110820082299863) with fermi= 1.0 adding state (1, -0.11111111678092336) with fermi= 1.0 adding state (2, -0.11111111114690274) with fermi= 1.0
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). Mathematically, this term is \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$).
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} In Hartree approximation, we have \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)
plot(R, ones(len(R))*2*Z)
show()
Poisson equation does not have the first order derivative, hence it can also be more efficiently solved by the Numerov algorithm.
Recall the Numerov algorithm as appropriate for the Poisson equation
\begin{eqnarray} x(h)+x(-h) = 2x(0)+h^2 (f(0)x(0)+u(0))+\frac{2}{4!}h^4 x^{(4)}(0)+O(h^6) \end{eqnarray}\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, dt):
x = zeros(len(U))
x[0] = x0
x[1] = dx*dt + x0
h2 = dt*dt
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
ux = -8*pi*R*rho
U2 = NumerovUP(ux, 0.0, 0.5, R[1]-R[0])
alpha2 = (2*Z-U2[-1])/R[-1]
U2 += alpha2 * R
plot(R,U0)
plot(R,U2);
For generic density the following routine will work:
def HartreeU(R, rho, Zatom):
ux = -8*pi*R*rho
U2 = NumerovUP(ux, 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);
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.)
Vxc = [2*exc.Vc(rs(rh)) + 2*exc.Vx(rs(rh)) for rh in rho]
Uks = U2 - 2*Z + Vxc*R
plot(R, -Uks/2);
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)=0. \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)=0. \end{equation}
@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]
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
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('Found bound state at E=%14.9f' % Ebound)
u0 = u1
return Ebnd
#R = linspace(1e-8,10,2**13+1)
#Zatom = 8
#E0=-1.2*Zatom**2
#Eshift=0.5 # sometimes energies can be positive!!!
#Esearch = -logspace(-4,log10(-E0+Eshift),200)[::-1] + Eshift
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.555322971 Found bound state at E= -0.206605857 Found bound state at E= -0.059520239 Found bound state at E= -0.004720150 Found bound state at E= -0.640986418 Found bound state at E= -0.243652142 Found bound state at E= -0.071991382 Found bound state at E= -0.007466669 Found bound state at E= -0.600054533 Found bound state at E= -0.217328796 Found bound state at E= -0.057702239 Found bound state at E= -0.558268879 Found bound state at E= -0.187465773
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.
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), 'with fermi=', ferm)
if N>=Zatom: break
return (rho,Ebs)
rho_new, Ebs = ChargeDensity(Bnd,R,Z,Uks)
adding state (1, -0.3204932088844849) with fermi= 1.0 adding state (2, -0.30002726657005885) with fermi= 1.0 adding state (0, -0.27766148533971746) with fermi= 1.0 adding state (3, -0.27913443951024725) with fermi= 0.7142857142857143
For oxygen, the total energy in this implementation is : -74.47303426133809 Hartree, while NIST shows -74.473077 Hartree. This is in excellent agreement. The small loss of accuracy is mainly due to the mesh used for R.