diff --git a/ast/arrays.py b/ast/arrays.py
index 5b13201105cbc977f8494cbec74f80c7e8391982..f1ffb86d8160d4df7d67b8d8e8df703082c88f35 100644
--- a/ast/arrays.py
+++ b/ast/arrays.py
@@ -27,9 +27,9 @@ class Arrays:
         return self.arrays
 
     def find(self, a_name):
-        arr = [a for a in self.arrays if a.name() == a_name]
-        if arr:
-            return arr[0]
+        array = [a for a in self.arrays if a.name() == a_name]
+        if array:
+            return array[0]
 
         return None
 
@@ -138,8 +138,7 @@ class ArrayAccess:
         return Expr(self.sim, other, self, '*')
 
     def __getitem__(self, index):
-        assert self.index is None, \
-            "Number of indexes higher than array dimension!"
+        assert self.index is None, "Number of indexes higher than array dimension!"
         self.indexes.append(as_lit_ast(self.sim, index))
         self.check_and_set_index()
         return self
@@ -171,7 +170,8 @@ class ArrayAccess:
         return self.sim.add_statement(Assign(self.sim, self, self + other))
 
     def type(self):
-        return self.array.type() if self.index is None else Type_Array
+        return self.array.type()
+        # return self.array.type() if self.index is None else Type_Array
 
     def scope(self):
         if self.index is None:
@@ -192,12 +192,10 @@ class ArrayAccess:
         agen = self.array.generate()
         igen = self.index.generate()
         if mem is False and self.generated is False:
-            self.sim.code_gen.generate_array_access(
-                self.acc_id, self.array.type(), agen, igen)
+            self.sim.code_gen.generate_array_access(self.acc_id, self.array.type(), agen, igen)
             self.generated = True
 
-        return self.sim.code_gen.generate_array_access_ref(
-            self.acc_id, agen, igen, mem)
+        return self.sim.code_gen.generate_array_access_ref(self.acc_id, agen, igen, mem)
 
     def transform(self, fn):
         self.array = self.array.transform(fn)
diff --git a/ast/cast.py b/ast/cast.py
index d069486df37437f6289e00a9251db5578408a9ed..970109e06a25556468d34530f076d30961a0cc3f 100644
--- a/ast/cast.py
+++ b/ast/cast.py
@@ -26,8 +26,7 @@ class Cast:
         return [self.expr]
 
     def generate(self, mem=False):
-        return self.sim.code_gen.generate_cast(
-            self.cast_type, self.expr.generate())
+        return self.sim.code_gen.generate_cast(self.cast_type, self.expr.generate())
 
     def transform(self, fn):
         self.expr = self.expr.transform(fn)
diff --git a/ast/expr.py b/ast/expr.py
index 728974f3e91715274d20c8aae93b8b87a0112ac8..ac889de3b5c27616b250834f34a3b8557495f326 100644
--- a/ast/expr.py
+++ b/ast/expr.py
@@ -86,8 +86,7 @@ class Expr:
         return self.__getitem__(2)
 
     def __getitem__(self, index):
-        assert self.lhs.type() == Type_Vector, \
-            "Cannot use operator [] on specified type!"
+        assert self.lhs.type() == Type_Vector, "Cannot use operator [] on specified type!"
         return ExprVec(self.sim, self, as_lit_ast(self.sim, index))
 
     def generated_vector_index(self, index):
@@ -126,6 +125,7 @@ class Expr:
         if lhs_type == Type_Float or rhs_type == Type_Float:
             return Type_Float
 
+        print(f"{lhs} ({lhs_type}) -- {rhs} ({rhs_type})\n")
         return None
 
     def type(self):
@@ -146,15 +146,13 @@ class Expr:
         lhs_expr = self.lhs.generate(mem)
         rhs_expr = self.rhs.generate()
         if self.op == '[]':
-            return self.sim.code_gen.generate_expr_access(
-                lhs_expr, rhs_expr, self.mem)
+            return self.sim.code_gen.generate_expr_access(lhs_expr, rhs_expr, self.mem)
 
         if self.generated is False:
             assert self.expr_type != Type_Vector, \
                 "Vector code must be generated through ExprVec class!"
 
-            self.sim.code_gen.generate_expr(
-                self.expr_id, self.expr_type, lhs_expr, rhs_expr, self.op)
+            self.sim.code_gen.generate_expr(self.expr_id, self.expr_type, lhs_expr, rhs_expr, self.op)
             self.generated = True
 
         return self.sim.code_gen.generate_expr_ref(self.expr_id)
@@ -162,19 +160,19 @@ class Expr:
     def generate_inline(self, mem=False, recursive=False):
         inline_lhs_expr = recursive and isinstance(self.lhs, Expr)
         inline_rhs_expr = recursive and isinstance(self.rhs, Expr)
-        lhs_expr = (self.lhs.generate_inline(recursive, mem) if inline_lhs_expr
-                    else self.lhs.generate(mem))
-        rhs_expr = (self.rhs.generate_inline(recursive) if inline_rhs_expr
-                    else self.rhs.generate())
+        lhs_expr = \
+            self.lhs.generate_inline(recursive, mem) if inline_lhs_expr \
+            else self.lhs.generate(mem)
+        rhs_expr = \
+            self.rhs.generate_inline(recursive) if inline_rhs_expr \
+            else self.rhs.generate()
 
         if self.op == '[]':
-            return self.sim.code_gen.generate_expr_access(
-                lhs_expr, rhs_expr, self.mem)
+            return self.sim.code_gen.generate_expr_access(lhs_expr, rhs_expr, self.mem)
 
         assert self.expr_type != Type_Vector, \
             "Vector code must be generated through ExprVec class!"
-        return self.sim.code_gen.generate_inline_expr(
-                lhs_expr, rhs_expr, self.op)
+        return self.sim.code_gen.generate_inline_expr(lhs_expr, rhs_expr, self.op)
 
     def transform(self, fn):
         self.lhs = self.lhs.transform(fn)
@@ -250,8 +248,7 @@ class ExprVec():
 
         index_expr = self.index.generate()
         if self.expr.op == '[]':
-            return self.sim.code_gen.generate_expr_access(
-                self.expr.generate(), index_expr, True)
+            return self.sim.code_gen.generate_expr_access(self.expr.generate(), index_expr, True)
 
         if self.expr.generated_vector_index(index_expr):
             self.sim.code_gen.generate_vec_expr(
@@ -264,8 +261,7 @@ class ExprVec():
 
             self.expr.vec_generated.append(index_expr)
 
-        return self.sim.code_gen.generate_vec_expr_ref(
-            self.expr.expr_id, index_expr, self.expr.mem)
+        return self.sim.code_gen.generate_vec_expr_ref(self.expr.expr_id, index_expr, self.expr.mem)
 
     def transform(self, fn):
         self.lhs = self.lhs.transform(fn)
diff --git a/ast/loops.py b/ast/loops.py
index 9cfcebe76258f198ee68c488f5b26b829f5914aa..f0e27c30b64151c8e24c8135eb116169c0b38909 100644
--- a/ast/loops.py
+++ b/ast/loops.py
@@ -110,8 +110,7 @@ class ParticleFor(For):
         return f"ParticleFor<>"
 
     def generate(self):
-        self.sim.code_gen.generate_for_preamble(
-            self.iterator.generate(), 0, self.sim.nparticles.generate())
+        self.sim.code_gen.generate_for_preamble(self.iterator.generate(), 0, self.sim.nlocal.generate())
         self.block.generate()
         self.sim.code_gen.generate_for_epilogue()
 
diff --git a/ast/properties.py b/ast/properties.py
index 4b061d544d55ab100240681120d6d83f3cab9763..02139fcd1afeb5e3dfe393689a5a226b4b9490c4 100644
--- a/ast/properties.py
+++ b/ast/properties.py
@@ -53,6 +53,9 @@ class Property:
     def layout(self):
         return self.prop_layout
 
+    def default(self):
+        return self.default_value
+
     def scope(self):
         return self.sim.global_scope
 
diff --git a/ast/transform.py b/ast/transform.py
index 1539ffc061759e207bd667b92ef7eb8692d46955..e2fa60e4368bf90f17715231ac69a754f9f079c6 100644
--- a/ast/transform.py
+++ b/ast/transform.py
@@ -35,7 +35,7 @@ class Transform:
 
                 elif layout == Layout_SoA:
                     flat_index = \
-                        ast.index * ast.expr.sim.nparticles + ast.expr.rhs
+                        ast.index * ast.expr.sim.particle_capacity + ast.expr.rhs
 
                 else:
                     raise Exception("Invalid property layout!")
diff --git a/code_gen/cgen.py b/code_gen/cgen.py
index 00b04075b0429716d270befed9b38ed5cea14577..e2642b0da2107409e290beeaf08e7e61aea951ef 100644
--- a/code_gen/cgen.py
+++ b/code_gen/cgen.py
@@ -74,8 +74,7 @@ class CGen:
         return f"sizeof({tkw})"
 
     def generate_for_preamble(iter_id, rmin, rmax):
-        printer.print(
-            f"for(int {iter_id} = {rmin}; {iter_id} < {rmax}; {iter_id}++) {{")
+        printer.print(f"for(int {iter_id} = {rmin}; {iter_id} < {rmax}; {iter_id}++) {{")
 
     def generate_for_epilogue():
         printer.print("}")
@@ -98,7 +97,7 @@ class CGen:
         return f"{lhs}[{rhs}]" if mem else f"{lhs}_{rhs}"
 
     def generate_vec_expr_ref(expr_id, index, mem):
-        return (f"e{expr_id}[{index}]" if mem else f"e{expr_id}_{index}")
+        return f"e{expr_id}[{index}]" if mem else f"e{expr_id}_{index}"
 
     def generate_vec_expr(expr_id, index, lhs, rhs, op, mem):
         ref = CGen.generate_vec_expr_ref(expr_id, index, mem)
diff --git a/particle.py b/particle.py
index 400529555ed81c2e27b2c4d9b0eb9b23d93c2cbc..cb4aea6977cb8bbd4d7e520d215f9716edf1a80c 100644
--- a/particle.py
+++ b/particle.py
@@ -10,7 +10,7 @@ sigma6 = sigma ** 6
 
 psim = pt.simulation()
 mass = psim.add_real_property('mass', 1.0)
-position = psim.add_vector_property('position', layout=Layout_SoA)
+position = psim.add_vector_property('position')
 velocity = psim.add_vector_property('velocity')
 force = psim.add_vector_property('force', vol=True)
 
diff --git a/sim/cell_lists.py b/sim/cell_lists.py
index 5f6bde9ffa909182b7aa33e758fac731f48a8598..7bb7efc058992263f8560b61b34b74c30cc492e4 100644
--- a/sim/cell_lists.py
+++ b/sim/cell_lists.py
@@ -30,7 +30,7 @@ class CellLists:
         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.stencil = self.sim.add_array('stencil', self.nstencil_max, Type_Int)
-        self.particle_cell = self.sim.add_array('particle_cell', self.sim.nparticles, Type_Int)
+        self.particle_cell = self.sim.add_array('particle_cell', self.sim.particle_capacity, Type_Int)
 
 
 class CellListsStencilBuild:
diff --git a/sim/lattice.py b/sim/lattice.py
index b4f6aea0715d7dd4b5a0ef64caf8ce61b6a364e3..4ea51d1183377af42d2642bfdab1d7ad086b5d58 100644
--- a/sim/lattice.py
+++ b/sim/lattice.py
@@ -1,3 +1,4 @@
+from ast.data_types import Type_Vector
 from ast.loops import For
 
 
@@ -20,15 +21,25 @@ class ParticleLattice():
                 n = int((d_max - d_min) / self.spacing[d] - 0.001) + 1
 
                 for d_idx in For(self.sim, 0, n):
-                    index = (d_idx if index is None else index * n + d_idx)
+                    # index = (d_idx if index is None else index * n + d_idx)
                     loop_indexes.append(d_idx)
 
                     if d == self.sim.dimensions - 1:
+                        index = self.sim.nlocal
+
                         for d_ in range(0, self.sim.dimensions):
-                            pos = self.grid.min(d_) + \
-                                  self.spacing[d_] * loop_indexes[d_]
+                            pos = self.grid.min(d_) + self.spacing[d_] * loop_indexes[d_]
                             self.positions[index][d_].set(pos)
 
-                        self.sim.nparticles.set(self.sim.nparticles + 1)
+                        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):
+                                    prop[index][d_].set(prop.default()[d_])
+
+                            else:
+                                prop[index].set(prop.default())
+
+                        self.sim.nlocal.set(self.sim.nlocal + 1)
 
         return self.sim.block
diff --git a/sim/particle_simulation.py b/sim/particle_simulation.py
index 9d093fec35355eca42adc100df7ff140853d8600..0b7ec054c089f1555bbfda19f33b56746f7ff4ea 100644
--- a/sim/particle_simulation.py
+++ b/sim/particle_simulation.py
@@ -26,7 +26,10 @@ class ParticleSimulation:
         self.properties = Properties(self)
         self.vars = Variables(self)
         self.arrays = Arrays(self)
-        self.nparticles = self.add_var('nparticles', Type_Int)
+        self.particle_capacity = self.add_var('particle_capacity', Type_Int, 100)
+        self.nlocal = self.add_var('nlocal', Type_Int)
+        self.nghost = self.add_var('nghost', Type_Int)
+        self.nparticles = self.nlocal + self.nghost
         self.grid = None
         self.cell_lists = None
         self.scope = []
@@ -159,10 +162,10 @@ class ParticleSimulation:
             CellListsStencilBuild(self.cell_lists).lower(),
             self.setups.lower(),
             Timestep(self, self.ntimesteps, [
-                (EnforcePBC(self.pbc).lower(), 20),
-                (SetupPBC(self.pbc).lower(), 20),
+                #(EnforcePBC(self.pbc).lower(), 20),
+                #(SetupPBC(self.pbc).lower(), 20),
                 (CellListsBuild(self.cell_lists).lower(), 20),
-                UpdatePBC(self.pbc).lower(),
+                #UpdatePBC(self.pbc).lower(),
                 PropertiesResetVolatile(self).lower(),
                 self.kernels.lower()
             ]).as_block()
diff --git a/sim/pbc.py b/sim/pbc.py
index 235fce79468ecceb0c882d8c549eb0950d6ed42a..f8c28fcb8a5c4701e9236ccdd2caca6baf613508 100644
--- a/sim/pbc.py
+++ b/sim/pbc.py
@@ -4,6 +4,7 @@ from ast.loops import For, ParticleFor
 from ast.select import Select
 from sim.resize import Resize
 
+
 class PBC:
     def __init__(self, sim, grid, cutneigh, pbc_flags=[1, 1, 1]):
         self.sim = sim
@@ -28,7 +29,7 @@ class UpdatePBC:
         pbc_map = self.pbc.pbc_map
         pbc_mult = self.pbc.pbc_mult
         positions = self.pbc.sim.property('position')
-        nlocal = self.pbc.sim.nparticles
+        nlocal = self.pbc.sim.nlocal
 
         sim.clear_block()
         for i in For(sim, 0, npbc):
@@ -77,7 +78,7 @@ class SetupPBC:
         pbc_map = self.pbc.pbc_map
         pbc_mult = self.pbc.pbc_mult
         positions = self.pbc.sim.property('position')
-        nlocal = self.pbc.sim.nparticles
+        nlocal = self.pbc.sim.nlocal
 
         sim.clear_block()
         for capacity, resize in Resize(sim, pbc_capacity, [pbc_map, pbc_mult]):
diff --git a/sim/properties.py b/sim/properties.py
index 9668761a400249fcd7397e380c85986ff7bb10f8..acb34feae2255c68c31206feb188b5cbc0dc5926 100644
--- a/sim/properties.py
+++ b/sim/properties.py
@@ -8,18 +8,18 @@ class PropertiesDecl:
         self.sim = sim
 
     def lower(self):
-        nparticles = self.sim.nparticles
+        particle_capacity = self.sim.particle_capacity
 
         self.sim.clear_block()
         for p in self.sim.properties.all():
             sizes = []
             if p.type() == Type_Float:
-                sizes = [nparticles]
+                sizes = [particle_capacity]
             elif p.type() == Type_Vector:
                 if p.flattened:
-                    sizes = [nparticles * self.sim.dimensions]
+                    sizes = [particle_capacity * self.sim.dimensions]
                 else:
-                    sizes = [nparticles, self.sim.dimensions]
+                    sizes = [particle_capacity, self.sim.dimensions]
             else:
                 raise Exception("Invalid property type!")