Background

At this stage, the only differences between your Schrödinger solver and a density functional solver are the energy and gradient functions, getE(), getgrad() and getH(), respectively.

To see this, we begin by noting that one important difference between the energy functions is that, now, each wave function $\psi_i$ has an associated occupancy, or ``filling factor,'' $f_i$. For simplicity we shall here always work in cases where all states have the same occupancy $f_i=f$. (Many physical systems have this property; usually $f=2$. The generalization to difference occupancies for different orbitals is conceptually straight forward but involves a lot of messy algebra.)

Under the conditions of the previous paragraph, three of the quantities from density-functional theory have very similar forms to those from the Schrödinger case. For instance, the density-functional form for the kinetic energy $(-1/2) \sum_i f_i \int \psi_i^* \nabla^2 \psi_i$ simply picks up an extra factor of $f$,

$$KE=-\frac{1}{2} \sum_i f_i \int \psi_i^* \nabla^2 \psi_i = f \cdot \frac{1}{2} {\mbox{Tr\,}}(W^\dagger L W N^{-1/2}).$$ (1)

The total electron density $\sum_i f_i \vert\psi_i\vert^2$ picks up a similar factor of $f$,

$$\vec n=f \cdot {\mbox{diag\,}}\left[ {\cal I}W N^{-1} \left({\cal I}W\right)^2 \right].$$ (2)

Finally, in terms of the vector of sample values $\vec n$, the formula for the electron-nuclear potential energy $PE=\int V(r) n(r)$ is identical to the corresponding formula in the Schrödinger case and thus leads to the same result,

$$PE=\tilde V^\dagger \vec n,$$ (3)

where

$$\tilde V \equiv {\cal J}^\dagger {\cal O}{\cal J}\vec V,$$ (4)

with $\vec V$ being the sample values of the potential.

The first new term appearing in density-functional theory is the electron-electron potential energy $U_{ee} = \frac{1}{2} \int n(r)^* \phi(r)$, where $\nabla^2 \phi(r) = - 4 \pi n(r)$. (As in a similar derivation in lecture, the ``*'' on $n(r)$ doesn't change the integral because $n(r)$ is real. We include it, however, so that we may use our currently defined operator ${\cal O}$.) Inserting the expansions $n(r) = \sum_\alpha b_\alpha(r) \hat n_\alpha$ and $\phi(r) = \sum_\alpha b_\alpha(r) \hat \phi_\alpha$ into the integral for $U_{ee}$, using our solution for Poisson's equation from Problem Set 7 and rearranging quantities within $\dagger$'s, we find

$$U_{ee} = \frac{1}{2} \vec n^\dagger {\cal J}^\dagger {\cal O}\hat \phi,$$ (5)

where the expansion coefficients for the solution to Poisson's equation are

$$\hat \phi=-4 \pi L^{-1} {\cal O}{\cal J}\vec n.$$

The final term we require is the exchange correlation energy $\int n_(r)^* \epsilon_{xc}(n(r))$. Defining the operator $\epsilon_{xc}(\vec n)$ as returning a vector each of whose components is the evaluation of the exchange-correlation energy for the corresponding component of $\vec n$ and inserting appropriate expansions gives the result,

$$E_{xc} = \vec n^\dagger {\cal J}^\dagger {\cal O}{\cal J}\epsilon_{xc}({\cal I}\vec n).$$ (6)

In summary, in terms of unconstrained wave function coefficients $W$, the total energy within density functional theory is

$$E_{LDA} = -f \cdot \frac{1}{2}{\mbox{Tr\,}}\left(W^\dagger L W N^{-1}\right) + \tilde V^\dagger n$$ (7)

$$+ \frac{1}{2} \vec n^\dagger {\cal J}^\dagger {\cal O}\hat \phi + \vec n^\dagger {\cal J}^\dagger {\cal O}{\cal J}\epsilon_{xc}(\vec n),$$

where

$$\begin{eqnarray*} N & \equiv & W^\dagger {\cal O}W \\ \vec n & \equiv & f \cdot... ...] \\ \hat \phi & \equiv & -4 \pi L^{-1} {\cal O}{\cal J}\vec n. \end{eqnarray*}$$

As with the Schrödinger case, two effects contribute to the gradient. First, there are contributions due to the constraints, which turn out to be of an identical form as we had previously. Including the appropriate factors for $f$, we thus have

$$\nabla_W E = f \left( H W - {\cal O}W N^{-1} W^\dagger H W \right) N^{-1}.$$ (8)

Finally, the operator $H()$ encapsulates the remaining contributions to the gradient, which come from the basic structure of the energy and has the form,

$$H \equiv -\frac{1}{2} L + {\cal I}^\dagger \left({\mbox{Diag\,}} V_{eff}\right) {\cal I},$$ (9)

where

$$V_{eff} \equiv \tilde V + {\cal J}^\dagger {\cal O}\hat \phi... ..._{xc}'(\vec n)\right] {\cal J}^\dagger {\cal O}{\cal J}\vec n,$$

with the operator $\epsilon_{xc}'(\vec n)$ defined similarly to $\epsilon_{xc}(\vec n)$ but now evaluating $\partial \epsilon_{xc}/\partial n$ on each component of $\vec n$.