New Page 2

Chapter 12 -THE SELF CONSISTENT FIELD METHOD

by Reinaldo Baretti Machín

e mail : reibaretti2004@yahoo.com

Lectures on Quantum Mechanics -Introduction

Chapter 1 - Some experimental facts

Chapter 2 - The Old Quantum Theory

Chapter 3 - The Postulates of Quantum Mechanics

Chapter 4- Potential Barriers

Chapter 5 - Potential Wells

Chapter 6 - Periodic Potentials

Chapter 7 - Angular Momentum

Chapter 8 - The Hydrogen Atom

Chapter 9 - The Electron Spin

Chapter 10 - The Two Electron Atom

Chapter 11- Ground State Energies of Lithium and Beryllium

Chapter 12- The Self Consistent Field Method

Chapter ** -The H₂ molecule

References :

1. The Hartree -Fock Method , http://en.wikipedia.org/wiki/Hartree-Fock

2. Principles of quantum mechanics; nonrelativistic wave mechanics with illustrative applications.

by W. V. Houston

3. Introduction to Quantum Mechanics with Applications to Chemistry by Linus Pauling and E. Bright Wilson Jr.

4. The Hartree-Fock Method for Atoms: A Numerical Approach by Charlotte Froese Fischer

5. Douglas Rayner Hartree - His Life in Science and Computing - Charlotte Froese Fischer

Charlotte Froese Fischer Douglas R. Hartree

Vladimir Aleksandrovich Fock.

Our objective is to solve the one electron equation in the case of helium ground state by numerical and analytical procedures,

showing by the example, the main difficulties of the self consistent method in atomic physics. The end result is called the Hartree-Fock energy of helium.

The one electron equation is ,

-(1/2)∆ φ_i⁽ⁿ⁾ -(Z/r) φ_i⁽ⁿ⁾ + [ ∫ ∑_j≠i(/r - r' / )^-1 {φ_i^(n-1) }² dτ'] φ_i⁽ⁿ⁾ = ε_i φ_i⁽ⁿ⁾ , (1)

-(1/2)∆ φ_i⁽ⁿ⁾ -(Z/r) φ_i⁽ⁿ⁾ + V(r) φ_i⁽ⁿ⁾ = ε_i φ_i⁽ⁿ⁾ , (2)

One has to solve (2) for N electrons. Actually the work can be halved halved since two electron with opposite spins share the same space wave function. For example in helium ground , the two electrons have the same φ_1s (r0 wave function.

in a many electron problem (neglecting exchange) .

In the case of helium

where V(r₁) = 1/r₁ - ( Z +1/r₁ ) exp(-2Zr₁) (3)

∆ V = - 4π ρ = - 4π φ ² (4)

In the case of spherical symmetry the differential equation reduces to

d²V/dr² + (2/r) dV/dr = - 4π ρ(r) . (5)

Approximating the derivatives as differences one gets

{V(r +∆r) - 2V(r) +V(r-∆r) }/∆r² + (2/r) {V(r) -V(r-∆r) }/∆r = - 4π ρ(r) . (6)

The finite difference solution of the equation is

V(r +∆r) = 2V(r) - V(r-∆r) -(2∆r) /r) {V(r) -V(r-∆r) }- (∆r)² 4π ρ(r) , (7)

where two initial conditions are needed.

Given that V(r) = ∫ ρ(r')/r - r' / )^-1 dτ' we have the initial condition

V(r=0) = ∫ ( ρ(r')/ r' ) dτ' . (10)

The initial condition on the derivative of V(r) at the origin is obtained from the DE (5). Multiplying (5) by r and letting r go to zero gives

2 (dV/dr)_r=0 = 0 , hence

(dV/dr)_r=0 = 0 . (11)

so V(0) is obtained from (10) and V(∆r) = V(0) because dV/dr)_r=0 = 0 .

As an example we find V(r) numerically for the case of helium where we assume φ_1s(r) = ( Z³ /π)^1/2 exp(-Zr) .The analytic answer is known to be (see ref 2. appendix V or Lecures on Quantum Mechanics, Chapter 10 - The Two Electron Atom )

V_analytic(r) = (1/r){1-exp(-2Zr) } - Z exp(-2Zr) . (12)

It is seen that lim_r→0 ( V_analytic(r) ) = Z .

var('r');Z=2;
v=(1/r)*(1-exp(-2*Z*r) ) - Z*exp(-2*Z*r);
y=plot(v,r,0.0001,5);
show(y)

Fig 1. Plot V(r) = (1/r){1-exp(-2Zr) } - Z exp(-2Zr) .

Results:

The table gives the numerical solution from eq (9) and the exact solution eq .(12).

r,phi(exact),phinum= 0.1000E-03 0.2000E+01 0.2000E+01
r,phi(exact),phinum= 0.1767E-02 0.2000E+01 0.2000E+01
r,phi(exact),phinum= 0.8510E-01 0.1967E+01 0.1967E+01
r,phi(exact),phinum= 0.1684E+00 0.1891E+01 0.1890E+01
r,phi(exact),phinum= 0.2518E+00 0.1790E+01 0.1789E+01
r,phi(exact),phinum= 0.3351E+00 0.1680E+01 0.1678E+01
r,phi(exact),phinum= 0.4184E+00 0.1567E+01 0.1565E+01
r,phi(exact),phinum= 0.5018E+00 0.1456E+01 0.1455E+01
r,phi(exact),phinum= 0.5851E+00 0.1352E+01 0.1350E+01
r,phi(exact),phinum= 0.6684E+00 0.1255E+01 0.1253E+01
r,phi(exact),phinum= 0.7518E+00 0.1166E+01 0.1164E+01
r,phi(exact),phinum= 0.8351E+00 0.1084E+01 0.1082E+01
r,phi(exact),phinum= 0.9184E+00 0.1010E+01 0.1009E+01
r,phi(exact),phinum= 0.1002E+01 0.9437E+00 0.9422E+00
r,phi(exact),phinum= 0.1085E+01 0.8835E+00 0.8820E+00
r,phi(exact),phinum= 0.1168E+01 0.8292E+00 0.8278E+00
r,phi(exact),phinum= 0.1252E+01 0.7802E+00 0.7788E+00
r,phi(exact),phinum= 0.1335E+01 0.7358E+00 0.7346E+00
r,phi(exact),phinum= 0.1418E+01 0.6957E+00 0.6945E+00
r,phi(exact),phinum= 0.1502E+01 0.6593E+00 0.6582E+00
r,phi(exact),phinum= 0.1585E+01 0.6262E+00 0.6251E+00
r,phi(exact),phinum= 0.1668E+01 0.5961E+00 0.5950E+00
r,phi(exact),phinum= 0.1752E+01 0.5685E+00 0.5675E+00
r,phi(exact),phinum= 0.1835E+01 0.5433E+00 0.5423E+00
r,phi(exact),phinum= 0.1918E+01 0.5201E+00 0.5192E+00
r,phi(exact),phinum= 0.2002E+01 0.4987E+00 0.4979E+00
r,phi(exact),phinum= 0.2085E+01 0.4790E+00 0.4782E+00
r,phi(exact),phinum= 0.2168E+01 0.4608E+00 0.4600E+00
r,phi(exact),phinum= 0.2252E+01 0.4438E+00 0.4431E+00
r,phi(exact),phinum= 0.2335E+01 0.4280E+00 0.4273E+00
r,phi(exact),phinum= 0.2418E+01 0.4133E+00 0.4127E+00
r,phi(exact),phinum= 0.2502E+01 0.3996E+00 0.3989E+00
r,phi(exact),phinum= 0.2585E+01 0.3868E+00 0.3861E+00
r,phi(exact),phinum= 0.2668E+01 0.3747E+00 0.3741E+00
r,phi(exact),phinum= 0.2752E+01 0.3634E+00 0.3628E+00
r,phi(exact),phinum= 0.2835E+01 0.3527E+00 0.3521E+00
r,phi(exact),phinum= 0.2918E+01 0.3426E+00 0.3421E+00
r,phi(exact),phinum= 0.3002E+01 0.3331E+00 0.3326E+00
r,phi(exact),phinum= 0.3085E+01 0.3241E+00 0.3236E+00
r,phi(exact),phinum= 0.3168E+01 0.3156E+00 0.3151E+00
r,phi(exact),phinum= 0.3252E+01 0.3075E+00 0.3070E+00
r,phi(exact),phinum= 0.3335E+01 0.2998E+00 0.2994E+00
r,phi(exact),phinum= 0.3418E+01 0.2925E+00 0.2921E+00
r,phi(exact),phinum= 0.3502E+01 0.2856E+00 0.2851E+00
r,phi(exact),phinum= 0.3585E+01 0.2789E+00 0.2785E+00
r,phi(exact),phinum= 0.3668E+01 0.2726E+00 0.2721E+00
r,phi(exact),phinum= 0.3752E+01 0.2665E+00 0.2661E+00
r,phi(exact),phinum= 0.3835E+01 0.2608E+00 0.2603E+00
r,phi(exact),phinum= 0.3918E+01 0.2552E+00 0.2548E+00
r,phi(exact),phinum= 0.4002E+01 0.2499E+00 0.2495E+00
r,phi(exact),phinum= 0.4085E+01 0.2448E+00 0.2444E+00
r,phi(exact),phinum= 0.4168E+01 0.2399E+00 0.2394E+00
r,phi(exact),phinum= 0.4252E+01 0.2352E+00 0.2347E+00
r,phi(exact),phinum= 0.4335E+01 0.2307E+00 0.2302E+00
r,phi(exact),phinum= 0.4418E+01 0.2263E+00 0.2258E+00
r,phi(exact),phinum= 0.4502E+01 0.2221E+00 0.2216E+00
r,phi(exact),phinum= 0.4585E+01 0.2181E+00 0.2176E+00
r,phi(exact),phinum= 0.4668E+01 0.2142E+00 0.2137E+00
r,phi(exact),phinum= 0.4752E+01 0.2105E+00 0.2099E+00
r,phi(exact),phinum= 0.4835E+01 0.2068E+00 0.2063E+00
r,phi(exact),phinum= 0.4918E+01 0.2033E+00 0.2028E+00
r,phi(exact),phinum= 0.4918E+01 0.2033E+00 0.1994E+00

Fig 2. V(r) exact and V(r) numerical obtained by solving Poisson's equation.

FORTRAN code for Poisson eq in atomic calculations

c potential through Poisson's equation may 13,2008
data nstep,a, epsi/3000, 2.,1.e-4/
c phinum is the numerical potential , Phi(r) is the anlytic solution
dimension phinum(0:5000)
psi1s(r)=sqrt(a**3/pi)*exp(-a*r)
phi(r)=(1./r)*(1.-exp(-2.*Z*r) ) - Z*exp(-2.*Z*r)
rho(r)=psi1s(r)**2
f1(r)=4.*pi*r*rho(r)
pi=2.*asin(1.)
Z=a
ri=epsi
rf=10./a
dr=(rf-ri)/float(nstep)
c calculates IC Phinum(r=0)
sum=0.
do 10 i=1,nstep
r=ri+dr*float(i)
sum=sum+(dr/2.)*(f1(r)+f1(r-dr))
10 continue
c initial conditions,the first derivative is zero; phinum(1)=phinum(0)
phinum(0)=sum
phinum(1)=phinum(0)
print*,'phinum(0)=',phinum(0)
c solution of Poisson eq by finite difference
do 20 i=2,nstep
r=ri+dr*float(i)
phinum(i)=2.*phinum(i-1)-phinum(i-2)+dr**2*( (-2./(r-dr))*
$ (phinum(i-1)-phinum(i-2))/dr -4.*pi*rho(r-dr) )
20 continue
c print phi theor and phi num
print 100 ,ri, phi(ri),phinum(0)
do 30 i=1,nstep,50
r=ri+dr*float(i)
print 100 ,r, phi(r),phinum(i)
30 continue
print 100 ,r, phi(r),phinum(nstep)
100 format(1x,'r,phi(exact),phinum=',3(3x,e11.4))
stop
end

Example of self consistent method

The first potential is taken to be (1./(r+epsi))*(1.-exp(-2.*a*r) ) - a*exp(-2.*a*r) with a=Z-5/16.

The next potential is obtained using Poisson's equation . Subroutine vr2 in the code , FORTRAN - SCF method , is employed. Then , using Mirosoft Excel a 3rd. degree polynomial is found to represent V(r).

(rlimit = 5.0 , nsteps = 4000 )

V(r) -e_i Ve(1s,1s) -E_total = - {2e_i -Ve(1s,1s)}

1. (1./(r+epsi))*(1.-exp(-2.*a*r) ) - a*exp(-2.*a*r) .9058 1.06870258 -2.88030243

(a= Z-5./16.)

2. -0.0342*r**3 + 0.3465*r**2 - 1.2342*r + 1.8519 .8936 1.08084798 -2.86804795

3.-0.0343*r**3 + 0.3471*r**2 - 1.2348*r + 1.8521 .8930 1.08147633 -2.86747646

In three runs we have reached agreement to three digits with the Hartree-Fock value given in ref. 4 as -2.81617 au.

The first ionization potential of helium is 24.580 V= 0. 9038 au which compares very well with the last -e_i = .8930 .

The various approximations used are:

a) finite difference method to solve all differential equations

b) integrations are done in the range 0 ≤ r ≤ 5.

c) after solving for V(r) employing Poisson's equation, it is approximated by a 3rd degree polynomial

FORTRAN - SCF method

dimension phinum(0:5000), psi(0:5000)
data nstep, epsi/4000,1.e-4/
data z,ei,ef,ne/2.,-.8946,-.90, 1/
psi1s(r)=sqrt(a**3/pi)*exp(-a*r)
pi=2.*asin(1.)
a=z-5./16.
200 e=ei
de=(ef-ei)/float(ne)
do 20 ie=1,ne
psi(0)=1.
psi(1)=psi(0)
rlimit=3.0
dr=rlimit/float(nstep)
c backward integration Psi(nstep) , psi(nstep-1)..
c psi(nstep)=exp(-sqrt(2.*abs(e))*rlimit)
c psi(nstep-1)=exp(-sqrt(2.*abs(e))*(rlimit-dr))
do 10 i=2,nstep
r=dr*float(i)
psi(i)=2.*psi(i-1)-psi(i-2) -(2.*dr/(r+dr))*(psi(i-1)-psi(i-2))
$-2.*dr**2*(z/(r-dr)-v(r-dr)+e)*psi(i-1)
10 continue
print 100,e,psi(nstep)
e=e+de
20 continue
if(ne.eq.1)call anorm(nstep,pi,rlimit,dr,psi)
print*,' '
if(ne.eq.1)call vr2(nstep,rlimit,dr,epsi,psi,phinum)
if(ne.eq.1)call aj1s1s(nstep,dr,pi,psi,ve1s1s)
print*,' '
print*,'ei,ve1s1s,etotal=',ei,ve1s1s ,2.*ei-ve1s1s
100 format('e,psi(nstep)/abs(psi)=',2(3x,e13.6))
stop
end

subroutine vr2(nstep,rlimit,dr,epsi,psi,phinum)
dimension phinum(0:5000),psi(0:5000)
rho(n)=psi(n)**2
f1(n,r)=4.*pi*r*rho(n)
h=abs(dr)
pi=2.*asin(1.)
c calculates IC Phinum(r=0)
sum=0.
do 10 i=1,nstep
r=dr*float(i)
sum=sum+(dr/2.)*(f1(i,r)+f1(i-1,r-dr))
10 continue
c initial conditions,the first derivative is zero; phinum(1)=phinum(0)
phinum(0)=sum
phinum(1)=phinum(0)
c print*,'phinum(0)=',phinum(0)
c solution of Poisson eq by finite difference
do 20 i=2,nstep
r=dr*float(i)
phinum(i)=2.*phinum(i-1)-phinum(i-2)+ dr**2*( (-2./(r-dr))*
$ (phinum(i-1)-phinum(i-2))/dr -4.*pi*rho(i-1) )
20 continue
do 30 i=0,nstep,70
r=dr*float(i)
print 100 ,r,phinum(i)
30 continue
print 100 ,dr*float(nstep),phinum(nstep)
100 format(1x,'r,phinum=',2(3x,e11.4))
return
end

subroutine anorm(nstep,pi,rlimit,dr,psi)
dimension psi (0:5000)
sum=0.
do 10 i=1,nstep
r=dr*float(i)
sum=sum+(dr/2.)*(psi(i)**2*4.*pi*r**2 +psi(i-1)**2
$ *4.*pi*(r-dr)**2 )
10 continue
factor=sqrt(sum)
do 20 i=0,nstep
psi(i)=psi(i)/factor
20 continue
sum=0.
do 30 i=1,nstep
r=dr*float(i)
sum=sum+(dr/2.)*(psi(i)**2*4.*pi*r**2 +psi(i-1)**2
$ *4.*pi*(r-dr)**2 )
30 continue
do 40 i=0,nstep,70
r=dr*float(i)
print 100,r,psi(i)
if(i.eq.nstep)print 100,r,psi(i)
40 continue
100 format('r,psi=',2(3x,e10.3))
c print*,'sum=',sum
return
end

subroutine aj1s1s(nstep,dr,pi,psi,ve1s1s)
dimension psi(0:5000)
f(r,n)= 4.*pi*r**2*v(r)*psi(n)**2
sum=0.
do 10 i=1,nstep
r=dr*float(i)
sum=sum+(dr/2.)*(f(r,i)+f(r-dr,i-1) )
10 continue
ve1s1s=sum
return
end

function V(r)
data z,epsi/2.,1.e-4 /
a=z-5./16.
c V = (1./(r+epsi))*(1.-exp(-2.*z*r) ) - z*exp(-2.*z*r)
c V = (1./(r+epsi))*(1.-exp(-2.*a*r) ) - a*exp(-2.*a*r)
c V = -0.0337*r**3 + 0.3432*r**2 - 1.2263*r + 1.8445
v = -0.0342*r**3 + 0.346*r**2 - 1.2314*r + 1.8483

return
end