diff --git a/examples/dem.py b/examples/dem.py
index ecd4889e00b5e0abb6d50040d3777187cbf661f4..60f86c0c3f899927b0375d71f9afe31a6d8f4a43 100644
--- a/examples/dem.py
+++ b/examples/dem.py
@@ -137,7 +137,6 @@ psim = pairs.simulation(
     double_prec=True,
     use_contact_history=True)
 
-#psim = pairs.simulation("dem", debug=True, timesteps=timeSteps)
 psim.add_position('position')
 psim.add_property('mass', pairs.real(), 1.0)
 psim.add_property('linear_velocity', pairs.vector())
@@ -177,6 +176,7 @@ psim.setup(update_mass_and_inertia, {'densityParticle_SI': densityParticle_SI,
                                      'pi': math.pi,
                                      'infinity': math.inf })
 
+#psim.compute_half()
 psim.build_neighbor_lists(linkedCellWidth + skin)
 psim.vtk_output(f"output/dem_{target}", frequency=visSpacing)
 
diff --git a/src/pairs/ir/apply.py b/src/pairs/ir/apply.py
index e2658e7298c22b845db4170e2ba343dc472f4ba3..43c21e446247ac3031818b7e0660625ffe9f18ec 100644
--- a/src/pairs/ir/apply.py
+++ b/src/pairs/ir/apply.py
@@ -7,6 +7,7 @@ from pairs.ir.properties import Property
 from pairs.ir.scalars import ScalarOp
 from pairs.ir.vectors import Vector
 from pairs.ir.types import Types
+from pairs.sim.flags import Flags
 from pairs.sim.lowerable import FinalLowerable, Lowerable
 
 
@@ -45,8 +46,8 @@ class Apply(Lowerable):
     def lower(self):
         if self.sim._compute_half:
             for _ in Filter(self.sim,
-                            ScalarOp.and_op(j < self.sim.nlocal,
-                                            ScalarOp.cmp(self.sim.particle_flags[j] & Flags.Fixed, 0))):
+                            ScalarOp.and_op(self._j < self.sim.nlocal,
+                                            ScalarOp.cmp(self.sim.particle_flags[self._j] & Flags.Fixed, 0))):
 
                 Assign(self.sim, self._prop[self._j], self._prop[self._j] - self._expr)
 
diff --git a/src/pairs/sim/interaction.py b/src/pairs/sim/interaction.py
index 3fcf52fb6b5e044edd242106353de56f8f459bf5..fb5beccf8aea6fb66245f5383829157998a1c8ac 100644
--- a/src/pairs/sim/interaction.py
+++ b/src/pairs/sim/interaction.py
@@ -59,6 +59,7 @@ class NeighborFor:
             ncells = self.cell_lists.ncells
             particle_cell = self.cell_lists.particle_cell
             cell_particles = self.cell_lists.cell_particles
+            particle_shape = self.sim.particle_shape
             nshapes = self.cell_lists.nshapes
 
             for shape in range(self.sim.max_shapes()):
@@ -67,17 +68,36 @@ class NeighborFor:
 
                     for _ in Filter(self.sim, ScalarOp.and_op(neigh_cell > 0, neigh_cell < ncells)):
                         start = sum([nshapes[neigh_cell][s] for s in range(shape)], 0)
-
                         for cell_particle in For(self.sim, start, start + nshapes[neigh_cell][shape]):
                             particle_id = cell_particles[neigh_cell][cell_particle]
-                            for _ in Filter(self.sim, ScalarOp.neq(particle_id, self.particle)):
+
+                            if self.sim._compute_half:
+                                shape_cond = particle_shape[particle_id] > shape
+                                condition = ScalarOp.or_op(shape_cond,
+                                                           ScalarOp.and_op(ScalarOp.not_op(shape_cond),
+                                                                           self.particle < particle_id))
+
+                            else:
+                                condition = ScalarOp.neq(particle_id, self.particle)
+
+                            for _ in Filter(self.sim, condition):
                                 yield Neighbor(self.sim, cell_particle, neigh_cell, particle_id, shape)
 
                 # Infinite particles
                 start = sum([nshapes[0][s] for s in range(shape)], 0)
                 for inf_id in For(self.sim, start, start + nshapes[0][shape]):
                     inf_particle = cell_particles[0][inf_id]
-                    for _ in Filter(self.sim, ScalarOp.neq(inf_particle, self.particle)):
+
+                    if self.sim._compute_half:
+                        shape_cond = particle_shape[inf_particle] > shape
+                        condition = ScalarOp.or_op(shape_cond,
+                                                   ScalarOp.and_op(ScalarOp.not_op(shape_cond),
+                                                                   self.particle < inf_particle))
+
+                    else:
+                        condition = ScalarOp.neq(inf_particle, self.particle)
+
+                    for _ in Filter(self.sim, condition):
                         yield Neighbor(self.sim, inf_id, 0, inf_particle, shape)
 
         else:
diff --git a/src/pairs/transformations/loops.py b/src/pairs/transformations/loops.py
index 829f6f7e11c6831303936e9062a12941ab289c2a..238effd920576d8f4933940d8066ae0c0c92fa73 100644
--- a/src/pairs/transformations/loops.py
+++ b/src/pairs/transformations/loops.py
@@ -42,7 +42,10 @@ class LICM(Mutator):
             Matrix,
             MatrixOp,
             PropertyAccess,
-            Quaternion,
+            # TODO: This is commented to avoid errors from the CUDA Explicit Euler kernel used
+            # for the DEM test cases. Allow lifting of Vector, Quaternion and Matrix until the
+            # kernel loops only (maybe other elements as well)
+            #Quaternion,
             QuaternionOp,
             ScalarOp,
             Select,