Using RateFilter to Search Rates in a Library¶

The Library class in pynucastro provides a high level interface for reading files containing one or more Reaclib rates and then filtering these rates based on user-specified criteria for the nuclei involved in the reactions. We can then use the resulting rates to build a network.

This example uses a Reaclib snapshot downloaded from:

https://groups.nscl.msu.edu/jina/reaclib/db/library.php?action=viewsnapshots.

Reading a rate snapshot¶

The Library class will look for the library file in the working directory or in the pynucastro/library subdirectory of the pynucastro package.

When the constructor is supplied a file name, pynucastro will read the contents of this file and interpret them as Reaclib rates in either the Reaclib 1 or 2 formats. The Library then stores the rates from the file as Rate objects.

[1]:

%matplotlib inline
import pynucastro as pyna

[2]:

library_file = '20180201ReaclibV2.22'

[3]:

mylibrary = pyna.rates.Library(library_file)

Filtering the Library¶

This example introduces the RateFilter class which allows us to define a set of reactants and products to search for in a Library object.

Inexact Filtering¶

Inexact filtering is like using wildcards: in the following example, the rate filter we define will match any rates in which \(\mathrm{^{12}C}\) is a reactant.

[4]:

c12_inexact_filter = pyna.rates.RateFilter(reactants=['c12'], exact=False)

Once we construct a RateFilter object, we can apply it to our Library by passing it to the Library.filter function.

Library.filter returns a new Library object containing the rates that match our RateFilter.

We can print a Library to see the rates it contains. In parentheses the rate identifier is printed, showing the Reaclib rate label as well as whether the rate is forward or reverse.

[5]:

c12_inexact_library = mylibrary.filter(c12_inexact_filter)
print(c12_inexact_library)

c12 --> n + c11    (c12 --> n + c11 <bb92_reaclib_reverse>)
c12 --> p + b11    (c12 --> p + b11 <nw00_reaclib_reverse>)
c12 --> he4 + he4 + he4    (c12 --> he4 + he4 + he4 <fy05_reaclib_reverse>)
c12 + n --> c13    (c12 + n --> c13 <ka02_reaclib>)
c12 + p --> n13    (c12 + p --> n13 <ls09_reaclib>)
c12 + he4 --> o16    (c12 + he4 --> o16 <nac2_reaclib>)
c12 + n --> p + b12    (c12 + n --> p + b12 <wag_reaclib_reverse>)
c12 + n --> he4 + be9    (c12 + n --> he4 + be9 <cf88_reaclib_reverse>)
c12 + he4 --> n + o15    (c12 + he4 --> n + o15 <cf88_reaclib_reverse>)
c12 + he4 --> p + n15    (c12 + he4 --> p + n15 <nacr_reaclib_reverse>)
c12 + c12 --> n + mg23    (c12 + c12 --> n + mg23 <cf88_reaclib_reverse>)
c12 + c12 --> p + na23    (c12 + c12 --> p + na23 <cf88_reaclib>)
c12 + c12 --> he4 + ne20    (c12 + c12 --> he4 + ne20 <cf88_reaclib>)
o16 + c12 --> n + si27    (o16 + c12 --> n + si27 <cf88_reaclib_reverse>)
o16 + c12 --> p + al27    (o16 + c12 --> p + al27 <cf88_reaclib>)
o16 + c12 --> he4 + mg24    (o16 + c12 --> he4 + mg24 <cf88_reaclib>)
ne20 + c12 --> n + s31    (ne20 + c12 --> n + s31 <rolf_reaclib>)
ne20 + c12 --> p + p31    (ne20 + c12 --> p + p31 <rolf_reaclib>)
ne20 + c12 --> he4 + si28    (ne20 + c12 --> he4 + si28 <rolf_reaclib>)

The rate identifiers above can be used to access individual Rate objects within a Library as follows:

[6]:

cago = c12_inexact_library.get_rate('c12 + he4 --> o16 <nac2_reaclib>')

[7]:

cago.plot()

_images/library-examples-filtering_12_0.png

Exact Filtering¶

Exact filtering is useful when you have a specific rate in mind or a specific combination of reactants or products. In the following example, we look for all rates of the form \(\mathrm{^{12}C + ^{12}C \rightarrow \ldots}\)

To use exact filtering, omit the exact keyword to the RateFilter constructor, as it is turned on by default.

Exact filtering does not mean all the nuclei involved in the rate must be specified, it means that all filtering options passed to the RateFilter constructor are strictly applied. In this case, the filter will return rates with exactly two reactants, both of which are \(\mathrm{^{12}C}\). However, the filter places no constraint on the products or number of products in the rate.

[8]:

c12_exact_filter = pyna.rates.RateFilter(reactants=['c12', 'c12'])
c12_exact_library = mylibrary.filter(c12_exact_filter)
print(c12_exact_library)

c12 + c12 --> n + mg23    (c12 + c12 --> n + mg23 <cf88_reaclib_reverse>)
c12 + c12 --> p + na23    (c12 + c12 --> p + na23 <cf88_reaclib>)
c12 + c12 --> he4 + ne20    (c12 + c12 --> he4 + ne20 <cf88_reaclib>)

Example: Building an Alpha Capture Network¶

In the next example, we use rate filtering to iteratively construct a Library containing the alpha capture rates linking \(\mathrm{^{12}C}\) to \(\mathrm{^{56}Ni}\).

After finding each successive link in the alpha capture chain, we call Library.heaviest() to find the heaviest nucleus in the filtered rates. This corresponds to the nucleus with the largest mass number, and in case of a tie between isobars, this returns the isobar with the smallest atomic number. We use this feature to find the reverse rate for each alpha capture reaction.

In the example below, we add each filtered library to our alpha capture library alpha_library, initialized as an empty Library. The Library class supports the addition operator by returning a new library containing the rates in the two libraries we added together.

This example also introduces the max_products keyword, which specifies we are looking for reactions producing at most max_products product nuclei.

Similarly, the RateFilter constructor supports the following keywords constraining the number of reactants and products:

min_reactants
max_reactants
min_products
max_products

Because we have omitted the argument exact=False, the filter constraints we apply are exact.

[9]:

alpha_library = pyna.rates.Library()

capture = pyna.rates.Nucleus('he4')
seed = pyna.rates.Nucleus('c12')

while True:
    ac_filter = pyna.rates.RateFilter(reactants=[capture, seed], max_products=1)
    ac_library = mylibrary.filter(ac_filter)
    alpha_library = alpha_library + ac_library

    heavy = ac_library.heaviest()
    ac_filter_inv = pyna.rates.RateFilter(reactants=[heavy], products=[capture, seed])
    ac_inv_library = mylibrary.filter(ac_filter_inv)
    alpha_library = alpha_library + ac_inv_library

    print(heavy)
    if heavy.A == 56:
        break
    else:
        seed = heavy

o16
ne20
mg24
si28
s32
ar36
ca40
ti44
cr48
fe52
ni56

We will next print out the library we constructed, seeing that we have both forward and reverse rates for the alpha chain.

Note that at this time pynucastro has not yet implemented nuclear partition functions, so these reverse rates are calculated only from detailed balance in the Reaclib library.

[10]:

print(alpha_library)

c12 + he4 --> o16    (c12 + he4 --> o16 <nac2_reaclib>)
o16 --> he4 + c12    (o16 --> he4 + c12 <nac2_reaclib_reverse>)
o16 + he4 --> ne20    (o16 + he4 --> ne20 <co10_reaclib>)
ne20 --> he4 + o16    (ne20 --> he4 + o16 <co10_reaclib_reverse>)
ne20 + he4 --> mg24    (ne20 + he4 --> mg24 <il10_reaclib>)
mg24 --> he4 + ne20    (mg24 --> he4 + ne20 <il10_reaclib_reverse>)
mg24 + he4 --> si28    (mg24 + he4 --> si28 <st08_reaclib>)
si28 --> he4 + mg24    (si28 --> he4 + mg24 <st08_reaclib_reverse>)
si28 + he4 --> s32    (si28 + he4 --> s32 <ths8_reaclib>)
s32 --> he4 + si28    (s32 --> he4 + si28 <ths8_reaclib_reverse>)
s32 + he4 --> ar36    (s32 + he4 --> ar36 <ths8_reaclib>)
ar36 --> he4 + s32    (ar36 --> he4 + s32 <ths8_reaclib_reverse>)
ar36 + he4 --> ca40    (ar36 + he4 --> ca40 <ths8_reaclib>)
ca40 --> he4 + ar36    (ca40 --> he4 + ar36 <ths8_reaclib_reverse>)
ca40 + he4 --> ti44    (ca40 + he4 --> ti44 <chw0_reaclib>)
ti44 --> he4 + ca40    (ti44 --> he4 + ca40 <chw0_reaclib_reverse>)
ti44 + he4 --> cr48    (ti44 + he4 --> cr48 <ths8_reaclib>)
cr48 --> he4 + ti44    (cr48 --> he4 + ti44 <ths8_reaclib_reverse>)
cr48 + he4 --> fe52    (cr48 + he4 --> fe52 <ths8_reaclib>)
fe52 --> he4 + cr48    (fe52 --> he4 + cr48 <ths8_reaclib_reverse>)
fe52 + he4 --> ni56    (fe52 + he4 --> ni56 <ths8_reaclib>)
ni56 --> he4 + fe52    (ni56 --> he4 + fe52 <ths8_reaclib_reverse>)

Next we can create a reaction network from our filtered alpha capture library by passing our library to a network constructor using the libraries keyword.

[11]:

alpha_network = pyna.networks.PythonNetwork(libraries=alpha_library)

And finally we can print an overview of the network as well as a Z-N plot of the nuclei linked via the reactions we selected.

[12]:

print(alpha_network.network_overview())

he4
  consumed by:
     c12 + he4 --> o16
     o16 + he4 --> ne20
     ne20 + he4 --> mg24
     mg24 + he4 --> si28
     si28 + he4 --> s32
     s32 + he4 --> ar36
     ar36 + he4 --> ca40
     ca40 + he4 --> ti44
     ti44 + he4 --> cr48
     cr48 + he4 --> fe52
     fe52 + he4 --> ni56
  produced by:
     o16 --> he4 + c12
     ne20 --> he4 + o16
     mg24 --> he4 + ne20
     si28 --> he4 + mg24
     s32 --> he4 + si28
     ar36 --> he4 + s32
     ca40 --> he4 + ar36
     ti44 --> he4 + ca40
     cr48 --> he4 + ti44
     fe52 --> he4 + cr48
     ni56 --> he4 + fe52

c12
  consumed by:
     c12 + he4 --> o16
  produced by:
     o16 --> he4 + c12

o16
  consumed by:
     o16 --> he4 + c12
     o16 + he4 --> ne20
  produced by:
     c12 + he4 --> o16
     ne20 --> he4 + o16

ne20
  consumed by:
     ne20 --> he4 + o16
     ne20 + he4 --> mg24
  produced by:
     o16 + he4 --> ne20
     mg24 --> he4 + ne20

mg24
  consumed by:
     mg24 --> he4 + ne20
     mg24 + he4 --> si28
  produced by:
     ne20 + he4 --> mg24
     si28 --> he4 + mg24

si28
  consumed by:
     si28 --> he4 + mg24
     si28 + he4 --> s32
  produced by:
     mg24 + he4 --> si28
     s32 --> he4 + si28

s32
  consumed by:
     s32 --> he4 + si28
     s32 + he4 --> ar36
  produced by:
     si28 + he4 --> s32
     ar36 --> he4 + s32

ar36
  consumed by:
     ar36 --> he4 + s32
     ar36 + he4 --> ca40
  produced by:
     s32 + he4 --> ar36
     ca40 --> he4 + ar36

ca40
  consumed by:
     ca40 --> he4 + ar36
     ca40 + he4 --> ti44
  produced by:
     ar36 + he4 --> ca40
     ti44 --> he4 + ca40

ti44
  consumed by:
     ti44 --> he4 + ca40
     ti44 + he4 --> cr48
  produced by:
     ca40 + he4 --> ti44
     cr48 --> he4 + ti44

cr48
  consumed by:
     cr48 --> he4 + ti44
     cr48 + he4 --> fe52
  produced by:
     ti44 + he4 --> cr48
     fe52 --> he4 + cr48

fe52
  consumed by:
     fe52 --> he4 + cr48
     fe52 + he4 --> ni56
  produced by:
     cr48 + he4 --> fe52
     ni56 --> he4 + fe52

ni56
  consumed by:
     ni56 --> he4 + fe52
  produced by:
     fe52 + he4 --> ni56

[13]:

alpha_network.plot()