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

Critical bug fixed: missing Kronecker deltas

This commit is contained in:
Nuwan Yapa 2024-04-02 18:48:41 -04:00
parent ccee9c38ea
commit fbab9c0316
2 changed files with 10 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
ho_basis.jl
View File

@ -54,9 +54,10 @@ function get_2p_basis(E_max)
return (Es, n1s, l1s, n2s, l2s) return (Es, n1s, l1s, n2s, l2s)
end end
function sp_T_matrix(ns, ls; ω=1.0, μ=1.0) function sp_T_matrix(ns, ls; mask=trues(length(ns),length(ns)), ω=1.0, μ=1.0)
mat = spzeros(length(ns), length(ns)) mat = spzeros(length(ns), length(ns))
for idx in CartesianIndices(mat) for idx in CartesianIndices(mat)
if !mask[idx]; continue; end
(i, j) = Tuple(idx) (i, j) = Tuple(idx)
if ls[i] == ls[j] if ls[i] == ls[j]
if ns[i] == ns[j] if ns[i] == ns[j]
@ -70,9 +71,10 @@ function sp_T_matrix(ns, ls; ω=1.0, μ=1.0)
return (ω / μ) .* mat return (ω / μ) .* mat
end end
function sp_V_matrix(V_l, ns, ls; dtype=Float64) function sp_V_matrix(V_l, ns, ls; mask=trues(length(ns),length(ns)), dtype=Float64)
mat = zeros(dtype, length(ns), length(ns)) mat = zeros(dtype, length(ns), length(ns))
Threads.@threads for idx in CartesianIndices(mat) Threads.@threads for idx in CartesianIndices(mat)
if !mask[idx]; continue; end
(i, j) = Tuple(idx) (i, j) = Tuple(idx)
if ls[i] == ls[j] if ls[i] == ls[j]
mat[idx] = V_l(ls[i], ns[i], ns[j]) mat[idx] = V_l(ls[i], ns[i], ns[j])

10
ho_basis_3body.jl
View File

@ -18,17 +18,19 @@ c2 = 2
println("No of threads = ", Threads.nthreads()) println("No of threads = ", Threads.nthreads())
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max) Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max)
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
mask2 = (n1s .== n1s') .&& (l1s .== l1s')
println("Basis size = ", length(Es)) println("Basis size = ", length(Es))
println("Constructing KE matrices") println("Constructing KE matrices")
@time "T1" T1 = sp_T_matrix(n1s, l1s; ω=ω, μ=μ1) @time "T1" T1 = sp_T_matrix(n1s, l1s; mask=mask1, ω=ω, μ=μ1)
@time "T2" T2 = sp_T_matrix(n2s, l2s; ω=ω, μ=c2^2 * μ2) @time "T2" T2 = sp_T_matrix(n2s, l2s; mask=mask2, ω=ω, μ=c2^2 * μ2)
println("Constructing PE matrices") println("Constructing PE matrices")
V1_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; ω=ω) V1_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; ω=ω)
V_relative_elem(l, n1, n2) = Va * V_Gaussian(Ra / c, l, n1, n2; ω=ω) V_relative_elem(l, n1, n2) = Va * V_Gaussian(Ra / c, l, n1, n2; ω=ω)
@time "V1" V1 = sp_V_matrix(V1_elem, n1s, l1s) @time "V1" V1 = sp_V_matrix(V1_elem, n1s, l1s; mask=mask1)
@time "V relative" V_relative = sp_V_matrix(V_relative_elem, n1s, l1s) + sp_V_matrix(V_relative_elem, n2s, l2s) @time "V relative" V_relative = sp_V_matrix(V_relative_elem, n1s, l1s; mask=mask1) + sp_V_matrix(V_relative_elem, n2s, l2s; mask=mask2)
@time "Moshinsky brackets" U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ) @time "Moshinsky brackets" U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
@time "V2" V2 = U' * V_relative * U @time "V2" V2 = U' * V_relative * U