Projector augmented-wave method

PHYSICAL REVIEW B VOLUME 50, NUMBER 24 Projector augmented-+rave method 15 DECEMBER 1994-II P. E. Blochl IBM Research Division, Zurich Research Laboratory, CH-8808 Ruschlikon, Switzerland (Received 13 June 1994; revised manuscript received 22 August 1994) An approach for electronic structure calculations is described that generalizes both the pseu- dopotential method and the linear augmented-plane-wave (LAPW) method in a natural way. The method allows high-quality first-principles molecular-dynamics calculations to be performed using the original fictitious Lagrangian approach of Car and Parrinello. Like the LAPW method it can be used to treat first-row and transition-metal elements with affordable effort and provides access to the full wave function. The augmentation procedure is generalized in that partial-wave expansions are not determined by the value and the derivative of the envelope function at some muKn-tin radius, but rather by the overlap with localized projector functions. The pseudopotential approach based on generalized separable pseudopotentials can be regained by a simple approximation. I. INTRODUCTION In the past few decades, electronic structure calcula- tions have made significant contributions to our under- standing of solid-state properties. The majority of such calculations are based on the local-density approximation (LDA) of the density-functional theory. ~ 2 The density- functional theory maps the ground state of an interacting electron gas onto the ground state of noninteracting elec- trons, which experience an effective potential. Numerous methods have been developed to solve the resulting one-particle Schrodinger equation of the LDA. The most widely used electronic structure methods can be divided into two classes: (i) the linear methodss de- veloped by Andersen &om the augmented-plane-wave (APW) method4's and the Korringa-Kohn-Rostocker methods' and (ii) the pseudopotential method based on norm-conserving O,b initio pseudopotentials invented by Hamann, Schluter, and Chiang. A third class, primar- ily employed in chemistry, uses Gaussian basis sets to expand the full wave functions. The linear methods can be subdivided into a vari- ety of methods ranging &om the most accurate linear augmented-plane-wave (LAPW) method to the linear muffin-tin orbital (LMTO) method, which, in a simpli- fied version, even allows some electronic structure calcu- lations to be performed with paper and pencil. The lin- ear methods deal with the full wave functions and treat all elements in the Periodic Table, i.e., 8-, p-, d-, and f-electron systems, on the same footing. The pseudopotential method, when used in combina- tion with a plane-wave basis set, on the other hand, has the advantage of formal simplicity. When applied to ei- ther first-row elements or systems with d or f electrons, even pseudopotentials become very "hard, " so that in practice either very large or complicated basis sets in- stead of plane waves have to be used. Similarly, treating semicore states as valence states, which is often necessary for early transition-metal elements and allmli and alkaline earth metals, results in hard pseudopotentials and sects the transferability of the pseudopotential. Vanderbilt's ultrasoft pseudopotentials9 ~o have improved this situa- tion significantly by relaxing the norm-conservation con- dition that is usually imposed on the pseudopotential ap- proach. This xnethod also allows first-rom and transition- metal elexnents to be dealt with in an economical way. Car and Parrinello have combined the density- functional theory with molecular-dynamics techniques. Here both the electronic structure problexn and the dy- namics of the atoms are solved simultaneously by a set of Newton's equations. In this way not only has the struc- ture deterxnination become a straightforward technique, but the fully dynaxnic time evolution of the atomic struc- ture has also become accessible. The Car-Parrinello method was first applied in the con- text of the plane-wave pseudopotential method. There is considerable interest in applying the same technique to all-electron (AE) methods, which allow one to deal efficiently with first-row and transition-xnetal elements and which supply inforxnation about the wave func- tion close to the nucleus probed by several experimen- tal techniques, but not provided by the pseudopoten- tial approach. These are, among many others, hyper- fine parameters and electric field gradients. ' Sev- eral features of the Car-Parrinello method have been im- plexnented into existing AE methods such as the com- bined minimization of electronic and nuclear degrees of &eedom. To my knowledge, however, no energy- conserving molecular-dynaxnics simulation has been per- forxned to date that can compare in terms of quality with simulations using the pseudopotential approach. This article describes an approach that combines the versatility of the LAPW method with the formal sim- plicity of the traditional plane-wave pseudopotential ap- proach. The method extends the augmented-wave meth- ods, such as the LAPW method, and the pseudopoten- tial method in a natural way. As an AE method it pro- vides the full wave functions that are not directly acces- sible with the pseudopotential approach, and the poten- tial is determined properly &oxn the full charge densi- 0163-1829/94/50(24)/17953(27)/$06. 00 50 17 953 Q~1994 The American Physical Society 17 954 P. E. BLOCHL ties. It will be demonstrated that the accuracy of the method described here compares well with the most ac- curate existing electronic structure methods based on the local-density approximation. The quality of first- principles molecular dynamics obtained with the present AE approach is in line with that of state-of-the-art Car- Parrinello calculations. Hence the first energy-conserving molecular-dynamics calculations based on the full wave functions were made possible. Finally, it can be imple- mented with relatively minor efFort into existing pseu- dopotential codes. The method has many similarities with both the ex- isting linear methods and the pseudopotential approach. We can therefore expect that this method will close the gap between the two. The LAPW method is a special case of the present method, and the pseudopotential for- malisrn is obtained by a well-defined approximation. This article is organized as follows. Section II estab- lishes the principles of the method. Section III describes which approximations are required in real calculations. Section IV derives the expressions for the Hamilton op- erator and forces. Section V describes the implemen- tation in a first-principles molecular-dynamics scheme. Section VI describes the basic ingredients used in the method, such as partial waves and projector functions. Section VII contains a detailed analysis of the errors in- troduced in Section III. Section VIII is devoted to nu- merical test calculations. Section IX shows the relation between the new method and existing approaches. II. FORMALISM A. Projector augmented-wave functions Wave functions of real materials have very different signatures in difFerent regions of space: in the bonding region the wave function is fairly smooth, whereas close to the nucleus the wave function oscillates rapidly owing to the large attractive potential of the nucleus. This is the source of the diKculty of electronic structure methods to describe the bonding region to a high degree of accuracy while accounting for the large variations in the atom cen- ter. The strategy of the augmented-wave methods has been to divide the wave function into parts, namely, a partial-wave expansion within an atom-centered sphere and envelope functions outside the spheres. The enve- lope function is expanded into plane waves or some other convenient basis set. Envelope function and partial-wave expansions are then matched with value and derivative at the sphere radius. Even though the present method has been inspired by the existing augmented-wave methods, I approach the problem in a somewhat difFerent way. The relation of my approach to the commonly used one described above will be described in Sec. IXB. Concerning the following derivation it is emphasized that the present method is, in a certain sense, the most general augmentation scheme. Let us consider the Hilbert space of all wave functions orthogonal to the core states. The physically relevant wave functions in this Hilbert space exhibit strong os- cillations, which make a numerical treatment cumber- some. Therefore, we transform the wave functions of this Hilbert space into a new, so-called pseudo (PS) Hilbert space. Mapping the physical valence wave functions onto the fictitious PS wave functions thus de6ned shall be a linear transformation and it shall transform the physi- cally relevant AE wave functions onto computationally convenient PS wave functions. The PS wave functions will be identified with the envelope functions of the jin- ear methods or the wave functions of the pseudopoten- tial approach. An AE wave function is a full one-electron Kohn-Sham wave function and is not to be confused with a many-electron wave function. All quantities related to the PS representation of the wave functions will hence- forth be indicated by a tilde. This transformation changes the representation of the wave functions in a way reminiscent of the change from a. Schrodinger to a Heisenberg picture. Knowing the trans- formation 7 from the PS wave function to the AE wave functions, we can obtain physical quantities, represented as the expectation value (A) of some operator A, frorii the PS wave functions ~4) either directly as (@~A~4) af- ter transformation to the true AE wave functions ~@) = 7 ~4) or as the expectation value (A) = (4~A~4) of a PS operator A = 7 tA7 in the Hilbert space of the PS wave functions. Similarly we can evaluate the total energy directly as a functional of the PS wave functions. The ground-state PS wave functions can be obtained from ~&71+)j, ~t~(@) Next, we choose a particular transformation. Since we will exploit the characteristics of particular atom types, we consider only transformations that differ &om identity by a sum of local, atom-centered contributions 7R such that 7 =&+).4. R Each local contribution 7~ acts only within some aug- mentation region OR enclosing the atom. This implies that AE and PS wave functions coincide outside the aug- mentation regions. The equivalent of the augmentation region in the linear methods is the mufBn-tin or atomic sphere. In the pseudopotential method the augmentation region corresponds to the so-called core region. The local terms 7R are defined for each augmenta- tion region individually by specifying the target func- tions ~P;) of the transformation 7 for set of initial func- tions ~P;) that is orthogonal to the core states and other- wise complete in the augmentation region, namely, by ~P;) = (1 + 7~)~P, ) within O~. I call the initial states ~P;) PS partial waves and the corresponding target func- tions ~P;) AE partial waves. A natural choice for these functions for the AE partial waves are solutions of the radial Schrodinger equation for the isolated atom, which are orthogonalized to the core states if necessary. Hence the index i refers to the atomic site R, the angular mo- mentum quantum numbers I = (I., m), and an additional 50 PROJECTOR AUGMENTED-WAVE METHOD 17 955 index n to label diferent partial waves for the same site and angular mornentuxn. For each such AE partial wave let us choose a PS partial wave denoted by 1$;). The PS partial waves must be identical to the corresponding AE partial waves outside the augmentation region and should themselves form a complete set of functions within the augmentation region. The rexnaining degree of &eedom in the choice of the PS partial waves will be exploited to map the physically relevant AE wave functions onto computationally convenient PS wave functions. In our case these are smooth functions. This forrnal definition must be turned into a closed expression for the transformation operator. We make use of the fact that, within the augmentation region, every PS wave function can be expanded into PS partial waves: I@) = ) 1$;)c; within OR. Since 1$,) = 7 1$;), the corresponding AE wave function is of the form 14) = 714') = ) Ip;)c; within OR, (4) with identical coefficients c, in both expansions. Hence we can express the AE wave function as where the expansion coefficients for the partial wave ex- pansions rexnain to be determined. Since we require the transformation 7 to be linear, the coefficients must be linear functionals of the PS wave functions. Hence the coefficients are scalar products (6) of the PS wave function with some fixed functions (p;I, which I will call projector functions. There is exactly one projector function for each PS partial wave. The projector functions must fulfill the condition P,. 1$,)(p;I = 1 within OR, so that the one-center ex- pansion P,. 1$;)(plilI) of a PS wave function is identical to the PS wave function 14') itself. This implies that The projector functions are localized in the augmentation region, even though more extended projector functions could in principle also be chosen. The most general form for the projector functions is (p;I = g.(((fi, lg~))),. (f&I, where the I fz) form an arbitrary, linearly independent set of functions. The projector functions are localized if the functions I f;) are localized. The reader interested at this point in a practical procedure to determine partial waves and projector functions might wish to jump to Sec. VI. In summary, a linear transformation between the valence wave functions and fictitious PS wave functions has been established. Using this trans- formation, the AE wave function can be obtained from the PS wave function by The three quantities that determine this transformation are (i) the AE partial waves 1$;) obtained by radially integrating the Schrodinger equation of the atomic en- ergy for a set of energies e~ and orthogonalization to the core states; (ii) one PS partial wave 1$;), which coincides with the corresponding AE partial wave outside some augmentation region for each AE partial wave; and (iii) one projector function Ip;) for each PS partial wave lo- calized within the augmentation region and which obeys the relation (p;1$~) = b;z. The partial waves are functions on a radial grid, mul- tiplied with spherical harxnonics. In our case the PS wave functions are expanded into plane waves, but other choices are equally possible. The projectors are also cal- culated as a radial function times spherical harxnonics, but are then transformed into the same representation as the PS wave functions, which, in our case, is a plane- wave representation. Since the projectors are tied to the atomic positions and since their shape is independent of the potential, their Fourier components are expressed as a product of a form factor and a structure factor. The core states 14") are decomposed in a way simi- lar to the valence wave functions. They are decomposed into three contributions: a PS core wave function I@'), which is identical to the true core state outside the aug- mentation region and a smooth continuation inside; an "AE core partial wave" IgP), which is identical to the AE core state 14") and is expressed as a radial function times spherical harmonics; and finally a "PS core partial wave" 1$'), which is identical to the PS core state 14'), but rep- resented as a radial function times spherical harmonics. The core state is therefore expressed as (IO) In contrast to the valence states, no projector functions need be defined for the core states, and the "coefficients" of the one-center terms are always unity. Furthermore, consistent with the &ozen-core approximation, the core states are imported &om an isolated atom. In prac- tice, a soft core scheme with core states that adjust to the instantaneous potential is also conceivable (see Sec. VII D 2), but has not been implemented. In the fol- lowing, the core states are implicitly included when sum- ming over energy states. Note that the corresponding coefficients are not defined via the scalar product with a PS wave function, even though, for the sake of simplicity, I will still use the symbol for all states. It should be noted that the frozen-core approximation allows certain nontrivial changes of the core wave func- tion during the self-consistency or molecular-dynamics simulation. The frozen-core approximation only restricts the variational degree of freedom to a simple unitary transformation among the core states (and occupied va- 17 956 P. E. BLOCHL 50 lence states). It does allow mixing among the core states due to changing potential. Therefore, to test the accu- racy of the frozen-core approximation one should never compare the core states of the isolated atom on a one-to- one basis with those obtained &om a relaxed-core calcu- lation in a crystal or molecule. At this point I will not discuss the components of the projector augmented-wave (PAW) method further. They are described in Sec. VI of this article. I will, however, continue to impose the condition that the AE and PS partial waves form complete sets of functions within the augmentation regions. In practical calculations the num- ber of partial waves and projectors needs to be truncated. The way to truncate the series and the errors involved are also described in detail in later sections of the paper. Here and in the following I will make extensive use of Dirac's bra and ket notation. A wave function in real space is written as (rl@) = 4(r); its complex conjugate function is (@lr& = 4'*(r). The Fourier components of the wave function are (Gl@) = 4(G) with a similar def- inition of its complex conjugate. A plane wave is of the form (rlG) = exp(iGr). I have adopted the convention for the Fourier transform that the forward transform of a function f is (7'If) = g&(rlG)(GI f) and the backward transform is of the form (Glf) = 1/V fv dr(GIr&&rI f&, where V is the volume of the unit cell. B. Operators Since in the PAW method the PS wave functions in- stead of the AE functions play the role of the variational parameters, we need to be able to obtain observable quantities as the expectation values of the PS wave func- tions. As the representation of the wave functions has been changed, we also need to transform our operators into new, so-called PS operators. Consider soine operator A: Its expectation value (A) = P„ f„(@„IAIDO'„), where n is the band index and f„ is the occupation of the state, can be obtained alternatively as (A) = g„f„&@„IAIDO'„). For quasilocal operators, such as the kinetic-energy operator —V'2/2 and the real-space projection operator Ir)(rl, which are needed to evaluate total energy and charge density, the PS operator has the form The general form of an operator is strongly reminiscent of generalized separable pseudopotentials. The PS op- erator contains three parts: The first part is an operator that directly acts on the PS wave function and is evalu- ated either in real or reciprocal space. The remaining two parts contain the projectors and the expectation value of the operator either between the AE or the PS partial waves, which can be easily obtained on radial grids using spherical harmonics and Clebsch-Gordan coefBcients. If the partial waves are unbound, the individual terms A' and A are not defined. However, since the PS and AE partial waves are identical outside the augmentation re- gion, these tails cancel exactly for each pair of partial waves. In practice, this problem is solved by truncat- ing the AE and PS partial waves somewhere outside the augmentation region in a completely identical way. There is an additional freedom to add a term of the form to the right-hand side of Eq. (11), where B is an ar- bitrary operator that is localized within t