diff --git a/src/run.jl b/src/run.jl
index 4daa814..570b245 100644
--- a/src/run.jl
+++ b/src/run.jl
@@ -1,6 +1,7 @@
 using Random: Random
 using JSON3: JSON3
 using JLD2: JLD2
+using StaticArrays: SVector
 
 function run_sim(
     dir::String;
@@ -40,6 +41,13 @@ function run_sim(
     bundles_dir = "$dir/bundles"
     bundle = JLD2.load_object("$bundles_dir/bundle_$n_bundles.jld2")
 
+    particles = generate_particles(bundle, sim_consts.N)
+
+    particles_c = [SVector(0.0, 0.0) for i in 1:(sim_consts.N)]
+    for i in 1:(sim_consts.N)
+        particles_c[i] = SVector(particles[i].c[1], particles[i].c[2])
+    end
+
     args = (
         v=sim_consts.v,
         skin_r=skin_r,
@@ -56,10 +64,12 @@ function run_sim(
         N=sim_consts.N,
         l=sim_consts.l,
         particle_diameter=sim_consts.particle_diameter,
-        particles=generate_particles(bundle, sim_consts.N),
+        particles=particles,
+        particles_c=particles_c,
         verlet_list=[PreVector{Int64}(sim_consts.N - 1) for i in 1:(sim_consts.N - 1)],
         n_bundle_snapshots=n_bundle_snapshots,
         bundle=Bundle(sim_consts.N, n_bundle_snapshots),
+        box=Box(SVector(2 * sim_consts.l, 2 * sim_consts.l), skin_r),
     )
 
     T0::Float64 = status.T
diff --git a/src/simulation.jl b/src/simulation.jl
index 7512ecb..e301c3a 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -1,30 +1,43 @@
+using Dates: Dates, now
 using ProgressMeter: @showprogress
 using LoopVectorization: @turbo
 using Distributions: Normal
-using Dates: Dates, now
+using CellListMap: Box, CellList, map_pairwise!, UpdateCellList!
+
 import Base: wait
 
 rand_normal01() = rand(Normal(0, 1))
 
-function update_verlet_list!(args)
+function push_to_verlet_list(verlet_list, i, j)
+    if i < j
+        push!(verlet_list[i], Int64(j))
+    else
+        push!(verlet_list[j], Int64(i))
+    end
+
+    return nothing
+end
+
+function update_verlet_list!(args, cl)
     @simd for pv in args.verlet_list
         reset!(pv)
     end
 
-    for i in 1:(args.N - 1)
-        p1 = args.particles[i]
-        for j in (i + 1):(args.N)
-            p2 = args.particles[j]
-
-            overlapping = are_overlapping(p1, p2, args.skin_r², args.l).overlapping
-
-            if overlapping
-                push!(args.verlet_list[i], j)
-            end
-        end
+    for i in 1:(args.N)
+        args.particles_c[i] = SVector(args.particles[i].c[1], args.particles[i].c[2])
     end
 
-    return nothing
+    cl = UpdateCellList!(args.particles_c, args.box, cl; parallel=false)
+
+    map_pairwise!(
+        (x, y, i, j, d2, output) -> push_to_verlet_list(args.verlet_list, i, j),
+        nothing,
+        args.box,
+        cl;
+        parallel=false,
+    )
+
+    return cl
 end
 
 function euler!(args)
@@ -84,11 +97,13 @@ function simulate(
 )
     bundle_snapshot_counter = 0
 
+    task::Union{Task,Nothing} = nothing
+
+    cl = CellList(args.particles_c, args.box; parallel=false)
+
     start_time = now()
     println("Started simulation at $start_time.")
 
-    task::Union{Task,Nothing} = nothing
-
     @showprogress 0.6 for (integration_step, t) in enumerate(T0:δt:T)
         if (integration_step % n_steps_before_snapshot == 0) && save_data
             wait(task)
@@ -106,7 +121,7 @@ function simulate(
         end
 
         if integration_step % n_steps_before_verlet_list_update == 0
-            update_verlet_list!(args)
+            cl = update_verlet_list!(args, cl)
         end
 
         euler!(args)