Robust determination of maximally-localized Wannier functions
We propose an algorithm to determine Maximally Localized Wannier Functions (MLWFs). This algorithm, based on recent theoretical developments, does not require any physical input such as initial guesses for the Wannier functions, unlike popular schemes based on the projection method. We discuss how the projection method can fail on fine grids when the initial guesses are too far from MLWFs. We demonstrate that our algorithm is able to find localized Wannier functions through tests on two-dimensional systems, simplified models of semiconductors, and realistic DFT systems by interfacing with the Wannier90 code. We also test our algorithm on the Haldane and Kane-Mele models to examine how it fails in the presence of topological obstructions.
Wannier functions are a well-established tool in the solid state physics community. In addition to providing intuition on chemical bonding, they are the theoretical and computational underpinning of many developments such as tight-binding approximations , interpolation of band structures , and the so-called modern theory of polarization or orbital magnetization [34, 36].
Wannier functions are however not uniquely defined, and the choice of the phase of the (quasi-)Bloch functions can have a dramatic impact on their spatial localization. As for a single isolated Bloch band, it was realized already in the sixties that in specific situations, appropriate choices of gauge lead to exponentially localized Wannier functions [20, 10, 9]. Progress was made in the eighties by Nenciu  as well as Helffer and Sjöstrand , who proved in the general case the existence of a Bloch gauge yielding exponentially localized Wannier functions.
Whenever Bloch bands intersect each other, as it generically happens in 3D crystals, Bloch functions cannot be smooth – and not even continuous – at the crossing points (excluding the exceptional case of 1D systems). In this situation, a multi-band approach is mandatory. As early realized [1, 10, 9], in insulators and semiconductors it is convenient to consider all the Bloch bands below the Fermi level as a whole, and to replace Bloch functions with quasi-Bloch functions, namely eigenfunctions of the projector on the occupied bands. Their preimages under the Bloch transform are called composite Wannier functions, since they refer to a composite family of energy bands. For the sake of a simpler terminology, we will skip the adjective “composite” in the next sections.
The existence of an orthonormal basis of well-localized Wannier functions is equivalent to the existence of an orthonormal frame of quasi-Bloch functions, spanning , which is both smooth and periodic as a function of . Such an existence problem has been a long-standing problem in theoretical solid-state physics. The solution for 1D systems was provided in [28, 30], while a solution for 2D and 3D systems required the crucial use of geometric ideas and methods [2, 31].
From a practical viewpoint, Marzari and Vanderbilt introduced an optimization procedure to minimize the spread of Wannier functions  and compute Maximally-Localized Wannier Functions (MLWFs). The corresponding localization functional is now known as the Marzari-Vanderbilt (MV) functional. The approach based on the minimization of the MV functional was developed before a full understanding of the theoretical criteria for localization was gained. It however gives very satisfactory results in many situations, and has become one of the standard tools of computational solid state physics. It was conjectured in  that global minimizers of the MV functional are exponentially localized, a fact which was proved recently in , thereby providing a firm and mathematically consistent ground for the MV optimization.
It has however been observed that, in some situations, the minimization of the MV functional could fail because the algorithm remained trapped in “false local minima” presenting unphysical oscillations [22, 25]. This issue is to be distinguished from physically relevant “real” local minima  and was associated with bad initial guesses obtained by the projection method and fine -point meshes. This problem has been considered recently in , where an alternative, more robust method is presented. This method is however still based on a projection and requires a physical input in the choice of basis functions, while the algorithm we present here does not require any input or parameter tuning. In that sense it is similar to the recently proposed SCDM algorithm [7, 8], which uses columns of the density matrix to provide localized Wannier functions. Compared to this method, our approach works directly in Bloch space, avoiding any representation of the density matrix in real space, and is readily implementable using the same input as standard MLWFs computation. It also does not depend on a potentially ill-conditioned matrix inversion.
We will show in this work that the issue of “false local minima” occurs when the initial guess corresponds to a Bloch gauge with vortex-like discontinuities, which prevent the convergence of the MV optimization algorithm on fine samplings of the Brillouin zone. To provide an admissible initial guess for the MV algorithm, we should really find a continuous Bloch gauge, a mathematically non-trivial task. This issue has been studied recently in  and , where the authors develop constructive algorithms which complement the abstract existence results of [2, 31]. Both methods are however not trivial to implement in practice, that of  requiring a perturbation argument that might yield to Bloch gauge with very sharp transitions, and that of  relying on cumbersome interpolation procedures.
In this work, we present an algorithm inspired by the theoretical works of  and  but more suited for actual implementation. This method fixes the gauge at the point in the Brillouin zone and then incrementally on its neighboring points, until the whole Brillouin zone is covered. Certain pseudo-periodicity conditions have then to be satisfied, which we enforce by the method of obstruction matrices developed in . This requires finding continuous and periodic logarithms of families of unitaries, a problem discussed in . We propose an algorithm to solve this problem. Although we are unable to prove that this step always produces a continuous Bloch gauge, we conjecture it to be the case, and support this claim by numerical tests.
Our algorithms only provides a continuous Bloch gauge, which yields algebraically-decaying Wannier functions. Therefore, we see this algorithm not as an end in itself but as a way to obtain good initial guesses for the MV localization procedure. We show by numerical examples on two- and three-dimensional systems that this two-step process always yields localized Wannier functions, which is not the case for state-of-the-art methods based on the projection method. We therefore advocate the use of our method when good initial guesses are not available. We however emphasize that our method only finds “real” local minima of the Marzari-Vanderbilt functional, which might or might not possess optimal spread (i. e. be a global minimum). Since they are devoid of any physical input, they might also not correspond to physically intuitive Wannier functions. Therefore, we do not see our method as a replacement of the traditional construction of Wannier functions for systems where physical insight are available, but rather as an automated procedure that will always produce exponentially-localized Wannier functions in the case when intuition is lacking.
This article is organized as follows. Section 2 introduces the notation and recalls some well-known facts on Wannier functions. We next present our algorithm in a one-dimensional case in Section 3, then its extensions to the two- and three-dimensional settings in Sections 4 and 5. The subtle issue of the appropriate choice of continuous and periodic logarithms is discussed in Section 6. We finally illustrate the approach with numerical results in Section 7, and present some perspectives in Section 8.
2 Notation and statement of the problem
We consider a -dimensional crystal. We denote by the corresponding Bravais lattice, by a unit cell centered on the origin, and by the reciprocal lattice, with . We choose the vectors as a basis for the reciprocal lattice. In these coordinates, the Brillouin zone is
where the symbol refers to the fact that hereafter we will identify the vector with its components with respect to the reciprocal basis.
2.1 Smooth families of projectors and their symmetries
The key input of our algorithm is a family of orthogonal projectors of rank , which is smooth with respect to the wave vector . Such projectors are obtained, in standard solid-state physics computations, by spectral projections of the effective one-body Schrödinger operator
where is a real-valued -periodic potential. More precisely, introduce the Bloch orbitals , with periodic parts , so that
The integer labels the band index, while denotes the quasi-momentum. For a fixed , the periodic functions form an orthogonal basis of consisting of solutions of the following eigenvalue problem:
where the are labeled by in increasing order. Our convention for the Bloch transform is
In order to construct Wannier functions from the Bloch orbitals, we identify a set of Bloch bands in the energy spectrum. We assume that they are isolated in the sense that
This gap condition is satisfied for instance by the occupied bands of an insulator or a semiconductor. It ensures that the spectral projector
The underlying symmetries of the Hamiltonian translate into corresponding symmetries of the family of projectors . We first note that, because the potential is real, satisfies the time-reversal property
where is the anti-unitary operator corresponding to complex conjugation, namely
We also note that, for any ,
with the following unitary translation operators
The convention we use here is slightly different from the one in : our definition of ensures that is an eigenvector of (so that is a translation by ). The family of operators is a unitary group representation of , in the sense that
and map functions defined for to functions defined for any by reflection and translation.
Similarly to , we do not use the specific structure of the model (1), so that the approach presented here can be applied to various other periodic models of quantum mechanics, such as tight-binding models or relativistic models described by a Dirac operator. We only require a continuous family of projectors satisfying (3) and (4), for some unitary operators and an anti-unitary operator satisfying
It is possible to further restrict the Brillouin zone under consideration when has more symmetries than the time-reversal and translation symmetries (3) and (4). The construction we present here does not take into account these possible additional symmetries, and therefore may produce Wannier functions that do not respect these symmetries. Our algorithm could be extended to work only on the irreducible Brillouin zone, but this is outside the scope of this paper.
2.2 Orthonormal frames
The fundamental element to construct well localized Wannier functions and the output of our algorithm is a Bloch frame depending smoothly on the wave vector . By definition, a Bloch frame (or, shortly, a frame) is a mapping from to an orthonormal basis of , such that each component satisfies the pseudo-periodicity property
The latter condition expresses the fact that such functions are compatible with the symmetries of the family of projectors given by (4).
We emphasize that, despite the similarity in the notation, we do not assume that the functions are (the periodic part of) Bloch functions. It is only assumed that
while does not hold true in general. Such functions are sometimes called quasi-Bloch functions in the literature.
As already noticed , the existence of a smooth Bloch frame is not trivial in view of the competition between the smoothness of and the property (6), which encodes a global topological constraint. For instance, it is known that in some models with broken time-reversal symmetry (e. g. in the presence of a magnetic field or in a Chern insulator) there cannot exist any such continuous frame, due to a topological obstruction [11, 16, 2, 24].
For two given frames and , we write for the matrix with entries
For a matrix , we write
Note that, if is a frame, then so is when the matrices are unitary for all values of . Similarly, for an operator , is obtained by applying to all components of independently. These conventions allow for a uniform notation whether are functions or coefficients in a basis set of size (in which case is a matrix, and the previous definitions are simply restatements of matrix multiplication rules).
2.3 Well localized Wannier functions
(Composite) Wannier functions are defined, for a given frame , by
and the translations of these functions:
The Wannier functions form a complete orthonormal basis of the subspace of associated with the chosen bands.
The localization properties of the Wannier functions are determined by the regularity of the frame . This can be seen from the fact that, because satisfies the pseudo-periodicity property (6), an integration by parts gives
Therefore, a frame yields Wannier functions that decay faster than any polynomial, and an analytic frame yields exponentially localized Wannier functions (see Section 2 of  for a more precise statement involving the Sobolev regularity of ). In practice, the frame straightforwardly coming out of the simulations, computed on the reduced Brillouin zone, is usually not smooth because of the arbitrary phases of the and of the possible crossings between the energy levels . Therefore, some correction through a family of unitary matrices has to be applied to the frame in order to transform it into a smooth frame . Although it is easy to fix phases to impose local smoothness, global smoothness is a more complicated issue because the frame, extended by reflection and translation to the full space , might not be continuous at the boundary of the reduced Brillouin zone. Certain compatibility or gluing conditions should be satisfied in order for this extension to be smooth.
Our problem is to compute real-valued localized Wannier functions, which, after Bloch transform, is equivalent to the following.
Find a smooth frame such that
The additional property (8) ensures that the Wannier functions are real-valued.
We discuss in the next section an algorithm for constructing continuous frames. We do not attempt to impose a higher degree of regularity, as done in [12, 4] via abstract procedures, but consider a more pragmatic approach where a subsequent smoothing of the continuous frame is obtained with the MV algorithm. Our hope is that the continuous frame we obtain at the end of our procedure is a sufficiently good initial guess for the minimization of the MV functional in order to actually converge towards a minimizer.
3 Algorithm in the one-dimensional case
We present here a practical algorithm to solve Problem 3 in dimension 1. It can be seen as a discretization of well-known procedures to construct smooth frames , dating back to at least , and discussed in  section IV.C.1. We however carefully present this simple case since, as in , it is the building block of the algorithm in higher dimensions.
We start from a given frame on the reduced Brillouin zone , which we then extend to the whole space by the relations (6) and (8). More precisely, we first extend to by , and then to any by where is such that . This procedure yields a globally continuous frame if and only if the frame on is continuous and satisfies the following compatibility conditions at and :
Step 1: Fix .
The first step is to choose a frame satisfying (P1). Since is real-valued, this can be done by choosing an orthonormal set of real-valued eigenfunctions of .
Step 2: Propagate from to .
Evolving the eigenvector at to can be done in various ways, for instance by the Sz.-Nagy intertwining unitary as done in . Here we describe a natural way of performing this operation on a mesh of -points, using a Löwdin orthogonalization procedure. Assume that a mesh of is given. We construct iteratively
Since is continuous, provided the mesh spacing is sufficiently small, the overlap matrix has its eigenvalues bounded away from 0, so that is a well-defined orthonormal basis of .
When the mesh spacing tends to zero, it can be shown that converges to the solution of the following ODE, with initial condition :
In particular, the function is smooth.
Step 3: Enforcing (P2) at .
The previous step yields a frame that is continuous, satisfies the compatibility condition (P1), but not (P2). Note however that , so that and are both orthonormal bases of the same space. There exists therefore a unitary matrix , which we call “obstruction matrix” following , such that
This matrix can be explicitly computed as . A simple computation also shows that :
We see that satisfies the compatibility condition (P2) if and only if , i. e. . Since , we can choose
With this choice, satisfies (P2). We define the fractional power of a normal matrix in the usual way, by applying the exponent to its eigenvalues. For matrices with complex eigenvalues such as , defining requires choosing a branch cut in the complex plane for the logarithm, which we fix by using the principal determination, i. e. for all .
Step 4: Correct the frame on .
The last step is to globally define a new frame on by interpolation:
Note that , and so satisfies (P1). In the limit of infinitely small mesh spacing, the mapping is continuous, and so is , which implies that is continuous on .
The compatibility conditions finally ensure that can be extended to a continuous frame for all .
We have only imposed continuity in our construction, but we actually get a smooth frame. This is because the frame we constructed is the same as the one we would obtain by defining on via the propagation (9), then . The latter frame is smooth because both and are.
4 Algorithm in the two-dimensional case
As in the one-dimensional case, we first construct a frame on the reduced Brillouin zone
and then extend it by symmetry to the whole space using (6) and (8). This yields a globally continuous frame if and only if the frame on the reduced Brillouin zone is continuous and the following compatibility conditions are satisfied:
Note that these conditions are now edge conditions, relating two points of the same edge (for (E1) and (E3)) or two points on opposite edges (for (E2)), in contrast to the point conditions of the one-dimensional case. These edge conditions imply compatibility conditions at the special points and , which are the fixed points of the symmetry group. All these conditions can be visualized on Figure 1.
We first use the one-dimensional construction to obtain a frame with , satisfying (E1) as well as
Then, for every , we propagate from to using the propagation procedure described in (9) in the one-dimensional case. This defines a frame on .
Step 2: Enforcing the compatibility conditions (E3) at the corners .
This is a point compatibility condition analogous to the one encountered in the one-dimensional case. To enforce it, we consider the unitary obstruction matrix
and define the modified frame
The modified frame satisfies (13) by construction. Note also that the transformation (14) preserves the compatibility conditions (E2) and (E1). Moreover, (E2) and (13) imply that also satisfies the corner condition at :
In order to simplify the notation, we drop the prime in the remainder of the algorithm, and simply denote by the modified frame obtained at the end of this step.
Step 3: Enforcing (E3) on the right edge.
We still need to modify to satisfy (E3) for any value of . We define to this end the following family of obstruction matrices:
These matrices satisfy by construction
At this point, we would like to define a new frame by
The frame satisfies all the compatibility conditions (E1)-(E2)-(E3), is continuous with respect to , but may fail to be continuous with respect to , because the eigenvalues of may pass through the branch cut of the logarithm in the negative real axis. We are therefore faced with the problem of the continuous periodic logarithm: find a matrix valued function , with hermitian, which is continuous and satisfies the following conditions for all :
If such a function exists, the frame
The behavior of the algorithm on a simple example can be visualized in Figure 2.
The algorithm produces a continuous frame. However, in contrast to the 1D case, it does not provide any additional regularity. Because the 1D construction guarantees a smooth frame, there are no discontinuities of derivatives in the vertical direction, i. e. at the points . However, the condition for the first derivative to be continuous at a point is
which is in general incompatible with (E3). This yields discontinuous derivatives at .
This means that the Wannier function will only decay algebraically in the direction, as can be seen in the right panel of Figure 2. This is however easily fixed by the application of the MV minimization algorithm, as will be shown on Figure 5 in Section 7. We note that the MV algorithm is essentially equivalent to other forms of smoothing, such as for instance the one employed in : the MV algorithm can be seen as the gradient flow of a Dirichlet-like energy functional, i. e. a heat flow, whose solution can be written as a convolution similar to the one used in .
After step 1, the frame satisfies the compatibility conditions related to vertical translations, and the associated Wannier function is localized vertically, but not horizontally (left panel). After step 2, the frame satisfies all the compatibility conditions at the corners (with ), but not on the vertical edges between these points. The associated Wannier function is still delocalized horizontally (middle panel). After step 3, the frame satisfies all the compatibility conditions and is globally continuous, although not smooth. The associated Wannier function is localized, but exhibits a slow algebraic horizontal decay (right panel). See Section 7 for the details of the model used.
5 Algorithm in the three-dimensional case
We use the same induction to pass from the two-dimensional to the three-dimensional case as the one to pass from the one-dimensional to the two-dimensional case. We define a frame on the reduced Brillouin zone
and extend it to the full reciprocal space by symmetry. The compatibility conditions to be satisfied by the frame are now
Step 2: Enforcing (F4) on the edges of the face .
As in the two-dimensional case, before enforcing (F4) on the whole face , we enforce it on the four edges of its boundary. In order to do this, we first fix the corner , for which the compatibility condition is
We define the obstruction matrix
By construction, satisfies (15). The compatibility conditions (F1)-(F2)-(F3) are also still valid for the modified frame . In addition, in view of (15) and (F2)-(F3), the frame satisfies the compatibility conditions at the other corners and .
We next enforce the condition (F4) on the edge , namely:
The corresponding unitary obstruction matrix is
Note that and . As in the two-dimensional case, provided we can solve the logarithm problem, we find a hermitian logarithm satisfying , and , and modify the frame as
The modified frame then satisfies (F4) on the two edges .
We repeat this construction on the edge , and introduce appropriate hermitian matrices such that
satisfies (F4) on the two edges .
To simplify the notation, we drop the primes in the remainder of the algorithm, and simply denote by the modified frame obtained at the end of this step.
Step 3: Enforcing (F4) on the whole face .
The compatibility condition (F4) is, on the whole face :
We therefore introduce the family of unitary obstruction matrices
Because of the previous step,
As in the two-dimensional case, we would like to define a new frame satisfying all the compatibility conditions as
We again face with the logarithm problem, which we discuss in Section 6.
As a mathematical curiosity, this construction immediately extends to higher dimensions, again provided the logarithm problem can be solved. By induction, we first build a frame on , then propagate it to . We next fix the obstruction on the boundary of , and conclude by finding a continuous logarithm of the obstruction matrix on to correct the frame.
6 The logarithm problem
We come back in this section to the so-called logarithm problem, summarized in Problem 8 below.
An element is denoted in this section. We also introduce the set of edges in dimension 3 (vertices in dimension 2)
With this notation, we can summarize the logarithm problem as follows.
Given a continuous mapping from to the set of unitary matrices, such that
find a continuous mapping from to the set of hermitian matrices, such that
and, for all ,
This is the problem as used by our construction in the previous two sections: (16) is satisfied because we fixed the obstruction matrix on the edges, while (17) stems from the structure of the conditions (E3) and (F4). Let us already note that this problem cannot be solved in general (see the counter-example in , which we recall in Section 6.3). We describe an algorithm that solves the problem provided that the type of eigenvalue collisions present in this counter-example does not occur. In our numerical tests, we did not encounter such difficulties and were always able to solve the logarithm problem. We refer to Section 6.3 for a discussion of this issue.
6.1 The two-dimensional case
A logarithm of can be found by starting at , and by following the eigenvalues of to ensure the continuity of . This procedure can be seen as a discrete counterpart of the phase following procedure used in [4, Lemma 2.13] for a single band.
Assume that is partitioned as . We set , and iteratively determine as a logarithm of , taking into account the phase information in . To this end, we diagonalize as , where is unitary and a diagonal matrix whose entries have modulus 1. We then set , where is the diagonal matrix given by