85 lines
2.8 KiB
Julia
85 lines
2.8 KiB
Julia
using Plots
|
|
|
|
training_c = [1.1, 0.9, 0.7, 0.5]
|
|
extrapolating_c = 0.0 : 0.2 : 1.2
|
|
|
|
training_ref = reverse([1.4750633616275919 - 0.0003021770706749637im
|
|
1.9567078295375822 - 0.0007646829108872369im
|
|
2.4351117758403076 - 0.001281037843108658im
|
|
2.9096543462392357 - 0.002962488527470604im])
|
|
|
|
exact_ref = reverse([4.076662025307587-0.012709842443350328im,
|
|
3.613318119833891-0.007335804709990623im,
|
|
3.1453431847006783-0.004030580410326795im,
|
|
2.672967129943755-0.00211498327461944im,
|
|
2.196542557810288-0.0010719835443437104im,
|
|
1.7164583929199813-0.0005455212208182736im,
|
|
1.233088227541505-0.0003070320106485624im])
|
|
|
|
include("../p_space_3body_resonance.jl")
|
|
H0 = H
|
|
|
|
# Vp = perturbation to make the state artificially bound
|
|
Vp_of_r(r) = -exp(-(r/3)^2)
|
|
Vp_l(j, k, kp) = Vl_mat_elem(Vp_of_r, j, k, kp; atol=atol, maxevals=maxevals, R_cutoff=R_cutoff)
|
|
|
|
@time "Vp block diagonal part" begin
|
|
Vpb_blocks = [kron_sum(get_V_matrix((k, kp) -> Vp_l(j1, k, kp), ks, ws), spzeros(length(ks), length(ks))) for (j1, _) in js]
|
|
Vpb = blockdiag(sparse.(Vpb_blocks)...)
|
|
end
|
|
|
|
@time "Vp2_HO" Vp2_HO = get_jacobi_V2_matrix(Vp_of_r, basis_ho, μω_global; atol=atol, maxevals=maxevals)
|
|
@time "Vp2" Vp2 = W_left * Vp2_HO * transpose(W_right)
|
|
@time "Vp" Vp = Vpb + Vp2
|
|
|
|
# free memory
|
|
basis = Hb_blocks = Hb = basis_ho = V2_HO = W_right = W_left = V2 = nothing
|
|
GC.gc()
|
|
|
|
exact = ComplexF64[]
|
|
training = ComplexF64[]
|
|
extrapolated = ComplexF64[]
|
|
training_vecs = Vector{ComplexF64}[]
|
|
|
|
for c in training_c
|
|
println("Training for c = $c")
|
|
global current_E = pop!(training_ref)
|
|
|
|
local H = H0 + c .* Vp
|
|
local evals, evecs = eigs(H, sigma=current_E, maxiter=5000, tol=1e-5, ritzvec=true, check=1)
|
|
|
|
global current_E = nearest(evals, current_E)
|
|
push!(training, current_E)
|
|
push!(training_vecs, evecs[:, nearestIndex(evals, current_E)])
|
|
end
|
|
|
|
weights = repeat(kron(ws, ws), jmax + 1)
|
|
weights_mat = spdiagm(weights)
|
|
|
|
println("Original EC dimensionality = $(length(training_vecs))")
|
|
@time "Gram-Schmidt" training_vecs = gram_schmidt!(training_vecs, weights; verbose=true) # orthonormalization
|
|
|
|
EC_basis = hcat(training_vecs...)
|
|
H0_EC = transpose(EC_basis) * weights_mat * H0 * EC_basis
|
|
Vp_EC = transpose(EC_basis) * weights_mat * Vp * EC_basis
|
|
|
|
for c in extrapolating_c
|
|
println("Extrapolating for c = $c")
|
|
global current_E = pop!(exact_ref)
|
|
|
|
local H = H0 + c .* Vp
|
|
local evals, _ = eigs(H, sigma=current_E, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
|
|
|
|
global current_E = nearest(evals, current_E)
|
|
push!(exact, current_E)
|
|
|
|
# extrapolation
|
|
H_EC = H0_EC + c .* Vp_EC
|
|
evals = eigvals(H_EC)
|
|
push!(extrapolated, nearest(evals, current_E))
|
|
end
|
|
|
|
scatter(real.(training),imag.(training), label="training")
|
|
scatter!(real.(exact),imag.(exact), label="exact")
|
|
scatter!(real.(extrapolated),imag.(extrapolated), label="extrapolated")
|
|
savefig("temp/Berggren_R2R.pdf") |