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?

  • 5
    $\begingroup$ Is your question how you write the total energy in equations? Because then I think wikipedia is a fine start. Or is it how to get the total energy in practical? Like inside the programs. $\endgroup$ Commented Aug 27, 2017 at 22:57

1 Answer 1


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.


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.