diff --git a/examples/dem.py b/examples/dem.py
index 661218ab6c56a93cd3db668beb6f93aaccb8aa79..f2a60d5b71728172186c87c9fc85cefa4784d002 100644
--- a/examples/dem.py
+++ b/examples/dem.py
@@ -11,7 +11,9 @@ def linear_spring_dashpot(i, j):
     k = radius[j] + 0.5 * penetration_depth
     contact_point = position[j] + contact_normal * k
 
-    rel_vel = -velocity_wf[i] - velocity_wf[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_t = rel_vel - rel_vel_n
 
@@ -57,8 +59,8 @@ def linear_spring_dashpot(i, j):
 
 
 def euler(i):
-    velocity[i] += dt * force[i] / mass[i]
-    position[i] += dt * velocity[i]
+    linear_velocity[i] += dt * force[i] / mass[i]
+    position[i] += dt * linear_velocity[i]
 
 
 def gravity(i):
@@ -90,8 +92,8 @@ gravity_SI = 1.0
 psim = pairs.simulation("dem", debug=True)
 psim.add_position('position')
 psim.add_property('mass', pairs.double(), 1.0)
-psim.add_property('velocity', pairs.vector())
-psim.add_property('velocity_wf', pairs.vector())
+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('radius', pairs.double(), 1.0)
 psim.add_feature('type', ntypes)
@@ -105,7 +107,7 @@ psim.add_contact_property('is_sticking', pairs.int32(), 0)
 psim.add_contact_property('tangential_spring_displacement', pairs.vector(), [0.0, 0.0, 0.0])
 psim.add_contact_property('impact_velocity_magnitude', pairs.double(), 0.0)
 
-psim.read_particle_data("data/fluidized_bed.input", ['mass', 'position', 'velocity'])
+psim.read_particle_data("data/fluidized_bed.input", ['mass', 'position', 'linear_velocity'])
 psim.build_neighbor_lists(cutoff_radius + skin)
 psim.vtk_output(f"output/test_{target}")
 psim.compute(linear_spring_dashpot, cutoff_radius, symbols={'dt': dt})
diff --git a/examples/lj.py b/examples/lj.py
index ef2ff66cf6d9b21aa35c88bb1de5d460de150356..dab50c7457308f902899efbd25b257eaa341186f 100644
--- a/examples/lj.py
+++ b/examples/lj.py
@@ -9,8 +9,8 @@ def lj(i, j):
 
 
 def euler(i):
-    velocity[i] += dt * force[i] / mass[i]
-    position[i] += dt * velocity[i]
+    linear_velocity[i] += dt * force[i] / mass[i]
+    position[i] += dt * linear_velocity[i]
 
 
 cmd = sys.argv[0]
@@ -29,12 +29,12 @@ sigma6 = sigma ** 6
 psim = pairs.simulation("lj", debug=True)
 psim.add_position('position')
 psim.add_property('mass', pairs.double(), 1.0)
-psim.add_property('velocity', pairs.vector())
+psim.add_property('linear_velocity', pairs.vector())
 psim.add_property('force', pairs.vector(), vol=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)])
-psim.read_particle_data("data/minimd_setup_32x32x32.input", ['type', 'mass', 'position', 'velocity', 'particle_flags'])
+psim.read_particle_data("data/minimd_setup_32x32x32.input", ['type', 'mass', 'position', 'linear_velocity', 'particle_flags'])
 psim.build_neighbor_lists(cutoff_radius + skin)
 psim.vtk_output(f"output/test_{target}")
 psim.compute(lj, cutoff_radius)
diff --git a/src/pairs/analysis/blocks.py b/src/pairs/analysis/blocks.py
index 14d36458eae1e0b4a42284e24416a2d4cb1a6304..611029df6e5f062ad2660ed45efd880d35d0ff65 100644
--- a/src/pairs/analysis/blocks.py
+++ b/src/pairs/analysis/blocks.py
@@ -189,6 +189,10 @@ class DetermineExpressionsOwnership(Visitor):
         self.visit_children(ast_node)
         self.update_ownership(ast_node)
 
+    def visit_Vector(self, ast_node):
+        self.visit_children(ast_node)
+        self.update_ownership(ast_node)
+
     def visit_VectorOp(self, ast_node):
         self.visit_children(ast_node)
         self.update_ownership(ast_node)
diff --git a/src/pairs/analysis/expressions.py b/src/pairs/analysis/expressions.py
index 96990e6557bd0f892c3fb57c07021b4d69139509..8f7b8d3e871e4573d7b8c56ca4dd9dfb9f04e313 100644
--- a/src/pairs/analysis/expressions.py
+++ b/src/pairs/analysis/expressions.py
@@ -46,6 +46,9 @@ class DetermineExpressionsTerminals(Visitor):
     def visit_Select(self, ast_node):
         self.traverse_expression(ast_node)
 
+    def visit_Vector(self, ast_node):
+        self.traverse_expression(ast_node)
+
     def visit_VectorOp(self, ast_node):
         self.traverse_expression(ast_node)
 
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index 2b6598dd76dc2c9c0d549c0dc9e60a74cc099090..2d2a2934af606a62d96de8a185925ff9a1cb6945 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -24,7 +24,7 @@ from pairs.ir.sizeof import Sizeof
 from pairs.ir.types import Types
 from pairs.ir.utils import Print
 from pairs.ir.variables import Var, VarDecl, Deref
-from pairs.ir.vectors import VectorAccess, VectorOp, ZeroVector
+from pairs.ir.vectors import Vector, VectorAccess, VectorOp, ZeroVector
 from pairs.sim.timestep import Timestep
 from pairs.code_gen.printer import Printer
 
@@ -365,6 +365,12 @@ class CGen:
                 tkw = Types.c_keyword(math_func.type())
                 self.print(f"const {tkw} {acc_ref} = {math_func.function_name()}({params});")
 
+            if isinstance(ast_node.elem, Vector):
+                vector = ast_node.elem
+                for dim in vector.indexes_to_generate():
+                    expr = self.generate_expression(vector.get_value(dim))
+                    self.print(f"const double v{vector.id()}_{dim} = {expr};")
+
             if isinstance(ast_node.elem, VectorOp):
                 vector_op = ast_node.elem
                 for dim in vector_op.indexes_to_generate():
@@ -861,6 +867,10 @@ class CGen:
         if isinstance(ast_node, VectorAccess):
             return self.generate_expression(ast_node.expr, mem, self.generate_expression(ast_node.index))
 
+        if isinstance(ast_node, Vector):
+            assert index is not None, "Index must be set for vector."
+            return f"v{ast_node.id()}_{index}"
+
         if isinstance(ast_node, VectorOp):
             assert index is not None, "Index must be set for vector operation."
             return f"e{ast_node.id()}_{index}"
diff --git a/src/pairs/ir/arrays.py b/src/pairs/ir/arrays.py
index ef785b4956e8c0ed5ac5af93b1883ecf76c1141a..309f87605257d329d10cbc2bbc88349b170d725f 100644
--- a/src/pairs/ir/arrays.py
+++ b/src/pairs/ir/arrays.py
@@ -1,5 +1,4 @@
 from functools import reduce
-from pairs.ir.assign import Assign
 from pairs.ir.ast_node import ASTNode
 from pairs.ir.ast_term import ASTTerm
 from pairs.ir.scalars import ScalarOp
@@ -188,12 +187,6 @@ class ArrayAccess(ASTTerm):
 
         return False
 
-    def set(self, other):
-        return self.sim.add_statement(Assign(self.sim, self, other))
-
-    def add(self, other):
-        return self.sim.add_statement(Assign(self.sim, self, self + other))
-
     def id(self):
         return self.acc_id
 
diff --git a/src/pairs/ir/assign.py b/src/pairs/ir/assign.py
index 9c492725aa3919dc0ce045f0c585ef4761066c62..30f1a1b5e85ffa7091e6438f69f6672142e78185 100644
--- a/src/pairs/ir/assign.py
+++ b/src/pairs/ir/assign.py
@@ -1,4 +1,5 @@
 from pairs.ir.ast_node import ASTNode
+from pairs.ir.ast_term import ASTTerm
 from pairs.ir.lit import Lit
 from pairs.ir.types import Types
 
@@ -11,13 +12,15 @@ class Assign(ASTNode):
 
         # When vector assignments occur, all indexes for the dest
         # and source terms must be generated
-        if self._dest.type() == Types.Vector:
+        if isinstance(self._dest, ASTTerm) and self._dest.type() == Types.Vector:
             for dim in range(sim.ndims()):
                 self._dest.add_index_to_generate(dim)
 
-                if self._src.type() == Types.Vector:
+                if isinstance(self._src, ASTTerm) and self._src.type() == Types.Vector:
                     self._src.add_index_to_generate(dim)
 
+        sim.add_statement(self)
+
     def __str__(self):
         return f"Assign<{self._dest, self._src}>"
 
diff --git a/src/pairs/ir/lit.py b/src/pairs/ir/lit.py
index d23367c82fd5218f1646c33eed9bf827b993e9c3..03484656fbe88d66f0eb49e5913eb4f29cbc8339 100644
--- a/src/pairs/ir/lit.py
+++ b/src/pairs/ir/lit.py
@@ -12,23 +12,16 @@ class Lit(ASTNode):
     def __init__(self, sim, value):
         super().__init__(sim)
         self.value = value
-        self.lit_type = Types.Invalid
 
-        if isinstance(value, int):
-            self.lit_type = Types.Int32
-
-        if isinstance(value, float):
-            self.lit_type = Types.Double
-
-        if isinstance(value, bool):
-            self.lit_type = Types.Boolean
-
-        if isinstance(value, str):
-            self.lit_type = Types.String
-
-        if isinstance(value, list):
-            self.lit_type = Types.Vector
+        type_mapping = {
+            int: Types.Int32,
+            float: Types.Double,
+            bool: Types.Boolean,
+            str: Types.String,
+            list: Types.Vector
+        }
 
+        self.lit_type = type_mapping.get(type(value), Types.Invalid)
         assert self.lit_type != Types.Invalid, "Invalid literal type!"
 
     def __str__(self):
diff --git a/src/pairs/ir/mutator.py b/src/pairs/ir/mutator.py
index 73c9bdce6a86a5279276fa942179cf68264303ad..ab54c9f23cec6530e245408fb6d72db638255f53 100644
--- a/src/pairs/ir/mutator.py
+++ b/src/pairs/ir/mutator.py
@@ -189,6 +189,10 @@ class Mutator:
         ast_node.block = self.mutate(ast_node.block)
         return ast_node
 
+    def mutate_Vector(self, ast_node):
+        ast_node._values = [self.mutate(v) for v in ast_node._values]
+        return ast_node
+
     def mutate_VectorAccess(self, ast_node):
         ast_node.expr = self.mutate(ast_node.expr)
         return ast_node
diff --git a/src/pairs/ir/properties.py b/src/pairs/ir/properties.py
index c36c8996cb79b8d057af1d6740f30905457d156a..712156ef1c7a0fbd06b50b0e0b259a9e2e5a2542 100644
--- a/src/pairs/ir/properties.py
+++ b/src/pairs/ir/properties.py
@@ -1,6 +1,5 @@
 from pairs.ir.ast_node import ASTNode
 from pairs.ir.ast_term import ASTTerm
-from pairs.ir.assign import Assign
 from pairs.ir.declaration import Decl
 from pairs.ir.layouts import Layouts
 from pairs.ir.lit import Lit
@@ -131,15 +130,6 @@ class PropertyAccess(ASTTerm):
         self.inlined = True
         return self
 
-    def set(self, other):
-        return self.sim.add_statement(Assign(self.sim, self, other))
-
-    def add(self, other):
-        return self.sim.add_statement(Assign(self.sim, self, self + other))
-
-    def sub(self, other):
-        return self.sim.add_statement(Assign(self.sim, self, self - other))
-
     def id(self):
         return self.acc_id
 
@@ -302,15 +292,6 @@ class ContactPropertyAccess(ASTTerm):
         self.inlined = True
         return self
 
-    def set(self, other):
-        return self.sim.add_statement(Assign(self.sim, self, other))
-
-    def add(self, other):
-        return self.sim.add_statement(Assign(self.sim, self, self + other))
-
-    def sub(self, other):
-        return self.sim.add_statement(Assign(self.sim, self, self - other))
-
     def id(self):
         return self.acc_id
 
diff --git a/src/pairs/ir/scalars.py b/src/pairs/ir/scalars.py
index 5848abc2a9e88491839d720cf8c690dc35669f5e..929a559d319f4b295889f021bc9e467329b18641 100644
--- a/src/pairs/ir/scalars.py
+++ b/src/pairs/ir/scalars.py
@@ -1,6 +1,5 @@
 from pairs.ir.ast_node import ASTNode
 from pairs.ir.ast_term import ASTTerm
-from pairs.ir.assign import Assign
 from pairs.ir.lit import Lit
 from pairs.ir.operators import Operators
 from pairs.ir.types import Types
@@ -64,18 +63,6 @@ class ScalarOp(ASTTerm):
     def z(self):
         return self.__getitem__(2)
 
-    def set(self, other):
-        assert self.mem is True, "Invalid assignment: lvalue expected!"
-        return self.sim.add_statement(Assign(self.sim, self, other))
-
-    def add(self, other):
-        assert self.mem is True, "Invalid assignment: lvalue expected!"
-        return self.sim.add_statement(Assign(self.sim, self, self + other))
-
-    def sub(self, other):
-        assert self.mem is True, "Invalid assignment: lvalue expected!"
-        return self.sim.add_statement(Assign(self.sim, self, self - other))
-
     def infer_type(lhs, rhs, op):
         lhs_type = lhs.type()
 
diff --git a/src/pairs/ir/variables.py b/src/pairs/ir/variables.py
index 17af8a5477eecde69a9d39a7b74d0978418fe9ac..b0ad8564fd405c3cac9c7e681f25b584d96819f8 100644
--- a/src/pairs/ir/variables.py
+++ b/src/pairs/ir/variables.py
@@ -26,7 +26,7 @@ class Variables:
         lit = Lit.cvt(self.sim, init)
         tmp_id = Variables.new_temp_id()
         tmp_var = Var(self.sim, f"tmp{tmp_id}", lit.type(), temp=True)
-        self.sim.add_statement(Assign(self.sim, tmp_var, lit))
+        Assign(self.sim, tmp_var, lit)
         return tmp_var
 
     def all(self):
@@ -58,12 +58,6 @@ class Var(ASTTerm):
         # Terminal copies are just themselves
         return self
 
-    def set(self, other):
-        return self.sim.add_statement(Assign(self.sim, self, other))
-
-    def add(self, other):
-        return self.sim.add_statement(Assign(self.sim, self, self + other))
-
     def name(self):
         return self.var_name
 
diff --git a/src/pairs/ir/vectors.py b/src/pairs/ir/vectors.py
index c11f28cc7bf852f462bf2fef28ba88a3443b0890..e6ed7c9b613f7dbb66fdd4288a499153cfe05134 100644
--- a/src/pairs/ir/vectors.py
+++ b/src/pairs/ir/vectors.py
@@ -1,4 +1,3 @@
-from pairs.ir.assign import Assign
 from pairs.ir.ast_node import ASTNode
 from pairs.ir.ast_term import ASTTerm
 from pairs.ir.scalars import ScalarOp
@@ -62,18 +61,6 @@ class VectorOp(ASTTerm):
     def z(self):
         return self.__getitem__(2)
 
-    def set(self, other):
-        assert self.mem is True, "Invalid assignment: lvalue expected!"
-        return self.sim.add_statement(Assign(self.sim, self, other))
-
-    def add(self, other):
-        assert self.mem is True, "Invalid assignment: lvalue expected!"
-        return self.sim.add_statement(Assign(self.sim, self, self + other))
-
-    def sub(self, other):
-        assert self.mem is True, "Invalid assignment: lvalue expected!"
-        return self.sim.add_statement(Assign(self.sim, self, self - other))
-
     def add_terminal(self, terminal):
         self.terminals.add(terminal)
 
@@ -94,44 +81,44 @@ class VectorAccess(ASTTerm):
     def type(self):
         return Types.Double
 
-    def set(self, other):
-        return self.sim.add_statement(Assign(self.sim, self, other))
-
-    def add(self, other):
-        return self.sim.add_statement(Assign(self.sim, self, self + other))
-
-    def sub(self, other):
-        return self.sim.add_statement(Assign(self.sim, self, self - other))
-
     def children(self):
         return [self.expr]
 
 
-class ConstVector(ASTTerm):
+class Vector(ASTTerm):
+    last_vector = 0
+
+    def new_id():
+        Vector.last_vector += 1
+        return Vector.last_vector - 1
+
     def __init__(self, sim, values):
-        assert isinstance(values, list) and len(values) == sim.ndims(), \
-            "ConstVector(): Given list is invalid!"
+        assert isinstance(values, list) and len(values) == sim.ndims(), "Vector(): Given list is invalid!"
         super().__init__(sim, VectorOp)
-        self.values = values
-        self.array = None
+        self._id = Vector.new_id()
+        self._values = values
+        self.terminals = set()
 
     def __str__(self):
-        return f"ConstVector<{self.values}>"
+        return f"Vector<{self.values}>"
 
-    def __getitem__(self, expr_ast):
-        if isinstance(index, Lit) or isinstance(index, int):
-            return self.values[index]
-
-        if self.array is None:
-            self.array = self.sim.add_static_array(f"cv{self.const_vector_id}", self.sim.ndims(), Types.Double)
+    def __getitem__(self, index):
+        return VectorAccess(self.sim, self, Lit.cvt(self.sim, index))
 
-        return self_array[index]
+    def id(self):
+        return self._id
 
     def type(self):
         return Types.Vector
 
+    def get_value(self, dimension):
+        return self._values[dimension]
+
+    def add_terminal(self, terminal):
+        self.terminals.add(terminal)
+
     def children(self):
-        return []
+        return self._values
 
 
 class ZeroVector(ASTTerm):
diff --git a/src/pairs/mapping/funcs.py b/src/pairs/mapping/funcs.py
index 4613748a573f3196b3c1dbb0f07c6c7cd733995d..f18b24365ca25d8a3771e209c996e1b2fefc4251 100644
--- a/src/pairs/mapping/funcs.py
+++ b/src/pairs/mapping/funcs.py
@@ -1,10 +1,11 @@
 import ast
 import inspect
-from pairs.ir.scalars import ScalarOp
+from pairs.ir.assign import Assign
 from pairs.ir.branches import Branch, Filter
 from pairs.ir.lit import Lit
 from pairs.ir.loops import ParticleFor
 from pairs.ir.operators import Operators
+from pairs.ir.scalars import ScalarOp
 from pairs.ir.types import Types
 from pairs.ir.vectors import VectorOp
 from pairs.mapping.keywords import Keywords
@@ -94,7 +95,7 @@ class BuildParticleIR(ast.NodeVisitor):
         if isinstance(lhs, UndefinedSymbol):
             self.add_symbols({lhs.symbol_id: rhs})
         else:
-            lhs.set(rhs)
+            Assign(self.sim, lhs, rhs)
 
     def visit_AugAssign(self, node):
         lhs = self.visit(node.target)
@@ -103,7 +104,7 @@ class BuildParticleIR(ast.NodeVisitor):
         if isinstance(lhs, UndefinedSymbol):
             self.add_symbols({lhs.symbol_id: rhs})
         else:
-            lhs.add(rhs)
+            Assign(self.sim, lhs, lhs + rhs)
 
     def visit_BinOp(self, node):
         #print(ast.dump(node))
diff --git a/src/pairs/mapping/keywords.py b/src/pairs/mapping/keywords.py
index d49b044da25e839e2b2810921207d96042af2a64..9cc16f1c12f11b58cddc0c404d81a53b84cc75ad 100644
--- a/src/pairs/mapping/keywords.py
+++ b/src/pairs/mapping/keywords.py
@@ -5,7 +5,7 @@ from pairs.ir.loops import Continue
 from pairs.ir.math import Sqrt
 from pairs.ir.select import Select
 from pairs.ir.types import Types
-from pairs.ir.vectors import ZeroVector
+from pairs.ir.vectors import Vector, ZeroVector
 
 
 class Keywords:
@@ -62,6 +62,16 @@ class Keywords:
         assert vector2.type() == Types.Vector, "dot(): Second argument must be a vector!"
         return sum([vector1[d] * vector2[d] for d in range(self.sim.ndims())])
 
+    def keyword_cross(self, args):
+        assert len(args) == 2, "cross() keyword requires two parameters!"
+        vector1 = args[0]
+        vector2 = args[1]
+        assert vector1.type() == Types.Vector, "cross(): First argument must be a vector!"
+        assert vector2.type() == Types.Vector, "cross(): Second argument must be a vector!"
+        return Vector(self.sim, [ vector1[1] * vector2[2] - vector1[2] * vector2[1],
+                                  vector1[2] * vector2[0] - vector1[0] * vector2[2],
+                                  vector1[0] * vector2[1] - vector1[1] * vector2[0] ])
+
     def keyword_normalized(self, args):
         assert len(args) == 1, "normalized() keyword requires one parameter!"
         vector = args[0]
diff --git a/src/pairs/sim/cell_lists.py b/src/pairs/sim/cell_lists.py
index 95dfc82de3dc5e2284c6631a492d82b1d99b3ac1..142eae8249189224ae97d5d6e11d607c58095af9 100644
--- a/src/pairs/sim/cell_lists.py
+++ b/src/pairs/sim/cell_lists.py
@@ -1,5 +1,6 @@
 from functools import reduce
 import math
+from pairs.ir.assign import Assign
 from pairs.ir.ast_term import ASTTerm
 from pairs.ir.atomic import AtomicAdd
 from pairs.ir.scalars import ScalarOp
@@ -51,19 +52,19 @@ class CellListsStencilBuild(Lowerable):
         sim.check_resize(cl.ncells_capacity, cl.ncells)
 
         for d in range(sim.ndims()):
-            cl.dim_ncells[d].set(Ceil(sim, (grid.max(d) - grid.min(d)) / cl.spacing[d]) + 2)
+            Assign(sim, cl.dim_ncells[d], Ceil(sim, (grid.max(d) - grid.min(d)) / cl.spacing[d]) + 2)
             ntotal_cells *= cl.dim_ncells[d]
 
-        cl.ncells.set(ntotal_cells + 1)
+        Assign(sim, cl.ncells, ntotal_cells + 1)
         for _ in sim.nest_mode():
-            cl.nstencil.set(0)
+            Assign(sim, cl.nstencil, 0)
             for d in range(sim.ndims()):
                 nneigh = cl.nneighbor_cells[d]
                 for d_idx in For(sim, -nneigh, nneigh + 1):
                     index = (d_idx if index is None else index * cl.dim_ncells[d - 1] + d_idx)
                     if d == sim.ndims() - 1:
-                        cl.stencil[cl.nstencil].set(index)
-                        cl.nstencil.set(cl.nstencil + 1)
+                        Assign(sim, cl.stencil[cl.nstencil], index)
+                        Assign(sim, cl.nstencil, cl.nstencil + 1)
 
 
 class CellListsBuild(Lowerable):
@@ -82,7 +83,7 @@ class CellListsBuild(Lowerable):
         sim.check_resize(cl.cell_capacity, cl.cell_sizes)
 
         for c in For(sim, 0, cl.ncells):
-            cl.cell_sizes[c].set(0)
+            Assign(sim, cl.cell_sizes[c], 0)
 
         for i in ParticleFor(sim, local_only=False):
             flat_index = sim.add_temp_var(0)
@@ -94,9 +95,9 @@ class CellListsBuild(Lowerable):
                 for d in range(sim.ndims()):
                     index_1d = (cell_index[d] if index_1d is None else index_1d * cl.dim_ncells[d] + cell_index[d])
 
-                flat_index.set(index_1d + 1)
+                Assign(sim, flat_index, index_1d + 1)
 
             for _ in Filter(sim, ScalarOp.and_op(flat_index >= 0, flat_index < cl.ncells)):
                 index_in_cell = AtomicAdd(sim, cl.cell_sizes[flat_index], 1)
-                cl.particle_cell[i].set(flat_index)
-                cl.cell_particles[flat_index][index_in_cell].set(i)
+                Assign(sim, cl.particle_cell[i], flat_index)
+                Assign(sim, cl.cell_particles[flat_index][index_in_cell], i)
diff --git a/src/pairs/sim/comm.py b/src/pairs/sim/comm.py
index a2d9f17cda2b5c84f6b0bd225804ce4644005cf7..55b3d6c0514eed7661d7d3356b7aa2cc7daf328e 100644
--- a/src/pairs/sim/comm.py
+++ b/src/pairs/sim/comm.py
@@ -1,3 +1,4 @@
+from pairs.ir.assign import Assign
 from pairs.ir.atomic import AtomicAdd
 from pairs.ir.scalars import ScalarOp
 from pairs.ir.block import pairs_device_block, pairs_host_block, pairs_inline
@@ -44,8 +45,9 @@ class Comm:
     @pairs_inline
     def borders(self):
         prop_list = [self.sim.property(p) for p in ['mass', 'position', 'particle_flags']]
-        self.nsend_all.set(0)
-        self.sim.nghost.set(0)
+        Assign(self.sim, self.nsend_all, 0)
+        Assign(self.sim, self.sim.nghost, 0)
+
         for step in range(self.dom_part.number_of_steps()):
             DetermineGhostParticles(self, step, self.sim.cell_spacing())
             CommunicateSizes(self, step)
@@ -53,20 +55,21 @@ class Comm:
             PackGhostParticles(self, step, prop_list)
             CommunicateData(self, step, prop_list)
             UnpackGhostParticles(self, step, prop_list)
-            self.sim.nghost.add(sum([self.nrecv[j] for j in self.dom_part.step_indexes(step)]))
+            Assign(self.sim, self.sim.nghost, self.sim.nghost + sum([self.nrecv[j] for j in self.dom_part.step_indexes(step)]))
 
     @pairs_inline
     def exchange(self):
-        prop_list = [self.sim.property(p) for p in ['mass', 'position', 'velocity', 'particle_flags']]
+        prop_list = [self.sim.property(p) for p in ['mass', 'position', 'linear_velocity', 'particle_flags']]
         for step in range(self.dom_part.number_of_steps()):
-            self.nsend_all.set(0)
-            self.sim.nghost.set(0)
+            Assign(self.sim, self.nsend_all, 0)
+            Assign(self.sim, self.sim.nghost, 0)
+
             for s in range(step):
                 for j in self.dom_part.step_indexes(s):
-                    self.nsend[j].set(0)
-                    self.nrecv[j].set(0)
-                    self.send_offsets[j].set(0)
-                    self.recv_offsets[j].set(0)
+                    Assign(self.sim, self.nsend[j], 0)
+                    Assign(self.sim, self.nrecv[j], 0)
+                    Assign(self.sim, self.send_offsets[j], 0)
+                    Assign(self.sim, self.recv_offsets[j], 0)
 
             DetermineGhostParticles(self, step, 0.0)
             CommunicateSizes(self, step)
@@ -131,24 +134,24 @@ class DetermineGhostParticles(Lowerable):
         #self.sim.check_resize(self.comm.send_capacity, nsend_all)
 
         for j in self.comm.dom_part.step_indexes(self.step):
-            nsend[j].set(0)
-            nrecv[j].set(0)
+            Assign(self.sim, nsend[j], 0)
+            Assign(self.sim, nrecv[j], 0)
 
         if is_exchange:
             for i in ParticleFor(self.sim):
-                exchg_flag[i].set(0)
+                Assign(self.sim, exchg_flag[i], 0)
 
         for i, j, _, pbc in self.comm.dom_part.ghost_particles(self.step, self.sim.position(), self.spacing):
             next_idx = AtomicAdd(self.sim, nsend_all, 1)
-            send_map[next_idx].set(i)
+            Assign(self.sim, send_map[next_idx], i)
 
             if is_exchange:
-                exchg_flag[i].set(1)
+                Assign(self.sim, exchg_flag[i], 1)
 
             for d in range(self.sim.ndims()):
-                send_mult[next_idx][d].set(pbc[d])
+                Assign(self.sim, send_mult[next_idx][d], pbc[d])
 
-            nsend[j].add(1)
+            Assign(self.sim, nsend[j], nsend[j] + 1)
 
 
 class SetCommunicationOffsets(Lowerable):
@@ -174,8 +177,8 @@ class SetCommunicationOffsets(Lowerable):
                 irecv += nrecv[j]
 
         for j in self.comm.dom_part.step_indexes(self.step):
-            send_offsets[j].set(isend)
-            recv_offsets[j].set(irecv)
+            Assign(self.sim, send_offsets[j], isend)
+            Assign(self.sim, recv_offsets[j], irecv)
             isend += nsend[j]
             irecv += nrecv[j]
 
@@ -212,13 +215,13 @@ class PackGhostParticles(Lowerable):
                         if p == self.sim.position():
                             src += send_mult[i][d] * self.sim.grid.length(d)
 
-                        send_buffer[i][p_offset + d].set(src)
+                        Assign(self.sim, send_buffer[i][p_offset + d], src)
 
                     p_offset += self.sim.ndims()
 
                 else:
                     cast_fn = lambda x: Cast(self.sim, x, Types.Double) if p.type() != Types.Double else x
-                    send_buffer[i][p_offset].set(cast_fn(p[m]))
+                    Assign(self.sim, send_buffer[i][p_offset], cast_fn(p[m]))
                     p_offset += 1
 
             
@@ -248,13 +251,13 @@ class UnpackGhostParticles(Lowerable):
             for p in self.prop_list:
                 if p.type() == Types.Vector:
                     for d in range(self.sim.ndims()):
-                        p[nlocal + i][d].set(recv_buffer[i][p_offset + d])
+                        Assign(self.sim, p[nlocal + i][d], recv_buffer[i][p_offset + d])
 
                     p_offset += self.sim.ndims()
 
                 else:
                     cast_fn = lambda x: Cast(self.sim, x, p.type()) if p.type() != Types.Double else x
-                    p[nlocal + i].set(cast_fn(recv_buffer[i][p_offset]))
+                    Assign(self.sim, p[nlocal + i], cast_fn(recv_buffer[i][p_offset]))
                     p_offset += 1
 
 
@@ -273,13 +276,13 @@ class RemoveExchangedParticles_part1(Lowerable):
             for need_copy in Branch(self.sim, particle_id < self.sim.nlocal - self.comm.nsend_all):
                 if need_copy:
                     for _ in While(self.sim, ScalarOp.cmp(self.comm.exchg_flag[send_pos], 1)):
-                        send_pos.set(send_pos - 1)
+                        Assign(self.sim, send_pos, send_pos - 1)
 
-                    self.comm.exchg_copy_to[i].set(send_pos)
-                    send_pos.set(send_pos - 1)
+                    Assign(self.sim, self.comm.exchg_copy_to[i], send_pos)
+                    Assign(self.sim, send_pos, send_pos - 1)
 
                 else:
-                    self.comm.exchg_copy_to[i].set(-1)
+                    Assign(self.sim, self.comm.exchg_copy_to[i], -1)
 
 
 class RemoveExchangedParticles_part2(Lowerable):
@@ -300,12 +303,12 @@ class RemoveExchangedParticles_part2(Lowerable):
                 for p in self.prop_list:
                     if p.type() == Types.Vector:
                         for d in range(self.sim.ndims()):
-                            p[dst][d].set(p[src][d])
+                            Assign(self.sim, p[dst][d], p[src][d])
 
                     else:
-                        p[dst].set(p[src])
+                        Assign(self.sim, p[dst], p[src])
 
-        self.sim.nlocal.set(self.sim.nlocal - self.comm.nsend_all)
+        Assign(self.sim, self.sim.nlocal, self.sim.nlocal - self.comm.nsend_all)
 
 
 class ChangeSizeAfterExchange(Lowerable):
@@ -317,7 +320,6 @@ class ChangeSizeAfterExchange(Lowerable):
 
     @pairs_host_block
     def lower(self):
-        sim = self.sim
-        sim.module_name(f"change_size_after_exchange{self.step}")
-        sim.check_resize(self.sim.particle_capacity, self.sim.nlocal)
-        self.sim.nlocal.set(sim.nlocal + sum([self.comm.nrecv[j] for j in self.comm.dom_part.step_indexes(self.step)]))
+        self.sim.module_name(f"change_size_after_exchange{self.step}")
+        self.sim.check_resize(self.sim.particle_capacity, self.sim.nlocal)
+        Assign(self.sim, self.sim.nlocal, self.sim.nlocal + sum([self.comm.nrecv[j] for j in self.comm.dom_part.step_indexes(self.step)]))
diff --git a/src/pairs/sim/contact_history.py b/src/pairs/sim/contact_history.py
index 818b5596ab5f667b907aa56d78ca8336fdca8edb..6541e71a8a5bf2cdb6774eb655c66f9df2a32bed 100644
--- a/src/pairs/sim/contact_history.py
+++ b/src/pairs/sim/contact_history.py
@@ -1,7 +1,8 @@
-from pairs.ir.scalars import ScalarOp
+from pairs.ir.assign import Assign
 from pairs.ir.block import pairs_device_block
 from pairs.ir.branches import Branch, Filter
 from pairs.ir.loops import ParticleFor, For
+from pairs.ir.scalars import ScalarOp
 from pairs.ir.types import Types
 from pairs.ir.utils import Print
 from pairs.sim.interaction import NeighborFor
@@ -36,33 +37,33 @@ class BuildContactHistory(Lowerable):
 
             for neigh in NeighborFor(self.sim, i, cell_lists, neighbor_lists):
                 j = neigh.particle_index()
-                contact_index.set(-1)
+                Assign(self.sim, contact_index, -1)
                 for k in For(self.sim, 0, num_contacts[i]):
                     for _ in Filter(self.sim, ScalarOp.cmp(contact_lists[i][k], j)):
-                        contact_index.set(k)
+                        Assign(self.sim, contact_index, k)
 
                 for _ in Filter(self.sim,
                                 ScalarOp.and_op(contact_index >= 0,
                                                 ScalarOp.neq(last_contact, contact_index))):
                     for contact_prop in self.sim.contact_properties:
                         tmp = self.sim.add_temp_var(contact_prop[i, last_contact])
-                        contact_prop[i, last_contact].set(contact_prop[i, contact_index])
-                        contact_prop[i, contact_index].set(tmp)
+                        Assign(self.sim, contact_prop[i, last_contact], contact_prop[i, contact_index])
+                        Assign(self.sim, contact_prop[i, contact_index], tmp)
 
-                    contact_lists[i][contact_index].set(contact_lists[i][last_contact])
-                    contact_lists[i][last_contact].set(j)
+                    Assign(self.sim, contact_lists[i][contact_index], contact_lists[i][last_contact])
+                    Assign(self.sim, contact_lists[i][last_contact], j)
 
-                last_contact.set(last_contact + 1)
+                Assign(self.sim, last_contact, last_contact + 1)
 
-            last_contact.set(0)
+            Assign(self.sim, last_contact, 0)
             for neigh in NeighborFor(self.sim, i, cell_lists, neighbor_lists):
                 j = neigh.particle_index()
                 for _ in Filter(self.sim, ScalarOp.neq(contact_lists[i][last_contact], j)):
                     for contact_prop in self.sim.contact_properties:
-                        contact_prop[i, last_contact].set(contact_prop.default())
+                        Assign(self.sim, contact_prop[i, last_contact], contact_prop.default())
 
-                    contact_lists[i][last_contact].set(j)
+                    Assign(self.sim, contact_lists[i][last_contact], j)
 
-                last_contact.set(last_contact + 1)
+                Assign(self.sim, last_contact, last_contact + 1)
 
-            num_contacts[i].set(last_contact)
+            Assign(self.sim, num_contacts[i], last_contact)
diff --git a/src/pairs/sim/lattice.py b/src/pairs/sim/lattice.py
index 9ab1f9b0549c0976211395c74a3e63a05edc545b..7e2cbec9131fdfdbdf3415eb66dca15281d8d739 100644
--- a/src/pairs/sim/lattice.py
+++ b/src/pairs/sim/lattice.py
@@ -1,3 +1,4 @@
+from pairs.ir.assign import Assign
 from pairs.ir.block import pairs_inline
 from pairs.ir.loops import For
 from pairs.ir.types import Types
@@ -32,17 +33,17 @@ class ParticleLattice(Lowerable):
 
                         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)
+                            Assign(self.sim, self.positions[index][d_], 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() == Types.Vector:
                                 for d_ in range(0, self.sim.ndims()):
-                                    prop[index][d_].set(prop.default()[d_])
+                                    Assign(self.sim, prop[index][d_], prop.default()[d_])
 
                             else:
-                                prop[index].set(prop.default())
+                                Assign(self.sim, prop[index], prop.default())
 
-                        self.sim.nlocal.set(self.sim.nlocal + 1)
+                        Assign(self.sim, self.sim.nlocal, self.sim.nlocal + 1)
 
         return self.sim.block
diff --git a/src/pairs/sim/neighbor_lists.py b/src/pairs/sim/neighbor_lists.py
index 04b48e533278aa6a792f9a1905fac736c4360a9e..f79a4d1a836a4579112acd7f436365220ed56d28 100644
--- a/src/pairs/sim/neighbor_lists.py
+++ b/src/pairs/sim/neighbor_lists.py
@@ -1,3 +1,4 @@
+from pairs.ir.assign import Assign
 from pairs.ir.block import pairs_device_block
 from pairs.ir.branches import Branch, Filter
 from pairs.ir.loops import ParticleFor
@@ -31,9 +32,9 @@ class NeighborListsBuild(Lowerable):
         sim.check_resize(sim.neighbor_capacity, neighbor_lists.numneighs)
 
         for i in ParticleFor(sim):
-            neighbor_lists.numneighs[i].set(0)
+            Assign(self.sim, neighbor_lists.numneighs[i], 0)
 
         for i, j in ParticleInteraction(sim, 2, cutoff_radius, bypass_neighbor_lists=True):
             numneighs = neighbor_lists.numneighs[i]
-            neighbor_lists.neighborlists[i][numneighs].set(j)
-            neighbor_lists.numneighs[i].set(numneighs + 1)
+            Assign(self.sim, neighbor_lists.neighborlists[i][numneighs], j)
+            Assign(self.sim, neighbor_lists.numneighs[i], numneighs + 1)
diff --git a/src/pairs/sim/pbc.py b/src/pairs/sim/pbc.py
index 47bb6676623c5d27d44ef27ba4ad0f56c740bc5c..277c18f3f373e65861e282c7979d87fe894c25ac 100644
--- a/src/pairs/sim/pbc.py
+++ b/src/pairs/sim/pbc.py
@@ -1,3 +1,4 @@
+from pairs.ir.assign import Assign
 from pairs.ir.block import pairs_device_block, pairs_host_block
 from pairs.ir.branches import Branch, Filter
 from pairs.ir.loops import For, ParticleFor
@@ -23,7 +24,7 @@ class EnforcePBC(Lowerable):
             # TODO: VecFilter?
             for d in range(0, ndims):
                 for _ in Filter(sim, positions[i][d] < grid.min(d)):
-                    positions[i][d].add(grid.length(d))
+                    Assign(sim, positions[i][d], positions[i][d] + grid.length(d)) 
 
                 for _ in Filter(sim, positions[i][d] > grid.max(d)):
-                    positions[i][d].sub(grid.length(d))
+                    Assign(sim, positions[i][d], positions[i][d] - grid.length(d)) 
diff --git a/src/pairs/sim/properties.py b/src/pairs/sim/properties.py
index ec4d7dd28acf85f210239e848326454dc836ba33..2f89336976b433c0a14182321355a13910243a96 100644
--- a/src/pairs/sim/properties.py
+++ b/src/pairs/sim/properties.py
@@ -1,3 +1,4 @@
+from pairs.ir.assign import Assign
 from pairs.ir.block import pairs_device_block, pairs_inline
 from pairs.ir.loops import ParticleFor
 from pairs.ir.memory import Malloc, Realloc
@@ -58,4 +59,4 @@ class PropertiesResetVolatile(Lowerable):
         self.sim.module_name("reset_volatile_properties")
         for i in ParticleFor(self.sim):
             for p in self.sim.properties.volatiles():
-                p[i].set(0.0)
+                Assign(self.sim, p[i], 0.0)
diff --git a/src/pairs/sim/read_from_file.py b/src/pairs/sim/read_from_file.py
index c447f8bef305d68c07736d3a381fd05a66d92430..ec2497b27c7746bface233c7610a24cf9dc5e99d 100644
--- a/src/pairs/sim/read_from_file.py
+++ b/src/pairs/sim/read_from_file.py
@@ -1,3 +1,4 @@
+from pairs.ir.assign import Assign
 from pairs.ir.block import pairs_inline
 from pairs.ir.functions import Call_Int, Call_Void
 from pairs.ir.particle_attributes import ParticleAttributeList
@@ -18,14 +19,14 @@ class ReadParticleData(Lowerable):
     def lower(self):
         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])
+            Assign(self.sim, self.grid.min(d), self.grid_buffer[d * 2 + 0])
+            Assign(self.sim, self.grid.max(d), 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.attrs, self.attrs.length()]))
+        Assign(self.sim, self.sim.nlocal, Call_Int(self.sim, "pairs::read_particle_data", [self.filename, self.attrs, self.attrs.length()]))
 
 
 class ReadFeatureData(Lowerable):
diff --git a/src/pairs/transformations/expressions.py b/src/pairs/transformations/expressions.py
index 8787ade5cbf5c16d2db97f60fef2059f79512d30..f82b539b5f7f677729a996495c6513d85af1f68d 100644
--- a/src/pairs/transformations/expressions.py
+++ b/src/pairs/transformations/expressions.py
@@ -219,6 +219,15 @@ class AddExpressionDeclarations(Mutator):
 
         return ast_node
 
+    def mutate_Vector(self, ast_node):
+        ast_node._values = [self.mutate(v) for v in ast_node._values]
+        vector_id = id(ast_node)
+        if vector_id not in self.declared_exprs and vector_id not in self.params:
+            self.push_decl(Decl(ast_node.sim, ast_node))
+            self.declared_exprs.append(vector_id)
+
+        return ast_node
+
     def mutate_VectorOp(self, ast_node):
         ast_node.lhs = self.mutate(ast_node.lhs)
         if not ast_node.operator().is_unary():
diff --git a/src/pairs/transformations/loops.py b/src/pairs/transformations/loops.py
index 2cf82a4d85b618de84bec3eddeebade34a7acb8d..c7ab965dfa9fcac1c39d9f2e8bb6050029d1c921 100644
--- a/src/pairs/transformations/loops.py
+++ b/src/pairs/transformations/loops.py
@@ -6,7 +6,7 @@ from pairs.ir.mutator import Mutator
 from pairs.ir.properties import PropertyAccess, ContactPropertyAccess
 from pairs.ir.scalars import ScalarOp
 from pairs.ir.select import Select
-from pairs.ir.vectors import VectorOp
+from pairs.ir.vectors import Vector, VectorOp
 
 
 class LICM(Mutator):
@@ -40,6 +40,7 @@ class LICM(Mutator):
             PropertyAccess,
             ScalarOp,
             Select,
+            Vector,
             VectorOp
         )
 
diff --git a/src/pairs/transformations/modules.py b/src/pairs/transformations/modules.py
index 57c26d6ee4dc149323cab23f347c3f328cf24529..87cd6c8538a1fc4c10b014636bcfde48db22d9e1 100644
--- a/src/pairs/transformations/modules.py
+++ b/src/pairs/transformations/modules.py
@@ -84,7 +84,7 @@ class AddResizeLogic(Mutator):
                 capacities = list(self.module_resizes[module].values())
                 resize_id = resizes[capacities.index(match_capacity)]
                 return Branch(ast_node.sim, src + 1 >= match_capacity,
-                              blk_if=Block(ast_node.sim, ast_node.sim.resizes[resize_id].set(src)),
+                              blk_if=Block(ast_node.sim, Assign(ast_node.sim, ast_node.sim.resizes[resize_id], src)),
                               blk_else=Block(ast_node.sim, ast_node))
 
         return ast_node