quantum chemistry - How does one actually get the energy from the Kohn–Sham equations?

This might be a silly question, but how does one acquire the energy of the system from the Kohn–Sham equations?

$\left[ -\frac{1}{2}\nabla^2 + V_{eN}(\vec{r}) + V_{ee}(\vec{r}) + V_{xc}(\vec{r}) \right]\phi_i(\vec{r})=\epsilon_i\phi_i(\vec{r})$

Where $\nabla$ is the Laplacian operator, $V_{eN}(\vec{r})$ is the electron-nuclei attraction, $V_{ee}(\vec{r})$ is the electron-electron repulsion, $V_{xc}(\vec{r})$ is the exchange-correlation potential, $\phi_i(\vec{r})$ are the one-electron wave functions and $\epsilon_i$ are the orbital energies.

Assuming that we run through this set of equations for each electron, we get a set of orbital energies: $\epsilon_i \cdots \epsilon_N$. How is the overall total energy of the system then calculated?

Answer

A common mistake is the thought that the total energy is the sum of all orbital energies $\{\epsilon_i\}$.

From Step #6 of Daniel Crawford's SCF programming project (modified slightly in some places):

The SCF electronic energy may be computed using the density matrix as:

$$E_{\text{elec}} = \sum_{\mu\nu}^{\text{AO}} D_{\mu\nu} (H_{\mu\nu}^{\text{core}} + F_{\mu\nu})$$

The total energy is the sum of the electronic energy and the nuclear repulsion energy:

$$E_{\text{total}} = E_{\text{elec}} + E_{\text{nuc}},$$

where the density matrix is defined as (Step #8)

$$D_{\mu\nu} = \sum_{m}^{\text{occ. MO}} C_{\mu m} C_{\nu m},$$

the Fock matrix as (Step #7)

$$\begin{align} F_{\mu\nu} &= H_{\mu\nu}^{\text{core}} + \sum_{\lambda\sigma}^{\text{AO}} D_{\lambda\sigma} \left[ 2(\mu\nu|\lambda\sigma) - (\mu\lambda|\nu\sigma) \right] \\ &= H_{\mu\nu}^{\text{core}} + 2 J_{\mu\nu} - K_{\mu\nu} , \end{align}$$

and the core Hamiltonian as (Step #2)

$$H_{\mu\nu}^{\text{core}} = T_{\mu\nu} + V_{\mu\nu}.$$

I've also introduced the definitions of the Coulomb matrix $J$ and the exchange matrix $K$:

$$\begin{align} J_{\mu\nu} &= \sum_{\lambda\sigma}^{\text{AO}} D_{\lambda\sigma} (\mu\nu|\lambda\sigma) \\ K_{\mu\nu} &= \sum_{\lambda\sigma}^{\text{AO}} D_{\lambda\sigma} (\mu\lambda|\nu\sigma) \\ \end{align}$$

---

Now, identify each of the terms in the Kohn-Sham equations with the terms from above.

$$\begin{align} \hat{T}_{e} &= -\frac{1}{2} \nabla^2 \rightarrow T_{\mu\nu} = \left< \chi_{\mu} \left| \hat{T} \right| \chi_{\nu} \right> \\ \hat{V}_{eN}(\vec{r}) &= \sum_{A}^{\text{nuclei}} \frac{Z_A}{|\vec{r} - \vec{R}_{A}|} \rightarrow V_{\mu\nu} = \left< \chi_{\mu} \left| \hat{V}_{eN} \right| \chi_{\nu} \right> \\ \hat{V}_{ee}(\vec{r}) &\stackrel{?}{\rightarrow} 2 \hat{J} \\ \hat{V}_{\text{XC}}(\vec{r}) &\stackrel{?}{\rightarrow} - \hat{K} \end{align}$$

This last part isn't quite correct though. Usually, when looking at the Kohn-Sham equations, one replaces the full electron-electron interaction $\hat{V}_{ee}$ with the sum of the Hartree potential $\hat{V}_{H}$, which gives the Coulomb energy, and the exchange-correlation potential $\hat{V}_{\text{XC}}$, which replaces the exact exchange $\hat{K}$ with a (currently approximate) expression for both the exchange term and the true electron-electron (correlated) interaction.

In terms of how the energy is actually calculated, all quantities from above are the same as in Hartree-Fock theory, except the calculation of the exact exchange integrals during the Fock build is replaced with calculating the exchange-correlation matrix $F^{\text{XC}}$, leading to

$$\begin{align} F_{\mu\nu}^{\alpha} &= H_{\mu\nu}^{\text{core}} + J_{\mu\nu} + F_{\mu\nu}^{\text{XC}\alpha} \\ F_{\mu\nu}^{\beta} &= H_{\mu\nu}^{\text{core}} + J_{\mu\nu} + F_{\mu\nu}^{\text{XC}\beta} \end{align}$$

For a density functional approximation (DFA) based on the generalized gradient approximation (GGA), where the functional is dependent on both the density $\rho(\mathbf{r})$ and its gradient $\nabla \rho(\mathbf{r})$,

$$\begin{align} E_{\text{XC}} &= \int f_{GGA}^{\text{DFA}}(\rho_{\alpha},\rho_{\beta},\gamma_{\alpha\alpha},\gamma_{\alpha\beta},\gamma_{\beta\beta}) \, \mathrm{d}\mathbf{r} \\ \gamma_{\alpha\alpha} &= |\nabla \rho_{\alpha}|^{2} \\ \gamma_{\beta\beta} &= |\nabla \rho_{\beta}|^{2} \\ \gamma_{\alpha\beta} &= \nabla \rho_{\alpha} \cdot \nabla \rho_{\beta} \\ \end{align}$$

The exchange-correlation parts of the Fock matrices are given by

$$F_{\mu\nu}^{\text{XC}\alpha} = \int \left[ \frac{\partial f}{\partial \rho_{\alpha}} \chi_{\mu}\chi_{\nu} + \left( 2\frac{\partial f}{\partial \gamma_{\alpha\alpha}} \nabla\rho_{\alpha} + \frac{\partial f}{\partial \gamma_{\alpha\beta}} \nabla\rho_{\beta} \right) \cdot \nabla(\chi_{\mu}\chi_{\nu}) \right] \mathrm{d}\mathbf{r}$$

$f^{\text{DFA}}$, $\frac{\partial f^{\text{DFA}}}{\partial \rho}$, and $\frac{f^{\text{DFA}}}{\partial \gamma}$ are unique closed-form expressions for each DFA, and are usually evaluated numerically on an atom-centered grid (ACG) such as a Lebedev grid. This generally requires mapping the set of AOs/basis functions $\{\chi\}$ onto this grid.

References

$\tiny{\text{As usual, sorry if I'm lazy with notation, being consistent is so difficult...}}$