Charge self-consistent many-body corrections using optimized projected localized orbitals
Malte Schüler, Oleg E. Peil, Gernot J. Kraberger, Ronald Pordzik, Martijn Marsman, Georg Kresse, Tim O. Wehling, Markus Aichhorn
CCharge self-consistent many-body corrections usingoptimized projected localized orbitals
M Sch¨uler , , O E Peil , , , G J Kraberger , R Pordzik , ,M Marsman , G Kresse , T O Wehling , and M Aichhorn Institut f¨ur Theoretische Physik, Universit¨at Bremen, Otto-Hahn-Allee 1,28359 Bremen, Germany Bremen Center for Computational Materials Science, Universit¨at Bremen, AmFallturm 1a, 28359 Bremen, Germany Centre de Physique Th´eorique, ´Ecole Polytechnique, CNRS, Universit´eParis-Saclay, 91128 Palaiseau, France Department of Quantum Matter Physics, University of Geneva, 24 QuaiErnest-Ansermet, 1211 Geneva 4, Switzerland Materials Center Leoben Forschung GmbH, A-8700 Leoben, Austria Institute of Theoretical and Computational Physics, Graz University ofTechnology, NAWI Graz, 8010 Graz, Austria Faculty of Physics and Center for Computational Materials Sciences,University of Vienna, Sensengasse 8/12, 1090 Wien, Austria
Abstract.
In order for methods combining ab-initio density-functional theoryand many-body techniques to become routinely used, a flexible, fast, and easy-to-use implementation is crucial. We present an implementation of a generalcharge self-consistent scheme based on projected localized orbitals in the ProjectorAugmented Wave framework in the Vienna Ab-Initio Simulation Package ( vasp ).We give a detailed description on how the projectors are optimally chosen and howthe total energy is calculated. We benchmark our implementation in combinationwith dynamical mean-field theory: First we study the charge-transfer insulatorNiO using a Hartree-Fock approach to solve the many-body Hamiltonian. Weaddress the advantages of the optimized against non-optimized projectors andfurthermore find that charge self-consistency decreases the dependence of thespectral function – especially the gap – on the double counting. Second, usingcontinuous-time quantum Monte Carlo we study a monolayer of SrVO , wherestrong orbital polarization occurs due to the reduced dimensionality. Using total-energy calculation for structure determination, we find that electronic correlationshave a non-negligible influence on the position of the apical oxygens, and thereforeon the thickness of the single SrVO layer. a r X i v : . [ c ond - m a t . s t r- e l ] N ov harge self-consistent many-body corrections using optimized PLOs
1. Introduction
The advances in the field of nanostructures, wherecontrol is now experimentally possible on an atom-by-atom scale, have given rise to a strong demandfor theoretical simulation tools that are capable ofsimulating complex correlated electron systems suchas heterostructures, clusters, or adatom arrays onsurfaces. There are several successful approaches thatcan be used for that. Among the most widely used areab-initio approaches, notably density functional theory(DFT) and many-body model Hamiltonians.Modern Kohn-Sham DFT, within the local densityapproximation (LDA) [1] or a generalized gradientapproximation (GGA) [2], yields various ground-state properties including crystal structures quitereliably for many materials, essentially without anyambiguous input parameters. The auxiliary single-particle energies can be seen as an estimate of theelectronic quasi-particle energies,[3] sometimes givinggood qualitative agreement. However, (semi)localfunctionals like LDA or GGA are well known tofail to capture the physics of strongly-correlatedmaterials such as Mott insulators, unconventionalsuperconductors, heavy Fermion systems, or Kondosystems.This is why, for treating these systems, modelHamiltonians such as the Hubbard model are usuallyemployed. Generally speaking, these approaches focuson the description of electron correlation phenomenain a minimal low-energy Hilbert space. They areaccessible by a broad set of many-body techniquesincluding dynamical mean-field theory (DMFT)[4]and generalizations thereof. These models naturallydepend on model parameters that are unknown apriori. Keeping the number of parameters low canlead to over-simplified models that fail to explainsufficiently well experimental findings. On the otherhand, keeping a large number of unknown parametersmakes the results ambiguous already from the outset.Therefore, to obtain the best of both worlds,it appears natural to combine the complementarystrengths of both approaches to model correlatedelectron systems in a realistic manner. To achievethis goal, there has been an ongoing effort in thedevelopment of approaches like DFT+DMFT[5, 6] formore than 20 years now. The class of materialsaddressed by DFT+DMFT is very wide and includesMott insulators, correlated metals, superconductors, and magnetic materials.On the DFT side, simulations using projectoraugmented wave (PAW) [7] basis sets turn outto provide a good compromise between accuracyand computational requirements in these systems.Thus, several DFT+DMFT approaches have beenimplemented with PAW basis sets on the DFT side.[8,9] In its most simple formulation, the DMFT isperformed on top of a converged DFT calculation(so-called one-shot DFT+DMFT). For many systems,especially those where the strong correlations cause arearrangement of charges, the success of the methodcan be significantly improved by performing a fullycharge self-consistent calculation. [10, 11, 12, 13, 14]There, the electron density obtained from the DMFTcalculation is fed back to the DFT code and theentire loop is self-consistently solved. There alreadyexist implementations of DFT+DMFT using the PAWformalism that include charge self-consistency. [15, 16,17] Here we describe a fully charge self-consistentimplementation based on optimized projected localorbitals in the Vienna Ab Initio Simulation Package( vasp ).[18, 19, 20] On the DFT side, it is flexibleand easy to use; vasp provides the data necessary toconstruct the low-energy model of the system and offersan interface so that the updated charge density can behanded back. This makes it possible to combine it, e.g.,with any kind of many-body correction beyond DFT insome correlated subspace defined via local projectionoperators without requiring any changes to the DFTcode. On the DMFT side, we present two codes makinguse of this extension of vasp .The paper is organized as follows. In section 2we explain the details of our approach, particularly,the definition of the optimized local projectors andthe way the charge feedback from the many-bodyto the DFT part is implemented. We then presentan application to the testbed material of NiO insection 3, where we first analyze the optimized versusnon-optimized projectors and second demonstrate howfull charge self-consistency affects simulations of theelectronic structure, particularly, in relation with theso-called double-counting problem. In section 4, wereport an application of our scheme to monolayers ofSrVO , where we compare our implementation, in the triqs [21] framework, to the already existing interfacebetween triqs and the DFT code wien2k .[22, 23]Furthermore, for that material, we demonstrate total harge self-consistent many-body corrections using optimized PLOs
2. Details of the implementation
To perform DFT+DMFT calculations one needs a wayto transform between the basis of the Kohn-Sham (KS)states and the localized basis of the subspace used todefine lattice models, e.g., a Hubbard model. Themost obvious way to do this is to perform a unitarytransformation of a subset of KS states to constructa set of local states. This method is at the heart ofvarious types of Wannier function based methods, suchas the maximally-localized Wannier functions.[24]In many cases, when KS bands with different char-acters are strongly entangled, it becomes however diffi-cult to construct a well-defined unitary transformation.In this case it is more advantageous to use projector op-erators which, acting on the KS Hilbert space, projectout only KS states with a desired character. This isthe basis of projected localized orbitals (PLO).[8]In the PLO formalism one starts by definingan orthonormal localized basis set | χ L (cid:105) associatedwith each correlated site, which is typically indexedby local quantum numbers, e.g., L = { l, m, σ, ... } with ( l, m ) and σ being orbital and spin angular-momentum quantum-numbers, respectively. {| χ L (cid:105)} spans a correlated subspace C at each correlated site.Assuming (cid:104) χ L | χ L (cid:48) (cid:105) = δ LL (cid:48) , any operator acting on thisspace can be constructed by projecting onto C using aprojector operator ˆ P C , i.e.,ˆ A imp = ˆ P C ˆ A ˆ P C , (1)ˆ P C = (cid:88) L | χ L (cid:105)(cid:104) χ L | . (2)An arbitrary vector | Ψ (cid:105) of the Hilbert space can bedecomposed in terms of local states by writing ˆ P C | Ψ (cid:105) = (cid:80) L | χ L (cid:105)(cid:104) χ L | Ψ (cid:105) . If we now consider a complete basis | Ψ µ (cid:105) , the projector operator is completely defined byspecifying PLO functions P L,µ ≡ (cid:104) χ L | Ψ µ (cid:105) .As described in detail in [8, 9], the PLO projectorin the PAW framework[7, 18] can be written as P R L,ν ( k ) = (cid:88) i (cid:104) χ R L | φ i (cid:105)(cid:104) ˜ p i | ˜Ψ ν k (cid:105) , (3)where | χ R L (cid:105) are localized basis functions associatedwith each correlated site R , | φ i (cid:105) are all-electron partialwaves, and | ˜ p i (cid:105) are the standard PAW projectors. Inthe following, we will omit the site index R unless itleads to a confusion. The index i stands for the PAWchannel n , the angular momentum quantum number l , and its magnetic quantum number m . | ˜Ψ ν k (cid:105) arepseudo-KS states. For the PAW potentials distributed with vasp , generally a minimum of two channels n foreach angular quantum number are used. For instancefor 3 d elements, the first 3 d channel is placed at theenergy of the bound 3 d state in the atom, and a secondchannel is added a few eV above the bound atomic3 d state. Inclusion of the second channel improvesthe transferability of the PAW potentials and thedescription of the scattering properties greatly. For s and p states in, e.g., 3 d transition metals, the firstchannel usually describes 3 s and 3 p semi-core statesand the second channel is placed somewhere in thevalence regime (4 s and 4 p states) or even above thevacuum level. Although a description of how theprojectors have been obtained is stored in more recentPAW potential files, it is not always straightforward forthe user to identify the best suitable PLO projectors.Specifically, various choices are possible for | χ L (cid:105) .[9] Conveniently, one might simply use the all-electron partial waves for a particular PAW channel n , | φ Ln (cid:105) , as the local basis. However, doing so in vasp can lead to a non-optimal projection. For instance, fortransition metals to the left of the periodic table, the d states in the solid might be more or less contractedthan those in the atom, so that the “best” | χ L (cid:105) is alinear combination of the two available 3 d channels.Likewise, projection onto the TM valence s states isoften difficult, because of the presence of semi-core s states in the PAW potential file. Ideally, the usershould not need to make a manual choice of the local | χ L (cid:105) functions.In the following we explain a protocol thatlargely resolves these issues. The first step is notstrictly required, but simplifies the subsequent codingsomewhat. We first construct a set of PAW projectorsand partial waves that are orthogonalized inside thePAW sphere. This can be achieved by diagonalizing theall-electron one-center overlap matrix O nn (cid:48) (for eachangular and magnetic quantum number L ), O nn (cid:48) = (cid:104) φ Ln | φ Ln (cid:48) (cid:105) , (4)which gives the eigenvalues λ n and the eigenvectormatrix U ,Λ = diag( λ , . . . , λ n ) , (5)Λ = U † OU. (6)The new partial waves and the corresponding PAWprojectors can now be defined as linear combinationsof the original ones: | ξ Ln (cid:105) = 1 √ λ n (cid:88) n (cid:48) U nn (cid:48) | φ Ln (cid:48) (cid:105) , (7) (cid:104) ˜ β Ln | = (cid:112) λ n (cid:88) n (cid:48) U † n (cid:48) n (cid:104) ˜ p Ln (cid:48) | . (8)This new set of projectors and partial waves have theimportant property that the all-electron partial waves | ξ Ln (cid:105) are orthogonal to each other and that they are harge self-consistent many-body corrections using optimized PLOs within a chosen (usersupplied) energy window [ ε P min , ε P max ]. Specifically, wediagonalize a matrix M nn (cid:48) = (cid:88) ν k (cid:104) ˜ β Ln | ˜Ψ ν k (cid:105)(cid:104) ˜Ψ ν k | ˜ β Ln (cid:48) (cid:105) , (9)where the sum is restricted to states with energies ε ν k ∈ [ ε P min , ε P max ]. We then pick the eigenvector υ n with the largest absolute eigenvalue to construct localbasis functions and corresponding PAW projectors, | χ L (cid:105) = (cid:88) n υ n | ξ Ln (cid:105) , (10) (cid:104) ˜ π L | = (cid:88) n υ ∗ n (cid:104) ˜ β Ln | , (11)with projector functions given by ‡ P L,ν ( k ) = (cid:104) ˜ π L | ˜Ψ ν k (cid:105) . (12)We note in passing that both steps could be combinedinto a single computational step, by diagonalization ofa generalized eigenvalue problem involving the overlapmatrix O and the matrix M evaluated using theoriginal set of projectors and partial waves.Note that since the all-electron partial waves donot form an orthonormal basis set between differentPAW spheres and because projection is performedfor a subset of KS states, the above procedure willusually produce non-normalized local states. However,normalized local states can be constructed if weorthonormalize the projector functions as follows, O LL (cid:48) = (cid:88) ν k P L,ν ( k ) P ∗ ν,L (cid:48) ( k ) , (13) P L,ν ( k ) ← (cid:88) L (cid:48) O − LL (cid:48) P L (cid:48) ,ν ( k ) . (14)On a technical level, the optimal projectors areconstructed internally in vasp according to (4-12) foreach sites R given an energy window [ ε P min , ε P max ]. Thelast orthonormalization step (14) is done externally aspost-processing. This allows one to choose a differentenergy window for the summation in (13) in orderto fine-tune the degree of localization of the impurityWannier functions.The projection scheme has been used quite widelyin the solid state community,[8, 9] even in combinationwith VASP. When using VASP, these calculationswere often not based on a solid mathematicalfoundation. VASP used to project onto all available ‡ The projectors P L,ν ( k ) according to (12) containing all phasefactors are written by vasp to a file called LOCPROJ. PAW projectors with a certain angular and momentumquantum number lm , summed the resulting densities,averaged the phase factors and intensities over all lm projectors, and finally wrote the results to afile ( PROCAR ). Such a prescription, in particular, theaveraging of the intensities and phase factors was done ad hoc without a proper mathematical prescription.This can result in unsatisfactory results if the twoprojectors span a completely different subspace, forinstance Ni 3 d and 4 d states, or transition metal semi-core p states and valence p states. We will demonstratethis issue for NiO in section 3.1. A general extension of DFT Kohn-Sham equations to amany-body problem based on the Hubbard model canbe formulated using the total-energy functional[25, 10] E [ n, ˆ G C ] = 1 β (cid:88) i ω n k Tr { ˆ G ( k , i ω n ) ˆ H KS ( k ) } + E H [ n ] + E xc [ n ] − (cid:90) d r ( v xc + v H ) n ( r )+ E corr [ ˆ G C ] − E dc [ ˆ G C ] , (15)where ˆ G ( k , i ω n ) is the Green’s function containingmany-body effects, and ˆ H KS ( k ) = (cid:80) ν | Ψ ν k (cid:105) ε ν k (cid:104) Ψ ν k | is the KS Hamiltonian. The charge density, n ( r ),and Green’s function in the correlated subspace,ˆ G C ( k , i ω n ), are both related to ˆ G ( k , i ω n ) by n ( r ) = 1 β Tr (cid:88) i ω n k (cid:104) r | ˆ G ( k , i ω n ) | r (cid:105) , (16)ˆ G C ( k , i ω n ) = ˆ P C ( k ) ˆ G ( k , i ω n ) ˆ P C ( k ) , (17)where ˆ P C k projects onto correlated localized states,see (2). E corr is the energy contribution of theinteraction term (Hubbard U -term) and E dc is thedouble-counting correction.The first four terms of the above expression formthe usual DFT total energy calculated for the densitymatrix and charge density of the interacting system,which is non-diagonal in this case. It is thus clear thatto get the correct value of the total energy the densitymatrix must be calculated in a self-consistent manner.In the framework of DFT+DMFT [4, 6] theinteracting Green’s function is defined as follows:ˆ G ( k , i ω n ) = (cid:104) (i ω n + µ )ˆ11 − ˆ H KS ( k ) − ˆΣ KS ( k , i ω n ) (cid:105) − , (18)where ˆΣ KS ( k , i ω n ) is obtained by up-folding the localself-energy,ˆΣ KS ( k , i ω n ) = (cid:88) νν (cid:48) | Ψ ν k (cid:105)(cid:104) Ψ ν (cid:48) k |· (cid:88) mm (cid:48) P ∗ ν,m ( k )Σ mm (cid:48) (i ω n ) P m (cid:48) ,ν (cid:48) ( k ) , (19) harge self-consistent many-body corrections using optimized PLOs | χ lm (cid:105) , areΣ mm (cid:48) (i ω n ) = Σ imp mm (cid:48) (i ω n ) − Σ dc mm (cid:48) , (20)which consist of the purely local impurity self-energyˆΣ imp (i ω n ), obtained from an impurity solver, and thedouble counting ˆΣ dc . Using the PLO functions we canwrite this Green’s function in terms of its KS matrixelements G νν (cid:48) ( k , i ω n ) = (cid:104) (i ω n + µ − ε ν k )ˆ11 − ˆΣ KS ( k , i ω n ) (cid:105) − νν (cid:48) . (21)If we define an energy window [ ε C min , ε C max ] whichselects a subset W C of KS states affected bycorrelations, the interacting charge density can bewritten as n ( r ) = (cid:88) k (cid:88) ν / ∈W C f ν k (cid:104) r | Ψ ν k (cid:105)(cid:104) Ψ ν k | r (cid:105) + (cid:88) k (cid:88) νν (cid:48) ∈W C (cid:104) r | Ψ ν k (cid:105) N νν (cid:48) ( k ) (cid:104) Ψ ν (cid:48) k | r (cid:105) , (22)where we use a non-diagonal density matrix N νν (cid:48) ( k ) = 1 β (cid:88) i ω n G νν (cid:48) ( k , i ω n ) , (23)and a diagonal density matrix formed by the KSoccupation numbers f ν k . Note that the chemicalpotential µ here is generally different from the KSchemical potential µ KS . The subset of KS states forthe optimization of the projectors W P and KS statesaffected by correlations W C can be chosen to be thesame but one does necessarily not have to do so.From a practical point of view it is convenientto split the charge density into a DFT part and acorrelation-induced part, n ( r ) = n DFT ( r ) + ∆ n ( r ) , (24) n DFT ( r ) = (cid:88) ν k f ν k (cid:104) r | Ψ ν k (cid:105)(cid:104) Ψ ν k | r (cid:105) (25)where ∆ n ( r ) is the correlation correction,∆ n ( r ) = (cid:88) k νν (cid:48) ∈W C (cid:104) r | Ψ ν k (cid:105) ∆ N νν (cid:48) ( k ) (cid:104) Ψ ν (cid:48) k | r (cid:105) , (26)∆ N νν (cid:48) ( k ) = N νν (cid:48) ( k ) − f ν k δ νν (cid:48) , (27)with the last quantity having the important property (cid:88) ν ∈W C (cid:88) k ∆ N νν ( k ) = 0 . (28)Such a definition of the new charge density is con-venient because it ensures charge neutrality betweenDFT iterations.In a similar way, the total energy can be splitinto a DFT part calculated using the correlated charge density given by (22) and a correlation part, E = E DFT [ n ]+ (cid:88) k (cid:88) ν ∈W C ∆ N νν ( k ) ε ν k + E corr [ ˆ G C ] − E dc [ ˆ G C ] . (29)Thus, in order to calculate the total energy andto obtain the charge density for the next KS iteration,only two quantities need to be calculated after a DMFTiteration: the density-matrix correction ∆ N νν ( k ) andthe interaction energy (including the double-countingterm) E corr − E dc .In this particular vasp implementation, weuse ∆ N νν ( k ) to obtain the natural orbitals by atransformation V given by diagonalizing the totalcorrelated density matrix, f (cid:48) ν k δ νν (cid:48) = (cid:88) µµ (cid:48) V νµ [ f µ k δ µµ (cid:48) +∆ N µµ (cid:48) ( k )] V ∗ µ (cid:48) ν (cid:48) , (30) | Ψ (cid:48) ν k (cid:105) = (cid:88) µ V νµ | Ψ µ k (cid:105) , (31)and the charge density and the one-electron energy are,then, obtained in the same way as in the normal KScycle, n ( r ) = (cid:88) ν k f (cid:48) ν k (cid:104) r | Ψ (cid:48) ν k (cid:105)(cid:104) Ψ (cid:48) ν k | r (cid:105) , (32)1 β (cid:88) i ω n k Tr { G ( k , i ω n ) H KS ( k ) } = (cid:88) ν k f (cid:48) ν k (cid:104) Ψ (cid:48) ν k | H KS | Ψ (cid:48) ν k (cid:105) . (33)Note that the DFT band energy is calculated usingthe original occupancies f ν k , and the second term in(29) represents essentially a band-energy correctionwhich takes into account the change in the densitymatrix induced by correlations. The occupancies f (cid:48) ν k will deviate from the usual Fermi-Dirac statistics as aresult of particle fluctuations from the occupied non-interacting KS states into unoccupied states.
3. Benchmark for NiO
NiO is a prototypical charge-transfer insulator wherethe electronic band gap on the order of 4 eV arisesfrom a combination of strong local Coulomb repulsion U within the Ni 3 d shell and the charge-transfer energy∆ = (cid:15) d − (cid:15) p , i.e., the difference in on-site energies ofO-2 p and Ni-3 d orbitals. [26] The uppermost valencestates are of hybrid Ni-3 d and O-2 p character, wherethe O-2 p contribution is dominant.DFT in (semi)local approximations like LDA orGGA fails to describe the insulating nature of NiO.Non-spin-polarized LDA and GGA yield a metal, whilespin-polarized calculations of the antiferromagneticallyordered phase of NiO yield an energy gap of ∼ . harge self-consistent many-body corrections using optimized PLOs X k / e V X k N i d w e i g h t opt, N b = 12old, N b = 12opt, N b = 24old, N b = 24 Figure 1.
Left: Ni d weight of the band shown as red dashed lineon the right for paramagnetic NiO. Results using non-optimizedprojectors are shown as dashed lines, those from the optimizedprojectors as solid lines. Orange and blue are results for 12 and24 bands, respectively. established gap.[27] The reason for this failure ofsemilocal DFT to describe the insulating state of NiOis well known to be the insufficient treatment of thestrong local Coulomb interaction in the Ni 3 d shell.NiO has thus become a testbed material for realisticcorrelated electron approaches such as DFT+U [28] orDFT+DMFT.[5, 29, 30, 31] For the example of NiO, we first investigate theproperties of the optimized projectors described insection 2. To this end, we analyze the Kohn-ShamHamiltonian transformed to a localized basis, whichwill later serve as a starting point for including localcorrelation effects on a Hartree-Fock level.In order to treat the Ni d states and ligand O p -states explicitly, we project the paramagnetic Kohn-Sham Hamiltonian of a two atomic unit cell using48 × × k points onto Ni d and O p statesaccording to (12). Therefor, we optimize the projectorsin an energy window around the Fermi level given by ε P min = − . ε P max = 1 . d states) and orthonormalize accordingto (14) for all available energies. We perform a Gram-Schmidt procedure to orthogonally complement the Ni d and O p states with additional states to end up withthe same number of localized states as Kohn-Shamstates.In figure 1 we analyze the orbital properties ofselected eigenstates of the resulting Hamiltonian. Theleft panel shows the Ni d weight of a single band(displayed as red dashed line in the right panel) acrossa k-path from the Γ to the X point. For a small numberof unoccupied bands, (total number of bands N b = 12),the non-optimized and optimized projectors give nearlyidentical results. However, by including additionalunoccupied bands ( N b = 24) the two schemes displayhuge differences. While the optimized projectors lead to only minor differences to the case of N b = 12 closeto the Γ point, the non-optimized projectors lead toconsiderably less overall Ni d weight and a stronglydiscontinuous behavior for part of the k pathway. Thisunderlines the advantage of the optimized projectors.The erroneous behavior of the non-optimizedprojectors roots in additional high energy bandshaving Ni-4 d character. Using the standard VASP“projectors”, these high lying valence orbitals showappreciable overlap with the 4 d partial waves (and4 d projectors). VASP used to lump the 4 d and3 d contributions into single values, which causes theissues we have just observed. Without going intomathematical details one can easily understand thisbehavior. The 3 d or 4 d states are automaticallyorthogonal, because the 4 d states possess an additionalnode inside the PAW sphere. By simply adding up thecontributions from the 3 d and 4 d projectors pretendingthat they correspond to a single main quantumnumber, and then performing an orthogonalization, thetotal charge in each d spin channel is one instead of two,explaining the large reduction of the d character usingthe old scheme.The failure of the unoptimized projectors could inprinciple be avoided by simply not including additionalempty bands in the construction of the Wannierfunctions. However, there are important cases wherethis is not possible. First of all, the situation of 4 d states close to the Fermi energy is realized in many latetransition metals. Second, one might be interested inquantities that include transitions to unoccupied statesover large energy scales, such as optical spectroscopy.In these cases, the optimized projectors give reliableresults. Furthermore, from a physical point of view, theprojectors should not depend strongly on the inclusionof unoccupied states. The DFT+U approach improves the LDA and GGAdescription of NiO by supplementing the Ni d statesby local Coulomb interaction which leads to a staticlocal self energy in (20).There are two ways to formulate DFT+U: First,the traditional one, where the Kohn-Sham potentialis directly augmented with a potential arising fromthe Hartree-Fock decoupling of the local Hubbardinteraction. Second, in terms of a charge self-consistentDFT+DMFT scheme, where the DMFT impurityproblem is solved in the Hartree-Fock approximationcalled DFT+DMFT(HF) in the following. Bothschemes are fully equivalent if the augmentation spacesare chosen to be strictly the same and the way spin-polarization is accounted for is the same.However, many DFT+U implementations, includ-ing the one in vasp , work with angular-momentum- harge self-consistent many-body corrections using optimized PLOs / eV A () a) LDA+DMFT(HF) total O , O p Ni d Ni d / eV b) VASP LDA+U Figure 2.
Partial density of states for spin up of Ni d and O p states from a) charge self-consistent DFT+DMFT(HF) with µ dc = 62 . vasp DFT+U with FLL double counting.Ni and Ni refers to the two different Ni atoms in the unitcell arising due to the anti-ferromagnetic ordering. The spin-updensity of the atom Ni is equal to the spin-down density of Ni and vice versa. The oxygen atoms in the unit cell show no spinpolarization, thus, O is equivalent to O . decomposed charge densities, which in general do notrelate to proper projectors in the definitions of the”+U” potentials.[32] Additionally, DFT+DMFT andDFT+U for magnetic materials can work with or with-out spin polarization in the DFT part.In the following, we compare DFT+U asimplemented in vasp to DFT+DMFT(HF) with andwithout charge self-consistency for the testbed materialNiO, and investigate the dependence of the electronicspectra on the double-counting.We fix the double counting in the vasp DFT+Uapproach according to a prescription known as “fullylocalized limit” (FLL)[33] which leads to an orbitalindependent and diagonal version of Σ dc mm (cid:48) in (20), inshort µ dc . For a meaningful comparison of spectralfeatures we choose the double counting potential in theDFT+DMFT(HF) approach such that we reproducethe size of the gap in vasp DFT+U, leading to µ dc =62 . × × U = 8 eV and J = 1 eV as obtained fromconstrained LDA calculations.[28] We perform spinpolarized calculations considering antiferromagneticordering of the two Ni atoms in the unit cell, where forboth, DFT+DMFT(HF) and vasp DFT+U, we onlyconsider spin polarization in the Hartree-Fock part butnot in the DFT part.The respective density of states for spin up arepresented in figure 2. In general, the DOS fromDFT+U and DFT+DMFT(HF) are very similar: agap of about 4 eV is opened between heavily spin-polarized Ni states in the conduction band and stronglyhybridized Ni-O states in the valence band. Onlydetails differ for the two approaches, such as a
DOS 0high56 58 60 62 64 DC / eV / e V a) fcscFLL
56 58 60 62 64 DC / eV b) one-shot A () FLLAMF
Figure 3.
Color-coded density of states from DFT+DMFT(HF)for double-counting potentials µ dc in steps of 0 . µ dc from FLL and AMF are depicted as dashed lines,where AMF in the case of charge self-consistency is 53 . slightly larger spin polarization of Ni states and alarger O contribution at the valence band edge forDFT+DMFT(HF).The double-counting problem poses a seri-ous problem when using many-body corrections inDFT+DMFT approaches predictively, since funda-mental properties such as the single particle gap de-pend on the double counting. Here, we performfully charge self-consistent (fcsc) as well as one-shot DFT+DMFT(HF) calculations for various fixeddouble-counting potentials µ dc to investigate its influ-ence on the spectra and on the gap.Strictly speaking, this approach does not representa proper double-counting scheme on first sight, sincethere seems to be no functional relation betweenthe double-counting potential µ dc and energy E dc .However, one can motivate this approach followingarguments in [34]. Instead of using the standard FLLform for the double-counting correction, one introducesan interaction parameter U (cid:48) in the double countingformulas, which is allowed to be slightly different from U used otherwise in the calculation. In that way, onecan again relate µ dc to a proper energy correction E dc .However, here we are not interested in totalenergies but in a systematic investigation of theinfluence of the double counting on the spectralproperties. That is why we refrain from an explicitintroduction of this parameter U (cid:48) . Furthermore,we want to circumvent here further complicatingambiguities in usual double-counting approaches (e.g.,if one should use the formal or self-consistentoccupation in formulas for E dc ).[30]The resulting total DOS for the fcsc and one-shot calculations are presented color coded in figure 3.In both cases the Ni states in the conduction band(dark blue around 5 eV) are shifted linearly towardsthe Fermi energy as µ dc increases. The O/Ni states harge self-consistent many-body corrections using optimized PLOs µ dc , i.e., in the oppositedirection than the Ni state in the conduction band.For µ dc (cid:46)
58 eV in case of fcsc and µ dc (cid:46)
57 eVin the case of one-shot calculations, these states crossthe Ni states. Then, the conduction band edge doesnot consist of Ni states, which is not in agreementwith experiment. Note that the “around-mean-field”AMF prescription lies in this regime both for one-shotand fcsc calculations and is thus not suitable for NiO.Similarly, the FLL prescription for fcsc calculations isvery close to this regime.For values of the double counting that give thecorrect conduction band character, the gap shrinkswith increasing µ dc . This dependence differs stronglyfor the one-shot and charge self-consistent calculations:the slope of the gap dE g /dµ dc in the physical regimediffers by a factor of nearly 2. Thus, the charge self-consistency does not cure but alleviates the double-counting problem by decreasing the influence of µ dc on the spectrum.
4. Benchmark for SrVO monolayer SrVO presents an example of a correlated transition-metal oxide experiencing a metal-insulator transition(MIT) driven by a dimensional crossover. [35] Inparticular, a monolayer of SrVO grown on a SrTiO substrate is an insulator with a Mott gap of around2 eV. In DFT, the compound is found to be metallicand one has to combine DFT with a many-bodytechnique to achieve the correct insulating solution.The related problem of a double layer of SrVO showsinsulating behavior in a non charge self-consistentDFT+DMFT treatment. [36]Here, we use a monolayer of SrVO to benchmarkthe fully charge-self-consistent (fcsc) DFT+DMFTimplementation of triqs / dfttools [23] within the vasp PLO formalism. The motivation for choosingthis particular system is that electronic correlationsinduce an appreciable charge redistribution, makingthe use of a fcsc DFT+DMFT scheme imperative.Importantly, using the triqs / dfttools framework,we can benchmark the presented vasp interface to theone that is based on the wien2k DFT package,[22] alsoimplemented in triqs / dfttools . Furthermore, wecompare our results to previously published fcsc data,also based on the wien2k code, but using MLWF as . ˚ A . ˚ A Figure 4.
The unit cell of the monolayer of SrVO . The Sratoms are green, the V atom gray and the O atoms red. Thelattice constant in the z -direction is 20 ˚A, i.e., there is about 16 ˚Aof vacuum between the periodic replica of the layers, effectivelygiving an isolated monolayer. correlated basis set.[13]The system is modeled by a free-standingmonolayer of SrVO with the in-plane lattice constantequal to that of SrTiO (3.92 ˚A), simulating thus theepitaxial geometry. The effect of the substrate isneglected. To isolate the individual layers in a periodicunit cell, a vacuum layer of about 16 ˚A is used in the vasp calculations.Based on geometry relaxations on the DFT level,in the out-of-plane direction, the V-O distance isreduced from 1 .
96 ˚A to 1 .
93 ˚A and the Sr-Sr distancefrom 3 .
92 ˚A to 3 .
52 ˚A. The Brillouin zone is sampledusing a 15 × × k -grid. The energy cutoff of the plane wave basisset is 400 eV, in accordance with the default valuefor the PAW potentials used. Using the proceduredescribed above, the projectors onto the V d -statesare calculated according to (12) by optimizing in anenergy window around the Fermi level given by ε P min = − ε P max = 1 . t g character) and orthonormalizing according to (14)in the same energy window. In the t g subspace, aHubbard-Kanamori interaction with U = 5 . J = 0 .
75 eV (in agreement with [13]) is added; theresulting double counting is estimated in the FLLscheme.[33, 37] The impurity problem is solved usingthe triqs / CTHYB
Quantum Monte Carlo[38] solverat an inverse temperature β = 40 eV − .The one-shot DFT+DMFT calculation results ina nearly complete polarization of the orbitals (seefillings in table 1), which is found to be equal in vasp and wien2k , and which is also in accordance withpublished data.[13] When using the charge feedback inthe fcsc framework, the empty orbitals become partlyrepopulated. This effect happens slightly stronger harge self-consistent many-body corrections using optimized PLOs Table 1.
Filling of the correlated orbitals in DFT+DMFT forone-shot and fully charge self-consistent calculations based on vasp and wien2k . one-shot fcsc d xy d xz , d yz d xy d xz , d yz vasp wien2k in vasp but the agreement of the two calculationsbased on the two different DFT codes is within theexpected difference between the two implementations.This re-population can also be seen in the spectralfunction (figure 5), which is obtained using analyticcontinuation of the local lattice Green’s function usingthe maximum entropy method.[39] Unlike the one-shot calculation (top panel of figure 5), the fcscscheme produces a lower Hubbard band with non-negligible spectral weight also for the degenerate d xz and d yz orbitals (bottom panel of figure 5). Thesame Mott gap (whose value is in good agreementwith experiment [35]) is found starting from both DFTcodes, with the peak positions being basically identical.The difference in peak height is compatible with thedifference of filling between the two methodologies.The calculated spectra are in excellent agreement withresults presented in [13].Full charge self-consistency allows us to calculatereliably the total energy as a function of structuralparameters and to determine the lowest-energystructure in DFT+DMFT. Here, as a proof ofprinciple, we calculate the total energy of thecompound as a function of the distance between Sr-O and V-O planes. We consider deviations ∆ z fromthe DFT-optimized structure. Positive ∆ z means thatthe upper Sr-O plane gets shifted upwards and thelower plane gets shifted downwards, thus increasingthe thickness of the slab in z direction. This changesthe splitting between the d xy , which is close to half-filling, and the degenerate d xz and d yz orbitals, whichare nearly empty. For more negative ∆ z (i.e., thinnerslabs), the d xy orbital approach even more half-filling,while the d xz and d yz orbitals are progressively emptied(see figure 6, bottom). Additionally, the bandwidthdecreases, enhancing thus the correlations and the sizeof the Mott gap (not shown here). The main resultof this total-energy calculation is that the minimumcalculated with DFT+DMFT is shifted significantlytowards lower ∆ z (Figure 6, top), indicating that thestructure with the lowest energy has a smaller slabwidth than obtained by DFT.The trend in the structural change can be roughlyexplained in terms of an additional energy gain byremoving the degeneracy between the in-plane d xy andout-of-plane d xz , d yz orbitals for smaller values of ∆ z . A () one-shot xz,yz VASPxy VASP xz,yz Wien2kxy Wien2k / eV A () fcsc Figure 5.
DFT+DMFT spectral function of the single layerof SrVO in a one-shot (top) and fully charge self-consistentcalculation (bottom). The calculations have been performedusing triqs / dfttools , once with wien2k (compare [13]) andonce with vasp as underlying DFT code. The resultingimaginary-time Green’s function was analytically continuedusing the Maximum Entropy Method.[39] E t o t ( e V ) DFTDMFT z (Å) f illi n g o f d xz , d yz Figure 6.
Top: Total-energy change of the single layer of SrVO when moving the upper and lower Sr-O-plane symmetrically by∆ z (from the DFT-optimized structure for ∆ z = 0) in DFT andin fully charge self-consistent DFT+DMFT. For convenience, theenergy for ∆ z = 0 is shifted to 0 in each case. Bottom: Fillingof the degenerate d xz and d yz orbitals per spin channel as afunction of ∆ z . The total filling of the impurity is 1 electron.The lines in both plots are guides to the eye. The values anderror bars of the DFT+DMFT calculations in both plots wereobtained by calculating the quantity for the four last iterationsand then determining the mean and the standard deviation. Formany data points the error bars are smaller than the markers.The total energy in DFT was converged to 10 − eV. harge self-consistent many-body corrections using optimized PLOs U − J . Once the apical oxygen is movedsufficiently close to the V ion, the anti-bonding orbitals d xz , d yz are pushed up in energy and only d xy remainsoccupied (half-filled). This removes the interorbitalenergy cost, lowering thus the total energy.
5. Conclusion
We have presented a charge self-consistent implemen-tation to combine DFT with many-body techniques inthe vasp package, based on optimized projector local-ized orbitals (PLO) in the PAW framework. The im-plemented optimization, seeking the partial wave withthe largest overlap with the relevant correlated sub-space, is crucial for concise projections and leads to astraight-forward connection between delocalized Kohn-Sham states and localized basis functions. As usual, inthe localized subspace, Hubbard-like Hamiltonians canbe used straightforwardly. In contrast to a maximally-localized Wannier projection, the projector formalismis very simple, easy to implement, preserves symme-try, and does not require any special precautions forstrongly entangled bands. Therefore, the projector for-malism is also well suited for the simulation of correla-tion effects in supercells with a large number of bands.We have exemplified the benefits of using optimizedprojectors for the case of NiO.We have benchmarked our fcsc implementationfor two cases. First, we have compared a standard vasp
DFT+U calculation with a charge self-consistentmean-field treatment of the DFT+DMFT Hamiltonianfor the case of NiO. We find only small deviations forthe DOS, which we relate to different projections usedin the standard vasp
DFT+U and the DFT+DMFTscheme. For NiO an important finding is that thedouble-counting problem is alleviated by the chargeself-consistency. With charge self-consistency theinfluence of the double-counting parameter on the bandgap is reduced by about a factor 2 compared to one-shot calculations.Second, we simulated a SrVO monolayer andfound strong orbital polarization, which is decreasedin charge self-consistency. This agrees with previ-ously published results using FLAPW+DMFT. Addi-tionally, as a proof of concept, we calculated the totalenergies from DFT+DMFT and found that correlationeffects lead to structural changes in SrVO monolayers,reducing the apical oxygen height in the single layer.The presented projector and fcsc scheme can beused to interface basically any many-body methodwith the vasp package. It offers a robust andconcise interface for materials studies as well as futuredevelopments of tools for strongly-correlated electron systems. Acknowledgments
We thank Olivier Parcollet and Michel Ferrero fordiscussions about the triqs interface. Funding bythe Austrian Science Fund (FWF) within ProjectNo. Y746 and within the F41 (SFB ViCoM)is gratefully acknowledged. O. P. is gratefulto A. Georges for discussions on the details ofcharge self-consistency. O. P. acknowledges supportfrom the Swiss National Science Foundation NCCRMARVEL and computing resources provided bythe Swiss National Supercomputing Centre (CSCS)under projects s575 and mr17. The supportfrom the Austrian federal government (in particularfrom Bundesministerium f¨ur Verkehr, Innovation undTechnologie and Bundesministerium f¨ur Wirtschaft,Familie und Jugend) represented by ¨OsterreichischeForschungsf¨orderungsgesellschaft mbH and the Styrianand the Tyrolean provincial government, representedby Steirische Wirtschaftsf¨orderungsgesellschaft mbHand Standortagentur Tirol, within the framework ofthe COMET Funding Programme is also gratefullyacknowledged. The work at the University of Bremenhas been supported by the Deutsche Forschungs-gemeinschaft (DFG) via FOR 1346 as well as theZentrale Forschungsf¨orderung of the University ofBremen.
References [1] Ceperley D M and Alder B J 1980
Phys. Rev. Lett. (7) 566–569 URL https://link.aps.org/doi/10.1103/PhysRevLett.45.566 [2] Perdew J P, Burke K and Ernzerhof M 1996 Phys. Rev.Lett. (18) 3865–3868[3] G¨orling A 1996 Phys. Rev. A (5) 3912–3915 URL https://link.aps.org/doi/10.1103/PhysRevA.54.3912 [4] Georges A, Kotliar G, Krauth W and Rozenberg M J 1996 Rev. Mod. Phys. (1) 13–125[5] Lichtenstein A I and Katsnelson M I 1998 Phys. Rev. B http://link.aps.org/doi/10.1103/PhysRevB.57.6884 [6] Kotliar G, Savrasov S Y, Haule K, Oudovenko V S,Parcollet O and Marianetti C A 2006 Reviews of ModernPhysics Phys. Rev. B (24) 17953–17979[8] Amadon B, Lechermann F, Georges A, Jollet F, WehlingT O and Lichtenstein A I 2008 Phys. Rev. B Journal of Physics: Condensed Matter http://iopscience.iop.org/0953-8984/23/8/085601 [10] Pourovskii L V, Amadon B, Biermann S and Georges A2007 Physical Review B Phys. Rev. B Phys. Rev. B (5) 054529 URL http://link.aps.org/doi/10.1103/PhysRevB.84.054529 harge self-consistent many-body corrections using optimized PLOs [13] Bhandary S, Assmann E, Aichhorn M and Held K 2016 Phys. Rev. B Comp. Mat. Sci. [15] Grieger D, Piefke C, Peil O E and Lechermann F 2012 Physical Review B https://link.aps.org/doi/10.1103/PhysRevB.86.155121 [16] Amadon B 2012 J. Phys. C Computer Physics Communications [18] Kresse G and Joubert D 1999
Phys. Rev. B (3) 1758–1775[19] Kresse G and Hafner J 1993 Phys. Rev. B (1) 558–561[20] Kresse G and Furthm¨uller J 1996 Phys. Rev. B (16)11169–11186[21] Parcollet O, Ferrero M, Ayral T, Hafermann H, KrivenkoI, Messio L and Seth P 2015 Computer PhysicsCommunications
WIEN2k, An augmented Plane Wave + LocalOrbitals Program for Calculating Crystal Properties (Techn. Universitaet Wien, Austria, ISBN 3-9501031-1-2.)[23] Aichhorn M, Pourovskii L, Seth P, Vildosola V, Zingl M,Peil O E, Deng X, Mravlje J, Kraberger G J, MartinsC, Ferrero M and Parcollet O 2016
Computer PhysicsCommunications
Rev. Mod. Phys. (4) 1419–1475[25] Savrasov S Y and Kotliar G 2004 Physical Review B Physical ReviewLetters https://link.aps.org/doi/10.1103/PhysRevLett.53.2339 [27] Bengone O, Alouani M, Bl¨ochl P and Hugel J 2000 PhysicalReview B https://link.aps.org/doi/10.1103/PhysRevB.62.16392 [28] Anisimov V I, Zaanen J and Andersen O K 1991 PhysicalReview B https://link.aps.org/doi/10.1103/PhysRevB.44.943 [29] Ren X, Leonov I, Keller G, Kollar M, Nekrasov I and Voll-hardt D 2006 Phys. Rev. B (19) 195114 URL https://link.aps.org/doi/10.1103/PhysRevB.74.195114 [30] Karolak M, Ulm G, Wehling T, Mazurenko V, PoteryaevA and Lichtenstein A 2010 Journal of ElectronSpectroscopy and Related Phenomena
11 – 15 ISSN0368-2048 proceedings of International Workshop onStrong Correlations and Angle-Resolved PhotoemissionSpectroscopy 2009 URL [31] Leonov I, Pourovskii L, Georges A and Abrikosov I A 2016
Phys. Rev. B (15) 155135 URL https://link.aps.org/doi/10.1103/PhysRevB.94.155135 [32] Amadon B 2012 Journal of Physics: Condensed Matter http://stacks.iop.org/0953-8984/24/i=7/a=075604 [33] Anisimov V I, Solovyev I V, Korotin M A, Czy˙zyk M T andSawatzky G A 1993 Physical Review B http://link.aps.org/doi/10.1103/PhysRevB.48.16929 [34] Park H, Millis A J and Marianetti C A 2014 Phys. Rev.B (24) 245133 URL https://link.aps.org/doi/10.1103/PhysRevB.89.245133 [35] Yoshimatsu K, Okabe T, Kumigashira H, Okamoto S,Aizaki S, Fujimori A and Oshima M 2010 Phys. Rev.Lett. (14) 147601 URL https://link.aps.org/doi/10.1103/PhysRevLett.104.147601 [36] Zhong Z, Wallerberger M, Tomczak J M, Taranto C,Parragh N, Toschi A, Sangiovanni G and Held K 2015
Physical Review Letters https://link.aps.org/doi/10.1103/PhysRevLett.114.246401 [37] Held K 2007
Adv. Phys. Computer Physics Communications