Self-consistency and the density

Department of Applied Physics, Chalmers, Göteborg
Göran Wahnström, Per Sundell, Martin Petisme, 2014-01-24
Self-consistency and the density-functional theory
In this exercise you are asked to determine the ground state energy for the helium atom using the Hartree
and the density functional method. You are recommended to use Python or .
e Helium Atom Hamiltonian
e Born-Oppenheimer Hamiltonian for the many-electron system reads
H=−
1 1 ∑ e2
ℏ2 ∑ 2
1 ∑ ∑ Zn e 2
∇i −
+
2me i
4πε 0 n i |rri − R n | 2 4πε 0
|rri − r j |
(1)
1 ∑ 2 ∑∑
1∑
Zn
1
∇i −
+
r
r
|r
−
R
|
|r
−
rj|
2 i
2
i
n
i
n
i
(2)
i̸=j
or in atomic units
H=−
i̸=j
For the helium atom it simpliﬁes to
1
2
2
1
H = − (∇21 + ∇22 ) − − +
r1
r2
r12
2
(3)
and where r 1 and r 2 denote the position coordinates of the two electrons and r12 = |rr1 − r 2 |.
Hartree–Fo for helium
e two-electron Hartree–Fock Ansatz, which satisﬁes the antisymmetry requirement that the wavefunction changes sign when two particles are exchanged, is
1
Ψ(xx1 , x 2 ) = √ [ψ 1 (xx1 )ψ 2 (xx2 ) − ψ 1 (xx2 )ψ 2 (xx1 )]
2
(4)
where the notation x i includes both the spatial and spin coordinate. For the ground state with paired
electrons we have ψ 1 (xx1 ) = φ 0 (rr1 )α(s1 ) and ψ 2 (xx2 ) = φ 0 (rr2 )β(s2 ) so that
1
Ψ(xx1 , x 2 ) = √ φ 0 (rr1 )φ 0 (rr2 )[α(s1 )β(s2 ) − α(s2 )β(s1 )]
2
(5)
where φ 0 is normalized and for the spin components we have |α(si )α(si )| = |β(si )β(si )| = 1 and α(si )β(si ) =
0.
Plugging our Ansatz into the Schrödinger equation and integrating our the r 2 dependence (see [1]),
we get the equation
[
]
∫
1 2
2
2 1
− ∇1 − + drr2 |φ(rr2 )|
φ(rr1 ) = E′ φ(rr1 ).
(6)
r1
r12
2
Density Functional eory
e density functional method is based on the observation by Hohenberg and Kohn (1964) that all the
ground-state properties of a many-body quantum-mechanical system of electrons may be obtained from
a knowledge of the electron density n(rr). ey proved that the total energy of a many-electron system
in an external potential is a unique functional of the electron density
∫
E[n] = F[n] + drr Vext (rr)n(rr)
(7)
and this functional is minimum and equal to the ground-state energy E0 when n(rr) is the ground-state
density n0 (rr),
E0 = E[n0 (rr)] = min E[n(rr)]
(8)
n(rr)
e notation
minn(rr) implies minimizing with respect to the density n(rr) keeping the number of electrons
∫
N = drr n(rr) constant.
1
e Kohn–Sham equations
A practical scheme was developed by Kohn and Sham (1965). ey wrote the functional on the form
∫
KS
E [n] = T0 [n] + EH [n] + Exc [n] + drr Vext (rr)n(rr)
(9)
where T0 [n] is the kinetic energy of a (ﬁctious) non-interacting electron system; EH [n] is the Hartree energy, the classical inter-electron repulsion; and Exc [n] is the exchange-correlation energy, the rest. e
laer includes all complicated many-body eﬀects and has to be approximated. Carrying out the minimization in equation (8) one arrives at the Kohn–Sham equations
[
]
1
− ∇2 + Veﬀ (rr) ψ i (rr) = ε i ψ i (rr)
(10)
2
with
Veﬀ (rr) = VH (rr) + Vxc (rr) + Vext (rr)
and
∫
VH (rr) =
Vxc (rr)
=
Vext (rr) =
with
n(rr) =
(11)
n(rr′ )
,
|rr − r′ |
δExc [n]
,
δn(rr)
∑ Zn
−
,
|rr − R n |
n
N
∑
drr′
(12)
(13)
(14)
|ψ i (rr)|2 .
(15)
i=1
e ground-state energy is then given by
E0 =
N
∑
i=1
εi −
1
2
∫
drr drr′
n(rr)n(rr′ )
+ Exc [n] −
|rr − r ′ |
∫
drr Vxc (rr)n(rr).
(16)
Notice that the ε i ’s are Lagrange multipliers introduced to handle the constraint that the number of electrons should be conserved,
∫
N=
drr n(rr) = const.
2
(17)
For the helium atom the Kohn–Sham equation reduces to
[
]
∫
1
2
|φ(rr′ )|2
r
− ∇2 − + 2 drr′
+
V
(r
)
φ(rr) = εφ(rr)
xc
r
|rr − r ′ |
2
with
1
E0 = 2ε −
2
∫
n(rr)n(rr′ )
drr drr
+ Exc [n] −
|rr − r′ |
′
and
(18)
∫
drr Vxc (rr)n(rr)
n(rr) = 2|φ(rr)|2 .
(19)
(20)
e Local Density Approximation
e exact form of the exchange-correlation functional is not known. It has to be approxiamted. e most
widely used is the Local Density Approximation (). In that approximation the exchange-correlation
functional Exc [n] is represented as
∫
Exc [n] = drr n(rr)ε hom
(21)
xc (n(rr ))
where ε hom
xc (n) is the exchange-correlation energy per electron of an homogeneous electron gas at density
n. e corresponding potential is given by
Vxc (rr) ≡
δExc [n]
d hom
= ε hom
ε (n)
xc (n) + n
δn
dn xc
(22)
e exchange-correlation energy can be separated into an exchange and a correlation contribution
hom
hom
ε hom
xc (n) = ε x (n) + ε c (n)
(23)
and correspondingly for the exchange-correlation potential. e exchange part can be evaluated exactly[2]
ε hom
x (n)
=−
3
4
(
3n
π
)1/3
(24)
and the correlation part has been evaluated essentially exactly by the numerical quantum Monte Carlo
technique.[3] e following analytical respresentation has been suggested[4]
ε hom
c (rs ) =
γ
, rs ≥ 1
√
1 + β 1 rs + β 2 rs
(25)
ε hom
c (rs ) =
A ln rs + B + Crs ln rs + Drs , rs < 1
(26)
with n = 3/(4πr3s ). e numerical values for the non-spinpolarized electron gas (in atomic units) are
A = 0.0311, B = −0.048, C = 0.0020, D = −0.0116, γ = −0.1423, β 1 = 1.0529, β 2 = 0.3334.
(27)
Numerical procedure
e Kohn–Sham equation (18) is non-linear in φ(rr) and has to be solved in an iterative manner. An initial
guess for the wave-function φ(rr) is made, the Hartree potential VH (r) and exchange-correlation potential
Vxc (r) are determined and the Kohn–Sham equation (18) can then be solved. A new wave-function φ(rr)
is obtained and the procedure is repeated until the ground-state energy E0 in equation (19) is converged.
3
e Hartree potential
e Hartree potential, the electrostatic potential generated by the charge distribution,
∫
n(rr′ )
VH (rr) ≡ drr′
|rr − r ′ |
(28)
which also can be wrien on diﬀerential form, the Poisson’s equation
∇2 VH (rr) = −4πn(rr)
(29)
For the helium atom the ground-state electron density is given by
n(rr) ≡ 2ns (rr) = 2|φ(rr)|2
(30)
where ns (rr) is the density of a single orbital. We can introduce the corresponding electro-static potential
according to
∇2 VsH (rr) = −4πns (rr)
(31)
Using the radial symmetry it takes the form
d2
U(r) = −4πrns (r)
dr2
with U(r) = rVsH (r). e density is normalized according to
∫
∫
drr ns (rr) = 4π dr r2 ns (r) = 1
and by introducing the variable
u(r) =
√
4πns (r) r =
√
4π rφ(r)
(32)
(33)
(34)
the following diﬀerential equation is obtained
d2
u2 (r)
U(r)
=
−
dr2
r
with the boundary conditions (see ijssen [1] for details)
U(0) = 0,
(35)
U(rmax ) = qmax = 1
By introducing the variable
U0 (r) ≡ U(r) −
r
rmax
the two-point boundary value problem is transformed into the standard form
d2
u2 (r)
U0 (r) = −
2
dr
r
(36)
(37)
with
U0 (0) = U0 (rmax ) = 0.
e radial Srödinger equation
In radial coordinates the Kohn–Sham equation takes the form
[
]
1 d2
2
−
− + VH (r) + Vx (r) + Vc (r) u(r) = εu(r)
r
2 dr2
(38)
with the boundary conditions
u(0) = u(rmax ) = 0
Finally, the ground-state energy, in the local-density approximation, is given by
]
[
∫
1
E0 = 2ε − 2 dr u2 (r) VH (r) + Vxc (r) − ε xc (r)
2
4
(39)
Boundary value problems
e boundary value problem for the diﬀerential equation (37) and eigenvalue equation (38) can be solved
with various numerical methods. e ﬁnite diﬀerence method converts the boundary value problem into
a system of algebraic equations. is is done by replacing the derivatives with ﬁnite diﬀerences.
Finite diﬀerences
In the ﬁnite diﬀerence method the space is discretized
xi = a + ih,
i = 0, . . . , n,
h = (b − a)/n,
(40)
and one seeks an approximate solution at the grid points
yi = y(xi ).
(41)
e derivatives are replaced by ﬁnite diﬀerences
y′ (xi ) =
y′′ (xi ) =
yi+1 − yi−1
2h
yi+1 − 2yi + yi−1
h2
and the resulting matrix equation can be solved using standard routines.
5
(42)
(43)
Task
1. Solve the problem in Section 4.3.2 in ijssen [1], i.e. write the wave function as
φ(rr) =
4
∑
Ci χ i (rr)
i=1
with
χ i (rr) = exp(−α i r2 )
and
α 1 = 0.297104,
α 2 = 1.236745,
α 3 = 5.749982,
α 4 = 38.216677.
en solve the equations resulting from plugging this into equation (6) self-consistently until your
energy change is less than 10−5 eV. Determine the ground state energy and wave function.
2. Consider Poisson’s equation. Solve it in radial coordinates (equation (37)) using ﬁnite diﬀerences.
Discretize the equation and approximate the derivates using a uniform grid. en, solve the resulting system of linear equations.
Test your program by using the ground state electron density for the hydrogen atom as input density.
You should obtain the Hartree potential
(
)
1
1
VH (r) = − 1 +
exp(−2r).
r
r
3. You also need to solve the Kohn–Sham equation (38). Write a program that solves the eigenvalue
problem using the ﬁnite diﬀerence method on a uniform grid.
Test your program by solving the radial Schrödinger equation for the hydrogen atom
(
)
1 d2
1
−
−
u(r) = Eu(r).
r
2 dr2
Compare your results with the analytical result for the ground state energy and wave function.
4. Use the above results to construct a program that determines the ground state energy in the Hartree
approximation. In this case, leave out the exchange-correlation terms, and remove the self-interaction
from the Hartree potential, i.e. use VH = VsH Make an initial guess for the density and iterate to
self-consistency. Stop the loop when the energy changes with less than 10−5 eV. Determine the
ground state energy, eigenvalue, and wave function. Make sure that your solution is converged
with respect to both the grid spacing and rmax . (Suggestion: increase rmax keeping the same grid
density until convergence, and then increase the grid density until convergence.) Compare your
results with Task 1.
5. Add exchange contributions to the previous task and solve it. (In this case, you should use the full
Hartree potential, i.e. VH = 2VsH .)
6. Finally, solve the task also including the correlation contributions.
References
[1] J.M. ijssen Computational Physics, Cambridge, (1999)
[2] N.W. Ashcro and N.D. Mermin Solid State Physics, Holt, Rinehart and Winston (1976)
[3] D.M. Ceperley and B.J. Alder, Ground State of the Electron Gas by a Stochastic Method,
Phys. Rev. Le. 45, 566 (1980)
[4] J.P. Perdew, A. Zunger, Self-interaction correction to density-functional approximations for manyelectron systems, Phys. Rev. B 23, 5048 (1981)
6

Download Report