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

Reorganize

This commit is contained in:
Nuwan Yapa 2024-09-19 18:16:37 -04:00
parent c61adf6f2b
commit 6e5e1d0132
1 changed files with 14 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
ho_basis.jl
View File

@ -11,6 +11,20 @@ N_lnk(l, n, k) = 1/sqrt_factorial(n) * binomial(n, k) * sqrt(gamma(n + 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) 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)]) 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)])
# 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
##############################################################
function get_sp_basis(E_max) function get_sp_basis(E_max)
Es = Int[] Es = Int[]
ns = Int[] ns = Int[]
@ -126,18 +140,6 @@ 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] return BRAC[1 + NP, 1 + n1′, 1 + MP, 1 + n1, 1 + n2, 1 + N, 1 + M, 1]
end 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
function get_jacobi_V_matrix(V_of_r, E_max, Λ, μ1ω1, μω_global; atol=10^-6, maxevals=10^5) function get_jacobi_V_matrix(V_of_r, E_max, Λ, μ1ω1, μω_global; atol=10^-6, maxevals=10^5)
_, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ) _, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
l_max = max(maximum(l1s), maximum(l2s)) l_max = max(maximum(l1s), maximum(l2s))