Reaction Network Calculator
   +                                                                         RNC:3.2 (Lynyrd Skynyrd version)
L.Fortunato*, G. Albertin, A.F. Lopez-Loaiza, M .Torresan
University of Padova & I.N.F.N. Padova - Italy
© All rights reserved. Distributed for study and non-profit research purpose only.
History :
v. 0.0: only proof of principle, originally written by L. Fortunato (2020)
v. 1.0: Jetter cycle reactions added by G. Albertin, constants revised (Nov 2022)
v. 2.0: Kinetic array and Rate array created. Equations rewritten. (May 2023)
v. 2.1: New reactions added. SpList introduced and automatized procedure for masses calc
v. 2.2: Post cycle reactions added by A.F Lopez-Loaiza.  Jetter-Post cycles merged. (Jul 2023) 
v. 2.3: Decays added; minor corrections and revisions; graphs improved; new data added by L. Fortunato. (Oct 2023)
v. 2.4: Stripped of unnecessary information, graphs improved. Energy distributions (Maxwell, Planck) added by M. Torresan, but not used in the present calc's. (Nov 2023)
v. 3.0: Bethe-Bloch Stopping power included in calculations. Eqns for Kinetic energy decrement added in additional file. Mistake on reaction n3 corrected. (Mar 2024)
v. 3.1: Extrapolations for cross-sections outside data are added (0 for E<Emin and  exponential for E>Emax). (Apr 2024)
v. 3.2: Conditions on positivity of NDSolve solutions added. (May2024)

(* Define the path to the directory where the files are stored  *)
path=NotebookDirectory[]
AppendTo[$Path,path ];
AppendTo[$Path,path <>"Data/"];

Cross Section Data Reading from external files in Data/ directory

(* Initial data reading, for each reaction of interest. Set the filename inside file list, taking in account that the standard units are mbarns for cross section and MeV for energy *)
filetoberead="Constants.nb"; (*Read file Constants*)
NotebookEvaluate[path<>"Data/"<>filetoberead,InsertResults->True];

(* filetoberead="Maxwell.nb";  (*Read file with Energy Ditributions*)
NotebookEvaluate[path<>"Data/"<>filetoberead,InsertResults->True]; *)

 filetoberead="Bethe-Bloch.nb";  (*Read file with Stoppin Power formulas *)
NotebookEvaluate[path<>"Data/"<>filetoberead,InsertResults->True]; 

files={"6Li(n,3H)4He.nb","3H(2H,n)4He.nb","2H(n,2n)p.nb","6Li(p,3He)4He.nb","3He(2H,p)4He.nb"};
nofi=Length[files]
For[ind=1,ind<=nofi,ind++,
filetoberead=files[[ind]];
NotebookEvaluate[path<>"Data/"<>filetoberead,InsertResults->True];
minE[ind]=data[[1,1]];  maxE[ind]=data[[-1,1]];
data0[ind]=data;
σ0[ind]=Interpolation[data0[ind]];
Bcoeff[ind_,En_]:= -(D[σ0[ind][dd],dd]/.dd->maxE[ind])/σ0[ind][maxE[ind]];
Acoeff[ind_,En_]:= σ0[ind][maxE[ind]]*Exp[Bcoeff[ind,En]*maxE[ind]];
σ[ind_][En_]:=If[En<minE[ind],0.,If[En>maxE[ind],Acoeff[ind,En]*Exp[-Bcoeff[ind,En]*En]  ,σ0[ind][En ]]]
]
(* Exponential extrapolation at high energy is rather arbitrary, it might be larger and misses any resonance *)

lplo[k_]:=LogLogPlot[σ[k][x],{x,0,maxE[k]},PlotLabel->files[[k]]]
Grid[{Table[lplo[k],{k,1,nofi}]}]
(*table with the list 1 and 2,.. and so on. As the number of reactions increase this list will grow*)
(*All cross-sections are set to 0. below the lowest datum, to avoid numerical problems later on, at most this is an underestimation for neutron cross-sections *)

Particle indexes:   1) n    2) 2H    3) 3H    4) 4He    5) 6Li    6)1H    7)3He ,  ... if new species are appended to the end of the list, the ordering is preserved in the equations below, no need for changes of indexes.

Use notation {Z,A} to conform to IsotopeData

SpList= {{0,1}  (* neutron *),  {1,2}   (* deuteron *),{1,3} (* tritium *),{2,4}  (* alpha *),  {3,6} (* 6Li *),{1,1}(* proton *), {2,3} (* 3He *)}  ;  

Do[
A[i]=SpList[[i,2]];
Z[i]=SpList[[i,1]];
M[i]=IsotopeData[SpList[[i]],"AtomicMass"][[1]]*amu-Z[i]*me
,{i,1,Length[SpList]}]

Decay constants (in Hz = 1/s)

(* Decay constants for radioactive species indicated with species index *) 
λ[1]= 1.13688*10^(-3)   (* neutron decay τ = 879.6 s *)
λ[3]= 3.63738*10^(-9)   (* tritium decay t1/2 = 12.32 y, remember 1 Julian year =365.25 days *)

Reactions  (parametrization of S-factors or cross-sections)

1)  6Li(n,T) 4He           2) D(T,n) 4He   warning: rate is written as if deuteron is moving       Jetter Cycle

3) n(D,2n)p

4) 6Li(p,3He) 4He     5) 3He(D,p) 4He         Post Cycle

(* Q-values of reactions in MeV *)
Q[1]:=4.7830; Q[2]:= 17.5892; Q[3]:=-2.2245; Q[4]:= 4.0197; Q[5]:= 18.353055 ; 

(*Velocity of incident particle , adimensional β=v/c  *) 
β[1][En_]:=Sqrt[2En/(M[1]) ]    (*Neutron Velocity*)
β[2][En_]:=Sqrt[2En/(M[2]) ]    (*Deuterium Velocity*) (*According to the cross section the deuteron is moving and the target is tritium *)
β[3][En_]:=Sqrt[2En/(M[1]) ]    (*Neutron Velocity*)
β[4][En_]:=Sqrt[2En/(M[6]) ]    (*Proton Velocity*)
β[5][En_]:=Sqrt[2En/(M[7]) ]    (*Helium-3 Velocity*)

(* These are not really rates, rather σβc are volumes over time σ*v products, while rate is mass/time *)
Do[σβc[i_][En_]:=NA*σ[k][En]*β[k][En]*c*10^(-31)/.k->i;
rate[i_][En_]:=densLiD * σβc[i][En]
,{i,1,nofi}
];
(*β in units of c, factor 10^(-31)to convert mb -> m^2, divided by Volmol to have mol/s*)

Reactions  rates

Do[plotσβc[k]=LogLogPlot[σβc[k][En] ,{En,minE[k],maxE[k]}, PlotRange->All, PlotStyle->{RandomColor[k]},PlotLabel->StringDrop[files[[k]],-3],AxesLabel->{"E [MeV]", "NA<σv> [m³/s mol]"},GridLines->Full] ,
{k,1,nofi}];
GraphicsGrid[{{plotσβc[1],plotσβc[2],plotσβc[3]},{plotσβc[4],plotσβc[5],Null}},ImageSize->Full]

(*I've found that the flux of neutron that a neutron beam facility can product is from 10^5 particle/cm^2/s to 10^14 particle/cm^2/s so ca 1 pmA of beam current *)

i = 10^(14)(* neutrons / cm² * s*)
i = i/(1.675*10^(-24)*NA)(* flow of neutrons mol/g*s *)

Energy  in MeV Array  Kin[ i ] [ j ] [time]     i: Species ;  j:Label for Energy  -> promoted to function of time, but still keep the neutrons from different reactions separated (different initial energy)

Rate Array [ i ] [ j ][time]      i :  Reaction index    ;     j : Label for Energy  also promoted to functions of time

Rate[1][1][t_]:= rate[1][Kin[1][1][t]]  (*Reaction 1 induced by neutron from the beam*)
Rate[1][2][t_]:= rate[1][Kin[1][2][t]]  (*Reaction 1 induced by neutron product in reaction 2*)

Rate[2][1][t_]:= rate[2][Kin[2][1][t]]   (*Reaction 2 induced by tritium product in reaction 1*)

Rate[3][1][t_]:= rate[3][Kin[1][1][t]]  (*Reaction 3 induced by neutron from the beam *)
Rate[3][2][t_]:= rate[3][Kin[1][2][t]]  (*Reaction 3 induced by neutron product in reaction 2*)

Rate[4][2][t_]:=rate[4][Kin[6][2][t]]  (*Reaction 4 induced by proton product in reaction 5*)

Rate[5][1][t_]:=rate[5][Kin[7][1][t]]  (* Reaction 5  induced by 3-He product in reaction 4*)

Set of equations & Boundary Conditions (these equations are in 1/seconds)

Jetter cycle + Post cycle

Distribution for endothermic 2H(n,2n)p

(*Read distribution file and do the interpolation*)
dis=Import["distribution_E_p.txt","Table"];  (*Reading data from Distribution.txt*)
P=Interpolation[dis];                                        (*Interpolation, units of 1/MeV *)
min=dis[[1,1]]
max=dis[[-1,1]]

(*Plot*)
plot1=Plot[P[x],{x,dis[[1,1]],dis[[-1,1]]},PlotRange->All,PlotStyle -> { Green,Thickness[0.010]}, PlotLegends-> Placed[{"Interpolation"},Right]];
plot2=ListPlot[dis,Filling->Axis,PlotStyle -> {Black},PlotLegends->Placed[{"Simulated Data"},Right]] ;
Show[plot1,plot2,AxesLabel->{"E [MeV]", "Probability"} ,ImageSize->500,LabelStyle->Directive[FontSize->Scaled[.03]],TicksStyle->Directive[FontSize->Scaled[.08]],AxesStyle->Black]
(* Integration: normalization check *)
NIntegrate[P[x],{x,dis[[1,1]],dis[[-1,1]]}]
dumkin=NIntegrate[P[x]x,{x,dis[[1,1]],dis[[-1,1]]}]

LogLinearPlot[σ[4][x],{x,min,max},AxesLabel->{"E [MeV]", "σ [mb]"} ,PlotStyle -> { Green,Thickness[0.005]},ImageSize->500,LabelStyle->Directive[FontSize->Scaled[.03]],TicksStyle->Directive[FontSize->Scaled[.08]],GridLines->Full,ImageSize->500,LabelStyle->Directive[FontSize->Scaled[.03]],TicksStyle->Directive[FontSize->Scaled[.08]],AxesStyle->Black]

(*v sigma just mediated on probability P *)
(* μ[En_]:=β[4][En]*σ[4][En]*P[En]*c*10^(-31)*NA  This was the ond one, next is simpler *)
μ[En_]:=σβc[4][En]*P[En]

LogLinearPlot[μ[x],{x,min,max},AxesLabel->{"E [MeV]", "μ(E) [m³/(mol*s*MeV)]"} ,PlotStyle -> { Blue,Thickness[0.005]},ImageSize->500,LabelStyle->Directive[FontSize->Scaled[.03]],TicksStyle->Directive[FontSize->Scaled[.08]],GridLines->Full,ImageSize->500,LabelStyle->Directive[FontSize->Scaled[.03]],TicksStyle->Directive[FontSize->Scaled[.08]],Filling->Axis,AxesStyle->Black]
NIntegrate[μ[x],{x,min,max}]  
Rate[4][1][t_]:=densLiD*μ[ Kin[6][1][t]] (*Reaction 4 induced by proton product in reaction 3 D-breakup*)

(* Index 0 means sink terms *)
Equ[1]:=D[n[1][1][t],t]==i-n[1][1][t]*n[5][1][t]*Rate[1][1][t] -λ[1]*n[1][1][t]
Equ[2]:=D[n[1][2][t],t]==n[2][1][t]*n[3][1][t]*Rate[2][1][t]-n[1][2][t]*n[5][1][t]*Rate[1][2][t] -n[1][2][t]*n[2][1][t]*Rate[3][2][t]-λ[1]*n[1][2][t]
Equ[3]:=D[n[1][3][t],t]==n[1][1][t]*n[2][1][t]*Rate[3][1][t]-λ[1]*n[1][3][t]     (* It has been corrected - L.F. Mar.2024 *)
Equ[4]:=D[n[1][4][t],t]==n[1][2][t]*n[2][1][t]*Rate[3][2][t]-λ[1]*n[1][4][t]
Equ[5]:=D[n[4][1][t],t]==n[1][1][t]*n[5][1][t]*Rate[1][1][t]+n[1][2][t]*n[5][1][t]*Rate[1][2][t]
Equ[6]:=D[n[4][2][t],t]==n[2][1][t]*n[3][1][t]*Rate[2][1][t]
Equ[7]:=D[n[4][3][t],t]==n[6][1][t]*n[5][1][t]*Rate[4][1][t]+n[6][2][t]*n[5][1][t]*Rate[4][2][t]
Equ[8]:=D[n[4][4][t],t]==n[2][1][t]*n[7][1][t]*Rate[5][1][t]
Equ[9]:=D[n[2][1][t],t]==-n[2][1][t]*n[3][1][t]*Rate[2][1][t]-n[2][1][t]*n[7][1][t]*Rate[5][1][t]-n[2][1][t]*n[1][2][t]*Rate[3][2][t]
Equ[10]:=D[n[3][1][t],t]==n[1][1][t]*n[5][1][t]*Rate[1][1][t]-n[2][1][t]*n[3][1][t]*Rate[2][1][t]+n[1][2][t]*n[5][1][t]*Rate[1][2][t]-λ[3]*n[3][1][t]
Equ[11]:=D[n[7][1][t],t]==n[6][2][t]*n[5][1][t]*Rate[4][2][t]+n[6][1][t]*n[5][1][t]*Rate[4][1][t]-n[2][1][t]*n[7][1][t]*Rate[5][1][t]
Equ[12]:=D[n[6][1][t],t]==n[2][1][t]*n[1][2][t]*Rate[3][2][t]-n[6][1][t]*n[5][1][t]*Rate[4][1][t]
Equ[13]:=D[n[6][2][t],t]==n[2][1][t]*n[7][1][t]*Rate[5][1][t]-n[6][2][t]*n[5][1][t]*Rate[4][2][t]
Equ[14]:=D[n[5][1][t],t]==-n[1][1][t]*n[5][1][t]*Rate[1][1][t]-n[1][2][t]*n[5][1][t]*Rate[1][2][t]-n[6][1][t]*n[5][1][t]*Rate[4][1][t]-n[6][2][t]*n[5][1][t]*Rate[4][2][t]
Equ[15]:= D[n[6][0][t],t]==λ[1]*n[1][1][t]+λ[1]*n[1][2][t]+λ[1]*n[1][3][t]+λ[1]*n[1][4][t]
Equ[16]:= D[n[7][0][t],t]==λ[3]*n[3][1][t]


BouCond:={n[1][1][0]==0.0,n[1][2][0]==0.0,n[1][3][0]==0.0,n[1][4][0]==0.0, n[2][1][0]==0.125, n[3][1][0]==0.0, n[4][1][0]==0.0,n[4][2][0]==0.0,n[4][3][0]==0.0,n[4][4][0]==0.0,n[5][1][0]==0.125,n[6][0][0]==0.0,n[6][1][0]==0.0,n[6][2][0]==0.0 ,n[7][1][0]==0.0, n[7][0][0]==0.0}

varia :={n[1][1][t],n[1][2][t],n[1][3][t],n[1][4][t] ,n[2][1][t], n[3][1][t], n[4][1][t],n[4][2][t],n[4][3][t],n[4][4][t],n[5][1][t],n[6][0][t] ,n[6][1][t],n[6][2][t] ,n[7][1][t],n[7][0][t]}

(* Slowing by Ionization of all charged particles - neutron deceleration not implemented *)
Slow[1]:= D[Kin[1][1][t],t]== 0.   (* neutrons do not decelerate *)
Slow[2]:= D[Kin[1][2][t],t]== 0.    (* neutrons do not decelerate *)
Slow[3]:= D[Kin[2][1][t],t] == -dEdt[Z[2],M[2],Kin[2][1][t]]   (* tritium from n-6Li *) 
Slow[4]:= D[Kin[4][1][t],t] == -dEdt[Z[4],M[4],Kin[4][1][t]]  (* alphas from n-6Li *)
Slow[5]:= D[Kin[4][2][t],t] == -dEdt[Z[4],M[4],Kin[4][2][t]]   (* alphas from D-T *)
Slow[6]:= D[Kin[4][3][t],t] == -dEdt[Z[4],M[4],Kin[4][3][t]]    (* alphas from p-6Li *)
Slow[7]:= D[Kin[4][4][t],t] == -dEdt[Z[4],M[4],Kin[4][4][t]]    (* alphas from d-3He *)
Slow[8]:= D[Kin[6][1][t],t] == -dEdt[Z[6],M[6],Kin[6][1][t]]    (* proton from d-breakup *)
Slow[9]:= D[Kin[6][2][t],t] == -dEdt[Z[6],M[6],Kin[6][2][t]]   (* proton from d-3He *)
Slow[10]:= D[Kin[7][1][t],t] == -dEdt[Z[7],M[7],Kin[7][1][t]]    (* helion from d-3He *)

newvaria:={Kin[1][1][t],Kin[1][2][t],Kin[2][1][t],Kin[4][1][t],Kin[4][2][t],Kin[4][3][t],Kin[4][4][t],Kin[6][1][t],Kin[6][2][t],Kin[7][1][t]}

KinCond:={Kin[1][1][0]==.24,Kin[1][2][0]== 14.0481,Kin[2][1][0]== 2.7273 * M[3]/M[2],Kin[4][1][0]==2.0557,Kin[4][2][0]==3.5411,Kin[4][3][0]==1.7188,Kin[4][4][0]==3.6908,Kin[6][1][0] == dumkin,Kin[6][2][0]== 14.6622,Kin[7][1][0]==2.2812};

abundancies={Equ[1],Equ[2],Equ[3],Equ[4],Equ[5],Equ[6],Equ[7],Equ[8],Equ[9],Equ[10],Equ[11],Equ[12],Equ[13],Equ[14],Equ[15],Equ[16],BouCond}; 
stopping={Slow[1],Slow[2],Slow[3],Slow[4],Slow[5],Slow[6],Slow[7],Slow[8],Slow[9],Slow[10], KinCond};
additionalconds={WhenEvent[n[1][1][t]<0,n[1][1][t]->0],WhenEvent[n[1][2][t]<0,n[1][2][t]->0],WhenEvent[n[1][3][t]<0,n[1][3][t]->0],WhenEvent[n[1][4][t]<0,n[1][4][t]->0],WhenEvent[n[2][1][t]<0,n[2][1][t]->0],WhenEvent[n[3][1][t]<0,n[3][1][t]->0],WhenEvent[n[4][1][t]<0,n[4][1][t]->0],WhenEvent[n[4][2][t]<0,n[4][2][t]->0],WhenEvent[n[4][3][t]<0,n[4][3][t]->0],WhenEvent[n[4][4][t]<0,n[4][4][t]->0],WhenEvent[n[5][1][t]<0,n[5][1][t]->0],WhenEvent[n[6][0][t]<0,n[6][0][t]->0],WhenEvent[n[6][1][t]<0,n[6][1][t]->0],WhenEvent[n[6][2][t]<0,n[6][2][t]->0],
WhenEvent[n[7][0][t]<0,n[7][0][t]->0],WhenEvent[n[7][1][t]<0,n[7][1][t]->0],
WhenEvent[Kin[1][1][t]<0,Kin[1][1][t]->0], WhenEvent[Kin[1][2][t]<0,Kin[1][2][t]->0],WhenEvent[Kin[2][1][t]<0,Kin[2][1][t]->0],WhenEvent[Kin[4][1][t]<0,Kin[4][1][t]->0],WhenEvent[Kin[4][2][t]<0,Kin[4][2][t]->0],WhenEvent[Kin[4][3][t]<0,Kin[4][3][t]->0],WhenEvent[Kin[4][4][t]<0,Kin[4][4][t]->0],WhenEvent[Kin[6][1][t]<0,Kin[6][1][t]->0],WhenEvent[Kin[6][2][t]<0,Kin[6][2][t]->0],WhenEvent[Kin[7][1][t]<0,Kin[7][1][t]->0]};
setofequations=Join[abundancies,stopping,additionalconds];
variables=Join[varia,newvaria];

tfinal=10*10^(-10);
sol=NDSolve[setofequations,variables,{t,0,tfinal}]

LogLinearPlot[Evaluate[{variables[[19]]}/.sol],{t,0,tfinal},PlotRange->All]

tke=Table[Evaluate[variables[[kz]]/.sol][[1]],{kz,17,26}];
LogLinearPlot[tke,{t,0,tfinal},PlotRange->{0,16}]

Y11[t_]=n[1][1][t]/. sol;
Y12[t_]=n[1][2][t]/. sol;
Y13[t_]=n[1][3][t]/. sol;
Y14[t_]=n[1][4][t]/. sol; 
Y21[t_]=n[2][1][t]/. sol;
Y31[t_]=n[3][1][t]/. sol;
Y41[t_]=n[4][1][t]/. sol;
Y42[t_]=n[4][2][t]/. sol;
Y43[t_]=n[4][3][t]/. sol;
Y44[t_]=n[4][4][t]/. sol;
Y51[t_]=n[5][1][t]/. sol;
Y60[t_]=n[6][0][t]/. sol;
Y61[t_]=n[6][1][t]/. sol;
Y62[t_]=n[6][2][t]/. sol;
Y70[t_]=n[7][0][t]/. sol;
Y71[t_]=n[7][1][t]/. sol; 

FinAbu={Y11[tfinal][[1]],Y12[tfinal][[1]],Y13[tfinal][[1]],Y14[tfinal][[1]],Y21[tfinal][[1]],Y31[tfinal][[1]],Y41[tfinal][[1]],Y42[tfinal][[1]],Y43[tfinal][[1]],Y44[tfinal][[1]],Y51[tfinal][[1]],Y60[tfinal][[1]],Y61[tfinal][[1]],Y62[tfinal][[1]],Y71[tfinal][[1]],Y70[tfinal][[1]]};
namz={"n_1","n_2","n_3","n_4","D","T","α_1","α_2","α_3","α_4","6Li","p_0","p_1","p_2","3He","3He-0"} ;
colorz={Purple,Brown,Gray,Black,Blue,Green,Yellow,Orange,Yellow,Orange,Magenta,Red,Red,Green,White,Pink,Purple};

BarChart[Log[FinAbu],ChartLabels->namz[[1;;-1]],ChartStyle->colorz[[1;;-1]],GridLines->All,AxesStyle->Black, PlotLabel->"Final"]

BarChart[Table[BouCond[[i,2]],{i,1,Length[BouCond]}],ChartLabels->namz[[2;;-1]],ChartStyle->colorz[[2;;-1]],GridLines->All,AxesStyle->Black, PlotLabel->"Initial"]

(* Result of including decays of the neutron -> negligible *)
Y60[tfinal][[1]]
Y70[tfinal][[1]]

Do[
evol[k]=LogLinearPlot[Evaluate[{variables[[k]]}/.sol],{t,10^(-20),tfinal},PlotStyle->{colorz[[k]]},PlotRange->All,PlotLegends-> Placed[ { "Y_"<>namz[[k]]<>"(t)"},Right], AxesLabel->{"t [s]", "Y(t) [mol/g]"},TicksStyle->Directive[FontSize->Scaled[.015]],GridLines->Full,AxesStyle->Black],
{k,1,Length[varia]}]

evol[1]
evol[2]
evol[3]
evol[4]

(* GraphicsGrid[Table[Table[evol[k+3j],{k,1,3}],{j,0,3}], ImageSize->Full] *)

LogLinearPlot[Evaluate[{n[1][2][t],n[1][3][t],n[1][4][t], n[2][1][t],n[3][1][t],n[4][1][t],n[4][2][t],n[4][3][t],n[4][4][t],n[5][1][t],n[6][1][t],n[6][2][t], n[7][1][t]}/.sol],{t,0,tfinal},PlotStyle->{ Brown,ColorData[4,"ColorList"],Black,Blue,Cyan,Yellow,ColorData[44,"ColorList"], Orange,ColorData[23,"ColorList"],Magenta,Red,ColorData[35,"ColorList"],Green}, PlotRange->{0,0.125},PlotLegends-> Placed[ { Subscript[Y,n2],Subscript[Y,n3],Subscript[Y,n4],Subscript[Y,D],Subscript[Y,T],Subscript[Y,α1],Subscript[Y,α2],Subscript[Y,α3],Subscript[Y,α4],Subscript[Y,6Li],Subscript[Y,p1],Subscript[Y,p2],Subscript[Y,3He]},Right], AxesLabel->{"t [s]", "Y(t) [mol/g]"},ImageSize->650,LabelStyle->Directive[FontSize->Scaled[0.02],Background->White],GridLines->Full,AxesStyle->Black,PlotLabel->Style["Molar fractions",15,Black],Background->White]

LogLogPlot[Evaluate[{n[1][2][t],n[1][3][t],n[1][4][t], n[2][1][t],n[3][1][t],n[4][1][t],n[4][2][t],n[4][3][t],n[4][4][t],n[5][1][t],n[6][1][t],n[6][2][t], n[7][1][t]}/.sol],{t,0,tfinal},PlotStyle->{ Brown,ColorData[4,"ColorList"],Black,Blue,Cyan,Yellow,ColorData[44,"ColorList"], Orange,ColorData[23,"ColorList"],Magenta,Red,ColorData[35,"ColorList"],Green}, PlotRange->{10^(-23),0.125},PlotLegends-> Placed[ { "Y_n_2(t)","Y_n_3(t)","Y_n_4(t)","Y_D(t)","Y_T(t)","Y_α_1(t)","Y_α_2(t)","Y_α_3(t)","Y_α_4(t)","Y_6Li(t)","Y_p_1(t)","Y_p_2(t)","Y_3He(t)"},Right], AxesLabel->{"t [s]", "Y(t) [mol/g]"},ImageSize->650,TicksStyle->Directive[FontSize->Scaled[.015]],GridLines->Full,AxesStyle->Black]

Conversion into alpha’s

Y41[tfinal][[1]] + Y42[tfinal][[1]] + Y43[tfinal][[1]] + Y44[tfinal][[1]]

(*P: power produced per unit volume [J/cm^3*s]]*)
(*pij:  i-> reaction , j-> label energy*)
factor=9.06*ech*10^(-9); 

(*Reaction 1 *)
p11=factor*densLiD*NA*Rate[1][1][t]*Q[1] ; (*beam neutrons*)
p12=factor*densLiD*NA*Rate[1][2][t]*Q[1];  (*fusion product neutrons*)

(*Reaction 2: Fusion.*)
p21=factor*densLiD*NA*Rate[2][1][t]*Q[2]; (*fission product tritium*)

(*Reaction 3: Deuterium-breakup.*)
p31=factor*densLiD*NA*Rate[3][2][t]*Q[3]; 

(*Reaction 4: Fission.*)
p41=factor*densLiD*NA*Rate[4][1][t]*Q[4];  (*protons from fusion product*)
p42=factor*densLiD*NA*Rate[4][2][t]*Q[4];  (*protons from fusion product*)

(*Reaction 5: Fusion.*)
p51=factor*densLiD*NA*Rate[5][1][t]*Q[5]; (*fission product tritium*)

(*Together*)
PowerJP={NIntegrate[Evaluate[{p11*n[1][1][t]*n[5][1][t]}/.sol],{t,0,tfinal}][[1]][[1]],NIntegrate[Evaluate[{p12*n[1][2][t]*n[5][1][t]}/.sol],{t,0,tfinal}][[1]][[1]],NIntegrate[Evaluate[{p21*n[2][1][t]*n[3][1][t]}/.sol],{t,0,tfinal}][[1]][[1]],NIntegrate[Evaluate[{p41*n[6][1][t]*n[5][1][t]}/.sol],{t,0,tfinal}][[1]][[1]],NIntegrate[Evaluate[{p42*n[6][2][t]*n[5][1][t]}/.sol],{t,0,tfinal}][[1]][[1]],NIntegrate[Evaluate[{p51*n[2][1][t]*n[7][1][t]}/.sol],{t,0,tfinal}][[1]][[1]],NIntegrate[Evaluate[{p31*n[2][1][t]*n[1][2][t]}/.sol],{t,0,tfinal}][[1]][[1]]}
Print["Integrated Energy: ",  Total[PowerJP], " GJ"]

LogLinearPlot[{
Callout[Evaluate[p11*n[1][1][t]*n[5][1][t]/.sol],"~"NumberForm[PowerJP[[1]],{4,1}]"GJ",Scaled[0.3],Appearance->Frame],Callout[Evaluate[p12*n[1][2][t]*n[5][1][t]/.sol],"~"NumberForm[PowerJP[[2]],{4,1}]"GJ",Scaled[0.4],Appearance->Frame],Callout[ Evaluate[p21*n[2][1][t]*n[3][1][t]/.sol],"~"NumberForm[PowerJP[[3]],{4,1}]"GJ",Scaled[0.5],Appearance->Frame], Callout[Evaluate[p41*n[6][1][t]*n[5][1][t]/.sol],"~"NumberForm[PowerJP[[4]],{4,1}]"GJ",Scaled[0.6],Appearance->Frame], Callout[Evaluate[p42*n[6][2][t]*n[5][1][t]/.sol],"~"NumberForm[PowerJP[[5]],{4,1}]"GJ",Scaled[0.7],Appearance->Frame], Callout[Evaluate[p51*n[2][1][t]*n[7][1][t]/.sol],"~"NumberForm[PowerJP[[6]],{4,1}]"GJ",Scaled[0.8],Appearance->Frame], Callout[Evaluate[p31*n[2][1][t]*n[1][2][t]/.sol],"~"NumberForm[PowerJP[[7]],{4,1}]"GJ",Below,Scaled[0.9],Appearance->Frame]},{t,0,tfinal},PlotStyle -> { Black,Black,Red,Black,Black,Red,Blue,Thickness[0.003]}, PlotRange->{All, All}, Filling->Axis,AxesLabel->{Style["t [s]",14], Style["P(t) [ GJ / s ]",14]},ImageSize->650,TicksStyle->Directive[FontSize->Scaled[.015]],GridLines->Full,AxesStyle->Black]

Beep[]