diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 7f1c8d127e1..8b060ca8147 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -174,6 +174,44 @@ jobs:
       - name: Test Cantera
         run: scons test
 
+  run-examples:
+    name: Run the Python examples using bash
+    runs-on: ubuntu-18.04
+    steps:
+      - uses: actions/checkout@v2
+        name: Checkout the repository
+        with:
+          submodules: recursive
+      - name: Setup Python
+        uses: actions/setup-python@v1
+        with:
+          python-version: '3.8'
+          architecture: x64
+      - name: Install Apt dependencies
+        run: sudo apt-get install libboost-dev gfortran graphviz liblapack-dev libblas-dev
+      - name: Upgrade pip
+        run: python3 -m pip install -U pip setuptools wheel
+      - name: Install Python dependencies
+        run: python3 -m pip install ruamel.yaml scons numpy cython h5py pandas matplotlib scipy
+      - name: Build Cantera
+        run: python3 `which scons` build -j2
+      - name: Run the examples
+        # See https://unix.stackexchange.com/a/392973 for an explanation of the -exec part
+        run: |
+          find interfaces/cython/cantera/examples -type f -iname "*.py" \
+          -exec sh -c 'for n; do echo "$n" | tee -a results.txt && python3 "$n" >> results.txt || exit 1; done' sh {} +
+        env:
+          PYTHONPATH: build/python
+          PYTHONWARNINGS: error
+          MPLBACKEND: Agg
+      - name: Save the results file for inspection
+        uses: actions/upload-artifact@v2
+        with:
+          path: results.txt
+        # Run this step if the job was successful or failed, but not if it was cancelled
+        # Using always() would run this step if the job was cancelled as well.
+        if: failure() || success()
+
   multiple-sundials:
     name: Sundials ${{ matrix.sundials-ver }}
     runs-on: ubuntu-latest
diff --git a/data/sofc.yaml b/data/sofc.yaml
index cbe0b95e7c9..2421ca57172 100644
--- a/data/sofc.yaml
+++ b/data/sofc.yaml
@@ -5,7 +5,7 @@ description: |-
   species thermochemistry ARE NOT REAL VALUES - they are chosen only
   for the purposes of this example.
 
-  Defines bulk (i.e., 3D) phases - a gas, a metal, and an oxide.
+  Defines bulk (that is, 3D) phases - a gas, a metal, and an oxide.
   The gas contains only the minimum number of species needed to model
   operation on hydrogen. The species definitions are imported from
   gri30.yaml. The initial composition is set to hydrogen + 5% water, but
@@ -26,14 +26,14 @@ description: |-
   charge balances on the metal, then we would require a more complex
   model. Note that there is no work function for this metal.
 
-  Note: the "const_cp" species thermo model is used throughout this
+  Note: the 'const_cp' species thermo model is used throughout this
   file (with the exception of the gaseous species, which use NASA
-  polynomials imported from gri30.cti). The const_cp model assumes a
+  polynomials imported from gri30.yaml). The const_cp model assumes a
   constant specific heat, which by default is zero. Parameters that
-  can be specified are cp0, t0, h0, and s0. If omitted, t0 = 300 K, h0
-  = 0, and s0 = 0. The thermo properties are computed as follows: h =
-  h0 + cp0*(t - t0), s = s0 + cp0*ln(t/t0). For work at a single
-  temperature, it is sufficient to specify only h0.
+  can be specified are cp0, t0, h0, and s0. If omitted, t0 = 300 K,
+  h0 = 0, and s0 = 0. The thermo properties are computed as follows:
+  h = h0 + cp0*(t - t0), s = s0 + cp0*ln(t/t0).
+  For work at a single temperature, it is sufficient to specify only h0.
 
   The 'oxide_bulk' phase is a very simple model for the bulk phase. We only
   consider the oxygen sublattice. The only species we define are a
@@ -41,7 +41,7 @@ description: |-
   required input, but is not used here, so may be set arbitrarily.
 
   The vacancy will be modeled as truly vacant - it contains no atoms,
-  has no charge, and has zero enthalpy and entropy.  This is different
+  has no charge, and has zero enthalpy and entropy. This is different
   from the usual convention in which the vacancy properties are are
   expressed relative to the perfect crystal lattice. For example, in
   the usual convention, an oxygen vacancy has charge +2. But the
@@ -65,7 +65,7 @@ description: |-
   'mol' and 'cm' were specified in the units directive below as the
   units for quantity and length, respectively.
   2. The 'reactions' field specifies that all reaction entries in this file
-  that have are in the 'metal_surface-reactions' field are reactions belonging
+  that are in the 'metal_surface-reactions' field are reactions belonging
   to this surface mechanism.
 
   On the oxide surface, we consider four species:
@@ -74,11 +74,11 @@ description: |-
       3. OH'(ox)  - a surface hydroxyl with charge -1
       4. H2O(ox)  - physisorbed neutral water
 
-  The 'tpb' phase is th etriple phase boundary between the metal, oxide, and gas. A
-  single species is specified, but it is not used, since all reactions
-  only involve species on either side of the tpb. Note that the site
-  density is in mol/cm. But since no reactions involve TPB species,
-  this parameter is unused.
+  The triple phase boundary (TPB) between the metal, oxide, and gas is
+  specified in the 'tpb' phase. A single species is specified, but it
+  is not used, since all reactions only involve species on either side
+  of the TPB. Note that the site density is in mol/cm. But since no
+  reactions involve TPB species, this parameter is unused.
 
 generator: cti2yaml
 cantera-version: 2.5.0a3
diff --git a/include/cantera/kinetics/RateCoeffMgr.h b/include/cantera/kinetics/RateCoeffMgr.h
index 8449e1b9601..4f89f464a74 100644
--- a/include/cantera/kinetics/RateCoeffMgr.h
+++ b/include/cantera/kinetics/RateCoeffMgr.h
@@ -43,10 +43,10 @@ class Rate1
 
     /**
      * Update the concentration-dependent parts of the rate coefficient, if any.
-     * Used by class SurfaceArrhenius to compute coverage-dependent *
+     * Used by class SurfaceArrhenius to compute coverage-dependent
      * modifications to the Arrhenius parameters. The array c should contain
      * whatever data the particular rate coefficient class needs to update its
-     * rates.  Note that this method does not return anything. To get the
+     * rates. Note that this method does not return anything. To get the
      * updated rates, method update must be called after the call to update_C.
      */
     void update_C(const doublereal* c) {
diff --git a/include/cantera/kinetics/RxnRates.h b/include/cantera/kinetics/RxnRates.h
index 8ade9c2fb37..a4ada5fbd11 100644
--- a/include/cantera/kinetics/RxnRates.h
+++ b/include/cantera/kinetics/RxnRates.h
@@ -162,7 +162,7 @@ class SurfaceArrhenius
     }
 
     /**
-     * Update the value the rate constant.
+     * Update the value of the rate constant.
      *
      * This function returns the actual value of the rate constant. It can be
      * safely called for negative values of the pre-exponential factor.
diff --git a/include/cantera/kinetics/StoichManager.h b/include/cantera/kinetics/StoichManager.h
index 9077ed11edd..66c2a138bae 100644
--- a/include/cantera/kinetics/StoichManager.h
+++ b/include/cantera/kinetics/StoichManager.h
@@ -40,7 +40,7 @@ namespace Cantera
  * \f]
  *
  * where \f$ \nu^{(p)_{k,i}} \f$ is the product-side stoichiometric
- * coefficient of species \a k in reaction \a i. This could be done be
+ * coefficient of species \a k in reaction \a i. This could be done by
  * straightforward matrix multiplication, but would be inefficient, since most
  * of the matrix elements of \f$ \nu^{(p)}_{k,i} \f$ are zero. We could do
  * better by using sparse-matrix algorithms to compute this product.
@@ -61,12 +61,12 @@ namespace Cantera
  * computing quantities that can be written as matrix multiplies.
  *
  * They are designed to explicitly unroll loops over species or reactions for
- * Operations on reactions that require knowing the reaction stoichiometry.
+ * operations on reactions that require knowing the reaction stoichiometry.
  *
  * This module consists of class StoichManager, and classes C1, C2, and C3.
  * Classes C1, C2, and C3 handle operations involving one, two, or three
  * species, respectively, in a reaction. Instances are instantiated with a
- * reaction number, and n species numbers (n = 1 for C1, etc.).  All three
+ * reaction number, and n species numbers (n = 1 for C1, etc.). All three
  * classes have the same interface.
  *
  * These classes are designed for use by StoichManager, and the operations
diff --git a/include/cantera/oneD/Boundary1D.h b/include/cantera/oneD/Boundary1D.h
index 82dc4e68062..f1a7807b92a 100644
--- a/include/cantera/oneD/Boundary1D.h
+++ b/include/cantera/oneD/Boundary1D.h
@@ -299,6 +299,10 @@ class ReactingSurf1D : public Boundary1D
         m_enabled = docov;
     }
 
+    bool coverageEnabled() {
+        return m_enabled;
+    }
+
     virtual std::string componentName(size_t n) const;
 
     virtual void init();
diff --git a/include/cantera/thermo/EdgePhase.h b/include/cantera/thermo/EdgePhase.h
index 22d982a48f0..2825fde0a58 100644
--- a/include/cantera/thermo/EdgePhase.h
+++ b/include/cantera/thermo/EdgePhase.h
@@ -22,8 +22,8 @@ namespace Cantera
  * thermodynamic object.
  *
  * All of the equations and formulations carry through from SurfPhase to this
- * EdgePhase object. It should be noted however, that dimensional object with
- * length dimensions, have their dimensions reduced by one.
+ * EdgePhase object. It should be noted that dimensional quantities with
+ * dimensions including a length have their dimensions reduced by one.
  *
  * @ingroup thermoprops
  */
diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd
index 9a8dcf5eeb8..26ff56dcd3c 100644
--- a/interfaces/cython/cantera/_cantera.pxd
+++ b/interfaces/cython/cantera/_cantera.pxd
@@ -757,6 +757,7 @@ cdef extern from "cantera/oneD/Boundary1D.h":
         CxxRreactingSurf1D()
         void setKineticsMgr(CxxInterfaceKinetics*) except +translate_exception
         void enableCoverageEquations(cbool) except +translate_exception
+        cbool coverageEnabled()
 
 
 cdef extern from "cantera/oneD/StFlow.h":
diff --git a/interfaces/cython/cantera/examples/kinetics/extract_submechanism.py b/interfaces/cython/cantera/examples/kinetics/extract_submechanism.py
index aa52aefea0d..19a7c9b0d81 100644
--- a/interfaces/cython/cantera/examples/kinetics/extract_submechanism.py
+++ b/interfaces/cython/cantera/examples/kinetics/extract_submechanism.py
@@ -14,7 +14,8 @@
 import cantera as ct
 import matplotlib.pyplot as plt
 
-all_species = ct.Species.listFromFile('gri30.yaml')
+input_file = 'gri30.yaml'
+all_species = ct.Species.listFromFile(input_file)
 species = []
 
 # Filter species
@@ -36,7 +37,8 @@
 print('Species: {0}'.format(', '.join(S.name for S in species)))
 
 # Filter reactions, keeping only those that only involve the selected species
-all_reactions = ct.Reaction.listFromFile('gri30.yaml')
+ref_phase = ct.Solution(thermo='ideal-gas', kinetics='gas', species=all_species)
+all_reactions = ct.Reaction.listFromFile(input_file, ref_phase)
 reactions = []
 
 print('\nReactions:')
@@ -51,8 +53,8 @@
     print(R.equation)
 print('\n')
 
-gas1 = ct.Solution('gri30.yaml')
-gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics',
+gas1 = ct.Solution(input_file)
+gas2 = ct.Solution(thermo='ideal-gas', kinetics='gas',
                    species=species, reactions=reactions)
 
 
diff --git a/interfaces/cython/cantera/examples/kinetics/reaction_path.py b/interfaces/cython/cantera/examples/kinetics/reaction_path.py
index 629ffda619e..544f84c9000 100644
--- a/interfaces/cython/cantera/examples/kinetics/reaction_path.py
+++ b/interfaces/cython/cantera/examples/kinetics/reaction_path.py
@@ -9,8 +9,9 @@
 Requires: cantera >= 2.5.0
 """
 
-import os
+from subprocess import run
 import sys
+from pathlib import Path
 
 import cantera as ct
 
@@ -33,14 +34,14 @@
 
 dot_file = 'rxnpath.dot'
 img_file = 'rxnpath.png'
-img_path = os.path.join(os.getcwd(), img_file)
+img_path = Path.cwd().joinpath(img_file)
 
 diagram.write_dot(dot_file)
 print(diagram.get_data())
 
-print("Wrote graphviz input file to '{0}'.".format(os.path.join(os.getcwd(), dot_file)))
+print("Wrote graphviz input file to '{0}'.".format(Path.cwd().joinpath(dot_file)))
 
-os.system('dot {0} -Tpng -o{1} -Gdpi=200'.format(dot_file, img_file))
+run('dot {0} -Tpng -o{1} -Gdpi=200'.format(dot_file, img_file).split())
 print("Wrote graphviz output file to '{0}'.".format(img_path))
 
 if "-view" in sys.argv:
diff --git a/interfaces/cython/cantera/examples/onedim/diffusion_flame.py b/interfaces/cython/cantera/examples/onedim/diffusion_flame.py
index e910e74633d..66906e9cfc1 100644
--- a/interfaces/cython/cantera/examples/onedim/diffusion_flame.py
+++ b/interfaces/cython/cantera/examples/onedim/diffusion_flame.py
@@ -41,7 +41,7 @@
 f.oxidizer_inlet.T = tin_o
 
 # Set the boundary emissivities
-f.set_boundary_emissivities(0.0, 0.0)
+f.boundary_emissivities = 0.0, 0.0
 # Turn radiation off
 f.radiation_enabled = False
 
diff --git a/interfaces/cython/cantera/examples/onedim/diffusion_flame_extinction.py b/interfaces/cython/cantera/examples/onedim/diffusion_flame_extinction.py
index 14b84313ad3..ac168b0d7b3 100644
--- a/interfaces/cython/cantera/examples/onedim/diffusion_flame_extinction.py
+++ b/interfaces/cython/cantera/examples/onedim/diffusion_flame_extinction.py
@@ -130,13 +130,12 @@
         f.solve(loglevel=0)
     except ct.CanteraError as e:
         print('Error: Did not converge at n =', n, e)
-    if np.max(f.T) > temperature_limit_extinction:
+    if not np.isclose(np.max(f.T), temperature_limit_extinction):
         # Flame is still burning, so proceed to next strain rate
         n_last_burning = n
         if hdf_output:
             group = 'extinction/{0:04d}'.format(n)
-            f.write_hdf(file_name, group=group, quiet=False,
-                        description='extinction iteration'.format(n))
+            f.write_hdf(file_name, group=group, quiet=True)
         else:
             file_name = 'extinction_{0:04d}.xml'.format(n)
             f.save(os.path.join(data_directory, file_name),
@@ -145,19 +144,32 @@
                    ', reaction mechanism ' + reaction_mechanism)
         T_max.append(np.max(f.T))
         a_max.append(np.max(np.abs(np.gradient(f.velocity) / np.gradient(f.grid))))
+        print('Flame burning at alpha = {:8.4F}. Proceeding to the next iteration, '
+              'with delta_alpha = {}'.format(alpha[-1], delta_alpha))
+    elif ((T_max[-2] - T_max[-1] < delta_T_min) and (delta_alpha < delta_alpha_min)):
         # If the temperature difference is too small and the minimum relative
-        # strain rate increase is reached, abort
-        if ((T_max[-2] - T_max[-1] < delta_T_min) &
-                (delta_alpha < delta_alpha_min)):
-            print('Flame extinguished at n = {0}.'.format(n),
-                  'Abortion criterion satisfied.')
-            break
+        # strain rate increase is reached, save the last, non-burning, solution
+        # to the output file and break the loop
+        T_max.append(np.max(f.T))
+        a_max.append(np.max(np.abs(np.gradient(f.velocity) / np.gradient(f.grid))))
+        if hdf_output:
+            group = 'extinction/{0:04d}'.format(n)
+            f.write_hdf(file_name, group=group, quiet=True)
+        else:
+            file_name = 'extinction_{0:04d}.xml'.format(n)
+            f.save(os.path.join(data_directory, file_name), name='solution', loglevel=0)
+        print('Flame extinguished at alpha = {0:8.4F}.'.format(alpha[-1]),
+              'Abortion criterion satisfied.')
+        break
     else:
         # Procedure if flame extinguished but abortion criterion is not satisfied
-        print('Flame extinguished at n = {0}. Restoring n = {1} with '
-              'alpha = {2}'.format(n, n_last_burning, alpha[n_last_burning]))
         # Reduce relative strain rate increase
         delta_alpha = delta_alpha / delta_alpha_factor
+
+        print('Flame extinguished at alpha = {0:8.4F}. Restoring alpha = {1:8.4F} and '
+              'trying delta_alpha = {2}'.format(
+                  alpha[-1], alpha[n_last_burning], delta_alpha))
+
         # Restore last burning solution
         if hdf_output:
             group = 'extinction/{0:04d}'.format(n_last_burning)
@@ -168,7 +180,15 @@
                       name='solution', loglevel=0)
 
 
-# Print some parameters at the extinction point
+# Print some parameters at the extinction point, after restoring the last burning
+# solution
+if hdf_output:
+    group = 'extinction/{0:04d}'.format(n_last_burning)
+    f.read_hdf(file_name, group=group)
+else:
+    file_name = 'extinction_{0:04d}.xml'.format(n_last_burning)
+    f.restore(os.path.join(data_directory, file_name),
+              name='solution', loglevel=0)
 print('----------------------------------------------------------------------')
 print('Parameters at the extinction point:')
 print('Pressure p={0} bar'.format(f.P / 1e5))
diff --git a/interfaces/cython/cantera/examples/onedim/flame_fixed_T.py b/interfaces/cython/cantera/examples/onedim/flame_fixed_T.py
index 7811da192a6..fd2bc54a7ed 100644
--- a/interfaces/cython/cantera/examples/onedim/flame_fixed_T.py
+++ b/interfaces/cython/cantera/examples/onedim/flame_fixed_T.py
@@ -7,6 +7,7 @@
 
 import cantera as ct
 import numpy as np
+from pathlib import Path
 
 ################################################################
 # parameter values
@@ -38,7 +39,9 @@
 
 # read temperature vs. position data from a file.
 # The file is assumed to have one z, T pair per line, separated by a comma.
-zloc, tvalues = np.genfromtxt('tdata.dat', delimiter=',', comments='#').T
+# The data file must be stored in the same folder as this script.
+data_file = Path(__file__).parent.joinpath('tdata.dat')
+zloc, tvalues = np.genfromtxt(str(data_file), delimiter=',', comments='#').T
 zloc /= max(zloc)
 
 # set the temperature profile to the values read in
diff --git a/interfaces/cython/cantera/examples/onedim/ion_burner_flame.py b/interfaces/cython/cantera/examples/onedim/ion_burner_flame.py
index cb7c744a559..e480dffc91c 100644
--- a/interfaces/cython/cantera/examples/onedim/ion_burner_flame.py
+++ b/interfaces/cython/cantera/examples/onedim/ion_burner_flame.py
@@ -9,7 +9,7 @@
 p = ct.one_atm
 tburner = 600.0
 reactants = 'CH4:1.0, O2:2.0, N2:7.52'  # premixed gas composition
-width = 0.5  # m
+width = 0.05  # m
 loglevel = 1  # amount of diagnostic output (0 to 5)
 
 gas = ct.Solution('gri30_ion.yaml')
diff --git a/interfaces/cython/cantera/examples/reactors/ic_engine.py b/interfaces/cython/cantera/examples/reactors/ic_engine.py
index 13b75c8fefb..c1260cfd0ee 100644
--- a/interfaces/cython/cantera/examples/reactors/ic_engine.py
+++ b/interfaces/cython/cantera/examples/reactors/ic_engine.py
@@ -191,6 +191,7 @@ def ca_ticks(t):
 t = states.t
 
 # pressure and temperature
+xticks = np.arange(0, 0.18, 0.02)
 fig, ax = plt.subplots(nrows=2)
 ax[0].plot(t, states.P / 1.e5)
 ax[0].set_ylabel('$p$ [bar]')
@@ -199,7 +200,8 @@ def ca_ticks(t):
 ax[1].plot(t, states.T)
 ax[1].set_ylabel('$T$ [K]')
 ax[1].set_xlabel(r'$\phi$ [deg]')
-ax[1].set_xticklabels(ca_ticks(ax[1].get_xticks()))
+ax[1].set_xticks(xticks)
+ax[1].set_xticklabels(ca_ticks(xticks))
 plt.show()
 
 # p-V diagram
@@ -224,7 +226,8 @@ def ca_ticks(t):
 ax.legend(loc=0)
 ax.set_ylabel('[kW]')
 ax.set_xlabel(r'$\phi$ [deg]')
-ax.set_xticklabels(ca_ticks(ax.get_xticks()))
+ax.set_xticks(xticks)
+ax.set_xticklabels(ca_ticks(xticks))
 plt.show()
 
 # gas composition
@@ -236,7 +239,8 @@ def ca_ticks(t):
 ax.legend(loc=0)
 ax.set_ylabel('$X_i$ [-]')
 ax.set_xlabel(r'$\phi$ [deg]')
-ax.set_xticklabels(ca_ticks(ax.get_xticks()))
+ax.set_xticks(xticks)
+ax.set_xticklabels(ca_ticks(xticks))
 plt.show()
 
 ######################################################################
@@ -245,22 +249,21 @@ def ca_ticks(t):
 
 # heat release
 Q = trapz(states.heat_release_rate * states.V, t)
-print('{:45s}{:>4.1f} kW'.format('Heat release rate per cylinder (estimate):',
-                                 Q / t[-1] / 1000.))
+output_str = '{:45s}{:>4.1f} {}'
+print(output_str.format('Heat release rate per cylinder (estimate):',
+                        Q / t[-1] / 1000., 'kW'))
 
 # expansion power
 W = trapz(states.dWv_dt, t)
-print('{:45s}{:>4.1f} kW'.format('Expansion power per cylinder (estimate):',
-                                 W / t[-1] / 1000.))
+print(output_str.format('Expansion power per cylinder (estimate):',
+                        W / t[-1] / 1000., 'kW'))
 
 # efficiency
 eta = W / Q
-print('{:45s}{:>4.1f} %'.format('Efficiency (estimate):',
-                                eta * 100.))
+print(output_str.format('Efficiency (estimate):', eta * 100., '%'))
 
 # CO emissions
 MW = states.mean_molecular_weight
 CO_emission = trapz(MW * states.mdot_out * states('CO').X[:, 0], t)
 CO_emission /= trapz(MW * states.mdot_out, t)
-print('{:45s}{:>4.1f} ppm'.format('CO emission (estimate):',
-                                  CO_emission * 1.e6))
+print(output_str.format('CO emission (estimate):', CO_emission * 1.e6, 'ppm'))
diff --git a/interfaces/cython/cantera/examples/reactors/reactor2.py b/interfaces/cython/cantera/examples/reactors/reactor2.py
index 0b33ffe9a56..168136016ab 100644
--- a/interfaces/cython/cantera/examples/reactors/reactor2.py
+++ b/interfaces/cython/cantera/examples/reactors/reactor2.py
@@ -27,8 +27,8 @@
 # First create each gas needed, and a reactor or reservoir for each one.
 
 # create an argon gas object and set its state
-ar = ct.Solution('argon.yaml')
-ar.TP = 1000.0, 20.0 * ct.one_atm
+ar = ct.Solution('air.yaml')
+ar.TPX = 1000.0, 20.0 * ct.one_atm, "AR:1"
 
 # create a reactor to represent the side of the cylinder filled with argon
 r1 = ct.IdealGasReactor(ar)
@@ -113,4 +113,4 @@
     plt.tight_layout()
     plt.show()
 else:
-    print("""To view a plot of these results, run this script with the option -plot""")
+    print("To view a plot of these results, run this script with the option --plot")
diff --git a/interfaces/cython/cantera/examples/surface_chemistry/catalytic_combustion.py b/interfaces/cython/cantera/examples/surface_chemistry/catalytic_combustion.py
index d75e20cabed..8b3566e2164 100644
--- a/interfaces/cython/cantera/examples/surface_chemistry/catalytic_combustion.py
+++ b/interfaces/cython/cantera/examples/surface_chemistry/catalytic_combustion.py
@@ -64,7 +64,7 @@
 # grid
 sim = ct.ImpingingJet(gas=gas, width=width, surface=surf_phase)
 
-# Objects of class StagnationFlow have members that represent the gas inlet
+# Objects of class ImpingingJet have members that represent the gas inlet
 # ('inlet') and the surface ('surface'). Set some parameters of these objects.
 sim.inlet.mdot = mdot
 sim.inlet.T = tinlet
@@ -74,7 +74,7 @@
 # Show the initial solution estimate
 sim.show_solution()
 
-# Solving problems with stiff chemistry coulpled to flow can require a
+# Solving problems with stiff chemistry coupled to flow can require a
 # sequential approach where solutions are first obtained for simpler problems
 # and used as the initial guess for more difficult problems.
 
diff --git a/interfaces/cython/cantera/examples/surface_chemistry/diamond_cvd.py b/interfaces/cython/cantera/examples/surface_chemistry/diamond_cvd.py
index 0b75c08d87e..ecf7d2301b4 100644
--- a/interfaces/cython/cantera/examples/surface_chemistry/diamond_cvd.py
+++ b/interfaces/cython/cantera/examples/surface_chemistry/diamond_cvd.py
@@ -9,7 +9,7 @@
 role in diamond CVD, and this example computes the growth rate and surface
 coverages as a function of [H] at the surface for fixed temperature and [CH3].
 
-Requires: cantera >= 2.5.0, pandas >= 0.21.0, matplotlib >= 2.0
+Requires: cantera >= 2.5.0, pandas >= 0.25.0, matplotlib >= 2.0
 """
 
 import csv
@@ -58,19 +58,13 @@
     import pandas as pd
     data = pd.read_csv('diamond.csv')
 
-    plt.figure()
-    plt.plot(data['H mole Fraction'], data['Growth Rate (microns/hour)'])
+    data.plot(x="H mole Fraction", y="Growth Rate (microns/hour)", legend=False)
     plt.xlabel('H Mole Fraction')
     plt.ylabel('Growth Rate (microns/hr)')
     plt.show()
 
-    plt.figure()
-    for name in data:
-        if name.startswith(('H mole', 'Growth')):
-            continue
-        plt.plot(data['H mole Fraction'], data[name], label=name)
-
-    plt.legend()
+    names = [name for name in data.columns if not name.startswith(('H mole', 'Growth'))]
+    data.plot(x='H mole Fraction', y=names, legend=True)
     plt.xlabel('H Mole Fraction')
     plt.ylabel('Coverage')
     plt.show()
diff --git a/interfaces/cython/cantera/examples/surface_chemistry/sofc.py b/interfaces/cython/cantera/examples/surface_chemistry/sofc.py
index 883dd074eca..5877ff1a0c3 100644
--- a/interfaces/cython/cantera/examples/surface_chemistry/sofc.py
+++ b/interfaces/cython/cantera/examples/surface_chemistry/sofc.py
@@ -1,7 +1,7 @@
 """
 A simple model of a solid oxide fuel cell.
 
-Unlike most SOFC models, this model does not use semi-empirical Butler- Volmer
+Unlike most SOFC models, this model does not use semi-empirical Butler-Volmer
 kinetics for the charge transfer reactions, but uses elementary, reversible
 reactions obeying mass-action kinetics for all reactions, including charge
 transfer. As this script will demonstrate, this approach allows computing the
@@ -120,14 +120,16 @@ def anode_curr(E):
     # get the species net production rates due to the anode-side TPB reaction
     # mechanism. The production rate array has the values for the neighbor
     # species in the order listed in the .yaml file, followed by the tpb phase.
-    # Since the first neighbor phase is the bulk metal, species 0 is the
-    # electron.
+    # The kinetics_species_index method finds the index of the species,
+    # accounting for the possibility of species being in different orders in the
+    # arrays.
     w = tpb_a.net_production_rates
+    electron_index = tpb_a.kinetics_species_index('electron')
 
     # the sign convention is that the current is positive when
-    # electrons are being delivered to the anode - i.e. it is positive
+    # electrons are being delivered to the anode - that is, it is positive
     # for fuel cell operation.
-    return ct.faraday * w[0] * TPB_length_per_area
+    return ct.faraday * w[electron_index] * TPB_length_per_area
 
 
 #####################################################################
@@ -167,13 +169,15 @@ def cathode_curr(E):
     # get the species net production rates due to the cathode-side TPB
     # reaction mechanism. The production rate array has the values for the
     # neighbor species in the order listed in the .yaml file, followed by the
-    # tpb phase. Since the first neighbor phase is the bulk metal, species 0
-    # is the electron.
+    # tpb phase. The kinetics_species_index method finds the index of the species,
+    # accounting for the possibility of species being in different orders in the
+    # arrays.
     w = tpb_c.net_production_rates
+    electron_index = tpb_c.kinetics_species_index('electron')
 
     # the sign convention is that the current is positive when electrons are
-    # being drawn from the cathode (i.e, negative production rate).
-    return -ct.faraday * w[0] * TPB_length_per_area
+    # being drawn from the cathode (that is, negative production rate).
+    return -ct.faraday * w[electron_index] * TPB_length_per_area
 
 
 # initialization
diff --git a/interfaces/cython/cantera/examples/transport/multiprocessing_viscosity.py b/interfaces/cython/cantera/examples/transport/multiprocessing_viscosity.py
index 4f23ae0b318..835a9178896 100644
--- a/interfaces/cython/cantera/examples/transport/multiprocessing_viscosity.py
+++ b/interfaces/cython/cantera/examples/transport/multiprocessing_viscosity.py
@@ -55,15 +55,14 @@ def parallel(mech, predicate, nProcs, nTemps):
     """
     P = ct.one_atm
     X = 'CH4:1.0, O2:1.0, N2:3.76'
-    pool = multiprocessing.Pool(processes=nProcs,
-                                initializer=init_process,
-                                initargs=(mech,))
-
-    y = pool.map(predicate,
-                 zip(itertools.repeat(mech),
-                     np.linspace(300, 900, nTemps),
-                     itertools.repeat(P),
-                     itertools.repeat(X)))
+    with multiprocessing.Pool(
+        processes=nProcs, initializer=init_process, initargs=(mech,)
+    ) as pool:
+        y = pool.map(predicate,
+                     zip(itertools.repeat(mech),
+                         np.linspace(300, 900, nTemps),
+                         itertools.repeat(P),
+                         itertools.repeat(X)))
     return y
 
 
diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py
index a5e2c2e3345..3232f4e1466 100644
--- a/interfaces/cython/cantera/onedim.py
+++ b/interfaces/cython/cantera/onedim.py
@@ -2,8 +2,8 @@
 # at https://cantera.org/license.txt for license and copyright information.
 
 from math import erf
-from os import path
 from email.utils import formatdate
+import warnings
 import numpy as np
 
 from ._cantera import *
@@ -51,7 +51,7 @@ def other_components(self, domain=None):
         if isinstance(dom, Inlet1D):
             return tuple([e for e in self._other
                           if e not in {'grid', 'lambda', 'eField'}])
-        elif isinstance(dom, IdealGasFlow):
+        elif isinstance(dom, (IdealGasFlow, IonFlow)):
             return self._other
         else:
             return ()
@@ -253,7 +253,27 @@ def radiation_enabled(self, enable):
         self.flame.radiation_enabled = enable
 
     def set_boundary_emissivities(self, e_left, e_right):
-        self.flame.set_boundary_emissivities(e_left, e_right)
+        """
+        .. deprecated:: 2.5
+
+             To be deprecated with version 2.5, and removed thereafter.
+             Replaced by property `boundary_emissivities`.
+        """
+        warnings.warn("Method 'set_boundary_emissivities' to be removed after "
+                      "Cantera 2.5. Replaced by property "
+                      "'boundary_emissivities'", DeprecationWarning)
+        self.flame.boundary_emissivities = e_left, e_right
+
+    @property
+    def boundary_emissivities(self):
+        """ Set/get boundary emissivities. """
+        return self.flame.boundary_emissivities
+
+    @boundary_emissivities.setter
+    def boundary_emissivities(self, epsilon):
+        if len(epsilon) != 2:
+            raise ValueError("Boundary emissivities must both be set at the same time.")
+        self.flame.boundary_emissivities = epsilon[0], epsilon[1]
 
     @property
     def grid(self):
diff --git a/interfaces/cython/cantera/onedim.pyx b/interfaces/cython/cantera/onedim.pyx
index 86f81dea246..97de224169a 100644
--- a/interfaces/cython/cantera/onedim.pyx
+++ b/interfaces/cython/cantera/onedim.pyx
@@ -432,6 +432,8 @@ cdef class ReactingSurface1D(Boundary1D):
         """Controls whether or not to solve the surface coverage equations."""
         def __set__(self, value):
             self.surf.enableCoverageEquations(<cbool>value)
+        def __get__(self):
+            return self.surf.coverageEnabled()
 
 
 cdef class _FlowBase(Domain1D):
@@ -578,10 +580,10 @@ cdef class _FlowBase(Domain1D):
         def __get__(self):
             return self.flow.leftEmissivity(), self.flow.rightEmissivity()
         def __set__(self, tuple epsilon):
-            if len(epsilon) == 2:
-                self.flow.setBoundaryEmissivities(epsilon[0], epsilon[1])
-            else:
-                raise ValueError("Setter requires tuple of length 2.")
+            if len(epsilon) != 2:
+                raise ValueError('Setting the boundary emissivities requires a '
+                                 'tuple of length 2.')
+            self.flow.setBoundaryEmissivities(epsilon[0], epsilon[1])
 
     property radiation_enabled:
         """ Determines whether or not to include radiative heat transfer """
diff --git a/interfaces/cython/cantera/reactionpath.pyx b/interfaces/cython/cantera/reactionpath.pyx
index afd1929f8a4..d23f1789843 100644
--- a/interfaces/cython/cantera/reactionpath.pyx
+++ b/interfaces/cython/cantera/reactionpath.pyx
@@ -1,6 +1,8 @@
 # This file is part of Cantera. See License.txt in the top-level directory or
 # at https://cantera.org/license.txt for license and copyright information.
 
+from pathlib import Path
+
 cdef class ReactionPathDiagram:
     def __cinit__(self, *args, **kwargs):
         self._log = new CxxStringStream()
@@ -158,7 +160,7 @@ cdef class ReactionPathDiagram:
         Write the reaction path diagram formatted for use by Graphviz's 'dot'
         program to the file named *filename*.
         """
-        open(filename, 'wb').write(self.get_dot().encode('utf-8'))
+        Path(filename).write_text(self.get_dot())
 
     def get_data(self):
         """
diff --git a/samples/f77/SConscript b/samples/f77/SConscript
index 0cd943c15d3..eb589c89b83 100644
--- a/samples/f77/SConscript
+++ b/samples/f77/SConscript
@@ -10,7 +10,8 @@ samples = [('ctlib', ['ctlib.f']),
            ('isentropic', ['isentropic.f'])]
 
 ftn_demo = localenv.Object('demo_ftnlib.cpp',
-                           CPPPATH=['#include', localenv['boost_inc_dir']])
+                           CPPPATH=['#include', localenv['boost_inc_dir'],
+                                    localenv['extra_inc_dirs']])
 for program_name, fortran_sources in samples:
     buildSample(localenv.Program, program_name,
                 fortran_sources + ftn_demo,
diff --git a/src/oneD/Boundary1D.cpp b/src/oneD/Boundary1D.cpp
index fdb336cc311..170956a8851 100644
--- a/src/oneD/Boundary1D.cpp
+++ b/src/oneD/Boundary1D.cpp
@@ -679,12 +679,10 @@ void ReactingSurf1D::eval(size_t jg, double* xg, double* rg,
     size_t ioffset = m_kin->kineticsSpeciesIndex(0, m_surfindex);
 
     if (m_enabled) {
-        double maxx = -1.0;
         for (size_t k = 0; k < m_nsp; k++) {
             r[k] = m_work[k + ioffset] * m_sphase->size(k) * rs0;
             r[k] -= rdt*(x[k] - prevSoln(k,0));
             diag[k] = 1;
-            maxx = std::max(x[k], maxx);
         }
         r[0] = 1.0 - sum;
         diag[0] = 0;
@@ -707,9 +705,17 @@ void ReactingSurf1D::eval(size_t jg, double* xg, double* rg,
         double* xb = x - nc;
         rb[c_offset_T] = xb[c_offset_T] - m_temp; // specified T
         size_t nSkip = m_flow_left->rightExcessSpecies();
+        size_t l_offset = 0;
+        ThermoPhase* left_thermo = &m_flow_left->phase();
+        for (size_t nth = 0; nth < m_kin->nPhases(); nth++) {
+            if (&m_kin->thermo(nth) == left_thermo) {
+                l_offset = m_kin->kineticsSpeciesIndex(0, nth);
+                break;
+            }
+        }
         for (size_t nl = 0; nl < m_left_nsp; nl++) {
             if (nl != nSkip) {
-                rb[c_offset_Y+nl] += m_work[nl]*mwleft[nl];
+                rb[c_offset_Y+nl] += m_work[nl + l_offset]*mwleft[nl];
             }
         }
     }
diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp
index 85fd70a3520..da3f985312b 100644
--- a/src/zeroD/Reactor.cpp
+++ b/src/zeroD/Reactor.cpp
@@ -296,9 +296,10 @@ double Reactor::evalSurfaces(double t, double* ydot)
         ydot[loc] = sum;
         loc += nk;
 
+        size_t bulkloc = kin->kineticsSpeciesIndex(m_thermo->speciesName(0));
         double wallarea = S->area();
         for (size_t k = 0; k < m_nsp; k++) {
-            m_sdot[k] += m_work[k]*wallarea;
+            m_sdot[k] += m_work[bulkloc + k] * wallarea;
             mdot_surf += m_sdot[k] * mw[k];
         }
     }
diff --git a/test/data/ch4_ion.yaml b/test/data/ch4_ion.yaml
index 617385a5ca0..eaadccf8499 100644
--- a/test/data/ch4_ion.yaml
+++ b/test/data/ch4_ion.yaml
@@ -1,102 +1,186 @@
+generator: cti2yaml
+cantera-version: 2.5.0a4
+date: Mon, 10 Aug 2020 12:55:00 +0000
+input-files: [ch4_ion.cti]
+
 units: {length: cm, quantity: mol, activation-energy: cal/mol}
 
 phases:
 - name: gas
   thermo: ideal-gas
+  elements: [O, H, C, N, E]
+  species:
+  - gri30.yaml/species: [H, O, OH, HO2, H2O2, C, CH, CH2, CH2(S), CH3, HCO,
+      CH2O, CH3O]
+  - species: [H2, O2, H2O, CH4, CO, CO2, N2, HCO+, H3O+, E, O2-]
+  skip-undeclared-third-bodies: true
+  kinetics: gas
+  reactions:
+  - gri30.yaml/reactions: declared-species
+  - reactions: declared-species
   transport: ionized-gas
-  species: [O2, H2O, CO2, N2, H3O+, E, O2-]
-  state: {T: 300, P: 1 atm}
+  state:
+    T: 300.0
+    P: 1.01325e+05
 
 species:
+- name: H2
+  composition: {H: 2}
+  thermo:
+    model: NASA7
+    temperature-ranges: [200.0, 1000.0, 3500.0]
+    data:
+    - [2.34433112, 7.98052075e-03, -1.9478151e-05, 2.01572094e-08, -7.37611761e-12,
+      -917.935173, 0.683010238]
+    - [3.3372792, -4.94024731e-05, 4.99456778e-07, -1.79566394e-10, 2.00255376e-14,
+      -950.158922, -3.20502331]
+  transport:
+    model: gas
+    geometry: linear
+    diameter: 2.92
+    well-depth: 38.0
+    polarizability: 0.455
+    rotational-relaxation: 280.0
+  note: TPIS78
 - name: O2
   composition: {O: 2}
   thermo:
     model: NASA7
-    temperature-ranges: [200, 1000, 3500]
+    temperature-ranges: [200.0, 1000.0, 3500.0]
     data:
-    - [3.78245636E+00, -2.99673416E-03, 9.84730201E-06, -9.68129509E-09,
-       3.24372837E-12, -1.06394356E+03, 3.65767573E+00]
-    - [3.28253784E+00, 1.48308754E-03, -7.57966669E-07, 2.09470555E-10,
-       -2.16717794E-14, -1.08845772E+03, 5.45323129E+00]
+    - [3.78245636, -2.99673416e-03, 9.84730201e-06, -9.68129509e-09, 3.24372837e-12,
+      -1063.94356, 3.65767573]
+    - [3.28253784, 1.48308754e-03, -7.57966669e-07, 2.09470555e-10, -2.16717794e-14,
+      -1088.45772, 5.45323129]
   transport:
     model: gas
     geometry: linear
     diameter: 3.46
-    well-depth: 107.40
+    well-depth: 107.4
     polarizability: 1.131
-    rotational-relaxation: 3.80
-    note: TPIS89
-
+    rotational-relaxation: 3.8
+  note: TPIS89
 - name: H2O
   composition: {H: 2, O: 1}
   thermo:
     model: NASA7
-    temperature-ranges: [200, 1000, 3500]
+    temperature-ranges: [200.0, 1000.0, 3500.0]
     data:
-    - [4.19864056E+00, -2.03643410E-03, 6.52040211E-06, -5.48797062E-09,
-       1.77197817E-12, -3.02937267E+04, -8.49032208E-01]
-    - [3.03399249E+00, 2.17691804E-03, -1.64072518E-07, -9.70419870E-11,
-       1.68200992E-14, -3.00042971E+04, 4.96677010E+00]
+    - [4.19864056, -2.0364341e-03, 6.52040211e-06, -5.48797062e-09, 1.77197817e-12,
+      -3.02937267e+04, -0.849032208]
+    - [3.03399249, 2.17691804e-03, -1.64072518e-07, -9.7041987e-11, 1.68200992e-14,
+      -3.00042971e+04, 4.9667701]
   transport:
     model: gas
     geometry: nonlinear
-    diameter: 2.60
-    well-depth: 572.40
+    diameter: 2.6
+    well-depth: 572.4
     dipole: 1.84
     polarizability: 1.053
-    rotational-relaxation: 4.00
-    note: L 8/89
-
+    rotational-relaxation: 4.0
+  note: L 8/89
+- name: CH4
+  composition: {C: 1, H: 4}
+  thermo:
+    model: NASA7
+    temperature-ranges: [200.0, 1000.0, 3500.0]
+    data:
+    - [5.14987613, -0.0136709788, 4.91800599e-05, -4.84743026e-08, 1.66693956e-11,
+      -1.02466476e+04, -4.64130376]
+    - [0.074851495, 0.0133909467, -5.73285809e-06, 1.22292535e-09, -1.0181523e-13,
+      -9468.34459, 18.437318]
+  transport:
+    model: gas
+    geometry: nonlinear
+    diameter: 3.75
+    well-depth: 141.4
+    polarizability: 2.6
+    rotational-relaxation: 13.0
+  note: L 8/88
+- name: CO
+  composition: {C: 1, O: 1}
+  thermo:
+    model: NASA7
+    temperature-ranges: [200.0, 1000.0, 3500.0]
+    data:
+    - [3.57953347, -6.1035368e-04, 1.01681433e-06, 9.07005884e-10, -9.04424499e-13,
+      -1.4344086e+04, 3.50840928]
+    - [2.71518561, 2.06252743e-03, -9.98825771e-07, 2.30053008e-10, -2.03647716e-14,
+      -1.41518724e+04, 7.81868772]
+  transport:
+    model: gas
+    geometry: linear
+    diameter: 3.65
+    well-depth: 98.1
+    polarizability: 1.95
+    rotational-relaxation: 1.8
+  note: TPIS79
 - name: CO2
   composition: {C: 1, O: 2}
   thermo:
     model: NASA7
-    temperature-ranges: [200, 1000, 3500]
+    temperature-ranges: [200.0, 1000.0, 3500.0]
     data:
-    - [2.35677352E+00, 8.98459677E-03, -7.12356269E-06, 2.45919022E-09,
-       -1.43699548E-13, -4.83719697E+04, 9.90105222E+00]
-    - [3.85746029E+00, 4.41437026E-03, -2.21481404E-06, 5.23490188E-10,
-       -4.72084164E-14, -4.87591660E+04, 2.27163806E+00]
+    - [2.35677352, 8.98459677e-03, -7.12356269e-06, 2.45919022e-09, -1.43699548e-13,
+      -4.83719697e+04, 9.90105222]
+    - [3.85746029, 4.41437026e-03, -2.21481404e-06, 5.23490188e-10, -4.72084164e-14,
+      -4.8759166e+04, 2.27163806]
   transport:
     model: gas
     geometry: linear
     diameter: 3.76
-    well-depth: 244.00
+    well-depth: 244.0
     polarizability: 2.65
-    rotational-relaxation: 2.10
-    note: L 7/88
-
+    rotational-relaxation: 2.1
+  note: L 7/88
 - name: N2
   composition: {N: 2}
   thermo:
     model: NASA7
-    temperature-ranges: [300, 1000, 5000]
+    temperature-ranges: [300.0, 1000.0, 5000.0]
     data:
-    - [3.29867700E+00, 1.40824040E-03, -3.96322200E-06, 5.64151500E-09,
-       -2.44485400E-12, -1.02089990E+03, 3.95037200E+00]
-    - [2.92664000E+00, 1.48797680E-03, -5.68476000E-07, 1.00970380E-10,
-       -6.75335100E-15, -9.22797700E+02, 5.98052800E+00]
+    - [3.298677, 1.4082404e-03, -3.963222e-06, 5.641515e-09, -2.444854e-12,
+      -1020.8999, 3.950372]
+    - [2.92664, 1.4879768e-03, -5.68476e-07, 1.0097038e-10, -6.753351e-15,
+      -922.7977, 5.980528]
   transport:
     model: gas
     geometry: linear
     diameter: 3.62
     well-depth: 97.53
     polarizability: 1.76
-    rotational-relaxation: 4.00
+    rotational-relaxation: 4.0
     dispersion-coefficient: 2.995
     quadrupole-polarizability: 3.602
   note: '121286'
-
+- name: HCO+
+  composition: {H: 1, C: 1, O: 1, E: -1}
+  thermo:
+    model: NASA7
+    temperature-ranges: [300.0, 1000.0, 5000.0]
+    data:
+    - [2.4739736, 8.671559e-03, -1.00315e-05, 6.7170527e-09, -1.7872674e-12,
+      9.9146608e+04, 8.17571187]
+    - [3.741188, 3.3441517e-03, -1.2397121e-06, 2.1189388e-10, -1.370415e-14,
+      9.8884078e+04, 2.07861357]
+  transport:
+    model: gas
+    geometry: linear
+    diameter: 3.59
+    well-depth: 498.0
+    polarizability: 1.356
+    dispersion-coefficient: 0.416
+  note: J12/70
 - name: H3O+
   composition: {H: 3, O: 1, E: -1}
   thermo:
     model: NASA7
-    temperature-ranges: [298.15, 1000, 6000]
+    temperature-ranges: [298.15, 1000.0, 6000.0]
     data:
-    - [3.79295270E+00, -9.10854000E-04, 1.16363549E-05, -1.21364887E-08,
-       4.26159663E-12, 7.07512401E+04, 1.47156856E+00]
-    - [2.49647716E+00, 5.72844920E-03, -1.83953281E-06, 2.73577439E-10,
-       -1.54093985E-14, 7.09729113E+04, 7.45850779E+00]
+    - [3.7929527, -9.10854e-04, 1.16363549e-05, -1.21364887e-08, 4.26159663e-12,
+      7.07512401e+04, 1.47156856]
+    - [2.49647716, 5.7284492e-03, -1.83953281e-06, 2.73577439e-10, -1.54093985e-14,
+      7.09729113e+04, 7.45850779]
   transport:
     model: gas
     geometry: nonlinear
@@ -104,39 +188,50 @@ species:
     well-depth: 106.2
     dipole: 1.417
     polarizability: 0.897
-    rotational-relaxation: 0.0
-    note: TPIS89
-
+  note: TPIS89
 - name: O2-
   composition: {E: 1, O: 2}
   thermo:
     model: NASA7
-    temperature-ranges: [298.15, 1000.00, 6000]
+    temperature-ranges: [298.15, 1000.0, 6000.0]
     data:
-    - [3.66442522E+00, -9.28741138E-04, 6.45477082E-06, -7.7470338E-09,
-       2.93332662E-12, -6.87076983E+03, 4.35140681E+00]
-    - [3.95666294E+00, 5.98141823E-04, -2.12133905E-07, 3.63267581E-11,
-       -2.24989228E-15, -7.06287229E+03, 2.27871017E+00]
+    - [3.66442522, -9.28741138e-04, 6.45477082e-06, -7.7470338e-09, 2.93332662e-12,
+      -6870.76983, 4.35140681]
+    - [3.95666294, 5.98141823e-04, -2.12133905e-07, 3.63267581e-11, -2.24989228e-15,
+      -7062.87229, 2.27871017]
   transport:
     model: gas
     geometry: linear
     diameter: 3.33
     well-depth: 136.5
     polarizability: 1.424
-    note: L4/89
-
+  note: L4/89
 - name: E
   composition: {E: 1}
   thermo:
     model: NASA7
-    temperature-ranges: [200, 3000]
+    temperature-ranges: [200.0, 1000.0, 6000.0]
     data:
-    - [2.5E+00, 0.0, 0.0, 0.0, 0.0, -7.45375000E+02, -1.17246902E+01]
+    - [2.5, 0.0, 0.0, 0.0, 0.0, -745.375, -11.7246902]
+    - [2.5, 0.0, 0.0, 0.0, 0.0, -745.375, -11.7246902]
   transport:
     model: gas
     geometry: atom
     diameter: 2.05
     well-depth: 145.0
     polarizability: 0.667
-    rotational-relaxation: 0.0
-    note: gas L10/92
+  note: gas L10/92
+
+reactions:
+- equation: CH + O => HCO+ + E  # Reaction 1
+  rate-constant: {A: 2.51e+11, b: 0.0, Ea: 1700}
+- equation: HCO+ + H2O => H3O+ + CO  # Reaction 2
+  rate-constant: {A: 1.51e+15, b: 0.0, Ea: 0.0}
+- equation: H3O+ + E => H2O + H  # Reaction 3
+  rate-constant: {A: 2.29e+18, b: -0.5, Ea: 0.0}
+- equation: H3O+ + E => OH + H + H  # Reaction 4
+  rate-constant: {A: 7.95e+21, b: -1.4, Ea: 0.0}
+- equation: H3O+ + E => H2 + OH  # Reaction 5
+  rate-constant: {A: 1.25e+19, b: -0.5, Ea: 0.0}
+- equation: H3O+ + E =>  O + H2 + H  # Reaction 6
+  rate-constant: {A: 6.0e+17, b: -0.3, Ea: 0.0}