@@ -4571,11 +4571,13 @@ def _make_rule(rr):
45714571 if n > 1 :
45724572 kin = arr ().fit_to_reactions (rs , recipe = recipe )
45734573 if n == 1 or kin .E0 .value_si < 0.0 :
4574+ # still run it through the averaging function when n=1 to standardize the units and run checks
45744575 kin = average_kinetics ([r .kinetics for r in rs ])
4575- #kin.comment = "Only one reaction or Arrhenius BM fit bad. Instead averaged from {} reactions.".format(n)
45764576 if n == 1 :
45774577 kin .uncertainty = RateUncertainty (mu = 0.0 , var = (np .log (fmax ) / 2.0 ) ** 2 , N = 1 , Tref = Tref , data_mean = data_mean , correlation = label )
4578+ kin .comment = f"Only one reaction rate: { rs [0 ]!s} "
45784579 else :
4580+ kin .comment = f"Blowers-Masel fit was bad (E0<0) so instead averaged from { n } reactions."
45794581 dlnks = np .array ([
45804582 np .log (
45814583 average_kinetics ([r .kinetics for r in rs [list (set (range (len (rs ))) - {i })]]).get_rate_coefficient (T = Tref ) / rxn .get_rate_coefficient (T = Tref )
@@ -4589,34 +4591,31 @@ def _make_rule(rr):
45894591 mu = np .dot (ws , dlnks ) / V1
45904592 s = np .sqrt (np .dot (ws , (dlnks - mu ) ** 2 ) / (V1 - V2 / V1 ))
45914593 kin .uncertainty = RateUncertainty (mu = mu , var = s ** 2 , N = n , Tref = Tref , data_mean = data_mean , correlation = label )
4592- else :
4593- if n == 1 :
4594- kin .uncertainty = RateUncertainty (mu = 0.0 , var = (np .log (fmax ) / 2.0 ) ** 2 , N = 1 , Tref = Tref , data_mean = data_mean , correlation = label )
4595- else :
4596- if isinstance (rs [0 ].kinetics , Arrhenius ):
4597- dlnks = np .array ([
4598- np .log (
4599- arr ().fit_to_reactions (rs [list (set (range (len (rs ))) - {i })], recipe = recipe )
4600- .to_arrhenius (rxn .get_enthalpy_of_reaction (Tref ))
4601- .get_rate_coefficient (T = Tref ) / rxn .get_rate_coefficient (T = Tref )
4602- ) for i , rxn in enumerate (rs )
4603- ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref
4604- else :
4605- dlnks = np .array ([
4606- np .log (
4607- arr ().fit_to_reactions (rs [list (set (range (len (rs ))) - {i })], recipe = recipe )
4608- .to_arrhenius_charge_transfer (rxn .get_enthalpy_of_reaction (Tref ))
4609- .get_rate_coefficient (T = Tref ) / rxn .get_rate_coefficient (T = Tref )
4610- ) for i , rxn in enumerate (rs )
4611- ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref
4612- varis = (np .array ([rank_accuracy_map [rxn .rank ].value_si for rxn in rs ]) / (2.0 * 8.314 * Tref )) ** 2
4613- # weighted average calculations
4614- ws = 1.0 / varis
4615- V1 = ws .sum ()
4616- V2 = (ws ** 2 ).sum ()
4617- mu = np .dot (ws , dlnks ) / V1
4618- s = np .sqrt (np .dot (ws , (dlnks - mu ) ** 2 ) / (V1 - V2 / V1 ))
4619- kin .uncertainty = RateUncertainty (mu = mu , var = s ** 2 , N = n , Tref = Tref , data_mean = data_mean , correlation = label )
4594+ else : # Blowers-Masel fit was good
4595+ if isinstance (rs [0 ].kinetics , Arrhenius ):
4596+ dlnks = np .array ([
4597+ np .log (
4598+ arr ().fit_to_reactions (rs [list (set (range (len (rs ))) - {i })], recipe = recipe )
4599+ .to_arrhenius (rxn .get_enthalpy_of_reaction (Tref ))
4600+ .get_rate_coefficient (T = Tref ) / rxn .get_rate_coefficient (T = Tref )
4601+ ) for i , rxn in enumerate (rs )
4602+ ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref
4603+ else : # SurfaceChargeTransfer or ArrheniusChargeTransfer
4604+ dlnks = np .array ([
4605+ np .log (
4606+ arr ().fit_to_reactions (rs [list (set (range (len (rs ))) - {i })], recipe = recipe )
4607+ .to_arrhenius_charge_transfer (rxn .get_enthalpy_of_reaction (Tref ))
4608+ .get_rate_coefficient (T = Tref ) / rxn .get_rate_coefficient (T = Tref )
4609+ ) for i , rxn in enumerate (rs )
4610+ ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref
4611+ varis = (np .array ([rank_accuracy_map [rxn .rank ].value_si for rxn in rs ]) / (2.0 * 8.314 * Tref )) ** 2
4612+ # weighted average calculations
4613+ ws = 1.0 / varis
4614+ V1 = ws .sum ()
4615+ V2 = (ws ** 2 ).sum ()
4616+ mu = np .dot (ws , dlnks ) / V1
4617+ s = np .sqrt (np .dot (ws , (dlnks - mu ) ** 2 ) / (V1 - V2 / V1 ))
4618+ kin .uncertainty = RateUncertainty (mu = mu , var = s ** 2 , N = n , Tref = Tref , data_mean = data_mean , correlation = label )
46204619
46214620 #site solute parameters
46224621 site_datas = [get_site_solute_data (rxn ) for rxn in rxns ]
0 commit comments