diff --git a/examples/lj.py b/examples/lj.py
index 97b66a45df42fede9c61a43d4fe73ba9ec8a3d22..214966e0634db3a81c9a1faf7988a23fc5751a73 100644
--- a/examples/lj.py
+++ b/examples/lj.py
@@ -5,7 +5,7 @@ import sys
 def lj(i, j):
     sr2 = 1.0 / squared_distance(i, j)
     sr6 = sr2 * sr2 * sr2 * sigma6[i, j]
-    apply(force, delta(i, j) * 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon[i, j])
+    apply(force, delta(i, j) * (48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon[i, j]))
 
 def euler(i):
     linear_velocity[i] += dt * force[i] / mass[i]
@@ -36,6 +36,7 @@ psim.add_feature_property('type', 'sigma6', pairs.real(), [epsilon for i in rang
 psim.set_domain([0.0, 0.0, 0.0, 53.747078, 53.747078, 53.747078])
 psim.read_particle_data("data/minimd_setup_32x32x32.input", ['type', 'mass', 'position', 'linear_velocity', 'flags'], pairs.point_mass())
 psim.reneighbor_every(20)
+#psim.compute_half()
 psim.build_neighbor_lists(cutoff_radius + skin)
 psim.vtk_output(f"output/lj_{target}")
 psim.compute(lj, cutoff_radius)
diff --git a/src/pairs/ir/apply.py b/src/pairs/ir/apply.py
index 43c21e446247ac3031818b7e0660625ffe9f18ec..6ae09dddf81dc77dc7d2c127a06de9ecbb79d3b0 100644
--- a/src/pairs/ir/apply.py
+++ b/src/pairs/ir/apply.py
@@ -44,11 +44,11 @@ class Apply(Lowerable):
 
     @pairs_inline
     def lower(self):
+        Assign(self.sim, self._reduction_variable, self._reduction_variable + self._expr)
+
         if self.sim._compute_half:
             for _ in Filter(self.sim,
                             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)
-
-        Assign(self.sim, self._reduction_variable, self._reduction_variable + self._expr)
diff --git a/src/pairs/sim/interaction.py b/src/pairs/sim/interaction.py
index fb5beccf8aea6fb66245f5383829157998a1c8ac..ff2f5c2dba828115b74d370528221d21f6d0fb8e 100644
--- a/src/pairs/sim/interaction.py
+++ b/src/pairs/sim/interaction.py
@@ -72,10 +72,8 @@ class NeighborFor:
                             particle_id = cell_particles[neigh_cell][cell_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))
+                                shape_cond = particle_shape[particle_id] > self.sim.get_shape_id(shape)
+                                condition = ScalarOp.or_op(shape_cond, self.particle < particle_id)
 
                             else:
                                 condition = ScalarOp.neq(particle_id, self.particle)
@@ -89,10 +87,8 @@ class NeighborFor:
                     inf_particle = cell_particles[0][inf_id]
 
                     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))
+                        shape_cond = particle_shape[inf_particle] > self.sim.get_shape_id(shape)
+                        condition = ScalarOp.or_op(shape_cond, self.particle < inf_particle)
 
                     else:
                         condition = ScalarOp.neq(inf_particle, self.particle)
@@ -266,7 +262,11 @@ class ParticleInteraction(Lowerable):
                             prop_reductions[prop] = prop_reductions[prop] + reduction
 
                     for prop, reduction in prop_reductions.items():
-                        Assign(self.sim, prop[i], reduction)
+                        if self.sim._compute_half:
+                            Assign(self.sim, prop[i], prop[i] + reduction)
+
+                        else:
+                            Assign(self.sim, prop[i], reduction)
 
         else:
             raise Exception("Interactions among more than two particles are currently not supported.")