diff --git a/openmc/data/function.py b/openmc/data/function.py
index e2efb4c10..bea6f5e9a 100644
--- a/openmc/data/function.py
+++ b/openmc/data/function.py
@@ -80,50 +80,46 @@ class Tabulated1D(object):
         # Get indices for interpolation
         idx = np.searchsorted(self.x, x, side='right') - 1
 
-        # Find lowest valid index
-        i_low = np.searchsorted(idx, 0)
-
+        # Loop over interpolation regions
         for k in range(len(self.breakpoints)):
-            # Determine which x values are within this interpolation range
-            i_high = np.searchsorted(idx, self.breakpoints[k] - 1)
+            # Get indices for the begining and ending of this region
+            i_begin = self.breakpoints[k-1] - 1 if k > 0 else 0
+            i_end = self.breakpoints[k] - 1
 
-            # Get x values and bounding (x,y) pairs
-            xk = x[i_low:i_high]
-            xi = self.x[idx[i_low:i_high]]
-            xi1 = self.x[idx[i_low:i_high] + 1]
-            yi = self.y[idx[i_low:i_high]]
-            yi1 = self.y[idx[i_low:i_high] + 1]
+            # Figure out which idx values lie within this region
+            contained = (idx >= i_begin) & (idx < i_end)
+
+            xk = x[contained]                 # x values in this region
+            xi = self.x[idx[contained]]       # low edge of corresponding bins
+            xi1 = self.x[idx[contained] + 1]  # high edge of corresponding bins
+            yi = self.y[idx[contained]]
+            yi1 = self.y[idx[contained] + 1]
 
             if self.interpolation[k] == 1:
                 # Histogram
-                y[i_low:i_high] = yi
+                y[contined] = yi
 
             elif self.interpolation[k] == 2:
                 # Linear-linear
-                y[i_low:i_high] = yi + (xk - xi)/(xi1 - xi)*(yi1 - yi)
+                y[contained] = yi + (xk - xi)/(xi1 - xi)*(yi1 - yi)
 
             elif self.interpolation[k] == 3:
                 # Linear-log
-                y[i_low:i_high] = yi + np.log(xk/xi)/np.log(xi1/xi)*(yi1 - yi)
+                y[contained] = yi + np.log(xk/xi)/np.log(xi1/xi)*(yi1 - yi)
 
             elif self.interpolation[k] == 4:
                 # Log-linear
-                y[i_low:i_high] = yi*np.exp((xk - xi)/(xi1 - xi)*np.log(yi1/yi))
+                y[contained] = yi*np.exp((xk - xi)/(xi1 - xi)*np.log(yi1/yi))
 
             elif self.interpolation[k] == 5:
                 # Log-log
-                y[i_low:i_high] = yi*np.exp(np.log(xk/xi)/np.log(xi1/xi)*np.log(yi1/yi))
+                y[contained] = (yi*np.exp(np.log(xk/xi)/np.log(xi1/xi)
+                                *np.log(yi1/yi)))
 
-            i_low = i_high
-
-        # In some cases, the first/last point of x may be less than the first
-        # value of self.x due only to precision, so we check if they're close
-        # and set them equal if so. Otherwise, the interpolated value might be
-        # out of range (and thus zero)
-        if np.isclose(x[0], self.x[0], 1e-8):
-            y[0] = self.y[0]
-        if np.isclose(x[-1], self.x[-1], 1e-8):
-            y[-1] = self.y[-1]
+        # In some cases, x values might be outside the tabulated region due only
+        # to precision, so we check if they're close and set them equal if so.
+        y[np.isclose(x, self.x[0], atol=1e-14)] = self.y[0]
+        y[np.isclose(x, self.x[-1], atol=1e-14)] = self.y[-1]
 
         return y if iterable else y[0]
 
diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py
index 63b2bab01..7cbcca8c7 100644
--- a/openmc/data/neutron.py
+++ b/openmc/data/neutron.py
@@ -326,6 +326,17 @@ class IncidentNeutron(object):
                     tgroup = group['total_nu']
                     rx.derived_products.append(Product.from_hdf5(tgroup))
 
+        # Build summed reactions.  Start from the highest MT number because high
+        # MTs never depend on lower MTs.
+        for mt_sum in sorted(SUM_RULES.keys())[::-1]:
+            if mt_sum not in data:
+                xs_components = [data[mt].xs for mt in SUM_RULES[mt_sum]
+                                 if mt in data]
+                if len(xs_components) > 0:
+                    rxn = Reaction(mt_sum)
+                    rxn.xs = Sum(xs_components)
+                    data.summed_reactions[mt_sum] = rxn
+
         # Read unresolved resonance probability tables
         if 'urr' in group:
             urr_group = group['urr']
diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py
index 62b02048b..ad707d276 100644
--- a/openmc/data/reaction.py
+++ b/openmc/data/reaction.py
@@ -400,7 +400,7 @@ class Reaction(object):
         # Read cross section
         if 'xs' in group:
             xs = group['xs'].value
-            rx.xs = Tabulated1D(energy, xs)
+            rx.xs = Tabulated1D(energy[rx.threshold_idx:], xs)
 
         # Determine number of products
         n_product = 0