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()

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()
