74 lines
2.6 KiB
Julia
74 lines
2.6 KiB
Julia
include("common.jl")
|
|
using TensorOperations, KrylovKit, LinearAlgebra
|
|
|
|
"A Hamiltonian that can be applied to a vector"
|
|
struct HOperator{T}
|
|
d::Int
|
|
n::Int
|
|
N::Int
|
|
L::T
|
|
μ::T
|
|
∂1::Matrix{Complex{T}}
|
|
Vs::Array{Complex{T}}
|
|
hermitian::Bool
|
|
function HOperator{T}(V_twobody::Function, d::Int, n::Int, N::Int, L::T, ϕ::T, μ::T, n_image::Int) where {T<:Float}
|
|
k = -N÷2:N÷2-1
|
|
∂1 = exp(-im * ϕ) .* ∂_1DOF.(L, N, k, k')
|
|
Vs = calculate_Vs(V_twobody, d, n, N, L, ϕ, n_image)
|
|
return new{T}(d, n, N, L, μ, ∂1, Vs, ϕ == 0.0)
|
|
end
|
|
end
|
|
|
|
Base.size(H::HOperator, i::Int)::Int = (i == 1 || i == 2) ? H.N^(H.d * (H.n - 1)) : throw(ArgumentError("HOperator only has 2 dimesions"))
|
|
Base.size(H::HOperator)::Dims{2} = (size(H, 1), size(H, 2))
|
|
|
|
"Dimensions of a vector to which H can be applied"
|
|
vectorDims(H::HOperator)::Dims = tuple(fill(H.N, H.d * (H.n - 1))...)
|
|
|
|
"Apply H on v and store the result in out"
|
|
function LinearAlgebra.mul!(out::Array{Complex{T}}, H::HOperator{T}, v::Array{Complex{T}})::Array{Complex{T}} where {T<:Float}
|
|
#LinearMaps.check_dim_mul(out,H,v) --- dimensions don't match
|
|
# apply V operator
|
|
@. out = H.Vs * v
|
|
# apply K opereator
|
|
coeff = -1 / (2 * H.μ)
|
|
coords = H.n - 1
|
|
nconList_v_template = -collect(1:H.d*(coords))
|
|
for dim = 1:H.d
|
|
for coord1 = 1:coords
|
|
for coord2 = 1:coord1
|
|
i1 = which_index(H.n, dim, coord1)
|
|
i2 = which_index(H.n, dim, coord2)
|
|
nconList_1 = [-i1, 1]
|
|
nconList_2 = [-i2, 2]
|
|
nconList_v = copy(nconList_v_template)
|
|
if i1 == i2
|
|
nconList_2[1] = 1
|
|
else
|
|
nconList_v[i1] = 1
|
|
end
|
|
nconList_v[i2] = 2
|
|
v_new = @ncon((H.∂1, H.∂1, v), (nconList_1, nconList_2, nconList_v))
|
|
out = axpy!(coeff, v_new, out)
|
|
end
|
|
end
|
|
end
|
|
return out
|
|
end
|
|
|
|
"Apply H on v and return the result"
|
|
function (H::HOperator{T})(v::Array{Complex{T}})::Array{Complex{T}} where {T<:Float}
|
|
out = similar(v)
|
|
return mul!(out, H, v)
|
|
end
|
|
|
|
tolerance = 1e-6
|
|
|
|
"Wrapper for KrylovKit.eigsolve"
|
|
function eig(H::HOperator{T}, levels::Int; resonances = !H.hermitian)::Tuple{Vector{Complex{T}},Any,Any} where {T<:Float}
|
|
x₀ = rand(Complex{T}, vectorDims(H))
|
|
evals, evecs, info = eigsolve(H, x₀, levels, resonances ? :LI : :SR; ishermitian = H.hermitian, tol = tolerance)
|
|
resonances || info.converged < levels && throw(error("Not enough convergence")) # don't check convergence for resonances
|
|
return evals, evecs, info
|
|
end
|