diff --git a/src/pairs/sim/read_from_file.py b/src/pairs/sim/read_from_file.py
index 7f37b5fcf4263e53cc713793d1b6ec516e4032b8..1d3d4edd5356cb71032cf5f488ba70e16539ed6e 100644
--- a/src/pairs/sim/read_from_file.py
+++ b/src/pairs/sim/read_from_file.py
@@ -1,5 +1,5 @@
 from pairs.ir.block import pairs_inline
-from pairs.ir.functions import Call_Int
+from pairs.ir.functions import Call_Int, Call_Void
 from pairs.ir.properties import PropertyList
 from pairs.ir.types import Types
 from pairs.sim.grid import MutableGrid
@@ -16,9 +16,13 @@ class ReadFromFile(Lowerable):
 
     @pairs_inline
     def lower(self):
-        self.sim.nlocal.set(Call_Int(self.sim, "pairs::read_particle_data",
-                            [self.filename, self.grid_buffer, self.props, self.props.length()]))
-
+        Call_Void(self.sim, "pairs::read_grid_data", [self.filename, self.grid_buffer])
         for d in range(self.sim.ndims()):
             self.grid.min(d).set(self.grid_buffer[d * 2 + 0])
             self.grid.max(d).set(self.grid_buffer[d * 2 + 1])
+
+        dom_part = self.sim.domain_partitioning()
+        grid_array = [[self.grid.min(d), self.grid.max(d)] for d in range(self.sim.ndims())]
+        Call_Void(self.sim, "pairs->initDomain", [param for delim in grid_array for param in delim]),
+        Call_Void(self.sim, "pairs->fillCommunicationArrays", [dom_part.neighbor_ranks, dom_part.pbc, dom_part.subdom])
+        self.sim.nlocal.set(Call_Int(self.sim, "pairs::read_particle_data", [self.filename, self.props, self.props.length()]))
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index ca7d56721a606af3de12514146854cbf4ea2275b..2ed480c5566493b5dbe109093bb1ba79ae4ffe2f 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -62,6 +62,7 @@ class Simulation:
         self.iter_id = 0
         self.vtk_file = None
         self._target = None
+        self._dom_part = DimensionRanges(self)
         self.nparticles = self.nlocal + self.nghost
         self.properties.add_capacity(self.particle_capacity)
 
@@ -236,20 +237,12 @@ class Simulation:
     def cell_spacing(self):
         return self.cell_lists.cutoff_radius
 
+    def domain_partitioning(self):
+        return self._dom_part
+
     def generate(self):
         assert self._target is not None, "Target not specified!"
-
-        dom_part = DimensionRanges(self)
-        comm = Comm(self, dom_part)
-
-        self.capture_statements(False)
-        grid_array = [[self.grid.min(d), self.grid.max(d)] for d in range(self.ndims())]
-        self.setups.add_statement([
-            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])
-        ])
-
-        self.capture_statements() # TODO: check if this is actually required
+        comm = Comm(self, self._dom_part)
 
         timestep = Timestep(self, self.ntimesteps, [
             (comm.exchange(), 20),