Band selection and disentanglement using maximally-localized Wannier functions: the cases of Co impurities in bulk copper and the Cu (111) surface

We have adapted the maximally-localized Wannier function approach of [I. Souza, N. Marzari and D. Vanderbilt, Phys. Rev. B 65, 035109 (2002)] to the density functional theory based Siesta method [J. M. Soler et al., J. Phys.: Cond. Mat. 14, 2745 (2002)] and applied it to the study of Co substitutional impurities in bulk copper as well as to the Cu (111) surface. In the Co impurity case, we have reduced the problem to the Co d-electrons and the Cu sp-band, permitting us to obtain an Anderson-like Hamiltonian from well defined density functional parameters in a fully orthonormal basis set. In order to test the quality of the Wannier approach to surfaces, we have studied the electronic structure of the Cu (111) surface by again transforming the density functional problem into the Wannier representation. An excellent description of the Shockley surface state is attained, permitting us to be confident in the application of this method to future studies of magnetic adsorbates in the presence of an extended surface state.


Introduction
The advent of density functional theory (DFT) in electronic structure calculations has revolutionized the fields of quantum chemistry and, more generally, of condensed matter physics [1]. However, due to necessary approximations of the unknown functionals (typically, the local density approximation or LDA, and semilocal approximations as the generalized gradient approximation or GGA), important limitations prevail in the description of many systems ranging from insulating materials to impurities in a metal host. The general approach to these problems is to go beyond DFT by proposing a simplified Hamiltonian that can be solved, at the expense of losing the parameter-free advantage of DFT. Recently, the situation is changing and new methods either stemming from DFT or making use of DFT for initial input are emerging [2,3]. One such example is the LDA+U approach [4,5] where the exchange-correlation potential for the local electron gas is complemented by missing strong localized correlations. The success of this approach has been considerable in explaining the opening of a band gap in systems with localized electronic states. Yet, the choice of the Coulomb intra-atomic energy U is somewhat arbitrary, depending on the choice of basis sets or other descriptors of the treated system [6].
Nakamura and co-workers [7] have developed a constrained LDA approach based on maximally localized Wannier functions (MLWF) [8]. The MLWF replace the linear muffin-tin orbitals (LMTO) that were initially used as a natural basis set to define U [4,9,10]. As Nakamura et al. emphasize [7], the properties of MLWF are particularly appealing for reducing the complicated DFT problem to a simplified Hamiltonian, where especial physics can be explored such as the localized electron correlations mentioned above. Wannier functions have been thoroughly explored in the work by Marzari and Vanderbilt [11], and an algorithmic approach to obtaining them from a DFT calculation is available. More recently, Souza, Marzari and Vanderbilt [8] have extended the approach to disentangle electronic bands, creating a compact local basis set that accurately reproduces the DFT electronic structure in a given energy window. This approach is independent of the actual implementation of the DFT calculation, yielding a natural way of describing an extended basis set calculation in terms of localized functions.
In the present work, we have interfaced the approach by Souza et al. [8] as implemented in Wannier90 [12] to the Siesta code [13]. Contrary to the case of Nakamura et al. [7], Siesta is an atomic basis set code, and it would seem natural to use the atomic orbitals to define U and associated model Hamiltonians. However, besides their optimized spread, MLWF have two important features: (i) they are a naturally orthogonal basis set, rendering tight-binding like approaches easy to use (ii) the extraordinary accuracy of MLWF in a defined energy window permits to disentangle electronic bands [8] and to have a simple tight-binding approach with DFT accuracy within the chosen energy window. Hence, our present implementation of MLWF permits us to translate the complex atomic orbital method into a simple orthogonal tight-binding approach that can be easily cast into an Anderson-Hamiltonian model [14] or a Hubbard one [15].
The manuscript is organized in four sections. The methods section deals with a small overview of the work of reference [8], and the technical details of the actual numerical implementation. The approach is then applied to two different systems: (i) a substitutional Co impurity in bulk Cu and (ii) the Cu (111) surface. Both systems are of uttermost interest. The case of Co in Cu is a classical example where strong correlations are in play leading to very large Kondo temperatures [16,17], hence it is interesting to understand the magnetic properties leading to the Kondo correlations in this system [18]. Besides, Co in Cu is one of the model giant magnetoresistance (GMR) systems [19,20,21]. The concentration of Co in Cu is an important parameter in its magnetic properties [21]. Indeed, when the Co to Cu atomic ratio is above 1 to 4, the system becomes a ferromagnet. Below this density, Fan et al. [21] experimentally find that the Co atoms are distant enough to show paramagnetism in agreement with the paradigmatic analysis by Goodenough [22]. We will concentrate in this paramagnetic phase, with two different densities, for 2 × 2 × 2 and 4 × 4 × 4 cubic supercells.
Our calculations are a first step in the study of the complex electronic structure of the Co-impurity problem that presents a high (∼ 500 K) Kondo temperature [17]. Hence, the evolution of the electronic structure as the Co concentration is reduced, permits us to study the effect of Co concentration and the use of MWLF gives us access to an Anderson-like Hamiltonian.
The study of the Cu(111) surface is also very interesting given the existence of the LL ′ gap (Γ in the surface Brillouin zone) and the associated Shockley state. Despite the locality of the MLWF, the calculations succeed in accounting for the surface Shockley state and in yielding the correct electronic structure about the Fermi energy.
This work shows that MLWF are accurate enough to study the magnetic properties of Co on Cu(111) and are a first step towards the study of the electronic structure [23] on surface systems as well as the Kondo effect of magnetic adsorbates that has recently received much experimental [24,25] and theoretical [26] attention.

Wannier functions and Wannier90
Wannier functions, w(r − R), are formally identical to Fourier coefficients of Bloch waves, ψ k (r), given by Here R denote Bravais vectors and Ω ⋆ is the volume of the Brillouine zone (BZ). This definition suffers from the indeterminacy of the phases of Bloch functions, ψ k (r). Furthermore, if a group of bands, {n} is considered, with corresponding Bloch functions, ψ kn (r), it is ambiguous to talk about an individual band and the above definition can be generalized to where U mn (k) is an arbitrary unitary matrix. The unitarity guarantees that the Wannier functions will be orthogonal. The arbitrariness of U mn (k) allows for tuning the phases of Bloch functions in the integral as well as the admixture of functions pertaining to different bands. Thus, there is a whole class of Wannier functions for the given band structure. Marzari and Vanderbilt devised a variational scheme that determines the Wannier functions with minimum total spread, Ω defined as where we have used the notation · m = w m | · |w m . For a given set of Bloch functions, the total spread, Ω, is a functional of the unitary matrices U mn (k). The Wannier functions so obtained are called maximally-localized Wannier functions (MLWF). The locality and the orthogonality of Wannier functions permit us to have a straightforward compact tight-binding representation of the Hamiltonian in systems where the group of bands of interest is separated from the rest of the electronic structure by a band gap. In metals, the absence of band gaps renders the separation of states more complicated. Souza et al. have devised a disentanglement procedure [8] by focusing on a certain energy window, hereinafter called outer energy window, and by selecting certain bands in it. Hence, the disentanglement procedure necessarily reduces the number of states inside the outer energy window, while exactly reproducing the electronic bands in a certain energy interval called inner energy window. The selection of bands proceeds via trial orbitals g n (r) that allow to define the character of the Wannier subspace of interest.
The numerical calculations of MLWF have been done with the Wannier90 [12] code. This package can be used as a post-processing tool with most first-principles codes. In practical simulations, the Bloch states, |ψ kn , are computed in a mesh of uniformly spaced k-points within the first-Brillouin zone. The basis set used to represent the wave functions changes from one electronic structure code to another one; however, the MLWF algorithm requires an input that is essentially basis-independent. In particular, the main ingredients are (i) the overlap matrix between the cell-periodic parts of the wavefunctions at neighboring k-points (ii) the Bloch energies ε kn on the regular grid of k-points and (iii) the coefficients of the above trial orbitals, g n (r), in the Bloch basis The trial functions we employ are of the form g n (r) = R n (r)Θ lmr (φ, θ), r = r, φ, θ, a product of nodeless hydrogenic-like radial part R n (r) and a real spherical harmonic Θ lmr (φ, θ) with angular momentum l and projection m r . From this, the Wannier90 code computes the final unitary transformation matrices U mn (k). With the aid of the matrices U mn (k) the Hamiltonian can be expressed in the Wannier basis and diagonalized at any k-point. These Wannier-based energy bands are further compared to the initial ab-initio bands in order to test the quality of the newly obtained Wannier basis set.
It is very useful to work with density of states projected onto a certain Wannier function of spin σ, w mσ . This projected density of states (PDOS) [27] is the spectral function, ρ mσ (ω), given by This Wannier PDOS or spectral function shows how the Wannier character is distributed in energy over the electronic band structure ‡. From (6) it is straightforward to define the Wannier function occupation at zero temperature as where µ is the Fermi level.

Pseudo-atomic orbital DFT calculations
Ab-initio DFT calculation presented in this work are based on strictly localized [28] numerical pseudo-atomic orbitals (PAO) [29] that are solutions to the atomic Kohn and Sham equation with norm-conserving pseudopotentials. In particular, the calculations are done using the Siesta [13] package. The matrix elements between basis functions are calculated by real-space integration, and the Hamiltonian eigenstates are labeled using the Bloch theorem, because periodic boundary conditions are imposed. The Bloch functions can be expanded into the PAO basis as follows where we assume that the unit cell contains the centers r µ of basis functions ϕ µ (r − r µ ) which are then repeated periodically to every other cell by the Bravais lattice vector R.
Finally, the complex numbers c µn (k) are expansion coefficients. The Brillouin zone is sampled uniformly. For the exchange and correlation potential entering the Kohn and Sham equation, we use the PBE generalized gradient approximation [30]. The chosen atomic basis set is an optimized double-ζ plus polarization for the valence states of Co and Cu, amounting to 15 basis functions per atom. All the parameters that define the shape and the range of the basis functions were obtained by a variational optimization of the enthalpy (energy plus a penalty for orbital volume increase) with a pressure of P = 0.1 GPa, following the recipe given in reference [31]. The substitutional Co impurity basis set was optimized in the Cu host, and the Cu basis set corresponds to the optimal one for bulk Cu. The question of optimal basis sets for surfaces is more intricate and has been discussed in [32]. We use the basis set obtained in reference [32], which is based on the same enthalpy minimization with an especial focus on the vacuum extension of the surface state density.

Implementation of Wannier functions
The implementation of maximally localized Wannier functions in Siesta consists in evaluation of (4) and (5). Expanding the Bloch functions according to (8), we obtain where and This reduces the computation of (4) and (5) to the calculation of a few matrix elements of localized functions, equations (11) and (12). The first integral is computed on the real space grid while the second integral makes use of the analytic angular dependence of the integrand in the same way as is done for the calculation of overlap matrices in reference [13]. This is an important difference with the implementation of [33] where an expansion on powers of the integrand is performed. Another difference of our implementation is that we write an interface for Wannier90, and hence the trial functions are the ones for Wannier90, while in reference [33] the original basis set is used.
Finally, we note that Brillouin zone sampling used to obtain the self-consistent electronic density in Siesta is essentially independent from the sampling used for the input data of Wannier90, i.e. Equations (4,5) and ǫ nk .

Ab-initio calculation of a cobalt impurity in bulk copper
The first system, the cobalt impurity in a bulk FCC copper host matrix, was simulated by a supercell where one host atom was replaced by a cobalt one. The lattice parameter was fixed to the theoretical value we found for pure copper (a lat = 3.690Å). We use two different unit cells to represent the system: (i) an eight-atom (2 × 2 × 2) cell § , where there is a 1/7 ratio of Co impurities, and (ii) a 64-atom one (4×4×4) where the cobalt concentration drops to 1/63. Besides its lower computational demand, the small cell permits us to represent the calculated bands and follow the Wannier disentanglement in a simpler way avoiding the larger band folding of the (4 × 4 × 4) supercell. According to the experimental data by Fan et al. [21] only above a concentration ratio of 1/4 does Co in Cu show ferromagnetic ordering; below this concentration, the system becomes paramagnetic. Our calculations correspond to impurity densities of the paramagnetic phase.
Interestingly, the measured saturation magnetic moment per atom is 1.5 µ B for pure Co and it drops to 0.4µ B /atom for a 1/9 concentration [21]. Our calculations yield 0.15µ B /atom for the 1/7 concentration, indicating a smaller spin-polarization of our DFT calculations. It is difficult to know the actual experimental error bar, but our calculations yield the right order of magnitude as well as the trend with decreasing Co concentration. Indeed, the supercell magnetization for the 2 × 2 × 2-case (1/7 concentration) is 1.03 µ B (0.15µ B /atom) while the supercell magnetization for the 4×4×4 is 1.57 µ B (0.025µ B /atom), indicating a saturation of the cell magnetization and thus a decrease of magnetization per atom, roughly linear with the number of electrons following the Slater-Pauling curve [21].
For the obtention of the MLWF, 55 bands in the 2 × 2 × 2 supercell and 460 in the case of a 4×4×4 supercell, were extracted for both spins, spanning the energy interval of (-14.0,9.1) eV with respect to the Fermi level (fully including the lowest valence bands). This energy interval is the outer window of Wannier90.
At the Fermi surface, the relevant states are both dispersive conduction bands and hybridized cobalt states originating from the Co incomplete d-shell. On the other hand, the closed copper d-shell gives rise to narrow bands, lying deeper below the Fermi energy (µ). In order to simplify the host electronic structure, we want these Cu d-bands to be projected out during the disentanglement process and keep a simpler sp-like band. This can be achieved by both (1) the choice of the inner window position that freezes the bands with dispersive character around µ and (2) the choice of trial orbitals. As for the latter, we used one spherically symmetric function per atom, located in one of the two interstitials of the fcc primitive cell. In addition, five cobalt centered orbitals with angular momentum l = 2 were used. The inner window choice is discussed in Section 3. Finally, for the spread, Ω, minimization a tolerance of 10 −10Å2 has been used.

Simulations of Cu (111)
The second system we studied is the Cu (111) surface. We used a 12-atom slab with a 1 × 1 surface unit cell. The k-point sampling was 14 × 14 × 1 and real space mesh cutoff of 800 Ry was used. Please, refer to [32] for more details. The Bloch functions corresponding to the lowest 120 bands were used to calculate the matrix M mn (k, b) using a k-point sampling 10 × 10 × 1 that largely suffices for convergence obtaining the MLWF. As in the previous case, we are interested in reproducing with the simplified scheme of MLWF the electronic structure about the Fermi energy, µ, hence we limit the inner energy window to the energy interval −1.0, 2.0 eV. The inner energy window must not contain more bands than the required number of Wannier functions, so the lower limit of the inner energy window is chosen somewhat above the narrow d-like bands and its upper edge of the inner window is limited by the higher-energy dispersive bands.
Among various configurations of trial orbitals, we found that the one described in figure 1 converges to good-quality MLWF. The figure shows the centers of the Wannier functions that are located both at the atomic positions of the slab (dots) and in between (crosses) to enhance the electronic structure description. Inside the slab, one spherical interstitial is sufficient to give a good description of the electronic structure close to the Fermi energy, µ. On the surface, two functions are necessary; we put them at the positions corresponding to the interstitials of the bulk geometry. The trial orbital set has been found stable: only the outermost Wannier functions move some 0.21Å outwards during the spread minimization. Table 1. Spread (Ω), on-site energies ( ǫ) and occupation (n) for the cobalt Wannier functions of a substitutional Co atom in the 2 × 2 × 2-Cu cell. These split into twoand three-fold degenerate states m = {e g , t 2g }. Majority (σ =↑) and minority (σ =↓) spins are given.

Cobalt impurity embedded in bulk copper
In this section, we present and discuss MLWF for the impurity system. The essence of what we are doing is a tight-binding representation of a conduction band (supported by one interstitial Wannier function per atom) hybridizing with impurity Wannier functions. This picture will be extended in the next section, where we attempt to separate the Coulomb interaction in the impurity in the spirit of Anderson Hamiltonian. Figure 2 shows the ab-initio band structure and the disentangled one for the 2×2×2cell. The color code denotes 100% overlap with a Co MLWF for red, and 0% for blue. The chosen inner energy window starts at −1.05 eV and extends up to +3.6 eV. This, in combination with trial orbital choice, efficiently disentangles the Cu sp-bands from the Cu d-bands, removes these last ones, and keeps the rest together with the full delectron structure of the Co impurity. As we described above, the spin polarization of the system is sizable due to the presence of Co atoms. This is clearly seen in the present graph, where the bands split according to their spin, in figure 2 (a), the majority spin is represented, and correspond to Co d-bands that coexist in the region of the Cu d-bands. The quality of these bands stemming from the MLWF calculation is considerably worse than the minority spin ones, figure 2 (b), because they lie outside the inner energy window. Nevertheless, the calculation contains information on the effect of the Cu dbands on the majority spin Co bands as we will discuss in the analysis of the PDOS. Moreover, the calculations dealing with the magnetic structure and other information close to the Fermi energy will be very accurate as can be seen in the excellent matching of the ab-initio bands and the MLWF ones near the Fermi level, figure 2. The calculation retains the Cu sp-bands that correspond to the colder-colored bands of figure 2 that strongly disperse. However, the hot-colored bands present flat features, revealing a small Co-Co intercell interaction, revealing that the dispersion largely comes from nearneighbor interactions, Co-Cu. For the 4 × 4 × 4 cell this interaction is even smaller, but the bands are not much flatter, what implies that for many Co-based properties our 2 × 2 × 2 calculations will suffice.
It can be seen that the MLWF disentanglement projects out the Cu d-bands but retains a realistic description of Co d-bands, as revealed by looking into the spin polarization of the system. Indeed, by evaluating (7), we can compute the electronic population of the different MLWF. For the 4 ×4 ×4 cell, the majority spin MLWF yields 4.59 electrons, while the minority spin gives a population of 2.65r;, if we add the spin polarization of the remaining Cu bands, we obtain 1.72 electrons, in good agreement with the 1.57 electrons of the full ab-initio calculation, discussed in section 2.4.
Valuable information is obtained by analyzing the PDOS in figures 3 and 4 and the MLWF occupancies, on-site energies and real-space spreads as shown in tables 1 and 2.
Due to the symmetry of the crystal field [22,34], we find two types of MLWF, the e ↑ g one at lower energy, which is twice degenerate according to the on-site energies (diagonal element of the MLWF Hamiltonian), and the one of t ↑ 2g symmetry, three-fold degenerate. We find that the minority spin MLWF's, e ↓ g , t ↓ 2g are actually different in spread from the majority spin ones, e ↑ g , t ↑ 2g . This is typical of unrestricted Hartree-Fock schemes, and it reflects the fact that both spins correspond to very different energy regions. The difference among MLWF's can be seen in the spread, Ω in tables 1 and 2 and in (3), that measures the extension of the MLWF. We see that e g MLWF's are more compact than the t 2g ones, and that the majority spin are more compact than the minority ones. It is interesting to note that the MLWF's are indeed very confined, in the present case, the more extended MLWF is basically zero beyond 2Å from its center.
The on-site energy is exactly the first moment of the spectral function (6), and hence can be directly correlated to figures 3 and 4. The occupancy is the integration of the spectral function, over occupied states, from −∞ to the Fermi energy, (7). Hence, these two quantities help us characterizing the spectral functions. As shown in figures 3 and 4, there are four sharp peaks that correspond to the above four types of MLWF. The presence of tails further from peak centers indicates that the MLWF's have an important hybridization with mainly the Cu host. These tails shift the on-site energy away from the peak energies. Most of the weight of the MLWF is, however, concentrated in a small energy window spanned by the peaks. Only in the case of the e ↑ g MLWF does the peak slightly lie outside the inner energy window used for the band disentanglement, and develops a tail in the Cu d-band region.
Small differences show up when going from the 2 × 2 × 2 cell to the 4 × 4 × 4 one. The main peaks lie at the same energy, and only the first momenta or on-site energies are slightly different (table 2) because the spectral function distribution is somewhat narrower for the smaller cell. This is due to the smaller content of Co atoms in the large cell, that leads to smaller Co-Co hoppings. These comparison permits us to conclude that the 2 × 2 × 2 cell is indeed a good approximation to the dilute (paramagnetic) regime, in agreement with the experimental findings by Fan et al. [21].
The fact that the spectral functions are rather localized in energies permit us to conclude that they are very good descriptors of the actual Co electronic structure. This is further seen in figure 5, where the MLWF spectral function is compared with the atomic basis one (3d Cobalt PAOs). Indeed, the PAO spectral function is more spread in energies, and reveal more mixing with the Cu-band structure, and, consequently, less Co character. It is also interesting to notice that the main peaks coincide, showing that the MLWF contain a lot of physics of the actual electronic structure.
From the spectral functions we see that the different electronic contributions of Co will have different physical properties. The t 2g electronic structure is basically fully occupied and slightly contributes to magnetism. However the e g states have a large spin polarization, carrying the leading rhôle in the magnetic properties of the Co impurities.

Model Hamiltonian
The results of the previous section show that the cobalt impurity physics can be qualitatively understood as an atomic d-orbital weakly hybridized with the Cu substrate. The hybridization determines the positions and widths of atomic levels, which in turn  fix the impurity occupancy. There is significant energy splitting regarding spin, which results in sizeable majority and minority spin populations. It is well understood that the origin of spin polarization lies in the intra-atomic Coulomb interaction [14]. However, the energy splitting for the e g levels is larger than for the t 2g ones, indicating that different t 2g and e g Coulomb matrix elements must be taken into account.
The substrate electrons (operators c † nkσ , c nkσ ) with Bloch energies ǫ nk hybridize with impurity given by H at with on-site energies ǫ m . Electrons on the impurity site are subject to Coulomb repulsion with matrix elements U mm ′ coupled to orbital occupancies n mσ = d † mσ d mσ . Interaction of equal spins is neglected by assuming the same value of direct and exchange integrals. The impurity and conduction band are connected by the hybridization term with matrix elements V nk,m .
Construction of the Hamiltonian (13) from first principles using MLWF for a cobalt impurity in copper substrate is now discussed. Orthogonality of Wannier functions allows us to divide the Kohn-Sham Hamiltonian into the blocks where H subs is the Hamiltonian in the subspace of substrate Wannier functions, H imp acts on the impurity Wannier functions (i. e. the cobalt e g and t 2g functions) while the remaining terms V, V † are off-diagonal. In order to be consistent with the mean-field character of the Kohn-Sham DFT we identify (15) with the mean-field solution of the Anderson Hamiltonian.
Firstly, H at is discussed. We impose that the on-site energies ǫ m obey restrictions of symmetry, i.e. have two values which will be denoted by ǫ e and ǫ t for e g and t 2g . In the correlation term, U mm ′ is a symmetric block matrix of the form A usual mean-field factorization leads to the simplified Hamiltonian where the modified atomic energies read The right hand sides are identified with on-site energies of impurity Wannier functions, for we take H MF at = H imp . The electron numbern eσ means the total mean occupancy of e g symmetry orbitals with spin σ and similarly for the t 2g states.
To build the Anderson Hamiltonian, the bare energies ǫ m and Coulomb integrals U mm ′ must be determined. Since the occupancies can be obtained from (7) we face a linear system of equations.
This system is underdetermined. We solve it by fixing the difference ∆ = ǫ t − ǫ e = 0.087eV (20) which is the crystal field splitting of atomic levels that is calculated as the t 2g − e g splitting of bulk copper d-orbitals . Taking the 4 × 4 × 4 case, we obtain: Eventhough we are not aware of any estimation of the U for cobalt in bulk copper, values reported in the literature for bulk Co are about 5 eV [7], and for cobalt adsorbed on the (111) surface of gold [37] are 2.8 eV. We expect our values to be smaller than the ones obtainable by the approach of [7]. The reason for this is that our approach gives a rule to obtain an Anderson Hamiltonian from the MLWF Hamiltonian which proceeds from our DFT calculation and hence has the known problems of the local and semilocal exchange and correlation functionals in determining U (see for example [4]). Namely, our U will correspond to the Hund's rule exchange matrix element rather than the electrostatic one appearing in LDA+U. Nevertheless, our values albeit smaller, are comparable to other U values for Co. Indeed, Antonides and co-workers measured 1.2 eV for the U of Co as an impurity [38] in excellent agreement with our calculation. We think that the main advantage of our approach is that it keeps the complete symmetry and multi-orbital character of the original DFT calculation.
The construction of substrate and hybridization terms in (13) is straightforward due to their one-particle form. Firstly, H subs is diagonalized, yielding the conduction In good agreement with the crystal field splitting in Ag and Au estimated in [ [35]] of about 0.15 eV, and in reference [36] of less than 0.1 eV. e g ↑ e g ↓ t 2g ↑ t 2g ↓ Figure 6. Hybridization function, Γ, calculated for the cobalt impurity in the 2 × 2 × 2 supercell of copper using a fine Brillouin zone sampling 40 × 40 × 40 and gaussian smearing 0.1 eV. The matrix elements corresponding to hybridization of e g states and t 2g states for both spins are given. Only diagonal contributions (see (24)) are shown because the non-diagonal hybridization functions are orders of magnitude smaller. Zero coincides with Fermi level.
band energies ǫ nk and Bloch states. In the next step, hoppings V and V † are transformed accordingly, leaving us with the matrix elements V nk,m between substrate Bloch states and impurity Wannier functions.
The essentially complete description of the impurity -substrate mixing is included in the hybridization function whose diagonal elements Γ mm give the spectral intensity of hybridization of a given impurity Wannier function m. The elements Γ mm (0) give (apart from a 2π factor) the inverse lifetime of an impurity electron. Off-diagonal terms provide substrate mediated intra-atomic hybridization intensity. The function Γ mm ′ (ω) is precisely − 1 π times the imaginary part of the retarded self-energy due to hybridization.
So far we have tacitly omitted spin polarization inherent in the matrices H subs , V and V † , coming from a spin polarized Kohn-Sham DFT. The hybridization function calculated for majority and minority spin directions is given in figure 6. Off-diagonal e g − t 2g matrix elements were omitted, for we found they are of the order of a few meV. The same holds for elements between different Wannier functions of the same symmetry. We see that all intensities have the same order of magnitude. The t 2g functions become much stronger for unoccupied substrate levels. The hybridization for minoritary spin is weaker than the majoritary-spin one. The most pronounced difference is between the e g ↑ and e g ↓ curves. The differences between hybridization functions for different spins partially reflect the polarization of the copper matrix by the magnetic moment of cobalt.  Figure 7 shows the comparison of the ab-initio band structure for the Cu (111) and the MLWF interpolated ones. The disentangling scheme has permitted us to retain the electronic structure with sp character about the Fermi energy. As a consequence neither the d-bands nor the upper limit of the LL' gap (Γ for the surface Brillouin zone) are described by the MLWF. However, the electronic structure in the inner energy window is perfectly reproduced.

Surface studies: the case of Cu (111)
In particular, we underline the excellent description of the Shockley surface state by the MLWF basis set. The band structure is basically indistinguishable from the ab-initio one which in turn is a very good description of the experimental one [32]. The minimization procedure leads to two types of differentiated MLWF sets, one set describing the bulk electronic structure with a spread of Ω = 3.3Å 2 and a surface MLWF with Ω = 5.7Å 2 and its original center displaced by 0.21Å into the vacuum region. Hence, the MWLF try to follow the behavior of the PAO basis [32] where the energy minimization was improved by using two distinct basis sets, one for the bulk electronic structure and one for the surface with diffuse orbitals. Indeed, not only does the surface MLWF reproduce the surface state dispersion (minimum at E 0 = −0.34 eV, and effective mass m * = 0.335 in electron masses, in good agreement with E 0 = −0.42 eV and m * = 0.37 of reference [32]) but it also shifts into the vacuum region in order to account for the surface spilling of charge.
The interpolated bands can be analyzed in terms of the MLWF character they have. Figure 7 depicts in a color code the character of the interpolated bands. In this way we find that the Shockley surface state is basically purely described by the surface MLWF nearΓ. As we move away fromΓ, the surface state has some bulk MLWF weight, until it reaches the gap edge and it becomes a surface resonance with a large bulk MLWF character.
Outside the inner energy window, the bulk MLWF has more weight. Indeed, the lowest and highest energy bands are described by the bulk MLWF. This signals that in this slab geometry, the electronic structure has mainly a surface character and bulk MLWF are the least indicated to describe the slab electronic structure.
The MLWF's correspond to s-like orbitals. A consequence of this is the disentanglement of the electronic structure as seen in figure 7. There, we see that the MLWF band structure spans from the bottom of the conduction band and crosses the d-band without mixing with it. Even outside the inner energy window, the MLWFbands excellently reproduce the sp-bands. The description considerably worsens outside the interpolating inner energy window.

Summary and conclusions
We have implemented an interface to the scheme of [12] to obtain a finite set of maximally-localized Wannier functions (MLWF's) from calculations using the Siesta code [13].
The method yields an efficient band disentanglement in the case of Co impurities in bulk Cu. We have been able to retain just an s-like band to represent the Cu electronic structure, while keeping most of the correct description of the Co electronic structure. As a consequence, we have analyzed the magnetic properties in the limit of small concentration of Co impurities, and trace back most of the Co properties to the ones described by a single MLWF, doubly degenerated with e g symmetry.
The reduction of the electronic structure to a few important elements is of great interest to obtain model Hamiltonians from full DFT calculations. We show that we can map the MLWF Hamiltonian into an Anderson model, having access to on-site energies, exact hybridization matrix elements and intra-atomic Coulomb matrix elements that reflect the correct symmetry of the problem. We show that by using this scheme, we obtain values of the intra-atomic Coulomb matrix that correspond to the semilocal exchange and correlation functional, and hence to the Hund's rule exchange, that are slightly smaller than the electrostatic ones obtained in LDA+U schemes [4,5,6]. These results are encouraging for future studies of strongly correlated systems using MLWF.
Even in the more stringent case of surface electronic structure, the description in terms of MLWF give very good results. We have applied the computational scheme to the Cu (111) surface and inside the chosen inner energy window the interpolation of the band structure is excellently reproducing the Shockley surface state with an accuracy comparable to the complete calculation of reference [32], and permitting us to obtained the sp-band disentangled from the Cu d-bands. In this way, the electronic structure problem is limited to the sp-electronic structure about the Fermi energy.
The surface case shows that not only is the MLWF method a mathematical trick to disentangle bands and reduce the problem to a smaller, tight-binding like, basis set. Indeed, one can analyze the obtained electronic structure in terms of the MLWF and have interesting insight. We have projected the electronic bands in the two types MLWF of the surface problem: the surface and the bulk MLWF's. In this way, we have verified the mainly surface character of the Shockley state near theΓ point, as well as its subsequent mixing with the bulk MLWF as the surface state energy increases, until it becomes a surface resonance.
In conclusion, we have implemented an interface to the scheme [12] to obtain a finite set of maximally-localized Wannier functions from calculations using the SIESTA code [13]. In this way, we can obtain DFT-based model Hamiltonians from ab-initio calculations with a very accurate description of the electronic structure inside a chosen energy window, that lends itself to the analysis and exploration of complex systems.