Conceptual and computational advances in multiple-scattering electronic-structure calculations

COMPUTATIONAL MATERIALS SCIENCE ELSEVIER

Computational Materials Science 10 (1998) 373-380

Conceptual and computational advances in multiple-scattering electronic-structure calculations Rudolf Zeller ’ Institutfir Festkiitpeflorschung,

Forschungszentrum Jiilich GmbH, D-S2425 Jiilich, Federal Republic of Germany

Abstract

A physically transparent transformation of the Korringa-Kohn-Rostoker (KKR or multiple-scattering) method into a tightbinding form is described. The transformation replaces the complicated, slowly decaying, traditional KKR structure constants by exponentially decaying “tight-binding” parameters. The main computational effort consists in the inversion of sparse matrices and scales for surfaces and interfaces, i.e. for systems with two-dimensional periodicity, linearly with the number of layers. This gives the opportunity to treat high-indexed surfaces as an approximation for almost isolated surface steps. Additional adatoms on surfaces and at steps can also be treated and it is discussed that reliable atomic forces and geometric arrangements can be obtained. Copyright 0 1998 Elsevier Science B.V. PACS: 70.15.-m;

71.20.-b

1. Introduction Since ab initio electronic-structure

calculations

require too much computing

power except for rather simple

modeling systems, the method of empirical potentials is widely used for the computational

simulation

and analysis

of materials. One step forward to bridge the gap between quantum-mechanical electronic-structure calculations and the method of empirical potentials is to increase the efficiency of density-functional calculations for larger systems. The search for efficient methods, particularly for ones, which need a computational effort scaling linearly with the system size [ 11, has recently received much attention. Most of these methods gain their efficiency by the assumption that the electronic wave functions decay fast enough in space so that a local region around each individual atom can be treated separately. This allows to treat insulating or semiconducting materials, but probably not metallic systems, where long-ranged electronic states at the Fermi energy are decisive for many physical properties. Close-packed metallic systems, on the other hand, are well treated by the multiple-scattering (KKR) method, which originated in the works of Lord Rayleigh [2], Kasterin [3] and Ewald [4] and was formulated for the Schriidinger equation by Korringa, Kohn and Rostoker [5]. The main historical disadvantages of the KKR method are the ’ Tel.: +49 2461-615192, fax: +49 2461-612620, e-mail: [email protected]. 0927-0256/98/$19,00 Copyright 0 1998 Elsevier Science B.V. All rights reserved PII SO927-0256(97)00105-5

314

R. iWler/Computational

Materials Science 10 (1998) 373-380

numerical difficulties connected with the long-range and complicated energy and wave vector dependencies of the so called KKK structure constants, and the original restriction to spherically symmetric, nonoverlapping potentials around each atom, the so-called muffin-tin approximation. The controversy whether space-filling potentials can be used has now been settled and it is commonly assumed that the KKK method can generally be applied [6,7]. This general applicability is also substantiated by calculations for atomic forces and displacements [8], which will be discussed in Section 6. The main issue of this paper concerns a recent development which bypasses the numerical complexity of the traditional multiple-scattering equations and allows to determine, for instance, the electronic structure of surfaces with an effort, which increases only linearly with the number of surface layers. This opens the way to treat highindexed surfaces and in this approximation isolated surface steps. An important feature of the KKK method is that it allows to determine the Green function for the considered material. Thus the treatment of local perturbations like impurities, adatoms on surfaces, at steps or even at kinks is straightforward [9,10]. The recent progress which drastically reduces the computational effort and algorithmic complexity of the traditional KKK method has its origin and motivation in the success of the tight-binding linear-muffin-tin-orbital (TB-LMTO) method of Andersen and Jepsen [ 111. The concept of the TB-LMTO screening constants was first extended to energy-dependent screening parameters by Andersen et al. [12] as a tool to obtain a TB form of the KKK method. Since the original determination of these screening parameters was numerically difficult, a simplified procedure was suggested by Szunyogh et al. [ 131,but the real breakthrough came by the discovery that the screening parameters can be obtained in a physically transparent and easy way. For the transparent determination and the subsequent use of the screening parameters, two different, but related techniques have been developed. One technique [ 141is based on the concept of unitary spherical waves, which are defined as solutions for a hard sphere solid, and is described in a wave function language. The other technique [ 151 is based on a repulsive reference system, which consists of constant, finite potentials inside nonoverlapping spheres, and uses a Green-function language. The Green-function technique will be presented here and was used to obtain the results discussed in Section 4.

2. Theory The theory has been described in previous publications [15-171 and only a short account is given here. The Green functions G(r, r’, E) and G’(r, r’, E) for two different Kohn-Sham effective potentials V(r) and V’(r) are connected by a Dyson equation: G(r, r’, E) = G’(r, r’, E) f

s

G’(r, r”, E)[V(r”) - V’(r”)]G(r”,

r’, E) dr”,

(1)

which by the multiple-scattering expression: G(r + R”, r’ + Rn’, E) = Pn’Gs(r + Rn, r’ + Rn, E) + c

RF(r, E)G’,$,(E)Rs(r’,

E)

(2)

LL’ can

be

transformed [6] into a system of linear equations:

G?;,(E)

= G;;’ (E) + F> r G;:;“(E) II” L”

~]t;&)

- tZ.~~,,,(E)lG”,‘~~,(E).

(3)

L”’

Cell-centered coordinates r and r’ are used in (2) and E and L = (e, m) denote energy and angular-momentum numbers. Ri and G, are wave functions and Green function [6] for a single potential restricted to the Voronoi

R. Zeller/Computational Materials Science 10 (1998) 373-380

315

cell around the atomic center R”. The r-matrices tzL, and t;F, are also single-site quantities [6] determined by the potentials V and V’ restricted to the cell at R”. The traditional reference system is empty space, denoted by 0 in the following. In empty space the potentiaf V” and consequently the t-matrix $L, vanish and the so-called structural Green-function matrix elements G:i: are analytically known as

withAR=R’-Rn’andC LLTL”= I& dS;Z,YL(r)YLT(r)YL” (r), where hr denotes spherical Hankel functions and YL spherical harmonics. Empty space seems to be the easiest reference system, but this is not true because it makes the calculations difficult and time consuming. For instance, in periodic geometries one needs Fourier sums like Gt,,(k,

E) = xexp(ikR” .’

- ikR”‘)G:‘:‘(E),

(5)

which converge only conditionally and require complicated Ewald procedures [4,18] for their evaluation as a consequence of the slow decaying Hankel functions in (4). This problem is connected with the well-known freeelectron singularities in the traditional KKR structure constants GiL, (k, E), which occur if the free-electron-band structure condition E = (k + Kh I2 is satisfied, where k denotes the wave vector and Kh a reciprocal lattice vector. The problem of the traditional reference system thus arises from the eigenstates of empty space. When the basic equations (1) and (3) are solved, these eigenstates are eliminated in the solution process and replaced by the eigenstates of the material, for which the Green function is calculated. It would be much more advantageous to start with a reference system, which has no eigenstates in the energy range, which is used by density-functional electronic-structure calculations and typically covers energies below 1 Ry. Such a reference system consists of an infinite array of repulsive, constant, finite potentials within nonoverlapping spheres around each center R” and a zero potential in the interstitial region between the spheres [ 151. From firstorder perturbation theory one expects that the eigenstates in the repulsive reference system have higher energies than in empty space and that the energy shift is given by the product of the potential strength and the volume filling. Calculated density-of-states for the reference system are shown in Fig. 1 and confirm the shift to higher energies. Contrary to first-order perturbation theory the shift is not uniform and also shows a saturation effect. The lowest eigenstate is at 1.35 Ry for potentials of 2 Ry height, at 2.25 Ry for potentials of 4 Ry height and just below 3 Ry for potentials of 20 Ry height. In Fig. 1 it is important to notice that the density of states vanishes below the energy of the lowest eigenstate. The consequences of this fact will be discussed in the next section. It is also important to notice that the vanishing density of states has been achieved by using an angular-momentum cut-off I,, = 3 in all equations. This means that the repulsive reference system already works if only low angular-momentum components are used in the calculations. Thus no convergence problems with respect to angular-momentum sums occur in the TB form of the KKK method.

3. Rasic facts Fig. 1 shows that for lower energies no eigenstates exist in the reference system and thus any local perturbation of the reference system must decay exponentially. This means that the Green function G’(r, r’, E) and the structural Green-function matrix elements G>’ (E) decrease exponentially with increasing distance Ar = ]r - r’] or AR = IR” - Rn’ 1.The exponential decay was verified by numerical calculations [ 12,151 which showed that the matrix

R. Zeller/Computational Materials Science 10 (1998) 373-380

376

0

1

2

Energy CRY > Fig. 1. Density-of-states (DOS) for a fee reference system. The height of the repulsive potentials is chosen as 0, 1, 2, 4 Ry and the corresponding DOS are plotted with dotted, solid, dashed and dash dotted curves. The Brillouin sampling was done with 11726special points in the irreducible part and smoothed by a temperatute broadening with T = 800 K.

elements G>Z’f (E) decrease by a factor of about 10 000, if the distance AR is increased by the amount of a lattice constant. Because of this fast decrease it was speculated [ 12,151 that the matrix elements can be neglected beyond a rather short coupling range. This speculation was verified by density-of-states calculations for empty space [ 16,171, by total-energy calculations for the fee metals Al, Cu and Pd [ 161 and by density-of-states surface [ 171. Some of these results will be discussed in Section 4. The short coupling

range has important

structural Green function

matrix elements

calculations

for the Cu

consequences for the numerical effort, both for the evaluation of the (E) of the reference system and for their use in (3). The matrix

G;”

elements can be evaluated from the ones of empty space, defined in (4), by an equation, which similar to (3) is given by G ;;‘(E)

= Go;;“‘(E)

+ r x G~;;:“(E)I;;““(E)G;“;‘(E). II” L”

(6)

Here Go, plays the role of the reference Green function and G’ is to be calculated. Because of the exponential decay of G;“L (E) the sum over n” can be restricted to a finite number of centers around n’. Thus (6) can be solved in real space for each center separately. This task is ideally suited for parallel computing since the double loop over the individual atoms and over the different values of E can be completely distributed over the computing nodes. The subsequent evaluation of (3) is simplified by the sparse matrix structure of G$( E), which is schematically given by

(7)

R. Zeller/Computational

Materials Science IO (1998) 373-380

377

where each x denotes a nonvanishing matrix block. For systems, which are periodic in two dimensions like surfaces, interfaces or multilayers, a two-dimensional lattice Fourier transform can be applied to (3) and the blocks in (7) are matrices in the angular-momentum indices. For general three-dimensional systems, each block x has again the structure shown in (7) with blocks x’, which again have the same structure with blocks x”, which are matrices in the angular-momentum indices. For surfaces, interfaces and multilayers the almost block tridiagonal structure (7) can efficiently be treated [ 171. The general case is more difficult, existing algorithms for general sparse matrices usually do not exploit the regular structure of (7) and are not tailored to obtain the blocks on the diagonal of the inverse of the sparse matrix. These blocks, essentially the on-site (n = n’) Green function matrix elements G”,n,,(E), are needed to determine the electronic density, the basic quantity of density-functional calculations. For that, a contour integration with about 20-40 complex energies is applied [21]. In this respect the TB-KKR method also differs from standard TB electronic-structure methods, where a TB Hamiltonian is used and the eigenvalue problem H W = E 9 must be solved. In the TB-KKR method the diagonal blocks of the inverse of the matrix M = 1 - G’(t - t’) are obtained by solving systems of linear equations. Another important difference concerns the coupling range. In usual TB methods a short coupling range is assumed in the considered material, whereas in the TB-KKR method the short range of TB parameters Gi’ (E) is an exact property of the reference system and not of the material, for which the electronic-structure calculations are intended. This makes the TB-KKR method equally well suited for insulating, semiconducting and metallic systems.

4. Accuracy

For the numerical accuracy of the TB-KKR method, it must be investigated how many TB parameters G:;‘(E) must evaluated in (6) and used in (3). Such investigations are described in previous publications [ 16,171and showed that finite real space clusters of repulsive potentials can be used to determine the TB parameters. The clusters were constructed around the central site n with a repulsive potential on this site and repulsive potentials on several shells of neighboring sites. Both indices n’ and n” in (6) were restricted to the sites in the cluster and the calculated TB parameters G;z’ (E) from the central site n to neighboring sites n’ were used in (3). If enough repulsive potentials were used, densities of states could be calculated very well almost up to the energy of the lowest eigenstates of the reference system. Already only 0.1 Ry below this energy the plotted densities of states [ 16,171cannot be distinguished from the exact ones obtained by standard KKR calculations. Results for equilibrium total energies, lattice constants and bulk moduli [ 161for the fee metals Al, Cu and Pd are shown in Fig. 2 as functions LL, (E). The convergence with respect of the number of repulsive potentials used to determine the TB parameters G’*““’ to the number of repulsive potentials is always rapid. Typical errors for total energies, lattice constants and bulk moduli are 1 mRy, 0.2 pm and 1 GPa, if 19 potentials on central, nearest and next-nearest neighbor sites are used, and 5 p Ry, 0.01 pm and 0.05 GPa, if 79 potentials are used. The calculations applied full-potential KKR codes [ 191 and used an angular-momentum expansion for single-site charge densities and potentials up to 1, = 6. They were nonrelativistic and the exchange-correlation potential as given by Vosko et al. [20] was used. The k point Brillouin sampling was facilitated by a temperature of T = 800 K as described in [ 16,2 11.

5. Effwiency for large-scale calculations

The main computational tasks of multiple-scattering theory consist in the calculation of single-site quantities like Ri(r, E) and &,(E) and of the structural Green-function matrix elements G?;,(E). All quantities can be

378

R. Zeller/Computational

0

50

Materials Science 10 (1998) 373-380

100

150

Fig. 2. Deviations dE, da, dB from standard KKR results for total energies (top), lattice constants (middle) and bulk moduli (bottom) for the fee metals Al, Cu and Pd as function of the number of repulsive potentials used to determine the TFJ parameters.

calculated in parallel for each energy E. A second parallelization strategy concerns the calculation of the singlesite quantities independently for each atom and the distribution of the inversion of the KKR matrix M onto the computing nodes by standard software [24]. The TB-KKR method uses matrices of economic size. With dimensions I6N, where N denotes the number of inequivalent atoms, good quantitative all-electron results are obtained even for transition metals. This matrix size should be compared to the one used in O(N) methods [ 11, where N denotes the number of electrons and where the prefactor is usually larger. Contrary to the traditional KKR matrix, the TB-KKR matrix is sparse for large systems and the TB-KKR matrix elements can be calculated much easier. Fourier transformations similar to (5) are straightforward, thus bypassing

R. iWler/Computational

any complicated

Ewald sums of the traditional

exploited for systems with two-dimensional layer it is achieved that only neighboring numerical

complexity

Materials Science 10 (1998) 373-380

379

KKR method. The sparsity of the TB-KKR

periodicity principal

like surfaces. By combining

matrix can directly be

several layers into a principal

layers interact in the reference system. As a consequence,

the

is O(N) for N different and O(log N) for N identical principal layers. Such fast calculations

allow to treat complicated

surfaces, e.g. high-indexed

surfaces with a periodic repetition of surface steps.

6. Usefulness for perturbed systems Since multiple-scattering material, the electronic Green function

method in its standard or TB form provides the full Green function for the considered

structure of locally perturbed

technique

[9]. Thus the electronic

adatoms on surfaces can be determined.

materials can straightforwardly

structure for impurities

be calculated

by the KKR-

in the bulk, in the surface layers or for

This is also true for adatoms and substitutional

atoms at steps and even at

kinks, which can be approximated by missing or additional chains of a few adatoms. For perturbed materials the atoms are displaced from the ideal positions of the underlying periodic lattice, and it is satisfying that the displacements can now be calculated by the KKR method [8]. Reliable forces and displacements around the perturbing atoms are obtained by an extension of the full-potential formalism and by the application of a modified, “ionic” HellmannFeynman theorem, which allows to treat the core electrons in spherical approximation. atoms is based on a reexpansion transformation

c = UGU-‘.

of the structural Green function

The transformation

The treatment of the displaced

matrix around the shifted positions by a matrix

matrix U depends on the shifts sn and is given by

where jl denotes spherical Bessel functions. The reliability of the calculated forces has been checked by calculating them in two different ways, either by the ionic Hellmann-Feynman theorem or by numerically differentiating total energies. It was found that these forces agreed very well, if the high lying core states are carefully treated. For instance, for a Ti impurity in Cu, a spherical, atomic-like treatment of the 3p core states overestimated the forces by 50%. Both the full potential anisotropy

and the hybridization

with the valence states were found to be important

effects for calculating reliable forces. A spherical-potential approximation around the atoms must completely fail as numerical tests have shown [22], which overestimated the forces typically by a factor of 10. The reliability of the calculated forces and resulting displacements, obtained by zero-force condition, the experimental information available from extended-X-ray-absorption-fine-structure impurities

in metals [23]. The measured nearest-neighbor

[S]. This agreement and the sensitivity full-potential

displacements

can nicely be checked against (EXAFS) measurements for

agree very well with the calculated

ones

of the calculated forces represent a stringent test for the accuracy of the used

KKR formalism.

7. Summary The traditional multiple-scattering (KKR) method can be transformed into a TB form by the concept of a reference system with repulsive, nonoverlapping, spherical potentials. The TB parameters can be calculated with finite clusters of repulsive potentials and used for accurate density-functional calculations. Contrary to the traditional KKR matrix, the TB-KKR matrix is sparse for large systems and can be easily calculated. The sparsity can be exploited and leads to O(N) algorithms for surface calculations which are expected to be of great use for complicated surfaces and surface steps. The TB-KKR method allows to determine the Green function of the considered materials and thus

380

R. Zeller/Computational

the treatment of local perturbations.

Materials

Recent improvements

Science 10 (1998) 373-380

in full-potential

calculations

show

within the HCM Network Contract No. ERBCHRXCT930369

and

that reliable atomic forces and resulting geometric atomic arrangements

The work has benefited from collaborations

multiple-scattering

can be calculated.

the TMR Network Contract No. EMRX-CT96-0089.

References [ l] M. Krajci and J. Hafner, Phys. Rev. Lett. 74 (1995) 5 100; Y. Wang, GM. Stocks, W.A. Shelton, D.M.C. Nicholson, Z. Szotek, and W. M. Temmerman, Phys. Rev. Lett. 75 (1995) 2867; LA. Abrikosov, A.M.N. Niklasson, S.I. Simak, B. Johansson, A.V. Ruban and H.L. Skriver, Phys. Rev. Lett. 76 (1996) 4203; S. Wei and M.Y. Chou, Phys. Rev. Lett. 76 (1996) 2650. T. Zhu, W. Pan and W. Yang, Phys. Rev. B 53 (1996) 12 7 13; S. Goedecker and L. Colombo, Phys. Rev. Lett. 73 (1994) 122; A.F. Voter, J.D. Kress and R.N. Silver, Phys. Rev. B 53 (1996) 12 733; A.P. Horsfield, A.M. Bratkovsky, M. Fearn, D.G. Pettifor and M. Aoki, Phys. Rev. B 53 (1996) 12694; S. Itoh, P. Ordejon, D.A. Drabold and R.M. Martin, Phys. Rev. B 53 (1996) 2132; E. Hemandez, M.J. Gillan and CM. Goringe, Phys. Rev. B 53 (1996) 7147; A. Canning, G. Galli, F. Mauri, A. De Vita and R. Car, Comp. Phys. Commun. 94 (1996) 89. [2] Lord Rayleigh, Philos. Mag. 34 (I 892) 48 I. [3] N.P. Kasterin, Diss. Moscow (1903). 141 P.P.Ewald,Ann.Physik49(1916) 1;49(1916) 117;64(1921)253; [5] J. Korringa, Physica 13 (1947) 392; W. Kohn and N. Rostoker. Phys. Rev. 94 (I 954) 1 I I I. [6] R. Zeller, J. Phys. C 20 (1987)2347. [7] R.G. Newton, Phys. Rev. Lett. 65 (1990) 203 I; J. Math. Phys. 33 ( 1992) 44; W.H. Butler, A. Gonis and X.-G. Zhang, Phys. Rev. B 45 (1992) 11 527; 48 (1993) 2 I 18: S. Bei der Kellen, Y. Oh, E. Badralexe and A.J. Freeman, Phys. Rev. B 51 (1995) 9560. [8] N. Papanikolaou, R. Zeher, P.H. Dederichs and N. Stefanou, Phys. Rev. B 55 (1997) 4157. [9] R. Zeller and P. H. Dederichs, Phys. Rev. Lett. 42 (1979) I7 13: P.J. Braspenning, R. Zeller, A. Lodder and P.H. Dederichs, Phys. Rev. B 29 (1984) 703. [IO] K. Wildberger, V.S. Stepanyuk, P. Lang, R. Zeller and P.H. Dederichs, Phys. Rev. Lett. 75 ( 1995) 509. [ 1I] O.K. Andersen and 0. Jepsen, Phys. Rev. Lett. 53 (1984) 2571. [ 121 O.K. Andersen, A.V. Postnikov and S. Yu Savrasov, in: Application of Multiple Scattering Theory to Materials Science, eds. W.H. Butler, P.H. Dederichs, A. Gonis and R.L. Weaver, MRS Symposia Proceedings No. 253 (Materials Research Society, Pittsburgh. 1992). [ 131 L. Szunyogh, B. Ujfalussy, P. Weinberger and J. Kollar, Phys. Rev. B 49 (1994) 272 I. [ 141 O.K. Andersen, 0. Jepsen and G. Krier, in: Lectures on Methods of Electronic Structure Calculations, eds. V. Kumar, O.K. Andersen and A. Mookerjee (World Scientific, Singapore, 1994). [ 151 R. Zeller, P.H. Dederichs, B. Ujfalussy. L. Szunyogh and P. Weinberger. Phys. Rev. B 52 (1995) 8807. [ 161 R. Zeller, Phys. Rev. B 55 (1997) 9400. [ 171 K. Wildberger, R. Zeller and P.H. Dederichs, Phys. Rev. B 55 ( 1997) IO 057. [18] ES. Ham and B. Segall, Phys. Rev. 124 (1961) 1786. [ 191 B. Drittler, M. Weinert, R. Zeller and P.H. Dederichs, Solid State Comm. 79 (1991) 3 I. [20] S.H. Vosko, L. Wilk and M. Nusair, Can. J. Phys. 58 (1980) 1200. [21] K. Wildberger, P. Lang, R. Zeller and P.H. Dederichs, Phys. Rev. B 52 (1995) II 502. [22] K. Abraham, B. Drittler, R. Zeller and PH. Dederichs, KFA-JUL-Report, Jtil-2451 (1991). [23] U. Scheuer and B. Lengeler, Phys. Rev. B 44 (1991) 9883. [24] J. Choi, J. Demmel, I. Dhillon, J. Dongarra, S. Ostrouchov, A. Petitet. K. Stanley, D. Walker and R.C. Whaley, Comp. Phys. Comm. 97 (1996) I.