close

Вход

Забыли?

вход по аккаунту

?

How to build an atom - Psi-k

код для вставки
How to build an atom
M.S.S. Brooks
September 29, 2009
Contents
1 Introduction
3
2 Elementary numerical techniques
5
2.1
A first attempt at the SchrВЁ
odinger equation for free electrons . . . . . . . . . . . . . . . . . . . .
5
2.2
Reduction of linear differential equations to a first order system . . . . . . . . . . . . . . . . . . .
6
2.3
A second attempt at the SchrВЁ
odinger equation for free electrons . . . . . . . . . . . . . . . . . . .
10
3 Hydrogen
12
3.1
Bound states in a Coulomb potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
3.2
Properties of the eigenstates of the spherical Coulomb potential . . . . . . . . . . . . . . . . . . .
13
3.3
The Pauli principle and the hydrogen model of the periodic system . . . . . . . . . . . . . . . . .
14
4 The real structure of the periodic table
15
4.1
Madelung rules for the periodic system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
4.2
Self-consistent charge densities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
16
5 Calculated properties of free atoms
18
5.1
Trends in ionicity and atom size
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18
5.2
Band masses of transition metals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
5.3
Exchange integrals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
6 Boundary conditions and series expansions for the SchrВЁ
odinger equation.
27
6.1
ScrВЁ
odinger equation for a spherical potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
6.2
Coupled first order radial equations - series expansion . . . . . . . . . . . . . . . . . . . . . . . .
28
6.3
Boundary conditions at large radius - the WKB approximation . . . . . . . . . . . . . . . . . . .
30
1
2
CONTENTS
7 Relativistic effects
32
7.1
The Dirac equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
32
7.2
Coupled first order radial Dirac equations in a central field . . . . . . . . . . . . . . . . . . . . .
33
7.3
The radial Pauli equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
7.4
Generic radial SchrВЁ
odinger, relativistic Pauli and Dirac equations . . . . . . . . . . . . . . . . . .
38
7.5
The spin polarised Dirac equation
40
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8 Virial theorem
42
8.1
Virial theorem from the SchrВЁ
odinger equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
42
8.2
Spherical charge density approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
44
8.3
Virial theorem from the Dirac equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
46
8.4
Kinetic energy from radial Dirac Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
47
9 Programs
50
9.1
One dimensional SchrВЁ
odinger equation - free electrons . . . . . . . . . . . . . . . . . . . . . . . .
50
9.2
Three dimensional SchrВЁ
odinger equation - free electrons . . . . . . . . . . . . . . . . . . . . . . .
50
9.3
Euler method for a system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
9.4
Runge-Kutta for a system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
52
9.5
Free electron SchrВЁ
odinger equation - revisited . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
52
9.6
Hydrogenic bound states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
54
9.7
Self-consistent field atom local density approximation
62
. . . . . . . . . . . . . . . . . . . . . . . .
3
1
Introduction
Not a great deal of attention is paid to atom cores or free atoms in solid state physics because there exist reliable
and much used programs that handle the core. Many of those programs were written when memory was scarce
and computers both slow and expensive. Every possible programming device was therefore used to achieve
accuracy whilst sparing memory and time. Although excellent, this meant that such programs could be hard to
follow. The modern personal computer has bucketfuls of memory and is very fast. There is therefore no need to
save time and memory when constructing free atom charge densities with a modern personal computer. There
are two features of the free atom problem that differentiate it from the physics of conduction electrons. Firstly,
one must find the discrete eigenvalues of bound states rather than a continuous energy spectrum. Secondly, there
are many branches (or principal quantum numbers) in a heavier atoms and we will need to find the eigenvalues
for them all.
In order to build a free atom it is going to be necessary to solve the Dirac equation in a self-consistent field
due to the other electrons. It is desirable to do this as simply and quickly as possible but in order to make the
process as comprehensible as possible and to sort out which properties of free atoms are due to single electron
solutions in a spherical potential, which are due to the self-consistent field, and which are due to the effects of
relativity, we will go through the process in a series of steps. So we shall begin at the very beginning.
Firstly, it will be necessary to solve some differential equations. Let’s take the one dimensional free particle
SchrВЁ
odinger equation as an example
2m
d2 П€
+Оµ 2 П€ =0
2
dx
(1.1)
The mass of the electron is m = 9.1 × 10−28 gm, it’s charge is e = 1.6 × 10−19 C and the reduced Planck constant
is = 6.58 × 10−16 eV.s. Let’s straight away, choose to use atomic Rydberg (Ry) units. Starting with Gaussian
units one sets
e2
= 2m =
=1
(1.2)
2
where m is the mass of an electron
and e its charge. The unit of angular momentum is , the unit of mass is 2m
в€љ
and the unit of charge is e/ 2. The unit of length is the Bohr radius
2
a0 в‰Ў
me2
= 0.5291 Г— 10в€’8 cm
The unit of energy is the Rydberg
ERy в‰Ў
e2
me4
=
= 2.18 × 10−11 erg ≈ 13.6eV
2a0
2 2
The unit of time is ratio of angular momentum to energy
t0 в‰Ў
ERy
= 4.84 Г— 10в€’17 s
The unit of velocity is the ratio of length to time
a0
e2
=
≈ 1.094 × 108 cm/s
t0
2
4
1
INTRODUCTION
and the velocity of light in vacuum, c, is1
c≈
2.998 Г— 1010 cm/s
≈ 274
1.094 Г— 108 cm/s
In atomic units, Equation (1.1) becomes
d2 П€
d2 П€
+ ОµП€ =
+ k2 П€ = 0
2
dx
dx2
(1.3)
в€љ
where k = Оµ. Of course we know that the solutions are oscillatory for Оµ > 0, exp(В±ikx), with the coefficients
A and B in
П€ = Aeikx + Beв€’ikx
determined by the boundary conditions. For example, if there is just a single boundary condition П€(0) = 0, the
solutions are в€ќ sin kx, for any value of k, with the remaining constant determined by choice of normalisation. A
second spatial boundary condition will restrict solutions to discrete values of k.
Let’s start by looking briefly at how second order differential equations like Equation (1.3) are solved numerically.
1 Alternatively,
one may use the relationship between c and the fine structure constant, О±,
α≈
1
e2
2
=
=
137
c
c
в†’
c=
2
≈ 274
О±
5
2
2.1
Elementary numerical techniques
A first attempt at the SchrВЁ
odinger equation for free electrons
Now the simplest, and least accurate, way to solve Equation (1.1) is to follow Euler, set up a linear grid with
step size, h, for the x-axis and make the approximations
dП€
П€i+1 в€’ П€i
≈
dx
h
and
d dП€
d2 П€
d П€i+1 в€’ П€i
П€i+1 + П€iв€’1 в€’ 2П€i
≈
=
≈
dx2
dx dx
dx
h
h2
Equation (1.3) is therefore approximated by
П€i+1 + П€iв€’1 в€’ 2П€i
+ П€i = 0
h2
or, if i is shifted to i + 1
П€i+2 = (2 в€’ h2 Оµ)П€i+1 в€’ П€i
(2.1)
This looks easy to solve if one has the initial values for П€1 and П€2 . For example the program listed in 9.1 took
a few minutes to write and produces the results shown on the left hand side of Fig. 2.1. What this bit of match
practice shows is that we have calculated unnormalised sin and cosin functions in a rather silly way. However one
does see clearly from this example firstly how the choice of initial values determines the solution and, secondly,
that the solutions are not normalised.
Figure 2.1: Approximate unnormalised solutions to the one dimensional SchrВЁodinger equation (left) and radial
SchrВЁ
odinger equation (right) for free electrons with an energy of 1Ry.
Now let’s move to three dimensions but stay with free electrons. The Scr¨odinger equation becomes
∇2 + εψ = 0
Since
∇2 =
1 ∂2 2
L2
r
в€’
2 r2
r2 ∂r2
(2.2)
(2.3)
Equation (2.2) becomes
в€’
L2
1 ∂2
r
в€’
в€’Оµ П€ =0
r ∂r2
2mr2
(2.4)
6
2
ELEMENTARY NUMERICAL TECHNIQUES
Since there is no potential, the Hamiltonian has at least spherical symmetry, L2 = L В· L = L2 is scalar, and the
SchrВЁ
odinger equation is separable with solutions of the form
П€(r) = il Rl (r)Ylm (Л†
r)
whence
(2.5)
l(l + 1)
d2
2 d
1 d2
l(l + 1)
r
в€’
+
Оµ
R
(r)
=
+
+ Оµ Rl (l) = 0
в€’
l
r dr2
r2
dr2
r dr
r2
(2.6)
and the three dimensional ScrВЁ
odinger equation has been reduced to a one-dimensional equation for the radial part
of the wave function, Rl (r). We may therefore solve Equation (2.5) in the same way that we solved Equation (1.3),
it’s just a little more complicated
П€i+1 2 1 +
П€i+2 =
h
ri+1
в€’ h2 Оµ в€’
1+
l(l + 1)
2
ri+1
в€’ П€i
(2.7)
2h
ri+1
leading to the program listed in 9.2 (where l has been set equal to zero), which produces the results shown on the
right hand side of Fig. 2.1. This time we have calculated unnormalised spherical Bessel and Neumann functions
in a silly way and not very accurately, but the fact is that second order differential equations are not solved in
this way since there are more convenient and accurate methods. It’s now time to look at the methods used in
practice.
2.2
Reduction of linear differential equations to a first order system
We want to solve the second-order differential equation
y = f (x, y, y ),
with initial values
y(x0 ) = y0 ,
y (x0 ) = z0
(2.8)
If we define z в‰Ў y , Equation (2.8) is equivalent to the first-order system
y
= z
z
= y = f (x, y, z)
(2.9)
with z(x0 ) = z0 2 . From now on we shall solve second order differential equations by first transforming them to
a first order coupled system.
We now need to address the question of accuracy. The best way to demonstrate how the Euler formula might be
improved is to use an example where we know the analytic result, such as the following. A simple second order
differential equation is
y =x+y
(2.10)
2 Of
course any n’th order equation may be reduced to a system of first order equations. Thus, for
y (n) = f (x, y, y , y (2) , . . . , y (nв€’1) )
we just define new functions
z1 = y ,
z2 = y (2) = z1
z3 = y (3) = z2
znв€’1 = y (nв€’1) = znв€’2
then
zn = y (n) = f (x, y, z1 , z2 , . . . , znв€’1 )
making a total of n first order coupled equations. For second order differential equations, n = 2, z1 = z and z2 = z .
2.2
Reduction of linear differential equations to a first order system
7
with the boundary conditions y(0) = 0 and y (0) = 1. Equation (2.10) reduces to the first-order system
y
= z
z
= x+y
(2.11)
with y(0) = 0 and z(0) = 1. The general analytic solution to the Equation (2.10) is
y = в€’x + c1 ex + c2 eв€’x
(2.12)
Due to the boundary condtion y(0) = 0, c1 = в€’c2 = c and
y = в€’x + c(ex в€’ eв€’x )
Due to the boundary condition y (0) = 1, c = 1 and
y = в€’x + (ex в€’ eв€’x )
and
y(x = 1) = в€’1 + e в€’ 1/e = 1.350402
(2.13)
So we know that y(x = 1) = 1.350402 and we can take a look at a couple of approximations for numerical
integration. The simplest numerical integration is the Euler method which is, for a step distance h,
yn+1
= yn + hyn = yn + hzn
zn+1
= zn + hzn = zn + h(xn + yn )
(2.14)
A program that uses the Euler method for a system is listed in 9.3. The results for increasing numbers of grid
points are:
Euler integration of z’=y+x and z=y’ system
grid points=
11 x=
1.000000 y=
1.245064
grid points= 211 x=
1.000000 y=
1.344836
grid points= 411 x=
1.000000 y=
1.347544
grid points= 611 x=
1.000000 y=
1.348479
grid points= 811 x=
1.000000 y=
1.348954
grid points= 1011 x=
1.000000 y=
1.349240
grid points= 1211 x=
1.000000 y=
1.349432
grid points= 1411 x=
1.000000 y=
1.349570
grid points= 1611 x=
1.000000 y=
1.349673
grid points= 1811 x=
1.000000 y=
1.349754
grid points= 2011 x=
1.000000 y=
1.349818
grid points= 2211 x=
1.000000 y=
1.349871
grid points= 2411 x=
1.000000 y=
1.349915
z=
z=
z=
z=
z=
z=
z=
z=
z=
z=
z=
z=
z=
1.942421
2.078840
2.082405
2.083635
2.084258
2.084635
2.084887
2.085067
2.085203
2.085309
2.085394
2.085463
2.085521
The convergence towards y = 1.3504 is tardy. In fact the Euler approximation is not very good. From the first
order differential equation y = f we know that
xi+m
yi+m = yi +
f (y, x)dx
(2.15)
xi
The Euler formula, which is the simplest and least accurate, corresponds to setting m = 1 and approximating
f (y, x) ≈ fi . Numerical accuracy depends upon improved approximations to the integral on the right of Equation (2.15), usually a Taylor polynomial approximation. The most used are the predictor corrector methods
where a simple formula is used to predict the next value of y and a more accurate corrector formula is used to
8
2
ELEMENTARY NUMERICAL TECHNIQUES
provide successive improvements. To see how this works let’s start with the simple Euler formula for a system
with y = f
yi+1 = yi + hyi = yi + hfi
which requires just yi and fi and predicts yi+1 . However the more accurate modified Euler formula is
yi+1 = yi +
h
h
y + yi+1 = yi + (fi + fi+1 )
2 i
2
which requires fi+1 = yi+1 as well. The predicted value of yi+1 may be used to estimate fi+1 = yi+1 so that the
modified Euler formula may be used to correct yi+1 . Successive corrections to yi+1 and yi+1 may then be made.
Two much used predictor corrector pairs are, the Milne pair
4h
2yiв€’2 в€’ yiв€’1 + 2yi
3
h
=yiв€’1 +
y
+ 4yi + yi+1
3 iв€’1
yi+1 =yiв€’3 +
(2.16a)
yi+1
(2.16b)
where Simpson’s rule will be recognised, and the Adams pair
h
в€’9yiв€’3 + 37yiв€’2 в€’ 59yiв€’1 + 55yi
24
h
=yi +
y
в€’ 5yiв€’1 + 19yi + 9yi+1
24 iв€’2
yi+1 =yi +
(2.17a)
yi+1
(2.17b)
Predictor corrector methods can be very accurate but may require a large number of initial values, which are
obtained from power series expansions.
An alternative is the Runge-Kutta method which avoids the computation of higher order derivatives and therefore
requires only one initial value. In place of the derivatives extra values of the function f must be provided to
duplicate the accuracy of the Taylor polynomial. To see how this works we will lay out the detail of the second
order Runge-Kutta method with y = f . The Taylor expansion of y up to second order is
y(x + h) = y(x) + hy (x) +
h2
h2
h2
y (x) = y(x) + hf (x) + f (x) = y(x) + hf (x) +
(fx + f fy ))
2
2
2
(2.18)
where we have used
∂f
∂y ∂f
df
=
+
= fx + f fy
dx
∂x ∂x ∂y
Another expansion of y up to second order is
y(x + h) =y(x) + a1 k1 + a2 k2
where
and
(2.19a)
k1 =hf (x, y)
(2.19b)
k2 =hf (x + bh, y + bk1 )
(2.19c)
with a1 , a2 and b to be evaluated. This is done by expanding k2 , Equation (2.19c), in a Taylor series
k2 = hf + h2 bfx + h2 bf fy = hf + h2 b(fx + f fy )
and substituting the result into Equation (2.19a), which becomes
y(x + h) = y(x) + h(a1 + a2 )f + h2 a2 b (fx + f fy )
(2.20)
The coefficients are evaluated by comparing Equation (2.18) with Equation (2.20), but they are underdetermined
and there is room for choice. One choice would be
a1 = a2 =
1
2
b=1
2.2
Reduction of linear differential equations to a first order system
9
from which it follows that
1
y(x + h) =y(x) + (k1 + k2 )
2
k1 =hf (x, y)
(2.21b)
k2 =hf (x + h, y + k1 )
(2.21c)
(2.21a)
Of course there are no free rides. The Runge-Kutte method has simulated a power series expansion, but does so
at the cost of requiring the function at intermediate values of the initial step-size. This will require interpolation
of a tabulated function if the analytic form is not known. In atomic calculations this requires interpolation of the
potential. Most frequently used is the fourth order Runge-Kutta method (popularly known as RK4) which is
obtained from a Taylor expansion of order four. For a system of two coupled first order equations it is, explicitly,
yn+1
=
zn+1
=
1
yn + (k1 + 2k2 + 2k3 + k4 )
6
1
zn + (l1 + 2l2 + 2l3 + l4 )
6
(2.22)
where
k1
= hyn = hf (xn , yn , zn )
l1
k4
= hzn = hg(xn , yn , zn )
1
1
1
= hf (xn + h, yn + k1 , zn + l1 )
2
2
2
1
1
1
= hg(xn + h, yn + k1 , zn + l1 )
2
2
2
1
1
1
= hf (xn + h, yn + k2 , zn + l2 )
2
2
2
1
1
1
= hg(xn + h, yn + k2 , zn + l2 )
2
2
2
= hf (xn + h, yn + k3 , zn + l3 )
l4
= hg(xn + h, yn + k3 , zn + l3 )
k2
l2
k3
l3
(2.23)
Thus, yn+1 is determined by yn plus the product of the size of the interval (h) and an estimated slope. The
estimated slope is a weighted average of slopes: k1 is the slope at the beginning of the interval, k2 is the slope
at the midpoint of the interval, using slope k1 to determine the value of y at the point xn + h/2 using Euler’s
method. k3 is also the slope at the midpoint using the slope k2 to determine the corrsponding value of y. k4
is the slope at the end of the interval, with the corresponding value of y determined using k3 . In averaging the
four slopes, greater weight is given to the slopes at the midpoint. The RK4 method is a fourth-order method,
meaning that the error per step is O(h5 ), with the accumulated error is O(h4 ). Since it requires initial values
only at the first grid point it provides good value for money. A program that integrates Equation (2.10) is listed
in 9.4 and the results are:
Runga-Kutta integration of z’=y+x and z=y’ system
grid points=
9 x=
1.000000 y=
1.350397 z=
grid points= 14 x=
1.000000 y=
1.350402 z=
grid points= 19 x=
1.000000 y=
1.350402 z=
grid points= 24 x=
1.000000 y=
1.350402 z=
grid points= 29 x=
1.000000 y=
1.350402 z=
2.086157
2.086161
2.086161
2.086161
2.086161
the improvement upon the Euler formula is so drastic, with the correct answer to 6 decimal places by 14 grid
points, that we should all be suitably impressed.
10
2
2.3
ELEMENTARY NUMERICAL TECHNIQUES
A second attempt at the SchrВЁ
odinger equation for free electrons
We now return to the radial free electron SchrВЇodinger equation to get a more accurate solution. Firstly, we will
replace the second order radial SchrВЇ
odinger equation by two coupled first order equations, and then add a couple
of new tricks. We also note that Equation (2.5) leads to the normalisation condition
dr |Ylm (Л†
r |2 Rl2 (r) =
П€|П€ = 1 =
dr r2 Rl2 (r)
(2.24)
and that Equation (2.6) may also be written
l(l + 1)
d2
l(l + 1)
d2
в€’
+ E rRl (r) =
в€’
+ k 2 rRl (r) = 0
2
2
2
dr
r
dr
r2
where k =
в€љ
(2.25)
Оµ. If we now let Pl = rRl and define Q through
Q=P в€’
(l + 1)
P
r
(2.26)
It is easy to show, using Equations (2.25) and (2.6), that
Q =P +
(l + 1)
(l + 1)
(l + 1)
Pв€’
P = в€’k 2 P в€’
Q
2
r
r
r
Consequently the coupled equations
(l + 1)
P +Q
r
(l + 1)
Q =в€’
Q в€’ k2 P
r
P =
(2.27a)
(2.27b)
are equivalent to Equation (2.25). So the first trick is to always use Pl = rRl , rather than Rl . The second trick,
which greatly improves accuracy, is to use a logarithmic mesh on the r-axis. Since most of the activity occurs
ar small radii (the kinetic energy is larger and the wave functions have greater curvature), we define x through
r в€ќ ex making the mesh of x dense near the origin. For numerical work with a logarithmic mesh one should
note that r = r0 ex therefore dr/dx = r, dr = rdx3 . Consequently, rdP/dr = dP/dx, rdQ/dr = dQ/dx, and
Equations (2.27) take the form, with y1 = P and y2 = Q,
dy1
dx
dy2
dx
= ay1 + by2
= в€’ay2 + dy1
(2.29)
with
a=l+1
b=r
d = в€’rОµ
(2.30)
The second equation, for r → 0, requires that if y1 (1) = 1 then y2 (1) = 0. The routine in 9.5 solves Equation (2.29) and normalises the wave functions with the function simf using Simpson’s rule. The results are
shown in Fig. 2.2 for l = 0, . . . , 6 where they may be recognized as spherical Bessel functions, apart from
normalization. The nodes of the Bessel functions are listed in table 2.3.
3 Integration;
R
R
dr φ = dx rφ, therefore, if the function to be integrated is simply multiplied by r, a normal integration routine for
a constant stepsize may be used. Differentiation; dφ/dr = (1/r)dφ/dx. If an atomic potential, φ, is be differentiated or interpolated,
rφ is slowly varying, approaching constant for a Coulomb potential. It is therefore useful to use
d
dφ
(rφ) = rφ + r2
в†’
dx
dr
A standard differentiation routine may then be used for d(rφ)/dx
dφ
=
dr
d
(rφ) −
dx
r2
rφ
(2.28)
2.3
A second attempt at the SchrВЁ
odinger equation for free electrons
n
1
2
3
4
nodes
0
1
2
3
l=0
3.14
6.28
9.42
12.56
l=1
4.49
7.73
10.9
14.1
l=2
5.76
9.10
12.3
l=3
6.99
10.42
13.7
11
l=4
8.18
11.7
l=5
9.36
12.96
Table 2.1: The nodes of the calculated spherical Bessel functions.
Figure 2.2: Radial solutions to the SchrВЁ
odinger equation in a spherical potential set to zero. The radial functions
are normalized by dr r2 jn2 (r) = 1 which differs from the normal normalization for Bessel functions (j0 (0) = 1).
As long as the energy is positive a potential, for example a Coulomb potential в€’2Z/r, may be added and the
radial equations remain the same except that d becomes d = (V в€’ Оµ)r and the boundary condition at r в†’ 0
becomes y2 = −zy1 /(l + 1). If the potential has finite range, it’s effect at large distances is just to phase shift
the solutions.
There are scattering solutions for all values of positive energy. The electrons making up free atoms, in contrast,
are all in bound states with negative energies. Consequently if one integrates the SchrВЁodinger equation outwards
a turning point will be reached where the electron energy becomes less than the potential. The wave functions
then change from oscillatory to exponential. Now there are two exponetial solutions, proportional to ear and
eв€’ar . One of them blows up at large r for an arbitrary value of the energy. The eigenvalues are at those energies
where the solution decreases exponentially at large radii. The way to find the eigenvalues therefore, is to integrate
outwards in the classically allowed region and to integrate inwards in the classically forbidden region, ensuring
the proper boundary conditions are imposed at both small and large radii. The inward and outward solutions
must then be matched at some point such that the entire solution is continuous and differentiable. This can
only occur at the eigenvalues. We therefore have to learn how to do this for a spherical potential and a simple
available example, for which we know the exact solutions, is the hydrogen atom with the Coulomb potential
equal to в€’2/r.
12
3
3.1
3
HYDROGEN
Hydrogen
Bound states in a Coulomb potential
Not only is hydrogen the universe’s most abundant element, but it’s electronic eigenstates are also known
analytically, which is more than be said for any other element. Since the nucleus is 103 times heavier than it’s
electron, we can treat the electron as a single particle moving in a potential в€’2/r. For all other atoms the
interaction between the electrons prevents analytic solutions. Later, we shall we shall be using density functional
theory (Hohenberg & Kohn 1964) to approximate the effect of electron interactions, with every element treated
similarly. That is, the relationship between potential and charge and spin density will be the same functional
of charge (and spin) density, regardless of Z. Curiously, the electronic structure of hydrogen then becomes far
more difficult to calculate - as we shall see. If one uses the more general theory and makes the self-interaction
correction for hydrogen, one returns to the single electron potential в€’2/r.
In the Coulomb potential, all the eigenstates with negative energy are bound states and the lowest bound state
will have an energy greater than в€’(2/n)2 , where n is the principle quantum number. The eigenvalues therefore
lie between в€’(2/n)2 and zero for a given branch, n. The program that deals with the hydrogenic atom is listed
in 9.6. The potential and the energy of a trial energy are plotted in Fig. 3.1. An outward integration, started
near r = 0, will be oscillatory up to X. It will have a number of nodes given by node = n в€’ l в€’ 1 which, for the
1s-state is zero, for the 2s-state is one, for the 2p-state is zero, etc. Beyond X there will be no nodes. This means
Figure 3.1: A trial energy in a Coulomb potential. For zero angular momentum the point X is the classical
turning point. For r < X the solutions are oscillatory and for r > X the physical solution falls off, becoming
exponential at large distances.
that we may determine the number of nodes by integrating the wave function out from r0 to X, labelled by ntp1,
3.2
Properties of the eigenstates of the spherical Coulomb potential
13
and counting the zero’s. We determine the energy range of the branch that we require by the shooting method.
The energy of integration is set in the subroutine getei g which calls the subroutine rkout to get the wave function
and count the nodes of the wave function. At the beginning the two bounding energies are ET = +40Ry4 and
EB = в€’(2Z/n)2 . A trial energy is set by (ET + EB)/2 and, if the number of nodes is too large, this becomes
ET . Otherwise it becomes EB. This process is continued until a range of 0.2Ry is found in the correct branch.
The subroutine gettp is used to find the turning point, X, and rkout is used to find the inner wave function up
to X. The subroutine rkin is used to integrate the wave function from a large value of the radius, labelled by
ntp2, inwards to X. Details the boundary conditions at small radius and X will be given later5 .
Now the outward and inward solutions don’t match at X. However, we are solving two coupled linear equations,
for P and Q and the boundary conditions only set the ratio of Q to P at r0 and X. We may therefore scale Pin
and Pout so that these are equal to unity at X, and scale Qin and Qout by the same amount. Qin and Qout are
then not continuous at X. Now it’s easy to see from the coupled equations that, if Qout > Qin , the eigenvalue
is too low, and vice versa for Qout < Qin . In fact it’s possible to predict the eigenvalue approximately using a
Taylor series expansion of the coupled equations, and this is usually done. But it’s not always entirely stable and
not worth it. We can simply shoot again. If Qout > Qin , one sets E = EB, and, if Qout < Qin , one sets E = ET
and tries a new energy at (ET + EB)/2. The process is repeated until the |ET в€’ EB| is less than some given
tolerance (in the program tol2 = 10в€’12 ). This method provides the prototype for heavier atoms than hydrogen.
3.2
Properties of the eigenstates of the spherical Coulomb potential
The one electron in hydrogen moves in the nuclear potential в€’1/r. Since the potential is spherically symmetric
the wave functions are products, Rnl Ylm , of radial functions and spherical harmonics. The characteristics of the
solutions are as follows.
(1) The principal quantum number, n, can take any integer value ≥ 1 and the energy, equal to −1/n2 , depends only upon this quantum number.
(2) The angular momentum quantum number, l, ranges from 0 to n в€’ 1, and the number of nodes in the
radial function is equal to n в€’ l в€’ 1. Since no wave function can have a negative number of nodes, the first
branch for any given angular momentum, l, has a principal quantum number n = l + 1 - e.g. the first d-states
(l = 2) are 3d. For any given angular momentum the number of nodes increases monotonically from zero with
increasing principal quantum number. The origin of this behaviour is the requirement that the solutions of
the wave equation be orthogonal. For a given angular momentum the orthogonality is not guaranteed by the
orthogonality of the spherical harmonics therefore the radial wave functions must be orthogonal. The radial
wave function of the higher energy state must therefore change sign one extra time. Additional nodes require
more curvature in the radial wave function corresponding to larger kinetic energy. They are also responsible for
an increase in the size of atoms as the electron density is pushed out.
4 This energy is, of course, far too large but it doesn’t matter. We shall later find that, for heavier atoms, it will be necessary to
guess an initial charge density which might be so bad that it makes the energies of occupied states positive. A restriction to negative
energies would then be a disaster.
5 The reason for starting the inward integration at r(ntp2) is that if one chooses a very large radius for the atom and starts the
inward integration very far out, an exponentially increasing solution, going inwards, will blow up before X is reached. To combat
this one starts integrating inwards from r(ntp2) and attaches an exponential tail, going outwards, beyond r(ntp2). This is all done
in rkin.
14
3
HYDROGEN
(3) The magnetic quantum number, m, ranges from в€’l to l and the degeneracy of the shell is6
nв€’1
l
nв€’1
(2l + 1) = n2 .
1=
l=0 m=в€’l
(3.1)
l=0
2
The probability densities, Pnl = r2 Rnl
, of the first three hydrogen wave functions are shown in Figs. 3.2 and 3.3.
The node in the 2s wave function has the effect of pushing the electron density away from the nucleus compared
with the 1s density. The 2p state has no nodes and would be closer to the nucleus than the 2s state except that
is has angular momentum and therefore a concomitant centrifugal potential, l(l + 1)/r2 , that repels the electron
density away from the nucleus. Compared to the 2s state, the inner part of the probability density of the 2p
state is further away from the nucleus (making it less penetrating) but the outer part of the probability density is
closer to the nucleus. The net effect is that the probability density (and therefore charge density) is compressed
increasingly into a ring as the angular momentum increases. Therefore the higher angular momentum states are
both less penetrating and more localised than lower angular momentum states with the same principle quantum
number.
2
Figure 3.2: The 1s, 2p and 2s radial probability densities Pnl = r2 Rnl
for hydrogen.
3.3
The Pauli principle and the hydrogen model of the periodic system
The Pauli principle7 states that no fermion state may contain more than one electron of a given spin and suggests
a rule for constructing the periodic system. It is postulated that heavier atoms may be built from hydrogen
by increasing the nuclear charge and the electron number such that they are always equal. The approximation
involved for the hydrogen model is that the potential for each electron is from the nuclear charge screened
by the other electrons leaving an effective nuclear charge of one. Then the the wave functions and energies are
approximately hydrogenic, and the energy of any shell defined by the principle quantum number, n, is E = в€’1/n2
with degeneracy equal to 2n2 when allowance is made for spin. The periodic table, constructed in this way, is
shown in Fig. 3.4.
The first row, for the 1s state, is correct. The second row is almost correct except that the 2s and 2p states
should not be degenerate as they are in the model. In reality the 2s, 2p degeneracy is lifted, the 2p states lying
6 Spin
7A
has not yet been introduced. Pauli degeneracy increases the total degeneracy to 2n2 .
thorough discussion of the relationship between spin and statistics is to be found in (Streater & Wightman 1978)
15
2
Figure 3.3: The 1s, 2s and 3s, 2p and 3d radial probability densities Pnl = r2 Rnl
for hydrogen.
higher in energy and becoming occupied later. The hydrogenic n = 2 shell is broken into 2s and 2p shells with
periods 2 and 6. The degeneracy of the hydrogenic model is also lifted in the third and subsequent rows. For
example, the 3d states have less binding energy the 3s and 3p states.
The origin of the removal of the degeneracy for a given branch is incomplete screening of the nuclear charge by
the electron charge density. The lower angular momentum states penetrate closer to the nucleus than higher
angular momentum states and therefore experience a potential that is less well screened by the charge density.
The binding energies of the lower angular momentum states are increased relative to those of the higher angular
momentum states, removing the degeneracy of the hydrogenic model. The periodic system therefore contains
periodicities of 2(2l + 1) (including spin) for each l rather than the full 2n2 degeneracy of the nth branch.
4
4.1
The real structure of the periodic table
Madelung rules for the periodic system
The rules from which the periodic table may be correctly built are due to Madelung:
1. First arrange the energy levels according to the hydrogenic approximation.
2. Look up the value of n + l as shown in Fig 4.1.
3. The degeneracy of the hydrogenic approximation is removed by filling the states according to the Pauli principle in order of increasing n + l then, for a given n + l, in order increasing n - subject to the requirement that
l < n.
For example, n + l = 1 allows only n = 1, l = 0 - the 1s state. If n + l = 2 only n = 2, l = 0 or 2s is allowed (
n = 1, l = 1 has n = l). If n + l = 3 (n = 1 l = 2 is not allowed) n can be 2 or 3 in which case l can be (a) 0
for n = 3, 3s state, (b) 1 for n = 2, 2p state and, following the third rule, the lowest value of n fills first, ie. 3p,
therefore the ordering so far is 1s, 2s, 2p, 3s. For n + l = 4, 3p and 4s are allowed and for n + l = 5, 3d, 4p and
5s are allowed and are filled in that order. The process is easily visualised through Fig 4.1. Comparison with
the real periodic table shows that the proposed ordering is correct for all known elements. We will be using the
Madelung rules to build an atom when we set up a table of the orbitals and their occupation numbers. .
16
4
THE REAL STRUCTURE OF THE PERIODIC TABLE
n = 5
g
n = 4
f
n = 3
d
n = 2
p
n = 1
s
1 s
2 p
3 d
2 s
3 s
1 s
2 p
3 d
2 s
3 s
Figure 3.4: The periodic table based upon the degeneracies of the hydrogen atom.
4.2
Self-consistent charge densities
Now we are must try to construct an atom with more than one electron and we will cut through the vast literature
on the subject of electron-electron interactions and use the simplest possible approximation. This is the local
density approximation (Kohn & Sham 1965) to density functional theory(Hohenberg & Kohn 1964)from which
there is no change to the SchrВЁ
odinger for hydrogen equation except that the potential for any electron is given
by
в€ћ
2Z
(r )
V (r) = в€’
+2
dr
+ Вµxc (r)
(4.1)
|r в€’ R|
|r
в€’r |
0
where is the electron charge density, r is the electron coordinate and R is the nuclear coordinate, which may
be taken to be the origin. Many approximations are available for Вµxc , which is a function of the electron density.
Typical is that due to von Barth & Hedin (1972), which we shall use here. An enormous simplification is made
by taking spherical averages in which case Equation (4.1) becomes
V (r) = в€’
2
2Z
+
r
r
в€ћ
r
(r )2 dr
0
(r )2 dr
(r ) + 2
r
(r )
+ Вµxc (r)
r
(4.2)
Now the problem is that we don’t know the electron density because we have not yet calculated it. Let’s go back
to hydrogen. If we had an electron density (the input density) we would just modify the hydrogen program,
listed in 9.6, by adding the extra three terms to the potential from Equation (4.2). We could then calculate
the eigenvalue and wave function of the 1s orbital. But we don’t have an input density so we shall have to
start with a guess or an estimate. The standard procedure is to estimate the electron density from the ThomasFermi approximation. The resulting eigenvalue and wave function are then incorrect, not least because the
Thomas-Fermi approximation for the density is poor (there is no shell structure). But at least we have a poor
approximation to the 1s wave function and we can recalculate the density from it. This density (the output
density) is also incorrect.
Long before formal density functional theory provided a rigorous basis for self-consistent calculations it was
known that if one could contrive to produce consistent input and output densities, the total electron density
4.2
Self-consistent charge densities
17
n=1
l=0 l=1 l=2 l=3 l=4 l=5 l=6 l=7
Filling sequence: n+l with l< n
1s
n+l
1
n=2
2s
2p
and increasing l
2,3
n=3
3s
3p 3d
n=4
4s
4p 4d
4f
n=5
5s
5p 5d
5f
5g
n=6
6s
6p 6d
6f
6g
6h
n=7
7s
7p 7d
7f
7g
7h
7i
n=8
8s
8p 8d
8f
8g
8h
8i
3,4,5
4,5,6,7
8j
Figure 4.1: Aufbau of the periodic table. The diagonals are constant n + l.
would be correct and the total energy at a variational minimum. The procedure is therefore to mix some of
output to the input density and try again, and to keep trying until what comes out is what goes in to some
given tolerance. In general this procedure is not stable but if the mixing ratio is kept low convergence is always
achieved for atoms. It is therefore a simple matter to add a self-consistent loop to the hydrogen program. There
are a multitude of schemes for accelerating convergence but, given the speed of modern personal computers,
convergence is so fast for free atoms with a fixed ratio that they hardly worthwhile. If one does this the results
for hydrogen are
EIGENVALUE=
TOTAL ENERGY=
-0.499109931604861
-0.923626197449808
which seems like a large step backwards. That the eigenvalue is not equal to the total energy is a property of
density functional theory. The eigenvalue does not have the same meaning in density functional theory as it
does for self-interaction corrected (SIC) hydrogen where it is also the total energy. More depressing is the value
of the total energy, which should be equal to в€’1Ry. The reason for the large error is that we have not allowed
any spin polarisation. The approximation that we have made in calculating Вµxc is that there is one half of an
electron with spin up and one half of an electron with spin down. This somewhat unphysical approximation for
an unpaired electron leads to a Coulomb interaction between two halves of an electron that appears in the total
energy. Later we will see that, if the electron is spin polarised, a single electron with spin up almost does not
interact with itself, and the total energy is almost в€’1Ry.
Nevertheless, we shall struggle on. Density functional theory states that the relationship between electron density
and potential has the same structure, independent of the value of Z. So we shall have a crack at many-electron
atoms, which requires that we introduce a table of orbitals to occupy with exactly Z electrons. We can use the
eigenstates of hydrogen to label these orbitals, 1s, 2s, 2p, 3s, 3p, 3d, . . .. We need to list the orbitals and to give
them names, which could be 1s, 2s, 2p, 3s, 3p, 3d, . . . and perhaps add a label for spin since in a spin polarised
calculation the spin up and spin down orbitals are no longer degenerate. However, we shall find that in relativistic
calculations the eigenvalues of states for a given value of l are split by spin-orbit interaction into states labelled
by total angular momentum j = l В± 1/2 therefore, with the wisdom of hindsight, a table of orbitals might look
like:
C TABLE OF ORBITALS
18
5
CALCULATED PROPERTIES OF FREE ATOMS
CHARACTER*8 NLJ(LDIM)
DATA NLJ/’ 1S1/2’,’ 2S1/2’,’ 2P1/2’,’ 2P3/2’,’
+ ’ 3P1/2 ’,’ 3P3/2 ’,’ 4S1/2 ’,’ 3D3/2 ’,’
+ ’ 4P1/2 ’,’ 4P3/2 ’,’ 5S1/2 ’,’ 4D3/2 ’,’
+ ’ 5P1/2 ’,’ 5P3/2 ’,’ 6S1/2 ’,’ 4F5/2 ’,’
+ ’ 5D3/2 ’,’ 5D5/2 ’,
+ ’ 6P1/2 ’,’ 6P3/2 ’,’ 7S1/2 ’,’ 5F5/2 ’,’
+ ’ 6D3/2 ’,’ 6D5/2 ’,’ 8S1/2 ’,’ 5G7/2 ’,’
+ ’ 6F5/2 ’,’ 6F7/2 ’,
+ ’ 7D3/2 ’,’ 7D5/2 ’,
+ ’ 7P1/2 ’,’ 7P3/2 ’,’ 6G7/2 ’,
+ ’ 6G9/2 ’,
+ ’ 6H9/2 ’,’ 6H11/2 ’,
+ ’ 6H9/2 ’,’ 6H11/2 ’,’ 6H11/2 ’,
+ ’ 7F5/2 ’,’ 7F7/2 ’,’ 7G7/2 ’,’ 7G9/2 ’,’
3S1/2’,
3D5/2 ’,
4D5/2 ’,
4F7/2 ’,
5F7/2
5G9/2
’,
’,
7H9/2
’/
We also need to set up the occupation numbers of the orbitals. For filled shells this is no problem but for outer
open shells a choice must be made. All of this is done in the subroutine setup which is part of the self-consistent
atomic program is listed in 9.78 . A table which associates the element name with the nuclear charge is also
included, so that information on the calulated eigenvalues looks like:
ATOM Fe
1S1/2
2S1/2
2P1/2
2P3/2
3S1/2
3P1/2
3P3/2
4S1/2
3D3/2
3D5/2
ATOMIC NUMBER=
26
eigenvalue=
-513.170714864160
eigenvalue=
-60.1251912919572
eigenvalue=
-51.9252187858152
eigenvalue=
-51.0102113290435
eigenvalue=
-6.90374797643669
eigenvalue=
-4.51245508603952
eigenvalue=
-4.39829228132764
eigenvalue= -0.434167420446594
eigenvalue= -0.609091963153407
eigenvalue= -0.597844758632259
occupation
occupation
occupation
occupation
occupation
occupation
occupation
occupation
occupation
occupation
=
=
=
=
=
=
=
=
=
=
2.00000000000000
2.00000000000000
2.00000000000000
4.00000000000000
2.00000000000000
2.00000000000000
4.00000000000000
2.00000000000000
2.40000000000000
3.60000000000000
Examples of the self-consistently calculated square of wave functions (or partial densities) for oxygen, iron and
uranium are shown in Figs. 4.2, 4.3 and 4.4. The partial densities for element number 136 are shown in Fig.4.5.
5
Calculated properties of free atoms
Since we have gone to the trouble of writing the program, it seems worthwhile to milk it for as many physical
properties as possible. Here are a few examples.
5.1
Trends in ionicity and atom size
The calculated outermost s and p energies of free atoms are shown in Figs. 5.1 and 5.2.
The outer electrons
8 Actually the program is the relativistic version but the SchrВЁ
odinger equation is recovered by setting the velocity of light to a
very large number.
5.1
Trends in ionicity and atom size
19
Figure 4.2: The partial radial electron densities for oxygen.
Figure 4.3: The partial radial electron densities for iron.
become less bound down a column of the periodic table, due to the addition of an orthogonality node for each row
down the column, which increases the kinetic energy. Both s and p energies sink across a row of the periodic table
with the s energies sinking more rapidly. These energy eigenvalues are closely related to the s and p ionisation
energies. The origin of this trend is the incomplete screening of electrons from the charge of the nucleus. As the
nuclear charge is increased by unity and an additional electron is added, the additional electron is not completely
screened by the other electrons and it is acted upon by an effective nuclear charge greater than one. The lack
of complete screening is greater for the lower angular momentum states, which penetrate closer to the nucleus,
and affects s electrons most. Consequently the s state energies sink more rapidly than p state energies.
Similar trends are discernable in the atomic sizes shown in Fig. 5.3. A good measure of effective atomic size is
the outermost node of the outermost valence electrons, which is related to the pseudo-potential. The radius of
outermost node is plotted in Fig. 5.3. The atoms increase in size down a column of the periodic table - again due
to the additional orthogonality node - and decrease in size across a row - due to incomplete screening. Contrary
to intuition, heavier atoms are not necessarily the largest because contraction across a row - due to incomplete
screening - is greater than expansion down a column - due to the additional orthogonality node. Thus the largest
atom in the periodic table is caesium.
20
5
CALCULATED PROPERTIES OF FREE ATOMS
Figure 4.4: The partial radial electron densities for uranium.
Figure 4.5: The partial radial electron densities for element 136.
A particularly interesting anomaly is that atoms following a transitions series9 , e.g. gallium, are smaller than
would be expected. This may be traced to the fact that the more penetrating p-states are not well screened by
the filled d-shell. Thus gallium is no larger than aluminium as shown in Fig. 5.4 where the radial densities of
the outermost wave functions of the two atoms are plotted and the outermost node is a measure of the effective
atom size.
Fig. 5.5 shows the s and d energy levels of transition metal atoms in the dnв€’1 s configuration. The d-level
sinks relative to the s-level across the row the effect again being due to the incomplete d в€’ d screening which is
responsible for the well known contraction of the size of the d-shell (for the f -shell of the rare earths this is known
as lanthanide contraction). Coulomb repulsion between the more localised 3d electrons tends to increase their
energies therefore the s в€’ d separation energy increase is greater for the 4d series. Curiously, the effect decreases
again for the 5d series for an entirely different reason. The heavier 5d atoms experience stronger relativistic effects
and the Darwin shift of the outer s-states increases their binding energies relative to the 5d-states, decreasing
the energy difference. Thus the ground state of the atoms Ni, Pd and Pt changes from 3d8 4s2 to 4d10 to 3d9 4s
9 For
a complete discussion, see Pettifor (1995).
5.2
Band masses of transition metals
21
Figure 5.1: Energies of the outermost s-states across the periodic table with transition metals, rare earths and
actinides omitted.
as the sp energy difference first increases and then decreases again.
5.2
Band masses of transition metals
If one knows the Wigner-Seitz radius of the atom in a solid, it is quite easy to estimate the band mass (O.K.Andersen,
H.L.Skriver, H.Nohl & B.Johansson 1979), m = 2/S 2 П†2 , from the self-consistent atom. The estimated bandwidths, obtained from the band masses are plotted in Fig.5.6.
5.3
Exchange integrals
The exchange integrals in the local density approximation may be used in Stoner theory. In terms of the exchange
integral matrices, the spin polarisation energy is, approximately(Brooks & Johansson 1983),
∆E =
ml Jll ml
ll
(5.1)
22
5
CALCULATED PROPERTIES OF FREE ATOMS
Figure 5.2: Energies of the outermost p-states across the periodic table with transition metals, rare earths and
actinides omitted.
The diagonal integrals are obtained from
Jll = в€’
2
3
r2 dr П†2l (r)П†2l (r)ars (r)/ (r)
(5.2)
where ars (r) is defined in the subroutine for the local density approximation. The diagonal exchange integrals
are calculated in the atom program upon request and the trends for transition metals are shown in Fig. 5.7.
5.3
Exchange integrals
23
Figure 5.3: Estimated atomic sizes for atoms with outermost s and p electrons. The radius of final node of the
outermost s-state is plotted.
24
5
CALCULATED PROPERTIES OF FREE ATOMS
Figure 5.4: The radial electron densities of Al, Ga and In. The gallium and Aluminium atoms are about the
same size.
5.3
Exchange integrals
25
Figure 5.5: The energies of s and d states of transition metal atoms.
Figure 5.6: Trends in the bandwidths of transition metals, rare earths and actinides, obtained from free atom
calculations.
26
5
CALCULATED PROPERTIES OF FREE ATOMS
Figure 5.7: Exchange integrals for the trnasition metal series calculated from the free atom densities.
27
6
Boundary conditions and series expansions for the SchrВЁ
odinger
equation.
In this section we look at the details of the radial ScrВЁodinger equation, the coupled radial equations, series
expansions at small radii and the WKB approximation at large radii.
6.1
ScrВЁ
odinger equation for a spherical potential
The ScrВЁ
odinger equation for a potential, V , is
2
в€’
2m
∇2 + V
ψ = εψ −→ ∇2 ψ = M (ε)ψ
where
M (Оµ) =
2m
2
(6.1)
(V в€’ Оµ)
(6.2)
is the motive energy. In one dimension ∇2 → d2 /dx2 and if M < 0 the solutions are oscillatory whereas if M > 0
the solutions are exponentials, with the classical turning points at M = 0. In a three dimensional spherical
potential
1 ∂2
L2
1 ∂2
L2
∂2
2 ∂
L2
∇2 = 2 2 r 2 − 2 2 =
r
в€’
=
+
в€’
(6.3)
2 r2
2 r2
r ∂r
r
r ∂r2
∂r2
r ∂r
the wave equation HП€ = ОµП€ becomes
1 ∂2
L2
r
в€’
+V П€
2m r ∂r2
2mr2
2
2
в€’
=
в€’
=
в€’
∂2
2 ∂
+
∂r2
r ∂r
2m
в€’
L2
+V П€
2mr2
p2r
L2
в€’
+ V П€ = ОµП€
2m 2mr2
(6.4)
where
pr = в€’i
∂
1
+
∂r
r
(6.5)
is the radial component of the momentum. Since L2 = L В· L = L2 , if the potential is spherical, the wave equation
is separable with solutions of the form
П€(r) = il Rl (r)Ylm (Л†
r)
(6.6)
whence
2
в€’
2
1 d2
l(l + 1)
r
+
+ V (r) в€’ Оµ Rl (r)
2m r dr2
2m r2
2
=
=
where k =
2
d2
l(l + 1)
+
+ V (r) в€’ Оµ rRl (r)
2m dr2
2m r2
d2
l(l + 1)
в€’
+ k 2 rRl (r) = 0
dr2
r2
в€’
(6.7)
2m(Оµ в€’ V ). Or
d2
rRl (r) = M (Оµ)rRl (r)
dr2
(6.8)
where
M (Оµ) =
2m
2
(V в€’ Оµ) +
l(l + 1)
r2
(6.9)
28
6
ВЁ
BOUNDARY CONDITIONS AND SERIES EXPANSIONS FOR THE SCHRODINGER
EQUATION.
Since the centrifugal potential is large and positive for small radii there are two classical turning points, M (Оµ) = 0,
for negative energies if the potential is attractive and Coulomb-like.
Another way to write Equation (6.7) is to use the dimensionless logarithmic derivative rR /R. Since
(rRl ) =
d
(Rl + rRl )
dr
(6.10)
Equation (6.7) becomes
=
d
r d
l(l + 1)
rRl =
(Rl + Rl Dl ) + k 2 в€’
(Rl + Rl Dl ) в€’ l(l + 1) + k 2 r2
2
dr
r
Rl dr
dDl
dDl
Dl + Dl2 + r
в€’ l(l + 1) + k 2 r2 (Dl в€’ l)(Dl + l + 1) + r
+ k2 r2 = 0
dr
dr
Hence
в€’r
dDl
= (Dl в€’ l)(Dl + l + 1) + k 2 r2
dr
(6.11)
In any region where the potential is flat Equation (6.7) reduces to the radial Helmholtz equation since k is
constant and, if the energy is also equal to the potential energy, to the radial Laplace equation since then k
is zero. The solutions to the single electron wave equation for a spherical ion core embedded in a region of
flat potential are obtained by joining the solutions to the Helmholtz equation for the interstitial region to the
solutions within the ion cores continuously and differentiably.
6.2
Coupled first order radial equations - series expansion
The ScrВЇ
odinger equation in a spherical symmetry, Equation (6.7), may be written, with P = rRl
where
P = ВµP
(6.12)
l(l + 1)
в€’ k2
r2
(6.13)
Q=P в€’
(l + 1)
P
r
(6.14)
P =Q+
(l + 1)
P
r
10
Вµ=
If one defines
or
в€љ
в€љ
development is as follows. If one writes P = rY (Y = rR) and lets r = exp(x), one finds d/dr = exp(в€’x)d/dx. If
Y = dY /dx one may evaluate dP/dr and d2 P/dr2 in terms of Y and Y . One finds
„
В«
dP
Y
= eв€’x/2
+Y
dr
2
10 Another
and
„
В«
d2 P
1
в€’3x/2
=
e
в€’
Y
+
Y
dr2
4
and again there are no terms in Y . Hence Equation (6.12) becomes
Y
= ОіY
where
1 2
)
2
For the first term to be negligible for l = 0 at small radii, we require that r0 = ex0
Оі = (V в€’ Оµ)e2x + (l +
1/8Z, since V в€ј в€’2Z/r at small radii.
6.2
Coupled first order radial equations - series expansion
29
then
Q
(l + 1)
(l + 1)
(l + 1)
(l + 1)
Pв€’
Pв€’
P = Вµ+
P
r2
r
r2
r
(l + 1)
(l + 1) (l + 1)2
в€’
Pв€’
Вµ+
Q
r2
r2
r
(l + 1)
l(l + 1)
(l + 1)
Pв€’
Вµв€’
Q = в€’k 2 P в€’
Q
r2
r
r
=
P +
=
=
Hence the coupled first order differential equations, with11 k 2
P
Q
2
= 2m(Оµ в€’ V )
(l + 1)
P
r
(l + 1)
= в€’k 2 P в€’
Q
r
= Q+
(6.15)
(6.16)
are equivalent to the SchrВЇ
odinger equation.
Equations (6.15) and (6.16) become (differentiation with respect to x)
P
Q
= Qex + (l + 1)P
2
(6.17)
x
= в€’k P e в€’ (l + 1)Q
(6.18)
In atomic (Rydberg) units k 2 = Оµ в€’ V . As r в†’ 0 x в†’ в€’в€ћ and V = в€’2Z/r = в€’2Z exp(в€’x), k 2 r в†’ 2Z hence
the limits of Equations (6.17) and (6.18) are
lim P
=
(l + 1)P
(6.19)
lim Q
=
в€’2ZP в€’ (l + 1)Q
(6.20)
r→0
r→0
which have solutions Q = A exp(ax), P = B exp(ax) with ОІ = A/B = Q/P = в€’Z/(l + 1) and a = l + 1. The
power series expansion is
Q = e(l+1)x Q0 + Q1 ex + Q2 e2x + В· В· В·
Q
= e(l+1)x (l + 1)Q0 + (l + 2)Q1 ex + (l + 3)Q2 e2x + В· В· В·
P
= e(l+1)x P0 + P1 ex + P2 e2x + В· В· В·
P
= e(l+1)x (l + 1)P0 + (l + 2)P1 ex + (l + 3)P2 e2x + В· В· В·
(6.21)
Substitution in Equation (6.17) and equating coefficients of exp(nx) yields
P1 = Q0
2P2 = Q1
3P3 = Q2
(6.22)
Similarly, Substitution in Equation (6.18) and equating coefficients of exp(nx) yields, with V = в€’2Z/r =
−2Z exp(−x) + ∆V
2(l + 1)Q0
=
в€’2ZP0
(l + 2)Q1 = −2ZP1 + (∆V − ε)P0 − (l + 1)Q1
(l + 3)Q2
=
−2ZP2 + (∆V − ε)P1 − (l + 1)Q2
An P0 is determined by normalisation.
11 In
atomic (Rydberg) units,
2 /2m
= 1 and k2 = (Оµ в€’ V ).
(6.23)
30
6
6.3
ВЁ
BOUNDARY CONDITIONS AND SERIES EXPANSIONS FOR THE SCHRODINGER
EQUATION.
Boundary conditions at large radius - the WKB approximation
We need the boundary condition for a radius, far beyond the classical turning point, to start the inward integration. We use a version of the WKB approximation. Consider the second order linear equation
p (x) + О±(x)p (x) + ОІ(x)p(x) = 0
(6.24)
which becomes, if one makes the substitution p(x) = v(x)w(x)
v + 2
w
w + О±w + ОІw
+О± v +
v=0
w
w
(6.25)
Therefore, if w is chosen such that 2w /w + О± = 0, the first derivative vanishes and, integrating,
w(x) = eв€’(1/2)
Then
1
w = в€’ О±w
2
and Equation (6.25) becomes
where
R
dx О±(x)
(6.26)
1
1
1
1
w = в€’ О± w в€’ О±w = в€’ О± w + О±2 w
2
2
2
4
v + Вµv = 0
(6.27)
1
1
Вµ(x) = ОІ в€’ О± в€’ О±2
2
4
(6.28)
and the solution to Equation (6.24) is
p(x) = v(x)eв€’(1/2)
R
dx О±(x)
(6.29)
If v(x) is slowly varying Equation (6.27) may be solved using the WKB approximation. If Вµ(x) were independent
в€љ
of x, the solutions to Equation (6.27) would be proportional to eiφ with φ2 = µ, φ = ± µ. If µ < 0, φ is
imaginary and the solutions to Equation (6.27) are exponentials whereas if Вµ > 0, П† is real and the solutions to
Equation (6.27) are oscillatory. One might therefore try
v = eiφ(x)
(6.30)
iφ − (φ )2 + µ = 0
(6.31)
in which case Equation (6.27) becomes
If П† is small the first approximation to Equation (6.31) is
в€љ
П† =В± Вµ
в†’ П†(x) = eВ±
in which case
П† =
R
dx
в€љ
Вµ(x)
(6.32)
1 Вµ
в€љ
2 Вµ
and, iterating Equation (6.31) by substituting this result for П† , one finds
(φ )2 ≈ µ ±
i Вµ
в€љ
2 Вµ
iВµ
в€љ
→φ ≈± µ+
4Вµ
(6.33)
Consequently
П†(x) = В±
dx
Вµ(x) +
i
ln Вµ(x)
4
(6.34)
6.3
Boundary conditions at large radius - the WKB approximation
Then
v(x) = eВ±i
R
dx
в€љ
Вµ(x)в€’ 14 ln Вµ
31
= Вµв€’1/4 eВ±i
R
dx
в€љ
Вµ(x)
The general solutions to Equation (6.27) is therefore, for Вµ > 0
в€љ
в€љ
R
R
v(x) = Вµв€’1/4 aei dx Вµ(x) + beв€’i dx Вµ(x)
(6.35)
(6.36)
and in a classically forbidden region, Вµ < 0,
v(x) = (в€’Вµ)в€’1/4 aeв€’
R
dx
в€љ
в€’Вµ(x)
R
+ be
dx
в€љ
в€’Вµ(x)
(6.37)
if the exponentially increasing solution is unphysical
v(x) = a(в€’Вµ)в€’1/4 eв€’
and
dv
=
dx
R
dx
в€љ
1Вµ
в€’ в€’Вµ +
4Вµ
в€љ
в€’Вµ(x)
v
(6.38)
(6.39)
From Equation (6.29)
dp
dx
=
=
=
R
dv в€’(1/2) R dx О± 1
e
в€’ О±veв€’(1/2) dx О±
dx
2
R
R
в€љ
1
1Вµ
veв€’(1/2) dx О± в€’ О±veв€’(1/2) dx О±
в€’ в€’Вµ +
4Вµ
2
в€љ
1
1Вµ
в€’ О± p
в€’ в€’Вµ +
4Вµ
2
(6.40)
32
7
7.1
7
RELATIVISTIC EFFECTS
Relativistic effects
The Dirac equation
Since the Dirac Hamiltonian is
H = cО± В· p + ОІmc2
(7.1)
The Dirac equation for the 4-component state, П€, is therefore
i
∂ψ
= HП€ = cО± В· p + ОІmc2 П€
∂t
(7.2)
In the limit v/c в†’ 0, the Dirac equation may be approximated by the Pauli equation, a scalar wave equation similar to the SchrВЇ
odinger equation with added spin degrees of freedom. Since the spin, s, is an angular
momentum it obeys the commutation rule s Г— s = i s. The magnetic moment associated with the spin is
Вµs = в€’(e /mc)s = в€’2Вµ0 s.
The operator
K = 1 + Пѓ В· l/
(7.3)
commutes with the spin-orbit interaction and is a constant of the motion with eigenvalues в€’Оє. For the two
possibilities j = l + 1/2 and j = l в€’ 1/2 one finds
1
2
1
j =lв€’
2
j =l+
K
= l+1
K
= в€’l
Оє = в€’l в€’ 1
Оє=l
(7.4)
The values of Оє for j = l+1/2 states are в€’1, в€’2, в€’3, . . ..For j = lв€’1/2 there are no l = 0 states and Оє = 1, 2, 3, . . ..
Therefore Оє takes the values
Оє = В±1, В±2, В· В· В·
(7.5)
The eigenstates of K are the eigenstates of j 2 with
j = |Оє| в€’
1
2
(7.6)
It is therefore common usage to write the spin angular functions as
П‡ВµОє =
ms
1
lm ms |jВµ YlВµв€’ms (Л†
r )П‡ms
2
(7.7)
with Оє replacing the quantum numbers j and l. It is conventional to also use the quantum numbers ВЇl and the
sign of Оє, SОє , defined by,
1
ВЇl = в€’Оє = l + 1
2
1
ВЇl = Оє в€’ 1 = l в€’ 1
j =lв€’
2
SОє = l в€’ ВЇl = 2(l в€’ j)
j =l+
(7.8a)
(7.8b)
(7.8c)
The parity of the eigenfunctions is determined by Ylm and is therefore
ПЂl = (в€’1)l = (в€’1)j+SОє /2
(7.9)
7.2
Coupled first order radial Dirac equations in a central field
33
Although it may seem that there are too many quantum numbers, it actually requires both l and j to label a
state uniquely whereas labelling by Оє is unique. For example, for principle quantum number equal to 2, the s1/2
and p1/2 both have j = 1/2, as does the j = 1/2 state with principle quantum number equal to 1. However Оє
changes sign, negative for j = l + 1/2 states and positive for j = l в€’ 1/2 states. The maximum value of l for
principle quantum number equal to n is n в€’ 1. It follows that Оє ranges from в€’1 в†’ в€’n for j = l + 1/2 states and
from 1 в†’ n в€’ 1 for j = l в€’ 1/2 states. For the 1s1/2 , Оє = в€’1, for the 2s1/2 , Оє = 1 and for the 3p1/2 , Оє = 1. The
quantum number Оє therefore contains information both upon the value of j and whether j = l В± 1/2.
The operator Пѓ r , defined by
rПѓ r = ПѓВ·r
(7.10)
(ПѓВ·r) = (r В· r) + iПѓВ·(r Г— r) = r2
(7.11)
has odd parity, since r has odd parity and
2
hence Пѓ 2r = 1. K and Пѓ r anti-commute. It follows from Equation (7.3) that
Пѓ r П‡Вµв€’Оє = в€’П‡ВµОє
Пѓ r П‡ВµОє = в€’П‡Вµв€’Оє
(Пѓ r )2 П‡ВµОє = П‡ВµОє
(7.12)
The operator σ r therefore converts the j = l ± 1/2 state to the j = l ∓ 1/2 state whilst keeping the value of j
constant by changing the value of l. If one starts with a j = l + 1/2 state, with angular momentum l, Оє = в€’l в€’ 1,
the effect of operating with Пѓ r is Оє в†’ в€’Оє or Оє = l + 1. Now Оє > 0 corresponding to a j = l в€’ 1/2 where
Оє = l = l + 1. The new angular momentum is therefore l = l + 1 = ВЇl = 2j в€’ l, from Equation (7.8). If one starts
with a j = l в€’ 1/2 state, with angular momentum l, Оє = l, the effect of operating with Пѓ r is again Оє в†’ в€’Оє or
Оє = в€’l. Now Оє < 0 corresponding to a j = l + 1/2 with Оє = в€’l for which the angular momentum is obtained
from Оє = в€’l в€’1. The new angular momentum is therefore l = l в€’1 = ВЇl = 2j в€’1, from Equation (7.8). Therefore
Пѓ r changes j = l + 1/2 states into j = l в€’ 1/2 states with l в†’ l + 1 and changes j = l в€’ 1/2 states into j = l + 1/2
states with l в†’ l в€’ 1, conserving j.
7.2
Coupled first order radial Dirac equations in a central field
The Dirac Hamiltonian for a spherical potential is
H = cО± В· p + ОІmc2 + V (r) = cОі 5 ОЈ В· p + ОІmc2 + V (r)
(7.13)
where
Оі5 =
0
1
1
0
and
ОЈ=
Пѓ
0
0
Пѓ
(7.14)
where ОЈ В· p is given by
ОЈВ·p=i
ОЈВ·r
∂
1
в€’
+ ОЈВ·l
r
∂r
r
(7.15)
The Dirac Hamiltonian for a spherical potential becomes
H = icОі 5
ОЈВ·r
r
в€’
∂
ОЈВ·l
+
∂r
r
+ ОІmc2 + V (r) = в€’icОі 5 ОЈr
∂
ОЈВ·l
в€’
∂r
r
+ ОІmc2 + V (r)
(7.16)
The four component Dirac wave functions are divided in two component upper and lower parts
П€=
П€u
П€l
(7.17)
34
7
RELATIVISTIC EFFECTS
where П€u and П€l are two-component spinors. The Dirac equation becomes
пЈ®
пЈ№пЈ® пЈ№
∂
1
П€
M1
в€’icПѓ r
в€’ Пѓ В· l пЈє пЈЇ uпЈє
пЈЇ
∂r r
пЈЇ
пЈєпЈЇ пЈє = 0
пЈ°
пЈ»пЈ° пЈ»
∂
1
в€’icПѓ r
в€’ ПѓВ·l
M2
П€l
∂r r
(7.18)
where
M1
= mc2 + V в€’ W
M2
= в€’mc2 + V в€’ W
(7.19)
In the non-relativistic limit, when only the large components are significant, the Pauli equation is recovered and
[l, H] = 0. One large component must have spin up and the other spin down and j 2 and jz must be constants
of the motion. The Dirac wave function
gОє П‡ОєВµ
О¦=
(7.20)
в€’ifОє Пѓr П‡ОєВµ
obeys the coupled first order differential equations
gОє
=
cfОє
=
Оє+1
gОє + 2M cfОє
r
Оєв€’1
cfОє + (V в€’ E)gОє
r
в€’
(7.21)
(7.22)
With the rest mass subtracted from the energy, Equation (7.19) becomes
M1 =mc2 + V в€’ W = V в€’ E
(7.23a)
M2 = в€’ mc2 + V в€’ W = в€’2mc2 + V в€’ E = в€’2mc2 1 +
where
M =m 1+
Eв€’V
2mc2
Eв€’V
2mc2
= в€’2M c2
(7.23b)
(7.24)
With P = rg and Q = rcf , Equation (7.21) becomes
QОє
Оє
= в€’ PОє + 2M QОє
r
Оє
=
QОє + (V в€’ E)PОє
r
PОє
в†’
QОє
в†’
PОє
(7.25)
(7.26)
Note that for c в†’ в€ћ 2M в†’ 1 and
Оє
в€’ PОє + QОє
r
Оє
QОє + (V в€’ E)PОє
r
(7.27)
(7.28)
which are equivalent to the coupled SchrВЇ
odinger equations if Оє = в€’l в€’ 1.
On a radial mesh r в€ќ exp x, Equations (7.25) and (7.26) become - with the prime now denoting differentiation
with respect to x
PОє
QОє
= в€’ОєPОє + 2M QОє ex
= ОєQОє + (V в€’ E)PОє e
(7.29)
x
(7.30)
7.3
The radial Pauli equation
35
For r в†’ 0, x в†’ в€’в€ћ and 2M в†’ 2Zeв€’x /c2 since V в†’ в€’2Z/r in Rydbergs
2Z
QОє
c2
lim QОє в†’ ОєQОє в€’ 2ZPОє ex
lim PОє в†’ в€’ОєPОє +
(7.31)
r→0
(7.32)
r→0
For solutions of the form P в€ј ArОі в€ј A exp(Оіx) and Q в€ј BrОі в€ј B exp(Оіx) one finds
(Оі + Оє)A в€’
2Z
B=0
c2
and
2ZA + (Оі в€’ Оє)B = 0
therefore
Оі 2 в€’ Оє2 +
2Z
c
2
=0
or
Hence
lim Оі = |Оє| в€’
c→∞
1
2|Оє|
Оє2 в€’
Оі=
2Z
c
2Z
c
2
(7.33)
2
(7.34)
The ratio Q/P for is
c2
2Z
B
QОє
в†’
(Оі + Оє) =
=c
=
r→0 Pκ
A
2Z
Оєв€’Оі
ОІ = lim
7.3
Оі+Оє
Оєв€’Оі
(7.35)
The radial Pauli equation
We now derive a second order radial equation for the major component from
Оє+1
gОє + 2M cfОє
r
Оєв€’1
cfОє (r) =
cfОє + (V в€’ E)gОє
r
gОє (r) = в€’
(7.36a)
(7.36b)
where, M is given by Equation (7.24). Equation (7.36a) may be solved for cfОє and, if it is differentiated with
respect to r, solved for fОє
cfОє (r) = {gОє (r) + (Оє + 1)gОє (r)/r} /2M
(7.37a)
2
cfОє (r) = gОє (r) + (Оє + 1)gОє (r)/r в€’ (Оє + 1)gОє (r)/r в€’ 2cfОє M
/2M
(7.37b)
Then, if Equations (7.37) are substituted into Equation (7.36b), one finds the following second order differential
equation for gОє
gОє (r) + 2
gОє (r) Оє(Оє + 1)
M
в€’
gОє (r) в€’
r
r2
M
gОє (r) +
Оє+1
g + 2M (E в€’ V )gОє (r) = 0
r
(7.38)
Equation (7.38) is equivalent to the Dirac radial coupled first order equations. Since Оє(Оє + 1) = l(l + 1), the
second term will be recognised to be the centrifugal potential and only one term, spin-orbit coupling, depends
upon Оє. The mass velocity and Darwin terms depend upon M and M . If the spin-orbit interaction is removed
by setting Оє + 1 = 0, we obtain the relativistic Pauli equation includes mass velocity and Darwin terms and
describes an electron with almost pure spin up and spin down states but with all other relativistic included, The
Оє dependence in Equations (7.36) then becomes a hindrance. It is therefore useful to replace the equation for
f with another for another function obeying a differential equation independent of Оє, for example one might
define, similar to Equation (7.36a),
1
g
(7.39)
П†=
2M c
36
7
RELATIVISTIC EFFECTS
Then it follows from Equation (7.38) that
2
l(l + 1) V в€’ E
П† =в€’ П†+
+
g
r
2M cr2
c
(7.40)
which has no Оє-dependence. Equations (7.39) and (7.40) are the new coupled radial equations, independent of
Оє where, from Equation (7.36b),
1 Оє+1
g
(7.41)
f =П†+
2M c r
An alternative to Equation (7.38) is obtained if one makes the substitution P = rg since then
g = в€’P/r2 + p /r
g = 2P/r3 в€’ 2P /r2 + P /r
g + 2g /r = P /r
and one finds
Pl + О±Pl = ВµPl
(7.42)
where Pl must obey the Dirac equation (7.25) and Оє(Оє + 1) = l(l + 1), Оє = в€’(1 + Пѓ В· l) and
О±
=
Вµ =
M
M
l(l + 1) M Оє
+
+ 2M (V в€’ E)
r2
M r
в€’
(7.43)
(7.44)
A comparison of Equations (7.42) with Equations (6.12) and a glance at Equation (6.14) suggests that, yet
another set of coupled Оє-independent first order radial equations may be obtained if we define QОє by
2M Ql = Pl в€’
l+1
Pl
r
(7.45)
If one sets П† = Q/cr, and uses the coupled Dirac equations one finds that
П†l =
1
1 l+1
1
1 l
(g + rg ) в€’
rg =
g в€’
g
2M cr
2M cr r
2M c
2M c r
(7.46)
Again, from the coupled Dirac equations
f=
1
1 Оє+1
1 Оє+l+1
g +
g =П†+
g
2M c
2M c r
2M c
r
(7.47)
Therefore the Dirac wave function, Equation (7.20), becomes, in terms of П†
О¦ОєВµ =
в€’i П†l +
gl П‡ОєВµ
1 Оє+1+1
gl
2M c
r
Пѓr П‡ОєВµ
(7.48)
With the radial functions now depending only upon l, linear combinations of the spin functions may be taken to
produce spin up and down states for the major component (the minor component still contains a spin mixture)
О¦lms =
в€’i П†l +
gl Ylm П‡s
Пѓ В· l+1 g
1
2M c
r
l
Пѓr Ylm П‡s
=
gl Ylm П‡s
1
Пѓ В· lg Пѓ Y П‡
в€’g
+
i 2M
l
r lm s
l
c
r
(7.49)
ortho-normalised according to
О¦l m s |О¦lms = Оґll ;mm ss Nl в€’ Оґll Sl Yl m П‡s |Пѓ В· l|Yl m П‡s
(7.50)
7.3
The radial Pauli equation
37
where
r2 dr gl2 +
Nl =
Sl =
1
4M 2 c2
1
4M 2 c2
r2 dr
gl2 +
2gl gl +
l(l + 1) 2
gl
r2
(7.51a)
gl2
r2
(7.51b)
The terms in 1/c2 may be neglected for valence states but not for core states.
The new coupled first order equations may be derived by first differentiating Equation (7.45), leading to
2M Ql + 2M Ql
l+1
l+1
Pl в€’
Pl
r2
r
l+1
l+1
Вµ+ 2
Pl в€’ О± +
r
r
l+1
l+1
Pl в€’ О± +
Вµ+ 2
r
r
= PОє +
=
=
Pl
2M Ql +
l+1
Pl
r
(7.52)
or, substituting from Equations (7.43) and (7.44)
2M
l+1
Ql + 2M Ql
r
=
l(l + 1)
l+1
в€’О±
Pl
2
r
r
M l+1
M Оє
2M (V в€’ E) +
Pl
M r
M r
Вµв€’
=
(7.53)
Hence we have the coupled first order differential equations
l+1
Pl
r
l+1
M Оє+l+1
+ (V в€’ E) Pl в€’
Ql =
Ql
2
2M
r
r
Pl =2M Ql +
(7.54a)
(7.54b)
Noting that Оє + l + 1 = l в€’ Пѓ В· l and one finds
l+1
Pl
r
M Оє+l+1
l+1
Ql =
+ (V в€’ E) Pl в€’
Ql
2
2M
r
r
Pl =2M Ql +
(7.55a)
(7.55b)
With x = ln r in the limit r в†’ 0
2M в†’ 2
Z в€’x
e
c2
M
в†’ в€’eв€’x
M
Z в€’2x
e
c2
M
c2
в†’
в€’
2M 2
2Z
2M в†’ в€’2
(7.56)
Then
Pl
=
Ql
=
2Z
l+1
Ql +
Pl
2
c
r
c2 Оє + l + 1 2Z
l+1
в€’
в€’
Pl в€’
Ql
2Z
r
r
r
(7.57)
(7.58)
38
7
RELATIVISTIC EFFECTS
where differentiation is with respect to r, or changing the variable to x = ln r
Pl
=
Ql
=
2Z
Q + (l + 1)Pl
c2 l
c2
в€’ (Оє + l + 1) в€’ 2Z Pl в€’ (l + 1)Ql
2Z
(7.59)
(7.60)
where differentiation is with respect to x. Now with solutions of the form P = p exp(Оіx) and Q = q exp(Оіx)
substitution in Equations (7.59) and (7.60) yields
c2
(Оє + l + 1) в€’ 2Z) p
2Z
(l + Оі + 1)q
=
в€’
(Оі в€’ l в€’ 1)p
=
2Z
q
c2
which may be solved for Оі
2Z
c2
Оі 2 = l(l + 1) в€’ Оє в€’
and for
2
(7.61)
(Оі в€’ l в€’ 1)
q
=
p
2Z/c2
(7.62)
Equation (7.38) may be rearranged and expressed in terms of logarithmic derivatives using
2
1 d dgОє
gОє
gОє + gОє = 2 r2
= 2 DОє + DОє2 + rDОє
r
r dr dr
r
in which case it becomes
в€’rDОє = (DОє в€’ 1)(DОє + l + 1) + 2M (E в€’ V )r2 +
1
V (DОє + Оє + 1)
2M c2
(7.63)
where M = в€’V /2c2 and the final term is due to spin-orbit coupling.
7.4
Generic radial SchrВЁ
odinger, relativistic Pauli and Dirac equations
SchrВЁ
odinger gave us a second order differential equation which we turned into to coupled first order equations.
Dirac gave is a fourth order differential equation, which we separated into two sets of coupled first order equations,
and then we derived a second order differential equation from the coupled first order equations to obtain the
relativistic Pauli equation. It is not surprising, therefore, that all three equations may be written in the generic
form
dP
=a11 P + a12 Q + a13 Qr
dr
dQ
=a21 P + a22 Q + a23 P r
r
dr
r
(7.64a)
(7.64b)
where, for the ScrВЁ
odinger equation
a11
=
l+1
a22
=
в€’(l + 1)
a12 = 0
a21 = 0
a13 = 1
a23 = (V в€’ E)
(7.65)
7.4
Generic radial SchrВЁ
odinger, relativistic Pauli and Dirac equations
39
and for the Pauli Equation
a11
=
a22
=
l+1
a12 = 2M r
a13 = 0
M
в€’(l + 1)
a21 =
(Оє + l + 1)
2M 2
a23 = (V в€’ E)
(7.66)
and for the Dirac equation
a11
=
в€’Оє
a22
=
Оє
a12 = в€’2SОє M r
a13 = 0
a23 = в€’SОє (V в€’ E)
a21 = 0
(7.67)
At the starting point for numerical outward integration, r в†’ 0, (in atomic Rydberg units)
V (r)r в†’ в€’2Z
2M r в†’
2Z
c2
M
c2
в†’
в€’
2M 2
2Z
(7.68)
therefore all of the amn are independent of r.
Equations (7.64) become, if P в†’ rОі p0 , Q в†’ rОі q0 ,
lim ОіP =a11 P + a12 Q
в†’
(Оі в€’ a11 )p0 = a12 q0
(7.69a)
lim ОіQ =a22 Q + a21 P
в†’
(Оі в€’ a22 )q0 = a21 p0
(7.69b)
r→0
r→0
Consequently
(Оі в€’ a11 )(Оі в€’ a22 ) = a12 a21
Оі в€’ a11
a21
q0
=
=
p0
a12
Оі в€’ a22
(7.70a)
(7.70b)
With p0 = 1, to be set later by normalisation, Оі and q0 may be obtained from Equations (7.70). A power series
expansion may then be started with
lim P (r) =rОі (p0 + p1 r + p2 r2 + p3 r3 + . . .)
(7.71a)
lim Q(r) =rОі (q0 + q1 r + q2 r2 + q3 r3 + . . .)
(7.71b)
r→0
r→0
Then Equations (7.64) become
rОів€’1 (Оіp0 + (Оі + 1)p1 r + . . .) =a11 rОів€’1 (p0 + p1 r + . . .) + a12 rОів€’1 (q0 + q1 r + . . .)
+a13 rОів€’1 (q0 r + q1 r2 + . . .)
r
Оів€’1
(Оіq0 + (Оі + 1)q1 r + . . .) =a22 r
Оів€’1
(q0 + q1 r + . . .) + a21 r
+a23 r
Оів€’1
2
(7.72a)
Оів€’1
(p0 + p1 r + . . .)
(p0 r + p1 r + . . .)
(7.72b)
(Оі + n в€’ a11 )pn в€’ a12 qn =a13 qnв€’1
(7.73a)
(Оі + n в€’ a22 )qn в€’ a21 pn =a23 pnв€’1
(7.73b)
Equating coefficients of rn Г— rОів€’1 one finds
hence
в€’a21
a + n в€’ a11
Оі + n в€’ a22
в€’a12
pn
qn
=
a23 pnв€’1
a13 qnв€’1
(7.74)
which may be solved for pn and qn
pn
qn
=
в€’a21
Оі + n в€’ a11
Оі + n в€’ a22
в€’a12
в€’1
в€’a12
в€’(Оі + n в€’ a22 )
в€’(Оі + n в€’ a11 )
в€’a21
a23 pnв€’1
a13 qnв€’1
Starting with p0 and q0 from the limit r в†’ 0, the wave functions may be iterated outwards.
(7.75)
40
7.5
7
RELATIVISTIC EFFECTS
The spin polarised Dirac equation
When an effective field along the z-axis acts on the spin an additional term Bz Пѓ z must be added to the Dirac
Hamiltonian. The diagonal matrix elements between the major components, for which the angular functions are
Ωjlµ are 2µ/(2l + 1) between the j = l + 1/2 states, for which spin and orbit are couple parallel, and −2µ/(2l + 1)
between the j = l в€’ 1/2 states for which spin and orbit are couple antiparallel. The minor components have
angular functions Ωj,2j−l,µ and, since 2j − l = l + 1 for j = l + 1/2 states, the diagonal components are,
2Вµ/(2l + 3). For the j = l в€’ 1/2 states 2j в€’ l = l в€’ 1 the matrix elements are в€’2Вµ/(2l в€’ 1). In addition, there
are matrix elements of Пѓ z between the j = l + 1/2 and j = l в€’ 1/2 states of the major components, equal to
1/2
GВµj j = в€’ (l + 1/2)2 в€’ Вµ2
/(l + 1/2), which couple the two pairs of Dirac equations. The matrix elements
of Пѓ z between the minor components are zero because the angular momenta are not equal and the spherical
harmonics are orthogonal. The Dirac coupled equations are therefore of order four
пЈ№
пЈ®
пЈ® пЈ№
∂
l+1
2Вµ
Вµ
в€’c
+
Bz Gj j
0
пЈє Pk
пЈЇV в€’ E + 2l + 1 Bz
∂r
r
пЈєпЈЇ пЈє
пЈЇ
пЈєпЈЇ пЈє
пЈЇ
2Вµ
2
пЈє пЈЇQk пЈє
 c ∂ − l+1
V в€’ E в€’ 2mc +
Bz
0
0
пЈєпЈЇ пЈє
пЈЇ
∂r
r
2l
+
3
пЈєпЈЇ пЈє = 0
пЈЇ
пЈєпЈЇ пЈє
пЈЇ
∂
l
2Вµ
Вµ
пЈє пЈЇRk пЈє
пЈЇ
Bz
c
в€’
Bz Gj j
0
V в€’E в€’
пЈєпЈЇ пЈє
пЈЇ
2l + 1
∂r r
пЈєпЈ° пЈ»
пЈЇ
пЈ»
пЈ°
∂
1
2Вµ
2
Sk
0
0
в€’c
+
V в€’ E в€’ 2mc в€’
Bz
∂r
r
2l в€’ 1
(7.76)
where
1/2
(l + 1/2)2 в€’ Вµ2
GВµj j = в€’
l + 1/2
where
пЈ№ пЈ®
пЈ№ пЈ®
пЈ№
rgk+ (r)Ωjlµ
rgk+ (r)Ωjlµ
Pk
Qk  rf + (r)Ωj,2j−l,µ  −irf + (r)σ r Ωjlµ 
k
пЈЇ пЈє=пЈЇ k в€’
пЈє пЈЇ
пЈє
Rk   rg (r)Ωjlµ  =  rg − (r)Ωjlµ 
k
k
Sk
rfk− (r)Ωj,2j−l,µ
−irfk− (r)σ r Ωjlµ
пЈ®
(7.77)
When the Dirac equation is spin polarised the 2j + 1 fold degeneracies are broken and all eigenstates are singlets.
For very large spin, the eigenstates are approximately spin up ans down states with the 2l + 1 degeneracy broken
by spin orbit coupling as illustrated in Fig. 7.1. When the shell is filled for 14 electrons, there can be no moment
and the states of a given j are degenerate. For a half-filled shell, a large spin can develop and states with
predominant spin up lie lower in energy than those with predominant spin down.
7.5
The spin polarised Dirac equation
41
Figure 7.1: Eigenvalues of an f-shell as a function of electron number, where spin polarisation and spin-orbit
coupling compete.
42
8
8
VIRIAL THEOREM
Virial theorem
8.1
Virial theorem from the SchrВЁ
odinger equation
We first of all operate with r · ∇ on the Schr¨odinger equation ∇2 ψ = (V − E)ψ and premultiply by ψ ∗ . One
finds
ψ ∗ (r · ∇)∇2 ψ = ψ ∗ (r · ∇)(V − E)ψ = ψ ∗ ψ(r · ∇)V + ψ ∗ (V − E)(r · ∇)ψ
But ∇2 ψ ∗ = (V − E)ψ ∗ therefore
ψ ∗ (r · ∇)∇2 ψ = ψ ∗ ψ(r · ∇)V + (∇2 ψ ∗ )(r · ∇)ψ
or
ψ ∗ (r · ∇)∇2 ψ − (∇2 ψ ∗ )(r · ∇)ψ = |ψ|2 (r · ∇)V
(8.1)
In components
пЈ®
ψ ∗ (r · ∇)∇2 ψ − (∇2 ψ ∗ )(r · ∇)ψ =
пЈ°П€ в€—
i
j
∂ 2 ψ∗
∂3ψ
в€’
xj
2
∂xj ∂xi
∂x2i
пЈ№
∂ψ 
xj
∂xj
j
(8.2)
Now consider the vector
which may be transformed to
G = ψ ∗ ∇ψ + ψ ∗ (r · ∇)∇ψ − (∇ψ ∗ )(r · ∇)ψ
(8.3)
G = ψ ∗ ∇(r · ∇ψ) − (∇ψ ∗ )(r · ∇)ψ
(8.4)
12
From Equation (8.3), the components of G are
Gi = (П€ в€— )2
Then
∂
∂xi
j
xj ∂ψ∂xj
П€в€—
= П€в€—
∂Gi
∂2ψ
= 2П€ в€— 2 + П€ в€—
∂xi
∂xi
and
i
∂2ψ
∂Gi
в€’ 2П€ в€— 2 =
∂xi
∂xi
∂ψ
+ П€в€—
∂xi
xj
j
пЈ°П€ в€—
j
which is just the RHS of equation (8.2). Hence, since
ψ ∗ (r · ∇)∇2 ψ − (∇2 ψ ∗ )(r · ∇)ψ =
i
j
∂2ψ
∂ψ ∗
в€’
∂xi ∂xj
∂xi
∂ 2 ψ∗
∂3ψ
в€’
2
∂xj ∂xi
∂x2i
пЈ®
i
xj
i
xj
j
∂3ψ
∂ 2 ψ∗
xj 2
в€’
∂xi ∂xj
∂x2i
xj
j
∂ψ
∂xj
∂ψ
∂xj
j
пЈ№
∂ψ 
xj
∂xj
(8.5)
(8.6)
(8.7)
∂Gi /∂xi = ∇ · G,
∂Gi
в€’ 2П€ в€—
∂xi
i
∂2ψ
= ∇ · G + 2ψ ∗ (−∇2 )ψ
∂x2i
(8.8)
and, from Equation (8.1),
2ψ ∗ (−∇2 )ψ + ∇ · G = |ψ|2 (r · ∇)V
12 The
vector identity
ψ ∗ ∇(r · ∇ψ) = ψ ∗ ∇ψ + ψ ∗ (r · ∇)∇ψ
is used.
(8.9)
8.1
Virial theorem from the SchrВЁ
odinger equation
43
If one integrates over space and the wave functions vanish at the boundaries of integration13 , the integral over
the second term vanishes and
2 dr ψ ∗ (−∇2 )ψ = dr |ψ|2 (r · ∇)V
(8.10)
or, if the eigenfunctions and eigenvalues are labelled by i,
dr 2ψi∗ (−∇2 )ψi =
2ti =
dr |ψi |2 (r · ∇)V
(8.11)
where ti is the kinetic energy of the eigenstate labelled by i. Twice the total kinetic energy is therefore, if no
orbitals have a surface contribution,
ni ti =
2T = 2
i
dr 2ψi∗ (−∇2 )ψi =
ni
i
ni
dr |ψi |2 (r · ∇)V
(8.12)
i
Let the electronic potential be composed of electron-nuclear, Veв€’n , electron-electron, Veв€’e , and exchangecorrelation, Вµxc , parts. The electron-nuclear potential is
Veв€’n = в€’
Z
|r в€’ R|
(8.13)
Now we want to differentiate Equation (8.13) WRT r
(r · ∇r )
Z
Z
Z
Z
=в€’
− (R · ∇R )
= Ve−n − (R · ∇R )
|r в€’ R|
|r в€’ R|
|r в€’ R|
|r в€’ R|
(8.14)
The contribution to Equation (8.12) is therefore
ni
dr |ψi |2 (r · ∇r )
i
в€’Z
ZQe
= −Ee−n + (R · ∇R )
|r в€’ R|
|r в€’ R|
(8.15)
where Qe is the electronic charge, and Eeв€’n the electron-nuclear potential energy. The electron-electron contribution to Equation (8.12) is
ni
dr |ψi |2 (r · ∇r )
i
dr
(r )
=
|r в€’ r |
dr (r)
dr (r )(r · ∇r )
1
|r в€’ r |
(8.16)
Due to the symmetry of the integral, one may interchange r and r on the RHS, in which case the RHS becomes
1
2
dr (r)
dr (r ) [(r · ∇r ) + (r ·∇r )]
Since
[(r · ∇r ) + (r ·∇r )]
1
|r в€’ r |
1
1
=в€’
|r в€’ r |
|r в€’ r |
equation (8.16) becomes
ni
dr |ψi |2 (r · ∇r )
dr
i
13 Periodic
boundary conditions may achieve the same.
(r )
1
=в€’
|r в€’ r |
2
dr dr
(r) (r )
= в€’Eeв€’e
|r в€’ r |
(8.17)
44
8
VIRIAL THEOREM
where Eeв€’e is the electron-electron potential energy. Finally, Equation (8.12) may be rearranged to
2T + Eeв€’n + Eeв€’e в€’
dr (r)(r · ∇r )µxc = Qe (R · ∇R )
Z
|r в€’ R|
(8.18)
where Вµxc is the exchange correlation potential. The total energy is
E = T + Eeв€’n + Eeв€’e + Exc
(8.19)
If the density is varied ОґE = 0. Hence, if the nuclear position is changed the changes in components of the energy
due to the resulting adjustment of charge density cancel (the Hellman-Feynman theorem) and
R · ∇E = R
dE
= R · ∇Ee−n = 0
dR
(8.20)
Finally, Equation (8.18) becomes
2T + Eeв€’n + Eeв€’e в€’
dr (r)(r · ∇r )µxc = 0
(8.21)
which is the virial theorem when there is an exchange correlation contribution to the total energy.
8.2
Spherical charge density approximation
The electron-electron interaction energy is
Eeв€’e =
1
2
dr
S
r>
2
r>
dr>
=
2
r<
dr<
0
0
S
r2 dr
(r)
Q(r)
r
(r )2 dr
(r )
=
0
(8.22)
r2 dr (r) = Q, with the spherical part of 1/|r в€’ R| =
If the charge density is spherical, normalised such that
1/r> the electron-electron interaction energy is
Eeв€’e
(r) (r )
|r в€’ r |
dr
(r> ) (r< )
=
r>
S
r2 dr
0
(r)
r
r
(r )2 dr
(r )
0
r
where
Q(r)
=
(8.23)
0
where Q is the charge inside r. The Coulomb potential energy is
S
r2 dr
Ecoulв€’tot = Eeв€’e + Eeв€’n =
0
(r)
[Q(r) в€’ Z]
r
(8.24)
The electron potential is
S
Z
(r )
+
dr
+ Вµxc (r)
|r в€’ R|
|r
в€’r |
0
which in the spherical approximation may be written in a number of ways
V (r) = в€’
V (r) = в€’
=в€’
=в€’
Z
1
+
r
r
r
r
r
S
(r )2 dr
(r )2 dr
(r )2 dr
(r ) +
0
S
0
(r )
+ Вµxc (r)
r
(r )2 dr
(r ) +
0
1
Zв€’
r
Z
+
r
S
(r )2 dr
r
(r )
+
r
r
(r )2 dr
0
(8.25)
(r )
(r )
+ Вµxc (r)
r
1
1
в€’
r
r
+ Вµxc (r)
(8.26a)
(8.26b)
(8.26c)
8.2
Spherical charge density approximation
45
If the SchrВЁ
odinger equation
2
∇2 ψi + V ψi = εi
2m
is premultiplied by П€iв€— and summed over the occupied orbitals, the kinetic energy is obtained
в€’
2
T =в€’
2m
dr ψi∗ ∇2 ψi =
ni
i
ni Оµi в€’
dr (r)V (r)
(8.27)
i
If Equation (8.26b) is differentiated WRT r, the contributions from the limits of the integrals cancel and
r
dV
1
dВµxc (r)
= [Z в€’ Q(r)] + r
dr
r
dr
(8.28)
Hence the integral
r2 dr (r)r
dV
dr
r2 dr
=
(r)
[Z в€’ Q(r)] +
r
r2 dr (r)r
= в€’Ecoulв€’tot +
r2 dr (r)r
dВµxc (r)
dr
dВµxc (r)
dr
(8.29)
rearranging
r2 dr (r)r
Ecoulв€’tot = в€’
dV
+
dr
r2 dr (r)r
dВµxc (r)
dr
(8.30)
The exchange correlation contribution may integrated
d(Вµxc ) = Вµxc d + dВµxc
and
d Оµxc
= Вµxc
d
в†’
d(Вµxc ) = d( Оµxc ) + dВµxc
в†’
(8.31)
Вµxc d = d( Оµxc )
(8.32)
whence
Then
r2 dr (r)r
where Оё =
S
0
r2 dr dВµxc = (Вµxc в€’ Оµxc )
dВµxc (r)
=
dr
dВµxc r3 (r) =
r3 dОё
(8.33)
(8.34)
dВµxc = (Вµxc в€’ Оµxc )|S0 . Partial integration yields
r2 dr (r)r
dВµxc (r)
dr
= r3 Оё|S0 в€’
3r2 dr Оё
S
= S 3 (S) [Вµxc (S) в€’ Оµxc (S)] в€’ 3
r2 dr (Вµxc в€’ Оµxc )
(8.35)
0
whence Equation (8.30) becomes, since the surface terms are zero for a free atom
Ecoulв€’tot = в€’
r2 dr (r)r
dV
+3
dr
S
r2 dr (Оµxc в€’ Вµxc )
(8.36)
0
Equation (8.11) may be used to insert the kinetic energy
S
r2 dr (Оµxc в€’ Вµxc )
Ecoulв€’tot = в€’2T + 3
0
(8.37)
46
8
VIRIAL THEOREM
Consequently the total energy, equal to Ecoulв€’tot + T + Exc , is
Etot = в€’T + 4Exc в€’ 3О“xc в€’ S 3 (S) [Оµxc (S) в€’ Вµxc (S)]
where
(8.38)
S
r2 dr Вµxc
О“xc =
(8.39)
0
There are, by definition, no surface terms for core states therefore, from Equation (8.11), the core kinetic energy
given by
dV
(8.40)
2Tc = r2 dr c (r)r
dr
and this may be used in Equation (8.30), leading to
S
r2 dr (Оµxc в€’ Вµxc ) = 0
2T + Ecoulв€’tot в€’ 3
(8.41)
0
8.3
Virial theorem from the Dirac equation
Now let’s do the same thing with the Dirac equation to derive an interesting result. The Dirac equation and its
conjugate are
−i cα · ∇ψ + βmc2 ψ =(E − V )ψ
в€—
в€—
2
i cα · ∇ψ + βmc ψ =(E − V )ψ
(8.42a)
в€—
(8.42b)
If one premultiplies the first of these equations by ψ ∗ r · ∇, one finds
−i (ψ ∗ r · ∇)(cα · ∇ψ + ψ ∗ βmc2 r · ∇ψ
=
ψ ∗ r · ∇(E − V )ψ
=
ψ ∗ (E − V )(r · ∇ψ) − ψ ∗ ψ(r · ∇V )
(8.43)
If the conjugate of the Dirac equation is used to substitute for П€ в€— (E в€’ V ) on the RHS, the terms involving ОІ
are found to cancel and one obtains
ψ ∗ (r · ∇) (cα · ∇ψ) + (cα · ∇ψ ∗ ) (r · ∇)ψ =
1 в€—
ψ ψ(r · ∇V )
i
(8.44)
Now consider
G = ∇· [ψ ∗ α(r · ∇)ψ] =
i
П€ в€— О±i xj
=
ij
∂
[ψ ∗ αi (r · ∇)ψ]
∂xi
∂ ψ
∂ψ ∗ ∂ψ
+ О±i xj
+
∂xi ∂xj
∂xi ∂xj
2
П€ в€— О±i
i
The LHS of Equation (8.44) may be written
П€ в€— О± i xj
ij
∂2ψ
∂ψ ∗ ∂ψ
+ О±i xj
∂xi ∂xj
∂xi ∂xj
which is just
П€ в€— О±i
Gв€’
i
∂ψ
П€
∂xi
or
G − ψ ∗ (α · ∇)ψ
∂ψ
П€
∂xi
(8.45)
8.4
Kinetic energy from radial Dirac Equation
47
Equation (8.44) may therefore be written
∇· [ψ ∗ cα(r · ∇)ψ] − ψ ∗ (cα · ∇)ψ =
1 в€—
ψ ψ(r · ∇V )
i
(8.46)
If one integrates over volume
dr ∇· [ψ ∗ α(r · ∇)ψ] −
dr ψ ∗ (α · ∇)ψ =
1
i c
dr ψ ∗ ψ(r · ∇V )
(8.47)
dS· [ψ ∗ cα(r · i ∇)ψ]
(8.48)
or, using the divergence theorem and rearranging
dr ψ ∗ (cα · i ∇)ψ +
dr |ψ|2 (r · ∇V ) =
The total free electron energy (including rest mass), T , from the Dirac equation is given by
dr П€ в€— (cО± В· p)П€ + mc2
dr П€ в€— ОІП€
dr П€ в€— (cО± В· p)П€ = T в€’ mc2
dr П€ в€— ОІП€
T =
whence
If the rest mass is subtracted from the first term on the RHS (the usual kinetic energy is T = T в€’ mc2 ) and
added to the second term on the RHS, one finds
dr П€ в€— (cО± В· p)П€ = T в€’ mc2
dr П€ в€— (ОІ в€’ 1)П€
(8.49)
and Equation (8.48) becomes the quantum relativistic virial theorem
dr П€ в€— (cО± В· p)П€ в€’
dr |ψ|2 (r · ∇V ) =
dSВ· [П€ в€— cО±(r В· p)П€]
(8.50)
T + 1в€’ОІ mc2
The first term is also equal to T + 1 в€’ ОІ mc2 , but not 2T as in the non-relativistic limit when 1 в€’ ОІ mc2 в†’ T .
The second term is equal to V for a Coulomb potential but not for an exchange-correlation potential (see
Equation (8.31)). Since the surface terms vanish for free atoms, the RHS of Equation (8.50) vanishes.
8.4
Kinetic energy from radial Dirac Equation
We start with the coupled radial Dirac equations from Equation (7.25)
Pi = в€’
Оє
Pi + 2Mi Qi
r
(8.51a)
Оє
Qi = Qi + (V в€’ Оµi )Pi
r
(8.51b)
where Mi is given by Equation (7.24) From Equation (8.51b)
dr (Оµi в€’ V )Pi2
=
в€’
=
в€’Qi Pi |S0 +
dr Qi Pi + Оє
dr
Pi Qi
r
dr Pi Qi + Оє
dr
Pi Qi
r
(8.52)
48
8
VIRIAL THEOREM
and if Equation (8.51a) is used to substitute for P in Equation (8.51b) one finds
dr (Оµi в€’ V )Pi2 = 2Mi
therefore
dr (Оµi в€’ V ) Pi2 +
Q2i
c2
dr Q2i в€’ Qi Pi |S0
dr Q2i m +
=2
or
Оµi в€’
r2 dr П€i2 V = 2
dr Q2i m +
Оµi в€’ V
c2
Оµi в€’ V
c2
(8.53)
в€’ Qi Pi |S0
(8.54)
в€’ Qi Pi |S0
(8.55)
Summing over all orbitals with occupation numbers, ni ,
ni Оµi в€’
dr (r)V (r) = 2
ni
dr Q2i m +
i
i
Оµi в€’ V
c2
ni Qi Pi |S0
в€’
(8.56)
i
Now
cО± В· p = (H в€’ V ) в€’ mc2 ОІ
(8.57)
Therefore, for the eigenstate labelled by i with total energy eigenvalue wi
П€i |cО± В· p|П€i
= wi в€’
r2 dr |П€i |2 V (r) в€’ mc2
= wi в€’
dr |П€i |2 V (r) в€’ mc2
= wi в€’
r2 dr |П€i |2 V (r) в€’ mc2
Q2i
c2
Q2
Pi2 + 2i
c
Pi2 в€’
dr
dr
Q2i
c2
+
2mc2
=
(wi в€’ mc2 ) в€’
dr
r2 dr (П€u )2i в€’ (П€l )2i
r2 dr |П€i |2 V (r) + 2m
dr Q2i
(8.58)
The eigenvalues with the rest mass subtracted are Оµi = wi в€’ mc2
П€i |cО± В· p|П€i = Оµi в€’
r2 dr |П€i |2 V (r) + 2m
dr Q2i
(8.59)
Similarly, the second contribution to the kinetic energy is
mc2 П€i |ОІ|П€i = mc2 в€’ 2m
dr Q2i
в†’
mc2 П€i |ОІ в€’ 1|П€i = в€’2m
dr Q2i
(8.60)
and the free electron energy of the state is
Ti = c П€i |О± В· p|П€i + mc2 П€i |ОІ|П€i = wi в€’
r2 dr |П€i |2 V (r)
(8.61)
However, the rest mass is just a constant and it is subtracted to yield the kinetic energy
Ti = c П€i |О± В· p|П€i + mc2 П€i |(ОІ в€’ 1)|П€i = Оµi в€’
r2 dr |П€i |2 V (r)
(8.62)
8.4
Kinetic energy from radial Dirac Equation
Das
-55,982
49
Snow
-56,116
Desclaux
-56,021
Equation (8.65) for T
-56,140
Table 8.1: The total electronic energy of a free uranium atom (Ry) from a variety of relativistic calculations.
and total kinetic energy
ni c П€i |О± В· p|П€i + mc2
T =
ni П€i |(ОІ в€’ 1)|П€i =
i
i
r2 dr (r)V (r)
ni Оµi в€’
KE-1
(8.63)
i
If Equation (8.59) is summed over all occupied orbitals, one finds
cО± В· p =
ni Оµi в€’
dr (r)V (r) + 2m
dr Q2i
ni
i
(8.64)
i
Equation (8.56) may be used to replace the first two terms on the RHS of Equation (8.64), in which case
cО± В· p = 2
dr Q2i 2m +
ni
i
Оµi в€’ V
c2
ni Qi Pi |S0
в€’
(8.65)
i
The non-relativistic limit c в†’ в€ћ is
cО± В· p = 4m
ni
dr Q2i в€’
i
ni Qi Pi |S0
ONLY FOR c в†’ в€ћ
i
Comparison, with Equation (8.64) yields, if the surface term is zero
T = 2m
ni
dr Q2i =
i
r2 dr (r)V (r) =
ni Оµi в€’
i
1
cО± В· p
2
ONLY FOR c в†’ в€ћ
(8.66)
Equations (8.60) and (8.65) may be combined to produce another expression for the kinetic energy
T = cО± В· p + mc2 ОІ в€’ 1 = 2
ni
dr Q2i m +
i
Оµi в€’ V
c2
ni Qi Pi |S0
в€’
KE-2
(8.67)
i
The composition of the relativistic kinetic energy is obtained from
T = cО± В· p + mc2 ОІ
в†’
T = T в€’ mc2 = cО± В· p + mc2 (ОІ в€’ 1)
where for c в†’ в€ћ, cО± В· p в†’ 2T and mc2 (1 в€’ ОІ) в†’ T . For a uranium atom, Equation (8.64) yields about
69, 000Ry for cО± В· p /214 , whereas Equation (8.67) yields about 62, 500Ry for the kinetic energy. Further, in
the latter case, the total energy is minus cО± В· p /2 for an isolated atom (surface terms absent) whereas it is not
equal to в€’T . A summary of the results of relativistic calculations (Das, Ramana, & Rajagopal 1980, Desclaux
1973, Snow, Canfield & Waber n.d., Kagawa 1975) is shown in Table 8.4.
Another expression for the relativistic kinetic energy may be obtained from Equations (8.50) and 8.60. When
the surface term is zero
ni П€i |cО± В· p|П€i + mc2
T =
i
where
14 The
ni П€i |ОІ в€’ 1|П€i =
i
r2 dr (r)r
dV
в€’
dr
ni
dr Q2i
KE-3
i
is the electron density.
exact amount depends, of course, upon other approximations and details, such as the approximation for Exc .
(8.68)
50
9
9.1
c
c
c
c
9.2
c
c
c
9
Programs
One dimensional SchrВЁ
odinger equation - free electrons
program onedsc
free electrons
integer i,jri
parameter (jri=500)
double precision x(jri),psi(jri,2),dx,e
parameter (e=1.D0)
parameter(dx = 0.05D0)
x-mesh
do i =1,jri
x(i) = dx*float(i-1)
enddo
boundary condition 1
psi(1,1) = 0.D0
psi(2,1) = 1.D-2
do i = 3,jri
psi(i,1) = 2.D0*psi(i-1,1) - psi(i-2,1) - dx*dx*e*psi(i-1,1)
enddo
boundary condition 2
psi(1,2) = 1.D0
psi(2,2) = 1.D0-1.D-20
do i = 3,jri
psi(i,2) = 2.D0*psi(i-1,2) - psi(i-2,2) - dx*dx*e*psi(i-1,2)
enddo
do i=1,jri
write(*,*) x(i),psi(i,1),psi(i,2)
enddo
end
Three dimensional SchrВЁ
odinger equation - free electrons
program radsch
free electrons
integer i,jri,l
parameter (jri=500)
double precision r(jri),psi(jri,2),dx,e
parameter (e=1.D0)
parameter(dx = 0.05D0)
l = 0
r-mesh
do i =1,jri
r(i) = dx*float(i-1)
enddo
boundary condition 1
psi(1,1) = -20.D0
PROGRAMS
9.3
c
9.3
Euler method for a system
psi(2,1) = -10.D0
do i = 3,jri
psi(i,1) = (2.D0*(1.D0+dx/r(i-1))*psi(i-1,1) - psi(i-2,1)
+ - dx*dx*( e-float( l*(l+1) )/r(i-1)**2 )*psi(i-1,1))/
+ (1.D0 + 2.D0*dx/r(i-1))
enddo
boundary condition 2
psi(1,2) = 1.D0
psi(2,2) = 1.D0-1.D-20
do i = 3,jri
psi(i,2) = (2.D0*(1.D0+dx/r(i-1))*psi(i-1,2) - psi(i-2,2)
+ - dx*dx*( e-float( l*(l+1) )/r(i-1)**2 )*psi(i-1,2))/
+ (1.D0 + 2.D0*dx/r(i-1))
enddo
do i=1,jri
write(*,*) r(i),psi(i,1),psi(i,2)
enddo
end
Euler method for a system
program Euler
parameter(nrd=2503)
integer nrd,n,i
double precision x(nrd),y(nrd),z(nrd),yp(nrd),zp(nrd)
double precision h
c integration of y’’-y=x with z=y’ ie z’=y+x and z=y’ system
write(*,*) "Euler integration of z’=y+x and z=y’ system"
c initialize
do n=11,2501,200
x(1)=0.D0
y(1)=0.D0
z(1)=1.D0
h=1.D0/float(n-1)
do i=1,n
yp(i) = z(i)
zp(i) = x(i)+y(i)
x(i+1)=x(i)+h
y(i+1)= y(i) + h*yp(i)
z(i+1)= z(i) + h*zp(i)
enddo
write(*,1) n,x(i-1),y(i-1),z(i-1)
1 format(’grid points= ’,i4,’ x= ’,f10.6,’ y= ’,f10.6,
+ ’ z= ’,f10.6)
enddo
end
51
52
9.4
9
Runge-Kutta for a system
program runge
parameter(nrd=303)
integer nrd,n,i
double precision x(nrd),y(nrd),z(nrd),yp(nrd),zp(nrd)
double precision xk(4),xl(4)
double precision h
c integration of y’’-y=x with z=y’ ie z’=y+x and z=y’ system
write(*,*) "Runga-Kutta integration of z’=y+x and z=y’ system"
do n=9,31,5
x(1)=0.D0
y(1)=0.D0
z(1)=1.D0
h=1.D0/float(n-1)
do i=1,n
xk(1) = h*z(i)
xl(1) = h*(x(i)+y(i))
xk(2) = h*(z(i)+0.5D0*xl(1))
xl(2) = h*(x(i)+0.5D0*h+y(i)+0.5D0*xk(1))
xk(3) = h*(z(i)+0.5D0*xl(2))
xl(3) = h*(x(i)+0.5D0*h+y(i)+0.5D0*xk(2))
xk(4) = h*(z(i)+xl(3))
xl(4) = h*(x(i)+h+y(i)+xk(3))
x(i+1)=x(i)+h
y(i+1)= y(i) + (xk(1)+2.D0*xk(2)+2.D0*xk(3)+xk(4))/6.D0
z(i+1)= z(I) + (xl(1)+2.D0*xl(2)+2.D0*xl(3)+xl(4))/6.D0
enddo
write(*,101) n,x(i-1),y(i-1),z(i-1)
101 format(’grid points= ’,i3,’ x= ’,f10.6,’ y= ’,f10.6,
+ ’ z= ’,f10.6)
enddo
end
9.5
Free electron SchrВЁ
odinger equation - revisited
PROGRAM FESCRO
DOUBLE PRECISION DX,E,R0
PARAMETER (R0=0.2D0,DX=0.005D0,E=0.25D0)
INTEGER NRD,L,I,M,N,NRDM,ID
PARAMETER (NRD=1000,N=2,M=4)
PARAMETER (NRDM=NRD-1)
DOUBLE PRECISION P(NRD),RAD(NRD),Y(N,NRD),Q(NRD),PS(0:5,NRD)
DOUBLE PRECISION XK(N,M),A,B,C,D,Y1,Y2,S
DOUBLE PRECISION G,SIMF
C---------------------------------------C
SET UP RADIAL MESH
PROGRAMS
9.5
Free electron SchrВЁ
odinger equation - revisited
DO I=1,NRD
RAD(I) = R0*DEXP(DX*FLOAT(I-1))
ENDDO
C
ANGULAR MOMENTUM
DO L = 0,5
C---------------------------------------C RUNGE-KUTTA
C---------------------------------------P(1) = 1.D0
Q(1) = 0.D0
Y(1,1)=1.D0
Y(2,1) = 0.D0
DO I = 1, NRDM
Y1 = Y(1,I)
Y2 = Y(2,I)
A = FLOAT(L+1)
B = RAD(I)
D = -E*RAD(I)
XK(1,1) = DX*G(1,Y1,Y2,A,B,D)
XK(2,1) = DX*G(2,Y1,Y2,A,B,D)
B = RAD(I)*DEXP(DX/2.D0)
D = -E*RAD(I)*DEXP(DX/2.D0)
XK(1,2)=DX*G(1,(Y1+0.5D0*XK(1,1)),(Y2+0.5D0*XK(2,1)),A,B,D)
XK(2,2)=DX*G(2,(Y1+0.5D0*XK(1,1)),(Y2+0.5D0*XK(2,1)),A,B,D)
XK(1,3)=DX*G(1,(Y1+0.5D0*XK(1,2)),(Y2+0.5D0*XK(2,2)),A,B,D)
XK(2,3)=DX*G(2,(Y1+0.5D0*XK(1,2)),(Y2+0.5D0*XK(2,2)),A,B,D)
B = RAD(I+1)
D = -E*RAD(I+1)
XK(1,4) = DX*G(1,(Y1+XK(1,3)),(Y2+XK(2,3)),A,B,D)
XK(2,4) = DX*G(2,(Y1+XK(1,3)),(Y2+XK(2,3)),A,B,D)
DO ID =1,2
Y(ID,I+1) = Y(ID,I) +
+ (XK(ID,1)+2.D0*(XK(ID,2)+XK(ID,3))+XK(ID,4))/6.D0
ENDDO
P(I+1) = Y(1,I+1)
Q(I+1) = Y(2,I+1)
ENDDO
DO I=1,NRD
Q(I)=P(I)*P(I)
ENDDO
S = SIMF(RAD,Q,NRD,DX)
DO I=1,NRD
53
54
9
P(I)=P(I)/(DSQRT(S)*RAD(I))
PS(L,I) = P(I)
ENDDO
ENDDO
DO I=1,NRD
WRITE(*,1) RAD(I)*DSQRT(E),(PS(L,I),L=0,5)
1 FORMAT(7F12.4)
ENDDO
END
DOUBLE PRECISION FUNCTION G(ID,Y1,Y2,A,B,D)
INTEGER ID
DOUBLE PRECISION Y1,Y2,A,B,D
IF(ID.EQ.1)THEN
G = A*Y1 + B*Y2
ELSEIF(ID.EQ.2)THEN
G = -A*Y2 + D*Y1
ENDIF
RETURN
END
DOUBLE PRECISION FUNCTION SIMF(X,YT,N,DX)
C SIMPSON’S RULE FOR LOG MESH
C----------------------------------------------------------INTEGER NRD,N,I
PARAMETER (NRD=1000)
DOUBLE PRECISION Y(NRD),X(NRD),YT(NRD),DX,SIMF
C----------------------------------------------------------C
MULTIPLIES BY X TO ACCOUNT FOR LOG MESH
DO I = 1,N
Y(I) = YT(I)*X(I)
ENDDO
C----------------------------------------------------------SIMF = 0.D0
DO I = 2,N-1,2
SIMF = SIMF + (Y(I-1)+4.D0*Y(I)+Y(I+1))
ENDDO
SIMF = DX*SIMF/3.D0
IF(MOD(N,2).EQ.0)THEN
SIMF = SIMF + DX*(5.D0*Y(N)+8.D0*Y(N-1)-Y(N-2))/12.D0
ENDIF
RETURN
END
9.6
Hydrogenic bound states
INCLUDE "./PARAMETERS"
DOUBLE PRECISION ZN,FN,XL,XJ,R0,TWOZ,RWS,CSQ,XKAP,CC
PROGRAMS
9.6
Hydrogenic bound states
INTEGER JRI,I,J,NZ,L,LIGHT
DOUBLE PRECISION RAD(NRD),V(NRD),DX
COMMON/RADIAL/RAD,V,DX,JRI
C
OPEN (UNIT=7,FILE=’H30.DAT’,FORM=’FORMATTED’)
C
REWIND 7
C---------------------------------------------------------------C SET UP CONSTANTS
C Z = NUCLEAR CHARGE, TWOZ = NUCLEAR POTENTIAL TIMES THE RADIUS IN ATOMIC UNITS
C C = VELOCITY OF LIGHT
C FN = PRINCIPLE QUANTUM NUMBER , NODEC = NUMBER OF NODES = FN -L -1
C L = ANGULAR MOMEMTUM
C---------------------------------------------------------------WRITE(*,*) ’INPUT ATOMIC NUMBER’
READ(*,*) ZN
TWOZ = 2.D0*ZN
WRITE(*,*) ’PRINCIPLE QUANTUM NUMBER’
READ(*,*) FN
WRITE(*,*) ’ANGULAR MOMENTUM QUANTUM NUMBER’
READ(*,*) XL
WRITE(*,*) ’XJ’
READ(*,*) XJ
C---------------------------------------------------------------C
SET UP LOGARITHMIC RADIAL MESH AND POTENTIAL: JRI= NUMBER OF MESH POINTS
C
DX = LOGARITHMIC INCREMENT
C
R0 = RADIUS OF FIRST MESH POINT
C---------------------------------------------------------------WRITE(*,*) ’INPUT NUMBER OF RADIAL GRID POINTS’
READ(*,*) JRI
IF (JRI.GT.NRD)THEN
WRITE(*,*) ’JRI GREATER THAN DIMENSIONS - INCREASE NRD’
STOP
ENDIF
WRITE(*,*) ’INPUT ATOM RADIUS IN ANGSTROM’
READ(*,*) RWS
WRITE(*,*) "WIGNER-SEITZ RADIUS IN ANGSTROM= ",RWS
RWS = RWS/AUC
C CONVERT TO ATOMIC UNITS
R0 = 1.D-9
DX = DLOG(RWS/R0)/FLOAT(JRI-1)
WRITE(*,*) "WIGNER-SEITZ RADIUS IN ATOMIC UNITS= ",RWS
WRITE(*,*) "
R0= ",R0," STEP SIZE, DX= ",DX
DO I=1,JRI
RAD(I) = R0*DEXP(DX*FLOAT(I-1))
ENDDO
C---------------------------------------------------------------DO I = 1,JRI
V(I) = -TWOZ/RAD(I)
ENDDO
C---------------------------------------------------------------WRITE(*,*) ’VELOCITY OF LIGHT= ’,C,’
IN ATOMIC (RYDBERG) UNITS’
55
56
WRITE(*,*) ’DO YOU WANT TO CHANGE IT?, NO=0 OR YES=1 ’
READ(*,*) LIGHT
IF(LIGHT.EQ.1)THEN
WRITE(*,*) ’ENTER VELOCITY OF LIGHT’
READ(*,*) CC
CSQ = CC*CC
ELSE
CSQ = C*C
ENDIF
L = INT(XL)
XKAP = (XJ+0.5D0)*2.D0*(XL-XJ)
CALL GETEIG(TWOZ,FN,L,XKAP,CSQ)
C--------------------------------------------------WRITE(*,*) ’
DONE
’
END
SUBROUTINE GETEIG(TWOZ,FN,L,XKAP,CSQ)
INCLUDE "./PARAMETERS"
C EIGENVALUE AND WAVE FUNCTION FOR GIVEN QUANTUM NUMBERS
DOUBLE PRECISION Y(NRD),XNORM
DOUBLE PRECISION E,ETOP,EBOT,EX,FN,TWOZ,DEL,TOL1,TOL2,GAM
DOUBLE PRECISION CSQ,DE,ERR,ET,EB,FAC,QOUT,QIN,XKAP
PARAMETER (TOL1=1.D-2,TOL2=1.D-12)
INTEGER JRI,L,NTP1,NTP2
INTEGER NODE,NODEC
DOUBLE PRECISION P(NRD),Q(NRD),RAD(NRD),V(NRD),DX
COMMON/RADIAL/RAD,V,DX,JRI
DOUBLE PRECISION SIMF
DE = 0.1D0
C----------------------------------------------------------C ETOP EBOT DEFINE ENERGY SCAN RANGE
PRINT*,’M ’,JRI
ETOP = 40.0D0
EBOT = - (TWOZ/FN)**2
NODEC = INT(FN)-L-1
GAM = DSQRT(XKAP**2 - TWOZ**2/CSQ)
NTP1 = JRI
NTP2 = JRI
NODE = 10
DEL = DABS(ETOP - EBOT)
C----------------------------------------------------------------------C GET NUMBER OF NODES AND SCAN RANGE WITH CORRECT NUMBER OF NODES
C----------------------------------------------------------------------C TOP OF RANGE
ET = ETOP
EB = EBOT
DO WHILE (DEL.GT.TOL1)
E = (ET + EB)/2.D0
C----------------------------------------------------------------------C SEARCH FOR TURNING POINT
9
PROGRAMS
9.6
Hydrogenic bound states
CALL GETTP(NTP1,NTP2,JRI,L,E,V,RAD)
CALL RKOUT(E,RAD,V,JRI,GAM,NODE,DX,P,Q,NTP1,XKAP,
+ CSQ,TWOZ,QOUT)
IF(NODE.LE.NODEC)THEN
EB = E
ELSEIF(NODE.GT.NODEC)THEN
ET = E
ENDIF
DEL = DABS(ET - EB)
ENDDO
C---------------------------------------------------------------------ETOP = E
C---------------------------------------------------------------------C BOTTOM OF RANGE
ET = ETOP
EB = EBOT
DEL = DABS(ETOP - EBOT)
DO WHILE (DEL.GT.TOL1)
E = (ET + EB)/2.D0
C----------------------------------------------------------------------C SEARCH FOR TURNING POINT
CALL GETTP(NTP1,NTP2,JRI,L,E,V,RAD)
CALL RKOUT(E,RAD,V,JRI,GAM,NODE,DX,P,Q,NTP1,XKAP,
+ CSQ,TWOZ,QOUT)
IF(NODE.GE.NODEC)THEN
ET = E
ELSEIF(NODE.LT.NODEC)THEN
EB = E
ENDIF
DEL = DABS(ET - EB)
ENDDO
EBOT = E
C---------------------------------------------------------------------C SCAN RANGE ACQUIRED - FROM EBOT TO ETOP
C---------------------------------------------------------------------DE = ETOP - EBOT
DO WHILE(DE.GT.TOL2)
E = (ETOP + EBOT)/2.D0
DEL = DABS(ETOP - EBOT)
C SEARCH FOR TURNING POINT
CALL GETTP(NTP1,NTP2,JRI,L,E,V,RAD)
C----------------------------------------------------------------------CALL RKOUT(E,RAD,V,JRI,GAM,NODE,DX,P,Q,NTP1,XKAP,
+ CSQ,TWOZ,QOUT)
CALL RKIN(E,RAD,V,JRI,L,DX,P,Q,NTP1,NTP2,
+ XKAP,CSQ,TWOZ,QIN)
ERR = QOUT - QIN
C PREDICTED ENERGY
IF(ERR.LT.0.D0)THEN
ETOP = E
57
58
ELSE
EBOT = E
ENDIF
DE = DABS(EBOT-ETOP)
ENDDO
C-------------------------------------------------------C NORMALIZE WAVE FUNCTION
C-------------------------------------------------------DO I =1,JRI
Y(I) = P(I)*P(I) + Q(I)*Q(I)/CSQ
ENDDO
XNORM = SIMF(RAD,Y,JRI,DX)
DO I =1,JRI
P(I) = P(I)/DSQRT(XNORM)
Q(I) = Q(I)/DSQRT(XNORM)
RAD(I) = RAD(I)*AUC
Y(I) = Y(I)/XNORM
WRITE(7,FMT=10) RAD(I),Y(I)
10 FORMAT(2F12.7)
ENDDO
WRITE(*,*) ’ EIGENVALUE= ’,E
RETURN
END
SUBROUTINE GETTP(NTP1,NTP2,JRI,L,E,V,RAD)
INCLUDE "./PARAMETERS"
C GET THE CLASSICAL TURNING POINT
INTEGER NTP1,NTP2,I,J,JRI,L
DOUBLE PRECISION EMV,E
DOUBLE PRECISION V(NRD),RAD(NRD)
EMV = V(JRI)+ FLOAT(L*(L+1))/RAD(JRI)**2
IF(E-EMV.GE.0.D0)THEN
NTP1 = JRI
NTP2 = JRI
RETURN
ENDIF
I = JRI
DO WHILE (E-EMV.LT.0.D0.AND.I.GT.1)
I = I - 1
EMV = V(I)+ FLOAT(L*(L+1))/RAD(I)**2
NTP1 = I + 1
ENDDO
EMV = V(NTP1)+ (FLOAT(L*(L+1))-EPS)/RAD(NTP1)**2
I = NTP1
DO WHILE (E-EMV.GT.0.D0.AND.I.LE.JRI-1)
I = I + 1
EMV = V(I)+ (FLOAT(L*(L+1))-EPS)/RAD(I)**2
NTP2 = I
ENDDO
RETURN
9
PROGRAMS
9.6
Hydrogenic bound states
END
DOUBLE PRECISION FUNCTION GKUTTA(ID,Y1,Y2,A11,A12,A22,A21)
INTEGER ID
DOUBLE PRECISION Y1,Y2,A11,A12,A21,A22
IF(ID.EQ.1)THEN
GKUTTA = A11*Y1 + A12*Y2
ELSEIF(ID.EQ.2)THEN
GKUTTA = A21*Y1 + A22*Y2
ENDIF
RETURN
END
SUBROUTINE RKIN(E,RAD,V,JRI,L,DX,P,Q,NTP1,NTP2,XKAP,
+ CSQ,TWOZ,QIN)
INCLUDE "./PARAMETERS"
INTEGER JRI,L
INTEGER M,N,ID,NTP1,NTP2
PARAMETER (N=2,M=4)
DOUBLE PRECISION E,Y(N,NRD),P(NRD),Q(NRD),RAD(NRD),V(NRD)
DOUBLE PRECISION XK(N,M),CSQ,TWOZ,XM,XMU,BETA,R2EMV
DOUBLE PRECISION A11,A12,A22,A21,EMV,Y1,Y2,DX,QIN,H,XKAP,XLD
PARAMETER (H = 0.5D0)
DOUBLE PRECISION GKUTTA
DX = -DX
XM = 0.5D0*(1.D0 + (E - V(NTP2))/CSQ)
XMU = - (FLOAT(L)+0.5D0)**2 - 2.D0*XM*(V(NTP2)-E)*RAD(NTP2)**2
IF(XMU.LT.0.D0.AND.NTP2.LT.JRI)THEN
C FREE ATOM BOUNDARY CONDITIONS
BETA = (XKAP - DSQRT(-XMU) +
+ 0.5D0*( (FLOAT(L)+0.5D0)**2 + XM*V(NTP2)*RAD(NTP2)**2 )/XMU)/
+ (2.D0*XM*RAD(NTP2))
ELSE
C WIGNER-SEITZ D=-L-1 OR D=L
XLD = FLOAT(-L-1)
BETA = (XLD+XKAP+1)/(2.D0*XM*RAD(NTP2))
ENDIF
P(NTP2)= 1.D-8
Q(NTP2) = BETA*P(NTP2)
Y(1,NTP2)= P(NTP2)
Y(2,NTP2) = Q(NTP2)
DO J = 1,NTP2-NTP1
I = NTP2 - J + 1
Y1 = Y(1,I)
Y2 = Y(2,I)
A11 = -XKAP
A12 = (1.D0 + (E - V(I))/CSQ)*RAD(I)
A22 = XKAP
A21 = (V(I)-E)*RAD(I)
XK(1,1) = DX*GKUTTA(1,Y1,Y2,A11,A12,A22,A21)
59
60
9
XK(2,1) = DX*GKUTTA(2,Y1,Y2,A11,A12,A22,A21)
A12 = (1.D0 + (E - 0.5D0*(V(I)+V(I-1)))/CSQ)*RAD(I)*DEXP(DX/2.D0)
A21 = (V(I-1)*RAD(I-1) - E*RAD(I)*DEXP(DX/2.D0))
XK(1,2)=DX*GKUTTA(1,(Y1+H*XK(1,1)),(Y2+H*XK(2,1)),A11,A12,A22,A21)
XK(2,2)=DX*GKUTTA(2,(Y1+H*XK(1,1)),(Y2+H*XK(2,1)),A11,A12,A22,A21)
XK(1,3)=DX*GKUTTA(1,(Y1+H*XK(1,2)),(Y2+H*XK(2,2)),A11,A12,A22,A21)
XK(2,3)=DX*GKUTTA(2,(Y1+H*XK(1,2)),(Y2+H*XK(2,2)),A11,A12,A22,A21)
A12 = (1.D0 + (E - V(I-1))/CSQ)*RAD(I-1)
A21 = (V(I-1)-E)*RAD(I-1)
XK(1,4) = DX*GKUTTA(1,(Y1+XK(1,3)),(Y2+XK(2,3)),A11,A12,A22,A21)
XK(2,4) = DX*GKUTTA(2,(Y1+XK(1,3)),(Y2+XK(2,3)),A11,A12,A22,A21)
DO ID =1,2
Y(ID,I-1) = Y(ID,I) +
+ (XK(ID,1)+2.D0*(XK(ID,2)+XK(ID,3))+XK(ID,4))/6.D0
ENDDO
P(I-1) = Y(1,I-1)
Q(I-1) = Y(2,I-1)
ENDDO
Y1 = P(NTP1)
DO I = NTP1,NTP2
P(I) = P(I)/Y1
Q(I) = Q(I)/Y1
ENDDO
QIN = Q(NTP1)
IF(NTP2.LT.JRI)THEN
DO I = NTP2,JRI
P(I) = P(NTP2)*DEXP(-DSQRT(-E)*(RAD(I)-RAD(NTP2)))
Q(I) = Q(NTP2)*DEXP(-DSQRT(-E)*(RAD(I)-RAD(NTP2)))
ENDDO
ENDIF
DX = -DX
RETURN
END
SUBROUTINE RKOUT(E,RAD,V,JRI,GAM,NODE,DX,P,Q,NTP1,
+ XKAP,CSQ,TWOZ,QOUT)
INCLUDE "./PARAMETERS"
C---------------------------------------C RUNGE-KUTTA
C---------------------------------------DOUBLE PRECISION DX,E,RATFG,S,H
PARAMETER (H = 0.5D0)
INTEGER JRI,NTP1,NTP1M,I,M,N,ID,NODE
DOUBLE PRECISION P(NRD),RAD(NRD)
DOUBLE PRECISION Y1,Y2,GAM,RNIE,S1,XKAP,QOUT
PARAMETER (N=2,M=4)
DOUBLE PRECISION Y(N,NRD),Q(NRD),V(NRD),D1,D2,D3
DOUBLE PRECISION XK(N,M),XP(5),GP(5)
DOUBLE PRECISION A11,A12,A22,A21,X,EMV,CSQ,TWOZ
DOUBLE PRECISION GKUTTA
PROGRAMS
9.6
Hydrogenic bound states
P(1) = 1.D-8
Q(1) = P(1)*(GAM+XKAP)*CSQ/TWOZ
Y(1,1)= P(1)
Y(2,1) = Q(1)
NTP1M = NTP1 - 1
NODE = 0
EMV = E - V(1)
DO I = 1,NTP1M
Y1 = Y(1,I)
Y2 = Y(2,I)
A11 = -XKAP
A12 = (1.D0 + (E - V(I))/CSQ)*RAD(I)
A22 = XKAP
A21 = (V(I)-E)*RAD(I)
XK(1,1) = DX*GKUTTA(1,Y1,Y2,A11,A12,A22,A21)
XK(2,1) = DX*GKUTTA(2,Y1,Y2,A11,A12,A22,A21)
A12 = (1.D0 + (E - 0.5D0*(V(I)+V(I+1)))/CSQ)*RAD(I)*DEXP(DX/2.D0)
A21 = ( V(I)*RAD(I) - E*RAD(I)*DEXP(DX/2.D0) )
XK(1,2)=DX*GKUTTA(1,(Y1+H*XK(1,1)),(Y2+H*XK(2,1)),A11,A12,A22,A21)
XK(2,2)=DX*GKUTTA(2,(Y1+H*XK(1,1)),(Y2+H*XK(2,1)),A11,A12,A22,A21)
XK(1,3)=DX*GKUTTA(1,(Y1+H*XK(1,2)),(Y2+H*XK(2,2)),A11,A12,A22,A21)
XK(2,3)=DX*GKUTTA(2,(Y1+H*XK(1,2)),(Y2+H*XK(2,2)),A11,A12,A22,A21)
A12 = (1.D0 + (E - V(I+1))/CSQ)*RAD(I+1)
A21 = (V(I+1)-E)*RAD(I+1)
XK(1,4) = DX*GKUTTA(1,(Y1+XK(1,3)),(Y2+XK(2,3)),A11,A12,A22,A21)
XK(2,4) = DX*GKUTTA(2,(Y1+XK(1,3)),(Y2+XK(2,3)),A11,A12,A22,A21)
DO ID = 1,2
Y(ID,I+1) = Y(ID,I) +
+ (XK(ID,1)+2.D0*(XK(ID,2)+XK(ID,3))+XK(ID,4))/6.D0
ENDDO
P(I+1) = Y(1,I+1)
Q(I+1) = Y(2,I+1)
NODE = NODE + DABS(DSIGN(1.D0,P(I+1))-DSIGN(1.D0,P(I)))/2.D0
ENDDO
IF(NTP1.EQ.JRI)THEN
NTP1M = NTP1
DO I = 1,NTP1M
P(I) = P(I)/P(NTP1)
Q(I) = Q(I)/P(NTP1)
ENDDO
QOUT = Q(NTP1)/P(NTP1)
ELSE
IF(NTP1.EQ.JRI)THEN
NTP1M = JRI
ENDIF
DO I = 1,NTP1M
P(I) = P(I)/P(NTP1)
Q(I) = Q(I)/P(NTP1)
ENDDO
QOUT = Q(NTP1)/P(NTP1)
61
62
9
ENDIF
RETURN
END
DOUBLE PRECISION FUNCTION SIMF(X,YT,N,DX)
C SIMPSON’S RULE FOR LOG MESH
C----------------------------------------------------------INCLUDE "./PARAMETERS"
INTEGER N,I
DOUBLE PRECISION Y(NRD),X(NRD),YT(NRD),DX,SIMF
C----------------------------------------------------------C
MULTIPLIES BY X TO ACCOUNT FOR LOG MESH
DO I = 1,N
Y(I) = YT(I)*X(I)
ENDDO
C----------------------------------------------------------SIMF = 0.D0
DO I = 2,N-1,2
SIMF = SIMF + (Y(I-1)+4.D0*Y(I)+Y(I+1))
ENDDO
SIMF = DX*SIMF/3.D0
C
END CORRECTION FO N EVEN
IF(MOD(N,2).EQ.0)THEN
SIMF = SIMF + DX*(5.D0*Y(N)+8.D0*Y(N-1)-Y(N-2))/12.D0
ENDIF
RETURN
END
9.7
C
Self-consistent field atom local density approximation
INCLUDE "../INCLUDE/PARAMETERS"
DOUBLE PRECISION ZN,XN(LDIM),XL(LDIM),XJ(LDIM),XZ(LDIM),XZS(LDIM)
DOUBLE PRECISION EIG(LDIM),SPINM,ORBM,IONRAD,STONH(NRD),STONI(NRD)
INTEGER JTOT,JRI,I,J,ITER,NZ,NC(LDIM),NV,LL
CHARACTER*4 NAM
CHARACTER*8 NLJ(LDIM)
.. COMMON BLOCKS ..
COMMON/SETUP1/ZN,XN,XL,XJ,XZ,XZS,NLJ,NAM,JTOT
DOUBLE PRECISION RAD(NRD),V(NRD),VR(NRD),VI(NRD),DVDR(NRD),DX,ABE
COMMON/RADIAL/RAD,V,VR,VI,DX
DOUBLE PRECISION RHOR2N(NRD),RHOR2O(NRD),
+ RHOR2(NRD),QOFR(NRD),VALRHO(NRD),VDUM(NRD)
COMMON/CHARGE/RHOR2N,RHOR2O,RHOR2,VALRHO
DOUBLE PRECISION BZ,VRS(NRD),RHOS(NRD)
COMMON/MAG1/BZ,RHOS,VRS
DOUBLE PRECISION PWAVE(LODIM,NRD),QWAVE(LODIM,NRD)
COMMON/MAG2/PWAVE,QWAVE
PROGRAMS
9.7
Self-consistent field atom local density approximation
DOUBLE PRECISION PHI,ALFA,RMS,TOL1,TOL2,DEL,OEIGS,TCH,XION,VOL3
DOUBLE PRECISION PRESS3,CADOTP,VEISUM,CSQ,RNODE(LDIM),STON
PARAMETER (TOL1 = 1.D-8,TOL2=1.D-8)
DOUBLE
DOUBLE
DOUBLE
DOUBLE
DOUBLE
PRECISION
PRECISION
PRECISION
PRECISION
PRECISION
S,R0,TWOZ,EIGSUM,TENEG1,TENEG2,TENEG3
DM2,DMM2,RWS,RWSB,VKE1,V3PV,VPRES1
ELPE,PE,KE1,KE2,RKE,KE3,VKE2,VALPE,VALVIR
VXCT,EXCT,VXCTS,EXCTS,VCT,VIR1,VIR2,VIR3
RHIN(NRD,110),RHOUT(NRD,110),FM(NRD,110),F(NRD)
DOUBLE PRECISION SIMF
DOUBLE PRECISION XKAP(LODIM),FN(LODIM),XMU(LODIM),DUMMY(LODIM)
DOUBLE PRECISION XJORB(LODIM),ORBEIG(LODIM),EIGORB(LODIM)
INTEGER NORB,IORB,KAPPA(LODIM),L(LODIM),JPAIR(LODIM),IORBE(LODIM)
DOUBLE PRECISION XKAPE(LODIM),FNE(LODIM),XMUE(LODIM)
INTEGER LE(LODIM),NITS,IBM,IBMQ
CHARACTER*8 NLO(LODIM),NLOE(LODIM)
DOUBLE PRECISION WFP(LDIM,NRD),WFQ(LDIM,NRD),BMASS
COMMON/WFN/WFP,WFQ
C-STATEMENT FUNCTIONS-------------------------------------------------------------DOUBLE PRECISION BZ1,BZ2
BZ1(LL) = 2.D0/FLOAT(2*LL+3)
BZ2(LL) = 2.D0/FLOAT(2*LL+1)
C CLEAR CHARGE AND SPIN DENSITY ARRAYS
BZ = 0.D0
DO I = 1,NRD
RHOR2O(I)=0.D0
RHOR2N(I)=0.D0
RHOR2(I)=0.D0
RHOS(I) = 0.D0
VRS(I) = 0.D0
ENDDO
C---------------------------------------------------------------C SET UP CONSTANTS
C Z = NUCLEAR CHARGE, TWOZ = NUCLEAR POTENTIAL TIMES THE RADIUS IN ATOMIC UNITS
C C = VELOCITY OF LIGHT
C FN = PRINCIPLE QUANTUM NUMBER , NODEC = NUMBER OF NODES = FN -L -1
C L = ANGULAR MOMEMTUM
C---------------------------------------------------------------CSQ = C*C
RMS = 1.D0
DEL = 1.D0
OEIGS = 0.D0
PHI = 0.20D0
ALFA = 1.2D0 - PHI
63
64
9
PRINT*,’INPUT NUCLEAR CHARGE’
READ*,ZN
TWOZ = 2.D0*ZN
XION = 0.D0
FN = 1.D0
L = 0
C---------------------------------------------------------------C
SET UP LOGARITHMIC RADIAL MESH AND POTENTIAL: JRI= NUMBER OF MESH POINTS
C
DX = LOGARITHMIC INCREMENT
C
R0 = RADIUS OF FIRST MESH POINT
C---------------------------------------------------------------JRI = 1000
IF (JRI.GT.NRD)THEN
WRITE(*,*) ’JRI GREATER THAN DIMENSIONS - INCREASE NRD’
STOP
ENDIF
RWS = 80.0D0
PRINT*, ’
WIGNER-SEITZ RADIUS IN ANGSTROM= ’,RWS
RWS = RWS/AUC
PRINT*, ’WIGNER-SEITZ RADIUS IN ATOMIC UNITS= ’,RWS
VOL3 = 4.D0*PI*RWS**3
PRINT*, ’
3 TIMES VOLUME IN AU= ’,VOL3
C CONVERT TO ATOMIC UNITS
R0 = 1.D-8
DX = DLOG(RWS/R0)/FLOAT(JRI-1)
PRINT*, ’
R0= ’,R0
PRINT*,’
STEP SIZE, DX= ’,DX
PRINT*,’OK?’
READ(*,*)
DO I=1,JRI
RAD(I) = R0*DEXP(DX*FLOAT(I-1))
ENDDO
C---------------------------------------------------------------CALL SETUP
C---------------------------------------------------------------S = 0.D0
WRITE(6,*) ’
J
L
XZ(J)’
WRITE(*,*) ’----------------------------------------------’
DO J=1,JTOT
S=S+XZ(J)
WRITE(6,*) J,NLJ(J),XL(J),XZ(J)
ENDDO
PRINT*, ’TOTAL NUMBER OF ELECTRONS= ’,S,’ IONICITY= ’,XION,’OK?’
READ(*,*)
C---------------------------------------------------------------C GET THOMAS-FERMI DENSITY AND POTENTIAL DUE TO ELECTRONS
PROGRAMS
9.7
Self-consistent field atom local density approximation
65
CALL TFPOT(JRI,ZN,XION,DX,RAD,VR)
CALL TFDENS(ZN,XION,JRI,RAD,RHOR2O,DX)
C CONSTRUCT TOTAL POTENTIAL TO START
DO I = 1,JRI
V(I) = (VR(I)-TWOZ)/RAD(I)
ENDDO
DO I =1,JRI
RHOR2(I)= RHOR2O(I)
ENDDO
CALL TPOT(JRI,EXCT,VXCT,EXCTS,VXCTS,STONH)
C---------------------------------------------------------------C SELF-CONSISTENCY LOOP
ITER = 0
DO WHILE (RMS.GT.TOL1.OR.DEL.GT.TOL2)
ITER = ITER + 1
WRITE(*,*) "
ITERATION NUMBER= ",ITER
IF(ITER.GT.99) WRITE(*,*) "TOO MANY ITERATIONS"
IF(ITER.GT.99) STOP
C---------------------------------------------------------------C INTERPOLATE POTENTIAL
CALL INTERP(JRI)
C--------------------------------------------------------------C GET CHARGE DENSITIES
CALL WAVEEQ(JRI,EIGSUM,NC,KE2,RKE,VKE2,CADOTP,VEISUM,EIG,RNODE)
C----------------------------------------------------------------C MIX CHARGE DENSITIES
ALFA = 0.50
RMS=0.D0
DO I = 1,JRI
RHOR2(I) = ALFA*RHOR2O(I) + (1.D0-ALFA)*RHOR2N(I)
RMS = RMS + (RHOR2O(I) - RHOR2N(I))**2
RHIN(I,ITER) = RHOR2O(I)
RHOUT(I,ITER) = RHOR2N(I)
FM(I,ITER) = RHOR2O(I) - RHOR2N(I)
ENDDO
RMS=RMS/FLOAT(JRI)
DEL = DABS((EIGSUM - OEIGS)/EIGSUM)
OEIGS = EIGSUM
WRITE(6,*) ’RMS ERROR FOR ATOM C.D.= ’,RMS,’
WRITE(6,*) ’EIGENVALUE SUM ERROR
= ’,DEL
C
TCH = SIMF(RAD,RHOR2,JRI,DX)
WRITE(*,*) ’TOTAL ELECTRONIC CHARGE SIMF
TCH = ADLINT(RAD,RHOR2,JRI,DX)
MIXING RATIO ’,ALFA
= ’,TCH
66
C
9
WRITE(*,*) ’TOTAL ELECTRONIC CHARGE ADLINT
= ’,TCH
C GET COULOMB POTENTIAL
CALL COUL(JRI,TWOZ,VCT)
C ADD EXC POTENTIAL
CALL TPOT(JRI,EXCT,VXCT,EXCTS,VXCTS,STONH)
C------------------------------------------------------------------C KE1
DO I=1,JRI
F(I) = - V(I)*RHOR2(I)
ENDDO
KE1 = EIGSUM + SIMF(RAD,F,JRI,DX)
C----------------------------------------------------------------C-----------------------------------------------------------------C KE3
CALL DIFFRP(V,DX,JRI,DVDR,RAD)
DO I = 1,JRI
F(I) = RHOR2(I)*RAD(I)*DVDR(I)
ENDDO
KE3 = SIMF(RAD,F,JRI,DX) - RKE
C----------------------------------------------------------------C ELECTRON-ELECTRON POTENTIAL ENERGY
C FOR SIMINT - LOG MESH
DO I = 1,JRI
F(I) = RHOR2(I)*RAD(I)
ENDDO
CALL SIMINT(DX,F,QOFR,JRI,NRD)
DO I = 1,JRI
F(I) = 2.D0*QOFR(I)*RHOR2(I)/RAD(I)
ENDDO
ELPE = SIMF(RAD,F,JRI,DX)
C-------------------------------------------------------------C TOTAL COULOMB PE
DO I = 1,JRI
F(I) = 2.D0*(QOFR(I)-ZN)*RHOR2(I)/RAD(I)
ENDDO
PE = SIMF(RAD,F,JRI,DX)
C---------------------------------------------------------C TE
TENEG1 = EIGSUM + VCT + EXCT - VXCT
TENEG2 = KE2 + PE + EXCT
TENEG3 = KE3 + PE + EXCT
C---------------------------------------------------------C VIRIAL CHECKS
CALL DIFFRP(V,DX,JRI,DVDR,RAD)
DO I = 1,JRI
F(I) = RHOR2(I)*RAD(I)*DVDR(I)
ENDDO
VIR1 = KE1 + RKE - SIMF(RAD,F,JRI,DX)
PROGRAMS
9.7
C
Self-consistent field atom local density approximation
VIR2 = CADOTP - SIMF(RAD,F,JRI,DX)
VIR2 = CADOTP + PE -3.D0*(EXCT - VXCT) + EXCTS - VXCTS
VIR3 = SIMF(RAD,F,JRI,DX) + PE -3.D0*(EXCT - VXCT) + EXCTS - VXCTS
C---------------------------------------------------------C REPLACE OLD CHARGE DENSITY WITH NEW FOR NEXT ITERATION
DO I = 1,JRI
RHOR2O(I) = RHOR2(I)
ENDDO
C---------------------------------------------------------C
END OF SC LOOP
ENDDO
C---------------------------------------------------------PRINT*,’ IONIC RADIUS’
IONRAD = -1.D0
DO J = 1,JTOT
IF(RNODE(J).GT.IONRAD)THEN
IONRAD = RNODE(J)
ENDIF
ENDDO
IONRAD = IONRAD*AUC
PRINT*,’ IONIC RADIUS ANGSTROM= ’,IONRAD,’
OK?’
READ(*,*)
WRITE(*,*) ’ NUMBER OF ITERATIONS= ’,ITER
C--------------------------------------------------NZ = ZN
WRITE(6,*) "ATOM ",NAM,"ATOMIC NUMBER= ",NZ
WRITE(*,*) "WIGNER-SEITZ RADIUS= ",RWS*AUC,"
R0= ",R0,
+ " STEP SIZE= ",DX
WRITE(*,*) ’
KINETIC ENERGY 1= ’,KE1
WRITE(*,*) ’
KINETIC ENERGY 2= ’,KE2
WRITE(*,*) ’
KINETIC ENERGY 3= ’,KE3
WRITE(*,*) ’
ELECTRON PE= ’,ELPE
WRITE(*,*) ’
COULOMB ENERGY= ’,VCT
WRITE(*,*) ’
POTENTIAL ENERGY= ’,PE
WRITE(*,*) ’
EIGENVALUE SUM= ’,EIGSUM
WRITE(*,*) ’EXCH-CORRELATION ENERGY= ’,EXCT
WRITE(*,*) ’
TOTAL ENERGY 1= ’,TENEG1
WRITE(*,*) ’
TOTAL ENERGY 2= ’,TENEG2
WRITE(*,*) ’
TOTAL ENERGY 3= ’,TENEG3
WRITE(*,*) ’
VIR1= ’,VIR1
WRITE(*,*) ’
VIR2= ’,VIR2
WRITE(*,*) ’
VIR3= ’,VIR3
WRITE(*,*) ’
INT OF Q^2= ’,RKE
WRITE(*,*) ’
DONE
’
67
68
9
PRINT*,’WANT THE BAND MASS?, ENTER 2 FOR YES, ELSE OTHER NUMBER’
READ(*,*) IBMQ
IF(IBMQ.EQ.2)THEN
PRINT*,’ENTER RADIUS FOR BAND MASS IN ANGSTROM’
READ(*,*) RWSB
RWSB = RWSB/AUC
IBM = INT( 1.D0 + DLOG(RWSB/R0)/DX )
DO J = 1,JTOT
BMASS = 2.D0/(RAD(IBM)*WFP(J,IBM)**2)
PRINT*,J,NLJ(J),BMASS,RAD(IBM)*AUC
ENDDO
ENDIF
PRINT*,’WANT EXCHANGE INTS?, ENTER 2 FOR YES, ELSE OTHER NUMBER’
READ(*,*) IBMQ
IF(IBMQ.EQ.2)THEN
DO J = 1,JTOT
DO I = 1,JRI
STONI(I) = STONH(I)*(WFP(J,I)**2+WFQ(J,I)**2/CSQ)**2
ENDDO
STON = - R2EV*SIMF(RAD,STONI,JRI,DX)
PRINT*,J,NLJ(J),’ EXCHANGE INTEGRAL= ’,STON
ENDDO
ENDIF
WRITE(*,*) ’
DONE
’
END
SUBROUTINE SETUP
INCLUDE "../INCLUDE/PARAMETERS"
INTEGER NORB1(10,10,0:6),NORB2(LDIM),XORB3(LDIM),NORB4(LDIM)
INTEGER NL,N,L,I,LL,NZ,ITOT,IJ,JTOT,J
DOUBLE PRECISION SUMM,EXTRA
DOUBLE PRECISION ZN,XN(LDIM),XL(LDIM),XJ(LDIM),XZ(LDIM),XZS(LDIM)
CHARACTER*4 NAM
CHARACTER*8 NLJ(LDIM)
C
.. COMMON BLOCKS ..
COMMON/SETUP1/ZN,XN,XL,XJ,XZ,XZS,NLJ,NAM,JTOT
C---------------------------------------------C TABLE OF ORBITALS IN CHARACTERS
DATA NLJ/’ 1S1/2’,’ 2S1/2’,’ 2P1/2’,’ 2P3/2’,’ 3S1/2’,
+ ’ 3P1/2 ’,’ 3P3/2 ’,’ 4S1/2 ’,’ 3D3/2 ’,’ 3D5/2 ’,
+ ’ 4P1/2 ’,’ 4P3/2 ’,’ 5S1/2 ’,’ 4D3/2 ’,’ 4D5/2 ’,
+ ’ 5P1/2 ’,’ 5P3/2 ’,’ 6S1/2 ’,’ 4F5/2 ’,’ 4F7/2 ’,
+ ’ 5D3/2 ’,’ 5D5/2 ’,
+ ’ 6P1/2 ’,’ 6P3/2 ’,’ 7S1/2 ’,’ 5F5/2 ’,’ 5F7/2 ’,
+ ’ 6D3/2 ’,’ 6D5/2 ’,’ 8S1/2 ’,’ 5G7/2 ’,’ 5G9/2 ’,
+ ’ 6F5/2 ’,’ 6F7/2 ’,
+ ’ 7D3/2 ’,’ 7D5/2 ’,
PROGRAMS
9.7
Self-consistent field atom local density approximation
+ ’ 7P1/2 ’,’ 7P3/2 ’,’ 6G7/2 ’,
+ ’ 6G9/2 ’,
+ ’ 6H9/2 ’,’ 6H11/2 ’,
+ ’ 6H9/2 ’,’ 6H11/2 ’,’ 6H11/2 ’,
+ ’ 7F5/2 ’,’ 7F7/2 ’,’ 7G7/2 ’,’ 7G9/2 ’,’ 7H9/2 ’/
C---------------------------------------------C TABLE OF ELEMENT NAMES
CHARACTER*2 NAMTAB(0:110)
DATA NAMTAB /’EV’,
+
’H ’,’HE’,’LI’,’BE’,’B ’, ’C ’,’N ’,’O ’,’F ’,’NE’,
+
’NA’,’MG’,’AL’,’SI’,’P ’, ’S ’,’CL’,’AR’,’K ’,’CA’,
+
’SC’,’TI’,’V ’,’CR’,’MN’, ’FE’,’CO’,’NI’,’CU’,’ZN’,
+
’GA’,’GE’,’AS’,’SE’,’BR’, ’KR’,’RB’,’SR’,’Y ’,’ZR’,
+
’NB’,’MO’,’TC’,’RU’,’RH’, ’PD’,’AG’,’CD’,’IN’,’SN’,
+
’SB’,’TE’,’I ’,’XE’,’CS’, ’BA’,’LA’,’CE’,’PR’,’ND’,
+
’PM’,’SM’,’EU’,’GD’,’TB’, ’DY’,’HO’,’ER’,’TM’,’YB’,
+
’LU’,’HF’,’TA’,’W ’,’RE’, ’OS’,’IR’,’PT’,’AU’,’HG’,
+
’TL’,’PB’,’BI’,’PO’,’AT’, ’RN’,’FR’,’RA’,’AC’,’TH’,
+
’PA’,’U ’,’NP’,’PU’,’AM’, ’CM’,’BK’,’CF’,’ES’,’FM’,
+
’MD’,’NO’,’LW’,’KT’,’??’, ’??’,’??’,’??’,’??’,’??’/
C-----------------------------------------------NZ = ZN
NAM=NAMTAB(NZ)
PRINT*,’THIS ATOM IS= ’, NAM,’ WITH NUCLEAR CHARGE= ’,NZ, ’OK?’
C
WRITE(*,’(" 1 ")’) !BP
READ(*,*)
C CLEAR ARRAY
DO I =1,LDIM
XORB3(I) = 0.D0
NORB4(I) = 0
NORB2(I) = 0
ENDDO
DO NL = 1,10
DO N = 1,10
DO L=0,6
NORB1(NL,N,L) = 0
ENDDO
ENDDO
ENDDO
C-------------------------------------------C SETS UP MADELUNG’S RULES IN NORB1
DO NL = 1,10
DO N = 1,10
DO L= 0,N-1
IF(NL.EQ.N+L)THEN
NORB1(NL,N,L) = NL
C
WRITE(*,*) N,L,NL,NORB1(NL,N,L)
ENDIF
ENDDO
ENDDO
69
70
9
PROGRAMS
ENDDO
C-------------------------------------------C SETS UP THE PRINCIPLE QUANTUM NUMBER AND ANGLUAR MOMENTUM ACCORDING TO MADELUNG’S RULES
C NORB2 AND NORB4. SETS UP SHELL OCCUPATION IN XORB3.
PRINT*,’
I
NORB2=N
NORB4=L
XORB3=SHELL-OCC’
I=0
PRINT*,’---------------------------------------------------------’
DO NL = 1,10
DO N = 1,10
L = NL - N
IF(L.LE.N-1.AND.NL.EQ.L+N.AND.L.GE.0)THEN
IF(NORB1(NL,N,L).NE.0)THEN
I = I + 1
NORB2(I) = N
NORB4(I) = L
XORB3(I) = FLOAT(2*(2*L+1))
PRINT*, I,NORB2(I),NORB4(I),XORB3(I)
ENDIF
ENDIF
ENDDO
ENDDO
PRINT*,’---------------------------------------------------------’
C-------------------------------------------C COUNTS NUMBER OF OCCUPIED ORBITALS IN FILLED SHELLS AND
PRINT*,’ I
SUMM
NORB2(I)
NORB4(I)
XORB3(I)’
PRINT*,’---------------------------------------------------------’
SUMM = 0.D0
I = 0
DO WHILE (SUMM.LE.FLOAT(NZ))
I = I + 1
SUMM = SUMM + XORB3(I)
PRINT*, I,SUMM,NORB2(I),NORB4(I),XORB3(I)
ENDDO
EXTRA = SUMM - FLOAT(NZ)
XORB3(I) = XORB3(I) - EXTRA
ITOT = I
PRINT*,’---------------------------------------------------------’
PRINT*, ’TOTAL NUMBER OF ORBITALS= ’,ITOT,
+’ EXTRA ELECTRONS= ’,EXTRA
PRINT*,’OUTERMOST ORBITAL OCCUPATION= ’,XORB3(I)
PRINT*,’---------------------------------------------------------’
SUMM = 0.D0
DO I = 1,ITOT
SUMM = SUMM + XORB3(I)
ENDDO
PRINT*,’---------------------------------------------------------’
PRINT*, ’CHECK ON TOTAL NUMBER OF ELECTRONS= ’,SUMM
9.7
Self-consistent field atom local density approximation
PRINT*,’---------------------------------------------------------’
C-------------------------------------------C SETS UP RELATIVISTIC QUANTUM NUMBERS
IJ=0
DO I = 1,ITOT
IJ = IJ + 1
L = NORB4(I)
IF(L.EQ.0)THEN
XJ(IJ) = NORB4 (I) + 0.5D0
XN(IJ) = NORB2 (I)
XL(IJ) = NORB4 (I)
XZ(IJ) = XORB3 (I)
ELSE
XJ(IJ) = NORB4 (I) - 0.5D0
XL(IJ) = NORB4 (I)
XN(IJ) = NORB2 (I)
XZ(IJ) = XORB3 (I)*FLOAT(L)/FLOAT((2*L+1))
IJ = IJ +1
XJ(IJ) = NORB4 (I) + 0.5D0
XL(IJ) = NORB4 (I)
XN(IJ) = NORB2 (I)
XZ(IJ) = XORB3 (I)*FLOAT(L+1)/FLOAT(2*L+1)
ENDIF
ENDDO
PRINT*,’---------------------------------------------------------’
JTOT = IJ
PRINT*,’ TOTAL NUMBER OF J-J COUPLED ORBITALS= ’,JTOT
PRINT*,’---------------------------------------------------------’
PRINT*, ’
XN(J)
XL(J)
+
XJ(J)
XZ(J)’
DO J = 1,JTOT
PRINT*, NLJ(J),XN(J),XL(J),XJ(J),XZ(J)
ENDDO
PRINT*,’OK?’
READ(*,*)
RETURN
END
SUBROUTINE TFDENS(ZN,XION,JRI,RAD,RHOR2O,DX)
INCLUDE "../INCLUDE/PARAMETERS"
DOUBLE PRECISION ZN,XION,DX
DOUBLE PRECISION RAD(NRD),RHOR2O(NRD)
DOUBLE PRECISION ALPHA,BETA,GAMMA,G3,CHARGE
INTEGER JRI,I
DOUBLE PRECISION SIMF
WRITE(*,*) ’ ENTERING THOMAS FERMI DENSITY’,ZN,XION,JRI
ALPHA=0.3058*(ZN-XION)**(1.0/3.0)
71
72
9
BETA=5.8632*DMAX1(0.D0,ZN-XION-2.D0)*ALPHA
GAMMA=DSQRT(4.0+3.2*XION)
G3=GAMMA**3
IF ((ZN-XION).LT.2.0) THEN
G3=0.5*(ZN-XION)*G3
ENDIF
DO I=1,JRI
X=ALPHA*RAD(I)
RHOR2O(I)=BETA*DSQRT(ALPHA*RAD(I))*DEXP(-3.0*ALPHA*RAD(I))
+ + G3*DEXP(-GAMMA*RAD(I))*RAD(I)**2
ENDDO
CHARGE = SIMF(RAD,RHOR2O,JRI,DX)
WRITE(*,*) "CHARGE CHECK IN THOMAS FERMI - Q= ",CHARGE
WRITE(*,*) ’ LEAVING THOMAS FERMI DENSITY’
RETURN
END
SUBROUTINE TFPOT(JRI,ZN,XION,DX,RAD,VR)
C
C **********************************************************************
C
C THOMAS FERMI POTENIAL
C
C **********************************************************************
C
C
.. PARAMETERS ..
INCLUDE "../INCLUDE/PARAMETERS"
C
..
C
.. SCALAR ARGUMENTS ..
DOUBLE PRECISION XION,ZN,DX
INTEGER JRI
C
..
C
.. ARRAY ARGUMENTS ..
DOUBLE PRECISION RAD(NRD),VR(NRD)
C
..
C
.. LOCAL SCALARS ..
DOUBLE PRECISION ELNO,WA,WC,WD,WE
INTEGER I
C
..
C
.. INTRINSIC FUNCTIONS ..
INTRINSIC SQRT
C
..
WRITE(*,*) ’ ENTERING THOMAS FERMI POTENTIAL’,ZN,XION,JRI
ELNO = ZN - XION
WA = ELNO - ZN - 1
C
PROGRAMS
9.7
C
C
C
Self-consistent field atom local density approximation
WRITE(6,8) JRI,ZN,ELNO,XION,WA
8 FORMAT(I5,4F10.4)
DO 10 I = 1,JRI
WC = SQRT((RAD(I)* (ZN+WA)** (1./3.))/0.8853)
WD = WC* (0.60112*WC+1.81061) + 1.
WE = WC* (WC* (WC* (WC* (0.04793*WC+0.21465)+0.77112)+
+
1.39515)+1.81061) + 1.
WC = (ZN+WA)* (WD/WE)**2 - WA
C
C RAD TIMES ELECTRONIC POTENTIAL
C
VR(I) =2.D0*(ZN -WC)
10 CONTINUE
WRITE(*,*) ’LEAVING THOMAS FERMI POTENTIAL’
C
C
C
C
C
C
DO I = 1,JRI
WRITE(6,7) VR(I),RAD(I)
ENDDO
7 FORMAT(2D14.5)
END
SUBROUTINE GETEIG(P,Q,TWOZ,FN,L,JRI,E,NODE,SKAP,XKAP,CSQ,XLOGD,
+ NTP1,NTP2,JNODE)
INCLUDE "../INCLUDE/PARAMETERS"
C EIGENVALUE AND WAVE FUNCTION FOR GIVEN QUANTUM NUMBERS
DOUBLE PRECISION Y(NRD),XNORM
DOUBLE PRECISION E,ETOP,EBOT,EX,FN,TWOZ,R0,DEL,TOL1,TOL2,GAM,XM
DOUBLE PRECISION CSQ,DE,ERR,ET,EB,FAC,QOUT,QIN,XKAP,SKAP,XLOGD
PARAMETER (TOL1=1.D-2,TOL2=1.D-10)
INTEGER JRI,L,NTP1,NTP2,JNODE
INTEGER NODE,NODEC
DOUBLE PRECISION P(NRD),Q(NRD)
DOUBLE PRECISION RAD(NRD),V(NRD),VR(NRD),VI(NRD),DX
COMMON/RADIAL/RAD,V,VR,VI,DX
DOUBLE PRECISION SIMF,ADLINT
DE = 0.1D0
C----------------------------------------------------------C ETOP EBOT DEFINE ENERGY SCAN RANGE
ETOP = 40.0D0
73
74
9
EBOT = - (TWOZ/FN)**2
NODEC = INT(FN+0.001D0)-L-1
GAM =
DSQRT(XKAP**2 - TWOZ**2/CSQ)
NTP1 = JRI
NTP2 = JRI
NODE = 10
DEL = DABS(ETOP - EBOT)
C----------------------------------------------------------------------C GET NUMBER OF NODES AND SCAN RANGE WITH CORRECT NUMBER OF NODES
C----------------------------------------------------------------------C TOP OF RANGE
ET = ETOP
EB = EBOT
DO WHILE (DEL.GT.TOL1)
E = (ET + EB)/2.D0
C----------------------------------------------------------------------C SEARCH FOR TURNING POINT
CALL GETTP(NTP1,NTP2,JRI,L,E,V,RAD)
CALL RKOUT(E,RAD,V,VI,JRI,L,GAM,NODE,DX,P,Q,NTP1,SKAP,XKAP,
+ CSQ,TWOZ,QOUT,XLOGD,JNODE)
IF(NODE.LE.NODEC)THEN
EB = E
ELSEIF(NODE.GT.NODEC)THEN
ET = E
ENDIF
DEL = DABS(ET - EB)
ENDDO
C---------------------------------------------------------------------ETOP = E
C
WRITE(*,*) ’TOP OF RANGE AND NODES ’,ETOP, NODE
C BOTTOM OF RANGE
ET = ETOP
EB = EBOT
DEL = DABS(ETOP - EBOT)
DO WHILE (DEL.GT.TOL1)
E = (ET + EB)/2.D0
C----------------------------------------------------------------------C SEARCH FOR TURNING POINT
PROGRAMS
9.7
Self-consistent field atom local density approximation
CALL GETTP(NTP1,NTP2,JRI,L,E,V,RAD)
CALL RKOUT(E,RAD,V,VI,JRI,L,GAM,NODE,DX,P,Q,NTP1,SKAP,XKAP,
+ CSQ,TWOZ,QOUT,XLOGD,JNODE)
IF(NODE.GE.NODEC)THEN
ET = E
ELSEIF(NODE.LT.NODEC)THEN
EB = E
ENDIF
DEL = DABS(ET - EB)
ENDDO
C---------------------------------------------------------------------EBOT = E
C
WRITE(*,*) ’BOTTOM OF RANGE AND NODES ’,EBOT, NODE
DE = ETOP - EBOT
DO WHILE(DE.GT.TOL2)
E = (ETOP + EBOT)/2.D0
DEL = DABS(ETOP - EBOT)
C SEARCH FOR TURNING POINT AGAIN
CALL GETTP(NTP1,NTP2,JRI,L,E,V,RAD)
C----------------------------------------------------------------------CALL RKOUT(E,RAD,V,VI,JRI,L,GAM,NODE,DX,P,Q,NTP1,SKAP,XKAP,
+ CSQ,TWOZ,QOUT,XLOGD,JNODE)
CALL RKIN(E,RAD,V,VI,JRI,L,DX,P,Q,NTP1,NTP2,
+ SKAP,XKAP,CSQ,TWOZ,QIN)
XM = 0.5D0*( 1.D0 + (E - V(NTP2))/CSQ )
XLOGD = - XKAP - 1.D0 - 2.D0*XM*SKAP*RAD(NTP2)*Q(NTP2)/P(NTP2)
ERR = QOUT - QIN
C PREDICTED ENERGY
IF(ERR.LT.0.D0)THEN
ETOP = E
ELSE
EBOT = E
ENDIF
DE = DABS(EBOT-ETOP)
ENDDO
C-------------------------------------------------------C NORMALIZE WAVE FUNCTION
C-------------------------------------------------------DO I =1,JRI
Y(I) = P(I)*P(I) + Q(I)*Q(I)/CSQ
ENDDO
75
76
9
XNORM = SIMF(RAD,Y,JRI,DX)
DO I =1,JRI
P(I) = P(I)/DSQRT(XNORM)
Q(I) = Q(I)/DSQRT(XNORM)
ENDDO
RETURN
END
SUBROUTINE GETTP(NTP1,NTP2,JRI,L,E,V,RAD)
INCLUDE "../INCLUDE/PARAMETERS"
C GET THE CLASSICAL TURNING POINT
INTEGER NTP1,NTP2,I,J,JRI,L
DOUBLE PRECISION EMV,E
DOUBLE PRECISION V(NRD),RAD(NRD)
EMV = V(JRI)+ FLOAT(L*(L+1))/RAD(JRI)**2
IF(E-EMV.GE.0.D0)THEN
NTP1 = JRI
NTP2 = JRI
RETURN
ENDIF
I = JRI
DO WHILE (E-EMV.LT.0.D0.AND.I.GT.1)
I = I - 1
EMV = V(I)+ FLOAT(L*(L+1))/RAD(I)**2
NTP1 = I + 1
ENDDO
EMV = V(NTP1)+ (FLOAT(L*(L+1))-EPS)/RAD(NTP1)**2
I = NTP1
DO WHILE (E-EMV.GT.0.D0.AND.I.LE.JRI-1)
I = I + 1
EMV = V(I)+ (FLOAT(L*(L+1))-EPS)/RAD(I)**2
NTP2 = I
ENDDO
RETURN
END
DOUBLE PRECISION FUNCTION GKUTTA(ID,Y1,Y2,A11,A12,A22,A21)
INTEGER ID
DOUBLE PRECISION Y1,Y2,A11,A12,A21,A22
IF(ID.EQ.1)THEN
GKUTTA = A11*Y1 + A12*Y2
ELSEIF(ID.EQ.2)THEN
GKUTTA = A21*Y1 + A22*Y2
ENDIF
PROGRAMS
9.7
Self-consistent field atom local density approximation
RETURN
END
SUBROUTINE RKIN(E,RAD,V,VI,JRI,L,DX,P,Q,NTP1,NTP2,SKAP,XKAP,
+ CSQ,TWOZ,QIN)
INCLUDE "../INCLUDE/PARAMETERS"
INTEGER JRI,L,LMAXH
PARAMETER (LMAXH = 4)
INTEGER M,N,ID,JRIM,NTP1,NTP2
PARAMETER (N=2,M=4)
DOUBLE PRECISION E,Y(N,NRD),P(NRD),Q(NRD),RAD(NRD),V(NRD),VI(NRD)
DOUBLE PRECISION XK(N,M),TWOM,CSQ,TWOZ,XM,XMU,BETA
DOUBLE PRECISION A11,A12,A22,A21,EMV,Y1,Y2,DX,QIN,H,SKAP,XKAP
DOUBLE PRECISION XLD,DVDR(NRD)
PARAMETER (H = 0.5D0)
DOUBLE PRECISION GKUTTA
EXTERNAL GKUTTA
SKAP = -1.
DX = -DX
XM = 0.5D0*(1.D0 + (E - V(NTP2))/CSQ)
XMU = - (FLOAT(L)+0.5D0)**2 - 2.D0*XM*(V(NTP2)-E)*RAD(NTP2)**2
IF(XMU.LT.0.D0.AND.NTP2.LT.JRI)THEN
CALL DIFFRP(V,DX,JRI,DVDR,RAD)
BETA = -SKAP*(XKAP - DSQRT(-XMU) +
+ 0.5D0*( (FLOAT(L)+0.5D0)**2 - XM*DVDR(NTP2)*RAD(NTP2)**3 )/XMU)/
+ (2.D0*XM*RAD(NTP2))
ELSE
IF(L.LE.1)THEN
XLD = FLOAT(L)
ELSE
XLD = - FLOAT(L+1)
ENDIF
BETA = - SKAP*(XLD+XKAP+1)/(2.D0*XM*RAD(NTP2))
ENDIF
P(NTP2)= 1.D-6
Q(NTP2) = BETA*P(NTP2)
Y(1,NTP2)= P(NTP2)
Y(2,NTP2) = Q(NTP2)
DO J = 1,NTP2-NTP1
I = NTP2 - J + 1
Y1 = Y(1,I)
Y2 = Y(2,I)
A11 = -XKAP
TWOM = 1.D0 + (E - V(I))/CSQ
A12 = -SKAP*TWOM*RAD(I)
77
78
9
A22 = XKAP
A21 = -SKAP*(V(I)-E)*RAD(I)
XK(1,1) = DX*GKUTTA(1,Y1,Y2,A11,A12,A22,A21)
XK(2,1) = DX*GKUTTA(2,Y1,Y2,A11,A12,A22,A21)
TWOM = 1.D0 + (E - 0.5D0*(V(I)+V(I-1)))/CSQ
A12 = -SKAP*TWOM*RAD(I)*DEXP(DX/2.D0)
A21 = -SKAP*( VI(I-1) - E*RAD(I)*DEXP(DX/2.D0) )
XK(1,2)=DX*GKUTTA(1,(Y1+H*XK(1,1)),(Y2+H*XK(2,1)),A11,A12,A22,A21)
XK(2,2)=DX*GKUTTA(2,(Y1+H*XK(1,1)),(Y2+H*XK(2,1)),A11,A12,A22,A21)
XK(1,3)=DX*GKUTTA(1,(Y1+H*XK(1,2)),(Y2+H*XK(2,2)),A11,A12,A22,A21)
XK(2,3)=DX*GKUTTA(2,(Y1+H*XK(1,2)),(Y2+H*XK(2,2)),A11,A12,A22,A21)
TWOM = 1.D0 + (E - V(I-1))/CSQ
A12 = -SKAP*TWOM*RAD(I-1)
A21 = -SKAP*(V(I-1)-E)*RAD(I-1)
XK(1,4) = DX*GKUTTA(1,(Y1+XK(1,3)),(Y2+XK(2,3)),A11,A12,A22,A21)
XK(2,4) = DX*GKUTTA(2,(Y1+XK(1,3)),(Y2+XK(2,3)),A11,A12,A22,A21)
DO ID =1,2
Y(ID,I-1) = Y(ID,I) +
+ (XK(ID,1)+2.D0*(XK(ID,2)+XK(ID,3))+XK(ID,4))/6.D0
ENDDO
P(I-1) = Y(1,I-1)
Q(I-1) = Y(2,I-1)
ENDDO
Y1 = P(NTP1)
DO I = NTP1,NTP2
P(I) = P(I)/Y1
Q(I) = Q(I)/Y1
ENDDO
QIN = Q(NTP1)
IF(NTP2.LT.JRI)THEN
DO I = NTP2,JRI
P(I) = P(NTP2)*DEXP(-DSQRT(-E-0.25*E**2/CSQ)*(RAD(I)-RAD(NTP2)))
Q(I) = Q(NTP2)*DEXP(-DSQRT(-E-0.25*E**2/CSQ)*(RAD(I)-RAD(NTP2)))
ENDDO
ENDIF
DX = -DX
RETURN
END
SUBROUTINE RKOUT(E,RAD,V,VI,JRI,L,GAM,NODE,DX,P,Q,NTP1,
+ SKAP,XKAP,CSQ,TWOZ,QOUT,XLOGD,JNODE)
PROGRAMS
9.7
Self-consistent field atom local density approximation
INCLUDE "../INCLUDE/PARAMETERS"
INTEGER NODE,L
DOUBLE PRECISION DX,E,RATFG,S,H
PARAMETER (H = 0.5D0)
INTEGER JRI,NTP1,NTP1M,I,M,N,ID,NOLD,JNODE
DOUBLE PRECISION P(NRD),RAD(NRD)
DOUBLE PRECISION Y1,Y2,GAM,RNIE,S1,SKAP,XKAP,QOUT,XM,XLOGD
PARAMETER (N=2,M=4)
DOUBLE PRECISION Y(N,NRD),Q(NRD),V(NRD),VI(NRD),D1,D2,D3
DOUBLE PRECISION XK(N,M)
DOUBLE PRECISION A11,A12,A22,A21,X,EMV,CSQ,TWOM,TWOZ
DOUBLE PRECISION GKUTTA
EXTERNAL GKUTTA
C---------------------------------------C RUNGE-KUTTA
C---------------------------------------JNODE = 0
XLOGD = 0.
SKAP = -1.
P(1) = 1.D-6
Q(1) = - P(1)*SKAP*(GAM+XKAP)*CSQ/TWOZ
Y(1,1)= P(1)
Y(2,1) = Q(1)
NTP1M = NTP1 - 1
NODE = 0
EMV = E - V(1)
DO I = 1,NTP1M
Y1 = Y(1,I)
Y2 = Y(2,I)
A11 = -XKAP
TWOM = 1.D0 + (E - V(I))/CSQ
A12 = -SKAP*TWOM*RAD(I)
A22 = XKAP
A21 = -SKAP*(V(I)-E)*RAD(I)
XK(1,1) = DX*GKUTTA(1,Y1,Y2,A11,A12,A22,A21)
XK(2,1) = DX*GKUTTA(2,Y1,Y2,A11,A12,A22,A21)
TWOM = 1.D0 + (E - 0.5D0*(V(I)+V(I+1)))/CSQ
A12 = -SKAP*TWOM*RAD(I)*DEXP(DX/2.D0)
79
80
9
A21 = -SKAP*(VI(I) - E*RAD(I)*DEXP(DX/2.D0))
XK(1,2)=DX*GKUTTA(1,(Y1+H*XK(1,1)),(Y2+H*XK(2,1)),A11,A12,A22,A21)
XK(2,2)=DX*GKUTTA(2,(Y1+H*XK(1,1)),(Y2+H*XK(2,1)),A11,A12,A22,A21)
XK(1,3)=DX*GKUTTA(1,(Y1+H*XK(1,2)),(Y2+H*XK(2,2)),A11,A12,A22,A21)
XK(2,3)=DX*GKUTTA(2,(Y1+H*XK(1,2)),(Y2+H*XK(2,2)),A11,A12,A22,A21)
TWOM = 1.D0 + (E - V(I+1))/CSQ
A12 = -SKAP*TWOM*RAD(I+1)
A21 = -SKAP*(V(I+1)-E)*RAD(I+1)
XK(1,4) = DX*GKUTTA(1,(Y1+XK(1,3)),(Y2+XK(2,3)),A11,A12,A22,A21)
XK(2,4) = DX*GKUTTA(2,(Y1+XK(1,3)),(Y2+XK(2,3)),A11,A12,A22,A21)
DO ID = 1,2
Y(ID,I+1) = Y(ID,I) +
+ (XK(ID,1)+2.D0*(XK(ID,2)+XK(ID,3))+XK(ID,4))/6.D0
ENDDO
P(I+1) = Y(1,I+1)
Q(I+1) = Y(2,I+1)
NOLD = NODE
NODE = NODE + DABS(DSIGN(1.D0,P(I+1))-DSIGN(1.D0,P(I)))/2.D0
IF(NODE.GT.NOLD)THEN
JNODE = I+1
ENDIF
ENDDO
IF(NTP1.EQ.JRI)THEN
NTP1M = NTP1
DO I = 1,NTP1M
P(I) = P(I)/P(NTP1)
Q(I) = Q(I)/P(NTP1)
ENDDO
XM = 0.5D0*(1.D0 + (E - V(NTP1)/CSQ))
XLOGD = - XKAP - 1.D0 - 2.D0*XM*SKAP*RAD(NTP1)*Q(NTP1)/P(NTP1)
QOUT = Q(NTP1)/P(NTP1)
ELSE
IF(NTP1.EQ.JRI)THEN
NTP1M = JRI
ENDIF
DO I = 1,NTP1M
P(I) = P(I)/P(NTP1)
Q(I) = Q(I)/P(NTP1)
ENDDO
QOUT = Q(NTP1)/P(NTP1)
ENDIF
PROGRAMS
9.7
Self-consistent field atom local density approximation
RETURN
END
SUBROUTINE WAVEEQ (JRI,EIGSUM,NC,KE2,RKE,VKE2,CADOTP,VEISUM,EIG,
+ RNODE)
INCLUDE "../INCLUDE/PARAMETERS"
DOUBLE PRECISION E,ETOP,EBOT,EX,FN,TWOZ,R0,DEL,BETA,TCH,TOL3
DOUBLE PRECISION CSQ,EMV,EMVM1,DE,ERR,DFDE,E1,E2,ET,EB,FAC,Y2
DOUBLE PRECISION EIGSUM,XKAP,SKAP,EIG(LDIM),VE1,VE4,KE2,RKE,VKE2
DOUBLE PRECISION CADOTP,VEISUM,RNODE(LDIM)
PARAMETER (TOL1=1.D-3)
INTEGER L,KAPPA,NTP1,NTP2,NTP11(LDIM),NTP12(LDIM),NODES(LDIM)
INTEGER NODE,NODEC,JNODE
DOUBLE PRECISION P(NRD),Q(NRD),F(NRD),XLOGD,XLD(LDIM)
DOUBLE PRECISION ZN,XN(LDIM),XL(LDIM),XJ(LDIM),XZ(LDIM),XZS(LDIM)
INTEGER JTOT,JRI,I,J,NC(LDIM)
CHARACTER*4 NAM
CHARACTER*8 NLJ(LDIM)
DOUBLE PRECISION WFP(LDIM,NRD),WFQ(LDIM,NRD)
COMMON/WFN/WFP,WFQ
COMMON/SETUP1/ZN,XN,XL,XJ,XZ,XZS,NLJ,NAM,JTOT
DOUBLE PRECISION RAD(NRD),V(NRD),VR(NRD),VI(NRD),DVDR(NRD),DX
COMMON/RADIAL/RAD,V,VR,VI,DX
DOUBLE PRECISION RHOR2N(NRD),RHOR2O(NRD),RHOR2(NRD),VALRHO(NRD)
COMMON/CHARGE/RHOR2N,RHOR2O,RHOR2,VALRHO
DOUBLE PRECISION SIMF,ADLINT
IF(JRI.GT.NRD) WRITE(*,*) ’ JRI GREATER THAN DIMENSIONS’,JRI,NRD
IF(JRI.GT.NRD) STOP
TWOZ = 2.D0*ZN
CSQ = C*C
DO I = 1,NRD
81
82
9
RHOR2N (I) = 0.D0
VALRHO (I) = 0.D0
ENDDO
VEISUM = 0.D0
EIGSUM = 0.D0
KE2 = 0.D0
RKE = 0.D0
VKE2 = 0.D0
CADOTP = 0.D0
DO J = 1,JTOT
FN = XN(J)
L = INT(XL(J))
SKAP = -1.D0
XKAP = (XJ(J)+0.5D0)*2.D0*(XL(J)-XJ(J))
KAPPA = INT(XKAP)
CALL GETEIG(P,Q,TWOZ,FN,L,JRI,E,NODE,SKAP,XKAP,CSQ,XLOGD,
+ NTP1,NTP2,JNODE)
EIG(J) = E
XLD(J) = XLOGD
NTP11(J) = NTP1
NTP12(J) = NTP2
NODES(J) = NODE
RNODE(J) = 0.D0
IF(JNODE.NE.0)THEN
RNODE(J) = RAD(JNODE)
ENDIF
EIGSUM = EIGSUM + XZ(J)*E
C--------------------------------------------------------------C ACCUMULATE CHARGE DENSITY
DO I =1,JRI
RHOR2N (I) = RHOR2N(I) + (P(I)*P(I)+Q(I)*Q(I)/CSQ)*XZ(J)
F(I) = Q(I)*Q(I)*(1.D0 +2.D0*(E-V(I))/CSQ)
WFP(J,I) = P(I)
WFQ(J,I) = Q(I)
ENDDO
C LOOK AT RKE
DO I =1,JRI
F(I) = Q(I)*Q(I)
ENDDO
RKE = RKE + XZ(J)*SIMF(RAD,F,JRI,DX)
C LOOK AT CADOTP
DO I =1,JRI
F(I) = 2.D0*Q(I)*Q(I)*(1.D0 +(E-V(I))/CSQ)
ENDDO
CADOTP = CADOTP + XZ(J)*SIMF(RAD,F,JRI,DX) +
+ XZ(J)*(Q(1)*P(1) - Q(JRI)*P(JRI))
PROGRAMS
9.7
Self-consistent field atom local density approximation
ENDDO
DO J = 1,JTOT
WRITE(*,*) NLJ(J),’ EIGENVALUE= ’,EIG(J),’ OCCUPATION = ’,XZ(J)
ENDDO
END
SUBROUTINE COUL(JRI,TWOZ,VCT)
INCLUDE "../INCLUDE/PARAMETERS"
DOUBLE PRECISION RAD(NRD),V(NRD),VR(NRD),VI(NRD),DX
COMMON/RADIAL/RAD,V,VR,VI,DX
DOUBLE PRECISION RHOR2N(NRD),RHOR2O(NRD),RHOR2(NRD),VALRHO(NRD)
COMMON/CHARGE/RHOR2N,RHOR2O,RHOR2,VALRHO
DOUBLE PRECISION Y(NRD),Z(NRD),F(NRD),TWOZ,VCT
INTEGER JRI,I
DOUBLE PRECISION SIMF
CALL SIMINT(DX,RHOR2,Y,JRI,NRD)
DO I =1,JRI
F(I) = RAD(I)*RHOR2(I)
ENDDO
CALL SIMINT(DX,F,Z,JRI,NRD)
DO I =1,JRI
V(I) = (Z(I)/RAD(I) + Y(JRI) - Y(I))*2.D0 - TWOZ/RAD(I)
ENDDO
C----------------------------------------------------------C COULOMB CONTRIBUTION TO TOTAL ENERGY
DO I=1,JRI
F(I) = - 0.5D0*(V(I)+TWOZ/RAD(I))*RHOR2(I)
ENDDO
VCT = SIMF(RAD,F,JRI,DX)
C----------------------------------------------------------RETURN
END
SUBROUTINE INTERP(JRI)
INCLUDE "../INCLUDE/PARAMETERS"
INTEGER N,J,I,A,K,M
PARAMETER (N=9,M=4)
C N=2M+1
DOUBLE PRECISION RAD(NRD),V(NRD),VR(NRD),VI(NRD),DX
COMMON/RADIAL/RAD,V,VR,VI,DX
DOUBLE PRECISION X(N),F(N),XI,FI,ERR
C FIRST M POINTS
DO I = 1,N
83
84
9
F(I) = V(I)*RAD(I)
X(I) = RAD(I)
ENDDO
DO I = 1,M
XI = X(I)*DEXP(DX/2.D0)
CALL AITKEN(N,X,F,XI,FI,ERR)
VI(I) = FI
ENDDO
C FINAL M POINTS
DO I = 1,N
J = JRI-N+I
F(I) = V(J)*RAD(J)
X(I) = RAD(J)
ENDDO
DO I = N-M+1,N
XI = X(I)*DEXP(DX/2.D0)
CALL AITKEN(N,X,F,XI,FI,ERR)
J = JRI-N+I
VI(J) = FI
ENDDO
C INTERIOR POINTS
DO K = M+1,JRI-M
J = K - M - 1
DO I = 1,N
F(I) = V(J+I)*RAD(J+I)
X(I) = RAD(J+I)
ENDDO
XI = X(M+1)*DEXP(DX/2.D0)
CALL AITKEN(N,X,F,XI,FI,ERR)
VI(K) = FI
ENDDO
DO I=1,JRI
ENDDO
RETURN
END
SUBROUTINE LOCEX(DX,R,CHAR1,CHAR2,EPSXC,VXC1,VXC2,JRI,STONH)
IMPLICIT REAL*8 (A-H,O-Z)
INCLUDE "../INCLUDE/PARAMETERS"
DIMENSION CHAR1(NRD),CHAR2(NRD),EPSXC(NRD),VXC1(NRD),VXC2(NRD),
C
*
ARS(NRD),R(NRD),DC1(NRD),DC2(NRD),C1(NRD),STONH(NRD)
*
ARS(NRD),R(NRD),STONH(NRD)
C ......................................................................
D43=4.D0/3.D0
PROGRAMS
9.7
Self-consistent field atom local density approximation
AA=.5D0**0.333333333333D0
D13 = 1.D0/3.D0
DO I = 1,JRI
CHAR=CHAR1(I)+CHAR2(I)
IF(CHAR.LE.1.E-8) THEN
EPSXC(I)=0.D0
VXC1(I)=0.D0
VXC2(I)=0.D0
ARS(I)=0.D0
ELSE
RCE=R(I)*R(I)
RS1=(CHAR/(3.D0*RCE))**0.333333333333333D0
RS=1.D0/RS1
EPSXP=-.91633059D0/RS
X=CHAR1(I)/CHAR
IF(X.LT..000001D0) X=.000001D0
IF(X.GT..999999D0) X=.999999D0
FX=(X**D43+(1.D0-X)**D43-AA)/(1.D0-AA)
RSF=RS/75.D0
RSF2=RSF*RSF
RSF3=RSF2*RSF
RSP=RS/30.D0
RSP2=RSP*RSP
RSP3=RSP2*RSP
FCF=(1.D0+RSF3)*DLOG(1.E0+1.E0/RSF)+0.5E0*RSF-RSF2-1.D0/3.D0
FCP=(1.D0+RSP3)*DLOG(1.E0+1.E0/RSP)+0.5D0*RSP-RSP2-1.D0/3.D0
EPSCP=-.0504D0*FCP
EPSCF=-.0254D0*FCF
CNY=5.1297628D0*(EPSCF-EPSCP)
EPSXC(I)=EPSXP+EPSCP+FX*(CNY+4.D0/3.D0*EPSXP)/5.1297628D0
ARS(I)=-1.22177412D0/RS+CNY
STONH(I) = 2.D0*ARS(I)/(3.D0*CHAR)
BRS=-0.0504D0*DLOG(1.E0+30.E0/RS)-CNY
TRX1=(2.D0*DABS(CHAR1(I))/CHAR)**0.333333333333333D0
TRX2=(2.D0*DABS(CHAR2(I))/CHAR)**0.333333333333333D0
VXR1 = 0.0D0
VXR2 = 0.0D0
VXC1(I)=ARS(I)*TRX1+BRS + VXR1
VXC2(I)=ARS(I)*TRX2+BRS + VXR2
ENDIF
ENDDO
RETURN
END
SUBROUTINE TPOT(JRI,EXCT,VXCT,EXCTS,VXCTS,STONH)
INCLUDE "../INCLUDE/PARAMETERS"
DOUBLE PRECISION RAD(NRD),V(NRD),VR(NRD),VI(NRD),DX
COMMON/RADIAL/RAD,V,VR,VI,DX
85
86
9
DOUBLE PRECISION RHOR2N(NRD),RHOR2O(NRD),RHOR2(NRD),VALRHO(NRD)
COMMON/CHARGE/RHOR2N,RHOR2O,RHOR2,VALRHO
DOUBLE PRECISION BZ,VRS(NRD),RHOS(NRD)
COMMON/MAG1/BZ,RHOS,VRS
DOUBLE PRECISION
DOUBLE PRECISION
DOUBLE PRECISION
INTEGER JRI,I
DOUBLE PRECISION
Y(NRD),Z(NRD),F(NRD),STONH(NRD)
EPSXC(NRD),VXUP(NRD),VXDO(NRD),EXCT,VXCT
EXCTS,VXCTS
SIMF
C-----------------------------------------------------C POTENTIAL
DO I =1,JRI
Y(I) = (RHOR2(I) + RHOS(I))/2.D0
Z(I) = (RHOR2(I) - RHOS(I))/2.D0
ENDDO
CALL LOCEX(DX,RAD,Y,Z,EPSXC,VXUP,VXDO,JRI,STONH)
DO I =1,JRI
V(I) = V(I) + (VXUP(I)+VXDO(I))/2.D0
VRS(I) = (VXUP(I)-VXDO(I))/2.D0
ENDDO
C----------------------------------------------------C EXCHANGE-CORRELATION CONTRIBUTION TO THE TOTAL ENERGY
DO I =1,JRI
F(I) = EPSXC(I)*RHOR2(I)
ENDDO
EXCT = SIMF(RAD,F,JRI,DX)
C----------------------------------------------------C VXCT FOR VIRIAL
DO I =1,JRI
F(I) =
VXUP(I)*Y(I) + VXDO(I)*Z(I)
ENDDO
VXCT = SIMF(RAD,F,JRI,DX)
C----------------------------------------------------C SURFACE TERMS FOR VIRIAL
EXCTS = RAD(JRI)*EPSXC(JRI)*VALRHO(JRI)
VXCTS = RAD(JRI)*VXUP(JRI)*VALRHO(JRI)/2.D0 +
+ VXDO(JRI)*VALRHO(JRI)/2.D0
RETURN
END
SUBROUTINE AITKEN(N,X,F,XI,FI,ERR)
INTEGER I,J,N,NMAX
PARAMETER (NMAX=15)
DOUBLE PRECISION X(N),F(N),FH(N),XI,FI,ERR,X1,X2,F1,F2
IF(N.GT.NMAX) THEN
PRINT*,’DIMENSIONS TOO SMALL IN AITKEN’
STOP
PROGRAMS
9.7
Self-consistent field atom local density approximation
ENDIF
DO I = 1,N
FH(I) = F(I)
ENDDO
DO I = 1,N-1
DO J = 1,N-I
X1 = X(J)
X2 = X(J+I)
F1 = FH(J)
F2 = FH(J+1)
FH(J) = (XI-X1)/(X2-X1)*F2 + (XI-X2)/(X1-X2)*F1
ENDDO
ENDDO
FI = FH(1)
ERR = (DABS(FI-F1) + DABS(FI-F2))/2.D0
RETURN
END
SUBROUTINE DIFFRP(F,DX,N,FP,RAD)
C-------------------------------------------------------------------C
FIVE-POINT DIFFERENTIATION FORMULA
C
INPUT F - CHANGED TO R*F TO DIFFERENTIATE
C
LOG MESH - USES D(RF)/DX = RF + R^2 DF/DR
C
DF/DR = (D(RF)/DX - RF)/R^2
C
OUTPUT FP = DF/DR
C-------------------------------------------------------------------C
.. PARAMETERS ..
INCLUDE "../INCLUDE/PARAMETERS"
DOUBLE PRECISION RAD(NRD),F(NRD),FP(NRD),RF(NRD),DX
INTEGER N,I
DO I = 1,N
RF(I) = F(I)*RAD(I)
ENDDO
FP(1) = RF(1)
FP(2) = (
+ 6.D0*RF(3)+ 20.D0*RF(5)/3.D0 + 6.D0*RF(7)/5.D0
+ -49.D0*RF(2)/20.D0 - 15.D0*RF(4)/2.D0 - 15.D0*RF(6)/4.D0
+ -RF(8)/6.D0
+ )/DX
DO I = 3,N-2
87
88
9
FP(I) = ((RF(I-2)+8.D0*RF(I+1))- (8.D0*RF(I-1)+RF(I+2)))/12.D0/DX
ENDDO
FP(N-1) = (RF(N)-RF(N-2))/2.D0/DX
FP(N) = (RF(N-2)/2.D0-2.D0*RF(N-1)+1.5D0*RF(N))/DX
DO I = 1,N
FP(I) = (FP(I)-RF(I))/RAD(I)**2
ENDDO
RETURN
END
DOUBLE PRECISION FUNCTION SIMF(X,YT,N,DX)
C SIMPSON’S RULE FOR LOG MESH
C----------------------------------------------------------INCLUDE "../INCLUDE/PARAMETERS"
INTEGER N,I
DOUBLE PRECISION Y(NRD),X(NRD),YT(NRD),DX,SIMF
C----------------------------------------------------------C
MULTIPLIES BY X TO ACCOUNT FOR LOG MESH
DO I = 1,N
Y(I) = YT(I)*X(I)
ENDDO
C----------------------------------------------------------SIMF = 0.D0
DO I = 2,N-1,2
SIMF = SIMF + (Y(I-1)+4.D0*Y(I)+Y(I+1))
ENDDO
SIMF = DX*SIMF/3.D0
IF(MOD(N,2).EQ.0)THEN
SIMF = SIMF + DX*(5.D0*Y(N)+8.D0*Y(N-1)-Y(N-2))/12.D0
ENDIF
RETURN
END
SUBROUTINE SIMINT(H,Y,Z,N,NDIM)
INTEGER N,NDIM,I
DOUBLE PRECISION Y(NDIM),Z(NDIM)
DOUBLE PRECISION H,SUM11,SUM21,H3,SUM12,SUM22
C
H3 = H/3.D0
IF (N.LT.6) THEN
WRITE(*,*) ’SIMINT - N TOO SMALL’,N
STOP
ELSE
SUM12 = H3* (Y(1)+4.D0*Y(2)+Y(3))
SUM11 = SUM12 + H3* (Y(3)+4.D0*Y(4)+Y(5))
SUM21 = H3* (Y(1)+31.D0*(Y(2)+Y(5))/8.D0 +
+ 21.D0*(Y(3)+Y(4))/8.D0+Y(6))
SUM22 = SUM21 - H3* (Y(4)+4.D0*Y(5)+Y(6))
PROGRAMS
9.7
Self-consistent field atom local density approximation
Z(1) = 0.D0
Z(2) = SUM22 - H3* (Y(2)+4.D0*Y(3)+Y(4))
Z(3) = SUM12
Z(4) = SUM22
IF (N-6.GT.0) THEN
C-----------------------------------------------------------C
INTEGRATION LOOP
C----------------------------------------------------------DO 10 I = 7,N,2
SUM12 = SUM11
SUM22 = SUM21
SUM11 = SUM12 + H3* (Y(I-2)+4.D0*Y(I-1)+Y(I))
Z(I-2) = SUM12
IF (I-N.GE.0) THEN
Z(N-1) = SUM22
Z(N) = SUM11
RETURN
ELSE
SUM21 = SUM22 + H3* (Y(I-1)+4.D0*Y(I)+Y(I+1))
Z(I-1) = SUM22
END IF
10
CONTINUE
END IF
Z(N-1) = SUM11
Z(N) = SUM21
END IF
END
89
90
REFERENCES
References
Brooks, M. & Johansson, B. (1983), �Exchange integral matrices and cohesive energies of transition metal atoms’,
J. Phys. F: Metal Phys. 13, L197–L202.
Das, M. P., Ramana, M. V., & Rajagopal, A. K. (1980), �Self-consistent relativistic density-functional theory:
Application to neutral uranium atom and some ions of lithium isoelectronic sequence’, Physical Review A
22, 9–13.
Desclaux, J. (1973), Atomic Data and Nuclear Data 12, 311.
Hohenberg, P. & Kohn, W. (1964), �Inhomogeneous electron gas’, Physical Review 136, 864–871.
Kagawa, T. (1975), �Relativistic hartree-fock-roothaan theory for open-shell atoms’, Physical Review A 12, 2245
– 2256.
Kohn, W. & Sham, L. J. (1965), �Self-consistent equations including exchange and correlation effects’, Physical
Review 140, A1133–A1138.
O.K.Andersen, H.L.Skriver, H.Nohl & B.Johansson (1979), Pure. Appl. Chem. 52, 93.
Pettifor, D. (1995), Bonding and Structure of Molecules and Solids, Clarendon Press, Oxford.
Snow, E. C., Canfield, J. M. & Waber, J. T. (n.d.).
Streater, R. & Wightman, A. (1978), PCT, Spin and Statistics and All That, Benjamin/Cummings, Reading,
MA.
von Barth, U. & Hedin, L. (1972), J. Phys. C: Solid State Phys. 5, 1629.
Документ
Категория
Без категории
Просмотров
12
Размер файла
7 545 Кб
Теги
1/--страниц
Пожаловаться на содержимое документа