diff --git a/examples/dem.py b/examples/dem.py
index 8ac9d8266a6a3625d025755b92e7d6f658609902..7d08a0f9cc62fd4bc74d58e8f31a1d28e6321211 100644
--- a/examples/dem.py
+++ b/examples/dem.py
@@ -4,27 +4,22 @@ import sys
 
 
 def linear_spring_dashpot(i, j):
-    penetration_depth = squared_distance(i, j) - radius[i] - radius[j]
-    skip_when(penetration_depth >= 0.0)
+    skip_when(penetration_depth(i, j) >= 0.0)
 
-    contact_normal = normalized(delta(i, j))
-    k = radius[j] + 0.5 * penetration_depth
-    contact_point = position[j] + contact_normal * k
+    velocity_wf_i = linear_velocity[i] + cross(angular_velocity[i], contact_point(i, j) - position[i])
+    velocity_wf_j = linear_velocity[j] + cross(angular_velocity[j], contact_point(i, j) - position[j])
 
-    velocity_wf_i = linear_velocity[i] + cross(angular_velocity[i], contact_point - position[i])
-    velocity_wf_j = linear_velocity[j] + cross(angular_velocity[j], contact_point - position[j])
     rel_vel = -velocity_wf_i - velocity_wf_j
-    rel_vel_n = dot(rel_vel, contact_normal) * contact_normal
+    rel_vel_n = dot(rel_vel, contact_normal(i, j)) * contact_normal(i, j)
     rel_vel_t = rel_vel - rel_vel_n
-
-    fN = stiffness_norm[i, j] * (-penetration_depth) * contact_normal + damping_norm[i, j] * rel_vel_n;
+    fN = stiffness_norm[i, j] * (-penetration_depth(i, j)) * contact_normal(i, j) + damping_norm[i, j] * rel_vel_n;
 
     tan_spring_disp = tangential_spring_displacement[i, j]
     impact_vel_magnitude = impact_velocity_magnitude[i, j]
     impact_magnitude = select(impact_vel_magnitude > 0.0, impact_vel_magnitude, length(rel_vel))
     sticking = is_sticking[i, j]
 
-    rotated_tan_disp = tan_spring_disp - contact_normal * (contact_normal * tan_spring_disp)
+    rotated_tan_disp = tan_spring_disp - contact_normal(i, j) * (contact_normal(i, j) * tan_spring_disp)
     new_tan_spring_disp = dt * rel_vel_t + \
                           select(squared_length(rotated_tan_disp) <= 0.0,
                                  zero_vector(),
@@ -44,8 +39,7 @@ def linear_spring_dashpot(i, j):
     n_sticking = select(cond1 or cond2 or fTLS_len < f_friction_abs_dynamic, 1, 0)
 
     if not cond1 and not cond2 and stiffness_tan[i, j] > 0.0:
-        tangential_spring_displacement[i, j] = \
-            (f_friction_abs * t - damping_tan[i, j] * rel_vel_t) / stiffness_tan[i, j]
+        tangential_spring_displacement[i, j] = (f_friction_abs * t - damping_tan[i, j] * rel_vel_t) / stiffness_tan[i, j]
 
     else:
         tangential_spring_displacement[i, j] = new_tan_spring_disp
@@ -94,8 +88,9 @@ psim.add_position('position')
 psim.add_property('mass', pairs.double(), 1.0)
 psim.add_property('linear_velocity', pairs.vector())
 psim.add_property('angular_velocity', pairs.vector())
-psim.add_property('force', pairs.vector(), vol=True)
+psim.add_property('force', pairs.vector(), volatile=True)
 psim.add_property('radius', pairs.double(), 1.0)
+psim.add_property('normal', pairs.vector())
 psim.add_feature('type', ntypes)
 psim.add_feature_property('type', 'stiffness_norm', pairs.double(), [stiffness_norm for i in range(ntypes * ntypes)])
 psim.add_feature_property('type', 'stiffness_tan', pairs.double(), [stiffness_tan for i in range(ntypes * ntypes)])
diff --git a/examples/lj.py b/examples/lj.py
index a9aa2ffe80844ad4b07233e776421211d58cd541..4fabc8a0b16d82597f8ed4cd206c2fbcd2749eb0 100644
--- a/examples/lj.py
+++ b/examples/lj.py
@@ -30,7 +30,7 @@ psim = pairs.simulation("lj", debug=True)
 psim.add_position('position')
 psim.add_property('mass', pairs.double(), 1.0)
 psim.add_property('linear_velocity', pairs.vector())
-psim.add_property('force', pairs.vector(), vol=True)
+psim.add_property('force', pairs.vector(), volatile=True)
 psim.add_feature('type', ntypes)
 psim.add_feature_property('type', 'epsilon', pairs.double(), [sigma for i in range(ntypes * ntypes)])
 psim.add_feature_property('type', 'sigma6', pairs.double(), [epsilon for i in range(ntypes * ntypes)])
diff --git a/src/pairs/ir/lit.py b/src/pairs/ir/lit.py
index 03484656fbe88d66f0eb49e5913eb4f29cbc8339..241000723f510d0ec1fde04ee509c1f79b8e085a 100644
--- a/src/pairs/ir/lit.py
+++ b/src/pairs/ir/lit.py
@@ -1,8 +1,8 @@
-from pairs.ir.ast_node import ASTNode
+from pairs.ir.ast_term import ASTTerm
 from pairs.ir.types import Types
 
 
-class Lit(ASTNode):
+class Lit(ASTTerm):
     def is_literal(a):
         return isinstance(a, (int, float, bool, str, list))
 
@@ -10,9 +10,6 @@ class Lit(ASTNode):
         return Lit(sim, a) if Lit.is_literal(a) else a
 
     def __init__(self, sim, value):
-        super().__init__(sim)
-        self.value = value
-
         type_mapping = {
             int: Types.Int32,
             float: Types.Double,
@@ -24,6 +21,11 @@ class Lit(ASTNode):
         self.lit_type = type_mapping.get(type(value), Types.Invalid)
         assert self.lit_type != Types.Invalid, "Invalid literal type!"
 
+        from pairs.ir.scalars import ScalarOp
+        from pairs.ir.vectors import VectorOp
+        super().__init__(sim, VectorOp if self.lit_type == Types.Vector else ScalarOp)
+        self.value = value
+
     def __str__(self):
         return f"Lit<{self.value}>"
 
diff --git a/src/pairs/ir/symbols.py b/src/pairs/ir/symbols.py
index afbaa93ce6a72db41089ddcd6c8da8f233c57522..6f2182d85c70374e1a84a3e1e55f1209c4012d05 100644
--- a/src/pairs/ir/symbols.py
+++ b/src/pairs/ir/symbols.py
@@ -22,4 +22,29 @@ class Symbol(ASTTerm):
         return self.sym_type
 
     def __getitem__(self, index):
-        return VectorAccess(self.sim, self, Lit.cvt(self.sim, index))
+        #return VectorAccess(self.sim, self, Lit.cvt(self.sim, index))
+        return SymbolAccess(self.sim, self, Lit.cvt(self.sim, index))
+
+
+class SymbolAccess(ASTTerm):
+    def __init__(self, sim, symbol, index):
+        assert symbol.type() == Types.Vector, "Only vector symbols can be indexed!"
+        super().__init__(sim, ScalarOp if symbol.type() != Types.Vector else VectorOp)
+        self._symbol = symbol
+        self._index = index
+        symbol.add_index_to_generate(index)
+
+    def __str__(self):
+        return f"SymbolAccess<{self._symbol}, {self._index}>"
+
+    def symbol(self):
+        return self._symbol
+
+    def index(self):
+        return self._index
+
+    def type(self):
+        if self._symbol.type() == Types.Vector:
+            return Types.Float
+
+        return self._symbol.type()
diff --git a/src/pairs/mapping/funcs.py b/src/pairs/mapping/funcs.py
index f18b24365ca25d8a3771e209c996e1b2fefc4251..333d90ca9fffcb2bc0e20896decb382b0927b7d7 100644
--- a/src/pairs/mapping/funcs.py
+++ b/src/pairs/mapping/funcs.py
@@ -78,9 +78,9 @@ class BuildParticleIR(ast.NodeVisitor):
 
         raise Exception("Invalid operator: {}".format(ast.dump(op)))
 
-    def __init__(self, sim, ctx_symbols={}, ctx_calls=[]):
+    def __init__(self, sim, ctx_symbols={}):
         self.sim = sim
-        self.ctx_symbols = ctx_symbols
+        self.ctx_symbols = ctx_symbols.copy()
         self.keywords = Keywords(sim)
 
     def add_symbols(self, symbols):
@@ -130,7 +130,7 @@ class BuildParticleIR(ast.NodeVisitor):
         if self.keywords.exists(func):
             return self.keywords(func, args)
 
-        if func == 'squared_distance' or func == 'delta':
+        if func in ['delta', 'squared_distance', 'penetration_depth', 'contact_point', 'contact_normal']:
             return self.ctx_symbols[f"__{func}__"]
 
         raise Exception(f"Undefined function called: {func}")
@@ -140,10 +140,7 @@ class BuildParticleIR(ast.NodeVisitor):
 
     def visit_Compare(self, node):
         #print(ast.dump(node))
-        valid_ops = (
-            ast.Eq, ast.NotEq, ast.Lt, ast.LtE, ast.Gt, ast.GtE, ast.Is, ast.IsNot, ast.In, ast.NotIn
-        )
-
+        valid_ops = (ast.Eq, ast.NotEq, ast.Lt, ast.LtE, ast.Gt, ast.GtE, ast.Is, ast.IsNot, ast.In, ast.NotIn)
         if len(node.ops) != 1 or not isinstance(node.ops[0], valid_ops):
             raise Exception(f"Chained comparisons or unsupported comparison found: {ast.dump(node)}")
 
@@ -231,29 +228,33 @@ def compute(sim, func, cutoff_radius=None, symbols={}):
     params = info.params()
     nparams = info.nparams()
 
+    # Compute functions must have parameters
+    assert nparams > 0, "Number of parameters from compute functions must be higher than zero!"
+
     # Convert literal symbols
     symbols = {symbol: Lit.cvt(sim, value) for symbol, value in symbols.items()}
 
-    # Start building IR
-    ir = BuildParticleIR(sim, symbols)
-    assert nparams > 0, "Number of parameters from compute functions must be higher than zero!"
-
     sim.init_block()
     sim.module_name(func.__name__)
 
     if nparams == 1:
         for i in ParticleFor(sim):
+            ir = BuildParticleIR(sim, symbols)
             ir.add_symbols({params[0]: i})
             ir.visit(tree)
 
     else:
-        pairs = ParticleInteraction(sim, nparams, cutoff_radius)
-        for i, j in pairs:
+        for interaction_data in ParticleInteraction(sim, nparams, cutoff_radius):
+            # Start building IR
+            ir = BuildParticleIR(sim, symbols)
             ir.add_symbols({
-                params[0]: i,
-                params[1]: j,
-                '__delta__': pairs.delta(),
-                '__squared_distance__': pairs.squared_distance()
+                params[0]: interaction_data.i(),
+                params[1]: interaction_data.j(),
+                '__delta__': interaction_data.delta(),
+                '__squared_distance__': interaction_data.squared_distance(),
+                '__penetration_depth__': interaction_data.penetration_depth(),
+                '__contact_point__': interaction_data.contact_point(),
+                '__contact_normal__': interaction_data.contact_normal()
             })
 
             ir.visit(tree)
diff --git a/src/pairs/sim/interaction.py b/src/pairs/sim/interaction.py
index 9b6eacb4790a1b005678355973bf521740cd8fde..6767960a7c7976a8ead4ba0503634b2e3170df3e 100644
--- a/src/pairs/sim/interaction.py
+++ b/src/pairs/sim/interaction.py
@@ -3,16 +3,20 @@ from pairs.ir.scalars import ScalarOp
 from pairs.ir.block import Block, pairs_inline
 from pairs.ir.branches import Filter
 from pairs.ir.loops import For, ParticleFor
+from pairs.ir.math import Sqrt
 from pairs.ir.types import Types
+from pairs.ir.vectors import Vector
 from pairs.sim.lowerable import Lowerable
+from pairs.sim.shapes import Shapes
 
 
 class Neighbor(ASTTerm):
-    def __init__(self, sim, neighbor_index, cell_id, particle_index):
+    def __init__(self, sim, neighbor_index, cell_id, particle_index, shape):
         super().__init__(sim, ScalarOp)
         self._neighbor_index = neighbor_index
         self._cell_id = cell_id
         self._particle_index = particle_index
+        self._shape = shape
 
     def __str__(self):
         return f"Neighbor<{self._neighbor_index}, {self._cell_id}>"
@@ -29,10 +33,13 @@ class Neighbor(ASTTerm):
     def particle_index(self):
         return self._particle_index
 
+    def shape(self):
+        return self._shape
+
 
 class NeighborFor:
-    def number_of_cases(neighbor_lists):
-        return 2 if neighbor_lists is None else 1
+    def number_of_cases(sim, use_cell_lists):
+        return 2 if sim.neighbor_lists is None or use_cell_lists else 1
 
     def __init__(self, sim, particle, cell_lists, neighbor_lists=None):
         self.sim = sim
@@ -45,80 +52,153 @@ class NeighborFor:
 
     def __iter__(self):
         if self.neighbor_lists is None:
-            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, ScalarOp.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, ScalarOp.neq(it, self.particle)):
-                            yield Neighbor(self.sim, nc, neigh_cell, it)
-
-            # Infinite particles
-            for inf_id in For(self.sim, 0, cl.cell_sizes[0]):
-                inf_particle = cl.cell_particles[0][inf_id]
-                for _ in Filter(self.sim, ScalarOp.neq(inf_particle, self.particle)):
-                    yield Neighbor(self.sim, inf_id, 0, inf_particle)
+            stencil = self.cell_lists.stencil
+            nstencil = self.cell_lists.nstencil
+            ncells = self.cell_lists.ncells
+            particle_cell = self.cell_lists.particle_cell
+            cell_particles = self.cell_lists.cell_particles
+            nshapes = self.cell_lists.nshapes
+
+            for shape in range(self.sim.max_shapes()):
+                for disp in For(self.sim, 0, self.cell_lists.nstencil):
+                    neigh_cell = particle_cell[self.particle] + stencil[disp]
+
+                    for _ in Filter(self.sim, ScalarOp.and_op(neigh_cell > 0, neigh_cell < ncells)):
+                        start = sum([nshapes[neigh_cell][s] for s in range(shape)], 0)
+
+                        for cell_particle in For(self.sim, start, start + nshapes[neigh_cell][shape]):
+                            particle_id = cell_particles[neigh_cell][cell_particle]
+                            for _ in Filter(self.sim, ScalarOp.neq(particle_id, self.particle)):
+                                yield Neighbor(self.sim, cell_particle, neigh_cell, particle_id, shape)
+
+                # Infinite particles
+                start = sum([nshapes[0][s] for s in range(shape)], 0)
+                for inf_id in For(self.sim, start, start + nshapes[0][shape]):
+                    inf_particle = cell_particles[0][inf_id]
+                    for _ in Filter(self.sim, ScalarOp.neq(inf_particle, self.particle)):
+                        yield Neighbor(self.sim, inf_id, 0, inf_particle, shape)
 
         else:
-            neighbor_lists = self.neighbor_lists
-            for k in For(self.sim, 0, neighbor_lists.numneighs[self.particle]):
-                yield Neighbor(self.sim, k, None, neighbor_lists.neighborlists[self.particle][k])
+            neighborlists = self.neighbor_lists.neighborlists
+            numneighs = self.neighbor_lists.numneighs
 
+            for shape in range(self.sim.max_shapes()):
+                for k in For(self.sim, 0, numneighs[self.particle][shape]):
+                    yield Neighbor(self.sim, k, None, neighborlists[self.particle][k], shape)
 
-class ParticleInteraction(Lowerable):
-    def __init__(self, sim, nbody, cutoff_radius, bypass_neighbor_lists=False):
-        super().__init__(sim)
-        self.nbody = nbody
-        self.cutoff_radius = cutoff_radius
-        self.bypass_neighbor_lists = bypass_neighbor_lists
-        self.i = sim.add_symbol(Types.Int32)
-        self.dp = sim.add_symbol(Types.Vector)
-        self.rsq = sim.add_symbol(Types.Double)
-        self.ncases = NeighborFor.number_of_cases(None if bypass_neighbor_lists else sim.neighbor_lists)
-        self.jlist = [sim.add_symbol(Types.Int32) for _ in range(self.ncases)]
-        self.jblocks = [Block(sim, []) for _ in range(self.ncases)]
-        self.current_j = -1
+
+class InteractionData:
+    def __init__(self, sim, shape):
+        self._i = sim.add_symbol(Types.Int32)
+        self._j = sim.add_symbol(Types.Int32)
+        self._delta = sim.add_symbol(Types.Vector)
+        self._squared_distance = sim.add_symbol(Types.Double)
+        self._penetration_depth = sim.add_symbol(Types.Double)
+        self._contact_point = sim.add_symbol(Types.Vector)
+        self._contact_normal = sim.add_symbol(Types.Vector)
+        self._shape = shape
+
+    def i(self):
+        return self._i
+
+    def j(self):
+        return self._j
 
     def delta(self):
-        return self.dp
+        return self._delta
 
     def squared_distance(self):
-        return self.rsq
+        return self._squared_distance
+
+    def penetration_depth(self):
+        return self._penetration_depth
+
+    def contact_point(self):
+        return self._contact_point
+
+    def contact_normal(self):
+        return self._contact_normal
+
+    def shape(self):
+        return self._shape
+
+
+class ParticleInteraction(Lowerable):
+    def __init__(self, sim, nbody, cutoff_radius, use_cell_lists=False):
+        super().__init__(sim)
+        self.nbody = nbody
+        self.cutoff_radius = cutoff_radius
+        self.use_cell_lists = use_cell_lists
+        self.ncases = NeighborFor.number_of_cases(sim, use_cell_lists)
+        self.interactions_data = []
+        self.blocks = [Block(sim, []) for _ in range(sim.max_shapes() * self.ncases)]
+        self.active_block = None
 
     def add_statement(self, stmt):
-        self.jblocks[self.current_j].add_statement(stmt)
+        self.active_block.add_statement(stmt)
 
     def __iter__(self):
         self.sim.add_statement(self)
         self.sim.enter(self)
 
         # Neighbors vary across iterations
-        for j in range(self.ncases):
-            self.current_j = j
-            yield self.i, self.jlist[j]
+        for shape in range(self.sim.max_shapes()):
+            for case in range(self.ncases):
+                self.active_block = self.blocks[shape * self.ncases + case]
+                self.interactions_data.append(InteractionData(self.sim, shape))
+                yield self.interactions_data[-1]
 
         self.sim.leave()
-        self.current_j = -1
+        self.active_block = None
 
     @pairs_inline
     def lower(self):
         if self.nbody == 2:
             position = self.sim.position()
+            radius = self.sim.property('radius')
+            normal = self.sim.property('normal')
             cell_lists = self.sim.cell_lists
-            neighbor_lists = None if self.bypass_neighbor_lists else self.sim.neighbor_lists
-            for i in ParticleFor(self.sim):
-                j = 0
-                self.i.assign(i)
+            neighbor_lists = None if self.use_cell_lists else self.sim.neighbor_lists
 
+            for i in ParticleFor(self.sim):
+                interaction = 0
                 for neigh in NeighborFor(self.sim, i, cell_lists, neighbor_lists):
-                    dp = position[i] - position[neigh.particle_index()]
-                    rsq = dp.x() * dp.x() + dp.y() * dp.y() + dp.z() * dp.z()
-                    self.jlist[j].assign(neigh.particle_index())
-                    self.dp.assign(dp)
-                    self.rsq.assign(rsq)
-                    self.sim.add_statement(Filter(self.sim, rsq < self.cutoff_radius, self.jblocks[j]))
-                    j += 1
+                    interaction_data = self.interactions_data[interaction]
+                    shape = interaction_data.shape()
+                    j = neigh.particle_index()
+
+                    if shape == Shapes.Sphere:
+                        delta = position[i] - position[j]
+                        squared_distance = delta.x() * delta.x() + delta.y() * delta.y() + delta.z() * delta.z()
+                        separation_dist = radius[i] + radius[j] + self.cutoff_radius
+                        cutoff_condition = squared_distance < separation_dist * separation_dist
+                        distance = Sqrt(self.sim, squared_distance)
+                        penetration_depth = distance - radius[i] - radius[j]
+                        contact_normal = delta * (1.0 / distance)
+                        k = radius[j] + 0.5 * penetration_depth
+                        contact_point = position[j] + contact_normal * k
+
+                    elif shape == Shapes.Halfspace:
+                        d = normal[j][0] * position[j][0] + normal[j][1] * position[j][1] + normal[j][2] * position[j][2]
+                        k = normal[j][0] * position[i][0] + normal[j][1] * position[i][1] + normal[j][2] * position[i][2]
+                        penetration_depth = k - radius[i] - d
+                        cutoff_condition = penetration_depth < self.cutoff_radius
+                        tmp = radius[i] + penetration_depth
+                        contact_normal = normal[j]
+                        contact_point = position[i] - Vector(self.sim, [tmp, tmp, tmp]) * normal[j]
+
+                    else:
+                        raise Exception("Invalid shape!")
+
+                    interaction_data.i().assign(i)
+                    interaction_data.j().assign(j)
+                    interaction_data.delta().assign(delta)
+                    interaction_data.squared_distance().assign(squared_distance)
+                    interaction_data.penetration_depth().assign(penetration_depth)
+                    interaction_data.contact_point().assign(contact_point)
+                    interaction_data.contact_normal().assign(contact_normal)
+                    self.sim.add_statement(Filter(self.sim, cutoff_condition, self.blocks[interaction]))
+                    interaction += 1
 
         else:
-            raise Exception("Interactions among more than two particles is currently not implemented!")
+            raise Exception("Interactions among more than two particles are currently not supported.")
diff --git a/src/pairs/sim/neighbor_lists.py b/src/pairs/sim/neighbor_lists.py
index bd39529b3f28fe2d5cf4d0ad98f5c90d5c4ff750..e54ed69b6d480dec0be7003a4e52ba6e1965a3fa 100644
--- a/src/pairs/sim/neighbor_lists.py
+++ b/src/pairs/sim/neighbor_lists.py
@@ -13,7 +13,7 @@ class NeighborLists:
         self.sim = cell_lists.sim
         self.cell_lists = cell_lists
         self.neighborlists = self.sim.add_array('neighborlists', [self.sim.particle_capacity, self.sim.neighbor_capacity], Types.Int32)
-        self.numneighs = self.sim.add_array('numneighs', self.sim.particle_capacity, Types.Int32)
+        self.numneighs = self.sim.add_array('numneighs', [self.sim.particle_capacity, self.sim.max_shapes()], Types.Int32)
 
 
 class BuildNeighborLists(Lowerable):
@@ -28,13 +28,18 @@ class BuildNeighborLists(Lowerable):
         cell_lists = neighbor_lists.cell_lists
         cutoff_radius = cell_lists.cutoff_radius
         position = sim.position()
-        sim.module_name("neighbor_lists_build")
+        sim.module_name("build_neighbor_lists")
         sim.check_resize(sim.neighbor_capacity, neighbor_lists.numneighs)
 
         for i in ParticleFor(sim):
-            Assign(self.sim, neighbor_lists.numneighs[i], 0)
+            for shape in range(sim.max_shapes()):
+                Assign(sim, neighbor_lists.numneighs[i][shape], 0)
 
-        for i, j in ParticleInteraction(sim, 2, cutoff_radius, bypass_neighbor_lists=True):
-            numneighs = neighbor_lists.numneighs[i]
-            Assign(self.sim, neighbor_lists.neighborlists[i][numneighs], j)
-            Assign(self.sim, neighbor_lists.numneighs[i], numneighs + 1)
+        for interaction_data in ParticleInteraction(sim, 2, cutoff_radius, use_cell_lists=True):
+            i = interaction_data.i()
+            j = interaction_data.j()
+            shape = interaction_data.shape()
+
+            numneighs = neighbor_lists.numneighs[i][shape]
+            Assign(sim, neighbor_lists.neighborlists[i][numneighs], j)
+            Assign(sim, neighbor_lists.numneighs[i][shape], numneighs + 1)
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index 56e30c330b2c9009589a2c02309eabb1fee3954f..328b9619599c530d1145ee9986b513811814d752 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -105,13 +105,13 @@ class Simulation:
     def ndims(self):
         return self.dims
 
-    def add_property(self, prop_name, prop_type, value=0.0, vol=False):
+    def add_property(self, prop_name, prop_type, value=0.0, volatile=False):
         assert self.property(prop_name) is None, f"Property already defined: {prop_name}"
-        return self.properties.add(prop_name, prop_type, value, vol)
+        return self.properties.add(prop_name, prop_type, value, volatile)
 
-    def add_position(self, prop_name, value=[0.0, 0.0, 0.0], vol=False, layout=Layouts.AoS):
+    def add_position(self, prop_name, value=[0.0, 0.0, 0.0], volatile=False, layout=Layouts.AoS):
         assert self.property(prop_name) is None, f"Property already defined: {prop_name}"
-        self.position_prop = self.properties.add(prop_name, Types.Vector, value, vol, layout)
+        self.position_prop = self.properties.add(prop_name, Types.Vector, value, volatile, layout)
         return self.position_prop
 
     def add_feature(self, feature_name, nkinds):
diff --git a/src/pairs/transformations/expressions.py b/src/pairs/transformations/expressions.py
index f82b539b5f7f677729a996495c6513d85af1f68d..c5d8c20bcc14f390c6420697ed21d6e6c65d4d1f 100644
--- a/src/pairs/transformations/expressions.py
+++ b/src/pairs/transformations/expressions.py
@@ -10,7 +10,16 @@ class ReplaceSymbols(Mutator):
         super().__init__(ast)
 
     def mutate_Symbol(self, ast_node):
-        return ast_node.assign_to
+        expr = self.mutate(ast_node.assign_to)
+        for i in ast_node.indexes_to_generate():
+            expr.add_index_to_generate(i)
+
+        return expr
+
+    def mutate_SymbolAccess(self, ast_node):
+        symbol = self.mutate(ast_node.symbol())
+        index = self.mutate(ast_node.index())
+        return symbol[index]
 
 
 class LowerNeighborIndexes(Mutator):