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

Proper mu and omega constants

This commit is contained in:
Nuwan Yapa 2024-04-17 16:25:06 -04:00
parent 1de8f45f84
commit 4c368ac903
5 changed files with 18 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
XZ_technique.jl
View File

@ -2,7 +2,7 @@ using LinearAlgebra, Plots
include("ho_basis.jl")
include("p_space.jl")
ω = 0.5 * exp(-1im * 0.47 * pi)
μω_gen = 0.5 * exp(-1im * 0.47 * pi)
μ = 0.5
l = 0
V1 = -5
@ -14,8 +14,8 @@ n_max = 15
ns = collect(0:n_max)
ls = fill(l, n_max + 1)
T = sp_T_matrix(ns, ls; ω=ω, μ=μ)
V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns); ω=ω) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns); ω=ω)
T = sp_T_matrix(ns, ls; μω_gen=μω_gen, μ=μ)
V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns); μω_gen=μω_gen) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns); μω_gen=μω_gen)
n_EC = 8
train_cs = (0.7 .+ 0.05 * randn(n_EC)) - 1im * (0.2 .+ 0.05 * randn(n_EC))

8
ho_basis.jl
View File

@ -7,8 +7,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.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)])
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[]
@ -54,7 +54,7 @@ function get_2p_basis(E_max)
return (Es, n1s, l1s, n2s, l2s)
end
function sp_T_matrix(ns, ls; mask=trues(length(ns),length(ns)), ω=1.0, μ=1.0)
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
@ -68,7 +68,7 @@ function sp_T_matrix(ns, ls; mask=trues(length(ns),length(ns)), ω=1.0, μ=1.0)
end
end
end
return (ω / μ) .* mat
return (μω_gen / μ) .* mat
end
function sp_V_matrix(V_l, ns, ls; mask=trues(length(ns),length(ns)), dtype=Float64)

6
ho_basis_2body.jl
View File

@ -3,7 +3,7 @@ include("ho_basis.jl")
include("p_space.jl")
E_max = 30
ω = 0.2
μω_gen = 0.2
Λ = 0
m = 1.0
@ -18,10 +18,10 @@ Es, n1s, l1s = get_sp_basis(E_max)
println("Basis size = ", length(Es))
println("Constructing KE matrices")
@time "T1" T1 = sp_T_matrix(n1s, l1s; ω=ω, μ=μ1)
@time "T1" T1 = sp_T_matrix(n1s, l1s; μω_gen=μω_gen, μ=μ1)
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; μω_gen=μω_gen)
@time "V1" V1 = sp_V_matrix(V1_elem, n1s, l1s)
println("Calculating spectrum")

10
ho_basis_3body.jl
View File

@ -3,7 +3,7 @@ include("ho_basis.jl")
include("p_space.jl")
E_max = 20
ω = 0.2
μω_gen = 0.2
Λ = 0
m = 1.0
@ -25,12 +25,12 @@ end
println("Basis size = ", length(Es))
println("Constructing KE matrices")
@time "T1" T1 = sp_T_matrix(n1s, l1s; mask=mask1, ω=ω, μ=μ1)
@time "T2" T2 = sp_T_matrix(n2s, l2s; mask=mask2, ω=ω, μ=μ2)
@time "T1" T1 = sp_T_matrix(n1s, l1s; mask=mask1, μω_gen=μω_gen, μ=μ1)
@time "T2" T2 = sp_T_matrix(n2s, l2s; mask=mask2, μω_gen=μω_gen, μ=μ2)
println("Constructing PE matrices")
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; ω=ω)
V1_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; μω_gen=μω_gen)
V_relative_elem(l, n1, n2) = Va * V_Gaussian(Ra / c, l, n1, n2; μω_gen=μω_gen)
@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; mask=mask1) + sp_V_matrix(V_relative_elem, n2s, l2s; mask=mask2)
@time "Moshinsky brackets" U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)

6
ho_basis_benchmark.jl
View File

@ -2,7 +2,7 @@ using LinearAlgebra, Plots
include("ho_basis.jl")
include("p_space.jl")
ω = 0.5 * exp(-1im * 0.45 * pi)
μω_gen = 0.5 * exp(-1im * 0.45 * pi)
μ = 0.5
l = 0
V1 = -5
@ -14,9 +14,9 @@ n_max = 15
ns = collect(0:n_max)
ls = fill(l, n_max + 1)
T = sp_T_matrix(ns, ls; ω=ω, μ=μ)
T = sp_T_matrix(ns, ls; μω_gen=μω_gen, μ=μ)
V_l(l, n1, n2) = V1 * V_Gaussian(R1, l, n1, n2; ω=ω) + V2 * V_Gaussian(R2, l, n1, n2; ω=ω)
V_l(l, n1, n2) = V1 * V_Gaussian(R1, l, n1, n2; μω_gen=μω_gen) + V2 * V_Gaussian(R2, l, n1, n2; μω_gen=μω_gen)
V = sp_V_matrix(V_l, ns, ls; dtype=ComplexF64)
cs = range(1.25, 0.25, 10)