Construction of maximally localized Wannier functions

Abstract

We present a general method for constructing maximally localized Wannier functions. It consists of three steps: (i) picking a localized trial wave function, (ii) performing a full band projection, and (iii) orthonormalizing with the Löwdin method. Our method is capable of producing maximally localized Wannier functions without further minimization, and it can be applied straightforwardly to random potentials without using supercells. The effectiveness of our method is demonstrated for both simple bands and composite bands.

Graphical abstract

Keywords

Wannier function / random potential / cold atomic gases

Wannier functions that are localized at lattice sites are a representation of electronic states in crystalline solids, alternative to Bloch waves [1]. They offer an insightful picture of chemical bonds, play a pivotal role in the modern theory of polarization [2, 3], and are the basis for efficient linear-scaling algorithms in electric-structure calculations [4, 5]. Wannier functions are also important in linking cold atom experiments in continuous light potentials with lattice Hamiltonians, such as the Bose-Hubbard model and Anderson random lattice [6-8].

Although there are explicit formulae that transform Bloch waves to Wannier functions, the construction of Wannier functions is far from trivial. The primary reason is that there are infinite Wannier functions for a given band because of the phase choices for Bloch waves, whereas in practice, the maximally localized Wannier functions (MLWFs) are sought and preferred [9]. Kohn showed a method to fix the Bloch wave phases for one-dimensional lattice to obtain these MLWFs [10]. In 1971, Teichler found a general method to construct Wannier functions, which is insensitive to the phases of Bloch waves [11]. However, this method does not guarantee maximal localization and depends on the initial trial function. A more sophisticated method involving numerical minimization of the spread of a Wannier function was developed to compute MLWFs [12]. There were also some new developments recently [13-15].

In this work, we present a general method for computing MLWFs. Our method is somewhat similar to Teichler’s method as it also involves projection and the Löwdin orthonormalization method [16]. Nevertheless there is a significant difference so that our method can produce MLWFs and is insensitive to the initial trial wave functions. Compared to the method in Ref. [12], our method does not need the minimization procedure. It uses the full band projection; therefore, it can be applied straightforwardly to random potentials without using supercells. The effectiveness of our method is demonstrated for both simple and composite bands.

We consider a simple band that is isolated from the other bands. Composite bands will be discussed later. Our method of constructing MLWFs consists of three steps:

1) Guess: Choose a set of trial wave functions |g_n\rangle, which are localized at lattice sites.

2) Projection: |{\xi }_n\rangle =P|g_n\rangle, where P={\displaystyle {\sum }_k}|{\psi }_k\rangle \langle {\psi }_k| with the summation over the whole first Brillouin zone for a given band. Note that P is a full band projection and different from the projections in Refs. [9, 11], which are at a given k point. As a result of this crucial difference, the projected function |{\xi }_n\rangle is localized at site n.

3) Orthonormalization: Use the Löwdin orthonormalization method [16] to transform |{\xi }_n\rangle into a set of MLWFs |w_n\rangle. If one uses other methods such as Kohn’s method [17] to orthonormalize |{\xi }_n\rangle, the resulting Wannier function is unlikely to be maximally localized.

Here is why our method is effective and capable of producing MLWFs. We re-write the full band projection operator in terms of Wannier functions, P={\displaystyle {\sum }_n}|w_n\rangle \langle w_n|, where the summation is over all the lattice sites. We thus have

\begin{array}{}\text{(1)}& |{\xi }_n\rangle ={\displaystyle \sum _m}|w_m\rangle \langle w_m|g_n\rangle \approx {\displaystyle \sum _{m=\langle n\rangle }}|w_m\rangle \langle w_m|g_n\rangle \mathrm{,}\end{array}

where \langle n\rangle indicates that the summation is only over site n and its nearest neighbors. It is clear that |{\xi }_n\rangle is localized at lattice site n. However, these projected functions |{\xi }_n\rangle’s are not orthonormal. Any orthonormalization of |{\xi }_n\rangle’s will give us a set of Wannier functions. For example, one may use Kohn’s method [17]. We choose the Löwdin method [16] as it can produce MLWFs. According to Ref. [18], the Löwdin orthogonalization uniquely minimizes the functional measuring the least squares distance between the given orbitals and the orthogonalized orbitals. In our case, for the set of projected orbitals |{\xi }_n\rangle, the Löwdin method produces Wannier functions |w_n\rangle that minimize

\begin{array}{}\text{(2)}& {\displaystyle \sum _n}{\displaystyle \int }\text{d}x|\langle x|{\xi }_n\rangle -\langle x|w_n\rangle {|}^2\mathrm{,}\end{array}

where the summation is over all the lattice sites. For a periodic potential, this is equivalent to maximizing \langle w_n|{\xi }_n\rangle =\langle w_n|g_n\rangle. Therefore, if |g_n\rangle is properly chosen, the resulting Wannier function is maximally localized. No further minimization such as the one in Ref. [12] is needed.

We consider a single particle Hamiltonian with the periodic potential V(x)=\cos (2\pi x/a), where a is the lattice constant, in order to illustrate our method. We construct the Wannier function for its lowest band. As each potential well is symmetric with respect to its center (the lowest point of the well), we expect that the MLWF is also symmetric with its highest peak at the center. Therefore, we choose \langle x|g_n\rangle ={\text{e}}^{-{(x-na)}^2/(2{\alpha }^2)}/\sqrt{2{\alpha }^2\pi }. It is clear that \langle w_n|g_n\rangle is the largest for the narrowest Wannier function allowed within the band or MLWF. The resulting Wannier function should be insensitive to the width of \langle x|g_n\rangle. This is exactly what we see in Fig. 1(a). With three Gaussians of different widths \alpha =a\mathrm{,0.5}a\mathrm{,0.1}a as trial functions, the resulting Wannier functions fall right on top of each other as seen in Fig. 1(a). We have computed the spread of the Wannier functions obtained with our method; they are almost identical for Gaussian trial functions of different widths as seen in Fig. 1(b).

We computed Wannier functions with two other methods to illustrate the effectiveness of our method further. One of them is the traditional method, where one fixes the phases of the Bloch waves according to the prescription in Ref. [10]. The other method is to orthonormalize |{\xi }_n\rangle’s with Kohn’s method [17]. The results obtained using these different methods are compared in Fig. 1(c), which shows that our result is in excellent agreement with that of the traditional method, whereas the results obtained with Kohn’s method are much worse.

The Löwdin method can be implemented in different ways. However, as long as it is implemented correctly, the method transforms a given set of non-orthogonal vectors to a unique set of orthonormal vectors. Nevertheless, we explicitly demonstrate its implementation here. We impose a periodic boundary condition with N unit cells. As a result, the crystal wave vector k takes N discrete values k_1\mathrm{,}{k}_2\mathrm{,}\cdots \mathrm{,}{k}_N. We let |{\psi }_j\rangle =|{\psi }_{k_j}\rangle and A_{nj}=\langle {\psi }_j|g_n\rangle. The Löwdin orthonormalization is then implemented as

\begin{array}{}\text{(3)}& |w_n\rangle ={\displaystyle \sum _{mj}}{(A{A}^{†})}_{nm}^{-1/2}{A}_{mj}|{\psi }_j\rangle \mathrm{.}\end{array}

If the trial function |g_n\rangle is translationally symmetric, \langle x|g_n\rangle =\langle x-r_n|g_0\rangle, we have A_{nj}={\text{e}}^{-\text{i}{k}_j{r}_n}{A}_{0j} and {(A{A}^{†})}_{nm}={\displaystyle {\sum }_j}{\text{e}}^{\text{i}{k}_j(r_m-r_n)}|\langle {\psi }_j|g_0\rangle {|}^2.

Our method is applicable to composite bands. In composite bands, one or more Wannier functions have nodes. Therefore, to achieve the largest \langle w_n|g_n\rangle for MLWFs, we need to choose |g_n\rangle’s such that they have nodes at proper positions. The node positions can be determined from the symmetries of the wells. In the worst case, we can determine these node positions by numerically computing the eigenstates of the local wells. Here we consider a two dimensional periodic potential V(x\mathrm{,}y)=V_0[\cos (2\pi x/a)+\cos (2\pi y/a)] and use its p-bands to illustrate the effectiveness of our method for composite bands. The two trial wave functions are chosen as g_{\mathrm{1,2}}=\sqrt{\frac{\pi }{2{\alpha }^4}}{\eta }_{\mathrm{1,2}}{\text{e}}^{-(x^2+y^2)/2{\alpha }^2} with {\eta }_1=x\mathrm{,}{\eta }_2=y, which are the first two excited states of a two dimensional harmonic oscillator. The results are plotted in Fig. 2. Fig. 2(a) shows a comparison between the trial function and the resulting Wannier function. Three Wannier functions obtained from different trial functions are shown in Figs. 2(c) and (d). Our calculations show that when the width of the trial function α is sufficiently narrow, the resulting Wannier functions are almost identical to each other. When they are plotted in the figure, they fall right on top of each other. Therefore, in Figs. 2(c) and (d) only three Wannier functions are plotted. The spread of a Wannier function \langle w|x^2+y^2|w\rangle is plotted as a function of α in Fig. 2(b), which shows that the Wannier function spread does not change when α is sufficiently narrow.

Our method is directly applicable to random potentials. The reason is that we use the full band projection. It can be constructed with all the energy eigenfunctions in a given band regardless of whether the eigenfunctions are Bloch waves. It is well known that even in random potentials, there exist eigen-energy “bands” that are isolated from other eigen-energies by gaps that are independent of the system size. Wannier functions exist for these bands [19]. Our method can be used to compute these Wannier functions by constructing the projection P using the energy eigenstates of a given random energy band. According to Eq. (2), the resulting Wannier functions maximize a sum, {\displaystyle {\sum }_n}\langle w_n|g_n\rangle. With a proper choice of |g_n\rangle’s, these Wannier functions can be regarded as maximally localized collectively. There are several methods to compute Wannier functions for random potentials. Our method has various advantages. The method in Ref. [20] is only applicable to cases in which the random potential is a perturbation to a periodic potential. Kivelson’s method [21] has difficulty for two or three dimensional systems. Our method is clearly more efficient than the one in Ref. [8]. We recently proposed a method to compute Wannier functions for random potentials [22]. Our current method is certainly 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\pi x/a)-1]. The well depth A_n=A[1+\eta \cdot {\Re }_n], where A is a constant, {\Re }_n denotes a sequence of random numbers between -0.5 and 0.5, and η denotes the relative strength of disorder. For instance, we refer to \eta =0.1 as a 10% disorder. We use A=5, a=1, and \eta =0.3 in the example shown in Fig. 3. In our computation, we choose the same Gaussian trial functions |g_n\rangle for different wells despite them being different. As shown in Fig. 3, our method successfully produces exponentially localized Wannier functions. It is clear that a narrower trial Gaussian function leads to more localized Wannier functions. Numerical results indicate that Gaussian trial functions with σ as small as 0.5a are sufficient to construct MLWFs.

Note that we have so far assumed that the system has time reversal symmetry. The energy eigenfunctions of these systems can always be made real. It follows that the projection P and all the matrices involved in the Löwdin method [Eq. (3)] can also be made real. Therefore, the final Wannier functions |w_n\rangle are also real as long as the trial functions are real. |w_n\rangle can be complex for systems in which the time reversal symmetry is broken. In this case, the Löwdin method maximizes \text{R}\text{e}\{\langle w_n|g_n\rangle \} and we may not obtain MLWFs. We shall leave it for future discussion.

In sum, we have presented a simple and general method for constructing MLWFs. In our method, the full band projection is used on localized trial wave functions. If the trial functions are properly chosen to respect the local potential configuration and have good node positions, the ensuing Löwdin method orthonormalizes them to MLWFs. No numerical minimization is needed. Our method can be directly applied to random potentials. The application of our method to real material is left for future.

References

[1] G. H. Wannier, The structure of electronic excitation levels in insulating crystals, Phys. Rev. 52, 191 (1937)

[2] R. Resta, Macroscopic polarization in crystalline dielectrics: The geometric phase approach, Rev. Mod. Phys. 66, 899 (1994)

[3] R. D. King-Smith and D. Vanderbilt, Theory of polarization of crystalline solids, Phys. Rev. B 47, 1651 (1993)

[4] S. Goedecker, Linear scaling electronic structure methods, Rev. Mod. Phys. 71, 1085 (1999)

[5] G. Galli, Linear scaling methods for electronic structure calculations and quantum molecular dynamics simulations, Current Opinion in Solid State and Materials Science 1(6), 864 (1996)

[6] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Cold bosonic atoms in optical lattices, Phys. Rev. Lett. 81, 3108 (1998)

[7] M. White, M. Pasienski, D. McKay, S. Q. Zhou, D. Ceperley, and B. DeMarco, Strongly interacting bosons in a disordered optical lattice, Phys. Rev. Lett. 102, 055301 (2009)

[8] S. Q. Zhou and D. M. Ceperley, Construction of localized wave functions for a disordered optical lattice and analysis of the resulting Hubbard model parameters, Phys. Rev. A 81, 013402 (2010)

[9] N. Marzari and D. Vanderbilt, Maximally localized generalized Wannier functions for composite energy bands, Phys. Rev. B 56, 12847 (1997)

[10] W. Kohn, Analytic properties of Bloch waves and Wannier functions, Phys. Rev. 115, 809 (1959)

[11] H. Teichler, Best Localized Symmetry-Adapted Wannier Functions of the Diamond Structure, Phys. Status Solidi B 43, 307 (1971)

[12] 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)

[13] H. D. Cornean, I. Herbst, and G. Nenciu, On the construction of composite Wannier functions, arXiv: 1506.07435 (2015)

[14] J. I. Mustafa, S. Coh, M. L. Cohen, and S. G. Louie, Automated construction of maximally localized Wannier functions: Optimized projection functions method, Phys. Rev. B 92, 165134 (2015), arXiv: 1508.04148 (2015)

[15] E. Cancès, A. Levitt, G. Panati, and G. Stoltz, Robust determination of maximally-localized Wannier functions, arXiv: 1605.07201 (2016)

[16] 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)

[17] W. Kohn, Construction of Wannier functions and applications to energy bands, Phys. Rev. B 7, 4388 (1973)

[18] J. G. Aiken, J. A. Erdos, and J. A. Goldstein, You have full text access to this content On Löwdin orthogonalization, Int. J. Quantum Chem. 18, 1101 (1980)

[19] A. Nenciu and G. Nenciu, Existence of exponentially localized Wannier functions for nonperiodic systems, Phys. Rev. B 47, 10112 (1993)

[20] W. Kohn and J. R. Onffroy, Wannier functions in a simple nonperiodic system, Phys. Rev. B 8, 2485 (1973)

[21] S. Kivelson, Wannier functions in one-dimensional disordered systems: Application to fractionally charged solitons, Phys. Rev. B 26, 4269 (1982)

[22] J. Zhu, Z. Chen, and B. Wu, Construction of Wannier functions in disordered systems, arXiv: 1512.02043 (2015)

Rights & permissions

Higher Education Press and Springer-Verlag Berlin Heidelberg 2017