From 788ea2d3951601f788a4d2255cd80f334cab4d6f Mon Sep 17 00:00:00 2001 From: Luna Pratali Maffei Date: Thu, 14 Nov 2024 15:59:33 +0100 Subject: [PATCH 1/4] alternative fits if default fails --- ratefit/fit/arr.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/ratefit/fit/arr.py b/ratefit/fit/arr.py index ac60ede..12d5ff2 100644 --- a/ratefit/fit/arr.py +++ b/ratefit/fit/arr.py @@ -204,10 +204,20 @@ def fit_doub_arr(temps, kts, sing_params, a_change, n_change, tref=1.0, # Perform a least-squares fit # note: previous version scipy.optimize.leastsq (unbounded): used method='lm' # same or better results obtained with x_scale='jac' for new cases tested - plsq = least_squares(_resid_func, init_guess, bounds=bounds, - args=(temps, kts, doub_tref), x_scale = 'jac', #method = 'lm', - ftol=1.0E-8, xtol=1.0E-8, max_nfev=100000) - + try: + plsq = least_squares(_resid_func, init_guess, bounds=bounds, + args=(temps, kts, doub_tref), x_scale = 'jac', #method = 'lm', + ftol=1.0E-8, xtol=1.0E-8, max_nfev=100000) + except ValueError: #test older method without bounds + try: + plsq = least_squares(_resid_func, init_guess, + args=(temps, kts, doub_tref), method='lm', + ftol=1.0E-8, xtol=1.0E-8, max_nfev=100000) + except ValueError: + plsq = least_squares(_resid_func, init_guess, + loss = 'arctan', + args=(temps, kts, doub_tref), + ftol=1.0E-8, xtol=1.0E-8, max_nfev=100000) # Retrieve the fit params and convert A back to the input tref raw_params = list(plsq.x) # list of length 6 raw_params[0] = raw_params[0] * (tref / doub_tref) ** raw_params[1] From e346f27e80717f85f7373f5eb5d72cde783acbaf Mon Sep 17 00:00:00 2001 From: Luna Pratali Maffei Date: Thu, 14 Nov 2024 16:00:36 +0100 Subject: [PATCH 2/4] fix test after new mess_io fix --- mechanalyzer/tests/test__prompt.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/mechanalyzer/tests/test__prompt.py b/mechanalyzer/tests/test__prompt.py index 2c5c74d..e3a1954 100644 --- a/mechanalyzer/tests/test__prompt.py +++ b/mechanalyzer/tests/test__prompt.py @@ -166,8 +166,8 @@ def test_rovib_dos(): ped_1800 = ped_df_frag1_dct[1.0][1800] ped_2000 = ped_df_frag1_dct[1.0][2000] - assert np.isclose((ped_1800.iloc[166]), 0.0173, atol=1e-3, rtol=1e-2) - assert np.isclose((ped_2000.iloc[166]), 0.0148, atol=1e-3, rtol=1e-2) + assert np.isclose((ped_1800.iloc[166]), 0.01985, atol=1e-3, rtol=1e-2) + assert np.isclose((ped_2000.iloc[166]), 0.01724, atol=1e-3, rtol=1e-2) assert np.isclose(np.trapz(ped_1800.values, x=ped_1800.index), 1) assert np.isclose(np.trapz(ped_2000.values, x=ped_2000.index), 1) @@ -230,7 +230,7 @@ def test_bf_from_fne(): _, _, _, _, fne_bf = _read_data() bf_tp_dct = bf.bf_tp_dct( 'fne', None, None, 0.1, fne=fne_bf['CH3CHCH3']) - print(bf_tp_dct) + assert np.allclose( bf_tp_dct['CH3CHCH3'][1.0][1][[0, 6, 14]], np.array([1., 0.99999987, 0.99798199]), atol=1e-3, rtol=1e-2) @@ -320,5 +320,5 @@ def _read_data(): test_thermal() test_bf_from_phi1a() test_bf_from_fne() - test_rovib_dos() test_new_ktp_dct() + test_rovib_dos() From 3cefbe6e5be18005f1b1a9d7fda46f3183448a81 Mon Sep 17 00:00:00 2001 From: Luna Pratali Maffei Date: Thu, 14 Nov 2024 16:01:43 +0100 Subject: [PATCH 3/4] BF reading for hoten fixed if -stable- spc not found --- mechanalyzer/calculator/bf.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/mechanalyzer/calculator/bf.py b/mechanalyzer/calculator/bf.py index 11d77bd..26a7202 100644 --- a/mechanalyzer/calculator/bf.py +++ b/mechanalyzer/calculator/bf.py @@ -87,8 +87,11 @@ def rename_ktp_dct(ktp_dct, label, hotsp_dct): renamed_ktp_dct = {} spc0 = list(set(hotsp_dct.keys()).intersection(label[1]))[0] - renamed_ktp_dct[label] = ktp_dct[spc0] - ktp_dct.pop(spc0) + if spc0 in ktp_dct.keys(): + renamed_ktp_dct[label] = ktp_dct[spc0] + ktp_dct.pop(spc0) + else: + print('*Warning: hot spc {} not found in dct. dissociates completely?'.format(spc0)) label1 = list(label[1]) label1.remove(spc0) @@ -173,7 +176,7 @@ def bf_tp_df_full(ped_df, hotbf_df): return bf_tp_df -def bf_tp_df_todct(bf_tp_df, bf_threshold = 1e-2, savefile=False, rxn='', model=''): +def bf_tp_df_todct(bf_tp_df, bf_threshold = 1e-3, savefile=False, rxn='', model=''): """ Converts the dataframe of hot branching fractions to dictionary and excludes invalid BFs @@ -294,8 +297,9 @@ def bf_df_fromktpdct(ktp_dct, reac_str, temps, pressures): # warning for neg vals - absolute value if any(bf_series.values < 0): - print('Warning: found negative BFs at {:1.0f} K and {:1.1e} atm' - .format(temp, pressure)) + # remove warning- too verbose + # print('Warning: found negative BFs at {:1.0f} K and {:1.1e} atm' + # .format(temp, pressure)) bf_series = abs(bf_series) bf_tp_df.at[temp, pressure] = bf_series/np.sum(bf_series.values) #bf_tp_df[pressure][temp] = bf_series/np.sum(bf_series.values) From 7353c9772de3043e1e477d8b4a6e7364394a4e41 Mon Sep 17 00:00:00 2001 From: Luna Pratali Maffei Date: Thu, 14 Nov 2024 16:02:39 +0100 Subject: [PATCH 4/4] improved stupid/unstable ped and dos merging --- mechanalyzer/calculator/ene_partition.py | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/mechanalyzer/calculator/ene_partition.py b/mechanalyzer/calculator/ene_partition.py index 211bad1..af8f038 100644 --- a/mechanalyzer/calculator/ene_partition.py +++ b/mechanalyzer/calculator/ene_partition.py @@ -335,23 +335,21 @@ def init_dos(pressure, temp): # calculate rho_non1(ene1_vect) rho_non1 = [] - for idx in self.ene1_vect: - # first iter should be 0 + #for idx in self.ene1_vect: + ###old, relied on the fact that energy vector had spacing of 1, however not always accurate + for idx, _ in enumerate(self.ene1_vect[:-1]): # the sum of the energies in rhovib_prod2 and - # rho_trasl is always ene1 - idx_ene_int = np.arange(0, idx+1, dtype=int) + # rho_trasl is always ene1 + # ene1 = ene1_vect_w0[idx_ene_int] + ene1_vect_w0[idx_ene_int[::-1]] + # old # idx_ene_int = np.arange(0, idx+1, dtype=int) + idx_ene_int = np.arange(0, idx+2, dtype=int) idx_ene_minus_ene_int = idx_ene_int[::-1] rho_non1_integrand = ( rho_rovib_prod2[idx_ene_int] * rho_trasl[idx_ene_minus_ene_int] ) - try: - rho_non1.append(np.trapz(rho_non1_integrand, - x=self.ene1_vect[idx_ene_int])) - except IndexError: - continue - #probably just missed 1 index - #print(pressure, temp, self.ene1_vect, idx, idx_ene_int) + rho_non1.append(np.trapz(rho_non1_integrand, + x=self.ene1_vect[idx_ene_int])) rho_non1 = np.array(rho_non1)