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

Fixed error regarding construction of W matrix

This commit is contained in:
Nuwan Yapa 2024-05-01 17:49:38 -04:00
parent d90551577d
commit 6b20372abd
1 changed files with 9 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
berggren_3body.jl
View File

@ -53,16 +53,16 @@ end
#############################################################################################
function get_W_matrix(basis, basis_HO, μ1ω1, μ2ω2=μ1ω1)
function get_W_matrix(basis, basis_HO, μ1ω1, μ2ω2=μ1ω1; weights=true)
Es, n1s, l1s, n2s, l2s = basis_HO
W = zeros(ComplexF64, length(basis), length(Es))
Threads.@threads for idx in CartesianIndices(W)
(i1, i2) = Tuple(idx)
((k1, w1), j1, (k2, w2), j2) = basis[i1]
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) * w1
elem2 = 1/sqrt(sqrt(μ2ω2)) * (-1)^n2s[i2] * ho_basis(j2, n2s[i2], 1/sqrt(μ2ω2) * k2) * w2
W[idx] = elem1 * elem2
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)
W[idx] = elem1 * elem2 * (weights ? w1 * w2 : 1)
end
end
return sparse(W)
@ -71,7 +71,7 @@ end
############################################################################################
Λ = 0
E_max = 24
E_max = 30
μω_global = 0.5
μ1ω1 = μω_global * 1/2
@ -90,8 +90,10 @@ V_relative_cache = LRU{Tuple{UInt8, UInt8, UInt8}, Float64}(maxsize=(1+l_max)*(1
@time "Moshinsky brackets" U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
@time "V2_HO" V2_HO = U' * V_relative * U
@time "W" W = get_W_matrix(basis, basis_HO, μ1ω1, μ2ω2)
@time "V2" V2 = W * V2_HO * transpose(W)
@time "W_right" W_right = get_W_matrix(basis, basis_HO, μ1ω1, μ2ω2; weights=true)
@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 "H" H = T1 + V1 + T2 + V2
@time "Eigenvalues" evals, _ = eigs(H, nev=3, ncv=24, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)