diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py
index 689546b86bf..d1d62841f5e 100644
--- a/openmc/deplete/chain.py
+++ b/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]
diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py
index 95881483474..636468b5bb0 100644
--- a/openmc/deplete/nuclide.py
+++ b/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')
diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py
index e90b6102241..ebb8599f44f 100644
--- a/tests/unit_tests/test_deplete_chain.py
+++ b/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 = """\
+
+
+
+
+
+ 1e-5 0.0253 1e3 1e5 1e6 1e7
+ 0.9 0.9 0.87 0.84 0.74 0.52
+ 0.1 0.1 0.13 0.16 0.26 0.48
+
+
+
+
+
+
+
+
+
+
+
+"""
+
+
+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])
diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py
index f2bb7d1b65e..750a1120d33 100644
--- a/tests/unit_tests/test_deplete_nuclide.py
+++ b/tests/unit_tests/test_deplete_nuclide.py
@@ -348,3 +348,68 @@ def test_deepcopy():
# Mutate the original and verify the copy remains intact
nuc *= 2
assert copied_nuc != nuc
+
+
+def test_isomeric_branching_xml_roundtrip():
+ """Test reading and writing energy-dependent isomeric branching data."""
+
+ # Build a nuclide with isomeric branching data
+ nuc = nuclide.Nuclide("Am241")
+ nuc.add_reaction("(n,gamma)", "Am242", 5537755.0, 1.0)
+
+ energies = np.array([1e-5, 0.0253, 1e3, 1e5, 1e6, 1e7])
+ products = {
+ "Am242": np.array([0.9, 0.9, 0.87, 0.84, 0.74, 0.52]),
+ "Am242_m1": np.array([0.1, 0.1, 0.13, 0.16, 0.26, 0.48]),
+ }
+ nuc.isomeric_branching["(n,gamma)"] = (energies, products)
+
+ # Write to XML
+ elem = nuc.to_xml_element()
+
+ # Check XML structure
+ iso_elems = elem.findall("isomeric_branching")
+ assert len(iso_elems) == 1
+ assert iso_elems[0].get("reaction") == "(n,gamma)"
+ assert iso_elems[0].find("energies") is not None
+ product_elems = iso_elems[0].findall("product")
+ assert len(product_elems) == 2
+ assert {p.get("nuclide") for p in product_elems} == {"Am242", "Am242_m1"}
+
+ # Read back from XML
+ nuc2 = nuclide.Nuclide.from_xml(elem)
+ assert "(n,gamma)" in nuc2.isomeric_branching
+ energies2, products2 = nuc2.isomeric_branching["(n,gamma)"]
+ assert np.allclose(energies, energies2)
+ for target in products:
+ assert target in products2
+ assert np.allclose(products[target], products2[target])
+
+
+def test_isomeric_branching_from_xml():
+ """Test reading isomeric branching data from XML string."""
+
+ data = """
+
+
+
+ 1e6 5e6 1e7 2e7
+ 0.6 0.55 0.5 0.45
+ 0.15 0.17 0.2 0.22
+ 0.25 0.28 0.3 0.33
+
+
+ """
+
+ elem = ET.fromstring(data)
+ nuc = nuclide.Nuclide.from_xml(elem)
+
+ assert nuc.name == "Ir191"
+ assert "(n,2n)" in nuc.isomeric_branching
+ energies, products = nuc.isomeric_branching["(n,2n)"]
+ assert len(energies) == 4
+ assert np.isclose(energies[0], 1e6)
+ assert set(products.keys()) == {"Ir190", "Ir190_m1", "Ir190_m2"}
+ assert np.allclose(products["Ir190"], [0.6, 0.55, 0.5, 0.45])
+ assert np.allclose(products["Ir190_m1"], [0.15, 0.17, 0.2, 0.22])
+ assert np.allclose(products["Ir190_m2"], [0.25, 0.28, 0.3, 0.33])