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

Simplification

This commit is contained in:
Nuwan Yapa 2024-12-05 18:12:55 -05:00
parent bdf35583aa
commit e181f1cdec
3 changed files with 12 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
berggren_3body_resonance.jl
View File

@ -47,18 +47,13 @@ end
U = blockdiag(sparse.(U_blocks)...) U = blockdiag(sparse.(U_blocks)...)
end end
@time "T" begin @time "Block diagonal part" begin
T_blocks = [kron_sum(to_berg_basis(get_T_matrix(ks, μ), j1), to_berg_basis(get_T_matrix(ks, μ), j2)) for (j1, j2) in js] Hb_blocks = [kron_sum(to_berg_basis(get_H_matrix((k, kp) -> V_l(j1, k, kp), ks, ws, μ), j1), to_berg_basis(get_H_matrix((k, kp) -> V_l(j2, k, kp), ks, ws, μ), j2)) for (j1, j2) in js]
T = blockdiag(sparse.(T_blocks)...) Hb = blockdiag(sparse.(Hb_blocks)...)
end end
@time "T_cross" T_cross = get_2p_p1p2_matrix(length(ks), js, Λ, berg_bases, berg_bases, ws) ./ (2*μ) @time "T_cross" T_cross = get_2p_p1p2_matrix(length(ks), js, Λ, berg_bases, berg_bases, ws) ./ (2*μ)
@time "V1 and V2" begin
V_blocks = [kron_sum(to_berg_basis(get_V_matrix((k, kp) -> V_l(j1, k, kp), ks, ws), j1), to_berg_basis(get_V_matrix((k, kp) -> V_l(j2, k, kp), ks, ws), j2)) for (j1, j2) in js]
V = blockdiag(sparse.(V_blocks)...)
end
E_max = 30 E_max = 30
μω_global = 0.5 μω_global = 0.5
# due to simple relative coordinates # due to simple relative coordinates
@ -72,7 +67,7 @@ E_max = 30
@time "V12_p" V12_p = W * V12_HO * transpose(W) @time "V12_p" V12_p = W * V12_HO * transpose(W)
@time "V12" V12 = transpose(U) * V12_p * U @time "V12" V12 = transpose(U) * V12_p * U
@time "H" H = T + T_cross + V + V12 @time "H" H = Hb + T_cross + V12
@time "Eigenvalues" evals, _ = eigs(H, sigma=target, maxiter=5000, tol=1e-5, ritzvec=false, check=1) @time "Eigenvalues" evals, _ = eigs(H, sigma=target, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
display(evals) display(evals)

13
p_space_3body.jl
View File

@ -32,14 +32,9 @@ basis_size = length(js) * length(ks)^2
@assert length(basis) == basis_size "Something wrong with the basis" @assert length(basis) == basis_size "Something wrong with the basis"
println("Basis size = $basis_size") println("Basis size = $basis_size")
@time "T" begin @time "Block diagonal part" begin
T_blocks = [kron_sum(get_T_matrix(ks, μ1), get_T_matrix(ks, μ2)) for _ in js] Hb_blocks = [kron_sum(get_H_matrix((k, kp) -> V_l(j1, k, kp), ks, ws, μ1), get_T_matrix(ks, μ2)) for (j1,_ ) in js]
T = blockdiag(sparse.(T_blocks)...) Hb = blockdiag(sparse.(Hb_blocks)...)
end
@time "V1" begin
V1_blocks = [kron(get_V_matrix((k, kp) -> V_l(j1, k, kp), ks, ws), I(length(ks))) for (j1, _) in js]
V1 = blockdiag(sparse.(V1_blocks)...)
end end
E_max = 30 E_max = 30
@ -53,7 +48,7 @@ E_max = 30
@time "W_left" W_left = get_W_matrix(basis, E_max, Λ, μ1ω1, μ2ω2; weights=false) @time "W_left" W_left = get_W_matrix(basis, E_max, Λ, μ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 = T + V1 + V2 @time "H" H = Hb + 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)

13
p_space_3body_resonance.jl
View File

@ -31,14 +31,9 @@ basis_size = length(js) * length(ks)^2
@assert length(basis) == basis_size "Something wrong with the basis" @assert length(basis) == basis_size "Something wrong with the basis"
println("Basis size = $basis_size") println("Basis size = $basis_size")
@time "T" begin @time "Block diagonal part" begin
T_blocks = [kron_sum(get_T_matrix(ks, μ1), get_T_matrix(ks, μ2)) for _ in js] Hb_blocks = [kron_sum(get_H_matrix((k, kp) -> V_l(j1, k, kp), ks, ws, μ1), get_T_matrix(ks, μ2)) for (j1,_ ) in js]
T = blockdiag(sparse.(T_blocks)...) Hb = blockdiag(sparse.(Hb_blocks)...)
end
@time "V1" begin
V1_blocks = [kron(get_V_matrix((k, kp) -> V_l(j1, k, kp), ks, ws), I(length(ks))) for (j1, _) in js]
V1 = blockdiag(sparse.(V1_blocks)...)
end end
E_max = 30 E_max = 30
@ -52,7 +47,7 @@ E_max = 30
@time "W_left" W_left = get_W_matrix(basis, E_max, Λ, μ1ω1, μ2ω2; weights=false) @time "W_left" W_left = get_W_matrix(basis, E_max, Λ, μ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 = T + V1 + V2 @time "H" H = Hb + V2
@time "Eigenvalues" evals, _ = eigs(H, sigma=target, maxiter=5000, tol=1e-5, ritzvec=false, check=1) @time "Eigenvalues" evals, _ = eigs(H, sigma=target, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
display(evals) display(evals)