diff --git a/docs/Project.toml b/docs/Project.toml
index 3b2bed0..aff1204 100644
--- a/docs/Project.toml
+++ b/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"
diff --git a/docs/make.jl b/docs/make.jl
index 7479707..fe61fc5 100644
--- a/docs/make.jl
+++ b/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
diff --git a/docs/src/FDM.md b/docs/src/FDM.md
index 5248bac..8f9f65f 100644
--- a/docs/src/FDM.md
+++ b/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
diff --git a/docs/src/GEM.md b/docs/src/GEM.md
new file mode 100644
index 0000000..6b9580d
--- /dev/null
+++ b/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
+```
diff --git a/docs/src/Rayleigh-Ritz.md b/docs/src/Rayleigh-Ritz.md
index 6901c04..dd70dbb 100644
--- a/docs/src/Rayleigh-Ritz.md
+++ b/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
diff --git a/docs/src/index.md b/docs/src/index.md
index 57f8f39..e39e0ad 100644
--- a/docs/src/index.md
+++ b/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
 ```
diff --git a/src/Basis.jl b/src/Basis.jl
index f1b4d2d..2d6e97e 100644
--- a/src/Basis.jl
+++ b/src/Basis.jl
@@ -1,4 +1,6 @@
-export BasisSet, Basis, GeometricBasisSet, PrimitiveBasis, ContractedBasis, SimpleGaussianBasis, GaussianBasis, PowerSlaterBasis, FC
+export BasisSet, Basis, GeometricBasisSet, PrimitiveBasis, ContractedBasis,
+       SimpleGaussianBasis, GaussianBasis, ComplexGaussianBasis,
+       ComplexGaussianBasisSet, PowerSlaterBasis, FC
 
 # type
 
@@ -56,6 +58,45 @@ Base.@kwdef struct GaussianBasis{T<:Real} <: PrimitiveBasis
   m::Int = 0
 end
 
+Base.@kwdef struct ComplexGaussianBasis{T<:Real,W<:Real} <: PrimitiveBasis
+  a::T = 1
+  ω::W = 1.0
+  component::Symbol = :cos
+  l::Int = 0
+  m::Int = 0
+end
+
+# GeometricBasisSet constructs primitives from one exponent at a time.
+GaussianBasis(a; l=0, m=0) = GaussianBasis(a, l, m)
+ComplexGaussianBasis(a, ω, component; l=0, m=0) =
+  ComplexGaussianBasis(a, ω, component, l, m)
+
+function ComplexGaussianBasisSet(r₁, rₙ, n::Int; ω=1.0, l=0, m=0)
+  n > 0 || throw(ArgumentError("n must be positive"))
+  primitives = ComplexGaussianBasis[]
+  for a in geometric(r₁, rₙ, n)
+    push!(primitives, ComplexGaussianBasis(a, ω, :cos, l, m))
+    push!(primitives, ComplexGaussianBasis(a, ω, :sin, l, m))
+  end
+  return BasisSet(primitives...)
+end
+
+@doc raw"""
+`ComplexGaussianBasisSet(r₁, rₙ, n; ω=1.0, l=0, m=0)`
+
+Construct the `2n` real, normalized complex-range Gaussian primitives
+
+```math
+r^l e^{-\nu_j r^2}\cos(\omega\nu_jr^2),\qquad
+r^l e^{-\nu_j r^2}\sin(\omega\nu_jr^2),
+```
+
+where `νⱼ = 1/rⱼ²` and the ranges from `r₁` through `rₙ` form a
+geometric progression. This is the real cos/sin form of the complex-range GEM
+basis introduced by Hiyama, Kino, and Kamimura (2003) and used for the
+highly excited hydrogen example of Hiyama and Kamimura (2018).
+""" ComplexGaussianBasisSet
+
 Base.@kwdef struct PowerSlaterBasis <: PrimitiveBasis
   n::Int = 0
   a::Real = 1
@@ -122,7 +163,6 @@ Base.length(::ContractedBasis{N}) where {N} = N
 
 φ(b::SimpleGaussianBasis, r) = exp(-b.a*r^2)
 φ(b::PowerSlaterBasis, r) = r^b.n * exp(-b.a*r)
-φ(b::GaussianBasis, r, θ, φ) = N(b.l) * r^b.l * exp(-b.a*r^2) * Y(b.l, b.m, θ, φ)
 @inline φ(b::ContractedBasis, coordinates...) = _contracted_value(
   getfield(b, :coefficients),
   getfield(b, :primitives),
@@ -201,6 +241,9 @@ function FC(hamiltonian::Hamiltonian, basisset::BasisSet; g::PowerSlaterBasis=Po
   return BasisSet(unique(complements)...)
 end
 
+_replace_exponent(b::ComplexGaussianBasis, a) =
+  ComplexGaussianBasis(a, b.ω, b.component, b.l, b.m)
+
 # function for testing
 
 function expikr(k,θk,φk,r,θr,φr)
@@ -297,14 +340,14 @@ julia> ν = TwoBody.geometric(0.1, 10.0, 5, nₘₐₓ = 10)
 @doc raw"""
 `GeometricBasisSet(basistype, r₁, rₙ, n; nₘᵢₙ=1, nₘₐₓ=n)`
 
-This is a basis set with exponentials generated by `geometric(r₁, rₙ, n; nₘₐₓ=n, nₘᵢₙ=1)`. You can define the same basis set as Table A2 in [E. Hiyama, M. Kamimura, Front. Phys. 13, 132106 (2018)](https://doi.org/10.1007/s11467-018-0828-5) like this:
+This is a basis set with exponents generated by `geometric(r₁, rₙ, n; nₘₐₓ=n, nₘᵢₙ=1)`. You can define the 20 real-range Gaussian basis functions used for the hydrogen atom in Appendix A.2 and Table VII of [E. Hiyama, M. Kamimura, Front. Phys. 13, 132106 (2018)](https://doi.org/10.1007/s11467-018-0828-5) like this:
 ```math
   r_1 = 0.1,
   r_{n_\mathrm{max}} = 80.0,
   n_\mathrm{max} = 20.
 ```
 ```@example
-BS = GeometricBasisSet(SimpleGaussianBasis, 0.1, 80.0, 20)
+BS = GeometricBasisSet(GaussianBasis, 0.1, 80.0, 20)
 ```
 """ GeometricBasisSet
 
@@ -489,10 +532,47 @@ For the hydrogen Hamiltonian, repeated application starting from
 @doc raw"""
 `GaussianBasis(a=1, l=0, m=0)`
 ```math
-\phi_i(r, θ, φ) = N _{il} r^l \exp(-a_i r^2) Y_l^m(θ, φ)
+\phi_{ilm}(r, θ, φ) = N _{il} r^l \exp(-a_i r^2) Y_l^m(θ, φ)
+```
+
+This normalized primitive is used by the Gaussian expansion method (GEM).
+`N_{il}` normalizes the primitive, while `l` and `m` specify its spherical
+harmonic. For a central Hamiltonian, combine primitives with the same `l` and
+`m` and pass the basis set to the existing Rayleigh–Ritz solver.
+
+For example, a ``p``-wave basis can be constructed from geometrically spaced
+exponents as follows:
+
+```julia
+exponents = TwoBody.geometric(0.1, 10.0, 20)
+BS = BasisSet((GaussianBasis(a; l=1, m=0) for a in exponents)...)
 ```
 """ GaussianBasis
 
+@doc raw"""
+`ComplexGaussianBasis(a=1, ω=1, component=:cos, l=0, m=0)`
+
+One normalized real component of a complex-range Gaussian pair:
+
+```math
+\phi^{\cos}_{nlm}=N^{\cos}_{nl}r^l e^{-a_n r^2}
+\cos(\omega a_n r^2)Y_l^m(\hat{\boldsymbol r}),
+```
+
+or the corresponding sine function when `component=:sin`. Equivalently, the
+pair is formed from exponents `(1±iω)aₙ`. Use
+`ComplexGaussianBasisSet(r₁, rₙ, n; ω, l, m)` to construct both members
+at every geometrically spaced range.
+
+The analytic overlap, kinetic, constant, power-law, Coulomb, linear, and
+Gaussian-potential matrix elements are supported by the ordinary
+Rayleigh–Ritz `solve` interface.
+
+References: E. Hiyama, Y. Kino, and M. Kamimura, *Prog. Part. Nucl. Phys.*
+**51**, 223–307 (2003), and E. Hiyama and M. Kamimura, *Front. Phys.* **13**,
+132106 (2018).
+""" ComplexGaussianBasis
+
 @doc raw"""
 `ContractedBasis([c1, c2, ...], [primitive1, primitive2, ...])`
 ```math
diff --git a/src/GEM.jl b/src/GEM.jl
new file mode 100644
index 0000000..c09e04a
--- /dev/null
+++ b/src/GEM.jl
@@ -0,0 +1,370 @@
+export φp, ψp
+
+import LinearAlgebra
+import SpecialFunctions
+
+# Generalized Gauss–Laguerre nodes and weights for x^α exp(-x). The only
+# non-analytic GEM matrix elements reuse these angular-momentum-specific grids.
+const _gausslaguerre_cache = Dict{Tuple{Float64,Int},Tuple{Vector{Float64},Vector{Float64}}}()
+const _gem_quadrature_order = 128
+
+function _gausslaguerre(α::Real, n::Int=_gem_quadrature_order)
+  α > -1 || throw(ArgumentError("Gauss–Laguerre quadrature requires α > -1"))
+  key = (Float64(α), n)
+  return get!(_gausslaguerre_cache, key) do
+    k = collect(1.0:n)
+    diagonal = 2 .* k .- 1 .+ α
+    indices = collect(1.0:n-1)
+    offdiagonal = sqrt.(indices .* (indices .+ α))
+    decomposition = LinearAlgebra.eigen(LinearAlgebra.SymTridiagonal(diagonal, offdiagonal))
+    weights = SpecialFunctions.gamma(α + 1) .* abs2.(decomposition.vectors[1, :])
+    (decomposition.values, weights)
+  end
+end
+
+function _check_gaussian_basis(b::GaussianBasis)
+  b.a > 0 || throw(ArgumentError("Gaussian exponent a must be positive"))
+  b.l isa Integer && b.l ≥ 0 || throw(ArgumentError("l must be a nonnegative integer"))
+  b.m isa Integer && abs(b.m) ≤ b.l || throw(ArgumentError("m must be an integer with |m| ≤ l"))
+  return b
+end
+
+function _check_complex_gaussian_basis(b::ComplexGaussianBasis)
+  b.a > 0 || throw(ArgumentError("Gaussian exponent a must be positive"))
+  b.ω > 0 || throw(ArgumentError("complex-range parameter ω must be positive"))
+  b.component in (:cos, :sin) ||
+    throw(ArgumentError("component must be :cos or :sin"))
+  b.l isa Integer && b.l ≥ 0 || throw(ArgumentError("l must be a nonnegative integer"))
+  b.m isa Integer && abs(b.m) ≤ b.l ||
+    throw(ArgumentError("m must be an integer with |m| ≤ l"))
+  return b
+end
+
+function _gaussian_normalization(a::Real, l::Integer)
+  a > 0 || throw(ArgumentError("Gaussian exponent a must be positive"))
+  l ≥ 0 || throw(ArgumentError("l must be nonnegative"))
+  q = l + 3/2
+  return sqrt(2 * (2 * a)^q / SpecialFunctions.gamma(q))
+end
+
+function _associated_legendre(l::Integer, m::Integer, x::Real)
+  0 ≤ m ≤ l || throw(ArgumentError("associated Legendre indices require 0 ≤ m ≤ l"))
+  pmm = one(float(x))
+  if m > 0
+    root = sqrt(max(zero(float(x)), one(float(x)) - x*x))
+    factor = one(float(x))
+    for _ in 1:m
+      pmm *= -factor * root
+      factor += 2
+    end
+  end
+  l == m && return pmm
+  pmmp1 = x * (2 * m + 1) * pmm
+  l == m + 1 && return pmmp1
+  previous, current = pmm, pmmp1
+  for degree in m+2:l
+    following = ((2 * degree - 1) * x * current - (degree + m - 1) * previous) / (degree - m)
+    previous, current = current, following
+  end
+  return current
+end
+
+function _spherical_harmonic(l::Integer, m::Integer, θ::Real, ϕ::Real)
+  l ≥ 0 || throw(ArgumentError("l must be nonnegative"))
+  abs(m) ≤ l || throw(ArgumentError("m must satisfy |m| ≤ l"))
+  if m < 0
+    return (-1)^(-m) * conj(_spherical_harmonic(l, -m, θ, ϕ))
+  end
+  ratio = SpecialFunctions.gamma(l - m + 1) / SpecialFunctions.gamma(l + m + 1)
+  normalization = sqrt((2 * l + 1) / (4π) * ratio)
+  return normalization * _associated_legendre(l, m, cos(θ)) * exp(im * m * ϕ)
+end
+
+"""Return the normalized radial part of a Gaussian GEM primitive."""
+function φ(b::GaussianBasis, r)
+  _check_gaussian_basis(b)
+  return _gaussian_normalization(b.a, b.l) * r^b.l * exp(-b.a * r^2)
+end
+
+"""Return a normalized Gaussian GEM primitive in position space."""
+φ(b::GaussianBasis, r, θ, ϕ) = φ(b, r) * _spherical_harmonic(b.l, b.m, θ, ϕ)
+
+function _complex_terms(b::ComplexGaussianBasis)
+  _check_complex_gaussian_basis(b)
+  exponents = (b.a * (1 + im * b.ω), b.a * (1 - im * b.ω))
+  coefficients = b.component === :cos ? (1/2, 1/2) : (im/2, -im/2)
+  return zip(exponents, coefficients)
+end
+
+function _complex_pair_sum(f, b1::ComplexGaussianBasis, b2::ComplexGaussianBasis)
+  return sum(
+    conj(c1) * c2 * f(conj(η1) + η2, conj(η1), η2)
+    for (η1, c1) in _complex_terms(b1), (η2, c2) in _complex_terms(b2)
+  )
+end
+
+function _complex_raw_moment(b1::ComplexGaussianBasis, b2::ComplexGaussianBasis, power)
+  q = b1.l + (power + 3) / 2
+  q > 0 || throw(ArgumentError("radial matrix element does not converge"))
+  return _complex_pair_sum((z, _, _) -> SpecialFunctions.gamma(q) / (2 * z^q), b1, b2)
+end
+
+function _complex_gaussian_normalization(b::ComplexGaussianBasis)
+  norm² = real(_complex_raw_moment(b, b, 0))
+  norm² > 0 || throw(ArgumentError("complex-range primitive has zero norm"))
+  return inv(sqrt(norm²))
+end
+
+"""Return the normalized radial part of a real complex-range Gaussian primitive."""
+function φ(b::ComplexGaussianBasis, r)
+  normalization = _complex_gaussian_normalization(b)
+  oscillation = b.component === :cos ? cos(b.ω * b.a * r^2) : sin(b.ω * b.a * r^2)
+  return normalization * r^b.l * exp(-b.a * r^2) * oscillation
+end
+
+φ(b::ComplexGaussianBasis, r, θ, ϕ) =
+  φ(b, r) * _spherical_harmonic(b.l, b.m, θ, ϕ)
+
+"""Return the radial momentum-space Gaussian primitive (unitary Fourier convention)."""
+function φp(b::GaussianBasis, p)
+  _check_gaussian_basis(b)
+  q = b.l + 3/2
+  return im^b.l * _gaussian_normalization(b.a, b.l) * p^b.l * exp(-p^2 / (4 * b.a)) / (2 * b.a)^q
+end
+
+"""Return the full momentum-space Gaussian primitive."""
+φp(b::GaussianBasis, p, θ, ϕ) = φp(b, p) * _spherical_harmonic(b.l, b.m, θ, ϕ)
+
+"""Evaluate a Rayleigh–Ritz eigenfunction in momentum space."""
+function ψp(result::ResultRayleighRitz, p; n::Int=1)
+  return sum(result.C[i, n] * φp(result.basisset[i], p) for i in 1:result.nₘₐₓ)
+end
+
+function ψp(result::ResultRayleighRitz, p, θ, ϕ; n::Int=1)
+  return sum(result.C[i, n] * φp(result.basisset[i], p, θ, ϕ) for i in 1:result.nₘₐₓ)
+end
+
+function ψ(result::ResultRayleighRitz, r, θ, ϕ; n::Int=1)
+  return sum(result.C[i, n] * φ(result.basisset[i], r, θ, ϕ) for i in 1:result.nₘₐₓ)
+end
+
+_same_channel(b1::GaussianBasis, b2::GaussianBasis) = b1.l == b2.l && b1.m == b2.m
+_same_channel(b1::ComplexGaussianBasis, b2::ComplexGaussianBasis) =
+  b1.l == b2.l && b1.m == b2.m
+
+function _normalization_product(b1::GaussianBasis, b2::GaussianBasis)
+  _check_gaussian_basis(b1)
+  _check_gaussian_basis(b2)
+  return _gaussian_normalization(b1.a, b1.l) * _gaussian_normalization(b2.a, b2.l)
+end
+
+function element(b1::GaussianBasis, b2::GaussianBasis)
+  _same_channel(b1, b2) || return 0.0
+  q = b1.l + 3/2
+  return _normalization_product(b1, b2) * SpecialFunctions.gamma(q) / (2 * (b1.a + b2.a)^q)
+end
+
+function _central_element(f, b1::GaussianBasis, b2::GaussianBasis)
+  _same_channel(b1, b2) || return 0.0
+  q = b1.l + 3/2
+  a = b1.a + b2.a
+  nodes, weights = _gausslaguerre(q - 1)
+  integral = sum(weights .* f.(sqrt.(nodes ./ a))) / (2 * a^q)
+  return _normalization_product(b1, b2) * integral
+end
+
+element(o::RestEnergy, b1::GaussianBasis, b2::GaussianBasis) = o.m * o.c^2 * element(b1, b2)
+
+function element(o::Laplacian, b1::GaussianBasis, b2::GaussianBasis)
+  _same_channel(b1, b2) || return 0.0
+  return -2 * o.coefficient * (2 * b1.l + 3) * b1.a * b2.a / (b1.a + b2.a) * element(b1, b2)
+end
+
+function element(o::Kinetic, b1::GaussianBasis, b2::GaussianBasis)
+  _same_channel(b1, b2) || return 0.0
+  return o.hbar^2 / o.m * (2 * b1.l + 3) * b1.a * b2.a / (b1.a + b2.a) * element(b1, b2)
+end
+
+element(o::Constant, b1::GaussianBasis, b2::GaussianBasis) = o.constant * element(b1, b2)
+
+function element(o::PowerLaw, b1::GaussianBasis, b2::GaussianBasis)
+  _same_channel(b1, b2) || return 0.0
+  q = b1.l + (o.exponent + 3) / 2
+  return o.coefficient * _normalization_product(b1, b2) * SpecialFunctions.gamma(q) / (2 * (b1.a + b2.a)^q)
+end
+
+element(o::Linear, b1::GaussianBasis, b2::GaussianBasis) = element(PowerLaw(o.coefficient, 1), b1, b2)
+element(o::Coulomb, b1::GaussianBasis, b2::GaussianBasis) = element(PowerLaw(o.coefficient, -1), b1, b2)
+
+function element(o::Gaussian, b1::GaussianBasis, b2::GaussianBasis)
+  _same_channel(b1, b2) || return 0.0
+  q = b1.l + 3/2
+  return o.coefficient * _normalization_product(b1, b2) * SpecialFunctions.gamma(q) / (2 * (b1.a + b2.a + o.exponent)^q)
+end
+
+function _exp_gaussian_integral(n::Integer, σ::Real, a::Real)
+  n ≥ 0 || throw(ArgumentError("integral power must be nonnegative"))
+  a > 0 || throw(ArgumentError("Gaussian exponent must be positive"))
+  i₀ = sqrt(π) * SpecialFunctions.erfcx(σ / (2 * sqrt(a))) / (2 * sqrt(a))
+  n == 0 && return i₀
+  i₁ = (1 - σ * i₀) / (2 * a)
+  n == 1 && return i₁
+  previous, current = i₀, i₁
+  for power in 2:n
+    following = ((power - 1) * previous - σ * current) / (2 * a)
+    previous, current = current, following
+  end
+  return current
+end
+
+function element(o::Exponential, b1::GaussianBasis, b2::GaussianBasis)
+  _same_channel(b1, b2) || return 0.0
+  radial = _exp_gaussian_integral(2 * b1.l + 2, o.exponent, b1.a + b2.a)
+  return o.coefficient * _normalization_product(b1, b2) * radial
+end
+
+
+function element(o::Yukawa, b1::GaussianBasis, b2::GaussianBasis)
+  _same_channel(b1, b2) || return 0.0
+  radial = _exp_gaussian_integral(2 * b1.l + 1, o.exponent, b1.a + b2.a)
+  return o.coefficient * _normalization_product(b1, b2) * radial
+end
+
+element(o::Custom, b1::GaussianBasis, b2::GaussianBasis) = _central_element(o.f, b1, b2)
+
+function element(o::Delta, b1::GaussianBasis, b2::GaussianBasis)
+  _same_channel(b1, b2) || return 0.0
+  b1.l == 0 || return 0.0
+  return o.coefficient * _normalization_product(b1, b2) / (4π)
+end
+
+function _momentum_moment(n::Integer, b1::GaussianBasis, b2::GaussianBasis)
+  _same_channel(b1, b2) || return 0.0
+  n ≥ 0 || throw(ArgumentError("moment order must be nonnegative"))
+  q = b1.l + 3/2
+  a = (inv(b1.a) + inv(b2.a)) / 4
+  prefactor = _normalization_product(b1, b2) / (4 * b1.a * b2.a)^q
+  return prefactor * SpecialFunctions.gamma(q + n) / (2 * a^(q + n))
+end
+
+function _sqrt_binomial(n::Integer)
+  n ≥ 0 || throw(ArgumentError("expansion order must be nonnegative"))
+  coefficient = 1.0
+  for k in 1:n
+    coefficient *= (3/2 - k) / k
+  end
+  return coefficient
+end
+
+function element(o::RelativisticCorrection, b1::GaussianBasis, b2::GaussianBasis)
+  o.n ≥ 1 || throw(ArgumentError("RelativisticCorrection requires n ≥ 1"))
+  coefficient = _sqrt_binomial(o.n) / (o.m^(2 * o.n - 1) * o.c^(2 * o.n - 2))
+  return coefficient * _momentum_moment(o.n, b1, b2)
+end
+
+function element(o::RelativisticKinetic, b1::GaussianBasis, b2::GaussianBasis)
+  _same_channel(b1, b2) || return 0.0
+  q = b1.l + 3/2
+  a = (inv(b1.a) + inv(b2.a)) / 4
+  nodes, weights = _gausslaguerre(q - 1)
+  p² = nodes ./ a
+  mc² = o.m * o.c^2
+  kinetic = p² .* o.c^2 ./ (sqrt.(o.m^2 * o.c^4 .+ p² .* o.c^2) .+ mc²)
+  radial = sum(weights .* kinetic) / (2 * a^q)
+  prefactor = _normalization_product(b1, b2) / (4 * b1.a * b2.a)^q
+  return prefactor * radial
+end
+
+function _complex_normalization_product(
+  b1::ComplexGaussianBasis,
+  b2::ComplexGaussianBasis,
+)
+  return _complex_gaussian_normalization(b1) * _complex_gaussian_normalization(b2)
+end
+
+function element(b1::ComplexGaussianBasis, b2::ComplexGaussianBasis)
+  _same_channel(b1, b2) || return 0.0
+  value = _complex_normalization_product(b1, b2) * _complex_raw_moment(b1, b2, 0)
+  return real(value)
+end
+
+element(o::RestEnergy, b1::ComplexGaussianBasis, b2::ComplexGaussianBasis) =
+  o.m * o.c^2 * element(b1, b2)
+
+function _complex_kinetic_element(b1::ComplexGaussianBasis, b2::ComplexGaussianBasis)
+  _same_channel(b1, b2) || return 0.0
+  q = b1.l + 3/2
+  raw = _complex_pair_sum(
+    (z, η1, η2) -> (2 * b1.l + 3) * η1 * η2 / z *
+                    SpecialFunctions.gamma(q) / (2 * z^q),
+    b1,
+    b2,
+  )
+  return real(_complex_normalization_product(b1, b2) * raw)
+end
+
+element(o::Laplacian, b1::ComplexGaussianBasis, b2::ComplexGaussianBasis) =
+  -2 * o.coefficient * _complex_kinetic_element(b1, b2)
+
+element(o::Kinetic, b1::ComplexGaussianBasis, b2::ComplexGaussianBasis) =
+  o.hbar^2 / o.m * _complex_kinetic_element(b1, b2)
+
+element(o::Constant, b1::ComplexGaussianBasis, b2::ComplexGaussianBasis) =
+  o.constant * element(b1, b2)
+
+function element(o::PowerLaw, b1::ComplexGaussianBasis, b2::ComplexGaussianBasis)
+  _same_channel(b1, b2) || return 0.0
+  value = o.coefficient * _complex_normalization_product(b1, b2) *
+          _complex_raw_moment(b1, b2, o.exponent)
+  return real(value)
+end
+
+element(o::Linear, b1::ComplexGaussianBasis, b2::ComplexGaussianBasis) =
+  element(PowerLaw(o.coefficient, 1), b1, b2)
+element(o::Coulomb, b1::ComplexGaussianBasis, b2::ComplexGaussianBasis) =
+  element(PowerLaw(o.coefficient, -1), b1, b2)
+
+function element(o::Gaussian, b1::ComplexGaussianBasis, b2::ComplexGaussianBasis)
+  _same_channel(b1, b2) || return 0.0
+  q = b1.l + 3/2
+  raw = _complex_pair_sum(
+    (z, _, _) -> SpecialFunctions.gamma(q) / (2 * (z + o.exponent)^q),
+    b1,
+    b2,
+  )
+  return real(o.coefficient * _complex_normalization_product(b1, b2) * raw)
+end
+
+@doc raw"""
+`φp(b::GaussianBasis, p)` returns the radial momentum-space primitive. With
+the unitary Fourier convention
+
+```math
+\widetilde\phi(\boldsymbol p)=(2\pi)^{-3/2}\int
+e^{i\boldsymbol p\cdot\boldsymbol r}\phi(\boldsymbol r)\,d^3r,
+```
+
+the normalized position-space primitive
+``N_l r^l e^{-a r^2}Y_l^m(\hat{\boldsymbol r})`` becomes
+
+```math
+\widetilde\phi_{alm}(\boldsymbol p)=
+i^l N_l\frac{p^l}{(2a)^{l+3/2}}e^{-p^2/(4a)}
+Y_l^m(\hat{\boldsymbol p}).
+```
+
+The two-argument method returns only its radial factor; the four-argument
+method `φp(b, p, θ, ϕ)` includes the spherical harmonic. Position and momentum
+must use reciprocal units (for example, GeV⁻¹ and GeV when `ℏ=c=1`).
+""" φp
+
+@doc raw"""
+`ψp(result, p; n=1)` evaluates a Gaussian-expansion eigenfunction in momentum
+space by summing `result.C[i,n] * φp(result.basisset[i], p)`. It therefore
+uses the same unitary Fourier convention as `φp`. The four-argument form
+`ψp(result, p, θ, ϕ; n=1)` also includes the spherical harmonic.
+
+This method currently applies to results built from real-range
+`GaussianBasis` primitives. Complex-range cos/sin primitives are intended for
+the position-space Rayleigh–Ritz calculation.
+""" ψp
diff --git a/src/TwoBody.jl b/src/TwoBody.jl
index 8e2607e..b8b5417 100644
--- a/src/TwoBody.jl
+++ b/src/TwoBody.jl
@@ -11,6 +11,7 @@ include("./Basis.jl")
 
 # Solvers
 include("./Rayleigh-Ritz.jl")
+include("./GEM.jl")
 include("./FDM.jl")
 include("./VNN.jl")
 include("./VMC.jl")
diff --git a/test/FDM.jl b/test/FDM.jl
index 9a29164..01b9834 100644
--- a/test/FDM.jl
+++ b/test/FDM.jl
@@ -38,13 +38,13 @@
   end
   
   # comparison with Antique.jl
-  HA = Antique.HydrogenAtom(Z=1, m_e=1.0, a_0=1.0, E_h=1.0, hbar=1.0)
+  HA = Antique.HydrogenAtom(Z=1, mₑ=1.0, a₀=1.0, Eₕ=1.0, ℏ=1.0)
 
   println("Energy")
   println("  i\tnumerical  \tanalytical")
   for i in 1:res.nₘₐₓ
     numerical  = res.E[i]
-    analytical = Antique.energy(HA, n=i)
+    analytical = Antique.energy(HA; n=i)
     error = iszero(analytical) ? abs(numerical-analytical) : abs((numerical-analytical)/analytical)
     acceptance = error < 1e-2
     @printf("%3d\t%.9f\t%.9f\t%s\n", i, numerical, analytical, acceptance ? "✔" :  "✗")
@@ -58,7 +58,7 @@
     for i in keys(res.method.R[begin:min(10,length(res.method.R))])
         r = res.method.R[i]
         numerical  = abs(res.ψ[i,n])
-        analytical = abs(Antique.wavefunction(HA, r, 0, 0, n=n))
+        analytical = abs(Antique.wavefunction(HA, r, 0, 0; n=n))
         error = iszero(analytical) ? abs(numerical-analytical) : abs((numerical-analytical)/analytical)
         acceptance = error < 5e-2
         @printf("%3d\t%.1f\t%.9f\t%.9f\t%s\n", i, r, numerical, analytical, acceptance ? "✔" :  "✗")
diff --git a/test/GEM.jl b/test/GEM.jl
new file mode 100644
index 0000000..1144f26
--- /dev/null
+++ b/test/GEM.jl
@@ -0,0 +1,146 @@
+@testset "GEM.jl" begin
+  @testset "normalized Gaussian basis" begin
+    for l in 0:3
+      basis = GaussianBasis(0.7, l, 0)
+      position_norm = quadgk(r -> r^2 * abs2(TwoBody.φ(basis, r)), 0, Inf, rtol=1e-11)[1]
+      momentum_norm = quadgk(p -> p^2 * abs2(φp(basis, p)), 0, Inf, rtol=1e-11)[1]
+      @test position_norm ≈ 1 atol=1e-10
+      @test momentum_norm ≈ 1 atol=1e-10
+    end
+
+    basis = GaussianBasis(0.7, 0, 0)
+    @test TwoBody.φ(basis, 0.4, 0.3, 0.2) ≈ TwoBody.φ(basis, 0.4) / sqrt(4π)
+    @test TwoBody._replace_exponent(GaussianBasis(0.7, 2, -1), 1.2) == GaussianBasis(1.2, 2, -1)
+    @test_throws ArgumentError TwoBody.φ(GaussianBasis(-1.0), 0.0)
+    @test_throws ArgumentError TwoBody.φ(GaussianBasis(1.0, 1, 2), 0.0)
+  end
+
+  @testset "matrix elements" begin
+    b1 = GaussianBasis(0.6, 1, 0)
+    b2 = GaussianBasis(1.1, 1, 0)
+    overlap = quadgk(r -> r^2 * TwoBody.φ(b1, r) * TwoBody.φ(b2, r), 0, Inf, rtol=1e-11)[1]
+    @test TwoBody.element(b1, b2) ≈ overlap rtol=1e-10
+    @test TwoBody.element(b1, GaussianBasis(1.1, 0, 0)) == 0
+
+    for operator in (
+      Constant(constant=0.3),
+      Linear(coefficient=0.2),
+      Coulomb(coefficient=-0.4),
+      PowerLaw(coefficient=0.7, exponent=2.5),
+      Gaussian(coefficient=-0.3, exponent=0.8),
+      Exponential(coefficient=0.2, exponent=0.5),
+      Yukawa(coefficient=-0.4, exponent=0.3),
+    )
+      numerical = quadgk(
+        r -> r^2 * TwoBody.φ(b1, r) * TwoBody.V(operator, r) * TwoBody.φ(b2, r),
+        0,
+        Inf,
+        rtol=1e-10,
+      )[1]
+      @test TwoBody.element(operator, b1, b2) ≈ numerical rtol=2e-7 atol=1e-10
+    end
+
+    kinetic = Kinetic(hbar=1.3, m=0.8)
+    radial_laplacian(r) = ForwardDiff.derivative(
+      x -> x^2 * ForwardDiff.derivative(y -> TwoBody.φ(b2, y), x),
+      r,
+    ) / r^2 - b2.l * (b2.l + 1) * TwoBody.φ(b2, r) / r^2
+    numerical_kinetic = quadgk(
+      r -> r^2 * TwoBody.φ(b1, r) * (-kinetic.hbar^2 / (2 * kinetic.m)) * radial_laplacian(r),
+      0,
+      Inf,
+      rtol=1e-9,
+    )[1]
+    @test TwoBody.element(kinetic, b1, b2) ≈ numerical_kinetic rtol=2e-8
+
+    correction = RelativisticCorrection(c=2.0, m=1.4, n=2)
+    numerical_correction = quadgk(
+      p -> p^2 * conj(φp(b1, p)) * (-p^4 / (8 * correction.m^3 * correction.c^2)) * φp(b2, p),
+      0,
+      Inf,
+      rtol=1e-10,
+    )[1]
+    @test TwoBody.element(correction, b1, b2) ≈ real(numerical_correction) rtol=2e-9
+
+    relativistic = RelativisticKinetic(c=1.0, m=1.4)
+    numerical_relativistic = quadgk(
+      p -> p^2 * conj(φp(b1, p)) * (sqrt(p^2 + relativistic.m^2) - relativistic.m) * φp(b2, p),
+      0,
+      Inf,
+      rtol=1e-10,
+    )[1]
+    @test TwoBody.element(relativistic, b1, b2) ≈ real(numerical_relativistic) rtol=2e-7
+  end
+
+  @testset "Rayleigh–Ritz integration" begin
+    hydrogen = Hamiltonian(Kinetic(hbar=1, m=1), Coulomb(coefficient=-1))
+    basisset = BasisSet(
+      GaussianBasis(13.00773),
+      GaussianBasis(1.962079),
+      GaussianBasis(0.444529),
+      GaussianBasis(0.1219492),
+    )
+    result = solve(hydrogen, basisset)
+    @test result.E[1] ≈ -0.499278 rtol=1e-5
+    @test result.C[:, 1]' * result.S * result.C[:, 1] ≈ 1 atol=1e-10
+    @test isfinite(ψp(result, 0.5))
+
+    geometric_basis = GeometricBasisSet(GaussianBasis, 0.1, 10.0, 5)
+    @test all(b -> b isa GaussianBasis && b.l == 0 && b.m == 0, geometric_basis.basis)
+  end
+
+  @testset "complex-range basis set" begin
+    basisset = ComplexGaussianBasisSet(0.015, 2000.0, 80; ω=1.5)
+    @test length(basisset) == 160
+
+    for basis in basisset.basis[1:2]
+      norm = quadgk(r -> r^2 * abs2(TwoBody.φ(basis, r)), 0, Inf, rtol=1e-10)[1]
+      @test norm ≈ 1 atol=2e-9
+    end
+  end
+
+  @testset "Hiyama real-range hydrogen benchmark" begin
+    hydrogen = Hamiltonian(Kinetic(hbar=1, m=1), Coulomb(coefficient=-1))
+    basisset = GeometricBasisSet(GaussianBasis, 0.1, 80.0, 20)
+    @test length(basisset) == 20
+
+    result = solve(hydrogen, basisset)
+    hiyama_table_vii = [
+      -0.499982,
+      -0.124998,
+      -0.055555,
+      -0.031249,
+      -0.019998,
+      -0.013883,
+      -0.010203,
+    ]
+    @test maximum(abs.(result.E[1:7] .- hiyama_table_vii)) ≤ 5e-7
+  end
+
+  @testset "Arifi et al. (2024) Hamiltonian" begin
+    ν = 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
+
+    hamiltonian = Hamiltonian(
+      RestEnergy(c=1, m=masses[1]),
+      RelativisticKinetic(c=1, m=masses[1]),
+      RestEnergy(c=1, m=masses[2]),
+      RelativisticKinetic(c=1, m=masses[2]),
+      Constant(constant=a),
+      Linear(coefficient=b),
+      Coulomb(coefficient=-4αs/3),
+      Gaussian(coefficient=hyperfine, exponent=λ^2),
+    )
+    result = solve(hamiltonian, GaussianBasis(ν))
+    # Regression value for the parameter set above.
+    @test result.E[1] ≈ 3.013183414 atol=5e-10
+  end
+end
diff --git a/test/Rayleigh-Ritz.jl b/test/Rayleigh-Ritz.jl
index 43528ff..6839826 100644
--- a/test/Rayleigh-Ritz.jl
+++ b/test/Rayleigh-Ritz.jl
@@ -155,7 +155,7 @@
   end
 
   # comparison with Antique.jl
-  HA = Antique.HydrogenAtom(Z=1, m_e=1.0, a_0=1.0, E_h=1.0, hbar=1.0)
+  HA = Antique.HydrogenAtom(Z=1, mₑ=1.0, a₀=1.0, Eₕ=1.0, ℏ=1.0)
 
   println("ψ(r)")
   println("  r\tnumerical  \tanalytical")
diff --git a/test/runtests.jl b/test/runtests.jl
index 2d84f92..ba7c80e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -12,6 +12,7 @@ using ForwardDiff
   include("Basis.jl")
   include("FC.jl")
   include("Rayleigh-Ritz.jl")
+  include("GEM.jl")
   include("FDM.jl")
   include("VNN.jl")
   include("VMC.jl")