@@ -4778,11 +4778,13 @@ def _make_rule(rr):
47784778 if n > 1 :
47794779 kin = arr ().fit_to_reactions (rs , recipe = recipe )
47804780 if n == 1 or kin .E0 .value_si < 0.0 :
4781+ # still run it through the averaging function when n=1 to standardize the units and run checks
47814782 kin = average_kinetics ([r .kinetics for r in rs ])
4782- #kin.comment = "Only one reaction or Arrhenius BM fit bad. Instead averaged from {} reactions.".format(n)
47834783 if n == 1 :
47844784 kin .uncertainty = RateUncertainty (mu = 0.0 , var = (np .log (fmax ) / 2.0 ) ** 2 , N = 1 , Tref = Tref , data_mean = data_mean , correlation = label )
4785+ kin .comment = f"Only one reaction rate: { rs [0 ]!s} "
47854786 else :
4787+ kin .comment = f"Blowers-Masel fit was bad (E0<0) so instead averaged from { n } reactions."
47864788 dlnks = np .array ([
47874789 np .log (
47884790 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 )
@@ -4796,34 +4798,31 @@ def _make_rule(rr):
47964798 mu = np .dot (ws , dlnks ) / V1
47974799 s = np .sqrt (np .dot (ws , (dlnks - mu ) ** 2 ) / (V1 - V2 / V1 ))
47984800 kin .uncertainty = RateUncertainty (mu = mu , var = s ** 2 , N = n , Tref = Tref , data_mean = data_mean , correlation = label )
4799- else :
4800- if n == 1 :
4801- kin .uncertainty = RateUncertainty (mu = 0.0 , var = (np .log (fmax ) / 2.0 ) ** 2 , N = 1 , Tref = Tref , data_mean = data_mean , correlation = label )
4802- else :
4803- if isinstance (rs [0 ].kinetics , Arrhenius ):
4804- dlnks = np .array ([
4805- np .log (
4806- arr ().fit_to_reactions (rs [list (set (range (len (rs ))) - {i })], recipe = recipe )
4807- .to_arrhenius (rxn .get_enthalpy_of_reaction (Tref ))
4808- .get_rate_coefficient (T = Tref ) / rxn .get_rate_coefficient (T = Tref )
4809- ) for i , rxn in enumerate (rs )
4810- ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref
4811- else :
4812- dlnks = np .array ([
4813- np .log (
4814- arr ().fit_to_reactions (rs [list (set (range (len (rs ))) - {i })], recipe = recipe )
4815- .to_arrhenius_charge_transfer (rxn .get_enthalpy_of_reaction (Tref ))
4816- .get_rate_coefficient (T = Tref ) / rxn .get_rate_coefficient (T = Tref )
4817- ) for i , rxn in enumerate (rs )
4818- ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref
4819- varis = (np .array ([rank_accuracy_map [rxn .rank ].value_si for rxn in rs ]) / (2.0 * 8.314 * Tref )) ** 2
4820- # weighted average calculations
4821- ws = 1.0 / varis
4822- V1 = ws .sum ()
4823- V2 = (ws ** 2 ).sum ()
4824- mu = np .dot (ws , dlnks ) / V1
4825- s = np .sqrt (np .dot (ws , (dlnks - mu ) ** 2 ) / (V1 - V2 / V1 ))
4826- kin .uncertainty = RateUncertainty (mu = mu , var = s ** 2 , N = n , Tref = Tref , data_mean = data_mean , correlation = label )
4801+ else : # Blowers-Masel fit was good
4802+ if isinstance (rs [0 ].kinetics , Arrhenius ):
4803+ dlnks = np .array ([
4804+ np .log (
4805+ arr ().fit_to_reactions (rs [list (set (range (len (rs ))) - {i })], recipe = recipe )
4806+ .to_arrhenius (rxn .get_enthalpy_of_reaction (Tref ))
4807+ .get_rate_coefficient (T = Tref ) / rxn .get_rate_coefficient (T = Tref )
4808+ ) for i , rxn in enumerate (rs )
4809+ ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref
4810+ else : # SurfaceChargeTransfer or ArrheniusChargeTransfer
4811+ dlnks = np .array ([
4812+ np .log (
4813+ arr ().fit_to_reactions (rs [list (set (range (len (rs ))) - {i })], recipe = recipe )
4814+ .to_arrhenius_charge_transfer (rxn .get_enthalpy_of_reaction (Tref ))
4815+ .get_rate_coefficient (T = Tref ) / rxn .get_rate_coefficient (T = Tref )
4816+ ) for i , rxn in enumerate (rs )
4817+ ]) # 1) fit to set of reactions without the current reaction (k) 2) compute log(kfit/kactual) at Tref
4818+ varis = (np .array ([rank_accuracy_map [rxn .rank ].value_si for rxn in rs ]) / (2.0 * 8.314 * Tref )) ** 2
4819+ # weighted average calculations
4820+ ws = 1.0 / varis
4821+ V1 = ws .sum ()
4822+ V2 = (ws ** 2 ).sum ()
4823+ mu = np .dot (ws , dlnks ) / V1
4824+ s = np .sqrt (np .dot (ws , (dlnks - mu ) ** 2 ) / (V1 - V2 / V1 ))
4825+ kin .uncertainty = RateUncertainty (mu = mu , var = s ** 2 , N = n , Tref = Tref , data_mean = data_mean , correlation = label )
48274826
48284827 #site solute parameters
48294828 site_datas = [get_site_solute_data (rxn ) for rxn in rxns ]
0 commit comments