Total energy

The electronic energy is decomposed as

\[ E[n] = T_s[n] + E_{ext}[n] + E_H[n] + E_{xc}[n], \]

with the external (electron-nucleus) energy \(E_{ext} = \int n(\mathbf{r})\, v_{ext}(\mathbf{r})\, d\mathbf{r}\).

Kinetic energy density functionals (OF-DFT)

  • Thomas-Fermi: \(T_{TF}[n] = C_{TF} \int n^{5/3}\, d\mathbf{r}\) with \(C_{TF} = \tfrac{3}{10}(3\pi^2)^{2/3}\).

  • von Weizsäcker: \(T_{vW}[n] = \tfrac{1}{8}\int \frac{|\nabla n|^2}{n}\, d\mathbf{r}\), exact for one- and two-electron systems.

  • TFvW: \(T[n] = c_{TF}\,T_{TF}[n] + \lambda\, T_{vW}[n]\). The gradient prefactor \(\lambda\) interpolates between the Thomas-Fermi limit (\(\lambda = 0\)) and the full von Weizsäcker correction.

The OF-DFT engine minimizes \(E[n]\) directly over the positivity-preserving variable \(w = \sqrt{n}\), constrained to the sphere \(\int w^2\, d\mathbf{r} = N\). The constrained gradient is \(g = 2w\,(v_{\mathrm eff} - \mu)\) with the chemical potential \(\mu = \tfrac{1}{N}\int w^2 v_{\mathrm eff}\, d\mathbf{r}\), which makes \(g\) tangent to the constraint sphere (\(\langle g, w\rangle = 0\)). Search directions are built with the Polak-Ribière(+) conjugate-gradient formula,

\[ d_k = -g_k + \beta_k\, d_{k-1}, \qquad \beta_k = \max\!\left(0,\ \frac{\langle g_k,\, g_k - g_{k-1}\rangle} {\langle g_{k-1},\, g_{k-1}\rangle}\right), \]

projected back onto the tangent space at \(w\) and followed by a backtracking line search with the normalization retraction \(w' = \mathrm{normalize}(w + \alpha d)\). The method restarts to steepest descent periodically and whenever a conjugate direction stops descending, guaranteeing a monotonically decreasing energy.

Exchange-correlation

The shared LDA functional combines Dirac (Slater) exchange, \(E_x = -C_x \int n^{4/3}\, d\mathbf{r}\) with \(C_x = \tfrac{3}{4}(3/\pi)^{1/3}\), and Perdew-Wang (1992) correlation. The same interface is reused by the OF-DFT and KS-DFT engines.

Kohn-Sham DFT

The Kohn-Sham Hamiltonian is

\[ \hat{H} = -\tfrac{1}{2}\nabla^2 + v_{ext} + v_H[n] + v_{xc}[n]. \]

The kinetic operator is applied spectrally (FFT), \(-\tfrac{1}{2}\nabla^2\psi = \tfrac{1}{2}\,\mathcal{F}^{-1}\!\big(|\mathbf{G}|^2\,\mathcal{F}\psi\big)\), and the occupied orbitals are obtained with a matrix-free Hermitian eigensolver inside a linearly-mixed fixed-point SCF loop. The total energy uses the band-energy decomposition

\[ E = \sum_i f_i \varepsilon_i - E_H[n] + E_{xc}[n] - \int v_{xc}[n]\, n\, d\mathbf{r}. \]

Frozen-Density Embedding

The supersystem density is partitioned as \(n_{tot} = \sum_I n_I\). Each subsystem is relaxed in the embedding potential generated by the others (Wesolowski-Warshel),

\[ v^I_{emb} = v^{others}_{ext} + v_H[n_{others}] + \big(v_{xc}[n_{tot}] - v_{xc}[n_I]\big) + \big(v_{Ts}[n_{tot}] - v_{Ts}[n_I]\big), \]

where the nonadditive kinetic potential uses an explicit KEDF. Subsystems may be solved with the OF-DFT or KS-DFT engine, enabling OF-in-OF, KS-in-KS, and mixed embedding. Freeze-and-thaw cycles alternate which subsystem is relaxed until the total density and embedded energy are self-consistent.