diff --git a/src/Particle.jl b/src/Particle.jl
index 206c1c4..79e1e72 100644
--- a/src/Particle.jl
+++ b/src/Particle.jl
@@ -10,11 +10,6 @@ mutable struct Particle
     end
 end
 
-get_c(p::Particle) = p.c
-get_c_x(p::Particle) = p.c[1]
-get_c_y(p::Particle) = p.c[2]
-get_φ(p::Particle) = p.φ
-
 function restrict_coordinate(value::Float64; l::Float64)
     if value < -l
         value += 2 * l
diff --git a/src/run.jl b/src/run.jl
index dc189c4..20cdd29 100644
--- a/src/run.jl
+++ b/src/run.jl
@@ -47,7 +47,7 @@ function run(;
         particles = generate_particles(grid_n, grid_box_width, l),
         skin_r = skin_r,
         skin_r² = skin_r^2,
-        verlet_list = [PreVector(Int64, N - 1) for i in 1:N],
+        verlet_list = [PreVector(Int64, N - 1) for i in 1:(N - 1)],
         n_frames = floor(Int64, integration_steps / n_steps_before_save) + 1,
     )
 
diff --git a/src/simulation.jl b/src/simulation.jl
index 4b0681e..cc74a18 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -4,7 +4,7 @@ function update_verlet_list!(args)
     @simd for pv in args.verlet_list
         reset!(pv)
     end
-
+    
     for i in 1:(args.N - 1)
         for j in (i + 1):args.N
             p1 = args.particles[i]
@@ -14,29 +14,35 @@ function update_verlet_list!(args)
 
             if overlapping
                 push!(args.verlet_list[i], j)
-                push!(args.verlet_list[j], i)
             end
         end
     end
 end
 
 function euler!(args)
-    Threads.@threads for p in args.particles
+    for i in 1:(args.N - 1)
+        p = args.particles[i]
         verlet_list = args.verlet_list[p.id]
 
-        for i in 1:verlet_list.last_ind
-            p2 = args.particles[verlet_list.v[i]]
+        for j in 1:verlet_list.last_ind
+            p2 = args.particles[verlet_list.v[j]]
 
             overlapping, r⃗₁₂, distance² = are_overlapping(p, p2, args.interaction_r², args.l)
 
             if overlapping
                 c = args.c₁ / (distance²^4) * (args.c₂ / (distance²^3) - 1)
-                @simd for j in 1:2
-                    p.tmp_c[j] -= c * r⃗₁₂[j]
+                
+                @simd for k in 1:2
+                    dck = c * r⃗₁₂[k]
+                    
+                    p.tmp_c[k] -= dck
+                    p2.tmp_c[k] += dck
                 end
             end
         end
+    end
 
+    @simd for p in args.particles
         e = SVector(cos(p.φ), sin(p.φ))
 
         @simd for i in 1:2
@@ -46,10 +52,8 @@ function euler!(args)
         p.φ += args.c₄ * rand_normal01()
 
         restrict_coordinates!(p; l=args.l)
-    end
 
-    @simd for particle in args.particles
-        update!(particle)
+        update!(p)
     end
 
     return nothing