From 12c09caa9562f3b742f07597b686a7e8fe596c7f Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Mon, 23 Dec 2024 17:03:55 -0500
Subject: [PATCH] Simplification

---
 calculations/3body_HO_B2R_EC.jl | 66 ++++++++++-----------------------
 1 file changed, 20 insertions(+), 46 deletions(-)

diff --git a/calculations/3body_HO_B2R_EC.jl b/calculations/3body_HO_B2R_EC.jl
index a10e2cc..54ac2ab 100644
--- a/calculations/3body_HO_B2R_EC.jl
+++ b/calculations/3body_HO_B2R_EC.jl
@@ -1,55 +1,29 @@
-using Arpack, SparseArrays, LRUCache
-using DelimitedFiles, Plots
-include("../ho_basis.jl")
+using Plots
 
-Λ = 0
-m = 1.0
-Va_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
-Vb_of_r(r) = -exp(-(r/3)^2)
-
-E_max = 40
-μω_global = 0.5 * exp(-2im * pi / 9)
-
-# due to Jacobi coordinates
-μ1ω1 = μω_global * 1/2
-μ2ω2 = μω_global * 2
-μ1 = m * 1/2
-μ2 = m * 2/3
-
-println("No of threads = ", Threads.nthreads())
-
-basis = ho_basis_2B(E_max, Λ)
-println("Basis size = ", basis.dim)
-
-@time "T1" T1 = get_sp_T_matrix(basis.n1s, basis.l1s, [basis.n2s, basis.l2s]; μω_gen=μ1ω1, μ=μ1)
-@time "T2" T2 = get_sp_T_matrix(basis.n2s, basis.l2s, [basis.n1s, basis.l1s]; μω_gen=μ2ω2, μ=μ2)
-
-@time "Va" Va = get_jacobi_V_matrix(Va_of_r, basis, μ1ω1, μω_global)
-@time "Vb" Vb = get_jacobi_V_matrix(Vb_of_r, basis, μ1ω1, μω_global)
-
-@time "Ha" Ha = T1 + T2 + Va
-@time "Eigenvalues" target_evals, _ = eigs(Ha, nev=5, ncv=50, which=:LI, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
-
-display(target_evals)
-
-# free memory
-basis = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
-GC.gc()
-
-current_E = -0.72763
+training_ref = -0.72763
+exact_ref = 4.0766890719636635 - 0.01275892774109674im
 
 training_c = [2.0, 1.9, 1.8]
 extrapolating_c = 0.0 : 0.2 : 1.2
 
+include("../ho_basis_3body_resonance.jl")
+H0 = H
+
+# Vp = perturbation to make the state artificially bound
+Vp_of_r(r) = -exp(-(r/3)^2)
+@time "Vp" Vp =  get_src_V_matrix(Vp_of_r, basis, μω, μω_global)
+
 exact = ComplexF64[]
 training = ComplexF64[]
 extrapolated = ComplexF64[]
 training_vecs = Vector{ComplexF64}[]
 
+current_E = training_ref
+
 for c in training_c
     println("Training for c = $c")
-    H = Ha + c .* Vb
-    evals, evecs = eigs(H, nev=3, ncv=24, which=:LI, maxiter=5000, tol=1e-5, ritzvec=true, check=1)
+    local H = H0 + c .* Vp
+    local evals, evecs = eigs(H, nev=3, ncv=24, which=:LI, maxiter=5000, tol=1e-5, ritzvec=true, check=1)
 
     global current_E = nearest(evals, current_E)
     push!(training, current_E)
@@ -61,21 +35,21 @@ training_vecs = vcat(training_vecs, conj(training_vecs))
 
 EC_basis = hcat(training_vecs...)
 N_EC = transpose(EC_basis) * EC_basis
-Ha_EC = transpose(EC_basis) * Ha * EC_basis
-Vb_EC = transpose(EC_basis) * Vb * EC_basis
+H0_EC = transpose(EC_basis) * H0 * EC_basis
+Vp_EC = transpose(EC_basis) * Vp * EC_basis
 
-current_E = 4.0766890719636635 - 0.01275892774109674im
+current_E = exact_ref
 
 for c in extrapolating_c
     println("Extrapolating for c = $c")
-    H = Ha + c .* Vb
-    evals, evecs = eigs(H, nev=3, ncv=24, which=:LI, maxiter=5000, tol=1e-5, ritzvec=true, check=1)
+    local H = H0 + c .* Vp
+    local evals, _ = eigs(H, nev=3, ncv=24, which=:LI, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
 
     global current_E = nearest(evals, current_E)
     push!(exact, current_E)
 
     # extrapolation
-    H_EC = Ha_EC + c .* Vb_EC
+    H_EC = H0_EC + c .* Vp_EC
     evals = eigvals(H_EC, N_EC)
     push!(extrapolated, nearest(evals, current_E))
 end