diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
index a8864b5256..9308097051 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
@@ -255,8 +255,7 @@ Foam::StandardChemistryModel<ThermoType>::tc() const
 
     const label nReaction = reactions_.size();
 
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
+    scalar omegaf, omegar;
 
     if (this->chemistry_)
     {
@@ -279,12 +278,11 @@ Foam::StandardChemistryModel<ThermoType>::tc() const
             forAll(reactions_, i)
             {
                 const Reaction<ThermoType>& R = reactions_[i];
-
-                R.omega(pi, Ti, c_, celli, pf, cf, lRef, pr, cr, rRef);
+                R.omega(pi, Ti, c_, celli, omegaf, omegar);
 
                 forAll(R.rhs(), s)
                 {
-                    tc[celli] += R.rhs()[s].stoichCoeff*pf*cf;
+                    tc[celli] += R.rhs()[s].stoichCoeff*omegaf;
                 }
             }
 
@@ -360,8 +358,7 @@ Foam::StandardChemistryModel<ThermoType>::calculateRR
 
     reactionEvaluationScope scope(*this);
 
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
+    scalar omegaf, omegar;
 
     forAll(rho, celli)
     {
@@ -376,16 +373,13 @@ Foam::StandardChemistryModel<ThermoType>::calculateRR
         }
 
         const Reaction<ThermoType>& R = reactions_[ri];
-        const scalar omegai = R.omega
-        (
-            pi, Ti, c_, celli, pf, cf, lRef, pr, cr, rRef
-        );
+        const scalar omegaI = R.omega(pi, Ti, c_, celli, omegaf, omegar);
 
         forAll(R.lhs(), s)
         {
             if (si == R.lhs()[s].index)
             {
-                RR[celli] -= R.lhs()[s].stoichCoeff*omegai;
+                RR[celli] -= R.lhs()[s].stoichCoeff*omegaI;
             }
         }
 
@@ -393,7 +387,7 @@ Foam::StandardChemistryModel<ThermoType>::calculateRR
         {
             if (si == R.rhs()[s].index)
             {
-                RR[celli] += R.rhs()[s].stoichCoeff*omegai;
+                RR[celli] += R.rhs()[s].stoichCoeff*omegaI;
             }
         }
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
index dd14ae8e29..4d51cf081b 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
@@ -135,8 +135,7 @@ void Foam::TDACChemistryModel<ThermoType>::omega
 {
     const bool reduced = mechRed_->active();
 
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
+    scalar omegaf, omegar;
 
     dcdt = Zero;
 
@@ -145,33 +144,26 @@ void Foam::TDACChemistryModel<ThermoType>::omega
         if (!reactionsDisabled_[i])
         {
             const Reaction<ThermoType>& R = this->reactions_[i];
-
-            scalar omegai = R.omega
-            (
-                p, T, c, li, pf, cf, lRef, pr, cr, rRef
-            );
+            const scalar omegaI = R.omega(p, T, c, li, omegaf, omegar);
 
             forAll(R.lhs(), s)
             {
-                label si = R.lhs()[s].index;
-                if (reduced)
-                {
-                    si = completeToSimplifiedIndex_[si];
-                }
-
+                const label si =
+                    reduced
+                  ? completeToSimplifiedIndex_[R.lhs()[s].index]
+                  : R.lhs()[s].index;
                 const scalar sl = R.lhs()[s].stoichCoeff;
-                dcdt[si] -= sl*omegai;
+                dcdt[si] -= sl*omegaI;
             }
+
             forAll(R.rhs(), s)
             {
-                label si = R.rhs()[s].index;
-                if (reduced)
-                {
-                    si = completeToSimplifiedIndex_[si];
-                }
-
+                const label si =
+                    reduced
+                  ? completeToSimplifiedIndex_[R.rhs()[s].index]
+                  : R.rhs()[s].index;
                 const scalar sr = R.rhs()[s].stoichCoeff;
-                dcdt[si] += sr*omegai;
+                dcdt[si] += sr*omegaI;
             }
         }
     }
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C
index 68f179e2c3..a4dc6fc309 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C
@@ -263,14 +263,13 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
     // Index of the other species involved in the rABNum
     RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
 
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
     forAll(this->chemistry_.reactions(), i)
     {
         const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
 
         // for each reaction compute omegai
-        scalar omegai = R.omega(p, T, c1, li, pf, cf, lRef, pr, cr, rRef);
+        scalar omegaf, omegar;
+        const scalar omegai = R.omega(p, T, c1, li, omegaf, omegar);
 
         // Then for each pair of species composing this reaction,
         // compute the rAB matrix (separate the numerator and
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C
index a807df2ca2..35011e41ef 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C
@@ -86,14 +86,13 @@ void Foam::chemistryReductionMethods::DRG<ThermoType>::reduceMechanism
     // Index of the other species involved in the rABNum
     RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
 
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
     forAll(this->chemistry_.reactions(), i)
     {
         const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
 
         // For each reaction compute omegai
-        scalar omegai = R.omega(p, T, c1, li, pf, cf, lRef, pr, cr, rRef);
+        scalar omegaf, omegar;
+        const scalar omegai = R.omega(p, T, c1, li, omegaf, omegar);
 
         // Then for each pair of species composing this reaction,
         // compute the rAB matrix (separate the numerator and
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C
index ce765a8e48..c4c45efc78 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C
@@ -131,15 +131,14 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
     // Index of the other species involved in the rABNum
     RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
 
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
     scalarField omegaV(this->chemistry_.reactions().size());
     forAll(this->chemistry_.reactions(), i)
     {
         const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
 
         // for each reaction compute omegai
-        scalar omegai = R.omega(p, T, c1, li, pf, cf, lRef, pr, cr, rRef);
+        scalar omegaf, omegar;
+        const scalar omegai = R.omega(p, T, c1, li, omegaf, omegar);
         omegaV[i] = omegai;
 
         // then for each pair of species composing this reaction,
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.C
index a42532b61b..709364751a 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.C
@@ -127,16 +127,15 @@ void Foam::chemistryReductionMethods::EFA<ThermoType>::reduceMechanism
     // Index of the other species involved in the rABNum
     RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
 
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
     forAll(this->chemistry_.reactions(), i)
     {
         const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
 
         // for each reaction compute omegai
-        R.omega(p, T, c1, li, pf, cf, lRef, pr, cr, rRef);
+        scalar omegaf, omegar;
+        R.omega(p, T, c1, li, omegaf, omegar);
 
-        scalar fr = mag(pf*cf)+mag(pr*cr);
+        scalar fr = mag(omegaf) + mag(omegar);
         scalar NCi(0.0),NHi(0.0),NOi(0.0),NNi(0.0);
         forAll(R.lhs(),s)
         {
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.C
index 7375d0edbb..5c0a4be2dc 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.C
@@ -88,14 +88,13 @@ void Foam::chemistryReductionMethods::PFA<ThermoType>::reduceMechanism
     // Index of the other species involved in the rABNum
     RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
 
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
     forAll(this->chemistry_.reactions(), i)
     {
         const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
 
         // for each reaction compute omegai
-        scalar omegai = R.omega(p, T, c1, li, pf, cf, lRef, pr, cr, rRef);
+        scalar omegaf, omegar;
+        const scalar omegai = R.omega(p, T, c1, li, omegaf, omegar);
 
         // then for each pair of species composing this reaction,
         // compute the rAB matrix (separate the numerator and
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
index ac6bdf0a19..bcf0b373cd 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
@@ -295,13 +295,9 @@ void Foam::Reaction<ReactionThermo>::omega
     scalarField& dcdt
 ) const
 {
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
+    scalar omegaf, omegar;
 
-    scalar omegaI = omega
-    (
-        p, T, c, li, pf, cf, lRef, pr, cr, rRef
-    );
+    const scalar omegaI = omega(p, T, c, li, omegaf, omegar);
 
     forAll(lhs(), i)
     {
@@ -325,111 +321,29 @@ Foam::scalar Foam::Reaction<ReactionThermo>::omega
     const scalar T,
     const scalarField& c,
     const label li,
-    scalar& pf,
-    scalar& cf,
-    label& lRef,
-    scalar& pr,
-    scalar& cr,
-    label& rRef
+    scalar& omegaf,
+    scalar& omegar
 ) const
 {
-
     scalar clippedT = min(max(T, this->Tlow()), this->Thigh());
 
-    const scalar kf = this->kf(p, clippedT, c, li);
-    const scalar kr = this->kr(kf, p, clippedT, c, li);
+    omegaf = this->kf(p, clippedT, c, li);
+    omegar = this->kr(omegaf, p, clippedT, c, li);
 
-    pf = 1;
-    pr = 1;
-
-    const label Nl = lhs().size();
-    const label Nr = rhs().size();
-
-    label slRef = 0;
-    lRef = lhs()[slRef].index;
-
-    pf = kf;
-    for (label s = 1; s < Nl; s++)
+    forAll(lhs(), i)
     {
-        const label si = lhs()[s].index;
-
-        if (c[si] < c[lRef])
-        {
-            const scalar exp = lhs()[slRef].exponent;
-            pf *= pow(max(c[lRef], 0), exp);
-            lRef = si;
-            slRef = s;
-        }
-        else
-        {
-            const scalar exp = lhs()[s].exponent;
-            pf *= pow(max(c[si], 0), exp);
-        }
+        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;
     }
-    cf = max(c[lRef], 0);
-
+    forAll(rhs(), i)
     {
-        const scalar exp = lhs()[slRef].exponent;
-        if (exp < 1)
-        {
-            if (cf > small)
-            {
-                pf *= pow(cf, exp - 1);
-            }
-            else
-            {
-                pf = 0;
-            }
-        }
-        else
-        {
-            pf *= pow(cf, exp - 1);
-        }
+        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;
     }
 
-    label srRef = 0;
-    rRef = rhs()[srRef].index;
-
-    // Find the matrix element and element position for the rhs
-    pr = kr;
-    for (label s = 1; s < Nr; s++)
-    {
-        const label si = rhs()[s].index;
-        if (c[si] < c[rRef])
-        {
-            const scalar exp = rhs()[srRef].exponent;
-            pr *= pow(max(c[rRef], 0), exp);
-            rRef = si;
-            srRef = s;
-        }
-        else
-        {
-            const scalar exp = rhs()[s].exponent;
-            pr *= pow(max(c[si], 0), exp);
-        }
-    }
-    cr = max(c[rRef], 0);
-
-    {
-        const scalar exp = rhs()[srRef].exponent;
-        if (exp < 1)
-        {
-            if (cr > small)
-            {
-                pr *= pow(cr, exp - 1);
-            }
-            else
-            {
-                pr = 0;
-            }
-        }
-        else
-        {
-            pr *= pow(cr, exp - 1);
-        }
-    }
-
-    return pf*cf - pr*cr;
+    return omegaf - omegar;
 }
 
 
@@ -449,10 +363,9 @@ void Foam::Reaction<ReactionThermo>::dwdc
     const List<label>& c2s
 ) const
 {
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
+    scalar omegaf, omegar;
 
-    omegaI = omega(p, T, c, li, pf, cf, lRef, pr, cr, rRef);
+    omegaI = omega(p, T, c, li, omegaf, omegar);
 
     forAll(lhs(), i)
     {
@@ -480,25 +393,17 @@ void Foam::Reaction<ReactionThermo>::dwdc
             const scalar el = lhs()[i].exponent;
             if (i == j)
             {
-                if (el < 1)
-                {
-                    if (c[si] > SMALL)
-                    {
-                        kf *= el*pow(c[si] + VSMALL, el - 1);
-                    }
-                    else
-                    {
-                        kf = 0;
-                    }
-                }
-                else
-                {
-                    kf *= el*pow(c[si], el - 1);
-                }
+                kf *=
+                    c[si] >= small || el >= 1
+                  ? el*pow(max(c[si], 0), el - 1)
+                  : 0;
             }
             else
             {
-                kf *= pow(c[si], el);
+                kf *=
+                    c[si] >= small || el >= 1
+                  ? pow(max(c[si], 0), el)
+                  : 0;
             }
         }
 
@@ -526,25 +431,17 @@ void Foam::Reaction<ReactionThermo>::dwdc
             const scalar er = rhs()[i].exponent;
             if (i == j)
             {
-                if (er < 1)
-                {
-                    if (c[si] > SMALL)
-                    {
-                        kr *= er*pow(c[si] + VSMALL, er - 1);
-                    }
-                    else
-                    {
-                        kr = 0;
-                    }
-                }
-                else
-                {
-                    kr *= er*pow(c[si], er - 1);
-                }
+                kr *=
+                    c[si] >= small || er >= 1
+                  ? er*pow(max(c[si], 0), er - 1)
+                  : 0;
             }
             else
             {
-                kr *= pow(c[si], er);
+                kr *=
+                    c[si] >= small || er >= 1
+                  ? pow(max(c[si], 0), er)
+                  : 0;
             }
         }
 
@@ -620,13 +517,13 @@ void Foam::Reaction<ReactionThermo>::dwdT
     {
         const label si = lhs()[i].index;
         const scalar el = lhs()[i].exponent;
-        dkfdT *= pow(c[si], el);
+        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 *= pow(c[si], er);
+        dkrdT *= c[si] >= small || er >= 1 ? pow(max(c[si], 0), er) : 0;
     }
 
     const scalar dqidT = dkfdT - dkrdT;
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
index a778592cae..734b55b18e 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
@@ -241,12 +241,8 @@ public:
                 const scalar T,
                 const scalarField& c,
                 const label li,
-                scalar& pf,
-                scalar& cf,
-                label& lRef,
-                scalar& pr,
-                scalar& cr,
-                label& rRef
+                scalar& omegaf,
+                scalar& omegar
             ) const;
 
         // Reaction rate coefficients