diff --git a/Makefile b/Makefile
index 18bad17bd043138c6056c4f667cd0c9c6df860de..29899a6abca070565ec5c04dc2cd36bee7b49ab0 100644
--- a/Makefile
+++ b/Makefile
@@ -17,7 +17,8 @@ lj_ns.cu:
 
 # Targets
 lj_cpu: lj_ns.cpp
-	g++ -o lj_cpu lj_ns.cpp
+#	g++ -o lj_cpu lj_ns.cpp -DDEBUG
+	mpic++ -o lj_cpu lj_ns.cpp -DDEBUG
 
 lj_gpu: lj_ns.cu
 	nvcc -o lj_gpu lj_ns.cu
diff --git a/runtime/comm.hpp b/runtime/comm.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d9a7abb240a866ee46cda4812408953f66f2e576
--- /dev/null
+++ b/runtime/comm.hpp
@@ -0,0 +1,31 @@
+#include "pairs.hpp"
+
+#pragma once
+
+namespace pairs {
+
+template<int ndims>
+void initDomain(PairsSimulation<ndims> *pairs, real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax) {
+    pairs->initDomain(xmin, xmax, ymin, ymax, zmin, zmax);
+}
+
+template<int ndims>
+void communicateSizes(PairsSimulation<ndims> *pairs, int dim, const int *send_sizes, int *recv_sizes) {
+    pairs->getDomainPartitioner()->communicateSizes(dim, send_sizes, recv_sizes);
+}
+
+template<int ndims>
+void communicateData(
+    PairsSimulation<ndims> *pairs, int dim, int elem_size,
+    const real_t *send_buf, const int *send_offsets, const int *nsend,
+    real_t *recv_buf, const int *recv_offsets, const int *nrecv) {
+
+    pairs->getDomainPartitioner()->communicateData(dim, elem_size, send_buf, send_offsets, nsend, recv_buf, recv_offsets, nrecv);
+}
+
+template<int ndims>
+void fillCommunicationArrays(PairsSimulation<ndims> *pairs, int neighbor_ranks[], int pbc[], real_t subdom[]) {
+    pairs->getDomainPartitioner()->fillArrays(neighbor_ranks, pbc, subdom);
+}
+
+}
diff --git a/runtime/domain/domain_partitioning.hpp b/runtime/domain/domain_partitioning.hpp
index bf43df21d2bfccc3be0e2409a5d1449142cf72e6..3936a506bdc702dda79d1263514c51ff92d667da 100644
--- a/runtime/domain/domain_partitioning.hpp
+++ b/runtime/domain/domain_partitioning.hpp
@@ -1,24 +1,24 @@
 #include <mpi.h>
-//---
-#include "pairs.hpp"
 
 #pragma once
 
+typedef double real_t;
+
 namespace pairs {
 
 template<int ndims>
-class DomainPartitioning {
+class DomainPartitioner {
 protected:
     int world_size, rank;
     real_t grid_min[ndims];
     real_t grid_max[ndims];
     
 public:
-    DomainPartitioning(real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax) {
-        static_assert(ndims == 3, "DomainPartitioning(): number of dimensions mismatching!");
-        PAIRS_ASSERT(xmax > xmin);
-        PAIRS_ASSERT(ymax > ymin);
-        PAIRS_ASSERT(zmax > zmin);
+    DomainPartitioner(real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax) {
+        static_assert(ndims == 3, "DomainPartitioner(): number of dimensions mismatching!");
+        //PAIRS_ASSERT(xmax > xmin);
+        //PAIRS_ASSERT(ymax > ymin);
+        //PAIRS_ASSERT(zmax > zmin);
 
         grid_min[0] = xmin;
         grid_max[0] = xmax;
@@ -29,13 +29,18 @@ public:
     }
 
     virtual void initialize();
-    virtual void communicateSizes(const real_t *nsend, const real_t *nrecv);
-    virtual void communicateData(const real_t *send_buf, size_t *nsend, const real_t *recv_buf, size_t *nrecv, size_t elem_size);
+    virtual void fillArrays(int neighbor_ranks[], int pbc[], real_t subdom[]);
+    virtual void communicateSizes(int dim, const int *nsend, int *nrecv);
+    virtual void communicateData(
+        int dim, int elem_size,
+        const real_t *send_buf, const int *send_offsets, const int *nsend,
+        real_t *recv_buf, const int *recv_offsets, const int *nrecv);
     virtual void finalize();
+    inline int getRank() { return rank; }
 };
 
 template<int ndims>
-class DimensionRanges : DomainPartitioning<ndims> {
+class DimensionRanges : DomainPartitioner<ndims> {
 protected:
     int nranks[ndims];
     int prev[ndims];
@@ -46,6 +51,9 @@ protected:
     real_t subdom_max[ndims];
 
 public:
+    DimensionRanges(real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax) :
+        DomainPartitioner<ndims>(xmin, xmax, ymin, ymax, zmin, zmax) {}
+
     void fillArrays(int neighbor_ranks[], int pbc[], real_t subdom[]) {
         for(int d = 0; d < ndims; d++) {
             neighbor_ranks[d * 2 + 0] = this->prev[d];
@@ -57,28 +65,28 @@ public:
         }
     }
 
-    void communicateSizes(int dim, const real_t *send_sizes, const real_t *recv_sizes) {
-        if(prev[dim] != rank) {
+    void communicateSizes(int dim, const int *send_sizes, int *recv_sizes) {
+        if(prev[dim] != this->getRank()) {
             MPI_Send(&send_sizes[dim * 2 + 0], 1, MPI_INT, prev[dim], 0, MPI_COMM_WORLD);
             MPI_Recv(&recv_sizes[dim * 2 + 0], 1, MPI_INT, next[dim], 0, MPI_COMM_WORLD);
         } else {
-            recv_sizes[d * 2 + 0] = send_sizes[dim * 2 + 0];
+            recv_sizes[dim * 2 + 0] = send_sizes[dim * 2 + 0];
         }
 
-        if(next[dim] != rank) {
+        if(next[dim] != this->getRank()) {
             MPI_Send(&send_sizes[dim * 2 + 1], 1, MPI_INT, next[dim], 0, MPI_COMM_WORLD);
             MPI_Recv(&recv_sizes[dim * 2 + 1], 1, MPI_INT, prev[dim], 0, MPI_COMM_WORLD);
         } else {
-            recv_sizes[dim * 2 + 1] = send_sizes[d * 2 + 1];
+            recv_sizes[dim * 2 + 1] = send_sizes[dim * 2 + 1];
         }
     }
 
     void communicateData(
-        int dim, size_t elem_size,
-        const real_t *send_buf, const size_t *send_offsets, const size_t *nsend,
-        const real_t *recv_buf, const size_t *recv_offsets, const size_t *nrecv) {
+        int dim, int elem_size,
+        const real_t *send_buf, const int *send_offsets, const int *nsend,
+        real_t *recv_buf, const int *recv_offsets, const int *nrecv) {
 
-        if(prev[dim] != rank) {
+        if(prev[dim] != this->getRank()) {
             MPI_Send(&send_buf[send_offsets[dim * 2 + 0]], nsend[dim * 2 + 0] * elem_size, MPI_DOUBLE, prev[dim], 0, MPI_COMM_WORLD);
             MPI_Recv(&recv_buf[recv_offsets[dim * 2 + 0]], nrecv[dim * 2 + 0] * elem_size, MPI_DOUBLE, next[dim], 0, MPI_COMM_WORLD);
         } else {
@@ -87,7 +95,7 @@ public:
             }
         }
 
-        if(next[dim] != rank) {
+        if(next[dim] != this->getRank()) {
             MPI_Send(&send_buf[send_offsets[dim * 2 + 1]], nsend[dim * 2 + 1] * elem_size, MPI_DOUBLE, next[dim], 0, MPI_COMM_WORLD);
             MPI_Recv(&recv_buf[recv_offsets[dim * 2 + 1]], nrecv[dim * 2 + 1] * elem_size, MPI_DOUBLE, prev[dim], 0, MPI_COMM_WORLD);
         } else {
@@ -98,6 +106,7 @@ public:
     }
 };
 
-class ListOfBoxes : DomainPartitioning {};
+template<int ndims>
+class ListOfBoxes : DomainPartitioner<ndims> {};
 
 }
diff --git a/runtime/pairs.hpp b/runtime/pairs.hpp
index a79fea8167dedb88342ea83017a282aa01d73e1e..e0d00da69d767f0adb58b0dbf6bd7b6e7cd2b366 100644
--- a/runtime/pairs.hpp
+++ b/runtime/pairs.hpp
@@ -10,15 +10,19 @@
 #   include "devices/dummy.hpp"
 #endif
 
+#include "domain/domain_partitioning.hpp"
+
 #pragma once
 
 #ifdef DEBUG
 #   include <assert.h>
 #   define PAIRS_DEBUG(...)     fprintf(stderr, __VA_ARGS__)
 #   define PAIRS_ASSERT(a)      assert(a)
+#   define PAIRS_EXCEPTION(a)
 #else
 #   define PAIRS_DEBUG(...)
 #   define PAIRS_ASSERT(a)
+#   define PAIRS_EXCEPTION(a)
 #endif
 
 #define PAIRS_ERROR(...)    fprintf(stderr, __VA_ARGS__)
@@ -42,6 +46,11 @@ enum DataLayout {
     SoA
 };
 
+enum DomainPartitioning {
+    DimRanges = 0,
+    BoxList,
+};
+
 template<typename real_t>
 class Vector3 {
 private:
@@ -243,23 +252,40 @@ public:
     }
 };
 
+template<int ndims>
 class PairsSimulation {
 private:
+    DomainPartitioner<ndims> *dom_part;
     std::vector<Property> properties;
     std::vector<Array> arrays;
     DeviceFlags *prop_flags, *array_flags;
+    DomainPartitioning dom_part_type;
     int nprops, narrays;
 public:
-    PairsSimulation(int nprops_, int narrays_) {
+    PairsSimulation(int nprops_, int narrays_, DomainPartitioning dom_part_type_) {
+        dom_part_type = dom_part_type_;
         prop_flags = new DeviceFlags(nprops_);
         array_flags = new DeviceFlags(narrays_);
     }
 
     ~PairsSimulation() {
+        dom_part->finalize();
         delete prop_flags;
         delete array_flags;
     }
 
+    void initDomain(real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax) {
+        if(dom_part_type == DimRanges) {
+            dom_part = new DimensionRanges<ndims>(xmin, xmax, ymin, ymax, zmin, zmax);
+        } else {
+            PAIRS_EXCEPTION("Domain partitioning type not implemented!\n");
+        }
+
+        dom_part->initialize();
+    }
+
+    DomainPartitioner<ndims> *getDomainPartitioner() { return dom_part; }
+
     template<typename T_ptr>
     void addArray(array_t id, std::string name, T_ptr **h_ptr, std::nullptr_t, size_t size) {
         PAIRS_ASSERT(size > 0);
diff --git a/runtime/read_from_file.hpp b/runtime/read_from_file.hpp
index 3376b132125b12c8e1f9b0d524d045ef760a7b8b..2404e9bb22105c1825c15efa838217d1f052fa4e 100644
--- a/runtime/read_from_file.hpp
+++ b/runtime/read_from_file.hpp
@@ -9,7 +9,8 @@
 
 namespace pairs {
 
-size_t read_particle_data(PairsSimulation *ps, const char *filename, double *grid_buffer, const property_t properties[], size_t nprops) {
+template<int ndims>
+size_t read_particle_data(PairsSimulation<ndims> *ps, const char *filename, double *grid_buffer, const property_t properties[], size_t nprops) {
     std::ifstream in_file(filename, std::ifstream::in);
     std::string line;
     size_t n = 0;
@@ -23,9 +24,7 @@ size_t read_particle_data(PairsSimulation *ps, const char *filename, double *gri
 
             while(std::getline(line_stream, in0, ',')) {
                 if(!read_grid_data) {
-                    // TODO: store this from PairsSimulation class
-                    const int num_dims = 3;
-                    PAIRS_ASSERT(i < num_dims * 2);
+                    PAIRS_ASSERT(i < ndims * 2);
                     grid_buffer[i] = std::stod(in0);
                 } else {
                     PAIRS_ASSERT(i < nprops);
diff --git a/runtime/vtk.hpp b/runtime/vtk.hpp
index 33d598af1bb63400851d5dff4f0ef3ddd1c76958..894fc119a6b57444ab2a076c66069f15643be788 100644
--- a/runtime/vtk.hpp
+++ b/runtime/vtk.hpp
@@ -8,7 +8,8 @@
 
 namespace pairs {
 
-void vtk_write_data(PairsSimulation *ps, const char *filename, int start, int end, int timestep) {
+template <int ndims>
+void vtk_write_data(PairsSimulation<ndims> *ps, const char *filename, int start, int end, int timestep) {
     std::string output_filename(filename);
     std::ostringstream filename_oss;
     filename_oss << filename << "_" << timestep << ".vtk";
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index e7e5fec49692bd3481ad5ed11675000738a2fe19..68bab60c5de2e45b9636d894aecc7c72f3ff333c 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -55,6 +55,7 @@ class CGen:
         self.print("#include <stdlib.h>")
         self.print("//---")
         self.print("#include \"runtime/pairs.hpp\"")
+        self.print("#include \"runtime/comm.hpp\"")
         self.print("#include \"runtime/read_from_file.hpp\"")
         self.print("#include \"runtime/vtk.hpp\"")
 
@@ -87,10 +88,11 @@ class CGen:
 
     def generate_module(self, module):
         if module.name == 'main':
+            ndims = module.sim.ndims()
             nprops = module.sim.properties.nprops()
             narrays = module.sim.arrays.narrays()
             self.print("int main() {")
-            self.print(f"    PairsSimulation *pairs = new PairsSimulation({nprops}, {narrays});")
+            self.print(f"    PairsSimulation<{ndims}> *pairs = new PairsSimulation<{ndims}>({nprops}, {narrays}, DimRanges);")
             self.generate_statement(module.block)
             self.print("    return 0;")
             self.print("}")
diff --git a/src/pairs/sim/comm.py b/src/pairs/sim/comm.py
index fc88339b03969e5e6e2678921b239db7dd0ca19f..8ce93a477acd74ada38c8d58d20ec31a081cd467 100644
--- a/src/pairs/sim/comm.py
+++ b/src/pairs/sim/comm.py
@@ -3,6 +3,7 @@ from pairs.ir.bin_op import BinOp
 from pairs.ir.block import pairs_device_block, pairs_host_block, pairs_inline
 from pairs.ir.branches import Branch, Filter
 from pairs.ir.cast import Cast
+from pairs.ir.functions import Call_Void
 from pairs.ir.loops import For, ParticleFor, While
 from pairs.ir.utils import Print
 from pairs.ir.select import Select
@@ -72,10 +73,11 @@ class CommunicateSizes(Lowerable):
         super().__init__(comm.sim)
         self.comm = comm
         self.step = step
+        self.sim.add_statement(self)
 
     @pairs_inline
     def lower(self):
-        Call_Void(self.sim, "pairs->communicateSizes", [self.step, self.comm.nsend, self.comm.nrecv])
+        Call_Void(self.sim, "pairs::communicateSizes", [self.step, self.comm.nsend, self.comm.nrecv])
 
 
 class CommunicateData(Lowerable):
@@ -84,11 +86,12 @@ class CommunicateData(Lowerable):
         self.comm = comm
         self.step = step
         self.prop_list = prop_list
+        self.sim.add_statement(self)
 
     @pairs_inline
     def lower(self):
         elem_size = sum([self.sim.ndims() if p.type() == Types.Vector else 1 for p in self.prop_list])
-        Call_Void(self.sim, "pairs->communicateData", [self.step, elem_size,
+        Call_Void(self.sim, "pairs::communicateData", [self.step, elem_size,
                                                        self.comm.send_buffer, self.comm.send_offsets, self.comm.nsend,
                                                        self.comm.recv_buffer, self.comm.recv_offsets, self.comm.nrecv])
 
diff --git a/src/pairs/sim/domain_partitioning.py b/src/pairs/sim/domain_partitioning.py
index d317c516f2a87ceb701a82742f1d55c386acc9bc..a26ea2b9a171aafb947ec922d28a39a325e5710a 100644
--- a/src/pairs/sim/domain_partitioning.py
+++ b/src/pairs/sim/domain_partitioning.py
@@ -12,7 +12,7 @@ class DimensionRanges:
         self.sim = sim
         self.neighbor_ranks = sim.add_static_array('neighbor_ranks', [sim.ndims() * 2], Types.Int32)
         self.pbc            = sim.add_static_array('pbc', [sim.ndims() * 2], Types.Int32)
-        self.subdom         = sim.add_static_array('subdom', [sim.ndims() * 2], Types.Int32)
+        self.subdom         = sim.add_static_array('subdom', [sim.ndims() * 2], Types.Double)
 
     def number_of_steps(self):
         return self.sim.ndims()
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index a7d9edf6659aad7138c4cc97312c297bb6582a0d..f03b9f653aff09471d253377d8f459e0bf70218c 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -1,6 +1,7 @@
 from pairs.ir.arrays import Arrays
 from pairs.ir.block import Block
 from pairs.ir.branches import Filter
+from pairs.ir.functions import Call_Void
 from pairs.ir.kernel import Kernel
 from pairs.ir.layouts import Layouts
 from pairs.ir.module import Module
@@ -153,6 +154,7 @@ class Simulation:
         self.setups.add_statement(read_object)
         self.grid = read_object.grid
 
+
     def build_cell_lists(self, spacing):
         self.cell_lists = CellLists(self, self.grid, spacing, spacing)
         return self.cell_lists
@@ -249,8 +251,11 @@ class Simulation:
         timestep.add(VTKWrite(self, self.vtk_file, timestep.timestep() + 1))
         self.leave()
 
+        grid_array = [[self.grid.min(d), self.grid.max(d)] for d in range(self.ndims())]
         body = Block.from_list(self, [
             self.setups,
+            Call_Void(self, "pairs::initDomain", [param for delim in grid_array for param in delim]),
+            Call_Void(self, "pairs::fillCommunicationArrays", [dom_part.neighbor_ranks, dom_part.pbc, dom_part.subdom]),
             CellListsStencilBuild(self, self.cell_lists),
             VTKWrite(self, self.vtk_file, 0),
             timestep.as_block()