diff --git a/Makefile b/Makefile
index 11fd99e99b063d7daf4f96484c4c76f8a29e6526..15cf2e1526b682db391fcab7e9a93c1d04b82e5f 100644
--- a/Makefile
+++ b/Makefile
@@ -1,5 +1,6 @@
 .PHONY: all build clean
 
+TESTCASE=dem
 PYCMD=python3
 CC=mpicxx
 NVCC=nvcc
@@ -7,11 +8,13 @@ NVCC_PATH:="$(shell which ${NVCC})"
 CUDA_BIN_PATH:="$(shell dirname ${NVCC_PATH})"
 CUDA_PATH:="$(shell dirname ${CUDA_BIN_PATH})"
 OBJ_PATH=obj
-CPU_SRC=lj.cpp
-GPU_SRC=lj.cu
+CPU_SRC="$(TESTCASE).cpp"
+CPU_BIN="$(TESTCASE)_cpu"
+GPU_SRC="$(TESTCASE).cu"
+GPU_BIN="$(TESTCASE)_gpu"
 DEBUG_FLAGS="-DDEBUG"
 
-all: clean build lj_cpu lj_gpu
+all: clean build $(CPU_BIN) $(GPU_BIN)
 	@echo "Everything was done!"
 
 build:
@@ -19,14 +22,14 @@ build:
 	$(PYCMD) setup.py build && $(PYCMD) setup.py install --user
 
 $(CPU_SRC):
-	@echo "Generating and compiling Lennard-Jones example for CPU..."
+	@echo "Generating and compiling $(TESTCASE) example for CPU..."
 	@mkdir -p $(OBJ_PATH)
-	$(PYCMD) examples/lj.py cpu
+	$(PYCMD) examples/$(TESTCASE).py cpu
 
 $(GPU_SRC):
-	@echo "Generating and compiling Lennard-Jones example for GPU..."
+	@echo "Generating and compiling $(TESTCASE) example for GPU..."
 	@mkdir -p $(OBJ_PATH)
-	$(PYCMD) examples/lj.py gpu
+	$(PYCMD) examples/$(TESTCASE).py gpu
 
 $(OBJ_PATH)/pairs.o: runtime/pairs.cpp
 	$(CC) -c -o $@ $< $(DEBUG_FLAGS)
@@ -38,14 +41,14 @@ $(OBJ_PATH)/dummy.o: runtime/devices/dummy.cpp
 	$(CC) -c -o $@ $< $(DEBUG_FLAGS)
 
 # Targets
-lj_cpu: $(CPU_SRC) $(OBJ_PATH)/pairs.o $(OBJ_PATH)/regular_6d_stencil.o $(OBJ_PATH)/dummy.o
-	$(CC) -O3 -o lj_cpu $(CPU_SRC) $(OBJ_PATH)/pairs.o $(OBJ_PATH)/regular_6d_stencil.o $(OBJ_PATH)/dummy.o -DDEBUG
+$(CPU_BIN): $(CPU_SRC) $(OBJ_PATH)/pairs.o $(OBJ_PATH)/regular_6d_stencil.o $(OBJ_PATH)/dummy.o
+	$(CC) -O3 -o $(CPU_BIN) $(CPU_SRC) $(OBJ_PATH)/pairs.o $(OBJ_PATH)/regular_6d_stencil.o $(OBJ_PATH)/dummy.o -DDEBUG
 
-lj_gpu: $(GPU_SRC) $(OBJ_PATH)/pairs.o $(OBJ_PATH)/regular_6d_stencil.o 
+$(GPU_BIN): $(GPU_SRC) $(OBJ_PATH)/pairs.o $(OBJ_PATH)/regular_6d_stencil.o 
 	$(NVCC) -c -o $(OBJ_PATH)/cuda_runtime.o runtime/devices/cuda.cu -DDEBUG
-	$(NVCC) -c -o $(OBJ_PATH)/lj_gpu.o $(GPU_SRC) -DDEBUG
-	$(CC) -o lj_gpu $(OBJ_PATH)/lj_gpu.o $(OBJ_PATH)/cuda_runtime.o $(OBJ_PATH)/pairs.o $(OBJ_PATH)/regular_6d_stencil.o -lcudart -L$(CUDA_PATH)/lib64
+	$(NVCC) -c -o $(OBJ_PATH)/$(GPU_BIN).o $(GPU_SRC) -DDEBUG
+	$(CC) -o $(GPU_BIN) $(OBJ_PATH)/$(GPU_BIN).o $(OBJ_PATH)/cuda_runtime.o $(OBJ_PATH)/pairs.o $(OBJ_PATH)/regular_6d_stencil.o -lcudart -L$(CUDA_PATH)/lib64
 
 clean:
 	@echo "Cleaning..."
-	rm -rf build lj_cpu lj_gpu lj.o $(CPU_SRC) $(GPU_SRC) dist pairs.egg-info functions functions.pdf $(OBJ_PATH)/pairs.o $(OBJ_PATH)/regular_6d_stencil.o $(OBJ_PATH)/cuda_runtime.o $(OBJ_PATH)/lj_gpu.o
+	rm -rf build $(CPU_BIN) $(GPU_BIN) lj.o $(CPU_SRC) $(GPU_SRC) dist pairs.egg-info functions functions.pdf $(OBJ_PATH)/pairs.o $(OBJ_PATH)/regular_6d_stencil.o $(OBJ_PATH)/cuda_runtime.o $(OBJ_PATH)/$(GPU_BIN).o
diff --git a/examples/dem.py b/examples/dem.py
index 62f26a7f746f3211afced06306b59ec060b383d6..731cfc6022cb0ef136e1cbb99671dab665860c20 100644
--- a/examples/dem.py
+++ b/examples/dem.py
@@ -3,21 +3,26 @@ import sys
 
 
 def linear_spring_dashpot(i, j):
-    delta = -penetration_depth
-    skip_if(delta < 0.0)
+    penetration_depth = rsq - radius[i] - radius[j]
+    skip_when(penetration_depth >= 0.0)
+
+    contact_normal = normalized(delta)
+    k = radius[j] + 0.5 * penetration_depth
+    contact_point = position[j] + contact_normal * k
 
     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
 
-    fN = stiffness_norm[i, j] * delta * contact_normal + damping_norm[i, j] * relVelN;
+    fN = stiffness_norm[i, j] * (-penetration_depth) * contact_normal + damping_norm[i, j] * relVelN;
 
     tan_spring_disp = tangential_spring_displacement[i, j]
     impact_vel_magnitude = impact_velocities_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)
-    new_tan_spring_disp = select(sq_len(rotated_tan_disp) <= 0.0,
+    new_tan_spring_disp = select(square_length(rotated_tan_disp) <= 0.0,
                                  0.0, 
                                  rotated_tan_disp * length(tan_spring_disp) / length(rotated_tan_disp))
     new_tan_spring_disp += dt * rel_vel_t
@@ -39,7 +44,7 @@ def linear_spring_dashpot(i, j):
                              new_tan_spring_disp2)
 
     tangential_spring_displacement[i, j] = n_T_spring_disp
-    impact_velocities_magnitude[i, j] = impact_vel_magnitude
+    impact_velocities_magnitude[i, j] = impact_magnitude
     is_sticking[i, j] = n_sticking
 
     fTabs = min(fTLS_len, f_friction_abs)
@@ -61,30 +66,36 @@ if target != 'cpu' and target != 'gpu':
 dt = 0.005
 cutoff_radius = 2.5
 skin = 0.3
-sigma = 1.0
-epsilon = 1.0
-sigma6 = sigma ** 6
+ntypes = 4
+stiffness_norm = 1.0
+stiffness_tan = 1.0
+damping_norm = 1.0
+damping_tan = 1.0
+friction_static = 1.0
+friction_dynamic = 1.0
 
 psim = pairs.simulation("dem", debug=True)
 psim.add_position('position')
-psim.add_property('mass', Types.Double, 1.0)
-psim.add_property('velocity', Types.Vector)
-psim.add_property('force', Types.Vector, vol=True)
-psim.add_feature('type')
-psim.add_feature_property('type', 'stiffness_norm', Types.Double)
-psim.add_feature_property('type', 'stiffness_tan', Types.Double)
-psim.add_feature_property('type', 'damping_norm', Types.Double)
-psim.add_feature_property('type', 'damping_tan', Types.Double)
-psim.add_feature_property('type', 'friction_static', Types.Double)
-psim.add_feature_property('type', 'friction_dynamic', Types.Double)
-psim.add_contact_history_property('is_sticking', Types.Bool, False)
-psim.add_contact_history_property('tangential_spring_displacement', Types.Vector, [0.0, 0.0, 0.0])
-psim.add_contact_history_property('impact_velocity_magnitude', Types.Double, 0.0)
-
-psim.from_file("data/minimd_setup_32x32x32.input", ['mass', 'position', 'velocity'])
+psim.add_property('mass', pairs.double(), 1.0)
+psim.add_property('velocity', pairs.vector())
+psim.add_property('velocity_wf', pairs.vector())
+psim.add_property('force', pairs.vector(), vol=True)
+psim.add_property('radius', pairs.double(), 1.0)
+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)])
+psim.add_feature_property('type', 'damping_norm', pairs.double(), [damping_norm for i in range(ntypes * ntypes)])
+psim.add_feature_property('type', 'damping_tan', pairs.double(), [damping_tan for i in range(ntypes * ntypes)])
+psim.add_feature_property('type', 'friction_static', pairs.double(), [friction_static for i in range(ntypes * ntypes)])
+psim.add_feature_property('type', 'friction_dynamic', pairs.double(), [friction_dynamic for i in range(ntypes * ntypes)])
+psim.add_contact_property('is_sticking', pairs.int32(), False)
+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.build_neighbor_lists(cutoff_radius + skin)
 psim.vtk_output(f"output/test_{target}")
-psim.compute(lj, cutoff_radius, {'sigma6': sigma6, 'epsilon': epsilon})
+psim.compute(linear_spring_dashpot, cutoff_radius)
 psim.compute(euler, symbols={'dt': dt})
 
 if target == 'gpu':
diff --git a/examples/lj.py b/examples/lj.py
index c7392daff3d0ecfbe37b9b1aa4d31aa4a13c8974..98b2b58c2b18acc130afbeb4f04b84951c9160ac 100644
--- a/examples/lj.py
+++ b/examples/lj.py
@@ -31,7 +31,7 @@ psim.add_position('position')
 psim.add_property('mass', pairs.double(), 1.0)
 psim.add_property('velocity', pairs.vector())
 psim.add_property('force', pairs.vector(), vol=True)
-psim.add_feature('type', 4)
+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'])
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index 7fd9480e77302e4b37e57fb9180f6286451a32cb..46091a6bc1b2c2fa3d17eca6d3eabd2f49ec3e71 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -12,7 +12,7 @@ from pairs.ir.functions import Call
 from pairs.ir.kernel import Kernel, KernelLaunch
 from pairs.ir.layouts import Layouts
 from pairs.ir.lit import Lit
-from pairs.ir.loops import For, Iter, ParticleFor, While
+from pairs.ir.loops import For, Iter, ParticleFor, While, Continue
 from pairs.ir.math import Ceil, Sqrt
 from pairs.ir.memory import Malloc, Realloc
 from pairs.ir.module import ModuleCall
@@ -233,6 +233,9 @@ class CGen:
                 self.generate_statement(stmt)
             self.print.add_indent(-4)
 
+        if isinstance(ast_node, Continue):
+            self.print("continue;")
+
         # TODO: Why there are Decls for other types?
         if isinstance(ast_node, Decl):
             if isinstance(ast_node.elem, BinOp):
@@ -243,13 +246,21 @@ class CGen:
                             lhs = self.generate_expression(bin_op.lhs, bin_op.mem, index=i)
                             rhs = self.generate_expression(bin_op.rhs, index=i)
                             operator = bin_op.operator()
-                            self.print(f"const double e{bin_op.id()}_{i} = {lhs} {operator.symbol()} {rhs};")
+
+                            if operator.is_unary():
+                                self.print(f"const double e{bin_op.id()}_{i} = {operator.symbol()}({lhs});")
+                            else:
+                                self.print(f"const double e{bin_op.id()}_{i} = {lhs} {operator.symbol()} {rhs};")
                     else:
                         lhs = self.generate_expression(bin_op.lhs, bin_op.mem)
                         rhs = self.generate_expression(bin_op.rhs)
                         operator = bin_op.operator()
                         tkw = Types.c_keyword(bin_op.type())
-                        self.print(f"const {tkw} e{bin_op.id()} = {lhs} {operator.symbol()} {rhs};")
+
+                        if operator.is_unary():
+                            self.print(f"const {tkw} e{bin_op.id()} = {operator.symbol()}({lhs});")
+                        else:
+                            self.print(f"const {tkw} e{bin_op.id()} = {lhs} {operator.symbol()} {rhs};")
 
             if isinstance(ast_node.elem, ArrayAccess):
                 array_access = ast_node.elem
@@ -259,6 +270,15 @@ class CGen:
                 acc_ref = f"a{array_access.id()}"
                 self.print(f"const {tkw} {acc_ref} = {array_name}[{acc_index}];")
 
+            if isinstance(ast_node.elem, ContactPropertyAccess):
+                contact_prop_access = ast_node.elem
+                contact_prop = contact_prop_access.contact_prop
+                prop_name = self.generate_expression(contact_prop)
+                tkw = Types.c_keyword(contact_prop_access.type())
+                acc_index = self.generate_expression(contact_prop_access.index)
+                acc_ref = f"cp{feature_prop_access.id()}"
+                self.print(f"const {tkw} {acc_ref} = {prop_name}[{acc_index}];")
+
             if isinstance(ast_node.elem, FeaturePropertyAccess):
                 feature_prop_access = ast_node.elem
                 feature_prop = feature_prop_access.feature_prop
@@ -524,6 +544,32 @@ class CGen:
 
             self.print(f"pairs->addProperty({p.id()}, \"{p.name()}\", &{ptr}, {d_ptr}, {ptype}, {playout}, {sizes});")
 
+        if isinstance(ast_node, RegisterContactProperty):
+            p = ast_node.property()
+            ptr = p.name()
+            d_ptr = f"d_{ptr}" if self.target.is_gpu() and p.device_flag else "nullptr"
+            tkw = Types.c_keyword(p.type())
+            ptype = "Prop_Integer"  if p.type() == Types.Int32 else \
+                    "Prop_Float"    if p.type() == Types.Double else \
+                    "Prop_Vector"   if p.type() == Types.Vector else \
+                    "Prop_Invalid"
+
+            assert ptype != "Prop_Invalid", "Invalid property type!"
+
+            playout = "AoS" if p.layout() == Layouts.AoS else \
+                      "SoA" if p.layout() == Layouts.SoA else \
+                      "Invalid"
+
+            sizes = ", ".join([str(self.generate_expression(BinOp.inline(size))) for size in ast_node.sizes()])
+
+            if self.target.is_gpu() and p.device_flag:
+                self.print(f"{tkw} *{ptr}, *{d_ptr};")
+                d_ptr = f"&{d_ptr}"
+            else:
+                self.print(f"{tkw} *{ptr};")
+
+            self.print(f"pairs->addContactProperty({p.id()}, \"{p.name()}\", &{ptr}, {d_ptr}, {ptype}, {playout}, {sizes});")
+
         if isinstance(ast_node, RegisterFeatureProperty):
             fp = ast_node.feature_property()
             ptr = fp.name()
@@ -601,7 +647,8 @@ class CGen:
                 lhs = self.generate_expression(ast_node.lhs, mem, index)
                 rhs = self.generate_expression(ast_node.rhs, index=index)
                 operator = ast_node.operator()
-                return f"({lhs} {operator.symbol()} {rhs})"
+                return f"({operator.symbol()}({lhs}))" if operator.is_unary() else \
+                       f"({lhs} {operator.symbol()} {rhs})"
 
             if ast_node.is_vector_kind():
                 assert index is not None, "Index must be set for vector reference!"
@@ -660,6 +707,17 @@ class CGen:
         if isinstance(ast_node, Property):
             return ast_node.name()
 
+        if isinstance(ast_node, ContactPropertyAccess):
+            assert not ast_node.is_vector_kind() or index is not None, "Index must be set for vector property access!"
+            prop_name = self.generate_expression(ast_node.contact_prop)
+
+            if mem or ast_node.inlined is True:
+                index_expr = ast_node.index if not ast_node.is_vector_kind() else ast_node.get_index_expression(index)
+                index_g = self.generate_expression(index_expr)
+                return f"{prop_name}[{index_g}]"
+
+            return f"cp{ast_node.id()}" + (f"_{index}" if ast_node.is_vector_kind() else "")
+
         if isinstance(ast_node, PropertyAccess):
             assert not ast_node.is_vector_kind() or index is not None, "Index must be set for vector property access!"
             prop_name = self.generate_expression(ast_node.prop)
diff --git a/src/pairs/ir/bin_op.py b/src/pairs/ir/bin_op.py
index f098e0b6c4cf9d86d98d79578c50e30ecedc1aae..66673b285f0c45588dd92e2c63bcb338e3134eda 100644
--- a/src/pairs/ir/bin_op.py
+++ b/src/pairs/ir/bin_op.py
@@ -83,6 +83,10 @@ class BinOp(VectorExpression):
 
     def infer_type(lhs, rhs, op):
         lhs_type = lhs.type()
+
+        if op.is_unary():
+            return lhs_type
+
         rhs_type = rhs.type()
 
         if op.is_conditional():
@@ -235,7 +239,7 @@ class ASTTerm(ASTNode):
         return BinOp(self.sim, self, other, Operators.BinXor)
 
     def __invert__(self):
-        return BinOp(self.sim, self, None, Operators.BinNeg)
+        return BinOp(self.sim, self, None, Operators.Invert)
 
     def and_op(self, other):
         return BinOp(self.sim, self, other, Operators.And)
diff --git a/src/pairs/ir/const_vector.py b/src/pairs/ir/const_vector.py
new file mode 100644
index 0000000000000000000000000000000000000000..adf94031a7a861b874b072bf0526ac223fb4746e
--- /dev/null
+++ b/src/pairs/ir/const_vector.py
@@ -0,0 +1,44 @@
+from pairs.ir.ast_node import ASTNode
+from pairs.ir.lit import Lit
+from pairs.ir.types import Types
+
+
+class ConstVector(ASTNode):
+    def __init__(self, sim, values):
+        assert isinstance(values, list) and len(values) == sim.ndims(), \
+            "ConstVector(): Given list is invalid!"
+        super().__init__(sim)
+        self.values = values
+        self.array = None
+
+    def __str__(self):
+        return f"ConstVector<{self.values}>"
+
+    def type(self):
+        return Types.Vector
+
+    def __getitem__(self, expr_ast):
+        if isinstance(expr_ast, Lit) or isinstance(expr_ast, int):
+            return self.values[expr_ast]
+
+        if self.array is None:
+            self.array = self.sim.add_static_array(f"cv{self.const_vector_id}", self.sim.ndims(), Types.Double)
+
+        return self_array[expr_ast]
+
+
+class ZeroVector(ASTNode):
+    def __init__(self, sim):
+        super().__init__(sim)
+
+    def __str__(self):
+        return f"ZeroVector<{self.values}>"
+
+    def type(self):
+        return Types.Vector
+
+    def __getitem__(self, expr_ast):
+        # This allows out-of-bound access to the vector, which may not be good.
+        # It is however difficult to evaluate this possibilty without performance
+        # loss when the expression value is only known at runtime
+        return Lit(self.sim, 0.0)
diff --git a/src/pairs/ir/loops.py b/src/pairs/ir/loops.py
index 553b5ad689ddeefd17e8aa376c189a8ea5d3229c..75a352ac0585259b26f86574322e21d435c11032 100644
--- a/src/pairs/ir/loops.py
+++ b/src/pairs/ir/loops.py
@@ -97,3 +97,8 @@ class While(ASTNode):
 
     def children(self):
         return [self.cond, self.block]
+
+
+class Continue(ASTNode):
+    def __init__(self, sim):
+        super().__init__(sim)
diff --git a/src/pairs/ir/math.py b/src/pairs/ir/math.py
index 3f74aea6491b5ae278d51e7df35d9607dae3881a..ed770f2ec370a18719e3d31b05f12bc0afb76728 100644
--- a/src/pairs/ir/math.py
+++ b/src/pairs/ir/math.py
@@ -3,7 +3,7 @@ from pairs.ir.types import Types
 
 
 class Sqrt(ASTTerm):
-    def __init__(self, sim, expr, cast_type):
+    def __init__(self, sim, expr):
         super().__init__(sim)
         self.expr = expr
 
diff --git a/src/pairs/ir/operators.py b/src/pairs/ir/operators.py
index 8e094e1cdaf65531744e6c643c258d9043637c55..e3bdf8ed6dcd77e50cc79dfdab57a77efe14ec17 100644
--- a/src/pairs/ir/operators.py
+++ b/src/pairs/ir/operators.py
@@ -28,7 +28,10 @@ class Operators:
     BinAnd  =   Op('&')
     BinOr   =   Op('|')
     BinXor  =   Op('^')
-    BinNeg  =   Op('~',  unary=True)
+    UAdd    =   Op('+',  unary=True)
+    USub    =   Op('-',  unary=True)
+    Invert  =   Op('~',  unary=True)
+    Not     =   Op('!',  unary=True)
     Eq      =   Op('==', conditional=True)
     Neq     =   Op('!=', conditional=True)
     Lt      =   Op('<',  conditional=True)
diff --git a/src/pairs/ir/properties.py b/src/pairs/ir/properties.py
index 8a63be2ad153269a029353a46ba3a99b14797e90..caa99e17d77dc1894f42c73a60396429d3fc7053 100644
--- a/src/pairs/ir/properties.py
+++ b/src/pairs/ir/properties.py
@@ -188,3 +188,157 @@ class ReallocProperty(ASTNode):
 
     def __str__(self):
         return f"ReallocProperty<{self.prop.name()}>"
+
+
+class ContactProperties:
+    def __init__(self, sim):
+        self.sim = sim
+        self.contact_properties = []
+
+    def add(self, cp_name, cp_type, cp_layout, cp_default):
+        cp = ContactProperty(self.sim, cp_name, cp_type, cp_layout, cp_default)
+        self.contact_properties.append(cp)
+        return cp
+
+    def find(self, cp_name):
+        contact_prop = [cp for cp in self.contact_properties if cp.name() == cp_name]
+        if contact_prop:
+            return contact_prop[0]
+
+        return None
+
+    def __iter__(self):
+        yield from self.contact_properties
+
+
+class ContactProperty(ASTTerm):
+    last_contact_prop_id = 0
+
+    def __init__(self, sim, name, dtype, layout, default):
+        super().__init__(sim)
+        self.contact_prop_id = ContactProperty.last_contact_prop_id
+        self.contact_prop_name = name
+        self.contact_prop_type = dtype
+        self.contact_prop_layout = layout
+        self.contact_prop_default = default
+        self.device_flag = False
+        ContactProperty.last_contact_prop_id += 1
+
+    def __str__(self):
+        return f"ContactProperty<{self.contact_prop_name}>"
+
+    def id(self):
+        return self.contact_prop_id
+
+    def name(self):
+        return self.contact_prop_name
+
+    def type(self):
+        return self.contact_prop_type
+
+    def layout(self):
+        return self.contact_prop_layout
+
+    def default(self):
+        return self.contact_prop_default
+
+    def ndims(self):
+        return 1 if self.contact_prop_type != Types.Vector else 2
+
+    def sizes(self):
+        neighbor_list_sizes = [self.sim.particle_capacity, self.sim.neighbor_capacity]
+        return neighbor_list_sizes if self.contact_prop_type != Types.Vector \
+               else [self.sim.ndims()] + neighbor_list_sizes
+
+    def __getitem__(self, expr):
+        return ContactPropertyAccess(self.sim, self, expr)
+
+
+class ContactPropertyAccess(ASTTerm, VectorExpression):
+    last_contact_prop_acc = 0
+
+    def new_id():
+        ContactPropertyAccess.last_contact_prop_acc += 1
+        return ContactPropertyAccess.last_contact_prop_acc - 1
+
+    def __init__(self, sim, contact_prop, index):
+        assert isinstance(index, tuple), "Two indexes must be used for contact property access!"
+        super().__init__(sim)
+        self.acc_id = ContactPropertyAccess.new_id()
+        self.contact_prop = contact_prop
+        self.index = Lit.cvt(sim, index)
+        self.inlined = False
+        self.terminals = set()
+
+    def __str__(self):
+        return f"ContactPropertyAccess<{self.contact_prop}, {self.index}>"
+
+    def copy(self):
+        return ContactPropertyAccess(self.sim, self.contact_prop, self.index)
+
+    def vector_index(self, v_index):
+        sizes = self.contact_prop.sizes()
+        layout = self.contact_prop.layout()
+        index = self.index * sizes[0] + v_index if layout == Layouts.AoS else \
+                v_index * sizes[1] + self.index if layout == Layouts.SoA else \
+                None
+
+        assert index is not None, "Invalid data layout"
+        return index
+
+    def inline_rec(self):
+        self.inlined = True
+        return self
+
+    def propagate_through(self):
+        return []
+
+    def id(self):
+        return self.acc_id
+
+    def type(self):
+        return self.contact_prop.type()
+
+    def add_terminal(self, terminal):
+        self.terminals.add(terminal)
+
+    def children(self):
+        return [self.contact_prop, self.index] + list(super().children())
+
+    def __getitem__(self, index):
+        super().__getitem__(index)
+        return VectorAccess(self.sim, self, Lit.cvt(self.sim, index))
+
+
+class RegisterContactProperty(ASTNode):
+    def __init__(self, sim, prop, sizes):
+        super().__init__(sim)
+        self.prop = prop
+        self.sizes_list = [Lit.cvt(sim, s) for s in sizes]
+        self.sim.add_statement(self)
+
+    def property(self):
+        return self.prop
+
+    def sizes(self):
+        return self.sizes_list
+
+    def __str__(self):
+        return f"RegisterContactProperty<{self.prop.name()}>"
+
+
+class ReallocContactProperty(ASTNode):
+    def __init__(self, sim, prop, sizes):
+        super().__init__(sim)
+        self.prop = prop
+        self.sizes_list = [Lit.cvt(sim, s) for s in sizes]
+        self.sim.add_statement(self)
+
+    def property(self):
+        return self.prop
+
+    def sizes(self):
+        return self.sizes_list
+
+    def __str__(self):
+        return f"ReallocContactProperty<{self.prop.name()}>"
diff --git a/src/pairs/ir/select.py b/src/pairs/ir/select.py
index c9a02604c1841bd266612fc5ecf15cd1f7a595cd..16fdb583194699ffc392500d29e3de7842d1c616 100644
--- a/src/pairs/ir/select.py
+++ b/src/pairs/ir/select.py
@@ -9,6 +9,10 @@ class Select(ASTNode):
         self.cond = Lit.cvt(sim, cond)
         self.expr_if = BinOp.inline(Lit.cvt(sim, expr_if))
         self.expr_else = BinOp.inline(Lit.cvt(sim, expr_else))
+        assert self.expr_if.type() == self.expr_else.type(), "Select: expressions must be of same type!"
+
+    def type(self):
+        return self.expr_if.type()
 
     def children(self):
         return [self.cond, self.expr_if, self.expr_else]
diff --git a/src/pairs/ir/symbols.py b/src/pairs/ir/symbols.py
index 556d073fb8d871a3ad024fda7e62b08b1b8a3d38..370136a6b4357a291573928273648fd83708dfd7 100644
--- a/src/pairs/ir/symbols.py
+++ b/src/pairs/ir/symbols.py
@@ -1,6 +1,7 @@
 from pairs.ir.ast_node import ASTNode
 from pairs.ir.assign import Assign
 from pairs.ir.bin_op import ASTTerm 
+from pairs.ir.types import Types
 
 
 class Symbol(ASTTerm):
@@ -10,10 +11,33 @@ class Symbol(ASTTerm):
         self.assign_to = None
 
     def __str__(self):
-        return f"Symbol<{self.var_type}>"
+        return f"Symbol<{self.sym_type}>"
 
     def assign(self, node):
         self.assign_to = node
 
     def type(self):
         return self.sym_type
+
+    def __getitem__(self, index):
+        return SymbolAccess(self.sim, self, index)
+
+
+class SymbolAccess(ASTTerm):
+    def __init__(self, sim, symbol, indexes):
+        super().__init__(sim)
+        self._symbol = symbol
+        self._indexes = indexes if isinstance(indexes, list) else [indexes]
+        assert symbol.type() == Types.Vector, "Only vector symbols can be indexed!"
+
+    def __str__(self):
+        return f"SymbolAccess<{self._symbol, self._indexes}>"
+
+    def symbol(self):
+        return self._symbol
+
+    def indexes(self):
+        return self._indexes
+
+    def type(self):
+        return self._symbol.type()
diff --git a/src/pairs/mapping/funcs.py b/src/pairs/mapping/funcs.py
index f79ad7de08e4de688c8e3cc92c80d01fe5c47fc6..74e23b687dd220905903c76704cb9aaf5db8103c 100644
--- a/src/pairs/mapping/funcs.py
+++ b/src/pairs/mapping/funcs.py
@@ -2,8 +2,15 @@ import ast
 import inspect
 from pairs.ir.assign import Assign
 from pairs.ir.bin_op import BinOp
-from pairs.ir.loops import ParticleFor
+from pairs.ir.block import Block
+from pairs.ir.branches import Filter
+from pairs.ir.const_vector import ZeroVector
+from pairs.ir.lit import Lit
+from pairs.ir.loops import ParticleFor, Continue
+from pairs.ir.math import Sqrt
 from pairs.ir.operators import Operators
+from pairs.ir.select import Select
+from pairs.ir.types import Types
 from pairs.sim.interaction import ParticleInteraction
 
 
@@ -26,8 +33,91 @@ class FetchParticleFuncInfo(ast.NodeVisitor):
         return self._params
 
 
+class Keywords:
+    def __init__(self, sim):
+        self.sim = sim
+
+    def get_method(self, method_name):
+        method = getattr(self, method_name, None)
+        return method if callable(method) else None
+
+    def __call__(self, keyword, args):
+        method = self.get_method(f"keyword_{keyword}")
+        assert method is not None, "Invalid keyword: {keyword}"
+        return method(args)
+
+    def exists(self, keyword):
+        method = self.get_method(f"keyword_{keyword}")
+        return method is not None
+
+    def keyword_select(self, args):
+        assert len(args) == 3, "select() keyword requires three parameter!"
+        return Select(self.sim, args[0], args[1], args[2])
+
+    def keyword_skip_when(self, args):
+        assert len(args) == 1, "skip_when() keyword requires one parameter!"
+        return Filter(self.sim, args[0], Block(self.sim, [Continue(self.sim)]))
+
+    def keyword_min(self, args):
+        e_min = args[0]
+        for a in args[1:]:
+            e_min = Min(self.sim, e_min, a)
+
+        return e_min
+
+    def keyword_max(self, args):
+        e_max = args[0]
+        for a in args[1:]:
+            e_max = Max(self.sim, e_max, a)
+
+        return e_max
+
+    def keyword_length(self, args):
+        assert len(args) == 1, "length() keyword requires one parameter!"
+        vector = args[0]
+        assert vector.type() == Types.Vector, "length(): Argument must be a vector!"
+        return Sqrt(self.sim, sum([vector[d] * vector[d] for d in range(self.sim.ndims())]))
+
+    def keyword_dot(self, args):
+        assert len(args) == 2, "dot() keyword requires two parameters!"
+        vector1 = args[0]
+        vector2 = args[1]
+        assert vector1.type() == Types.Vector, "dot(): First argument must be a vector!"
+        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_normalized(self, args):
+        assert len(args) == 1, "normalized() keyword requires one parameter!"
+        vector = args[0]
+        assert vector.type() == Types.Vector, "length(): Argument must be a vector!"
+        length = self.keyword_length([vector])
+        inv_length = Lit(self.sim, 1.0) / length
+        return Select(self.sim, length > Lit(self.sim, 0.0), vector * inv_length, ZeroVector(self.sim))
+
+    def keyword_square_length(self, args):
+        assert len(args) == 1, "length() keyword requires one parameter!"
+        vector = args[0]
+        assert vector.type() == Types.Vector, "length(): Argument must be a vector!"
+        return sum([vector[d] * vector[d] for d in range(self.sim.ndims())])
+
+
 class BuildParticleIR(ast.NodeVisitor):
-    def get_op(op):
+    def get_unary_op(op):
+        if isinstance(op, ast.UAdd):
+            return Operators.UAdd
+
+        if isinstance(op, ast.USub):
+            return Operators.USub
+
+        if isinstance(op, ast.Invert):
+            return Operators.Invert
+
+        if isinstance(op, ast.Not):
+            return Operators.Not
+
+        raise Exception("Invalid operator: {}".format(ast.dump(op)))
+
+    def get_binary_op(op):
         if isinstance(op, ast.Add):
             return Operators.Add
 
@@ -49,11 +139,13 @@ class BuildParticleIR(ast.NodeVisitor):
         self.sim = sim
         self.ctx_symbols = ctx_symbols
         self.ctx_calls = ctx_calls
+        self.keywords = Keywords(sim)
 
     def add_symbols(self, symbols):
         self.ctx_symbols.update(symbols)
 
     def visit_Assign(self, node):
+        print(ast.dump(node))
         assert len(node.targets) == 1, "Only one target is allowed on assignments!"
         lhs = self.visit(node.targets[0])
         rhs = self.visit(node.value)
@@ -78,12 +170,15 @@ class BuildParticleIR(ast.NodeVisitor):
         assert not isinstance(lhs, UndefinedSymbol), f"Undefined lhs used in BinOp: {lhs.symbol_id}"
         rhs = self.visit(node.right)
         assert not isinstance(rhs, UndefinedSymbol), f"Undefined rhs used in BinOp: {rhs.symbol_id}"
-        return BinOp(self.sim, lhs, rhs, BuildParticleIR.get_op(node.op))
+        return BinOp(self.sim, lhs, rhs, BuildParticleIR.get_binary_op(node.op))
 
     def visit_Call(self, node):
-        func = self.visit(node.func)
+        func = self.visit(node.func).symbol_id
         args = [self.visit(a) for a in node.args]
 
+        if self.keywords.exists(func):
+            return self.keywords(func, args)
+
         for c in self.ctx_calls:
             if c['func'] == func and len(c['args']) == len(args) and all([c['args'][a] == args[a] for a in range(0, len(args))]):
                 return c['value']
@@ -112,6 +207,10 @@ class BuildParticleIR(ast.NodeVisitor):
         if as_feature_prop is not None:
             return as_feature_prop
 
+        as_contact_prop = self.sim.contact_property(node.id)
+        if as_contact_prop is not None:
+            return as_contact_prop
+
         as_var = self.sim.var(node.id)
         if as_var is not None:
             return as_var
@@ -129,6 +228,14 @@ class BuildParticleIR(ast.NodeVisitor):
         #print(ast.dump(node))
         return tuple(self.visit(v) for v in node.elts)
 
+    def visit_UnaryOp(self, node):
+        #print(ast.dump(node))
+        operand = self.visit(node.operand)
+        assert not isinstance(operand, UndefinedSymbol), \
+            f"Undefined operand used in UnaryOp: {operand.symbol_id}"
+        return BinOp(self.sim, operand, None, BuildParticleIR.get_unary_op(node.op))
+
+
 def compute(sim, func, cutoff_radius=None, symbols={}):
     src = inspect.getsource(func)
     tree = ast.parse(src, mode='exec')
diff --git a/src/pairs/sim/contact_history.py b/src/pairs/sim/contact_history.py
new file mode 100644
index 0000000000000000000000000000000000000000..55e4b0fd7109024ae14eed226a8c5bfa63071671
--- /dev/null
+++ b/src/pairs/sim/contact_history.py
@@ -0,0 +1,67 @@
+from pairs.ir.block import pairs_device_block
+from pairs.ir.branches import Branch, Filter
+from pairs.ir.loops import ParticleFor
+from pairs.ir.types import Types
+from pairs.ir.utils import Print
+from pairs.sim.interaction import ParticleInteraction
+from pairs.sim.lowerable import Lowerable
+
+
+class ContactHistory:
+    def __init__(self, cell_lists):
+        self.sim = cell_lists.sim
+        self.cell_lists = cell_lists
+        self.contact_lists = self.sim.add_array('contact_lists', [self.sim.particle_capacity, self.sim.neighbor_capacity], Types.Int32)
+        self.num_contacts = self.sim.add_array('num_contacts', self.sim.particle_capacity, Types.Int32)
+
+
+class BuildContactHistory(Lowerable):
+    def __init__(self, sim, contact_history):
+        super().__init__(sim)
+        self.contact_history = contact_history
+
+    @pairs_device_block
+    def lower(self):
+        contact_history = self.contact_history
+        cell_lists = self.contact_history.cell_lists
+        neighbor_lists = self.contact_history.neighbor_lists
+        contact_lists = self.contact_history.contact_lists
+        num_contacts = self.contact_history.num_contacts
+        sim.module_name("build_contact_history")
+        last_contact_id = self.sim.add_temp_var(0)
+        neighbor_contact = self.sim.add_temp_var(0)
+
+        for i in ParticleFor(self.sim):
+            last_contact_id.set(0)
+            for j in NeighborFor(self.sim, i, cell_lists, neighbor_lists):
+                neighbor_contact.set(-1)
+                for k in For(self.sim, 0, num_contacts[i]):
+                    if contact_lists[i][k] == j:
+                        neighbor_contact.set(k)
+
+                for _ in Filter(self.sim, neighbor_contact >= 0):
+                    contact_lists[i][k].set(contact_lists[i][last_contact_id])
+
+                    for contact_prop in self.sim.contact_properties():
+                        tmp = self.sim.add_temp_var(contact_prop[i][last_contact_id])
+                        contact_prop[i][last_contact_id].set(contact_prop[i][k])
+                        contact_prop[i][k].set(tmp)
+
+                    contact_lists[i][last_contact_id].set(j)
+
+                last_contact_id.set(last_contact_id + 1)
+
+            last_contact_id.set(0)
+            for j in NeighborFor(self.sim, i, cell_lists, neighbor_lists):
+                neighbor_contact.set(-1)
+                for k in For(self.sim, 0, num_contacts[i]):
+                    if contact_lists[i][k] == j:
+                        neighbor_contact.set(k)
+
+                for _ in Filter(self.sim, neighbor_contact < 0):
+                    for contact_prop in self.sim.contact_properties():
+                        contact_prop[i][last_contact_id].set(contact_prop.default())
+
+                    contact_lists[i][last_contact_id].set(j)
+
+                last_contact_id.set(last_contact_id + 1)
diff --git a/src/pairs/sim/neighbor_lists.py b/src/pairs/sim/neighbor_lists.py
index eda0018e18c104e4d374d1a5b3b2f23fff5af94f..8042e778904a379e89925053fac924f7adbb3823 100644
--- a/src/pairs/sim/neighbor_lists.py
+++ b/src/pairs/sim/neighbor_lists.py
@@ -11,8 +11,7 @@ class NeighborLists:
     def __init__(self, cell_lists):
         self.sim = cell_lists.sim
         self.cell_lists = cell_lists
-        self.capacity = self.sim.add_var('neighborlist_capacity', Types.Int32, 32)
-        self.neighborlists = self.sim.add_array('neighborlists', [self.sim.particle_capacity, self.capacity], Types.Int32)
+        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)
 
 
diff --git a/src/pairs/sim/properties.py b/src/pairs/sim/properties.py
index a0b40571b164b06d3df3c6b454d6f648ef2d84bb..37b28d26854b356dbe0056c40a9711f0e9613009 100644
--- a/src/pairs/sim/properties.py
+++ b/src/pairs/sim/properties.py
@@ -1,7 +1,7 @@
 from pairs.ir.block import pairs_device_block, pairs_inline
 from pairs.ir.loops import ParticleFor
 from pairs.ir.memory import Malloc, Realloc
-from pairs.ir.properties import RegisterProperty
+from pairs.ir.properties import RegisterProperty, RegisterContactProperty
 from pairs.ir.types import Types
 from pairs.ir.utils import Print
 from pairs.sim.lowerable import Lowerable, FinalLowerable
@@ -32,6 +32,25 @@ class PropertiesAlloc(FinalLowerable):
                 RegisterProperty(self.sim, p, sizes)
 
 
+class ContactPropertiesAlloc(FinalLowerable):
+    def __init__(self, sim, realloc=False):
+        self.sim = sim
+
+    @pairs_inline
+    def lower(self):
+        capacity = sum(self.sim.contact_properties.capacities)
+        for p in self.sim.contact_properties.all():
+            sizes = []
+            if Types.is_real(p.type()) or Types.is_integer(p.type()):
+                sizes = [capacity]
+            elif p.type() == Types.Vector:
+                sizes = [capacity, self.sim.ndims()]
+            else:
+                raise Exception("Invalid contact property type!")
+
+            RegisterContactProperty(self.sim, p, sizes)
+
+
 class PropertiesResetVolatile(Lowerable):
     def __init__(self, sim):
         self.sim = sim
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index b35b0b19cd1b63c46b59d565da1d1e2afa12f59e..62e3a9db9ac8fb61f3db0f62f453693079a6f23c 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -6,7 +6,7 @@ from pairs.ir.functions import Call_Void
 from pairs.ir.kernel import Kernel
 from pairs.ir.layouts import Layouts
 from pairs.ir.module import Module
-from pairs.ir.properties import Properties
+from pairs.ir.properties import Properties, ContactProperties
 from pairs.ir.symbols import Symbol
 from pairs.ir.types import Types
 from pairs.ir.variables import Variables
@@ -39,7 +39,9 @@ class Simulation:
         self.arrays = Arrays(self)
         self.features = Features(self)
         self.feature_properties = FeatureProperties(self)
+        self.contact_properties = ContactProperties(self)
         self.particle_capacity = self.add_var('particle_capacity', Types.Int32, particle_capacity)
+        self.neighbor_capacity = self.add_var('neighbor_capacity', Types.Int32, particle_capacity)
         self.nlocal = self.add_var('nlocal', Types.Int32)
         self.nghost = self.add_var('nghost', Types.Int32)
         self.resizes = self.add_array('resizes', 3, Types.Int32, arr_sync=False)
@@ -120,6 +122,10 @@ class Simulation:
         assert self.property(prop_name) is None, f"Property already defined: {prop_name}"
         return self.feature_properties.add(feature, prop_name, prop_type, prop_data)
 
+    def add_contact_property(self, prop_name, prop_type, prop_default, layout=Layouts.AoS):
+        assert self.property(prop_name) is None, f"Property already defined: {prop_name}"
+        return self.contact_properties.add(prop_name, prop_type, layout, prop_default)
+
     def property(self, prop_name):
         return self.properties.find(prop_name)
 
@@ -132,6 +138,9 @@ class Simulation:
     def feature_property(self, feature_prop_name):
         return self.feature_properties.find(feature_prop_name)
 
+    def contact_property(self, contact_prop_name):
+        return self.contact_properties.find(contact_prop_name)
+
     def add_array(self, arr_name, arr_sizes, arr_type, arr_layout=Layouts.AoS, arr_sync=True):
         assert self.array(arr_name) is None, f"Array already defined: {arr_name}"
         return self.arrays.add(arr_name, arr_sizes, arr_type, arr_layout, arr_sync)
@@ -269,12 +278,14 @@ class Simulation:
     def generate(self):
         assert self._target is not None, "Target not specified!"
         comm = Comm(self, self._dom_part)
+        contact_history = ContactHistory(self, self.cell_lists)
 
         timestep = Timestep(self, self.ntimesteps, [
             (comm.exchange(), 20),
             (comm.borders(), comm.synchronize(), 20),
             (CellListsBuild(self, self.cell_lists), 20),
             (NeighborListsBuild(self, self.neighbor_lists), 20),
+            (BuildContactHistory(self, contact_history), 20),
             PropertiesResetVolatile(self),
             self.functions
         ])
@@ -294,6 +305,7 @@ class Simulation:
             VariablesDecl(self),
             ArraysDecl(self),
             PropertiesAlloc(self),
+            ContactPropertiesAlloc(self),
             FeaturePropertiesAlloc(self)
         ])
 
diff --git a/src/pairs/transformations/expressions.py b/src/pairs/transformations/expressions.py
index 15dc68843d63b9d00f3d68b426eea210622e96b0..768ad2768eb8f81e510bc8bbc5a7c28d4e0c7d49 100644
--- a/src/pairs/transformations/expressions.py
+++ b/src/pairs/transformations/expressions.py
@@ -12,6 +12,13 @@ class ReplaceSymbols(Mutator):
     def mutate_Symbol(self, ast_node):
         return ast_node.assign_to
 
+    def mutate_SymbolAccess(self, ast_node):
+        access = ast_node.symbol()
+        for i in ast_node.indexes():
+            access = access[i]
+
+        return access
+
 
 class SimplifyExpressions(Mutator):
     def __init__(self, ast=None):