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

Simplify Hamiltonian struct

This commit is contained in:
ysyapa 2023-04-20 21:04:08 -04:00
parent f002b36820
commit 247d564071
1 changed files with 13 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
Hamiltonian.jl
View File

@ -10,26 +10,27 @@ struct Hamiltonian{T}
N::Int
L::T
μ::T
∂1 # Matrix{Complex{T}} or Nothing
K_diag # CuTensor{Complex{T}} or Nothing
K_mixed # CuTensor{Complex{T}} or Nothing
Vs # Array{Complex{T}} or CuArray{Complex{T}}
K_partial::Matrix{Complex{T}}
K_diag::Union{CuTensor{Complex{T}},Nothing}
K_mixed::Union{CuTensor{Complex{T}},Nothing}
Vs::Union{Array{Complex{T}},CuArray{Complex{T}}}
hermitian::Bool
mode::Hamiltonian_backend
function Hamiltonian{T}(V_twobody::Function, d::Int, n::Int, N::Int, L::T, ϕ::T, μ::T, n_image::Int, mode::Hamiltonian_backend) where {T<:Float}
@assert mode != gpu_cutensor || CUDA.functional() && CUDA.has_cuda() && CUDA.has_cuda_gpu() "CUDA not available"
k = -N÷2:N÷2-1
Vs = calculate_Vs(V_twobody, d, n, N, L, ϕ, n_image)
hermitian = ϕ == 0.0
if mode == cpu_tensor
∂1 = exp(-im * ϕ) .* ∂_1DOF.(L, N, k, k')
return new{T}(d, n, N, L, μ, ∂1, nothing, nothing, Vs, hermitian, mode)
elseif mode == gpu_cutensor
K_partial = (exp(-im * ϕ) * im / sqrt(2 * μ)) .* ∂_1DOF.(L, N, k, k')
K_partial = (exp(-im * ϕ) * im / sqrt(2 * μ)) .* ∂_1DOF.(L, N, k, k')
K_diag = nothing
K_mixed = nothing
if mode == gpu_cutensor
K_diag = CuTensor(CuArray(K_partial * K_partial), ['a', 'A'])
K_mixed = CuTensor(CuArray(K_partial), ['a', 'A']) * CuTensor(CuArray(K_partial), ['b', 'B'])
return new{T}(d, n, N, L, μ, nothing, K_diag, K_mixed, CuArray(Vs), hermitian, mode)
Vs = CuArray(Vs)
end
return new{T}(d, n, N, L, μ, K_partial, K_diag, K_mixed, Vs, hermitian, mode)
end
end
@ -45,7 +46,6 @@ function LinearAlgebra.mul!(out::Array{Complex{T}}, H::Hamiltonian{T}, v::Array{
# 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
@ -62,8 +62,8 @@ function LinearAlgebra.mul!(out::Array{Complex{T}}, H::Hamiltonian{T}, v::Array{
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)
v_new = @ncon((H.K_partial, H.K_partial, v), (nconList_1, nconList_2, nconList_v))
out = axpy!(1, v_new, out)
end
end
end