Adding isomeric branching to chain

openmc/deplete/chain.py

@@ -120,6 +120,187 @@
 __all__ = ["Chain", "REACTIONS"]
 
 
+def _read_activation_products(evaluation):
+    """Read MF=9/MF=10 data from an ENDF evaluation to get isomeric products.
+
+    For each reaction (MT) that has MF=9 or MF=10 data, returns the product
+    states with their energy-dependent yields (branching fractions).
+
+    Parameters
+    ----------
+    evaluation : openmc.data.endf.Evaluation
+        ENDF evaluation for a nuclide
+
+    Returns
+    -------
+    dict
+        ``{mt: {(Z, A, m): yield_func}}`` where yield_func is a
+        :class:`openmc.data.Tabulated1D` giving the yield as a function of
+        energy. m is the isomeric state number (0=ground, 1=first metastable,
+        etc.), derived by sequentially numbering the excited states in order
+        of increasing LFS.
+    """
+    results = {}
+
+    # Find which MTs have MF=8 (which signals MF=9/MF=10 data)
+    mf8_mts = set()
+    for mf, mt, nc, mod in evaluation.reaction_list:
+        if mf == 8 and mt not in (454, 457, 459):
+            mf8_mts.add(mt)
+
+    for mt in sorted(mf8_mts):
+        if (8, mt) not in evaluation.section:
+            continue
+
+        file_obj = StringIO(evaluation.section[8, mt])
+        items = openmc.data.endf.get_head_record(file_obj)
+        n_states = items[4]
+        decay_sublib = (items[5] == 1)
+
+        # Determine which MF files (9 or 10) are present
+        present = {9: False, 10: False}
+        for _ in range(n_states):
+            if decay_sublib:
+                items = openmc.data.endf.get_cont_record(file_obj)
+            else:
+                items, _ = openmc.data.endf.get_list_record(file_obj)
+            lmf = int(items[2])
+            if lmf in present:
+                present[lmf] = True
+
+        products = {}
+
+        # Read MF=3 cross section for this MT (needed for MF=10 yield calc)
+        neutron_xs = None
+        if present[10] and (3, mt) in evaluation.section:
+            file_obj_3 = StringIO(evaluation.section[3, mt])
+            openmc.data.endf.get_head_record(file_obj_3)
+            _, neutron_xs = openmc.data.endf.get_tab1_record(file_obj_3)
+
+        for mf in (9, 10):
+            if not present[mf]:
+                continue
+            if (mf, mt) not in evaluation.section:
+                continue
+
+            file_obj = StringIO(evaluation.section[mf, mt])
+            items = openmc.data.endf.get_head_record(file_obj)
+            n_states = items[4]
+
+            for _ in range(n_states):
+                items, xs = openmc.data.endf.get_tab1_record(file_obj)
+                Z, A = divmod(int(items[2]), 1000)
+                lfs = int(items[3])
+
+                if mf == 9:
+                    # MF=9 gives multiplicity (yield) directly
+                    products[(Z, A, lfs)] = xs
+                elif neutron_xs is not None:
+                    # MF=10 gives production XS; convert to yield by
+                    # dividing by the reaction cross section
+                    energy = np.union1d(xs.x, neutron_xs.x)
+                    prod_xs = xs(energy)
+                    nxs = neutron_xs(energy)
+                    idx = np.where(nxs > 0)
+                    yield_ = np.zeros_like(energy)
+                    yield_[idx] = prod_xs[idx] / nxs[idx]
+                    products[(Z, A, lfs)] = \
+                        openmc.data.Tabulated1D(energy, yield_)
+
+        if products:
+            # Remap LFS level numbers to sequential isomeric state numbers.
+            # LFS in MF9/MF10 refers to the energy level, not the isomeric
+            # state number.  Isomeric states are numbered sequentially:
+            # ground (LFS=0) -> m=0, then m=1, m=2, ... for excited states
+            # in order of increasing LFS.
+            za_groups: dict[tuple[int, int], list[int]] = {}
+            for (Z, A, lfs) in products:
+                za_groups.setdefault((Z, A), []).append(lfs)
+
+            remapped = {}
+            for (Z, A), lfs_values in za_groups.items():
+                lfs_values.sort()
+                for m, lfs in enumerate(lfs_values):
+                    remapped[(Z, A, m)] = products[(Z, A, lfs)]
+
+            results[mt] = remapped
+
+    return results
+
+
+def _add_isomeric_reactions(nuclide, rx_name, q_value, iso_products,
+                            decay_data, missing_rx_product, parent):
+    """Add reaction entries for each isomeric product state.
+
+    Uses MF=9/MF=10 data to create separate reaction entries for ground
+    and metastable product states with appropriate branching ratios.
+    Also stores the full energy-dependent yield data on the nuclide
+    for use in spectrum-dependent branching ratio collapse.
+
+    Parameters
+    ----------
+    nuclide : openmc.deplete.Nuclide
+        Nuclide to add reactions to
+    rx_name : str
+        Reaction name, e.g. '(n,gamma)'
+    q_value : float
+        Q value of the reaction in [eV]
+    iso_products : dict
+        ``{(Z, A, m): yield_func}`` from :func:`_read_activation_products`
+    decay_data : dict
+        Decay data keyed by nuclide name
+    missing_rx_product : list
+        List to append missing products to
+    parent : str
+        Parent nuclide name
+    """
+    # Map product states to daughter names, filtering missing nuclides
+    yield_funcs = {}  # {daughter_name: yield_func}
+    ground_state = None
+    for (Z, A, m), yield_func in iso_products.items():
+        daughter = gnds_name(Z, A, m)
+
+        if daughter not in decay_data:
+            orig_daughter = daughter
+            daughter = replace_missing(daughter, decay_data)
+            if daughter is None:
+                missing_rx_product.append((parent, rx_name, orig_daughter))
+                continue
+
+        # If multiple states map to same daughter, keep the one with
+        # larger max yield (shouldn't normally happen)
+        if daughter in yield_funcs:
+            if np.max(yield_func.y) > np.max(yield_funcs[daughter].y):
+                yield_funcs[daughter] = yield_func
+        else:
+            yield_funcs[daughter] = yield_func
+
+        if m == 0:
+            ground_state = daughter
+
+    if not yield_funcs:
+        return
+
+    # Build common energy grid from union of all product grids
+    common_energies = yield_funcs[next(iter(yield_funcs))].x
+    for yf in yield_funcs.values():
+        common_energies = np.union1d(common_energies, yf.x)
+
+    # Interpolate all yields onto common grid
+    yield_arrays = {}
+    for daughter, yf in yield_funcs.items():
+        yield_arrays[daughter] = np.clip(yf(common_energies), 0.0, None)
+
+    # Store energy-dependent data on the nuclide
+    nuclide.isomeric_branching[rx_name] = (common_energies, yield_arrays)
+
+    # Add a single reaction entry for the ground state target.
+    # The isomeric_branching data is the source of truth for
+    # the actual products and energy-dependent branching.
+    target = ground_state or next(iter(yield_funcs))
+    nuclide.add_reaction(rx_name, target, q_value, 1.0)
+
+
 def replace_missing(product, decay_data):
     """Replace missing product with suitable decay daughter.
 
@@ -361,6 +542,7 @@ def from_endf(cls, decay_files, fpy_files, neutron_files,
         if progress:
             print('Processing neutron sub-library files...')
         reactions = {}
+        activation_products = {}
         for f in neutron_files:
             evaluation = openmc.data.endf.Evaluation(f)
             name = evaluation.gnds_name
@@ -372,6 +554,11 @@ def from_endf(cls, decay_files, fpy_files, neutron_files,
                     q_value = openmc.data.endf.get_cont_record(file_obj)[1]
                     reactions[name][mt] = q_value
 
+            # Read MF=9/MF=10 activation products to get isomeric
+            # branching data for each reaction
+            activation_products[name] = _read_activation_products(
+                evaluation)
+
         # Determine what decay and FPY nuclides are available
         if progress:
             print('Processing decay sub-library files...')
@@ -438,19 +625,11 @@ def from_endf(cls, decay_files, fpy_files, neutron_files,
             fissionable = False
             if parent in reactions:
                 reactions_available = set(reactions[parent].keys())
+                parent_products = activation_products.get(parent, {})
                 for name in transmutation_reactions:
                     mts = REACTIONS[name].mts
                     delta_A, delta_Z = openmc.data.DADZ[name]
                     if mts & reactions_available:
-                        A = data.nuclide['mass_number'] + delta_A
-                        Z = data.nuclide['atomic_number'] + delta_Z
-                        daughter = f'{openmc.data.ATOMIC_SYMBOL[Z]}{A}'
-
-                        if daughter not in decay_data:
-                            daughter = replace_missing(daughter, decay_data)
-                            if daughter is None:
-                                missing_rx_product.append((parent, name, daughter))
-
                         # Store Q value
                         for mt in sorted(mts):
                             if mt in reactions[parent]:
@@ -459,7 +638,36 @@ def from_endf(cls, decay_files, fpy_files, neutron_files,
                         else:
                             q_value = 0.0
 
-                        nuclide.add_reaction(name, daughter, q_value, 1.0)
+                        # Check if MF=9/MF=10 provides isomeric product
+                        # states for this reaction
+                        iso_products = {}
+                        for mt in mts:
+                            if mt in parent_products:
+                                iso_products = parent_products[mt]
+                                break
+
+                        if iso_products:
+                            # Use MF=9/MF=10 data to resolve isomeric
+                            # products with branching ratios
+                            _add_isomeric_reactions(
+                                nuclide, name, q_value, iso_products,
+                                decay_data, missing_rx_product, parent)
+                        else:
+                            # Fall back to ground state only
+                            A = data.nuclide['mass_number'] + delta_A
+                            Z = data.nuclide['atomic_number'] + delta_Z
+                            daughter = \
+                                f'{openmc.data.ATOMIC_SYMBOL[Z]}{A}'
+
+                            if daughter not in decay_data:
+                                daughter = replace_missing(
+                                    daughter, decay_data)
+                                if daughter is None:
+                                    missing_rx_product.append(
+                                        (parent, name, daughter))
+
+                            nuclide.add_reaction(
+                                name, daughter, q_value, 1.0)
 
                 if any(mt in reactions_available for mt in openmc.data.FISSION_MTS):
                     q_value = reactions[parent][18]

openmc/deplete/nuclide.py

@@ -122,6 +122,11 @@ def __init__(self, name=None):
         # Reaction paths
         self.reactions = []
 
+        # Energy-dependent isomeric branching data.
+        # {rx_type: (energies, {target: yields})} where energies is a 1D
+        # array and yields are 1D arrays on the same grid.
+        self.isomeric_branching = {}
+
         # Decay sources
         self.sources = {}
 
@@ -270,6 +275,17 @@ def from_xml(cls, element, root=None, fission_q=None):
             nuc.reactions.append(ReactionTuple(
                 r_type, target, Q, branching_ratio))
 
+        # Read energy-dependent isomeric branching data
+        for iso_elem in element.iter('isomeric_branching'):
+            rx_type = get_text(iso_elem, 'reaction')
+            energies = np.fromstring(iso_elem.find('energies').text, sep=' ')
+            products = {}
+            for p_elem in iso_elem.iter('product'):
+                target = get_text(p_elem, 'nuclide')
+                yields = np.fromstring(p_elem.text, sep=' ')
+                products[target] = yields
+            nuc.isomeric_branching[rx_type] = (energies, products)
+
         fpy_elem = element.find('neutron_fission_yields')
         if fpy_elem is not None:
             # Check for use of FPY from other nuclide
@@ -331,6 +347,18 @@ def to_xml_element(self):
             if br != 1.0:
                 rx_elem.set('branching_ratio', str(br))
 
+        # Write energy-dependent isomeric branching data
+        if self.isomeric_branching:
+            for rx_type, (energies, products) in self.isomeric_branching.items():
+                iso_elem = ET.SubElement(elem, 'isomeric_branching')
+                iso_elem.set('reaction', rx_type)
+                e_elem = ET.SubElement(iso_elem, 'energies')
+                e_elem.text = ' '.join(str(e) for e in energies)
+                for target, yields in products.items():
+                    p_elem = ET.SubElement(iso_elem, 'product')
+                    p_elem.set('nuclide', target)
+                    p_elem.text = ' '.join(str(y) for y in yields)
+
         if self.yield_data:
             fpy_elem = ET.SubElement(elem, 'neutron_fission_yields')
 

tests/unit_tests/test_deplete_chain.py

@@ -587,3 +587,63 @@ def test_reduce(gnd_simple_chain, endf_chain):
     reduced_chain = endf_chain.reduce(['U235'])
     assert 'H1' in reduced_chain
     assert 'H2' in reduced_chain
+
+
+_ISOMERIC_CHAIN = """\
+<depletion_chain>
+  <nuclide name="Am241" half_life="13651800000.0" decay_modes="1"
+           decay_energy="5627985.0" reactions="1">
+    <decay type="alpha" target="Np237" branching_ratio="1.0"/>
+    <reaction type="(n,gamma)" Q="5537755.0" target="Am242"/>
+    <isomeric_branching reaction="(n,gamma)">
+      <energies>1e-5 0.0253 1e3 1e5 1e6 1e7</energies>
+      <product nuclide="Am242">0.9 0.9 0.87 0.84 0.74 0.52</product>
+      <product nuclide="Am242_m1">0.1 0.1 0.13 0.16 0.26 0.48</product>
+    </isomeric_branching>
+  </nuclide>
+  <nuclide name="Am242" half_life="57672.0" decay_modes="1"
+           decay_energy="0.0" reactions="0">
+    <decay type="beta-" target="Cm242" branching_ratio="1.0"/>
+  </nuclide>
+  <nuclide name="Am242_m1" half_life="4449312000.0" decay_modes="1"
+           decay_energy="0.0" reactions="0">
+    <decay type="IT" target="Am242" branching_ratio="1.0"/>
+  </nuclide>
+  <nuclide name="Np237" reactions="0"/>
+  <nuclide name="Cm242" reactions="0"/>
+</depletion_chain>
+"""
+
+
+def test_isomeric_branching_chain_roundtrip():
+    """Test that a chain with isomeric branching data can be read, written,
+    and re-read with energy-dependent yields preserved."""
+    with cdtemp():
+        with open('chain_iso.xml', 'w') as fh:
+            fh.write(_ISOMERIC_CHAIN)
+
+        chain = Chain.from_xml('chain_iso.xml')
+        am241 = chain['Am241']
+
+        # Verify isomeric branching was read
+        assert '(n,gamma)' in am241.isomeric_branching
+        energies, products = am241.isomeric_branching['(n,gamma)']
+        assert len(energies) == 6
+        assert 'Am242' in products
+        assert 'Am242_m1' in products
+
+        # Check energy-dependent values
+        assert np.isclose(products['Am242'][0], 0.9)
+        assert np.isclose(products['Am242'][-1], 0.52)
+        assert np.isclose(products['Am242_m1'][0], 0.1)
+        assert np.isclose(products['Am242_m1'][-1], 0.48)
+
+        # Write and re-read
+        chain.export_to_xml('chain_iso_rt.xml')
+        chain2 = Chain.from_xml('chain_iso_rt.xml')
+        am241_2 = chain2['Am241']
+
+        energies2, products2 = am241_2.isomeric_branching['(n,gamma)']
+        assert np.allclose(energies, energies2)
+        for target in products:
+            assert np.allclose(products[target], products2[target])