Simple reduction (still wrong)

common.jl

@@ -15,17 +15,18 @@ struct system{T}
     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 = calculate_all_data(N)
+            unique_i, unique_point, multiplicity, labels = calculate_all_data(N)
         elseif sym == A1
-            unique_i, unique_point, multiplicity = calculate_A1_data(N)
+            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)
+        return new{T}(d, n, N, convert(T, L), convert(T, μ), sym, unique_i, unique_point, multiplicity, labels)
     end
 end
 
@@ -90,7 +91,6 @@ function calculate_Vs(s::system{T}, V_twobody::Function, ϕ::T, n_image::Int)::A
                 end
             end
         end
-        Vs[i] *= s.multiplicity[i[1]]
     end
     return Vs
 end

irrep.jl

@@ -3,12 +3,13 @@ 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
+    return unique_i, unique_point, multiplicity, labels
 end
 
 function calculate_A1_data(N::Int)
@@ -38,7 +39,7 @@ function calculate_A1_data(N::Int)
     multiplicity = [count(labels.==i) for i in unique_i]
     unique_point = transpose(lattice[:, unique_i])
 
-    return unique_i, unique_point, multiplicity
+    return unique_i, unique_point, multiplicity, labels
 end
 
 function sym_reduce(s, K_partial)
@@ -47,13 +48,12 @@ function sym_reduce(s, K_partial)
     K_partial_y = kron(kron(I, K_partial), I)
     K_partial_z = kron(kron(I, I), K_partial)
 
-    for (i, j) in enumerate(s.unique_i)
-        K_partial_x[j, :] *= sqrt(s.multiplicity[i])
-        K_partial_x[:, j] *= sqrt(s.multiplicity[i])
-        K_partial_y[j, :] *= sqrt(s.multiplicity[i])
-        K_partial_y[:, j] *= sqrt(s.multiplicity[i])
-        K_partial_z[j, :] *= sqrt(s.multiplicity[i])
-        K_partial_z[:, j] *= sqrt(s.multiplicity[i])
+    for (i, label) in enumerate(s.labels)
+        if i != label
+            K_partial_x[:, i] += K_partial_x[:, label]
+            K_partial_y[:, i] += K_partial_y[:, label]
+            K_partial_z[:, i] += K_partial_z[:, label]
+        end
     end
 
     K_partial_x = K_partial_x[s.unique_i, s.unique_i]

testing.jl

@@ -12,7 +12,7 @@ n = 2
 N = 16
 println("\n$n-body system with N=$N")
 
-for L::T in 5.0:9.0
+for L::T in [16]
     println("L=$L")
     println("Constructing Hamiltonian")
     s=system{T}(3,n,N,L,0.5,A1)