# Symmetry-adapted Wannier functions in the maximal localization procedure

###### Abstract

A procedure to construct symmetry-adapted Wannier functions in the framework of the maximally-localized Wannier function approach[Marzari and Vanderbilt, Phys. Rev. B 56, 12847 (1997); Souza, Marzari, and Vanderbilt, ibid. 65, 035109 (2001)] is presented. In this scheme the minimization of the spread functional of the Wannier functions is performed with constraints that are derived from symmetry properties of the specified set of the Wannier functions and the Bloch functions used to construct them, therefore one can obtain a solution that does not necessarily yield the global minimum of the spread functional. As a test of this approach, results of atom-centered Wannier functions for GaAs and Cu are presented.

###### pacs:

71.15.Ap## I introduction

After its proposal by Marzari and VanderbiltMarzari97 (), the maximally-localized Wannier function approachMarzari97 (); Souza01 (); Marzari12 () has been widely used as a convenient tool to construct localized orthonormal functions in crystals. These Wannier functions are obtained from unitary transformation of the Bloch functions, whose phase factors are chosen so that the spatial spreads of the Wannier functionsBlount () are minimized. The maximally-localized Wannier functions have been employed for a number of applications, such as analysis of chemical bonding and basis functions for linear-scaling calculations or model Hamiltonians in strongly-correlated systemsMarzari12 ().

As found by Souza et al.Souza01 () in the case of the -like Wannier function of Cu, the Wannier functions obtained in the maximal-localization procedure do not necessarily reflect the spatial symmetry of the system; the center of the 4-like Wannier function of Cu is not at the Cu atom but at the tetrahedral interstitial site. A similar result was reported by Thygesen et al.Thygesen05 (). This feature of the maximal-localization approach has some drawbacks; the center of the maximally-localized Wannier functions is not necessarily on atom positions or other high-symmetry points, which sometimes makes the interpretation of the obtained Wannier functions difficult. Furthermore, due to the lack of definite symmetry in the Wannier functions, one has to calculate the transformation matrix from the Bloch functions to the Wannier functions for all -points in the first Brillouin zone, not only inside the irreducible part of it.

The connection between the symmetry of the crystal and the properties of Wannier functions was first discussed by des CloizeauxCloizeaux63 () from the view point of group theory. His basic idea is that the Wannier functions can be chosen to be the basis of the irreducible representations of a subgroup of the full symmetry group of the system whose elements do not change the given point in the unit cell, and he derived the relation between these Wannier functions and the eigenfunctions of the one-particle Hamiltonian of the crystal(i.e. Bloch states). This idea is based on the site-symmetry group and the theory of the induced representationsEvarestov97 (), and there have been a number of works considering these symmetry-adapted Wannier functionsKohn73 (); Boehm79 (); Kruger87 (); Sporkmann94 (); Evarestov97 (); Smirnov01 (); Smirnov05 (); Posternak02 (); Casassa06 (). Based on this idea, in this work we propose a procedure to construct symmetry-adapted Wannier functions in the framework of the maximally-localized Wannier function approach. Considering the symmetry properties of the specified set of the Wannier functions and the Bloch functions used, we derive a formula that the transformation matrix follows, and perform the minimization of the spread functional with this symmetry constraint. Our procedure enables one to control the symmetry and center of the Wannier functions, and it also enables one to generate the transformation matrix for a general -point from its symmetry-equivalent point inside the irreducible Brillouin zone(IBZ), which simplifies the minimization of the spread functional. As a test of our approach, we consider Wannier functions of GaAs and Cu.

## Ii method

### ii.1 Symmetry-adapted Wannier functions

#### ii.1.1 Site symmetry group and symmetry-adapted Wannier functions

In this section we summarize the main points of site symmetry group and symmetry-adapted Wannier functions following Ref. Evarestov97, . We refer to Refs. Cloizeaux63, ; Evarestov97, for details of the theory.

The starting point of this idea is to specify a set of positions in real space(“sites”) in which one or more Wannier functions will be centered. These sites can be either atomic positions, at chemical bonds, or interstitial sites, depending on the case of interest. The site symmetry group of a given point , denoted by , is a subgroup of the full symmetry group of the crystal whose elements leave unchanged: namely, satisfies

(1) |

where are the rotation and the translation part of the symmetry operation with a lattice translation vector. The full symmetry group can be decomposed into left cosets of the subgroup as

(2) |

where

(3) |

In the above equation, is one of the symmetry operations that maps to its symmetry-equivalent point as

(4) |

Here corresponds to the original point (i.e. and , where denotes the identity operation). These points constitute a crystallographic orbit whose multiplicity is given by , where is the number of symmetry operations in the full crystal group without pure translations, and is the number of elements in . The vector is chosen so that transforms to the point which lies in the unit cell.

From the site symmetry group for a given point , the symmetry-adapted Wannier functions centered at are defined as the basis functions of the irreducible representations of the site symmetry group ; these Wannier functions are represented as , where labels the irreducible representations and runs over the basis functions of the irreducible representation , and is the dimension of the irreducible representation . For these Wannier functions transform as

(5) |

From we can generate Wannier functions centered at as

(6) |

Therefore the symmetry-adapted Wannier functions can be specified by one representative point of their centers (i.e. Wyckoff position) and the irreducible representations of the corresponding site-symmetry group. If the irreducible representation is real, the corresponding Wannier functions can be chosen to be realCloizeaux63 (). From these symmetry-adapted Wannier functions , one can construct the Bloch functions as

(7) |

(8) |

Here is the number of -points.

To understand how these Wannier functions transform with respect to the operations in the full symmetry group of the system, the theory of induced representations is used; it can be shown that, for a given and any element of the full symmetry group, , there exist one pair of and that satisfy the following equationEvarestov97 ():

(9) |

where

(10) | ||||

(11) |

Using Eq. (9), it can be shown that these symmetry-adapted Wannier functions and the Bloch functions transform asEvarestov97 ()

(12) |

(13) |

In the above equations the index or the symmetry operation is determined by and according to Eq. (9). By writing , Eq. (13) can be rewritten as

(14) |

where is a block diagonal matrix

(15) |

Due to its block-diagonal form, the matrix can contain blocks corresponding to nonequivalent Wannier centers(different Wyckoff positions), and the number of the blocks in is given as the sum of the number of the irreducible representations considered for a given set of Wannier centers; when there are more than one set of Wannier functions belonging to the same irreducible representation, contains the same number of identical blocks as the number of these multiple sets.

#### ii.1.2 Construction of the Bloch functions from the eigenstates of the Hamiltonian

Our objective is to construct the Bloch wavefunctions defined in Eq. (7), which are related to the symmetry-adapted Wannier functions through unitary transformation and transform according to Eq. (14), from the linear combination of the eigenfunctions of some one-particle Hamiltonian that is invariant under the full symmetry operations of the system. In this work we use the Kohn-Sham Hamiltonian of density functional theoryDFT (). By using the Kohn-Sham wavefunctions , we construct the orthonormal Bloch functions as

(16) |

Since the Kohn-Sham wavefunctions form the basis of the irreducible representations of the full symmetry group of the system , they transform as

(17) |

for . From Eqs. (14), (16), (17), one can obtain the following relation between the transformation matrices and :

(18) |

Therefore can be calculated from by providing and . In this work we do not consider time-reversal symmetry, but Eq. (18) can be generalized to include it. For the symmetry operations that transform to itself, namely, for the operations in the little group of denoted by , Eq. (18) yields the condition that has to fulfill:

(19) |

where unitarity of is used.

Equations (12), (13), (18), and (19) are the central equations of this work. We force the transformation matrix to follow Eq. (18) for all , which guarantees that the resulting Wannier functions transform according to Eq. (12). Since and are related by Eq. (18), we need to calculate the transformation matrix only for -points inside the IBZ, which reduces the computational cost. Note that since in practice the Wannier functions are calculated using a limited subspace spanned by a finite number of the Kohn-Sham states inside a chosen “energy window”, it is not possible to construct for any desired irreducible representations. If a given irreducible representation is not compatible with the symmetry of the Kohn-Sham states inside the energy window, Eq. (19) cannot be fulfilled.

### ii.2 Maximally-localized Wannier functions

In the maximally-localized Wannier function approachMarzari97 (); Souza01 (), the Wannier functions are obtained by minimizing the spread functional

(20) |

where is the Wannier function whose center is in the cell . The matrix elements and are calculated asMarzari97 (); Blount ()

(21) | ||||

(22) |

where is the cell-periodic part of the Bloch function (Eq. (16)) and is the vector that connects a given -point with its neighbors and is its weight. The spread functional given by Eq. (20) can conveniently be decomposed as , whereMarzari97 ()

(23) | ||||

(24) | ||||

(25) |

It can be shown that is invariant under the unitary transformation of the Bloch functionsMarzari97 (). The algorithm to minimize the spread functional is given in Refs. Marzari97, and Souza01, for both the cases where the bands of interest are isolated from other bands and are entangled with other bands. The minimization can be done as a postprocess to density-functional calculations, and the necessary input data are , the overlap matrix elements between the states at and from which the spread functional is calculated via Eqs. (21) and (22), and the initial guess of the transformation matrix, which is obtained by orthonormalizing the following matrixMarzari97 ()

(26) |

where is an initial guess of the Wannier function .

### ii.3 Minimization of the spread functional under symmetry constraint

#### ii.3.1 Input data

To perform minimization of the spread functional under the symmetry constraint (Eqs. (18) and (19)), in addition to the overlap matrix elements and the initial guess of the transformation matrices(Eq. (26)), one needs the matrix representation of the symmetry operations in the basis of the Bloch functions defined in Eq. (7) and the Kohn-Sham states, (Eq. (15)) and (Eq. (17)). The former is in many cases obtained by specifying the center and the character of the Wannier functions (ex. , , ), and calculating the rotation matrix for the elements of the corresponding site-symmetry group expressed in the basis of these functions. The latter is calculated from the Kohn-Sham wavefunctions as

(27) |

Similar to the overlap matrices , in Eq. (27) can be calculated with any basis set, and after all data are calculated the procedure is basis independent, as in the original maximally-localized Wannier function approach. We also note that this procedure does not require any specific phase factor relation between and ; the Kohn-Sham wavefunctions at can be calculated independently from those at , or can be generated from the wavefunctions at by performing symmetry operations. It is also important to include all degenerate states in the calculation of , inside the specified energy window.

The initial transformation matrix has to follow Eq. (19). We can construct that fulfills this requirement iteratively as follows; starting from that is calculated from the initial guess of the Wannier functions (Eq. (26)), we first calculate

(28) |

and in the next step this is orthonormalized by, e.g. using singular value decomposition. This cycle is repeated until we get converged . For a limited energy window, it is not always possible to construct for a given set of the irreducible representations. A measure to check the convergence of can be

(29) |

which is zero if fulfills Eq. (19).

#### ii.3.2 Isolated set of bands

In the case where we construct Wannier functions from bands that are separated from all other bands, since any unitary transformation of the Bloch states does not change , we only have to consider the variation of with respect to the change

(30) |

where is an infinitesimal antiunitary matrixMarzari97 (). Using the relation Eq. (18), for the gradient of the spread functional is calculated as

(31) |

where is the number of the symmetry operations that leave unchanged, and is the gradient of with respect to , whose explicit form is given in Ref. Marzari97, . It can be shown the new set of optimized along this direction also satisfy Eq. (19). The transformation matrices for -points not inside the IBZ are obtained via Eq. (18).

#### ii.3.3 Entangled bands

When we construct the Wannier functions from the states that are entangled with other bands, generally the number of the Bloch states inside a given energy window is larger than the number of the Wannier functions . Following Ref. Souza01, , in this case we minimize the spread functional using the two-step procedure; first we determine the optimal subspace inside the specified energy window spanned by orthonormal states that minimizes , and in the second step the remaining part of the spread functional, , is minimized within the chosen subspace. In the first step, we search for the optimal wavefunctions

(32) |

which minimize and also transform according to Eq. (13). This set of wavefunctions are obtained by the variation of

(33) |

where is a Lagrange multiplier and is the periodic part of . In this work we calculate by using the steepest descend method; in each iteration the wavefunctions are minimized along the direction

(34) |

where is the Hermitian operator defined as

(35) |

Here is the projection operator defined asSouza01 ()

(36) |

and is calculated as . In practice, for each state , we diagonalize in the subspace spanned by and , and we construct a new from the eigenvector with the larger eigenvalue of this matrix. In each iteration, after all are updated, we orthonormalize and impose the condition Eq. (19) by using the method described above. This point is an important difference between the current scheme and the usual maximally-localized Wannier function approach; in the latter is chosen to be the eigenvectors of the largest eigenvalues of (Eq. (36)). The optimal subspace chosen in the conventional approach does not necessarily match the subspace spanned by the desired symmetry-adapted Wannier functions.

After the wavefunctions are obtained, we calculate the transformation matrix that yields the Bloch wavefunctions(Eq. (7)) as a linear combination of ,

(37) |

which yields the minimum of in the chosen subspace, in the same way as in the isolated band case. Since and both transform according to Eq. (13), the relation of the transformation matrices and is modified from Eqs. (18) and (19) as follows:

(38) |

(39) | ||||

### ii.4 Computational details

In this work we perform calculations using the plane-wave DFT code TAPPYamauchi96 () with norm-conserving Troullier-Martins type pseudopotentialsTroullier91 (). We employ the generalized gradient approximationPerdew96 () for the exchange-correlation functional. For the minimization of the spread functional, the routines in Wannier90 librarywannier90 () are used. All calculations are done using experimental lattice constants, that are and for GaAs and Cu, respectively, and we use and k-point sampling including (the point). Energy cutoffs of the plane-wave basis are 25 Ry and 64 Ry for GaAs and Cu, respectively. Spin-orbit coupling is not included in the calculations.

## Iii results

### iii.1 GaAs

First we consider constructing four Wannier functions from the four valence bands in GaAs, whose band structure is shown in the solid lines in Fig. 1. As shown in Ref. Marzari97, , in this system the maximal localization procedure yields four localized functions centered on four covalent bonds. From group-theoretical view, those correspond to the irreducible representation of site-symmetry group of the Wyckoff position , and one can obtain the same results with our symmetry-constrained minimization procedure. In this system any point along the bond yields the same set of the matrices (Eq. (15)), therefore starting from the initial Wannier functions centered at an arbitrary point along the bond and its symmetry-equivalent three points, after the minimization their centers are moved to the points which yield the minimum of the spread functional, that are around away from Ga atom.

Another set of symmetry-adopted Wannier functions that are compatible with the symmetry of these four valence states are -like and -like functions centered at the anion(As) atom, which correspond to the irreducible representations and , respectively, of site-symmetry group of the Wyckoff position . In the usual minimization without symmetry constraint, these atom-centered Wannier functions are a stationary point of the spread functional but not the global minimum of it, and therefore they are not a stable solution as discussed by Marzari and VanderbiltMarzari97 (). In our approach, these atom-centered Wannier functions are easily obtained by providing corresponding matrix . In Table 1 we compare the spreads of the bond-centered and atom-centered Wannier functions. Following Ref. Marzari97, , in the table we also show the results obtained by combining two independent calculations for - and -like Wannier functions; namely, in this calculation we first perform two calculations to obtain the -like and -like Wannier functions separately, from the lowest band and higher three bands, respectively. By using these and transformation matrices we construct the transformation matrix in the block-diagonal form without further optimization. As anticipated, compared to this separate result, the atom-centered Wannier functions constructed with four valence bands are more localized. This is mainly due to the reduction in the off-diagonal contribution of the spread functional (Eq. (24)).

In this system it is also possible to construct -like and - functions centered at the cation(Ga) atom from the four valence bands. These - and -like functions correspond to the irreducible representations and , respectively, of site-symmetry group of Wyckoff position . The spreads of these cation-centered functions are also shown in Table 1, and as anticipated, these Wannier functions are more delocalized compared to bond-centered or anion-centered ones. Unlike the As-centered case, it is not possible to construct these - and -like Ga-centered functions separately from the lowest band and other three bands; Equation (19) cannot be fulfilled separately for and unitary transformation matrices for the lowest band and higher three bands, respectively, but it can be fulfilled if we use the four valence bands together. The reason for this becomes clear from Fig. 1, where we plot the Wannier-interpolated band-structureSouza01 () calculated separately from the -like Wannier function and from the three -like Ga-centered Wannier functions. One can see the Ga-centered -like Wannier function is connected to and states, not the lowest and states which are constructed from the -like functions. This shows a close connection between the symmetry of the Wannier functions and the band structure; the correspondence between irreducible representations of a given site symmetry group and the Bloch functions at high-symmetry -points can be found in the tables in Ref. Evarestov97, . At the point, there are two states belonging to the same irreducible representation () that contribute to both the -like and -like Wannier functions, therefore at this point the two interpolated bands deviate from the original ones. We finally note that the original band structure of the four valence bands can be reproduced by using these -like and -like Ga-centered Wannier functions together in the interpolation.

bond-centered | ||||||
---|---|---|---|---|---|---|

6.124 | 0.006 | 0.630 | 6.760 | 1.690 | ||

7.870 | 0.006 | 0.566 | 8.442 | 2.110 | ||

centered on As | ||||||

6.124 | 0.012 | 3.502 | 9.639 | 1.450 | 2.730 | |

7.870 | 0.012 | 3.826 | 11.708 | 1.510 | 3.399 | |

*centered on As | ||||||

6.124 | 0.064 | 4.388 | 10.576 | 1.828 | 2.916 | |

7.870 | 0.069 | 4.943 | 12.882 | 2.032 | 3.617 | |

centered on Ga | ||||||

6.124 | 0.151 | 7.648 | 13.924 | 2.448 | 3.825 | |

7.870 | 0.112 | 9.028 | 17.011 | 2.615 | 4.798 |

### iii.2 Cu

Energy window [-10 eV :10 eV] | |||||||
---|---|---|---|---|---|---|---|

atom-centered | |||||||

3.901 | 0.000 | 0.248 | 4.149 | 1.959 | 0.439 | 0.437 | |

5.613 | 0.000 | 0.208 | 5.820 | 3.474 | 0.466 | 0.475 | |

centered at | |||||||

3.555 | 0.000 | 0.407 | 3.962 | 1.706 | 0.447 | 0.457 | |

4.815 | 0.000 | 0.561 | 5.376 | 2.855 | 0.488 | 0.528 | |

centered at | |||||||

3.279 | 0.167 | 0.500 | 3.946 | 1.599 | 0.489 | 0.441 | |

3.968 | 0.107 | 0.511 | 4.587 | 2.042 | 0.534 | 0.471 | |

Energy window [-10 eV :20 eV] | |||||||

atom-centered | |||||||

3.342 | 0.000 | 0.041 | 3.382 | 1.463 | 0.386 | 0.381 | |

3.685 | 0.000 | 0.016 | 3.701 | 1.705 | 0.400 | 0.397 | |

centered at | |||||||

3.059 | 0.000 | 0.164 | 3.222 | 1.247 | 0.393 | 0.398 | |

3.347 | 0.000 | 0.176 | 3.523 | 1.439 | 0.412 | 0.424 | |

centered at | |||||||

2.947 | 0.011 | 0.250 | 3.208 | 1.164 | 0.422 | 0.390 | |

3.168 | 0.008 | 0.252 | 3.428 | 1.284 | 0.443 | 0.407 |

Next we consider constructing six Wannier functions from one -like and five -like states for bulk copper in fcc structure. Souza et al.Souza01 () showed that in this system one obtains five -like Wannier functions centered on Cu atom that are split into and states and one -like Wannier function whose center is not on Cu atom but at the tetrahedral interstitial site . In the six-band case, this tetrahedrally-centered Wannier function is not regarded as a symmetry-adapted Wannier function, as due to the inversion symmetry this site (Wyckoff position ) is equivalent to and thus one needs one additional -like Wannier function centered at the latter site to make them the basis functions of the full symmetry group, resulting a seven-band model as discussed by Souza et alSouza01 ().

The possible -like symmetry-adapted Wannier function in this six-band case is (i) irreducible representation centered on Cu atom(Wyckoff position ) and (ii) irreducible representation centered at (Wyckoff position ). In Table 2, we compare the spreads of these sets of the six Wannier functions obtained by using the energy window of [-10eV:+10eV] and [-10eV:+20eV]. In these calculations, the centers of the and Wannier functions are on the Cu atom. As expected, the two symmetry-adapted solutions yield larger spreads than the tetrahedrally-centered Wannier functions obtained via the unconstrained minimization, while the spreads of -like Wannier functions do not vary very much in these three cases. The -like Wannier function centered at is found to be more localized than the atom-centered Wannier function, which may be traced back to the fact that the -like band in Cu is very extended and it has a larger weight in the interstitial region.

The gauge-invariant part of the spread functional (, Eq. (23)) is also different in the three cases, which indicates that the optimal subspace(Eq. (32)) chosen in the first step of the minimization procedure is different in these three cases as a result of the symmetry constraint. As reported previously by Souza et al.Souza01 () and Thygesen et al.Thygesen05 (), in the case of the energy window [-10eV:+10eV], we find that without symmetry constraint the atom-centered -like Wannier function is not a stable solution. In the case of the larger energy window of [-10eV:+20eV], we get the atom-centered -like Wannier function without symmetry constraint by using atom-centered gaussian functions as initial trial functions, however we find that without symmetry constraint this solution is unstable against a small perturbation of the initial states, and this clearly shows the importance of the symmetry constraint when constructing Wannier functions from extended bands. In Fig. 2 we plot these three -like Wannier functions. Since the two symmetry-adapted functions (Fig. 2 (a) and (b)) belong to the irreducible representation, they are invariant with respect to transformations of their site symmetry group. It can be seen that the most localized tetrahedrally-centered solution (Fig. 2 (c)), which is obtained without symmetry constraint, is also symmetric with respect to the rotation of around the [111] axis, which indicates this function also reflects some site-symmetry properties of the tetrahedral site.

In Fig. 3 we plot Wannier-interpolated band structures calculated with different sets of the Wannier functions. For completeness, in Fig. 3(d) we also show the interpolated band structure calculated with seven Wannier functions, namely, five -like Wannier functions and two equivalent -like Wannier functions centered at (). As in the case of GaAs and as also discussed by Souza et al.Souza01 (), for high-symmetry points in the Brillouin zone one can predict which Bloch states can be formed from a given set of the symmetry-adapted Wannier functions; as seen in Fig. 3 the atom-centered -like Wannier function (Fig. 3(a)) is connected to , , and states, while from the Wannier function centered at , , , and states are formed. The low-lying and states are formed with the tetrahedrally-centered Wannier functions (Figs. 3(c) and (d)), as discussed by Souza et alSouza01 ().

As can be seen in Fig. 3(a), the , , and states which are formed by the atom-centered -like Wannier function are located in a relatively high energy region, and this is why the atom-centered -like Wannier function is unstable in the conventional maximal localization approach when a small energy window is chosen. Indeed, in our calculation, with a smaller choice of the energy window, we cannot satisfy the relation Eq. (19) for the atom-centered -like Wannier function. This shows the importance of selecting the energy window properly, as the symmetry of the Wannier functions are determined by the symmetry properties of the Bloch functions inside the energy window through Eqs. (18) and (27).

## Iv conclusions

In this paper we have presented a systematic procedure to generate symmetry-adapted Wannier functions based on the theory of site-symmetry and induction group combined with the maximally-localized Wannier function approach. This scheme can easily be implemented in the existing maximally-localized Wannier function calculation code, and it allows one to calculate localized functions of a specified symmetry which do not necessarily yield the global minimum of the spread functional. It also provides the relation between the unitary transformation matrices for symmetry-equivalent -points, which simplifies the minimization process and also improves accuracy of the calculation.

The results for GaAs and Cu show that the calculated Wannier functions are indeed localized and have the specified symmetry properties, and they reflect the symmetry of the Bloch functions inside the energy window used in the calculation.

These symmetry-adapted Wannier functions are suitable for symmetry analysis of the band-structure of the system and for accurate basis functions of the tight-binding model. Generalizations of the present method, such as including spin-orbit coupling, would be interesting subjects to be investigated.

###### Acknowledgements.

We thank K. Nakamura for providing us with his pseudopotential data, and F. Aryasetiawan and F. Nilsson for helpful discussions. This work was supported by Swedish Research Council and the Scandinavia-Japan Sasakawa Foundation.## References

- (1) N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
- (2) I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 65, 035109 (2001).
- (3) N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012).
- (4) E. I. Blount, Solid State Phys. 13, 305 (1962).
- (5) K. S. Thygesen, L. B. Hansen, and K. W. Jacobsen, Phys. Rev. B 72, 125119 (2005).
- (6) J. des Cloizeaux, Phys. Rev. 129, 554 (1963).
- (7) W. Kohn, Phys. Rev. B 7, 4388 (1973).
- (8) J. von Boehm and J.-L. Calais, J. Phys. C: Solid State Phys. 12, 3661 (1979).
- (9) E. Krüger, Phys. Rev. B 36, 2263 (1987).
- (10) B. Sporkmann and H. Bross, Phys. Rev. B 49, 10869 (1994).
- (11) R. A. Evarestov and V. P. Smirnov, Site Symmetry in Crystals: Theory and Applications, 2nd ed., Springer Series in Solid State Sciences (Springer-Verlag, New York, 1997).
- (12) V. P. Smirnov and D. E. Usvyat, Phys. Rev. B 64, 245108 (2001).
- (13) V. P. Smirnov and R. A. Evarestov, Phys. Rev. B 72, 075138 (2005).
- (14) M. Posternak, A. Baldereschi, S. Massidda, and N. Marzari, Phys. Rev. B 65, 184422 (2002).
- (15) S. Casassa, C. M. Zicovich-Wilson, and C. Pisani, Theor. Chem. Acc. 116, 726 (2006).
- (16) R. M. Dreizler and E. K. U. Gross, Density Functional Theory (Springer-Verlag, Berlin, 1990).
- (17) J. Yamauchi, M. Tsukada, S. Watanabe, and O. Sugino, Phys. Rev. B 54, 5586 (1996).
- (18) N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
- (19) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- (20) A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Comput. Phys. Commun. 178, 685 (2008).