Make compatible with #5

.gitignore

@@ -3,4 +3,5 @@ src/test.py
 src/test.jl
 .DS_Store
 .vscode
-docs/build/
+docs/build/
+Manifest.toml

Examples/HydrogenAnion.jl

@@ -12,44 +12,19 @@ K_transformed = psys.J * K * psys.J'
 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]
 
-let
-    n_basis = 50
-    b1 = default_b0(psys.scale)
-    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())
-        A = generate_A_matrix(bij, w_raw)
-        push!(basis_fns, Rank0Gaussian(A))
-
-        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)
-        E₀_step = minimum(vals)
-
-        push!(E₀_list, E₀_step)
-        println("Step $i: E₀ = $E₀_step")
-
-    end
-
-    E₀ = minimum(E₀_list)
-    Eᵗʰ = -0.527751016523
-    Δ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"
-    )
-
-end
+
+ops = Operator[
+    KineticOperator(K_transformed);
+    (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))...
+]
+
+A = solve_ECG(ops, psys, 100)
+
+E = -0.527751016523
+ΔE = abs(A.ground_state - E)
+@info "Energy difference" ΔE
+
+n, E = convergence(A)
+plot(n, E)

Examples/Hydrogen_p-wave.jl

@@ -21,15 +21,23 @@ E₀_list = Float64[]
 
 a_vec = [1.0]
 
+ops = Operator[
+    KineticOperator(K_transformed);
+    (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))...
+]
+
+A = solve_ECG(ops, psys, 50)
+
 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]))
+    A = _generate_A_matrix(bij, w_raw)
+    # pass a plain vector for the rank-1 displacement
+    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))...
+        KineticOperator(K_transformed);
+        (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))...
     ]
 
     H = build_hamiltonian_matrix(basis, ops)
@@ -52,7 +60,7 @@ for i in 1:n_basis
     println("Step $i: E₀ = $E₀")
 end
 
-E_exact = -0.125  #
+E_exact = -0.125
 E_min = minimum(E₀_list)
 @show ΔE = abs(E_min - E_exact)
 

Examples/Hydrogen_s-wave.jl

@@ -3,6 +3,8 @@ using LinearAlgebra
 using Plots
 using QuasiMonteCarlo
 
+import FewBodyECG: _generate_A_matrix
+
 masses = [1.0e15, 1.0]  # proton, electron
 psys = ParticleSystem(masses)
 
@@ -19,16 +21,18 @@ b1 = 1.5
 basis_fns = GaussianBase[]
 E₀_list = Float64[]
 
+ops = Operator[
+    KineticOperator(K_transformed);
+    (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))...
+]
+
+
 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)
+    A = _generate_A_matrix(bij, w_raw)
     push!(basis_fns, Rank0Gaussian(A))
 
     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)

Examples/Positronium.jl

@@ -14,6 +14,13 @@ 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)
@@ -23,13 +30,13 @@ let
 
     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)
+        A = _generate_A_matrix(bij, w_raw)
         push!(basis_fns, Rank0Gaussian(A))
 
         basis = BasisSet(basis_fns)
         ops = Operator[
-            KineticEnergy(K_transformed);
-            (CoulombPotential(c, w) for (c, w) in zip(coeffs, w_raw))...
+            KineticOperator(K_transformed);
+            (CoulombOperator(c, w) for (c, w) in zip(coeffs, w_raw))...
         ]
 
         H = build_hamiltonian_matrix(basis, ops)
@@ -39,13 +46,13 @@ let
         E₀ = minimum(vals)
 
         push!(E₀_list, E₀)
-        println("Step $i: E₀ = $E₀")
+        @info "Step $i" E₀ = E₀
     end
 
     E₀ = minimum(E₀_list)
     Eᵗʰ = -0.2620050702328
     ΔE = abs(E₀ - Eᵗʰ)
-    @show Δ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]