A fusion of the LAPW and the LMTO methods: the augmented plane wave plus muffin-tin orbital (PMT) method
We present a new full-potential method to solve the one-body problem, for example, in the local density approximation. The method uses the augmented plane waves (APWs) and the generalized muffin-tin orbitals (MTOs) together as basis sets to represent the eigenfunctions. Since the MTOs can efficiently describe localized orbitals, e.g, transition metal 3 orbitals, the total energy convergence with basis size is drastically improved in comparison with the linearized APW method. Required parameters to specify MTOs are given by atomic calculations in advance. Thus the robustness, reliability, easy-of-use, and efficiency at this method can be superior to the linearized APW and MTO methods. We show how it works in typical examples, Cu, Fe, Li, SrTiO, and GaAs.
pacs:71.15.Ap, 71.15.Fv 71.15.-m
There are several kinds of all-electron full potential (FP) methods to obtain numerically-accurate solutions in the local density approximation to the density functional theory Kohn and Sham (1965). Among such FP methods, most popular ones are the linearized augmented plane wave (LAPW) method, and the projector augmented-wave (PAW) method Andersen (1975); Soler and Williams (1989); Blochl (1994); Kresse and Joubert (1999). They both use plane waves (PWs) to expand the eigenfunctions in the interstitial regions. However, PWs do not efficiently describe the localized character of eigenfunctions just around the muffin-tins (MTs). For example, oxygen (denoted as O() below) and transition metal’s 3 orbitals are well localized and atomic-like even in solids, thus we need to superpose many PWs to represent these orbitals. For example, as shown in Ref.Liu et al., 1994 (and below), the energy cutoff for the augmented PW(APW) 15Ry is required in fcc Cu to obtain 1 mRy convergence for total energy in LAPW. In contrast, such orbitals can be rather effectively represented by localized basis in real space. In fact, it is already accomplished in the linearized muffin-tin orbital (LMTO) method, which differs from the LAPW method in that envelope functions consists of atomic-like localized orbitals Mesfessel and van Schilfgaarde (); Methfessel et al. (2000) instead of PWs. Such a localized augmented waves are called as MT orbital (MTO).
To circumvent the inefficiency in the LAPW method, we have implemented a new method named as linearized augmented plane wave plus muffin-tin orbital (PMT) method. The PMT method includes not only the APWs but also MTOs in its basis set. Our implementation becomes LAPW in the no MTO limit. As we show later, we can very effectively reduce the number of basis set by including MTOs; we see the rapid convergence of the total energy as a function of the number of APWs (or energy cutoff ). As far as we tested, APWs with Ry in addition to minimum MTOs will be reasonably good enough for usual applications; e.g. for 1 mRy convergence of total energy for Cu. Even in comparison with the LMTO method of Ref.Methfessel et al. (2000), the PMT method is quite advantageous in its simplicity. The parameters to specify these minimum MTOs ( and for each channel. See next paragraph.) are automatically determined by atomic calculations in advance. This is a great advantage practically because optimization of these parameters is a highly non-linear problem Methfessel et al. (2000) which makes the LMTO difficult to use. Thus the PMT method can satisfy the requirements for latest first-principle methods, reliability, speed, easy-of-use, and robustness very well. In what follows, we explain points in our method, and then we show how it works.
We adapted a variant of the LMTO method developed by Methfessel, van Schilfgaarde, and their collaborators Mesfessel and van Schilfgaarde (); Methfessel et al. (2000). This method uses the “smooth Hankel function” invented by Methfessel as a modification of the usual Hankel functions so as to mimic the atomic orbitals Methfessel et al. (2000); Bott et al. (1998). It contains an additional parameter called as the smoothing radius in addition to the usual energy parameter which specifies the damping behavior. By choosing , we can control bending of the function just outside of MT. By using this degree of freedom, we can reproduce eigenvalues of atoms very well even if we substitute the eigenfunction outside of MT with such a smooth Hankel function. This is an important feature of the function in comparison with others like Gaussians,which are not directly fit to the atomic orbitals. Analytic properties of the smooth Hankel are analyzed in detail by Bott et al, and all the required operations, e.g, Bloch sum, to perform electronic structure calculations are well established Bott et al. (1998). The augmentation procedure requires the one-center expansion of the envelops functions. In our method, the one-center expansion is given as the linear combinations of the generalized Gaussian (Gaussians polynomials), where the expansion coefficients are determined by a projection procedure as described in Sec.XII in Ref. Bott et al., 1998. Then the generalized Gaussians in each angular momentum -channel are replaced by the linear combinations of a radial function and its energy derivative in the usual manner of augmentation Methfessel et al. (2000). When we use high energy cutoff Ry, we needed to use generalized Gaussians for the one-center expansion; however, is good enough for practical usage with Ry.
Another key point in our method is the smoothed charge density treatment introduced by Soler and Williams Soler and Williams (1989) to the LAPW method. The charge density is divided into the smooth part, the true density on MTs, and the counterpart on MTs. The smooth part is tabulated on regular uniform mesh in real space, and the others are tabulated on radial mesh in the spherical harmonics expansion in each MT sphere. The PAW methods Blochl (1994); Kresse and Joubert (1999) also use this treatment. It enables low-angular momentum cutoff for augmentation and makes the calculation very efficient. As for the regular mesh in real space, we usually need to use the spatial resolution corresponding to the cutoff energy Ry, which is determined so as to reproduce the smooth Hankel function well in real space. Note that we still have some inefficiency, e.g. an atom in a large supercell requires a fine mesh everywhere only in order to describe the density around the atom. This problem is common to any method which uses an uniform mesh for density.
Though the LMTO formalism shown in Mesfessel and van Schilfgaarde (); Methfessel et al. (2000) is intended for such MTOs constructed from the smooth Hankel functions, it is essentially applicable to any kind of envelope functions. As we explained above, our formalism is not so different from LAPW/PAW formalisms shown in Kresse and Joubert (1999); Blochl (1994) except in the augmentation(projection) procedure. The atomic forces are calculated Methfessel et al. (2000); Methfessel and van Schilfgaarde (1993) in the same manner as in PAW Kresse and Joubert (1999). For deep cores, we usually use a frozen core approximation which allows the extension of core densities outside of MT (but no augmentation) Mesfessel and van Schilfgaarde (). Further, we use the local orbital (lo) method E. Sjostedt (2000) in some cases; for example, to treat 3 semicore for Ga (denoted as Ga([lo]) below), or to reproduce high-lying bands for Cu by Cu([lo]).
Results: In Fig.1, we plot the total energies as functions of the inverse of the number of basis at the point in order to observe its convergence as (the number of basis is controlled by as shown on Fig.1 together). Here we includes minimum MTOs whose parameters and are fixed by the atomic calculations (here “minimum” means only the atomic bound states). The lo’s are included as explained in Fig.1.
There is a problem of linear dependency in the basis set when we use large . For example, in Li, we could not include MTOs of Li() (Li() and Li()) as basis at , because then they are well expanded by PWs. This occurs also for other cases; when we use high enough to expand a MTO, the rank of the overlap matrix of basis set is reduced by one. Thus possible ways to use large are: (1) keep only well localized MTOs which are not yet expanded by given , or (2) remove a subspace of basis through the diagonalization procedure of the overlap matrix before solving the secular equation. For (2), we need to introduce some threshold to judge the linear dependency. This can cause an artificial discontinuity when changing lattice constants and so on. Thus (1) should be safer, but here we use the procedure (2) with careful check so that such discontinuity do not occur. We use the number of basis after the procedure (2) for to plot Fig. 1. Even for LAPW cases, we applied the procedure (2), e.g, to the case for 20Ry in SrTiO; then we reduce the dimension of Hamiltonian from 606 to (data point at the left end of SrTiO in Fig.1).
As is seen in all the cases, the PMT method shows smooth and rapid convergence for the total energy at . On the other hand, the convergence in LAPW (no MTO limit in our PMT implementation) shows a characteristic feature; it is way off until it reaches to some , e.g, Ry in Cu. This is consistent with previous calculations Liu et al. (1994). is required to reproduce the localized orbitals. We also show the convergence behaviors for band gap and magnetic moments by dotted lines (right y-axis); they are quite satisfactory.
Fig.2 shows the energy bands up to 100eV above the Fermi energy. Though the Cu panel in Fig.1 shows the little effects of local orbital to the total energy, it affects energy bands above eV. Calculations including [lo] gives good agreements each other. This means that we have no artificial bands (ghost bands). The MTO(high) panel is by the pure LMTO method where we use 34+5 basis (+[lo]). With some careful choice of and , the LMTO method can be very efficient and accurate.
Fig.3 shows the total energy as function of the lattice constant. All lines looks very parallel. This shows that PMT do not include some systematic errors. This is true also for other materials (not shown). Table 1 shows the obtained lattice constants and related parameters. We also showed values for = 5Ry. Even though we still need other extensive tests, we think that such low cutoff is reliable enough for practical applications. For = 5Ry, we only need basis (when we do not include lo) for Cu.
|Lattice constant.||Bulk Modulus||Cohesive|
|LAPW111Ref. Khein et al. (1995)||6.65||192||—|
|expt.111Ref. Khein et al. (1995)||6.81||131||—|
|PAW222Ref. Kresse and Joubert (1999)||5.20||247||—|
|LAPW333Ref. Stixrude et al. (1994)||5.210||245||—|
|expt.333Ref. Stixrude et al. (1994)||5.417||172||—|
|bcc Li||6.347(6.341)||15.3(15.5)||-.124(-.123) 777The difference from PAW may be because it uses the non-polarized Li atom as reference. If we do so, we have -.150(-.150).|
|PAW222Ref. Kresse and Joubert (1999)||6.355||15.0||-.149|
|PP 444Ref. Liborio et al. (2005). PP denotes a pseudopotential method.||7.31||203||—|
|expt.555Ref. Filippi et al. (1994)||7.39||184||—|
|LAPW555Ref. Filippi et al. (1994)||10.62||74||-.587|
|expt.555Ref. Filippi et al. (1994)||10.68||76||-.479|
Fig.4 shows the total energies with different MTO sets combined with different numbers of APWs. Curve ”A” includes just MTOs for O(), sufficient for a crude representation of the valence bands. Also included are Sr([lo]) and Ti([lo]). A large number of APWs is needed to get a good total energy: APWs are needed to converge energy to within 50 mRy of the converged result. ”B” corresponds to an extreme tight-binding basis, consisting of Sr() and Ti() in addition to “A”. (Note that the conduction band is mainly Ti(), and O()). The total energy of the MTO basis alone (no APWs) is rather crude — more than 200 mRy underbound. However, only 25 orbitals (plus 6 for lo’s) are included in this basis. The energy drops rapidly as low-energy APWs are included: adding 40 APWs is sufficient to converge energy to 50 mRy. As more APWs are added, the gain in energy becomes more gradual; convergence is very slow for large . ”C” differs from ”B” only in that Sr() orbital was added. With the addition of these 5 orbitals, the MTO-only basis is already rather reasonable. This would be the smallest acceptable MTO-basis. As in the ”B” basis, there is initially a rapid gain in energy as the first few APWs are added, followed by a progressively slower gain in energy as more APWs are added. ”D” is a standard LMTO minimum basis: orbitals on all atoms. Comparing ”C” or ”E” to ”A” shows that the MTO basis is vastly more efficient than the APW basis in converging the total energy. This is true until a minimum basis is reached. Beyond this point, the gain APWs and more MTOs improve the total energy with approximately the same efficiency, as the next tests show. ”E” is a standard LMTO larger basis: orbitals on Sr and Ti, and on O. Comparing ”D” and ’F” shows that the efficiency of any one orbital added to to the standard MTO minimum basis is rather similar in the APW and MTO cases. Thus, increasing the MTO basis from 51 to 81 orbitals in the MTO basis lowers the energy by 33 mRy; adding 33 APWs to the minimum basis (”D”) lowers the energy by 36 mRy. ”F” enlarges the MTO basis still more, with Sr: , Ti: , O: . Also local orbitals are used to represent the high-lying states Ti([lo]) and O([lo],[lo]). For occupied states, these orbitals have little effect for total energy as in the case of Cu.
In conclusion, we have implemented the PMT method whose basis set consists of the APWs together with the MTOs which are localized in real space. This idea is consistent with the nature of the eigenfunctions in solids; they can be localized as atoms or extended as PW. This method combines advantages of LAPW/PAW and LMTO. We have implemented force calculations, but no results are shown here. One of the advantage is in its flexibility. As an example, in order to treat Cu impurity in Si, it will be possible to choose very low just responsible for empty regions because MTOs for Cu and MTOs for Si span its neighbors already very well. Convergence is easily checked by changing (much simpler than LMTO). In future, we can use the PMT method to make a natural division of the Kohn-Sham Hamiltonian into localized blocks and extended blocks, instead of the energy windows method for the maximally localized Wanner functions Souza et al. (2001). The problem of large should be solved to make PMT more efficient.
This work was supported by ONR contract N00014-7-1-0479. We are also indebted to the Ira A. Fulton High Performance Computing Initiative.
- Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
- Andersen (1975) O. K. Andersen, Phys. Rev. B 12, 3060 (1975).
- Soler and Williams (1989) J. M. Soler and A. R. Williams, Phys. Rev. B 40, 1560 (1989).
- Blochl (1994) P. E. Blochl, Phys. Rev. B 50, 17953 (1994).
- Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
- Liu et al. (1994) A. Y. Liu, D. J. Singh, and H. Krakauer, Phys. Rev. B 49, 17424 (1994).
- (7) M. Mesfessel and M. van Schilfgaarde, ’NFP manual 1.01 Oct 10,1997’. NFP is previous to the current LMTO package lmf maintained by M. van Schilfgaarde.
- Methfessel et al. (2000) M. Methfessel, M. van Schilfgaarde, and R. A. Casali, in Lecture Notes in Physics, edited by H. Dreysse (Springer-Verlag, Berlin, 2000), vol. 535.
- Bott et al. (1998) E. Bott, M. Methfessel, W. Krabs, and P. C. Schmidt, J. Math. Phys. 39, 3393 (1998).
- Methfessel and van Schilfgaarde (1993) M. Methfessel and M. van Schilfgaarde, Phys. Rev. B 48, 4937 (1993).
- E. Sjostedt (2000) D. J. S. E. Sjostedt, L. Nordstrom, Solid State Communications 114, 15 (2000).
- von Barth and Hedin (1972) U. von Barth and L. Hedin, J. Phys. C 5, 1692 (1972).
- S. H. Vosko and Nusair (1980) L. W. S. H. Vosko and M. Nusair, Can. J. Phys. 58, 1200 (1980).
- Khein et al. (1995) A. Khein, D. J. Singh, and C. J. Umrigar, Phys. Rev. B 51, 4105 (1995).
- Stixrude et al. (1994) L. Stixrude, R. E. Cohen, and D. J. Singh, Phys. Rev. B 50, 6442 (1994).
- Liborio et al. (2005) L. M. Liborio, C. G. Sanchez, A. T. Paxton, and M. W. Finnis, Journal of Physics: Condensed Matter 17, L223 (2005), URL http://stacks.iop.org/0953-8984/17/L223.
- Filippi et al. (1994) C. Filippi, D. J. Singh, and C. J. Umrigar, Phys. Rev. B 50, 14947 (1994).
- Souza et al. (2001) I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 65, 035109 (2001).