diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp
index 1776fa63724..35957c6d436 100644
--- a/src/core/particle_data.cpp
+++ b/src/core/particle_data.cpp
@@ -821,42 +821,85 @@ int remove_particle(int p_id) {
   return ES_OK;
 }
 
-void local_remove_particle(int part) {
-  int ind, c;
-  Particle *p = local_particles[part];
-  ParticleList *pl = nullptr, *tmp;
+/**
+ * @brief Extract an indexed particle from a list.
+ *
+ * Removes a particle from a particle list and
+ * from the particle index.
+ *
+ * @param i Index of particle to remove,
+ *          needs to be valid.
+ * @param sl List to remove the particle from,
+ *           needs to be non-empty.
+ * @return The extracted particle.
+ */
+Particle extract_indexed_particle(ParticleList *sl, int i) {
+  assert(sl->n > 0);
+  assert(i < sl->n);
+  Particle *src = &sl->part[i];
+  Particle *end = &sl->part[sl->n - 1];
 
-  /* the tricky - say ugly - part: determine
-     the cell the particle is located in by checking
-     whether the particle address is inside the array */
-  for (c = 0; c < local_cells.n; c++) {
-    tmp = local_cells.cell[c];
-    ind = p - tmp->part;
-    if (ind >= 0 && ind < tmp->n) {
-      pl = tmp;
-      break;
-    }
+  Particle p = std::move(*src);
+
+  assert(p.p.identity <= max_seen_particle);
+  local_particles[p.p.identity] = nullptr;
+
+  if (src != end) {
+    new (src) Particle(std::move(*end));
   }
-  if (!pl) {
-    fprintf(stderr,
-            "%d: INTERNAL ERROR: could not find cell of particle %d, exiting\n",
-            this_node, part);
-    errexit();
+
+  if (realloc_particlelist(sl, --sl->n)) {
+    update_local_particles(sl);
+  } else if (src != end) {
+    local_particles[src->p.identity] = src;
   }
+  return p;
+}
 
-  free_particle(p);
+namespace {
+std::pair<Cell *, size_t> find_particle(Particle *p, Cell *c) {
+  for (int i = 0; i < c->n; ++i) {
+    if ((c->part + i) == p) {
+      return {c, i};
+    }
+  }
+  return {nullptr, 0};
+}
+
+std::pair<Cell *, size_t> find_particle(Particle *p, CellPList cells) {
+  for (auto &c : cells) {
+    auto res = find_particle(p, c);
+    if (res.first) {
+      return res;
+    }
+  }
 
-  /* remove local_particles entry */
-  local_particles[p->p.identity] = nullptr;
+  return {nullptr, 0};
+}
+} // namespace
 
-  if (&pl->part[pl->n - 1] != p) {
-    /* move last particle to free position */
-    *p = pl->part[pl->n - 1];
+void local_remove_particle(int part) {
+  Particle *p = local_particles[part];
+  assert(p);
+  assert(not p->l.ghost);
+
+  /* If the particles are sorted we can use the
+   * cell system to find the cell containing the
+   * particle. Otherwise we do a brute force search
+   * of the cells. */
+  Cell *cell = nullptr;
+  size_t n = 0;
+  if (Cells::RESORT_NONE == get_resort_particles()) {
+    std::tie(cell, n) = find_particle(p, find_current_cell(*p));
+  }
 
-    /* update the local_particles array for the moved particle */
-    local_particles[p->p.identity] = p;
+  if (not cell) {
+    std::tie(cell, n) = find_particle(p, local_cells);
   }
-  pl->n--;
+
+  assert(cell && cell->part && (n < cell->n) && ((cell->part + n) == p));
+
+  Particle p_destroy = extract_indexed_particle(cell, n);
 }
 
 void local_place_particle(int part, const double p[3], int _new) {
diff --git a/testsuite/python/particle.py b/testsuite/python/particle.py
index c899798db6c..f544109ac48 100644
--- a/testsuite/python/particle.py
+++ b/testsuite/python/particle.py
@@ -248,7 +248,7 @@ def test_parallel_property_setters(self):
             "If this test hangs, there is an mpi deadlock in a particle property setter.")
         for p in espressomd.particle_data.particle_attributes:
             # Uncomment to identify guilty property
-            #print( p)
+            # print( p)
 
             if not hasattr(s.part[0], p):
                 raise Exception(
@@ -260,7 +260,19 @@ def test_parallel_property_setters(self):
             # Cause a differtn mpi callback to uncover deadlock immediately
             x = getattr(s.part[:], p)
 
+    def test_zz_remove_all(self):
+        for id in self.system.part[:].id:
+            self.system.part[id].remove()
+        self.system.part.add(pos=np.random.random(
+            (100, 3)) * self.system.box_l, id=np.arange(100, dtype=int))
+        ids = self.system.part[:].id
+        np.random.shuffle(ids)
+        for id in ids:
+            self.system.part[id].remove()
+        with self.assertRaises(Exception):
+            self.system.part[17].remove()
+
 
 if __name__ == "__main__":
-    #print("Features: ", espressomd.features())
+    # print("Features: ", espressomd.features())
     ut.main()