DFT calculations of solids with LAPW and WIEN2k

ry 17 s w z S echn ed fo on t ic bo ility, pond roble s. O plo uch c for systems containing about 100 atoms per unit cell. A selection of representative examples and the references to the original r 2003 Elsevier Inc. All rights reserved. comes to the nanometer (nm) scale or atomic dimen- of the key advances. Other applications are found in the area of magnetic recording or optical storage media. A experimental data such as equilibrium geometries, bulk understood, one often must rely on ab initio calcula- tions. They are more demanding in terms of computer ARTICLE IN PRESS typical example in chemistry is heterogeneous catalysis, in which one likes to understand the details of catalytic resources and thus allow only the treatment of cell units smaller than force-field calculations. The advantage of first-principle (ab initio) methods lies in the fact that they can be carried out without knowing any experi- mental data of the system. The following presentation �Fax: +43-1-58801-15698. E-mail address: kschwarz@theochem.tuwien.ac.at. 0022-4596/03/$ - see front matter r 2003 Elsevier Inc. All rights reserved. doi:10.1016 sions measured in A˚, the properties are determined by the electronic structure of the solid. In the development of modern materials an understanding of the atomic scale is frequently essential in order to replace trial and error procedures by a systematic materials design. Modern devices in the electronic industry provide such an example, where the increased miniaturization is one moduli or special vibrational frequencies (phonons). Such schemes have reached a high level of sophistication and are useful in their range of applicability, namely within a given class of materials provided good parameters are already known from closely related systems. If, however, such parameters are not available, or if a system shows unusual phenomena that are not yet Keywords: DFT; LAPW; Energy bands; WIEN2k; Electronic structure 1. Introduction Solid materials are of great technological interest. Different materials are governed by very different length and time scales. They may differ by many orders of magnitude depending on their applications. Let us focus on the length scale, where from meters (m) to micro- meters (mm) classical mechanics and continuum models are the dominating concepts to investigate the properties of the corresponding materials. However, when one processes between molecules interacting with a solid surface, e.g., in zeolites. One way of studying complex systems is to perform computer simulations. Calculation of solids in general (metals, insulators, semiconductors, minerals, etc.) can be performed with a variety of methods from classical to quantum mechanical (QM) approaches. The former are mainly force field schemes, in which the forces that determine the interactions between the atoms are parameterized in order to reproduce a series of literature is given. Journal of Solid State Chemist DFT calculations of solid Karlhein Institute of Materials Chemistry, Vienna University of T Received 16 January 2003; received in revis Abstract In solids one often starts with an ideal crystal that is studied several atoms (at certain positions) and is repeated with period structure that is responsible for properties such as relative stab electrical, mechanical, optical or magnetic behavior, etc. Corres functional theory (DFT), according to which the many-body p one-electron equations, the so-called Kohn–Sham (KS) equation the linearized-augmented-plane-wave (LAPW) method that is em properties on the atomic scale (see www.wien2k.at). Nowadays s /S0022-4596(03)00213-5 6 (2003) 319–328 ith LAPW and WIEN2k chwarz� ology, Getreidemarkt 9/165-TC, A1060 Vienna, Austria rm 28 March 2003; accepted 10 April 2003 he atomic scale at zero temperature. The unit cell may contain undary conditions. Quantum mechanics governs the electronic chemical bonding, relaxation of the atoms, phase transitions, ing first principles calculations are mainly done within density m of interacting electrons and nuclei is mapped to a series of ne among the most precise schemes to solve the KS equations is yed for example in the computer code WIEN2k to study crystal alculations can be done—on sufficiently powerful computers— treated exactly but correlation effects are ignored by ARTICLE IN PRESS tate C definition. The latter can be included by more sophis- ticated approaches (e.g., a linear combination of Slater determinants) such as the configuration interaction (CI) scheme that progressively requires more computer time with a scaling as bad as N7 when the system size (N) grows. As a consequence it is only feasible to study small systems, which contain a few atoms. An alternative scheme is the DFT that is commonly used to calculate the electronic structure of complex systems containing many atoms such as large molecules or solids [1,2]. It is based on the electron density rather than on the wave functions and treats exchange and correlation, but both approximately. Since it became the method of choice for solids it will be described in more detail below. The ideal crystal is defined by the unit cell, which may contain many atoms and is repeated indefinitely according to translational symmetry. Periodic boundary conditions are used to describe the infinite crystal by knowing the properties in one unit cell. The additional symmetry operations (inversion, rotation, mirror planes, etc.) that leave the ideal crystal invariant allow both providing symmetry labels and simplifying the calcula- tions. Each ideal structure belongs to one of the 230 space groups that are characterized in the International Tables [3]. Only the electronic ground state is studied. When thermal vibrations are considered in the form of phonon calculations [4] an adiabatic (Born Oppenhei- mer) approximation in which the electronic degrees of freedom are decoupled from the nuclear motions is used. 2. Density functional theory (DFT) 2.1. The Kohn–Sham equations The well-established scheme to calculate electronic properties of solids is based on the DFT, for which Walter Kohn has received the Nobel Prize in chemistry in 1998. DFT is a universal approach to the quantum mechanical many-body problem, where the system of interacting electrons is mapped in a unique manner onto an effective non-interacting system that has the same will be restricted to ab initio methods whose main characteristics shall be briefly sketched. The fact that electrons are indistinguishable and are fermions requires that their wave functions be anti- symmetric when two electrons are interchanged. This situation leads to the phenomenon of exchange. Starting with molecules there are mainly two types of approaches for a full quantum mechanical treatment, HF and DFT. The traditional scheme is (or was) the Hartree–Fock (HF) method which is based on a wave function description (with one Slater determinant). Exchange is K. Schwarz / Journal of Solid S320 total density. Hohenberg and Kohn [1] have shown that the ground state electron density r (in atoms, molecules or solids) uniquely defines the total energy E; i.e., E½r� must be a functional of the density. Thus one does not need to know the many-body wave function. The non- interacting particles of this auxiliary system move in an effective local one-particle potential, which consists of a classical mean-field (Hartree) part and an exchange- correlation part Vxc (due to quantum mechanics) that, in principle, incorporates all correlation effects exactly. According to the variational principle a set of effective one-particle Schro¨dinger equations, the so-called Kohn– Sham (KS) equations [2], must be solved. Its form is ½�r2 þ Vextð~rrÞ þ VC½rð~rrÞ� þ Vxc½rð~rrÞ��Fið~rrÞ ¼ eiFið~rrÞ ð1Þ when written in Rydberg atomic units for an atom with the obvious generalization to molecules and solids. The four terms represent the kinetic energy operator, the external potential from the nucleus, the Coulomb-, and exchange-correlation potential, VC and Vxc: The KS equations must be solved iteratively till self-consistency is reached. The iteration cycles are needed because of the interdependence between orbitals and potential. In the KS scheme the electron density is obtained by summing over all occupied states, i.e., by filling the KS orbitals (with increasing energy) according to the aufbau principle. rð~rrÞ ¼ Xocc i ½fið~rrÞ�2: ð2Þ From the electron density the VC and Vxc potentials for the next iteration can be calculated, which define the KS orbitals. This closes the SCF loop. The exact functional form of the potential Vxc is not known and thus one needs to make approximations. Early applications were done by using results from quantum Monte Carlo calculations for the homogeneous electron gas, for which the problem of exchange and correlation can be solved exactly, leading to the original local density approximation (LDA). LDA works reasonably well but has some shortcomings mostly due to the tendency of overbinding, which cause e.g., too small lattice con- stants. Modern versions of DFT, especially those using the generalized gradient approximation (GGA), im- proved the LDA by adding gradient terms of the electron density and reached (almost) chemical accu- racy, as for example the version by Perdew, Burke, Ernzerhof (PBE) [5]. In the study of large systems the strategy differs for schemes based on HF or DFT. In HF based methods the Hamiltonian is well defined but can be solved only approximately (e.g., due to limited basis sets). In DFT, however, one must first choose the functional that is used to represent the exchange and correlation effects hemistry 176 (2003) 319–328 (or approximations to them) but then one can solve this rather than a pseudo-potential approach with unphysi- ARTICLE IN PRESS tate C effective Hamiltonian almost exactly, i.e., with very high precision. Thus in both cases an approximation enters (either in the first or second step) but the sequence is reversed. This perspective illustrates the importance in DFT calculations of improving the functional, since this defines the quality of the calculation. 2.2. Solving the DFT equation, the choice of basis sets and wave functions Many computer programs that can solve the DFT equations are available but they differ in the basis sets. Many use an LCAO (linear combination of atomic orbitals) scheme in one form or another. Some use Gaussian or Slater type orbitals (GTOs or STOs), others use plane wave (PW) basis sets with or without augmentations, and some others make use of muffin tin orbitals (MTOs) as in linear combination of MTOs (LMTO) or augmented spherical wave (ASW). In the former schemes the basis functions are given in analytic form, but in the latter the radial wave functions are obtained by numerically integrating the radial Schro¨- dinger equation, whereas the angular dependence is given analytically. Closely related to the basis set used is the explicit form of the wave functions, which can be well represented by them. These can be nodeless pseudo-wave functions or all-electron wave functions including the complete radial nodal structure and a proper description close to the nucleus. 2.3. The form of the potential In the muffin tin or the atomic sphere approximation (MTA or ASA) an atomic sphere, in which the potential (and charge density) is assumed to be spherically symmetric, surrounds each atom in the crystal. While these schemes work reasonably well in highly coordi- nated, closely packed systems (as for example face centered cubic metals) they become very approximate in all non-isotropic cases (e.g., layered compounds, semi- conductors, or open structures). Schemes that make no shape approximation in the form of the potential are termed full-potential schemes (see Section 3.2). With a proper choice of pseudo-potential one can focus on the valence electrons, which are relevant for chemical bonding, and replace the inner part of their wave functions by a nodeless pseudo-function that can be expanded in PWs with good convergence. 2.4. Relativistic effects If a solid contains only light elements, non-relativistic calculations are well justified, but as soon as a system of interest contains heavier elements, relativistic effects can K. Schwarz / Journal of Solid S no longer be neglected. In the medium range of atomic cal wave functions near the nucleus. On the other hand for an efficient optimization of a structure, in which the shape (and symmetry) of the unit cell changes, it is very helpful to know the corresponding stress tensor. These tensors are much easier to obtain in pseudo-potential schemes and thus are available there. In augmentation schemes, however, such algorithms become more tedious and consequently are often not implemented. On the other hand all-electron methods do not depend on choices of pseudo-potentials and contain the full wave function information. Thus, the choice of method for a particular application depends on the properties of interest and may affect the accuracy, ease or difficulty to calculate them. 3. The full-potential linearized augmented plane wave (LAPW) method One among the most precise schemes for solving the Kohn–Sham equations is the full-potential linearized augmented plane wave (FP-LAPW) method (see e.g., [8]). There are several programs employing this method such as FLAPW (Freeman’s group), FLEUR (Blu¨gel’s group), D. Singh’s code and others. Here we focus on the WIEN code that has been developed in our group during the last two decades and is used worldwide by more than 500 groups coming from universities and industrial laboratories. The original version (WIEN) was the first LAPW code that was published [9] and thus was made available for other users. 3.1. The LAPW method In the LAPW method [8] the unit cell is partitioned numbers (up to about 54) the so-called scalar relativistic schemes [6] are often used, which describe the main contraction or expansion of various orbitals (due to the Darwin s-shift or the mass-velocity term) but omit spin– orbit splitting. This version is computationally easy and thus is highly recommended for all systems. The spin– orbit part can be included in a second-variational treatment [7]. For very heavy elements it may be necessary to add p1 2 orbitals [41] or to solve Dirac’s equation, which has all these terms included. 2.5. Method of choice and properties As a consequence of the aspects described above, different methods have their advantages or disadvan- tages when it comes to computing various quantities. For example, properties, which rely on the knowledge of the density close to the nucleus (hyperfine fields, electric field gradients, etc.), require an all-electron description hemistry 176 (2003) 319–328 321 into (non-overlapping) atomic spheres centered on the ARTICLE IN PRESS tate C atomic sites (I) and an interstitial region (II). For the construction of basis functions—and only for that— the muffin tin approximation (MTA) is used according to which the potential is assumed to be spherically symmetric within the atomic spheres but constant outside. Atomic-like functions are used in region I but plane waves (PW) in region II. Each PW is augmented by a corresponding atomic solution inside every atomic sphere. Three schemes of augmentation (APW, LAPW, APW+lo) have been suggested over the years and illustrate the progress in this development of APW-type calculations that was discussed in a recent paper [10]. Here only a brief summary will be given. The energy dependence of the atomic radial functions ucðr;EÞ can be treated in different ways. In Slater’s APW [11] this was done by choosing a fixed energy E; which leads to a non-linear eigenvalue problem, since the basis functions become energy dependent. In LAPW, suggested by Andersen [12], a linearization of this energy dependence is used by solving the radial Schro¨dinger equation for a fixed linearization energy Ec but adding an energy derivative of this function to increase the variational flexibility. Inside sphere a the atomic function is given by a sum of partial waves (radial functions times spherical harmonics), where L labels the quantum numbers (c;m). X L ½aaKL uacðr0Þ þ baKL ’uacðr0Þ�YLð Kr0Þ: ð3Þ The corresponding two coefficients a and b (weight for function and derivative) can be chosen such as to match each plane wave (characterized by K) to the atomic solution in value and slope at the sphere boundary (for details see e.g., [8,10]). In the APW plus local orbitals (APW+lo) method by Sjo¨stedt et al. [13] the matching is again (as in APW) only done in value. The crystalline wave functions (of Bloch type) are expanded in these APWs leading (in the latter two cases of LAPW or APW+lo) to a general eigenvalue problem. The size of the matrix is mainly given by the number of plane waves (PWs) but is increased slightly by the additional local orbitals that are used. As a rule one needs about 50–100 PWs for every atom in the unit cell in order to achieve good convergence. APW+lo leads—on the one hand—to a significant speedup (by an order of magnitude) and—on the other hand—to a comparable high accuracy with respect to LAPW [14]. The historical development and the details of this latest development, which is the basis for the new WIEN2k program [15], is given in Refs. [10,14]. The new version combines the best features of all APW-based methods. It was known that LAPW converges some- what slower than APW due to the constraint of having differential basis functions and thus it was advantageous to go back to APW. However, the energy-independent K. Schwarz / Journal of Solid S322 basis introduced in LAPW is crucial, since it avoids the non-linear eigenvalue problem of APW, and thus is kept. The local orbitals provide the necessary variational flexibility that make the new scheme efficient [10,13,14]. 3.2. The muffin tin approximation and the full potential The MTA was frequently used in the 1970s and works reasonably well in highly coordinated (closed packed) systems. However, for covalently bonded solids, open or layered structures, MTA is a poor approximation and leads to serious discrepancies with experiment. In all these cases a full-potential treatment is essential. In the full-potential schemes both, the potential and charge density, are expanded into lattice harmonics inside each atomic sphere: X LM VLMðrÞYLMðrÞ ð4Þ and as a Fourier series in the interstitial region: X K VKe i~KK:~rr: ð5Þ Thus, their form (shown for the potential in Eqs. (4) and (5)) is completely general so that such a scheme is termed full-potential calculation. The foundation was laid by the pioneering work of the Freeman group leading to the FLAPW [16,17]. In order to have the smallest number of LM values in the lattice harmonics expansion (Eq. (4)) a local coordinate system for each atomic sphere is defined according to the point group symmetry of the corresponding atom. A rotation matrix relates the local to the global coordinate system of the unit cell. In addition to reducing the number of LM terms in Eq. (4) the local coordinate system also provides orbitals that are properly oriented with respect to the ligands, which may help the interpretation. The choice of sphere radii is not very critical in full- potential calculations in contrast to MTA, in which one would, e.g., obtain different radii as optimum choice depending on whether one looks at the potential (maximum between two adjacent atoms) or the charge density (minimum between two adjacent atoms). There- fore in MTA one must make a compromise between these two criteria which are both reasonable. In full- potential calculations one can efficiently handle this problem and is rather insensitive to the choice of atomic sphere radii. 3.3. Computational considerations In the newest version WIEN2k [15] the alternative basis set (APW+lo) is used inside the ato