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

Introduced omega

This commit is contained in:
Nuwan Yapa 2024-02-27 13:45:43 -05:00
parent 14b574500c
commit b4e7897cb7
2 changed files with 7 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
ho_basis.jl
View File

@ -6,8 +6,8 @@ include("helper.jl")
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)^(k1 + k2) * (1 + 1/R^2)^-(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2)
V_Gaussian(R, l, n1, n2) = (-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)])
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[]
@ -55,7 +55,7 @@ end
get_V_matrix(V_l, ls, ns) = throw("unimplemented")
function sp_T_matrix(ns, ls)
function sp_T_matrix(ns, ls; ω=1.0)
mat = zeros(length(ns), length(ns))
for idx in CartesianIndices(mat)
(i, j) = Tuple(idx)
@ -70,7 +70,7 @@ function sp_T_matrix(ns, ls)
end
end
end
return mat
return ω .* mat
end
get_H_matrix(V_l, ns, ls) = get_T_matrix(ns, ls) + get_V_matrix(V_l, ns, ls)

5
ho_basis_benchmark.jl
View File

@ -2,6 +2,7 @@ using LinearAlgebra, Plots
include("ho_basis.jl")
include("p_space.jl")
ω = 0.5
μ = 0.5
l = 0
V1 = -5
@ -13,8 +14,8 @@ n_max = 15
ns = collect(0:n_max)
ls = fill(l, n_max + 1)
T = 1/μ * sp_T_matrix(ns, ls)
V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns)) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns))
T = 1/μ * sp_T_matrix(ns, ls; ω=ω)
V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns); ω=ω) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns); ω=ω)
cs = range(1.35, 0.9, 5)