99 lines
3.2 KiB
Julia
99 lines
3.2 KiB
Julia
using LinearAlgebra, SparseArrays, Arpack
|
|
include("helper.jl")
|
|
include("p_space.jl")
|
|
include("ho_basis.jl")
|
|
|
|
println("No of threads = ", Threads.nthreads())
|
|
|
|
atol = 10^-5
|
|
maxevals = 10^5
|
|
|
|
# E_target = -0.3919
|
|
|
|
m = 1.0
|
|
μ1 = m * 1/2
|
|
μ2 = m * 2/3
|
|
|
|
Va = -2
|
|
Ra = 2
|
|
V_of_r(r) = Va * exp(-r^2 / Ra^2)
|
|
V_l(j, k, kp) = Vl_mat_elem(V_of_r, j, k, kp; atol=atol, maxevals=maxevals, R_cutoff=16)
|
|
|
|
vertices = [0, 0.5 - 0.3im, 1, 4]
|
|
subdivisions = [10, 10, 10]
|
|
ks, ws = get_mesh(vertices, subdivisions)
|
|
js = collect(0:4)
|
|
|
|
sp_basis_size = length(ks) * length(js)
|
|
sp_ws = spdiagm(vcat(fill(ws, length(js))...))
|
|
|
|
basis = collect(Iterators.product(zip(ks, ws), js, zip(ks, ws), js))[:]
|
|
for ((k1, w1), j1, (k2, w2), j2) in basis
|
|
@assert k1 ∈ ks && k2 ∈ ks && w1 ∈ ws && w2 ∈ ws && j1 ∈ js && j2 ∈ js "Something wrong with the basis ordering"
|
|
end
|
|
println("Basis size = $(length(basis))")
|
|
|
|
@time "T1" begin
|
|
sp_T1_j = [get_T_matrix(ks, μ1) for j in js]
|
|
sp_T1 = blockdiag(sparse.(sp_T1_j)...)
|
|
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
|
|
|
|
@time "V1" begin
|
|
sp_V1_j = [get_V_matrix((k, kp) -> V_l(j, k, kp), ks, ws) for j in js]
|
|
sp_V1 = blockdiag(sparse.(sp_V1_j)...)
|
|
V1 = kron(sp_V1, I(sp_basis_size))
|
|
end
|
|
|
|
#############################################################################################
|
|
|
|
function get_W_matrix(basis, basis_HO, μ1ω1, μ2ω2=μ1ω1)
|
|
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
|
|
end
|
|
end
|
|
return sparse(W)
|
|
end
|
|
|
|
############################################################################################
|
|
|
|
Λ = 0
|
|
E_max = 24
|
|
μω_global = 0.5
|
|
|
|
μ1ω1 = μω_global * 1/2
|
|
μ2ω2 = μω_global * 2
|
|
|
|
basis_HO = get_2p_basis(E_max)
|
|
Es, n1s, l1s, n2s, l2s = basis_HO
|
|
l_max = max(maximum(l1s), maximum(l2s))
|
|
n_max = max(maximum(n1s), maximum(n2s))
|
|
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
|
|
mask2 = (n1s .== n1s') .&& (l1s .== l1s')
|
|
|
|
V_relative_elem(l, n1, n2) = V_numerical(V_of_r, l, n1, n2; μω_gen=μω_global, atol=atol, maxevals=maxevals)
|
|
V_relative_cache = LRU{Tuple{UInt8, UInt8, UInt8}, Float64}(maxsize=(1+l_max)*(1+n_max)^2)
|
|
@time "V relative" V_relative = sp_V_matrix(V_relative_elem, n1s, l1s; mask=mask1, dtype=Float64, cache=V_relative_cache) + sp_V_matrix(V_relative_elem, n2s, l2s; mask=mask2, dtype=Float64, cache=V_relative_cache)
|
|
@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 "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)
|
|
display(evals)
|