Compare commits
No commits in common. "main" and "complex-scaling" have entirely different histories.
main
...
complex-sc
|
|
@ -1,15 +1,3 @@
|
|||
# VSCode
|
||||
.vscode/
|
||||
|
||||
# HPC scripts and logs
|
||||
hpc/
|
||||
|
||||
# Calculation outputs
|
||||
*.dat
|
||||
*.csv
|
||||
*.hdf5
|
||||
*.out
|
||||
|
||||
# Temporary and scratch files
|
||||
temp/
|
||||
scratch/
|
||||
|
|
|
|||
|
|
@ -0,0 +1,73 @@
|
|||
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
|
||||
|
|
@ -0,0 +1,98 @@
|
|||
include("common.jl")
|
||||
using KrylovKit, LinearAlgebra, CUDA, CUDA.CUTENSOR
|
||||
|
||||
@assert CUDA.functional() && CUDA.has_cuda() && CUDA.has_cuda_gpu() "CUDA not available"
|
||||
|
||||
"A Hamiltonian that can be applied to a vector"
|
||||
struct HOperator{T}
|
||||
d::Int
|
||||
n::Int
|
||||
N::Int
|
||||
K_diag::CuTensor{Complex{T}}
|
||||
K_mixed::CuTensor{Complex{T}}
|
||||
Vs::CuArray{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
|
||||
K_partial = (exp(-im * ϕ) * im / sqrt(2 * μ)) .* ∂_1DOF.(L, N, k, k')
|
||||
K_diag = CuTensor(CuArray(K_partial * K_partial), ['a', 'A'])
|
||||
K_mixed = CuTensor(CuArray(K_partial), ['a', 'A']) * CuTensor(CuArray(K_partial), ['b', 'B'])
|
||||
Vs = calculate_Vs(V_twobody, d, n, N, L, ϕ, n_image)
|
||||
return new{T}(d, n, N, K_diag, K_mixed, 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))...)
|
||||
|
||||
"cuTENSOR contraction and accumulation (C = A * B + C)"
|
||||
function contract_accumulate!(C::CuTensor, A::CuTensor, B::CuTensor)::CuTensor
|
||||
CUTENSOR.contraction!(one(eltype(C)), A.data, A.inds, CUTENSOR.CUTENSOR_OP_IDENTITY, B.data, B.inds, CUTENSOR.CUTENSOR_OP_IDENTITY,
|
||||
one(eltype(C)), C.data, C.inds, CUTENSOR.CUTENSOR_OP_IDENTITY, CUTENSOR.CUTENSOR_OP_IDENTITY)
|
||||
return C
|
||||
end
|
||||
|
||||
"Apply H on v and store the result in out"
|
||||
function LinearAlgebra.mul!(out::CuArray{Complex{T}}, H::HOperator{T}, v::CuArray{Complex{T}})::CuArray{Complex{T}} where {T<:Float}
|
||||
#LinearMaps.check_dim_mul(out,H,v) --- dimensions don't match
|
||||
ctx = context()
|
||||
# apply V operator
|
||||
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))
|
||||
v_t = CuTensor(v, copy(inds_template))
|
||||
out_t = CuTensor(out, copy(inds_template))
|
||||
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)
|
||||
@assert v_t.inds == inds_template "v indices permuted"
|
||||
if i1 == i2
|
||||
@assert H.K_diag.inds[2] == 'A' "K_diag indices permuted"
|
||||
H.K_diag.inds[1] = 'a' - 1 + i1
|
||||
v_t.inds[i1] = 'A'
|
||||
#synchronize(ctx)
|
||||
NVTX.@range "K-diag" out_t = contract_accumulate!(out_t, H.K_diag, v_t)
|
||||
v_t.inds[i1] = 'a' - 1 + i1
|
||||
else
|
||||
@assert H.K_mixed.inds[2] == 'A' && H.K_mixed.inds[4] == 'B' "K_mixed indices permuted"
|
||||
H.K_mixed.inds[1] = 'a' - 1 + i1
|
||||
H.K_mixed.inds[3] = 'a' - 1 + i2
|
||||
# OPTIMIZE: A and B can be swapped
|
||||
v_t.inds[i1] = 'A'
|
||||
v_t.inds[i2] = 'B'
|
||||
#synchronize(ctx)
|
||||
NVTX.@range "K-mixed" out_t = contract_accumulate!(out_t, H.K_mixed, v_t)
|
||||
v_t.inds[i1] = 'a' - 1 + i1
|
||||
v_t.inds[i2] = 'a' - 1 + i2
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
@assert out_t.inds == inds_template "out indices permuted"
|
||||
synchronize(ctx)
|
||||
return out_t.data
|
||||
end
|
||||
|
||||
"Apply H on v and return the result"
|
||||
function (H::HOperator{T})(v::CuArray{Complex{T}})::CuArray{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₀ = CUDA.rand(Complex{T}, vectorDims(H)...) # ... added
|
||||
synchronize()
|
||||
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
|
||||
145
Hamiltonian.jl
145
Hamiltonian.jl
|
|
@ -1,145 +0,0 @@
|
|||
include("common.jl")
|
||||
using TensorOperations, KrylovKit, LinearAlgebra, CUDA, cuTENSOR, NVTX
|
||||
|
||||
@enum Hamiltonian_backend cpu_tensor gpu_cutensor
|
||||
|
||||
"A Hamiltonian that can be applied to a vector"
|
||||
struct Hamiltonian{T}
|
||||
s::system{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}(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 = -s.N÷2:s.N÷2-1
|
||||
Vs = calculate_Vs(s, V_twobody, convert(T, ϕ), n_image)
|
||||
hermitian = ϕ == 0.0
|
||||
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
|
||||
K_diag = CuTensor(CuArray(K_partial * K_partial), ['a', 'A'])
|
||||
K_mixed = CuTensor(CuArray(K_partial), ['a', 'A']) * CuTensor(CuArray(K_partial), ['b', 'B'])
|
||||
Vs = CuArray(Vs)
|
||||
end
|
||||
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.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.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}
|
||||
#LinearMaps.check_dim_mul(out,H,v) --- dimensions don't match
|
||||
# apply V operator
|
||||
@. out = H.Vs * v
|
||||
# apply K opereator
|
||||
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.s, dim, coord1)
|
||||
i2 = which_index(H.s, 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.K_partial, H.K_partial, v), (nconList_1, nconList_2, nconList_v))
|
||||
out = axpy!(1, v_new, out)
|
||||
end
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
"cuTENSOR contraction and accumulation (C = A * B + C)"
|
||||
function contract_accumulate!(C::CuTensor, A::CuTensor, B::CuTensor)::CuTensor
|
||||
cuTENSOR.contraction!(one(eltype(C)), A.data, A.inds, cuTENSOR.CUTENSOR_OP_IDENTITY, B.data, B.inds, cuTENSOR.CUTENSOR_OP_IDENTITY,
|
||||
one(eltype(C)), C.data, C.inds, cuTENSOR.CUTENSOR_OP_IDENTITY, cuTENSOR.CUTENSOR_OP_IDENTITY)
|
||||
return C
|
||||
end
|
||||
|
||||
"Apply 'H' on 'v' and store the result in 'out' using the 'gpu_cutensor' backend"
|
||||
function LinearAlgebra.mul!(out::CuArray{Complex{T}}, H::Hamiltonian{T}, v::CuArray{Complex{T}})::CuArray{Complex{T}} where {T<:Float}
|
||||
#LinearMaps.check_dim_mul(out,H,v) --- dimensions don't match
|
||||
ctx = context()
|
||||
# apply V operator
|
||||
NVTX.@range "V" @. out = H.Vs * v
|
||||
synchronize(ctx)
|
||||
# apply K opereator
|
||||
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.s.d
|
||||
for coord1 = 1:coords
|
||||
for coord2 = 1:coord1
|
||||
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"
|
||||
H.K_diag.inds[1] = 'a' - 1 + i1
|
||||
v_t.inds[i1] = 'A'
|
||||
#synchronize(ctx)
|
||||
NVTX.@range "K-diag" out_t = contract_accumulate!(out_t, H.K_diag, v_t)
|
||||
v_t.inds[i1] = 'a' - 1 + i1
|
||||
else
|
||||
@assert H.K_mixed.inds[2] == 'A' && H.K_mixed.inds[4] == 'B' "K_mixed indices permuted"
|
||||
H.K_mixed.inds[1] = 'a' - 1 + i1
|
||||
H.K_mixed.inds[3] = 'a' - 1 + i2
|
||||
# OPTIMIZE: A and B can be swapped
|
||||
v_t.inds[i1] = 'A'
|
||||
v_t.inds[i2] = 'B'
|
||||
#synchronize(ctx)
|
||||
NVTX.@range "K-mixed" out_t = contract_accumulate!(out_t, H.K_mixed, v_t)
|
||||
v_t.inds[i1] = 'a' - 1 + i1
|
||||
v_t.inds[i2] = 'a' - 1 + i2
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
@assert out_t.inds == inds_template "out indices permuted"
|
||||
synchronize(ctx)
|
||||
return out_t.data
|
||||
end
|
||||
|
||||
"Apply 'H' on 'v' and return the result"
|
||||
function (H::Hamiltonian)(v)
|
||||
out = similar(v)
|
||||
return mul!(out, H, v)
|
||||
end
|
||||
|
||||
tolerance = 1e-6
|
||||
|
||||
"Wrapper for KrylovKit.eigsolve"
|
||||
function eig(H::Hamiltonian{T}, levels::Int; resonances = !H.hermitian)::Tuple{Vector,Vector,KrylovKit.ConvergenceInfo} where {T<:Float}
|
||||
if H.mode == cpu_tensor
|
||||
x₀ = rand(Complex{T}, vectorDims(H)...)
|
||||
elseif H.mode == gpu_cutensor
|
||||
x₀ = CUDA.rand(Complex{T}, vectorDims(H)...)
|
||||
synchronize()
|
||||
end
|
||||
evals, evecs, info = eigsolve(H, x₀, levels, resonances ? :LI : :SR; ishermitian = H.hermitian, tol = tolerance, krylovdim = levels * 8)
|
||||
info.converged < levels && throw(error("Not enough convergence"))
|
||||
if H.hermitian evals = real.(evals) end
|
||||
if H.mode == gpu_cutensor # to avoid possible GPU memory leak
|
||||
CUDA.reclaim()
|
||||
GC.gc(true)
|
||||
end
|
||||
return evals, evecs, info
|
||||
end
|
||||
|
|
@ -1,2 +0,0 @@
|
|||
[TensorOperations]
|
||||
precompile_workload = true
|
||||
|
|
@ -1,7 +0,0 @@
|
|||
[deps]
|
||||
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
|
||||
KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
|
||||
NVTX = "5da4648a-3479-48b8-97b9-01cb529c0a1f"
|
||||
Preferences = "21216c6a-2e73-6563-6e65-726566657250"
|
||||
TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
|
||||
cuTENSOR = "011b41b2-24ef-40a8-b3eb-fa098493e9e1"
|
||||
27
README.md
27
README.md
|
|
@ -1,27 +0,0 @@
|
|||
# DVR-jl
|
||||
|
||||
Solves the quantum $n$-body problem in finite volume (lattice) with periodic boundary conditions. Uses discrete variable representation (DVR) with optional support for complex scaling to study resonances. All details can be found in [H. Yu, N. Yapa, and S. König, Complex scaling in finite volume, Phys. Rev. C 109, 014316 (2024)](https://doi.org/10.1103/PhysRevC.109.014316).
|
||||
|
||||
Written in Julia with optional CUDA GPU acceleration (experimental).
|
||||
|
||||
## Installation
|
||||
|
||||
Make sure you have Julia installed. Required packages can be installed with a single command:
|
||||
```bash
|
||||
julia --project=. -e 'import Pkg; Pkg.instantiate()'
|
||||
```
|
||||
|
||||
## Usage
|
||||
|
||||
See `calculations/3b_bound.jl` for an example on a 3-body bound state.
|
||||
See `calculations/3b_res_from_paper.jl` for an example of a 3-body resonance via complex scaling.
|
||||
|
||||
## Planned features
|
||||
|
||||
- [ ] Spin and isospin degrees of freedom for nuclear calculations
|
||||
- [ ] Multi-node HPC support
|
||||
- [ ] Parity and cubic symmetries ($S_4$)
|
||||
|
||||
## Acknowledgments
|
||||
|
||||
The author gratefully acknowledges the guidance from Sebastian König.
|
||||
14
benchmark.jl
14
benchmark.jl
|
|
@ -1,21 +1,21 @@
|
|||
include("Hamiltonian.jl")
|
||||
using CUDA
|
||||
|
||||
GPU_mode = !("CPU" in ARGS) && CUDA.functional() && CUDA.has_cuda() && CUDA.has_cuda_gpu()
|
||||
|
||||
println("Running with ",Threads.nthreads()," thread(s)")
|
||||
|
||||
if GPU_mode
|
||||
mode=gpu_cutensor
|
||||
include("GPU.jl")
|
||||
println("Available GPUs:")
|
||||
print(" ")
|
||||
println.(name.(devices()))
|
||||
else
|
||||
mode=cpu_tensor
|
||||
include("CPU.jl")
|
||||
end
|
||||
|
||||
T=Float32
|
||||
|
||||
function V_test(r2)
|
||||
function V_test(r2::T)::T
|
||||
return -4*exp(-r2/4)
|
||||
end
|
||||
|
||||
|
|
@ -27,11 +27,11 @@ end
|
|||
|
||||
N=10
|
||||
n_image=1
|
||||
μ=0.5
|
||||
|
||||
for L in 5.0:14.0
|
||||
for L::T in 5.0:14.0
|
||||
println("Constructing H operator...")
|
||||
s=system{T}(3,3,N,L)
|
||||
@time H=Hamiltonian{T}(s,V_test,0,n_image,mode)
|
||||
@time H=HOperator{T}(V_test,3,3,N,L,convert(T,μ),n_image)
|
||||
println("Applying H 1000 times...")
|
||||
if GPU_mode
|
||||
v=CUDA.rand(Complex{T},vectorDims(H)...)
|
||||
|
|
|
|||
|
|
@ -1,19 +0,0 @@
|
|||
include("../Hamiltonian.jl")
|
||||
mode = cpu_tensor
|
||||
T = Float32
|
||||
|
||||
V_gauss(r2) = -2 * exp(-r2 / 4)
|
||||
|
||||
d = 3
|
||||
n = 3
|
||||
N = 20
|
||||
L = 15
|
||||
n_imag = 1
|
||||
ϕ = 0
|
||||
|
||||
s = system{T}(d, n, N, L)
|
||||
H = Hamiltonian{T}(s, V_gauss, ϕ, n_imag, mode)
|
||||
@time evals, _, info = eig(H, 5)
|
||||
|
||||
print(info.numops, " operations")
|
||||
display(evals)
|
||||
|
|
@ -1,36 +0,0 @@
|
|||
# 10.1007/s00601-020-01550-8
|
||||
# Fig. 7
|
||||
# E_R = 4.18(8)
|
||||
|
||||
#./En.run -d 3 -n 3 -N 16 -c pot=v_shifted_gauss,v0=2.0,r=1.5,a=3.0 -c n_eig=20 -c which=li -c tol=1e-6 -L 16 -c phi=0.3 -v
|
||||
|
||||
include("../Hamiltonian.jl")
|
||||
mode = cpu_tensor
|
||||
T = Float32 # single-precision mode
|
||||
|
||||
using Plots
|
||||
|
||||
V_gauss(r2) =
|
||||
2 * exp(-((sqrt(r2) - 3) / 1.5) ^ 2)
|
||||
|
||||
d = 3
|
||||
n = 3
|
||||
N = 16
|
||||
L = 16
|
||||
n_imag = 0
|
||||
|
||||
for ϕ::T in 0.2:0.05:0.4
|
||||
s = system{T}(d, n, N, L)
|
||||
H = Hamiltonian{T}(s, V_gauss, ϕ, n_imag, mode)
|
||||
@time evals, _, info = eig(H, 20)
|
||||
|
||||
print(info.numops, " operations")
|
||||
display(evals)
|
||||
|
||||
scatter(real.(evals), imag.(evals); legend=false)
|
||||
xlabel!("Re E")
|
||||
ylabel!("Im E")
|
||||
xlims!(0, 6)
|
||||
ylims!(-0.6, 0)
|
||||
savefig("temp/phi$(Int(round(ϕ * 100))).png")
|
||||
end
|
||||
|
|
@ -1,24 +0,0 @@
|
|||
include("../Hamiltonian.jl")
|
||||
mode = cpu_tensor
|
||||
T = Float32 # single-precision mode
|
||||
|
||||
V_gauss(r2) =
|
||||
-10 * exp(-(sqrt(r2)) ^ 2)
|
||||
|
||||
d = 3
|
||||
n = 2
|
||||
N = 96
|
||||
ϕ = pi/6
|
||||
n_imag = 1
|
||||
|
||||
open("ComplexScaling-FV-P-res.dat", "w") do f
|
||||
for L = range(20, 35, length=16)
|
||||
println("Calculating L=", L)
|
||||
s = system{T}(d, n, N, L)
|
||||
H = Hamiltonian{T}(s, V_gauss, ϕ, n_imag, mode)
|
||||
@time evals, _, info = eig(H, 40)
|
||||
|
||||
dataline = vcat([L], hcat(real.(evals), imag.(evals))'[:])
|
||||
println(f, join(dataline, '\t'))
|
||||
end
|
||||
end
|
||||
|
|
@ -1,24 +0,0 @@
|
|||
include("../Hamiltonian.jl")
|
||||
mode = cpu_tensor
|
||||
T = Float32 # single-precision mode
|
||||
|
||||
V_gauss(r2) =
|
||||
-10 * exp(-(sqrt(r2)) ^ 2)
|
||||
|
||||
d = 3
|
||||
n = 2
|
||||
N = 30
|
||||
L = 6
|
||||
n_imag = 1
|
||||
|
||||
open("ComplexScaling-FV-S-bound-phi.dat", "w") do f
|
||||
for ϕ = range(0.0, 0.5, length=11)
|
||||
println("Calculating ϕ=", ϕ)
|
||||
s = system{T}(d, n, N, L)
|
||||
H = Hamiltonian{T}(s, V_gauss, ϕ, n_imag, mode)
|
||||
@time evals, _, info = eig(H, 10, resonances = false)
|
||||
|
||||
dataline = vcat([ϕ], hcat(real.(evals), imag.(evals))'[:])
|
||||
println(f, join(dataline, '\t'))
|
||||
end
|
||||
end
|
||||
|
|
@ -1,24 +0,0 @@
|
|||
include("../Hamiltonian.jl")
|
||||
mode = cpu_tensor
|
||||
T = Float32 # single-precision mode
|
||||
|
||||
V_gauss(r2) =
|
||||
2 * exp(- ((sqrt(r2)-3)/1.5) ^ 2)
|
||||
|
||||
d = 3
|
||||
n = 2
|
||||
N = 96
|
||||
L = 30
|
||||
n_imag = 1
|
||||
|
||||
open("ComplexScaling-FV-S-res-phi.dat", "w") do f
|
||||
for ϕ = range(0.1, 0.6, length=26)
|
||||
println("Calculating ϕ=", ϕ)
|
||||
s = system{T}(d, n, N, L)
|
||||
H = Hamiltonian{T}(s, V_gauss, ϕ, n_imag, mode)
|
||||
@time evals, _, info = eig(H, 40, resonances = true)
|
||||
|
||||
dataline = vcat([ϕ], hcat(real.(evals), imag.(evals))'[:])
|
||||
println(f, join(dataline, '\t'))
|
||||
end
|
||||
end
|
||||
|
|
@ -1,67 +0,0 @@
|
|||
using Plots, Arpack
|
||||
|
||||
include("../helper.jl")
|
||||
include("../Hamiltonian.jl")
|
||||
|
||||
mode = cpu_tensor
|
||||
T = Float32 # single-precision mode
|
||||
|
||||
V_r2(c) = r2 -> c * (-5 * exp(-r2/3) + 2 * exp(-r2/10))
|
||||
|
||||
d = 3
|
||||
n = 2
|
||||
N = 48
|
||||
L = 30
|
||||
ϕ = pi/6
|
||||
n_imag = 1
|
||||
s = system{T}(d, n, N, L)
|
||||
|
||||
train_cs = range(0.78, 0.45, length=5)
|
||||
train_ref = reverse([0.05387926313545913-0.008900278182520881im,
|
||||
0.11254295298924327-0.020515067379548786im,
|
||||
0.16060154707503538-0.03716539208626717im,
|
||||
0.19741353362674618-0.05994519982799412im,
|
||||
0.2219100763497223-0.08959449893439568im])
|
||||
|
||||
extrapolate_cs = range(0.38, 0.22, length=5)
|
||||
extrapolate_ref = reverse([0.23165109150003316-0.12052751440975719im,
|
||||
0.23190549514995962-0.1406687118589838im,
|
||||
0.22763660218046278-0.1626190970863793im,
|
||||
0.21807104244164865-0.18635600686249373im,
|
||||
0.2020979906072586-0.21180157628258728im])
|
||||
|
||||
training_E = ComplexF64[]
|
||||
training_vec = Array[]
|
||||
exact_E = ComplexF64[]
|
||||
extrapolated_E = ComplexF64[]
|
||||
|
||||
for c in train_cs
|
||||
println("Training c=", c)
|
||||
H = Hamiltonian{T}(s, V_r2(c), ϕ, n_imag, mode)
|
||||
@time evals, evecs, info = eig(H, 20, resonances = true)
|
||||
i = nearestIndex(evals, pop!(train_ref))
|
||||
push!(training_E, evals[i])
|
||||
push!(training_vec, evecs[i])
|
||||
end
|
||||
|
||||
N_EC = [sum(x .* y) for (x, y) in Iterators.product(training_vec, training_vec)]
|
||||
|
||||
for c in extrapolate_cs
|
||||
println("Extrapolating c=", c)
|
||||
H = Hamiltonian{T}(s, V_r2(c), ϕ, n_imag, mode)
|
||||
@time evals, _, info = eig(H, 40, resonances = true)
|
||||
nearestE = nearest(evals, pop!(extrapolate_ref))
|
||||
push!(exact_E, nearestE)
|
||||
|
||||
# EC extrapolation
|
||||
H_training_vec = H.(training_vec)
|
||||
H_EC = [sum(x .* y) for (x, y) in Iterators.product(training_vec, H_training_vec)]
|
||||
|
||||
evals = eigvals(H_EC, N_EC)
|
||||
push!(extrapolated_E, nearestE)
|
||||
end
|
||||
|
||||
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")
|
||||
savefig("temp/EC-R2R-S.pdf")
|
||||
59
common.jl
59
common.jl
|
|
@ -1,64 +1,53 @@
|
|||
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
|
||||
|
||||
system{T}(d::Int, n::Int, N::Int, L::Real, μ::Real=0.5) where {T<:Float} = new{T}(d, n, N, convert(T, L), convert(T, μ))
|
||||
end
|
||||
|
||||
norm_square(x::Array{Int})::Int = sum(x .* x)
|
||||
|
||||
"Eq (46): Partial derivative matrix element for 1 degree of freedom"
|
||||
function ∂_1DOF(s::system{T}, k::Int, l::Int)::Complex{T} where {T<:Float}
|
||||
function ∂_1DOF(L::T, N::Int, k::Int, l::Int)::Complex{T} where {T<:Float}
|
||||
if k == l
|
||||
return -im * (π / s.L)
|
||||
return -im * (π / L)
|
||||
else
|
||||
return (π / s.L) * (-1)^(k - l) * exp(-im * π * (k - l) / s.N) / sin(π * (k - l) / s.N)
|
||||
return (π / L) * (-1)^(k - l) * exp(-im * π * (k - l) / N) / sin(π * (k - l) / N)
|
||||
end
|
||||
end
|
||||
|
||||
"Which index (dimension of the multidimensional array) corresponds to spatial dimension 'dim' and particle 'p'?"
|
||||
which_index(s::system, dim::Int, p::Int)::Int = (dim - 1) * (s.n - 1) + p
|
||||
which_index(n::Int, dim::Int, p::Int)::Int = (dim - 1) * (n - 1) + p
|
||||
|
||||
"Δk (distance in terms of lattice paramter) between two particles along the given dimension"
|
||||
function get_Δk(s::system, i::CartesianIndex, dim::Int, p1::Int, p2::Int)::Int
|
||||
function get_Δk(n::Int, N::Int, i::CartesianIndex, dim::Int, p1::Int, p2::Int)::Int
|
||||
if p1 == p2
|
||||
return 0
|
||||
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
|
||||
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
|
||||
else
|
||||
return i[which_index(s, dim, p1)] - i[which_index(s, dim, p2)]
|
||||
return i[which_index(n, dim, p1)] - i[which_index(n, dim, p2)]
|
||||
end
|
||||
end
|
||||
|
||||
"Calculate diagonal elements of the V matrix"
|
||||
function calculate_Vs(s::system{T}, V_twobody::Function, ϕ::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))...)
|
||||
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))...)
|
||||
Threads.@threads for i in CartesianIndices(Vs)
|
||||
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
|
||||
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
|
||||
else
|
||||
min_Δk[dim] = Δk
|
||||
end
|
||||
end
|
||||
for image in images
|
||||
Δk² = norm_square(min_Δk .- (s.N .* image))
|
||||
Δk² = norm_square(min_Δk .- (N .* image))
|
||||
Vs[i] += V_twobody(Δk² * coeff²)
|
||||
end
|
||||
end
|
||||
|
|
|
|||
|
|
@ -6,9 +6,8 @@
|
|||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# prerequisite packages: KrylovKit, TensorOperations, LinearAlgebra, CUDA#tb/cutensor, Plots\n",
|
||||
"include(\"Hamiltonian.jl\")\n",
|
||||
"mode = cpu_tensor # using CPU mode\n",
|
||||
"# prerequisite packages: KrylovKit, TensorOperations, LinearAlgebra, CUDA#tb/cutensor\n",
|
||||
"include(\"CPU.jl\") # using CPU mode\n",
|
||||
"T = Float32 # single-precision mode"
|
||||
]
|
||||
},
|
||||
|
|
@ -24,12 +23,12 @@
|
|||
"d = 3\n",
|
||||
"n = 3\n",
|
||||
"N = 6\n",
|
||||
"L = 12\n",
|
||||
"ϕ = 0.0\n",
|
||||
"L::T = 12\n",
|
||||
"ϕ::T = 0.0\n",
|
||||
"mu::T = 0.5\n",
|
||||
"n_imag = 1\n",
|
||||
"\n",
|
||||
"s = system{T}(d, n, N, L)\n",
|
||||
"H = Hamiltonian{T}(s, V_gauss, ϕ, n_imag, mode)\n",
|
||||
"H = HOperator{T}(V_gauss, d, n, N, L, ϕ, mu, n_imag)\n",
|
||||
"@time evals, evecs, info = eig(H, 5)\n",
|
||||
"print(info.numops, \" operations : \")\n",
|
||||
"println(evals)"
|
||||
|
|
@ -49,12 +48,12 @@
|
|||
"d = 3\n",
|
||||
"n = 2\n",
|
||||
"N = 32\n",
|
||||
"L = 16\n",
|
||||
"ϕ = 0.5\n",
|
||||
"L::T = 16\n",
|
||||
"ϕ::T = 0.5\n",
|
||||
"mu::T = 0.5\n",
|
||||
"n_imag = 0\n",
|
||||
"\n",
|
||||
"s = system{T}(d, n, N, L)\n",
|
||||
"H = Hamiltonian{T}(s, V_gauss, ϕ, n_imag, mode)\n",
|
||||
"H = HOperator{T}(V_gauss, d, n, N, L, ϕ, mu, n_imag)\n",
|
||||
"@time evals, evecs, info = eig(H, 20)\n",
|
||||
"print(info.numops, \" operations : \")\n",
|
||||
"print(evals)\n",
|
||||
|
|
|
|||
|
|
@ -1,5 +0,0 @@
|
|||
"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)]
|
||||
|
|
@ -1,71 +0,0 @@
|
|||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"include(\"Hamiltonian.jl\")\n",
|
||||
"\n",
|
||||
"println(\"Running with \",Threads.nthreads(),\" thread(s)\")\n",
|
||||
"println(\"Available GPUs:\")\n",
|
||||
"println.(name.(devices()))\n",
|
||||
"\n",
|
||||
"T=Float32\n",
|
||||
"\n",
|
||||
"function V_test(r2)\n",
|
||||
" return -4*exp(-r2/4)\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"function test(mode)\n",
|
||||
" for (n,N) in [(2,16),(3,8)]\n",
|
||||
" println(\"\\n$n-body system with N=$N\")\n",
|
||||
" n_image=0\n",
|
||||
" for L::T in 5.0:9.0\n",
|
||||
" print(\"L=$L\")\n",
|
||||
" s=system{T}(3,n,N,L)\n",
|
||||
" H=Hamiltonian{T}(s,V_test,0.0,n_image,mode)\n",
|
||||
" evals,_,_ = eig(H,5)\n",
|
||||
" println(real.(evals))\n",
|
||||
" end\n",
|
||||
" end\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"test(cpu_tensor)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"test(gpu_cutensor)"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Julia 1.8.5",
|
||||
"language": "julia",
|
||||
"name": "julia-1.8"
|
||||
},
|
||||
"language_info": {
|
||||
"file_extension": ".jl",
|
||||
"mimetype": "application/julia",
|
||||
"name": "julia",
|
||||
"version": "1.8.5"
|
||||
},
|
||||
"orig_nbformat": 4
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 2
|
||||
}
|
||||
Loading…
Reference in New Issue