Variational method

Examples/Helium.jl

@@ -0,0 +1,66 @@
+using FewBodyECG
+using Antique
+using Plots
+using QuasiMonteCarlo
+
+masses = [1e15, 1.0, 1.0]
+
+os = Operators(masses, [+2, -1, -1])   # nucleus (Z=2), e₁, e₂
+os += "Kinetic"
+os += "Coulomb"   # auto: nucleus-e₁ (-2), nucleus-e₂ (-2), e₁-e₂ (+1)
+
+E_gs_exact = -2.9037242   # 1s² ¹S  (ground state)
+E_ex_exact = -2.17523     # 1s2s ¹S (first excited singlet)
+
+println("Ground state (1s²)...")
+result_gs = solve_ECG(os, 250; scale = 1.0, verbose = false)
+ΔE_gs = result_gs.ground_state - E_gs_exact
+println("  E = $(round(result_gs.ground_state; digits=6))  (exact: $E_gs_exact,  error: $(round(ΔE_gs; digits=6)))")
+
+println("\nFirst excited ¹S state (1s2s)...")
+result_ex = solve_ECG(os, 200; scale = 1.0, verbose = false, state = 2)
+ΔE_ex = result_ex.ground_state - E_ex_exact
+println("  E = $(round(result_ex.ground_state; digits=6))  (exact: $E_ex_exact,  error: $(round(ΔE_ex; digits=6)))")
+
+n_gs, E_gs = convergence(result_gs)
+n_ex, E_ex = convergence(result_ex)
+
+p1 = plot(n_gs, E_gs,
+    label = "1s² (ground)", lw = 2,
+    xlabel = "Basis size", ylabel = "E (Ha)",
+    title = "Helium convergence")
+plot!(p1, n_ex, E_ex, label = "1s2s (excited)", lw = 2, ls = :dash)
+hline!(p1, [E_gs_exact, E_ex_exact], ls = :dot, color = :gray, label = "Exact")
+display(p1)
+
+r_gs, ρ_gs = correlation_function(result_gs; rmax = 5.0)
+r_ex, ρ_ex = correlation_function(result_ex; rmax = 5.0)
+
+p2 = plot(r_gs, ρ_gs,
+    label = "1s² (ground)", lw = 2,
+    xlabel = "r (a.u.)", ylabel = "r²|ψ(r)|²",
+    title = "Helium radial correlation")
+    plot!(p2, r_ex, ρ_ex, label = "1s2s (excited)", lw = 2, ls = :dash)
+display(p2)
+
+HeP = HydrogenAtom(Z = 2)   # atomic units: Eₕ=1, a₀=1, mₑ=1, ℏ=1 (defaults)
+E_hep_exact = Antique.E(HeP; n = 1)   # = -2.0 Ha
+
+os_hep = Operators([1e15, 1.0], [+2, -1])
+os_hep += "Kinetic"
+os_hep += "Coulomb"
+
+result_hep = solve_ECG(os_hep, 250; scale = 0.5, verbose = false)
+ΔE_hep = result_hep.ground_state - E_hep_exact
+println("  E (ECG)    = $(round(result_hep.ground_state; digits=8))")
+println("  E (Antique)= $(round(E_hep_exact; digits=8))   error: $(round(ΔE_hep; sigdigits=3))")
+
+r_hep, ρ_ecg = correlation_function(result_hep; rmax = 3.0, npoints = 300)
+ρ_antique = [r^2 * Antique.R(HeP, r; n = 1, l = 0)^2 for r in r_hep]
+
+p3 = plot(r_hep, ρ_ecg,
+    label = "ECG", lw = 2,
+    xlabel = "r (a.u.)", ylabel = "r²|ψ(r)|²",
+    title = "He⁺ 1s radial density: ECG vs Antique.jl")
+plot!(p3, r_hep, ρ_antique, label = "Antique (exact)", lw = 2, ls = :dash)
+display(p3)

Examples/HydrogenAnion.jl

@@ -5,25 +5,26 @@ using QuasiMonteCarlo
 
 masses = [1.0e15, 1.0, 1.0]
 
-Λmat = Λ(masses)
-kin = KineticOperator(Λmat)
-J, U = _jacobi_transform(masses)
+os = Operators(masses, [+1, -1, -1])   # proton, e₁, e₂
+os += "Kinetic"
+os += "Coulomb"
 
-w_list = [[1, -1, 0], [1, 0, -1], [0, 1, -1]]
+result = solve_ECG(os, 250, scale = 1.0, verbose=false)
 
-w_raw = [U' * w for w in w_list]
-coeffs = [-1.0, -1.0, +1.0]
-
-ops = Operator[
-    kin;
-    (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))...
-]
-
-result = solve_ECG(ops, 250, scale = 1.0)
-
-E = -0.527751016523
-ΔE = abs(result.ground_state - E)
+E_exact = -0.527751016523
+ΔE = abs(result.ground_state - E_exact)
 @info "Energy difference" ΔE
 
-n, E = convergence(result)
-plot(n, E)
+n_conv, E_conv = convergence(result)
+p1 = plot(n_conv, E_conv,
+    xlabel = "Basis size", ylabel = "E (Ha)",
+    label = "Ground state", lw = 2,
+    title = "Hydrogen anion convergence")
+display(p1)
+
+r_grid, ρ = correlation_function(result; rmax = 10.0)
+p2 = plot(r_grid, ρ,
+    xlabel = "r (a.u.)", ylabel = "r²|ψ(r)|²",
+    label = "Hydrogen anion", lw = 2,
+    title = "Hydrogen radial correlation")
+display(p2)

Examples/HydrogenStates.jl

@@ -0,0 +1,163 @@
+using FewBodyECG
+import Antique
+using LinearAlgebra
+using Plots
+
+atom = Antique.HydrogenAtom()
+
+masses = [1.0e15, 1.0]
+Λmat   = Λ(masses)
+_, U   = _jacobi_transform(masses)
+w      = U' * Float64.([1, -1])   # electron-proton separation in Jacobi coords
+
+ops = Operator[KineticOperator(Λmat); CoulombOperator(-1.0, w)]
+
+sr_1s = solve_ECG_sequential(ops, 16;
+    n_candidates = 8, scale = 1.0, max_iterations_step = 120, verbose = true)
+
+E_1s = sr_1s.ground_state
+println("\n  E(1s)       = $(round(E_1s; digits = 8)) Ha")
+println("  E_exact(1s) = $(round(Antique.E(atom; n = 1); digits = 8)) Ha")
+println("  |ΔE|        = $(round(abs(E_1s - Antique.E(atom; n = 1)); sigdigits = 2))\n")
+
+a_p    = [1.0]
+s_zero = [0.0]
+
+# Log-spaced widths covering the spatial extent of the 2p orbital (peak ~4 a₀)
+alphas_p = exp10.(range(log10(0.003), log10(3.0), length = 16))
+basis_p  = GaussianBase[]
+E_2p_conv = Float64[]
+vecs_2p   = Matrix{Float64}(undef, 0, 0)
+
+for (k, α) in enumerate(alphas_p)
+    push!(basis_p, Rank1Gaussian([α;;], a_p, s_zero))
+    bset = BasisSet(basis_p)
+    H = build_hamiltonian_matrix(bset, ops)
+    S = build_overlap_matrix(bset)
+    vals, vecs = solve_generalized_eigenproblem(H, S)
+    E_k = minimum(vals)
+    push!(E_2p_conv, E_k)
+    global vecs_2p = vecs
+    println("  step $(lpad(k, 2))  α = $(rpad(round(α; digits = 4), 7))" *
+            "  E = $(round(E_k; digits = 8))")
+end
+
+E_2p = minimum(E_2p_conv)
+c_2p = vecs_2p[:, 1]
+println("\n  E(2p)       = $(round(E_2p; digits = 8)) Ha")
+println("  E_exact(2p) = $(round(Antique.E(atom; n = 2); digits = 8)) Ha")
+println("  |ΔE|        = $(round(abs(E_2p - Antique.E(atom; n = 2)); sigdigits = 2))\n")
+
+# Orthogonal polarization vectors define a pure d-wave channel.
+a_d = reshape([1.0, 0.0, 0.0], 1, 3)
+b_d = reshape([0.0, 0.0, 1.0], 1, 3)
+alphas_d = exp10.(range(log10(0.002), log10(0.8), length = 24))
+basis_d   = GaussianBase[]
+E_3d_conv = Float64[]
+vecs_3d   = Matrix{Float64}(undef, 0, 0)
+
+for (k, α) in enumerate(alphas_d)
+    push!(basis_d, Rank2Gaussian([α;;], a_d, b_d, s_zero))
+    bset = BasisSet(basis_d)
+    H = build_hamiltonian_matrix(bset, ops)
+    S = build_overlap_matrix(bset)
+    vals, vecs = solve_generalized_eigenproblem(H, S)
+    E_k = minimum(vals)
+    push!(E_3d_conv, E_k)
+    global vecs_3d = vecs
+    println("  step $(lpad(k, 2))  α = $(rpad(round(α; digits = 4), 7))" *
+            "  E = $(round(E_k; digits = 8))")
+end
+
+E_rank2 = minimum(E_3d_conv)
+c_rank2 = vecs_3d[:, 1]
+println("\n  E(3d, Rank2 pure d) = $(round(E_rank2; digits = 8)) Ha")
+println("  E_exact(3d)         = $(round(Antique.E(atom; n = 3); digits = 8)) Ha")
+println("  |ΔE|                = $(round(abs(E_rank2 - Antique.E(atom; n = 3)); sigdigits = 2))\n")
+
+function ψ_rank1(rval, c, bfs)
+    r_vec = [rval]
+    return sum(
+        c[i] * dot(bfs[i].a, r_vec) *
+            exp(-dot(r_vec, parent(bfs[i].A) * r_vec) + dot(bfs[i].s, r_vec))
+            for i in eachindex(bfs)
+    )
+end
+
+function ψ_rank2(rval, c, bfs)
+    # Directional profile for the pure d-wave basis: r̂ = (x+z)/√2.
+    θ_d, φ_d = π / 4, 0.0
+    r_cart = rval .* [sin(θ_d) * cos(φ_d), sin(θ_d) * sin(φ_d), cos(θ_d)]
+    r_vec = [rval]
+    return sum(
+        c[i] * (
+            bfs[i].a isa AbstractMatrix ?
+            dot(vec(bfs[i].a), r_cart) * dot(vec(bfs[i].b), r_cart) :
+            dot(bfs[i].a, r_vec) * dot(bfs[i].b, r_vec)
+        ) * exp(-dot(r_vec, parent(bfs[i].A) * r_vec) + dot(bfs[i].s, r_vec))
+            for i in eachindex(bfs)
+    )
+end
+
+function normalise(ρ, grid)
+    dr = step(grid)
+    return ρ ./ (sum(ρ) * dr)
+end
+
+r_1s = range(0.01, 12.0, length = 600)
+r_2p = range(0.01, 22.0, length = 600)
+r_3d = range(0.01, 35.0, length = 600)
+
+ρ_ecg_1s   = normalise([r^2 * abs2(ψ₀([r], sr_1s)) for r in r_1s], r_1s)
+ρ_ecg_2p   = normalise([r^2 * abs2(ψ_rank1(r, c_2p, basis_p)) for r in r_2p], r_2p)
+ρ_ecg_rank2 = normalise([r^2 * abs2(ψ_rank2(r, c_rank2, basis_d)) for r in r_3d], r_3d)
+
+θ_p, φ_p = 0.0, 0.0
+θ_d, φ_d = π / 4, 0.0
+
+ρ_exact_1s = normalise(
+    [r^2 * abs2(Antique.ψ(atom, r, 0.0, 0.0; n = 1, l = 0, m = 0)) for r in r_1s],
+    r_1s,
+)
+ρ_exact_2p = normalise(
+    [r^2 * abs2(Antique.ψ(atom, r, θ_p, φ_p; n = 2, l = 1, m = 0)) for r in r_2p],
+    r_2p,
+)
+ρ_exact_3d = normalise(
+    [r^2 * abs2(Antique.ψ(atom, r, θ_d, φ_d; n = 3, l = 2, m = 1)) for r in r_3d],
+    r_3d,
+)
+
+p = plot(
+    layout = (3, 1),
+    size   = (720, 1000),
+    left_margin   = 6Plots.mm,
+    bottom_margin = 4Plots.mm,
+    legend = :topright,
+)
+
+plot!(p[1], collect(r_1s), ρ_ecg_1s;
+    label = "ECG  Rank0 (16 fn.)", lw = 2.5, color = :steelblue)
+plot!(p[1], collect(r_1s), ρ_exact_1s;
+    label = "Antique 1s", lw = 1.8, ls = :dash, color = :black)
+xlabel!(p[1], "r (a.u.)")
+ylabel!(p[1], "r²|ψ(r)|²")
+title!(p[1], "1s  (L=0)  E = $(round(E_1s; digits=6)) Ha  |  exact = $(round(Antique.E(atom; n = 1); digits = 6))")
+
+plot!(p[2], collect(r_2p), ρ_ecg_2p;
+    label = "ECG  Rank1 (16 fn.)", lw = 2.5, color = :tomato)
+plot!(p[2], collect(r_2p), ρ_exact_2p;
+    label = "Antique 2p (θ=0)", lw = 1.8, ls = :dash, color = :black)
+xlabel!(p[2], "r (a.u.)")
+ylabel!(p[2], "r²|ψ(r)|²")
+title!(p[2], "2p  (L=1)  E = $(round(E_2p; digits=6)) Ha  |  exact = $(round(Antique.E(atom; n = 2); digits = 6))")
+
+plot!(p[3], collect(r_3d), ρ_ecg_rank2;
+    label = "ECG  Rank2 (12 fn.)", lw = 2.5, color = :seagreen)
+plot!(p[3], collect(r_3d), ρ_exact_3d;
+    label = "Antique 3d (θ=π/4, m=1)", lw = 1.8, ls = :dash, color = :black)
+xlabel!(p[3], "r (a.u.)")
+ylabel!(p[3], "r²|ψ(r)|²")
+title!(p[3], "pure d-wave  Rank2  E = $(round(E_rank2; digits=6)) Ha  |  exact 3d = $(round(Antique.E(atom; n = 3); digits = 6))")
+
+display(p)

Examples/Hydrogen_p-wave.jl

@@ -2,27 +2,15 @@ using FewBodyECG
 using LinearAlgebra
 using Plots
 
-masses = [1.0e15, 1.0]
+masses = [1e12, 1.0]
 
-Λmat = Λ(masses)
-kin = KineticOperator(Λmat)
+os = Operators(masses)
+os += "Kinetic"
+os += "Coulomb", 1, 2, -1.0   # p-e (attraction)
 
-J, U = _jacobi_transform(masses)
-
-w_raw = [U' * [1, -1]]
-coeffs = [-1.0]
-
-ops = Operator[
-    kin;
-    (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))...
-]
-
-# Polarization vector for p-wave (selects one spatial direction)
 a_vec = [1.0]
 s_zero = [0.0]
 
-# Use a range of Gaussian widths that spans the spatial extent of the 2p orbital.
-# The 2p state is more diffuse than 1s, so we need wider Gaussians (smaller α).
 alphas = [0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0]
 
 basis_fns = GaussianBase[]
@@ -33,7 +21,7 @@ for (i, α) in enumerate(alphas)
 
     basis = BasisSet(basis_fns)
 
-    H = build_hamiltonian_matrix(basis, ops)
+    H = build_hamiltonian_matrix(basis, os)
     S = build_overlap_matrix(basis)
 
     global vals, vecs = solve_generalized_eigenproblem(H, S)

Examples/Positronium.jl

@@ -5,18 +5,26 @@ import FewBodyECG: default_scale, convergence
 
 masses = [1.0, 1.0, 1.0]
 
-Λmat = Λ(masses)
-kin = KineticOperator(Λmat)
+os = Operators(masses, [+1, -1, -1])   # e⁺, e⁻, e⁻
+os += "Kinetic"
+os += "Coulomb"
 
-J, U = _jacobi_transform(masses)
-w_list = [[1, -1, 0], [1, 0, -1], [0, 1, -1]]
-w_raw = [U' * w for w in w_list]
-
-coeffs = [-1.0, -1.0, +1.0]
-coulomb_ops = [CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw)]
-
-ops = Operator[kin; coulomb_ops...]
 scale = default_scale(masses)
 
-result = solve_ECG(ops, 300, sampler = SobolSample(); scale = scale)
+# Ps⁻ has only one bound state; for excited-state examples see Helium.jl.
+result = solve_ECG(os, 300, sampler = SobolSample(); scale = scale, verbose=false, state = 1)
 println("E ≈ ", result.ground_state)
+
+n_conv, E_conv = convergence(result)
+p1 = plot(n_conv, E_conv,
+    xlabel = "Basis size", ylabel = "E (Ha)",
+    label = "Ground state", lw = 2,
+    title = "Positronium convergence")
+display(p1)
+
+r_grid, ρ = correlation_function(result; rmax = 15.0)
+p2 = plot(r_grid, ρ,
+    xlabel = "r (a.u.)", ylabel = "r²|ψ(r)|²",
+    label = "Positronium", lw = 2,
+    title = "Positronium radial correlation")
+display(p2)