Formulation of the problem

A system of noninteracting particles in a potential $V$ is described by

$$\hat{H}\psi_{n}(\mathbf{r}) = \left[ -\frac{\hbar^{2}}{2m}\... ...ght] \psi_{n}(\mathbf{r}) = \epsilon_{n} \psi_{n}(\mathbf{r}),$$ (1)

where $\hat{H}$ is the single-particle Hamiltonian of the system, with energy eigenvalues $\epsilon_{n}$ and corresponding eigenstates $\psi_{n}(\mathbf{r})$. The eigenstates satisfy the orthogonality constraints given by

$$\int \psi^{\ast}_{m}(\mathbf{r}) \psi^{\ }_{n}(\mathbf{r}) \mathrm{d} \mathbf{r} = \delta_{mn},$$ (2)

for all $m$ and $n$. For instance, within the Kohn-Sham scheme of density-functional theory [22,23,24], $\hat{H}$ is the Kohn-Sham Hamiltonian and $V$ is the effective potential.

The total band-structure energy is given by

$$E = \sum_{n} f^{\ }_{n} \epsilon^{\ }_{n} = \sum_{n} f^{\ }_... ...f{r}) \hat{H} \psi^{\ }_{n}(\mathbf{r}) \mathrm{d} \mathbf{r},$$ (3)

where $f_{n}$ is the occupancy of state $\psi_{n}(\mathbf{r})$ [25]: at the energy minimum, all states below and above the Fermi-level have occupancy unity and zero, respectively.

In the case of linear-scaling calculations, the $\mathcal{N}$ lowest extended eigenstates $\psi_{n}(\mathbf{r})$ ( $n \in \{1,\ldots,\mathcal{N} \}$) are expressed in terms of a set of $\mathcal{N}$ localized functions $\phi_{\alpha}(\mathbf{r})$ ( $\alpha \in \{1,\ldots,\mathcal{N} \}$) that are generally nonorthogonal:

$$\psi^{\ }_{n}(\mathbf{r}) = \sum_{\alpha} \phi^{\ }_{\alpha}(\mathbf{r}) M^{\alpha}_{\ n},$$ (4)

where $\mathbf{M}$ is a square ($\mathcal{N}$ by $\mathcal{N}$), nonsingular matrix of coefficients, and $\mathcal{N}$ can be equal to or greater than the number of occupied eigenstates. The overlap matrix $S_{\alpha\beta}$ of the localized functions $\phi_{\alpha}(\mathbf{r})$ is

$$S^{\ }_{\alpha\beta} = \int \phi^{\ast}_{\alpha}(\mathbf{r}) \phi^{\ }_{\beta}(\mathbf{r}) \mathrm{d} \mathbf{r},$$ (5)

and on substitution of Eq. (4) into the orthogonality relation given by Eq. (2) we find that $S_{\alpha\beta}$ satisfies

$$(M^{\dagger})_{n}^{\ \alpha} S^{\ }_{\alpha \beta} M^{\beta}_{\ m} = \delta^{\ }_{nm},$$ (6)

where a distinction has been made between contravariant and covariant quantities [26,27] through the use of superscript and subscript Greek suffixes, respectively.

Substituting Eq. (4) into the energy expression of Eq. (3), and defining

$$K^{\alpha\beta} = \sum_{n} M^{\alpha}_{\ n} f^{\ }_{n} (M^{\dagger})_{n}^{\ \beta},$$ (7)

$$H^{\ }_{\alpha\beta} = \int \phi^{\ast}_{\alpha}(\mathbf{r}) \hat{H} \phi^{\ }_{\beta}(\mathbf{r}) \mathrm{d} \mathbf{r},$$ (8)

the band-structure energy becomes

$$E = \sum_{\alpha\beta} H_{\alpha\beta} K^{\beta\alpha},$$ (9)

where $K^{\alpha\beta}$ is referred to as the density kernel [28].

We consider the localized functions $\phi_{\alpha}(\mathbf{r})$ to be represented in terms of a basis $D_{\mu}(\mathbf{r})$ as follows:

$$\phi^{\ }_{\alpha}(\mathbf{r}) = \sum_{\mu} D^{\ }_{\mu}(\mathbf{r}) c^{\mu}_{\ \alpha},$$ (10)

for some coefficients $c^{\mu}_{\ \alpha}$. As the basis functions $D_{\mu}(\mathbf{r})$ may be in general nonorthogonal, the tensor properties must be taken into account through the use of superscript and subscript Greek suffixes.

Defining

$$h^{\ }_{\mu\nu} = \int D^{\ast}_{\mu}(\mathbf{r}) \hat{H} D^{\ }_{\nu}(\mathbf{r}) \mathrm{d} \mathbf{r},$$ (11)

and using Eqs. (7)-(10), the energy may be written as

$$E = (c^{\dagger})_{\alpha}^{\ \mu} h^{\ }_{\mu\nu} c^{\nu}_... ...}^{\ \mu} h^{\ }_{\mu\nu} c^{\nu}_{\ \beta} M^{ \beta}_{\ n}.$$ (12)

Suffixes $\alpha$ and $\beta$ run over the localized functions $\{ \phi \}$, $\mu$ and $\nu$ run over the basis functions $\{ D \}$ and $n$ runs over the extended orthogonal orbitals $\{ \psi \}$. We have adopted the Einstein summation convention for all repeated Greek suffixes, and continue to do so from here on.

It is both convenient and physically meaningful to perform the minimization of the energy functional in two nested loops, as in the ensemble density-functional method of Marzari et al. [29]: in the inner loop we minimize the energy with respect to the elements of the density kernel $K^{\alpha\beta}$ using one of a number of methods [30,31,32] to impose the constraint that the ground state density matrix be idempotent and give the correct number of electrons; in the outer loop we optimize the localized functions $\phi_{\alpha}(\mathbf{r})$ with respect to their coefficients $c^{\mu}_{\ \alpha}$ in the basis $D_{\mu}(\mathbf{r})$ [16].