Construction and solution of a Wannier-functions based Hamiltonian in the pseudopotential plane-wave framework for strongly correlated materials
Ab initio determination of model Hamiltonian parameters for strongly correlated materials is a key issue in applying many-particle theoretical tools to real narrow-band materials. We propose a self-contained calculation scheme to construct, with an ab initio approach, and solve such a Hamiltonian. The scheme uses a Wannier-function-basis set, with the Coulomb interaction parameter obtained specifically for these Wannier functions via constrained Density functional theory (DFT) calculations. The Hamiltonian is solved by Dynamical Mean-Field Theory (DMFT) with the effective impurity problem treated by the Quantum Monte Carlo (QMC) method. Our scheme is based on the pseudopotential plane-wave method, which makes it suitable for developments addressing the challenging problem of crystal structural relaxations and transformations due to correlation effects. We have applied our scheme to the “charge transfer insulator” material nickel oxide and demonstrate a good agreement with the experimental photoemission spectra.
The problem of inter-electron Coulomb interaction in solids is a many-particle problem which cannot be solved without major approximations. The most successful and widely used approach is electronic density functional theory (DFT) Kohn-Sham (); jones () —within the local density approximation (LDA) or generalized gradient approximation (GGA), where all electrons feel the same one-particle potential defined by the electronic density distribution in the system. For wide-band materials, where the kinetic energy term dominates and the inter-electron Coulomb interaction can be treated in an average way with an energy (time) independent potential, this approach works very well. However for narrow-band materials this is not the case. The Coulomb interaction dominates over the kinetic energy, leading to strong correlations between electrons. Hence a static one-electron potential, as in DFT, becomes a bad approximation. Nevertheless, DFT can still give good results in such systems for some static integral properties. Moreover, DFT calculations can serve as a starting point for more sophisticated approaches designed to treat strongly correlated systems.
The physics of strongly correlated systems was historically studied via solution of model Hamiltonians, such as the Hubbard Hubbard () and Anderson models And (). In such Hamiltonians, the Coulomb interaction term is defined using a set of localized atomic-like orbitals centered on atomic sites. While the kinetic energy is invariant with the choice of the wave-functions-basis set of the model, the Coulomb interaction term significantly depends on the specific form of the atomic-like orbitals.
Wannier functions Wannier () are defined as Fourier transforms of Bloch functions, from wavevector to real space. They are considered nowaday as an optimal choice of basis set to construct model Hamiltonians, Wannier-DMFT1 (); Wannier-DMFT2 () because they have a form of atomic centered localized orbitals and represent a complete basis set for the Bloch functions Hilbert space. However, Wannier functions are not uniquely defined and one needs to impose some additional conditions to make them unique.
Another source of uncertainty in the construction of the Coulomb interaction term is the value of the onsite Coulomb repulsion parameter . Attempts to simply determine it via the integral of the Coulomb potential multiplied by the squares of the orbital wave functions Solovyev (); Ferdi () gave values 2-3 times larger than the experimental estimates. This is due to the neglect of strong relaxation and screening effects. The latter effects can be calculated using perturbation theory, but the results strongly depend on the choice of the screening channels and the number of higher and lower lying states included in the expansion series.
DFT can provide not only good data for the kinetic energy terms in the model Hamiltonians, but also gives a practical alternative way to calculate the Coulomb parameter using constrained DFT calculationsAG (); ConstrainPickett (); Coco (); nakamura (). In the constrained calculations, the DFT equations are solved with a fixed occupancy of the localized orbitals, and is defined as a derivative of the orbital energy with respect to its occupancy.
After its construction, the Hamiltonian for strongly correlated electrons needs to be solved in a way that should be as close to exact as possible. Recently Dynamical Mean-Field Theory (DMFT) vollha93 (); pruschke (); georges96 (); kotliar () became a very popular tool to describe strongly correlated materials, especially when the numerically exact Quantum Monte Carlo (QMC)Hirsh () method is used as a solver for the effective impurity problem. While neglecting, in its standard version, the inter-site correlations, DMFT fully accounts for the local dynamics. When the system is not very close to an ordering-disordering transition, single-site DMFT usually provides a satisfactory agreement between calculated and experimental spectra.
Various methods are available to carry out DFT electronic structure calculations, which are based on different approximations for the wave-functions expansion. Any one of those methods may be used, in principle, to construct and solve the Hamiltonian for strongly correlated materials, if only the electronic spectral functions are needed as the results of the calculations. However, Coulomb correlation effects can lead to strong renormalization of the electron-lattice coupling and hence to complicated phase transitions due to lattice distortions. To describe such effects, lattice relaxation should be taken into account explicitly.
The pseudopotential plane-wave method for first-principles DFT calculations is well suited and widely used to determine lattice relaxations in solidspseudo () – in recent years this has also been extended to DFT plus onsite Coulomb interaction (DFT + ) calculations.Coco (); Fabris05 () It would be desirable to use this method also as a basis to develop a self-contained scheme to determine, via DMFT and Wannier functions, the properties of strongly correlated materials. So far, DMFT computations with Wannier functions have been implemented either using Muffin-Tin-Orbital (LMTO) methods, with Wannier functions constructed using the N-order muffin-tin-orbital scheme,Pavarini04 () or more recently using some mixed-basis methods, with maximally localized Wannier functions,Wannier-DMFT2 () and some an adjustable parameter value.
In the present work, we propose a calculation scheme where all the terms of the Hamiltonian are generated, in a consistent way, within the ab initio pseudopotential plane-wave framework. Starting from the DFT pseudopotential-plane-waves method, we construct atomic-centered Wannier functions and produce the Hamiltonian kinetic-energy term in the basis of the Wannier functions. Then the value of the Coulomb parameter for electrons in these Wannier functions is calculated via constrained DFT calculations. At the last stage, the DMFT-QMC method is used to solve the Hamiltonian. Recently, the electronic structure of nickel oxide was successfully described within the DMFT methodKunes07 (). In the current article NiO plays the role of a well known and already well described test system. The goal is to demonstrate the ability of our new method to reproduce the electronic structure of real system with a level of agreement comparable to previous LMTO-based DMFT calculations. Calculated spectral functions are compared with recent NiO DMFT-calculation resultsKunes07 () and with the experimental photoemission spectra.
ii.1 Wannier functions
Wannier functions (WFs) are defined as Fourier transforms of Bloch functions : Wannier ()
where is the lattice translation vector, the band number and the reciprocal lattice vector. WFs are not uniquely defined because, in the single band case, there is a freedom of choice of the phases of the Bloch functions, , as a function of , and in the multiband case, any set of orthogonal linear combination of Bloch functions could be used in (1). The uncertainty in the WF’s definition corresponds to a freedom of choice for a unitary transformation matrix
One of the most commonly used approach to generate WFs was proposed by N. Marzari and D. Vanderbilt MarzariVanderbilt (). They use a condition of maximum localization of the WFs, that results in a variational procedure for the matrix . As an initial step before the variational process, a set of trial localized orbitals in the form of atomic orbitals was chosen and projected onto the subspace of Bloch functions. Later KuRosner () it was shown that this initial guess for the WFs of transition-metal oxides is usually so good that the variational procedure can be dropped and the projection of the trial orbitals onto the subspace of Bloch functions can be used to define the unitary transformation matrix .
In the present work, we employ the pseudopotential method and a plane wave basis set. Hence, site centered pseudoatomic orbitals were chosen as a set of trial orbitals.
Nonorthogonalized approximations to the WFs in the direct and reciprocal space are calculated as projection of the pseudoatomic orbitals onto a subspace of Bloch functions that is defined by setting an energy interval or some band numbers :
In the plane waves basis, Bloch functions and Bloch sums of pseudoatomic orbitals can be decomposed as:
where is a combined index ( is the index for the atom, are the orbital and magnetic quantum numbers, respectively, and is the spin projection).
Using Eq. (4), together with these decompositions, one obtains:
To generate orthogonalized WFs, one should define the overlap matrix:
Orthogonalized Wannier functions are then obtained as:
For practical use—which includes generating the hamiltonian matrix for the kinetic energy term in the model Hamiltonian and performing constrained LDA/GGA calculation, one needs to compute the matrix elements of a given operator in the basis of the WFs. This can be conveniently done in reciprocal space.
The matrix elements of the one-electron Hamiltonian in reciprocal space are defined as:
where is the eigenvalue of the one-electron Hamiltonian for band .
In real space, the Hamiltonian matrix reads:
The Wannier functions occupancy matrix is given by:
where is the step function and is the Fermi energy.
The transformation from plane waves to the Wannier functions basis is determined by Eqs. (7-II.1) and for the matrix elements by Eqs. (13) and (14). It can be convenient also to back transform from WFs to the plane-waves basis. For example, if some constrained potential has a diagonal form in the WFs basis, , it can be written in the plane-waves basis as:
The localized Wannier functions are used as a basis set to define the Hamiltonian for strongly correlated materials:
Here is the creation (annihilation) operator for an electron in the state defined by the Wannier function , is the occupancy operator for this state, and is a correction term to avoid double-counting for the Coulomb interaction that is already taken into account in DFT in an averaged way.
The first Hamiltonian term, on the right-hand side of Eq.(16), corresponds to the kinetic energy, and the matrix elements are defined by Eq.(13). The Coulomb interaction is described by the second term on the right-hand side of Eq.(16). The Coulomb matrix can be nonzero, in the most general case, for all orbitals. However, normally one takes into account only interactions between electrons in WFs having the symmetry of the transition-metal-ion or orbitals. The corresponding matrix elements can be obtained using only two parameters: the direct Coulomb interaction parameter and the exchange Coulomb interaction (Hund) parameter (see for details Ref. ani97, ). The double-counting correction is taken as:
Here we focus on a -electron system and use for the double-counting potential:
and the fact that one-electron energies in DFT are derivatives of the total energy with respect to the orbitals occupancy
ii.3 Constrained DFT calculations for the Coulomb interaction parameter
In a rigorous way, the Coulomb interaction parameter for electrons in a state described by Wannier function should be calculated via the integral:
where the screened Coulomb interaction is defined via the operator equation:
with the bare Coulomb interaction, , and the polarization operator :
However various attemptsSolovyev (); Ferdi () to use Eqs. (21-23) to calculate the Coulomb interaction parameter gave large variations in the resulting values. This is due to the many possibilities in choosing the channels of the screening via the definition of a set of occupied and unoccupied states in Eq.(23). It either gives a strong underestimation of the values, when one includes into the screening channels transitions between -states (that is clearly unphysical because electrons cannot screen themselves), or too large values of , when only states other than -states are included in the sums in Eq.(23), but the summation over higher and lower energy states is not extended far enough.
Another way to determine the Coulomb interaction parameter is to use constrained DFT calculations,AG (); ConstrainPickett (); Coco (); nakamura () where the screening and relaxation effects are taken into account explicitly in a self-consistent procedure. If one uses the assumption that the contribution to the DFT total energy from the interaction between -electrons is given by Eq.(19), then can be calculated as a second derivative of the DFT total energy with respect to the -orbital occupancy:
using Eq.(20) this can be expressed via the first derivative of the one-electron energy :
To use Eq.(25), one needs to perform DFT calculations with a constraint fixing the occupancy of -orbitals to certain values. In practice, this is done via an auxiliary potential in the form of a projection operator acting on -symmetry WFs :
(in reciprocal space this equation will take the form of Eq.(II.1)). One-electron energy can be calculated as diagonal elements of the Hamiltonian matrix, in Eq.(13), and the corresponding occupancy as diagonal elements of the occupation matrix, in Eq.(14):
and the derivative (Eq.(25)) is computed then numerically.
ii.4 Dynamical Mean-Field Theory
The problem defined by the Hamiltonian (16) can be solved by any of the methods developed to treat many-body effects. In the present work we have used Dynamical Mean-Field Theory (DMFT) vollha93 (); pruschke (); georges96 () which was recently found to be a powerful tool to numerically solve multiband Hubbard models.
In DMFT, the lattice problem becomes an effective single-site problem, which has to be solved self-consistently for the matrix self-energy and the local matrix Green function in the WFs basis set:
where is the chemical potential, is the noninteracting part of the Hamiltonian (16):
(in reciprocal space, the matrix elements for can be calculated using Eq.(12)) and is the self-energy in the Wannier functions basis:
The DMFT single-site problem may be viewed as a self-consistent single-impurity Anderson model georges96 (). The corresponding local one-particle matrix Green function can be written as a functional integral georges96 () involving an action where the Hamiltonian of the correlation problem under investigation, including the interaction term with the Hubbard interaction, enters LDADMFT (). The action depends on the bath matrix Green function through
To solve the functional integral of the effective single-impurity Anderson problem, various methods can be used: quantum Monte Carlo (QMC), numerical renormalization group (NRG), exact diagonalization (ED), noncrossing approximation (NCA), etc. (for a brief overview of the methods see Ref. LDADMFT, ). In the present work, the QMC Hirsh () method is used to solve the impurity problem. The real-frequencies single-particle spectral functions are computed using the maximum entropy method.Jarrell96 ()
Iii Calculations for NiO
The calculation scheme described in Section II was applied in the present work to the problem of the classical Mott insulator nickel oxide, NiO. A DMFT study of the NiO electronic structure was done recentlyKunes07 () with LMTO-based Wannier functions and here we use previous work as a reference. Actually, in the Zaanen-Sawatzky-Allen classification,ZSA () NiO is normally considered as a ’charge transfer insulator’ where the minimal energy excitation across the gap happens between occupied oxygen -states and unoccupied transition metal -states. This fact makes it crucial to include in the calculation not only partially filled -states (as is usually done for strongly correlated materials), but also explicitly take into account oxygen -orbitals.
The pseudopotential plane-waves method, as implemented in the Quantum-ESPRESSO ESPRESSO () package, was used in the present work. The calculations were performed within the local density approximation to density-functional theory, using Vanderbilt ultrasoft pseudopotential Vanderbilt () in the Rappe-Rabe-Kaxiras-Joannopoulos form.RRKJ () For the lattice parameter of NiO (rock-salt structure), we used the experimental value Å. A kinetic-energy cutoff of 40 Ry was employed for the plane-wave expansion of the electronic states. The integrations in reciprocal space were performed using a (4,4,4) Monkhorst-Pack MonkhorstPack () k-point grid.
The Wannier functions, introduced in Section II, are defined by the choice of the Bloch functions Hilbert space and by a set of trial localized orbitals. To show how this choice influences the results, we performed two calculations for the WFs. One with only bands formed by the Ni states included in the Hilbert functional space and the second one where this space was extended by inclusion also of oxygen -bands. In both cases, the projection was done on pseudoatomic wave functions.
The results obtained without and with the -bands are presented in Fig.1 and Fig.2, respectively. In both cases, the -bands obtained in the pseudopotential calculations are exactly reproduced by the bands obtained diagonalizing the k-space Hamiltonian matrix (Eq.(12)). However the spatial distributions of the corresponding Wannier functions are different in the two cases. To demonstrate that, spatial distribution of -like Wannier function charge density was plotted (see Fig.1 and Fig.2). Omission of the oxygen bands in the first calculation results in the appearance of a contribution from orbitals on neighboring oxygen atoms in the Wannier function obtained by projection of the Ni-atom orbitals (see Fig.1). Extending the Hilbert space with explicit inclusion of oxygen bands results in a Wannier function that is nearly a pure -orbital (see Fig.2). In the following, we will use WF obtained with the full sets of bands and pseudo-atomic orbitals including all Ni and O states.
The partial densities of states calculated from the WF hamiltonian are shown in Fig.3. One can distinguish the occupied oxygen -band and, at higher energy, the partially occupied -band. The sub-band is nearly filled, while the Fermi level is located essentially in the middle of the sub-band.
The calculation of Coulomb interaction parameter was done using straightforward constrained DFT procedureFerdi (). We performed a self-consistent calculation in a supercell containing two formula units. The constrained potential (Eq.(26)) on the first Ni atom was positive and equal to 0.1 Ry, and for the second Ni atom it had the same magnitude and opposite sign. The Wannier functions occupations and energies were calculated and the value was then evaluated as described in Sec. II, using Eqs. (25-27). The resulting value for the Coulomb interaction parameter was found to be 6.6 eV. While this is smaller than the value of 8 eV obtained from constrained LDA calculations for linear-muffin-tin orbitals,AZA () it is larger than the value of of 4.5 eV obtained from linear-response-theory calculations for pseudo-atomic orbitals.Coco ()
The interacting Hamiltonian, constructed from the pseudopotential plane-waves calculations, was solved by DMFT-QMC simulations. In the QMC simulations, the inverse temperature value was eV (), and we used 80 time slices and Monte Carlo sweeps. Other DMFT+QMC-calculation details are the same as in Kunes07 (). The DMFT calculations result in strong changes in the spectral functions (see Fig.4) compared to the LDA densities of states (Fig.3), consistent with recent DMFT calculations based on the LMTO method.Ren06 (); Kunes07 () First of all, the correlated spectra show a large-gap insulator in contrast to the metallic LDA solution. The half filled band is now split into an empty upper-Hubbard band and an occupied lower-Hubbard band. In addition, with explicit inclusion of the oxygen band, the lower Hubbard band is found to largely overlap in energy with the oxygen -band and to strongly hybridize with it. Furthermore, although the band remains occupied, it is shifted down in energy so that it also overlaps with the oxygen band.
In Fig.4, we present a comparison of the experimental photoemission (direct and inverse) spectrasaw84 () with the spectral functions obtained from our DMFT pseudopotential calculations. The photoemission spectra measured with 120 eV photons reflect more Ni states and those with 66 eV have more oxygen contributions, due to the difference in the cross section values. One can observe a good agreement between the experimental and calculated spectra. The calculations reproduce a large band gap as well as remarkable redistribution of spectral weight from the top of the valence band to lower energies, going from the 120 eV to the 66 eV spectrum. Comparison with the theoretical orbital-resolved result in Fig.4, indicate that this is consistent with a decreasing -states spectral density contribution and an increasing -states contribution to the spectrum.
We have presented a computational scheme for strongly correlated materials where, in the framework of the pseudopotential plane-waves method, atomic-centered Wannier functions are calculated and the Hamiltonian matrix elements are evaluated in the Wannier-functions basis. Then, via constrain DFT calculations, the Coulomb interaction parameter is obtained for electrons in these Wannier functions, allowing us to generate in a consistent way all necessary components of the Hamiltonian. This Hamiltonian is solved by Dynamical-Mean-Field Theory with the numerically exact Quantum Monte-Carlo method as effective-impurity-problem solver. We have applied this scheme to nickel oxide and demonstrate good agreement with the experimental photoemission data.
This work was supported by the Russian Foundation for Basic Research under the grant RFFI 07-02-00041. We are grateful to Jan Kunes and Dieter Vollhardt for helpful discussions and for supplying the QMC-DMFT computer code. We also acknowledge support for this work by the Light Source Theory Network, LighTnet, of the EU.
- (1) W. Kohn and L.J. Sham, Phys. Rev. 140, A1133 (1965).
- (2) R. O. Jones, O. Gunnarsson, Rev. Mod. Phys. 61, 689 (1989).
- (3) J. Hubbard, Proc. Roy. Soc. (London) 276, 238 (1963).
- (4) P. W. Anderson, Phys. Rev. 124, 41 (1961).
- (5) G. H. Wannier, Phys. Rev. 52, 191 (1937).
- (6) V. I. Anisimov, D. E. Kondakov, A. V. Kozhevnikov, I. A. Nekrasov, Z. V. Pchelkina, J. W. Allen, S.-K. Mo, H.-D. Kim, P. Metcalf, S. Suga, A. Sekiyama, G. Keller, I. Leonov, X. Ren, and D. Vollhardt, Phys. Rev. B 71, 125119 (2005).
- (7) F. Lechermann, A. Georges, A. Poteryaev, S. Biermann, M. Posternak, A. Yamasaki, and O. K. Andersen, Phys. Rev. B. 74, 125120 (2006).
- (8) I. V. Solovyev and M. Imada, Phys. Rev. B. 71, 045103 (2005).
- (9) F. Aryasetiawan, K. Karlsson, O. Jepsen, and U. Schonberger, Phys. Rev. B. 74, 125106 (2006).
- (10) See, e.g., S. Fabris, G. Vicario, G. Balducci, S. De Gironcoli, and S. Baroni, J. Phys. Chem. B 109, 22860 (2005); G. Trimarchi and N. Binggeli, Phys. Rev. B 71, 035101 (2005); and references therein.
- (11) V. I. Anisimov and O. Gunnarsson, Phys. Rev. B. 43, 7570 (1991).
- (12) W. E. Pickett, S. C. Erwin, and E. C. Ethridge, Phys. Rev. B. 58,1201 (1998).
- (13) M. Cococcioni and S. de Gironcoli, Phys. Rev. B. 71, 035105 (2005).
- (14) K. Nakamura, R. Arita, Y. Yoshimoto, and Sh. Tsuneyuki, Phys. Rev. B. 74, 235113 (2006).
- (15) D. Vollhardt, 1993, in Correlated Electron Systems edited by V. J. Emery (Singapore: World Scientific), p. 57.
- (16) Th. Pruschke, M. Jarrell, and J. K. Freericks, Adv. in Phys. 44, 187 (1995).
- (17) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
- (18) G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet and C. A. Marianetti, Rev. Mod. Phys. 78, 865 (2006).
- (19) R. M. Fye and J. E. Hirsch, Phys. Rev. B. 38, 433 (1988).
- (20) M. Jarrell and J. E. Gubernatis, Phy. Rep. 269, 133 (1996).
- (21) Richard M. Martin, Electronic structure, Basic Theory and practical methods, Cambridge University press 2004; D. Marx and J. Hutter, in Modern Methods and Algorithms of Quantum Chemistry, pp. 301 V449 (1rst edition, paperback) or pp. 329 V477 (second edition, hardcover), Editor: J. Grotendorst (John von Neumann Institute for Computing, Forschungszentrum Julich 2000); M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias, and J. D. Joannopoulos, Rev. Mod. Phys. 64, 1045 (1992).
- (22) E. Pavarini, S. Biermann, A. Poteryaev, A. I. Lichtenstein, A. Georges, and O. K. Andersen, Phys. Rev. Lett. 92, 176403 (2004).
- (23) N. Marzari and D. Vanderbilt, Phys. Rev. B. 56, 12847 (1997).
- (24) Wei Ku, H. Rosner, W. E. Pickett, and R. T. Scalettar, Phys. Rev. Lett. 89, 167204 (2002).
- (25) V. I. Anisimov, F. Aryasetiawan, and A. Lichtenstein, J. Phys.: Condens. Matter 9, 767 (1997); in correlated transition-metal compounds in general, a good approximation for the J of the transition metal is eV.
- (26) By electrons we imply electrons in Wannier orbitals of symmetry.
- (27) K. Held, I. A. Nekrasov, N. Blümer, V. I. Anisimov, and D. Vollhardt, Int. J. Mod. Phys. B 15, 2611 (2001).
- (28) J. Zaanen, G. A. Sawatzky, J. W. Allen, Phys. Rev. Lett. 55, 418 (1985).
- (29) S. Baroni, A. Dal Corso, S. de Gironcoli, P. Giannozzi, C. Cavazzoni, G. Ballabio, S. Scandolo, G. Chiarotti, P. Focher, A. Pasquarello, K. Laasonen, A. Trave, R. Car, N. Marzari, and A. Kokalj, http://www.pwscf.org/.
- (30) A. M. Rappe, K. M. Rabe, E. Kaxiras, and J. D. Joannopoulos, Phys. Rev. B 41, 1227 (1990).
- (31) D. Vanderbilt, Phys. Rev. B 41, 7892 (1990).
- (32) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
- (33) V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys. Rev. B. 44, 943 (1991).
- (34) G. A. Sawatzky and J. W. Allen, Phys. Rev. Lett. 53, 2339 (1984).
- (35) A. Kokalj, Comp. Mater. Sci., 28, 155 (2003). Code available from http://www.xcrysden.org/.
- (36) X. Ren, I. Leonov, G. Keller, M. Kollar, I. Nekrasov, and D. Vollhardt, Phys. Rev. B 74, 195114 (2006).
- (37) J. Kunes, V. I. Anisimov, A. V. Lukoyanov, and D. Vollhardt, Phys. Rev. B 75, 165115 (2007).