diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py
index ae8ea8191f..e60878d601 100644
--- a/openmc/data/__init__.py
+++ b/openmc/data/__init__.py
@@ -14,3 +14,4 @@ from .nbody import *
 from .thermal import *
 from .urr import *
 from .library import *
+from .fission_energy import *
diff --git a/openmc/data/endf_utils.py b/openmc/data/endf_utils.py
new file mode 100644
index 0000000000..6db3c611cf
--- /dev/null
+++ b/openmc/data/endf_utils.py
@@ -0,0 +1,43 @@
+"""This module contains a few utility functions for reading ENDF_ data. It is by
+no means enough to read an entire ENDF file.  For a more complete ENDF reader,
+see Pyne_.
+
+.. _ENDF: http://www.nndc.bnl.gov/endf
+.. _Pyne: http://www.pyne.io
+
+"""
+
+import re
+
+def read_float(float_string):
+    """Parse ENDF 6E11.0 formatted string into a float."""
+    assert len(float_string) == 11
+    pattern = '([\s\\-]\d+\\.\d+)([\\+\\-]\d+)'
+    mantissa, exponent = re.match(pattern, float_string).groups()
+    return float(mantissa + 'e' + exponent)
+
+
+def read_CONT_line(line):
+    """Parse 80-column line from ENDF CONT record into floats and ints."""
+    return (read_float(line[0:11]), read_float(line[11:22]), int(line[22:33]),
+            int(line[33:44]), int(line[44:55]), int(line[55:66]),
+            int(line[66:70]), int(line[70:72]), int(line[72:75]),
+            int(line[75:80]))
+
+def identify_nuclide(fname):
+    """Read the header of an ENDF file and extract identifying information."""
+    with open(fname, 'r') as fh:
+        # Skip the tape id (TPID).
+        line = fh.readline()
+
+        # Read the first HEAD and CONT info.
+        line = fh.readline()
+        ZA, AW, LRP, LFI, NLIB, NMOD, MAT, MF, MT, NS = read_CONT_line(line)
+        line = fh.readline()
+        ELIS, STA, LIS, LISO, junk, NFOR, MAT, MF, MT, NS = read_CONT_line(line)
+
+    # Return dictionary of the most important identifying information.
+    return {'Z': int(ZA) // 1000,
+            'A': int(ZA) % 1000,
+            'LIS': LIS,
+            'LISO': LISO}
diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py
new file mode 100644
index 0000000000..5716a3e858
--- /dev/null
+++ b/openmc/data/fission_energy.py
@@ -0,0 +1,285 @@
+from collections import Callable
+import sys
+#from warnings import warn
+
+import numpy as np
+from numpy.polynomial.polynomial import Polynomial
+
+from .function import Tabulated1D, Sum
+from .endf_utils import read_float, read_CONT_line, identify_nuclide
+import openmc.checkvalue as cv
+
+if sys.version_info[0] >= 3:
+    basestring = str
+
+
+class FissionEnergyRelease(object):
+    def __init__(self):
+        self._fragments = None
+        self._prompt_neutrons = None
+        self._delayed_neutrons = None
+        self._prompt_photons = None
+        self._delayed_photons = None
+        self._betas = None
+        self._neutrinos = None
+        self._form = None
+
+    @property
+    def fragments(self):
+        return self._fragments
+
+    @property
+    def prompt_neutrons(self):
+        return self._prompt_neutrons
+
+    @property
+    def delayed_neutrons(self):
+        return self._delayed_neutrons
+
+    @property
+    def prompt_photons(self):
+        return self._prompt_photons
+
+    @property
+    def delayed_photons(self):
+        return self._delayed_photons
+
+    @property
+    def betas(self):
+        return self._betas
+
+    @property
+    def neutrinos(self):
+        return self._neutrinos
+
+    @property
+    def recoverable(self):
+        return Sum([self.fragments, self.prompt_neutrons, self.delayed_neutrons,
+                    self.prompt_photons, self.delayed_photons, self.betas])
+
+    @property
+    def total(self):
+        return Sum([self.fragments, self.prompt_neutrons, self.delayed_neutrons,
+                    self.prompt_photons, self.delayed_photons, self.betas,
+                    self.neutrinos])
+    
+    @property
+    def form(self):
+        return self._form
+
+    @fragments.setter
+    def fragments(self, energy_release):
+        cv.check_type('fragments', energy_release, Callable)
+        self._fragments = energy_release
+
+    @prompt_neutrons.setter
+    def prompt_neutrons(self, energy_release):
+        cv.check_type('prompt_neutrons', energy_release, Callable)
+        self._prompt_neutrons = energy_release
+
+    @delayed_neutrons.setter
+    def delayed_neutrons(self, energy_release):
+        cv.check_type('delayed_neutrons', energy_release, Callable)
+        self._delayed_neutrons = energy_release
+
+    @prompt_photons.setter
+    def prompt_photons(self, energy_release):
+        cv.check_type('prompt_photons', energy_release, Callable)
+        self._prompt_photons = energy_release
+
+    @delayed_photons.setter
+    def delayed_photons(self, energy_release):
+        cv.check_type('delayed_photons', energy_release, Callable)
+        self._delayed_photons = energy_release
+
+    @betas.setter
+    def betas(self, energy_release):
+        cv.check_type('betas', energy_release, Callable)
+        self._betas = energy_release
+
+    @neutrinos.setter
+    def neutrinos(self, energy_release):
+        cv.check_type('neutrinos', energy_release, Callable)
+        self._neutrinos = energy_release
+
+    @form.setter
+    def form(self, form):
+        cv.check_value('format', form, ('Madland', 'Sher-Beck'))
+        self._form = form
+
+    @classmethod
+    def from_endf(cls, filename, incident_neutron):
+        """Generate fission energy release data from an ENDF file.
+
+        Parameters
+        ----------
+        filename : str
+            Name of the ENDF file containing fission energy release data
+
+        incident_neutron : openmc.data.IncidentNeutron
+            Corresponding incident neutron dataset
+
+        Returns
+        -------
+        openmc.data.FissionEnergyRelease
+            Fission energy release data
+
+        """
+
+        # Check to make sure this ENDF file matches the expected isomer.
+        ident = identify_nuclide(filename)
+        if ident['Z'] != incident_neutron.atomic_number:
+            pass
+        if ident['A'] != incident_neutron.mass_number:
+            pass
+        if ident['LISO'] != incident_neutron.metastable:
+            pass
+
+        # Extract the MF=1, MT=458 section.
+        lines = []
+        with open(filename, 'r') as fh:
+            line = fh.readline()
+            while line != '':
+                if line[70:75] == ' 1458':
+                    lines.append(line)
+                line = fh.readline()
+
+        # Read the number of coefficients in this LIST record.
+        NPL = read_CONT_line(lines[1])[4]
+
+        # Parse the ENDF LIST into an array.
+        data = []
+        for i in range(NPL):
+            row, column = divmod(i, 6)
+            data.append(read_float(lines[2 + row][11*column:11*(column+1)]))
+
+        # Declare the coefficient names and the order they are given in.  The
+        # LIST contains a value followed immediately by an uncertainty for each
+        # of these components, times the polynomial order + 1.  If we only find
+        # one value for each of these components, then we need to use the
+        # Sher-Beck formula for energy dependence.  Otherwise, it is a
+        # polynomial.
+        labels = ('EFR', 'ENP', 'END', 'EGP', 'EGD', 'EB', 'ENU', 'ER', 'ET')
+
+        # Associate each set of values and uncertainties with its label.
+        value = dict()
+        uncertainty = dict()
+        for i in range(len(labels)):
+            value[labels[i]] = data[2*i::18]
+            uncertainty[labels[i]] = data[2*i + 1::18]
+
+        # In ENDF/B-7.1, data for 2nd-order coefficients were mistakenly not
+        # converted from MeV to eV.  Check for this error and fix it if present.
+        n_coeffs = len(value['EFR'])
+        if n_coeffs == 3:  # Only check 2nd-order data.
+            # Check each energy component for the error.  If a 1 MeV neutron
+            # causes a change of more than 100 MeV, we know something is wrong.
+            error_present = False
+            for coeffs in value.values():
+                second_order = coeffs[2]
+                if abs(second_order) * 1e12 > 1e8:
+                    error_present = True
+                    break
+
+            # If we found the error, reduce all 2nd-order coeffs by 10**6.
+            if error_present:
+                for coeffs in value.values(): coeffs[2] *= 1e-6
+                for coeffs in uncertainty.values(): coeffs[2] *= 1e-6
+
+            # Perform the sanity check again... just in case.
+            for coeffs in value.values():
+                second_order = coeffs[2]
+                if abs(second_order) * 1e12 > 1e8:
+                    raise ValueError("Encountered a ludicrously large second-"
+                                     "order polynomial coefficient.")
+
+        # Convert eV to MeV.
+        for coeffs in value.values():
+            for i in range(len(coeffs)):
+                coeffs[i] *= 10**(-6 + 6*i)
+        for coeffs in uncertainty.values():
+            for i in range(len(coeffs)):
+                coeffs[i] *= 10**(-6 + 6*i)
+
+        out = cls()
+        if n_coeffs > 1:
+            out.form = 'Madland'
+            out.fragments = Polynomial(value['EFR'])
+            out.prompt_neutrons = Polynomial(value['ENP'])
+            out.delayed_neutrons = Polynomial(value['END'])
+            out.prompt_photons = Polynomial(value['EGP'])
+            out.delayed_photons = Polynomial(value['EGD'])
+            out.betas = Polynomial(value['EB'])
+            out.neutrinos = Polynomial(value['ENU'])
+        else:
+            out.form = 'Sher-Beck'
+            raise NotImplemented
+
+        return out
+
+    @classmethod
+    def from_hdf5(cls, group):
+        """Generate fission energy release data from an HDF5 group.
+
+        Parameters
+        ----------
+        group : h5py.Group
+            HDF5 group to read from
+
+        Returns
+        -------
+        openmc.data.FissionEnergyRelease
+            Fission energy release data
+
+        """
+
+        obj = cls()
+        if group.attrs['format'] == 'Madland':
+            obj.fragments = Polynomial(group['fragments'].value)
+            obj.prompt_neutrons = Polynomial(group['prompt_neutrons'].value)
+            obj.delayed_neutrons = Polynomial(group['delayed_neutrons'].value)
+            obj.prompt_photons = Polynomial(group['prompt_photons'].value)
+            obj.delayed_photons = Polynomial(group['delayed_photons'].value)
+            obj.betas = Polynomial(group['betas'].value)
+            obj.neutrinos = Polynomial(group['neutrinos'].value)
+        elif group.attrs['format'] == 'Sher-Beck':
+            raise NotImplemented
+        else:
+            raise ValueError('Unrecognized energy release format')
+
+        return obj
+
+    def to_hdf5(self, group):
+        """Write energy release data to an HDF5 group
+
+        Parameters
+        ----------
+        group : h5py.Group
+            HDF5 group to write to
+
+        """
+
+        if self.form == 'Madland':
+            group.attrs['format'] = np.string_('Madland')
+            group.create_dataset('fragments', data=self.fragments.coef)
+            group.create_dataset('prompt_neutrons',
+                                 data=self.prompt_neutrons.coef)
+            group.create_dataset('delayed_neutrons',
+                                 data=self.delayed_neutrons.coef)
+            group.create_dataset('prompt_photons',
+                                 data=self.prompt_photons.coef)
+            group.create_dataset('delayed_photons',
+                                 data=self.delayed_photons.coef)
+            group.create_dataset('betas', data=self.betas.coef)
+            group.create_dataset('neutrinos', data=self.neutrinos.coef)
+        elif self.form == 'Sher-Beck':
+            group.attrs['format'] = np.string_('Sher-Beck')
+            self.fragments.to_hdf5(group, 'fragments')
+            self.prompt_neutrons.to_hdf5(group, 'prompt_neutrons')
+            self.delayed_neutrons.to_hdf5(group, 'delayed_neutrons')
+            self.prompt_photons.to_hdf5(group, 'prompt_photons')
+            self.delayed_photons.to_hdf5(group, 'delayed_photons')
+            self.betas.to_hdf5(group, 'betas')
+            self.neutrinos.to_hdf5(group, 'neutrinos')
+        else:
+            raise ValueError('Unrecognized energy release format')
diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py
index 63b2bab01b..84d9bcb2c6 100644
--- a/openmc/data/neutron.py
+++ b/openmc/data/neutron.py
@@ -9,6 +9,7 @@ import h5py
 
 from .data import ATOMIC_SYMBOL, SUM_RULES
 from .ace import Table, get_table
+from .fission_energy import FissionEnergyRelease
 from .function import Tabulated1D, Sum
 from .product import Product
 from .reaction import Reaction, _get_photon_products
@@ -81,6 +82,7 @@ class IncidentNeutron(object):
         self.temperature = temperature
 
         self._energy = None
+        self._fission_energy = None
         self.reactions = OrderedDict()
         self.summed_reactions = OrderedDict()
         self.urr = None
@@ -126,6 +128,10 @@ class IncidentNeutron(object):
     def energy(self):
         return self._energy
 
+    @property
+    def fission_energy(self):
+        return self._fission_energy
+
     @property
     def temperature(self):
         return self._temperature
@@ -186,6 +192,12 @@ class IncidentNeutron(object):
         cv.check_type('energy grid', energy, Iterable, Real)
         self._energy = energy
 
+    @fission_energy.setter
+    def fission_energy(self, fission_energy):
+        cv.check_type('fission energy release', fission_energy,
+                      FissionEnergyRelease)
+        self._fission_energy = fission_energy
+
     @reactions.setter
     def reactions(self, reactions):
         cv.check_type('reactions', reactions, Mapping)
@@ -276,6 +288,11 @@ class IncidentNeutron(object):
             urr_group = g.create_group('urr')
             self.urr.to_hdf5(urr_group)
 
+        # Write fission energy release data
+        if self.fission_energy is not None:
+            fer_group = g.create_group('fission_energy_release')
+            self.fission_energy.to_hdf5(fer_group)
+
         f.close()
 
     @classmethod
@@ -331,6 +348,11 @@ class IncidentNeutron(object):
             urr_group = group['urr']
             data.urr = ProbabilityTables.from_hdf5(urr_group)
 
+        # Read fission energy release data
+        if 'fission_energy_release' in group:
+            fer_group = group['fission_energy_release']
+            data.fission_energy = FissionEnergyRelease.from_hdf5(fer_group)
+
         return data
 
     @classmethod