vibe-qc/Documentation

Molecular restricted open-shell Kohn-Sham DFT

In vibe-qc’s molecular ROKS driver, the two spins share spatial orbitals but have separate densities and physical KS operators. The UKS result supplies the fixed-grid energy derivative with respect to those densities. The ROHF result explains the single effective operator used to update the orbitals. These results coincide as stationarity conditions for ordinary integer shell occupations; for a fractional frontier they need not.

Spin densities, energy and shell projectors

Take real orbitals \(C\) in a fixed AO basis with positive definite overlap \(S\), \(C^TSC=I\), at fixed nuclei. For spin \(\sigma\), let \(D_\sigma=C\,\mathrm{diag}(n^\sigma)\,C^T\). A closed orbital has occupations \((1,1)\) and a virtual orbital \((0,0)\). For a frontier of size \(m>0\), total occupation \(q_i\) on orbital \(i\) and spin excess \(\Delta=N_\alpha-N_\beta\), the implemented split is

\[ n_i^\alpha=(q_i+\Delta/m)/2,\qquad n_i^\beta=(q_i-\Delta/m)/2. \]

Each pair sums to \(q_i\), and the frontier carries total spin excess \(\Delta\). Admissibility (\(0\le n_i^\sigma\le1\)) and the correct total electron count are input conditions. For example, two orbitals each with \(q_i=1.5\) carry one excess electron as \((1,1/2)\) on each orbital.

Set \(D=D_\alpha+D_\beta\) and use the real AO pairing \(\langle A,B\rangle=\sum_{\mu\nu} A_{\mu\nu}B_{\mu\nu}\). For pure LDA/GGA on a fixed grid, the electronic energy and physical spin operators are

\[ E=\langle h,D\rangle+\tfrac12\langle J[D],D\rangle+E_{\rm xc}, \qquad F_\sigma=h+J[D]+V_{\rm xc}^\sigma . \]

Coulomb self-adjointness gives the full \(J[D]\) in the derivative; \(V_{\rm xc}^\sigma\) is the derivative of the same integrated XC energy. The UKS derivation gives \(\delta E=\sum_\sigma\langle F_\sigma,\delta D_\sigma\rangle\) at a differentiable point. Substituting the two physical operators into \(\tfrac12\sum_\sigma\langle h+F_\sigma,D_\sigma\rangle\), subtracting \(\tfrac12\sum_\sigma\langle V_{\rm xc}^\sigma,D_\sigma\rangle\) and adding \(E_{\rm xc}\) recovers the electronic energy above. Nuclear repulsion is added once to report the total.

For a shell \(s\), its selector is \(P_s=C M_s C^T S\) for a binary diagonal mask \(M_s\) (use \(P_v=I-P_c-P_o\) for the complement). Acting on \(C\) returns \(C M_s\). A spin-weighted mask instead returns \(C\,\mathrm{diag}(n^\sigma)\); replacing the open projector by a fractional density would scale some columns. The explicit partition builder uses the binary masks even for fractional densities. With \(\bar F=(F_\alpha+F_\beta)/2\), its MO blocks are:

Shell pair	Operator between the corresponding MO columns
closed/open	\(F_\beta\)
open/virtual	\(F_\alpha\)
closed/virtual; all same-shell blocks	\(\bar F\)

The reverse blocks are transposes. This follows by inserting the ROHF projectors on the left and right of each block; it does not use an energy derivative. The effective matrix selects orbital updates and is not substituted for either physical spin operator in the energy.

What does orbital diagonality establish?

Rotate two columns by \(C_i(t)=C_i\cos t+C_j\sin t\) and \(C_j(t)=C_j\cos t-C_i\sin t\). Holding occupations fixed, expansion of the two outer products gives

\[ D'_\sigma(0)=(n_i^\sigma-n_j^\sigma)(C_j C_i^T+C_i C_j^T). \]

Apply the UKS derivative to this path and use symmetry of each real \(F_\sigma\). Writing \(F^\sigma_{ij}=C_i^T F_\sigma C_j\) gives

\[ E'(0)=2\big[(n_i^\alpha-n_j^\alpha)F^\alpha_{ij} +(n_i^\beta-n_j^\beta)F^\beta_{ij}\big]. \]

For integer closed/open, open/virtual and closed/virtual rotations, the slopes are respectively \(2F^\beta_{ij}\), \(2F^\alpha_{ij}\) and \(4\bar F_{ij}\). Same-shell slopes vanish because the occupations agree. Thus zero intershell effective couplings imply zero first derivatives for all these shared rotations. This is a first-order condition, not a minimum or an SCF convergence proof.

For a \((1,1/2)\) frontier orbital rotating into an empty virtual, the slope is \(2F^\alpha_{ij}+F^\beta_{ij}\), whereas the effective coupling is still \(F^\alpha_{ij}\). For example, matrix elements \(0\) and \(1\) give zero effective coupling and unit slope. Other constraints or molecular symmetry may make the extra term vanish in a particular state. A fractional-shell stationarity check must evaluate the occupation-weighted slopes for the allowed rotations; a small effective-Fock residual alone is insufficient.

Proofs and implementation checks

The pinned Lean project adds ten real finite-dimensional theorems in ROKS.lean:

Theorems (prefix roks_)	Statement
frontier_spin, weighted_action	Spin split and action of a weighted mask in the overlap metric
open_virtual	ROHF block selection applied to the physical spin KS matrices
rotated_entry, turn_derivative	Per-entry rotation identity and its derivative at zero
rotation_derivative, cross_contraction	Chain rule given the physical energy differential; contraction of the symmetric pair into an MO element
integer_slopes, fractional_open_virtual	Occupation factors for the three integer intershell rotations and the fractional counterexample formula
half_trace	Electronic energy reconstruction over spin and AO entries

python formal/check_proofs.py builds the project and audits the transitive axioms. For the numerical bridge run python formal/check_roks.py with a native vibe-qc installation supporting molecular PBE. The check uses the production spin partition and effective Fock from python/vibeqc/rohf.py and the physical Fock and energy builder from python/vibeqc/roks.py. It compares both spin Focks with direct native XC and Coulomb evaluations and reconstructs the energy and all sixteen MO block entries. A nonsymmetric four-H geometry in STO-3G supplies the nonorthogonal basis; the four-point grid and prescribed densities are algebra tests, not molecular predictions.

Fixed occupation fixture	Maximum block error (\(E_h\))	Energy ledger error (\(E_h\))	Maximum slope error, six centered differences (\(E_h\)/rad)
Integer	2.22e-16	0	4.91e-10
Fractional	2.22e-16	1.33e-15	5.22e-10

The rotation step was \(2\times10^{-5}\) radians and the slope tolerance was \(2\times10^{-8}\) Hartree/radian. Substituting twice the effective coupling for the physical slope on the selected open/virtual rotation misses by 0.189929 Hartree/radian in the fractional fixture, versus \(1.82\times10^{-10}\) for integer occupations. This distinguishes the two conventions without presuming that a particular converged molecule has a residual.

Existing full-driver calculations with the same Python ROKS and ROHF module contents provide separate integration checks:

Case	vibe-qc total (\(E_h\))	Comparison
Closed-shell water, PBE/STO-3G	-75.22152819233663	ORCA 6.1.1 RKS, explicit libxc PBE, NoRI: absolute difference 6.51e-8 \(E_h\)
H doublet, LDA/6-31G	-0.47604447049266124	ORCA UKS, libxc IDs 1 and 7: absolute difference 3.66e-12 \(E_h\)
OH doublet, PBE/def2-SVP, 0.97 angstrom	-75.57307171512876	12 iterations, occupations 1.5/1.5 on the frontier, spin populations 5/4; no matched external ensemble reference

The first two are limits where RKS or UKS is a matched comparison; both pass the original \(10^{-5}\) Hartree tolerance. The OH calculation exercises the fractional restricted route; its reported convergence does not independently establish variational stationarity. See Lehtola et al., Sec. 9, Eqs. 44-46 for the spin orbital-rotation background; the shared-orbital and fractional occupation factors above follow from the displayed density path.

The proofs assume real arithmetic, fixed occupations, fixed geometry and grid, symmetric physical operators and a differentiable energy with the stated derivative. They do not cover occupation switches, grid screening thresholds, quadrature error, floating-point bounds, hybrids, nuclear gradients or state-specific excited-singlet ROKS. The numerical checks exercise selected production paths; the Lean model is not a verification of the whole engine.