diff --git a/src/CabanaPD_Constants.hpp b/src/CabanaPD_Constants.hpp
index 6e9d11d..4ed7404 100644
--- a/src/CabanaPD_Constants.hpp
+++ b/src/CabanaPD_Constants.hpp
@@ -12,13 +12,12 @@
 #ifndef CONSTANTS_H
 #define CONSTANTS_H
 
-#include "mpi.h"
+#include "Kokkos_Core.hpp"
 
 namespace CabanaPD
 {
 
-// FIXME: Kokkos::numbers::pi_v<double> when minimum 4.0 is required.
-constexpr double pi = M_PI;
+constexpr double pi = Kokkos::numbers::pi_v<double>;
 
 } // namespace CabanaPD
 #endif
diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp
index ea7d87b..0e87bd5 100644
--- a/src/CabanaPD_Force.hpp
+++ b/src/CabanaPD_Force.hpp
@@ -60,7 +60,7 @@
 #ifndef FORCE_H
 #define FORCE_H
 
-#include <cmath>
+#include <Kokkos_Core.hpp>
 
 #include <CabanaPD_ForceModels.hpp>
 #include <CabanaPD_Particles.hpp>
@@ -86,8 +86,8 @@ getDistanceComponents( const PosType& x, const PosType& u, const int i,
     rx = xi_x + eta_u;
     ry = xi_y + eta_v;
     rz = xi_z + eta_w;
-    r = sqrt( rx * rx + ry * ry + rz * rz );
-    xi = sqrt( xi_x * xi_x + xi_y * xi_y + xi_z * xi_z );
+    r = Kokkos::sqrt( rx * rx + ry * ry + rz * rz );
+    xi = Kokkos::sqrt( xi_x * xi_x + xi_y * xi_y + xi_z * xi_z );
     s = ( r - xi ) / xi;
 }
 
@@ -112,7 +112,7 @@ KOKKOS_INLINE_FUNCTION void getLinearizedDistanceComponents(
     const double eta_v = u( j, 1 ) - u( i, 1 );
     xi_z = x( j, 2 ) - x( i, 2 );
     const double eta_w = u( j, 2 ) - u( i, 2 );
-    xi = sqrt( xi_x * xi_x + xi_y * xi_y + xi_z * xi_z );
+    xi = Kokkos::sqrt( xi_x * xi_x + xi_y * xi_y + xi_z * xi_z );
     s = ( xi_x * eta_u + xi_y * eta_v + xi_z * eta_w ) / ( xi * xi );
 }
 
diff --git a/src/CabanaPD_Prenotch.hpp b/src/CabanaPD_Prenotch.hpp
index b147513..fb700e3 100644
--- a/src/CabanaPD_Prenotch.hpp
+++ b/src/CabanaPD_Prenotch.hpp
@@ -20,35 +20,6 @@
 
 namespace CabanaPD
 {
-KOKKOS_INLINE_FUNCTION
-double abs( const double a )
-{
-#if KOKKOS_VERSION >= 30700
-    return Kokkos::abs( a );
-#else
-    return Kokkos::Experimental::abs( a );
-#endif
-}
-
-KOKKOS_INLINE_FUNCTION
-double fmin( const double a, const double b )
-{
-#if KOKKOS_VERSION >= 30700
-    return Kokkos::fmin( a, b );
-#else
-    return Kokkos::Experimental::fmin( a, b );
-#endif
-}
-
-KOKKOS_INLINE_FUNCTION
-double fmax( const double a, const double b )
-{
-#if KOKKOS_VERSION >= 30700
-    return Kokkos::fmax( a, b );
-#else
-    return Kokkos::Experimental::fmax( a, b );
-#endif
-}
 
 KOKKOS_INLINE_FUNCTION
 double dot( Kokkos::Array<double, 3> a, Kokkos::Array<double, 3> b )
@@ -60,7 +31,10 @@ double dot( Kokkos::Array<double, 3> a, Kokkos::Array<double, 3> b )
 }
 
 KOKKOS_INLINE_FUNCTION
-double norm( Kokkos::Array<double, 3> a ) { return sqrt( dot( a, a ) ); }
+double norm( Kokkos::Array<double, 3> a )
+{
+    return Kokkos::sqrt( dot( a, a ) );
+}
 
 KOKKOS_INLINE_FUNCTION
 Kokkos::Array<double, 3> cross( Kokkos::Array<double, 3> a,
@@ -108,9 +82,9 @@ int linePlaneIntersection( const Kokkos::Array<double, 3> p0,
                            const Kokkos::Array<double, 3> l,
                            const double tol = 1e-10 )
 {
-    if ( abs( dot( l, n ) ) < tol )
+    if ( Kokkos::abs( dot( l, n ) ) < tol )
     {
-        if ( abs( dot( diff( p0, l0 ), n ) ) < tol )
+        if ( Kokkos::abs( dot( diff( p0, l0 ), n ) ) < tol )
             return 1;
         else
             return 2;
@@ -131,7 +105,7 @@ int bondPrenotchIntersection( const Kokkos::Array<double, 3> v1,
     double norm_cross_v1_v2 = norm( cross_v1_v2 );
 
     // Check if v1 and v2 are parallel.
-    assert( abs( norm( cross( v1, v2 ) ) ) > tol );
+    assert( Kokkos::abs( norm( cross( v1, v2 ) ) ) > tol );
 
     // Define plane normal.
     auto n = scale( cross_v1_v2, 1.0 / norm_cross_v1_v2 );
@@ -159,8 +133,10 @@ int bondPrenotchIntersection( const Kokkos::Array<double, 3> v1,
         double lj2 = -dot( cross( diff( x_j, p0 ), v1 ), cross_v1_v2 ) /
                      norm2_cross_v1_v2;
 
-        if ( !( fmin( li1, lj1 ) > 1 + tol || fmax( li1, lj1 ) < -tol ||
-                fmin( li2, lj2 ) > 1 + tol || fmax( li2, lj2 ) < -tol ) )
+        if ( !( Kokkos::fmin( li1, lj1 ) > 1 + tol ||
+                Kokkos::fmax( li1, lj1 ) < -tol ||
+                Kokkos::fmin( li2, lj2 ) > 1 + tol ||
+                Kokkos::fmax( li2, lj2 ) < -tol ) )
             keep_bond = 0;
     }
     // Case II: no intersection.
@@ -168,7 +144,7 @@ int bondPrenotchIntersection( const Kokkos::Array<double, 3> v1,
     // Case III: single point intersection.
     else if ( case_flag == 3 )
     {
-        assert( abs( dot( l, n ) ) > tol );
+        assert( Kokkos::abs( dot( l, n ) ) > tol );
 
         // Check if intersection point belongs to the bond.
         auto d = dot( diff( p0, l0 ), n ) / dot( l, n );
@@ -185,7 +161,8 @@ int bondPrenotchIntersection( const Kokkos::Array<double, 3> v1,
             double l2 = -dot( cross( diff( p, p0 ), v1 ), cross_v1_v2 ) /
                         norm2_cross_v1_v2;
 
-            if ( -tol < fmin( l1, l2 ) && fmax( l1, l2 ) < 1 + tol )
+            if ( -tol < Kokkos::fmin( l1, l2 ) &&
+                 Kokkos::fmax( l1, l2 ) < 1 + tol )
                 // Intersection.
                 keep_bond = 0;
         }
diff --git a/src/force/CabanaPD_ForceModels_LPS.hpp b/src/force/CabanaPD_ForceModels_LPS.hpp
index 8f2999f..7cd947c 100644
--- a/src/force/CabanaPD_ForceModels_LPS.hpp
+++ b/src/force/CabanaPD_ForceModels_LPS.hpp
@@ -12,6 +12,8 @@
 #ifndef FORCE_MODELS_LPS_H
 #define FORCE_MODELS_LPS_H
 
+#include <Kokkos_Core.hpp>
+
 #include <CabanaPD_ForceModels.hpp>
 #include <CabanaPD_Types.hpp>
 
@@ -95,11 +97,11 @@ struct ForceModel<LPS, Fracture> : public ForceModel<LPS, Elastic>
         G0 = _G0;
         if ( influence_type == 1 )
         {
-            s0 = sqrt( 5.0 * G0 / 9.0 / K / delta ); // 1/xi
+            s0 = Kokkos::sqrt( 5.0 * G0 / 9.0 / K / delta ); // 1/xi
         }
         else
         {
-            s0 = sqrt( 8.0 * G0 / 15.0 / K / delta ); // 1
+            s0 = Kokkos::sqrt( 8.0 * G0 / 15.0 / K / delta ); // 1
         }
         bond_break_coeff = ( 1.0 + s0 ) * ( 1.0 + s0 );
     }
diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp
index f0200ec..b488685 100644
--- a/src/force/CabanaPD_ForceModels_PMB.hpp
+++ b/src/force/CabanaPD_ForceModels_PMB.hpp
@@ -12,6 +12,8 @@
 #ifndef FORCE_MODELS_PMB_H
 #define FORCE_MODELS_PMB_H
 
+#include <Kokkos_Core.hpp>
+
 #include <CabanaPD_Constants.hpp>
 #include <CabanaPD_ForceModels.hpp>
 #include <CabanaPD_Types.hpp>
@@ -88,7 +90,7 @@ struct ForceModel<PMB, Fracture, TemperatureIndependent>
     {
         base_type::set_param( _delta, _K );
         G0 = _G0;
-        s0 = sqrt( 5.0 * G0 / 9.0 / K / delta );
+        s0 = Kokkos::sqrt( 5.0 * G0 / 9.0 / K / delta );
         bond_break_coeff = ( 1.0 + s0 ) * ( 1.0 + s0 );
     }
 
diff --git a/src/force/CabanaPD_Force_HertzianContact.hpp b/src/force/CabanaPD_Force_HertzianContact.hpp
index 7c2d918..688636b 100644
--- a/src/force/CabanaPD_Force_HertzianContact.hpp
+++ b/src/force/CabanaPD_Force_HertzianContact.hpp
@@ -12,7 +12,7 @@
 #ifndef CONTACT_HERTZIAN_H
 #define CONTACT_HERTZIAN_H
 
-#include <cmath>
+#include <Kokkos_Core.hpp>
 
 #include <CabanaPD_Force.hpp>
 #include <CabanaPD_Input.hpp>
@@ -62,8 +62,8 @@ class Force<MemorySpace, HertzianModel>
         const int n_frozen = particles.frozenOffset();
         const int n_local = particles.localOffset();
 
-        const double coeff_h_n = 4.0 / 3.0 * Es * std::sqrt( Rs );
-        const double coeff_h_d = -2.0 * sqrt( 5.0 / 6.0 ) * beta;
+        const double coeff_h_n = 4.0 / 3.0 * Es * Kokkos::sqrt( Rs );
+        const double coeff_h_d = -2.0 * Kokkos::sqrt( 5.0 / 6.0 ) * beta;
 
         const auto vol = particles.sliceVolume();
         const auto rho = particles.sliceDensity();
@@ -76,11 +76,6 @@ class Force<MemorySpace, HertzianModel>
 
         auto contact_full = KOKKOS_LAMBDA( const int i, const int j )
         {
-            using Kokkos::abs;
-            using Kokkos::min;
-            using Kokkos::pow;
-            using Kokkos::sqrt;
-
             double fcx_i = 0.0;
             double fcy_i = 0.0;
             double fcz_i = 0.0;
@@ -97,9 +92,10 @@ class Force<MemorySpace, HertzianModel>
             double Sn = 0.0;
             if ( delta_n < 0.0 )
             {
-                coeff =
-                    min( 0.0, -coeff_h_n * pow( abs( delta_n ), 3.0 / 2.0 ) );
-                Sn = 2.0 * Es * sqrt( Rs * abs( delta_n ) );
+                coeff = Kokkos::min(
+                    0.0, -coeff_h_n *
+                             Kokkos::pow( Kokkos::abs( delta_n ), 3.0 / 2.0 ) );
+                Sn = 2.0 * Es * Kokkos::sqrt( Rs * Kokkos::abs( delta_n ) );
             }
 
             coeff /= vol( i );
@@ -118,7 +114,7 @@ class Force<MemorySpace, HertzianModel>
                                                  vy, vz, vn );
 
             double ms = ( rho( i ) * vol( i ) ) / 2.0;
-            double fnd = coeff_h_d * sqrt( Sn * ms ) * vn / vol( i );
+            double fnd = coeff_h_d * Kokkos::sqrt( Sn * ms ) * vn / vol( i );
 
             fcx_i = fnd * rx / r;
             fcy_i = fnd * ry / r;
diff --git a/src/force/CabanaPD_HertzianContact.hpp b/src/force/CabanaPD_HertzianContact.hpp
index 56435da..cd5a8d3 100644
--- a/src/force/CabanaPD_HertzianContact.hpp
+++ b/src/force/CabanaPD_HertzianContact.hpp
@@ -12,7 +12,7 @@
 #ifndef CONTACTMODEL_HERTZIAN_H
 #define CONTACTMODEL_HERTZIAN_H
 
-#include <cmath>
+#include <Kokkos_Core.hpp>
 
 #include <CabanaPD_Force.hpp>
 #include <CabanaPD_Input.hpp>
@@ -46,17 +46,14 @@ struct HertzianModel : public ContactModel
     void set_param( const double _radius, const double _nu, const double _E,
                     const double _e )
     {
-        using Kokkos::log;
-        using Kokkos::pow;
-        using Kokkos::sqrt;
-
         nu = _nu;
         radius = _radius;
         Rs = 0.5 * radius;
-        Es = _E / ( 2.0 * pow( 1.0 - nu, 2.0 ) );
+        Es = _E / ( 2.0 * Kokkos::pow( 1.0 - nu, 2.0 ) );
         e = _e;
-        double ln_e = log( e );
-        beta = -ln_e / sqrt( pow( ln_e, 2.0 ) + pow( pi, 2.0 ) );
+        double ln_e = Kokkos::log( e );
+        beta = -ln_e / Kokkos::sqrt( Kokkos::pow( ln_e, 2.0 ) +
+                                     Kokkos::pow( pi, 2.0 ) );
     }
 };