diff --git a/arkane/encorr/bac.py b/arkane/encorr/bac.py
index c511a9f47b4..c2e2ccbda07 100644
--- a/arkane/encorr/bac.py
+++ b/arkane/encorr/bac.py
@@ -1092,10 +1092,13 @@ def get_confidence_intervals(x: np.ndarray,
         weights = np.eye(n)
 
     e = y - ypred  # Residuals
-    sigma2 = e.T @ weights @ e / (n - p)  # MSE
+    dof = n - p
+    if dof <= 0:
+        return np.full(p, np.nan), np.full((p, p), np.nan)
+    sigma2 = e.T @ weights @ e / dof  # MSE
     cov = sigma2 * np.linalg.inv(x.T @ weights @ x)  # covariance matrix
     se = np.sqrt(np.diag(cov))  # standard error
-    tdist = distributions.t.ppf(1 - alpha / 2, n - p)  # student-t
+    tdist = distributions.t.ppf(1 - alpha / 2, dof)  # student-t
     ci = tdist * se  # confidence interval half-width
     return ci, cov
 
@@ -1103,6 +1106,7 @@ def get_confidence_intervals(x: np.ndarray,
 def _covariance_to_correlation(cov: np.ndarray) -> np.ndarray:
     """Convert (unscaled) covariance matrix to correlation matrix"""
     v = np.sqrt(np.diag(cov))
-    corr = cov / np.outer(v, v)
+    with np.errstate(divide='ignore', invalid='ignore'):
+        corr = cov / np.outer(v, v)
     corr[cov == 0] = 0
     return corr
diff --git a/arkane/output.py b/arkane/output.py
index 7657834c004..1fc8bd1f690 100644
--- a/arkane/output.py
+++ b/arkane/output.py
@@ -56,6 +56,20 @@ def __init__(self, level=0, indent=4):
         self.level = level
         self.indent = indent
 
+    @staticmethod
+    def _is_string(node):
+        return isinstance(node, ast.Constant) and isinstance(node.value, str)
+
+    @staticmethod
+    def _is_number(node):
+        return isinstance(node, ast.Constant) and isinstance(node.value, (int, float, complex)) and not isinstance(node.value, bool)
+
+    @classmethod
+    def _is_simple_constant(cls, node):
+        return cls._is_string(node) or cls._is_number(node) or (
+            isinstance(node, ast.Constant) and isinstance(node.value, bool)
+        )
+
     def visit_Call(self, node):
         """
         Return a pretty representation of the class or function call represented by `node`.
@@ -84,7 +98,7 @@ def visit_List(self, node):
         """
         Return a pretty representation of the list represented by `node`.
         """
-        if any([not isinstance(e, (ast.Str, ast.Num, ast.UnaryOp)) for e in node.elts]):
+        if any([not (self._is_simple_constant(e) or isinstance(e, ast.UnaryOp)) for e in node.elts]):
             # Split elements onto multiple lines
             result = '[\n'
             self.level += 1
@@ -105,11 +119,11 @@ def visit_Tuple(self, node):
         """
         # If the tuple represents a quantity, keep it on one line
         is_quantity = True
-        if len(node.elts) == 0 or not isinstance(node.elts[0], (ast.Num, ast.List)) or (
+        if len(node.elts) == 0 or not (self._is_number(node.elts[0]) or isinstance(node.elts[0], ast.List)) or (
                 isinstance(node.elts[0], ast.List) and
-                any([not isinstance(e, (ast.Num, ast.UnaryOp)) for e in node.elts[0].elts])):
+                any([not (self._is_number(e) or isinstance(e, ast.UnaryOp)) for e in node.elts[0].elts])):
             is_quantity = False
-        elif len(node.elts) < 2 or not isinstance(node.elts[1], ast.Str):
+        elif len(node.elts) < 2 or not self._is_string(node.elts[1]):
             is_quantity = False
 
         if not is_quantity:
@@ -131,8 +145,8 @@ def visit_Dict(self, node):
         """
         Return a pretty representation of the dict represented by `node`.
         """
-        if (any([not isinstance(e, (ast.Str, ast.Num)) for e in node.keys])
-                or any([not isinstance(e, (ast.Str, ast.Num)) for e in node.values])):
+        if (any([not self._is_simple_constant(e) for e in node.keys])
+                or any([not self._is_simple_constant(e) for e in node.values])):
             # Split elements onto multiple lines
             result = '{\n'
             self.level += 1
@@ -149,20 +163,18 @@ def visit_Dict(self, node):
             self.string = result
             return result
 
-    def visit_Str(self, node):
-        """
-        Return a pretty representation of the string represented by `node`.
-        """
-        result = repr(node.s)
-        self.string = result
-        return result
-
-    def visit_Num(self, node):
+    def visit_Constant(self, node):
         """
-        Return a pretty representation of the number represented by `node`.
+        Return a pretty representation of the constant represented by `node`.
         """
-        result = '{0:g}'.format(node.n)
-        # result = repr(node.n)
+        if isinstance(node.value, str):
+            result = repr(node.value)
+        elif isinstance(node.value, bool):
+            result = repr(node.value)
+        elif isinstance(node.value, (int, float, complex)) and not isinstance(node.value, bool):
+            result = '{0:g}'.format(node.value)
+        else:
+            return None
         self.string = result
         return result
 
diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index 0b43a7134cd..685d92cd778 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -43,6 +43,7 @@
 from copy import deepcopy
 
 import numpy as np
+from scipy.optimize import OptimizeWarning
 from sklearn.model_selection import KFold
 
 from rmgpy import settings
@@ -4654,13 +4655,20 @@ def _make_rule(rr):
             kin.uncertainty = RateUncertainty(mu=0.0, var=(np.log(fmax) / 2.0) ** 2, N=1, Tref=Tref, data_mean=data_mean, correlation=label)
         else:
             if isinstance(rs[0].kinetics, Arrhenius):
-                dlnks = np.array([
-                    np.log(
-                        arr().fit_to_reactions(np.delete(rs, i), recipe=recipe)
-                .to_arrhenius(rxn.get_enthalpy_of_reaction(Tref))
-                .get_rate_coefficient(T=Tref) / rxn.get_rate_coefficient(T=Tref)
-            ) for i, rxn in enumerate(rs)
-        ])  # 1) fit to set of reactions without the current reaction (k)  2) compute log(kfit/kactual) at Tref
+                with warnings.catch_warnings():
+                    # Jackknife refits use sparse samples and only consume fitted rates, not covariance estimates.
+                    warnings.filterwarnings(
+                        "ignore",
+                        message="Covariance of the parameters could not be estimated",
+                        category=OptimizeWarning,
+                    )
+                    dlnks = np.array([
+                        np.log(
+                            arr().fit_to_reactions(np.delete(rs, i), recipe=recipe)
+                            .to_arrhenius(rxn.get_enthalpy_of_reaction(Tref))
+                            .get_rate_coefficient(T=Tref) / rxn.get_rate_coefficient(T=Tref)
+                        ) for i, rxn in enumerate(rs)
+                    ])  # 1) fit to set of reactions without the current reaction (k)  2) compute log(kfit/kactual) at Tref
             else:
                 dlnks = np.array([
                     np.log(
diff --git a/rmgpy/data/solvation.py b/rmgpy/data/solvation.py
index bfa1fcb9661..585cbd2e284 100644
--- a/rmgpy/data/solvation.py
+++ b/rmgpy/data/solvation.py
@@ -2381,7 +2381,7 @@ def get_Kfactor_parameters(self, delG298, delH298, delS298, solvent_name, T_tran
 
         # Generate Amatrix and bvector for Ax = b
         Amatrix = np.zeros((4, 4))
-        bvec = np.zeros((4, 1))
+        bvec = np.zeros(4)
         # 1. Tr*ln(K-factor) value at T = 298 K
         rho_g_298 = get_gas_saturation_density(solvent_name, 298)
         rho_l_298 = get_liquid_saturation_density(solvent_name, 298)
diff --git a/rmgpy/rmg/pdep.py b/rmgpy/rmg/pdep.py
index 237cdada2a4..f6daa5292c1 100644
--- a/rmgpy/rmg/pdep.py
+++ b/rmgpy/rmg/pdep.py
@@ -432,9 +432,14 @@ def solve_ss_network(self, T, P):
             ind = isomer_spcs.index(self.source[0])
             b[ind] = -1.0  # flux at source
         else:
-            b = -b / b.sum()  # 1.0 flux from source
+            total_source_flux = b.sum()
+            if total_source_flux == 0:
+                return None
+            b = -b / total_source_flux  # 1.0 flux from source
 
         if len(b) == 1:
+            if A[0, 0] == 0:
+                return None
             return np.array([b[0] / A[0, 0]])
 
         con = np.linalg.cond(A)
diff --git a/rmgpy/statmech/conformer.pyx b/rmgpy/statmech/conformer.pyx
index a391484fc35..e3e911248db 100644
--- a/rmgpy/statmech/conformer.pyx
+++ b/rmgpy/statmech/conformer.pyx
@@ -541,6 +541,15 @@ cpdef double phi(double beta, int k, double E, logQ) except -10000000:
     T = 1.0 / (constants.R * beta)
     return logQ(T) - k * log(beta) + beta * E
 
+
+def _phi_for_optimizer(beta, int k, double E, logQ):
+    """
+    Evaluate :func:`phi` from SciPy optimizers, which pass the scalar search
+    coordinate as a one-element array.
+    """
+    return phi(float(np.asarray(beta).item()), k, E, logQ)
+
+
 @cython.boundscheck(False)
 @cython.wraparound(False)
 def get_density_of_states_forst(np.ndarray[np.float64_t, ndim=1] e_list, logQ, int order=1):
@@ -557,6 +566,7 @@ def get_density_of_states_forst(np.ndarray[np.float64_t, ndim=1] e_list, logQ, i
 
     cdef np.ndarray[np.float64_t, ndim=1] dens_states, sum_states
     cdef double x, dx, v, E, dE
+    cdef object opt_result
     cdef int i, k
 
     if order != 1 and order != 2:
@@ -580,10 +590,10 @@ def get_density_of_states_forst(np.ndarray[np.float64_t, ndim=1] e_list, logQ, i
 
         # Find minimum of phi  func x0  arg     xtol  ftol maxi  maxf fullout  disp retall  callback
         try:
-            x = scipy.optimize.fmin(phi, x, (k, E, logQ), 1e-8, 1e-8, 100, 1000, False, False, False, None)
+            opt_result = scipy.optimize.fmin(_phi_for_optimizer, np.array([x]), (k, E, logQ), 1e-8, 1e-8, 100, 1000, False, False, False, None)
         except ValueError:
             break
-        x = float(x)
+        x = float(opt_result.item())
         dx = 1e-2 * x
 
         # Evaluate derivatives needed for steepest descents approximation numerically
diff --git a/test/arkane/arkaneOutputTest.py b/test/arkane/arkaneOutputTest.py
index adde2c6873a..fa1347765d8 100644
--- a/test/arkane/arkaneOutputTest.py
+++ b/test/arkane/arkaneOutputTest.py
@@ -146,8 +146,8 @@ def test_prettify(self):
                 ],
                 'kJ/mol',
             ),
-            quantum = None,
-            semiclassical = None,
+            quantum = True,
+            semiclassical = False,
         ),
     ],
     spin_multiplicity = 2,
diff --git a/test/arkane/ess/arkaneGaussianTest.py b/test/arkane/ess/arkaneGaussianTest.py
index 8109b708a96..6c2a0bdc632 100644
--- a/test/arkane/ess/arkaneGaussianTest.py
+++ b/test/arkane/ess/arkaneGaussianTest.py
@@ -93,10 +93,10 @@ def test_load_ethylene_from_gaussian_log_cbsqb3(self):
         trans = [mode for mode in conformer.modes if isinstance(mode, IdealGasTranslation)][0]
         rot = [mode for mode in conformer.modes if isinstance(mode, NonlinearRotor)][0]
         vib = [mode for mode in conformer.modes if isinstance(mode, HarmonicOscillator)][0]
-        t_list = np.array([298.15], float)
-        assert abs(trans.get_partition_function(t_list) - 5.83338e6) < 1e1
-        assert abs(rot.get_partition_function(t_list) - 2.59622e3) < 1e-2
-        assert abs(vib.get_partition_function(t_list) - 1.0481e0) < 1e-4
+        t = 298.15
+        assert abs(trans.get_partition_function(t) - 5.83338e6) < 1e1
+        assert abs(rot.get_partition_function(t) - 2.59622e3) < 1e-2
+        assert abs(vib.get_partition_function(t) - 1.0481e0) < 1e-4
 
         assert round(abs(e0 / constants.Na / constants.E_h - -78.467452), 4) == 0
         assert conformer.spin_multiplicity == 1
@@ -152,10 +152,10 @@ def test_load_oxygen_from_gaussian_log(self):
         trans = [mode for mode in conformer.modes if isinstance(mode, IdealGasTranslation)][0]
         rot = [mode for mode in conformer.modes if isinstance(mode, LinearRotor)][0]
         vib = [mode for mode in conformer.modes if isinstance(mode, HarmonicOscillator)][0]
-        t_list = np.array([298.15], float)
-        assert abs(trans.get_partition_function(t_list) - 7.11169e6) < 1e1
-        assert abs(rot.get_partition_function(t_list) - 7.13316e1) < 1e-4
-        assert abs(vib.get_partition_function(t_list) - 1.00037e0) < 1e-4
+        t = 298.15
+        assert abs(trans.get_partition_function(t) - 7.11169e6) < 1e1
+        assert abs(rot.get_partition_function(t) - 7.13316e1) < 1e-4
+        assert abs(vib.get_partition_function(t) - 1.00037e0) < 1e-4
 
         assert round(abs(e0 / constants.Na / constants.E_h - -150.3784877), 4) == 0
         assert conformer.spin_multiplicity == 3
@@ -180,11 +180,10 @@ def test_load_ethylene_from_gaussian_log_g3(self):
         trans = [mode for mode in conformer.modes if isinstance(mode, IdealGasTranslation)][0]
         rot = [mode for mode in conformer.modes if isinstance(mode, NonlinearRotor)][0]
         vib = [mode for mode in conformer.modes if isinstance(mode, HarmonicOscillator)][0]
-        t_list = np.array([298.15], float)
-
-        assert abs(trans.get_partition_function(t_list) - 5.83338e6) < 1e1
-        assert abs(rot.get_partition_function(t_list) - 2.53410e3) < 1e-2
-        assert abs(vib.get_partition_function(t_list) - 1.0304e0) < 1e-4
+        t = 298.15
+        assert abs(trans.get_partition_function(t) - 5.83338e6) < 1e1
+        assert abs(rot.get_partition_function(t) - 2.53410e3) < 1e-2
+        assert abs(vib.get_partition_function(t) - 1.0304e0) < 1e-4
 
         assert round(abs(e0 / constants.Na / constants.E_h - -78.562189), 4) == 0
         assert conformer.spin_multiplicity == 1
diff --git a/test/arkane/ess/molproTest.py b/test/arkane/ess/molproTest.py
index 36a121f6238..afeb864fce2 100644
--- a/test/arkane/ess/molproTest.py
+++ b/test/arkane/ess/molproTest.py
@@ -34,8 +34,6 @@
 import os
 
 
-import numpy as np
-
 import rmgpy.constants as constants
 from rmgpy.statmech import (
     IdealGasTranslation,
@@ -129,11 +127,10 @@ def test_load_hosi_from_molpro_log(self):
         trans = [mode for mode in conformer.modes if isinstance(mode, IdealGasTranslation)][0]
         rot = [mode for mode in conformer.modes if isinstance(mode, NonlinearRotor)][0]
         vib = [mode for mode in conformer.modes if isinstance(mode, HarmonicOscillator)][0]
-        t_list = np.array([298.15], float)
-
-        assert abs(trans.get_partition_function(t_list) - 9.175364e7) < 1e1
-        assert abs(rot.get_partition_function(t_list) - 1.00005557e5) < 1e-2
-        assert abs(vib.get_partition_function(t_list) - 1.9734989e0) < 1e-4
+        t = 298.15
+        assert abs(trans.get_partition_function(t) - 9.175364e7) < 1e1
+        assert abs(rot.get_partition_function(t) - 1.00005557e5) < 1e-2
+        assert abs(vib.get_partition_function(t) - 1.9734989e0) < 1e-4
 
         assert round(abs(e0 / constants.Na / constants.E_h - -768.275662), 4) == 0
         assert conformer.spin_multiplicity == 1