diff --git a/.github/workflows/runic.yml b/.github/workflows/runic.yml
index 4f0136d..1c754c0 100644
--- a/.github/workflows/runic.yml
+++ b/.github/workflows/runic.yml
@@ -12,6 +12,12 @@ jobs:
   runic:
     name: Runic formatting
     runs-on: ubuntu-latest
+    # Permissions needed for reviewdog/action-suggester to post comments
+    permissions:
+      contents: read
+      checks: write
+      issues: write
+      pull-requests: write
     steps:
       - uses: actions/checkout@v4
       # - uses: julia-actions/setup-julia@v2
@@ -20,4 +26,12 @@ jobs:
       # - uses: julia-actions/cache@v2
       - uses: fredrikekre/runic-action@v1
         with:
-          version: '1'
\ No newline at end of file
+          version: '1'
+          format_files: true
+        # Fail on next step instead
+        continue-on-error: ${{ github.event_name == 'pull_request' }}
+      - uses: reviewdog/action-suggester@v1
+        if: github.event_name == 'pull_request'
+        with:
+          tool_name: Runic
+          fail_level: warning
\ No newline at end of file
diff --git a/Examples/HydrogenAnion.jl b/Examples/HydrogenAnion.jl
index 4f0b62a..350f56e 100644
--- a/Examples/HydrogenAnion.jl
+++ b/Examples/HydrogenAnion.jl
@@ -3,26 +3,26 @@ using LinearAlgebra
 using Plots
 using QuasiMonteCarlo
 
-masses = [1e15, 1.0, 1.0] 
+masses = [1.0e15, 1.0, 1.0]
 psys = ParticleSystem(masses)
 
-K = Diagonal([0.0, 1/2, 1/2])
+K = Diagonal([0.0, 1 / 2, 1 / 2])
 K_transformed = psys.J * K * psys.J'
 
-w_list = [ [1, -1, 0], [1, 0, -1], [0, 1, -1] ]
+w_list = [[1, -1, 0], [1, 0, -1], [0, 1, -1]]
 
 w_raw = [psys.U' * w for w in w_list]
 
 let
     n_basis = 50
     b1 = default_b0(psys.scale)
-    method = :quasirandom  
+    method = :quasirandom
     basis_fns = GaussianBase[]
     E₀_list = Float64[]
     coeffs = [-1.0, -1.0, +1.0]
 
     for i in 1:n_basis
-        bij = generate_bij(method, i, length(w_raw), b1; qmc_sampler=SobolSample())
+        bij = generate_bij(method, i, length(w_raw), b1; qmc_sampler = SobolSample())
         A = generate_A_matrix(bij, w_raw)
         push!(basis_fns, Rank0Gaussian(A))
 
@@ -31,7 +31,7 @@ let
             KineticEnergy(K_transformed);
             (CoulombPotential(c, w) for (c, w) in zip(coeffs, w_raw))...
         ]
-                
+
         H = build_hamiltonian_matrix(basis, ops)
         S = build_overlap_matrix(basis)
         vals, _ = solve_generalized_eigenproblem(H, S)
@@ -47,8 +47,9 @@ let
     ΔE = abs(E₀ - Eᵗʰ)
     @show ΔE
 
-    plot(1:n_basis, E₀_list, xlabel="Number of Gaussians", ylabel="E₀ [Hartree]",
-         lw=2, label="Ground state energy", title="Hydrogen Anion Convergence")
+    plot(
+        1:n_basis, E₀_list, xlabel = "Number of Gaussians", ylabel = "E₀ [Hartree]",
+        lw = 2, label = "Ground state energy", title = "Hydrogen Anion Convergence"
+    )
 
 end
-
diff --git a/Examples/Hydrogen_p-wave.jl b/Examples/Hydrogen_p-wave.jl
index 870be01..54a2293 100644
--- a/Examples/Hydrogen_p-wave.jl
+++ b/Examples/Hydrogen_p-wave.jl
@@ -1,61 +1,66 @@
 using FewBodyECG
 using LinearAlgebra
 using Plots
-using QuasiMonteCarlo
 
-masses = [1e15, 1.0]
-psys = ParticleSystem(masses)
+particles = [Particle(1.0e15, 1.0, nothing), Particle(1.0, -1.0, nothing)]
+sys = System(particles, true)
 
-K = Diagonal([0.0, 0.5])
-K_transformed = psys.J * K * psys.J'
+ops = Operator[
+    Kinetic(1, 1.0, 1.0),
+    Coulomb(1, 2, -1.0),
+]
 
-w_raw = [psys.U' * [1, -1]]
-coeffs = [-1.0]
+overlap(bra::Rank0Gaussian, ket::Rank0Gaussian) = begin
+    A, B = bra.A, ket.A
+    R = inv(A + B)
+    (π^size(R, 1) / det(A + B))^(3 / 2)
+end
 
 n_basis = 25
-method = :quasirandom  
-b1 = 1.4
-
-basis_fns = GaussianBase[]
-E₀_list = Float64[]
+basis_fns = Vector{GaussianBase}(undef, 0)
+E0_list = Float64[]
 
-a_vec = [1.0]  
-    
 for i in 1:n_basis
-    bij = generate_bij(method, i, length(w_raw), b1; qmc_sampler=HaltonSample())
-    A = generate_A_matrix(bij, w_raw)
-    push!(basis_fns, Rank1Gaussian(A, [a_vec]))
-
-    basis = BasisSet(basis_fns)
-    ops = Operator[
-        KineticEnergy(K_transformed);
-        (CoulombPotential(c, w) for (c, w) in zip(coeffs, w_raw))...
-    ]
-
-    H = build_hamiltonian_matrix(basis, ops)
-    S = build_overlap_matrix(basis)
-
-    vals, _ = solve_generalized_eigenproblem(H, S)
-    valid = vals .> 1e-12
-    S⁻¹₂ = Diagonal(1 ./ sqrt.(vals[valid]))
-    H̃ = S⁻¹₂ * H[valid, valid] * S⁻¹₂
-    eigvals = eigen(H̃).values
-
-    real_eigvals = eigvals[abs.(imag.(eigvals)) .< 1e-10]
-    if !isempty(real_eigvals)
-        E₀ = minimum(real_eigvals)
-    else
-        E₀ = NaN
-        @warn "All eigenvalues are complex at step $i"
+    dim = 1
+    α = 0.2 + 0.1 * i
+    A = α * I(dim)
+    push!(basis_fns, Rank0Gaussian(A))
+
+    basis = ECGBasis(basis_fns)
+
+    N = length(basis.functions)
+    H = zeros(Float64, N, N)
+    S = zeros(Float64, N, N)
+
+    for p in 1:N, q in 1:N
+        bra = basis.functions[p]
+        ket = basis.functions[q]
+        S[p, q] = overlap(bra, ket)
+
+        for op in ops
+            if op isa Kinetic
+                H[p, q] += compute_matrix_element(bra, ket, op, K)
+            elseif op isa Coulomb
+                H[p, q] += compute_matrix_element(bra, ket, op, w)
+            end
+        end
     end
-    push!(E₀_list, E₀)
-    println("Step $i: E₀ = $E₀")
+
+
+    F = eigen(H, S)
+    vals = real(F.values)
+    E0 = minimum(vals)
+    push!(E0_list, E0)
+    @info "step $i" E0
 end
 
-E_exact = -0.125  #
-E_min = minimum(E₀_list)
-@show ΔE = abs(E_min - E_exact)
+E_exact = -0.125
+E_min = minimum(E0_list)
+ΔE = abs(E_min - E_exact)
+@show ΔE
 
-plot(1:n_basis, E₀_list, xlabel="Number of Gaussians", ylabel="E₀ [Hartree]",
-     lw=2, label="E₀ estimate", title="p-wave Hydrogen Convergence")
-hline!([E_exact], label="Exact: -0.125", linestyle=:dash)
+plot(
+    1:n_basis, E0_list, xlabel = "Number of Gaussians", ylabel = "E₀ [Hartree]",
+    lw = 2, label = "E₀ estimate", title = "S-wave Hydrogen (minimal demo)"
+)
+hline!([E_exact], label = "Exact: -0.125", linestyle = :dash)
diff --git a/Examples/Hydrogen_s-wave.jl b/Examples/Hydrogen_s-wave.jl
index 9b0d742..4eec94e 100644
--- a/Examples/Hydrogen_s-wave.jl
+++ b/Examples/Hydrogen_s-wave.jl
@@ -3,7 +3,7 @@ using LinearAlgebra
 using Plots
 using QuasiMonteCarlo
 
-masses = [1e15, 1.0]  # proton, electron
+masses = [1.0e15, 1.0]  # proton, electron
 psys = ParticleSystem(masses)
 
 K = Diagonal([0.0, 0.5])
@@ -13,14 +13,14 @@ w_raw = [psys.U' * [1, -1]]  # r₁ - r₂
 coeffs = [-1.0]  # Coulomb attraction
 
 n_basis = 25
-method = :quasirandom 
+method = :quasirandom
 b1 = 1.5
 
 basis_fns = GaussianBase[]
 E₀_list = Float64[]
 
 for i in 1:n_basis
-    bij = generate_bij(method, i, length(w_raw), b1; qmc_sampler=SobolSample())
+    bij = generate_bij(method, i, length(w_raw), b1; qmc_sampler = SobolSample())
     A = generate_A_matrix(bij, w_raw)
     push!(basis_fns, Rank0Gaussian(A))
 
@@ -34,20 +34,22 @@ for i in 1:n_basis
     S = build_overlap_matrix(basis)
 
     λs, Us = eigen(S)
-    keep = λs .> 1e-10
+    keep = λs .> 1.0e-10
     S⁻¹₂ = Us[:, keep] * Diagonal(1 ./ sqrt.(λs[keep])) * Us[:, keep]'
     H̃ = Symmetric(S⁻¹₂ * H * S⁻¹₂)
     E₀ = minimum(eigen(H̃).values)
-    
+
     push!(E₀_list, E₀)
-    
-    
+
+
 end
 
 E_exact = -0.5
 E_min = minimum(E₀_list)
 @show ΔE = abs(E_min - E_exact)
 
-plot(1:n_basis, E₀_list, xlabel="Number of Gaussians", ylabel="E₀ [Hartree]",
-     lw=2, label="E₀ estimate", title="s-wave Hydrogen Convergence")
-hline!([E_exact], label="Exact: -0.5", linestyle=:dash)
\ No newline at end of file
+plot(
+    1:n_basis, E₀_list, xlabel = "Number of Gaussians", ylabel = "E₀ [Hartree]",
+    lw = 2, label = "E₀ estimate", title = "s-wave Hydrogen Convergence"
+)
+hline!([E_exact], label = "Exact: -0.5", linestyle = :dash)
diff --git a/Examples/Positronium.jl b/Examples/Positronium.jl
index c532602..daa6f93 100644
--- a/Examples/Positronium.jl
+++ b/Examples/Positronium.jl
@@ -6,7 +6,7 @@ using QuasiMonteCarlo
 masses = [1.0, 1.0, 1.0]
 psys = ParticleSystem(masses)
 
-K = Diagonal([1/2, 1/2, 1/2])
+K = Diagonal([1 / 2, 1 / 2, 1 / 2])
 K_transformed = psys.J * K * psys.J'
 
 w_list = [[1, -1, 0], [1, 0, -1], [0, 1, -1]]
@@ -22,7 +22,7 @@ let
     E₀_list = Float64[]
 
     for i in 1:n_basis
-        bij = generate_bij(method, i, length(w_raw), b1; qmc_sampler=SobolSample())
+        bij = generate_bij(method, i, length(w_raw), b1; qmc_sampler = SobolSample())
         A = generate_A_matrix(bij, w_raw)
         push!(basis_fns, Rank0Gaussian(A))
 
@@ -47,17 +47,20 @@ let
     ΔE = abs(E₀ - Eᵗʰ)
     @show ΔE
 
-    r = range(0.01, 14.0, length=400)
+    r = range(0.01, 14.0, length = 400)
     ρ_r = [rval^2 * abs2(ψ₀([rval, 0.0], c₀, basis_fns)) for rval in r]
-    
-    p1 = plot(r, ρ_r, xlabel="r (a.u.)", ylabel="r²|ψ₀(r)|²",
-              lw=2, label="r²C(r)", title="Electron-Positron Correlation Function")
-    
 
-    p2 = plot(1:n_basis, E₀_list, xlabel="Number of Gaussians", ylabel="E₀ [Hartree]",
-    lw=2, label="Ground state energy", title="Positronium Convergence")
+    p1 = plot(
+        r, ρ_r, xlabel = "r (a.u.)", ylabel = "r²|ψ₀(r)|²",
+        lw = 2, label = "r²C(r)", title = "Electron-Positron Correlation Function"
+    )
 
-    plot(p1, p2, layout=(2, 1))
 
-end
+    p2 = plot(
+        1:n_basis, E₀_list, xlabel = "Number of Gaussians", ylabel = "E₀ [Hartree]",
+        lw = 2, label = "Ground state energy", title = "Positronium Convergence"
+    )
+
+    plot(p1, p2, layout = (2, 1))
 
+end
diff --git a/Manifest.toml b/Manifest.toml
index c30d4d5..5b1aa9e 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,22 +2,7 @@
 
 julia_version = "1.11.6"
 manifest_format = "2.0"
-project_hash = "635cf368533e28beaba4d8c3c9a432477d04656c"
-
-[[deps.ADTypes]]
-git-tree-sha1 = "60665b326b75db6517939d0e1875850bc4a54368"
-uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
-version = "1.17.0"
-
-    [deps.ADTypes.extensions]
-    ADTypesChainRulesCoreExt = "ChainRulesCore"
-    ADTypesConstructionBaseExt = "ConstructionBase"
-    ADTypesEnzymeCoreExt = "EnzymeCore"
-
-    [deps.ADTypes.weakdeps]
-    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-    ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
-    EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869"
+project_hash = "cb3b6f52648369280bb1a7256f54569c38c07503"
 
 [[deps.Accessors]]
 deps = ["CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "MacroTools"]
@@ -43,80 +28,16 @@ version = "0.1.42"
     Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
     Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
-[[deps.Adapt]]
-deps = ["LinearAlgebra", "Requires"]
-git-tree-sha1 = "f7817e2e585aa6d924fd714df1e2a84be7896c60"
-uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "4.3.0"
-
-    [deps.Adapt.extensions]
-    AdaptSparseArraysExt = "SparseArrays"
-    AdaptStaticArraysExt = "StaticArrays"
-
-    [deps.Adapt.weakdeps]
-    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-
 [[deps.AliasTables]]
 deps = ["PtrArrays", "Random"]
 git-tree-sha1 = "9876e1e164b144ca45e9e3198d0b689cadfed9ff"
 uuid = "66dad0bd-aa9a-41b7-9441-69ab47430ed8"
 version = "1.1.3"
 
-[[deps.ArrayInterface]]
-deps = ["Adapt", "LinearAlgebra"]
-git-tree-sha1 = "dbd8c3bbbdbb5c2778f85f4422c39960eac65a42"
-uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.20.0"
-
-    [deps.ArrayInterface.extensions]
-    ArrayInterfaceBandedMatricesExt = "BandedMatrices"
-    ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices"
-    ArrayInterfaceCUDAExt = "CUDA"
-    ArrayInterfaceCUDSSExt = "CUDSS"
-    ArrayInterfaceChainRulesCoreExt = "ChainRulesCore"
-    ArrayInterfaceChainRulesExt = "ChainRules"
-    ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
-    ArrayInterfaceMetalExt = "Metal"
-    ArrayInterfaceReverseDiffExt = "ReverseDiff"
-    ArrayInterfaceSparseArraysExt = "SparseArrays"
-    ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore"
-    ArrayInterfaceTrackerExt = "Tracker"
-
-    [deps.ArrayInterface.weakdeps]
-    BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
-    BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
-    CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
-    CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e"
-    ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
-    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-    GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
-    Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
-    ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
-    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-    StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
-    Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
-
 [[deps.Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 version = "1.11.0"
 
-[[deps.CommonSubexpressions]]
-deps = ["MacroTools"]
-git-tree-sha1 = "cda2cfaebb4be89c9084adaca7dd7333369715c5"
-uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
-version = "0.3.1"
-
-[[deps.Compat]]
-deps = ["TOML", "UUIDs"]
-git-tree-sha1 = "0037835448781bb46feb39866934e243886d756a"
-uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.18.0"
-weakdeps = ["Dates", "LinearAlgebra"]
-
-    [deps.Compat.extensions]
-    CompatLinearAlgebraExt = "LinearAlgebra"
-
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
@@ -173,133 +94,17 @@ git-tree-sha1 = "9e2f36d3c96a820c678f2f1f1782582fcf685bae"
 uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 version = "1.9.1"
 
-[[deps.DiffResults]]
-deps = ["StaticArraysCore"]
-git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621"
-uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
-version = "1.1.0"
-
-[[deps.DiffRules]]
-deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
-git-tree-sha1 = "23163d55f885173722d1e4cf0f6110cdbaf7e272"
-uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
-version = "1.15.1"
-
-[[deps.DifferentiationInterface]]
-deps = ["ADTypes", "LinearAlgebra"]
-git-tree-sha1 = "16946a4d305607c3a4af54ff35d56f0e9444ed0e"
-uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
-version = "0.7.7"
-
-    [deps.DifferentiationInterface.extensions]
-    DifferentiationInterfaceChainRulesCoreExt = "ChainRulesCore"
-    DifferentiationInterfaceDiffractorExt = "Diffractor"
-    DifferentiationInterfaceEnzymeExt = ["EnzymeCore", "Enzyme"]
-    DifferentiationInterfaceFastDifferentiationExt = "FastDifferentiation"
-    DifferentiationInterfaceFiniteDiffExt = "FiniteDiff"
-    DifferentiationInterfaceFiniteDifferencesExt = "FiniteDifferences"
-    DifferentiationInterfaceForwardDiffExt = ["ForwardDiff", "DiffResults"]
-    DifferentiationInterfaceGPUArraysCoreExt = "GPUArraysCore"
-    DifferentiationInterfaceGTPSAExt = "GTPSA"
-    DifferentiationInterfaceMooncakeExt = "Mooncake"
-    DifferentiationInterfacePolyesterForwardDiffExt = ["PolyesterForwardDiff", "ForwardDiff", "DiffResults"]
-    DifferentiationInterfaceReverseDiffExt = ["ReverseDiff", "DiffResults"]
-    DifferentiationInterfaceSparseArraysExt = "SparseArrays"
-    DifferentiationInterfaceSparseConnectivityTracerExt = "SparseConnectivityTracer"
-    DifferentiationInterfaceSparseMatrixColoringsExt = "SparseMatrixColorings"
-    DifferentiationInterfaceStaticArraysExt = "StaticArrays"
-    DifferentiationInterfaceSymbolicsExt = "Symbolics"
-    DifferentiationInterfaceTrackerExt = "Tracker"
-    DifferentiationInterfaceZygoteExt = ["Zygote", "ForwardDiff"]
-
-    [deps.DifferentiationInterface.weakdeps]
-    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-    DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
-    Diffractor = "9f5e2b26-1114-432f-b630-d3fe2085c51c"
-    Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
-    EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869"
-    FastDifferentiation = "eb9bf01b-bf85-4b60-bf87-ee5de06c00be"
-    FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
-    FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
-    ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-    GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
-    GTPSA = "b27dd330-f138-47c5-815b-40db9dd9b6e8"
-    Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6"
-    PolyesterForwardDiff = "98d1487c-24ca-40b6-b7ab-df2af84e126b"
-    ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
-    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-    SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
-    SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
-    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-    Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
-    Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
-    Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
-
-[[deps.Distributed]]
-deps = ["Random", "Serialization", "Sockets"]
-uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
-version = "1.11.0"
-
 [[deps.DocStringExtensions]]
 git-tree-sha1 = "7442a5dfe1ebb773c29cc2962a8980f47221d76c"
 uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 version = "0.9.5"
 
-[[deps.EnumX]]
-git-tree-sha1 = "bddad79635af6aec424f53ed8aad5d7555dc6f00"
-uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
-version = "1.0.5"
-
-[[deps.FillArrays]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "6a70198746448456524cb442b8af316927ff3e1a"
-uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "1.13.0"
-
-    [deps.FillArrays.extensions]
-    FillArraysPDMatsExt = "PDMats"
-    FillArraysSparseArraysExt = "SparseArrays"
-    FillArraysStatisticsExt = "Statistics"
-
-    [deps.FillArrays.weakdeps]
-    PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
-    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-    Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-
-[[deps.FiniteDiff]]
-deps = ["ArrayInterface", "LinearAlgebra", "Setfield"]
-git-tree-sha1 = "31fd32af86234b6b71add76229d53129aa1b87a9"
-uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
-version = "2.28.1"
-
-    [deps.FiniteDiff.extensions]
-    FiniteDiffBandedMatricesExt = "BandedMatrices"
-    FiniteDiffBlockBandedMatricesExt = "BlockBandedMatrices"
-    FiniteDiffSparseArraysExt = "SparseArrays"
-    FiniteDiffStaticArraysExt = "StaticArrays"
-
-    [deps.FiniteDiff.weakdeps]
-    BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
-    BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
-    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-
-[[deps.ForwardDiff]]
-deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"]
-git-tree-sha1 = "ce15956960057e9ff7f1f535400ffa14c92429a4"
-uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "1.1.0"
-
-    [deps.ForwardDiff.extensions]
-    ForwardDiffStaticArraysExt = "StaticArrays"
-
-    [deps.ForwardDiff.weakdeps]
-    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-
-[[deps.Future]]
-deps = ["Random"]
-uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
-version = "1.11.0"
+[[deps.FewBodyHamiltonians]]
+git-tree-sha1 = "8d67c1859d3682d75b2419633cc94edabb9dbbca"
+repo-rev = "main"
+repo-url = "https://github.com/JuliaFewBody/FewBodyHamiltonians.jl"
+uuid = "3a126c26-e5d7-4a95-83c3-3b69f8a11ded"
+version = "0.0.1"
 
 [[deps.IntegerMathUtils]]
 git-tree-sha1 = "4c1acff2dc6b6967e7e750633c50bc3b8d83e617"
@@ -324,12 +129,6 @@ git-tree-sha1 = "e2222959fbc6c19554dc15174c81bf7bf3aa691c"
 uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
 version = "0.2.4"
 
-[[deps.JLLWrappers]]
-deps = ["Artifacts", "Preferences"]
-git-tree-sha1 = "0533e564aae234aff59ab625543145446d8b6ec2"
-uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.7.1"
-
 [[deps.LatticeRules]]
 deps = ["Random"]
 git-tree-sha1 = "7f5b02258a3ca0221a6a9710b0a0a2e8fb4957fe"
@@ -340,12 +139,6 @@ version = "0.0.1"
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 version = "1.11.0"
 
-[[deps.LineSearches]]
-deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"]
-git-tree-sha1 = "4adee99b7262ad2a1a4bbbc59d993d24e55ea96f"
-uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
-version = "7.4.0"
-
 [[deps.LinearAlgebra]]
 deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -382,69 +175,16 @@ version = "1.2.0"
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 version = "1.11.0"
 
-[[deps.NLSolversBase]]
-deps = ["ADTypes", "DifferentiationInterface", "Distributed", "FiniteDiff", "ForwardDiff"]
-git-tree-sha1 = "25a6638571a902ecfb1ae2a18fc1575f86b1d4df"
-uuid = "d41bc354-129a-5804-8e4c-c37616107c6c"
-version = "7.10.0"
-
-[[deps.NaNMath]]
-deps = ["OpenLibm_jll"]
-git-tree-sha1 = "9b8215b1ee9e78a293f99797cd31375471b2bcae"
-uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
-version = "1.1.3"
-
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
 version = "0.3.27+1"
 
-[[deps.OpenLibm_jll]]
-deps = ["Artifacts", "Libdl"]
-uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
-version = "0.8.5+0"
-
-[[deps.OpenSpecFun_jll]]
-deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "1346c9208249809840c91b26703912dff463d335"
-uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
-version = "0.5.6+0"
-
-[[deps.Optim]]
-deps = ["Compat", "EnumX", "FillArrays", "ForwardDiff", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"]
-git-tree-sha1 = "61942645c38dd2b5b78e2082c9b51ab315315d10"
-uuid = "429524aa-4258-5aef-a3af-852621145aeb"
-version = "1.13.2"
-
-    [deps.Optim.extensions]
-    OptimMOIExt = "MathOptInterface"
-
-    [deps.Optim.weakdeps]
-    MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
-
 [[deps.OrderedCollections]]
 git-tree-sha1 = "05868e21324cede2207c6f0f466b4bfef6d5e7ee"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 version = "1.8.1"
 
-[[deps.Parameters]]
-deps = ["OrderedCollections", "UnPack"]
-git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
-uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
-version = "0.12.3"
-
-[[deps.PositiveFactorizations]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "17275485f373e6673f7e7f97051f703ed5b15b20"
-uuid = "85a6dd25-e78a-55b7-8502-1745935b8125"
-version = "0.2.4"
-
-[[deps.Preferences]]
-deps = ["TOML"]
-git-tree-sha1 = "0f27480397253da18fe2c12a4ba4eb9eb208bf3d"
-uuid = "21216c6a-2e73-6563-6e65-726566657250"
-version = "1.5.0"
-
 [[deps.Primes]]
 deps = ["IntegerMathUtils"]
 git-tree-sha1 = "25cdd1d20cd005b52fc12cb6be3f75faaf59bb9b"
@@ -492,22 +232,12 @@ version = "0.7.0"
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 version = "1.11.0"
 
-[[deps.Setfield]]
-deps = ["ConstructionBase", "Future", "MacroTools", "StaticArraysCore"]
-git-tree-sha1 = "c5391c6ace3bc430ca630251d02ea9687169ca68"
-uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
-version = "1.1.2"
-
 [[deps.Sobol]]
 deps = ["DelimitedFiles", "Random"]
 git-tree-sha1 = "5a74ac22a9daef23705f010f72c81d6925b19df8"
 uuid = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4"
 version = "1.5.0"
 
-[[deps.Sockets]]
-uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
-version = "1.11.0"
-
 [[deps.SortingAlgorithms]]
 deps = ["DataStructures"]
 git-tree-sha1 = "64d974c2e6fdf07f8155b5b2ca2ffa9069b608d9"
@@ -519,23 +249,6 @@ deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 version = "1.11.0"
 
-[[deps.SpecialFunctions]]
-deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "41852b8679f78c8d8961eeadc8f62cef861a52e3"
-uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "2.5.1"
-
-    [deps.SpecialFunctions.extensions]
-    SpecialFunctionsChainRulesCoreExt = "ChainRulesCore"
-
-    [deps.SpecialFunctions.weakdeps]
-    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-
-[[deps.StaticArraysCore]]
-git-tree-sha1 = "192954ef1208c7019899fbf8049e717f92959682"
-uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
-version = "1.4.3"
-
 [[deps.Statistics]]
 deps = ["LinearAlgebra"]
 git-tree-sha1 = "ae3bb1eb3bba077cd276bc5cfc337cc65c3075c0"
@@ -563,21 +276,11 @@ deps = ["Artifacts", "Libdl", "libblastrampoline_jll"]
 uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
 version = "7.7.0+0"
 
-[[deps.TOML]]
-deps = ["Dates"]
-uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
-version = "1.0.3"
-
 [[deps.UUIDs]]
 deps = ["Random", "SHA"]
 uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 version = "1.11.0"
 
-[[deps.UnPack]]
-git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
-uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
-version = "1.0.2"
-
 [[deps.Unicode]]
 uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 version = "1.11.0"
diff --git a/Project.toml b/Project.toml
index 46732f3..8b0f5da 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,17 +4,15 @@ authors = ["Shuhei Ohno", "Martin Mikkelsen"]
 version = "1.0.5"
 
 [deps]
+FewBodyHamiltonians = "3a126c26-e5d7-4a95-83c3-3b69f8a11ded"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
 Aqua = "0.8.13"
+FewBodyHamiltonians = "0.0.1"
 LinearAlgebra = "1.7.3"
-Optim = "1.7.8"
 QuasiMonteCarlo = "0.3.3"
-SpecialFunctions = "2.5.0"
 Test = "1.11.0"
 julia = "1.7"
 
diff --git a/src/FewBodyECG.jl b/src/FewBodyECG.jl
index 1de026e..36f1cb2 100644
--- a/src/FewBodyECG.jl
+++ b/src/FewBodyECG.jl
@@ -5,24 +5,6 @@ include("coordinates.jl")
 include("matrix_elements.jl")
 include("hamiltonian.jl")
 include("sampling.jl")
-include("optimization.jl")
 include("utils.jl")
 
-using .Types
-using .Coordinates
-using .MatrixElements
-using .Hamiltonian
-using .Sampling
-using .Optimization
-using .Utils
-
-export generate_A_matrix, generate_bij, default_b0, jacobi_transform, transform_list, transform_coordinates, inverse_transform_coordinates, ParticleSystem, Particle, GaussianBase, Rank0Gaussian, Rank1Gaussian, Rank2Gaussian, BasisSet, Operator, KineticEnergy, CoulombPotential, FewBodyHamiltonian, MatrixElementResult
-
-export compute_matrix_element, build_overlap_matrix, build_operator_matrix,
-    build_hamiltonian_matrix, solve_generalized_eigenproblem,
-    generate_basis, compute_ground_state_energy, optimize_ground_state_energy
-
-export ψ₀, plot_wavefunction, plot_density
-
-
 end
diff --git a/src/coordinates.jl b/src/coordinates.jl
index 449ebd6..30acf66 100644
--- a/src/coordinates.jl
+++ b/src/coordinates.jl
@@ -182,4 +182,4 @@ function inverse_transform_coordinates(U::Matrix{Float64}, x::Vector{Float64})::
     return U * x
 end
 
-end 
+end
diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl
index eb13eb3..e4dc292 100644
--- a/src/hamiltonian.jl
+++ b/src/hamiltonian.jl
@@ -3,18 +3,17 @@ module Hamiltonian
 using LinearAlgebra
 using ..Types
 using ..MatrixElements
+using FewBodyHamiltonians
 
 export build_overlap_matrix, build_operator_matrix, build_hamiltonian_matrix, solve_generalized_eigenproblem
 
-struct IdentityOperator <: Operator end
-
 function compute_overlap_element(bra::GaussianBase, ket::GaussianBase)
     A, B = bra.A, ket.A
     R = inv(A + B)
     return (π^length(R) / det(A + B))^(3 / 2)
 end
 
-function build_overlap_matrix(basis::BasisSet)
+function build_overlap_matrix(basis::ECGBasis)
     n = length(basis.functions)
     S = zeros(n, n)
     for i in 1:n, j in 1:i
@@ -24,7 +23,7 @@ function build_overlap_matrix(basis::BasisSet)
     return S
 end
 
-function build_operator_matrix(basis::BasisSet, op::Operator)
+function build_operator_matrix(basis::ECGBasis, op::Operator)
     n = length(basis.functions)
     H = zeros(n, n)
     for i in 1:n, j in 1:i
@@ -34,7 +33,7 @@ function build_operator_matrix(basis::BasisSet, op::Operator)
     return H
 end
 
-function build_hamiltonian_matrix(basis::BasisSet, operators::Vector{Operator})
+function build_hamiltonian_matrix(basis::ECGBasis, operators::Vector{Operator})
     H = zeros(length(basis.functions), length(basis.functions))
     for op in operators
         H .+= build_operator_matrix(basis, op)
@@ -47,4 +46,4 @@ function solve_generalized_eigenproblem(H::Matrix{Float64}, S::Matrix{Float64})
     return real(vals), real(vecs)
 end
 
-end # module
+end
diff --git a/src/matrix_elements.jl b/src/matrix_elements.jl
index d745870..b6bcf43 100644
--- a/src/matrix_elements.jl
+++ b/src/matrix_elements.jl
@@ -1,85 +1,72 @@
 module MatrixElements
 
 using LinearAlgebra
+using FewBodyHamiltonians
 using ..Types
-using SpecialFunctions: erf
 
 export compute_matrix_element
 
-"""
-compute_matrix_element(bra, ket, op)
 
-Compute the matrix element ⟨bra|op|ket⟩ using analytic expressions.
-"""
+flatten1(x) = x isa AbstractVector{<:AbstractVector} ? vec(x[1]) : vec(x)
 
-function compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian, op::KineticEnergy)
+function compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian, op::Kinetic, K::AbstractMatrix)
     A, B = bra.A, ket.A
-    K = op.K
     R = inv(A + B)
-    M0 = (π^length(R) / det(A + B))^(3 / 2)
+    M0 = (π^size(R, 1) / det(A + B))^(3 / 2)
     return 6 * tr(B * K * A * R) * M0
 end
 
-function compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian, op::CoulombPotential)
-    A, B, w = bra.A, ket.A, op.w
+function compute_matrix_element(bra::Rank0Gaussian, ket::Rank0Gaussian, op::Coulomb, w::AbstractVector)
+    A, B = bra.A, ket.A
     R = inv(A + B)
     β = 1 / (dot(w, R * w))
-    M0 = (π^length(R) / det(A + B))^(3 / 2)
+    M0 = (π^size(R, 1) / det(A + B))^(3 / 2)
     return op.coefficient * 2 * sqrt(β / π) * M0
 end
 
-function compute_matrix_element(bra::Rank1Gaussian, ket::Rank1Gaussian, op::CoulombPotential)
-    A, B, a, b, w = bra.A, ket.A, bra.a, ket.a, op.w
+function compute_matrix_element(bra::Rank1Gaussian, ket::Rank1Gaussian, op::Coulomb, w::AbstractVector)
+    A, B = bra.A, ket.A
+    a, b = flatten1(bra.a), flatten1(ket.a)
     R = inv(A + B)
     β = 1 / (dot(w, R * w))
-    M0 = (π^length(R) / det(A + B))^(3 / 2)
+    M0 = (π^size(R, 1) / det(A + B))^(3 / 2)
     M1 = 0.5 * dot(b, R * a) * M0
     q2 = 0.25 * dot(a .+ b, R * (w * w') * (a .+ b))
     return 2 * sqrt(β / π) * M1 - sqrt(β^3 / π) / 3 * q2 * M0
 end
 
-function compute_matrix_element(bra::Rank1Gaussian, ket::Rank1Gaussian, op::KineticEnergy)
-    a = bra.a isa AbstractVector{<:AbstractVector} ? vec(bra.a[1]) : vec(bra.a)
-    b = ket.a isa AbstractVector{<:AbstractVector} ? vec(ket.a[1]) : vec(ket.a)
-
+function compute_matrix_element(bra::Rank1Gaussian, ket::Rank1Gaussian, op::Kinetic, K::AbstractMatrix)
     A, B = bra.A, ket.A
-    K = op.K
+    a, b = flatten1(bra.a), flatten1(ket.a)
     R = inv(A + B)
-    M0 = (π^length(R) / det(A + B))^(3 / 2)
+    M0 = (π^size(R, 1) / det(A + B))^(3 / 2)
     M1 = 0.5 * dot(b, R * a) * M0
-
     T1 = 6 * tr(B * K * A * R) * M1
     T2 = dot(b, a) * M0
     T3 = dot(a, R * B * A * R * b) * M0
     T4 = dot(b, R * B * a) * M0
     T5 = dot(a, R * A * b) * M0
-
     return T1 + T2 + T3 - T4 - T5
 end
 
-
-function compute_matrix_element(bra::Rank2Gaussian, ket::Rank2Gaussian, op::KineticEnergy)
+function compute_matrix_element(bra::Rank2Gaussian, ket::Rank2Gaussian, op::Kinetic, K::AbstractMatrix)
     A, B = bra.A, ket.A
-    a, b, c, d = bra.a, bra.b, ket.a, ket.b
-    K = op.K
+    a, b = flatten1(bra.a), flatten1(bra.b)
+    c, d = flatten1(ket.a), flatten1(ket.b)
     R = inv(A + B)
-    M0 = (π^length(R) / det(A + B))^(3 / 2)
-
+    M0 = (π^size(R, 1) / det(A + B))^(3 / 2)
     M2 = 0.25 * (
         dot(a, R * b) * dot(c, R * d) +
             dot(a, R * c) * dot(b, R * d) +
             dot(a, R * d) * dot(b, R * c)
     ) * M0
-
     T1 = 6 * tr(B * K * A * R) * M2
-
     T2 = 0.5 * (
         dot(a, K * c) * dot(b, R * d) +
             dot(a, K * d) * dot(b, R * c) +
             dot(b, K * c) * dot(a, R * d) +
             dot(b, K * d) * dot(a, R * c)
     ) * M0
-
     T3 = 0.5 * (
         dot(a, R * B * K * A * R * b) * dot(c, R * d) +
             dot(a, R * B * K * A * R * c) * dot(b, R * d) +
@@ -94,7 +81,6 @@ function compute_matrix_element(bra::Rank2Gaussian, ket::Rank2Gaussian, op::Kine
             dot(d, R * B * K * A * R * b) * dot(a, R * c) +
             dot(d, R * B * K * A * R * c) * dot(a, R * b)
     ) * M0
-
     T4 = -0.5 * (
         dot(a, R * B * K * b) * dot(c, R * d) +
             dot(b, R * B * K * a) * dot(c, R * d) +
@@ -103,7 +89,6 @@ function compute_matrix_element(bra::Rank2Gaussian, ket::Rank2Gaussian, op::Kine
             dot(d, R * B * K * a) * dot(b, R * c) +
             dot(d, R * B * K * b) * dot(a, R * c)
     ) * M0
-
     T5 = -0.5 * (
         dot(c, K * A * a) * dot(b, R * d) +
             dot(c, K * A * b) * dot(a, R * d) +
@@ -112,49 +97,35 @@ function compute_matrix_element(bra::Rank2Gaussian, ket::Rank2Gaussian, op::Kine
             dot(d, K * A * b) * dot(a, R * c) +
             dot(d, K * A * c) * dot(a, R * b)
     ) * M0
-
     return T1 + T2 + T3 + T4 + T5
 end
 
-function compute_matrix_element(bra::Rank2Gaussian, ket::Rank2Gaussian, op::CoulombPotential)
+function compute_matrix_element(bra::Rank2Gaussian, ket::Rank2Gaussian, op::Coulomb, w::AbstractVector)
     A, B = bra.A, ket.A
-    a, b, c, d = bra.a, bra.b, ket.a, ket.b
-    w = op.w
+    a, b = flatten1(bra.a), flatten1(bra.b)
+    c, d = flatten1(ket.a), flatten1(ket.b)
     R = inv(A + B)
-    M0 = (π^length(R) / det(A + B))^(3 / 2)
-
+    M0 = (π^size(R, 1) / det(A + B))^(3 / 2)
     β = 1 / dot(w, R * w)
     q = 0.5 * dot(w, R * (a + b + c + d))
-
     M2 = 0.25 * (
         dot(a, R * b) * dot(c, R * d) +
             dot(a, R * c) * dot(b, R * d) +
             dot(a, R * d) * dot(b, R * c)
     ) * M0
-
     term1 = 2 * sqrt(β / π) * M2
-
     q2_1 = dot(a, R * (w * w') * b) * dot(c, R * d)
     q2_2 = dot(a, R * (w * w') * c) * dot(b, R * d)
     q2_3 = dot(a, R * (w * w') * d) * dot(b, R * c)
     q2_4 = dot(b, R * (w * w') * c) * dot(a, R * d)
     q2_5 = dot(b, R * (w * w') * d) * dot(a, R * c)
     q2_6 = dot(c, R * (w * w') * d) * dot(a, R * b)
-
     q4_1 = dot(a, R * (w * w') * b) * dot(c, R * (w * w') * d)
     q4_2 = dot(a, R * (w * w') * c) * dot(b, R * (w * w') * d)
     q4_3 = dot(a, R * (w * w') * d) * dot(b, R * (w * w') * c)
-
-    term2 = -2 * sqrt(β / π) * β / 3 * 0.25 * (
-        q2_1 + q2_2 + q2_3 + q2_4 + q2_5 + q2_6
-    ) * M0
-
-    term3 = 2 * sqrt(β / π) * β^2 / 10 * 0.5 * (
-        q4_1 + q4_2 + q4_3
-    ) * M0
-
+    term2 = -2 * sqrt(β / π) * β / 3 * 0.25 * (q2_1 + q2_2 + q2_3 + q2_4 + q2_5 + q2_6) * M0
+    term3 = 2 * sqrt(β / π) * β^2 / 10 * 0.5 * (q4_1 + q4_2 + q4_3) * M0
     return op.coefficient * (term1 + term2 + term3)
 end
 
-
-end 
+end
diff --git a/src/optimization.jl b/src/optimization.jl
deleted file mode 100644
index 013a4fe..0000000
--- a/src/optimization.jl
+++ /dev/null
@@ -1,30 +0,0 @@
-module Optimization
-
-using Optim
-using ..Types
-using ..Coordinates
-using ..Hamiltonian
-
-export optimize_ground_state_energy
-
-"""
-optimize_ground_state_energy(init_widths::Vector{Matrix{Float64}}, ops::Vector{Operator}; max_iter=100)
-
-Uses local optimization (Nelder-Mead) on Gaussian widths to minimize ground state energy.
-"""
-function optimize_ground_state_energy(init_widths::Vector{Matrix{Float64}}, ops::Vector{Operator}; max_iter::Int = 100)
-    vecdim = size(init_widths[1], 1)
-    nfuncs = length(init_widths)
-    flat_init = vcat([vec(A) for A in init_widths]...)
-
-    function objective(x)
-        widths = [reshape(x[((i - 1) * vecdim^2 + 1):(i * vecdim^2)], vecdim, vecdim) for i in 1:nfuncs]
-        basis = generate_basis(widths)
-        return compute_ground_state_energy(basis, ops)
-    end
-
-    result = optimize(objective, flat_init, NelderMead(); iterations = max_iter)
-    return Optim.minimum(result)
-end
-
-end 
diff --git a/src/sampling.jl b/src/sampling.jl
index 76d604f..e825a7f 100644
--- a/src/sampling.jl
+++ b/src/sampling.jl
@@ -6,12 +6,13 @@ using ..Hamiltonian
 using LinearAlgebra
 using QuasiMonteCarlo
 
-export generate_basis, compute_ground_state_energy, generate_bij
+export generate_basis, generate_bij
+
 
 """
 generate_basis(widths::Vector{Matrix{Float64}}, rank::Int=0)
 
-Construct a `BasisSet` from a list of correlation matrices and optional rank.
+Construct a `ECGBasis` from a list of correlation matrices and optional rank.
 """
 function generate_basis(widths::Vector{Matrix{Float64}}, rank::Int = 0)
     if rank == 0
@@ -19,7 +20,7 @@ function generate_basis(widths::Vector{Matrix{Float64}}, rank::Int = 0)
     else
         error("Only Rank0Gaussian implemented in generate_basis")
     end
-    return BasisSet(funcs)
+    return ECGBasis(funcs)
 end
 
 """
@@ -37,16 +38,4 @@ function generate_bij(method::Symbol, i::Int, n_terms::Int, b1::Float64; qmc_sam
     end
 end
 
-"""
-compute_ground_state_energy(basis::BasisSet, ops::Vector{Operator})
-
-Construct the Hamiltonian and overlap matrices and return the lowest eigenvalue.
-"""
-function compute_ground_state_energy(basis::BasisSet, ops::Vector{Operator})
-    H = build_hamiltonian_matrix(basis, ops)
-    S = build_overlap_matrix(basis)
-    vals, _ = solve_generalized_eigenproblem(H, S)
-    return minimum(vals)
 end
-
-end 
diff --git a/src/types.jl b/src/types.jl
index c080692..fc131ba 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -1,60 +1,46 @@
 module Types
 
-export Particle,
-    GaussianBase, Rank0Gaussian, Rank1Gaussian, Rank2Gaussian,
-    BasisSet,
-    Operator, KineticEnergy, CoulombPotential,
-    FewBodyHamiltonian, MatrixElementResult
+using FewBodyHamiltonians
+
+export Operator, Hamiltonian, Kinetic, Coulomb, Particle, System, GaussianBase, Rank0Gaussian, Rank1Gaussian, Rank2Gaussian, ECGBasis, MatrixElementResult
 
 abstract type GaussianBase end
 
 struct Particle
-    mass::Float64
-    charge::Float64
-    label::Symbol
-end
-
-struct Rank0Gaussian <: GaussianBase
-    A::Matrix{Float64}  # Correlation matrix
+    mass::Real
+    charge::Real
+    spin::Union{Nothing, Real}
 end
 
-struct Rank1Gaussian <: GaussianBase
-    A::Matrix{Float64}                 # Correlation matrix
-    a::Vector{Vector{Float64}}         # Polarization vectors
+struct System
+    particles::Vector{Particle}
+    remove_com::Bool
 end
 
-struct Rank2Gaussian <: GaussianBase
-    A::Matrix{Float64}
-    a::Vector{Vector{Float64}}         # First polarization vector set
-    b::Vector{Vector{Float64}}         # Second polarization vector set
+struct Rank0Gaussian{T <: Real, M <: AbstractMatrix{T}} <: GaussianBase
+    A::M
 end
 
-struct BasisSet
-    functions::Vector{GaussianBase}
+struct Rank1Gaussian{T <: Real, M <: AbstractMatrix{T}, V <: AbstractVector{<:AbstractVector{T}}} <: GaussianBase
+    A::M
+    a::V
 end
 
-abstract type Operator end
-
-struct KineticEnergy <: Operator
-    K::Matrix{Float64}  # required for ⟨T⟩
+struct Rank2Gaussian{T <: Real, M <: AbstractMatrix{T}, V <: AbstractVector{<:AbstractVector{T}}} <: GaussianBase
+    A::M
+    a::V
+    b::V
 end
 
-struct CoulombPotential <: Operator
-    coefficient::Float64
-    w::Vector{Float64}
+struct ECGBasis{F <: GaussianBase}
+    functions::Vector{F}
 end
 
-
-struct FewBodyHamiltonian
-    basis::BasisSet
-    operators::Vector{Operator}
+struct MatrixElementResult{B <: GaussianBase, O <: Operator, T <: Real}
+    bra::B
+    ket::B
+    operator::O
+    value::T
 end
 
-struct MatrixElementResult
-    bra::GaussianBase
-    ket::GaussianBase
-    operator::Operator
-    value::Float64
 end
-
-end # module