From fbab9c0316012b1852018f7566747bbeae25b67d Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Tue, 2 Apr 2024 18:48:41 -0400
Subject: [PATCH] Critical bug fixed: missing Kronecker deltas

---
 ho_basis.jl       |  6 ++++--
 ho_basis_3body.jl | 10 ++++++----
 2 files changed, 10 insertions(+), 6 deletions(-)

diff --git a/ho_basis.jl b/ho_basis.jl
index 8547e5a..c1f8292 100644
--- a/ho_basis.jl
+++ b/ho_basis.jl
@@ -54,9 +54,10 @@ function get_2p_basis(E_max)
     return (Es, n1s, l1s, n2s, l2s)
 end
 
-function sp_T_matrix(ns, ls; ω=1.0, μ=1.0)
+function sp_T_matrix(ns, ls; mask=trues(length(ns),length(ns)), ω=1.0, μ=1.0)
     mat = spzeros(length(ns), length(ns))
     for idx in CartesianIndices(mat)
+        if !mask[idx]; continue; end
         (i, j) = Tuple(idx)
         if ls[i] == ls[j]
             if ns[i] == ns[j]
@@ -70,9 +71,10 @@ function sp_T_matrix(ns, ls; ω=1.0, μ=1.0)
     return (ω / μ) .* mat
 end
 
-function sp_V_matrix(V_l, ns, ls; dtype=Float64)
+function sp_V_matrix(V_l, ns, ls; mask=trues(length(ns),length(ns)), dtype=Float64)
     mat = zeros(dtype, length(ns), length(ns))
     Threads.@threads for idx in CartesianIndices(mat)
+        if !mask[idx]; continue; end
         (i, j) = Tuple(idx)
         if ls[i] == ls[j]
             mat[idx] = V_l(ls[i], ns[i], ns[j])
diff --git a/ho_basis_3body.jl b/ho_basis_3body.jl
index 5cac046..7ec5c13 100644
--- a/ho_basis_3body.jl
+++ b/ho_basis_3body.jl
@@ -18,17 +18,19 @@ c2 = 2
 println("No of threads = ", Threads.nthreads())
 
 Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max)
+mask1 = (n2s .== n2s') .&& (l2s .== l2s')
+mask2 = (n1s .== n1s') .&& (l1s .== l1s')
 println("Basis size = ", length(Es))
 
 println("Constructing KE matrices")
-@time "T1" T1 = sp_T_matrix(n1s, l1s; ω=ω, μ=μ1)
-@time "T2" T2 = sp_T_matrix(n2s, l2s; ω=ω, μ=c2^2 * μ2)
+@time "T1" T1 = sp_T_matrix(n1s, l1s; mask=mask1, ω=ω, μ=μ1)
+@time "T2" T2 = sp_T_matrix(n2s, l2s; mask=mask2, ω=ω, μ=c2^2 * μ2)
 
 println("Constructing PE matrices")
 V1_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; ω=ω)
 V_relative_elem(l, n1, n2) = Va * V_Gaussian(Ra / c, l, n1, n2; ω=ω)
-@time "V1" V1 = sp_V_matrix(V1_elem, n1s, l1s)
-@time "V relative" V_relative = sp_V_matrix(V_relative_elem, n1s, l1s) + sp_V_matrix(V_relative_elem, n2s, l2s)
+@time "V1" V1 = sp_V_matrix(V1_elem, n1s, l1s; mask=mask1)
+@time "V relative" V_relative = sp_V_matrix(V_relative_elem, n1s, l1s; mask=mask1) + sp_V_matrix(V_relative_elem, n2s, l2s; mask=mask2)
 @time "Moshinsky brackets" U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
 @time "V2" V2 =  U' * V_relative * U