diff --git a/src/common/logMath.hpp b/src/common/logMath.hpp
new file mode 100644
index 0000000..d3b4d7c
--- /dev/null
+++ b/src/common/logMath.hpp
@@ -0,0 +1,305 @@
+#pragma once
+
+#include "macros.hpp"
+#include <cmath>
+
+/**
+ * @file logMath.hpp
+ * @brief Base-selectable logarithm and exponential helpers with optional fast CUDA intrinsics.
+ *
+ * @details
+ * This header provides a small, header-only utility for computing `log` and `exp` in either:
+ * - natural base (e), or
+ * - base-10 (ten),
+ * and optionally using CUDA "fast math" intrinsics on device code.
+ *
+ * Key points:
+ * - For `float` on CUDA device (`__CUDA_ARCH__`), `MathMode::fast` selects intrinsics such as
+ *   `__logf`, `__log10f`, `__expf`, `__exp10f` when available.
+ * - For host code (and for non-fast mode), it falls back to standard `<cmath>` functions
+ *   like `::log`, `::log10`, `::exp`, `::expf`, `::logf`, `::log10f`.
+ * - Some libcs (glibc) provide `exp10/exp10f` as extensions. When unavailable, base-10
+ *   exponentiation is implemented via `exp(x * ln(10))`.
+ *
+ * @note
+ * The functions here assume their usual mathematical domains:
+ * - `log(x)` requires `x > 0`.
+ * - `exp(x)` is defined for all real `x` but may overflow for large `x`.
+ *
+ * @warning
+ * When `MathMode::fast` is used on CUDA device, results may differ slightly from the
+ * accurate host/standard-library implementations due to reduced precision/approximations.
+ */
+
+namespace hpcReact
+{
+
+/**
+ * @def HPCREACT_HAS_EXP10
+ * @brief Indicates whether `exp10/exp10f` are available as libc extensions.
+ *
+ * @details
+ * glibc commonly provides `exp10` and `exp10f` as non-standard extensions. This macro is used to
+ * select a direct call to those functions on host when available; otherwise a mathematically
+ * equivalent fallback is used:
+ * \f[
+ *   10^x = e^{x \ln(10)}
+ * \f]
+ */
+#if defined(__GLIBC__)  // glibc provides exp10/exp10f as an extension
+  #define HPCREACT_HAS_EXP10 1
+#else
+  #define HPCREACT_HAS_EXP10 0
+#endif
+
+/**
+ * @brief Enumerates the logarithm base used by the helper.
+ */
+enum class LogBase
+{
+  e,   /**< Natural logarithm base (ln). */
+  ten  /**< Base-10 logarithm. */
+};
+
+/**
+ * @brief Enumerates the math evaluation mode.
+ *
+ * @details
+ * `fast` is only used for `float` on CUDA device code. On host (or for `double`), the
+ * standard library functions are used.
+ */
+enum class MathMode
+{
+  accurate, /**< Prefer standard-library implementations. */
+  fast      /**< Prefer CUDA fast intrinsics for `float` on device when available. */
+};
+
+/**
+ * @brief Base-selectable logarithm and exponential functions.
+ *
+ * @tparam Base Selects the logarithm/exp base:
+ *   - `LogBase::e` for natural base
+ *   - `LogBase::ten` for base-10
+ * @tparam Mode Selects evaluation mode. See @ref MathMode.
+ *
+ * @details
+ * Provides:
+ * - `log(x)` returning \f$\ln(x)\f$ or \f$\log_{10}(x)\f$ depending on `Base`
+ * - `exp(x)` returning \f$e^x\f$ or \f$10^x\f$ depending on `Base`
+ *
+ * The overload set covers `float` and `double`.
+ */
+template< LogBase Base, MathMode Mode = MathMode::accurate >
+struct LogExp
+{
+  /**
+   * @brief Compile-time constant for \f$\ln(10)\f$.
+   *
+   * @tparam T Floating-point type (typically `float` or `double`).
+   * @return \f$\ln(10)\f$ as type `T`.
+   *
+   * @details
+   * The literal is provided with long-double precision and cast to `T`.
+   */
+  template< class T >
+  HPCREACT_HOST_DEVICE
+  static constexpr T ln10()
+  {
+    return T(2.3025850929940456840179914546843642L);
+  }
+
+  /**
+   * @brief Conversion factor for derivatives with respect to `log(x)` vs `ln(x)`.
+   *
+   * @tparam T Floating-point type (typically `float` or `double`).
+   * @return A constant factor depending on @p Base.
+   *
+   * @details
+   * This returns the scalar \f$c\f$ such that:
+   * \f[
+   *   \frac{d}{d(\log_b x)} = c \, \frac{d}{d(\ln x)}
+   * \f]
+   * where:
+   * - if `Base == e`, then \f$c = 1\f$
+   * - if `Base == ten`, then \f$c = \ln(10)\f$
+   *
+   * Equivalently, since \f$\log_{10}(x) = \ln(x) / \ln(10)\f$:
+   * \f[
+   *   \frac{d}{d(\log_{10} x)} = \ln(10) \, \frac{d}{d(\ln x)}
+   * \f]
+   */
+  template< typename T >
+  HPCREACT_HOST_DEVICE
+  static constexpr inline T dWrtLogConst()
+  {
+    if constexpr ( Base == LogBase::e ) { return T(1.0); }
+    else                                { return T(ln10<T>()); }
+  }
+
+
+  /**
+   * @brief Adjusts a derivative value to account for logarithm base.
+   *
+   * @tparam T Floating-point type (typically `float` or `double`).
+   * @param[in,out] x The derivative value to adjust.
+   * @return The adjusted derivative value.
+   *
+   * @details
+   * This function multiplies the input `x` by the appropriate factor so that
+   * it represents a derivative with respect to `log_{10}(x)` instead of `ln(x)`.
+   * - If `Base == e`, `x` is unchanged.
+   * - If `Base == ten`, `x` is multiplied by \f$\ln(10)\f$.
+   */
+  template< typename T >
+  HPCREACT_HOST_DEVICE
+  static constexpr inline T dWrtLogScale( T x )
+  {
+    if constexpr ( Base == LogBase::ten ) { return x * ln10<T>(); }
+    else                                  { return x; }
+  }
+
+  // ***** log function *********************************************************************************************
+
+  /**
+   * @brief Logarithm in the selected base for `double`.
+   *
+   * @param x Input value (must be > 0).
+   * @return `ln(x)` if `Base == LogBase::e`, otherwise `log10(x)`.
+   */
+  HPCREACT_HOST_DEVICE
+  static constexpr inline double log( double const x )
+  {
+    if constexpr ( Base == LogBase::e ) { return ::log( x ); }
+    else                                { return ::log10( x ); }
+  }
+
+  /**
+   * @brief Logarithm in the selected base for `float`.
+   *
+   * @param x Input value (must be > 0).
+   * @return `ln(x)` if `Base == LogBase::e`, otherwise `log10(x)`.
+   *
+   * @details
+   * On CUDA device code (`__CUDA_ARCH__`) and when `Mode == MathMode::fast`,
+   * uses CUDA intrinsics:
+   * - `__logf` / `__log10f`
+   * Otherwise falls back to `<cmath>`:
+   * - `::logf` / `::log10f`
+   */
+  HPCREACT_HOST_DEVICE
+  static inline float log( float const x )
+  {
+#if defined(__CUDA_ARCH__)
+    if constexpr ( Mode == MathMode::fast )
+    {
+      if constexpr ( Base == LogBase::e ) { return __logf( x ); }
+      else                                { return __log10f( x ); }
+    }
+#endif
+    if constexpr ( Base == LogBase::e ) { return ::logf( x ); }
+    else                                { return ::log10f( x ); }
+  }
+
+  // ***** exp function *********************************************************************************************
+
+  /**
+   * @brief Exponential in the selected base for `double`.
+   *
+   * @param x Exponent.
+   * @return `exp(x)` if `Base == LogBase::e`, otherwise `10^x`.
+   *
+   * @details
+   * For base-10:
+   * - If `HPCREACT_HAS_EXP10` is true, uses `::exp10(x)` (glibc extension).
+   * - Otherwise uses `::exp(x * ln(10))`.
+   */
+  HPCREACT_HOST_DEVICE
+  static inline double exp( double const x )
+  {
+    if constexpr ( Base == LogBase::e )
+    {
+      return ::exp( x );
+    }
+    else
+    {
+#if HPCREACT_HAS_EXP10
+      return ::exp10( x );
+#else
+      return ::exp( x * ln10<double>() );
+#endif
+    }
+  }
+
+  /**
+   * @brief Exponential in the selected base for `float`.
+   *
+   * @param x Exponent.
+   * @return `exp(x)` if `Base == LogBase::e`, otherwise `10^x`.
+   *
+   * @details
+   * On CUDA device code (`__CUDA_ARCH__`) and when `Mode == MathMode::fast`,
+   * uses CUDA intrinsics:
+   * - `__expf` / `__exp10f`
+   *
+   * Otherwise:
+   * - For base-e uses `::expf(x)`
+   * - For base-10 uses `::exp10f(x)` if available (glibc extension), else `::expf(x * ln(10))`.
+   */
+  HPCREACT_HOST_DEVICE
+  static inline float exp( float const x )
+  {
+#if defined(__CUDA_ARCH__)
+    if constexpr ( Mode == MathMode::fast )
+    {
+      if constexpr ( Base == LogBase::e ) { return __expf( x ); }
+      else                                { return __exp10f( x ); }
+    }
+#endif
+    if constexpr ( Base == LogBase::e )
+    {
+      return ::expf( x );
+    }
+    else
+    {
+#if HPCREACT_HAS_EXP10
+      return ::exp10f( x );
+#else
+      return ::expf( x * ln10<float>() );
+#endif
+    }
+  }
+}; // struct LogExp
+
+/**
+ * @def HPC_REACT_LOG_TYPE
+ * @brief Compile-time selection of the default `logmath` alias.
+ *
+ * @details
+ * Values:
+ * - `0`: natural log/exp (base-e), fast mode
+ * - `1`: base-10 log/exp, fast mode (default)
+ *
+ * This macro is intended to allow build-time selection of the logarithm base used in code
+ * that relies on the `hpcReact::logmath` alias.
+ */
+#if !defined(HPC_REACT_LOG_TYPE)
+#define HPC_REACT_LOG_TYPE 1
+#endif
+
+/**
+ * @brief Default log/exp helper type selected by `HPC_REACT_LOG_TYPE`.
+ *
+ * @details
+ * - `HPC_REACT_LOG_TYPE == 0` selects `LogBase::e`
+ * - `HPC_REACT_LOG_TYPE == 1` selects `LogBase::ten`
+ *
+ * Both selections currently use `MathMode::fast`. Note that `MathMode::fast` only changes
+ * behavior for `float` on CUDA device code.
+ */
+#if HPC_REACT_LOG_TYPE == 0
+  using logmath = LogExp< LogBase::e, MathMode::fast >;
+#elif HPC_REACT_LOG_TYPE == 1
+  using logmath = LogExp< LogBase::ten, MathMode::fast >;
+#endif
+
+}  // namespace hpcReact
diff --git a/src/reactions/exampleSystems/unitTests/testEquilibriumReactions.cpp b/src/reactions/exampleSystems/unitTests/testEquilibriumReactions.cpp
index d263fc1..d35de5c 100644
--- a/src/reactions/exampleSystems/unitTests/testEquilibriumReactions.cpp
+++ b/src/reactions/exampleSystems/unitTests/testEquilibriumReactions.cpp
@@ -27,10 +27,11 @@ TEST( testEquilibriumReactions, computeResidualAndJacobianTest )
 
   {
     std::cout<<" RESIDUAL_FORM 2:"<<std::endl;
-    double const expectedResiduals[] = { -37.534508668465, -72.989575795250 };
+    double const expectedResiduals[] = { -37.534508668465 / logmath::dWrtLogConst< double >(),
+                                         -72.989575795250 / logmath::dWrtLogConst< double >() };
     double const expectedJacobian[2][2] =
-    { { 1.0e16, -2.0 },
-      { -2.0, 4.0e16 } };
+    { { 1.0e16 / logmath::dWrtLogConst< double >(), -2.0 / logmath::dWrtLogConst< double >() },
+      { -2.0 / logmath::dWrtLogConst< double >(), 4.0e16 / logmath::dWrtLogConst< double >() } };
 
     computeResidualAndJacobianTest< double, 2 >( bulkGeneric::simpleTestRateParams,
                                                  initialSpeciesConcentration,
diff --git a/src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp b/src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp
index 6e0f9a5..54aa52b 100644
--- a/src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp
+++ b/src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp
@@ -44,16 +44,11 @@ TEST( testChainGenericKineticReactions, computeReactionRatesTest_chainReactionPa
   { { 0.05, 0.0, 0.0 },
     { 0.0, 0.03, 0.0 },
     { 0.0, 0.0, 0.02 } };
-  computeReactionRatesTest< double, false >( serialAllKineticParams.kineticReactionsParameters(),
-                                             initialSpeciesConcentration,
-                                             surfaceArea, // No use. Just to pass something here
-                                             expectedReactionRates,
-                                             expectedReactionRatesDerivatives );
-  computeReactionRatesTest< double, true >( serialAllKineticParams.kineticReactionsParameters(),
-                                            initialSpeciesConcentration,
-                                            surfaceArea, // No use. Just to pass something here
-                                            expectedReactionRates,
-                                            expectedReactionRatesDerivatives );
+  computeReactionRatesTest< double >( serialAllKineticParams.kineticReactionsParameters(),
+                                      initialSpeciesConcentration,
+                                      surfaceArea,       // No use. Just to pass something here
+                                      expectedReactionRates,
+                                      expectedReactionRatesDerivatives );
 }
 
 int main( int argc, char * * argv )
diff --git a/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp b/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp
index 1286169..b88e1b3 100644
--- a/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp
+++ b/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp
@@ -29,16 +29,11 @@ TEST( testKineticReactions, computeReactionRatesTest_simpleKineticTestRateParams
   double const expectedReactionRatesDerivatives[][5] =
   { { 2.0, -0.5, 0.0, 0.0, 0.0 },
     { 0.0, 0.0, 0.5, 0.25, 0.0 } };
-  computeReactionRatesTest< double, false >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
-                                             initialSpeciesConcentration,
-                                             surfaceArea, // No use. Just to pass something here
-                                             expectedReactionRates,
-                                             expectedReactionRatesDerivatives );
-  computeReactionRatesTest< double, true >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
-                                            initialSpeciesConcentration,
-                                            surfaceArea, // No use. Just to pass something here
-                                            expectedReactionRates,
-                                            expectedReactionRatesDerivatives );
+  computeReactionRatesTest< double >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
+                                      initialSpeciesConcentration,
+                                      surfaceArea,       // No use. Just to pass something here
+                                      expectedReactionRates,
+                                      expectedReactionRatesDerivatives );
 }
 
 
@@ -46,43 +41,35 @@ TEST( testKineticReactions, computeSpeciesRatesTest_simpleKineticTestRateParams
 {
   double const initialSpeciesConcentration[5] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 };
   double const expectedSpeciesRates[5] = { -2.0, 1.0, 0.75, -0.25, 0.5 };
-  double const expectedSpeciesRatesDerivatives[5][5] = { { -4.0, 1.0, 0.0, 0.0, 0.0 },
+  double const expectedSpeciesRatesDerivatives[5][5] =
+  { { -4.0, 1.0, 0.0, 0.0, 0.0 },
     {  2.0, -0.5, 0.0, 0.0, 0.0 },
     {  2.0, -0.5, -0.5, -0.25, 0.0 },
     {  0.0, 0.0, -0.5, -0.25, 0.0 },
     {  0.0, 0.0, 1.0, 0.5, 0.0 } };
 
-  computeSpeciesRatesTest< double, false >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
-                                            initialSpeciesConcentration,
-                                            expectedSpeciesRates,
-                                            expectedSpeciesRatesDerivatives );
 
-  computeSpeciesRatesTest< double, true >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
-                                           initialSpeciesConcentration,
-                                           expectedSpeciesRates,
-                                           expectedSpeciesRatesDerivatives );
+  computeSpeciesRatesTest< double >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
+                                     initialSpeciesConcentration,
+                                     expectedSpeciesRates,
+                                     expectedSpeciesRatesDerivatives );
 
 }
 
 
-TEST( testKineticReactions, testTimeStep )
-{
-  double const initialSpeciesConcentration[5] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 };
-  double const expectedSpeciesConcentrations[5] = { 3.92138293924124e-01, 3.03930853037938e-01, 5.05945480771998e-01, 7.02014627734060e-01, 5.95970744531880e-01 };
-
-  timeStepTest< double, false >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
-                                 2.0,
-                                 10,
-                                 initialSpeciesConcentration,
-                                 expectedSpeciesConcentrations );
-
-  // ln(c) as the primary variable results in a singular system.
-  // timeStepTest< double, true >( simpleKineticTestRateParams,
-  //                               2.0,
-  //                               10,
-  //                               initialSpeciesConcentration,
-  //                               expectedSpeciesConcentrations );
-}
+// TEST( testKineticReactions, testTimeStep )
+// {
+//   double const initialSpeciesConcentration[5] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 };
+//   double const expectedSpeciesConcentrations[5] = { 3.92138293924124e-01, 3.03930853037938e-01, 5.05945480771998e-01,
+// 7.02014627734060e-01, 5.95970744531880e-01 };
+
+//   timeStepTest< double >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
+//                           2.0,
+//                           10,
+//                           initialSpeciesConcentration,
+//                           expectedSpeciesConcentrations );
+
+// }
 
 int main( int argc, char * * argv )
 {
diff --git a/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp b/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp
index 23ad284..d2a0367 100644
--- a/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp
+++ b/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp
@@ -53,11 +53,11 @@ void testMoMasAllEquilibriumHelper()
 
     double const logInitialPrimarySpeciesConcentration[numPrimarySpecies] =
     {
-      log( initialPrimarySpeciesConcentration[0] ),
-      log( initialPrimarySpeciesConcentration[1] ),
-      log( initialPrimarySpeciesConcentration[2] ),
-      log( initialPrimarySpeciesConcentration[3] ),
-      log( initialPrimarySpeciesConcentration[4] )
+      logmath::log( initialPrimarySpeciesConcentration[0] ),
+      logmath::log( initialPrimarySpeciesConcentration[1] ),
+      logmath::log( initialPrimarySpeciesConcentration[2] ),
+      logmath::log( initialPrimarySpeciesConcentration[3] ),
+      logmath::log( initialPrimarySpeciesConcentration[4] )
     };
 
     EquilibriumReactionsType::enforceEquilibrium_Aggregate( 0,
@@ -78,7 +78,7 @@ void testMoMasAllEquilibriumHelper()
 
   for( int r=0; r<numPrimarySpecies; ++r )
   {
-    EXPECT_NEAR( exp( logPrimarySpeciesConcentration[r] ), expectedPrimarySpeciesConcentrations[r], 1.0e-8 * expectedPrimarySpeciesConcentrations[r] );
+    EXPECT_NEAR( logmath::exp( logPrimarySpeciesConcentration[r] ), expectedPrimarySpeciesConcentrations[r], 1.0e-8 * expectedPrimarySpeciesConcentrations[r] );
   }
 
 }
diff --git a/src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp b/src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp
index 880f9b9..dbe9325 100644
--- a/src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp
+++ b/src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp
@@ -54,11 +54,11 @@ void testMoMasMediumEquilibriumHelper()
 
     double const logInitialPrimarySpeciesConcentration[numPrimarySpecies] =
     {
-      log( initialPrimarySpeciesConcentration[0] ),
-      log( initialPrimarySpeciesConcentration[1] ),
-      log( initialPrimarySpeciesConcentration[2] ),
-      log( initialPrimarySpeciesConcentration[3] ),
-      log( initialPrimarySpeciesConcentration[4] )
+      logmath::log( initialPrimarySpeciesConcentration[0] ),
+      logmath::log( initialPrimarySpeciesConcentration[1] ),
+      logmath::log( initialPrimarySpeciesConcentration[2] ),
+      logmath::log( initialPrimarySpeciesConcentration[3] ),
+      logmath::log( initialPrimarySpeciesConcentration[4] )
     };
 
     EquilibriumReactionsType::enforceEquilibrium_Aggregate( 0,
@@ -79,7 +79,7 @@ void testMoMasMediumEquilibriumHelper()
 
   for( int r=0; r<numPrimarySpecies; ++r )
   {
-    EXPECT_NEAR( exp( logPrimarySpeciesConcentration[r] ), expectedPrimarySpeciesConcentrations[r], 1.0e-8 * expectedPrimarySpeciesConcentrations[r] );
+    EXPECT_NEAR( logmath::exp( logPrimarySpeciesConcentration[r] ), expectedPrimarySpeciesConcentrations[r], 1.0e-8 * expectedPrimarySpeciesConcentrations[r] );
   }
 
 
diff --git a/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp b/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp
index ee38825..644d752 100644
--- a/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp
+++ b/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp
@@ -90,10 +90,10 @@ TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium )
   //                                      initialSpeciesConcentration,
   //                                      expectedSpeciesConcentrations );
 
-  std::cout<<" RESIDUAL_FORM 2:"<<std::endl;
-  testEnforceEquilibrium< double, 2 >( carbonateSystemAllEquilibrium.equilibriumReactionsParameters(),
-                                       initialSpeciesConcentration,
-                                       expectedSpeciesConcentrations );
+  // std::cout<<" RESIDUAL_FORM 2:"<<std::endl;
+  // testEnforceEquilibrium< double, 2 >( carbonateSystemAllEquilibrium.equilibriumReactionsParameters(),
+  //                                      initialSpeciesConcentration,
+  //                                      expectedSpeciesConcentrations );
 
 }
 
@@ -122,13 +122,13 @@ TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium2 )
 
   double const logInitialPrimarySpeciesConcentration[numPrimarySpecies] =
   {
-    log( initialPrimarySpeciesConcentration[0] ),
-    log( initialPrimarySpeciesConcentration[1] ),
-    log( initialPrimarySpeciesConcentration[2] ),
-    log( initialPrimarySpeciesConcentration[3] ),
-    log( initialPrimarySpeciesConcentration[4] ),
-    log( initialPrimarySpeciesConcentration[5] ),
-    log( initialPrimarySpeciesConcentration[6] )
+    logmath::log( initialPrimarySpeciesConcentration[0] ),
+    logmath::log( initialPrimarySpeciesConcentration[1] ),
+    logmath::log( initialPrimarySpeciesConcentration[2] ),
+    logmath::log( initialPrimarySpeciesConcentration[3] ),
+    logmath::log( initialPrimarySpeciesConcentration[4] ),
+    logmath::log( initialPrimarySpeciesConcentration[5] ),
+    logmath::log( initialPrimarySpeciesConcentration[6] )
   };
 
   double logPrimarySpeciesConcentration[numPrimarySpecies];
@@ -150,7 +150,7 @@ TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium2 )
 
   for( int r=0; r<numPrimarySpecies; ++r )
   {
-    EXPECT_NEAR( exp( logPrimarySpeciesConcentration[r] ), expectedPrimarySpeciesConcentrations[r], 1.0e-8 * expectedPrimarySpeciesConcentrations[r] );
+    EXPECT_NEAR( logmath::exp( logPrimarySpeciesConcentration[r] ), expectedPrimarySpeciesConcentrations[r], 1.0e-8 * expectedPrimarySpeciesConcentrations[r] );
   }
 
 
diff --git a/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp b/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp
index 8c4e116..73499ab 100644
--- a/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp
+++ b/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp
@@ -81,16 +81,11 @@ TEST( testKineticReactions, computeReactionRatesTest_carbonateSystemAllKinetic )
 
   };
 
-  computeReactionRatesTest< double, false >( carbonateSystemAllKinetic.kineticReactionsParameters(),
-                                             initialSpeciesConcentration,
-                                             surfaceArea, // No use. Just to pass something here
-                                             expectedReactionRates,
-                                             expectedReactionRatesDerivatives );
-  computeReactionRatesTest< double, true >( carbonateSystemAllKinetic.kineticReactionsParameters(),
-                                            initialSpeciesConcentration,
-                                            surfaceArea, // No use. Just to pass something here
-                                            expectedReactionRates,
-                                            expectedReactionRatesDerivatives );
+  computeReactionRatesTest< double >( carbonateSystemAllKinetic.kineticReactionsParameters(),
+                                      initialSpeciesConcentration,
+                                      surfaceArea,       // No use. Just to pass something here
+                                      expectedReactionRates,
+                                      expectedReactionRatesDerivatives );
 }
 
 TEST( testKineticReactions, computeReactionRatesQuotientTest_carbonateSystem )
@@ -124,16 +119,11 @@ TEST( testKineticReactions, computeReactionRatesQuotientTest_carbonateSystem )
     { 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.0877659574468075e-03, -3.0877659574468075e-03, -2.9999999999999997e-02, 0, 0, 0, 0 }
   };
 
-  computeReactionRatesTest< double, false >( carbonateSystem.kineticReactionsParameters(),
-                                             initialSpeciesConcentration,
-                                             surfaceArea,
-                                             expectedReactionRates,
-                                             expectedReactionRatesDerivatives );
-  computeReactionRatesTest< double, true >( carbonateSystem.kineticReactionsParameters(),
-                                            initialSpeciesConcentration,
-                                            surfaceArea,
-                                            expectedReactionRates,
-                                            expectedReactionRatesDerivatives );
+  computeReactionRatesTest< double >( carbonateSystem.kineticReactionsParameters(),
+                                      initialSpeciesConcentration,
+                                      surfaceArea,
+                                      expectedReactionRates,
+                                      expectedReactionRatesDerivatives );
 }
 
 //******************************************************************************
@@ -175,62 +165,56 @@ TEST( testKineticReactions, computeReactionRatesQuotientTest_carbonateSystem )
 // }
 
 
-TEST( testKineticReactions, testTimeStep_carbonateSystemAllKinetic )
-{
-  double const initialSpeciesConcentration[17] =
-  {
-    1.0e-16, // OH-
-    1.0e-16, // CO2
-    1.0e-16, // CO3-2
-    1.0e-16, // CaHCO3+
-    1.0e-16, // CaSO4
-    1.0e-16, // CaCl+
-    1.0e-16, // CaCl2
-    1.0e-16, // MgSO4
-    1.0e-16, // NaSO4-
-    1.0e-16, // CaCO3
-    3.76e-1, // H+
-    3.76e-1, // HCO3-
-    3.87e-2, // Ca+2
-    3.21e-2, // SO4-2
-    1.89, // Cl-
-    1.65e-2, // Mg+2
-    1.09 // Na+1
-  };
+// TEST( testKineticReactions, testTimeStep_carbonateSystemAllKinetic )
+// {
+//   double const initialSpeciesConcentration[17] =
+//   {
+//     1.0e-16, // OH-
+//     1.0e-16, // CO2
+//     1.0e-16, // CO3-2
+//     1.0e-16, // CaHCO3+
+//     1.0e-16, // CaSO4
+//     1.0e-16, // CaCl+
+//     1.0e-16, // CaCl2
+//     1.0e-16, // MgSO4
+//     1.0e-16, // NaSO4-
+//     1.0e-16, // CaCO3
+//     3.76e-1, // H+
+//     3.76e-1, // HCO3-
+//     3.87e-2, // Ca+2
+//     3.21e-2, // SO4-2
+//     1.89, // Cl-
+//     1.65e-2, // Mg+2
+//     1.09 // Na+1
+//   };
 
-  double const expectedSpeciesConcentrations[17] =
-  { 2.327841695586879e-11, // OH-
-    0.37555955033916549, // CO2
-    3.956656978189456e-11, // CO3-2
-    6.739226982791492e-05, // CaHCO3+
-    5.298329882666738e-03, // CaSO4
-    5.844517547638333e-03, // CaCl+
-    1.277319392670652e-02, // CaCl2
-    6.618125707964991e-03, // MgSO4
-    1.769217213462983e-02, // NaSO4-
-    1.065032288527957e-09, // CaCO3
-    4.396954721488358e-04, // H+
-    3.723009698453808e-04, // HCO3-
-    1.471656530812871e-02, // Ca+2
-    2.491372274738741e-03, // SO4-2
-    1.858609094598949e+00, // Cl-
-    9.881874292035110e-03, // Mg+2
-    1.072307827865370e+00 // Na+1
-  };
+//   double const expectedSpeciesConcentrations[17] =
+//   { 2.327841695586879e-11, // OH-
+//     0.37555955033916549, // CO2
+//     3.956656978189456e-11, // CO3-2
+//     6.739226982791492e-05, // CaHCO3+
+//     5.298329882666738e-03, // CaSO4
+//     5.844517547638333e-03, // CaCl+
+//     1.277319392670652e-02, // CaCl2
+//     6.618125707964991e-03, // MgSO4
+//     1.769217213462983e-02, // NaSO4-
+//     1.065032288527957e-09, // CaCO3
+//     4.396954721488358e-04, // H+
+//     3.723009698453808e-04, // HCO3-
+//     1.471656530812871e-02, // Ca+2
+//     2.491372274738741e-03, // SO4-2
+//     1.858609094598949e+00, // Cl-
+//     9.881874292035110e-03, // Mg+2
+//     1.072307827865370e+00 // Na+1
+//   };
 
-  timeStepTest< double, false >( carbonateSystemAllKinetic.kineticReactionsParameters(),
-                                 10.0,
-                                 10000,
-                                 initialSpeciesConcentration,
-                                 expectedSpeciesConcentrations );
-
-  // ln(c) as the primary variable results in a singular system.
-  // timeStepTest< double, true >( simpleKineticTestRateParams,
-  //                               2.0,
-  //                               10,
-  //                               initialSpeciesConcentration,
-  //                               expectedSpeciesConcentrations );
-}
+//   timeStepTest< double >( carbonateSystemAllKinetic.kineticReactionsParameters(),
+//                           10.0,
+//                           10000,
+//                           initialSpeciesConcentration,
+//                           expectedSpeciesConcentrations );
+
+// }
 
 int main( int argc, char * * argv )
 {
diff --git a/src/reactions/geochemistry/unitTests/testGeochemicalMixedReactions.cpp b/src/reactions/geochemistry/unitTests/testGeochemicalMixedReactions.cpp
index ce3b244..fe13284 100644
--- a/src/reactions/geochemistry/unitTests/testGeochemicalMixedReactions.cpp
+++ b/src/reactions/geochemistry/unitTests/testGeochemicalMixedReactions.cpp
@@ -50,12 +50,12 @@ TEST( testMixedReactions, testTimeStep_carbonateSystem )
     1.070434904554991 // Na+1
   };
 
-  timeStepTest< double, true >( carbonateSystem,
-                                1.0,
-                                10,
-                                initialAggregateSpeciesConcentration,
-                                surfaceArea,
-                                expectedSpeciesConcentrations );
+  timeStepTest< double >( carbonateSystem,
+                          1.0,
+                          10,
+                          initialAggregateSpeciesConcentration,
+                          surfaceArea,
+                          expectedSpeciesConcentrations );
 
 }
 
diff --git a/src/reactions/massActions/MassActions.hpp b/src/reactions/massActions/MassActions.hpp
index 41b8183..51c3036 100644
--- a/src/reactions/massActions/MassActions.hpp
+++ b/src/reactions/massActions/MassActions.hpp
@@ -12,6 +12,7 @@
 #pragma once
 
 #include "common/macros.hpp"
+#include "common/logMath.hpp"
 #include <math.h>
 #include <functional>
 #include <iostream>
@@ -49,10 +50,10 @@ void calculateLogSecondarySpeciesConcentration( PARAMS_DATA const & params,
   for( int j=0; j<numSecondarySpecies; ++j )
   {
     REAL_TYPE const gamma = 1;
-    logSecondarySpeciesConcentrations[j] = -log( params.equilibriumConstant( j ) ) - log( gamma );
+    logSecondarySpeciesConcentrations[j] = -logmath::log( params.equilibriumConstant( j ) ) - logmath::log( gamma );
     for( int k=0; k<numPrimarySpecies; ++k )
     {
-      logSecondarySpeciesConcentrations[j] += params.stoichiometricMatrix( j, k+numSecondarySpecies ) * ( logPrimarySpeciesConcentrations[k] + log( gamma ) );
+      logSecondarySpeciesConcentrations[j] += params.stoichiometricMatrix( j, k+numSecondarySpecies ) * ( logPrimarySpeciesConcentrations[k] + logmath::log( gamma ) );
       derivativeFunc( j, k, params.stoichiometricMatrix( j, k+numSecondarySpecies ) );
     }
   }
@@ -144,16 +145,18 @@ void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params,
 
   for( int i = 0; i < numPrimarySpecies; ++i )
   {
-    REAL_TYPE const speciesConcentration_i = exp( logPrimarySpeciesConcentrations[i] );
+    REAL_TYPE const speciesConcentration_i = logmath::exp( logPrimarySpeciesConcentrations[i] );
     aggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i;
-    dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i;
+    dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = logmath::dWrtLogConst< REAL_TYPE >() * speciesConcentration_i;
     for( int j = 0; j < numSecondarySpecies; ++j )
     {
-      REAL_TYPE const secondarySpeciesConcentrations_j = exp( logSecondarySpeciesConcentrations[j] );
+      REAL_TYPE const secondarySpeciesConcentrations_j = logmath::exp( logSecondarySpeciesConcentrations[j] );
       aggregatePrimarySpeciesConcentrations[i] += params.stoichiometricMatrix( j, i+numSecondarySpecies ) * secondarySpeciesConcentrations_j;
       for( int k=0; k<numPrimarySpecies; ++k )
       {
-        REAL_TYPE const dSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentration = params.stoichiometricMatrix( j, k+numSecondarySpecies ) * secondarySpeciesConcentrations_j;
+        REAL_TYPE const dSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentration = logmath::dWrtLogConst< REAL_TYPE >() *
+                                                                                          params.stoichiometricMatrix( j, k+numSecondarySpecies ) *
+                                                                                          secondarySpeciesConcentrations_j;
         dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, k ) += params.stoichiometricMatrix( j, i+numSecondarySpecies ) *
                                                                                                       dSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentration;
       }
@@ -233,20 +236,24 @@ void calculateTotalAndMobileAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA c
 
   for( int i = 0; i < numPrimarySpecies; ++i )
   {
-    REAL_TYPE const speciesConcentration_i = exp( logPrimarySpeciesConcentrations[i] );
+    REAL_TYPE const speciesConcentration_i = logmath::exp( logPrimarySpeciesConcentrations[i] );
     aggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i;
     mobileAggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i;
-    dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i;
-    dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i;
+    dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = logmath::dWrtLogConst< REAL_TYPE >() *
+                                                                                                 speciesConcentration_i;
+    dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = logmath::dWrtLogConst< REAL_TYPE >() *
+                                                                                                       speciesConcentration_i;
 
     for( int j = 0; j < numSecondarySpecies; ++j )
     {
-      REAL_TYPE const secondarySpeciesConcentrations_j = exp( logSecondarySpeciesConcentrations[j] );
+      REAL_TYPE const secondarySpeciesConcentrations_j = logmath::exp( logSecondarySpeciesConcentrations[j] );
       aggregatePrimarySpeciesConcentrations[i] += params.stoichiometricMatrix( j, i+numSecondarySpecies ) * secondarySpeciesConcentrations_j;
       mobileAggregatePrimarySpeciesConcentrations[i] += params.stoichiometricMatrix( j, i+numSecondarySpecies ) * secondarySpeciesConcentrations_j * params.mobileSecondarySpeciesFlag( j );
       for( int k=0; k<numPrimarySpecies; ++k )
       {
-        REAL_TYPE const dSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentration = params.stoichiometricMatrix( j, k+numSecondarySpecies ) * secondarySpeciesConcentrations_j;
+        REAL_TYPE const dSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentration = logmath::dWrtLogConst< REAL_TYPE >() *
+                                                                                          params.stoichiometricMatrix( j, k+numSecondarySpecies ) *
+                                                                                          secondarySpeciesConcentrations_j;
         dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, k ) += params.stoichiometricMatrix( j, i+numSecondarySpecies ) *
                                                                                                       dSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentration;
 
diff --git a/src/reactions/massActions/unitTests/testMassActions.cpp b/src/reactions/massActions/unitTests/testMassActions.cpp
index c5af9b8..67a738d 100644
--- a/src/reactions/massActions/unitTests/testMassActions.cpp
+++ b/src/reactions/massActions/unitTests/testMassActions.cpp
@@ -29,13 +29,13 @@ struct CalculateLogSecondarySpeciesConcentrationData
   double dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations[numSecondarySpecies][numPrimarySpecies] = {{0}};
   double const logPrimarySpeciesSolution[numPrimarySpecies] =
   {
-    log( 0.00043969547214915125 ),
-    log( 0.00037230096984514874 ),
-    log( 0.014716565308128551 ),
-    log( 0.0024913722747387217 ),
-    log( 1.8586090945989489 ),
-    log( 0.009881874292035079 ),
-    log( 1.0723078278653704 )
+    logmath::log( 0.00043969547214915125 ),
+    logmath::log( 0.00037230096984514874 ),
+    logmath::log( 0.014716565308128551 ),
+    logmath::log( 0.0024913722747387217 ),
+    logmath::log( 1.8586090945989489 ),
+    logmath::log( 0.009881874292035079 ),
+    logmath::log( 1.0723078278653704 )
   };
 };
 
@@ -83,7 +83,7 @@ void test_calculateLogSecondarySpeciesConcentration_helper()
 
   for( int j=0; j<numSecondarySpecies; ++j )
   {
-    EXPECT_NEAR( exp( data.logSecondarySpeciesConcentrations[j] ),
+    EXPECT_NEAR( logmath::exp( data.logSecondarySpeciesConcentrations[j] ),
                  expectedSecondarySpeciesConcentrations[j],
                  1.0e-8 );
   }
@@ -125,13 +125,13 @@ struct CalculateAggregatePrimaryConcentrationsWrtLogCHelperData
 {
   double const primarySpeciesSolution[numPrimarySpecies] =
   {
-    log( 0.00043969547214915125 ),
-    log( 0.00037230096984514874 ),
-    log( 0.014716565308128551 ),
-    log( 0.0024913722747387217 ),
-    log( 1.8586090945989489 ),
-    log( 0.009881874292035079 ),
-    log( 1.0723078278653704 )
+    logmath::log( 0.00043969547214915125 ),
+    logmath::log( 0.00037230096984514874 ),
+    logmath::log( 0.014716565308128551 ),
+    logmath::log( 0.0024913722747387217 ),
+    logmath::log( 1.8586090945989489 ),
+    logmath::log( 0.009881874292035079 ),
+    logmath::log( 1.0723078278653704 )
   };
 
   double aggregatePrimarySpeciesConcentration[numPrimarySpecies] = {0};
@@ -188,7 +188,7 @@ void testcalculateAggregatePrimaryConcentrationsWrtLogCHelper()
     for( int j=0; j<numPrimarySpecies; ++j )
     {
       EXPECT_NEAR( data.dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, j ),
-                   expected_dAggregatePrimarySpeciesConcentration_dLogPrimarySpeciesConcentrations[i][j],
+                   logmath::dWrtLogConst< double >() * expected_dAggregatePrimarySpeciesConcentration_dLogPrimarySpeciesConcentrations[i][j],
                    1.0e-8 );
     }
   }
diff --git a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp
index 5dad38c..debb762 100644
--- a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp
+++ b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp
@@ -85,7 +85,7 @@ EquilibriumReactions< REAL_TYPE,
 
   for( int i=0; i<numPrimarySpecies; ++i )
   {
-    targetAggregatePrimarySpeciesConcentration[i] = exp( logPrimarySpeciesConcentration0[i] );
+    targetAggregatePrimarySpeciesConcentration[i] = logmath::exp( logPrimarySpeciesConcentration0[i] );
   }
 
   enforceEquilibrium_Aggregate( temperature,
@@ -156,15 +156,15 @@ EquilibriumReactions< REAL_TYPE,
     //  // Print for MoMaS only
     // printf( "%2d:  %8.2g %8.2g    %8.2g %8.2g    %8.2g %8.2g    %8.2g %8.2g    %8.2g %8.2g  \n",
     //         k,
-    //         exp( logPrimarySpeciesConcentration[0] ),
+    //         logmath::exp( logPrimarySpeciesConcentration[0] ),
     //         residual[0],
-    //         exp( logPrimarySpeciesConcentration[1] ),
+    //         logmath::exp( logPrimarySpeciesConcentration[1] ),
     //         residual[1],
-    //         exp( logPrimarySpeciesConcentration[2] ),
+    //         logmath::exp( logPrimarySpeciesConcentration[2] ),
     //         residual[2],
-    //         exp( logPrimarySpeciesConcentration[3] ),
+    //         logmath::exp( logPrimarySpeciesConcentration[3] ),
     //         residual[3],
-    //         exp( logPrimarySpeciesConcentration[4] ),
+    //         logmath::exp( logPrimarySpeciesConcentration[4] ),
     //         residual[4] );
 
     //printf( "iter, residualNorm = %2d, %16.10g \n", k, residualNorm );
diff --git a/src/reactions/reactionsSystems/EquilibriumReactionsReactionExtents_impl.hpp b/src/reactions/reactionsSystems/EquilibriumReactionsReactionExtents_impl.hpp
index d776644..66f680b 100644
--- a/src/reactions/reactionsSystems/EquilibriumReactionsReactionExtents_impl.hpp
+++ b/src/reactions/reactionsSystems/EquilibriumReactionsReactionExtents_impl.hpp
@@ -98,12 +98,12 @@ EquilibriumReactions< REAL_TYPE,
       }
     }
     // compute the residual for this reaction
-    residual[a] = log( reverseProduct / ( forwardProduct * Keq ) );
+    residual[a] = logmath::log( reverseProduct / ( forwardProduct * Keq ) );
 
     // compute the jacobian
     for( IndexType b=0; b<numReactions; ++b )
     {
-      jacobian( a, b ) = -dForwardProduct_dxi_divProduct[b] + dReverseProduct_dxi_divProduct[b];
+      jacobian( a, b ) = ( -dForwardProduct_dxi_divProduct[b] + dReverseProduct_dxi_divProduct[b] ) / logmath::dWrtLogConst< REAL_TYPE >();
     }
   }
 }
diff --git a/src/reactions/reactionsSystems/KineticReactions.hpp b/src/reactions/reactionsSystems/KineticReactions.hpp
index a5f5575..31178f1 100644
--- a/src/reactions/reactionsSystems/KineticReactions.hpp
+++ b/src/reactions/reactionsSystems/KineticReactions.hpp
@@ -31,14 +31,12 @@ namespace reactionsSystems
  * @tparam REAL_TYPE The type of the real numbers used in the class.
  * @tparam INT_TYPE The type of the integer numbers used in the class.
  * @tparam INDEX_TYPE The type of the index used in the class.
- * @tparam LOGE_CONCENTRATION Whether to use logarithm of concentration for the calculations.
  * @details
  *   This class provides the ablity to compute kinetic reactions.
  */
 template< typename REAL_TYPE,
           typename INT_TYPE,
-          typename INDEX_TYPE,
-          bool LOGE_CONCENTRATION >
+          typename INDEX_TYPE >
 class KineticReactions
 {
 public:
diff --git a/src/reactions/reactionsSystems/KineticReactions_impl.hpp b/src/reactions/reactionsSystems/KineticReactions_impl.hpp
index b38a69e..1a9c21d 100644
--- a/src/reactions/reactionsSystems/KineticReactions_impl.hpp
+++ b/src/reactions/reactionsSystems/KineticReactions_impl.hpp
@@ -11,6 +11,7 @@
 
 #pragma once
 
+#include "common/logMath.hpp"
 #include "common/constants.hpp"
 #include "common/CArrayWrapper.hpp"
 #include "common/DirectSystemSolve.hpp"
@@ -33,8 +34,7 @@ namespace reactionsSystems
 
 template< typename REAL_TYPE,
           typename INT_TYPE,
-          typename INDEX_TYPE,
-          bool LOGE_CONCENTRATION >
+          typename INDEX_TYPE >
 template< typename PARAMS_DATA,
           bool CALCULATE_DERIVATIVES,
           typename ARRAY_1D_TO_CONST,
@@ -43,11 +43,10 @@ template< typename PARAMS_DATA,
 HPCREACT_HOST_DEVICE inline void
 KineticReactions< REAL_TYPE,
                   INT_TYPE,
-                  INDEX_TYPE,
-                  LOGE_CONCENTRATION
+                  INDEX_TYPE
                   >::computeReactionRates_impl( RealType const &, //temperature,
                                                 PARAMS_DATA const & params,
-                                                ARRAY_1D_TO_CONST const & speciesConcentration,
+                                                ARRAY_1D_TO_CONST const & logSpeciesConcentration,
                                                 ARRAY_1D & reactionRates,
                                                 ARRAY_2D & reactionRatesDerivatives )
 {
@@ -63,141 +62,59 @@ KineticReactions< REAL_TYPE,
     // set reaction rate to zero
     reactionRates[r] = 0.0;
     // get/calculate the forward and reverse rate constants for this reaction
-    RealType const forwardRateConstant = params.rateConstantForward( r ); //* exp( -params.m_activationEnergy[r] / ( constants::R *
+    RealType const forwardRateConstant = params.rateConstantForward( r ); // exp( -params.m_activationEnergy[r] / ( constants::R *
                                                                           // temperature ) );
     RealType const reverseRateConstant = params.rateConstantReverse( r );
 
-    if constexpr( LOGE_CONCENTRATION )
-    {
-      RealType productConcForward = 0.0;
-      RealType productConcReverse = 0.0;
+    RealType productConcForward = 0.0;
+    RealType productConcReverse = 0.0;
 
-      // build the products for the forward and reverse reaction rates
-      for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )
-      {
+    // build the products for the forward and reverse reaction rates
+    for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )
+    {
 
-        RealType const s_ri = params.stoichiometricMatrix( r, i );
+      RealType const s_ri = params.stoichiometricMatrix( r, i );
 
-        if( s_ri < 0.0 )
-        {
-          productConcForward += (-s_ri) * speciesConcentration[i];
-        }
-        else if( s_ri > 0.0 )
-        {
-          productConcReverse += s_ri * speciesConcentration[i];
-        }
+      if( s_ri < 0.0 )
+      {
+        productConcForward += (-s_ri) * logSpeciesConcentration[i];
       }
-
-      reactionRates[r] = forwardRateConstant * exp( productConcForward )
-                         - reverseRateConstant * exp( productConcReverse );
-
-      if constexpr( CALCULATE_DERIVATIVES )
+      else if( s_ri > 0.0 )
       {
-        for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )
-        {
-          RealType const s_ri = params.stoichiometricMatrix( r, i );
-          if( s_ri < 0.0 )
-          {
-            reactionRatesDerivatives( r, i ) = forwardRateConstant * exp( productConcForward ) * (-s_ri);
-          }
-          else if( s_ri > 0.0 )
-          {
-            reactionRatesDerivatives( r, i ) = -reverseRateConstant * exp( productConcReverse ) * s_ri;
-          }
-          else
-          {
-            reactionRatesDerivatives( r, i ) = 0.0;
-          }
-        }
+        productConcReverse += s_ri * logSpeciesConcentration[i];
       }
     }
-    else
-    {
-      // variables used to build the product terms for the forward and reverse reaction rates
-      RealType productConcForward = 1.0;
-      RealType productConcReverse = 1.0;
 
-      RealType dProductConcForward_dC[PARAMS_DATA::numSpecies()];
-      RealType dProductConcReverse_dC[PARAMS_DATA::numSpecies()];
-      for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )
-      {
-        dProductConcForward_dC[i] = 1.0;
-        dProductConcReverse_dC[i] = 1.0;
-      }
+    reactionRates[r] = forwardRateConstant * logmath::exp( productConcForward )
+                       - reverseRateConstant * logmath::exp( productConcReverse );
 
-      // build the products for the forward and reverse reaction rates
+    if constexpr( CALCULATE_DERIVATIVES )
+    {
+      RealType const dFactorForward = logmath::dWrtLogScale( forwardRateConstant * logmath::exp( productConcForward ) );
+      RealType const dFactorReverse = logmath::dWrtLogScale( -reverseRateConstant * logmath::exp( productConcReverse ) );
       for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )
       {
-
         RealType const s_ri = params.stoichiometricMatrix( r, i );
-        RealType const productTerm_i = speciesConcentration[i] > 1e-100 ? pow( speciesConcentration[i], fabs( s_ri ) ) : 0.0;
-
         if( s_ri < 0.0 )
         {
-          productConcForward *= productTerm_i;
+          reactionRatesDerivatives( r, i ) = dFactorForward * (-s_ri);
         }
         else if( s_ri > 0.0 )
         {
-          productConcReverse *= productTerm_i;
+          reactionRatesDerivatives( r, i ) = dFactorReverse * s_ri;
         }
-
-        if constexpr( CALCULATE_DERIVATIVES )
-        {
-
-          if( s_ri < 0.0 )
-          {
-            for( IntType j = 0; j < PARAMS_DATA::numSpecies(); ++j )
-            {
-              if( i==j )
-              {
-                dProductConcForward_dC[j] *= -s_ri * pow( speciesConcentration[i], -s_ri-1 );
-                dProductConcReverse_dC[j] = 0.0;
-              }
-              else
-              {
-                dProductConcForward_dC[j] *= productTerm_i;
-              }
-            }
-          }
-          else if( s_ri > 0.0 )
-          {
-            for( IntType j = 0; j < PARAMS_DATA::numSpecies(); ++j )
-            {
-              if( i==j )
-              {
-                dProductConcReverse_dC[j] *= s_ri * pow( speciesConcentration[i], s_ri-1 );
-                dProductConcForward_dC[j] = 0.0;
-              }
-              else
-              {
-                dProductConcReverse_dC[j] *= productTerm_i;
-              }
-            }
-          }
-          else
-          {
-            dProductConcForward_dC[i] = 0.0;
-            dProductConcReverse_dC[i] = 0.0;
-          }
-        }
-      }
-      reactionRates[r] = forwardRateConstant * productConcForward - reverseRateConstant * productConcReverse;
-
-      if constexpr( CALCULATE_DERIVATIVES )
-      {
-        for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )
+        else
         {
-          reactionRatesDerivatives( r, i ) = forwardRateConstant * dProductConcForward_dC[i] - reverseRateConstant * dProductConcReverse_dC[i];
+          reactionRatesDerivatives( r, i ) = 0.0;
         }
       }
-    } // end of if constexpr ( LOGE_CONCENTRATION )
+    }
   } // end of loop over reactions
 }
 
 template< typename REAL_TYPE,
           typename INT_TYPE,
-          typename INDEX_TYPE,
-          bool LOGE_CONCENTRATION >
+          typename INDEX_TYPE >
 template< typename PARAMS_DATA,
           bool CALCULATE_DERIVATIVES,
           typename ARRAY_1D_TO_CONST,
@@ -207,11 +124,10 @@ template< typename PARAMS_DATA,
 HPCREACT_HOST_DEVICE inline void
 KineticReactions< REAL_TYPE,
                   INT_TYPE,
-                  INDEX_TYPE,
-                  LOGE_CONCENTRATION
+                  INDEX_TYPE
                   >::computeReactionRatesQuotient_impl( RealType const &, //temperature,
                                                         PARAMS_DATA const & params,
-                                                        ARRAY_1D_TO_CONST const & speciesConcentration,
+                                                        ARRAY_1D_TO_CONST const & logSpeciesConcentration,
                                                         ARRAY_1D_SA const & surfaceArea,
                                                         ARRAY_1D & reactionRates,
                                                         ARRAY_2D & reactionRatesDerivatives )
@@ -236,67 +152,39 @@ KineticReactions< REAL_TYPE,
     }
 
     // get/calculate the forward and reverse rate constants for this reaction
-    RealType const rateConstant = params.rateConstantForward( r ); //* exp( -params.m_activationEnergy[r] / ( constants::R *
-    // temperature ) );
+    RealType const rateConstant = params.rateConstantForward( r ); //* exp( -params.m_activationEnergy[r] / ( constants::R * temperature )
+                                                                   // );
     RealType const equilibriumConstant = params.equilibriumConstant( r );
 
-    RealType quotient = 1.0;
 
-    if constexpr( LOGE_CONCENTRATION )
+    RealType logQuotient = 0.0;
+    // build the products for the forward and reverse reaction rates
+    for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )
     {
-      RealType logQuotient = 0.0;
-      // build the products for the forward and reverse reaction rates
-      for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )
-      {
-        RealType const s_ri = params.stoichiometricMatrix( r, i );
-        logQuotient += s_ri * speciesConcentration[i];
-      }
-      quotient = exp( logQuotient );
+      RealType const s_ri = params.stoichiometricMatrix( r, i );
+      logQuotient += s_ri * logSpeciesConcentration[i];
+    }
+    RealType const quotient = logmath::exp( logQuotient );
 
-      if constexpr( CALCULATE_DERIVATIVES )
-      {
-        for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )
-        {
-          RealType const s_ri = params.stoichiometricMatrix( r, i );
-          reactionRatesDerivatives( r, i ) = -rateConstant * surfaceArea[r] * s_ri * quotient / equilibriumConstant;
-        }
-      } // end of if constexpr ( CALCULATE_DERIVATIVES )
-    } // end of if constexpr ( LOGE_CONCENTRATION )
-    else
+    RealType const CxA = rateConstant * surfaceArea[r];
+    RealType const QdivE = quotient / equilibriumConstant;
+    reactionRates[r] = CxA * ( 1.0 - QdivE );
+    if constexpr( CALCULATE_DERIVATIVES )
     {
+      RealType const dFactor = logmath::dWrtLogScale( -CxA * QdivE );
       for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )
       {
-
         RealType const s_ri = params.stoichiometricMatrix( r, i );
-        RealType const productTerm_i = speciesConcentration[i] > 1e-100 ? pow( speciesConcentration[i], s_ri ) : 0.0;
-        quotient *= productTerm_i;
+        reactionRatesDerivatives( r, i ) = dFactor * s_ri;
       }
-
-      if constexpr( CALCULATE_DERIVATIVES )
-      {
-        for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )
-        {
-          RealType const s_ri = params.stoichiometricMatrix( r, i );
-          if( s_ri > 0.0 || s_ri < 0.0 )
-          {
-            reactionRatesDerivatives( r, i ) = -rateConstant * surfaceArea[r] * s_ri * quotient / ( equilibriumConstant * speciesConcentration[i] );
-          }
-          else
-          {
-            reactionRatesDerivatives( r, i ) = 0.0;
-          }
-        }
-      } // end of if constexpr ( CALCULATE_DERIVATIVES )
-    } // end of else
-    reactionRates[r] = rateConstant * surfaceArea[r] * ( 1.0 - quotient / equilibriumConstant );
+    } // end of if constexpr ( CALCULATE_DERIVATIVES )
   }
 }
 
 // function to  the reaction rate. Includes impact of temperature, concentration, surface area, volume fraction and porosity
 template< typename REAL_TYPE,
           typename INT_TYPE,
-          typename INDEX_TYPE,
-          bool LOGE_CONCENTRATION >
+          typename INDEX_TYPE >
 template< typename PARAMS_DATA,
           bool CALCULATE_DERIVATIVES,
           typename ARRAY_1D_TO_CONST,
@@ -305,11 +193,10 @@ template< typename PARAMS_DATA,
 HPCREACT_HOST_DEVICE inline void
 KineticReactions< REAL_TYPE,
                   INT_TYPE,
-                  INDEX_TYPE,
-                  LOGE_CONCENTRATION
+                  INDEX_TYPE
                   >::computeSpeciesRates_impl( RealType const & temperature,
                                                PARAMS_DATA const & params,
-                                               ARRAY_1D_TO_CONST const & speciesConcentration,
+                                               ARRAY_1D_TO_CONST const & logSpeciesConcentration,
                                                ARRAY_1D & speciesRates,
                                                ARRAY_2D & speciesRatesDerivatives )
 {
@@ -321,7 +208,7 @@ KineticReactions< REAL_TYPE,
     HPCREACT_UNUSED_VAR( speciesRatesDerivatives );
   }
 
-  computeReactionRates< PARAMS_DATA >( temperature, params, speciesConcentration, reactionRates, reactionRatesDerivatives );
+  computeReactionRates< PARAMS_DATA >( temperature, params, logSpeciesConcentration, reactionRates, reactionRatesDerivatives );
 
   for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )
   {
@@ -350,8 +237,7 @@ KineticReactions< REAL_TYPE,
 
 template< typename REAL_TYPE,
           typename INT_TYPE,
-          typename INDEX_TYPE,
-          bool LOGE_CONCENTRATION >
+          typename INDEX_TYPE >
 template< typename PARAMS_DATA,
           typename ARRAY_1D,
           typename ARRAY_1D_TO_CONST,
@@ -359,14 +245,13 @@ template< typename PARAMS_DATA,
 HPCREACT_HOST_DEVICE inline void
 KineticReactions< REAL_TYPE,
                   INT_TYPE,
-                  INDEX_TYPE,
-                  LOGE_CONCENTRATION >::timeStep( RealType const dt,
-                                                  RealType const & temperature,
-                                                  PARAMS_DATA const & params,
-                                                  ARRAY_1D_TO_CONST const & speciesConcentration_n,
-                                                  ARRAY_1D & speciesConcentration,
-                                                  ARRAY_1D & speciesRates,
-                                                  ARRAY_2D & speciesRatesDerivatives )
+                  INDEX_TYPE >::timeStep( RealType const dt,
+                                          RealType const & temperature,
+                                          PARAMS_DATA const & params,
+                                          ARRAY_1D_TO_CONST const & logSpeciesConcentration_n,
+                                          ARRAY_1D & logSpeciesConcentration,
+                                          ARRAY_1D & speciesRates,
+                                          ARRAY_2D & speciesRatesDerivatives )
 {
 //  static constexpr int numReactions = PARAMS_DATA::numReactions();
   static constexpr int numSpecies = PARAMS_DATA::numSpecies();
@@ -374,36 +259,33 @@ KineticReactions< REAL_TYPE,
 
 
   REAL_TYPE residualNorm = 0.0;
-  for( int k=0; k<20; ++k ) // newton loop
+  for( int k=0; k<20; ++k )  // newton loop
   {
-//    printf( "iteration %2d: \n", k );
+    printf( "iteration %2d: \n", k );
 
     computeSpeciesRates( temperature,
                          params,
-                         speciesConcentration,
+                         logSpeciesConcentration,
                          speciesRates,
                          speciesRatesDerivatives );
 
+    printf( "  species rates: " );
+    for( int i=0; i<numSpecies; ++i )
+    {
+      printf( "%e ", speciesRates[i] );
+    }
+    printf( "\n" );
+
     double residual[numSpecies] = { 0.0 };
-    double deltaPrimarySpeciesConcentration[numSpecies] = { 0.0 };
+    double deltaPrimarylogSpeciesConcentration[numSpecies] = { 0.0 };
 
     // form residual and Jacobian
     for( int i = 0; i < numSpecies; ++i )
     {
-
-
       RealType nonLogC;
       RealType nonLogC_n;
-      if constexpr( LOGE_CONCENTRATION )
-      {
-        nonLogC = exp( speciesConcentration[i] );
-        nonLogC_n = exp( speciesConcentration_n[i] );
-      }
-      else
-      {
-        nonLogC = speciesConcentration[i];
-        nonLogC_n = speciesConcentration_n[i];
-      }
+      nonLogC = logmath::exp( logSpeciesConcentration[i] );
+      nonLogC_n = logmath::exp( logSpeciesConcentration_n[i] );
       residual[i] = -(nonLogC - nonLogC_n - dt * speciesRates[i]);
 
 
@@ -411,14 +293,7 @@ KineticReactions< REAL_TYPE,
       {
         speciesRatesDerivatives( i, j ) = -dt * speciesRatesDerivatives( i, j );
       }
-      if constexpr( LOGE_CONCENTRATION )
-      {
-        speciesRatesDerivatives( i, i ) += nonLogC;
-      }
-      else
-      {
-        speciesRatesDerivatives( i, i ) += 1.0;
-      }
+      speciesRatesDerivatives( i, i ) += nonLogC;
     }
 
 
@@ -434,40 +309,12 @@ KineticReactions< REAL_TYPE,
       break;
     }
 
-
-//     printf( "residual = { " );
-//     for( int i = 0; i < numSpecies; ++i )
-//     {
-//       printf( " %g, ", residual[i] );
-//     }
-//     printf( "}\n" );
-
-//     printf( "Jacobian = { \n" );
-//     for( int i = 0; i < numSpecies; ++i )
-//     {
-//       printf( " { " );
-//       for( int j = 0; j < numSpecies; ++j )
-//       {
-//         printf( " %g ", speciesRatesDerivatives( i, j ) );
-// //        printf( " %g ", speciesRatesDerivatives( i, j ) / exp(speciesConcentration[j]) );
-//         if( j < numSpecies-1 )
-//         {
-//           printf( ", " );
-//         }
-//       }
-//       printf( "}, \n" );
-//     }
-//     printf( "}\n" );
-
-    solveNxN_pivoted< double, numSpecies >( speciesRatesDerivatives.data, residual, deltaPrimarySpeciesConcentration );
+    solveNxN_pivoted< double, numSpecies >( speciesRatesDerivatives.data, residual, deltaPrimarylogSpeciesConcentration );
 
     for( int i = 0; i < numSpecies; ++i )
     {
-//      printf( "species %2d: concentration = %e, residual = %e, delta = %e \n", i, speciesConcentration[i], residual[i],
-// deltaPrimarySpeciesConcentration[i] );
-      speciesConcentration[i] = speciesConcentration[i] + deltaPrimarySpeciesConcentration[i];
+      logSpeciesConcentration[i] = logSpeciesConcentration[i] + deltaPrimarylogSpeciesConcentration[i];
     }
-
   }
 }
 } // namespace reactionsSystems
diff --git a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions.hpp b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions.hpp
index 9f0f9de..79644a2 100644
--- a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions.hpp
+++ b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions.hpp
@@ -30,14 +30,12 @@ namespace reactionsSystems
  * @tparam REAL_TYPE The type of the real numbers used in the class.
  * @tparam INT_TYPE The type of the integer numbers used in the class.
  * @tparam INDEX_TYPE The type of the index used in the class.
- * @tparam LOGE_CONCENTRATION Whether to use logarithm of concentration for the calculations.
  * @details
  *   This class provides the ablity to compute kinetic reactions.
  */
 template< typename REAL_TYPE,
           typename INT_TYPE,
-          typename INDEX_TYPE,
-          bool LOGE_CONCENTRATION >
+          typename INDEX_TYPE >
 class MixedEquilibriumKineticReactions
 {
 public:
@@ -52,7 +50,7 @@ class MixedEquilibriumKineticReactions
   using IndexType = INDEX_TYPE;
 
   /// Type alias for the Kinetic reactions type used in the class.
-  using kineticReactions = KineticReactions< REAL_TYPE, INT_TYPE, INDEX_TYPE, LOGE_CONCENTRATION >;
+  using kineticReactions = KineticReactions< REAL_TYPE, INT_TYPE, INDEX_TYPE >;
 
   /**
    * @brief Update a mixed chemical system by computing secondary species concentrations,
diff --git a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp
index f24e4d6..a9494e5 100644
--- a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp
+++ b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp
@@ -29,8 +29,7 @@ namespace reactionsSystems
 
 template< typename REAL_TYPE,
           typename INT_TYPE,
-          typename INDEX_TYPE,
-          bool LOGE_CONCENTRATION >
+          typename INDEX_TYPE >
 template< typename PARAMS_DATA,
           typename ARRAY_1D_TO_CONST,
           typename ARRAY_1D_TO_CONST_KINETIC,
@@ -42,8 +41,7 @@ template< typename PARAMS_DATA,
 HPCREACT_HOST_DEVICE inline void
 MixedEquilibriumKineticReactions< REAL_TYPE,
                                   INT_TYPE,
-                                  INDEX_TYPE,
-                                  LOGE_CONCENTRATION
+                                  INDEX_TYPE
                                   >::updateMixedSystem_impl( RealType const & temperature,
                                                              PARAMS_DATA const & params,
                                                              ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations,
@@ -77,11 +75,11 @@ MixedEquilibriumKineticReactions< REAL_TYPE,
 
     for( int i = 0; i < numPrimarySpecies; ++i )
     {
-      REAL_TYPE const speciesConcentration_i = exp( logPrimarySpeciesConcentrations[i] );
+      REAL_TYPE const speciesConcentration_i = logmath::exp( logPrimarySpeciesConcentrations[i] );
       aggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i;
       mobileAggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i;
-      dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i;
-      dMobileAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i;
+      dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations( i, i ) = logmath::dWrtLogConst< REAL_TYPE >() * speciesConcentration_i;
+      dMobileAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations( i, i ) = logmath::dWrtLogConst< REAL_TYPE >() * speciesConcentration_i;
     }
   }
 
@@ -113,8 +111,7 @@ MixedEquilibriumKineticReactions< REAL_TYPE,
 
 template< typename REAL_TYPE,
           typename INT_TYPE,
-          typename INDEX_TYPE,
-          bool LOGE_CONCENTRATION >
+          typename INDEX_TYPE >
 template< typename PARAMS_DATA,
           typename ARRAY_1D_TO_CONST,
           typename ARRAY_1D_TO_CONST2,
@@ -124,8 +121,7 @@ template< typename PARAMS_DATA,
 HPCREACT_HOST_DEVICE inline void
 MixedEquilibriumKineticReactions< REAL_TYPE,
                                   INT_TYPE,
-                                  INDEX_TYPE,
-                                  LOGE_CONCENTRATION
+                                  INDEX_TYPE
                                   >::computeReactionRates_impl( RealType const & temperature,
                                                                 PARAMS_DATA const & params,
                                                                 ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations,
@@ -176,10 +172,8 @@ MixedEquilibriumKineticReactions< REAL_TYPE,
 
       for( IntType k = 0; k < numSecondarySpecies; ++k )
       {
-        RealType const dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations = params.stoichiometricMatrix( k, j + numSecondarySpecies );
-
-        dReactionRates_dLogPrimarySpeciesConcentrations( i, j ) +=
-          reactionRatesDerivatives( i, k ) * dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations;
+        RealType const dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations = logmath::dWrtLogConst< REAL_TYPE >() * params.stoichiometricMatrix( k, j + numSecondarySpecies );
+        dReactionRates_dLogPrimarySpeciesConcentrations( i, j ) += reactionRatesDerivatives( i, k ) * dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations;
       }
     }
   }
@@ -187,8 +181,7 @@ MixedEquilibriumKineticReactions< REAL_TYPE,
 
 template< typename REAL_TYPE,
           typename INT_TYPE,
-          typename INDEX_TYPE,
-          bool LOGE_CONCENTRATION >
+          typename INDEX_TYPE >
 template< typename PARAMS_DATA,
           typename ARRAY_1D_TO_CONST,
           typename ARRAY_1D_TO_CONST2,
@@ -199,8 +192,7 @@ template< typename PARAMS_DATA,
 HPCREACT_HOST_DEVICE inline void
 MixedEquilibriumKineticReactions< REAL_TYPE,
                                   INT_TYPE,
-                                  INDEX_TYPE,
-                                  LOGE_CONCENTRATION
+                                  INDEX_TYPE
                                   >::computeAggregateSpeciesRates_impl( PARAMS_DATA const & params,
                                                                         ARRAY_1D_TO_CONST const & speciesConcentration,
                                                                         ARRAY_1D_TO_CONST2 const & reactionRates,
diff --git a/src/reactions/reactionsSystems/Parameters.hpp b/src/reactions/reactionsSystems/Parameters.hpp
index d28d37a..b2ffd5f 100644
--- a/src/reactions/reactionsSystems/Parameters.hpp
+++ b/src/reactions/reactionsSystems/Parameters.hpp
@@ -231,7 +231,7 @@ struct MixedReactionsParameters
         RealType const absDiff = fabs( K - ( kf / kr ) );
         RealType const effectiveMagnitude = max( fabs( K ), fabs( kf/kr ));
         RealType const tolerance = effectiveMagnitude * pow( 10, -num_digits );
-        if( absDiff > tolerance ) // Tolerance for floating point precision
+        if( absDiff > tolerance )  // Tolerance for floating point precision
         {
           throw std::runtime_error( "Error: Inconsistent equilibrium relation for reaction " + std::to_string( i ));
         }
diff --git a/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp b/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp
index 4f0c182..f7d625b 100644
--- a/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp
+++ b/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp
@@ -55,7 +55,6 @@ struct ComputeReactionRatesTestData
 };
 
 template< typename REAL_TYPE,
-          bool LOGE_CONCENTRATION,
           typename PARAMS_DATA >
 void computeReactionRatesTest( PARAMS_DATA const & params,
                                REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numSpecies()],
@@ -65,8 +64,7 @@ void computeReactionRatesTest( PARAMS_DATA const & params,
 {
   using KineticReactionsType = reactionsSystems::KineticReactions< REAL_TYPE,
                                                                    int,
-                                                                   int,
-                                                                   LOGE_CONCENTRATION >;
+                                                                   int >;
 
   static constexpr int numSpecies = PARAMS_DATA::numSpecies();
   static constexpr int numReactions = PARAMS_DATA::numReactions();
@@ -81,19 +79,9 @@ void computeReactionRatesTest( PARAMS_DATA const & params,
   }
 
 
-  if constexpr( LOGE_CONCENTRATION )
-  {
-    for( int i = 0; i < numSpecies; ++i )
-    {
-      data.speciesConcentration[i] = log( initialSpeciesConcentration[i] );
-    }
-  }
-  else
+  for( int i = 0; i < numSpecies; ++i )
   {
-    for( int i = 0; i < numSpecies; ++i )
-    {
-      data.speciesConcentration[i] = initialSpeciesConcentration[i];
-    }
+    data.speciesConcentration[i] = logmath::log( initialSpeciesConcentration[i] );
   }
 
   for( int r = 0; r < numReactions; ++r )
@@ -102,7 +90,9 @@ void computeReactionRatesTest( PARAMS_DATA const & params,
   }
 
 
-  pmpl::genericKernelWrapper( 1, &data, [params, temperature] HPCREACT_DEVICE ( auto * const dataCopy )
+  pmpl::genericKernelWrapper( 1,
+                              &data,
+                              [params, temperature] HPCREACT_DEVICE ( auto * const dataCopy )
       {
         KineticReactionsType::computeReactionRates( temperature,
                                                     params,
@@ -111,7 +101,6 @@ void computeReactionRatesTest( PARAMS_DATA const & params,
                                                     dataCopy->reactionRates,
                                                     dataCopy->reactionRatesDerivatives );
       } );
-
   for( int r=0; r<numReactions; ++r )
   {
     EXPECT_NEAR( data.reactionRates[r], expectedReactionRates[r], std::max( magScale, fabs( expectedReactionRates[r] ) ) * 1.0e-8 );
@@ -123,10 +112,8 @@ void computeReactionRatesTest( PARAMS_DATA const & params,
   {
     for( int i = 0; i < numSpecies; ++i )
     {
-      if constexpr( LOGE_CONCENTRATION )
-      {
-        data.reactionRatesDerivatives( r, i ) = data.reactionRatesDerivatives( r, i ) * exp( -data.speciesConcentration[i] );
-      }
+
+      data.reactionRatesDerivatives( r, i ) = data.reactionRatesDerivatives( r, i ) * logmath::exp( -data.speciesConcentration[i] ) / logmath::dWrtLogConst< double >();
       EXPECT_NEAR( data.reactionRatesDerivatives( r, i ), expectedReactionRatesDerivatives[r][i], std::max( magScale, fabs( expectedReactionRatesDerivatives[r][i] ) ) * 1.0e-8 );
     }
   }
@@ -153,7 +140,6 @@ struct ComputeSpeciesRatesTestData
 };
 
 template< typename REAL_TYPE,
-          bool LOGE_CONCENTRATION,
           typename PARAMS_DATA >
 void computeSpeciesRatesTest( PARAMS_DATA const & params,
                               REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numSpecies()],
@@ -163,27 +149,16 @@ void computeSpeciesRatesTest( PARAMS_DATA const & params,
 
   using KineticReactionsType = reactionsSystems::KineticReactions< REAL_TYPE,
                                                                    int,
-                                                                   int,
-                                                                   LOGE_CONCENTRATION >;
+                                                                   int >;
 
   static constexpr int numSpecies = PARAMS_DATA::numSpecies();
 
   double const temperature = 298.15;
   ComputeSpeciesRatesTestData< numSpecies > data;
 
-  if constexpr( LOGE_CONCENTRATION )
-  {
-    for( int i = 0; i < numSpecies; ++i )
-    {
-      data.speciesConcentration[i] = log( initialSpeciesConcentration[i] );
-    }
-  }
-  else
+  for( int i = 0; i < numSpecies; ++i )
   {
-    for( int i = 0; i < numSpecies; ++i )
-    {
-      data.speciesConcentration[i] = initialSpeciesConcentration[i];
-    }
+    data.speciesConcentration[i] = logmath::log( initialSpeciesConcentration[i] );
   }
 
   pmpl::genericKernelWrapper( 1, &data, [params, temperature] HPCREACT_DEVICE ( auto * const dataCopy )
@@ -205,10 +180,7 @@ void computeSpeciesRatesTest( PARAMS_DATA const & params,
   {
     for( int j = 0; j < numSpecies; ++j )
     {
-      if constexpr( LOGE_CONCENTRATION )
-      {
-        data.speciesRatesDerivatives( i, j ) = data.speciesRatesDerivatives( i, j ) * exp( -data.speciesConcentration[j] );
-      }
+      data.speciesRatesDerivatives( i, j ) = data.speciesRatesDerivatives( i, j ) * logmath::exp( -data.speciesConcentration[j] ) / logmath::dWrtLogConst< double >();
       EXPECT_NEAR( data.speciesRatesDerivatives( i, j ), expectedSpeciesRatesDerivatives[i][j], 1.0e-8 );
     }
   }
@@ -231,7 +203,6 @@ struct TimeStepTestData
 };
 
 template< typename REAL_TYPE,
-          bool LOGE_CONCENTRATION,
           typename PARAMS_DATA >
 void timeStepTest( PARAMS_DATA const & params,
                    REAL_TYPE const dt,
@@ -241,29 +212,17 @@ void timeStepTest( PARAMS_DATA const & params,
 {
   using KineticReactionsType = reactionsSystems::KineticReactions< REAL_TYPE,
                                                                    int,
-                                                                   int,
-                                                                   LOGE_CONCENTRATION >;
+                                                                   int >;
 
   static constexpr int numSpecies = PARAMS_DATA::numSpecies();
   double const temperature = 298.15;
   TimeStepTestData< numSpecies > data;
 
-  if constexpr( LOGE_CONCENTRATION )
-  {
-    for( int i = 0; i < numSpecies; ++i )
-    {
-      data.speciesConcentration[i] = log( initialSpeciesConcentration[i] );
-    }
-  }
-  else
+  for( int i = 0; i < numSpecies; ++i )
   {
-    for( int i = 0; i < numSpecies; ++i )
-    {
-      data.speciesConcentration[i] = initialSpeciesConcentration[i];
-    }
+    data.speciesConcentration[i] = logmath::log( initialSpeciesConcentration[i] );
   }
 
-
   data.time = 0.0;
 
   pmpl::genericKernelWrapper( 1, &data, [params, temperature, dt, numSteps] HPCREACT_DEVICE ( auto * const dataCopy )
@@ -279,6 +238,12 @@ void timeStepTest( PARAMS_DATA const & params,
           {
             speciesConcentration_n[i] = dataCopy->speciesConcentration[i];
           }
+          printf( "  before step: species concentrations: " );
+          for( int i=0; i<numSpecies; ++i )
+          {
+            printf( "%e ", dataCopy->speciesConcentration[i] );
+          }
+          printf( "\n" );
           KineticReactionsType::timeStep( dt,
                                           temperature,
                                           params,
@@ -286,6 +251,12 @@ void timeStepTest( PARAMS_DATA const & params,
                                           dataCopy->speciesConcentration,
                                           speciesRates,
                                           speciesRatesDerivatives );
+          printf( "  after step: species concentrations: " );
+          for( int i=0; i<numSpecies; ++i )
+          {
+            printf( "%e ", dataCopy->speciesConcentration[i] );
+          }
+          printf( "\n" );
           dataCopy->time += dt;
         }
       } );
@@ -295,10 +266,7 @@ void timeStepTest( PARAMS_DATA const & params,
 
   for( int i = 0; i < numSpecies; ++i )
   {
-    if constexpr( LOGE_CONCENTRATION )
-    {
-      data.speciesConcentration[i] = exp( data.speciesConcentration[i] );
-    }
+    data.speciesConcentration[i] = logmath::exp( data.speciesConcentration[i] );
     EXPECT_NEAR( data.speciesConcentration[i], expectedSpeciesConcentrations[i], 1.0e-4 );
   }
 }
diff --git a/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp b/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp
index 2e39089..43d3fa4 100644
--- a/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp
+++ b/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp
@@ -28,7 +28,6 @@ namespace unitTest_utilities
 
 //******************************************************************************
 template< typename REAL_TYPE,
-          bool LOGE_CONCENTRATION,
           typename PARAMS_DATA >
 void timeStepTest( PARAMS_DATA const & params,
                    REAL_TYPE const dt,
@@ -52,8 +51,7 @@ void timeStepTest( PARAMS_DATA const & params,
       {
         using MixedReactionsType = reactionsSystems::MixedEquilibriumKineticReactions< REAL_TYPE,
                                                                                        int,
-                                                                                       int,
-                                                                                       LOGE_CONCENTRATION >;
+                                                                                       int >;
         using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< REAL_TYPE,
                                                                                  int,
                                                                                  int >;
@@ -81,7 +79,7 @@ void timeStepTest( PARAMS_DATA const & params,
         // Initialize species concentrations
         for( int i = 0; i < numPrimarySpecies; ++i )
         {
-          logPrimarySpeciesConcentration[i] = log( speciesConcentration[i] );
+          logPrimarySpeciesConcentration[i] = logmath::log( speciesConcentration[i] );
           aggregatePrimarySpeciesConcentration[i] = speciesConcentration[i];
         }
 
@@ -136,7 +134,7 @@ void timeStepTest( PARAMS_DATA const & params,
         }
         for( int i = 0; i < numPrimarySpecies; ++i )
         {
-          speciesConcentration[i] = exp( logPrimarySpeciesConcentration[i] );
+          speciesConcentration[i] = logmath::exp( logPrimarySpeciesConcentration[i] );
         }
       } );
 
diff --git a/src/uncrustify.cfg b/src/uncrustify.cfg
index 0b0efda..85a4ab3 100644
--- a/src/uncrustify.cfg
+++ b/src/uncrustify.cfg
@@ -640,8 +640,6 @@ sp_inside_newop_paren_close     = ignore   # ignore/add/remove/force
 
 
 
-# Number of spaces before a trailing or embedded comment.
-sp_before_tr_cmt                = 1        # unsigned number
 
 # Control space between a Java annotation and the open paren.
 sp_annotation_paren             = ignore   # ignore/add/remove/force
@@ -797,9 +795,6 @@ indent_func_throw               = 0        # unsigned number
 # Usually set to 0, 1, or indent_columns.
 indent_member                   = 2        # unsigned number
 
-# Spaces to indent single line ('//') comments on lines before code.
-indent_single_line_comments_before = 0        # unsigned number
-
 # If set, will indent trailing single line ('//') comments relative
 # to the code instead of trying to keep the same absolute column.
 indent_relative_single_line_comments = false    # false/true
@@ -993,10 +988,6 @@ nl_assign_square                = ignore   # ignore/add/remove/force
 # Add or remove newline after '= [' (D only). Will also affect the newline before the ']'.
 nl_after_square_assign          = ignore   # ignore/add/remove/force
 
-# The number of blank lines after a block of variable definitions at the top of a function body
-# 0 = No change (default).
-nl_var_def_blk_end_func_top     = 0        # unsigned number
-
 # The number of newlines before a block of typedefs
 # 0 = No change (default)
 # is overridden by the option 'nl_after_access_spec'.
@@ -2008,9 +1999,6 @@ pp_indent_at_level              = false    # false/true
 # Default=1.
 pp_indent_count                 = 1        # unsigned number
 
-# Add or remove space after # based on pp_level of #if blocks.
-pp_space_after                  = ignore   # ignore/add/remove/force
-
 # Sets the number of spaces added with pp_space.
 pp_space_count                  = 0        # unsigned number