From 4dc35c6810d2bbebbeb4111318ea5bb56afecffd Mon Sep 17 00:00:00 2001
From: Henry Weller <http://openfoam.org>
Date: Fri, 15 Jun 2018 12:26:59 +0100
Subject: [PATCH] thermophysicalModels: Implementation of the full algebraic
 Jacobian

including third-body and pressure dependent derivatives, and derivative of the
temperature term.  The complete Jacobian is more robust than the incomplete and
partially approximate form used previously and improves the efficiency of the
stiff ODE solvers which rely on the Jacobian.

Reaction rate evaluation moved from the chemistryModel to specie library to
simplfy support for alternative reaction rate expressions and associated
Jacobian terms.

Temperature clipping included in the Reaction class. This is inactive by default
but for most cases it is advised to provide temperature limits (high and
low). These are provided in the foamChemistryFile with the keywords Thigh and
Tlow. When using chemkinToFoam these values are set to the limits of the Janaf
thermodynamic data.  With the new Jacobian this temperature clipping has proved
very beneficial for stability and for some cases essential.

Improvement of the TDAC MRU list better integrated in add and grow functions.

To get the most out of this significant development it is important to re-tune
the ODE integration tolerances, in particular the absTol in the odeCoeffs
sub-dictionary of the chemistryProperties dictionary:

odeCoeffs
{
    solver          seulex;
    absTol          1e-12;
    relTol          0.01;
}

Typically absTol can now be set to 1e-8 and relTol to 0.1 except for ignition
time problems, and with theses settings the integration is still robust but for
many cases a lot faster than previously.

Code development and integration undertaken by
Francesco Contino
Henry G. Weller, CFD Direct
---
 .../chemkinToFoam/chemkinToFoam.C             |  18 +-
 .../StandardChemistryModel.C                  | 319 +++-----------
 .../StandardChemistryModel.H                  |  20 +-
 .../TDACChemistryModel/TDACChemistryModel.C   | 271 +++++-------
 .../TDACChemistryModel/TDACChemistryModel.H   |  12 +-
 .../TDACChemistryModel/tabulation/ISAT/ISAT.C |  84 ++--
 .../radiationCoupledBase.C                    |  21 +-
 .../chemkinReader/chemkinLexer.L              |  12 +
 .../chemkinReader/chemkinReader.C             |  20 +-
 .../Reactions/solidReaction/solidReaction.H   |  31 +-
 .../solidArrheniusReactionRate.H              |  33 +-
 .../solidArrheniusReactionRateI.H             |  52 ++-
 .../IrreversibleReaction.C                    | 156 ++++++-
 .../IrreversibleReaction.H                    |  70 ++-
 .../NonEquilibriumReversibleReaction.C        | 113 ++++-
 .../NonEquilibriumReversibleReaction.H        |  47 +-
 .../reaction/Reactions/Reaction/Reaction.C    | 417 +++++++++++++++++-
 .../reaction/Reactions/Reaction/Reaction.H    | 161 ++++++-
 .../reaction/Reactions/Reaction/ReactionI.H   |  21 +
 .../Reactions/ReactionList/ReactionList.C     |   9 +-
 .../Reactions/ReactionProxy/ReactionProxy.C   | 207 +++++++++
 .../Reactions/ReactionProxy/ReactionProxy.H   | 189 ++++++++
 .../ReversibleReaction/ReversibleReaction.C   | 111 +++++
 .../ReversibleReaction/ReversibleReaction.H   |  47 +-
 .../ArrheniusReactionRate.H                   |  33 +-
 .../ArrheniusReactionRateI.H                  |  51 +++
 .../ChemicallyActivatedReactionRate.H         |  37 +-
 .../ChemicallyActivatedReactionRateI.H        | 124 +++++-
 .../FallOffReactionRate/FallOffReactionRate.H |  35 +-
 .../FallOffReactionRateI.H                    | 113 ++++-
 .../JanevReactionRate/JanevReactionRate.H     |  31 +-
 .../JanevReactionRate/JanevReactionRateI.H    |  71 ++-
 .../LandauTellerReactionRate.H                |  31 +-
 .../LandauTellerReactionRateI.H               |  77 +++-
 .../LangmuirHinshelwoodReactionRate.H         |  31 +-
 .../LangmuirHinshelwoodReactionRateI.H        |  73 ++-
 .../LindemannFallOffFunction.H                |  18 +-
 .../LindemannFallOffFunctionI.H               |  26 +-
 .../SRIFallOffFunction/SRIFallOffFunction.H   |  18 +-
 .../SRIFallOffFunction/SRIFallOffFunctionI.H  |  39 +-
 .../TroeFallOffFunction/TroeFallOffFunction.H |  18 +-
 .../TroeFallOffFunctionI.H                    |  97 +++-
 .../infiniteReactionRate.H                    |  31 +-
 .../infiniteReactionRateI.H                   |  40 +-
 .../powerSeries/powerSeriesReactionRate.H     |  31 +-
 .../powerSeries/powerSeriesReactionRateI.H    |  60 ++-
 .../thirdBodyArrheniusReactionRate.H          |  33 ++
 .../thirdBodyArrheniusReactionRateI.H         |  63 ++-
 .../thirdBodyEfficienciesI.H                  |   4 +-
 .../specie/thermo/eConst/eConstThermo.H       |  12 +-
 .../specie/thermo/eConst/eConstThermoI.H      |  22 +
 .../specie/thermo/hConst/hConstThermo.H       |  11 +-
 .../specie/thermo/hConst/hConstThermoI.H      |  20 +
 .../thermo/hPolynomial/hPolynomialThermo.H    |  11 +-
 .../thermo/hPolynomial/hPolynomialThermoI.H   |  30 ++
 .../specie/thermo/hPower/hPowerThermo.H       |  11 +-
 .../specie/thermo/hPower/hPowerThermoI.H      |  24 +
 .../specie/thermo/hRefConst/hRefConstThermo.H |  11 +-
 .../thermo/hRefConst/hRefConstThermoI.H       |  20 +
 .../specie/thermo/janaf/janafThermo.H         |  11 +-
 .../specie/thermo/janaf/janafThermoI.H        |  25 ++
 .../specie/thermo/thermo/thermo.H             |  13 +-
 .../specie/thermo/thermo/thermoI.H            |  29 ++
 .../nc7h16/constant/chemistryProperties       |   2 +-
 .../constant/panelRegion/chemistryProperties  |   2 +-
 .../DLR_A_LTS/constant/chemistryProperties    |  19 +-
 .../constant/chemistryProperties.test         | 142 ------
 .../RAS/DLR_A_LTS/constant/reactionsGRI       |   3 +
 .../reactingFoam/RAS/SandiaD_LTS/Allrun       |   4 +-
 .../SandiaD_LTS/constant/chemistryProperties  |   6 +-
 .../RAS/SandiaD_LTS/constant/reactionsGRI     |   3 +
 .../RAS/SandiaD_LTS/system/controlDict        |   8 +-
 .../constant/chemistryProperties              |   2 +-
 .../constant/chemistryProperties              |   2 +-
 .../constant/chemistryProperties              |   7 +-
 .../constant/reactionsGRI                     |   3 +
 .../system/controlDict                        |   2 +-
 .../constant/chemistryProperties              |  21 +-
 .../constant/reactionsGRI                     |  85 ++--
 .../constant/chemistryProperties              |   8 +-
 .../constant/chemistryProperties.new          |   2 +-
 .../constant/reactionsGRI                     |  85 ++--
 .../constant/chemistryProperties              |   2 +-
 .../constant/chemistryProperties.gas          |   2 +-
 84 files changed, 3380 insertions(+), 906 deletions(-)
 create mode 100644 src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C
 create mode 100644 src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H
 delete mode 100644 tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/chemistryProperties.test

diff --git a/applications/utilities/thermophysical/chemkinToFoam/chemkinToFoam.C b/applications/utilities/thermophysical/chemkinToFoam/chemkinToFoam.C
index cb219fb83..e094ee28a 100644
--- a/applications/utilities/thermophysical/chemkinToFoam/chemkinToFoam.C
+++ b/applications/utilities/thermophysical/chemkinToFoam/chemkinToFoam.C
@@ -67,15 +67,13 @@ int main(int argc, char *argv[])
 
     chemkinReader cr(species, args[1], args[3], args[2], newFormat);
 
-
     OFstream reactionsFile(args[4]);
-    reactionsFile
-        << "elements" << cr.elementNames() << token::END_STATEMENT << nl << nl;
-    reactionsFile
-        << "species" << cr.species() << token::END_STATEMENT << nl << nl;
+    reactionsFile.writeKeyword("elements")
+        << cr.elementNames() << token::END_STATEMENT << nl << nl;
+    reactionsFile.writeKeyword("species")
+        << cr.species() << token::END_STATEMENT << nl << nl;
     cr.reactions().write(reactionsFile);
 
-
     // Temporary hack to splice the specie composition data into the thermo file
     // pending complete integration into the thermodynamics structure
 
@@ -103,6 +101,14 @@ int main(int argc, char *argv[])
 
     thermoDict.write(OFstream(args[5])(), false);
 
+    reactionsFile << nl;
+
+    reactionsFile.writeKeyword("Tlow")
+        << Reaction<gasHThermoPhysics>::TlowDefault
+        << token::END_STATEMENT << nl;
+    reactionsFile.writeKeyword("Thigh")
+        << Reaction<gasHThermoPhysics>::ThighDefault
+        << token::END_STATEMENT << nl << nl;
 
     Info<< "End\n" << endl;
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
index ad2c75577..adc979f62 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
@@ -56,7 +56,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::StandardChemistryModel
         BasicChemistryModel<ReactionThermo>::template lookupOrDefault<scalar>
         (
             "Treact",
-            0.0
+            0
         )
     ),
     RR_(nSpecie_),
@@ -80,7 +80,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::StandardChemistryModel
                     IOobject::NO_WRITE
                 ),
                 thermo.p().mesh(),
-                dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+                dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0)
             )
         );
     }
@@ -109,8 +109,6 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omega
     scalarField& dcdt
 ) const
 {
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
 
     dcdt = Zero;
 
@@ -118,24 +116,7 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omega
     {
         const Reaction<ThermoType>& R = reactions_[i];
 
-        const scalar omegai = omega
-        (
-            R, c, T, p, pf, cf, lRef, pr, cr, rRef
-        );
-
-        forAll(R.lhs(), s)
-        {
-            const label si = R.lhs()[s].index;
-            const scalar sl = R.lhs()[s].stoichCoeff;
-            dcdt[si] -= sl*omegai;
-        }
-
-        forAll(R.rhs(), s)
-        {
-            const label si = R.rhs()[s].index;
-            const scalar sr = R.rhs()[s].stoichCoeff;
-            dcdt[si] += sr*omegai;
-        }
+        R.omega(p, T, c, dcdt);
     }
 }
 
@@ -156,123 +137,11 @@ Foam::scalar Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omegaI
 ) const
 {
     const Reaction<ThermoType>& R = reactions_[index];
-    scalar w = omega(R, c, T, p, pf, cf, lRef, pr, cr, rRef);
+    scalar w = R.omega(p, T, c, pf, cf, lRef, pr, cr, rRef);
     return(w);
 }
 
 
-template<class ReactionThermo, class ThermoType>
-Foam::scalar Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omega
-(
-    const Reaction<ThermoType>& R,
-    const scalarField& c,
-    const scalar T,
-    const scalar p,
-    scalar& pf,
-    scalar& cf,
-    label& lRef,
-    scalar& pr,
-    scalar& cr,
-    label& rRef
-) const
-{
-    const scalar kf = R.kf(p, T, c);
-    const scalar kr = R.kr(kf, p, T, c);
-
-    pf = 1.0;
-    pr = 1.0;
-
-    const label Nl = R.lhs().size();
-    const label Nr = R.rhs().size();
-
-    label slRef = 0;
-    lRef = R.lhs()[slRef].index;
-
-    pf = kf;
-    for (label s = 1; s < Nl; s++)
-    {
-        const label si = R.lhs()[s].index;
-
-        if (c[si] < c[lRef])
-        {
-            const scalar exp = R.lhs()[slRef].exponent;
-            pf *= pow(max(c[lRef], 0.0), exp);
-            lRef = si;
-            slRef = s;
-        }
-        else
-        {
-            const scalar exp = R.lhs()[s].exponent;
-            pf *= pow(max(c[si], 0.0), exp);
-        }
-    }
-    cf = max(c[lRef], 0.0);
-
-    {
-        const scalar exp = R.lhs()[slRef].exponent;
-        if (exp < 1.0)
-        {
-            if (cf > small)
-            {
-                pf *= pow(cf, exp - 1.0);
-            }
-            else
-            {
-                pf = 0.0;
-            }
-        }
-        else
-        {
-            pf *= pow(cf, exp - 1.0);
-        }
-    }
-
-    label srRef = 0;
-    rRef = R.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 = R.rhs()[s].index;
-        if (c[si] < c[rRef])
-        {
-            const scalar exp = R.rhs()[srRef].exponent;
-            pr *= pow(max(c[rRef], 0.0), exp);
-            rRef = si;
-            srRef = s;
-        }
-        else
-        {
-            const scalar exp = R.rhs()[s].exponent;
-            pr *= pow(max(c[si], 0.0), exp);
-        }
-    }
-    cr = max(c[rRef], 0.0);
-
-    {
-        const scalar exp = R.rhs()[srRef].exponent;
-        if (exp < 1.0)
-        {
-            if (cr > small)
-            {
-                pr *= pow(cr, exp - 1.0);
-            }
-            else
-            {
-                pr = 0.0;
-            }
-        }
-        else
-        {
-            pr *= pow(cr, exp - 1.0);
-        }
-    }
-
-    return pf*cf - pr*cr;
-}
-
-
 template<class ReactionThermo, class ThermoType>
 void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
 (
@@ -286,29 +155,29 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
 
     forAll(c_, i)
     {
-        c_[i] = max(c[i], 0.0);
+        c_[i] = max(c[i], 0);
     }
 
     omega(c_, T, p, dcdt);
 
     // Constant pressure
     // dT/dt = ...
-    scalar rho = 0.0;
-    scalar cSum = 0.0;
+    scalar rho = 0;
+    scalar cSum = 0;
     for (label i = 0; i < nSpecie_; i++)
     {
         const scalar W = specieThermo_[i].W();
         cSum += c_[i];
         rho += W*c_[i];
     }
-    scalar cp = 0.0;
+    scalar cp = 0;
     for (label i=0; i<nSpecie_; i++)
     {
         cp += c_[i]*specieThermo_[i].cp(p, T);
     }
     cp /= rho;
 
-    scalar dT = 0.0;
+    scalar dT = 0;
     for (label i = 0; i < nSpecie_; i++)
     {
         const scalar hi = specieThermo_[i].ha(p, T);
@@ -319,7 +188,7 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
     dcdt[nSpecie_] = -dT;
 
     // dp/dt = ...
-    dcdt[nSpecie_ + 1] = 0.0;
+    dcdt[nSpecie_ + 1] = 0;
 }
 
 
@@ -329,7 +198,7 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::jacobian
     const scalar t,
     const scalarField& c,
     scalarField& dcdt,
-    scalarSquareMatrix& dfdc
+    scalarSquareMatrix& J
 ) const
 {
     const scalar T = c[nSpecie_];
@@ -337,131 +206,67 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::jacobian
 
     forAll(c_, i)
     {
-        c_[i] = max(c[i], 0.0);
+        c_[i] = max(c[i], 0);
     }
 
-    dfdc = Zero;
-
-    // Length of the first argument must be nSpecie_
-    omega(c_, T, p, dcdt);
+    J = Zero;
+    dcdt = Zero;
 
+    // To compute the species derivatives of the temperature term,
+    // the enthalpies of the individual species is needed
+    scalarField hi(nSpecie_);
+    scalarField cpi(nSpecie_);
+    for (label i = 0; i < nSpecie_; i++)
+    {
+        hi[i] = specieThermo_[i].ha(p, T);
+        cpi[i] = specieThermo_[i].cp(p, T);
+    }
+    scalar omegaI = 0;
+    List<label> dummy;
     forAll(reactions_, ri)
     {
         const Reaction<ThermoType>& R = reactions_[ri];
-
-        const scalar kf0 = R.kf(p, T, c_);
-        const scalar kr0 = R.kr(kf0, p, T, c_);
-
-        forAll(R.lhs(), j)
-        {
-            const label sj = R.lhs()[j].index;
-            scalar kf = kf0;
-            forAll(R.lhs(), i)
-            {
-                const label si = R.lhs()[i].index;
-                const scalar el = R.lhs()[i].exponent;
-                if (i == j)
-                {
-                    if (el < 1.0)
-                    {
-                        if (c_[si] > small)
-                        {
-                            kf *= el*pow(c_[si], el - 1.0);
-                        }
-                        else
-                        {
-                            kf = 0.0;
-                        }
-                    }
-                    else
-                    {
-                        kf *= el*pow(c_[si], el - 1.0);
-                    }
-                }
-                else
-                {
-                    kf *= pow(c_[si], el);
-                }
-            }
-
-            forAll(R.lhs(), i)
-            {
-                const label si = R.lhs()[i].index;
-                const scalar sl = R.lhs()[i].stoichCoeff;
-                dfdc(si, sj) -= sl*kf;
-            }
-            forAll(R.rhs(), i)
-            {
-                const label si = R.rhs()[i].index;
-                const scalar sr = R.rhs()[i].stoichCoeff;
-                dfdc(si, sj) += sr*kf;
-            }
-        }
-
-        forAll(R.rhs(), j)
-        {
-            const label sj = R.rhs()[j].index;
-            scalar kr = kr0;
-            forAll(R.rhs(), i)
-            {
-                const label si = R.rhs()[i].index;
-                const scalar er = R.rhs()[i].exponent;
-                if (i == j)
-                {
-                    if (er < 1.0)
-                    {
-                        if (c_[si] > small)
-                        {
-                            kr *= er*pow(c_[si], er - 1.0);
-                        }
-                        else
-                        {
-                            kr = 0.0;
-                        }
-                    }
-                    else
-                    {
-                        kr *= er*pow(c_[si], er - 1.0);
-                    }
-                }
-                else
-                {
-                    kr *= pow(c_[si], er);
-                }
-            }
-
-            forAll(R.lhs(), i)
-            {
-                const label si = R.lhs()[i].index;
-                const scalar sl = R.lhs()[i].stoichCoeff;
-                dfdc(si, sj) += sl*kr;
-            }
-            forAll(R.rhs(), i)
-            {
-                const label si = R.rhs()[i].index;
-                const scalar sr = R.rhs()[i].stoichCoeff;
-                dfdc(si, sj) -= sr*kr;
-            }
-        }
+        scalar kfwd, kbwd;
+        R.dwdc(p, T, c_, J, dcdt, omegaI, kfwd, kbwd, false, dummy);
+        R.dwdT(p, T, c_, omegaI, kfwd, kbwd, J, false, dummy, nSpecie_);
     }
 
-    // Calculate the dcdT elements numerically
-    const scalar delta = 1.0e-3;
-
-    omega(c_, T + delta, p, dcdt_);
+    // The species derivatives of the temperature term are partially computed
+    // while computing dwdc, they are completed hereunder:
+    scalar cpMean = 0;
+    scalar dcpdTMean = 0;
     for (label i=0; i<nSpecie_; i++)
     {
-        dfdc(i, nSpecie_) = dcdt_[i];
+        cpMean += c_[i]*cpi[i]; // J/(m3.K)
+        dcpdTMean += c_[i]*specieThermo_[i].dcpdT(p, T);
     }
-
-    omega(c_, T - delta, p, dcdt_);
+    scalar dTdt = 0.0;
     for (label i=0; i<nSpecie_; i++)
     {
-        dfdc(i, nSpecie_) = 0.5*(dfdc(i, nSpecie_) - dcdt_[i])/delta;
+        dTdt += hi[i]*dcdt[i]; // J/(m3.s)
+    }
+    dTdt /= -cpMean; // K/s
+
+    for (label i = 0; i < nSpecie_; i++)
+    {
+        J(nSpecie_, i) = 0;
+        for (label j = 0; j < nSpecie_; j++)
+        {
+            J(nSpecie_, i) += hi[j]*J(j, i);
+        }
+        J(nSpecie_, i) += cpi[i]*dTdt; // J/(mol.s)
+        J(nSpecie_, i) /= -cpMean;    // K/s / (mol/m3)
     }
 
-    dfdc(nSpecie_, nSpecie_) = 0;
-    dfdc(nSpecie_ + 1, nSpecie_) = 0;
+    // ddT of dTdt
+    J(nSpecie_, nSpecie_) = 0;
+    for (label i = 0; i < nSpecie_; i++)
+    {
+        J(nSpecie_, nSpecie_) += cpi[i]*dcdt[i] + hi[i]*J(i, nSpecie_);
+    }
+    J(nSpecie_, nSpecie_) += dTdt*dcpdTMean;
+    J(nSpecie_, nSpecie_) /= -cpMean;
+    J(nSpecie_, nSpecie_) += dTdt/T;
 }
 
 
@@ -509,7 +314,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::tc() const
             const scalar Ti = T[celli];
             const scalar pi = p[celli];
 
-            scalar cSum = 0.0;
+            scalar cSum = 0;
 
             for (label i=0; i<nSpecie_; i++)
             {
@@ -521,7 +326,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::tc() const
             {
                 const Reaction<ThermoType>& R = reactions_[i];
 
-                omega(R, c_, Ti, pi, pf, cf, lRef, pr, cr, rRef);
+                R.omega(pi, Ti, c_, pf, cf, lRef, pr, cr, rRef);
 
                 forAll(R.rhs(), s)
                 {
@@ -557,7 +362,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::Qdot() const
                 false
             ),
             this->mesh_,
-            dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0.0)
+            dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0)
         )
     );
 
@@ -600,7 +405,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::calculateRR
                 IOobject::NO_WRITE
             ),
             this->mesh(),
-            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0)
         )
     );
 
@@ -628,7 +433,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::calculateRR
         }
 
         const Reaction<ThermoType>& R = reactions_[ri];
-        const scalar omegai = omega(R, c_, Ti, pi, pf, cf, lRef, pr, cr, rRef);
+        const scalar omegai = R.omega(pi, Ti, c_, pf, cf, lRef, pr, cr, rRef);
 
         forAll(R.lhs(), s)
         {
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H
index 1a75627a0..c7f1b5005 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -163,22 +163,6 @@ public:
             scalarField& dcdt
         ) const;
 
-        //- Return the reaction rate for reaction r and the reference
-        //  species and charateristic times
-        virtual scalar omega
-        (
-            const Reaction<ThermoType>& r,
-            const scalarField& c,
-            const scalar T,
-            const scalar p,
-            scalar& pf,
-            scalar& cf,
-            label& lRef,
-            scalar& pr,
-            scalar& cr,
-            label& rRef
-        ) const;
-
 
         //- Return the reaction rate for iReaction and the reference
         //  species and charateristic times
@@ -254,7 +238,7 @@ public:
                 const scalar t,
                 const scalarField& c,
                 scalarField& dcdt,
-                scalarSquareMatrix& dfdc
+                scalarSquareMatrix& J
             ) const;
 
             virtual void solve
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
index 7caf8bf55..777fa225d 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
@@ -168,9 +168,9 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::omega
         {
             const Reaction<ThermoType>& R = this->reactions_[i];
 
-            scalar omegai = omega
+            scalar omegai = R.omega
             (
-                R, c, T, p, pf, cf, lRef, pr, cr, rRef
+                p, T, c, pf, cf, lRef, pr, cr, rRef
             );
 
             forAll(R.lhs(), s)
@@ -232,17 +232,17 @@ Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::omega
         if (c[si] < c[lRef])
         {
             const scalar exp = R.lhs()[slRef].exponent;
-            pf *= pow(max(c[lRef], 0.0), exp);
+            pf *= pow(max(c[lRef], 0), exp);
             lRef = si;
             slRef = s;
         }
         else
         {
             const scalar exp = R.lhs()[s].exponent;
-            pf *= pow(max(c[si], 0.0), exp);
+            pf *= pow(max(c[si], 0), exp);
         }
     }
-    cf = max(c[lRef], 0.0);
+    cf = max(c[lRef], 0);
 
     {
         const scalar exp = R.lhs()[slRef].exponent;
@@ -274,17 +274,17 @@ Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::omega
         if (c[si] < c[rRef])
         {
             const scalar exp = R.rhs()[srRef].exponent;
-            pr *= pow(max(c[rRef], 0.0), exp);
+            pr *= pow(max(c[rRef], 0), exp);
             rRef = si;
             srRef = s;
         }
         else
         {
             const scalar exp = R.rhs()[s].exponent;
-            pr *= pow(max(c[si], 0.0), exp);
+            pr *= pow(max(c[si], 0), exp);
         }
     }
-    cr = max(c[rRef], 0.0);
+    cr = max(c[rRef], 0);
 
     {
         const scalar exp = R.rhs()[srRef].exponent;
@@ -334,14 +334,14 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::derivatives
         // efficiencies
         for (label i=0; i<NsDAC_; i++)
         {
-            this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0.0);
+            this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0);
         }
     }
     else
     {
         for (label i=0; i<this->nSpecie(); i++)
         {
-            this->c_[i] = max(c[i], 0.0);
+            this->c_[i] = max(c[i], 0);
         }
     }
 
@@ -398,7 +398,8 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::jacobian
 (
     const scalar t,
     const scalarField& c,
-    scalarSquareMatrix& dfdc
+    scalarField& dcdt,
+    scalarSquareMatrix& J
 ) const
 {
     const bool reduced = mechRed_->active();
@@ -416,183 +417,119 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::jacobian
         this->c_ = completeC_;
         for (label i=0; i<NsDAC_; i++)
         {
-            this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0.0);
+            this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0);
         }
     }
     else
     {
         forAll(this->c_, i)
         {
-            this->c_[i] = max(c[i], 0.0);
+            this->c_[i] = max(c[i], 0);
         }
     }
 
-    dfdc = Zero;
+    J = Zero;
+    dcdt = Zero;
+    scalarField hi(this->c_.size());
+    scalarField cpi(this->c_.size());
+    forAll(hi, i)
+    {
+        hi[i] = this->specieThermo_[i].ha(p, T);
+        cpi[i] = this->specieThermo_[i].cp(p, T);
+    }
+
+    scalar omegaI = 0;
 
     forAll(this->reactions_, ri)
     {
         if (!reactionsDisabled_[ri])
         {
             const Reaction<ThermoType>& R = this->reactions_[ri];
-
-            const scalar kf0 = R.kf(p, T, this->c_);
-            const scalar kr0 = R.kr(kf0, p, T, this->c_);
-
-            forAll(R.lhs(), j)
-            {
-                label sj = R.lhs()[j].index;
-                if (reduced)
-                {
-                    sj = completeToSimplifiedIndex_[sj];
-                }
-                scalar kf = kf0;
-                forAll(R.lhs(), i)
-                {
-                    const label si = R.lhs()[i].index;
-                    const scalar el = R.lhs()[i].exponent;
-                    if (i == j)
-                    {
-                        if (el < 1)
-                        {
-                            if (this->c_[si] > small)
-                            {
-                                kf *= el*pow(this->c_[si], el - 1);
-                            }
-                            else
-                            {
-                                kf = 0;
-                            }
-                        }
-                        else
-                        {
-                            kf *= el*pow(this->c_[si], el - 1);
-                        }
-                    }
-                    else
-                    {
-                        kf *= pow(this->c_[si], el);
-                    }
-                }
-
-                forAll(R.lhs(), i)
-                {
-                    label si = R.lhs()[i].index;
-                    if (reduced)
-                    {
-                        si = completeToSimplifiedIndex_[si];
-                    }
-                    const scalar sl = R.lhs()[i].stoichCoeff;
-                    dfdc(si, sj) -= sl*kf;
-                }
-                forAll(R.rhs(), i)
-                {
-                    label si = R.rhs()[i].index;
-                    if (reduced)
-                    {
-                        si = completeToSimplifiedIndex_[si];
-                    }
-                    const scalar sr = R.rhs()[i].stoichCoeff;
-                    dfdc(si, sj) += sr*kf;
-                }
-            }
-
-            forAll(R.rhs(), j)
-            {
-                label sj = R.rhs()[j].index;
-                if (reduced)
-                {
-                    sj = completeToSimplifiedIndex_[sj];
-                }
-                scalar kr = kr0;
-                forAll(R.rhs(), i)
-                {
-                    const label si = R.rhs()[i].index;
-                    const scalar er = R.rhs()[i].exponent;
-                    if (i == j)
-                    {
-                        if (er < 1)
-                        {
-                            if (this->c_[si] > small)
-                            {
-                                kr *= er*pow(this->c_[si], er - 1);
-                            }
-                            else
-                            {
-                                kr = 0;
-                            }
-                        }
-                        else
-                        {
-                            kr *= er*pow(this->c_[si], er - 1);
-                        }
-                    }
-                    else
-                    {
-                        kr *= pow(this->c_[si], er);
-                    }
-                }
-
-                forAll(R.lhs(), i)
-                {
-                    label si = R.lhs()[i].index;
-                    if (reduced)
-                    {
-                        si = completeToSimplifiedIndex_[si];
-                    }
-                    const scalar sl = R.lhs()[i].stoichCoeff;
-                    dfdc(si, sj) += sl*kr;
-                }
-                forAll(R.rhs(), i)
-                {
-                    label si = R.rhs()[i].index;
-                    if (reduced)
-                    {
-                        si = completeToSimplifiedIndex_[si];
-                    }
-                    const scalar sr = R.rhs()[i].stoichCoeff;
-                    dfdc(si, sj) -= sr*kr;
-                }
-            }
+            scalar kfwd, kbwd;
+            R.dwdc
+            (
+                p,
+                T,
+                this->c_,
+                J,
+                dcdt,
+                omegaI,
+                kfwd,
+                kbwd,
+                reduced,
+                completeToSimplifiedIndex_
+            );
+            R.dwdT
+            (
+                p,
+                T,
+                this->c_,
+                omegaI,
+                kfwd,
+                kbwd,
+                J,
+                reduced,
+                completeToSimplifiedIndex_,
+                this->nSpecie_
+            );
         }
     }
 
-    // Calculate the dcdT elements numerically
-    const scalar delta = 1e-3;
-
-    omega(this->c_, T + delta, p, this->dcdt_);
-    for (label i=0; i<this->nSpecie_; i++)
+    // The species derivatives of the temperature term are partially computed
+    // while computing dwdc, they are completed hereunder:
+    scalar cpMean = 0;
+    scalar dcpdTMean = 0;
+    forAll(this->c_, i)
     {
-        dfdc(i, this->nSpecie_) = this->dcdt_[i];
+        cpMean += this->c_[i]*cpi[i]; // J/(m3.K)
+        // Already multiplied by rho
+        dcpdTMean += this->c_[i]*this->specieThermo_[i].dcpdT(p, T);
     }
 
-    omega(this->c_, T - delta, p, this->dcdt_);
-    for (label i=0; i<this->nSpecie_; i++)
+    scalar dTdt = 0;
+    forAll(hi, i)
     {
-        dfdc(i, this->nSpecie_) =
-            0.5*(dfdc(i, this->nSpecie_) - this->dcdt_[i])/delta;
+        if (reduced)
+        {
+            const label si = completeToSimplifiedIndex_[i];
+            if (si != -1)
+            {
+                dTdt += hi[i]*dcdt[si]; // J/(m3.s)
+            }
+        }
+        else
+        {
+            dTdt += hi[i]*dcdt[i]; // J/(m3.s)
+        }
+    }
+    dTdt /= -cpMean; // K/s
+    dcdt[this->nSpecie_] = dTdt;
+
+    for (label i = 0; i < this->nSpecie_; i++)
+    {
+        J(this->nSpecie_, i) = 0;
+        for (label j = 0; j < this->nSpecie_; j++)
+        {
+            const label sj = reduced ? simplifiedToCompleteIndex_[j] : j;
+            J(this->nSpecie_, i) += hi[sj]*J(j, i);
+        }
+        const label si = reduced ? simplifiedToCompleteIndex_[i] : i;
+        J(this->nSpecie_, i) += cpi[si]*dTdt; // J/(mol.s)
+        J(this->nSpecie_, i) /= -cpMean;    // K/s / (mol/m3)
     }
 
-    dfdc(this->nSpecie_, this->nSpecie_) = 0;
-    dfdc(this->nSpecie_ + 1, this->nSpecie_) = 0;
-}
-
-
-template<class ReactionThermo, class ThermoType>
-void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::jacobian
-(
-    const scalar t,
-    const scalarField& c,
-    scalarField& dcdt,
-    scalarSquareMatrix& dfdc
-) const
-{
-    jacobian(t, c, dfdc);
-
-    const scalar T = c[this->nSpecie_];
-    const scalar p = c[this->nSpecie_ + 1];
-
-    // Note: Uses the c_ field initialized by the call to jacobian above
-    omega(this->c_, T, p, dcdt);
+    // ddT of dTdt
+    J(this->nSpecie_, this->nSpecie_) = 0;
+    for (label i = 0; i < this->nSpecie_; i++)
+    {
+        const label si = reduced ? simplifiedToCompleteIndex_[i] : i;
+        J(this->nSpecie_, this->nSpecie_) +=
+            cpi[si]*dcdt[i]
+          + hi[si]*J(i, this->nSpecie_);
+    }
+    J(this->nSpecie_, this->nSpecie_) += dTdt*dcpdTMean;
+    J(this->nSpecie_, this->nSpecie_) /= -cpMean;
+    J(this->nSpecie_, this->nSpecie_) += dTdt/T;
 }
 
 
@@ -915,7 +852,7 @@ setTabulationResultsAdd
     const label celli
 )
 {
-    tabulationResults_[celli] = 0.0;
+    tabulationResults_[celli] = 0;
 }
 
 
@@ -923,7 +860,7 @@ template<class ReactionThermo, class ThermoType>
 void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::
 setTabulationResultsGrow(const label celli)
 {
-    tabulationResults_[celli] = 1.0;
+    tabulationResults_[celli] = 1;
 }
 
 
@@ -934,7 +871,7 @@ setTabulationResultsRetrieve
     const label celli
 )
 {
-    tabulationResults_[celli] = 2.0;
+    tabulationResults_[celli] = 2;
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H
index 87cc80eb6..f2cc41a68 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -221,20 +221,12 @@ public:
                 scalarField& dcdt
             ) const;
 
-            //- Pure jacobian function for tabulation
-            void jacobian
-            (
-                const scalar t,
-                const scalarField& c,
-                scalarSquareMatrix& dfdc
-            ) const;
-
             virtual void jacobian
             (
                 const scalar t,
                 const scalarField& c,
                 scalarField& dcdt,
-                scalarSquareMatrix& dfdc
+                scalarSquareMatrix& J
             ) const;
 
             virtual void solve
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C
index f8d957707..9b3d86233 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C
@@ -65,7 +65,7 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::ISAT
     (
         this->coeffsDict_.lookupOrDefault
         (
-            "minBalanceThreshold",0.1*chemisTree_.maxNLeafs()
+            "minBalanceThreshold", 0.1*chemisTree_.maxNLeafs()
         )
     ),
     MRURetrieve_(this->coeffsDict_.lookupOrDefault("MRURetrieve", false)),
@@ -169,17 +169,8 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::addToMRU
         {
             if (MRUList_.size() == maxMRUSize_)
             {
-                if (iter() == MRUList_.last())
-                {
-                    MRUList_.remove(iter);
-                    MRUList_.insert(phi0);
-                }
-                else
-                {
-                    FatalErrorInFunction
-                        << "wrong MRUList construction"
-                        << exit(FatalError);
-                }
+                MRUList_.remove(iter);
+                MRUList_.insert(phi0);
             }
             else
             {
@@ -238,13 +229,13 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::calcNewC
 
                 // As we use an approximation of A, Rphiq should be checked for
                 // negative values
-                Rphiq[i] = max(0.0,Rphiq[i]);
+                Rphiq[i] = max(0, Rphiq[i]);
             }
             // The species is not active A(i, j) = I(i, j)
             else
             {
                 Rphiq[i] += dphi[i];
-                Rphiq[i] = max(0.0,Rphiq[i]);
+                Rphiq[i] = max(0, Rphiq[i]);
             }
         }
         else // Mechanism reduction is not active
@@ -255,7 +246,7 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::calcNewC
             }
             // As we use a first order gradient matrix, Rphiq should be checked
             // for negative values
-            Rphiq[i] = max(0.0,Rphiq[i]);
+            Rphiq[i] = max(0, Rphiq[i]);
         }
     }
 }
@@ -322,6 +313,9 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::cleanAndBalance()
         }
         x = xtmp;
     }
+
+    MRUList_.clear();
+
     // Check if the tree should be balanced according to criterion:
     //  -the depth of the tree bigger than a*log2(size), log2(size) being the
     //      ideal depth (e.g. 4 leafs can be stored in a tree of depth 2)
@@ -333,7 +327,6 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::cleanAndBalance()
     )
     {
         chemisTree_.balance();
-        MRUList_.clear();
         treeModified = true;
     }
 
@@ -365,7 +358,7 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
         Rcq[i] = rhoi*Rphiq[s2c]/this->chemistry_.specieThermo()[s2c].W();
     }
     Rcq[speciesNumber] = Rphiq[Rphiq.size() - nAdditionalEqns_];
-    Rcq[speciesNumber+1] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 1];
+    Rcq[speciesNumber + 1] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 1];
     if (this->variableTimeStep())
     {
         Rcq[speciesNumber + 2] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 2];
@@ -373,15 +366,15 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
 
     // Aaa is computed implicitly,
     // A is given by A = C(psi0, t0+dt), where C is obtained through solving
-    // d/dt C(psi0,t) = J(psi(t))C(psi0,t)
+    // d/dt C(psi0, t) = J(psi(t))C(psi0, t)
     // If we solve it implicitly:
-    // (C(psi0, t0+dt) - C(psi0,t0))/dt = J(psi(t0+dt))C(psi0,t0+dt)
+    // (C(psi0, t0+dt) - C(psi0, t0))/dt = J(psi(t0+dt))C(psi0, t0+dt)
     // The Jacobian is thus computed according to the mapping
     // 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
-
-    this->chemistry_.jacobian(runTime_.value(), Rcq, A);
+    scalarField dcdt(speciesNumber + 2, Zero);
+    this->chemistry_.jacobian(runTime_.value(), Rcq, dcdt, A);
 
     // The jacobian is computed according to the molar concentration
     // the following conversion allows the code to use A with mass fraction
@@ -410,21 +403,48 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
         // Columns for pressure and temperature
         A(i, speciesNumber) *=
             -dt*this->chemistry_.specieThermo()[si].W()/rhoi;
-        A(i, speciesNumber+1) *=
+        A(i, speciesNumber + 1) *=
             -dt*this->chemistry_.specieThermo()[si].W()/rhoi;
     }
 
-    // For temperature and pressure, only unity on the diagonal
-    A(speciesNumber, speciesNumber) = 1;
-    A(speciesNumber + 1, speciesNumber + 1) = 1;
+    // For the temperature and pressure lines, ddc(dTdt)
+    // should be converted in ddY(dTdt)
+    for (label i=0; i<speciesNumber; i++)
+    {
+        label si = i;
+        if (mechRedActive)
+        {
+            si = this->chemistry_.simplifiedToCompleteIndex()[i];
+        }
+
+        A(speciesNumber, i) *=
+            -dt*rhoi/this->chemistry_.specieThermo()[si].W();
+        A(speciesNumber + 1, i) *=
+            -dt*rhoi/this->chemistry_.specieThermo()[si].W();
+    }
+
+    A(speciesNumber, speciesNumber) = -dt*A(speciesNumber, speciesNumber) + 1;
+
+    A(speciesNumber + 1, speciesNumber + 1) =
+        -dt*A(speciesNumber + 1, speciesNumber + 1) + 1;
+
     if (this->variableTimeStep())
     {
-        A[speciesNumber + 2][speciesNumber + 2] = 1;
+        A(speciesNumber + 2, speciesNumber + 2) = 1;
     }
 
     // Inverse of (I-dt*J(psi(t0+dt)))
     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++)
+    {
+        A(speciesNumber, i) = 0;
+        A(speciesNumber + 1, i) = 0;
+    }
 }
 
 
@@ -498,7 +518,7 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::retrieve
         }
         lastSearch_->lastTimeUsed() = this->chemistry_.timeSteps();
         addToMRU(phi0);
-        calcNewC(phi0,phiq, Rphiq);
+        calcNewC(phi0, phiq, Rphiq);
         nRetrieved_++;
         return true;
     }
@@ -525,11 +545,12 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
     // option is on, the code first tries to grow the point hold by lastSearch_
     if (lastSearch_ && growPoints_)
     {
-        if (grow(lastSearch_,phiq, Rphiq))
+        if (grow(lastSearch_, phiq, Rphiq))
         {
             nGrowth_++;
             growthOrAddFlag = 0;
-            // the structure of the tree is not modified, return false
+            addToMRU(lastSearch_);
+            //the structure of the tree is not modified, return false
             return growthOrAddFlag;
         }
     }
@@ -605,7 +626,10 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
         scaleFactor_.size(),
         lastSearch_ // lastSearch_ may be nullptr (handled by binaryTree)
     );
-
+    if (lastSearch_ != nullptr)
+    {
+        addToMRU(lastSearch_);
+    }
     nAdd_++;
 
     return growthOrAddFlag;
diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/radiationCoupledBase/radiationCoupledBase.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/radiationCoupledBase/radiationCoupledBase.C
index bc41bce27..3b410cd38 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/radiationCoupledBase/radiationCoupledBase.C
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/radiationCoupledBase/radiationCoupledBase.C
@@ -28,7 +28,7 @@ License
 #include "mappedPatchBase.H"
 #include "fvPatchFieldMapper.H"
 #include "radiationModel.H"
-#include "noRadiation.H"
+#include "opaqueSolid.H"
 #include "absorptionEmissionModel.H"
 
 // * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * * //
@@ -151,29 +151,20 @@ Foam::scalarField Foam::radiationCoupledBase::emissivity() const
 
             const polyMesh& nbrMesh = mpp.sampleMesh();
 
-            const radiation::radiationModel& radiation =
-                nbrMesh.lookupObject<radiation::radiationModel>
+            const radiation::opaqueSolid& radiation =
+                nbrMesh.lookupObject<radiation::opaqueSolid>
                 (
                     "radiationProperties"
                 );
 
-            if (isType<radiation::noRadiation>(radiation))
-            {
-                FatalErrorInFunction
-                    << "No radiation model defined for region "
-                    << nbrMesh.name() << nl
-                    << "    required for option 'solidRadiation' "
-                       "of patchField type " << type()
-                    << " on patch " << patch_.name()
-                    << abort(FatalError);
-            }
-
             const fvMesh& nbrFvMesh = refCast<const fvMesh>(nbrMesh);
 
             const fvPatch& nbrPatch =
                 nbrFvMesh.boundary()[mpp.samplePolyPatch().index()];
 
-
+            // NOTE: for an opaqueSolid the absorptionEmission model returns the
+            // emissivity of the surface rather than the emission coefficient
+            // and the input specification MUST correspond to this.
             scalarField emissivity
             (
                 radiation.absorptionEmission().e()().boundaryField()
diff --git a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L
index 552d759d7..73cc72a70 100644
--- a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L
+++ b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L
@@ -651,6 +651,18 @@ bool finishReaction = false;
     }
 
 <readThermoSpecieName>{end} {
+        Reaction<gasHThermoPhysics>::TlowDefault = max
+        (
+            Reaction<gasHThermoPhysics>::TlowDefault,
+            currentLowT
+        );
+
+        Reaction<gasHThermoPhysics>::ThighDefault = min
+        (
+            Reaction<gasHThermoPhysics>::ThighDefault,
+            currentHighT
+        );
+
         BEGIN(INITIAL);
     }
 
diff --git a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C
index 029956ba8..710192740 100644
--- a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C
+++ b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C
@@ -26,6 +26,7 @@ License
 #include "chemkinReader.H"
 #include <fstream>
 #include "atomicWeights.H"
+#include "ReactionProxy.H"
 #include "IrreversibleReaction.H"
 #include "ReversibleReaction.H"
 #include "NonEquilibriumReversibleReaction.H"
@@ -182,7 +183,7 @@ void Foam::chemkinReader::addReactionType
                 new IrreversibleReaction
                 <Reaction, gasHThermoPhysics, ReactionRateType>
                 (
-                    Reaction<gasHThermoPhysics>
+                    ReactionProxy<gasHThermoPhysics>
                     (
                         speciesTable_,
                         lhs.shrink(),
@@ -202,7 +203,7 @@ void Foam::chemkinReader::addReactionType
                 new ReversibleReaction
                 <Reaction, gasHThermoPhysics, ReactionRateType>
                 (
-                    Reaction<gasHThermoPhysics>
+                    ReactionProxy<gasHThermoPhysics>
                     (
                         speciesTable_,
                         lhs.shrink(),
@@ -490,7 +491,7 @@ void Foam::chemkinReader::addReaction
                     new NonEquilibriumReversibleReaction
                         <Reaction, gasHThermoPhysics, ArrheniusReactionRate>
                     (
-                        Reaction<gasHThermoPhysics>
+                        ReactionProxy<gasHThermoPhysics>
                         (
                             speciesTable_,
                             lhs.shrink(),
@@ -546,7 +547,7 @@ void Foam::chemkinReader::addReaction
                         thirdBodyArrheniusReactionRate
                     >
                     (
-                        Reaction<gasHThermoPhysics>
+                        ReactionProxy<gasHThermoPhysics>
                         (
                             speciesTable_,
                             lhs.shrink(),
@@ -646,9 +647,13 @@ void Foam::chemkinReader::addReaction
                 reactions_.append
                 (
                     new NonEquilibriumReversibleReaction
-                        <Reaction, gasHThermoPhysics, LandauTellerReactionRate>
+                    <
+                        Reaction,
+                        gasHThermoPhysics,
+                        LandauTellerReactionRate
+                    >
                     (
-                        Reaction<gasHThermoPhysics>
+                        ReactionProxy<gasHThermoPhysics>
                         (
                             speciesTable_,
                             lhs.shrink(),
@@ -780,6 +785,9 @@ void Foam::chemkinReader::read
     const fileName& transportFileName
 )
 {
+    Reaction<gasHThermoPhysics>::TlowDefault = 0;
+    Reaction<gasHThermoPhysics>::ThighDefault = great;
+
     transportDict_.read(IFstream(transportFileName)());
 
     if (thermoFileName != fileName::null)
diff --git a/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.H b/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.H
index 84e4f61f2..a98126696 100644
--- a/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.H
+++ b/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.H
@@ -59,7 +59,7 @@ inline Ostream& operator<<(Ostream&, const solidReaction<ReactionThermo>&);
 template<class ReactionThermo>
 class solidReaction
 :
-    public  Reaction<ReactionThermo>
+    public Reaction<ReactionThermo>
 {
 
 private:
@@ -129,28 +129,6 @@ public:
         );
 
 
-        //- Construct and return a clone
-        virtual autoPtr<Reaction<ReactionThermo>> clone() const
-        {
-            return autoPtr<Reaction<ReactionThermo>>
-            (
-                new solidReaction<ReactionThermo>(*this)
-            );
-        }
-
-        //- Construct and return a clone with new speciesTable
-        virtual autoPtr<Reaction<ReactionThermo>> clone
-        (
-            const speciesTable& species
-        ) const
-        {
-            return autoPtr<Reaction<ReactionThermo>>
-            (
-                new solidReaction<ReactionThermo>(*this, species)
-            );
-        }
-
-
     //- Destructor
     virtual ~solidReaction()
     {}
@@ -160,12 +138,13 @@ public:
 
         // Access
 
-
-            // - Access to gas components of the reaction
+            //- Access to the gas components of the left hand side
             virtual const List<specieCoeffs>& grhs() const;
+
+            //- Access to the gas components of the right hand side
             virtual const List<specieCoeffs>& glhs() const;
 
-            // - Access to gas specie list
+            //- Access to gas specie list
             virtual const speciesTable& gasSpecies() const;
 
 
diff --git a/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRate.H b/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRate.H
index c5642e627..86ff4639d 100644
--- a/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRate.H
+++ b/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,6 +104,37 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) 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;
+
+        //- Species concentration derivative of the pressure dependent term
+        //  By default this value is 1 as it multiplies the third-body term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            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;
+
 
          //- Write to stream
         inline void write(Ostream& os) const;
diff --git a/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRateI.H b/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRateI.H
index cd64c2f36..bf32298db 100644
--- a/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRateI.H
+++ b/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRateI.H
@@ -60,18 +60,60 @@ inline Foam::scalar Foam::solidArrheniusReactionRate::operator()
     const scalarField&
 ) const
 {
-    scalar ak = A_;
-
     if (T < Tcrit_)
     {
-        ak *= 0.0;
+        return 0;
     }
     else
     {
-        ak *= exp(-Ta_/T);
+        return A_*exp(-Ta_/T);
     }
+}
 
-    return ak;
+
+inline Foam::scalar Foam::solidArrheniusReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField&
+) const
+{
+    if (T < Tcrit_)
+    {
+        return 0;
+    }
+    else
+    {
+        return A_*exp(-Ta_/T)*Ta_/sqr(T);
+    }
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::solidArrheniusReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::solidArrheniusReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::solidArrheniusReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C b/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C
index 32a301fab..25b843724 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -106,6 +106,160 @@ Foam::scalar Foam::IrreversibleReaction
 }
 
 
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::kr
+(
+    const scalar kfwd,
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::kr
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkfdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return k_.ddT(p, T, c);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkrdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const scalar dkfdT,
+    const scalar kr
+) const
+{
+    return 0;
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::beta() const
+{
+    return k_.beta();
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+void Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    k_.dcidc(p, T, c, dcidc);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return k_.dcidT(p, T, c);
+}
+
+
 template
 <
     template<class> class ReactionType,
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H b/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H
index 3c3aa7ead..12771b748 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,8 +25,7 @@ Class
     Foam::IrreversibleReaction
 
 Description
-    Simple extension of Reaction to handle reversible reactions using
-    equilibrium thermodynamics.
+    Simple extension of Reaction to handle irreversible reactions
 
 SourceFiles
     IrreversibleReaction.C
@@ -164,6 +163,71 @@ public:
                 const scalarField& c
             ) const;
 
+            //- Reverse rate constant from the given forward rate constant
+            //  Returns 0
+            virtual scalar kr
+            (
+                const scalar kfwd,
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Reverse rate constant
+            //  Returns 0
+            virtual scalar kr
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+
+        // IrreversibleReaction Jacobian functions
+
+            //- Temperature derivative of forward rate
+            virtual scalar dkfdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Temperature derivative of reverse rate
+            //  Returns 0
+            virtual scalar dkrdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const scalar dkfdT,
+                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;
+
+            //- Species concentration derivative of the pressure dependent term
+            //  By default this value is 1 as it multiplies the third-body term
+            virtual void dcidc
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcidc
+            ) const;
+
+            //- Temperature derivative of the pressure dependent term
+            //  By default this value is 0 since ddT of molecularity is approx.0
+            virtual scalar dcidT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
 
         //- Write
         virtual void write(Ostream&) const;
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C
index 255cf2b93..d3b43b8c6 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -178,6 +178,117 @@ Foam::NonEquilibriumReversibleReaction
 }
 
 
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar
+Foam::NonEquilibriumReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkfdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return fk_.ddT(p, T, c);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar
+Foam::NonEquilibriumReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkrdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const scalar dkfdT,
+    const scalar kr
+) const
+{
+    return rk_.ddT(p, T, c);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::NonEquilibriumReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::beta() const
+{
+    return rk_.beta();
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+void Foam::NonEquilibriumReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    fk_.dcidc(p, T, c, dcidc);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::NonEquilibriumReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return fk_.dcidT(p, T, c);
+}
+
+
 template
 <
     template<class> class ReactionType,
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H
index 8e76691d3..939c0e680 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -167,6 +167,51 @@ public:
             ) const;
 
 
+        // ReversibleReaction Jacobian functions
+
+            //- Temperature derivative of forward rate
+            virtual scalar dkfdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Temperature derivative of backward rate
+            virtual scalar dkrdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const scalar dkfdT,
+                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;
+
+            //- Species concentration derivative of the pressure dependent term
+            //  By default this value is 1 as it multiplies the third-body term
+            virtual void dcidc
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcidc
+            ) const;
+
+            //- Temperature derivative of the pressure dependent term
+            //  By default this value is 0 since ddT of molecularity is approx.0
+            virtual scalar dcidT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+
         //- Write
         virtual void write(Ostream&) const;
 };
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
index dc2392eca..f17036f37 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
@@ -29,7 +29,13 @@ License
 // * * * * * * * * * * * * * * * * Static Data * * * * * * * * * * * * * * * //
 
 template<class ReactionThermo>
-Foam::label Foam::Reaction<ReactionThermo>::nUnNamedReactions = 0;
+Foam::label Foam::Reaction<ReactionThermo>::nUnNamedReactions(0);
+
+template<class ReactionThermo>
+Foam::scalar Foam::Reaction<ReactionThermo>::TlowDefault(0);
+
+template<class ReactionThermo>
+Foam::scalar Foam::Reaction<ReactionThermo>::ThighDefault(great);
 
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
@@ -160,6 +166,8 @@ Foam::Reaction<ReactionThermo>::Reaction
     ReactionThermo::thermoType(*thermoDatabase[species[0]]),
     name_("un-named-reaction-" + Foam::name(getNewReactionID())),
     species_(species),
+    Tlow_(TlowDefault),
+    Thigh_(ThighDefault),
     lhs_(lhs),
     rhs_(rhs)
 {
@@ -177,6 +185,8 @@ Foam::Reaction<ReactionThermo>::Reaction
     ReactionThermo::thermoType(r),
     name_(r.name() + "Copy"),
     species_(species),
+    Tlow_(r.Tlow()),
+    Thigh_(r.Thigh()),
     lhs_(r.lhs_),
     rhs_(r.rhs_)
 {}
@@ -197,7 +207,7 @@ Foam::Reaction<ReactionThermo>::specieCoeffs::specieCoeffs
     }
     else
     {
-        stoichCoeff = 1.0;
+        stoichCoeff = 1;
     }
 
     exponent = stoichCoeff;
@@ -330,7 +340,9 @@ Foam::Reaction<ReactionThermo>::Reaction
 :
     ReactionThermo::thermoType(*thermoDatabase[species[0]]),
     name_(dict.dictName()),
-    species_(species)
+    species_(species),
+    Tlow_(dict.lookupOrDefault<scalar>("Tlow", TlowDefault)),
+    Thigh_(dict.lookupOrDefault<scalar>("Thigh", ThighDefault))
 {
     setLRhs
     (
@@ -388,39 +400,402 @@ void Foam::Reaction<ReactionThermo>::write(Ostream& os) const
 
 
 template<class ReactionThermo>
-Foam::scalar Foam::Reaction<ReactionThermo>::kf
+void Foam::Reaction<ReactionThermo>::ddot
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    scalarField& d
 ) const
 {
-    return 0.0;
 }
 
 
 template<class ReactionThermo>
-Foam::scalar Foam::Reaction<ReactionThermo>::kr
+void Foam::Reaction<ReactionThermo>::fdot
 (
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& f
+) const
+{
+}
+
+
+template<class ReactionThermo>
+void Foam::Reaction<ReactionThermo>::omega
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcdt
+) const
+{
+    scalar pf, cf, pr, cr;
+    label lRef, rRef;
+
+    scalar omegaI = omega
+    (
+        p, T, c, pf, cf, lRef, pr, cr, rRef
+    );
+
+    forAll(lhs_, i)
+    {
+        const label si = lhs_[i].index;
+        const scalar sl = lhs_[i].stoichCoeff;
+        dcdt[si] -= sl*omegaI;
+    }
+    forAll(rhs_, i)
+    {
+        const label si = rhs_[i].index;
+        const scalar sr = rhs_[i].stoichCoeff;
+        dcdt[si] += sr*omegaI;
+    }
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::Reaction<ReactionThermo>::omega
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalar& pf,
+    scalar& cf,
+    label& lRef,
+    scalar& pr,
+    scalar& cr,
+    label& rRef
+) const
+{
+
+    scalar clippedT = min(max(T, this->Tlow()), this->Thigh());
+
+    const scalar kf = this->kf(p, clippedT, c);
+    const scalar kr = this->kr(kf, p, clippedT, c);
+
+    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++)
+    {
+        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);
+        }
+    }
+    cf = max(c[lRef], 0);
+
+    {
+        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);
+        }
+    }
+
+    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;
+}
+
+
+template<class ReactionThermo>
+void Foam::Reaction<ReactionThermo>::dwdc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarSquareMatrix& J,
+    scalarField& dcdt,
+    scalar& omegaI,
+    scalar& kfwd,
+    scalar& kbwd,
+    const bool reduced,
+    const List<label>& c2s
+) const
+{
+
+    scalar pf, cf, pr, cr;
+    label lRef, rRef;
+
+    omegaI = omega(p, T, c, pf, cf, lRef, pr, cr, rRef);
+
+    forAll(lhs_, i)
+    {
+        const label si = reduced ? c2s[lhs_[i].index] : lhs_[i].index;
+        const scalar sl = lhs_[i].stoichCoeff;
+        dcdt[si] -= sl*omegaI;
+    }
+    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);
+    kbwd = this->kr(kfwd, p, T, c);
+
+    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)
+            {
+                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);
+                }
+            }
+            else
+            {
+                kf *= pow(c[si], el);
+            }
+        }
+
+        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)
+            {
+                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);
+                }
+            }
+            else
+            {
+                kr *= pow(c[si], er);
+            }
+        }
+
+        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
+    if (this->beta().size() > 0)
+    {
+        scalarField dcidc(this->beta().size());
+        this->dcidc(p, T, c, dcidc);
+        forAll(this->beta(), j)
+        {
+            label sj = this-> 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;
+                }
+            }
+        }
+    }
+}
+
+
+template<class ReactionThermo>
+void Foam::Reaction<ReactionThermo>::dwdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const scalar omegaI,
     const scalar kfwd,
-    const scalar p,
-    const scalar T,
-    const scalarField& c
+    const scalar kbwd,
+    scalarSquareMatrix& J,
+    const bool reduced,
+    const List<label>& c2s,
+    const label indexT
 ) const
 {
-    return 0.0;
-}
+    scalar kf = kfwd;
+    scalar kr = kbwd;
 
+    scalar dkfdT = this->dkfdT(p, T, c);
+    scalar dkrdT = this->dkrdT(p, T, c, dkfdT, kr);
 
-template<class ReactionThermo>
-Foam::scalar Foam::Reaction<ReactionThermo>::kr
-(
-    const scalar p,
-    const scalar T,
-    const scalarField& c
-) const
-{
-    return 0.0;
+    scalar sumExp = 0.0;
+    forAll(lhs_, i)
+    {
+        const label si = lhs_[i].index;
+        const scalar el = lhs_[i].exponent;
+        const scalar cExp = pow(c[si], el);
+        dkfdT *= cExp;
+        kf *= cExp;
+        sumExp += el;
+    }
+    kf *= -sumExp/T;
+
+    sumExp = 0.0;
+    forAll(rhs_, i)
+    {
+        const label si = rhs_[i].index;
+        const scalar er = rhs_[i].exponent;
+        const scalar cExp = pow(c[si], er);
+        dkrdT *= cExp;
+        kr *= cExp;
+        sumExp += er;
+    }
+    kr *= -sumExp/T;
+
+    // dqidT includes the third-body (or pressure dependent) effect
+    scalar dqidT = dkfdT - dkrdT + kf - kr;
+
+    // For reactions including third-body efficiencies or pressure dependent
+    // reaction, an additional term is needed
+    scalar dcidT = this->dcidT(p, T, c);
+    dcidT *= omegaI;
+
+    // J(i, indexT) = sum_reactions nu_i dqdT
+    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);
+    }
+    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);
+    }
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
index 13f6d9a53..dce5e7688 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
@@ -40,6 +40,8 @@ SourceFiles
 #include "speciesTable.H"
 #include "HashPtrTable.H"
 #include "scalarField.H"
+#include "simpleMatrix.H"
+#include "Tuple2.H"
 #include "typeInfo.H"
 #include "runTimeSelectionTables.H"
 
@@ -84,6 +86,8 @@ public:
         //- Number of un-named reactions
         static label nUnNamedReactions;
 
+        //- Default temperature limits of applicability of reaction rates
+        static scalar TlowDefault, ThighDefault;
 
     // Public data types
 
@@ -134,6 +138,9 @@ private:
         //- List of specie names present in reaction system
         const speciesTable& species_;
 
+        //- Temperature limits of applicability of reaction rates
+        scalar Tlow_, Thigh_;
+
         //- Specie info for the left-hand-side of the reaction
         List<specieCoeffs> lhs_;
 
@@ -201,25 +208,13 @@ public:
         );
 
         //- Construct and return a clone
-        virtual autoPtr<Reaction<ReactionThermo>> clone() const //
-        {
-            return autoPtr<Reaction<ReactionThermo>>
-            (
-                new Reaction<ReactionThermo>(*this)
-            );
-        }
+        virtual autoPtr<Reaction<ReactionThermo>> clone() const = 0;
 
         //- Construct and return a clone with new speciesTable
         virtual autoPtr<Reaction<ReactionThermo>> clone
         (
             const speciesTable& species
-        ) const
-        {
-            return autoPtr<Reaction<ReactionThermo>>
-            (
-                new Reaction<ReactionThermo>(*this, species)
-            );
-        }
+        ) const = 0;
 
 
     // Selectors
@@ -242,14 +237,22 @@ public:
 
         // Access
 
+            inline word& name();
             inline const word& name() const;
 
-            // - Access to basic components of the reaction
+            inline scalar Tlow() const;
+            inline scalar Thigh() const;
+
+            //- Access to the components of the left hand side
             inline const List<specieCoeffs>& lhs() const;
+
+            //- Access to the components of the right hand side
             inline const List<specieCoeffs>& rhs() const;
 
-            // - Access to gas components of the reaction
+            //- Access to the gas components of the left hand side
             virtual const List<specieCoeffs>& grhs() const;
+
+            //- Access to the gas components of the right hand side
             virtual const List<specieCoeffs>& glhs() const;
 
             //- Access to specie list
@@ -268,6 +271,49 @@ public:
             );
 
 
+        // Reaction rate coefficients
+
+            //- Forward reaction rate
+            void ddot
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& d
+            ) const;
+
+            //- Backward reaction rate
+            void fdot
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& f
+            ) const;
+
+            //- Net reaction rate for individual species
+            void omega
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcdt
+            ) const;
+
+            //- Net reaction rate
+            scalar omega
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalar& pf,
+                scalar& cf,
+                label& lRef,
+                scalar& pr,
+                scalar& cr,
+                label& rRef
+            ) const;
+
         // Reaction rate coefficients
 
             //- Forward rate constant
@@ -276,7 +322,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c
-            ) const;
+            ) const = 0;
 
             //- Reverse rate constant from the given forward rate constant
             virtual scalar kr
@@ -285,18 +331,91 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c
-            ) const;
+            ) const = 0;
 
-            //- Reverse rate constant.
-            //  Note this evaluates the forward rate constant and divides by the
-            //  equilibrium constant
+            //- Reverse rate constant
             virtual scalar kr
             (
                 const scalar p,
                 const scalar T,
                 const scalarField& c
+            ) const = 0;
+
+
+        // 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,
+                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 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
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const = 0;
+
+            //- Temperature derivative of reverse rate
+            virtual scalar dkrdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const scalar dkfdT,
+                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;
+
+            //- Species concentration derivative of the pressure dependent term
+            virtual void dcidc
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcidc
+            ) const = 0;
+
+            //- Temperature derivative of the pressure dependent term
+            virtual scalar dcidT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const = 0;
+
 
         //- Write
         virtual void write(Ostream&) const;
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/ReactionI.H b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/ReactionI.H
index eba046f64..b52565c63 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/ReactionI.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/ReactionI.H
@@ -32,6 +32,13 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class ReactionThermo>
+inline word& Reaction<ReactionThermo>::name()
+{
+    return name_;
+}
+
+
 template<class ReactionThermo>
 inline const word& Reaction<ReactionThermo>::name() const
 {
@@ -39,6 +46,20 @@ inline const word& Reaction<ReactionThermo>::name() const
 }
 
 
+template<class ReactionThermo>
+inline scalar Reaction<ReactionThermo>::Tlow() const
+{
+    return Tlow_;
+}
+
+
+template<class ReactionThermo>
+inline scalar Reaction<ReactionThermo>::Thigh() const
+{
+    return Thigh_;
+}
+
+
 template<class ReactionThermo>
 inline const List<typename Reaction<ReactionThermo>::specieCoeffs>&
 Reaction<ReactionThermo>::lhs() const
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReactionList/ReactionList.C b/src/thermophysicalModels/specie/reaction/Reactions/ReactionList/ReactionList.C
index 6d3f2539f..cfae96819 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/ReactionList/ReactionList.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/ReactionList/ReactionList.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -84,6 +84,13 @@ bool Foam::ReactionList<ThermoType>::readReactionDict()
 {
     const dictionary& reactions(dict_.subDict("reactions"));
 
+    // Set general temperature limits from the dictionary
+    Reaction<ThermoType>::TlowDefault =
+        dict_.lookupOrDefault<scalar>("Tlow", 0);
+
+    Reaction<ThermoType>::ThighDefault =
+        dict_.lookupOrDefault<scalar>("Thigh", great);
+
     forAllConstIter(dictionary, reactions, iter)
     {
         const word reactionName = iter().keyword();
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C b/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C
new file mode 100644
index 000000000..992521fd5
--- /dev/null
+++ b/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C
@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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 "ReactionProxy.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ReactionThermo>
+Foam::ReactionProxy<ReactionThermo>::ReactionProxy
+(
+    const speciesTable& species,
+    const List<typename Reaction<ReactionThermo>::specieCoeffs>& lhs,
+    const List<typename Reaction<ReactionThermo>::specieCoeffs>& rhs,
+    const HashPtrTable<ReactionThermo>& thermoDatabase
+)
+:
+    Reaction<ReactionThermo>
+    (
+        species,
+        lhs,
+        rhs,
+        thermoDatabase
+    )
+{}
+
+
+template<class ReactionThermo>
+Foam::ReactionProxy<ReactionThermo>::ReactionProxy
+(
+    const Reaction<ReactionThermo>& r,
+    const speciesTable& species
+)
+:
+    Reaction<ReactionThermo>
+    (
+        r,
+        species
+    )
+{}
+
+
+template<class ReactionThermo>
+Foam::ReactionProxy<ReactionThermo>::ReactionProxy
+(
+    const speciesTable& species,
+    const HashPtrTable<ReactionThermo>& thermoDatabase,
+    const dictionary& dict
+)
+:
+    Reaction<ReactionThermo>
+    (
+        species,
+        thermoDatabase,
+        dict
+    )
+{}
+
+
+template<class ReactionThermo>
+Foam::autoPtr<Foam::Reaction<ReactionThermo>>
+Foam::ReactionProxy<ReactionThermo>::clone() const
+{
+    NotImplemented;
+    return autoPtr<Foam::Reaction<ReactionThermo>>();
+}
+
+
+template<class ReactionThermo>
+Foam::autoPtr<Foam::Reaction<ReactionThermo>>
+Foam::ReactionProxy<ReactionThermo>::clone
+(
+    const speciesTable& species
+) const
+{
+    NotImplemented;
+    return autoPtr<Reaction<ReactionThermo>>();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::kf
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::kr
+(
+    const scalar kfwd,
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::kr
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::dkfdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::dkrdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const scalar dkfdT,
+    const scalar kr
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+template<class ReactionThermo>
+const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::ReactionProxy<ReactionThermo>::beta() const
+{
+    NotImplemented;
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+template<class ReactionThermo>
+void Foam::ReactionProxy<ReactionThermo>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    NotImplemented;
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H b/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H
new file mode 100644
index 000000000..2349dd27c
--- /dev/null
+++ b/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H
@@ -0,0 +1,189 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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::ReactionProxy
+
+Description
+    Proxy version of Reaction which provides dummy implementations of the
+    abstract virtual functions.
+
+    Used for read-construction and format conversion.
+
+SourceFiles
+    ReactionProxy.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ReactionProxy_H
+#define ReactionProxy_H
+
+#include "Reaction.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class ReactionProxy Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ReactionThermo>
+class ReactionProxy
+:
+    public Reaction<ReactionThermo>
+{
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        ReactionProxy
+        (
+            const speciesTable& species,
+            const List<typename Reaction<ReactionThermo>::specieCoeffs>& lhs,
+            const List<typename Reaction<ReactionThermo>::specieCoeffs>& rhs,
+            const HashPtrTable<ReactionThermo>& thermoDatabase
+        );
+
+        //- Construct as copy given new speciesTable
+        ReactionProxy
+        (
+            const Reaction<ReactionThermo>&,
+            const speciesTable& species
+        );
+
+        //- Construct from dictionary
+        ReactionProxy
+        (
+            const speciesTable& species,
+            const HashPtrTable<ReactionThermo>& thermoDatabase,
+            const dictionary& dict
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<Reaction<ReactionThermo>> clone() const;
+
+        //- Construct and return a clone with new speciesTable
+        virtual autoPtr<Reaction<ReactionThermo>> clone
+        (
+            const speciesTable& species
+        ) const;
+
+
+    //- Destructor
+    virtual ~ReactionProxy()
+    {}
+
+
+    // Member Functions
+
+        // Reaction rate coefficients
+
+            //- Forward rate constant
+            virtual scalar kf
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Reverse rate constant from the given forward rate constant
+            virtual scalar kr
+            (
+                const scalar kfwd,
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Reverse rate constant
+            virtual scalar kr
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+
+        // Jacobian coefficients
+
+            //- Temperature derivative of forward rate
+            virtual scalar dkfdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Temperature derivative of reverse rate
+            virtual scalar dkrdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const scalar dkfdT,
+                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;
+
+            //- Species concentration derivative of the pressure dependent term
+            virtual void dcidc
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcidc
+            ) const;
+
+            //- Temperature derivative of the pressure dependent term
+            virtual scalar dcidT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "ReactionProxy.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C
index 2ccc2fd88..58964a9ba 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C
@@ -151,6 +151,117 @@ Foam::scalar Foam::ReversibleReaction
 }
 
 
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::ReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkfdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return k_.ddT(p, T, c);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::ReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkrdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const scalar dkfdT,
+    const scalar kr
+) const
+{
+    scalar Kc = max(this->Kc(p, T), rootSmall);
+
+    return dkfdT/Kc - kr*this->dKcdTbyKc(p, T);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::ReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::beta() const
+{
+    return k_.beta();
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+void Foam::ReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    k_.dcidc(p, T, c, dcidc);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::ReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return k_.dcidT(p, T, c);
+}
+
+
 template
 <
     template<class> class ReactionType,
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H
index 798b75759..9cb1c4d64 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -180,6 +180,51 @@ public:
             ) const;
 
 
+        // ReversibleReaction Jacobian functions
+
+            //- Temperature derivative of forward rate
+            virtual scalar dkfdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Temperature derivative of backward rate
+            virtual scalar dkrdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const scalar dkfdT,
+                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;
+
+            //- Species concentration derivative of the pressure dependent term
+            //  By default this value is 1 as it multiplies the third-body term
+            virtual void dcidc
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcidc
+            ) const;
+
+            //- Temperature derivative of the pressure dependent term
+            //  By default this value is 0 since ddT of molecularity is approx.0
+            virtual scalar dcidT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+
         //- Write
         virtual void write(Ostream&) const;
 };
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H
index fbe89154b..604161773 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -100,6 +100,37 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) 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;
+
+        //- Species concentration derivative of the pressure dependent term
+        //  By default this value is 1 as it multiplies the third-body term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            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;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H
index a856eb911..61bb2e1b6 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H
@@ -75,6 +75,57 @@ inline Foam::scalar Foam::ArrheniusReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::ArrheniusReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField&
+) const
+{
+    scalar ak = A_;
+
+    if (mag(beta_) > vSmall)
+    {
+        ak *= pow(T, beta_);
+    }
+
+    if (mag(Ta_) > vSmall)
+    {
+        ak *= exp(-Ta_/T);
+    }
+
+    return ak*(beta_+Ta_/T)/T;
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::ArrheniusReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::ArrheniusReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::ArrheniusReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::ArrheniusReactionRate::write(Ostream& os) const
 {
     os.writeKeyword("A") << A_ << token::END_STATEMENT << nl;
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H
index 87eab0d40..e70a3977b 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,6 +69,7 @@ class ChemicallyActivatedReactionRate
         ReactionRate kInf_;
         ChemicallyActivationFunction F_;
         thirdBodyEfficiencies thirdBodyEfficiencies_;
+        List<Tuple2<label, scalar>> beta_;
 
 
 public:
@@ -109,6 +110,40 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) 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_;
+        }
+
+        //- Species concentration derivative of the pressure dependent term
+        //  By default this value is 1 as it multiplies the third-body term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            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;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H
index 3e54f9a6a..d85a8c5bb 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,7 +42,12 @@ 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>
@@ -60,7 +65,19 @@ inline Foam::ChemicallyActivatedReactionRate
     kInf_(species, dict),
     F_(dict),
     thirdBodyEfficiencies_(species, dict)
-{}
+{
+    forAll(thirdBodyEfficiencies_, i)
+    {
+        beta_.append
+        (
+            Tuple2<label, scalar>
+            (
+                i,
+                thirdBodyEfficiencies_[i]
+            )
+        );
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -77,15 +94,108 @@ inline Foam::scalar Foam::ChemicallyActivatedReactionRate
     const scalarField& c
 ) const
 {
-    scalar k0 = k0_(p, T, c);
-    scalar kInf = kInf_(p, T, c);
-
-    scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+    const scalar k0 = k0_(p, T, c);
+    const scalar kInf = kInf_(p, T, c);
+    const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
 
     return k0*(1/(1 + Pr))*F_(T, Pr);
 }
 
 
+template<class ReactionRate, class ChemicallyActivationFunction>
+inline Foam::scalar Foam::ChemicallyActivatedReactionRate
+<
+    ReactionRate,
+    ChemicallyActivationFunction
+>::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    const scalar k0 = k0_(p, T, c);
+    const scalar kInf = kInf_(p, T, c);
+    const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+
+    return (1/(1 + Pr))*F_(T, Pr)*k0_.ddT(p, T, c);
+}
+
+
+template<class ReactionRate, class ChemicallyActivationFunction>
+inline void Foam::ChemicallyActivatedReactionRate
+<
+    ReactionRate,
+    ChemicallyActivationFunction
+>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    const scalar M = thirdBodyEfficiencies_.M(c);
+
+    if (M > small)
+    {
+        const scalar k0 = k0_(p, T, c);
+        const scalar kInf = kInf_(p, T, c);
+
+        const scalar Pr = k0*M/kInf;
+        const scalar F = F_(T, Pr);
+
+        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;
+        }
+    }
+}
+
+
+template<class ReactionRate, class ChemicallyActivationFunction>
+inline Foam::scalar Foam::ChemicallyActivatedReactionRate
+<
+    ReactionRate,
+    ChemicallyActivationFunction
+>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    const scalar M = thirdBodyEfficiencies_.M(c);
+
+    if (M > small)
+    {
+        const scalar k0 = k0_(p, T, c);
+        const scalar kInf = kInf_(p, T, c);
+
+        const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+        const scalar F = F_(T, Pr);
+        const scalar dPrdT =
+            Pr*(k0_.ddT(p, T, c)/k0 - kInf_.ddT(p, T, c)/kInf - 1/T);
+        const scalar dFdT = F_.ddT(Pr, F, dPrdT, T);
+
+        return (-dPrdT/(1 + Pr) + dFdT/F);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
 template<class ReactionRate, class ChemicallyActivationFunction>
 inline void Foam::ChemicallyActivatedReactionRate
 <
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H
index e89aa8949..eb2498f60 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,6 +68,7 @@ class FallOffReactionRate
         ReactionRate kInf_;
         FallOffFunction F_;
         thirdBodyEfficiencies thirdBodyEfficiencies_;
+        List<Tuple2<label, scalar>> beta_;
 
 
 public:
@@ -106,6 +107,38 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) 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_;
+        }
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H
index fff348d48..228776e61 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,7 +39,12 @@ FallOffReactionRate
     kInf_(kInf),
     F_(F),
     thirdBodyEfficiencies_(tbes)
-{}
+{
+    forAll(tbes, i)
+    {
+        beta_.append(Tuple2<label, scalar>(i, tbes[i]));
+    }
+}
 
 
 template<class ReactionRate, class FallOffFunction>
@@ -54,7 +59,19 @@ 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  * * * * * * * * * * * * * //
@@ -68,15 +85,97 @@ Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::operator()
     const scalarField& c
 ) const
 {
-    scalar k0 = k0_(p, T, c);
-    scalar kInf = kInf_(p, T, c);
-
-    scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+    const scalar k0 = k0_(p, T, c);
+    const scalar kInf = kInf_(p, T, c);
+    const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
 
     return kInf*(Pr/(1 + Pr))*F_(T, Pr);
 }
 
 
+template<class ReactionRate, class FallOffFunction>
+inline Foam::scalar
+Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    const scalar k0 = k0_(p, T, c);
+    const scalar kInf = kInf_(p, T, c);
+    const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+
+    return (Pr/(1 + Pr))*F_(T, Pr)*kInf_.ddT(p, T, c);
+}
+
+
+template<class ReactionRate, class FallOffFunction>
+inline void Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    const scalar M = thirdBodyEfficiencies_.M(c);
+
+    if (M > small)
+    {
+        const scalar k0 = k0_(p, T, c);
+        const scalar kInf = kInf_(p, T, c);
+        const scalar Pr = k0*M/kInf;
+        const scalar F = F_(T, Pr);
+
+        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;
+        }
+    }
+}
+
+
+template<class ReactionRate, class FallOffFunction>
+inline Foam::scalar
+Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    const scalar M = thirdBodyEfficiencies_.M(c);
+
+    if (M > small)
+    {
+        const scalar k0 = k0_(p, T, c);
+        const scalar kInf = kInf_(p, T, c);
+
+        const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+        const scalar F = F_(T, Pr);
+        const scalar dPrdT =
+            Pr*(k0_.ddT(p, T, c)/k0 - kInf_.ddT(p, T, c)/kInf - 1/T);
+        const scalar dFdT = F_.ddT(Pr, F, dPrdT, T);
+
+        return (dPrdT/(Pr*(1 + Pr)) + dFdT/F);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
 template<class ReactionRate, class FallOffFunction>
 inline void Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::write
 (
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H
index ab37e054f..a1b813ab4 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -103,6 +103,35 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) 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;
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H
index 636c3b0c2..cdfb3dc73 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H
@@ -69,7 +69,7 @@ inline Foam::scalar Foam::JanevReactionRate::operator()
         lta *= pow(T, beta_);
     }
 
-    scalar expArg = 0.0;
+    scalar expArg = 0;
 
     if (mag(Ta_) > vSmall)
     {
@@ -89,6 +89,75 @@ inline Foam::scalar Foam::JanevReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::JanevReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField&
+) const
+{
+    scalar lta = A_;
+
+    if (mag(beta_) > vSmall)
+    {
+        lta *= pow(T, beta_);
+    }
+
+    scalar expArg = 0;
+
+    if (mag(Ta_) > vSmall)
+    {
+        expArg -= Ta_/T;
+    }
+
+    scalar lnT = log(T);
+
+    for (int n=0; n<nb_; n++)
+    {
+        expArg += b_[n]*pow(lnT, n);
+    }
+
+    scalar deriv = b_[1];
+
+    for (int n=2; n<nb_; n++)
+    {
+        deriv += n*b_[n]*pow(lnT, n-1);
+    }
+
+    lta *= exp(expArg);
+
+    return lta*(beta_+Ta_/T+deriv)/T;
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::JanevReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::JanevReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::JanevReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::JanevReactionRate::write(Ostream& os) const
 {
     os.writeKeyword("A") << A_ << nl;
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H
index c353475be..8b1c3572e 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,6 +102,35 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) 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;
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H
index 4d8e92b55..913dbd215 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H
@@ -72,7 +72,7 @@ inline Foam::scalar Foam::LandauTellerReactionRate::operator()
         lta *= pow(T, beta_);
     }
 
-    scalar expArg = 0.0;
+    scalar expArg = 0;
 
     if (mag(Ta_) > vSmall)
     {
@@ -98,6 +98,81 @@ inline Foam::scalar Foam::LandauTellerReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::LandauTellerReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField&
+) const
+{
+    scalar lta = A_;
+
+    if (mag(beta_) > vSmall)
+    {
+        lta *= pow(T, beta_);
+    }
+
+    scalar expArg = 0;
+    scalar deriv = 0;
+
+    if (mag(Ta_) > vSmall)
+    {
+        scalar TaT = Ta_/T;
+        expArg -= TaT;
+        deriv += TaT;
+    }
+
+    if (mag(B_) > vSmall)
+    {
+        scalar BT = B_/cbrt(T);
+        expArg += BT;
+        deriv -= BT/3;
+    }
+
+    if (mag(C_) > vSmall)
+    {
+        scalar CT = C_/pow(T, 2.0/3.0);
+        expArg += CT;
+        deriv -= 2*CT/3;
+    }
+
+    if (mag(expArg) > vSmall)
+    {
+        lta *= exp(expArg);
+    }
+
+    return lta*(beta_+deriv)/T;
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::LandauTellerReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::LandauTellerReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::LandauTellerReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::LandauTellerReactionRate::write(Ostream& os) const
 {
     os.writeKeyword("A") << A_ << token::END_STATEMENT << nl;
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H
index 1d1c5dd33..4b38fc897 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,6 +104,35 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) 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;
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H
index d0b13587b..7429b2d73 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -89,6 +89,77 @@ inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    scalar den =
+    (
+        T
+       *sqr(1 + A_[1]*exp(-Ta_[1]/T)*c[co_] + A_[2]*exp(-Ta_[2]/T)*c[c3h6_])
+       *(1 + A_[3]*exp(-Ta_[3]/T)*sqr(c[co_])*sqr(c[c3h6_]))
+       *(1 + A_[4]*exp(-Ta_[4]/T)*pow(c[no_], 0.7))
+    );
+    scalar rate = A_[0]*exp(-Ta_[0]/T)/ den;
+
+    scalar derivDen =
+    (
+        sqr(1 + A_[1]*exp(-Ta_[1]/T)*c[co_] + A_[2]*exp(-Ta_[2]/T)*c[c3h6_])
+       *(1 + A_[3]*exp(-Ta_[3]/T)*sqr(c[co_])*sqr(c[c3h6_]))
+       *(1 + A_[4]*exp(-Ta_[4]/T)*pow(c[no_], 0.7))
+      + 2*T*(1 + A_[1]*exp(-Ta_[1]/T)*c[co_] + A_[2]*exp(-Ta_[2]/T)*c[c3h6_])
+       *
+        (
+            A_[1]*exp(-Ta_[1]/T)*c[co_]*Ta_[1]/sqr(T)
+          + A_[2]*exp(-Ta_[2]/T)*c[c3h6_]*Ta_[2]/sqr(T)
+        )
+       *(1 + A_[3]*exp(-Ta_[3]/T)*sqr(c[co_])*sqr(c[c3h6_]))
+       *(1 + A_[4]*exp(-Ta_[4]/T)*pow(c[no_], 0.7))
+      + T
+       *sqr(1 + A_[1]*exp(-Ta_[1]/T)*c[co_] + A_[2]*exp(-Ta_[2]/T)*c[c3h6_])
+       *(A_[3]*exp(-Ta_[3]/T)*Ta_[3]*sqr(c[co_])*sqr(c[c3h6_])/sqr(T))
+       *(1 + A_[4]*exp(-Ta_[4]/T)*pow(c[no_], 0.7))
+      + T
+       *sqr(1 + A_[1]*exp(-Ta_[1]/T)*c[co_] + A_[2]*exp(-Ta_[2]/T)*c[c3h6_])
+       *(1 + A_[3]*exp(-Ta_[3]/T)*sqr(c[co_])*sqr(c[c3h6_]))
+       *(A_[4]*exp(-Ta_[4]/T)*Ta_[4]*pow(c[no_], 0.7))/sqr(T)
+    );
+
+    return rate*(Ta_[0]/sqr(T) - derivDen/den);
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::LangmuirHinshelwoodReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::LangmuirHinshelwoodReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::LangmuirHinshelwoodReactionRate::write(Ostream& os) const
 {
     FixedList<Tuple2<scalar, scalar>, n_> coeffs;
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunction.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunction.H
index 76c2f3acc..23409455b 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     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,6 +80,22 @@ public:
             const scalar Pr
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdT,
+            const scalar T
+        ) const;
+
+        inline scalar ddc
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdc,
+            const scalar T
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunctionI.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunctionI.H
index 8421898b6..3ab557fed 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     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,6 +48,30 @@ inline Foam::scalar Foam::LindemannFallOffFunction::operator()
 }
 
 
+inline Foam::scalar Foam::LindemannFallOffFunction::ddT
+(
+    const scalar,
+    const scalar,
+    const scalar,
+    const scalar
+) const
+{
+    return 0;
+}
+
+
+inline Foam::scalar Foam::LindemannFallOffFunction::ddc
+(
+    const scalar,
+    const scalar,
+    const scalar,
+    const scalar
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::LindemannFallOffFunction::write(Ostream& os) const
 {}
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunction.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunction.H
index 2a73db4be..7e099ea0d 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     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,6 +91,22 @@ public:
             const scalar Pr
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdT,
+            const scalar T
+        ) const;
+
+        inline scalar ddc
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdc,
+            const scalar T
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunctionI.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunctionI.H
index 82fdf78c8..58ac126ff 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunctionI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunctionI.H
@@ -60,11 +60,48 @@ inline Foam::scalar Foam::SRIFallOffFunction::operator()
     const scalar Pr
 ) const
 {
-    scalar X = 1.0/(1.0 + sqr(log10(max(Pr, small))));
+    scalar X = 1.0/(1 + sqr(log10(max(Pr, small))));
     return d_*pow(a_*exp(-b_/T) + exp(-T/c_), X)*pow(T, e_);
 }
 
 
+inline Foam::scalar Foam::SRIFallOffFunction::ddT
+(
+    const scalar Pr,
+    const scalar F,
+    const scalar dPrdT,
+    const scalar T
+) 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_))
+        );
+}
+
+
+inline Foam::scalar Foam::SRIFallOffFunction::ddc
+(
+    const scalar Pr,
+    const scalar F,
+    const scalar dPrdc,
+    const scalar T
+) 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_));
+}
+
+
 inline void Foam::SRIFallOffFunction::write(Ostream& os) const
 {
     os.writeKeyword("a") << a_ << token::END_STATEMENT << nl;
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunction.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunction.H
index 18b39038f..4a9025fa4 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     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,6 +91,22 @@ public:
             const scalar Pr
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdT,
+            const scalar T
+        ) const;
+
+        inline scalar ddc
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdc,
+            const scalar T
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunctionI.H b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunctionI.H
index 3f009702c..79cc84ed0 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunctionI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunctionI.H
@@ -71,7 +71,102 @@ inline Foam::scalar Foam::TroeFallOffFunction::operator()
     scalar n = 0.75 - 1.27*logFcent;
 
     scalar logPr = log10(max(Pr, small));
-    return pow(10.0, logFcent/(1.0 + sqr((logPr + c)/(n - d*(logPr + c)))));
+    return pow(10, logFcent/(1 + sqr((logPr + c)/(n - d*(logPr + c)))));
+}
+
+
+inline Foam::scalar Foam::TroeFallOffFunction::ddT
+(
+    const scalar Pr,
+    const scalar F,
+    const scalar dPrdT,
+    const scalar T
+) 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);
+
+    scalar dFcentdT =
+    (
+        (alpha_ - 1)*exp(-T/Tsss_)/Tsss_
+      - alpha_*exp(-T/Ts_)/Ts_
+      + Tss_*exp(-Tss_/T)/sqr(T)
+    );
+
+    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;
+
+    scalar dlogPrdT = dPrdT/Pr/logTen;
+
+    scalar dParentdT =
+        2.0*(logPr + c)/sqr(n - d*(logPr + c))
+       *(
+            (dlogPrdT + dcdT)
+          - (logPr + c)*(dndT - d*(dlogPrdT + dcdT))/(n - d*(logPr + c))
+        );
+
+    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))))
+        )
+    );
+}
+
+
+inline Foam::scalar Foam::TroeFallOffFunction::ddc
+(
+    const scalar Pr,
+    const scalar F,
+    const scalar dPrdc,
+    const scalar T
+) 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
+        )
+    );
+
+    scalar dlogPrdc = dPrdc/Pr/logTen;
+    scalar d = 0.14;
+    scalar c = -0.4 - 0.67*logFcent;
+    scalar n = 0.75 - 1.27*logFcent;
+
+    scalar dParentdc =
+        2.0*(logPr + c)/sqr(n - d*(logPr + c))
+       *(
+            (dlogPrdc)
+          - (logPr + c)*(-d*(dlogPrdc))/(n - d*(logPr + c))
+        );
+
+    return
+    (
+        F*logTen
+       *(
+          - logFcent*dParentdc/sqr(1.0 + sqr((logPr + c)/(n - d*(logPr + c))))
+        )
+    );
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRate.H
index ec8abe2cd..0dc9f13b2 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -87,6 +87,35 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) 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;
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRateI.H
index 918603451..137cf0ec9 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRateI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,6 +53,44 @@ inline Foam::scalar Foam::infiniteReactionRate::operator()
     return (1);
 }
 
+inline Foam::scalar Foam::infiniteReactionRate::ddT
+(
+    const scalar p,
+    const scalar,
+    const scalarField&
+) const
+{
+    return (0);
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::infiniteReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::infiniteReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::infiniteReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
 
 inline Foam::Ostream& Foam::operator<<
 (
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H
index f12d3ba0d..4187365aa 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -103,6 +103,35 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) 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;
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H
index f727c8866..446b29886 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H
@@ -69,7 +69,7 @@ inline Foam::scalar Foam::powerSeriesReactionRate::operator()
         lta *= pow(T, beta_);
     }
 
-    scalar expArg = 0.0;
+    scalar expArg = 0;
 
     forAll(coeffs_, n)
     {
@@ -82,6 +82,64 @@ inline Foam::scalar Foam::powerSeriesReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::powerSeriesReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField&
+) const
+{
+    scalar lta = A_;
+
+    if (mag(beta_) > vSmall)
+    {
+        lta *= pow(T, beta_);
+    }
+
+    scalar expArg = 0;
+    scalar deriv = 0;
+
+    forAll(coeffs_, n)
+    {
+        scalar cT = coeffs_[n]/pow(T, n + 1);
+        expArg += cT;
+        deriv -= (n + 1)*cT;
+    }
+
+    lta *= exp(expArg);
+
+    return lta*(beta_+deriv)/T;
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::powerSeriesReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::powerSeriesReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::powerSeriesReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::powerSeriesReactionRate::write(Ostream& os) const
 {
     os.writeKeyword("A") << A_ << token::END_STATEMENT << nl;
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H
index 3c7375d73..9e7bae406 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H
@@ -61,6 +61,7 @@ class thirdBodyArrheniusReactionRate
     // Private data
 
         thirdBodyEfficiencies thirdBodyEfficiencies_;
+        List<Tuple2<label, scalar>> beta_;
 
 
 public:
@@ -99,6 +100,38 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) 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_;
+        }
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H
index d9b39b61a..7061aeea4 100644
--- a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H
+++ b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,12 @@ 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
@@ -50,7 +55,19 @@ inline Foam::thirdBodyArrheniusReactionRate::thirdBodyArrheniusReactionRate
         dict
     ),
     thirdBodyEfficiencies_(species, dict)
-{}
+{
+    forAll(thirdBodyEfficiencies_, i)
+    {
+        beta_.append
+        (
+            Tuple2<label, scalar>
+            (
+                i,
+                thirdBodyEfficiencies_[i]
+            )
+        );
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -68,6 +85,46 @@ inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return
+        thirdBodyEfficiencies_.M(c)
+       *ArrheniusReactionRate::ddT(p, T, c);
+}
+
+
+inline void Foam::thirdBodyArrheniusReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    scalar M = thirdBodyEfficiencies_.M(c);
+    forAll(beta_, i)
+    {
+        dcidc[i] = beta_[i].second()/max(M, small);
+    }
+}
+
+
+inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return -1.0/T;
+}
+
+
 inline void Foam::thirdBodyArrheniusReactionRate::write(Ostream& os) const
 {
     ArrheniusReactionRate::write(os);
diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyEfficiencies/thirdBodyEfficienciesI.H b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyEfficiencies/thirdBodyEfficienciesI.H
index 97a0e8eee..a1aa8265b 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     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -83,7 +83,7 @@ inline Foam::thirdBodyEfficiencies::thirdBodyEfficiencies
 
 inline Foam::scalar Foam::thirdBodyEfficiencies::M(const scalarList& c) const
 {
-    scalar M = 0.0;
+    scalar M = 0;
     forAll(*this, i)
     {
         M += operator[](i)*c[i];
diff --git a/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H b/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H
index 2a9c599ff..b212d93ec 100644
--- a/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H
+++ b/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -155,6 +155,16 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
+
         // I-O
 
             //- Write to Ostream
diff --git a/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H b/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H
index f913d48f9..c2522b34f 100644
--- a/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H
@@ -139,6 +139,28 @@ inline Foam::scalar Foam::eConstThermo<EquationOfState>::S
 }
 
 
+template<class EquationOfState>
+inline Foam::scalar Foam::eConstThermo<EquationOfState>::dGdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return 0;
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::eConstThermo<EquationOfState>::dCpdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return 0;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState>
diff --git a/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H b/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H
index dadbb0712..5a5c7626e 100644
--- a/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H
+++ b/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -148,6 +148,15 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
         // I-O
 
             //- Write to Ostream
diff --git a/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H b/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H
index cf75c27e2..5c48335b3 100644
--- a/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H
@@ -136,6 +136,26 @@ inline Foam::scalar Foam::hConstThermo<EquationOfState>::S
 }
 
 
+template<class EquationOfState>
+inline Foam::scalar Foam::hConstThermo<EquationOfState>::dGdT
+(
+    const scalar p, const scalar T
+) const
+{
+    return 0;
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::hConstThermo<EquationOfState>::dCpdT
+(
+    const scalar p, const scalar T
+) const
+{
+    return 0;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState>
diff --git a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H
index f9b73bd93..c9b5adb60 100644
--- a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H
+++ b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -198,6 +198,15 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
         // I-O
 
             //- Write to Ostream
diff --git a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
index 370b3502e..84a96ed19 100644
--- a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
@@ -126,6 +126,36 @@ inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::S
 }
 
 
+template<class EquationOfState, int PolySize>
+inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::dGdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return
+    (
+        hCoeffs_.derivative(T)
+      - T*sCoeffs_.derivative(T)
+      - sCoeffs_.value(T)
+    );
+}
+
+
+template<class EquationOfState, int PolySize>
+inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::dCpdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return
+    (
+        CpCoeffs_.derivative(T)
+    );
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState, int PolySize>
diff --git a/src/thermophysicalModels/specie/thermo/hPower/hPowerThermo.H b/src/thermophysicalModels/specie/thermo/hPower/hPowerThermo.H
index 7ca3b2bbd..cf19e37e4 100644
--- a/src/thermophysicalModels/specie/thermo/hPower/hPowerThermo.H
+++ b/src/thermophysicalModels/specie/thermo/hPower/hPowerThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -169,6 +169,15 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
     // Member operators
 
         inline void operator+=(const hPowerThermo&);
diff --git a/src/thermophysicalModels/specie/thermo/hPower/hPowerThermoI.H b/src/thermophysicalModels/specie/thermo/hPower/hPowerThermoI.H
index 2b2f61016..115a9af7f 100644
--- a/src/thermophysicalModels/specie/thermo/hPower/hPowerThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/hPower/hPowerThermoI.H
@@ -164,6 +164,30 @@ inline Foam::scalar Foam::hPowerThermo<EquationOfState>::S
 }
 
 
+template<class EquationOfState>
+inline Foam::scalar Foam::hPowerThermo<EquationOfState>::dGdT
+(
+    const scalar p, const scalar T
+) const
+{
+    // To be implemented
+    NotImplemented;
+    return 0;
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::hPowerThermo<EquationOfState>::dCpdT
+(
+    const scalar p, const scalar T
+) const
+{
+    // To be implemented
+    NotImplemented;
+    return 0;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState>
diff --git a/src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermo.H b/src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermo.H
index 5af026540..a139fa564 100644
--- a/src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermo.H
+++ b/src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -152,6 +152,15 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
         // I-O
 
             //- Write to Ostream
diff --git a/src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermoI.H b/src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermoI.H
index e5ebfd954..15e7988ca 100644
--- a/src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermoI.H
@@ -142,6 +142,26 @@ inline Foam::scalar Foam::hRefConstThermo<EquationOfState>::S
 }
 
 
+template<class EquationOfState>
+inline Foam::scalar Foam::hRefConstThermo<EquationOfState>::dGdT
+(
+    const scalar p, const scalar T
+) const
+{
+    return 0;
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::hRefConstThermo<EquationOfState>::dCpdT
+(
+    const scalar p, const scalar T
+) const
+{
+    return 0;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState>
diff --git a/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H b/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H
index f16b76ebf..8b487c38c 100644
--- a/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H
+++ b/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -187,6 +187,15 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
         // I-O
 
             //- Write to Ostream
diff --git a/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H b/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
index ad8b1a9fa..7391c4b30 100644
--- a/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
@@ -237,6 +237,31 @@ inline Foam::scalar Foam::janafThermo<EquationOfState>::S
 }
 
 
+template<class EquationOfState>
+inline Foam::scalar Foam::janafThermo<EquationOfState>::dGdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    const coeffArray& a = coeffs(T);
+    return -((a[0] + a[5]/T)/T + a[1]/2 + T*(a[2]/3 + T*(a[3]/4 + T*a[4]/5)));
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::janafThermo<EquationOfState>::dCpdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    const coeffArray& a = coeffs(T);
+    return
+        (((4*a[4]*T + 3*a[3])*T + 2*a[2])*T + a[1]);
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState>
diff --git a/src/thermophysicalModels/specie/thermo/thermo/thermo.H b/src/thermophysicalModels/specie/thermo/thermo/thermo.H
index a5cdba6ec..3a644d5d8 100644
--- a/src/thermophysicalModels/specie/thermo/thermo/thermo.H
+++ b/src/thermophysicalModels/specie/thermo/thermo/thermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::thermo
+    Foam::species::thermo
 
 Description
     Basic thermodynamics type based on the use of fitting functions for
@@ -318,6 +318,15 @@ public:
             ) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of B (acooding to Niemeyer et al.) w.r.t. temperature
+            inline scalar dKcdTbyKc(const scalar p, const scalar T) const;
+
+            //- Derivative of cp w.r.t. temperature
+            inline scalar dcpdT(const scalar p, const scalar T) const;
+
+
         // I-O
 
             //- Write to Ostream
diff --git a/src/thermophysicalModels/specie/thermo/thermo/thermoI.H b/src/thermophysicalModels/specie/thermo/thermo/thermoI.H
index 791042fe7..0fced0cee 100644
--- a/src/thermophysicalModels/specie/thermo/thermo/thermoI.H
+++ b/src/thermophysicalModels/specie/thermo/thermo/thermoI.H
@@ -444,6 +444,35 @@ inline Foam::scalar Foam::species::thermo<Thermo, Type>::TEa
 }
 
 
+template<class Thermo, template<class> class Type>
+inline Foam::scalar
+Foam::species::thermo<Thermo, Type>::dKcdTbyKc
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    const scalar nm = this->Y()/this->W();
+
+    if (equal(nm, small))
+    {
+        return -this->dGdT(Pstd, T)*this->Y()/RR;
+    }
+    else
+    {
+        return -(nm/T + this->dGdT(Pstd, T)*this->Y()/RR);
+    }
+}
+
+
+template<class Thermo, template<class> class Type>
+inline Foam::scalar
+Foam::species::thermo<Thermo, Type>::dcpdT(const scalar p, const scalar T) const
+{
+    return this->dCpdT(p, T)*this->W();;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class Thermo, template<class> class Type>
diff --git a/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties b/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties
index 2dc50a724..9c75d7dff 100644
--- a/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties
+++ b/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties
@@ -33,7 +33,7 @@ EulerImplicitCoeffs
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-14;
+    absTol          1e-12;
     relTol          1e-1;
 }
 
diff --git a/tutorials/combustion/fireFoam/LES/oppositeBurningPanels/constant/panelRegion/chemistryProperties b/tutorials/combustion/fireFoam/LES/oppositeBurningPanels/constant/panelRegion/chemistryProperties
index df407c07f..46997d9d2 100644
--- a/tutorials/combustion/fireFoam/LES/oppositeBurningPanels/constant/panelRegion/chemistryProperties
+++ b/tutorials/combustion/fireFoam/LES/oppositeBurningPanels/constant/panelRegion/chemistryProperties
@@ -28,7 +28,7 @@ initialChemicalTimeStep 1e-07;
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          0.01;
 }
 
diff --git a/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/chemistryProperties b/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/chemistryProperties
index c250d3841..aafbdf700 100644
--- a/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/chemistryProperties
@@ -35,9 +35,9 @@ initialChemicalTimeStep 1e-07;
 
 odeCoeffs
 {
-    solver          Rosenbrock34;  // Rosenbrock34, seulex or rodas23
-    absTol          1e-12;
-    relTol          1e-7;
+    solver          seulex;
+    absTol          1e-08;
+    relTol          0.1;
 }
 
 reduction
@@ -86,7 +86,7 @@ tabulation
     printNumRetrieve   off;
 
     // Tolerance used for retrieve and grow
-    tolerance   1e-4;
+    tolerance   0.003;
 
     // ISAT is the only method currently available
     method    ISAT;
@@ -97,24 +97,24 @@ tabulation
         otherSpecies 1;
         Temperature  1000;
         Pressure     1e15;
-        deltaT       0.5;
+        deltaT       1;
     }
 
     // Maximum number of leafs stored in the binary tree
-    maxNLeafs  2000;
+    maxNLeafs  5000;
 
     // Maximum life time of the leafs (in time steps) used in unsteady
     // simulations to force renewal of the stored chemPoints and keep the tree
     // small
-    chPMaxLifeTime 100;
+    chPMaxLifeTime 1000;
 
     // Maximum number of growth allowed on a chemPoint to avoid distorted
     // chemPoints
-    maxGrowth  10;
+    maxGrowth  100;
 
     // Number of time steps between analysis of the tree to remove old
     // chemPoints or try to balance it
-    checkEntireTreeInterval  5;
+    checkEntireTreeInterval  500;
 
     // Parameters used to decide whether to balance or not if the tree's depth
     // is larger than maxDepthFactor*log2(nLeafs) then balance the tree
@@ -123,7 +123,6 @@ tabulation
     // Try to balance the tree only if the size of the tree is greater
     minBalanceThreshold 30;
 
-
     // Activate the use of a MRU (most recently used) list
     MRURetrieve false;
 
diff --git a/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/chemistryProperties.test b/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/chemistryProperties.test
deleted file mode 100644
index cb5c5b610..000000000
--- a/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/chemistryProperties.test
+++ /dev/null
@@ -1,142 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      chemistryProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-chemistryType
-{
-    solver            ode;
-    method            TDAC;
-}
-
-chemistry           on;
-
-importantSpecies
-{
-    CH4;
-    H2O;
-    O2;
-    CO2;
-}
-
-initialChemicalTimeStep 1e-07;
-
-odeCoeffs
-{
-    solver          Rosenbrock34;  // Rosenbrock34, seulex or rodas23
-    absTol          1e-12;
-    relTol          1e-7;
-}
-
-reduction
-{
-    // Activate reduction
-    active  on;
-
-    // Switch logging of the reduction statistics and performance
-    log         on;
-
-    // Tolerance depends on the reduction method, see details for each method
-    tolerance   1e-4;
-
-    // Available methods: DRG, DAC, DRGEP, PFA, EFA
-    method DAC;
-
-    // Search initiating set (SIS) of species, needed for most methods
-    initialSet
-    {
-        CO;
-        CH4;
-        HO2;
-    }
-
-    // For DAC, option to automatically change the SIS switch from HO2 to H2O
-    // and CO to CO2, + disable fuel
-    automaticSIS    off;
-
-    // When automaticSIS, the method needs to know the fuel
-    fuelSpecies
-    {
-        CH4 1;
-    }
-}
-
-tabulation
-{
-    // Activate tabulation
-    active      on;
-
-    // Switch logging of the tabulation statistics and performance
-    log         on;
-
-    printProportion    off;
-
-    printNumRetrieve   off;
-
-    // Tolerance used for retrieve and grow
-    tolerance   0.003;
-
-    // ISAT is the only method currently available
-    method    ISAT;
-
-    // Scale factors used in the definition of the ellipsoid of accuracy
-    scaleFactor
-    {
-        otherSpecies 1;
-        Temperature  10000;
-        Pressure     1e15;
-        deltaT       1;
-    }
-
-    // Maximum number of leafs stored in the binary tree
-    maxNLeafs  5000;
-
-    // Maximum life time of the leafs (in time steps) used in unsteady
-    // simulations to force renewal of the stored chemPoints and keep the tree
-    // small
-    chPMaxLifeTime 1000;
-
-    // Maximum number of growth allowed on a chemPoint to avoid distorted
-    // chemPoints
-    maxGrowth  100;
-
-    // Number of time steps between analysis of the tree to remove old
-    // chemPoints or try to balance it
-    checkEntireTreeInterval  500;
-
-    // Parameters used to decide whether to balance or not if the tree's depth
-    // is larger than maxDepthFactor*log2(nLeafs) then balance the tree
-    maxDepthFactor   2;
-
-    // Try to balance the tree only if the size of the tree is greater
-    minBalanceThreshold 30;
-
-    // Activate the use of a MRU (most recently used) list
-    MRURetrieve false;
-
-    // Maximum size of the MRU list
-    maxMRUSize 0;
-
-    // Allow to grow points
-    growPoints  true;
-
-    // When mechanism reduction is used, new dimensions might be added
-    // maxNumNewDim set the maximum number of new dimensions added during a
-    // growth
-    maxNumNewDim 10;
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/reactionsGRI b/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/reactionsGRI
index 8bd6be358..9e89ab5d9 100644
--- a/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/reactionsGRI
+++ b/tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/reactionsGRI
@@ -3608,3 +3608,6 @@ reactions
         Ta              0;
     }
 }
+
+Tlow            250;
+Thigh           5000;
diff --git a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/Allrun b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/Allrun
index 2769119da..d8137ce0a 100755
--- a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/Allrun
+++ b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/Allrun
@@ -14,16 +14,16 @@ runApplication chemkinToFoam \
 runApplication blockMesh
 runApplication setFields
 
-
 # Run the application without chemistry until 1500 to let the flow field develop
+foamDictionary -entry "startTime" -set "0" system/controlDict
 foamDictionary -entry "writeInterval" -set "1500" system/controlDict
 foamDictionary -entry "endTime" -set "1500" system/controlDict
 foamDictionary -entry "chemistry" -set "off" constant/chemistryProperties
 
 runApplication $application
 
-
 # Run with chemistry until flame reach its full size
+foamDictionary -entry "startTime" -set "1500" system/controlDict
 foamDictionary -entry "writeInterval" -set "100" system/controlDict
 foamDictionary -entry "endTime" -set "5000" system/controlDict
 foamDictionary -entry "chemistry" -set "on" constant/chemistryProperties
diff --git a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/chemistryProperties b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/chemistryProperties
index f8c1f4c3d..cf015f487 100644
--- a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/chemistryProperties
@@ -20,7 +20,7 @@ chemistryType
     method          TDAC;
 }
 
-chemistry       off;
+chemistry       on;
 
 importantSpecies
 {
@@ -35,8 +35,8 @@ initialChemicalTimeStep 1e-07;
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
-    relTol          1e-07;
+    absTol          1e-08;
+    relTol          0.1;
 }
 
 reduction
diff --git a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/reactionsGRI b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/reactionsGRI
index 8bd6be358..9e89ab5d9 100644
--- a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/reactionsGRI
+++ b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/reactionsGRI
@@ -3608,3 +3608,6 @@ reactions
         Ta              0;
     }
 }
+
+Tlow            250;
+Thigh           5000;
diff --git a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/system/controlDict b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/system/controlDict
index bf9887f7b..774eb04d0 100644
--- a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/system/controlDict
+++ b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/system/controlDict
@@ -16,19 +16,19 @@ FoamFile
 
 application     reactingFoam;
 
-startFrom       latestTime;
+startFrom       startTime;
 
-startTime       0;
+startTime       1500;
 
 stopAt          endTime;
 
-endTime         1500;
+endTime         5000;
 
 deltaT          1;
 
 writeControl    runTime;
 
-writeInterval   1500;
+writeInterval   100;
 
 purgeWrite      0;
 
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties
index 3812e6131..f49edd526 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties
@@ -33,7 +33,7 @@ EulerImplicitCoeffs
 odeCoeffs
 {
     solver          Rosenbrock43;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          0.01;
 }
 
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS/constant/chemistryProperties b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS/constant/chemistryProperties
index 3812e6131..f49edd526 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS/constant/chemistryProperties
@@ -33,7 +33,7 @@ EulerImplicitCoeffs
 odeCoeffs
 {
     solver          Rosenbrock43;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          0.01;
 }
 
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/chemistryProperties b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/chemistryProperties
index 2604f1fd8..9fd9f1912 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/chemistryProperties
@@ -24,11 +24,12 @@ chemistryType
 chemistry       on;
 
 initialChemicalTimeStep 1e-7;
+//maxChemicalTimeStep 1E-3;
 
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          1e-1;
 }
 
@@ -78,7 +79,7 @@ tabulation
     printNumRetrieve   off;
 
     // Tolerance used for retrieve and grow
-    tolerance   1e-3;
+    tolerance   3e-3;
 
     // ISAT is the only method currently available
     method    ISAT;
@@ -87,7 +88,7 @@ tabulation
     scaleFactor
     {
         otherSpecies 1;
-        Temperature  25000;
+        Temperature  10000;
         Pressure     1e15;
         deltaT       1;
     }
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/reactionsGRI b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/reactionsGRI
index de8ff0a51..694ea19b6 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/reactionsGRI
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/reactionsGRI
@@ -5588,3 +5588,6 @@ reactions
         }
     }
 }
+
+Thigh   5000;
+Tlow    200;
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/system/controlDict b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/system/controlDict
index 792dbf249..0aff953b0 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/system/controlDict
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/system/controlDict
@@ -23,7 +23,7 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         1000;
+endTime         1500;
 
 deltaT          1;
 
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/chemistryProperties b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/chemistryProperties
index d63f983b7..678d5e60f 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/chemistryProperties
@@ -7,29 +7,28 @@
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      chemistryProperties;
+    version         2;
+    format          ascii;
+    class           dictionary;
+    location        "constant";
+    object          chemistryProperties;
 }
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 chemistryType
 {
-    solver            ode;
+    solver          ode;
 }
 
 chemistry       on;
 
-initialChemicalTimeStep 1e-7;
-maxChemicalTimeStep     1e-4;
+initialChemicalTimeStep 1e-07;
 
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
-    relTol          1e-1;
+    absTol          1e-8;
+    relTol          0.1;
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/reactionsGRI b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/reactionsGRI
index 276111048..694ea19b6 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/reactionsGRI
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/reactionsGRI
@@ -2349,7 +2349,7 @@ reactions
         A               1.2e+11;
         beta            -1;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2415,7 +2415,7 @@ reactions
         A               5e+11;
         beta            -1;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2481,7 +2481,7 @@ reactions
         A               2.8e+12;
         beta            -0.86;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 1)
@@ -2547,7 +2547,7 @@ reactions
         A               1e+12;
         beta            -1;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2613,7 +2613,7 @@ reactions
         A               2.2e+16;
         beta            -2;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2679,7 +2679,7 @@ reactions
         A               1.87e+14;
         beta            -1;
         Ta              8554.25;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2745,7 +2745,7 @@ reactions
         A               1.06e+14;
         beta            -1.41;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2811,7 +2811,7 @@ reactions
         A               1.3e+11;
         beta            -0.11;
         Ta              2505.89;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2877,7 +2877,7 @@ reactions
         A               4.48e+13;
         beta            -1.32;
         Ta              372.362;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2943,7 +2943,7 @@ reactions
         A               3.1e+11;
         beta            0;
         Ta              27197.5;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -3009,7 +3009,7 @@ reactions
         A               1.04e+26;
         beta            -3.3;
         Ta              63704;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -3075,7 +3075,7 @@ reactions
         A               1.18e+13;
         beta            0;
         Ta              42630.4;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -3155,7 +3155,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3240,7 +3240,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3325,7 +3325,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 3)
@@ -3410,7 +3410,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3495,7 +3495,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3580,7 +3580,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3665,7 +3665,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3750,7 +3750,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3835,7 +3835,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3920,7 +3920,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4005,7 +4005,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4090,7 +4090,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4175,7 +4175,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4260,7 +4260,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4345,7 +4345,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4430,7 +4430,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4515,7 +4515,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4600,7 +4600,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4685,7 +4685,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4770,7 +4770,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4851,7 +4851,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4932,7 +4932,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5017,7 +5017,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5102,7 +5102,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5187,7 +5187,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5272,7 +5272,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5357,7 +5357,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5442,7 +5442,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5527,7 +5527,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5588,3 +5588,6 @@ reactions
         }
     }
 }
+
+Thigh   5000;
+Tlow    200;
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties
index 3a230eb6f..9fd9f1912 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties
@@ -24,12 +24,12 @@ chemistryType
 chemistry       on;
 
 initialChemicalTimeStep 1e-7;
-maxChemicalTimeStep     1e-4;
+//maxChemicalTimeStep 1E-3;
 
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          1e-1;
 }
 
@@ -79,7 +79,7 @@ tabulation
     printNumRetrieve   off;
 
     // Tolerance used for retrieve and grow
-    tolerance   1e-4;
+    tolerance   3e-3;
 
     // ISAT is the only method currently available
     method    ISAT;
@@ -88,7 +88,7 @@ tabulation
     scaleFactor
     {
         otherSpecies 1;
-        Temperature  2500;
+        Temperature  10000;
         Pressure     1e15;
         deltaT       1;
     }
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties.new b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties.new
index 0b4586f51..330d01023 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties.new
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties.new
@@ -28,7 +28,7 @@ initialChemicalTimeStep 1e-7;
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          1e-1;
 }
 
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/reactionsGRI b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/reactionsGRI
index 276111048..694ea19b6 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/reactionsGRI
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/reactionsGRI
@@ -2349,7 +2349,7 @@ reactions
         A               1.2e+11;
         beta            -1;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2415,7 +2415,7 @@ reactions
         A               5e+11;
         beta            -1;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2481,7 +2481,7 @@ reactions
         A               2.8e+12;
         beta            -0.86;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 1)
@@ -2547,7 +2547,7 @@ reactions
         A               1e+12;
         beta            -1;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2613,7 +2613,7 @@ reactions
         A               2.2e+16;
         beta            -2;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2679,7 +2679,7 @@ reactions
         A               1.87e+14;
         beta            -1;
         Ta              8554.25;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2745,7 +2745,7 @@ reactions
         A               1.06e+14;
         beta            -1.41;
         Ta              0;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2811,7 +2811,7 @@ reactions
         A               1.3e+11;
         beta            -0.11;
         Ta              2505.89;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2877,7 +2877,7 @@ reactions
         A               4.48e+13;
         beta            -1.32;
         Ta              372.362;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -2943,7 +2943,7 @@ reactions
         A               3.1e+11;
         beta            0;
         Ta              27197.5;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -3009,7 +3009,7 @@ reactions
         A               1.04e+26;
         beta            -3.3;
         Ta              63704;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -3075,7 +3075,7 @@ reactions
         A               1.18e+13;
         beta            0;
         Ta              42630.4;
-        coeffs          
+        coeffs
 53
 (
 (CH4 2)
@@ -3155,7 +3155,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3240,7 +3240,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3325,7 +3325,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 3)
@@ -3410,7 +3410,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3495,7 +3495,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3580,7 +3580,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3665,7 +3665,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3750,7 +3750,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3835,7 +3835,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -3920,7 +3920,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4005,7 +4005,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4090,7 +4090,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4175,7 +4175,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4260,7 +4260,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4345,7 +4345,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4430,7 +4430,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4515,7 +4515,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4600,7 +4600,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4685,7 +4685,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4770,7 +4770,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4851,7 +4851,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -4932,7 +4932,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5017,7 +5017,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5102,7 +5102,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5187,7 +5187,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5272,7 +5272,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5357,7 +5357,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5442,7 +5442,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5527,7 +5527,7 @@ reactions
         }
         thirdBodyEfficiencies
         {
-            coeffs          
+            coeffs
 53
 (
 (CH4 2)
@@ -5588,3 +5588,6 @@ reactions
         }
     }
 }
+
+Thigh   5000;
+Tlow    200;
diff --git a/tutorials/lagrangian/reactingParcelFoam/counterFlowFlame2DLTS/constant/chemistryProperties b/tutorials/lagrangian/reactingParcelFoam/counterFlowFlame2DLTS/constant/chemistryProperties
index 3812e6131..f49edd526 100644
--- a/tutorials/lagrangian/reactingParcelFoam/counterFlowFlame2DLTS/constant/chemistryProperties
+++ b/tutorials/lagrangian/reactingParcelFoam/counterFlowFlame2DLTS/constant/chemistryProperties
@@ -33,7 +33,7 @@ EulerImplicitCoeffs
 odeCoeffs
 {
     solver          Rosenbrock43;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          0.01;
 }
 
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnEvaporatingReacting/constant/chemistryProperties.gas b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnEvaporatingReacting/constant/chemistryProperties.gas
index d328270d6..d27415ed4 100644
--- a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnEvaporatingReacting/constant/chemistryProperties.gas
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnEvaporatingReacting/constant/chemistryProperties.gas
@@ -33,7 +33,7 @@ EulerImplicitCoeffs
 odeCoeffs
 {
     solver          Rosenbrock43;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          0.01;
 }