From 3bdf6285e27786b7601d9e9ce3c8f4e211242c37 Mon Sep 17 00:00:00 2001
From: ysyapa <ysyapa@ncsu.edu>
Date: Wed, 1 May 2024 17:46:08 +0000
Subject: [PATCH] V2 not working

---
 berggren_3body.jl | 98 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 98 insertions(+)
 create mode 100644 berggren_3body.jl

diff --git a/berggren_3body.jl b/berggren_3body.jl
new file mode 100644
index 0000000..5c5a283
--- /dev/null
+++ b/berggren_3body.jl
@@ -0,0 +1,98 @@
+using LinearAlgebra, SparseArrays, Arpack
+include("helper.jl")
+include("p_space.jl")
+include("ho_basis.jl")
+
+println("No of threads = ", Threads.nthreads())
+
+atol = 10^-5
+maxevals = 10^5
+
+# E_target = -0.3919
+
+m = 1.0
+μ1 = m * 1/2
+μ2 = m * 2/3
+
+Va = -2
+Ra = 2
+V_of_r(r) = Va * exp(-r^2 / Ra^2)
+V_l(j, k, kp) = Vl_mat_elem(V_of_r, j, k, kp; atol=atol, maxevals=maxevals, R_cutoff=16)
+
+vertices = [0, 0.5 - 0.3im, 1, 4]
+subdivisions = [10, 10, 10]
+ks, ws = get_mesh(vertices, subdivisions)
+js = collect(0:4)
+
+sp_basis_size = length(ks) * length(js)
+sp_ws = spdiagm(vcat(fill(ws, length(js))...))
+
+basis = collect(Iterators.product(zip(ks, ws), js, zip(ks, ws), js))[:]
+for ((k1, w1), j1, (k2, w2), j2) in basis
+    @assert k1 ∈ ks && k2 ∈ ks && w1 ∈ ws && w2 ∈ ws && j1 ∈ js && j2 ∈ js "Something wrong with the basis ordering"
+end
+println("Basis size = $(length(basis))")
+
+@time "T1" begin
+    sp_T1_j = [get_T_matrix(ks, μ1) for j in js]
+    sp_T1 = blockdiag(sparse.(sp_T1_j)...)
+    T1 = kron(sp_T1, I(sp_basis_size))
+end
+
+@time "T2" begin
+    sp_T2_j = [get_T_matrix(ks, μ2) for j in js]
+    sp_T2 = blockdiag(sparse.(sp_T2_j)...)
+    T2 = kron(I(sp_basis_size), sp_T2)
+end
+
+@time "V1" begin
+    sp_V1_j = [get_V_matrix((k, kp) -> V_l(j, k, kp), ks, ws) for j in js]
+    sp_V1 = blockdiag(sparse.(sp_V1_j)...)
+    V1 = kron(sp_V1, I(sp_basis_size))
+end
+
+#############################################################################################
+
+function get_W_matrix(basis, basis_HO, μ1ω1, μ2ω2=μ1ω1)
+    Es, n1s, l1s, n2s, l2s = basis_HO
+    W = zeros(ComplexF64, length(basis), length(Es))
+    Threads.@threads for idx in CartesianIndices(W)
+        (i1, i2) = Tuple(idx)
+        ((k1, w1), j1, (k2, w2), j2) = basis[i1]
+        if j1 == l1s[i2] && j2 == l2s[i2]
+            elem1 = 1/sqrt(sqrt(μ1ω1)) * ho_basis(j1, n1s[i2], 1/sqrt(μ1ω1) * k1) * w1
+            elem2 = 1/sqrt(sqrt(μ2ω2)) * ho_basis(j2, n2s[i2], 1/sqrt(μ2ω2) * k2) * w2
+            W[idx] = elem1 * elem2
+        end
+    end
+    return sparse(W)
+end
+
+############################################################################################
+
+Λ = 0
+E_max = 24
+μω_global = 0.5
+
+μ1ω1 = μω_global * 1/2
+μ2ω2 = μω_global * 2
+
+basis_HO = get_2p_basis(E_max)
+Es, n1s, l1s, n2s, l2s = basis_HO
+l_max = max(maximum(l1s), maximum(l2s))
+n_max = max(maximum(n1s), maximum(n2s))
+mask1 = (n2s .== n2s') .&& (l2s .== l2s')
+mask2 = (n1s .== n1s') .&& (l1s .== l1s')
+
+V_relative_elem(l, n1, n2) = V_numerical(V_of_r, l, n1, n2; μω_gen=μω_global, atol=atol, maxevals=maxevals)
+V_relative_cache = LRU{Tuple{UInt8, UInt8, UInt8}, Float64}(maxsize=(1+l_max)*(1+n_max)^2)
+@time "V relative" V_relative = sp_V_matrix(V_relative_elem, n1s, l1s; mask=mask1, dtype=Float64, cache=V_relative_cache) + sp_V_matrix(V_relative_elem, n2s, l2s; mask=mask2, dtype=Float64, cache=V_relative_cache)
+@time "Moshinsky brackets" U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
+@time "V2_HO" V2_HO =  U' * V_relative * U
+
+@time "W" W = get_W_matrix(basis, basis_HO, μ1ω1, μ2ω2)
+@time "V2" V2 = W * V2_HO * transpose(W)
+@time "H" H = T1 + V1 + T2 + V2
+
+@time "Eigenvalues" evals, _ = eigs(H, nev=3, ncv=24, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
+display(evals)