Added 1D system for testing

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

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