diff --git a/berggren_3body_resonance.jl b/berggren_3body_resonance.jl index a6c46b9..2c98cc4 100644 --- a/berggren_3body_resonance.jl +++ b/berggren_3body_resonance.jl @@ -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 diff --git a/calculations/2in3body_HO_B2R_EC.jl b/calculations/2in3body_HO_B2R_EC.jl index ac1691d..2e103a3 100644 --- a/calculations/2in3body_HO_B2R_EC.jl +++ b/calculations/2in3body_HO_B2R_EC.jl @@ -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 diff --git a/calculations/3body_Berggren_B2R_EC.jl b/calculations/3body_Berggren_B2R_EC.jl index 9824e42..23675f5 100644 --- a/calculations/3body_Berggren_B2R_EC.jl +++ b/calculations/3body_Berggren_B2R_EC.jl @@ -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 diff --git a/calculations/3body_Berggren_R2R_EC.jl b/calculations/3body_Berggren_R2R_EC.jl index d147e4b..1b2c6bc 100644 --- a/calculations/3body_Berggren_R2R_EC.jl +++ b/calculations/3body_Berggren_R2R_EC.jl @@ -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[] diff --git a/calculations/3body_HO_B2R_EC.jl b/calculations/3body_HO_B2R_EC.jl index f08406b..246a221 100644 --- a/calculations/3body_HO_B2R_EC.jl +++ b/calculations/3body_HO_B2R_EC.jl @@ -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 diff --git a/calculations/3body_HO_R2R_EC.jl b/calculations/3body_HO_R2R_EC.jl index d71e3bd..f560042 100644 --- a/calculations/3body_HO_R2R_EC.jl +++ b/calculations/3body_HO_R2R_EC.jl @@ -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[] diff --git a/calculations/3body_dis_HO_EC.jl b/calculations/3body_dis_HO_EC.jl index d123813..90d24c1 100644 --- a/calculations/3body_dis_HO_EC.jl +++ b/calculations/3body_dis_HO_EC.jl @@ -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] diff --git a/ho_basis.jl b/ho_basis.jl index 5e4bd99..9873407 100644 --- a/ho_basis.jl +++ b/ho_basis.jl @@ -4,6 +4,73 @@ using LRUCache include("helper.jl") include("math.jl") +"1-body HO basis" +struct ho_basis_1B + dim::Int # dimensionality of the basis + Es::Vector{Int} + ns::Vector{Int} + ls::Vector{Int} + + function ho_basis_1B(E_max) + Es = Int[] + ns = Int[] + ls = Int[] + + # E = 2*n + l + for E in 0 : E_max + for n in 0 : E ÷ 2 + l = E - 2*n + push!(Es, E) + push!(ns, n) + push!(ls, l) + end + end + + return new(length(Es), Es, ns, ls) + end +end + +"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[] + n2s = Int[] + l2s = Int[] + + # E = 2*n1 + l1 + 2*n2 + l2 + for E in E_max : -2 : 0 # same parity states only + for n1 in 0 : E ÷ 2 + for n2 in 0 : (E - 2*n1) ÷ 2 + for l1 in 0 : (E - 2*n1 - 2*n2) + l2 = E - 2*n1 - 2*n2 - l1 + if Λ≥0 && !triangle_ineq(l1, l2, Λ); continue; end + push!(Es, E) + push!(n1s, n1) + push!(l1s, l1) + push!(n2s, n2) + push!(l2s, l2) + end + end + end + end + + 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) @@ -11,51 +78,6 @@ function V_numerical(V_of_r, l, n1, n2; μω_gen=1.0, atol=0, maxevals=10^7) return const_part * integral end -function get_sp_basis(E_max) - Es = Int[] - ns = Int[] - ls = Int[] - - # E = 2*n + l - for E in 0 : E_max - for n in 0 : E ÷ 2 - l = E - 2*n - push!(Es, E) - push!(ns, n) - push!(ls, l) - end - end - - return (Es, ns, ls) -end - -function get_2p_basis(E_max, Λ=-1) - Es = Int[] - n1s = Int[] - l1s = Int[] - n2s = Int[] - l2s = Int[] - - # E = 2*n1 + l1 + 2*n2 + l2 - for E in E_max : -2 : 0 # same parity states only - for n1 in 0 : E ÷ 2 - for n2 in 0 : (E - 2*n1) ÷ 2 - for l1 in 0 : (E - 2*n1 - 2*n2) - l2 = E - 2*n1 - 2*n2 - l1 - if Λ≥0 && !triangle_ineq(l1, l2, Λ); continue; end - push!(Es, E) - push!(n1s, n1) - push!(l1s, l1) - push!(n2s, n2) - push!(l2s, l2) - end - end - end - end - - return (Es, n1s, l1s, n2s, l2s) -end - function get_sp_T_matrix(ns, ls; mask=trues(length(ns),length(ns)), μω_gen=1.0, μ=1.0) mat = spzeros(length(ns), length(ns)) for idx in CartesianIndices(mat) @@ -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 diff --git a/ho_basis_2body.jl b/ho_basis_2body.jl index f8e704f..80b9d62 100644 --- a/ho_basis_2body.jl +++ b/ho_basis_2body.jl @@ -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 diff --git a/ho_basis_3body.jl b/ho_basis_3body.jl index 15d39d3..9053573 100644 --- a/ho_basis_3body.jl +++ b/ho_basis_3body.jl @@ -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") diff --git a/ho_basis_3body_resonance.jl b/ho_basis_3body_resonance.jl index 8cf1066..1c775c2 100644 --- a/ho_basis_3body_resonance.jl +++ b/ho_basis_3body_resonance.jl @@ -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 diff --git a/p_space_3body.jl b/p_space_3body.jl index 1dfba2d..f92eccc 100644 --- a/p_space_3body.jl +++ b/p_space_3body.jl @@ -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 diff --git a/p_space_3body_resonance.jl b/p_space_3body_resonance.jl index 7c9cf37..81264d9 100644 --- a/p_space_3body_resonance.jl +++ b/p_space_3body_resonance.jl @@ -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 diff --git a/test/moshinsky.jl b/test/moshinsky.jl index 515a108..85019dc 100644 --- a/test/moshinsky.jl +++ b/test/moshinsky.jl @@ -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