Assemble Hamiltonian by operator

docs/src/Rayleigh-Ritz.md

@@ -292,6 +292,10 @@ println("  Reference: ", -0.495010)
 
 The results agree with those in Table 1 of the Supporting Information. The small differences between the calculated and reference values arise from rounding in the published parameters.
 
+## Acknowledgments
+
+We thank the participants of the [JuliaHEP Workshop 2025](https://indico.cern.ch/e/juliahep2025) for feedback on matrix assembly.
+
 ## API reference
 
 ### Solver

src/Rayleigh-Ritz.jl

@@ -365,11 +365,12 @@ end
 
 function matrix(hamiltonian::Hamiltonian, basisset::BasisSet)
   nₘₐₓ = length(basisset.basis)
-  # H = [element(hamiltonian, basisset.basis[i], basisset.basis[j]) for i=1:nₘₐₓ, j=1:nₘₐₓ]
-  H = Array{Float64}(undef, nₘₐₓ, nₘₐₓ)
-  for j in 1:nₘₐₓ
-    for i in 1:j
-      H[i,j] = element(hamiltonian, basisset.basis[i], basisset.basis[j])
+  H = zeros(Float64, nₘₐₓ, nₘₐₓ)
+  for operator in hamiltonian.terms
+    for j in 1:nₘₐₓ
+      for i in 1:j
+        H[i,j] += element(operator, basisset.basis[i], basisset.basis[j])
+      end
     end
   end
   return LinearAlgebra.Symmetric(H)

test/Rayleigh-Ritz.jl

@@ -13,6 +13,8 @@
     SimpleGaussianBasis(0.1219492),
   )
   res = solve(H, BS)
+
+  @test TwoBody.matrix(H, BS) ≈ sum(TwoBody.matrix(term, BS) for term in H.terms)
 
   println("4π×∫|ψ(r)|²r²dr = 1")
   println("  i\tnumerical  \tanalytical")