vibe-qc/Documentation

Molecular Kohn-Sham DFT on a fixed quadrature

For real molecular RKS and collinear UKS, the XC potential stored in the AO Fock matrix represents the derivative of the discretized XC energy with respect to the density matrix. The same finite-sum argument covers LDA and GGA. A GGA also differentiates the density gradients. In UKS this introduces a mixed-spin term whose coefficient differs from the same-spin coefficient. The PBE energy depends on both spin densities and their gradients, as in equations (2) and (3) of Perdew, Burke and Ernzerhof (1996).

Densities, gradients and the XC energy

Fix real AO values \(\chi_\mu(g)\), their three spatial derivatives and weights \(w_g\) at finitely many grid points. RKS uses the spin-summed density \(D\). UKS uses \(D_\alpha,D_\beta\) separately, with total density \(D_\alpha+D_\beta\). For each density matrix supplied to a channel, define

\[ \rho_g = \sum_{\mu\nu}D_{\mu\nu}\chi_\mu\chi_\nu,\qquad \partial_k\rho_g=\sum_{\mu\nu}D_{\mu\nu} [(\partial_k\chi_\mu)\chi_\nu+\chi_\mu(\partial_k\chi_\nu)] . \]

All AO values on the right are at \(g\). The sums use ordered AO pairs. For a symmetric density, interchanging the two indices shows that the two gradient terms are equal after summation. This gives the factor two in the production density-gradient contraction. Both fields vary linearly along a density perturbation \(D+tH\). Lean proves these identities in ks_density_gradient and ks_field_line.

Write \(e\) for energy per unit volume, including the total density factor when a functional supplies energy per electron. The molecular functional wrapper returns this volume integrand as exc. Both drivers compute \(E_{\rm xc}=\sum_g w_g e_g\). For RKS PBE the local inputs are \(\rho\) and \(\sigma=\nabla\rho\cdot\nabla\rho\). For UKS PBE they are \(\rho_\alpha\), \(\rho_\beta\), \(\sigma_{\alpha\alpha}\), \(\sigma_{\alpha\beta}\) and \(\sigma_{\beta\beta}\), where \(\sigma_{st}=\nabla\rho_s\cdot\nabla\rho_t\).

Why the projected AO potential is the derivative

Let \(a=\nabla\rho_\alpha\) and \(b=\nabla\rho_\beta\). For density perturbations that induce \(da,db\), the mixed invariant has derivative \(da\cdot b+a\cdot db\). A same-spin invariant has derivative \(2a\cdot da\). The Lean theorem ks_sigma_deriv differentiates the finite dot product along the affine perturbation. If \(v_{st}\) is the partial derivative of \(e\) with respect to its single input \(\sigma_{st}\), collecting the UKS terms gives

\[ u_\alpha=2v_{\alpha\alpha}a+v_{\alpha\beta}b,\qquad u_\beta=2v_{\beta\beta}b+v_{\alpha\beta}a . \]

For RKS there is one vector \(u=2v_\sigma\nabla\rho\). These coefficient identities are ks_uks_gradient_flow and ks_rks_gradient_flow. In LDA the gradient terms drop out. For each spin channel \(s\), or the single total-density channel of RKS, the energy derivative is represented by

\[ V^s_{\mu\nu}=\sum_g w_g\left[ v_{\rho_s}\chi_\mu\chi_\nu +u_s\cdot\{(\nabla\chi_\mu)\chi_\nu+\chi_\mu\nabla\chi_\nu\}\right]. \]

To prove the projection once for arbitrary finite grids, flatten the ordered AO pair and spin into an index \(i\). Write \(K_{gfi}\) for the AO kernel of field \(f\) (density or a Cartesian gradient component of a named channel). The field variation is \(\sum_i K_{gfi}H_i\). At the actual local fields, assume the volume integrand is differentiable and that the supplied coefficients \(v_{gf}\) represent its derivative in these coordinates. Applying the chain rule and interchanging finite sums gives

\[ \left.\frac{d}{dt}E_{\rm xc}(D+tH)\right|_{t=0} =\sum_i H_i \sum_g w_g\sum_f v_{gf}K_{gfi}. \]

ks_grid_adjoint proves the finite sum interchange. ks_xc_directional gives the derivative as a Lean HasDerivAt statement. The local functional derivative is an explicit hypothesis. On the equal-spin path \(D_\alpha=D_\beta=D/2\), the two perturbations are each \(H/2\). Their contractions give \((V^\alpha+V^\beta)/2\) against \(H\) (ks_equal_spin_pullback). For an exchange-symmetric spin functional the two potentials coincide there. The native equal-spin check tests both potential matrices as well as the energy.

Total energy accounting

For pure molecular KS, \(F_s=h+J[D_{\rm total}]+V^s_{\rm xc}\), where \(h\) is the core Hamiltonian. Write \(\langle D,F\rangle\) for the ordered AO contraction, summed over the applicable channels. The total energy is

\[ E=\langle D,h\rangle+\tfrac12\langle D,J\rangle+E_{\rm xc}+E_{\rm nuc}. \]

Starting with \(\langle D,F\rangle\) requires subtracting half the Coulomb contraction and the full XC potential contraction, then adding the XC functional and nuclear energies. ks_energy_ledger proves the equality of these two expressions. Production sums the individual energy components. The theorem assumes one evaluated density throughout the ledger.

Implementation and numerical checks

The RKS and UKS grid builders in cpp/src/rks.cpp and cpp/src/uks.cpp supply densities and invariants to the functional, then project its derivatives. The checker calls the corresponding native entry points evaluate_rks_xc_potential and evaluate_uks_xc_potential. It uses two hydrogen STO-3G AOs, four fixed points with unequal weights, and distinct, positive-definite alpha and beta density matrices. It independently builds the AO pair products and their gradients. This reference shares AO values and pointwise libxc evaluations with the native builder. A separate check differences the native XC energies, perturbing every independent symmetric AO entry in each channel with step \(10^{-5}\). An off-diagonal perturbation changes both entries, so its derivative includes both weights.

Native calculation	Maximum absolute AO potential discrepancy	Maximum absolute energy-derivative discrepancy
RKS LDA	2.78e-17	1.88e-12
UKS LDA	8.33e-17	4.32e-12
RKS PBE	2.78e-17	4.30e-12
UKS PBE	5.56e-17	3.17e-12

The projection and centered-difference tolerances were \(10^{-12}\) and \(10^{-9}\) respectively. In the open-spin PBE fixture, doubling the mixed coefficient changed the alpha AO potential by \(4.46\times10^{-3}\). Swapping the spin potentials changed it by \(2.96\times10^{-2}\). These controls test conventions that an equal-spin calculation alone cannot distinguish.

The complete energy checks reuse converged, all-electron, unfitted PBE/STO-3G calculations. ORCA 6.1.1 used explicit libxc PBE exchange and correlation, NoRI, DefGrid3 and VeryTightSCF. The compared vibe-qc grid used 150 radial points and Lebedev order 59. Both codes used the same geometries and spin.

Method and system	vibe-qc energy (Hartree)	ORCA energy (Hartree)	Absolute difference (Hartree)
RKS, neutral singlet water	-75.2215282935867	-75.22152812721895	1.66e-7
UKS, neutral doublet OH	-74.56983247538918	-74.5698324769221	1.53e-9

In bohr the water positions were O at the origin, H at \((0,0,1.81)\) and H at \((1.71,0,-0.45)\). OH uses the first O-H pair. Both differences pass the original \(10^{-5}\) Hartree tolerance. The energy components also reconstruct both native totals within \(2\times10^{-14}\) Hartree on both recorded grids (99/41 and 150/59). The RKS and UKS density-field functions, XC slice builders and fixed-density entry points had identical source text in the fixture runtime’s source revision, the SCF runtime’s source revision and the main revision examined here. The energy differences include the effects of using different quadratures.

Run python formal/check_proofs.py with the pinned Lean environment and python formal/check_ks.py with a native vibe-qc installation that exposes both XC evaluators. The nine mathematical claims apply to exact real sums at fixed geometry and grid. Pointwise smoothness and the accuracy of the functional derivatives remain assumptions. Density cutoffs, numerical integration error, floating-point execution, hybrid exchange, meta-GGA, nonlocal functionals, ROKS and nuclear gradients require separate evidence.