Sanity Checks

examples/basic_example/results.txt

@@ -52,3 +52,12 @@ Occupation number =
             2p : ε = -1.006028995427625e+00        n = 6.0
             3s : ε = -7.701608721789199e-02        n = 0.999999999991612
             3p : ε = -1.111637731682302e-02        n = 8.387956995647983e-12
+
+SANITY CHECKS
+-------------
+Sanity checks = 
+            energy balance (Etot - ΣE)          = 0.0e+0
+            electron count error (Σn - N)       = 0.0e+0
+            density norm error (∫4πr²ρdr - N)   = 3.962675521052006e-09
+            orthonormality error (‖UᵀM₀U - I‖∞) = 1.285377138060539e-09
+            virial ratio (Ekin / -Etot)         = 9.999999999576582e-01

src/AtomicKohnSham.jl

@@ -141,6 +141,8 @@ module AtomicKohnSham
     export KSESolution
     export eval_orbital, eval_density, eval_hartree, eval_nuclear,
         eval_kinetic_potential, eval_vxc, eval_effective_potential
+    export SanityChecks
+    include("solution/sanity.jl")
     include("solution/types.jl")
     include("solution/show.jl")
     include("solution/postprocess.jl")

src/discretization/assemble.jl

@@ -38,6 +38,11 @@ function assemble_hartree_pot!(discretization::KSEDiscretization,
     @unpack A, F = discretization.femops
     @unpack Hartree = discretization.ksham
     @unpack B, W = cache.hartw
+    # HartreeWorkspace allocates B/W with length 0 when model.hartree == 0
+    # (see discretization.jl); A\B below would otherwise throw a
+    # DimensionMismatch. W is already empty in that case, so there's nothing
+    # left to do.
+    isempty(B) && return nothing
     if nspin == 1
         tensor_matrix_dict!(B, D, F)
         W .= A\B

src/discretization/energies.jl

@@ -97,6 +97,9 @@ function compute_hartree_energy(discretization::KSEDiscretization, D::AbstractAr
     @unpack Rmax, cache, nspin, N = discretization
     @unpack A, F = discretization.femops
     @unpack B, W = cache.hartw
+    # HartreeWorkspace allocates B/W with length 0 when model.hartree == 0 (see
+    # discretization.jl); A\B below would otherwise throw a DimensionMismatch.
+    isempty(B) && return zero(eltype(D))
     if nspin == 1
         tensor_matrix_dict!(B, D, F)
         W .= A\B
@@ -118,6 +121,8 @@ function compute_hartree_mix_energy(discretization::KSEDiscretization,
     @unpack Rmax, cache, nspin, N = discretization
     @unpack A, F = discretization.femops
     @unpack B, W = cache.hartw
+    # See compute_hartree_energy: B/W are empty when model.hartree == 0.
+    isempty(B) && return zero(eltype(D0))
     if nspin == 1
         tensor_matrix_dict!(B, D0, F)
         W .= A\B

src/solution/report.jl

@@ -172,6 +172,11 @@ function write_report(sol::KSESolution{T}, path::AbstractString) where T
         println(io, "OCCUPIED ORBITALS")
         println(io, "-----------------")
         _print_occupied(io, sol; styled = false)
+        println(io)
+
+        println(io, "SANITY CHECKS")
+        println(io, "-------------")
+        _print_sanity(io, sol; styled = false)
     end
     nothing
 end

src/solution/sanity.jl

@@ -0,0 +1,96 @@
+"""
+    SanityChecks{T}
+
+Physical consistency diagnostics computed from a converged (or stopped)
+Kohn–Sham solution, independent of any comparison to a previous run's output.
+
+Each field checks an identity or conservation law that must hold regardless of
+the specific system/functional/discretization. A nonzero value here indicates
+an actual computational problem (in FEM assembly, the eigensolver, the aufbau
+scheme, or energy bookkeeping) even when the SCF loop reports `SUCCESS` and no
+exception is thrown along the way — the kind of bug a golden-value regression
+test on the final energy alone cannot distinguish from "the physics is right".
+
+# Fields
+- `energy_balance::T`: `Etot - (Ekin+Ecou+Ehar+Eexc)`. Zero by definition;
+  any deviation means the energy components don't actually add up to the
+  reported total.
+- `electron_count_error::T`: `∑ᵢnᵢ - N`. Zero: occupation numbers must sum to
+  the model's electron count.
+- `density_norm_error::T`: `∫4πr²ρ(r)dr - N`, evaluated by numerical
+  quadrature independently of the occupation numbers. Zero: cross-checks that
+  the density matrix is consistent with the occupations, through a different
+  code path than `electron_count_error`.
+- `orthonormality_error::T`: `maxₗ,σ ‖Uₗ,σᵀ M₀ Uₗ,σ - I‖∞`, where `M₀` is the
+  FEM mass matrix. Zero: the Kohn–Sham orbitals within each angular
+  momentum/spin channel must be orthonormal.
+- `virial_ratio::T`: `Ekin / (-Etot)`. By the quantum virial theorem this goes
+  to exactly `1` whenever the only potentials present are homogeneous of
+  degree -1 in `r` — which includes *both* the nuclear Coulomb term and the
+  Hartree term (any `hartree` coefficient), just not exchange–correlation
+  (LDA-type functionals scale differently, e.g. as `ρ^(4/3)`, so no XC means
+  no closed-form target). Unlike the other fields above, this is *not* a
+  fixed-tolerance bookkeeping identity: in a finite radial domain `[0,Rmax]`
+  it converges to `1` only as `Rmax → ∞` (the deviation is a domain-truncation
+  effect — empirically it shrinks rapidly with `Rmax`, consistent with cutting
+  off the exponentially decaying tail of a bound state, and is largely
+  insensitive to basis/polynomial order at fixed `Rmax`), so treat it as a
+  convergence diagnostic with a loose tolerance, not an exact check.
+
+Any field is `NaN` if the underlying quantity (`D`, `U`, or `n`) is
+unavailable from the algorithm cache used to solve the problem.
+"""
+struct SanityChecks{T <: Real}
+    energy_balance::T
+    electron_count_error::T
+    density_norm_error::T
+    orthonormality_error::T
+    virial_ratio::T
+end
+
+function _density_norm_error(discretization::KSEDiscretization{T}, D) where {T}
+    isnothing(D) && return T(NaN)
+    @unpack fem_integration_method, nspin, N = discretization
+    @unpack y, wy = fem_integration_method
+    ρ = similar(y)
+    if nspin == 1
+        eval_density!(ρ, discretization, D, y)
+    else
+        @views DUP = D[:, :, 1]
+        @views DDOWN = D[:, :, 2]
+        ρup = similar(y)
+        ρdown = similar(y)
+        eval_density!(ρup, discretization, DUP, y)
+        eval_density!(ρdown, discretization, DDOWN, y)
+        @. ρ = ρup + ρdown
+    end
+    @. ρ *= y^2
+    4π * dot(wy, ρ) - N
+end
+
+function _orthonormality_error(discretization::KSEDiscretization{T}, U) where {T}
+    isnothing(U) && return T(NaN)
+    @unpack lₕ, nₕ, nspin = discretization
+    M₀ = discretization.femops.M₀
+    err = zero(T)
+    for σ in 1:nspin
+        for l in 1:(lₕ + 1)
+            Ulσ = nspin == 1 ? (@view U[:, :, l]) : (@view U[:, :, l, σ])
+            gram = Ulσ' * M₀ * Ulσ
+            err = max(err, maximum(abs, gram - I))
+        end
+    end
+    err
+end
+
+function _sanity_checks(discretization::KSEDiscretization{T}, model::KSEModel,
+        D, U, n, energies::Energies) where {T}
+    energy_balance = energies.Etot -
+                      (energies.Ekin + energies.Ecou + energies.Ehar + energies.Eexc)
+    electron_count_error = isnothing(n) ? T(NaN) : sum(n) - model.N
+    density_norm_error = _density_norm_error(discretization, D)
+    orthonormality_error = _orthonormality_error(discretization, U)
+    virial_ratio = iszero(energies.Etot) ? T(NaN) : energies.Ekin / (-energies.Etot)
+    SanityChecks{T}(energy_balance, electron_count_error, density_norm_error,
+        orthonormality_error, virial_ratio)
+end

src/solution/show.jl

@@ -51,6 +51,28 @@ function _print_occupied(io::IO, sol::KSESolution{T}; styled::Bool = true) where
     end
 end
 
+"""
+Print the "Sanity checks" block (energy balance, electron count, density
+normalization, orbital orthonormality, virial ratio) shared by
+`Base.show(io, sol)` and `write_report`. See [`SanityChecks`](@ref) for what
+each value means and why a nonzero error indicates a real problem.
+"""
+function _print_sanity(io::IO, sol::KSESolution{T}; styled::Bool = true) where {T}
+    if styled
+        printstyled(io, "Sanity checks = \n"; bold = true, color = :blue)
+    else
+        println(io, "Sanity checks = ")
+    end
+    sc = sol.sanity
+    _write_kv_block(io, 12, [
+        ("energy balance (Etot - ΣE)", _sci_str(sc.energy_balance)),
+        ("electron count error (Σn - N)", _sci_str(sc.electron_count_error)),
+        ("density norm error (∫4πr²ρdr - N)", _sci_str(sc.density_norm_error)),
+        ("orthonormality error (‖UᵀM₀U - I‖∞)", _sci_str(sc.orthonormality_error)),
+        ("virial ratio (Ekin / -Etot)", _sci_str(sc.virial_ratio)),
+    ])
+end
+
 function Base.show(io::IO, sol::KSESolution{T}) where {T}
     printstyled(io, "Name : $(sol.name)\n"; bold = true)
     printstyled(io, "Success = "; bold = true)
@@ -62,4 +84,6 @@ function Base.show(io::IO, sol::KSESolution{T}) where {T}
     println(io, sol.stopping_criteria)
 
     _print_occupied(io, sol)
+    println(io)
+    _print_sanity(io, sol)
 end

src/solution/types.jl

@@ -62,6 +62,9 @@ solution.
 - `occupied`: Summary of occupied orbitals, sorted by orbital energy. Each entry
   contains the shell label, orbital energy, and occupation number.
 - `W`: Coefficients of the Hartree potential.
+- `sanity::SanityChecks`: Physical consistency diagnostics (energy bookkeeping,
+  electron count, density normalization, orbital orthonormality, virial
+  ratio) computed independently of any comparison to a previous run's output.
 - `logbook::LogBook`: Logbook containing recorded quantities during the SCF
   iterations.
 - `name::String`: Name of the solution, used for identification or labeling.
@@ -85,6 +88,8 @@ struct KSESolution{T <: Real, TD, TU, TE, TN, TO, TW, logbookType <: LogBook, C<
 
     W::TW                               # Coefficients of the Hartree potential
 
+    sanity::SanityChecks{T}             # Physical consistency diagnostics
+
     logbook::logbookType                # LogBook of all recorded quantities
 
     name::String                        # Name of the solution
@@ -108,6 +113,8 @@ struct KSESolution{T <: Real, TD, TU, TE, TN, TO, TW, logbookType <: LogBook, C<
 
         W = discretization.cache.hartw.W
 
+        sanity = _sanity_checks(discretization, model, D, U, n, solver.energies)
+
         # Context
         context = KSEContext(model, alg, lₕ, nₕ, basis, fem_integration_method)
 
@@ -125,6 +132,7 @@ struct KSESolution{T <: Real, TD, TU, TE, TN, TO, TW, logbookType <: LogBook, C<
                              solver.energies,
                              D, U, ϵ, n, occupied,
                              W,
+                             sanity,
                              solver.logbook,
                              name,
                              context)

test/hydrogen.jl

@@ -22,6 +22,35 @@ function run_solve(ex::TEX, ec::TEC,
     return sol
 end
 
+# ===================================================================
+#                             SANITY CHECKS
+# ===================================================================
+"""
+    check_sanity(sol; tol=1e-6, virial_tol=nothing)
+
+Assert the Tier-1 physical consistency checks in `sol.sanity` (energy
+balance, electron count, density normalization, orbital orthonormality) --
+identities/conservation laws that must hold regardless of the system or
+functional, so a violation means a real computational problem even though the
+SCF loop reported `SUCCESS`.
+
+If `virial_tol` is given, also assert the quantum virial theorem
+`Ekin ≈ -Etot` (Tier 2b). This only holds when there is no exchange–
+correlation functional (Hartree and nuclear-Coulomb are both homogeneous of
+degree -1 in `r`, XC generally is not), and only up to the domain-truncation
+defect that vanishes as `Rmax → ∞` -- so it must never be asserted for a
+functional-bearing model, and needs a looser tolerance than the Tier-1 checks.
+"""
+function check_sanity(sol; tol::Real = 1e-6, virial_tol::Union{Real, Nothing} = nothing)
+    sc = sol.sanity
+    @test abs(sc.energy_balance) < tol
+    @test abs(sc.electron_count_error) < tol
+    @test abs(sc.density_norm_error) < tol
+    @test abs(sc.orthonormality_error) < tol
+    if !isnothing(virial_tol)
+        @test abs(sc.virial_ratio - 1) < virial_tol
+    end
+end
 
 @testset "Hydrogen RHF Spin Nopoloarized" begin
 
@@ -48,6 +77,10 @@ end
     @test occ[1] == "1s"
     @test abs(occ[2] + 0.046222) < 1e-5
     @test occ[3] == 1.0
+
+    # No exchange-correlation here (only nuclear Coulomb + Hartree, both ~1/r):
+    # the virial theorem Ekin = -Etot applies, up to the Rmax-truncation defect.
+    check_sanity(sol; virial_tol = 1e-5)
 end
 
 
@@ -77,6 +110,9 @@ end
     @test occ[1] == "1s"
     @test abs(occ[2] + 0.194250) < 1e-6
     @test occ[3] == 1.0
+
+    # Exchange functional is present here, so no virial identity to check.
+    check_sanity(sol)
 end
 
 
@@ -95,4 +131,41 @@ end
 
     @test_nowarn run_solve(ex, ec, Rmax, Nmesh, ordermax, maxiter, scftol)
     sol = run_solve(ex, ec, Rmax, Nmesh, ordermax, maxiter, scftol)
+
+    check_sanity(sol)
+end
+
+@testset "Hydrogen Bare Coulomb (exact)" begin
+    # No Hartree (hartree=0) and no exchange-correlation: the bare
+    # one-electron Coulomb problem, which is exactly solvable in closed form.
+    # Unlike every other testset here, this is not a regression check against
+    # a value this code happened to produce in the past -- it's checked
+    # against E_nl = -Z^2/(2n^2), a fact about hydrogen that holds
+    # independently of this codebase.
+    ex = NoFunctional(n_spin = 1)
+    ec = NoFunctional(n_spin = 1)
+    model = KSEModel(; Z = 1, N = 1, hartree = 0, ex = ex, ec = ec)
+    mesh = expmesh(0, 1000, 30; s = 1.2)
+    basis = P1IntLegendreBasis(mesh; ordermax = 10)
+    dis = KSEDiscretization(basis, model; lh = 1, nh = 2)
+    aufbau = OptimizedAufbau(max_degen = 1)
+    alg = ODA(tinit = 0.6, aufbau = aufbau, scftol = 1e-12)
+
+    @test_nowarn groundstate(model, dis, alg; maxiter = 200)
+    sol = groundstate(model, dis, alg; maxiter = 200)
+    @test sol.success == "SUCCESS"
+
+    # E_1s = -Z^2/2 exactly.
+    @test isapprox(sol.energies.Etot, -0.5; atol = 1e-8)
+
+    # sol.ϵ[l+1, k] = E_{n=k+l, l}; 2s and 2p must match the accidental
+    # l-degeneracy of the pure Coulomb problem exactly.
+    @test isapprox(sol.ϵ[1, 1], -1/2;  atol = 1e-8)   # 1s (n=1)
+    @test isapprox(sol.ϵ[1, 2], -1/8;  atol = 1e-8)   # 2s (n=2)
+    @test isapprox(sol.ϵ[2, 1], -1/8;  atol = 1e-8)   # 2p (n=2)
+    @test isapprox(sol.ϵ[2, 2], -1/18; atol = 1e-8)   # 3p (n=3)
+
+    # Pure Coulomb: the virial theorem Ekin = -Etot must hold almost exactly
+    # (Rmax=1000 makes the domain-truncation defect negligible here).
+    check_sanity(sol; tol = 1e-8, virial_tol = 1e-8)
 end

test/scandium.jl

@@ -22,6 +22,18 @@ function run_solve(ex::TEX, ec::TEC,
     return sol
 end
 
+# See hydrogen.jl for the documented version of this helper.
+function check_sanity(sol; tol::Real = 1e-6, virial_tol::Union{Real, Nothing} = nothing)
+    sc = sol.sanity
+    @test abs(sc.energy_balance) < tol
+    @test abs(sc.electron_count_error) < tol
+    @test abs(sc.density_norm_error) < tol
+    @test abs(sc.orthonormality_error) < tol
+    if !isnothing(virial_tol)
+        @test abs(sc.virial_ratio - 1) < virial_tol
+    end
+end
+
 
 @testset "Scandium RHF Spin Nopoloarized" begin
 
@@ -84,4 +96,8 @@ end
     @test occ3d[1] == "3d"
     @test abs(occ3d[2] + 0.00262) < 1e-5
     @test abs(occ3d[3]/5 - 0.0056) < 1e-4
+
+    # No exchange-correlation here (only nuclear Coulomb + Hartree, both ~1/r):
+    # the virial theorem Ekin = -Etot applies, up to the Rmax-truncation defect.
+    check_sanity(sol; tol = 1e-5, virial_tol = 1e-5)
 end