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...}}$
