112 lines
3.1 KiB
Julia
112 lines
3.1 KiB
Julia
using SparseArrays
|
|
using SpecialFunctions
|
|
include("helper.jl")
|
|
include("osbrackets.jl")
|
|
|
|
# Gaussian potentials in HO space
|
|
inv_factorial(n) = Iterators.prod(inv.(1:n))
|
|
sqrt_factorial(n) = Iterators.prod(sqrt.(n:-1:1))
|
|
N_lnk(l, n, k) = 1/sqrt_factorial(n) * binomial(n, k) * sqrt(gamma(n + l + 3/2)) / gamma(k + l + 3/2)
|
|
Talmi(l, R, k1, k2; ω=1.0) = (-1)^(k1 + k2) * (1 + 1/(ω * R^2))^-(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2)
|
|
V_Gaussian(R, l, n1, n2; ω=1.0) = (-1)^(n1 + n2) * better_sum([N_lnk(l, n1, k1) * N_lnk(l, n2, k2) * Talmi(l, R, k1, k2; ω=ω) for (k1, k2) in Iterators.product(0:n1, 0:n2)])
|
|
|
|
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)
|
|
Es = Int[]
|
|
n1s = Int[]
|
|
l1s = Int[]
|
|
n2s = Int[]
|
|
l2s = Int[]
|
|
|
|
# E = 2*n1 + l1 + 2*n2 + l2
|
|
for E in 2*E_max : -2 : 0
|
|
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
|
|
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 sp_T_matrix(ns, ls; mask=trues(length(ns),length(ns)), ω=1.0, μ=1.0)
|
|
mat = spzeros(length(ns), length(ns))
|
|
for idx in CartesianIndices(mat)
|
|
if !mask[idx]; continue; end
|
|
(i, j) = Tuple(idx)
|
|
if ls[i] == ls[j]
|
|
if ns[i] == ns[j]
|
|
mat[idx] = ns[j] + ls[i]/2 + 3/4
|
|
elseif abs(ns[i]-ns[j]) == 1
|
|
n_max = max(ns[i], ns[j])
|
|
mat[idx] = -(1/2) * sqrt(n_max * (n_max + ls[i] + 1/2))
|
|
end
|
|
end
|
|
end
|
|
return (ω / μ) .* mat
|
|
end
|
|
|
|
function sp_V_matrix(V_l, ns, ls; mask=trues(length(ns),length(ns)), dtype=Float64)
|
|
mat = zeros(dtype, length(ns), length(ns))
|
|
Threads.@threads for idx in CartesianIndices(mat)
|
|
if !mask[idx]; continue; end
|
|
(i, j) = Tuple(idx)
|
|
if ls[i] == ls[j]
|
|
mat[idx] = V_l(ls[i], ns[i], ns[j])
|
|
end
|
|
end
|
|
return sparse(mat)
|
|
end
|
|
|
|
function Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
|
|
Emax = maximum(Es) ÷ 2 # TODO: Too many steps. Simplify.
|
|
|
|
ul1s = unique(l1s)
|
|
ul2s = unique(l2s)
|
|
BRACs = Matrix{Array}(undef, length(ul1s), length(ul2s))
|
|
for l1 in ul1s
|
|
for l2 in ul2s
|
|
BRACs[l1, l2] = cal_BRAC(Emax, Λ, l1, l2)
|
|
end
|
|
end
|
|
|
|
mat = spzeros(length(Es), length(Es))
|
|
s = hcat(Es, n1s, l1s, n2s, l2s)
|
|
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
|
|
phase = (-1)^(N + n + n1 + n2)
|
|
mat[i, j] = phase * get_bracket(BRACs[l1, l2], Emax, Λ, N, L, n, l, n1, n2)
|
|
end
|
|
end
|
|
|
|
return mat
|
|
end
|