How does one actually get the energy from the Kohn–Sham equations?
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} &amp;= H_{\mu\nu}^{\text{core}} + \sum_{\lambda\sigma}^{\text{AO}} D_{\lambda\sigma} \left[ 2(\mu\nu|\lambda\sigma) - (\mu\lambda|\nu\sigma) \right] \\
&amp;= 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} &amp;= \sum_{\lambda\sigma}^{\text{AO}} D_{\lambda\sigma} (\mu\nu|\lambda\sigma) \\
K_{\mu\nu} &amp;= \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} &amp;= -\frac{1}{2} \nabla^2 \rightarrow T_{\mu\nu} = \left&lt; \chi_{\mu} \left| \hat{T} \right| \chi_{\nu} \right&gt; \\
\hat{V}_{eN}(\vec{r}) &amp;= \sum_{A}^{\text{nuclei}} \frac{Z_A}{|\vec{r} - \vec{R}_{A}|} \rightarrow V_{\mu\nu} = \left&lt; \chi_{\mu} \left| \hat{V}_{eN} \right| \chi_{\nu} \right&gt; \\
\hat{V}_{ee}(\vec{r}) &amp;\stackrel{?}{\rightarrow} 2 \hat{J} \\
\hat{V}_{\text{XC}}(\vec{r}) &amp;\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} &amp;= H_{\mu\nu}^{\text{core}} + J_{\mu\nu} + F_{\mu\nu}^{\text{XC}\alpha} \\
F_{\mu\nu}^{\beta} &amp;= 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}} &amp;= \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} &amp;= |\nabla \rho_{\alpha}|^{2} \\
\gamma_{\beta\beta} &amp;= |\nabla \rho_{\beta}|^{2} \\
\gamma_{\alpha\beta} &amp;= \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


John A. Pople; Peter M. W. Gill; Benny G. Johnson; Kohn–Sham density-functional theory within a finite basis set. Chem. Phys. Lett. 1992, 199 557-560


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