Updated Gaussians

Examples/HydrogenAnion.jl

@@ -4,27 +4,26 @@ using Plots
 using QuasiMonteCarlo
 
 masses = [1.0e15, 1.0, 1.0]
-psys = ParticleSystem(masses)
 
-K = Diagonal([0.0, 1 / 2, 1 / 2])
-K_transformed = psys.J * K * psys.J'
+Λmat = Λ(masses)
+kin = KineticOperator(Λmat)
+J, U = _jacobi_transform(masses)
 
 w_list = [[1, -1, 0], [1, 0, -1], [0, 1, -1]]
 
-w_raw = [psys.U' * w for w in w_list]
+w_raw = [U' * w for w in w_list]
 coeffs = [-1.0, -1.0, +1.0]
 
-
 ops = Operator[
-    KineticOperator(K_transformed);
+    kin;
     (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))...
 ]
 
-A = solve_ECG(ops, psys, 100)
+result = solve_ECG(ops, 250, scale = 1.0)
 
 E = -0.527751016523
-ΔE = abs(A.ground_state - E)
+ΔE = abs(result.ground_state - E)
 @info "Energy difference" ΔE
 
-n, E = convergence(A)
+n, E = convergence(result)
 plot(n, E)

Examples/Positronium.jl

@@ -1,73 +1,22 @@
-using FewBodyECG
-using LinearAlgebra
-using Plots
+using FewBodyECG, LinearAlgebra
 using QuasiMonteCarlo
+using Plots
+import FewBodyECG: default_scale, convergence
 
 masses = [1.0, 1.0, 1.0]
-psys = ParticleSystem(masses)
 
-K = Diagonal([1 / 2, 1 / 2, 1 / 2])
-K_transformed = psys.J * K * psys.J'
+Λmat = Λ(masses)
+kin = KineticOperator(Λmat)
 
+J, U = _jacobi_transform(masses)
 w_list = [[1, -1, 0], [1, 0, -1], [0, 1, -1]]
-w_raw = [psys.U' * w for w in w_list]
-
-coeffs = [+1.0, -1.0, -1.0]
-
-ops = Operator[
-    KineticOperator(K_transformed);
-    (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))...
-]
-
-A = solve_ECG(ops, psys, 100)
-
-let
-    n_basis = 100
-    b1 = default_b0(psys.scale)
-    method = :quasirandom
-    basis_fns = GaussianBase[]
-    E₀_list = Float64[]
-
-    for i in 1:n_basis
-        bij = generate_bij(method, i, length(w_raw), b1; qmc_sampler = SobolSample())
-        A = _generate_A_matrix(bij, w_raw)
-        push!(basis_fns, Rank0Gaussian(A))
-
-        basis = BasisSet(basis_fns)
-        ops = Operator[
-            KineticOperator(K_transformed);
-            (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))...
-        ]
-
-        H = build_hamiltonian_matrix(basis, ops)
-        S = build_overlap_matrix(basis)
-        vals, vecs = solve_generalized_eigenproblem(H, S)
-        global c₀ = vecs[:, 1]
-        E₀ = minimum(vals)
-
-        push!(E₀_list, E₀)
-        @info "Step $i" E₀ = E₀
-    end
-
-    E₀ = minimum(E₀_list)
-    Eᵗʰ = -0.2620050702328
-    ΔE = abs(E₀ - Eᵗʰ)
-    @info "Energy difference" ΔE
-
-    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"
-    )
-
+w_raw = [U' * w for w in w_list]
 
-    p2 = plot(
-        1:n_basis, E₀_list, xlabel = "Number of Gaussians", ylabel = "E₀ [Hartree]",
-        lw = 2, label = "Ground state energy", title = "Positronium Convergence"
-    )
+coeffs = [-1.0, -1.0, +1.0]
+coulomb_ops = [CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw)]
 
-    plot(p1, p2, layout = (2, 1))
+ops = Operator[kin; coulomb_ops...]
+scale = default_scale(masses)
 
-end
+result = solve_ECG(ops, 300, sampler = SobolSample(); scale = scale)
+println("E ≈ ", result.ground_state)

Examples/tdμ.jl

@@ -0,0 +1,25 @@
+using FewBodyECG, LinearAlgebra
+using QuasiMonteCarlo
+import FewBodyECG: default_scale, convergence
+using Plots
+
+masses = [5496.918, 3670.481, 206.7686]
+
+Λmat = Λ(masses)
+kin = KineticOperator(Λmat)
+
+J, U = _jacobi_transform(masses)
+base_w = [[1, -1, 0], [1, 0, -1], [0, 1, -1]]
+coeffs = [+1.0, -1.0, -1.0]
+w_list = [c .* w for (c, w) in zip(coeffs, base_w)]
+w_raw = [U' * w for w in w_list]
+
+coulomb_ops = [CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw)]
+
+ops = Operator[kin; coulomb_ops...]
+scale = default_scale(masses)
+result = solve_ECG(ops, 550, sampler = HaltonSample(); scale = scale)
+println("E ≈ ", result.ground_state)
+
+a, b = correlation_function(result)
+conv = FewBodyECG.convergence(result)

Project.toml

@@ -1,7 +1,7 @@
 name = "FewBodyECG"
 uuid = "083b1810-24a1-4a79-9a41-145bb2bb8ceb"
-authors = ["Shuhei Ohno", "Martin Mikkelsen"]
 version = "1.0.5"
+authors = ["Shuhei Ohno", "Martin Mikkelsen"]
 
 [deps]
 FewBodyHamiltonians = "3a126c26-e5d7-4a95-83c3-3b69f8a11ded"