Update README.md

.gitignore

@@ -1,3 +1,9 @@
+# VSCode
+.vscode/
+
+# HPC scripts and logs
+hpc/
+
 # Calculation outputs
 *.dat
 *.csv

Hamiltonian.jl

@@ -6,23 +6,27 @@ using TensorOperations, KrylovKit, LinearAlgebra, CUDA, cuTENSOR, NVTX
 "A Hamiltonian that can be applied to a vector"
 struct Hamiltonian{T}
     s::system{T}
-    K::Union{CuTensor{Complex{T}}, Matrix{Complex{T}}}
-    Vs::Union{Array{Complex{T}}, 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}(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"
-        hermitian = ϕ == 0.0
-        Vs = calculate_Vs(s, V_twobody, convert(T, ϕ), n_image)
         k = -s.N÷2:s.N÷2-1
-        ∂ = ∂_1DOF.(Ref(s), k, k')
-        K = exp(-2im * convert(T, ϕ)) .* (∂ * ∂)    # TODO: Calculate K matrix elements directly
+        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 = CuTensor(CuArray(K), ['a', 'A'])
+            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, Vs, hermitian, mode)
+        return new{T}(s, K_partial, K_diag, K_mixed, Vs, hermitian, mode)
     end
 end
 
@@ -41,22 +45,30 @@ function LinearAlgebra.mul!(out::Array{Complex{T}}, H::Hamiltonian{T}, v::Array{
     coords = H.s.n - 1
     nconList_v_template = -collect(1:H.s.d*(coords))
     for dim = 1:H.s.d
-        for coord = 1:coords
-            i = which_index(H.s, dim, coord)
-            nconList_K = [-i, 1]
-            nconList_v = copy(nconList_v_template)
-            nconList_v[i] = 1
-            v_new = @ncon((H.K, v), (nconList_K, nconList_v))
-            coeff = -1 / (2 * H.s.μs[coord])
-            out = axpy!(coeff, v_new, out)
+        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!(alpha::Number, C::CuTensor, A::CuTensor, B::CuTensor)::CuTensor
-    cuTENSOR.contraction!(alpha, A.data, A.inds, cuTENSOR.CUTENSOR_OP_IDENTITY, B.data, B.inds, cuTENSOR.CUTENSOR_OP_IDENTITY,
+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
@@ -74,16 +86,31 @@ function LinearAlgebra.mul!(out::CuArray{Complex{T}}, H::Hamiltonian{T}, v::CuAr
     v_t = CuTensor(v, copy(inds_template))
     out_t = CuTensor(out, copy(inds_template))
     for dim = 1:H.s.d
-        for coord = 1:coords
-            i = which_index(H.s, dim, coord)
-            @assert v_t.inds == inds_template "v indices permuted"
-            @assert H.K.inds[2] == 'A' "K_diag indices permuted"
-            H.K.inds[1] = 'a' - 1 + i
-            v_t.inds[i] = 'A'
-            #synchronize(ctx)
-            coeff = -1 / (2 * H.s.μs[coord])
-            NVTX.@range "K" out_t = contract_accumulate!(coeff, out_t, H.K, v_t)
-            v_t.inds[i] = 'a' - 1 + i
+        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"

LocalPreferences.toml

@@ -0,0 +1,2 @@
+[TensorOperations]
+precompile_workload = true

Project.toml

@@ -0,0 +1,7 @@
+[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"

README.md

@@ -0,0 +1,27 @@
+# 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.

calculations/3b_bound.jl

@@ -0,0 +1,19 @@
+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)

calculations/3b_res_from_paper.jl

@@ -0,0 +1,36 @@
+# 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

calculations/ComplexScaling-FV-P-res.jl

@@ -0,0 +1,24 @@
+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

calculations/ComplexScaling-FV-S-bound-phi.jl

@@ -0,0 +1,24 @@
+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

calculations/ComplexScaling-FV-S-res-phi.jl

@@ -0,0 +1,24 @@
+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

calculations/EC-R2R-S.jl

@@ -0,0 +1,67 @@
+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")

common.jl

@@ -1,38 +1,18 @@
 Float = Union{Float32,Float64}
 
-norm_square(x) = sum(x .* x)
-reducedMass(m1, m2) = 1 / (1/m1 + 1/m2)
-
 "A few-body system defined by its physical parameters"
 struct system{T}
     d::Int
     n::Int
     N::Int
     L::T
-    μs::Vector{T}
-    invU::Matrix{T}
+    μ::T
 
-    function system{T}(d::Int, n::Int, N::Int, L::Real) where {T<:Float}
-        μs = [1/((coord + 1)^2 * reducedMass(coord, 1)) for coord in 1:(n - 1)]
-
-        # TODO: Optimize
-        U = Matrix{T}(undef, n, n)
-        for i in CartesianIndices(U)
-            if i[1] + 1 == i[2]
-                U[i] = -i[1]
-            elseif i[1] >= i[2]
-                U[i] = 1
-            else
-                U[i] = 0
-            end
-        end
-        U[n, :] .= 1/n
-        invU = inv(U)[:, 1:(n - 1)]
-        
-        return new{T}(d, n, N, convert(T, L), μs, invU)
-    end
+    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}
     if k == l
@@ -42,20 +22,19 @@ function ∂_1DOF(s::system{T}, k::Int, l::Int)::Complex{T} where {T<:Float}
     end
 end
 
-"Which index (dimension of the multidimensional array) corresponds to spatial dimension 'dim' of coordinate 'coord'?"
-which_index(s::system, dim::Int, coord::Int)::Int = (dim - 1) * (s.n - 1) + coord
+"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
 
-"Get the distance to the nearest image of the particle"
-function nearest(s::system, Δk)
-    # TODO: Optimize
-    while true
-        if Δk >= s.N ÷ 2
-            Δk -= s.N
-        elseif Δk < -s.N ÷ 2
-            Δk += s.N
-        else
-            return Δk
-        end            
+"Δ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
+    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
+    else
+        return i[which_index(s, dim, p1)] - i[which_index(s, dim, p2)]
     end
 end
 
@@ -65,17 +44,21 @@ function calculate_Vs(s::system{T}, V_twobody::Function, ϕ::T, n_image::Int)::A
     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)
-        xs = reshape(collect(Tuple(i)), s.n - 1, s.d) .- (s.N ÷ 2 + 1)
-        rs = s.invU * xs
         for p1 in 1:s.n
-            for p2 in 1:(p1 - 1)
-                Δk = Array{T}(undef, s.d)
+            for p2 in (p1 + 1):s.n
+                min_Δk = Array{Int}(undef, s.d)
                 for dim in 1:s.d
-                    Δk_temp = rs[p1, dim] - rs[p2, dim]
-                    Δk[dim] = nearest(s, Δk_temp)
+                    Δ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(Δk .- (s.N .* image))
+                    Δk² = norm_square(min_Δk .- (s.N .* image))
                     Vs[i] += V_twobody(Δk² * coeff²)
                 end
             end

helper.jl

@@ -0,0 +1,5 @@
+"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)]

test_jacobi.jl

@@ -1,21 +0,0 @@
-include("Hamiltonian.jl")
-
-T=Float32
-
-function V_test(r2)
-    return -4*exp(-r2/4)
-end
-
-for (n,N) in [(2,16), (3,8)]
-    println("\n$n-body system with N=$N")
-    n_image=0
-    for L::T in 5.0:9.0
-        print("L=$L: ")
-        s=system{T}(3,n,N,L)
-        print("Constructing H...")
-        H=Hamiltonian{T}(s,V_test,0.0,n_image,cpu_tensor)
-        print("Diagonalizing...")
-        evals,_,_ = eig(H,5)
-        println(real.(evals))
-    end
-end