Compare commits
1 Commits
main
...
k_moshinsk
| Author | SHA1 | Date |
|---|---|---|
 Nuwan Yapa
|
8b3eec01e8 |
173
EC.jl
173
EC.jl
|
|
@ -1,173 +0,0 @@
|
|||
using Statistics, SparseArrays, LinearAlgebra, Arpack, Plots
|
||||
include("common.jl")
|
||||
|
||||
"EC model for a Hamiltonian family H(c) = H0 + c * H1"
|
||||
mutable struct affine_EC
|
||||
H0::AbstractMatrix{ComplexF64}
|
||||
H1::AbstractMatrix{ComplexF64}
|
||||
weights::Vector{ComplexF64}
|
||||
trained::Bool
|
||||
|
||||
H0_EC
|
||||
H1_EC
|
||||
N_EC
|
||||
|
||||
ensemble_size::Int
|
||||
H0_EC_ensemble
|
||||
H1_EC_ensemble
|
||||
N_EC_ensemble
|
||||
|
||||
training_E::Vector{ComplexF64}
|
||||
exact_E::Vector{ComplexF64}
|
||||
extrapolated_E::Vector{ComplexF64}
|
||||
extrapolated_CI::Vector{ComplexF64}
|
||||
|
||||
affine_EC(H0::AbstractMatrix{ComplexF64}, H1::AbstractMatrix{ComplexF64}, weights::Vector{ComplexF64}=ones(ComplexF64, size(H0, 1)); ensemble_size=0) =
|
||||
new(H0, H1, weights, false, nothing, nothing, nothing, ensemble_size, Matrix[], Matrix[], Matrix[], ComplexF64[], ComplexF64[], ComplexF64[], ComplexF64[])
|
||||
end
|
||||
|
||||
"Train an 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 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.
|
||||
If orthonormalize_threshold > 0, Gram-Schmidt orthonormalization is performed, using this value as the threshold for dropping redundant vectors."
|
||||
function train!(EC::affine_EC, c_vals; ref_eval=-10.0, CAEC=false, gram_schmidt_threshold=0, verbose=true, tol=1e-5)
|
||||
training_vecs = Vector{ComplexF64}[]
|
||||
|
||||
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 = popfirst!(ref_eval)
|
||||
end
|
||||
|
||||
H = EC.H0 + c .* EC.H1
|
||||
|
||||
evals, evecs = eigs(H, sigma=current_E, maxiter=5000, tol=tol, ritzvec=true, check=1)
|
||||
current_E = nearest(evals, current_E)
|
||||
|
||||
push!(EC.training_E, current_E)
|
||||
push!(training_vecs, evecs[:, nearestIndex(evals, current_E)])
|
||||
end
|
||||
|
||||
CAEC && append!(training_vecs, conj.(training_vecs))
|
||||
|
||||
(EC.H0_EC, EC.H1_EC, EC.N_EC) = get_reduced_matrices(EC, training_vecs, gram_schmidt_threshold; verbose=true)
|
||||
|
||||
for _ in 1:EC.ensemble_size
|
||||
subsample = resample(length(training_vecs))
|
||||
if gram_schmidt_threshold > 0
|
||||
(H0_EC, H1_EC, N_EC) = get_reduced_matrices(EC, training_vecs, gram_schmidt_threshold, subsample; verbose=false)
|
||||
push!(EC.H0_EC_ensemble, H0_EC)
|
||||
push!(EC.H1_EC_ensemble, H1_EC)
|
||||
push!(EC.N_EC_ensemble, N_EC)
|
||||
else
|
||||
push!(EC.H0_EC_ensemble, EC.H0_EC[subsample, subsample])
|
||||
push!(EC.H1_EC_ensemble, EC.H1_EC[subsample, subsample])
|
||||
push!(EC.N_EC_ensemble, EC.N_EC[subsample, subsample])
|
||||
end
|
||||
end
|
||||
|
||||
EC.trained = true
|
||||
end
|
||||
|
||||
function get_reduced_matrices(EC::affine_EC, training_vecs, gram_schmidt_threshold, subsample=1:length(training_vecs); verbose=false)
|
||||
vecs = deepcopy(training_vecs[subsample])
|
||||
|
||||
if gram_schmidt_threshold > 0; vecs = gram_schmidt!(vecs, EC.weights, gram_schmidt_threshold; verbose=verbose); end
|
||||
|
||||
EC_basis = hcat(vecs...)
|
||||
weights_mat = spdiagm(EC.weights)
|
||||
H0_EC = transpose(EC_basis) * weights_mat * EC.H0 * EC_basis
|
||||
H1_EC = transpose(EC_basis) * weights_mat * EC.H1 * EC_basis
|
||||
N_EC = transpose(EC_basis) * weights_mat * EC_basis
|
||||
|
||||
return (H0_EC, H1_EC, N_EC)
|
||||
end
|
||||
|
||||
resample(n::Int) = rand(1:n, n) |> unique |> sort
|
||||
|
||||
"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.
|
||||
If precalculated_exact_E is provided, ref_eval is ignored.
|
||||
If pseudo_inv_tol > 0, the GEVP is avoided using Moore-Penrose psuedoinverse, using this value as the relative tolerance for dropping redundant vectors."
|
||||
function extrapolate!(EC::affine_EC, c_vals; ref_eval=EC.training_E[end], pseudo_inv_tol=0, verbose=true, tol=1e-5, precalculated_exact_E=nothing)
|
||||
@assert EC.trained "EC model must be trained using train() before extrapolation"
|
||||
|
||||
for c in c_vals
|
||||
global current_E
|
||||
|
||||
if isnothing(precalculated_exact_E)
|
||||
if ref_eval isa Number
|
||||
current_E = ref_eval
|
||||
ref_eval = nothing
|
||||
elseif !isnothing(ref_eval)
|
||||
current_E = popfirst!(ref_eval)
|
||||
end
|
||||
|
||||
verbose && println("Extact solution for c = $c")
|
||||
|
||||
H = EC.H0 + c .* EC.H1
|
||||
evals, _ = eigs(H, sigma=current_E, maxiter=5000, tol=tol, ritzvec=false, check=1)
|
||||
current_E = nearest(evals, current_E)
|
||||
else
|
||||
current_E = popfirst!(precalculated_exact_E)
|
||||
end
|
||||
|
||||
push!(EC.exact_E, current_E)
|
||||
|
||||
verbose && println("Extrapolating for c = $c")
|
||||
|
||||
evals = get_extrapolated_evals(EC.H0_EC, EC.H1_EC, EC.N_EC, c, pseudo_inv_tol)
|
||||
push!(EC.extrapolated_E, nearest(evals, current_E))
|
||||
|
||||
if EC.ensemble_size > 0
|
||||
E_ensemble = zeros(ComplexF64, EC.ensemble_size)
|
||||
for i in 1:EC.ensemble_size
|
||||
evals = get_extrapolated_evals(EC.H0_EC_ensemble[i], EC.H1_EC_ensemble[i], EC.N_EC_ensemble[i], c, pseudo_inv_tol)
|
||||
E_ensemble[i] = nearest(evals, current_E)
|
||||
end
|
||||
re_CI = std(real.(E_ensemble))
|
||||
im_CI = std(imag.(E_ensemble))
|
||||
push!(EC.extrapolated_CI, complex(re_CI, im_CI))
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
"Solve the GEVP with or without Moore-Penrose psuedoinverse"
|
||||
function get_extrapolated_evals(H0_EC, H1_EC, N_EC, c, pseudo_inv_tol)
|
||||
H_EC = H0_EC + c .* H1_EC
|
||||
if pseudo_inv_tol > 0
|
||||
inv_N_EC = pinv(N_EC; atol=pseudo_inv_tol)
|
||||
H_EC = inv_N_EC * H_EC
|
||||
return eigvals(H_EC)
|
||||
else
|
||||
return eigvals(H_EC, N_EC)
|
||||
end
|
||||
end
|
||||
|
||||
"Export EC data as CSV file"
|
||||
exportCSV(EC::affine_EC, filename) = exportCSV(filename, (EC.training_E, EC.exact_E, EC.extrapolated_E), ("training", "exact", "extrapolated"))
|
||||
|
||||
"Plot EC data and optionally save figure to a file"
|
||||
function plot(EC::affine_EC, save_fig_filename=nothing; basis_points=nothing, basis_contour=nothing, xlims=nothing, ylims=nothing)
|
||||
scatter(real.(EC.training_E), imag.(EC.training_E), label="training")
|
||||
scatter!(real.(EC.exact_E), imag.(EC.exact_E), label="exact", markercolor=:white)
|
||||
if EC.ensemble_size > 0
|
||||
scatter!(real.(EC.extrapolated_E), imag.(EC.extrapolated_E), xerror=real.(EC.extrapolated_CI), yerror=imag.(EC.extrapolated_CI), label="extrapolated", m=:x)
|
||||
else
|
||||
scatter!(real.(EC.extrapolated_E), imag.(EC.extrapolated_E), label="extrapolated", m=:x)
|
||||
end
|
||||
|
||||
isnothing(basis_points) || scatter!(real.(basis_points), imag.(basis_points), m=:x, label="basis")
|
||||
isnothing(basis_contour) || plot!(real.(basis_contour), imag.(basis_contour), label="contour")
|
||||
|
||||
isnothing(xlims) || xlims!(xlims...)
|
||||
isnothing(ylims) || ylims!(ylims...)
|
||||
|
||||
isnothing(save_fig_filename) || savefig(save_fig_filename)
|
||||
end
|
||||
|
|
@ -0,0 +1,100 @@
|
|||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"using Plots\n",
|
||||
"include(\"p_space.jl\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"vertices = [0, 0.4 - 0.15im, 0.8, 6]\n",
|
||||
"subdivisions = [128, 128, 128]\n",
|
||||
"p, w = get_mesh(vertices, subdivisions)\n",
|
||||
"mesh_E = p.*p ./ (2*0.5)\n",
|
||||
"\n",
|
||||
"# ResonanceEC: Eq. (20)\n",
|
||||
"V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"training_points = range(0.75, 0.45, 5)\n",
|
||||
"training_E = Vector{ComplexF64}(undef, length(training_points))\n",
|
||||
"EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points))\n",
|
||||
"\n",
|
||||
"for (j, c) in enumerate(training_points)\n",
|
||||
" evals, evecs = eigen(get_H_matrix(V_system(c), p, w))\n",
|
||||
" i = identify_pole_i(p, evals)\n",
|
||||
" training_E[j] = evals[i]\n",
|
||||
" EC_basis[:, j] = evecs[:, i]\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"scatter(real.(training_E), imag.(training_E), label=\"training\")\n",
|
||||
"plot!(real.(mesh_E), imag.(mesh_E), label=\"contour\")\n",
|
||||
"xlims!(0,1)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"extrapolate_points = range(0.40, 0.25, 5)\n",
|
||||
"ref_E = 0.2 - 0.1im\n",
|
||||
"\n",
|
||||
"exact_E = Vector{ComplexF64}(undef, length(extrapolate_points))\n",
|
||||
"extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points))\n",
|
||||
"\n",
|
||||
"EC_basis_w = EC_basis .* w\n",
|
||||
"N_EC = transpose(EC_basis_w) * EC_basis\n",
|
||||
"\n",
|
||||
"for (j, c) in enumerate(extrapolate_points)\n",
|
||||
" exact_E[j] = quick_pole_E(V_system(c))\n",
|
||||
"\n",
|
||||
" EC_basis_w = EC_basis .* w\n",
|
||||
" H = get_H_matrix(V_system(c), p, w)\n",
|
||||
" H_EC = transpose(EC_basis_w) * H * EC_basis\n",
|
||||
" evals = eigvals(H_EC, N_EC)\n",
|
||||
" i = argmin(abs.(evals .- ref_E))\n",
|
||||
" ref_E = evals[i]\n",
|
||||
" extrapolate_E[j] = evals[i]\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"scatter(real.(training_E), imag.(training_E), label=\"training\")\n",
|
||||
"scatter!(real.(exact_E), imag.(exact_E), label=\"exact\")\n",
|
||||
"scatter!(real.(extrapolate_E), imag.(extrapolate_E), label=\"extrapolated\")\n",
|
||||
"plot!(real.(mesh_E), imag.(mesh_E), label=\"contour\")\n",
|
||||
"xlims!(0,1)"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Julia 1.10.2",
|
||||
"language": "julia",
|
||||
"name": "julia-1.10"
|
||||
},
|
||||
"language_info": {
|
||||
"file_extension": ".jl",
|
||||
"mimetype": "application/julia",
|
||||
"name": "julia",
|
||||
"version": "1.10.2"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 2
|
||||
}
|
||||
|
|
@ -1,7 +1,5 @@
|
|||
[deps]
|
||||
Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
|
||||
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
|
||||
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
|
||||
DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
|
||||
FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
|
||||
HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
|
||||
|
|
@ -9,7 +7,6 @@ LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
|
|||
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
|
||||
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
|
||||
QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
|
||||
Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
|
||||
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
|
||||
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
|
||||
WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b"
|
||||
|
|
|
|||
44
berggren.jl
44
berggren.jl
|
|
@ -1,22 +1,30 @@
|
|||
using SparseArrays
|
||||
include("math.jl")
|
||||
include("ho_basis.jl")
|
||||
|
||||
"berg_bases1/2 are lists of (1+l_max) matrices containing the eigenbases corresponding to 1st and 2nd DOFs respectively,
|
||||
js are a list of tuples (j1, j2) corresponding to 1st and 2nd DOFs respectively,
|
||||
and ws are the weights needed to evaluate the inner products"
|
||||
function get_2p_p1p2_matrix(mesh_size, js, Λ, berg_bases1::Vector{Matrix{ComplexF64}}, berg_bases2::Vector{Matrix{ComplexF64}}, ks::Vector{ComplexF64}, ws::Vector{ComplexF64})
|
||||
# TODO: Cache / precalculate
|
||||
integral1(np, lp, n, l) = sum(ks .* berg_bases1[1 + lp][:, np] .* ws .* berg_bases1[1 + l][:, n])
|
||||
integral2(np, lp, n, l) = sum(ks .* berg_bases2[1 + lp][:, np] .* ws .* berg_bases2[1 + l][:, n])
|
||||
|
||||
basis = iter_prod(js, 1:mesh_size, 1:mesh_size)
|
||||
mat = zeros(ComplexF64, length(basis), length(basis))
|
||||
Threads.@threads for idx in CartesianIndices(mat)
|
||||
(ip, i) = Tuple(idx)
|
||||
((j1p, j2p), n1p, n2p) = basis[ip]
|
||||
((j1, j2), n1, n2) = basis[i]
|
||||
val = racahs_reduction_formula(n1p, j1p, n2p, j2p, n1, j1, n2, j2, Λ, integral1, integral2)
|
||||
if !(val ≈ 0); mat[idx] = val; end
|
||||
"Basis transformation from HO to momentum space"
|
||||
function get_W_matrix(basis_p, E_max, Λ, μ1ω1, μ2ω2=μ1ω1; weights=true)
|
||||
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
|
||||
W = zeros(ComplexF64, length(basis_p), length(Es))
|
||||
Threads.@threads for idx in CartesianIndices(W)
|
||||
(i1, i2) = Tuple(idx)
|
||||
((j1, j2), (k1, w1), (k2, w2)) = basis_p[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)
|
||||
elem2 = 1/sqrt(sqrt(μ2ω2)) * (-1)^n2s[i2] * ho_basis(j2, n2s[i2], 1/sqrt(μ2ω2) * k2)
|
||||
W[idx] = elem1 * elem2 * (weights ? w1 * w2 : 1)
|
||||
end
|
||||
end
|
||||
return sparse(mat)
|
||||
return sparse(W)
|
||||
end
|
||||
|
||||
"Moshinsky transform in k-space"
|
||||
function get_Moshinsky_in_kspace(E_max, basis_p, μ1ω1, μ2ω2, Λ)
|
||||
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
|
||||
U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
|
||||
|
||||
W_right = get_W_matrix(basis_p, E_max, Λ, μ1ω1, μ2ω2; weights=false)
|
||||
W_left = get_W_matrix(basis_p, E_max, Λ, μ1ω1, μ2ω2; weights=false)
|
||||
|
||||
out = W_left * U * transpose(W_right)
|
||||
return out
|
||||
end
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
using SparseArrays, Arpack
|
||||
include("common.jl")
|
||||
include("helper.jl")
|
||||
include("p_space.jl")
|
||||
|
||||
E_target = -0.3919
|
||||
|
|
@ -0,0 +1,63 @@
|
|||
using LinearAlgebra, SparseArrays, Arpack
|
||||
include("helper.jl")
|
||||
include("p_space.jl")
|
||||
include("berggren.jl")
|
||||
|
||||
println("No of threads = ", Threads.nthreads())
|
||||
|
||||
atol = 10^-5
|
||||
maxevals = 10^5
|
||||
R_cutoff = 16
|
||||
|
||||
Λ = 0
|
||||
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=R_cutoff)
|
||||
|
||||
vertices = [0, 0.5 - 0.3im, 1, 4]
|
||||
subdivisions = [10, 10, 10]
|
||||
ks, ws = get_mesh(vertices, subdivisions)
|
||||
|
||||
jmax = 4
|
||||
tri((j1, j2)) = triangle_ineq(j1, j2, Λ)
|
||||
js = collect(Iterators.filter(tri, iter_prod(0:jmax, 0:jmax)))
|
||||
|
||||
basis = iter_prod(js, zip(ks, ws), zip(ks, ws)) # basis = ((j1, j2), (k1, w1), (k2, w2))
|
||||
basis_size = length(js) * length(ks)^2
|
||||
@assert length(basis) == basis_size "Something wrong with the basis"
|
||||
println("Basis size = $basis_size")
|
||||
|
||||
@time "T" begin
|
||||
T_blocks = [kron_sum(get_T_matrix(ks, μ1), get_T_matrix(ks, μ2)) for _ in js]
|
||||
T = blockdiag(sparse.(T_blocks)...)
|
||||
end
|
||||
|
||||
@time "V1" begin
|
||||
V1_blocks = [kron(get_V_matrix((k, kp) -> V_l(j1, k, kp), ks, ws), I(length(ks))) for (j1, _) in js]
|
||||
V1 = blockdiag(sparse.(V1_blocks)...)
|
||||
end
|
||||
|
||||
@time "V_relative" begin
|
||||
V_relative_blocks = [kron(get_V_matrix((k, kp) -> V_l(j1, k, kp), ks, ws), get_V_matrix((k, kp) -> V_l(j2, k, kp), ks, ws)) for (j1, j2) in js]
|
||||
V_relative = blockdiag(sparse.(V_relative_blocks)...)
|
||||
end
|
||||
|
||||
@time "weights" weights = spdiagm(vcat(fill(kron(ws, ws), length(js))...))
|
||||
|
||||
# Moshinsky transform
|
||||
E_max = 40
|
||||
μω_global = 0.5
|
||||
μ1ω1 = μω_global * 1/2
|
||||
μ2ω2 = μω_global * 2
|
||||
@time "M" M = get_Moshinsky_in_kspace(E_max, basis, μ1ω1, μ2ω2, Λ)
|
||||
@time "V2" V2 = M' * weights * V_relative * M * weights
|
||||
|
||||
@time "H" H = T + V1 + V2
|
||||
|
||||
@time "Eigenvalues" evals, _ = eigs(H, nev=3, ncv=24, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
|
||||
display(evals)
|
||||
|
|
@ -1,7 +1,6 @@
|
|||
using LinearAlgebra, SparseArrays, Arpack
|
||||
include("common.jl")
|
||||
include("helper.jl")
|
||||
include("p_space.jl")
|
||||
include("ho_basis.jl")
|
||||
include("berggren.jl")
|
||||
|
||||
println("No of threads = ", Threads.nthreads())
|
||||
|
|
@ -11,7 +10,8 @@ R_cutoff = 16
|
|||
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
μ = m/2 # due to simple relative coordinates
|
||||
μ1 = m * 1/2
|
||||
μ2 = m * 2/3
|
||||
|
||||
target = 4.0766890719636875 - 0.012758927741074495im
|
||||
|
||||
|
|
@ -31,45 +31,28 @@ basis_size = length(js) * length(ks)^2
|
|||
@assert length(basis) == basis_size "Something wrong with the basis"
|
||||
println("Basis size = $basis_size")
|
||||
|
||||
# generate Berggren bases
|
||||
@time "berg_bases" begin
|
||||
berg_bases = Vector{Matrix{ComplexF64}}(undef, jmax + 1)
|
||||
for j in 0:jmax
|
||||
_, berg_basis = eigen(get_H_matrix((k, kp) -> V_l(j, k, kp), ks, ws, μ); permute=false, scale=false)
|
||||
N_berg = sum(berg_basis.^2 .* ws, dims=1)
|
||||
berg_bases[1 + j] = berg_basis ./ transpose(sqrt.(N_berg))
|
||||
end
|
||||
to_berg_basis(mat, j) = transpose(berg_bases[1 + j] .* ws) * mat * berg_bases[1 + j]
|
||||
@time "T" begin
|
||||
T_blocks = [kron_sum(get_T_matrix(ks, μ1), get_T_matrix(ks, μ2)) for _ in js]
|
||||
T = blockdiag(sparse.(T_blocks)...)
|
||||
end
|
||||
|
||||
@time "U_berggren" begin
|
||||
U_blocks = [kron(berg_bases[1 + j1], berg_bases[1 + j2]) for (j1, j2) in js]
|
||||
U = blockdiag(sparse.(U_blocks)...)
|
||||
@time "V1" begin
|
||||
V1_blocks = [kron(get_V_matrix((k, kp) -> V_l(j1, k, kp), ks, ws), I(length(ks))) for (j1, _) in js]
|
||||
V1 = blockdiag(sparse.(V1_blocks)...)
|
||||
end
|
||||
|
||||
@time "Block diagonal part" begin
|
||||
Hb_blocks = [kron_sum(to_berg_basis(get_H_matrix((k, kp) -> V_l(j1, k, kp), ks, ws, μ), j1), to_berg_basis(get_H_matrix((k, kp) -> V_l(j2, k, kp), ks, ws, μ), j2)) for (j1, j2) in js]
|
||||
Hb = blockdiag(sparse.(Hb_blocks)...)
|
||||
end
|
||||
|
||||
@time "T_cross" T_cross = get_2p_p1p2_matrix(length(ks), js, Λ, berg_bases, berg_bases, ks, ws) ./ (2*μ)
|
||||
|
||||
E_max = 30
|
||||
μω_global = 0.5
|
||||
# due to simple relative coordinates
|
||||
μω = μω_global * 2
|
||||
μ = m/2
|
||||
μ1ω1 = μω_global * 1/2
|
||||
μ2ω2 = μω_global * 2
|
||||
|
||||
basis_ho = ho_basis_2B(E_max, Λ)
|
||||
@time "V2_HO" V2_HO = get_jacobi_V2_matrix(V_of_r, E_max, Λ, μω_global)
|
||||
|
||||
@time "V12_HO" V12_HO = get_src_V12_matrix(V_of_r, basis_ho, μω_global; atol=10^-6, maxevals=10^5)
|
||||
@time "W_right" W_right = get_W_matrix(basis, E_max, Λ, μ1ω1, μ2ω2; weights=true)
|
||||
@time "W_left" W_left = get_W_matrix(basis, E_max, Λ, μ1ω1, μ2ω2; weights=false)
|
||||
|
||||
@time "W" W = get_W_matrix(basis, basis_ho, μω, μω; weights=true)
|
||||
|
||||
@time "V12_p" V12_p = W * V12_HO * transpose(W)
|
||||
@time "V12" V12 = transpose(U) * V12_p * U
|
||||
|
||||
@time "H" H = Hb + T_cross + V12
|
||||
@time "V2" V2 = W_left * V2_HO * transpose(W_right)
|
||||
@time "H" H = T + V1 + V2
|
||||
|
||||
@time "Eigenvalues" evals, _ = eigs(H, sigma=target, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
|
||||
display(evals)
|
||||
|
|
|
|||
|
|
@ -1,33 +0,0 @@
|
|||
using Roots, DelimitedFiles
|
||||
|
||||
include("../EC.jl")
|
||||
include("../common.jl")
|
||||
include("../p_space.jl")
|
||||
|
||||
μ = 0.5
|
||||
V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20)
|
||||
|
||||
# determining c0 with EC
|
||||
temp_c = range(1.1, 0.9, 3)
|
||||
p, w = get_mesh([0, 8], [256])
|
||||
H0 = get_T_matrix(p, μ)
|
||||
V = get_V_matrix(V_system(1), p, w)
|
||||
EC = affine_EC(H0, V, w)
|
||||
train!(EC, temp_c; ref_eval=-0.2, CAEC=false, verbose=false)
|
||||
quick_extrapolate(c) = minimum(abs2, get_extrapolated_evals(EC.H0_EC, EC.H1_EC, EC.N_EC, c, 0))
|
||||
c0 = find_zero(quick_extrapolate, 0.85)
|
||||
|
||||
training_c = range(1.2, 0.9, 9) # original: range(1.35, 0.9, 5)
|
||||
extrapolating_c = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
|
||||
|
||||
data_c = vcat(training_c, extrapolating_c)
|
||||
data_E = [quick_pole_E(V_system(c)) for c in data_c]
|
||||
|
||||
# export to CSV
|
||||
file = "temp/2body_data.csv"
|
||||
delim = ','
|
||||
open(file, "w") do f
|
||||
writedlm(f, ["c" "re_E" "im_E"], delim)
|
||||
writedlm(f, [c0 0 0], delim) # first entry for the threshold
|
||||
writedlm(f, hcat(data_c, real.(data_E), imag.(data_E)), delim)
|
||||
end
|
||||
|
|
@ -0,0 +1,88 @@
|
|||
using Arpack, SparseArrays, LRUCache
|
||||
using DelimitedFiles, Plots
|
||||
include("../ho_basis.jl")
|
||||
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
V_of_r(r) = -5 * exp(-r^2 / 3) + 2 * exp(-r^2 / 10)
|
||||
|
||||
E_max = 40
|
||||
μω_global = 0.2 * exp(-2im * pi / 9)
|
||||
|
||||
# due to Jacobi coordinates
|
||||
μ1ω1 = μω_global * 1/2
|
||||
μ2ω2 = μω_global * 2
|
||||
μ1 = m * 1/2
|
||||
μ2 = m * 2/3
|
||||
|
||||
println("No of threads = ", Threads.nthreads())
|
||||
|
||||
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
|
||||
l_max = max(maximum(l1s), maximum(l2s))
|
||||
n_max = max(maximum(n1s), maximum(n2s))
|
||||
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
|
||||
mask2 = (n1s .== n1s') .&& (l1s .== l1s')
|
||||
|
||||
println("Basis size = ", length(Es))
|
||||
|
||||
@time "T1" T1 = get_sp_T_matrix(n1s, l1s; mask=mask1, μω_gen=μ1ω1, μ=μ1)
|
||||
@time "T2" T2 = get_sp_T_matrix(n2s, l2s; mask=mask2, μω_gen=μ2ω2, μ=μ2)
|
||||
|
||||
@time "V" Vb = get_jacobi_V_matrix(V_of_r, E_max, Λ, μ1ω1, μω_global)
|
||||
|
||||
@time "H0" Ha = T1 + T2
|
||||
|
||||
# free memory
|
||||
Es = n1s = l1s = n2s = l2s = mask1 = mask2 = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
|
||||
GC.gc()
|
||||
|
||||
current_E = -0.26141959851000807
|
||||
|
||||
training_c = range(1.1, 0.9, 5)
|
||||
extrapolating_c = range(0.78, 0.45, 7)
|
||||
|
||||
exact = ComplexF64[]
|
||||
training = ComplexF64[]
|
||||
extrapolated = ComplexF64[]
|
||||
training_vecs = Vector{ComplexF64}[]
|
||||
|
||||
for c in training_c
|
||||
println("Training for c = $c")
|
||||
H = Ha + c .* Vb
|
||||
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
|
||||
|
||||
# CA-EC
|
||||
training_vecs = vcat(training_vecs, conj(training_vecs))
|
||||
|
||||
EC_basis = hcat(training_vecs...)
|
||||
N_EC = transpose(EC_basis) * EC_basis
|
||||
Ha_EC = transpose(EC_basis) * Ha * EC_basis
|
||||
Vb_EC = transpose(EC_basis) * Vb * EC_basis
|
||||
|
||||
current_E = 0.05387926316280764 - 0.008900278183544552im
|
||||
|
||||
for c in extrapolating_c
|
||||
println("Extrapolating for c = $c")
|
||||
H = Ha + c .* Vb
|
||||
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 = Ha_EC + c .* Vb_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")
|
||||
xlims!(-0.3, 0.3)
|
||||
ylims!(-0.2, 0)
|
||||
savefig("temp/2in3_HO_B2R.pdf")
|
||||
|
|
@ -1,40 +1,136 @@
|
|||
using Arpack, SparseArrays, LRUCache
|
||||
using DelimitedFiles, Plots
|
||||
include("../ho_basis.jl")
|
||||
include("../p_space.jl")
|
||||
include("../EC.jl")
|
||||
include("../berggren.jl")
|
||||
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
|
||||
|
||||
ϕ = 0.1
|
||||
vertices = [0, 4 * exp(-1im * ϕ), 5, 6]
|
||||
subdivisions = [40, 12, 15]
|
||||
jmax = 6
|
||||
|
||||
E_max = 40
|
||||
μω_global = 0.5
|
||||
|
||||
H0, weights = get_3b_H_matrix(jacobi, V_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m)
|
||||
|
||||
# Vp = perturbation to make the state artificially bound
|
||||
Vp_of_r(r) = -exp(-(r/3)^2)
|
||||
Vp, _ = get_3b_H_matrix(jacobi, Vp_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m, false, true)
|
||||
println("No of threads = ", Threads.nthreads())
|
||||
|
||||
training_c = [2.6, 2.4, 2.2, 2.0, 1.8]
|
||||
extrapolating_c = 0.0 : 0.2 : 1.2
|
||||
|
||||
training_ref = -2.22 # complete list not needed because identification is simple
|
||||
|
||||
extrapolating_ref = [4.076662025307587-0.012709842443350328im,
|
||||
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]
|
||||
1.233088227541505-0.0003070320106485624im])
|
||||
|
||||
EC = affine_EC(H0, Vp, weights)
|
||||
train!(EC, training_c; ref_eval=training_ref, CAEC=true)
|
||||
extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref)
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
Va_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
|
||||
Vb_of_r(r) = -exp(-(r/3)^2)
|
||||
|
||||
exportCSV(EC, "temp/Berggren_B2R.csv")
|
||||
plot(EC, "temp/Berggren_B2R.pdf")
|
||||
atol = 10^-5
|
||||
maxevals = 10^5
|
||||
R_cutoff = 16
|
||||
|
||||
# due to Jacobi coordinates
|
||||
μ1 = m * 1/2
|
||||
μ2 = m * 2/3
|
||||
|
||||
vertices = [0, 2 - 0.2im, 3, 4]
|
||||
subdivisions = [16, 10, 10]
|
||||
ks, ws = get_mesh(vertices, subdivisions)
|
||||
|
||||
jmax = 4
|
||||
tri((j1, j2)) = triangle_ineq(j1, j2, Λ)
|
||||
js = collect(Iterators.filter(tri, iter_prod(0:jmax, 0:jmax)))
|
||||
|
||||
basis = iter_prod(js, zip(ks, ws), zip(ks, ws)) # basis = ((j1, j2), (k1, w1), (k2, w2))
|
||||
basis_size = length(js) * length(ks)^2
|
||||
weights_mat = spdiagm(repeat(kron(ws, ws), jmax + 1))
|
||||
@assert length(basis) == basis_size "Something wrong with the basis"
|
||||
println("Basis size = $basis_size")
|
||||
|
||||
@time "T" begin
|
||||
T_blocks = [kron_sum(get_T_matrix(ks, μ1), get_T_matrix(ks, μ2)) for _ in js]
|
||||
T = blockdiag(sparse.(T_blocks)...)
|
||||
end
|
||||
|
||||
@time "Va1" begin
|
||||
Va_l(j, k, kp) = Vl_mat_elem(Va_of_r, j, k, kp; atol=atol, maxevals=maxevals, R_cutoff=R_cutoff)
|
||||
Va1_blocks = [kron(get_V_matrix((k, kp) -> Va_l(j1, k, kp), ks, ws), I(length(ks))) for (j1, _) in js]
|
||||
Va1 = blockdiag(sparse.(Va1_blocks)...)
|
||||
end
|
||||
|
||||
@time "Vb1" begin
|
||||
Vb_l(j, k, kp) = Vl_mat_elem(Vb_of_r, j, k, kp; atol=atol, maxevals=maxevals, R_cutoff=R_cutoff)
|
||||
Vb1_blocks = [kron(get_V_matrix((k, kp) -> Vb_l(j1, k, kp), ks, ws), I(length(ks))) for (j1, _) in js]
|
||||
Vb1 = blockdiag(sparse.(Vb1_blocks)...)
|
||||
end
|
||||
|
||||
E_max = 40
|
||||
μω_global = 0.5
|
||||
μ1ω1 = μω_global * 1/2
|
||||
μ2ω2 = μω_global * 2
|
||||
|
||||
@time "Va2_HO" Va2_HO = get_jacobi_V2_matrix(Va_of_r, E_max, Λ, μω_global; atol=atol, maxevals=maxevals)
|
||||
@time "Vb2_HO" Vb2_HO = get_jacobi_V2_matrix(Vb_of_r, E_max, Λ, μω_global; atol=atol, maxevals=maxevals)
|
||||
|
||||
@time "W_right" W_right = get_W_matrix(basis, E_max, Λ, μ1ω1, μ2ω2; weights=true)
|
||||
@time "W_left" W_left = get_W_matrix(basis, E_max, Λ, μ1ω1, μ2ω2; weights=false)
|
||||
|
||||
@time "Va2" Va2 = W_left * Va2_HO * transpose(W_right)
|
||||
@time "Vb2" Vb2 = W_left * Vb2_HO * transpose(W_right)
|
||||
|
||||
@time "Ha" Ha = T + Va1 + Va2
|
||||
@time "Vb" Vb = Vb1 + Vb2
|
||||
@time "Eigenvalues" test_evals, _ = eigs(Ha, sigma=exact_ref[end], maxiter=5000, tol=1e-5, ritzvec=false, check=1)
|
||||
|
||||
display(test_evals)
|
||||
|
||||
# free memory
|
||||
Es = n1s = l1s = n2s = l2s = mask1 = mask2 = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
|
||||
GC.gc()
|
||||
|
||||
current_E = training_ref
|
||||
|
||||
exact = ComplexF64[]
|
||||
training = ComplexF64[]
|
||||
extrapolated = ComplexF64[]
|
||||
training_vecs = Vector{ComplexF64}[]
|
||||
|
||||
for c in training_c
|
||||
println("Training for c = $c")
|
||||
H = Ha + c .* Vb
|
||||
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
|
||||
|
||||
# CA-EC
|
||||
training_vecs = vcat(training_vecs, conj(training_vecs))
|
||||
|
||||
EC_basis = hcat(training_vecs...)
|
||||
N_EC = transpose(EC_basis) * weights_mat * EC_basis
|
||||
Ha_EC = transpose(EC_basis) * weights_mat * Ha * EC_basis
|
||||
Vb_EC = transpose(EC_basis) * weights_mat * Vb * EC_basis
|
||||
|
||||
for c in extrapolating_c
|
||||
println("Extrapolating for c = $c")
|
||||
global current_E = pop!(exact_ref)
|
||||
|
||||
H = Ha + c .* Vb
|
||||
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 = Ha_EC + c .* Vb_EC
|
||||
evals = eigvals(H_EC, N_EC)
|
||||
push!(extrapolated, nearest(evals, current_E))
|
||||
end
|
||||
|
||||
exportCSV("temp/Berggren_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/Berggren_B2R.pdf")
|
||||
|
|
@ -1,42 +1,134 @@
|
|||
using Arpack, SparseArrays, LRUCache
|
||||
using DelimitedFiles, Plots
|
||||
include("../ho_basis.jl")
|
||||
include("../p_space.jl")
|
||||
include("../EC.jl")
|
||||
include("../berggren.jl")
|
||||
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
|
||||
|
||||
vertices = [0, 2 - 0.2im, 3, 4]
|
||||
subdivisions = [15, 10, 10]
|
||||
jmax = 4
|
||||
|
||||
E_max = 40
|
||||
μω_global = 0.5
|
||||
|
||||
H0, weights = get_3b_H_matrix(jacobi, V_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m)
|
||||
|
||||
# Vp = perturbation to make the state artificially bound
|
||||
Vp_of_r(r) = -exp(-(r/3)^2)
|
||||
Vp, _ = get_3b_H_matrix(jacobi, Vp_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m, false, true)
|
||||
println("No of threads = ", Threads.nthreads())
|
||||
|
||||
training_c = [1.1, 0.9, 0.7, 0.5]
|
||||
extrapolating_c = 0.0 : 0.2 : 1.2
|
||||
|
||||
training_ref = [1.4750633616275919 - 0.0003021770706749637im
|
||||
training_ref = reverse([1.4750633616275919 - 0.0003021770706749637im
|
||||
1.9567078295375822 - 0.0007646829108872369im
|
||||
2.4351117758403076 - 0.001281037843108658im
|
||||
2.9096543462392357 - 0.002962488527470604im]
|
||||
2.9096543462392357 - 0.002962488527470604im])
|
||||
|
||||
extrapolating_ref = [4.076662025307587-0.012709842443350328im,
|
||||
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]
|
||||
1.233088227541505-0.0003070320106485624im])
|
||||
|
||||
EC = affine_EC(H0, Vp, weights)
|
||||
train!(EC, training_c; ref_eval=training_ref, CAEC=false)
|
||||
extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref)
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
Va_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
|
||||
Vb_of_r(r) = -exp(-(r/3)^2)
|
||||
|
||||
exportCSV(EC, "temp/Berggren_R2R.csv")
|
||||
plot(EC, "temp/Berggren_R2R.pdf")
|
||||
atol = 10^-5
|
||||
maxevals = 10^5
|
||||
R_cutoff = 16
|
||||
|
||||
# due to Jacobi coordinates
|
||||
μ1 = m * 1/2
|
||||
μ2 = m * 2/3
|
||||
|
||||
vertices = [0, 2 - 0.1im, 3, 4]
|
||||
subdivisions = [16, 10, 10]
|
||||
ks, ws = get_mesh(vertices, subdivisions)
|
||||
|
||||
jmax = 4
|
||||
tri((j1, j2)) = triangle_ineq(j1, j2, Λ)
|
||||
js = collect(Iterators.filter(tri, iter_prod(0:jmax, 0:jmax)))
|
||||
|
||||
basis = iter_prod(js, zip(ks, ws), zip(ks, ws)) # basis = ((j1, j2), (k1, w1), (k2, w2))
|
||||
basis_size = length(js) * length(ks)^2
|
||||
weights_mat = spdiagm(repeat(kron(ws, ws), jmax + 1))
|
||||
@assert length(basis) == basis_size "Something wrong with the basis"
|
||||
println("Basis size = $basis_size")
|
||||
|
||||
@time "T" begin
|
||||
T_blocks = [kron_sum(get_T_matrix(ks, μ1), get_T_matrix(ks, μ2)) for _ in js]
|
||||
T = blockdiag(sparse.(T_blocks)...)
|
||||
end
|
||||
|
||||
@time "Va1" begin
|
||||
Va_l(j, k, kp) = Vl_mat_elem(Va_of_r, j, k, kp; atol=atol, maxevals=maxevals, R_cutoff=R_cutoff)
|
||||
Va1_blocks = [kron(get_V_matrix((k, kp) -> Va_l(j1, k, kp), ks, ws), I(length(ks))) for (j1, _) in js]
|
||||
Va1 = blockdiag(sparse.(Va1_blocks)...)
|
||||
end
|
||||
|
||||
@time "Vb1" begin
|
||||
Vb_l(j, k, kp) = Vl_mat_elem(Vb_of_r, j, k, kp; atol=atol, maxevals=maxevals, R_cutoff=R_cutoff)
|
||||
Vb1_blocks = [kron(get_V_matrix((k, kp) -> Vb_l(j1, k, kp), ks, ws), I(length(ks))) for (j1, _) in js]
|
||||
Vb1 = blockdiag(sparse.(Vb1_blocks)...)
|
||||
end
|
||||
|
||||
E_max = 40
|
||||
μω_global = 0.5
|
||||
μ1ω1 = μω_global * 1/2
|
||||
μ2ω2 = μω_global * 2
|
||||
|
||||
@time "Va2_HO" Va2_HO = get_jacobi_V2_matrix(Va_of_r, E_max, Λ, μω_global; atol=atol, maxevals=maxevals)
|
||||
@time "Vb2_HO" Vb2_HO = get_jacobi_V2_matrix(Vb_of_r, E_max, Λ, μω_global; atol=atol, maxevals=maxevals)
|
||||
|
||||
@time "W_right" W_right = get_W_matrix(basis, E_max, Λ, μ1ω1, μ2ω2; weights=true)
|
||||
@time "W_left" W_left = get_W_matrix(basis, E_max, Λ, μ1ω1, μ2ω2; weights=false)
|
||||
|
||||
@time "Va2" Va2 = W_left * Va2_HO * transpose(W_right)
|
||||
@time "Vb2" Vb2 = W_left * Vb2_HO * transpose(W_right)
|
||||
|
||||
@time "Ha" Ha = T + Va1 + Va2
|
||||
@time "Vb" Vb = Vb1 + Vb2
|
||||
@time "Eigenvalues" test_evals, _ = eigs(Ha, sigma=exact_ref[end], maxiter=5000, tol=1e-5, ritzvec=false, check=1)
|
||||
|
||||
display(test_evals)
|
||||
|
||||
# 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")
|
||||
global current_E = pop!(training_ref)
|
||||
|
||||
H = Ha + c .* Vb
|
||||
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
|
||||
|
||||
EC_basis = hcat(training_vecs...)
|
||||
N_EC = transpose(EC_basis) * weights_mat * EC_basis
|
||||
Ha_EC = transpose(EC_basis) * weights_mat * Ha * EC_basis
|
||||
Vb_EC = transpose(EC_basis) * weights_mat * Vb * EC_basis
|
||||
|
||||
for c in extrapolating_c
|
||||
println("Extrapolating for c = $c")
|
||||
global current_E = pop!(exact_ref)
|
||||
|
||||
H = Ha + c .* Vb
|
||||
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 = Ha_EC + c .* Vb_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/Berggren_R2R.pdf")
|
||||
|
|
@ -1,40 +0,0 @@
|
|||
include("../p_space.jl")
|
||||
include("../EC.jl")
|
||||
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
|
||||
|
||||
ϕ = 0.1
|
||||
vertices = [0, 6 * exp(-1im * ϕ)]
|
||||
subdivisions = [50]
|
||||
jmax = 4
|
||||
|
||||
E_max = 40
|
||||
μω_global = 0.5
|
||||
|
||||
H0, weights = get_3b_H_matrix(jacobi, V_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m)
|
||||
|
||||
# Vp = perturbation to make the state artificially bound
|
||||
Vp_of_r(r) = -exp(-(r/3)^2)
|
||||
Vp, _ = get_3b_H_matrix(jacobi, Vp_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m, false, true)
|
||||
|
||||
training_c = [2.6, 2.4, 2.2, 2.0, 1.8]
|
||||
extrapolating_c = 0.0 : 0.2 : 1.2
|
||||
|
||||
training_ref = -2.22 # complete list not needed because identification is simple
|
||||
|
||||
exact_E = [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]
|
||||
|
||||
EC = affine_EC(H0, Vp, weights)
|
||||
train!(EC, training_c; ref_eval=training_ref, CAEC=true)
|
||||
extrapolate!(EC, extrapolating_c; precalculated_exact_E=exact_E)
|
||||
|
||||
exportCSV(EC, "temp/CSM_B2R.csv")
|
||||
plot(EC, "temp/CSM_B2R.pdf")
|
||||
|
|
@ -1,68 +0,0 @@
|
|||
using Roots, LinearAlgebra, Plots
|
||||
|
||||
include("../EC.jl")
|
||||
include("../common.jl")
|
||||
include("../ho_basis.jl")
|
||||
|
||||
V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
|
||||
ϕ = 0.1
|
||||
μω_global = 0.5 * exp(-2im * ϕ)
|
||||
E_max = 40
|
||||
|
||||
H0 = get_3b_H_matrix(jacobi, V_of_r, μω_global, E_max, Λ, m, true, true)
|
||||
|
||||
# Vp = perturbation to make the state artificially bound
|
||||
Vp_of_r(r) = -exp(-(r/3)^2)
|
||||
@time "Vp" Vp = get_3b_H_matrix(jacobi, Vp_of_r, μω_global, E_max, Λ, m, false, true)
|
||||
|
||||
training_ref = -2.22
|
||||
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]
|
||||
|
||||
training_c = range(2.8, 1.8, 5)
|
||||
extrapolating_c = 0.0 : 0.2 : 1.2
|
||||
|
||||
EC = affine_EC(H0, Vp)
|
||||
train!(EC, training_c; ref_eval=training_ref, CAEC=true)
|
||||
extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref)
|
||||
|
||||
# determining c0 with EC
|
||||
approx_c0 = 1.5
|
||||
quick_extrapolate(c) = minimum(abs2, get_extrapolated_evals(EC.H0_EC, EC.H1_EC, EC.N_EC, c, 1e-14))
|
||||
c0 = find_zero(quick_extrapolate, approx_c0)
|
||||
|
||||
order::Int = ceil((length(training_c) - 1) / 2) # order of the Pade approximant
|
||||
|
||||
# Solve coefficients as a linear system
|
||||
training_k = alt_sqrt.(EC.training_E)
|
||||
M_left_element(c, i) = alt_sqrt(c - c0)^i
|
||||
M_left = M_left_element.(training_c, (0:order)')
|
||||
M_right = -training_k .* M_left[:, 2:end] # remove the first column
|
||||
M = hcat(M_left, M_right) # M = [M_left | M_right]
|
||||
sol = M \ training_k
|
||||
a = sol[1:order+1]
|
||||
b = [1; sol[order+2:end]]
|
||||
|
||||
# Pade approximant
|
||||
polynomial(a, c) = sum(i -> a[i+1] * alt_sqrt(c - c0)^i, 0:order)
|
||||
pade_approx(c) = polynomial(a, c) / polynomial(b, c)
|
||||
|
||||
# Extrapolate
|
||||
extrapolated_k = pade_approx.([extrapolating_c; training_c])
|
||||
extrapolated_E = extrapolated_k .^ 2
|
||||
|
||||
# Plotting
|
||||
scatter(real.(EC.training_E), imag.(EC.training_E), label="training")
|
||||
scatter!(real.(EC.exact_E), imag.(EC.exact_E), label="exact")
|
||||
scatter!(real.(EC.extrapolated_E), imag.(EC.extrapolated_E), label="CAEC", m=:x)
|
||||
scatter!(real.(extrapolated_E), imag.(extrapolated_E), label="ACCC", m=:+)
|
||||
title!("3-body extrapolation with $(length(training_c)) training points")
|
||||
savefig("temp/3body_HO_B2R_ACCC-$(length(training_c)).pdf")
|
||||
|
|
@ -1,36 +1,93 @@
|
|||
using Arpack, SparseArrays, LRUCache
|
||||
using DelimitedFiles, Plots
|
||||
include("../ho_basis.jl")
|
||||
include("../EC.jl")
|
||||
|
||||
V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
Va_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
|
||||
Vb_of_r(r) = -exp(-(r/3)^2)
|
||||
|
||||
ϕ = 0.1
|
||||
μω_global = 0.5 * exp(-2im * ϕ)
|
||||
E_max = 40
|
||||
μω_global = 0.5 * exp(-2im * pi / 9)
|
||||
|
||||
H0 = get_3b_H_matrix(jacobi, V_of_r, μω_global, E_max, Λ, m, true, true)
|
||||
# due to Jacobi coordinates
|
||||
μ1ω1 = μω_global * 1/2
|
||||
μ2ω2 = μω_global * 2
|
||||
μ1 = m * 1/2
|
||||
μ2 = m * 2/3
|
||||
|
||||
# Vp = perturbation to make the state artificially bound
|
||||
Vp_of_r(r) = -exp(-(r/3)^2)
|
||||
@time "Vp" Vp = get_3b_H_matrix(jacobi, Vp_of_r, μω_global, E_max, Λ, m, false, true)
|
||||
println("No of threads = ", Threads.nthreads())
|
||||
|
||||
training_ref = -2.22
|
||||
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
|
||||
l_max = max(maximum(l1s), maximum(l2s))
|
||||
n_max = max(maximum(n1s), maximum(n2s))
|
||||
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
|
||||
mask2 = (n1s .== n1s') .&& (l1s .== l1s')
|
||||
|
||||
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]
|
||||
println("Basis size = ", length(Es))
|
||||
|
||||
training_c = [2.6, 2.4, 2.2, 2.0, 1.8]
|
||||
@time "T1" T1 = get_sp_T_matrix(n1s, l1s; mask=mask1, μω_gen=μ1ω1, μ=μ1)
|
||||
@time "T2" T2 = get_sp_T_matrix(n2s, l2s; mask=mask2, μω_gen=μ2ω2, μ=μ2)
|
||||
|
||||
@time "Va" Va = get_jacobi_V_matrix(Va_of_r, E_max, Λ, μ1ω1, μω_global)
|
||||
@time "Vb" Vb = get_jacobi_V_matrix(Vb_of_r, E_max, Λ, μ1ω1, μω_global)
|
||||
|
||||
@time "Ha" Ha = T1 + T2 + Va
|
||||
@time "Eigenvalues" target_evals, _ = eigs(Ha, nev=5, ncv=50, which=:LI, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
|
||||
|
||||
display(target_evals)
|
||||
|
||||
# free memory
|
||||
Es = n1s = l1s = n2s = l2s = mask1 = mask2 = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
|
||||
GC.gc()
|
||||
|
||||
current_E = -0.72763
|
||||
|
||||
training_c = [2.0, 1.9, 1.8]
|
||||
extrapolating_c = 0.0 : 0.2 : 1.2
|
||||
|
||||
EC = affine_EC(H0, Vp)
|
||||
train!(EC, training_c; ref_eval=training_ref, CAEC=true)
|
||||
extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref)
|
||||
exact = ComplexF64[]
|
||||
training = ComplexF64[]
|
||||
extrapolated = ComplexF64[]
|
||||
training_vecs = Vector{ComplexF64}[]
|
||||
|
||||
exportCSV(EC, "temp/HO_B2R.csv")
|
||||
plot(EC, "temp/HO_B2R.pdf")
|
||||
for c in training_c
|
||||
println("Training for c = $c")
|
||||
H = Ha + c .* Vb
|
||||
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
|
||||
|
||||
# CA-EC
|
||||
training_vecs = vcat(training_vecs, conj(training_vecs))
|
||||
|
||||
EC_basis = hcat(training_vecs...)
|
||||
N_EC = transpose(EC_basis) * EC_basis
|
||||
Ha_EC = transpose(EC_basis) * Ha * EC_basis
|
||||
Vb_EC = transpose(EC_basis) * Vb * EC_basis
|
||||
|
||||
current_E = 4.0766890719636635 - 0.01275892774109674im
|
||||
|
||||
for c in extrapolating_c
|
||||
println("Extrapolating for c = $c")
|
||||
H = Ha + c .* Vb
|
||||
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 = Ha_EC + c .* Vb_EC
|
||||
evals = eigvals(H_EC, N_EC)
|
||||
push!(extrapolated, nearest(evals, current_E))
|
||||
end
|
||||
|
||||
exportCSV("temp/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/HO_B2R.pdf")
|
||||
|
|
@ -1,24 +1,54 @@
|
|||
include("../ho_basis.jl")
|
||||
include("../EC.jl")
|
||||
using DelimitedFiles, Plots
|
||||
include("../ho_basis_3body_resonance.jl")
|
||||
|
||||
V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
|
||||
μω_global = 0.5 * exp(-2im * pi / 9)
|
||||
E_max = 40
|
||||
|
||||
T = get_3b_H_matrix(src, V_of_r, μω_global, E_max, Λ, m, true, false)
|
||||
V = get_3b_H_matrix(src, V_of_r, μω_global, E_max, Λ, m, false, true)
|
||||
|
||||
ref_E = 5.9673 - 0.0006im
|
||||
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]
|
||||
|
||||
EC = affine_EC(T, V)
|
||||
train!(EC, training_c; ref_eval=ref_E, CAEC=false)
|
||||
extrapolate!(EC, extrapolating_c)
|
||||
@time "H0" H0 = T1 + T2
|
||||
|
||||
exportCSV(EC, "temp/HO_R2R.csv")
|
||||
plot(EC, "temp/HO_R2R.pdf")
|
||||
# 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")
|
||||
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
|
||||
|
||||
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")
|
||||
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, N_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")
|
||||
|
|
@ -1,77 +0,0 @@
|
|||
include("../p_space.jl")
|
||||
include("../EC.jl")
|
||||
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
|
||||
# Distinguishable particles: V12 = bound and V13 & V23 = resonant
|
||||
Vsubsystem_of_r(r) = -2 * exp(-r^2/4)
|
||||
Vdecay_of_r(r) = -exp(-r^2 / 3) + exp(-r^2 / 10)
|
||||
|
||||
ϕ = 0.1
|
||||
vertices = [0, 6 * exp(-1im * ϕ)]
|
||||
subdivisions = [50]
|
||||
jmax = 4
|
||||
|
||||
E_max = 40
|
||||
μω_global = 0.4
|
||||
|
||||
# Jacobi coordinates
|
||||
μ1ω1 = μω_global * 1/2
|
||||
μ2ω2 = μω_global * 2
|
||||
μ1 = m * 1/2
|
||||
μ2 = m * 2/3
|
||||
|
||||
atol=10^-5; maxevals=10^5; R_cutoff=16; verbose=true;
|
||||
|
||||
###########
|
||||
|
||||
verbose && println("No of threads = ", Threads.nthreads())
|
||||
|
||||
V_l(j, k, kp) = Vl_mat_elem(Vsubsystem_of_r, j, k, kp; atol=atol, maxevals=maxevals, R_cutoff=R_cutoff)
|
||||
|
||||
ks, ws = get_mesh(vertices, subdivisions)
|
||||
weights = repeat(kron(ws, ws), jmax + 1)
|
||||
block_size = length(ks)
|
||||
tri((j1, j2)) = triangle_ineq(j1, j2, Λ)
|
||||
js = collect(Iterators.filter(tri, iter_prod(0:jmax, 0:jmax)))
|
||||
basis = iter_prod(js, zip(ks, ws), zip(ks, ws)) # basis = ((j1, j2), (k1, w1), (k2, w2))
|
||||
basis_size = length(js) * length(ks)^2
|
||||
@assert length(basis) == basis_size "Something wrong with the basis"
|
||||
verbose && println("Basis size = $basis_size")
|
||||
|
||||
@time "Block diagonal part" begin
|
||||
blocks = [kron_sum(get_H_matrix((k, kp) -> V_l(j1, k, kp), ks, ws, μ1), get_T_matrix(ks, μ2)) for (j1, _) in js]
|
||||
Ha = blockdiag(sparse.(blocks)...)
|
||||
end
|
||||
|
||||
basis_ho = ho_basis_2B(E_max, Λ)
|
||||
verbose && println("HO basis size = ", basis_ho.dim)
|
||||
|
||||
@time "V2_HO" V2_HO = get_jacobi_V2_matrix(Vdecay_of_r, basis_ho, μω_global)
|
||||
|
||||
@time "W_right" W_right = get_W_matrix(basis, basis_ho, μ1ω1, μ2ω2; weights=true)
|
||||
@time "W_left" W_left = get_W_matrix(basis, basis_ho, μ1ω1, μ2ω2; weights=false)
|
||||
|
||||
@time "V2" Vb = W_left * V2_HO * transpose(W_right)
|
||||
|
||||
###########
|
||||
|
||||
training_c = [-0.55, -0.7, -0.85, -1, -1.2]
|
||||
extrapolating_c = [0.2, 0.1, 0.0, -0.1, -0.2, -0.3]
|
||||
|
||||
ref_E = -0.5173809356244544
|
||||
|
||||
exact_ref = [-0.31360746615280954-0.07689284936870341im
|
||||
-0.3233372403877718-0.06011323914565665im
|
||||
-0.339615582074795-0.04239442037174759im
|
||||
-0.36333816534241997-0.02648721825958402im
|
||||
-0.39376650561322885-0.014382935339817332im
|
||||
-0.4299479825535172-0.006510710745123606im]
|
||||
|
||||
EC = affine_EC(Ha, Vb)
|
||||
train!(EC, training_c; ref_eval=ref_E, CAEC=true)
|
||||
extrapolate!(EC, extrapolating_c; precalculated_exact_E=exact_ref)
|
||||
|
||||
exportCSV(EC, "temp/dis_CSM_B2R.csv")
|
||||
plot(EC, "temp/dis_CSM_B2R.pdf")
|
||||
|
|
@ -1,5 +1,6 @@
|
|||
using Arpack, SparseArrays, LRUCache
|
||||
using DelimitedFiles, Plots
|
||||
include("../ho_basis.jl")
|
||||
include("../EC.jl")
|
||||
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
|
|
@ -19,34 +20,78 @@ E_max = 40
|
|||
|
||||
println("No of threads = ", Threads.nthreads())
|
||||
|
||||
basis = ho_basis_2B(E_max, Λ)
|
||||
println("Basis size = ", basis.dim)
|
||||
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
|
||||
l_max = max(maximum(l1s), maximum(l2s))
|
||||
n_max = max(maximum(n1s), maximum(n2s))
|
||||
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
|
||||
mask2 = (n1s .== n1s') .&& (l1s .== l1s')
|
||||
|
||||
@time "T1" T1 = get_sp_T_matrix(basis.n1s, basis.l1s, [basis.n2s, basis.l2s]; μω_gen=μ1ω1, μ=μ1)
|
||||
@time "T2" T2 = get_sp_T_matrix(basis.n2s, basis.l2s, [basis.n1s, basis.l1s]; μω_gen=μ2ω2, μ=μ2)
|
||||
println("Basis size = ", length(Es))
|
||||
|
||||
@time "Va" Va = get_jacobi_V1_matrix(Va_of_r, basis, μ1ω1)
|
||||
@time "Vb" Vb = get_jacobi_V2_matrix(Vb_of_r, basis, μω_global)
|
||||
@time "T1" T1 = get_sp_T_matrix(n1s, l1s; mask=mask1, μω_gen=μ1ω1, μ=μ1)
|
||||
@time "T2" T2 = get_sp_T_matrix(n2s, l2s; mask=mask2, μω_gen=μ2ω2, μ=μ2)
|
||||
|
||||
@time "Va" Va = get_jacobi_V1_matrix(Va_of_r, E_max, Λ, μ1ω1)
|
||||
@time "Vb" Vb = get_jacobi_V2_matrix(Vb_of_r, E_max, Λ, μω_global)
|
||||
|
||||
@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)
|
||||
|
||||
training_c = [-0.55, -0.7, -0.85, -1, -1.2]
|
||||
extrapolating_c = [0.2, 0.1, 0.0, -0.1, -0.2, -0.3]
|
||||
# free memory
|
||||
Es = n1s = l1s = n2s = l2s = mask1 = mask2 = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing
|
||||
GC.gc()
|
||||
|
||||
ref_E = -0.5173809356244544
|
||||
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]
|
||||
|
||||
exact_ref = [-0.3136074661528041-0.07689284936868852im
|
||||
-0.323337240387771-0.06011323914564878im
|
||||
-0.33961558207479553-0.04239442037174764im
|
||||
-0.3633381653424224-0.026487218259589693im
|
||||
-0.393766505613234-0.014382935339825854im
|
||||
-0.42994798255352606-0.006510710745131777im]
|
||||
current_E = -0.5173809356244544
|
||||
|
||||
EC = affine_EC(Ha, Vb)
|
||||
train!(EC, training_c; ref_eval=ref_E, CAEC=true) # try CAEC=false !!!
|
||||
extrapolate!(EC, extrapolating_c, precalculated_exact_E = exact_ref)
|
||||
exact = ComplexF64[]
|
||||
training = ComplexF64[]
|
||||
extrapolated = ComplexF64[]
|
||||
training_vecs = Vector{ComplexF64}[]
|
||||
|
||||
exportCSV(EC, "temp/dis_HO_B2R.csv")
|
||||
plot(EC, "temp/dis_HO_B2R.pdf")
|
||||
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
|
||||
|
||||
# CA-EC
|
||||
training_vecs = vcat(training_vecs, conj(training_vecs))
|
||||
|
||||
EC_basis = hcat(training_vecs...)
|
||||
N_EC = transpose(EC_basis) * EC_basis
|
||||
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, N_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")
|
||||
|
|
@ -1,50 +0,0 @@
|
|||
using Roots, LinearAlgebra, Plots
|
||||
|
||||
include("../EC.jl")
|
||||
include("../common.jl")
|
||||
include("../p_space.jl")
|
||||
|
||||
μ = 0.5
|
||||
V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20)
|
||||
|
||||
# determining c0 with EC
|
||||
temp_c = range(1.1, 0.9, 3)
|
||||
p, w = get_mesh([0, 8], [256])
|
||||
H0 = get_T_matrix(p, μ)
|
||||
V = get_V_matrix(V_system(1), p, w)
|
||||
EC = affine_EC(H0, V, w)
|
||||
train!(EC, temp_c; ref_eval=-0.2, CAEC=false, verbose=false)
|
||||
quick_extrapolate(c) = minimum(abs2, get_extrapolated_evals(EC.H0_EC, EC.H1_EC, EC.N_EC, c, 0))
|
||||
c0 = find_zero(quick_extrapolate, 0.85)
|
||||
|
||||
# Calculation of training and extrapolating E
|
||||
training_c = range(1.2, 0.9, 9) # original: range(1.35, 0.9, 5)
|
||||
training_E = [quick_pole_E(V_system(c)) for c in training_c]
|
||||
training_k = alt_sqrt.(2μ .* training_E)
|
||||
|
||||
extrapolating_c = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
|
||||
exact_E = [quick_pole_E(V_system(c)) for c in extrapolating_c]
|
||||
|
||||
order::Int = ceil((length(training_c) - 1) / 2) # order of the Pade approximant
|
||||
|
||||
# Solve coefficients as a linear system
|
||||
M_left_element(c, i) = alt_sqrt(c - c0)^i
|
||||
M_left = M_left_element.(training_c, (0:order)')
|
||||
M_right = -training_k .* M_left[:, 2:end] # remove the first column
|
||||
M = hcat(M_left, M_right) # M = [M_left | M_right]
|
||||
sol = M \ training_k
|
||||
a = sol[1:order+1]
|
||||
b = [1; sol[order+2:end]]
|
||||
|
||||
# Pade approximant
|
||||
polynomial(a, c) = sum(i -> a[i+1] * alt_sqrt(c - c0)^i, 0:order)
|
||||
pade_approx(c) = polynomial(a, c) / polynomial(b, c)
|
||||
|
||||
# Extrapolate
|
||||
extrapolated_k = pade_approx.([training_c; extrapolating_c])
|
||||
extrapolated_E = (extrapolated_k .^ 2) / (2μ)
|
||||
|
||||
# Plotting
|
||||
scatter(real.(training_E), imag.(training_E), label="training")
|
||||
scatter!(real.(exact_E), imag.(exact_E), label="exact")
|
||||
scatter!(real.(extrapolated_E), imag.(extrapolated_E), label="extrapolated", m=:star5)
|
||||
|
|
@ -1,12 +1,12 @@
|
|||
include("../EC.jl")
|
||||
include("../common.jl")
|
||||
using Plots
|
||||
include("../helper.jl")
|
||||
include("../p_space.jl")
|
||||
|
||||
# contour
|
||||
p, w = get_mesh([0, 0.4 - 0.15im, 0.8, 6], [128, 128, 128])
|
||||
|
||||
μ = 0.5
|
||||
V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20)
|
||||
# ResonanceEC: Eq. (20)
|
||||
V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))
|
||||
|
||||
# generating a Berggren basis with a pole using the same system
|
||||
basis_c = 0.6
|
||||
|
|
@ -16,18 +16,48 @@ N_berg = sqrt.(diag(transpose(berg_basis .* w) * berg_basis))
|
|||
berg_basis = berg_basis ./ transpose(N_berg)
|
||||
berg_basis_w = berg_basis .* w
|
||||
|
||||
H0 = transpose(berg_basis_w) * get_T_matrix(p, μ) * berg_basis
|
||||
V = transpose(berg_basis_w) * get_V_matrix(V_system(1), p, w) * berg_basis
|
||||
training_points = range(1.1, 0.9, 5) # original: range(1.35, 0.9, 5)
|
||||
|
||||
training_c = range(1.1, 0.9, 5) # original: range(1.35, 0.9, 5)
|
||||
extrapolating_c = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
|
||||
training_E = Vector{ComplexF64}(undef, length(training_points))
|
||||
EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points))
|
||||
|
||||
training_ref = -0.26
|
||||
extrapolating_ref = [quick_pole_E(V_system(c)) for c in extrapolating_c]
|
||||
# training
|
||||
for (j, c) in enumerate(training_points)
|
||||
H = get_H_matrix(V_system(c), p, w)
|
||||
H_berg = transpose(berg_basis_w) * H * berg_basis
|
||||
evals, evecs = eigen(H_berg)
|
||||
i = argmin(real.(evals))
|
||||
# i = identify_pole_i(basis_p, evals)
|
||||
training_E[j] = evals[i]
|
||||
EC_basis[:, j] = evecs[:, i]
|
||||
end
|
||||
|
||||
EC = affine_EC(H0, V)
|
||||
train!(EC, training_c; ref_eval=training_ref, CAEC=true)
|
||||
extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref)
|
||||
EC_basis = hcat(EC_basis, conj.(EC_basis)) # CA-EC
|
||||
N_EC = transpose(EC_basis) * EC_basis
|
||||
|
||||
exportCSV(EC, "temp/2b_GSM_B2R.csv")
|
||||
plot(EC, "temp/2b_GSM_B2R.pdf"; basis_points=basis_E, xlims=(0, 0.3), ylims=(-0.120, 0.020))
|
||||
extrapolate_points = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
|
||||
|
||||
exact_E = Vector{ComplexF64}(undef, length(extrapolate_points))
|
||||
extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points))
|
||||
|
||||
# extrapolating
|
||||
for (j, c) in enumerate(extrapolate_points)
|
||||
exact_E[j] = quick_pole_E(V_system(c))
|
||||
|
||||
H = get_H_matrix(V_system(c), p, w)
|
||||
H_berg = transpose(berg_basis_w) * H * berg_basis
|
||||
H_EC = transpose(EC_basis) * H_berg * EC_basis
|
||||
evals = eigvals(H_EC, N_EC)
|
||||
i = argmin(abs.(evals .- exact_E[j]))
|
||||
extrapolate_E[j] = evals[i]
|
||||
end
|
||||
|
||||
exportCSV("temp/2b_GSM_B2R.csv", (training_E, exact_E, extrapolate_E), ("training", "exact", "extrapolated"))
|
||||
|
||||
scatter(real.(training_E), imag.(training_E), label="training")
|
||||
scatter!(real.(exact_E), imag.(exact_E), label="exact")
|
||||
scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated")
|
||||
scatter!(real.(basis_E), imag.(basis_E), m=:x, label="Berggren basis")
|
||||
xlims!(0,0.3)
|
||||
ylims!(-0.120,0.020)
|
||||
savefig("temp/2b_GSM_B2R.pdf")
|
||||
|
|
|
|||
|
|
@ -1,30 +1,51 @@
|
|||
include("../EC.jl")
|
||||
include("../common.jl")
|
||||
using Plots
|
||||
include("../p_space.jl")
|
||||
|
||||
berggren_mesh = get_mesh([0, 0.4 - 0.15im, 0.8, 6], [128, 128, 128])
|
||||
csm_mesh = get_mesh([0, 8 - 3im], [512])
|
||||
|
||||
for (mesh, name) in zip((berggren_mesh, csm_mesh), ("beggren", "csm"))
|
||||
for mesh in (berggren_mesh, csm_mesh)
|
||||
p, w = mesh
|
||||
mesh_E = p.*p ./ (2*0.5)
|
||||
|
||||
μ = 0.5
|
||||
V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20)
|
||||
|
||||
H0 = get_T_matrix(p, μ)
|
||||
V = get_V_matrix(V_system(1), p, w)
|
||||
|
||||
training_c = range(1.1, 0.9, 5) # original: range(1.35, 0.9, 5)
|
||||
extrapolating_c = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
|
||||
# ResonanceEC: Eq. (20)
|
||||
V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))
|
||||
|
||||
training_ref = [quick_pole_E(V_system(c)) for c in training_c]
|
||||
exact_E = [quick_pole_E(V_system(c)) for c in extrapolating_c]
|
||||
training_points = range(1.1, 0.9, 5) # original: range(1.35, 0.9, 5)
|
||||
training_E = Vector{ComplexF64}(undef, length(training_points))
|
||||
EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points))
|
||||
|
||||
EC = affine_EC(H0, V, w)
|
||||
train!(EC, training_c; ref_eval=training_ref, CAEC=true)
|
||||
extrapolate!(EC, extrapolating_c; precalculated_exact_E=exact_E)
|
||||
|
||||
#exportCSV(EC, "temp/2b_comparison_$name.csv")
|
||||
plot(EC, "temp/2b_comparison_$name.pdf"; basis_contour=mesh_E, xlims=(-0.3,0.3), ylims=(-0.120,0.020))
|
||||
for (j, c) in enumerate(training_points)
|
||||
evals, evecs = eigen(get_H_matrix(V_system(c), p, w))
|
||||
i = identify_pole_i(p, evals)
|
||||
training_E[j] = evals[i]
|
||||
EC_basis[:, j] = evecs[:, i]
|
||||
end
|
||||
|
||||
EC_basis = hcat(EC_basis, conj.(EC_basis)) # CA-EC
|
||||
EC_basis_w = EC_basis .* w
|
||||
N_EC = transpose(EC_basis_w) * EC_basis
|
||||
|
||||
extrapolate_points = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
|
||||
|
||||
exact_E = Vector{ComplexF64}(undef, length(extrapolate_points))
|
||||
extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points))
|
||||
|
||||
for (j, c) in enumerate(extrapolate_points)
|
||||
exact_E[j] = quick_pole_E(V_system(c))
|
||||
|
||||
H = get_H_matrix(V_system(c), p, w)
|
||||
H_EC = transpose(EC_basis_w) * H * EC_basis
|
||||
evals = eigvals(H_EC, N_EC)
|
||||
i = argmin(abs.(evals .- exact_E[j]))
|
||||
extrapolate_E[j] = evals[i]
|
||||
end
|
||||
|
||||
scatter(real.(training_E), imag.(training_E), label="training")
|
||||
scatter!(real.(exact_E), imag.(exact_E), label="exact")
|
||||
scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated")
|
||||
plot!(real.(mesh_E), imag.(mesh_E), label="contour")
|
||||
xlims!(-0.3,0.3)
|
||||
ylims!(-0.120,0.020)
|
||||
savefig("temp/" * string(rand(UInt16)) * ".pdf")
|
||||
end
|
||||
|
|
@ -1,81 +0,0 @@
|
|||
#%%
|
||||
import pandas as pd
|
||||
import torch
|
||||
import numpy as np
|
||||
|
||||
#%%
|
||||
df = pd.read_csv('../temp/2body_data.csv').sort_values(by='c')
|
||||
df['E'] = df['re_E'] + 1j * df['im_E']
|
||||
train_data = df[df['re_E'] < 0]
|
||||
target_data = df[df['re_E'] > 0]
|
||||
|
||||
train_cs = train_data['c'].to_numpy()
|
||||
train_Es = torch.tensor(train_data['E'].to_numpy(), dtype=torch.complex128)
|
||||
|
||||
#%%
|
||||
# hyperparameters
|
||||
N = 9
|
||||
|
||||
# initialize random Hamiltonians
|
||||
H0 = torch.randn(N, N, dtype=torch.complex128)
|
||||
H0 = (H0 + torch.transpose(H0, 0, 1)).requires_grad_() # symmetric
|
||||
H1 = torch.randn(N, N, dtype=torch.complex128)
|
||||
H1 = (H1 + torch.transpose(H1, 0, 1)).requires_grad_() # symmetric
|
||||
|
||||
#%%
|
||||
# training
|
||||
|
||||
# generate a set of c values to follow by subdividing the training cs
|
||||
subdivisions = 3
|
||||
c_steps = np.concatenate([np.linspace(start, stop, subdivisions, endpoint=False) for (start, stop) in zip(train_cs, train_cs[1:])])
|
||||
c_steps = np.append(c_steps, train_cs[-1])
|
||||
|
||||
lr = 0.05
|
||||
epochs = 100000
|
||||
for epoch in range(epochs):
|
||||
Es = torch.empty(len(train_data), dtype=torch.complex128)
|
||||
current_E = 0.0 # start at the threshold
|
||||
for c in c_steps:
|
||||
H = H0 + c * H1
|
||||
evals = torch.linalg.eigvals(H)
|
||||
current_E = evals[torch.argmin(torch.abs(evals - current_E))]
|
||||
if np.any(c == train_cs):
|
||||
index = np.where(c == train_cs)[0][0]
|
||||
Es[index] = current_E
|
||||
|
||||
loss = ((Es - train_Es).abs() ** 2).sum()
|
||||
|
||||
if epoch % 1000 == 0:
|
||||
print(f"Training {(epoch+1)/epochs:.1%} \t Loss: {loss}")
|
||||
|
||||
if H0.grad is not None:
|
||||
H0.grad.zero_()
|
||||
if H1.grad is not None:
|
||||
H1.grad.zero_()
|
||||
loss.backward()
|
||||
|
||||
with torch.no_grad():
|
||||
H0 -= lr * H0.grad
|
||||
H1 -= lr * H1.grad
|
||||
|
||||
# %%
|
||||
# evaluate for all points
|
||||
all_c = torch.tensor(df['c'].values, dtype=torch.float64)
|
||||
exact_E = torch.tensor(df['E'].values, dtype=torch.complex128)
|
||||
pred_Es = torch.empty(len(df), dtype=torch.complex128)
|
||||
with torch.no_grad():
|
||||
for (index, (c, E)) in enumerate(zip(all_c, exact_E)):
|
||||
H = H0 + c * H1
|
||||
evals = torch.linalg.eigvals(H)
|
||||
i = torch.argmin(torch.abs(evals - E)) # TODO: more robust way to identify the eigenvector
|
||||
pred_Es[index]= evals[i]
|
||||
|
||||
# %%
|
||||
# plot the results
|
||||
import matplotlib.pyplot as plt
|
||||
plt.scatter(train_data['re_E'], train_data['im_E'], label='training')
|
||||
plt.scatter(target_data['re_E'], target_data['im_E'], label='target')
|
||||
plt.scatter(pred_Es.real, pred_Es.imag, marker='x', label='predicted')
|
||||
plt.legend()
|
||||
|
||||
# %%
|
||||
|
|
@ -1,28 +0,0 @@
|
|||
include("../EC.jl")
|
||||
include("../common.jl")
|
||||
include("../p_space.jl")
|
||||
|
||||
# contour
|
||||
vertices = [0, 0.4 - 0.15im, 0.8, 6]
|
||||
subdivisions = [128, 128, 128]
|
||||
p, w = get_mesh(vertices, subdivisions)
|
||||
mesh_E = p.*p ./ (2*0.5)
|
||||
|
||||
μ = 0.5
|
||||
V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20)
|
||||
|
||||
H0 = get_T_matrix(p, μ)
|
||||
V = get_V_matrix(V_system(1), p, w)
|
||||
|
||||
training_c = range(0.75, 0.45, 5)
|
||||
extrapolating_c = range(0.40, 0.25, 5)
|
||||
|
||||
training_ref = [quick_pole_E(V_system(c)) for c in training_c]
|
||||
exact_E = [quick_pole_E(V_system(c)) for c in extrapolating_c]
|
||||
|
||||
EC = affine_EC(H0, V, w)
|
||||
train!(EC, training_c; ref_eval=training_ref, CAEC=false)
|
||||
extrapolate!(EC, extrapolating_c; precalculated_exact_E=exact_E)
|
||||
|
||||
#exportCSV(EC, "temp/2b_R2R.csv")
|
||||
plot(EC, "temp/2b_R2R.pdf"; basis_contour=mesh_E, xlims=(0, 1))
|
||||
|
|
@ -1,33 +1,66 @@
|
|||
include("../EC.jl")
|
||||
include("../common.jl")
|
||||
using Plots
|
||||
include("../helper.jl")
|
||||
include("../p_space.jl")
|
||||
|
||||
# contour
|
||||
p, w = get_mesh([0, 0.4 - 0.15im, 0.8, 6], [128, 128, 128])
|
||||
p, w = get_mesh([0, 0.4 - 0.08im, 0.8, 6], [128, 128, 128])
|
||||
contour_E = p.^2
|
||||
|
||||
μ = 0.5
|
||||
V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20)
|
||||
# ResonanceEC: Eq. (20)
|
||||
V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))
|
||||
|
||||
# generating a Berggren basis with a pole using the same system
|
||||
basis_c = 0.6
|
||||
basis_E, berg_basis = eigen(get_H_matrix(V_system(basis_c), p, w); permute=false, scale=false)
|
||||
pole_E = quick_pole_E(V_system(basis_c)) # basis only has 1 pole
|
||||
basis_p = sqrt.(basis_E)
|
||||
N_berg = sqrt.(diag(transpose(berg_basis .* w) * berg_basis))
|
||||
berg_basis = berg_basis ./ transpose(N_berg)
|
||||
berg_basis_w = berg_basis .* w
|
||||
|
||||
H0 = transpose(berg_basis_w) * get_T_matrix(p, μ) * berg_basis
|
||||
V = transpose(berg_basis_w) * get_V_matrix(V_system(1), p, w) * berg_basis
|
||||
training_points = range(0.79, 0.66, 4) # original: range(1.35, 0.9, 5)
|
||||
|
||||
training_c = range(0.79, 0.66, 4) # original: range(1.35, 0.9, 5)
|
||||
extrapolating_c = range(0.62, 0.40, 6) # original: range(0.75, 0.40, 8)
|
||||
training_E = Vector{ComplexF64}(undef, length(training_points))
|
||||
EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points))
|
||||
|
||||
training_ref = [quick_pole_E(V_system(c)) for c in training_c]
|
||||
extrapolating_ref = [quick_pole_E(V_system(c)) for c in extrapolating_c]
|
||||
# training
|
||||
for (j, c) in enumerate(training_points)
|
||||
H = get_H_matrix(V_system(c), p, w)
|
||||
H_berg = transpose(berg_basis_w) * H * berg_basis
|
||||
evals, evecs = eigen(H_berg)
|
||||
i = identify_pole_i(basis_p, evals)
|
||||
training_E[j] = evals[i]
|
||||
EC_basis[:, j] = evecs[:, i]
|
||||
end
|
||||
|
||||
EC = affine_EC(H0, V)
|
||||
train!(EC, training_c; ref_eval=training_ref, CAEC=false)
|
||||
extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref)
|
||||
N_EC = transpose(EC_basis) * EC_basis
|
||||
|
||||
exportCSV(EC, "temp/2b_GSM_R2R.csv")
|
||||
plot(EC, "temp/2b_GSM_R2R.pdf"; basis_points=basis_E, xlims=(0, 0.3), ylims=(-0.120, 0.020))
|
||||
extrapolate_points = range(0.62, 0.40, 6) # original: range(0.75, 0.40, 8)
|
||||
|
||||
exact_E = Vector{ComplexF64}(undef, length(extrapolate_points))
|
||||
extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points))
|
||||
|
||||
# extrapolating
|
||||
for (j, c) in enumerate(extrapolate_points)
|
||||
exact_E[j] = quick_pole_E(V_system(c))
|
||||
|
||||
H = get_H_matrix(V_system(c), p, w)
|
||||
H_berg = transpose(berg_basis_w) * H * berg_basis
|
||||
H_EC = transpose(EC_basis) * H_berg * EC_basis
|
||||
evals = eigvals(H_EC, N_EC)
|
||||
i = argmin(abs.(evals .- exact_E[j]))
|
||||
extrapolate_E[j] = evals[i]
|
||||
end
|
||||
|
||||
exportCSV("temp/2b_GSM_R2R.csv", (training_E, exact_E, extrapolate_E, [pole_E]), ("training", "exact", "extrapolated", "basis"))
|
||||
|
||||
contour_E_export = contour_E[real.(contour_E) .< 1] # to trim unnecessary data outside axis limits
|
||||
exportCSV("temp/2b_GSM_R2R_contour.csv", contour_E_export)
|
||||
|
||||
scatter(real.(training_E), imag.(training_E), label="training")
|
||||
scatter!(real.(exact_E), imag.(exact_E), label="exact")
|
||||
scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated")
|
||||
scatter!(real.(basis_E), imag.(basis_E), m=:x, label="Berggren basis")
|
||||
xlims!(0,0.3)
|
||||
ylims!(-0.120,0.020)
|
||||
savefig("temp/2b_GSM_R2R.pdf")
|
||||
|
|
|
|||
|
|
@ -1,39 +0,0 @@
|
|||
using Plots
|
||||
|
||||
include("../../EC.jl")
|
||||
include("../../ho_basis.jl")
|
||||
include("../../p_space.jl")
|
||||
|
||||
angle = 0.25 * pi # DOESN'T WORK WITHOUT ROTATION
|
||||
μω_gen = 0.5 * exp(-2im * angle)
|
||||
μ = 0.5
|
||||
l = 0
|
||||
V1 = -5
|
||||
R1 = sqrt(3)
|
||||
V2 = 2
|
||||
R2 = sqrt(10)
|
||||
n_max = 15
|
||||
|
||||
ns = collect(0:n_max)
|
||||
ls = fill(l, n_max + 1)
|
||||
|
||||
T = get_sp_T_matrix(ns, ls; μω_gen=μω_gen, μ=μ)
|
||||
V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns); μω_gen=μω_gen) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns); μω_gen=μω_gen)
|
||||
|
||||
n_EC = 8
|
||||
train_cs = (0.7 .+ 0.05 * randn(n_EC)) - 1im * (0.2 .+ 0.05 * randn(n_EC))
|
||||
target_cs = [0.5]
|
||||
|
||||
near_E = 0.2 + 0.2im
|
||||
exact_E = [0.20845136860234303 - 0.07100640993695649im]
|
||||
|
||||
EC = affine_EC(T, V; ensemble_size=32)
|
||||
train!(EC, train_cs; ref_eval=near_E, CAEC=false)
|
||||
extrapolate!(EC, target_cs; precalculated_exact_E=exact_E)
|
||||
|
||||
plot(EC; xlims=(0,0.3), ylims=(-0.3,0.3))
|
||||
hline!([0], color=:red, label="continuum")
|
||||
xlabel!("Re(E)")
|
||||
ylabel!("Im(E)")
|
||||
plot!(legend=:bottomleft)
|
||||
savefig("temp/2b_HO_XZ.pdf")
|
||||
|
|
@ -1,50 +0,0 @@
|
|||
using Plots
|
||||
|
||||
include("../../EC.jl")
|
||||
include("../../ho_basis.jl")
|
||||
include("../../p_space.jl")
|
||||
|
||||
# paramters of the system
|
||||
|
||||
angle = 0.0
|
||||
μ = 0.5
|
||||
l = 0
|
||||
V1 = -5
|
||||
R1 = sqrt(3)
|
||||
V2 = 2
|
||||
R2 = sqrt(10)
|
||||
|
||||
n_EC = 8
|
||||
train_cs = (0.7 .+ 0.03 * randn(n_EC)) - 1im * (0.2 .+ 0.03 * randn(n_EC))
|
||||
near_E = 0.2 + 0.2im
|
||||
|
||||
target_c = 0.5
|
||||
exact_E = 0.20845136860234303 - 0.07100640993695649im
|
||||
|
||||
vertices = [0, 4 * exp(-1im * angle)]
|
||||
subdivisions = [256]
|
||||
ks, ws = get_mesh(vertices, subdivisions)
|
||||
|
||||
V_of_r(r) = V1 * exp(-r^2 / R1^2) + V2 * exp(-r^2 / R2^2)
|
||||
V_mat_elem(k, kp) = Vl_mat_elem(V_of_r, l, k, kp; atol=10^-5, maxevals=10^5, R_cutoff=16)
|
||||
V = get_V_matrix(V_mat_elem, ks, ws)
|
||||
T = get_T_matrix(ks, μ)
|
||||
|
||||
EC_p_space = affine_EC(T, V)
|
||||
train!(EC_p_space, train_cs; ref_eval=near_E, CAEC=false)
|
||||
extrapolate!(EC_p_space, [target_c]; precalculated_exact_E=[exact_E])
|
||||
|
||||
# Plotting
|
||||
|
||||
theme(:dark) # Set the global theme to dark
|
||||
|
||||
scatter([real(exact_E)], [imag(exact_E)], label="exact", marker=:circle, markercolor=:white, bg = :black) # black background
|
||||
scatter!(real.(EC_p_space.training_E), imag.(EC_p_space.training_E), label="training", marker=:circle, color=:blue)
|
||||
scatter!(real.(EC_p_space.extrapolated_E), imag.(EC_p_space.extrapolated_E), label="extrapolated", marker=:x, color=:green)
|
||||
hline!([0], color=:red, label="continuum")
|
||||
plot!(legend=:bottomleft)
|
||||
xlabel!("Re(E)")
|
||||
ylabel!("Im(E)")
|
||||
xlims!(0, 0.3)
|
||||
ylims!(-0.3, 0.3)
|
||||
savefig("temp/2body_p_space.pdf")
|
||||
|
|
@ -1,36 +0,0 @@
|
|||
include("../../ho_basis.jl")
|
||||
include("../../EC.jl")
|
||||
|
||||
V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
|
||||
ϕ = 0.1 # DOESN'T WORK WITHOUT ROTATION
|
||||
μω_global = 0.5 * exp(-2im * ϕ)
|
||||
E_max = 40
|
||||
|
||||
H0 = get_3b_H_matrix(jacobi, V_of_r, μω_global, E_max, Λ, m, true, true)
|
||||
|
||||
# Vp = perturbation to make the state artificially bound
|
||||
Vp_of_r(r) = -exp(-(r/3)^2)
|
||||
@time "Vp" Vp = get_3b_H_matrix(jacobi, Vp_of_r, μω_global, E_max, Λ, m, false, true)
|
||||
|
||||
training_ref = 2 + 0.5im
|
||||
|
||||
exact_E = [4.076642792419057-0.012998408352259658im,
|
||||
3.6129849325287-0.007397677539402868im,
|
||||
3.145212908643357-0.0038660337822150753im,
|
||||
2.6729225739451596-0.0021090370393881063im,
|
||||
2.196385760253282-0.0010430088245526555im,
|
||||
1.7162659936896967-0.0004515351140200029,
|
||||
1.2329926791785895-0.00017698044022813525im]
|
||||
|
||||
training_c = [0.6 - 0.16im] .+ 0.04 .* (randn(8) .+ 0.5im * randn(8))
|
||||
extrapolating_c = 0.0 : 0.2 : 1.2
|
||||
|
||||
EC = affine_EC(H0, Vp)
|
||||
train!(EC, training_c; ref_eval=training_ref, CAEC=false)
|
||||
extrapolate!(EC, extrapolating_c; precalculated_exact_E=exact_E)
|
||||
|
||||
exportCSV(EC, "temp/3b_HO_XZ.csv")
|
||||
plot(EC, "temp/3b_HO_XZ.pdf")
|
||||
|
|
@ -1,45 +0,0 @@
|
|||
include("../../p_space.jl")
|
||||
include("../../EC.jl")
|
||||
|
||||
using Arpack
|
||||
|
||||
# target = 4.0766890719636875 - 0.012758927741074495im
|
||||
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
|
||||
|
||||
vertices = [0, 2 - 0.2im, 3, 4] # TODO: real contour instead of Berggren basis
|
||||
subdivisions = [15, 10, 10]
|
||||
jmax = 4
|
||||
|
||||
E_max = 40
|
||||
μω_global = 0.5
|
||||
|
||||
@time "H0" H0, _ = get_3b_H_matrix(jacobi, V_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m, true, true)
|
||||
|
||||
# Vp = perturbation to make the state artificially bound
|
||||
Vp_of_r(r) = -exp(-(r/3)^2)
|
||||
@time "Vp" Vp, _ = get_3b_H_matrix(jacobi, Vp_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m, false, true)
|
||||
|
||||
training_ref = 2 + 0.5im
|
||||
|
||||
exact_E = [4.076642792419057-0.012998408352259658im,
|
||||
3.6129849325287-0.007397677539402868im,
|
||||
3.145212908643357-0.0038660337822150753im,
|
||||
2.6729225739451596-0.0021090370393881063im,
|
||||
2.196385760253282-0.0010430088245526555im,
|
||||
1.7162659936896967-0.0004515351140200029,
|
||||
1.2329926791785895-0.00017698044022813525im]
|
||||
|
||||
training_c = [-1.5 - 0.5im] .+ (randn(8) .+ 0.05im * randn(8))
|
||||
extrapolating_c = 0.0 : 0.2 : 1.2
|
||||
|
||||
EC = affine_EC(H0, Vp)
|
||||
train!(EC, training_c; ref_eval=training_ref, CAEC=false)
|
||||
extrapolate!(EC, extrapolating_c; precalculated_exact_E=exact_E)
|
||||
|
||||
exportCSV(EC, "temp/3b_p_space_XZ.csv")
|
||||
plot(EC, "temp/3b_p_space_XZ.pdf")
|
||||
|
||||
# Results: training points are all over the place, and extrapolated values are garbage.
|
||||
|
|
@ -1,71 +0,0 @@
|
|||
using Plots
|
||||
|
||||
include("../../EC.jl")
|
||||
include("../../ho_basis.jl")
|
||||
include("../../p_space.jl")
|
||||
|
||||
# paramters of the system
|
||||
|
||||
angle = 0.0 * pi
|
||||
μ = 0.5
|
||||
l = 0
|
||||
V1 = -5
|
||||
R1 = sqrt(3)
|
||||
V2 = 2
|
||||
R2 = sqrt(10)
|
||||
|
||||
n_EC = 8
|
||||
train_cs = (0.7 .+ 0.03 * randn(n_EC)) - 1im * (0.2 .+ 0.03 * randn(n_EC))
|
||||
near_E = 0.2 + 0.2im
|
||||
|
||||
target_c = 0.5
|
||||
exact_E = 0.20845136860234303 - 0.07100640993695649im
|
||||
|
||||
# HO basis
|
||||
|
||||
global EC_HO
|
||||
begin
|
||||
println("HO basis calculation")
|
||||
μω_gen = 0.5 * exp(-1im * angle)
|
||||
n_max = 40
|
||||
ns = collect(0:n_max)
|
||||
ls = fill(l, n_max + 1)
|
||||
|
||||
T = get_sp_T_matrix(ns, ls; μω_gen=μω_gen, μ=μ)
|
||||
V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns); μω_gen=μω_gen) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns); μω_gen=μω_gen)
|
||||
|
||||
global EC_HO = affine_EC(T, V)
|
||||
train!(EC_HO, train_cs; ref_eval=near_E, CAEC=false)
|
||||
extrapolate!(EC_HO, [target_c]; precalculated_exact_E=[exact_E])
|
||||
end
|
||||
|
||||
# p-space
|
||||
|
||||
global EC_p_space
|
||||
begin
|
||||
println("p-space calculation")
|
||||
vertices = [0, 4 * exp(-1im * angle)]
|
||||
subdivisions = [256]
|
||||
ks, ws = get_mesh(vertices, subdivisions)
|
||||
|
||||
V_of_r(r) = V1 * exp(-r^2 / R1^2) + V2 * exp(-r^2 / R2^2)
|
||||
V_mat_elem(k, kp) = Vl_mat_elem(V_of_r, l, k, kp; atol=10^-5, maxevals=10^5, R_cutoff=16)
|
||||
V = get_V_matrix(V_mat_elem, ks, ws)
|
||||
T = get_T_matrix(ks, μ)
|
||||
|
||||
global EC_p_space = affine_EC(T, V)
|
||||
train!(EC_p_space, train_cs; ref_eval=near_E, CAEC=false)
|
||||
extrapolate!(EC_p_space, [target_c]; precalculated_exact_E=[exact_E])
|
||||
end
|
||||
|
||||
# Plotting
|
||||
|
||||
scatter([real(exact_E)], [imag(exact_E)], label="Exact", marker=:circle, markercolor=:white)
|
||||
scatter!(real.(EC_HO.training_E), imag.(EC_HO.training_E), label="HO basis training", marker=:circle, color=:blue)
|
||||
scatter!(real.(EC_HO.extrapolated_E), imag.(EC_HO.extrapolated_E), label="HO basis extrapolated", marker=:x, color=:blue)
|
||||
scatter!(real.(EC_p_space.training_E), imag.(EC_p_space.training_E), label="p-space training", marker=:circle, color=:red)
|
||||
scatter!(real.(EC_p_space.extrapolated_E), imag.(EC_p_space.extrapolated_E), label="p-space extrapolated", marker=:x, color=:red)
|
||||
plot!(legend=:bottomleft)
|
||||
xlabel!("Re(E)")
|
||||
ylabel!("Im(E)")
|
||||
savefig("temp/2b_p_space_vs_HO.pdf")
|
||||
|
|
@ -0,0 +1,55 @@
|
|||
using LinearAlgebra, Plots
|
||||
include("../ho_basis.jl")
|
||||
include("../p_space.jl")
|
||||
|
||||
μω_gen = 0.5 * exp(-1im * 0.47 * pi)
|
||||
μ = 0.5
|
||||
l = 0
|
||||
V1 = -5
|
||||
R1 = sqrt(3)
|
||||
V2 = 2
|
||||
R2 = sqrt(10)
|
||||
n_max = 15
|
||||
|
||||
ns = collect(0:n_max)
|
||||
ls = fill(l, n_max + 1)
|
||||
|
||||
T = get_sp_T_matrix(ns, ls; μω_gen=μω_gen, μ=μ)
|
||||
V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns); μω_gen=μω_gen) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns); μω_gen=μω_gen)
|
||||
|
||||
n_EC = 8
|
||||
train_cs = (0.7 .+ 0.05 * randn(n_EC)) - 1im * (0.2 .+ 0.05 * randn(n_EC))
|
||||
target_cs = range(0.77, 0.22, 6)
|
||||
|
||||
train_E = zeros(ComplexF64, n_EC)
|
||||
EC_basis = zeros(ComplexF64, (n_max + 1, length(train_cs)))
|
||||
exact_E = zeros(ComplexF64, length(target_cs))
|
||||
extrapolate_E = similar(exact_E)
|
||||
|
||||
near_E = 0.2 + 0.2im
|
||||
|
||||
for (j, c) in enumerate(train_cs)
|
||||
H = T + c .* V
|
||||
evals, evecs = eigen(H)
|
||||
i = argmin(abs.(evals .- near_E))
|
||||
train_E[j] = evals[i]
|
||||
EC_basis[:, j] = evecs[:, i]
|
||||
end
|
||||
|
||||
N_EC = transpose(EC_basis) * EC_basis
|
||||
|
||||
for (j, c) in enumerate(target_cs)
|
||||
exact_E[j] = quick_pole_E((p, q) -> c*(V1*g0(R1, p, q) + V2*g0(R2, p, q)), μ; cs_angle=0.5)
|
||||
|
||||
H = T + c .* V
|
||||
H_EC = transpose(EC_basis) * H * EC_basis
|
||||
evals = eigvals(H_EC, N_EC)
|
||||
i = argmin(abs.(evals .- exact_E[j]))
|
||||
extrapolate_E[j] = evals[i]
|
||||
end
|
||||
|
||||
scatter(real.(train_E), imag.(train_E), label="training")
|
||||
scatter!(real.(exact_E), imag.(exact_E), label="exact")
|
||||
scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated")
|
||||
xlims!(-0.2,0.3)
|
||||
ylims!(-0.3,0.3)
|
||||
106
common.jl
106
common.jl
|
|
@ -1,106 +0,0 @@
|
|||
using LinearAlgebra, DelimitedFiles, SparseArrays
|
||||
|
||||
@enum coordinate_system jacobi src
|
||||
|
||||
"Square root function with the branch cut along the postive imaginary axis"
|
||||
alt_sqrt(x::Number)::ComplexF64 = sqrt(im * x) / sqrt(im)
|
||||
|
||||
"Sum over array while minimizing catastrophic cancellation as much as possible"
|
||||
function better_sum(arr::Array{T}) where T<:Real
|
||||
pos_arr = arr[arr .> 0]
|
||||
neg_arr = arr[arr .< 0]
|
||||
|
||||
sort!(pos_arr)
|
||||
sort!(neg_arr, rev=true)
|
||||
|
||||
return sum(pos_arr) + sum(neg_arr)
|
||||
end
|
||||
|
||||
better_sum(arr::Array{ComplexF64}) = better_sum(real.(arr)) + 1im * better_sum(imag.(arr))
|
||||
|
||||
"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[nearestIndex(list, ref)]
|
||||
|
||||
"Simple implementation of the Kronecker sum"
|
||||
function kron_sum(A::AbstractMatrix, B::AbstractMatrix)
|
||||
@assert size(A, 1) == size(A, 2) && size(B, 1) == size(B, 2) "Matrices should be square"
|
||||
return kron(A, I(size(B, 1))) + kron(I(size(A, 1)), B)
|
||||
end
|
||||
|
||||
"Flattened vector version of Iterators.product(...) with index hierachy reversed -- leftmost index has the highest hierachy"
|
||||
iter_prod(args...) = reverse.(collect(Iterators.product(reverse(args)...))[:])
|
||||
|
||||
"Export CSV data for a scatter plot taking in data as a list of complex vectors (x=Re, y=Im), and a list of corresponding labels, typically used for EC results"
|
||||
function exportCSV(file::String, data, labels=nothing)
|
||||
columns = ["x" "y" "label"]
|
||||
|
||||
if isnothing(labels)
|
||||
columns[end] = ""
|
||||
labels = fill("", length(data))
|
||||
end
|
||||
|
||||
open(file, "w") do f
|
||||
writedlm(f, columns)
|
||||
for (d, label) in zip(data, labels)
|
||||
l = fill(label, length(d))
|
||||
writedlm(f, hcat(reim(d)..., l))
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
"In-place c-orthogonalization via (modified) Gram-Schmidt. Only significant vectors are returned (c-normalized).
|
||||
The number of significant vectors to return are determined by the original singular values (compared to the threshold)."
|
||||
function gram_schmidt!(vecs::Vector{Vector{T}}, ws=ones(length(vecs[1])), threshold=1e-5; verbose=false) where T<: Number
|
||||
c_product(i, j) = sum(vecs[i] .* ws .* vecs[j])
|
||||
norm(i) = c_product(i, i)
|
||||
proj(i, j) = (c_product(i, j) / norm(i)) .* vecs[i] # component of vec[i] in vec[j]
|
||||
|
||||
# initial normalization
|
||||
for i in eachindex(vecs)
|
||||
vecs[i] ./= sqrt(norm(i))
|
||||
end
|
||||
|
||||
verbose && println("Absolute singular values = $(round.(c_singular_values(vecs, ws); sigdigits=1))")
|
||||
target_dim = c_rank(vecs, ws, threshold)
|
||||
verbose && println("Target dimensionality = $target_dim")
|
||||
|
||||
selected_vecs_i = Integer[]
|
||||
|
||||
while length(selected_vecs_i) < target_dim
|
||||
i = argmax(i -> abs(norm(i)), setdiff(eachindex(vecs), selected_vecs_i)) # find the largest vector from the remaining
|
||||
push!(selected_vecs_i, i)
|
||||
for j in setdiff(eachindex(vecs), selected_vecs_i)
|
||||
vecs[j] .-= proj(i, j)
|
||||
end
|
||||
end
|
||||
|
||||
verbose && println("Absolute norms of selected vectors = $(round.(selected_vecs_i .|> norm .|> abs; sigdigits=1))")
|
||||
verbose && println("Absolute norms of dropped vectors = $(round.(setdiff(eachindex(vecs), selected_vecs_i) .|> norm .|> abs; sigdigits=1))")
|
||||
|
||||
# final normalization
|
||||
for i in selected_vecs_i
|
||||
vecs[i] ./= sqrt(norm(i))
|
||||
end
|
||||
|
||||
return vecs[selected_vecs_i]
|
||||
end
|
||||
|
||||
"Same as above but the basis is provided as a matrix instead of an array of vectors"
|
||||
gram_schmidt!(vecs::Matrix{T}, ws=ones(size(vecs, 1)), threshold=1e-5; verbose=false) where T<: Number = hcat(gram_schmidt!([vecs[:, i] for i in axes(vecs, 2)], ws, threshold; verbose=verbose)...)
|
||||
|
||||
"Rank of a basis set (pre-normalized) under c-product"
|
||||
c_rank(vecs, ws, threshold=1e-8) = count(c_singular_values(vecs, ws) .> threshold)
|
||||
|
||||
"Singular values (magnitudes) of a basis set (pre-normalized) under c-product"
|
||||
function c_singular_values(vecs, ws)
|
||||
basis = hcat(vecs...)
|
||||
N = transpose(basis) * spdiagm(ws) * basis
|
||||
singular_values = eigvals(N) .|> abs .|> sqrt
|
||||
return singular_values
|
||||
end
|
||||
|
|
@ -0,0 +1,50 @@
|
|||
using LinearAlgebra, DelimitedFiles
|
||||
|
||||
"Sum over array while minimizing catastrophic cancellation as much as possible"
|
||||
function better_sum(arr::Array{T}) where T<:Real
|
||||
pos_arr = arr[arr .> 0]
|
||||
neg_arr = arr[arr .< 0]
|
||||
|
||||
sort!(pos_arr)
|
||||
sort!(neg_arr, rev=true)
|
||||
|
||||
return sum(pos_arr) + sum(neg_arr)
|
||||
end
|
||||
|
||||
better_sum(arr::Array{ComplexF64}) = better_sum(real.(arr)) + 1im * better_sum(imag.(arr))
|
||||
|
||||
"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[nearestIndex(list, ref)]
|
||||
|
||||
"Simple implementation of the Kronecker sum"
|
||||
function kron_sum(A::AbstractMatrix, B::AbstractMatrix)
|
||||
@assert size(A, 1) == size(A, 2) && size(B, 1) == size(B, 2) "Matrices should be square"
|
||||
return kron(A, I(size(B, 1))) + kron(I(size(A, 1)), B)
|
||||
end
|
||||
|
||||
"Flattened vector version of Iterators.product(...) with index hierachy reversed -- leftmost index has the highest hierachy"
|
||||
iter_prod(args...) = reverse.(collect(Iterators.product(reverse(args)...))[:])
|
||||
|
||||
"Export CSV data for a scatter plot taking in data as a list of complex vectors (x=Re, y=Im), and a list of corresponding labels, typically used for EC results"
|
||||
function exportCSV(file::String, data, labels=nothing)
|
||||
columns = ["x" "y" "label"]
|
||||
|
||||
if isnothing(labels)
|
||||
columns[end] = ""
|
||||
labels = fill("", length(data))
|
||||
end
|
||||
|
||||
open(file, "w") do f
|
||||
writedlm(f, columns)
|
||||
for (d, label) in zip(data, labels)
|
||||
l = fill(label, length(d))
|
||||
writedlm(f, hcat(reim(d)..., l))
|
||||
end
|
||||
end
|
||||
end
|
||||
256
ho_basis.jl
256
ho_basis.jl
|
|
@ -1,60 +1,9 @@
|
|||
using SparseArrays
|
||||
using QuadGK
|
||||
using LRUCache
|
||||
include("common.jl")
|
||||
include("helper.jl")
|
||||
include("math.jl")
|
||||
|
||||
"2-body HO basis (used for 3-body systems without the CM DOF)"
|
||||
struct ho_basis_2B
|
||||
E_max::Int
|
||||
Λ::Int
|
||||
dim::Int # dimensionality of the basis
|
||||
Es::Vector{Int}
|
||||
n1s::Vector{Int}
|
||||
l1s::Vector{Int}
|
||||
n2s::Vector{Int}
|
||||
l2s::Vector{Int}
|
||||
|
||||
function ho_basis_2B(E_max, Λ=-1)
|
||||
Es = Int[]
|
||||
n1s = Int[]
|
||||
l1s = Int[]
|
||||
n2s = Int[]
|
||||
l2s = Int[]
|
||||
|
||||
# E = 2*n1 + l1 + 2*n2 + l2
|
||||
for E in E_max : -2 : 0 # same parity states only
|
||||
for n1 in 0 : E ÷ 2
|
||||
for n2 in 0 : (E - 2*n1) ÷ 2
|
||||
for l1 in 0 : (E - 2*n1 - 2*n2)
|
||||
l2 = E - 2*n1 - 2*n2 - l1
|
||||
if Λ≥0 && !triangle_ineq(l1, l2, Λ); continue; end
|
||||
push!(Es, E)
|
||||
push!(n1s, n1)
|
||||
push!(l1s, l1)
|
||||
push!(n2s, n2)
|
||||
push!(l2s, l2)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
return new(E_max, Λ, length(Es), Es, n1s, l1s, n2s, l2s)
|
||||
end
|
||||
end
|
||||
|
||||
"Number of possible distinct matrix elements for a given l."
|
||||
function cache_size_l(E_max::Int, l::Int)::Int
|
||||
n_max = (E_max - l) ÷ 2
|
||||
return (n_max * (n_max + 1)) ÷ 2
|
||||
end
|
||||
|
||||
"Number of possible distinct matrix elements for all l."
|
||||
cache_size(E_max::Int)::Int = sum(l -> cache_size_l(E_max, l), 0 : E_max)
|
||||
|
||||
"Preallocation of cache for PE matrix elements."
|
||||
prealloc_V_cache(E_max::Int, dtype::DataType=Float64) = LRU{Tuple{UInt8, UInt8, UInt8}, dtype}(maxsize=cache_size(E_max))
|
||||
|
||||
function V_numerical(V_of_r, l, n1, n2; μω_gen=1.0, atol=0, maxevals=10^7)
|
||||
const_part = sqrt(μω_gen) * ho_basis_const(l, n1) * ho_basis_const(l, n2)
|
||||
integrand(r) = ho_basis_func(l, n1, sqrt(μω_gen) * r) * ho_basis_func(l, n2, sqrt(μω_gen) * r) * V_of_r(r)
|
||||
|
|
@ -62,13 +11,56 @@ function V_numerical(V_of_r, l, n1, n2; μω_gen=1.0, atol=0, maxevals=10^7)
|
|||
return const_part * integral
|
||||
end
|
||||
|
||||
"KE matrix for a single DOF. Set kron_deltas=[other quantum numbers] for other DOFs which the operator does not act on.
|
||||
E.g. get_sp_T_matrix(n1s, l1s, kron_deltas=[n2s l2s])"
|
||||
function get_sp_T_matrix(ns, ls, kron_deltas=[]; μω_gen=1.0, μ=1.0)
|
||||
function get_sp_basis(E_max)
|
||||
Es = Int[]
|
||||
ns = Int[]
|
||||
ls = Int[]
|
||||
|
||||
# E = 2*n + l
|
||||
for E in 0 : E_max
|
||||
for n in 0 : E ÷ 2
|
||||
l = E - 2*n
|
||||
push!(Es, E)
|
||||
push!(ns, n)
|
||||
push!(ls, l)
|
||||
end
|
||||
end
|
||||
|
||||
return (Es, ns, ls)
|
||||
end
|
||||
|
||||
function get_2p_basis(E_max, Λ=-1)
|
||||
Es = Int[]
|
||||
n1s = Int[]
|
||||
l1s = Int[]
|
||||
n2s = Int[]
|
||||
l2s = Int[]
|
||||
|
||||
# E = 2*n1 + l1 + 2*n2 + l2
|
||||
for E in E_max : -2 : 0 # same parity states only
|
||||
for n1 in 0 : E ÷ 2
|
||||
for n2 in 0 : (E - 2*n1) ÷ 2
|
||||
for l1 in 0 : (E - 2*n1 - 2*n2)
|
||||
l2 = E - 2*n1 - 2*n2 - l1
|
||||
if Λ≥0 && !triangle_ineq(l1, l2, Λ); continue; end
|
||||
push!(Es, E)
|
||||
push!(n1s, n1)
|
||||
push!(l1s, l1)
|
||||
push!(n2s, n2)
|
||||
push!(l2s, l2)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
return (Es, n1s, l1s, n2s, l2s)
|
||||
end
|
||||
|
||||
function get_sp_T_matrix(ns, ls; mask=trues(length(ns),length(ns)), μω_gen=1.0, μ=1.0)
|
||||
mat = spzeros(length(ns), length(ns))
|
||||
for idx in CartesianIndices(mat)
|
||||
if !mask[idx]; continue; end
|
||||
(i, j) = Tuple(idx)
|
||||
all(arr -> arr[i]==arr[j], kron_deltas) || continue # check if all Kronecker deltas are non-zero
|
||||
if ls[i] == ls[j]
|
||||
if ns[i] == ns[j]
|
||||
mat[idx] = ns[j] + ls[i]/2 + 3/4
|
||||
|
|
@ -81,14 +73,11 @@ function get_sp_T_matrix(ns, ls, kron_deltas=[]; μω_gen=1.0, μ=1.0)
|
|||
return (μω_gen / μ) .* mat
|
||||
end
|
||||
|
||||
"PE matrix for a single DOF. Set kron_deltas=[other quantum numbers] for other DOFs which the operator does not act on.
|
||||
E.g. get_sp_V_matrix(n1s, l1s, kron_deltas=[n2s l2s])
|
||||
Providing a preallocated cache is optional. Otherwise, provided E_max value is used to initialize one (defaults to 100)."
|
||||
function get_sp_V_matrix(V_l, ns, ls, kron_deltas=[]; dtype=Float64, E_max=100, cache=prealloc_V_cache(E_max, dtype))
|
||||
function get_sp_V_matrix(V_l, ns, ls; mask=trues(length(ns),length(ns)), dtype=Float64, cache=LRU{Tuple{UInt8, UInt8, UInt8}, dtype}(maxsize=(1+maximum(ns))^2))
|
||||
mat = zeros(dtype, length(ns), length(ns))
|
||||
Threads.@threads for idx in CartesianIndices(mat)
|
||||
if !mask[idx]; continue; end
|
||||
(i, j) = Tuple(idx)
|
||||
all(arr -> arr[i]==arr[j], kron_deltas) || continue # check if all Kronecker deltas are non-zero
|
||||
if ls[i] == ls[j]
|
||||
l = UInt8(ls[i])
|
||||
n1, n2 = UInt8.(minmax(ns[i], ns[j])) # assuming transpose symmetry
|
||||
|
|
@ -98,12 +87,12 @@ function get_sp_V_matrix(V_l, ns, ls, kron_deltas=[]; dtype=Float64, E_max=100,
|
|||
return sparse(mat)
|
||||
end
|
||||
|
||||
function Moshinsky_transform(basis::ho_basis_2B)
|
||||
NQMAX = maximum(basis.Es)
|
||||
@assert all(mod.(basis.Es, 2) .== mod(NQMAX, 2)) "Can only admit basis states with same parity"
|
||||
function Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
|
||||
NQMAX = maximum(Es)
|
||||
@assert all(mod.(Es, 2) .== mod(NQMAX, 2)) "Can only admit basis states with same parity"
|
||||
|
||||
LMIN = basis.Λ
|
||||
LMAX = basis.Λ
|
||||
LMIN = Λ
|
||||
LMAX = Λ
|
||||
CO = 1/sqrt(2)
|
||||
SI = 1/sqrt(2)
|
||||
|
||||
|
|
@ -111,15 +100,15 @@ function Moshinsky_transform(basis::ho_basis_2B)
|
|||
BRAC = zeros(Float64, 1 + LMAX, 1 + (NQMAX - LMIN) ÷ 2, 1 + (NQMAX - LMIN) ÷ 2, 1 + (NQMAX - LMIN) ÷ 2, 1 + (NQMAX - LMIN) ÷ 2, 1 + LMAX, 1 + (NQMAX-LMIN) ÷ 2, 1 + LMAX-LMIN)
|
||||
@ccall "../OSBRACKETS/allosbrac.so".allosbrac_(NQMAX::Ref{Int32},LMIN::Ref{Int32},LMAX::Ref{Int32},CO::Ref{Float64},SI::Ref{Float64},BRAC::Ptr{Array{Float64}})::Cvoid
|
||||
|
||||
mat = zeros(basis.dim, basis.dim)
|
||||
mat = zeros(length(Es), length(Es))
|
||||
|
||||
s = hcat(basis.Es, basis.n1s, basis.l1s, basis.n2s, basis.l2s)
|
||||
s = hcat(Es, n1s, l1s, n2s, l2s)
|
||||
Threads.@threads for idx in CartesianIndices(mat)
|
||||
(i, j) = Tuple(idx)
|
||||
(Elhs, N, L, n, l) = s[i, :]
|
||||
(Erhs, n1, l1, n2, l2) = s[j, :]
|
||||
if Elhs == Erhs && triangle_ineq(L, l, basis.Λ) && triangle_ineq(l1, l2, basis.Λ)
|
||||
mat[i, j] = (-1)^(n1 + n2 + N + n) * pick_Moshinsky_bracket(BRAC, n1, l1, n2, l2, N, L, n, l, basis.Λ)
|
||||
if Elhs == Erhs && triangle_ineq(L, l, Λ) && triangle_ineq(l1, l2, Λ)
|
||||
mat[i, j] = (-1)^(n1 + n2 + N + n) * pick_Moshinsky_bracket(BRAC, n1, l1, n2, l2, N, L, n, l, Λ)
|
||||
end
|
||||
end
|
||||
|
||||
|
|
@ -137,115 +126,82 @@ function pick_Moshinsky_bracket(BRAC, n1′, l1′, n2′, l2′, n1, l1, n2, l2
|
|||
return BRAC[1 + NP, 1 + n1′, 1 + MP, 1 + n1, 1 + n2, 1 + N, 1 + M, 1]
|
||||
end
|
||||
|
||||
function get_jacobi_V_matrix(V_of_r, basis::ho_basis_2B, μ1ω1, μω_global; atol=10^-6, maxevals=10^5)
|
||||
V1 = get_jacobi_V1_matrix(V_of_r, basis, μ1ω1; atol=atol, maxevals=maxevals)
|
||||
V2 = get_jacobi_V2_matrix(V_of_r, basis, μω_global; atol=atol, maxevals=maxevals)
|
||||
function get_jacobi_V_matrix(V_of_r, E_max, Λ, μ1ω1, μω_global; atol=10^-6, maxevals=10^5)
|
||||
V1 = get_jacobi_V1_matrix(V_of_r, E_max, Λ, μ1ω1; atol=atol, maxevals=maxevals)
|
||||
V2 = get_jacobi_V2_matrix(V_of_r, E_max, Λ, μω_global; atol=atol, maxevals=maxevals)
|
||||
return V1 + V2
|
||||
end
|
||||
|
||||
function get_jacobi_V1_matrix(V_of_r, basis::ho_basis_2B, μ1ω1; atol=10^-6, maxevals=10^5)
|
||||
function get_jacobi_V1_matrix(V_of_r, E_max, Λ, μ1ω1; atol=10^-6, maxevals=10^5)
|
||||
_, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
|
||||
l_max = max(maximum(l1s), maximum(l2s))
|
||||
n_max = max(maximum(n1s), maximum(n2s))
|
||||
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
|
||||
|
||||
V1_elem(l, n1, n2) = V_numerical(V_of_r, l, n1, n2; μω_gen=μ1ω1, atol=atol, maxevals=maxevals)
|
||||
V1 = get_sp_V_matrix(V1_elem, basis.n1s, basis.l1s, [basis.n2s, basis.l2s]; dtype=ComplexF64, E_max=basis.E_max)
|
||||
V1_cache = LRU{Tuple{UInt8, UInt8, UInt8}, ComplexF64}(maxsize=(1+l_max)*(1+n_max)^2)
|
||||
V1 = get_sp_V_matrix(V1_elem, n1s, l1s; mask=mask1, dtype=ComplexF64, cache=V1_cache)
|
||||
|
||||
return V1
|
||||
end
|
||||
|
||||
function get_jacobi_V2_matrix(V_of_r, basis::ho_basis_2B, μω_global; atol=10^-6, maxevals=10^5)
|
||||
function get_jacobi_V2_matrix(V_of_r, E_max, Λ, μω_global; atol=10^-6, maxevals=10^5)
|
||||
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
|
||||
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 = prealloc_V_cache(basis.E_max, ComplexF64)
|
||||
V_relative_cache = LRU{Tuple{UInt8, UInt8, UInt8}, ComplexF64}(maxsize=(1+l_max)*(1+n_max)^2)
|
||||
|
||||
V_relative = get_sp_V_matrix(V_relative_elem, basis.n1s, basis.l1s, [basis.n2s, basis.l2s]; dtype=ComplexF64, cache=V_relative_cache) + get_sp_V_matrix(V_relative_elem, basis.n2s, basis.l2s, [basis.n1s, basis.l1s]; dtype=ComplexF64, cache=V_relative_cache)
|
||||
U = Moshinsky_transform(basis)
|
||||
V2 = transpose(U) * V_relative * U
|
||||
V_relative = get_sp_V_matrix(V_relative_elem, n1s, l1s; mask=mask1, dtype=ComplexF64, cache=V_relative_cache) + get_sp_V_matrix(V_relative_elem, n2s, l2s; mask=mask2, dtype=ComplexF64, cache=V_relative_cache)
|
||||
U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
|
||||
V2 = U' * V_relative * U
|
||||
|
||||
return V2
|
||||
end
|
||||
|
||||
function get_2p_p1p2_matrix(basis::ho_basis_2B, μ1ω1, μ2ω2; dtype=Float64)
|
||||
# TODO: Cache for integrals
|
||||
integral1(np, lp, n, l) = integral_HO(np, lp, n, l, μ1ω1)
|
||||
integral2(np, lp, n, l) = integral_HO(np, lp, n, l, μ2ω2)
|
||||
|
||||
mat = zeros(dtype, basis.dim, basis.dim)
|
||||
function get_2p_p1p2_matrix(n1s, l1s, n2s, l2s, Λ, μ1ω1, μ2ω2; dtype=Float64)
|
||||
mat = zeros(dtype, length(n1s), length(n1s))
|
||||
Threads.@threads for idx in CartesianIndices(mat)
|
||||
(i, j) = Tuple(idx)
|
||||
val = racahs_reduction_formula(basis.n1s[i], basis.l1s[i], basis.n2s[i], basis.l2s[i], basis.n1s[j], basis.l1s[j], basis.n2s[j], basis.l2s[j], basis.Λ, integral1, integral2)
|
||||
val = racahs_reduction_formula(n1s[i], l1s[i], n2s[i], l2s[i], n1s[j], l1s[j], n2s[j], l2s[j], Λ, μ1ω1, μ2ω2)
|
||||
if !(val ≈ 0); mat[idx] = val; end
|
||||
end
|
||||
return sparse(mat)
|
||||
end
|
||||
|
||||
function get_src_V_matrix(V_of_r, basis::ho_basis_2B, μω, μω_global; atol=10^-6, maxevals=10^5)
|
||||
function get_src_V_matrix(V_of_r, E_max, Λ, μω, μω_global; atol=10^-6, maxevals=10^5)
|
||||
_, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
|
||||
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_elem(l, n1, n2) = V_numerical(V_of_r, l, n1, n2; μω_gen=μω, atol=atol, maxevals=maxevals)
|
||||
V_cache = prealloc_V_cache(basis.E_max, ComplexF64)
|
||||
V_cache = LRU{Tuple{UInt8, UInt8, UInt8}, ComplexF64}(maxsize=(1+l_max)*(1+n_max)^2)
|
||||
|
||||
V1 = get_sp_V_matrix(V_elem, basis.n1s, basis.l1s, [basis.n2s, basis.l2s]; dtype=ComplexF64, cache=V_cache)
|
||||
V2 = get_sp_V_matrix(V_elem, basis.n2s, basis.l2s, [basis.n1s, basis.l1s]; dtype=ComplexF64, cache=V_cache)
|
||||
V1 = get_sp_V_matrix(V_elem, n1s, l1s; mask=mask1, dtype=ComplexF64, cache=V_cache)
|
||||
V2 = get_sp_V_matrix(V_elem, n2s, l2s; mask=mask2, dtype=ComplexF64, cache=V_cache)
|
||||
|
||||
V12 = get_src_V12_matrix(V_of_r, basis, μω_global; atol=atol, maxevals=maxevals)
|
||||
V12 = get_src_V12_matrix(V_of_r, E_max, Λ, μω_global; atol=atol, maxevals=maxevals)
|
||||
|
||||
return V1 + V2 + V12
|
||||
end
|
||||
|
||||
function get_src_V12_matrix(V_of_r, basis::ho_basis_2B, μω_global; atol=10^-6, maxevals=10^5)
|
||||
function get_src_V12_matrix(V_of_r, E_max, Λ, μω_global; atol=10^-6, maxevals=10^5)
|
||||
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
|
||||
l_max = max(maximum(l1s), maximum(l2s))
|
||||
n_max = max(maximum(n1s), maximum(n2s))
|
||||
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
|
||||
|
||||
V_relative_elem(l, n1, n2) = V_numerical(V_of_r, l, n1, n2; μω_gen=μω_global, atol=atol, maxevals=maxevals)
|
||||
V_relative = get_sp_V_matrix(V_relative_elem, basis.n1s, basis.l1s, [basis.n2s, basis.l2s]; dtype=ComplexF64, E_max=basis.E_max)
|
||||
V_relative_cache = LRU{Tuple{UInt8, UInt8, UInt8}, ComplexF64}(maxsize=(1+l_max)*(1+n_max)^2)
|
||||
|
||||
U = Moshinsky_transform(basis)
|
||||
V12 = transpose(U) * V_relative * U
|
||||
V_relative = get_sp_V_matrix(V_relative_elem, n1s, l1s; mask=mask1, dtype=ComplexF64, cache=V_relative_cache)
|
||||
U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
|
||||
V12 = U' * V_relative * U
|
||||
|
||||
return V12
|
||||
end
|
||||
|
||||
"Basis transformation from HO to momentum space"
|
||||
function get_W_matrix(basis_p, basis::ho_basis_2B, μ1ω1, μ2ω2=μ1ω1; weights=true)
|
||||
W = zeros(ComplexF64, length(basis_p), basis.dim)
|
||||
Threads.@threads for idx in CartesianIndices(W)
|
||||
(i1, i2) = Tuple(idx)
|
||||
((j1, j2), (k1, w1), (k2, w2)) = basis_p[i1]
|
||||
if j1 == basis.l1s[i2] && j2 == basis.l2s[i2]
|
||||
elem1 = 1/sqrt(sqrt(μ1ω1)) * (-1)^basis.n1s[i2] * ho_basis(j1, basis.n1s[i2], 1/sqrt(μ1ω1) * k1)
|
||||
elem2 = 1/sqrt(sqrt(μ2ω2)) * (-1)^basis.n2s[i2] * ho_basis(j2, basis.n2s[i2], 1/sqrt(μ2ω2) * k2)
|
||||
W[idx] = elem1 * elem2 * (weights ? w1 * w2 : 1)
|
||||
end
|
||||
end
|
||||
return sparse(W)
|
||||
end
|
||||
|
||||
function get_3b_H_matrix(coord_system::coordinate_system, V_of_r, μω_global, E_max, Λ, m=1.0, kinetic_part=true, potential_part=true; atol=10^-5, maxevals=10^5, verbose=true)
|
||||
if coord_system == jacobi
|
||||
μ1ω1 = μω_global * 1/2
|
||||
μ2ω2 = μω_global * 2
|
||||
μ1 = m * 1/2
|
||||
μ2 = m * 2/3
|
||||
elseif coord_system == src
|
||||
μ1ω1 = μ2ω2 = μω = μω_global * 2
|
||||
μ1 = μ2 = μ = m/2
|
||||
end
|
||||
|
||||
verbose && println("No of threads = ", Threads.nthreads())
|
||||
|
||||
@time "Basis" basis = ho_basis_2B(E_max, Λ)
|
||||
verbose && println("Basis size = ", basis.dim)
|
||||
|
||||
out = spzeros(basis.dim, basis.dim)
|
||||
|
||||
if kinetic_part
|
||||
verbose && println("Constructing KE matrices")
|
||||
@time "T1" out += get_sp_T_matrix(basis.n1s, basis.l1s, [basis.n2s, basis.l2s]; μω_gen=μ1ω1, μ=μ1)
|
||||
@time "T2" out += get_sp_T_matrix(basis.n2s, basis.l2s, [basis.n1s, basis.l1s]; μω_gen=μ2ω2, μ=μ2)
|
||||
if coord_system == src
|
||||
@time "T_cross" out += get_2p_p1p2_matrix(basis, μ1ω1, μ2ω2; dtype=ComplexF64) ./ (2*μ)
|
||||
end
|
||||
end
|
||||
|
||||
if potential_part
|
||||
verbose && println("Constructing PE matrices")
|
||||
if coord_system == jacobi
|
||||
@time "V" out += get_jacobi_V_matrix(V_of_r, basis, μ1ω1, μω_global; atol=atol, maxevals=maxevals)
|
||||
elseif coord_system == src
|
||||
@time "V" out += get_src_V_matrix(V_of_r, basis, μω, μω_global; atol=atol, maxevals=maxevals)
|
||||
end
|
||||
end
|
||||
|
||||
return out
|
||||
end
|
||||
|
|
@ -0,0 +1,30 @@
|
|||
using Arpack, SparseArrays
|
||||
include("ho_basis.jl")
|
||||
|
||||
E_max = 30
|
||||
μω_gen = 0.2
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
|
||||
Va = -2
|
||||
Ra = 2
|
||||
|
||||
μ1 = m * 1/2
|
||||
|
||||
println("No of threads = ", Threads.nthreads())
|
||||
|
||||
Es, n1s, l1s = get_sp_basis(E_max)
|
||||
println("Basis size = ", length(Es))
|
||||
|
||||
println("Constructing KE matrices")
|
||||
@time "T1" T1 = get_sp_T_matrix(n1s, l1s; μω_gen=μω_gen, μ=μ1)
|
||||
|
||||
println("Constructing PE matrices")
|
||||
V1_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; μω_gen=μω_gen)
|
||||
@time "V1" V1 = get_sp_V_matrix(V1_elem, n1s, l1s)
|
||||
|
||||
println("Calculating spectrum")
|
||||
@time "H" H = T1 + V1
|
||||
@time "Eigenvalues" evals, _ = eigs(H, nev=3, ncv=30, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
|
||||
|
||||
display(evals)
|
||||
|
|
@ -1,14 +1,44 @@
|
|||
using Arpack
|
||||
using Arpack, SparseArrays
|
||||
include("ho_basis.jl")
|
||||
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
V_of_r(r) = -2 * exp(-r^2 / 4)
|
||||
Va = -2
|
||||
Ra = 2
|
||||
|
||||
E_max = 40
|
||||
μω_global = 0.3
|
||||
|
||||
H = get_3b_H_matrix(src, V_of_r, μω_global, E_max, Λ, m)
|
||||
# due to simple relative coordinates
|
||||
μω = μω_global * 2
|
||||
μ = m/2
|
||||
|
||||
println("No of threads = ", Threads.nthreads())
|
||||
|
||||
@time "Basis" begin
|
||||
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
|
||||
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
|
||||
mask2 = (n1s .== n1s') .&& (l1s .== l1s')
|
||||
end
|
||||
|
||||
println("Basis size = ", length(Es))
|
||||
|
||||
println("Constructing KE matrices")
|
||||
@time "T1" T1 = get_sp_T_matrix(n1s, l1s; mask=mask1, μω_gen=μω, μ=μ)
|
||||
@time "T2" T2 = get_sp_T_matrix(n2s, l2s; mask=mask2, μω_gen=μω, μ=μ)
|
||||
@time "T_cross" T_cross = get_2p_p1p2_matrix(n1s, l1s, n2s, l2s, Λ, μω, μω) ./ (2*μ)
|
||||
|
||||
println("Constructing PE matrices")
|
||||
V_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; μω_gen=μω)
|
||||
V_relative_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; μω_gen=μω_global)
|
||||
@time "V1" V1 = get_sp_V_matrix(V_elem, n1s, l1s; mask=mask1)
|
||||
@time "V2" V2 = get_sp_V_matrix(V_elem, n2s, l2s; mask=mask2)
|
||||
@time "V relative" V_relative = get_sp_V_matrix(V_relative_elem, n1s, l1s; mask=mask1)
|
||||
@time "Moshinsky brackets" U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
|
||||
@time "V12" V12 = U' * V_relative * U
|
||||
|
||||
println("Calculating spectrum")
|
||||
@time "H" H = T1 + T2 + T_cross + V1 + V2 + V12
|
||||
@time "Eigenvalues" evals, _ = eigs(H, nev=3, ncv=30, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
|
||||
|
||||
display(evals)
|
||||
|
|
@ -1,15 +1,41 @@
|
|||
using Arpack
|
||||
using Arpack, SparseArrays, LRUCache
|
||||
include("ho_basis.jl")
|
||||
|
||||
target_E = 4.07656088827514 - 0.012743522750966718im
|
||||
V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
|
||||
|
||||
μω_global = 0.5 * exp(-2im * pi / 9)
|
||||
E_max = 30
|
||||
μω_global = 0.5 * exp(-2im * pi / 9)
|
||||
|
||||
H = get_3b_H_matrix(src, V_of_r, μω_global, E_max, Λ, m)
|
||||
# due to simple relative coordinates
|
||||
μω = μω_global * 2
|
||||
μ = m/2
|
||||
|
||||
println("No of threads = ", Threads.nthreads())
|
||||
|
||||
@time "Basis" begin
|
||||
Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
|
||||
l_max = max(maximum(l1s), maximum(l2s))
|
||||
n_max = max(maximum(n1s), maximum(n2s))
|
||||
mask1 = (n2s .== n2s') .&& (l2s .== l2s')
|
||||
mask2 = (n1s .== n1s') .&& (l1s .== l1s')
|
||||
end
|
||||
|
||||
println("Basis size = ", length(Es))
|
||||
|
||||
println("Constructing KE matrices")
|
||||
|
||||
@time "T1" T1 = get_sp_T_matrix(n1s, l1s; mask=mask1, μω_gen=μω, μ=μ)
|
||||
@time "T2" T2 = get_sp_T_matrix(n2s, l2s; mask=mask2, μω_gen=μω, μ=μ)
|
||||
@time "T_cross" T_cross = get_2p_p1p2_matrix(n1s, l1s, n2s, l2s, Λ, μω, μω; dtype=ComplexF64) ./ (2*μ)
|
||||
|
||||
println("Constructing PE matrices")
|
||||
|
||||
@time "V" V = get_src_V_matrix(V_of_r, E_max, Λ, μω, μω_global)
|
||||
|
||||
println("Calculating spectrum")
|
||||
@time "H" H = T1 + T2 + T_cross + V
|
||||
@time "Eigenvalues" evals, _ = eigs(H, nev=5, ncv=50, which=:LI, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
|
||||
|
||||
@time "Eigenvalues" evals, _ = eigs(H, nev=5, ncv=50, sigma=target_E, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
|
||||
display(evals)
|
||||
|
|
@ -0,0 +1,30 @@
|
|||
using HDF5
|
||||
|
||||
include("ho_basis_3body.jl")
|
||||
|
||||
@time "Eigenvectors" evals, evecs = eigs(H, nev=1, ncv=10, which=:SR, maxiter=5000, tol=1e-8, ritzvec=true, check=1)
|
||||
|
||||
idx = argmin(real.(evals))
|
||||
|
||||
E = evals[idx]
|
||||
@assert imag(E) ≈ 0 "Energy is not real"
|
||||
E = real(E)
|
||||
println("Exporting E = $E")
|
||||
|
||||
wf = real.(evecs[:, idx])
|
||||
@assert all(imag.(wf) .≈ 0) "Wave function is not real"
|
||||
wf = real.(wf)
|
||||
|
||||
h5open("temp/HO_3b.h5", "w") do fid
|
||||
write_attribute(fid, "E", E)
|
||||
write_attribute(fid, "Λ", Λ)
|
||||
write_attribute(fid, "μ1ω1", μ1ω1)
|
||||
write_attribute(fid, "μ2ω2", μ2ω2)
|
||||
write_attribute(fid, "μ1", μ1)
|
||||
write_attribute(fid, "μ2", μ2)
|
||||
write_dataset(fid, "n1", n1s)
|
||||
write_dataset(fid, "l1", l1s)
|
||||
write_dataset(fid, "n2", n2s)
|
||||
write_dataset(fid, "l2", l2s)
|
||||
write_dataset(fid, "wf", wf)
|
||||
end
|
||||
|
|
@ -1,5 +1,5 @@
|
|||
using LinearAlgebra
|
||||
include("../ho_basis.jl")
|
||||
include("ho_basis.jl")
|
||||
|
||||
l = 0
|
||||
V0 = -10
|
||||
25
math.jl
25
math.jl
|
|
@ -1,5 +1,5 @@
|
|||
using SpecialFunctions, WignerSymbols
|
||||
include("common.jl")
|
||||
include("helper.jl")
|
||||
|
||||
# Gaussian potentials in HO space
|
||||
inv_factorial(n) = Iterators.prod(inv.(1:n))
|
||||
|
|
@ -22,26 +22,19 @@ double_factorial(n::Int) = Iterators.prod(big, n:-2:1)
|
|||
"Gaussian integral for n ∈ Integers (Ref: Wolfram MathWorld + simplifications)"
|
||||
gauss_int(a, n) = double_factorial(n - 1) / (2.0 * a)^((n + 1)/2) * (iseven(n) ? sqrt(π / 2) : 1)
|
||||
|
||||
"Gives ∫dp p u' u where u' and u are HO functions with different l (Ref: worked out in Mathematica)"
|
||||
function integral_HO(np, lp, n, l, μω)
|
||||
"Gives ∫dp p u' u where u' and u' may have different l (Ref: worked out in Mathematica)"
|
||||
function integral(np, lp, n, l, μω)
|
||||
s = [(-1)^(m + mp) * gauss_int(1, 2 * m + 2 * mp + l + lp + 3) * N_lnk(l, n, m) * N_lnk(lp, np, mp) for (m, mp) in Iterators.product(0:n, 0:np)]
|
||||
return 2 * sqrt(μω) * better_sum(s)
|
||||
end
|
||||
|
||||
"Gives <n' l'|| p ||n l> for the HO basis, where integral(np, lp, n, l) is a function that returns ∫dp p u' u"
|
||||
function reduced_matrix_element(np, lp, n, l, integral)::ComplexF64
|
||||
wig::Float64 = wigner3j(Float64, lp, 1, l, 0, 0, 0)
|
||||
if wig == 0
|
||||
return 0
|
||||
else
|
||||
return (-1)^lp * sqrt(2*lp + 1) * sqrt(2*l + 1) * wig * integral(np, lp, n, l)
|
||||
end
|
||||
end
|
||||
"Gives <n' l'|| p ||n l> for the HO basis"
|
||||
reduced_matrix_element(np, lp, n, l, μω) = (-1)^lp * sqrt(2*lp + 1) * sqrt(2*l + 1) * wigner3j(Float64, lp, 1, l, 0, 0, 0) * integral(np, lp, n, l, μω)
|
||||
|
||||
"Matrix element <n1p l1p n2p l2p| p1⋅p2 |n1 l1 n2 l2> (Ref: de-Shalit & Talmi, Eq 15.5), where integral(np, lp, n, l) is a function that returns ∫dp p u' u"
|
||||
function racahs_reduction_formula(n1p, l1p, n2p, l2p, n1, l1, n2, l2, Λ, integral1, integral2)
|
||||
"Matrix element <n1p l1p n2p l2p| p1⋅p2 |n1 l1 n2 l2> (Ref: de-Shalit & Talmi, Eq 15.5)"
|
||||
function racahs_reduction_formula(n1p, l1p, n2p, l2p, n1, l1, n2, l2, Λ, μ1ω1, μ2ω2)
|
||||
val = wigner6j(Float64, l1p, l2p, Λ, l2, l1, 1)
|
||||
if val != 0; val *= reduced_matrix_element(n1p, l1p, n1, l1, integral1); end
|
||||
if val != 0; val *= reduced_matrix_element(n2p, l2p, n2, l2, integral2); end
|
||||
if val != 0; val *= reduced_matrix_element(n1p, l1p, n1, l1, μ1ω1); end
|
||||
if val != 0; val *= reduced_matrix_element(n2p, l2p, n2, l2, μ2ω2); end
|
||||
return (-1)^(l1 + l2p + Λ) * val
|
||||
end
|
||||
54
p_space.jl
54
p_space.jl
|
|
@ -1,6 +1,5 @@
|
|||
using LinearAlgebra
|
||||
using SpecialFunctions, FastGaussQuadrature, QuadGK
|
||||
include("ho_basis.jl")
|
||||
|
||||
function gausslegendre_shifted(a, b, n)
|
||||
scale = (b - a) / 2
|
||||
|
|
@ -60,56 +59,3 @@ function Vl_mat_elem(V_of_r, l, p, q; atol=0, maxevals=10^7, R_cutoff=Inf)
|
|||
(integral, _) = quadgk(integrand, 0, R_cutoff; atol=atol, maxevals=maxevals)
|
||||
return (2 / pi) * integral
|
||||
end
|
||||
|
||||
"Return the Hamiltonian matrix (and the array of weights) for the given system."
|
||||
function get_3b_H_matrix(coord_system::coordinate_system, V_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m=1.0, kinetic_part=true, potential_part=true; atol=10^-5, maxevals=10^5, R_cutoff=16, verbose=true)
|
||||
if coord_system == jacobi
|
||||
μ1ω1 = μω_global * 1/2
|
||||
μ2ω2 = μω_global * 2
|
||||
μ1 = m * 1/2
|
||||
μ2 = m * 2/3
|
||||
else
|
||||
error("Only Jacobi coordinates are implemented")
|
||||
end
|
||||
|
||||
verbose && println("No of threads = ", Threads.nthreads())
|
||||
|
||||
V_l(j, k, kp) = Vl_mat_elem(V_of_r, j, k, kp; atol=atol, maxevals=maxevals, R_cutoff=R_cutoff)
|
||||
|
||||
ks, ws = get_mesh(vertices, subdivisions)
|
||||
weights = repeat(kron(ws, ws), jmax + 1)
|
||||
block_size = length(ks)
|
||||
tri((j1, j2)) = triangle_ineq(j1, j2, Λ)
|
||||
js = collect(Iterators.filter(tri, iter_prod(0:jmax, 0:jmax)))
|
||||
basis = iter_prod(js, zip(ks, ws), zip(ks, ws)) # basis = ((j1, j2), (k1, w1), (k2, w2))
|
||||
basis_size = length(js) * length(ks)^2
|
||||
@assert length(basis) == basis_size "Something wrong with the basis"
|
||||
verbose && println("Basis size = $basis_size")
|
||||
|
||||
out = spzeros(basis_size, basis_size)
|
||||
|
||||
@time "Block diagonal part" begin
|
||||
if kinetic_part & potential_part
|
||||
Hb_blocks = [kron_sum(get_H_matrix((k, kp) -> V_l(j1, k, kp), ks, ws, μ1), get_T_matrix(ks, μ2)) for (j1, _) in js]
|
||||
elseif kinetic_part
|
||||
Hb_blocks = [kron_sum(get_T_matrix(ks, μ1), get_T_matrix(ks, μ2)) for _ in js]
|
||||
elseif potential_part
|
||||
Hb_blocks = [kron_sum(get_V_matrix((k, kp) -> V_l(j1, k, kp), ks, ws), spzeros(block_size, block_size)) for (j1, _) in js]
|
||||
end
|
||||
out += blockdiag(sparse.(Hb_blocks)...)
|
||||
end
|
||||
|
||||
if potential_part
|
||||
basis_ho = ho_basis_2B(E_max, Λ)
|
||||
verbose && println("HO basis size = ", basis_ho.dim)
|
||||
|
||||
@time "V2_HO" V2_HO = get_jacobi_V2_matrix(V_of_r, basis_ho, μω_global)
|
||||
|
||||
@time "W_right" W_right = get_W_matrix(basis, basis_ho, μ1ω1, μ2ω2; weights=true)
|
||||
@time "W_left" W_left = get_W_matrix(basis, basis_ho, μ1ω1, μ2ω2; weights=false)
|
||||
|
||||
@time "V2" out += W_left * V2_HO * transpose(W_right)
|
||||
end
|
||||
|
||||
return (out, weights)
|
||||
end
|
||||
|
|
@ -1,18 +0,0 @@
|
|||
using Arpack
|
||||
include("p_space.jl")
|
||||
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
V_of_r(r) = -2 * exp(-r^2 / 4)
|
||||
|
||||
vertices = [0, 0.5 - 0.3im, 1, 4]
|
||||
subdivisions = [10, 10, 10]
|
||||
jmax = 4
|
||||
|
||||
E_max = 30
|
||||
μω_global = 0.5
|
||||
|
||||
H, _ = get_3b_H_matrix(jacobi, V_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m)
|
||||
|
||||
@time "Eigenvalues" evals, _ = eigs(H, nev=3, ncv=24, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
|
||||
display(evals)
|
||||
|
|
@ -1,20 +0,0 @@
|
|||
using Arpack
|
||||
include("p_space.jl")
|
||||
|
||||
target = 4.0766890719636875 - 0.012758927741074495im
|
||||
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
|
||||
|
||||
vertices = [0, 2 - 0.2im, 3, 4]
|
||||
subdivisions = [15, 10, 10]
|
||||
jmax = 4
|
||||
|
||||
E_max = 40
|
||||
μω_global = 0.5
|
||||
|
||||
H, _ = get_3b_H_matrix(jacobi, V_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m)
|
||||
|
||||
@time "Eigenvalues" evals, _ = eigs(H, sigma=target, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
|
||||
display(evals)
|
||||
|
|
@ -0,0 +1,62 @@
|
|||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"using Plots\n",
|
||||
"include(\"helper.jl\")\n",
|
||||
"include(\"p_space.jl\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"vertices = [0, 1 - 0.5im, 2, 6]\n",
|
||||
"subdivisions = [64, 64, 64]\n",
|
||||
"p, w = get_mesh(vertices, subdivisions)\n",
|
||||
"\n",
|
||||
"scatter(real.(p), imag.(p))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# ComplexScaling-FV: Eq. (54)\n",
|
||||
"V_pq(p, q) = -10 * g1(1, p, q)\n",
|
||||
"\n",
|
||||
"H = get_H_matrix(V_pq, p, w)\n",
|
||||
"evals = eigen(H).values\n",
|
||||
"\n",
|
||||
"E_target = 0.258 - 0.164im\n",
|
||||
"E = nearest(evals, E_target)\n",
|
||||
"\n",
|
||||
"print(\"E = $E\")\n",
|
||||
"scatter(real.(evals), imag.(evals), xlim = (0,1))"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Julia 1.9.0",
|
||||
"language": "julia",
|
||||
"name": "julia-1.9"
|
||||
},
|
||||
"language_info": {
|
||||
"file_extension": ".jl",
|
||||
"mimetype": "application/julia",
|
||||
"name": "julia",
|
||||
"version": "1.9.0"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 2
|
||||
}
|
||||
|
|
@ -0,0 +1,70 @@
|
|||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"using Plots\n",
|
||||
"include(\"helper.jl\")\n",
|
||||
"include(\"p_space.jl\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"vertices = [0, 0.5 - 0.3im, 1, 6]\n",
|
||||
"subdivisions = [64, 64, 64]\n",
|
||||
"p, w = get_mesh(vertices, subdivisions)\n",
|
||||
"\n",
|
||||
"# resonance example from my thesis\n",
|
||||
"V_basis(p, q) = 2*g0(4, p, q) - 3*g0(2, p, q)\n",
|
||||
"\n",
|
||||
"basis_eig = eigen(get_H_matrix(V_basis, p, w))\n",
|
||||
"basis = basis_eig.vectors .* w\n",
|
||||
"\n",
|
||||
"evals = basis_eig.values\n",
|
||||
"E_target = 0.7\n",
|
||||
"E = include(\"helper.jl\")\n",
|
||||
"\n",
|
||||
"print(\"E = $E\")\n",
|
||||
"scatter(real.(evals), imag.(evals), xlim = (0,2))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# ResonanceEC: Eq. (20)\n",
|
||||
"V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))\n",
|
||||
"\n",
|
||||
"H = get_H_matrix(V_system(0.45), p, w)\n",
|
||||
"H_berggren = transpose(basis) * H * basis\n",
|
||||
"\n",
|
||||
"evals = eigvals(H)\n",
|
||||
"scatter(real.(evals), imag.(evals), xlim = (0, 0.5))"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Julia 1.10.2",
|
||||
"language": "julia",
|
||||
"name": "julia-1.10"
|
||||
},
|
||||
"language_info": {
|
||||
"file_extension": ".jl",
|
||||
"mimetype": "application/julia",
|
||||
"name": "julia",
|
||||
"version": "1.10.2"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 2
|
||||
}
|
||||
|
|
@ -1,4 +1,4 @@
|
|||
include("../p_space_3body_resonance.jl")
|
||||
include("../berggren_3body_resonance.jl")
|
||||
|
||||
@time "Eigenvectors" evals, evecs = eigs(H, sigma=target, maxiter=5000, tol=1e-5, ritzvec=true, check=1)
|
||||
|
||||
65
test/misc.jl
65
test/misc.jl
|
|
@ -1,65 +0,0 @@
|
|||
using SparseArrays, LinearAlgebra
|
||||
include("../common.jl")
|
||||
|
||||
#### gram_schmidt ####
|
||||
|
||||
n = 64
|
||||
N = 8
|
||||
|
||||
vecs = [rand(n) + 1im .* rand(n) for _ in 1:N]
|
||||
|
||||
ws = rand(n)
|
||||
ws_mat = spdiagm(ws)
|
||||
|
||||
basis = hcat(vecs...)
|
||||
|
||||
H = rand(n, n)
|
||||
H += transpose(H) # complex symmetric
|
||||
|
||||
H_EC = transpose(basis) * ws_mat * H * ws_mat * basis
|
||||
N_EC = transpose(basis) * ws_mat * basis
|
||||
|
||||
evals = eigvals(H_EC, N_EC)
|
||||
println("Eigenvalues with GEVP:")
|
||||
display(evals)
|
||||
|
||||
ortho_basis = hcat(gram_schmidt!(vecs, ws)...)
|
||||
N_ortho = transpose(ortho_basis) * ws_mat * ortho_basis
|
||||
println("Norm matrix after Gram-Schmidt:")
|
||||
display(round.(N_ortho; digits=2)) # should be ≈I
|
||||
@assert N_ortho ≈ I(N) "Gram-Schmidt did not yield an orthogonal basis"
|
||||
|
||||
H_EC_ortho = transpose(ortho_basis) * ws_mat * H * ws_mat * ortho_basis
|
||||
|
||||
evals_ortho = eigvals(H_EC_ortho)
|
||||
println("Eigenvalues after Gram-Schmidt:")
|
||||
display(evals_ortho)
|
||||
|
||||
@assert evals ≈ evals_ortho "Gram-Schmidt did not preserve the eigenvalues"
|
||||
|
||||
############
|
||||
println("\nRepeat with a redundant basis\n")
|
||||
|
||||
println("Original dimensionality = $(length(vecs))")
|
||||
|
||||
for pow in 6:9
|
||||
noise = rand(n) ./ 10^pow
|
||||
new_vec = vecs[1] + noise
|
||||
push!(vecs, new_vec)
|
||||
end
|
||||
|
||||
println("Dimensionality before Gram-Schmidt = $(length(vecs))")
|
||||
|
||||
ortho_vecs = gram_schmidt!(vecs, ws; verbose=true)
|
||||
ortho_basis = hcat(ortho_vecs...)
|
||||
println("Dimensionality after Gram-Schmidt = $(length(ortho_vecs))")
|
||||
|
||||
H_EC_ortho = transpose(ortho_basis) * ws_mat * H * ws_mat * ortho_basis
|
||||
|
||||
evals_ortho = eigvals(H_EC_ortho)
|
||||
println("Eigenvalues after Gram-Schmidt:")
|
||||
display(evals_ortho)
|
||||
|
||||
@assert isapprox(evals, evals_ortho; atol=1e-3) "Gram-Schmidt did not approximately preserve the eigenvalues"
|
||||
|
||||
######################
|
||||
|
|
@ -1,4 +1,4 @@
|
|||
println("### Test: transpose(U) * U == identity")
|
||||
println("### Test: U' * U == identity")
|
||||
|
||||
using LinearAlgebra
|
||||
|
||||
|
|
@ -9,13 +9,13 @@ E_max = 30
|
|||
|
||||
println("No of threads = ", Threads.nthreads())
|
||||
|
||||
@time "Basis" basis = ho_basis_2B(E_max, Λ)
|
||||
@time "Basis" Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
|
||||
|
||||
println("Basis size = ", basis.dim)
|
||||
println("Basis size = ", length(Es))
|
||||
|
||||
@time "Moshinsky brackets" U = Moshinsky_transform(basis)
|
||||
@time "Moshinsky brackets" U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
|
||||
|
||||
check = transpose(U) * U - I
|
||||
check = U' * U - I
|
||||
|
||||
println("Maximum = ", maximum(abs.(check)))
|
||||
println("Norm = ", sum(check .* conj(check)))
|
||||
|
|
|
|||
|
|
@ -1,57 +0,0 @@
|
|||
using LinearAlgebra, SparseArrays, Plots
|
||||
include("../p_space.jl")
|
||||
include("../ho_basis.jl")
|
||||
include("../berggren.jl")
|
||||
|
||||
println("No of threads = ", Threads.nthreads())
|
||||
atol = 10^-5
|
||||
maxevals = 10^5
|
||||
R_cutoff = 16
|
||||
|
||||
Λ = 0
|
||||
m = 1.0
|
||||
μ = m/2 # due to simple relative coordinates
|
||||
|
||||
vertices = [0, 2 - 0.2im, 3, 4]
|
||||
subdivisions = [15, 10, 10]
|
||||
ks, ws = get_mesh(vertices, subdivisions)
|
||||
|
||||
jmax = 4
|
||||
tri((j1, j2)) = triangle_ineq(j1, j2, Λ)
|
||||
js = collect(Iterators.filter(tri, iter_prod(0:jmax, 0:jmax)))
|
||||
|
||||
basis = iter_prod(js, zip(ks, ws), zip(ks, ws)) # basis = ((j1, j2), (k1, w1), (k2, w2))
|
||||
basis_size = length(js) * length(ks)^2
|
||||
@assert length(basis) == basis_size "Something wrong with the basis"
|
||||
println("Basis size = $basis_size")
|
||||
|
||||
V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
|
||||
V_l(j, k, kp) = Vl_mat_elem(V_of_r, j, k, kp; atol=atol, maxevals=maxevals, R_cutoff=R_cutoff)
|
||||
|
||||
# generate p-space bases (actually identity matrices)
|
||||
@time "p-space bases" ps_bases = [Matrix(spdiagm(1 ./ sqrt.(ws))) for _ in 0:jmax]
|
||||
|
||||
# generate Berggren bases
|
||||
@time "Berggren bases" begin
|
||||
berg_bases = Vector{Matrix{ComplexF64}}(undef, jmax + 1)
|
||||
for j in 0:jmax
|
||||
_, berg_basis = eigen(get_H_matrix((k, kp) -> V_l(j, k, kp), ks, ws, μ); permute=false, scale=false)
|
||||
N_berg = sum(berg_basis.^2 .* ws, dims=1)
|
||||
berg_bases[1 + j] = berg_basis ./ transpose(sqrt.(N_berg))
|
||||
end
|
||||
end
|
||||
|
||||
@time "BB tranform matrix" begin
|
||||
U_blocks = [kron(berg_bases[1 + j1] .* ws, berg_bases[1 + j2] .* ws) for (j1, j2) in js]
|
||||
U = blockdiag(sparse.(U_blocks)...)
|
||||
end
|
||||
|
||||
@time "In p-space" T_cross_PS = get_2p_p1p2_matrix(length(ks), js, Λ, ps_bases, ps_bases, ws)
|
||||
@time "In BB" T_cross_BB = get_2p_p1p2_matrix(length(ks), js, Λ, berg_bases, berg_bases, ws)
|
||||
|
||||
@time "Basis transform" T_cross_transformed = transpose(U) * T_cross_PS * U
|
||||
|
||||
diff = abs.(T_cross_transformed .- T_cross_BB)
|
||||
println("Max error = $(maximum(diff))")
|
||||
|
||||
#@time "In in HO" T_cross_HO = get_2p_p1p2_matrix(n1s, l1s, n2s, l2s, Λ, μω, μω; dtype=ComplexF64)
|
||||
|
|
@ -1,21 +0,0 @@
|
|||
using Plots
|
||||
include("../common.jl")
|
||||
include("../p_space.jl")
|
||||
|
||||
vertices = [0, 1 - 0.5im, 2, 6]
|
||||
subdivisions = [64, 64, 64]
|
||||
p, w = get_mesh(vertices, subdivisions)
|
||||
|
||||
scatter(real.(p), imag.(p))
|
||||
|
||||
# ComplexScaling-FV: Eq. (54)
|
||||
V_pq(p, q) = -10 * g1(1, p, q)
|
||||
|
||||
H = get_H_matrix(V_pq, p, w)
|
||||
evals = eigen(H).values
|
||||
|
||||
E_target = 0.258 - 0.164im
|
||||
E = nearest(evals, E_target)
|
||||
|
||||
print("E = $E")
|
||||
scatter(real.(evals), imag.(evals), xlim = (0,1))
|
||||
|
|
@ -1,30 +0,0 @@
|
|||
using Plots
|
||||
include("../common.jl")
|
||||
include("../p_space.jl")
|
||||
|
||||
vertices = [0, 0.5 - 0.3im, 1, 6]
|
||||
subdivisions = [64, 64, 64]
|
||||
p, w = get_mesh(vertices, subdivisions)
|
||||
|
||||
# resonance example from my thesis
|
||||
V_basis(p, q) = 2*g0(4, p, q) - 3*g0(2, p, q)
|
||||
|
||||
basis_eig = eigen(get_H_matrix(V_basis, p, w))
|
||||
basis = basis_eig.vectors .* w
|
||||
|
||||
basis_evals = basis_eig.values
|
||||
E_target = 0.7
|
||||
E = nearest(basis_evals, E_target)
|
||||
|
||||
print("pole E = $E")
|
||||
|
||||
# ResonanceEC: Eq. (20)
|
||||
V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))
|
||||
|
||||
H = get_H_matrix(V_system(0.45), p, w)
|
||||
H_berggren = transpose(basis) * H * basis
|
||||
|
||||
evals = eigvals(H)
|
||||
scatter(real.(evals), imag.(evals), label="E")
|
||||
scatter!(real.(basis_evals), imag.(basis_evals), label="basis", markershape=:x)
|
||||
xlims!((0, 1))
|
||||
Loading…
Reference in New Issue