diff --git a/CADETMatch/scores/fractionation.py b/CADETMatch/scores/fractionation.py
index 4897227..6287e86 100644
--- a/CADETMatch/scores/fractionation.py
+++ b/CADETMatch/scores/fractionation.py
@@ -32,7 +32,7 @@ def run(sim_data, feature):
         selected = (times >= start) & (times <= stop)
 
         local_times = times[selected]
-        local_values = simulation.root.output.solution.unit_001["solution_outlet_comp_%03d" % component][selected]
+        local_values = simulation.root.output.solution[feature['unit']["solution_outlet_comp_%03d" % component][selected]
 
         sim_value = numpy.trapz(local_values, local_times)
 
@@ -87,6 +87,7 @@ def setup(sim, feature, selectedTimes, selectedValues, CV_time, abstol):
 
             funcs.append( (start, stop, int(component), value, func) )
     temp['funcs'] = funcs
+    temp['unit'] = feature['unit_name']
     settings.count = len(funcs)
     return temp
 
diff --git a/CADETMatch/scores/fractionationCombine.py b/CADETMatch/scores/fractionationCombine.py
index c5dfeb2..8c12717 100644
--- a/CADETMatch/scores/fractionationCombine.py
+++ b/CADETMatch/scores/fractionationCombine.py
@@ -40,7 +40,7 @@ def run(sim_data, feature):
         selected = (times >= start) & (times <= stop)
 
         local_times = times[selected]
-        local_values = simulation.root.output.solution.unit_001["solution_outlet_comp_%03d" % component][selected]
+        local_values = simulation.root.output.solution[feature['unit']["solution_outlet_comp_%03d" % component][selected]
 
         sim_value = numpy.trapz(local_values, local_times)
 
@@ -104,6 +104,7 @@ def setup(sim, feature, selectedTimes, selectedValues, CV_time, abstol):
     temp['funcs'] = funcs
     temp['components'] = [int(i) for i in headers[2:]]
     temp['samplesPerComponent'] = rows
+    temp['unit'] = feature['unit_name']
     return temp
 
 def headers(experimentName, feature):
diff --git a/CADETMatch/scores/fractionationMeanVariance.py b/CADETMatch/scores/fractionationMeanVariance.py
index bb3b170..715cd02 100644
--- a/CADETMatch/scores/fractionationMeanVariance.py
+++ b/CADETMatch/scores/fractionationMeanVariance.py
@@ -35,7 +35,7 @@ def run(sim_data, feature):
     for (start, stop, component, values, func_mean_time, func_variance_time, func_mean_value, func_variance_value) in funcs:
         time_center = (start + stop)/2.0
                 
-        sim_values = util.fractionate(start, stop, times, simulation.root.output.solution.unit_001["solution_outlet_comp_%03d" % component])
+        sim_values = util.fractionate(start, stop, times, simulation.root.output.solution[feature['unit']["solution_outlet_comp_%03d" % component])
        
         mean_sim_time, variance_sim_time, skew_sim_time, mean_sim_value, variance_sim_value, skew_sim_value = util.fracStat(time_center, sim_values)
 
@@ -93,6 +93,7 @@ def setup(sim, feature, selectedTimes, selectedValues, CV_time, abstol):
     temp['funcs'] = funcs
     temp['components'] = [int(i) for i in headers[2:]]
     temp['samplesPerComponent'] = rows
+    temp['unit'] = feature['unit_name']
     return temp
 
 def headers(experimentName, feature):
diff --git a/CADETMatch/scores/fractionationMoment.py b/CADETMatch/scores/fractionationMoment.py
index f8aef65..8991a48 100644
--- a/CADETMatch/scores/fractionationMoment.py
+++ b/CADETMatch/scores/fractionationMoment.py
@@ -35,7 +35,7 @@ def run(sim_data, feature):
     for (start, stop, component, values, func_mean_time, func_variance_time, func_skew_time, func_mean_value, func_variance_value, func_skew_value) in funcs:
         time_center = (start + stop)/2.0
                 
-        sim_values = util.fractionate(start, stop, times, simulation.root.output.solution.unit_001["solution_outlet_comp_%03d" % component])
+        sim_values = util.fractionate(start, stop, times, simulation.root.output.solution[feature['unit']["solution_outlet_comp_%03d" % component])
        
         mean_sim_time, variance_sim_time, skew_sim_time, mean_sim_value, variance_sim_value, skew_sim_value = util.fracStat(time_center, sim_values)
 
@@ -99,6 +99,7 @@ def setup(sim, feature, selectedTimes, selectedValues, CV_time, abstol):
     temp['funcs'] = funcs
     temp['components'] = [int(i) for i in headers[2:]]
     temp['samplesPerComponent'] = rows
+    temp['unit'] = feature['unit_name']
     return temp
 
 def headers(experimentName, feature):
diff --git a/CADETMatch/scores/fractionationSSE.py b/CADETMatch/scores/fractionationSSE.py
index aae85e6..3e5c58a 100644
--- a/CADETMatch/scores/fractionationSSE.py
+++ b/CADETMatch/scores/fractionationSSE.py
@@ -37,7 +37,7 @@ def run(sim_data, feature):
         selected = (times >= start) & (times <= stop)
 
         local_times = times[selected]
-        local_values = simulation.root.output.solution.unit_001["solution_outlet_comp_%03d" % component][selected]
+        local_values = simulation.root.output.solution[feature['unit']]["solution_outlet_comp_%03d" % component][selected]
 
         sim_value = numpy.trapz(local_values, local_times)
 
@@ -76,6 +76,7 @@ def setup(sim, feature, selectedTimes, selectedValues, CV_time, abstol):
 
     temp['start'] = start
     temp['stop'] = stop
+    temp['unit'] = feature['unit_name']
 
     smallestTime = min(data['Stop'] - data['Start'])
     abstolFraction = abstol * smallestTime
diff --git a/CADETMatch/scores/fractionationSlide.py b/CADETMatch/scores/fractionationSlide.py
index b92866e..00d144c 100644
--- a/CADETMatch/scores/fractionationSlide.py
+++ b/CADETMatch/scores/fractionationSlide.py
@@ -75,7 +75,7 @@ def run(sim_data, feature):
 
     for component, value_func in funcs:
         exp_values = numpy.array(data[str(component)])
-        sim_value = simulation.root.output.solution.unit_001["solution_outlet_comp_%03d" % component]
+        sim_value = simulation.root.output.solution[feature['unit']["solution_outlet_comp_%03d" % component]
 
         rollLeft, rollRight, searchMax = rollRange(times, sim_value, searchIndexStart, searchIndexStop)
 
@@ -148,6 +148,7 @@ def setup(sim, feature, selectedTimes, selectedValues, CV_time, abstol):
     temp['samplesPerComponent'] = rows
     temp['CV_time'] = CV_time
     temp['funcs'] = funcs
+    temp['unit'] = feature['unit_name']
     return temp
 
 def headers(experimentName, feature):
diff --git a/CADETMatch/util.py b/CADETMatch/util.py
index df0d958..a122031 100644
--- a/CADETMatch/util.py
+++ b/CADETMatch/util.py
@@ -1401,12 +1401,16 @@ def setupSimulation(sim, times):
 
     sim.root.input.solver.user_solution_times = times
     sim.root.input.solver.sections.section_times[-1] = times[-1]
-    sim.root.input['return'].unit_001.write_solution_particle = 0
-    sim.root.input['return'].unit_001.write_solution_column_inlet = 1
-    sim.root.input['return'].unit_001.write_solution_column_outlet = 1
-    sim.root.input['return'].unit_001.write_solution_inlet = 1
-    sim.root.input['return'].unit_001.write_solution_outlet = 1
-    sim.root.input['return'].unit_001.split_components_data = 0
+
+    for i in range(sim.root.input.model.nunits):
+        unit = 'unit_%03d' % i
+
+        sim.root.input['return'][unit].write_solution_particle = 0
+        sim.root.input['return'][unit].write_solution_column_inlet = 1
+        sim.root.input['return'][unit].write_solution_column_outlet = 1
+        sim.root.input['return'][unit].write_solution_inlet = 1
+        sim.root.input['return'][unit].write_solution_outlet = 1
+        sim.root.input['return'][unit].split_components_data = 0
     sim.root.input.solver.nthreads = 1
 
 def graph_corner_process(cache, last=False, interval=1200):