added spin-orbit example

docs/make.jl

@@ -23,8 +23,13 @@ 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",
         ],
         "API" => "API.md",

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

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

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)
 ````
 
 ---

docs/src/examples/harmonium.md

@@ -0,0 +1,107 @@
+```@meta
+EditURL = "../../../examples/harmonium.jl"
+```
+
+````@example harmonium
+using FewBodyECG
+using Plots
+
+ω = 0.5
+masses = [1.0e15, 1.0, 1.0]
+
+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
+````
+
+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...,
+)
+````
+
+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
+
+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,
+)
+````
+
+---
+
+*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
+

docs/src/examples/helium.md

@@ -16,15 +16,15 @@ 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])
 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)

docs/src/examples/muonic_ions.md

@@ -0,0 +1,52 @@
+```@meta
+EditURL = "../../../examples/muonic_ions.jl"
+```
+
+# 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.
+
+````@example muonic_ions
+using FewBodyECG
+
+mμ, md, mt = 206.7686, 3670.481, 5496.918
+
+systems = [
+    (
+        "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, svm200, best_ref) in systems
+    ops = Operators(masses, [+1.0, +1.0, -1.0])
+    ops += "Kinetic"
+    ops += "Coulomb"
+
+    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
+````
+
+---
+
+*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
+

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).*
+