Extension of the basis set of linearized augmented plane wave method (LAPW) by using supplemented tight binding basis functions
In order to increase the accuracy of the linearized augmented plane wave method (LAPW) we present a new approach where the plane wave basis function is augmented by two different atomic radial components constructed at two different linearization energies corresponding to two different electron bands (or energy windows). We demonstrate that this case can be reduced to the standard treatment within the LAPW paradigm where the usual basis set is enriched by the basis functions of the tight binding type, which go to zero with zero derivative at the sphere boundary. We show that the task is closely related with the problem of extended core states which is currently solved by applying the LAPW method with local orbitals (LAPW+LO). In comparison with LAPW+LO, the number of supplemented basis functions in our approach is doubled, which opens up a new channel for the extension of the LAPW and LAPW+LO basis sets. The appearance of new supplemented basis functions absent in the LAPW+LO treatment is closely related with the existence of the component in the canonical LAPW method. We discuss properties of additional tight binding basis functions and apply the extended basis set for computation of electron energy bands of lanthanum (face and body centered structures) and hexagonal close packed lattice of cadmium. We demonstrate that the new treatment gives lower total energies in comparison with both canonical LAPW and LAPW+LO, with the energy difference more pronounced for intermediate and poor LAPW basis sets.
pacs:71.15.-m; 71.15.Ap; 71.20.-b
The choice of a basis set which may first appear as “the black art” Sza () is being constantly debated within the quantum chemistry community. It is well known that in molecular calculations there are two main groups of molecular basis sets introduced by Pople and collaborators (see Sza () and references therein) and more recently by Dunning Dun (). Both groups supply us with the whole hierarchy of basis sets, where at each step we can enrich the main set by polarization functions of high orbital type or by diffuse functions. For example, one has to add polarization functions if polarization effects are expected to be important, or diffuse functions if we want to refine the description of extended molecular states. Not surprisingly, the actual choice of a basis set depends on the task to be solved and is considered as a difficult problem. For heavy and laborous calculations the choice of a basis set is crucial since on one hand we want to obtain a reliable result and on the other hand minimize the computer time to achieve the goal.
In contrast to this complicated hierarchy of molecular basis sets, the choice of bases in electronic band structure calculations and here we imply mainly the linear augmented plane wave (LAPW) method And (); Koe (); blapw (), seems rather simple. The number of augmented plane waves is commonly determined by the parameter , where is the smallest muffin-tin (MT) radius, and is the maximal value of the plane wave vector. The value implies that the kinetic energy cut off is (in atomic units). However, in practice some band structure calculations can not be carried out without so called local orbitals (LO) wien2k (). Such situations occur in systems with semicore electron states which cannot be fully confined within the MT-spheres. The problem of extended core states and its relation to our approach is discussed in detail in the next section. What concerns us here is that the introduction of local orbitals represents an extension of the canonical LAPW basis set albeit the form of new basis states (local orbitals) is very different from the standard augmented plane wave basis function. The LAPW+LO scheme proposed by Singh Sin-2 (), wien2k () is practical, but the way it has been introduced is not fully satisfactory. The form of the local orbital functions is not derived from a general approach, and arguments for adding LO basis states are purely variational.
In the present paper we show that the appearance of supplemented basis functions can be understood as a result of refinement of the LAPW band scheme when in an effort to increase its accuracy we use two linearization energies (corresponding to two electron bands). We will demonstrate that new basis states are of two types. The first group gives the local orbitals in the form suggested by Sing Sin-2 (). however, the basis functions of second type have different form which is not used in the LAPW+LO method. Therefore, the canonical LAPW basis set and also the LAPW+LO basis set can be extended to a more complete basis set. New basis functions and consequences of their introduction are closely examined in the present work.
The paper is organized as follows. We start with revisiting the problem of extended core states which gives rise to the LAPW+LO scheme and formulate our initial statement for the refinement of the LAPW method, Sec. II. In Sec. III we present our method which results in adding tight binding basis functions to the canonical LAPW basis set. In Sec. IV we apply the method to electron band structure calculations of the face centered and body centered lattice of lanthanum and the hexagonal close packed lattice of cadmium. Our conclusions are summarized in Sec. V.
Ii The problem of extended core states
The linear augmented plane wave (LAPW) method And (); Koe (); blapw () is probably the most precise method for electronic band structure calculations and is widely used for the calculation of materials properties wien2k ().
In the LAPW method And (); Koe (); blapw (); wien2k () space is partitioned in the region inside the nonoverlapping muffin-tin () spheres and the interstitial region . The basis functions where are given by
with radial parts
Here the index refers to the type of atom (or -sphere) in the unit cell, the radius is counted from the center of the sphere (i.e. ), is the volume of the unit cell. Radial functions are solutions of the Schrödinger equation in the spherically averaged crystal potential computed at the linearization energy , and is the derivative of with respect to at . The coefficients and are found from the condition that the basis function is continuous with continuous derivative at the sphere boundary, ( is the radius of the -sphere ). In the following for compactness we omit the index and restore it when needed. Linearization energies are chosen close to average values of corresponding band energies or to the Fermi level. The extended electron basis states defined by Eq. (1) as a rule are orthogonal to the core states. This is a consequence of the relation
applied to a core state with orbital quantum numbers and the radial wave function , and a partial radial function of valence state, , Eq. (2), with the same angular dependence. Notice that Eq. (3) ensures the orthogonality between the extended and core states if two conditions at the sphere boundary are satisfied for each of the core states,
Semicore states leaking out of the -regions should be treated as extended states. This in turn requires that the linearization energy is chosen near the energy of the semicore level, , because the LAPW basis describes only states near well. However, as is quite far from the Fermi energy and the valence band energy, the choice inevitably gives poor description for partial valence states. On the other hand, the option is not satisfactory for the semicore states situated substantially deeper in energy. As discussed in Refs. Sin-1 (); Sin-2 (); Bla-1 (); Sin-3 (); Sin-4 (); Goe () there is no simple solution to this dilemma. Even worse, in many cases the attempt to use a single value of for both valence and semicore states leads to the appearance of so called “ghost bands” Sin-1 (); Sin-2 () giving false band energy positions. As a remedy one can divide the energy spectrum in two windows (energy panels) and use two different sets of for calculations of semicore and valence states, respectively blapw (). This technique however is also not fully satisfactory because now there is no single Hamiltonian matrix for the problem and strict orthogonality between electron states belonging to different energy windows is not guaranteed. Ideally, in the MT-region there should be two different types of radial components with the same angular dependence . That is, in Eq. (1)
|refers to the semicore states with , and|
|refers to the valence states with . Both states, Eqs. (5b) and (5c), should merge to a single wave component of the plane wave|
at the surface of MT sphere. Now, however the boundary problem becomes ill-defined, because for the component there are four coefficients, , , and for only two boundary conditions.
In Ref. Sin-2 () Singh has proposed to increase the number of boundary conditions to four by matching the value of the basis function and its first three radial derivatives at the sphere surface. This gives rise to super-linearized APW method denoted as SLAPW-4 Sin-2 () because four functions, four coefficients and four boundary conditions are involved. In a simpler super-linearized modification called SLAPW-3 Sin-2 () the first radial part , Eq. (5b), is supplemented by only one function (instead of , Eq. (5c)). The three coefficients (, , ) are determined by requiring continuity of the basis function and its two derivatives.
In comparison with standard LAPW method in both SLAPW modifications there are additional requirements for the plane wave convergence. Indeed, the plane wave expansion of the interstitial region must converge either to the correct second and third derivative (SLAPW-4) or to the second derivative (SLAPW-3). Because of that more plane waves are needed to satisfy these additional requirements, the plane wave energy cutoff parameter should be increased and calculations become much more costly Sin-2 (),blapw ().
To circumvent the problem and improve the LAPW efficiency Singh put forward a third approach based on local orbitals (LAPW+LO) Sin-2 (). In the LAPW+LO approach the same three radial functions as in SLAPW-3 are used (i.e. , and ), but the coefficient of is fixed (say, ) and the two remaining coefficients (, ) are found from the conditions that the local orbital goes to zero with zero derivative at the sphere boundary. Nowadays, LAPW+LO is widely used for band structure calculations of solids with semicore states blapw (); wien2k (). However, conceptually the LAPW+LO method is understood as a procedure giving additional variational freedom through an increase of the number of basis functions. It is not clear why additional basis functions should include these particular components (i.e. , and ). The proposed zero boundary conditions for local functions are not derived from a general physical statement.
Inspired by the LAPW+LO method Sin-2 () in the present study we formulate a more general approach to the problem. Unlike the LAPW+LO approach which uses variational arguments for its foundation, we will derive supplemented basis states from the initial requirement that two different radial functions (i.e. and , Eqs. (5b), (5c)), having the same angular part merge in a single plane wave function , Eq. (5d), in the interstitial region. Unlike SLAPW-4 or SLAPW-3 we retain only two joining conditions across the -sphere boundary. As a result, we will obtain two types of supplementary tight-binding basis functions (see Eqs. (15b) and (15c) below), satisfying Bloch’s theorem.
Iii Description of the method
As discussed in Sec. II, in the case of semicore states we have two types of radial solutions in the MT-region with the same angular dependence but different linearization energies and : , Eq. (5b), and , Eq. (5c). One of the radial functions can refer to extended states, i.e. , while the other can refer to supplementary angular states . As we will see later in Sec. IV in practice we describe the semicore states as extended states with and valence states with the same as supplementary states for which . (For metals one can take .) Since in the interstitial -region both types of solutions are represented by the plane wave function , Eq. (5d), they become indistinguishable there. In the LAPW method there are two matching conditions (for the function and its derivative) on the sphere boundary. Therefore, in our case we have
Here we adopt short notations , , , , , , and have used the Rayleigh expansion of the plane wave on the sphere surface. Since there are four coefficients (, , and ) and only two equations, it is clear that the general solution to Eqs. (6a), (6b), forms a two dimensional linear space with two linear independent basis vectors.
Further, introducing standard LAPW quantities , , where
Notice that the standard LAPW solution for and without supplementary states, i.e. when , , can be found from the following system
Defining auxiliary quantities and
where and are arbitrary numbers. The full radial component of the basis function inside the -sphere , Eq. (5a), is written as
(Here notations , etc. refer to radial functions.)
Notice that since the coefficients and are arbitrary, they should be found from the standard variational procedure by requiring the minimization of the LAPW ground state energy. Furthermore, the form (14) suggests considering three linear independent radial parts (i.e. , , ) instead of the single function . Explicitly,
Here , etc. are corresponding radial functions and
The first function, Eq. (15a), is in fact the standard radial part of the type, , Eq. (2), entering the usual LAPW basis function , Eq. (1). Its coefficients and are given by the LAPW boundary relations, Eqs. (9a), (9b), Two other functions however are very different from and should be included to the LAPW basis set as extra basis states,
where . The important thing is that their coefficients , are found from Eq. (12a), while coefficients , from Eq. (12b). In respect to two new functions , Eqs. (12a), (12b) impose the following boundary conditions
These relations have a simple interpretation: new supplementary basis functions are required to be orthogonal to the standard LAPW radial functions. We want to stress here, that the conditions (18a) and (18b) are not assumed or introduced at our will. They are derived from the initial equations (6a), (6b) [or equivalently from Eqs. (8a), (8b)] and are used to obtain the general solution, Eqs. (13a)-(13d).
From Eq. (15b), (15c) and (16) it follows that the supplementary basis states in principle depend on the index , i.e. . (We recall that , where is a vector of the reciprocal lattice.) However, since all functions with different index have the same radial part, for , or for , they are simply proportional to each other, . Therefore, to avoid the linear dependence we should choose only one set of functions corresponding to a single index . The obvious choice is to use the function with , and the reciprocal lattice vector . In that case the coefficient [compare with Eq. (16)] can be further rationalized by omitting the multiplier [or equivalently, including it in factors and , Eq. (10a), (10b)]. Thus, we substitute with
(In principle, since the local function is not orthonormal, we can simply put , but the form (19) being similar to the constant coefficient for the standard LAPW basis function, simplifies some expressions for programming.)
To study the transformational properties of supplementary basis functions , Eq. (17), we first rewrite them in the following form,
where for each site we have introduced two local functions () of the -type,
Notice that each local function is strictly confined inside the -sphere , because both and satisfy the boundary conditions (18a), (18b). We can then extend the function to the interstitial region () by requiring . For the whole crystal we thus have
This is a clear manifestation of the tight binding wave function. The multiplier in Eqs. (22) and (20) ensures that the supplementary wave functions obey the Bloch theorem. It is worth noting that usually the tight-binding description is spoiled by the presence of overlap between tails of wave function centered at neighboring sites. In the present method the tight-binding functions, Eqs. (22) and (20), are free from this disadvantage because the local functions (and their first derivatives) go to zero at the sphere boundary and the overlap is absent. Thus, supplementary tight-binding functions can be considered as additional basis states orthogonal to the standard LAPW basis set.
All matrix elements between the supplementary basis functions in the spherically symmetric potential are quoted explicitly in Appendix B, and all matrix elements between and standard LAPW basis functions are listed in Appendix C. For briefness we do not quote here the partial charges and electron density associated with supplementary basis states. They are tightly connected with the overlap matrix elements given by Eqs. (28a), (30a) and (31a) of Appendix B, and Eqs. (36a), (37a) of Appendix C. Concerning the full potential expressions for the extended basis set it is worth noting that after some algebraic transformations the equations can be obtained by selecting in standard FLAPW equations the contributions with the orbital indices referring to the components of supplemented states and combining them together according to Eqs. (15b), (15c), (17).
The tight binding basis functions have a very important and practical property: they work even in the case when their expansion energy lies not far from the LAPW linear expansion energy . (We recall that the whole procedure is designed to treat the complicated case of semicore states when is supposed to be separated from by at least 10 eV.) The limiting case is considered in detail in Appendix D, and also discussed in calculations of Cd in Sec. IV.3.
Iv Practical Implementation
iv.1 Face centered cubic structure of La
We have applied the method developed in Sec. III to full potential electron band structure calculations of face centered cubic (fcc) structure of lanthanum. Atomic lanthanum has completely filled 5 semicore electron shell lying at -22.12 eV which can slightly mix with valence states (, , ) at energies from -3 to -2 eV. For lanthanum here and below we use the Perdew-Burke-Ernzerhof (PBE) PBE () variant of the generalized gradient approximation (GGA) which gives rather accurate lattice constants for our systems.
We have employed our original version of LAPW code with the potential of general form Nik (). Integration in the irreducible part of the Brillouin zone has been performed over 240 special points. Angular expansions for the electron density and wave function inside MT-sphere have been done up to . The number of basis functions has been limited by the condition , resulting in 65 basis states. In addition to the standard LAPW basis functions we have considered 6 supplementary tight binding basis functions with the angular dependence, which are located strictly inside the MT-sphere, Sec. III.
The semicore states have been treated as band states with the LAPW linear expansion energy lying 0.5 eV above the band bottom energy (which is eV for the equilibrium lattice constant Å). The linear expansion energy for supplementary tight binding states has been fixed at 1 eV below the Fermi energy, eV. Under these conditions two supplementary radial functions shown in Fig. 1 are given by
|where , , and|
where , .
Notice that the number of nodes for both radial functions is three (excluding points with and ) which allows us to consider these functions as “compressed” basis states (i.e. with the principal quantum number ) strictly confined within MT-sphere. The iteration procedure with supplemented tight binding functions has been stable converging to a self consistent solution without additional difficulties. Our PBE-GGA calculations result in equilibrium fcc lattice constant Å which compares well with the experimental value, Å, Ref. fcc-La, .
To compare our treatment with the standard (LAPW+LO) method which uses only the first local function , Eq. (23a), we have performed a series of calculations, the results of which are summarized in Tables 1, 2, 3. In all cases the present method gives lower values of the total energy, Table 1. However, the energy difference which amounts to 0.145 eV for the poor basis set () becomes smaller for the intermediate basis sets and decreases to a small value (0.004 eV) for the best basis set (). Nevertheless, even this difference is clearly noticeable in energy band characteristics. In particular, band energy spectrum demonstrates that the energy difference is of 0.07 eV for the best basis set, Tables 2, 3. It is worth noting that for all basis sets the present method results gives smaller band energies, Table 2.
iv.2 Body centered cubic structure of La
For the body centered cubic (bcc) phase of lanthanum we have used a plane wave cut off parameter (79 basis states plus 6 additional tight binding states), for the expansion of the electron density and wave functions inside MT-spheres, 285 special points in the irreducible part of the Brillouin zone during the self-consistent procedure.
As for fcc-La the LAPW linear expansion energy for the extended states of bcc-La has been has been chosen at 0.5 eV above the band bottom energy (i.e. eV for the equilibrium lattice constant Å). The linear expansion energy for the supplemented tight binding states has been fixed at 1.0 eV below the Fermi energy ( eV). Two supplemented radial functions, quoted in Eqs. (23a) and (23b), are defined by the coefficients , for , and , for .
The equilibrium lattice constant found for bcc-La, Å, is in good correspondence with the experimental value Å bcc-La (). The calculated band structure of bcc lattice of lanthanum (PBE exchange and correlation) is plotted in Fig. 2.
As for fcc-La, the present scheme gives lower values of total energy of bcc-La for all basis sets. The total energy difference increases with worsening of the basis quality, Table 4. For the best basis set () the present method gives energy spectrum shifted downwards by eV in comparison with FLAPW+LO band energy values, Tables 6, 5. The difference between energy bands parameters reaches 0.2 eV for poor basis set (), Table 5. For the best basis set the total energy of bcc-La (phase) is eV higher than the total energy of fcc-La (phase), Tables 1 and 4, which is in agreement with the fact that at normal pressure La exists only at high temperatures (1138 K) in narrow temperature range (53 K) bcc-La ().
iv.3 Hexagonal close packed structure of Cd
Hexagonal close packed (hcp) cadmium is a special case because its states are not separated from the valence states by an energy gap, Fig. 3.
In fact there is a small overlap between the top of band and the bottom of band. This implies that if the linear expansion energy of extended states is properly chosen (for example, at 0.5 eV above the band bottom energy, eV) the electronic band structure of cadmium can be carried out without supplemented tight-binding basis states of type. Thus, both calculations i.e. with and without supplemented states, can be directly compared with each other. The second peculiarity is that because of the gapless energy spectrum the radial distribution of states does not change much throughout the valence band. If we chose the linear expansion energy of supplemented states at 1.0 eV below the Fermi energy ( eV), then the energy difference eV is a relatively small value. In that case one can expect that the extended states and the supplemented states are close to being linearly dependent, and we can test the scenario described in Appendix D. (Two supplementary radial functions, Eqs. (23a) and (23b), are specified by the coefficients , for and , for .)
The technical parameters of FLAPW calculations were the following: the plane wave cut off parameter (149 basis states and 10 supplemented tight binding states), for the expansion of electron density and wave functions inside the MT-spheres, 216 special points in the irreducible part of the Brillouin zone during the self-consistent procedure, and the PBE PBE () form of the exchange correlation potential.
We observe that for all basis sets the present approach gives lower total energy values, Table 7. For the best basis set () the obtained total energy difference with the LAPW+LO value, 0.02 eV, is larger than for fcc-La or bcc-La. In the present approach band energies computed with the best basis set are lowered by 0.04-0.09 eV in comparison with LAPW+LO and LAPW values, Tables 8 and 9. These energy shifts are also typical for the other points of the Brillouin zone.
We have presented a new method for the improvement of the LAPW description of the electronic band structure by using two linearization energies for the same partial component. Starting with two LAPW radial functions, Eqs. (5b) and (5c), having the same angular dependence but different linearization energies ( and ) inside MT-spheres, we have demonstrated that their augmentation to the basis plane wave can be performed by constructing additional basis functions () in the form of Eq. (17) [two functions, Eqs. (15b) and (15c), for each -component]. The supplementary basis functions have zero values and slopes on the sphere surface, Eqs. (18a), (18b), and are linear independent of the usual LAPW basis states. The constructed basis functions are of the tight-binding type, Eq. (22), and obey Bloch’s law.
In contrast to the LAPW+LO method with only one supplemented function, Eq. (15b), in our treatment for each -component there are two supplemented functions [Eq. (15b) and (15c)]. The second basis function (absent in LAPW+LO) closely examined in this work, owes its appearance to the function in the canonical LAPW method. Thus, the basis sets of LAPW and LAPW+LO methods can be extended further by adding supplemented functions of the tight-binding type, Eq. (15c).
In Sec. IV, the present method with extended basis set has been applied to the study of the face centered and body centered phases of lanthanum (La and La) with the semicore shell separated by a gap of forbidden states from the valence states and to the hexagonal close packed structure of cadmium, where the semicore states overlap with the valence states. In all cases we have observed a systematic improvement in the values of total energy in comparison with the standard LAPW+LO treatment, Tables 1, 4, 7. The difference with LAPW+LO total energy is only 0.003-0.004 eV for La and 0.019 eV for Cd for the best basis set () but significantly increases in going to intermediate ( or 8) and poor ( or 7) basis sets.
Acknowledgements.A.V.N. acknowledges useful discussions with B. Verberck, E.V. Tkalya, A.V. Bibikov.
The matrix elements for the overlap and Hamiltonian operator between supplementary states,
are partitioned in three different blocks, when (block ), (block ), and , or , (block ).
For the first block () we have