From 4f0dfc3bdf09ae6cf9caf0b806af6fabf3002186 Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Wed, 15 Sep 2021 12:06:07 +0100
Subject: [PATCH] chemistryModel: Change of state variable to mass fraction

The chemistry model now solves a system of mass fractions, temperature
and pressure, rather than a system of concentrations, temperature and
pressure.

The new form now accounts for the change in reaction rate associated
with thermal expansion. Thermal expansion (or contraction) can dilute
(or concentrate) the species concentrations, thereby reducing (or
increasing) the reaction rates. Previously it was not possible to
include this term because it was not computationally feasible to
evaluate it in a system in which the state variable was concentration.

The reaction rate interface has been simplified with respect to the
generation of derivatives. Reactions are defined as k*C, where k is the
reaction rate and C is the product of concentrations (raised to their
stoichiomentric coefficients and/or specified powers). Reaction rate
classes now provide two logical functions governing the derivatives of
k; i.e., ddT (derivative w.r.t. temperature) and ddc (derivative w.r.t.
concentration). Previously the reaction rate interface was closely
related to the form of third-body reactions, which made it inconvenient
to implement rates that were not very third-body-like.

It is now possible to verify the implementations of the jacobian methods
by comparison with finite-difference based evaluations of the rate
methods. This has been done and a number of bugs have been found and
fixed in the reaction rate classes.
---
 src/ODE/ODESystem/ODESystem.C                 |  34 +-
 .../basicChemistryModel/basicChemistryModel.C |  21 +-
 .../basicChemistryModel/basicChemistryModel.H |  16 +
 .../chemistryModel/chemistryModel.C           | 481 +++++++++++-------
 .../chemistryModel/chemistryModel.H           |  54 +-
 .../chemistryModel/chemistryModelI.H          |  20 +-
 .../chemistryModel/tabulation/ISAT/ISAT.C     |  70 +--
 .../chemistryModel/tabulation/ISAT/ISAT.H     |   6 +-
 .../ISAT/chemPointISAT/chemPointISAT.H        |   2 +-
 .../chemistryTabulationMethod.C               |   2 +-
 .../chemistryTabulationMethod.H               |   6 +-
 .../chemistryTabulationMethodNew.C            |   2 +-
 .../noChemistryTabulation.C                   |   2 +-
 .../noChemistryTabulation.H                   |   2 +-
 .../chemistryModel/chemistrySolver/ode/ode.C  |   4 +-
 ...uxLimitedLangmuirHinshelwoodReactionRate.H |  19 +-
 ...xLimitedLangmuirHinshelwoodReactionRateI.H |  21 +-
 .../IrreversibleReaction.C                    |  20 +-
 .../IrreversibleReaction.H                    |  28 +-
 .../NonEquilibriumReversibleReaction.C        |  24 +-
 .../NonEquilibriumReversibleReaction.H        |  25 +-
 .../reaction/Reactions/Reaction/Reaction.C    | 354 ++++++-------
 .../reaction/Reactions/Reaction/Reaction.H    |  96 ++--
 .../Reactions/ReactionProxy/ReactionProxy.C   |  17 +-
 .../Reactions/ReactionProxy/ReactionProxy.H   |  21 +-
 .../ReversibleReaction/ReversibleReaction.C   |  26 +-
 .../ReversibleReaction/ReversibleReaction.H   |  25 +-
 .../ArrheniusReactionRate.H                   |  25 +-
 .../ArrheniusReactionRateI.H                  |  21 +-
 .../ChemicallyActivatedReactionRate.H         |  29 +-
 .../ChemicallyActivatedReactionRateI.H        | 170 ++++---
 .../FallOffReactionRate/FallOffReactionRate.H |  30 +-
 .../FallOffReactionRateI.H                    | 132 +++--
 .../JanevReactionRate/JanevReactionRate.H     |  23 +-
 .../JanevReactionRate/JanevReactionRateI.H    |  21 +-
 .../LandauTellerReactionRate.H                |  23 +-
 .../LandauTellerReactionRateI.H               |  21 +-
 .../LangmuirHinshelwoodReactionRate.H         |  23 +-
 .../LangmuirHinshelwoodReactionRateI.H        |  73 +--
 .../MichaelisMentenReactionRate.H             |  25 +-
 .../MichaelisMentenReactionRateI.H            |  28 +-
 .../LindemannFallOffFunction.H                |  14 +-
 .../LindemannFallOffFunctionI.H               |   6 +-
 .../SRIFallOffFunction/SRIFallOffFunction.H   |  14 +-
 .../SRIFallOffFunction/SRIFallOffFunctionI.H  |  57 ++-
 .../TroeFallOffFunction/TroeFallOffFunction.H |  14 +-
 .../TroeFallOffFunctionI.H                    | 157 +++---
 .../powerSeries/powerSeriesReactionRate.H     |  23 +-
 .../powerSeries/powerSeriesReactionRateI.H    |  21 +-
 .../thirdBodyArrheniusReactionRate.H          |  27 +-
 .../thirdBodyArrheniusReactionRateI.H         |  50 +-
 .../thirdBodyEfficiencies.H                   |  12 +-
 .../thirdBodyEfficienciesI.H                  |  17 +-
 .../reaction/specieCoeffs/specieCoeffs.C      |   4 +-
 .../reaction/specieCoeffs/specieCoeffs.H      | 121 +++--
 .../reaction/specieExponent/specieExponent.H  | 154 ++++++
 .../reaction/specieExponent/specieExponentI.H | 174 +++++++
 57 files changed, 1595 insertions(+), 1312 deletions(-)
 create mode 100644 src/thermophysicalModels/specie/reaction/specieExponent/specieExponent.H
 create mode 100644 src/thermophysicalModels/specie/reaction/specieExponent/specieExponentI.H

diff --git a/src/ODE/ODESystem/ODESystem.C b/src/ODE/ODESystem/ODESystem.C
index 54ee308b96..0e4e465baf 100644
--- a/src/ODE/ODESystem/ODESystem.C
+++ b/src/ODE/ODESystem/ODESystem.C
@@ -57,8 +57,15 @@ void Foam::ODESystem::check
     jacobian(x, y, li, dfdx1, d2fdxdyAnalytic);
 
     // Compare derivatives
-    Info<< "[derivatives] dfdx = " << dfdx0 << nl;
-    Info<< "[   jacobian] dfdx = " << dfdx1 << nl;
+    Info<< "[derivatives] dfdx = ( ";
+    forAll(dfdx0, i) { Info<< dfdx0[i] << ' '; }
+    Info<< ")" << nl;
+    Info<< "[   jacobian] dfdx = ( ";
+    forAll(dfdx1, i) { Info<< dfdx1[i] << ' '; }
+    Info<< ")" << nl;
+    Info<< "[      ratio] dfdx = ( ";
+    forAll(dfdx1, i) { Info<< dfdx1[i]/stabilise(dfdx0[i], rootVSmall) << ' '; }
+    Info<< ")" << nl;
 
     // Construct a Jacobian using the finite differences and the derivatives
     // method
@@ -83,10 +90,25 @@ void Foam::ODESystem::check
 
     for (label i = 0; i < nEqns(); ++ i)
     {
-        Info<< "[derivatives] d2fdxdy[" << i << "] = "
-            << UList<scalar>(d2fdxdyFiniteDifference[i], nEqns()) << nl;
-        Info<< "[   jacobian] d2fdxdy[" << i << "] = "
-            << UList<scalar>(d2fdxdyAnalytic[i], nEqns()) << nl;
+        UList<scalar> FD(d2fdxdyFiniteDifference[i], nEqns());
+        Info<< "[derivatives] d2fdxdy[" << i << "] = ( ";
+        forAll(FD, i) { Info<< FD[i] << ' '; }
+        Info<< ")" << nl;
+    }
+    for (label i = 0; i < nEqns(); ++ i)
+    {
+        UList<scalar> A(d2fdxdyAnalytic[i], nEqns());
+        Info<< "[   jacobian] d2fdxdy[" << i << "] = ( ";
+        forAll(A, i) { Info<< A[i] << ' '; }
+        Info<< ")" << nl;
+    }
+    for (label i = 0; i < nEqns(); ++ i)
+    {
+        UList<scalar> FD(d2fdxdyFiniteDifference[i], nEqns());
+        UList<scalar> A(d2fdxdyAnalytic[i], nEqns());
+        Info<< "[      ratio] d2fdxdy[" << i << "] = ( ";
+        forAll(A, i) { Info<< A[i]/stabilise(FD[i], rootVSmall) << ' '; }
+        Info<< ")" << nl;
     }
 }
 
diff --git a/src/thermophysicalModels/chemistryModel/basicChemistryModel/basicChemistryModel.C b/src/thermophysicalModels/chemistryModel/basicChemistryModel/basicChemistryModel.C
index 12a9b113bf..b704bb790c 100644
--- a/src/thermophysicalModels/chemistryModel/basicChemistryModel/basicChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/basicChemistryModel/basicChemistryModel.C
@@ -25,7 +25,25 @@ License
 
 #include "basicChemistryModel.H"
 
-/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
+/* * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * */
+
+namespace Foam
+{
+    template<>
+    const char* NamedEnum<basicChemistryModel::jacobianType, 2>::names[] =
+    {
+        "fast",
+        "exact"
+    };
+}
+
+
+const Foam::NamedEnum
+<
+    Foam::basicChemistryModel::jacobianType,
+    2
+> Foam::basicChemistryModel::jacobianTypeNames_;
+
 
 namespace Foam
 {
@@ -33,6 +51,7 @@ namespace Foam
     defineRunTimeSelectionTable(basicChemistryModel, thermo);
 }
 
+
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
 void Foam::basicChemistryModel::correct()
diff --git a/src/thermophysicalModels/chemistryModel/basicChemistryModel/basicChemistryModel.H b/src/thermophysicalModels/chemistryModel/basicChemistryModel/basicChemistryModel.H
index c87cb1e9ff..7550254be2 100644
--- a/src/thermophysicalModels/chemistryModel/basicChemistryModel/basicChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/basicChemistryModel/basicChemistryModel.H
@@ -54,6 +54,22 @@ class basicChemistryModel
 :
     public IOdictionary
 {
+public:
+
+    // Public Enumerations
+
+        //- Enumeration for the type of Jacobian to be calculated by the
+        //  chemistry model
+        enum class jacobianType
+        {
+            fast,
+            exact
+        };
+
+        //- Jacobian type names
+        static const NamedEnum<jacobianType, 2> jacobianTypeNames_;
+
+
 protected:
 
     // Protected data
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C
index 0d680b50ca..92e0f95329 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C
@@ -39,16 +39,22 @@ Foam::chemistryModel<ThermoType>::chemistryModel
     basicChemistryModel(thermo),
     ODESystem(),
     log_(this->lookupOrDefault("log", false)),
-    Y_(this->thermo().composition().Y()),
+    jacobianType_
+    (
+        this->found("jacobian")
+      ? jacobianTypeNames_.read(this->lookup("jacobian"))
+      : jacobianType::fast
+    ),
+    Yvf_(this->thermo().composition().Y()),
     mixture_(refCast<const multiComponentMixture<ThermoType>>(this->thermo())),
     specieThermos_(mixture_.specieThermos()),
     reactions_(mixture_.species(), specieThermos_, this->mesh(), *this),
-    nSpecie_(Y_.size()),
-    nReaction_(reactions_.size()),
-    Treact_(basicChemistryModel::template lookupOrDefault<scalar>("Treact", 0)),
+    nSpecie_(Yvf_.size()),
     RR_(nSpecie_),
+    Y_(nSpecie_),
     c_(nSpecie_),
-    dcdt_(nSpecie_),
+    YTpWork_(scalarField(nSpecie_ + 2)),
+    YTpYTpWork_(scalarSquareMatrix(nSpecie_ + 2)),
     cTos_(nSpecie_, -1),
     sToc_(nSpecie_),
     mechRedPtr_
@@ -81,7 +87,7 @@ Foam::chemistryModel<ThermoType>::chemistryModel
             (
                 IOobject
                 (
-                    "RR." + Y_[fieldi].name(),
+                    "RR." + Yvf_[fieldi].name(),
                     this->mesh().time().timeName(),
                     this->mesh(),
                     IOobject::NO_READ,
@@ -94,7 +100,7 @@ Foam::chemistryModel<ThermoType>::chemistryModel
     }
 
     Info<< "chemistryModel: Number of species = " << nSpecie_
-        << " and reactions = " << nReaction_ << endl;
+        << " and reactions = " << nReaction() << endl;
 
     // When the mechanism reduction method is used, the 'active' flag for every
     // species should be initialised (by default 'active' is true)
@@ -102,11 +108,11 @@ Foam::chemistryModel<ThermoType>::chemistryModel
     {
         const basicSpecieMixture& composition = this->thermo().composition();
 
-        forAll(Y_, i)
+        forAll(Yvf_, i)
         {
             typeIOobject<volScalarField> header
             (
-                Y_[i].name(),
+                Yvf_[i].name(),
                 this->mesh().time().timeName(),
                 this->mesh(),
                 IOobject::NO_READ
@@ -137,101 +143,95 @@ Foam::chemistryModel<ThermoType>::~chemistryModel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class ThermoType>
-void Foam::chemistryModel<ThermoType>::omega
-(
-    const scalar p,
-    const scalar T,
-    const scalarField& c, // Contains all species even when mechRed is active
-    const label li,
-    scalarField& dcdt
-) const
-{
-    if (mechRedActive_)
-    {
-        forAll(sToc_, si)
-        {
-            dcdt_[sToc_[si]] = 0;
-        }
-
-        forAll(reactions_, i)
-        {
-            if (!mechRed_.reactionDisabled(i))
-            {
-                const Reaction<ThermoType>& R = reactions_[i];
-
-                R.omega(p, T, c, li, dcdt_);
-            }
-        }
-
-        forAll(sToc_, si)
-        {
-            dcdt[si] = dcdt_[sToc_[si]];
-        }
-    }
-    else
-    {
-        dcdt = Zero;
-
-        forAll(reactions_, i)
-        {
-            const Reaction<ThermoType>& R = reactions_[i];
-            R.omega(p, T, c, li, dcdt);
-        }
-    }
-}
-
-
 template<class ThermoType>
 void Foam::chemistryModel<ThermoType>::derivatives
 (
     const scalar time,
-    const scalarField& c,
+    const scalarField& YTp,
     const label li,
-    scalarField& dcdt
+    scalarField& dYTpdt
 ) const
 {
-    const scalar T = c[nSpecie_];
-    const scalar p = c[nSpecie_ + 1];
-
     if (mechRedActive_)
     {
         forAll(sToc_, i)
         {
-            c_[sToc_[i]] = max(c[i], 0);
+            Y_[sToc_[i]] = max(YTp[i], 0);
         }
     }
     else
     {
-        forAll(c_, i)
+        forAll(Y_, i)
         {
-            c_[i] = max(c[i], 0);
+            Y_[i] = max(YTp[i], 0);
         }
     }
 
+    const scalar T = YTp[nSpecie_];
+    const scalar p = YTp[nSpecie_ + 1];
+
+    // Evaluate the mixture density
+    scalar rhoM = 0;
+    for (label i=0; i<Y_.size(); i++)
+    {
+        rhoM += Y_[i]/specieThermos_[i].rho(p, T);
+    }
+    rhoM = 1/rhoM;
+
+    // Evaluate the concentrations
+    for (label i=0; i<Y_.size(); i ++)
+    {
+        c_[i] = rhoM/specieThermos_[i].W()*Y_[i];
+    }
+
     // Evaluate contributions from reactions
-    omega(p, T, c_, li, dcdt);
+    dYTpdt = Zero;
+    forAll(reactions_, ri)
+    {
+        if (!mechRed_.reactionDisabled(ri))
+        {
+            reactions_[ri].dNdtByV
+            (
+                p,
+                T,
+                c_,
+                li,
+                dYTpdt,
+                mechRedActive_,
+                cTos_,
+                0
+            );
+        }
+    }
+
+    // Reactions return dNdtByV, so we need to convert the result to dYdt
+    for (label i=0; i<nSpecie_; i++)
+    {
+        const scalar WiByrhoM = specieThermos_[sToc(i)].W()/rhoM;
+        scalar& dYidt = dYTpdt[i];
+        dYidt *= WiByrhoM;
+    }
 
     // Evaluate the effect on the thermodynamic system ...
 
-    // c*Cp
-    scalar ccp = 0;
-    for (label i=0; i<c_.size(); i++)
+    // Evaluate the mixture Cp
+    scalar CpM = 0;
+    for (label i=0; i<Y_.size(); i++)
     {
-        ccp += c_[i]*specieThermos_[i].cp(p, T);
+        CpM += Y_[i]*specieThermos_[i].Cp(p, T);
     }
 
     // dT/dt
-    scalar& dTdt = dcdt[nSpecie_];
-    dTdt = 0;
+    scalar& dTdt = dYTpdt[nSpecie_];
     for (label i=0; i<nSpecie_; i++)
     {
-        dTdt -= dcdt[i]*specieThermos_[sToc(i)].ha(p, T);
+        dTdt -= dYTpdt[i]*specieThermos_[sToc(i)].Ha(p, T);
     }
-    dTdt /= ccp;
+    dTdt /= CpM;
 
     // dp/dt = 0 (pressure is assumed constant)
-    dcdt[nSpecie_ + 1] = 0;
+    scalar& dpdt = dYTpdt[nSpecie_ + 1];
+    dpdt = 0;
 }
 
 
@@ -239,125 +239,217 @@ template<class ThermoType>
 void Foam::chemistryModel<ThermoType>::jacobian
 (
     const scalar t,
-    const scalarField& c,
+    const scalarField& YTp,
     const label li,
-    scalarField& dcdt,
+    scalarField& dYTpdt,
     scalarSquareMatrix& J
 ) const
 {
-    // If the mechanism reduction is active, the computed Jacobian
-    // is compact (size of the reduced set of species)
-    // but according to the information of the complete set
-    // (i.e. for the third-body efficiencies)
-
     if (mechRedActive_)
     {
         forAll(sToc_, i)
         {
-            c_[sToc_[i]] = max(c[i], 0);
+            Y_[sToc_[i]] = max(YTp[i], 0);
         }
     }
     else
     {
         forAll(c_, i)
         {
-            c_[i] = max(c[i], 0);
+            Y_[i] = max(YTp[i], 0);
         }
     }
 
-    const scalar T = c[nSpecie_];
-    const scalar p = c[nSpecie_ + 1];
+    const scalar T = YTp[nSpecie_];
+    const scalar p = YTp[nSpecie_ + 1];
 
-    dcdt = Zero;
-    J = Zero;
+    // Evaluate the specific volumes and mixture density
+    scalarField& v = YTpWork_[0];
+    for (label i=0; i<Y_.size(); i++)
+    {
+        v[i] = 1/specieThermos_[i].rho(p, T);
+    }
+    scalar rhoM = 0;
+    for (label i=0; i<Y_.size(); i++)
+    {
+        rhoM += Y_[i]*v[i];
+    }
+    rhoM = 1/rhoM;
+
+    // Evaluate the concentrations
+    for (label i=0; i<Y_.size(); i ++)
+    {
+        c_[i] = rhoM/specieThermos_[i].W()*Y_[i];
+    }
+
+    // Evaluate the derivatives of concentration w.r.t. mass fraction
+    scalarSquareMatrix& dcdY = YTpYTpWork_[0];
+    for (label i=0; i<nSpecie_; i++)
+    {
+        const scalar rhoMByWi = rhoM/specieThermos_[sToc(i)].W();
+        switch (jacobianType_)
+        {
+            case jacobianType::fast:
+                {
+                    dcdY(i, i) = rhoMByWi;
+                }
+                break;
+            case jacobianType::exact:
+                for (label j=0; j<nSpecie_; j++)
+                {
+                    dcdY(i, j) =
+                        rhoMByWi*((i == j) - rhoM*v[sToc(j)]*Y_[sToc(i)]);
+                }
+                break;
+        }
+    }
+
+    // Evaluate the mixture thermal expansion coefficient
+    scalar alphavM = 0;
+    for (label i=0; i<Y_.size(); i++)
+    {
+        alphavM += Y_[i]*rhoM*v[i]*specieThermos_[i].alphav(p, T);
+    }
 
     // Evaluate contributions from reactions
+    dYTpdt = Zero;
+    scalarSquareMatrix& ddNdtByVdcTp = YTpYTpWork_[1];
+    for (label i=0; i<nSpecie_ + 2; i++)
+    {
+        for (label j=0; j<nSpecie_ + 2; j++)
+        {
+            ddNdtByVdcTp[i][j] = 0;
+        }
+    }
     forAll(reactions_, ri)
     {
         if (!mechRed_.reactionDisabled(ri))
         {
-            const Reaction<ThermoType>& R = reactions_[ri];
-            scalar omegaI, kfwd, kbwd;
-            R.dwdc
+            reactions_[ri].ddNdtByVdcTp
             (
                 p,
                 T,
                 c_,
                 li,
-                J,
-                dcdt,
-                omegaI,
-                kfwd,
-                kbwd,
-                mechRedActive_,
-                cTos_
-            );
-            R.dwdT
-            (
-                p,
-                T,
-                c_,
-                li,
-                omegaI,
-                kfwd,
-                kbwd,
-                J,
+                dYTpdt,
+                ddNdtByVdcTp,
                 mechRedActive_,
                 cTos_,
-                nSpecie_
+                0,
+                nSpecie_,
+                YTpWork_[1],
+                YTpWork_[2]
             );
         }
     }
 
+    // Reactions return dNdtByV, so we need to convert the result to dYdt
+    for (label i=0; i<nSpecie_; i++)
+    {
+        const scalar WiByrhoM = specieThermos_[sToc(i)].W()/rhoM;
+        scalar& dYidt = dYTpdt[i];
+        dYidt *= WiByrhoM;
+
+        for (label j=0; j<nSpecie_; j++)
+        {
+            scalar ddNidtByVdYj = 0;
+            switch (jacobianType_)
+            {
+                case jacobianType::fast:
+                    {
+                        const scalar ddNidtByVdcj = ddNdtByVdcTp(i, j);
+                        ddNidtByVdYj = ddNidtByVdcj*dcdY(j, j);
+                    }
+                    break;
+                case jacobianType::exact:
+                    for (label k=0; k<nSpecie_; k++)
+                    {
+                        const scalar ddNidtByVdck = ddNdtByVdcTp(i, k);
+                        ddNidtByVdYj += ddNidtByVdck*dcdY(k, j);
+                    }
+                    break;
+            }
+
+            scalar& ddYidtdYj = J(i, j);
+            ddYidtdYj = WiByrhoM*ddNidtByVdYj + rhoM*v[sToc(j)]*dYidt;
+        }
+
+        scalar ddNidtByVdT = ddNdtByVdcTp(i, nSpecie_);
+        for (label j=0; j<nSpecie_; j++)
+        {
+            const scalar ddNidtByVdcj = ddNdtByVdcTp(i, j);
+            ddNidtByVdT -= ddNidtByVdcj*c_[sToc(j)]*alphavM;
+        }
+
+        scalar& ddYidtdT = J(i, nSpecie_);
+        ddYidtdT = WiByrhoM*ddNidtByVdT + alphavM*dYidt;
+
+        scalar& ddYidtdp = J(i, nSpecie_ + 1);
+        ddYidtdp = 0;
+    }
+
     // Evaluate the effect on the thermodynamic system ...
 
-    // c*Cp
-    scalar ccp = 0, dccpdT = 0;
-    forAll(c_, i)
+    // Evaluate the mixture Cp and its derivative
+    scalarField& Cp = YTpWork_[3];
+    scalar CpM = 0, dCpMdT = 0;
+    for (label i=0; i<Y_.size(); i++)
     {
-        ccp += c_[i]*specieThermos_[i].cp(p, T);
-        dccpdT += c_[i]*specieThermos_[i].dcpdT(p, T);
+        Cp[i] = specieThermos_[i].Cp(p, T);
+        CpM += Y_[i]*Cp[i];
+        dCpMdT += Y_[i]*specieThermos_[i].dCpdT(p, T);
     }
 
     // dT/dt
-    scalar& dTdt = dcdt[nSpecie_];
+    scalarField& Ha = YTpWork_[4];
+    scalar& dTdt = dYTpdt[nSpecie_];
     for (label i=0; i<nSpecie_; i++)
     {
-        dTdt -= dcdt[i]*specieThermos_[sToc(i)].ha(p, T);
+        Ha[sToc(i)] = specieThermos_[sToc(i)].Ha(p, T);
+        dTdt -= dYTpdt[i]*Ha[sToc(i)];
     }
-    dTdt /= ccp;
+    dTdt /= CpM;
 
     // dp/dt = 0 (pressure is assumed constant)
+    scalar& dpdt = dYTpdt[nSpecie_ + 1];
+    dpdt = 0;
 
-    // d(dTdt)/dc
-    for (label i = 0; i < nSpecie_; i++)
+    // d(dTdt)/dY
+    for (label i=0; i<nSpecie_; i++)
     {
-        scalar& d2Tdtdci = J(nSpecie_, i);
-        for (label j = 0; j < nSpecie_; j++)
+        scalar& ddTdtdYi = J(nSpecie_, i);
+        ddTdtdYi = 0;
+        for (label j=0; j<nSpecie_; j++)
         {
-            const scalar d2cjdtdci = J(j, i);
-            d2Tdtdci -= d2cjdtdci*specieThermos_[sToc(j)].ha(p, T);
+            const scalar ddYjdtdYi = J(j, i);
+            ddTdtdYi -= ddYjdtdYi*Ha[sToc(j)];
         }
-        d2Tdtdci -= specieThermos_[sToc(i)].cp(p, T)*dTdt;
-        d2Tdtdci /= ccp;
+        ddTdtdYi -= Cp[sToc(i)]*dTdt;
+        ddTdtdYi /= CpM;
     }
 
     // d(dTdt)/dT
-    scalar& d2TdtdT = J(nSpecie_, nSpecie_);
-    for (label i = 0; i < nSpecie_; i++)
+    scalar& ddTdtdT = J(nSpecie_, nSpecie_);
+    ddTdtdT = 0;
+    for (label i=0; i<nSpecie_; i++)
     {
-        const scalar d2cidtdT = J(i, nSpecie_);
-        const label si = sToc(i);
-        d2TdtdT -=
-            dcdt[i]*specieThermos_[si].cp(p, T)
-          + d2cidtdT*specieThermos_[si].ha(p, T);
+        const scalar dYidt = dYTpdt[i];
+        const scalar ddYidtdT = J(i, nSpecie_);
+        ddTdtdT -= dYidt*Cp[sToc(i)] + ddYidtdT*Ha[sToc(i)];
     }
-    d2TdtdT -= dTdt*dccpdT;
-    d2TdtdT /= ccp;
+    ddTdtdT -= dTdt*dCpMdT;
+    ddTdtdT /= CpM;
 
-    // d(dpdt)/dc = 0 (pressure is assumed constant)
+    // d(dTdt)/dp = 0 (pressure is assumed constant)
+    scalar& ddTdtdp = J(nSpecie_, nSpecie_ + 1);
+    ddTdtdp = 0;
 
-    // d(dpdt)/dT = 0 (pressure is assumed constant)
+    // d(dpdt)/dYiTp = 0 (pressure is assumed constant)
+    for (label i=0; i<nSpecie_ + 2; i++)
+    {
+        scalar& ddpdtdYiTp = J(nSpecie_ + 1, i);
+        ddpdtdYiTp = 0;
+    }
 }
 
 
@@ -395,7 +487,7 @@ Foam::chemistryModel<ThermoType>::tc() const
 
             for (label i=0; i<nSpecie_; i++)
             {
-                c_[i] = rhoi*Y_[i][celli]/specieThermos_[i].W();
+                c_[i] = rhoi*Yvf_[i][celli]/specieThermos_[i].W();
             }
 
             // A reaction's rate scale is calculated as it's molar
@@ -470,7 +562,7 @@ Foam::chemistryModel<ThermoType>::Qdot() const
 
         scalarField& Qdot = tQdot.ref();
 
-        forAll(Y_, i)
+        forAll(Yvf_, i)
         {
             forAll(Qdot, celli)
             {
@@ -522,7 +614,7 @@ Foam::chemistryModel<ThermoType>::calculateRR
 
         for (label i=0; i<nSpecie_; i++)
         {
-            const scalar Yi = Y_[i][celli];
+            const scalar Yi = Yvf_[i][celli];
             c_[i] = rhoi*Yi/specieThermos_[i].W();
         }
 
@@ -566,6 +658,8 @@ void Foam::chemistryModel<ThermoType>::calculate()
     const scalarField& T = this->thermo().T();
     const scalarField& p = this->thermo().p();
 
+    scalarField& dNdtByV = YTpWork_[0];
+
     reactionEvaluationScope scope(*this);
 
     forAll(rho, celli)
@@ -576,15 +670,33 @@ void Foam::chemistryModel<ThermoType>::calculate()
 
         for (label i=0; i<nSpecie_; i++)
         {
-            const scalar Yi = Y_[i][celli];
+            const scalar Yi = Yvf_[i][celli];
             c_[i] = rhoi*Yi/specieThermos_[i].W();
         }
 
-        omega(pi, Ti, c_, celli, dcdt_);
+        dNdtByV = Zero;
 
-        for (label i=0; i<nSpecie_; i++)
+        forAll(reactions_, ri)
         {
-            RR_[i][celli] = dcdt_[i]*specieThermos_[i].W();
+            if (!mechRed_.reactionDisabled(ri))
+            {
+                reactions_[ri].dNdtByV
+                (
+                    pi,
+                    Ti,
+                    c_,
+                    celli,
+                    dNdtByV,
+                    mechRedActive_,
+                    cTos_,
+                    0
+                );
+            }
+        }
+
+        for (label i=0; i<mechRed_.nActiveSpecies(); i++)
+        {
+            RR_[sToc(i)][celli] = dNdtByV[i]*specieThermos_[sToc(i)].W();
         }
     }
 }
@@ -615,35 +727,41 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
         return deltaTMin;
     }
 
-    tmp<volScalarField> trho(this->thermo().rho0());
-    const scalarField& rho = trho();
+    tmp<volScalarField> trhovf(this->thermo().rho());
+    const volScalarField& rhovf = trhovf();
+    tmp<volScalarField> trho0vf(this->thermo().rho0());
+    const volScalarField& rho0vf = trho0vf();
 
-    const scalarField& T = this->thermo().T().oldTime();
-    const scalarField& p = this->thermo().p().oldTime();
+    const volScalarField& T0vf = this->thermo().T().oldTime();
+    const volScalarField& p0vf = this->thermo().p().oldTime();
 
     reactionEvaluationScope scope(*this);
 
-    scalarField c0(nSpecie_);
+    scalarField Y0(nSpecie_);
 
     // Composition vector (Yi, T, p, deltaT)
     scalarField phiq(nEqns() + 1);
     scalarField Rphiq(nEqns() + 1);
 
-    forAll(rho, celli)
+    forAll(rho0vf, celli)
     {
-        const scalar rhoi = rho[celli];
-        scalar pi = p[celli];
-        scalar Ti = T[celli];
+        const scalar rho = rhovf[celli];
+        const scalar rho0 = rho0vf[celli];
+
+        scalar p = p0vf[celli];
+        scalar T = T0vf[celli];
 
         for (label i=0; i<nSpecie_; i++)
         {
-            c_[i] =
-                rhoi*Y_[i].oldTime()[celli]/specieThermos_[i].W();
-            c0[i] = c_[i];
-            phiq[i] = Y_[i].oldTime()[celli];
+            Y_[i] = Y0[i] = Yvf_[i].oldTime()[celli];
         }
-        phiq[nSpecie()] = Ti;
-        phiq[nSpecie() + 1] = pi;
+
+        for (label i=0; i<nSpecie_; i++)
+        {
+            phiq[i] = Yvf_[i].oldTime()[celli];
+        }
+        phiq[nSpecie()] = T;
+        phiq[nSpecie() + 1] = p;
         phiq[nSpecie() + 2] = deltaT[celli];
 
         // Initialise time progress
@@ -660,8 +778,10 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
             // Retrieved solution stored in Rphiq
             for (label i=0; i<nSpecie(); i++)
             {
-                c_[i] = rhoi*Rphiq[i]/specieThermos_[i].W();
+                Y_[i] = Rphiq[i];
             }
+            T = Rphiq[nSpecie()];
+            p = Rphiq[nSpecie() + 1];
         }
         // This position is reached when tabulation is not used OR
         // if the solution is not retrieved.
@@ -671,17 +791,29 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
         {
             if (mechRedActive_)
             {
+                for (label i=0; i<nSpecie_; i++)
+                {
+                    const scalar Yi = Yvf_[i][celli];
+                    c_[i] = rho0*Yi/specieThermos_[i].W();
+                }
+
                 // Reduce mechanism change the number of species (only active)
                 mechRed_.reduceMechanism
                 (
-                    pi,
-                    Ti,
+                    p,
+                    T,
                     c_,
                     sc_,
                     cTos_,
                     sToc_,
                     celli
                 );
+
+                sY_.setSize(nSpecie_);
+                for (label i=0; i<nSpecie_; i++)
+                {
+                    sY_[i] = specieThermos_[sToc(i)].W()*sc_[i]/rho0;
+                }
             }
 
             if (log_)
@@ -699,9 +831,9 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
                     // Solve the reduced set of ODE
                     solve
                     (
-                        pi,
-                        Ti,
-                        sc_,
+                        p,
+                        T,
+                        sY_,
                         celli,
                         dt,
                         deltaTChem_[celli]
@@ -709,13 +841,12 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
 
                     for (label i=0; i<mechRed_.nActiveSpecies(); i++)
                     {
-                        c_[sToc_[i]] = sc_[i];
+                        Y_[sToc_[i]] = sY_[i];
                     }
                 }
                 else
                 {
-                    solve
-                    (pi, Ti, c_, celli, dt, deltaTChem_[celli]);
+                    solve(p, T, Y_, celli, dt, deltaTChem_[celli]);
                 }
                 timeLeft -= dt;
             }
@@ -729,15 +860,15 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
             // the stored points (either expand or add)
             if (tabulation_.tabulates())
             {
-                forAll(c_, i)
+                forAll(Y_, i)
                 {
-                    Rphiq[i] = c_[i]/rhoi*specieThermos_[i].W();
+                    Rphiq[i] = Y_[i];
                 }
-                Rphiq[Rphiq.size()-3] = Ti;
-                Rphiq[Rphiq.size()-2] = pi;
+                Rphiq[Rphiq.size()-3] = T;
+                Rphiq[Rphiq.size()-2] = p;
                 Rphiq[Rphiq.size()-1] = deltaT[celli];
 
-                tabulation_.add(phiq, Rphiq, celli, rhoi, deltaT[celli]);
+                tabulation_.add(phiq, Rphiq, celli, rho0, deltaT[celli]);
             }
 
             // When operations are done and if mechanism reduction is active,
@@ -745,19 +876,17 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
             // to the total number of species (stored in the mechRed object)
             if (mechRedActive_)
             {
-                nSpecie_ = mechRed_.nSpecie();
+                setNSpecie(mechRed_.nSpecie());
             }
-            deltaTMin = min(deltaTChem_[celli], deltaTMin);
 
-            deltaTChem_[celli] =
-                min(deltaTChem_[celli], deltaTChemMax_);
+            deltaTMin = min(deltaTChem_[celli], deltaTMin);
+            deltaTChem_[celli] = min(deltaTChem_[celli], deltaTChemMax_);
         }
 
         // Set the RR vector (used in the solver)
         for (label i=0; i<nSpecie_; i++)
         {
-            RR_[i][celli] =
-                (c_[i] - c0[i])*specieThermos_[i].W()/deltaT[celli];
+            RR_[i][celli] = (Y_[i]*rho - Y0[i]*rho0)/deltaT[celli];
         }
     }
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.H
index 5ca6de3029..c5c5a9a470 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.H
@@ -130,8 +130,11 @@ class chemistryModel
         //- Switch to select performance logging
         Switch log_;
 
+        //- Type of the Jacobian to be calculated
+        const jacobianType jacobianType_;
+
         //- Reference to the field of specie mass fractions
-        const PtrList<volScalarField>& Y_;
+        const PtrList<volScalarField>& Yvf_;
 
         //- Reference to the multi component mixture
         const multiComponentMixture<ThermoType>& mixture_;
@@ -145,24 +148,27 @@ class chemistryModel
         //- Number of species
         label nSpecie_;
 
-        //- Number of reactions
-        label nReaction_;
-
-        //- Temperature below which the reaction rates are assumed 0
-        scalar Treact_;
-
         //- List of reaction rate per specie [kg/m^3/s]
         PtrList<volScalarField::Internal> RR_;
 
+        //- Temporary mass fraction field
+        mutable scalarField Y_;
+
+        //- Temporary simplified mechanism mass fraction field
+        DynamicField<scalar> sY_;
+
         //- Temporary concentration field
         mutable scalarField c_;
 
-        //- Temporary rate-of-change of concentration field
-        mutable scalarField dcdt_;
-
         //- Temporary simplified mechanism concentration field
         DynamicField<scalar> sc_;
 
+        //- Specie-temperature-pressure workspace fields
+        mutable FixedList<scalarField, 5> YTpWork_;
+
+        //- Specie-temperature-pressure workspace matrices
+        mutable FixedList<scalarSquareMatrix, 2> YTpYTpWork_;
+
         //- Temporary map from complete to simplified concentration fields
         //  c -> sc
         List<label> cTos_;
@@ -241,22 +247,6 @@ public:
         //- The number of reactions
         virtual inline label nReaction() const;
 
-        //- Temperature below which the reaction rates are assumed 0
-        inline scalar Treact() const;
-
-        //- Temperature below which the reaction rates are assumed 0
-        inline scalar& Treact();
-
-        //- dc/dt = omega, rate of change in concentration, for each species
-        virtual void omega
-        (
-            const scalar p,
-            const scalar T,
-            const scalarField& c,
-            const label li,
-            scalarField& dcdt
-        ) const;
-
         //- Calculates the reaction rates
         virtual void calculate();
 
@@ -306,17 +296,17 @@ public:
             virtual void derivatives
             (
                 const scalar t,
-                const scalarField& c,
+                const scalarField& YTp,
                 const label li,
-                scalarField& dcdt
+                scalarField& dYTpdt
             ) const;
 
             virtual void jacobian
             (
                 const scalar t,
-                const scalarField& c,
+                const scalarField& YTp,
                 const label li,
-                scalarField& dcdt,
+                scalarField& dYTpdt,
                 scalarSquareMatrix& J
             ) const;
 
@@ -324,7 +314,7 @@ public:
             (
                 scalar& p,
                 scalar& T,
-                scalarField& c,
+                scalarField& Y,
                 const label li,
                 scalar& deltaT,
                 scalar& subDeltaT
@@ -332,7 +322,7 @@ public:
 
 
         //- Return a reference to the list of mass fraction fields
-        inline const PtrList<volScalarField>& Y();
+        inline const PtrList<volScalarField>& Y() const;
 
 
         // Mechanism reduction access functions
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModelI.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModelI.H
index 883c02b838..61e5c58b93 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModelI.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModelI.H
@@ -91,21 +91,7 @@ inline Foam::label Foam::chemistryModel<ThermoType>::nSpecie() const
 template<class ThermoType>
 inline Foam::label Foam::chemistryModel<ThermoType>::nReaction() const
 {
-    return nReaction_;
-}
-
-
-template<class ThermoType>
-inline Foam::scalar Foam::chemistryModel<ThermoType>::Treact() const
-{
-    return Treact_;
-}
-
-
-template<class ThermoType>
-inline Foam::scalar& Foam::chemistryModel<ThermoType>::Treact()
-{
-    return Treact_;
+    return reactions_.size();
 }
 
 
@@ -126,9 +112,9 @@ Foam::chemistryModel<ThermoType>::RR(const label i)
 
 template<class ThermoType>
 inline const Foam::PtrList<Foam::volScalarField>&
-Foam::chemistryModel<ThermoType>::Y()
+Foam::chemistryModel<ThermoType>::Y() const
 {
-    return this->Y_;
+    return this->Yvf_;
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/ISAT.C b/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/ISAT.C
index 9ccd4f76fe..0ead5424b2 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/ISAT.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/ISAT.C
@@ -32,7 +32,7 @@ template<class ThermoType>
 Foam::chemistryTabulationMethods::ISAT<ThermoType>::ISAT
 (
     const dictionary& chemistryProperties,
-    chemistryModel<ThermoType>& chemistry
+    const chemistryModel<ThermoType>& chemistry
 )
 :
     chemistryTabulationMethod<ThermoType>
@@ -352,16 +352,17 @@ void Foam::chemistryTabulationMethods::ISAT<ThermoType>::computeA
     const scalar dt
 )
 {
-    const label speciesNumber = chemistry_.nSpecie();
-    scalarField Rcq(chemistry_.nEqns() + 1);
-    for (label i=0; i<speciesNumber; i++)
+    const label nSpecie = chemistry_.nSpecie();
+
+    scalarField Rphiqs(chemistry_.nEqns() + 1);
+    for (label i=0; i<nSpecie; i++)
     {
         const label si = chemistry_.sToc(i);
-        Rcq[i] = rhoi*Rphiq[si]/chemistry_.specieThermos()[si].W();
+        Rphiqs[i] = Rphiq[si];
     }
-    Rcq[speciesNumber] = Rphiq[Rphiq.size() - 3];
-    Rcq[speciesNumber + 1] = Rphiq[Rphiq.size() - 2];
-    Rcq[speciesNumber + 2] = Rphiq[Rphiq.size() - 1];
+    Rphiqs[nSpecie] = Rphiq[Rphiq.size() - 3];
+    Rphiqs[nSpecie + 1] = Rphiq[Rphiq.size() - 2];
+    Rphiqs[nSpecie + 2] = Rphiq[Rphiq.size() - 1];
 
     // Aaa is computed implicitly,
     // A is given by A = C(psi0, t0+dt), where C is obtained through solving
@@ -372,62 +373,29 @@ void Foam::chemistryTabulationMethods::ISAT<ThermoType>::computeA
     // C(psi0,t0+dt)*(I-dt*J(psi(t0+dt))) = C(psi0, t0)
     // A = C(psi0,t0)/(I-dt*J(psi(t0+dt)))
     // where C(psi0,t0) = I
-    scalarField dcdt(speciesNumber + 2, Zero);
-    chemistry_.jacobian(runTime_.value(), Rcq, li, dcdt, A);
+    scalarField dYTpdt(nSpecie + 2, Zero);
+    chemistry_.jacobian(runTime_.value(), Rphiqs, li, dYTpdt, A);
 
-    // The jacobian is computed according to the molar concentration
-    // the following conversion allows the code to use A with mass fraction
-    for (label i=0; i<speciesNumber; i++)
+    // Inverse of I - dt*J(psi(t0 + dt))
+    for (label i=0; i<nSpecie + 2; i++)
     {
-        const label si = chemistry_.sToc(i);
-
-        for (label j=0; j<speciesNumber; j++)
+        for (label j=0; j<nSpecie + 2; j++)
         {
-            const label sj = chemistry_.sToc(j);
-
-            A(i, j) *=
-              -dt*chemistry_.specieThermos()[si].W()
-               /chemistry_.specieThermos()[sj].W();
+            A(i, j) *= -dt;
         }
-
         A(i, i) += 1;
-        // Columns for pressure and temperature
-        A(i, speciesNumber) *=
-            -dt*chemistry_.specieThermos()[si].W()/rhoi;
-        A(i, speciesNumber + 1) *=
-            -dt*chemistry_.specieThermos()[si].W()/rhoi;
     }
-
-    // For the temperature and pressure lines, ddc(dTdt)
-    // should be converted in ddY(dTdt)
-    for (label i=0; i<speciesNumber; i++)
-    {
-        const label si = chemistry_.sToc(i);
-
-        A(speciesNumber, i) *=
-            -dt*rhoi/chemistry_.specieThermos()[si].W();
-        A(speciesNumber + 1, i) *=
-            -dt*rhoi/chemistry_.specieThermos()[si].W();
-    }
-
-    A(speciesNumber, speciesNumber) = -dt*A(speciesNumber, speciesNumber) + 1;
-
-    A(speciesNumber + 1, speciesNumber + 1) =
-        -dt*A(speciesNumber + 1, speciesNumber + 1) + 1;
-
-    A(speciesNumber + 2, speciesNumber + 2) = 1;
-
-    // Inverse of (I-dt*J(psi(t0+dt)))
+    A(nSpecie + 2, nSpecie + 2) = 1;
     LUscalarMatrix LUA(A);
     LUA.inv(A);
 
     // After inversion, lines of p and T are set to 0 except diagonal.  This
     // avoid skewness of the ellipsoid of accuracy and potential issues in the
     // binary tree.
-    for (label i=0; i<speciesNumber; i++)
+    for (label i=0; i<nSpecie; i++)
     {
-        A(speciesNumber, i) = 0;
-        A(speciesNumber + 1, i) = 0;
+        A(nSpecie, i) = 0;
+        A(nSpecie + 1, i) = 0;
     }
 }
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/ISAT.H b/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/ISAT.H
index 067818e160..18ad424ee0 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/ISAT.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/ISAT.H
@@ -63,7 +63,7 @@ class ISAT
 
         const dictionary coeffsDict_;
 
-        chemistryModel<ThermoType>& chemistry_;
+        const chemistryModel<ThermoType>& chemistry_;
 
         //- Switch to select performance logging
         Switch log_;
@@ -214,7 +214,7 @@ public:
         ISAT
         (
             const dictionary& chemistryProperties,
-            chemistryModel<ThermoType>& chemistry
+            const chemistryModel<ThermoType>& chemistry
         );
 
         //- Disallow default bitwise copy construction
@@ -233,7 +233,7 @@ public:
             return true;
         }
 
-        chemistryModel<ThermoType>& chemistry()
+        const chemistryModel<ThermoType>& chemistry()
         {
             return chemistry_;
         }
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H b/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H
index fd6a30be05..ab982e6753 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H
@@ -276,7 +276,7 @@ public:
     // Member Functions
 
         //- Access to the chemistryModel
-        inline chemistryModel<ThermoType>& chemistry()
+        inline const chemistryModel<ThermoType>& chemistry()
         {
             return table_.chemistry();
         }
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.C b/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.C
index 1247f7e335..3e99a2ad8f 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.C
@@ -32,7 +32,7 @@ template<class ThermoType>
 Foam::chemistryTabulationMethod<ThermoType>::chemistryTabulationMethod
 (
     const dictionary& dict,
-    chemistryModel<ThermoType>& chemistry
+    const chemistryModel<ThermoType>& chemistry
 )
 {}
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.H b/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.H
index 8a64b950c8..bcea0c4bb6 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.H
@@ -71,7 +71,7 @@ public:
         dictionary,
         (
             const dictionary& dict,
-            chemistryModel<ThermoType>& chemistry
+            const chemistryModel<ThermoType>& chemistry
         ),
         (dict, chemistry)
     );
@@ -83,7 +83,7 @@ public:
         chemistryTabulationMethod
         (
             const dictionary& dict,
-            chemistryModel<ThermoType>& chemistry
+            const chemistryModel<ThermoType>& chemistry
         );
 
 
@@ -92,7 +92,7 @@ public:
         static autoPtr<chemistryTabulationMethod> New
         (
             const IOdictionary& dict,
-            chemistryModel<ThermoType>& chemistry
+            const chemistryModel<ThermoType>& chemistry
         );
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethodNew.C b/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethodNew.C
index 24b69385c8..53f6c115fb 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethodNew.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethodNew.C
@@ -34,7 +34,7 @@ Foam::autoPtr<Foam::chemistryTabulationMethod<ThermoType>>
 Foam::chemistryTabulationMethod<ThermoType>::New
 (
     const IOdictionary& dict,
-    chemistryModel<ThermoType>& chemistry
+    const chemistryModel<ThermoType>& chemistry
 )
 {
     if (dict.found("tabulation"))
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.C b/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.C
index 712082a966..9909479276 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.C
@@ -31,7 +31,7 @@ template<class ThermoType>
 Foam::chemistryTabulationMethods::none<ThermoType>::none
 (
     const dictionary& chemistryProperties,
-    chemistryModel<ThermoType>& chemistry
+    const chemistryModel<ThermoType>& chemistry
 )
 :
     chemistryTabulationMethod<ThermoType>
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.H b/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.H
index 17c1475cce..5320a28563 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.H
@@ -65,7 +65,7 @@ public:
         none
         (
             const dictionary& chemistryProperties,
-            chemistryModel<ThermoType>& chemistry
+            const chemistryModel<ThermoType>& chemistry
         );
 
         //- Disallow default bitwise copy construction
diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C
index 4cfca81803..0e3bbc5110 100644
--- a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C
+++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C
@@ -77,8 +77,8 @@ void Foam::ode<ChemistryModel>::solve
     if (debug)
     {
         scalarField dcTp(this->nEqns(), rootSmall);
-        dcTp[nSpecie] = 300*rootSmall;
-        dcTp[nSpecie+1] = 1e5*rootSmall;
+        dcTp[nSpecie] = T*rootSmall;
+        dcTp[nSpecie+1] = p*rootSmall;
         this->check(0, cTp_, dcTp, li);
     }
 
diff --git a/src/thermophysicalModels/chemistryModel/reactions/fluxLimitedLangmuirHinshelwood/fluxLimitedLangmuirHinshelwoodReactionRate.H b/src/thermophysicalModels/chemistryModel/reactions/fluxLimitedLangmuirHinshelwood/fluxLimitedLangmuirHinshelwoodReactionRate.H
index a8eabfcf3e..271f140781 100644
--- a/src/thermophysicalModels/chemistryModel/reactions/fluxLimitedLangmuirHinshelwood/fluxLimitedLangmuirHinshelwoodReactionRate.H
+++ b/src/thermophysicalModels/chemistryModel/reactions/fluxLimitedLangmuirHinshelwood/fluxLimitedLangmuirHinshelwoodReactionRate.H
@@ -149,28 +149,15 @@ public:
             const label li
         ) const;
 
-        //- Third-body efficiencies (beta = 1-alpha)
-        //  non-empty only for third-body reactions
-        //  with enhanced molecularity (alpha != 1)
-        inline const List<Tuple2<label, scalar>>& beta() const;
+        inline bool hasDdc() const;
 
-        //- Species concentration derivative of the pressure dependent term
-        inline void dcidc
+        inline void ddc
         (
             const scalar p,
             const scalar T,
             const scalarField& c,
             const label li,
-            scalarField& dcidc
-        ) const;
-
-        //- Temperature derivative of the pressure dependent term
-        inline scalar dcidT
-        (
-            const scalar p,
-            const scalar T,
-            const scalarField& c,
-            const label li
+            scalarField& ddc
         ) const;
 
         //- Write to stream
diff --git a/src/thermophysicalModels/chemistryModel/reactions/fluxLimitedLangmuirHinshelwood/fluxLimitedLangmuirHinshelwoodReactionRateI.H b/src/thermophysicalModels/chemistryModel/reactions/fluxLimitedLangmuirHinshelwood/fluxLimitedLangmuirHinshelwoodReactionRateI.H
index ed8aabe214..d8b5795b02 100644
--- a/src/thermophysicalModels/chemistryModel/reactions/fluxLimitedLangmuirHinshelwood/fluxLimitedLangmuirHinshelwoodReactionRateI.H
+++ b/src/thermophysicalModels/chemistryModel/reactions/fluxLimitedLangmuirHinshelwood/fluxLimitedLangmuirHinshelwoodReactionRateI.H
@@ -339,33 +339,22 @@ inline Foam::scalar Foam::fluxLimitedLangmuirHinshelwoodReactionRate::ddT
 }
 
 
-inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
-Foam::fluxLimitedLangmuirHinshelwoodReactionRate::beta() const
+inline bool Foam::fluxLimitedLangmuirHinshelwoodReactionRate::hasDdc() const
 {
-    return NullObjectRef<List<Tuple2<label, scalar>>>();
+    return false;
 }
 
 
-inline void Foam::fluxLimitedLangmuirHinshelwoodReactionRate::dcidc
+inline void Foam::fluxLimitedLangmuirHinshelwoodReactionRate::ddc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
     const label li,
-    scalarField& dcidc
-) const
-{}
-
-
-inline Foam::scalar Foam::fluxLimitedLangmuirHinshelwoodReactionRate::dcidT
-(
-    const scalar p,
-    const scalar T,
-    const scalarField& c,
-    const label li
+    scalarField& ddc
 ) const
 {
-    return 0;
+    ddc = 0;
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C b/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C
index 17024bc4f1..2e33ac2614 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C
@@ -169,37 +169,39 @@ Foam::scalar Foam::IrreversibleReaction<ReactionThermo, ReactionRate>::dkrdT
 
 
 template<class ReactionThermo, class ReactionRate>
-const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
-Foam::IrreversibleReaction<ReactionThermo, ReactionRate>::beta() const
+bool Foam::IrreversibleReaction<ReactionThermo, ReactionRate>::hasDkdc() const
 {
-    return k_.beta();
+    return k_.hasDdc();
 }
 
 
 template<class ReactionThermo, class ReactionRate>
-void Foam::IrreversibleReaction<ReactionThermo, ReactionRate>::dcidc
+void Foam::IrreversibleReaction<ReactionThermo, ReactionRate>::dkfdc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
     const label li,
-    scalarField& dcidc
+    scalarField& dkfdc
 ) const
 {
-    k_.dcidc(p, T, c, li, dcidc);
+    k_.ddc(p, T, c, li, dkfdc);
 }
 
 
 template<class ReactionThermo, class ReactionRate>
-Foam::scalar Foam::IrreversibleReaction<ReactionThermo, ReactionRate>::dcidT
+void Foam::IrreversibleReaction<ReactionThermo, ReactionRate>::dkrdc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
-    const label li
+    const label li,
+    const scalarField& dkfdc,
+    const scalar kr,
+    scalarField& dkrdc
 ) const
 {
-    return k_.dcidT(p, T, c, li);
+    dkrdc = 0;
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H b/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H
index 800b341548..168a977011 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H
@@ -181,8 +181,7 @@ public:
                 const label li
             ) const;
 
-            //- Temperature derivative of reverse rate
-            //  Returns 0
+            //- Temperature derivative of reverse rate. Returns zero.
             virtual scalar dkrdT
             (
                 const scalar p,
@@ -193,35 +192,34 @@ public:
                 const scalar kr
             ) const;
 
-            //- Third-body efficiencies (beta = 1-alpha)
-            //  non-empty only for third-body reactions
-            //  with enhanced molecularity (alpha != 1)
-            virtual const List<Tuple2<label, scalar>>& beta() const;
+            //- Does this reaction have concentration-dependent rate constants?
+            virtual bool hasDkdc() const;
 
-            //- Species concentration derivative of the pressure dependent term
-            //  By default this value is 1 as it multiplies the third-body term
-            virtual void dcidc
+            //- Concentration derivative of forward rate
+            void dkfdc
             (
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
                 const label li,
-                scalarField& dcidc
+                scalarField& dkfdc
             ) const;
 
-            //- Temperature derivative of the pressure dependent term
-            //  By default this value is 0 since ddT of molecularity is approx.0
-            virtual scalar dcidT
+            //- Concentration derivative of reverse rate. Returns zero.
+            void dkrdc
             (
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
-                const label li
+                const label li,
+                const scalarField& dkfdc,
+                const scalar kr,
+                scalarField& dkrdc
             ) const;
 
 
         //- Write
-        virtual void write(Ostream&) const;    // Private Member Functions
+        virtual void write(Ostream&) const;
 
 
     // Member Operators
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C
index 3656063b92..81207a2452 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C
@@ -181,40 +181,40 @@ Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::dkrdT
 
 
 template<class ReactionThermo, class ReactionRate>
-const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
-Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::
-beta() const
+bool Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::
+hasDkdc() const
 {
-    return fk_.beta();
+    return fk_.hasDdc() || rk_.hasDdc();
 }
 
 
 template<class ReactionThermo, class ReactionRate>
-void
-Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::dcidc
+void Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::dkfdc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
     const label li,
-    scalarField& dcidc
+    scalarField& dkfdc
 ) const
 {
-    fk_.dcidc(p, T, c, li, dcidc);
+    fk_.ddc(p, T, c, li, dkfdc);
 }
 
 
 template<class ReactionThermo, class ReactionRate>
-Foam::scalar
-Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::dcidT
+void Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::dkrdc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
-    const label li
+    const label li,
+    const scalarField& dkfdc,
+    const scalar kr,
+    scalarField& dkrdc
 ) const
 {
-    return fk_.dcidT(p, T, c, li);
+    rk_.ddc(p, T, c, li, dkrdc);
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H
index 946d85dade..15dbbf9b3c 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H
@@ -163,7 +163,7 @@ public:
                 const label li
             ) const;
 
-            //- Reverse rate constant from the given formard rate constant
+            //- Reverse rate constant from the given forward rate constant
             virtual scalar kr
             (
                 const scalar kfwd,
@@ -207,30 +207,29 @@ public:
                 const scalar kr
             ) const;
 
-            //- Third-body efficiencies (beta = 1-alpha)
-            //  non-empty only for third-body reactions
-            //  with enhanced molecularity (alpha != 1)
-            virtual const List<Tuple2<label, scalar>>& beta() const;
+            //- Does this reaction have concentration-dependent rate constants?
+            virtual bool hasDkdc() const;
 
-            //- Species concentration derivative of the pressure dependent term
-            //  By default this value is 1 as it multiplies the third-body term
-            virtual void dcidc
+            //- Concentration derivative of forward rate
+            void dkfdc
             (
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
                 const label li,
-                scalarField& dcidc
+                scalarField& dkfdc
             ) const;
 
-            //- Temperature derivative of the pressure dependent term
-            //  By default this value is 0 since ddT of molecularity is approx.0
-            virtual scalar dcidT
+            //- Concentration derivative of reverse rate
+            void dkrdc
             (
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
-                const label li
+                const label li,
+                const scalarField& dkfdc,
+                const scalar kr,
+                scalarField& dkrdc
             ) const;
 
 
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
index bcf0b373cd..c861d30887 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
@@ -286,30 +286,30 @@ void Foam::Reaction<ReactionThermo>::write(Ostream& os) const
 
 
 template<class ReactionThermo>
-void Foam::Reaction<ReactionThermo>::omega
+void Foam::Reaction<ReactionThermo>::C
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
     const label li,
-    scalarField& dcdt
+    scalar& Cf,
+    scalar& Cr
 ) const
 {
-    scalar omegaf, omegar;
-
-    const scalar omegaI = omega(p, T, c, li, omegaf, omegar);
+    Cf = Cr = 1;
 
     forAll(lhs(), i)
     {
         const label si = lhs()[i].index;
-        const scalar sl = lhs()[i].stoichCoeff;
-        dcdt[si] -= sl*omegaI;
+        const specieExponent& el = lhs()[i].exponent;
+        Cf *= c[si] >= small || el >= 1 ? pow(max(c[si], 0), el) : 0;
     }
+
     forAll(rhs(), i)
     {
         const label si = rhs()[i].index;
-        const scalar sr = rhs()[i].stoichCoeff;
-        dcdt[si] += sr*omegaI;
+        const specieExponent& er = rhs()[i].exponent;
+        Cr *= c[si] >= small || er >= 1 ? pow(max(c[si], 0), er) : 0;
     }
 }
 
@@ -325,225 +325,229 @@ Foam::scalar Foam::Reaction<ReactionThermo>::omega
     scalar& omegar
 ) const
 {
-    scalar clippedT = min(max(T, this->Tlow()), this->Thigh());
+    const scalar clippedT = min(max(T, this->Tlow()), this->Thigh());
 
-    omegaf = this->kf(p, clippedT, c, li);
-    omegar = this->kr(omegaf, p, clippedT, c, li);
+    // Rate constants
+    const scalar kf = this->kf(p, clippedT, c, li);
+    const scalar kr = this->kr(kf, p, clippedT, c, li);
 
-    forAll(lhs(), i)
-    {
-        const label si = lhs()[i].index;
-        const scalar el = lhs()[i].exponent;
-        omegaf *= c[si] >= small || el >= 1 ? pow(max(c[si], 0), el) : 0;
-    }
-    forAll(rhs(), i)
-    {
-        const label si = rhs()[i].index;
-        const scalar er = rhs()[i].exponent;
-        omegar *= c[si] >= small || er >= 1 ? pow(max(c[si], 0), er) : 0;
-    }
+    // Concentration products
+    scalar Cf, Cr;
+    this->C(p, T, c, li, Cf, Cr);
+
+    omegaf = kf*Cf;
+    omegar = kr*Cr;
 
     return omegaf - omegar;
 }
 
 
 template<class ReactionThermo>
-void Foam::Reaction<ReactionThermo>::dwdc
+void Foam::Reaction<ReactionThermo>::dNdtByV
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
     const label li,
-    scalarSquareMatrix& J,
-    scalarField& dcdt,
-    scalar& omegaI,
-    scalar& kfwd,
-    scalar& kbwd,
+    scalarField& dNdtByV,
     const bool reduced,
-    const List<label>& c2s
+    const List<label>& c2s,
+    const label Nsi0
 ) const
 {
     scalar omegaf, omegar;
-
-    omegaI = omega(p, T, c, li, omegaf, omegar);
+    const scalar omega = this->omega(p, T, c, li, omegaf, omegar);
 
     forAll(lhs(), i)
     {
         const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
         const scalar sl = lhs()[i].stoichCoeff;
-        dcdt[si] -= sl*omegaI;
+        dNdtByV[Nsi0 + si] -= sl*omega;
     }
     forAll(rhs(), i)
     {
         const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
         const scalar sr = rhs()[i].stoichCoeff;
-        dcdt[si] += sr*omegaI;
-    }
-
-    kfwd = this->kf(p, T, c, li);
-    kbwd = this->kr(kfwd, p, T, c, li);
-
-    forAll(lhs(), j)
-    {
-        const label sj = reduced ? c2s[lhs()[j].index] : lhs()[j].index;
-        scalar kf = kfwd;
-        forAll(lhs(), i)
-        {
-            const label si = lhs()[i].index;
-            const scalar el = lhs()[i].exponent;
-            if (i == j)
-            {
-                kf *=
-                    c[si] >= small || el >= 1
-                  ? el*pow(max(c[si], 0), el - 1)
-                  : 0;
-            }
-            else
-            {
-                kf *=
-                    c[si] >= small || el >= 1
-                  ? pow(max(c[si], 0), el)
-                  : 0;
-            }
-        }
-
-        forAll(lhs(), i)
-        {
-            const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
-            const scalar sl = lhs()[i].stoichCoeff;
-            J(si, sj) -= sl*kf;
-        }
-        forAll(rhs(), i)
-        {
-            const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
-            const scalar sr = rhs()[i].stoichCoeff;
-            J(si, sj) += sr*kf;
-        }
-    }
-
-    forAll(rhs(), j)
-    {
-        const label sj = reduced ? c2s[rhs()[j].index] : rhs()[j].index;
-        scalar kr = kbwd;
-        forAll(rhs(), i)
-        {
-            const label si = rhs()[i].index;
-            const scalar er = rhs()[i].exponent;
-            if (i == j)
-            {
-                kr *=
-                    c[si] >= small || er >= 1
-                  ? er*pow(max(c[si], 0), er - 1)
-                  : 0;
-            }
-            else
-            {
-                kr *=
-                    c[si] >= small || er >= 1
-                  ? pow(max(c[si], 0), er)
-                  : 0;
-            }
-        }
-
-        forAll(lhs(), i)
-        {
-            const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
-            const scalar sl = lhs()[i].stoichCoeff;
-            J(si, sj) += sl*kr;
-        }
-        forAll(rhs(), i)
-        {
-            const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
-            const scalar sr = rhs()[i].stoichCoeff;
-            J(si, sj) -= sr*kr;
-        }
-    }
-
-    // When third-body species are involved, additional terms are added
-    // beta function returns an empty list when third-body are not involved
-    const List<Tuple2<label, scalar>>& beta = this->beta();
-    if (notNull(beta))
-    {
-        // This temporary array needs to be cached for efficiency
-        scalarField dcidc(beta.size());
-        this->dcidc(p, T, c, li, dcidc);
-
-        forAll(beta, j)
-        {
-            label sj = beta[j].first();
-            sj = reduced ? c2s[sj] : sj;
-            if (sj != -1)
-            {
-                forAll(lhs(), i)
-                {
-                    const label si =
-                        reduced ? c2s[lhs()[i].index] : lhs()[i].index;
-                    const scalar sl = lhs()[i].stoichCoeff;
-                    J(si, sj) -= sl*dcidc[j]*omegaI;
-                }
-                forAll(rhs(), i)
-                {
-                    const label si =
-                        reduced ? c2s[rhs()[i].index] : rhs()[i].index;
-                    const scalar sr = rhs()[i].stoichCoeff;
-                    J(si, sj) += sr*dcidc[j]*omegaI;
-                }
-            }
-        }
+        dNdtByV[Nsi0 + si] += sr*omega;
     }
 }
 
 
 template<class ReactionThermo>
-void Foam::Reaction<ReactionThermo>::dwdT
+void Foam::Reaction<ReactionThermo>::ddNdtByVdcTp
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
     const label li,
-    const scalar omegaI,
-    const scalar kfwd,
-    const scalar kbwd,
-    scalarSquareMatrix& J,
+    scalarField& dNdtByV,
+    scalarSquareMatrix& ddNdtByVdcTp,
     const bool reduced,
     const List<label>& c2s,
-    const label indexT
+    const label Nsi0,
+    const label Tsi,
+    scalarField& cTpWork0,
+    scalarField& cTpWork1
 ) const
 {
-    scalar dkfdT = this->dkfdT(p, T, c, li);
-    scalar dkrdT = this->dkrdT(p, T, c, li, dkfdT, kbwd);
+    // Rate constants
+    const scalar kf = this->kf(p, T, c, li);
+    const scalar kr = this->kr(kf, p, T, c, li);
 
-    forAll(lhs(), i)
-    {
-        const label si = lhs()[i].index;
-        const scalar el = lhs()[i].exponent;
-        dkfdT *= c[si] >= small || el >= 1 ? pow(max(c[si], 0), el) : 0;
-    }
-    forAll(rhs(), i)
-    {
-        const label si = rhs()[i].index;
-        const scalar er = rhs()[i].exponent;
-        dkrdT *= c[si] >= small || er >= 1 ? pow(max(c[si], 0), er) : 0;
-    }
+    // Concentration products
+    scalar Cf, Cr;
+    this->C(p, T, c, li, Cf, Cr);
 
-    const scalar dqidT = dkfdT - dkrdT;
+    // Overall reaction rate
+    const scalar omega = kf*Cf - kr*Cr;
 
-    // For reactions including third-body efficiencies or pressure dependent
-    // reaction, an additional term is needed
-    const scalar dcidT = omegaI*this->dcidT(p, T, c, li);
-
-    // J(i, indexT) = sum_reactions nu_i dqdT
+    // Specie reaction rates
     forAll(lhs(), i)
     {
         const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
         const scalar sl = lhs()[i].stoichCoeff;
-        J(si, indexT) -= sl*(dqidT + dcidT);
+        dNdtByV[Nsi0 + si] -= sl*omega;
     }
     forAll(rhs(), i)
     {
         const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
         const scalar sr = rhs()[i].stoichCoeff;
-        J(si, indexT) += sr*(dqidT + dcidT);
+        dNdtByV[Nsi0 + si] += sr*omega;
+    }
+
+    // Jacobian contributions from the derivative of the concentration products
+    // w.r.t. concentration
+    {
+        forAll(lhs(), j)
+        {
+            const label sj = reduced ? c2s[lhs()[j].index] : lhs()[j].index;
+
+            scalar dCfdcj = 1;
+            forAll(lhs(), i)
+            {
+                const label si = lhs()[i].index;
+                const specieExponent& el = lhs()[i].exponent;
+                if (i == j)
+                {
+                    dCfdcj *=
+                        c[si] >= small || el >= 1
+                      ? el*pow(max(c[si], 0), el - specieExponent(1))
+                      : 0;
+                }
+                else
+                {
+                    dCfdcj *=
+                        c[si] >= small || el >= 1
+                      ? pow(max(c[si], 0), el)
+                      : 0;
+                }
+            }
+
+            forAll(lhs(), i)
+            {
+                const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
+                const scalar sl = lhs()[i].stoichCoeff;
+                ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) -= sl*kf*dCfdcj;
+            }
+            forAll(rhs(), i)
+            {
+                const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
+                const scalar sr = rhs()[i].stoichCoeff;
+                ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) += sr*kf*dCfdcj;
+            }
+        }
+
+        forAll(rhs(), j)
+        {
+            const label sj = reduced ? c2s[rhs()[j].index] : rhs()[j].index;
+
+            scalar dCrcj = 1;
+            forAll(rhs(), i)
+            {
+                const label si = rhs()[i].index;
+                const specieExponent& er = rhs()[i].exponent;
+                if (i == j)
+                {
+                    dCrcj *=
+                        c[si] >= small || er >= 1
+                      ? er*pow(max(c[si], 0), er - specieExponent(1))
+                      : 0;
+                }
+                else
+                {
+                    dCrcj *=
+                        c[si] >= small || er >= 1
+                      ? pow(max(c[si], 0), er)
+                      : 0;
+                }
+            }
+
+            forAll(lhs(), i)
+            {
+                const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
+                const scalar sl = lhs()[i].stoichCoeff;
+                ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) += sl*kr*dCrcj;
+            }
+            forAll(rhs(), i)
+            {
+                const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
+                const scalar sr = rhs()[i].stoichCoeff;
+                ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) -= sr*kr*dCrcj;
+            }
+        }
+    }
+
+    // Jacobian contributions from the derivative of the rate constants
+    // w.r.t. temperature
+    {
+        const scalar dkfdT = this->dkfdT(p, T, c, li);
+        const scalar dkrdT = this->dkrdT(p, T, c, li, dkfdT, kr);
+
+        const scalar dwdT = dkfdT*Cf - dkrdT*Cr;
+        forAll(lhs(), i)
+        {
+            const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
+            const scalar sl = lhs()[i].stoichCoeff;
+            ddNdtByVdcTp(Nsi0 + si, Tsi) -= sl*dwdT;
+        }
+        forAll(rhs(), i)
+        {
+            const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
+            const scalar sr = rhs()[i].stoichCoeff;
+            ddNdtByVdcTp(Nsi0 + si, Tsi) += sr*dwdT;
+        }
+    }
+
+    // Jacobian contributions from the derivative of the rate constants
+    // w.r.t. concentration
+    if (hasDkdc())
+    {
+        scalarField& dkfdc = cTpWork0;
+        scalarField& dkrdc = cTpWork1;
+
+        this->dkfdc(p, T, c, li, dkfdc);
+        this->dkrdc(p, T, c, li, dkfdc, kr, dkrdc);
+
+        forAll(c, j)
+        {
+            const label sj = reduced ? c2s[j] : j;
+
+            if (sj == -1) continue;
+
+            const scalar dwdc = dkfdc[j]*Cf - dkrdc[j]*Cr;
+            forAll(lhs(), i)
+            {
+                const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
+                const scalar sl = lhs()[i].stoichCoeff;
+                ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) -= sl*dwdc;
+            }
+            forAll(rhs(), i)
+            {
+                const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
+                const scalar sr = rhs()[i].stoichCoeff;
+                ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) += sr*dwdc;
+            }
+        }
     }
 }
 
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
index 734b55b18e..11311075bf 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
@@ -224,14 +224,15 @@ public:
 
         // Reaction rate coefficients
 
-            //- Net reaction rate for individual species
-            void omega
+            //- Concentration powers
+            void C
             (
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
                 const label li,
-                scalarField& dcdt
+                scalar& Cf,
+                scalar& Cr
             ) const;
 
             //- Net reaction rate
@@ -245,6 +246,20 @@ public:
                 scalar& omegar
             ) const;
 
+            //- The net reaction rate for each species involved
+            void dNdtByV
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const label li,
+                scalarField& dNdtByV,
+                const bool reduced,
+                const List<label>& c2s,
+                const label Nsi0
+            ) const;
+
+
         // Reaction rate coefficients
 
             //- Forward rate constant
@@ -259,7 +274,7 @@ public:
             //- Reverse rate constant from the given forward rate constant
             virtual scalar kr
             (
-                const scalar kfwd,
+                const scalar kf,
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
@@ -278,40 +293,6 @@ public:
 
         // Jacobian coefficients
 
-            //- Derivative of the net reaction rate for each species involved
-            //  w.r.t. the species concentration
-            void dwdc
-            (
-                const scalar p,
-                const scalar T,
-                const scalarField& c,
-                const label li,
-                scalarSquareMatrix& J,
-                scalarField& dcdt,
-                scalar& omegaI,
-                scalar& kfwd,
-                scalar& kbwd,
-                const bool reduced,
-                const List<label>& c2s
-            ) const;
-
-            //- Derivative of the net reaction rate for each species involved
-            //  w.r.t. the temperature
-            void dwdT
-            (
-                const scalar p,
-                const scalar T,
-                const scalarField& c,
-                const label li,
-                const scalar omegaI,
-                const scalar kfwd,
-                const scalar kbwd,
-                scalarSquareMatrix& J,
-                const bool reduced,
-                const List<label>& c2s,
-                const label indexT
-            ) const;
-
             //- Temperature derivative of forward rate
             virtual scalar dkfdT
             (
@@ -332,30 +313,49 @@ public:
                 const scalar kr
             ) const = 0;
 
-            //- Third-body efficiencies (beta = 1-alpha)
-            //  non-empty only for third-body reactions
-            //  with enhanced molecularity (alpha != 1)
-            virtual const List<Tuple2<label, scalar>>& beta() const = 0;
+            //- Does this reaction have concentration-dependent rate constants?
+            virtual bool hasDkdc() const = 0;
 
-            //- Species concentration derivative of the pressure dependent term
-            virtual void dcidc
+            //- Concentration derivative of forward rate
+            virtual void dkfdc
             (
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
                 const label li,
-                scalarField& dcidc
+                scalarField& dkfdc
             ) const = 0;
 
-            //- Temperature derivative of the pressure dependent term
-            virtual scalar dcidT
+            //- Concentration derivative of reverse rate
+            virtual void dkrdc
             (
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
-                const label li
+                const label li,
+                const scalarField& dkfdc,
+                const scalar kr,
+                scalarField& dkrdc
             ) const = 0;
 
+            //- Derivative of the net reaction rate for each species involved
+            //  w.r.t. the concentration and temperature
+            void ddNdtByVdcTp
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const label li,
+                scalarField& dNdtByV,
+                scalarSquareMatrix& ddNdtByVdcTp,
+                const bool reduced,
+                const List<label>& c2s,
+                const label csi0,
+                const label Tsi,
+                scalarField& cTpWork0,
+                scalarField& cTpWork1
+            ) const;
+
 
         //- Write
         virtual void write(Ostream&) const;
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C b/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C
index 7610898bc3..f603f3fa13 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C
@@ -185,22 +185,21 @@ Foam::scalar Foam::ReactionProxy<ReactionThermo>::dkrdT
 
 
 template<class ReactionThermo>
-const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
-Foam::ReactionProxy<ReactionThermo>::beta() const
+bool Foam::ReactionProxy<ReactionThermo>::hasDkdc() const
 {
     NotImplemented;
-    return NullObjectRef<List<Tuple2<label, scalar>>>();
+    return false;
 }
 
 
 template<class ReactionThermo>
-void Foam::ReactionProxy<ReactionThermo>::dcidc
+void Foam::ReactionProxy<ReactionThermo>::dkfdc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
     const label li,
-    scalarField& dcidc
+    scalarField& dkfdc
 ) const
 {
     NotImplemented;
@@ -208,16 +207,18 @@ void Foam::ReactionProxy<ReactionThermo>::dcidc
 
 
 template<class ReactionThermo>
-Foam::scalar Foam::ReactionProxy<ReactionThermo>::dcidT
+void Foam::ReactionProxy<ReactionThermo>::dkrdc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
-    const label li
+    const label li,
+    const scalarField& dkfdc,
+    const scalar kr,
+    scalarField& dkrdc
 ) const
 {
     NotImplemented;
-    return 0;
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H b/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H
index 477f7482bf..5ea4e8ed9d 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H
@@ -162,28 +162,29 @@ public:
                 const scalar kr
             ) const;
 
-            //- Third-body efficiencies (beta = 1-alpha)
-            //  non-empty only for third-body reactions
-            //  with enhanced molecularity (alpha != 1)
-            virtual const List<Tuple2<label, scalar>>& beta() const;
+            //- Does this reaction have concentration-dependent rate constants?
+            virtual bool hasDkdc() const;
 
-            //- Species concentration derivative of the pressure dependent term
-            virtual void dcidc
+            //- Concentration derivative of forward rate
+            void dkfdc
             (
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
                 const label li,
-                scalarField& dcidc
+                scalarField& dkfdc
             ) const;
 
-            //- Temperature derivative of the pressure dependent term
-            virtual scalar dcidT
+            //- Concentration derivative of reverse rate
+            void dkrdc
             (
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
-                const label li
+                const label li,
+                const scalarField& dkfdc,
+                const scalar kr,
+                scalarField& dkrdc
             ) const;
 };
 
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C
index f9b465accc..5b1363e531 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C
@@ -164,44 +164,48 @@ Foam::scalar Foam::ReversibleReaction<ReactionThermo, ReactionRate>::dkrdT
     const scalar kr
 ) const
 {
-    scalar Kc = max(this->Kc(p, T), rootSmall);
+    const scalar Kc = max(this->Kc(p, T), rootSmall);
 
-    return dkfdT/Kc - kr*this->dKcdTbyKc(p, T);
+    return dkfdT/Kc - (Kc > rootSmall ? kr*this->dKcdTbyKc(p, T) : 0);
 }
 
 
 template<class ReactionThermo, class ReactionRate>
-const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
-Foam::ReversibleReaction<ReactionThermo, ReactionRate>::beta() const
+bool Foam::ReversibleReaction<ReactionThermo, ReactionRate>::hasDkdc() const
 {
-    return k_.beta();
+    return k_.hasDdc();
 }
 
 
 template<class ReactionThermo, class ReactionRate>
-void Foam::ReversibleReaction<ReactionThermo, ReactionRate>::dcidc
+void Foam::ReversibleReaction<ReactionThermo, ReactionRate>::dkfdc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
     const label li,
-    scalarField& dcidc
+    scalarField& dkfdc
 ) const
 {
-    k_.dcidc(p, T, c, li, dcidc);
+    k_.ddc(p, T, c, li, dkfdc);
 }
 
 
 template<class ReactionThermo, class ReactionRate>
-Foam::scalar Foam::ReversibleReaction<ReactionThermo, ReactionRate>::dcidT
+void Foam::ReversibleReaction<ReactionThermo, ReactionRate>::dkrdc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
-    const label li
+    const label li,
+    const scalarField& dkfdc,
+    const scalar kr,
+    scalarField& dkrdc
 ) const
 {
-    return k_.dcidT(p, T, c, li);
+    const scalar Kc = max(this->Kc(p, T), rootSmall);
+
+    dkrdc = dkfdc/Kc;
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H
index 517f3868c0..dded56aea2 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H
@@ -149,7 +149,7 @@ public:
                 const label li
             ) const;
 
-            //- Reverse rate constant from the given formard rate constant
+            //- Reverse rate constant from the given forward rate constant
             virtual scalar kr
             (
                 const scalar kfwd,
@@ -193,30 +193,29 @@ public:
                 const scalar kr
             ) const;
 
-            //- Third-body efficiencies (beta = 1-alpha)
-            //  non-empty only for third-body reactions
-            //  with enhanced molecularity (alpha != 1)
-            virtual const List<Tuple2<label, scalar>>& beta() const;
+            //- Does this reaction have concentration-dependent rate constants?
+            virtual bool hasDkdc() const;
 
-            //- Species concentration derivative of the pressure dependent term
-            //  By default this value is 1 as it multiplies the third-body term
-            virtual void dcidc
+            //- Concentration derivative of forward rate
+            void dkfdc
             (
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
                 const label li,
-                scalarField& dcidc
+                scalarField& dkfdc
             ) const;
 
-            //- Temperature derivative of the pressure dependent term
-            //  By default this value is 0 since ddT of molecularity is approx.0
-            virtual scalar dcidT
+            //- Concentration derivative of reverse rate
+            void dkrdc
             (
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
-                const label li
+                const label li,
+                const scalarField& dkfdc,
+                const scalar kr,
+                scalarField& dkrdc
             ) const;
 
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H
index 4e64520d97..a4c6876984 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H
@@ -99,6 +99,7 @@ public:
         //- Post-evaluation hook
         inline void postEvaluate() const;
 
+        //- Return the rate
         inline scalar operator()
         (
             const scalar p,
@@ -107,6 +108,7 @@ public:
             const label li
         ) const;
 
+        //- The derivative of the rate w.r.t. temperature
         inline scalar ddT
         (
             const scalar p,
@@ -115,30 +117,17 @@ public:
             const label li
         ) const;
 
-        //- Third-body efficiencies (beta = 1-alpha)
-        //  non-empty only for third-body reactions
-        //  with enhanced molecularity (alpha != 1)
-        inline const List<Tuple2<label, scalar>>& beta() const;
+        //- Is the rate a function of concentration?
+        inline bool hasDdc() const;
 
-        //- Species concentration derivative of the pressure dependent term
-        //  By default this value is 1 as it multiplies the third-body term
-        inline void dcidc
+        //- The derivative of the rate w.r.t. concentration
+        inline void ddc
         (
             const scalar p,
             const scalar T,
             const scalarField& c,
             const label li,
-            scalarField& dcidc
-        ) const;
-
-        //- Temperature derivative of the pressure dependent term
-        //  By default this value is 0 since ddT of molecularity is approx.0
-        inline scalar dcidT
-        (
-            const scalar p,
-            const scalar T,
-            const scalarField& c,
-            const label li
+            scalarField& ddc
         ) const;
 
         //- Write to stream
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H
index 850f0cf177..0133774e3b 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H
@@ -108,33 +108,22 @@ inline Foam::scalar Foam::ArrheniusReactionRate::ddT
 }
 
 
-inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
-Foam::ArrheniusReactionRate::beta() const
+inline bool Foam::ArrheniusReactionRate::hasDdc() const
 {
-    return NullObjectRef<List<Tuple2<label, scalar>>>();
+    return false;
 }
 
 
-inline void Foam::ArrheniusReactionRate::dcidc
+inline void Foam::ArrheniusReactionRate::ddc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
     const label li,
-    scalarField& dcidc
-) const
-{}
-
-
-inline Foam::scalar Foam::ArrheniusReactionRate::dcidT
-(
-    const scalar p,
-    const scalar T,
-    const scalarField& c,
-    const label li
+    scalarField& ddc
 ) const
 {
-    return 0;
+    ddc = 0;
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H
index adc2cdb450..fbbb467cb4 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H
@@ -69,7 +69,6 @@ class ChemicallyActivatedReactionRate
         ReactionRate kInf_;
         ChemicallyActivationFunction F_;
         thirdBodyEfficiencies thirdBodyEfficiencies_;
-        List<Tuple2<label, scalar>> beta_;
 
 
 public:
@@ -109,6 +108,7 @@ public:
         //- Post-evaluation hook
         inline void postEvaluate() const;
 
+        //- Return the rate
         inline scalar operator()
         (
             const scalar p,
@@ -117,6 +117,7 @@ public:
             const label li
         ) const;
 
+        //- The derivative of the rate w.r.t. temperature
         inline scalar ddT
         (
             const scalar p,
@@ -125,33 +126,17 @@ public:
             const label li
         ) const;
 
-        //- Third-body efficiencies (beta = 1-alpha)
-        //  non-empty only for third-body reactions
-        //  with enhanced molecularity (alpha != 1)
-        inline const List<Tuple2<label, scalar>>& beta() const
-        {
-            return beta_;
-        }
+        //- Is the rate a function of concentration?
+        inline bool hasDdc() const;
 
-        //- Species concentration derivative of the pressure dependent term
-        //  By default this value is 1 as it multiplies the third-body term
-        inline void dcidc
+        //- The derivative of the rate w.r.t. concentration
+        inline void ddc
         (
             const scalar p,
             const scalar T,
             const scalarField& c,
             const label li,
-            scalarField& dcidc
-        ) const;
-
-        //- Temperature derivative of the pressure dependent term
-        //  By default this value is 0 since ddT of molecularity is approx.0
-        inline scalar dcidT
-        (
-            const scalar p,
-            const scalar T,
-            const scalarField& c,
-            const label li
+            scalarField& ddc
         ) const;
 
         //- Write to stream
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H
index 8d104b8080..44a9e89308 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H
@@ -23,6 +23,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#define CHEMICALLY_ACTIVATION_FUNCTION_JACOBIAN false
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class ReactionRate, class ChemicallyActivationFunction>
@@ -42,12 +44,7 @@ inline Foam::ChemicallyActivatedReactionRate
     kInf_(kInf),
     F_(F),
     thirdBodyEfficiencies_(tbes)
-{
-    forAll(tbes, i)
-    {
-        beta_.append(Tuple2<label, scalar>(i, tbes[i]));
-    }
-}
+{}
 
 
 template<class ReactionRate, class ChemicallyActivationFunction>
@@ -61,23 +58,11 @@ inline Foam::ChemicallyActivatedReactionRate
     const dictionary& dict
 )
 :
-    k0_(species, dict),
-    kInf_(species, dict),
-    F_(dict),
-    thirdBodyEfficiencies_(species, dict)
-{
-    forAll(thirdBodyEfficiencies_, i)
-    {
-        beta_.append
-        (
-            Tuple2<label, scalar>
-            (
-                i,
-                thirdBodyEfficiencies_[i]
-            )
-        );
-    }
-}
+    k0_(species, dict.subDict("k0")),
+    kInf_(species, dict.subDict("kInf")),
+    F_(dict.subDict("F")),
+    thirdBodyEfficiencies_(species, dict.subDict("thirdBodyEfficiencies"))
+{}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -121,9 +106,11 @@ inline Foam::scalar Foam::ChemicallyActivatedReactionRate
 {
     const scalar k0 = k0_(p, T, c, li);
     const scalar kInf = kInf_(p, T, c, li);
-    const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+    const scalar M = thirdBodyEfficiencies_.M(c);
+    const scalar Pr = k0/kInf*M;
+    const scalar F = F_(T, Pr);
 
-    return k0*(1/(1 + Pr))*F_(T, Pr);
+    return k0/(1 + Pr)*F;
 }
 
 
@@ -142,9 +129,39 @@ inline Foam::scalar Foam::ChemicallyActivatedReactionRate
 {
     const scalar k0 = k0_(p, T, c, li);
     const scalar kInf = kInf_(p, T, c, li);
-    const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+    const scalar M = thirdBodyEfficiencies_.M(c);
+    const scalar Pr = k0/kInf*M;
+    const scalar F = F_(T, Pr);
 
-    return (1/(1 + Pr))*F_(T, Pr)*k0_.ddT(p, T, c, li);
+    const scalar dk0dT = k0_.ddT(p, T, c, li);
+
+    #if CHEMICALLY_ACTIVATION_FUNCTION_JACOBIAN
+
+    const scalar dkInfdT = kInf_.ddT(p, T, c, li);
+    const scalar dPrdT = (M*dk0dT - Pr*dkInfdT)/kInf;
+    const scalar dFdT = F_.ddT(T, Pr, F) + F_.ddPr(T, Pr, F)*dPrdT;
+
+    return
+        dk0dT/(1 + Pr)*F
+      - k0*dPrdT/sqr(1 + Pr)*F
+      + k0/(1 + Pr)*dFdT;
+
+    #else
+
+    return dk0dT/(1 + Pr)*F;
+
+    #endif
+}
+
+
+template<class ReactionRate, class ChemicallyActivationFunction>
+inline bool Foam::ChemicallyActivatedReactionRate
+<
+    ReactionRate,
+    ChemicallyActivationFunction
+>::hasDdc() const
+{
+    return true;
 }
 
 
@@ -153,74 +170,44 @@ inline void Foam::ChemicallyActivatedReactionRate
 <
     ReactionRate,
     ChemicallyActivationFunction
->::dcidc
+>::ddc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
     const label li,
-    scalarField& dcidc
+    scalarField& ddc
 ) const
 {
+    const scalar k0 = k0_(p, T, c, li);
+    const scalar kInf = kInf_(p, T, c, li);
     const scalar M = thirdBodyEfficiencies_.M(c);
+    const scalar Pr = k0/kInf*M;
+    const scalar F = F_(T, Pr);
 
-    if (M > small)
-    {
-        const scalar k0 = k0_(p, T, c, li);
-        const scalar kInf = kInf_(p, T, c, li);
+    #if CHEMICALLY_ACTIVATION_FUNCTION_JACOBIAN
 
-        const scalar Pr = k0*M/kInf;
-        const scalar F = F_(T, Pr);
+    scalarField dk0dc(c.size(), 0);
+    k0_.ddc(p, T, c, li, dk0dc);
+    scalarField dkInfdc(c.size(), 0);
+    kInf_.ddc(p, T, c, li, dkInfdc);
+    tmp<scalarField> tdMdc = thirdBodyEfficiencies_.dMdc(c);
+    const scalarField& dMdc = tdMdc();
+    scalarField dPrdc((M*dk0dc - Pr*dkInfdc + k0*dMdc)/kInf);
+    const scalar dFdPr = F_.ddPr(T, Pr, F);
 
-        forAll(beta_, i)
-        {
-            const scalar dPrdci = -beta_[i].second()*k0/kInf;
-            const scalar dFdci = F_.ddc(Pr, F, dPrdci, T);
-            dcidc[i] = (-dPrdci/(1 + Pr) + dFdci/F);
-        }
-    }
-    else
-    {
-        forAll(beta_, i)
-        {
-            dcidc[i] = 0;
-        }
-    }
-}
+    ddc =
+        dk0dc/(1 + Pr)*F
+      - k0*dPrdc/sqr(1 + Pr)*F
+      + k0/(1 + Pr)*dFdPr*dPrdc;
 
+    #else
 
-template<class ReactionRate, class ChemicallyActivationFunction>
-inline Foam::scalar Foam::ChemicallyActivatedReactionRate
-<
-    ReactionRate,
-    ChemicallyActivationFunction
->::dcidT
-(
-    const scalar p,
-    const scalar T,
-    const scalarField& c,
-    const label li
-) const
-{
-    const scalar M = thirdBodyEfficiencies_.M(c);
+    k0_.ddc(p, T, c, li, ddc);
 
-    if (M > small)
-    {
-        const scalar k0 = k0_(p, T, c, li);
-        const scalar kInf = kInf_(p, T, c, li);
+    ddc *= 1/(1 + Pr)*F;
 
-        const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
-        const scalar F = F_(T, Pr);
-        const scalar dPrdT =
-            Pr*(k0_.ddT(p, T, c, li)/k0 - kInf_.ddT(p, T, c, li)/kInf - 1/T);
-        const scalar dFdT = F_.ddT(Pr, F, dPrdT, T);
-
-        return (-dPrdT/(1 + Pr) + dFdT/F);
-    }
-    else
-    {
-        return 0;
-    }
+    #endif
 }
 
 
@@ -231,10 +218,33 @@ inline void Foam::ChemicallyActivatedReactionRate
     ChemicallyActivationFunction
 >::write(Ostream& os) const
 {
+    os  << indent << "k0" << nl;
+    os  << indent << token::BEGIN_BLOCK << nl;
+    os  << incrIndent;
     k0_.write(os);
+    os  << decrIndent;
+    os  << indent << token::END_BLOCK << nl;
+
+    os  << indent << "kInf" << nl;
+    os  << indent << token::BEGIN_BLOCK << nl;
+    os  << incrIndent;
     kInf_.write(os);
+    os  << decrIndent;
+    os  << indent << token::END_BLOCK << nl;
+
+    os  << indent << "F" << nl;
+    os  << indent << token::BEGIN_BLOCK << nl;
+    os  << incrIndent;
     F_.write(os);
+    os  << decrIndent;
+    os  << indent << token::END_BLOCK << nl;
+
+    os  << indent << "thirdBodyEfficiencies" << nl;
+    os  << indent << token::BEGIN_BLOCK << nl;
+    os  << incrIndent;
     thirdBodyEfficiencies_.write(os);
+    os  << decrIndent;
+    os  << indent << token::END_BLOCK << nl;
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H
index 7534200cdb..7d7aafaf0b 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H
@@ -65,10 +65,12 @@ class FallOffReactionRate
     // Private Data
 
         ReactionRate k0_;
+
         ReactionRate kInf_;
+
         FallOffFunction F_;
+
         thirdBodyEfficiencies thirdBodyEfficiencies_;
-        List<Tuple2<label, scalar>> beta_;
 
 
 public:
@@ -106,6 +108,7 @@ public:
         //- Post-evaluation hook
         inline void postEvaluate() const;
 
+        //- Return the rate
         inline scalar operator()
         (
             const scalar p,
@@ -114,6 +117,7 @@ public:
             const label li
         ) const;
 
+        //- The derivative of the rate w.r.t. temperature
         inline scalar ddT
         (
             const scalar p,
@@ -122,31 +126,17 @@ public:
             const label li
         ) const;
 
-        //- Third-body efficiencies (beta = 1-alpha)
-        //  non-empty only for third-body reactions
-        //  with enhanced molecularity (alpha != 1)
-        inline const List<Tuple2<label, scalar>>& beta() const
-        {
-            return beta_;
-        }
+        //- Is the rate a function of concentration?
+        inline bool hasDdc() const;
 
-        //- Species concentration derivative of the pressure dependent term
-        inline void dcidc
+        //- The derivative of the rate w.r.t. concentration
+        inline void ddc
         (
             const scalar p,
             const scalar T,
             const scalarField& c,
             const label li,
-            scalarField& dcidc
-        ) const;
-
-        //- Temperature derivative of the pressure dependent term
-        inline scalar dcidT
-        (
-            const scalar p,
-            const scalar T,
-            const scalarField& c,
-            const label li
+            scalarField& ddc
         ) const;
 
         //- Write to stream
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H
index 2300527234..8863ba2c2c 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H
@@ -23,6 +23,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#define FALL_OFF_FUNCTION_JACOBIAN false
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class ReactionRate, class FallOffFunction>
@@ -39,12 +41,7 @@ FallOffReactionRate
     kInf_(kInf),
     F_(F),
     thirdBodyEfficiencies_(tbes)
-{
-    forAll(tbes, i)
-    {
-        beta_.append(Tuple2<label, scalar>(i, tbes[i]));
-    }
-}
+{}
 
 
 template<class ReactionRate, class FallOffFunction>
@@ -59,19 +56,7 @@ FallOffReactionRate
     kInf_(species, dict.subDict("kInf")),
     F_(dict.subDict("F")),
     thirdBodyEfficiencies_(species, dict.subDict("thirdBodyEfficiencies"))
-{
-    forAll(thirdBodyEfficiencies_, i)
-    {
-        beta_.append
-        (
-            Tuple2<label, scalar>
-            (
-                i,
-                thirdBodyEfficiencies_[i]
-            )
-        );
-    }
-}
+{}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -106,9 +91,11 @@ Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::operator()
 {
     const scalar k0 = k0_(p, T, c, li);
     const scalar kInf = kInf_(p, T, c, li);
-    const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+    const scalar M = thirdBodyEfficiencies_.M(c);
+    const scalar Pr = k0/kInf*M;
+    const scalar F = F_(T, Pr);
 
-    return kInf*(Pr/(1 + Pr))*F_(T, Pr);
+    return kInf*(Pr/(1 + Pr))*F;
 }
 
 
@@ -124,77 +111,78 @@ Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::ddT
 {
     const scalar k0 = k0_(p, T, c, li);
     const scalar kInf = kInf_(p, T, c, li);
-    const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+    const scalar M = thirdBodyEfficiencies_.M(c);
+    const scalar Pr = k0/kInf*M;
+    const scalar F = F_(T, Pr);
 
-    return (Pr/(1 + Pr))*F_(T, Pr)*kInf_.ddT(p, T, c, li);
+    const scalar dkInfdT = kInf_.ddT(p, T, c, li);
+
+    #if FALL_OFF_FUNCTION_JACOBIAN
+
+    const scalar dk0dT = k0_.ddT(p, T, c, li);
+    const scalar dPrdT = (M*dk0dT - Pr*dkInfdT)/kInf;
+    const scalar dFdT = F_.ddT(T, Pr, F) + F_.ddPr(T, Pr, F)*dPrdT;
+
+    return
+        dkInfdT*(Pr/(1 + Pr))*F
+      + kInf*dPrdT/sqr(1 + Pr)*F
+      + kInf*(Pr/(1 + Pr))*dFdT;
+
+    #else
+
+    return dkInfdT*(Pr/(1 + Pr))*F;
+
+    #endif
 }
 
 
 template<class ReactionRate, class FallOffFunction>
-inline void Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::dcidc
+inline bool
+Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::hasDdc() const
+{
+    return true;
+}
+
+
+template<class ReactionRate, class FallOffFunction>
+inline void Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::ddc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
     const label li,
-    scalarField& dcidc
+    scalarField& ddc
 ) const
 {
+    const scalar k0 = k0_(p, T, c, li);
+    const scalar kInf = kInf_(p, T, c, li);
     const scalar M = thirdBodyEfficiencies_.M(c);
+    const scalar Pr = k0/kInf*M;
+    const scalar F = F_(T, Pr);
 
-    if (M > small)
-    {
-        const scalar k0 = k0_(p, T, c, li);
-        const scalar kInf = kInf_(p, T, c, li);
-        const scalar Pr = k0*M/kInf;
-        const scalar F = F_(T, Pr);
+    #if FALL_OFF_FUNCTION_JACOBIAN
 
-        forAll(beta_, i)
-        {
-            const scalar dPrdci = -beta_[i].second()*k0/kInf;
-            const scalar dFdci = F_.ddc(Pr, F, dPrdci, T);
-            dcidc[i] = (dPrdci/(Pr*(1 + Pr)) + dFdci/F);
-        }
-    }
-    else
-    {
-        forAll(beta_, i)
-        {
-            dcidc[i] = 0;
-        }
-    }
-}
+    scalarField dk0dc(c.size(), 0);
+    k0_.ddc(p, T, c, li, dk0dc);
+    scalarField dkInfdc(c.size(), 0);
+    kInf_.ddc(p, T, c, li, dkInfdc);
+    tmp<scalarField> tdMdc = thirdBodyEfficiencies_.dMdc(c);
+    const scalarField& dMdc = tdMdc();
+    scalarField dPrdc((M*dk0dc - Pr*dkInfdc + k0*dMdc)/kInf);
+    const scalar dFdPr = F_.ddPr(T, Pr, F);
 
+    ddc =
+        dkInfdc*(Pr/(1 + Pr))*F
+      + kInf*dPrdc/sqr(1 + Pr)*F
+      + kInf*(Pr/(1 + Pr))*dFdPr*dPrdc;
 
-template<class ReactionRate, class FallOffFunction>
-inline Foam::scalar
-Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::dcidT
-(
-    const scalar p,
-    const scalar T,
-    const scalarField& c,
-    const label li
-) const
-{
-    const scalar M = thirdBodyEfficiencies_.M(c);
+    #else
 
-    if (M > small)
-    {
-        const scalar k0 = k0_(p, T, c, li);
-        const scalar kInf = kInf_(p, T, c, li);
+    kInf_.ddc(p, T, c, li, ddc);
 
-        const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
-        const scalar F = F_(T, Pr);
-        const scalar dPrdT =
-            Pr*(k0_.ddT(p, T, c, li)/k0 - kInf_.ddT(p, T, c, li)/kInf - 1/T);
-        const scalar dFdT = F_.ddT(Pr, F, dPrdT, T);
+    ddc *= (Pr/(1 + Pr))*F;
 
-        return (dPrdT/(Pr*(1 + Pr)) + dFdT/F);
-    }
-    else
-    {
-        return 0;
-    }
+    #endif
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H
index 9c9e291ab2..76d3cfaa01 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H
@@ -102,6 +102,7 @@ public:
         //- Post-evaluation hook
         inline void postEvaluate() const;
 
+        //- Return the rate
         inline scalar operator()
         (
             const scalar p,
@@ -110,6 +111,7 @@ public:
             const label li
         ) const;
 
+        //- The derivative of the rate w.r.t. temperature
         inline scalar ddT
         (
             const scalar p,
@@ -118,28 +120,17 @@ public:
             const label li
         ) const;
 
-        //- Third-body efficiencies (beta = 1-alpha)
-        //  non-empty only for third-body reactions
-        //  with enhanced molecularity (alpha != 1)
-        inline const List<Tuple2<label, scalar>>& beta() const;
+        //- Is the rate a function of concentration?
+        inline bool hasDdc() const;
 
-        //- Species concentration derivative of the pressure dependent term
-        inline void dcidc
+        //- The derivative of the rate w.r.t. concentration
+        inline void ddc
         (
             const scalar p,
             const scalar T,
             const scalarField& c,
             const label li,
-            scalarField& dcidc
-        ) const;
-
-        //- Temperature derivative of the pressure dependent term
-        inline scalar dcidT
-        (
-            const scalar p,
-            const scalar T,
-            const scalarField& c,
-            const label li
+            scalarField& ddc
         ) const;
 
         //- Write to stream
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H
index ec9aefe00e..739cd8650f 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H
@@ -140,33 +140,22 @@ inline Foam::scalar Foam::JanevReactionRate::ddT
 }
 
 
-inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
-Foam::JanevReactionRate::beta() const
+inline bool Foam::JanevReactionRate::hasDdc() const
 {
-    return NullObjectRef<List<Tuple2<label, scalar>>>();
+    return false;
 }
 
 
-inline void Foam::JanevReactionRate::dcidc
+inline void Foam::JanevReactionRate::ddc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
     const label li,
-    scalarField& dcidc
-) const
-{}
-
-
-inline Foam::scalar Foam::JanevReactionRate::dcidT
-(
-    const scalar p,
-    const scalar T,
-    const scalarField& c,
-    const label li
+    scalarField& ddc
 ) const
 {
-    return 0;
+    ddc = 0;
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H
index 221e240a5b..6036443385 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H
@@ -101,6 +101,7 @@ public:
         //- Post-evaluation hook
         inline void postEvaluate() const;
 
+        //- Return the rate
         inline scalar operator()
         (
             const scalar p,
@@ -109,6 +110,7 @@ public:
             const label li
         ) const;
 
+        //- The derivative of the rate w.r.t. temperature
         inline scalar ddT
         (
             const scalar p,
@@ -117,28 +119,17 @@ public:
             const label li
         ) const;
 
-        //- Third-body efficiencies (beta = 1-alpha)
-        //  non-empty only for third-body reactions
-        //  with enhanced molecularity (alpha != 1)
-        inline const List<Tuple2<label, scalar>>& beta() const;
+        //- Is the rate a function of concentration?
+        inline bool hasDdc() const;
 
-        //- Species concentration derivative of the pressure dependent term
-        inline void dcidc
+        //- The derivative of the rate w.r.t. concentration
+        inline void ddc
         (
             const scalar p,
             const scalar T,
             const scalarField& c,
             const label li,
-            scalarField& dcidc
-        ) const;
-
-        //- Temperature derivative of the pressure dependent term
-        inline scalar dcidT
-        (
-            const scalar p,
-            const scalar T,
-            const scalarField& c,
-            const label li
+            scalarField& ddc
         ) const;
 
         //- Write to stream
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H
index 884d2b89c9..eaa6e82335 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H
@@ -155,33 +155,22 @@ inline Foam::scalar Foam::LandauTellerReactionRate::ddT
 }
 
 
-inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
-Foam::LandauTellerReactionRate::beta() const
+inline bool Foam::LandauTellerReactionRate::hasDdc() const
 {
-    return NullObjectRef<List<Tuple2<label, scalar>>>();
+    return false;
 }
 
 
-inline void Foam::LandauTellerReactionRate::dcidc
+inline void Foam::LandauTellerReactionRate::ddc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
     const label li,
-    scalarField& dcidc
-) const
-{}
-
-
-inline Foam::scalar Foam::LandauTellerReactionRate::dcidT
-(
-    const scalar p,
-    const scalar T,
-    const scalarField& c,
-    const label li
+    scalarField& ddc
 ) const
 {
-    return 0;
+    ddc = 0;
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H
index 72f26e14a4..414c39824d 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H
@@ -104,6 +104,7 @@ public:
         //- Post-evaluation hook
         inline void postEvaluate() const;
 
+        //- Return the rate
         inline scalar operator()
         (
             const scalar p,
@@ -112,6 +113,7 @@ public:
             const label li
         ) const;
 
+        //- The derivative of the rate w.r.t. temperature
         inline scalar ddT
         (
             const scalar p,
@@ -120,28 +122,17 @@ public:
             const label li
         ) const;
 
-        //- Third-body efficiencies (beta = 1-alpha)
-        //  non-empty only for third-body reactions
-        //  with enhanced molecularity (alpha != 1)
-        inline const List<Tuple2<label, scalar>>& beta() const;
+        //- Is the rate a function of concentration?
+        inline bool hasDdc() const;
 
-        //- Species concentration derivative of the pressure dependent term
-        inline void dcidc
+        //- The derivative of the rate w.r.t. concentration
+        inline void ddc
         (
             const scalar p,
             const scalar T,
             const scalarField& c,
             const label li,
-            scalarField& dcidc
-        ) const;
-
-        //- Temperature derivative of the pressure dependent term
-        inline scalar dcidT
-        (
-            const scalar p,
-            const scalar T,
-            const scalarField& c,
-            const label li
+            scalarField& ddc
         ) const;
 
         //- Write to stream
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H
index b0439e5753..c2a10eda3a 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H
@@ -83,15 +83,13 @@ inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::operator()
     const scalar c1m = pow(c[r_[0]], m_[1]);
     const scalar c2m = pow(c[r_[1]], m_[2]);
 
-    const scalar TaByT0 = Ta_[0]/T;
-    const scalar TaByT1 = Ta_[1]/T;
-    const scalar TaByT2 = Ta_[2]/T;
+    const scalar k0 = A_[0]*pow(T, beta_[0])*exp(-Ta_[0]/T);
+    const scalar k1 = A_[1]*pow(T, beta_[1])*exp(-Ta_[1]/T);
+    const scalar k2 = A_[2]*pow(T, beta_[2])*exp(-Ta_[2]/T);
 
-    const scalar k0 = A_[0]*pow(T, beta_[0])*exp(-TaByT0);
-    const scalar k1 = A_[1]*pow(T, beta_[1])*exp(-TaByT1);
-    const scalar k2 = A_[2]*pow(T, beta_[2])*exp(-TaByT2);
+    const scalar b = a_ + k1*c1m + k2*c2m;
 
-    return k0/pow(a_ + k1*c1m + k2*c2m, m_[0]);
+    return k0/pow(b, m_[0]);
 }
 
 
@@ -106,50 +104,55 @@ inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::ddT
     const scalar c1m = pow(c[r_[0]], m_[1]);
     const scalar c2m = pow(c[r_[1]], m_[2]);
 
-    const scalar TaByT0 = Ta_[0]/T;
-    const scalar TaByT1 = Ta_[1]/T;
-    const scalar TaByT2 = Ta_[2]/T;
+    const scalar k0 = A_[0]*pow(T, beta_[0])*exp(-Ta_[0]/T);
+    const scalar k1 = A_[1]*pow(T, beta_[1])*exp(-Ta_[1]/T);
+    const scalar k2 = A_[2]*pow(T, beta_[2])*exp(-Ta_[2]/T);
 
-    const scalar k0 = A_[0]*pow(T, beta_[0])*exp(-TaByT0);
-    const scalar k1 = A_[1]*pow(T, beta_[1])*exp(-TaByT1);
-    const scalar k2 = A_[2]*pow(T, beta_[2])*exp(-TaByT2);
+    const scalar dk0dT = k0/T*(beta_[0] + Ta_[0]/T);
+    const scalar dk1dT = k1/T*(beta_[1] + Ta_[1]/T);
+    const scalar dk2dT = k2/T*(beta_[2] + Ta_[2]/T);
 
-    return
-       (
-           (beta_[0] + TaByT0)*k0
-         - m_[0]*k0*((beta_[1] + TaByT1)*k1*c1m + (beta_[2] + TaByT2)*k2*c2m)
-          /(a_ + k1*c1m + k2*c2m)
-       )/(pow(a_ + k1*c1m + k2*c2m, m_[0])*T);
+    const scalar b = a_ + k1*c1m + k2*c2m;
+    const scalar dbdT = dk1dT*c1m + dk2dT*c2m;
+
+    return (dk0dT - k0*m_[0]*dbdT/b)/pow(b, m_[0]);
 }
 
 
-inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
-Foam::LangmuirHinshelwoodReactionRate::beta() const
+inline bool Foam::LangmuirHinshelwoodReactionRate::hasDdc() const
 {
-    return NullObjectRef<List<Tuple2<label, scalar>>>();
+    return true;
 }
 
 
-inline void Foam::LangmuirHinshelwoodReactionRate::dcidc
+inline void Foam::LangmuirHinshelwoodReactionRate::ddc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
     const label li,
-    scalarField& dcidc
-) const
-{}
-
-
-inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::dcidT
-(
-    const scalar p,
-    const scalar T,
-    const scalarField& c,
-    const label li
+    scalarField& ddc
 ) const
 {
-    return 0;
+    const scalar c1m = pow(c[r_[0]], m_[1]);
+    const scalar c2m = pow(c[r_[1]], m_[2]);
+
+    const scalar dc1mdc1 = m_[1]*c1m/c[r_[0]];
+    const scalar dc2mdc2 = m_[2]*c2m/c[r_[1]];
+
+    const scalar k0 = A_[0]*pow(T, beta_[0])*exp(-Ta_[0]/T);
+    const scalar k1 = A_[1]*pow(T, beta_[1])*exp(-Ta_[1]/T);
+    const scalar k2 = A_[2]*pow(T, beta_[2])*exp(-Ta_[2]/T);
+
+    const scalar b = a_ + k1*c1m + k2*c2m;
+    const scalar dbdc1 = k1*dc1mdc1;
+    const scalar dbdc2 = k2*dc2mdc2;
+
+    const scalar k = k0/pow(b, m_[0]);
+
+    ddc = 0;
+    ddc[r_[0]] = - k*m_[0]/b*dbdc1;
+    ddc[r_[1]] = - k*m_[1]/b*dbdc2;
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRate.H
index 65f0f484c5..52f1c79f07 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRate.H
@@ -77,8 +77,6 @@ class MichaelisMentenReactionRate
         //- The substrate specie index
         label s_;
 
-        List<Tuple2<label, scalar>> beta_;
-
 
 public:
 
@@ -106,6 +104,7 @@ public:
         //- Post-evaluation hook
         inline void postEvaluate() const;
 
+        //- Return the rate
         inline scalar operator()
         (
             const scalar p,
@@ -114,6 +113,7 @@ public:
             const label li
         ) const;
 
+        //- The derivative of the rate w.r.t. temperature
         inline scalar ddT
         (
             const scalar p,
@@ -122,28 +122,17 @@ public:
             const label li
         ) const;
 
-        //- Third-body efficiencies (beta = 1-alpha)
-        //  non-empty only for third-body reactions
-        //  with enhanced molecularity (alpha != 1)
-        inline const List<Tuple2<label, scalar>>& beta() const;
+        //- Is the rate a function of concentration?
+        inline bool hasDdc() const;
 
-        //- Species concentration derivative of the pressure dependent term
-        inline void dcidc
+        //- The derivative of the rate w.r.t. concentration
+        inline void ddc
         (
             const scalar p,
             const scalar T,
             const scalarField& c,
             const label li,
-            scalarField& dcidc
-        ) const;
-
-        //- Temperature derivative of the pressure dependent term
-        inline scalar dcidT
-        (
-            const scalar p,
-            const scalar T,
-            const scalarField& c,
-            const label li
+            scalarField& ddc
         ) const;
 
         //- Write to stream
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRateI.H
index 6e2b2697a7..e28e713914 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRateI.H
@@ -35,9 +35,7 @@ inline Foam::MichaelisMentenReactionRate::MichaelisMentenReactionRate
     Vmax_(dict.lookup<scalar>("Vmax")),
     Km_(dict.lookup<scalar>("Km")),
     s_(st[dict.lookup("S")])
-{
-    beta_.append(Tuple2<label, scalar>(s_, 1.0));
-}
+{}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -74,35 +72,23 @@ inline Foam::scalar Foam::MichaelisMentenReactionRate::ddT
 }
 
 
-inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
-Foam::MichaelisMentenReactionRate::beta() const
+inline bool Foam::MichaelisMentenReactionRate::hasDdc() const
 {
-    return beta_;
+    return true;
 }
 
 
-inline void Foam::MichaelisMentenReactionRate::dcidc
+inline void Foam::MichaelisMentenReactionRate::ddc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
     const label li,
-    scalarField& dcidc
+    scalarField& ddc
 ) const
 {
-    dcidc[0] = -1.0/(Km_ + c[s_]);
-}
-
-
-inline Foam::scalar Foam::MichaelisMentenReactionRate::dcidT
-(
-    const scalar p,
-    const scalar T,
-    const scalarField& c,
-    const label li
-) const
-{
-    return 0;
+    ddc = 0;
+    ddc[s_] = - Vmax_/sqr(Km_ + c[s_]);
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunction.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunction.H
index 4e28931ee3..caab63bb3b 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunction.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunction.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -82,18 +82,16 @@ public:
 
         inline scalar ddT
         (
+            const scalar T,
             const scalar Pr,
-            const scalar F,
-            const scalar dPrdT,
-            const scalar T
+            const scalar F
         ) const;
 
-        inline scalar ddc
+        inline scalar ddPr
         (
+            const scalar T,
             const scalar Pr,
-            const scalar F,
-            const scalar dPrdc,
-            const scalar T
+            const scalar F
         ) const;
 
         //- Write to stream
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunctionI.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunctionI.H
index 85231025dc..7702096f73 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunctionI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunctionI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -50,7 +50,6 @@ inline Foam::scalar Foam::LindemannFallOffFunction::operator()
 
 inline Foam::scalar Foam::LindemannFallOffFunction::ddT
 (
-    const scalar,
     const scalar,
     const scalar,
     const scalar
@@ -60,9 +59,8 @@ inline Foam::scalar Foam::LindemannFallOffFunction::ddT
 }
 
 
-inline Foam::scalar Foam::LindemannFallOffFunction::ddc
+inline Foam::scalar Foam::LindemannFallOffFunction::ddPr
 (
-    const scalar,
     const scalar,
     const scalar,
     const scalar
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunction.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunction.H
index 65dd2d6376..58daa57db4 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunction.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunction.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -93,18 +93,16 @@ public:
 
         inline scalar ddT
         (
+            const scalar T,
             const scalar Pr,
-            const scalar F,
-            const scalar dPrdT,
-            const scalar T
+            const scalar F
         ) const;
 
-        inline scalar ddc
+        inline scalar ddPr
         (
+            const scalar T,
             const scalar Pr,
-            const scalar F,
-            const scalar dPrdc,
-            const scalar T
+            const scalar F
         ) const;
 
         //- Write to stream
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunctionI.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunctionI.H
index 9222f919f9..ff4ab77851 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunctionI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunctionI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -60,45 +60,52 @@ inline Foam::scalar Foam::SRIFallOffFunction::operator()
     const scalar Pr
 ) const
 {
-    scalar X = 1.0/(1 + sqr(log10(max(Pr, small))));
-    return d_*pow(a_*exp(-b_/T) + exp(-T/c_), X)*pow(T, e_);
+    const scalar logPr = log10(max(Pr, small));
+
+    const scalar X = 1/(1 + sqr(logPr));
+
+    const scalar psi = a_*exp(-b_/T) + exp(-T/c_);
+
+    return d_*pow(psi, X)*pow(T, e_);
 }
 
 
 inline Foam::scalar Foam::SRIFallOffFunction::ddT
 (
+    const scalar T,
     const scalar Pr,
-    const scalar F,
-    const scalar dPrdT,
-    const scalar T
+    const scalar F
 ) const
 {
-    scalar X = 1.0/(1 + sqr(log10(max(Pr, small))));
-    scalar dXdPr = -X*X*2*log10(Pr)/Pr/log(10.0);
-    return
-        F
-       *(
-            e_/T
-          + X
-           *(a_*b_*exp(-b_/T)/sqr(T) - exp(-T/c_)/c_)
-           /(a_*exp(-b_/T) + exp(-T/c_))
-          + dXdPr*dPrdT*log(a_*exp(-b_/T) + exp(-T/c_))
-        );
+    const scalar logPr = log10(max(Pr, small));
+
+    const scalar X = 1/(1 + sqr(logPr));
+
+    const scalar psi = a_*exp(-b_/T) + exp(-T/c_);
+    const scalar dpsidT = a_*b_/sqr(T)*exp(-b_/T) - 1/c_*exp(-T/c_);
+
+    return F*(X/psi*dpsidT + e_/T);
 }
 
 
-inline Foam::scalar Foam::SRIFallOffFunction::ddc
+inline Foam::scalar Foam::SRIFallOffFunction::ddPr
 (
+    const scalar T,
     const scalar Pr,
-    const scalar F,
-    const scalar dPrdc,
-    const scalar T
+    const scalar F
 ) const
 {
-    scalar X = 1.0/(1 + sqr(log10(max(Pr, small))));
-    scalar dXdPr = -X*X*2*log10(Pr)/Pr/log(10.0);
-    scalar dXdc = dXdPr*dPrdc;
-    return F*dXdc*log(a_*exp(-b_/T) + exp(-T/c_));
+    static const scalar logTen = log(scalar(10));
+
+    const scalar logPr = log10(max(Pr, small));
+    const scalar dlogPrdPr = Pr >= small ? 1/(logTen*Pr) : 0;
+
+    const scalar X = 1/(1 + sqr(logPr));
+    const scalar dXdPr = -sqr(X)*2*logPr*dlogPrdPr;
+
+    const scalar psi = a_*exp(-b_/T) + exp(-T/c_);
+
+    return F*log(psi)*dXdPr;
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunction.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunction.H
index b1edefea50..11ec779101 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunction.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunction.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -93,18 +93,16 @@ public:
 
         inline scalar ddT
         (
+            const scalar T,
             const scalar Pr,
-            const scalar F,
-            const scalar dPrdT,
-            const scalar T
+            const scalar F
         ) const;
 
-        inline scalar ddc
+        inline scalar ddPr
         (
+            const scalar T,
             const scalar Pr,
-            const scalar F,
-            const scalar dPrdc,
-            const scalar T
+            const scalar F
         ) const;
 
         //- Write to stream
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunctionI.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunctionI.H
index 1fad3a8f76..229ae5c189 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunctionI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunctionI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,116 +57,103 @@ inline Foam::scalar Foam::TroeFallOffFunction::operator()
     const scalar Pr
 ) const
 {
-    scalar logFcent = log10
-    (
-        max
-        (
-            (1 - alpha_)*exp(-T/Tsss_) + alpha_*exp(-T/Ts_) + exp(-Tss_/T),
-            small
-        )
-    );
+    const scalar logPr = log10(max(Pr, small));
 
-    scalar c = -0.4 - 0.67*logFcent;
+    const scalar Fcent =
+        (1 - alpha_)*exp(-T/Tsss_)
+      + alpha_*exp(-T/Ts_)
+      + exp(-Tss_/T);
+
+    const scalar logFcent = log10(max(Fcent, small));
+
+    const scalar c = -0.4 - 0.67*logFcent;
     static const scalar d = 0.14;
-    scalar n = 0.75 - 1.27*logFcent;
+    const scalar n = 0.75 - 1.27*logFcent;
 
-    scalar logPr = log10(max(Pr, small));
-    return pow(10, logFcent/(1 + sqr((logPr + c)/(n - d*(logPr + c)))));
+    const scalar x1 = n - d*(logPr + c);
+    const scalar x2 = (logPr + c)/x1;
+    const scalar x3 = 1 + sqr(x2);
+    const scalar x4 = logFcent/x3;
+
+    return pow(10, x4);
 }
 
 
 inline Foam::scalar Foam::TroeFallOffFunction::ddT
 (
+    const scalar T,
     const scalar Pr,
-    const scalar F,
-    const scalar dPrdT,
-    const scalar T
+    const scalar F
 ) const
 {
-    scalar logPr = log10(max(Pr, small));
-    scalar logTen = log(10.0);
-    scalar Fcent =
-    (
-        max
-        (
-            (1 - alpha_)*exp(-T/Tsss_) + alpha_*exp(-T/Ts_) + exp(-Tss_/T),
-            small
-        )
-    );
-    scalar logFcent = log10(Fcent);
+    static const scalar logTen = log(scalar(10));
 
-    scalar dFcentdT =
-    (
-        (alpha_ - 1)*exp(-T/Tsss_)/Tsss_
-      - alpha_*exp(-T/Ts_)/Ts_
-      + Tss_*exp(-Tss_/T)/sqr(T)
-    );
+    const scalar logPr = log10(max(Pr, small));
 
-    scalar d = 0.14;
-    scalar dlogFcentdT = dFcentdT/Fcent/logTen;
-    scalar c = -0.4 - 0.67*logFcent;
-    scalar dcdT = -0.67*dlogFcentdT;
-    scalar n = 0.75 - 1.27*logFcent;
-    scalar dndT = -1.27*dlogFcentdT;
+    const scalar Fcent =
+        (1 - alpha_)*exp(-T/Tsss_)
+      + alpha_*exp(-T/Ts_)
+      + exp(-Tss_/T);
+    const scalar dFcentdT =
+      - (1 - alpha_)/Tsss_*exp(-T/Tsss_)
+      - alpha_/Ts_*exp(-T/Ts_)
+      + Tss_/sqr(T)*exp(-Tss_/T);
 
-    scalar dlogPrdT = dPrdT/Pr/logTen;
+    const scalar logFcent = log10(max(Fcent, small));
+    const scalar dlogFcentdT = Fcent >= small ? dFcentdT/Fcent/logTen : 0;
 
-    scalar dParentdT =
-        2.0*(logPr + c)/sqr(n - d*(logPr + c))
-       *(
-            (dlogPrdT + dcdT)
-          - (logPr + c)*(dndT - d*(dlogPrdT + dcdT))/(n - d*(logPr + c))
-        );
+    const scalar c = -0.4 - 0.67*logFcent;
+    const scalar dcdT = -0.67*dlogFcentdT;
+    static const scalar d = 0.14;
+    const scalar n = 0.75 - 1.27*logFcent;
+    const scalar dndT = - 1.27*dlogFcentdT;
 
-    return
-    (
-        F*logTen
-       *(
-            dlogFcentdT/(1.0 + sqr((logPr + c)/(n - d*(logPr + c))))
-          - logFcent*dParentdT/sqr(1.0 + sqr((logPr + c)/(n - d*(logPr + c))))
-        )
-    );
+    const scalar x1 = n - d*(logPr + c);
+    const scalar dx1dT = dndT - d*dcdT;
+    const scalar x2 = (logPr + c)/x1;
+    const scalar dx2dT = (dcdT - x2*dx1dT)/x1;
+    const scalar x3 = 1 + sqr(x2);
+    const scalar dx3dT = 2*x2*dx2dT;
+    const scalar x4 = logFcent/x3;
+    const scalar dx4dT = (dlogFcentdT - x4*dx3dT)/x3;
+
+    return logTen*F*dx4dT;
 }
 
 
-inline Foam::scalar Foam::TroeFallOffFunction::ddc
+inline Foam::scalar Foam::TroeFallOffFunction::ddPr
 (
+    const scalar T,
     const scalar Pr,
-    const scalar F,
-    const scalar dPrdc,
-    const scalar T
+    const scalar F
 ) const
 {
-    scalar logPr = log10(max(Pr, small));
-    scalar logTen = log(10.0);
-    scalar logFcent = log10
-    (
-        max
-        (
-            (1 - alpha_)*exp(-T/Tsss_) + alpha_*exp(-T/Ts_) + exp(-Tss_/T),
-            small
-        )
-    );
+    static const scalar logTen = log(scalar(10));
 
-    scalar dlogPrdc = dPrdc/Pr/logTen;
-    scalar d = 0.14;
-    scalar c = -0.4 - 0.67*logFcent;
-    scalar n = 0.75 - 1.27*logFcent;
+    const scalar logPr = log10(max(Pr, small));
+    const scalar dlogPrdPr = Pr >= small ? 1/(logTen*Pr) : 0;
 
-    scalar dParentdc =
-        2.0*(logPr + c)/sqr(n - d*(logPr + c))
-       *(
-            (dlogPrdc)
-          - (logPr + c)*(-d*(dlogPrdc))/(n - d*(logPr + c))
-        );
+    const scalar Fcent =
+        (1 - alpha_)*exp(-T/Tsss_)
+      + alpha_*exp(-T/Ts_)
+      + exp(-Tss_/T);
 
-    return
-    (
-        F*logTen
-       *(
-          - logFcent*dParentdc/sqr(1.0 + sqr((logPr + c)/(n - d*(logPr + c))))
-        )
-    );
+    const scalar logFcent = log10(max(Fcent, small));
+
+    const scalar c = -0.4 - 0.67*logFcent;
+    static const scalar d = 0.14;
+    const scalar n = 0.75 - 1.27*logFcent;
+
+    const scalar x1 = n - d*(logPr + c);
+    const scalar dx1dPr = -d*dlogPrdPr;
+    const scalar x2 = (logPr + c)/x1;
+    const scalar dx2dPr = (dlogPrdPr - x2*dx1dPr)/x1;
+    const scalar x3 = 1 + sqr(x2);
+    const scalar dx3dPr = 2*x2*dx2dPr;
+    const scalar x4 = logFcent/x3;
+    const scalar dx4dPr = -x4*dx3dPr/x3;
+
+    return logTen*F*dx4dPr;
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H
index 6e903200cf..b8f8abd81f 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H
@@ -102,6 +102,7 @@ public:
         //- Post-evaluation hook
         inline void postEvaluate() const;
 
+        //- Return the rate
         inline scalar operator()
         (
             const scalar p,
@@ -110,6 +111,7 @@ public:
             const label li
         ) const;
 
+        //- The derivative of the rate w.r.t. temperature
         inline scalar ddT
         (
             const scalar p,
@@ -118,28 +120,17 @@ public:
             const label li
         ) const;
 
-        //- Third-body efficiencies (beta = 1-alpha)
-        //  non-empty only for third-body reactions
-        //  with enhanced molecularity (alpha != 1)
-        inline const List<Tuple2<label, scalar>>& beta() const;
+        //- Is the rate a function of concentration?
+        inline bool hasDdc() const;
 
-        //- Species concentration derivative of the pressure dependent term
-        inline void dcidc
+        //- The derivative of the rate w.r.t. concentration
+        inline void ddc
         (
             const scalar p,
             const scalar T,
             const scalarField& c,
             const label li,
-            scalarField& dcidc
-        ) const;
-
-        //- Temperature derivative of the pressure dependent term
-        inline scalar dcidT
-        (
-            const scalar p,
-            const scalar T,
-            const scalarField& c,
-            const label li
+            scalarField& ddc
         ) const;
 
         //- Write to stream
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H
index 14ee608f74..6f174cc05b 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H
@@ -122,33 +122,22 @@ inline Foam::scalar Foam::powerSeriesReactionRate::ddT
 }
 
 
-inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
-Foam::powerSeriesReactionRate::beta() const
+inline bool Foam::powerSeriesReactionRate::hasDdc() const
 {
-    return NullObjectRef<List<Tuple2<label, scalar>>>();
+    return false;
 }
 
 
-inline void Foam::powerSeriesReactionRate::dcidc
+inline void Foam::powerSeriesReactionRate::ddc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
     const label li,
-    scalarField& dcidc
-) const
-{}
-
-
-inline Foam::scalar Foam::powerSeriesReactionRate::dcidT
-(
-    const scalar p,
-    const scalar T,
-    const scalarField& c,
-    const label li
+    scalarField& ddc
 ) const
 {
-    return 0;
+    ddc = 0;
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H
index 596ff47922..2b067f9a23 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H
@@ -61,7 +61,6 @@ class thirdBodyArrheniusReactionRate
     // Private Data
 
         thirdBodyEfficiencies thirdBodyEfficiencies_;
-        List<Tuple2<label, scalar>> beta_;
 
 
 public:
@@ -99,6 +98,7 @@ public:
         //- Post-evaluation hook
         inline void postEvaluate() const;
 
+        //- Return the rate
         inline scalar operator()
         (
             const scalar p,
@@ -107,6 +107,7 @@ public:
             const label li
         ) const;
 
+        //- The derivative of the rate w.r.t. temperature
         inline scalar ddT
         (
             const scalar p,
@@ -115,31 +116,17 @@ public:
             const label li
         ) const;
 
-        //- Third-body efficiencies (beta = 1-alpha)
-        //  non-empty only for third-body reactions
-        //  with enhanced molecularity (alpha != 1)
-        inline const List<Tuple2<label, scalar>>& beta() const
-        {
-            return beta_;
-        }
+        //- Is the rate a function of concentration?
+        inline bool hasDdc() const;
 
-        //- Species concentration derivative of the pressure dependent term
-        inline void dcidc
+        //- The derivative of the rate w.r.t. concentration
+        inline void ddc
         (
             const scalar p,
             const scalar T,
             const scalarField& c,
             const label li,
-            scalarField& dcidc
-        ) const;
-
-        //- Temperature derivative of the pressure dependent term
-        inline scalar dcidT
-        (
-            const scalar p,
-            const scalar T,
-            const scalarField& c,
-            const label li
+            scalarField& ddc
         ) const;
 
         //- Write to stream
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H
index 2fabaf78c3..23a1d3fd76 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H
@@ -35,12 +35,7 @@ inline Foam::thirdBodyArrheniusReactionRate::thirdBodyArrheniusReactionRate
 :
     ArrheniusReactionRate(A, beta, Ta),
     thirdBodyEfficiencies_(tbes)
-{
-    forAll(tbes, i)
-    {
-        beta_.append(Tuple2<label, scalar>(i, tbes[i]));
-    }
-}
+{}
 
 
 inline Foam::thirdBodyArrheniusReactionRate::thirdBodyArrheniusReactionRate
@@ -55,19 +50,7 @@ inline Foam::thirdBodyArrheniusReactionRate::thirdBodyArrheniusReactionRate
         dict
     ),
     thirdBodyEfficiencies_(species, dict)
-{
-    forAll(thirdBodyEfficiencies_, i)
-    {
-        beta_.append
-        (
-            Tuple2<label, scalar>
-            (
-                i,
-                thirdBodyEfficiencies_[i]
-            )
-        );
-    }
-}
+{}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -112,32 +95,25 @@ inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::ddT
 }
 
 
-inline void Foam::thirdBodyArrheniusReactionRate::dcidc
+inline bool Foam::thirdBodyArrheniusReactionRate::hasDdc() const
+{
+    return true;
+}
+
+
+inline void Foam::thirdBodyArrheniusReactionRate::ddc
 (
     const scalar p,
     const scalar T,
     const scalarField& c,
     const label li,
-    scalarField& dcidc
+    scalarField& ddc
 ) const
 {
-    scalar M = thirdBodyEfficiencies_.M(c);
-    forAll(beta_, i)
-    {
-        dcidc[i] = beta_[i].second()/max(M, small);
-    }
-}
+    const scalar k = ArrheniusReactionRate::operator()(p, T, c, li);
 
-
-inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::dcidT
-(
-    const scalar p,
-    const scalar T,
-    const scalarField& c,
-    const label li
-) const
-{
-    return -1.0/T;
+    ddc = thirdBodyEfficiencies_.dMdc(c);
+    ddc *= k;
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyEfficiencies/thirdBodyEfficiencies.H b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyEfficiencies/thirdBodyEfficiencies.H
index d2eab2cfa9..6b177697a4 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyEfficiencies/thirdBodyEfficiencies.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyEfficiencies/thirdBodyEfficiencies.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,7 @@ SourceFiles
 #ifndef thirdBodyEfficiencies_H
 #define thirdBodyEfficiencies_H
 
-#include "scalarList.H"
+#include "scalarField.H"
 #include "speciesTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -55,7 +55,7 @@ Ostream& operator<<(Ostream&, const thirdBodyEfficiencies&);
 
 class thirdBodyEfficiencies
 :
-    public scalarList
+    public scalarField
 {
     // Private Data
 
@@ -84,7 +84,11 @@ public:
     // Member Functions
 
         //- Calculate and return M, the concentration of the third-bodies
-        inline scalar M(const scalarList& c) const;
+        inline scalar M(const scalarField& c) const;
+
+        //- Calculate and return the derivative of M, w.r.t. the species
+        //  concentrations
+        inline tmp<scalarField> dMdc(const scalarField& c) const;
 
         //- Write to stream
         inline void write(Ostream& os) const;
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyEfficiencies/thirdBodyEfficienciesI.H b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyEfficiencies/thirdBodyEfficienciesI.H
index 32312cd67c..a5c6f7d678 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyEfficiencies/thirdBodyEfficienciesI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyEfficiencies/thirdBodyEfficienciesI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,7 +33,7 @@ inline Foam::thirdBodyEfficiencies::thirdBodyEfficiencies
     const scalarList& efficiencies
 )
 :
-    scalarList(efficiencies),
+    scalarField(efficiencies),
     species_(species)
 {
     if (size() != species_.size())
@@ -52,7 +52,7 @@ inline Foam::thirdBodyEfficiencies::thirdBodyEfficiencies
     const dictionary& dict
 )
 :
-    scalarList(species.size()),
+    scalarField(species.size()),
     species_(species)
 {
     if (dict.found("coeffs"))
@@ -74,14 +74,14 @@ inline Foam::thirdBodyEfficiencies::thirdBodyEfficiencies
     else
     {
         scalar defaultEff = dict.lookup<scalar>("defaultEfficiency");
-        scalarList::operator=(defaultEff);
+        scalarField::operator=(defaultEff);
     }
 }
 
 
 // * * * * * * * * * * * * * * * Member functions  * * * * * * * * * * * * * //
 
-inline Foam::scalar Foam::thirdBodyEfficiencies::M(const scalarList& c) const
+inline Foam::scalar Foam::thirdBodyEfficiencies::M(const scalarField& c) const
 {
     scalar M = 0;
     forAll(*this, i)
@@ -93,6 +93,13 @@ inline Foam::scalar Foam::thirdBodyEfficiencies::M(const scalarList& c) const
 }
 
 
+inline Foam::tmp<Foam::scalarField>
+Foam::thirdBodyEfficiencies::dMdc(const scalarField& c) const
+{
+    return *this;
+}
+
+
 inline void Foam::thirdBodyEfficiencies::write(Ostream& os) const
 {
     List<Tuple2<word, scalar>> coeffs(species_.size());
diff --git a/src/thermophysicalModels/specie/reaction/specieCoeffs/specieCoeffs.C b/src/thermophysicalModels/specie/reaction/specieCoeffs/specieCoeffs.C
index 5198d2a53a..072529d157 100644
--- a/src/thermophysicalModels/specie/reaction/specieCoeffs/specieCoeffs.C
+++ b/src/thermophysicalModels/specie/reaction/specieCoeffs/specieCoeffs.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2019-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2019-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -191,7 +191,7 @@ void Foam::specieCoeffs::reactionStr
             reaction << scs[i].stoichCoeff;
         }
         reaction << species[scs[i].index];
-        if (mag(scs[i].exponent - scs[i].stoichCoeff) > small)
+        if (mag(scs[i].exponent.operator scalar() - scs[i].stoichCoeff) > small)
         {
             reaction << "^" << scs[i].exponent;
         }
diff --git a/src/thermophysicalModels/specie/reaction/specieCoeffs/specieCoeffs.H b/src/thermophysicalModels/specie/reaction/specieCoeffs/specieCoeffs.H
index 2be13b6fdf..b5601fd5a2 100644
--- a/src/thermophysicalModels/specie/reaction/specieCoeffs/specieCoeffs.H
+++ b/src/thermophysicalModels/specie/reaction/specieCoeffs/specieCoeffs.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2019-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2019-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,6 +36,7 @@ SourceFiles
 #ifndef specieCoeffs_H
 #define specieCoeffs_H
 
+#include "specieExponent.H"
 #include "speciesTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -47,7 +48,6 @@ namespace Foam
 class specieCoeffs;
 Ostream& operator<<(Ostream&, const specieCoeffs&);
 
-
 /*---------------------------------------------------------------------------*\
                         Class specieCoeffs Declaration
 \*---------------------------------------------------------------------------*/
@@ -56,62 +56,83 @@ class specieCoeffs
 {
 public:
 
-    label index;
-    scalar stoichCoeff;
-    scalar exponent;
+    // Public Data
 
-    specieCoeffs()
-    :
-        index(-1),
-        stoichCoeff(0),
-        exponent(1)
-    {}
+        //- Index of the specie
+        label index;
 
-    specieCoeffs(const speciesTable& species, Istream& is);
+        //- Stoichiometric coefficient
+        scalar stoichCoeff;
 
-    bool operator==(const specieCoeffs& sc) const
-    {
-        return index == sc.index;
-    }
+        //- Exponent of the specie concentration
+        specieExponent exponent;
 
-    bool operator!=(const specieCoeffs& sc) const
-    {
-        return index != sc.index;
-    }
 
-    friend Ostream& operator<<(Ostream& os, const specieCoeffs& sc)
-    {
-        os  << sc.index << token::SPACE
-            << sc.stoichCoeff << token::SPACE
-            << sc.exponent;
-        return os;
-    }
+    // Constructors
 
-    //- Construct the left- and right-hand-side reaction coefficients
-    static void setLRhs
-    (
-        Istream&,
-        const speciesTable&,
-        List<specieCoeffs>& lhs,
-        List<specieCoeffs>& rhs
-    );
+        //- Construct null
+        specieCoeffs()
+        :
+            index(-labelMax),
+            stoichCoeff(NaN),
+            exponent()
+        {}
 
-    //- Write the string representation of the specieCoeffs list
-    static void reactionStr
-    (
-        OStringStream& reaction,
-        const speciesTable&,
-        const List<specieCoeffs>& scs
-    );
+        //- Construct from species table and input stream
+        specieCoeffs(const speciesTable& species, Istream& is);
 
-    //- Return string representation of reaction
-    static string reactionStr
-    (
-        OStringStream& reaction,
-        const speciesTable&,
-        const List<specieCoeffs>& lhs,
-        const List<specieCoeffs>& rhs
-    );
+
+    // Member Functions
+
+        //- Construct the left- and right-hand-side reaction coefficients
+        static void setLRhs
+        (
+            Istream&,
+            const speciesTable&,
+            List<specieCoeffs>& lhs,
+            List<specieCoeffs>& rhs
+        );
+
+        //- Write the string representation of the specieCoeffs list
+        static void reactionStr
+        (
+            OStringStream& reaction,
+            const speciesTable&,
+            const List<specieCoeffs>& scs
+        );
+
+        //- Return string representation of reaction
+        static string reactionStr
+        (
+            OStringStream& reaction,
+            const speciesTable&,
+            const List<specieCoeffs>& lhs,
+            const List<specieCoeffs>& rhs
+        );
+
+
+    // Member Operators
+
+        //- Equality comparison
+        bool operator==(const specieCoeffs& sc) const
+        {
+            return index == sc.index;
+        }
+
+        //- Inequality comparison
+        bool operator!=(const specieCoeffs& sc) const
+        {
+            return index != sc.index;
+        }
+
+        //- Write to output stream
+        friend Ostream& operator<<(Ostream& os, const specieCoeffs& sc)
+        {
+            os  << sc.index << token::SPACE
+                << sc.stoichCoeff << token::SPACE
+                << sc.exponent;
+            return os;
+        }
 };
 
 
diff --git a/src/thermophysicalModels/specie/reaction/specieExponent/specieExponent.H b/src/thermophysicalModels/specie/reaction/specieExponent/specieExponent.H
new file mode 100644
index 0000000000..8878b15cb7
--- /dev/null
+++ b/src/thermophysicalModels/specie/reaction/specieExponent/specieExponent.H
@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::specieExponent
+
+Description
+
+SourceFiles
+    specieExponentI.H
+    specieExponent.C
+    specieExponentIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef specieExponent_H
+#define specieExponent_H
+
+#include "label.H"
+#include "scalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class Ostream;
+
+// Forward declaration of friend functions and operators
+class specieExponent;
+scalar pow(scalar x, const specieExponent& e);
+specieExponent operator+(const specieExponent& a, const specieExponent& b);
+specieExponent operator-(const specieExponent& a, const specieExponent& b);
+Ostream& operator<<(Ostream&, const specieExponent&);
+
+
+/*---------------------------------------------------------------------------*\
+                          Class specieExponent Declaration
+\*---------------------------------------------------------------------------*/
+
+class specieExponent
+{
+    // Private Data
+
+        //- Marker value to indicate that there is no possible integer
+        //  representation of this exponent
+        static const label noIntegerExponent_ = labelMax;
+
+        //- Integer exponent
+        label integerExponent_;
+
+        //- Scalar exponent
+        scalar scalarExponent_;
+
+
+    // Private Member Functions
+
+        //- Return whether or not this exponent has an integer representation
+        inline bool hasIntegerExponent() const;
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline specieExponent();
+
+        //- Construct from integer
+        inline specieExponent(const label integerExponent);
+
+        //- Construct from scalar
+        inline specieExponent(const scalar scalarExponent);
+
+
+    // Member Operators
+
+        //- Cast to scalar
+        inline operator scalar() const;
+
+        //- Assign to integer
+        inline specieExponent& operator=(const label integerExponent);
+
+        //- Assign to scalar
+        inline specieExponent& operator=(const scalar scalarExponent);
+
+        //- Negate a specie exponent
+        inline specieExponent operator-() const;
+
+
+    // Friend Functions
+
+        //- Compute the power of a number to a specie exponent
+        inline friend scalar pow(const scalar x, const specieExponent& e);
+
+
+    // Friend Operators
+
+        //- Sum two specie exponents
+        inline friend specieExponent operator+
+        (
+            const specieExponent& a,
+            const specieExponent& b
+        );
+
+        //- Subtract two specie exponents
+        inline friend specieExponent operator-
+        (
+            const specieExponent& a,
+            const specieExponent& b
+        );
+
+
+    // IOstream Operators
+
+        //- Write to output stream
+        inline friend Ostream& operator<<(Ostream& os, const specieExponent& e);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "specieExponentI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/reaction/specieExponent/specieExponentI.H b/src/thermophysicalModels/specie/reaction/specieExponent/specieExponentI.H
new file mode 100644
index 0000000000..eb9068e982
--- /dev/null
+++ b/src/thermophysicalModels/specie/reaction/specieExponent/specieExponentI.H
@@ -0,0 +1,174 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "specieExponent.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+inline bool Foam::specieExponent::hasIntegerExponent() const
+{
+    return integerExponent_ != noIntegerExponent_;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::specieExponent::specieExponent()
+:
+    integerExponent_(noIntegerExponent_),
+    scalarExponent_(NaN)
+{}
+
+
+inline Foam::specieExponent::specieExponent(const label integerExponent)
+:
+    integerExponent_(integerExponent),
+    scalarExponent_(integerExponent)
+{}
+
+
+inline Foam::specieExponent::specieExponent(const scalar scalarExponent)
+:
+    integerExponent_(labelMax),
+    scalarExponent_(scalarExponent)
+{
+    const label integerExponent = floor(scalarExponent);
+    if (integerExponent == scalarExponent)
+    {
+        integerExponent_ = integerExponent;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline Foam::specieExponent::operator scalar() const
+{
+    return scalarExponent_;
+}
+
+
+inline Foam::specieExponent&
+Foam::specieExponent::operator=(const label integerExponent)
+{
+    return *this = specieExponent(integerExponent);
+}
+
+
+inline Foam::specieExponent&
+Foam::specieExponent::operator=(const scalar scalarExponent)
+{
+    return *this = specieExponent(scalarExponent);
+}
+
+
+inline Foam::specieExponent Foam::specieExponent::operator-() const
+{
+    if (hasIntegerExponent())
+    {
+        return -integerExponent_;
+    }
+    else
+    {
+        return -scalarExponent_;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
+
+inline Foam::scalar Foam::pow(const scalar x, const specieExponent& e)
+{
+    if (e.hasIntegerExponent())
+    {
+        if (e.integerExponent_ == 0)
+        {
+            return 1;
+        }
+        else
+        {
+            scalar xx = e.integerExponent_ > 0 ? x : 1/x;
+
+            scalar y = 1;
+
+            for (label i = e.integerExponent_; i != 0; i /= 2)
+            {
+                if (i % 2)
+                {
+                    y *= xx;
+                }
+
+                xx *= xx;
+            }
+
+            return y;
+        }
+    }
+    else
+    {
+        return pow(x, e.scalarExponent_);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+inline Foam::specieExponent Foam::operator+
+(
+    const specieExponent& a,
+    const specieExponent& b
+)
+{
+    if (a.hasIntegerExponent() && b.hasIntegerExponent())
+    {
+        return a.integerExponent_ + b.integerExponent_;
+    }
+    else
+    {
+        return a.scalarExponent_ + b.scalarExponent_;
+    }
+}
+
+
+inline Foam::specieExponent Foam::operator-
+(
+    const specieExponent& a,
+    const specieExponent& b
+)
+{
+    return a + (-b);
+}
+
+
+// * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
+
+inline Foam::Ostream& Foam::operator<<(Ostream& os, const specieExponent& e)
+{
+    os  << e.operator scalar();
+    return os;
+}
+
+
+// ************************************************************************* //