- Email: [email protected]
Parametrization of the band structure of FCC crystals
Computer Physics Communications 30 (1983) 271—275 North-Holland Publishing Company
271
PARAMETRIZATION OF THE BAND STRUCTURE OF FCC CRYSTALS Carlo SALUSTRI
*
Max - Planck - Inst itut für Metallforschung, Inst itut für Physik, Heisenbergstrasse 1, 7000 Stuttgart 80, Fed. Rep. Germany
Received 1 July 1982; in revised form 4 May 1983
PROGRAM SUMMARY Title of the program: BAPAR
No. of lines in combined program and test deck: 445
Catalogue number: ACKU
Keyword: band structure interpolation
Program obtainable from: CPC Program Library, Queen’s Uni-
Nature of the physical problem
versity of Belfast N. Ireland (see application form in this issue)
The program fits the energy bands of fcc crystals.
Computer: CYBER 170-720; Installation: Ecole Polytechnique
Method of solution
Federal, Lausanne
The method consists in treating the integrals contained in the Hamiltonian matrix as independent parameters. With a starting Set of approximated values of these integrals, the program evaluates the eigenvalues of the Hamiltonian matrix and calcu-
Operating system: NOS.BE Programming language used: FORTRAN
lates the mean value of the difference between them and the energy values of the bands to be fitted. The user has to give a
High -speed storage required: 14000 words No. of bits in a word: 60
limit to this difference; the program modifies the parameters until this limit is reached, finding at the ened of the iterations the “best estimate” set of parameters which fits the true band.
Peripherals used: card reader, line printer
60—70 s.
Typical running time
*
Present address: Istituto di Elettronica, dello Stato Solido (CNR), Via Cineto Romano 42, 00156 Roma, Italy.
001 0-4655/83/0000—0000/$03.00
©
1983 North-Holland
272
C. Salustri
/
Band structure of fcc crystals
LONG WRITE-UP 1. Introduction
thogonal between the same and between different sites [2]:
Calculations of the energy bands in solids have been carried out extensively in the last few years, using different methods, e.g., the OPW or APW approach, to solve the one-electron Schrodinger equation for the crystal. The main difficulty in handling the results of such calculations is the restricted number of k points at which the energy values are available. The present paper follows the scheme proposed by Slater and Koster to interpolate the results of these calculations by considering the energy matrix elements between localized functions as fitting parameters to be adapted to the true energy band [1].The number of these parameters is determined by the number of bands and the range of the interatomic interactions one takes into account. The method yields analytical expressions of the energy as a function of the wave vector k, starting from the energy values calculated at restricted symmetry points of the Brillouin zone (normally 89 for fcc crystals).
)
fw~(r — R, wM( r — R, )dr
=
~
(4)
Putting eq. (2) into eq. (1), we see that the E,(k) and the coefficients b,~(k)satisfy the follow-
ing secular equation: ~A~(k)b(k)
=
E~(k)h~(k),
(5)
where A~(k) =
~
exp[ik.(R,
—
R
1)]
R1
x fw~(r—R,)H0w~(r—R,)dr.
(6)
We see in eq. (6) that the integrals
I~(R1, R,) =fw~(r_ R~)H0w~(r— R1)dr
(7)
are in practice independent parameters. We want to determine a “best estimated” set of 2. Theory The eigen-functions of the one-electron Hamiltonian H0 are Bloch functions p~,(k,r) which satisfy the SchrOdinger equation:
parameters that approximate the I~(R,,R1). In other words, we have to find those parameters that, eq. (5), provide eigenvalues ç,(k, P) such put that into the difference 2(P) ~ [r~(k, P) E,,(k)}2 (8) =
—
E~(k)~1(k,r), (I) where the solutions E~(k) represent the energy
~
bands of the crystal. They will be in general degenerate at certain k points, We describe the Bloch functions as linear combinations of orthogonal wave functions:
is minimum [3]. In eq. (8), P is a vector whose components are the parameters to find and M is the number of ab initio calculated eigenvalues. In practice, instead of minimizing eq. (8), one imposes the condition
H0(r)~(k, r)
~
r)
=
=
~b,~(k)u~(k,
n
.
(2)
r),
(9)
I)
where u,,(k, r)
where 6 is a small number. If a first set of approximated parameters does not satisfy condition =
~
exp(ik-R 1)w,,(r
—
R1).
(3)
(9), one must modify them with the condition
R1
In eq. (3) the w~(r—R1)have the property of being localized on the lattice site R, and are or-
(10) Assuming that ,,(k, P) is linear in P, and
C. Salustri
/
Band structure offcc crystals
defining: D,1
=
which the difference ~ P) E~(k) remains greater than a value y. Unfortunately, this number —
~ d~~(k, F) n,k dP1
=
~
[~(~)
n,k
(11)
th,,(k, P) p(0)
dP1
p(0)
E~(k,p(0))] d~~(k, P) d~
p0’
(12) where dc,1(k, dP, F)
273
1 ‘c’ b’ dA —~
=
(13)
b,~,
it can be shown that a new set of parameters, obtained from the first one by:
oscillates, but, to minimize it, the user can let the computer continue the iterations after the 6 value is reached. In the original version, the program stops when the 6 value is reached, whatever the value of POUT is. In this version two EISPACK subroutines (TRED2 and IMTQL2) for theDdiagonalization of for theA inversion of matrix are used. IMGC These matrix and one EPFL library subroutine should if necessary be replaced by locally-available equivalents, and their actions are described in the listing.
p 1
=
(14)
p,~ + ~(D~)11l1
3.3. Input data
I
satisfy condition (10).
The input data are approximated starting values of the 32 parameters and the values of 6 and y. The values resulting from our test run for ferro-
3. Program structure
magnetic nickel can be used as starting values for many fcc crystals.
3.1. Application to ferromagnetic nickel We applied the procedure described above in fitting the first six bands of ferromagnetic nickel (spin up) obtained with self-consistent calculations by Skriver [4]. The electronic structure of ferromagnetic nickel can be adequately described by a 9 X 9 Hamiltonian matrix, corresponding to one s-, five d- and three p-functions with a total of 32 independent parameters, taking into account first and second
‘~
Test run
In our test, 6 0.005 Ry and y 0.008 Ry. At the last iteration POUT was 104. A better parametrization could be reached with double precision calculations and more computer time. The values of the 32 parameters we obtained =
=
0.2-
7
0.0-
neighbours. 3.2. Program flow
-0.2-
The program solves eq. (5) with an initial set of parameters, calculates eigenvalues and eigenvectors of the matrix A0~(k,P) and modifies the parameters according to eq. (14) to reach convergence of the mean difference ~(P). During each computer iteration the difference )((P) is calculated at 89 points of the Brillouin zone (M 534). The parameters are then modified until condition (9) is satisfied. During each iteration the program calculates also the number (POUT) of k points at =
~ ~ ~ Lu
-
0.6__________
0.8 I
Fig. I.
t~
________
x
—r——
Z
W
Q
L
A
I
K
Interpolated bands of ferromagnetic nickel in the case of
spin up. The continuous line represents the results of the test
run, while the dots are the energy values obtained by Skriver [4]
274
C. Salustri
/
Band structure of fcc crystals
Table I Values of the 32 best interpolated parameters in the 9 > 9 Hamiltonian matrix in the case of ferromagnetic nickel (spin up). The indices are as in eq. (7) with R = 0 P,,
(000)=0.35l491 (llO)= —0.114658
~ P~
P~,
(200) = 0.038 169 (ll0)= —0.097740 (220) 0.007 102 (llO)=0.066683 (110) = —0.038515 (002) = 0.029531 (000)=0.613431 (llo)=0.125575 (011) = 0.034001 (200)=0.082787 (200)=).024894 (011) = ).064696 (Oll)=0.071194 (002) = 0.001775
~
~ P~3~,z
~
P. P~3~2~2
are listed in table 1, while fig. 1 shows the very good quality of the fit.
References [lJJ.C. Slater and G.F. Koster, Phys. Rev. 94 (1954) 1498. [2]P.O. Lowdin, J. Mol. Spectry. 13 (1964) 326.
(000)=0.016726 0~
P3_3.~_~ P3~2_,.23,2_r2
P3~3~~
P~
(ll0)=0.022061 (011) = 0.023220 (200)= 0.034963 (002) = —0.197050 (Oll)= —0.026851 (110) = 0.008032 (000) = —0.014354 (llO)=0.00l 189 (llO)= —0.002756 (002) = — 0.007864 (002)= —0.274037 (llO)= —0.021915 (110) = 0.033434 (0Il)= —0.003759 (020) = 0.004464
[3] J.W.D. Connolly, in: Electronic Density of States. Nat!. Bur. Stand. Special PubI. No. 323 (US GPO, Washington, DC, 1971) p. 27. [4]H. Skriver, private communication. See also H. Skriver and 0. Andersen, Inst. Phys. Conf. No. 39. cap. 1(1978)100.
C. Salustri
/
Band structure offcc ctystals
TEST RUN OUTPUT • )341354795 PJUT=
51t~
SLIMNK PtJUT=
1~4
•
PC
THE ~r~J !rs1’E,~PJLATE[i PARAIETE~S ARE: 1)= • 35L4~1 2l~ —.114058 31= .L3~ 4I~ —. J~174( 51~ 81= .Cu~d3 71= —.03a51.5 -31= .2’4531
P(
fl=
PC PC j)(
PC PC Pt
PC
P(1Ll=
.813431 • 125575 •OJ41l
P(1,2)=
•
P(tQl=
P(13)= PC~.5l= Pt18)=
.3711.9’,
PC17)=
•‘~1b726
P(i3)= P(191= Pt~~ ))
.~)b1 •(-~3,~2(~ —.
19735~~
P(22)= PC23) 354 .O~-L1.a9 P1281= Pt~’7)=
—.O;7364
—.27’,)37 —.C2L~15 • :33~,34
P(31)
—~C.33759
.C-(os4~4
$S8
275
271
PARAMETRIZATION OF THE BAND STRUCTURE OF FCC CRYSTALS Carlo SALUSTRI
*
Max - Planck - Inst itut für Metallforschung, Inst itut für Physik, Heisenbergstrasse 1, 7000 Stuttgart 80, Fed. Rep. Germany
Received 1 July 1982; in revised form 4 May 1983
PROGRAM SUMMARY Title of the program: BAPAR
No. of lines in combined program and test deck: 445
Catalogue number: ACKU
Keyword: band structure interpolation
Program obtainable from: CPC Program Library, Queen’s Uni-
Nature of the physical problem
versity of Belfast N. Ireland (see application form in this issue)
The program fits the energy bands of fcc crystals.
Computer: CYBER 170-720; Installation: Ecole Polytechnique
Method of solution
Federal, Lausanne
The method consists in treating the integrals contained in the Hamiltonian matrix as independent parameters. With a starting Set of approximated values of these integrals, the program evaluates the eigenvalues of the Hamiltonian matrix and calcu-
Operating system: NOS.BE Programming language used: FORTRAN
lates the mean value of the difference between them and the energy values of the bands to be fitted. The user has to give a
High -speed storage required: 14000 words No. of bits in a word: 60
limit to this difference; the program modifies the parameters until this limit is reached, finding at the ened of the iterations the “best estimate” set of parameters which fits the true band.
Peripherals used: card reader, line printer
60—70 s.
Typical running time
*
Present address: Istituto di Elettronica, dello Stato Solido (CNR), Via Cineto Romano 42, 00156 Roma, Italy.
001 0-4655/83/0000—0000/$03.00
©
1983 North-Holland
272
C. Salustri
/
Band structure of fcc crystals
LONG WRITE-UP 1. Introduction
thogonal between the same and between different sites [2]:
Calculations of the energy bands in solids have been carried out extensively in the last few years, using different methods, e.g., the OPW or APW approach, to solve the one-electron Schrodinger equation for the crystal. The main difficulty in handling the results of such calculations is the restricted number of k points at which the energy values are available. The present paper follows the scheme proposed by Slater and Koster to interpolate the results of these calculations by considering the energy matrix elements between localized functions as fitting parameters to be adapted to the true energy band [1].The number of these parameters is determined by the number of bands and the range of the interatomic interactions one takes into account. The method yields analytical expressions of the energy as a function of the wave vector k, starting from the energy values calculated at restricted symmetry points of the Brillouin zone (normally 89 for fcc crystals).
)
fw~(r — R, wM( r — R, )dr
=
~
(4)
Putting eq. (2) into eq. (1), we see that the E,(k) and the coefficients b,~(k)satisfy the follow-
ing secular equation: ~A~(k)b(k)
=
E~(k)h~(k),
(5)
where A~(k) =
~
exp[ik.(R,
—
R
1)]
R1
x fw~(r—R,)H0w~(r—R,)dr.
(6)
We see in eq. (6) that the integrals
I~(R1, R,) =fw~(r_ R~)H0w~(r— R1)dr
(7)
are in practice independent parameters. We want to determine a “best estimated” set of 2. Theory The eigen-functions of the one-electron Hamiltonian H0 are Bloch functions p~,(k,r) which satisfy the SchrOdinger equation:
parameters that approximate the I~(R,,R1). In other words, we have to find those parameters that, eq. (5), provide eigenvalues ç,(k, P) such put that into the difference 2(P) ~ [r~(k, P) E,,(k)}2 (8) =
—
E~(k)~1(k,r), (I) where the solutions E~(k) represent the energy
~
bands of the crystal. They will be in general degenerate at certain k points, We describe the Bloch functions as linear combinations of orthogonal wave functions:
is minimum [3]. In eq. (8), P is a vector whose components are the parameters to find and M is the number of ab initio calculated eigenvalues. In practice, instead of minimizing eq. (8), one imposes the condition
H0(r)~(k, r)
~
r)
=
=
~b,~(k)u~(k,
n
.
(2)
r),
(9)
I)
where u,,(k, r)
where 6 is a small number. If a first set of approximated parameters does not satisfy condition =
~
exp(ik-R 1)w,,(r
—
R1).
(3)
(9), one must modify them with the condition
R1
In eq. (3) the w~(r—R1)have the property of being localized on the lattice site R, and are or-
(10) Assuming that ,,(k, P) is linear in P, and
C. Salustri
/
Band structure offcc crystals
defining: D,1
=
which the difference ~ P) E~(k) remains greater than a value y. Unfortunately, this number —
~ d~~(k, F) n,k dP1
=
~
[~(~)
n,k
(11)
th,,(k, P) p(0)
dP1
p(0)
E~(k,p(0))] d~~(k, P) d~
p0’
(12) where dc,1(k, dP, F)
273
1 ‘c’ b’ dA —~
=
(13)
b,~,
it can be shown that a new set of parameters, obtained from the first one by:
oscillates, but, to minimize it, the user can let the computer continue the iterations after the 6 value is reached. In the original version, the program stops when the 6 value is reached, whatever the value of POUT is. In this version two EISPACK subroutines (TRED2 and IMTQL2) for theDdiagonalization of for theA inversion of matrix are used. IMGC These matrix and one EPFL library subroutine should if necessary be replaced by locally-available equivalents, and their actions are described in the listing.
p 1
=
(14)
p,~ + ~(D~)11l1
3.3. Input data
I
satisfy condition (10).
The input data are approximated starting values of the 32 parameters and the values of 6 and y. The values resulting from our test run for ferro-
3. Program structure
magnetic nickel can be used as starting values for many fcc crystals.
3.1. Application to ferromagnetic nickel We applied the procedure described above in fitting the first six bands of ferromagnetic nickel (spin up) obtained with self-consistent calculations by Skriver [4]. The electronic structure of ferromagnetic nickel can be adequately described by a 9 X 9 Hamiltonian matrix, corresponding to one s-, five d- and three p-functions with a total of 32 independent parameters, taking into account first and second
‘~
Test run
In our test, 6 0.005 Ry and y 0.008 Ry. At the last iteration POUT was 104. A better parametrization could be reached with double precision calculations and more computer time. The values of the 32 parameters we obtained =
=
0.2-
7
0.0-
neighbours. 3.2. Program flow
-0.2-
The program solves eq. (5) with an initial set of parameters, calculates eigenvalues and eigenvectors of the matrix A0~(k,P) and modifies the parameters according to eq. (14) to reach convergence of the mean difference ~(P). During each computer iteration the difference )((P) is calculated at 89 points of the Brillouin zone (M 534). The parameters are then modified until condition (9) is satisfied. During each iteration the program calculates also the number (POUT) of k points at =
~ ~ ~ Lu
-
0.6__________
0.8 I
Fig. I.
t~
________
x
—r——
Z
W
Q
L
A
I
K
Interpolated bands of ferromagnetic nickel in the case of
spin up. The continuous line represents the results of the test
run, while the dots are the energy values obtained by Skriver [4]
274
C. Salustri
/
Band structure of fcc crystals
Table I Values of the 32 best interpolated parameters in the 9 > 9 Hamiltonian matrix in the case of ferromagnetic nickel (spin up). The indices are as in eq. (7) with R = 0 P,,
(000)=0.35l491 (llO)= —0.114658
~ P~
P~,
(200) = 0.038 169 (ll0)= —0.097740 (220) 0.007 102 (llO)=0.066683 (110) = —0.038515 (002) = 0.029531 (000)=0.613431 (llo)=0.125575 (011) = 0.034001 (200)=0.082787 (200)=).024894 (011) = ).064696 (Oll)=0.071194 (002) = 0.001775
~
~ P~3~,z
~
P. P~3~2~2
are listed in table 1, while fig. 1 shows the very good quality of the fit.
References [lJJ.C. Slater and G.F. Koster, Phys. Rev. 94 (1954) 1498. [2]P.O. Lowdin, J. Mol. Spectry. 13 (1964) 326.
(000)=0.016726 0~
P3_3.~_~ P3~2_,.23,2_r2
P3~3~~
P~
(ll0)=0.022061 (011) = 0.023220 (200)= 0.034963 (002) = —0.197050 (Oll)= —0.026851 (110) = 0.008032 (000) = —0.014354 (llO)=0.00l 189 (llO)= —0.002756 (002) = — 0.007864 (002)= —0.274037 (llO)= —0.021915 (110) = 0.033434 (0Il)= —0.003759 (020) = 0.004464
[3] J.W.D. Connolly, in: Electronic Density of States. Nat!. Bur. Stand. Special PubI. No. 323 (US GPO, Washington, DC, 1971) p. 27. [4]H. Skriver, private communication. See also H. Skriver and 0. Andersen, Inst. Phys. Conf. No. 39. cap. 1(1978)100.
C. Salustri
/
Band structure offcc ctystals
TEST RUN OUTPUT • )341354795 PJUT=
51t~
SLIMNK PtJUT=
1~4
•
PC
THE ~r~J !rs1’E,~PJLATE[i PARAIETE~S ARE: 1)= • 35L4~1 2l~ —.114058 31= .L3~ 4I~ —. J~174( 51~ 81= .Cu~d3 71= —.03a51.5 -31= .2’4531
P(
fl=
PC PC j)(
PC PC Pt
PC
P(1Ll=
.813431 • 125575 •OJ41l
P(1,2)=
•
P(tQl=
P(13)= PC~.5l= Pt18)=
.3711.9’,
PC17)=
•‘~1b726
P(i3)= P(191= Pt~~ ))
.~)b1 •(-~3,~2(~ —.
19735~~
P(22)= PC23) 354 .O~-L1.a9 P1281= Pt~’7)=
—.O;7364
—.27’,)37 —.C2L~15 • :33~,34
P(31)
—~C.33759
.C-(os4~4
$S8
275