Naive n-body extension

Hamiltonian.jl

@@ -9,6 +9,9 @@ struct Hamiltonian{T}
     K_partial::Matrix{Complex{T}}
     K_diag::Union{CuTensor{Complex{T}},Nothing}
     K_mixed::Union{CuTensor{Complex{T}},Nothing}
+    K_partial_1::Union{Tuple,Nothing}
+    K_partial_2::Union{Tuple,Nothing}
+    K_partial_c::Union{Tuple,Nothing}
     Vs::Union{Array{Complex{T}},CuArray{Complex{T}}}
     hermitian::Bool
     mode::Hamiltonian_backend
@@ -16,25 +19,27 @@ struct Hamiltonian{T}
     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
+        K_partial_1, K_partial_2, K_partial_c = sym_reduce(s, K_partial)
+        Vs = calculate_Vs(s, V_twobody, convert(T, ϕ), n_image)
         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)
+        else
+            K_diag = nothing
+            K_mixed = nothing    
         end
-        return new{T}(s, K_partial, K_diag, K_mixed, Vs, hermitian, mode)
+        return new{T}(s, K_partial, K_diag, K_mixed, K_partial_1, K_partial_2, K_partial_c, 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.s.N^(H.s.d * (H.s.n - 2)) * length(H.s.unique_i) : 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(length(H.s.unique_i), fill(H.s.N, H.s.d * (H.s.n - 2))...)
 
 "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,10 +47,9 @@ 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))
+    nconList_v_template = -collect(1:(H.s.d * (H.s.n - 2) + 1))
     for dim = 1:H.s.d
-        for coord1 = 1:coords
+        for coord1 = 1:(H.s.n - 1)
             for coord2 = 1:coord1
                 i1 = which_index(H.s, dim, coord1)
                 i2 = which_index(H.s, dim, coord2)
@@ -58,7 +62,17 @@ 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))
+
+                if coord1 == 1 && coord2 == 1
+                    tensor1 = H.K_partial_1[dim]
+                    tensor2 = H.K_partial_2[dim]
+                else
+                    tensor1 = coord1 == 1 ? H.K_partial_c[dim] : H.K_partial
+                    tensor2 = coord2 == 1 ? H.K_partial_c[dim] : H.K_partial
+                end
+
+                v_new = @ncon((tensor1, tensor2, v), (nconList_1, nconList_2, nconList_v))
+                
                 out = axpy!(1, v_new, out)
             end
         end
@@ -134,7 +148,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 * 4)
+    KrylovKit_hermitian = H.hermitian && H.s.sym == all
+    evals, evecs, info = eigsolve(H, x₀, levels, resonances ? :LI : :SR; ishermitian = KrylovKit_hermitian, tol = tolerance, krylovdim = levels * 4)
     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

common.jl

@@ -1,4 +1,7 @@
+include("irrep.jl")
+
 Float = Union{Float32,Float64}
+@enum rep all A1
 
 "A few-body system defined by its physical parameters"
 struct system{T}
@@ -8,7 +11,23 @@ struct system{T}
     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, μ))
+    sym::rep
+    unique_i::Array{Int}
+    unique_point::Array{Int}
+    multiplicity::Array{Int}
+    labels::Array{Int}
+
+    function system{T}(d::Int, n::Int, N::Int, L::Real, μ::Real=0.5, sym::rep=all) where {T<:Float}
+        @assert d == 3 "Only supports 3D"
+        if sym == all
+            unique_i, unique_point, multiplicity, labels = calculate_all_data(N)
+        elseif sym == A1
+            unique_i, unique_point, multiplicity, labels = calculate_A1_data(N)
+        else
+            throw(ArgumentError("Symmetry not yet implemented"))
+        end
+        return new{T}(d, n, N, convert(T, L), convert(T, μ), sym, unique_i, unique_point, multiplicity, labels)
+    end
 end
 
 norm_square(x::Array{Int})::Int = sum(x .* x)
@@ -23,18 +42,27 @@ function ∂_1DOF(s::system{T}, k::Int, l::Int)::Complex{T} where {T<:Float}
 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(s::system, dim::Int, p::Int)::Int = p == 1 ? 1 : (dim - 1) * (s.n - 2) + p
+
+"Δk (distance in terms of lattice paramter) between two particles along the given dimension"
+function get_k(s::system, i::CartesianIndex, dim::Int, p::Int)::Int
+    if p == 1
+        s.unique_point[i[1], dim]
+    else
+        return i[which_index(s, dim, p)] - s.N ÷ 2 - 1
+    end
+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)
+        return -get_k(s, i, dim, p2)
     elseif p2 == s.n
-        return i[which_index(s, dim, p1)] - s.N ÷ 2 - 1
+        return get_k(s, i, dim, p1)
     else
-        return i[which_index(s, dim, p1)] - i[which_index(s, dim, p2)]
+        return get_k(s, i, dim, p1) - get_k(s, i, dim, p2)
     end
 end
 
@@ -42,7 +70,7 @@ end
 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))...)
+    Vs = zeros(Complex{T}, length(s.unique_i), fill(s.N, s.d * (s.n - 2))...)
     Threads.@threads for i in CartesianIndices(Vs)
         for p1 in 1:s.n
             for p2 in (p1 + 1):s.n

irrep.jl

@@ -0,0 +1,77 @@
+using DelimitedFiles, LinearAlgebra
+
+function calculate_all_data(N::Int)
+    ks = -N÷2:N÷2-1
+    lattice = hcat((collect.(Iterators.product(ks, ks, ks)))...)
+    labels = reshape(collect(1:N^3), (N, N, N))
+
+    unique_i = collect(1:N^3)
+    multiplicity = fill(1, length(unique_i))
+    unique_point = transpose(lattice)
+
+    return unique_i, unique_point, multiplicity, labels
+end
+
+function calculate_A1_data(N::Int)
+    rotations = readdlm("rotations.mat", ',', Int, '\n')
+    rotations = reshape(rotations, (24, 3, 3))
+    
+    ks = -N÷2:N÷2-1
+    lattice = hcat((collect.(Iterators.product(ks, ks, ks)))...)
+    labels = reshape(collect(1:N^3), (N, N, N))
+    
+    for r in 1:24
+        rotated_lattice = Matrix(rotations[r, :, :]) * lattice
+        for index in 1:N^3
+            rotated_lattice_point = rotated_lattice[:, index]
+            (i, j, k) = mod1.(rotated_lattice_point .+ (N÷2 + 1), N)
+            old_label = max(labels[index], labels[i, j, k])
+            new_label = min(labels[index], labels[i, j, k])
+            if old_label != new_label
+                for o in findall(isequal(old_label), labels)
+                    labels[o] = new_label
+                end
+            end
+        end
+    end
+
+    unique_i = unique(labels)
+    multiplicity = [count(labels.==i) for i in unique_i]
+    unique_point = transpose(lattice[:, unique_i])
+
+    return unique_i, unique_point, multiplicity, labels
+end
+
+function sym_reduce(s, K_partial)
+    I = one(K_partial)
+
+    K_partial_x1 = kron(kron(K_partial, I), I)
+    K_partial_y1 = kron(kron(I, K_partial), I)
+    K_partial_z1 = kron(kron(I, I), K_partial)
+
+    K_partial_x1 = K_partial_x1[s.unique_i, :]
+    K_partial_y1 = K_partial_y1[s.unique_i, :]
+    K_partial_z1 = K_partial_z1[s.unique_i, :]
+
+    K_partial_x2 = kron(kron(K_partial, I), I)
+    K_partial_y2 = kron(kron(I, K_partial), I)
+    K_partial_z2 = kron(kron(I, I), K_partial)
+
+    for (i, label) in enumerate(s.labels)
+        if i != label
+            K_partial_x2[:, label] .+= K_partial_x2[:, i]
+            K_partial_y2[:, label] .+= K_partial_y2[:, i]
+            K_partial_z2[:, label] .+= K_partial_z2[:, i]
+        end
+    end
+
+    K_partial_x2 = K_partial_x2[:, s.unique_i]
+    K_partial_y2 = K_partial_y2[:, s.unique_i]
+    K_partial_z2 = K_partial_z2[:, s.unique_i]
+
+    K_partial_xc = K_partial_x2[s.unique_i, :]
+    K_partial_yc = K_partial_y2[s.unique_i, :]
+    K_partial_zc = K_partial_z2[s.unique_i, :]
+
+    return (K_partial_x1, K_partial_y1, K_partial_z1), (K_partial_x2, K_partial_y2, K_partial_z2), (K_partial_xc, K_partial_yc, K_partial_zc)
+end

rotations.mat

@@ -0,0 +1,24 @@
+1,0,0,0,1,0,0,0,1
+1,0,0,0,0,-1,0,1,0
+0,0,1,0,1,0,-1,0,0
+0,-1,0,1,0,0,0,0,1
+1,0,0,0,-1,0,0,0,-1
+-1,0,0,0,1,0,0,0,-1
+-1,0,0,0,-1,0,0,0,1
+1,0,0,0,0,1,0,-1,0
+0,0,-1,0,1,0,1,0,0
+0,1,0,-1,0,0,0,0,1
+0,0,-1,1,0,0,0,-1,0
+0,0,1,-1,0,0,0,-1,0
+0,0,-1,-1,0,0,0,1,0
+0,0,1,1,0,0,0,1,0
+0,1,0,0,0,-1,-1,0,0
+0,-1,0,0,0,-1,1,0,0
+0,-1,0,0,0,1,-1,0,0
+0,1,0,0,0,1,1,0,0
+0,1,0,1,0,0,0,0,-1
+0,-1,0,-1,0,0,0,0,-1
+0,0,1,0,-1,0,1,0,0
+0,0,-1,0,-1,0,-1,0,0
+-1,0,0,0,0,1,0,1,0
+-1,0,0,0,0,-1,0,-1,0

testing.jl

@@ -0,0 +1,23 @@
+include("Hamiltonian.jl")
+
+println("Running with ",Threads.nthreads()," thread(s)")
+
+T=Float32
+
+function V_test(r2)
+    return -4*exp(-r2/4)
+end
+
+n = 3
+N = 8
+println("\n$n-body system with N=$N")
+
+for L::T in [16]
+    println("L=$L")
+    println("Constructing Hamiltonian")
+    s=system{T}(3,n,N,L,0.5,A1)
+    @time H=Hamiltonian{T}(s,V_test,0,0,cpu_tensor)
+    println("Solving eigenvalues")
+    @time evals,_,_ = eig(H,5)
+    println(evals)
+end