Add benchmark Hamiltonian database

README.md

@@ -16,17 +16,15 @@ config:
 ---
 flowchart TD
   A["Hamiltonian.jl"]
-  B["Basis.jl"]
   C["Rayleigh-Ritz.jl"]
-  D["GEM.jl"]
-  E["FiniteDifferenceMatrices.jl"]
   F["FDM.jl"]
   G["VMC.jl"]
+  H["DB.jl"]
   Z["TwoBody.jl"]
-  A --> C & D & F & G
-  B --> C & D
-  E --> F
-  C & D & F & G --> Z
+  A --> H
+  A --> C & F & G
+  H --> C & F & G
+  C & F & G --> Z
 ```
 
 ## Developer's Guide

docs/make.jl

@@ -18,6 +18,7 @@ makedocs(;
     pages=[
         "Home" => "index.md",
         "Hamiltonian" => "Hamiltonian.md",
+        "Database" => "DB.md",
         "Rayleigh-Ritz Method" => "Rayleigh-Ritz.md",
         "Finite Difference Method" => "FDM.md",
         "API reference" => "API.md",

docs/src/DB.md

@@ -0,0 +1,53 @@
+```@meta
+CurrentModule = TwoBody
+```
+
+# Database
+
+The internal database provides benchmark Hamiltonians and reference energies
+for testing solvers. Each Hamiltonian can be passed to a solver and its result
+compared with the reference energy. The database is not exported because it is
+intended for package development.
+
+## Usage
+
+Retrieve a benchmark using its key:
+
+```julia-repl
+julia> entry = TwoBody.db(:hydrogen)
+
+julia> entry.hamiltonian
+
+julia> entry.energy
+```
+
+Add a benchmark using `TwoBody.put!`:
+
+```julia-repl
+julia> hamiltonian = Hamiltonian(
+           NonRelativisticKinetic(ℏ = 1.0, m = 1.0),
+           CoulombPotential(coefficient = -1.0),
+       )
+
+julia> TwoBody.put!(:example, hamiltonian, -0.5)
+```
+
+Duplicate keys are rejected, and Hamiltonians are copied on registration and
+lookup to protect stored benchmarks.
+
+## Data
+
+| Key | System | Reference energy |
+|:--|:--|--:|
+| `:hydrogen` | Hydrogen ground state | `-0.5` |
+| `:positronium` | Positronium ground state | `-0.25` |
+| `:harmonic_oscillator` | Three-dimensional harmonic oscillator ground state | `1.5` |
+
+## API
+
+```@docs; canonical=false
+TwoBody.DatabaseEntry
+TwoBody.put!
+TwoBody.db
+TwoBody.dbkeys
+```

docs/src/FDM.md

@@ -94,7 +94,7 @@ println("------------------------------")
 println(" n     numerical    analytical")
 println("------------------------------")
 for n in 1:4
-  @printf("%2d  %+.9f  %+.9f\n", n, res.E[n], Antique.E(HA,n=n))
+  @printf("%2d  %+.9f  %+.9f\n", n, res.E[n], Antique.energy(HA,n=n))
 end
 
 # wave function
@@ -119,7 +119,7 @@ for n in 1:4
   X = res.method.R
   Y = 4π * X .^2 .* res.ψ[:,n] .^ 2
   scatter!(axis, X, Y, label="TwoBody.jl", markersize=6)
-  lines!(axis, 0..50, r -> 4π * r^2 * abs(Antique.ψ(HA,r,0,0,n=n))^2, label="Antique.jl", color=:black)
+  lines!(axis, 0..50, r -> 4π * r^2 * abs(Antique.wavefunction(HA,r,0,0,n=n))^2, label="Antique.jl", color=:black)
   axislegend(axis, "n = $n", position=:rt, framevisible=false)
 end
 save("assets/FDM_HA.svg", fig) # hide
@@ -149,7 +149,7 @@ println("------------------------------")
 println(" n     numerical    analytical")
 println("------------------------------")
 for n in 1:4
-  @printf("%2d  %+.9f  %+.9f\n", n-1, res.E[n], Antique.E(SO,n=n-1))
+  @printf("%2d  %+.9f  %+.9f\n", n-1, res.E[n], Antique.energy(SO,n=n-1))
 end
 
 # wave function
@@ -174,7 +174,7 @@ for n in 1:4
   X = res.method.R
   Y = 4π * X .^2 .* res.ψ[:,n] .^ 2
   scatter!(axis, X, Y, label="TwoBody.jl", markersize=6)
-  lines!(axis, 0..50, r -> 4π * r^2 * abs(Antique.ψ(SO,r,0,0,n=n-1))^2, label="Antique.jl", color=:black)
+  lines!(axis, 0..50, r -> 4π * r^2 * abs(Antique.wavefunction(SO,r,0,0,n=n-1))^2, label="Antique.jl", color=:black)
   axislegend(axis, "n = $(n-1)", position=:rt, framevisible=false)
 end
 fig

docs/src/Rayleigh-Ritz.md

@@ -93,7 +93,7 @@ println("------------------------------")
 println(" n     numerical    analytical")
 println("------------------------------")
 for n in 1:4
-  @printf("%2d  %+.9f  %+.9f\n", n, res.E[n], Antique.E(HA,n=n))
+  @printf("%2d  %+.9f  %+.9f\n", n, res.E[n], Antique.energy(HA,n=n))
 end
 
 # wave function
@@ -116,7 +116,7 @@ for n in 1:4
     )
   )
   lines!(axis, 0..50, r -> 4π * r^2 * abs(TwoBody.ψ(res,r,n=n))^2, label="TwoBody.jl")
-  lines!(axis, 0..50, r -> 4π * r^2 * abs(Antique.ψ(HA,r,0,0,n=n))^2, label="Antique.jl", color=:black, linestyle=:dash)
+  lines!(axis, 0..50, r -> 4π * r^2 * abs(Antique.wavefunction(HA,r,0,0,n=n))^2, label="Antique.jl", color=:black, linestyle=:dash)
   axislegend(axis, "n = $n", position=:rt, framevisible=false)
 end
 fig
@@ -147,7 +147,7 @@ println("------------------------------")
 println(" n     numerical    analytical")
 println("------------------------------")
 for n in 1:4
-  @printf("%2d  %+.9f  %+.9f\n", n-1, res.E[n], Antique.E(SO,n=n-1))
+  @printf("%2d  %+.9f  %+.9f\n", n-1, res.E[n], Antique.energy(SO,n=n-1))
 end
 
 # wave function
@@ -170,7 +170,7 @@ for n in 1:4
     )
   )
   lines!(axis, 0..50, r -> 4π * r^2 * abs(TwoBody.ψ(res,r,n=n))^2, label="TwoBody.jl")
-  lines!(axis, 0..50, r -> 4π * r^2 * abs(Antique.ψ(SO,r,0,0,n=n-1))^2, label="Antique.jl", color=:black, linestyle=:dash)
+  lines!(axis, 0..50, r -> 4π * r^2 * abs(Antique.wavefunction(SO,r,0,0,n=n-1))^2, label="Antique.jl", color=:black, linestyle=:dash)
   axislegend(axis, "n = $(n-1)", position=:rt, framevisible=false)
 end
 fig

docs/src/index.md

@@ -81,7 +81,7 @@ using CairoMakie
 fig = Figure(size=(420,300), fontsize=11, backgroundcolor=:transparent)
 axis = Axis(fig[1,1], xlabel=L"$r / a_0$", ylabel=L"$\psi(r) / a_0^{-3/2}$", ylabelsize=16.5, xlabelsize=16.5, limits=(0,4,0,1.1/sqrt(π)))
 lines!(axis, 0..5, r -> abs(TwoBody.ψ(res,r)), label="TwoBody.jl")
-lines!(axis, 0..5, r -> abs(Antique.ψ(HA,r,0,0)), linestyle=:dash, color=:black, label="Antique.jl")
+lines!(axis, 0..5, r -> abs(Antique.wavefunction(HA,r,0,0)), linestyle=:dash, color=:black, label="Antique.jl")
 axislegend(axis, position=:rt, framevisible=false)
 fig
 ```

src/DB.jl

@@ -0,0 +1,99 @@
+"""
+    DatabaseEntry(hamiltonian, energy)
+
+A benchmark problem stored in the database. `hamiltonian` is ready to be
+passed to a solver and `energy` is its reference energy.
+"""
+struct DatabaseEntry{T<:Real}
+    hamiltonian::Hamiltonian
+    energy::T
+end
+
+import Base: put!
+
+# Keeping the registry private makes `db` the single lookup boundary and leaves
+# room for validation, lazy loading, and provenance.
+const _DATABASE = Dict{Symbol,DatabaseEntry}()
+
+"""
+    put!(key, hamiltonian, energy)
+
+Add a benchmark problem to the database. `key` may be a `Symbol` or string,
+`hamiltonian` must be a [`Hamiltonian`](@ref), and `energy` must be real.
+Registering the same key twice throws an `ArgumentError`.
+"""
+function put!(key::Symbol, hamiltonian::Hamiltonian, energy::T) where {T<:Real}
+    haskey(_DATABASE, key) &&
+        throw(ArgumentError("database key $(repr(key)) is already registered"))
+
+    entry = DatabaseEntry(deepcopy(hamiltonian), energy)
+    _DATABASE[key] = entry
+    return entry
+end
+
+put!(key::AbstractString, hamiltonian::Hamiltonian, energy::Real) =
+    put!(Symbol(key), hamiltonian, energy)
+
+# PoC data in atomic units.
+put!(
+    :hydrogen,
+    Hamiltonian(
+        NonRelativisticKinetic(ℏ = 1.0, m = 1.0),
+        CoulombPotential(coefficient = -1.0),
+    ),
+    -0.5,
+)
+
+put!(
+    :positronium,
+    Hamiltonian(
+        NonRelativisticKinetic(ℏ = 1.0, m = 0.5),
+        CoulombPotential(coefficient = -1.0),
+    ),
+    -0.25,
+)
+
+put!(
+    :harmonic_oscillator,
+    Hamiltonian(
+        NonRelativisticKinetic(ℏ = 1.0, m = 1.0),
+        PowerLawPotential(coefficient = 0.5, exponent = 2.0),
+    ),
+    1.5,
+)
+
+"""
+    db(key::Union{Symbol,AbstractString}) -> DatabaseEntry
+
+Return the benchmark Hamiltonian and reference energy associated with `key`.
+The returned Hamiltonian is independent of the stored value and can safely be
+modified by callers.
+
+# Examples
+
+```julia
+entry = db(:hydrogen)
+result = solve(entry.hamiltonian, method)
+isapprox(result.values[1], entry.energy)
+```
+"""
+function db(key::Symbol)
+    haskey(_DATABASE, key) || throw(
+        ArgumentError(
+            "unknown database key $(repr(key)); available keys: " *
+            join(repr.(dbkeys()), ", "),
+        ),
+    )
+
+    entry = _DATABASE[key]
+    return DatabaseEntry(deepcopy(entry.hamiltonian), entry.energy)
+end
+
+db(key::AbstractString) = db(Symbol(key))
+
+"""
+    dbkeys() -> Vector{Symbol}
+
+Return the available database keys in deterministic order.
+"""
+dbkeys() = sort!(collect(keys(_DATABASE)); by = string)

src/TwoBody.jl

@@ -3,6 +3,9 @@ module TwoBody
 # Hamiltonian
 include("./Hamiltonian.jl")
 
+# Database
+include("./DB.jl")
+
 # Basis
 include("./Basis.jl")
 

test/DB.jl

@@ -0,0 +1,38 @@
+@testset "Database" begin
+    hydrogen = TwoBody.db(:hydrogen)
+
+    @test hydrogen isa TwoBody.DatabaseEntry{Float64}
+    @test hydrogen.hamiltonian isa Hamiltonian
+    @test hydrogen.energy == -0.5
+    @test TwoBody.db("hydrogen").energy == hydrogen.energy
+    @test TwoBody.dbkeys() == [:harmonic_oscillator, :hydrogen, :positronium]
+    @test_throws ArgumentError TwoBody.db(:unknown)
+
+    positronium = TwoBody.db(:positronium)
+    @test positronium.energy == -0.25
+    @test positronium.hamiltonian[1].m == 0.5
+
+    # A caller may modify its Hamiltonian without corrupting the registry.
+    pop!(hydrogen.hamiltonian.terms)
+    @test length(hydrogen.hamiltonian) == 1
+    @test length(TwoBody.db(:hydrogen).hamiltonian) == 2
+
+    key = :temporary_test_problem
+    hamiltonian = Hamiltonian(
+        NonRelativisticKinetic(ℏ = 1.0, m = 1.0),
+        CoulombPotential(coefficient = -2.0),
+    )
+
+    try
+        registered = TwoBody.put!(key, hamiltonian, -2.0)
+        @test registered isa TwoBody.DatabaseEntry{Float64}
+        @test TwoBody.db(key).energy == -2.0
+        @test_throws ArgumentError TwoBody.put!(key, hamiltonian, -2.0)
+
+        # Registration also isolates the stored Hamiltonian from the caller.
+        pop!(hamiltonian.terms)
+        @test length(TwoBody.db(key).hamiltonian) == 2
+    finally
+        delete!(TwoBody._DATABASE, key)
+    end
+end

test/runtests.jl

@@ -7,7 +7,8 @@ using SpecialFunctions
 using ForwardDiff
 
 @testset verbose = true "TwoBody.jl" begin
+  include("DB.jl")
   include("Basis.jl")
   include("Rayleigh-Ritz.jl")
   include("FDM.jl")
-end
+end