Effective Coulomb interactions in solids under pressure
Correlated materials are extremely sensitive to external stimuli, such as temperature or pressure. Describing the electronic properties of such systems often requires applying many-body techniques to effective low energy problems in the spirit of the Hubbard model, or extensions thereof. While the effect of pressure on structures and bands has been investigated extensively within density-functional based methods, the pressure dependence of electron-electron interactions has so far received little attention. As a step toward ab initio pressure studies for realistic systems within a setup of maximally localized Wannier functions and the constrained random phase approximation, we examine in this paper the paradigmatic pressure dependence of Coulomb interactions. While compression commonly causes the “extension” of Wannier functions, and thus transfer elements, to grow, we find the – seemingly counter-intuitive – tendency that the bare Coulomb interaction increases under compression as well. We reconcile these behaviors by appealing to a semi-analytical tight-binding model. We moreover argue that, for this model, the requirement of maximal Wannier localization is equivalent to maximizing the Coulomb interaction matrix elements. We then apply the above first principles techniques to fcc hydrogen under pressure. While we find our comprehension of the bare Coulomb interaction confirmed, the induced changes in screening strengths lead to an effective one band model with a Hubbard interaction that is non-monotonous under pressure.
The panoply of structural, orbital and spin degrees of freedom and the joint presence of important electronic Coulomb interactions cause correlated materials  to be the realm neither of band-theory nor of model many-body physics each on their own. Therefore, in recent years, much ingenuity was invested into finding ways to merge the “realism” of the former with the accurate description of correlation effects of the latter.
With the exception of the GW approximation [2, 3, 4] to Hedin’s equations , the electronic structure approach GW+DMFT  that combines GW with the dynamical mean-field theory (DMFT) , and a recent proposal for a self-energy downfolding , this combining commonly amounts to extracting a low-energy one-particle Hamiltonian from density-functional-theory-based methods, such as the local-density approximation (LDA), [8, 9] supplementing it with interaction terms [e.g. of the Hubbard-Hund (U and J) type] and to solving the resulting “realistic model” with a chosen many-body technique.
Hence, in this approach there are two intertwined preparative tasks that generate the many-body problem. The deducing of the low-energy one-particle part can be achieved e.g. by a tight binding fit of relevant bands, the downfolding [10, 11] procedure within e.g. muffin-tin-based methods, [11, 12, 13] or by the generation of (maximally localized [14, 15]) Wannier functions.
The interaction parameters of realistic many-body models, in turn, are often chosen rather empirically than from a solid first principles basis, a fact that has caused many objections in the past. In particular, when tracking properties as a function of an external parameter – pressure, in our case – the evolution of the interaction has mostly been discarded.
This and the quest for going beyond mere qualitative results towards, eventually, the quantitative design of materials, highlights the need for accurate ways to determine all ingredients of realistic models in an ab initio fashion. Nowadays, the most popular methods for the computation of interaction matrix elements are the constrained LDA technique , and the constrained random phase approximation (cRPA) .
A recent, and promising approach is the use of Wannier functions within the cRPA setup , which allows for a deducing of the one-particle and two-particle parts of the Hamiltonian on the same footing. Moreover, working in a localized Wannier type of basis is often a requirement of many-body approaches such as the DMFT . As to the interaction matrix elements, the cRPA technique allows for a precise elimination of the screening channels of the chosen orbital subspace that constitutes the effective model .
While these techniques have already been applied for the setting up of many-body models of some complex materials [19, 20, 21], a basic understanding of the effects of pressure on the Coulomb interactions within a Wannier setup is lacking. This is the aim of the current work.
In a first part, we investigate a semi-analytical tight binding model of a one-band solid in one dimension, track transfer matrix elements, the bare (i.e. unscreened) Coulomb interaction, and the spread of the maximally localized Wannier orbitals as a function of lattice spacing. Being able to access the decomposition of the maximally localized Wannier functions onto the tight-binding basis will allow to understand the surprising finding that under pressure the Coulomb interaction matrix element augments, while, at the same time, transfer elements describing the delocalization of electrons grow as well. As a more realistic example, we, second, apply the fully ab initio approach of the cRPA within maximally localized Wannier functions  to fcc hydrogen, which is found to exhibit the explained generic behavior of the bare Coulomb interaction matrix elements. However, the partially screened Coulomb interaction – the Hubbard of an effective low energy model for the half-filled 1s orbital – actually shows a non-monotonous trend – a consequence of two opposing effects in screening processes.
The method of using the Wannier orbital construction in conjunction with the cRPA technique has been presented in detail in Ref. . While not fully reviewing the approach, we will discuss some issues relevant for the understanding of our results and introduce some notation.
For the one-particle band-structure, a density functional calculation is performed. For the realistic case of fcc hydrogen, we will employ the LDA  in the full-potential linear muffin-tin orbital (FP-LMTO)  realization. For obtaining the random phase approximation (RPA) polarizations we employ the code of Ref.  with the maximally localized Wannier extension of Ref. , and construct an effective system for the isolated 1s orbital. That is, we introduce the sub-Hilbert space , spanned by the 1s Kohn-Sham wave function. Since the aim is to work within a localized basis, the extraction of the low-energy part is done by a construction of Wannier functions  for , as described in Refs. [14, 15, 19]. The corresponding effective interactions are then computed within the constrained RPA  formalism. This amounts to screening the matrix elements of the bare Coulomb interaction , which in the Wannier basis are given by
with a partial RPA polarization
where and are the polarizations of the full and the sub-Hilbert spaces, respectively. The latter , when using the Kohn-Sham orbitals, can be expressed as
i.e. by transitions restricted to the effective sub-system, in our case . Within this notation, the strengths of screening channels are influenced by two effects : the matrix elements [the overlap integrals of wave functions that occur when calculating matrix elements of , in analogy to Eq. (1)], and the energy differences of the Kohn-Sham excitations, , that appear in the denominators. The virtue used for the constraining is the fact that the screening contributions are additive . Indeed, when using the above decomposition , the fully screened interaction (of the GW formalism [2, 3, 4]) can be given in terms of the partially screened interaction for the effective model of the one-band orbital, , by . The Hubbard of the 1s sub-system is given by the on-site Wannier function matrix element of .
The major observation in this context is that the construction of Wannier functions, and thus also the determination of interaction matrix elements are not unique . Indeed a unitary transformation applied to the periodic part of the wave functions, while preserving the Bloch functions, changes the Wannier orbitals. This gauge freedom can be used to choose the Wannier basis that is most suitable for the final purpose. The aim of the low-energy system is to correct for local (Hubbard-type) interaction effects. To this end, there exist at least two proposals on how to choose an optimal Wannier basis set : in the maximally localized Wannier approach [14, 15] the extension (“spread”) of the Wannier functions is minimized : denoting the Wannier states by kets, , this spread can be chosen as 
The minimization of this is of course only one of the possible options to fix the Wannier functions, yet a very natural one, since it can e.g. be shown [24, 14] that, in one dimension, the resulting Wannier functions are eigenfunctions of the subset projected position operator, an intuitive criterion for real space localization 111In higher dimensions, the non-commutation of the components of the position operator impedes this property [24, 14]..
In the second scheme, the screened local Coulomb interaction matrix element – the Hubbard , as given by the on-site component of from above – is maximized [25, 19] to determine the Wannier functions. For the case of SrVO, it has been shown  that both approaches yield very similar results. Indeed, appealing to e.g. the equation of the bare Coulomb interaction matrix element in the Wannier basis, Eq. (1), it seems plausible that a greater localization of Wannier functions (smaller ) results in an increased interaction matrix element . For the simple model in one dimension that we discuss in the following, we in fact motivate the equivalence of the maximally localized Wannier functions and the basis in which the (bare) Coulomb interaction matrix element is maximal.
However, we stress that the intuitive correspondence between a stronger localization of Wannier functions, in the sense of the spread , and a larger interaction matrix element in this basis does not hold in general. Here, one has to distinguish between the changes in the Wannier localization that are induced by a modification of the Wannier gauge from those caused by modifications of the lattice, caused e.g. by pressure. The fact that for the discussed model, both methods to fix the Wannier gauge are equivalent, states that for a given pressure the maximally localized Wannier functions yield the maximal Coulomb interaction matrix element. Yet, as we shall see, a pressure induced increase in the converged minimal Wannier spread is actually quite naturally concomitant with a greater bare interaction.
Iii Results and discussion
External pressure, or structural changes in general, provides an impetus to alter not only the one-particle band-structure of a material, but also the Wannier functions. Therefore, when investigating the pressure dependence of effective, i.e. screened, interaction matrix elements, one has to distinguish between influences of the former, which enter via a modification of screening channels, and of the latter that not only affects the polarization, but, on a more fundamental level, modifies already the bare Coulomb interaction matrix elements. The fact that Wannier functions of a solid are not eigenfunctions of the system, may result in a nonstraightforward evolution when parameters such as external pressure change.
iii.1 One-band tight-binding model in one dimension
As a first model system, we investigate a tight-binding parametrization of a one-band solid in one dimension. In that case the maximally localized Wannier function is naturally given by the Fourier transform of the Bloch functions when inversion symmetry is verified [24, 14]. Moreover, since no higher energy orbitals are present, no partial screening can occur [in Eq. (2) : ] and the Hubbard equals the on-site matrix element of the bare Coulomb interaction, .
iii.1.1 Bloch function in tight binding
As building blocks of the tight binding basis functions we opt for a hydrogen-like 1s orbital in one dimension
with the Bohr radius . This orbital is a solution to the time-independent Schrödinger equation with a single binding delta potential : It represents an eigenstate of the “atomic” Hamilton operator
with the eigenvalue . While the tight-binding approach with this orbital can in principle be used to approximate the Bloch eigenfunction for a Hamiltonian with any potential, its use is obviously most justified for a Kronig-Penney type of model [26, 27] with a Dirac-comb potential.
The Bloch function is written as :
Here, the factors assure the orthonormality of the Bloch function , and is determined to be :
where, with the lattice constant , and the integer , denotes the distance to the th neighboring site. Further, denotes the overlap integral between the atomic function at the origin and its th neighbor, and is given by :
iii.1.2 Wannier function
where we defined
or, with , and the reference ,
Demanding , the normalization constant becomes :
The quantity is real for all for symmetry reasons. Figure 1 shows its behavior for different lattice constants : in the limit of large atomic separation (), the overlaps are negligible, and the Wannier function , Eq. (11), will equal the atomic orbital , Eq. (5). Thus the distribution picks up a single mode of the array, for the representative site “0”. When pressure is applied, and the lattice constant shrinks, finite overlaps of the (non-orthogonal) hydrogen orbitals entails contributions from neighboring sites to mix in, and the distribution broadens (see Figure 1). Since 222This can be seen (hand-wavingly) as follows. Consider only the first term in the sum. Then is defined for , which is verified for . For odd , the relevant (=even) part of the periodic nominator, is always negative at the boundaries . On the other hand, is always positive (), and becomes largest at the boundary. Thus the momentum integral that would vanish without the denominator, becomes negative, since the modulation of the denominator gives the negative values a higher weight. Since the argument will remain valid when considering the complete sum over . , and , the normalization of the Wannier function, Eq. (14), causes the coefficient of the atomic orbital at the origin to become larger than 1 : . This results in a greater probability density around the site origins, , when pressure is applied.
The corresponding Wannier functions of the above cases are shown in Figure 2. As anticipated from the above discussion, more weight is accumulated at the origin : a harbinger for a larger on-site Coulomb interaction. On the other hand, contrary to the atomic limit, the tails of the Wannier functions extend over several lattice constants, before the exponential decay sets in. This behavior at larger distance points to an increase in the Wannier spread and the growing of transfer integrals.
iii.1.3 Wannier spread
Since by symmetry, , the spread of the Wannier function [ as given in Eq. (4)] reduces to :
While the first term is always positive, the second will be negative for odd and positive for even neighbors (see the form of ). Given the decay behavior of , the second term will thus be negative in total, yet small. Indeed the major contribution to the spread comes from the first term, making it plausible that the spread, as defined by , increases with a more widely distributed , as shown in Figure 3.
Although helpful for the understanding of the current model study, the spread is not too good a quantity for gaining quantitative insights from pressure studies of realistic systems, as we will discuss later.
iii.1.4 Transfer integral
For the transfer integral , we need to explicitly specify the non-interacting Hamiltonian, , and, for simplicity, we shall choose a Kronig-Penny-type model [26, 27] with an ionic Dirac-comb potential, i.e.
Then the nearest neighbor hopping can be expressed as :
As can be inferred from the dependence of the coefficients , Eq. (LABEL:An), and the overlap , Eq. (9), the pressure induced delocalization increases the transfer integral, as expected. Figure 4 displays the hopping as a function of the lattice constant.
iii.1.5 Electron-electron interaction
All of the above were concerned with the one-particle picture, i.e. the non-interacting Hamiltonian . The on-site electron-electron interaction matrix element , that governs the two-particle term in the final many-body model, reads in the Wannier basis :
when choosing a hardcore (h.c.) or Coulombic (C.) interaction :
As shown in Figure 5 and anticipated before, the sharpening of the Wannier function at the origin causes greater interaction matrix elements, when pressure is applied, irrespective of the chosen electron-electron interaction.
As previously stated, the on-site electron-electron interaction from above equals the Hubbard , since higher energy orbitals, and thus screening effects, are absent by construction.
iii.1.6 Maximally localized Wannier functions versus maximal Hubbard interaction
In Sec. I we mentioned that another technique to choose the Wannier function gauge is given by the request to maximize the static local, partially screened Coulomb interaction [25, 19] – the Hubbard .
While not actually performing this approach for our model, we give evidence that in this simple case, both techniques are equivalent. As said before, the Hubbard equals the bare Coulomb interaction since we consider only a single band, so, contrary to the general case , the argument does not involve any screening related effects.
The quest is thus to find a Wannier gauge, meaning an additional factor in Eq. (10), that yields the greatest interaction element as given by Eq. (III.1.5). The choice of gauge can be absorbed into the distribution , and we define
Having seen the correspondence between the value of and the interaction strength, one might endeavor to solve the functional derivative for . Yet the dependence of the normalization constant, , given in Eq. (14), makes this a rather tedious task analytically. Instead, we shall argue for a specific example for the choice of Wannier gauge, and make some general comments. Consider a gauge field that is linear in momentum, [see also Ref. ]. Figure 6 displays the Wannier function (upper panel) and the resulting on-site interaction (bottom panel) for the hardcore case (see Eq. (19)), for different gauge parameters , and for a fixed lattice spacing . With the inversion symmetry of the Wannier function with respect to the site is lost. This was the requirement for maximally localized Wannier functions in one dimension [24, 14], and as seen in Figure 6, for the tail of the Wannier function augments, and the Coulomb interaction decreases. While is a mere lattice translation, amounts to a case, where the Wannier function of site has equal positive weight at and , corresponding to . As discussed above, an increase in is caused by the mixing in of negative components to from neighboring sites, leading to a smaller normalization factor in Eq. (14). Owing to the symmetry, in the current example, , negative contributions will come only from the sites for and for . Yet the overall gain in renormalization is distributed over the two equivalent positions , and for which . As a result of this shifting of weight to the site , and thus the on-site Coulomb interaction are much lower than in the case with inversion symmetry around that site. As seen in the bottom panel of Figure 6, it is indeed the maximally localized Wannier functions () that yield the greatest possible Coulomb interaction matrix element for our simple model.
For gauge fields depending on the momentum to the power , the argument is geometrically less obvious, but still true as verified numerically. Indeed, only the integrand in the coefficient of Eq. (LABEL:An2) is always positive for as a function of . Therefore any modulation in , with , will decrease the corresponding integral to a greater extent than for , in which case the integrand changes sign with already for . As a consequence, the decrease in the contribution to the overall normalization will be greater on a relative scale than for , and thus decreases with any non-constant .
This can be taken as a further indication that also for realistic systems, the maximally localized Wannier functions and the maximal Hubbard approach are generally giving the same results.
iii.2 fcc hydrogen under pressure
Towards a more realistic application of the gained insight, we apply the fully ab initio approach  of the cRPA within maximally localized Wannier functions to the “simplest” realistic system, namely solid hydrogen.
While at low pressure, solid hydrogen forms a molecular crystal. It was conjectured, already in the 1930s, that at high pressure hydrogen should become an isotropic metal . Here, however, we shall not be interested with the precise phase diagram of solid hydrogen. For the sake of simplicity, we assume throughout the discussion a face-centered cubic crystal structure with variable lattice constant. Moreover, we are well aware that for the current case the problem of self-interaction  within the LDA formalism is a particularly severe issue. However, here, we are not concerned with the accuracy of the LDA band-structure but with compression induced trends in a 3d multi-orbital setup.
For the one-particle band-structure, we employ in this work the LDA in the full-potential FP-LMTO  realization. In the LMTO basis, we include orbitals up to the 4f, using local orbitals  for multiple orbitals per l-channel, and use a Brillouin zone discretization with up to points. As described in Ref. , maximally localized Wannier functions are then constructed for the hydrogen 1s band, which is entirely isolated from all other Kohn-Sham excitations. In other words, the effective model consists of a single half-filled orbital. We stress that while the Wannier functions of different sites are orthogonal by constructions, the LMTO basis functions – in analogy to the tight-binding parametrization in the preceding section – are not. As a consequence, the same prototypical response to pressure as discussed above can be expected in the current case.
Figure 7 shows the LDA band-structure for the extremal lattice constants that we consider. As expected, under growing compression, the dispersions increase and unoccupied bands are shifted upwards333The fact that the energy difference between the 1s and 2s,2p bands at a.u. is already smaller than the atomic ionisation energy is owing to the aforementioned spurious self-interaction in the LDA. .
In line with this is the behavior of the hopping – the nearest neighbor transfer matrix element in basis of the maximally localized Wannier function for the subspace of the lowest Kohn-Sham excitation. As shown in Figure 9, it augments with decreasing lattice constant , accounting for the greater delocalization. In Figure 8 is shown the maximally localized Wannier function (it is real, cf. Ref. ) of the hydrogen atom at the origin, as a function of the (scaled) distance towards the nearest neighbor at e.g. . When the lattice constant shrinks, clearly witnessed is both, a growing weight at the neighboring position (), as well as an increase in the Wannier function at the origin (). This is in complete analogy to what was shown in Figure 2 for the model considered above444We note that the issue of self-interaction inherent in the LDA treatment impedes the Wannier function to converge towards the well-known 1s atomic eigenfunction . This is why the values of for large separation are below teh indicated . .
It is thus expected that the on-site matrix element of the bare Coulomb interaction grows under compression as before. And, indeed this is the case, as can be inferred from Figure 9 in the second panel from the bottom. Yet, does that entail for the Wannier spread the same behavior as witnessed in the one dimensional model ? In three dimensions, the gain in spread by effects of hybridizations with neighboring sites might turn out less prominent, since the region occupied by nearest neighbor atoms is relatively small. Hence, the angular integral in , Eq. (4), even for the radius of the distance to the 12 neighboring atoms (fcc) will run very much over a sphere on which the Wannier function is mostly zero. And, indeed, as displayed in the second panel from the top in Figure 9, the Wannier spread actually decreases under compression555We note that the true value of the spread at infinite separation should be – the atomic value. Again this is due to the inadequateness of the LDA. .
A direct interpretation of the spread in pressure studies has thus to be taken with caution. First, depending on the crystal-structure (dimension, number and position of neighboring atoms) the spread does not necessarily reflect the increased delocalization of charge carrier but can be dominated by the more isotropic concentration of weight at the origin (see Figure 8). Second, as seen in Figure 9 (and already in the seminal work, Ref. , Table I), the momentum space convergence of is poor. Also, a change in lattice constant upon compression, changes the accuracy of the spread function when the number of k-points is kept constant. While the spread is of course the entity that is minimized in order to obtain the maximally localized Wannier functions for a given pressure, the spread itself is not a reliable measure for trends in the Wannier functions upon compression. Instead, one should either plot the functions, or resort to the bare Coulomb interaction matrix elements, which – owing to the different powers of the position operator – converges well for a moderate k-point sampling.
In the current case, while still constructing an effective one-band model, the initial system contains higher energy orbitals. As a result, and contrary to the simple tight-binding model, there is a non-zero partial polarization, of Eq. (2), that screens the bare Coulomb interaction . The on-site part of the screened interaction – the Hubbard – within cRPA is displayed in the lowest panel of Figure 9.
Interestingly, this quantity, in contrast to the bare interaction, shows a non-monotonic behavior, that reflects the struggle of two opposing effects in the polarization . As can be inferred from Eq. (2), and the equation for , Eq. (II), the changes in the polarization originate from modifications of transitions from the Wannier “1s” into orbitals at higher energy. These can be altered by two ingredients (see again Eq. (II)) : the transition matrix elements, and the (Kohn-Sham) transition energies , i.e. the band-structure. The effect of pressure will be different for these two mechanisms. Indeed, the increasing compression pushes the bands further apart, as seen in Figure 7, diminishing the polarization. The increased overlaps/hybridizations of orbitals, on the other hand, tend to make the polarization grow. In order to separate the influence of the two contributions, we computed the partially screened Coulomb interaction, , as a function of lattice spacing, albeit while keeping the band-structure fixed at the values obtained for the highest pressure, . Thus the corresponding screened interaction will contain only the effect of changes in the Bloch functions. As indicated by the symbols in the figure of the Hubbard U, the latter decreases under compression, proving the above conjectured opposition in trend to the influence of the band-structure.
This effect is surely very material specific. In systems with more electrons, where pressure will e.g. cause the enhancement of bonding/anti-bonding splittings, it can be expected that the changes in the band-structure often prevail such as to reduce the polarization under compression.
Iv Conclusions and Perspectives
In conclusion, we have studied the influence of external pressure onto the construction of effective low energy many-body systems. Using maximally localized Wannier functions for a one-band tight binding model, we rationalized the counter-intuitive, yet prototypical behavior of the bare Coulomb interaction, namely that its matrix elements augment upon compression, as a consequence of the delocalization of the Wannier functions. This we understood to be caused by increased admixtures of non-orthogonal nearest neighbor tight-binding functions when the lattice spacing shrinks. As a more realistic system, we investigated fcc hydrogen under pressure, and constructed an effective model for the half-filled 1s orbital using maximally localized Wannier functions. For the transfer integrals and the bare Coulomb interaction, we witnessed the same tendencies as in the model case. Yet, the Hubbard , calculated as the on-site screened interaction within the constraint RPA technique, exhibited a non-monotonous trend under compression. This we traced back to a struggle between two opposing effects in the strength of screening. All these highlight the intricacy of mechanisms that influence effective models, and emphasizes the need for reliable ab initio techniques for their construction.
This work was in part supported by the G-COE program of MEXT(G-03) and Grant-in-Aid for Scientific Research from MEXT Japan (Grants No. 19019013 and No. 19051016).
-  M. Imada, A. Fujimori, and Y. Tokura. Metal-insulator transitions. Rev. Mod. Phys., 70(4):1039–1263, Oct 1998.
-  L. Hedin. New method for calculating the one-particle green’s function with application to the electron-gas problem. Phys. Rev., 139(3A):A796–A823, Aug 1965.
-  F Aryasetiawan and O Gunnarsson. The gw method. Rep. Prog. Phys., 61(3):237–312, 1998.
-  Giovanni Onida, Lucia Reining, and Angel Rubio. Electronic excitations: density-functional versus many-body green’s-function approaches. Rev. Mod. Phys., 74(2):601–659, Jun 2002.
-  S. Biermann, F. Aryasetiawan, and A. Georges. First-principles approach to the electronic structure of strongly correlated systems: Combining the approximation and dynamical mean-field theory. Phys. Rev. Lett., 90(8):086402, Feb 2003.
-  A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg. Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions. Rev. Mod. Phys., 68(1):13, Jan 1996.
-  F. Aryasetiawan, J. M. Tomczak, T. Miyake, and R. Sakuma. Downfolded self-energy of many-electron systems. Phys. Rev. Lett., 102(17):176402, May 2009. preprint: arXiv:0806.3373.
-  W. Kohn. Nobel lecture: Electronic structure of matter-wave functions and density functionals. Rev. Mod. Phys., 71(5):1253–1266, Oct 1999.
-  R. O. Jones and O. Gunnarsson. The density functional formalism, its applications and prospects. Rev. Mod. Phys., 61(3):689–746, Jul 1989.
-  Per-Olov Löwdin. A note on the quantum-mechanical perturbation theory. J. Chem. Phys., 19(11):1396, Nov 1951.
-  Walter R. L. Lambrecht and O. K. Andersen. Minimal basis sets in the linear muffin-tin orbital method: Application to the diamond-structure crystals c, si, and ge. Phys. Rev. B, 34(4):2439–2449, Aug 1986.
-  O. K. Andersen. Linear methods in band theory. Phys. Rev. B, 12(8):3060–3083, Oct 1975.
-  O. K. Andersen and T. Saha-Dasgupta. Muffin-tin orbitals of arbitrary order. Phys. Rev. B, 62(24):R16219–R16222, Dec 2000.
-  Nicola Marzari and David Vanderbilt. Maximally localized generalized wannier functions for composite energy bands. Phys. Rev. B, 56(20):12847–12865, Nov 1997.
-  Ivo Souza, Nicola Marzari, and David Vanderbilt. Maximally localized wannier functions for entangled energy bands. Phys. Rev. B, 65(3):035109, Dec 2001.
-  Gregory H. Wannier. Dynamics of band electrons in electric and magnetic fields. Rev. Mod. Phys., 34(4):645–655, Oct 1962.
-  V. I. Anisimov and O. Gunnarsson. Density-functional calculation of effective coulomb interactions in metals. Phys. Rev. B, 43(10):7570–7574, Apr 1991.
-  F. Aryasetiawan, M. Imada, A. Georges, G. Kotliar, S. Biermann, and A. I. Lichtenstein. Frequency-dependent local interactions and low-energy effective models from electronic structure calculations. Phys. Rev. B, 70(19):195104, Nov 2004.
-  Takashi Miyake and F. Aryasetiawan. Screened coulomb interaction in the maximally localized wannier basis. Phys. Rev. B, 77(8):085122, 2008.
-  Kazuma Nakamura, Ryotaro Arita, and Masatoshi Imada. Ab initio derivation of low-energy model for iron-based superconductors lafeaso and lafepo. J. Phys. Soc. Jpn., 77:093711, 2008.
-  T. Miyake, L. Pourovskii, V. Vildosola, S. Biermann, and A. Georges. d- and f-orbital correlations in the refeaso compounds. J. Phys. Soc. Jpn., 77SC(Supplement C):99–102, 2008.
-  M. Methfessel, Mark van Schilfgaarde, and R.I. Casali. A full-potential lmto method based on smooth hankel functions. in Electronic Structure and Physical Properties of Solids: The Uses of the LMTO Method, Lecture Notes in Physics. H. Dreysse, ed., 535:114–147, 2000.
-  Takao Kotani, Mark van Schilfgaarde, and Sergey V. Faleev. Quasiparticle self-consistent gw method: A basis for the independent-particle approximation. Phys. Rev. B, 76(16):165106, 2007.
-  S. Kivelson. Wannier functions in one-dimensional disordered systems: Application to fractionally charged solitons. Phys. Rev. B, 26(8):4269–4277, Oct 1982.
-  Clyde Edmiston and Klaus Ruedenberg. Localized atomic and molecular orbitals. Rev. Mod. Phys., 35(3):457–464, Jul 1963.
-  R. de L. Kronig and W. G. Penney. Quantum mechanics of electrons in crystal lattices. Proc. R. Soc. London Ser. A, 130(814):499–513, 1931.
-  F. B. Pedersen, G. T. Einevoll, and P. C. Hemmer. Wannier functions for the kronig-penney model. Phys. Rev. B, 44(11):5470–5475, Sep 1991.
-  E. Wigner and H. B. Huntington. On the possibility of a metallic modification of hydrogen. The Journal of Chemical Physics, 3(12):764–770, 1935.
-  J. P. Perdew and Alex Zunger. Self-interaction correction to density-functional approximations for many-electron systems. Phys. Rev. B, 23(10):5048–5079, May 1981.