Moshinsky in k-space implemented (not working)

berggren.jl

@@ -16,3 +16,15 @@ function get_W_matrix(basis_p, E_max, Λ, μ1ω1, μ2ω2=μ1ω1; weights=true)
     end
     return sparse(W)
 end
+
+"Moshinsky transform in k-space"
+function get_Moshinsky_in_kspace(E_max, basis_p, μ1ω1, μ2ω2, Λ)
+    Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)    
+    U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
+    
+    W_right = get_W_matrix(basis_p, E_max, Λ, μ1ω1, μ2ω2; weights=false)
+    W_left = get_W_matrix(basis_p, E_max, Λ, μ1ω1, μ2ω2; weights=false)    
+
+    out = W_left * U * transpose(W_right)
+    return out
+end

berggren_3body.jl

@@ -42,17 +42,21 @@ end
     V1 = blockdiag(sparse.(V1_blocks)...)
 end
 
-E_max = 30
+@time "V_relative" begin
+    V_relative_blocks = [kron(get_V_matrix((k, kp) -> V_l(j1, k, kp), ks, ws), get_V_matrix((k, kp) -> V_l(j2, k, kp), ks, ws)) for (j1, j2) in js]
+    V_relative = blockdiag(sparse.(V_relative_blocks)...)
+end
+
+@time "weights" weights = spdiagm(vcat(fill(kron(ws, ws), length(js))...))
+
+# Moshinsky transform
+E_max = 40
 μω_global = 0.5
 μ1ω1 = μω_global * 1/2
 μ2ω2 = μω_global * 2
+@time "M" M = get_Moshinsky_in_kspace(E_max, basis, μ1ω1, μ2ω2, Λ)
+@time "V2" V2 = M' * weights * V_relative * M * weights
 
-@time "V2_HO" V2_HO =  get_jacobi_V2_matrix(V_of_r, E_max, Λ, μω_global)
-
-@time "W_right" W_right = get_W_matrix(basis, E_max, Λ, μ1ω1, μ2ω2; weights=true)
-@time "W_left" W_left = get_W_matrix(basis, E_max, Λ, μ1ω1, μ2ω2; weights=false)
-
-@time "V2" V2 = W_left * V2_HO * transpose(W_right)
 @time "H" H = T + V1 + V2
 
 @time "Eigenvalues" evals, _ = eigs(H, nev=3, ncv=24, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)