Merge branch 'master' into debugging

common.jl

@@ -1,5 +1,7 @@
 Float = Union{Float32,Float64}
 
+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}
     if k == l
@@ -9,35 +11,43 @@ function ∂_1DOF(L::T, N::Int, k::Int, l::Int)::Complex{T} where {T<:Float}
     end
 end
 
-"Which index (dimension of the multidimensional array) corresponds to this dimension and coordinate?"
-which_index(n::Int, dim::Int, coord::Int)::Int = (dim - 1) * (n - 1) + coord
+"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
 
-"k value of the given degree of freedom at the corresponding index, with coord=0 always returning 0"
-get_k(n::Int, N::Int, i::CartesianIndex, dim::Int, coord::Int)::Int =
-    coord == 0 ? 0 : i[which_index(n, dim, coord)] - N ÷ 2 - 1
-
-"k value of the DOF at the specified cubic image"
-get_shifted_k(n::Int, N::Int, i::CartesianIndex, dim::Int, coord::Int, image::Vector{Int})::Int =
-    get_k(n, N, i, dim, coord) + N * image[dim]
-
-"Difference of k values between two particles"
-get_Δk(n::Int, N::Int, i::CartesianIndex, dim::Int, coord1::Int, coord2::Int, image::Vector{Int})::Int =
-    get_k(n, N, i, dim, coord1) - get_shifted_k(n, N, i, dim, coord2, image)
+"Δ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
+    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
+    else
+        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(V_twobody::Function, d::Int, n::Int, N::Int, L::T, n_image::Int)::Array{T} where {T<:Float}
     L²_over_N² = (L / N)^2
-    coords = n - 1
     images = collect.(Iterators.product(fill(-n_image:n_image, d)...)) # TODO: Learn how to use tuples instead of vectors
-    Vs = zeros(T, fill(N, d * coords)...)
-    for image in images
-        Threads.@threads for i in CartesianIndices(Vs)
-            for coord1 in 1:coords
-                for coord2 in 0:coord1-1
-                    Δk² = 0
-                    for dim in 1:d
-                        Δk² += get_Δk(n, N, i, dim, coord1, coord2, image)^2
+    Vs = zeros(T, fill(N, d * (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
+                    else
+                        min_Δk[dim] = Δk
                     end
+                end
+                for image in images
+                    Δk² = norm_square(min_Δk .- (N .* image))
                     Vs[i] += V_twobody(Δk² * L²_over_N²)
                 end
             end

example.ipynb

@@ -0,0 +1,47 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# prerequisites: KrylovKit, TensorOperations, LinearAlgebra, CUDA#tb/cutensor (for GPU mode)\n",
+    "\n",
+    "include(\"CPU.jl\") # using CPU mode\n",
+    "T = Float32\n",
+    "\n",
+    "V_test(r2::T)::T =\n",
+    "    -4 * exp(-r2 / 4)\n",
+    "\n",
+    "d = 3\n",
+    "n = 3\n",
+    "N = 6\n",
+    "L::T = 12.0\n",
+    "mu::T = 0.5\n",
+    "n_imag = 1\n",
+    "\n",
+    "H = HOperator{T}(V_test, 3, 3, N, L, 0.5f0, 1)\n",
+    "@time evals, evecs, info = eig(H, 5)\n",
+    "print(info.numops, \" operations : \")\n",
+    "println(evals)"
+   ]
+  }
+ ],
+ "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
+}

tester.jl

@@ -1,43 +0,0 @@
-# ./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
-
-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
-    include("GPU.jl")
-    println("Available GPUs:")
-    print("    ")
-    println.(name.(devices()))
-else
-    include("CPU.jl")
-end
-
-T=Float32
-
-function V_zero(r2::T)::T
-    return 0.0
-end
-
-function V_test(r2::T)::T
-    return -4*exp(-r2/4)
-end
-
-N=6
-n_image=1
-μ=0.5
-
-levels=5
-
-for L::T in 5.0:14.0
-    H=HOperator{T}(V_test,3,3,N,L,convert(T,μ),n_image)
-    if GPU_mode
-        CUDA.@time evals,evecs,info=eig(H,levels)
-    else
-        @time evals,evecs,info=eig(H,levels)
-    end
-    print(info.numops," operations : ")
-    println(evals)
-end