From 3986d19ea7b19469f20eac9ca19581fc279f58cc Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Fri, 14 Aug 2026 11:24:35 +0900
Subject: [PATCH] Assemble Hamiltonian by operator

---
 docs/src/Rayleigh-Ritz.md |  4 ++++
 src/Rayleigh-Ritz.jl      | 11 ++++++-----
 test/Rayleigh-Ritz.jl     |  2 ++
 3 files changed, 12 insertions(+), 5 deletions(-)

diff --git a/docs/src/Rayleigh-Ritz.md b/docs/src/Rayleigh-Ritz.md
index dd70dbb..56caccc 100644
--- a/docs/src/Rayleigh-Ritz.md
+++ b/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
diff --git a/src/Rayleigh-Ritz.jl b/src/Rayleigh-Ritz.jl
index 4ac8c17..73f492b 100644
--- a/src/Rayleigh-Ritz.jl
+++ b/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)
diff --git a/test/Rayleigh-Ritz.jl b/test/Rayleigh-Ritz.jl
index 6839826..e1ae146 100644
--- a/test/Rayleigh-Ritz.jl
+++ b/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")