142 lines
5.1 KiB
Julia
142 lines
5.1 KiB
Julia
using SparseArrays
|
||
using SpecialFunctions
|
||
using QuadGK
|
||
include("helper.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; μω_gen=1.0) = (-1)^(k1 + k2) * (1 + 1/(μω_gen * R^2))^-(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2)
|
||
V_Gaussian(R, l, n1, n2; μω_gen=1.0) = (-1)^(n1 + n2) * better_sum([N_lnk(l, n1, k1) * N_lnk(l, n2, k2) * Talmi(l, R, k1, k2; μω_gen=μω_gen) 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 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
|
||
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)), μω_gen=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 (μω_gen / μ) .* mat
|
||
end
|
||
|
||
function sp_V_matrix(V_l, ns, ls; mask=trues(length(ns),length(ns)), dtype=Float64, cache=fill(convert(dtype, NaN), 1+maximum(ls), 1+maximum(ns), 1+maximum(ns)))
|
||
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]
|
||
l = ls[i]
|
||
n1, n2 = minmax(ns[i], ns[j]) # assuming transpose symmetry
|
||
if isnan(cache[1+l, 1+n1, 1+n2])
|
||
cache[1+l, 1+n1, 1+n2] = V_l(l, n1, n2) # hopefully no race condition
|
||
@assert !isnan(cache[1+l, 1+n1, 1+n2]) "V matrix element returned NaN"
|
||
end
|
||
mat[idx] = cache[1+l, 1+n1, 1+n2]
|
||
end
|
||
end
|
||
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"
|
||
|
||
LMIN = Λ
|
||
LMAX = Λ
|
||
CO = 1/sqrt(2)
|
||
SI = 1/sqrt(2)
|
||
|
||
# dimensions BRAC(0:LMAX,0:(NQMAX-LMIN)/2,0:(NQMAX-LMIN)/2,0:(NQMAX-LMIN)/2,0:(NQMAX-LMIN)/2,0:LMAX,0:(NQMAX-LMIN)/2,LMIN:LMAX)
|
||
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))
|
||
|
||
s = hcat(Es, n1s, l1s, n2s, 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, Λ)
|
||
end
|
||
end
|
||
|
||
return sparse(mat)
|
||
end
|
||
|
||
function pick_Moshinsky_bracket(BRAC, n1′, l1′, n2′, l2′, n1, l1, n2, l2, Λ) # Efros notation -- don't confuse
|
||
ϵ = (l1 + l2 - Λ) % 2
|
||
NP = (l1′ - l2′ + Λ - ϵ) ÷ 2
|
||
MP = (l1′ + l2′ - Λ - ϵ) ÷ 2
|
||
N = (l1 - l2 + Λ - ϵ) ÷ 2
|
||
M = (l1 + l2 - Λ - ϵ) ÷ 2
|
||
|
||
# BRAC(NP,N1P,MP,N1,N2,N,M,L)
|
||
return BRAC[1 + NP, 1 + n1′, 1 + MP, 1 + n1, 1 + n2, 1 + N, 1 + M, 1]
|
||
end
|
||
|
||
# numerical evaluation of V matrix elements
|
||
sqrt_double_factorial(n) = Iterators.prod(sqrt.(n:-2:1))
|
||
sqrt_sqrt_pi = sqrt(sqrt(pi))
|
||
laguerre(l, n, x) = gamma(n + l + 3/2) * better_sum([(-x * x)^k / gamma(k + l + 3/2) * inv_factorial(n - k) * inv_factorial(k) for k in 0:n])
|
||
ho_basis(l, n, x) = (-1)^n / sqrt_sqrt_pi * 2^((n + l + 2) / 2) * sqrt_factorial(n) / sqrt_double_factorial(2*n + 2*l + 1) * x^(l + 1) * exp(-x^2 / 2) * laguerre(l, n, x)
|
||
|
||
function V_numerical(V_of_r, l, n1, n2; μω_gen=1.0, atol=0, maxevals=10^7)
|
||
integrand(r) = sqrt(μω_gen) * ho_basis(l, n1, sqrt(μω_gen) * r) * ho_basis(l, n2, sqrt(μω_gen) * r) * V_of_r(r)
|
||
(integral, _) = quadgk(integrand, 0, Inf; atol=atol, maxevals=maxevals)
|
||
return integral
|
||
end
|