Merge pull request #317 from lpratalimaffei/dev

Dev

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)

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)
 

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

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]