z.ab_patch

diff --git a/INSTALL.sh b/INSTALL.sh
new file mode 100644
index 000000000..0e158dff6
--- /dev/null
+++ b/INSTALL.sh
@@ -0,0 +1,20 @@
+CXX=g++
+CC=gcc
+cxx_flags="-std=c++17"
+prefix=${HOME}/opt/
+python_package=n
+f90_interface=n
+system_eigen=n
+extra_inc_dirs="/opt/include /usr/include/yaml-cpp"  # Include both directories
+system_blas_lapack=n
+boost_inc_dir=$(pwd)/ext/cliboost
+system_sundials=n
+system_yamlcpp=n
+
+scons build CXX=${CXX} CC=${CC} cxx_flags="${cxx_flags}" prefix=${prefix} \
+ python_package=${python_package} f90_interface=${f90_interface} \
+ system_eigen=${system_eigen} extra_inc_dirs="${extra_inc_dirs}" \
+ system_blas_lapack=${system_blas_lapack} boost_inc_dir=${boost_inc_dir} system_sundials=${system_sundials} \
+ system_yamlcpp=${system_yamlcpp} -j8
+
+scons install
diff --git a/include/cantera/kinetics/Arrhenius.h b/include/cantera/kinetics/Arrhenius.h
index 2791faca8..473cfb262 100644
--- a/include/cantera/kinetics/Arrhenius.h
+++ b/include/cantera/kinetics/Arrhenius.h
@@ -1,212 +1,215 @@
-/**
- * @file Arrhenius.h
- * Header for reaction rates that involve Arrhenius-type kinetics.
- */
-
-// This file is part of Cantera. See License.txt in the top-level directory or
-// at https://cantera.org/license.txt for license and copyright information.
-
-#ifndef CT_ARRHENIUS_H
-#define CT_ARRHENIUS_H
-
-#include "cantera/base/ct_defs.h"
-#include "cantera/base/Units.h"
-#include "cantera/kinetics/ReactionData.h"
-#include "ReactionRate.h"
-#include "MultiRate.h"
-
-namespace Cantera
-{
-
-class AnyValue;
-class AnyMap;
-
-//! Data container holding shared data specific to ArrheniusRate
-/**
- * The data container `ArrheniusData` holds precalculated data common to
- * all `ArrheniusRate` objects.
- */
-struct ArrheniusData : public ReactionData
-{
-    bool update(const ThermoPhase& phase, const Kinetics& kin) override;
-    using ReactionData::update;
-};
-
-
-//! Base class for Arrhenius-type Parameterizations
-/*!
- * This base class provides a minimally functional interface that allows for parameter
- * access from derived classes as well as classes that use Arrhenius-type expressions
- * internally, for example FalloffRate and PlogRate.
- * @ingroup arrheniusGroup
- */
-class ArrheniusBase : public ReactionRate
-{
-public:
-    //! Default constructor.
-    ArrheniusBase() {}
-
-    //! Constructor.
-    /*!
-     *  @param A  Pre-exponential factor. The unit system is (kmol, m, s); actual units
-     *      depend on the reaction order and the dimensionality (surface or bulk).
-     *  @param b  Temperature exponent (non-dimensional)
-     *  @param Ea  Activation energy in energy units [J/kmol]
-     */
-    ArrheniusBase(double A, double b, double Ea);
-
-    //! Constructor based on AnyValue content
-    ArrheniusBase(const AnyValue& rate, const UnitSystem& units,
-                  const UnitStack& rate_units);
-
-    explicit ArrheniusBase(const AnyMap& node, const UnitStack& rate_units={});
-
-    //! Perform object setup based on AnyValue node information
-    /*!
-     *  Used to set parameters from a child of the reaction node, which may have
-     *  different names for different rate parameterizations, such as falloff rates.
-     *
-     *  @param rate  Child of the reaction node containing Arrhenius rate parameters.
-     *      For example, the `rate-coefficient` node for a standard Arrhenius reaction.
-     *  @param units  Unit system
-     *  @param rate_units  Unit definitions specific to rate information
-     */
-    void setRateParameters(const AnyValue& rate,
-                           const UnitSystem& units,
-                           const UnitStack& rate_units);
-
-    //! Get Arrhenius parameters used to populate the `rate-coefficient` or
-    //! equivalent field
-    void getRateParameters(AnyMap& node) const;
-
-    void setParameters(const AnyMap& node, const UnitStack& rate_units) override;
-
-    void getParameters(AnyMap& node) const override;
-
-    //! Check rate expression
-    void check(const string& equation) override;
-
-    void validate(const string& equation, const Kinetics& kin) override;
-
-    //! Return the pre-exponential factor *A* (in m, kmol, s to powers depending
-    //! on the reaction order)
-    /*!
-     * Class specializations may provide alternate definitions that describe
-     * an effective pre-exponential factor that depends on the thermodynamic state.
-     */
-    virtual double preExponentialFactor() const {
-        return m_A;
-    }
-
-    //! Return the temperature exponent *b*
-    /*!
-     * Class specializations may provide alternate definitions that describe
-     * an effective temperature exponent that depends on the thermodynamic state.
-     */
-    virtual double temperatureExponent() const {
-        return m_b;
-    }
-
-    //! Return the activation energy *Ea* [J/kmol]
-    //! The value corresponds to the constant specified by input parameters;
-    /*!
-     * Class specializations may provide alternate definitions that describe
-     * an effective activation energy that depends on the thermodynamic state.
-     */
-    virtual double activationEnergy() const {
-        return m_Ea_R * GasConstant;
-    }
-
-    //! Return reaction order associated with the reaction rate
-    double order() const {
-        return m_order;
-    }
-
-    //! Set units of the reaction rate expression
-    void setRateUnits(const UnitStack& rate_units) override {
-        ReactionRate::setRateUnits(rate_units);
-        if (rate_units.size() > 1) {
-            m_order = 1 - rate_units.product().dimension("quantity");
-        } else {
-            m_order = NAN;
-        }
-    }
-
-    //! Get flag indicating whether negative A values are permitted
-    bool allowNegativePreExponentialFactor() const {
-        return m_negativeA_ok;
-    }
-
-    //! Set flag indicating whether negative A values are permitted
-    void setAllowNegativePreExponentialFactor(bool value) {
-        m_negativeA_ok = value;
-    }
-
-protected:
-    bool m_negativeA_ok = false; //!< Permissible negative A values
-    double m_A = NAN; //!< Pre-exponential factor
-    double m_b = NAN; //!< Temperature exponent
-    double m_Ea_R = 0.; //!< Activation energy (in temperature units)
-    double m_E4_R = 0.; //!< Optional 4th energy parameter (in temperature units)
-    double m_logA = NAN; //!< Logarithm of pre-exponential factor
-    double m_order = NAN; //!< Reaction order
-    string m_A_str = "A"; //!< The string for the pre-exponential factor
-    string m_b_str = "b"; //!< The string for temperature exponent
-    string m_Ea_str = "Ea"; //!< The string for activation energy
-    string m_E4_str = ""; //!< The string for an optional 4th parameter
-};
-
-//! Arrhenius reaction rate type depends only on temperature
-/*!
- * A reaction rate coefficient of the following form.
- *
- *   @f[
- *        k_f =  A T^b \exp (-Ea/RT)
- *   @f]
- *
- * @ingroup arrheniusGroup
- */
-class ArrheniusRate : public ArrheniusBase
-{
-public:
-    using ArrheniusBase::ArrheniusBase; // inherit constructors
-
-    unique_ptr<MultiRateBase> newMultiRate() const override {
-        return make_unique<MultiRate<ArrheniusRate, ArrheniusData>>();
-    }
-
-    const string type() const override {
-        return "Arrhenius";
-    }
-
-    //! Evaluate reaction rate
-    double evalRate(double logT, double recipT) const {
-        return m_A * std::exp(m_b * logT - m_Ea_R * recipT);
-    }
-
-    //! Evaluate natural logarithm of the rate constant.
-    double evalLog(double logT, double recipT) const {
-        return m_logA + m_b * logT - m_Ea_R * recipT;
-    }
-
-    //! Evaluate reaction rate
-    /*!
-     *  @param shared_data  data shared by all reactions of a given type
-     */
-    double evalFromStruct(const ArrheniusData& shared_data) const {
-        return m_A * std::exp(m_b * shared_data.logT - m_Ea_R * shared_data.recipT);
-    }
-
-    //! Evaluate derivative of reaction rate with respect to temperature
-    //! divided by reaction rate
-    /*!
-     *  @param shared_data  data shared by all reactions of a given type
-     */
-    double ddTScaledFromStruct(const ArrheniusData& shared_data) const {
-        return (m_Ea_R * shared_data.recipT + m_b) * shared_data.recipT;
-    }
-};
-
-}
-
-#endif
+/**
+ * @file Arrhenius.h
+ * Header for reaction rates that involve Arrhenius-type kinetics.
+ */
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#ifndef CT_ARRHENIUS_H
+#define CT_ARRHENIUS_H
+
+#include "cantera/base/ct_defs.h"
+#include "cantera/base/Units.h"
+#include "cantera/kinetics/ReactionData.h"
+#include "ReactionRate.h"
+#include "MultiRate.h"
+
+namespace Cantera
+{
+
+class AnyValue;
+class AnyMap;
+
+//! Data container holding shared data specific to ArrheniusRate
+/**
+ * The data container `ArrheniusData` holds precalculated data common to
+ * all `ArrheniusRate` objects.
+ */
+struct ArrheniusData : public ReactionData
+{
+    bool update(const ThermoPhase& phase, const Kinetics& kin) override;
+    using ReactionData::update;
+};
+
+
+//! Base class for Arrhenius-type Parameterizations
+/*!
+ * This base class provides a minimally functional interface that allows for parameter
+ * access from derived classes as well as classes that use Arrhenius-type expressions
+ * internally, for example FalloffRate and PlogRate.
+ * @ingroup arrheniusGroup
+ */
+class ArrheniusBase : public ReactionRate
+{
+public:
+    //! Default constructor.
+    ArrheniusBase() {}
+
+    //! Constructor.
+    /*!
+     *  @param A  Pre-exponential factor. The unit system is (kmol, m, s); actual units
+     *      depend on the reaction order and the dimensionality (surface or bulk).
+     *  @param T0 Reference temperature (K)
+     *  @param b  Temperature exponent (non-dimensional)
+     *  @param Ea  Activation energy in energy units [J/kmol]
+     */
+    ArrheniusBase(double A, double b,double T0,double Ea);
+
+    //! Constructor based on AnyValue content
+    ArrheniusBase(const AnyValue& rate, const UnitSystem& units,
+                  const UnitStack& rate_units);
+
+    explicit ArrheniusBase(const AnyMap& node, const UnitStack& rate_units={});
+
+    //! Perform object setup based on AnyValue node information
+    /*!
+     *  Used to set parameters from a child of the reaction node, which may have
+     *  different names for different rate parameterizations, such as falloff rates.
+     *
+     *  @param rate  Child of the reaction node containing Arrhenius rate parameters.
+     *      For example, the `rate-coefficient` node for a standard Arrhenius reaction.
+     *  @param units  Unit system
+     *  @param rate_units  Unit definitions specific to rate information
+     */
+    void setRateParameters(const AnyValue& rate,
+                           const UnitSystem& units,
+                           const UnitStack& rate_units);
+
+    //! Get Arrhenius parameters used to populate the `rate-coefficient` or
+    //! equivalent field
+    void getRateParameters(AnyMap& node) const;
+
+    void setParameters(const AnyMap& node, const UnitStack& rate_units) override;
+
+    void getParameters(AnyMap& node) const override;
+
+    //! Check rate expression
+    void check(const string& equation) override;
+
+    void validate(const string& equation, const Kinetics& kin) override;
+
+    //! Return the pre-exponential factor *A* (in m, kmol, s to powers depending
+    //! on the reaction order)
+    /*!
+     * Class specializations may provide alternate definitions that describe
+     * an effective pre-exponential factor that depends on the thermodynamic state.
+     */
+    virtual double preExponentialFactor() const {
+        return m_A;
+    }
+
+    //! Return the temperature exponent *b*
+    /*!
+     * Class specializations may provide alternate definitions that describe
+     * an effective temperature exponent that depends on the thermodynamic state.
+     */
+    virtual double temperatureExponent() const {
+        return m_b;
+    }
+
+    //! Return the activation energy *Ea* [J/kmol]
+    //! The value corresponds to the constant specified by input parameters;
+    /*!
+     * Class specializations may provide alternate definitions that describe
+     * an effective activation energy that depends on the thermodynamic state.
+     */
+    virtual double activationEnergy() const {
+        return m_Ea_R * GasConstant;
+    }
+
+    //! Return reaction order associated with the reaction rate
+    double order() const {
+        return m_order;
+    }
+
+    //! Set units of the reaction rate expression
+    void setRateUnits(const UnitStack& rate_units) override {
+        ReactionRate::setRateUnits(rate_units);
+        if (rate_units.size() > 1) {
+            m_order = 1 - rate_units.product().dimension("quantity");
+        } else {
+            m_order = NAN;
+        }
+    }
+
+    //! Get flag indicating whether negative A values are permitted
+    bool allowNegativePreExponentialFactor() const {
+        return m_negativeA_ok;
+    }
+
+    //! Set flag indicating whether negative A values are permitted
+    void setAllowNegativePreExponentialFactor(bool value) {
+        m_negativeA_ok = value;
+    }
+
+protected:
+    bool m_negativeA_ok = false; //!< Permissible negative A values
+    double m_A = NAN; //!< Pre-exponential factor
+    double m_b = NAN; //!< Temperature exponent
+    double m_T0 = 1.0; //!< Reference temperature
+    double m_Ea_R = 0.; //!< Activation energy (in temperature units)
+    double m_E4_R = 0.; //!< Optional 4th energy parameter (in temperature units)
+    double m_logA = NAN; //!< Logarithm of pre-exponential factor
+    double m_order = NAN; //!< Reaction order
+    string m_A_str = "A"; //!< The string for the pre-exponential factor
+    string m_b_str = "b"; //!< The string for temperature exponent
+    string m_T0_str = "T0"; //!< The string for reference temperature
+    string m_Ea_str = "Ea"; //!< The string for activation energy
+    string m_E4_str = ""; //!< The string for an optional 4th parameter
+};
+
+//! Arrhenius reaction rate type depends only on temperature
+/*!
+ * A reaction rate coefficient of the following form.
+ *
+ *   @f[
+ *        k_f =  A (T/T0)^b \exp (-Ea/RT)
+ *   @f]
+ *
+ * @ingroup arrheniusGroup
+ */
+class ArrheniusRate : public ArrheniusBase
+{
+public:
+    using ArrheniusBase::ArrheniusBase; // inherit constructors
+
+    unique_ptr<MultiRateBase> newMultiRate() const override {
+        return make_unique<MultiRate<ArrheniusRate, ArrheniusData>>();
+    }
+
+    const string type() const override {
+        return "T-Arrhenius";
+    }
+
+    //! Evaluate reaction rate
+    double evalRate(double logT, double recipT) const {
+        return m_A * std::exp((m_b * (logT - std::log(m_T0))) - (m_Ea_R * recipT));
+    }
+
+    //! Evaluate natural logarithm of the rate constant.
+    double evalLog(double logT, double recipT) const {
+        return m_logA + (m_b * (logT - std::log(m_T0))) - (m_Ea_R * recipT);
+    }
+
+    //! Evaluate reaction rate
+    /*!
+     *  @param shared_data  data shared by all reactions of a given type
+     */
+    double evalFromStruct(const ArrheniusData& shared_data) const {
+        return m_A * std::exp((m_b * (shared_data.logT - std::log(m_T0))) - (m_Ea_R * shared_data.recipT));
+    }
+
+    //! Evaluate derivative of reaction rate with respect to temperature
+    //! divided by reaction rate
+    /*!
+     *  @param shared_data  data shared by all reactions of a given type
+     */
+    double ddTScaledFromStruct(const ArrheniusData& shared_data) const {
+        return (m_Ea_R * shared_data.recipT + m_b) * shared_data.recipT;
+    }
+};
+
+}
+
+#endif
diff --git a/include/cantera/kinetics/BlowersMaselRate.h b/include/cantera/kinetics/BlowersMaselRate.h
index d10db0b72..cbef45a7a 100644
--- a/include/cantera/kinetics/BlowersMaselRate.h
+++ b/include/cantera/kinetics/BlowersMaselRate.h
@@ -81,7 +81,7 @@ public:
      *  @param w  Average bond dissociation energy of the bond being formed and
      *      broken in the reaction, in energy units [J/kmol]
      */
-    BlowersMaselRate(double A, double b, double Ea0, double w);
+    BlowersMaselRate(double A, double b, double T0, double Ea0, double w);
 
     explicit BlowersMaselRate(const AnyMap& node,
                               const UnitStack& rate_units={});
diff --git a/include/cantera/kinetics/Photolysis.h b/include/cantera/kinetics/Photolysis.h
index ee9b4c0cf..358e9f5b3 100644
--- a/include/cantera/kinetics/Photolysis.h
+++ b/include/cantera/kinetics/Photolysis.h
@@ -1,4 +1,6 @@
-//! @file Photolysis.h
+/** @file Photolysis.h
+ * Header for reaction rates that involve Photochemical reactions
+ */
 
 #ifndef CT_PHOTOLYSIS_H
 #define CT_PHOTOLYSIS_H
@@ -69,12 +71,16 @@ class PhotolysisBase : public ReactionRate {
   //! @param branch_map Map of branch names to branch indices
   void setRateParameters(const AnyValue& rate, map<string, int> const& branch_map);
 
+  //! Get the parameters corresponding to node rate-constant
   void getParameters(AnyMap& node) const override;
 
+  //! Get the parameters for a given node with flow style output
   void getRateParameters(AnyMap& node) const;
 
+  //! Checks for temperature range, and wavelength data
   void check(string const& equation) override;
 
+  //! Checks for valid species, stoichiometric balance, and consistency with photolysis branches
   void validate(const string& equation, const Kinetics& kin) override;
 
   vector<double> getCrossSection(double temp, double wavelength) const;
@@ -122,14 +128,21 @@ class PhotolysisRate : public PhotolysisBase {
     return make_unique<MultiRate<PhotolysisRate, PhotolysisData>>();
   }
 
+  //! reaction string type for photolysis reactions
   const string type() const override {
     return "Photolysis";
   }
-
+  
+  //! net stoichiometric coefficients of photolysis products
   Composition const& photoProducts() const override {
     return m_net_products;
   }
 
+/**
+ * @brief Calculates the photolysis reaction rate and updates stoichiometric concentration of products
+ * 
+ * @return total photolysis rate from all the branches
+ */
   double evalFromStruct(PhotolysisData const& data);
 
  protected:
diff --git a/include/cantera/kinetics/TwoTempPlasmaRate.h b/include/cantera/kinetics/TwoTempPlasmaRate.h
index 7f17e9de7..f09f11b55 100644
--- a/include/cantera/kinetics/TwoTempPlasmaRate.h
+++ b/include/cantera/kinetics/TwoTempPlasmaRate.h
@@ -67,10 +67,11 @@ public:
      *  @param A  Pre-exponential factor. The unit system is (kmol, m, s); actual units
      *      depend on the reaction order and the dimensionality (surface or bulk).
      *  @param b  Temperature exponent (non-dimensional)
+     *  @param T0 Reference temperature (K)
      *  @param Ea  Activation energy in energy units [J/kmol]
      *  @param EE  Activation electron energy in energy units [J/kmol]
      */
-    TwoTempPlasmaRate(double A, double b, double Ea=0.0, double EE=0.0);
+    TwoTempPlasmaRate(double A, double b, double T0, double Ea=0.0, double EE=0.0);
 
     TwoTempPlasmaRate(const AnyMap& node, const UnitStack& rate_units={});
 
diff --git a/src/kinetics/Arrhenius.cpp b/src/kinetics/Arrhenius.cpp
index ecd89b203..42f82d618 100644
--- a/src/kinetics/Arrhenius.cpp
+++ b/src/kinetics/Arrhenius.cpp
@@ -1,156 +1,160 @@
-//! @file Arrhenius.cpp
-
-// This file is part of Cantera. See License.txt in the top-level directory or
-// at https://cantera.org/license.txt for license and copyright information.
-
-#include "cantera/kinetics/Arrhenius.h"
-#include "cantera/thermo/ThermoPhase.h"
-
-namespace Cantera
-{
-
-ArrheniusBase::ArrheniusBase(double A, double b, double Ea)
-    : m_A(A)
-    , m_b(b)
-    , m_Ea_R(Ea / GasConstant)
-{
-    if (m_A > 0.0) {
-        m_logA = std::log(m_A);
-    }
-    m_valid = true;
-}
-
-ArrheniusBase::ArrheniusBase(const AnyValue& rate, const UnitSystem& units,
-                             const UnitStack& rate_units)
-{
-    setRateUnits(rate_units);
-    setRateParameters(rate, units, rate_units);
-}
-
-ArrheniusBase::ArrheniusBase(const AnyMap& node, const UnitStack& rate_units)
-{
-    setParameters(node, rate_units);
-}
-
-void ArrheniusBase::setRateParameters(
-    const AnyValue& rate, const UnitSystem& units, const UnitStack& rate_units)
-{
-    m_Ea_R = 0.; // assume zero if not provided
-    m_E4_R = 0.; // assume zero if not provided
-    if (rate.empty()) {
-        m_A = NAN;
-        m_b = NAN;
-        m_logA = NAN;
-        setRateUnits(Units(0.));
-        return;
-    }
-
-    if (rate.is<AnyMap>()) {
-
-        auto& rate_map = rate.as<AnyMap>();
-        m_A = units.convertRateCoeff(rate_map[m_A_str], conversionUnits());
-        m_b = rate_map[m_b_str].asDouble();
-        if (rate_map.hasKey(m_Ea_str)) {
-            m_Ea_R = units.convertActivationEnergy(rate_map[m_Ea_str], "K");
-        }
-        if (rate_map.hasKey(m_E4_str)) {
-            m_E4_R = units.convertActivationEnergy(rate_map[m_E4_str], "K");
-        }
-    } else {
-        auto& rate_vec = rate.asVector<AnyValue>(2, 4);
-        m_A = units.convertRateCoeff(rate_vec[0], conversionUnits());
-        m_b = rate_vec[1].asDouble();
-        if (rate_vec.size() > 2) {
-            m_Ea_R = units.convertActivationEnergy(rate_vec[2], "K");
-        }
-        if (rate_vec.size() > 3) {
-            m_E4_R = units.convertActivationEnergy(rate_vec[3], "K");
-        }
-    }
-    if (m_A > 0.0) {
-        m_logA = std::log(m_A);
-    }
-    m_valid = true;
-}
-
-void ArrheniusBase::getRateParameters(AnyMap& node) const
-{
-    if (!valid()) {
-        // Return empty/unmodified AnyMap
-        return;
-    }
-
-    if (conversionUnits().factor() != 0.0) {
-        node[m_A_str].setQuantity(m_A, conversionUnits());
-    } else {
-        node[m_A_str] = m_A;
-        // This can't be converted to a different unit system because the dimensions of
-        // the rate constant were not set. Can occur if the reaction was created outside
-        // the context of a Kinetics object and never added to a Kinetics object.
-        node["__unconvertible__"] = true;
-    }
-    node[m_b_str] = m_b;
-    node[m_Ea_str].setQuantity(m_Ea_R, "K", true);
-    if (m_E4_str != "") {
-        node[m_E4_str].setQuantity(m_E4_R, "K", true);
-    }
-    node.setFlowStyle();
-}
-
-void ArrheniusBase::setParameters(const AnyMap& node, const UnitStack& rate_units)
-{
-    ReactionRate::setParameters(node, rate_units);
-    m_negativeA_ok = node.getBool("negative-A", false);
-    if (!node.hasKey("rate-constant")) {
-        setRateParameters(AnyValue(), node.units(), rate_units);
-        return;
-    }
-    setRateParameters(node["rate-constant"], node.units(), rate_units);
-}
-
-void ArrheniusBase::getParameters(AnyMap& node) const {
-    if (m_negativeA_ok) {
-        node["negative-A"] = true;
-    }
-    AnyMap rateNode;
-    getRateParameters(rateNode);
-    if (!rateNode.empty()) {
-        // RateType object is configured
-        node["rate-constant"] = std::move(rateNode);
-    }
-}
-
-void ArrheniusBase::check(const string& equation)
-{
-    if (!m_negativeA_ok && m_A < 0) {
-        if (equation == "") {
-            throw CanteraError("ArrheniusBase::check",
-                "Detected negative pre-exponential factor (A={}).\n"
-                "Enable 'allowNegativePreExponentialFactor' to suppress "
-                "this message.", m_A);
-        }
-        throw InputFileError("ArrheniusBase::check", m_input,
-            "Undeclared negative pre-exponential factor found in reaction '{}'",
-            equation);
-    }
-}
-
-void ArrheniusBase::validate(const string& equation, const Kinetics& kin)
-{
-    if (!valid()) {
-        throw InputFileError("ArrheniusBase::validate", m_input,
-            "Rate object for reaction '{}' is not configured.", equation);
-    }
-}
-
-bool ArrheniusData::update(const ThermoPhase& phase, const Kinetics& kin)
-{
-    double T = phase.temperature();
-    if (T == temperature) {
-        return false;
-    }
-    update(T);
-    return true;
-}
-
-}
+//! @file Arrhenius.cpp
+
+// This file is part of Cantera. See License.txt in the top-level directory or
+// at https://cantera.org/license.txt for license and copyright information.
+
+#include "cantera/kinetics/Arrhenius.h"
+#include "cantera/thermo/ThermoPhase.h"
+
+namespace Cantera
+{
+
+ArrheniusBase::ArrheniusBase(double A, double b, double T0, double Ea)
+    : m_A(A)
+    , m_b(b)
+    , m_T0(T0)
+    , m_Ea_R(Ea / GasConstant)
+{
+    if (m_A > 0.0) {
+        m_logA = std::log(m_A);
+    }
+    m_valid = true;
+}
+
+ArrheniusBase::ArrheniusBase(const AnyValue& rate, const UnitSystem& units,
+                             const UnitStack& rate_units)
+{
+    setRateUnits(rate_units);
+    setRateParameters(rate, units, rate_units);
+}
+
+ArrheniusBase::ArrheniusBase(const AnyMap& node, const UnitStack& rate_units)
+{
+    setParameters(node, rate_units);
+}
+
+void ArrheniusBase::setRateParameters(
+    const AnyValue& rate, const UnitSystem& units, const UnitStack& rate_units)
+{
+    m_Ea_R = 0.; // assume zero if not provided
+    m_E4_R = 0.; // assume zero if not provided
+    if (rate.empty()) {
+        m_A = NAN;
+        m_b = NAN;
+	m_T0 = 1.0;
+        m_logA = NAN;
+        setRateUnits(Units(0.));
+        return;
+    }
+
+    if (rate.is<AnyMap>()) {
+
+        auto& rate_map = rate.as<AnyMap>();
+        m_A = units.convertRateCoeff(rate_map[m_A_str], conversionUnits());
+        m_b = rate_map[m_b_str].asDouble();
+	m_T0 = rate_map[m_T0_str].asDouble();
+        if (rate_map.hasKey(m_Ea_str)) {
+            m_Ea_R = units.convertActivationEnergy(rate_map[m_Ea_str], "K");
+        }
+        if (rate_map.hasKey(m_E4_str)) {
+            m_E4_R = units.convertActivationEnergy(rate_map[m_E4_str], "K");
+        }
+    } else {
+        auto& rate_vec = rate.asVector<AnyValue>(2, 4);
+        m_A = units.convertRateCoeff(rate_vec[0], conversionUnits());
+        m_b = rate_vec[1].asDouble();
+        if (rate_vec.size() > 2) {
+            m_Ea_R = units.convertActivationEnergy(rate_vec[2], "K");
+        }
+        if (rate_vec.size() > 3) {
+            m_E4_R = units.convertActivationEnergy(rate_vec[3], "K");
+        }
+    }
+    if (m_A > 0.0) {
+        m_logA = std::log(m_A);
+    }
+    m_valid = true;
+}
+
+void ArrheniusBase::getRateParameters(AnyMap& node) const
+{
+    if (!valid()) {
+        // Return empty/unmodified AnyMap
+        return;
+    }
+
+    if (conversionUnits().factor() != 0.0) {
+        node[m_A_str].setQuantity(m_A, conversionUnits());
+    } else {
+        node[m_A_str] = m_A;
+        // This can't be converted to a different unit system because the dimensions of
+        // the rate constant were not set. Can occur if the reaction was created outside
+        // the context of a Kinetics object and never added to a Kinetics object.
+        node["__unconvertible__"] = true;
+    }
+    node[m_b_str] = m_b;
+    node[m_T0_str] = m_T0;
+    node[m_Ea_str].setQuantity(m_Ea_R, "K", true);
+    if (m_E4_str != "") {
+        node[m_E4_str].setQuantity(m_E4_R, "K", true);
+    }
+    node.setFlowStyle();
+}
+
+void ArrheniusBase::setParameters(const AnyMap& node, const UnitStack& rate_units)
+{
+    ReactionRate::setParameters(node, rate_units);
+    m_negativeA_ok = node.getBool("negative-A", false);
+    if (!node.hasKey("rate-constant")) {
+        setRateParameters(AnyValue(), node.units(), rate_units);
+        return;
+    }
+    setRateParameters(node["rate-constant"], node.units(), rate_units);
+}
+
+void ArrheniusBase::getParameters(AnyMap& node) const {
+    if (m_negativeA_ok) {
+        node["negative-A"] = true;
+    }
+    AnyMap rateNode;
+    getRateParameters(rateNode);
+    if (!rateNode.empty()) {
+        // RateType object is configured
+        node["rate-constant"] = std::move(rateNode);
+    }
+}
+
+void ArrheniusBase::check(const string& equation)
+{
+    if (!m_negativeA_ok && m_A < 0) {
+        if (equation == "") {
+            throw CanteraError("ArrheniusBase::check",
+                "Detected negative pre-exponential factor (A={}).\n"
+                "Enable 'allowNegativePreExponentialFactor' to suppress "
+                "this message.", m_A);
+        }
+        throw InputFileError("ArrheniusBase::check", m_input,
+            "Undeclared negative pre-exponential factor found in reaction '{}'",
+            equation);
+    }
+}
+
+void ArrheniusBase::validate(const string& equation, const Kinetics& kin)
+{
+    if (!valid()) {
+        throw InputFileError("ArrheniusBase::validate", m_input,
+            "Rate object for reaction '{}' is not configured.", equation);
+    }
+}
+
+bool ArrheniusData::update(const ThermoPhase& phase, const Kinetics& kin)
+{
+    double T = phase.temperature();
+    if (T == temperature) {
+        return false;
+    }
+    update(T);
+    return true;
+}
+
+}
diff --git a/src/kinetics/BlowersMaselRate.cpp b/src/kinetics/BlowersMaselRate.cpp
index c51d8d718..528d17a13 100644
--- a/src/kinetics/BlowersMaselRate.cpp
+++ b/src/kinetics/BlowersMaselRate.cpp
@@ -42,8 +42,8 @@ BlowersMaselRate::BlowersMaselRate()
     m_E4_str = "w";
 }
 
-BlowersMaselRate::BlowersMaselRate(double A, double b, double Ea0, double w)
-    : ArrheniusBase(A, b, Ea0)
+BlowersMaselRate::BlowersMaselRate(double A, double b, double T0, double Ea0, double w)
+    : ArrheniusBase(A, b, T0, Ea0)
 {
     m_Ea_str = "Ea0";
     m_E4_str = "w";
diff --git a/src/kinetics/Photolysis.cpp b/src/kinetics/Photolysis.cpp
index 67d2ccb87..a87a8802d 100644
--- a/src/kinetics/Photolysis.cpp
+++ b/src/kinetics/Photolysis.cpp
@@ -44,11 +44,13 @@ bool PhotolysisData::check() const
                            "Wavelength grid must have at least two points.");
     }
 
+    // Check that wavelength grid values are positive
     if (wavelength[0] <= 0.0) {
         throw CanteraError("PhotolysisData::update",
                            "Wavelength grid must be positive.");
     }
 
+    // Check that wavelength grid values are monotonic and increasing
     for (size_t i = 1; i < wavelength.size(); i++) {
         if (wavelength[i] <= wavelength[i-1]) {
             throw CanteraError("PhotolysisData::update",
@@ -61,12 +63,14 @@ bool PhotolysisData::check() const
         throw CanteraError("PhotolysisData::update",
                            "Actinic flux is empty.");
     }
-
+    
+    // Check that actinic flux grid should have same size as wavelength grid
     if (actinicFlux.size() != wavelength.size()) {
         throw CanteraError("PhotolysisData::update",
                            "Actinic flux must have the same size as the wavelength grid.");
     }
 
+    // Check that actinic flux values are positive
     for (size_t i = 0; i < actinicFlux.size(); i++) {
         if (actinicFlux[i] < 0.0) {
             throw CanteraError("PhotolysisData::update",
@@ -87,6 +91,7 @@ PhotolysisBase::PhotolysisBase(
   m_ntemp = temp.size();
   m_nwave = wavelength.size();
 
+ // Grid for temperature and wavelength
   m_temp_wave_grid.resize(m_ntemp + m_nwave);
   for (size_t i = 0; i < m_ntemp; i++) {
     m_temp_wave_grid[i] = temp[i];
@@ -100,6 +105,7 @@ PhotolysisBase::PhotolysisBase(
     m_branch.push_back(parseCompString(branch));
   }
 
+  // Check if cross-section data size 
   if (m_ntemp * m_nwave * branches.size() != m_crossSection.size()) {
     throw CanteraError("PhotolysisBase::PhotolysisBase",
                        "Cross-section data size does not match the temperature, "
@@ -182,16 +188,19 @@ void PhotolysisBase::setParameters(AnyMap const& node, UnitStack const& rate_uni
   } else if (rtmp.products != rtmp.reactants) { // this is not photoabsorption
     m_branch.push_back(rtmp.products);
   }
-
+   
   if (node.hasKey("cross-section")) {
     for (auto const& data: node["cross-section"].asVector<AnyMap>()) {
       auto format = data["format"].asString();
       auto temp = data["temperature-range"].asVector<double>(2, 2);
+      
+      //Check temperature range to be monotonically increasing
       if (temp[0] >= temp[1]) {
         throw CanteraError("PhotolysisBase::setParameters",
                            "Temperature range must be strictly increasing.");
       }
-
+      
+      //Check for gaps in temperature range
       if (temperature.empty()) {
         temperature = temp;
       } else {
@@ -209,9 +218,12 @@ void PhotolysisBase::setParameters(AnyMap const& node, UnitStack const& rate_uni
           result.first.push_back(entry[0]);
           result.second.push_back(entry[1]);
         }
+      
+      //Read file names for VULCAN photochemistry database
       } else if (format == "VULCAN") {
         auto files = data["filenames"].asVector<string>();
         result = load_xsection_vulcan(files, m_branch);
+      //Read file name for KINETICS photochemistry database
       } else if (format == "KINETICS7") {
         auto files = data["filenames"].asVector<string>();
         result = load_xsection_kinetics7(files, m_branch);
@@ -267,11 +279,13 @@ void PhotolysisBase::setParameters(AnyMap const& node, UnitStack const& rate_uni
   m_valid = true;
 }
 
+ //Set the rate parameters to flow style ?? (TBD)
 void PhotolysisBase::getRateParameters(AnyMap& node) const
 {
   node.setFlowStyle();
 }
 
+ //Get rate parameters for photolysis reaction
 void PhotolysisBase::getParameters(AnyMap& node) const
 {
   AnyMap rateNode;
@@ -282,6 +296,7 @@ void PhotolysisBase::getParameters(AnyMap& node) const
   }
 }
 
+//Check temperature range, and wavelength data for photolysis reactions
 void PhotolysisBase::check(string const& equation)
 {
   if (m_ntemp < 2) {
@@ -296,6 +311,7 @@ void PhotolysisBase::check(string const& equation)
   }
 }
 
+//Check rate coefficient 
 void PhotolysisBase::validate(string const& equation, Kinetics const& kin)
 {
   if (!valid()) {
@@ -330,6 +346,7 @@ void PhotolysisBase::validate(string const& equation, Kinetics const& kin)
     }
   }
 
+  // Check for consistency in species in reaction string, and photolysis branches
   if (species_from_equation != species_from_branches) {
     throw InputFileError("PhotolysisBase::validate", m_input,
                        "Reaction '{}' has different products than the photolysis branches.", equation);
@@ -353,6 +370,7 @@ vector<double> PhotolysisBase::getCrossSection(double temp, double wavelength) c
   return cross;
 }
 
+// Evaluate the photolysis rate, and update the stoichiometric coefficient for different branches
 double PhotolysisRate::evalFromStruct(PhotolysisData const& data) {
     double wmin = m_temp_wave_grid[m_ntemp];
     double wmax = m_temp_wave_grid.back();
@@ -382,10 +400,9 @@ double PhotolysisRate::evalFromStruct(PhotolysisData const& data) {
 
     double coord[2] = {data.temperature, data.wavelength[0]};
     size_t len[2] = {m_ntemp, m_nwave};
+ 
 
-    // debug
-    //std::cout << "coord = " << coord[0] << " " << coord[1] << std::endl;
-
+    // N-space interpolation to determine photolysis cross section
     interpn(cross1, coord, m_crossSection.data(), m_temp_wave_grid.data(),
         len, 2, m_branch.size());
 
@@ -417,6 +434,7 @@ double PhotolysisRate::evalFromStruct(PhotolysisData const& data) {
 
       // photodissociation only
       for (size_t n = 1; n < m_branch.size(); n++) {
+        //Photochemical rate constant
         double rate = 0.5 * (data.wavelength[i+1] - data.wavelength[i])
           * (cross1[n] * data.actinicFlux[i] + cross2[n] * data.actinicFlux[i+1]);
 
diff --git a/src/kinetics/TwoTempPlasmaRate.cpp b/src/kinetics/TwoTempPlasmaRate.cpp
index ce3f46949..1a930ada5 100644
--- a/src/kinetics/TwoTempPlasmaRate.cpp
+++ b/src/kinetics/TwoTempPlasmaRate.cpp
@@ -51,8 +51,8 @@ TwoTempPlasmaRate::TwoTempPlasmaRate()
     m_E4_str = "Ea-electron";
 }
 
-TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double Ea, double EE)
-    : ArrheniusBase(A, b, Ea)
+TwoTempPlasmaRate::TwoTempPlasmaRate(double A, double b, double T0, double Ea, double EE)
+    : ArrheniusBase(A, b, T0, Ea)
 {
     m_Ea_str = "Ea-gas";
     m_E4_str = "Ea-electron";