Add Gaussian expansion method

docs/Project.toml

@@ -2,6 +2,7 @@
 Antique = "be6e5d0e-34a5-4c8f-af83-e1b5389203d8"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Lux = "b2108857-7c20-44ae-9111-449ecde12c47"
 Optimisers = "3bd65402-5787-11e9-1adc-39752487f4e2"
 TwoBody = "a92d7657-722c-45a6-9d18-9da4c8a753b6"

docs/make.jl

@@ -3,33 +3,48 @@ using Documenter
 
 DocMeta.setdocmeta!(TwoBody, :DocTestSetup, :(using TwoBody); recursive=true)
 
+local_build = get(ENV, "DOCUMENTER_LOCAL", "false") == "true"
+gem_only = get(ENV, "DOCUMENTER_GEM_ONLY", "false") == "true"
+remote_options = local_build ?
+    (repo = "", remotes = nothing) :
+    (repo = "https://github.com/JuliaFewBody/TwoBody.jl/blob/{commit}{path}#{line}",)
+
+pages = gem_only ?
+    ["Gaussian Expansion Method" => "GEM.md"] :
+    [
+        "Home" => "index.md",
+        "Hamiltonian" => "Hamiltonian.md",
+        "Database" => "DB.md",
+        "Rayleigh-Ritz Method" => "Rayleigh-Ritz.md",
+        "Gaussian Expansion Method" => "GEM.md",
+        "Free Complement Method" => "Free-Complement.md",
+        "Finite Difference Method" => "FDM.md",
+        "Variational Neural Network" => "VNN.md",
+        "Variational Monte Carlo" => "VMC.md",
+        "API reference" => "API.md",
+    ]
+
 makedocs(;
     modules=[TwoBody],
     authors="Shuhei Ohno",
-    repo="https://github.com/JuliaFewBody/TwoBody.jl/blob/{commit}{path}#{line}",
     sitename="TwoBody.jl",
     format=Documenter.HTML(;
         prettyurls=get(ENV, "CI", "false") == "true",
         canonical="https://juliafewbody.github.io/TwoBody.jl/",
+        edit_link=local_build ? nothing : :commit,
         assets=String[
             "./assets/logo.ico",
         ],
     ),
-    pages=[
-        "Home" => "index.md",
-        "Hamiltonian" => "Hamiltonian.md",
-        "Database" => "DB.md",
-        "Rayleigh-Ritz Method" => "Rayleigh-Ritz.md",
-        "Free Complement Method" => "Free-Complement.md",
-        "Finite Difference Method" => "FDM.md",
-        "Variational Neural Network" => "VNN.md",
-        "Variational Monte Carlo" => "VMC.md",
-        "API reference" => "API.md",
-    ],
+    pages=pages,
+    pagesonly=gem_only,
+    checkdocs=gem_only ? :none : :all,
+    remote_options...,
 )
 
 if get(ENV, "CI", "false") == "true"
     deploydocs(;
         repo="github.com/JuliaFewBody/TwoBody.jl",
+        push_preview=true,
     )
 end

docs/src/FDM.md

@@ -85,7 +85,7 @@ res = solve(H, FDM, info=0, nₘₐₓ=4)
 
 # benchmark
 import Antique
-HA = Antique.HydrogenAtom(Z=1, E_h=1.0, a_0=1.0, m_e=1.0, hbar=1.0)
+HA = Antique.HydrogenAtom(Z=1, Eₕ=1.0, a₀=1.0, mₑ=1.0, ℏ=1.0)
 
 # energy
 using Printf
@@ -94,7 +94,7 @@ println("------------------------------")
 println(" n     numerical    analytical")
 println("------------------------------")
 for n in 1:4
-  @printf("%2d  %+.9f  %+.9f\n", n, res.E[n], Antique.energy(HA,n=n))
+  @printf("%2d  %+.9f  %+.9f\n", n, res.E[n], Antique.energy(HA; n=n))
 end
 
 # wave function
@@ -119,7 +119,7 @@ for n in 1:4
   X = res.method.R
   Y = 4π * X .^2 .* res.ψ[:,n] .^ 2
   scatter!(axis, X, Y, label="TwoBody.jl", markersize=6)
-  lines!(axis, 0..50, r -> 4π * r^2 * abs(Antique.wavefunction(HA,r,0,0,n=n))^2, label="Antique.jl", color=:black)
+  lines!(axis, 0..50, r -> 4π * r^2 * abs(Antique.wavefunction(HA, r, 0, 0; n=n))^2, label="Antique.jl", color=:black)
   axislegend(axis, "n = $n", position=:rt, framevisible=false)
 end
 save("assets/FDM_HA.svg", fig) # hide
@@ -140,7 +140,7 @@ res = solve(H, FDM, info=0, nₘₐₓ=4)
 
 # benchmark
 import Antique
-SO = Antique.SphericalOscillator(k=1.0, mu=1.0, hbar=1.0)
+SO = Antique.SphericalOscillator(k=1.0, μ=1.0, ℏ=1.0)
 
 # energy
 using Printf
@@ -149,7 +149,7 @@ println("------------------------------")
 println(" n     numerical    analytical")
 println("------------------------------")
 for n in 1:4
-  @printf("%2d  %+.9f  %+.9f\n", n-1, res.E[n], Antique.energy(SO,n=n-1))
+  @printf("%2d  %+.9f  %+.9f\n", n-1, res.E[n], Antique.energy(SO; n=n-1))
 end
 
 # wave function
@@ -174,7 +174,7 @@ for n in 1:4
   X = res.method.R
   Y = 4π * X .^2 .* res.ψ[:,n] .^ 2
   scatter!(axis, X, Y, label="TwoBody.jl", markersize=6)
-  lines!(axis, 0..50, r -> 4π * r^2 * abs(Antique.wavefunction(SO,r,0,0,n=n-1))^2, label="Antique.jl", color=:black)
+  lines!(axis, 0..50, r -> 4π * r^2 * abs(Antique.wavefunction(SO, r, 0, 0; n=n-1))^2, label="Antique.jl", color=:black)
   axislegend(axis, "n = $(n-1)", position=:rt, framevisible=false)
 end
 fig

docs/src/GEM.md

@@ -0,0 +1,154 @@
+```@meta
+CurrentModule = TwoBody
+```
+
+# Gaussian Expansion Method
+
+The Gaussian expansion method (GEM) uses the existing Rayleigh–Ritz solver
+with normalized Gaussian primitives,
+
+```math
+\psi_{lm}(\boldsymbol r) = \sum_{n=1}^{n_\mathrm{max}} c_n
+N_{nl} r^l e^{-\nu_n r^2}Y_{lm}(\hat{\boldsymbol r}).
+```
+
+The ranges are normally placed in a geometric progression so that one basis
+set covers both short- and long-distance behavior. The implementation follows
+the formulation reviewed by [Hiyama, Kino, and Kamimura
+(2003)](https://www.sciencedirect.com/science/article/abs/pii/S0146641003900159?via%3Dihub)
+and uses TwoBody.jl's ordinary `solve` function for the generalized
+Rayleigh–Ritz eigenvalue problem.
+
+## Usage
+
+GEM uses the same Hamiltonian construction and `solve` interface as the
+[Rayleigh–Ritz method](@ref "Rayleigh-Ritz Method"). Only the basis set is
+changed below.
+
+Run the following code before each use.
+
+```@example gem
+using TwoBody
+```
+
+Define the Hamiltonian. This example uses the non-relativistic hydrogen atom
+in atomic units,
+
+```math
+\hat{H} = -\frac{1}{2}\nabla^2 - \frac{1}{r}.
+```
+
+```@example gem
+H = Hamiltonian(
+  Kinetic(hbar = 1, m = 1),
+  Coulomb(coefficient = -1),
+)
+nothing # hide
+```
+
+Define the Gaussian basis set.
+
+```@example gem
+BS = GeometricBasisSet(GaussianBasis, 0.1, 10.0, 20)
+nothing # hide
+```
+
+Solve the generalized eigenvalue problem with the Rayleigh–Ritz solver.
+
+```@example gem
+result = solve(H, BS)
+result.E[1]
+```
+
+The position- and momentum-space definitions and unit conventions are given
+in the `GaussianBasis`, `φp`, and `ψp` docstrings in the API reference below.
+
+## Example of Hydrogen Atom
+
+Appendix A.2 and Table VII of [Hiyama and Kamimura
+(2018)](https://link.springer.com/article/10.1007/s11467-018-0828-5)
+calculate the lowest seven ``l=0`` states of the hydrogen atom with 20
+real-range Gaussians,
+
+```math
+\phi_n(r) = N_n e^{-\nu_n r^2}.
+```
+
+The Gaussian ranges are placed in a geometric progression with
+``n_{\max}=20``, ``r_1=0.1`` a.u., and ``r_{20}=80`` a.u.
+
+```@example gem
+H = Hamiltonian(
+  Kinetic(hbar = 1, m = 1),
+  Coulomb(coefficient = -1),
+)
+BS = GeometricBasisSet(GaussianBasis, 0.1, 80.0, 20)
+solve(H, BS)
+```
+
+| principal level ``n`` | TwoBody.jl (hartree) | Hiyama Table VII (hartree) | exact (hartree) |
+|:--:|--:|--:|--:|
+| 1 | ``-0.499~982`` | ``-0.499~982`` | ``-0.500~000`` |
+| 2 | ``-0.124~998`` | ``-0.124~998`` | ``-0.125~000`` |
+| 3 | ``-0.055~555`` | ``-0.055~555`` | ``-0.055~556`` |
+| 4 | ``-0.031~249`` | ``-0.031~249`` | ``-0.031~250`` |
+| 5 | ``-0.019~998`` | ``-0.019~998`` | ``-0.020~000`` |
+| 6 | ``-0.013~883`` | ``-0.013~883`` | ``-0.013~889`` |
+| 7 | ``-0.010~203`` | ``-0.010~203`` | ``-0.010~204`` |
+
+The TwoBody.jl results agree with all seven values in Table VII at the six
+decimal places reported there. The complex-range hydrogen calculation in
+Appendix A.6.2 is a separate example for highly excited states and is not the
+calculation reproduced here.
+
+## Example of Charmonium
+
+The following calculation uses one Gaussian with ``\nu=0.2443`` and the SGA
+parameters of [Arifi et al.
+(2024)](https://arxiv.org/abs/2401.07933). Natural units are used, so energies
+and masses are in GeV and lengths are in GeV``^{-1}``.
+
+```@example gem
+ν = 0.2443
+masses = (1.6324, 1.6324)
+a = -0.4235
+b = 0.1655
+αs = 0.4410
+Λ = 0.9639
+spin = -3/4
+
+reduced_mass = inv(inv(masses[1]) + inv(masses[2]))
+λ = Λ * sqrt(reduced_mass)
+hyperfine = 32π * αs * (λ / sqrt(π))^3 /
+            (9 * masses[1] * masses[2]) * spin
+
+H = Hamiltonian(
+  RestEnergy(m=masses[1]), RelativisticKinetic(m=masses[1]),
+  RestEnergy(m=masses[2]), RelativisticKinetic(m=masses[2]),
+  Constant(constant=a),
+  Linear(coefficient=b),
+  Coulomb(coefficient=-4αs/3),
+  Gaussian(coefficient=hyperfine, exponent=λ^2),
+)
+
+eta_c = solve(H, GaussianBasis(ν))
+round(eta_c.E[1] * 1000; digits=6)
+```
+
+| calculation | ``\eta_c`` mass (MeV) |
+|:--|--:|
+| TwoBody.jl with the parameters above | 3013.183414 |
+| Arifi et al. SGA, Table 2 | 3012 |
+
+The 1.18 MeV difference from the paper table is consistent with using the
+rounded parameters and rounded variational exponent shown above.
+
+## API reference
+
+```@docs; canonical=false
+TwoBody.GaussianBasis
+TwoBody.ComplexGaussianBasis
+TwoBody.ComplexGaussianBasisSet
+TwoBody.φp
+TwoBody.ψp
+```

docs/src/Rayleigh-Ritz.md

@@ -159,7 +159,7 @@ res = solve(H, BS)
 
 # benchmark
 import Antique
-HA = Antique.HydrogenAtom(Z=1, E_h=1.0, a_0=1.0, m_e=1.0, hbar=1.0)
+HA = Antique.HydrogenAtom(Z=1, Eₕ=1.0, a₀=1.0, mₑ=1.0, ℏ=1.0)
 
 # energy
 using Printf
@@ -168,7 +168,7 @@ println("------------------------------")
 println(" n     numerical    analytical")
 println("------------------------------")
 for n in 1:4
-  @printf("%2d  %+.9f  %+.9f\n", n, res.E[n], Antique.energy(HA,n=n))
+  @printf("%2d  %+.9f  %+.9f\n", n, res.E[n], Antique.energy(HA; n=n))
 end
 
 # wave function
@@ -191,7 +191,7 @@ for n in 1:4
     )
   )
   lines!(axis, 0..50, r -> 4π * r^2 * abs(TwoBody.ψ(res,r,n=n))^2, label="TwoBody.jl")
-  lines!(axis, 0..50, r -> 4π * r^2 * abs(Antique.wavefunction(HA,r,0,0,n=n))^2, label="Antique.jl", color=:black, linestyle=:dash)
+  lines!(axis, 0..50, r -> 4π * r^2 * abs(Antique.wavefunction(HA, r, 0, 0; n=n))^2, label="Antique.jl", color=:black, linestyle=:dash)
   axislegend(axis, "n = $n", position=:rt, framevisible=false)
 end
 fig
@@ -213,7 +213,7 @@ res = solve(H, BS)
 
 # benchmark
 import Antique
-SO = Antique.SphericalOscillator(k=1.0, mu=1.0, hbar=1.0)
+SO = Antique.SphericalOscillator(k=1.0, μ=1.0, ℏ=1.0)
 
 # energy
 using Printf
@@ -222,7 +222,7 @@ println("------------------------------")
 println(" n     numerical    analytical")
 println("------------------------------")
 for n in 1:4
-  @printf("%2d  %+.9f  %+.9f\n", n-1, res.E[n], Antique.energy(SO,n=n-1))
+  @printf("%2d  %+.9f  %+.9f\n", n-1, res.E[n], Antique.energy(SO; n=n-1))
 end
 
 # wave function
@@ -245,7 +245,7 @@ for n in 1:4
     )
   )
   lines!(axis, 0..50, r -> 4π * r^2 * abs(TwoBody.ψ(res,r,n=n))^2, label="TwoBody.jl")
-  lines!(axis, 0..50, r -> 4π * r^2 * abs(Antique.wavefunction(SO,r,0,0,n=n-1))^2, label="Antique.jl", color=:black, linestyle=:dash)
+  lines!(axis, 0..50, r -> 4π * r^2 * abs(Antique.wavefunction(SO, r, 0, 0; n=n-1))^2, label="Antique.jl", color=:black, linestyle=:dash)
   axislegend(axis, "n = $(n-1)", position=:rt, framevisible=false)
 end
 fig

docs/src/index.md

@@ -74,14 +74,14 @@ res = solve(H, BS)
 
 # benchmark
 import Antique
-HA = Antique.HydrogenAtom(Z=1, E_h=1.0, a_0=1.0, m_e=1.0, hbar=1.0)
+HA = Antique.HydrogenAtom(Z=1, Eₕ=1.0, a₀=1.0, mₑ=1.0, ℏ=1.0)
 
 # plot
 using CairoMakie
 fig = Figure(size=(420,300), fontsize=11, backgroundcolor=:transparent)
 axis = Axis(fig[1,1], xlabel=L"$r / a_0$", ylabel=L"$\psi(r) / a_0^{-3/2}$", ylabelsize=16.5, xlabelsize=16.5, limits=(0,4,0,1.1/sqrt(π)))
 lines!(axis, 0..5, r -> abs(TwoBody.ψ(res,r)), label="TwoBody.jl")
-lines!(axis, 0..5, r -> abs(Antique.wavefunction(HA,r,0,0)), linestyle=:dash, color=:black, label="Antique.jl")
+lines!(axis, 0..5, r -> abs(Antique.wavefunction(HA, r, 0, 0)), linestyle=:dash, color=:black, label="Antique.jl")
 axislegend(axis, position=:rt, framevisible=false)
 fig
 ```