Band structure calculation using KKR-Green’s function methodby Kazunori Sato (June 1998)
modiﬁed by Masako Ogura and Hisazumi Akai (August 2008)

1

Introduction

In this document, how to use Machikaneyama2002 (AkaiKKR), a KKR-CPA-LDA package, isbrieﬂy guided. The KKR method is one of the methods of electronic structure calculationand is also called“Green’s function method.” KKR indicates the initials of Korringa, Kohn,
and Rostoker, who invented this method. CPA means the coherent potential approximation,
which can deal with random systems. Green’s function might be diﬃcult to understand if you
are trained with ordinary band structure calculations, which solve an eigenvalue problem by
diagonalization. However, you can apply KKR-CPA to much wider range of situation thanordinary band structure calculations can do. For example, it can deal with ﬁnite temperaturemagnetism and partial disorder systems since CPA can treat not only periodic systems but also
random systems.

2

Get started

The KKR band structure code has been developed by H. Akai since 1979. The KKR package isa set of subroutines written in FORTRAN77. It is not so diﬃcult for you to perform some small
calculations using this package even if you are not familiar with the computational physics andnumerical techniques used in the method.

2.1

Get the program

Usually, the downloaded KKR program package is archived and compressed. Thus compressedﬁle is named as ’cpa2002v00x.tar.gz’. First you put ’cpa2002v00x.tar.gz’ under your home
directly. Next, you must decompress it and get the original form. The following commands
should be executed.
~> tar xvfz cpa2002v00x.tar.gz
Now, you ﬁnd the following ﬁles in the directory, /cpa2002v00x/.
~>cd cpa2002v00x~/cpa2002v00x>ls
data
gpdos.f
gpd
inmakefile
readmesourceutil

Here, the ﬁle ’makeﬁle’ is used to compile the sources. New directories, ‘source’, ’in’, ’out’
and ’data’, are created. Under the directory ’source’ you ﬁnd all of the source codes, and ’in’,’out’ and ’data’ may be used to save input, output and data ﬁles, respectively. You can delete’cpa2002v00x.tar.gz’ which is not needed anymore.

2.2

Compile the program

Since the package contains many subroutines, it is convenient to use the ’make’ command. The’make’ command refers to the ’makeﬁle’ and compiles the ﬁles whose names are written in the’makeﬁle’. Moreover, the ’make’ command ﬁnds out the ﬁles that have been modiﬁed since thelast time, and compiles those ﬁles.
1^/cpa2002v00x> make
Just key in, and you will get an executable ﬁle ’specx’.

2.3

Execute the program

To start a calculation, an input ﬁle in which lattice parameters and some numerical parameters
are written is needed. How to write the input ﬁle is explained in the next section. If you want
to execute a KKR calculation with your input ﬁle ’input’ and save a result in a ﬁle ’output’, key
in the following.
~/cpa2002v006x> specx < in/input > out/output &The last character ’&’ means that this job should be carried out as a background job.

3

Execution of the calculation

3.1

Input File

The lattice parameters, such as the Bravais type, the lattice constants and atomic numbers, andso on, are written in an input ﬁle and it is read in by the program when the kkr is executed. As
an example, the input ﬁle used in the calculation of fcc NiFe alloy is the following.
c----------------------- input data ------------------------c
go
file
go
data/fenic ----------------------------------------------------------c
brvtypac/ab/a alpha beta gamma
bcc
5.3056 ,,,,,,c ----------------------------------------------------------c
edelt
ewidthreltypsdftypmagtyprecord
0.001
1.2
nrl
mjw
mag
init
c ----------------------------------------------------------c
outtypbzqltymaxitrpmix
update480
0.024
c ----------------------------------------------------------cntyp
1c ----------------------------------------------------------c
type
ncmp
rmt
field
l_max
anclr
conc
FeNi
2
0
0
2
26
60
28
40
c ----------------------------------------------------------cnatm
1c ----------------------------------------------------------c
atmicx
atmtyp000FeNic -----------------------------------------------------------The lines which begin with the letter ’c’ is neglected. Names of the ﬁelds are also commented
in the sample input ﬁle. The delimiter is a space ’ ’ or a comma ’,’. In table 1, what is written
in each ﬁeld is summarized.
In the following, some comments and notes are given.
• [go] In ’dsp’ and ’dos’ modes, values of record, outtyp and maxitr are automatically set
to be ‘2nd’, ‘quit’ and ‘1’, respectively.
2

Table 1: Parameters written in the input ﬁle.
Field
go

ﬁle
brvtyp

Value
go
dos
dsp
spc
ﬁle name
fcc
bcc
hcp(hex)
sc
bct
st
fco
bco
bso
so
bsm
sm
trc
rhb(trg)
fct

a
c/a
b/a
alpha
beta
gamma
edelt
ewidthreltypsdftypmagtyprecordouttypbzqltymaxitr
pmix
ntyp
type
ncmp
rmt
ﬁeld
l max
anclr
conc
natm
atmicx
atmtyp∼0.001—nrl
sra
mjw, vbh, vwn
mag
nmag-mag(rvrs)
kick
init
2nd
1st
updatequit
0, 1, 2,...
t, s, m, l, h, u
0.01-0.03
type namecoordinates
type name

Meaning
perform a band structure calculation.
calculate a density of states.
display a previous result.
calculate Bloch spectral function (dispersion relation).
ﬁle name in which a potential data is written.
which type of the Bravais lattices.
face centered cubic
body centered cubic
hexagonal close packed
simple cubic
body centered tetragonal
simple tetragonal
face centered orthorhombic
body centered orthorhombic
base centered orthorhombic
simple orthorhombic
base centered monoclinic
simple monoclinic
triclinic
rhombohedral (trigonal)
face centered tetragonallattice constant in atomic unit.
c/a ratios of lattice constants.
b/a ratios of lattice constants.
α in degrees.
β in degrees.
γ in degrees.
Imaginary part at the Fermi level in Ry (see Fig. 1).
Width of the energy contour in Ry (see Fig. 1).
no relativistic treatment.
scaler relativistic approximation.
which exchange correlation potential is used.
magnetic.
non-magnetic.
change the sign of the magnetization.
transfer to the magnetic state artiﬁcially.
initially start a calculation.
continues the last calculation.
continues the second last calculation.
potential data is updated.
potential data is not updated.
The bigger, the ﬁner mesh in the Brillouin zone.
see nfqlty.f
the maximum number of the iteration loop.
a parameter used in mixing V in and V out .
how many types of atom exist in a unit cell.
names of the respective types of site.
how many components exist in the site.
muﬃn tin radius in a.
external magnetic ﬁeld at each site in Ry.
the maximum angular momentum taken into account.
atomic number.
concentration of the components at the site.
number of atoms in a unit cell.
where each atom is (in a).
which type of the site at the lattice point.

3

Im E

energy contour C

core state

valence stateedelt

ewidthRe E

Figure 1: The energy contour.
• [ﬁle] You can specify a ﬁle name of a potential data using a relative pass.
• [brvtyp] If c/a, b/a, α, β and γ are evident from ’brvtyp’ , you can omit to write them.
• [record] The previous potential data is stored together with the latest potential datain the data, in case that the latest calculation terns out rubbish. ’2nd’ means that thecalculation starts with the latest data ﬁle, and ’1st’ means that the calculation starts withthe data of the last but one.
• [rmt] If the muﬃn-tin spheres whose radii are given by ’rmt’ conﬂict each other theirradii are reduced. If ’rmt’ = 0, muﬃn tin radii are set so that their ratio is equal to theratio of the atomic radii.
• [anclr] ’anclr’=0 means that a vacancy is put on the lattice point.
• [conc] Concentrations are normalized in each site. It is not necessary that the sum ofthem equal to 1.
• [atmicx] The lattice point at which an atom is put is speciﬁed with the Cartesian coordinate (e.g., written as 0.5, 0.5, 0.5) or with the primitive vector (e.g., written as 0.5a, 0.5b,0.5c). They can have a fractional form such as 1/3, 2/3a, and 1d0/3d0x.
• If the input data for several diﬀerent systems are contained in a single ﬁle, they are
performed sequentially.

3.2

Output File

A result of a calculation is shown on a display or is redirected to an output ﬁle. As an example,
the output of the calculation for NiFe alloy is shown, for which the input ﬁle has already beenshown.
2-Feb-2005
meshr mse
400
35

ng

mxl
15

3

data read in
go=go
file=data/feni
brvtyp=bcc a= 5.30560 c/a=1.00000 b/a=1.00000
alpha= 0.0 beta= 0.0 gamma= 0.0
edelt= 0.0010 ewidth= 1.200 reltyp=nrl
sdftyp=mjw
magtyp=mag
record=init outtyp=update bzqlty=4
maxitr= 80 pmix=0.02400
ntyp= 1 natm= 1 ncmpx= 2
complex energy mesh1( -1.2000, 0.0000)
4( -1.1970, 0.0121)
7( -1.1710, 0.0368)
10( -1.0562, 0.0779)

2( -1.1998, 0.0030)
5( -1.1928, 0.0186)
8( -1.1478, 0.0489)
11( -0.9792, 0.0930)

3( -1.1990, 0.0070)
6( -1.1849, 0.0267)
9( -1.1111, 0.0629)
12( -0.8787, 0.1063)

413( -0.7577, 0.1158)
16( -0.3646, 0.1104)
19( -0.1156, 0.0708)
22( -0.0276, 0.0360)
25( -0.0057, 0.0165)
28( -0.0011, 0.0073)
31( -0.0002, 0.0032)
34( 0.0000, 0.0014)

14( -0.6240, 0.1199)
17( -0.2591, 0.0988)
20( -0.0735, 0.0575)
23( -0.0165, 0.0280)
26( -0.0033, 0.0126)
29( -0.0006, 0.0055)
32( -0.0001, 0.0024)
35( 0.0000, 0.0010)

15( -0.4891, 0.1179)
18( -0.1764, 0.0850)
21( -0.0455, 0.0458)
24( -0.0098, 0.0216)
27( -0.0019, 0.0096)
30( -0.0004, 0.0042)
33( -0.0001, 0.0018)

file to be accessed=data/feni
createdlattice constantbravais=bcc
a= 5.30560
c/a=1.0000
b/a=1.0000
alpha= 90.00
beta= 90.00
gamma= 90.00
primitive translation vectors
a=( -0.50000 0.50000 0.50000)
b=( 0.50000 -0.50000 0.50000)
c=( 0.50000 0.50000 -0.50000)
type of site
type=FeNi

rmt=0.43301 field= 0.000
lmxtyp= 2
component= 1 anclr= 26. conc= 0.6000
component= 2 anclr= 28. conc= 0.4000

atoms in the unit cell
position=
0.00000000

0.00000000

0.00000000

type=FeNi***msg in spmain...new ew, ez generated
ew=
0.09998 ez=
0.80300
preta= 0.35542 eta= 0.35542
isymop= 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
last= 243

np= 15nt= 141

nrpt= 169

nk=

29

itr= 1rms error = -1.872
itr= 2rms error = -2.390
itr= 3rms error = -2.879
itr= 4rms error = -3.357
itr= 5rms error = -3.890
itr= 6rms error = -3.879
itr= 7rms error = -4.878
itr= 8rms error = -4.850
itr= 9rms error = -5.574
itr= 10rms error = -6.235
interval= 10
cpu time=
0.00 secnl
cnf
energy----------------------------------1s
2.000-508.5203
2s
2.000-59.20742p
6.000-51.18063s
2.000-6.80273p
6.000-4.45633d
6.000-0.66964s
2.000-0.4930
itr= 1rms error = -1.724
itr= 2rms error = -2.239
itr= 3rms error = -2.657
itr= 4rms error = -2.967
itr= 5rms error = -3.645
itr= 6rms error = -3.695
itr= 7rms error = -4.268
itr= 8rms error = -5.061
itr= 9rms error = -5.356
itr= 10rms error = -5.470
itr= 11rms error = -6.321
interval= 11
cpu time=
0.01 secnl
cnf
energy-----------------------------------

5

nd=

1

1s
2s
2p
3s
3p
3d
4s

2.000
2.000
6.000
2.000
6.000
8.000
2.000

-595.8142-70.7050-61.8163-7.9863-5.2730-0.7806-0.5252

record 1 will be overlaied by input andrecord 2 will be replaced by new output.
core configuration for Z= 26
state 1s 2s 2p 3s 3p 3d 4s 4p 4d 5s 5p 4f 5d 6s 6p 5f 6d 7sup
1 1 3 1 3 0 0 0 0 0 0 0 0 0 0 0 0 0
down
1 1 3 1 3 0 0 0 0 0 0 0 0 0 0 0 0 0
core configuration for Z= 28
state 1s 2s 2p 3s 3p 3d 4s 4p 4d 5s 5p 4f 5d 6s 6p 5f 6d 7sup
1 1 3 1 3 0 0 0 0 0 0 0 0 0 0 0 0 0
down
1 1 3 1 3 0 0 0 0 0 0 0 0 0 0 0 0 0***** self-consistent iteration starts *****
FeNi
itr= 1 neu= -1.4601 moment= 0.0905 te= -2718.3225808 err= 0.119
itr= 2 neu= -0.8991 moment= 0.1151 te= -2718.3221114 err= -0.277
itr= 3 neu= -0.2974 moment= 0.1168 te= -2718.3157206 err= -0.149
itr= 4 neu= 0.1787 moment= 0.1296 te= -2718.3197036 err= -0.627
itr= 5 neu= 0.5050 moment= 0.1933 te= -2718.3214574 err= -0.407
itr= 6 neu= 0.6526 moment= 0.2873 te= -2718.3261849 err= -0.483
itr= 7 neu= 0.5415 moment= 0.4203 te= -2718.3258145 err= -0.530
itr= 8 neu= 0.3129 moment= 0.5384 te= -2718.3258235 err= -1.056
itr= 9 neu= 0.1143 moment= 0.7051 te= -2718.3265769 err= -0.801
itr= 10 neu= -0.0203 moment= 0.8951 te= -2718.3295190 err= -0.873
itr= 11 neu= -0.1221 moment= 1.0540 te= -2718.3323905 err= -0.967
itr= 12 neu= -0.1746 moment= 1.2526 te= -2718.3352375 err= -0.844
itr= 13 neu= -0.1626 moment= 1.4605 te= -2718.3369449 err= -0.608
itr= 14 neu= -0.1454 moment= 1.6426 te= -2718.3395758 err= -0.754
itr= 15 neu= -0.1518 moment= 1.6880 te= -2718.3406928 err= -0.650
itr= 16 neu= -0.1663 moment= 1.7401 te= -2718.3410295 err= -0.607
itr= 17 neu= -0.1522 moment= 1.7988 te= -2718.3415185 err= -1.255
itr= 18 neu= -0.1303 moment= 1.7861 te= -2718.3414079 err= -1.326
itr= 19 neu= -0.0791 moment= 1.7878 te= -2718.3410749 err= -1.206
itr= 20 neu= -0.0246 moment= 1.7774 te= -2718.3410051 err= -1.553
itr= 21 neu= 0.0213 moment= 1.7734 te= -2718.3409834 err= -1.224
itr= 22 neu= 0.0494 moment= 1.7695 te= -2718.3410544 err= -1.282
itr= 23 neu= 0.0583 moment= 1.7667 te= -2718.3410975 err= -1.829
itr= 24 neu= 0.0518 moment= 1.7646 te= -2718.3410481 err= -1.769
itr= 25 neu= 0.0361 moment= 1.7609 te= -2718.3409914 err= -1.879
itr= 26 neu= 0.0172 moment= 1.7587 te= -2718.3409593 err= -2.014
itr= 27 neu= 0.0005 moment= 1.7597 te= -2718.3409526 err= -1.849
itr= 28 neu= -0.0111 moment= 1.7612 te= -2718.3409625 err= -2.108
itr= 29 neu= -0.0169 moment= 1.7625 te= -2718.3409716 err= -2.097
itr= 30 neu= -0.0176 moment= 1.7634 te= -2718.3409749 err= -2.114
itr= 31 neu= -0.0147 moment= 1.7639 te= -2718.3409743 err= -2.287
itr= 32 neu= -0.0102 moment= 1.7640 te= -2718.3409723 err= -2.403
itr= 33 neu= -0.0055 moment= 1.7640 te= -2718.3409707 err= -2.315
itr= 34 neu= -0.0017 moment= 1.7638 te= -2718.3409703 err= -2.499
itr= 35 neu= 0.0009 moment= 1.7637 te= -2718.3409703 err= -2.727
itr= 36 neu= 0.0024 moment= 1.7635 te= -2718.3409704 err= -2.832
itr= 37 neu= 0.0029 moment= 1.7633 te= -2718.3409703 err= -2.844
itr= 38 neu= 0.0026 moment= 1.7632 te= -2718.3409701 err= -2.749
itr= 39 neu= 0.0019 moment= 1.7631 te= -2718.3409699 err= -2.813
itr= 40 neu= 0.0010 moment= 1.7631 te= -2718.3409698 err= -2.999
itr= 41 neu= 0.0002 moment= 1.7631 te= -2718.3409698 err= -3.004
interval= 41
cpu time=
17.79 sec
sdftyp=mjw
reltyp=nrl
dmpc=0.024
FeNi
itr= 41 neu 0.0002 chr,spn 8.8000 1.7631 intc,ints 0.9632 -0.0234
rms err= -3.164 -3.274 -3.004 -3.177
ef=
0.7260630
0.7343370 def=
2.6962356 12.2016065
total energy= -2718.3409698

6

*** type-FeNi
Fe (z= 26.0) ***core charge in the muffin-tin sphere =17.9787766
valence charge in the cell (spin up ) =0.19826(s)0.19326(p)4.36075(d)valence charge in the cell (spin down) =0.19758(s)0.20985(p)1.84235(d)total charge= 24.98082
valence charge (up/down)=
4.75227
2.24978
spin moment=
2.50249 orbital moment=
0.00000
core level (spin up )-507.1266873 Ry(1s)-5.5302084 Ry(3s)core level (spin down)-507.1176964 Ry(1s)-5.3171341 Ry(3s)-57.8685202 Ry(2s)-3.1878588 Ry(3p)-49.8296484 Ry(2p)-57.7567489 Ry(2s)-2.9808461 Ry(3p)-49.7425563 Ry(2p)hyperfine field of FeNi
-239.803 KG (core= -250.980 KG valence=
11.177 KG )core contribution-21.109 KG(1s)-548.149 KG(2s)318.277 KG(3s)charge density at the nucleus
11820.1880 (core=
11814.5387 valence=core contribution
10701.4478(1s)972.7303(2s)5.6493 )140.3607(3s)*** type-FeNi
Ni (z= 28.0) ***core charge in the muffin-tin sphere =17.9914848
valence charge in the cell (spin up ) =0.21450(s)0.20389(p)4.50139(d)valence charge in the cell (spin down) =0.22638(s)0.23660(p)3.74662(d)total charge= 27.12086
valence charge (up/down)=
4.91978
4.20959
spin moment=
0.71019 orbital moment=
0.00000
core level (spin up )-594.3771108 Ry(1s)-6.5907736 Ry(3s)core level (spin down)-594.3686783 Ry(1s)-6.5162732 Ry(3s)-69.2884479 Ry(2s)-3.8813868 Ry(3p)-60.3956628 Ry(2p)-69.2422357 Ry(2s)-3.8088998 Ry(3p)-60.3587004 Ry(2p)hyperfine field of FeNi
-166.611 KG (core=-77.215 KG valence=-89.396 KG )core contribution-9.105 KG(1s)-204.777 KG(2s)136.667 KG(3s)charge density at the nucleus
14827.9567 (core=14822.0421 valence=core contribution
13396.9377(1s)1243.1999(2s)cpu used

5.9145 )181.9046(3s)17.94 sec

First, a date, ’meshr’, ’mse’, ’ng’ and ’mxl’ are output. The next block shows meshes onan energy contour and the contents of the input ﬁle. If muﬃn-tin radii are modiﬁed, the new
radii are shown. The next information shows how many lattice points (nrpt), reciprocal lattice
points (ngpt) and k-points (nk) are needed. When ’record=init’, the LDA calculation on each
atom is performed and an initial potential is generated. Then, conﬁgurations of core states arelisted. The core conﬁgurations are set up in subroutine ’corcnf.f’. It is noticed that the corestates must not come into the energy contour to avoid double counting of the core electrons.
After these output of initial conditions, the record of the self-consistent iteration is output.
For each iteration loop, charge neutrality, spin moment, total energy and error (the average
diﬀerence between input and output potentials) are shown to conﬁrm a convergence. If theerror becomes under a tolerance which is set in the subroutine ’specx.f’, the iteration loop stops,
and a result is printed out. If the number of loops exceeds ’maxitr’, the iteration loop also stops.
Lastly, from the obtained electronic structure, the hyperﬁne ﬁeld, and so on, are calculated.

7

Table 2: Parameters used in ’specx.f’.
Parameter
natmmxncnpmxmsizemx
mxlmx
nk1mx,nk3mx
msex
ngmx
nrpmx
ngpmx
msr
mse0
tol

3.3

Meaning
maximum number of atoms in a unit cell.
maximum number of types of atoms in a unit cell.
∑natm
maximum size of the KKR matrix. ≥ i=1 (lmax,i + 1)2
l ≥mxlmx are truncated.
nk1mx+nk3mx is maximum number of k-points in the Brillouin zone.
mesh points on the energy contour for go=dos, spc.
The chebyshev expansion is performed up to ngmx-th order.
maximum number of lattice points used in the Ewalt’s sum.
maximum number of reciprocal lattice points used in the Ewalt’s sum.
radial mesh points.
(data) mesh points on the energy contour for go=go.
(data) tolerance of convergence

Parameters which deﬁne the size of arrays

By the way, size of an array used in the KKR program depends on a system to be calculated.
It is necessary to reserve enough area. The size of the arrays in the KKR program is declaredin the subroutine ’specx.f’. The parameters used in the ’specx.f’ are summarized in table 2.
If unsuitable parameters are used in the input ﬁle or specx.f, the program will stop with anerror message. Then, you can recognize what is wrong.

4

Exercises
1. Calculate the electronic structure of bcc Fe. Use the lattice constant a = 5.27 bohr andtake magnetism into account. Does the obtained result show that Fe is ferromagnetic ?
How large is the magnetic moment ?
2. Calculate the density of states of bcc Fe.
3. Estimate the lattice constant of Fe. First, calculate the electronic structure for various
lattice constants. Then compare the total energy. The lattice constant which gives thelowest energy is realizable.
4. Calculate the electronic structure of fcc Ni. Use the lattice constant a = 6.60 bohr. Checkthe magnetic moment, density of states, lattice constant and so on.
5. (Impurity problem) Calculate bcc Fe with a single impurity which substitute a Fe atom.
Plot the calculated hyperﬁne ﬁeld as a function of the impurity’s atomic number. Discuss
the behavior of the hyperﬁne ﬁeld.
6. (Random alloy) Calculate fcc Ni-Fe random alloy. Use the lattice constant of Ni. Discuss
the density of states. What happens when the concentration of Fe is high ?
7. (Random alloy) Calculate various kinds of 3d transition metal alloys (Change combination of elements, concentration and structure). Plot the obtained magnetic moments asa function of the averaged atomic number. Thus obtained curve is called Slater-Paulingcurve.
8. (Local moment disorder) Calculate the electronic structure of bcc Fe in the local
moment disorder (LMD) state. We can treat the LMD state as a random alloy of Fewith up and down spin moments. Prepare the potential data ﬁle using /util/fmg.f8

9. (Local moment disorder) In the mean-ﬁeld theory, the Curie temperature TC is evaluated as follows:
kB TC =

2
(ELMD − EF ) ,3

where kB is Boltsmann constant, EF and ELMD are total energy in the ferromagnetic and
LMD states, respectively. Calculate Curie temperature from the obtained total energy.
Use the values, kB = 6.3336 × 10−6 Ry/K, if you need. Compare the results with theexperimental value, TC (exp.) = 1043 K.
10. (Local moment disorder) Calculate the electronic structure of fcc Ni in the LMD state.
Does the calculation works ? If not, discuss the reason.
11. Calculate the electronic structure of zinc-blende GaAs. The zinc-blende structure is constructed by the fcc structure. How is the input ﬁle ?
12. Put some defects in zinc-blende GaAs. Discuss the results.
13. Put some donors or accepters in zinc-blende GaAs. Discuss the results.
14. Put some magnetic impurities in zinc-blende GaAs. Discuss the results.

Figure 2: Zinc-blende, half-Heusler and full-Heusler structures.

9More in Details
5

Density Functional Theory

5.1

Hohenberg-Kohn Theorem

The Hohenberg-Kohn theorem on which the density functional theory is established is explained
according to the original paper.
An N -electron system under an external potential v(r) is described with the following Hamiltonian H.
H = T + V + U,

(1)

where T is the kinematic energy, V is the potential energy and U is the Coulomb energy between
electrons, i.e.,
T

=

∑

−∇2i

(2)

v(ri )(3)

2,|r − rj |
i<j i

(4)

i

V

=

∑
i

U

=

∑

where the atomic unit (h̄2 = 1, e2 = 2, me = 1/2) is used. The N -electron wave function Φ
satisﬁes the Schrödinger equation written as
HΦ = EΦ.

(5)

In the framework of the quantum theory, the ground state of the system, Ψ, under the given
external potential v is obtained by the variational principle
HΨ = E0 Ψ(6)

E0 [Ψ] = Min〈Ψ|H|Ψ〉(7)〈Ψ|Ψ〉 = N.

(8)Once the N -electron wave function Ψ is obtained, the density distribution ρ and any expectation
values are obtained from Ψ. That is Ψ, ρ and expectation energy E0 are the functionals of v.
Ψ = Ψ[v](9)ρ = ρ[v](10)E0 = E0 [v](11)

In this procedure, we must deal with Ψ which has 3N degrees of freedom and must take thevariation of the total energy functional with respect to Ψ.
If we can regard the density distribution ρ as a basic function in spite of Ψ, it will make theproblem easier because ρ has only 3 degrees of freedom. It is shown that v is a functional of ρ byTable 3: atomic unit
lengthenergyunith̄2= 1(bohr) = 0.529(Å)
me e2
4
me e
= 1(Ry) = 13.6(eV)
2h̄2

10

the reductio ad absurdum. If the same ρ gives two diﬀerent v and v ′ , we get the diﬀerent wave
functions Ψ, Ψ′ and the ground state energy E, E ′ according to the two diﬀerent Hamiltonians
H, H ′ associated with v, v ′ respectively. This leads the following two inequalities.
E ′ = 〈Ψ′ |H ′ |Ψ′ 〉< 〈Ψ|H ′ |Ψ〉= 〈Ψ|H − v + v ′ |Ψ〉∫(v ′ − v)ρdv

= E+

(12)

E = 〈Ψ|H|Ψ〉< 〈Ψ′ |H|Ψ′ 〉= 〈Ψ′ |H ′ − v ′ + v|Ψ′ 〉
∫

′

= E +

(v − v ′ )ρdv

(13)

The sum of these two inequalities gives the inconsistent result E + E ′ < E ′ + E. Therefore, v isa unique functional of ρ. In turn, the N -electron ground state wave function Ψ and all groundstate properties are functionals of ρ.
v = v[ρ](14)Ψ = Ψ[ρ](15)〈A〉 = 〈Ψ[ρ]|Â|Ψ[ρ]〉= A[ρ](16)

Next, the variational principle for Ev [ρ] is derived. When Ψ and Ψ′ are the ground statesfor given v and v ′ respectively, the following inequality is satisﬁed because of the variational
principle of Ev [Ψ] for Ψ.
Ev [Ψ] ≤ Ev [Ψ′ ](17)Now, we introduce a new energy functional Ev deﬁned as∫

Ev [ρ] =ρvdv + F [ρ]F [ρ] = T [ρ] + U [ρ]

(18)
(19)

The variational principle for Ev [ρ] for ρ is derived as
Ev [ρ] ≤ Ev [ρ′ ].

(20)

Because of the universality of the functional F [ρ], once we know it, the ground state propertiesof any system can be determined.

5.2

Kohn-Sham Equations

Kohn and Sham developed an exact way to treat an inhomogeneous electron system from thetheory of Hohenberg and Kohn. It has already shown that the total energy Ev of ground state
for a given external potential v is written as,
Ev = Ev [ρ]∫=vρdv +

12∫

ρ(r)ρ(r′ )drdr′ + G[ρ],|r − r′ |(21)

where G[ρ] was the universal functional of the density. G may be wrote as
G[ρ] = Ts [ρ] + EXC [ρ],
11

(22)

where Ts is the kinetic energy of a non-interacting system of electrons. All of the many bodyeﬀects are included into EXC . To make a discussion clear, an eﬀective potential Veﬀ is introduced,
and Ev is divided into two terms such as
Ev [ρ] = EB [ρ] + ED [ρ],(23)EB [ρ] = Ts +

(24)

∫

ED [ρ] = −Veﬀ ρdv,∫∫

Veﬀ ρdv +vρdv + EH + EXC .

Therefore, the variational equation given by Hohenberg and Kohn, i.e.
δEv
= 0,
δρ
is divided into a couple of the simultaneous equations:δEB
δED= 0, and= 0.
δρδρThen the following variational equation called Kohn-Sham equations are derived.
−∇2i ψi (r) + Vef f (r)ψi (r) = ϵi ψi (r)∫Vef f (r) = v(r) +ρ(r) =

∑

(25)

(26)

(27)

(28)
ρ(r′ )|r − r′ |dv ′ + VXC (r)(29)|ψi |2(30)

δ ẼXC [ρ]
δρ

(31)

i

VXC (r) =5.3

Local density approximation

For a homogeneous electron gas, the exchange-correlation energy is accurately calculated.
ẼXC [ρ]/N= εXC (ρ)(32)

In the local density approximation (LDA), we regard an inhomogeneous electron system as asystem locally homogeneous, i.e.,∫

ẼXC [ρ] =ρ(r)εXC (ρ(r)) dv.

(33)

Therefore, the exchange-correlation potential is given by
VXC (r) = εXC (ρ(r)) + ρ(r)5.4

dεXC (ρ(r)).
ρ(r)

Iteration

The self consistent Kohn-Sham equations are solved by an iteration procedure.
• Give a trial Veﬀ .
• Solve the single-electron Schrödinger equation.
• Calculate ρ(r) by summing∑2i |ψi | up to the Fermi level.

• Solve the Poisson equation for ρ(r) and construct the Hartree potential.
• Calculate Veﬀ using LDA.
• Compare the new and old Veﬀ and improve Veﬀ .
Repeat this procedure until the new and old potential become the same.

12

(34)

6

KKR-Green’s function method

6.1

Single-Site problem

Let us consider the scattering due to a single potential at the origin. This potential vanishes at
a region far from the origin.
6.1.1

Schrödinger equation

The Schrödinger equation which describes stationary states isHψ(r) = [−∇2 + V (r)]ψ = Eψ.

(35)

When the potential is spherical, the solutions are separable using the spherical coordinates:
ψ(r) = Rl (r)Ylm (θ, φ),
[

(36)

]

1 d
d
l(l + 1)
− 2 r2 ++ V (r) − E Rl = 0,r dr dr
r2

[

(37)

]

1
∂
∂
1
∂2
l(l + 1)sinθ+Ylm (θ, φ) = 0.
+r2 sin θ ∂θ∂θ r2 sin2 θ ∂φ2
r2

(38)

The solution of the angular part is spherical harmonics. Using Rl (r) = ul (r)/r, the radial
equation is written as[]d2
l(l + 1)
− 2+
+ V (r) − E ul (r) = 0.
dr
r2

(39)

The radial wave function is normalized as∫ RM T0

6.1.2[rRl (r; E)]2 dr = 1.

(40)

Asymptotic form

The exact radial wave function outside the potential is√√
Pl (r; E) = Al (E)jl ( Er) + Bl (E)nl ( Er).

(41)

jl and nl are the spherical Bessel and Neumann functions, which are regular and irregular atthe origin, respectively. Outside the potential, whole wave function satisfying any particular∑bounding conditions has a form ψ(r) = L CL (E)Pl (r; E)YL (θ, φ), where L = (l, m). As theboundary condition, the wave function which represents scattering should have the following
asymptotic form:√

ψ(r) → ei Ez +f (θ) i√Ere.
r

Therefore, Al (E)and Bl (E) can be expressed with a parameter ηl (E) as√√
Pl (r; E) = eiηl [cos ηl jl ( Er) + sin ηl nl ( Er)]√lπ
+ ηl ).
→ eiηl sin( Er −
2

13

(42)

(43)

ηl is the phase shift of the wave function doe to the scattering of the single potential. Using the
Hankel function of the ﬁrst kind,(1)hl= jl + inl ,√√(1)
Pl (r; E) = jl ( Er) + ieηl sin ηl hl ( Er)
√
√
√
(1)
= jl ( Er) − i Etl (E)hl ( Er).

(44)

(45)

As a results, the scattering problem becomes the problem to obtain the t-matrix deﬁned by1
tl (E) = − √ eiηl sin ηl
E
6.1.3

(46)

WronskianLet us consider to construct the solution of the Schrödinger equation outside the potential range
from the Wronskian relation. From the radial Schrödinger equations,[]l(l + 1)d2
+ V (r) − E rPl (r; E) = 0,
− 2+
dr
r2
] ( √
[
)
jl ( Er)l(l + 1)
d2
√
−E r
− 2+
= 0,
dr
r2
nl ( Er)the following integrals are evaluated.
∫ r√√rjl ( Er)V (r)rPl (r; E)dr = W [rjl ( Er), rPl (r, E)]r ,
0
∫ ∞
r

√
√
1
rnl ( Er)V (r)rPl (r; E)dr = − √ eiηl cos ηl − W [rnl ( Er), rPl (r, E)]r .
E

(47)
(48)

(49)
(50)

The ﬁrst term in Eq. (50) is obtained from the asymptotic form of Pl (r; E). Wronskian is deﬁned
asdb(r) da(r)
−
b(r).
(51)
dr
dr
√
√
√
Multipling Eq. (49) by nl ( E) and jl ( E) to Eq. (50) by nl ( E), and summing these equations,
we obtain√√ ∫ ∞ √√iηl
Pl (r; E) = e cos ηl jl ( Er) − E
jl ( Er< )nl ( Er> )V (r′ )Pl (r′ ; E)r′2 dr′
W [a(r), b(r)] = a(r)

0

∫ ∞ √
√
√
√
→ e cos ηl jl ( Er) − Enl ( Er)jl ( Er′ )V (r′ )Pl (r′ ; E)r′2 dr′0√√= eiηl [cos ηl jl ( Er) − sin ηl nl ( Er)]√√= Cl (E)jl ( Er) − Sl (E)nl ( Er)√√(1)= [Cl (E) − iSl (E)] jl ( Er) + iSl (E)hl ( Er),
iηl

(52)

where
Cl (E) = eiηl cos ηl√√=
EW [rPl (r; E), rnl ( Er)]r>RM T ,

(53)

iηl

Sl (E) = esin ηl√√= − EW [rPl (r; E), rjl ( Er)]r>RM T .

(54)

As a results, t-matrix is calculated from Wronskian:1
Sl (E)
tl (E) = − √.
E Cl (E) − iSl (E)
14

(55)

6.1.4

Green’s function in the free space

Green’s function in the free space is deﬁned as[∇2 + E]g(r, r ′) = δ(r − r ′).

(56)

The solution is obtained using the Fourier transformation
g(r, r ′) = g(r − r ′)√′e(i E|r −r |)= −.
4π|r − r ′|In addition, it is expanded into the partial waves as
√ ∑ √
√
(1)
g(r, r ′) = −i E
jl ( Er< )YL (r<)hl ( Er> )YL (r>).

(57)

(58)

L

Consider an inhomogeneous diﬀerential equation,[∇2 + E]φL (r) = V (r)φL (r).

(59)

Its general solution is expressed as a sum of its particular solution and the general solution to
the homogeneous diﬀerential equation,[∇2 + E]φ0L (r) = 0.

(60)

Since a particular solution can be obtained from Green’s function, the general solution is∫φL (r) = ZL (E)φ0L (r) +g(r − r ′)V (r′ )φL (r ′)dr ′,(61)

where ZL (E) is a normalization constant.

6.2
6.2.1

Traditional KKR
Muﬃn-tin potential

Usually, the KKR method is formulated using the muﬃn-tin potential approximation. Under
this approximation, the Schrödinger equation can be solved exactly.
The potentials are approximated as(v(r) =

v(r) for r < RM T
0 for r > RM T ,

)

(62)

where the origin is taken at the center of the atomic sphere. RM T is the radius of the spherewhich does not overlapped each other.
In the muﬃn-tin approximation, a constant potential is assumed in the interstitial region.
The origin of the energy is chosen such that potential at the interstitial region is zero as shownin Fig. 3.
6.2.2

Cell-centered expantion

Let us consider the Schrödinger equation,[−∇2 + V (r)]ψ(r) = Eψ(r).

15

(63)

rm
m
Rm
O

n
Rn

R

rn

MT

v(r)

r

Figure 3: Muﬃn-tin potential
Using the free space Green’s function, ψ(r) can be written as∫

ψ(r) =∫=

g(r, r ′)V (r ′)ψ(r ′)dv ′′)g(Rm + rm, Rn + rn

∑∫

=

∑

′ )ψ(R + r ′ )dv ′v(Rn + rnnnnn′ )v(r ′ )ψ(r ′ )dv ′ .
g(rm, Rn − Rm + rn
n
n
n

(64)

n

Here, ψ(r) is expanded into the atomic wave functions in the m-th sphere ψLm as
ψ(r) =∑

CLm φmL (rm).

(65)

L

φm
L satisﬁes the Shcrödinger equation in the m-th sphere and expressed by using g,m[−∇2 + v(r)]φm
L (r) = EφL (r),

∫

m
0
φm
L (r) = ZL (E)φL (r) +

(66)
′
′
g(r, r ′)v(r ′)φm
L (r )dv .

(67)

Substituting these expressions, we obtain∑∫

CLm [ZLm (E)φ0L (rm) +′ )v(r ′ )φm (r ′ )dv ′ ]g(rm, rm
m
m L
m

L

=

∑∫

′ )v(r ′ )g(rm, Rn − Rm + rn
n

n

6.2.3

∑

′ )dv ′ .
CLn φnL (rnnL

Structure constant′ ) is expanded into the partial waves asg(rm, Rn − Rm + rn
√ ∑ √′ ) = −i E
jl ( Erm )YL (rm)g(rm, Rn − Rm + rn
L

16

(68)

√
(1)
′ + R − R |)Y (r ′ + R − R )×hl ( E|rn
n
m
n
m
L n
√ ∑ √
jl ( Erm )YL (rm)= −i EL∑

×4π

L′ ,L′′√′′′(1)il−l +l CLL′ L′′ hl′′ ( E|Rn − Rm|)√′)×YL′′ (Rn − Rm)jl′ ( Ern′ )YL′ (rn
∑

mn 0
′
φ0L (rm)gLL′ φL′ (rn),(69)√ ∑ l−l′ +l′′√(1)
= −i4π E
i
CLL′ L′′ hl′′ ( E|Rn − Rm|)YL′′ (Rn − Rm).

(70)

=

LL′

mn
gLL
′

L′′

mn is called the structure constant. It depends only on the lattice structure, not on the potential
gLL′
at the lattice point.

6.2.4

KKR-matrix

From the partial wave expansion of Green’s function, we can obtain the equation that the
expansion coeﬃcients of the wave function should satisfy:∑∫

CLm [ZLm (E)φ0L (rm) +L=′ )v(r ′ )φm (r ′ )dv ′ ]g(rm, rm
m
m L
m

∑ ∫ ∑
n̸=m

∫

CLm ZLm (E)φ0L (rm) =

L

∑

=

′ )dv ′
CLn φnL (rnnLφ0L (rm)

∑ ∑

∫

mn
gLL
′

n̸=m L′L∑∑′ )v(r ′ )φm (r ′ )dv ′ ,g(rm, rm
m
m L
m

+
∑

LL′

mn 0
′
′
φ0L (rm)gLL′ φL′ (rn)v(rn)φ0L (rm)

∑ ∑

(71)

′ )v(r ′ )φn (r ′ )dv ′ C n
φ0L′ (rn
n L′
n L n

mn n
n
gLL
′ YL′ (E)CL′ .

(72)n̸=m L′

L

This leads to the secular equation
∑

mn n
n
[δLL′ δmn ZLm (E) − gLL′ YL′ (E)] CL′ = 0.

(73)

L′ n

Here, ZLm (E) and YLm′ (E) are
ZLm (E) = Clm (E) − iSlm (E),(74)

YLm′ (E)(75)=

Slm (E).

Therefore,mn n
det [δLL′ δmn − gLL′ tL′ (E)] = 0(76)

must be satisﬁed.

6.3
6.3.1

Crystal Green’s Function
Integral equation

In the previous section, we constructed the wave functions of scattered electrons using the freespace Green’s function. In this section, the Green’s function of the system G(r, r ′) is constructed.
[∇2 + E − V (r)]G(r, r ′) = δ(r − r ′)
17

(77)

Transposing the potential term to r.h.s. as[∇2 + E]G(r, r ′) = V (r)G(r, r ′) + δ(r − r ′)= f (r, r ′),we ﬁnd that G is expressed as
G(r, r ′)∫=g(r, r ′′)f (r ′′, r ′)dv ′′= g(r, r ′) +
6.3.2

(78)

∫

g(r, r ′′)V (r′′ )G(r ′′, r ′)dv ′′ .

(79)

Dyson equation

The free space Green’s function is expanded asg(r + Rm, r ′ + Rn) = δmn g(r, r ′) +∑mn′
JL (r)gLL′ JL′ (r ),

(80)

LL′

√ ∑
g(r, r ′) = −i E
JL (r<)HL (r>).

(81)

L

where

√
JL (r) = jl ( Er)YL (r),√(1)
HL (r) = hl ( Er)YL (r).

(82)
(83)

Similarly G(r, r′ ) is expressed as′
G(r + Rm, r ′ + Rn) = δmn Gms (r, r ) +∑′JL (r)Gmn
LL′ JL′ (r ),

(84)

LL′

√ ∑
′ ) = −i E
Gm(r,r
JL (r<)HL (r>),
s

(85)

L

√
JL (r) = Pl ( Er)YL (r),√HL (r) = Ql ( Er)YL (r).

(86)
(87)

The direct insertion of this formula to the integral equation yields the expansion coeﬃcients
Gmn
LL
mn
Gmn
LL′ = gLL′ +

∑

ml l
ln
gLL
′′ tL′′ GL′′ L′ .

(88)

L′′ l

6.3.3

Periodic boundary condition

mn
When the scattering centers are aligned periodically such as in a crystal, Gmn
LL′ and gLL′ are
Fourier transformed as

∫

Gmn
LL′

=

mn
gLL
′

=

1
GLL′ (k)e−ik(Rm−Rn) dk,
τ 1stBZ∫
1gLL′ (k)e−ik(Rm−Rn) dk.
τ 1stBZ(89)(90)

Substituting them into Eq. (88) yields
GLL′ (k) = gLL′ (k) +∑gLL′′ (k)tL′′ GL′′ L′ (k).

(91)

[δL′′ L − gL′′ L (k)tL ] gL′′ L′ (k).

(92)

L′′

We ﬁnally obtain the equation,
GLL′ (k) =

∑
L′′

18

6.3.4

Density of States

The electron density distribution is directly obtained from the Green’s function of the system.
It is shown from the eigen function expansion of the Green’s function.
The eigenfunction φn of the Schrödinger equation Hφ = Eφ which associated with the
eigenvalue En and the Green’s function of the system satisfy
Hφn = En φn ,(93)

[E − H]G(r, r ′) = δ(r − r ′).

(94)

Expanding G(r, r′ ) into φn ’s with the expansion coeﬃcients Gn (r′ ),∑G(r, r ′) =Gn (r ′)φn (r),

(95)

n

and multiplying it by (E − H) from the left, we obtain∑(E − H)G(r, r ′) =Gn (r ′)(E − H)φn (r),n= δ(r − r ′).

(96)

Multiplying Eq. (96) by φ∗m (r) and the volume integration leads toGn (r ′) = φ∗ (r ′)/(E − En ).

(97)

Therefore,∑ φ∗ (r ′ )φn (r)G(r, r ′) =nn

E − En(98)

is obtained. With an identity
1
1= P. − iπδ(x),x + iϵxG(r, r ′) =(99)

∑ φ∗ (r ′)φn (r)nn E + iϵ − En

= P.

∑ φ∗ (r ′ )φn (r)nn

E − En− iπ∑δ(E − En )φ∗n (r ′)φn (r)

(100)

n

is obtained. This gives the expression for the electron density distribution ρρ(r, E) =∑δ(E − En )φ∗n (r)φn (r)n1= − ImG(r, r).
π

(101)

Therefore, once the Green’s function of the crystal is known, we obtain the electron density
directly from it.

19

Contour Integration of the Green's Function
Im E

Im E

Ef

Ef

C

A
Re E

Re E

δ (E-En)Im G(z)2Γ
z=E+i ΓE

Figure 4: Contour integration of the Green’s Function
6.3.5

Contour integration

Now, the density distribution is calculated from the crystal Green’s function. To obtain thedensity distribution of electrons, it is needed to integrate the Green’s function with respect
to the energy accurately. However, it is rather diﬃcult because of the many-peaked structureof the density of states. To overcome this diﬃculty the contour integration technique is often
employed. First, the Green’s function is continuated to the complex energy plane.
G(r, r; z) =∑ φ∗ (r)φn (r)
n

z−E

n

=

∫ ∞ ∑

δ(E ′ − En )

−∞ n
∫
1 ∞

= −

π

−∞

φ∗n (r)φn (r) ′dEz − E′ImG(r, r ′; E ′ )/(z − E ′ )dE ′(102)

Here, G(z) is analytic in the whole complex plane except on the real axis. So we can deform theintegration path arbitrary in the complex plane.
n(r) =∑|φn (r)|2

n

∫ EF

=

−∞

ρ(r, E)dE∫

1 EF
ImG(r, r; E)dE
π −∞∫1
ImG(r, r; z)dz
= −
π A
∫
1
= −
ImG(r, r; z)dz
π C
= −

(103)

The contour integration can be performed accurately, because the peaked structure of the DOS
is smeared out. At z = E + iΓ, G(z) is
ImG(r, r; z) = ImG(r, r; E + iΓ)
20

1
π

=

∫ ∞
−∞

ImG(E ′ )ΓdE ′ .
(E ′ − E)2 − Γ2(104)This means the delta function on the real axis is broadened into the Lorentzian shape with thehalf-width of Γ (Fig. 4).

6.4

Coherent potential approximation (CPA)

Let us consider a random alloy system which consists of n components, A1 , A2 , ... An andtheir concentrations are x1 , x2 , ...xn . Suppose that the atom Ai is at the origin in the eﬀective
medium. Green’s function which starts from the origin and comes back to the origin is
GiLL′ =∑[

G̃LL′′ 1 − (ti − t̃)G̃

L′′

]−1
L′′ L′.

(105)

G̃LL′ and t̃ are the Green’s function and t-matrix (coherent t-matrix) of the eﬀective medium,
respectively. CPA is an eﬃcient approximation to determine t̃. We use the self-consistentequation,n∑xi GiLL′ = G̃LL′ .

(106)

i=1

This equation means that we determine the Green’s function of the medium by taking weighted
average of the Green’s function where the component atom is placed at the origin in the eﬀective
medium (see Fig. (5)).

xa

A

xb

B

Figure 5: Idea of CPA.

21

START

22

Figure 6: The ﬂow chart of the band structure calculation.

Construct starting potenpotenv.f
tial.

gsdatp.f atmicv.f

Atomic calculation.

No

Starting potential
exists?

gtchst.f

Structural constants.

Configure lattice vectors
and the reciprocal lattice
vectors.
genlpt.f madlng.fchklat.f prmvec.f
bzmesh.fWhich System ?

specx.f readin.f

Yes

Check neutrality andneutrl.f
shift EF.

gensm.f chrdnc.f rrcnv.f

Energy integral of theGreen function and get
charge and spin densities.

spckkr.f kkrsed.f
cinvrx.f bzmsmb.f

Back scattering.

Radial wave functions
and the atomic phase
phasea.f radial.f
shifts.

Define complex energy
drvmesh.f
contour.

Self-consistent
iteration starts.

Fix the core state conficorcnf.f
gulation.

No

END

dsenum.f hypera.f

Display results.

Iteration loop ends.

Yes

erranc.f

Check convergence.

banden.f totalw.f

Total energy calculation.

potenv.f

Solve the Poissn equation
and add the exchangecorrelation potential toconstruct new potential.

7
Subroutines

The KKR package consists of over 100 subroutines. Here, to get a survey of the KKR package,a brief ﬂowchart is given in Fig. 6. A comment on each subroutine is also given in table 4.

Table 4: Subroutines in the KKR package.
Subroutine
specx.freadin.f
spmain.fchklat.fprmvec.fbzmesh.fgenlpt.fmadlng.f
gtchst.fgsdatp.f
atmicv.f
corcnf.f
drvmsh.fphasea.fradial.ffczero.f
cstate.f
spckkr.fkkrsed.fcinvrx.f
bzmsmb.fchrdnc.fneutrl.fpotenv.fbanden.f
totalw.ferranc.fdsenum.f
hypera.fComment
Some parameters and arrays are deﬁned and spmain is called.
A given input ﬁle is read in.
Main routine of the KKR.
Check muﬃn tin radii and modify them if they conﬂict.
The primitive lattice vector and the reciprocal lattice vector are generated.
k-points in the Brillouin zone are set up.
Lattice points in the real and reciprocal space are generated.
The madelung constants are calculated.
The KKR structure constants are calculated.
Initial potentials are constructed from the atomic calculations.
The atomic calculation is performed.
The conﬁguration of the core states is conﬁrmed.
The energy contour is deﬁned.
The atomic t-matrix is calculated.
The radial Schrödinger equation is solved.
Single site Green’s function.
The radial Schrödinger equation is solved for the core states.
The crystal Green’s function is constructed.
The structure constant are interpolated for each energy point.
(1 − gt)−1 is calculated.
Integration in k-space.
The charge density distribution is calculated.
The charge neutrality is calculated and ϵf is shifted.
The Poisson equation is solved. The exchange-correlation potential is calculated.
The energy eigenvalues are summed up.
The total energy of the system is calculated.
The convergence is checked, and new input potential is constructed.
The results are printed out.
The hyperﬁne ﬁeld is calculated.

23

8

References• Scattering theory– J. J. Sakurai Modern Quantum Mechanics– L. I. Schiﬀ Quantum Mechanics– J. Alfken Mathematical Methods for Physicists• KKR method– W. Kohn, N. Rostoker Phys. Rev. 94(1954)1111.
– K. H. Johnson J. Chem. Phys. 45(1966)3085.
– F. Herman, A. R. Williams, K. H. Johnson J. Chem. Phys. 61(1974)3508.
– R. Podloucky, R. Zeller, P. H. Dederichs Phys. Rev. 22(1980)5777.
– R. Zeller, J. Deutz, P. H. Dederichs Solid State Comm. 44(1982)993.
– H. Akai’s note on Korringa-Kohn-Rostoker Method (1993).
– T. Kotani, H. Akai Phys. Rev. B54(1996)16502.
• Density functional theory– P. Hohenberg, W. Kohn Phys. Rev. 136(1964)B864.
– W. Kohn, L. J. Sham Phys. Rev. 140(1965)A1133.
• Package (Machikaneyama2002, AkaiKKR)– http://sham.phys.sci.osaka-u.ac.jp/˜kkr/

24