diff --git a/common.jl b/common.jl
index d187f99..6a63689 100644
--- a/common.jl
+++ b/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