Introduction

Several methods have been proposed (Bekas 2005)

Theory

We define a set of orthonormal functions $$\left\{\varsigma(\vec{r})_j\right\}_j$$ that span the same basis as the orbitals $$\left\{\varphi(\vec{r})_j\right\}_j$$. The are related through a orthogonal transformation $${\overleftrightarrow{R}}$$ such that $\varphi_i(\vec{r}) = \sum_jR_{ij}\varsigma_j(\vec{r}) \ .$ \begin{aligned} n(\vec{r}) &= \sum_if_i\varphi_i^*(\vec{r})\varphi_i(\vec{r}) \\ &=\sum_{ijk}f_iR^*_{ij}\varsigma_j^*(\vec{r})R_{ik}\varsigma_k(\vec{r})\\ &=\sum_{jk}\varsigma_j^*(\vec{r})\left(\sum_iR^*_{ij}f_iR_{ik}\right)\varsigma_k(\vec{r})\\ &=\sum_{jk}\varsigma_j^*(\vec{r})F_{jk}\varsigma_k(\vec{r})\end{aligned}

$F_{jk}=\sum_iR^*_{ij}f_iR_{ik} \label{eq:occupationmatrix}$

\begin{aligned} \mathrm{Tr}({\overleftrightarrow{F}})&=\sum_jF_{jj}\\ &=\sum_{ij}R^*_{ij}f_iR_{ij}\\ &=\sum_{j}f_j\\ &=N\end{aligned}

\begin{aligned} \sum_if_i\epsilon_i &= \sum_if_i\int\mathrm{d}\vec{r}\varphi_i^*(\vec{r})\int\mathrm{d}\vec{r}'\mathrm{H}(\vec{r},\vec{r}')\varphi_i(\vec{r}') \\ &=\sum_{ijk}f_i\int\mathrm{d}\vec{r}R^*_{ij}\varsigma_j^*(\vec{r})\int\mathrm{d}\vec{r}'\mathrm{H}(\vec{r},\vec{r}')R_{ik}\varsigma_k(\vec{r}')\\ &=\sum_{jk}F_{jk}\int\mathrm{d}\vec{r}\varsigma_j^*(\vec{r})\int\mathrm{d}\vec{r}'\mathrm{H}(\vec{r},\vec{r}')\varsigma_k(\vec{r}')\\ &=\sum_{jk}F_{jk}{\cal{E}}_{jk}\\ &=\mathrm{Tr}({\overleftrightarrow{F}}^\intercal {\overleftrightarrow{{\cal{E}}}})\end{aligned}

${\cal{E}}_{jk}=\int\mathrm{d}\vec{r}\mathrm{d}\vec{r}'\varsigma_j^*(\vec{r})\mathrm{H}(\vec{r},\vec{r}')\varsigma_k(\vec{r}')$

In temperature-dependent DFT, the occupations are obtained from the Fermi-Dirac distribution $f_i=\frac1{1+e^{(\epsilon_i-\mu)/kT}} \label{eq:fermidirac}$ where the chemical potential $$\mu$$ is obtained from imposing the condition that the occupations must sum to the total number of electrons. This procedure establishes a relation between the set of KS eigenvalues $$\epsilon_i$$ and the occupations $$f_i$$. For our approach, we must find a relation between the subspace Hamiltonian $${\overleftrightarrow{h}}$$ and the occupation matrix $${\overleftrightarrow{F}}$$. Formally, we can diagonalize $${\overleftrightarrow{h}}$$ to obtain both $${\overleftrightarrow{R}}$$ and $$\epsilon_i$$, from there we can then obtain $${\overleftrightarrow{F}}$$ from eqs. \ref{eq:fermidirac} and \ref{eq:occupationmatrix} as $F_{jk}=\sum_iR^*_{ij}\frac1{1+e^{(\epsilon_i-\mu)/kT}}R_{ik}\ . \label{eq:fermidiracdiagonal}$ Of course this involves the diagonalization of $${\overleftrightarrow{{\cal{E}}}}$$ which is what we wanted to avoid from the beginning. But since $${\overleftrightarrow{R}}$$ and $$\epsilon_i$$ are the eigenvectors and eigenvalues of $${\overleftrightarrow{{\cal{E}}}}$$ we can rewrite eq. \ref{eq:fermidiracdiagonal} as a matrix formula ${\overleftrightarrow{F}}=\left[I + \exp\left(\frac1{kT}\left({\overleftrightarrow{{\cal{E}}}}-\mu\right)\right)\right]^{-1} \ . \label{eq:fermidiracmatrix}$ with the condition that the value of $$\mu$$ must be such that $$\mathrm{Tr}({\overleftrightarrow{F}})=N$$.

We have converted the problem from the diagonalization of $${\cal{E}}$$ to the evaluation of the function of a matrix for which several methods exists. So the question is now, can we find a method to evaluate \ref{eq:fermidiracmatrix} that is more efficient for large matrices than diagonalization.

While most approaches are focused on sparse matrices (Sidje 2010)

We can use a Padé approximation (Moler 2003) $\exp\left({\overleftrightarrow{A}}\right)\sim{\overleftrightarrow{R}}_n({\overleftrightarrow{A}})=\left[{\overleftrightarrow{P}}_n({\overleftrightarrow{A}})\right]^{-1}\left[{\overleftrightarrow{Q}}_n({\overleftrightarrow{A}})\right]$ where $P_n({\overleftrightarrow{A}})=\sum_{j=0}^n\frac{(2n-j)!\,n!}{(2n)!\,j!\,(n-j)!}{\overleftrightarrow{A}}^j$ and $Q_n({\overleftrightarrow{A}})=\sum_{j=0}^n\left(-1\right)^j\frac{(2n-j)!\,n!}{(2n)!\,j!\,(n-j)!}{\overleftrightarrow{A}}^j\ .$

$\frac{\partial{\overleftrightarrow{F}}}{\partial \mu} = \frac1T\exp\left(\frac1{T}\left({\overleftrightarrow{{\cal{E}}}}-\mu\right)\right){\overleftrightarrow{F}}^{-2}$

(CHECK ORDER, NON-COMMUTATIVE MULTIPLICATION)