From 6b20372abdb01fd8b631e0ac97f60626860fb046 Mon Sep 17 00:00:00 2001 From: Nuwan Yapa Date: Wed, 1 May 2024 17:49:38 -0400 Subject: [PATCH] Fixed error regarding construction of W matrix --- berggren_3body.jl | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/berggren_3body.jl b/berggren_3body.jl index 1d777ae..5ce9900 100644 --- a/berggren_3body.jl +++ b/berggren_3body.jl @@ -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)