yapanuwan
/
BergEC-jl
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

struct for HO basis

This commit is contained in:
Nuwan Yapa 2024-12-10 19:41:15 -05:00
parent c0e2b0c910
commit 9f363d2ff1
14 changed files with 184 additions and 179 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
berggren_3body_resonance.jl
View File

@ -60,9 +60,11 @@ E_max = 30
μω = μω_global * 2
μ = m/2
@time "V12_HO" V12_HO = get_src_V12_matrix(V_of_r, E_max, Λ, μω_global; atol=10^-6, maxevals=10^5)
basis_ho = ho_basis_2B(E_max, Λ)
@time "W" W = get_W_matrix(basis, E_max, Λ, μω, μω; weights=true)
@time "V12_HO" V12_HO = get_src_V12_matrix(V_of_r, basis_ho, μω_global; atol=10^-6, maxevals=10^5)
@time "W" W = get_W_matrix(basis, basis_ho, μω, μω; weights=true)
@time "V12_p" V12_p = W * V12_HO * transpose(W)
@time "V12" V12 = transpose(U) * V12_p * U

18
calculations/2in3body_HO_B2R_EC.jl
View File

@ -17,23 +17,21 @@ E_max = 40
println("No of threads = ", Threads.nthreads())
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
l_max = max(maximum(l1s), maximum(l2s))
n_max = max(maximum(n1s), maximum(n2s))
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
mask2 = (n1s .== n1s') .&& (l1s .== l1s')
basis = ho_basis_2B(E_max, Λ)
l_max = max(maximum(basis.l1s), maximum(basis.l2s))
n_max = max(maximum(basis.n1s), maximum(basis.n2s))
println("Basis size = ", length(Es))
println("Basis size = ", basis.dim)
@time "T1" T1 = get_sp_T_matrix(n1s, l1s; mask=mask1, μω_gen=μ1ω1, μ=μ1)
@time "T2" T2 = get_sp_T_matrix(n2s, l2s; mask=mask2, μω_gen=μ2ω2, μ=μ2)
@time "T1" T1 = get_sp_T_matrix(basis.n1s, basis.l1s; mask=mask1(basis), μω_gen=μ1ω1, μ=μ1)
@time "T2" T2 = get_sp_T_matrix(basis.n2s, basis.l2s; mask=mask2(basis), μω_gen=μ2ω2, μ=μ2)
@time "V" Vb = get_jacobi_V_matrix(V_of_r, E_max, Λ, μ1ω1, μω_global)
@time "V" Vb = get_jacobi_V_matrix(V_of_r, basis, μ1ω1, μω_global)
@time "H0" Ha = T1 + T2
# free memory
Es = n1s = l1s = n2s = l2s = mask1 = mask2 = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
basis = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
GC.gc()
current_E = -0.26141959851000807

12
calculations/3body_Berggren_B2R_EC.jl
View File

@ -67,11 +67,13 @@ E_max = 40
μ1ω1 = μω_global * 1/2
μ2ω2 = μω_global * 2
@time "Va2_HO" Va2_HO = get_jacobi_V2_matrix(Va_of_r, E_max, Λ, μω_global; atol=atol, maxevals=maxevals)
@time "Vb2_HO" Vb2_HO = get_jacobi_V2_matrix(Vb_of_r, E_max, Λ, μω_global; atol=atol, maxevals=maxevals)
basis_ho = ho_basis_2B(E_max, Λ)
@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 "Va2_HO" Va2_HO = get_jacobi_V2_matrix(Va_of_r, basis_ho, μω_global; atol=atol, maxevals=maxevals)
@time "Vb2_HO" Vb2_HO = get_jacobi_V2_matrix(Vb_of_r, basis_ho, μω_global; atol=atol, maxevals=maxevals)
@time "W_right" W_right = get_W_matrix(basis, basis_ho, μ1ω1, μ2ω2; weights=true)
@time "W_left" W_left = get_W_matrix(basis, basis_ho, μ1ω1, μ2ω2; weights=false)
@time "Va2" Va2 = W_left * Va2_HO * transpose(W_right)
@time "Vb2" Vb2 = W_left * Vb2_HO * transpose(W_right)
@ -83,7 +85,7 @@ E_max = 40
display(test_evals)
# free memory
Es = n1s = l1s = n2s = l2s = mask1 = mask2 = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
basis = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
GC.gc()
current_E = training_ref

12
calculations/3body_Berggren_R2R_EC.jl
View File

@ -70,11 +70,13 @@ E_max = 40
μ1ω1 = μω_global * 1/2
μ2ω2 = μω_global * 2
@time "Va2_HO" Va2_HO = get_jacobi_V2_matrix(Va_of_r, E_max, Λ, μω_global; atol=atol, maxevals=maxevals)
@time "Vb2_HO" Vb2_HO = get_jacobi_V2_matrix(Vb_of_r, E_max, Λ, μω_global; atol=atol, maxevals=maxevals)
basis_ho = ho_basis_2B(E_max, Λ)
@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 "Va2_HO" Va2_HO = get_jacobi_V2_matrix(Va_of_r, basis_ho, μω_global; atol=atol, maxevals=maxevals)
@time "Vb2_HO" Vb2_HO = get_jacobi_V2_matrix(Vb_of_r, basis_ho, μω_global; atol=atol, maxevals=maxevals)
@time "W_right" W_right = get_W_matrix(basis, basis_ho, μ1ω1, μ2ω2; weights=true)
@time "W_left" W_left = get_W_matrix(basis, basis_ho, μ1ω1, μ2ω2; weights=false)
@time "Va2" Va2 = W_left * Va2_HO * transpose(W_right)
@time "Vb2" Vb2 = W_left * Vb2_HO * transpose(W_right)
@ -86,7 +88,7 @@ E_max = 40
display(test_evals)
# free memory
Es = n1s = l1s = n2s = l2s = mask1 = mask2 = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
basis = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
GC.gc()
exact = ComplexF64[]

20
calculations/3body_HO_B2R_EC.jl
View File

@ -18,19 +18,17 @@ E_max = 40
println("No of threads = ", Threads.nthreads())
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
l_max = max(maximum(l1s), maximum(l2s))
n_max = max(maximum(n1s), maximum(n2s))
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
mask2 = (n1s .== n1s') .&& (l1s .== l1s')
basis = ho_basis_2B(E_max, Λ)
l_max = max(maximum(basis.l1s), maximum(basis.l2s))
n_max = max(maximum(basis.n1s), maximum(basis.n2s))
println("Basis size = ", length(Es))
println("Basis size = ", basis.dim)
@time "T1" T1 = get_sp_T_matrix(n1s, l1s; mask=mask1, μω_gen=μ1ω1, μ=μ1)
@time "T2" T2 = get_sp_T_matrix(n2s, l2s; mask=mask2, μω_gen=μ2ω2, μ=μ2)
@time "T1" T1 = get_sp_T_matrix(basis.n1s, basis.l1s; mask=mask1(basis), μω_gen=μ1ω1, μ=μ1)
@time "T2" T2 = get_sp_T_matrix(basis.n2s, basis.l2s; mask=mask2(basis), μω_gen=μ2ω2, μ=μ2)
@time "Va" Va = get_jacobi_V_matrix(Va_of_r, E_max, Λ, μ1ω1, μω_global)
@time "Vb" Vb = get_jacobi_V_matrix(Vb_of_r, E_max, Λ, μ1ω1, μω_global)
@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)
@ -38,7 +36,7 @@ println("Basis size = ", length(Es))
display(target_evals)
# free memory
Es = n1s = l1s = n2s = l2s = mask1 = mask2 = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
basis = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
GC.gc()
current_E = -0.72763

2
calculations/3body_HO_R2R_EC.jl
View File

@ -9,7 +9,7 @@ extrapolating_c = 1.05 .- [0.0 : 0.1 : 0.4; 0.45 : 0.05 : 0.60]
@time "H0" H0 = T1 + T2
# free memory
Es = n1s = l1s = n2s = l2s = mask1 = mask2 = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
basis = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
GC.gc()
exact = ComplexF64[]

20
calculations/3body_dis_HO_EC.jl
View File

@ -20,19 +20,17 @@ E_max = 40
println("No of threads = ", Threads.nthreads())
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
l_max = max(maximum(l1s), maximum(l2s))
n_max = max(maximum(n1s), maximum(n2s))
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
mask2 = (n1s .== n1s') .&& (l1s .== l1s')
basis = ho_basis_2B(E_max, Λ)
l_max = max(maximum(basis.l1s), maximum(basis.l2s))
n_max = max(maximum(basis.n1s), maximum(basis.n2s))
println("Basis size = ", length(Es))
println("Basis size = ", basis.dim)
@time "T1" T1 = get_sp_T_matrix(n1s, l1s; mask=mask1, μω_gen=μ1ω1, μ=μ1)
@time "T2" T2 = get_sp_T_matrix(n2s, l2s; mask=mask2, μω_gen=μ2ω2, μ=μ2)
@time "T1" T1 = get_sp_T_matrix(basis.n1s, basis.l1s; mask=mask1(basis), μω_gen=μ1ω1, μ=μ1)
@time "T2" T2 = get_sp_T_matrix(basis.n2s, basis.l2s; mask=mask2(basis), μω_gen=μ2ω2, μ=μ2)
@time "Va" Va = get_jacobi_V1_matrix(Va_of_r, E_max, Λ, μ1ω1)
@time "Vb" Vb = get_jacobi_V2_matrix(Vb_of_r, E_max, Λ, μω_global)
@time "Va" Va = get_jacobi_V1_matrix(Va_of_r, basis, μ1ω1)
@time "Vb" Vb = get_jacobi_V2_matrix(Vb_of_r, basis, μω_global)
@time "Ha" Ha = T1 + T2 + Va
@time "Eigenvalues" target_evals, _ = eigs(Ha, nev=5, ncv=50, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
@ -40,7 +38,7 @@ println("Basis size = ", length(Es))
display(target_evals)
# free memory
Es = n1s = l1s = n2s = l2s = mask1 = mask2 = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
basis = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
GC.gc()
training_c = [-0.5, -0.65, -0.8, -1, -1.2]

133
ho_basis.jl
View File

@ -4,14 +4,14 @@ using LRUCache
include("helper.jl")
include("math.jl")
function V_numerical(V_of_r, l, n1, n2; μω_gen=1.0, atol=0, maxevals=10^7)
const_part = sqrt(μω_gen) * ho_basis_const(l, n1) * ho_basis_const(l, n2)
integrand(r) = ho_basis_func(l, n1, sqrt(μω_gen) * r) * ho_basis_func(l, n2, sqrt(μω_gen) * r) * V_of_r(r)
(integral, _) = quadgk(integrand, 0, Inf; atol=atol, maxevals=maxevals)
return const_part * integral
end
"1-body HO basis"
struct ho_basis_1B
dim::Int # dimensionality of the basis
Es::Vector{Int}
ns::Vector{Int}
ls::Vector{Int}
function get_sp_basis(E_max)
function ho_basis_1B(E_max)
Es = Int[]
ns = Int[]
ls = Int[]
@ -26,10 +26,21 @@ function get_sp_basis(E_max)
end
end
return (Es, ns, ls)
return new(length(Es), Es, ns, ls)
end
end
function get_2p_basis(E_max, Λ=-1)
"2-body HO basis"
struct ho_basis_2B
Λ::Int
dim::Int # dimensionality of the basis
Es::Vector{Int}
n1s::Vector{Int}
l1s::Vector{Int}
n2s::Vector{Int}
l2s::Vector{Int}
function ho_basis_2B(E_max, Λ=-1)
Es = Int[]
n1s = Int[]
l1s = Int[]
@ -53,7 +64,18 @@ function get_2p_basis(E_max, Λ=-1)
end
end
return (Es, n1s, l1s, n2s, l2s)
return new(Λ, length(Es), Es, n1s, l1s, n2s, l2s)
end
end
mask1(basis::ho_basis_2B) = (basis.n2s .== basis.n2s') .&& (basis.l2s .== basis.l2s')
mask2(basis::ho_basis_2B) = (basis.n1s .== basis.n1s') .&& (basis.l1s .== basis.l1s')
function V_numerical(V_of_r, l, n1, n2; μω_gen=1.0, atol=0, maxevals=10^7)
const_part = sqrt(μω_gen) * ho_basis_const(l, n1) * ho_basis_const(l, n2)
integrand(r) = ho_basis_func(l, n1, sqrt(μω_gen) * r) * ho_basis_func(l, n2, sqrt(μω_gen) * r) * V_of_r(r)
(integral, _) = quadgk(integrand, 0, Inf; atol=atol, maxevals=maxevals)
return const_part * integral
end
function get_sp_T_matrix(ns, ls; mask=trues(length(ns),length(ns)), μω_gen=1.0, μ=1.0)
@ -87,12 +109,12 @@ function get_sp_V_matrix(V_l, ns, ls; mask=trues(length(ns),length(ns)), dtype=F
return sparse(mat)
end
function Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
NQMAX = maximum(Es)
@assert all(mod.(Es, 2) .== mod(NQMAX, 2)) "Can only admit basis states with same parity"
function Moshinsky_transform(basis::ho_basis_2B)
NQMAX = maximum(basis.Es)
@assert all(mod.(basis.Es, 2) .== mod(NQMAX, 2)) "Can only admit basis states with same parity"
LMIN = Λ
LMAX = Λ
LMIN = basis.Λ
LMAX = basis.Λ
CO = 1/sqrt(2)
SI = 1/sqrt(2)
@ -100,15 +122,15 @@ function Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
BRAC = zeros(Float64, 1 + LMAX, 1 + (NQMAX - LMIN) ÷ 2, 1 + (NQMAX - LMIN) ÷ 2, 1 + (NQMAX - LMIN) ÷ 2, 1 + (NQMAX - LMIN) ÷ 2, 1 + LMAX, 1 + (NQMAX-LMIN) ÷ 2, 1 + LMAX-LMIN)
@ccall "../OSBRACKETS/allosbrac.so".allosbrac_(NQMAX::Ref{Int32},LMIN::Ref{Int32},LMAX::Ref{Int32},CO::Ref{Float64},SI::Ref{Float64},BRAC::Ptr{Array{Float64}})::Cvoid
mat = zeros(length(Es), length(Es))
mat = zeros(basis.dim, basis.dim)
s = hcat(Es, n1s, l1s, n2s, l2s)
s = hcat(basis.Es, basis.n1s, basis.l1s, basis.n2s, basis.l2s)
Threads.@threads for idx in CartesianIndices(mat)
(i, j) = Tuple(idx)
(Elhs, N, L, n, l) = s[i, :]
(Erhs, n1, l1, n2, l2) = s[j, :]
if Elhs == Erhs && triangle_ineq(L, l, Λ) && triangle_ineq(l1, l2, Λ)
mat[i, j] = (-1)^(n1 + n2 + N + n) * pick_Moshinsky_bracket(BRAC, n1, l1, n2, l2, N, L, n, l, Λ)
if Elhs == Erhs && triangle_ineq(L, l, basis.Λ) && triangle_ineq(l1, l2, basis.Λ)
mat[i, j] = (-1)^(n1 + n2 + N + n) * pick_Moshinsky_bracket(BRAC, n1, l1, n2, l2, N, L, n, l, basis.Λ)
end
end
@ -126,100 +148,89 @@ function pick_Moshinsky_bracket(BRAC, n1′, l1′, n2′, l2′, n1, l1, n2, l2
return BRAC[1 + NP, 1 + n1′, 1 + MP, 1 + n1, 1 + n2, 1 + N, 1 + M, 1]
end
function get_jacobi_V_matrix(V_of_r, E_max, Λ, μ1ω1, μω_global; atol=10^-6, maxevals=10^5)
V1 = get_jacobi_V1_matrix(V_of_r, E_max, Λ, μ1ω1; atol=atol, maxevals=maxevals)
V2 = get_jacobi_V2_matrix(V_of_r, E_max, Λ, μω_global; atol=atol, maxevals=maxevals)
function get_jacobi_V_matrix(V_of_r, basis::ho_basis_2B, μ1ω1, μω_global; atol=10^-6, maxevals=10^5)
V1 = get_jacobi_V1_matrix(V_of_r, basis, μ1ω1; atol=atol, maxevals=maxevals)
V2 = get_jacobi_V2_matrix(V_of_r, basis, μω_global; atol=atol, maxevals=maxevals)
return V1 + V2
end
function get_jacobi_V1_matrix(V_of_r, E_max, Λ, μ1ω1; atol=10^-6, maxevals=10^5)
_, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
l_max = max(maximum(l1s), maximum(l2s))
n_max = max(maximum(n1s), maximum(n2s))
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
function get_jacobi_V1_matrix(V_of_r, basis::ho_basis_2B, μ1ω1; atol=10^-6, maxevals=10^5)
l_max = max(maximum(basis.l1s), maximum(basis.l2s))
n_max = max(maximum(basis.n1s), maximum(basis.n2s))
V1_elem(l, n1, n2) = V_numerical(V_of_r, l, n1, n2; μω_gen=μ1ω1, atol=atol, maxevals=maxevals)
V1_cache = LRU{Tuple{UInt8, UInt8, UInt8}, ComplexF64}(maxsize=(1+l_max)*(1+n_max)^2)
V1 = get_sp_V_matrix(V1_elem, n1s, l1s; mask=mask1, dtype=ComplexF64, cache=V1_cache)
V1 = get_sp_V_matrix(V1_elem, basis.n1s, basis.l1s; mask=mask1(basis), dtype=ComplexF64, cache=V1_cache)
return V1
end
function get_jacobi_V2_matrix(V_of_r, E_max, Λ, μω_global; atol=10^-6, maxevals=10^5)
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
l_max = max(maximum(l1s), maximum(l2s))
n_max = max(maximum(n1s), maximum(n2s))
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
mask2 = (n1s .== n1s') .&& (l1s .== l1s')
function get_jacobi_V2_matrix(V_of_r, basis::ho_basis_2B, μω_global; atol=10^-6, maxevals=10^5)
l_max = max(maximum(basis.l1s), maximum(basis.l2s))
n_max = max(maximum(basis.n1s), maximum(basis.n2s))
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}, ComplexF64}(maxsize=(1+l_max)*(1+n_max)^2)
V_relative = get_sp_V_matrix(V_relative_elem, n1s, l1s; mask=mask1, dtype=ComplexF64, cache=V_relative_cache) + get_sp_V_matrix(V_relative_elem, n2s, l2s; mask=mask2, dtype=ComplexF64, cache=V_relative_cache)
U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
V_relative = get_sp_V_matrix(V_relative_elem, basis.n1s, basis.l1s; mask=mask1(basis), dtype=ComplexF64, cache=V_relative_cache) + get_sp_V_matrix(V_relative_elem, basis.n2s, basis.l2s; mask=mask2(basis), dtype=ComplexF64, cache=V_relative_cache)
U = Moshinsky_transform(basis)
V2 = U' * V_relative * U
return V2
end
function get_2p_p1p2_matrix(n1s, l1s, n2s, l2s, Λ, μ1ω1, μ2ω2; dtype=Float64)
function get_2p_p1p2_matrix(basis::ho_basis_2B, μ1ω1, μ2ω2; dtype=Float64)
# TODO: Cache for integrals
integral1(np, lp, n, l) = integral_HO(np, lp, n, l, μ1ω1)
integral2(np, lp, n, l) = integral_HO(np, lp, n, l, μ2ω2)
mat = zeros(dtype, length(n1s), length(n1s))
mat = zeros(dtype, basis.dim, basis.dim)
Threads.@threads for idx in CartesianIndices(mat)
(i, j) = Tuple(idx)
val = racahs_reduction_formula(n1s[i], l1s[i], n2s[i], l2s[i], n1s[j], l1s[j], n2s[j], l2s[j], Λ, integral1, integral2)
val = racahs_reduction_formula(basis.n1s[i], basis.l1s[i], basis.n2s[i], basis.l2s[i], basis.n1s[j], basis.l1s[j], basis.n2s[j], basis.l2s[j], basis.Λ, integral1, integral2)
if !(val ≈ 0); mat[idx] = val; end
end
return sparse(mat)
end
function get_src_V_matrix(V_of_r, E_max, Λ, μω, μω_global; atol=10^-6, maxevals=10^5)
_, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
l_max = max(maximum(l1s), maximum(l2s))
n_max = max(maximum(n1s), maximum(n2s))
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
mask2 = (n1s .== n1s') .&& (l1s .== l1s')
function get_src_V_matrix(V_of_r, basis::ho_basis_2B, μω, μω_global; atol=10^-6, maxevals=10^5)
l_max = max(maximum(basis.l1s), maximum(basis.l2s))
n_max = max(maximum(basis.n1s), maximum(basis.n2s))
V_elem(l, n1, n2) = V_numerical(V_of_r, l, n1, n2; μω_gen=μω, atol=atol, maxevals=maxevals)
V_cache = LRU{Tuple{UInt8, UInt8, UInt8}, ComplexF64}(maxsize=(1+l_max)*(1+n_max)^2)
V1 = get_sp_V_matrix(V_elem, n1s, l1s; mask=mask1, dtype=ComplexF64, cache=V_cache)
V2 = get_sp_V_matrix(V_elem, n2s, l2s; mask=mask2, dtype=ComplexF64, cache=V_cache)
V1 = get_sp_V_matrix(V_elem, basis.n1s, basis.l1s; mask=mask1(basis), dtype=ComplexF64, cache=V_cache)
V2 = get_sp_V_matrix(V_elem, basis.n2s, basis.l2s; mask=mask2(basis), dtype=ComplexF64, cache=V_cache)
V12 = get_src_V12_matrix(V_of_r, E_max, Λ, μω_global; atol=atol, maxevals=maxevals)
V12 = get_src_V12_matrix(V_of_r, basis, μω_global; atol=atol, maxevals=maxevals)
return V1 + V2 + V12
end
function get_src_V12_matrix(V_of_r, E_max, Λ, μω_global; atol=10^-6, maxevals=10^5)
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
l_max = max(maximum(l1s), maximum(l2s))
n_max = max(maximum(n1s), maximum(n2s))
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
function get_src_V12_matrix(V_of_r, basis::ho_basis_2B, μω_global; atol=10^-6, maxevals=10^5)
l_max = max(maximum(basis.l1s), maximum(basis.l2s))
n_max = max(maximum(basis.n1s), maximum(basis.n2s))
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}, ComplexF64}(maxsize=(1+l_max)*(1+n_max)^2)
V_relative = get_sp_V_matrix(V_relative_elem, n1s, l1s; mask=mask1, dtype=ComplexF64, cache=V_relative_cache)
U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
V_relative = get_sp_V_matrix(V_relative_elem, basis.n1s, basis.l1s; mask=mask1(basis), dtype=ComplexF64, cache=V_relative_cache)
U = Moshinsky_transform(basis)
V12 = U' * V_relative * U
return V12
end
"Basis transformation from HO to momentum space"
function get_W_matrix(basis_p, E_max, Λ, μ1ω1, μ2ω2=μ1ω1; weights=true)
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
W = zeros(ComplexF64, length(basis_p), length(Es))
function get_W_matrix(basis_p, basis::ho_basis_2B, μ1ω1, μ2ω2=μ1ω1; weights=true)
W = zeros(ComplexF64, length(basis_p), basis.dim)
Threads.@threads for idx in CartesianIndices(W)
(i1, i2) = Tuple(idx)
((j1, j2), (k1, w1), (k2, w2)) = basis_p[i1]
if j1 == l1s[i2] && j2 == l2s[i2]
elem1 = 1/sqrt(sqrt(μ1ω1)) * (-1)^n1s[i2] * ho_basis(j1, n1s[i2], 1/sqrt(μ1ω1) * k1)
elem2 = 1/sqrt(sqrt(μ2ω2)) * (-1)^n2s[i2] * ho_basis(j2, n2s[i2], 1/sqrt(μ2ω2) * k2)
if j1 == basis.l1s[i2] && j2 == basis.l2s[i2]
elem1 = 1/sqrt(sqrt(μ1ω1)) * (-1)^basis.n1s[i2] * ho_basis(j1, basis.n1s[i2], 1/sqrt(μ1ω1) * k1)
elem2 = 1/sqrt(sqrt(μ2ω2)) * (-1)^basis.n2s[i2] * ho_basis(j2, basis.n2s[i2], 1/sqrt(μ2ω2) * k2)
W[idx] = elem1 * elem2 * (weights ? w1 * w2 : 1)
end
end

8
ho_basis_2body.jl
View File

@ -13,15 +13,15 @@ Ra = 2
println("No of threads = ", Threads.nthreads())
Es, n1s, l1s = get_sp_basis(E_max)
println("Basis size = ", length(Es))
basis = ho_basis_1B(E_max)
println("Basis size = ", basis.dim)
println("Constructing KE matrices")
@time "T1" T1 = get_sp_T_matrix(n1s, l1s; μω_gen=μω_gen, μ=μ1)
@time "T1" T1 = get_sp_T_matrix(basis.ns, basis.ls; μω_gen=μω_gen, μ=μ1)
println("Constructing PE matrices")
V1_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; μω_gen=μω_gen)
@time "V1" V1 = get_sp_V_matrix(V1_elem, n1s, l1s)
@time "V1" V1 = get_sp_V_matrix(V1_elem, basis.ns, basis.ls)
println("Calculating spectrum")
@time "H" H = T1 + V1

22
ho_basis_3body.jl
View File

@ -15,26 +15,22 @@ E_max = 40
println("No of threads = ", Threads.nthreads())
@time "Basis" begin
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
mask2 = (n1s .== n1s') .&& (l1s .== l1s')
end
@time "Basis" basis = ho_basis_2B(E_max, Λ)
println("Basis size = ", length(Es))
println("Basis size = ", basis.dim)
println("Constructing KE matrices")
@time "T1" T1 = get_sp_T_matrix(n1s, l1s; mask=mask1, μω_gen=μω, μ=μ)
@time "T2" T2 = get_sp_T_matrix(n2s, l2s; mask=mask2, μω_gen=μω, μ=μ)
@time "T_cross" T_cross = get_2p_p1p2_matrix(n1s, l1s, n2s, l2s, Λ, μω, μω) ./ (2*μ)
@time "T1" T1 = get_sp_T_matrix(basis.n1s, basis.l1s; mask=mask1(basis), μω_gen=μω, μ=μ)
@time "T2" T2 = get_sp_T_matrix(basis.n2s, basis.l2s; mask=mask2(basis), μω_gen=μω, μ=μ)
@time "T_cross" T_cross = get_2p_p1p2_matrix(basis, μω, μω) ./ (2*μ)
println("Constructing PE matrices")
V_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; μω_gen=μω)
V_relative_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; μω_gen=μω_global)
@time "V1" V1 = get_sp_V_matrix(V_elem, n1s, l1s; mask=mask1)
@time "V2" V2 = get_sp_V_matrix(V_elem, n2s, l2s; mask=mask2)
@time "V relative" V_relative = get_sp_V_matrix(V_relative_elem, n1s, l1s; mask=mask1)
@time "Moshinsky brackets" U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
@time "V1" V1 = get_sp_V_matrix(V_elem, basis.n1s, basis.l1s; mask=mask1(basis))
@time "V2" V2 = get_sp_V_matrix(V_elem, basis.n2s, basis.l2s; mask=mask2(basis))
@time "V relative" V_relative = get_sp_V_matrix(V_relative_elem, basis.n1s, basis.l1s; mask=mask1(basis))
@time "Moshinsky brackets" U = Moshinsky_transform(basis)
@time "V12" V12 = U' * V_relative * U
println("Calculating spectrum")

18
ho_basis_3body_resonance.jl
View File

@ -14,25 +14,19 @@ E_max = 30
println("No of threads = ", Threads.nthreads())
@time "Basis" begin
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
l_max = max(maximum(l1s), maximum(l2s))
n_max = max(maximum(n1s), maximum(n2s))
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
mask2 = (n1s .== n1s') .&& (l1s .== l1s')
end
@time "Basis" basis = ho_basis_2B(E_max, Λ)
println("Basis size = ", length(Es))
println("Basis size = ", basis.dim)
println("Constructing KE matrices")
@time "T1" T1 = get_sp_T_matrix(n1s, l1s; mask=mask1, μω_gen=μω, μ=μ)
@time "T2" T2 = get_sp_T_matrix(n2s, l2s; mask=mask2, μω_gen=μω, μ=μ)
@time "T_cross" T_cross = get_2p_p1p2_matrix(n1s, l1s, n2s, l2s, Λ, μω, μω; dtype=ComplexF64) ./ (2*μ)
@time "T1" T1 = get_sp_T_matrix(basis.n1s, basis.l1s; mask=mask1(basis), μω_gen=μω, μ=μ)
@time "T2" T2 = get_sp_T_matrix(basis.n2s, basis.l2s; mask=mask2(basis), μω_gen=μω, μ=μ)
@time "T_cross" T_cross = get_2p_p1p2_matrix(basis, μω, μω; dtype=ComplexF64) ./ (2*μ)
println("Constructing PE matrices")
@time "V" V = get_src_V_matrix(V_of_r, E_max, Λ, μω, μω_global)
@time "V" V = get_src_V_matrix(V_of_r, basis, μω, μω_global)
println("Calculating spectrum")
@time "H" H = T1 + T2 + T_cross + V

8
p_space_3body.jl
View File

@ -42,10 +42,12 @@ E_max = 30
μ1ω1 = μω_global * 1/2
μ2ω2 = μω_global * 2
@time "V2_HO" V2_HO = get_jacobi_V2_matrix(V_of_r, E_max, Λ, μω_global)
basis_ho = ho_basis_2B(E_max, Λ)
@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_HO" V2_HO = get_jacobi_V2_matrix(V_of_r, basis_ho, μω_global)
@time "W_right" W_right = get_W_matrix(basis, basis_ho, μ1ω1, μ2ω2; weights=true)
@time "W_left" W_left = get_W_matrix(basis, basis_ho, μ1ω1, μ2ω2; weights=false)
@time "V2" V2 = W_left * V2_HO * transpose(W_right)
@time "H" H = Hb + V2

8
p_space_3body_resonance.jl
View File

@ -41,10 +41,12 @@ E_max = 30
μ1ω1 = μω_global * 1/2
μ2ω2 = μω_global * 2
@time "V2_HO" V2_HO = get_jacobi_V2_matrix(V_of_r, E_max, Λ, μω_global)
basis_ho = ho_basis_2B(E_max, Λ)
@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_HO" V2_HO = get_jacobi_V2_matrix(V_of_r, basis_ho, μω_global)
@time "W_right" W_right = get_W_matrix(basis, basis_ho, μ1ω1, μ2ω2; weights=true)
@time "W_left" W_left = get_W_matrix(basis, basis_ho, μ1ω1, μ2ω2; weights=false)
@time "V2" V2 = W_left * V2_HO * transpose(W_right)
@time "H" H = Hb + V2

6
test/moshinsky.jl
View File

@ -9,11 +9,11 @@ E_max = 30
println("No of threads = ", Threads.nthreads())
@time "Basis" Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
@time "Basis" basis = ho_basis_2B(E_max, Λ)
println("Basis size = ", length(Es))
println("Basis size = ", basis.dim)
@time "Moshinsky brackets" U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
@time "Moshinsky brackets" U = Moshinsky_transform(basis)
check = U' * U - I