Hartree Potential and the Poisson Equation

Every self-consistent DFT calculation has to answer the same question at each SCF step: given the current electron density $n(\mathbf{r})$, what classical electrostatic potential does it create? That potential — the Hartree potential $v_\mathrm{H}(\mathbf{r})$ — is the dominant mean-field contribution to the effective potential felt by each electron. It is also the term where a plane-wave basis earns its keep most visibly: what is a three-dimensional integral equation in real space becomes a one-line pointwise division in reciprocal space.

All equations below use Rydberg atomic units ($\hbar = 1$, $m_e = 1/2$, $e^2 = 2$, $a_0 = 1$), matching the KRONOS source. The key unit-system consequences: bare Coulomb repulsion between two unit charges is $2/r$ (not $1/r$), and the Hartree potential prefactor is $8\pi$ (not $4\pi$). Where a Hartree-unit equivalent is shown for comparison, it is labeled explicitly.

The classical Coulomb integral

The Hartree potential is the electrostatic potential generated by the electron density treated as a classical charge distribution. Each volume element $d^3r'$ at position $\mathbf{r}'$ carries charge $-n(\mathbf{r}') d^3r'$ (electrons carry charge $-1$ in atomic units) and creates a Coulomb potential at $\mathbf{r}$. In Rydberg units the Coulomb kernel between two unit charges is $2/|\mathbf{r}-\mathbf{r}'|$, so the total Hartree potential is

$v_\mathrm{H}(\mathbf{r}) = \int \frac{2\, n(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} \, d^3r'$

In Hartree atomic units the same expression reads $\int n(\mathbf{r}')/|\mathbf{r}-\mathbf{r}'| \, d^3r'$ (Coulomb kernel is $1/r$). The factor of two in the Rydberg form comes directly from $e^2 = 2$ in those units; the physics is identical, only the prefactor differs.

$v_\mathrm{H}(\mathbf{r})$ is a purely classical, mean-field potential. It contains no exchange (no Pauli exclusion between electrons of the same spin) and no correlation (no correlated motion beyond the average field). Those effects are captured entirely by $v_\mathrm{xc}$. The Hartree term is the dominant electron-electron contribution in most systems — orders of magnitude larger than $v_\mathrm{xc}$ in typical valence electron densities — which makes solving the Poisson equation accurately at every SCF step non-negotiable.

The Poisson equation

The Hartree potential satisfies the Poisson equation. To see why, apply the Laplacian $\nabla^2$ to the convolution integral above. Using the distributional identity $\nabla^2(1/r) = -4\pi\,\delta^3(\mathbf{r})$ (the Green's function of $-\nabla^2$ is $1/4\pi r$), the differentiation collapses the convolution:

$\nabla^2 v_\mathrm{H}(\mathbf{r}) = 2 \int n(\mathbf{r}') \nabla^2 \frac{1}{|\mathbf{r}-\mathbf{r}'|} \, d^3r' = 2 \int n(\mathbf{r}') \bigl(-4\pi\bigr)\delta^3(\mathbf{r}-\mathbf{r}') \, d^3r'$

$\boxed{\nabla^2 v_\mathrm{H}(\mathbf{r}) = -8\pi\, n(\mathbf{r})} \quad \text{(Rydberg units)}$

The Hartree-unit version is $\nabla^2 v_\mathrm{H} = -4\pi n$ — the standard textbook form from Jackson's Classical Electrodynamics. The factor of two difference is again $e^2 = 2$ vs $e^2 = 1$.

Physically this is Gauss's law in differential form: $\nabla \cdot \mathbf{E} = 4\pi \rho_\mathrm{charge}$ (in Gaussian units), rewritten for the scalar potential $\mathbf{E} = -\nabla v_\mathrm{H}$. The electron density $n(\mathbf{r})$ plays the role of the charge source.

Diagonal in G-space

Solving the Poisson equation directly in real space for a periodic system would require inverting the full Laplacian — a dense $N_\mathrm{grid}^3 \times N_\mathrm{grid}^3$ operator. In G-space the Laplacian is diagonal: $\nabla^2 \to -|\mathbf{G}|^2$. Taking the Fourier transform of both sides of the Poisson equation,

$-|\mathbf{G}|^2 \, v_\mathrm{H}(\mathbf{G}) = -8\pi\, n(\mathbf{G})$

which rearranges instantly to

$\boxed{v_\mathrm{H}(\mathbf{G}) = \frac{8\pi\, n(\mathbf{G})}{|\mathbf{G}|^2}} \quad \text{(Rydberg units)}$

The Hartree-unit equivalent, $v_\mathrm{H}(\mathbf{G}) = 4\pi n(\mathbf{G})/|\mathbf{G}|^2$, is the form found in most electronic structure textbooks. In both cases the Poisson solve reduces to a single pointwise division over the G-vector list — $O(N_\mathrm{G})$ work with a tiny prefactor and no FFT required. This is why the Hartree potential update is never a bottleneck even for large supercells.

The density $n(\mathbf{G})$ used here lives on the full density FFT grid (controlled by ecutrho), not the smaller wavefunction grid (ecutwfc). This is essential: the density contains Fourier components up to $2\sqrt{E_\mathrm{cut}^\mathrm{wfc}}$ because it is formed from products of wavefunctions, so the Poisson solve must operate on the larger grid. In KRONOS, ecutrho $\geq 4 \times$ ecutwfc for norm-conserving pseudopotentials (PAW requires $\geq 12\times$).

The G=0 divergence and charge neutrality

The formula $v_\mathrm{H}(\mathbf{G}) = 8\pi n(\mathbf{G})/|\mathbf{G}|^2$ diverges at $\mathbf{G} = 0$. This is not a numerical accident — it reflects a genuine physical fact. For an infinite periodic solid, the $\mathbf{G}=0$ component of the Coulomb potential measures the energy cost of the net charge in the simulation cell interacting with its own infinite periodic images. That cost is infinite unless the cell is charge-neutral.

Physical resolution: charge neutrality. A real crystal is electrically neutral: the total electronic charge exactly cancels the ionic charge,

$N_\mathrm{electrons} + N_\mathrm{ions}^\mathrm{eff} = 0 \quad \Longleftrightarrow \quad n(\mathbf{G}=0) = \frac{N_e}{\Omega} = \bar{n}$

where $N_e$ is the number of valence electrons, $\Omega$ is the unit cell volume, and $\bar{n}$ is the mean density. KRONOS enforces this as a hard requirement: if the input does not produce a charge-neutral cell (electrons plus pseudopotential valence charges), the calculation aborts.

Numerical convention: set $v_\mathrm{H}(\mathbf{G}=0) = 0$. When the cell is neutral, the infinite $1/|\mathbf{G}|^2$ divergence at $\mathbf{G}=0$ from the Hartree term cancels exactly against the identical divergence from the ionic local potential $v_\mathrm{loc}(\mathbf{G}=0)$ and the Ewald reciprocal-space sum. These three terms together produce a finite, well-defined total energy. In KRONOS — as in Quantum ESPRESSO and most plane-wave codes — all three contributions are set to zero at $\mathbf{G}=0$ individually, and their combined finite remainder is carried through the Ewald correction. Setting $v_\mathrm{H}(\mathbf{G}=0) = 0$ is equivalent to choosing the arbitrary zero of electrostatic energy, which has no observable consequence for energy differences or forces.

In hartree.cpp, the implementation checks g2 < 1e-12 and assigns zero explicitly, bypassing the division:

if (g2 < 1.0e-12) {
    vhartree[ig] = complex_t{0.0, 0.0};
} else {
    vhartree[ig] = prefactor * density_g[ig] / g2;
}

This threshold-based skip also prevents any finite-precision issues from near-zero $|\mathbf{G}|^2$ values near the $\Gamma$ point.

Hartree energy

The Hartree energy is the classical electrostatic self-energy of the electron charge distribution. It is defined as half the interaction energy of the density with its own Hartree potential (the factor of $1/2$ avoids double-counting each pair of interacting volume elements):

$E_\mathrm{H} = \frac{1}{2} \int v_\mathrm{H}(\mathbf{r}) \, n(\mathbf{r}) \, d^3r$

Using Parseval's theorem to convert to G-space (the $\Omega$ factor comes from the cell volume normalization of the Fourier transform):

$E_\mathrm{H} = \frac{\Omega}{2} \sum_{\mathbf{G} \neq 0} v_\mathrm{H}^*(\mathbf{G}) \, n(\mathbf{G})$

The $\mathbf{G}=0$ term is excluded from the sum because $v_\mathrm{H}(\mathbf{G}=0) = 0$ by convention. Substituting $v_\mathrm{H}(\mathbf{G}) = 8\pi n(\mathbf{G})/|\mathbf{G}|^2$, this can also be written

$E_\mathrm{H} = \frac{\Omega}{2} \sum_{\mathbf{G} \neq 0} \frac{8\pi\, |n(\mathbf{G})|^2}{|\mathbf{G}|^2}$

which is explicitly real and non-negative — the Hartree energy is always positive, reflecting the repulsive character of electron-electron Coulomb interaction.

FFT normalization note. In the KRONOS FFT convention both density_g and vhartree_g carry a factor of $N_\mathrm{grid}$ relative to the physics convention (FFTW3 does not normalize its transforms). The product $v_\mathrm{H}^*(\mathbf{G}) \, n(\mathbf{G})$ is therefore $N_\mathrm{grid}^2$ too large, and HartreeSolver::energy() divides by num_grid * num_grid to recover the correct Rydberg value:

const double n2 = static_cast<double>(num_grid) * static_cast<double>(num_grid);
return 0.5 * volume * esum / n2;

Double-counting in the DFT total energy

The Kohn-Sham band energy $E_\mathrm{band} = \sum_{n\mathbf{k}} f_{n\mathbf{k}} \varepsilon_{n\mathbf{k}}$ is computed by summing the KS eigenvalues. Each eigenvalue $\varepsilon_{n\mathbf{k}}$ is the expectation value of the full effective-potential Hamiltonian — kinetic plus $v_\mathrm{H}$ plus $v_\mathrm{xc}$ plus external. The band energy sum therefore already includes the Hartree potential contribution $\int v_\mathrm{H}(\mathbf{r}) \, n(\mathbf{r}) \, d\mathbf{r}$, which is twice $E_\mathrm{H}$ (since $E_\mathrm{H} = \frac{1}{2}\int v_\mathrm{H} n$). Similarly, $\int v_\mathrm{xc} \, n$ is already inside the eigenvalue sum but is not $E_\mathrm{xc}$ (the latter is a different functional of $n$).

To recover the correct total energy from the band sum, both overcounts must be removed. In Rydberg units:

$\boxed{E_\mathrm{tot} = E_\mathrm{band} - E_\mathrm{H}[n] + E_\mathrm{xc}[n] - \int v_\mathrm{xc}(\mathbf{r}) \, n(\mathbf{r}) \, d\mathbf{r} + E_\mathrm{ion-ion}}$

Breaking down each term:

• $E_\mathrm{band}$: sum of KS eigenvalues, includes kinetic + $\int v_\mathrm{H} n$ + $\int v_\mathrm{xc} n$ + external.
• $-E_\mathrm{H}$: removes the overcounted half of the Hartree self-energy ($E_\mathrm{band}$ has $\int v_\mathrm{H} n = 2E_\mathrm{H}$, so subtracting $E_\mathrm{H}$ leaves $+E_\mathrm{H}$, the correct contribution).
• $+E_\mathrm{xc}[n]$: the true XC energy, a functional of the density.
• $-\int v_\mathrm{xc} n$: removes the XC potential energy that was included in the eigenvalues.
• $E_\mathrm{ion-ion}$: classical ion-ion interaction, computed by Ewald summation. See Ewald Summation in the legacy Physics Notes for the derivation; a dedicated page is on the roadmap.

This decomposition is implemented in src/solver/scf.cpp. The Hartree energy is computed after each SCF iteration by HartreeSolver::energy(), and the double-counting correction is applied immediately before writing the total energy to the JSON output.

How KRONOS implements this

The Hartree solver lives in src/potential/hartree.hpp and src/potential/hartree.cpp. The public interface exposes two methods:

• HartreeSolver::compute(density_g) — takes the electron density in G-space (indexed on the full density grid, i.e. ecutrho G-vectors), applies $v_\mathrm{H}(\mathbf{G}) = 8\pi n(\mathbf{G})/|\mathbf{G}|^2$ with the G=0 guard, and returns $v_\mathrm{H}(\mathbf{G})$ as a CVec.

• HartreeSolver::energy(density_g, vhartree_g, volume, num_grid) — computes $E_\mathrm{H} = (\Omega/2)\sum_{\mathbf{G}\neq 0} v_\mathrm{H}^*(\mathbf{G})\,n(\mathbf{G})$ with the $N_\mathrm{grid}^2$ normalization correction described above. The energy loop starts at ig = 0 but the G=0 term contributes nothing because vhartree_g[0] = 0.

The returned $v_\mathrm{H}(\mathbf{G})$ is added to the local pseudopotential and XC contributions to form the total effective potential $v_\mathrm{eff}(\mathbf{G})$, which is then inverse-FFTed to real space for the local Hamiltonian application. Because the Hartree solve is a single pass through $N_G \sim \mathcal{O}(10^5)$ G-vectors with no FFT, it is typically less than 0.1% of the total wall time even for large cells.

References

• Martin, R. M. Electronic Structure: Basic Theory and Practical Methods, Cambridge University Press, Ch. 7 (2004) — Hartree energy, double-counting, and the total energy expression.
• Pickett, W. E. "Pseudopotential methods in condensed matter applications", Comput. Phys. Rep. 9, 115 (1989) — comprehensive treatment of the Poisson solve and G=0 convention in plane-wave DFT.
• Jackson, J. D. Classical Electrodynamics, 3rd ed., Wiley, Ch. 1 (1999) — Poisson equation derivation from Gauss's law and the Green's function of the Laplacian.
• Payne, M. C., Teter, M. P., Allan, D. C., Arias, T. A., Joannopoulos, J. D. "Iterative minimization techniques for ab initio total-energy calculations", Rev. Mod. Phys. 64, 1045 (1992) — practical plane-wave implementation of Hartree and double-counting.