diff --git a/docs/source/io_formats/mgxs_library.rst b/docs/source/io_formats/mgxs_library.rst
index 9262a234e..2a1e63d43 100644
--- a/docs/source/io_formats/mgxs_library.rst
+++ b/docs/source/io_formats/mgxs_library.rst
@@ -38,10 +38,10 @@ MGXS Library Specification
 The data within <library name> contains the temperature-dependent multi-group
 data for the nuclide or material that it represents.
 
-:Attributes: - **atomic_weight_ratio** (*double*) -- The atomic weight ratio (optional,
-               i.e. it is not meaningful for material-wise data)
-             - **fissionable** (*int*) -- Whether the dataset is fissionable
-               (1) or not (0).
+:Attributes: - **atomic_weight_ratio** (*double*) -- The atomic weight ratio
+               (optional, i.e. it is not meaningful for material-wise data).
+             - **fissionable** (*bool*) -- Whether the dataset is fissionable
+               (True) or not (False).
              - **representation** (*char[]*) -- The method used to generate and
                represent the multi-group cross sections.  That is, whether they
                were generated with scalar flux weighting (or reduced to a
@@ -65,6 +65,15 @@ data for the nuclide or material that it represents.
                or number of points (depending on the value of `scatter_format`)
                used to describe the angular distribution associated with each
                group-to-group transfer probability.
+             - **scatter_shape** (*char[]*) -- The shape of the provided
+               scatter and multiplicity matrix. The values provided are strings
+               describing the ordering the scattering array is provided in
+               row-major (i.e., C/C++ and Python) indexing. Valid values are
+               "[Order][G][G']" or "[Order][G'][G]" where "G'" denotes the
+               secondary/outgoing energy groups, "G" denotes the incoming
+               energy groups, and "Order" is the angular distribution index.
+               This value is not required; if not the default value of
+               "[Order][G][G']" will be assumed.
 
 **/<library name>/kTs/**
 
@@ -79,29 +88,29 @@ Temperature-dependent data, provided for temperature <TTT>K.
 :Datasets: - **total** (*double[]* or *double[][][]*) -- Total cross section.
              This is a 1-D vector if `representation` is "isotropic", or a 3-D
              vector if `representation` is "angle" with dimensions of
-             [groups, azimuthal, polar].
+             [groups][azimuthal][polar].
            - **absorption** (*double[]* or *double[][][]*) -- Absorption
              cross section.
              This is a 1-D vector if `representation` is "isotropic", or a 3-D
              vector if `representation` is "angle" with dimensions of
-             [groups, azimuthal, polar].
+             [groups][azimuthal][polar].
            - **fission** (*double[]* or *double[][][]*) -- Fission
              cross section.
              This is a 1-D vector if `representation` is "isotropic", or a 3-D
              vector if `representation` is "angle" with dimensions of
-             [groups, azimuthal, polar].  This is only required if the dataset
+             [groups][azimuthal][polar].  This is only required if the dataset
              is fissionable and fission-tallies are expected to be used.
            - **kappa-fission** (*double[]* or *double[][][]*) -- Kappa-Fission
              (energy-release from fission) cross section.
              This is a 1-D vector if `representation` is "isotropic", or a 3-D
              vector if `representation` is "angle" with dimensions of
-             [groups, azimuthal, polar].  This is only required if the dataset
+             [groups][azimuthal][polar].  This is only required if the dataset
              is fissionable and fission-tallies are expected to be used.
            - **chi** (*double[]* or *double[][][]*) -- Fission neutron energy
              spectra.
              This is a 1-D vector if `representation` is "isotropic", or a 3-D
              vector if `representation` is "angle" with dimensions of
-             [groups, azimuthal, polar].  This is only required if the dataset
+             [groups][azimuthal][polar].  This is only required if the dataset
              is fissionable and fission-tallies are expected to be used.
            - **nu-fission** (*double[]* to *double[][][][]*) -- Nu-Fission
              cross section.
@@ -119,34 +128,30 @@ Temperature-dependent data, provided for temperature <TTT>K.
 Data specific to neutron scattering for the temperature <TTT>K
 
 :Datasets: - **g_min** (*int[]* or *int[][][]*) --
-             Minimum (most energetic) outgoing groups with non-zero values of
-             the scattering matrix. These group numbers use the standard
+             Minimum (most energetic) groups with non-zero values of
+             the scattering matrix provided.  If `scatter_shape` is
+             "[Order][G][G']" then `g_min` will describe the minimum values
+             of "G'" for each "G"; if `scatter_shape` is "[Order][G'][G]"
+             then `g_min` will describe the minimum values of "G" for each "G'".
+             These group numbers use the standard
              ordering where the fastest neutron energy group is group 1 while
              the slowest neutron energy group is group G.
              The dimensionality of `g_min` is:
-             `g_min[g_in]`, or `g_min[num_polar][num_azimuthal][g_in]`.
+             `g_min[g]`, or `g_min[num_polar][num_azimuthal][g]`.
              The former is used when `representation` is "isotropic", and the
              latter when `representation` is "angle".
            - **g_max** (*int[]* or *int[][][]*) --
-             Maximum (least energetic) outgoing groups with non-zero values of
-             the scattering matrix. These group numbers use the standard
-             ordering where the fastest neutron energy group is group 1 while
-             the slowest neutron energy group is group G.
-             The dimensionality of `g_max` is:
-             `g_max[g_in]`, or `g_max[num_polar][num_azimuthal][g_in]`.
-             The former is used when `representation` is "isotropic", and the
-             latter when `representation` is "angle".
+             Similar to `g_min`, except this dataset describes the maximum
+             (least energetic) groups with non-zero values of
+             the scattering matrix.
            - **scatter_matrix** (*double[]*) -- Flattened representation of the
-             scattering moment matrices. The pre-flattened array is shaped as
-             follows (in row-major format):
-             `scatter_matrix[order(+1)][g_in][g_out]`, or
-             `scatter_matrix[num_polar][num_azimuthal][order(+1)][g_in][g_out]`
-             The former is used when `representation` is "isotropic", and the
-             latter when `representation` is "angle".  Note that if the value of
-             `scatter_format` is "legendre", the order dimension will be one
-             larger than the value of `order`, otherwise it will match `order`.
-             Finally, the g_out dimension has a dimensionality of
-             `g_min` to `g_max`.
+             scattering moment matrices. The pre-flattened array corresponds to
+             the shape provied in `scatter_shape`, but if `representation` is
+             "angle" the dimensionality in `scatter_shape` is prepended by
+             "[num_polar][num_azimuthal]" dimensions. The right-most energy
+             group dimension will only include the entries between `g_min` and
+             `g_max`.
+             dimension has a dimensionality of `g_min` to `g_max`.
            - **multiplicity_matrix** (*double[]*) -- Flattened representation of
              the scattering moment matrices. This dataset provides the code with
              a scaling factor to account for neutrons being produced in (n,xn)
@@ -154,9 +159,6 @@ Data specific to neutron scattering for the temperature <TTT>K
              for every Legendre moment or histogram/tabular bin. This dataset is
              optional, if it is not provided no multiplication (i.e., values of
              1.0) will be assumed.
-             The pre-flattened array is shaped as follows (in row-major format):
-             `multiplicity_matrix[g_in][g_out]`, or
-             `multiplicity_matrix[num_polar][num_azimuthal][g_in][g_out]`
-             The former is used when `representation` is "isotropic", and the
-             latter when `representation` is "angle". Finally, the g_out
-             dimension has a dimensionality of `g_min` to `g_max`.
\ No newline at end of file
+             The pre-flattened array is shapes consistent with `scatter_matrix`
+             except the "[Order]" dimension in `scatter_shape` is ignored since
+             this data is assumed isotropic.
\ No newline at end of file
diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py
index 64210606d..c03aa7877 100644
--- a/openmc/mgxs/library.py
+++ b/openmc/mgxs/library.py
@@ -186,14 +186,14 @@ class Library(object):
 
     @property
     def domains(self):
-        if self._domains is 'all':
-            if self.domain_type is 'material':
+        if self._domains == 'all':
+            if self.domain_type == 'material':
                 return self.openmc_geometry.get_all_materials()
             elif self.domain_type in ['cell', 'distribcell']:
                 return self.openmc_geometry.get_all_material_cells()
-            elif self.domain_type is 'universe':
+            elif self.domain_type == 'universe':
                 return self.openmc_geometry.get_all_universes()
-            elif self.domain_type is 'mesh':
+            elif self.domain_type == 'mesh':
                 raise ValueError('Unable to get domains for Mesh domain type')
             else:
                 raise ValueError('Unable to get domains without a domain type')
@@ -265,7 +265,7 @@ class Library(object):
     @mgxs_types.setter
     def mgxs_types(self, mgxs_types):
         all_mgxs_types = openmc.mgxs.MGXS_TYPES + openmc.mgxs.MDGXS_TYPES
-        if mgxs_types is 'all':
+        if mgxs_types == 'all':
             self._mgxs_types = all_mgxs_types
         else:
             cv.check_iterable_type('mgxs_types', mgxs_types, basestring)
@@ -277,7 +277,7 @@ class Library(object):
     def by_nuclide(self, by_nuclide):
         cv.check_type('by_nuclide', by_nuclide, bool)
 
-        if by_nuclide == True and self.domain_type is 'mesh':
+        if by_nuclide == True and self.domain_type == 'mesh':
             raise ValueError('Unable to create MGXS library by nuclide with '
                              'mesh domain')
 
@@ -287,7 +287,7 @@ class Library(object):
     def domain_type(self, domain_type):
         cv.check_value('domain type', domain_type, openmc.mgxs.DOMAIN_TYPES)
 
-        if self.by_nuclide == True and domain_type is 'mesh':
+        if self.by_nuclide == True and domain_type == 'mesh':
             raise ValueError('Unable to create MGXS library by nuclide with '
                              'mesh domain')
 
@@ -297,21 +297,21 @@ class Library(object):
     def domains(self, domains):
 
         # Use all materials, cells or universes in the geometry as domains
-        if domains is 'all':
+        if domains == 'all':
             self._domains = domains
 
         # User specified a list of material, cell or universe domains
         else:
-            if self.domain_type is 'material':
+            if self.domain_type == 'material':
                 cv.check_iterable_type('domain', domains, openmc.Material)
                 all_domains = self.openmc_geometry.get_all_materials()
             elif self.domain_type in ['cell', 'distribcell']:
                 cv.check_iterable_type('domain', domains, openmc.Cell)
                 all_domains = self.openmc_geometry.get_all_material_cells()
-            elif self.domain_type is 'universe':
+            elif self.domain_type == 'universe':
                 cv.check_iterable_type('domain', domains, openmc.Universe)
                 all_domains = self.openmc_geometry.get_all_universes()
-            elif self.domain_type is 'mesh':
+            elif self.domain_type == 'mesh':
                 cv.check_iterable_type('domain', domains, openmc.Mesh)
 
                 # The mesh and geometry are independent, so set all_domains
@@ -355,7 +355,7 @@ class Library(object):
     def correction(self, correction):
         cv.check_value('correction', correction, ('P0', None))
 
-        if correction is 'P0' and self.legendre_order > 0:
+        if correction == 'P0' and self.legendre_order > 0:
             warn('The P0 correction will be ignored since the scattering '
                  'order "{}" is greater than zero'.format(self.legendre_order))
 
@@ -367,7 +367,7 @@ class Library(object):
         cv.check_greater_than('legendre_order', legendre_order, 0, equality=True)
         cv.check_less_than('legendre_order', legendre_order, 10, equality=True)
 
-        if self.correction is 'P0' and legendre_order > 0:
+        if self.correction == 'P0' and legendre_order > 0:
             msg = 'The P0 correction will be ignored since the scattering ' \
                   'order {} is greater than zero'.format(self.legendre_order)
             warn(msg, RuntimeWarning)
@@ -505,7 +505,7 @@ class Library(object):
         self._openmc_geometry = statepoint.summary.openmc_geometry
         self._nuclides = statepoint.summary.nuclides
 
-        if statepoint.run_mode is 'k-eigenvalue':
+        if statepoint.run_mode == 'k-eigenvalue':
             self._keff = statepoint.k_combined[0]
 
         # Load tallies for each MGXS for each domain and mgxs type
@@ -543,13 +543,13 @@ class Library(object):
 
         """
 
-        if self.domain_type is 'material':
+        if self.domain_type == 'material':
             cv.check_type('domain', domain, (openmc.Material, Integral))
-        elif self.domain_type is 'cell' or self.domain_type is 'distribcell':
+        elif self.domain_type == 'cell' or self.domain_type == 'distribcell':
             cv.check_type('domain', domain, (openmc.Cell, Integral))
-        elif self.domain_type is 'universe':
+        elif self.domain_type == 'universe':
             cv.check_type('domain', domain, (openmc.Universe, Integral))
-        elif self.domain_type is 'mesh':
+        elif self.domain_type == 'mesh':
             cv.check_type('domain', domain, (openmc.Mesh, Integral))
 
         # Check that requested domain is included in library
@@ -662,7 +662,7 @@ class Library(object):
         # Clone this Library to initialize the subdomain-averaged version
         subdomain_avg_library = copy.deepcopy(self)
 
-        if subdomain_avg_library.domain_type is 'distribcell':
+        if subdomain_avg_library.domain_type == 'distribcell':
             subdomain_avg_library.domain_type = 'cell'
         else:
             return subdomain_avg_library
@@ -671,7 +671,7 @@ class Library(object):
         for domain in self.domains:
             for mgxs_type in self.mgxs_types:
                 mgxs = subdomain_avg_library.get_mgxs(domain, mgxs_type)
-                if mgxs.domain_type is 'distribcell':
+                if mgxs.domain_type == 'distribcell':
                     avg_mgxs = mgxs.get_subdomain_avg_xs()
                     subdomain_avg_library.all_mgxs[domain.id][mgxs_type] = avg_mgxs
 
@@ -751,7 +751,7 @@ class Library(object):
             for mgxs_type in self.mgxs_types:
                 mgxs = self.all_mgxs[domain.id][mgxs_type]
 
-                if subdomains is 'avg':
+                if subdomains == 'avg':
                     mgxs = mgxs.get_subdomain_avg_xs()
 
                 mgxs.build_hdf5_store(filename, directory, xs_type=xs_type,
@@ -896,7 +896,7 @@ class Library(object):
 
         # Build & add metadata to XSdata object
         name = xsdata_name
-        if nuclide is not 'total':
+        if nuclide != 'total':
             name += '_' + nuclide
         xsdata = openmc.XSdata(name, self.energy_groups)
 
@@ -912,7 +912,7 @@ class Library(object):
         self.scatter_format = 'legendre'
         self.representation = 'isotropic'
 
-        if nuclide is not 'total':
+        if nuclide != 'total':
             xsdata.atomic_weight_ratio = self._nuclides[nuclide][1]
 
         if subdomain is None:
@@ -921,7 +921,7 @@ class Library(object):
             subdomain = [subdomain]
 
         # Now get xs data itself
-        if 'nu-transport' in self.mgxs_types and self.correction is 'P0':
+        if 'nu-transport' in self.mgxs_types and self.correction == 'P0':
             mymgxs = self.get_mgxs(domain, 'nu-transport')
             xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
                                   subdomain=subdomain)
@@ -1056,7 +1056,7 @@ class Library(object):
         # Initialize file
         mgxs_file = openmc.MGXSLibrary(self.energy_groups)
 
-        if self.domain_type is 'mesh':
+        if self.domain_type == 'mesh':
             # Create the xsdata objects and add to the mgxs_file
             i = 0
             for domain in self.domains:
@@ -1089,7 +1089,7 @@ class Library(object):
                         xsdata_name = 'set' + str(i + 1)
                     else:
                         xsdata_name = xsdata_names[i]
-                    if nuclide is not 'total':
+                    if nuclide != 'total':
                         xsdata_name += '_' + nuclide
 
                     xsdata = self.get_xsdata(domain, xsdata_name,
@@ -1159,14 +1159,14 @@ class Library(object):
         # the multiple meshes could be overlapping or in disparate regions
         # of the continuous energy model. The next step makes sure there is
         # only one before continuing.
-        if self.domain_type is 'mesh':
+        if self.domain_type == 'mesh':
             cv.check_length("domains", self.domains, 1, 1)
 
         # Get the MGXS File Data
         mgxs_file = self.create_mg_library('macro', xsdata_names)
 
         # Now move on the creating the geometry and assigning materials
-        if self.domain_type is 'mesh':
+        if self.domain_type == 'mesh':
             root = openmc.Universe(name='root', universe_id=0)
 
             # Add cells representative of the mesh with reflective BC
@@ -1212,13 +1212,13 @@ class Library(object):
                 materials.append(material)
 
                 # Differentiate Geometry with new Material
-                if self.domain_type is 'material':
+                if self.domain_type == 'material':
                     # Fill all appropriate Cells with new Material
                     for cell in all_cells:
                         if cell.fill.id == domain.id:
                             cell.fill = material
 
-                elif self.domain_type is 'cell':
+                elif self.domain_type == 'cell':
                     for cell in all_cells:
                         if cell.id == domain.id:
                             cell.fill = material
@@ -1283,7 +1283,7 @@ class Library(object):
              ('nu-transport' not in self.mgxs_types))):
             error_flag = True
             warn('A "nu-transport" MGXS type is required since a "P0" '
-                 'correction is applied, but a "nu-transport" MGXS is '
+                 'correction is applied, but a "nu-transport" MGXS == '
                  'not provided.')
         elif (((self.correction is None) and
                ('total' not in self.mgxs_types))):
diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py
index edddd0bfa..35dff1223 100644
--- a/openmc/mgxs_library.py
+++ b/openmc/mgxs_library.py
@@ -8,7 +8,7 @@ import h5py
 import openmc
 import openmc.mgxs
 from openmc.checkvalue import check_type, check_value, check_greater_than, \
-    check_less_than, check_iterable_type
+    check_iterable_type
 
 if sys.version_info[0] >= 3:
     basestring = str
@@ -16,6 +16,7 @@ if sys.version_info[0] >= 3:
 # Supported incoming particle MGXS angular treatment representations
 _REPRESENTATIONS = ['isotropic', 'angle']
 _SCATTER_TYPES = ['tabular', 'legendre', 'histogram']
+_SCATTER_SHAPES = ["[Order][G][G']"]
 
 
 class XSdata(object):
@@ -51,6 +52,8 @@ class XSdata(object):
         Whether or not this is a fissionable data set.
     scatter_format : {'legendre', 'histogram', or 'tabular'}
         Angular distribution representation (legendre, histogram, or tabular)
+    scatter_shapes : {"[Order][G][G']"}
+        Dimensionality of the scattering and multiplicity matrices
     order : int
         Either the Legendre order, number of bins, or number of points used to
         describe the angular distribution associated with each group-to-group
@@ -106,8 +109,8 @@ class XSdata(object):
         :attr:`XSdata.nu_fission` attribute instead.
     nu_fission : dict of numpy.ndarray
         Group-wise fission production cross section vector (i.e., if
-        :attr:`XSdata.chi` is provided), or is the group-wise fission production
-        matrix.
+        :attr:`XSdata.chi` is provided), or is the group-wise fission
+        production matrix.
     inverse_velocities : dict of numpy.ndarray
         Inverse of velocities, in units of sec/cm.
 
@@ -142,6 +145,7 @@ class XSdata(object):
         self._atomic_weight_ratio = None
         self._fissionable = False
         self._scatter_format = 'legendre'
+        self._scatter_shape = "[Order][G][G']"
         self._order = None
         self._num_polar = None
         self._num_azimuthal = None
@@ -184,6 +188,10 @@ class XSdata(object):
     def scatter_format(self):
         return self._scatter_format
 
+    @property
+    def scatter_shape(self):
+        return self._scatter_shape
+
     @property
     def order(self):
         return self._order
@@ -242,28 +250,28 @@ class XSdata(object):
 
     @property
     def vector_shape(self):
-        if self.representation is 'isotropic':
+        if self.representation == 'isotropic':
             return (self.energy_groups.num_groups,)
-        elif self.representation is 'angle':
+        elif self.representation == 'angle':
             return (self.num_polar, self.num_azimuthal,
                     self.energy_groups.num_groups)
 
     @property
     def matrix_shape(self):
-        if self.representation is 'isotropic':
+        if self.representation == 'isotropic':
             return (self.energy_groups.num_groups,
                     self.energy_groups.num_groups)
-        elif self.representation is 'angle':
+        elif self.representation == 'angle':
             return (self.num_polar, self.num_azimuthal,
                     self.energy_groups.num_groups,
                     self.energy_groups.num_groups)
 
     @property
     def pn_matrix_shape(self):
-        if self.representation is 'isotropic':
+        if self.representation == 'isotropic':
             return (self.num_orders, self.energy_groups.num_groups,
                     self.energy_groups.num_groups)
-        elif self.representation is 'angle':
+        elif self.representation == 'angle':
             return (self.num_polar, self.num_azimuthal, self.num_orders,
                     self.energy_groups.num_groups,
                     self.energy_groups.num_groups)
@@ -309,6 +317,12 @@ class XSdata(object):
         check_value('scatter_format', scatter_format, _SCATTER_TYPES)
         self._scatter_format = scatter_format
 
+    @scatter_shape.setter
+    def scatter_shape(self, scatter_shape):
+        # check to see it is of a valid type and value
+        check_value('scatter_shape', scatter_shape, _SCATTER_SHAPES)
+        self._scatter_shape = scatter_shape
+
     @order.setter
     def order(self, order):
         # Check type and value
@@ -694,10 +708,10 @@ class XSdata(object):
         check_value('temperature', temperature, self.temperatures)
 
         i = np.where(self.temperatures == temperature)[0][0]
-        if self.representation is 'isotropic':
+        if self.representation == 'isotropic':
             self._total[i] = total.get_xs(nuclides=nuclide, xs_type=xs_type,
                                           subdomains=subdomain)
-        elif self.representation is 'angle':
+        elif self.representation == 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
             raise ValueError(msg)
 
@@ -740,11 +754,11 @@ class XSdata(object):
         check_value('temperature', temperature, self.temperatures)
 
         i = np.where(self.temperatures == temperature)[0][0]
-        if self.representation is 'isotropic':
+        if self.representation == 'isotropic':
             self._absorption[i] = absorption.get_xs(nuclides=nuclide,
                                                     xs_type=xs_type,
                                                     subdomains=subdomain)
-        elif self.representation is 'angle':
+        elif self.representation == 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
             raise ValueError(msg)
 
@@ -787,11 +801,11 @@ class XSdata(object):
         check_value('temperature', temperature, self.temperatures)
 
         i = np.where(self.temperatures == temperature)[0][0]
-        if self.representation is 'isotropic':
+        if self.representation == 'isotropic':
             self._fission[i] = fission.get_xs(nuclides=nuclide,
                                               xs_type=xs_type,
                                               subdomains=subdomain)
-        elif self.representation is 'angle':
+        elif self.representation == 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
             raise ValueError(msg)
 
@@ -835,11 +849,11 @@ class XSdata(object):
         check_value('temperature', temperature, self.temperatures)
 
         i = np.where(self.temperatures == temperature)[0][0]
-        if self.representation is 'isotropic':
+        if self.representation == 'isotropic':
             self._nu_fission[i] = nu_fission.get_xs(nuclides=nuclide,
                                                     xs_type=xs_type,
                                                     subdomains=subdomain)
-        elif self.representation is 'angle':
+        elif self.representation == 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
             raise ValueError(msg)
 
@@ -892,11 +906,11 @@ class XSdata(object):
         check_value('temperature', temperature, self.temperatures)
 
         i = np.where(self.temperatures == temperature)[0][0]
-        if self.representation is 'isotropic':
+        if self.representation == 'isotropic':
             self._kappa_fission[i] = k_fission.get_xs(nuclides=nuclide,
                                                       xs_type=xs_type,
                                                       subdomains=subdomain)
-        elif self.representation is 'angle':
+        elif self.representation == 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
             raise ValueError(msg)
 
@@ -943,10 +957,10 @@ class XSdata(object):
         check_value('temperature', temperature, self.temperatures)
 
         i = np.where(self.temperatures == temperature)[0][0]
-        if self.representation is 'isotropic':
+        if self.representation == 'isotropic':
             self._chi[i] = chi.get_xs(nuclides=nuclide,
                                       xs_type=xs_type, subdomains=subdomain)
-        elif self.representation is 'angle':
+        elif self.representation == 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
             raise ValueError(msg)
 
@@ -994,7 +1008,7 @@ class XSdata(object):
         check_type('temperature', temperature, Real)
         check_value('temperature', temperature, self.temperatures)
 
-        if self.scatter_format is not 'legendre':
+        if self.scatter_format != 'legendre':
             msg = 'Anisotropic scattering representations other than ' \
                   'Legendre expansions have not yet been implemented in ' \
                   'openmc.mgxs.'
@@ -1011,7 +1025,7 @@ class XSdata(object):
                         [self.order])
 
         i = np.where(self.temperatures == temperature)[0][0]
-        if self.representation is 'isotropic':
+        if self.representation == 'isotropic':
             # Get the scattering orders in the outermost dimension
             self._scatter_matrix[i] = np.zeros((self.num_orders,
                                                 self.energy_groups.num_groups,
@@ -1021,7 +1035,7 @@ class XSdata(object):
                     scatter.get_xs(nuclides=nuclide, xs_type=xs_type,
                                    moment=moment, subdomains=subdomain)
 
-        elif self.representation is 'angle':
+        elif self.representation == 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
             raise ValueError(msg)
 
@@ -1089,7 +1103,7 @@ class XSdata(object):
                         ['universe', 'cell', 'material', 'mesh'])
 
         i = np.where(self.temperatures == temperature)[0][0]
-        if self.representation is 'isotropic':
+        if self.representation == 'isotropic':
             nuscatt = nuscatter.get_xs(nuclides=nuclide,
                                        xs_type=xs_type, moment=0,
                                        subdomains=subdomain)
@@ -1100,7 +1114,7 @@ class XSdata(object):
                                        xs_type=xs_type, moment=0,
                                        subdomains=subdomain)
                 self._multiplicity_matrix[i] = np.divide(nuscatt, scatt)
-        elif self.representation is 'angle':
+        elif self.representation == 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
             raise ValueError(msg)
         self._multiplicity_matrix[i] = \
@@ -1123,7 +1137,7 @@ class XSdata(object):
         if self.representation is not None:
             grp.attrs['representation'] = np.array(self.representation,
                                                    dtype='S')
-            if self.representation is 'angle':
+            if self.representation == 'angle':
                 if self.num_azimuthal is not None:
                     grp.attrs['num_azimuthal'] = self.num_azimuthal
                 if self.num_polar is not None:
@@ -1131,6 +1145,9 @@ class XSdata(object):
         if self.scatter_format is not None:
             grp.attrs['scatter_format'] = np.array(self.scatter_format,
                                                    dtype='S')
+        if self.scatter_shape is not None:
+            grp.attrs['scatter_shape'] = np.array(self.scatter_shape,
+                                                  dtype='S')
         if self.order is not None:
             grp.attrs['order'] = self.order
 
@@ -1170,7 +1187,7 @@ class XSdata(object):
 
             # Get the sparse scattering data to print to the library
             G = self.energy_groups.num_groups
-            if self.representation is 'isotropic':
+            if self.representation == 'isotropic':
                 g_out_bounds = np.zeros((G, 2), dtype=np.int)
                 for g_in in range(G):
                     nz = np.nonzero(self._scatter_matrix[i][0, g_in, :])
@@ -1206,7 +1223,7 @@ class XSdata(object):
                 scatt_grp.create_dataset("g_min", data=g_out_bounds[:, 0])
                 scatt_grp.create_dataset("g_max", data=g_out_bounds[:, 1])
 
-            elif self.representation is 'angle':
+            elif self.representation == 'angle':
                 Np = self.num_polar
                 Na = self.num_azimuthal
                 g_out_bounds = np.zeros((Np, Na, G, 2), dtype=np.int)
diff --git a/scripts/openmc-update-mgxs b/scripts/openmc-update-mgxs
index 6122ff626..658be9fef 100755
--- a/scripts/openmc-update-mgxs
+++ b/scripts/openmc-update-mgxs
@@ -47,14 +47,12 @@ def parse_args():
                         help='input XML file')
     parser.add_argument('-o', '--output', nargs='?', default='',
                         help='output file, in HDF5 format')
-    parser.add_argument('-c', '--compression', type=int,
-                        help='HDF5 Compression Level')
     args = vars(parser.parse_args())
 
-    if args['output'] is '':
+    if args['output'] == '':
         filename = args['input'].name
         extension = filenameos.path.splitext()
-        if extension is '.xml':
+        if extension == '.xml':
             filename = filename[:filename.rfind('.')] + '.h5'
         args['output'] = filename
 
@@ -75,7 +73,7 @@ def get_data(element, entry):
     return value
 
 
-if __name__ is '__main__':
+if __name__ == '__main__':
     args = parse_args()
 
     # Parse the XML data.
@@ -117,7 +115,7 @@ if __name__ is '__main__':
         representation = get_data(xsdata_elem, 'representation')
         if representation is None:
             representation = 'isotropic'
-        if representation is 'angle':
+        if representation == 'angle':
             n_azi = int(get_data(xsdata_elem, 'num_azimuthal'))
             n_pol = int(get_data(xsdata_elem, 'num_polar'))
 
@@ -146,14 +144,14 @@ if __name__ is '__main__':
                                      representation=representation))
             if awr is not None:
                 xsd[-1].atomic_weight_ratio = awr
-            if representation is 'angle':
+            if representation == 'angle':
                 xsd[-1].num_azimuthal = n_azi
                 xsd[-1].num_polar = n_pol
             xsd[-1].scatter_format = scatter_format
             xsd[-1].order = order
             names.append(name)
 
-        if scatter_format is 'legendre':
+        if scatter_format == 'legendre':
             order_dim = order + 1
         else:
             order_dim = order
diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90
index ddce50be3..cbdf470b0 100644
--- a/src/mgxs_header.F90
+++ b/src/mgxs_header.F90
@@ -338,6 +338,13 @@ module mgxs_header
       else
         this % scatter_format = ANGLE_LEGENDRE
       end if
+      if (attribute_exists(xs_id, "scatter_shape")) then
+        call read_attribute(temp_str, xs_id, "scatter_shape")
+        temp_str = trim(temp_str)
+        if (to_lower(temp_str) /= "[order][g][g']") then
+          call fatal_error("Invalid scatter_shape option!")
+        end if
+      end if
       if (attribute_exists(xs_id, "fissionable")) then
         call read_attribute(this % fissionable, xs_id, "fissionable")
       else
diff --git a/tests/1d_mgxs.h5 b/tests/1d_mgxs.h5
index 401fe02fc..3befba194 100644
Binary files a/tests/1d_mgxs.h5 and b/tests/1d_mgxs.h5 differ