Quantum Chemistry of Solids / 13-Symmetry and Localization of Crystalline Orbitals

3.3 Symmetry of Localized Crystalline Orbitals. Wannier Functions

97

where Uγmµ , ij (k) ≡ Uσ , τ (k) is a unitary matrix, which in the case of a nondegenerate band is reduced to a phase factor exp(iα(k)). It follows from (3.110) and (3.111)

Wij(β)(r − an) = N − 21
			Uγmµ, ij (k)· exp(−ikan) · ϕmµ(γ) (k, r)	(3.112)
		k γmµ
or
Wt(r) = N − 21			with Ust = Uγmµ, ij (k)· exp(−ikan)
	Ust · ϕs(r)			(3.112a)

s

If the functions ψij(β)(k, r) satisfy the conditions of orthonormality

(2π)3			(ψij	(k, r)) ψi j (k , r)dr = δii δjj δ(k − k )	(3.113)
	Va		(β)	(β)
the functions Wij(β)(r − an) form an orthonormal system:
(Wij(β)(r − an), Wi(βj) (r − an )) = δii δjj δnn or (Wt(r), Wt (r)) = δtt					(3.114)

Let Q(q,β) be the space of a simple indrep (q, β) (for simplicity). The space Q(q,β) is spanned by both the set of orthonormal functions Wij(β)(r − an) ≡ Wt(β)(r) and the set

of Bloch functions ϕ(mµγ) (k, r), = ϕs(r), see (3.112). The orthonormal functions Wt(β)(r) in (3.112) can be chosen real (if irrep β is real) and transform according to the irrep β of the site groups Lqj,n of the points qj,n. Their localization depends on the choice of the matrix Uγmµ, ij (k) in (3.112). The existence of Wannier functions decreasing exponentially at infinity (for the model of an infinite crystal) is established in many cases. The uniqueness of these functions for nondegenerate bands in crystals with centers of inversion was proved [44,57]. In this case Wannier functions correspond to a

special choice of phase factors Uγmµ, ij (k) = exp(iα(k)) of Bloch orbitals ϕ(mµγ) (k, r) = ϕs(r) in (3.112). Any other choice of phase factors destroys either the symmetry properties of the Wannier functions, or their reality, or both [42]. Obviously these functions are as well localized as possible (fall o exponentially in the model of infinite crystal). If the choice of phase factors is not correct the Wannier functions lose the exponential character of their decrease and, therefore, are not maximally localized according to any reasonable criterion of localization. Unfortunately, the uniqueness of Wannier functions is not yet proved for the more general case of the degenerate bands in crystals with centers of inversion, where it apparently exists.

As a criterion of localization for a localized function W (r) one uses the value of the integral over the whole space of the crystal [42]

I = ρ(r)|W (r)|2dr

(3.115)

with the weight function ρ(r) ≥0, which is supposed to be invariant under the operations from the site-symmetry group Gq. Particular choices of the weight function are:

98	3 Symmetry and Localization of Crystalline Orbitals
	1.
	ρ(r) = (r − q)2	Boys localization criteria		(3.115a)
	2.	1,	if r ∆;
	ρ(r) =			(3.115b)
		0,
			if r / ∆

where ∆ is some region surrounding the point q of symmetry localization of the
function W	(β)	(r);
3.	0
				(r − q)2
		ρ(r) = (πr2)−3/2 exp(	−		)	(3.115c)
		o		ro2

which underlines the contribution in (3.115) of the values of the function |W0(β)(r)|2

inside the sphere of radius ro.		As a special		case of (3.115c) for ro		→ 0 one has
ρ(r) = δ(r − q).				(β)	(r − an) in the
One searches for the set of nonorthogonal localized functions Vij
form
Vij(β)(r − an) = N − 21	k	Cγmµ, ij (k)· exp(−ikan) · ϕmµ(γ) (k, r)				(3.116)
		γmµ
or
	Vt(β)(r) = N − 21			Cst · ϕs(r)		(3.116a)
			s

The system of functions Vt(β)(r) ≡ Vij(β)(r − an) can be obtained from the function V11(β)(r) ≡ V0(β)(r) in the same way as the functions Wij(β)(r − an) from the function

W	(β)	(r) (see above). Therefore, it is su cient to find only one function, for example
	11
		V0(β)(r) = N − 21
			Cs0 · ϕs(r)	(3.117)

s

The coe cients Cs0 can be found from the following variational problem: to find the coe cients Cs0 in (3.117), maximizing (or minimizing) the functional (3.115) and satisfying the supplementary condition:

|V0(β)(r)|2dr =1

(3.118)

This variational problem is equivalent to the eigenvalues and eigenvectors problem

	for the matrix:	1
	Ass =			ρ(r)ϕs (r)ϕs (r)dr	(3.119)
		N

The eigenvalues of the matrix A are stationary values of the localization criterion (3.115), and the eigenvectors corresponding to these values are required coe cients of the expansion (3.117). In our case, it is necessary to search for the eigenvector corresponding to the highest eigenvalue for the choices (3.115b) and (3.115c) of the weight function ρ(r) and to the lowest one for the case (3.115a). Let us note that it is su cient to use the variational procedure in the subspace of the first basis vectors of

3.3 Symmetry of Localized Crystalline Orbitals. Wannier Functions

99

the irrep β of the site group Lq instead of the whole space of the canonical orbitals of the energy band under consideration.

Though the set of the functions Vt(r) is not orthogonal, these functions are close to the accurate localized Wannier functions Wt(β)(r) ≡ Wij(β)(r − an). They can be chosen real (for real irrep β of the site group Lq) and satisfy all the symmetry requirements for the functions Wt(β)(r). The orthonormal system Wt(β)(r) is generated

	(β)				orthogonalization procedure for
	from the functions Vt (r) by a suitable symmetrical
						!
	periodic systems:	!	t	(S− 21 )t tVt(β)(r)
		Wt(β)(r) =				(3.120)

where S is the overlapping matrix of the functions Vt(β)(r):

	Stt ≡ (Vt(β)(r),Vt(β)(r))
Combining (3.116a) and (3.120) we get
!	(r) = N − 21	t		(S− 21 )
Wt(β)		Cst	·		·	ϕs(r)
				t t

s

(3.121)

(3.122)

As the symmetrical orthogonalization procedure (3.122) leaves unchanged the reality and symmetry properties of the functions, the set of orthonormalized functions

!(β)

Wt (r) satisfy all the requirements to the localized Wannier functions (reality, symmetry requirements and orthonormality) and (in the case when this set of functions is unique) has to coincide with the latter:

W (β)	(β)	(r)	(3.123)
!t	(r) = Wt

The weight function ρ(r) in the functional I, see (3.115), can be varied. In particular, one can choose for the region ∆ in (3.115b) a mu n-tin sphere, some part of a Wigner–Seitz cell (even very small),etc . The functions Vt(r) depend on the choice of ρ(r), but, according to calculations, after the procedure of symmetric orthogonalization we always have the same result even for degenerated bands in crystals with a center of inversion. Apparently, the proposed method gives, in these cases, the orthonormal set of maximally localized Wannier functions. Thus, the numerical calculations imply constructing matrix Ass (3.119), diagonalization of this matrix, then obtaining the overlap matrix Stt (3.121), taking the matrix square root of it and other operations of linear algebra.

To demonstrate the reliability of the proposed variational method let us consider two examples of its applications given in [42] – the Wannier-function generation in silicon and MgO crystals.

In accordance with the theory of induced (band) representations the corresponding Wannier functions in the silicon crystal (four per unit cell) are centered at the middle of the bonds between the nearest Si atoms (Wycko position c with site group Gc = D3d) and transform according to the irrep a1g of the site group D3d, see Sect. 3.3.1.

For the variational procedure two sets of Bloch functions were used (both obtained with the help of CRYSTAL code [23]). The first set (S1) corresponds to the full electron Hartree–Fock LCAO calculations, the second one (S2) – to the pseudopotential

100 3 Symmetry and Localization of Crystalline Orbitals

Hartree–Fock LCAO method. A model of a finite crystal with periodical boundary conditions (cyclic model) with the main region composed of 4 × 4 × 4 = 64 primitive cells was adopted. The first basis consists of 13 s- and p- atomic-like functions per atom, the pseudopotential basis consists of 2 s-, 6 p- and 5 d-functions per atom. The weight function ρ(r) = δ(r − q) has been taken in (3.115).

Table 3.17. Localized states for the upper valence band of Si crystal along [1,1,1] direction

√

(the origin is taken in the middle of the bond, x in a 3/96 units, a being the conventional lattice constant), [42]

x	W1(x)	W8(x)	W64(x)	W64(no)(x)
0	-0.2300	-0.2036	-0.1982	-0.1988
3	-0.2190	-0.1975	-0.1926	-0.1930
6	-0.1566	-0.1465	-0.1436	-0.1437
9	0.1204	0.1216	0.1199	0.1196
11	0.4116	0.3840	0.3725	0.3747
12	-0.5831	-0.7333	-0.7579	-0.7406
15	-0.0596	-0.0401	-0.0354	-0.0373
18	0.0497	0.0336	0.0288	0.0307
22	0.0620	0.0501	0.0453	0.0471
27	0.0250	0.0327	0.0317	0.0321
33	-0.0099	0.0097	0.0127	0.0121
39	-0.0098	-0.0002	0.0040	0.0030
45	0.0045	-0.0017	0.0011	0.0004
51	0.0045	-0.0002	0.0006	0.0003
57	-0.0098	0.0005	0.0003	0.0003
63	-0.0099	0.0000	0.0000	-0.0001
69	0.0250	-0.0006	-0.0002	-0.0002
75	0.0649	-0.0015	-0.0002	-0.0001
81	-0.0596	0.0008	0.0000	-0.0001
84	-0.5831	-0.0082	-0.0005	0.0000
87	0.1204	-0.0001	0.0000	0.0000
90	-0.1566	-0.0002	-0.0001	-0.0002
93	-0.2190	-0.0006	-0.0002	-0.0004
96	-0.2300	-0.0008	-0.0003	-0.0006

The bigger the cyclic model of a crystal, the closer are its Wannier functions to those of the model of an infinite crystal. In [42] the convergence of the Wannierfunctions W (r) (for the set S1) for di erent sizes of cubic cells defining the cyclic model of a crystal was studied. The smallest cubic cell consists of four primitive cells, the larger ones consist of 8 and 64 smallest cubic cells. The total number of wavevectors used in BZ summation is four times larger than the number of smallest cubic cells in the supercel. The localized functions W1(x), W8(x), W64(x) were calculated according to (3.112) as a linear combinations of 4 (supercell coincides with a smallest cubic

3.3 Symmetry of Localized Crystalline Orbitals. Wannier Functions

101

cell), 32 (2 × 2 × 2 = 8 cubic cells in a supercell) and 256 (4 × 4 × 4 = 64 cubic cells in a supercell) Bloch states. Table 3.17 gives the values of the Wannier functions W1(x), W8(x), W64(x) and W64(no)(x) (the latter is nonorthogonalized) in some points along the [1, 1, 1] direction. As is seen from the table the function W64(no)(x) is very close to W64(x). So the orthogonalization procedure changes the Wannier function insignificantly. When comparing the functions W1(x), W8(x), W64(x) it is necessary to take into account that these functions are normalized di erently, in the volume of

1, 8 and 64 cubic cells, respectively. The fact that their values are relatively close to

√

each other from x = 0 to x = 48 (in units of a 3/96, the translational period in the

√

[1, 1, 1] direction is equal to a 3) shows a good localization of the Wannier function under consideration.

Figure 3.9 gives the W64(x) and pseudopotential(pp) Wannier function W64(pp)(x) in the [111] direction.

The functions di er significantly at the atom cores. Outside the cores the behavior of the functions is alike. This is quite natural, since the pseudopotential Wannier functions are constructed from smooth pseudowave Bloch functions, that are nonorthogonal to the core states. Orthogonalization of these smooth Wannier functions to the localized core functions would lead to an increase of oscillations at atom cores observed in the Wannier functions W64(x). The di erences between the Wannier functions in the interatomic space are due to the normalization.

In this case, the method reveals a very good stability with respect to the choice of the weight function ρ(r) form. The computations give the same resulting orthogonalized Wannier functions, whereas the intermediate nonorthogonalized functions Vt(r) turn out to be di erent. An amazing feature was noted. Even if one uses the weight functions ρ(r) centered at any point in the localization region of the Wannier function and so the functions Vt(r) do not have the symmetry compatible with the site group Gq, the same orthogonalized Wannier functions arise after the symmetrical orthogonalization procedure. The latter not only conserves the symmetry of the localized orbitals, but reconstructs(!) it up to an appropriate level. The cause of such flexibility is the fact that the most localized Wannier functions corresponding to a certain energy band are unique and just these functions arise without fail as the result of the proposed variational procedure [42]).

To construct the Wannier functions for the upper valence bands of the perfect MgO crystal a set of MgO Bloch functions has been applied [42], that were obtained in pseudopotential Hartree–Fock LCAO calculations with the CRYSTAL code [23]. The 4 × 4 × 4 cyclic system has been used in the calculations. The valence band of MgO represents two separated bands, and thus it is possible to construct two independent sets of Wannier functions (for each of the bands). The method of induced representations gives all the Wannier functions being centered on O atoms (Wycko position b with site group Gc = Oh) and transforming – according to the irrep ag of the site group Oh for the lower band (one Wannier function per unit cell) and according to the irrep t1u – for the upper band (three Wannier functions per unit cell). The weight function ρ(r) = δ(r − q) was taken in (3.115) located some distance away from the centering point of the corresponding Wannier function – as the Wannier functions of the upper band are antisymmetric and thus equal zero in their centering points, while, as stated above, the symmetrical centering of ρ(r) function is not necessarily needed.

Figure 3.10 shows the Wannier functions for both the bands in the [111] direction.

102 3 Symmetry and Localization of Crystalline Orbitals

D

D

6L 6L

> @ D[LV

E

D

6L 6L

> @ D[LV

Fig. 3.9. Wannier functions for the silicon upper valence band along the [111] direction: a) full electron calculation; b) pseudopotential calculation

The function corresponding to the upper band (dashed curve) is determined up to a linear combination of the three functions transforming via t1u irrep. One can see that both Wannier functions are almost completely localized around one of the oxygen atoms, which confirms the ionic character of MgO compound.

The character of the localization of Wannier functions depends on the analytical properties of Bloch states (as a function of the wavevector) that are essentially determined by the nature of the system under consideration. One can arbitrary change only the form of an unitary transformation of Bloch functions. It is just this arbitrariness that is used in the variational approach [42] to assure the best localization of Wannier functions. The accuracy of the Wannier functions obtained by the proposed method is determined solely by the accuracy of the Bloch functions and the size of the supercell used. As the calculations have shown, the proposed method is reliable and useful in the problem of generation of the localized Wannier functions. In the two examples

3.3 Symmetry of Localized Crystalline Orbitals. Wannier Functions			103
D	0J	2

> @ D[LV

Fig. 3.10. Wannier functions for the MgO upper valence band along [111] direction: solid line for WF of a1g symmetry, dotted line – for one of the three WFs of t1u symmetry

considered the Bloch states of only occupied bands were used for WF generation. Inclusion of the vacant band states of the appropriate symmetry allows Wannier-type atomic orbitals (WTAO) to be generated. WTAO are used for the chemical-bonding analysis in crystals. We consider WTAO applications in Chap. 9.