The very old and successful density-functional technique of half-occupation is revisited [J. C. Slater, Adv. Quant. Chem. 6, 1 (1972)]. We use it together with the modern exchange-correlation approximations to calculate atomic ionization energies and band gaps in semiconductors [L. G. Ferreira et al., Phys. Rev. B 78, 125116 (2008)]. Here we enlarge the results of the previous paper, add to its understandability, and show when the technique might fail. Even in this latter circumstance, the calculated band gaps are far better than those of simple LDA or GGA. As before, the difference between the Kohn-Sham ground state one-particle eigenvalues and the half-occupation eigenvalues is simply interpreted as the self-energy (not self-interaction) of the particle excitation. In both cases, that of atomic ionization energies and semiconductor band gaps, the technique is proven to be very worthy, because not only the results can be very precise but the calculations are fast and very simple.

The very old half-occupation technique1–3 and its formalization through the Janak theorem4 have been continuously used since the seventies. Even now we find many papers citing those authors and using their techniques in many computational problems.5–13 The Janak theorem refers to the derivative of the total energy with respect to the occupation, which equals the Kohn-Sham (KS) eigenvalue.14 If the derivative is linear with the occupation then the half-occupation technique follows. Part of the recurring publications on the subject refers to verifying the linearity,15 and this since the early days of the technique.16 The first applications of the half-occupation technique were for atoms, thus we revisit this problem, now using the modern exchange-correlation approximations.17,18 In this respect, our work runs parallel to that of Kraisler et al.,11 among others, though they did not use the half-occupation technique, which is a matter of necessity in our case, as we will see in the case of extended systems (crystals).

After atoms, molecules and clusters, the half-occupation technique must face the problem of extended systems such as crystals. In this respect, aside from our work, being reviewed, we only know the work of Liberman.19 Of course, removing 1/2 electron from an infinite system is no perturbation, so that a special handling of the half-occupation technique is required. Related to this special technique is the definition of “self-energy”, or the “self-interaction”, known to be missing in the Kohn-Sham eigenvalues.20 The self-interaction correction (SIC) also has a long history. It was first defined by Perdew and Zunger17 and frequently proposed as the difference between the ionization energy (IE) and the KS eigenvalue.20,21 It happens that whenever the Janak theorem is valid, namely when the energy functional can be differentiated with respect to the occupations, one can prove that the difference between the IE and the KS eigenvalue, is about 1/2 the SIC defined by Ref. 17. This factor of 1/2, apparently first discovered by Trickey,22 has also been numerically verified by Refs. 23. To avoid confusion we are naming the difference between the IE and the KS eigenvalue, following from Janak's theorem, as “self-energy” because its expression is indeed that of a self-energy.

The technique was proposed by Slater's group in the early seventies.1–3 At that time the density-functional (DF) approximations were unknown and researchers were using the X − α exchange to obtain the best results available at the time. The half-ion technique gained its formalization much latter with Janak4 which proved

(1)

In this Eq. fα is the occupation of the Kohn-Sham (KS) eigenstate,14 not necessarily an integer, eα is the KS eigenvalue, dependent on the occupation fα, and E(fα) is the total energy of the atom (or ion). Since the seminal works of Slater and Janak, the half-ion technique is being intensively used to this day. The proof of Eq. (1) is certainly exact for the local-density approximation (LDA) to DF, and it is likely to be also valid for the generalized-gradient approximation (GGA), though GGA may come in very different formulations.

The half-ion technique begins by integrating Eq. (1) from the ion to the atom

(2)

The last equality assumes that the KS eigenvalue is linear with the occupation. The linearity was verified very early in the development of the half-ion technique15,16 and is now being confirmed by a very large set of atoms. Of course the half-ion technique cannot be as precise as the difference of total energies, yet its precision is typically better than 0.1eV and, sometimes it is an important tool of research.24,25 In the case of extended systems, as crystals, it is unavoidable. Instead of the approximation of Eq. (2), ref. 13 among other authors, uses

which cannot be interpreted as the energy of a transition state and, as Numerical Analysis says, it is a slightly poorer approximation.

Apparently the authors of the half-ion technique failed to give an explanation of why it works and the physical meaning for the difference between the half-ion eigenvalue eα( − 1/2) and the eigenvalue eα(0) at full occupation. An understanding may be reached starting from

(3)

where

$n_{\alpha }=\psi ^*_{\alpha }\psi _{\alpha }$
nα=ψα*ψα is the number density of the KS state α, n = ∑ifini is the number density of the whole electron cloud, and KE[] is the one-electron kinetic energy functional. Taking the derivative of Eq. (3) with respect to fα obtain

(4)
(5)

The first term of Sα is typically the self-energy of the charge distribution nα. In many instances this first term will be the most important and that is why we are naming Sα as self-energy. A numerical comparison between the first term and the whole self-energy is in the footnote.26 The first two terms of Eq. (5) coincide with Trickey's post-hoc correction,22 except for the unit systems being used. The self-energy is always positive and it is larger if the number-density nα is denser (core states) and smaller if it is diffuse (valence states). In deriving Eq. (3) to obtain Eq. (5) we apparently miss the derivative of nα with respect to fα, except when it is contained in the total density n. The reason is that the eigenvalue is stationary with respect to variations in the normalized wavefunctions ψα, thus the derivative is zero (Hellman-Feynman theorem).27 It must be stressed that the self-energy will never be calculated with Eq. (5) but in a much simpler way, even simpler than by using any known self-interaction expression.17,22,28 The simplicity makes the method very useful.

Integrating Eq. (4) obtain

(6)

the last equality again depending on linearity, which again may be assumed because it is usually verified. Truly, the self-energy tends to be independent of the occupations.

Now inserting Eq. (6) into Eq. (2) obtain

(7)

Meaning that the ionization of electron α requires not only the negative of the eigenvalue energy but the self-energy as well.

We calculated many ionization energies (IE) with the half-occupation technique, using both the LDA approximation of Perdew and Zunger17 to the results of Ceperley and Alder29 and the GGA of Perdew, Burke, and Ernzerhof.18 Our calculations were all with spin polarization,30 thus for the heavier atoms we had to make a perturbation correction to include relativistic effects.

Our results are compared with experiment and with calculated results by other methods. The NIST site31 presents very large tables of IE data. Of course we cannot compare our results with coupled-clusters. Coupled-cluster is known to reach a precision of sub-meV but only for very light atoms,32 for heavier atoms their precision is not better than ours.31 We are mainly interested in extended systems, specially in semiconductors, where heavier atoms are important. Thus we are comparing our results calculated with spin-polarization, radial grid integration, LDA and GGA half-ion technique with basis functions results, Hartree-Fock (HF) and Density-Functional (DF) results presented in the NIST site. In the NIST site we chose the results for the largest basis (aug-cc-pVQZ, when available), and for the honorable B3LYP functional. The comparison is in Fig. 1. As expected, only HF has important deviations from experimental values. The results presented in Ref. 11 also seem to coincide with ours and NIST’s, barring HF, though a detailed comparison is not possible because those authors do not tabulate their results. Ref. 11 also obtains the ionization energies as the difference between the total energies of atom and ion, not from half-occupation. Looking more closely, at the rms of the difference between calculated and experimental values, our GGA-1/218 calculation is outstanding. Its superiority with respect to LDA-1/217 only proves that, in searching for better exchange- correlation approximations, scientists have not been idle. The apparent superiority of GGA-1/2 with respect to B3LYP might be a problem of method (numerical against bases), not necessarily a question of quality of the approximation.

FIG. 1.

(color online) Ionization energies: Comparison of calculated results with experiment for the atoms tabulated in the figure. The experimental data, HF, and B3LYP are taken from the NIST site.31 The LDA-1/2 uses the exchange-correlation approximation of ref. 17. The GGA-1/2 uses ref. 18.

FIG. 1.

(color online) Ionization energies: Comparison of calculated results with experiment for the atoms tabulated in the figure. The experimental data, HF, and B3LYP are taken from the NIST site.31 The LDA-1/2 uses the exchange-correlation approximation of ref. 17. The GGA-1/2 uses ref. 18.

Close modal

One important advantage of the half-ion technique is its ability to calculate heavier atoms and higher order IE. To show this ability we plot in Fig. 2 the first and second IE's for the atoms tabulated in reference 31 and some heavier atoms. In the calculation of the second IE one must delete 3/2 electrons from the atom. The inclusion of the heavier atoms increased the rms error to some extent. The second IE has higher rms error than the first. In the case of Cs and Tl we had to apply a perturbation scheme to include the relativistic corrections into the spin-polarized calculation.

FIG. 2.

(color online) Ionization energies: Comparison of calculated results for the first and second IE with experiment for the atoms tabulated in the Fig. The exchange-correlation used was GGA of 18.

FIG. 2.

(color online) Ionization energies: Comparison of calculated results for the first and second IE with experiment for the atoms tabulated in the Fig. The exchange-correlation used was GGA of 18.

Close modal

Following Eq. (6), the IE, which equals minus the eigenvalue at half occupation, is obtained if we subtract the self-energy Sα from the eigenvalue at full occupation. The self-energy is the quantum-mechanical average

(8)

of the “self-energy potential” VS defined as follows

(9)

The self-energy potential keeps a simple relationship with the KS potential for the half ion and the atom. Indeed consider

(10)

Assuming linearity and comparing with Eq, (9) obtain

(11)

Therefore, instead obtaining IE from the eigenvalue of the half ion, we may as well subtract the self-energy potential from the potential of the atom and run a single non self-consistent iteration to determine the eigenvalue at half occupation. After the subtraction the potential is V( − 1/2, r), no longer V(0, r), and the eigenvalue of a single iteration is eα( − 1/2). This seems an undue complication in the case of atoms, but it is important in the case of many-atom systems. In those cases we will want to use the self-energy potential defined in the atom and subtract it from the potentials of the atoms (of the same species) composing the system.

First we recall that the KS eigenstates in a crystal are Bloch states whose wavefunctions extend through the infinite crystal. As such, the Bloch state has an exactly null self-energy. Now, when removing or exciting an electron from an occupied band, say the topmost valence band, the excitation or hole does behave as being in a localized wavefunction, therefore with non null self-energy. In fact, it would be impossible to remove an electron infinitely delocalized in any finite time. Further, theories of excitons, potential barriers, etc, would be in trouble if holes were Bloch states. Though these arguments seem to be very simple, there is no calculation of a localized hole in a solid that we know of, specially in the case when such localized hole could scatter into the continuum of band states, or, in other words, when the localized hole is degenerate with the band continuum.33 

Thus, to proceed with the half-occupation for crystals, we have to assume that holes in a valence band, and electrons in a conduction band are localized with non-null self-energy. Generalizing Eqs. (2) and (6), to describe the excitation of an electron from the valence to the conduction band, we write

(12)

This last equation is frequently writen as

where Δxc is known as the “derivative discontinuity” that is caused by an exchange-correlation beyond LDA or GGA.34–36 The half-occupation technique allows an entirely different, and perhaps simpler, interpretation of this discontinuity: It is a self-energy caused by the localization of the particle excitation and it is not the effect of exotic exchange-correlation functionals.

Since we do not know how to calculate localized holes and electrons in solids, which are not even stationary states, we assume that their self-energies are again quantum mechanical averages of the atomic self-energy potential VS. In a crystal such as NaCl we would add VS(Cl) to the chlorine atoms and VS(Na) to the sodiums. We add to all atoms of the crystal because the calculation with just one or a few perturbed atoms is much too complicate.

There is a problem with adding VS to all the atoms of an infinite crystal: the potential will diverge. A VS is a potential of an excess charge of 1/2 proton (see Eq. (9)) and has a tail of 0.5/r that cannot be summed in an infinite lattice. Therefore the tail has to be trimmed by a step function, and we have been using

(13)

Thus we are making

(14)

In this Eq., the exponent m is as large as possible and it is limited only by the band-calculation code. We have been using always m = 8. Thus we are introducing a single parameter in the calculation, CUT, which is determined variationally, as explained in ref. 37.

Thus our recipe, named LDA-1/2, to determine band gaps, consists of adding the atomic self-energy potential to the atoms of the crystal. Now, in all our calculations with crystals, and we made many, we never found an instance when the self-energy of the electron in the conduction band had any importance in Eq. (12). The valence band is usually made of the valence state of the anion, while the conduction band is usually a mixture of many atomic states, there included the valence state of the cation. Then, usually the self-energy of the cation could be neglected, or, which is the same, the extremizing solution usually was CUTcation = 0. In this respect we found the calculation of the work functions of Al, Cu, Na very symptomatic. A work function is a band gap between the Fermi level at the conduction band and the free electron band at infinity. In those cases we always found CUT = 0, that is the LDA-1/2 method reducing to standard LDA.

Fig. 3 shows the situation of the exceptional case ZnO. This material is exceptional in two ways: first the pure LDA or GGA calculations give a band gap smaller than 1.0 eV, while our LDA-1/2 is able to bring the gap to a value close to the experimental 3.4 eV. A major shift and a major success of our technique. Secondly, it is exceptional because the self-energy of the cation Zn cannot be neglected. Fig. 3 explains how the calculation was made. We began by including the self-energy potential for the anion

(15)

and by maximizing the gap with respect to CUTO. Then with the maximizing CUTO fixed, we included the self-energy potential for Zn

(16)

A new maximization was made, resulting in a gap near the experimental value.

FIG. 3.

(color online) The figure depicts the choice of CUT for O and Zn in the case of ZnO. First we maximize the gap varying CUTO, then we vary CUTZn to reach a band gap value near the correct value.

FIG. 3.

(color online) The figure depicts the choice of CUT for O and Zn in the case of ZnO. First we maximize the gap varying CUTO, then we vary CUTZn to reach a band gap value near the correct value.

Close modal

We have been naming our method37 of calculating band gaps as the LDA-1/2, though sometimes we use the Generalized Gradient Approximation, for instance of Ref. 18, instead of the Local-Density Approximation of Refs. 17 or.38 The NSCF version is that explained in the last paragraph of section III. In other words,

  1. The (model) number-density is chosen that of the ground state unperturbed by VS.

  2. For a single interaction we choose the potential as the sum of the ground state potential plus the CUT-trimmed self-energy potential.

  3. Solve the Schroedinger equation.

  4. Then CUT is chosen so as to make the band gap extreme.

The procedure works very well for highly ionic crystals, as shown in Fig 4. There we compare calculated band gaps with measured numbers. The pink balls are calculations according to the present procedure. They fall very near the straight line of perfection. Simple LDA, without any correction, gives band gaps smaller, as shown by the blue crosses. Then, using the procedure of the next subsection, that is the SCF LDA-1/2 we obtain the red squares that are half way between the pure LDA results and the NSCF LDA-1/2 results.

FIG. 4.

(color online) Comparison of calculated band gaps with experiment. The best results were obtained with the NSCF LDA-1/2 technique (pink balls). The worst were obtained with pure LDA (crosses). The standard LDA-1/2 technique37 (red squares), here being renamed SCF LDA-1/2, gives results between the two. The band structure calculations were made with the codes VASP39,40 and WIEN2k.41 

FIG. 4.

(color online) Comparison of calculated band gaps with experiment. The best results were obtained with the NSCF LDA-1/2 technique (pink balls). The worst were obtained with pure LDA (crosses). The standard LDA-1/2 technique37 (red squares), here being renamed SCF LDA-1/2, gives results between the two. The band structure calculations were made with the codes VASP39,40 and WIEN2k.41 

Close modal

This is the original LDA-1/2 method.37 The atomic self-energy potential, trimmed by the function Θ, Eq. (14), is added to the atoms of the crystal, as an external potential. Then self-consistent iterations are run. The parameter CUT is chosen according to make the gap extreme. Since the 2008 paper, many new semiconductors were calculated. The new plot of Calculated Band Gaps versus Experimental Band Gaps (Fig. 5) is very impressive, probably the best set of calculated band gaps available in the literature.

FIG. 5.

(color online) Comparison of calculated band gaps with experiment. The red square are the SCF LDA-1/2 (standard LDA-1/2). The crosses are standard LDA. The small gap semiconductors are metals (negative gaps), when calculated with LDA. LDA-1/2 corrects the situation. The band structure calculations were made with the codes VASP39,40 and WIEN2k.41 

FIG. 5.

(color online) Comparison of calculated band gaps with experiment. The red square are the SCF LDA-1/2 (standard LDA-1/2). The crosses are standard LDA. The small gap semiconductors are metals (negative gaps), when calculated with LDA. LDA-1/2 corrects the situation. The band structure calculations were made with the codes VASP39,40 and WIEN2k.41 

Close modal

It is very difficult to calculate the self-energy in the crystal itself and we have been recurring to atomic physics to calculate an approximation. We found that for most semiconductors the self-energy potential of the atom was transferable to crystals after trimming its 0.5/r tail. For the highly ionic crystal we found that no self-consistent loops could be made. For the crystals in the Fig. 6 we are finding that neither procedure, self-consistent or non self-consistent is very good. In the case of these oxides we found that a working procedure was to calculate the self-energy potential of Oxygen as

(17)

instead of Eq. (11) for the neutral atoms.

FIG. 6.

(color online) Comparison of calculated band gaps with experiment. The red squares are the SCF LDA-1/2 (standard LDA-1/2). The crosses are standard LDA, with no self-energy potential. The pink balls are NSCF LDA-1/2 calculations. The band structure calculations were made with the codes VASP39,40 and WIEN2k.41 

FIG. 6.

(color online) Comparison of calculated band gaps with experiment. The red squares are the SCF LDA-1/2 (standard LDA-1/2). The crosses are standard LDA, with no self-energy potential. The pink balls are NSCF LDA-1/2 calculations. The band structure calculations were made with the codes VASP39,40 and WIEN2k.41 

Close modal

We are paying a price for using an atomic self-energy potential instead of a self-energy potential defined in the crystal: we do not know a priori whether to run self-consistent iterations or NSCF. On the other hand, it is easy to verify that all covalent semiconductors and insulators require full SCF calculations, while the highly ionic insulators require treating the self-energy potential in a single NSCF iteration. Further, the localized valence band hole energy, evSv, falls within the band continuum for those cases when we make the standard SCF LDA-1/2, and falls outside for the cases when we use NSCF. Thus there is a simple explanation for the difference between the two procedures: when the hole energy falls within the band continuum the atomic self-energy potential is screened by the other band electrons and the self-consistent iterations must be completed to account for that screening. Contrarily, when the hole energy falls outside the band, no screening by the electrons of the other atoms must be considered and the situation is all similar to that of atoms.

In this paper we show the power of Slater's half-ion technique (transition state technique). Aside from its known usefulness to finite systems, such as atoms, molecules, and cluster, we now extend its use to crystals and semiconductors. The extension is attained by taking a second derivative of the energy functional, one step beyond Janak's first derivative, which allowed us to define a self-energy, means to calculate ionization energies and band gaps. The set of band gaps we are presenting have a GW precision,42 though it is obtained with minimal computational effort.

This paper resulted from our experience with the half-occupation technique after its proposal in Ref. [37]. In the case of atoms we compared our ionization energies with those of other methods and verified that we are performing better, except for the coupled cluster results of the first line of the Periodic Table. In the case of crystals we now included highly ionic insulators, for which our technique had to be modified. We are also presenting band gap results for many new covalent semiconductors. And to be complete, we found some solids (oxides) for which the technique did not perform very well.

The authors are thankful to the Brazilian Funding Agencies FAPESP and CNPq for continuous support. The authors are indebted to Profs. Marcio H. F. Bettega, Maria C. Santos and Mauro F. S. Ribeiro Jr. for their suggestions.

1.
J. C.
Slater
and
K. H.
Johnson
,
Phys. Rev. B
5
,
844
(
1972
).
2.
J. C.
Slater
,
Adv. Quantum Chem.
6
,
1
(
1972
).
3.
J. C.
Slater
and
J. H.
Wood
,
Int. J. Quant. Chem. Suppl.
4
,
3
(
1971
).
4.
J. F.
Janak
,
Phys. Rev. B
18
,
7165
(
1978
).
5.
Y.
Imamura
,
R.
Kobayashi
, and
H.
Nakai
,
J. Chem. Phys.
134
,
124113
(
2011
).
6.
D. H. E.
abd
E. R.
Johnson
,
X.
Hu
, and
W.
Yang
,
J. Phys. Chem. A
115
,
76
(
2011
).
7.
S.
Sanna
,
W. G.
Schmidt
,
T.
Frauenheim
, and
U.
Gerstmann
,
Phys. Rev. B
80
,
104120
(
2009
).
8.
J.-W.
Song
,
M. A.
Watson
, and
K.
Hirao
,
J. Chem. Phys.
131
,
144108
(
2009
).
9.
S.
Sanna
,
T.
Frauenheim
, and
U.
Gerstmann
,
Phys. Rev. B
78
,
085201
(
2008
).
10.
T.
Körzdörfer
and
S.
Kümmel
,
Phys. Rev. B
82
,
155206
(
2010
).
11.
E.
Kraisler
,
G.
Makov
, and
I.
Kelson
,
Phys. Rev. A
82
,
042516
(
2010
).
12.
E. R.
Johnson
,
W.
Yang
, and
E. R.
Davidson
,
J. Chem. Phys.
133
,
164107
(
2010
).
13.
F.
Gallino
,
G.
Pacchioni
, and
C. D.
Valentin
,
J. Chem. Phys.
133
,
144512
(
2010
).
14.
W.
Kohn
and
L. J.
Sham
,
Phys. Rev.
140
,
A1133
(
1965
).
15.
C.
Göransson
,
W.
Olovsson
, and
I. A.
Abrikosov
,
Phys. Rev. B
72
,
134203
(
2005
).
16.
J. R.
Leite
and
L. G.
Ferreira
,
Phys. Rev. A
3
,
1224
(
1971
).
17.
J. P.
Perdew
and
A.
Zunger
,
Phys. Rev. B
23
,
5048
(
1981
).
18.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
19.
D. A.
Liberman
,
Phys. Rev. B
62
,
6851
(
2000
).
20.
T.
Körzdörfer
,
J. Chem. Phys.
134
,
094111
(
2011
).
21.
I.
Dabo
,
A.
Ferretti
,
N.
Poilvert
,
Y.
Li
,
N.
Marzari
, and
M.
Cococcioni
,
Phys. Rev. B
82
,
115121
(
2010
).
22.
S. B.
Trickey
,
Phys. Rev. Letters
56
,
881
(
1986
).
23.
C. D.
Pemmaraju
,
T.
Archer
,
D.
Sánchez-Portal
, and
S.
Sanvito
,
Phys. Rev. B
75
,
045101
(
2007
).
24.
J. M.
García-Lastra
,
P. L.
Cook
,
F. J.
Himpsel
, and
A.
Rubio
,
J. Chem. Phys.
133
,
151103
(
2010
).
25.
R. R.
Zope
,
T.
Baruah
,
S. L.
Richardson
,
M. R.
Pederson
, and
B. I.
Dunlap
,
J. Chem. Phys.
133
,
034301
(
2010
).
26.
In the table below the Electrostatic self-energy is the first term of Eq. (5), Self-energy is calculated according to Eq. (6). Entries are in Rydbergs.
27.
We rewrite Eq. (3) as
where H is the part of the Eq. operating on nα. Then, as ψα is a normalized eigenfuntion of H, the derivative with respect to any of its parameters is null.
28.
T.
Körzdörfer
,
S.
Kümmel
,
N.
Marom
, and
L.
Kronik
,
Phys. Rev. B
79
,
201205
R
(
2009
).
29.
D. M.
Ceperley
and
B. J.
Alder
,
Phys. Rev. Lett.
45
,
566
(
1980
).
30.
S.
Froyen
,
N.
Troullier
,
J. L.
Martins
, and
A.
Garcia
, Code ATOM inside the SIESTA package (
2006
).
31.
cccbdb.nist.gov, Calculated data D.6, Ionization Energies.
32.
W.
Klopper
,
R. A.
Bachorz
,
D. P.
Tew
, and
C.
Hättig
,
Phys. Rev. A
81
,
022503
(
2010
).
33.
The reader should not confuse localized quasi-stationary states in crystals, such as the holes we are talking about, with localized basis functions as the Wannier functions.
34.
J. P.
Perdew
and
M.
Levy
,
Phys. Rev. Lett.
51
,
1884
(
1983
).
35.
L. J.
Sham
and
M.
Schlüter
,
Phys. Rev. Lett.
51
,
1888
(
1983
).
36.
N.
Helbig
,
N. N.
Lathiotakis
, and
E. K. U.
Gross
,
Phys. Rev. A
79
,
022504
(
2009
).
37.
L. G.
Ferreira
,
M.
Marques
, and
L. K.
Teles
,
Phys. Rev. B
78
,
125116
(
2008
).
38.
J. P.
Perdew
and
Y.
Wang
,
Phys. Rev. B
45
,
13244
(
1992
).
39.
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
47
,
R558
(
1993
).
40.
G.
Kresse
and
J.
Furthmüller
,
Comput. Mater. Sci.
6
,
15
(
1996
).
41.
P.
Blaha
,
K.
Schwarz
,
P.
Sorantin
, and
S. B.
Trickey
,
Comput. Phys. Commun.
59
,
399
(
1990
), see www.wien2k.at.
42.
G.
Onida
,
L.
Reining
, and
A.
Rubio
,
Rev. Mod. Phys.
74
,
601
(
2002
).
This article is distributed under a Creative Commons Attribution 3.0 Unported License.