up


Closest Wannier functions to a given set of localized orbitals
The original has been published in Phys. Rev. B 110, 125115 (2024).

Taisuke Ozaki, ISSP, Univ. of Tokyo
Abstract

A non-iterative 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 TTF-TCNQ molecular crystal, and a topological insulator of Bi2Se3. 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 large-scale systems with complicated electronic structures [3]. To overcome this challenge, methods aimed at automated high-throughput 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 ill-conditioned 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 {ϕ}, which can be obtained by solving the Kohn-Sham (KS) equation [25] within the density functional theory (DFT) [4], normalized as

ϕ𝐤1μ|ϕ𝐤2ν=NBCδ𝐤1𝐤2δμν, (1)

where 𝐤 and μ are the k-vector and band index, respectively, and NBC is the number of primitive cells in the Born-von Karman (BvK) boundary condition. We consider the projection of a localized guiding function Q onto the Bloch functions {ϕ} weighted with a window function w(ε) as

|L𝐑p=1NBC𝐤μe-i𝐤𝐑|ϕ𝐤μa𝐤μ,p (2)

with

a𝐤μ,p=w(ε𝐤μ)ϕ𝐤μ|Q𝟎p, (3)

where 𝐑 and p are the translational lattice vector and index of the localized orbital, respectively, ε𝐤μ is the eigenvalue of the KS equation, and Q𝟎pQ𝐑p (𝐑=𝟎). In the summation of Eq. (2), the numbers of k-points and the Bloch functions to be included per k-point are NBC and Nband, 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𝐑p}, such as atomic orbitals, hybrid orbitals, or molecular orbitals (MOs), are assumed to localize in the unit cell specified by 𝐑. When the window function w(ε) is taken to be unity, Eq. (2) is nothing but an expansion of Q𝐑p with Bloch functions from the identity relation 1NBC𝐤μ|ϕ𝐤μϕ𝐤μ|=I^. To introduce an expansion of Q𝐑p by the Bloch functions {ϕ𝐤μ} in a subspace, we use the following window function w(ε):

w(ε)=1-exp(x0+x1)(1+exp(x0))(1+exp(x1))+δ (4)

with the definition of x0β(ε0-ε) and x1β(ε-ε1), where β=1kBT with the Boltzmann constant of kB 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+exp(x)). It is also possible to introduce two temperatures T0 and T1. The usefulness of the two temperature scheme will be demonstrated for the case of GaAs in Sec. IV. The last term δ in Eq. (4) is a small constant, e.g., 10-12, which is introduced to avoid the ill-conditioning of the matrix consisting of a𝐤μ,p by Eq. (3) as will be discussed later. As shown in Fig. 1, the parameters of ε0 and ε1 (ε0<ε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 ε0 and ε1 in w(ε). 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.


Figure 1: Window functions of Eq. (4) with kBT=0.1, 1.0, and 3.0 eV. ε0 and ε1 were set to be -5.0 and 5.0 eV, respectively, and δ=10-12 was used.

The function L𝐑p defined by Eq. (2) must be similar to the localized function Q𝐑p, but they are generally non-orthogonal to each other. We now consider generating a set of closest Wannier functions (CWFs) {W𝐑p} to a set of localized functions {L𝐑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𝐑p=1NBC𝐤μe-i𝐤𝐑|ϕ𝐤μb𝐤μ,p, (5)

and defining the residual function R:

|R𝐑p=|L𝐑p-|W𝐑p, (6)

the distance measure (DM) function F we minimize is given by

F[B] = pR𝟎p|R𝟎p, (7)
= 1NBC𝐤X[B,𝐤]

with

X[B,𝐤]=tr[(A(𝐤)-B(𝐤))(A(𝐤)-B(𝐤))], (8)

where the elements of the matrices A(𝐤) and B(𝐤) are given by a𝐤μ,p and b𝐤μ,p, respectively. The second line of Eq. (7) is obtained by considering the orthonormality of Eq. (1), and X[B,𝐤] of Eq. (8) is regarded as the squared Frobenius norm of (A(𝐤)-B(𝐤)) [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 NCWF per unit cell, being equivalent to the number of the guiding functions per cell, and NCWF should be smaller than or equal to Nband to guarantee the linear independency of the subspace spanned by the CWFs, resulting in size of Nband×NCWF for the matrices A(𝐤) and B(𝐤). Assuming B(𝐤)B(𝐤)=I, where I is of size NCWF×NCWF, 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(𝐤) is regarded as a part of a unitary matrix of size Nband×Nband, 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(𝐤) as U(𝐤) which is calculated by the polar decomposition of A(𝐤)=U(𝐤)P(𝐤), where U(𝐤) and P(𝐤) are unitary and hermitian, respectively [27]. Let us prove the statement below. The polar decomposition of A(𝐤) is obtained via SVD of A(𝐤) as

A=WΣV=WVVΣV=UP, (9)

where we dropped the dependency on 𝐤 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 Nband×NCWF and NCWF×NCWF, respectively, and Σ is the singular value diagonal matrix in size of NCWF×NCWF. Note that UWV and PVΣV. It is worth mentioning that W and U are partial unitary matrices, and hold WW=I (WWI) and UU=I (UUI), and that V is a full unitary matrix holding VV=I (VV=I). We evaluate X of Eq. (8) for each 𝐤 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] = 2tr[12(AB+BA)-P], (10)
= 2tr[ΣD-Σ],
= 2nσn(dnn-1)

with

D=12(VUBV+VBUV). (11)

The second line of Eq. (10) is derived by using the polar decomposition of A and P=VΣV, and σn and dnn in the third line are diagonal elements of Σ 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 VU (VB) and BV (UV). Another proof for the upper bound is also given based on the Cauchy-Schwarz inequality in the appendix. By considering the upper bound of the diagonal elements of D and 0σn, the third line of Eq. (10) leads to the following inequality:

X[U]X[B]. (12)

The equality of Eq. (12) holds if B=U. If some of σn is zero, the corresponding dnn 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 σn are positive, all the dnns 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 dnn are unity happens when

UB+BU=2I, (13)

which is derived by equating D=I and multiplying V and V 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 σn are positive, since F[B] of Eq. (7) consists of the sum of X[B,𝐤] over k. The uniqueness of U itself is related to that of the SVD for the matrix A of Eq. (3). If A has NCWF positive singular values which are non-degenerate, the SVD is uniquely determined except for the non-unique 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=WKKV=WV with KK=I. Further, noting that the non-unique phases in W and V are canceled out when the matrix U is computed as WV, we conclude that U is uniquely determined if A has the NCWF positive singular values. The violation from the positive definiteness of the singular values can be avoided by the small constant of δ in Eq. (4). On the other hand, if δ in Eq. (4) is set to be 0 and a narrow window, including fewer than NCWF eigenstates for a given 𝐤, is used with a small kBT, the matrix A has NPSV positive singular values, where NPSV<NCWF. In this case, a subspace of the (NCWF-NPSV) dimension is arbitrary chosen to form W and V of the dimension of NCWF in addition to the subspace of the NPSV dimension, which breaks the symmetry of CWFs. The symmetry breaking is due to the non-uniqueness of the SVD when some of the singular values are zero. The situation can be avoided by introducing a small constant δ 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 non-unique 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]=1NBC𝐤p(σ𝐤p-1)2, (14)

where σ𝐤p is a singular value of the matrix A(𝐤). From Eq. (14), we find that the mean squared distance between L and W is related to the deviation of σ 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

𝐑ei𝐤𝐑L𝟎p|L𝐑q=μa𝐤μ,p*a𝐤μ,q, (15)

which is obtained by noting that 1NBC𝐑ei(𝐤-𝐤)𝐑=δ𝐤𝐤. By writing Eq. (15) as S(𝐤)=A(𝐤)A(𝐤) in a matrix form, and using Eq. (9), we have the following relation:

S(𝐤)=A(𝐤)A(𝐤)=V(𝐤)Σ2(𝐤)V(𝐤). (16)

Comparing Eq. (16) with P(𝐤)=V(𝐤)Σ(𝐤)V(𝐤), one obtains a relation P(𝐤)=S1/2(𝐤). If P(𝐤) is invertible, the matrix U(𝐤) is given with Eq. (9) by

U(𝐤)=A(𝐤)S-1/2(𝐤). (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 ill-conditioned. 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 Nband(NCWF) Bloch functions per k-point are included in Eqs. (2) and (5), it should be emphasized that the proposed method always disentangles Nband states to generate NCWF Wannier functions. A couple of examples will be shown in Sec. IV, including valence and low-lying 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(𝐤) are obtained by the polar decomposition, tight-binding (TB) parameters are calculated as expectation values of the KS Hamiltonian H^KS by summing contributions over 𝐤 and μ as

t𝟎p,𝐑q = W𝟎p|H^KS|W𝐑q, (18)
= 1NBC𝐤με𝐤μu𝐤μ,p*u𝐤μ,qe-i𝐤𝐑,

where u𝐤μ,q is the matrix element of U(𝐤). 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. 1.

    Determining ε0, ε1, and kBT 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. 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. 3.

    Calculation of A(𝐤) by Eq. (3). The calculation of the overlap ϕ𝐤μ|Q𝟎p depends on the implementation of the KS method.

  4. 4.

    Performing the SVD of A(𝐤). A proper numerical library might be used.

  5. 5.

    Calculation of U(𝐤). Using the left and right singular matrices, W(𝐤) and V(𝐤), of A(𝐤), U(𝐤) is calculated as W(𝐤)V(𝐤).

  6. 6.

    Calculations of the CWFs and TB parameters. The CWFs can be calculated using Eq. (5) with B=U obtained by the step 5. The calculation is expressed by the sum of the matrix-matrix product over 𝐤. Also, the TB parameters are obtained by Eq. (18).

It should be stressed that the minimization of F can be efficiently performed in a non-iterative 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(NBCNbandNCWFNbasis), O(NBCNband3), O(NBCNbandNCWF2), or O(NBCNbandNCWFNbasis) for the step 3, 4, 5, or 6, respectively, where Nbasis 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(NBCNbandNCWF).

For isolated systems, the Brillouin zone sampling is limited to only the Γ 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 DFT-KS 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 norm-conserving pseudopotentials (PPs) [31, 32] and optimized pseudo-atomic 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.0-s2p2d1, 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 Slater-type 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 exchange-correlation functional [36]. An electronic temperature of 300 K was used to count the number of electrons by the Fermi-Dirac function for all the systems we considered. For all the calculations, δ=10-12 was used in the window function of Eq. (4), and Nband was set to be equivalent to the number of basis functions in both the summation of Eqs. (2) and (5).

Table 1: Basis functions and valence states included in PPs. *For the calculations of atomic charges, the basis functions listed in Table III were used.

Element Basis functions Valence states in PP
H *H7.0-s2p2d1 1s
C *C6.0-s3p2d2 2s, 2p
N *N6.0-s3p2d2 2s, 2p
Na *Na9.0-s3p3d2f2 2s, 2p, 3s
Si Si7.0-s2p2d2f1 3s, 3p
Cl *Cl9.0-s3p3d2f2 3s, 3p
S S7.0-s3p2d2f1 3s, 3p
Cu Cu6.0-s3p3d3f1 3s, 3p, 3d, 4s
Ga Ga7.0-s3p2d2f1 3d, 4s, 4p
As As7.0-s3p2d2f1 3d, 4s, 4p
Se Se7.0-s3p2d2 4s, 4p
Bi Bi8.0-s3p2d2f1 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 ϕ𝐤μ is expanded by PAOs χ [33, 34] as

|ϕ𝐤μ=𝐑ei𝐤𝐑iαc𝐤μ,iα|χ𝐑iα, (19)

where i and α are atomic and orbital indices, and c is a linear combination of PAO (LCPAO) coefficient. Also note that 𝐫|χ𝐑iαχiα(𝐫-τi-𝐑), where τ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

P^=1NBC𝐤μ|ϕ𝐤μf(ε𝐤μ)ϕ𝐤μ|, (20)

where f is the Fermi-Dirac function. Considering that {χ} are a set of non-orthogonal basis set [37], where the overlap integral is denoted by s𝐑iα,𝐑jβχ𝐑iα|χ𝐑jβ, the density matrix is calculated with the projection operator by

ρ𝐑iα,𝐑jβ = 𝐤μχ~𝐑iα|P^|χ~𝐑jβ, (21)
= 1NBC𝐤μei𝐤(𝐑-𝐑)f(ε𝐤μ)c𝐤μ,iαc𝐤μ,jβ*,

where χ~ is the dual orbital defined by

|χ~𝐑iα = 𝐑jβ|χ𝐑jβs𝐑jβ,𝐑iα-1 (22)

holding the following orthonormality:

χ𝐑iα|χ~𝐑jβ = χ~𝐑iα|χ𝐑jβ=δ𝐑𝐑δijδαβ. (23)

By setting 𝐑=𝐑=𝟎 and i=j in Eq. (21), and diagonalizing the diagonal block element consisting of ρ𝟎iα,𝟎iβ, 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]

Nele = 2tr[P^], (24)
= 2𝐑iαχ~𝐑iα|P^|χ𝐑iα,

In the second line of Eq. (24) we used the following identity operator:

I^=𝐑iα|χ𝐑iαχ~𝐑iα|=𝐑iα|χ~𝐑iαχ𝐑iα|. (25)

Since Nele=NBCNele(0) with Nele(0)=2iαχ~𝟎iα|P^|χ𝟎iα, it is enough to consider Nele(0). We utilize the formula for Nele(0) to calculate embedded MOs in molecules and bulks, and rewrite it the sum of partial traces as

Nele(0)=2gtrg[(χ~𝟎g|P^|χ𝟎g)]=2gtrg[Λg], (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 |χ𝟎g) stands for a subset of PAOs as

|χ𝟎g)=(,|χ𝟎i1,|χ𝟎i2,,), (27)

where PAOs on all the atoms in the group g are included in the subset. The notation of |χ~𝟎g) is the counterpart for the dual orbitals. Using the identity operator of Eq. (25), Λg in Eq. (26) is given by

Λg=𝐑gρ𝟎g,𝐑gs𝐑g,𝟎g (28)

with definition of block elements:

ρ𝐑g,𝐑g = (χ~𝐑g|P^|χ~𝐑g), (29)
s𝐑g,𝐑g = (χ𝐑g|χ𝐑g). (30)

We now introduce the embedded MOs orbitals in molecules and bulks by a right-singular vectors of Λg [38]. Since the elements of Λg are real in case of the collinear DFT, the right-singular vectors can be obtained with real components. However, the SVD of Λg tends to produce the right-singular vectors with complex components, which results in the complex CWFs making analysis complicated. To obtain the real right-singular vectors, we perform the eigendecomposition of ΛgΛg as

ΛgΛg=YgΩg2Yg, (31)

where Ωg2 is the diagonal matrix consisting of the eigenvalues of ΛgΛg, and Yg is a unitary matrix consisting of the corresponding eigenvectors {yg,ν}. By noting that the SVD of Λg is given by ZgΩgYg, we see that the right-singular vectors Yg can be obtained by diagonalizing ΛgΛg as Yg. The square roots of the eigenvalues of ΛgΛg are singular values of Λg, and can be approximately regarded as the occupation of the corresponding eigenvector yg,ν. If Λg is a real matrix, Yg obtained by the diagonalization of ΛgΛg is guaranteed to be a real unitary matrix. We further normalize the eigenvector yg,ν by considering the overlap integrals of Eq. (30) as

|y¯g,ν=|yg,νyg,ν|s𝟎g,𝟎g|yg,ν, (32)

and calculate an expectation value of the KS Hamiltonian H^KS with y¯g,ν as

ϵg,ν=y¯g,ν|H^KS|y¯g,ν. (33)

After {y¯g,ν} are ordered based on the expectation values, we employ selected ones among {y¯g,ν}, 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, {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 projection-only 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 px-, py-, and pz-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 ε0 and ε1, respectively, which covers the energy range of valence bands. On the other hand, two kBT 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 kBT=3.0 eV. We also note that the interpolated bands by the projection-only 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 kBT of 0.01 and 3.0 eV, respectively. For kBT=0.01 eV, the first four singular values at each k point are close to 1, while the remaining four approach δ(=10-12). Conversely, for kBT=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 kBT=0.01 eV is larger than that for kBT=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 kBT=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 atomic-like orthogonal orbitals, which well span the subspace by the valence bands, and to calculate an effective charge allocated to each atom using the minimal atomic-like 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.


Figure 2: Interpolated bands (red circles) of silicon in the diamond structure, calculated by the TB Hamiltonian derived from CWFs with (a) kBT=0.01 and (b) 3.0 eV in the window function of Eq. (4), respectively. ε0 and ε1 relative to the chemical potential were set to be -15.0 and 0.0 eV, respectively. The solid line is the original one directly calculated by the conventional scheme. The result with the projection-only WFs is shown in (c), calculated by the projection method in Ref. [8] with the outer window (-14.0 to 11.0 eV) and the inner window (-14.0 to 5.0 eV). The number of k-points for the Brillouin zone sampling was 13 × 13 × 13. The experimental lattice constant of 5.43 Å was used. The values of the DM function per CWF are 0.571 and 0.444 for (a) and (b), respectively. (d) a CWF in the case (b), where isovalues of ±0.04 (orange:0.04, blue:-0.04) were used for drawing the isosurfaces using OpenMX Viewer [40]. The same isovalues were used for the cases of Figs. 3(d), 4(c), and 4(d). (e) TB parameters t (eV) for the case (a) as a function of distance.


Figure 3: Interpolated bands (red circles) of copper in the FCC structure, calculated by the TB Hamiltonian derived from CWFs with (a) ε0=-5.5, ε1=-1.0 eV, kBT=1.0 eV, and hybrid 3d-orbitals as {Q}, and (b) ε0=-11.0, ε1=24.0 eV, kBT=3.0 eV, and hybrid 3d-, 4s-, 4p-orbitals as {Q}. The solid line is the original one directly calculated by the conventional scheme. The result with the projection-only WFs calculated with the same guiding functions in (b) is shown in (c), for which the projection method in Ref. [8] was used with the outer window (-10.0 to 0.0 eV) and the inner window (-4.0 to -2.0 eV). The number of k-points for the Brillouin zone sampling was 21 × 21 × 21. The experimental lattice constant of 3.61 Å was used. The values of the DM function per CWF are 0.081 and 0.211 for (a) and (b), respectively. (d) a CWF in the case (a). (e) TB parameters t (eV) for the case (a) as a function of distance.

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 projection-only WFs is also shown in Fig. 3 (c). By selecting the parameters ε0 and ε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 dz2-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 projection-only Wannier functions (WFs) demonstrate good agreement with the reference, as shown in Fig. 3 (c). This agreement is achieved through the minimization of ΩI, which is originally called gauge invariant part of the spread function Ω, as discussed in Ref. [8]. Without minimizing ΩI, however, the bands obtained by the projection-only 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).


Figure 4: Interpolated bands, red lines in (a) and red circles in (b), of TTF-TCNQ, calculated by the TB Hamiltonian derived from CWFs with (a) ε0=-22.0 and ε1=1.0 eV relative to the chemical potential, kBT=2.0 eV, and hybrid atomic valence orbitals as {Q}, and (b) ε0=-1.0 and ε1=1.0 eV relative to the chemical potential, kBT=0.01 eV, and the 26th MO and the 37th MO for TTF and TCNQ as {Q}. The solid line is the original one directly calculated by the conventional scheme. The number of k-points for the Brillouin zone sampling was 5 × 19 × 3. An experimental crystal structure was used [41]. The values of the DM function per CWF are 0.441 and 0.022 for (a) and (b), respectively. CWFs for (c) TTF in TTF-TCNQ, and (d) TCNQ in TTF-TCNQ, obtained from the case (b).

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 ω of Λg are listed in Table 2. It can be seen that the singular values largely change from ω26 to ω27 and from ω37 to ω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.

Table 2: Expectation values ϵ (eV) of the KS Hamiltonian calculated by Eq. (33), which is relative to the chemical potential, and singular values ω of Λg for embedded MOs in the TTF-TCNQ molecular crystal.

TTF TCNQ
index ϵ ω index ϵ ω
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 non-collinear DFT [42, 43] with spin-orbit interaction (SOI) [32]. The KS orbitals are expressed by two-component spinor, and the SOI is introduced by fully relativistic j-dependent PPs [28]. The theoretical framework is readily extended to the non-collinear 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 Bi2Se3 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.


Figure 5: Interpolated bands (red circles) of Bi2Se3 (a) without and (b) with the spin-orbit interaction, calculated by the TB Hamiltonian derived from CWFs. The solid line is the original ones directly calculated by the conventional scheme. ε0=-6.0 and ε1=3.0 eV relative to the chemical potential, kBT=1.0 eV, and hybrid Bi 6p-orbitals and Se 4p-orbitals as {Q}, were used in both the calculations. The number of k-points for the Brillouin zone sampling was 7 × 7 × 7. An experimental crystal structure was used [44]. The values of the DM function per CWF are 0.163 and 0.184 for (a) and (b), respectively.

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 projection-only 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 kBT0=kBT1 = 0.001 eV, interpolated bands are obtained as shown in Fig. 6 (b). Other the other hand, including outer bands with kBT0=kBT1 = 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 k-point, 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 kBT0 = 0.001 and kBT1 = 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.


Figure 6: Interpolated bands (red circles) of GaAs in the zinc blende structure, calculated by the TB Hamiltonian derived from (a) projection-only WFs calculated by the projection method in Ref. [8] with the outer window (-10.0 to 12.0 eV) and the inner window (-8.0 to 6.0 eV), CWFs with (b) kBT0=kBT1 = 0.001 eV, (c) kBT0=kBT1 = 1.0 eV, and (d) kBT0 = 0.001 and kBT1 = 1.0 eV, respectively. For (b), (c), and (d), ε0 and ε1 relative to the chemical potential were set to be -9.0 and 0.0 eV, respectively. Ga-4s, 4p, and As-4p orbitals were used as the guiding functions in all the cases. The solid line is the original one directly calculated by the conventional scheme. The number of k-points for the Brillouin zone sampling was 13 × 13 × 13. The experimental lattice constant of 5.65 Å was used. The value of the DM function per CWF is 0.578 in the case (d).

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 atomic-like orbitals. Since the CWFs are a set of orthonormal functions, the effective atomic charge κi of atom i is easily calculated without arbitrariness as

κi = Ni(val)-2patom_iW𝟎p|P^|W𝟎p, (34)
= Ni(val)-2NBCpatom_i𝐤,μf(ε𝐤μ)u𝐤μ,p*u𝐤μ,p,

where the factor of 2 is due to the spin degeneracy, P^ is the projection operator defined by Eq. (20), and Ni(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 projection-only WFs (POWF) [8], the Mulliken population analysis [45], and the Bader analysis [46, 47].

Table 3: Effective atomic charges in a HCN molecule and the NaCl bulk, calculated by the proposed method (CWF), the MLWFs (MLWF), the projection-only WFs (POWF) [8], the Mulliken population analysis (MP), and the Bader analysis (Bader) with a variety of basis functions, where A represents the constituting atoms. Valence atomic orbitals were used as the guiding functions {Q}. For the window function, ε0=-55.0 and ε1=0.0 eV for HCN and ε0=-35.0 and ε1=0.0 eV for NaCl, relative to the chemical potential, and kBT=0.1 eV are used. The same number of MLWFs and POWFs as the cases of CWF were generated starting from the same valence atomic orbitals as the initial guess, where the outer and inner windows, relative to the chemical potential, were set to be -22.0 to 14.0 (-30.0 to 37.0) and -22.0 to 0.0 (-29.0 to 0.0) eV for HCN (NaCl), respectively. The number of k-points for the Brillouin zone sampling is 9 × 9 × 9 for the NaCl bulk with the lattice constant of 5.63 Å.

HCN CWF MLWF POWF MP Bader
Basis H C N H C N H C N H C N H C N
A6.0-s1p1 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
A6.0-s2p2 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
A6.0-s2p2d1 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
A6.0-s3p3d2 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
A6.0-s3p3d2f2 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
A9.0-s2p1 0.391 -0.391 0.760 -0.760 0.627 -0.627 0.716 -0.716 0.858 -0.858
A9.0-s3p2 0.422 -0.422 0.628 -0.628 0.698 -0.698 0.595 -0.595 0.865 -0.865
A9.0-s3p3d2 0.421 -0.421 0.876 -0.876 0.697 -0.697 0.158 -0.158 0.853 -0.853
A9.0-s3p3d2f2 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 non-ionic 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 σ bonds between H-C and C-N are contributed by the MLWFs (l) and (p), respectively, which originate from C-2px and N-2px orbitals, respectively. Thus, the MLWF (j), originating from H-1s orbital, does not largely contribute to the σ bonds between H-C, reducing the population on the H atom. Similarly, the population on the N atom increases, since the σ bond between C-N is formed by the MLWF (p), originating from N-2px. 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 bond-centered. 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.


Figure 7: Wannier functions of a HCN molecule calculated by (a)-(i) the CWF and (j)-(r) the MLWF methods. The number that follows is the occupation number considering the degree of spin degeneracy. The white, silver, and blue spheres correspond to hydrogen, carbon, and nitrogen atoms, respectively. In all the cases, isovalues of ±0.08 (orange:0.08, blue:-0.08) are used for drawing the isosurfaces using OpenMX Viewer [40]. The computational conditions are the same as for the cases of A6.0-s3p3d2f2 in Table III.

For comparison, we also calculated the effective charges using the projection-only 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-, px-, py-, pz-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 ε0=-22.0 eV, ε1=8.0 eV, and kBT0=kBT1=0.1 eV in the CWF method, and calculating the effective charges, contributions from low-lying unoccupied orbitals are included equally. We confirmed that this calculation results in effective charges nearly equivalent to those obtained by the projection-only method (not shown). It is evident that even when using the inner window scheme with the projection-only 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 real-space 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 non-iterative 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 non-iterative 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 Nband Bloch functions to generate NCWF CWFs, where NCWFNband. 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 matrix-matrix product and SVD. We have implemented the proposed method into the OpenMX code, which is based on the density functional theory (DFT), numerical pseudo-atomic orbitals (PAOs), and norm-conserving 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 tight-binding (TB) models derived from the CWFs reproduce well the targeted conventional bands of a wide variety of systems including Si, Cu, the TTF-TCNQ molecular crystal, and a topological insulator of Bi2Se3. 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 Nband×NCWF. Since D is hermitian, its eigenvalues are real. Let λ be an eigenvalue of D and x be the corresponding eigenvector. Then, the following equation holds:

12(VUBV+VBUV)|x=λ|x. (35)

Defining |y=UV|x and |z=BV|x, and operating x| from the left side of Eq. (35), we have

y|z+z|y=2λx|x. (36)

Noting y|z+z|y2|y|z|, and using the Cauchy-Schwarz inequality |y|z|y|yz|z, we obtain

y|z+z|y2y|yz|z. (37)

Since y|y=z|z=x|x, combining Eq. (36) and Eq. (37) results in

λx|xx|x. (38)

Considering x|x0, the upper bound of the eigenvalues of D is found to be

λ1. (39)

Noting that the matrix D can be written by x and λ as

D=ν|xνλνxν|, (40)

the diagonal elements dnn of the matrix D is given by

dnn=ν|n|xν|2λν, (41)

where n|xν is the n-th element in the vector xν. Since the upper bound of λ is unity, and {x} forms a unitary matrix, we obtain the upper bound of the diagonal elements of the matrix D as

dnn1. (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 mean-field 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 high-throughput 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 La4Ba2Cu2O10: 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 tight-binding analysis, Phys. Rev. B 78, 245112 (2008).
  • [15] P.-O. Löwdin, On the non-orthogonality 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 maximally-localised 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, donor-acceptor 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, Electron-phonon interactions using the projector augmented-wave 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 maximally-localised 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, Zero-point 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, Self-consistent 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, pseudo-atomic basis functions, and pseudopotentials are available under terms of the GNU-GPL on a website (http://www.openmx-square.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 three-dimensional domain decomposition method for large-scale DFT electronic structure calculations, Comput. Phys. Commun. 185, 777 (2014).
  • [31] I. Morrison, D.M. Bylander, and L. Kleinman, Nonlocal hermitian norm-conserving vanderbilt pseudopotential, Phys. Rev. B 47, 6728 (1993).
  • [32] G. Theurich and N.A. Hill, Self-consistent treatment of spin-orbit coupling in solids using relativistic fully separable ab initio pseudopotentials, Phys. Rev. B 64, 073106 (2001).
  • [33] T. Ozaki, Variationally optimized atomic orbitals for large-scale 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. Flores-Livas, 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() divide-conquer method with localized single-particle natural orbitals, Phys. Rev. B 98, 245137 (2018).
  • [38] As discussed in Ref. [37], strictly speaking the right eigenvectors obtained by the eigendecomposition of Λ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.openmx-square.org/cwf/).
  • [40] Y.-T. Lee and T. Ozaki, OpenMX Viewer: A web-based crystalline and molecular graphical user interface program, J. Mol. Graph. Model. 89, 192 (2019); https://www.openmx-square.org/viewer/
  • [41] T.J. Kistenmacher, T.E. Phillips, and D.O. Cowan, The crystal structure of the 1:1 radical cation-radical anion salt of 2,2’-bis-l,3-dithiole (TTF) and 7,7,8,8-tetracyanoquinodimethane (TCNQ), Acta Cryst. B30, 763 (1974).
  • [42] U. von Barth and L. Hedin, A local exchange-correlation 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 non-collinear magnetism, J. Phys. F: Met. Phys. 18, 469 (1988).
  • [44] S. Nakajima, The crystal structure of Bi2Te3-xSex, 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 non-contact 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).