diff --git a/code_gen/cgen.py b/code_gen/cgen.py
index b39bc2f1548a87744b6bd605edec3a109f44043e..f1a913b649402b2791a47686b384875ebab8ee7b 100644
--- a/code_gen/cgen.py
+++ b/code_gen/cgen.py
@@ -1,15 +1,16 @@
 from ir.assign import Assign
-from ir.arrays import ArrayAccess, ArrayDecl
+from ir.arrays import Array, ArrayAccess, ArrayDecl
 from ir.block import Block
 from ir.branches import Branch
 from ir.cast import Cast
 from ir.bin_op import BinOp, BinOpDef
-from ir.data_types import Type_Int, Type_Float, Type_Vector
+from ir.data_types import Type_Int, Type_Float, Type_String, Type_Vector
+from ir.functions import Call
 from ir.lit import Lit
 from ir.loops import For, Iter, ParticleFor, While
 from ir.math import Sqrt
 from ir.memory import Malloc, Realloc
-from ir.properties import Property
+from ir.properties import Property, PropertyList
 from ir.select import Select
 from ir.sizeof import Sizeof
 from ir.utils import Print
@@ -20,6 +21,8 @@ from code_gen.printer import Printer
 
 
 class CGen:
+    temp_id = 0
+
     def type2keyword(type_):
         return (
             'double' if type_ == Type_Float or type_ == Type_Vector
@@ -100,6 +103,10 @@ class CGen:
 
             self.print("}") 
 
+        if isinstance(ast_node, Call):
+            call = self.generate_expression(ast_node)
+            self.print("{call};")
+
         if isinstance(ast_node, For):
             iterator = self.generate_expression(ast_node.iterator)
             lower_range = None
@@ -161,6 +168,9 @@ class CGen:
             self.print("}")
 
     def generate_expression(self, ast_node, mem=False, index=None):
+        if isinstance(ast_node, Array):
+            return ast_node.name()
+
         if isinstance(ast_node, ArrayAccess):
             index = self.generate_expression(ast_node.index)
             array_name = ast_node.array.name()
@@ -202,6 +212,10 @@ class CGen:
 
             return f"e{ast_node.id()}"
 
+        if isinstance(ast_node, Call):
+            params = ", ".join(["ps"] + [str(self.generate_expression(p)) for p in ast_node.parameters()])
+            return f"{ast_node.name()}({params})"
+
         if isinstance(ast_node, Cast):
             tkw = CGen.type2keyword(ast_node.cast_type)
             expr = self.generate_expression(ast_node.expr)
@@ -213,11 +227,23 @@ class CGen:
 
         if isinstance(ast_node, Lit):
             assert mem is False, "Literal is not lvalue!"
+
+            if ast_node.type() == Type_String:
+                return f"\"{ast_node.value}\""
+
             return ast_node.value
 
         if isinstance(ast_node, Property):
             return ast_node.name()
 
+        if isinstance(ast_node, PropertyList):
+            tid = CGen.temp_id
+            list_ref = f"prop_list_{tid}"
+            list_def = ", ".join(str(p.id) for p in ast_node)
+            self.print(f"const int {list_ref}[] = {{{list_def}}};")
+            CGen.temp_id += 1
+            return list_ref
+
         if isinstance(ast_node, Sizeof):
             assert mem is False, "Sizeof expression is not lvalue!"
             tkw = CGen.type2keyword(ast_node.data_type)
diff --git a/ir/assign.py b/ir/assign.py
index 388e79c41d29286853a10def0b636cd054804154..997cd4c8bb40546b68708ec517737f52efeee3e1 100644
--- a/ir/assign.py
+++ b/ir/assign.py
@@ -13,7 +13,7 @@ class Assign(ASTNode):
         if dest.type() == Type_Vector:
             self.assignments = []
 
-            for i in range(0, sim.dimensions):
+            for i in range(0, sim.ndims()):
                 from ir.bin_op import BinOp
                 dim_src = src if not isinstance(src, BinOp) or src.type() != Type_Vector else src[i]
                 self.assignments.append((dest[i], dim_src))
diff --git a/ir/data_types.py b/ir/data_types.py
index d7b2a4c07884cc3add6daf28092cd2fa319d7780..e3dc0ad6eef799b5d7165c86b032aeaad83eb8a8 100644
--- a/ir/data_types.py
+++ b/ir/data_types.py
@@ -2,5 +2,6 @@ Type_Invalid = -1
 Type_Int = 0
 Type_Float = 1
 Type_Bool = 2
-Type_Vector = 3
-Type_Array = 4
+Type_String = 3
+Type_Vector = 4
+Type_Array = 5
diff --git a/ir/functions.py b/ir/functions.py
new file mode 100644
index 0000000000000000000000000000000000000000..4ae220deff22b9ec30326a0899cbec166a4f1f44
--- /dev/null
+++ b/ir/functions.py
@@ -0,0 +1,26 @@
+from ir.bin_op import ASTTerm
+from ir.data_types import Type_Int
+from ir.lit import as_lit_ast
+
+
+class Call(ASTTerm):
+    def __init__(self, sim, func_name, params, return_type):
+        super().__init__(sim)
+        self.func_name = func_name
+        self.params = [as_lit_ast(sim, p) for p in params]
+        self.return_type = return_type
+
+    def name(self):
+        return self.func_name
+
+    def parameters(self):
+        return self.params
+
+    def type(self):
+        return self.return_type
+
+
+class Call_Int(Call):
+    def __init__(self, sim, func_name, parameters):
+        super().__init__(sim, func_name, parameters, Type_Int)
+        
diff --git a/ir/lit.py b/ir/lit.py
index 89a55feb0e194e1d4112b5e8f355a162297c6dfb..ab6215714fa3e2f287d36c7c48a7398e52fe4b91 100644
--- a/ir/lit.py
+++ b/ir/lit.py
@@ -1,9 +1,9 @@
 from ir.ast_node import ASTNode
-from ir.data_types import Type_Invalid, Type_Int, Type_Float, Type_Bool, Type_Vector
+from ir.data_types import Type_Invalid, Type_Int, Type_Float, Type_Bool, Type_String, Type_Vector
 
 
 def is_literal(a):
-    return isinstance(a, (int, float, bool, list))
+    return isinstance(a, (int, float, bool, str, list))
 
 
 def as_lit_ast(sim, a):
@@ -25,6 +25,9 @@ class Lit(ASTNode):
         if isinstance(value, bool):
             self.lit_type = Type_Bool
 
+        if isinstance(value, str):
+            self.lit_type = Type_String
+
         if isinstance(value, list):
             self.lit_type = Type_Vector
 
diff --git a/ir/loops.py b/ir/loops.py
index f035afa708c6d45797e69bfde4d2aea55185ca18..0243dcf4134035d9bb543aee7ba6f420cc0c456b 100644
--- a/ir/loops.py
+++ b/ir/loops.py
@@ -116,9 +116,7 @@ class NeighborFor():
             cl = self.cell_lists
             for s in For(self.sim, 0, cl.nstencil):
                 neigh_cell = cl.particle_cell[self.particle] + cl.stencil[s]
-                for _ in Filter(self.sim,
-                                BinOp.and_op(neigh_cell >= 0,
-                                             neigh_cell <= cl.ncells_all)):
+                for _ in Filter(self.sim, BinOp.and_op(neigh_cell >= 0, neigh_cell <= cl.ncells)):
                     for nc in For(self.sim, 0, cl.cell_sizes[neigh_cell]):
                         it = cl.cell_particles[neigh_cell][nc]
                         for _ in Filter(self.sim, BinOp.neq(it, self.particle)):
diff --git a/ir/properties.py b/ir/properties.py
index 3adc4e28bccd82aff8bca4e8ccf6e36b350c49d6..3b5635c4a1bb9bc3cb97bdfe4444077409b2b07e 100644
--- a/ir/properties.py
+++ b/ir/properties.py
@@ -40,17 +40,24 @@ class Properties:
 
 
 class Property(ASTNode):
+    last_prop_id = 0
+
     def __init__(self, sim, name, dtype, default, volatile, layout=Layout_AoS):
         super().__init__(sim)
+        self.id = Property.last_prop_id
         self.prop_name = name
         self.prop_type = dtype
         self.prop_layout = layout
         self.default_value = default
         self.volatile = volatile
+        Property.last_prop_id += 1
 
     def __str__(self):
         return f"Property<{self.prop_name}>"
 
+    def id(self):
+        return self.id
+
     def name(self):
         return self.prop_name
 
@@ -69,3 +76,21 @@ class Property(ASTNode):
     def __getitem__(self, expr):
         from ir.bin_op import BinOp
         return BinOp(self.sim, self, expr, '[]', True)
+
+
+class PropertyList(ASTNode):
+    def __init__(self, sim, properties_list):
+        super().__init__(sim)
+        self.list = []
+        for p in properties_list:
+            if isinstance(p, Property):
+                self.list.append(p)
+
+            if isinstance(p, str):
+                self.list.append(sim.prop(p))
+
+    def __iter__(self):
+        yield from self.list
+
+    def length(self):
+        return len(self.list)
diff --git a/ir/transform.py b/ir/transform.py
index a9c5eb00b4b1359e28df55e1f3c3147ae911c6fc..ba738b6070bcbad0ed4090fa5a285c4e6f75725e 100644
--- a/ir/transform.py
+++ b/ir/transform.py
@@ -23,7 +23,7 @@ class Transform:
                     flat_index = None
 
                     if layout == Layout_AoS:
-                        flat_index = ast.rhs * ast.sim.dimensions + i
+                        flat_index = ast.rhs * ast.sim.ndims() + i
 
                     elif layout == Layout_SoA:
                         flat_index = i * ast.sim.particle_capacity + ast.rhs
diff --git a/runtime/pairs.hpp b/runtime/pairs.hpp
index 6a356cb30566e4a284774d51726ff3c0bdf63c48..57c346cb65ab24d5e9e62c891a29298245f721c6 100644
--- a/runtime/pairs.hpp
+++ b/runtime/pairs.hpp
@@ -2,57 +2,107 @@
 #include <vector>
 
 #define PAIRS_ASSERT(a)
+#define PAIRS_ERROR(a)
 
 namespace pairs {
 
 typedef int property_t;
 typedef int layout_t;
 
+enum PropertyType {
+    Prop_Invalid = -1,
+    Prop_Integer = 0,
+    Prop_Float,
+    Prop_Vector
+}
+
+enum DataLayout {
+    Invalid = -1,
+    AoS = 0,
+    SoA
+}
+
 class PairsSim {
     PairsSim() : {}
 
 public:
-    Property getProperty(property_t id) {
-        return std::find_if(properties.begin(), properties.end(), [id](Property p) { return p.getId() == id; });
-    }
+    void addProperty(Property prop) { properties.push_back(prop); }
 
-    Property getPropertyByName(std::string name) {
-        return std::find_if(properties.begin(), properties.end(), [name](Property p) { return p.getName() == name; });
-    }
+    Property *getProperty(property_t id) {
+        auto p = std::find_if(properties.begin(), properties.end(), [id](Property p) {
+            return p.getId() == id;
+        });
 
-    VectorPtr getFloatPropertyMutablePtr(property_t property) {
-        return static_cast<VectorPtr>(getProperty(property).getPointer());
+        return (p != std::end(properties)) ? p : nullptr;
     }
 
-    int *getIntegerPropertyMutablePtr(property_t property) {
-        return static_cast<int *>(getProperty(property).getPointer());
-    }
+    Property *getPropertyByName(std::string name) {
+        auto p = std::find_if(properties.begin(), properties.end(), [name](Property p) {
+            return p.getName() == name;
+        });
 
-    double *getFloatPropertyMutablePtr(property_t property) {
-        return static_cast<double *>(getProperty(property).getPointer());
+        return (p != std::end(properties)) ? p : nullptr;
     }
 
+    VectorPtr getVectorProperty(property_t property) { return static_cast<VectorPtr>(getProperty(property)->getPointer()); }
+    int *getIntegerPropertyPtr(property_t property) { return static_cast<int *>(getProperty(property)->getPointer()); }
+    double *getFloatPropertyPtr(property_t property) { return static_cast<double *>(getProperty(property)->getPointer()); }
+
 private:
     std::vector<Property> properties;
 };
 
 class Property {
-    Property(property_t id_, std::string name_) : id(id_), name(name_), layout(-1), ptr(nullptr) {}
-    Property(property_t id_, std::string name_, layout_t layout_) : id(id_), name(name_), layout(layout_), ptr(nullptr) {}
-    Property(property_t id_, std::string name_, void *ptr_) : id(id_), name(name_), layout(-1), ptr(ptr_) {}
-    Property(property_t id_, std::string name_, layout_t layout_, void *ptr_) : id(id_), name(name_), layout(layout_), ptr(ptr_) {}
+    Property(property_t id_, std::string name_, void *ptr_, PropertyType type_) :
+        id(id_),
+        name(name_),
+        ptr(ptr_),
+        type(type_),
+        layout(Invalid) {}
+
+    Property(property_t id_, std::string name_, void *ptr_, PropertyType type_, layout_t layout_, size_t sx_ size_t sy) :
+        id(id_),
+        name(name_),
+        type(type_),
+        layout(layout_),
+        sx(sx_),
+        sy(sy_) {
+
+        PAIRS_ASSERT(type != Prop_Invalid && layout_ != Invalid && sx_ > 0 && sy_ > 0);
+        ptr = static_cast<void *>(new VectorPtr(ptr_, layout_, sx_, sy_));
+    }
 
 public:
     property_t getId() { return id; }
     std::string getName() { return name; }
-    void *getPointer() { return ptr; }
-    void setPointer(void *ptr_) { ptr = ptr_; }
+    void *getPointer() { return ptr.get(); }
+    void setPointer(void *ptr_) { ptr = (layout == Invalid) ? ptr_ : static_cast<void *>(new VectorPtr(ptr_, layout, sx, sy)); }
 
 private:
     property_t id;
-    layout_t layout;
     std::string name;
-    void *ptr;
+    std::shared_ptr<void *> ptr;
+    PropertyType type;
+    layout_t layout;
+    size_t sx, sy;
+};
+
+class VectorPtr {
+    VectorPtr(void *ptr_, layout_t layout_, size_t sx_, size_t sy_) : ptr(ptr_), layout(layout_), sx(sx_), sy(sy_) {}
+
+public:
+    double &operator[](int i, int j) {
+        double *dptr = static_cast<double *>(ptr.get());
+        if(layout == AoS) { return dptr[i * sy + j]; }
+        if(layout == SoA) { return dptr[j * sx + i]; }
+        PAIRS_ERROR("VectorPtr::operator[]: Invalid data layout!");
+        return 0.0;
+    }
+
+private:
+    std::shared_ptr<void *> ptr;
+    layout_t layout;
+    size_t sx, sy;
 };
 
 template<typename real_t>
diff --git a/runtime/read_from_file.hpp b/runtime/read_from_file.hpp
index b1bb7e7745d03e9ea297195bdecc2619839a4b94..d523479020495c894bd7b5825b921d1c2314ae9c 100644
--- a/runtime/read_from_file.hpp
+++ b/runtime/read_from_file.hpp
@@ -4,46 +4,52 @@
 
 namespace pairs {
 
-size_t read_particle_data(PairsSim *ps, const char *filename, property_t properties[], size_t nprops) {
+size_t read_particle_data(PairsSim *ps, const char *filename, double *grid_buffer, property_t properties[], size_t nprops) {
     std::ifstream in_file(filename);
     std::string line;
     size_t n = 0;
 
     if(in_file.is_open()) {
         while(getline(in_file, line)) {
-            int p = 0;
+            int i = 0;
             char *in0 = strtok(line, ",");
             while(in0 != NULL) {
-                PAIRS_ASSERT(p < nprops);
-                prop_type = ps->getPropertyType(p);
-
-                switch(prop_type) {
-                    case PROPERTY_TYPE_VECTOR:
-                        auto vector_ptr = ps->getVectorPropertyMutablePtr(p);
-                        char *in1 = strtok(NULL, ",");
-                        PAIRS_ASSERT(in1 != NULL);
-                        char *in2 = strtok(NULL, ",");
-                        PAIRS_ASSERT(in2 != NULL);
-                        vector_ptr[n] = Vector3(std::stod(in0), std::stod(in1), std::stod(in2));
-                        break;
-
-                    case PROPERTY_TYPE_INTEGER:
-                        auto int_ptr = ps->getIntegerPropertyMutablePtr(p);
-                        int_ptr[n] = std::stoi(in0);
-                        break;
-
-                    case PROPERTY_TYPE_FLOAT:
-                        auto float_ptr = ps->getFloatPropertyMutablePtr(p);
-                        float_ptr[n] = std::stod(in0);
-                        break;
-
-                    default:
-                        fprintf(stderr, "read_particle_data(): Invalid property type!");
-                        return -1;
+                if(grid_data_unread) {
+                    PAIRS_ASSERT(i < ps->getNumDims() * 2);
+                    grid_buffer[i] = std::stod(in0);
+                } else {
+                    PAIRS_ASSERT(i < nprops);
+                    property_t p = properties[i];
+                    auto prop_type = ps->getProperty(p)->getType();
+
+                    switch(prop_type) {
+                        case Prop_Vector:
+                            auto vector_ptr = ps->getVectorPropertyPtr(p);
+                            char *in1 = strtok(NULL, ",");
+                            PAIRS_ASSERT(in1 != NULL);
+                            char *in2 = strtok(NULL, ",");
+                            PAIRS_ASSERT(in2 != NULL);
+                            vector_ptr[n] = Vector3(std::stod(in0), std::stod(in1), std::stod(in2));
+                            break;
+
+                        case Prop_Integer:
+                            auto int_ptr = ps->getIntegerPropertyPtr(p);
+                            int_ptr[n] = std::stoi(in0);
+                            break;
+
+                        case Prop_Float:
+                            auto float_ptr = ps->getFloatPropertyPtr(p);
+                            float_ptr[n] = std::stod(in0);
+                            break;
+
+                        default:
+                            fprintf(stderr, "read_particle_data(): Invalid property type!");
+                            return -1;
+                    }
                 }
 
                 in0 = strtok(NULL, ",");
-                p++;
+                i++;
             }
 
             n++;
diff --git a/sim/cell_lists.py b/sim/cell_lists.py
index 4c1400a70e38d128e547e36ecedfca079c4fd24f..76d7224a726cd55e576618a536ed5d485c9a6b98 100644
--- a/sim/cell_lists.py
+++ b/sim/cell_lists.py
@@ -17,20 +17,17 @@ class CellLists:
         self.cutoff_radius = cutoff_radius
 
         self.nneighbor_cells = [
-            math.ceil(cutoff_radius / (
-                spacing if not isinstance(spacing, list)
-                else spacing[d]))
-            for d in range(0, sim.dimensions)]
+            math.ceil(cutoff_radius / (spacing if not isinstance(spacing, list) else spacing[d])) for d in range(sim.ndims())
+        ]
 
         self.nstencil = self.sim.add_var('nstencil', Type_Int)
-        self.nstencil_max = reduce((lambda x, y: x * y), [
-            self.nneighbor_cells[d] * 2 + 1 for d in range(0, sim.dimensions)])
-
-        self.ncells_all = self.sim.add_var('ncells_all', Type_Int)
+        self.nstencil_max = reduce((lambda x, y: x * y), [self.nneighbor_cells[d] * 2 + 1 for d in range(sim.ndims())])
+        self.ncells = self.sim.add_var('ncells', Type_Int, 1)
+        self.ncells_capacity = self.sim.add_var('ncells_capacity', Type_Int, 100)
+        self.dim_ncells = self.sim.add_static_array('dim_cells', self.sim.ndims(), Type_Int)
         self.cell_capacity = self.sim.add_var('cell_capacity', Type_Int, 20)
-        self.ncells = self.sim.add_static_array('ncells', self.sim.dimensions, Type_Int)
-        self.cell_particles = self.sim.add_array('cell_particles', [self.ncells_all, self.cell_capacity], Type_Int)
-        self.cell_sizes = self.sim.add_array('cell_sizes', self.ncells_all, Type_Int)
+        self.cell_particles = self.sim.add_array('cell_particles', [self.ncells_capacity, self.cell_capacity], Type_Int)
+        self.cell_sizes = self.sim.add_array('cell_sizes', self.ncells_capacity, Type_Int)
         self.stencil = self.sim.add_array('stencil', self.nstencil_max, Type_Int)
         self.particle_cell = self.sim.add_array('particle_cell', self.sim.particle_capacity, Type_Int)
 
@@ -41,28 +38,29 @@ class CellListsStencilBuild:
 
     def lower(self):
         cl = self.cell_lists
-        ncells = cl.sim.grid.get_ncells_for_spacing(cl.spacing)
+        grid = cl.sim.grid
         index = None
+        nall = 1
 
         cl.sim.clear_block()
         cl.sim.add_statement(Print(cl.sim, "CellListsStencilBuild"))
 
-        nall = 1
-        for d in range(0, cl.sim.dimensions):
-            cl.ncells[d].set(ncells[d])
-            nall *= ncells[d]
+        for d in range(cl.sim.ndims()):
+            cl.dim_ncells[d].set(Cast.int(cl.sim, (grid.max(d) - grid.min(d)) / cl.spacing))
+            nall *= cl.dim_ncells[d]
 
-        cl.ncells_all.set_initial_value(nall)
+        cl.ncells.set(nall)
+        for resize in Resize(cl.sim, cl.ncells_capacity):
+            for _ in Filter(cl.sim, cl.ncells >= cl.ncells_capacity):
+                resize.set(cl.ncells)
 
         for _ in cl.sim.nest_mode():
             cl.nstencil.set(0)
-            for d in range(0, cl.sim.dimensions):
+            for d in range(cl.sim.ndims()):
                 nneigh = cl.nneighbor_cells[d]
                 for d_idx in For(cl.sim, -nneigh, nneigh + 1):
-                    index = (d_idx if index is None
-                             else index * cl.ncells[d - 1] + d_idx)
-
-                    if d == cl.sim.dimensions - 1:
+                    index = (d_idx if index is None else index * cl.dim_ncells[d - 1] + d_idx)
+                    if d == cl.sim.ndims() - 1:
                         cl.stencil[cl.nstencil].set(index)
                         cl.nstencil.set(cl.nstencil + 1)
 
@@ -76,27 +74,26 @@ class CellListsBuild:
     def lower(self):
         cl = self.cell_lists
         grid = cl.sim.grid
-        spc = cl.spacing
         positions = cl.sim.property('position')
 
         cl.sim.clear_block()
         cl.sim.add_statement(Print(cl.sim, "CellListsBuild"))
         for resize in Resize(cl.sim, cl.cell_capacity):
-            for c in For(cl.sim, 0, cl.ncells_all):
+            for c in For(cl.sim, 0, cl.ncells):
                 cl.cell_sizes[c].set(0)
 
             for i in ParticleFor(cl.sim, local_only=False):
                 cell_index = [
-                    Cast.int(cl.sim, (positions[i][d] - grid.min(d)) / spc)
-                    for d in range(0, cl.sim.dimensions)]
+                    Cast.int(cl.sim, (positions[i][d] - grid.min(d)) / cl.spacing)
+                    for d in range(0, cl.sim.ndims())]
 
                 flat_idx = None
-                for d in range(0, cl.sim.dimensions):
+                for d in range(0, cl.sim.ndims()):
                     flat_idx = (cell_index[d] if flat_idx is None
-                                else flat_idx * cl.ncells[d] + cell_index[d])
+                                else flat_idx * cl.dim_ncells[d] + cell_index[d])
 
                 cell_size = cl.cell_sizes[flat_idx]
-                for _ in Filter(cl.sim, BinOp.and_op(flat_idx >= 0, flat_idx <= cl.ncells_all)):
+                for _ in Filter(cl.sim, BinOp.and_op(flat_idx >= 0, flat_idx <= cl.ncells)):
                     for cond in Branch(cl.sim, cell_size >= cl.cell_capacity):
                         if cond:
                             resize.set(cell_size)
diff --git a/sim/grid.py b/sim/grid.py
index a563a5093b6ad10dba08fd9f0b19416f0712c57b..0eea7ed01ac93abaa86715b0fe107ed980479193 100644
--- a/sim/grid.py
+++ b/sim/grid.py
@@ -1,3 +1,6 @@
+from ir.data_types import Type_Float
+
+
 class Grid:
     def __init__(self, sim, config):
         self.sim = sim
@@ -34,10 +37,6 @@ class Grid:
     def length(self, dim):
         return self.max(dim) - self.min(dim)
 
-    def get_ncells_for_spacing(self, spacing):
-        return [int(self.max(d) - self.min(d) / spacing)
-                for d in range(0, self.ndims)]
-
 
 class Grid2D(Grid):
     def __init__(self, sim, xmin, xmax, ymin, ymax):
@@ -49,3 +48,14 @@ class Grid3D(Grid):
     def __init__(self, sim, xmin, xmax, ymin, ymax, zmin, zmax):
         config = [[xmin, xmax], [ymin, ymax], [zmin, zmax]]
         super().__init__(sim, config)
+
+
+class MutableGrid(Grid):
+    last_id = 0
+
+    def __init__(self, sim, ndims):
+        self.id = MutableGrid.last_id
+        prefix = f"grid{self.id}_"
+        config = [[sim.add_var(f"{prefix}d{d}_min", Type_Float), sim.add_var(f"{prefix}d{d}_max", Type_Float)] for d in range(ndims)]
+        super().__init__(sim, config)
+        MutableGrid.last_id += 1
diff --git a/sim/lattice.py b/sim/lattice.py
index 70671416e3ade33250cbd65dc878db1e32efbcf1..b9985e670a729e35a235afe2cf7ddec362c3b47f 100644
--- a/sim/lattice.py
+++ b/sim/lattice.py
@@ -16,7 +16,7 @@ class ParticleLattice():
 
         self.sim.clear_block()
         for _ in self.sim.nest_mode():
-            for d in range(0, self.sim.dimensions):
+            for d in range(0, self.sim.ndims()):
                 d_min, d_max = self.grid.range(d)
                 n = int((d_max - d_min) / self.spacing[d] - 0.001) + 1
 
@@ -24,17 +24,17 @@ class ParticleLattice():
                     # index = (d_idx if index is None else index * n + d_idx)
                     loop_indexes.append(d_idx)
 
-                    if d == self.sim.dimensions - 1:
+                    if d == self.sim.ndims() - 1:
                         index = self.sim.nlocal
 
-                        for d_ in range(0, self.sim.dimensions):
+                        for d_ in range(0, self.sim.ndims()):
                             pos = self.grid.min(d_) + self.spacing[d_] * loop_indexes[d_]
                             self.positions[index][d_].set(pos)
 
                         for prop in [p for p in self.sim.properties.all()
                                      if p.volatile is False and p.name() != self.positions.name()]:
                             if prop.type() == Type_Vector:
-                                for d_ in range(0, self.sim.dimensions):
+                                for d_ in range(0, self.sim.ndims()):
                                     prop[index][d_].set(prop.default()[d_])
 
                             else:
diff --git a/sim/particle_simulation.py b/sim/particle_simulation.py
index 5b82cd4ce611a1e5bfeeb2f124e52f8970fa7c4c..870985fc7f7348366a1efc5e028064b08dd2c712 100644
--- a/sim/particle_simulation.py
+++ b/sim/particle_simulation.py
@@ -26,6 +26,7 @@ from transformations.set_used_bin_ops import set_used_bin_ops
 from transformations.simplify import simplify_expressions
 from transformations.LICM import move_loop_invariant_code
 
+
 class ParticleSimulation:
     def __init__(self, code_gen, dims=3, timesteps=100):
         self.code_gen = code_gen
@@ -48,7 +49,7 @@ class ParticleSimulation:
         self.block = Block(self, [])
         self.setups = SetupWrapper()
         self.kernels = KernelWrapper()
-        self.dimensions = dims
+        self.dims = dims
         self.ntimesteps = timesteps
         self.expr_id = 0
         self.iter_id = 0
@@ -56,6 +57,9 @@ class ParticleSimulation:
         self.nparticles = self.nlocal + self.nghost
         self.properties.add_capacity(self.particle_capacity)
 
+    def ndims(self):
+        return self.dims
+
     def add_real_property(self, prop_name, value=0.0, vol=False):
         assert self.property(prop_name) is None, f"Property already defined: {prop_name}"
         return self.properties.add(prop_name, Type_Float, value, vol)
@@ -191,8 +195,8 @@ class ParticleSimulation:
         timestep.add(Block(self, VTKWrite(self, self.vtk_file, timestep.timestep() + 1)))
 
         body = Block.from_list(self, [
-            CellListsStencilBuild(self.cell_lists).lower(),
             self.setups.lower(),
+            CellListsStencilBuild(self.cell_lists).lower(),
             Block(self, VTKWrite(self, self.vtk_file, 0)),
             timestep.as_block()
         ])
diff --git a/sim/pbc.py b/sim/pbc.py
index e99694134b36d1819406cfcf95fe5a4908502f49..1f066be079ec89809b1ec3c8f4eb370823f6506c 100644
--- a/sim/pbc.py
+++ b/sim/pbc.py
@@ -15,7 +15,7 @@ class PBC:
         self.npbc = sim.add_var('npbc', Type_Int)
         self.pbc_capacity = sim.add_var('pbc_capacity', Type_Int, 100)
         self.pbc_map = sim.add_array('pbc_map', [self.pbc_capacity], Type_Int)
-        self.pbc_mult = sim.add_array('pbc_mult', [self.pbc_capacity, sim.dimensions], Type_Int)
+        self.pbc_mult = sim.add_array('pbc_mult', [self.pbc_capacity, sim.ndims()], Type_Int)
 
 
 class UpdatePBC:
@@ -24,7 +24,7 @@ class UpdatePBC:
 
     def lower(self):
         sim = self.pbc.sim
-        ndims = sim.dimensions
+        ndims = sim.ndims()
         grid = self.pbc.grid
         npbc = self.pbc.npbc
         pbc_map = self.pbc.pbc_map
@@ -49,7 +49,7 @@ class EnforcePBC:
 
     def lower(self):
         sim = self.pbc.sim
-        ndims = sim.dimensions
+        ndims = sim.ndims()
         grid = self.pbc.grid
         positions = sim.property('position')
 
@@ -73,7 +73,7 @@ class SetupPBC:
 
     def lower(self):
         sim = self.pbc.sim
-        ndims = sim.dimensions
+        ndims = sim.ndims()
         grid = self.pbc.grid
         cutneigh = self.pbc.cutneigh
         npbc = self.pbc.npbc
diff --git a/sim/properties.py b/sim/properties.py
index b9a87e91d6c3ba241ed1bba63c7631c33173d9b8..f99d85778cad1d499e02600e880c843f299977f1 100644
--- a/sim/properties.py
+++ b/sim/properties.py
@@ -18,7 +18,7 @@ class PropertiesAlloc:
             if p.type() == Type_Float:
                 sizes = [capacity]
             elif p.type() == Type_Vector:
-                sizes = [capacity * self.sim.dimensions]
+                sizes = [capacity * self.sim.ndims()]
             else:
                 raise Exception("Invalid property type!")
 
diff --git a/sim/read_from_file.py b/sim/read_from_file.py
index 546966cd9fc6d532d641c10be7a368e866ec8b69..6a597a2a50175f3bb2bd0468843e9bafe5bdb345 100644
--- a/sim/read_from_file.py
+++ b/sim/read_from_file.py
@@ -1,49 +1,24 @@
-from ir.data_types import Type_Int, Type_Float, Type_Vector
-from ir.loops import For
-from sim.grid import Grid
+from ir.data_types import Type_Float
+from ir.functions import Call_Int
+from ir.properties import PropertyList
+from sim.grid import MutableGrid
+
 
 class ReadFromFile():
     def __init__(self, sim, filename, props):
         self.sim = sim
         self.filename = filename
-        self.props = props
-        self.grid = None
+        self.props = PropertyList(sim, props)
+        self.grid = MutableGrid(sim, sim.ndims())
+        self.grid_buffer = self.sim.add_static_array("grid_buffer", [self.sim.ndims() * 2], Type_Float)
 
     def lower(self):
-        ndims = None
-        nlocal = self.sim.nlocal
-        line_number = 0
-
         self.sim.clear_block()
-        with open(self.filename, "r") as fp:
-            for line in fp:
-                current_data = line.split(',')
-                if line_number == 0:
-                    assert len(current_data) % 2 == 0, "Number of dimensions in header must be even!"
-                    ndims = int(len(current_data) / 2)
-                    config = [[float(current_data[d * 2]), float(current_data[d * 2 + 1])] for d in range(0, ndims)]
-                    self.grid = Grid(self.sim, config)
-                else:
-                    i = 0
-                    for p in self.props:
-                        if p.type() == Type_Vector:
-                            for d in range(0, ndims):
-                                p[nlocal][d].set(float(current_data[i + d]))
-
-                            i += ndims
-
-                        else:
-                            if p.type() == Type_Int:
-                                p[nlocal].set(int(current_data[i]))
-                            elif p.type() == Type_Float:
-                                p[nlocal].set(float(current_data[i]))
-                            else:
-                                raise Exception(f"Invalid property type at line {line_number}!")
-
-                            i += 1
-
-                    nlocal.set(nlocal + 1)
+        self.sim.nlocal.set(Call_Int(self.sim, "pairs::read_particle_data",
+                            [self.filename, self.grid_buffer, self.props, self.props.length()]))
 
-                line_number += 1
+        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])
 
         return self.sim.block
diff --git a/sim/resize.py b/sim/resize.py
index afeb9cd99785ce51bc5c5860a5a172a70f7e3053..4c354fa86399b3f7902d6a47e5cae296e8985e44 100644
--- a/sim/resize.py
+++ b/sim/resize.py
@@ -27,7 +27,7 @@ class Resize:
                     capacity = sum(self.sim.properties.capacities)
                     for p in properties.all():
                         if p.type() == Type_Vector:
-                            sizes = capacity * self.sim.dimensions
+                            sizes = capacity * self.sim.ndims()
                         else:
                             sizes = capacity
 
diff --git a/transformations/flatten.py b/transformations/flatten.py
index d673c190419d8d1fbab362ac576b80d7d9f0f3d9..7b0ac8e739d394ee0d782bef84ae7db3ffecc10a 100644
--- a/transformations/flatten.py
+++ b/transformations/flatten.py
@@ -17,7 +17,7 @@ class FlattenPropertyAccesses(Mutator):
                 flat_index = None
 
                 if layout == Layout_AoS:
-                    flat_index = ast_node.rhs * ast_node.sim.dimensions + i
+                    flat_index = ast_node.rhs * ast_node.sim.ndims() + i
 
                 elif layout == Layout_SoA:
                     flat_index = i * ast_node.sim.particle_capacity + ast_node.rhs