Construction of Maximally Localized Wannier Functions
CConstruction of Maximally Localized Wannier Functions
Junbo Zhu, Zhu Chen, and Biao Wu
1, 3, 4, 5, ∗ International Center for Quantum Materials, School of Physics, Peking University, Beijing 100871, China Institute of Applied Physics and Computational Mathematics, Beijing 100088, China Collaborative Innovation Center of Quantum Matter, Beijing 100871, China Wilczek Quantum Center, College of Science, Zhejiang University of Technology, Hangzhou 310014, China Synergetic Innovation Center for Quantum Effects and Applications (SICQEA),Hunan Normal University, Changsha 410081, China (Dated: September 30, 2016)We present a general method of constructing maximally localized Wannier functions. It consistsof three steps: (1) picking a localized trial wave function, (2) performing a full band projection,and (3) orthonormalizing with the L¨owdin method. Our method is capable of producing maximallylocalized Wannier functions without further minimization, and it can be applied straightforwardlyto random potentials without using supercells. The effectiveness of our method is demonstrated forboth simple bands and composite bands.
PACS numbers: 71.15.-m, 71.15.Ap, 72.80.Ng, 37.10.Jk
Wannier functions that are localized at lattice sites arean alternative representation of electronic states in crys-talline solid to the Bloch wave representation [1]. Theyoffer an insightful picture of chemical bonds, play a piv-otal role in the modern theory of polarization [2, 3], andare the basis for an efficient linear-scaling algorithmsin electric-structure calculations [4, 5]. Wannier func-tions are also important in linking cold atom experimentsin continuous light potentials with lattice Hamiltonians,such as the Bose-Hubbard model and Anderson randomlattice [6–8].Although there are explicit formula that transformBloch waves to Wannier functions, the construction ofWannier functions is far from trivial. The primary rea-son is that there are infinite number of Wannier func-tions for a given band due to phase choices for Blochwaves whereas in practice the maximally localized Wan-nier functions (MLWFs) are sought and preferred [9].For one dimensional lattice, Kohn showed how to fix theBloch wave phases to obtain these MLWFs [10]. In 1971,Teichler found a general method to construct Wannierfunctions, which is insensitive to the phases of Blochwaves [11]. However, this method does not guaranteemaximal localization and depends on initial trial func-tion. A more sophisticated method involving numericalminimization of the spread of a Wannier function hasbeen developed to compute MLWFs [12]. There are alsosome new developments recently [13–15].In this work we present a general method for com-puting MLWFs. Our method is somehow similar to Te-ichler’s method as it also involves projection and theL¨owdin orthonormalization method [16]. Neverthelessthere is a significant difference so that our method canproduce MLWFs and is insensitive to the initial trial wavefunctions. Compared to the method in Ref. [12], ourmethod does not need the minimization procedure. Asour method uses the full band projection, it can be ap- plied straightforwardly to random potential without us-ing supercells. The effectiveness of our method is demon-strated for both simple bands and composite bands.We consider a simple band that is isolated from otherbands. Composite bands will be discussed later. Ourmethod of constructing MLWFs consists of three steps:1. Guess: choose a set of trial wave functions | g n (cid:105) ,which are localized at lattice sites.2. Projection: | ξ n (cid:105) = P | g n (cid:105) , where P = (cid:80) k | ψ k (cid:105) (cid:104) ψ k | with the summation over the whole first Brillouinzone. Note that P is a full band projection and isdifferent from the projection in Refs. [9, 11], whichis at a given k point. As a result of this crucialdifference, the projected function | ξ n (cid:105) is localizedat site n .3. Orthonormalization: use the L¨owdin orthonormal-ization method [16] to transform | ξ n (cid:105) into a set ofMLWFs | w n (cid:105) . If one uses other methods such asKohn’s method [17] to orthonormalize | ξ n (cid:105) , the re-sulted Wannier function is unlikely maximally lo-calized.Here is why our method is effective and capable of pro-ducing MLWFs. We re-write the full band projection op-erator in terms of Wannier functions, P = (cid:80) n | w n (cid:105) (cid:104) w n | ,where the summation is over all lattice sites. We thushave | ξ n (cid:105) = (cid:88) m | w m (cid:105) (cid:104) w m | g n (cid:105) ≈ (cid:88) m = (cid:104) n (cid:105) | w m (cid:105) (cid:104) w m | g n (cid:105) , (1)where (cid:104) n (cid:105) indicates that the summation is only over site n and its nearest neighbors. It is clear that | ξ n (cid:105) is localizedat lattice site n . However, these projected functions | ξ n (cid:105) ’sare not orthonormal. Any orthonormalization of | ξ n (cid:105) ’swill give us a set of Wannier functions. For example, one a r X i v : . [ c ond - m a t . o t h e r] S e p may use Kohn’s method [17]. We choose the L¨owdinmethod [16] as it can produce the MLWFs. Accord-ing to Ref. [18], the L¨owdin orthogonalization uniquelyminimizes the functional measuring the least squares dis-tance between the given orbitals and the orthogonalizedorbitals. In our case, for the set of projected orbitals | ξ n (cid:105) , the L¨owdin method produces Wannier functions | w n (cid:105) that minimize (cid:88) n (cid:90) dx | (cid:104) x | ξ n (cid:105) − (cid:104) x | w n (cid:105) | , (2)where the summation is over all lattice sites. For a peri-odic potential, this is equivalent to maximizing (cid:104) w n | ξ n (cid:105) = (cid:104) w n | g n (cid:105) . Therefore, if | g n (cid:105) is properly chosen, the re-sulted Wannier function is maximally localized. No fur-ther minimization such as the one in Ref. [12] is needed.To illustrate our method, we consider a single par-ticle Hamiltonian with the periodic potential V ( x ) =cos(2 πx/a ), where a is the lattice constant. We con-struct the Wannier function for its lowest band. As eachpotential well is symmetric with respect to its center (thelowest point of the well), we expect that the MLWF isalso symmetric with its highest peak at the center. So,we choose (cid:104) x | g n (cid:105) = e − ( x − na ) / α / √ α π . It is clearthat (cid:104) w n | g n (cid:105) is the largest for a narrowest Wannier func-tion allowed within the band or MLWF. The resultedWannier function should be insensitive to the width of (cid:104) x | g n (cid:105) . This is exactly what we see in Fig. 1(a). Withthree Gaussians of different widths α = a, . a, . a astrial functions, the resulted Wannier functions fall righton top of each other as seen in Fig. 1(a). We have com-puted the spread of the Wannier functions obtained withour method; they are almost identical for Gaussian trialfunctions of different widths as seen Fig. 1(b).To further illustrate the effectiveness of our method, wehave computed Wannier functions with two other meth-ods. One is the traditional method, where one fixes thephases of the Bloch waves according to the prescriptionin Ref. [10]. The other method is to orthonormalize | ξ n (cid:105) ’s with Kohn’s method [17]. They are compared toour results in Fig. 1(c), where we see that our result isin excellent agreement with the traditional method whilethe one obtained with Kohn’s method is much worse.The L¨owdin method can be implemented differently.However, as long as it is implemented correctly, themethod transforms a given set of non-orthogonal vectorsto a unique set of orthonormal vectors. Nevertheless, weshow here explicitly how we implement it. We imposea periodic boundary condition with N unit cells. As aresult, the crystal wave vector k takes N discrete values k , k , · · · , k N . We let | ψ j (cid:105) = | ψ k j (cid:105) and A nj = (cid:104) ψ j | g n (cid:105) .The L¨owdin orthonormalization is then implemented as | w n (cid:105) = (cid:88) mj ( AA † ) − / nm A mj | ψ j (cid:105) . (3) FIG. 1: (color online) Wannier functions for one dimensionalperiodic potential V ( x ) = cos(2 πx/a ). (a) The Wannier func-tions obtained with our method with trial Gaussian functionsof different widths α . (b) The spread of the Wannier func-tion obtained with our method as a function of the width ofthe Gaussian function α . (c) The Wannier functions obtainedwith three different methods. The black line is our method;the blue line is obtained with the traditional method, whereone fixes the phases of Bloch wave according to the prescrip-tion given in Ref. [10]; the red one is obtained by orthonor-malizing | ξ n (cid:105) ’s with Kohn’s method [17]. If the trial function | g n (cid:105) is translationally symmetric, (cid:104) x | g n (cid:105) = (cid:104) x − r n | g (cid:105) , we have A nj = e − ik j r n A j and( AA † ) nm = (cid:80) j e ik j ( r m − r n ) | (cid:104) ψ j | g (cid:105) | .Our method is applicable for composite bands. Incomposite bands, one or more Wannier functions havenodes. Therefore, to have largest (cid:104) w n | g n (cid:105) for MLWFs,we need to choose such | g n (cid:105) ’s that they have nodesat proper positions. The node positions can be deter-mined by symmetries of the wells. In the worst case,we can determine these node positions by numericallycomputing the eigenstates of the local wells. Here weconsider a two dimensional periodic potential V ( x, y ) = V [cos(2 πx/a ) + cos(2 πy/a )] and use its p-bands to il-lustrate the effectiveness of our method for compos-ite bands. The two trial wave functions are chosen as g , = (cid:112) π α η , e − ( x + y ) / α with η = x, η = y , whichare the two first excited states of a two dimensional har-monic oscillator. The results are plotted in Fig.2. Shownin Fig.2(a) is a comparison between the trial functionand the resulted Wannier function. Three Wannier func-tions obtained from different trial functions are shown inFig.2(c) and (d). Our numerical calculations show thatwhen the width of the trial function α is narrow enough,the resulted Wannier functions are almost identical toeach other. When they are plotted in the figure, theyfall right on top of each other. So, in Fig.2(c) and (d)only three Wannier functions are plotted. The spread of aWannier function (cid:104) w | x + y | w (cid:105) is plotted as the functionof α in Fig.2(b), which shows that the Wannier functionspread no longer changes when α is narrow enough.Our method is directly applicable to random poten-tials. The reason is that we use the full band projection.It can be constructed with all the energy eigenfunctionsin a given band no matter the eigenfunctions are Blochwaves or not. It is well known that even in random po-tentials, there exist eigen-energy “band” that are isolatedfrom other eigen-energies by gaps that are independentof the system size. For these bands, there exist Wan-nier functions [19]. Our method can be used to com-pute these Wannier functions by constructing the pro-jection P with the energy eigenstates of a given randomenergy band. According to Eq.(2), the resulted Wannierfunctions maximize a sum, (cid:80) n (cid:104) w n | g n (cid:105) . With a properchoice of | g n (cid:105) ’s, these Wannier functions can be regardedas maximally localized collectively. There already ex-ist several methods to compute Wannier functions forrandom potentials. Our method has various advantages.The method in Ref. [20] is only applicable to the casewhere the random potential is a perturbation to a peri-odic potential. Kivelson’s method [21] has difficulty fortwo or three dimensional systems. Our method is clearlymore efficient than the one in Ref. [8]. We ourselves haverecently proposed a method to compute Wannier func-tions for random potentials [22]. Our current method iscertainly superior.We now illustrate our method in disordered systems.We choose a disordered potential as a series of cosine-type wells of random depths, V n ( x ) = A n [cos(2 πx/a ) − A n = A [1+ η · R n ], where A is a constant and R n denotes a sequence of random numbers between -0.5and 0.5, and η denotes the relative strength of disorder.For instance, we refer to η = 0 . A = 5 and a = 1, and η = 0 . | g n (cid:105) for different wells despite the wells aredifferent. As shown in Fig. 3, our method produces suc- FIG. 2: (color online) Wannier functions of the composite p -bands in a two dimensional lattice V ( x, y ) = V [cos(2 πx/a ) +cos(2 πy/a )] with V = −
30. (a) The left is the trial wavefunction g with α = 0 . a ; the right is the Wannier functionobtained with this trial function. The white dashed line marksthe unit cell. (b) The spread (cid:104) w | x + y | w (cid:105) of a Wannier func-tion with respect to varying α . The convergence to the blackline is obvious. The black line is the width of the Wannierfunction w ( x, y ) = w ( x ) w ( y ), where w ( x ) as the s -bandWannier function of cos(2 πx/a ) and w ( y ) as the p -band Wan-nier function of cos(2 πy/a ) are obtained with the traditionalmethod. (c) Logarithmic plot of Wannier functions at x = 0along y -axis for different values of α : α = 0 . a (red), 0 . a (green), 0 . a (blue). Note that the latter two almost over-lap. (d) Logarithmic plot of Wannier functions at y along the x -axis for different choices of α : α = 0 . a (red), 0 . a (green),0 . a (blue). y is the highest peak position of the Wannierfunctions. cessfully exponentially localized Wannier functions. Itis clear that a narrower trial Gaussian function leads tomore localized Wannier functions. Numerical results in-dicate that Gaussian trial functions with σ as small as0.5 a is enough to construct MLWFs.For bands of non-trivial topology, the existing meth- FIG. 3: (color online) Wannier functions constructed by ourmethod for one dimensional 30% disordered potential, whichis plotted as the black line at the top of (a). The initialtrial Gaussian functions are the same for different wells. (a)The Wannier functions obtained with Gaussian functions ofdifferent widths α . (b) The spread of the Wannier functionsobtained with our method as a function of the width of thetrial Gaussian function. ods [9, 11] face inherent singularity at special k points.One has to find ingenious ways to circumvent the sin-gularity to obtain exponentially localized Wannier func-tions [23]. Our method has the potential to circumventthe singularity without any special re-design as we usethe full band projection instead of the projection at in-dividual k points.Note that we have so far assumed that the system hastime reversal symmetry. For these systems, its energyeigenfunctions can always be made real. It follows thatthe projection P and all the matrices involved in theL¨owdin method (Eq.3) can also be made real. Therefore,the final Wannier functions | w n (cid:105) are also real as long asthe trial functions are real. For systems where the timereversal symmetry is broken, | w n (cid:105) can be complex. Inthis case, the L¨owdin method maximizes Re {(cid:104) w n | g n (cid:105)} and we may not obtain MLWFs. We shall leave it forfuture discussion.In sum, we have presented a simple and general method for constructing MLWFs. In our method, the full bandprojection is used on localized trial wave functions. Ifthe trial functions are properly chosen to respect the lo-cal potential configuration and have good node positions,the ensuing L¨owdin method orthonormalize them to ML-WFs. No numerical minimization is needed. Our methodcan be directly applied to random potentials. The appli-cation of our method to a real material is left for future.We thank Ji Feng and Xianqing Lin for helpful dis-cussion. This work is supported by the National BasicResearch Program of China (Grants No. 2013CB921903and No. 2012CB921300) and the National Natural Sci-ence Foundation of China (Grants No. 11274024, No.11334001, and No. 11429402). ∗ Electronic address: [email protected][1] G. H. Wannier, Phys. Rev. , 191 (1937).[2] R. Resta, Rev. Mod. Phys. , 899 (1994).[3] R. D. King-Smith and D. Vanderbilt, Phys. Rev. B ,1651 (1993).[4] S. Goedecker, Rev. Mod. Phys. , 1085 (1999).[5] G. Galli, Current Opinion in Solid State and MaterialsScience , 864 (1996).[6] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, andP. Zoller, Phys. Rev. Lett. , 3108 (1998).[7] M. White, M. Pasienski, D. McKay, S. Q. Zhou,D. Ceperley, and B. DeMarco, Phys. Rev. Lett. ,055301 (2009).[8] S. Q. Zhou and D. M. Ceperley, Phys. Rev. A , 013402(2010).[9] N. Marzari and D. Vanderbilt, Phys. Rev. B , 12847(1997).[10] W. Kohn, Phys. Rev. , 809 (1959).[11] H. Teichler, Phys. Status Solidi B , 307 (1971).[12] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, andD. Vanderbilt, Rev. Mod. Phys. , 1419 (2012).[13] H. D. Cornean, I. Herbst, and G. Nenciu, AnnalesHenri Poincar´e (2016) doi:10.1007/s00023-016-0489-2.(arXiv:1506.07435).[14] J. I. Mustafa, S. Coh, M. L. Cohen, and S. G. Louie,Phys. Rev. B , 165134 (2015).[15] E. Canc`es, A. Levitt, G. Panati, and G. Stoltz,arXiv:160507201 (2016).[16] P. L¨owdin, The Journal of Chemical Physics , 365(1950).[17] W. Kohn, Phys. Rev. B , 4388 (1973).[18] J. G. Aiken, J. A. Erdos, and J. A. Goldstein, Interna-tional Journal of Quantum Chemistry , 1101 (1980).[19] A. Nenciu and G. Nenciu, Phys. Rev. B , 10112 (1993).[20] W. Kohn and J. R. Onffroy, Phys. Rev. B , 2485 (1973).[21] S. Kivelson, Phys. Rev. B , 4269 (1982).[22] J. Zhu, Z. Chen, and B. Wu, arXiv:1512.02043 (2015).[23] G. W. Winkler, A. A. Soluyanov, and M. Troyer, Phys.Rev. B93