diff --git a/Hamiltonian.jl b/Hamiltonian.jl
index db5c174..b0954d9 100644
--- a/Hamiltonian.jl
+++ b/Hamiltonian.jl
@@ -11,6 +11,7 @@ struct Hamiltonian{T}
     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
@@ -20,7 +21,7 @@ struct Hamiltonian{T}
         k = -s.N÷2:s.N÷2-1
         hermitian = ϕ == 0.0
         K_partial = (exp(-im * convert(T, ϕ)) * im / sqrt(2 * s.μ)) .* ∂_1DOF.(Ref(s), k, k')
-        K_partial_1, K_partial_2 = sym_reduce(s, K_partial)
+        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'])
@@ -30,7 +31,7 @@ struct Hamiltonian{T}
             K_diag = nothing
             K_mixed = nothing    
         end
-        return new{T}(s, K_partial, K_diag, K_mixed, K_partial_1, K_partial_2, 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
 
@@ -62,8 +63,14 @@ function LinearAlgebra.mul!(out::Array{Complex{T}}, H::Hamiltonian{T}, v::Array{
                 end
                 nconList_v[i2] = 2
 
-                tensor1 = coord1 == 1 ? H.K_partial_1[dim] : H.K_partial
-                tensor2 = coord2 == 1 ? H.K_partial_2[dim] : H.K_partial
+                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)
diff --git a/irrep.jl b/irrep.jl
index 293b1fd..fe922fa 100644
--- a/irrep.jl
+++ b/irrep.jl
@@ -49,6 +49,10 @@ function sym_reduce(s, K_partial)
     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)
@@ -61,13 +65,13 @@ function sym_reduce(s, K_partial)
         end
     end
 
-    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 = K_partial_x2[:, s.unique_i]
     K_partial_y2 = K_partial_y2[:, s.unique_i]
     K_partial_z2 = 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_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
\ No newline at end of file
diff --git a/testing.jl b/testing.jl
index 416f5bb..b951ece 100644
--- a/testing.jl
+++ b/testing.jl
@@ -8,8 +8,8 @@ function V_test(r2)
     return -4*exp(-r2/4)
 end
 
-n = 2
-N = 16
+n = 3
+N = 8
 println("\n$n-body system with N=$N")
 
 for L::T in [16]