Including metastable states as reactions in chainfiles

[shimwell]

Bug Description

Currently we can produce chain files with reactions and decays. The decay options can contain nuclides ending with _m* as I guess this comes from the ENDF file directly.

However I don't think we can include the production of metastable states from reactions. For example production of Co58m is not possible and I think this is quite an important one for steel activation in fusion

It looks like the chain.REACTIONS currently includes reactions that change the number of protons and neutrons but not the meta_stable state

I've tried calling openmc.deplete.Chain.from_endf with '(n,n1)' reaction and as it is not in the REACTIONS dictionary I get a key error.

Would be be acceptable if we include some more reactions in the chain.REACTIONS dictionary

I notice that data.REACTION_NAME has a line of code to add lots of n,n* reactions

openmc/openmc/data/reaction.py

Line 59 in 847e9b5

REACTION_NAME.update({i: '(n,n{})'.format(i - 50) for i in range(51, 91)})

Perhaps we can add something like the line above to the chain.REACTIONS?

I think if we want to accommodated reactions that result in a metastable then the ReactionInfo namedtuple might need another term. It is currently defined with mts, dadz and secondaries but perhaps it should have another term for meta.
To accommodate this reaction we might what to change from dAdZ to dAdZdM in the named tuple

Steps to Reproduce

chain = openmc.deplete.Chain.from_endf(
        decay_files=decay_files,
        fpy_files=fpy_files,
        neutron_files=neutron_files,
        reactions=['(n,n1)']
    )

Environment

Ubuntu 18.04, python 3.11, current openmc-dev (soon to be 0.14.0)

[paulromano]

Production of metastable states from reactions is complicated. When we have a reaction like (n,n5), that means that it is producing the 5th excited state. However, that doesn't mean it is a metastable state and in fact most excited states are very short-lived. To give a concrete example, Am242 has excited states at 44 keV and and 48.6 keV. The first excited state decays instantly but the second excited state is metastable, i.e., it is the first metastable state (Am242_m1). So, the reaction Am242(n,n2) should produce Am242_m1. Higher MTs like (n,n3) could also produce Am242_m1 if the higher excited states decay into the second excited state, but now you need to know the probability of that transition.

Hypothetically, if you have the full level information from, say, ENSDF, you could do this systematically for inelastic scattering reactions. However, 1) we don't have ENSDF data available to us and 2) ENDF contains some of this information but not all, and 3) even if you have this information, it doesn't help you for other reactions like (n,γ) because you don't naturally know the probability of producing a given excited state of the compound nucleus. ENDF also allows evaluators to put radionuclide production cross sections in MF=9/MF=10 so that one doesn't really need to think about level diagrams and transition probabilities, but the use of MF=9/MF=10 is sporadic in the major libraries (I believe special purpose activation libraries rely entirely on MF=9/MF=10).

tl;dr The change you proposed in #2763 is not going to solve this problem, and as far as I can tell there is no "easy" fix.

[shimwell]

Thanks for all the info, hope to find a nice solution one day but sounds like we need to start will the nuclear data

[shimwell]

Could we "hard code" a few key nuclides like Co58_m if we know their specifc production route

[paulromano]

Yes, hard-coding is how we handle capture branching ratios (see depletion/add_branching_ratios.py on the data repo). The deplete.Chain class has a set_branch_ratios method that sets capture branching ratios, so the first step is probably to extend that to handles reactions other than (n,γ).

[shimwell]

Thanks for the info Paul, this is starting to make sense.

There part I'm least sure about is where do the contents of the json files in the data repo come from. Wondering how this file and this file were created? Is there a python script that grabs the data from somewhere

[shimwell]

I have made a PR to the data repo attempting to make the adding_branch_ratios flexible so that other (non n,gamma) reactions can be supported

[jbae11]

Hello @shimwell and @paulromano from the past. I recently ran into this issue too with some activation validation studies, where OpenMC was not matching with the results of FISPACT and ORIGEN, mostly due to the metastable (and their daughter) inventories. I looked into this issue on my own time and found a potential solution, and I wanted to run by y'all to see how sane it sounds:

Issue

1. Metastable generation paths are not properly implemented in openmc.deplete

• (e.g., Am241 -> (n,gamma) -> Am242m, or Ag107 -> (n,gamma) -> Ag108m)
• This is because one reaction only has one target, as listed in the chain file:
  • <reaction type="(n,gamma)" Q="5539101.0" target="Am242"/>
• This yield is also energy-dependent, which further complicates issue.
• In other words, (n,gamma) reaction of Am241, should yield either Am242 or Am242m, depending on the incident neutron spectrum.

How others approach

1. ORIGEN implements an energy-dependent partial cross sections (https://www.ornl.gov/file/origen-data-0/display#page=10). So there's technically two (n,gamma) reactions for Am241 that are total (n,gamma) cross sections multiplied by the energy-dependent branching ratio, as such:

• 

Solution Implementation

• Please note the implementation is (hopefully) optimized to minimize codebase perturbation and ensure backwards compatibility

1. Additional data in Chainfile - leveraging neutron_fission_yields

Currently, neutron_fission_yields are defined in the chainfile as a set of:

• energy points
• list of product and yield for each energy bin

A similar dataset can be provided in the Chainfile for energy-dependent-reaction-yields, as such:

 # fission yield looks like this:
 <neutron_fission_yields>
      <energies>0.0253 500000.0 14000000.0</energies>
      <fission_yields energy="0.0253">
        <products>Ag103 Ag105 ...</products>
        <data>0.0 0.0 0.0 0.0 ...</data>
      </fission_yields>
      <fission_yields energy="500000.0">
        <products>Ag103 Ag105 ...</products>
        <data>0.0 0.0 ...</data>
      </fission_yields>
      <fission_yields energy="14000000.0">
        <products>Ag103 Ag105 ...</products>
        <data>0.0 0.0 ...</data>
      </fission_yields>
  </neutron_fission_yields>

# for Ag107
# the default target can be overwritten, or used as 1 outside energy bounds
<reaction type="(n,gamma)" Q="7269690.0" target="Ag108"/>
      <yield_data type="(n,gamma)">
        <energies>0.0253 500000.0 14000000.0</energies>
        <yields energy="0.0253">
          <products>Ag108 Ag108_m1</products>
          <data>0.1 0.9 </data>
        </yields>
        <yields energy="500000.0">
          <products>Ag108 Ag108_m1</products>
          <data>0.1 0.9 </data>
        </yields>
        <yields energy="14000000.0">
          <products>Ag108 Ag108_m1</products>
          <data>0.1 0.9 </data>
        </yields>
      </yield_data>
      

2. Implementation to openmc.deplete

First, we'd have to implement getting the neutron-flux-spectrum-weighted yield for each split reaction, similar to fission_yields calculated in openmc.deplete.OpenMCOperator.

Now with that data, all you'd have to do is to change the term on the transmutation matrix that's created by openmc.deplete.Chain.form_matrix(), which is, again, similar to how the fission yields are implemented.

Obivously,

Devil's in the details. I ran into the following issues:

1. This branching data is weird..?

• For most branching issues like Am241 (n,gamma) and Ag107 (n,gamma) splits, I was able to find them in the TENDL2019 ENDF libraries by using @paulromano endf module, something like:

    for f in files:
        mat = endf.Material(f)
        evaluation = openmc.data.endf.Evaluation(f)
        iso = evaluation.gnds_name
        if iso == 'Es205': continue
        section = mat.section_data
        sub_dict = {}
        for k, v in section.items():
            if k[0] == 9: # this is the MF, and MF=9 is production
                # product level
                mt = k[1]
                assert((10, mt) not in section.keys())
                if mt in [4, 18]: # production of neutron
                    continue
                name = mt_to_name[mt]
                levels = section[k]['levels']
                if len(levels) > 1:
                    print(iso, name)
                    sub_dict[name] = levels

• but the issue is if MF=10, also listed as production, is in a weird format and I'm not sure how to tackle it.
  • I found the crucial Ag109 -> (n,2n) -> Ag108 or Ag108m reaction here, which is also included in ORIGEN

2. The flux-spectrum-normalization of the yield data can be challenging.

• the yield data in MF=9 has varying energy bins, even within a reaction:
  • 
• this could be an easy fix using smart interpolation and checking to ensure they sum up to 1 at every point.

TL;DR

Leveraging the existing fission_yields logic in openmc.deplete.Chain may make implementation of 'one-reaction-multiple-targets' capability. But nuclear data is weird and hard to automatically gracefully parse through.

Here's a small result for the FNS activation benchmark for silver, where I manually changed the transmutation matrix to have the multiple paths using the one-group-collapsed yield using the ENDF python tool above:

[paulromano]

Thanks for your attention to this @jbae11! Very interesting (and creative) solution there with using the fission yields logic. That being said, I don't think I would be in favor of adopting that as a permanent/preferred solution. Perhaps in a pinch you could try to force things into that abstraction but it just doesn't feel natural.

To respond to a few particular comments:

This is because one reaction only has one target, as listed in the chain file

This is sort of true and sort of not. It's true that each reaction only has one target, but you can list the same reaction twice with different targets and associated branching ratios. We use this currently to treat capture branching ratios for select reactions. This could be done for reactions other than (n,gamma) too but we just haven't put in data yet. Of course, as you've pointed out, those branching ratios are spectrum dependent. Sometimes, the ratio between the levels is mostly flat (e.g., Ni58(n,p)) in which case a constant value is fine, but in other cases (e.g., Ag108(n,2n)) the ratio does vary so then you are forced to pick a value and hope it's correct.

This branching data is weird..?

Totally! There are at least three different ways that evaluators have given to specify cross sections to different levels of the compound nucleus, which means it's really hard to know what the right procedure is to try to automatically determine them:

1. Evaluator gives cross sections to different levels in the MF=9 or MF=10 sections (common in activation libraries like EAF)
2. Evaluator gives ground and metastable states as different "products" using MF=6
3. Evaluator uses MT=600+ to specify cross sections to different excited states of the residual nucleus (one of which may be metastable)

Oh the joys of nuclear data!

[jbae11]

joys of nuclear data indeed. Thanks for the discussion and insight. I think given the currently existing infrastructure of the branching ratios, maybe we can add a function to modify Chain object to tack on the branching ratios from a multigroup flux. Can I ask how the branching_ratios_pwr.json was generated?

I'm working on taking the energy-dependent branching ratios from ORIGEN data and testing it for the FNS benchmark

[paulromano]

Those PWR branching ratios come from here. Further commentary on depletion chain data can also be found here.

[shimwell]

adding support for energy-dependent branching ratios in the chain xml files and openmc looks like a necessary feature that helps us get a more accurate/better answer for these depletion simulations.

Is this feature something we can architect/designed and then afterwards look for reliable methods of populating the branching ratio data from nuclear data?

Separately would it be worth making DT fusion specific branching ratios file for fusion, something equivalent to the sfr and pwr files but perhaps using 14MeV to determine the ratios. This would give us something we can use with the existing openmc workflow and get better results in the short term.

I think @eepeterson and group are also interested in this topic