Runic PR

.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'
+          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

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
-

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)

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)
+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)

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