EC.jl implemented for all 3-body systems
This commit is contained in:
parent
9084820ddc
commit
29bbceac03
15
EC.jl
15
EC.jl
|
|
@ -23,18 +23,18 @@ end
|
|||
If a list is provided for ref_eval, they are used as reference values for picking the closest eigenvalues at each sampling point.
|
||||
If a single number is provided for ref_eval, it is used as a reference for the first point, and the previous eigenvalue is used as the reference for each successive point.
|
||||
Set orthonormalize_threshold=0 to skip Gram-Schmidt orthonormalization and use GEVP. Otherwise this value is used as the threshold for dropping redundant vectors."
|
||||
function train!(EC::affine_EC, c_vals; ref_eval=nothing, orthonormalize_threshold=1e-5, CAEC=false, verbose=true, tol=1e-5)
|
||||
function train!(EC::affine_EC, c_vals; ref_eval=-10.0, orthonormalize_threshold=1e-5, CAEC=false, verbose=true, tol=1e-5)
|
||||
training_vecs = Vector{ComplexF64}[]
|
||||
current_E = -10.0 # randomly chosen reference value
|
||||
|
||||
for c in c_vals
|
||||
verbose && println("Training for c = $c")
|
||||
|
||||
global current_E
|
||||
if ref_eval isa Number
|
||||
current_E = ref_eval
|
||||
ref_eval = nothing
|
||||
elseif !isnothing(ref_eval)
|
||||
current_E = pop!(ref_eval)
|
||||
current_E = popfirst!(ref_eval)
|
||||
end
|
||||
|
||||
H = EC.H0 + c .* EC.H1
|
||||
|
|
@ -64,22 +64,21 @@ function train!(EC::affine_EC, c_vals; ref_eval=nothing, orthonormalize_threshol
|
|||
EC.trained = true
|
||||
end
|
||||
|
||||
"Extrapolate using a train EC model for a given range of c values
|
||||
"Extrapolate using a trained EC model for a given range of c values
|
||||
If a list is provided for ref_eval, they are used as reference values for picking the closest eigenvalues at each point.
|
||||
If a single number is provided for ref_eval, it is used as a reference for the first point, and the previous eigenvalue is used as the reference for each successive point."
|
||||
function extrapolate!(EC::affine_EC, c_vals; ref_eval=nothing, verbose=true, tol=1e-5)
|
||||
function extrapolate!(EC::affine_EC, c_vals; ref_eval=EC.training_E[end], verbose=true, tol=1e-5)
|
||||
@assert EC.trained "EC model must be trained using train() before extrapolation"
|
||||
|
||||
current_E = -10.0 # randomly chosen reference value
|
||||
|
||||
for c in c_vals
|
||||
verbose && println("Extact solution for c = $c")
|
||||
|
||||
global current_E
|
||||
if ref_eval isa Number
|
||||
current_E = ref_eval
|
||||
ref_eval = nothing
|
||||
elseif !isnothing(ref_eval)
|
||||
current_E = pop!(ref_eval)
|
||||
current_E = popfirst!(ref_eval)
|
||||
end
|
||||
|
||||
H = EC.H0 + c .* EC.H1
|
||||
|
|
|
|||
|
|
@ -1,20 +1,20 @@
|
|||
using Plots
|
||||
include("../EC.jl")
|
||||
|
||||
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
|
||||
training_ref = [1.4750633616275919 - 0.0003021770706749637im
|
||||
1.9567078295375822 - 0.0007646829108872369im
|
||||
2.4351117758403076 - 0.001281037843108658im
|
||||
2.9096543462392357 - 0.002962488527470604im])
|
||||
2.9096543462392357 - 0.002962488527470604im]
|
||||
|
||||
exact_ref = reverse([4.076662025307587-0.012709842443350328im,
|
||||
extrapolating_ref = [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])
|
||||
1.233088227541505-0.0003070320106485624im]
|
||||
|
||||
include("../p_space_3body_resonance.jl")
|
||||
H0 = H
|
||||
|
|
@ -32,54 +32,11 @@ end
|
|||
@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 = affine_EC(H0, Vp, weights)
|
||||
train!(EC, training_c; ref_eval=training_ref, CAEC=false)
|
||||
extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref)
|
||||
|
||||
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")
|
||||
exportCSV(EC, "temp/Berggren_R2R.csv")
|
||||
plot(EC, "temp/Berggren_R2R.pdf")
|
||||
|
|
@ -1,7 +1,7 @@
|
|||
include("../EC.jl")
|
||||
|
||||
training_ref = -0.72763
|
||||
exact_ref = 4.0766890719636635 - 0.01275892774109674im
|
||||
extrapolating_ref = 4.0766890719636635 - 0.01275892774109674im
|
||||
|
||||
training_c = [2.0, 1.9, 1.8]
|
||||
extrapolating_c = 0.0 : 0.2 : 1.2
|
||||
|
|
@ -15,7 +15,7 @@ Vp_of_r(r) = -exp(-(r/3)^2)
|
|||
|
||||
EC = affine_EC(H0, Vp)
|
||||
train!(EC, training_c; ref_eval=training_ref, CAEC=true)
|
||||
extrapolate!(EC, extrapolating_c; ref_eval=exact_ref)
|
||||
extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref)
|
||||
|
||||
exportCSV(EC, "temp/HO_B2R.csv")
|
||||
plot(EC, "temp/HO_B2R.pdf")
|
||||
|
|
|
|||
|
|
@ -1,56 +1,16 @@
|
|||
using DelimitedFiles, Plots
|
||||
include("../ho_basis_3body_resonance.jl")
|
||||
include("../EC.jl")
|
||||
|
||||
current_E = 5.9673 - 0.0006im
|
||||
ref_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
|
||||
include("../ho_basis_3body_resonance.jl")
|
||||
@time "H0" H0 = T1 + T2 + T_cross
|
||||
|
||||
# free memory
|
||||
basis = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
|
||||
GC.gc()
|
||||
EC = affine_EC(H0, V)
|
||||
train!(EC, training_c; ref_eval=ref_E, CAEC=false)
|
||||
extrapolate!(EC, extrapolating_c)
|
||||
|
||||
exact = ComplexF64[]
|
||||
training = ComplexF64[]
|
||||
extrapolated = ComplexF64[]
|
||||
training_vecs = Vector{ComplexF64}[]
|
||||
|
||||
for c in training_c
|
||||
println("Training for c = $c")
|
||||
local H = H0 + c .* V
|
||||
local 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
|
||||
|
||||
println("Original EC dimensionality = $(length(training_vecs))")
|
||||
@time "Gram-Schmidt" training_vecs = gram_schmidt!(training_vecs; verbose=true) # orthonormalization
|
||||
|
||||
EC_basis = hcat(training_vecs...)
|
||||
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")
|
||||
local H = H0 + c .* V
|
||||
local 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)
|
||||
push!(extrapolated, nearest(evals, current_E))
|
||||
end
|
||||
|
||||
exportCSV("temp/NCSM.csv", (training, exact, extrapolated), ("training", "exact", "extrapolated"))
|
||||
|
||||
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")
|
||||
exportCSV(EC, "temp/HO_R2R.csv")
|
||||
plot(EC, "temp/HO_R2R.pdf")
|
||||
|
|
|
|||
|
|
@ -1,6 +1,5 @@
|
|||
using Arpack, SparseArrays, LRUCache
|
||||
using DelimitedFiles, Plots
|
||||
include("../ho_basis.jl")
|
||||
include("../EC.jl")
|
||||
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
|
|
@ -31,62 +30,16 @@ println("Basis size = ", basis.dim)
|
|||
|
||||
@time "Ha" Ha = T1 + T2 + Va
|
||||
@time "Eigenvalues" target_evals, _ = eigs(Ha, nev=5, ncv=50, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
|
||||
|
||||
display(target_evals)
|
||||
|
||||
# free memory
|
||||
basis = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
|
||||
GC.gc()
|
||||
|
||||
training_c = [-0.5, -0.65, -0.8, -1, -1.2]
|
||||
extrapolating_c = [0.8, 0.6, 0.4, 0.2, 0.1, 0.0, -0.1, -0.2, -0.3]
|
||||
|
||||
current_E = -0.5173809356244544
|
||||
ref_E = -0.5173809356244544
|
||||
|
||||
exact = ComplexF64[]
|
||||
training = ComplexF64[]
|
||||
extrapolated = ComplexF64[]
|
||||
training_vecs = Vector{ComplexF64}[]
|
||||
EC = affine_EC(Ha, Vb)
|
||||
train!(EC, training_c; ref_eval=ref_E, CAEC=true) # try CAEC=false !!!
|
||||
extrapolate!(EC, extrapolating_c)
|
||||
|
||||
for c in training_c
|
||||
print("Training for c = $c: ")
|
||||
H = Ha + c .* Vb
|
||||
evals, evecs = eigs(H, nev=3, ncv=24, which=:SR, maxiter=5000, tol=1e-5, ritzvec=true, check=1)
|
||||
|
||||
global current_E = nearest(evals, current_E)
|
||||
println(current_E)
|
||||
push!(training, current_E)
|
||||
push!(training_vecs, evecs[:, nearestIndex(evals, current_E)])
|
||||
end
|
||||
|
||||
training_vecs = vcat(training_vecs, conj(training_vecs)) # CA-EC
|
||||
println("Original EC dimensionality = $(length(training_vecs))")
|
||||
@time "Gram-Schmidt" training_vecs = gram_schmidt!(training_vecs; verbose=true) # orthonormalization
|
||||
|
||||
EC_basis = hcat(training_vecs...)
|
||||
Ha_EC = transpose(EC_basis) * Ha * EC_basis
|
||||
Vb_EC = transpose(EC_basis) * Vb * EC_basis
|
||||
|
||||
current_E = -0.3005521915662689 - 0.13612069020686351im
|
||||
|
||||
for c in extrapolating_c
|
||||
print("Extrapolating for c = $c: ")
|
||||
H = Ha + c .* Vb
|
||||
evals, evecs = eigs(H, nev=3, ncv=24, which=:SR, maxiter=5000, tol=1e-5, ritzvec=true, check=1)
|
||||
|
||||
global current_E = nearest(evals, current_E)
|
||||
println(current_E)
|
||||
push!(exact, current_E)
|
||||
|
||||
# extrapolation
|
||||
H_EC = Ha_EC + c .* Vb_EC
|
||||
evals = eigvals(H_EC)
|
||||
push!(extrapolated, nearest(evals, current_E))
|
||||
end
|
||||
|
||||
# exportCSV("temp/dis_HO_B2R.csv", (training, exact, extrapolated), ("training", "exact", "extrapolated"))
|
||||
|
||||
scatter(real.(training),imag.(training), label="training")
|
||||
scatter!(real.(exact),imag.(exact), label="exact")
|
||||
scatter!(real.(extrapolated),imag.(extrapolated), label="extrapolated")
|
||||
savefig("temp/dis_HO_B2R.pdf")
|
||||
exportCSV(EC, "temp/dis_HO_B2R.csv")
|
||||
plot(EC, "temp/dis_HO_B2R.pdf")
|
||||
Loading…
Reference in New Issue