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

Base constructed according to triangle inequality

This commit is contained in:
ysyapa 2024-06-04 19:44:39 +00:00
parent 59f74b0263
commit d230476589
1 changed files with 16 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
berggren_3body.jl
View File

@ -8,6 +8,7 @@ println("No of threads = ", Threads.nthreads())
atol = 10^-5 atol = 10^-5
maxevals = 10^5 maxevals = 10^5
Λ = 0
m = 1.0 m = 1.0
μ1 = m * 1/2 μ1 = m * 1/2
μ2 = m * 2/3 μ2 = m * 2/3
@ -20,33 +21,24 @@ V_l(j, k, kp) = Vl_mat_elem(V_of_r, j, k, kp; atol=atol, maxevals=maxevals, R_cu
vertices = [0, 0.5 - 0.3im, 1, 4] vertices = [0, 0.5 - 0.3im, 1, 4]
subdivisions = [10, 10, 10] subdivisions = [10, 10, 10]
ks, ws = get_mesh(vertices, subdivisions) ks, ws = get_mesh(vertices, subdivisions)
js = collect(0:4)
sp_basis_size = length(ks) * length(js) jmax = 4
sp_ws = spdiagm(vcat(fill(ws, length(js))...)) tri((j1, j2)) = triangle_ineq(j1, j2, Λ)
js = collect(Iterators.filter(tri, Iterators.product(0:jmax, 0:jmax)))
basis = collect(Iterators.product(zip(ks, ws), js, zip(ks, ws), js))[:] basis = collect(Iterators.product(zip(ks, ws), zip(ks, ws), js))[:]
for ((k1, w1), j1, (k2, w2), j2) in basis basis_size = length(ks)^2 * length(js)
@assert k1 ∈ ks && k2 ∈ ks && w1 ∈ ws && w2 ∈ ws && j1 ∈ js && j2 ∈ js "Something wrong with the basis ordering" @assert length(basis) == basis_size "Something wrong with the basis"
end println("Basis size = $basis_size")
println("Basis size = $(length(basis))")
@time "T1" begin @time "T" begin
sp_T1_j = [get_T_matrix(ks, μ1) for j in js] T_blocks = [kron_sum(get_T_matrix(ks, μ1), get_T_matrix(ks, μ2)) for _ in js]
sp_T1 = blockdiag(sparse.(sp_T1_j)...) T = blockdiag(sparse.(T_blocks)...)
T1 = kron(sp_T1, I(sp_basis_size))
end
@time "T2" begin
sp_T2_j = [get_T_matrix(ks, μ2) for j in js]
sp_T2 = blockdiag(sparse.(sp_T2_j)...)
T2 = kron(I(sp_basis_size), sp_T2)
end end
@time "V1" begin @time "V1" begin
sp_V1_j = [get_V_matrix((k, kp) -> V_l(j, k, kp), ks, ws) for j in js] V1_blocks = [kron(get_V_matrix((k, kp) -> V_l(j1, k, kp), ks, ws), I(length(ks))) for (j1, _) in js]
sp_V1 = blockdiag(sparse.(sp_V1_j)...) V1 = blockdiag(sparse.(V1_blocks)...)
V1 = kron(sp_V1, I(sp_basis_size))
end end
############################################################################################# #############################################################################################
@ -56,7 +48,7 @@ function get_W_matrix(basis, basis_HO, μ1ω1, μ2ω2=μ1ω1; weights=true)
W = zeros(ComplexF64, length(basis), length(Es)) W = zeros(ComplexF64, length(basis), length(Es))
Threads.@threads for idx in CartesianIndices(W) Threads.@threads for idx in CartesianIndices(W)
(i1, i2) = Tuple(idx) (i1, i2) = Tuple(idx)
((k1, w1), j1, (k2, w2), j2) = basis[i1] ((k1, w1), (k2, w2), (j1, j2)) = basis[i1]
if j1 == l1s[i2] && j2 == l2s[i2] if j1 == l1s[i2] && j2 == l2s[i2]
elem1 = 1/sqrt(sqrt(μ1ω1)) * (-1)^n1s[i2] * ho_basis(j1, n1s[i2], 1/sqrt(μ1ω1) * k1) elem1 = 1/sqrt(sqrt(μ1ω1)) * (-1)^n1s[i2] * ho_basis(j1, n1s[i2], 1/sqrt(μ1ω1) * k1)
elem2 = 1/sqrt(sqrt(μ2ω2)) * (-1)^n2s[i2] * ho_basis(j2, n2s[i2], 1/sqrt(μ2ω2) * k2) elem2 = 1/sqrt(sqrt(μ2ω2)) * (-1)^n2s[i2] * ho_basis(j2, n2s[i2], 1/sqrt(μ2ω2) * k2)
@ -68,14 +60,13 @@ end
############################################################################################ ############################################################################################
Λ = 0
E_max = 30 E_max = 30
μω_global = 0.5 μω_global = 0.5
μ1ω1 = μω_global * 1/2 μ1ω1 = μω_global * 1/2
μ2ω2 = μω_global * 2 μ2ω2 = μω_global * 2
basis_HO = get_2p_basis(E_max) basis_HO = get_2p_basis(E_max, Λ)
Es, n1s, l1s, n2s, l2s = basis_HO Es, n1s, l1s, n2s, l2s = basis_HO
l_max = max(maximum(l1s), maximum(l2s)) l_max = max(maximum(l1s), maximum(l2s))
n_max = max(maximum(n1s), maximum(n2s)) n_max = max(maximum(n1s), maximum(n2s))
@ -92,7 +83,7 @@ V_relative_cache = LRU{Tuple{UInt8, UInt8, UInt8}, Float64}(maxsize=(1+l_max)*(1
@time "W_left" W_left = get_W_matrix(basis, basis_HO, μ1ω1, μ2ω2; weights=false) @time "W_left" W_left = get_W_matrix(basis, basis_HO, μ1ω1, μ2ω2; weights=false)
@time "V2" V2 = W_left * V2_HO * transpose(W_right) @time "V2" V2 = W_left * V2_HO * transpose(W_right)
@time "H" H = T1 + V1 + T2 + V2 @time "H" H = T + V1 + V2
@time "Eigenvalues" evals, _ = eigs(H, nev=3, ncv=24, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1) @time "Eigenvalues" evals, _ = eigs(H, nev=3, ncv=24, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
display(evals) display(evals)