WIP - Update electrochemical reaction type input for BV and Marcus kinetics

include/cantera/kinetics/InterfaceRate.h

@@ -57,6 +57,7 @@ struct InterfaceData : public BlowersMaselData
         electricPotentials.resize(nPhases, 0.);
         standardChemPotentials.resize(nSpecies, 0.);
         standardConcentrations.resize(nSpecies, 0.);
+        chemPotentials.resize(nSpecies, 0.);
         ready = true;
     }
 
@@ -67,6 +68,7 @@ struct InterfaceData : public BlowersMaselData
     vector_fp electricPotentials; //!< electric potentials of phases
     vector_fp standardChemPotentials; //!< standard state chemical potentials
     vector_fp standardConcentrations; //!< standard state concentrations
+    vector_fp chemPotentials;
 };
 
 
@@ -123,6 +125,11 @@ class InterfaceRateBase
         return m_exchangeCurrentDensityFormulation;
     }
 
+    // String storing the input electrochemical kinetics form used if non-standard
+    std::string eChemForm() {
+        return m_eChemForm;
+    }
+
     //! Build rate-specific parameters based on Reaction and Kinetics context
     //! @param rxn  Reaction associated with rate parameterization
     //! @param kin  Kinetics object associated with rate parameterization
@@ -161,23 +168,34 @@ class InterfaceRateBase
         // Calculate reaction rate correction. Only modify those with a non-zero
         // activation energy.
         double correction = 1.;
+        double beta_tmp = m_beta;
         if (m_deltaPotential_RT != 0.) {
             // Comments preserved from previous implementation:
             // Below we decrease the activation energy below zero.
             // NOTE, there is some discussion about this point. Should we decrease the
             // activation energy below zero? I don't think this has been decided in any
             // definitive way. The treatment below is numerically more stable, however.
-            correction = exp(-m_beta * m_deltaPotential_RT);
+            if (m_eChemForm == "Marcus") {
+                // If Marcus kinetics is being used, adjust beta by the reorganization 
+                // energy m_lambdaMarcus
+                beta_tmp += (m_etaF / 4. / m_lambdaMarcus); //m_etaF / 4 / m_lambdaMarcus;
+            }
+            correction = exp(-beta_tmp * m_deltaPotential_RT);
         }
 
         // Update correction if exchange current density formulation format is used.
-        if (m_exchangeCurrentDensityFormulation) {
+        if (m_eChemForm == "Butler-Volmer" || m_exchangeCurrentDensityFormulation) {
             // Comment preserved from previous implementation:
             // We need to have the straight chemical reaction rate constant to
             // come out of this calculation.
             double tmp = exp(-m_beta * m_deltaGibbs0_RT);
             tmp /= m_prodStandardConcentrations * Faraday;
             correction *= tmp;
+        } else if (m_eChemForm == "Marcus") {
+            // If Marcus kinetics is being used, adjust the correction value.
+            double tmp = exp(-m_lambda_RT / 4.) * exp(-beta_tmp * m_deltaGibbs_RT);
+            tmp /= m_prodStandardConcentrations * Faraday;
+            correction *= tmp;
         }
         return correction;
     }
@@ -204,6 +222,31 @@ class InterfaceRateBase
         return NAN;
     }
 
+    //! Return the Marcus theory reorganization energy parameter
+    double lambdaMarcus() const {
+        if (m_chargeTransfer) {
+            return m_lambdaMarcus;
+        }
+        return NAN;
+    }
+
+    //! Return the overpotential times Faraday's constant parameter
+    double etaF() {
+        if (m_chargeTransfer) {
+            return m_etaF;
+        }
+        return NAN;
+    }
+
+    //! Return Marcus theory reorganization energy parameter
+    //! Divided by the Gas Constant and temperature
+    double lambda_RT() {
+        if (m_chargeTransfer) {
+            return m_lambda_RT;
+        }
+        return NAN;
+    }
+
     //! Return site density [kmol/m^2]
     /*!
      *  @warning  This method is an experimental part of the %Cantera API and
@@ -233,6 +276,11 @@ class InterfaceRateBase
     double m_mcov; //!< Coverage term in reaction rate
     bool m_chargeTransfer; //!< Boolean indicating use of electrochemistry
     bool m_exchangeCurrentDensityFormulation; //! Electrochemistry only
+    std::string m_eChemForm; //! String storing echem kinetics form
+    double m_lambdaMarcus; //! Marcus theory reorganization energy
+    double m_etaF; //! Overpotential times Faraday
+    double m_lambda_RT; //! Marcus theory reorganization energy over RT
+    double m_deltaGibbs_RT; //! Normalized Gibbs free energy change 
     double m_beta; //!< Forward value of apparent electrochemical transfer coefficient
     double m_deltaPotential_RT; //!< Normalized electric potential energy change
     double m_deltaGibbs0_RT; //!< Normalized standard state Gibbs free energy change

src/kinetics/InterfaceRate.cpp

@@ -77,6 +77,7 @@ bool InterfaceData::update(const ThermoPhase& phase, const Kinetics& kin)
             electricPotentials[n] = ph.electricPotential();
             ph.getPartialMolarEnthalpies(partialMolarEnthalpies.data() + start);
             ph.getStandardChemPotentials(standardChemPotentials.data() + start);
+            ph.getChemPotentials(chemPotentials.data() + start);
             size_t nsp = ph.nSpecies();
             for (size_t k = 0; k < nsp; k++) {
                 // only used for exchange current density formulation
@@ -102,6 +103,10 @@ InterfaceRateBase::InterfaceRateBase()
     , m_chargeTransfer(false)
     , m_exchangeCurrentDensityFormulation(false)
     , m_beta(0.5)
+    , m_lambdaMarcus(0.)
+    , m_etaF(NAN)
+    , m_lambda_RT(NAN)
+    , m_deltaGibbs_RT(NAN)
     , m_deltaPotential_RT(NAN)
     , m_deltaGibbs0_RT(NAN)
     , m_prodStandardConcentrations(NAN)
@@ -117,6 +122,12 @@ void InterfaceRateBase::setParameters(const AnyMap& node)
     if (node.hasKey("beta")) {
         m_beta = node["beta"].asDouble();
     }
+    if (node.hasKey("echem-kinetics-form")) {
+        m_eChemForm = node["echem-kinetics-form"].asString();
+        if (m_eChemForm == "Marcus" && node.hasKey("lambda")) {
+            m_lambdaMarcus = node["lambda"].asDouble();
+        }
+    }
     m_exchangeCurrentDensityFormulation = node.getBool(
         "exchange-current-density-formulation", false);
 }
@@ -132,6 +143,12 @@ void InterfaceRateBase::getParameters(AnyMap& node) const
         if (m_beta != 0.5) {
             node["beta"] = m_beta;
         }
+        if (m_eChemForm != "") {
+            node["echem-kinetics-form"] = m_eChemForm;
+            if (m_eChemForm == "Marcus") {
+                node["lambda"] = m_lambdaMarcus;
+            }
+        }
         if (m_exchangeCurrentDensityFormulation) {
             node["exchange-current-density-formulation"] = true;
         }
@@ -238,26 +255,45 @@ void InterfaceRateBase::updateFromStruct(const InterfaceData& shared_data) {
     // Update change in electrical potential energy
     if (m_chargeTransfer) {
         m_deltaPotential_RT = 0.;
+        m_etaF = 0.;
         for (const auto& ch : m_netCharges) {
             m_deltaPotential_RT +=
                 shared_data.electricPotentials[ch.first] * ch.second;
         }
+        m_etaF += m_deltaPotential_RT;
         m_deltaPotential_RT /= GasConstant * shared_data.temperature;
     }
 
     // Update quantities used for exchange current density formulation
-    if (m_exchangeCurrentDensityFormulation) {
-        m_deltaGibbs0_RT = 0.;
+    if (m_eChemForm != "") {
         m_prodStandardConcentrations = 1.;
-        for (const auto& item : m_stoichCoeffs) {
-            m_deltaGibbs0_RT +=
-                shared_data.standardChemPotentials[item.first] * item.second;
-            if (item.second > 0.) {
-                m_prodStandardConcentrations *=
-                    shared_data.standardConcentrations[item.first];
+        if (m_eChemForm == "Butler-Volmer" || m_exchangeCurrentDensityFormulation) {
+            m_deltaGibbs0_RT = 0.;
+            for (const auto& item : m_stoichCoeffs) {
+                m_deltaGibbs0_RT +=
+                    shared_data.standardChemPotentials[item.first] * item.second;
+                if (item.second > 0.) {
+                    m_prodStandardConcentrations *=
+                        shared_data.standardConcentrations[item.first];
+                }
+            }
+            m_deltaGibbs0_RT /= GasConstant * shared_data.temperature;
+        }
+        else if (m_eChemForm == "Marcus") {
+            m_deltaGibbs_RT = 0.;
+            m_lambda_RT = 0.;
+            for (const auto& item : m_stoichCoeffs) {
+                m_deltaGibbs_RT +=
+                    shared_data.chemPotentials[item.first] * item.second;
+                if (item.second > 0.) {
+                    m_prodStandardConcentrations *=
+                        shared_data.standardConcentrations[item.first];
+                }
             }
+            m_etaF += m_deltaGibbs_RT;
+            m_deltaGibbs_RT /= GasConstant * shared_data.temperature;
+            m_lambda_RT = m_lambdaMarcus / GasConstant / shared_data.temperature;
         }
-        m_deltaGibbs0_RT /= GasConstant * shared_data.temperature;
     }
 }