Localized Orbitals, Wannier Functions

The representation of the electronic ground state in terms of localized Wannier orbitals [84] provides a powerful tool in the study of periodic solids. Recent advances in the formulation of a theory of electronic polarization [78] and the development of linear-scaling methods [60] have rejuvenated the use of Wannier functions as an analysis tool. Namely, Wannier functions afford an insightful picture to the nature of chemical bonding and aid in the understanding of classical chemical concepts (e.g. nonbonding electron pairs or valency) in terms of quantum mechanics.

Wannier functions (WF) are defined in terms of a unitary transformation performed on the occupied Bloch orbitals (BO) [84]. One major problem in a practical calculation is their non-uniqueness. This is a result of the indeterminacy of the BO's, which are, in the case of a single band, only determined up to a phase factor, in the multi-band case, up to an arbitrary unitary transformation among all occupied orbitals at every point in the Brillouin zone. As proposed recently by Marzari and Vanderbilt [85], one can resolve this non-uniqueness by requiring that the total spread of the localized function be minimal. This criterion is in close analogy with the Boys-Foster method [86] for finite systems, here one uses the spread defined through the conventional position operator. The new technique has been successfully applied to crystal systems and to small molecules within a general k-point scheme[85,87]. An extension to disordered systems within the ${\bf\Gamma}$ -point approximation was recently performed[88]. This is of particular interest when one would like a localized orbital picture within the framework of Car-Parrinello molecular dynamics (CPMD). Here we examine the problem focusing on the ${\bf\Gamma}$ -point approximation only. Upon minimization of the spread functional the appropriate unitary transformation to the localized orbitals can be calculated. With explicit knowledge of the spread functional we can derive the complete expressions required to implement the iterative minimization procedure.

We begin by reviewing the work of Resta [89]. In his treatment, the fundamental object for studying localization of an electronic state within Born-Von Karman boundary conditions is the dimensionless complex number,

$$z=\int_{L}dx \;\exp(i2\pi x/L)\;\vert\psi(x)\vert^2 .$$ (156)

Here, $L$ is the linear dimension, and $\psi(x)$ denotes the wavefunction. By considering the definition of the spread of the wavefunction to be $\Omega=\langle x^2\rangle-\langle x\rangle^2$ , where $\langle\cdots\rangle$ denotes an expectation value, Resta has shown that to $O(1/L^2)$ the functional for the spread in one-dimension to be,

$$\Omega={1\over (2\pi)^2}\ln \vert z\vert^2.$$ (157)

One goal of this study is to generalize Eq. (156) to three-dimensions and obtain the appropriate generalization of Eq. (157). Thus, we choose to study the following dimensionless complex number within Born-Von Karman boundary conditions,

$$z_I=\int_{V}d{\bf r}\;\exp(i{\bf G}_I\cdot {\bf r})\;\vert\psi({\bf r})\vert^2 .$$ (158)

Here, $I$ labels a general reciprocal lattice vector, ${\bf G}_I=l_I{\bf b}_1+m_I{\bf b}_2+n_I{\bf b}_3$ , where ${\bf b}_{\alpha}$ are the primitive reciprocal lattice vectors, the integers $l$ , $m$ , and $n$ are the Miller indices, $V$ is the volume of the supercell, and $\psi({\bf r})$ denotes the wavefunction. We must find an appropriate function of the $z_I$ 's that gives the three dimensional spread in the case of an arbitrary simulation cell. We proceed by noting that in a molecular dynamics simulation the cell parameters (primitive lattice vectors) to describe systems of general symmetry are given by ${\bf a}_1$ , ${\bf a}_2$ and ${\bf a}_3$ . It is convenient to form a matrix of these cell parameters, $\stackrel{\leftrightarrow}{\bf h}=({\bf a}_1, {\bf a}_2, {\bf a}_3)$ where the volume $V$ of the simulation cell is given by the determinant of $\stackrel{\leftrightarrow}{\bf h}$ . It is also very useful to define scaled coordinates, ${\bf s}=\stackrel{\leftrightarrow}{\bf h}^{-1}\cdot {\bf r}$ that lie in the unit cube. In molecular dynamics simulations, this allows one to perform periodic boundary conditions for systems with general symmetry by first transforming to the unit cube, performing cubic periodic boundary conditions, and transforming back to the general cell with the action of $\stackrel{\leftrightarrow}{\bf h}$ One can also compute the reciprocal space vectors for systems of general symmetry with knowledge of the matrix of cell parameters. Thus, the $I$ -th reciprocal lattice vector,

$${\bf G}_I=2\pi\left(\stackrel{\leftrightarrow}{\bf h}^{-1}\right)^{\rm T}\cdot \hat{\bf g}_I .$$ (159)

Here, the superscript $\rm T$ denotes transposition, and $\hat{\bf g}_I=(l_I,m_I,n_I)$ is the $I$ -th Miller index. We then substitute this expression into Eq. (158) and use the definition of ${\bf r}$ to obtain,

$$z_I={\rm det}\stackrel{\leftrightarrow}{\bf h}\int^1_0d{\bf s}\;\... ...\right)\;\vert\psi(\stackrel{\leftrightarrow}{\bf h}\cdot {\bf s})\vert^2 .$$ (160)

Note that the exponential in Eq. (160) is independent of any coordinate system. Following Resta [89] we can write the electron density in terms of a superposition of localized density and its periodic images, $\vert\psi(\stackrel{\leftrightarrow}{\bf h}\cdot {\bf s})\vert^2=\sum\limits... ...}{\bf h}\cdot {\bf s}_0-\stackrel{\leftrightarrow}{\bf h}\cdot \hat{\bf m})$ . Here $\hat{\bf m}$ is a vector of integers and $\stackrel{\leftrightarrow}{\bf h}\cdot {\bf s}_0$ is the center of the distribution such that $\int_{-\infty}^{\infty}d{\bf s}\;\stackrel{\leftrightarrow}{\bf h}\cdot {\bf s}\;n_{\rm loc}(\stackrel{\leftrightarrow}{\bf h}\cdot {\bf s})=0$ . Using the Poisson summation formula [90], we rewrite Eq. (160),

$$z_I=\exp\left(i 2\pi\hat{\bf g}^{\rm T}_I\cdot {\bf s}_0\right)... ...}(-2\pi\hat{\bf g}^{\rm T}_I\cdot \stackrel{\leftrightarrow}{\bf h}^{-1}) ,$$ (161)

where $\hat{n}_{\rm loc}$ denotes the Fourier transform of $n_{\rm loc}$ . Furthermore, since we are considering $n_{\rm loc}$ to be localized, its Fourier transform is smooth over reciprocal distances and we can be assured that it is well represented about $\hat{g}_I=0$ . We expand $\hat{n}_{\rm loc}(-2\pi\hat{\bf g}^{\rm T}_I\cdot \stackrel{\leftrightarrow}{\bf h}^{-1})$ to second order, obtaining,

$$\hat{n}_{\rm loc}(-2\pi\hat{\bf g}^{\rm T}_I\cdot \stackrel{\le... ...al\hat{g}_{\alpha,I} \partial\hat{g}_{\beta,I}}\vert _{\hat{g}_I=0}+\ldots .$$ (162)

The second term in Eq.(162) is zero given our imposed condition $\langle\stackrel{\leftrightarrow}{\bf h}\cdot {\bf s}\rangle=0$ . Thus, we are left with,

$$\hat{n}_{\rm loc}(-2\pi\hat{\bf g}^{\rm T}_I\cdot \stackrel{\le... ...lpha s_\beta\;n_{\rm loc}(\stackrel{\leftrightarrow}{\bf h}\cdot {\bf s}) .$$ (163)

Combining Eq. (163) and Eq. (161), we obtain,

$$1-\vert z_I\vert=V {(2\pi)^2\over 2}\sum_{\alpha,\beta}\hat{g}_{... ...lpha s_\beta\;n_{\rm loc}(\stackrel{\leftrightarrow}{\bf h}\cdot {\bf s}) .$$ (164)

Keeping in mind that $\int_{-\infty}^{\infty}d{\bf s}\;\stackrel{\leftrightarrow}{\bf h}\cdot {\bf s}\;n_{\rm loc}(\stackrel{\leftrightarrow}{\bf h}\cdot {\bf s})=0$ , one can define the spread of the electronic distribution for the case of a general box through,

$$\langle r^2\rangle - \langle r\rangle^2 = \langle\left(\stackrel... ...}s_{\beta} \;n_{\rm loc}(\stackrel{\leftrightarrow}{\bf h}\cdot {\bf s}) .$$ (165)

Here, $g_{\alpha\beta}=\sum_{\mu}\stackrel{\leftrightarrow}{h}^{\rm T}_{\alpha\mu}\stackrel{\leftrightarrow}{h}_{\mu\beta}$ can be thought of as a metric tensor to describe the corresponding distances in the unit cube. Eq. (165) shows us exactly how the length scales are built into the spread through the metric tensor. From direct comparison of Eq. (164) and Eq. (165) we see that for supercells of general symmetry we need to choose linear combinations of $\hat{g}_{\alpha,I}\hat{g}_{\beta,I}$ that reproduce the metric tensor, $g_{\alpha\beta}$ . However, as stated earlier, $\hat{g}_{\alpha,I}$ are dimensionless numbers. Thus, an appropriate generalization takes the form of a sum rule,

$$g_{\alpha\beta}=\sum_I\omega_I\hat{g}_{\alpha,I}\hat{g}_{\beta,I} .$$ (166)

Here, $\omega_I$ are the ``weights'' with the appropriate dimensions to be determined later. Thus, it should also be clear that $g_{\alpha\beta}$ will have at most six independent entries (for triclinic symmetry) and thus a maximum of six weights are needed. It is interesting to note that by multiplying Eq. (166) on the left and right hand sides by $\stackrel{\leftrightarrow}{\bf h}^{-1}$ and using the definition of ${\bf G}_I$ , one will recover the rule used by Silvestrelli [91] and by Marzari and Vanderbilt [85]. Finally, we generalize to more than one state, $\vert\psi\rangle\rightarrow \vert\psi_n\rangle$ and the desired expression for the spread, $\Omega$ in a supercell of general symmetry is,

$$\Omega = {2\over (2\pi)^2}\sum_n^{\rm N_{states}}\sum_I\omega_I (1- \vert ... ... O(2\pi\hat{{\bf g}}^{\rm T}_I\cdot \stackrel{\leftrightarrow}{\bf h}^{-1})^2$$

$$z_{I,n} = \int_{V}d{\bf r}\;\exp(i{\bf G}_I\cdot {\bf r})\;\vert\psi_n({\bf r})\vert^2 ,$$ (167)

where Eq. (166) determine the ${\bf G}_I$ .

At this point it is useful to make contact with other spread formulas that are present in the current literature. Following Resta's derivation one finds the formula [89], that in our notation reads,

$$\Omega= - {1\over (2\pi)^2}\sum_n^{\rm N_{states}}\sum_I\omega_I\log \vert z_{I,n}\vert^ 2 ,$$ (168)

with $z_{I,n}$ defined as above. Eq. (168) is obtained by inserting Eq. (163) into Eq. (161), taking the $\log$ of the absolute value and expanding to consistent order.

Silvestrelli[91] on the other hand uses (again, in our notation),

$$\Omega= {1\over (2\pi)^2}\sum_n^{\rm N_{states}}\sum_I\omega_I (1- \vert z_{I,n}\vert^2) ,$$ (169)

with a similar definition for $z_{I,n}$ . Obviously Eq. (169) is obtained from Eq. (168) by an expansion of the log.

At first glance, it seems confusing that there are different definitions for the spread. Admittedly, one has to keep in mind that all formulae are only valid up to the order given in Eq. (167). Thus, although different, they are consistent and there is no fundamental reason to choose one definition of the spread over another.

One can also derive a general expression for the expectation value of the periodic position operator for computing the center of the localized function. Recall, that for a cubic simulation supercell the expectation value of the position operator is given as,

$$r_{\alpha ,n} = -{L\over 2\pi}{\rm Im}\;\log z_{\alpha,n}$$

$$z_{\alpha,n} = \int_{V}d{\bf r}\;\exp(i{\bf g}_{\alpha}\cdot {\bf r})\;\vert\psi_n({\bf r})\vert^2 ,$$ (170)

where $\hat{\bf g}_1=(1,0,0)$ , $\hat{\bf g}_2=(0,1,0)$ , and $\hat{\bf g}_3=(0,0,1)$ , and ${\rm Im}$ denotes the imaginary part. Again, the salient feature of Eq.(170) is that the expectation value of the exponential is invariant with respect to the choice of cell. Thus, a general equation for the expectation value of the position operator in supercells of arbitrary symmetry is,

$$r_{\alpha ,n}=-\sum_\beta{\stackrel{\leftrightarrow}{h}_{\alpha\beta}\over 2\pi} {\rm Im}\;\log z_{\alpha,n} .$$ (171)

We now proceed to determine the weights $\omega_I$ as defined in the sum rule Eq. (166) for supercells of general symmetry. Recall that the metric, $\stackrel{\leftrightarrow}{\bf g}$ will contain at most six independent entries as defined by the case of least symmetry, triclinic. Thus, Eq. (166) is a linear set of six equations with six unknowns. We have freedom to choose the six Miller indices, $\hat{{\bf g}}^I$ of which we are to take the linear combinations of. For computational convenience of computing $z_i$ we choose the first six indices that take you from one to the next point in the Brillouin zone. Namely, $\hat{{\bf g}}^1=(1,0,0)$ , $\hat{{\bf g}}^2=(0,1,0)$ , $\hat{{\bf g}}^3=(0,0,1)$ , $\hat{{\bf g}}^4=(1,1,0)$ , $\hat{{\bf g}}^5=(1,0,1)$ , $\hat{{\bf g}}^6=(0,1,1)$ . With this choice of $\hat{{\bf g}}^i$ the explicit system of equations based on Eq. (166) takes the following simple form,

$$\left(\begin{array}{cccccc} 1&0&0&1&1&0\ 0&0&0&1&0&0\ 0&0&0&0... ...ay}{c} g_{11}\ g_{12}\ g_{13}\ g_{22}\ g_{23}\ g_{33} \end{array}\right)$$ (172)

Thus, the solution to Eq. (172) yields the following set of general weights,

$$\omega_1 = g_{11}-g_{12}-g_{13}$$

$$\omega_2 = g_{22}-g_{12}-g_{23}$$

$$\omega_3 = g_{33}-g_{13}-g_{23}$$

$$\omega_4 = g_{12}$$

$$\omega_5 = g_{13}$$

$$\omega_6 = g_{23}$$ (173)

Eq. (173) indeed reduces to the specific cases computed in Ref.[91]. However, here, the case for triclinic symmetry is also included. Thus, with knowledge of the cell parameters, in conjunction with Eq. (167) allows one to compute the maximally localized WF.