diff --git a/examples/modular/sphere_box_global.cpp b/examples/modular/sphere_box_global.cpp
index 44390624e809ebfaf1abcd156222ea38d2f37816..0edb395355cd15b75772f4cce95ad84a24a814fe 100644
--- a/examples/modular/sphere_box_global.cpp
+++ b/examples/modular/sphere_box_global.cpp
@@ -4,7 +4,7 @@
 
 #include "sphere_box_global.hpp"
 
-// cmake -DINPUT_SCRIPT=../examples/modular/sphere_box_global.py -DWALBERLA_DIR=../../walberla -DBUILD_APP=ON -DUSER_SOURCE_FILES=../examples/modular/sphere_box_global.cpp -DCOMPILE_CUDA=ON ..
+// cmake -DINPUT_SCRIPT=../examples/modular/sphere_box_global.py DUSER_SOURCE_FILES=../examples/modular/sphere_box_global.cpp -DUSE_WALBERLA=ON -DCOMPILE_CUDA=ON ..
 
 void set_feature_properties(std::shared_ptr<PairsAccessor> &ac){
     ac->setTypeStiffness(0,0, 1e6);
diff --git a/src/pairs/sim/global_interaction.py b/src/pairs/sim/global_interaction.py
index 3bb6594f48f4ff38180f9fba61563b884177d9d1..cfbc030c8c35f03e315cb24c7f5c436b506d92e6 100644
--- a/src/pairs/sim/global_interaction.py
+++ b/src/pairs/sim/global_interaction.py
@@ -23,11 +23,6 @@ class GlobalLocalInteraction(ParticleInteraction):
     @pairs_device_block
     def lower(self):
         self.sim.module_name(f"{self.module_name}_global_local_interactions")
-        if self.sim._target.is_gpu():
-            first_cell_bytes = self.sim.add_temp_var(0)
-            Assign(self.sim, first_cell_bytes, self.cell_lists.cell_capacity * Sizeof(self.sim, Types.Int32))
-            CopyArray(self.sim, self.cell_lists.cell_sizes, Contexts.Host, Actions.ReadOnly, first_cell_bytes)
-        
         for ishape in range(self.maxs): # shape of globals
             if self.include_shape(ishape):
                 for jshape in range(self.maxs): # shape of locals
@@ -47,7 +42,7 @@ class GlobalLocalInteraction(ParticleInteraction):
                                         ScalarOp.cmp(self.sim.particle_shape[j], self.sim.get_shape_id(jshape)),
                                         ScalarOp.not_op(self.sim.particle_flags[j] & (Flags.Infinite | Flags.Global)))):
                                         for _ in Filter(self.sim, ScalarOp.neq(i, j)):
-                                            self.compute_interaction(i, j, ishape, jshape)
+                                            self.compute_interaction(i, j, ishape, jshape, atomic=True)
 
 
 class GlobalGlobalInteraction(ParticleInteraction):
diff --git a/src/pairs/sim/interaction.py b/src/pairs/sim/interaction.py
index 49aff87d50c1f6b7d508e54b3df3347d1004869d..1c350cc04174f9c60580865a84e5bf3a48038dc6 100644
--- a/src/pairs/sim/interaction.py
+++ b/src/pairs/sim/interaction.py
@@ -9,6 +9,7 @@ from pairs.ir.select import Select
 from pairs.ir.types import Types
 from pairs.ir.vectors import Vector
 from pairs.ir.print import Print
+from pairs.ir.atomic import AtomicInc
 from pairs.sim.flags import Flags
 from pairs.sim.lowerable import Lowerable
 from pairs.sim.shapes import Shapes
@@ -374,7 +375,7 @@ class ParticleInteraction(Lowerable):
         self.sim.leave()
         self.active_block = None
 
-    def apply_reductions(self, i, ishape, jshape):
+    def apply_reductions(self, i, ishape, jshape, atomic=False):
         prop_reductions = {}
         for app in self.apply_list[ishape*self.maxs + jshape]:
             prop = app.prop()
@@ -385,9 +386,18 @@ class ParticleInteraction(Lowerable):
                 prop_reductions[prop] = prop_reductions[prop] + reduction
 
         for prop, reduction in prop_reductions.items():
-            Assign(self.sim, prop[i], prop[i] + reduction)
+            if atomic:
+                if Types.is_scalar(prop.type()):
+                    AtomicInc(self.sim, prop[i], reduction)
+                    
+                else:
+                    for d in range(Types.number_of_elements(self.sim, prop.type())):
+                        AtomicInc(self.sim, prop[i][d], reduction[d])
+
+            else:
+                Assign(self.sim, prop[i], prop[i] + reduction)
 
-    def compute_interaction(self, i, j, ishape, jshape):
+    def compute_interaction(self, i, j, ishape, jshape, atomic=False):
         interaction_data = self.interactions_data[ishape*self.maxs + jshape]
         interaction_data.i().assign(i)
         interaction_data.j().assign(j)
@@ -416,7 +426,7 @@ class ParticleInteraction(Lowerable):
 
         # The i-j block is executed only if the cutoff_condition of the i-j interaction is met
         self.sim.add_statement(Filter(self.sim, interaction_data.cutoff_condition, self.blocks[ishape*self.maxs + jshape]))
-        self.apply_reductions(i, ishape, jshape)
+        self.apply_reductions(i, ishape, jshape, atomic)
 
     @pairs_block
     def lower(self):