Projectability disentanglement for accurate and automated electronic-structure Hamiltonians

Maximally-localized Wannier functions (MLWFs) are broadly used to characterize the electronic structure of materials. Generally, one can construct MLWFs describing isolated bands (e.g. valence bands of insulators) or entangled bands (e.g. valence and conduction bands of insulators, or metals). Obtaining accurate and compact MLWFs often requires chemical intuition and trial and error, a challenging step even for experienced researchers and a roadblock for high-throughput calculations. Here, we present an automated approach, projectability-disentangled Wannier functions (PDWFs), that constructs MLWFs spanning the occupied bands and their complement for the empty states, providing a tight-binding picture of optimized atomic orbitals in crystals. Key to the algorithm is a projectability measure for each Bloch state onto atomic orbitals, determining if that state should be kept identically, discarded, or mixed into the disentanglement. We showcase the accuracy on a test set of 200 materials, and the reliability by constructing 21,737 Wannier Hamiltonians.


I. INTRODUCTION
In periodic crystals, the electronic structure is usually described using one-particle Bloch wavefunctions.While choosing a basis set that is also periodic to describe these wavefunctions can often be beneficial, an alternative approach is to adopt localized orbitals in real space.One such choice of orbitals are Wannier functions (WFs), that can be obtained by Fourier transforming the periodic wavefunctions from reciprocal to real space.WFs are not unique, as they depend on the choice of the gauge (i.e., the choice of the phases of the wavefunctions) at each point in the Brillouin zone (BZ).Maximallylocalized Wannier functions (MLWFs) [1][2][3][4] are obtained by a gauge choice that is optimized to provide the most localized set of WFs, i.e., those that minimize the sum of their quadratic spread in real space [1].Having a very localized representation of the electronic structure not only provides an insightful analysis of chemical bonding in solids, but also brings a formal connection between the MLWF centers and the modern theory of electric polarization [5].Moreover, the real-space locality of MLWF allows for accurate and fast interpolation of physical operators [6], enabling calculations of material properties that require dense samplings of the BZ, such as Fermi * junfeng.qiao@epfl.chsurface, orbital magnetization [7], anomalous Hall conductivity [8,9], and spin Hall conductivity [10], to name a few.Practically, one obtains MLWFs starting from a set of Bloch wavefunctions, calculated e.g., from densityfunctional theory (DFT).Often, these Bloch states are projected onto some localized orbitals (usually chosen by the user) to generate initial guesses for MLWFs.In an insulator, by minimizing the spread functional [1] which measures localization, one can obtain a set of MLWFs, i.e., "Wannierize" a material.The Wannierization contains an additional disentanglement step [2] if the target Bloch states are not isolated from other band manifolds.For such entangled bands-metals or the conduction bands of insulators-one needs to first identify the relevant Bloch states that will be used to construct ML-WFs, and then mix or "disentangle" these from all the Bloch states [2].Practically, the choices for the initial projections and states to be disentangled substantially influence the shape and the quality of the final MLWFs.
In recent years, a lot of effort has been devoted to obtaining high-quality MLWFs and automate the Wannierization procedure.Focus of the research can be categorized into the following classes: (a) Novel minimization algorithms, such as: the symmetry-adapted WF method that adds constraints to impose the symmetries of the resulting WFs [11]; the simultaneous diagonalization algorithm that directly minimizes the spread functional for an isolated (or "Γ-only") system [12]; the partlyoccupied WF method, where the total spread is directly minimized in one step [13,14], rather than performing a two-step minimization for its gauge-invariant and gaugedependent parts as in the standard procedure [2]; or the variational formulation, that combines single-step optimization with manifold optimization to make the minimization algorithm more robust [15]; (b) new forms for the spread functional, such as the selectively localized WFs (SLWFs) for which only a subset of WFs of interest are localized and a penalty term is added to constrain the position of the WF centers [16], or the spreadbalanced WF method, that adds a penalty term to distribute the spread as uniformly as possible among all WFs [17]; (c) targeting a subset of orbitals, e.g.SLWF for a subset of MLWFs [16] or the optimized projection functions method where starting projections for the Wannierization are generated from a larger group of initial ones [18]; (d) matrix manifold algorithms instead of projection methods to construct a smooth gauge in a noniterative way [19,20]; (e) basis-vector decomposition of the density matrix, e.g. the selected columns of the density matrix (SCDM) algorithm [21,22], that starts from the density matrix of the system and uses QR decomposition with column pivoting (QRCP) to automatically generate an optimal set of basis vectors from the columns of the density matrix.
At the same time, high-throughput (HT) calculations have become increasingly popular for materials discovery and design.Calculations and results managed by workflow engines are collected into databases of original calculations, such as the Materials Project [23], AFLOW [24], OQMD [25], CMR [26], and the Materials Cloud [27], or aggregated, as in NOMAD [28].Thanks to recent research advances on Wannierization algorithms, it starts now to be possible to run HT Wannierizations for many materials and generate tight-binding (TB) models that reliably describe their physics.So far, several attempts have been made in this direction.Gresch et al. [29] gathered 195 Wannier TB Hamiltonians and applied postprocessing symmetrization to study strained III-V semiconductor materials.Vitale et al. [30] implemented the SCDM algorithm and designed a protocol to determine automatically the remaining free parameters of the algorithm; this protocol, implemented into automated workflows, was verified to work well for band interpolations on a set of 200 structures (metals, or valence and conduction bands of insulators) and 81 insulators (valence bands only).Garrity and Choudhary [31] accumulated a Wannier TB Hamiltonian database of 1771 materials using the standard hydrogenic orbital projections.However, there are still several challenges for an accurate and automated HT Wannierization, some of which might be more relevant depending on the research goal and the specific property to compute: MLWFs should be able to faithfully represent the original band structure, often (e.g., for transport properties) at least for those bands close to the Fermi energy; MLWFs should resemble the physically intuitive atomic orbitals for solids that would enter into Bloch sums; the algorithm should be fully and reliably automated and the implementation should be efficient for HT calculations.
To overcome the challenges mentioned above, in this paper we present a new methodology for automated Wannierization.First, we choose physically-inspired orbitals as initial projectors for MLWFs, that is, the pseudo-atomic orbitals (PAOs) from pseudopotentials [32].Then, for each state |nk⟩ (n is the band index, k is the Bloch quasi-momentum) we decide if it should be dropped, kept identically, or thrown into the disentanglement algorithm depending on the value of its projectability onto the chosen set of PAOs, replacing the standard disentanglement and frozen manifolds based only on energy windows.This approach naturally and powerfully targets the TB picture of atomic orbitals in crystals, as it will also become apparent from our results.Moreover, we fully automate this approach and implement it in the form of open-source AiiDA [33][34][35] workflows.To assess its effectiveness and precision, we compare the quality of the band interpolation and the locality of the Wannier Hamiltonians generated with the present approach, which we name as projectability-disentangled Wannier functions (PDWFs), with the results from the SCDM algorithm [30].Statistics from 200 materials demonstrate that PDWFs are more localized and more atomic-like, and the band interpolation is accurate at the meV scale.Furthermore, to demonstrate the reliability and automation of our method and workflows, we carry out a largescale high-throughput Wannierization of 21,737 materials from the Materials Cloud [27,36].
To set the context for the following paragraphs, here we briefly summarize the notations for WFs; a detailed description can be found in Refs.[1][2][3].WFs |w nR ⟩ are unitary transformations of Bloch wavefunctions |ψ mk ⟩, given by where k and R are the Bloch quasi-momentum in the BZ and a real-space lattice vector, respectively; m is the band index, and n is the Wannier-function index (running from 1 to the number of WFs J).For an isolated group of bands, J is equal to the number of bands, and the U mnk are unitary matrices; for entangled bands, the number of bands considered at each k−point is J k ≥ J, and the U mnk are semi-unitary rectangular matrices.MLWFs are the minimizers of the quadratic spread functional [1] Since Eq. ( 2) is a minimization problem with multiple local minima, initial guesses for U mnk substantially influence the optimization path and the final minimum obtained.In order to target the most localized and chemically appealing solution, Marzari and Vanderbilt [1] used hydrogenic wavefunctions |g n ⟩ (i.e., analytic solutions of the isolated hydrogenic Schrödinger equation) to provide a set of sensible initial guesses |ϕ nk ⟩, after projection on the space defined by the relevant Bloch states: The projection matrices A mnk = ⟨ψ mk |g n ⟩, after Löwdin orthonormalization [37], form the initial guesses for U mnk .We underline that while the gauge of Bloch wavefunctions |ψ mk ⟩ is arbitrary, Eq. ( 3) is invariant to such gauge freedom: suppose |ψ ′ ik ⟩ are also solutions of the electronic structure problem, then |ψ ′ ik ⟩ are related to |ψ mk ⟩ by some unitary matrices does not depend on the gauge of Bloch wavefunctions, where superscript * denotes conjugate transpose.For entangled bands, the "standard" disentanglement approach [2] uses energy windows to choose the disentanglement and frozen manifolds: (a) an (outer) disentanglement window that includes a large set of Bloch states, which can be mixed together to obtain a smaller disentangled manifold; (b) an (inner) frozen window that specifies a smaller set of Bloch states (often states around Fermi energy) which are kept unchanged in the final disentangled manifold.
Since in the following sections the present results are compared with SCDM, we also summarize the SCDM procedure here.The SCDM method [21] starts from the real-space density matrix ⟨r|P k |r ′ ⟩ where P k = J k m=1 |ψ mk ⟩ ⟨ψ mk |, and uses QR factorization with column pivoting (QRCP) to decompose ⟨r|P k |r ′ ⟩ into a set of localized real-space orbitals, thanks to the nearsightedness principle [38,39] stating that the matrix elements ⟨r|P k |r ′ ⟩ decay exponentially with the distance between two points r and r ′ in insulating systems.While storing the full ⟨r|P k |r ′ ⟩ is memory intensive (it has size N r × N r , where N r is the number of real-space grid points), one can equivalently decompose the matrix formed by the real-space Bloch wavefunctions Ψ * k = [ψ 1k , . . ., ψ J k k ] * , which has a smaller size J k × N r .For periodic systems, often the choice of columns in the QRCP algorithm can be performed using the wavefunctions at the Γ point only (Ψ Γ ) [40], and the same column selection is then used for all other kpoints.For entangled bands, since the density matrix is not continuous across the k-points, one can construct a quasi-density matrix (or equivalently a matrix of wavefunctions) , where f (ε mk ) is a smooth function of the energy eigenvalues ε mk , specifying the target energy window for the constructed MLWFs.Often the complementary error function 1  2 erfc( ε−µ σ ) is chosen as f (ε), and the choice of µ and σ determines the shape of MLWFs, as well as band-interpolation quality.Using projectability, defined later in Eq. ( 5), µ and σ can be automatically chosen, thus automating the Wannier-ization process [30].

A. Pseudo-atomic-orbital projections
In addition to the hydrogenic orbitals discussed above, alternative starting guesses for the Wannierization can be used.For instance, in pseudopotential plane-wave methods, PAOs are localized orbitals originating from the pseudopotential generation procedure [32].In this procedure, for each element, atomic wavefunctions of an isolated atom are pseudized to remove the radial nodes and are localized functions around the atom; spherical harmonics with well-defined angular-momentum character (s, p, d, or f ) are chosen for their angular dependency.Then, the PAOs are summed over lattice points with appropriate phases to obtain Bloch sums, Fourier transformed to a plane-wave basis, Löwdin-orthonormalized, and finally taken as the projectors for initial projections.PAOs are commonly used for analyzing the orbital contributions to band structures, as the basis set for noniterative construction of TB Hamiltonians [32], or as projectors in DFT+Hubbard calculations [41].
In order to understand the contribution of each orbital |g n ⟩ to a Bloch state |ψ mk ⟩, we define a measure of projectability as the square of the inner product between |ψ mk ⟩ and |g n ⟩: the projectability of |ψ mk ⟩ onto all PAOs is then defined as If the projectors |g n ⟩ are complete for |ψ mk ⟩, then p mk = 1.The band projectability is a very useful criterion to identify the orbital character of the bands; this is exemplified in Fig. 1a, where we show the projectability of the bands of graphene onto 2s and 2p PAOs for carbon.It is immediately apparent how one can easily identify states in the conduction manifold that have a strong 2p and 2s component.
Compared with the hydrogenic projections, which is the method used by default in Wannier90 [4] and its interface code to Quantum ESPRESSO [42] (called pw2wannier90.x),PAOs are better adapted to each element since they come exactly from the pseudopotential used in the actual solid-state calculation.Moreover, in pseudopotentials with semicore states, the PAOs for semicores are nodeless and those for valence wavefunctions have at least one radial node (so as to be orthogonal to the semicore states with same angular momentum); thus band projectability can clearly differentiate semicore from valence, making PAOs more convenient than the hydrogenic orbitals, for which the user would need to manually set the correct radial functions for both semicore and valence projectors.For these reasons, we use in this work the PAOs as initial and more accurate projections.If needed, higher energy orbitals not included in the pseudopotential file can be constructed, for example, using solutions of Schrödinger equation under confinement potential [43,44] (see also discussion in Section II F).

B. Projectability disentanglement
As mentioned, the standard disentanglement approach selects the disentanglement and frozen manifolds via two energy windows [2].We refer to this as energy disentanglement (ED).However, since bands have dispersions across the BZ, a fixed window for all k-points might not be an optimal choice.Taking the graphene band structure (Fig. 1a) as an example, the bands with large projectability are mixed with many free-electron bands with zero projectability (grey bands in the conduction region).In this case, one is faced with several options for the outer and inner energy windows, each with different shortcomings: (a) If the inner window includes freeelectron bands, the final MLWFs are mixtures of 2s, 2p atomic orbitals and free-electron bands, delocalizing the resulting MLWFs; (b) if the outer window excludes both the free-electron bands and the atomic-orbital states inside free-electron bands, the WFs lack the anti-bonding part of the bonding/anti-bonding closure [13], again degrading the localization of WF; (c) if the upper bound of the inner window is set to its maximal allowed value, i.e. the blue dashed line positioned at the minimum of freeelectron bands in Fig. 1b, and all the DFT eigenstates are included in the outer window, the disentanglement algorithm [2] will extract an optimally smooth manifold, at the expense of decreasing the chemical representability of the atomic-orbital bands in the free-electron region; in other words, the MLWFs obtained lose the information of the TB atomic orbitals in this chemical environment (see Fig. 1b).
The graphene case highlights the limitations of the standard ED.Instead, we propose here to select the disentanglement and frozen manifolds based on the projectability p mk of each state on the chosen PAOs (i.e., states are selected irrespective of their energy, but rather based on their chemical representativeness).Specifically, we select states based on two thresholds p min and p max : (a) If p mk < p min , the state ψ mk is discarded.(b) If p mk ≥ p max , the state ψ mk is kept identically.Crucially, all states for which p min ≤ p mk < p max are thrown in the disentanglement algorithm.Optimal numerical values for p min and p max are discussed later.In the case of graphene, p max identifies the fully atomic-orbital states inside the free-electron bands, while p min removes the fully free-electron bands from the disentanglement process, preventing the mixing of atomic and free-electron states.The two thresholds p min and p max constitute the parameters of the disentanglement process, replacing the four defining energy windows (the lower and upper bounds of the outer and inner energy windows).We note that projectability disentanglement is different from partly-occupied WF [13,14] in that the latter uses an energy window to select frozen states and minimizes the total spread functional directly, while projectability disentanglement selects the localized states using projectability instead of a constant energy window across k-points.In fact, one can combine projectability disentanglement with a variational formulation [15] to construct MLWFs by minimizing directly the total spread functional.
Ideally, if PAOs were always a complete set to describe valence and near-Fermi-energy conduction bands, the PD would select the most relevant Bloch states and accurately interpolate these DFT bands.However, since the PAOs are fixed orbitals from isolated single-atom calculations for each element, if the chemical environment in the crystal structure is significantly different from that of pseudopotential generation, then the total projectability p mk might be smaller than 1 for bands around the conduction band minimum (CBM) or even for valence bands.In such cases, one solution is to increase the number of PAOs, i.e., adding more projectors with higher angular momentum, as we will discuss in Section II F. However, since one almost always wants to correctly reproduce valence bands (plus possibly the bottom of the conduction) but at the same time keep the Wannier Hamiltonian small for computational reasons, we suggest to additionally freeze all the states that sit below the Fermi energy in metals (or below the CBM for insulators) and also those a few eV above (typically, 2 eV or so).Such a combination of PD+ED gives accurate interpolation of bands below and around the Fermi energy (or band edges for insulators), as well as maximally restoring the atomic-orbital picture.
We stress here that, even if we call the resulting Wannier functions PDWFs for clarity, our optimal suggestion is to always also freeze the states in the energy window mentioned above, as we discuss in the next sections.

C. Comparison
We choose four prototypical materials to discuss the present method: graphene, silicon, copper, and strontium vanadate (SrVO 3 ).Graphene is a difficult case where atomic-orbital states highly mix with free-electron bands; silicon tests the Wannierization of both valence and conduction bands of an insulator; copper is a test on a metal; and SrVO 3 represents the class of (metallic) perovskites.We compare the shapes, centers, and spreads of the resulting MLWFs using the five methods mentioned earlier: hydrogenic projection with ED (i.e., the standard approach), SCDM, PAO projection with ED, PAO projection with PD, and PAO projection with PD+ED.

Graphene
The original and interpolated band structures for the five methods discussed are shown in Figs.1b to 1f.The blue dashed lines in Figs.1b, 1d and 1f indicate the top of the inner energy window, which is set optimally (and manually) to just below the free-electron bands, to freeze as much as possible the atomic-orbital states but exclude any free-electron state.For PD and PD+ED, we choose p max = 0.85 and p min = 0.02 (we will discuss later on the choice of these thresholds).Comparing Fig. 1d and Fig. 1b, one sees that ED produces similar bands irrespective of using hydrogenic or PAO projection.However, as shown in Fig. 2 (first and third row), the ML-WFs for the two cases fall into slightly different minima: MLWFs from hydrogenic projection with ED are p z and hybridized s ± p orbitals pointing towards the center of the hexagon, while MLWFs from PAO with ED are p z , p x , and s ± p y .This is due to the fact that the PAO projections guide the minimization towards spherical har- FIG. 2. Graphene MLWFs: shapes, centers, and spreads obtained using different methods.dWFC is the distance of the WF center from the nearest-neighbor atom, and ΩWF is the MLWF spread.The multiplicity is the number of equivalent MLWFs, i.e. having the same dWFC, ΩWF, and shape, but different orientations.
monics, while the hydrogenic projections are farther away from such local minimum and the optimization algorithm happens to escape and converge to a better minimum.
A possible future work is to introduce more advanced optimization algorithms to improve the convergence of maximal localization.Both the PAO with PD and PAO with PD+ED cases reach the same set of MLWFs, p z , p x , and s ± p y , but with larger spreads than the PAO with ED, since the PD and PD+ED freeze more states, giving thus less freedom for maximal localization.Nevertheless, the interpolated bands of the PAO with PD and PAO with PD+ED cases can much better reproduce the atomic-orbital states inside the free-electron bands.Finally, compared to other cases, SCDM includes some freeelectron bands, some of which can be even reproduced by the Wannier interpolation.However, in order to follow those free-electron bands, abrupt changes of character and band derivative are needed in the conduction band.
As required by Nyquist-Shannon sampling theorem [45], this results in a denser k-space sampling needed to obtain a good interpolation quality.Moreover, the MLWFs are much more delocalized and do not resemble atomic orbitals: as shown in Fig. 2, the last two MLWFs for SCDM are floating away from the graphene 2D lattice, blurring the TB picture of atomic orbitals in solids.

Silicon
The SCDM method obtains four front-bonding and four back-bonding MLWFs, while all other cases lead to atom-centered s and p MLWFs, as shown in Fig. S3.While overall the SCDM bands (Fig. 3c) seem to reproduce relatively better the higher conduction bands, they fail to correctly reproduce the bottom of the conduction band near the X point, induce more wiggles around X and W, and have much larger spreads.Due to the low projectability of Bloch states around X (p mk around 0.83), the CBM is not correctly reproduced in the PAO with PD, as these are not frozen in PD with the current choice of p max = 0.95 and p min = 0.01.To explicitly freeze the CBM, p max would need to be lowered below 0.83.However, such kind of decrease will also result in freezing some high-energy conduction bands, degrading the localization.PD+ED overcomes this by explicitly freezing the near-Fermi-energy and low-projectability states at the CBM, but still only freezing those atomic-orbital states in the high-energy conduction bands that possess high projectability (see Fig. 3f), thus improving band interpolation.We note that the lower projectability of silicon CBM is intrinsic to the material-its CBM also includes 3d character.Therefore, by adding d PAOs, the CBM projectability increases (from 0.83 to 0.99) and one can restore a high-quality band-structure interpolation within the PD method: as shown in Fig. 3e, the low-energy conduction bands are correctly reproduced once we regenerate a silicon pseudopotential including 3d PAOs.Therefore, PD is sufficient to obtain an accurate band interpolation if enough PAOs are included (we will also discuss this later in Section II F).For completeness, we show the SCDM interpolation using the regenerated pseudopotential in Fig. 3c: the added d PAOs help select a larger manifold thanks to the increased projectability, enabling SCDM to reproduce higher conduction bands, as well as fixing the wrong interpolation at the W point.Moreover, additional PAOs can also benefit ED, since the frozen window can be enlarged to reproduce more states.In general, adding more PAOs improves interpolation quality in cases where the target bands have low projectability, at the price of increased computational cost.PD+ED is a better option for reaching a good interpolation accuracy while keeping the size of the corresponding TB model small.

Copper and SrVO3
Results for copper and SrVO 3 are only shown in the SI (Figs.S4, S6, S7 and S9), since the conclusions are the same: PD+ED consistently provides the best interpolation quality among all methods we consider, while not requiring to increase the size of the Hamiltonian model, and results in WFs that resemble atomic orbitals or their hybridization.

D. High-throughput verification
In this section we discuss the applicability of the present PDWF method to obtain, in a fully automated way and without user input, WFs for any material.In order to assess quantitatively its performance, we compare it to SCDM, that can also be fully automated (see Ref. [30]).
In all results that follow, we exclude semicore orbitals in both methods, since these low-energy states correspond to almost flat bands and do not play any role in the chemistry of the materials.We compare quantitatively the band interpolation quality between the two methods and the corresponding WF centers and spreads on the 200-structure set used in Ref. [30] for both oc-cupied and unoccupied bands, totalling 6818 MLWFs for each method.In accordance with Refs.[30,46], the band interpolation quality is measured by the average band distance, and the max band distance, where fnk = f DFT nk (E F + ν, σ)f Wan nk (E F + ν, σ) and f (E F +ν, σ) is the Fermi-Dirac distribution.Here E F +ν and σ are fictitious Fermi levels and smearing widths which we choose for comparing a specific range of bands.Since the Wannier TB model describes the low-energy valence electrons, it is expected that the band interpolation deviates from the original in the higher conduction band region.Therefore, the higher ν is, the larger η ν is expected to be.In the following paragraphs, we will use η 0 and η 2 to compare bands below E F and E F + 2 eV, respectively; σ is always fixed at 0.1 eV.
In the supplementary information Section S8, we provide comparisons between the Wannier-interpolated bands and the DFT bands for both PDWF and SCDM, their respective band distances, and the Hamiltonian decay plots for each of the 200 materials.We discuss these properties in the following.

Projectability thresholds and automation
For PDWF, we set the maximum of the inner window to the Fermi energy + 2 eV for metals, or to the CBM + 2 eV for insulators, to fully reproduce states around Fermi energy or the band edges.We also specify the two additional parameters p min and p max .From our tests, in most cases p max = 0.95 and p min = 0.01 already produce very good results.However, since chemical environments vary across different crystal structures, the two parameters are not universal and influence the quality of band interpolation.Figure 4 shows the variation of band distances w.r.t.p min and p max for several materials.For Al 3 V (Figs. 4a and 4b), η 0 and η 2 reach a minimum at two different sets of parameters, i.e., p max = 0.99, p min = 0.01 and p max = 0.97, p min = 0.01, respectively.In some cases, the variation of η w.r.t.p max and p min can be non-monotonic and display multiple local minima: For instance, in Au 2 Ti (Fig. 4c) at p min = 0.01, η 2 decreases from p max = 0.90 to 0.95 but increases from p max = 0.95 to 0.98 and finally reaches a local minimum at p max = 0.99.In other cases, η can be quite stable and largely independent of the parameters: e.g., for Ba 6 Ge 10 (Fig. 4d), η 2 reaches the same minimum for p max = 0.99 to 0.88.
Therefore, we implement an iterative optimization workflow to automatically find the optimal values for p max and p min , in order to fully automate the Wannierization procedure.The workflow is released as part of the aiida-wannier90-workflows package [47].First, we run a QE band structure workflow to get the reference DFT bands for calculating η 2 ; in addition, the DFT bands are also used to calculate the band gap of the material.Second, we run an optimization workflow with the following settings: The maximum of the inner window is set to Fermi energy + 2 eV for metals and CBM + 2 eV for insulators, respectively; p max and p min are set to the defaults of 0.95 and 0.01, respectively.Third, if the average band distance η 2 is less than a threshold (set to 10 meV here), the workflow stops; otherwise, the workflow iterates on a mesh of p max and p min , i.e. p max decreasing from 0.99 to 0.80 with step size -0.01, and p min = 0.01 or 0.02, until η 2 ≤ threshold.If η 2 is still larger than the threshold after exhausting all the parameter combinations, the workflow will output the minimum-η 2 calculation.

Band distance
To compare quantitatively the band interpolation quality of SCDM and PDWF, we Wannierize the 200 structures mentioned earlier and calculate their band distances with respect to the corresponding DFT bands.We choose η 2 and η max 2 to compare near-Fermi-energy bands.The histograms of the band distances for the 200 structures are shown in Fig. 5.To directly compare SCDM and PDWF, the mean and median value of η of the 200 calculations are shown as vertical lines in each panel.For PDWF, the mean η 2 is 4.231 meV, to be compared with 11.201 meV for SCDM.For η max 2 (that is a more stringent test of the quality of interpolation) the PDWF method also performs better, with a η max 2 = 36.743meV vs. 84.011meV for SCDM.We can also observe this trend in Fig. 5: For η 2 and η max 2 , the PDWF histogram bins are much more clustered towards η = 0. Note that in the cumulative histograms of η 2 , at η = 20 meV, the PDWF cumulative count is closer to the total number of calculations (200).This indicates that the PDWF has a higher success rate in reducing the interpolation error below 20 meV.Similarly, for η max 2 , PDWF has a higher success rate in reducing the interpolation error under 100 meV (to get a better overview of η and η max , we further show the same histograms of η in a wider range 0 meV to 100 meV, and η max in range 0 meV to 500 meV, in Figs.S11 and S12).To reduce the effect of major outliers, we can also compare the interpolation accuracy of successful calculations, i.e., excluding the outlier calculations which have significantly large band distances.As shown in Table S1, the η ≤20 2 , i.e., the average of all the calculations for which η 2 ≤ 20 meV, indicates that PDWF (2.922 meV) is twice as good as SCDM (5.280 meV), and also has a higher success rate: for η ≤20  S1.
In summary, PDWF provides more accurate and robust interpolations, especially for bands around the Fermi energy or the band gap edges, which are the most relevant bands for many applications.Last but not least, a higher energy range can be accurately interpolated by increasing the number of PAOs (see Section II F).

MLWF centers
Since we are aiming at restoring a tight-binding atomic-orbital picture with PDWF, we compare the distance of the WF centers from the nearest-neighboring (NN) and next-nearest-neighboring (NNN) atoms, again both for SCDM and PDWF.For each method, we compute d NN and d NNN , i.e., the average distance of all the 6818 MLWFs from the respective NN and NNN atoms.If d NN is 0, then the atomic-orbital picture is strictly preserved.However, this is unlikely to happen since there is no constraint on the WF centers during both the disentanglement and the localization, and the final PDWFs, resembling atomic orbitals, are optimized according to the chemical environment.Still, if a WF center is much closer to the NN atom than to the NNN atom, then one can still assign it to the NN atom, preserving the atomic-orbital picture.Figure 6   This can also be observed in Fig. 6: The overlap of the d NN and d NNN histograms is smaller for PDWF than for SCDM.To further understand the overlaps, we plot the histogram of the ratio d NN /d NNN of each MLWF in the insets of Fig. 6.For a MLWF, if d NN /d NNN = 1, then the MLWF is a bonding orbital centered between two atoms; while if d NN /d NNN ≪ 1, then it can be regarded as an (almost) atomic orbital.The histogram of the ratio of SCDM has a long tail extending towards 1.0, i.e., there are a large number of SCDM MLWFs sitting close to bond centers; on the contrary, the vast majority of the PDWF MLWFs are closer to the NN atom.
We can further compare the effect of maximal localization on the WF centers.The WFs from the projection matrices A mnk are strictly atom-centered, i.e. d NN = 0.The inset of Fig. S13a shows the histogram of the initial WFs, i.e., after disentanglement and before maximal localization, and the final MLWFs, i.e., after maximal localization, for PDWF.If one chooses d NN ≤ 0.1 Å as the criterion for atom-centered MLWFs, then 5594/6818 = 82.0% of the initial WFs and 2045/6818 = 30.0% of the final MLWFs are atom-centered.The disentanglement and maximal localization improve the band interpolation, but since there is no constraint on the WF center in the spread functional Eq. ( 2), many of the final MLWF centers are not atom-centered.As a comparison, for SCDM, 955/6818 = 14.0% of the initial WFs and 1823/6818 = 26.7% of the final MLWFs are atom-centered.For completeness, the statistics and histograms of initial and final d NN , d NNN , and d NN /d NNN are shown in Table S2 and Fig. S13.
In summary, for PDWF, most of the initial WFs (after disentanglement and before maximal localization) are atom-centered; many drift a bit away from atom centers during the localization, but the MLWFs are still much closer to the NN than to NNN atoms.For SCDM, most of the initial WFs are away from atom centers, and maximal localization pushes some of the WFs back to atoms, but there is still a large number of MLWFs for which an atom representing the WF center cannot be clearly identified.To exactly fix the MLWFs to atomic positions, one needs to add constraints to the spread functional [16], at the cost of potentially having worse interpolators.However, this is beyond the scope of the current work, and here we rely on the atom-centered PAO projectors to guide the MLWFs towards the atomic positions, so that the final MLWFs are optimally localized and atom-centered.

MLWF spreads
Next, we investigate the spread distributions of SCDM and PDWF.Usually, we want localized MLWFs to restore the TB atomic orbitals.Figure 7 shows the histograms of the spread distributions for the two methods.The SCDM spreads have a long tail extending over 10 Å2 in Fig. 7b, due to its inclusion of free-electron states in the density matrix, thus resulting in more delocalized MLWFs as discussed earlier (see e.g.Fig. 2).On the contrary, the PDWF selects and freezes atomic-orbital states from the remaining bands, leading to much more localized MLWFs, thus much more clustered in a narrow range of 0 Å2 to 4 Å2 , and already at 5 Å2 the cumulative histogram almost reaches the total number of MLWFs (see Fig. 7a).This can be interpreted as follows: The PAO initial projections guide the spread minimization toward the (local) minimum resembling spherical harmonics, whereas the SCDM-decomposed basis vectors are designed to be mathematical objects spanning as much as possible the density matrix, but result in WFs for which it is harder to assign definite orbital characters.
We can further compare the average initial (after disentanglement but before maximal localization) and final (after disentanglement and maximal localization) spreads between the two methods, as shown in Table S3 and corresponding histograms in Fig. S14.Maximal localization is needed to bring SCDM spreads, from the initial Ω i = 30.82Å2 to the final Ω f = 3.54 Å2 ; For PDWF, the initial Ω i = 2.72 Å2 is already excellent, and much better than the final Ω f for SCDM; localization then brings it to an optimal Ω f = 1.41 Å2 .

Hamiltonian decay
Finally, we compare the decay length of the Wannier gauge Hamiltonian between the two methods in Fig. 8. Thanks to the localization of MLWFs, the expectation values of quantum mechanical operators in the MLWF basis, such as the Hamiltonian H(R), decay rapidly with respect to the lattice vector R (exponentially in insulators [48,49] and properly disentangled metals).To compare this decay for the Hamiltonian matrix elements, we approximate the Frobenius norm of the Hamiltonian as where τ measures the decay length.Then τ is fitted by least squares to the calculated ∥H(R)∥; as shown in Fig. 8a, the Hamiltonian of PDWF decays faster than SCDM for Br 2 Ti, which is selected here to represent the general trend between PDWF and SCDM Hamiltonians.Fig. 8b shows the histogram of τ for the 200 materials; the mean τ are 2.266 Å for PDWF and 2.659 Å for SCDM, respectively, indicating that the PDWF Hamiltonian decays faster than SCDM, consistent with the better band interpolation of PDWF discussed in Fig. 5.

E. High-throughput Wannierization
Based on the above verification, we run a HT Wannierization using PDWF for 21,737 materials, selected from the non-magnetic materials of the MC3D database [36].Figure 9 shows the band distance histograms for η 2 and η max 2 .Overall, the statistics follow the same trend as the 200 materials set in Fig. 5: the average η 2 and average η max 2 are 3.685 meV and 42.768 meV, respectively.Note in Fig. 9a the η 2 is not truncated at 10 meV, but rather due to the automated optimization workflow: results that have η 2 larger than a threshold (10 meV) are further optimized with respect to p min and p max , thus improving the average band distance η 2 .In Table S4 we show several other statistics for the band distances.
The excellent interpolation quality of PDWF can be assessed, for instance, from the number of systems with η 2 ≤ 20 meV, that are ≈ 97.8% of all the calculations (21259/21737); the corresponding bands distance calculated on these 21259 calculations is η ≤20 2 = 2.118 meV.This remarkable result show how automated and reliable Wannierizations can now be deployed automatically both for individual calculation and for HT application.

F. Additional PAOs for high-energy high-accuracy interpolation
Based on the HT Wannierization results, one can identify cases where the interpolation quality can be further improved by increasing the number of PAOs.Typically, the number of PAOs is determined during pseudopotential generation, and they are usually the orbitals describing low-energy valence electrons.In some cases, the bonding/anti-bonding combinations of these PAOs are not sufficient to span the space of target conduction bands, leading to a loss of interpolation quality.We use silicon as an example to illustrate the difficulties of accurately describing its CBM [50], which is not located at any high-symmetry k−point, but along the Γ−X line.The common choice of one s and three p hydrogenic or PAOs projectors per atom results in oscillations in the Wannier-interpolated bands at the meV level.To remedy this, one can use a larger set of PAOs, e.g., by regenerating a silicon pseudopotential including d PAOs as discussed in Section II C 2. However, generating a new pseudopotential requires extensive testing and validation, therefore another solution could be using a set of PAOs different from the pseudopotential ones.To compare this second approach, we test here also PAOs obtained from the OpenMX code [44], and Wannierize silicon using one s, three p, and five d PAOs per atom using ED.This provides a much better description of the CBM, as shown in Fig. S17 Moreover, the additional d orbitals allow to raise the inner energy window and better reproduce a larger number of conduction bands, as shown in Fig. S18, which might be beneficial for some applications.For completeness, we also show the WF spreads and shapes of d orbitals in Fig. S19.However, there are some caveats to this approach.When using external PAOs, ideally one should generate them using the same pseudization scheme as the pseudopotentials used in the DFT calculations.The PAOs from OpenMX are instead generated using a different scheme, resulting in lower projectabilities (smaller than one even for the valence bands, as shown in Fig. S21).In such case, PD cannot reproduce the original bands (see Fig. S20b), thus ED (with a higher inner energy window) is needed to obtain accurate interpolation (see Fig. S18d).In comparison, the pseudopotential PAOs which we regenerated with 3d orbitals (as discussed in Section II C 2) are better projectors for the wavefunctions.Indeed, the first 12 bands have projectabilities almost equal to 1, and as a consequence PD itself already provides accurate band interpolation (all the low-energy conduction states are frozen since their projectabilities are high, see Fig. S20a).Moreover, we mention that when adding additional projectors one needs to make sure that they have the correct number of radial nodes: e.g., the gold pseudopotential from SSSP [46] contains 5s + 5p semicore states, and 6s + 5d orbitals for valence electrons.If one wants to add an additional 6p orbital, it is important to ensure that the 6p orbital has one radial node, such that it is orthogonal to the nodeless 5p semicore state; Otherwise, the Bloch wavefunctions would project onto the 5p semicore state, and PD would only disentangle the 5p semicore states instead of the 6p orbitals contributing to bands above the Fermi energy.In summary, including more projectors can further improve the interpolation quality, but at the expense of increasing the number of orbitals in the model.The combination of PD and ED enables to improve the interpolation quality of low-projectability states while keeping the TB model size small.Automatic checks could be implemented in the future in the AiiDA workflows to detect whether the projectability drops below a certain threshold, and in that case either raise a warning or automatically add more projectors.

III. CONCLUSIONS
We present an automated method for the automated, robust, and reliable construction of tight-binding models based on MLWFs.The approach applies equally well to metals, insulators and semiconductors, providing in all cases atomic-like orbitals that span both the occupied states, and the empty ones whose character remains orbital-like and and not free-electron-like.The method is based on the band projectability onto pseudo-atomic orbitals to select which states are kept identically, dropped, or passed on to the established disentanglement procedure.We augment such projectability-based selection with an additional energy window to guarantee that all states around the Fermi level or the conduction band edge are well reproduced, showing that such a combination enables accurate interpolation even when minimal sets of initial atomic orbitals are chosen.This results in compact Wannier tight-binding models that provide accurate band interpolations while preserving the picture of atomic orbitals in crystals.We refer to the method collectively as projectability-disentangled Wannier functions (PDWF).
The Wannierization process is implemented as fully automated AiiDA workflows.We compare PDWFs with the other method that is also fully automated, namely SCDM.We show with a detailed study of 200 structures that PDWFs lead to more accurate band interpolations (with errors with respect to the original bands at the meV scale), and are more atom-centered and more localized than those originating from SCDM.The high accuracy in band interpolations, the target atomic orbitals obtained, and the low computational cost make PDWFs an ideal choice for automated or high-throughput Wannierization, which we demonstrate by performing the Wannierization of 21,737 non-magnetic structures from the Materials Cloud MC3D database.

IV. METHODS
We implement the PAO projection in the pw2wannier90.xexecutable inside Quantum ESPRESSO (QE) [42,51]; the PD and PD+ED methods are implemented on top of the Wannier90 code [4].In terms of the practical implementation, computing PAO projections is more efficient in both computational time and memory than the SCDM QR decomposition with column pivoting (QRCP) algorithm, since the A mnk matrices (i.e., the inner products of Bloch wavefunctions with PAOs) can be evaluated in the plane-wave G vector space, rather than requiring a Fourier transform and decomposition of very large real-space wavefunction matrices.Furthermore, since the HT Wannierization can be computationally intensive, we implement a "k−pool parallelization strategy" inside pw2wannier90.x,similarly to the main pw.x code of QE, to efficiently utilize many-core architectures by parallelizing over "pools" of processors for the almost trivially-parallel computations at each k−point.Test results show that k−pool parallelization significantly improves the efficiency of pw2wannier90.x(benchmarks are shown in Fig. S10).
The DFT calculations are carried out using QE, with the SSSP efficiency (version 1.1, PBE functional) library [46] for pseudopotentials and its recommended energy cutoffs.The HT calculations are managed with the AiiDA infrastructure [33][34][35] which submits QE and Wannier90 calculations to remote clusters, parses, and stores the results into a database, while also orchestrating all sequences of simulations and workflows.The automated AiiDA workflows are open-source and hosted on GitHub [47].The workflows accept a crystal structure as input and provide the Wannier-interpolated band structure, the real-space MLWFs, and a number of additional quantities as output.Semicore states from pseudopotentials are automatically detected and excluded from the Wannierizations, except for a few cases where some semicore states overlap with valence states; in such cases, all the semicore states are Wannierized, otherwise the band interpolation quality would be degraded, especially for SCDM.A regular k−point mesh is used for the Wannier calculations, with a k−point spacing of 0.2 Å−1 , as selected by the protocol in Vitale et al. [30].MLWFs are rendered with VESTA [52].Figures are generated by matplotlib [53].
V. DATA AVAILABILITY All data generated for this work can be obtained from the Materials Cloud Archive (https://doi.org/10.24435/materialscloud:v4-e9).
The modifications to the codes mentioned above implemented for this work will become available in the next releases of Quantum ESPRESSO (pw2wannier90.x) and Wannier90.

FIG. 1 .
FIG. 1.Comparisons of graphene band structures interpolated using different methods.(a) DFT band structure, shown as grey lines.The colored dots represent the projectabilities onto carbon 2s (green) and 2p (red) orbitals.The size of each dot is proportional to the total projectability p mk of the band m at k-point k; see Eq. (5).For a detailed plot of total projectability, see Fig.S1.Comparisons of the original and the Wannier-interpolated bands for (b) hydrogenic projections with energy disentanglement (ED), (c) SCDM, (d) PAO with ED, (e) PAO with projectability disentanglement (PD), and (f) PAO with PD+ED.The Fermi energy EF (horizontal black dashed line) is at zero; the horizontal blue dashed line denotes the top of the inner energy window, where applicable.

FIG. 3 .
FIG.3.Comparisons of silicon band structures interpolated using different methods.(a) DFT band structure, shown as grey lines.The colored dots represent the projectabilities of silicon 3s (green) and 3p (red) orbitals.The size of the dot is proportional to the total projectability p mk of the band m at k-point k.For a detailed plot of total projectability, see Fig.S2.Comparisons of the original and the Wannier-interpolated bands for (b) hydrogenic projections with ED, (c) SCDM, (d) PAO with ED, (e) PAO with PD, and (f) PAO with PD+ED.The CBM (horizontal black dashed line) is at zero; the horizontal blue dashed line denotes the top of the inner energy window, i.e., CBM + 2 eV, where applicable.Note in (c), (e), and (f), the cyan lines with circle markers show the interpolated bands obtained including also 3d orbitals, and consequently increasing the dimensionality of the disentangled manifold.These additional states are beneficial because of the presence of an intrinsic d component at the bottom of the conduction manifold, and lead to more accurate band interpolations.

2 ,
193/200 = 96.5% of the structures have η 2 ≤ 20 meV, while for SCDM it is 183/200 = 91.5%.More details are listed in Table shows the histograms for d NN and d NNN for the two methods.The PDWF average d NN = 0.43 Å is smaller than the SCDM d NN = 0.53 Å, and correspondingly the PDWF d NNN = 2.19 Å is instead larger than the SCDM d NNN = 2.11 Å.

FIG. 5 .
FIG. 5. Histogram (red) and cumulative histogram (blue) of the band distances η2 and η max 2 for 200 reference structures.(a) η2 of PDWF, (b) η2 of SCDM, (c) η max 2 of PDWF, and (d) η max 2 of SCDM.The orange (green) vertical line is the mean (median) of the band distance for the 200 structures; their values are shown in the right of each panel; PDWF provides approximately an improvement by a factor of 3.

2 d 29 FIG. 6 .
FIG. 6. Histogram of the distances of the WF centers from the NN atom (red, d NN ) and NNN atom (green, d NNN ), for 200 reference structures.(a) PDWF and (b) SCDM.The inset of each panel shows the histogram of the ratio of dNN/dNNN.The numbers in the lower right of each inset are the averages over all the 6818 MLWFs; PDWF provides MLWFs that are both closer to the NN atom and further away from the NNN atom.

FIG. 7 .
FIG. 7. Histogram (red) and cumulative histogram (blue) of WF spreads for 200 reference structures.(a) PDWF and (b) SCDM.The orange (green) vertical line is the mean (median) spread of the 6818 MLWFs, their values are shown in the right of each panel.The long tail of MLWF spreads obtained with SCDM is absent in PDWF.

FIG. 8 .
FIG. 8. Exponential decay of the Hamiltonian H(R) in the basis of MLWFs.(a) Exponential-form fitting of Frobenius norm of the Hamiltonian ∥H(R)∥ w.r.t. to the 2norm of lattice vector ∥R∥ for the case of Br2Ti, for PDWF (red) and SCDM (blue).The τ reported are the fitted decay lengths of the PDWF and SCDM Hamiltonians, respectively.(b) Histogram of decay lengths τ for the 200 reference materials, obtained using PDWF (red) and SCDM (blue).The vertical lines indicate the mean τ of PDWF and SCDM, respectively.

FIG. 9 . 2 .
FIG. 9. Histogram (red) and cumulative histogram (blue) of the PDWF band distances for 21,737 nonmagnetic structures obtained from the materials cloud MC3D database [36].(a) Average band distance η2 and (b) max band distance η max 2 .The orange (green) vertical line is the mean (median) of the band distance for the 21,737 structures; their values are shown in the right of each panel.