From e02765d0e4244365635d79f8de65704d6cc61357 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Fri, 10 Jul 2026 12:48:41 +0200 Subject: [PATCH 1/5] added spin-orbit example --- docs/Manifest.toml | 2 +- docs/make.jl | 3 + docs/src/API.md | 11 + docs/src/examples/gaussian_well.md | 5 - docs/src/examples/harmonium.md | 53 +++++ docs/src/examples/helium.md | 4 +- docs/src/examples/spin_orbit.md | 100 +++++++++ docs/src/examples/tdmu.md | 10 +- docs/src/examples/three_body_gaussian.md | 54 +++++ examples/harmonium.jl | 42 ++++ examples/helium.jl | 4 +- examples/spin_orbit.jl | 87 ++++++++ examples/tdmu.jl | 4 +- examples/three_body_gaussian.jl | 43 ++++ src/FewBodyECG.jl | 11 +- src/gradient.jl | 29 +-- src/linalg.jl | 45 ++-- src/matrix_elements.jl | 151 +++++++++++-- src/operators.jl | 16 +- src/sampling.jl | 8 +- src/solution.jl | 2 +- src/solve.jl | 4 +- src/state.jl | 7 +- src/types.jl | 202 ++++++++++++++++-- src/utils/convergence.jl | 11 + src/utils/observables.jl | 34 +++ src/{recipes.jl => utils/plotting.jl} | 15 +- src/{observables.jl => utils/wavefunction.jl} | 15 +- test/runtests.jl | 3 + test/test_hamiltonian.jl | 20 +- test/test_harmonic_oscillator.jl | 35 +++ test/test_harmonium.jl | 20 ++ test/test_matrix_elements.jl | 68 +++++- test/test_observables.jl | 41 +++- test/test_operators.jl | 16 ++ test/test_sampling.jl | 2 +- test/test_solve.jl | 5 +- test/test_spin_orbit.jl | 62 ++++++ test/test_types.jl | 12 +- test/test_variational.jl | 20 +- 40 files changed, 1140 insertions(+), 136 deletions(-) create mode 100644 docs/src/examples/harmonium.md create mode 100644 docs/src/examples/spin_orbit.md create mode 100644 docs/src/examples/three_body_gaussian.md create mode 100644 examples/harmonium.jl create mode 100644 examples/spin_orbit.jl create mode 100644 examples/three_body_gaussian.jl create mode 100644 src/utils/convergence.jl create mode 100644 src/utils/observables.jl rename src/{recipes.jl => utils/plotting.jl} (71%) rename src/{observables.jl => utils/wavefunction.jl} (58%) create mode 100644 test/test_harmonic_oscillator.jl create mode 100644 test/test_harmonium.jl create mode 100644 test/test_spin_orbit.jl diff --git a/docs/Manifest.toml b/docs/Manifest.toml index cc0c7ae..77d1c3a 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -389,7 +389,7 @@ version = "8.0.1+0" [[deps.FewBodyECG]] deps = ["Antique", "FewBodyHamiltonians", "ForwardDiff", "LinearAlgebra", "OptimKit", "QuasiMonteCarlo", "RecipesBase", "SpecialFunctions"] -path = ".." +path = "/Users/pzb464/Documents/FewBodyECG.jl" uuid = "083b1810-24a1-4a79-9a41-145bb2bb8ceb" version = "2.0.0" diff --git a/docs/make.jl b/docs/make.jl index 4932119..ba9d636 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -25,6 +25,9 @@ makedocs( "tdmu" => "examples/tdmu.md", "H2+ (non-BO)" => "examples/h2plus.md", "Gaussian wells" => "examples/gaussian_well.md", + "Hooke's atom" => "examples/harmonium.md", + "Three-body Gaussian" => "examples/three_body_gaussian.md", + "Spin-orbit doublet" => "examples/spin_orbit.md", "Workflow" => "examples/workflow.md", ], "API" => "API.md", diff --git a/docs/src/API.md b/docs/src/API.md index 1517ca0..b3b5d8e 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -9,6 +9,15 @@ Operator KineticOperator CoulombOperator GaussianOperator +OscillatorOperator +ManyBodyGaussianOperator +GaussianTensorOperator +GaussianSpinOrbitOperator +SpinProjection +up +down +SpinState +SpinGaussian GaussianBase Rank0Gaussian Rank1Gaussian @@ -37,8 +46,10 @@ ConvergenceReport StageResult converged energies +convergence wavefunction Wavefunction +radial_profile ``` ## Power-user layer diff --git a/docs/src/examples/gaussian_well.md b/docs/src/examples/gaussian_well.md index b3bfa17..0da2f93 100644 --- a/docs/src/examples/gaussian_well.md +++ b/docs/src/examples/gaussian_well.md @@ -2,11 +2,6 @@ EditURL = "../../../examples/gaussian_well.jl" ``` -# Gaussian well - -A finite-range Gaussian attraction, V(r) = -V0 exp(-gamma r^2), is useful for -short-range model interactions and has a binding threshold unlike Coulomb. - ````@example gaussian_well using FewBodyECG using Plots diff --git a/docs/src/examples/harmonium.md b/docs/src/examples/harmonium.md new file mode 100644 index 0000000..90bf0a1 --- /dev/null +++ b/docs/src/examples/harmonium.md @@ -0,0 +1,53 @@ +```@meta +EditURL = "../../../examples/harmonium.jl" +``` + +# Hooke's atom (harmonium) + +Two electrons in a common harmonic trap, repelling through the Coulomb +interaction: + +```math +H = -\tfrac12\nabla_1^2 - \tfrac12\nabla_2^2 + + \tfrac12\omega^2 (r_1^2 + r_2^2) + \frac{1}{r_{12}}. +``` + +For the trap frequency ``\omega = \tfrac12`` the ground-state energy is known +in closed form (Taut, *Phys. Rev. A* **48**, 3561 (1993)): **exactly 2 Ha**. +The trap is supplied by an `OscillatorOperator` between a fixed heavy centre +and each electron (coefficient ``\tfrac12\omega^2``); the only Coulomb term is +the electron–electron repulsion. The system is spherically symmetric and real, +so the plain stochastic `SVM` solver applies. + +````@example harmonium +using FewBodyECG +using Plots + +ω = 0.5 +masses = [1.0e15, 1.0, 1.0] # heavy trap centre + two electrons + +ops = Operators(masses) +ops += "Kinetic" +ops += ("Oscillator", 1, 2, 0.5 * ω^2) # ½ω² r₁² (trap on electron 1) +ops += ("Oscillator", 1, 3, 0.5 * ω^2) # ½ω² r₂² (trap on electron 2) +ops += ("Coulomb", 2, 3, 1.0) # +1/r₁₂ electron–electron repulsion + +exact = 2.0 # Taut 1993, ω = 1/2 + +sol = solve(ops, SVM(basis = 80, candidates = 40, scale = 2.0)) +println("Hooke's atom E₀ = ", sol.E₀, " Ha (Taut exact ", exact, ", Δ = ", sol.E₀ - exact, ")") +println("variational upper bound respected: ", sol.E₀ ≥ exact) +sol + +# Convergence toward the exact energy +plot(sol, exact) + +# Ground-state radial density of one electron +plot(wavefunction(sol); coord = 1, rmax = 8.0) +plot(convergence(sol)) +```` + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + diff --git a/docs/src/examples/helium.md b/docs/src/examples/helium.md index a6fc284..e164c6d 100644 --- a/docs/src/examples/helium.md +++ b/docs/src/examples/helium.md @@ -16,7 +16,7 @@ helium += "Kinetic" helium += "Coulomb" he_ref = -2.9037 -he = solve(helium, SVM(basis = 35, candidates = 25, scale = 1.0)) +he = solve(helium, SVM(basis = 50, candidates = 25, scale = 1.0)) println("Helium E0 = ", he.E₀, " Ha (reference ", he_ref, ", Δ = ", he.E₀ - he_ref, ")") hminus = Operators([1.0e15, 1.0, 1.0], [+1.0, -1.0, -1.0]) @@ -24,7 +24,7 @@ hminus += "Kinetic" hminus += "Coulomb" hm_ref = -0.52775 -hm = solve(hminus, SVM(basis = 30, candidates = 20, scale = 1.0)) +hm = solve(hminus, SVM(basis = 50, candidates = 20, scale = 1.0)) println("H- E0 = ", hm.E₀, " Ha (reference ", hm_ref, ", Δ = ", hm.E₀ - hm_ref, ")") plot(he, he_ref) diff --git a/docs/src/examples/spin_orbit.md b/docs/src/examples/spin_orbit.md new file mode 100644 index 0000000..979b4f0 --- /dev/null +++ b/docs/src/examples/spin_orbit.md @@ -0,0 +1,100 @@ +```@meta +EditURL = "../../../examples/spin_orbit.jl" +``` + +# Spin-orbit doublet — the complex Hermitian path + +Spin-dependent interactions make the Hamiltonian **complex Hermitian**, and +FewBodyECG assembles and diagonalises it directly. Here two spin-½ particles +share a p-wave-like spatial manifold — Gaussians shifted along `x`, `y`, `z`, +each dressed with the spin-triplet state `|↑↑⟩` — bound by a central Gaussian +attraction. A [`GaussianSpinOrbitOperator`](@ref) couples the relative +orbital motion to the total spin `S₁ + S₂`. + +With the coupling off, the two upper orbital states are degenerate. Turning it +on lifts the degeneracy: the `m_ℓ = ±1` combinations shift in opposite +directions, giving a doublet whose splitting grows **linearly** with the +coupling strength — a model fine-structure splitting. + +````@example spin_orbit +using FewBodyECG +using LinearAlgebra +using Plots + +masses = [1.0, 1.0] +ops = Operators(masses) +ops += "Kinetic" +ops += ("Gaussian", 1, 2, -4.0, 0.4) # central binding well +terms = ops.terms +w = only(op.w for op in terms if op isa GaussianOperator) + +# p-wave-like manifold: unit-shifted Gaussians along x, y, z with |↑↑⟩ spin +shift(dir) = (v = zeros(1, 3); v[1, dir] = 0.8; v) +basis = BasisSet( + [ + SpinGaussian(Rank0Gaussian([0.4;;], shift(d)), SpinState([up, up])) for d in 1:3 + ] +) + +H₀ = build_hamiltonian_matrix(basis, terms) # real, SO-free +Hₛₒ = build_hamiltonian_matrix(basis, FewBodyECG.Operator[GaussianSpinOrbitOperator(1.0, 0.4, w, 1, 2)]) +S = build_overlap_matrix(basis) + +println("H₀ eltype = ", eltype(H₀), " (spin-orbit-free, real)") +println("Hₛₒ eltype = ", eltype(Hₛₒ), " Hermitian? ", isapprox(Hₛₒ, Hₛₒ')) + +# Sweep the spin-orbit strength λ and collect the (real) spectrum +λs = range(0, 2, length = 21) +levels = map(λs) do λ + E, _ = solve_generalized_eigenproblem(H₀ .+ λ .* Hₛₒ, S) + sort(real.(E)) +end +spectrum = reduce(hcat, levels)' # rows = λ, columns = levels + +# The upper two levels form the spin-orbit doublet; its gap is linear in λ. +gap = spectrum[:, 3] .- spectrum[:, 2] +println("doublet gap at λ = 0, 1, 2: ", round.(gap[[1, 11, 21]], digits = 4)) + +plot( + λs, spectrum; + xlabel = "spin-orbit coupling λ", ylabel = "E (Ha)", + label = ["ground" "doublet⁻" "doublet⁺"], linewidth = 2, + title = "Spin-orbit doublet splitting", +) +```` + +## Convergence of the manifold + +The three-function manifold above is the minimal one. Enriching each Cartesian +direction with a geometric ladder of widths gives a sequence of **nested** +bases, so the ground-state energy (here at fixed coupling `λ = 1`) descends +monotonically toward the variational limit. + +````@example spin_orbit +αs = exp10.(range(-0.8, 0.6, length = 8)) # fixed ladder ⇒ nested bases +sg(α, d) = SpinGaussian(Rank0Gaussian([α;;], shift(d)), SpinState([up, up])) +λ = 1.0 +so = GaussianSpinOrbitOperator(1.0, 0.4, w, 1, 2) + +sizes = Int[] +ground = Float64[] +for k in 1:length(αs) + b = BasisSet([sg(α, d) for α in αs[1:k] for d in 1:3]) + H = build_hamiltonian_matrix(b, terms) .+ λ .* build_hamiltonian_matrix(b, FewBodyECG.Operator[so]) + E, _ = solve_generalized_eigenproblem(H, build_overlap_matrix(b)) + push!(sizes, 3k) + push!(ground, minimum(real.(E))) +end +println("converged ground-state energy: ", ground[end], " Ha") + +plot( + sizes, ground; + xlabel = "basis size", ylabel = "E₀ (Ha)", label = "ground state (λ = 1)", + marker = :circle, linewidth = 2, title = "Basis-size convergence", +) +```` + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + diff --git a/docs/src/examples/tdmu.md b/docs/src/examples/tdmu.md index d8fc39f..215284d 100644 --- a/docs/src/examples/tdmu.md +++ b/docs/src/examples/tdmu.md @@ -2,12 +2,6 @@ EditURL = "../../../examples/tdmu.jl" ``` -# tdmu muonic molecular ion - -The tdmu ion is deeply bound because the muon is much heavier than an -electron. A loose stochastic run is enough to show the energy scale; use a -tighter `tol` and larger basis for production values near -111.36444 Ha. - ````@example tdmu using FewBodyECG using Plots @@ -18,8 +12,8 @@ ops += "Coulomb" sol = solve( ops, - SVM(basis = 40, candidates = 25, scale = 0.03); - tol = 1.0e-2, + SVM(basis = 100, candidates = 25, scale = 0.03); + tol = 1.0e-3, window = 10, ) sol diff --git a/docs/src/examples/three_body_gaussian.md b/docs/src/examples/three_body_gaussian.md new file mode 100644 index 0000000..f342e41 --- /dev/null +++ b/docs/src/examples/three_body_gaussian.md @@ -0,0 +1,54 @@ +```@meta +EditURL = "../../../examples/three_body_gaussian.jl" +``` + +# Three-body Gaussian model with a many-body regulator + +A heavy centre binds two identical light particles through pairwise Gaussian +attractions. A repulsive many-body Gaussian `exp(-rᵀWr)` acts as a regulator +that lifts the energy — a physically useful interaction that stays compatible +with the rank-0 stochastic solver. + +````@example three_body_gaussian +using FewBodyECG +using Plots + +masses = [1.0e15, 1.0, 1.0] + +build() = begin + o = Operators(masses) + o += "Kinetic" + o += ("Gaussian", 1, 2, -4.0, 0.5) # centre–particle-2 attraction + o += ("Gaussian", 1, 3, -4.0, 0.5) # centre–particle-3 attraction + o +end + +# Without the regulator (λ = 0) +ops₀ = build() +sol₀ = solve(ops₀, SVM(basis = 30, candidates = 20, scale = 1.0)) + +# With a repulsive many-body Gaussian regulator (λ > 0, W positive-definite) +W = [0.6 0.1; 0.1 0.6] +λ = 3.0 +ops₁ = build() +ops₁ += ManyBodyGaussianOperator(λ, W) +sol₁ = solve(ops₁, SVM(basis = 30, candidates = 20, scale = 1.0)) + +ΔE = sol₁.E₀ - sol₀.E₀ +println("E₀ (λ = 0) = ", sol₀.E₀, " Ha") +println("E₀ (λ = ", λ, ") = ", sol₁.E₀, " Ha") +println("regulator-induced shift ΔE = ", ΔE, " Ha") + +# Convergence and radial profile through the utilities +steps, history = convergence(sol₁) +pc = plot(steps, history; xlabel = "step", ylabel = "E (Ha)", label = "with regulator") +plot!(pc, convergence(sol₀)...; label = "no regulator") +pc + +plot(wavefunction(sol₁); coord = 1, rmax = 6.0) +```` + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + diff --git a/examples/harmonium.jl b/examples/harmonium.jl new file mode 100644 index 0000000..0f9aef8 --- /dev/null +++ b/examples/harmonium.jl @@ -0,0 +1,42 @@ +# # Hooke's atom (harmonium) +# +# Two electrons in a common harmonic trap, repelling through the Coulomb +# interaction: +# +# ```math +# H = -\tfrac12\nabla_1^2 - \tfrac12\nabla_2^2 +# + \tfrac12\omega^2 (r_1^2 + r_2^2) + \frac{1}{r_{12}}. +# ``` +# +# For the trap frequency ``\omega = \tfrac12`` the ground-state energy is known +# in closed form (Taut, *Phys. Rev. A* **48**, 3561 (1993)): **exactly 2 Ha**. +# The trap is supplied by an `OscillatorOperator` between a fixed heavy centre +# and each electron (coefficient ``\tfrac12\omega^2``); the only Coulomb term is +# the electron–electron repulsion. The system is spherically symmetric and real, +# so the plain stochastic `SVM` solver applies. + +using FewBodyECG +using Plots + +ω = 0.5 +masses = [1.0e15, 1.0, 1.0] # heavy trap centre + two electrons + +ops = Operators(masses) +ops += "Kinetic" +ops += ("Oscillator", 1, 2, 0.5 * ω^2) # ½ω² r₁² (trap on electron 1) +ops += ("Oscillator", 1, 3, 0.5 * ω^2) # ½ω² r₂² (trap on electron 2) +ops += ("Coulomb", 2, 3, 1.0) # +1/r₁₂ electron–electron repulsion + +exact = 2.0 # Taut 1993, ω = 1/2 + +sol = solve(ops, SVM(basis = 80, candidates = 40, scale = 2.0)) +println("Hooke's atom E₀ = ", sol.E₀, " Ha (Taut exact ", exact, ", Δ = ", sol.E₀ - exact, ")") +println("variational upper bound respected: ", sol.E₀ ≥ exact) +sol + +## Convergence toward the exact energy +plot(sol, exact) + +## Ground-state radial density of one electron +plot(wavefunction(sol); coord = 1, rmax = 8.0) +plot(convergence(sol)) \ No newline at end of file diff --git a/examples/helium.jl b/examples/helium.jl index b2eb5ff..69efd65 100644 --- a/examples/helium.jl +++ b/examples/helium.jl @@ -11,7 +11,7 @@ helium += "Kinetic" helium += "Coulomb" he_ref = -2.9037 -he = solve(helium, SVM(basis = 35, candidates = 25, scale = 1.0)) +he = solve(helium, SVM(basis = 50, candidates = 25, scale = 1.0)) println("Helium E0 = ", he.E₀, " Ha (reference ", he_ref, ", Δ = ", he.E₀ - he_ref, ")") hminus = Operators([1.0e15, 1.0, 1.0], [+1.0, -1.0, -1.0]) @@ -19,7 +19,7 @@ hminus += "Kinetic" hminus += "Coulomb" hm_ref = -0.52775 -hm = solve(hminus, SVM(basis = 30, candidates = 20, scale = 1.0)) +hm = solve(hminus, SVM(basis = 50, candidates = 20, scale = 1.0)) println("H- E0 = ", hm.E₀, " Ha (reference ", hm_ref, ", Δ = ", hm.E₀ - hm_ref, ")") plot(he, he_ref) diff --git a/examples/spin_orbit.jl b/examples/spin_orbit.jl new file mode 100644 index 0000000..99c592d --- /dev/null +++ b/examples/spin_orbit.jl @@ -0,0 +1,87 @@ +# # Spin-orbit doublet — the complex Hermitian path +# +# Spin-dependent interactions make the Hamiltonian **complex Hermitian**, and +# FewBodyECG assembles and diagonalises it directly. Here two spin-½ particles +# share a p-wave-like spatial manifold — Gaussians shifted along `x`, `y`, `z`, +# each dressed with the spin-triplet state `|↑↑⟩` — bound by a central Gaussian +# attraction. A [`GaussianSpinOrbitOperator`](@ref) couples the relative +# orbital motion to the total spin `S₁ + S₂`. +# +# With the coupling off, the two upper orbital states are degenerate. Turning it +# on lifts the degeneracy: the `m_ℓ = ±1` combinations shift in opposite +# directions, giving a doublet whose splitting grows **linearly** with the +# coupling strength — a model fine-structure splitting. + +using FewBodyECG +using LinearAlgebra +using Plots + +masses = [1.0, 1.0] +ops = Operators(masses) +ops += "Kinetic" +ops += ("Gaussian", 1, 2, -4.0, 0.4) # central binding well +terms = ops.terms +w = only(op.w for op in terms if op isa GaussianOperator) + +## p-wave-like manifold: unit-shifted Gaussians along x, y, z with |↑↑⟩ spin +shift(dir) = (v = zeros(1, 3); v[1, dir] = 0.8; v) +basis = BasisSet( + [ + SpinGaussian(Rank0Gaussian([0.4;;], shift(d)), SpinState([up, up])) for d in 1:3 + ] +) + +H₀ = build_hamiltonian_matrix(basis, terms) # real, SO-free +Hₛₒ = build_hamiltonian_matrix(basis, FewBodyECG.Operator[GaussianSpinOrbitOperator(1.0, 0.4, w, 1, 2)]) +S = build_overlap_matrix(basis) + +println("H₀ eltype = ", eltype(H₀), " (spin-orbit-free, real)") +println("Hₛₒ eltype = ", eltype(Hₛₒ), " Hermitian? ", isapprox(Hₛₒ, Hₛₒ')) + +## Sweep the spin-orbit strength λ and collect the (real) spectrum +λs = range(0, 2, length = 21) +levels = map(λs) do λ + E, _ = solve_generalized_eigenproblem(H₀ .+ λ .* Hₛₒ, S) + sort(real.(E)) +end +spectrum = reduce(hcat, levels)' # rows = λ, columns = levels + +## The upper two levels form the spin-orbit doublet; its gap is linear in λ. +gap = spectrum[:, 3] .- spectrum[:, 2] +println("doublet gap at λ = 0, 1, 2: ", round.(gap[[1, 11, 21]], digits = 4)) + +plot( + λs, spectrum; + xlabel = "spin-orbit coupling λ", ylabel = "E (Ha)", + label = ["ground" "doublet⁻" "doublet⁺"], linewidth = 2, + title = "Spin-orbit doublet splitting", +) + +# ## Convergence of the manifold +# +# The three-function manifold above is the minimal one. Enriching each Cartesian +# direction with a geometric ladder of widths gives a sequence of **nested** +# bases, so the ground-state energy (here at fixed coupling `λ = 1`) descends +# monotonically toward the variational limit. + +αs = exp10.(range(-0.8, 0.6, length = 8)) # fixed ladder ⇒ nested bases +sg(α, d) = SpinGaussian(Rank0Gaussian([α;;], shift(d)), SpinState([up, up])) +λ = 1.0 +so = GaussianSpinOrbitOperator(1.0, 0.4, w, 1, 2) + +sizes = Int[] +ground = Float64[] +for k in 1:length(αs) + b = BasisSet([sg(α, d) for α in αs[1:k] for d in 1:3]) + H = build_hamiltonian_matrix(b, terms) .+ λ .* build_hamiltonian_matrix(b, FewBodyECG.Operator[so]) + E, _ = solve_generalized_eigenproblem(H, build_overlap_matrix(b)) + push!(sizes, 3k) + push!(ground, minimum(real.(E))) +end +println("converged ground-state energy: ", ground[end], " Ha") + +plot( + sizes, ground; + xlabel = "basis size", ylabel = "E₀ (Ha)", label = "ground state (λ = 1)", + marker = :circle, linewidth = 2, title = "Basis-size convergence", +) diff --git a/examples/tdmu.jl b/examples/tdmu.jl index 0e9d8d1..1ec0fb9 100644 --- a/examples/tdmu.jl +++ b/examples/tdmu.jl @@ -7,8 +7,8 @@ ops += "Coulomb" sol = solve( ops, - SVM(basis = 40, candidates = 25, scale = 0.03); - tol = 1.0e-2, + SVM(basis = 100, candidates = 25, scale = 0.03); + tol = 1.0e-3, window = 10, ) sol diff --git a/examples/three_body_gaussian.jl b/examples/three_body_gaussian.jl new file mode 100644 index 0000000..7fad632 --- /dev/null +++ b/examples/three_body_gaussian.jl @@ -0,0 +1,43 @@ +# # Three-body Gaussian model with a many-body regulator +# +# A heavy centre binds two identical light particles through pairwise Gaussian +# attractions. A repulsive many-body Gaussian `exp(-rᵀWr)` acts as a regulator +# that lifts the energy — a physically useful interaction that stays compatible +# with the rank-0 stochastic solver. + +using FewBodyECG +using Plots + +masses = [1.0e15, 1.0, 1.0] + +build() = begin + o = Operators(masses) + o += "Kinetic" + o += ("Gaussian", 1, 2, -4.0, 0.5) # centre–particle-2 attraction + o += ("Gaussian", 1, 3, -4.0, 0.5) # centre–particle-3 attraction + o +end + +## Without the regulator (λ = 0) +ops₀ = build() +sol₀ = solve(ops₀, SVM(basis = 30, candidates = 20, scale = 1.0)) + +## With a repulsive many-body Gaussian regulator (λ > 0, W positive-definite) +W = [0.6 0.1; 0.1 0.6] +λ = 3.0 +ops₁ = build() +ops₁ += ManyBodyGaussianOperator(λ, W) +sol₁ = solve(ops₁, SVM(basis = 30, candidates = 20, scale = 1.0)) + +ΔE = sol₁.E₀ - sol₀.E₀ +println("E₀ (λ = 0) = ", sol₀.E₀, " Ha") +println("E₀ (λ = ", λ, ") = ", sol₁.E₀, " Ha") +println("regulator-induced shift ΔE = ", ΔE, " Ha") + +## Convergence and radial profile through the utilities +steps, history = convergence(sol₁) +pc = plot(steps, history; xlabel = "step", ylabel = "E (Ha)", label = "with regulator") +plot!(pc, convergence(sol₀)...; label = "no regulator") +pc + +plot(wavefunction(sol₁); coord = 1, rmax = 6.0) diff --git a/src/FewBodyECG.jl b/src/FewBodyECG.jl index bfc1f9b..dd0a27e 100644 --- a/src/FewBodyECG.jl +++ b/src/FewBodyECG.jl @@ -15,12 +15,15 @@ const Operator = FewBodyHamiltonians.Operator # system building export Operators, coulomb_weights, Operator, KineticOperator, CoulombOperator, GaussianOperator, + OscillatorOperator, ManyBodyGaussianOperator, + GaussianTensorOperator, GaussianSpinOrbitOperator, + SpinProjection, up, down, SpinState, SpinGaussian, GaussianBase, Rank0Gaussian, Rank1Gaussian, Rank2Gaussian, BasisSet # solving export solve, SVM, Refine, Variational, GrowVariational, Pipeline, →, AutoDiff # results export Solution, ConvergenceReport, StageResult, converged, energies -export wavefunction, Wavefunction +export wavefunction, Wavefunction, convergence, radial_profile # power-user layer export build_hamiltonian_matrix, build_overlap_matrix, solve_generalized_eigenproblem, Λ, jacobi_transform, default_scale @@ -37,8 +40,10 @@ include("solution.jl") include("state.jl") include("solve.jl") include("gradient.jl") -include("observables.jl") -include("recipes.jl") +include("utils/wavefunction.jl") +include("utils/convergence.jl") +include("utils/observables.jl") +include("utils/plotting.jl") end diff --git a/src/gradient.jl b/src/gradient.jl index 8e32b51..41de911 100644 --- a/src/gradient.jl +++ b/src/gradient.jl @@ -31,23 +31,24 @@ function _encode_basis(basis::BasisSet{<:Rank0Gaussian}) for g in basis.functions C = cholesky(Symmetric(Matrix(g.A))) append!(params, _chol_to_params(Matrix(C.L))) - append!(params, Float64.(g.s)) # shift vector (unconstrained) + append!(params, Float64.(vec(parent(g.s)))) # N×3 shift, column-major end return params end # Decode a flat parameter vector back into a BasisSet{Rank0Gaussian}. -# Layout per Gaussian: [n_chol Cholesky params | n_dim shift params]. +# Layout per Gaussian: [n_chol Cholesky params | 3·n_dim shift params (N×3)]. function _decode_basis(θ::AbstractVector, n_basis::Int, n_dim::Int) T = eltype(θ) n_chol = n_dim * (n_dim + 1) ÷ 2 - n_per = n_chol + n_dim - fns = Vector{Rank0Gaussian{T, Matrix{T}, Vector{T}}}(undef, n_basis) + n_shift = 3 * n_dim + n_per = n_chol + n_shift + fns = Vector{Rank0Gaussian{T, Matrix{T}, Matrix{T}}}(undef, n_basis) for i in 1:n_basis start = (i - 1) * n_per + 1 A = _params_to_matrix(θ[start:(start + n_chol - 1)], n_dim) - s = θ[(start + n_chol):(start + n_per - 1)] - fns[i] = Rank0Gaussian(Matrix(A), Vector(s)) + s = reshape(θ[(start + n_chol):(start + n_per - 1)], n_dim, 3) + fns[i] = Rank0Gaussian(Matrix(A), Matrix(s)) end return BasisSet(fns) end @@ -57,11 +58,11 @@ end function _variational_engine( terms, n::Int, θ0, scale::Float64, maxiter::Int, gtol::Float64, verbose::Bool; - shift_init::Symbol = :qmc + shift_init::Symbol = :zeros ) n_dim = size(first(op for op in terms if op isa KineticOperator).K, 1) n_chol = n_dim * (n_dim + 1) ÷ 2 # Cholesky params per Gaussian - n_per = n_chol + n_dim # total params per Gaussian (A + shift) + n_per = n_chol + 3 * n_dim # total params per Gaussian (A + N×3 shift) regularization = 1.0e-10 if θ0 === nothing @@ -70,7 +71,7 @@ function _variational_engine( for i in 1:n bij = generate_bij(:quasirandom, i, length(w_list), scale) A = _generate_A_matrix(bij, w_list) - s = shift_init === :zeros ? zeros(n_dim) : generate_shift(:quasirandom, i, n_dim, scale) + s = shift_init === :zeros ? zeros(n_dim, 3) : generate_shift(:quasirandom, i, n_dim, scale) push!(fns, Rank0Gaussian(A, s)) end θ0 = _encode_basis(BasisSet(fns)) @@ -116,6 +117,9 @@ function _variational_engine( catch zeros(Float64, length(θ)) end + # ForwardDiff through a near-singular overlap can yield silent NaNs; + # treat a non-finite gradient as a zero-gradient barrier. + all(isfinite, G) || (G = zeros(Float64, length(θ))) return val, G end @@ -140,11 +144,11 @@ end function _sequential_engine( terms, n::Int, θ0, scale::Float64, candidates::Int, maxiter_step::Int, gtol::Float64, verbose::Bool; - shift_init::Symbol = :qmc + shift_init::Symbol = :zeros ) n_dim = size(first(op for op in terms if op isa KineticOperator).K, 1) n_chol = n_dim * (n_dim + 1) ÷ 2 - n_per = n_chol + n_dim + n_per = n_chol + 3 * n_dim w_list = [op.w for op in terms if op isa CoulombOperator] regularization = 1.0e-10 @@ -169,7 +173,7 @@ function _sequential_engine( attempt = (step - 1) * candidates + c bij = generate_bij(:quasirandom, attempt, length(w_list), scale) A = _generate_A_matrix(bij, w_list) - s = shift_init === :zeros ? zeros(n_dim) : generate_shift(:quasirandom, attempt, n_dim, scale) + s = shift_init === :zeros ? zeros(n_dim, 3) : generate_shift(:quasirandom, attempt, n_dim, scale) cand = Rank0Gaussian(A, s) θ_c = _encode_basis(BasisSet([cand])) θ_t = [θ_running; θ_c] @@ -235,6 +239,7 @@ function _sequential_engine( catch zeros(Float64, length(θ)) end + all(isfinite, G) || (G = zeros(Float64, length(θ))) return val, G end diff --git a/src/linalg.jl b/src/linalg.jl index a30088e..9a66634 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -12,24 +12,26 @@ Return the ECG overlap matrix `S` with entries `` for a `BasisSet`. """ function build_overlap_matrix(basis::BasisSet{<:GaussianBase}) n = length(basis.functions) - T = eltype(parent(first(basis.functions).A)) - S = Matrix{T}(undef, n, n) + # Infer the scalar type from an actual matrix element so complex overlaps + # (e.g. spin-orbit / spinor bases) allocate a complex matrix. + v11 = _compute_overlap_element(basis.functions[1], basis.functions[1]) + S = Matrix{typeof(v11)}(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] = val - S[j, i] = val + S[j, i] = conj(val) end return S end function _build_operator_matrix(basis::BasisSet{<:GaussianBase}, op::FewBodyHamiltonians.Operator) n = length(basis.functions) - T = eltype(parent(first(basis.functions).A)) - H = Matrix{T}(undef, n, n) + v11 = _compute_matrix_element(basis.functions[1], basis.functions[1], op) + H = Matrix{typeof(v11)}(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] = val - H[j, i] = val + H[j, i] = conj(val) end return H end @@ -43,7 +45,10 @@ terms. """ function build_hamiltonian_matrix(basis::BasisSet{<:GaussianBase}, operators::AbstractVector{<:FewBodyHamiltonians.Operator}) n = length(basis.functions) - T = eltype(parent(first(basis.functions).A)) + g1 = first(basis.functions) + # Infer the scalar type from matrix elements (spin wrappers expose no `A`). + T = isempty(operators) ? typeof(_compute_overlap_element(g1, g1)) : + mapreduce(op -> typeof(_compute_matrix_element(g1, g1, op)), promote_type, operators) H = zeros(T, n, n) for op in operators H .+= _build_operator_matrix(basis, op) @@ -58,8 +63,8 @@ Solve the symmetric generalized eigenproblem `H*c = E*S*c`, returning eigenvalues and `S`-orthonormal eigenvectors. """ function solve_generalized_eigenproblem( - H::AbstractMatrix{<:Real}, - S::AbstractMatrix{<:Real}; + H::AbstractMatrix{<:Number}, + S::AbstractMatrix{<:Number}; max_condition::Real = 1.0e12, regularization::Real = 0.0 ) @@ -71,8 +76,11 @@ function solve_generalized_eigenproblem( error("Overlap matrix S contains NaN or Inf values") end - H_sym = Symmetric((H + H') / 2) - S_sym = Symmetric((S + S') / 2) + # Hermitian symmetrisation works for both real (→ Symmetric) and complex + # (→ conjugate-symmetric) input; eigenvalues are real, eigenvectors keep + # their (possibly complex) scalar type. + H_sym = Hermitian((H + H') / 2) + S_sym = Hermitian((S + S') / 2) cond_S = cond(S_sym) if cond_S > max_condition @@ -82,28 +90,29 @@ function solve_generalized_eigenproblem( end if regularization > 0 - S_sym = Symmetric(Matrix(S_sym) + regularization * I) + S_sym = Hermitian(Matrix(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 = Symmetric(Matrix(S_sym) + ε * I) + S_sym = Hermitian(Matrix(S_sym) + ε * I) if !isposdef(S_sym) error("Overlap matrix not positive definite even after regularization") end end - # Solve the generalised symmetric eigenvalue problem H c = λ S c via - # LAPACK's divide-and-conquer driver (dsygvd). This is more reliable than - # manually factorising S and back-transforming, and returns eigenvectors - # normalised so that vᵀ S v = I. + # Solve the generalised Hermitian eigenvalue problem H c = λ S c via + # LAPACK's divide-and-conquer driver (dsygvd / zhegvd). This is more + # reliable than manually factorising S and back-transforming, and returns + # eigenvectors normalised so that cᴴ S c = I. Eigenvalues are real; + # eigenvectors are kept complex when the problem is complex. local evals, vecs try F = eigen(H_sym, S_sym) evals = real.(F.values) - vecs = real.(F.vectors) + vecs = F.vectors catch e @error "Generalised eigenvalue decomposition failed" exception = e rethrow(e) diff --git a/src/matrix_elements.jl b/src/matrix_elements.jl index 80f569c..1bced05 100644 --- a/src/matrix_elements.jl +++ b/src/matrix_elements.jl @@ -6,6 +6,10 @@ compute_matrix_element(bra, ket, op) Compute the matrix element ⟨bra|op|ket⟩ using analytic expressions. """ +# Supervector contraction tr(xᵀ M y): reduces to xᵀ M y for length-N vectors, +# and correctly sums the three Cartesian channels for N×3 shifts. +_superdot(x, M, y) = tr(transpose(x) * M * y) + function _compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian) A, B = parent(bra.A), parent(ket.A) a, b = bra.s, ket.s @@ -13,7 +17,8 @@ function _compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian) R = inv(S) n = size(S, 1) M0 = (π^n / det(S))^(3 / 2) - return exp(0.25 * (a + b)' * R * (a + b)) * M0 + v = a + b + return exp(0.25 * _superdot(v, R, v)) * M0 end function _compute_matrix_element(bra::Rank1Gaussian, ket::Rank1Gaussian) @@ -42,9 +47,11 @@ function _compute_matrix_element(bra::Rank2Gaussian, ket::Rank2Gaussian) _check_polarization_compat(ket.a, ket.b) _check_polarization_compat(bra.a, ket.a) if (any(!iszero, bra.s) || any(!iszero, ket.s)) && _pol_ncomp(bra.a) > 1 - throw(ArgumentError( - "Rank2 overlap with nonzero shifts currently requires single-component polarizations" - )) + throw( + ArgumentError( + "Rank2 overlap with nonzero shifts currently requires single-component polarizations" + ) + ) end A, B = bra.A, ket.A @@ -94,11 +101,12 @@ function _compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian, op::Kin S = A + B R = inv(S) M = _compute_matrix_element(bra, ket) + v = a + b 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) + _superdot(b, K, a) + + _superdot(v, R * B * K * A * R, v) - + _superdot(v, R * B * K, a) - + _superdot(b, K * A * R, v) return term * M end @@ -120,7 +128,7 @@ function _compute_matrix_element(bra::Rank1Gaussian, ket::Rank1Gaussian, op::Kin T2 = _polar_contract(a, K, b) * M0 T3 = ( _polar_contract(b, R * A * K * B * R, a) + - _polar_contract(a, R * A * K * B * R, b) + _polar_contract(a, R * A * K * B * R, b) ) * M0 T4 = -_polar_contract(a, R * A * K, b) * M0 T5 = -_polar_contract(b, R * B * K, a) * M0 @@ -235,8 +243,12 @@ function _compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian, op::Cou S = A + B R = inv(S) M = _compute_matrix_element(bra, ket) + v = a + b β = 1 / (w' * R * w) - q = 0.5 * (w' * R * (a + b)) + # With 3D shifts the mean displacement wᵀR v is a Cartesian 3-vector; the + # shifted Coulomb kernel depends on its magnitude q = ‖wᵀR v‖/2. + qvec = 0.5 * (transpose(w) * R * v) + q = norm(qvec) # Use the limiting form 2√(β/π) when the scaled argument x = √β·q is # small so that erf(x)/q = √β·erf(x)/x → 2√(β/π) accurately. # Thresholding on x (not q alone) correctly handles all β values. @@ -260,8 +272,8 @@ function _compute_matrix_element(bra::Rank1Gaussian, ket::Rank1Gaussian, op::Cou return op.coefficient * ( 2 * sqrt(β / π) * M1 - - sqrt(β / π) * β / 3 * - _polar_project_dot(a, Rw, b, Rw) * M0 + sqrt(β / π) * β / 3 * + _polar_project_dot(a, Rw, b, Rw) * M0 ) end @@ -325,13 +337,112 @@ function _compute_matrix_element(bra::Rank2Gaussian, ket::Rank2Gaussian, op::Cou return op.coefficient * (term1 + term2 + term3) end +# Overlap prefactor (π^n / det B)^(3/2) for an n-coordinate, 3D ECG. +_overlap_prefactor(B) = (π^size(B, 1) / det(B))^(3 / 2) + +# Central rank-0 datum after shifting the combined exponent by `W`: the inverse +# exponent `R`, the combined shift `v = bra.s + ket.s`, and the scalar overlap +# `M` of the W-shifted Gaussians. Central operators reduce to `coefficient·M`. +function _updated_rank0_data(bra::Rank0Gaussian, ket::Rank0Gaussian, W) + Bp = parent(bra.A) + parent(ket.A) + W + R = inv(Bp) + v = bra.s + ket.s + M = exp(_superdot(v, R, v) / 4) * _overlap_prefactor(Bp) + return (; R, v, M) +end + function _compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian, op::GaussianOperator) - A, B = parent(bra.A), parent(ket.A) - a, b = bra.s, ket.s - γ, w = op.γ, op.w - n = size(A, 1) - # V(r_ij) = exp(-γ (w'r)²) shifts the exponent matrix S → S' = S + γ ww' - S_prime = Symmetric(A + B + γ * (w * w')) - R_prime = inv(S_prime) - return op.coefficient * exp(0.25 * (a + b)' * R_prime * (a + b)) * (π^n / det(S_prime))^(3 / 2) + # V(r_ij) = coefficient·exp(-γ (w'r)²) shifts the exponent by γ w wᵀ. + data = _updated_rank0_data(bra, ket, op.γ * (op.w * op.w')) + return op.coefficient * data.M +end + +function _compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian, op::ManyBodyGaussianOperator) + # V = coefficient·exp(-rᵀ W r) shifts the exponent by the full matrix W. + data = _updated_rank0_data(bra, ket, op.W) + return op.coefficient * data.M +end + +function _compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian, op::OscillatorOperator) + # V = coefficient·|wᵀr|² is a second radial moment of the (unshifted-exponent) + # overlap: ⟨(wᵀr)²⟩ = 3·(wᵀ R w)/2 (3D variance) + ‖wᵀ R v / 2‖² (mean²). + w = op.w + data = _updated_rank0_data(bra, ket, zero(parent(bra.A))) + R, v, M = data.R, data.v, data.M + mean = w' * R * v / 2 + moment = 3 * (w' * R * w) / 2 + dot(vec(mean), vec(mean)) + return op.coefficient * moment * M +end + +# ── Spin-½ local matrix elements ───────────────────────────────────────────── +# ⟨br|Sᶜ|kt⟩ for a single site; c ∈ (1,2,3) = (x,y,z), in units of ℏ. +function _spin_element(c::Int, br::SpinProjection, kt::SpinProjection) + if c == 3 + return br == kt ? (br == up ? 0.5 + 0.0im : -0.5 + 0.0im) : 0.0 + 0.0im + elseif c == 1 + return br == kt ? (0.0 + 0.0im) : (0.5 + 0.0im) + else + return (br == up && kt == down) ? -0.5im : + (br == down && kt == up) ? 0.5im : 0.0 + 0.0im + end +end + +_spin_vec(br::SpinProjection, kt::SpinProjection) = + ComplexF64[_spin_element(c, br, kt) for c in 1:3] + +# Product of spectator-site δ's; sites in `exclude` are handled by spin operators. +function _spectator(bra::SpinState, ket::SpinState, exclude) + length(bra.projections) == length(ket.projections) || return 0.0 + for k in eachindex(bra.projections) + k in exclude && continue + bra.projections[k] == ket.projections[k] || return 0.0 + end + return 1.0 +end + +_spin_overlap(bra::SpinState, ket::SpinState) = _spectator(bra, ket, ()) + +# Non-conjugating contraction Σₖ x[k] y[k]. +_cdot(x, y) = sum(x[k] * y[k] for k in eachindex(x)) + +# Spin overlap and central operators factor through the spin overlap. +_compute_matrix_element(bra::SpinGaussian, ket::SpinGaussian) = + _compute_matrix_element(bra.orbital, ket.orbital) * _spin_overlap(bra.spin, ket.spin) + +_compute_matrix_element(bra::SpinGaussian, ket::SpinGaussian, op::FewBodyHamiltonians.Operator) = + _compute_matrix_element(bra.orbital, ket.orbital, op) * _spin_overlap(bra.spin, ket.spin) + +function _compute_matrix_element(bra::SpinGaussian, ket::SpinGaussian, op::GaussianTensorOperator) + i, j, w = op.i, op.j, op.w + spectator = _spectator(bra.spin, ket.spin, (i, j)) + spectator == 0 && return 0.0 + 0.0im + data = _updated_rank0_data(bra.orbital, ket.orbital, op.γ * (op.w * op.w')) + R, v, M = data.R, data.v, data.M + q = vec(w' * R * v / 2) # 3-vector mean displacement + var = (w' * R * w) / 2 # per-Cartesian variance + Si = _spin_vec(bra.spin.projections[i], ket.spin.projections[i]) + Sj = _spin_vec(bra.spin.projections[j], ket.spin.projections[j]) + # Σ_ab ⟨rₐr_b⟩ Siₐ Sj_b = var (Si·Sj) + (Si·q)(Sj·q) + spatial = var * _cdot(Si, Sj) + _cdot(Si, q) * _cdot(Sj, q) + result = spatial + if op.traceless + trace = 3 * var + dot(q, q) # ⟨r²⟩ + result -= trace / 3 * _cdot(Si, Sj) + end + return op.coefficient * M * spectator * result +end + +function _compute_matrix_element(bra::SpinGaussian, ket::SpinGaussian, op::GaussianSpinOrbitOperator) + i, j, w = op.i, op.j, op.w + data = _updated_rank0_data(bra.orbital, ket.orbital, op.γ * (op.w * op.w')) + R, M = data.R, data.M + left = vec(w' * R * parent(bra.orbital.s)) # 3-vector from bra shift + right = vec(w' * R * parent(ket.orbital.s)) # 3-vector from ket shift + Lvec = -im .* cross(left, right) .* (M / 4) # complex orbital 3-vector + spec_i = _spectator(bra.spin, ket.spin, (i,)) + spec_j = _spectator(bra.spin, ket.spin, (j,)) + Si = _spin_vec(bra.spin.projections[i], ket.spin.projections[i]) + Sj = _spin_vec(bra.spin.projections[j], ket.spin.projections[j]) + Sfull = spec_i .* Si .+ spec_j .* Sj # ⟨Sᵢ + Sⱼ⟩ (3-vector) + return op.coefficient * _cdot(Lvec, Sfull) end diff --git a/src/operators.jl b/src/operators.jl index b438d17..ded6989 100644 --- a/src/operators.jl +++ b/src/operators.jl @@ -137,12 +137,12 @@ end function Base.:+(ops::Operators, term::Tuple{<:AbstractString, <:Integer, <:Integer, <:Real}) name, i, j, coeff = term - name == "Coulomb" || - throw(ArgumentError("Unknown operator \"$name\". Supported: \"Coulomb\".")) + name in ("Coulomb", "Oscillator") || + throw(ArgumentError("Unknown operator \"$name\". Supported: \"Coulomb\", \"Oscillator\".")) ops.masses !== nothing || throw( ArgumentError( - "String-based \"Coulomb\" requires Operators(masses)." + "String-based \"$name\" requires Operators(masses)." ) ) i != j || throw(ArgumentError("Particle indices must be distinct, got i = j = $i.")) @@ -153,7 +153,11 @@ function Base.:+(ops::Operators, term::Tuple{<:AbstractString, <:Integer, <:Inte e_ij[i] = 1.0 e_ij[j] = -1.0 w = ops._U' * e_ij - push!(ops.terms, CoulombOperator(Float64(coeff), w)) + if name == "Coulomb" + push!(ops.terms, CoulombOperator(Float64(coeff), w)) + else + push!(ops.terms, OscillatorOperator(Float64(coeff), w)) + end return ops end @@ -180,6 +184,10 @@ function Base.show(io::IO, ops::Operators) println(io, " + $(op.coefficient) × Coulomb(w = $(round.(op.w; digits = 3)))") elseif op isa GaussianOperator println(io, " + $(op.coefficient) × Gaussian(γ = $(round(op.γ; digits = 3)), w = $(round.(op.w; digits = 3)))") + elseif op isa OscillatorOperator + println(io, " + $(op.coefficient) × Oscillator(w = $(round.(op.w; digits = 3)))") + elseif op isa ManyBodyGaussianOperator + println(io, " + $(op.coefficient) × ManyBodyGaussian(W = $(round.(op.W; digits = 3)))") else println(io, " + $(typeof(op))") end diff --git a/src/sampling.jl b/src/sampling.jl index 0341393..bd5da3f 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -17,10 +17,12 @@ function generate_bij( end 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) : + # Sample a full three-dimensional shift for each of the `dim` Jacobi + # coordinates, returned as the paper's `dim × 3` supervector. + u = method === :quasirandom ? _qmc_point(i, 3 * dim; sampler = qmc_sampler) : + method === :random ? rand(3 * dim) : error("Unsupported method $method") - return scale .* (2.0 .* u .- 1.0) + return reshape(scale .* (2.0 .* u .- 1.0), dim, 3) end function _generate_A_matrix(bij::AbstractVector{<:Real}, w_list::AbstractVector{<:AbstractVector{<:Real}}) diff --git a/src/solution.jl b/src/solution.jl index 51c8219..e540889 100644 --- a/src/solution.jl +++ b/src/solution.jl @@ -51,7 +51,7 @@ ascending), `basis::BasisSet`, `coefficients` (generalized eigenvectors, struct Solution E::Vector{Float64} basis::BasisSet - coefficients::Matrix{Float64} + coefficients::AbstractMatrix{<:Number} operators::Vector{FewBodyHamiltonians.Operator} state::Int stages::Vector{StageResult} diff --git a/src/solve.jl b/src/solve.jl index 09e9e92..d4d242f 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -14,11 +14,11 @@ function _ctx(terms, masses; state, tol, window, verbose) state ≥ 1 || throw(ArgumentError("state must be ≥ 1, got $state")) w_list = Vector{Float64}[ op.w for op in terms - if op isa Union{CoulombOperator, GaussianOperator} + if op isa Union{CoulombOperator, GaussianOperator, OscillatorOperator} ] isempty(w_list) && throw( ArgumentError( - "need at least one pairwise potential term (Coulomb/Gaussian) " * + "need at least one pairwise potential term (Coulomb/Gaussian/Oscillator) " * "to define the candidate geometry" ) ) diff --git a/src/state.jl b/src/state.jl index 80f365d..250eab2 100644 --- a/src/state.jl +++ b/src/state.jl @@ -36,12 +36,15 @@ function _candidate_columns(cand, basis, operators) end # Draw the next quasi-random Rank0 candidate; advances the stream counter. +# Stochastic candidates are unshifted (s = 0): the standard SVM correlated +# Gaussian for spatially symmetric (L = 0) ground states, where a nonzero shift +# only breaks the symmetry and degrades convergence. Shifted bases are reached +# through the gradient methods and manual construction (both N×3-aware). function _draw_candidate!(st::BasisState, scale::Float64, sampler, w_list, d) st.draw += 1 bij = generate_bij(:quasirandom, st.draw, length(w_list), scale; qmc_sampler = sampler) A = _generate_A_matrix(bij, w_list) - s = generate_shift(:quasirandom, st.draw, d, scale; qmc_sampler = sampler) - return Rank0Gaussian(A, s) + return Rank0Gaussian(A, zeros(d, 3)) end # Append `cand` (whose columns are `cols`): update eigensolver + S/H caches. diff --git a/src/types.jl b/src/types.jl index 9618556..3c421ef 100644 --- a/src/types.jl +++ b/src/types.jl @@ -10,24 +10,53 @@ Concrete subtypes differ by the rank of the polynomial prefactor: """ abstract type GaussianBase end +""" + _shift_matrix(s::AbstractVector) -> Matrix + +Map a legacy length-`N` shift onto the paper's three-dimensional `N × 3` +supervector by placing it in the `z` Cartesian component (third column). +A matrix argument is returned unchanged. +""" +function _shift_matrix(s::AbstractVector{T}) where {T <: Real} + result = zeros(T, length(s), 3) + result[:, 3] .= s + return result +end +_shift_matrix(s::AbstractMatrix{<:Real}) = s + +""" + _gaussian_data(A, s) -> (A::Matrix, s::Matrix) + +Validate and normalise the exponent matrix `A` (square) and the `N × 3` shift +`s`, promoting both to a common element type. +""" +function _gaussian_data(A::AbstractMatrix{<:Real}, s::AbstractMatrix{<:Real}) + size(A, 1) == size(A, 2) || throw(ArgumentError("A must be square")) + size(s) == (size(A, 1), 3) || + throw(ArgumentError("s must have size (size(A,1), 3)")) + T = promote_type(eltype(A), eltype(s)) + return Matrix{T}(A), Matrix{T}(s) +end + """ Rank0Gaussian(A, s) -Basis function ``g(\\mathbf{r}) = \\exp(-\\mathbf{r}^T A\\,\\mathbf{r} + \\mathbf{s}^T\\mathbf{r})``. +Basis function ``g(\\mathbf{r}) = \\exp(-\\mathbf{r}^T A\\,\\mathbf{r} + \\operatorname{tr}(s^T \\mathbf{r}))``. # Fields - `A` : symmetric positive-definite ``n_{\\text{dim}} \\times n_{\\text{dim}}`` matrix controlling the Gaussian width and correlations. -- `s` : shift vector ``\\mathbf{s} \\in \\mathbb{R}^{n_{\\text{dim}}}``; controls the location of the Gaussian maximum. +- `s` : shift supervector of size ``n_{\\text{dim}} \\times 3``; row `i` is the Cartesian shift of Jacobi coordinate `i`. A length-`N` vector is accepted for compatibility and mapped to the `z` component. """ -struct Rank0Gaussian{T <: Real, M <: AbstractMatrix{T}, V <: AbstractVector{T}} <: GaussianBase +struct Rank0Gaussian{T <: Real, M <: AbstractMatrix{T}, S <: AbstractMatrix{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)")) - return new{T, typeof(A), typeof(s)}(Symmetric(A), s) + s::S + function Rank0Gaussian(A::AbstractMatrix{<:Real}, s::AbstractMatrix{<:Real}) + Ad, sd = _gaussian_data(A, s) + return new{eltype(Ad), typeof(Ad), typeof(sd)}(Symmetric(Ad), sd) end end +Rank0Gaussian(A::AbstractMatrix{<:Real}, s::AbstractVector{<:Real}) = + Rank0Gaussian(A, _shift_matrix(s)) const Polarization{T} = Union{AbstractVector{T}, AbstractMatrix{T}} @@ -84,11 +113,11 @@ Rank-1 (p-wave-like) ECG basis function with linear prefactor. - a matrix of size `size(A,1) × ncomp` (multi-component polarization). """ struct Rank1Gaussian{ - T <: Real, - M <: AbstractMatrix{T}, - P <: Polarization{T}, - V <: AbstractVector{T}, -} <: GaussianBase + T <: Real, + M <: AbstractMatrix{T}, + P <: Polarization{T}, + V <: AbstractVector{T}, + } <: GaussianBase A::Symmetric{T, M} a::P s::V @@ -112,12 +141,12 @@ match `size(A,1)`. For matrix polarizations, `a` and `b` must have the same number of columns (`ncomp`), enabling multi-component pure d-wave channels. """ struct Rank2Gaussian{ - T <: Real, - M <: AbstractMatrix{T}, - P <: Polarization{T}, - Q <: Polarization{T}, - V <: AbstractVector{T}, -} <: GaussianBase + T <: Real, + M <: AbstractMatrix{T}, + P <: Polarization{T}, + Q <: Polarization{T}, + V <: AbstractVector{T}, + } <: GaussianBase A::Symmetric{T, M} a::P b::Q @@ -204,6 +233,141 @@ struct GaussianOperator{T <: Real} <: FewBodyHamiltonians.PotentialTerm w::AbstractVector{T} end +""" + OscillatorOperator(coefficient, w) + +Harmonic (oscillator) two-body potential ``V = \\text{coefficient}\\cdot|w^T\\mathbf{r}|^2``, +where ``w^T\\mathbf{r}`` is the inter-particle coordinate selected by the Jacobi +weight vector `w`. + +# Fields +- `coefficient` : coupling constant. +- `w` : weight vector in Jacobi coordinates selecting the pair. +""" +struct OscillatorOperator{T <: Real} <: FewBodyHamiltonians.PotentialTerm + coefficient::T + w::AbstractVector{T} + function OscillatorOperator(coefficient::Real, w::AbstractVector{<:Real}) + T = promote_type(typeof(coefficient), eltype(w)) + return new{T}(T(coefficient), Vector{T}(w)) + end +end + +""" + ManyBodyGaussianOperator(coefficient, W) + +Many-body Gaussian interaction ``V = \\text{coefficient}\\cdot\\exp(-\\mathbf{r}^T W\\,\\mathbf{r})`` +with a symmetric positive-definite exponent matrix `W` acting on all Jacobi +coordinates at once (e.g. a repulsive regulator). + +# Fields +- `coefficient` : coupling constant. +- `W` : symmetric positive-definite exponent matrix. +""" +struct ManyBodyGaussianOperator{T <: Real} <: FewBodyHamiltonians.PotentialTerm + coefficient::T + W::AbstractMatrix{T} + function ManyBodyGaussianOperator(coefficient::Real, W::AbstractMatrix{<:Real}) + issymmetric(W) || throw(ArgumentError("W must be symmetric")) + isposdef(W) || throw(ArgumentError("W must be positive definite")) + T = promote_type(typeof(coefficient), eltype(W)) + return new{T}(T(coefficient), Matrix{T}(W)) + end +end + +""" + SpinProjection + +Spin-½ projection eigenstates: [`up`](@ref) (+½) and [`down`](@ref) (−½). +""" +@enum SpinProjection down = -1 up = 1 + +""" + up :: SpinProjection + +Spin-½ projection eigenstate with eigenvalue +½. +""" +up + +""" + down :: SpinProjection + +Spin-½ projection eigenstate with eigenvalue −½. +""" +down + +""" + SpinState(projections) + +Direct-product spin-½ state: one [`SpinProjection`](@ref) per particle site. +""" +struct SpinState + projections::Vector{SpinProjection} +end +SpinState(projections::AbstractVector{SpinProjection}) = SpinState(collect(projections)) + +""" + SpinGaussian(orbital, spin) + +An explicitly correlated Gaussian with an attached direct-product spin state. +`orbital` is a [`Rank0Gaussian`](@ref); `spin` is a [`SpinState`](@ref). Only +introduced to support the tensor and spin-orbit interactions; central operators +factor through the spin overlap. +""" +struct SpinGaussian{G <: Rank0Gaussian} <: GaussianBase + orbital::G + spin::SpinState +end + +""" + GaussianTensorOperator(coefficient, γ, w, i, j; traceless = true) + +Gaussian-form tensor interaction coupling the coordinate `wᵀr` (range `γ`) to +the spins on sites `i` and `j`. With `traceless = true` the rank-2 spatial +tensor `rₐr_b − ⅓r²δₐ_b` is used. +""" +struct GaussianTensorOperator{T <: Real} <: FewBodyHamiltonians.PotentialTerm + coefficient::T + γ::T + w::AbstractVector{T} + i::Int + j::Int + traceless::Bool + function GaussianTensorOperator( + coefficient::Real, γ::Real, w::AbstractVector{<:Real}, + i::Integer, j::Integer; traceless::Bool = true + ) + γ > 0 || throw(ArgumentError("γ must be positive, got γ = $γ")) + i != j || throw(ArgumentError("sites i and j must differ, got i = j = $i")) + T = promote_type(typeof(coefficient), typeof(γ), eltype(w)) + return new{T}(T(coefficient), T(γ), Vector{T}(w), Int(i), Int(j), traceless) + end +end + +""" + GaussianSpinOrbitOperator(coefficient, γ, w, i, j) + +Gaussian-form spin-orbit interaction coupling the orbital motion in the +coordinate `wᵀr` (range `γ`) to the total spin `Sᵢ + Sⱼ`. Produces complex +Hermitian matrix elements for shifted Gaussians. +""" +struct GaussianSpinOrbitOperator{T <: Real} <: FewBodyHamiltonians.PotentialTerm + coefficient::T + γ::T + w::AbstractVector{T} + i::Int + j::Int + function GaussianSpinOrbitOperator( + coefficient::Real, γ::Real, w::AbstractVector{<:Real}, + i::Integer, j::Integer + ) + γ > 0 || throw(ArgumentError("γ must be positive, got γ = $γ")) + i != j || throw(ArgumentError("sites i and j must differ, got i = j = $i")) + T = promote_type(typeof(coefficient), typeof(γ), eltype(w)) + return new{T}(T(coefficient), T(γ), Vector{T}(w), Int(i), Int(j)) + end +end + struct ECG{G <: GaussianBase, O} basis::BasisSet{G} operators::Vector{O} diff --git a/src/utils/convergence.jl b/src/utils/convergence.jl new file mode 100644 index 0000000..bbd3d25 --- /dev/null +++ b/src/utils/convergence.jl @@ -0,0 +1,11 @@ +""" + convergence(sol::Solution) -> (steps, history) + +Return the cumulative solver-step indices `1:length(energies(sol))` together +with the per-step target-state energy `history = energies(sol)`, ready for +plotting a convergence curve. See also [`energies`](@ref) and `plot(sol)`. +""" +function convergence(sol::Solution) + history = energies(sol) + return 1:length(history), history +end diff --git a/src/utils/observables.jl b/src/utils/observables.jl new file mode 100644 index 0000000..00fda8a --- /dev/null +++ b/src/utils/observables.jl @@ -0,0 +1,34 @@ +""" + radial_profile(ψ::Wavefunction; coord = 1, rmax = 10.0, npoints = 400, normalize = true) + +Sample the radial density `r²|ψ(r)|²` along Jacobi coordinate `coord` on the +physical half-line `r ≥ 0` (the other coordinates held at zero), returning +`(r, density)`. When `normalize = true` the density is scaled so that its +trapezoidal integral over `[0, rmax]` equals 1. + +Because `r²|ψ|²` is defined only for non-negative radial distance, no mirrored +negative-`r` branch is produced. +""" +function radial_profile( + ψ::Wavefunction; + coord::Int = 1, rmax::Real = 10.0, npoints::Int = 400, normalize::Bool = true + ) + d = size(first(ψ.basis.functions).A, 1) + 1 ≤ coord ≤ d || throw(ArgumentError("coord must be in 1:$d")) + r = collect(range(0.0, float(rmax), length = npoints)) + density = similar(r) + v = zeros(d) + for k in eachindex(r) + fill!(v, 0.0) + v[coord] = r[k] + density[k] = r[k]^2 * abs2(ψ(v)) + end + if normalize + area = zero(eltype(density)) + for i in 1:(npoints - 1) + area += (density[i] + density[i + 1]) * (r[i + 1] - r[i]) / 2 + end + area > 0 && (density ./= area) + end + return r, density +end diff --git a/src/recipes.jl b/src/utils/plotting.jl similarity index 71% rename from src/recipes.jl rename to src/utils/plotting.jl index 819dcdf..07c2d21 100644 --- a/src/recipes.jl +++ b/src/utils/plotting.jl @@ -27,20 +27,13 @@ using RecipesBase end end -# plot(ψ; coord = 1, rmax = 10.0, npoints = 400): radial profile r²|ψ|² -# along one Jacobi coordinate (others fixed at 0). +# plot(ψ; coord = 1, rmax = 10.0, npoints = 400): half-line radial density +# r²|ψ|² along one Jacobi coordinate, delegated to `radial_profile`. @recipe function f(ψ::Wavefunction; coord = 1, rmax = 10.0, npoints = 400) - d = length(first(ψ.basis.functions).s) - 1 ≤ coord ≤ d || throw(ArgumentError("coord must be in 1:$d")) - rs = range(1.0e-3, rmax, length = npoints) - ys = map(rs) do r - v = zeros(d) - v[coord] = r - r^2 * abs2(ψ(v)) - end + r, density = radial_profile(ψ; coord, rmax, npoints) xguide --> "r (Jacobi coordinate $coord, mass-weighted)" yguide --> "r²|ψ(r)|²" label --> "|ψ|²" linewidth --> 2 - collect(rs), ys + r, density end diff --git a/src/observables.jl b/src/utils/wavefunction.jl similarity index 58% rename from src/observables.jl rename to src/utils/wavefunction.jl index 3e96a02..b85199e 100644 --- a/src/observables.jl +++ b/src/utils/wavefunction.jl @@ -4,14 +4,20 @@ Callable variational wavefunction `ψ(r) = Σᵢ cᵢ gᵢ(r)` in **Jacobi coordinates** (mass-weighted: the package's Jacobi transform normalises each relative coordinate by √μ — see `jacobi_transform`). Obtained from -[`wavefunction`](@ref); plot with `plot(ψ; coord = i)`. +[`wavefunction`](@ref); plot with `plot(ψ; coord = i)` or sample with +[`radial_profile`](@ref). """ struct Wavefunction basis::BasisSet - c::Vector{Float64} + c::AbstractVector{<:Number} end -_gauss(g, r) = exp(-(r' * g.A * r) + g.s' * r) +# `r` holds the (1D) amplitude of each Jacobi coordinate along the z axis; with +# an isotropic A this reproduces the radial Gaussian, and only the z component +# of the N×3 shift couples to it. +_shift_z(g::Rank0Gaussian) = @view parent(g.s)[:, 3] +_shift_z(g) = g.s +_gauss(g, r) = exp(-(r' * g.A * r) + _shift_z(g)' * r) _eval(g::Rank0Gaussian, r) = _gauss(g, r) _eval(g::Rank1Gaussian, r) = sum(_polar_projection(g.a, r)) * _gauss(g, r) _eval(g::Rank2Gaussian, r) = @@ -22,6 +28,9 @@ _eval(g::Rank2Gaussian, r) = """ wavefunction(sol::Solution; state = sol.state) -> Wavefunction + +Build the callable [`Wavefunction`](@ref) for the given `state` from a +[`Solution`](@ref)'s basis and generalized-eigenvector coefficients. """ wavefunction(sol::Solution; state::Int = sol.state) = Wavefunction(getfield(sol, :basis), getfield(sol, :coefficients)[:, state]) diff --git a/test/runtests.jl b/test/runtests.jl index 271d32a..96c85f5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,6 +19,9 @@ using FewBodyECG include("test_hamiltonian.jl") include("test_svm_eigen.jl") include("test_hydrogen.jl") + include("test_harmonic_oscillator.jl") + include("test_harmonium.jl") + include("test_spin_orbit.jl") include("test_utils.jl") include("test_types.jl") include("test_variational.jl") diff --git a/test/test_hamiltonian.jl b/test/test_hamiltonian.jl index ee54bfc..c23edb5 100644 --- a/test/test_hamiltonian.jl +++ b/test/test_hamiltonian.jl @@ -81,6 +81,15 @@ end @testset "solve_generalized_eigenproblem" begin + @testset "Complex Hermitian problem keeps complex eigenvectors" begin + # off-diagonal √2·i makes the eigenvalues exactly 0 and 3 + H = ComplexF64[1 sqrt(2) * im; -sqrt(2) * im 2] + S = Matrix{Float64}(I, 2, 2) + E, C = solve_generalized_eigenproblem(H, S) + @test E ≈ [0.0, 3.0] + @test eltype(C) <: Complex + end + @testset "Basic solve" begin H = [2.0 0.5; 0.5 3.0] S = [1.0 0.1; 0.1 1.0] @@ -281,6 +290,15 @@ end @testset "solve_generalized_eigenproblem" begin + @testset "Complex Hermitian problem keeps complex eigenvectors" begin + # off-diagonal √2·i makes the eigenvalues exactly 0 and 3 + H = ComplexF64[1 sqrt(2) * im; -sqrt(2) * im 2] + S = Matrix{Float64}(I, 2, 2) + E, C = solve_generalized_eigenproblem(H, S) + @test E ≈ [0.0, 3.0] + @test eltype(C) <: Complex + end + @testset "Basic solve" begin H = [2.0 0.5; 0.5 3.0] S = [1.0 0.1; 0.1 1.0] @@ -488,7 +506,7 @@ 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 + @test size(s) == (dim, 3) end end end diff --git a/test/test_harmonic_oscillator.jl b/test/test_harmonic_oscillator.jl new file mode 100644 index 0000000..6f0a1f1 --- /dev/null +++ b/test/test_harmonic_oscillator.jl @@ -0,0 +1,35 @@ +using Test +using LinearAlgebra +using FewBodyECG +import Antique + +# Isotropic 3D harmonic oscillator (k = m = ℏ = 1) modelled as a unit-mass +# particle bound to a fixed centre by an OscillatorOperator. Its l = 0 states +# (nᵣ = 0,1,2) coincide with Antique's odd 1D states n = 1,3,5 at E = 3/2, 7/2, +# 11/2. A tempered geometric basis of unshifted Gaussians resolves all three +# to spectroscopic accuracy; the sqrt(2) factor converts Antique's full-line +# odd-state normalization to the half-line reduced-radial normalization +# ∫₀^∞ u² dr = 1. +@testset "Harmonic oscillator vs Antique" begin + ops = Operators([1.0e15, 1.0]) + ops += "Kinetic" + ops += ("Oscillator", 1, 2, 0.5) + + w = only(op.w for op in ops.terms if op isa OscillatorOperator) + αs = exp10.(range(-1.0, 1.0, length = 16)) + basis = BasisSet([Rank0Gaussian(α .* (w * transpose(w)), zeros(1, 3)) for α in αs]) + H = build_hamiltonian_matrix(basis, ops) + S = build_overlap_matrix(basis) + ecg_energies, coeffs = solve_generalized_eigenproblem(H, S) + + HO = Antique.HarmonicOscillator(k = 1.0, m = 1.0, ℏ = 1.0) + @test ecg_energies[1:3] ≈ [Antique.E(HO; n = n) for n in (1, 3, 5)] atol = 1.0e-5 + + radius = range(1.0e-4, 6, length = 800) + ψ = Wavefunction(basis, coeffs[:, 1]) + ecg_wavefunction = [ψ([r]) for r in radius] + u_ecg = sqrt(4π) .* radius .* ecg_wavefunction + u_ref = sqrt(2) .* [Antique.ψ(HO, r; n = 1) for r in radius] + u_ecg .*= sign(dot(u_ecg, u_ref)) + @test maximum(abs.(u_ecg - u_ref)) < 1.0e-3 +end diff --git a/test/test_harmonium.jl b/test/test_harmonium.jl new file mode 100644 index 0000000..0e4c4fc --- /dev/null +++ b/test/test_harmonium.jl @@ -0,0 +1,20 @@ +using Test +using FewBodyECG + +# Hooke's atom (harmonium): two electrons in a harmonic trap with Coulomb +# repulsion. For ω = 1/2 the ground-state energy is exactly 2 Ha +# (Taut, Phys. Rev. A 48, 3561 (1993)). Exercises the OscillatorOperator +# together with Coulomb repulsion on the plain stochastic solver. +@testset "Hooke's atom vs Taut (ω = 1/2)" begin + ω = 0.5 + ops = Operators([1.0e15, 1.0, 1.0]) + ops += "Kinetic" + ops += ("Oscillator", 1, 2, 0.5 * ω^2) + ops += ("Oscillator", 1, 3, 0.5 * ω^2) + ops += ("Coulomb", 2, 3, 1.0) + + sol = solve(ops, SVM(basis = 80, candidates = 40, scale = 2.0)) + + @test sol.E₀ ≈ 2.0 atol = 1.0e-3 # Taut's exact energy + @test sol.E₀ > 2.0 - 1.0e-9 # variational upper bound never violated +end diff --git a/test/test_matrix_elements.jl b/test/test_matrix_elements.jl index 905eb16..71c34b4 100644 --- a/test/test_matrix_elements.jl +++ b/test/test_matrix_elements.jl @@ -478,7 +478,7 @@ end g = Rank0Gaussian(A, s) w = [1.0] op_narrow = GaussianOperator(1.0, 10.0, w) - op_wide = GaussianOperator(1.0, 0.1, w) + op_wide = GaussianOperator(1.0, 0.1, w) @test _compute_matrix_element(g, g, op_narrow) < _compute_matrix_element(g, g, op_wide) end @@ -498,3 +498,69 @@ end @test_throws ArgumentError ops += ("Gaussian", 1, 2, -1.0, -1.0) end end + +@testset "Scalar Gaussian, oscillator, and many-body elements (3D shift)" begin + A, B = [1.2;;], [0.8;;] + bra = Rank0Gaussian(B, reshape([-0.1, 0.4, 0.2], 1, 3)) + ket = Rank0Gaussian(A, reshape([0.2, -0.1, 0.3], 1, 3)) + w, γ = [1.0], 0.7 + + _M(Bmat, v) = exp(tr(v' * inv(Bmat) * v) / 4) * (π / det(Bmat))^(3 / 2) + v = parent(bra.s) + parent(ket.s) + + # Gaussian: exponent shifted by γ w wᵀ + Bg = A + B + γ * (w * w') + expected_gaussian = 1.0 * _M(Bg, v) + @test _compute_matrix_element(bra, ket, GaussianOperator(1.0, γ, w)) ≈ expected_gaussian + + # Oscillator: second radial moment of the plain overlap + Bo = A + B + R = inv(Bo) + mean = w' * R * v / 2 + moment = 3 * (w' * R * w) / 2 + dot(vec(mean), vec(mean)) + expected_oscillator = 1.0 * moment * _M(Bo, v) + @test _compute_matrix_element(bra, ket, OscillatorOperator(1.0, w)) ≈ expected_oscillator + + # Many-body: exponent shifted by the full matrix W + W = [0.4;;] + expected_manybody = 1.0 * _M(A + B + W, v) + @test _compute_matrix_element(bra, ket, ManyBodyGaussianOperator(1.0, W)) ≈ expected_manybody + + # Central operators factor through the plain overlap when W → 0 limit checked + @test _compute_matrix_element(bra, ket, ManyBodyGaussianOperator(2.0, [1.0e-9;;])) ≈ + 2.0 * _compute_matrix_element(bra, ket) rtol = 1.0e-6 +end + +@testset "Spin tensor and spin-orbit interactions" begin + # tensor Hermiticity (zero-shift s-wave orbital) + orbital = Rank0Gaussian([1.0;;], zeros(1, 3)) + updown = SpinGaussian(orbital, SpinState([up, down])) + downup = SpinGaussian(orbital, SpinState([down, up])) + tensor = GaussianTensorOperator(1.0, 0.4, [1.0], 1, 2) + @test _compute_matrix_element(updown, downup, tensor) ≈ + conj(_compute_matrix_element(downup, updown, tensor)) + + # central operators factor through the spin overlap + K = KineticOperator([0.5;;]) + @test _compute_matrix_element(updown, updown, K) ≈ + _compute_matrix_element(orbital, orbital, K) + @test _compute_matrix_element(updown, downup, K) == 0 # orthogonal spins + + # spin-orbit: nonzero complex element with shifted, non-parallel orbitals + orb1 = Rank0Gaussian([1.0;;], reshape([0.3, 0.0, 0.1], 1, 3)) + orb2 = Rank0Gaussian([1.2;;], reshape([0.0, 0.4, 0.2], 1, 3)) + sg1 = SpinGaussian(orb1, SpinState([up, up])) + sg2 = SpinGaussian(orb2, SpinState([up, up])) + so = GaussianSpinOrbitOperator(1.0, 0.5, [1.0], 1, 2) + el = _compute_matrix_element(sg1, sg2, so) + @test el isa Complex + @test abs(el) > 0 + @test _compute_matrix_element(sg1, sg2, so) ≈ conj(_compute_matrix_element(sg2, sg1, so)) + + # assembled Hamiltonian is complex Hermitian + H = build_hamiltonian_matrix(BasisSet([sg1, sg2]), FewBodyECG.Operator[so]) + @test eltype(H) <: Complex + @test H ≈ H' + E, C = solve_generalized_eigenproblem(H, build_overlap_matrix(BasisSet([sg1, sg2]))) + @test all(isfinite, E) +end diff --git a/test/test_observables.jl b/test/test_observables.jl index d3d01f3..eda973e 100644 --- a/test/test_observables.jl +++ b/test/test_observables.jl @@ -16,7 +16,7 @@ sol = solve(ops, SVM(basis = 15, candidates = 15, scale = 1.0)) c = sol.coefficients[:, 1] fns = sol.basis.functions ref = sum( - c[i] * exp(-([0.5]' * fns[i].A * [0.5]) + fns[i].s' * [0.5]) + c[i] * exp(-([0.5]' * fns[i].A * [0.5]) + (@view parent(fns[i].s)[:, 3])' * [0.5]) for i in eachindex(fns) ) @test ψ([0.5]) ≈ ref rtol = 1.0e-12 @@ -30,6 +30,45 @@ sol = solve(ops, SVM(basis = 15, candidates = 15, scale = 1.0)) @test ψ2([0.7]) ≈ 0.7 * 0.7 * exp(-0.49) rtol = 1.0e-12 end +@testset "convergence and radial_profile utilities" begin + steps, history = convergence(sol) + @test steps == 1:length(history) + @test history == energies(sol) + + r, density = radial_profile(wavefunction(sol); rmax = 4, npoints = 200) + @test first(r) ≥ 0 + @test all(≥(0), density) + @test isapprox( + sum((density[i] + density[i + 1]) * (r[i + 1] - r[i]) / 2 for i in 1:(length(r) - 1)), + 1; atol = 1.0e-8 + ) + + # unnormalized profile is the bare r²|ψ|² + r2, d2 = radial_profile(wavefunction(sol); rmax = 4, npoints = 200, normalize = false) + @test all(≥(0), d2) + + # rank-1 and rank-2 wavefunction profiles are supported and half-line normalized + g1 = Rank1Gaussian([1.0;;], [1.0], [0.0]) + ψ1 = Wavefunction(BasisSet([g1]), [1.0]) + r1, dens1 = radial_profile(ψ1; rmax = 6, npoints = 300) + @test first(r1) ≥ 0 && all(≥(0), dens1) + @test isapprox( + sum((dens1[i] + dens1[i + 1]) * (r1[i + 1] - r1[i]) / 2 for i in 1:(length(r1) - 1)), + 1; atol = 1.0e-8 + ) + + g2 = Rank2Gaussian([1.0;;], [1.0], [1.0], [0.0]) + ψ2 = Wavefunction(BasisSet([g2]), [1.0]) + r2b, dens2 = radial_profile(ψ2; rmax = 6, npoints = 300) + @test all(≥(0), dens2) + @test isapprox( + sum((dens2[i] + dens2[i + 1]) * (r2b[i + 1] - r2b[i]) / 2 for i in 1:(length(r2b) - 1)), + 1; atol = 1.0e-8 + ) + + @test_throws ArgumentError radial_profile(wavefunction(sol); coord = 5) +end + @testset "Recipes" begin # convergence recipe plots = RecipesBase.apply_recipe(Dict{Symbol, Any}(), sol) diff --git a/test/test_operators.jl b/test/test_operators.jl index 448df0e..2b1f7d4 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -407,4 +407,20 @@ import FewBodyECG: jacobi_transform, Λ @test sol_vec.E₀ < -0.46 end end + + @testset "Oscillator, Gaussian, and many-body builders/validation" begin + ops = Operators([1.0e15, 1.0]) + ops += ("Oscillator", 1, 2, 0.5) + @test ops[1] isa OscillatorOperator + ops += ("Gaussian", 1, 2, -5.0, 1.0) + @test ops[2] isa GaussianOperator + + # many-body exponent must be symmetric positive-definite + @test_throws ArgumentError ManyBodyGaussianOperator(1.0, [1.0 1.0; 0.0 1.0]) + @test_throws ArgumentError ManyBodyGaussianOperator(1.0, [-1.0;;]) + @test ManyBodyGaussianOperator(1.0, [0.4;;]) isa ManyBodyGaussianOperator + + # unknown 4-tuple operator name is rejected + @test_throws ArgumentError (Operators([1.0e15, 1.0]) + ("Nope", 1, 2, 0.5)) + end end diff --git a/test/test_sampling.jl b/test/test_sampling.jl index 22f6ae1..2c2b0bc 100644 --- a/test/test_sampling.jl +++ b/test/test_sampling.jl @@ -26,7 +26,7 @@ import FewBodyECG: _generate_A_matrix, build_rank0, generate_bij, generate_shift @testset "Full" begin s = generate_shift(:quasirandom, 1, 3, 2.0) - @test length(s) == 3 + @test size(s) == (3, 3) # N×3 supervector shift @test all(abs.(s) .<= 2.0 .+ eps()) s_q1 = generate_shift(:quasirandom, 5, 3, 1.5) diff --git a/test/test_solve.jl b/test/test_solve.jl index 2132041..cbeadd2 100644 --- a/test/test_solve.jl +++ b/test/test_solve.jl @@ -36,10 +36,11 @@ ops = Operators([1.0e15, 1.0], [+1.0, -1.0]); ops += "Kinetic"; ops += "Coulomb" sol2 = solve(ops, SVM(basis = 25, candidates = 20, scale = 1.0); state = 2) @test sol2.state == 2 && sol2.E₀ == sol2.E[2] && sol2.E₀ > sol2.E[1] - # warm start grows an existing basis + # warm start grows an existing basis (unshifted candidates may occasionally + # be linearly dependent and skipped — the report documents this honestly) small = solve(ops, SVM(basis = 5, candidates = 10, scale = 1.0)) bigger = solve(ops, SVM(basis = 10, candidates = 10, scale = 1.0); init = small) - @test length(bigger.basis.functions) == 15 + @test 5 < length(bigger.basis.functions) <= 15 @test bigger.E₀ <= small.E₀ + 1.0e-12 # early stop: an impossible independence floor rejects every candidate, diff --git a/test/test_spin_orbit.jl b/test/test_spin_orbit.jl new file mode 100644 index 0000000..1d1db01 --- /dev/null +++ b/test/test_spin_orbit.jl @@ -0,0 +1,62 @@ +using Test +using LinearAlgebra +using FewBodyECG + +# Model spin-orbit doublet: two spin-½ particles in a p-wave-like manifold of +# shifted Gaussians bound by a central well. A GaussianSpinOrbitOperator makes +# the Hamiltonian complex Hermitian and splits the degenerate upper doublet +# linearly in the coupling strength. Exercises the complex generalized +# eigensolver end-to-end. +@testset "Spin-orbit doublet (complex Hermitian)" begin + ops = Operators([1.0, 1.0]) + ops += "Kinetic" + ops += ("Gaussian", 1, 2, -4.0, 0.4) + terms = ops.terms + w = only(op.w for op in terms if op isa GaussianOperator) + + shift(dir) = (v = zeros(1, 3); v[1, dir] = 0.8; v) + basis = BasisSet( + [ + SpinGaussian(Rank0Gaussian([0.4;;], shift(d)), SpinState([up, up])) for d in 1:3 + ] + ) + + H₀ = build_hamiltonian_matrix(basis, terms) + Hₛₒ = build_hamiltonian_matrix( + basis, FewBodyECG.Operator[GaussianSpinOrbitOperator(1.0, 0.4, w, 1, 2)] + ) + S = build_overlap_matrix(basis) + + # spin-orbit-free block is real; spin-orbit block is genuinely complex Hermitian + @test eltype(H₀) <: Real + @test eltype(Hₛₒ) <: Complex + @test Hₛₒ ≈ Hₛₒ' + @test !isapprox(Hₛₒ, real.(Hₛₒ)) # actually complex, not a real matrix in disguise + + gap(λ) = begin + E, C = solve_generalized_eigenproblem(H₀ .+ λ .* Hₛₒ, S) + @test maximum(abs, imag.(E)) < 1.0e-10 # eigenvalues are real + @test eltype(C) <: Complex # eigenvectors stay complex + Es = sort(real.(E)) + Es[3] - Es[2] + end + + # degenerate at λ = 0, splits, and the gap is linear in λ + @test gap(0.0) < 1.0e-8 + g1 = gap(1.0) + g2 = gap(2.0) + @test g1 > 1.0e-3 + @test g2 ≈ 2 * g1 rtol = 1.0e-3 # linear splitting + + # nested-basis enrichment converges the ground state monotonically + αs = exp10.(range(-0.8, 0.6, length = 8)) + sg(α, d) = SpinGaussian(Rank0Gaussian([α;;], shift(d)), SpinState([up, up])) + ground = map(1:length(αs)) do k + b = BasisSet([sg(α, d) for α in αs[1:k] for d in 1:3]) + Hk = build_hamiltonian_matrix(b, terms) .+ + build_hamiltonian_matrix(b, FewBodyECG.Operator[GaussianSpinOrbitOperator(1.0, 0.4, w, 1, 2)]) + minimum(real.(solve_generalized_eigenproblem(Hk, build_overlap_matrix(b))[1])) + end + @test all(diff(ground) .≤ 1.0e-9) # monotone non-increasing (variational) + @test ground[end] < ground[1] # actually improves +end diff --git a/test/test_types.jl b/test/test_types.jl index 5489210..90fccfe 100644 --- a/test/test_types.jl +++ b/test/test_types.jl @@ -11,7 +11,7 @@ using FewBodyHamiltonians @test isa(g, Rank0Gaussian) @test isa(g.A, Symmetric) - @test g.s == s + @test g.s == FewBodyECG._shift_matrix(s) # N×3, legacy vector → z component A_ns = rand(2, 3) @test_throws ArgumentError Rank0Gaussian(A_ns, [1.0, 2.0]) @@ -62,7 +62,7 @@ end @test isa(g, Rank0Gaussian) @test isa(g.A, Symmetric) - @test g.s == s + @test g.s == FewBodyECG._shift_matrix(s) # N×3, legacy vector → z component A_ns = rand(2, 3) @test_throws ArgumentError Rank0Gaussian(A_ns, [1.0, 2.0]) @@ -155,6 +155,14 @@ end @test_throws LinearAlgebra.PosDefException validate!(g1_indef) end +@testset "Rank0Gaussian N×3 shift representation" begin + g = Rank0Gaussian([1.0;;], reshape([0.1, -0.2, 0.3], 1, 3)) + @test size(g.s) == (1, 3) + @test Rank0Gaussian([1.0;;], [0.2]).s == reshape([0.0, 0.0, 0.2], 1, 3) + # mismatched shift shape is rejected + @test_throws ArgumentError Rank0Gaussian([1.0;;], reshape([0.1, 0.2], 1, 2)) +end + @testset "BasisSet, KineticOperator, CoulombOperator, and ECG composition" begin A = [2.0 0.0; 0.0 3.0] s1 = [1.0, 0.0] diff --git a/test/test_variational.jl b/test/test_variational.jl index 919f01f..3874733 100644 --- a/test/test_variational.jl +++ b/test/test_variational.jl @@ -50,7 +50,7 @@ end basis = BasisSet([g1, g2]) θ = _encode_basis(basis) - @test length(θ) == 4 # 2 functions × (n_chol=1 + n_dim=1) + @test length(θ) == 8 # 2 functions × (n_chol=1 + 3·n_dim=3) basis2 = _decode_basis(θ, 2, 1) @test length(basis2.functions) == 2 @@ -68,7 +68,7 @@ end basis = BasisSet([g]) θ = _encode_basis(basis) - @test length(θ) == 5 # 1 function × (n_chol=3 + n_dim=2) + @test length(θ) == 9 # 1 function × (n_chol=3 + 3·n_dim=6) basis2 = _decode_basis(θ, 1, 2) @test Matrix(basis2.functions[1].A) ≈ A rtol = 1.0e-8 @@ -78,12 +78,12 @@ end @testset "_encode_basis / _decode_basis round-trip with non-zero shifts" begin # Verify shift vectors are correctly preserved through the encode/decode cycle. A = [3.0 0.5; 0.5 2.0] - s = [0.3, -0.1] + s = [0.3 -0.2 0.5; -0.1 0.4 -0.3] # full N×3 shift g = Rank0Gaussian(A, s) basis = BasisSet([g]) θ = _encode_basis(basis) - @test length(θ) == 5 + @test length(θ) == 9 basis2 = _decode_basis(θ, 1, 2) @test Matrix(basis2.functions[1].A) ≈ A rtol = 1.0e-8 @@ -200,15 +200,15 @@ end # --------------------------------------------------------------------------- @testset "shift vectors are included in optimised parameters" begin - # _encode_basis should pack n_chol + n_dim params per Gaussian. - # For a 1-D (hydrogen) basis: n_chol=1, n_dim=1 → 2 params per function. - # The second param is the shift; _decode_basis should round-trip it. + # _encode_basis packs n_chol + 3·n_dim params per Gaussian (N×3 shift). + # For a 1-D (hydrogen) basis: n_chol=1, 3·n_dim=3 → 4 params per function. + # The shift is optimised; _decode_basis round-trips it. ops = _hydrogen_ops() sol = solve(ops, Variational(basis = 4, scale = 1.0, maxiter = 50)) - # Each basis function has a 1-D shift vector stored in s. + # Each basis function has a 1×3 shift supervector stored in s. for g in sol.basis.functions - @test length(g.s) == 1 - @test isfinite(g.s[1]) + @test size(g.s) == (1, 3) + @test all(isfinite, g.s) end end From 4afc5d994b96ac584f6c38db2c46dd9044cdee04 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Fri, 10 Jul 2026 15:13:57 +0200 Subject: [PATCH 2/5] updated hookes model --- docs/src/examples/h2plus.md | 2 +- docs/src/examples/harmonium.md | 100 +++++++++++++++++++++++++-------- examples/h2plus.jl | 2 +- examples/harmonium.jl | 88 +++++++++++++++++++++-------- test/test_harmonium.jl | 13 +++++ 5 files changed, 157 insertions(+), 48 deletions(-) diff --git a/docs/src/examples/h2plus.md b/docs/src/examples/h2plus.md index fd59ce3..94149d1 100644 --- a/docs/src/examples/h2plus.md +++ b/docs/src/examples/h2plus.md @@ -29,7 +29,7 @@ println("H2+ E0 = ", sol.E₀, " Ha (reference ", h2p_ref, ", Δ = ", sol.E₀ println("bound below H + p+ threshold? ", sol.E₀ < -0.5) plot(sol, h2p_ref) -plot(wavefunction(sol); coord = 1, rmax = 80.0) +plot(wavefunction(sol); coord = 1, rmax = 100.0) ```` --- diff --git a/docs/src/examples/harmonium.md b/docs/src/examples/harmonium.md index 90bf0a1..d95e732 100644 --- a/docs/src/examples/harmonium.md +++ b/docs/src/examples/harmonium.md @@ -2,29 +2,12 @@ EditURL = "../../../examples/harmonium.jl" ``` -# Hooke's atom (harmonium) - -Two electrons in a common harmonic trap, repelling through the Coulomb -interaction: - -```math -H = -\tfrac12\nabla_1^2 - \tfrac12\nabla_2^2 - + \tfrac12\omega^2 (r_1^2 + r_2^2) + \frac{1}{r_{12}}. -``` - -For the trap frequency ``\omega = \tfrac12`` the ground-state energy is known -in closed form (Taut, *Phys. Rev. A* **48**, 3561 (1993)): **exactly 2 Ha**. -The trap is supplied by an `OscillatorOperator` between a fixed heavy centre -and each electron (coefficient ``\tfrac12\omega^2``); the only Coulomb term is -the electron–electron repulsion. The system is spherically symmetric and real, -so the plain stochastic `SVM` solver applies. - ````@example harmonium using FewBodyECG using Plots ω = 0.5 -masses = [1.0e15, 1.0, 1.0] # heavy trap centre + two electrons +masses = [1.0e15, 1.0, 1.0] ops = Operators(masses) ops += "Kinetic" @@ -38,13 +21,84 @@ sol = solve(ops, SVM(basis = 80, candidates = 40, scale = 2.0)) println("Hooke's atom E₀ = ", sol.E₀, " Ha (Taut exact ", exact, ", Δ = ", sol.E₀ - exact, ")") println("variational upper bound respected: ", sol.E₀ ≥ exact) sol +```` + +Shared styling so the three figures read as one consistent set: a magma-family +magenta for the ECG (computed) data, a neutral slate for the reference/exact. + +````@example harmonium +FIG = ( + titlefontsize = 12, guidefontsize = 10, tickfontsize = 9, legendfontsize = 9, + framestyle = :box, grid = true, gridalpha = 0.12, size = (640, 420), dpi = 200, + left_margin = 4Plots.mm, bottom_margin = 3Plots.mm, +) +ecg_color = RGB(0.72, 0.16, 0.42) # computed (ECG) +ref_color = RGB(0.45, 0.47, 0.52) # reference / exact + +plot( + sol, exact; + palette = [ecg_color, ref_color], linewidth = 2, + title = "Convergence to the exact energy", FIG..., +) +```` -# Convergence toward the exact energy -plot(sol, exact) +Back to laboratory coordinates: `jacobi_transform` gives J mapping physical +particle positions to the (mass-weighted) Jacobi coordinates the solution uses, +so ψ(J·r) evaluates the wavefunction directly in coordinate space. + +````@example harmonium +ψ = wavefunction(sol) +J, _ = jacobi_transform(masses) +Ψ(z₁, z₂) = ψ(J * [0.0, z₁, z₂]) # electrons at z₁, z₂ on the axis; centre at 0 +```` + +Relative wavefunction χ(u): electrons at ±u/2 (centre of mass at the trap origin) +vs Taut's closed form χ(u) ∝ (1 + u/2) e^{-u²/8}. + +````@example harmonium +u = range(0, 10, length = 200) +χ_ecg = [Ψ(x / 2, -x / 2) for x in u] +χ_exact = [(1 + x / 2) * exp(-x^2 / 8) for x in u] +χ_ecg ./= maximum(abs, χ_ecg) # peak-normalise for shape +χ_exact ./= maximum(χ_exact) +χ_ecg .*= sign(sum(χ_ecg .* χ_exact)) +println("max |Δχ| (shape) = ", maximum(abs, χ_ecg .- χ_exact)) + +plot(u, χ_exact; label = "exact (Taut 1993)", color = ref_color, linewidth = 3, FIG...) +plot!( + u, χ_ecg; label = "ECG", color = ecg_color, linestyle = :dash, linewidth = 2, + xlabel = "interelectronic distance u", ylabel = "χ(u) (peak-normalised)", + title = "Relative wavefunction vs Taut's closed form", legend = :topright, +) +```` + +Full two-electron density in coordinate space. The reduced amplitude along the +diagonal z₁ = z₂ is the Coulomb hole — the electrons avoid each other. + +````@example harmonium +lim = 4.5 +zs = range(-lim, lim, length = 251) +density = [abs2(Ψ(z₁, z₂)) for z₂ in zs, z₁ in zs] +density ./= maximum(density) # peak-normalise → colorbar 0…1 -# Ground-state radial density of one electron -plot(wavefunction(sol); coord = 1, rmax = 8.0) -plot(convergence(sol)) +heatmap( + zs, zs, density; + c = :magma, clims = (0, 1), + xlims = (-lim, lim), ylims = (-lim, lim), aspect_ratio = :equal, + xlabel = "electron 1 position z₁", ylabel = "electron 2 position z₂", + colorbar_title = "\n|Ψ|² (normalised)", colorbar_titlefontsize = 9, + title = "Hooke's atom — two-electron density", + titlefontsize = 12, guidefontsize = 10, tickfontsize = 9, + grid = false, framestyle = :box, widen = false, + size = (600, 500), dpi = 200, left_margin = 3Plots.mm, +) +# dashed line marks where the electrons coincide — the trough is the Coulomb hole +plot!( + [-lim, lim], [-lim, lim]; + color = :white, alpha = 0.4, linestyle = :dash, linewidth = 1.5, + label = "z₁ = z₂ (Coulomb hole)", legend = :topleft, foreground_color_legend = nothing, + background_color_legend = RGBA(0, 0, 0, 0.35), legendfontcolor = :white, legendfontsize = 8, +) ```` --- diff --git a/examples/h2plus.jl b/examples/h2plus.jl index ac0e1b4..216cbb0 100644 --- a/examples/h2plus.jl +++ b/examples/h2plus.jl @@ -24,4 +24,4 @@ println("H2+ E0 = ", sol.E₀, " Ha (reference ", h2p_ref, ", Δ = ", sol.E₀ println("bound below H + p+ threshold? ", sol.E₀ < -0.5) plot(sol, h2p_ref) -plot(wavefunction(sol); coord = 1, rmax = 80.0) +plot(wavefunction(sol); coord = 1, rmax = 100.0) diff --git a/examples/harmonium.jl b/examples/harmonium.jl index 0f9aef8..c759d5c 100644 --- a/examples/harmonium.jl +++ b/examples/harmonium.jl @@ -1,25 +1,8 @@ -# # Hooke's atom (harmonium) -# -# Two electrons in a common harmonic trap, repelling through the Coulomb -# interaction: -# -# ```math -# H = -\tfrac12\nabla_1^2 - \tfrac12\nabla_2^2 -# + \tfrac12\omega^2 (r_1^2 + r_2^2) + \frac{1}{r_{12}}. -# ``` -# -# For the trap frequency ``\omega = \tfrac12`` the ground-state energy is known -# in closed form (Taut, *Phys. Rev. A* **48**, 3561 (1993)): **exactly 2 Ha**. -# The trap is supplied by an `OscillatorOperator` between a fixed heavy centre -# and each electron (coefficient ``\tfrac12\omega^2``); the only Coulomb term is -# the electron–electron repulsion. The system is spherically symmetric and real, -# so the plain stochastic `SVM` solver applies. - using FewBodyECG using Plots ω = 0.5 -masses = [1.0e15, 1.0, 1.0] # heavy trap centre + two electrons +masses = [1.0e15, 1.0, 1.0] ops = Operators(masses) ops += "Kinetic" @@ -34,9 +17,68 @@ println("Hooke's atom E₀ = ", sol.E₀, " Ha (Taut exact ", exact, ", Δ = " println("variational upper bound respected: ", sol.E₀ ≥ exact) sol -## Convergence toward the exact energy -plot(sol, exact) +# Shared styling so the three figures read as one consistent set: a magma-family +# magenta for the ECG (computed) data, a neutral slate for the reference/exact. +FIG = ( + titlefontsize = 12, guidefontsize = 10, tickfontsize = 9, legendfontsize = 9, + framestyle = :box, grid = true, gridalpha = 0.12, size = (640, 420), dpi = 200, + left_margin = 4Plots.mm, bottom_margin = 3Plots.mm, +) +ecg_color = RGB(0.72, 0.16, 0.42) # computed (ECG) +ref_color = RGB(0.45, 0.47, 0.52) # reference / exact + +plot( + sol, exact; + palette = [ecg_color, ref_color], linewidth = 2, + title = "Convergence to the exact energy", FIG..., +) + +# Back to laboratory coordinates: `jacobi_transform` gives J mapping physical +# particle positions to the (mass-weighted) Jacobi coordinates the solution uses, +# so ψ(J·r) evaluates the wavefunction directly in coordinate space. +ψ = wavefunction(sol) +J, _ = jacobi_transform(masses) +Ψ(z₁, z₂) = ψ(J * [0.0, z₁, z₂]) # electrons at z₁, z₂ on the axis; centre at 0 + +# Relative wavefunction χ(u): electrons at ±u/2 (centre of mass at the trap origin) +# vs Taut's closed form χ(u) ∝ (1 + u/2) e^{-u²/8}. +u = range(0, 10, length = 200) +χ_ecg = [Ψ(x / 2, -x / 2) for x in u] +χ_exact = [(1 + x / 2) * exp(-x^2 / 8) for x in u] +χ_ecg ./= maximum(abs, χ_ecg) # peak-normalise for shape +χ_exact ./= maximum(χ_exact) +χ_ecg .*= sign(sum(χ_ecg .* χ_exact)) +println("max |Δχ| (shape) = ", maximum(abs, χ_ecg .- χ_exact)) + +plot(u, χ_exact; label = "exact (Taut 1993)", color = ref_color, linewidth = 3, FIG...) +plot!( + u, χ_ecg; label = "ECG", color = ecg_color, linestyle = :dash, linewidth = 2, + xlabel = "interelectronic distance u", ylabel = "χ(u) (peak-normalised)", + title = "Relative wavefunction vs Taut's closed form", legend = :topright, +) + +# Full two-electron density in coordinate space. The reduced amplitude along the +# diagonal z₁ = z₂ is the Coulomb hole — the electrons avoid each other. +lim = 4.5 +zs = range(-lim, lim, length = 251) +density = [abs2(Ψ(z₁, z₂)) for z₂ in zs, z₁ in zs] +density ./= maximum(density) # peak-normalise → colorbar 0…1 -## Ground-state radial density of one electron -plot(wavefunction(sol); coord = 1, rmax = 8.0) -plot(convergence(sol)) \ No newline at end of file +heatmap( + zs, zs, density; + c = :magma, clims = (0, 1), + xlims = (-lim, lim), ylims = (-lim, lim), aspect_ratio = :equal, + xlabel = "electron 1 position z₁", ylabel = "electron 2 position z₂", + colorbar_title = "\n|Ψ|² (normalised)", colorbar_titlefontsize = 9, + title = "Hooke's atom — two-electron density", + titlefontsize = 12, guidefontsize = 10, tickfontsize = 9, + grid = false, framestyle = :box, widen = false, + size = (600, 500), dpi = 200, left_margin = 3Plots.mm, +) +## dashed line marks where the electrons coincide — the trough is the Coulomb hole +plot!( + [-lim, lim], [-lim, lim]; + color = :white, alpha = 0.4, linestyle = :dash, linewidth = 1.5, + label = "z₁ = z₂ (Coulomb hole)", legend = :topleft, foreground_color_legend = nothing, + background_color_legend = RGBA(0, 0, 0, 0.35), legendfontcolor = :white, legendfontsize = 8, +) diff --git a/test/test_harmonium.jl b/test/test_harmonium.jl index 0e4c4fc..721a528 100644 --- a/test/test_harmonium.jl +++ b/test/test_harmonium.jl @@ -17,4 +17,17 @@ using FewBodyECG @test sol.E₀ ≈ 2.0 atol = 1.0e-3 # Taut's exact energy @test sol.E₀ > 2.0 - 1.0e-9 # variational upper bound never violated + + # Relative wavefunction vs Taut's closed form χ(u) ∝ (1 + u/2) e^{-u²/8}. + # Transform physical electron positions back to Jacobi coordinates via J. + ψ = wavefunction(sol) + J, _ = jacobi_transform([1.0e15, 1.0, 1.0]) + Ψ(z₁, z₂) = ψ(J * [0.0, z₁, z₂]) + u = range(0, 10, length = 200) + χ_ecg = [Ψ(x / 2, -x / 2) for x in u] + χ_exact = [(1 + x / 2) * exp(-x^2 / 8) for x in u] + χ_ecg ./= maximum(abs, χ_ecg) + χ_exact ./= maximum(χ_exact) + χ_ecg .*= sign(sum(χ_ecg .* χ_exact)) + @test maximum(abs, χ_ecg .- χ_exact) < 0.02 end From 8cd9df6404eea0bf3c6f883ff8a4370437089c2e Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Fri, 10 Jul 2026 17:26:00 +0200 Subject: [PATCH 3/5] updated examples --- docs/make.jl | 2 ++ docs/src/examples/muonic_ions.md | 40 ++++++++++++++++++++++ docs/src/examples/nuclear.md | 58 ++++++++++++++++++++++++++++++++ examples/muonic_ions.jl | 29 ++++++++++++++++ examples/nuclear.jl | 47 ++++++++++++++++++++++++++ test/runtests.jl | 2 ++ test/test_muonic.jl | 26 ++++++++++++++ test/test_nuclear.jl | 37 ++++++++++++++++++++ 8 files changed, 241 insertions(+) create mode 100644 docs/src/examples/muonic_ions.md create mode 100644 docs/src/examples/nuclear.md create mode 100644 examples/muonic_ions.jl create mode 100644 examples/nuclear.jl create mode 100644 test/test_muonic.jl create mode 100644 test/test_nuclear.jl diff --git a/docs/make.jl b/docs/make.jl index ba9d636..aaa42e0 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -23,9 +23,11 @@ makedocs( "Positronium" => "examples/positronium.md", "Helium and H-" => "examples/helium.md", "tdmu" => "examples/tdmu.md", + "Muonic ions" => "examples/muonic_ions.md", "H2+ (non-BO)" => "examples/h2plus.md", "Gaussian wells" => "examples/gaussian_well.md", "Hooke's atom" => "examples/harmonium.md", + "Nuclear (Minnesota/Volkov)" => "examples/nuclear.md", "Three-body Gaussian" => "examples/three_body_gaussian.md", "Spin-orbit doublet" => "examples/spin_orbit.md", "Workflow" => "examples/workflow.md", diff --git a/docs/src/examples/muonic_ions.md b/docs/src/examples/muonic_ions.md new file mode 100644 index 0000000..5d936e2 --- /dev/null +++ b/docs/src/examples/muonic_ions.md @@ -0,0 +1,40 @@ +```@meta +EditURL = "../../../examples/muonic_ions.jl" +``` + +# Muonic molecular ions (ddμ, dtμ, ttμ) + +A negative muon binds two hydrogen-isotope nuclei into a tiny molecular ion — +the mechanism behind muon-catalyzed fusion. Each is a three-body Coulomb +problem `Operators(masses, charges)` away, differing only in the nuclear masses +(electron-mass atomic units). The heavy, deeply-bound, two-scale nature makes +these among the hardest ECG problems; the Suzuki–Varga K = 200 benchmarks are +reproduced here to a few mHa. + +````@example muonic_ions +using FewBodyECG + +mμ, md, mt = 206.7686, 3670.481, 5496.918 +systems = [ + ("ddμ", [md, md, mμ], NaN), + ("dtμ", [md, mt, mμ], -111.36444), # Suzuki–Varga Table 8.1 + ("ttμ", [mt, mt, mμ], -112.973), # Suzuki–Varga Table 8.1 +] + +for (name, masses, ref) in systems + ops = Operators(masses, [+1.0, +1.0, -1.0]) + ops += "Kinetic" + ops += "Coulomb" + sol = solve(ops, SVM(basis = 300, candidates = 40, scale = 0.02); tol = 1.0e-4, window = 15) + msg = isnan(ref) ? "(prediction)" : "ref $ref Δ = $(round(sol.E₀ - ref, digits = 5))" + println(rpad(name, 5), " E₀ = ", round(sol.E₀, digits = 5), " Ha ", msg) +end +```` + +Heavier nuclei localise the muon more tightly, so the binding deepens along +ddμ → dtμ → ttμ. + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + diff --git a/docs/src/examples/nuclear.md b/docs/src/examples/nuclear.md new file mode 100644 index 0000000..6a4380e --- /dev/null +++ b/docs/src/examples/nuclear.md @@ -0,0 +1,58 @@ +```@meta +EditURL = "../../../examples/nuclear.jl" +``` + +# Nuclear few-body with Gaussian potentials + +Realistic soft-core nucleon–nucleon potentials are **sums of Gaussians**, so +[`GaussianOperator`](@ref) represents them exactly — opening the nuclear domain +with no new machinery. We work in nuclear units (energies in MeV, lengths in +fm); a particle's package "mass" is `mc² / (ħc)²` so that `ħ²/2m → (ħc)²/2mc²` +with `ħc = 197.327 MeV·fm` (giving `ħ²/2mₙ ≈ 20.7 MeV·fm²`). + +Two classic central-force benchmarks: +* the **deuteron** with the Minnesota potential (which reproduces the two-body + binding), and +* the **triton** with the Volkov V1 potential (a standard three-body central + benchmark). + +````@example nuclear +using FewBodyECG + +const ħc = 197.3269804 +mpkg(mc²) = mc² / ħc^2 +mp, mn = mpkg(938.272), mpkg(939.565) + +# Deuteron — Minnesota triplet-even central potential +# V(r) = 200 e^{-1.487 r²} − 178 e^{-0.639 r²} (MeV, r in fm) +deut = Operators([mp, mn]) +deut += "Kinetic" +deut += ("Gaussian", 1, 2, 200.0, 1.487) # repulsive core +deut += ("Gaussian", 1, 2, -178.0, 0.639) # triplet attraction + +sol_d = solve(deut, SVM(basis = 40, candidates = 25, scale = 3.0)) +Ed_ref = -2.202 +println("deuteron E = ", round(sol_d.E₀, digits = 4), " MeV (Minnesota ", Ed_ref, ")") + +# Triton — Volkov V1 central potential on all three pairs +# V(r) = 144.86 e^{-(r/0.82)²} − 83.34 e^{-(r/1.60)²} (MeV, r in fm) +γR, γA = 1 / 0.82^2, 1 / 1.6^2 +trit = Operators([mn, mn, mp]) +trit += "Kinetic" +trit += ("Gaussian", 1, 2, 144.86, γR); trit += ("Gaussian", 1, 2, -83.34, γA) +trit += ("Gaussian", 1, 3, 144.86, γR); trit += ("Gaussian", 1, 3, -83.34, γA) +trit += ("Gaussian", 2, 3, 144.86, γR); trit += ("Gaussian", 2, 3, -83.34, γA) + +sol_t = solve(trit, SVM(basis = 120, candidates = 30, scale = 3.5)) +Et_ref = -8.46 +println("triton E = ", round(sol_t.E₀, digits = 4), " MeV (Volkov V1 ", Et_ref, ")") +```` + +The triton's central-force binding is ~4× the deuteron's — the classic +few-body sequence. (The physical triton also needs spin-isospin and tensor +forces, which a purely central model omits.) + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + diff --git a/examples/muonic_ions.jl b/examples/muonic_ions.jl new file mode 100644 index 0000000..a0b14a8 --- /dev/null +++ b/examples/muonic_ions.jl @@ -0,0 +1,29 @@ +# # Muonic molecular ions (ddμ, dtμ, ttμ) +# +# A negative muon binds two hydrogen-isotope nuclei into a tiny molecular ion — +# the mechanism behind muon-catalyzed fusion. Each is a three-body Coulomb +# problem `Operators(masses, charges)` away, differing only in the nuclear masses +# (electron-mass atomic units). The heavy, deeply-bound, two-scale nature makes +# these among the hardest ECG problems; the Suzuki–Varga K = 200 benchmarks are +# reproduced here to a few mHa. + +using FewBodyECG + +mμ, md, mt = 206.7686, 3670.481, 5496.918 +systems = [ + ("ddμ", [md, md, mμ], NaN), + ("dtμ", [md, mt, mμ], -111.36444), # Suzuki–Varga Table 8.1 + ("ttμ", [mt, mt, mμ], -112.973), # Suzuki–Varga Table 8.1 +] + +for (name, masses, ref) in systems + ops = Operators(masses, [+1.0, +1.0, -1.0]) + ops += "Kinetic" + ops += "Coulomb" + sol = solve(ops, SVM(basis = 300, candidates = 40, scale = 0.02); tol = 1.0e-4, window = 15) + msg = isnan(ref) ? "(prediction)" : "ref $ref Δ = $(round(sol.E₀ - ref, digits = 5))" + println(rpad(name, 5), " E₀ = ", round(sol.E₀, digits = 5), " Ha ", msg) +end + +# Heavier nuclei localise the muon more tightly, so the binding deepens along +# ddμ → dtμ → ttμ. diff --git a/examples/nuclear.jl b/examples/nuclear.jl new file mode 100644 index 0000000..eede6fb --- /dev/null +++ b/examples/nuclear.jl @@ -0,0 +1,47 @@ +# # Nuclear few-body with Gaussian potentials +# +# Realistic soft-core nucleon–nucleon potentials are **sums of Gaussians**, so +# [`GaussianOperator`](@ref) represents them exactly — opening the nuclear domain +# with no new machinery. We work in nuclear units (energies in MeV, lengths in +# fm); a particle's package "mass" is `mc² / (ħc)²` so that `ħ²/2m → (ħc)²/2mc²` +# with `ħc = 197.327 MeV·fm` (giving `ħ²/2mₙ ≈ 20.7 MeV·fm²`). +# +# Two classic central-force benchmarks: +# * the **deuteron** with the Minnesota potential (which reproduces the two-body +# binding), and +# * the **triton** with the Volkov V1 potential (a standard three-body central +# benchmark). + +using FewBodyECG + +const ħc = 197.3269804 +mpkg(mc²) = mc² / ħc^2 +mp, mn = mpkg(938.272), mpkg(939.565) + +## Deuteron — Minnesota triplet-even central potential +## V(r) = 200 e^{-1.487 r²} − 178 e^{-0.639 r²} (MeV, r in fm) +deut = Operators([mp, mn]) +deut += "Kinetic" +deut += ("Gaussian", 1, 2, 200.0, 1.487) # repulsive core +deut += ("Gaussian", 1, 2, -178.0, 0.639) # triplet attraction + +sol_d = solve(deut, SVM(basis = 40, candidates = 25, scale = 3.0)) +Ed_ref = -2.202 +println("deuteron E = ", round(sol_d.E₀, digits = 4), " MeV (Minnesota ", Ed_ref, ")") + +## Triton — Volkov V1 central potential on all three pairs +## V(r) = 144.86 e^{-(r/0.82)²} − 83.34 e^{-(r/1.60)²} (MeV, r in fm) +γR, γA = 1 / 0.82^2, 1 / 1.6^2 +trit = Operators([mn, mn, mp]) +trit += "Kinetic" +trit += ("Gaussian", 1, 2, 144.86, γR); trit += ("Gaussian", 1, 2, -83.34, γA) +trit += ("Gaussian", 1, 3, 144.86, γR); trit += ("Gaussian", 1, 3, -83.34, γA) +trit += ("Gaussian", 2, 3, 144.86, γR); trit += ("Gaussian", 2, 3, -83.34, γA) + +sol_t = solve(trit, SVM(basis = 120, candidates = 30, scale = 3.5)) +Et_ref = -8.46 +println("triton E = ", round(sol_t.E₀, digits = 4), " MeV (Volkov V1 ", Et_ref, ")") + +# The triton's central-force binding is ~4× the deuteron's — the classic +# few-body sequence. (The physical triton also needs spin-isospin and tensor +# forces, which a purely central model omits.) diff --git a/test/runtests.jl b/test/runtests.jl index 96c85f5..9ab2c1a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -22,6 +22,8 @@ using FewBodyECG include("test_harmonic_oscillator.jl") include("test_harmonium.jl") include("test_spin_orbit.jl") + include("test_nuclear.jl") + include("test_muonic.jl") include("test_utils.jl") include("test_types.jl") include("test_variational.jl") diff --git a/test/test_muonic.jl b/test/test_muonic.jl new file mode 100644 index 0000000..ea5c6a7 --- /dev/null +++ b/test/test_muonic.jl @@ -0,0 +1,26 @@ +using Test +using FewBodyECG + +# Muonic molecular ions: three-body Coulomb systems reproducing the Suzuki–Varga +# Table 8.1 K = 200 benchmark energies. A moderate basis is used here to keep +# the suite fast (the example shows the few-mHa agreement at basis = 300). +@testset "Muonic molecular ions vs Suzuki–Varga" begin + mμ, md, mt = 206.7686, 3670.481, 5496.918 + + solve_ion(masses) = solve( + (o = Operators(masses, [+1.0, +1.0, -1.0]); o += "Kinetic"; o += "Coulomb"; o), + SVM(basis = 200, candidates = 40, scale = 0.02); tol = 1.0e-4, window = 15, + ).E₀ + + E_dt = solve_ion([md, mt, mμ]) + E_tt = solve_ion([mt, mt, mμ]) + + # reproduce the benchmarks from above (variational upper bound) to <40 mHa + @test E_dt ≈ -111.36444 atol = 0.04 + @test E_dt > -111.36444 # variational upper bound + @test E_tt ≈ -112.973 atol = 0.04 + @test E_tt > -112.973 # variational upper bound + + # heavier nuclei bind deeper: ttμ below dtμ + @test E_tt < E_dt +end diff --git a/test/test_nuclear.jl b/test/test_nuclear.jl new file mode 100644 index 0000000..f27db39 --- /dev/null +++ b/test/test_nuclear.jl @@ -0,0 +1,37 @@ +using Test +using FewBodyECG + +# Nuclear few-body benchmarks with Gaussian NN potentials (energies in MeV, +# lengths in fm; particle "mass" = mc²/(ħc)² so ħ²/2m → (ħc)²/2mc²). +@testset "Nuclear Gaussian potentials (deuteron, triton)" begin + ħc = 197.3269804 + mpkg(mc²) = mc² / ħc^2 + mp, mn = mpkg(938.272), mpkg(939.565) + + # ħ²/2mₙ ≈ 20.7 MeV·fm² for a nucleon + @test 1 / (2 * mpkg(938.918)) ≈ 20.736 atol = 1.0e-2 + + # Deuteron — Minnesota triplet-even central: benchmark ≈ -2.202 MeV + deut = Operators([mp, mn]) + deut += "Kinetic" + deut += ("Gaussian", 1, 2, 200.0, 1.487) + deut += ("Gaussian", 1, 2, -178.0, 0.639) + sol_d = solve(deut, SVM(basis = 40, candidates = 25, scale = 3.0)) + @test sol_d.E₀ ≈ -2.202 atol = 5.0e-3 + @test sol_d.E₀ > -2.202 - 1.0e-3 # variational upper bound + + # Triton — Volkov V1 central on all pairs: benchmark ≈ -8.46 MeV + γR, γA = 1 / 0.82^2, 1 / 1.6^2 + trit = Operators([mn, mn, mp]) + trit += "Kinetic" + for (i, j) in ((1, 2), (1, 3), (2, 3)) + trit += ("Gaussian", i, j, 144.86, γR) + trit += ("Gaussian", i, j, -83.34, γA) + end + sol_t = solve(trit, SVM(basis = 100, candidates = 30, scale = 3.5)) + @test sol_t.E₀ ≈ -8.46 atol = 2.0e-2 + @test sol_t.E₀ > -8.48 # variational upper bound + + # the three-body system is bound well below the two-body deuteron + @test sol_t.E₀ < sol_d.E₀ - 5.0 +end From 0b9e683fc2bbc2431b7566019aa3e77f688d6d10 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Sat, 11 Jul 2026 13:05:30 +0200 Subject: [PATCH 4/5] added 4 particle example --- examples/hps.jl | 27 +++++++++++++++++++++++ examples/muonic_ions.jl | 47 ++++++++++++++++++++++------------------- examples/ps2.jl | 24 +++++++++++++++++++++ 3 files changed, 76 insertions(+), 22 deletions(-) create mode 100644 examples/hps.jl create mode 100644 examples/ps2.jl diff --git a/examples/hps.jl b/examples/hps.jl new file mode 100644 index 0000000..d9b68d6 --- /dev/null +++ b/examples/hps.jl @@ -0,0 +1,27 @@ +# # Positronium hydride (HPs⁺) +# +# A four-body Coulomb system consisting of a fixed proton, two electrons, and +# a positron. The particle order is proton, electron, electron, positron. +# This rank-0 example targets the symmetric spatial ground sector of the two +# identical electrons. The reference is the infinite-proton-mass SVM value +# from Suzuki–Varga Table 8.6 (K = 1200); this run uses a smaller basis as a +# practical example rather than claiming to reproduce K = 1200. + +using FewBodyECG +using Plots + +masses = [1.0e15, 1.0, 1.0, 1.0] +charges = [+1.0, -1.0, -1.0, +1.0] + +ops = Operators(masses, charges) +ops += "Kinetic" +ops += "Coulomb" + +hps_ref = -0.7891964 +sol = solve(ops, GrowVariational(basis = 50, candidates = 20, scale = 1.0)) + +println("HPs⁺ E₀ = ", sol.E₀, " Ha") +println("reference (SVM, K = 1200) = ", hps_ref, " Ha Δ = ", sol.E₀ - hps_ref) +println("variational upper bound respected: ", sol.E₀ ≥ hps_ref) + +plot(sol, hps_ref; title = "Positronium-hydride convergence") diff --git a/examples/muonic_ions.jl b/examples/muonic_ions.jl index a0b14a8..a5bdf56 100644 --- a/examples/muonic_ions.jl +++ b/examples/muonic_ions.jl @@ -1,29 +1,32 @@ -# # Muonic molecular ions (ddμ, dtμ, ttμ) -# -# A negative muon binds two hydrogen-isotope nuclei into a tiny molecular ion — -# the mechanism behind muon-catalyzed fusion. Each is a three-body Coulomb -# problem `Operators(masses, charges)` away, differing only in the nuclear masses -# (electron-mass atomic units). The heavy, deeply-bound, two-scale nature makes -# these among the hardest ECG problems; the Suzuki–Varga K = 200 benchmarks are -# reproduced here to a few mHa. - -using FewBodyECG - -mμ, md, mt = 206.7686, 3670.481, 5496.918 systems = [ - ("ddμ", [md, md, mμ], NaN), - ("dtμ", [md, mt, mμ], -111.36444), # Suzuki–Varga Table 8.1 - ("ttμ", [mt, mt, mμ], -112.973), # Suzuki–Varga Table 8.1 + ( + "dtμ", + [md, mt, mμ], + -111.36444, # Suzuki–Varga SVM, K = 200 + -111.364511474, # other method, K = 1400 + ), + ( + "ttμ", + [mt, mt, mμ], + -112.97300, # Suzuki–Varga SVM, K = 200 + -112.9730179, # other method, K = 500 + ), ] -for (name, masses, ref) in systems +for (name, masses, svm200, best_ref) in systems ops = Operators(masses, [+1.0, +1.0, -1.0]) ops += "Kinetic" ops += "Coulomb" - sol = solve(ops, SVM(basis = 300, candidates = 40, scale = 0.02); tol = 1.0e-4, window = 15) - msg = isnan(ref) ? "(prediction)" : "ref $ref Δ = $(round(sol.E₀ - ref, digits = 5))" - println(rpad(name, 5), " E₀ = ", round(sol.E₀, digits = 5), " Ha ", msg) -end -# Heavier nuclei localise the muon more tightly, so the binding deepens along -# ddμ → dtμ → ttμ. + sol = solve( + ops, + SVM(basis = 200, candidates = 40, scale = 0.02); + tol = 1.0e-4, + window = 15, + ) + + println(name) + println(" ECG E₀ = ", sol.E₀, " Ha") + println(" Δ vs SVM K=200 = ", sol.E₀ - svm200, " Ha") + println(" Δ vs high-K result = ", sol.E₀ - best_ref, " Ha") +end \ No newline at end of file diff --git a/examples/ps2.jl b/examples/ps2.jl new file mode 100644 index 0000000..1eeb0cb --- /dev/null +++ b/examples/ps2.jl @@ -0,0 +1,24 @@ +using FewBodyECG +using Plots + +masses = [1.0, 1.0, 1.0, 1.0] +charges = [+1.0, -1.0, +1.0, -1.0] + +ops = Operators(masses, charges) +ops += "Kinetic" +ops += "Coulomb" + +ps2_ref = -0.516003778 +sol = solve( + ops, + SVM(basis = 150, candidates = 30, scale = 1.0); + tol = 1.0e-4, + window = 15, +) + +println("Ps₂ E₀ = ", sol.E₀, " Ha") +println("reference (SVM, K = 800) = ", ps2_ref, " Ha Δ = ", sol.E₀ - ps2_ref) +println("variational upper bound respected: ", sol.E₀ ≥ ps2_ref) + +plot(sol, ps2_ref; title = "Dipositronium convergence") +plot(wavefunction(sol)) \ No newline at end of file From 440f49898184e968c8941d9c69d974cce9c9747a Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Sat, 11 Jul 2026 15:49:08 +0200 Subject: [PATCH 5/5] updated tests --- docs/src/examples/muonic_ions.md | 46 +++++++++++++++++++----------- examples/hps.jl | 10 +------ examples/muonic_ions.jl | 11 +++++++- test/test_operators.jl | 48 ++++++++++++++++++++++++++++++++ 4 files changed, 88 insertions(+), 27 deletions(-) diff --git a/docs/src/examples/muonic_ions.md b/docs/src/examples/muonic_ions.md index 5d936e2..8a04275 100644 --- a/docs/src/examples/muonic_ions.md +++ b/docs/src/examples/muonic_ions.md @@ -2,38 +2,50 @@ EditURL = "../../../examples/muonic_ions.jl" ``` -# Muonic molecular ions (ddμ, dtμ, ttμ) +# Muonic molecular ions (dtμ, ttμ) -A negative muon binds two hydrogen-isotope nuclei into a tiny molecular ion — -the mechanism behind muon-catalyzed fusion. Each is a three-body Coulomb -problem `Operators(masses, charges)` away, differing only in the nuclear masses -(electron-mass atomic units). The heavy, deeply-bound, two-scale nature makes -these among the hardest ECG problems; the Suzuki–Varga K = 200 benchmarks are -reproduced here to a few mHa. +Compare the ECG results with both the Suzuki–Varga K = 200 SVM values and +the higher-basis values listed in Table 8.1. ````@example muonic_ions using FewBodyECG mμ, md, mt = 206.7686, 3670.481, 5496.918 + systems = [ - ("ddμ", [md, md, mμ], NaN), - ("dtμ", [md, mt, mμ], -111.36444), # Suzuki–Varga Table 8.1 - ("ttμ", [mt, mt, mμ], -112.973), # Suzuki–Varga Table 8.1 + ( + "dtμ", + [md, mt, mμ], + -111.36444, # Suzuki–Varga SVM, K = 200 + -111.364511474, # other method, K = 1400 + ), + ( + "ttμ", + [mt, mt, mμ], + -112.97300, # Suzuki–Varga SVM, K = 200 + -112.9730179, # other method, K = 500 + ), ] -for (name, masses, ref) in systems +for (name, masses, svm200, best_ref) in systems ops = Operators(masses, [+1.0, +1.0, -1.0]) ops += "Kinetic" ops += "Coulomb" - sol = solve(ops, SVM(basis = 300, candidates = 40, scale = 0.02); tol = 1.0e-4, window = 15) - msg = isnan(ref) ? "(prediction)" : "ref $ref Δ = $(round(sol.E₀ - ref, digits = 5))" - println(rpad(name, 5), " E₀ = ", round(sol.E₀, digits = 5), " Ha ", msg) + + sol = solve( + ops, + SVM(basis = 200, candidates = 40, scale = 0.02); + tol = 1.0e-4, + window = 15, + ) + + println(name) + println(" ECG E₀ = ", sol.E₀, " Ha") + println(" Δ vs SVM K=200 = ", sol.E₀ - svm200, " Ha") + println(" Δ vs high-K result = ", sol.E₀ - best_ref, " Ha") end ```` -Heavier nuclei localise the muon more tightly, so the binding deepens along -ddμ → dtμ → ttμ. - --- *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* diff --git a/examples/hps.jl b/examples/hps.jl index d9b68d6..45d6d88 100644 --- a/examples/hps.jl +++ b/examples/hps.jl @@ -1,12 +1,3 @@ -# # Positronium hydride (HPs⁺) -# -# A four-body Coulomb system consisting of a fixed proton, two electrons, and -# a positron. The particle order is proton, electron, electron, positron. -# This rank-0 example targets the symmetric spatial ground sector of the two -# identical electrons. The reference is the infinite-proton-mass SVM value -# from Suzuki–Varga Table 8.6 (K = 1200); this run uses a smaller basis as a -# practical example rather than claiming to reproduce K = 1200. - using FewBodyECG using Plots @@ -25,3 +16,4 @@ println("reference (SVM, K = 1200) = ", hps_ref, " Ha Δ = ", sol.E₀ - hps_r println("variational upper bound respected: ", sol.E₀ ≥ hps_ref) plot(sol, hps_ref; title = "Positronium-hydride convergence") +plot(wavefunction(sol)) \ No newline at end of file diff --git a/examples/muonic_ions.jl b/examples/muonic_ions.jl index a5bdf56..96feb1d 100644 --- a/examples/muonic_ions.jl +++ b/examples/muonic_ions.jl @@ -1,3 +1,12 @@ +# # Muonic molecular ions (dtμ, ttμ) +# +# Compare the ECG results with both the Suzuki–Varga K = 200 SVM values and +# the higher-basis values listed in Table 8.1. + +using FewBodyECG + +mμ, md, mt = 206.7686, 3670.481, 5496.918 + systems = [ ( "dtμ", @@ -29,4 +38,4 @@ for (name, masses, svm200, best_ref) in systems println(" ECG E₀ = ", sol.E₀, " Ha") println(" Δ vs SVM K=200 = ", sol.E₀ - svm200, " Ha") println(" Δ vs high-K result = ", sol.E₀ - best_ref, " Ha") -end \ No newline at end of file +end diff --git a/test/test_operators.jl b/test/test_operators.jl index 2b1f7d4..e1aef06 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -257,6 +257,54 @@ import FewBodyECG: jacobi_transform, Λ @test occursin("1 term", str) @test !occursin("1 terms", str) end + + @testset "GaussianOperator branch" begin + ops = Operators([1.0e15, 1.0]) + ops += ("Gaussian", 1, 2, -5.0, 1.0) + str = sprint(show, ops) + @test occursin("Gaussian(γ =", str) + @test occursin("-5.0", str) + @test occursin("w =", str) + end + + @testset "OscillatorOperator branch" begin + ops = Operators([1.0e15, 1.0]) + ops += ("Oscillator", 1, 2, 0.5) + str = sprint(show, ops) + @test occursin("Oscillator(w =", str) + @test occursin("0.5", str) + end + + @testset "ManyBodyGaussianOperator branch" begin + ops = Operators() + ops += ManyBodyGaussianOperator(2.0, [0.4;;]) + str = sprint(show, ops) + @test occursin("ManyBodyGaussian(W =", str) + @test occursin("2.0", str) + end + + @testset "Fallback branch for other operator types" begin + ops = Operators() + ops += GaussianTensorOperator(1.0, 1.0, [1.0], 1, 2) + str = sprint(show, ops) + @test occursin("GaussianTensorOperator", str) + end + + @testset "All branches together" begin + ops = Operators([1.0e15, 1.0], [+1, -1]) + ops += "Kinetic" + ops += "Coulomb" + ops += ("Gaussian", 1, 2, -5.0, 1.0) + ops += ("Oscillator", 1, 2, 0.5) + ops += ManyBodyGaussianOperator(2.0, [0.4;;]) + str = sprint(show, ops) + @test occursin("5 terms", str) + @test occursin("Kinetic", str) + @test occursin("Coulomb", str) + @test occursin("Gaussian(γ =", str) + @test occursin("Oscillator(w =", str) + @test occursin("ManyBodyGaussian(W =", str) + end end @testset "coulomb_weights" begin