diff --git a/thermo/bulk.py b/thermo/bulk.py
index 38b0417c..6bf598ef 100644
--- a/thermo/bulk.py
+++ b/thermo/bulk.py
@@ -372,8 +372,8 @@ class BulkSettings:
                 'mu_VL', 'mu_VL_power_exponent', 'k_LL', 'k_LL_power_exponent', 'k_VL', 'k_VL_power_exponent',
                 'sigma_LL', 'sigma_LL_power_exponent', 'T_liquid_volume_ref', 'T_normal', 'P_normal', 'T_standard',
                 'P_standard', 'T_gas_ref', 'P_gas_ref', 'speed_of_sound', 'kappa', 'isobaric_expansion',
-                'Joule_Thomson', 'VL_ID', 'VL_ID_settings', 'S_ID', 'S_ID_settings', 'solid_sort_method', 
-                'liquid_sort_method', 'liquid_sort_cmps', 'solid_sort_cmps', 'liquid_sort_cmps_neg', 'solid_sort_cmps_neg', 
+                'Joule_Thomson', 'VL_ID', 'VL_ID_settings', 'S_ID', 'S_ID_settings', 'solid_sort_method',
+                'liquid_sort_method', 'liquid_sort_cmps', 'solid_sort_cmps', 'liquid_sort_cmps_neg', 'solid_sort_cmps_neg',
                 'liquid_sort_prop', 'solid_sort_prop', 'phase_sort_higher_first', 'water_sort', 'equilibrium_perturbation')
 
     def as_json(self):
diff --git a/thermo/chemical_package.py b/thermo/chemical_package.py
index 06aacba8..3a1367d5 100644
--- a/thermo/chemical_package.py
+++ b/thermo/chemical_package.py
@@ -462,7 +462,7 @@ def compound_index(self, CAS=None, name=None, smiles=None, InChI=None,
 
     @property
     def CAS_to_index(self):
-        r'''Dictionary of CAS: index, used for efficiently 
+        r'''Dictionary of CAS: index, used for efficiently
         looking up which index a specified CAS number is in.
         This can save a lot of time, but does take a little more memory.
 
diff --git a/thermo/chemical_utils.py b/thermo/chemical_utils.py
index c72c2eaf..5231ccf3 100644
--- a/thermo/chemical_utils.py
+++ b/thermo/chemical_utils.py
@@ -165,7 +165,7 @@ def S0_basis_converter(c, S0_liq=None, S0_gas=None, T_ref=298.15):
     # 'S': {'Tm': 388.360, 'Tb': 882.117, 'Hfus': 1721, 'Hvap': 53326,  'solid_T_trans': [368.300], 'solid_H_trans': [401]}, # TODO
     # 'Be': {'Tm': 1560.000, 'Tb': 2741.437, 'Hfus': 7895, 'Hvap': 291572, 'solid_T_trans': [1527.000], 'solid_H_trans': [6849]}, # TODO
     # 'Fe': {'Tm': 3133.345, 'Hfus': 349585, 'solid_T_trans': [1184.000, 1665.000], 'solid_H_trans': [900, 837]}, # TODO
-    
+
     }
 
 shomate_gas_elements = ('H', 'O', 'N', 'F', 'P', 'Cl', 'Br', 'I', 'Mg', 'B', 'Pb', 'Li', 'Na', 'Al', 'K', 'V', 'Cr')
@@ -206,12 +206,12 @@ def element_HeatCapacitySolid_cache(CAS):
 def standard_state_ideal_gas_formation(c, T, Hf=None, Sf=None, T_ref=298.15):
     r'''This function calculates the standard state ideal-gas heat of formation
     of a compound at a specified from first principles. The entropy change and
-    Gibbs free energy change of formation are also returned. 
+    Gibbs free energy change of formation are also returned.
     This special condition is usually tabulated in thermodynamic tables, and
     this function is intended to be consistent with the standard conventions.
 
     The values returned depend on:
-    
+
     * The heat of formation as an ideal gas of the compound at 298.15 K (must be measured experimentally)
     * The ideal gas absolute entropy of the compound at 298.15 K (calculated from heat capacity and phase transitions of the compound and its various phases)
     * The chosen reference states of the elements (convention, but hardcoded)
@@ -260,7 +260,7 @@ def standard_state_ideal_gas_formation(c, T, Hf=None, Sf=None, T_ref=298.15):
     Hf_ref = Hf if Hf is not None else c.Hfgm
     Sf_ref = Sf if Sf is not None else c.Sfgm
     atoms = c.atoms
-    return _standard_state_ideal_gas_formation_direct(T=T, Hf_ref=Hf_ref, Sf_ref=Sf_ref, 
+    return _standard_state_ideal_gas_formation_direct(T=T, Hf_ref=Hf_ref, Sf_ref=Sf_ref,
             atoms=atoms, gas_Cp=c.HeatCapacityGas, T_ref=T_ref)
 
 _standard_formation_reaction_cache = {}
@@ -282,14 +282,14 @@ def _standard_state_ideal_gas_formation_direct(T, Hf_ref, Sf_ref, atoms, gas_Cp,
 
     dH_compound = gas_Cp.T_dependent_property_integral(T_ref, T)
     dS_compound = gas_Cp.T_dependent_property_integral_over_T(T_ref, T)
-    
+
     H_calc = reactant_coeff*Hf_ref + reactant_coeff*dH_compound
     S_calc = reactant_coeff*Sf_ref + reactant_coeff*dS_compound
     # if the compound is an element it will need special handling to go from solid liquid to gas if needed
 
     solid_ele = set(['C'])
     liquid_ele = set([''])
-    
+
     for coeff, ele_data in zip(elemental_counts, elemental_composition):
         ele = list(ele_data.keys())[0]
         element_obj = periodic_table[ele]
@@ -334,7 +334,7 @@ def _standard_state_ideal_gas_formation_direct(T, Hf_ref, Sf_ref, atoms, gas_Cp,
         elif ele == 'P':
             # White phosphorus is the basis here
             T_alpha_beta_P = 195.400
-            Htrans_alpha_beta_P = 521.0 # 525.5104 reported in 
+            Htrans_alpha_beta_P = 521.0 # 525.5104 reported in
             # The thermodynamic properties of elementary phosphorus The heat capacities of two crystalline modifications of red phosphorus, of α and β white phosphorus, and of black phosphorus from 15 to 300 K
             Tm_P = 317.300
             Hfus_P = 659
@@ -389,7 +389,7 @@ def _standard_state_ideal_gas_formation_direct(T, Hf_ref, Sf_ref, atoms, gas_Cp,
         H_calc -= coeff*dH_ele
         S_calc -= coeff*dS_ele
     G_calc = H_calc - T*S_calc
-    
+
     H_calc, S_calc, G_calc = H_calc/reactant_coeff, S_calc/reactant_coeff, G_calc/reactant_coeff
-    
+
     return H_calc, S_calc, G_calc
diff --git a/thermo/equilibrium.py b/thermo/equilibrium.py
index 7038ab52..9600c3eb 100644
--- a/thermo/equilibrium.py
+++ b/thermo/equilibrium.py
@@ -3395,7 +3395,7 @@ def get_bulk_property(self):
                'Gs_ideal_gas_standard_state', 'S_ideal_gas_standard_state', 'Ss_ideal_gas_standard_state',
 
                 'concentrations_mass_gas', 'concentrations_mass_gas_normal', 'concentrations_mass_gas_standard',
-                'concentrations_gas_standard', 'concentrations_gas_normal', 'concentrations_gas'              
+                'concentrations_gas_standard', 'concentrations_gas_normal', 'concentrations_gas'
  ]
 
 
diff --git a/thermo/fitting.py b/thermo/fitting.py
index 306c358a..315f2eae 100644
--- a/thermo/fitting.py
+++ b/thermo/fitting.py
@@ -68,17 +68,17 @@ def split_data(x, y, folds=5, seed=42):
     if pts != len(y):
         raise ValueError("Wrong size")
     z = list(range(pts))
-    
+
     Random(seed).shuffle(z)
     # go one by one and assign data to each group
-    
+
     fold_xs = [[] for _ in range(folds)]
     fold_ys = [[] for _ in range(folds)]
     for i in range(pts):
         l = i%folds
         fold_xs[l].append(x[i])
         fold_ys[l].append(y[i])
-    
+
     return fold_xs, fold_ys
 
 def assemble_fit_test_groups(x_groups, y_groups):
@@ -98,7 +98,7 @@ def assemble_fit_test_groups(x_groups, y_groups):
             if j != i:
                 x_train.extend(x_groups[j])
                 y_train.extend(y_groups[j])
-        
+
         train_x_groups.append(x_train)
         train_y_groups.append(y_train)
     return (train_x_groups, test_x_groups, train_y_groups, test_y_groups)
@@ -851,7 +851,7 @@ def jac_wrapped_for_leastsq(params):
         statistics['min_ratio'] = float(stats[2])
         statistics['max_ratio'] = float(stats[3])
         try:
-            pcov = float(pcov) 
+            pcov = float(pcov)
         except:
             # ndarray or None
             pass
diff --git a/thermo/flash/flash_base.py b/thermo/flash/flash_base.py
index 4ac7b86f..40c597cb 100644
--- a/thermo/flash/flash_base.py
+++ b/thermo/flash/flash_base.py
@@ -832,7 +832,7 @@ def _finish_initialization_base(self):
         self.scalar = scalar_statuses.pop()
 
         self.supports_lnphis_args = all(p.supports_lnphis_args for p in self.phases)
-        
+
         # Make the phases aware of the constants and properties
         constants = self.constants
         correlations = self.correlations
diff --git a/thermo/flash/flash_utils.py b/thermo/flash/flash_utils.py
index f75c875f..ab4ce7f2 100644
--- a/thermo/flash/flash_utils.py
+++ b/thermo/flash/flash_utils.py
@@ -70,7 +70,7 @@
     'generate_phase_boundaries_naive',
     'water_wet_bulb_temperature',
     'water_dew_point_from_humidity',
-    'solve_water_wet_bulb_temperature_nested', 
+    'solve_water_wet_bulb_temperature_nested',
     'solve_water_wet_bulb_temperature_direct',
 ]
 
@@ -5280,13 +5280,13 @@ def water_dew_point_from_humidity(flasher, T, P, humidity, zs_air, zs_added):
     except:
         return None
     N = flasher.N
-    ratio = sat.flash_convergence['mixing_factor'] 
+    ratio = sat.flash_convergence['mixing_factor']
     feed_ns = zs_air
     zs_for_factor = flash_mixing_remove_overlap(zs_air, zs_added)
     n_factor = sum([feed_ns[i] if zs_for_factor[i] > 0 else 0 for i in range(N)])
     ns_out = [n_factor*(zs_for_factor[i] + zs_added[i]*ratio*humidity) for i in range(N)]
     zs_my_flash = normalize(ns_out)
-    
+
     return water_wet_bulb_temperature(flasher, zs=zs_my_flash, T=T, P=P)
 
 
@@ -5413,10 +5413,10 @@ def solve_water_wet_bulb_temperature_direct(flasher, zs, T, P, T_wet_bulb):
     H_water_product = water_product.H()
     flashes2 = []
 
-    '''Known: T_wet_bulb. Product T. Product P.
+    """Known: T_wet_bulb. Product T. Product P.
 
     Unknown: x_w out, how many moles water to add to cause saturation of water, the temperature of the end flash
-    '''
+    """
     def wet_bulb_T_error_direct(guess):
         # TODO make a standalone function for this
         moles_water_1, moles_water_2 = guess
@@ -5439,7 +5439,7 @@ def wet_bulb_T_error_direct(guess):
         zs_out = normalize(ns_out_2)
 
         # As best as I can tell, the VF=1 flash is required otherwise the newton solver will
-        # go nuts with no jacobian. 
+        # go nuts with no jacobian.
         flash_out = flasher.flash(VF=1, P=P, zs=zs_out)
         flashes2[:] = [flash_out]
         energy_err = H_out_balanced - flash_out.H()
diff --git a/thermo/functional_groups.py b/thermo/functional_groups.py
index 87bf98b3..323ec997 100644
--- a/thermo/functional_groups.py
+++ b/thermo/functional_groups.py
@@ -2829,7 +2829,7 @@ def count_rings_attatched_to_rings(mol, allow_neighbors=True, atom_rings=None):
 benzene_smarts = 'c1ccccc1'
 
 def benene_rings(mol):
-    r'''Given a `rdkit.Chem.rdchem.Mol` object, returns the number of benzene rings 
+    r'''Given a `rdkit.Chem.rdchem.Mol` object, returns the number of benzene rings
     in the molecule.
 
     Parameters
@@ -2914,9 +2914,9 @@ def is_radionuclide(mol):
 
     References
     ----------
-    .. [1] Kondev, F. G., M. Wang, W. J. Huang, S. Naimi, and G. Audi. 
-       "The NUBASE2020 Evaluation of Nuclear Physics Properties *." 
-       Chinese Physics C 45, no. 3 (March 2021): 030001. 
+    .. [1] Kondev, F. G., M. Wang, W. J. Huang, S. Naimi, and G. Audi.
+       "The NUBASE2020 Evaluation of Nuclear Physics Properties *."
+       Chinese Physics C 45, no. 3 (March 2021): 030001.
        https://doi.org/10.1088/1674-1137/abddae.
     '''
     for atom in mol.GetAtoms():
diff --git a/thermo/heat_capacity.py b/thermo/heat_capacity.py
index 9fb01188..9ee461f2 100644
--- a/thermo/heat_capacity.py
+++ b/thermo/heat_capacity.py
@@ -1551,7 +1551,7 @@ class HeatCapacityLiquidMixture(MixtureProperty):
     ranked_methods = [LALIBERTE, LINEAR]
     pure_references = ('HeatCapacityLiquids',)
     pure_reference_types = (HeatCapacityLiquid,)
-    
+
     pure_constants = ('MWs', )
     custom_args = pure_constants
 
diff --git a/thermo/phases/gibbs_eos.py b/thermo/phases/gibbs_eos.py
index cc2be38d..918f8970 100644
--- a/thermo/phases/gibbs_eos.py
+++ b/thermo/phases/gibbs_eos.py
@@ -27,7 +27,7 @@
 
 class GibbsEOS(Phase):
     # http://www.iapws.org/relguide/Advise3.pdf is useful
-    def V(self):    
+    def V(self):
         return self.dG_dP()
         # return self._V
 
@@ -41,7 +41,7 @@ def dV_dP(self):
         return self.d2G_dP2()
 
     dV_dP_T = dV_dP
-    
+
     def dP_dV(self):
         return 1.0/self.d2G_dP2()
 
@@ -55,12 +55,12 @@ def dV_dT(self):
 
     def dT_dV(self):
         return 1.0/self.d2G_dTdP()
-    
+
     dT_dV_P = dT_dV
 
     def d2V_dT2(self):
         return self.d3G_dT2dP()
-    
+
     def d2V_dTdP(self):
         return self.d3G_dTdP2()
 
@@ -72,7 +72,7 @@ def dP_dT(self):
 
     def dT_dP(self):
         return 1.0/self.dP_dT()
-    
+
     dT_dP_V = dT_dP
 
 
@@ -81,10 +81,10 @@ def H(self):
 
     def S(self):
         return -self.dG_dT()
-    
+
     def Cp(self):
         return -self.T*self.d2G_dT2()
-    
+
     def U(self):
         return self.G() - self.T*self.dG_dT() - self.P*self.dG_dP()
 
@@ -109,6 +109,6 @@ def d2S_dP2(self):
 
     def d2S_dTdP(self):
         return -self.d3G_dT2dP()
-    
 
-    
+
+
diff --git a/thermo/phases/iapws_phase.py b/thermo/phases/iapws_phase.py
index 1379b600..c79c0de2 100644
--- a/thermo/phases/iapws_phase.py
+++ b/thermo/phases/iapws_phase.py
@@ -47,7 +47,7 @@ class IAPWS95(HelmholtzEOS):
     Hfs = [-241822.0]
     Sfs = [-44.5]
     Gfs = [-228554.325]
-    
+
     N = 1
     rhoc_inv = rho_to_Vm(rhoc_mass, _MW)
     rhoc = 1.0/rhoc_inv
@@ -756,17 +756,17 @@ def _set_core(self):
 
         tau2 = tau*tau
         v0 = pi - pi0
-        r2 = (v0*(v0*(2.34801409215913e-11 - 2.85651142904972e-11j) 
+        r2 = (v0*(v0*(2.34801409215913e-11 - 2.85651142904972e-11j)
                      + (-5.57107698030123e-5 + 4.64578634580806e-5j))
                      + (-72.597457432922 - 78.100842711287j))
 
-        r2p = (v0*(7.677551608692879e-14-9.340239477516712e-14j) 
+        r2p = (v0*(7.677551608692879e-14-9.340239477516712e-14j)
                + (-9.108171704568459e-08+7.595411065038183e-08j))
-        g0 = (v0*(v0*(v0*(-5.56464869058991e-22*v0 + 3.39746123271053e-15) 
+        g0 = (v0*(v0*(v0*(-5.56464869058991e-22*v0 + 3.39746123271053e-15)
                - 1.89369929326131e-8) + 0.655022213658955) - 632020.233335886)
         g0p = (v0*(v0*(-3.6390648292032365e-24*v0 + 1.666356094695489e-17)
                 - 6.192030151739651e-11) + 0.0010708979275295713)
-        g0pp = ((-1.7848556441943294e-26*v0 + 5.44866189611331e-20)*v0 
+        g0pp = ((-1.7848556441943294e-26*v0 + 5.44866189611331e-20)*v0
                 - 1.0123370045204503e-13)
 
 
@@ -774,17 +774,17 @@ def _set_core(self):
         log_t1_n_tau = logc(t1 - tau)
         log_t2_tau = logc(t2 + tau)
         log_t2_n_tau = logc(t2 - tau)
-        
+
         log_t1 = logc(t1)
         log_t2 = logc(t2)
-        
+
         t2_inv = 1.0/t2
         t1_inv = 1.0/t1
-        
+
         tau_t2inv = tau*t2_inv
         t2_log_t2 = t2*log_t2
         t2_log_t2_2 = t2_log_t2 + t2_log_t2
-        
+
         g_real_sum = (r1*((t1-tau)*log_t1_n_tau + (t1+tau)*log_t1_tau - 2.0*t1*log_t1 - tau2*t1_inv)
                    + (r2*((t2-tau)*log_t2_n_tau + (t2+tau)*log_t2_tau - t2_log_t2_2 - tau2*t2_inv))).real
         g = g0 - s0*Tt*tau + Tt*g_real_sum
@@ -798,7 +798,7 @@ def _set_core(self):
         g_TT = 0.0036608581051398447*g_TT_real_sum #1.0/Tt*g_TT_real_sum
 
         x0 = (t2-tau)*log_t2_n_tau + (t2+tau)*log_t2_tau - t2_log_t2_2 - tau2*t2_inv
-        
+
         g_P_real_sum = (r2p*(x0)).real
         g_P = g0p + Tt*g_P_real_sum
         g_PP_real_sum = (self.r2pp*(x0)).real
@@ -828,19 +828,19 @@ def G(self):
 
     def dG_dT(self):
         return self._dG_dT
-    
+
     def d2G_dT2(self):
         return self._d2G_dT2
 
     def dG_dP(self):
         return self._dG_dP
-    
+
     def d2G_dP2(self):
         return self._d2G_dP2
-    
+
     def d2G_dTdP(self):
         return self._d2G_dTdP
-        
+
     d2G_dPdT = d2G_dTdP
 
     def P_sub(self):
diff --git a/thermo/phases/phase.py b/thermo/phases/phase.py
index 240ea939..b03b748b 100644
--- a/thermo/phases/phase.py
+++ b/thermo/phases/phase.py
@@ -358,12 +358,12 @@ def is_same_model(self, other_phase, ignore_phase=False):
         if identical_model_attribute_ids:
             if ignore_phase:
                 return True
-            return self.__class__.__name__ == other_phase.__class__.__name__ 
+            return self.__class__.__name__ == other_phase.__class__.__name__
 
         # Using identities only we could not confirm if the phase was the same or not.
         # The values may still be the same if the identities are not.
         return self.model_hash(ignore_phase) == other_phase.model_hash(ignore_phase)
-    
+
     def model_hash(self, ignore_phase=False):
         r'''Method to compute a hash of a phase.
 
@@ -6363,7 +6363,7 @@ def _der(self, property=a, differentiate_by=b, at_constant=c):
              # These will probably need their doc fixed
              'd2G_dP2', 'd2G_dT2',
              'd2G_dPdT', 'd2G_dTdP',
-             
+
              ]
 derivatives_thermodynamic_mass = []
 
diff --git a/thermo/redlich_kister.py b/thermo/redlich_kister.py
index 3875134a..48e2788d 100644
--- a/thermo/redlich_kister.py
+++ b/thermo/redlich_kister.py
@@ -107,7 +107,7 @@ def redlich_kister_excess_inner(N, N_terms, a_tensor, xs):
 def redlich_kister_build_structure(N, shape, data, indexes):
     r'''Builds a redlich-kister compatible
     tensor (data structure) from pairs of indexes, and
-    the coefficients associated with those indexes. 
+    the coefficients associated with those indexes.
     This is especially important because of the asymmetry
     of the model.
 
@@ -136,7 +136,7 @@ def redlich_kister_build_structure(N, shape, data, indexes):
     and should be cached
     '''
     out = []
-    
+
     if len(shape) == 1:
         one_d = True
         two_d = False
@@ -145,12 +145,12 @@ def redlich_kister_build_structure(N, shape, data, indexes):
         two_d = True
     else:
         raise ValueError("Shape must be provided")
-        
+
     if len(indexes) != len(data):
         raise ValueError("Index and data length must be the same")
-        
+
     data_dict = {idx: d for idx, d in zip(indexes, data)}
-        
+
     for i in range(N):
         l = []
         for j in range(N):
@@ -159,7 +159,7 @@ def redlich_kister_build_structure(N, shape, data, indexes):
             else:
                 l.append([[0.0]*shape[1] for _ in range(shape[0])])
         out.append(l)
-        
+
     for i in range(N):
         for j in range(N):
             if (i, j) in data_dict:
@@ -248,7 +248,7 @@ def redlich_kister_excess_inner_binary(ais, xs):
     system. This calculation is optimized. This works
     with the same values of coefficients as `redlich_kister_excess_inner`
     but without the excess dimensionality of the input
-    data. 
+    data.
 
     Parameters
     ----------
@@ -288,7 +288,7 @@ def redlich_kister_excess_binary(coefficients, x0, T, N_T, N_terms):
     factor = 1.0
     for i in range(N_terms):
         offset = i*N_T
-        ai = coefficients[offset] 
+        ai = coefficients[offset]
         if N_T >= 2:
             ai +=  coefficients[offset+1]*Tinv
         if N_T >= 3:
diff --git a/thermo/stream.py b/thermo/stream.py
index 02c419ed..e91089fc 100644
--- a/thermo/stream.py
+++ b/thermo/stream.py
@@ -2000,7 +2000,7 @@ class EquilibriumStream(EquilibriumState):
     '''Creates an EquilibriumStream object, built off :obj:`EquilibriumState`
     to contain flow rate amounts, making mass and energy balances easier.
 
-    EquilibriumStreams can have their flow rate, state, and composition 
+    EquilibriumStreams can have their flow rate, state, and composition
     defined using any sufficient set of the following. Note that not all
     sets of specs will have a solution, or a unique solution, or an
     algorithm to solve the problem implemented.
@@ -2040,14 +2040,14 @@ class EquilibriumStream(EquilibriumState):
 
     Note that the liquid flow rates `Ql` and `Qls` will by default use the
     pure component liquid standard molar densities, but the temperature and
-    pressure used to calculate the liquid molar densities can be set with 
-    the two-tuple `Vf_TP`. 
+    pressure used to calculate the liquid molar densities can be set with
+    the two-tuple `Vf_TP`.
     See :obj:`EquilibriumState.V_liquids_ref` for details.
 
     Parameters
     ----------
-    flasher : One of :obj:`thermo.flash.FlashPureVLS`, :obj:`thermo.flash.FlashVL`, :obj:`thermo.flash.FlashVLN`, 
-        The configured flash object which can perform flashes for the 
+    flasher : One of :obj:`thermo.flash.FlashPureVLS`, :obj:`thermo.flash.FlashVL`, :obj:`thermo.flash.FlashVLN`,
+        The configured flash object which can perform flashes for the
         configured components, [-]
     zs : list, optional
         Mole fractions of all components [-]
@@ -2893,10 +2893,10 @@ def energy_balance(inlets, outlets, reactive=False, use_mass=False):
 
         return True
 
-    elif (in_unknown_count==1 and out_unknown_count == 1 and use_mass 
-        and isinstance(inlets[in_unknown_idx], StreamArgs) and isinstance(outlets[out_unknown_idx], StreamArgs) 
+    elif (in_unknown_count==1 and out_unknown_count == 1 and use_mass
+        and isinstance(inlets[in_unknown_idx], StreamArgs) and isinstance(outlets[out_unknown_idx], StreamArgs)
         and inlets[in_unknown_idx].state_specified and outlets[out_unknown_idx].state_specified):
-        '''
+        """
         from sympy import *
         m_in_known, m_in_unknown, m_out_known, m_out_unknown = symbols('m_in_known, m_in_unknown, m_out_known, m_out_unknown')
         e_in_known, e_out_known = symbols('e_in_known, e_out_known')
@@ -2907,7 +2907,7 @@ def energy_balance(inlets, outlets, reactive=False, use_mass=False):
 
         Eq0 = Eq(e_in_known+e_in_unknown, e_out_unknown+e_out_known)
         Eq1 = Eq(m_in_known+ m_in_unknown, m_out_known+m_out_unknown)
-        solve([Eq0, Eq1], [m_in_unknown, m_out_unknown])'''
+        solve([Eq0, Eq1], [m_in_unknown, m_out_unknown])"""
         unknown_in_state = inlets[in_unknown_idx].flash_state()
         unknown_out_state = outlets[out_unknown_idx].flash_state()
         H_mass_in_unknown = unknown_in_state.H_mass() if not reactive else unknown_in_state.H_reactive_mass()
@@ -2926,7 +2926,7 @@ def energy_balance(inlets, outlets, reactive=False, use_mass=False):
             if isinstance(f, StreamArgs) and f.state_specified:
                 unknown_inlet_idxs.append(i)
         if len(unknown_inlet_idxs) == 2:
-            '''
+            """
             from sympy import *
             m_in_known, m_in_unknown0, m_in_unknown1, m_out_known = symbols('m_in_known, m_in_unknown0, m_in_unknown1, m_out_known')
             e_in_known, e_known = symbols('e_in_known, e_known')
@@ -2938,7 +2938,7 @@ def energy_balance(inlets, outlets, reactive=False, use_mass=False):
             Eq0 = Eq(e_known, e_in_known + e_unkown1+e_unkown0 )
             Eq1 = Eq(m_in_known+ m_in_unknown0+ m_in_unknown1, m_out_known)
             solve([Eq0, Eq1], [m_in_unknown0, m_in_unknown1])
-            '''
+            """
             in_unknown_idx_0, in_unknown_idx_1 = unknown_inlet_idxs
             unknown_state_0 = inlets[in_unknown_idx_0].flash_state()
             unknown_state_1 = inlets[in_unknown_idx_1].flash_state()
@@ -2948,7 +2948,7 @@ def energy_balance(inlets, outlets, reactive=False, use_mass=False):
             m_in_known = sum(v.m for i, v in enumerate(inlets) if (i != in_unknown_idx_0 and i != in_unknown_idx_1 and not isinstance(v, EnergyStream)))
             m_out_known = sum(v.m for i, v in enumerate(outlets) if not isinstance(v, EnergyStream))
 
-            inlets[in_unknown_idx_0].m = (H_mass_in_unknown_1*m_in_known - H_mass_in_unknown_1*m_out_known - energy_in_known + outlet_energy)/(H_mass_in_unknown_0 - H_mass_in_unknown_1) 
+            inlets[in_unknown_idx_0].m = (H_mass_in_unknown_1*m_in_known - H_mass_in_unknown_1*m_out_known - energy_in_known + outlet_energy)/(H_mass_in_unknown_0 - H_mass_in_unknown_1)
             inlets[in_unknown_idx_1].m = (-H_mass_in_unknown_0*m_in_known + H_mass_in_unknown_0*m_out_known + energy_in_known - outlet_energy)/(H_mass_in_unknown_0 - H_mass_in_unknown_1)
             return True
     elif out_unknown_count == 2 and in_unknown_count == 0 and use_mass:
@@ -2967,7 +2967,7 @@ def energy_balance(inlets, outlets, reactive=False, use_mass=False):
             m_out_known = sum(v.m for i, v in enumerate(outlets) if (i != out_unknown_idx_0 and i != out_unknown_idx_1 and not isinstance(v, EnergyStream)))
             m_in_known = sum(v.m for i, v in enumerate(inlets) if not isinstance(v, EnergyStream))
 
-            outlets[out_unknown_idx_0].m = (H_mass_out_unknown_1*m_out_known - H_mass_out_unknown_1*m_in_known - energy_out_known + inlet_energy)/(H_mass_out_unknown_0 - H_mass_out_unknown_1) 
+            outlets[out_unknown_idx_0].m = (H_mass_out_unknown_1*m_out_known - H_mass_out_unknown_1*m_in_known - energy_out_known + inlet_energy)/(H_mass_out_unknown_0 - H_mass_out_unknown_1)
             outlets[out_unknown_idx_1].m = (-H_mass_out_unknown_0*m_out_known + H_mass_out_unknown_0*m_in_known + energy_out_known - inlet_energy)/(H_mass_out_unknown_0 - H_mass_out_unknown_1)
             return True
     return False
diff --git a/thermo/utils/mixture_property.py b/thermo/utils/mixture_property.py
index 88bafb33..a5084f01 100644
--- a/thermo/utils/mixture_property.py
+++ b/thermo/utils/mixture_property.py
@@ -935,7 +935,7 @@ def plot_isobaric_isothermal(self, T, P, methods=[], pts=50, only_valid=True,
             combined = list(product(iter_Ts, iter_Ps))
             iter_Ts = [T for T, _ in combined]
             iter_Ps = [P for _, P in combined]
-        
+
 
         xs = linspace(0, 1, pts)
         for T, P in zip(iter_Ts, iter_Ps):
@@ -1084,7 +1084,7 @@ def plot_binary(self, P=None, T=None, pts=30, Tmin=None, Tmax=None, Pmin=1E5,
         from matplotlib.ticker import FormatStrFormatter
         from numpy import ma
 
-        
+
 
         if T is None:
             # vary T
diff --git a/thermo/utils/t_dependent_property.py b/thermo/utils/t_dependent_property.py
index e2a891ba..8b7f4739 100644
--- a/thermo/utils/t_dependent_property.py
+++ b/thermo/utils/t_dependent_property.py
@@ -3328,7 +3328,7 @@ def add_piecewise_method(self, name, method_names, T_ranges):
         method_names : list[str]
             Each of the existing methods, [-]
         T_ranges : list[float]
-            A list of temperatures consisting of [Tmin, T_transition1, ..., Tmax]; 
+            A list of temperatures consisting of [Tmin, T_transition1, ..., Tmax];
             this is size len(method_names) + 1, [K]
 
         Notes
@@ -4737,7 +4737,7 @@ def __init__(self, extrapolation, **kwargs):
             piecewise_methods = kwargs['piecewise_methods']
             for name, data in piecewise_methods.items():
                 self.add_piecewise_method(name=name, method_names=data['methods'], T_ranges=data['T_ranges'])
-        
+
         try:
             method = kwargs['method']
         except:
diff --git a/thermo/utils/tp_dependent_property.py b/thermo/utils/tp_dependent_property.py
index 34db912d..b32641ae 100644
--- a/thermo/utils/tp_dependent_property.py
+++ b/thermo/utils/tp_dependent_property.py
@@ -107,7 +107,7 @@ def __init__(self, extrapolation, **kwargs):
         super().__init__(extrapolation, **kwargs)
 
         self.P_limits = {}
-        '''Pressure limits on a per-component basis. Not currently used.'''
+        """Pressure limits on a per-component basis. Not currently used."""
 
         if kwargs:
             P_correlation_keys_to_parameters = self.P_correlation_keys_to_parameters
diff --git a/thermo/vapor_pressure.py b/thermo/vapor_pressure.py
index 46afb0b0..3faae807 100644
--- a/thermo/vapor_pressure.py
+++ b/thermo/vapor_pressure.py
@@ -301,11 +301,11 @@ class VaporPressure(TDependentProperty):
        Equations for the Metallic Elements: 298-2500K." Canadian Metallurgical
        Quarterly 23, no. 3 (July 1, 1984): 309-13.
        https://doi.org/10.1179/cmq.1984.23.3.309.
-    .. [8] Hall, K. R. Vapor Pressure and Antoine Constants for Hydrocarbons, 
+    .. [8] Hall, K. R. Vapor Pressure and Antoine Constants for Hydrocarbons,
        and S, Se, Te, and Halogen Containing Organic Compounds. Springer, 1999.
     .. [9] Dykyj, J., and K. R. Hall. "Vapor Pressure and Antoine Constants for
        Oxygen Containing Organic Compounds". 2000.
-    .. [10] Hall, K. R. Vapor Pressure and Antoine Constants for Nitrogen 
+    .. [10] Hall, K. R. Vapor Pressure and Antoine Constants for Nitrogen
        Containing Organic Compounds. Springer, 2001.
     '''
 
@@ -718,11 +718,11 @@ class SublimationPressure(TDependentProperty):
        Equations for the Metallic Elements: 298-2500K." Canadian Metallurgical
        Quarterly 23, no. 3 (July 1, 1984): 309-13.
        https://doi.org/10.1179/cmq.1984.23.3.309.
-    .. [3] Hall, K. R. Vapor Pressure and Antoine Constants for Hydrocarbons, 
+    .. [3] Hall, K. R. Vapor Pressure and Antoine Constants for Hydrocarbons,
        and S, Se, Te, and Halogen Containing Organic Compounds. Springer, 1999.
     .. [4] Dykyj, J., and K. R. Hall. "Vapor Pressure and Antoine Constants for
        Oxygen Containing Organic Compounds". 2000.
-    .. [5] Hall, K. R. Vapor Pressure and Antoine Constants for Nitrogen 
+    .. [5] Hall, K. R. Vapor Pressure and Antoine Constants for Nitrogen
        Containing Organic Compounds. Springer, 2001.
     '''