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;
+}
+
+
+// ************************************************************************* //