diff --git a/README.md b/README.md index 226e993..dc904c4 100644 --- a/README.md +++ b/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 diff --git a/docs/make.jl b/docs/make.jl index 57353b7..f830c53 100644 --- a/docs/make.jl +++ b/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", diff --git a/docs/src/DB.md b/docs/src/DB.md new file mode 100644 index 0000000..092fdf5 --- /dev/null +++ b/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 +``` diff --git a/docs/src/FDM.md b/docs/src/FDM.md index 16352c9..adc6ce2 100644 --- a/docs/src/FDM.md +++ b/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 diff --git a/docs/src/Rayleigh-Ritz.md b/docs/src/Rayleigh-Ritz.md index bd4c47f..a60508a 100644 --- a/docs/src/Rayleigh-Ritz.md +++ b/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 diff --git a/docs/src/index.md b/docs/src/index.md index ed36578..bbc3b1b 100644 --- a/docs/src/index.md +++ b/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 ``` diff --git a/src/DB.jl b/src/DB.jl new file mode 100644 index 0000000..808ca1a --- /dev/null +++ b/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) diff --git a/src/TwoBody.jl b/src/TwoBody.jl index a455922..37adc95 100644 --- a/src/TwoBody.jl +++ b/src/TwoBody.jl @@ -3,6 +3,9 @@ module TwoBody # Hamiltonian include("./Hamiltonian.jl") +# Database +include("./DB.jl") + # Basis include("./Basis.jl") diff --git a/test/DB.jl b/test/DB.jl new file mode 100644 index 0000000..9ea22f2 --- /dev/null +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 89abff9..ad24e51 100644 --- a/test/runtests.jl +++ b/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 \ No newline at end of file +end