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

EC for 3 body resonance

This commit is contained in:
Nuwan Yapa 2024-05-04 21:40:43 -04:00
parent 6b20372abd
commit 3200b7e041
2 changed files with 57 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

53
3body_resonance_EC.jl Normal file
View File

@ -0,0 +1,53 @@
using Plots
include("3body_resonance.jl")
current_E = 5.9673 - 0.0006im
training_c = 2.0 : -0.2 : 1.2
extrapolating_c = 1.05 .- [0.0 : 0.1 : 0.4; 0.45 : 0.05 : 0.60]
@time "H0" H0 = T1 + T2
@time "V" V = V1 + V2
# free memory
Es = n1s = l1s = n2s = l2s = mask1 = mask2 = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
GC.gc()
exact = ComplexF64[]
training = ComplexF64[]
extrapolated = ComplexF64[]
training_vecs = Vector{ComplexF64}[]
for c in training_c
println("Training for c = $c")
H = H0 + c .* V
evals, evecs = eigs(H, nev=3, ncv=24, which=:LI, 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
EC_basis = hcat(training_vecs...)
N_EC = transpose(EC_basis) * EC_basis
H0_EC = transpose(EC_basis) * H0 * EC_basis
V_EC = transpose(EC_basis) * V * EC_basis
for c in extrapolating_c
println("Extrapolating for c = $c")
H = H0 + c .* V
evals, evecs = eigs(H, nev=3, ncv=24, which=:LI, maxiter=5000, tol=1e-5, ritzvec=true, check=1)
global current_E = nearest(evals, current_E)
push!(exact, current_E)
# extrapolation
H_EC = H0_EC + c .* V_EC
evals = eigvals(H_EC, N_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/NCSM.pdf")

5
helper.jl
View File

@ -16,8 +16,11 @@ better_sum(arr::Array{ComplexF64}) = better_sum(real.(arr)) + 1im * better_sum(i
"The triangle inequality for angular momenta"
triangle_ineq(l1, l2, L) = abs(l1 - l2) ≤ L && L ≤ (l1 + l2)
"Index of the nearest value in a list to a given reference point"
nearestIndex(list::Array, ref) = argmin(norm.(list .- ref))
"Nearest value in a list to a given reference point"
nearest(list::Array, ref) = list[argmin(norm.(list .- ref))]
nearest(list::Array, ref) = list[nearestIndex(list, ref)]
"Simple implementation of the Kronecker sum"
function kron_sum(A::AbstractMatrix, B::AbstractMatrix)