diff --git a/CMakeLists.txt b/CMakeLists.txt
index e37fe7311..dc2113b2b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -68,6 +68,8 @@ set(HERMES_SOURCES
     src/component.cxx
     src/component_scheduler.cxx
     src/cx_reaction.cxx
+    src/elastic_collisions.cxx
+    src/hydrogen_molecule_reactions.cxx
     src/ionisation.cxx
     src/izn_rec_reaction.cxx
     src/div_ops.cxx
@@ -132,6 +134,7 @@ set(HERMES_SOURCES
     include/component_scheduler.hxx
     include/diamagnetic_drift.hxx
     include/div_ops.hxx
+    include/elastic_collisions.hxx
     include/electromagnetic.hxx
     include/electron_force_balance.hxx
     include/evolve_density.hxx
@@ -144,10 +147,12 @@ set(HERMES_SOURCES
     include/guarded_options.hxx
     include/neutral_full_velocity.hxx
     include/hermes_utils.hxx
+    include/hydrogen_molecule_reactions.hxx
     include/solkit_hydrogen_charge_exchange.hxx
     include/integrate.hxx
     include/ionisation.hxx
     include/izn_rec_reaction.hxx
+    include/molecular_reactions.hxx
     include/neutral_boundary.hxx
     include/neutral_mixed.hxx
     include/neutral_parallel_diffusion.hxx
diff --git a/hermes-3.cxx b/hermes-3.cxx
index 5a06834b9..baef922cf 100644
--- a/hermes-3.cxx
+++ b/hermes-3.cxx
@@ -51,6 +51,7 @@
 #include "include/fixed_fraction_radiation.hxx"
 #include "include/fixed_temperature.hxx"
 #include "include/fixed_velocity.hxx"
+#include "include/hydrogen_molecule_reactions.hxx"
 #include "include/ionisation.hxx"
 #include "include/isothermal.hxx"
 #include "include/izn_rec_reaction.hxx"
@@ -96,8 +97,8 @@ BOUT_OVERRIDE_DEFAULT_OPTION("mesh:calcParallelSlices_on_communicate", false);
 class DecayLengthBoundary : public BoundaryOp {
 public:
   DecayLengthBoundary() : gen(nullptr) {}
-  DecayLengthBoundary(BoundaryRegion* region,  std::shared_ptr<FieldGenerator> g)
-    : BoundaryOp(region), gen(std::move(g)) {}
+  DecayLengthBoundary(BoundaryRegion* region, std::shared_ptr<FieldGenerator> g)
+      : BoundaryOp(region), gen(std::move(g)) {}
 
   using BoundaryOp::clone;
   /// Create a copy of this boundary condition
@@ -118,10 +119,11 @@ class DecayLengthBoundary : public BoundaryOp {
     ASSERT1(mesh == f.getMesh());
 
     // Get cell radial length
-    Coordinates *coord = mesh->getCoordinates();
+    Coordinates* coord = mesh->getCoordinates();
     Field2D dx = coord->dx;
     Field2D g11 = coord->g11;
-    Field2D dr = dx / sqrt(g11); // cell radial length. dr = dx/(Bpol * R) and g11 = (Bpol*R)**2 
+    Field2D dr =
+        dx / sqrt(g11); // cell radial length. dr = dx/(Bpol * R) and g11 = (Bpol*R)**2
 
     // Only implemented for cell centre quantities
     ASSERT1(f.getLocation() == CELL_CENTRE);
@@ -129,11 +131,13 @@ class DecayLengthBoundary : public BoundaryOp {
     // This loop goes over the first row of boundary cells (in X and Y)
     for (bndry->first(); !bndry->isDone(); bndry->next1d()) {
       for (int zk = 0; zk < mesh->LocalNz; zk++) { // Loop over Z points
-        BoutReal decay_length = 3; // Default decay length is 3 normalised units (usually ~3mm)
+        BoutReal decay_length =
+            3; // Default decay length is 3 normalised units (usually ~3mm)
         if (gen) {
           // Pick up the boundary condition setting from the input file
           // Must be specified in normalised units like the other BCs inputs
-          decay_length = gen->generate(bout::generator::Context(bndry, zk, CELL_CENTRE, 0.0, mesh));
+          decay_length =
+              gen->generate(bout::generator::Context(bndry, zk, CELL_CENTRE, 0.0, mesh));
         }
         // Set value in inner guard cell f(bndry->x, bndry->y, zk)
         // using the final domain cell value f(bndry->x - bndry->bx, bndry->y - bndry->by, zk)
@@ -142,19 +146,21 @@ class DecayLengthBoundary : public BoundaryOp {
         //    (-1, 0) X inner boundary (Core or PF)
         //    (0, 1)  Y upper boundary (outer lower target)
         //    (0, -1) Y lower boundary (inner lower target)
-        
+
         // Distance between final cell centre and inner guard cell centre in normalised units
-        BoutReal distance = 0.5 * (dr(bndry->x, bndry->y) + 
-          dr(bndry->x - bndry->bx, bndry->y - bndry->by));
+        BoutReal distance =
+            0.5
+            * (dr(bndry->x, bndry->y) + dr(bndry->x - bndry->bx, bndry->y - bndry->by));
 
-        // Exponential decay 
-        f(bndry->x, bndry->y, zk) =
-          f(bndry->x - bndry->bx, bndry->y - bndry->by, zk) * exp(-1 * distance / decay_length);
+        // Exponential decay
+        f(bndry->x, bndry->y, zk) = f(bndry->x - bndry->bx, bndry->y - bndry->by, zk)
+                                    * exp(-1 * distance / decay_length);
 
         // Set any remaining guard cells (i.e. the outer guards) to the same value
         // Should the outer guards have the decay continue, or just copy what the inners have?
         for (int i = 1; i < bndry->width; i++) {
-          f(bndry->x + i * bndry->bx, bndry->y + i * bndry->by, zk) = f(bndry->x, bndry->y, zk);
+          f(bndry->x + i * bndry->bx, bndry->y + i * bndry->by, zk) =
+              f(bndry->x, bndry->y, zk);
         }
       }
     }
@@ -170,8 +176,8 @@ class DecayLengthBoundary : public BoundaryOp {
 
 int Hermes::init(bool restarting) {
 
-  auto &options = Options::root()["hermes"];
-  
+  auto& options = Options::root()["hermes"];
+
   output.write("\nGit Version of Hermes: {:s}\n", hermes::version::revision);
   options["revision"] = hermes::version::revision;
   options["revision"].setConditionallyUsed();
@@ -193,7 +199,7 @@ int Hermes::init(bool restarting) {
   units["inv_meters_cubed"] = Nnorm;
   units["eV"] = Tnorm;
   units["Tesla"] = Bnorm;
-  units["seconds"] = 1./Omega_ci;
+  units["seconds"] = 1. / Omega_ci;
   units["meters"] = rho_s0;
 
   // Put into the options tree, so quantities can be normalised
@@ -212,11 +218,13 @@ int Hermes::init(bool restarting) {
   TRACE("Loading metric tensor");
 
   if (options.isSet("loadmetric")) {
-    throw BoutException("Error: The loadmetric option has been renamed to recalculate_metric.\n"
-                        "Note: The default (true/false) is inverted for the new option.\n"
-                        "Setting recalculate_metric=false (the default) uses the metric in the grid file.\n"
-                        "Setting recalculate_metric=true loads Rxy, Bpxy etc from the grid file.\n"
-                        "  This assumes an orthogonal coordinate system. See manual for details.");
+    throw BoutException(
+        "Error: The loadmetric option has been renamed to recalculate_metric.\n"
+        "Note: The default (true/false) is inverted for the new option.\n"
+        "Setting recalculate_metric=false (the default) uses the metric in the grid "
+        "file.\n"
+        "Setting recalculate_metric=true loads Rxy, Bpxy etc from the grid file.\n"
+        "  This assumes an orthogonal coordinate system. See manual for details.");
   }
 
   if (options["recalculate_metric"]
@@ -227,28 +235,23 @@ int Hermes::init(bool restarting) {
   } else {
     // Check that the grid file contains at least one metric tensor component
     // Note: Older grid files did not, so would silently default to the identity metric
-    if (!(mesh->sourceHasVar("J") or
-          mesh->sourceHasVar("g11") or
-          mesh->sourceHasVar("g22") or
-          mesh->sourceHasVar("g33") or
-          mesh->sourceHasVar("g12") or
-          mesh->sourceHasVar("g23") or
-          mesh->sourceHasVar("g13") or
-          mesh->sourceHasVar("g_11") or
-          mesh->sourceHasVar("g_22") or
-          mesh->sourceHasVar("g_33") or
-          mesh->sourceHasVar("g_12") or
-          mesh->sourceHasVar("g_23") or
-          mesh->sourceHasVar("g_13"))) {
-      throw BoutException("Grid input does not contain any metric components (J, g11, g_22 etc).\n"
-                          "Set hermes:recalculate_metric=true to calculate from Rxy, Bpxy etc.\n"
-                          "If the default identity metric is intended then set e.g. mesh:J=1\n");
+    if (!(mesh->sourceHasVar("J") or mesh->sourceHasVar("g11")
+          or mesh->sourceHasVar("g22") or mesh->sourceHasVar("g33")
+          or mesh->sourceHasVar("g12") or mesh->sourceHasVar("g23")
+          or mesh->sourceHasVar("g13") or mesh->sourceHasVar("g_11")
+          or mesh->sourceHasVar("g_22") or mesh->sourceHasVar("g_33")
+          or mesh->sourceHasVar("g_12") or mesh->sourceHasVar("g_23")
+          or mesh->sourceHasVar("g_13"))) {
+      throw BoutException(
+          "Grid input does not contain any metric components (J, g11, g_22 etc).\n"
+          "Set hermes:recalculate_metric=true to calculate from Rxy, Bpxy etc.\n"
+          "If the default identity metric is intended then set e.g. mesh:J=1\n");
     }
 
     if (options["normalise_metric"]
-             .doc("Normalise input metric tensor? (assumes input is in SI units)")
-             .withDefault<bool>(true)) {
-      Coordinates *coord = mesh->getCoordinates();
+            .doc("Normalise input metric tensor? (assumes input is in SI units)")
+            .withDefault<bool>(true)) {
+      Coordinates* coord = mesh->getCoordinates();
       // To use non-orthogonal metric
       // Normalise
       coord->dx /= rho_s0 * rho_s0 * Bnorm;
@@ -296,7 +299,7 @@ int Hermes::init(bool restarting) {
 int Hermes::rhs(BoutReal time, bool linear) {
   // Need to reset the state, since fields may be modified in transform steps
   state = Options();
-  
+
   set(state["time"], time);
   state["units"] = units.copy();
   set(state["linear"], linear);
@@ -330,82 +333,70 @@ void Hermes::outputVars(Options& options) {
   // Save normalisation quantities. These may be used by components
   // to calculate conversion factors to SI units
 
-  set_with_attrs(options["Tnorm"], Tnorm, {
-      {"units", "eV"},
-      {"conversion", 1}, // Already in SI units
-      {"standard_name", "temperature normalisation"},
-      {"long_name", "temperature normalisation"}
-    });
-  set_with_attrs(options["Nnorm"], Nnorm, {
-      {"units", "m^-3"},
-      {"conversion", 1},
-      {"standard_name", "density normalisation"},
-      {"long_name", "Number density normalisation"}
-    });
-  set_with_attrs(options["Bnorm"], Bnorm, {
-      {"units", "T"},
-      {"conversion", 1},
-      {"standard_name", "magnetic field normalisation"},
-      {"long_name", "Magnetic field normalisation"}
-    });
-  set_with_attrs(options["Cs0"], Cs0, {
-      {"units", "m/s"},
-      {"conversion", 1},
-      {"standard_name", "velocity normalisation"},
-      {"long_name", "Sound speed normalisation"}
-    });
-  set_with_attrs(options["Omega_ci"], Omega_ci, {
-      {"units", "s^-1"},
-      {"conversion", 1},
-      {"standard_name", "frequency normalisation"},
-      {"long_name", "Cyclotron frequency normalisation"}
-    });
-  set_with_attrs(options["rho_s0"], rho_s0, {
-      {"units", "m"},
-      {"conversion", 1},
-      {"standard_name", "length normalisation"},
-      {"long_name", "Gyro-radius length normalisation"}
-    });
+  set_with_attrs(options["Tnorm"], Tnorm,
+                 {{"units", "eV"},
+                  {"conversion", 1}, // Already in SI units
+                  {"standard_name", "temperature normalisation"},
+                  {"long_name", "temperature normalisation"}});
+  set_with_attrs(options["Nnorm"], Nnorm,
+                 {{"units", "m^-3"},
+                  {"conversion", 1},
+                  {"standard_name", "density normalisation"},
+                  {"long_name", "Number density normalisation"}});
+  set_with_attrs(options["Bnorm"], Bnorm,
+                 {{"units", "T"},
+                  {"conversion", 1},
+                  {"standard_name", "magnetic field normalisation"},
+                  {"long_name", "Magnetic field normalisation"}});
+  set_with_attrs(options["Cs0"], Cs0,
+                 {{"units", "m/s"},
+                  {"conversion", 1},
+                  {"standard_name", "velocity normalisation"},
+                  {"long_name", "Sound speed normalisation"}});
+  set_with_attrs(options["Omega_ci"], Omega_ci,
+                 {{"units", "s^-1"},
+                  {"conversion", 1},
+                  {"standard_name", "frequency normalisation"},
+                  {"long_name", "Cyclotron frequency normalisation"}});
+  set_with_attrs(options["rho_s0"], rho_s0,
+                 {{"units", "m"},
+                  {"conversion", 1},
+                  {"standard_name", "length normalisation"},
+                  {"long_name", "Gyro-radius length normalisation"}});
   scheduler->outputVars(options);
 }
 
 void Hermes::restartVars(Options& options) {
-  set_with_attrs(options["Tnorm"], Tnorm, {
-      {"units", "eV"},
-      {"conversion", 1}, // Already in SI units
-      {"standard_name", "temperature normalisation"},
-      {"long_name", "temperature normalisation"}
-    });
-  set_with_attrs(options["Nnorm"], Nnorm, {
-      {"units", "m^-3"},
-      {"conversion", 1},
-      {"standard_name", "density normalisation"},
-      {"long_name", "Number density normalisation"}
-    });
-  set_with_attrs(options["Bnorm"], Bnorm, {
-      {"units", "T"},
-      {"conversion", 1},
-      {"standard_name", "magnetic field normalisation"},
-      {"long_name", "Magnetic field normalisation"}
-    });
-  set_with_attrs(options["Cs0"], Cs0, {
-      {"units", "m/s"},
-      {"conversion", 1},
-      {"standard_name", "velocity normalisation"},
-      {"long_name", "Sound speed normalisation"}
-    });
-  set_with_attrs(options["Omega_ci"], Omega_ci, {
-      {"units", "s^-1"},
-      {"conversion", 1},
-      {"standard_name", "frequency normalisation"},
-      {"long_name", "Cyclotron frequency normalisation"}
-    });
-  set_with_attrs(options["rho_s0"], rho_s0, {
-      {"units", "m"},
-      {"conversion", 1},
-      {"standard_name", "length normalisation"},
-      {"long_name", "Gyro-radius length normalisation"}
-    });
+  set_with_attrs(options["Tnorm"], Tnorm,
+                 {{"units", "eV"},
+                  {"conversion", 1}, // Already in SI units
+                  {"standard_name", "temperature normalisation"},
+                  {"long_name", "temperature normalisation"}});
+  set_with_attrs(options["Nnorm"], Nnorm,
+                 {{"units", "m^-3"},
+                  {"conversion", 1},
+                  {"standard_name", "density normalisation"},
+                  {"long_name", "Number density normalisation"}});
+  set_with_attrs(options["Bnorm"], Bnorm,
+                 {{"units", "T"},
+                  {"conversion", 1},
+                  {"standard_name", "magnetic field normalisation"},
+                  {"long_name", "Magnetic field normalisation"}});
+  set_with_attrs(options["Cs0"], Cs0,
+                 {{"units", "m/s"},
+                  {"conversion", 1},
+                  {"standard_name", "velocity normalisation"},
+                  {"long_name", "Sound speed normalisation"}});
+  set_with_attrs(options["Omega_ci"], Omega_ci,
+                 {{"units", "s^-1"},
+                  {"conversion", 1},
+                  {"standard_name", "frequency normalisation"},
+                  {"long_name", "Cyclotron frequency normalisation"}});
+  set_with_attrs(options["rho_s0"], rho_s0,
+                 {{"units", "m"},
+                  {"conversion", 1},
+                  {"standard_name", "length normalisation"},
+                  {"long_name", "Gyro-radius length normalisation"}});
   scheduler->restartVars(options);
 }
 
diff --git a/include/elastic_collisions.hxx b/include/elastic_collisions.hxx
new file mode 100644
index 000000000..ced9f5f38
--- /dev/null
+++ b/include/elastic_collisions.hxx
@@ -0,0 +1,51 @@
+#pragma once
+#ifndef ELASTIC_COLLISIONS_H
+#define ELASTIC_COLLISIONS_H
+
+#include "reaction.hxx"
+#include <string>
+
+namespace hermes {
+
+/**
+ * @brief Base class for elastic collision 'reactions'.
+ */
+struct ElasticCollision : public Reaction {
+  /**
+   * @brief Main constructor for ElasticCollision.
+   *
+   * @details Passes `add_pop_change_sources`=`false` to the Reaction
+   * constructor to disable standard density, momentum, energy sources due to population
+   * changes.
+   *
+   * @param name
+   * @param options  The options object
+   */
+  ElasticCollision(std::string name, Options& options);
+
+  /**
+   * @brief ElasticCollision constructor used by component factory.
+   *
+   * @param name
+   * @param options  The options object
+   * @param solver  The solver object for the simulation (discarded by this class)
+   */
+  ElasticCollision(std::string name, Options& options, Solver*);
+
+  /**
+   * @brief Perform additional transformations specific to elastic collisions.
+   *
+   * @param state
+   * @param rate_calc_results
+   */
+  void transform_additional([[maybe_unused]] GuardedOptions& state,
+                            [[maybe_unused]] const RateData& rate_calc_results) final;
+
+private:
+  /// @brief Names of the two reactant species (== the product species)
+  std::string r1, r2;
+};
+
+} // namespace hermes
+
+#endif // ELASTIC_COLLISIONS_H
diff --git a/include/hydrogen_molecule_reactions.hxx b/include/hydrogen_molecule_reactions.hxx
new file mode 100644
index 000000000..76c16b1bf
--- /dev/null
+++ b/include/hydrogen_molecule_reactions.hxx
@@ -0,0 +1,262 @@
+#pragma once
+#ifndef HYDROGEN_MOLECULE_REACTIONS_H
+#define HYDROGEN_MOLECULE_REACTIONS_H
+
+#include "cx_reaction.hxx"
+#include "elastic_collisions.hxx"
+#include "molecular_reactions.hxx"
+#include <string>
+
+namespace hermes {
+
+enum class HIsotope { h, d, t };
+
+/**
+ * @brief Class to handle charge exchange between Hydrogen (isotope) ions and molecules.
+ *
+ * @details Based on CXReaction.
+ */
+struct MolHCX : public CXReaction {
+  /**
+   * @brief Main constructor for MolHCX.
+   *
+   * @param name
+   * @param options The options object
+   */
+  MolHCX(std::string name, Options& options);
+
+  /**
+   * @brief Constructor used by component factory.
+   *
+   * @param name
+   * @param options The options object
+   * @param solver The solver object for the simulation
+   */
+  MolHCX(std::string name, Options& options, Solver*);
+};
+
+/**
+ * @brief Class to handle dissociation of molecular Hydrogen (and its isotopes).
+ */
+template <HIsotope isotope>
+struct MolHDissociation : public Dissociation {
+  /**
+   * @brief Constructor for MolHDissociation.
+   *
+   * @param name
+   * @param options The options object
+   */
+  MolHDissociation(std::string name, Options& options) : Dissociation(name, options) {}
+
+  /**
+   * @brief Constructor used by component factory.
+   *
+   * @param name
+   * @param options The options object
+   * @param solver The solver object for the simulation
+   */
+  MolHDissociation(std::string name, Options& options, Solver* solver)
+      : Dissociation(name, options, solver) {}
+};
+
+/**
+ * @brief Class to handle dissociative excitation of molecular Hydrogen (and its
+ * isotopes).
+ */
+template <HIsotope isotope>
+struct MolHDissociativeExc : public DissociativeExc {
+  /**
+   * @brief Constructor for MolHDissociativeExc.
+   *
+   * @param name
+   * @param options The options object
+   */
+  MolHDissociativeExc(std::string name, Options& options)
+      : DissociativeExc(name, options) {}
+
+  /**
+   * @brief Constructor used by component factory.
+   *
+   * @param name
+   * @param options The options object
+   * @param solver The solver object for the simulation
+   */
+  MolHDissociativeExc(std::string name, Options& options, Solver* solver)
+      : DissociativeExc(name, options, solver) {}
+};
+
+/**
+ * @brief Class to handle dissociative ionisation of molecular Hydrogen (and its
+ * isotopes).
+ */
+template <HIsotope isotope>
+struct MolHDissociativeIzn : public DissociativeIzn {
+  /**
+   * @brief Constructor for MolHDissociativeIzn.
+   *
+   * @param name
+   * @param options The options object
+   */
+  MolHDissociativeIzn(std::string name, Options& options)
+      : DissociativeIzn(name, options) {}
+
+  /**
+   * @brief Constructor used by component factory.
+   *
+   * @param name
+   * @param options The options object
+   * @param solver The solver object for the simulation
+   */
+  MolHDissociativeIzn(std::string name, Options& options, Solver* solver)
+      : DissociativeIzn(name, options, solver) {}
+};
+
+/**
+ * @brief Class to handle non-dissociative ionisation of molecular Hydrogen (and its
+ * isotopes).
+ */
+template <HIsotope isotope>
+struct MolHNonDissociativeIzn : public NonDissociativeIzn {
+  /**
+   * @brief Constructor for MolHNonDissociativeIzn.
+   *
+   * @param name
+   * @param options The options object
+   */
+  MolHNonDissociativeIzn(std::string name, Options& options)
+      : NonDissociativeIzn(name, options) {}
+
+  /**
+   * @brief Constructor used by component factory.
+   *
+   * @param name
+   * @param options The options object
+   * @param solver The solver object for the simulation
+   */
+  MolHNonDissociativeIzn(std::string name, Options& options, Solver* solver)
+      : NonDissociativeIzn(name, options, solver) {}
+};
+
+/**
+ * @brief Class to handle dissociative recombination of molecular Hydrogen (and its
+ * isotopes).
+ */
+template <HIsotope isotope>
+struct MolHDissociativeRec : public DissociativeRec {
+  /**
+   * @brief Constructor for MolHDissociativeRec.
+   *
+   * @param name
+   * @param options The options object
+   */
+  MolHDissociativeRec(std::string name, Options& options)
+      : DissociativeRec(name, options) {}
+
+  /**
+   * @brief Constructor used by component factory.
+   *
+   * @param name
+   * @param options The options object
+   * @param solver The solver object for the simulation
+   */
+  MolHDissociativeRec(std::string name, Options& options, Solver* solver)
+      : DissociativeRec(name, options, solver) {}
+};
+
+/**
+ * @brief Templated hydrogen molecule elastic collision reaction.
+ *
+ * @details Child class of ElasticCollision, templated on hydrogen isotope.
+ */
+template <HIsotope isotope>
+struct MolHElasticCollision : public ElasticCollision {
+  /**
+   * @brief Constructor for MolHElasticCollision.
+   *
+   * @param name The name of the reaction
+   * @param options The options object
+   */
+  MolHElasticCollision(std::string name, Options& options)
+      : ElasticCollision(name, options) {}
+
+  /**
+   * @brief Constructor used by component factory.
+   *
+   * @param name The name of the reaction
+   * @param options The options object
+   * @param solver The solver object for the simulation
+   */
+  MolHElasticCollision(std::string name, Options& options,
+                       [[maybe_unused]] Solver* solver)
+      : MolHElasticCollision(name, options) {}
+};
+
+} // namespace hermes
+
+// Register molecular Hydrogen reactions
+namespace {
+
+// Non-dissociative ionisation of H isotope molecules
+RegisterComponent<hermes::MolHNonDissociativeIzn<hermes::HIsotope::h>>
+    register_h2_nondiss_izn_h("h2 + e -> h2+ + 2e");
+RegisterComponent<hermes::MolHNonDissociativeIzn<hermes::HIsotope::d>>
+    register_d2_nondiss_izn_d("d2 + e -> d2+ + 2e");
+RegisterComponent<hermes::MolHNonDissociativeIzn<hermes::HIsotope::t>>
+    register_t2_nondiss_izn_t("t2 + e -> t2+ + 2e");
+
+// Dissociation of H isotope molecules
+RegisterComponent<hermes::MolHDissociation<hermes::HIsotope::h>>
+    register_h2_diss("h2 + e -> 2h + e");
+RegisterComponent<hermes::MolHDissociation<hermes::HIsotope::d>>
+    register_d2_diss("d2 + e -> 2d + e");
+RegisterComponent<hermes::MolHDissociation<hermes::HIsotope::t>>
+    register_t2_diss("t2 + e -> 2t + e");
+
+// Dissociative ionisation of H isotope molecules
+RegisterComponent<hermes::MolHDissociativeIzn<hermes::HIsotope::h>>
+    register_h2_diss_izn_h("h2 + e -> h + h+ + 2e");
+RegisterComponent<hermes::MolHDissociativeIzn<hermes::HIsotope::d>>
+    register_d2_diss_izn_d("d2 + e -> d + d+ + 2e");
+RegisterComponent<hermes::MolHDissociativeIzn<hermes::HIsotope::t>>
+    register_t2_diss_izn_t("t2 + e -> t + t+ + 2e");
+
+// Charge exchange between H isotope molecules and H isotope ions
+RegisterComponent<hermes::MolHCX> register_h2_h_cx("h2 + h+ -> h2+ + h");
+RegisterComponent<hermes::MolHCX> register_d2_d_cx("d2 + d+ -> d2+ + d");
+RegisterComponent<hermes::MolHCX> register_t2_t_cx("t2 + t+ -> t2+ + t");
+
+// Dissociative recombination of H isotope molecular ions
+RegisterComponent<hermes::MolHDissociativeRec<hermes::HIsotope::h>>
+    register_h2_plus_rec_h("h2+ + e -> 2h");
+RegisterComponent<hermes::MolHDissociativeRec<hermes::HIsotope::d>>
+    register_d2_plus_rec_d("d2+ + e -> 2d");
+RegisterComponent<hermes::MolHDissociativeRec<hermes::HIsotope::t>>
+    register_t2_plus_rec_t("t2+ + e -> 2t");
+
+// Dissociative excitation of H isotope molecular ions
+RegisterComponent<hermes::MolHDissociativeExc<hermes::HIsotope::h>>
+    register_h2_plus_diss_exc_h("h2+ + e -> h + h+ + e");
+RegisterComponent<hermes::MolHDissociativeExc<hermes::HIsotope::d>>
+    register_d2_plus_diss_exc_d("d2+ + e -> d + d+ + e");
+RegisterComponent<hermes::MolHDissociativeExc<hermes::HIsotope::t>>
+    register_t2_plus_diss_exc_t("t2+ + e -> t + t+ + e");
+
+// Dissociative ionisation of H isotope molecular ions
+RegisterComponent<hermes::MolHDissociativeIzn<hermes::HIsotope::h>>
+    register_h2_plus_diss_izn_h("h2+ + e -> 2h+ + 2e");
+RegisterComponent<hermes::MolHDissociativeIzn<hermes::HIsotope::d>>
+    register_d2_plus_diss_izn_d("d2+ + e -> 2d+ + 2e");
+RegisterComponent<hermes::MolHDissociativeIzn<hermes::HIsotope::t>>
+    register_t2_plus_diss_izn_t("t2+ + e -> 2t+ + 2e");
+
+// Elastic collisions between H isotope molecules and H isotope ions
+RegisterComponent<hermes::MolHElasticCollision<hermes::HIsotope::h>>
+    register_h2_h_elastic("h2 + h+ -> h2+ + h+");
+RegisterComponent<hermes::MolHElasticCollision<hermes::HIsotope::d>>
+    register_d2_d_elastic("d2 + d+ -> d2+ + d+");
+RegisterComponent<hermes::MolHElasticCollision<hermes::HIsotope::t>>
+    register_t2_t_elastic("t2 + t+ -> t2+ + t+");
+
+} // namespace
+
+#endif // HYDROGEN_MOLECULE_REACTIONS_H
diff --git a/include/izn_rec_reaction.hxx b/include/izn_rec_reaction.hxx
index daec56344..086aca8b4 100644
--- a/include/izn_rec_reaction.hxx
+++ b/include/izn_rec_reaction.hxx
@@ -51,6 +51,9 @@ protected:
   /// Name of the (heavy) species with which collision freqs. are associated in the state
   std::string heavy_collfreq_species;
 
+  /// Radiation rate multiplier, extracted from input options
+  BoutReal radiation_multiplier;
+
 private:
   /// Short reaction type string used in diagnostic names ("iz" or "rec")
   const std::string short_reaction_type;
@@ -81,7 +84,7 @@ struct IznReaction : public IznRecReaction {
    * @param solver  The solver object for the simulation (discarded by this class)
    */
   IznReaction(std::string name, Options& options, [[maybe_unused]] Solver* solver)
-      : IznReaction(name, options){};
+      : IznReaction(name, options) {};
 };
 
 /**
@@ -106,7 +109,7 @@ struct RecReaction : public IznRecReaction {
    * @param solver  The solver object for the simulation (discarded by this class)
    */
   RecReaction(std::string name, Options& options, [[maybe_unused]] Solver* solver)
-      : RecReaction(name, options){};
+      : RecReaction(name, options) {};
 };
 
 } // namespace hermes
diff --git a/include/molecular_reactions.hxx b/include/molecular_reactions.hxx
new file mode 100644
index 000000000..61d576def
--- /dev/null
+++ b/include/molecular_reactions.hxx
@@ -0,0 +1,130 @@
+#pragma once
+#ifndef MOLECULAR_REACTIONS_H
+#define MOLECULAR_REACTIONS_H
+
+#include "cx_reaction.hxx"
+#include "reaction.hxx"
+#include <string>
+
+namespace hermes {
+
+/**
+ * @brief Base class for dissociation reactions involving molecules.
+ *
+ * @details Inherits from Reaction to handle dissociation of molecular species.
+ */
+struct Dissociation : public Reaction {
+  /**
+   * @brief Main constructor for Dissociation.
+   *
+   * @param name
+   * @param options The options object
+   */
+  Dissociation(std::string name, Options& options) : Reaction(name, options) {};
+
+  /**
+   * @brief Constructor used by component factory.
+   *
+   * @param name
+   * @param options The options object
+   * @param solver The solver object for the simulation
+   */
+  Dissociation(std::string name, Options& options, [[maybe_unused]] Solver* solver)
+      : Dissociation(name, options) {};
+};
+
+/**
+ * @brief Base class for dissociative excitation reactions involving molecules.
+ */
+struct DissociativeExc : public Reaction {
+  /**
+   * @brief Main constructor for DissociativeExc.
+   *
+   * @param name
+   * @param options The options object
+   */
+  DissociativeExc(std::string name, Options& options) : Reaction(name, options) {};
+
+  /**
+   * @brief Constructor used by component factory.
+   *
+   * @param name
+   * @param options The options object
+   * @param solver The solver object for the simulation
+   */
+  DissociativeExc(std::string name, Options& options, [[maybe_unused]] Solver* solver)
+      : DissociativeExc(name, options) {};
+};
+
+/**
+ * @brief Base class for dissociative ionization reactions involving molecules.
+ */
+struct DissociativeIzn : public Reaction {
+  /**
+   * @brief Main constructor for DissociativeIzn.
+   *
+   * @param name
+   * @param options The options object
+   */
+  DissociativeIzn(std::string name, Options& options) : Reaction(name, options) {};
+
+  /**
+   * @brief Constructor used by component factory.
+   *
+   * @param name
+   * @param options The options object
+   * @param solver The solver object for the simulation
+   */
+  DissociativeIzn(std::string name, Options& options, [[maybe_unused]] Solver* solver)
+      : DissociativeIzn(name, options) {};
+};
+
+/**
+ * @brief Base class for non-dissociative ionization reactions involving molecules.
+ */
+struct NonDissociativeIzn : public Reaction {
+  /**
+   * @brief Main constructor for NonDissociativeIzn.
+   *
+   * @param name
+   * @param options The options object
+   */
+  NonDissociativeIzn(std::string name, Options& options) : Reaction(name, options) {};
+
+  /**
+   * @brief Constructor used by component factory.
+   *
+   * @param name
+   * @param options The options object
+   * @param solver The solver object for the simulation
+   */
+  NonDissociativeIzn(std::string name, Options& options, [[maybe_unused]] Solver* solver)
+      : NonDissociativeIzn(name, options) {};
+};
+
+/**
+ * @brief Base class for dissociative recombination reactions involving molecules.
+ */
+struct DissociativeRec : public Reaction {
+  /**
+   * @brief Main constructor for DissociativeRec.
+   *
+   * @param name
+   * @param options The options object
+   */
+  DissociativeRec(std::string name, Options& options) : Reaction(name, options) {};
+
+  /**
+   * @brief Constructor used by component factory.
+   *
+   * @param name
+   * @param options The options object
+   * @param solver The solver object for the simulation
+   */
+  DissociativeRec(std::string name, Options& options, [[maybe_unused]] Solver* solver)
+      : DissociativeRec(name, options) {};
+};
+
+} // namespace hermes
+
+#endif // MOLECULAR_REACTIONS_H
diff --git a/include/reaction.hxx b/include/reaction.hxx
index 9f9a36460..5adb4193e 100644
--- a/include/reaction.hxx
+++ b/include/reaction.hxx
@@ -60,8 +60,10 @@ struct Reaction : public ReactionBase {
    *
    * @param name
    * @param options Options object
+   * @param add_pop_change_sources Whether to add population change sources (default:
+   * true)
    */
-  Reaction(std::string name, Options& options);
+  Reaction(std::string name, Options& options, bool add_pop_change_sources = true);
 
   /**
    * @brief Copy all diagnostics into the output, setting the appropriate metadata at the
@@ -83,8 +85,8 @@ protected:
   Options& units;
   BoutReal Tnorm, Nnorm, FreqNorm;
 
-  /// Rate multipliers, extracted from input options
-  BoutReal rate_multiplier, radiation_multiplier;
+  /// Rate multiplier. Defaults to 1, but subclasses may read an option and overwrite.
+  BoutReal rate_multiplier;
 
   /// Output diagnostics?
   bool diagnose;
@@ -210,6 +212,8 @@ protected:
   }
 
 private:
+  const bool add_pop_change_sources;
+
   // Channels to determine how momentum and energy are distributed to product species
   std::map<std::string, std::map<std::string, BoutReal>> energy_channels;
   std::map<std::string, std::map<std::string, BoutReal>> momentum_channels;
diff --git a/include/reaction_settings.hxx b/include/reaction_settings.hxx
index 2f2cd3eee..8eece0357 100644
--- a/include/reaction_settings.hxx
+++ b/include/reaction_settings.hxx
@@ -58,31 +58,72 @@ static inline void add_default_id(ReactionDataTypes data_type,
 static inline void generate_default_data_ids_map() {
   constexpr char H_isotopes[] = {'h', 'd', 't'};
 
-  // H isotope izn
+  // Various atomic and molecular reactions involving H isotopes
+  std::string reaction_str;
   for (const char& isotope : H_isotopes) {
-    std::string reaction_str = fmt::format("{} + e -> {}+ + 2e", isotope, isotope);
+    reaction_str = fmt::format("{} + e -> {}+ + 2e", isotope, isotope);
+    // Ionisation of atoms
     add_default_id(ReactionDataTypes::Amjuel, reaction_str, ReactionCoeffTypes::sigma_v,
                    "H.4_2.1.5");
     add_default_id(ReactionDataTypes::Amjuel, reaction_str, ReactionCoeffTypes::sigma_v_E,
                    "H.10_2.1.5");
-  }
 
-  // H isotope rec
-  for (const char& isotope : H_isotopes) {
-    std::string reaction_str = fmt::format("{}+ + e -> {}", isotope, isotope);
+    // Recombination of ions
+    reaction_str = fmt::format("{}+ + e -> {}", isotope, isotope);
     add_default_id(ReactionDataTypes::Amjuel, reaction_str, ReactionCoeffTypes::sigma_v,
                    "H.4_2.1.8");
     add_default_id(ReactionDataTypes::Amjuel, reaction_str, ReactionCoeffTypes::sigma_v_E,
                    "H.10_2.1.8");
+
+    // Dissociation of neutral molecules
+    reaction_str = fmt::format("{}2 + e -> 2{} + e", isotope, isotope);
+    add_default_id(ReactionDataTypes::Amjuel, reaction_str, ReactionCoeffTypes::sigma_v,
+                   "H.4_2.2.5g");
+
+    // Non-dissociative ionisation of neutral molecules
+    reaction_str = fmt::format("{}2 + e -> {}2+ + 2e", isotope, isotope);
+    add_default_id(ReactionDataTypes::Amjuel, reaction_str, ReactionCoeffTypes::sigma_v,
+                   "H.4_2.2.9");
+
+    // Dissociative ionisation of neutral molecules
+    reaction_str = fmt::format("{}2 + e -> {} + {}+ + 2e", isotope, isotope, isotope);
+    add_default_id(ReactionDataTypes::Amjuel, reaction_str, ReactionCoeffTypes::sigma_v,
+                   "H.4_2.2.10");
+
+    // Dissociative ionisation of (singly-) ionised molecules
+    reaction_str = fmt::format("{}2+ + e -> 2{}+ + 2e", isotope, isotope, isotope);
+    add_default_id(ReactionDataTypes::Amjuel, reaction_str, ReactionCoeffTypes::sigma_v,
+                   "H.4_2.2.11");
+
+    // Dissociative excitation of (singly-) ionised molecules
+    reaction_str = fmt::format("{}2+ + e -> {} + {}+ + e", isotope, isotope, isotope);
+    add_default_id(ReactionDataTypes::Amjuel, reaction_str, ReactionCoeffTypes::sigma_v,
+                   "H.4_2.2.12");
+
+    // Dissociative recombination of (singly-) ionised molecules
+    reaction_str = fmt::format("{}2+ + e -> 2{}", isotope, isotope);
+    add_default_id(ReactionDataTypes::Amjuel, reaction_str, ReactionCoeffTypes::sigma_v,
+                   "H.4_2.2.14");
   }
 
-  // H isotope CX
+  // CX and elastic collisions involving H isotopes
   for (const char& isotope1 : H_isotopes) {
     for (const char& isotope2 : H_isotopes) {
+      // Atom-ion CX
       std::string reaction_str =
           fmt::format("{} + {}+ -> {}+ + {}", isotope1, isotope2, isotope1, isotope2);
       add_default_id(ReactionDataTypes::Amjuel, reaction_str, ReactionCoeffTypes::sigma_v,
                      "H.2_3.1.8");
+      // Molecule-ion CX
+      reaction_str =
+          fmt::format("{}2 + {}+ -> {}2+ + {}", isotope1, isotope2, isotope1, isotope2);
+      add_default_id(ReactionDataTypes::Amjuel, reaction_str, ReactionCoeffTypes::sigma_v,
+                     "H.2_3.2.3");
+      // Molecule-ion elastic collisions
+      reaction_str =
+          fmt::format("{}2 + {}+ -> {}2 + {}+", isotope1, isotope2, isotope1, isotope2);
+      add_default_id(ReactionDataTypes::Amjuel, reaction_str, ReactionCoeffTypes::sigma_v,
+                     "H.2_0.3T");
     }
   }
 
diff --git a/include/species_parser.hxx b/include/species_parser.hxx
index ac8a11338..79d2a2a36 100644
--- a/include/species_parser.hxx
+++ b/include/species_parser.hxx
@@ -11,11 +11,12 @@
 class SpeciesParser {
 public:
   /**
-   * @brief Construct a new SpeciesParser object by extracting the element and charge from
-   * a species string.
+   * @brief Construct a new SpeciesParser object by extracting the base species and charge
+   * from a species string.
    * @details Valid string requirements:
-   * - Element name must be 1 or 2 letters and is stored in lower case.
-   * - Any integer can precede the element name, but is discarded
+   * - Base species (atom, molecule, or electron) consists of 1 or 2 letters optionally
+   * followed by an integer, and is stored in lower case. e.g. "he" or "h2"
+   * - Any integer can *precede* the base species, but is discarded
    *     e.g. "2he" is parsed as "he"
    * - A "+" or a "-", optionally followed by an integer is interpreted as the charge.
    *     e.g. "he-1" (-1), "H" (0) or "ne+8" (+8)
@@ -27,15 +28,14 @@ public:
   SpeciesParser(const std::string& species_str);
 
   int get_charge() const { return this->charge; }
-  std::string get_element() const { return this->element; }
+  std::string get_base_species() const { return this->base_species; }
   std::string get_str() const { return this->species_str; }
 
   /**
-   * @brief Get an appropriate long name for the element in this species (e.g. "hydrogen"
-   * for "h").
+   * @brief Get an appropriate long name for the base species (e.g. "hydrogen" for "h").
    *
-   * @return std::string The long name of the species, or the element name if no long name
-   * is defined
+   * @return std::string The long name of the species, or the base species name if no long
+   * name is defined
    */
   std::string long_name() const;
 
@@ -54,13 +54,13 @@ public:
   SpeciesParser recombined();
 
 private:
-  SpeciesParser(const std::string element, int charge);
+  SpeciesParser(const std::string base_species, int charge);
 
   /// Species string
   std::string species_str;
 
-  /// Atomic element of the species (or 'e' for electrons)
-  std::string element;
+  /// Base species form: atom, molecule, or electron (unneutralized)
+  std::string base_species;
 
   /// Charge of the species
   int charge;
diff --git a/json_database/AMJUEL_H.2_0.3T.json b/json_database/AMJUEL_H.2_0.3T.json
new file mode 100644
index 000000000..45f41d3c2
--- /dev/null
+++ b/json_database/AMJUEL_H.2_0.3T.json
@@ -0,0 +1,21 @@
+{
+    "info": {
+        "fit_type": "T",
+        "ref": "Amjuel H.2 0.3T",
+        "reaction_lbl": "p + H2 -> p + H2 (total rate coefficient)",
+        "notes": "Total rate coefficient for elastic collisions between H2 and H+. Maxwellian rate coefficient vs. T_p, with H2 at rest, obtained by taking the corresponding Beam-Maxw. rate coefficient at Eb=0.06 eV and verification by independent integration of cross-section"
+    },
+    "coeffs": [
+        [
+            -1.857812000000E+01,
+            2.411708000000E-01,
+            1.046088000000E-02,
+            -1.203649000000E-02,
+            -3.679626000000E-03,
+            2.895358000000E-04,
+            1.354441000000E-04,
+            -2.712317000000E-06,
+            -1.356528000000E-06
+        ]
+    ]
+}
diff --git a/json_database/AMJUEL_H.2_3.2.3.json b/json_database/AMJUEL_H.2_3.2.3.json
new file mode 100644
index 000000000..184f072ff
--- /dev/null
+++ b/json_database/AMJUEL_H.2_3.2.3.json
@@ -0,0 +1,21 @@
+{
+    "info": {
+        "fit_type": "T",
+        "ref": "Amjuel H.2 3.2.3",
+        "reaction_lbl": "p + H2 -> H + H2+",
+        "notes": "H2/H+ charge exchange. Effective ion conversion rate (charge exchange on H2) v eff = v H2(v=0) + v=114 v H2(v) pH2(v) Same vibrational distribution (as function of T_e) as above. Therefore: single parameter fit vs. T_e, since vibrational distribution does not depend upon density, E_0 is fixed (0.1 eV) and Tp = Te = T."
+    },
+    "coeffs": [
+        [
+            -2.163099643422E+01,
+            3.206843053514E+00,
+            -3.369939911269E+00,
+            1.290238400703E+00,
+            -3.988189754178E-01,
+            1.462287796966E-01,
+            -3.524154596754E-02,
+            4.146324082808E-03,
+            -1.846022446828E-04
+        ]
+    ]
+}
diff --git a/json_database/AMJUEL_H.4_2.2.10.json b/json_database/AMJUEL_H.4_2.2.10.json
new file mode 100644
index 000000000..535a6c77a
--- /dev/null
+++ b/json_database/AMJUEL_H.4_2.2.10.json
@@ -0,0 +1,109 @@
+{
+    "info": {
+        "fit_type": "nT",
+        "ref": "Amjuel H.4 2.2.10",
+        "reaction_lbl": "e + H2 -> 2e + H + H+",
+        "notes": "Ionization and dissociation of H2"
+    },
+    "coeffs": [
+        [
+            -3.793749300315E+01,
+            -3.333162972531E-01,
+            1.849601203843E-01,
+            -8.803945197107E-02,
+            2.205180180735E-02,
+            -2.852568161901E-03,
+            1.942314738448E-04,
+            -6.597388255594E-06,
+            8.798544848606E-08
+        ],
+        [
+            1.280249398154E+01,
+            1.028969438485E+00,
+            -3.271855492638E-01,
+            1.305597441611E-01,
+            -3.408439821910E-02,
+            4.591924060066E-03,
+            -3.167471002157E-04,
+            1.070920193931E-05,
+            -1.408139742113E-07
+        ],
+        [
+            -3.778148553140E+00,
+            -1.415561059533E+00,
+            2.928509524911E-01,
+            -7.425165688158E-02,
+            2.028424685287E-02,
+            -3.042376564749E-03,
+            2.279124955373E-04,
+            -8.197224564797E-06,
+            1.130682076163E-07
+        ],
+        [
+            2.499987501522E-01,
+            1.032922656537E+00,
+            -1.580288004759E-01,
+            9.934702707539E-03,
+            -2.450845732158E-03,
+            5.716646876513E-04,
+            -5.339115778704E-05,
+            2.135848413694E-06,
+            -3.072223247387E-08
+        ],
+        [
+            2.480574522949E-01,
+            -4.372934216955E-01,
+            6.448433196301E-02,
+            1.229222932630E-03,
+            -9.281410519553E-04,
+            5.946235618034E-05,
+            -8.758032156912E-08,
+            -7.270955072707E-08,
+            1.100087131523E-09
+        ],
+        [
+            -9.960628182831E-02,
+            1.092652428162E-01,
+            -1.782307798975E-02,
+            1.192181214757E-04,
+            2.310636556641E-04,
+            -2.492990725967E-05,
+            1.217600444191E-06,
+            -3.624263301602E-08,
+            6.139167092128E-10
+        ],
+        [
+            1.709129400742E-02,
+            -1.574889001363E-02,
+            2.865310743302E-03,
+            -1.700396064727E-04,
+            -1.502644504654E-06,
+            3.297869416435E-07,
+            6.572135289627E-10,
+            4.269190108005E-10,
+            -3.666090917669E-11
+        ],
+        [
+            -1.435304503973E-03,
+            1.203823111704E-03,
+            -2.350465388313E-04,
+            2.507288189894E-05,
+            -3.077975735212E-06,
+            3.748299687254E-07,
+            -2.613600078122E-08,
+            8.263175463927E-10,
+            -8.509179497022E-12
+        ],
+        [
+            4.808639828229E-05,
+            -3.761591649539E-05,
+            7.490531472388E-06,
+            -1.077314971617E-06,
+            1.950247963978E-07,
+            -2.569729600929E-08,
+            1.804377780165E-09,
+            -6.031847199601E-11,
+            7.416020205748E-13
+        ]
+    ]
+}
diff --git a/json_database/AMJUEL_H.4_2.2.11.json b/json_database/AMJUEL_H.4_2.2.11.json
new file mode 100644
index 000000000..b400510cc
--- /dev/null
+++ b/json_database/AMJUEL_H.4_2.2.11.json
@@ -0,0 +1,109 @@
+{
+    "info": {
+        "fit_type": "nT",
+        "ref": "Amjuel H.4 2.2.11",
+        "reaction_lbl": "e + H2+ -> 2e + H+ + H+",
+        "notes": "Dissociative ionization of H2+"
+    },
+    "coeffs": [
+        [
+            -3.708803769397E+01,
+            9.784233987341E-02,
+            -7.200361272130E-03,
+            6.496843022778E-03,
+            -1.420590818760E-03,
+            1.703620321164E-04,
+            -1.160738946400E-05,
+            4.148222302162E-07,
+            -6.007853385325E-09
+        ],
+        [
+            1.561780529774E+01,
+            -1.673256230592E-02,
+            2.743322772895E-02,
+            -1.026956102747E-02,
+            1.999561527383E-03,
+            -2.043607814503E-04,
+            1.084177127603E-05,
+            -2.671800995803E-07,
+            2.093182411476E-09
+        ],
+        [
+            -6.874406034117E+00,
+            -7.782929961315E-03,
+            -6.888773684846E-03,
+            2.306107197863E-03,
+            -4.029222834436E-04,
+            3.932152471491E-05,
+            -2.094907364150E-06,
+            5.682907060010E-08,
+            -6.320752545610E-10
+        ],
+        [
+            2.010540060675E+00,
+            -3.226785148562E-03,
+            -6.181192193854E-03,
+            2.388146990238E-03,
+            -5.018901320009E-04,
+            5.520233512352E-05,
+            -3.080798536641E-06,
+            7.864770315002E-08,
+            -6.357395371638E-10
+        ],
+        [
+            -3.614768906120E-01,
+            3.710098881765E-03,
+            2.045814599796E-03,
+            -8.523935993991E-04,
+            1.751295192861E-04,
+            -1.944203941844E-05,
+            1.138888354831E-06,
+            -3.256303793266E-08,
+            3.501794038444E-10
+        ],
+        [
+            2.956861321735E-02,
+            -5.524443504504E-04,
+            -2.457951062112E-05,
+            3.433179945503E-05,
+            -1.450208898992E-06,
+            -2.447566480782E-07,
+            1.375679100044E-08,
+            4.863880510459E-10,
+            -3.004374374556E-11
+        ],
+        [
+            9.662490252868E-04,
+            -1.548556801431E-04,
+            1.417215042439E-05,
+            -6.444863591678E-06,
+            -1.566028729499E-06,
+            4.152486680818E-07,
+            -2.855068942744E-08,
+            6.081804811000E-10,
+            9.512865901179E-13
+        ],
+        [
+            -3.543571865464E-04,
+            4.662969089421E-05,
+            -1.471117766355E-05,
+            5.235585096328E-06,
+            -5.779667826854E-07,
+            2.139729421817E-08,
+            -3.656048425230E-10,
+            3.759866326965E-11,
+            -1.486151370215E-12
+        ],
+        [
+            1.827109843671E-05,
+            -3.179895716088E-06,
+            1.432429412413E-06,
+            -5.141065080107E-07,
+            7.734387173369E-08,
+            -6.163336831045E-09,
+            3.128313515842E-10,
+            -1.061842444216E-11,
+            1.771099769640E-13
+        ]
+    ]
+}
diff --git a/json_database/AMJUEL_H.4_2.2.12.json b/json_database/AMJUEL_H.4_2.2.12.json
new file mode 100644
index 000000000..e70f226cd
--- /dev/null
+++ b/json_database/AMJUEL_H.4_2.2.12.json
@@ -0,0 +1,109 @@
+{
+    "info": {
+        "fit_type": "nT",
+        "ref": "Amjuel H.4 2.2.12",
+        "reaction_lbl": "e + H2+ -> e + H + H+",
+        "notes": "Dissociative excitation of H2+. Effective dissociative excitation of H2+ to H and H+, including the component e+H2+ -> H + H* -> H + H+ from dissociative recombination of H2+ with excited products."
+    },
+    "coeffs": [
+        [
+            -1.793443274600E+01,
+            -4.932783688604E-02,
+            1.039088280849E-01,
+            -4.375935166008E-02,
+            9.196691651936E-03,
+            -1.043378648769E-03,
+            6.600342421838E-05,
+            -2.198466460165E-06,
+            3.004145701249E-08
+        ],
+        [
+            2.236108757681E+00,
+            -2.545406018621E-02,
+            -1.160421006835E-01,
+            4.407846563362E-02,
+            -8.192521304984E-03,
+            8.200277386433E-04,
+            -4.508284363534E-05,
+            1.282824614809E-06,
+            -1.474719350236E-08
+        ],
+        [
+            -3.620018994703E-01,
+            6.721527680150E-02,
+            1.564387124002E-02,
+            -4.939045440424E-03,
+            4.263195867947E-04,
+            1.034216805418E-05,
+            -3.975028601900E-06,
+            2.322116289258E-07,
+            -4.381217154470E-09
+        ],
+        [
+            -4.353922258965E-01,
+            -3.051033606589E-02,
+            3.512861172521E-02,
+            -1.179504564265E-02,
+            2.091772760029E-03,
+            -1.991100044575E-04,
+            1.018080238045E-05,
+            -2.597941866088E-07,
+            2.524118386011E-09
+        ],
+        [
+            1.580381801957E-01,
+            2.493654957203E-03,
+            -1.601970998119E-02,
+            5.346709597939E-03,
+            -8.711870134835E-04,
+            7.542066727545E-05,
+            -3.410778344979E-06,
+            7.120460603822E-08,
+            -4.412295474522E-10
+        ],
+        [
+            1.697880687685E-02,
+            2.106675963900E-03,
+            4.521983358170E-04,
+            -3.017151690655E-04,
+            6.209239389357E-05,
+            -7.598119096817E-06,
+            5.523273241689E-07,
+            -2.130508249251E-08,
+            3.319099650589E-10
+        ],
+        [
+            -1.521914651109E-02,
+            -7.527862162788E-04,
+            9.095551479381E-04,
+            -2.372576223034E-04,
+            3.018561480848E-05,
+            -1.365255868731E-06,
+            -4.604769733903E-08,
+            5.867910270430E-09,
+            -1.357779142836E-10
+        ],
+        [
+            2.406276368070E-03,
+            9.971361856278E-05,
+            -1.760978402353E-04,
+            4.877659148871E-05,
+            -6.477358351729E-06,
+            3.541106430252E-07,
+            1.309772899670E-09,
+            -8.072907334230E-10,
+            2.074669430611E-11
+        ],
+        [
+            -1.219469579955E-04,
+            -4.785505675232E-06,
+            9.858840337511E-06,
+            -2.779210878533E-06,
+            3.720379996058E-07,
+            -2.110289928486E-08,
+            3.753875073646E-11,
+            4.024906665497E-11,
+            -1.075990572574E-12
+        ]
+    ]
+}
diff --git a/json_database/AMJUEL_H.4_2.2.14.json b/json_database/AMJUEL_H.4_2.2.14.json
new file mode 100644
index 000000000..b0947320c
--- /dev/null
+++ b/json_database/AMJUEL_H.4_2.2.14.json
@@ -0,0 +1,109 @@
+{
+    "info": {
+        "fit_type": "nT",
+        "ref": "Amjuel H.4 2.2.14",
+        "reaction_lbl": "e + H2+ -> H + H",
+        "notes": "Dissociative recombination of H2+. Effective dissociative recombination of H2+ to H and H, subtracting the ionising component e+H2+ -> H + H* -> H + H+ from dissociative recombination of H2+ with excited products."
+    },
+    "coeffs": [
+        [
+            -1.664335253647E+01,
+            8.953780953631E-02,
+            -1.056411030518E-01,
+            4.477000808690E-02,
+            -9.729945434357E-03,
+            1.174456882002E-03,
+            -7.987743820637E-05,
+            2.842957892768E-06,
+            -4.104508608435E-08
+        ],
+        [
+            -6.005444031657E-01,
+            4.063933992726E-02,
+            -4.753947846841E-02,
+            2.188304031377E-02,
+            -5.201085606791E-03,
+            6.866340394051E-04,
+            -5.059940013116E-05,
+            1.930213882205E-06,
+            -2.963966822809E-08
+        ],
+        [
+            4.494812032769E-04,
+            7.884508616595E-05,
+            3.688007562485E-04,
+            -4.659255785539E-04,
+            1.907115980400E-04,
+            -3.434324710145E-05,
+            3.067651560323E-06,
+            -1.325689465590E-07,
+            2.212493073620E-09
+        ],
+        [
+            1.632894866655E-04,
+            3.108116177617E-04,
+            -3.521552580917E-04,
+            -2.233169775063E-04,
+            1.869415236037E-04,
+            -4.329991211511E-05,
+            4.465256901322E-06,
+            -2.136296167564E-07,
+            3.873085368404E-09
+        ],
+        [
+            -7.234142549752E-05,
+            -1.316311320262E-03,
+            1.643509328764E-03,
+            -6.412764282779E-04,
+            1.048891053765E-04,
+            -7.018555173322E-06,
+            4.776213235854E-08,
+            1.380537343974E-08,
+            -4.199397846492E-10
+        ],
+        [
+            -1.504085050039E-05,
+            1.315865970237E-04,
+            -1.025653773999E-04,
+            5.310324781249E-05,
+            -1.831888048039E-05,
+            3.423755373077E-06,
+            -3.303384352061E-07,
+            1.551627097700E-08,
+            -2.809391819541E-10
+        ],
+        [
+            1.113923667684E-05,
+            2.711411525392E-05,
+            -8.495922363727E-05,
+            4.026487801017E-05,
+            -6.289324474240E-06,
+            1.911447036702E-07,
+            3.638198230235E-08,
+            -3.235540606394E-09,
+            7.605442050634E-11
+        ],
+        [
+            -1.843926162250E-06,
+            -1.663674537499E-06,
+            1.308069926896E-05,
+            -7.324021449032E-06,
+            1.431739868187E-06,
+            -1.085644779665E-07,
+            1.143164983367E-09,
+            2.151595003971E-10,
+            -7.052562220005E-12
+        ],
+        [
+            9.864173150662E-08,
+            -2.212261708468E-07,
+            -4.431749501051E-07,
+            3.270530731011E-07,
+            -7.282085521177E-08,
+            6.578253567957E-09,
+            -1.925258267827E-10,
+            -4.217474167519E-12,
+            2.364754029318E-13
+        ]
+    ]
+}
diff --git a/json_database/AMJUEL_H.4_2.2.5g.json b/json_database/AMJUEL_H.4_2.2.5g.json
new file mode 100644
index 000000000..89ee2a89b
--- /dev/null
+++ b/json_database/AMJUEL_H.4_2.2.5g.json
@@ -0,0 +1,109 @@
+{
+    "info": {
+        "fit_type": "nT",
+        "ref": "Amjuel H.4 2.2.5g",
+        "reaction_lbl": "e + H2 -> e + H + H",
+        "notes": "Dissociation of H2 with vibrational distribution. H2 multi-step model, data: Sawada/Fujimoto/Greenland. H(1), H2, H2+, H+ transported (slow species). The H2(v) are in vibrational equilibrium (depends only upon T_e), and the electronic levels in the molecules as discussed in (Sawada and Fujimoto, 1995) are taken into account. CX losses from vibr. distribution are computed assuming Te = Ti, ne = np, and an energy of H2-Beam = 0.1 eV."
+    },
+    "coeffs": [
+        [
+            -2.702372540584E+01,
+            -3.152103191633E-03,
+            5.990692171729E-03,
+            -3.151252835426E-03,
+            7.457309144890E-04,
+            -9.238664007853E-05,
+            6.222557542845E-06,
+            -2.160024578659E-07,
+            3.028755759836E-09
+        ],
+        [
+            1.081756417479E+01,
+            -1.487216964825E-02,
+            1.417396532101E-02,
+            -4.689911797083E-03,
+            7.180338663163E-04,
+            -5.502798587526E-05,
+            1.983066081752E-06,
+            -2.207639762507E-08,
+            -2.116339335271E-10
+        ],
+        [
+            -5.368872027676E+00,
+            5.419787589654E-03,
+            -1.747268613395E-02,
+            9.532963297450E-03,
+            -2.196705622859E-03,
+            2.611447288152E-04,
+            -1.695536960581E-05,
+            5.737375510694E-07,
+            -7.940900078995E-09
+        ],
+        [
+            1.340684229143E+00,
+            1.058157580038E-02,
+            -3.446019122786E-03,
+            -7.032769815599E-04,
+            4.427959286553E-04,
+            -7.370484189164E-05,
+            5.746786010618E-06,
+            -2.182085196303E-07,
+            3.264045809897E-09
+        ],
+        [
+            -1.561644923145E-01,
+            -3.847438570333E-03,
+            3.571477356851E-03,
+            -1.103305795473E-03,
+            1.476712517858E-04,
+            -8.461162952132E-06,
+            9.757111870171E-08,
+            8.130014050833E-09,
+            -2.234996157750E-10
+        ],
+        [
+            -1.444731533894E-04,
+            -3.194532513126E-04,
+            -2.987368098475E-04,
+            2.092094838648E-04,
+            -4.339352509941E-05,
+            4.009328699469E-06,
+            -1.762651912129E-07,
+            3.357860444624E-09,
+            -1.857322587267E-11
+        ],
+        [
+            2.117693926546E-03,
+            2.679309814780E-04,
+            -1.037559373832E-04,
+            7.297053580368E-06,
+            1.454171585421E-06,
+            -2.251616910293E-07,
+            9.191700327811E-09,
+            -2.052366968228E-11,
+            -3.567738654108E-12
+        ],
+        [
+            -2.143738340207E-04,
+            -3.539232757385E-05,
+            1.909399233821E-05,
+            -3.819368125069E-06,
+            3.754063159414E-07,
+            -2.441872829462E-08,
+            1.437490161488E-09,
+            -6.172308568891E-11,
+            1.104905484620E-12
+        ],
+        [
+            6.979740947331E-06,
+            1.462031952352E-06,
+            -8.858634506391E-07,
+            2.099830142707E-07,
+            -2.606862169776E-08,
+            2.039813579349E-09,
+            -1.113483084607E-10,
+            3.859777100010E-12,
+            -5.909099891913E-14
+        ]
+    ]
+}
diff --git a/json_database/AMJUEL_H.4_2.2.9.json b/json_database/AMJUEL_H.4_2.2.9.json
new file mode 100644
index 000000000..f90d36974
--- /dev/null
+++ b/json_database/AMJUEL_H.4_2.2.9.json
@@ -0,0 +1,109 @@
+{
+    "info": {
+        "fit_type": "nT",
+        "ref": "Amjuel H.4 2.2.9",
+        "reaction_lbl": "e + H2 -> 2e + H2+",
+        "notes": "Ionization of H2"
+    },
+    "coeffs": [
+        [
+            -3.574773783577E+01,
+            3.470247049909E-01,
+            -9.683166540937E-02,
+            1.959576276250E-03,
+            2.479361119190E-03,
+            -1.196632952666E-04,
+            -1.862956119592E-05,
+            1.669867158509E-06,
+            -3.673736278200E-08
+        ],
+        [
+            1.769208985507E+01,
+            -1.311169841222E+00,
+            4.700486215943E-01,
+            -5.521175478827E-02,
+            -2.689651616933E-03,
+            7.308915874002E-04,
+            -2.920560755694E-05,
+            -3.148831240316E-07,
+            2.514856386324E-08
+        ],
+        [
+            -8.291764008409E+00,
+            1.591701525694E+00,
+            -5.814996025336E-01,
+            9.160898084105E-02,
+            -4.770789631868E-03,
+            1.994775632224E-05,
+            -7.511552245648E-06,
+            1.089689676313E-06,
+            -2.920863498031E-08
+        ],
+        [
+            2.555712347240E+00,
+            -8.625268584825E-01,
+            2.612076696684E-01,
+            -3.686525285376E-02,
+            1.945480608139E-03,
+            -3.690918356665E-05,
+            4.836340453567E-06,
+            -4.165748666929E-07,
+            9.265898224345E-09
+        ],
+        [
+            -5.370404654062E-01,
+            2.375816996323E-01,
+            -4.165908778170E-02,
+            1.732469114063E-03,
+            3.693513203529E-04,
+            -4.931268184607E-05,
+            2.727501534044E-06,
+            -1.081027384449E-07,
+            2.420509440644E-09
+        ],
+        [
+            7.443307905391E-02,
+            -3.322214182214E-02,
+            -2.351235556666E-03,
+            1.723053881691E-03,
+            -2.096625925098E-04,
+            1.358575558294E-05,
+            -1.041586202167E-06,
+            6.928574330531E-08,
+            -1.746656185835E-09
+        ],
+        [
+            -6.391785721973E-03,
+            1.862554278190E-03,
+            1.540632467396E-03,
+            -3.547150770477E-04,
+            1.392157055273E-05,
+            1.047463944093E-06,
+            1.513510667993E-08,
+            -9.915499708242E-09,
+            3.298173891188E-10
+        ],
+        [
+            3.001729098239E-04,
+            3.497202259366E-05,
+            -1.742029226138E-04,
+            2.296551698214E-05,
+            2.357520372192E-06,
+            -5.306085513950E-07,
+            2.223137028418E-08,
+            3.340169309800E-10,
+            -2.560542889504E-11
+        ],
+        [
+            -5.607182991432E-06,
+            -5.779550092391E-06,
+            6.495742927455E-06,
+            -3.040011333889E-07,
+            -2.361542565281E-07,
+            3.655056080262E-08,
+            -1.771478792301E-09,
+            1.334615260635E-11,
+            6.831564719957E-13
+        ]
+    ]
+}
diff --git a/src/elastic_collisions.cxx b/src/elastic_collisions.cxx
new file mode 100644
index 000000000..b5759afa5
--- /dev/null
+++ b/src/elastic_collisions.cxx
@@ -0,0 +1,66 @@
+#include "elastic_collisions.hxx"
+
+namespace hermes {
+
+///
+ElasticCollision::ElasticCollision(std::string name, Options& options)
+    : Reaction(name, options, false) {
+  std::vector<std::string> reactants =
+      this->parser->get_species(species_filter::reactants);
+  std::vector<std::string> products = this->parser->get_species(species_filter::products);
+
+  // Must have 2 reactants
+  if (reactants.size() != 2) {
+    throw BoutException(fmt::format("Expected exactly two reactant species in "
+                                    "elastic collision reaction, but found {}: {}",
+                                    reactants.size(), fmt::join(reactants, ", ")));
+  }
+
+  // Check that reactants and products match (order-independent)
+  std::sort(reactants.begin(), reactants.end());
+  std::sort(products.begin(), products.end());
+  ASSERT2(reactants == products);
+  if (reactants != products) {
+    throw BoutException(fmt::format("Reactant and product species don't match in "
+                                    "elastic collision reaction: reactants are {} but "
+                                    "products are {}",
+                                    fmt::join(reactants, ", "),
+                                    fmt::join(products, ", ")));
+  }
+
+  // Store reactant names for convenience
+  this->r1 = reactants[0];
+  this->r2 = reactants[1];
+
+  // Diagnostics
+  if (this->diagnose) {
+    add_diagnostic(
+        this->r1, fmt::format("F{}{}", this->r1, this->r2),
+        fmt::format("Momentum exchange due to elastic collisions between {} and {}",
+                    this->r1, this->r2),
+        ReactionDiagnosticType::momentum_src, this->rate_data->src_str());
+  }
+}
+
+///
+void ElasticCollision::transform_additional(
+    [[maybe_unused]] GuardedOptions& state,
+    [[maybe_unused]] const RateData& rate_calc_results) {
+
+  // Frictional force between the two species
+  const BoutReal& m1 = state["species"][this->r1]["AA"].GetRef<BoutReal>();
+  const BoutReal& m2 = state["species"][this->r2]["AA"].GetRef<BoutReal>();
+  const Field3D& v1 = state["species"][this->r1]["velocity"].GetRef<Field3D>();
+  const Field3D& v2 = state["species"][this->r2]["velocity"].GetRef<Field3D>();
+
+  Field3D mu = (m1 * m2) / (m1 + m2);
+  Field3D momentum_exchange = rate_calc_results.rate * mu * (v1 - v2);
+
+  // Subtract from r1 and diagnostic
+  update_source<subtract<Field3D>>(state, this->r1, ReactionDiagnosticType::momentum_src,
+                                   momentum_exchange);
+  // Add to r2
+  add(state["species"][this->r2]["momentum_source"], momentum_exchange);
+}
+
+} // namespace hermes
diff --git a/src/hydrogen_molecule_reactions.cxx b/src/hydrogen_molecule_reactions.cxx
new file mode 100644
index 000000000..ac6bd5937
--- /dev/null
+++ b/src/hydrogen_molecule_reactions.cxx
@@ -0,0 +1,11 @@
+#include "hydrogen_molecule_reactions.hxx"
+
+namespace hermes {
+
+// MolHCX implementation
+MolHCX::MolHCX(std::string name, Options& options) : CXReaction(name, options) {}
+
+MolHCX::MolHCX(std::string name, Options& options, [[maybe_unused]] Solver* solver)
+    : MolHCX(name, options) {}
+
+} // namespace hermes
diff --git a/src/izn_rec_reaction.cxx b/src/izn_rec_reaction.cxx
index 580cd94c6..c809ce661 100644
--- a/src/izn_rec_reaction.cxx
+++ b/src/izn_rec_reaction.cxx
@@ -18,6 +18,9 @@ IznRecReaction::IznRecReaction(std::string short_reaction_type, std::string name
   this->heavy_product =
       this->parser->get_single_species(species_filter::heavy, species_filter::products);
 
+  // Multiplier defaults to 1, subclasses may read appropriate an option and overwrite
+  this->radiation_multiplier = 1;
+
   /*
   Data for calculating the electron energy loss rate.
   Data type is assumed to be the same as the izn/rec rate data for now and the data ID is
@@ -226,4 +229,4 @@ RecReaction::RecReaction(std::string name, Options& options)
           .withDefault<BoutReal>(1.0);
 }
 
-} // namespace hermes
\ No newline at end of file
+} // namespace hermes
diff --git a/src/reaction.cxx b/src/reaction.cxx
index ebebf9aa7..0d0a2d95e 100644
--- a/src/reaction.cxx
+++ b/src/reaction.cxx
@@ -14,10 +14,11 @@
 namespace hermes {
 
 ///
-Reaction::Reaction(std::string name, Options& options)
+Reaction::Reaction(std::string name, Options& options, bool add_pop_change_sources)
     : ReactionBase({readOnly("species:{sp}:{r_val}"), readOnly("species:e:{e_val}"),
                     readWrite("species:{sp}:{w_val}")}),
-      units(options["units"]), name(name) {
+      units(options["units"]), add_pop_change_sources(add_pop_change_sources),
+      name(name) {
 
   // Extract some relevant options, units to member vars for readability
   this->Tnorm = get<BoutReal>(this->units["eV"]);
@@ -28,6 +29,9 @@ Reaction::Reaction(std::string name, Options& options)
                        .doc("Output additional diagnostics?")
                        .withDefault<bool>(false);
 
+  // Multiplier defaults to 1, subclasses may read appropriate an option and overwrite
+  this->rate_multiplier = 1.0;
+
   std::string reaction_str, data_src_id;
   ReactionDataTypes data_src_type;
 
@@ -189,12 +193,21 @@ void Reaction::init_channel_weights(GuardedOptions& state) {
     double total_energy_weight = std::accumulate(
         this->energy_channels[reactant].begin(), this->energy_channels[reactant].end(),
         0.0, [](double sum, const auto& pair) { return sum + pair.second; });
-    ASSERT0(total_energy_weight >= 0 && total_energy_weight <= 1);
+    if (total_energy_weight < 0 || total_energy_weight > 1.0) {
+      throw BoutException(fmt::format("Total energy channel weight for reactant '{}' is "
+                                      "{} (must satisfy 0 <= weight <= 1).",
+                                      reactant, total_energy_weight));
+    }
     double total_momentum_weight =
         std::accumulate(this->momentum_channels[reactant].begin(),
                         this->momentum_channels[reactant].end(), 0.0,
                         [](double sum, const auto& pair) { return sum + pair.second; });
-    ASSERT0(total_momentum_weight >= 0 && total_momentum_weight <= 1);
+    if (total_momentum_weight < 0 || total_momentum_weight > 1.0) {
+      throw BoutException(
+          fmt::format("Total momentum channel weight for reactant '{}' is "
+                      "{} (must satisfy 0 <= weight <= 1).",
+                      reactant, total_momentum_weight));
+    }
   }
 }
 
@@ -277,63 +290,67 @@ void Reaction::transform_impl(GuardedOptions& state) {
   // Subclasses perform any additional transform tasks
   transform_additional(state, rate_calc_results);
 
-  // Use the stoichiometric values to set density sources for all species
-  Field3D density_source(0.0);
-  for (const auto& sp_name : this->parser->get_species()) {
-    int pop_change = this->parser->pop_change(sp_name);
-    if (pop_change != 0) {
-      // Density sources
-      density_source = pfactors.at(sp_name) * pop_change * rate_calc_results.rate;
-      update_source<add<Field3D>>(state, sp_name, ReactionDiagnosticType::density_src,
-                                  density_source);
+  if (this->add_pop_change_sources) {
+
+    // Use the stoichiometric values to set density sources for all species
+    Field3D density_source(0.0);
+    for (const auto& sp_name : this->parser->get_species()) {
+      int pop_change = this->parser->pop_change(sp_name);
+      if (pop_change != 0) {
+        // Density sources
+        density_source = pfactors.at(sp_name) * pop_change * rate_calc_results.rate;
+        update_source<add<Field3D>>(state, sp_name, ReactionDiagnosticType::density_src,
+                                    density_source);
+      }
     }
-  }
 
-  // Population change-driven sources for all species other than electrons
-  init_channel_weights(state);
-  Field3D momentum_source, energy_source;
-  for (const auto& [sp_name, pop_change_s] : this->parser->get_mom_energy_pop_changes()) {
-    // No momentum, energy source for electrons due to pop change
-    if (sp_name.compare("e") == 0) {
-      continue;
-    }
-    momentum_source = 0.0;
-    energy_source = 0.0;
-    if (pop_change_s < 0) {
-      // For species with net loss, sources follows directly from pop change
-      momentum_source = pop_change_s * rate_calc_results.rate
-                        * get<BoutReal>(state["species"][sp_name]["AA"])
-                        * get<Field3D>(state["species"][sp_name]["velocity"]);
-      energy_source = pop_change_s * rate_calc_results.rate * (3. / 2)
-                      * get<Field3D>(state["species"][sp_name]["temperature"]);
-    } else if (pop_change_s > 0) {
-      // Species with net gain receive a proportion of the momentum and energy lost by
-      // consumed reactants. See init_channel_weights() for default splitting factors.
-      for (auto& rsp_name : heavy_reactant_species) {
-        // All consumed (net loss) reactants can contribute
-        int pop_change_r = this->parser->pop_change_reactant(rsp_name);
-        if (pop_change_r < 0) {
-          momentum_source += -pop_change_r * pfactors.at(rsp_name)
-                             * this->momentum_channels[rsp_name][sp_name]
-                             * rate_calc_results.rate
-                             * get<BoutReal>(state["species"][rsp_name]["AA"])
-                             * get<Field3D>(state["species"][rsp_name]["velocity"]);
-          energy_source += -pop_change_r * pfactors.at(rsp_name)
-                           * this->energy_channels[rsp_name][sp_name]
-                           * rate_calc_results.rate * (3. / 2)
-                           * get<Field3D>(state["species"][rsp_name]["temperature"]);
+    // Population change-driven sources for all species other than electrons
+    init_channel_weights(state);
+    Field3D momentum_source, energy_source;
+    for (const auto& [sp_name, pop_change_s] :
+         this->parser->get_mom_energy_pop_changes()) {
+      // No momentum, energy source for electrons due to pop change
+      if (sp_name.compare("e") == 0) {
+        continue;
+      }
+      momentum_source = 0.0;
+      energy_source = 0.0;
+      if (pop_change_s < 0) {
+        // For species with net loss, sources follows directly from pop change
+        momentum_source = pop_change_s * rate_calc_results.rate
+                          * get<BoutReal>(state["species"][sp_name]["AA"])
+                          * get<Field3D>(state["species"][sp_name]["velocity"]);
+        energy_source = pop_change_s * rate_calc_results.rate * (3. / 2)
+                        * get<Field3D>(state["species"][sp_name]["temperature"]);
+      } else if (pop_change_s > 0) {
+        // Species with net gain receive a proportion of the momentum and energy lost by
+        // consumed reactants. See init_channel_weights() for default splitting factors.
+        for (auto& rsp_name : heavy_reactant_species) {
+          // All consumed (net loss) reactants can contribute
+          int pop_change_r = this->parser->pop_change_reactant(rsp_name);
+          if (pop_change_r < 0) {
+            momentum_source += -pop_change_r * pfactors.at(rsp_name)
+                               * this->momentum_channels[rsp_name][sp_name]
+                               * rate_calc_results.rate
+                               * get<BoutReal>(state["species"][rsp_name]["AA"])
+                               * get<Field3D>(state["species"][rsp_name]["velocity"]);
+            energy_source += -pop_change_r * pfactors.at(rsp_name)
+                             * this->energy_channels[rsp_name][sp_name]
+                             * rate_calc_results.rate * (3. / 2)
+                             * get<Field3D>(state["species"][rsp_name]["temperature"]);
+          }
         }
+      } else {
+        // No pop change
+        continue;
       }
-    } else {
-      // No pop change
-      continue;
-    }
 
-    // Update sources
-    update_source<add<Field3D>>(state, sp_name, ReactionDiagnosticType::momentum_src,
-                                momentum_source);
-    update_source<add<Field3D>>(state, sp_name, ReactionDiagnosticType::energy_src,
-                                energy_source);
+      // Update sources
+      update_source<add<Field3D>>(state, sp_name, ReactionDiagnosticType::momentum_src,
+                                  momentum_source);
+      update_source<add<Field3D>>(state, sp_name, ReactionDiagnosticType::energy_src,
+                                  energy_source);
+    }
   }
 }
 
diff --git a/src/species_parser.cxx b/src/species_parser.cxx
index 430e75648..56146448d 100644
--- a/src/species_parser.cxx
+++ b/src/species_parser.cxx
@@ -10,8 +10,8 @@ static std::map<std::string, std::string> long_names = {
     {"t", "tritium"},  {"he", "helium"},
 };
 
-std::string construct_species_str(std::string element, int charge) {
-  if (element == "e") {
+std::string construct_species_str(std::string base_species, int charge) {
+  if (base_species == "e") {
     // Special case for electrons
     if (charge != -1) {
       throw BoutException(
@@ -38,7 +38,7 @@ std::string construct_species_str(std::string element, int charge) {
         charge_str = "-" + std::to_string(-charge);
       }
     }
-    return element + charge_str;
+    return base_species + charge_str;
   }
 }
 
@@ -47,21 +47,21 @@ std::string construct_species_str(std::string element, int charge) {
 ///
 SpeciesParser::SpeciesParser(const std::string& species_str) {
 
-  // Extract element name, charge and ionisation state with regex
-  // Any number preceding the element is discarded
-  std::regex pattern("^([0-9]*)([a-zA-Z]{1,2})([\\+|\\-]?)([0-9]*)$");
+  // Extract base species, charge and ionisation state with regex
+  // Any number preceding the base species is discarded
+  std::regex pattern("^([0-9]*)([a-zA-Z]{1,2}[0-9]*)([\\+|\\-]?)([0-9]*)$");
   std::smatch matches;
   bool has_matches = std::regex_search(species_str, matches, pattern);
-  // String must provide at least an element name
+  // String must provide at least a base species
   if (!has_matches || !matches[1].matched) {
     throw BoutException(
         fmt::format("Unable to extract charge from species name {}", species_str));
   }
 
-  // Store element name; always lower case
-  this->element = matches[2];
-  std::transform(this->element.begin(), this->element.end(), this->element.begin(),
-                 ::tolower);
+  // Store base species; always lower case
+  this->base_species = matches[2];
+  std::transform(this->base_species.begin(), this->base_species.end(),
+                 this->base_species.begin(), ::tolower);
 
   // Extract charge, electron is a special case
   if (species_str == "e") {
@@ -73,39 +73,39 @@ SpeciesParser::SpeciesParser(const std::string& species_str) {
   }
 
   // Stored species string discards any leading number and is always lower case
-  this->species_str = construct_species_str(this->element, this->charge);
+  this->species_str = construct_species_str(this->base_species, this->charge);
 }
 
 ///
-SpeciesParser::SpeciesParser(const std::string element, int charge)
-    : element(element), charge(charge) {
-  this->species_str = construct_species_str(element, charge);
+SpeciesParser::SpeciesParser(const std::string base_species, int charge)
+    : base_species(base_species), charge(charge) {
+  this->species_str = construct_species_str(base_species, charge);
 }
 
 ///
 std::string SpeciesParser::long_name() const {
-  auto it = long_names.find(this->element);
+  auto it = long_names.find(this->base_species);
   if (it != long_names.end()) {
     return it->second;
   } else {
-    return this->element;
+    return this->base_species;
   }
 }
 
 ///
 SpeciesParser SpeciesParser::ionised() {
-  if (this->element == "e") {
+  if (this->base_species == "e") {
     throw BoutException("Cannot change electron charge!");
   }
   int new_charge = this->charge + 1;
-  return SpeciesParser(this->element, new_charge);
+  return SpeciesParser(this->base_species, new_charge);
 }
 
 ///
 SpeciesParser SpeciesParser::recombined() {
-  if (this->element == "e") {
+  if (this->base_species == "e") {
     throw BoutException("Cannot change electron charge!");
   }
   int new_charge = this->charge - 1;
-  return SpeciesParser(this->element, new_charge);
+  return SpeciesParser(this->base_species, new_charge);
 }
diff --git a/tests/unit/reactions/D2DissIzn.nc b/tests/unit/reactions/D2DissIzn.nc
new file mode 100644
index 000000000..72b25ba3c
Binary files /dev/null and b/tests/unit/reactions/D2DissIzn.nc differ
diff --git a/tests/unit/reactions/D2DpCX.nc b/tests/unit/reactions/D2DpCX.nc
new file mode 100644
index 000000000..f5f1c7446
Binary files /dev/null and b/tests/unit/reactions/D2DpCX.nc differ
diff --git a/tests/unit/reactions/H2Diss.nc b/tests/unit/reactions/H2Diss.nc
new file mode 100644
index 000000000..41dca7f98
Binary files /dev/null and b/tests/unit/reactions/H2Diss.nc differ
diff --git a/tests/unit/reactions/H2HpElastic.nc b/tests/unit/reactions/H2HpElastic.nc
new file mode 100644
index 000000000..6889a7176
Binary files /dev/null and b/tests/unit/reactions/H2HpElastic.nc differ
diff --git a/tests/unit/reactions/H2NonDissIzn.nc b/tests/unit/reactions/H2NonDissIzn.nc
new file mode 100644
index 000000000..51c79cfa4
Binary files /dev/null and b/tests/unit/reactions/H2NonDissIzn.nc differ
diff --git a/tests/unit/reactions/T2pDissExc.nc b/tests/unit/reactions/T2pDissExc.nc
new file mode 100644
index 000000000..2d026f855
Binary files /dev/null and b/tests/unit/reactions/T2pDissExc.nc differ
diff --git a/tests/unit/reactions/T2pDissRec.nc b/tests/unit/reactions/T2pDissRec.nc
new file mode 100644
index 000000000..cfe6da9b4
Binary files /dev/null and b/tests/unit/reactions/T2pDissRec.nc differ
diff --git a/tests/unit/test_molecule_reactions.hxx b/tests/unit/test_molecule_reactions.hxx
new file mode 100644
index 000000000..69d797fc3
--- /dev/null
+++ b/tests/unit/test_molecule_reactions.hxx
@@ -0,0 +1,154 @@
+#pragma once
+#ifndef TEST_MOLECULE_REACTIONS_H
+#define TEST_MOLECULE_REACTIONS_H
+
+#include "hydrogen_molecule_reactions.hxx"
+
+namespace hermes {
+
+/**
+ * @brief Class to test dissociative reactions.
+ *
+ */
+
+template <typename RTYPE>
+class DissociativeReactionTest : public ReactionTest<RTYPE> {
+
+  static_assert(std::is_base_of<Reaction, RTYPE>(),
+                "Template arg to DissociativeReactionTest must derive from Reaction");
+
+protected:
+  DissociativeReactionTest(std::string lbl, std::string reaction_str)
+      : ReactionTest<RTYPE>(lbl, reaction_str) {
+
+    // Generate sorted lists of heavy reactants and products
+    this->heavy_reactants =
+        this->parser.get_species(species_filter::reactants, species_filter::heavy);
+    std::sort(this->heavy_reactants.begin(), this->heavy_reactants.end());
+    this->heavy_products =
+        this->parser.get_species(species_filter::products, species_filter::heavy);
+    std::sort(this->heavy_products.begin(), this->heavy_products.end());
+    this->heavy_species = this->parser.get_species(species_filter::heavy);
+    std::sort(this->heavy_species.begin(), this->heavy_species.end());
+  };
+
+  /// Sorted lists of heavy reactants, products, all species, for use in generating the
+  /// test state.
+  std::vector<std::string> heavy_reactants;
+  std::vector<std::string> heavy_products;
+  std::vector<std::string> heavy_species;
+
+  /**
+   * @brief Generated state should work with any number of heavy reactants/products.
+   *
+   * @return Options
+   */
+  Options generate_state() final {
+    // Must have at least one heavy reactant and one heavy product.
+    if (heavy_reactants.size() == 0 || heavy_products.size() == 0) {
+      throw BoutException(
+          fmt::format("DissociativeReactionTest::generate_state expects reactions with "
+                      "at least one heavy reactant and one heavy product."));
+    }
+
+    std::string comp_name("test" + this->lbl);
+    Field3D e_vel(1.0);
+    Options state{{comp_name,
+                   {
+                       {"type", this->parser.get_reaction_str()},
+                       {"diagnose", true},
+                   }},
+                  {"units", {{"eV", 1.0}, {"inv_meters_cubed", 1.0}, {"seconds", 1.0}}},
+                  {"species", {{"e", {{"AA", 1. / 1836}}}}}};
+
+    // Linear functions for various fields that are inputs to the reaction transforms
+
+    // Cycle through axes when setting fields, starting with x
+    // (species names are sorted, so the state generated should be order independent)
+    linfunc_axis axis = linfunc_axis::x;
+
+    // Densities and temperatures for electrons
+    state["species"]["e"]["density"] = FieldFactory::get()->create3D(
+        this->gen_lin_field_str(this->logn_min, this->logn_max, axis++), &state, mesh);
+    state["species"]["e"]["temperature"] = FieldFactory::get()->create3D(
+        this->gen_lin_field_str(this->logT_min, this->logT_max, axis++), &state, mesh);
+
+    // Densities and temperatures for all heavy reactants
+    for (const auto& reactant : this->heavy_reactants) {
+      state["species"][reactant]["density"] = FieldFactory::get()->create3D(
+          this->gen_lin_field_str(this->logn_min, this->logn_max, axis++), &state, mesh);
+
+      state["species"][reactant]["temperature"] = FieldFactory::get()->create3D(
+          this->gen_lin_field_str(this->logT_min, this->logT_max, axis++), &state, mesh);
+    }
+
+    // Masses, velocities for all heavy species
+    for (const auto& sp_name : this->heavy_species) {
+      state["species"][sp_name]["AA"] = 1.0;
+      state["species"][sp_name]["velocity"] = FieldFactory::get()->create3D(
+          this->gen_lin_field_str(this->logv_min, this->logv_max, axis++), &state, mesh);
+    }
+    return state;
+  }
+};
+
+// Dissociation of H2
+class H2DissTest : public DissociativeReactionTest<MolHDissociation<HIsotope::h>> {
+public:
+  H2DissTest()
+      : DissociativeReactionTest<MolHDissociation<HIsotope::h>>("H2Diss",
+                                                                "h2 + e -> 2h + e") {}
+};
+
+// Dissociative excitation of T2
+class T2pDissExcTest : public DissociativeReactionTest<MolHDissociativeExc<HIsotope::t>> {
+public:
+  T2pDissExcTest()
+      : DissociativeReactionTest<MolHDissociativeExc<HIsotope::t>>(
+            "T2pDissExc", "t2+ + e -> t + t+ + e") {}
+};
+
+// Dissociative ionisation of D2
+class D2DissIznTest : public DissociativeReactionTest<MolHDissociativeIzn<HIsotope::d>> {
+public:
+  D2DissIznTest()
+      : DissociativeReactionTest<MolHDissociativeIzn<HIsotope::d>>(
+            "D2DissIzn", "d2 + e -> d + d+ + 2e") {}
+};
+
+// Non-dissociative ionisation of H2
+class H2NonDissIznTest
+    : public DissociativeReactionTest<MolHNonDissociativeIzn<HIsotope::h>> {
+public:
+  H2NonDissIznTest()
+      : DissociativeReactionTest<MolHNonDissociativeIzn<HIsotope::h>>(
+            "H2NonDissIzn", "h2 + e -> h2+ + 2e") {}
+};
+
+// Dissociative recombination of T2+
+class T2pDissRecTest : public DissociativeReactionTest<MolHDissociativeRec<HIsotope::t>> {
+public:
+  T2pDissRecTest()
+      : DissociativeReactionTest<MolHDissociativeRec<HIsotope::t>>("T2pDissRec",
+                                                                   "t2+ + e -> 2t") {}
+};
+
+// CX involving molecules
+class D2DpCXTest : public CXReactionTest<CXReaction> {
+public:
+  D2DpCXTest() : CXReactionTest<CXReaction>("D2DpCX", "d2 + d+ -> d2+ + d") {}
+};
+
+// Elastic collision between H2 and H+
+// Use DissociativeReactionTest fixture for now
+class H2HpElasticCollisionTest
+    : public DissociativeReactionTest<MolHElasticCollision<HIsotope::h>> {
+public:
+  H2HpElasticCollisionTest()
+      : DissociativeReactionTest<MolHElasticCollision<HIsotope::h>>(
+            "H2HpElastic", "h2 + h+ -> h2 + h+") {}
+};
+
+} // namespace hermes
+
+#endif // TEST_MOLECULE_REACTIONS_H
diff --git a/tests/unit/test_reactions.cxx b/tests/unit/test_reactions.cxx
index ba886a230..56a472486 100644
--- a/tests/unit/test_reactions.cxx
+++ b/tests/unit/test_reactions.cxx
@@ -1,13 +1,14 @@
 #include "test_adas_reactions.hxx"
 #include "test_cx_reactions.hxx"
 #include "test_izn_rec_reactions.hxx"
+#include "test_molecule_reactions.hxx"
 
 namespace hermes {
 
-//======================== General reaction class tests =======================
+// ======================= General reaction class tests =======================
 /// @brief Test parsing of various input optionsReactionBase constructor should throw if
-/// the reaction type string is not
-TEST(ReactionTest, InputOptions) {
+/// the reaction type string is not valid
+TEST(ReactionConfigTest, InputOptions) {
   const std::string comp_name = "test";
 
   // Base input with two reaction strings
@@ -48,7 +49,23 @@ TEST(ReactionTest, InputOptions) {
   ASSERT_THROW(CXReaction(comp_name, invalid_input2), BoutException);
 }
 
-//======================== CX reaction class tests =======================
+/// @brief Values < 0 or > 1 for momentum/energy channels should throw
+TEST(ReactionConfigTest, InvalidChannelWeightOptions) {
+  Options base_options{
+      {"test", {{"type", "(h2 + d+ -> d+ + h2)"}, {"data_ids", "H.2_0.3T"}}},
+      {"units", {{"eV", 1.0}, {"inv_meters_cubed", 1.0}, {"seconds", 1.0}}}};
+
+  ReactionBase::reset_instance_counter();
+  Options options1 = base_options.copy();
+  auto invalid_mom_weights_reaction = InvalidMomWeightsReaction("test", options1);
+  ASSERT_THROW(invalid_mom_weights_reaction.transform(options1), BoutException);
+
+  ReactionBase::reset_instance_counter();
+  Options options2 = base_options.copy();
+  auto invalid_energy_weights_reaction = InvalidEnergyWeightsReaction("test", options2);
+  ASSERT_THROW(invalid_energy_weights_reaction.transform(options2), BoutException);
+}
+// ========================== CX reaction class tests =========================
 
 /// @brief CXReaction constructor should throw for strings that aren't valid CX
 /// reactions
@@ -97,7 +114,36 @@ TEST(CXReactionTest, OrderIndependentReactionStrs) {
   }
 }
 
-//====================== Reaction source regression tests =====================
+// =================== ElasticCollision reaction class tests ==================
+/// @brief Test validity of various reaction strings for elastic collisions
+TEST(ElasticCollisionTest, ReactionStrings) {
+  Options base_options{
+      {"test", {{"data_ids", "H.2_0.3T"}}},
+      {"units", {{"eV", 1.0}, {"inv_meters_cubed", 1.0}, {"seconds", 1.0}}}};
+
+  // Swapped order should work
+  std::string valid_reaction_str = "h2 + d+ -> d+ + h2 ";
+  ReactionBase::reset_instance_counter();
+  Options options = base_options.copy();
+  options["test"]["type"] = valid_reaction_str;
+  ASSERT_NO_THROW(ElasticCollision("test", options));
+
+  // Invalid elastic collision reaction strings
+  std::string too_few_species = "h -> h";
+  std::string too_many_species = "h + d+ + t -> h + d+ + t";
+  std::string reactants_products_mismatch = "h2 + d+ -> h2 +d";
+
+  // Test that constructor throws for each invalid reaction string
+  for (const auto& invalid_reaction_str :
+       {too_few_species, too_many_species, reactants_products_mismatch}) {
+    ReactionBase::reset_instance_counter();
+    Options options = base_options.copy();
+    options["test"]["type"] = invalid_reaction_str;
+    ASSERT_THROW(ElasticCollision("test", options), BoutException);
+  }
+}
+
+// ===================== Reaction source regression tests =====================
 
 // H isotopes ionization
 TEST_F(HIznTest, SourcesRegression) { sources_regression_test(); }
@@ -120,6 +166,21 @@ TEST_F(DTpCXTest, SourcesRegression) { sources_regression_test(); }
 // H/H+ CX with neutral momentum gain turned off
 TEST_F(HHpCXTest_noNeutralMomGain, SourcesRegression) { sources_regression_test(); }
 
+/* Various dissociation reactions involving neutral/ionic molecules of H
+ * isotopes(non-exhaustive)
+ */
+TEST_F(H2DissTest, SourcesRegression) { sources_regression_test(); }
+TEST_F(T2pDissExcTest, SourcesRegression) { sources_regression_test(); }
+TEST_F(D2DissIznTest, SourcesRegression) { sources_regression_test(); }
+TEST_F(H2NonDissIznTest, SourcesRegression) { sources_regression_test(); }
+TEST_F(T2pDissRecTest, SourcesRegression) { sources_regression_test(); }
+
+// CX involving molecules of H isotopes (non-exhaustive)
+TEST_F(D2DpCXTest, SourcesRegression) { sources_regression_test(); }
+
+// Elastic collision between H2 and H+
+TEST_F(H2HpElasticCollisionTest, SourcesRegression) { sources_regression_test(); }
+
 // He ionization
 TEST_F(HeIzn01Test, SourcesRegression) { sources_regression_test(); }
 
diff --git a/tests/unit/test_reactions.hxx b/tests/unit/test_reactions.hxx
index 3fb65361a..2927a66a4 100644
--- a/tests/unit/test_reactions.hxx
+++ b/tests/unit/test_reactions.hxx
@@ -22,6 +22,31 @@
 
 enum class linfunc_axis { x, y, z };
 
+/**
+ * @brief Prefix increment operator for linfunc_axis.
+ *
+ * @param current_val The current value of the enum.
+ * @return linfunc_axis& The incremented value of the enum.
+ */
+inline linfunc_axis& operator++(linfunc_axis& current) {
+  current = (current == linfunc_axis::z)
+                ? linfunc_axis::x
+                : static_cast<linfunc_axis>(static_cast<int>(current) + 1);
+  return current;
+}
+
+/**
+ * @brief Postfix increment operator for linfunc_axis.
+ *
+ * @param current The current value of the enum.
+ * @return linfunc_axis The old value of the enum.
+ */
+inline linfunc_axis operator++(linfunc_axis& current, int) {
+  linfunc_axis old = current;
+  ++current;
+  return old;
+}
+
 /// Global mesh
 namespace bout::globals {
 extern Mesh* mesh;
@@ -47,7 +72,7 @@ class ReactionTest : public FakeMeshFixture_tmpl<8, 8, 8> {
 
 protected:
   ReactionTest(std::string lbl, std::string reaction_str)
-      : lbl(lbl), parser(reaction_str){};
+      : lbl(lbl), parser(reaction_str) {};
 
   std::string lbl;
   ReactionParser parser;
@@ -216,6 +241,29 @@ protected:
   }
 };
 
+/// An invalid reaction ( > 1 total energy weight for the first reactant)
+struct InvalidEnergyWeightsReaction : public Reaction {
+  InvalidEnergyWeightsReaction(const std::string& name, Options& options)
+      : Reaction(name, options) {
+    std::string r1 = this->parser->get_reactant_by_position(1);
+    for (const std::string& product :
+         this->parser->get_species(species_filter::products)) {
+      this->set_energy_channel_weight(r1, product, 1.1);
+    }
+  }
+};
+/// An invalid reaction ( < 0 total momentum weight for the first reactant)
+struct InvalidMomWeightsReaction : public Reaction {
+  InvalidMomWeightsReaction(const std::string& name, Options& options)
+      : Reaction(name, options) {
+    std::string r1 = this->parser->get_reactant_by_position(1);
+    for (const std::string& product :
+         this->parser->get_species(species_filter::products)) {
+      this->set_momentum_channel_weight(r1, product, -0.1);
+    }
+  }
+};
+
 } // namespace hermes
 
 #endif
diff --git a/tests/unit/test_species_parser.cxx b/tests/unit/test_species_parser.cxx
index 1f69f8786..f8712b114 100644
--- a/tests/unit/test_species_parser.cxx
+++ b/tests/unit/test_species_parser.cxx
@@ -9,7 +9,7 @@ class SpeciesParserTest : public ::testing::Test {};
 /// @brief Check parsing of a neutral species with a prefix > 1
 TEST_F(SpeciesParserTest, ParseNeutral) {
   SpeciesParser parser("2he");
-  EXPECT_EQ(parser.get_element(), "he");
+  EXPECT_EQ(parser.get_base_species(), "he");
   EXPECT_EQ(parser.get_charge(), 0);
   // Species string should not include the prefix number and should be lower case
   EXPECT_EQ(parser.get_str(), "he");
@@ -18,17 +18,17 @@ TEST_F(SpeciesParserTest, ParseNeutral) {
 /// @brief Check parsing of various positively charged species.
 TEST_F(SpeciesParserTest, ParsePositivelyCharged) {
   SpeciesParser parser1("h+");
-  EXPECT_EQ(parser1.get_element(), "h");
+  EXPECT_EQ(parser1.get_base_species(), "h");
   EXPECT_EQ(parser1.get_charge(), 1);
   EXPECT_EQ(parser1.get_str(), "h+");
 
   SpeciesParser parser2("He+2");
-  EXPECT_EQ(parser2.get_element(), "he");
+  EXPECT_EQ(parser2.get_base_species(), "he");
   EXPECT_EQ(parser2.get_charge(), 2);
   EXPECT_EQ(parser2.get_str(), "he+2");
 
   SpeciesParser parser3("ne+8");
-  EXPECT_EQ(parser3.get_element(), "ne");
+  EXPECT_EQ(parser3.get_base_species(), "ne");
   EXPECT_EQ(parser3.get_charge(), 8);
   EXPECT_EQ(parser3.get_str(), "ne+8");
 }
@@ -37,17 +37,17 @@ TEST_F(SpeciesParserTest, ParsePositivelyCharged) {
 TEST_F(SpeciesParserTest, ParseNegativelyCharged) {
 
   SpeciesParser parser1("He-1");
-  EXPECT_EQ(parser1.get_element(), "he");
+  EXPECT_EQ(parser1.get_base_species(), "he");
   EXPECT_EQ(parser1.get_charge(), -1);
   EXPECT_EQ(parser1.get_str(), "he-");
 
   SpeciesParser parser2("li-2");
-  EXPECT_EQ(parser2.get_element(), "li");
+  EXPECT_EQ(parser2.get_base_species(), "li");
   EXPECT_EQ(parser2.get_charge(), -2);
   EXPECT_EQ(parser2.get_str(), "li-2");
 
   SpeciesParser parser3("ne-5");
-  EXPECT_EQ(parser3.get_element(), "ne");
+  EXPECT_EQ(parser3.get_base_species(), "ne");
   EXPECT_EQ(parser3.get_charge(), -5);
   EXPECT_EQ(parser3.get_str(), "ne-5");
 }
@@ -57,13 +57,13 @@ TEST_F(SpeciesParserTest, ParseElectron) {
   // "e" parsed as electron
   SpeciesParser parser1("e");
   EXPECT_EQ(parser1.get_charge(), -1);
-  EXPECT_EQ(parser1.get_element(), "e");
+  EXPECT_EQ(parser1.get_base_species(), "e");
   EXPECT_EQ(parser1.get_str(), "e");
 
   // "e-" parsed as electron
   SpeciesParser parser2("e-");
   EXPECT_EQ(parser2.get_charge(), -1);
-  EXPECT_EQ(parser2.get_element(), "e");
+  EXPECT_EQ(parser2.get_base_species(), "e");
   EXPECT_EQ(parser2.get_str(), "e");
 
   // Charges other than -1 throw for electrons
@@ -85,7 +85,7 @@ TEST_F(SpeciesParserTest, LongNames) {
   SpeciesParser parser5("t");
   ASSERT_EQ(parser5.long_name(), "tritium");
 
-  // Species without preset long names should return the element name
+  // Species without preset long names should return the base species name
   SpeciesParser parser6("li");
   ASSERT_EQ(parser6.long_name(), "li");
 }
@@ -94,7 +94,7 @@ TEST_F(SpeciesParserTest, LongNames) {
 TEST_F(SpeciesParserTest, IoniseNeutral) {
   SpeciesParser neutral("h");
   SpeciesParser ionised = neutral.ionised();
-  EXPECT_EQ(ionised.get_element(), "h");
+  EXPECT_EQ(ionised.get_base_species(), "h");
   EXPECT_EQ(ionised.get_charge(), 1);
   EXPECT_EQ(ionised.get_str(), "h+");
 }
@@ -103,7 +103,7 @@ TEST_F(SpeciesParserTest, IoniseNeutral) {
 TEST_F(SpeciesParserTest, IoniseSinglyIonised) {
   SpeciesParser singly("li+");
   SpeciesParser doubly = singly.ionised();
-  EXPECT_EQ(doubly.get_element(), "li");
+  EXPECT_EQ(doubly.get_base_species(), "li");
   EXPECT_EQ(doubly.get_charge(), 2);
   EXPECT_EQ(doubly.get_str(), "li+2");
 }
@@ -112,7 +112,7 @@ TEST_F(SpeciesParserTest, IoniseSinglyIonised) {
 TEST_F(SpeciesParserTest, IoniseHighlyIonised) {
   SpeciesParser highly("ne+8");
   SpeciesParser next = highly.ionised();
-  EXPECT_EQ(next.get_element(), "ne");
+  EXPECT_EQ(next.get_base_species(), "ne");
   EXPECT_EQ(next.get_charge(), 9);
   EXPECT_EQ(next.get_str(), "ne+9");
 }
@@ -125,12 +125,12 @@ TEST_F(SpeciesParserTest, MultipleIonisations) {
 
   SpeciesParser hep1 = he.ionised();
   EXPECT_EQ(hep1.get_charge(), 1);
-  EXPECT_EQ(hep1.get_element(), "he");
+  EXPECT_EQ(hep1.get_base_species(), "he");
   EXPECT_EQ(hep1.get_str(), "he+");
 
   SpeciesParser hep2 = hep1.ionised();
   EXPECT_EQ(hep2.get_charge(), 2);
-  EXPECT_EQ(hep2.get_element(), "he");
+  EXPECT_EQ(hep2.get_base_species(), "he");
   EXPECT_EQ(hep2.get_str(), "he+2");
 }
 
@@ -138,7 +138,7 @@ TEST_F(SpeciesParserTest, MultipleIonisations) {
 TEST_F(SpeciesParserTest, RecombineSinglyIonised) {
   SpeciesParser singly("h+");
   SpeciesParser neutral = singly.recombined();
-  EXPECT_EQ(neutral.get_element(), "h");
+  EXPECT_EQ(neutral.get_base_species(), "h");
   EXPECT_EQ(neutral.get_charge(), 0);
   EXPECT_EQ(neutral.get_str(), "h");
 }
@@ -147,7 +147,7 @@ TEST_F(SpeciesParserTest, RecombineSinglyIonised) {
 TEST_F(SpeciesParserTest, RecombineHighlyIonised) {
   SpeciesParser highly("ne+9");
   SpeciesParser next = highly.recombined();
-  EXPECT_EQ(next.get_element(), "ne");
+  EXPECT_EQ(next.get_base_species(), "ne");
   EXPECT_EQ(next.get_charge(), 8);
   EXPECT_EQ(next.get_str(), "ne+8");
 }
@@ -159,12 +159,12 @@ TEST_F(SpeciesParserTest, MultipleRecombinations) {
 
   SpeciesParser hep1 = hep2.recombined();
   EXPECT_EQ(hep1.get_charge(), 1);
-  EXPECT_EQ(hep1.get_element(), "he");
+  EXPECT_EQ(hep1.get_base_species(), "he");
   EXPECT_EQ(hep1.get_str(), "he+");
 
   SpeciesParser he0 = hep1.recombined();
   EXPECT_EQ(he0.get_charge(), 0);
-  EXPECT_EQ(he0.get_element(), "he");
+  EXPECT_EQ(he0.get_base_species(), "he");
   EXPECT_EQ(he0.get_str(), "he");
 }
 
@@ -187,4 +187,4 @@ TEST_F(SpeciesParserTest, InvalidStrings) {
   ASSERT_THROW(SpeciesParser("%&?!"), BoutException);
   ASSERT_THROW(SpeciesParser("he +"), BoutException);
   ASSERT_THROW(SpeciesParser("H -"), BoutException);
-}
\ No newline at end of file
+}