Closest Wannier functions to a given set of localized orbitals
The original has been published
in Phys. Rev. B 110, 125115 (2024).
Abstract
A noniterative method is presented to calculate the closest Wannier functions (CWFs) to a given set of localized guiding functions, such as atomic orbitals, hybrid atomic orbitals, and molecular orbitals, based on minimization of a distance measure function. It is shown that the minimization is directly achieved by a polar decomposition of a projection matrix via singular value decomposition, making iterative calculations and complications arising from the choice of the gauge irrelevant. The disentanglement of bands is inherently addressed by introducing a smoothly varying window function and a greater number of Bloch functions, even for isolated bands. In addition to atomic and hybrid atomic orbitals, we introduce embedded molecular orbitals in molecules and bulks as the guiding functions, and demonstrate that the Wannier interpolated bands accurately reproduce the targeted conventional bands of a wide variety of systems including Si, Cu, the TTFTCNQ molecular crystal, and a topological insulator of Bi${}_{2}$Se${}_{3}$. We further show the usefulness of the proposed method for calculating effective atomic charges. These numerical results not only establish our proposed method as an efficient alternative for calculating WFs, but also suggest that the concept of CWF can serve as a foundation for developing novel methods to analyze electronic structures and calculate physical properties.
I INTRODUCTION
Wannier functions (WFs) [1, 2, 3] play a central role in analyzing electronic structures of real materials and advancing electronic structure methods along with density functional theory (DFT) [4] and the other electronic structure theories [5, 6]. A widely adopted method for calculating WFs involves maximizing their localization, which can be reformulated as minimizing the spread function characterizing the variance of the WFs in real space [7, 8]. The concept of maximal localization leads to an elegant formulation, resulting in maximally localized Wannier functions (MLWFs), and the following compact representation of the Hamiltonian allows us to analyze the electronic structures of real materials using a localized orbital picture. However, it is also recognized that the minimization of the spread function often encounters local minima, particularly in largescale systems with complicated electronic structures [3]. To overcome this challenge, methods aimed at automated highthroughput Wannierisation have been developed [9, 10, 11]. In addition, other methods for generating (nonorthogonal) WFs have been proposed based on projection methods [12, 13, 14], which produce compact atomic like orbitals. The orthogonalization of such nonorthogonal orbitals obtained by the projection can be achieved using the Löwdin orthogonalization procedure [15, 16]. These orthogonalized orbitals are utilized as an initial guess of the WFs in the minimization of the spread function to obtain MLWFs [7, 17]. Apart from the role of the initial guess, it is worth noting that the Löwdin orthogonalized orbitals possess an important variational property, which in particular minimizes the sum of the squared distances between the orthogonal and original nonorthogonal orbitals in the Hilbert space [18]. As long as the original nonorthogonal orbitals are localized, the localization of the Löwdin orthogonalized orbitals in real space is guaranteed in a sense of minimization of the distance. In quantum chemistry, this property is leveraged to develop methods for generating localized orthogonal orbitals. These methods, which include the natural atomic and bond orbital method [19] and the intrinsic atomic orbital method (IAO) [20], involve obtaining orthogonalized orbitals through Löwdin type orthogonalization, giving heavy weight to occupied states via the density matrix. On the other hand, in solid state physics, the variational property inherent in the Löwdin orthogonalization has not yet been fully explored. This is because the Löwdin orthogonalization requires the calculation of ${S}^{1/2}$ for the overlap matrix $S$, which causes numerical instability for nearly illconditioned matrices. Owing to this difficulty, advancements in methodologies along these lines appear to be slow. It should be mentioned that avoiding the calculation of ${S}^{1/2}$ has already been proposed based on singular value decomposition (SVD) to calculate an initial guess of WFs for the maximal localization [8]. However, the relationship between the WFs obtained through this procedure and the variational principle has not been sufficiently discussed so far.
Our study aims to develop an efficient and robust method for generating the closest Wannier functions (CWFs) to a given set of localized guiding functions, based on a variational principle and polar decomposition of the projection matrix. We will demonstrate that exploiting the variational property leads to a versatile method that eliminates the need for iterative optimization and avoids complications related to the gauge choice. The projection method based on the polar decomposition has already been employed as a workaround [21, 22] with successful applications [23, 24]. However, to the best of our knowledge, no study has thoroughly discussed the mathematical and physical significance of this method based on the closest principle. Our study establishes the theoretical foundation by formulating the closest principle as a variational problem, and clarifies the relationship between the Löwdin orthogonalization and the polar decomposition. We will also introduce a smoothly varying window function that enables us to flexibly control targeted states to be Wannierized in a similar way to that proposed in Ref. [10]. The present study offers a more rigorous justification for these practices. With the theoretical advancement, the concept of CWF can serve as a foundation for developing novel methods to analyze electronic structures and calculate physical properties as exemplified in the analysis of effective atomic charges in Sec. IV.
The structure of this paper is as follows: In Sec. II, we present the theory for calculating CWFs. Sec. III provides a detailed discussion on the implementation of the method. A series of benchmark calculations are presented in Sec. IV. Finally, in Sec. V, we summarize the theory of CWFs and suggest the potential role of CWFs as a foundation for developing novel methods.
II THEORY
Let us start by introducing Bloch functions $\{\varphi \}$, which can be obtained by solving the KohnSham (KS) equation [25] within the density functional theory (DFT) [4], normalized as
$\u27e8{\varphi}_{{\mathbf{k}}_{1}\mu}{\varphi}_{{\mathbf{k}}_{2}\nu}\u27e9={N}_{\mathrm{BC}}{\delta}_{{\mathbf{k}}_{1}{\mathbf{k}}_{2}}{\delta}_{\mu \nu},$  (1) 
where $\mathbf{k}$ and $\mu $ are the kvector and band index, respectively, and ${N}_{\mathrm{BC}}$ is the number of primitive cells in the Bornvon Karman (BvK) boundary condition. We consider the projection of a localized guiding function $Q$ onto the Bloch functions $\{\varphi \}$ weighted with a window function $w(\epsilon )$ as
${L}_{\mathbf{R}p}\u27e9={\displaystyle \frac{1}{{N}_{\mathrm{BC}}}}{\displaystyle \sum _{\mathbf{k}\mu}}{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \mathbf{R}}{\varphi}_{\mathbf{k}\mu}\u27e9{a}_{\mathbf{k}\mu ,p}$  (2) 
with
${a}_{\mathbf{k}\mu ,p}=w({\epsilon}_{\mathbf{k}\mu})\u27e8{\varphi}_{\mathbf{k}\mu}{Q}_{\mathrm{\U0001d7ce}p}\u27e9,$  (3) 
where $\mathbf{R}$ and $p$ are the translational lattice vector and index of the localized orbital, respectively, ${\epsilon}_{\mathbf{k}\mu}$ is the eigenvalue of the KS equation, and ${Q}_{\mathrm{\U0001d7ce}p}\equiv {Q}_{\mathbf{R}p}$ $(\mathbf{R}=\mathrm{\U0001d7ce})$. In the summation of Eq. (2), the numbers of kpoints and the Bloch functions to be included per kpoint are ${N}_{\mathrm{BC}}$ and ${N}_{\mathrm{band}}$, respectively. Throughout this paper, we do not consider the spin dependency on the formulation for the sake of simplicity, but the generalization is straightforward. The localized functions $\{{Q}_{\mathbf{R}p}\}$, such as atomic orbitals, hybrid orbitals, or molecular orbitals (MOs), are assumed to localize in the unit cell specified by $\mathbf{R}$. When the window function $w(\epsilon )$ is taken to be unity, Eq. (2) is nothing but an expansion of ${Q}_{\mathbf{R}p}$ with Bloch functions from the identity relation $\frac{1}{{N}_{\mathrm{BC}}}{\sum}_{\mathbf{k}\mu}{\varphi}_{\mathbf{k}\mu}\u27e9\u27e8{\varphi}_{\mathbf{k}\mu}=\widehat{I}$. To introduce an expansion of ${Q}_{\mathbf{R}p}$ by the Bloch functions $\{{\varphi}_{\mathbf{k}\mu}\}$ in a subspace, we use the following window function $w(\epsilon )$:
$w(\epsilon )={\displaystyle \frac{1\mathrm{exp}({x}_{0}+{x}_{1})}{(1+\mathrm{exp}({x}_{0}))(1+\mathrm{exp}({x}_{1}))}}+\delta $  (4) 
with the definition of ${x}_{0}\equiv \beta ({\epsilon}_{0}\epsilon )$ and ${x}_{1}\equiv \beta (\epsilon {\epsilon}_{1})$, where $\beta =\frac{1}{{k}_{\mathrm{B}}T}$ with the Boltzmann constant of ${k}_{\mathrm{B}}$ and a temperature of $T$. The window function of Eq. (4) is obtained by subtracting 1 from the sum of the two functions of $1/(1+\mathrm{exp}(x))$. It is also possible to introduce two temperatures ${T}_{0}$ and ${T}_{1}$. The usefulness of the two temperature scheme will be demonstrated for the case of GaAs in Sec. IV. The last term $\delta $ in Eq. (4) is a small constant, e.g., ${10}^{12}$, which is introduced to avoid the illconditioning of the matrix consisting of ${a}_{\mathbf{k}\mu ,p}$ by Eq. (3) as will be discussed later. As shown in Fig. 1, the parameters of ${\epsilon}_{0}$ and ${\epsilon}_{1}$ ($$) determine the energy range where the Bloch states are included in the expansion of Eq. (2) with a large weight, and the temperature $T$ provides the degree of smearing around ${\epsilon}_{0}$ and ${\epsilon}_{1}$ in $w(\epsilon )$. To control the degree of smearing, we introduce the temperature of $T$ so that one can intuitively understand the degree of smearing. It should be noted that the choice of the window function by Eq. (4) is not unique, and the other choices should give almost equivalent results as long as they are chosen properly. One may also consider that the window function has no special advantages over the outer and inner window scheme [8] and merely functions in a similar manner. However, the window function defined by Eq. (4) exhibits superior performance in the disentanglement of bands and calculation of effective charges. The advantages will be clearly demonstrated in Secs. IV A and IV B, respectively.
The function ${L}_{\mathbf{R}p}$ defined by Eq. (2) must be similar to the localized function ${Q}_{\mathbf{R}p}$, but they are generally nonorthogonal to each other. We now consider generating a set of closest Wannier functions (CWFs) $\{{W}_{\mathbf{R}p}\}$ to a set of localized functions $\{{L}_{\mathbf{R}p}\}$ in the sense that the sum of the squared distance between $L$ and $W$ in the Hilbert space is minimized. As well as Eq. (2), noting that the CWFs can be expressed [3] as
${W}_{\mathbf{R}p}\u27e9={\displaystyle \frac{1}{{N}_{\mathrm{BC}}}}{\displaystyle \sum _{\mathbf{k}\mu}}{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \mathbf{R}}{\varphi}_{\mathbf{k}\mu}\u27e9{b}_{\mathbf{k}\mu ,p},$  (5) 
and defining the residual function $R$:
${R}_{\mathbf{R}p}\u27e9={L}_{\mathbf{R}p}\u27e9{W}_{\mathbf{R}p}\u27e9,$  (6) 
the distance measure (DM) function $F$ we minimize is given by
$F[B]$  $=$  $\sum _{p}}\u27e8{R}_{\mathrm{\U0001d7ce}p}{R}_{\mathrm{\U0001d7ce}p}\u27e9,$  (7)  
$=$  $\frac{1}{{N}_{\mathrm{BC}}}}{\displaystyle \sum _{\mathbf{k}}}X[B,\mathbf{k}]$ 
with
$X[B,\mathbf{k}]=\mathrm{tr}\left[\left({A}^{\u2020}(\mathbf{k}){B}^{\u2020}(\mathbf{k})\right)\left(A(\mathbf{k})B(\mathbf{k})\right)\right],$  (8) 
where the elements of the matrices $A(\mathbf{k})$ and $B(\mathbf{k})$ are given by ${a}_{\mathbf{k}\mu ,p}$ and ${b}_{\mathbf{k}\mu ,p}$, respectively. The second line of Eq. (7) is obtained by considering the orthonormality of Eq. (1), and $X[B,\mathbf{k}]$ of Eq. (8) is regarded as the squared Frobenius norm of $\left(A(\mathbf{k})B(\mathbf{k})\right)$ [26]. The same Bloch functions as in Eq. (2) are included in the summation of Eq. (5). The number of CWFs to be generated is ${N}_{\mathrm{CWF}}$ per unit cell, being equivalent to the number of the guiding functions per cell, and ${N}_{\mathrm{CWF}}$ should be smaller than or equal to ${N}_{\mathrm{band}}$ to guarantee the linear independency of the subspace spanned by the CWFs, resulting in size of ${N}_{\mathrm{band}}\times {N}_{\mathrm{CWF}}$ for the matrices $A(\mathbf{k})$ and $B(\mathbf{k})$. Assuming ${B}^{\u2020}(\mathbf{k})B(\mathbf{k})=I$, where $I$ is of size ${N}_{\mathrm{CWF}}\times {N}_{\mathrm{CWF}}$, the CWFs are ensured to form a set of orthonormal functions. Under this constraint, we consider the optimization of the DM function $F$. Furthermore, the matrix $B(\mathbf{k})$ is regarded as a part of a unitary matrix of size ${N}_{\mathrm{band}}\times {N}_{\mathrm{band}}$, and will be referred to as a partial unitary matrix in subsequent discussions.
The minimum of the DM function of Eq. (7) is obtained by choosing $B(\mathbf{k})$ as $U(\mathbf{k})$ which is calculated by the polar decomposition of $A(\mathbf{k})=U(\mathbf{k})P(\mathbf{k})$, where $U(\mathbf{k})$ and $P(\mathbf{k})$ are unitary and hermitian, respectively [27]. Let us prove the statement below. The polar decomposition of $A(\mathbf{k})$ is obtained via SVD of $A(\mathbf{k})$ as
$A=W\mathrm{\Sigma}{V}^{\u2020}=W{V}^{\u2020}V\mathrm{\Sigma}{V}^{\u2020}=UP,$  (9) 
where we dropped the dependency on $\mathbf{k}$ for simplicity of the notation, and hereafter we will denote the dependency if necessary. $W$ and $V$ are the left and right singular matrices in size of ${N}_{\mathrm{band}}\times {N}_{\mathrm{CWF}}$ and ${N}_{\mathrm{CWF}}\times {N}_{\mathrm{CWF}}$, respectively, and $\mathrm{\Sigma}$ is the singular value diagonal matrix in size of ${N}_{\mathrm{CWF}}\times {N}_{\mathrm{CWF}}$. Note that $U\equiv W{V}^{\u2020}$ and $P\equiv V\mathrm{\Sigma}{V}^{\u2020}$. It is worth mentioning that $W$ and $U$ are partial unitary matrices, and hold ${W}^{\u2020}W=I$ ($W{W}^{\u2020}\ne I$) and ${U}^{\u2020}U=I$ ($U{U}^{\u2020}\ne I$), and that $V$ is a full unitary matrix holding ${V}^{\u2020}V=I$ ($V{V}^{\u2020}=I$). We evaluate $X$ of Eq. (8) for each $\mathbf{k}$ with both the matrix $U$ obtained by the polar decomposition and an arbitrary partial unitary matrix $B$, and calculate the difference as
$X[U]X[B]$  $=$  $2\mathrm{t}\mathrm{r}\left[{\displaystyle \frac{1}{2}}\left({A}^{\u2020}B+{B}^{\u2020}A\right)P\right],$  (10)  
$=$  $2\mathrm{t}\mathrm{r}\left[\mathrm{\Sigma}D\mathrm{\Sigma}\right],$  
$=$  $2{\displaystyle \sum _{n}}{\sigma}_{n}\left({d}_{nn}1\right)$ 
with
$D={\displaystyle \frac{1}{2}}\left({V}^{\u2020}{U}^{\u2020}BV+{V}^{\u2020}{B}^{\u2020}UV\right).$  (11) 
The second line of Eq. (10) is derived by using the polar decomposition of $A$ and $P=V\mathrm{\Sigma}{V}^{\u2020}$, and ${\sigma}_{n}$ and ${d}_{nn}$ in the third line are diagonal elements of $\mathrm{\Sigma}$ and $D$, respectively. The upper bound of the diagonal elements of the hermitian matrix $D$ is found to be unity, since the matrix $D$ consists of the sum of the product of partial unitary matrices of ${V}^{\u2020}{U}^{\u2020}$ (${V}^{\u2020}{B}^{\u2020}$) and $BV$ ($UV$). Another proof for the upper bound is also given based on the CauchySchwarz inequality in the appendix. By considering the upper bound of the diagonal elements of $D$ and $0\le {\sigma}_{n}$, the third line of Eq. (10) leads to the following inequality:
$X[U]\le X[B].$  (12) 
The equality of Eq. (12) holds if $B=U$. If some of ${\sigma}_{n}$ is zero, the corresponding ${d}_{nn}$ in Eq. (10) can be arbitrarily chosen. So, we see from Eq. (11) that $U$ is not uniquely determined. If all the singular values of ${\sigma}_{n}$ are positive, all the ${d}_{nn}$s should be unity when $X[U]=X[B]$. The uniqueness of the solution can be confirmed as follows: Since the matrix $D$ consists of the products of the partial unitary matrices as discussed above, the case that all the diagonal elements ${d}_{nn}$ are unity happens when
${U}^{\u2020}B+{B}^{\u2020}U=2I,$  (13) 
which is derived by equating $D=I$ and multiplying $V$ and ${V}^{\u2020}$ from the left and right of the equation, respectively. By noting again that $B$ and $U$ in Eq. (13) are the partial unitary matrices, it is found that Eq. (13) holds if and only if $B=U$. Thus, we have proven the statement that the minimum of the DM function of Eq. (7) is uniquely determined when $B=U$ as long as all the singular values of ${\sigma}_{n}$ are positive, since $F[B]$ of Eq. (7) consists of the sum of $X[B,\mathbf{k}]$ over k. The uniqueness of $U$ itself is related to that of the SVD for the matrix $A$ of Eq. (3). If $A$ has ${N}_{\mathrm{CWF}}$ positive singular values which are nondegenerate, the SVD is uniquely determined except for the nonunique phases in $W$ and $V$. Even if there are degenerate singular values, the matrix $U$ is invariant, since the freedom of the unitary transformation $K$ for the degenerate left and right singular vectors is canceled out as $U=WK{K}^{\u2020}{V}^{\u2020}=W{V}^{\u2020}$ with $K{K}^{\u2020}=I$. Further, noting that the nonunique phases in $W$ and $V$ are canceled out when the matrix $U$ is computed as $W{V}^{\u2020}$, we conclude that $U$ is uniquely determined if $A$ has the ${N}_{\mathrm{CWF}}$ positive singular values. The violation from the positive definiteness of the singular values can be avoided by the small constant of $\delta $ in Eq. (4). On the other hand, if $\delta $ in Eq. (4) is set to be 0 and a narrow window, including fewer than ${N}_{\mathrm{CWF}}$ eigenstates for a given $\mathbf{k}$, is used with a small ${k}_{\mathrm{B}}T$, the matrix $A$ has ${N}_{\mathrm{PSV}}$ positive singular values, where $$. In this case, a subspace of the $\left({N}_{\mathrm{CWF}}{N}_{\mathrm{PSV}}\right)$ dimension is arbitrary chosen to form $W$ and $V$ of the dimension of ${N}_{\mathrm{CWF}}$ in addition to the subspace of the ${N}_{\mathrm{PSV}}$ dimension, which breaks the symmetry of CWFs. The symmetry breaking is due to the nonuniqueness of the SVD when some of the singular values are zero. The situation can be avoided by introducing a small constant $\delta $ in the window function of Eq. (4). The simple treatment restores the symmetry of CWFs, allowing us to calculate symmetry preserving CWFs for a wide variety of choices in the energy range of the window function. It is also noted that the nonunique phase of the Bloch functions is canceled out via the polar decomposition of $A$ and Eq. (5). We have the arbitrariness of phase both in calculating the Bloch functions and singular vectors $W$ and $V$, and unitary freedom $K$ for $W$ and $V$ when singular values are degenerate. However they are all canceled out as discussed above. Thus, we see that the proposed method is free from complications arising from the choice of gauge.
Using Eqs. (8) and (9), the minimum of the DM function $F$ at $B=U$ is given by
$F[U]={\displaystyle \frac{1}{{N}_{\mathrm{BC}}}}{\displaystyle \sum _{\mathbf{k}p}}{\left({\sigma}_{\mathbf{k}p}1\right)}^{2},$  (14) 
where ${\sigma}_{\mathbf{k}p}$ is a singular value of the matrix $A(\mathbf{k})$. From Eq. (14), we find that the mean squared distance between $L$ and $W$ is related to the deviation of $\sigma $ from unity.
The proposed method based on the polar decomposition of $A$ is closely related to the Löwdin orthogonalization [15] as shown below. The Fourier transform of the overlap integrals for $\{L\}$ is given by
$\sum _{\mathbf{R}}}{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \mathbf{R}}\u27e8{L}_{\mathrm{\U0001d7ce}p}{L}_{\mathbf{R}q}\u27e9={\displaystyle \sum _{\mu}}{a}_{\mathbf{k}\mu ,p}^{*}{a}_{\mathbf{k}\mu ,q},$  (15) 
which is obtained by noting that $\frac{1}{{N}_{\mathrm{BC}}}{\sum}_{\mathbf{R}}{\mathrm{e}}^{\mathrm{i}(\mathbf{k}{\mathbf{k}}^{\prime})\cdot \mathbf{R}}={\delta}_{{\mathrm{\mathbf{k}\mathbf{k}}}^{\prime}}$. By writing Eq. (15) as $S(\mathbf{k})={A}^{\u2020}(\mathbf{k})A(\mathbf{k})$ in a matrix form, and using Eq. (9), we have the following relation:
$S(\mathbf{k})={A}^{\u2020}(\mathbf{k})A(\mathbf{k})={V}^{\u2020}(\mathbf{k}){\mathrm{\Sigma}}^{2}(\mathbf{k})V(\mathbf{k}).$  (16) 
Comparing Eq. (16) with $P(\mathbf{k})=V(\mathbf{k})\mathrm{\Sigma}(\mathbf{k}){V}^{\u2020}(\mathbf{k})$, one obtains a relation $P(\mathbf{k})={S}^{1/2}(\mathbf{k})$. If $P(\mathbf{k})$ is invertible, the matrix $U(\mathbf{k})$ is given with Eq. (9) by
$U(\mathbf{k})=A(\mathbf{k}){S}^{1/2}(\mathbf{k}).$  (17) 
The matrix $U$ obtained by Eq. (17) is exactly equivalent to that by the Löwdin orthogonalization [15], and the closest property of the Löwdin orthogonalized orbitals to a given set of orbitals was proven in Ref. [18]. On the other hand, the proposed method does not require the calculation of ${S}^{1/2}$, and can be applied in a numerical stable manner even to the case that the matrix $A$ is nearly illconditioned. The definition of $A$ by Eq. (3) with the window function allows us to calculate CWFs in various respects, e.g., heavily weighting to occupied states and disentangling of localized states embedded in the wider bands as demonstrated in Sec. IV. In this sense, the method we propose can be regarded as a generalization of the Löwdin orthogonalization. We also notice that the projection method in Ref. [8] is equivalent to the proposed method in the sense that $U$ in Eq. (17) is calculated via the polar decomposition, and that the WFs obtained through the procedure [8] can be regarded as CWFs under the constraints of how the matrix $A$ is constructed.
The disentanglement of bands can be easily performed by properly selecting the window function of Eq. (3). Since the ${N}_{\mathrm{band}}\phantom{\rule{veryverythickmathspace}{0ex}}(\ge {N}_{\mathrm{CWF}})$ Bloch functions per kpoint are included in Eqs. (2) and (5), it should be emphasized that the proposed method always disentangles ${N}_{\mathrm{band}}$ states to generate ${N}_{\mathrm{CWF}}$ Wannier functions. A couple of examples will be shown in Sec. IV, including valence and lowlying conduction bands for the diamond Si and narrow $3d$bands embedded in the wider $4s$band for the face centered cubic (FCC) Cu.
Once the $k$dependent $U(\mathbf{k})$ are obtained by the polar decomposition, tightbinding (TB) parameters are calculated as expectation values of the KS Hamiltonian ${\widehat{H}}_{\mathrm{KS}}$ by summing contributions over $\mathbf{k}$ and $\mu $ as
${t}_{\mathrm{\U0001d7ce}p,\mathbf{R}q}$  $=$  $\u27e8{W}_{\mathrm{\U0001d7ce}p}{\widehat{H}}_{\mathrm{KS}}{W}_{\mathbf{R}q}\u27e9,$  (18)  
$=$  $\frac{1}{{N}_{\mathrm{BC}}}}{\displaystyle \sum _{\mathbf{k}\mu}}{\epsilon}_{\mathbf{k}\mu}{u}_{\mathbf{k}\mu ,p}^{*}{u}_{\mathbf{k}\mu ,q}{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \mathbf{R}},$ 
where ${u}_{\mathbf{k}\mu ,q}$ is the matrix element of $U(\mathbf{k})$. Like with MLWFs, the TB parameters can be used for the Wannier interpolation in calculations of band structures and physical quantities.
The computational procedure to generate the CWFs is summarized as follows:

1.
Determining ${\epsilon}_{0}$, ${\epsilon}_{1}$, and ${k}_{\mathrm{B}}T$ in Eq. (4) by checking the band dispersion of a system of interest. The parameters should be chosen properly so that the targeted eigenstates can be included.

2.
Choosing a set of localized orbitals {$Q$}. Atomic orbitals, hybrid orbitals, and MOs might be possible choices. Depending on the symmetry of the targeted eigenstates, a proper set of localized orbitals needs to be chosen, e.g., $d$orbitals need to be employed in the generation of CWFs for the $d$bands in FCC Cu as shown later on.

3.
Calculation of $A(\mathbf{k})$ by Eq. (3). The calculation of the overlap $\u27e8{\varphi}_{\mathbf{k}\mu}{Q}_{\mathrm{\U0001d7ce}p}\u27e9$ depends on the implementation of the KS method.

4.
Performing the SVD of $A(\mathbf{k})$. A proper numerical library might be used.

5.
Calculation of $U(\mathbf{k})$. Using the left and right singular matrices, $W(\mathbf{k})$ and $V(\mathbf{k})$, of $A(\mathbf{k})$, $U(\mathbf{k})$ is calculated as $W(\mathbf{k}){V}^{\u2020}(\mathbf{k})$.
 6.
It should be stressed that the minimization of $F$ can be efficiently performed in a noniterative manner by the steps 1 to 6 using the polar decomposition via the SVD. The computational cost of each step is estimated to be $O({N}_{\mathrm{BC}}{N}_{\mathrm{band}}{N}_{\mathrm{CWF}}{N}_{\mathrm{basis}})$, $O({N}_{\mathrm{BC}}{N}_{\mathrm{band}}^{3})$, $O({N}_{\mathrm{BC}}{N}_{\mathrm{band}}{N}_{\mathrm{CWF}}^{2})$, or $O({N}_{\mathrm{BC}}{N}_{\mathrm{band}}{N}_{\mathrm{CWF}}{N}_{\mathrm{basis}})$ for the step 3, 4, 5, or 6, respectively, where ${N}_{\mathrm{basis}}$ is the number of basis functions to expand the KS orbitals. If the localized basis functions are used, the computational cost of the step 3 becomes $O({N}_{\mathrm{BC}}{N}_{\mathrm{band}}{N}_{\mathrm{CWF}})$.
For isolated systems, the Brillouin zone sampling is limited to only the $\mathrm{\Gamma}$ point. No modification of the proposed method is needed to calculate the CWFs.
The theoretical framework we have discussed so far is general for any implementation of the DFTKS method. However, the choice of the the localized functions $\{Q\}$ in Eq. (3) may depend on the implementation. So, we will discuss how the localized functions $\{Q\}$ can be properly generated in Sec. III.
III IMPLEMENTATION
III.1 General
We have implemented the proposed method into the OpenMX DFT software package [28, 29, 30] which is based on normconserving pseudopotentials (PPs) [31, 32] and optimized pseudoatomic orbitals (PAOs) [33, 34] as basis set. All the benchmark calculations were performed with a computational condition of a production level. The basis functions used are listed in Table 1. In the abbreviation of basis functions such as H7.0s2p2d1, H stands for the atomic symbol, 7.0 the cutoff radius (Bohr) in the generation by the confinement scheme [33, 34], and s2p2d1 means the employment of two, two, and one optimized radial functions for the $s$, $p$, and $d$orbitals, respectively. The radial functions were optimized by a variational optimization method [33]. These basis functions we used can be regarded as at least double zeta plus double polarization basis sets if we follow the terminology of Gaussian or Slatertype basis functions. Valence states included in the PPs are listed in Table 1. All the PPs and PAOs we used in the study were taken from the database (2019) in the OpenMX website [28], which were benchmarked by the delta gauge method [35]. Real space grid techniques were used for the numerical integrations and the solution of the Poisson equation using fast Fourier transform (FFT) with the energy cutoff of 250 to 500 Ryd [29]. We used a generalized gradient approximation (GGA) proposed by Perdew, Burke, and Ernzerhof to the exchangecorrelation functional [36]. An electronic temperature of 300 K was used to count the number of electrons by the FermiDirac function for all the systems we considered. For all the calculations, $\delta ={10}^{12}$ was used in the window function of Eq. (4), and ${N}_{\mathrm{band}}$ was set to be equivalent to the number of basis functions in both the summation of Eqs. (2) and (5).
Element  Basis functions  Valence states in PP 

H  ${}^{*}$H7.0s2p2d1  $1s$ 
C  ${}^{*}$C6.0s3p2d2  $2s$, $2p$ 
N  ${}^{*}$N6.0s3p2d2  $2s$, $2p$ 
Na  ${}^{*}$Na9.0s3p3d2f2  $2s$, $2p$, $3s$ 
Si  Si7.0s2p2d2f1  $3s$, $3p$ 
Cl  ${}^{*}$Cl9.0s3p3d2f2  $3s$, $3p$ 
S  S7.0s3p2d2f1  $3s$, $3p$ 
Cu  Cu6.0s3p3d3f1  $3s$, $3p$, $3d$, $4s$ 
Ga  Ga7.0s3p2d2f1  $3d$, $4s$, $4p$ 
As  As7.0s3p2d2f1  $3d$, $4s$, $4p$ 
Se  Se7.0s3p2d2  $4s$, $4p$ 
Bi  Bi8.0s3p2d2f1  $5d$, $6s$, $6p$ 
III.2 Choice of Guiding Functions $\{Q\}$
Depending on the targeted eigenstates in the window function of Eq. (4), one can choose either atomic orbitals, hybrid atomic orbitals, or MOs as the guiding functions $\{Q\}$ in Eq. (3). In this subsection we discuss the three kinds of the localized guiding functions $\{Q\}$, i.e., atomic orbitals, hybrid atomic orbitals, and embedded MOs in molecules and bulks, and how they can be generated in our implementation.
In our implementation, the Bloch function ${\varphi}_{\mathbf{k}\mu}$ is expanded by PAOs $\chi $ [33, 34] as
${\varphi}_{\mathbf{k}\mu}\u27e9={\displaystyle \sum _{\mathbf{R}}}{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \mathbf{R}}{\displaystyle \sum _{i\alpha}}{c}_{\mathbf{k}\mu ,i\alpha}{\chi}_{\mathbf{R}i\alpha}\u27e9,$  (19) 
where $i$ and $\alpha $ are atomic and orbital indices, and $c$ is a linear combination of PAO (LCPAO) coefficient. Also note that $\u27e8\mathbf{r}{\chi}_{\mathbf{R}i\alpha}\u27e9\equiv {\chi}_{i\alpha}(\mathbf{r}{\tau}_{i}\mathbf{R})$, where ${\tau}_{i}$ is the position of atom $i$. We use PAOs as the guiding functions $\{Q\}$ of atomic orbitals, corresponding to valence orbitals or specific orbitals such as localized $d$orbitals in a narrow energy window.
The choice of the atomic orbitals gives a good guiding function, however, the resultant CWFs may not recover the symmetry of CWFs respecting the bond direction to the neighboring atoms unlike hybrid atomic orbitals due to the closest property of CWFs. In this case, it would be better to use the hybrid atomic orbitals rather than the atomic orbitals as explained below. Let us introduce a projection operator for the occupied space by
$\widehat{P}={\displaystyle \frac{1}{{N}_{\mathrm{BC}}}}{\displaystyle \sum _{\mathbf{k}\mu}}{\varphi}_{\mathbf{k}\mu}\u27e9f({\epsilon}_{\mathbf{k}\mu})\u27e8{\varphi}_{\mathbf{k}\mu},$  (20) 
where $f$ is the FermiDirac function. Considering that $\{\chi \}$ are a set of nonorthogonal basis set [37], where the overlap integral is denoted by ${s}_{\mathbf{R}i\alpha ,{\mathbf{R}}^{\prime}j\beta}\equiv \u27e8{\chi}_{\mathbf{R}i\alpha}{\chi}_{{\mathbf{R}}^{\prime}j\beta}\u27e9$, the density matrix is calculated with the projection operator by
${\rho}_{\mathbf{R}i\alpha ,{\mathbf{R}}^{\prime}j\beta}$  $=$  $\sum _{\mathbf{k}\mu}}\u27e8{\stackrel{~}{\chi}}_{\mathbf{R}i\alpha}\widehat{P}{\stackrel{~}{\chi}}_{{\mathbf{R}}^{\prime}j\beta}\u27e9,$  (21)  
$=$  $\frac{1}{{N}_{\mathrm{BC}}}}{\displaystyle \sum _{\mathbf{k}\mu}}{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \left(\mathbf{R}{\mathbf{R}}^{\prime}\right)}f({\epsilon}_{\mathbf{k}\mu}){c}_{\mathbf{k}\mu ,i\alpha}{c}_{\mathbf{k}\mu ,j\beta}^{*},$ 
where $\stackrel{~}{\chi}$ is the dual orbital defined by
${\stackrel{~}{\chi}}_{\mathbf{R}i\alpha}\u27e9$  $=$  $\sum _{{\mathbf{R}}^{\prime}j\beta}}{\chi}_{{\mathbf{R}}^{\prime}j\beta}\u27e9{s}_{{\mathbf{R}}^{\prime}j\beta ,\mathbf{R}i\alpha}^{1$  (22) 
holding the following orthonormality:
$\u27e8{\chi}_{\mathbf{R}i\alpha}{\stackrel{~}{\chi}}_{{\mathbf{R}}^{\prime}j\beta}\u27e9$  $=$  $\u27e8{\stackrel{~}{\chi}}_{\mathbf{R}i\alpha}{\chi}_{{\mathbf{R}}^{\prime}j\beta}\u27e9={\delta}_{{\mathrm{\mathbf{R}\mathbf{R}}}^{\prime}}{\delta}_{ij}{\delta}_{\alpha \beta}.$  (23) 
By setting $\mathbf{R}={\mathbf{R}}^{\prime}=\mathrm{\U0001d7ce}$ and $i=j$ in Eq. (21), and diagonalizing the diagonal block element consisting of ${\rho}_{\mathrm{\U0001d7ce}i\alpha ,\mathrm{\U0001d7ce}i\beta}$, which is associated with selected atomic orbitals on the atomic site $i$, we obtain the hybrid atomic orbitals respecting the symmetry of the bond direction to the neighboring atoms. We employ the hybrid atomic orbitals as the localized functions $\{Q\}$ in our implementation. When the same atomic orbitals are chosen to form the diagonal block element as for the case of atomic orbitals $\{Q\}$, the resultant CWFs span the same subspace in the Hilbert space.
For some systems the electronic structures can be rather easily understood by employing MOs as building blocks. A molecular crystal is such a case, where the eigenstates near the chemical potential can be approximately constructed by a linear combination of the MOs of constituting molecules. Another example is the bond in molecules and bulks. The bond between two atoms embedded in the system might be analyzed by MOs associated with the two atoms. Here we show how embedded MOs in molecules and bulks can be calculated in a simple procedure. Let us start by noting that using Eq. (20) and assuming that the Bloch function is expressed by Eq. (19) the total number of electrons for a spin degenerate case is given by [37]
${N}_{\mathrm{ele}}$  $=$  $2\mathrm{t}\mathrm{r}\left[\widehat{P}\right],$  (24)  
$=$  $2{\displaystyle \sum _{\mathbf{R}i\alpha}}\u27e8{\stackrel{~}{\chi}}_{\mathbf{R}i\alpha}\widehat{P}{\chi}_{\mathbf{R}i\alpha}\u27e9,$ 
In the second line of Eq. (24) we used the following identity operator:
$\widehat{I}={\displaystyle \sum _{\mathbf{R}i\alpha}}{\chi}_{\mathbf{R}i\alpha}\u27e9\u27e8{\stackrel{~}{\chi}}_{\mathbf{R}i\alpha}={\displaystyle \sum _{\mathbf{R}i\alpha}}{\stackrel{~}{\chi}}_{\mathbf{R}i\alpha}\u27e9\u27e8{\chi}_{\mathbf{R}i\alpha}.$  (25) 
Since ${N}_{\mathrm{ele}}={N}_{\mathrm{BC}}{N}_{\mathrm{ele}}^{(0)}$ with ${N}_{\mathrm{ele}}^{(0)}=2{\sum}_{i\alpha}\u27e8{\stackrel{~}{\chi}}_{\mathrm{\U0001d7ce}i\alpha}\widehat{P}{\chi}_{\mathrm{\U0001d7ce}i\alpha}\u27e9$, it is enough to consider ${N}_{\mathrm{ele}}^{(0)}$. We utilize the formula for ${N}_{\mathrm{ele}}^{(0)}$ to calculate embedded MOs in molecules and bulks, and rewrite it the sum of partial traces as
${N}_{\mathrm{ele}}^{(0)}=2{\displaystyle \sum _{g}}{\mathrm{tr}}_{g}\left[({\stackrel{~}{\chi}}_{\mathrm{\U0001d7ce}g}\widehat{P}{\chi}_{\mathrm{\U0001d7ce}g})\right]=2{\displaystyle \sum _{g}}{\mathrm{tr}}_{g}\left[{\mathrm{\Lambda}}_{g}\right],$  (26) 
where $g$ is the index of group including PAOs on grouped atoms, e.g., a set of PAOs on a molecule. The notation of ${\chi}_{\mathrm{\U0001d7ce}g})$ stands for a subset of PAOs as
${\chi}_{\mathrm{\U0001d7ce}g})=(\mathrm{\cdots},{\chi}_{\mathrm{\U0001d7ce}i1}\u27e9,{\chi}_{\mathrm{\U0001d7ce}i2}\u27e9,\mathrm{\cdots},),$  (27) 
where PAOs on all the atoms in the group $g$ are included in the subset. The notation of ${\stackrel{~}{\chi}}_{\mathrm{\U0001d7ce}g})$ is the counterpart for the dual orbitals. Using the identity operator of Eq. (25), ${\mathrm{\Lambda}}_{g}$ in Eq. (26) is given by
${\mathrm{\Lambda}}_{g}={\displaystyle \sum _{\mathbf{R}{g}^{\prime}}}{\rho}_{\mathrm{\U0001d7ce}g,\mathbf{R}{g}^{\prime}}{s}_{\mathbf{R}{g}^{\prime},\mathrm{\U0001d7ce}g}$  (28) 
with definition of block elements:
${\rho}_{\mathbf{R}g,{\mathbf{R}}^{\prime}{g}^{\prime}}$  $=$  $({\stackrel{~}{\chi}}_{\mathbf{R}g}\widehat{P}{\stackrel{~}{\chi}}_{{\mathbf{R}}^{\prime}{g}^{\prime}}),$  (29)  
${s}_{\mathbf{R}g,{\mathbf{R}}^{\prime}{g}^{\prime}}$  $=$  $({\chi}_{\mathbf{R}g}{\chi}_{{\mathbf{R}}^{\prime}{g}^{\prime}}).$  (30) 
We now introduce the embedded MOs orbitals in molecules and bulks by a rightsingular vectors of ${\mathrm{\Lambda}}_{g}$ [38]. Since the elements of ${\mathrm{\Lambda}}_{g}$ are real in case of the collinear DFT, the rightsingular vectors can be obtained with real components. However, the SVD of ${\mathrm{\Lambda}}_{g}$ tends to produce the rightsingular vectors with complex components, which results in the complex CWFs making analysis complicated. To obtain the real rightsingular vectors, we perform the eigendecomposition of ${\mathrm{\Lambda}}_{g}^{\u2020}{\mathrm{\Lambda}}_{g}$ as
${\mathrm{\Lambda}}_{g}^{\u2020}{\mathrm{\Lambda}}_{g}={Y}_{g}{\mathrm{\Omega}}_{g}^{2}{Y}_{g}^{\u2020},$  (31) 
where ${\mathrm{\Omega}}_{g}^{2}$ is the diagonal matrix consisting of the eigenvalues of ${\mathrm{\Lambda}}_{g}^{\u2020}{\mathrm{\Lambda}}_{g}$, and ${Y}_{g}$ is a unitary matrix consisting of the corresponding eigenvectors $\{{y}_{g,\nu}\}$. By noting that the SVD of ${\mathrm{\Lambda}}_{g}$ is given by ${Z}_{g}{\mathrm{\Omega}}_{g}{Y}_{g}^{\u2020}$, we see that the rightsingular vectors ${Y}_{g}$ can be obtained by diagonalizing ${\mathrm{\Lambda}}_{g}^{\u2020}{\mathrm{\Lambda}}_{g}$ as ${Y}_{g}$. The square roots of the eigenvalues of ${\mathrm{\Lambda}}_{g}^{\u2020}{\mathrm{\Lambda}}_{g}$ are singular values of ${\mathrm{\Lambda}}_{g}$, and can be approximately regarded as the occupation of the corresponding eigenvector ${y}_{g,\nu}$. If ${\mathrm{\Lambda}}_{g}$ is a real matrix, ${Y}_{g}$ obtained by the diagonalization of ${\mathrm{\Lambda}}_{g}^{\u2020}{\mathrm{\Lambda}}_{g}$ is guaranteed to be a real unitary matrix. We further normalize the eigenvector ${y}_{g,\nu}$ by considering the overlap integrals of Eq. (30) as
${\overline{y}}_{g,\nu}\u27e9={\displaystyle \frac{{y}_{g,\nu}\u27e9}{\sqrt{\u27e8{y}_{g,\nu}{s}_{\mathrm{\U0001d7ce}g,\mathrm{\U0001d7ce}g}{y}_{g,\nu}\u27e9}}},$  (32) 
and calculate an expectation value of the KS Hamiltonian ${\widehat{H}}_{\mathrm{KS}}$ with ${\overline{y}}_{g,\nu}$ as
${\u03f5}_{g,\nu}=\u27e8{\overline{y}}_{g,\nu}{\widehat{H}}_{\mathrm{KS}}{\overline{y}}_{g,\nu}\u27e9.$  (33) 
After $\{{\overline{y}}_{g,\nu}\}$ are ordered based on the expectation values, we employ selected ones among $\{{\overline{y}}_{g,\nu}\}$, e.g., ones near the chemical potential, as the guiding functions $\{Q\}$ in Eq. (3) by monitoring the expectation values and the corresponding singular values. They are what we call embedded MOs in molecules and bulks. The reason why the expectation values are used to select the embedded MOs in addition to the singular values is that the singular values are found to be almost similar to each other as long as they are associated with occupied states. As demonstrated in Sec. IV, $\{\overline{y}\}$ work as the good guiding functions to capture the electronic structure of a molecular crystal.
IV NUMERICAL RESULTS
IV.1 Wannier Interpolated Bands
We present five numerical examples to demonstrate the broad applicability of the proposed method across various systems. In addition, a series of benchmark calculations for 30 systems are shown as supplementary information on a website [39]. In Figs. 2 (a) and (b), the interpolated bands of Si in the diamond structure, calculated by the TB Hamiltonian derived from CWFs, are compared with those by the conventional calculation. For comparison, the result with the projectiononly WFs is also shown in Fig. 2 (c), calculated by the projection method in Ref. [8]. The hybrid atomic orbitals consisting of valence minimal orbitals, which are one $s$ and a set of ${p}_{x}$, ${p}_{y}$, and ${p}_{z}$orbitals per Si, were used as the guiding functions $\{Q\}$. For the window function of Eq. (4) we use 15.0 and 0.0, relative to the chemical potential, in eV for ${\epsilon}_{0}$ and ${\epsilon}_{1}$, respectively, which covers the energy range of valence bands. On the other hand, two ${k}_{\mathrm{B}}T$ of 0.01 and 3.0 eV were used to check how conduction bands are reproduced depending on the effect of smearing. As clearly seen, the two cases well reproduce the valence bands, and the conduction bands are also reproduced well in the case (b) of ${k}_{\mathrm{B}}T=3.0$ eV. We also note that the interpolated bands by the projectiononly WFs reproduce well the reference as shown in Fig. 2 (c). Evaluating the DM function per CWF using Eq. (14), we obtain values of 0.571 and 0.444 for ${k}_{\mathrm{B}}T$ of 0.01 and 3.0 eV, respectively. For ${k}_{\mathrm{B}}T=0.01$ eV, the first four singular values at each $k$ point are close to 1, while the remaining four approach $\delta \phantom{\rule{veryverythickmathspace}{0ex}}(={10}^{12})$. Conversely, for ${k}_{\mathrm{B}}T=3.0$ eV, the eight singular values at each $k$ point distribute ranging from 1.2 to 0.1. These differences in the distribution of singular values explain why the value of the DM function for ${k}_{\mathrm{B}}T=0.01$ eV is larger than that for ${k}_{\mathrm{B}}T=3.0$ eV. One of the obtained CWFs is shown in Fig. 2 (d), which represents a $p$like CWF. Although the proposed method does not impose the localization of WFs, the decaying property of TB parameters $t$ exhibits a subexponential behavior as shown in Fig. 2 (e). In case of ${k}_{\mathrm{B}}T=0.01$ eV, the conduction bands are treated less importantly due to heavily weighting to the valence bands, resulting in the poor reproduction of the conduction bands. At first glance, one may consider that this is an undesirable feature in such a treatment. However, the feature enables us to generate the minimal atomiclike orthogonal orbitals, which well span the subspace by the valence bands, and to calculate an effective charge allocated to each atom using the minimal atomiclike orthogonal orbitals in an unbiased way. We will demonstrate the calculation of the effective atomic charges and stress the usefulness of the method later on.
The disentanglement of bands in metals is demonstrated for Cu in the FCC structure as shown in Figs. 3 (a) and (b). For comparison, the result with the projectiononly WFs is also shown in Fig. 3 (c). By selecting the parameters ${\epsilon}_{0}$ and ${\epsilon}_{1}$ so that the $3d$ bands are included, and using only five $d$orbitals as the guiding functions $\{Q\}$, the five $d$bands are reproduced as shown in Fig. 3 (a). We see that the $3d$bands are properly disentangled from the dispersive $4s$ band, and one of the obtained CWFs preserves the shape of ${d}_{{z}^{2}}$orbital as shown in Fig. 3 (d) in accordance with the subexponential decaying property of the TB parameters as shown in Fig. 3 (e). When $3d$, $4s$, and $4p$orbitals are used as the guiding functions $\{Q\}$, a broad range of bands are reproduced as shown in Fig. 3 (b). Since our PP of Cu includes $3s$, $3p$, $3d$, and $4s$ states, the original band structure has the deep $3s$ and $3p$bands. We set the parameters in the window function so as to discard the $3s$ and $3p$bands, and use the $3d$, $4s$, and $4p$orbitals, not $3s$ and $3p$orbitals, as the guiding functions $\{Q\}$. So, the $3s$ and $3p$bands are not included in the interpolated bands (not shown). The examples show that the CWFs can be flexibly and easily constructed in accordance with the purpose of study without numerical difficulties. On the other hand, the interpolated bands by the projectiononly Wannier functions (WFs) demonstrate good agreement with the reference, as shown in Fig. 3 (c). This agreement is achieved through the minimization of ${\mathrm{\Omega}}_{\mathrm{I}}$, which is originally called gauge invariant part of the spread function $\mathrm{\Omega}$, as discussed in Ref. [8]. Without minimizing ${\mathrm{\Omega}}_{\mathrm{I}}$, however, the bands obtained by the projectiononly WFs lose smoothness and exhibit slight deviations from the reference (not shown) as discussed in Ref. [3]. Our method achieves good agreement with the reference without any iterative calculations, as illustrated in Fig. 3 (a). This remarkable property can be attributed to the use of the smoothly varying window function of Eq. (4).
The interpolated bands for a molecular crystal, consisting of tetrathiafulvalene (TTF) and tetracyanoquinodimethane (TCNQ), are shown in Figs. 4 (a) and (b). By employing hybrid atomic orbitals as the guiding functions $\{Q\}$, the wide range of bands are reproduced as shown in Fig. 4 (a). On the other hand, as shown in Fig. 4 (b), only four bands near the chemical potential are reproduced using the embedded MOs in the bulk as the guiding functions $\{Q\}$. The MOs were calculated by grouping the TTF and TCNQ molecules separately as explained in Sec. III B. The expectation values calculated by Eq. (33) and the singular values $\omega $ of ${\mathrm{\Lambda}}_{g}$ are listed in Table 2. It can be seen that the singular values largely change from ${\omega}_{26}$ to ${\omega}_{27}$ and from ${\omega}_{37}$ to ${\omega}_{38}$ for TTF and TCNQ, respectively. So, the 26th MO and the 37th MO for TTF and TCNQ, respectively, were used as the guiding functions $\{Q\}$. Since there are two TTF molecules and two TCNQ molecules in the unit cell, we have the four embedded MOs as $\{Q\}$, resulting in the four bands. In Figs. 4 (c) and (d) two of the resultant CWFs are shown, which localize in TTF and TCNQ molecules, respectively, as expected. The value of the DM function per CWF is found to be 0.022, which implies that the guiding MOs are very close to the resultant CWFs. The example demonstrates that the proposed method provides a direct way to calculate CWFs for targeted bands in molecular crystals.
TTF  TCNQ  

index  $\u03f5$  $\omega $  index  $\u03f5$  $\omega $ 
23  0.069  1.359  34  0.135  1.082 
24  0.719  1.346  35  0.375  1.329 
25  0.761  1.449  36  0.477  1.336 
26  3.408  0.864  37  4.928  0.535 
27  5.947  0.047  38  9.389  0.036 
28  6.197  0.047  39  9.795  0.019 
We extend the proposed method to the noncollinear DFT [42, 43] with spinorbit interaction (SOI) [32]. The KS orbitals are expressed by twocomponent spinor, and the SOI is introduced by fully relativistic $j$dependent PPs [28]. The theoretical framework is readily extended to the noncollinear DFT with the SOI without any difficulty. However, it should be noted that the hybrid atomic orbitals and embedded MOs in molecules and bulks as the guiding functions $\{Q\}$ can be complex functions in the case. In Figs. 5 (a) and (b), we show the Wannier interpolated bands of Bi${}_{2}$Se${}_{3}$ which is known to be a topological insulator when the SOI is included in the DFT calculation. For the cases without and with the SOI, corresponding to Figs. 5 (a) and (b), respectively, it is confirmed that the interpolated bands accurately reproduce the original ones regardless of existence of the band inversion, demonstrating a wide variety of applicability of the proposed method.
Finally, we show results for GaAs in Figs. 6 (a)(d) to demonstrate how the window function of Eq. (4) with the two temperature scheme plays a significant role in reproducing the targeted bands. The interpolated bands calculated with the projectiononly WFs calculated by the projection method in Ref. [8] reproduce well the reference ones as shown in Fig. 6 (a). When the four bands located just below the chemical potential, which mainly consists of Ga$4s$ and As$4p$ orbitals, are focused by discarding the other bands with ${k}_{\mathrm{B}}{T}_{0}={k}_{\mathrm{B}}{T}_{1}$ = 0.001 eV, interpolated bands are obtained as shown in Fig. 6 (b). Other the other hand, including outer bands with ${k}_{\mathrm{B}}{T}_{0}={k}_{\mathrm{B}}{T}_{1}$ = 1.0 eV leads to an erratic behavior as shown in Fig. 6 (c), where the CWFs attempt to reproduce the lower band located around 12 eV and low lying unoccupied bands depending on the kpoint, resulting in the erratic behavior. Noting that the three bands in between 7 and 0 eV consists of mainly Ga$4s$ and As$4p$ orbitals and the low lying unoccupied bands are largely contributed by Ga$4p$ orbitals, the band around 12 eV needs to be completely discarded and the low lying unoccupied bands should be partly included when Ga$4s$, $4p$, and As$4p$ orbitals are used as the guiding functions. The result in Fig. 6 (d) is the case, where ${k}_{\mathrm{B}}{T}_{0}$ = 0.001 and ${k}_{\mathrm{B}}{T}_{1}$ = 1.0 eV were used. We see that the interpolated bands fully reproduce the reference ones as expected, which clearly demonstrates the usefulness of the two temperature scheme.
IV.2 Effective Atomic Charges
From the construction, one can calculate CWFs closest to (hybrid) valence atomic orbitals, while fully respecting the occupied states. The feature of CWFs provides us a unique way to evaluate effective atomic charges based on CWFs being orthogonal atomiclike orbitals. Since the CWFs are a set of orthonormal functions, the effective atomic charge ${\kappa}_{i}$ of atom $i$ is easily calculated without arbitrariness as
${\kappa}_{i}$  $=$  ${N}_{i}^{(\mathrm{val})}2{\displaystyle \sum _{p\in \mathrm{atom}\mathrm{\_}i}}\u27e8{W}_{\mathrm{\U0001d7ce}p}\widehat{P}{W}_{\mathrm{\U0001d7ce}p}\u27e9,$  (34)  
$=$  ${N}_{i}^{(\mathrm{val})}{\displaystyle \frac{2}{{N}_{\mathrm{BC}}}}{\displaystyle \sum _{p\in \mathrm{atom}\mathrm{\_}i}}{\displaystyle \sum _{\mathbf{k},\mu}}f({\epsilon}_{\mathbf{k}\mu}){u}_{\mathbf{k}\mu ,p}^{*}{u}_{\mathbf{k}\mu ,p},$ 
where the factor of 2 is due to the spin degeneracy, $\widehat{P}$ is the projection operator defined by Eq. (20), and ${N}_{i}^{(\mathrm{val})}$ is the number of valence electrons in the PP of atom $i$. Though the summation over orbitals $p$ belonging to atom $i$ is taken into account in Eq. (34), one can analyze the orbital resolved charges as well. Table III shows effective atomic charges in a HCN molecule and the NaCl bulk, calculated by Eq. (34) together with effective charges calculated by MLWFs [3, 7], the projectiononly WFs (POWF) [8], the Mulliken population analysis [45], and the Bader analysis [46, 47].
HCN  CWF  MLWF  POWF  MP  Bader  
Basis  H  C  N  H  C  N  H  C  N  H  C  N  H  C  N 
$A$6.0s1p1  0.077  0.052  0.025  0.665  0.303  0.969  0.262  0.215  0.047  0.384  0.164  0.221  1.000  1.746  2.746 
$A$6.0s2p2  0.069  0.003  0.073  0.793  0.067  0.860  0.312  0.009  0.303  0.070  0.321  0.252  0.201  2.439  2.639 
$A$6.0s2p2d1  0.066  0.009  0.076  0.793  0.102  0.907  0.307  0.010  0.317  0.008  0.393  0.385  0.183  2.523  2.706 
$A$6.0s3p3d2  0.066  0.008  0.074  0.699  0.343  1.042  0.314  0.001  0.315  0.110  0.298  0.408  0.182  2.539  2.721 
$A$6.0s3p3d2f2  0.066  0.008  0.074  0.721  0.344  1.065  0.308  0.007  0.315  0.167  0.297  0.464  0.182  2.532  2.715 
NaCl  CWF  MLWF  POWF  MP  Bader  
Basis  Na  Cl  Na  Cl  Na  Cl  Na  Cl  Na  Cl  
$A$9.0s2p1  0.391  0.391  0.760  0.760  0.627  0.627  0.716  0.716  0.858  0.858  
$A$9.0s3p2  0.422  0.422  0.628  0.628  0.698  0.698  0.595  0.595  0.865  0.865  
$A$9.0s3p3d2  0.421  0.421  0.876  0.876  0.697  0.697  0.158  0.158  0.853  0.853  
$A$9.0s3p3d2f2  0.421  0.421  0.877  0.877  0.697  0.697  0.062  0.062  0.850  0.850 
Both the systems are known to be notorious cases because of the difficulty in calculating the effective charges [20]. We see that the effective charges by the CWFs quickly converge as a function of basis functions, while those by the Mulliken population are highly dependent on the choice of the basis function [48]. The nonionic nature of the CN bond in the HCN molecule, obtained by the CWF method, is consistent with the picture given by the intrinsic atomic orbital method (IAO) [20]. On the other hand, it turns out that the effective charges calculated by MLWFs are relatively sensitive to the choice of basis functions compared to the cases of CWF, and exhibit large charge transfer in HCN and NaCl. The resultant ionic CN bond in the MLWF method differs from the picture derived by the IAO and CWF methods. The overestimation of the charge transfer in the MLWF method can be understood by comparing the WFs as shown in Fig. 7. It is confirmed that the CWFs (a)(i) well preserve the shape of the atomic orbitals as expected. By contrast, the MLWFs (j), (k), (l), (o), and (p) deform significantly from the shape of atomic orbitals as a result of the maximally localization of the spread function. The MLWFs (j) and (k), originating from H$1s$ and C$2s$ orbitals, respectively, reduces the occupations by the deformation, while the occupations of the MLWFs (l), (o), and (p) increase due to the deformation. The $\sigma $ bonds between HC and CN are contributed by the MLWFs (l) and (p), respectively, which originate from C2${p}_{x}$ and N2${p}_{x}$ orbitals, respectively. Thus, the MLWF (j), originating from H$1s$ orbital, does not largely contribute to the $\sigma $ bonds between HC, reducing the population on the H atom. Similarly, the population on the N atom increases, since the $\sigma $ bond between CN is formed by the MLWF (p), originating from N$2{p}_{x}$. It is apparent that attributing the occupation of the MLWFs (l) and (p) to a single atom is not justified due to the fact that they are bondcentered. This implies that MLWFs should not be used to decompose electron population to each atom as demonstrated for the HCN molecule. On the other hand, CWFs are closest to atomic orbitals when the atomic orbitals are used as the guiding functions, and the population calculated by the CWFs can be naturally attributed to each atom, which allows us to calculate the atomic populations in a physically convincing manner.
For comparison, we also calculated the effective charges using the projectiononly Wannier functions (POWF) [8] and found that the obtained values are relatively sensitive to the choice of basis functions and lie between the results obtained by the CWF and MLWF methods. Here, we analyze why the results from the POWF deviate from those by the CWF even when the inner window scheme [8] was used to focus on the valence states. In the HCN molecule, the projection of the $s$orbital for hydrogen and the $s$, ${p}_{x}$, ${p}_{y}$, ${p}_{z}$orbitals for carbon and nitrogen results in nine Wannier functions. Since there are five occupied orbitals, applying the inner window scheme to the energy range of these orbitals first selects a subspace spanned by them. Subsequently, four vectors are chosen from the subspace consisting of unoccupied orbitals so that they can be orthogonal to the former subspace and have maximal overlap with the subspace obtained by projection of guiding functions onto a subspace spanned by states included in the outer window. Once the latter subspace is selected, each subspace is given equal weight in constructing the projection matrix. As a result, the nine Wannier functions generated are not specifically designed to primarily accommodate the occupied orbitals. This statement can be verified using the CWF method: By setting the window function with ${\epsilon}_{0}=22.0$ eV, ${\epsilon}_{1}=8.0$ eV, and ${k}_{\mathrm{B}}{T}_{0}={k}_{\mathrm{B}}{T}_{1}=0.1$ eV in the CWF method, and calculating the effective charges, contributions from lowlying unoccupied orbitals are included equally. We confirmed that this calculation results in effective charges nearly equivalent to those obtained by the projectiononly method (not shown). It is evident that even when using the inner window scheme with the projectiononly approach, properly computing effective charges remains a challenge. To calculate the effective charges assigned to atoms, it is reasonable to use atomic orbitals as projection orbitals. Typically, the number of atomic orbitals used exceeds the number of occupied orbitals, leading to similar challenges as observed in the case of NaCl. A simple solution to resolve this problem is to introduce a smoothly varying window function even in the inner window scheme when constructing the projection matrix. In fact, we confirmed that introducing a smoothly varying window function in the inner window scheme almost reproduces those by the CWF method (not shown). However, one should note that the use of the smoothly varying window function in the CWF is a much simpler method.
The Bader analysis shows erratic charge transfer in the CN bond, implying the limited applicability of the Bader analysis based on a realspace decomposition to the triple bond [20], while the effective atomic charges in NaCl are similar to those by the MLWF method.
It is considered from the comparison that the effective atomic charges obtained from the CWFs could serve as a valuable tool for analyzing electronic structures in a manner related to the approach discussed in Ref. [20]. Further investigation in this direction will be conducted in future studies.
V CONCLUSIONS
We presented a noniterative method to calculate the closest Wannier functions (CWFs) to a given set of localized guiding functions such as atomic orbitals, hybrid atomic orbitals, and molecular orbitals (MOs). We defined the distance measure (DM) function $F$ by the sum of the squared distance between the projection functions $L$ and Wannier functions $W$ in the Hilbert space, and considered minimizing the DM function. It was shown that the minimization of the DM function is achieved by the polar decomposition of the projection matrix $A$ with a window function via the singular value decomposition (SVD) in a noniterative manner. The CWFs can be uniquely constructed as long as the projection matrix $A$ is positive definite. It was also shown that the method is free from subtle choice of the gauge. The disentanglement of bands is naturally taken into account by introducing a smoothly varying window function, and including ${N}_{\mathrm{band}}$ Bloch functions to generate ${N}_{\mathrm{CWF}}$ CWFs, where ${N}_{\mathrm{CWF}}\le {N}_{\mathrm{band}}$. Even for isolated bands, the method always performs the disentanglement of bands without any iterative calculation and difficulty. The advantage of the smoothly varying window function over the inner window scheme [8] was clearly demonstrated in the disentanglement of $d$bands in copper and in the calculations of the effective charges. Also, it is noted that the implementation of the method is relatively simple, which requires mainly matrixmatrix product and SVD. We have implemented the proposed method into the OpenMX code, which is based on the density functional theory (DFT), numerical pseudoatomic orbitals (PAOs), and normconserving pseudopotentials (PPs), and introduced three types of guiding functions, i.e., atomic orbitals, hybrid atomic orbitals, and embedded MOs in molecules and bulks. The first two are easily employed from the PAOs. For the last one, we developed a method to calculate embedded MOs in molecules and bulks, which focuses on a partial trace formula of the projection operator for the occupied space and applies SVD to the partial matrix for the projection operator. The interpolated bands by tightbinding (TB) models derived from the CWFs reproduce well the targeted conventional bands of a wide variety of systems including Si, Cu, the TTFTCNQ molecular crystal, and a topological insulator of Bi${}_{2}$Se${}_{3}$. These successful reproduction of targeted bands clearly demonstrates a wide variety of applicability of the proposed method. We further show the usefulness of the proposed method in calculating effective atomic charges, implying that the CWFs closest to atomic orbitals can be used as a measure to analyze electronic structures from one system to the others. Thus, we conclude that the proposed method is an alternative way in efficiently calculating WFs, and the concept of CWFs will provide a basis for development of novel methods of analyzing electronic structures and calculating physical properties [49, 50].
Acknowledgements.
The author wishes to express gratitude to Prof. Yoshiaki Sugimoto for inspiring the development of CWFs through collaborative research. Thanks are also extended to Mr. Ryotaro Koshoji for his insightful comments on the theoretical aspect. Part of the computation in this study was carried out using the computational facility of the Institute for Solid State Physics at the University of Tokyo.Appendix A: The upper bound of the diagonal elements of the matrix $D$
A proof for the upper bound of the diagonal elements of the matrix $D$ is given in the appendix. Note that $UV$ and $BV$ in Eq. (11) are partial unitary matrices in size of ${N}_{\mathrm{band}}\times {N}_{\mathrm{CWF}}$. Since $D$ is hermitian, its eigenvalues are real. Let $\lambda $ be an eigenvalue of $D$ and $x$ be the corresponding eigenvector. Then, the following equation holds:
$\frac{1}{2}}\left({V}^{\u2020}{U}^{\u2020}BV+{V}^{\u2020}{B}^{\u2020}UV\right)x\u27e9=\lambda x\u27e9.$  (35) 
Defining $y\u27e9=UVx\u27e9$ and $z\u27e9=BVx\u27e9$, and operating $\u27e8x$ from the left side of Eq. (35), we have
$\u27e8yz\u27e9+\u27e8zy\u27e9=2\lambda \u27e8xx\u27e9.$  (36) 
Noting $\u27e8yz\u27e9+\u27e8zy\u27e9\le 2\u27e8yz\u27e9$, and using the CauchySchwarz inequality $\u27e8yz\u27e9\le \sqrt{\u27e8yy\u27e9\u27e8zz\u27e9}$, we obtain
$\u27e8yz\u27e9+\u27e8zy\u27e9\le 2\sqrt{\u27e8yy\u27e9\u27e8zz\u27e9}.$  (37) 
Since $\u27e8yy\u27e9=\u27e8zz\u27e9=\u27e8xx\u27e9$, combining Eq. (36) and Eq. (37) results in
$\lambda \u27e8xx\u27e9\le \u27e8xx\u27e9.$  (38) 
Considering $\u27e8xx\u27e9\ne 0$, the upper bound of the eigenvalues of $D$ is found to be
$\lambda \le 1.$  (39) 
Noting that the matrix $D$ can be written by $x$ and $\lambda $ as
$D={\displaystyle \sum _{\nu}}{x}_{\nu}\u27e9{\lambda}_{\nu}\u27e8{x}_{\nu},$  (40) 
the diagonal elements ${d}_{nn}$ of the matrix $D$ is given by
${d}_{nn}={\displaystyle \sum _{\nu}}{\u27e8n{x}_{\nu}\u27e9}^{2}{\lambda}_{\nu},$  (41) 
where $\u27e8n{x}_{\nu}\u27e9$ is the $n$th element in the vector ${x}_{\nu}$. Since the upper bound of $\lambda $ is unity, and $\{x\}$ forms a unitary matrix, we obtain the upper bound of the diagonal elements of the matrix $D$ as
${d}_{nn}\le 1.$  (42) 
References
 [1] G.H. Wannier, The structure of electronic excitation levels in insulating crystals, Phys. Rev. 52, 191 (1937).
 [2] W. Kohn, Analytic properties of bloch waves and wannier functions, Phys. Rev. 115, 809 (1959).
 [3] N. Marzari, A.A. Mostofi, J.R. Yates, I. Souza, and D. Vanderbilt, Maximally localized wannier functions: Theory and applications Rev. Mod. Phys. 84, 1419 (2012).
 [4] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136, B864 (1964).
 [5] D.R. Hamann and D. Vanderbilt, Maximally localized wannier functions for GW quasiparticles, Phys. Rev. B 79, 045109 (2009).
 [6] F. Lechermann, A. Georges, A. Poteryaev, S. Biermann, M. Posternak, A. Yamasaki, and O.K. Andersen, Dynamical meanfield theory using Wannier functions: A flexible route to electronic structure calculations of strongly correlated materials, Phys. Rev. B 74, 125120 (2006).
 [7] N. Marzari and D. Vanderbilt, Maximally localized generalized wannier functions for composite energy bands, Phys. Rev. B 56, 12847 (1997).
 [8] I. Souza, N. Marzari, and D. Vanderbilt, Maximally localized wannier functions for entangled energy bands, Phys. Rev. B 65, 035109 (2001).
 [9] V. Vitale, G. Pizzi, A. Marrazzo, J.R. Yates, N. Marzari, and A.A. Mostofi, Automated highthroughput wannierisation, npj Computational Materials 6, 66 (2020).
 [10] A. Damle and L. Lin, Disentanglement via entanglement: A unified method for wannier localization, Multiscale Model. Simul. 16, 1392 (2018).
 [11] J. Qiao, G. Pizzi, and N. Marzari, Automated mixing of maximally localized wannier functions into target manifolds, npj Computational Materials 9, 206 (2023).
 [12] W. Ku, H. Rosner, W.E. Pickett, and R.T. Scalettar, Insulating ferromagnetism in La${}_{4}$Ba${}_{2}$Cu${}_{2}$O${}_{10}$: An ab initio wannier function analysis, Phys. Rev. Lett. 89, 167204 (2002).
 [13] W.C. Lu, C.Z. Wang, T.L. Chan, K. Ruedenberg, and K.M. Ho, Representation of electronic structures in crystals in terms of highly localized quasiatomic minimal basis orbitals, Phys. Rev. B 70, 041101(R) (2004).
 [14] X. Qian, J. Li, L. Qi, C.Z. Wang, T.L. Chan, Y.X. Yao, K.M. Ho, and S. Yip, Quasiatomic orbitals for ab initio tightbinding analysis, Phys. Rev. B 78, 245112 (2008).
 [15] P.O. Löwdin, On the nonorthogonality problem connected with the use of atomic wave functions in the theory of molecules and crystals, J. Chem. Phys. 18, 365 (1950).
 [16] J. Des Cloizeaux, Energy bands and projection operators in a crystal: Analytic and asymptotic properties, Phys. Rev. 135, A685 (1964).
 [17] A.A. Mostofi, J.R. Yates, Y.S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Wannier90: A tool for obtaining maximallylocalised wannier functions, Comput. Phys. Commun. 178, 685 (2008).
 [18] B.C. Carlson and J.M. Keller, Orthogonalization procedures and the localization of wannier functions, Phys. Rev. 105, 102 (1957).
 [19] A.E. Reed, L.A. Curtiss, and F. Weinhold, Intermolecular interactions from a natural bond orbital, donoracceptor viewpoint, Chem. Rev. 88, 899 (1988).
 [20] G. Knizia, Intrinsic atomic orbitals: An unbiased bridge between quantum theory and chemical concepts, J. Chem. Theory Comput. 9, 4834 (2013).
 [21] M. Engel, M. Marsman, C. Franchini, and G. Kresse, Electronphonon interactions using the projector augmentedwave method and wannier functions, Phys. Rev. B 101, 184302 (2020).
 [22] A.A. Mostofi, J.R. Yates, G. Pizzi, Y.S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, An updated version of wannier90: A tool for obtaining maximallylocalised wannier functions, Comput. Phys. Commun. 185, 2309 (2014).
 [23] T. Schäfer, F. Libisch, G. Kresse, A. Grüneis, Local embedding of coupled cluster theory into the random phase approximation using plane waves, J. Chem. Phys. 154, 011101 (2021).
 [24] M. Engel, H. Miranda, L. Chaput, A. Togo, C. Verdi, M. Marsman, and G. Kresse, Zeropoint renormalization of the band gap of semiconductors and insulators using the projector augmented wave method, Phys. Rev. B 106, 094316 (2022).
 [25] W. Kohn and L. J. Sham, Selfconsistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965).
 [26] Matrix Computations (3rd ed.). Baltimore: The Johns Hopkins University Press (1996), G. Golub and C.F Van Loan.
 [27] K. Fan and A. J. Hoffman, Some metric inequalities in the space of matrices, Proc. Amer. Math. Soc. 6, 111 (1955).
 [28] The code, OpenMX, pseudoatomic basis functions, and pseudopotentials are available under terms of the GNUGPL on a website (http://www.openmxsquare.org/).
 [29] T. Ozaki and H. Kino, Efficient projector expansion for the ab initio LCAO method, Phys. Rev. B 72, 045121 (2005).
 [30] T.V.T. Duy and T. Ozaki, A threedimensional domain decomposition method for largescale DFT electronic structure calculations, Comput. Phys. Commun. 185, 777 (2014).
 [31] I. Morrison, D.M. Bylander, and L. Kleinman, Nonlocal hermitian normconserving vanderbilt pseudopotential, Phys. Rev. B 47, 6728 (1993).
 [32] G. Theurich and N.A. Hill, Selfconsistent treatment of spinorbit coupling in solids using relativistic fully separable ab initio pseudopotentials, Phys. Rev. B 64, 073106 (2001).
 [33] T. Ozaki, Variationally optimized atomic orbitals for largescale electronic structures, Phys. Rev. B. 67, 155108, (2003).
 [34] T. Ozaki and H. Kino, Numerical atomic basis orbitals from H to Kr, Phys. Rev. B 69, 195113 (2004).
 [35] K. Lejaeghere, G. Bihlmayer, T. Björkman, P. Blaha, S. Blügel, V. Blum, D. Caliste, I.E. Castelli, S.J. Clark, A. Dal Corso, S. de Gironcoli, T. Deutsch, J.K. Dewhurst, I. Di Marco, C. Draxl, M. Dułak, O. Eriksson, J.A. FloresLivas, K.F. Garrity, L. Genovese, P. Giannozzi, M. Giantomassi, S. Goedecker, X. Gonze, O. Grånäs, E.K. Gross, A. Gulans, F. Gygi, D.R. Hamann, P.J. Hasnip, N.A. Holzwarth, D. Iu̧san, D.B. Jochym, F. Jollet, D. Jones, G. Kresse, K. Koepernik, E. Küçü kbenli, Y.O. Kvashnin, I.L. Locht, S. Lubeck, M. Marsman, N. Marzari, U. Nitzsche, L. Nordströ m, T. Ozaki, L. Paulatto, C.J. Pickard, W. Poelmans, M.I. Probert, K. Refson, M. Richter, G.M. Rignanese, S. Saha, M. Scheffler, M. Schlipf, K. Schwarz, S. Sharma, F. Tavazza, P. Thunströ m, A. Tkatchenko, M. Torrent, D. Vanderbilt, M.J. van Setten, V. Van Speybroeck, J.M. Wills, J.R. Yates, G.X. Zhang, and S. Cottenier, Reproducibility in density functional theory calculations of solids, Science 351, aad3000 (2016).
 [36] J.P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
 [37] T. Ozaki, M. Fukuda, and G. Jiang, Efficient O() divideconquer method with localized singleparticle natural orbitals, Phys. Rev. B 98, 245137 (2018).
 [38] As discussed in Ref. [37], strictly speaking the right eigenvectors obtained by the eigendecomposition of ${\mathrm{\Lambda}}_{g}$ correspond to embedded MOs, since the eigenvalues are directly related to the occupations. However, the right eigenvectors can be complex even for the collinear case. Thus, we use the SVD which allows us to obtain real vectors at least within the collinear case. We confirmed from numerical experiments that the right singular vectors are qualitatively similar to the right eigenvectors.
 [39] A series of benchmark results are available on a website (http://www.openmxsquare.org/cwf/).
 [40] Y.T. Lee and T. Ozaki, OpenMX Viewer: A webbased crystalline and molecular graphical user interface program, J. Mol. Graph. Model. 89, 192 (2019); https://www.openmxsquare.org/viewer/
 [41] T.J. Kistenmacher, T.E. Phillips, and D.O. Cowan, The crystal structure of the 1:1 radical cationradical anion salt of 2,2’bisl,3dithiole (TTF) and 7,7,8,8tetracyanoquinodimethane (TCNQ), Acta Cryst. B30, 763 (1974).
 [42] U. von Barth and L. Hedin, A local exchangecorrelation potential for the spin polarized case. i, J. Phys. C: Solid State Phys. 5, 1629 (1972).
 [43] J. Kubler, K.H. Hock, J. Sticht, and A. R. Williams, Density functional theory of noncollinear magnetism, J. Phys. F: Met. Phys. 18, 469 (1988).
 [44] S. Nakajima, The crystal structure of Bi${}_{2}$Te${}_{3x}$Se${}_{x}$, J. Phys. Chem. Solids 24, 479 (1963).
 [45] R.S. Mulliken, Electronic population analysis on LCAO–MO molecular wave functions. I, J. Chem. Phys. 23, 1833 (1955).
 [46] Atoms in Molecules: A Quantum Theory. Oxford University Press, New York (1990), R.F.W. Bader.
 [47] G. Henkelman, A. Arnaldsson, and H. Jónsson, A fast and robust algorithm for bader decomposition of charge density, Comput. Mater. Sci. 36, 354 (2006).
 [48] R.S. Mulliken, Iodine revisited, J. Chem. Phys. 55, 288 (1971).
 [49] R. Zhang, Y. Yasui, M. Fukuda, T. Ozaki, M. Ogura, T. Makino, D. Takeuchi, and Y. Sugimoto, Atomic observation on diamond (001) surfaces with noncontact atomic force microscopy, arXiv:2403.17427.
 [50] A. Inda, R. Oiwa, S. Hayami, H.M. Yamamoto, and H. Kusunose, Quantification of chirality based on electric toroidal monopole, J. Chem. Phys. 160, 184117 (2024).