From 27a3310e20078b03ab3dacb50c9fdefd051f493b Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Tue, 28 Oct 2025 18:15:51 +0100 Subject: [PATCH 01/15] updated exports --- Examples/HydrogenAnion.jl | 2 +- src/FewBodyECG.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Examples/HydrogenAnion.jl b/Examples/HydrogenAnion.jl index 78ec4c5..0d0240f 100644 --- a/Examples/HydrogenAnion.jl +++ b/Examples/HydrogenAnion.jl @@ -27,4 +27,4 @@ E = -0.527751016523 @info "Energy difference" ΔE n, E = convergence(A) -plot(n, E) +plot(n, E) \ No newline at end of file diff --git a/src/FewBodyECG.jl b/src/FewBodyECG.jl index 49de1de..c537865 100644 --- a/src/FewBodyECG.jl +++ b/src/FewBodyECG.jl @@ -15,7 +15,7 @@ export Operator export build_hamiltonian_matrix, build_overlap_matrix, solve_generalized_eigenproblem, solve_ECG, convergence -export ψ₀ +export ψ₀, SolverResults, convergence include("types.jl") include("coordinates.jl") From 86278161e0d81f2dff2969f84731151c60d5d9c4 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Wed, 29 Oct 2025 07:17:15 +0100 Subject: [PATCH 02/15] updated types and tests for shifts --- src/types.jl | 17 ++++++------ test/test_hamiltonian.jl | 52 ++++++++++++++++++++++++++---------- test/test_matrix_elements.jl | 14 +++++++--- 3 files changed, 57 insertions(+), 26 deletions(-) diff --git a/src/types.jl b/src/types.jl index 31fc39f..e848bcc 100644 --- a/src/types.jl +++ b/src/types.jl @@ -3,18 +3,19 @@ using FewBodyHamiltonians abstract type GaussianBase end struct Rank0Gaussian <: GaussianBase - A::Matrix{Float64} + A::Matrix{Number} + s::Vector{Number} end struct Rank1Gaussian <: GaussianBase - A::Matrix{Float64} - a::Vector{Float64} + A::Matrix{Number} + a::Vector{Number} end struct Rank2Gaussian <: GaussianBase - A::Matrix{Float64} - a::Vector{Float64} - b::Vector{Float64} + A::Matrix{Number} + a::Vector{Number} + b::Vector{Number} end struct BasisSet @@ -22,12 +23,12 @@ struct BasisSet end struct KineticOperator <: FewBodyHamiltonians.KineticTerm - K::Matrix{Float64} + K::Matrix{Number} end struct CoulombOperator <: FewBodyHamiltonians.PotentialTerm coefficient::Float64 - w::Vector{Float64} + w::Vector{Number} end struct ECG diff --git a/test/test_hamiltonian.jl b/test/test_hamiltonian.jl index 50a3225..2762a68 100644 --- a/test/test_hamiltonian.jl +++ b/test/test_hamiltonian.jl @@ -10,8 +10,11 @@ import FewBodyECG: _compute_overlap_element, _build_operator_matrix, _compute_ma A = [1.0 0.2; 0.2 1.5] B = [0.9 0.1; 0.1 1.2] - bra = Rank0Gaussian(A) - ket = Rank0Gaussian(B) + s₁ = randn(2) + s₂ = randn(2) + + bra = Rank0Gaussian(A,s₁) + ket = Rank0Gaussian(B,s₂) val = _compute_overlap_element(bra, ket) @@ -25,36 +28,46 @@ import FewBodyECG: _compute_overlap_element, _build_operator_matrix, _compute_ma end -spd(M) = 0.5 * (M + M') + (size(M, 1) == 1 ? 1.0e-12 : 0.0)I - @testset "_compute_overlap_element basic properties" begin A = [1.0;;] B = [1.0;;] - bra = Rank0Gaussian(A) - ket = Rank0Gaussian(B) + + s₁ = randn(2) + s₂ = randn(2) + + bra = Rank0Gaussian(A, s₁) + ket = Rank0Gaussian(B, s₂) val = _compute_overlap_element(bra, ket) @test isapprox(val, (π / 2)^(3 / 2); atol = 1.0e-12) @test isapprox( - _compute_overlap_element(Rank0Gaussian(A), Rank0Gaussian(B)), - _compute_overlap_element(Rank0Gaussian(B), Rank0Gaussian(A)); + _compute_overlap_element(Rank0Gaussian(A, s₁), Rank0Gaussian(B, s₂)), + _compute_overlap_element(Rank0Gaussian(B, s₂), Rank0Gaussian(A, s₁)); atol = 1.0e-12 ) A2 = [1.0 0.2; 0.2 1.5] B2 = [0.9 0.1; 0.1 1.2] - val2 = _compute_overlap_element(Rank0Gaussian(A2), Rank0Gaussian(B2)) + + val2 = _compute_overlap_element(Rank0Gaussian(A2, s₁), Rank0Gaussian(B2, s₂)) @test val2 > 0 end +spd(M) = 0.5 * (M + M') + (size(M, 1) == 1 ? 1.0e-12 : 0.0)I + @testset "build_overlap_matrix structure & values" begin A = spd([1.0 0.2; 0.2 1.5]) B = spd([0.9 0.1; 0.1 1.2]) C = spd([1.3 0.0; 0.0 0.8]) - g1 = Rank0Gaussian(A) - g2 = Rank0Gaussian(B) - g3 = Rank0Gaussian(C) + + s₁ = randn(2) + s₂ = randn(2) + s₃ = randn(2) + + g1 = Rank0Gaussian(A, s₁) + g2 = Rank0Gaussian(B, s₂) + g3 = Rank0Gaussian(C, s₃) basis = BasisSet([g1, g2, g3]) S = build_overlap_matrix(basis) @@ -74,7 +87,12 @@ end A = spd([1.0 0.2; 0.2 1.5]) B = spd([0.9 0.1; 0.1 1.2]) C = spd([1.3 0.0; 0.0 0.8]) - g = [Rank0Gaussian(A), Rank0Gaussian(B), Rank0Gaussian(C)] + + s₁ = randn(2) + s₂ = randn(2) + s₃ = randn(2) + + g = [Rank0Gaussian(A, s₁), Rank0Gaussian(B, s₂), Rank0Gaussian(C, s₃)] basis = BasisSet(g) K = [2.0 0.1; 0.1 2.0] @@ -102,7 +120,13 @@ end A = spd([1.0 0.2; 0.2 1.5]) B = spd([0.9 0.1; 0.1 1.2]) C = spd([1.3 0.0; 0.0 0.8]) - g = [Rank0Gaussian(A), Rank0Gaussian(B), Rank0Gaussian(C)] + + s₁ = randn(2) + s₂ = randn(2) + s₃ = randn(2) + + + g = [Rank0Gaussian(A, s₁), Rank0Gaussian(B, s₂), Rank0Gaussian(C, s₃)] basis = BasisSet(g) kop = KineticOperator([2.0 0.0; 0.0 2.0]) diff --git a/test/test_matrix_elements.jl b/test/test_matrix_elements.jl index 1db5897..2a201e9 100644 --- a/test/test_matrix_elements.jl +++ b/test/test_matrix_elements.jl @@ -9,8 +9,11 @@ import FewBodyECG: _compute_matrix_element B = [0.9 0.1; 0.1 1.2] K = rand(2, 2) - bra = Rank0Gaussian(A) - ket = Rank0Gaussian(B) + s₁ = randn(2) + s₂ = randn(2) + + bra = Rank0Gaussian(A, s₁) + ket = Rank0Gaussian(B, s₂) op = KineticOperator(K) result = _compute_matrix_element(bra, ket, op) @@ -27,8 +30,11 @@ end w = [1.0, -1.0] coefficient = 1.5 - bra = Rank0Gaussian(A) - ket = Rank0Gaussian(B) + s₁ = randn(2) + s₂ = randn(2) + + bra = Rank0Gaussian(A, s₁) + ket = Rank0Gaussian(B, s₂) op = CoulombOperator(coefficient, w) result = _compute_matrix_element(bra, ket, op) From 6f1447657078d0eed5737ce7e786f48fc57fcaac Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Thu, 30 Oct 2025 09:51:05 +0100 Subject: [PATCH 03/15] updated structs --- Examples/Positronium.jl | 5 +---- Project.toml | 4 +++- src/FewBodyECG.jl | 2 +- src/coordinates.jl | 47 +++++++++++++++-------------------------- src/matrix_elements.jl | 15 ++++++++++--- src/types.jl | 4 ++-- 6 files changed, 36 insertions(+), 41 deletions(-) diff --git a/Examples/Positronium.jl b/Examples/Positronium.jl index 30c4033..f434b84 100644 --- a/Examples/Positronium.jl +++ b/Examples/Positronium.jl @@ -6,16 +6,13 @@ using QuasiMonteCarlo masses = [1.0, 1.0, 1.0] psys = ParticleSystem(masses) -K = Diagonal([1 / 2, 1 / 2, 1 / 2]) -K_transformed = psys.J * K * psys.J' - w_list = [[1, -1, 0], [1, 0, -1], [0, 1, -1]] w_raw = [psys.U' * w for w in w_list] coeffs = [+1.0, -1.0, -1.0] ops = Operator[ - KineticOperator(K_transformed); + KineticOperator(masses); (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))... ] diff --git a/Project.toml b/Project.toml index 85cf46a..6c2f3b2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,11 +1,12 @@ name = "FewBodyECG" uuid = "083b1810-24a1-4a79-9a41-145bb2bb8ceb" -authors = ["Shuhei Ohno", "Martin Mikkelsen"] version = "1.0.5" +authors = ["Shuhei Ohno", "Martin Mikkelsen"] [deps] FewBodyHamiltonians = "3a126c26-e5d7-4a95-83c3-3b69f8a11ded" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" @@ -13,6 +14,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Aqua = "0.8.13" FewBodyHamiltonians = "0.0.2" LinearAlgebra = "1.7.3" +Plots = "1.41.1" QuasiMonteCarlo = "0.3.3" SpecialFunctions = "2.5.0" Test = "1.11.0" diff --git a/src/FewBodyECG.jl b/src/FewBodyECG.jl index c537865..6ab05a1 100644 --- a/src/FewBodyECG.jl +++ b/src/FewBodyECG.jl @@ -3,7 +3,7 @@ module FewBodyECG using LinearAlgebra using FewBodyHamiltonians -export ParticleSystem, default_b0 +export ParticleSystem, default_b0, Λ export generate_bij, _generate_A_matrix diff --git a/src/coordinates.jl b/src/coordinates.jl index cdbfa3c..c1040cd 100644 --- a/src/coordinates.jl +++ b/src/coordinates.jl @@ -1,18 +1,4 @@ -""" - struct ParticleSystem - -A structure representing a system of particles with associated masses and coordinate transformations. - -# Fields -- `masses::Vector{Float64}`: A vector containing the masses of the particles in the system. Must contain at least two elements. -- `J::Matrix{Float64}`: The Jacobi transformation matrix for the particle system, computed based on the masses. -- `U::Matrix{Float64}`: An auxiliary transformation matrix for the particle system, computed based on the masses. -- `scale::Union{Symbol,Nothing}`: An optional symbol indicating the scale of the system (e.g., `:atomic`, `:molecular`, `:nuclear`). Defaults to `nothing`. -# Constructor -- `ParticleSystem(masses::Vector{Float64}; scale::Union{Symbol,Nothing}=nothing)`: - Creates a new `ParticleSystem` instance. The `masses` vector must contain at least two elements. The `scale` parameter is optional and can be used to specify the scale of the system. The Jacobi and auxiliary transformation matrices (`J` and `U`) are computed internally using the `_jacobi_transform` function. -""" struct ParticleSystem masses::Vector{Float64} J::Matrix{Float64} @@ -50,6 +36,17 @@ function _jacobi_transform(masses::Vector{Float64})::Tuple{Matrix{Float64}, Matr return J, U end +function Λ(masses::Vector{<:Number}) + J, _ = _jacobi_transform(masses) + Minv = Diagonal(0.5 ./ masses) + Λ = Symmetric(J * Minv * J') + return Λ +end + +function KineticOperator(masses::Vector{<:Number}) + return KineticOperator(Λ(masses)) +end + """ default_b0(scale::Union{Symbol,Nothing}) -> Float64 @@ -70,14 +67,13 @@ function default_b0(scale::Union{Symbol, Nothing}) error("Unknown scale: $scale") end -function _generate_A_matrix(bij::Vector{Float64}, w_list::Vector{Vector{Float64}})::Matrix{Float64} - @assert length(bij) == length(w_list) "Length of `bij` and `w_list` must be equal." - dim = length(w_list[1]) +function _generate_A_matrix(bij::AbstractVector{<:Number}, w_list::AbstractVector{<:AbstractVector{<:Number}})::Matrix{Float64} + bijf = Float64.(bij) + w_listf = [Float64.(w) for w in w_list] + dim = length(w_listf[1]) A = zeros(Float64, dim, dim) - for i in 1:length(bij) - w = w_list[i] - @assert length(w) == dim "All weight vectors must have the same dimension." - A += (w * w') / (bij[i]^2) + for i in 1:length(bijf) + A .+= (w_listf[i] * w_listf[i]') / (bijf[i]^2) end return A end @@ -97,15 +93,6 @@ function _shift_vectors(a::Matrix{Float64}, b::Matrix{Float64}, mat::Union{Nothi return sum_val end -function _generate_weight_vector(dim::Int, i::Int, j::Int)::Vector{Int} - @assert 1 ≤ i ≤ dim "Index `i` must be between 1 and $dim." - @assert 1 ≤ j ≤ dim "Index `j` must be between 1 and $dim." - w = zeros(Int, dim) - w[i] = 1 - w[j] = -1 - return w -end - function _transform_coordinates(J::Matrix{Float64}, r::Vector{Float64})::Vector{Float64} @assert size(J, 2) == length(r) "Matrix `J` columns must match length of vector `r`." return J * r diff --git a/src/matrix_elements.jl b/src/matrix_elements.jl index b5d5689..dd7b477 100644 --- a/src/matrix_elements.jl +++ b/src/matrix_elements.jl @@ -6,12 +6,21 @@ compute_matrix_element(bra, ket, op) Compute the matrix element ⟨bra|op|ket⟩ using analytic expressions. """ +function _compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian) + B = bra.A + ket.A + v = bra.s + ket.s + N = size(B, 1) + return exp(0.25 * transpose(v) * inv(B) * v) * ((π^N) / det(B))^(3 / 2) +end + function _compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian, op::KineticOperator) A, B = bra.A, ket.A - K = op.K + v = bra.s + ket.s + K = op.K R = inv(A + B) - M0 = (π^length(R) / det(A + B))^(3 / 2) - return 6 * tr(B * K * A * R) * M0 + M0 = _compute_matrix_element(bra, ket) + u = 0.5 * inv(B) * v + return (6 * tr(B * K * A * R) + transpose(bra.s - 2 * bra.A * u) * K * (s - 2 * ket.A * u)) * M0 end function _compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian, op::CoulombOperator) diff --git a/src/types.jl b/src/types.jl index e848bcc..122e96f 100644 --- a/src/types.jl +++ b/src/types.jl @@ -22,8 +22,8 @@ struct BasisSet functions::Vector{GaussianBase} end -struct KineticOperator <: FewBodyHamiltonians.KineticTerm - K::Matrix{Number} +struct KineticOperator{T<:Number} <: FewBodyHamiltonians.KineticTerm + K::AbstractMatrix{T} end struct CoulombOperator <: FewBodyHamiltonians.PotentialTerm From c56becb8ce31adfa2e6e177e895d9fc17f65148c Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Fri, 31 Oct 2025 14:56:20 +0100 Subject: [PATCH 04/15] updated structs and solver --- Examples/HydrogenAnion.jl | 17 +++--- Examples/Positronium.jl | 74 ++++---------------------- src/FewBodyECG.jl | 2 +- src/coordinates.jl | 66 ++--------------------- src/hamiltonian.jl | 109 +++++++++++++++++++++----------------- src/matrix_elements.jl | 43 +++++++++------ src/sampling.jl | 43 ++++++++++----- src/types.jl | 61 ++++++++++++++------- 8 files changed, 184 insertions(+), 231 deletions(-) diff --git a/Examples/HydrogenAnion.jl b/Examples/HydrogenAnion.jl index 0d0240f..7a7c08a 100644 --- a/Examples/HydrogenAnion.jl +++ b/Examples/HydrogenAnion.jl @@ -4,27 +4,26 @@ using Plots using QuasiMonteCarlo masses = [1.0e15, 1.0, 1.0] -psys = ParticleSystem(masses) -K = Diagonal([0.0, 1 / 2, 1 / 2]) -K_transformed = psys.J * K * psys.J' +Λmat = Λ(masses) +kin = KineticOperator(Λmat) +J, U = _jacobi_transform(masses) w_list = [[1, -1, 0], [1, 0, -1], [0, 1, -1]] -w_raw = [psys.U' * w for w in w_list] +w_raw = [U' * w for w in w_list] coeffs = [-1.0, -1.0, +1.0] - ops = Operator[ - KineticOperator(K_transformed); + kin; (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))... ] -A = solve_ECG(ops, psys, 100) +result = solve_ECG(ops, masses, 100) E = -0.527751016523 -ΔE = abs(A.ground_state - E) +ΔE = abs(result.ground_state - E) @info "Energy difference" ΔE -n, E = convergence(A) +n, E = convergence(result) plot(n, E) \ No newline at end of file diff --git a/Examples/Positronium.jl b/Examples/Positronium.jl index f434b84..0337a53 100644 --- a/Examples/Positronium.jl +++ b/Examples/Positronium.jl @@ -1,70 +1,18 @@ -using FewBodyECG -using LinearAlgebra -using Plots -using QuasiMonteCarlo +using FewBodyECG, LinearAlgebra masses = [1.0, 1.0, 1.0] -psys = ParticleSystem(masses) -w_list = [[1, -1, 0], [1, 0, -1], [0, 1, -1]] -w_raw = [psys.U' * w for w in w_list] - -coeffs = [+1.0, -1.0, -1.0] - -ops = Operator[ - KineticOperator(masses); - (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))... -] - -A = solve_ECG(ops, psys, 100) - -let - n_basis = 100 - b1 = default_b0(psys.scale) - method = :quasirandom - basis_fns = GaussianBase[] - E₀_list = Float64[] - - for i in 1:n_basis - bij = generate_bij(method, i, length(w_raw), b1; qmc_sampler = SobolSample()) - A = _generate_A_matrix(bij, w_raw) - push!(basis_fns, Rank0Gaussian(A)) - - basis = BasisSet(basis_fns) - ops = Operator[ - KineticOperator(K_transformed); - (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))... - ] - - H = build_hamiltonian_matrix(basis, ops) - S = build_overlap_matrix(basis) - vals, vecs = solve_generalized_eigenproblem(H, S) - global c₀ = vecs[:, 1] - E₀ = minimum(vals) - - push!(E₀_list, E₀) - @info "Step $i" E₀ = E₀ - end - - E₀ = minimum(E₀_list) - Eᵗʰ = -0.2620050702328 - ΔE = abs(E₀ - Eᵗʰ) - @info "Energy difference" ΔE - - r = range(0.01, 14.0, length = 400) - ρ_r = [rval^2 * abs2(ψ₀([rval, 0.0], c₀, basis_fns)) for rval in r] - - p1 = plot( - r, ρ_r, xlabel = "r (a.u.)", ylabel = "r²|ψ₀(r)|²", - lw = 2, label = "r²C(r)", title = "Electron-Positron Correlation Function" - ) +Λmat = Λ(masses) +kin = KineticOperator(Λmat) +J, U = _jacobi_transform(masses) +w_list = [[1, -1, 0], [1, 0, -1], [0, 1, -1]] +w_raw = [U' * w for w in w_list] - p2 = plot( - 1:n_basis, E₀_list, xlabel = "Number of Gaussians", ylabel = "E₀ [Hartree]", - lw = 2, label = "Ground state energy", title = "Positronium Convergence" - ) +coeffs = [-1.0, -1.0, +1.0] +coulomb_ops = [CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw)] - plot(p1, p2, layout = (2, 1)) +ops = Operator[kin; coulomb_ops...] -end +result = solve_ECG(ops, masses, 150; scale=1) +println("E ≈ ", result.ground_state) \ No newline at end of file diff --git a/src/FewBodyECG.jl b/src/FewBodyECG.jl index 6ab05a1..c596b34 100644 --- a/src/FewBodyECG.jl +++ b/src/FewBodyECG.jl @@ -3,7 +3,7 @@ module FewBodyECG using LinearAlgebra using FewBodyHamiltonians -export ParticleSystem, default_b0, Λ +export Λ, _jacobi_transform export generate_bij, _generate_A_matrix diff --git a/src/coordinates.jl b/src/coordinates.jl index c1040cd..fcdb365 100644 --- a/src/coordinates.jl +++ b/src/coordinates.jl @@ -1,17 +1,3 @@ - -struct ParticleSystem - masses::Vector{Float64} - J::Matrix{Float64} - U::Matrix{Float64} - scale::Union{Symbol, Nothing} # :atomic, :molecular, :nuclear, etc. - - function ParticleSystem(masses::Vector{Float64}; scale::Union{Symbol, Nothing} = nothing) - @assert length(masses) ≥ 2 "At least two masses are required for a particle system." - J, U = _jacobi_transform(masses) - return new(masses, J, U, scale) - end -end - function _jacobi_transform(masses::Vector{Float64})::Tuple{Matrix{Float64}, Matrix{Float64}} N = length(masses) @assert N ≥ 2 "At least two masses are required for Jacobi transformation." @@ -36,63 +22,17 @@ function _jacobi_transform(masses::Vector{Float64})::Tuple{Matrix{Float64}, Matr return J, U end -function Λ(masses::Vector{<:Number}) +function Λ(masses::Vector{<:Real}) J, _ = _jacobi_transform(masses) - Minv = Diagonal(0.5 ./ masses) + Minv = Diagonal(1.9 ./ masses) Λ = Symmetric(J * Minv * J') return Λ end -function KineticOperator(masses::Vector{<:Number}) +function KineticOperator(masses::Vector{<:Real}) return KineticOperator(Λ(masses)) end -""" - default_b0(scale::Union{Symbol,Nothing}) -> Float64 - -Returns a default value for the parameter `b0` based on the provided `scale`. - -# Arguments -- `scale::Union{Symbol,Nothing}`: A symbol representing the scale type or `nothing`. - - `:atomic`: Returns `1.0`, corresponding to the Bohr radius in atomic units. - - `:molecular`: Returns `3.0`, representing a typical molecular bond length. - - `:nuclear`: Returns `0.03`, approximately 1 femtometer in atomic units. - - `nothing`: Returns `1.0` as a fallback default. -""" -function default_b0(scale::Union{Symbol, Nothing}) - scale === :atomic && return 1.0 - scale === :molecular && return 3.0 - scale === :nuclear && return 0.03 - scale === nothing && return 10.0 - error("Unknown scale: $scale") -end - -function _generate_A_matrix(bij::AbstractVector{<:Number}, w_list::AbstractVector{<:AbstractVector{<:Number}})::Matrix{Float64} - bijf = Float64.(bij) - w_listf = [Float64.(w) for w in w_list] - dim = length(w_listf[1]) - A = zeros(Float64, dim, dim) - for i in 1:length(bijf) - A .+= (w_listf[i] * w_listf[i]') / (bijf[i]^2) - end - return A -end - -function _shift_vectors(a::Matrix{Float64}, b::Matrix{Float64}, mat::Union{Nothing, Matrix{Float64}} = nothing)::Float64 - n = size(a, 2) - @assert n == size(b, 2) "Matrices `a` and `b` must have the same number of columns." - mat = mat === nothing ? I(n) : mat - @assert size(mat) == (n, n) "Matrix `mat` must be square with size equal to number of vectors." - - sum_val = 0.0 - for i in 1:n - for j in 1:n - sum_val += mat[i, j] * dot(view(a, :, i), view(b, :, j)) - end - end - return sum_val -end - function _transform_coordinates(J::Matrix{Float64}, r::Vector{Float64})::Vector{Float64} @assert size(J, 2) == length(r) "Matrix `J` columns must match length of vector `r`." return J * r diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl index b8695ac..47692d1 100644 --- a/src/hamiltonian.jl +++ b/src/hamiltonian.jl @@ -1,81 +1,92 @@ using FewBodyHamiltonians -using QuasiMonteCarlo +using LinearAlgebra -function _compute_overlap_element(bra::GaussianBase, ket::GaussianBase) - A, B = bra.A, ket.A - R = inv(A + B) - return (π^length(R) / det(A + B))^(3 / 2) +function _compute_overlap_element(bra::Rank0Gaussian, ket::Rank0Gaussian) + _compute_matrix_element(bra, ket) end -function build_overlap_matrix(basis::BasisSet) +function build_overlap_matrix(basis::BasisSet{<:GaussianBase}) n = length(basis.functions) - S = zeros(n, n) + S = Matrix{Float64}(undef, n, n) for i in 1:n, j in 1:i - val = _compute_overlap_element(basis.functions[i], basis.functions[j]) - S[i, j] = S[j, i] = val + val = _compute_overlap_element(basis.functions[i]::Rank0Gaussian, basis.functions[j]::Rank0Gaussian) + S[i, j] = val + S[j, i] = val end - return S + S end -function _build_operator_matrix(basis::BasisSet, op::FewBodyHamiltonians.Operator) +function _build_operator_matrix(basis::BasisSet{<:GaussianBase}, op::FewBodyHamiltonians.Operator) n = length(basis.functions) - H = zeros(n, n) + H = Matrix{Float64}(undef, n, n) for i in 1:n, j in 1:i - val = _compute_matrix_element(basis.functions[i], basis.functions[j], op) - H[i, j] = H[j, i] = val + val = _compute_matrix_element(basis.functions[i]::Rank0Gaussian, basis.functions[j]::Rank0Gaussian, op) + H[i, j] = val + H[j, i] = val end - return H + H end -function build_hamiltonian_matrix(basis::BasisSet, operators::AbstractVector{<:FewBodyHamiltonians.Operator}) - H = zeros(length(basis.functions), length(basis.functions)) +function build_hamiltonian_matrix(basis::BasisSet{<:GaussianBase}, operators::AbstractVector{<:FewBodyHamiltonians.Operator}) + n = length(basis.functions) + H = zeros(Float64, n, n) for op in operators H .+= _build_operator_matrix(basis, op) end - return H + H end -function solve_generalized_eigenproblem(H::Matrix{Float64}, S::Matrix{Float64}) - try - F = cholesky(S; check = true) - L = F.L - Linv = inv(L) - A = Linv * H * Linv' - vals, vecs = eigen(Symmetric(A)) - vecs_orig = Linv' * vecs - return real(vals), real(vecs_orig) - catch err - @warn "Cholesky on S failed, falling back to generalized eigen solver" exception = (err, catch_backtrace()) - vals, vecs = eigen(H, S) - return real(vals), real(vecs) - end +function solve_generalized_eigenproblem(H::AbstractMatrix{<:Real}, S::AbstractMatrix{<:Real}) + F = cholesky(Symmetric(S); check = true) + L = F.L + A = (L \ H) / L' + evals, evecs = eigen(Symmetric(A)) + vecs = L' \ evecs + real(evals), real(vecs) end -function solve_ECG(operators::Vector{FewBodyHamiltonians.Operator}, system::ParticleSystem, n::Int = 50; sampler = SobolSample(), method::Symbol = :quasirandom, verbose::Bool = true) - b₁ = default_b0(system.scale) - basis_fns = GaussianBase[] - E₀ = Float64[] - coulomb_length = count(x -> x isa CoulombOperator, operators) +function diagonalize(basis::BasisSet{<:GaussianBase}, operators::AbstractVector{<:FewBodyHamiltonians.Operator}) + H = build_hamiltonian_matrix(basis, operators) + S = build_overlap_matrix(basis) + solve_generalized_eigenproblem(H, S) +end + +function solve_ECG(operators::Vector{<:FewBodyHamiltonians.Operator}, + n::Int=50; + sampler=HaltonSample(), + method::Symbol=:quasirandom, + scale::Real=0.2, + sscale::Real=0.1, + verbose::Bool=true) + + b₁ = float(scale) + basis_fns = Rank0Gaussian[] + E_hist = Float64[] + vecs_list = Any[] + w_list = [op.w for op in operators if op isa CoulombOperator] - E₀_list = Float64[] - vecs_list = [] - for i in 1:n - bij = generate_bij(method, i, coulomb_length, b₁; qmc_sampler = sampler) - A = _generate_A_matrix(bij, w_list) - push!(basis_fns, Rank0Gaussian(A)) + n_pairs = length(w_list) + d = length(w_list[1]) - basis = BasisSet(basis_fns) + for i in 1:n + bij = generate_bij(method, i, n_pairs, b₁; qmc_sampler=sampler) + A, _ = _generate_A_matrix(bij, w_list) + s = generate_shift(method, i, d, sscale; qmc_sampler=sampler) # ▼ use nonzero shifts + push!(basis_fns, Rank0Gaussian(A, s)) + basis = BasisSet{Rank0Gaussian}(basis_fns) H = build_hamiltonian_matrix(basis, operators) S = build_overlap_matrix(basis) λs, Us = solve_generalized_eigenproblem(H, S) - E₀ = minimum(λs) + E0 = minimum(λs) - push!(E₀_list, E₀) + push!(E_hist, E0) push!(vecs_list, Us) - verbose && @info "Step $i" E₀ = E₀ + verbose && @info "Step $i" E₀=E0 end - @info "Minimum found" E₀ - return SolverResults(basis_fns, n, operators, :quasirandom, sampler, b₁, E₀, E₀_list, vecs_list) + + Emin = last(E_hist) + @info "Minimum found" E₀=Emin + return SolverResults(basis_fns, n, operators, method, sampler, b₁, Emin, E_hist, vecs_list) end diff --git a/src/matrix_elements.jl b/src/matrix_elements.jl index dd7b477..808c3db 100644 --- a/src/matrix_elements.jl +++ b/src/matrix_elements.jl @@ -7,28 +7,41 @@ Compute the matrix element ⟨bra|op|ket⟩ using analytic expressions. """ function _compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian) - B = bra.A + ket.A - v = bra.s + ket.s - N = size(B, 1) - return exp(0.25 * transpose(v) * inv(B) * v) * ((π^N) / det(B))^(3 / 2) + A, B = bra.A, ket.A + a, b = bra.s, ket.s + S = A + B + R = inv(S) + n = size(S, 1) + M0 = (π^n / det(S))^(3/2) + return exp(0.25 * (a + b)' * R * (a + b)) * M0 end function _compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian, op::KineticOperator) A, B = bra.A, ket.A - v = bra.s + ket.s - K = op.K - R = inv(A + B) - M0 = _compute_matrix_element(bra, ket) - u = 0.5 * inv(B) * v - return (6 * tr(B * K * A * R) + transpose(bra.s - 2 * bra.A * u) * K * (s - 2 * ket.A * u)) * M0 + a, b = bra.s, ket.s + K = op.K + S = A + B + R = inv(S) + M = _compute_matrix_element(bra, ket) + term = 6 * tr(B * K * A * R) + + b' * K * a + + (a + b)' * R * B * K * A * R * (a + b) - + (a + b)' * R * B * K * a - + b' * K * A * R * (a + b) + return term * M end function _compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian, op::CoulombOperator) - A, B, w = bra.A, ket.A, op.w - R = inv(A + B) - β = 1 / (dot(w, R * w)) - M0 = (π^length(R) / det(A + B))^(3 / 2) - return op.coefficient * 2 * sqrt(β / π) * M0 + A, B = bra.A, ket.A + a, b = bra.s, ket.s + w = op.w + S = A + B + R = inv(S) + M = _compute_matrix_element(bra, ket) + β = 1 / (w' * R * w) + q = 0.5 * (w' * R * (a + b)) + f = abs(q) < 1e-12 ? (2 * sqrt(β / π)) : (erf(sqrt(β) * q) / q) + return op.coefficient * f * M end function _compute_matrix_element(bra::Rank1Gaussian, ket::Rank1Gaussian, op::CoulombOperator) diff --git a/src/sampling.jl b/src/sampling.jl index 54a3f06..9281fba 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -1,18 +1,37 @@ using QuasiMonteCarlo -export generate_bij +export generate_bij, generate_shift, _generate_A_matrix, build_rank0 + +function _qmc_point(i::Int, d::Int; sampler = HaltonSample()) + QuasiMonteCarlo.sample(i + 1, d, sampler)[:, end] +end + +function generate_bij(method::Symbol, i::Int, n_terms::Int, b1::Float64; + qmc_sampler=HaltonSample(), bmin=0.02*b1, bmax=20b1) + u = method === :quasirandom ? QuasiMonteCarlo.sample(i + 1, n_terms, qmc_sampler)[:, end] : + method === :random ? rand(n_terms) : + error("Unsupported method $method") + bmin .* (bmax/bmin) .^ u +end + +function generate_shift(method::Symbol, i::Int, dim::Int, sscale::Real; + qmc_sampler=HaltonSample()) + u = method === :quasirandom ? QuasiMonteCarlo.sample(i + 1, dim, qmc_sampler)[:, end] : + method === :random ? rand(dim) : + error("Unsupported method $method") + sscale .* (2u .- 1) +end -""" - generate_bij(method::Symbol, i::Int, n_terms::Int, b1::Float64) -> Vector{Float64} -Generate a bij vector using the specified sampling method. -""" -function generate_bij(method::Symbol, i::Int, n_terms::Int, b1::Float64; qmc_sampler = HaltonSample()) - if method == :quasirandom - return QuasiMonteCarlo.sample(i + 1, n_terms, qmc_sampler)[:, end] * b1 - elseif method == :random - return rand(n_terms) * b1 - else - error("Unsupported sampling method: $method") +function _generate_A_matrix(bij::AbstractVector{<:Real}, w_list::AbstractVector{<:AbstractVector{<:Real}}) + b = Float64.(bij) + W = [Float64.(w) for w in w_list] + length(b) == length(W) || throw(ArgumentError("Length of bij must equal number of w vectors")) + d = length(W[1]) + A = zeros(Float64, d, d) + for k in eachindex(b) + A .+= (W[k] * W[k]') / (b[k]^2) end + s = zeros(Float64, d) + A, s end diff --git a/src/types.jl b/src/types.jl index 122e96f..4373a3c 100644 --- a/src/types.jl +++ b/src/types.jl @@ -2,36 +2,59 @@ using FewBodyHamiltonians abstract type GaussianBase end -struct Rank0Gaussian <: GaussianBase - A::Matrix{Number} - s::Vector{Number} +struct Rank0Gaussian{T<:Real, M<:AbstractMatrix{T}, V<:AbstractVector{T}} <: GaussianBase + A::Symmetric{T,M} + s::V + function Rank0Gaussian(A::AbstractMatrix{T}, s::AbstractVector{T}) where {T<:Real} + size(A,1) == size(A,2) || throw(ArgumentError("A must be square")) + length(s) == size(A,1) || throw(ArgumentError("length(s) != size(A,1)")) + new{T, typeof(A), typeof(s)}(Symmetric(A), s) + end end -struct Rank1Gaussian <: GaussianBase - A::Matrix{Number} - a::Vector{Number} +struct Rank1Gaussian{T<:Real, M<:AbstractMatrix{T}, V<:AbstractVector{T}} <: GaussianBase + A::Symmetric{T,M} + a::V + s::V + function Rank1Gaussian(A::AbstractMatrix{T}, a::AbstractVector{T}, s::AbstractVector{T}) where {T<:Real} + size(A,1) == size(A,2) || throw(ArgumentError("A must be square")) + (length(a) == size(A,1) && length(s) == size(A,1)) || + throw(ArgumentError("length(a) and length(s) must equal size(A,1)")) + new{T, typeof(A), typeof(a)}(Symmetric(A), a, s) + end end -struct Rank2Gaussian <: GaussianBase - A::Matrix{Number} - a::Vector{Number} - b::Vector{Number} +struct Rank2Gaussian{T<:Real, M<:AbstractMatrix{T}, V<:AbstractVector{T}} <: GaussianBase + A::Symmetric{T,M} + a::V + b::V + s::V + function Rank2Gaussian(A::AbstractMatrix{T}, a::AbstractVector{T}, b::AbstractVector{T}, s::AbstractVector{T}) where {T<:Real} + size(A,1) == size(A,2) || throw(ArgumentError("A must be square")) + (length(a) == size(A,1) && length(b) == size(A,1) && length(s) == size(A,1)) || + throw(ArgumentError("length(a), length(b), length(s) must equal size(A,1)")) + new{T, typeof(A), typeof(a)}(Symmetric(A), a, b, s) + end end -struct BasisSet - functions::Vector{GaussianBase} +struct BasisSet{G<:GaussianBase} + functions::Vector{G} end -struct KineticOperator{T<:Number} <: FewBodyHamiltonians.KineticTerm +struct KineticOperator{T<:Real} <: FewBodyHamiltonians.KineticTerm K::AbstractMatrix{T} end -struct CoulombOperator <: FewBodyHamiltonians.PotentialTerm - coefficient::Float64 - w::Vector{Number} +struct CoulombOperator{T<:Real} <: FewBodyHamiltonians.PotentialTerm + coefficient::T + w::AbstractVector{T} end -struct ECG - basis::BasisSet - operators::Vector{Operator} +struct ECG{G<:GaussianBase,O} + basis::BasisSet{G} + operators::Vector{O} end + +validate!(g::Rank0Gaussian) = (cholesky(g.A); g) +validate!(g::Rank1Gaussian) = (cholesky(g.A); g) +validate!(g::Rank2Gaussian) = (cholesky(g.A); g) \ No newline at end of file From 3accb6e8eff7326003f2595ba7b71f21eb1ebb0d Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Tue, 4 Nov 2025 14:21:37 +0100 Subject: [PATCH 05/15] updated sampling and shifts --- Examples/HydrogenAnion.jl | 2 +- Examples/Positronium.jl | 4 ++-- src/coordinates.jl | 2 +- src/hamiltonian.jl | 4 ++-- src/sampling.jl | 50 ++++++++++++++++++++++++++------------- 5 files changed, 40 insertions(+), 22 deletions(-) diff --git a/Examples/HydrogenAnion.jl b/Examples/HydrogenAnion.jl index 7a7c08a..813950f 100644 --- a/Examples/HydrogenAnion.jl +++ b/Examples/HydrogenAnion.jl @@ -19,7 +19,7 @@ ops = Operator[ (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))... ] -result = solve_ECG(ops, masses, 100) +result = solve_ECG(ops, 500, scale=1.5) E = -0.527751016523 ΔE = abs(result.ground_state - E) diff --git a/Examples/Positronium.jl b/Examples/Positronium.jl index 0337a53..4393614 100644 --- a/Examples/Positronium.jl +++ b/Examples/Positronium.jl @@ -9,10 +9,10 @@ J, U = _jacobi_transform(masses) w_list = [[1, -1, 0], [1, 0, -1], [0, 1, -1]] w_raw = [U' * w for w in w_list] -coeffs = [-1.0, -1.0, +1.0] +coeffs = [+1.0, -1.0, -1.0] coulomb_ops = [CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw)] ops = Operator[kin; coulomb_ops...] -result = solve_ECG(ops, masses, 150; scale=1) +result = solve_ECG(ops, 150; scale=10.0) println("E ≈ ", result.ground_state) \ No newline at end of file diff --git a/src/coordinates.jl b/src/coordinates.jl index fcdb365..4616c71 100644 --- a/src/coordinates.jl +++ b/src/coordinates.jl @@ -24,7 +24,7 @@ end function Λ(masses::Vector{<:Real}) J, _ = _jacobi_transform(masses) - Minv = Diagonal(1.9 ./ masses) + Minv = Diagonal(0.5 ./ masses) Λ = Symmetric(J * Minv * J') return Λ end diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl index 47692d1..8fa8a45 100644 --- a/src/hamiltonian.jl +++ b/src/hamiltonian.jl @@ -70,8 +70,8 @@ function solve_ECG(operators::Vector{<:FewBodyHamiltonians.Operator}, for i in 1:n bij = generate_bij(method, i, n_pairs, b₁; qmc_sampler=sampler) - A, _ = _generate_A_matrix(bij, w_list) - s = generate_shift(method, i, d, sscale; qmc_sampler=sampler) # ▼ use nonzero shifts + A = _generate_A_matrix(bij, w_list) + s = generate_shift(method, i, length(w_list[1]), sscale; qmc_sampler=sampler) push!(basis_fns, Rank0Gaussian(A, s)) basis = BasisSet{Rank0Gaussian}(basis_fns) diff --git a/src/sampling.jl b/src/sampling.jl index 9281fba..f69f472 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -7,31 +7,49 @@ function _qmc_point(i::Int, d::Int; sampler = HaltonSample()) end function generate_bij(method::Symbol, i::Int, n_terms::Int, b1::Float64; - qmc_sampler=HaltonSample(), bmin=0.02*b1, bmax=20b1) - u = method === :quasirandom ? QuasiMonteCarlo.sample(i + 1, n_terms, qmc_sampler)[:, end] : + qmc_sampler = HaltonSample(), + bmin = 0.02*b1, bmax = 20b1) + bmin > 0 || throw(ArgumentError("bmin must be > 0")) + bmax > bmin || throw(ArgumentError("bmax must be > bmin")) + u = method === :quasirandom ? _qmc_point(i, n_terms; sampler=qmc_sampler) : method === :random ? rand(n_terms) : error("Unsupported method $method") - bmin .* (bmax/bmin) .^ u + bmin .* (bmax / bmin) .^ u end +# symmetric shifts in [-sscale, +sscale] function generate_shift(method::Symbol, i::Int, dim::Int, sscale::Real; - qmc_sampler=HaltonSample()) - u = method === :quasirandom ? QuasiMonteCarlo.sample(i + 1, dim, qmc_sampler)[:, end] : + qmc_sampler = HaltonSample()) + u = method === :quasirandom ? _qmc_point(i, dim; sampler=qmc_sampler) : method === :random ? rand(dim) : error("Unsupported method $method") - sscale .* (2u .- 1) + sscale .* (2u .- 1) end - -function _generate_A_matrix(bij::AbstractVector{<:Real}, w_list::AbstractVector{<:AbstractVector{<:Real}}) +# A = W * Diag(1/b.^2) * W' +function _generate_A_matrix(bij::AbstractVector{<:Real}, + w_list::AbstractVector{<:AbstractVector{<:Real}}) b = Float64.(bij) - W = [Float64.(w) for w in w_list] - length(b) == length(W) || throw(ArgumentError("Length of bij must equal number of w vectors")) - d = length(W[1]) - A = zeros(Float64, d, d) - for k in eachindex(b) - A .+= (W[k] * W[k]') / (b[k]^2) + m = length(b) + m == length(w_list) || throw(ArgumentError("Length(bij) must equal number of w vectors")) + + d = length(w_list[1]) + @inbounds for k in 2:m + length(w_list[k]) == d || throw(ArgumentError("All w vectors must have same length")) + end + + W = Matrix{Float64}(undef, d, m) + @inbounds for k in 1:m + @views W[:, k] = Float64.(w_list[k]) end - s = zeros(Float64, d) - A, s + + A = W * Diagonal(1.0 ./ (b .^ 2)) * W' + A end + +function build_rank0(bij::AbstractVector{<:Real}, + w_list::AbstractVector{<:AbstractVector{<:Real}}, + s::AbstractVector{<:Real}) + A = _generate_A_matrix(bij, w_list) + Rank0Gaussian(A, Float64.(s)) +end \ No newline at end of file From 904cbbb96e05cd7b042ea42a8cb8d2b52527ae14 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Thu, 6 Nov 2025 08:32:20 +0100 Subject: [PATCH 06/15] updated examples --- Examples/HydrogenAnion.jl | 3 ++- Examples/Positronium.jl | 5 +++-- "Examples/td\316\274.jl" | 19 +++++++++++++++++++ 3 files changed, 24 insertions(+), 3 deletions(-) create mode 100644 "Examples/td\316\274.jl" diff --git a/Examples/HydrogenAnion.jl b/Examples/HydrogenAnion.jl index 813950f..0f9d4e8 100644 --- a/Examples/HydrogenAnion.jl +++ b/Examples/HydrogenAnion.jl @@ -26,4 +26,5 @@ E = -0.527751016523 @info "Energy difference" ΔE n, E = convergence(result) -plot(n, E) \ No newline at end of file +plot(n, E) + diff --git a/Examples/Positronium.jl b/Examples/Positronium.jl index 4393614..3dac99e 100644 --- a/Examples/Positronium.jl +++ b/Examples/Positronium.jl @@ -14,5 +14,6 @@ coulomb_ops = [CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw)] ops = Operator[kin; coulomb_ops...] -result = solve_ECG(ops, 150; scale=10.0) -println("E ≈ ", result.ground_state) \ No newline at end of file +result = solve_ECG(ops, 250,sampler=SobolSample(); scale=10.0) +println("E ≈ ", result.ground_state) + diff --git "a/Examples/td\316\274.jl" "b/Examples/td\316\274.jl" new file mode 100644 index 0000000..19eff8e --- /dev/null +++ "b/Examples/td\316\274.jl" @@ -0,0 +1,19 @@ +using FewBodyECG, LinearAlgebra + +masses = [5496.918, 3670.481, 206.7686] + +Λmat = Λ(masses) +kin = KineticOperator(Λmat) + +J, U = _jacobi_transform(masses) +base_w = [[1, -1, 0], [1, 0, -1], [0, 1, -1]] +coeffs = [+1.0, +1.0, -1.0] +w_list = [c .* w for (c, w) in zip(coeffs, base_w)] +w_raw = [U' * w for w in w_list] + +coulomb_ops = [CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw)] + +ops = Operator[kin; coulomb_ops...] + +result = solve_ECG(ops, 250, sampler=SobolSample(); scale=0.025) +println("E ≈ ", result.ground_state) From 8a478f89d139df5c8ea63e994d70ec367b80cb57 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Wed, 19 Nov 2025 10:18:11 +0100 Subject: [PATCH 07/15] added plot utilities and matrix elements --- Examples/HydrogenAnion.jl | 4 +- Examples/Hydrogen_p-wave.jl | 71 ------ Examples/Hydrogen_s-wave.jl | 59 ----- Examples/Positronium.jl | 7 +- "Examples/td\316\274.jl" | 6 +- Project.toml | 2 - src/coordinates.jl | 6 +- src/hamiltonian.jl | 100 +++++++-- src/matrix_elements.jl | 12 +- src/sampling.jl | 36 ++- src/types.jl | 46 ++-- src/utils.jl | 62 +++++- test/test_coordinates.jl | 149 ++++++++----- test/test_hamiltonian.jl | 178 ++------------- test/test_matrix_elements.jl | 411 ++++++++++++++++++++--------------- test/test_sampling.jl | 18 -- 16 files changed, 532 insertions(+), 635 deletions(-) delete mode 100644 Examples/Hydrogen_p-wave.jl delete mode 100644 Examples/Hydrogen_s-wave.jl diff --git a/Examples/HydrogenAnion.jl b/Examples/HydrogenAnion.jl index 0f9d4e8..77b6d98 100644 --- a/Examples/HydrogenAnion.jl +++ b/Examples/HydrogenAnion.jl @@ -2,6 +2,7 @@ using FewBodyECG using LinearAlgebra using Plots using QuasiMonteCarlo +using FewBodyDB masses = [1.0e15, 1.0, 1.0] @@ -19,7 +20,7 @@ ops = Operator[ (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))... ] -result = solve_ECG(ops, 500, scale=1.5) +result = solve_ECG(ops, 250, scale = 1.0) E = -0.527751016523 ΔE = abs(result.ground_state - E) @@ -27,4 +28,3 @@ E = -0.527751016523 n, E = convergence(result) plot(n, E) - diff --git a/Examples/Hydrogen_p-wave.jl b/Examples/Hydrogen_p-wave.jl deleted file mode 100644 index 006456f..0000000 --- a/Examples/Hydrogen_p-wave.jl +++ /dev/null @@ -1,71 +0,0 @@ -using FewBodyECG -using LinearAlgebra -using Plots -using QuasiMonteCarlo - -masses = [1.0e15, 1.0] -psys = ParticleSystem(masses) - -K = Diagonal([0.0, 0.5]) -K_transformed = psys.J * K * psys.J' - -w_raw = [psys.U' * [1, -1]] -coeffs = [-1.0] - -n_basis = 25 -method = :quasirandom -b1 = 1.4 - -basis_fns = GaussianBase[] -E₀_list = Float64[] - -a_vec = [1.0] - -ops = Operator[ - KineticOperator(K_transformed); - (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))... -] - -A = solve_ECG(ops, psys, 50) - -for i in 1:n_basis - bij = generate_bij(method, i, length(w_raw), b1; qmc_sampler = HaltonSample()) - A = _generate_A_matrix(bij, w_raw) - # pass a plain vector for the rank-1 displacement - push!(basis_fns, Rank1Gaussian(A, a_vec)) - - basis = BasisSet(basis_fns) - ops = Operator[ - KineticOperator(K_transformed); - (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))... - ] - - H = build_hamiltonian_matrix(basis, ops) - S = build_overlap_matrix(basis) - - vals, _ = solve_generalized_eigenproblem(H, S) - valid = vals .> 1.0e-12 - S⁻¹₂ = Diagonal(1 ./ sqrt.(vals[valid])) - H̃ = S⁻¹₂ * H[valid, valid] * S⁻¹₂ - eigvals = eigen(H̃).values - - real_eigvals = eigvals[abs.(imag.(eigvals)) .< 1.0e-10] - if !isempty(real_eigvals) - E₀ = minimum(real_eigvals) - else - E₀ = NaN - @warn "All eigenvalues are complex at step $i" - end - push!(E₀_list, E₀) - println("Step $i: E₀ = $E₀") -end - -E_exact = -0.125 -E_min = minimum(E₀_list) -@show ΔE = abs(E_min - E_exact) - -plot( - 1:n_basis, E₀_list, xlabel = "Number of Gaussians", ylabel = "E₀ [Hartree]", - lw = 2, label = "E₀ estimate", title = "p-wave Hydrogen Convergence" -) -hline!([E_exact], label = "Exact: -0.125", linestyle = :dash) diff --git a/Examples/Hydrogen_s-wave.jl b/Examples/Hydrogen_s-wave.jl deleted file mode 100644 index 003ea5a..0000000 --- a/Examples/Hydrogen_s-wave.jl +++ /dev/null @@ -1,59 +0,0 @@ -using FewBodyECG -using LinearAlgebra -using Plots -using QuasiMonteCarlo - -import FewBodyECG: _generate_A_matrix - -masses = [1.0e15, 1.0] # proton, electron -psys = ParticleSystem(masses) - -K = Diagonal([0.0, 0.5]) -K_transformed = psys.J * K * psys.J' - -w_raw = [psys.U' * [1, -1]] # r₁ - r₂ -coeffs = [-1.0] # Coulomb attraction - -n_basis = 25 -method = :quasirandom -b1 = 1.5 - -basis_fns = GaussianBase[] -E₀_list = Float64[] - -ops = Operator[ - KineticOperator(K_transformed); - (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))... -] - - -for i in 1:n_basis - bij = generate_bij(method, i, length(w_raw), b1; qmc_sampler = SobolSample()) - A = _generate_A_matrix(bij, w_raw) - push!(basis_fns, Rank0Gaussian(A)) - - basis = BasisSet(basis_fns) - - H = build_hamiltonian_matrix(basis, ops) - S = build_overlap_matrix(basis) - - λs, Us = eigen(S) - keep = λs .> 1.0e-10 - S⁻¹₂ = Us[:, keep] * Diagonal(1 ./ sqrt.(λs[keep])) * Us[:, keep]' - H̃ = Symmetric(S⁻¹₂ * H * S⁻¹₂) - E₀ = minimum(eigen(H̃).values) - - push!(E₀_list, E₀) - - -end - -E_exact = -0.5 -E_min = minimum(E₀_list) -@show ΔE = abs(E_min - E_exact) - -plot( - 1:n_basis, E₀_list, xlabel = "Number of Gaussians", ylabel = "E₀ [Hartree]", - lw = 2, label = "E₀ estimate", title = "s-wave Hydrogen Convergence" -) -hline!([E_exact], label = "Exact: -0.5", linestyle = :dash) diff --git a/Examples/Positronium.jl b/Examples/Positronium.jl index 3dac99e..1606b60 100644 --- a/Examples/Positronium.jl +++ b/Examples/Positronium.jl @@ -7,13 +7,12 @@ kin = KineticOperator(Λmat) J, U = _jacobi_transform(masses) w_list = [[1, -1, 0], [1, 0, -1], [0, 1, -1]] -w_raw = [U' * w for w in w_list] +w_raw = [U' * w for w in w_list] -coeffs = [+1.0, -1.0, -1.0] +coeffs = [+1.0, -1.0, -1.0] coulomb_ops = [CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw)] ops = Operator[kin; coulomb_ops...] -result = solve_ECG(ops, 250,sampler=SobolSample(); scale=10.0) +result = solve_ECG(ops, 250, sampler = SobolSample(); scale = 0.2) println("E ≈ ", result.ground_state) - diff --git "a/Examples/td\316\274.jl" "b/Examples/td\316\274.jl" index 19eff8e..3d9a4fe 100644 --- "a/Examples/td\316\274.jl" +++ "b/Examples/td\316\274.jl" @@ -9,11 +9,11 @@ J, U = _jacobi_transform(masses) base_w = [[1, -1, 0], [1, 0, -1], [0, 1, -1]] coeffs = [+1.0, +1.0, -1.0] w_list = [c .* w for (c, w) in zip(coeffs, base_w)] -w_raw = [U' * w for w in w_list] - +w_raw = [U' * w for w in w_list] + coulomb_ops = [CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw)] ops = Operator[kin; coulomb_ops...] -result = solve_ECG(ops, 250, sampler=SobolSample(); scale=0.025) +result = solve_ECG(ops, 250, sampler = HaltonSample(); scale = 0.025) println("E ≈ ", result.ground_state) diff --git a/Project.toml b/Project.toml index 6c2f3b2..a871bcc 100644 --- a/Project.toml +++ b/Project.toml @@ -6,7 +6,6 @@ authors = ["Shuhei Ohno", "Martin Mikkelsen"] [deps] FewBodyHamiltonians = "3a126c26-e5d7-4a95-83c3-3b69f8a11ded" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" @@ -14,7 +13,6 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Aqua = "0.8.13" FewBodyHamiltonians = "0.0.2" LinearAlgebra = "1.7.3" -Plots = "1.41.1" QuasiMonteCarlo = "0.3.3" SpecialFunctions = "2.5.0" Test = "1.11.0" diff --git a/src/coordinates.jl b/src/coordinates.jl index 4616c71..e4948e0 100644 --- a/src/coordinates.jl +++ b/src/coordinates.jl @@ -23,9 +23,9 @@ function _jacobi_transform(masses::Vector{Float64})::Tuple{Matrix{Float64}, Matr end function Λ(masses::Vector{<:Real}) - J, _ = _jacobi_transform(masses) - Minv = Diagonal(0.5 ./ masses) - Λ = Symmetric(J * Minv * J') + J, _ = _jacobi_transform(masses) + Minv = Diagonal(0.5 ./ masses) + Λ = Symmetric(J * Minv * J') return Λ end diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl index 8fa8a45..421086f 100644 --- a/src/hamiltonian.jl +++ b/src/hamiltonian.jl @@ -2,7 +2,7 @@ using FewBodyHamiltonians using LinearAlgebra function _compute_overlap_element(bra::Rank0Gaussian, ket::Rank0Gaussian) - _compute_matrix_element(bra, ket) + return _compute_matrix_element(bra, ket) end function build_overlap_matrix(basis::BasisSet{<:GaussianBase}) @@ -13,7 +13,7 @@ function build_overlap_matrix(basis::BasisSet{<:GaussianBase}) S[i, j] = val S[j, i] = val end - S + return S end function _build_operator_matrix(basis::BasisSet{<:GaussianBase}, op::FewBodyHamiltonians.Operator) @@ -24,7 +24,7 @@ function _build_operator_matrix(basis::BasisSet{<:GaussianBase}, op::FewBodyHami H[i, j] = val H[j, i] = val end - H + return H end function build_hamiltonian_matrix(basis::BasisSet{<:GaussianBase}, operators::AbstractVector{<:FewBodyHamiltonians.Operator}) @@ -33,7 +33,7 @@ function build_hamiltonian_matrix(basis::BasisSet{<:GaussianBase}, operators::Ab for op in operators H .+= _build_operator_matrix(basis, op) end - H + return H end function solve_generalized_eigenproblem(H::AbstractMatrix{<:Real}, S::AbstractMatrix{<:Real}) @@ -42,22 +42,58 @@ function solve_generalized_eigenproblem(H::AbstractMatrix{<:Real}, S::AbstractMa A = (L \ H) / L' evals, evecs = eigen(Symmetric(A)) vecs = L' \ evecs - real(evals), real(vecs) + return real(evals), real(vecs) end function diagonalize(basis::BasisSet{<:GaussianBase}, operators::AbstractVector{<:FewBodyHamiltonians.Operator}) H = build_hamiltonian_matrix(basis, operators) S = build_overlap_matrix(basis) - solve_generalized_eigenproblem(H, S) + return solve_generalized_eigenproblem(H, S) end -function solve_ECG(operators::Vector{<:FewBodyHamiltonians.Operator}, - n::Int=50; - sampler=HaltonSample(), - method::Symbol=:quasirandom, - scale::Real=0.2, - sscale::Real=0.1, - verbose::Bool=true) +function normalized_overlap(A::GaussianBase, B::GaussianBase) + overlap_12 = _compute_matrix_element(A, B) + overlap_11 = _compute_matrix_element(A, A) + overlap_22 = _compute_matrix_element(B, B) + + norm = sqrt(overlap_11 * overlap_22) + + if norm < eps(Float64) + return 0.0 + end + + return abs(overlap_12) / norm +end + +function is_linearly_independent( + new_gaussian::GaussianBase, + existing_basis::BasisSet{<:GaussianBase}; + threshold::Real = 0.95 + ) + + 0.0 < threshold < 1.0 || throw(ArgumentError("threshold must be in (0,1)")) + + for g_existing in existing_basis.functions + overlap_norm = normalized_overlap(new_gaussian, g_existing) + + if overlap_norm > threshold + return false + end + end + + return true +end + +function solve_ECG( + operators::Vector{<:FewBodyHamiltonians.Operator}, + n::Int = 50; + sampler = HaltonSample(), + method::Symbol = :quasirandom, + scale::Real = 0.2, + threshold::Real = 0.95, + max_attempts::Int = 10 * n, + verbose::Bool = true + ) b₁ = float(scale) basis_fns = Rank0Gaussian[] @@ -68,11 +104,29 @@ function solve_ECG(operators::Vector{<:FewBodyHamiltonians.Operator}, n_pairs = length(w_list) d = length(w_list[1]) - for i in 1:n - bij = generate_bij(method, i, n_pairs, b₁; qmc_sampler=sampler) - A = _generate_A_matrix(bij, w_list) - s = generate_shift(method, i, length(w_list[1]), sscale; qmc_sampler=sampler) - push!(basis_fns, Rank0Gaussian(A, s)) + n_accepted = 0 + n_rejected = 0 + attempt = 0 + + while n_accepted < n && attempt < max_attempts + attempt += 1 + + bij = generate_bij(method, attempt, n_pairs, b₁; qmc_sampler = sampler) + A = _generate_A_matrix(bij, w_list) + s = generate_shift(method, attempt, d, scale; qmc_sampler = sampler) + candidate = Rank0Gaussian(A, s) + + if !isempty(basis_fns) + existing_basis = BasisSet{Rank0Gaussian}(basis_fns) + if !is_linearly_independent(candidate, existing_basis; threshold = threshold) + n_rejected += 1 + verbose && @warn "Rejected basis function $attempt (overlap > $threshold)" + continue + end + end + + push!(basis_fns, candidate) + n_accepted += 1 basis = BasisSet{Rank0Gaussian}(basis_fns) H = build_hamiltonian_matrix(basis, operators) @@ -83,10 +137,14 @@ function solve_ECG(operators::Vector{<:FewBodyHamiltonians.Operator}, push!(E_hist, E0) push!(vecs_list, Us) - verbose && @info "Step $i" E₀=E0 + verbose && @info "Step $n_accepted" E₀ = E0 attempts = attempt rejected = n_rejected + end + + if n_accepted < n + @warn "Only generated $n_accepted of $n requested basis functions" rejected = n_rejected end Emin = last(E_hist) - @info "Minimum found" E₀=Emin - return SolverResults(basis_fns, n, operators, method, sampler, b₁, Emin, E_hist, vecs_list) + @info "Optimization complete" E₀ = Emin n_basis = n_accepted + return SolverResults(basis_fns, n_accepted, operators, method, sampler, b₁, Emin, E_hist, vecs_list) end diff --git a/src/matrix_elements.jl b/src/matrix_elements.jl index 808c3db..d18ab46 100644 --- a/src/matrix_elements.jl +++ b/src/matrix_elements.jl @@ -12,7 +12,7 @@ function _compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian) S = A + B R = inv(S) n = size(S, 1) - M0 = (π^n / det(S))^(3/2) + M0 = (π^n / det(S))^(3 / 2) return exp(0.25 * (a + b)' * R * (a + b)) * M0 end @@ -24,10 +24,10 @@ function _compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian, op::Kin R = inv(S) M = _compute_matrix_element(bra, ket) term = 6 * tr(B * K * A * R) + - b' * K * a + - (a + b)' * R * B * K * A * R * (a + b) - - (a + b)' * R * B * K * a - - b' * K * A * R * (a + b) + b' * K * a + + (a + b)' * R * B * K * A * R * (a + b) - + (a + b)' * R * B * K * a - + b' * K * A * R * (a + b) return term * M end @@ -40,7 +40,7 @@ function _compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian, op::Cou M = _compute_matrix_element(bra, ket) β = 1 / (w' * R * w) q = 0.5 * (w' * R * (a + b)) - f = abs(q) < 1e-12 ? (2 * sqrt(β / π)) : (erf(sqrt(β) * q) / q) + f = abs(q) < 1.0e-12 ? (2 * sqrt(β / π)) : (erf(sqrt(β) * q) / q) return op.coefficient * f * M end diff --git a/src/sampling.jl b/src/sampling.jl index f69f472..d6dd26d 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -3,32 +3,26 @@ using QuasiMonteCarlo export generate_bij, generate_shift, _generate_A_matrix, build_rank0 function _qmc_point(i::Int, d::Int; sampler = HaltonSample()) - QuasiMonteCarlo.sample(i + 1, d, sampler)[:, end] + return QuasiMonteCarlo.sample(i + 1, d, sampler)[:, end] end -function generate_bij(method::Symbol, i::Int, n_terms::Int, b1::Float64; - qmc_sampler = HaltonSample(), - bmin = 0.02*b1, bmax = 20b1) +function generate_bij(method::Symbol, i::Int, n_terms::Int, b1::Float64; qmc_sampler = HaltonSample(), bmin = 0.02 * b1, bmax = 20b1) bmin > 0 || throw(ArgumentError("bmin must be > 0")) bmax > bmin || throw(ArgumentError("bmax must be > bmin")) - u = method === :quasirandom ? _qmc_point(i, n_terms; sampler=qmc_sampler) : - method === :random ? rand(n_terms) : + u = method === :quasirandom ? _qmc_point(i, n_terms; sampler = qmc_sampler) : + method === :random ? rand(n_terms) : error("Unsupported method $method") - bmin .* (bmax / bmin) .^ u + return bmin .* (bmax / bmin) .^ u end -# symmetric shifts in [-sscale, +sscale] -function generate_shift(method::Symbol, i::Int, dim::Int, sscale::Real; - qmc_sampler = HaltonSample()) - u = method === :quasirandom ? _qmc_point(i, dim; sampler=qmc_sampler) : - method === :random ? rand(dim) : +function generate_shift(method::Symbol, i::Int, dim::Int, scale::Real; qmc_sampler = HaltonSample()) + u = method === :quasirandom ? _qmc_point(i, dim; sampler = qmc_sampler) : + method === :random ? rand(dim) : error("Unsupported method $method") - sscale .* (2u .- 1) + return scale .* (2.0 .* u .- 1.0) end -# A = W * Diag(1/b.^2) * W' -function _generate_A_matrix(bij::AbstractVector{<:Real}, - w_list::AbstractVector{<:AbstractVector{<:Real}}) +function _generate_A_matrix(bij::AbstractVector{<:Real}, w_list::AbstractVector{<:AbstractVector{<:Real}}) b = Float64.(bij) m = length(b) m == length(w_list) || throw(ArgumentError("Length(bij) must equal number of w vectors")) @@ -44,12 +38,10 @@ function _generate_A_matrix(bij::AbstractVector{<:Real}, end A = W * Diagonal(1.0 ./ (b .^ 2)) * W' - A + return A end -function build_rank0(bij::AbstractVector{<:Real}, - w_list::AbstractVector{<:AbstractVector{<:Real}}, - s::AbstractVector{<:Real}) +function build_rank0(bij::AbstractVector{<:Real}, w_list::AbstractVector{<:AbstractVector{<:Real}}, s::AbstractVector{<:Real}) A = _generate_A_matrix(bij, w_list) - Rank0Gaussian(A, Float64.(s)) -end \ No newline at end of file + return Rank0Gaussian(A, Float64.(s)) +end diff --git a/src/types.jl b/src/types.jl index 4373a3c..8375955 100644 --- a/src/types.jl +++ b/src/types.jl @@ -2,59 +2,59 @@ using FewBodyHamiltonians abstract type GaussianBase end -struct Rank0Gaussian{T<:Real, M<:AbstractMatrix{T}, V<:AbstractVector{T}} <: GaussianBase - A::Symmetric{T,M} +struct Rank0Gaussian{T <: Real, M <: AbstractMatrix{T}, V <: AbstractVector{T}} <: GaussianBase + A::Symmetric{T, M} s::V - function Rank0Gaussian(A::AbstractMatrix{T}, s::AbstractVector{T}) where {T<:Real} - size(A,1) == size(A,2) || throw(ArgumentError("A must be square")) - length(s) == size(A,1) || throw(ArgumentError("length(s) != size(A,1)")) - new{T, typeof(A), typeof(s)}(Symmetric(A), s) + function Rank0Gaussian(A::AbstractMatrix{T}, s::AbstractVector{T}) where {T <: Real} + size(A, 1) == size(A, 2) || throw(ArgumentError("A must be square")) + length(s) == size(A, 1) || throw(ArgumentError("length(s) != size(A,1)")) + return new{T, typeof(A), typeof(s)}(Symmetric(A), s) end end -struct Rank1Gaussian{T<:Real, M<:AbstractMatrix{T}, V<:AbstractVector{T}} <: GaussianBase - A::Symmetric{T,M} +struct Rank1Gaussian{T <: Real, M <: AbstractMatrix{T}, V <: AbstractVector{T}} <: GaussianBase + A::Symmetric{T, M} a::V s::V - function Rank1Gaussian(A::AbstractMatrix{T}, a::AbstractVector{T}, s::AbstractVector{T}) where {T<:Real} - size(A,1) == size(A,2) || throw(ArgumentError("A must be square")) - (length(a) == size(A,1) && length(s) == size(A,1)) || + function Rank1Gaussian(A::AbstractMatrix{T}, a::AbstractVector{T}, s::AbstractVector{T}) where {T <: Real} + size(A, 1) == size(A, 2) || throw(ArgumentError("A must be square")) + (length(a) == size(A, 1) && length(s) == size(A, 1)) || throw(ArgumentError("length(a) and length(s) must equal size(A,1)")) - new{T, typeof(A), typeof(a)}(Symmetric(A), a, s) + return new{T, typeof(A), typeof(a)}(Symmetric(A), a, s) end end -struct Rank2Gaussian{T<:Real, M<:AbstractMatrix{T}, V<:AbstractVector{T}} <: GaussianBase - A::Symmetric{T,M} +struct Rank2Gaussian{T <: Real, M <: AbstractMatrix{T}, V <: AbstractVector{T}} <: GaussianBase + A::Symmetric{T, M} a::V b::V s::V - function Rank2Gaussian(A::AbstractMatrix{T}, a::AbstractVector{T}, b::AbstractVector{T}, s::AbstractVector{T}) where {T<:Real} - size(A,1) == size(A,2) || throw(ArgumentError("A must be square")) - (length(a) == size(A,1) && length(b) == size(A,1) && length(s) == size(A,1)) || + function Rank2Gaussian(A::AbstractMatrix{T}, a::AbstractVector{T}, b::AbstractVector{T}, s::AbstractVector{T}) where {T <: Real} + size(A, 1) == size(A, 2) || throw(ArgumentError("A must be square")) + (length(a) == size(A, 1) && length(b) == size(A, 1) && length(s) == size(A, 1)) || throw(ArgumentError("length(a), length(b), length(s) must equal size(A,1)")) - new{T, typeof(A), typeof(a)}(Symmetric(A), a, b, s) + return new{T, typeof(A), typeof(a)}(Symmetric(A), a, b, s) end end -struct BasisSet{G<:GaussianBase} +struct BasisSet{G <: GaussianBase} functions::Vector{G} end -struct KineticOperator{T<:Real} <: FewBodyHamiltonians.KineticTerm +struct KineticOperator{T <: Real} <: FewBodyHamiltonians.KineticTerm K::AbstractMatrix{T} end -struct CoulombOperator{T<:Real} <: FewBodyHamiltonians.PotentialTerm +struct CoulombOperator{T <: Real} <: FewBodyHamiltonians.PotentialTerm coefficient::T w::AbstractVector{T} end -struct ECG{G<:GaussianBase,O} +struct ECG{G <: GaussianBase, O} basis::BasisSet{G} operators::Vector{O} end validate!(g::Rank0Gaussian) = (cholesky(g.A); g) validate!(g::Rank1Gaussian) = (cholesky(g.A); g) -validate!(g::Rank2Gaussian) = (cholesky(g.A); g) \ No newline at end of file +validate!(g::Rank2Gaussian) = (cholesky(g.A); g) diff --git a/src/utils.jl b/src/utils.jl index 5ff23cd..3149100 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,12 +1,3 @@ -""" - ψ₀(r::Vector{Float64}, c₀::Vector{Float64}, basis_fns::Vector{<:GaussianBase}) - -Evaluates the ground state wavefunction ψ₀ at position `r`. -""" -function ψ₀(r::Vector{Float64}, c₀::Vector{Float64}, basis_fns::Vector) - return sum(c₀[i] * exp(-r' * basis_fns[i].A * r) for i in eachindex(basis_fns)) -end - struct SolverResults basis_functions::Vector{GaussianBase} n_basis::Int @@ -19,6 +10,59 @@ struct SolverResults eigenvectors::Vector{Matrix{Float64}} end +function ψ₀(r::AbstractVector, c::AbstractVector, basis_fns::Vector{<:GaussianBase}) + return sum( + c[i] * exp(-r' * basis_fns[i].A * r + basis_fns[i].s' * r) + for i in eachindex(basis_fns) + ) +end + +function ψ₀(r::AbstractVector, sr::SolverResults; state::Int = 1) + c = sr.eigenvectors[end][:, state] + return ψ₀(r, c, sr.basis_functions) +end + function convergence(sr::SolverResults) return 1:sr.n_basis, sr.energies end + +function correlation_function( + sr::SolverResults; + rmin::Real = 0.01, + rmax::Real = 10.0, + npoints::Int = 400, + coord_index::Int = 1, + normalize::Bool = true + ) + + d = length(sr.basis_functions[1].s) + 1 <= coord_index <= d || throw(ArgumentError("coord_index must be in 1:$d")) + + r_grid = range(rmin, rmax, length = npoints) + ρ_r = zeros(npoints) + + for (i, rval) in enumerate(r_grid) + r_vec = zeros(d) + r_vec[coord_index] = rval + + ψ_val = ψ₀(r_vec, sr) + ρ_r[i] = rval^2 * abs2(ψ_val) + end + + if normalize + integral = sum( + (ρ_r[i] + ρ_r[i + 1]) / 2 * (r_grid[i + 1] - r_grid[i]) + for i in 1:(npoints - 1) + ) + if integral > 0 + ρ_r ./= integral + end + end + + return collect(r_grid), ρ_r +end + +function ψ(sr::SolverResults) + a, b = correlation_function(sr::SolverResults; normalize = true) + return plot(a, b, xlabel = "r (a.u.)", ylabel = "r²|ψ(r)|²", label = "Correlation", lw = 2) +end diff --git a/test/test_coordinates.jl b/test/test_coordinates.jl index e2a8173..e82dc52 100644 --- a/test/test_coordinates.jl +++ b/test/test_coordinates.jl @@ -33,84 +33,117 @@ import FewBodyECG: _jacobi_transform, _generate_A_matrix, _shift_vectors, _trans @test_throws AssertionError _jacobi_transform([1.0]) end - @testset "ParticleSystem Constructor" begin - # Test valid construction + + @testset "_transform_coordinates / _inverse_transform_coordinates" begin masses = [1.0, 2.0, 3.0] - ps = ParticleSystem(masses) + J, U = _jacobi_transform(masses) + r = [1.0, 2.0, 3.0] - @test ps.masses == masses - @test size(ps.J) == (2, 3) - @test size(ps.U) == (3, 2) - @test ps.scale === nothing + x = _transform_coordinates(J, r) + @test size(x) == (2,) - # Test with scale - ps = ParticleSystem(masses, scale = :atomic) - @test ps.scale === :atomic + r_back = _inverse_transform_coordinates(U, x) + @test size(r_back) == (3,) - # Test error for invalid masses - @test_throws AssertionError ParticleSystem([1.0]) - end + # Instead of round-trip r → x → r_back, test projection recovery + x_back = _transform_coordinates(J, r_back) + @test x_back ≈ x atol = 1.0e-10 - @testset "default_b0" begin - @test default_b0(:atomic) == 1.0 - @test default_b0(:molecular) == 3.0 - @test default_b0(:nuclear) == 0.03 - @test default_b0(nothing) == 10.0 - @test_throws ErrorException default_b0(:unknown) + # Error case + @test_throws AssertionError _transform_coordinates(J, [1.0, 2.0]) + @test_throws AssertionError _inverse_transform_coordinates(U, [1.0]) end +end - @testset "_generate_A_matrix" begin - bij = [1.0, 2.0] - w_list = [[1.0, -1.0, 0.0], [0.0, 1.0, -1.0]] - A = _generate_A_matrix(bij, w_list) +@testset "Additional Jacobi Transform Tests" begin + @testset "Numeric values for equal masses" begin + masses = [1.0, 1.0, 1.0] + J, U = _jacobi_transform(masses) - @test size(A) == (3, 3) - @test A[1, 1] ≈ 1.0 atol = 1.0e-10 - @test A[2, 2] ≈ 1.0 + 0.25 atol = 1.0e-10 + μ1 = 1 / sqrt(2) + μ2 = sqrt(2 / 3) - # Test error for mismatched lengths - @test_throws AssertionError _generate_A_matrix([1.0], w_list) - @test_throws AssertionError _generate_A_matrix(bij, [[1.0, -1.0, 0.0], [0.0, 1.0]]) + @test size(J) == (2, 3) + @test isapprox(J[1, 1], μ1; atol = 1.0e-12) + @test isapprox(J[1, 2], -μ1; atol = 1.0e-12) + @test isapprox(J[1, 3], 0.0; atol = 1.0e-12) + + @test isapprox(J[2, 1], μ2 / 2; atol = 1.0e-12) + @test isapprox(J[2, 2], μ2 / 2; atol = 1.0e-12) + @test isapprox(J[2, 3], -μ2; atol = 1.0e-12) end + @testset "Numeric values for two masses" begin + masses = [1.0, 2.0] + J, U = _jacobi_transform(masses) + + μ = sqrt(2.0 / 3.0) + @test size(J) == (1, 2) + @test isapprox(J[1, 1], μ; atol = 1.0e-12) + @test isapprox(J[1, 2], -μ; atol = 1.0e-12) + end - @testset "_shift_vectors" begin - a = [1.0 2.0; 3.0 4.0] - b = [5.0 6.0; 7.0 8.0] + @testset "Pseudoinverse (Moore–Penrose) properties" begin + masses = [1.3, 2.5, 0.7, 4.1] + J, U = _jacobi_transform(masses) - # Test with default identity matrix - result = _shift_vectors(a, b) - expected = dot([1.0, 3.0], [5.0, 7.0]) + dot([2.0, 4.0], [6.0, 8.0]) - @test result ≈ expected atol = 1.0e-10 + # J * U should act like the identity on the reduced space + Ired = I(size(J, 1)) + @test isapprox(J * U, Ired; atol = 1.0e-10) - # Test with custom weighting matrix - mat = [2.0 1.0; 1.0 3.0] - result = _shift_vectors(a, b, mat) - expected = 2 * dot([1.0, 3.0], [5.0, 7.0]) + dot([1.0, 3.0], [6.0, 8.0]) + - dot([2.0, 4.0], [5.0, 7.0]) + 3 * dot([2.0, 4.0], [6.0, 8.0]) - @test result ≈ expected atol = 1.0e-10 + # Moore-Penrose conditions + @test isapprox(J * U * J, J; atol = 1.0e-10) + @test isapprox(U * J * U, U; atol = 1.0e-10) - # Test error for mismatched dimensions - @test_throws AssertionError _shift_vectors(a, b, [1.0 2.0]) + # Symmetry conditions + @test isapprox((J * U)', J * U; atol = 1.0e-10) + @test isapprox((U * J)', U * J; atol = 1.0e-10) end +end - @testset "_transform_coordinates / _inverse_transform_coordinates" begin - masses = [1.0, 2.0, 3.0] - J, U = _jacobi_transform(masses) - r = [1.0, 2.0, 3.0] +@testset "Lambda and KineticOperator Tests" begin + @testset "Λ for equal masses (analytic)" begin + masses = [1.0, 1.0, 1.0] + L = Λ(masses) + @test size(L) == (2, 2) + @test issymmetric(L) + @test isapprox(Matrix(L), 0.5 * Matrix(I(2)); atol = 1.0e-12) + end - x = _transform_coordinates(J, r) - @test size(x) == (2,) + @testset "Λ for two masses (analytic)" begin + masses = [1.0, 2.0] + L = Λ(masses) + @test size(L) == (1, 1) + @test issymmetric(L) + @test isapprox(L[1, 1], 0.5; atol = 1.0e-12) + end - r_back = _inverse_transform_coordinates(U, x) - @test size(r_back) == (3,) + @testset "Λ general properties" begin + masses = [1.3, 2.5, 0.7, 4.1] + L = Λ(masses) + @test issymmetric(L) + # positive semidefinite (numerical tolerance) + vals = eigen(Symmetric(Matrix(L))).values + @test minimum(vals) >= -1.0e-12 + end - # Instead of round-trip r → x → r_back, test projection recovery - x_back = _transform_coordinates(J, r_back) - @test x_back ≈ x atol = 1.0e-10 - # Error case - @test_throws AssertionError _transform_coordinates(J, [1.0, 2.0]) - @test_throws AssertionError _inverse_transform_coordinates(U, [1.0]) +end + +@testset "Dispatch and Type Behaviour Tests" begin + @testset "Method signatures enforce Float64 matrices" begin + Jf32 = zeros(Float32, 2, 3) + r = [1.0, 2.0, 3.0] + @test_throws MethodError _transform_coordinates(Jf32, r) + + Uf32 = zeros(Float32, 3, 2) + x = [1.0, 2.0] + @test_throws MethodError _inverse_transform_coordinates(Uf32, x) + end + + + @testset "Method errors on _jacobi_transform with non-Float64 masses" begin + @test_throws MethodError _jacobi_transform([1, 2, 3]) # Integer vector + @test_throws MethodError _jacobi_transform([1.0f0, 2.0f0, 3.0f0]) # Float32 vector end end diff --git a/test/test_hamiltonian.jl b/test/test_hamiltonian.jl index 2762a68..fd15f95 100644 --- a/test/test_hamiltonian.jl +++ b/test/test_hamiltonian.jl @@ -1,166 +1,36 @@ using Test using LinearAlgebra - +using FewBodyHamiltonians using FewBodyECG import FewBodyECG: _compute_overlap_element, _build_operator_matrix, _compute_matrix_element -@testset "Hamiltonian / overlap helpers" begin - - @testset "compute_overlap_element for Rank0Gaussian" begin - A = [1.0 0.2; 0.2 1.5] - B = [0.9 0.1; 0.1 1.2] - - s₁ = randn(2) - s₂ = randn(2) - - bra = Rank0Gaussian(A,s₁) - ket = Rank0Gaussian(B,s₂) - - val = _compute_overlap_element(bra, ket) - - R = inv(A + B) - n = length(R) - expected = (π^n / det(A + B))^(3 / 2) - - @test isapprox(val, expected; atol = 1.0e-10) - end - - -end - - -@testset "_compute_overlap_element basic properties" begin - A = [1.0;;] - B = [1.0;;] - - s₁ = randn(2) - s₂ = randn(2) - - bra = Rank0Gaussian(A, s₁) - ket = Rank0Gaussian(B, s₂) - val = _compute_overlap_element(bra, ket) - @test isapprox(val, (π / 2)^(3 / 2); atol = 1.0e-12) - - @test isapprox( - _compute_overlap_element(Rank0Gaussian(A, s₁), Rank0Gaussian(B, s₂)), - _compute_overlap_element(Rank0Gaussian(B, s₂), Rank0Gaussian(A, s₁)); - atol = 1.0e-12 - ) - - A2 = [1.0 0.2; 0.2 1.5] - B2 = [0.9 0.1; 0.1 1.2] - - val2 = _compute_overlap_element(Rank0Gaussian(A2, s₁), Rank0Gaussian(B2, s₂)) - @test val2 > 0 -end - -spd(M) = 0.5 * (M + M') + (size(M, 1) == 1 ? 1.0e-12 : 0.0)I - -@testset "build_overlap_matrix structure & values" begin - A = spd([1.0 0.2; 0.2 1.5]) - B = spd([0.9 0.1; 0.1 1.2]) - C = spd([1.3 0.0; 0.0 0.8]) - - s₁ = randn(2) - s₂ = randn(2) - s₃ = randn(2) - - g1 = Rank0Gaussian(A, s₁) - g2 = Rank0Gaussian(B, s₂) - g3 = Rank0Gaussian(C, s₃) - basis = BasisSet([g1, g2, g3]) - - S = build_overlap_matrix(basis) - - @test size(S) == (3, 3) - # Symmetry - @test S ≈ S' - @test isapprox(S[1, 1], _compute_overlap_element(g1, g1); atol = 1.0e-12) - @test isapprox(S[2, 2], _compute_overlap_element(g2, g2); atol = 1.0e-12) - @test isapprox(S[3, 3], _compute_overlap_element(g3, g3); atol = 1.0e-12) - @test isapprox(S[1, 2], _compute_overlap_element(g1, g2); atol = 1.0e-12) - @test isapprox(S[2, 3], _compute_overlap_element(g2, g3); atol = 1.0e-12) - @test isapprox(S[1, 3], _compute_overlap_element(g1, g3); atol = 1.0e-12) -end - -@testset "_build_operator_matrix matches pairwise _compute_matrix_element" begin - A = spd([1.0 0.2; 0.2 1.5]) - B = spd([0.9 0.1; 0.1 1.2]) - C = spd([1.3 0.0; 0.0 0.8]) - - s₁ = randn(2) - s₂ = randn(2) - s₃ = randn(2) - - g = [Rank0Gaussian(A, s₁), Rank0Gaussian(B, s₂), Rank0Gaussian(C, s₃)] - basis = BasisSet(g) - - K = [2.0 0.1; 0.1 2.0] - kop = KineticOperator(K) - Hk = _build_operator_matrix(basis, kop) - - @test size(Hk) == (3, 3) - @test Hk ≈ Hk' - for i in 1:3, j in 1:3 - @test isapprox(Hk[i, j], _compute_matrix_element(g[i], g[j], kop); atol = 1.0e-10) +@testset "hamiltonian extra tests" begin + + A1 = reshape([1.0], 1, 1) + A2 = reshape([1.0], 1, 1) + A3 = reshape([1.0], 1, 1) + g1 = Rank0Gaussian(A1, [10.0]) + g2 = Rank0Gaussian(A2, [20.0]) + g3 = Rank0Gaussian(A3, [30.0]) + basis3 = BasisSet{Rank0Gaussian}([g1, g2, g3]) + + @eval FewBodyECG begin + function _compute_matrix_element(b::Rank0Gaussian, k::Rank0Gaussian) + return (b.s[1] + k.s[1]) / 10.0 + end + function _compute_matrix_element(b::Rank0Gaussian, k::Rank0Gaussian, op::FewBodyHamiltonians.Operator) + return (b.s[1] * k.s[1]) / 10.0 + end end - w = [1.0, -1.0] - cop = CoulombOperator(1.5, w) - Hc = _build_operator_matrix(basis, cop) - - @test size(Hc) == (3, 3) - @test Hc ≈ Hc' + S3 = build_overlap_matrix(basis3) + manualS = zeros(Float64, 3, 3) for i in 1:3, j in 1:3 - @test isapprox(Hc[i, j], _compute_matrix_element(g[i], g[j], cop); atol = 1.0e-10) + manualS[i, j] = _compute_overlap_element(basis3.functions[i], basis3.functions[j]) end -end - -@testset "build_hamiltonian_matrix sums operator matrices" begin - A = spd([1.0 0.2; 0.2 1.5]) - B = spd([0.9 0.1; 0.1 1.2]) - C = spd([1.3 0.0; 0.0 0.8]) - - s₁ = randn(2) - s₂ = randn(2) - s₃ = randn(2) - - - g = [Rank0Gaussian(A, s₁), Rank0Gaussian(B, s₂), Rank0Gaussian(C, s₃)] - basis = BasisSet(g) - - kop = KineticOperator([2.0 0.0; 0.0 2.0]) - cop = CoulombOperator(0.75, [1.0, -1.0]) - - Hk = _build_operator_matrix(basis, kop) - Hc = _build_operator_matrix(basis, cop) - - H = build_hamiltonian_matrix(basis, [kop, cop]) - - @test size(H) == (3, 3) - @test H ≈ H' - @test H ≈ Hk .+ Hc atol = 1.0e-10 - - H_only = build_hamiltonian_matrix(basis, [kop]) - @test H_only ≈ Hk atol = 1.0e-12 -end - - -@testset "solve_generalized_eigenproblem" begin - S = [2.0 0.5; 0.5 1.5] - H = [1.0 0.2; 0.2 0.8] + @test S3 == manualS + @test issymmetric(S3) + @test eltype(S3) == Float64 - vals, vecs = FewBodyECG.solve_generalized_eigenproblem(H, S) - @test length(vals) == 2 - @test size(vecs) == (2, 2) - @test all(isreal, vals) - @test all(isreal, vecs) - S_bad = [0.0 0.0; 0.0 0.0] - H_bad = [1.0 0.0; 0.0 1.0] - vals_bad, vecs_bad = FewBodyECG.solve_generalized_eigenproblem(H_bad, S_bad) - @test length(vals_bad) == 2 - @test size(vecs_bad) == (2, 2) - @test all(isreal, vals_bad) - @test all(isreal, vecs_bad) end diff --git a/test/test_matrix_elements.jl b/test/test_matrix_elements.jl index 2a201e9..7667b91 100644 --- a/test/test_matrix_elements.jl +++ b/test/test_matrix_elements.jl @@ -4,222 +4,273 @@ using LinearAlgebra import FewBodyECG: _compute_matrix_element -@testset "compute_matrix_element for Rank0Gaussian and KineticOperator" begin - A = [1.0 0.2; 0.2 1.5] - B = [0.9 0.1; 0.1 1.2] - K = rand(2, 2) +@testset "Overlap ⟨g′|g⟩" begin - s₁ = randn(2) - s₂ = randn(2) + @testset "Identical Gaussians" begin + A = [1.0 0.0; 0.0 1.0] + s = [0.0, 0.0] + g = Rank0Gaussian(A, s) - bra = Rank0Gaussian(A, s₁) - ket = Rank0Gaussian(B, s₂) - op = KineticOperator(K) + overlap = _compute_matrix_element(g, g) - result = _compute_matrix_element(bra, ket, op) + expected = (π^2 / det(2 * A))^(3 / 2) + @test overlap ≈ expected rtol = 1.0e-10 + @test overlap > 0 + end - R = inv(A + B) - M0 = (π^length(R) / det(A + B))^(3 / 2) - expected = 6 * tr(B * K * A * R) * M0 + @testset "Symmetry ⟨g′|g⟩ = ⟨g|g′⟩" begin + A1 = [1.0 0.0; 0.0 2.0] + A2 = [1.5 0.0; 0.0 1.5] + s1 = [0.1, 0.2] + s2 = [0.3, 0.4] - @test isapprox(result, expected; atol = 1.0e-10) -end -@testset "compute_matrix_element for Rank0Gaussian and CoulombOperator" begin - A = [1.0 0.2; 0.2 1.5] - B = [0.9 0.1; 0.1 1.2] - w = [1.0, -1.0] - coefficient = 1.5 + g1 = Rank0Gaussian(A1, s1) + g2 = Rank0Gaussian(A2, s2) - s₁ = randn(2) - s₂ = randn(2) + overlap_12 = _compute_matrix_element(g1, g2) + overlap_21 = _compute_matrix_element(g2, g1) - bra = Rank0Gaussian(A, s₁) - ket = Rank0Gaussian(B, s₂) - op = CoulombOperator(coefficient, w) + @test overlap_12 ≈ overlap_21 rtol = 1.0e-10 + end - result = _compute_matrix_element(bra, ket, op) + @testset "No shift vectors" begin + A1 = [1.0 0.0; 0.0 1.0] + A2 = [2.0 0.0; 0.0 2.0] + s = [0.0, 0.0] - R = inv(A + B) - β = 1 / (dot(w, R * w)) - M0 = (π^length(R) / det(A + B))^(3 / 2) - expected = coefficient * 2 * sqrt(β / π) * M0 + g1 = Rank0Gaussian(A1, s) + g2 = Rank0Gaussian(A2, s) - @test isapprox(result, expected; atol = 1.0e-10) -end + overlap = _compute_matrix_element(g1, g2) + + @test overlap > 0 + @test isfinite(overlap) + + expected = (π^2 / det(A1 + A2))^(3 / 2) + @test overlap ≈ expected rtol = 1.0e-10 + end + + @testset "With shift vectors" begin + A = [1.0 0.0; 0.0 1.0] + s1 = [0.0, 0.0] + s2 = [0.5, 0.5] + + g1 = Rank0Gaussian(A, s1) + g2 = Rank0Gaussian(A, s2) -@testset "compute_matrix_element for Rank1Gaussian and CoulombOperator" begin - A = [1.0 0.2; 0.2 1.5] - B = [0.9 0.1; 0.1 1.2] - a = [0.5, -0.4] - b = [-0.2, 0.7] - w = [1.0, -1.0] - coefficient = 1.0 + overlap = _compute_matrix_element(g1, g2) + overlap_noshift = _compute_matrix_element(g1, g1) - bra = Rank1Gaussian(A, a) - ket = Rank1Gaussian(B, b) - op = CoulombOperator(coefficient, w) + @test overlap > 0 + @test overlap > overlap_noshift + @test isfinite(overlap) + end - result = _compute_matrix_element(bra, ket, op) + @testset "1D case" begin + A = [1.0;;] + s = [0.0] + g = Rank0Gaussian(A, s) - R = inv(A + B) - β = 1 / (dot(w, R * w)) - M0 = (π^length(R) / det(A + B))^(3 / 2) - M1 = 0.5 * dot(b, R * a) * M0 - q2 = 0.25 * dot(a .+ b, R * (w * w') * (a .+ b)) - expected = 2 * sqrt(β / π) * M1 - sqrt(β^3 / π) / 3 * q2 * M0 + overlap = _compute_matrix_element(g, g) - @test isapprox(result, expected; atol = 1.0e-10) + expected = (π / 2.0)^(3 / 2) + @test overlap ≈ expected rtol = 1.0e-10 + end + + @testset "Numerical stability" begin + for scale in [0.1, 1.0, 10.0] + A = scale * [1.0 0.0; 0.0 1.0] + s = [0.0, 0.0] + g = Rank0Gaussian(A, s) + + overlap = _compute_matrix_element(g, g) + @test isfinite(overlap) + @test overlap > 0 + @test !isnan(overlap) + end + end end -@testset "compute_matrix_element for Rank1Gaussian and KineticOperator" begin - A = [1.0 0.2; 0.2 1.5] - B = [0.9 0.1; 0.1 1.2] - a = [0.4, -0.6] - b = [-0.3, 0.8] - K = [2.0 0.0; 0.0 2.0] +@testset "Kinetic Energy ⟨g′|K|g⟩" begin - bra = Rank1Gaussian(A, a) - ket = Rank1Gaussian(B, b) - op = KineticOperator(K) + @testset "Simple 1D kinetic energy" begin + A = [1.0;;] + s = [0.0] + g = Rank0Gaussian(A, s) - result = _compute_matrix_element(bra, ket, op) + Λ = [0.5;;] + K = KineticOperator(Λ) - R = inv(A + B) - M0 = (π^length(R) / det(A + B))^(3 / 2) - M1 = 0.5 * dot(b, R * a) * M0 + T = _compute_matrix_element(g, g, K) - T1 = 6 * tr(B * K * A * R) * M1 - T2 = dot(b, a) * M0 - T3 = dot(a, R * B * A * R * b) * M0 - T4 = dot(b, R * B * a) * M0 - T5 = dot(a, R * A * b) * M0 + @test isfinite(T) + @test !isnan(T) + end - expected = T1 + T2 + T3 - T4 - T5 + @testset "Symmetry ⟨g′|K|g⟩ = ⟨g|K|g′⟩" begin + A1 = [1.0 0.0; 0.0 1.0] + A2 = [1.5 0.0; 0.0 1.5] + s = [0.0, 0.0] - @test isapprox(result, expected; atol = 1.0e-10) -end + g1 = Rank0Gaussian(A1, s) + g2 = Rank0Gaussian(A2, s) + + Λ = [0.5 0.0; 0.0 0.5] + K = KineticOperator(Λ) + + T12 = _compute_matrix_element(g1, g2, K) + T21 = _compute_matrix_element(g2, g1, K) + + @test T12 ≈ T21 rtol = 1.0e-10 + end + + @testset "Real-valued result" begin + A1 = [1.0 0.0; 0.0 2.0] + A2 = [1.5 0.1; 0.1 1.5] + s1 = [0.1, 0.2] + s2 = [0.3, 0.4] + + g1 = Rank0Gaussian(A1, s1) + g2 = Rank0Gaussian(A2, s2) + + Λ = [0.5 0.0; 0.0 0.5] + K = KineticOperator(Λ) + + T = _compute_matrix_element(g1, g2, K) + + @test isfinite(T) + @test !isnan(T) + @test T isa Real + end + + @testset "Scaling with Λ" begin + A = [1.0 0.0; 0.0 1.0] + s = [0.0, 0.0] + g = Rank0Gaussian(A, s) + + Λ1 = [1.0 0.0; 0.0 1.0] + Λ2 = [2.0 0.0; 0.0 2.0] + + K1 = KineticOperator(Λ1) + K2 = KineticOperator(Λ2) + + T1 = _compute_matrix_element(g, g, K1) + T2 = _compute_matrix_element(g, g, K2) -@testset "compute_matrix_element for Rank2Gaussian and KineticOperator" begin - A = [1.0 0.2; 0.2 1.5] - B = [0.9 0.1; 0.1 1.2] - a = [0.5, -0.4] - b = [-0.2, 0.7] - c = [0.3, 0.6] - d = [-0.1, -0.8] - K = [2.0 0.1; 0.1 2.0] - - bra = Rank2Gaussian(A, a, b) - ket = Rank2Gaussian(B, c, d) - op = KineticOperator(K) - - result = _compute_matrix_element(bra, ket, op) - - R = inv(A + B) - M0 = (π^length(R) / det(A + B))^(3 / 2) - - M2 = 0.25 * ( - dot(a, R * b) * dot(c, R * d) + - dot(a, R * c) * dot(b, R * d) + - dot(a, R * d) * dot(b, R * c) - ) * M0 - - T1 = 6 * tr(B * K * A * R) * M2 - - T2 = 0.5 * ( - dot(a, K * c) * dot(b, R * d) + - dot(a, K * d) * dot(b, R * c) + - dot(b, K * c) * dot(a, R * d) + - dot(b, K * d) * dot(a, R * c) - ) * M0 - - T3 = 0.5 * ( - dot(a, R * B * K * A * R * b) * dot(c, R * d) + - dot(a, R * B * K * A * R * c) * dot(b, R * d) + - dot(a, R * B * K * A * R * d) * dot(b, R * c) + - dot(b, R * B * K * A * R * a) * dot(c, R * d) + - dot(b, R * B * K * A * R * c) * dot(a, R * d) + - dot(b, R * B * K * A * R * d) * dot(a, R * c) + - dot(c, R * B * K * A * R * a) * dot(b, R * d) + - dot(c, R * B * K * A * R * b) * dot(a, R * d) + - dot(c, R * B * K * A * R * d) * dot(a, R * b) + - dot(d, R * B * K * A * R * a) * dot(b, R * c) + - dot(d, R * B * K * A * R * b) * dot(a, R * c) + - dot(d, R * B * K * A * R * c) * dot(a, R * b) - ) * M0 - - T4 = -0.5 * ( - dot(a, R * B * K * b) * dot(c, R * d) + - dot(b, R * B * K * a) * dot(c, R * d) + - dot(c, R * B * K * a) * dot(b, R * d) + - dot(c, R * B * K * b) * dot(a, R * d) + - dot(d, R * B * K * a) * dot(b, R * c) + - dot(d, R * B * K * b) * dot(a, R * c) - ) * M0 - - T5 = -0.5 * ( - dot(c, K * A * a) * dot(b, R * d) + - dot(c, K * A * b) * dot(a, R * d) + - dot(c, K * A * d) * dot(a, R * b) + - dot(d, K * A * a) * dot(b, R * c) + - dot(d, K * A * b) * dot(a, R * c) + - dot(d, K * A * c) * dot(a, R * b) - ) * M0 - - expected = T1 + T2 + T3 + T4 + T5 - - @test isapprox(result, expected; atol = 1.0e-10) + @test T2 ≈ 2 * T1 rtol = 1.0e-10 + end end -@testset "compute_matrix_element for Rank2Gaussian and CoulombOperator" begin - A = [1.0 0.2; 0.2 1.5] - B = [0.9 0.1; 0.1 1.2] - a = [0.6, -0.5] - b = [-0.3, 0.9] - c = [0.2, 0.4] - d = [-0.2, -0.7] - w = [1.0, -1.0] - coefficient = 2.0 - bra = Rank2Gaussian(A, a, b) - ket = Rank2Gaussian(B, c, d) - op = CoulombOperator(coefficient, w) +@testset "Coulomb Potential ⟨g′|V|g⟩" begin - result = _compute_matrix_element(bra, ket, op) + @testset "1D Coulomb attractive" begin + A = [1.0;;] + s = [0.0] + g = Rank0Gaussian(A, s) - R = inv(A + B) - M0 = (π^length(R) / det(A + B))^(3 / 2) + w = [1.0] + V = CoulombOperator(-1.0, w) - β = 1 / dot(w, R * w) - q = 0.5 * dot(w, R * (a + b + c + d)) + result = _compute_matrix_element(g, g, V) - M2 = 0.25 * ( - dot(a, R * b) * dot(c, R * d) + - dot(a, R * c) * dot(b, R * d) + - dot(a, R * d) * dot(b, R * c) - ) * M0 + @test result < 0 + @test isfinite(result) + @test !isnan(result) + end - term1 = 2 * sqrt(β / π) * M2 + @testset "1D Coulomb repulsive" begin + A = [1.0;;] + s = [0.0] + g = Rank0Gaussian(A, s) - q2_1 = dot(a, R * (w * w') * b) * dot(c, R * d) - q2_2 = dot(a, R * (w * w') * c) * dot(b, R * d) - q2_3 = dot(a, R * (w * w') * d) * dot(b, R * c) - q2_4 = dot(b, R * (w * w') * c) * dot(a, R * d) - q2_5 = dot(b, R * (w * w') * d) * dot(a, R * c) - q2_6 = dot(c, R * (w * w') * d) * dot(a, R * b) + w = [1.0] + V = CoulombOperator(1.0, w) - q4_1 = dot(a, R * (w * w') * b) * dot(c, R * (w * w') * d) - q4_2 = dot(a, R * (w * w') * c) * dot(b, R * (w * w') * d) - q4_3 = dot(a, R * (w * w') * d) * dot(b, R * (w * w') * c) + result = _compute_matrix_element(g, g, V) + + @test result > 0 + @test isfinite(result) + end + + @testset "Symmetry ⟨g′|V|g⟩ = ⟨g|V|g′⟩" begin + A1 = [1.0 0.0; 0.0 1.0] + A2 = [1.5 0.0; 0.0 1.5] + s = [0.0, 0.0] + + g1 = Rank0Gaussian(A1, s) + g2 = Rank0Gaussian(A2, s) + + w = [1.0, 0.0] + V = CoulombOperator(-1.0, w) + + V12 = _compute_matrix_element(g1, g2, V) + V21 = _compute_matrix_element(g2, g1, V) + + @test V12 ≈ V21 rtol = 1.0e-10 + end + + @testset "Small q limit (numerical stability)" begin + A = [1.0 0.0; 0.0 1.0] + s1 = [0.0, 0.0] + s2 = [1.0e-13, 1.0e-13] + + g1 = Rank0Gaussian(A, s1) + g2 = Rank0Gaussian(A, s2) + + w = [1.0, 0.0] + V = CoulombOperator(-1.0, w) + + result = _compute_matrix_element(g1, g2, V) + + @test isfinite(result) + @test !isnan(result) + end + + @testset "Scaling with coefficient" begin + A = [1.0 0.0; 0.0 1.0] + s = [0.0, 0.0] + g = Rank0Gaussian(A, s) + + w = [1.0, 0.0] + V1 = CoulombOperator(-1.0, w) + V2 = CoulombOperator(-2.0, w) + + result1 = _compute_matrix_element(g, g, V1) + result2 = _compute_matrix_element(g, g, V2) + + @test result2 ≈ 2 * result1 rtol = 1.0e-10 + end + + @testset "Different w vectors" begin + A = [1.0 0.0; 0.0 1.0] + s = [0.0, 0.0] + g = Rank0Gaussian(A, s) + + w1 = [1.0, 0.0] + w2 = [0.0, 1.0] + + V1 = CoulombOperator(-1.0, w1) + V2 = CoulombOperator(-1.0, w2) + + result1 = _compute_matrix_element(g, g, V1) + result2 = _compute_matrix_element(g, g, V2) + + @test result1 ≈ result2 rtol = 1.0e-10 + end +end - term2 = -2 * sqrt(β / π) * β / 3 * 0.25 * ( - q2_1 + q2_2 + q2_3 + q2_4 + q2_5 + q2_6 - ) * M0 +@testset "Error Cases" begin - term3 = 2 * sqrt(β / π) * β^2 / 10 * 0.5 * ( - q4_1 + q4_2 + q4_3 - ) * M0 + @testset "Dimension mismatch" begin + A1 = [1.0;;] # 1D + A2 = [1.0 0.0; 0.0 1.0] # 2D + s1 = [0.0] + s2 = [0.0, 0.0] - expected = coefficient * (term1 + term2 + term3) + g1 = Rank0Gaussian(A1, s1) + g2 = Rank0Gaussian(A2, s2) - @test isapprox(result, expected; atol = 1.0e-10) + @test_throws Exception _compute_matrix_element(g1, g2) + end end diff --git a/test/test_sampling.jl b/test/test_sampling.jl index f8c4041..ba61021 100644 --- a/test/test_sampling.jl +++ b/test/test_sampling.jl @@ -3,23 +3,5 @@ using FewBodyECG @testset "Sampling Module Tests" begin - @testset "generate_bij function" begin - b1 = 1.5 - n_terms = 5 - i = 2 - - bij_quasi = generate_bij(:quasirandom, i, n_terms, b1) - @test length(bij_quasi) == n_terms - @test all(0 .<= bij_quasi .<= b1) - using QuasiMonteCarlo - expected_quasi = QuasiMonteCarlo.sample(i + 1, n_terms, HaltonSample())[:, end] * b1 - @test isapprox(bij_quasi, expected_quasi; atol = 1.0e-12) - - bij_rand = generate_bij(:random, i, n_terms, b1) - @test length(bij_rand) == n_terms - @test all(0 .<= bij_rand .<= b1) - - @test_throws ErrorException generate_bij(:foo, i, n_terms, b1) - end end From 1c8a7eba1fa328670ebd645be6abc762f888fbcc Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Wed, 19 Nov 2025 10:39:17 +0100 Subject: [PATCH 08/15] updated documentation --- Examples/HydrogenAnion.jl | 1 - docs/Manifest.toml | 288 +++++++++++++++++++------------------- docs/Project.toml | 1 - docs/src/examples.md | 56 ++------ docs/src/index.md | 62 ++------ 5 files changed, 171 insertions(+), 237 deletions(-) diff --git a/Examples/HydrogenAnion.jl b/Examples/HydrogenAnion.jl index 77b6d98..56127c4 100644 --- a/Examples/HydrogenAnion.jl +++ b/Examples/HydrogenAnion.jl @@ -2,7 +2,6 @@ using FewBodyECG using LinearAlgebra using Plots using QuasiMonteCarlo -using FewBodyDB masses = [1.0e15, 1.0, 1.0] diff --git a/docs/Manifest.toml b/docs/Manifest.toml index 2e5e8bc..8489179 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -1,13 +1,13 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.11.7" +julia_version = "1.12.1" manifest_format = "2.0" -project_hash = "7bc105c61ab9f45572602f700a9e0c9169cae73d" +project_hash = "c701bdbde2b173cbc9ce030ccf07a0026a315dfc" [[deps.ADTypes]] -git-tree-sha1 = "27cecae79e5cc9935255f90c53bb831cc3c870d7" +git-tree-sha1 = "8be2ae325471fc20b11c27bb34b518541d07dd3a" uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b" -version = "1.18.0" +version = "1.19.0" [deps.ADTypes.extensions] ADTypesChainRulesCoreExt = "ChainRulesCore" @@ -55,9 +55,9 @@ version = "0.1.42" [[deps.Adapt]] deps = ["LinearAlgebra", "Requires"] -git-tree-sha1 = "f7817e2e585aa6d924fd714df1e2a84be7896c60" +git-tree-sha1 = "7e35fca2bdfba44d797c53dfe63a51fabf39bfc0" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "4.3.0" +version = "4.4.0" [deps.Adapt.extensions] AdaptSparseArraysExt = "SparseArrays" @@ -79,15 +79,15 @@ version = "1.1.2" [[deps.ArrayInterface]] deps = ["Adapt", "LinearAlgebra"] -git-tree-sha1 = "dbd8c3bbbdbb5c2778f85f4422c39960eac65a42" +git-tree-sha1 = "d81ae5489e13bc03567d4fbbb06c546a5e53c857" uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "7.20.0" +version = "7.22.0" [deps.ArrayInterface.extensions] ArrayInterfaceBandedMatricesExt = "BandedMatrices" ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices" ArrayInterfaceCUDAExt = "CUDA" - ArrayInterfaceCUDSSExt = "CUDSS" + ArrayInterfaceCUDSSExt = ["CUDSS", "CUDA"] ArrayInterfaceChainRulesCoreExt = "ChainRulesCore" ArrayInterfaceChainRulesExt = "ChainRules" ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore" @@ -182,9 +182,9 @@ version = "0.3.1" [[deps.Compat]] deps = ["TOML", "UUIDs"] -git-tree-sha1 = "0037835448781bb46feb39866934e243886d756a" +git-tree-sha1 = "9d8a54ce4b17aa5bdce0ea5c34bc5e7c340d16ad" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.18.0" +version = "4.18.1" weakdeps = ["Dates", "LinearAlgebra"] [deps.Compat.extensions] @@ -193,7 +193,7 @@ weakdeps = ["Dates", "LinearAlgebra"] [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.1.1+0" +version = "1.3.0+1" [[deps.CompositionsBase]] git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad" @@ -242,9 +242,9 @@ version = "1.16.0" [[deps.DataStructures]] deps = ["OrderedCollections"] -git-tree-sha1 = "6c72198e6a101cccdd4c9731d3985e904ba26037" +git-tree-sha1 = "e357641bb3e0638d353c4b29ea0e40ea644066a6" uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.19.1" +version = "0.19.3" [[deps.Dates]] deps = ["Printf"] @@ -277,9 +277,9 @@ version = "1.15.1" [[deps.DifferentiationInterface]] deps = ["ADTypes", "LinearAlgebra"] -git-tree-sha1 = "16946a4d305607c3a4af54ff35d56f0e9444ed0e" +git-tree-sha1 = "80bd15222b3e8d0bc70d921d2201aa0084810ce5" uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" -version = "0.7.7" +version = "0.7.12" [deps.DifferentiationInterface.extensions] DifferentiationInterfaceChainRulesCoreExt = "ChainRulesCore" @@ -337,9 +337,9 @@ version = "0.9.5" [[deps.Documenter]] deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "CodecZlib", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "TOML", "Test", "Unicode"] -git-tree-sha1 = "47ffb8f27ffc01e2e57e7ae5365ae5ceef87b03d" +git-tree-sha1 = "70c521ca3a23c576e12655d15963977c9766c26b" uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -version = "1.14.1" +version = "1.16.0" [[deps.Downloads]] deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] @@ -365,27 +365,21 @@ version = "0.1.11" [[deps.Expat_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "7bb1361afdb33c7f2b085aa49ea8fe1b0fb14e58" +git-tree-sha1 = "27af30de8b5445644e8ffe3bcb0d72049c089cf1" uuid = "2e619515-83b5-522b-bb60-26c02a35a201" -version = "2.7.1+0" +version = "2.7.3+0" [[deps.FFMPEG]] deps = ["FFMPEG_jll"] -git-tree-sha1 = "83dc665d0312b41367b7263e8a4d172eac1897f4" +git-tree-sha1 = "95ecf07c2eea562b5adbd0696af6db62c0f52560" uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" -version = "0.4.4" +version = "0.4.5" [[deps.FFMPEG_jll]] deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"] -git-tree-sha1 = "3a948313e7a41eb1db7a1e733e6335f17b4ab3c4" +git-tree-sha1 = "ccc81ba5e42497f4e76553a5545665eed577a663" uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5" -version = "7.1.1+0" - -[[deps.FewBodyECG]] -deps = ["FewBodyHamiltonians", "LinearAlgebra", "QuasiMonteCarlo", "SpecialFunctions"] -path = ".." -uuid = "083b1810-24a1-4a79-9a41-145bb2bb8ceb" -version = "1.0.5" +version = "8.0.0+0" [[deps.FewBodyHamiltonians]] git-tree-sha1 = "3cf35661914b5eb9bce1f2a5ced704b4133ee1d4" @@ -398,9 +392,9 @@ version = "1.11.0" [[deps.FillArrays]] deps = ["LinearAlgebra"] -git-tree-sha1 = "173e4d8f14230a7523ae11b9a3fa9edb3e0efd78" +git-tree-sha1 = "5bfcd42851cf2f1b303f51525a54dc5e98d408a3" uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" -version = "1.14.0" +version = "1.15.0" [deps.FillArrays.extensions] FillArraysPDMatsExt = "PDMats" @@ -414,9 +408,9 @@ version = "1.14.0" [[deps.FiniteDiff]] deps = ["ArrayInterface", "LinearAlgebra", "Setfield"] -git-tree-sha1 = "31fd32af86234b6b71add76229d53129aa1b87a9" +git-tree-sha1 = "9340ca07ca27093ff68418b7558ca37b05f8aeb1" uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" -version = "2.28.1" +version = "2.29.0" [deps.FiniteDiff.extensions] FiniteDiffBandedMatricesExt = "BandedMatrices" @@ -449,9 +443,9 @@ version = "1.3.7" [[deps.ForwardDiff]] deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"] -git-tree-sha1 = "dc41303865a16274ecb8450c220021ce1e0cf05f" +git-tree-sha1 = "cd33c7538e68650bd0ddbb3f5bd50a4a0fa95b50" uuid = "f6369f11-7733-5829-9624-2563aa707210" -version = "1.2.1" +version = "1.3.0" [deps.ForwardDiff.extensions] ForwardDiffStaticArraysExt = "StaticArrays" @@ -484,15 +478,21 @@ version = "3.4.0+2" [[deps.GR]] deps = ["Artifacts", "Base64", "DelimitedFiles", "Downloads", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Preferences", "Printf", "Qt6Wayland_jll", "Random", "Serialization", "Sockets", "TOML", "Tar", "Test", "p7zip_jll"] -git-tree-sha1 = "1828eb7275491981fa5f1752a5e126e8f26f8741" +git-tree-sha1 = "f305bdb91e1f3fcc687944c97f2ede40585b1bd5" uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" -version = "0.73.17" +version = "0.73.19" + + [deps.GR.extensions] + GRIJuliaExt = "IJulia" + + [deps.GR.weakdeps] + IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" [[deps.GR_jll]] deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "FreeType2_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Qt6Base_jll", "Zlib_jll", "libpng_jll"] -git-tree-sha1 = "27299071cc29e409488ada41ec7643e0ab19091f" +git-tree-sha1 = "de439fbc02b9dc0e639e67d7c5bd5811ff3b6f06" uuid = "d2c73de3-f751-5644-a686-071e5b155ba9" -version = "0.73.17+0" +version = "0.73.19+1" [[deps.GettextRuntime_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll"] @@ -520,9 +520,9 @@ version = "3.7.0+0" [[deps.Git_jll]] deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"] -git-tree-sha1 = "e2aef26f7d273f1e5b1daba56837c47b49b4388f" +git-tree-sha1 = "b6a684587ebe896d9f68ae777f648205940f0f70" uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb" -version = "2.51.1+0" +version = "2.51.3+0" [[deps.Glib_jll]] deps = ["Artifacts", "GettextRuntime_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Zlib_jll"] @@ -543,9 +543,9 @@ version = "1.0.2" [[deps.HTTP]] deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "PrecompileTools", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"] -git-tree-sha1 = "ed5e9c58612c4e081aecdb6e1a479e18462e041e" +git-tree-sha1 = "5e6fe50ae7f23d171f44e311c2960294aaa0beb5" uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" -version = "1.10.17" +version = "1.10.19" [[deps.HarfBuzz_jll]] deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll"] @@ -555,9 +555,9 @@ version = "8.5.1+0" [[deps.IOCapture]] deps = ["Logging", "Random"] -git-tree-sha1 = "b6d6bfdd7ce25b0f9b2f6b3dd56b2673a66c8770" +git-tree-sha1 = "0ee181ec08df7d7c911901ea38baf16f755114dc" uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89" -version = "0.2.5" +version = "1.0.0" [[deps.IntegerMathUtils]] git-tree-sha1 = "4c1acff2dc6b6967e7e750633c50bc3b8d83e617" @@ -580,9 +580,9 @@ weakdeps = ["Dates", "Test"] InverseFunctionsTestExt = "Test" [[deps.IrrationalConstants]] -git-tree-sha1 = "e2222959fbc6c19554dc15174c81bf7bf3aa691c" +git-tree-sha1 = "b2d91fe939cae05960e760110b328288867b5758" uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" -version = "0.2.4" +version = "0.2.6" [[deps.JLFzf]] deps = ["REPL", "Random", "fzf_jll"] @@ -597,10 +597,16 @@ uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" version = "1.7.1" [[deps.JSON]] -deps = ["Dates", "Mmap", "Parsers", "Unicode"] -git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a" +deps = ["Dates", "Logging", "Parsers", "PrecompileTools", "StructUtils", "UUIDs", "Unicode"] +git-tree-sha1 = "5b6bb73f555bc753a6153deec3717b8904f5551c" uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -version = "0.21.4" +version = "1.3.0" + + [deps.JSON.extensions] + JSONArrowExt = ["ArrowTypes"] + + [deps.JSON.weakdeps] + ArrowTypes = "31f734f8-188a-4ce0-8406-c8a06bd891cd" [[deps.JpegTurbo_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -608,6 +614,11 @@ git-tree-sha1 = "4255f0032eafd6451d707a51d5f0248b8a165e4d" uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8" version = "3.1.3+0" +[[deps.JuliaSyntaxHighlighting]] +deps = ["StyledStrings"] +uuid = "ac6e5ff7-fb65-4e79-a425-ec3bc9c03011" +version = "1.12.0" + [[deps.LAME_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] git-tree-sha1 = "059aabebaa7c82ccb853dd4a0ee9d17796f7e1bc" @@ -672,24 +683,24 @@ uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" version = "0.6.4" [[deps.LibCURL_jll]] -deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "OpenSSL_jll", "Zlib_jll", "nghttp2_jll"] uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" -version = "8.6.0+0" +version = "8.11.1+1" [[deps.LibGit2]] -deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"] +deps = ["LibGit2_jll", "NetworkOptions", "Printf", "SHA"] uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" version = "1.11.0" [[deps.LibGit2_jll]] -deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "OpenSSL_jll"] uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5" -version = "1.7.2+0" +version = "1.9.0+0" [[deps.LibSSH2_jll]] -deps = ["Artifacts", "Libdl", "MbedTLS_jll"] +deps = ["Artifacts", "Libdl", "OpenSSL_jll"] uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" -version = "1.11.0+1" +version = "1.11.3+1" [[deps.Libdl]] uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" @@ -715,32 +726,32 @@ version = "1.18.0+0" [[deps.Libmount_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "706dfd3c0dd56ca090e86884db6eda70fa7dd4af" +git-tree-sha1 = "3acf07f130a76f87c041cfb2ff7d7284ca67b072" uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9" -version = "2.41.1+0" +version = "2.41.2+0" [[deps.Libtiff_jll]] deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "LERC_jll", "Libdl", "XZ_jll", "Zlib_jll", "Zstd_jll"] -git-tree-sha1 = "4ab7581296671007fc33f07a721631b8855f4b1d" +git-tree-sha1 = "f04133fe05eff1667d2054c53d59f9122383fe05" uuid = "89763e89-9b03-5906-acba-b20f662cd828" -version = "4.7.1+0" +version = "4.7.2+0" [[deps.Libuuid_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "d3c8af829abaeba27181db4acb485b18d15d89c6" +git-tree-sha1 = "2a7a12fc0a4e7fb773450d17975322aa77142106" uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700" -version = "2.41.1+0" +version = "2.41.2+0" [[deps.LineSearches]] deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"] -git-tree-sha1 = "4adee99b7262ad2a1a4bbbc59d993d24e55ea96f" +git-tree-sha1 = "a8b1215fb05581a1f9e403bec46a1333e7eb1ffb" uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" -version = "7.4.0" +version = "7.4.1" [[deps.LinearAlgebra]] deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -version = "1.11.0" +version = "1.12.0" [[deps.LogExpFunctions]] deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] @@ -764,9 +775,9 @@ version = "1.11.0" [[deps.LoggingExtras]] deps = ["Dates", "Logging"] -git-tree-sha1 = "f02b56007b064fbfddb4c9cd60161b6dd0f40df3" +git-tree-sha1 = "f00544d95982ea270145636c181ceda21c4e2575" uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36" -version = "1.1.0" +version = "1.2.0" [[deps.MacroTools]] git-tree-sha1 = "1e0228a030642014fe5cfe68c2c0a818f9e3f522" @@ -774,7 +785,7 @@ uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" version = "0.5.16" [[deps.Markdown]] -deps = ["Base64"] +deps = ["Base64", "JuliaSyntaxHighlighting", "StyledStrings"] uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" version = "1.11.0" @@ -791,14 +802,15 @@ uuid = "739be429-bea8-5141-9913-cc70e7f3736d" version = "1.1.9" [[deps.MbedTLS_jll]] -deps = ["Artifacts", "Libdl"] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "3cce3511ca2c6f87b19c34ffc623417ed2798cbd" uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.28.6+0" +version = "2.28.10+0" [[deps.Measures]] -git-tree-sha1 = "c13304c81eec1ed3af7fc20e75fb6b26092a1102" +git-tree-sha1 = "b513cedd20d9c914783d8ad83d08120702bf2c77" uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e" -version = "0.3.2" +version = "0.3.3" [[deps.Missings]] deps = ["DataAPI"] @@ -812,7 +824,7 @@ version = "1.11.0" [[deps.MozillaCACerts_jll]] uuid = "14a3606d-f60d-562e-9121-12d972cd8159" -version = "2023.12.12" +version = "2025.5.20" [[deps.NLSolversBase]] deps = ["ADTypes", "DifferentiationInterface", "Distributed", "FiniteDiff", "ForwardDiff"] @@ -828,7 +840,7 @@ version = "1.1.3" [[deps.NetworkOptions]] uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" -version = "1.2.0" +version = "1.3.0" [[deps.Ogg_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -839,30 +851,29 @@ version = "1.3.6+0" [[deps.OpenBLAS_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" -version = "0.3.27+1" +version = "0.3.29+0" [[deps.OpenLibm_jll]] deps = ["Artifacts", "Libdl"] uuid = "05823500-19ac-5b8b-9628-191a04bc5112" -version = "0.8.5+0" +version = "0.8.7+0" [[deps.OpenSSH_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "OpenSSL_jll", "Zlib_jll"] -git-tree-sha1 = "1f2f0911e1c02f28a390bb720f97f3349c4dcefb" +git-tree-sha1 = "301412a644646fdc0ad67d0a87487466b491e53d" uuid = "9bd350c2-7e96-507f-8002-3f2e150b4e1b" -version = "10.0.2+0" +version = "10.2.1+0" [[deps.OpenSSL]] -deps = ["BitFlags", "Dates", "MozillaCACerts_jll", "OpenSSL_jll", "Sockets"] -git-tree-sha1 = "f1a7e086c677df53e064e0fdd2c9d0b0833e3f6e" +deps = ["BitFlags", "Dates", "MozillaCACerts_jll", "NetworkOptions", "OpenSSL_jll", "Sockets"] +git-tree-sha1 = "386b47442468acfb1add94bf2d85365dea10cbab" uuid = "4d8831e6-92b7-49fb-bdf8-b643e874388c" -version = "1.5.0" +version = "1.6.0" [[deps.OpenSSL_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "2ae7d4ddec2e13ad3bddf5c0796f7547cf682391" +deps = ["Artifacts", "Libdl"] uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" -version = "3.5.2+0" +version = "3.5.1+0" [[deps.OpenSpecFun_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"] @@ -896,13 +907,13 @@ version = "1.8.1" [[deps.PCRE2_jll]] deps = ["Artifacts", "Libdl"] uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15" -version = "10.42.0+1" +version = "10.44.0+1" [[deps.Pango_jll]] deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "FriBidi_jll", "Glib_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl"] -git-tree-sha1 = "1f7f9bbd5f7a2e5a9f7d96e51c9754454ea7f60b" +git-tree-sha1 = "0662b083e11420952f2e62e17eddae7fc07d5997" uuid = "36c8627f-9965-5494-a995-c6b170f724f3" -version = "1.56.4+0" +version = "1.57.0+0" [[deps.Parameters]] deps = ["OrderedCollections", "UnPack"] @@ -925,7 +936,7 @@ version = "0.44.2+0" [[deps.Pkg]] deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "Random", "SHA", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.11.0" +version = "1.12.0" weakdeps = ["REPL"] [deps.Pkg.extensions] @@ -939,15 +950,15 @@ version = "3.3.0" [[deps.PlotUtils]] deps = ["ColorSchemes", "Colors", "Dates", "PrecompileTools", "Printf", "Random", "Reexport", "StableRNGs", "Statistics"] -git-tree-sha1 = "3ca9a356cd2e113c420f2c13bea19f8d3fb1cb18" +git-tree-sha1 = "26ca162858917496748aad52bb5d3be4d26a228a" uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043" -version = "1.4.3" +version = "1.4.4" [[deps.Plots]] -deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "JLFzf", "JSON", "LaTeXStrings", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "Pkg", "PlotThemes", "PlotUtils", "PrecompileTools", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "RelocatableFolders", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "TOML", "UUIDs", "UnicodeFun", "UnitfulLatexify", "Unzip"] -git-tree-sha1 = "bfe839e9668f0c58367fb62d8757315c0eac8777" +deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "JLFzf", "JSON", "LaTeXStrings", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "Pkg", "PlotThemes", "PlotUtils", "PrecompileTools", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "RelocatableFolders", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "TOML", "UUIDs", "UnicodeFun", "Unzip"] +git-tree-sha1 = "12ce661880f8e309569074a61d3767e5756a199f" uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -version = "1.40.20" +version = "1.41.1" [deps.Plots.extensions] FileIOExt = "FileIO" @@ -971,9 +982,9 @@ version = "0.2.4" [[deps.PrecompileTools]] deps = ["Preferences"] -git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f" +git-tree-sha1 = "07a921781cab75691315adc645096ed5e370cb77" uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" -version = "1.2.1" +version = "1.3.3" [[deps.Preferences]] deps = ["TOML"] @@ -999,9 +1010,9 @@ version = "1.3.0" [[deps.Qt6Base_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Vulkan_Loader_jll", "Xorg_libSM_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_cursor_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "libinput_jll", "xkbcommon_jll"] -git-tree-sha1 = "eb38d376097f47316fe089fc62cb7c6d85383a52" +git-tree-sha1 = "34f7e5d2861083ec7596af8b8c092531facf2192" uuid = "c0090381-4147-56d7-9ebc-da0b1113ec56" -version = "6.8.2+1" +version = "6.8.2+2" [[deps.Qt6Declarative_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Qt6Base_jll", "Qt6ShaderTools_jll"] @@ -1017,9 +1028,9 @@ version = "6.8.2+1" [[deps.Qt6Wayland_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Qt6Base_jll", "Qt6Declarative_jll"] -git-tree-sha1 = "e1d5e16d0f65762396f9ca4644a5f4ddab8d452b" +git-tree-sha1 = "8f528b0851b5b7025032818eb5abbeb8a736f853" uuid = "e99dba38-086e-5de3-a5b1-6e4c66e897c3" -version = "6.8.2+1" +version = "6.8.2+2" [[deps.QuasiMonteCarlo]] deps = ["Accessors", "ConcreteStructs", "LatticeRules", "LinearAlgebra", "Primes", "Random", "Requires", "Sobol", "StatsBase"] @@ -1034,7 +1045,7 @@ version = "0.3.3" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" [[deps.REPL]] -deps = ["InteractiveUtils", "Markdown", "Sockets", "StyledStrings", "Unicode"] +deps = ["InteractiveUtils", "JuliaSyntaxHighlighting", "Markdown", "Sockets", "StyledStrings", "Unicode"] uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" version = "1.11.0" @@ -1128,13 +1139,13 @@ version = "1.2.2" [[deps.SparseArrays]] deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -version = "1.11.0" +version = "1.12.0" [[deps.SpecialFunctions]] deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] -git-tree-sha1 = "41852b8679f78c8d8961eeadc8f62cef861a52e3" +git-tree-sha1 = "f2685b435df2613e25fc10ad8c26dddb8640f547" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "2.5.1" +version = "2.6.1" [deps.SpecialFunctions.extensions] SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" @@ -1144,14 +1155,14 @@ version = "2.5.1" [[deps.StableRNGs]] deps = ["Random"] -git-tree-sha1 = "95af145932c2ed859b63329952ce8d633719f091" +git-tree-sha1 = "4f96c596b8c8258cc7d3b19797854d368f243ddc" uuid = "860ef19b-820b-49d6-a774-d7a799459cd3" -version = "1.0.3" +version = "1.0.4" [[deps.StaticArraysCore]] -git-tree-sha1 = "192954ef1208c7019899fbf8049e717f92959682" +git-tree-sha1 = "6ab403037779dae8c514bad259f32a447262455a" uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" -version = "1.4.3" +version = "1.4.4" [[deps.Statistics]] deps = ["LinearAlgebra"] @@ -1171,9 +1182,23 @@ version = "1.7.1" [[deps.StatsBase]] deps = ["AliasTables", "DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] -git-tree-sha1 = "2c962245732371acd51700dbb268af311bddd719" +git-tree-sha1 = "064b532283c97daae49e544bb9cb413c26511f8c" uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -version = "0.34.6" +version = "0.34.8" + +[[deps.StructUtils]] +deps = ["Dates", "UUIDs"] +git-tree-sha1 = "79529b493a44927dd5b13dde1c7ce957c2d049e4" +uuid = "ec057cc2-7a8d-4b58-b3b3-92acb9f63b42" +version = "2.6.0" + + [deps.StructUtils.extensions] + StructUtilsMeasurementsExt = ["Measurements"] + StructUtilsTablesExt = ["Tables"] + + [deps.StructUtils.weakdeps] + Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" + Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" [[deps.StyledStrings]] uuid = "f489334b-da3d-4c2e-b8f0-e476e12c162b" @@ -1182,7 +1207,7 @@ version = "1.11.0" [[deps.SuiteSparse_jll]] deps = ["Artifacts", "Libdl", "libblastrampoline_jll"] uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" -version = "7.7.0+0" +version = "7.8.3+2" [[deps.TOML]] deps = ["Dates"] @@ -1235,25 +1260,6 @@ git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf" uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1" version = "0.4.1" -[[deps.Unitful]] -deps = ["Dates", "LinearAlgebra", "Random"] -git-tree-sha1 = "6258d453843c466d84c17a58732dda5deeb8d3af" -uuid = "1986cc42-f94f-5a68-af5c-568840ba703d" -version = "1.24.0" -weakdeps = ["ConstructionBase", "ForwardDiff", "InverseFunctions", "Printf"] - - [deps.Unitful.extensions] - ConstructionBaseUnitfulExt = "ConstructionBase" - ForwardDiffExt = "ForwardDiff" - InverseFunctionsUnitfulExt = "InverseFunctions" - PrintfExt = "Printf" - -[[deps.UnitfulLatexify]] -deps = ["LaTeXStrings", "Latexify", "Unitful"] -git-tree-sha1 = "af305cc62419f9bd61b6644d19170a4d258c7967" -uuid = "45397f5d-5981-4c77-b2b3-fc36d6e9b728" -version = "1.7.0" - [[deps.Unzip]] git-tree-sha1 = "ca0969166a028236229f63514992fc073799bb78" uuid = "41fe7b60-77ed-43a1-b4f0-825fd5a5650d" @@ -1321,9 +1327,9 @@ version = "1.3.7+0" [[deps.Xorg_libXfixes_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll"] -git-tree-sha1 = "9caba99d38404b285db8801d5c45ef4f4f425a6d" +git-tree-sha1 = "75e00946e43621e09d431d9b95818ee751e6b2ef" uuid = "d091e8ba-531a-589c-9de9-94069b037ed8" -version = "6.0.1+0" +version = "6.0.2+0" [[deps.Xorg_libXi_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libXext_jll", "Xorg_libXfixes_jll"] @@ -1363,9 +1369,9 @@ version = "1.1.3+0" [[deps.Xorg_xcb_util_cursor_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_jll", "Xorg_xcb_util_renderutil_jll"] -git-tree-sha1 = "c5bf2dad6a03dfef57ea0a170a1fe493601603f2" +git-tree-sha1 = "9750dc53819eba4e9a20be42349a6d3b86c7cdf8" uuid = "e920d4aa-a673-5f3a-b3d7-f755a4d47c43" -version = "0.1.5+0" +version = "0.1.6+0" [[deps.Xorg_xcb_util_image_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_xcb_util_jll"] @@ -1418,7 +1424,7 @@ version = "1.6.0+0" [[deps.Zlib_jll]] deps = ["Libdl"] uuid = "83775a58-1f1d-513f-b197-d71354ab007a" -version = "1.2.13+1" +version = "1.3.1+2" [[deps.Zstd_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -1440,9 +1446,9 @@ version = "0.61.1+0" [[deps.libaom_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "4bba74fa59ab0755167ad24f98800fe5d727175b" +git-tree-sha1 = "371cc681c00a3ccc3fbc5c0fb91f58ba9bec1ecf" uuid = "a4ae2306-e953-59d6-aa16-d00cac43593b" -version = "3.12.1+0" +version = "3.13.1+0" [[deps.libass_jll]] deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Zlib_jll"] @@ -1453,7 +1459,7 @@ version = "0.17.4+0" [[deps.libblastrampoline_jll]] deps = ["Artifacts", "Libdl"] uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" -version = "5.11.0+0" +version = "5.15.0+0" [[deps.libdecor_jll]] deps = ["Artifacts", "Dbus_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pango_jll", "Wayland_jll", "xkbcommon_jll"] @@ -1500,12 +1506,12 @@ version = "1.1.7+0" [[deps.nghttp2_jll]] deps = ["Artifacts", "Libdl"] uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" -version = "1.59.0+0" +version = "1.64.0+1" [[deps.p7zip_jll]] deps = ["Artifacts", "Libdl"] uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" -version = "17.4.0+2" +version = "17.5.0+2" [[deps.x264_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] diff --git a/docs/Project.toml b/docs/Project.toml index 2e161fa..17b7dec 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,6 +1,5 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -FewBodyECG = "083b1810-24a1-4a79-9a41-145bb2bb8ceb" FewBodyHamiltonians = "3a126c26-e5d7-4a95-83c3-3b69f8a11ded" Optim = "429524aa-4258-5aef-a3af-852621145aeb" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" diff --git a/docs/src/examples.md b/docs/src/examples.md index 0561683..226d16b 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -9,53 +9,27 @@ using Plots using QuasiMonteCarlo masses = [1.0e15, 1.0, 1.0] -psys = ParticleSystem(masses) -K = Diagonal([0.0, 1 / 2, 1 / 2]) -K_transformed = psys.J * K * psys.J' +Λmat = Λ(masses) +kin = KineticOperator(Λmat) +J, U = _jacobi_transform(masses) w_list = [[1, -1, 0], [1, 0, -1], [0, 1, -1]] -w_raw = [psys.U' * w for w in w_list] +w_raw = [U' * w for w in w_list] +coeffs = [-1.0, -1.0, +1.0] -let - n_basis = 50 - b1 = default_b0(psys.scale) - method = :quasirandom - basis_fns = GaussianBase[] - E₀_list = Float64[] - coeffs = [-1.0, -1.0, +1.0] +ops = Operator[ + kin; + (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))... +] - for i in 1:n_basis - bij = generate_bij(method, i, length(w_raw), b1; qmc_sampler = SobolSample()) - A = _generate_A_matrix(bij, w_raw) - push!(basis_fns, Rank0Gaussian(A)) +result = solve_ECG(ops, 250, scale = 1.0) - basis = BasisSet(basis_fns) - ops = Operator[ - KineticOperator(K_transformed); - (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))... - ] +E = -0.527751016523 +ΔE = abs(result.ground_state - E) +@info "Energy difference" ΔE - H = build_hamiltonian_matrix(basis, ops) - S = build_overlap_matrix(basis) - vals, _ = solve_generalized_eigenproblem(H, S) - E₀_step = minimum(vals) - - push!(E₀_list, E₀_step) - println("Step $i: E₀ = $E₀_step") - - end - - E₀ = minimum(E₀_list) - Eᵗʰ = -0.527751016523 - ΔE = abs(E₀ - Eᵗʰ) - @show ΔE - - plot( - 1:n_basis, E₀_list, xlabel = "Number of Gaussians", ylabel = "E₀ [Hartree]", - lw = 2, label = "Ground state energy", title = "Hydrogen Anion Convergence" - ) - -end +n, E = convergence(result) +plot(n, E) ``` \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md index 7dd2cad..a529dac 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -26,64 +26,20 @@ using Plots using QuasiMonteCarlo masses = [1.0, 1.0, 1.0] -psys = ParticleSystem(masses) -K = Diagonal([1 / 2, 1 / 2, 1 / 2]) -K_transformed = psys.J * K * psys.J' +Λmat = Λ(masses) +kin = KineticOperator(Λmat) +J, U = _jacobi_transform(masses) w_list = [[1, -1, 0], [1, 0, -1], [0, 1, -1]] -w_raw = [psys.U' * w for w in w_list] +w_raw = [U' * w for w in w_list] coeffs = [+1.0, -1.0, -1.0] +coulomb_ops = [CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw)] -let - n_basis = 100 - b1 = default_b0(psys.scale) - method = :quasirandom - basis_fns = GaussianBase[] - E₀_list = Float64[] +ops = Operator[kin; coulomb_ops...] - for i in 1:n_basis - bij = generate_bij(method, i, length(w_raw), b1; qmc_sampler = SobolSample()) - A = _generate_A_matrix(bij, w_raw) - push!(basis_fns, Rank0Gaussian(A)) - - basis = BasisSet(basis_fns) - ops = Operator[ - KineticOperator(K_transformed); - (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))... - ] - - H = build_hamiltonian_matrix(basis, ops) - S = build_overlap_matrix(basis) - vals, vecs = solve_generalized_eigenproblem(H, S) - global c₀ = vecs[:, 1] - E₀ = minimum(vals) - - push!(E₀_list, E₀) - println("Step $i: E₀ = $E₀") - end - - E₀ = minimum(E₀_list) - Eᵗʰ = -0.2620050702328 - ΔE = abs(E₀ - Eᵗʰ) - @show ΔE - - r = range(0.01, 14.0, length = 400) - ρ_r = [rval^2 * abs2(ψ₀([rval, 0.0], c₀, basis_fns)) for rval in r] - - p1 = plot( - r, ρ_r, xlabel = "r (a.u.)", ylabel = "r²|ψ₀(r)|²", - lw = 2, label = "r²C(r)", title = "Electron-Positron Correlation Function" - ) - - - p2 = plot( - 1:n_basis, E₀_list, xlabel = "Number of Gaussians", ylabel = "E₀ [Hartree]", - lw = 2, label = "Ground state energy", title = "Positronium Convergence" - ) - - plot(p1, p2, layout = (2, 1)) - -end +result = solve_ECG(ops, 250, sampler = SobolSample(); scale = 0.2) +n, E = convergence(result) +plot(n, E) ``` \ No newline at end of file From cecb2b7f9d5b3b740787a9ef0b681b879bb94eb5 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Wed, 19 Nov 2025 13:44:26 +0100 Subject: [PATCH 09/15] updated tests --- src/hamiltonian.jl | 6 -- test/test_coordinates.jl | 8 -- test/test_hamiltonian.jl | 203 +++++++++++++++++++++++++++++++---- test/test_matrix_elements.jl | 68 ------------ test/test_sampling.jl | 43 ++++++++ 5 files changed, 223 insertions(+), 105 deletions(-) diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl index 421086f..7cf8e22 100644 --- a/src/hamiltonian.jl +++ b/src/hamiltonian.jl @@ -45,12 +45,6 @@ function solve_generalized_eigenproblem(H::AbstractMatrix{<:Real}, S::AbstractMa return real(evals), real(vecs) end -function diagonalize(basis::BasisSet{<:GaussianBase}, operators::AbstractVector{<:FewBodyHamiltonians.Operator}) - H = build_hamiltonian_matrix(basis, operators) - S = build_overlap_matrix(basis) - return solve_generalized_eigenproblem(H, S) -end - function normalized_overlap(A::GaussianBase, B::GaussianBase) overlap_12 = _compute_matrix_element(A, B) overlap_11 = _compute_matrix_element(A, A) diff --git a/test/test_coordinates.jl b/test/test_coordinates.jl index e82dc52..a2d7000 100644 --- a/test/test_coordinates.jl +++ b/test/test_coordinates.jl @@ -5,7 +5,6 @@ import FewBodyECG: _jacobi_transform, _generate_A_matrix, _shift_vectors, _trans @testset "Coordinates Module Tests" begin @testset "_jacobi_transform" begin - # Test with simple equal masses masses = [1.0, 1.0, 1.0] J, U = _jacobi_transform(masses) @@ -13,7 +12,6 @@ import FewBodyECG: _jacobi_transform, _generate_A_matrix, _shift_vectors, _trans @test size(U) == (3, 2) @test J * U ≈ I(2) atol = 1.0e-10 - # Test with different masses masses = [1.0, 2.0, 3.0] J, U = _jacobi_transform(masses) @@ -21,7 +19,6 @@ import FewBodyECG: _jacobi_transform, _generate_A_matrix, _shift_vectors, _trans @test size(U) == (3, 2) @test J * U ≈ I(2) atol = 1.0e-10 - # Test with exactly two masses masses = [1.0, 2.0] J, U = _jacobi_transform(masses) @@ -29,7 +26,6 @@ import FewBodyECG: _jacobi_transform, _generate_A_matrix, _shift_vectors, _trans @test size(U) == (2, 1) @test J * U ≈ [1.0] atol = 1.0e-10 - # Test error for fewer than two masses @test_throws AssertionError _jacobi_transform([1.0]) end @@ -45,11 +41,9 @@ import FewBodyECG: _jacobi_transform, _generate_A_matrix, _shift_vectors, _trans r_back = _inverse_transform_coordinates(U, x) @test size(r_back) == (3,) - # Instead of round-trip r → x → r_back, test projection recovery x_back = _transform_coordinates(J, r_back) @test x_back ≈ x atol = 1.0e-10 - # Error case @test_throws AssertionError _transform_coordinates(J, [1.0, 2.0]) @test_throws AssertionError _inverse_transform_coordinates(U, [1.0]) end @@ -87,7 +81,6 @@ end masses = [1.3, 2.5, 0.7, 4.1] J, U = _jacobi_transform(masses) - # J * U should act like the identity on the reduced space Ired = I(size(J, 1)) @test isapprox(J * U, Ired; atol = 1.0e-10) @@ -95,7 +88,6 @@ end @test isapprox(J * U * J, J; atol = 1.0e-10) @test isapprox(U * J * U, U; atol = 1.0e-10) - # Symmetry conditions @test isapprox((J * U)', J * U; atol = 1.0e-10) @test isapprox((U * J)', U * J; atol = 1.0e-10) end diff --git a/test/test_hamiltonian.jl b/test/test_hamiltonian.jl index fd15f95..0f6a44f 100644 --- a/test/test_hamiltonian.jl +++ b/test/test_hamiltonian.jl @@ -2,35 +2,192 @@ using Test using LinearAlgebra using FewBodyHamiltonians using FewBodyECG -import FewBodyECG: _compute_overlap_element, _build_operator_matrix, _compute_matrix_element +import FewBodyECG: _compute_overlap_element, _build_operator_matrix, _compute_matrix_element, normalized_overlap, is_linearly_independent -@testset "hamiltonian extra tests" begin - A1 = reshape([1.0], 1, 1) - A2 = reshape([1.0], 1, 1) - A3 = reshape([1.0], 1, 1) - g1 = Rank0Gaussian(A1, [10.0]) - g2 = Rank0Gaussian(A2, [20.0]) - g3 = Rank0Gaussian(A3, [30.0]) - basis3 = BasisSet{Rank0Gaussian}([g1, g2, g3]) +@testset "build_overlap_matrix" begin + + @testset "Size and symmetry" begin + g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + g2 = Rank0Gaussian([1.5 0.0; 0.0 1.5], [0.0, 0.0]) + g3 = Rank0Gaussian([2.0 0.0; 0.0 2.0], [0.1, 0.1]) + + basis = BasisSet([g1, g2, g3]) + S = build_overlap_matrix(basis) + + @test size(S) == (3, 3) + @test issymmetric(S) + @test all(isfinite, S) + end + + + @testset "Consistency" begin + g1 = Rank0Gaussian([1.0;;], [0.0]) + g2 = Rank0Gaussian([2.0;;], [0.0]) + + basis = BasisSet([g1, g2]) + S = build_overlap_matrix(basis) + + S11 = _compute_matrix_element(g1, g1) + S12 = _compute_matrix_element(g1, g2) + S22 = _compute_matrix_element(g2, g2) + + @test S[1,1] ≈ S11 rtol=1e-12 + @test S[1,2] ≈ S12 rtol=1e-12 + @test S[2,2] ≈ S22 rtol=1e-12 + end +end - @eval FewBodyECG begin - function _compute_matrix_element(b::Rank0Gaussian, k::Rank0Gaussian) - return (b.s[1] + k.s[1]) / 10.0 - end - function _compute_matrix_element(b::Rank0Gaussian, k::Rank0Gaussian, op::FewBodyHamiltonians.Operator) - return (b.s[1] * k.s[1]) / 10.0 - end +@testset "build_hamiltonian_matrix" begin + + @testset "Size and symmetry" begin + g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + g2 = Rank0Gaussian([2.0 0.0; 0.0 2.0], [0.0, 0.0]) + basis = BasisSet([g1, g2]) + + K = KineticOperator([0.5 0.0; 0.0 0.5]) + H = build_hamiltonian_matrix(basis, [K]) + + @test size(H) == (2, 2) + @test issymmetric(H) + @test all(isfinite, H) + end + + @testset "Additivity" begin + # H(K+V) should equal H(K) + H(V) + g = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + basis = BasisSet([g]) + + K = KineticOperator([0.5 0.0; 0.0 0.5]) + V = CoulombOperator(-1.0, [1.0, 0.0]) + + H_K = build_hamiltonian_matrix(basis, [K]) + H_V = build_hamiltonian_matrix(basis, [V]) + H_both = build_hamiltonian_matrix(basis, [K, V]) + + @test H_both ≈ H_K + H_V rtol=1e-12 end + + @testset "Empty operators" begin + g = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + basis = BasisSet([g]) + + H = build_hamiltonian_matrix(basis, FewBodyHamiltonians.Operator[]) + + @test H == zeros(1, 1) + end +end + - S3 = build_overlap_matrix(basis3) - manualS = zeros(Float64, 3, 3) - for i in 1:3, j in 1:3 - manualS[i, j] = _compute_overlap_element(basis3.functions[i], basis3.functions[j]) +@testset "solve_generalized_eigenproblem" begin + + @testset "Basic solve" begin + H = [2.0 0.5; 0.5 3.0] + S = [1.0 0.1; 0.1 1.0] + + evals, evecs = solve_generalized_eigenproblem(H, S) + + @test length(evals) == 2 + @test size(evecs) == (2, 2) + @test all(isfinite, evals) + @test all(isfinite, evecs) + end + + @testset "Eigenvalue equation H*v = λ*S*v" begin + H = [3.0 0.5; 0.5 2.0] + S = [1.0 0.2; 0.2 1.0] + + evals, evecs = solve_generalized_eigenproblem(H, S) + + # Check each eigenpair + for i in 1:2 + λ = evals[i] + v = evecs[:, i] + residual = H * v - λ * S * v + @test norm(residual) < 1e-9 + end end - @test S3 == manualS - @test issymmetric(S3) - @test eltype(S3) == Float64 + + @testset "Requires positive definite S" begin + H = [1.0 0.0; 0.0 1.0] + S_bad = [-1.0 0.0; 0.0 1.0] + + @test_throws Exception solve_generalized_eigenproblem(H, S_bad) + end +end +@testset "normalized_overlap" begin + + @testset "Symmetry" begin + g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + g2 = Rank0Gaussian([2.0 0.0; 0.0 2.0], [0.5, 0.5]) + + overlap_12 = normalized_overlap(g1, g2) + overlap_21 = normalized_overlap(g2, g1) + + @test overlap_12 ≈ overlap_21 rtol=1e-12 + end + + @testset "Range [0, 1]" begin + g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + g2 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.5, 0.5]) + g3 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [5.0, 5.0]) + + overlap_12 = normalized_overlap(g1, g2) + overlap_13 = normalized_overlap(g1, g3) + + @test 0.0 <= overlap_12 <= 1.0 + @test 0.0 <= overlap_13 <= 1.0 + end + + @testset "Identical Gaussians" begin + A = [1.5 0.0; 0.0 1.5] + s = [0.3, 0.4] + g1 = Rank0Gaussian(A, s) + g2 = Rank0Gaussian(A, s) + + overlap = normalized_overlap(g1, g2) + @test overlap ≈ 1.0 rtol=1e-10 + end +end + +@testset "is_linearly_independent" begin + + @testset "Empty basis (always independent)" begin + g = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + empty_basis = BasisSet(Rank0Gaussian[]) + + @test is_linearly_independent(g, empty_basis; threshold=0.95) + end + + + @testset "Well-separated (should accept)" begin + g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + g2 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [5.0, 5.0]) + + basis = BasisSet([g1]) + + @test is_linearly_independent(g2, basis; threshold=0.95) + end + + @testset "Threshold behavior" begin + g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + g2 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.3, 0.3]) + + basis = BasisSet([g1]) + overlap = normalized_overlap(g2, g1) + + @test is_linearly_independent(g2, basis; threshold=overlap + 0.01) + end + + @testset "Invalid threshold" begin + g = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + basis = BasisSet([g]) + + @test_throws ArgumentError is_linearly_independent(g, basis; threshold=0.0) + @test_throws ArgumentError is_linearly_independent(g, basis; threshold=1.0) + @test_throws ArgumentError is_linearly_independent(g, basis; threshold=-0.5) + @test_throws ArgumentError is_linearly_independent(g, basis; threshold=1.5) + end end diff --git a/test/test_matrix_elements.jl b/test/test_matrix_elements.jl index 7667b91..ec5dd13 100644 --- a/test/test_matrix_elements.jl +++ b/test/test_matrix_elements.jl @@ -6,18 +6,6 @@ import FewBodyECG: _compute_matrix_element @testset "Overlap ⟨g′|g⟩" begin - @testset "Identical Gaussians" begin - A = [1.0 0.0; 0.0 1.0] - s = [0.0, 0.0] - g = Rank0Gaussian(A, s) - - overlap = _compute_matrix_element(g, g) - - expected = (π^2 / det(2 * A))^(3 / 2) - @test overlap ≈ expected rtol = 1.0e-10 - @test overlap > 0 - end - @testset "Symmetry ⟨g′|g⟩ = ⟨g|g′⟩" begin A1 = [1.0 0.0; 0.0 2.0] A2 = [1.5 0.0; 0.0 1.5] @@ -33,22 +21,6 @@ import FewBodyECG: _compute_matrix_element @test overlap_12 ≈ overlap_21 rtol = 1.0e-10 end - @testset "No shift vectors" begin - A1 = [1.0 0.0; 0.0 1.0] - A2 = [2.0 0.0; 0.0 2.0] - s = [0.0, 0.0] - - g1 = Rank0Gaussian(A1, s) - g2 = Rank0Gaussian(A2, s) - - overlap = _compute_matrix_element(g1, g2) - - @test overlap > 0 - @test isfinite(overlap) - - expected = (π^2 / det(A1 + A2))^(3 / 2) - @test overlap ≈ expected rtol = 1.0e-10 - end @testset "With shift vectors" begin A = [1.0 0.0; 0.0 1.0] @@ -66,29 +38,6 @@ import FewBodyECG: _compute_matrix_element @test isfinite(overlap) end - @testset "1D case" begin - A = [1.0;;] - s = [0.0] - g = Rank0Gaussian(A, s) - - overlap = _compute_matrix_element(g, g) - - expected = (π / 2.0)^(3 / 2) - @test overlap ≈ expected rtol = 1.0e-10 - end - - @testset "Numerical stability" begin - for scale in [0.1, 1.0, 10.0] - A = scale * [1.0 0.0; 0.0 1.0] - s = [0.0, 0.0] - g = Rank0Gaussian(A, s) - - overlap = _compute_matrix_element(g, g) - @test isfinite(overlap) - @test overlap > 0 - @test !isnan(overlap) - end - end end @testset "Kinetic Energy ⟨g′|K|g⟩" begin @@ -174,7 +123,6 @@ end result = _compute_matrix_element(g, g, V) - @test result < 0 @test isfinite(result) @test !isnan(result) end @@ -189,7 +137,6 @@ end result = _compute_matrix_element(g, g, V) - @test result > 0 @test isfinite(result) end @@ -259,18 +206,3 @@ end @test result1 ≈ result2 rtol = 1.0e-10 end end - -@testset "Error Cases" begin - - @testset "Dimension mismatch" begin - A1 = [1.0;;] # 1D - A2 = [1.0 0.0; 0.0 1.0] # 2D - s1 = [0.0] - s2 = [0.0, 0.0] - - g1 = Rank0Gaussian(A1, s1) - g2 = Rank0Gaussian(A2, s2) - - @test_throws Exception _compute_matrix_element(g1, g2) - end -end diff --git a/test/test_sampling.jl b/test/test_sampling.jl index ba61021..3d225b8 100644 --- a/test/test_sampling.jl +++ b/test/test_sampling.jl @@ -1,7 +1,50 @@ using Test using FewBodyECG +using LinearAlgebra @testset "Sampling Module Tests" begin + @testset "generate bij" begin + b = generate_bij(:quasirandom, 1, 5, 2.0; bmin = 0.5, bmax = 4.0) + @test length(b) == 5 + @test all(x -> x >= 0.5 && x <= 4.0, b) + end + @testset "Similar" begin + + bq1 = generate_bij(:quasirandom, 2, 4, 1.0) + bq2 = generate_bij(:quasirandom, 2, 4, 1.0) + @test bq1 ≈ bq2 + end + + @testset "Errors" begin + @test_throws ErrorException generate_bij(:unsupported, 1, 3, 1.0) + end + + @testset "Full" begin + + s = generate_shift(:quasirandom, 1, 3, 2.0) + @test length(s) == 3 + @test all(abs.(s) .<= 2.0 .+ eps()) + + s_q1 = generate_shift(:quasirandom, 5, 3, 1.5) + s_q2 = generate_shift(:quasirandom, 5, 3, 1.5) + @test s_q1 ≈ s_q2 + + @test_throws ErrorException generate_shift(:nope, 1, 2, 1.0) + + bij = [1.0, 2.0] + w1 = [1.0, 0.0, 0.0] + w2 = [0.0, 1.0, 0.0] + A = _generate_A_matrix(bij, [w1, w2]) + expected = Matrix(Diagonal([1.0, 1.0 / 4.0, 0.0])) + @test A ≈ expected + + @test_throws ArgumentError _generate_A_matrix([1.0, 2.0], [w1]) # length mismatch + @test_throws ArgumentError _generate_A_matrix([1.0, 2.0], [w1, [1.0, 0.0]]) # differing dimensions + + svec = [0.1, 0.2, 0.3] + rg = build_rank0([1.0, 2.0], [w1, w2], svec) + @test isa(rg, FewBodyECG.Rank0Gaussian) || isa(rg, Rank0Gaussian) + end end From 21d98db1158581ce5761684434456a0ac8df7e4d Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Wed, 19 Nov 2025 13:56:01 +0100 Subject: [PATCH 10/15] updated tests --- test/runtests.jl | 3 +- test/test_hamiltonian.jl | 122 +++++++++++++++++++-------------------- test/test_sampling.jl | 2 +- test/test_types.jl | 39 +++++++++++++ 4 files changed, 103 insertions(+), 63 deletions(-) create mode 100644 test/test_types.jl diff --git a/test/runtests.jl b/test/runtests.jl index 3bea91f..8d4efc7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,5 +9,6 @@ using FewBodyECG include("test_coordinates.jl") include("test_matrix_elements.jl") include("test_hamiltonian.jl") - + include("test_types.jl") + end diff --git a/test/test_hamiltonian.jl b/test/test_hamiltonian.jl index 0f6a44f..c3add76 100644 --- a/test/test_hamiltonian.jl +++ b/test/test_hamiltonian.jl @@ -6,188 +6,188 @@ import FewBodyECG: _compute_overlap_element, _build_operator_matrix, _compute_ma @testset "build_overlap_matrix" begin - + @testset "Size and symmetry" begin g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) g2 = Rank0Gaussian([1.5 0.0; 0.0 1.5], [0.0, 0.0]) g3 = Rank0Gaussian([2.0 0.0; 0.0 2.0], [0.1, 0.1]) - + basis = BasisSet([g1, g2, g3]) S = build_overlap_matrix(basis) - + @test size(S) == (3, 3) @test issymmetric(S) @test all(isfinite, S) end - - + + @testset "Consistency" begin g1 = Rank0Gaussian([1.0;;], [0.0]) g2 = Rank0Gaussian([2.0;;], [0.0]) - + basis = BasisSet([g1, g2]) S = build_overlap_matrix(basis) - + S11 = _compute_matrix_element(g1, g1) S12 = _compute_matrix_element(g1, g2) S22 = _compute_matrix_element(g2, g2) - - @test S[1,1] ≈ S11 rtol=1e-12 - @test S[1,2] ≈ S12 rtol=1e-12 - @test S[2,2] ≈ S22 rtol=1e-12 + + @test S[1, 1] ≈ S11 rtol = 1.0e-12 + @test S[1, 2] ≈ S12 rtol = 1.0e-12 + @test S[2, 2] ≈ S22 rtol = 1.0e-12 end end @testset "build_hamiltonian_matrix" begin - + @testset "Size and symmetry" begin g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) g2 = Rank0Gaussian([2.0 0.0; 0.0 2.0], [0.0, 0.0]) basis = BasisSet([g1, g2]) - + K = KineticOperator([0.5 0.0; 0.0 0.5]) H = build_hamiltonian_matrix(basis, [K]) - + @test size(H) == (2, 2) @test issymmetric(H) @test all(isfinite, H) end - + @testset "Additivity" begin # H(K+V) should equal H(K) + H(V) g = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) basis = BasisSet([g]) - + K = KineticOperator([0.5 0.0; 0.0 0.5]) V = CoulombOperator(-1.0, [1.0, 0.0]) - + H_K = build_hamiltonian_matrix(basis, [K]) H_V = build_hamiltonian_matrix(basis, [V]) H_both = build_hamiltonian_matrix(basis, [K, V]) - - @test H_both ≈ H_K + H_V rtol=1e-12 + + @test H_both ≈ H_K + H_V rtol = 1.0e-12 end - + @testset "Empty operators" begin g = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) basis = BasisSet([g]) - + H = build_hamiltonian_matrix(basis, FewBodyHamiltonians.Operator[]) - + @test H == zeros(1, 1) end end @testset "solve_generalized_eigenproblem" begin - + @testset "Basic solve" begin H = [2.0 0.5; 0.5 3.0] S = [1.0 0.1; 0.1 1.0] - + evals, evecs = solve_generalized_eigenproblem(H, S) - + @test length(evals) == 2 @test size(evecs) == (2, 2) @test all(isfinite, evals) @test all(isfinite, evecs) end - + @testset "Eigenvalue equation H*v = λ*S*v" begin H = [3.0 0.5; 0.5 2.0] S = [1.0 0.2; 0.2 1.0] - + evals, evecs = solve_generalized_eigenproblem(H, S) - + # Check each eigenpair for i in 1:2 λ = evals[i] v = evecs[:, i] residual = H * v - λ * S * v - @test norm(residual) < 1e-9 + @test norm(residual) < 1.0e-9 end end - + @testset "Requires positive definite S" begin H = [1.0 0.0; 0.0 1.0] S_bad = [-1.0 0.0; 0.0 1.0] - + @test_throws Exception solve_generalized_eigenproblem(H, S_bad) end end @testset "normalized_overlap" begin - + @testset "Symmetry" begin g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) g2 = Rank0Gaussian([2.0 0.0; 0.0 2.0], [0.5, 0.5]) - + overlap_12 = normalized_overlap(g1, g2) overlap_21 = normalized_overlap(g2, g1) - - @test overlap_12 ≈ overlap_21 rtol=1e-12 + + @test overlap_12 ≈ overlap_21 rtol = 1.0e-12 end - + @testset "Range [0, 1]" begin g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) g2 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.5, 0.5]) g3 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [5.0, 5.0]) - + overlap_12 = normalized_overlap(g1, g2) overlap_13 = normalized_overlap(g1, g3) - - @test 0.0 <= overlap_12 <= 1.0 - @test 0.0 <= overlap_13 <= 1.0 + + @test 0.0 <= overlap_12 <= 1.0 + @test 0.0 <= overlap_13 <= 1.0 end - + @testset "Identical Gaussians" begin A = [1.5 0.0; 0.0 1.5] s = [0.3, 0.4] g1 = Rank0Gaussian(A, s) g2 = Rank0Gaussian(A, s) - + overlap = normalized_overlap(g1, g2) - @test overlap ≈ 1.0 rtol=1e-10 + @test overlap ≈ 1.0 rtol = 1.0e-10 end end @testset "is_linearly_independent" begin - + @testset "Empty basis (always independent)" begin g = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) empty_basis = BasisSet(Rank0Gaussian[]) - - @test is_linearly_independent(g, empty_basis; threshold=0.95) + + @test is_linearly_independent(g, empty_basis; threshold = 0.95) end - - + + @testset "Well-separated (should accept)" begin g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) g2 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [5.0, 5.0]) - + basis = BasisSet([g1]) - - @test is_linearly_independent(g2, basis; threshold=0.95) + + @test is_linearly_independent(g2, basis; threshold = 0.95) end - + @testset "Threshold behavior" begin g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) g2 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.3, 0.3]) - + basis = BasisSet([g1]) overlap = normalized_overlap(g2, g1) - - @test is_linearly_independent(g2, basis; threshold=overlap + 0.01) + + @test is_linearly_independent(g2, basis; threshold = overlap + 0.01) end - + @testset "Invalid threshold" begin g = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) basis = BasisSet([g]) - - @test_throws ArgumentError is_linearly_independent(g, basis; threshold=0.0) - @test_throws ArgumentError is_linearly_independent(g, basis; threshold=1.0) - @test_throws ArgumentError is_linearly_independent(g, basis; threshold=-0.5) - @test_throws ArgumentError is_linearly_independent(g, basis; threshold=1.5) + + @test_throws ArgumentError is_linearly_independent(g, basis; threshold = 0.0) + @test_throws ArgumentError is_linearly_independent(g, basis; threshold = 1.0) + @test_throws ArgumentError is_linearly_independent(g, basis; threshold = -0.5) + @test_throws ArgumentError is_linearly_independent(g, basis; threshold = 1.5) end end diff --git a/test/test_sampling.jl b/test/test_sampling.jl index 3d225b8..de48972 100644 --- a/test/test_sampling.jl +++ b/test/test_sampling.jl @@ -46,5 +46,5 @@ using LinearAlgebra svec = [0.1, 0.2, 0.3] rg = build_rank0([1.0, 2.0], [w1, w2], svec) @test isa(rg, FewBodyECG.Rank0Gaussian) || isa(rg, Rank0Gaussian) - end + end end diff --git a/test/test_types.jl b/test/test_types.jl new file mode 100644 index 0000000..69cbae7 --- /dev/null +++ b/test/test_types.jl @@ -0,0 +1,39 @@ +using Test +using LinearAlgebra +using FewBodyECG + +@testset "Rank0Gaussian constructor and validate!" begin + A = [2.0 0.0; 0.0 3.0] + s = [1.0, 2.0] + g = Rank0Gaussian(A, s) + + @test isa(g, Rank0Gaussian) + @test isa(g.A, Symmetric) + @test g.s == s + + A_ns = rand(2, 3) + @test_throws ArgumentError Rank0Gaussian(A_ns, [1.0, 2.0]) + + @test_throws ArgumentError Rank0Gaussian(A, [1.0]) + + A_indef = [0.0 -1.0; -1.0 0.0] + g_indef = Rank0Gaussian(A_indef, s) +end + +@testset "Rank1Gaussian constructor and validate!" begin + A = [3.0 0.0; 0.0 4.0] + a = [0.1, 0.2] + s = [1.0, 2.0] + g1 = Rank1Gaussian(A, a, s) + + @test isa(g1, Rank1Gaussian) + @test isa(g1.A, Symmetric) + @test g1.a == a + @test g1.s == s + + @test_throws ArgumentError Rank1Gaussian(rand(2, 3), a, s) + + @test_throws ArgumentError Rank1Gaussian(A, [1.0], s) + @test_throws ArgumentError Rank1Gaussian(A, a, [1.0]) + +end From 967d6000f8c6f60d7e01f6cdd5d9b4173b84ad60 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Wed, 19 Nov 2025 14:02:30 +0100 Subject: [PATCH 11/15] updated matrix elements and tests --- src/matrix_elements.jl | 100 +---------------------------------- test/test_types.jl | 115 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 116 insertions(+), 99 deletions(-) diff --git a/src/matrix_elements.jl b/src/matrix_elements.jl index d18ab46..ac55e26 100644 --- a/src/matrix_elements.jl +++ b/src/matrix_elements.jl @@ -71,102 +71,4 @@ function _compute_matrix_element(bra::Rank1Gaussian, ket::Rank1Gaussian, op::Kin T5 = dot(a, R * A * b) * M0 return T1 + T2 + T3 - T4 - T5 -end - -function _compute_matrix_element(bra::Rank2Gaussian, ket::Rank2Gaussian, op::KineticOperator) - A, B = bra.A, ket.A - a, b, c, d = bra.a, bra.b, ket.a, ket.b - K = op.K - R = inv(A + B) - M0 = (π^length(R) / det(A + B))^(3 / 2) - - M2 = 0.25 * ( - dot(a, R * b) * dot(c, R * d) + - dot(a, R * c) * dot(b, R * d) + - dot(a, R * d) * dot(b, R * c) - ) * M0 - - T1 = 6 * tr(B * K * A * R) * M2 - - T2 = 0.5 * ( - dot(a, K * c) * dot(b, R * d) + - dot(a, K * d) * dot(b, R * c) + - dot(b, K * c) * dot(a, R * d) + - dot(b, K * d) * dot(a, R * c) - ) * M0 - - T3 = 0.5 * ( - dot(a, R * B * K * A * R * b) * dot(c, R * d) + - dot(a, R * B * K * A * R * c) * dot(b, R * d) + - dot(a, R * B * K * A * R * d) * dot(b, R * c) + - dot(b, R * B * K * A * R * a) * dot(c, R * d) + - dot(b, R * B * K * A * R * c) * dot(a, R * d) + - dot(b, R * B * K * A * R * d) * dot(a, R * c) + - dot(c, R * B * K * A * R * a) * dot(b, R * d) + - dot(c, R * B * K * A * R * b) * dot(a, R * d) + - dot(c, R * B * K * A * R * d) * dot(a, R * b) + - dot(d, R * B * K * A * R * a) * dot(b, R * c) + - dot(d, R * B * K * A * R * b) * dot(a, R * c) + - dot(d, R * B * K * A * R * c) * dot(a, R * b) - ) * M0 - - T4 = -0.5 * ( - dot(a, R * B * K * b) * dot(c, R * d) + - dot(b, R * B * K * a) * dot(c, R * d) + - dot(c, R * B * K * a) * dot(b, R * d) + - dot(c, R * B * K * b) * dot(a, R * d) + - dot(d, R * B * K * a) * dot(b, R * c) + - dot(d, R * B * K * b) * dot(a, R * c) - ) * M0 - - T5 = -0.5 * ( - dot(c, K * A * a) * dot(b, R * d) + - dot(c, K * A * b) * dot(a, R * d) + - dot(c, K * A * d) * dot(a, R * b) + - dot(d, K * A * a) * dot(b, R * c) + - dot(d, K * A * b) * dot(a, R * c) + - dot(d, K * A * c) * dot(a, R * b) - ) * M0 - - return T1 + T2 + T3 + T4 + T5 -end - -function _compute_matrix_element(bra::Rank2Gaussian, ket::Rank2Gaussian, op::CoulombOperator) - A, B = bra.A, ket.A - a, b, c, d = bra.a, bra.b, ket.a, ket.b - w = op.w - R = inv(A + B) - M0 = (π^length(R) / det(A + B))^(3 / 2) - - β = 1 / dot(w, R * w) - q = 0.5 * dot(w, R * (a + b + c + d)) - - M2 = 0.25 * ( - dot(a, R * b) * dot(c, R * d) + - dot(a, R * c) * dot(b, R * d) + - dot(a, R * d) * dot(b, R * c) - ) * M0 - - term1 = 2 * sqrt(β / π) * M2 - - q2_1 = dot(a, R * (w * w') * b) * dot(c, R * d) - q2_2 = dot(a, R * (w * w') * c) * dot(b, R * d) - q2_3 = dot(a, R * (w * w') * d) * dot(b, R * c) - q2_4 = dot(b, R * (w * w') * c) * dot(a, R * d) - q2_5 = dot(b, R * (w * w') * d) * dot(a, R * c) - q2_6 = dot(c, R * (w * w') * d) * dot(a, R * b) - - q4_1 = dot(a, R * (w * w') * b) * dot(c, R * (w * w') * d) - q4_2 = dot(a, R * (w * w') * c) * dot(b, R * (w * w') * d) - q4_3 = dot(a, R * (w * w') * d) * dot(b, R * (w * w') * c) - - term2 = -2 * sqrt(β / π) * β / 3 * 0.25 * ( - q2_1 + q2_2 + q2_3 + q2_4 + q2_5 + q2_6 - ) * M0 - - term3 = 2 * sqrt(β / π) * β^2 / 10 * 0.5 * ( - q4_1 + q4_2 + q4_3 - ) * M0 - - return op.coefficient * (term1 + term2 + term3) -end +end \ No newline at end of file diff --git a/test/test_types.jl b/test/test_types.jl index 69cbae7..7db3d60 100644 --- a/test/test_types.jl +++ b/test/test_types.jl @@ -1,6 +1,8 @@ using Test using LinearAlgebra using FewBodyECG +import FewBodyECG: validate! +using FewBodyHamiltonians @testset "Rank0Gaussian constructor and validate!" begin A = [2.0 0.0; 0.0 3.0] @@ -37,3 +39,116 @@ end @test_throws ArgumentError Rank1Gaussian(A, a, [1.0]) end + + +@testset "Rank0Gaussian constructor and validate!" begin + A = [2.0 0.0; 0.0 3.0] + s = [1.0, 2.0] + g = Rank0Gaussian(A, s) + + @test isa(g, Rank0Gaussian) + @test isa(g.A, Symmetric) + @test g.s == s + + A_ns = rand(2, 3) + @test_throws ArgumentError Rank0Gaussian(A_ns, [1.0, 2.0]) + + @test_throws ArgumentError Rank0Gaussian(A, [1.0]) + + A_indef = [0.0 -1.0; -1.0 0.0] + g_indef = Rank0Gaussian(A_indef, s) +end + +@testset "Rank1Gaussian constructor and validate!" begin + A = [3.0 0.0; 0.0 4.0] + a = [0.1, 0.2] + s = [1.0, 2.0] + g1 = Rank1Gaussian(A, a, s) + + @test isa(g1, Rank1Gaussian) + @test isa(g1.A, Symmetric) + @test g1.a == a + @test g1.s == s + + @test_throws ArgumentError Rank1Gaussian(rand(2, 3), a, s) + + @test_throws ArgumentError Rank1Gaussian(A, [1.0], s) + @test_throws ArgumentError Rank1Gaussian(A, a, [1.0]) + +end + +@testset "Rank2Gaussian constructor and validate!" begin + A = [4.0 0.0; 0.0 5.0] + a = [0.1, 0.2] + b = [0.3, 0.4] + s = [1.0, 2.0] + + g2 = Rank2Gaussian(A, a, b, s) + @test isa(g2, Rank2Gaussian) + @test isa(g2.A, Symmetric) + @test g2.a == a + @test g2.b == b + @test g2.s == s + + @test_throws ArgumentError Rank2Gaussian(rand(2,3), a, b, s) + @test_throws ArgumentError Rank2Gaussian(A, [1.0], b, s) + @test_throws ArgumentError Rank2Gaussian(A, a, [1.0], s) + @test_throws ArgumentError Rank2Gaussian(A, a, b, [1.0]) + + # validate! succeeds for positive-definite A + @test validate!(g2) === g2 + + # validate! throws for indefinite A + A_indef = [0.0 -1.0; -1.0 0.0] + g2_indef = Rank2Gaussian(A_indef, a, b, s) + @test_throws LinearAlgebra.PosDefException validate!(g2_indef) +end + +@testset "validate! for Rank0 and Rank1 positive/negative-definite" begin + A_pd = [2.0 0.0; 0.0 2.0] + s = [0.0, 0.0] + g0 = Rank0Gaussian(A_pd, s) + @test validate!(g0) === g0 + + A_indef = [0.0 -1.0; -1.0 0.0] + g0_indef = Rank0Gaussian(A_indef, s) + @test_throws LinearAlgebra.PosDefException validate!(g0_indef) + + A1 = [1.0 0.0; 0.0 1.0] + a = [0.0, 0.0] + g1 = Rank1Gaussian(A1, a, s) + @test validate!(g1) === g1 + + g1_indef = Rank1Gaussian(A_indef, a, s) + @test_throws LinearAlgebra.PosDefException validate!(g1_indef) +end + +@testset "BasisSet, KineticOperator, CoulombOperator, and ECG composition" begin + A = [2.0 0.0; 0.0 3.0] + s1 = [1.0, 0.0] + s2 = [0.0, 1.0] + g1 = Rank0Gaussian(A, s1) + g2 = Rank0Gaussian(A, s2) + + bset = BasisSet([g1, g2]) + @test isa(bset, BasisSet) + @test length(bset.functions) == 2 + @test bset.functions[1] == g1 + @test bset.functions[2] == g2 + + K = [1.0 0.0; 0.0 1.0] + kop = KineticOperator(K) + @test kop.K == K + @test kop isa FewBodyHamiltonians.KineticTerm + + coeff = 2.5 + w = [1.0, -1.0] + cop = CoulombOperator(coeff, w) + @test cop.coefficient == coeff + @test cop.w == w + @test cop isa FewBodyHamiltonians.PotentialTerm + + ecg = ECG(bset, [kop, cop]) + @test ecg.basis === bset + @test ecg.operators == [kop, cop] +end From 4ccd9acb5cc645fdace2c24b7724f7eda15409a1 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Wed, 19 Nov 2025 14:10:03 +0100 Subject: [PATCH 12/15] Removed rank1 matrix elements --- src/matrix_elements.jl | 29 ----------------------------- test/runtests.jl | 2 +- test/test_types.jl | 2 +- 3 files changed, 2 insertions(+), 31 deletions(-) diff --git a/src/matrix_elements.jl b/src/matrix_elements.jl index ac55e26..b78b811 100644 --- a/src/matrix_elements.jl +++ b/src/matrix_elements.jl @@ -43,32 +43,3 @@ function _compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian, op::Cou f = abs(q) < 1.0e-12 ? (2 * sqrt(β / π)) : (erf(sqrt(β) * q) / q) return op.coefficient * f * M end - -function _compute_matrix_element(bra::Rank1Gaussian, ket::Rank1Gaussian, op::CoulombOperator) - A, B, a, b, w = bra.A, ket.A, bra.a, ket.a, op.w - R = inv(A + B) - β = 1 / (dot(w, R * w)) - M0 = (π^length(R) / det(A + B))^(3 / 2) - M1 = 0.5 * dot(b, R * a) * M0 - q2 = 0.25 * dot(a .+ b, R * (w * w') * (a .+ b)) - return 2 * sqrt(β / π) * M1 - sqrt(β^3 / π) / 3 * q2 * M0 -end - -function _compute_matrix_element(bra::Rank1Gaussian, ket::Rank1Gaussian, op::KineticOperator) - a = bra.a isa AbstractVector{<:AbstractVector} ? vec(bra.a[1]) : vec(bra.a) - b = ket.a isa AbstractVector{<:AbstractVector} ? vec(ket.a[1]) : vec(ket.a) - - A, B = bra.A, ket.A - K = op.K - R = inv(A + B) - M0 = (π^length(R) / det(A + B))^(3 / 2) - M1 = 0.5 * dot(b, R * a) * M0 - - T1 = 6 * tr(B * K * A * R) * M1 - T2 = dot(b, a) * M0 - T3 = dot(a, R * B * A * R * b) * M0 - T4 = dot(b, R * B * a) * M0 - T5 = dot(a, R * A * b) * M0 - - return T1 + T2 + T3 - T4 - T5 -end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 8d4efc7..d43d85d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,5 +10,5 @@ using FewBodyECG include("test_matrix_elements.jl") include("test_hamiltonian.jl") include("test_types.jl") - + end diff --git a/test/test_types.jl b/test/test_types.jl index 7db3d60..2190e2d 100644 --- a/test/test_types.jl +++ b/test/test_types.jl @@ -90,7 +90,7 @@ end @test g2.b == b @test g2.s == s - @test_throws ArgumentError Rank2Gaussian(rand(2,3), a, b, s) + @test_throws ArgumentError Rank2Gaussian(rand(2, 3), a, b, s) @test_throws ArgumentError Rank2Gaussian(A, [1.0], b, s) @test_throws ArgumentError Rank2Gaussian(A, a, [1.0], s) @test_throws ArgumentError Rank2Gaussian(A, a, b, [1.0]) From 872539fabead966980325bb49ea3bf36f8cb6660 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Thu, 11 Dec 2025 12:40:54 +0100 Subject: [PATCH 13/15] updated sampling --- Examples/Positronium.jl | 8 ++- "Examples/td\316\274.jl" | 12 +++- src/FewBodyECG.jl | 2 +- src/hamiltonian.jl | 135 +++++++++++++++++++++++++++++++++++---- src/sampling.jl | 5 +- test/test_hamiltonian.jl | 4 +- 6 files changed, 145 insertions(+), 21 deletions(-) diff --git a/Examples/Positronium.jl b/Examples/Positronium.jl index 1606b60..145afcd 100644 --- a/Examples/Positronium.jl +++ b/Examples/Positronium.jl @@ -1,4 +1,7 @@ using FewBodyECG, LinearAlgebra +using QuasiMonteCarlo +using Plots +import FewBodyECG: default_scale, convergence masses = [1.0, 1.0, 1.0] @@ -9,10 +12,11 @@ J, U = _jacobi_transform(masses) w_list = [[1, -1, 0], [1, 0, -1], [0, 1, -1]] w_raw = [U' * w for w in w_list] -coeffs = [+1.0, -1.0, -1.0] +coeffs = [-1.0, -1.0, +1.0] coulomb_ops = [CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw)] ops = Operator[kin; coulomb_ops...] +scale = default_scale(masses) -result = solve_ECG(ops, 250, sampler = SobolSample(); scale = 0.2) +result = solve_ECG(ops, 300, sampler = SobolSample(); scale = scale) println("E ≈ ", result.ground_state) diff --git "a/Examples/td\316\274.jl" "b/Examples/td\316\274.jl" index 3d9a4fe..c388d98 100644 --- "a/Examples/td\316\274.jl" +++ "b/Examples/td\316\274.jl" @@ -1,4 +1,7 @@ using FewBodyECG, LinearAlgebra +using QuasiMonteCarlo +import FewBodyECG: default_scale, convergence +using Plots masses = [5496.918, 3670.481, 206.7686] @@ -7,13 +10,16 @@ kin = KineticOperator(Λmat) J, U = _jacobi_transform(masses) base_w = [[1, -1, 0], [1, 0, -1], [0, 1, -1]] -coeffs = [+1.0, +1.0, -1.0] +coeffs = [+1.0, -1.0, -1.0] w_list = [c .* w for (c, w) in zip(coeffs, base_w)] w_raw = [U' * w for w in w_list] coulomb_ops = [CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw)] ops = Operator[kin; coulomb_ops...] - -result = solve_ECG(ops, 250, sampler = HaltonSample(); scale = 0.025) +scale = default_scale(masses) +result = solve_ECG(ops, 550, sampler = HaltonSample(); scale = scale) println("E ≈ ", result.ground_state) + +a,b = correlation_function(result) +conv = FewBodyECG.convergence(result) \ No newline at end of file diff --git a/src/FewBodyECG.jl b/src/FewBodyECG.jl index c596b34..90e84e2 100644 --- a/src/FewBodyECG.jl +++ b/src/FewBodyECG.jl @@ -15,7 +15,7 @@ export Operator export build_hamiltonian_matrix, build_overlap_matrix, solve_generalized_eigenproblem, solve_ECG, convergence -export ψ₀, SolverResults, convergence +export ψ₀, SolverResults, convergence, correlation_function, ψ include("types.jl") include("coordinates.jl") diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl index 7cf8e22..44e0a34 100644 --- a/src/hamiltonian.jl +++ b/src/hamiltonian.jl @@ -36,13 +36,84 @@ function build_hamiltonian_matrix(basis::BasisSet{<:GaussianBase}, operators::Ab return H end -function solve_generalized_eigenproblem(H::AbstractMatrix{<:Real}, S::AbstractMatrix{<:Real}) - F = cholesky(Symmetric(S); check = true) +function solve_generalized_eigenproblem( + H::AbstractMatrix{<:Real}, + S::AbstractMatrix{<:Real}; + max_condition::Real = 1.0e12, + regularization::Real = 0.0 + ) + + if any(!isfinite, H) + error("Hamiltonian matrix H contains NaN or Inf values") + end + if any(!isfinite, S) + error("Overlap matrix S contains NaN or Inf values") + end + + H_sym = Symmetric((H + H') / 2) + S_sym = Symmetric((S + S') / 2) + + cond_S = cond(S_sym) + if cond_S > max_condition + @warn "Overlap matrix poorly conditioned (κ=$cond_S), adding regularization" + if regularization == 0.0 + regularization = maximum(abs.(diag(S_sym))) * 1.0e-10 + end + end + + if regularization > 0 + S_sym = S_sym + regularization * I + end + + if !isposdef(S_sym) + @warn "Overlap matrix not positive definite, adding regularization" + ε = maximum(abs.(diag(S_sym))) * 1.0e-8 + S_sym = S_sym + ε * I + + if !isposdef(S_sym) + error("Overlap matrix not positive definite even after regularization") + end + end + + # Cholesky decomposition with error handling + local F + try + F = cholesky(S_sym) + catch e + @error "Cholesky decomposition failed" exception = e + @error "Overlap matrix info" condition = cond(S_sym) min_eigval = minimum(eigvals(S_sym)) + rethrow(e) + end + L = F.L - A = (L \ H) / L' - evals, evecs = eigen(Symmetric(A)) + + # Transform to standard eigenvalue problem + A = (L \ Matrix(H_sym)) / L' + A_sym = Symmetric((A + A') / 2) + + # Check for NaN/Inf after transformation + if any(!isfinite, A_sym) + error("Transformed matrix contains NaN or Inf after Cholesky transformation") + end + + # Solve standard eigenvalue problem + local evals, evecs + try + evals, evecs = eigen(A_sym) + catch e + @error "Eigenvalue decomposition failed" exception = e + @error "Transformed matrix info" condition = cond(A_sym) + rethrow(e) + end + + # Transform eigenvectors back vecs = L' \ evecs - return real(evals), real(vecs) + + # Ensure real + evals_real = real.(evals) + vecs_real = real.(vecs) + + return evals_real, vecs_real end function normalized_overlap(A::GaussianBase, B::GaussianBase) @@ -78,6 +149,11 @@ function is_linearly_independent( return true end +function default_scale(masses::Vector{<:Real}) + μ = minimum(masses[masses .< 1e10]) + return 1 / sqrt(μ) +end + function solve_ECG( operators::Vector{<:FewBodyHamiltonians.Operator}, n::Int = 50; @@ -86,6 +162,7 @@ function solve_ECG( scale::Real = 0.2, threshold::Real = 0.95, max_attempts::Int = 10 * n, + max_condition::Real = 1e12, verbose::Bool = true ) @@ -110,6 +187,7 @@ function solve_ECG( s = generate_shift(method, attempt, d, scale; qmc_sampler = sampler) candidate = Rank0Gaussian(A, s) + # Check linear independence if !isempty(basis_fns) existing_basis = BasisSet{Rank0Gaussian}(basis_fns) if !is_linearly_independent(candidate, existing_basis; threshold = threshold) @@ -120,13 +198,48 @@ function solve_ECG( end push!(basis_fns, candidate) - n_accepted += 1 - basis = BasisSet{Rank0Gaussian}(basis_fns) - H = build_hamiltonian_matrix(basis, operators) - S = build_overlap_matrix(basis) + local H, S, λs, Us + try + basis = BasisSet{Rank0Gaussian}(basis_fns) + H = build_hamiltonian_matrix(basis, operators) + S = build_overlap_matrix(basis) + + # Check for NaN/Inf BEFORE eigensolve + if any(!isfinite, H) + @warn "Hamiltonian contains NaN/Inf at step $(n_accepted+1), rejecting basis function" + pop!(basis_fns) + n_rejected += 1 + continue + end + + if any(!isfinite, S) + @warn "Overlap contains NaN/Inf at step $(n_accepted+1), rejecting basis function" + pop!(basis_fns) + n_rejected += 1 + continue + end + + # Check condition number + cond_S = cond(Symmetric(S)) + if cond_S > max_condition + @warn "Overlap poorly conditioned (κ=$cond_S) at step $(n_accepted+1), rejecting" + pop!(basis_fns) + n_rejected += 1 + continue + end - λs, Us = solve_generalized_eigenproblem(H, S) + # Try to solve + λs, Us = solve_generalized_eigenproblem(H, S; max_condition) + + catch e + @warn "Failed at step $(n_accepted+1): $e" + pop!(basis_fns) # Remove problematic function + n_rejected += 1 + continue + end + + n_accepted += 1 E0 = minimum(λs) push!(E_hist, E0) @@ -141,4 +254,4 @@ function solve_ECG( Emin = last(E_hist) @info "Optimization complete" E₀ = Emin n_basis = n_accepted return SolverResults(basis_fns, n_accepted, operators, method, sampler, b₁, Emin, E_hist, vecs_list) -end +end \ No newline at end of file diff --git a/src/sampling.jl b/src/sampling.jl index d6dd26d..248f421 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -6,13 +6,14 @@ function _qmc_point(i::Int, d::Int; sampler = HaltonSample()) return QuasiMonteCarlo.sample(i + 1, d, sampler)[:, end] end -function generate_bij(method::Symbol, i::Int, n_terms::Int, b1::Float64; qmc_sampler = HaltonSample(), bmin = 0.02 * b1, bmax = 20b1) +function generate_bij(method::Symbol, i::Int, n_terms::Int, b1::Float64; + qmc_sampler = HaltonSample(), bmin = 0.02 * b1, bmax = 5.0 * b1) bmin > 0 || throw(ArgumentError("bmin must be > 0")) bmax > bmin || throw(ArgumentError("bmax must be > bmin")) u = method === :quasirandom ? _qmc_point(i, n_terms; sampler = qmc_sampler) : method === :random ? rand(n_terms) : error("Unsupported method $method") - return bmin .* (bmax / bmin) .^ u + return bmin .+ (bmax - bmin) .* u # Now properly uses the range! end function generate_shift(method::Symbol, i::Int, dim::Int, scale::Real; qmc_sampler = HaltonSample()) diff --git a/test/test_hamiltonian.jl b/test/test_hamiltonian.jl index c3add76..aab1487 100644 --- a/test/test_hamiltonian.jl +++ b/test/test_hamiltonian.jl @@ -3,7 +3,7 @@ using LinearAlgebra using FewBodyHamiltonians using FewBodyECG import FewBodyECG: _compute_overlap_element, _build_operator_matrix, _compute_matrix_element, normalized_overlap, is_linearly_independent - +using QuasiMonteCarlo @testset "build_overlap_matrix" begin @@ -190,4 +190,4 @@ end @test_throws ArgumentError is_linearly_independent(g, basis; threshold = -0.5) @test_throws ArgumentError is_linearly_independent(g, basis; threshold = 1.5) end -end +end \ No newline at end of file From da766b0a84d7988f288fda7e2668ce51d046bbd5 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Thu, 11 Dec 2025 12:51:45 +0100 Subject: [PATCH 14/15] updated tests --- test/test_hamiltonian.jl | 423 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 422 insertions(+), 1 deletion(-) diff --git a/test/test_hamiltonian.jl b/test/test_hamiltonian.jl index aab1487..5c9dff5 100644 --- a/test/test_hamiltonian.jl +++ b/test/test_hamiltonian.jl @@ -190,4 +190,425 @@ end @test_throws ArgumentError is_linearly_independent(g, basis; threshold = -0.5) @test_throws ArgumentError is_linearly_independent(g, basis; threshold = 1.5) end -end \ No newline at end of file +end + +using Test +using LinearAlgebra +using FewBodyHamiltonians +using FewBodyECG +import FewBodyECG: _compute_overlap_element, _build_operator_matrix, _compute_matrix_element +import FewBodyECG: normalized_overlap, is_linearly_independent, default_scale +import FewBodyECG: _jacobi_transform, _generate_A_matrix, generate_bij, generate_shift +using QuasiMonteCarlo + +# ============================================================================= +# Existing tests (preserved) +# ============================================================================= + +@testset "build_overlap_matrix" begin + + @testset "Size and symmetry" begin + g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + g2 = Rank0Gaussian([1.5 0.0; 0.0 1.5], [0.0, 0.0]) + g3 = Rank0Gaussian([2.0 0.0; 0.0 2.0], [0.1, 0.1]) + + basis = BasisSet([g1, g2, g3]) + S = build_overlap_matrix(basis) + + @test size(S) == (3, 3) + @test issymmetric(S) + @test all(isfinite, S) + end + + + @testset "Consistency" begin + g1 = Rank0Gaussian([1.0;;], [0.0]) + g2 = Rank0Gaussian([2.0;;], [0.0]) + + basis = BasisSet([g1, g2]) + S = build_overlap_matrix(basis) + + S11 = _compute_matrix_element(g1, g1) + S12 = _compute_matrix_element(g1, g2) + S22 = _compute_matrix_element(g2, g2) + + @test S[1, 1] ≈ S11 rtol = 1.0e-12 + @test S[1, 2] ≈ S12 rtol = 1.0e-12 + @test S[2, 2] ≈ S22 rtol = 1.0e-12 + end +end + +@testset "build_hamiltonian_matrix" begin + + @testset "Size and symmetry" begin + g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + g2 = Rank0Gaussian([2.0 0.0; 0.0 2.0], [0.0, 0.0]) + basis = BasisSet([g1, g2]) + + K = KineticOperator([0.5 0.0; 0.0 0.5]) + H = build_hamiltonian_matrix(basis, [K]) + + @test size(H) == (2, 2) + @test issymmetric(H) + @test all(isfinite, H) + end + + @testset "Additivity" begin + # H(K+V) should equal H(K) + H(V) + g = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + basis = BasisSet([g]) + + K = KineticOperator([0.5 0.0; 0.0 0.5]) + V = CoulombOperator(-1.0, [1.0, 0.0]) + + H_K = build_hamiltonian_matrix(basis, [K]) + H_V = build_hamiltonian_matrix(basis, [V]) + H_both = build_hamiltonian_matrix(basis, [K, V]) + + @test H_both ≈ H_K + H_V rtol = 1.0e-12 + end + + @testset "Empty operators" begin + g = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + basis = BasisSet([g]) + + H = build_hamiltonian_matrix(basis, FewBodyHamiltonians.Operator[]) + + @test H == zeros(1, 1) + end +end + + +@testset "solve_generalized_eigenproblem" begin + + @testset "Basic solve" begin + H = [2.0 0.5; 0.5 3.0] + S = [1.0 0.1; 0.1 1.0] + + evals, evecs = solve_generalized_eigenproblem(H, S) + + @test length(evals) == 2 + @test size(evecs) == (2, 2) + @test all(isfinite, evals) + @test all(isfinite, evecs) + end + + @testset "Eigenvalue equation H*v = λ*S*v" begin + H = [3.0 0.5; 0.5 2.0] + S = [1.0 0.2; 0.2 1.0] + + evals, evecs = solve_generalized_eigenproblem(H, S) + + # Check each eigenpair + for i in 1:2 + λ = evals[i] + v = evecs[:, i] + residual = H * v - λ * S * v + @test norm(residual) < 1.0e-9 + end + end + + @testset "Requires positive definite S" begin + H = [1.0 0.0; 0.0 1.0] + S_bad = [-1.0 0.0; 0.0 1.0] + + @test_throws Exception solve_generalized_eigenproblem(H, S_bad) + end +end + + +@testset "normalized_overlap" begin + + @testset "Symmetry" begin + g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + g2 = Rank0Gaussian([2.0 0.0; 0.0 2.0], [0.5, 0.5]) + + overlap_12 = normalized_overlap(g1, g2) + overlap_21 = normalized_overlap(g2, g1) + + @test overlap_12 ≈ overlap_21 rtol = 1.0e-12 + end + + @testset "Range [0, 1]" begin + g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + g2 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.5, 0.5]) + g3 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [5.0, 5.0]) + + overlap_12 = normalized_overlap(g1, g2) + overlap_13 = normalized_overlap(g1, g3) + + @test 0.0 <= overlap_12 <= 1.0 + @test 0.0 <= overlap_13 <= 1.0 + end + + @testset "Identical Gaussians" begin + A = [1.5 0.0; 0.0 1.5] + s = [0.3, 0.4] + g1 = Rank0Gaussian(A, s) + g2 = Rank0Gaussian(A, s) + + overlap = normalized_overlap(g1, g2) + @test overlap ≈ 1.0 rtol = 1.0e-10 + end +end + +@testset "is_linearly_independent" begin + + @testset "Empty basis (always independent)" begin + g = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + empty_basis = BasisSet(Rank0Gaussian[]) + + @test is_linearly_independent(g, empty_basis; threshold = 0.95) + end + + + @testset "Well-separated (should accept)" begin + g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + g2 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [5.0, 5.0]) + + basis = BasisSet([g1]) + + @test is_linearly_independent(g2, basis; threshold = 0.95) + end + + @testset "Threshold behavior" begin + g1 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + g2 = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.3, 0.3]) + + basis = BasisSet([g1]) + overlap = normalized_overlap(g2, g1) + + @test is_linearly_independent(g2, basis; threshold = overlap + 0.01) + end + + @testset "Invalid threshold" begin + g = Rank0Gaussian([1.0 0.0; 0.0 1.0], [0.0, 0.0]) + basis = BasisSet([g]) + + @test_throws ArgumentError is_linearly_independent(g, basis; threshold = 0.0) + @test_throws ArgumentError is_linearly_independent(g, basis; threshold = 1.0) + @test_throws ArgumentError is_linearly_independent(g, basis; threshold = -0.5) + @test_throws ArgumentError is_linearly_independent(g, basis; threshold = 1.5) + end +end + +# ============================================================================= +# New tests for default_scale +# ============================================================================= + +@testset "default_scale" begin + + @testset "Basic functionality" begin + # Equal masses + masses_equal = [1.0, 1.0, 1.0] + scale = default_scale(masses_equal) + @test scale ≈ 1.0 rtol = 1.0e-10 + @test scale > 0 + end + + @testset "Muonic systems (heavy nuclei)" begin + # Muon mass ≈ 206.77 mₑ, so scale ~ 1/√207 ≈ 0.07 + masses_muonic = [5000.0, 5000.0, 206.77] + scale = default_scale(masses_muonic) + @test scale ≈ 1 / sqrt(206.77) rtol = 1.0e-6 + @test scale < 0.1 # Should be small for muonic systems + end + + @testset "Ignores infinite masses" begin + # Infinite mass (fixed nucleus) should be ignored + masses_fixed = [1.0e15, 1.0, 1.0] + scale = default_scale(masses_fixed) + @test scale ≈ 1.0 rtol = 1.0e-10 # Should use minimum of finite masses + end + + @testset "Single light particle" begin + masses = [1000.0, 1000.0, 1.0] + scale = default_scale(masses) + @test scale ≈ 1.0 rtol = 1.0e-10 + end + + @testset "Hydrogen-like (electron mass = 1)" begin + masses_hydrogen = [1836.15, 1.0] # proton, electron + scale = default_scale(masses_hydrogen) + @test scale ≈ 1.0 rtol = 1.0e-10 + end +end + +# ============================================================================= +# Tests for generate_bij +# ============================================================================= + +@testset "generate_bij" begin + + @testset "Output range with bmin/bmax" begin + b1 = 1.0 + bmin = 0.02 * b1 + bmax = 5.0 * b1 + + for i in 1:100 + bij = generate_bij(:quasirandom, i, 3, b1; qmc_sampler = HaltonSample()) + @test all(bij .>= bmin - 1e-10) + @test all(bij .<= bmax + 1e-10) + end + end + + @testset "No zero values (would cause singularity)" begin + b1 = 1.0 + for i in 1:100 + bij = generate_bij(:quasirandom, i, 3, b1; qmc_sampler = HaltonSample()) + @test all(bij .> 0) # Critical: b=0 would give A → ∞ + end + end + + @testset "Reproducibility of quasirandom" begin + b1 = 2.0 + bij1 = generate_bij(:quasirandom, 42, 3, b1; qmc_sampler = SobolSample()) + bij2 = generate_bij(:quasirandom, 42, 3, b1; qmc_sampler = SobolSample()) + @test bij1 ≈ bij2 + end + + @testset "Different indices give different values" begin + b1 = 1.0 + bij1 = generate_bij(:quasirandom, 1, 3, b1; qmc_sampler = HaltonSample()) + bij2 = generate_bij(:quasirandom, 2, 3, b1; qmc_sampler = HaltonSample()) + @test bij1 != bij2 + end +end + +@testset "generate_shift" begin + + @testset "Output range" begin + scale = 2.0 + for i in 1:100 + s = generate_shift(:quasirandom, i, 3, scale; qmc_sampler = HaltonSample()) + @test all(abs.(s) .<= scale + 1e-10) + end + end + + @testset "Correct dimension" begin + for dim in [1, 2, 3, 5] + s = generate_shift(:quasirandom, 1, dim, 1.0; qmc_sampler = HaltonSample()) + @test length(s) == dim + end + end +end + + +@testset "_generate_A_matrix" begin + + @testset "Symmetry" begin + w_list = [[1.0, -1.0], [1.0, 0.0]] + bij = [1.0, 2.0] + A = _generate_A_matrix(bij, w_list) + @test issymmetric(A) + end + + + @testset "Correct size" begin + w_list = [[1.0, -1.0], [1.0, 0.0]] + bij = [1.0, 2.0] + A = _generate_A_matrix(bij, w_list) + @test size(A) == (2, 2) + end + + @testset "Small bij gives large A elements" begin + w_list = [[1.0, -1.0]] + A_large_b = _generate_A_matrix([10.0], w_list) + A_small_b = _generate_A_matrix([0.1], w_list) + @test maximum(abs.(A_small_b)) > maximum(abs.(A_large_b)) + end +end + + +@testset "Physics: Variational principle" begin + # Energy should decrease monotonically as basis size increases + + masses = [1.0e15, 1.0] + Λmat = Λ(masses) + kin = KineticOperator(Λmat) + J, U = _jacobi_transform(masses) + w_raw = [U' * [1.0, -1.0]] + coulomb = CoulombOperator(-1.0, w_raw[1]) + ops = Operator[kin, coulomb] + + result = solve_ECG(ops, 20; scale = 1.5, verbose = false) + + # Check monotonic decrease (with some tolerance for numerical noise) + for i in 2:length(result.energies) + @test result.energies[i] <= result.energies[i-1] + 1e-10 + end +end + +@testset "Physics: Scale sensitivity" begin + # Demonstrate that scale matters for convergence + + masses = [1.0e15, 1.0] + Λmat = Λ(masses) + kin = KineticOperator(Λmat) + J, U = _jacobi_transform(masses) + w_raw = [U' * [1.0, -1.0]] + coulomb = CoulombOperator(-1.0, w_raw[1]) + ops = Operator[kin, coulomb] + + # Good scale for hydrogen + result_good = solve_ECG(ops, 15; scale = 1.5, verbose = false) + + # Bad scale (too small - Gaussians too narrow) + result_bad = solve_ECG(ops, 15; scale = 0.05, verbose = false) + + # Good scale should give better (lower) energy + @test result_good.ground_state < result_bad.ground_state +end + + +@testset "solve_ECG" begin + + @testset "Returns correct structure" begin + masses = [1.0e15, 1.0] + Λmat = Λ(masses) + kin = KineticOperator(Λmat) + J, U = _jacobi_transform(masses) + w_raw = [U' * [1.0, -1.0]] + coulomb = CoulombOperator(-1.0, w_raw[1]) + ops = Operator[kin, coulomb] + + result = solve_ECG(ops, 10; scale = 1.0, verbose = false) + + @test length(result.basis_functions) == result.n_basis + @test length(result.energies) == result.n_basis + @test result.ground_state == last(result.energies) + @test result.ground_state == minimum(result.energies) + end + + @testset "Respects max_attempts" begin + masses = [1.0e15, 1.0] + Λmat = Λ(masses) + kin = KineticOperator(Λmat) + J, U = _jacobi_transform(masses) + w_raw = [U' * [1.0, -1.0]] + coulomb = CoulombOperator(-1.0, w_raw[1]) + ops = Operator[kin, coulomb] + + # Request many basis functions but limit attempts + result = solve_ECG(ops, 1000; scale = 1.0, max_attempts = 50, verbose = false) + + @test result.n_basis <= 50 + end + + @testset "Handles linear dependence rejection" begin + masses = [1.0e15, 1.0] + Λmat = Λ(masses) + kin = KineticOperator(Λmat) + J, U = _jacobi_transform(masses) + w_raw = [U' * [1.0, -1.0]] + coulomb = CoulombOperator(-1.0, w_raw[1]) + ops = Operator[kin, coulomb] + + # Very strict threshold should cause rejections + result = solve_ECG(ops, 10; scale = 1.0, threshold = 0.5, verbose = false) + + # Should still produce valid results + @test result.n_basis >= 1 + @test isfinite(result.ground_state) + end +end From 1d646a310298ff1b9081859775c7696bc8b78342 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Thu, 11 Dec 2025 13:07:51 +0100 Subject: [PATCH 15/15] updated tests --- "Examples/td\316\274.jl" | 4 +- src/hamiltonian.jl | 16 +- src/sampling.jl | 6 +- test/runtests.jl | 2 +- test/test_hamiltonian.jl | 36 ++-- test/test_utils.jl | 419 +++++++++++++++++++++++++++++++++++++++ 6 files changed, 452 insertions(+), 31 deletions(-) create mode 100644 test/test_utils.jl diff --git "a/Examples/td\316\274.jl" "b/Examples/td\316\274.jl" index c388d98..b9b18fa 100644 --- "a/Examples/td\316\274.jl" +++ "b/Examples/td\316\274.jl" @@ -21,5 +21,5 @@ scale = default_scale(masses) result = solve_ECG(ops, 550, sampler = HaltonSample(); scale = scale) println("E ≈ ", result.ground_state) -a,b = correlation_function(result) -conv = FewBodyECG.convergence(result) \ No newline at end of file +a, b = correlation_function(result) +conv = FewBodyECG.convergence(result) diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl index 44e0a34..561ac71 100644 --- a/src/hamiltonian.jl +++ b/src/hamiltonian.jl @@ -150,8 +150,8 @@ function is_linearly_independent( end function default_scale(masses::Vector{<:Real}) - μ = minimum(masses[masses .< 1e10]) - return 1 / sqrt(μ) + μ = minimum(masses[masses .< 1.0e10]) + return 1 / sqrt(μ) end function solve_ECG( @@ -162,7 +162,7 @@ function solve_ECG( scale::Real = 0.2, threshold::Real = 0.95, max_attempts::Int = 10 * n, - max_condition::Real = 1e12, + max_condition::Real = 1.0e12, verbose::Bool = true ) @@ -207,14 +207,14 @@ function solve_ECG( # Check for NaN/Inf BEFORE eigensolve if any(!isfinite, H) - @warn "Hamiltonian contains NaN/Inf at step $(n_accepted+1), rejecting basis function" + @warn "Hamiltonian contains NaN/Inf at step $(n_accepted + 1), rejecting basis function" pop!(basis_fns) n_rejected += 1 continue end if any(!isfinite, S) - @warn "Overlap contains NaN/Inf at step $(n_accepted+1), rejecting basis function" + @warn "Overlap contains NaN/Inf at step $(n_accepted + 1), rejecting basis function" pop!(basis_fns) n_rejected += 1 continue @@ -223,7 +223,7 @@ function solve_ECG( # Check condition number cond_S = cond(Symmetric(S)) if cond_S > max_condition - @warn "Overlap poorly conditioned (κ=$cond_S) at step $(n_accepted+1), rejecting" + @warn "Overlap poorly conditioned (κ=$cond_S) at step $(n_accepted + 1), rejecting" pop!(basis_fns) n_rejected += 1 continue @@ -233,7 +233,7 @@ function solve_ECG( λs, Us = solve_generalized_eigenproblem(H, S; max_condition) catch e - @warn "Failed at step $(n_accepted+1): $e" + @warn "Failed at step $(n_accepted + 1): $e" pop!(basis_fns) # Remove problematic function n_rejected += 1 continue @@ -254,4 +254,4 @@ function solve_ECG( Emin = last(E_hist) @info "Optimization complete" E₀ = Emin n_basis = n_accepted return SolverResults(basis_fns, n_accepted, operators, method, sampler, b₁, Emin, E_hist, vecs_list) -end \ No newline at end of file +end diff --git a/src/sampling.jl b/src/sampling.jl index 248f421..060b6e7 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -6,8 +6,10 @@ function _qmc_point(i::Int, d::Int; sampler = HaltonSample()) return QuasiMonteCarlo.sample(i + 1, d, sampler)[:, end] end -function generate_bij(method::Symbol, i::Int, n_terms::Int, b1::Float64; - qmc_sampler = HaltonSample(), bmin = 0.02 * b1, bmax = 5.0 * b1) +function generate_bij( + method::Symbol, i::Int, n_terms::Int, b1::Float64; + qmc_sampler = HaltonSample(), bmin = 0.02 * b1, bmax = 5.0 * b1 + ) bmin > 0 || throw(ArgumentError("bmin must be > 0")) bmax > bmin || throw(ArgumentError("bmax must be > bmin")) u = method === :quasirandom ? _qmc_point(i, n_terms; sampler = qmc_sampler) : diff --git a/test/runtests.jl b/test/runtests.jl index d43d85d..0fa9fa8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,5 +10,5 @@ using FewBodyECG include("test_matrix_elements.jl") include("test_hamiltonian.jl") include("test_types.jl") - + include("test_utils.jl") end diff --git a/test/test_hamiltonian.jl b/test/test_hamiltonian.jl index 5c9dff5..29d51ff 100644 --- a/test/test_hamiltonian.jl +++ b/test/test_hamiltonian.jl @@ -444,11 +444,11 @@ end b1 = 1.0 bmin = 0.02 * b1 bmax = 5.0 * b1 - + for i in 1:100 bij = generate_bij(:quasirandom, i, 3, b1; qmc_sampler = HaltonSample()) - @test all(bij .>= bmin - 1e-10) - @test all(bij .<= bmax + 1e-10) + @test all(bij .>= bmin - 1.0e-10) + @test all(bij .<= bmax + 1.0e-10) end end @@ -481,7 +481,7 @@ end scale = 2.0 for i in 1:100 s = generate_shift(:quasirandom, i, 3, scale; qmc_sampler = HaltonSample()) - @test all(abs.(s) .<= scale + 1e-10) + @test all(abs.(s) .<= scale + 1.0e-10) end end @@ -522,7 +522,7 @@ end @testset "Physics: Variational principle" begin # Energy should decrease monotonically as basis size increases - + masses = [1.0e15, 1.0] Λmat = Λ(masses) kin = KineticOperator(Λmat) @@ -530,18 +530,18 @@ end w_raw = [U' * [1.0, -1.0]] coulomb = CoulombOperator(-1.0, w_raw[1]) ops = Operator[kin, coulomb] - + result = solve_ECG(ops, 20; scale = 1.5, verbose = false) - + # Check monotonic decrease (with some tolerance for numerical noise) for i in 2:length(result.energies) - @test result.energies[i] <= result.energies[i-1] + 1e-10 + @test result.energies[i] <= result.energies[i - 1] + 1.0e-10 end end @testset "Physics: Scale sensitivity" begin # Demonstrate that scale matters for convergence - + masses = [1.0e15, 1.0] Λmat = Λ(masses) kin = KineticOperator(Λmat) @@ -549,13 +549,13 @@ end w_raw = [U' * [1.0, -1.0]] coulomb = CoulombOperator(-1.0, w_raw[1]) ops = Operator[kin, coulomb] - + # Good scale for hydrogen result_good = solve_ECG(ops, 15; scale = 1.5, verbose = false) - + # Bad scale (too small - Gaussians too narrow) result_bad = solve_ECG(ops, 15; scale = 0.05, verbose = false) - + # Good scale should give better (lower) energy @test result_good.ground_state < result_bad.ground_state end @@ -571,9 +571,9 @@ end w_raw = [U' * [1.0, -1.0]] coulomb = CoulombOperator(-1.0, w_raw[1]) ops = Operator[kin, coulomb] - + result = solve_ECG(ops, 10; scale = 1.0, verbose = false) - + @test length(result.basis_functions) == result.n_basis @test length(result.energies) == result.n_basis @test result.ground_state == last(result.energies) @@ -588,10 +588,10 @@ end w_raw = [U' * [1.0, -1.0]] coulomb = CoulombOperator(-1.0, w_raw[1]) ops = Operator[kin, coulomb] - + # Request many basis functions but limit attempts result = solve_ECG(ops, 1000; scale = 1.0, max_attempts = 50, verbose = false) - + @test result.n_basis <= 50 end @@ -603,10 +603,10 @@ end w_raw = [U' * [1.0, -1.0]] coulomb = CoulombOperator(-1.0, w_raw[1]) ops = Operator[kin, coulomb] - + # Very strict threshold should cause rejections result = solve_ECG(ops, 10; scale = 1.0, threshold = 0.5, verbose = false) - + # Should still produce valid results @test result.n_basis >= 1 @test isfinite(result.ground_state) diff --git a/test/test_utils.jl b/test/test_utils.jl new file mode 100644 index 0000000..3811122 --- /dev/null +++ b/test/test_utils.jl @@ -0,0 +1,419 @@ +using Test +using LinearAlgebra +using FewBodyHamiltonians +using FewBodyECG +import FewBodyECG: _jacobi_transform, _generate_A_matrix, generate_bij, generate_shift +import FewBodyECG: ψ₀, convergence, correlation_function, SolverResults +using QuasiMonteCarlo + +function create_mock_solver_results(; + n_basis::Int = 5, + dim::Int = 2, + scale::Float64 = 1.0 + ) + # Create simple basis functions + basis_fns = GaussianBase[] + for i in 1:n_basis + A = Diagonal(fill(0.5 * i, dim)) + s = zeros(dim) + push!(basis_fns, Rank0Gaussian(Matrix(A), s)) + end + + # Create mock operators + K = KineticOperator(Diagonal(fill(0.5, dim))) + V = CoulombOperator(-1.0, [1.0; zeros(dim - 1)]) + operators = Operator[K, V] + + # Create mock energies (decreasing sequence) + energies = [-0.1 * i for i in 1:n_basis] + + # Create mock eigenvectors + eigenvectors = [randn(i, i) for i in 1:n_basis] + # Normalize columns + for i in 1:n_basis + for j in 1:i + eigenvectors[i][:, j] ./= norm(eigenvectors[i][:, j]) + end + end + + return SolverResults( + basis_fns, + n_basis, + operators, + :quasirandom, + HaltonSample(), + scale, + energies[end], + energies, + eigenvectors + ) +end + +@testset "SolverResults" begin + + @testset "Construction" begin + sr = create_mock_solver_results(n_basis = 3, dim = 2) + + @test sr.n_basis == 3 + @test length(sr.basis_functions) == 3 + @test length(sr.energies) == 3 + @test length(sr.eigenvectors) == 3 + @test sr.method == :quasirandom + @test sr.length_scale == 1.0 + end + + @testset "Ground state is last energy" begin + sr = create_mock_solver_results(n_basis = 5) + @test sr.ground_state == sr.energies[end] + end + + @testset "Operators stored correctly" begin + sr = create_mock_solver_results() + @test length(sr.operators) == 2 + @test sr.operators[1] isa KineticOperator + @test sr.operators[2] isa CoulombOperator + end + + @testset "Different samplers" begin + basis_fns = [Rank0Gaussian([1.0;;], [0.0])] + ops = Operator[KineticOperator([0.5;;])] + + sr_halton = SolverResults(basis_fns, 1, ops, :quasirandom, HaltonSample(), 1.0, -0.5, [-0.5], [ones(1, 1)]) + sr_sobol = SolverResults(basis_fns, 1, ops, :quasirandom, SobolSample(), 1.0, -0.5, [-0.5], [ones(1, 1)]) + + @test sr_halton.sampler isa HaltonSample + @test sr_sobol.sampler isa SobolSample + end +end + +@testset "ψ₀" begin + + @testset "Basic evaluation with coefficients" begin + # Single Gaussian: ψ = c * exp(-r'Ar + s'r) + A = [1.0 0.0; 0.0 1.0] + s = [0.0, 0.0] + g = Rank0Gaussian(A, s) + basis_fns = [g] + c = [1.0] + + r = [0.0, 0.0] + ψ_val = ψ₀(r, c, basis_fns) + + # At origin with s=0: exp(-0 + 0) = 1 + @test ψ_val ≈ 1.0 rtol = 1.0e-10 + end + + @testset "Gaussian decay" begin + A = [1.0 0.0; 0.0 1.0] + s = [0.0, 0.0] + g = Rank0Gaussian(A, s) + basis_fns = [g] + c = [1.0] + + ψ_origin = ψ₀([0.0, 0.0], c, basis_fns) + ψ_far = ψ₀([3.0, 3.0], c, basis_fns) + + # Should decay away from origin + @test abs(ψ_far) < abs(ψ_origin) + @test ψ_far ≈ exp(-18.0) rtol = 1.0e-10 # exp(-(3² + 3²)) + end + + @testset "Shift vector effect" begin + A = [1.0 0.0; 0.0 1.0] + s = [2.0, 0.0] # Shift in x-direction + g = Rank0Gaussian(A, s) + basis_fns = [g] + c = [1.0] + + # Maximum should be shifted + ψ_origin = ψ₀([0.0, 0.0], c, basis_fns) + ψ_shifted = ψ₀([1.0, 0.0], c, basis_fns) # Closer to maximum + + @test abs(ψ_shifted) > abs(ψ_origin) + end + + @testset "Linear combination" begin + A1 = [1.0 0.0; 0.0 1.0] + A2 = [2.0 0.0; 0.0 2.0] + s = [0.0, 0.0] + g1 = Rank0Gaussian(A1, s) + g2 = Rank0Gaussian(A2, s) + basis_fns = [g1, g2] + + c = [0.5, 0.5] + r = [0.0, 0.0] + + ψ_val = ψ₀(r, c, basis_fns) + + # At origin: 0.5 * 1 + 0.5 * 1 = 1 + @test ψ_val ≈ 1.0 rtol = 1.0e-10 + end + + @testset "Negative coefficients" begin + A = [1.0 0.0; 0.0 1.0] + s = [0.0, 0.0] + g = Rank0Gaussian(A, s) + basis_fns = [g] + + c_pos = [1.0] + c_neg = [-1.0] + r = [0.0, 0.0] + + @test ψ₀(r, c_pos, basis_fns) ≈ -ψ₀(r, c_neg, basis_fns) rtol = 1.0e-10 + end + + @testset "With SolverResults" begin + sr = create_mock_solver_results(n_basis = 3, dim = 2) + + # Should not throw + r = [0.5, 0.5] + ψ_val = ψ₀(r, sr; state = 1) + + @test isfinite(ψ_val) + end + + @testset "Different states" begin + sr = create_mock_solver_results(n_basis = 5, dim = 2) + r = [0.1, 0.1] + + # Different states should generally give different values + ψ_1 = ψ₀(r, sr; state = 1) + ψ_2 = ψ₀(r, sr; state = 2) + + @test isfinite(ψ_1) + @test isfinite(ψ_2) + # They might be equal by chance, but usually won't be + end + + @testset "1D case" begin + A = [2.0;;] + s = [0.0] + g = Rank0Gaussian(A, s) + basis_fns = [g] + c = [1.0] + + r = [1.0] + ψ_val = ψ₀(r, c, basis_fns) + + @test ψ_val ≈ exp(-2.0) rtol = 1.0e-10 + end +end + +@testset "convergence" begin + + @testset "Returns correct range and energies" begin + sr = create_mock_solver_results(n_basis = 10) + + indices, energies = convergence(sr) + + @test indices == 1:10 + @test energies == sr.energies + @test length(indices) == length(energies) + end + + @testset "Single basis function" begin + sr = create_mock_solver_results(n_basis = 1) + + indices, energies = convergence(sr) + + @test indices == 1:1 + @test length(energies) == 1 + end + + @testset "Energies are same object" begin + sr = create_mock_solver_results(n_basis = 5) + + _, energies = convergence(sr) + + # Should be the same array (not a copy) + @test energies === sr.energies + end +end + +@testset "correlation_function" begin + + @testset "Output dimensions" begin + sr = create_mock_solver_results(n_basis = 5, dim = 2) + + r_grid, ρ = correlation_function(sr; npoints = 100) + + @test length(r_grid) == 100 + @test length(ρ) == 100 + end + + @testset "Grid range" begin + sr = create_mock_solver_results(n_basis = 5, dim = 2) + + rmin, rmax = 0.5, 5.0 + r_grid, _ = correlation_function(sr; rmin = rmin, rmax = rmax, npoints = 50) + + @test first(r_grid) ≈ rmin + @test last(r_grid) ≈ rmax + end + + @testset "Non-negative density" begin + sr = create_mock_solver_results(n_basis = 5, dim = 2) + + _, ρ = correlation_function(sr; normalize = false) + + # r²|ψ|² should always be non-negative + @test all(ρ .>= 0) + end + + @testset "Normalization" begin + sr = create_mock_solver_results(n_basis = 5, dim = 2) + + r_grid, ρ_norm = correlation_function(sr; normalize = true, npoints = 500) + _, ρ_unnorm = correlation_function(sr; normalize = false, npoints = 500) + + # Normalized should integrate to ~1 + dr = r_grid[2] - r_grid[1] + integral_norm = sum(ρ_norm) * dr + + # Check that normalization changed something (unless already normalized) + if sum(ρ_unnorm) * dr > 1.0e-10 + @test integral_norm ≈ 1.0 rtol = 0.1 # Rough due to trapezoidal rule + end + end + + @testset "coord_index selection" begin + sr = create_mock_solver_results(n_basis = 3, dim = 3) + + # Should work for all valid indices + for idx in 1:3 + r_grid, ρ = correlation_function(sr; coord_index = idx) + @test length(r_grid) > 0 + @test all(isfinite, ρ) + end + end + + @testset "Invalid coord_index" begin + sr = create_mock_solver_results(n_basis = 3, dim = 2) + + @test_throws ArgumentError correlation_function(sr; coord_index = 0) + @test_throws ArgumentError correlation_function(sr; coord_index = 3) + @test_throws ArgumentError correlation_function(sr; coord_index = -1) + end + + @testset "Returns Vector not Range" begin + sr = create_mock_solver_results(n_basis = 3, dim = 2) + + r_grid, _ = correlation_function(sr) + + @test r_grid isa Vector + end + + @testset "Finite values" begin + sr = create_mock_solver_results(n_basis = 5, dim = 2) + + r_grid, ρ = correlation_function(sr) + + @test all(isfinite, r_grid) + @test all(isfinite, ρ) + end + + @testset "Custom npoints" begin + sr = create_mock_solver_results(n_basis = 3, dim = 2) + + for np in [10, 100, 1000] + r_grid, ρ = correlation_function(sr; npoints = np) + @test length(r_grid) == np + @test length(ρ) == np + end + end +end + + +@testset "Integration: Utils with real solver" begin + + @testset "Hydrogen atom utilities" begin + # Set up hydrogen atom + masses = [1.0e15, 1.0] + Λmat = Λ(masses) + kin = KineticOperator(Λmat) + J, U = _jacobi_transform(masses) + w_raw = [U' * [1.0, -1.0]] + coulomb = CoulombOperator(-1.0, w_raw[1]) + ops = Operator[kin, coulomb] + + result = solve_ECG(ops, 15; scale = 1.5, verbose = false) + + # Test convergence + indices, energies = convergence(result) + @test length(indices) == result.n_basis + @test energies[end] == result.ground_state + + # Test wavefunction evaluation + r = [0.5] + ψ_val = ψ₀(r, result; state = 1) + @test isfinite(ψ_val) + + r_grid, ρ = correlation_function(result; npoints = 100) + @test all(ρ .>= 0) + @test length(r_grid) == 100 + end + + @testset "Three-body utilities" begin + # Set up three-body system + masses = [1000.0, 1000.0, 1.0] + Λmat = Λ(masses) + kin = KineticOperator(Λmat) + J, U = _jacobi_transform(masses) + + w_list = [[1, -1, 0], [1, 0, -1], [0, 1, -1]] + w_raw = [U' * w for w in w_list] + coeffs = [+1.0, -1.0, -1.0] + coulomb_ops = [CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw)] + ops = Operator[kin; coulomb_ops...] + + result = solve_ECG(ops, 10; scale = 1.0, verbose = false) + + _, energies = convergence(result) + for i in 2:length(energies) + @test energies[i] <= energies[i - 1] + 1.0e-10 + end + + for coord_idx in 1:2 + r_grid, ρ = correlation_function(result; coord_index = coord_idx) + @test all(isfinite, ρ) + end + end +end + +@testset "Edge cases" begin + + @testset "Very small basis" begin + basis_fns = [Rank0Gaussian([1.0;;], [0.0])] + ops = Operator[KineticOperator([0.5;;])] + eigvecs = [ones(1, 1)] + + sr = SolverResults( + basis_fns, 1, ops, :quasirandom, HaltonSample(), + 1.0, -0.5, [-0.5], eigvecs + ) + + # All utilities should work + @test ψ₀([0.0], sr) ≈ 1.0 + @test convergence(sr) == (1:1, [-0.5]) + + r_grid, ρ = correlation_function(sr) + @test length(r_grid) > 0 + end + + @testset "Large coordinates" begin + sr = create_mock_solver_results(n_basis = 3, dim = 2) + + r_large = [100.0, 100.0] + ψ_val = ψ₀(r_large, sr) + + @test isfinite(ψ_val) + @test abs(ψ_val) < 1.0e-10 + end + + @testset "Zero at correlation function boundaries" begin + sr = create_mock_solver_results(n_basis = 5, dim = 2) + + r_grid, ρ = correlation_function(sr; rmin = 1.0e-6, rmax = 1.0) + @test ρ[1] ≈ 0 atol = 1.0e-10 # r²|ψ|² → 0 as r → 0 + end +end