32 lines
1.1 KiB
Julia
32 lines
1.1 KiB
Julia
using Arpack
|
|
include("ho_basis.jl")
|
|
|
|
coord_system = src
|
|
atol = 10 ^ -5
|
|
maxevals = 10 ^ 5
|
|
|
|
Λ = 0
|
|
m = 1.0
|
|
V_of_r(r) = -2 * exp(-r^2 / 4)
|
|
|
|
# three-body potential V3b = V3b_of_r(r12) * V3b_of_r(r23) * V3b_of_r(r31)
|
|
V3b_of_r(r) = -1 * exp(-r^2 / 4)
|
|
|
|
E_max = 40
|
|
μω_global = 0.3
|
|
|
|
# compute three-body potential matrix elements
|
|
μω = μω_global * 2
|
|
basis = ho_basis_2B(E_max, Λ)
|
|
V3b_elem(l, n1, n2) = V_numerical(V3b_of_r, l, n1, n2; μω_gen=μω, atol=atol, maxevals=maxevals)
|
|
V_cache = prealloc_V_cache(basis.E_max, ComplexF64)
|
|
@time "V1" V1 = get_sp_V_matrix(V3b_elem, basis.n1s, basis.l1s, [basis.n2s, basis.l2s]; dtype=ComplexF64, cache=V_cache)
|
|
@time "V2" V2 = get_sp_V_matrix(V3b_elem, basis.n2s, basis.l2s, [basis.n1s, basis.l1s]; dtype=ComplexF64, cache=V_cache)
|
|
@time "V12" V12 = get_src_V12_matrix(V3b_of_r, basis, μω_global; atol=atol, maxevals=maxevals)
|
|
@time "V3b" V3b = (V1 * V2) * V12 # commutative?
|
|
|
|
H = get_3b_H_matrix(coord_system, V_of_r, μω_global, E_max, Λ, m) + V3b
|
|
|
|
@time "Eigenvalues" evals, _ = eigs(H, nev=3, ncv=30, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
|
|
display(evals)
|