Automated high-throughput wannierisation

Maximally-localised Wannier functions (MLWFs) are routinely used to compute from first-principles advanced materials properties that require very dense Brillouin zone integration and to build accurate tight-binding models for scale-bridging simulations. At the same time, high-thoughput (HT) computational materials design is an emergent field that promises to accelerate the reliable and cost-effective design and optimisation of new materials with target properties. The use of MLWFs in HT workflows has been hampered by the fact that generating MLWFs automatically and robustly without any user intervention and for arbitrary materials is, in general, very challenging. We address this problem directly by proposing a procedure for automatically generating MLWFs for HT frameworks. Our approach is based on the selected columns of the density matrix method and we present the details of its implementation in an AiiDA workflow. We apply our approach to a dataset of 200 bulk crystalline materials that span a wide structural and chemical space. We assess the quality of our MLWFs in terms of the accuracy of the band-structure interpolation that they provide as compared to the band-structure obtained via full first-principles calculations. Finally, we provide a downloadable virtual machine that allows to reproduce the results of this paper, including all first-principles and atomistic simulations as well as the computational workflows.


I. INTRODUCTION
The combination of modern high-performance computing, robust and scalable software for first-principles electronic structure calculations, and the development of computational workflow management platforms, has the potential to accelerate the design and discovery of materials with tailored properties using first-principles highthroughput (HT) calculations. [2][3][4][5] Wannier functions (WFs) play a key role in contemporary state-of-the-art first-principles electronic structure calculations. First, they provide a means by which to bridge lengthscales by enabling the transfer of information from the atomic scale (e.g., density-functional theory and many-body perturbation theory calculations) to mesoscopic scales at the level of functional nano-devices (e.g., tight-binding calculations with a first-principlesderived WF basis). 6,7 Second, the compact WF representation provides a means by which advanced materials properties that require very fine sampling of electronic states in the Brillouin zone (BZ) may be computed at much lower computational cost, yet without any loss of accuracy, via Wannier interpolation. 8 Among several variants of WFs, 9 maximally-localised Wannier functions (MLWFs), based on the minimisation of the Marzari-Vanderbilt quadratic spread functional Ω, are those most employed in actual calculations in the solid state. 9 One ingredient in the canonical minimisation procedure is the specification of a set of initial guesses for the MLWFs. These are typically trial functions localised in real-space that are specified by the users, based on their experience and chemical intuition. As shall be described in more detail later, in the case of an isolated manifold of bands, the final result for the MLWFs is almost always found to be independent of the choice of initial guess. 10 In the case of entangled bands, 11 however, this tends not to be the case and the choice of initial guess strongly affects the quality of the final MLWFs, presenting a challenge to the development of a general-purpose approach to generating MLWFs automatically without user intervention.
Several approaches have been put forward to remove the necessity for user-intervention in generating MLWFs, including the iterative projection method of Mustafa et al., 12 the smooth orthonormal Bloch frames of Levitt et al., 13 and the automated construction of pseudoatomic orbitals rather than WFs as the local basis to represent the target space, as described by Agapito et al.. [14][15][16] In addition, some ad hoc solutions have been proposed, whose range of applicability is focused onto specific classes of materials. [17][18][19][20] A recently proposed algorithm by Damle et al., 1,21 known as the selected columns of the density matrix (SCDM) method, has shown great promise in avoiding the need for user intervention in obtaining MLWFs. Based on QR factorisation with column pivoting (QRCP) of the reduced single-particle density matrix, SCDM can be used without the need for an initial guess, making the approach ideally suited for HT calculations. The method is robust, being based on standard linear-algebra routines rather than on iterative minimisation. Moreover, the authors have proposed an efficient algorithm for the QRCP factorisation that operates on a smaller and numerically more tractable matrix than the full density matrix. Finally, SCDM is parameter-free for an isolated set of composite bands, and requires only two parameters in the case of entangled bands together with the choice of the target dimensionality for the disentangled subspace (i.e., the number of MLWFs required). We emphasize here that the SCDM method can be seen as an extension to solid-state periodic systems of the Cholesky orbitals approach of Aquilante et al., 22 that has been developed from a quantum-chemistry molecular perspective for finite systems. SCDM focuses instead on periodic systems, and it is based on a real-space grid discretisation of the wavefunctions. We discuss in more detail this equivalence in Sec. II B 1 and Appendix E.
In this article, we present a fully-automated protocol based on the SCDM algorithm for the construction of MLWFs, in which the two free parameters are determined automatically (in our HT approach the dimensionality of the disentangled space is fixed by the total number of states used to generate the pseudopotentials in the DFT calculations). We have implemented the SCDM algorithm in the pw2wannier90 interface code between the Quantum ESPRESSO software package 23 and the Wannier90 code. 24 We have used our implementation as the basis for a complete computational workflow for obtaining MLWFs and electronic properties based on Wannier interpolation of the BZ, starting only from the specification of the initial crystal structure. We have implemented our workflow within the AiiDA 25 materials informatics platform, and we used it to perform a HT study on a dataset of 200 materials.
We anticipate here that the SCDM method works extremely well for band-structure interpolations both for insulating and metallic systems, but is less suitable for other applications where, for instance, a specific symmetry character of the WFs is required.
The manuscript is organised as follows. A summary of the background theory is presented in Sec. II, starting with MLWFs for isolated and entangled bands in Sec. II A followed by the SCDM algorithm in Sec. II B, where we focus in particular on providing a physical interpretation of the method. In Sec. III we provide a preliminary comparison, for a few well-known materials, between MLWFs obtained via the conventional method (i.e., with userdefined initial guesses) and those obtained from SCDM. Sec. IV contains the validation of the SCDM method and our workflow for the valence bands of 81 insulating materials. In Sec. V we then discuss our automated protocol to determine the free parameters in the case of entangled bands and validate it on a dataset of 200 semiconducting and metallic materials. Finally, details on the implementation of the SCDM method in pw2wannier90 and of the AiiDA workflow are presented in Sec. VI.

II. BACKGROUND THEORY
We summarise in this section the main concepts and notations related to maximally-localised Wannier functions that will be useful in the rest of the paper, following the notation of Ref. [9].
A Wannier function associated to a band n can be obtained via a unitary transformation of the Bloch state |ψ nk , known as Wannier transform 26 where V is the real-space primitive cell volume, R is a Bravais lattice vector, and the integral is over the first BZ. For clarity of notation, we assume spin-degeneracy unless otherwise specified. The gauge freedom of the Bloch state under multiplication by a k-dependent phase e iϕn(k) results in a non-uniqueness in the definition of the Wannier function. Maximally-localised Wannier functions represent the choice of gauge in which the real-space quadratic spread of the Wannier function is minimised. 9,10 In order to obtain a minimal TB basis set it is therefore beneficial to select the optimal phases that minimise the total spread, so that overlaps and Hamiltonian matrix elements between different Wannier functions decay rapidly to zero as a function of the distance between their centres. Since the integral transformation in Eq. (1) is still a unitary transformation, the resulting {|w Rn } span the same Hilbert space as the original Bloch states {|ψ nk }. Moreover, from the orthogonality of the |ψ nk readily follows the orthogonality of the |w Rn , since unitary transformations preserve inner products. Finally, two WFs |w Rn and |w R n transform into each other under translation by the Bravais lattice vector R − R . 27 A. Maximally-localised Wannier functions (MLWFs)

Isolated bands
For an isolated set of J bands describing, e.g., the valence bands of a semiconductor, the most general phase choice for a Wannier transform can be written as where U (k) is a unitary matrix that, at each wave vector k, mixes Bloch states belonging to different bands, giving as a result a set of J composite WFs. The localisation of the WFs may be improved by choosing the unitary matrices U (k) such that | ψ nk = m |ψ mk U (k) mn in Eq. (2) is as smooth as possible, i.e., analytic with respect to k (see, e.g., Duffin 28 ). Different approaches have been put forward [29][30][31][32][33] to generate well-localised WFs. In the where · n ≡ w n0 | · |w n0 and r n = r n = w n0 |r|w n0 is the centre of the n-th Wannier function. The resulting WFs are known as maximally-localised Wannier functions (MLWFs), and are the solid-state equivalent of the Foster-Boys molecular orbitals [34][35][36] in quantum chemistry.
It can be shown that 9,10 Ω I is gauge invariant, whereas Ω depends on the particular choice of the gauge (i.e., on the choice of U (k) ). For an isolated group of bands, therefore, Ω I is evaluated once and for all in the initial gauge and minimising the total spread Ω is equivalent to minimising only the gauge-dependent part Ω. For crystalline solids with translational symmetry, it is natural to work in reciprocal space, henceforth referred as k-space. Applying Blount's identities 27 for the representation of the position operator r and r 2 in k-space and discretising in k (on a uniform grid) gives 10 , and where the vectors {b} connect a BZ mesh point k to its nearest neighbours k+b, the associated weights w b come from the finite difference representation of the gradient operator in k-space (a result of the change of representation r → i/ ∇ k ), and M (k,b) is given by Since the gradient of Ω with respect to the U (k,b) mn degrees of freedom can be expressed analytically as function of the M (k,b) mn , the minimisation of the spread functional may be obtained, for instance, by steepest-descent or conjugate-gradient methods (see Refs. [9 and 10]).

Entangled bands
In many applications, the group of bands of interest are "entangled", i.e., are not separated by an energy gap from other bands throughout the whole Brillouin zone.
Souza, Marzari and Vanderbilt 11 (SMV) proposed a "disentanglement" strategy that involves two steps. In the first step, one defines an energy window that encompasses the states of interest and which contains J win k bands at each k. This defines a local Hilbert space F(k) at each k-point, which is spanned by the J win k states. Then, for a given number J ≤ min k J win k of target Wannier functions, one finds the optimal set of J-dimensional subspaces {S(k)}, with S(k) ⊆ F(k), that have maximum intrinsic smoothness over the BZ, where the intrinsic smoothness of the Hilbert space is measured by Ω I . Heuristically, Ω I represents the "change of character" of the states across the Brillouin zone. (For a rigorous derivation see Ref. [10].) The subspaces S(k) are defined as the span of {|u opt nk }, which are obtained via a unitary transformation on the |u nk that span F(k): Note that here the U dis(k) are rectangular J win k × J matrices, and are unitary in the sense that (U dis(k) ) † U dis(k) = 1 J (with 1 J being the J × J identity matrix), ensuring that {|u opt nk } form an orthonormal set.
Maximum intrinsic smoothness is achieved by choosing U dis(k) to minimise Ω I , which, as discussed earlier, is a measure of the "spillage" between neighbouring subspaces S(k). 11 In the second step, having defined a J-dimensional subspace |u opt nk at each k, one proceeds by minimising Ω following the same recipe described in the previous section for the case of an isolated manifold of bands. Further details on the disentanglement procedure can be found in Refs. [9] and [11].

Initial projections
The iterative minimisation of Ω I starts with an initial guess for the subspaces S(k). However, the spread functional is non-convex and the minimisation may get trapped in a local minimum, often resulting in complexvalued WFs 10 (in the absence of spin-orbit coupling, the WFs at the global spread minimum are expected to be real 37 ). For gradient-based minimisation methods, thus, the ability to reach the global minimum strongly depends on the choice of an appropriate starting point, sufficiently close to the final solution. To this aim, if one has a chemical intuition of the target J Wannier functions, an initial guess of J trial localised functions g n (r) can be defined. These are then projected at every k onto the J win k Bloch states inside the target energy window (for isolated bands, J win k = J, ∀ k), yielding: mn , (10) where, at every k, A (k) mn = ψ mk | g n is a J × J square matrix in the case of an isolated manifold of bands and a J win k × J rectangular matrix in the case of entangled bands. The initial unitary matrix U dis(k) can then be obtained by orthonormalising the projected guess orbitals |φ nk through a Löwdin orthogonalisation of A (k) : One possible choice, for instance, is to start from the Bloch states themselves as the projection functions (g n (r) = ψ nk (r)), so that the elements of A (k) are the (random) phases of the Bloch states that are computed by the ab initio code. In the case of isolated bands, even a poor initial choice such as this is often sufficient to reach the global minimum of the spread functional (with enough iterations of the minimisation algorithm). Conversely, in the case of entangled bands, the two-step "disentanglement" procedure is usually unable to reach the global minimum of the spread functional unless the initial trial orbitals are already quite close to the final solution. This strong dependence of the SMV minimisation algorithm on the initial trial functions, and hence on the user's intuition and intervention, has been the main obstruction in the development of fully-automated workflows for generating MLWFs for high-throughput applications.

B. The SCDM algorithm and its physical interpretation
An alternative method to the SMV approach described in II A 2 has recently been proposed by Damle, Lin and Ying 1,21 in the form of the aforementioned selected columns of the density matrix (SCDM) algorithm. The method uses a QR factorisation with column pivoting (QRCP) 38 of the single-particle density matrix (DM), to fix the gauge freedom in a single step, without the need for an iterative minimisation algorithm. In this section, we outline the core concepts of the SCDM method, focusing mainly on the aspects needed to provide a physical interpretation and facilitate its understanding. We refer to the original publications 1,21 for additional details.

SCDM for insulators sampled at Γ only
For clarity, we start by considering a system sampled at a single k-point, e.g. Γ, and so we drop the index k from the DM and other quantities; the extension to multiple k-points is given in the next subsection. We start by considering systems with a finite band-gap between the J valence bands and the conduction bands, e.g., insulators and semiconductors.
Let us first recall that P = J n=1 |ψ n ψ n | is gauge-invariant and it is a projector on the space S spanned by the J valence wavefunctions {|ψ n }. Moreover, in the insulating case, the real-space representation P (r, r ) ≡ r|P |r of the DM decays exponentially with the distance between two points r and r : P (r, r ) ∼ e −γ|r−r | . This is the well-known nearsightedness principle. [39][40][41] In particular, this means that for a given fixed r = r 0 , the function ϕ r0 (r) ≡ P (r, r = r 0 ) = dr P (r, r )δ(r − r 0 ) (13) represents the projection on the subspace S of a delta function centred at r 0 , and that this projection is an exponentially-localised orbital.
To understand the numerical implementation of the method, we consider from now on the real-space discretised version of the DM. The J valence wavefunctions (or, in the case of periodic systems, the periodic part u nk (r) of the J valence Bloch states) can be stored on a grid of n G points in real space r 1 , r 2 , . . . , r n G . We can then define the following n G × J matrix Ψ that contains the values of the J wavefunctions on the grid points: With this definition, the orthonormality condition is written as Ψ † Ψ = 1 J , while the density matrix (which in discretised form is an n G × n G matrix) can be written as P = ΨΨ † , i.e., P ij = J n=1 ψ n (r i )ψ * n (r j ). We can now interpret the j-th column C j of the DM, C j i ≡ P ij , as the projection on the valence subspace S of a test orbital φ j that is zero everywhere except at the j-th grid position (i.e., at position r j ). This statement is the discretised version of the projection of a delta function in Eq. (13), i.e., apart from normalisation, φ j is the discretised version of δ(r − r j ). Therefore, thanks to the near-sightedness principle, the orbitals represented by the columns of the DM are localised.
This statement is at the core of the SCDM method. In fact, when searching for Wannier functions, we are looking for a complete and orthogonal basis set of J localised functions that span the subspace S. In our case, the set of all columns C j clearly spans the whole subspace S (since the P operator is the projector on S). However, in essentially all practical situations, J n G and the set of all these n G orbitals is redundant. In addition, these orbitals are not orthogonal-intuitively, projecting on delta functions centred at two neighbouring points will typically result in a large overlap between the projected orbitals-and not normalised (e.g., in the limiting case of a delta function centred at a position in space where there is no charge density, the resulting projection will have zero norm). Selecting any set of J linearly-independent columns would form a basis for S, and an initial guess for the Wannier functions could be obtained by orthonormalising these J columns, e.g., with a Löwdin symmetric orthogonalisation. However, if these J columns are not already almost orthogonal, the orthogonalisation will be numerically unstable and, most importantly, will mix them and thereby degrade their localisation. Therefore, the goal of the SCDM method is to select the "most representative" J columns, i.e., the columns that possess the largest norm and that are as orthogonal to each other as possible, i.e. the most "well-conditioned subset", so that the Löwdin orthogonalisation will mix these orbitals as little as possible (Löwdin orthogonalisation minimises the squared difference between the original and orthogonalised functions 42 ). Equivalently, as every column is the projection of a delta-like test orbital centred at r i , we can say that the SCDM algorithm selects J points, from among the original n G grid points, that define the "most representative" localised projected orbitals. To achieve this goal, SCDM uses the standard linear algebra QRCP method, 38 which factorises a matrix P as P Π = QR, where Q is a matrix with orthonormal columns, R is a upper-triangular matrix, and Π is a permutation matrix that swaps the columns of P so that the the diagonal elements of R are in order of decreasing magnitude |R 11 | ≥ |R 22 | ≥ · · · ≥ |R n G n G | (see Appendix B for more details). The relevant output of the algorithm is the Π permutation matrix, or more specifically the indexes of the first J columns chosen by the algorithm: these are the "most representative" columns discussed above and, after orthonormalisation, they provide the best guess for the localised Wannier functions of the system. With a slight abuse of notation, in the following we will use the symbol Π also to identify the vector of indexes of the permutation matrix, such that Π(i) = j has the following meaning: Π ij = 1, and all of the other elements in the j-th column are equal to zero.
QRCP (a greedy algorithm) selects columns as follows: since R is triangular (and Q has orthonormal columns), the norm of the first selected column C Π(1) of P is |R 11 | 2 and must be the largest possible, therefore the algorithm will choose the column with the largest norm. The second column C Π(2) is chosen to maximise |R 22 | 2 that, due to the properties of Q and R, is the component of C Π(2) orthogonal to C Π(1) , as shown in Appendix D. So, the QRCP algorithm will select as the second vector the one with the largest orthogonal component to the first, and in general will select the k-th vector as the one with the largest orthogonal component to the subspace spanned by the previous (k − 1) columns (to be more precise the actual selection process is a heuristic for trying to keep principal sub-matrices of R as wellconditioned as possible). It is worth mentioning that this approach is related to the Cholesky orbitals approach of Aquilante et al., 22 that applies to finite (non-periodic) systems and for a different basis set (a basis of atomic orbitals rather than a real-space grid discretisation). In particular, the Cholesky algorithm used in Ref. [22] is a refined version of the original Cholesky decomposition specifically adapted for positive semi-definite matrices, i.e., Cholesky decomposition with full column pivoting (CholCP) Π T P Π = L † L, where L is an upper triangular matrix and Π is a permutation matrix. In Appendix E we demonstrate that the selection of the columns in CholCP is the same as in QRCP, at least for the first J = rank(P ) columns, i.e., (P Π) :,1:J = (PΠ) :,1:J . This is due to well-known connections between QR factorizations and Cholesky factorizations. 38 Finally, the two methods use undoubtedly related ideas but they are not direct analogues since there are multiple "variants" of SCDM when using localised orbitals.
For an effective practical implementation of the method, a final step is required. In fact, the P matrix can be extremely large, since n G can be of the order of 100 000 or more (while J is often of the order of 10-100). Therefore, applying the QRCP algorithm directly to P is impractical, both for the memory required to store it (O(n 2 G )), and for the time needed to compute the result (O(J × n 2 G )). Instead, using the fact that P = ΨΨ † and that the original columns of Ψ are orthonormal, one can prove (see Appendix C) that the same permutation matrix Π can be obtained applying the QRCP algorithm directly to the much smaller matrix Ψ † (of size J × n G ), with a computational cost that scales as O(J 2 × n G ). Moreover, the matrix obtained from the first J columns of (Ψ † Π) may be used as the A mn projection matrix of Eq. (10) as a starting point for the usual Wannierisation procedure in order to obtain MLWFs.
Finally, it is worth noting the connection with the "canonical" approach of user-defined initial guesses (e.g., atomic-like orbitals at specified centres): the SCDM method may be thought of as using as initial guesses a set of extremely localised s-like "orbitals" (actually, δ functions), whose centres (located at the points of the real-space grid) are optimally chosen by the SCDM algorithm via the QRCP factorisation.

SCDM for periodic systems: SCDM-k
We now extend the discussion to the case of k-point sampling with more than one k-point (i.e., not only at Γ), still considering an isolated manifold (e.g., the valence bands). The DM P (k) = k P k = n,k |ψ nk ψ nk | is an analytic function of k, 37,43 and it is also proven that WFs with an exponential decay exist; 44 numerical studies for the specific case of MLWFs have confirmed this claim for several materials, 44,45 and recently there has been a formal proof for 2D and 3D time-reversal-invariant insulators. 37 The SCDM method has been extended also to the case of k-sampling 1 and named in this case "SCDM-k". In summary, the goal is now to select a common set of columns for all the k-dependent density matrices P k . Ref. [1] discusses extensively how the method can be extended to a k-point sampling with more than one k point and it shows detailed results of the convergence as a function of the number of k points used in the columnselection algorithm. The final conclusion of the authors is that it is typically sufficient to select the columns using a single "anchor" k point (typically chosen to be Γ), i.e., it is sufficient to compute the permutation matrix Π using a QRCP on P k=Γ only. Then, this selection of columns can be used for all other k-points.

Extension to entangled bands
Finally, the extension to the entangled case (e.g., for metals or when considering also the conduction bands of insulators and semiconductors) has been proposed in Ref. [1]. In this case, a so-called quasi-density matrix is defined, where f ( nk ) is an occupancy function. The isolatedbands case can be recovered by setting f ( nk ) = 1 for energy values nk within the energy range of the isolated bands, and zero elsewhere. For the typical cases of interest of this work (metals, and valence bands and low-energy conduction bands in semiconductors and insulators), one needs bands up to a given energy (typically slightly above the Fermi energy). Then, as suggested in Ref. [1], f ( ) can be chosen as the complementary error function: This function depends on two free parameters µ and σ, whose choice is critical to tune the algorithm and obtain a set of Wannier functions that correctly interpolate the low-energy electronic bands of a given material. In Sec. V A we describe our protocol to choose the values of µ and σ based on the electronic structure of the material, allowing us to implement a fully automated workflow to construct its Wannier functions via the SCDM method. The algorithm then proceeds as in the case for isolated bands, computing the QRCP factorisation on the quasi-density-matrix or, in practice, on the matrix F k Ψ † k at the k = Γ anchor point, with F k a diagonal matrix with matrix elements {f ( 1,k ), . . . , f ( J win k ,k )}. This approach, therefore, constitutes an alternative to the SMV disentanglement procedure described in Sec. II A 2: the matrices obtained from the first J selected columns of F k Ψ † k at each k form the projection matrices A (k) , and the U dis(k) matrices of Eq. (9) are obtained using the Löwdin transformation of Eq. (11).

SCDM and MLWFs
The SCDM algorithm is able to robustly generate welllocalised functions without the need for an initial guess. Whilst this makes the algorithm well-suited for direct integration within HT frameworks, the selection of the columns cannot be controlled by external parameters (at least for isolated bands), and therefore it is not possible to enforce constraints that might be desirable, such as point symmetries. On the contrary, when explicitly specifying atomic-like initial projections, these (if appopriately chosen) provide at least some degree of chemical and symmetry information. In Sec. III we discuss how this affects the WFs obtained by the algorithm. Our aim is to leverage on the ability of SCDM to automatically generate a good set of localised functions, and to use these to seed the MV algorithm for the minimisation of the total spread functional, which will give in turn an automated protocol to generate MLWFs. Being able to automatically generate MLWFs will also allow users to seamlessly exploit the set of computational tools that have been developed in recent years for MLWFs and implemented in various codes, such as Wannier90. In practice, this entails employing the SCDM algorithm to compute the A (k) matrices of Eq. 10 as follows: where the J points r n are obtained from the first J columns of the permutation matrix Π, computed at Γ, i.e., representing the reduced matrix formed by the first J columns of Π Γ .

SCDM and "disentanglement"
It is worth noting that the SCDM method can be also combined with the SMV disentanglement procedure of Sec. II A 2, as a means of seeding the initial subspace projection. However, this introduces two additional parameters associated with the SMV approach, namely ε outer , and ε inner , giving a total of four parameters (together with µ and σ). ε outer defines the upper limit of the socalled "outer" energy window discussed in Sec. II A 2, and ε inner defines the upper limit of a smaller energy window contained within the outer energy window. This inner window is used to "freeze" the Bloch states within during the minimisation of Ω I , such that they are fully preserved within the selected subspaces {S(k)} (see Ref. [11] for a comprehensive description of the outer and inner energy windows). Each additional parameter makes it increasingly difficult to find a robust and automated protocol for obtaining MLWFs. Consequently, when combining SCDM with SMV disentanglement, an optimal selection of all the parameters can be achieved only in an ad hoc, non-automatic fashion (hence only for few materials). As shown in Sec. II B 3, SCDM employs a generalised form of the density matrix Eq. (15), which implicitly defines an energy window via the function f (ε) and selects a smooth manifold by construction. Intuitively, this suggests that SCDM can be used in lieu of the SMV disentanglement procedure. In general, we have found that for the sole purpose of interpolating the energy bands up to a given energy, performing SMV disentanglement step on top of SCDM has at best a marginal improvement on the quality of the interpolation (see Sec. V), and in some cases can even be detrimental due to the case-by-case sensitivity on the choice of energy windows. For this reason, in Sec. V we focus exclusively on a protocol for the automatic selection of the free parameters in SCDM, i.e., µ and σ, without considering any additional SMV disentanglement.

III. SCDM VS MLWFS IN WELL-KNOWN MATERIALS
As a precursor to the fully-automated high-throughput study on a set of 200 materials that focuses on automatic Wannierisation and band interpolation from SCDM projections and which will be presented in Sec. V, in this section we consider in greater depth and detail the performance of the SCDM method on a small set of simple systems with well-known Wannier representations of the electronic structure. Specifically, we compare quadratic spreads, centres and symmetries of the WFs computed from the SCDM gauge (as described in Sec. II B) with the ones computed from carefully chosen initial projections. Comparative studies between SCDM localised functions and MLWFs on well-known materials have recently appeared in the literature. 1,46 However, here we expand on different aspects, focusing in particular on the combination of the SCDM and the MV approaches (SCDM+MLWFs), to better assess its range of applicability, for instance for beyond-DFT methods, e.g., ab initio tight-binding, 47,48 DFT+U 49-51 and DMFT, 52,53 where the symmetries of the Wannier functions are important.
All DFT calculations have been carried out with Quantum ESPRESSO, using the PBE exchangecorrelation functional and Vanderbilt ultrasoft pseudopotentials. 54 MLWFs are generated from Bloch states calculated on a 10 × 10 × 10 Monkhorst-Pack grid of k-points. The SCDM method has been implemented in the pw2wannier90 code, which interfaces Quantum ESPRESSO with the Wannier90 code, 24,55 as explained in Sec. VI. Wannier90 is used throughout this work to generate the WFs on a real-space grid and to perform the interpolation of band structures in reciprocal space.
We consider four different schemes for generating Wan-nier functions: (1) Full minimisation of Ω using the SMV disentanglement algorithm to minimise Ω I and the MV algorithm to minimise Ω (DIS+MLWF); (2) Minimisation of Ω I only, using the SMV algorithm (DIS); (3) Minimisation of Ω only, using the MV algorithm (MLWF); and (4) No minimisation of Ω (proj-ONLY). In each case, the initial J-dimensional subspace at each k is determined in one of two ways, either by the SCDM method or by projection onto specific atomic-like localised orbitals (Eq. (10)).

A. Silicon
We start by studying the Wannierisation of a manifold of bands consisting of the four valence bands plus the four low-lying conduction bands in silicon, the latter being entangled with bands at higher energies. For the SCDM method, we use σ = 2 eV and µ = 10 eV. This choice is equivalent to that of Ref. [1], taking into account a shift in the absolute energy scale, which shifts the value of µ. The outer and inner energy windows (described in Sec. II A 2), obtained through convergence tests, are set to ε outer = 17.0 eV and ε inner = 6.5 eV.
When using initial projections onto atomic-like orbitals, we find that the spread functional Ω has three minima that are very close to each other and each of which gives eight real MLWFs. The global minimum corresponds to four sp 3 -type MLWFs per Si atom in the twoatom unit cell, oriented in a back-bonding (BB) configuration, i.e., with the major lobes of the sp 3 -type MLWFs pointing towards the tetrahedral interstitial sites. A representative example of one such BB MLWF is shown in the isosurface plots in the first row of Fig. 1. Intuitively, from an atomic orbital perspective, one might instead expect the sp 3 -type MLWFs to be in a front-bonding (FB) configuration, i.e., with the major lobes pointing towards the vertices of the tetrahedra centred on the two non-equivalent Si atoms, as shown in the isosurface plots in the second row of Fig. 1. However, this FB configuration corresponds to a slightly larger value of the total spread Ω and, therefore, constitutes a local minimum of the spread. A third (intermediate) local minimum gives four sp 3 -type MLWFs that are in the BB configuration on one Si atom in the unit cell and four sp 3 -type in the FB configuration on the other Si atom. At variance with what is stated in Ref. [46], all these cases can be found by specifying as initial projections four appropriately oriented sp 3 -type orbitals on each Si atom in the unit cell. 56 With these initial projections, the four different minimisation options described earlier give the same qualitative results. Going from the DIS+MLWF case to DIS to MLWF to proj-ONLY, the spreads of the MLWFs increase, as expected, but the FB/BB character is consistently present (see the top two rows of Fig. 1 Figure 1: WFs obtained by wannierising the four valence bands plus the four low-lying conduction bands in silicon with the four different options described in the text. First row: the initial subspace is defined by projecting the Bloch states ψ nk (r) on eight appropriately oriented sp 3 -type orbitals giving back-bonding (BB) MLWFs in all cases.
Second row: as above but with different orientations for the sp 3 -type orbitals, resulting in front-bonding (FB) MLWFs in all cases. Third row: the initial subspace is obtained from the SCDM method. Here, the eight sp 3 -type WFs are in the BB configuration only when a full minimisation is performed. In all other cases a mixture of configurations is obtained instead. The values below each WF isosurface (isovalue=±0.45Å −3/2 ) is the value of the individual spread inÅ 2 .
26.54Å 2 to 20.06Å 2 in both the FB and BB cases, showing that the initial and final selected subspaces from the two different choices of projection have the same intrinsic smoothness. Instead, starting from SCDM to define the initial subspace, we obtain different qualitative results for the four different minimisation schemes. Wannier functions in the BB configuration are found when a full minimisation is performed (i.e., SCDM followed by SMV and MV minimisation). A representative example of one such WF is shown in the third row and first column of Fig. 1. SCDM selects a less smooth initial subspace (Ω I = 27.54Å 2 ) than specifying atomic orbital initial projections (26.54Å 2 ), but the final spreads are the same as in the equivalent BB case with atomic orbital initial projections. We also observed that in the case of SCDM, the minimisation of both Ω I and Ω required more it-erations to achieve the same level of convergence, perhaps reflecting the fact that the initial subspace is less smooth. When using the other minimisation schemes, we find functions of both FB and BB character, all with slightly different individual spreads. Representative isosurfaces are shown in the last three columns of the row labelled "SCDM" in Fig. 1. It is clear that the tetrahedral site symmetry is not preserved in the resulting WFs. Moreover, there is no clear pattern in the individual spreads going from the DIS case to the proj-ONLY case.
When looking at the interpolated band structure, however, a different picture emerges. In the case of choosing atomic orbital projections, the interpolation is very poor if no SMV disentanglement step is included in the minimisation. This shows the importance of disentangling the correct manifold and it is in agreement with what has been previously reported in the literature. 9 On the other hand, in the case of an SCDM-generated initial subspace, the interpolation is only marginally affected by the minimisation scheme employed (see Fig. S1 in the Supplementary Materials).
To summarise, in silicon SCDM performs very well when combined with full spread minimisation, both in terms of the symmetries of the WFs and band interpolation (see Fig. S1). When SCDM is used in isolation, the individual spreads of the resulting WFs are larger than WFs generated from user-defined atomic orbital projections; the quality of band structure interpolation, however, is almost independent of whether or not subsequent spread minimisation is carried out.

B. Copper
Copper presents a paradigmatic case of a transition metal where a set of bands (e.g., of d-orbital character) cross and mix in a narrow energy window around the Fermi energy with a set of broad, nearly-free-electron bands. In this case, the SMV algorithm turns out to be very sensitive to the choice of the initial gauge and a good Wannier representation of the band structure can be achieved only by a careful choice of both initial projections and energy windows. Consequently, the possibility of bypassing these user-intensive steps makes the SCDM an attractive approach. This is particularly important for methodologies such as ab initio tight binding, 48 DFT+U 50 and DMFT, 53 which deal with strong correlation in a local subspace, e.g., the subspace spanned by d orbitals (for transition metals or transition-metal oxides) or f orbitals (for rare-earth or actinide intermetallics). For copper, as suggested by Souza et al., 11 in order to generate a faithful representation of the band structure around the Fermi level, we work with a manifold of dimension J = 7, which contains one more function than the conventional minimal basis usually employed in tightbinding models. For this system, we focus only on the full minimisation scheme (DIS+MLWF), as it is the most rep-resentative when comparing the symmetries of the WFs, as shown in the previous section. For the disentanglement step we set ε outer = 38.0 eV and ε inner = 19.0 eV. For SCDM, we set µ = 11.40 eV and σ = 2.0 eV. The Fermi energy in our calculation is at 12.18 eV. As shown in Ref. [11], appropriately selected initial projections are five d-type orbitals centred on the Cu atom and two stype orbitals, each centred on one of the two tetrahedral interstitial sites. The resulting seven MLWFs respect the symmetries one would expect from group theory. In fact, the five d-like functions give a representation of dimension 3+2 of the O h point group (which is isomorphic to the site-symmetry group of the origin), with the usual t 2g and e g character (see Fig. 2(a) and Fig. 2(b)). The two s-like functions give each a one-dimensional representation (a 1 ) of T d (which is the site-symmetry group of the tetrahedral interstitial sites), as shown in Fig. 2(c).
When using SCDM projections, the symmetries of the d-type MLWFs are not fully recovered. This can clearly be seen in Figs. 2(d) and 2(e), where the d-type functions show mixed t 2g /e g character (this is a feature of all five d-type functions).

IV. ISOLATED BANDS
Until here, we have looked into the details of the Wannier functions that can be obtained from SCDM projections, by focusing on the paradigmatic examples of silicon and copper (Sec. III). We focussed on comparing Wannier functions as obtained by adopting different initial projections, given that good atomic-like projections can often be easily identified through chemical intuition. Now we take a complementary perspective, by considering any given crystal structure, where we face the problem of finding good initial projections without any prior chemical knowledge of the system. This is particularly relevant for high-throughput studies, where crystal-structure databases are systematically screened with first-principles simulations. In order to produce high-throughput Wannier functions, it is fundamental to provide an algorithm that does not require human interaction in the choice of the initial projections. In addition, such an algorithm must be able to use only information that is either contained in the crystal structure and the pseudopotential, or that can be computed by a simple first-principles simulation, such as the projected density of states. To this aim, human-specified atomic-like projections are not suitable, and we propose the SCDM method as the workhorse for the automated choice of the initial projections.
In order to ascertain the effectiveness of the SCDM method in generating well-localised Wannier functions in an automated way, we start by testing the algorithm for isolated manifolds. We compare Wannier interpolations and direct DFT calculations for the band structure of the valence bands of a set of 81 insulating bulk crystalline materials spanning a wide range of chemical and Figure 2: MLWFs obtained by performing a full minimisation of Ω for the s-d complex, formed by seven bands (J = 7) around the Fermi level in copper. First column: three representative MLWFs obtained from using atomic orbital projections to define the initial subspace (see main text for description). Panel (a) shows one of the three MLWFs with t 2g character; panel (b) shows one of the two MLWFs with e g character; panel (c) shows one of the two broad s-like orbitals centred on a tetrahedral-interstitial site. Second column: three representative MLWFs obtained from using SCDM to define the initial subspace. Panel (d) and (e) show two of the five MLWFs with mixed t 2g /e g character; panel (f) shows one of the two broad s-like orbitals centred on an tetrahedral-interstitial site. Below each function its individual spread inÅ 2 is reported. Isosurfaces are plotted with an isovalue of ±0.45Å −3/2 . structural space, for the full list the Reader is referred to Ref. [57]. We quantify the differences between two band structures by introducing a simple metric that is inspired by the so-called "bands distance" introduced in Ref. [58]. Here we define the distance between DFT and Wannier-interpolated bands as: where ε DFT nk and ε Wan nk are respectively the DFT and Wannier-interpolated band structures, and the summation runs over the occupied bands only. Later in Sec. V, we will introduce a finite smearing to deal with conduction-band states and metallic systems. As in Ref. [58], to take into account the possibility that significant differences between band structures may occur only in sub-regions of the Brillouin zone or in small energy ranges, we also compute (19) where, essentially, we select the point (nk) with the worst interpolation, which is responsible for the largest contribution to η. We use η and η max to assess the effect of iteratively minimising the spread Ω to obtain maximally-localised Wannier functions ("SCDM+MLWF"), compared to the one-shot Wannier orbitals that are obtained by using the SCDM projections only ("SCDM-only"). We note that in the following MLWF might refer either to a maximally-localised WF or to the maximal localisation procedure itself, the meaning being always clear from the context.
For each of the 81 structures of the benchmark set, we first perform a variable-cell optimisation and we then compute the band structure on a high-symmetry path using DFT. The cell and the path are standardised using seekpath according to the prescription of Ref. [59]. The ground-state charge density is obtained using a kpoint density of 0.2Å −1 in the irreducible Brillouin zone (unless otherwise stated). Band structures are then calculated using the charge density frozen from the earlier calculation and sampling the high-symmetry path with a density of 0.01Å −1 . Then we compute the WFs and the real-space Hamiltonian with Wannier90, starting from a non-self-consistent field (NSCF) DFT calculation performed on a possibly different k-point grid on the full BZ and employing the ground-state charge density computed earlier. At this point, the bands distance is then calculated by diagonalising the Wannier Hamiltonian using the TBmodels code 60 on the same k-points used in the DFT bands calculation.
All DFT calculations are carried out using the Quantum ESPRESSO distribution, 23 employing the PBE functional 61 and a beta version 62 of the SSSP v1.0 efficiency pseudopotential library. 58,63-67 In Fig. 3 we report histograms of η and η max for four different k-point densities, namely ρ k = 0.15, 0.2, 0.3 and 0.4Å −1 , used in the NSCF step to construct Wannier functions. We stress that for an isolated set of bands, such as for the valence bands of an insulator, the SCDM method involves no free parameters and the only parameter to set is the density ρ k of a uniform k-point grid that is used to diagonalise the Hamiltonian. Hence it is fundamental to elaborate a strategy for the choice of ρ k , as this finally removes every free parameter from the construction of Wannier functions for isolated bands.
The SCDM method is found to work well for all of the 81 systems studied, with the exception of two that have very poor interpolation. Notably, these two structures (three if we consider the SCDM-only method) are the ones that exhibit the highest initial spread Ω per Wannier function. Although a large initial spread does not necessarily imply poor interpolation, it certainly correlates with a potential risk of poor Wannierisation and it could be used as a marker for triggering a check on the quality of bands interpolation within the calculation workflow.
Those systems with η > 20 meV or, in other words, interpolated bands that are significantly less accurate with respect to the majority of the sample, are considered to be outliers. In Table I, we report the number of the outliers for the four different k-point densities, both in the case of SCDM-only and SCDM+MLWF. Clearly, increasing the k-point density produces fewer outliers and, in this respect, the SCDM+MLWF seems to converge slightly faster than SCDM-only, in agreement with the results shown in Fig. 3. As we will discuss shortly, the superior performance of SCDM+MLWF is linked with the increased localisation associated with the MLWF procedure. As mentioned before, localisation is also related to the poor interpolation of the outliers: at all k-point densities, outliers are among the systems with the largest initial spreads. On one hand, a larger initial spread signals a potential problem with the SCDM projections, on the other hand it requires a denser k-point density for convergence (the less localised the Wannier functions are, the more long-range the Wannier Hamiltonian is). Fig. 3 also shows that, when considering valence bands only, the MLWF procedure moderately improves the quality of band interpolation with respect to SCDMonly, resulting in narrower η and η max distributions, although band interpolation is often already excellent using an SCDM-only approach. We emphasise, however, that it is known that for the valence bands of gapped systems, a set of randomly-centred Gaussian functions can be often used as starting projections leading to good MLWFs. We compare, therefore, the performance of SCDM projections versus randomly-centred Gaussian orbital projections as a starting point for the MLWF procedure (which we refer to as the "random+MLWF" scheme), assessing their comparative robustness and accuracy of band interpolation. Fig. 5 reports the distribution of η and η max at the density ρ k = 0.2Å −1 .
We now elaborate on the differences between random and SCDM initial projections. First, random projections typically generate a much higher initial spread (7.5Å 2 per WF) compared to SCDM (1.0Å 2 per WF). We find that the MLWF procedure is often sufficient to localise Wannier functions even in the case of large initial spreads: for 63 out of 81 materials the MLWF procedure brings both the random projections and the SCDM projections cases to the same minimum spread value. Notably, it never happens that the spread is similar and the quality of the interpolation is very different, while the opposite happens only in the case of He, a pathological case (1 atom and 2 electrons per cell) where random projections give a poorly localised Wannier function while still being able to provide a very good interpolation. For 15 materials (16 if we include He), random projections provide a very poor starting point and the MLWF procedure remains trapped in a local minimum with large spread. In these cases, instead, SCDM projections are a good starting point with low spread and the MLWF procedure further reduces it and a higher-quality interpolation is achieved, as demonstrated by the lower η values. Finally, there are two materials for which both SCDM-only and SCDM+MLWF do not perform well, but where random+MLWF happens to perform better than SCDM+MLWF. For one of these cases, Al 2 Os, we have checked that excluding the semi-core states greatly improves the performance and the quality of the interpolated bands. We believe that the reason lies in the fact that, if semi-core states are present, then there are some projections, centred on the same site, that possess the same symmetry character, e.g., p-like projections with different principal quantum numbers (for instance 1pand 2p-like). With a relatively low plane-wave energy cutoff, the real-space grid is too coarse and there are not enough degrees of freedom for the column selection in the QRCP step to distinguish or describe sufficiently well these same-symmetry-character states.
In the other case, Se 2 Sn, there are no semi-core states. Here instead, some SCDM projections show an initial value of Ω D -the sum of the diagonal elements of Ω in (5)-that is not zero or very close to zero (Ω D > 0.5Å 2 ), which could be used as a diagnostic indicator for problematic systems. In particular, SCDM+MLWF seems to get trapped in a state in which there are a number of well-localised WFs and two that are diffuse and spread over multiple sites. This set of WF are real with a total spread of 28Å 2 and Ω D of 2Å 2 . We found that a possible solution to recover a good interpolation is to add some noise (adding small random numbers to the search direction components, as implemented in Wannier90) during the minimisation to help the algorithm escape from the unwanted local minimum.
We propose some technical solutions that could be easily added to a workflow: • Automatically detect and exclude semi-core states (if any). This is generally a safe choice as these states are not physically interesting for most applications. Alternatively, one could retain the semicore states and increase the cutoff energy (or equivalently the density of the real-space grid).
• If the problem is not in describing semi-core states, then check the value of Ω D , if it is above a given threshold (e.g., > 1.0Å 2 ) for one or more initial projections, introduce some noise in the minimisation.
• If none of the above work, switch to ran-dom+MLWF projections, which may give a better final result.
To study now more in detail the effect of minimising the spread, we start by comparing the total spread Ω obtained using SCDM+MLWF and SCDM-only, by computing: where Ω SCDM and Ω MLWF are the total spreads obtained with SCDM-only and SCDM+MLWF, respectively. As reported in Fig. 6, the SCDM-only Wannier functions are already well localised and ∆Ω Ω MLWF is less than 10% for 68% (55/81) of systems, and less than 20% for 88% of them (71/81).
An interesting question is whether the difference in spread due to the MLWF procedure correlates with the difference in the quality of the interpolation. To assess this, we compute the quantity  Figure 6: Histogram of the relative variation of the total quadratic spread Ω before and after the MLWF procedure for the valence bands of our set of 81 insulators, obtained for ρ k = 0.2Å −1 . The SCDM+MLWF procedure provides Wannier functions that are moderately more localised with respect to SCDM-only, with a relative variation within 10 − 20% for most materials.
where η SCDM and η MLWF are the band distances obtained with SCDM-only and SCDM+MLWF respectively. Fig. 7 shows a scatter plot of ∆η vs. ∆Ω/Ω MLWF , showing that a reduction in the spread typically implies an improvement in the quality of the interpolation (∆η < 0). These findings highlight that SCDM-only Wannier functions are already sufficiently localised and represent well the valence manifold, and the subsequent MLWF procedure (starting from a very good guess) safely refines the initial choice of SCDM, improving the accuracy of the Wannier Hamiltonian by increasing localisation. In general, the greatest benefit from the MLWF procedure is visible in the interpolation of the almost-flat semi-core states. In fact often, when using SCDM-only Wannier functions for the interpolation of these states, the interpolated bands show an oscillatory behaviour, with the maximum absolute difference with respect to the DFT bands of the order of a few meV (comparable to the spread of those bands). From our results, a smoother and more accurate interpolation is usually recovered after a MLWF procedure. Before discussing the case of entangled bands, we summarise here the main conclusions that can be drawn for isolated bands. All the results we obtained, displayed in interpolation is very high for 98% of the structures, with only 2 (out of 81) cases showing a poor interpolation. Although SCDM-only Wannier functions are shown to provide already accurate band structures, the MLWF procedure appears to improve both the quality of interpolation (lower η) and localisation (lower spread). Hence, we suggest the SCDM+MLWF method with ρ k = 0.2Å −1 as the standard protocol for producing accurate and efficient Wannier Hamiltonians describing the valence bands of bulk insulating crystals.

V. ENTANGLED BANDS
We now consider the case of entangled bands. With the intent of describing a fully automatic protocol, we limit ourselves to the case of Wannier interpolation of all states up to a given energy (excluding, if appropriate, manifolds of low-lying semicore states that are isolated in energy from the rest of the band structure) and we do not consider the case of computing Wannier functions for a manifold of bands of given symmetry within a narrow energy window (e.g., d states in copper or t 2g /e g states in a transition-metal oxide, see Sec. III B) that is entangled with bands above and below in energy.
In the case of entangled bands, the SCDM method demands the choice of three free parameters: µ and σ, as described at the end of Sec. II B, as well as J, the target number of Wannier functions. These parameters play a fundamental role in the selection of the columns of the quasi-DM and hence greatly affect the overall quality of the subspace selection and, consequently, the bands interpolation. In particular, since there is no equivalent definition of an inner energy window 11 in the SCDM method, it is not guaranteed that a subspace that includes the physically-relevant lowest-lying bands will be selected because the greedy QRCP algorithm, owing to an inappropriate choice of µ and σ, might favour states that are higher in energy. It is, therefore, key to the success of the automation process to have a protocol that automatically chooses these parameters in a robust and systematic way. We will now describe such a protocol, and in Sec. V B we show its effectiveness on a large set of chemically diverse materials.

A. Protocol
To identify appropriate values of µ, σ and J, we first compute the "projectability" p nk , which measures how well each Bloch state |ψ nk is represented in a Hilbert space A defined by a given set of localised functions. Indeed, in the entangled case, WFs contain contributions from the valence states plus specific conduction states, typically corresponding to the anti-bonding partners of the valence states. The selection of these specific conduction states-out of the very many-can be challenging, because they are not necessarily the lowest energy ones. This idea motivates the use of projectability as a measure to see which conduction states might be more important.
Similarly to Agapito et al., 15 we choose as our localised functions the set of N PAO pseudo-atomic orbitals (PAO) φ Ilm (r) employed in the generation of the pseudopotentials, where I is an index running over the atoms in the cell and lm define the usual angular momentum quantum numbers. We then construct Bloch sums φ µk (r) = 1 Nµ R e −ik·R φ µ (r − R), where µ = {Ilm} and N µ is the number of lattice vectors R contained in the Born-von Karman cell (which is equal to the number of k-points sampled in the BZ). Finally, a Hilbert space A k at each k-point in the BZ is defined as the space spanned by the Löwdin-orthogonalised functions φ µk (r) = ν (S k −1/2 ) µν φ νk (r), with S k µν = φ µk (r)|φ νk (r) , and A is given by the direct sum A = k A k . The projectability of each Bloch state onto A is then defined as where 0 ≤ p nk ≤ 1. The projections ψ nk |φ k Ilm are computed straightforwardly using the projwfc.x code from Quantum ESPRESSO. In particular, for the pseudopotentials considered in this work, the number of valence electrons and the atomic orbitals included in the pseudopotential files may be found in Table S1 of the Supplementary Materials.
As the first step of our protocol, we choose J as the total number of projections N PAO considered in the sum of Eq. (22). Since we aim to interpolate the bands up to a given energy above the Fermi level, fixing J = N PAO is a conservative choice, as the number of PAOs is usually greater or equal to the number of valence bands plus few conduction bands. In addition, we stress that the bands that correspond to localised states, within the full band-structure of a solid, are those that naturally emerge from bonding and anti-bonding combinations of atomic orbitals, with higher energy states having more free-electron character and thus less amenable to Wannierisation.
We then use the values of the projectability to inform the choice of µ and σ. First, we plot the projectability for all Bloch states as a function of the corresponding band energy nk , as shown in Fig. 8 (to illustrate the procedure, we show plots for one prototypical material, namely crystalline tungsten (W), but similar plots and trends also hold for the other materials considered in this work). The general trend is that p nk ∼ 1 for low-energy states, which are well-represented by the chosen pseudo-atomic orbitals, and p nk ∼ 0 for high-energy states that originate either from free-electron-like states or from localised states with an orbital character that is not included in the set of Supplementary Table S1, e.g., atomic orbitals with principal quantum number n > 3 (i.e., more than two radial nodes). We then fit this plot to a complementary error function as in Eq. (16), extracting the two parameters µ fit and σ fit . The core of our protocol lies on the actual choice of the µ and σ parameters used as input for the SCDM method by setting Let us now motivate this choice. We observe that σ fit measures the typical energy spread of the bands originating from states within A, and therefore is a good physical guess also for σ. The naive choice µ = µ fit , however, produces extremely poor interpolation of the bands for most of the materials that we have tested in Sec. V B. The reason is that it gives too great a weight in Eq. (15) to states that have relatively small projectability (p nk < 1). As a consequence the SCDM algorithm might select columns representing better these states rather than those with projectability close to 1 at low energy, that are essential and physically relevant to include. In these cases, the corresponding band interpolation shows large oscillations and has large errors with respect to the DFT band structure in large portions of the BZ. We need therefore to choose a smaller value µ < µ fit . On the other hand, however, we note that the weight of states much above µ becomes numerically zero in Eq. (15), i.e., these states become completely unknown to the algorithm. Therefore, by choosing a too low value of µ, i.e., discarding too many relevant states, the SCDM algorithm will fail because it will have to choose J columns within a matrix of smaller rank. We need, therefore, a general and automatic recipe for choosing an appropriate, intermediate value of µ. Our choice µ = µ fit − κσ fit is guided by the consideration that states that start to have a significant component of their character outside A should be weighted in SCDM by Eq. (16) with a small weight, that is still though not exactly zero, giving the algorithm some freedom to pick up some of their character (for instance, states at energy ≥ µ fit have more than 50% of their character outside A and are weighted in SCDM with a factor ≤ 1 2 erfc(κ) (e.g., κ = 3 gives 1 2 erfc(3) ≈ 10 −5 ). In order to explain better our specific choice of κ = 3, we consider again the case of tungsten for the SCDM+MLWF case and we report in Fig. 9 the final total spread Ω (left-hand side) and the band distance η (right-hand side) as a function of a range of values of µ and σ. In particular, in the case of entangled bands, we generalise the definition of η by introducing a smearing, as we have mentioned in the previous section. More specifically, we extend the definition of the distance between DFT and Wannier-interpolated bands to: wheref and f DFT(Wan) nk (ν, τ ) is the Fermi-Dirac distribution for the state at energy ε DFT(Wan) nk , ν is a fictitious chemical potential and τ is a smearing width computed on the direct (ε DFT nk ) and Wannier-interpolated (ε Wan nk ) band structures. As in Sec. IV, we take into account the possibility that significant differences between band structures may occur only in sub-regions of the Brillouin zone or in small energy ranges, so we also compute In particular, the value of ν in f nk (ν, τ ) is set to 1 eV above the Fermi energy and the smearing width τ is 0.1 eV. In this way, only states up to slightly more than 1 eV above the Fermi level have a weight significantly different from zero when comparing band structures. In both panels of Fig. 9, we also show the line representing µ = µ fit − 3σ to discuss our choice of κ = 3, as well as the point (µ fit − 3σ fit , σ fit ) on this line. Our target is to have η as small as possible, indicating a good interpolation of the band structure. As visible in Fig. 9, and as mentioned in the previous two paragraphs, large values of µ and σ degrade significantly the quality of the band interpolation: in this case there are many states at high energy with a non-negligible weight and the QRCP, being a greedy algorithm, might select a subspace that better represents these states rather than the lowest energy states. It can also be seen that a larger µ, which results in more states with higher weight, gives the SCDM algorithm more freedom in the choice of the subspace, which in turn results in a lower total spread Ω (at the expenses of a potentially worse interpolation).
On the other hand, also moving to the region of small µ and σ is detrimental for the quality of the band interpolation (and partially also for the value of Ω). Even if the values of η in this region are not so large as in the region of large µ and σ, the quality of the interpolation is much less robust and both η and Ω depend strongly on the precise values of the two parameters. In this case, we are discarding relevant states from the initial space used for the column selection of F k Ψ k , therefore removing important information needed by the method for a good interpolation.
Our choice of κ = 3, thus, together with σ = σ fit , allows us to locate our choice of (µ, σ) in the intermediate region where η is small and both η and Ω are relatively insensitive to small variations of the two parameters. Ultimately, this specific choice for κ will be justified and validated in our high-throughput study of Sec. V B, where we show that the automated algorithm resulting from this choice is robust when tested on 200 chemically and structurally different materials, whose full list is available in Ref. [57].
We also emphasise here that the choice of µ and σ plays two different roles: the first is to give a relative weight to the states at the anchor point, namely Γ, that are used for the SCDM column selection; the second is to have a smooth dependence of the subspace as a function of k, therefore resulting in a small Ω I .  Figure 9: Bands distance η (left panels) and total position spread Ω (right panel) using SCDM+MLWF for tungsten (W), as a function of the SCDM input parameters µ and σ. The blue line represent µ = µ fit − 3σ where the red dot corresponds to the choice dictated by our protocol µ = µ fit − 3σ fit . The smearing function to compute η has smearing τ = 0.1 eV and ν is set to 1 eV above the Fermi energy.

B. High-throughput verification
In this section we present the results of the highthroughput calculations for the general case of 200 materials that have been chosen to cover a large region of structural (12 different Bravais lattices) and chemical (67 different elements) space. The free parameters in the SCDM method have been chosen by the automatic procedure outlined in the previous section. The structure of this section parallels the one for isolated bands; in particular, we make use of the bands distance η introduced in Eq. (24) to quantitatively assess the Wannier interpolation. In the case of metals, we also need to appropriately select the value of the fictitious chemical potential ν and of the smearing width τ in the distribution f nk (ν, τ ) of Eq. (25) (the final values used in this work are reported in the previous section), in order for η and η max to be reliable measures for the interpolation quality of the bands of physical interest. Indeed, the Wannier-interpolated bands are not expected to reproduce accurately the dispersion of the DFT bands at high energies; and the energy up to which the Wannier-interpolated bands may be deemed to be accurate depends mainly on the number of target WFs J which, in turn, is determined in our procedure by the number of PAOs in the pseudopotentials. In most applications, however, the high-energy bands are not of interest; therefore, ν and τ should be chosen so as to define a bands distance that only takes into account the relevant low-energy bands. For most practical applications, this means for states up to a small amount (usually a few eV) above the Fermi energy.
To verify up to which energy the interpolation is accurate (for the number of PAOs in the pseudopotentials chosen in this work, see Supplementary Table S1) we show in Fig. 10 the distribution of band distances for different values of ν = ε F + ∆, with ∆ = 1, 2, 3, 4, 5 eV, and τ fixed at 0.1 eV in order to have a smooth but sharpedged Fermi-Dirac distribution. When ν is set at 4 eV or more above the Fermi energy (∆ ≥ 4 eV, bottom panels in Fig. 10), the distribution is very broad and with a long tail. In this case states much above the Fermi energy, where the Wannier interpolation does not reproduce any more the DFT band structure, are given a non-negligible weight f nk which significantly increases the value of the band distance. The distribution becomes much more narrow and closer to η = 0 eV for ∆ ≤ 3 eV; in particular, for ν = ε F + 1.0 eV, 98% of the materials have η < 50 meV. Since for many applications having a good interpolation up to 1 eV above the Fermi energy is sufficient, in the rest of this work we choose ν = ε F + 1.0 eV (for entangled bands) as a reliable measure of the quality of the interpolation in the energy region of interest.
As in the case of isolated bands, the first step is to study the effect of the k-point grid density on the interpolation, to fix the last free parameter in the calculations. As shown in Fig. 11, a grid with density ρ k = 0.2Å −1 is typically sufficient to provide accurate interpolated band structures: in particular, 94% of the materials  Fig. 3. We therefore set ρ k to 0.2Å −1 for further analysis in this section. Fig. 12a shows the Wannier-interpolated bands (red lines) for tungsten (W), a metallic system, and Fig. 12b shows the Wannier-interpolated valence bands plus few conduction bands (in red) for the insulator C 3 Mg 2 (and these can be compared with Fig. 4b for the interpolation of the valence bands only). Unlike the case of isolated bands, for entangled bands the MLWF procedure substantially increases the localisation of the resulting Wannier functions from SCDM projections, giving for instance a ∆Ω Ω MLWF between 20 − 60% for 75% (149/200) of materials, with 30 materials showing a 60% or more increase in ∆Ω Ω MLWF , see Fig. 13. We now look at how the difference in spread due to the MLWF procedure correlates with the difference in the quality of the interpolated band structures. Although the correlation is not as strong as in the case of isolated bands, it can be seen (Fig. 14) that the trend is almost reversed: reducing the spread tends to worsen the quality of the band interpolation. In fact, the majority of systems (71%, 142/200) show a positive change in ∆η, meaning that SCDM-only provides better interpolation. The main reason behind this effect is that, in the selection of the optimal manifold S(k), the SCDM algorithm might include contributions from higher energy states. The subsequent MLWF step does not use information on the target band structure. Therefore, while mixing the states via the U matrix to minimise the spread, such spurious contributions can be distributed on the lower-energy states and, as a consequence, worsen the interpolation quality. However, we emphasise that in most cases, even when the MLWF algorithm increases the value of η, it does so only marginally: in 182 out of 200 systems (91%) the MLWF scheme either increases η by less than 5 meV or reduces it. More in detail, 163 out of these 182 materials show a variation |∆η| within only 5.0 meV, and only one system among these exhibits η MLWF > 20 meV. Moreover, for the remaining 19 (out of 182) systems the MLWF procedure improves the bands interpolation, notably yielding η MLWF < 20 meV for all of them. Finally, for the remaining 18 systems (9%), the MLWF scheme worsens the results with |∆η| > 5 meV and only in 6 cases the interpolation quality is quite poor (η MLWF > 20 meV). In all these cases, a possible reason for failure might be related to the choice of columns in the SCDM algorithm, which is performed only at Γ (see discussion in Sec. II B 2), for materials where the relative order of electronic states at Γ and at the BZ boundary is inverted. In this situation, spurious contributions might enter into the QR decomposition as discussed above.

VI. AIIDA WORKFLOWS
AiiDA 25 is a python materials' informatics platform to automate, manage and coordinate simulations and workflows, and to encourage sharing of both the resulting data and the workflow codes used to generate them. While general in its design, its plugins cover many materials science codes, including Quantum ESPRESSO 68 and Wannier90. 69 Our implementation of the SCDM method inside the open-source code Quantum ESPRESSO makes it available to any researcher. Moreover, our protocol for the choice of the SCDM parameters discussed in Sec. V A describes an effective procedure to automatically compute the Wannier functions of any material. However, the actual computation starting only from the crystal coordinates is non-trivial. The choice of numerical parameters (cutoffs, k−point grid density, convergence parameters) requires some prior knowledge and experience. Moreover, the full simulation for each material involves a complex sequence of steps, requiring a user to run over 10 different executables. Therefore, we have implemented the full procedure as AiiDA workflows, making it thus possible to repeat seamlessly the calculations for many different materials with minimal effort.
Furthermore, AiiDA keeps track of the provenance of the data generated in the simulations in a fully automated way, in the form of a directed graph (see Fig. 15 for an example of the provenance tracked for one material), where nodes can be calculations, workflows or data. This means that any researcher accessing the Ai-iDA database can inspect not only the final data, but also explore which calculation generated it, its relevant (raw and parsed) outputs and the complete set of its input parameters, and see how these input data were, in turn, obtained as output of previous calculations, traversing the graph up to the original input crystal structure.
The AiiDA workflows that we have written start by calling existing subworkflows available in the AiiDAquantumespresso 68 plug-in that, given a crystal structure, perform a variable-cell atomic relaxation to obtain the converged DFT charge density. These workflows also contain useful heuristics and recovery mechanisms to reach convergence in case of common problems (e.g., by changing the diagonalisation algorithm) as well as automatic selection of parameters, including pseudopotentials and cutoffs from the SSSP library. 58 Once the charge density is computed, the workflow first standardises the cell using the symmetry-detection library spglib 70 and the seekpath 59 library that, in addition, provide a standardised band-structure path. Then, it proceeds along two parallel branches: on one side, it computes the DFT band structure along the suggested path. In parallel, it computes the Wannier functions: if first computes wavefunctions on a full uniform grid using a non-selfconsistent Quantum ESPRESSO calculation, and then computes the PDOS, the projectabilities, and fits them to obtain the µ and σ parameters for the SCDM. Using these The SCDM+MLWF procedure provides Wannier functions that are substantially more localised with respect to SCDM-only, with a relative variation between 20 − 60% for most materials.
data, it prepares the Wannier90 input file and runs it in pre-processing mode to generate the input file needed by the code interfacing Quantum ESPRESSO with Wannier90 (pw2wannier90.x). The latter is then run to compute quantities needed by Wannier90, including the A (k) matrices obtained with the SCDM method. Finally, the workflow drives the execution of Wannier90 to compute the (maximally-localised) Wannier functions and produce the output quantities of interest (spreads, interpolated band structure on the same path of the DFT code, plots of the Wannier functions, etc.).
In an effort to improve the verification and dissemination of computational results, and in order to make the present work available to all, we are distributing all codes and workflows discussed here within a preconfigured virtual machine (VM) 57 based on the Quantum Mobile VM available on the Materials Cloud. 71 The relevant quantum codes (Quantum ESPRESSO, Wannier90) and the informatics' platform AiiDA come pre-installed and configured in the VM, ready to run through the workflows described above. A simple README file guides new users in the installation of the VM and in the execution of the workflow, to compute-with essentially no user intervention-the interpolated band structure of a material of choice.

VII. CONCLUSIONS
We have presented an approach to generate a set of maximally localised Wannier functions in an automated way that has the advantage of being simple, robust and applicable also in the more general case of so-called entangled bands. The high sensitivity of iterative minimisation algorithms to the initial conditions, which was a long-standing problem in particular for the entangledband case, is overcome by employing the selected columns of the density matrix 1,21 (SCDM) algorithm to automatically choose the initial subspace. For the Wannierisation of isolated bands, SCDM is a parameter-free method, whereas for entangled bands two real numbers µ and σ must be specified, whose appropriate choice is critical for the success of the method, in addition to the target dimensionality of the manifold to be described (i.e., the number of Wannier functions). We have proposed and validated a protocol to choose these parameters by leveraging information encoded in the projectability of the Bloch states on pseudo-atomic orbitals. We found that the SCDM method works very well for band-structure interpolations, but does not perform as well for other kind of applications where, for instance, a specific symmetry character of the WFs is desirable.  Figure 15: The provenance graph automatically generated by AiiDA when running a Wannier90 calculation using Quantum ESPRESSO as the input code for an InSe crystal. Red arrows represent caller-called relationships between a workflow and a subworkflow or a calculation; continuous lines connect calculations to their inputs and to the outputs they create, while dotted lines connect workflows to the data they return. In the top-right part of the graph, a set of workflows drive variable-cell relaxations of the initial structure via Quantum ESPRESSO; the central part contains the self-consistent, non-self-consistent and band-structure Quantum ESPRESSO calculations; in the bottom-left part are locate the calculations computing the projection of the wavefunctions on a localised atomic basis set, and in the bottom-right part of the graph we can find the Wannier90 calculation, producing a set of output nodes that includes the Wannier-interpolated band structure.
To make the method available to any researcher, we have implemented the SCDM algorithm in pw2wannier90, part of the open-source Quantum ESPRESSO distribution, and added corresponding functionality to the open-source Wannier90 code. We have also discussed how the full procedure is implemented as AiiDA 25 workflows, encoding the knowledge that is needed to perform all steps (DFT simulations, selection of the parameters, Wannierisation) into an automated software. This enables MLWFs to be obtained and used to calculate material properties by providing the crystal structure of a material as the only input. Furthermore, we are distributing publicly and freely all codes and workflows discussed in this work within a virtual machine 57 preconfigured with the open source codes AiiDA, Quantum ESPRESSO and Wannier90. This VM allows anyone to explore and reproduce straightforwardly the present results without the need to install or configure anything, and without the need of implementing again workflows and algorithms, in the true spirit of Open Science. In addition, interested researchers are not constrained to re-run the calculations performed in this work, but can perform their own simulations, either with different parameters or on new materials. To the best of our knowledge, this is the first time that such level of reproducibility is offered accompanying a scientific paper in the field of DFT simulations.
We have demonstrated the robustness of the present approach by carrying out high-throughput calculations on a dataset of 200 bulk crystalline materials, of which 81 are insulators, spanning a wide chemical and structural space. The main metric we used to assess the results is the so-called band distance, 58 quantifying the difference between the Wannier-interpolated band structures and the corresponding direct DFT band structures. In particular, we obtain excellent interpolations: for entangled bands, 97% of the materials show an average bands distance η < 20 meV and 72% show η < 5 meV. For the insulating subset, when limiting to valence bands only, 93% show η < 2 meV.
We believe that this work is a significant step forward towards completely automated high-throughput calculations of advanced materials properties exploiting Wan-nier functions.

Appendix A: SCDM implementation in Quantum ESPRESSO
To implement the SCDM method one needs the wavefunctions from the ab initio code represented on a real space grid, see Eq. (14). Since these are not directly accessible to Wannier90, we decided to implement the method in one of the open-source DFT-to-Wannier90 interface packages available. In particular, we have chosen the pw2wannier90 FORTRAN code, distributed with the open-source Quantum ESPRESSO suite. 23 Our SCDM implementation is available since the v6.3 release of Quantum ESPRESSO. It includes the extension of the method to k−points and to entangled bands, and it is parallelised using the Message Passing Interface (MPI).
To compute the A (k) mn projection matrices using SCDM, the auto projections keyword must be set to .true. in the Wannier90 input file. In addition, the following keywords should be defined in the pw2wannier90.f90 input file: scdm proj, scdm entanglement, scdm mu and scdm sigma. In particular, scdm proj is a boolean flag to enable the SCDM method. scdm entanglement is a string defining the functional form of the f (ε) function in Eq. (15). In the cases described in this paper, the value is either isolated (isolated bands) or erfc (entangled bands with the f (ε) of Eq. (16)). An additional choice we implemented is gaussian, see Ref. [1] for its functional form. Finally, scdm mu and scdm sigma (not needed if scdm entanglement is isolated) define, respectively, the values of µ and σ in Eq. (16) (in eV).
In pw2wannier90.x, the QRCP factorisation of the Ψ † k=Γ matrix is obtained through the LAPACK routine ZGEQP3. Presently the factorisation is performed on a single MPI process (since ZGEQP3 is not available in the parallel ScaLAPACK routines) and the resulting permutation matrix Π is broadcast to all processes. After the computation of the non-orthogonal SCDM functions, a Löwdin orthogonalisation is performed. This step is not needed when providing the A (k) mn matrices to Wannier90, since the same orthogonalisation is performed by the code before the start of the minimisation. However, having the orthogonalisation step also in the pw2wannier90.x interface allows users to directly employ the SCDM functions without further processing, if needed.
As a final note, we emphasise that when ultrasoft pseudopotentials are employed, the ψ nk (r) wavefunctions satisfy a generalised orthogonality condition with a non-trivial metricŜ being a function of the core augmentation charges. 54 In this case the u nk stored by Quantum ESPRESSO are not orthonormal, resulting in Ψ being non unitary. However, in practice this usually has only a marginal effect on the results. Indeed, as we have shown, the algorithm manages to find good Wannier functions also when employing ultrasoft pseudopotentials and therefore no adaptation has been applied for the ultrasoft case.
where the additional columnsQ of Q are chosen to complete the columns of Q to an orthonormal basis of R n G (always possible) and R extends R with (n G − J) additional rows of zeros.
We want now to prove that P Π = Q R is a QRCP decomposition of P . Indeed, by multiplying by blocks the two matrices Q and R , we get Q R = Q R + Q0 = Q R = P Π by virtue of Eq. (C3). Moreover, Q is a unitary matrix by construction, and R is clearly an upper-triangular matrix since R is according to point 3 of Appendix B, and the diagonal elements are still sorted in decreasing magnitude order since the additional elements are all zero. Therefore, we have shown that the same permutation matrix Π obtained by applying the QRCP to Ψ † is a valid QRCP permutation matrix also for P .
A different, equivalent approach to show the same result is to observe that the (complex) scalar product v 1 · v 2 ≡ (v * 1 ) T v 2 between columns of P is the same as the scalar product of the columns of Ψ † . Indeed, we first note that, as it can be easily proven from its explicit expression Eq. (12), P is a projector and it holds P 2 = P and P † = P . Therefore, we have P † P = P . But the elements of P † P are nothing else than the scalar products of the columns of P , and therefore P ij = p i · p j , with p i indicating the i−th column of P . At the same time, from the definition of P = ΨΨ † = (Ψ † ) † Ψ † we immediately notice that the elements of P are also the scalar products of the columns of Ψ † , i.e. the complex conjugate ψ * i of the wavefunctions of the system, proving our statement that P ij = p i · p j = ( ψ i |ψ j ) * . (C6) Appendix D: Geometrical interpretation of the column selection in the QRCP algorithm QRCP is a greedy algorithm, where the Π matrix is constructed by picking the columns one by one to obtain the condition |R 11 | ≥ |R 22 | ≥ . . . ≥ |R JJ |. In the case of the QRCP decomposition of a P matrix, the first column to be picked (p 1 ) i ≡ (P Π) i1 is chosen as the one with largest norm. This can be easily proven by noting that because of the triangular form of R. Moreover, since the columns of Q have norm one, p 1 = |R 11 |, that by construction (see point 3 of Appendix B) is the largest possible.
More generally, the j−th column p j is chosen to maximise the norm of the component p ⊥ j orthogonal to the subspace S j−1 spanned by the previous (j − 1) columns. To prove this, let us first write p j = p j + p ⊥ j , where p j is the projection of p j within S j−1 . We first note (again due to the triangular form of R) that in general the first j columns of Q also span the space S j and, moreover, they are a orthonormal basis set for S j since the Q columns are orthonormal. Furthermore p j is, by definition, in the S j subspace. Therefore, we can write the j−th column in this basis set of S j as (p j ) i = (P Π) ij = (QR) ij = j m=1 Q im R mj , and, thanks to the orthonormality of the {q m } basis (q i being the i−th column of Q), we have or equivalently in vector form p ⊥ j = q j R jj . Therefore, the norm of this orthogonal component is simply p ⊥ j = |R jj | that, again, is chosen by the algorithm to have maximal value (in order to have decreasing diagonal elements of R), therefore proving our intuitive explanation of the QRCP column selection.
To give a more physical interpretation of the column selection in terms of the charge density or wavefunctions, we observe that from Eq. (C6) we know that the square modulus of the i−th column of P is p i 2 = P ii , and the diagonal element of P is simply ρ(r i ), i.e. the charge density at the discretised grid point r i . Therefore, the algorithm will choose the first column p Π(1) as the one corresponding to the point in space r i with maximal charge density (i.e., the projection of a delta-like function centred on r i ).
1. Define an initial matrix P (0) = P being the density matrix of the system.
2. At every step k ≥ 1, choose Π(k) as the index of the column where the matrix P (k−1) has maximum diagonal element. Also, we define the k−th Cholesky vector c k as the Π(k)−th column of P (k−1) , rescaled by the inverse square root of the corresponding diagonal element: 3. Define the matrix P (k) for the next iteration as follows: (where c k · c † k indicates a matrix product).
4. Iterate the previous two points until the needed number of selected columns is obtained.
We can now show that this approach is equivalent to the selection of columns of the QRCP algorithm. In particular, in the first step, the Cholesky approach selects the column corresponding to the largest diagonal element of P , which is exactly the same choice as the QRCP algorithm, as discussed in Appendix D.
At the second step (k = 1), substituting Eq. (E1) in Eq. (E2) and using P (0) = P , we have In particular, we can notice now that the diagonal elements P (1) jj of P (1) can be written as (where we have used P † = P ) and therefore the choice of j = Π(2) based on the largest diagonal element of P (1) , as prescribed by the Cholesky algorithm, is equivalent to the QRCP choice maximising Eq. (D2). Finally, we note that the Π(1)−th column of P (1) is composed only by zeros (and analogously for the Π(1)−th row since the P (i) matrices are Hermitian), since P (1) iΠ(1) = P iΠ(1) − P iΠ(1) P * Π(1)Π(1) P Π(1)Π (1) = P iΠ(1) − P iΠ(1) = 0 (E5) (where we have used the fact that the diagonal elements of P are real). This fact, beside allowing to give numerical stability to the Cholesky-orbital algorithm by forcing these elements to be numerically zero, allows us to "remove" the zero row and column from P (1) and repeat the reasoning by induction for all subsequent Cholesky vectors, working with smaller and smaller matrices. Equivalently, one could understand more intuitively the result by noting that the Cholesky vectors of Eq. (E1) are normalised to 1 because of Eq. (C6). Therefore, Eq. (E2) constructs a new projection operator P (k) projecting on the subspace of the span of P (k−1) that is also orthogonal to c k , and then the Cholesky algorithm selects the largest vector in this subspace, that is exactly what the QRCP algorithm also does, as discussed in Appendix D. . The DFT band-structure is also shown for reference (dotted black). It is worth clarifying that regardless of the initial projections, after a full minimisation-solid red line in all three panels-the Wannier interpolation is extremely good, particularly for the valence manifold. In fact, the three methods give almost indistinguishable results. Table S1: List of number of valence electrons included in the pseudopotential (Z val ) and atomic pseudoorbitals included in the pseudopotential file (1s refers to an atomic s pseudo-orbital without radial nodes, 2p to an atomic p pseudo-orbital with one radial node, . .