diff --git a/src/core/AtomDecomposition.cpp b/src/core/AtomDecomposition.cpp
index 9398eb35e15..ef519856387 100644
--- a/src/core/AtomDecomposition.cpp
+++ b/src/core/AtomDecomposition.cpp
@@ -147,6 +147,8 @@ AtomDecomposition::AtomDecomposition(const boost::mpi::communicator &comm,
   mark_cells();
 }
 
-Utils::Vector3d AtomDecomposition::max_range() const {
+Utils::Vector3d AtomDecomposition::max_cutoff() const {
   return Utils::Vector3d::broadcast(std::numeric_limits<double>::infinity());
 }
+
+Utils::Vector3d AtomDecomposition::max_range() const { return max_cutoff(); }
\ No newline at end of file
diff --git a/src/core/AtomDecomposition.hpp b/src/core/AtomDecomposition.hpp
index 0a1863b43e6..1c36e7f5828 100644
--- a/src/core/AtomDecomposition.hpp
+++ b/src/core/AtomDecomposition.hpp
@@ -92,7 +92,9 @@ class AtomDecomposition : public ParticleDecomposition {
     return id_to_cell(p.identity());
   }
 
+  Utils::Vector3d max_cutoff() const override;
   Utils::Vector3d max_range() const override;
+
   /* Return true if minimum image convention is
    * needed for distance calculation. */
   boost::optional<BoxGeometry> minimum_image_distance() const override {
diff --git a/src/core/CellStructure.cpp b/src/core/CellStructure.cpp
index a96c1c43f66..4448318c698 100644
--- a/src/core/CellStructure.cpp
+++ b/src/core/CellStructure.cpp
@@ -179,6 +179,10 @@ ParticleRange CellStructure::ghost_particles() {
   return Cells::particles(decomposition().ghost_cells());
 }
 
+Utils::Vector3d CellStructure::max_cutoff() const {
+  return decomposition().max_cutoff();
+}
+
 Utils::Vector3d CellStructure::max_range() const {
   return decomposition().max_range();
 }
diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp
index 492dae3235e..b3688a956fe 100644
--- a/src/core/CellStructure.hpp
+++ b/src/core/CellStructure.hpp
@@ -127,13 +127,12 @@ struct CellStructure {
   /** One of @ref Cells::Resort, announces the level of resort needed.
    */
   unsigned m_resort_particles = Cells::RESORT_NONE;
+  bool m_rebuild_verlet_list = true;
+  std::vector<std::pair<Particle *, Particle *>> m_verlet_list;
 
 public:
   bool use_verlet_list = true;
 
-  bool m_rebuild_verlet_list = true;
-  std::vector<std::pair<Particle *, Particle *>> m_verlet_list;
-
   /**
    * @brief Update local particle index.
    *
@@ -240,6 +239,9 @@ struct CellStructure {
 public:
   int decomposition_type() const { return m_type; }
 
+  /** Maximal cutoff supported by current cell system. */
+  Utils::Vector3d max_cutoff() const;
+
   /** Maximal pair range supported by current cell system. */
   Utils::Vector3d max_range() const;
 
@@ -328,6 +330,7 @@ struct CellStructure {
     return assert(m_decomposition), *m_decomposition;
   }
 
+private:
   ParticleDecomposition &decomposition() {
     return assert(m_decomposition), *m_decomposition;
   }
@@ -445,12 +448,11 @@ struct CellStructure {
       std::unique_ptr<ParticleDecomposition> &&decomposition) {
     clear_particle_index();
 
-    auto local_parts = local_particles();
-    std::vector<Particle> particles(local_parts.begin(), local_parts.end());
-
-    m_decomposition = std::move(decomposition);
+    /* Swap in new cell system */
+    std::swap(m_decomposition, decomposition);
 
-    for (auto &p : particles) {
+    /* Add particles to new system */
+    for (auto &p : Cells::particles(decomposition->local_cells())) {
       add_particle(std::move(p));
     }
   }
@@ -510,31 +512,21 @@ struct CellStructure {
     }
   }
 
-public:
-  /** Non-bonded pair loop with potential use
-   * of verlet lists.
-   * @param pair_kernel Kernel to apply
-   */
-  template <class PairKernel> void non_bonded_loop(PairKernel pair_kernel) {
-    link_cell(pair_kernel);
-  }
-
-  /** Non-bonded pair loop with potential use
-   * of verlet lists.
+  /** Non-bonded pair loop with verlet lists.
+   *
    * @param pair_kernel Kernel to apply
    * @param verlet_criterion Filter for verlet lists.
    */
   template <class PairKernel, class VerletCriterion>
-  void non_bonded_loop(PairKernel &&pair_kernel,
-                       const VerletCriterion &verlet_criterion) {
+  void verlet_list_loop(PairKernel pair_kernel,
+                        const VerletCriterion &verlet_criterion) {
     /* In this case the verlet list update is attached to
      * the pair kernel, and the verlet list is rebuilt as
      * we go. */
-    if (use_verlet_list && m_rebuild_verlet_list) {
+    if (m_rebuild_verlet_list) {
       m_verlet_list.clear();
 
-      link_cell([&pair_kernel, &verlet_criterion,
-                 this](Particle &p1, Particle &p2, Distance const &d) {
+      link_cell([&](Particle &p1, Particle &p2, Distance const &d) {
         if (verlet_criterion(p1, p2, d)) {
           m_verlet_list.emplace_back(&p1, &p2);
           pair_kernel(p1, p2, d);
@@ -542,7 +534,7 @@ struct CellStructure {
       });
 
       m_rebuild_verlet_list = false;
-    } else if (use_verlet_list && not m_rebuild_verlet_list) {
+    } else {
       auto const maybe_box = decomposition().minimum_image_distance();
       /* In this case the pair kernel is just run over the verlet list. */
       if (maybe_box) {
@@ -558,6 +550,27 @@ struct CellStructure {
                       distance_function(*pair.first, *pair.second));
         }
       }
+    }
+  }
+
+public:
+  /** Non-bonded pair loop.
+   * @param pair_kernel Kernel to apply
+   */
+  template <class PairKernel> void non_bonded_loop(PairKernel pair_kernel) {
+    link_cell(pair_kernel);
+  }
+
+  /** Non-bonded pair loop with potential use
+   * of verlet lists.
+   * @param pair_kernel Kernel to apply
+   * @param verlet_criterion Filter for verlet lists.
+   */
+  template <class PairKernel, class VerletCriterion>
+  void non_bonded_loop(PairKernel pair_kernel,
+                       const VerletCriterion &verlet_criterion) {
+    if (use_verlet_list) {
+      verlet_list_loop(pair_kernel, verlet_criterion);
     } else {
       /* No verlet lists, just run the kernel with pairs from the cells. */
       link_cell(pair_kernel);
diff --git a/src/core/DomainDecomposition.cpp b/src/core/DomainDecomposition.cpp
index 1465b926783..89a14f784f1 100644
--- a/src/core/DomainDecomposition.cpp
+++ b/src/core/DomainDecomposition.cpp
@@ -262,13 +262,15 @@ void DomainDecomposition::fill_comm_cell_lists(ParticleList **part_lists,
         *part_lists++ = &(cells.at(i).particles());
       }
 }
-Utils::Vector3d DomainDecomposition::max_range() const {
+Utils::Vector3d DomainDecomposition::max_cutoff() const {
   auto dir_max_range = [this](int i) {
     return std::min(0.5 * m_box.length()[i], m_local_box.length()[i]);
   };
 
   return {dir_max_range(0), dir_max_range(1), dir_max_range(2)};
 }
+
+Utils::Vector3d DomainDecomposition::max_range() const { return cell_size; }
 int DomainDecomposition::calc_processor_min_num_cells() const {
   /* the minimal number of cells can be lower if there are at least two nodes
      serving a direction,
diff --git a/src/core/DomainDecomposition.hpp b/src/core/DomainDecomposition.hpp
index c8c2cf80623..0c6fe81252e 100644
--- a/src/core/DomainDecomposition.hpp
+++ b/src/core/DomainDecomposition.hpp
@@ -112,6 +112,7 @@ struct DomainDecomposition : public ParticleDecomposition {
   }
 
   void resort(bool global, std::vector<ParticleChange> &diff) override;
+  Utils::Vector3d max_cutoff() const override;
   Utils::Vector3d max_range() const override;
 
   boost::optional<BoxGeometry> minimum_image_distance() const override {
diff --git a/src/core/ParticleDecomposition.hpp b/src/core/ParticleDecomposition.hpp
index 0a6e853e4ce..8bc637c6d69 100644
--- a/src/core/ParticleDecomposition.hpp
+++ b/src/core/ParticleDecomposition.hpp
@@ -115,7 +115,13 @@ class ParticleDecomposition {
   /**
    * @brief Maximum supported cutoff.
    */
+  virtual Utils::Vector3d max_cutoff() const = 0;
+
+  /**
+   * @brief Range in which calculations are performed.
+   */
   virtual Utils::Vector3d max_range() const = 0;
+
   /**
    * @brief Return the box geometry needed for distance calculation
    *        if minimum image convention should be used needed for
diff --git a/src/core/cells.cpp b/src/core/cells.cpp
index 969492971b2..3ec3554b5fc 100644
--- a/src/core/cells.cpp
+++ b/src/core/cells.cpp
@@ -25,8 +25,10 @@
  *  Implementation of cells.hpp.
  */
 #include "cells.hpp"
+
 #include "Particle.hpp"
 #include "communication.hpp"
+#include "errorhandling.hpp"
 #include "event.hpp"
 #include "grid.hpp"
 #include "integrate.hpp"
@@ -40,6 +42,7 @@
 #include <utils/mpi/gather_buffer.hpp>
 
 #include <algorithm>
+#include <boost/range/algorithm/min_element.hpp>
 #include <functional>
 #include <stdexcept>
 #include <vector>
@@ -48,23 +51,24 @@
 CellStructure cell_structure;
 
 /**
- * @brief Get pairs closer than distance from the cells.
+ * @brief Get pairs of particles that are closer than a distance and fulfill a
+ * filter criterion.
  *
- * This is mostly for testing purposes and uses link_cell
- * to get pairs out of the cellsystem by a simple distance
- * criterion.
+ * It uses link_cell to get pairs out of the cellsystem
+ * by a simple distance criterion and
  *
  * Pairs are sorted so that first.id < second.id
  */
-std::vector<std::pair<int, int>> get_pairs(double distance) {
+template <class Filter>
+std::vector<std::pair<int, int>> get_pairs_filtered(double const distance,
+                                                    Filter filter) {
   std::vector<std::pair<int, int>> ret;
+  on_observable_calc();
   auto const cutoff2 = distance * distance;
-
-  cells_update_ghosts(Cells::DATA_PART_POSITION | Cells::DATA_PART_PROPERTIES);
-
-  auto pair_kernel = [&ret, &cutoff2](Particle const &p1, Particle const &p2,
-                                      Distance const &d) {
-    if (d.dist2 < cutoff2)
+  auto pair_kernel = [&ret, &cutoff2, &filter](Particle const &p1,
+                                               Particle const &p2,
+                                               Distance const &d) {
+    if (d.dist2 < cutoff2 and filter(p1) and filter(p2))
       ret.emplace_back(p1.p.identity, p2.p.identity);
   };
 
@@ -79,29 +83,61 @@ std::vector<std::pair<int, int>> get_pairs(double distance) {
   return ret;
 }
 
-void mpi_get_pairs_local(int, int) {
-  on_observable_calc();
+std::vector<std::pair<int, int>> get_pairs(double const distance) {
+  return get_pairs_filtered(distance, [](Particle const &p) { return true; });
+}
 
-  double distance;
-  boost::mpi::broadcast(comm_cart, distance, 0);
+std::vector<std::pair<int, int>>
+get_pairs_of_types(double const distance, std::vector<int> const &types) {
+  return get_pairs_filtered(distance, [types](Particle const &p) {
+    return std::any_of(types.begin(), types.end(),
+                       [p](int const type) { return p.p.type == type; });
+  });
+}
 
+void get_pairs_local(double const distance) {
   auto local_pairs = get_pairs(distance);
-
   Utils::Mpi::gather_buffer(local_pairs, comm_cart);
 }
 
-REGISTER_CALLBACK(mpi_get_pairs_local)
+REGISTER_CALLBACK(get_pairs_local)
 
-std::vector<std::pair<int, int>> mpi_get_pairs(double distance) {
-  mpi_call(mpi_get_pairs_local, 0, 0);
-  on_observable_calc();
+void get_pairs_of_types_local(double const distance,
+                              std::vector<int> const &types) {
+  auto local_pairs = get_pairs_of_types(distance, types);
+  Utils::Mpi::gather_buffer(local_pairs, comm_cart);
+}
 
-  boost::mpi::broadcast(comm_cart, distance, 0);
+REGISTER_CALLBACK(get_pairs_of_types_local)
+
+namespace detail {
+void search_distance_sanity_check(double const distance) {
+  /** get_pairs_filtered() finds pairs via the non_bonded_loop. The maximum
+   *finding range is therefore limited by the decomposition that is used
+   **/
+  auto range = *boost::min_element(cell_structure.max_range());
+  if (distance > range) {
+    runtimeErrorMsg() << "pair search distance " << distance
+                      << " bigger than the decomposition range " << range
+                      << ".\n";
+  }
+}
+} // namespace detail
 
+std::vector<std::pair<int, int>> mpi_get_pairs(double const distance) {
+  detail::search_distance_sanity_check(distance);
+  mpi_call(get_pairs_local, distance);
   auto pairs = get_pairs(distance);
-
   Utils::Mpi::gather_buffer(pairs, comm_cart);
+  return pairs;
+}
 
+std::vector<std::pair<int, int>>
+mpi_get_pairs_of_types(double const distance, std::vector<int> const &types) {
+  detail::search_distance_sanity_check(distance);
+  mpi_call(get_pairs_of_types_local, distance, types);
+  auto pairs = get_pairs_of_types(distance, types);
+  Utils::Mpi::gather_buffer(pairs, comm_cart);
   return pairs;
 }
 
diff --git a/src/core/cells.hpp b/src/core/cells.hpp
index e0c2ad5c52c..7e9db3f52db 100644
--- a/src/core/cells.hpp
+++ b/src/core/cells.hpp
@@ -88,6 +88,14 @@ void cells_update_ghosts(unsigned data_parts);
  */
 std::vector<std::pair<int, int>> mpi_get_pairs(double distance);
 
+/**
+ * @brief Get pairs closer than @p distance if both their types are in @p types
+ *
+ * Pairs are sorted so that first.id < second.id
+ */
+std::vector<std::pair<int, int>>
+mpi_get_pairs_of_types(double distance, std::vector<int> const &types);
+
 /** Check if a particle resorting is required. */
 void check_resort_particles();
 
diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp
index de1bc09c83a..8a1d0735f09 100644
--- a/src/core/integrate.cpp
+++ b/src/core/integrate.cpp
@@ -362,7 +362,7 @@ int python_integrate(int n_steps, bool recalc_forces, bool reuse_forces_par) {
     /* maximal skin that can be used without resorting is the maximal
      * range of the cell system minus what is needed for interactions. */
     skin = std::min(0.4 * max_cut,
-                    *boost::min_element(cell_structure.max_range()) - max_cut);
+                    *boost::min_element(cell_structure.max_cutoff()) - max_cut);
     mpi_bcast_parameter(FIELD_SKIN);
   }
 
diff --git a/src/core/short_range_loop.hpp b/src/core/short_range_loop.hpp
index 2b903b27cc4..b11131d60aa 100644
--- a/src/core/short_range_loop.hpp
+++ b/src/core/short_range_loop.hpp
@@ -38,7 +38,7 @@ struct True {
 template <class BondKernel, class PairKernel,
           class VerletCriterion = detail::True>
 void short_range_loop(BondKernel bond_kernel, PairKernel pair_kernel,
-                      double distance_cutoff = 1.,
+                      double distance_cutoff,
                       const VerletCriterion &verlet_criterion = {}) {
   ESPRESSO_PROFILER_CXX_MARK_FUNCTION;
 
diff --git a/src/core/tuning.cpp b/src/core/tuning.cpp
index b26da82d4d8..538b3d3f984 100644
--- a/src/core/tuning.cpp
+++ b/src/core/tuning.cpp
@@ -105,7 +105,7 @@ void tune_skin(double min_skin, double max_skin, double tol, int int_steps,
    * the maximal range that can be supported by the cell system, but
    * never larger than half the box size. */
   double const max_permissible_skin = std::min(
-      *boost::min_element(cell_structure.max_range()) - maximal_cutoff(),
+      *boost::min_element(cell_structure.max_cutoff()) - maximal_cutoff(),
       0.5 * *boost::max_element(box_geo.length()));
 
   if (adjust_max_skin and max_skin > max_permissible_skin)
diff --git a/src/python/espressomd/cellsystem.pxd b/src/python/espressomd/cellsystem.pxd
index 798f777f89a..7057ed7ed2b 100644
--- a/src/python/espressomd/cellsystem.pxd
+++ b/src/python/espressomd/cellsystem.pxd
@@ -37,7 +37,8 @@ cdef extern from "cells.hpp":
 
     const DomainDecomposition * get_domain_decomposition()
 
-    vector[pair[int, int]] mpi_get_pairs(double distance)
+    vector[pair[int, int]] mpi_get_pairs(double distance) except +
+    vector[pair[int, int]] mpi_get_pairs_of_types(double distance, vector[int] types) except +
     vector[int] mpi_resort_particles(int global_flag)
     void mpi_bcast_cell_structure(int cs)
     void mpi_set_use_verlet_lists(bool use_verlet_lists)
diff --git a/src/python/espressomd/cellsystem.pyx b/src/python/espressomd/cellsystem.pyx
index 48eadba1185..cb7eccfe5d9 100644
--- a/src/python/espressomd/cellsystem.pyx
+++ b/src/python/espressomd/cellsystem.pyx
@@ -23,10 +23,10 @@ from . cimport integrate
 from .globals cimport FIELD_SKIN, FIELD_NODEGRID
 from .globals cimport verlet_reuse, skin
 from .globals cimport mpi_bcast_parameter
+from libcpp.vector cimport vector
 from .cellsystem cimport cell_structure
-from .utils cimport handle_errors
-from .utils import is_valid_type
-from .utils cimport Vector3i
+from .utils cimport handle_errors, Vector3i, check_type_or_throw_except
+
 
 cdef class CellSystem:
     def set_domain_decomposition(self, use_verlet_lists=True):
@@ -109,8 +109,47 @@ cdef class CellSystem:
         self.skin = d['skin']
         self.node_grid = d['node_grid']
 
-    def get_pairs_(self, distance):
-        return mpi_get_pairs(distance)
+    def get_pairs(self, distance, types='all'):
+        """
+        Get pairs of particles closer than threshold value
+
+        Parameters
+        ----------
+        distance : :obj:`float`
+            Pairs of particles closer than ``distance`` are found.
+        types (optional) : list of :obj:`int` or ``'all'``
+            Restrict the pair search to the specified types. Defaults to ``'all'``, in which case all particles are considered.
+
+        Returns
+        -------
+        list of tuples of :obj:`int`
+            The particle pairs identified by their index
+
+        Raises
+        ------
+        Exception
+            If the pair search distance is greater than the cell size
+        """
+        pairs = None
+        if types == 'all':
+            pairs = mpi_get_pairs(distance)
+        else:
+            pairs = self._get_pairs_of_types(distance, types)
+        handle_errors("")
+        return pairs
+
+    def _get_pairs_of_types(self, distance, types):
+        """
+        This function needs to be separated from ``self.get_pairs()`` because ``cdef``
+        cannot be inside an ``else`` block
+        """
+
+        check_type_or_throw_except(
+            types, len(types), int, 'types must be a list of int')
+        cdef vector[int] types_c
+        for type in types:
+            types_c.push_back(type)
+        return mpi_get_pairs_of_types(distance, types_c)
 
     def resort(self, global_flag=True):
         """
diff --git a/testsuite/python/pairs.py b/testsuite/python/pairs.py
index aa9d462b7e5..35459fb1652 100644
--- a/testsuite/python/pairs.py
+++ b/testsuite/python/pairs.py
@@ -15,40 +15,47 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 import unittest as ut
-import unittest_decorators as utx
 import espressomd
 
 
-@utx.skipIfMissingFeatures(["LENNARD_JONES"])
 class PairTest(ut.TestCase):
+    """
+    Tests the particle pair finder of the cell system.
+    It checks that particles are found if their distance is below the threshold, 
+    no matter the type of cell system, periodicity and the image box they are in.
+    Also tests that the ``types`` argument works as expected and an exception is raised
+    when the distance threshold is larger than the cell size.
+    """
+
     s = espressomd.System(box_l=[1.0, 1.0, 1.0])
 
     def setUp(self):
-        self.s.time_step = 0.1
-        self.s.thermostat.turn_off()
-
-        self.s.part.clear()
+        self.s.time_step = 0.1     
         self.s.box_l = 3 * [10.]
         self.s.cell_system.skin = 0.3
 
         # Force an appropriate cell grid
-        self.s.non_bonded_inter[0, 0].lennard_jones.set_params(
-            epsilon=0.0, sigma=1.0, cutoff=1.5, shift=0.0)
+        self.s.min_global_cut = 1.6
 
         vel = [1., 2., 3.]
 
-        self.s.part.add(id=0, pos=[5., 4.5, 5.], v=vel)
-        self.s.part.add(id=1, pos=[5., 5.5, 5.], v=vel)
-        self.s.part.add(id=2, pos=[9.5, 5., 5.], v=vel)
-        self.s.part.add(id=3, pos=[0.5, 5., 5.], v=vel)
-        self.s.part.add(id=4, pos=[5., 5., 9.5], v=vel)
-        self.s.part.add(id=5, pos=[5., 5., 0.5], v=vel)
-        self.s.part.add(id=6, pos=[5., 9.5, 5.], v=vel)
-        self.s.part.add(id=7, pos=[5., 0.5, 5.], v=vel)
-        self.s.part.add(id=8, pos=[5., 9.5, 9.5], v=vel)
-        self.s.part.add(id=9, pos=[5., 0.5, 0.5], v=vel)
-        self.s.part.add(id=10, pos=[1., 1., 1.], v=vel)
-        self.s.part.add(id=11, pos=[9., 9., 9.], v=vel)
+        self.s.part.add(id=0, pos=[5., 4.5, 5.], v=vel, type=0)
+        self.s.part.add(id=1, pos=[5., 5.5, 5.], v=vel, type=1)
+        self.s.part.add(id=2, pos=[9.5, 5., 5.], v=vel, type=2)
+        self.s.part.add(id=3, pos=[0.5, 5., 5.], v=vel, type=3)
+        self.s.part.add(id=4, pos=[5., 5., 9.5], v=vel, type=4)
+        self.s.part.add(id=5, pos=[5., 5., 0.5], v=vel, type=5)
+        self.s.part.add(id=6, pos=[5., 9.5, 5.], v=vel, type=6)
+        self.s.part.add(id=7, pos=[5., 0.5, 5.], v=vel, type=7)
+        self.s.part.add(id=8, pos=[5., 9.5, 9.5], v=vel, type=8)
+        self.s.part.add(id=9, pos=[5., 0.5, 0.5], v=vel, type=9)
+        self.s.part.add(id=10, pos=[1., 1., 1.], v=vel, type=10)
+        self.s.part.add(id=11, pos=[9., 9., 9.], v=vel, type=11)
+
+        self.types_to_get_pairs = [0, 1, 4, 5, 6]
+
+    def tearDown(self):
+        self.s.part.clear()
 
     def expected_pairs(self, periodicity):
         if all(periodicity == (1, 1, 1)):
@@ -56,50 +63,53 @@ def expected_pairs(self, periodicity):
         elif all(periodicity == (1, 1, 0)):
             return [(0, 1), (2, 3), (6, 7)]
 
-    def check(self):
-        pairs = self.s.cell_system.get_pairs_(1.5)
+    def expected_pairs_with_types(self, periodicity):
+        if all(periodicity == (1, 1, 1)):
+            return [(0, 1), (4, 5)]
+        elif all(periodicity == (1, 1, 0)):
+            return [(0, 1)]
+
+    def check_pairs(self):
+        pairs = self.s.cell_system.get_pairs(1.5)
         epairs = self.expected_pairs(self.s.periodicity)
-        print(pairs)
+        self.assertSetEqual(set(pairs), set(epairs))
+
+        pairs_by_type = self.s.cell_system.get_pairs(
+            1.5, types=self.types_to_get_pairs)
+        epairs_by_type = self.expected_pairs_with_types(self.s.periodicity)
+        self.assertSetEqual(set(pairs_by_type), set(epairs_by_type))
+
+    def check_exception(self):
+        with self.assertRaises(Exception):
+            self.s.cell_system.get_pairs(3.)
 
-        self.assertEqual(len(pairs), len(epairs))
-        for p in pairs:
-            self.assertIn(p, epairs)
+    def run_and_check(self, n_steps=100):
+        self.s.integrator.run(0)
+        self.check_pairs()
+        self.s.integrator.run(n_steps)
+        self.check_pairs()
 
     def test_nsquare(self):
         self.s.cell_system.set_n_square()
         self.s.periodicity = [1, 1, 1]
-
-        self.s.integrator.run(0)
-        self.check()
-        self.s.integrator.run(100)
-        self.check()
+        self.run_and_check()
 
     def test_nsquare_partial_z(self):
         self.s.cell_system.set_n_square()
         self.s.periodicity = [1, 1, 0]
-
-        self.s.integrator.run(0)
-        self.check()
-        self.s.integrator.run(100)
-        self.check()
+        self.run_and_check()
 
     def test_dd(self):
         self.s.cell_system.set_domain_decomposition()
         self.s.periodicity = [1, 1, 1]
-
-        self.s.integrator.run(0)
-        self.check()
-        self.s.integrator.run(100)
-        self.check()
+        self.run_and_check()
+        self.check_exception()
 
     def test_dd_partial_z(self):
         self.s.cell_system.set_domain_decomposition()
         self.s.periodicity = [1, 1, 0]
-
-        self.s.integrator.run(0)
-        self.check()
-        self.s.integrator.run(100)
-        self.check()
+        self.run_and_check()
+        self.check_exception()
 
 
 if __name__ == "__main__":
diff --git a/testsuite/python/random_pairs.py b/testsuite/python/random_pairs.py
index 748d5dff3be..450d027d6a6 100644
--- a/testsuite/python/random_pairs.py
+++ b/testsuite/python/random_pairs.py
@@ -58,7 +58,7 @@ def tearDown(self):
 
     def pairs_n2(self, dist):
         # Go through list of all possible pairs for full periodicity
-        # and skip those that ar not within the desired distance
+        # and skip those that are not within the desired distance
         # for the current periodicity
 
         pairs = []
@@ -73,7 +73,7 @@ def check_duplicates(self, l):
             self.assertEqual(e, 1)
 
     def check_pairs(self, n2_pairs):
-        cs_pairs = self.system.cell_system.get_pairs_(1.5)
+        cs_pairs = self.system.cell_system.get_pairs(1.5)
         self.check_duplicates(cs_pairs)
         self.assertGreater(len(cs_pairs), 0)
         self.assertEqual(n2_pairs ^ set(cs_pairs), set())