Separate struct for system parameters
This commit is contained in:
ysyapa
parent
247d564071
commit
85b0478095
|
|
@ -5,11 +5,7 @@ using TensorOperations, KrylovKit, LinearAlgebra, CUDA, CUDA.CUTENSOR
|
|||
|
||||
"A Hamiltonian that can be applied to a vector"
|
||||
struct Hamiltonian{T}
|
||||
d::Int
|
||||
n::Int
|
||||
N::Int
|
||||
L::T
|
||||
μ::T
|
||||
s::system{T}
|
||||
K_partial::Matrix{Complex{T}}
|
||||
K_diag::Union{CuTensor{Complex{T}},Nothing}
|
||||
K_mixed::Union{CuTensor{Complex{T}},Nothing}
|
||||
|
|
@ -17,12 +13,12 @@ struct Hamiltonian{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}
|
||||
function Hamiltonian{T}(s::system{T}, V_twobody::Function, ϕ::Real, 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)
|
||||
k = -s.N÷2:s.N÷2-1
|
||||
Vs = calculate_Vs(V_twobody, s, convert(T, ϕ), n_image)
|
||||
hermitian = ϕ == 0.0
|
||||
K_partial = (exp(-im * ϕ) * im / sqrt(2 * μ)) .* ∂_1DOF.(L, N, k, k')
|
||||
K_partial = (exp(-im * convert(T, ϕ)) * im / sqrt(2 * s.μ)) .* ∂_1DOF.(Ref(s), k, k')
|
||||
K_diag = nothing
|
||||
K_mixed = nothing
|
||||
if mode == gpu_cutensor
|
||||
|
|
@ -30,15 +26,15 @@ struct Hamiltonian{T}
|
|||
K_mixed = CuTensor(CuArray(K_partial), ['a', 'A']) * CuTensor(CuArray(K_partial), ['b', 'B'])
|
||||
Vs = CuArray(Vs)
|
||||
end
|
||||
return new{T}(d, n, N, L, μ, K_partial, K_diag, K_mixed, Vs, hermitian, mode)
|
||||
return new{T}(s, K_partial, K_diag, K_mixed, Vs, hermitian, mode)
|
||||
end
|
||||
end
|
||||
|
||||
Base.size(H::Hamiltonian, i::Int)::Int = (i == 1 || i == 2) ? H.N^(H.d * (H.n - 1)) : throw(ArgumentError("Hamiltonian only has 2 dimesions"))
|
||||
Base.size(H::Hamiltonian, i::Int)::Int = (i == 1 || i == 2) ? H.s.N^(H.s.d * (H.s.n - 1)) : throw(ArgumentError("Hamiltonian only has 2 dimesions"))
|
||||
Base.size(H::Hamiltonian)::Dims{2} = (size(H, 1), size(H, 2))
|
||||
|
||||
"Dimensions of a vector to which 'H' can be applied"
|
||||
vectorDims(H::Hamiltonian)::Dims = tuple(fill(H.N, H.d * (H.n - 1))...)
|
||||
vectorDims(H::Hamiltonian)::Dims = tuple(fill(H.s.N, H.s.d * (H.s.n - 1))...)
|
||||
|
||||
"Apply 'H' on 'v' and store the result in 'out' using the 'cpu_tensor' backend"
|
||||
function LinearAlgebra.mul!(out::Array{Complex{T}}, H::Hamiltonian{T}, v::Array{Complex{T}})::Array{Complex{T}} where {T<:Float}
|
||||
|
|
@ -46,13 +42,13 @@ function LinearAlgebra.mul!(out::Array{Complex{T}}, H::Hamiltonian{T}, v::Array{
|
|||
# apply V operator
|
||||
@. out = H.Vs * v
|
||||
# apply K opereator
|
||||
coords = H.n - 1
|
||||
nconList_v_template = -collect(1:H.d*(coords))
|
||||
for dim = 1:H.d
|
||||
coords = H.s.n - 1
|
||||
nconList_v_template = -collect(1:H.s.d*(coords))
|
||||
for dim = 1:H.s.d
|
||||
for coord1 = 1:coords
|
||||
for coord2 = 1:coord1
|
||||
i1 = which_index(H.n, dim, coord1)
|
||||
i2 = which_index(H.n, dim, coord2)
|
||||
i1 = which_index(H.s, dim, coord1)
|
||||
i2 = which_index(H.s, dim, coord2)
|
||||
nconList_1 = [-i1, 1]
|
||||
nconList_2 = [-i2, 2]
|
||||
nconList_v = copy(nconList_v_template)
|
||||
|
|
@ -85,15 +81,15 @@ function LinearAlgebra.mul!(out::CuArray{Complex{T}}, H::Hamiltonian{T}, v::CuAr
|
|||
NVTX.@range "V" @. out = H.Vs * v
|
||||
synchronize(ctx)
|
||||
# apply K opereator
|
||||
coords = H.n - 1
|
||||
inds_template = ('a' - 1) .+ collect(1:H.d*(coords))
|
||||
coords = H.s.n - 1
|
||||
inds_template = ('a' - 1) .+ collect(1:H.s.d*(coords))
|
||||
v_t = CuTensor(v, copy(inds_template))
|
||||
out_t = CuTensor(out, copy(inds_template))
|
||||
for dim = 1:H.d
|
||||
for dim = 1:H.s.d
|
||||
for coord1 = 1:coords
|
||||
for coord2 = 1:coord1
|
||||
i1 = which_index(H.n, dim, coord1)
|
||||
i2 = which_index(H.n, dim, coord2)
|
||||
i1 = which_index(H.s, dim, coord1)
|
||||
i2 = which_index(H.s, dim, coord2)
|
||||
@assert v_t.inds == inds_template "v indices permuted"
|
||||
if i1 == i2
|
||||
@assert H.K_diag.inds[2] == 'A' "K_diag indices permuted"
|
||||
|
|
|
|||
|
|
@ -29,9 +29,10 @@ N=10
|
|||
n_image=1
|
||||
μ=0.5
|
||||
|
||||
for L::T in 5.0:14.0
|
||||
for L in 5.0:14.0
|
||||
println("Constructing H operator...")
|
||||
@time H=Hamiltonian{T}(V_test,3,3,N,L,convert(T,0),convert(T,μ),n_image,mode)
|
||||
s=system{T}(3,3,N,L,μ)
|
||||
@time H=Hamiltonian{T}(s,V_test,0,n_image,mode)
|
||||
println("Applying H 1000 times...")
|
||||
if GPU_mode
|
||||
v=CUDA.rand(Complex{T},vectorDims(H)...)
|
||||
|
|
|
|||
61
common.jl
61
common.jl
|
|
@ -1,53 +1,66 @@
|
|||
Float = Union{Float32,Float64}
|
||||
|
||||
"A few-body system defined by its physical parameters"
|
||||
struct system{T}
|
||||
d::Int
|
||||
n::Int
|
||||
N::Int
|
||||
L::T
|
||||
μ::T
|
||||
|
||||
function system{T}(d::Int, n::Int, N::Int, L::Real, μ::Real) where {T<:Float}
|
||||
return new{T}(d, n, N, convert(T, L), convert(T, μ))
|
||||
end
|
||||
end
|
||||
|
||||
norm_square(x::Array{Int})::Int = sum(x .* x)
|
||||
|
||||
"Eq (46): Partial derivative matrix element for 1 degree of freedom"
|
||||
function ∂_1DOF(L::T, N::Int, k::Int, l::Int)::Complex{T} where {T<:Float}
|
||||
function ∂_1DOF(s::system{T}, k::Int, l::Int)::Complex{T} where {T<:Float}
|
||||
if k == l
|
||||
return -im * (π / L)
|
||||
return -im * (π / s.L)
|
||||
else
|
||||
return (π / L) * (-1)^(k - l) * exp(-im * π * (k - l) / N) / sin(π * (k - l) / N)
|
||||
return (π / s.L) * (-1)^(k - l) * exp(-im * π * (k - l) / s.N) / sin(π * (k - l) / s.N)
|
||||
end
|
||||
end
|
||||
|
||||
"Which index (dimension of the multidimensional array) corresponds to spatial dimension 'dim' and particle 'p'?"
|
||||
which_index(n::Int, dim::Int, p::Int)::Int = (dim - 1) * (n - 1) + p
|
||||
which_index(s::system, dim::Int, p::Int)::Int = (dim - 1) * (s.n - 1) + p
|
||||
|
||||
"Δk (distance in terms of lattice paramter) between two particles along the given dimension"
|
||||
function get_Δk(n::Int, N::Int, i::CartesianIndex, dim::Int, p1::Int, p2::Int)::Int
|
||||
function get_Δk(s::system, i::CartesianIndex, dim::Int, p1::Int, p2::Int)::Int
|
||||
if p1 == p2
|
||||
return 0
|
||||
elseif p1 == n
|
||||
return -(i[which_index(n, dim, p2)] - N ÷ 2 - 1)
|
||||
elseif p2 == n
|
||||
return i[which_index(n, dim, p1)] - N ÷ 2 - 1
|
||||
elseif p1 == s.n
|
||||
return -(i[which_index(s, dim, p2)] - s.N ÷ 2 - 1)
|
||||
elseif p2 == s.n
|
||||
return i[which_index(s, dim, p1)] - s.N ÷ 2 - 1
|
||||
else
|
||||
return i[which_index(n, dim, p1)] - i[which_index(n, dim, p2)]
|
||||
return i[which_index(s, dim, p1)] - i[which_index(s, dim, p2)]
|
||||
end
|
||||
end
|
||||
|
||||
"Calculate diagonal elements of the V matrix"
|
||||
function calculate_Vs(V_twobody::Function, d::Int, n::Int, N::Int, L::T, ϕ::T, n_image::Int)::Array{Complex{T}} where {T<:Float}
|
||||
coeff² = (exp(im * ϕ) * L / N)^2
|
||||
images = collect.(Iterators.product(fill(-n_image:n_image, d)...)) # TODO: Learn how to use tuples instead of vectors
|
||||
Vs = zeros(Complex{T}, fill(N, d * (n - 1))...)
|
||||
function calculate_Vs(V_twobody::Function, s::system{T}, ϕ::T, n_image::Int)::Array{Complex{T}} where {T<:Float}
|
||||
coeff² = (exp(im * ϕ) * s.L / s.N)^2
|
||||
images = collect.(Iterators.product(fill(-n_image:n_image, s.d)...)) # TODO: Learn how to use tuples instead of vectors
|
||||
Vs = zeros(Complex{T}, fill(s.N, s.d * (s.n - 1))...)
|
||||
Threads.@threads for i in CartesianIndices(Vs)
|
||||
for p1 in 1:n
|
||||
for p2 in (p1 + 1):n
|
||||
min_Δk = Array{Int}(undef, d)
|
||||
for dim in 1:d
|
||||
Δk = get_Δk(n, N, i, dim, p1, p2)
|
||||
if Δk > N ÷ 2
|
||||
min_Δk[dim] = Δk - N
|
||||
elseif Δk < -N ÷ 2
|
||||
min_Δk[dim] = Δk + N
|
||||
for p1 in 1:s.n
|
||||
for p2 in (p1 + 1):s.n
|
||||
min_Δk = Array{Int}(undef, s.d)
|
||||
for dim in 1:s.d
|
||||
Δk = get_Δk(s, i, dim, p1, p2)
|
||||
if Δk > s.N ÷ 2
|
||||
min_Δk[dim] = Δk - s.N
|
||||
elseif Δk < -s.N ÷ 2
|
||||
min_Δk[dim] = Δk + s.N
|
||||
else
|
||||
min_Δk[dim] = Δk
|
||||
end
|
||||
end
|
||||
for image in images
|
||||
Δk² = norm_square(min_Δk .- (N .* image))
|
||||
Δk² = norm_square(min_Δk .- (s.N .* image))
|
||||
Vs[i] += V_twobody(Δk² * coeff²)
|
||||
end
|
||||
end
|
||||
|
|
|
|||
|
|
@ -24,12 +24,13 @@
|
|||
"d = 3\n",
|
||||
"n = 3\n",
|
||||
"N = 6\n",
|
||||
"L::T = 12\n",
|
||||
"ϕ::T = 0.0\n",
|
||||
"μ::T = 0.5\n",
|
||||
"L = 12\n",
|
||||
"ϕ = 0.0\n",
|
||||
"μ = 0.5\n",
|
||||
"n_imag = 1\n",
|
||||
"\n",
|
||||
"H = Hamiltonian{T}(V_gauss, d, n, N, L, ϕ, μ, n_imag, mode)\n",
|
||||
"s = system{T}(d, n, N, L, μ)\n",
|
||||
"H = Hamiltonian{T}(s, V_gauss, ϕ, n_imag, mode)\n",
|
||||
"@time evals, evecs, info = eig(H, 5)\n",
|
||||
"print(info.numops, \" operations : \")\n",
|
||||
"println(evals)"
|
||||
|
|
@ -49,12 +50,13 @@
|
|||
"d = 3\n",
|
||||
"n = 2\n",
|
||||
"N = 32\n",
|
||||
"L::T = 16\n",
|
||||
"ϕ::T = 0.5\n",
|
||||
"μ::T = 0.5\n",
|
||||
"L = 16\n",
|
||||
"ϕ = 0.5\n",
|
||||
"μ = 0.5\n",
|
||||
"n_imag = 0\n",
|
||||
"\n",
|
||||
"H = Hamiltonian{T}(V_gauss, d, n, N, L, ϕ, μ, n_imag, mode)\n",
|
||||
"s = system{T}(d, n, N, L, μ)\n",
|
||||
"H = Hamiltonian{T}(s, V_gauss, ϕ, n_imag, mode)\n",
|
||||
"@time evals, evecs, info = eig(H, 20)\n",
|
||||
"print(info.numops, \" operations : \")\n",
|
||||
"print(evals)\n",
|
||||
|
|
|
|||
Loading…
Reference in New Issue