Added 1D system for testing

.gitignore

@@ -1,15 +1,3 @@
-# VSCode
-.vscode/
-
-# HPC scripts and logs
-hpc/
-
-# Calculation outputs
-*.dat
-*.csv
-*.hdf5
-*.out
-
 # Temporary and scratch files
 temp/
 scratch/

Hamiltonian.jl

@@ -1,40 +1,43 @@
 include("common.jl")
-using TensorOperations, KrylovKit, LinearAlgebra, CUDA, cuTENSOR, NVTX
+using TensorOperations, KrylovKit, LinearAlgebra, CUDA, CUDA.CUTENSOR
 
 @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}}}
+    d::Int
+    n::Int
+    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}}
     hermitian::Bool
     mode::Hamiltonian_backend
-
-    function Hamiltonian{T}(s::system{T}, V_twobody::Function, ϕ::Real, n_image::Int, mode::Hamiltonian_backend) where {T<:Float}
+    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 = -s.N÷2:s.N÷2-1
-        Vs = calculate_Vs(s, V_twobody, convert(T, ϕ), n_image)
+        k = -N÷2:N÷2-1
+        Vs = calculate_Vs(V_twobody, d, n, N, L, ϕ, 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
+        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_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)
+            return new{T}(d, n, N, L, μ, nothing, K_diag, K_mixed, CuArray(Vs), hermitian, mode)
         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, 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)::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))...)
+vectorDims(H::Hamiltonian)::Dims = tuple(fill(H.N, H.d * (H.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}
@@ -42,13 +45,14 @@ function LinearAlgebra.mul!(out::Array{Complex{T}}, H::Hamiltonian{T}, v::Array{
     # 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
+    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.s, dim, coord1)
-                i2 = which_index(H.s, dim, coord2)
+                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)
@@ -58,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.K_partial, H.K_partial, v), (nconList_1, nconList_2, nconList_v))
-                out = axpy!(1, v_new, out)
+                v_new = @ncon((H.∂1, H.∂1, v), (nconList_1, nconList_2, nconList_v))
+                out = axpy!(coeff, v_new, out)
             end
         end
     end
@@ -68,8 +72,8 @@ 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)
+    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
 
@@ -81,15 +85,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.s.n - 1
-    inds_template = ('a' - 1) .+ collect(1:H.s.d*(coords))
+    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.s.d
+    for dim = 1:H.d
         for coord1 = 1:coords
             for coord2 = 1:coord1
-                i1 = which_index(H.s, dim, coord1)
-                i2 = which_index(H.s, dim, coord2)
+                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"
@@ -134,12 +138,8 @@ function eig(H::Hamiltonian{T}, levels::Int; resonances = !H.hermitian)::Tuple{V
         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"))
+    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
     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

LocalPreferences.toml

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

Project.toml

@@ -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"

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.

V_dump.ipynb

@@ -0,0 +1,51 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# ./En.run -d 3 -n 3 -e 5 -c eps=0 -c pot=v_gauss,v0=-4,r=2 -N 6 -L 5:14 -c n_imag=1\n",
+    "\n",
+    "include(\"common.jl\")\n",
+    "\n",
+    "T=Float32\n",
+    "\n",
+    "function V_test(r2::T)::T\n",
+    "    return -4*exp(-r2/4)\n",
+    "end\n",
+    "\n",
+    "N=6\n",
+    "L::T=14.0\n",
+    "n_image=0\n",
+    "\n",
+    "V=calculate_Vs(V_test, 3, 3, N, L, n_image)\n",
+    "\n",
+    "outfile = \"temp/V_vals.dat\"\n",
+    "\n",
+    "open(outfile, \"w\") do f\n",
+    "    for i in V\n",
+    "        write(f, i)\n",
+    "    end\n",
+    "end"
+   ]
+  }
+ ],
+ "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
+}

V_verify.ipynb

@@ -0,0 +1,87 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from itertools import chain\n",
+    "import numpy as np\n",
+    "import math\n",
+    "from scipy import sparse\n",
+    "from scipy.sparse.linalg import eigsh\n",
+    "\n",
+    "#sample potential\n",
+    "def V_gauss(dr_sqr):\n",
+    "    return -4*math.exp(-dr_sqr/4)\n",
+    "\n",
+    "n=3 # no of particles\n",
+    "L=14\n",
+    "N=6 # no of lattice points\n",
+    "mu=1/2 # reduced mass\n",
+    "\n",
+    "DOF=(n-1)*3 # degrees of freedom after excluding CM\n",
+    "\n",
+    "s=np.arange(N**DOF) # matrix index\n",
+    "\n",
+    "# k index for each particle and each dimension\n",
+    "k=np.empty((n-1,3),dtype=np.dtype)\n",
+    "for dof in range(DOF):\n",
+    "    k[dof//3,dof%3]=(s%N**(DOF-dof))//N**(DOF-1-dof)-N//2\n",
+    "\n",
+    "x=k*(L/N) # x coordinate from k index\n",
+    "\n",
+    "# adding up all non-local interactions\n",
+    "V_local=np.zeros(N**DOF)\n",
+    "\n",
+    "# 2-body interactions\n",
+    "V2=np.vectorize(V_gauss)\n",
+    "dxs=chain((x[i,:] for i in range(n-1)), (x[i,:]-x[j,:] for (i,j) in np.ndindex((n-1,n-1)) if i<j)) # with last particle + with each other\n",
+    "for dx in dxs:\n",
+    "    dr_sqr=np.sum(dx*dx)\n",
+    "    V_local+=V2(dr_sqr)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "filename = \"temp/V_vals.dat\"\n",
+    "\n",
+    "with open(filename, 'br') as f:\n",
+    "    buffer = f.read()\n",
+    "\n",
+    "V_julia = np.frombuffer(buffer, dtype=np.float32)\n",
+    "\n",
+    "abs_diff = np.abs(V_local-V_julia)\n",
+    "s = np.mean(abs_diff)\n",
+    "print(s)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "base",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.13"
+  },
+  "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

benchmark.jl

@@ -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=Hamiltonian{T}(V_test,3,3,N,L,convert(T,0),convert(T,μ),n_image,mode)
     println("Applying H 1000 times...")
     if GPU_mode
         v=CUDA.rand(Complex{T},vectorDims(H)...)

calculations/3b_bound.jl

@@ -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)

calculations/3b_res_from_paper.jl

@@ -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

calculations/ComplexScaling-FV-P-res.jl

@@ -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

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

@@ -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

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

@@ -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

calculations/EC-R2R-S.jl

@@ -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")

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

example.ipynb

@@ -24,12 +24,12 @@
     "d = 3\n",
     "n = 3\n",
     "N = 6\n",
-    "L = 12\n",
-    "ϕ = 0.0\n",
+    "L::T = 12\n",
+    "ϕ::T = 0.0\n",
+    "μ::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 = Hamiltonian{T}(V_gauss, d, n, N, L, ϕ, μ, n_imag, mode)\n",
     "@time evals, evecs, info = eig(H, 5)\n",
     "print(info.numops, \" operations : \")\n",
     "println(evals)"
@@ -49,12 +49,12 @@
     "d = 3\n",
     "n = 2\n",
     "N = 32\n",
-    "L = 16\n",
-    "ϕ = 0.5\n",
+    "L::T = 16\n",
+    "ϕ::T = 0.5\n",
+    "μ::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 = Hamiltonian{T}(V_gauss, d, n, N, L, ϕ, μ, n_imag, mode)\n",
     "@time evals, evecs, info = eig(H, 20)\n",
     "print(info.numops, \" operations : \")\n",
     "print(evals)\n",

helper.jl

@@ -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)]

mat_eig.jl

@@ -0,0 +1,42 @@
+include("Hamiltonian.jl")
+
+T=Float32
+
+function V_test(r2)
+    return -4 * exp(-r2 / 4)
+end
+
+d = 1
+n = 3
+N = 8
+L::T = 16
+n_imag = 0
+μ::T = 0.5
+
+H = HOperator{T}(V_test, d, n, N, L, μ, n_imag, cpu_tensor)
+dim = N ^ (d * (n - 1))
+H_mat = zeros(Complex{T}, dim, dim)
+iter = CartesianIndices(vectorDims(H))
+
+open("temp/mat_dump.csv", "w") do f
+    # this can be heavily optimized by getting rid of 'bi' vector
+    for i in 1 : dim
+        bi = zeros(Complex{T}, vectorDims(H)...)
+        bi[iter[i]] = 1
+
+        for j in 1 : dim
+            bj = zeros(Complex{T}, vectorDims(H)...)
+            bj[iter[j]] = 1
+
+            Hbj = similar(bj)
+            Hbj = mul!(Hbj, H, bj)
+
+            H_mat[i, j] = dot(bi, Hbj)
+        end
+    end
+end
+
+evals, _ = eigen(H_mat)
+evals = real.(evals)
+
+print(evals)

tester.ipynb

@@ -0,0 +1,122 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# ./En.run -d 3 -n 3 -e 5 -c pot=v_gauss,v0=-4,r=2 -N 8 -L 4 -c n_imag=0\n",
+    "# Calculating...\n",
+    "# -> N = 8\n",
+    "# -> L = 4\n",
+    "# -> Spectrum = {-6.07632,-3.81486,-3.71969,-3.71968,-3.38263}...\n",
+    "# Done.\n",
+    "# -> Time used = 12s\n",
+    "\n",
+    "include(\"Hamiltonian.jl\")\n",
+    "\n",
+    "T=Float32\n",
+    "\n",
+    "function V_test(r2)\n",
+    "    return -4*exp(-r2/4)\n",
+    "end\n",
+    "\n",
+    "N=8\n",
+    "L::T=4\n",
+    "n_imag=0\n",
+    "\n",
+    "H=Hamiltonian{T}(V_test,3,3,N,L,convert(T,0),convert(T,0.5),n_imag,cpu_tensor)\n",
+    "@time evals,evecs,info=eig(H,5)\n",
+    "print(info.numops,\" operations : \")\n",
+    "println(evals)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# ./En.run -d 3 -n 3 -e 5 -c pot=v_gauss,v0=-4,r=2 -N 6 -L 14 -c n_imag=0\n",
+    "# Calculating...\n",
+    "# -> N = 6\n",
+    "# -> L = 14\n",
+    "# -> Spectrum = {-8.49538,-3.35492,-3.34356,-3.32830,-3.07909}...\n",
+    "# Done.\n",
+    "# -> Time used = 0.229049s\n",
+    "\n",
+    "include(\"Hamiltonian.jl\")\n",
+    "\n",
+    "T=Float32\n",
+    "\n",
+    "function V_test(r2)\n",
+    "    return -4*exp(-r2/4)\n",
+    "end\n",
+    "\n",
+    "N=6\n",
+    "L::T=14\n",
+    "n_imag=0\n",
+    "\n",
+    "H=Hamiltonian{T}(V_test,3,3,N,L,convert(T,0),convert(T,0.5),n_imag,cpu_tensor)\n",
+    "@time evals,evecs,info=eig(H,5)\n",
+    "print(info.numops,\" operations : \")\n",
+    "println(evals)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# ./En.run -d 3 -n 3 -e 5 -c pot=v_gauss,v0=-4,r=2 -N 6 -L 14 -c n_imag=0\n",
+    "# Calculating...\n",
+    "# -> N = 6\n",
+    "# -> L = 14\n",
+    "# -> Spectrum = {-8.49538,-3.35492,-3.34356,-3.32830,-3.07909}...\n",
+    "# Done.\n",
+    "# -> Time used = 0.229049s\n",
+    "\n",
+    "include(\"Hamiltonian.jl\")\n",
+    "tolerance=1e-10\n",
+    "T=Float64\n",
+    "\n",
+    "function V_test(r2)\n",
+    "    return -4*exp(-r2/4)\n",
+    "end\n",
+    "\n",
+    "N=4\n",
+    "L::T=16\n",
+    "n_imag=0\n",
+    "\n",
+    "H=Hamiltonian{T}(V_test,1,3,N,L,convert(T,0),convert(T,0.5),n_imag,cpu_tensor)\n",
+    "evals,evecs,info=eig(H,16)\n",
+    "display(evals)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "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
+}

testing.ipynb

@@ -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
-}