Tests

Project.toml

@@ -16,6 +16,7 @@ Libxc = "66e17ffc-8502-11e9-23b5-c9248d0eb96d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
@@ -38,6 +39,7 @@ KrylovKit = "0.9.5"
 Libxc = "0.3.19"
 Optim = "1.11.0"
 Printf = "1.11.0"
+QuadGK = "2.11.3"
 julia = "1.9"
 
 [extras]

examples/atomic density comparison/run.jl

@@ -0,0 +1,83 @@
+cd(@__DIR__)
+using Pkg
+Pkg.activate(joinpath(@__DIR__, ".."))
+using AtomicKohnSham
+using CairoMakie
+using Libxc
+
+# Ground-state radial densities for a handful of atoms, overlaid on one
+# plot. No log.txt/results.txt here -- just solve and plot.
+
+# ===================== STEP 1 : PARAMETERS =====================#
+Zs = [3, 16, 17, 18, 20, 26]  # Li, S, Cl, Ar, Ca, Fe
+
+# DISCRETIZATION : P1IntLegendreBasis + expmesh (same for every atom)
+Rmax        = 1000
+Nmesh       = 30
+ordermax    = 10
+
+# ALGORITHM : ODA + Optimized Aufbau
+# scftol=1e-10 with the default adaptive line-search mixing (tinit=0.4,
+# frozen_t=false) is enough on its own to bring the density tail down to
+# its Float64-precision floor (~1e-31, i.e. eps(Float64)^2 -- checked
+# directly: eval_density(argon, r=50) with these exact settings reaches
+# 4.4e-31 in just 26 iterations). A *looser* scftol leaks numerical noise
+# from the density matrix straight into the quadratic form eval_density
+# evaluates: at scftol=1e-8, rho(50) floors around 6e-19 and then
+# *increases* with r past that, which isn't physics.
+#
+# Iron alone needs different handling: its 3d/4s near-degeneracy makes this
+# adaptive mixing oscillate indefinitely (checked: stopping criterion
+# bounces between 1e-9 and 1e-7 for hundreds of iterations, regardless of
+# scftol). frozen_t=true with a fixed, fairly strong tinit tames that. Only
+# apply that fix to iron, though -- tried it on every atom first, and a
+# frozen mixing throttles each iteration's step enough that the *easy*
+# atoms can satisfy the stopping criterion (successive-iterate difference
+# below scftol) well before their density is actually converged, a false
+# positive that looked identical to real convergence until checked against
+# the density values directly.
+default_alg(aufbau) = ODA(tinit = 0.4, aufbau = aufbau, scftol = 1e-10)
+iron_alg(aufbau) = ODA(tinit = 0.15, frozen_t = true, aufbau = aufbau, scftol = 1e-12)
+
+default_maxiter = 100
+# Iron plateaus around stopping~2e-9 rather than cleanly reaching scftol,
+# so it reports MAXITERS below -- reported honestly rather than loosened
+# away, since even without satisfying the criterion exactly its density
+# still reaches ~1e-30, far better than a falsely-"successful" run would.
+iron_maxiter = 800
+
+# PLOT
+X = AtomicKohnSham.exprange(1e-3, 99.9, 10000; s = 1.5)
+
+# ===================== STEP 2 : BUILD & SOLVE =====================#
+"""
+Angular momentum cutoff and orbitals-per-channel, scaled with Z so heavier
+atoms (which need d orbitals, more radial nodes) get a big enough basis
+without over-paying for light ones.
+"""
+basis_size(Z) = (lh = min(ceil(Int, Z / 20), 3), nh = clamp(Z, 2, 10))
+
+sols = KSESolution[]
+for Z in Zs
+    lh, nh = basis_size(Z)
+
+    ex = Functional(:lda_x; n_spin = 1)
+    ec = NoFunctional(1)
+    model = KSEModel(; Z = Z, N = Z, ex = ex, ec = ec)
+
+    mesh = expmesh(0, Rmax, Nmesh; s = 1.2)
+    basis = P1IntLegendreBasis(mesh; ordermax = ordermax)
+    dis = KSEDiscretization(basis, model; lh = lh, nh = nh,
+        fem_integration_method = GaussLegendre(basis, 1000))
+
+    aufbau = OptimizedAufbau(tol = 1e-2)
+    alg = Z == 26 ? iron_alg(aufbau) : default_alg(aufbau)
+    maxiter = Z == 26 ? iron_maxiter : default_maxiter
+
+    sol = groundstate(model, dis, alg; maxiter = maxiter)
+    println(sol.name, ": ", sol.success, " in ", sol.niter, " iterations")
+    push!(sols, sol)
+end
+
+# ===================== STEP 3 : PLOT =====================#
+save("densities.pdf", plot_density(sols, X))

examples/basic_example_doublefloat/description.md

@@ -0,0 +1,79 @@
+# Description
+
+`examples/basic_example` redone with all FEM/SCF arithmetic carried out in
+`Double64` (~32 significant digits, via `DoubleFloats.jl`) instead of
+`Float64`: sodium (`Z=11, N=11`), `Rmax=500`, `Nmesh=20`, `ordermax=10`,
+`lh=2`, `nh=10`, `OptimizedAufbau(max_degen=2, tol=0.01)`,
+`ODA(tinit=0.6, scftol=1e-20)`, `maxiter=100`, same output pipeline
+(`log.txt`, `results.txt`, `density.pdf`, `orbitals.pdf`, `potentials.pdf`,
+`convergence.pdf`, `energy_breakdown.pdf`).
+
+**One parameter differs from `basic_example` on purpose: no exchange
+functional** (`ex = ec = NoFunctional(1)` instead of Slater exchange), for a
+reason worth spelling out, because it wasn't the original plan.
+
+## Why: Gauss-Legendre quadrature is Float64-only
+
+The first version of this example kept Slater exchange (via
+`BuiltinFunctional(:lda_x)`, since `Libxc.jl`'s `Functional` is a C-library
+wrapper that only accepts `Float64` — that fix is still in
+`BuiltinFunctional.jl` and still correct). But the total energy came out
+capped at Float64 precision regardless: tightening `scftol` from `1e-9` to
+`1e-20` barely moved the last few digits of `Etot`.
+
+The cause: `Ekin`/`Ecou`/`Ehar` are computed by **exact polynomial
+integration** and are genuinely `Double64`-precise. `Eexc`, though, is
+computed via `compute_exc_energy`, which integrates using quadrature nodes
+from `GaussLegendre` (`src/fem/integration methods.jl`). Those nodes come
+from `FastGaussQuadrature.gausslegendre(npoints)`, which has exactly one
+method and always returns `Float64` — there's no generic-precision variant.
+Worse, the failure is easy to miss: intermediate quantities like `Pgenx`
+*do* get computed at full `Double64` precision inside `GaussLegendre`'s
+constructor (confirmed by evaluating the same expression in isolation), but
+the struct's type parameter is pinned to `eltype(x)` — and `x` is the
+`Float64` array straight out of `gausslegendre`. The `Double64`-precise
+`Pgenx` gets silently truncated back to `Float64` the moment it's stored
+into the `GaussLegendre{Float64}` struct. So any model with an
+exchange-correlation functional is presently ceiling-capped at Float64
+precision no matter what type the rest of the pipeline uses. Fixing that for
+real means a generic-precision Gauss-Legendre node/weight generator — doable
+without a new dependency (`Legendre(n,T)` in `fem/legendre polynomial.jl` is
+already generic, so nodes could be its roots via Newton's method), but real
+numerical work, out of scope here.
+
+So: `hartree=1` is kept (the Hartree term is exact integration too — the `F`
+tensor, not quadrature — so it has no such problem), but there's no
+exchange-correlation functional. `GaussLegendre` is still constructed in
+`run.jl`, matching `basic_example`'s structure, but with no XC functional
+it's never actually touched by the computation.
+
+## Cross-validation
+
+Running the *same* Hartree-only, no-XC configuration in both precisions:
+
+|      | Float64            | Double64                          |
+|------|---------------------|------------------------------------|
+| Etot | -148.38279166622362 | -148.382791666223127793970294844... |
+
+These agree to essentially the full precision `Float64` can express (~16
+significant digits, right at `eps(Float64)`) — not partway through and then
+diverging, like the Slater-exchange version did. That's the direct
+confirmation: with the quadrature-dependent term removed, `Double64`
+actually delivers ~30 correct digits, `Float64` reproduces as much of that
+as its own precision allows, and the earlier, XC-including version's
+mismatch really was the `GaussLegendre` cap and nothing else.
+
+One remaining wrinkle, for anyone reading `results.txt` closely: the
+"density norm error" sanity check (`∫4πr²ρdr - N`) still shows `~2.7e-9`,
+not machine-`Double64`-small like the other four sanity checks. That check
+independently re-integrates the density via `GaussLegendre`'s (`Float64`)
+quadrature nodes purely for verification — it's a limitation of that
+diagnostic itself, not of the actual SCF computation, which the other four
+checks (energy balance, electron count, orthonormality, and a stopping
+criterion that reaches `1.6e-30`) confirm is genuinely full-precision.
+
+`run.jl` activates the shared `examples/` environment itself
+(`Pkg.activate(joinpath(@__DIR__, ".."))`), so just run:
+```julia
+julia run.jl
+```