diff --git a/examples/lj.py b/examples/lj.py
index e703122d5469e96458dd831aedbfa4cf9c456090..fef1614371f4f45782cee3f8e6bb086b326ffb91 100644
--- a/examples/lj.py
+++ b/examples/lj.py
@@ -4,8 +4,8 @@ import sys
 
 def lj(i, j):
     sr2 = 1.0 / rsq
-    sr6 = sr2 * sr2 * sr2 * sigma6
-    force[i] += delta * 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon
+    sr6 = sr2 * sr2 * sr2 * sigma6[i, j]
+    force[i] += delta * 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon[i, j]
 
 
 def euler(i):
@@ -22,19 +22,23 @@ if target != 'cpu' and target != 'gpu':
 dt = 0.005
 cutoff_radius = 2.5
 skin = 0.3
+ntypes = 4
 sigma = 1.0
 epsilon = 1.0
 sigma6 = sigma ** 6
 
 psim = pairs.simulation("lj", debug=True)
-psim.add_real_property('mass', 1.0)
 psim.add_position('position')
-psim.add_vector_property('velocity')
-psim.add_vector_property('force', vol=True)
-psim.from_file("data/minimd_setup_32x32x32.input", ['mass', 'position', 'velocity'])
+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', 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.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}")
-psim.compute(lj, cutoff_radius, {'sigma6': sigma6, 'epsilon': epsilon})
+psim.compute(lj, cutoff_radius)
 psim.compute(euler, symbols={'dt': dt})
 
 if target == 'gpu':
diff --git a/examples/lj_onetype.py b/examples/lj_onetype.py
new file mode 100644
index 0000000000000000000000000000000000000000..e703122d5469e96458dd831aedbfa4cf9c456090
--- /dev/null
+++ b/examples/lj_onetype.py
@@ -0,0 +1,45 @@
+import pairs
+import sys
+
+
+def lj(i, j):
+    sr2 = 1.0 / rsq
+    sr6 = sr2 * sr2 * sr2 * sigma6
+    force[i] += delta * 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon
+
+
+def euler(i):
+    velocity[i] += dt * force[i] / mass[i]
+    position[i] += dt * velocity[i]
+
+
+cmd = sys.argv[0]
+target = sys.argv[1] if len(sys.argv[1]) > 1 else "none"
+if target != 'cpu' and target != 'gpu':
+    print(f"Invalid target, use {cmd} <cpu/gpu>")
+
+
+dt = 0.005
+cutoff_radius = 2.5
+skin = 0.3
+sigma = 1.0
+epsilon = 1.0
+sigma6 = sigma ** 6
+
+psim = pairs.simulation("lj", debug=True)
+psim.add_real_property('mass', 1.0)
+psim.add_position('position')
+psim.add_vector_property('velocity')
+psim.add_vector_property('force', vol=True)
+psim.from_file("data/minimd_setup_32x32x32.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(euler, symbols={'dt': dt})
+
+if target == 'gpu':
+    psim.target(pairs.target_gpu())
+else:
+    psim.target(pairs.target_cpu())
+
+psim.generate()
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index e4b1fde9f494c486188194e5b0d7293e39205cb4..da6a9c97be46ab0c5633993e2a3426cb09b44882 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -15,7 +15,9 @@ 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.properties import Property, PropertyAccess, PropertyList, RegisterProperty, ReallocProperty
+from pairs.ir.features import RegisterFeature
+from pairs.ir.particle_attributes import ParticleAttributeList
+from pairs.ir.properties import Property, PropertyAccess, RegisterProperty, ReallocProperty
 from pairs.ir.select import Select
 from pairs.ir.sizeof import Sizeof
 from pairs.ir.types import Types
@@ -75,6 +77,13 @@ class CGen:
                     size = self.generate_expression(BinOp.inline(array.alloc_size()))
                     self.print(f"__constant__ {tkw} d_{array.name()}[{size}];")
 
+            for feature_prop in self.sim.feature_properties():
+                if feature_prop.device_flag:
+                    t = feature_prop.type()
+                    tkw = Types.c_keyword(t)
+                    size = self.generate_expression(BinOp.inline(feature_prop.array_size()))
+                    self.print(f"__constant__ {tkw} d_{feature_prop.name()}[{size}];")
+
         self.print("")
 
         for kernel in self.sim.kernels():
@@ -236,6 +245,14 @@ 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, FeatureAccess):
+                feature_access = ast_node.elem
+                feature_name = self.generate_expression(feature_access.feature.name())
+                tkw = Types.c_keyword(feature_access.type())
+                acc_index = self.generate_expression(feature_access.index)
+                acc_ref = f"f{array_access.id()}"
+                self.print(f"const {tkw} {acc_ref} = {feature_name}[{acc_index}];")
+
             if isinstance(ast_node.elem, PropertyAccess):
                 prop_access = ast_node.elem
                 prop_name = self.generate_expression(prop_access.prop)
@@ -482,6 +499,16 @@ class CGen:
 
             self.print(f"pairs->addProperty({p.id()}, \"{p.name()}\", &{ptr}, {d_ptr}, {ptype}, {playout}, {sizes});")
 
+        if isinstance(ast_node, RegisterFeatureProperty):
+            fp = ast_node.feature_property()
+            ptr = fp.name()
+            d_ptr = f"d_{ptr}" if self.target.is_gpu() and p.device_flag else "nullptr"
+            array_size = fp.array_size()
+            nkinds = fp.feature().nkinds()
+            self.print(f"{tkw} {ptr}[{array_size}];")
+            self.print(f"pairs->addFeatureProperty({fp.id()}, \"{fp.name()}\", &{ptr}, {d_ptr}, {nkinds}, {array_size});")
+
+
         if isinstance(ast_node, Timestep):
             self.generate_statement(ast_node.block)
 
@@ -599,10 +626,18 @@ class CGen:
 
             return f"p{ast_node.id()}" + (f"_{index}" if ast_node.is_vector_kind() else "")
 
-        if isinstance(ast_node, PropertyList):
+        if isinstance(ast_node, FeatureAccess):
+            feature_name = self.generate_expression(ast_node.feature.name())
+            if mem or ast_node.inlined is True:
+                index = self.generate_expression(ast_node.index)
+                return f"{feature_name}[{index}]"
+
+            return f"f{ast_node.id()}"
+
+        if isinstance(ast_node, ParticleAttributeList):
             tid = CGen.temp_id
-            list_ref = f"prop_list_{tid}"
-            list_def = ", ".join(str(p.id()) for p in ast_node)
+            list_ref = f"attr_list_{tid}"
+            list_def = ", ".join(str([a.id() for p in ast_node]))
             self.print(f"const int {list_ref}[] = {{{list_def}}};")
             CGen.temp_id += 1
             return list_ref
diff --git a/src/pairs/ir/features.py b/src/pairs/ir/features.py
index 5434cff3adf0879324a49cbe5944f2204f828d44..ae50dd156b1c4ff71399572256415a45fb1a4eda 100644
--- a/src/pairs/ir/features.py
+++ b/src/pairs/ir/features.py
@@ -11,7 +11,6 @@ class Features:
     def __init__(self, sim):
         self.sim = sim
         self.features = []
-        self.feature_properties = []
 
     def add(self, f_name):
         f = Feature(self.sim, f_name)
@@ -35,11 +34,12 @@ class Features:
 class Feature(ASTNode):
     last_feature_id = 0
 
-    def __init__(self, sim, name):
+    def __init__(self, sim, name, nkinds):
         super().__init__(sim)
         self.feature_id = Feature.last_feature_id
         self.feature_name = name
-        self.feature_count = self.sim.add_var(f"count_{self.feature_name}", Types.Int32)
+        self.feature_prop = self.sim.add_property(self.feature_name, Types.Int32)
+        self.feature_nkinds = nkinds
         Feature.last_feature_id += 1
 
     def __str__(self):
@@ -51,8 +51,143 @@ class Feature(ASTNode):
     def name(self):
         return self.feature_name
 
-    def count(self):
-        return self.feature_count
+    def prop(self):
+        return self.feature_prop
+
+    def nkinds(self):
+        return self.feature_nkinds
+
+    def __getitem__(self, expr):
+        return FeatureAccess(self.sim, self, expr)
+
+
+class FeatureProperties:
+    def __init__(self, sim):
+        self.sim = sim
+        self.feature_properties = []
+
+    def add(self, fp_feature, fp_name, fp_type, fp_data, fp_layout=Layouts.AoS):
+        fp = FeatureProperty(self.sim, fp_feature, fp_name, fp_type, fp_data, fp_layout)
+        self.feature_properties.append(fp)
+        return fp
+
+    def nfeatures(self):
+        return len(self.features)
+
+    def find(self, fp_name):
+        feature_prop = [fp for fp in self.feature_proeprties if fp.name() == fp_name]
+        if feature_prop:
+            return feature_prop[0]
+
+        return None
+
+    def __iter__(self):
+        yield from self.feature_properties
+
+
+class FeatureProperty(ASTTerm):
+    last_feature_prop_id = 0
+
+    def __init__(self, sim, feature, name, dtype, data, layout=Layouts.AoS):
+        super().__init__(sim)
+        self.feature_prop_id = FeatureProperty.last_feature_prop_id
+        self.feature_prop_feature = feature
+        self.feature_prop_name = name
+        self.feature_prop_type = dtype
+        self.feature_prop_data = data
+        self.feature_prop_layout = layout
+        self.device_flag = False
+        FeatureProperty.last_feature_prop_id += 1
+
+    def __str__(self):
+        return f"FeatureProperty<{self.feature_prop_name}>"
+
+    def id(self):
+        return self.feature_prop_id
+
+    def feature(self):
+        return self.feature_prop_feature
+
+    def name(self):
+        return self.feature_prop_name
+
+    def type(self):
+        return self.feature_prop_type
+
+    def layout(self):
+        return self.feature_prop_layout
+
+    def ndims(self):
+        return 1 if self.feature_prop_type != Types.Vector else 2
+
+    def sizes(self):
+        return [self.feature_prop_feature.nkinds()] if self.feature_prop_type != Types.Vector \
+               else [self.sim.ndims(), self.feature_prop_feature.nkinds()]
+
+    def array_size(self):
+        nelems = self.feature_prop.feature.nkinds() * \
+                 (1 if self.feature_prop_type != Types.Vector else self.sim.ndims())
+        return nelems * nelems
+
+    def __getitem__(self, expr):
+        return FeaturePropertyAccess(self.sim, self, expr)
+
+
+class FeaturePropertyAccess(ASTTerm, VectorExpression):
+    last_feature_prop_acc = 0
+
+    def new_id():
+        PropertyAccess.last_feature_prop_acc += 1
+        return PropertyAccess.last_feature_prop_acc - 1
+
+    def __init__(self, sim, feature_prop, index):
+        assert isinstance(index, tuple), "Two indexes must be used for feature property access!"
+        super().__init__(sim)
+        self.acc_id = FeaturePropertyAccess.new_id()
+        self.feature_prop = feature_prop
+        feature = self.feature_prop.feature()
+        self.index = Lit.cvt(sim, feature[index[0]] * feature.nkinds() + feature[index[1]])
+        self.inlined = False
+        self.terminals = set()
+
+    def __str__(self):
+        return f"FeaturePropertyAccess<{self.feature_prop}, {self.index}>"
+
+    def copy(self):
+        return FeaturePropertyAccess(self.sim, self.feature_prop, self.index)
+
+    def vector_index(self, v_index):
+        sizes = self.prop.sizes()
+        layout = self.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.feature_prop.type()
+
+    def add_terminal(self, terminal):
+        self.terminals.add(terminal)
+
+    def children(self):
+        return [self.feature_prop, self.index] + list(super().children())
+
+    def __getitem__(self, index):
+        super().__getitem__(index)
+        return VectorAccess(self.sim, self, Lit.cvt(self.sim, index))
 
 
 class FeatureAccess(ASTTerm):
@@ -96,10 +231,23 @@ class FeatureAccess(ASTTerm):
         return self.acc_id
 
     def type(self):
-        return self.prop.type()
+        return Types.Int32
 
     def add_terminal(self, terminal):
         self.terminals.add(terminal)
 
     def children(self):
         return [self.prop, self.index]
+
+
+class RegisterFeatureProperty(ASTNode):
+    def __init__(self, sim, feature_prop):
+        super().__init__(sim)
+        self.feature_prop = feature_prop
+        self.sim.add_statement(self)
+
+    def feature_property(self):
+        return self.feature_prop
+
+    def __str__(self):
+        return f"RegisterFeatureProperty<{self.feature_prop.name()}>"
diff --git a/src/pairs/ir/particle_attributes.py b/src/pairs/ir/particle_attributes.py
new file mode 100644
index 0000000000000000000000000000000000000000..0917d4fa831724196332b4f45880882d03bda96c
--- /dev/null
+++ b/src/pairs/ir/particle_attributes.py
@@ -0,0 +1,30 @@
+from pairs.ir.ast_node import ASTNode
+from pairs.ir.properties import Property
+from pairs.ir.features import Feature
+
+
+class ParticleAttributeList(ASTNode):
+    def __init__(self, sim, items):
+        super().__init__(sim)
+        self.list = []
+
+        for i in items:
+            if isinstance(i, Property):
+                self.list.append(i)
+
+            if isinstance(i, Feature):
+                self.list.append(i.prop())
+
+            if isinstance(i, str):
+                prop = sim.property(i)
+                if prop is None:
+                    prop = sim.feature(i).prop()
+
+                assert prop is not None, f"Attribute not found: {i}"
+                self.list.append(prop)
+
+    def __iter__(self):
+        yield from self.list
+
+    def length(self):
+        return len(self.list)
diff --git a/src/pairs/ir/properties.py b/src/pairs/ir/properties.py
index e43d38642fc69bcf2efffe6dc2a2e7ac702f72ae..8a63be2ad153269a029353a46ba3a99b14797e90 100644
--- a/src/pairs/ir/properties.py
+++ b/src/pairs/ir/properties.py
@@ -14,8 +14,8 @@ class Properties:
         self.capacities = []
         self.defs = {}
 
-    def add(self, p_name, p_type, p_value, p_volatile, p_layout=Layouts.AoS, p_feature=None):
-        p = Property(self.sim, p_name, p_type, p_value, p_volatile, p_layout, p_feature)
+    def add(self, p_name, p_type, p_value, p_volatile, p_layout=Layouts.AoS):
+        p = Property(self.sim, p_name, p_type, p_value, p_volatile, p_layout)
         self.props.append(p)
         self.defs[p_name] = p_value
         return p
@@ -52,13 +52,12 @@ class Properties:
 class Property(ASTNode):
     last_prop_id = 0
 
-    def __init__(self, sim, name, dtype, default, volatile, layout=Layouts.AoS, feature=None):
+    def __init__(self, sim, name, dtype, default, volatile, layout=Layouts.AoS):
         super().__init__(sim)
         self.prop_id = Property.last_prop_id
         self.prop_name = name
         self.prop_type = dtype
         self.prop_layout = layout
-        self.prop_feature = feature
         self.default_value = default
         self.volatile = volatile
         self.device_flag = False
@@ -79,9 +78,6 @@ class Property(ASTNode):
     def layout(self):
         return self.prop_layout
 
-    def feature(self):
-        return self.prop_feature
-
     def default(self):
         return self.default_value
 
@@ -107,16 +103,7 @@ class PropertyAccess(ASTTerm, VectorExpression):
         super().__init__(sim)
         self.acc_id = PropertyAccess.new_id()
         self.prop = prop
-
-        if prop.feature() == None:
-            assert isinstance(index, int), "Only one index must be used for feature property!"
-            self.index = Lit.cvt(sim, index)
-
-        else:
-            assert isinstance(index, tuple), "Two indexes must be used for feature property!"
-            feature = self.prop.feature()
-            self.index = Lit.cvt(sim, feature[index[0]] * feature.count() + feature[index[1]])
-
+        self.index = Lit.cvt(sim, index)
         self.inlined = False
         self.terminals = set()
 
@@ -169,24 +156,6 @@ class PropertyAccess(ASTTerm, VectorExpression):
         return VectorAccess(self.sim, self, Lit.cvt(self.sim, index))
 
 
-class PropertyList(ASTNode):
-    def __init__(self, sim, properties_list):
-        super().__init__(sim)
-        self.list = []
-        for p in properties_list:
-            if isinstance(p, Property):
-                self.list.append(p)
-
-            if isinstance(p, str):
-                self.list.append(sim.prop(p))
-
-    def __iter__(self):
-        yield from self.list
-
-    def length(self):
-        return len(self.list)
-
-
 class RegisterProperty(ASTNode):
     def __init__(self, sim, prop, sizes):
         super().__init__(sim)
diff --git a/src/pairs/sim/features.py b/src/pairs/sim/features.py
new file mode 100644
index 0000000000000000000000000000000000000000..4de75518a16b564f9f235da8881f9b076fa7578e
--- /dev/null
+++ b/src/pairs/sim/features.py
@@ -0,0 +1,12 @@
+from pairs.ir.features import RegisterFeature
+from pairs.sim.lowerable import FinalLowerable
+
+
+class FeaturePropertiesAlloc(FinalLowerable):
+    def __init__(self, sim):
+        self.sim = sim
+
+    @pairs_inline
+    def lower(self):
+        for fp in self.sim.feature_properties:
+            RegisterFeatureProperty(self.sim, fp)
diff --git a/src/pairs/sim/read_from_file.py b/src/pairs/sim/read_from_file.py
index 1d3d4edd5356cb71032cf5f488ba70e16539ed6e..c447f8bef305d68c07736d3a381fd05a66d92430 100644
--- a/src/pairs/sim/read_from_file.py
+++ b/src/pairs/sim/read_from_file.py
@@ -1,16 +1,16 @@
 from pairs.ir.block import pairs_inline
 from pairs.ir.functions import Call_Int, Call_Void
-from pairs.ir.properties import PropertyList
+from pairs.ir.particle_attributes import ParticleAttributeList
 from pairs.ir.types import Types
 from pairs.sim.grid import MutableGrid
 from pairs.sim.lowerable import Lowerable
 
 
-class ReadFromFile(Lowerable):
-    def __init__(self, sim, filename, props):
+class ReadParticleData(Lowerable):
+    def __init__(self, sim, filename, items):
         super().__init__(sim)
         self.filename = filename
-        self.props = PropertyList(sim, props)
+        self.attrs = ParticleAttributeList(sim, items)
         self.grid = MutableGrid(sim, sim.ndims())
         self.grid_buffer = self.sim.add_static_array("grid_buffer", [self.sim.ndims() * 2], Types.Double)
 
@@ -25,4 +25,16 @@ class ReadFromFile(Lowerable):
         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.props, self.props.length()]))
+        self.sim.nlocal.set(Call_Int(self.sim, "pairs::read_particle_data", [self.filename, self.attrs, self.attrs.length()]))
+
+
+class ReadFeatureData(Lowerable):
+    def __init__(self, sim, filename, feature, items):
+        super().__init__(sim)
+        self.filename = filename
+        self.feature = feature
+        self.attrs = ParticleAttributeList(sim, items)
+
+    @pairs_inline
+    def lower(self):
+        Call_Int(self.sim, "pairs::read_feature_data", [self.filename, self.feature.id(), self.attrs, self.attrs.length()])
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index d6ccee4c9c2f2d14ff7963eccb36cc050937868c..6db3b091bbb7e0c81acb8b8c5257ae2cc2f12e92 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -1,6 +1,7 @@
 from pairs.ir.arrays import Arrays
 from pairs.ir.block import Block
 from pairs.ir.branches import Filter
+from pairs.ir.features import Features, FeatureProperties
 from pairs.ir.functions import Call_Void
 from pairs.ir.kernel import Kernel
 from pairs.ir.layouts import Layouts
@@ -15,12 +16,13 @@ from pairs.sim.arrays import ArraysDecl
 from pairs.sim.cell_lists import CellLists, CellListsBuild, CellListsStencilBuild
 from pairs.sim.comm import Comm
 from pairs.sim.domain_partitioning import DimensionRanges
+from pairs.sim.features import RegisterFeatureProperty
 from pairs.sim.grid import Grid2D, Grid3D
 from pairs.sim.lattice import ParticleLattice
 from pairs.sim.neighbor_lists import NeighborLists, NeighborListsBuild
 from pairs.sim.pbc import EnforcePBC
 from pairs.sim.properties import PropertiesAlloc, PropertiesResetVolatile
-from pairs.sim.read_from_file import ReadFromFile
+from pairs.sim.read_from_file import ReadParticleData
 from pairs.sim.timestep import Timestep
 from pairs.sim.variables import VariablesDecl
 from pairs.sim.vtk import VTKWrite
@@ -35,6 +37,8 @@ class Simulation:
         self.properties = Properties(self)
         self.vars = Variables(self)
         self.arrays = Arrays(self)
+        self.features = Features(self)
+        self.feature_properties = FeatureProperties(self)
         self.particle_capacity = self.add_var('particle_capacity', Types.Int32, particle_capacity)
         self.nlocal = self.add_var('nlocal', Types.Int32)
         self.nghost = self.add_var('nghost', Types.Int32)
@@ -110,11 +114,11 @@ class Simulation:
         assert self.feature(feature_name) is None, f"Feature already defined: {feature_name}"
         return self.features.add(feature_name)
 
-    def add_feature_property(self, feature_name, prop_name, prop_type):
+    def add_feature_property(self, feature_name, prop_name, prop_type, prop_data):
         feature = self.feature(feature_name)
         assert feature is not None, f"Feature not found: {feature_name}"
         assert self.property(prop_name) is None, f"Property already defined: {prop_name}"
-        return self.properties.add(prop_name, prop_type, value, vol, feature=feature)
+        return self.feature_properties.add(feature, prop_name, prop_type, prop_data)
 
     def property(self, prop_name):
         return self.properties.find(prop_name)
@@ -122,6 +126,12 @@ class Simulation:
     def position(self):
         return self.position_prop
 
+    def feature(self, feature_name):
+        return self.feature.find(feature_name)
+
+    def feature_property(self, feature_prop_name):
+        return self.feature_properties.find(feature_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)
@@ -161,12 +171,17 @@ class Simulation:
         lattice = ParticleLattice(self, grid, spacing, props, positions)
         self.setups.add_statement(lattice)
 
-    def from_file(self, filename, prop_names):
+    def read_particle_data(self, filename, prop_names):
         props = [self.property(prop_name) for prop_name in prop_names]
-        read_object = ReadFromFile(self, filename, props)
+        read_object = ReadParticleData(self, filename, props)
         self.setups.add_statement(read_object)
         self.grid = read_object.grid
 
+    #def read_feature_data(self, filename, feature_name, feature_props):
+    #    feature = self.feature(feature_name)
+    #    props = [self.property(prop_name) for prop_name in feature_props]
+    #    read_object = ReadFeatureData(self, filename, feature, props)
+    #    self.setups.add_statement(read_object)
 
     def build_cell_lists(self, spacing):
         self.cell_lists = CellLists(self, self.grid, spacing, spacing)
@@ -279,6 +294,7 @@ class Simulation:
             VariablesDecl(self),
             ArraysDecl(self),
             PropertiesAlloc(self),
+            FeaturePropertiesAlloc(self)
         ])
 
         program = Module(self, name='main', block=Block.merge_blocks(decls, body))