diff --git a/src/boxgrid.jl b/src/boxgrid.jl
index 8118365394281ae870f3e98c6a447b2b03d38b6c..ec9cfa9e0d1dd0110916c8be66901b5d041bb708 100644
--- a/src/boxgrid.jl
+++ b/src/boxgrid.jl
@@ -158,3 +158,119 @@ function relevant_boxes(so::AbstractObstacle, bg::BoxGrid, w::AbstractDomain{3};
     @warn "Automatic box detection for 3D obstacles is not supported. Override `relevant_boxes` for your obstacle to manually select boxes" maxlog=1
     return vec(CartesianIndices(bg.boxes))
 end
+
+
+
+## BoxGrid pairwise interaction computation
+
+@forcedef pairwise_interactions(interaction::Function, sim::Simulation, bg::BoxGrid) =
+    pairwise_interactions(interaction, sim, bg, Val{isthreaded(sim)}())
+
+
+"""
+    pairwise_interactions(interaction::Function, sim::Simulation, bg::BoxGrid, ::Val{false})
+Uses a [`BoxGrid`](@ref) to efficiently compute pairwise interactions.
+`BoxGrid` structure is mutated in the sense that particles are sorted into their boxes.
+"""
+@forcedef function pairwise_interactions(interaction::Function, sim::Simulation, bg::BoxGrid, ::Val{false})
+    # For execution order preservation compared to multithreaded simulations the
+    # boxes are iterated in the same striped manner.
+    boxparticles!(bg, particles(sim))
+
+    numstripes, numrows... = size(bg.boxes)
+    rowit = CartesianIndices(numrows)
+
+    for i = 1:2:numstripes
+        for j = rowit
+            I = CartesianIndex(i,j)
+            isempty(bg.boxes[I]) || @forcefun pairwise_interactions(interaction, sim, bg.boxes[I])
+            for J ∈ bg.neighborlist[I]
+                if !isempty(bg.boxes[J])
+                    @forcefun pairwise_interactions(interaction, sim, bg.boxes[I], bg.boxes[J])
+                end
+            end # for J ∈ neighbors
+        end
+    end
+    for i = 2:2:numstripes
+        for j = rowit
+            I = CartesianIndex(i,j)
+            isempty(bg.boxes[I]) || @forcefun pairwise_interactions(interaction, sim, bg.boxes[I])
+            for J ∈ bg.neighborlist[I]
+                if !isempty(bg.boxes[J])
+                    @forcefun pairwise_interactions(interaction, sim, bg.boxes[I], bg.boxes[J])
+                end
+            end # for J ∈ neighbors
+        end
+    end
+end
+
+
+"""
+    pairwise_interactions(interaction::Function, sim::Simulation, bg::BoxGrid, ::Val{true})
+Multi-threaded version of [`pairwise_interactions!_serial`](@ref).
+
+## Extended help
+
+Multithreading is implemented by dividing the system into two sets of interleaved stripes (slabs in 3D)
+of boxes, each stripe being one box wide.
+This division ensures that each stripe is independent of all other stripes in same set.
+
+Therefore, we can first integrate all the stripes in the first set in parallel, synchronise, then
+integrate the second set.
+"""
+@forcedef function pairwise_interactions(interaction::Function, sim::Simulation, bg::BoxGrid, ::Val{true})
+    boxparticles!(bg, particles(sim))
+
+    numstripes, numrows... = size(bg.boxes)
+    rowit = CartesianIndices(numrows)
+
+    Threads.@threads for i = 1:2:numstripes
+        for j = rowit
+            I = CartesianIndex(i,j)
+            isempty(bg.boxes[I]) || @forcefun pairwise_interactions(interaction, sim, bg.boxes[I])
+            for J ∈ bg.neighborlist[I]
+                if !isempty(bg.boxes[J])
+                    @forcefun pairwise_interactions(interaction, sim, bg.boxes[I], bg.boxes[J])
+                end
+            end # for J ∈ neighbors
+        end
+    end
+    Threads.@threads for i = 2:2:numstripes
+        for j = rowit
+            I = CartesianIndex(i,j)
+            isempty(bg.boxes[I]) || @forcefun pairwise_interactions(interaction, sim, bg.boxes[I])
+            for J ∈ bg.neighborlist[I]
+                if !isempty(bg.boxes[J])
+                    @forcefun pairwise_interactions(interaction, sim, bg.boxes[I], bg.boxes[J])
+                end
+            end # for J ∈ neighbors
+        end
+    end
+end
+
+
+@forcedef function pairwise_interactions(interaction::Function, sim::Simulation, ::NoBoxing)
+    @forcefun pairwise_interactions(interaction, sim, particles(sim))
+end
+
+@forcedef function pairwise_interactions(interaction::Function, sim::Simulation,
+                                         particles::AbstractVector{<: AbstractParticle})
+    N = length(particles)
+    for i = 1:N
+        c1 = particles[i]
+        for j = i+1:N
+            @forcefun interaction(c1, particles[j], sim)
+        end
+    end
+end
+
+@forcedef function pairwise_interactions(interaction::Function, sim::Simulation,
+                                         particles1::AbstractVector{P},
+                                         particles2::AbstractVector{P}) where {P <: AbstractParticle}
+    for i in eachindex(particles1)
+        c1 = particles1[i]
+        for j in eachindex(particles2)
+            @forcefun interaction(c1, particles2[j], sim)
+        end
+    end
+end
diff --git a/src/evolve.jl b/src/evolve.jl
index 71de2d25f7591b285b36a776da25bc541a44e0c0..63210135aba6667a89b44d4fc40842ef4cb65245 100644
--- a/src/evolve.jl
+++ b/src/evolve.jl
@@ -153,52 +153,4 @@ end
     computeforces!(sim)
 Computes the external and internal forces in the simulation.
 """
-@forcedef function computeforces!(sim::Simulation)
-    if isthreaded(sim)
-        @forcefun threaded_computeforces!(sim, sim.particles)
-    else
-        @forcefun serial_computeforces!(sim, sim.particles)
-    end
-end
-
-# fallback for non-parallelised PCs
-@forcedef threaded_computeforces!(sim, pc) = @forcefun serial_computeforces!(sim, pc)
-
-
-"""
-    computeforces!(particles, sim, ::NoBoxing)
-Computes all external and internal forces of `particles` naively, computing the
-interactions for all pairs of particles.
-"""
-@forcedef function computeforces!(particles::AbstractVector{T},
-                             sim::Simulation,
-                             ::NoBoxing) where {T <: AbstractParticle}
-    N = length(particles)
-    for i = 1:N
-        c1 = particles[i]
-        @forcefun internalforces!(c1, sim)
-        for j = i+1:N
-            @forcefun externalforces!(c1, particles[j], sim)
-        end
-    end
-end
-
-
-
-"""
-    computeforces!(particles1, particles2, sim)
-Computes the external forces arising due to interaction of particles from `particles1`
-with particles from `particles2`.
-
-The function assumes that `particles1` and `particles2` do not share elements.
-This condition is only checked within an [`@assert`](@ref), so there might be no
-signs of things going wrong depending on the optimization level!
-"""
-@forcedef function computeforces!(particles1::AbstractVector{T}, particles2::AbstractVector{T},
-                             sim::Simulation) where {T <: AbstractParticle}
-
-    for c1 ∈ particles1, c2 ∈ particles2
-        @assert c1 != c2
-        @forcefun externalforces!(c1, c2, sim)
-    end
-end
+@forcedef computeforces!(sim::Simulation) = @forcefun computeforces!(sim, sim.particles)
diff --git a/src/particlecontainers/particlehandling_common.jl b/src/particlecontainers/particlehandling_common.jl
index 3ad94602e0729c67f68462a528c5b399b30af405..87d33064a3bf0686ebadb1095988f7e666c1ab47 100644
--- a/src/particlecontainers/particlehandling_common.jl
+++ b/src/particlecontainers/particlehandling_common.jl
@@ -127,125 +127,3 @@ function findbyid(particles::AbstractVector{<:AbstractParticle}, id::UInt64)
 end
 
 findbyid(particles::AbstractVector{<:AbstractParticle}, ids::AbstractSet{UInt64}) = sort!([findbyid(particles, id) for id ∈ ids])
-
-
-
-function splitbytype(particles::AbstractVector{T}) where {T<:AbstractParticle}
-    if !(T isa Union)
-        return Dict(T => particles)
-    end
-
-    ret = Dict{UInt64, Vector{<:AbstractParticle}}()
-    for p ∈ particles
-        S = typeof(p)
-        Ssymb = hash(S)
-        haskey(ret, Ssymb) || (ret[Ssymb] = S[])
-        push!(ret[Ssymb], p)
-    end
-    return values(ret)
-end
-
-function maskbytype(particles::AbstractVector{U}) where {U}
-    if !(U isa Union)
-        return U, trues(length(particles))
-    end
-
-    l = length(particles)
-    ut = Base.uniontypes(U)
-    utl = length(ut)
-    masks = [falses(l) for i ∈ 1:utl]
-
-    for i ∈ 1:utl
-        for j ∈ 1:l
-            masks[i][j] = isa(particles[j], ut[i])
-        end
-    end
-
-    return ut, masks
-end
-
-
-
-## Common implementations for BoxGrid based particles containers
-
-"""
-    boxing_based_serial_computeforces(sim, pc::AbstractParticleContainer)
-Uses a [`BoxGrid`](@ref) to efficiently compute pairwise interactions.
-"""
-@forcedef function boxing_based_serial_computeforces!(sim::Simulation, pc::AbstractParticleContainer)
-    # For execution order preservation compared to multithreaded simulations the
-    # boxes are iterated in the same striped manner.
-    bg = pc.boxgrid
-    boxparticles!(bg, particles(pc))
-
-    numstripes, numrows... = size(bg.boxes)
-    rowit = CartesianIndices(numrows)
-
-    for i = 1:2:numstripes
-        for j = rowit
-            I = CartesianIndex(i,j)
-            isempty(bg.boxes[I]) || @forcefun computeforces!(bg.boxes[I], sim, NoBoxing())
-            for J ∈ bg.neighborlist[I]
-                if !isempty(bg.boxes[J])
-                    @forcefun computeforces!(bg.boxes[I], bg.boxes[J], sim)
-                end
-            end # for J ∈ neighbors
-        end
-    end
-    for i = 2:2:numstripes
-        for j = rowit
-            I = CartesianIndex(i,j)
-            isempty(bg.boxes[I]) || @forcefun computeforces!(bg.boxes[I], sim, NoBoxing())
-            for J ∈ bg.neighborlist[I]
-                if !isempty(bg.boxes[J])
-                    @forcefun computeforces!(bg.boxes[I], bg.boxes[J], sim)
-                end
-            end # for J ∈ neighbors
-        end
-    end
-end
-
-
-"""
-    boxing_based_threaded_computeforces(sim, pc::AbstractParticleContainer)
-Multi-threaded version of [`boxing_based_serial_computeforces`](@ref).
-
-## Extended help
-
-Multithreading is implemented by dividing the system into two sets of interleaved stripes (slabs in 3D)
-of boxes, each stripe being one box wide.
-This division ensures that each stripe is independent of all other stripes in same set.
-
-Therefore, we can first integrate all the stripes in the first set in parallel, synchronise, then
-integrate the second set.
-"""
-@forcedef function boxing_based_threaded_computeforces!(sim::Simulation, pc::AbstractParticleContainer)
-    bg = pc.boxgrid
-    boxparticles!(bg, particles(pc))
-
-    numstripes, numrows... = size(bg.boxes)
-    rowit = CartesianIndices(numrows)
-
-    Threads.@threads for i = 1:2:numstripes
-        for j = rowit
-            I = CartesianIndex(i,j)
-            isempty(bg.boxes[I]) || @forcefun computeforces!(bg.boxes[I], sim, NoBoxing())
-            for J ∈ bg.neighborlist[I]
-                if !isempty(bg.boxes[J])
-                    @forcefun computeforces!(bg.boxes[I], bg.boxes[J], sim)
-                end
-            end # for J ∈ neighbors
-        end
-    end
-    Threads.@threads for i = 2:2:numstripes
-        for j = rowit
-            I = CartesianIndex(i,j)
-            isempty(bg.boxes[I]) || @forcefun computeforces!(bg.boxes[I], sim, NoBoxing())
-            for J ∈ bg.neighborlist[I]
-                if !isempty(bg.boxes[J])
-                    @forcefun computeforces!(bg.boxes[I], bg.boxes[J], sim)
-                end
-            end # for J ∈ neighbors
-        end
-    end
-end
diff --git a/src/particlecontainers/particleobstaclecontainer.jl b/src/particlecontainers/particleobstaclecontainer.jl
index fd47eef2876606de7139bc7cc84e67478a0de431..ce31d27b7041213c4fef6b941d4a6cb157162a84 100644
--- a/src/particlecontainers/particleobstaclecontainer.jl
+++ b/src/particlecontainers/particleobstaclecontainer.jl
@@ -58,47 +58,30 @@ addobstacles!(ptc::ParticleObstacleContainer, obs) = append!(ptc.sov, obs)
 setparticles!(pc::ParticleObstacleContainer, ps) = (empty!(pc.pv); append!(pc.pv, ps))
 alldynamicobjects(sim::Simulation{<:AbstractParticleContainer}) = particles(sim)
 
-
-@forcedef function serial_computeforces!(sim::Simulation, pc::ParticleObstacleContainer{N, PT, SO, BG}) where {N, PT, SO, BG<:NoBoxing}
-    @forcefun computeforces!(pc.pv, sim, pc.boxgrid)
-    for n in axes(pc.sov, 1)
-        so = pc.sov[n]
-        @forcefun obstacleforces!(so, pc.pv, sim)
-    end
-end
-
-
-@forcedef function serial_computeforces!(sim::Simulation, pc::ParticleObstacleContainer{N, PT, SO, BG}) where {N, PT, SO, BG<:BoxGrid}
-    @forcefun boxing_based_serial_computeforces!(sim, pc)
-
+@forcedef function computeforces!(sim::Simulation, pc::ParticleObstacleContainer)
     bg = pc.boxgrid
-    for n in axes(pc.sov, 1)
-        so = pc.sov[n]
-        for idx in pc.so_boxes[n]
-            @forcefun obstacleforces!(so, bg.boxes[idx], sim)
-        end
+    @forcefun pairwise_interactions(externalforces!, sim, bg)
+    for p in particles(sim)
+        internalforces!(p, sim)
     end
-end
-
-
-@forcedef function threaded_computeforces!(sim::Simulation, pc::ParticleObstacleContainer{N, PT, SO, BG}) where {N, PT, SO, BG<:BoxGrid}
-    @forcefun boxing_based_threaded_computeforces!(sim, pc)
-
-    bg = pc.boxgrid
-    for n in axes(pc.sov, 1)
-        so = pc.sov[n]
-        for idx in pc.so_boxes[n]
-            @forcefun obstacleforces!(so, bg.boxes[idx], sim)
+    if pc.boxgrid isa BoxGrid
+        for n in axes(pc.sov, 1)
+            so = pc.sov[n]
+            for idx in pc.so_boxes[n]
+                # particles are sorted into boxes, within pairwise_interactions
+                @forcefun obstacleforces!(so, bg.boxes[idx], sim)
+            end
+        end
+    else
+        for n in axes(pc.sov, 1)
+            so = pc.sov[n]
+            @forcefun obstacleforces!(so, particles(sim), sim)
         end
     end
 end
 
-
-
 obstacles(pc::ParticleObstacleContainer) = pc.sov
 
-
-
 function Base.show(io::IO, ::MIME"text/plain", pc::ParticleObstacleContainer{N, PT, SO, BG}) where {N, PT, SO, BG}
     print(io,
 """ParticleObstacleContainer{$N, $PT, $SO, $BG}
diff --git a/src/particlecontainers/polocontainer.jl b/src/particlecontainers/polocontainer.jl
index 44202a2c775f3283c568d5cd2f361aef4cd5aeef..aa6d7c4fcf85780c0edf15dc076bbd95e80fee30 100644
--- a/src/particlecontainers/polocontainer.jl
+++ b/src/particlecontainers/polocontainer.jl
@@ -60,33 +60,17 @@ setparticles!(pc::POLOContainer, ps) = (empty!(pc.pv); append!(pc.pv, ps))
 alldynamicobjects(sim::Simulation{<:POLOContainer}) = vcat(sim.particles.pv, sim.particles.lov)
 
 
-@forcedef function serial_computeforces!(sim::Simulation, pc::POLOContainer)
-    @forcefun boxing_based_serial_computeforces!(sim, pc)
-
-    bg = pc.boxgrid
-    for n in axes(pc.sov, 1)
-        so = pc.sov[n]
-        for idx in pc.so_boxes[n]
-            @forcefun obstacleforces!(so, bg.boxes[idx], sim)
-        end
-    end
-
-    for n in axes(pc.lov, 1)
-        lo = pc.lov[n]
-        for cell in particles(sim)
-            @forcefun externalforces!(lo, cell, sim)
-        end
+@forcedef function computeforces!(sim::Simulation, pc::POLOContainer)
+    @forcefun pairwise_interactions(externalforces!, sim, pc.boxgrid)
+    for p in particles(sim)
+        internalforces!(p, sim)
     end
-end
-
-
-@forcedef function threaded_computeforces!(sim::Simulation, pc::POLOContainer)
-    @forcefun boxing_based_threaded_computeforces!(sim, pc)
 
     bg = pc.boxgrid
     for n in axes(pc.sov, 1)
         so = pc.sov[n]
         for idx in pc.so_boxes[n]
+            # particles are sorted into boxes, within pairwise_interactions
             @forcefun obstacleforces!(so, bg.boxes[idx], sim)
         end
     end
@@ -99,7 +83,6 @@ end
     end
 end
 
-
 obstacles(pc::POLOContainer) = pc.sov
 
 
diff --git a/src/particlecontainers/simpleparticlecontainer.jl b/src/particlecontainers/simpleparticlecontainer.jl
index d582f6e200b9e46564aea215fba1d52a15761c92..cd66c9870b21ce5e329ff00d6b5948f4a27182c6 100644
--- a/src/particlecontainers/simpleparticlecontainer.jl
+++ b/src/particlecontainers/simpleparticlecontainer.jl
@@ -47,22 +47,13 @@ alldynamicobjects(sim::Simulation{<:SimpleParticleContainer}) = sim.particles.pa
 
 ## FORCE COMPUTATION ###########################################################
 
-@forcedef function serial_computeforces!(sim::Simulation, pc::SimpleParticleContainer{N, PT, BG}) where {N, PT, BG<:NoBoxing}
-    @forcefun computeforces!(pc.particles, sim, pc.boxgrid)
+@forcedef function computeforces!(sim::Simulation, pc::SimpleParticleContainer)
+    @forcefun pairwise_interactions(externalforces!, sim, pc.boxgrid)
+    for p in particles(sim)
+        internalforces!(p, sim)
+    end
 end
 
-
-@forcedef function serial_computeforces!(sim::Simulation, pc::SimpleParticleContainer{N, PT, BG}) where {N, PT, BG<:BoxGrid}
-    @forcefun boxing_based_serial_computeforces!(sim, pc)
-end
-
-@forcedef function threaded_computeforces!(sim::Simulation, pc::SimpleParticleContainer{N, PT, BG}) where {N, PT, BG<:BoxGrid}
-    @forcefun boxing_based_threaded_computeforces!(sim, pc)
-end
-
-
-
-
 function Base.show(io::IO, ::MIME"text/plain", pc::SimpleParticleContainer{N, PT, BG}) where {N, PT, BG}
     print(io,
 """SimpleParticleContainer{$N, $PT, $BG}