diff --git a/graphics/radial_distribution.jl b/graphics/radial_distribution.jl
index c09afe8..1e1c63d 100644
--- a/graphics/radial_distribution.jl
+++ b/graphics/radial_distribution.jl
@@ -8,12 +8,12 @@ function gen_rdf_graphics()
     Random.seed!(1)
 
     box_length = 100
-    box_height = 100
+    box_length = 100
 
     graphics_export_dir = "exports/graphics"
     mkpath(graphics_export_dir)
 
-    Drawing(box_length, box_height, "$graphics_export_dir/rdf_shells.pdf")
+    Drawing(box_length, box_length, "$graphics_export_dir/rdf_shells.pdf")
     origin()
 
     particle_radius = 6
@@ -23,7 +23,6 @@ function gen_rdf_graphics()
     circle(Point(0, 0), particle_radius, :fill)
 
     N = 50
-    twice_max_particle_coordinate = box_length
 
     selected_shell_ind = 3
     selected_lower_radius = (selected_shell_ind - 1) * Δr
@@ -34,8 +33,8 @@ function gen_rdf_graphics()
 
     for p1_ind in 1:N
         while true
-            x = (rand() - 0.5) * twice_max_particle_coordinate
-            y = (rand() - 0.5) * twice_max_particle_coordinate
+            x = (rand() - 0.5) * box_length
+            y = (rand() - 0.5) * box_length
 
             p1_c = SVector(x, y)
 
@@ -43,7 +42,7 @@ function gen_rdf_graphics()
 
             no_collision = true
 
-            if distance > 2 * particle_radius
+            if distance >= 2 * particle_radius
                 for p2_ind in 1:(p1_ind - 1)
                     if ReCo.norm2d(p1_c - particle_cs[p2_ind]) < 2 * particle_radius
                         no_collision = false
diff --git a/graphics/verlet_and_cell_lists.jl b/graphics/verlet_and_cell_lists.jl
new file mode 100644
index 0000000..33b11b8
--- /dev/null
+++ b/graphics/verlet_and_cell_lists.jl
@@ -0,0 +1,94 @@
+using Luxor
+using Random: Random
+using StaticArrays: SVector
+
+using ReCo: ReCo
+
+function gen_verlet_and_cell_lists_graphics()
+    Random.seed!(2)
+
+    box_length = 100
+
+    graphics_export_dir = "exports/graphics"
+    mkpath(graphics_export_dir)
+
+    Drawing(box_length, box_length, "$graphics_export_dir/verlet_and_cell_lists.pdf")
+    origin()
+
+    R_particle = 4.2
+    σ = 2 * R_particle
+    R_interaction = 2^(1 / 6) * σ
+    R_skin = 2 * R_interaction
+
+    setcolor("blue")
+    reference_particle_x = 0.25 * R_skin
+    reference_particle_y = -0.1 * R_skin
+    circle(Point(reference_particle_x, reference_particle_y), R_particle, :fill)
+    reference_particle_c = SVector(reference_particle_x, reference_particle_y)
+
+    N = 92
+
+    particle_cs = Vector{SVector{2,Float64}}(undef, N)
+
+    x = y = distance = 0.0
+
+    for p1_ind in 1:N
+        while true
+            x = (rand() - 0.5) * box_length
+            y = (rand() - 0.5) * box_length
+
+            p1_c = SVector(x, y)
+
+            distance = ReCo.norm2d(p1_c - reference_particle_c)
+
+            no_collision = true
+
+            if distance >= 2 * R_particle
+                for p2_ind in 1:(p1_ind - 1)
+                    if ReCo.norm2d(p1_c - particle_cs[p2_ind]) < 2 * R_particle
+                        no_collision = false
+
+                        break
+                    end
+                end
+
+                if no_collision
+                    particle_cs[p1_ind] = p1_c
+                    break
+                end
+            end
+        end
+
+        if 2 * R_particle <= distance < R_interaction
+            setcolor("red")
+        elseif R_interaction <= distance < R_skin
+            setcolor("green")
+        else
+            setcolor("orange")
+        end
+
+        circle(x, y, R_particle, :fill)
+    end
+
+    setcolor("black")
+    setline(0.28)
+    for R in (R_interaction, R_skin)
+        circle(Point(reference_particle_x, reference_particle_y), R, :stroke)
+    end
+
+    half_n_lines = floor(Int64, box_length / (2 * R_skin) + 0.5)
+    for i in 1:half_n_lines
+        coordinate = (i - 0.5) * R_skin
+        line(Point(coordinate, -box_length), Point(coordinate, box_length), :stroke)
+        line(Point(-coordinate, -box_length), Point(-coordinate, box_length), :stroke)
+        line(Point(-box_length, coordinate), Point(box_length, coordinate), :stroke)
+        line(Point(-box_length, -coordinate), Point(box_length, -coordinate), :stroke)
+    end
+
+    setcolor((0.2, 0.0, 1.0, 0.3))
+    rect(-1.5 * R_skin, -1.5 * R_skin, 3 * R_skin, 3 * R_skin, :fill)
+
+    finish()
+
+    return nothing
+end
\ No newline at end of file
diff --git a/src/setup.jl b/src/setup.jl
index c2b040b..66386c3 100644
--- a/src/setup.jl
+++ b/src/setup.jl
@@ -67,21 +67,20 @@ function gen_sim_consts(
     ϵ = DEFAULT_ϵ
     interaction_r = 2^(1 / 6) * σ
 
-    buffer = 2.5
-    max_approach_after_one_integration_step = buffer * (2 * v₀ * δt) / interaction_r
-    @assert skin_to_interaction_r_ratio >= 1 + max_approach_after_one_integration_step
+    skin_r = skin_to_interaction_r_ratio * interaction_r
+
+    buffer = 2.2
+    max_approach_after_one_integration_step = buffer * (2 * v₀ * δt)
+    @assert skin_r >= interaction_r + max_approach_after_one_integration_step
 
     if v₀ != 0.0
-        n_steps_before_verlet_list_update = round(
-            Int64,
-            (skin_to_interaction_r_ratio - 1) / max_approach_after_one_integration_step,
+        n_steps_before_verlet_list_update = floor(
+            Int64, (skin_r - interaction_r) / max_approach_after_one_integration_step
         )
     else
         n_steps_before_verlet_list_update = 100
     end
 
-    skin_r = skin_to_interaction_r_ratio * interaction_r
-
     grid_n = round(Int64, ceil(sqrt(n_particles)))
 
     n_particles = grid_n^2
diff --git a/src/simulation.jl b/src/simulation.jl
index 3380694..0ddbad4 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -57,11 +57,11 @@ function euler!(
             state_update_helper_hook!(env_helper, id1, id2, r⃗₁₂)
 
             if overlapping
-                factor = args.ϵσ⁶δtμₜ24 / (distance²^4) * (args.σ⁶2 / (distance²^3) - 1.0)
-                dc = factor * r⃗₁₂
+                factor = args.ϵσ⁶δtμₜ24 / (distance²^4) * (1.0 - args.σ⁶2 / (distance²^3))
+                μₜF⃗₁₂ = factor * r⃗₁₂ # Force acting on 1 from 2 multiplied with μₜ
 
-                p1.tmp_c -= dc
-                p2.tmp_c += dc
+                p1.tmp_c += μₜF⃗₁₂
+                p2.tmp_c -= μₜF⃗₁₂
             end
         end
     end