diff --git a/examples/lj.py b/examples/lj.py
index fef1614371f4f45782cee3f8e6bb086b326ffb91..5d4a5d5010e784702e5520b6951946482f15efb5 100644
--- a/examples/lj.py
+++ b/examples/lj.py
@@ -29,12 +29,12 @@ sigma6 = sigma ** 6
 
 psim = pairs.simulation("lj", 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_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_property('type', 'epsilon', Types.Double, [sigma for i in range(ntypes * ntypes)])
-psim.add_feature_property('type', 'sigma6', Types.Double, [epsilon for i in range(ntypes * 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'])
 psim.build_neighbor_lists(cutoff_radius + skin)
 psim.vtk_output(f"output/test_{target}")
diff --git a/src/pairs/__init__.py b/src/pairs/__init__.py
index d0e878d7d3fc05e7d9e80a22782a53e002a18dcc..2807d54642f7f47fc2be74edab2b30b5a5320089 100644
--- a/src/pairs/__init__.py
+++ b/src/pairs/__init__.py
@@ -1,3 +1,4 @@
+from pairs.ir.types import Types
 from pairs.code_gen.cgen import CGen
 from pairs.code_gen.target import Target
 from pairs.sim.simulation import Simulation
@@ -11,3 +12,12 @@ def target_cpu():
 
 def target_gpu():
     return Target(Target.Backend_CUDA, Target.Feature_GPU)
+
+def int32():
+    return Types.Int32
+
+def double():
+    return Types.Double
+
+def vector():
+    return Types.Vector
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index da6a9c97be46ab0c5633993e2a3426cb09b44882..c33c32cdb5e0e2de3f0cf71d0e815fc109bc4d4f 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -15,7 +15,7 @@ from pairs.ir.loops import For, Iter, ParticleFor, While
 from pairs.ir.math import Ceil, Sqrt
 from pairs.ir.memory import Malloc, Realloc
 from pairs.ir.module import ModuleCall
-from pairs.ir.features import RegisterFeature
+from pairs.ir.features import RegisterFeatureProperty
 from pairs.ir.particle_attributes import ParticleAttributeList
 from pairs.ir.properties import Property, PropertyAccess, RegisterProperty, ReallocProperty
 from pairs.ir.select import Select
diff --git a/src/pairs/ir/features.py b/src/pairs/ir/features.py
index ae50dd156b1c4ff71399572256415a45fb1a4eda..e202ff3e821ffab13e032097bec17bd98baa9dff 100644
--- a/src/pairs/ir/features.py
+++ b/src/pairs/ir/features.py
@@ -12,10 +12,10 @@ class Features:
         self.sim = sim
         self.features = []
 
-    def add(self, f_name):
-        f = Feature(self.sim, f_name)
+    def add(self, f_name, f_nkinds):
+        f = Feature(self.sim, f_name, f_nkinds)
         self.features.append(f)
-        return p
+        return f
 
     def nfeatures(self):
         return len(self.features)
@@ -75,7 +75,7 @@ class FeatureProperties:
         return len(self.features)
 
     def find(self, fp_name):
-        feature_prop = [fp for fp in self.feature_proeprties if fp.name() == fp_name]
+        feature_prop = [fp for fp in self.feature_properties if fp.name() == fp_name]
         if feature_prop:
             return feature_prop[0]
 
diff --git a/src/pairs/mapping/funcs.py b/src/pairs/mapping/funcs.py
index adad3d223da2a78a7e5feac4246464e7b3eb0c70..50ca74ab6e06bf75d6d528180d40e005235914ad 100644
--- a/src/pairs/mapping/funcs.py
+++ b/src/pairs/mapping/funcs.py
@@ -108,6 +108,10 @@ class BuildParticleIR(ast.NodeVisitor):
         if as_prop is not None:
             return as_prop
 
+        as_feature_prop = self.sim.feature_property(node.id)
+        if as_feature_prop is not None:
+            return as_feature_prop
+
         as_var = self.sim.var(node.id)
         if as_var is not None:
             return as_var
diff --git a/src/pairs/sim/features.py b/src/pairs/sim/features.py
index 4de75518a16b564f9f235da8881f9b076fa7578e..70f7d44bbb2c7be0fac8ff63f256cbd9c6219848 100644
--- a/src/pairs/sim/features.py
+++ b/src/pairs/sim/features.py
@@ -1,4 +1,5 @@
-from pairs.ir.features import RegisterFeature
+from pairs.ir.block import pairs_device_block, pairs_inline
+from pairs.ir.features import RegisterFeatureProperty
 from pairs.sim.lowerable import FinalLowerable
 
 
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index 6db3b091bbb7e0c81acb8b8c5257ae2cc2f12e92..da35c410b13bda7d23b13cd718e76456d064916c 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -110,9 +110,9 @@ class Simulation:
         self.position_prop = self.properties.add(prop_name, Types.Vector, value, vol, layout)
         return self.position_prop
 
-    def add_feature(self, feature_name):
+    def add_feature(self, feature_name, nkinds):
         assert self.feature(feature_name) is None, f"Feature already defined: {feature_name}"
-        return self.features.add(feature_name)
+        return self.features.add(feature_name, nkinds)
 
     def add_feature_property(self, feature_name, prop_name, prop_type, prop_data):
         feature = self.feature(feature_name)
@@ -127,7 +127,7 @@ class Simulation:
         return self.position_prop
 
     def feature(self, feature_name):
-        return self.feature.find(feature_name)
+        return self.features.find(feature_name)
 
     def feature_property(self, feature_prop_name):
         return self.feature_properties.find(feature_prop_name)