diff --git a/libsymmetrix/CMakeLists.txt b/libsymmetrix/CMakeLists.txt
index c3fa384..3e97ce0 100644
--- a/libsymmetrix/CMakeLists.txt
+++ b/libsymmetrix/CMakeLists.txt
@@ -11,6 +11,7 @@ set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 
 find_package(BLAS REQUIRED)
 add_subdirectory(external/sphericart)
+option(SYMMETRIX_KOKKOS "Enable Kokkos mode in Symmetrix." ON)
 if(SYMMETRIX_KOKKOS)
     if(TARGET Kokkos::kokkos)
         message(STATUS "Symmetrix: Found existing Kokkos target.")
diff --git a/libsymmetrix/source/mace_kokkos.cpp b/libsymmetrix/source/mace_kokkos.cpp
index a0ee77f..80ba9bd 100644
--- a/libsymmetrix/source/mace_kokkos.cpp
+++ b/libsymmetrix/source/mace_kokkos.cpp
@@ -24,27 +24,30 @@ using Kokkos::TeamVectorRange;
 using Kokkos::TeamVectorMDRange;
 using Kokkos::View;
 
-MACEKokkos::MACEKokkos(std::string filename)
+template <typename Precision>
+MACEKokkos<Precision>::MACEKokkos(std::string filename)
 {
     load_from_json(filename);
 }
 
-MACEKokkos::~MACEKokkos()
+template <typename Precision>
+MACEKokkos<Precision>::~MACEKokkos()
 {
     Kokkos::fence();
 
     M0_monomials = Kokkos::View<Kokkos::View<int**,Kokkos::LayoutRight>*,Kokkos::SharedSpace>();
-    M0_weights = Kokkos::View<Kokkos::View<double***,Kokkos::LayoutRight>*,Kokkos::SharedSpace>();
+    M0_weights = Kokkos::View<Kokkos::View<Precision***,Kokkos::LayoutRight>*,Kokkos::SharedSpace>();
     M0_poly_spec = Kokkos::View<Kokkos::View<int**,Kokkos::LayoutRight>*,Kokkos::SharedSpace>();
-    M0_poly_coeff = Kokkos::View<Kokkos::View<double***,Kokkos::LayoutRight>*,Kokkos::SharedSpace>();
-    M0_poly_values = Kokkos::View<Kokkos::View<double***,Kokkos::LayoutRight>*,Kokkos::SharedSpace>();
-    M0_poly_adjoints = Kokkos::View<Kokkos::View<double***,Kokkos::LayoutRight>*,Kokkos::SharedSpace>();
+    M0_poly_coeff = Kokkos::View<Kokkos::View<Precision***,Kokkos::LayoutRight>*,Kokkos::SharedSpace>();
+    M0_poly_values = Kokkos::View<Kokkos::View<Precision***,Kokkos::LayoutRight>*,Kokkos::SharedSpace>();
+    M0_poly_adjoints = Kokkos::View<Kokkos::View<Precision***,Kokkos::LayoutRight>*,Kokkos::SharedSpace>();
 
-    A1_weights = Kokkos::View<Kokkos::View<double**,Kokkos::LayoutRight>*,Kokkos::SharedSpace>();
-    A1_weights_trans = Kokkos::View<Kokkos::View<double**,Kokkos::LayoutRight>*,Kokkos::SharedSpace>();
+    A1_weights = Kokkos::View<Kokkos::View<Precision**,Kokkos::LayoutRight>*,Kokkos::SharedSpace>();
+    A1_weights_trans = Kokkos::View<Kokkos::View<Precision**,Kokkos::LayoutRight>*,Kokkos::SharedSpace>();
 }
 
-void MACEKokkos::compute_node_energies_forces(
+template <typename Precision>
+void MACEKokkos<Precision>::compute_node_energies_forces(
     const int num_nodes,
     Kokkos::View<const int*> node_types,
     Kokkos::View<const int*> num_neigh,
@@ -94,7 +97,8 @@ void MACEKokkos::compute_node_energies_forces(
     reverse_A0(num_nodes, node_types, num_neigh, neigh_types, xyz, r);
 }
 
-void MACEKokkos::compute_R0(
+template <typename Precision>
+void MACEKokkos<Precision>::compute_R0(
     const int num_nodes,
     Kokkos::View<const int*> node_types,
     Kokkos::View<const int*> num_neigh,
@@ -168,7 +172,8 @@ void MACEKokkos::compute_R0(
     Kokkos::fence();
 }
 
-void MACEKokkos::compute_R1(
+template <typename Precision>
+void MACEKokkos<Precision>::compute_R1(
     const int num_nodes,
     Kokkos::View<const int*> node_types,
     Kokkos::View<const int*> num_neigh,
@@ -183,7 +188,8 @@ void MACEKokkos::compute_R1(
     Kokkos::fence();
 }
 
-void MACEKokkos::compute_Y(Kokkos::View<const double*> xyz) {
+template <typename Precision>
+void MACEKokkos<Precision>::compute_Y(Kokkos::View<const double*> xyz) {
 
 #ifndef SYMMETRIX_SPHERICART_CUDA
 
@@ -213,7 +219,7 @@ void MACEKokkos::compute_Y(Kokkos::View<const double*> xyz) {
     auto h_xyz = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), xyz_shuffled);
     auto h_Y = Kokkos::create_mirror_view(Y);
     auto h_Y_grad = Kokkos::create_mirror_view(Y_grad);
-    sphericart::SphericalHarmonics<double> sphericart(l_max);
+    sphericart::SphericalHarmonics<Precision> sphericart(l_max);
     sphericart.compute_array_with_gradients(
         h_xyz.data(), 3*num, h_Y.data(), 3*num*num_lm, h_Y_grad.data(), 3*num*num_lm),
     Kokkos::deep_copy(Y, h_Y);
@@ -250,15 +256,43 @@ void MACEKokkos::compute_Y(Kokkos::View<const double*> xyz) {
         Kokkos::realloc(Y, num*num_lm);
         Kokkos::realloc(Y_grad, 3*num*num_lm);
     }
-    sphericart::cuda::SphericalHarmonics<double> sphericart(l_max);
-    sphericart.compute_with_gradients(xyz.data(), num, Y.data(), Y_grad.data());
+    auto Y = this->Y;
+    auto Y_grad = this->Y_grad;
+
+    // shuffle to match e3nn
+    if (xyz_shuffled.extent(0) < 3*num)
+        Kokkos::realloc(xyz_shuffled, 3*num);
+    auto xyz_shuffled = this->xyz_shuffled;
+    Kokkos::parallel_for("shuffle_xyz", num, KOKKOS_LAMBDA (int i) {
+        xyz_shuffled(3*i) = xyz(3*i+2);
+        xyz_shuffled(3*i+1) = xyz(3*i);
+        xyz_shuffled(3*i+2) = xyz(3*i+1);
+    });
+    Kokkos::fence();
+
+    // call sphericart
+    sphericart::cuda::SphericalHarmonics<Precision> sphericart(l_max);
+    sphericart.compute_with_gradients(xyz_shuffled.data(), num, Y.data(), Y_grad.data());
+
+    // unshuffle gradient
+    if (Y_grad_shuffled.extent(0) < 3*num*num_lm)
+        Kokkos::realloc(Y_grad_shuffled, 3*num*num_lm);
+    Kokkos::deep_copy(Y_grad_shuffled, Y_grad);
+    auto Y_grad_shuffled = this->Y_grad_shuffled;
+    Kokkos::parallel_for("unshuffle_Y_grad", num, KOKKOS_LAMBDA (int i) {
+        for (int lm=0; lm<num_lm; ++lm) {
+            Y_grad(3*i*num_lm+0*num_lm+lm) = Y_grad_shuffled(3*i*num_lm+1*num_lm+lm);
+            Y_grad(3*i*num_lm+1*num_lm+lm) = Y_grad_shuffled(3*i*num_lm+2*num_lm+lm);
+            Y_grad(3*i*num_lm+2*num_lm+lm) = Y_grad_shuffled(3*i*num_lm+0*num_lm+lm);
+        }
+    });
+    Kokkos::fence();
+
     // normalize to match e3nn conventions
     const double normalization_factor = 2 * std::sqrt(M_PI);
-    auto Y = this->Y;
     Kokkos::parallel_for("normalize_Y", num*num_lm, KOKKOS_LAMBDA (int i) {
         Y(i) *= normalization_factor;
     });
-    auto Y_grad = this->Y_grad;
     Kokkos::parallel_for("normalize_Y_grad", 3*num*num_lm, KOKKOS_LAMBDA (int i) {
         Y_grad(i) *= normalization_factor;
     });
@@ -268,7 +302,8 @@ void MACEKokkos::compute_Y(Kokkos::View<const double*> xyz) {
     Kokkos::fence();
 }
 
-void MACEKokkos::compute_A0(
+template <typename Precision>
+void MACEKokkos<Precision>::compute_A0(
     const int num_nodes,
     View<const int*> node_types,
     View<const int*> num_neigh,
@@ -315,7 +350,8 @@ void MACEKokkos::compute_A0(
     Kokkos::fence();
 }
 
-void MACEKokkos::reverse_A0(
+template <typename Precision>
+void MACEKokkos<Precision>::reverse_A0(
     const int num_nodes,
     Kokkos::View<const int*> node_types,
     Kokkos::View<const int*> num_neigh,
@@ -354,8 +390,8 @@ void MACEKokkos::reverse_A0(
                 const double x_ij = xyz(3*ij) / r_ij;
                 const double y_ij = xyz(3*ij+1) / r_ij;
                 const double z_ij = xyz(3*ij+2) / r_ij;
-                const double* Y_ij = &Y(ij*num_lm);
-                const double* Y_grad_ij = &Y_grad(3*ij*num_lm);
+                const Precision* Y_ij = &Y(ij*num_lm);
+                const Precision* Y_grad_ij = &Y_grad(3*ij*num_lm);
                 double f_x, f_y, f_z;
                 Kokkos::parallel_reduce(
                     Kokkos::TeamThreadRange(team_member, num_lm),
@@ -382,7 +418,8 @@ void MACEKokkos::reverse_A0(
     Kokkos::fence();
 }
 
-void MACEKokkos::compute_A0_scaled(
+template <typename Precision>
+void MACEKokkos<Precision>::compute_A0_scaled(
     const int num_nodes,
     Kokkos::View<const int*> node_types,
     Kokkos::View<const int*> num_neigh,
@@ -437,7 +474,8 @@ void MACEKokkos::compute_A0_scaled(
     Kokkos::fence();
 }
 
-void MACEKokkos::reverse_A0_scaled(
+template <typename Precision>
+void MACEKokkos<Precision>::reverse_A0_scaled(
     const int num_nodes,
     Kokkos::View<const int*> node_types,
     Kokkos::View<const int*> num_neigh,
@@ -510,7 +548,8 @@ void MACEKokkos::reverse_A0_scaled(
 }
 
 #if 0
-void MACEKokkos::compute_M0(
+template <typename Precision>
+void MACEKokkos<Precision>::compute_M0(
     const int num_nodes,
     Kokkos::View<const int*> node_types)
 {
@@ -545,14 +584,15 @@ void MACEKokkos::compute_M0(
 #endif
 
 //#if 0
-void MACEKokkos::compute_M0(int num_nodes, Kokkos::View<const int*> node_types)
+template <typename Precision>
+void MACEKokkos<Precision>::compute_M0(int num_nodes, Kokkos::View<const int*> node_types)
 {
     if (M0.extent(0) < num_nodes)
         Kokkos::realloc(M0, num_nodes, num_LM, num_channels);
     Kokkos::deep_copy(M0, 0.0);
     for (int LM=0; LM<num_LM; ++LM) {
         if (M0_poly_values(LM).extent(0) < num_nodes)
-            M0_poly_values(LM) = Kokkos::View<double***,Kokkos::LayoutRight>(
+            M0_poly_values(LM) = Kokkos::View<Precision***,Kokkos::LayoutRight>(
                 Kokkos::view_alloc(std::string("M0_poly_values_")+std::to_string(LM),Kokkos::WithoutInitializing),
                 num_nodes, M0_poly_coeff(LM).extent(1), num_channels);
     }
@@ -599,7 +639,8 @@ void MACEKokkos::compute_M0(int num_nodes, Kokkos::View<const int*> node_types)
 //#endif
 
 #if 0
-void MACEKokkos::reverse_M0(
+template <typename Precision>
+void MACEKokkos<Precision>::reverse_M0(
     const int num_nodes,
     Kokkos::View<const int*> node_types)
 {
@@ -640,14 +681,15 @@ void MACEKokkos::reverse_M0(
 #endif
 
 //#if 0
-void MACEKokkos::reverse_M0(int num_nodes, Kokkos::View<const int*> node_types)
+template <typename Precision>
+void MACEKokkos<Precision>::reverse_M0(int num_nodes, Kokkos::View<const int*> node_types)
 {
     if (A0_adj.extent(0) < num_nodes)
         Kokkos::realloc(A0_adj, A0.extent(0), A0.extent(1), A0.extent(2));
     Kokkos::deep_copy(A0_adj, 0.0);
     for (int LM=0; LM<num_LM; ++LM) {
         if (M0_poly_adjoints(LM).extent(0) < num_nodes)
-            M0_poly_adjoints(LM) = Kokkos::View<double***,Kokkos::LayoutRight>(
+            M0_poly_adjoints(LM) = Kokkos::View<Precision***,Kokkos::LayoutRight>(
                 Kokkos::view_alloc(std::string("M0_poly_adjoints_")+std::to_string(LM),Kokkos::WithoutInitializing),
                 num_nodes, M0_poly_coeff(LM).extent(1), num_channels);
     }
@@ -702,7 +744,8 @@ void MACEKokkos::reverse_M0(int num_nodes, Kokkos::View<const int*> node_types)
 }
 //#endif
 
-void MACEKokkos::compute_H1(
+template <typename Precision>
+void MACEKokkos<Precision>::compute_H1(
     const int num_nodes)
 {
     if (H1.extent(0) < M0.extent(0))
@@ -730,7 +773,8 @@ void MACEKokkos::compute_H1(
     Kokkos::fence();
 }
 
-void MACEKokkos::reverse_H1(
+template <typename Precision>
+void MACEKokkos<Precision>::reverse_H1(
     const int num_nodes)
 {
     if (M0_adj.extent(0) < M0.extent(0))
@@ -758,7 +802,8 @@ void MACEKokkos::reverse_H1(
     Kokkos::fence();
 }
 
-void MACEKokkos::compute_Phi1(
+template <typename Precision>
+void MACEKokkos<Precision>::compute_Phi1(
     const int num_nodes,
     Kokkos::View<const int*> num_neigh,
     Kokkos::View<const int*> neigh_indices)
@@ -873,7 +918,8 @@ void MACEKokkos::compute_Phi1(
     Kokkos::fence();
 }
 
-void MACEKokkos::reverse_Phi1(
+template <typename Precision>
+void MACEKokkos<Precision>::reverse_Phi1(
     const int num_nodes,
     Kokkos::View<const int*> num_neigh,
     Kokkos::View<const int*> neigh_indices,
@@ -978,7 +1024,8 @@ void MACEKokkos::reverse_Phi1(
     Kokkos::fence();
 }
 
-void MACEKokkos::compute_A1(int num_nodes)
+template <typename Precision>
+void MACEKokkos<Precision>::compute_A1(int num_nodes)
 {
     // The core matrix multiplication is:
     //         [A1_il]_mk = \sum_(ek') [Phi1_il]_m(ek') [W_il]_(ek')k
@@ -1006,7 +1053,7 @@ void MACEKokkos::compute_A1(int num_nodes)
                 if (ll == l)
                     num_eta += 1;
             }
-            auto Phi1_il = Kokkos::View<double**,Kokkos::LayoutRight,Kokkos::MemoryUnmanaged>(
+            auto Phi1_il = Kokkos::View<Precision**,Kokkos::LayoutRight,Kokkos::MemoryUnmanaged>(
                 &Phi1(i,lme,0), 2*l+1, num_eta*num_channels);
             auto A1_il = Kokkos::subview(A1, i, Kokkos::make_pair(l*l,l*(l+2)+1), Kokkos::ALL);
             KokkosBatched::TeamGemm<Kokkos::TeamPolicy<>::member_type,
@@ -1018,7 +1065,8 @@ void MACEKokkos::compute_A1(int num_nodes)
     Kokkos::fence();
 }
 
-void MACEKokkos::reverse_A1(int num_nodes)
+template <typename Precision>
+void MACEKokkos<Precision>::reverse_A1(int num_nodes)
 {
     // The core matrix multiplication is:
     //         [dE/dPhi1_il]_m(ek) = \sum_k' [dE/dA1_il]_mk' [trans(W_il)]_k'(ek)
@@ -1047,7 +1095,7 @@ void MACEKokkos::reverse_A1(int num_nodes)
                     num_eta += 1;
             }
             auto dA1_il = Kokkos::subview(A1_adj, i, Kokkos::make_pair(l*l,l*l+2*l+1), Kokkos::ALL);
-            auto dPhi1_il = Kokkos::View<double**,Kokkos::LayoutRight,Kokkos::MemoryUnmanaged>(
+            auto dPhi1_il = Kokkos::View<Precision**,Kokkos::LayoutRight,Kokkos::MemoryUnmanaged>(
                 &dPhi1(i,lme,0), 2*l+1, num_eta*num_channels);
             KokkosBatched::TeamGemm<Kokkos::TeamPolicy<>::member_type,
                                     KokkosBatched::Trans::NoTranspose,
@@ -1058,7 +1106,8 @@ void MACEKokkos::reverse_A1(int num_nodes)
     Kokkos::fence();
 }
 
-void MACEKokkos::compute_A1_scaled(
+template <typename Precision>
+void MACEKokkos<Precision>::compute_A1_scaled(
     const int num_nodes,
     Kokkos::View<const int*> node_types,
     Kokkos::View<const int*> num_neigh,
@@ -1115,7 +1164,8 @@ void MACEKokkos::compute_A1_scaled(
     Kokkos::fence();
 }
 
-void MACEKokkos::reverse_A1_scaled(
+template <typename Precision>
+void MACEKokkos<Precision>::reverse_A1_scaled(
     const int num_nodes,
     Kokkos::View<const int*> node_types,
     Kokkos::View<const int*> num_neigh,
@@ -1189,7 +1239,8 @@ void MACEKokkos::reverse_A1_scaled(
 }
 
 #if 0
-void MACEKokkos::compute_M1(int num_nodes, Kokkos::View<const int*> node_types)
+template <typename Precision>
+void MACEKokkos<Precision>::compute_M1(int num_nodes, Kokkos::View<const int*> node_types)
 {
     Kokkos::realloc(M1, num_nodes, num_channels);
 
@@ -1220,7 +1271,8 @@ void MACEKokkos::compute_M1(int num_nodes, Kokkos::View<const int*> node_types)
 #endif
 
 //#if 0
-void MACEKokkos::compute_M1(int num_nodes, Kokkos::View<const int*> node_types)
+template <typename Precision>
+void MACEKokkos<Precision>::compute_M1(int num_nodes, Kokkos::View<const int*> node_types)
 {
     if (M1.extent(0) < num_nodes)
         Kokkos::realloc(M1, num_nodes, num_channels);
@@ -1269,7 +1321,8 @@ void MACEKokkos::compute_M1(int num_nodes, Kokkos::View<const int*> node_types)
 //#endif
 
 #if 0
-void MACEKokkos::reverse_M1(int num_nodes, Kokkos::View<const int*> node_types)
+template <typename Precision>
+void MACEKokkos<Precision>::reverse_M1(int num_nodes, Kokkos::View<const int*> node_types)
 {
     Kokkos::realloc(A1_adj, A1.extent(0), A1.extent(1), A1.extent(2));
     Kokkos::deep_copy(A1_adj, 0.0);
@@ -1305,7 +1358,8 @@ void MACEKokkos::reverse_M1(int num_nodes, Kokkos::View<const int*> node_types)
 }
 #endif
 
-void MACEKokkos::reverse_M1(int num_nodes, Kokkos::View<const int*> node_types)
+template <typename Precision>
+void MACEKokkos<Precision>::reverse_M1(int num_nodes, Kokkos::View<const int*> node_types)
 {
     if (A1_adj.extent(0) < num_nodes)
         Kokkos::realloc(A1_adj, A1.extent(0), A1.extent(1), A1.extent(2));
@@ -1362,7 +1416,8 @@ void MACEKokkos::reverse_M1(int num_nodes, Kokkos::View<const int*> node_types)
     Kokkos::fence();
 }
 
-void MACEKokkos::compute_H2(int num_nodes, Kokkos::View<const int*> node_types)
+template <typename Precision>
+void MACEKokkos<Precision>::compute_H2(int num_nodes, Kokkos::View<const int*> node_types)
 {
     if (H2.extent(0) < num_nodes or H2.extent(1) != num_channels)
         Kokkos::realloc(H2, num_nodes, num_channels);
@@ -1399,7 +1454,8 @@ void MACEKokkos::compute_H2(int num_nodes, Kokkos::View<const int*> node_types)
     Kokkos::fence();
 }
 
-void MACEKokkos::reverse_H2(int num_nodes, Kokkos::View<const int*> node_types, bool zero_H1_adj)
+template <typename Precision>
+void MACEKokkos<Precision>::reverse_H2(int num_nodes, Kokkos::View<const int*> node_types, bool zero_H1_adj)
 {
     if (H1_adj.extent(0) < H1.extent(0))
         Kokkos::resize(H1_adj, H1.extent(0), H1.extent(1), H1.extent(2));
@@ -1431,7 +1487,8 @@ void MACEKokkos::reverse_H2(int num_nodes, Kokkos::View<const int*> node_types,
     Kokkos::fence();
 }
 
-double MACEKokkos::compute_readouts(int num_nodes, const Kokkos::View<const int*> node_types)
+template <typename Precision>
+double MACEKokkos<Precision>::compute_readouts(int num_nodes, const Kokkos::View<const int*> node_types)
 {
     if (H1_adj.extent(0) < H1.extent(0))
         Kokkos::realloc(H1_adj, H1.extent(0), H1.extent(1), H1.extent(2));
@@ -1484,7 +1541,8 @@ double MACEKokkos::compute_readouts(int num_nodes, const Kokkos::View<const int*
     return energy;
 }
 
-void MACEKokkos::load_from_json(std::string filename)
+template <typename Precision>
+void MACEKokkos<Precision>::load_from_json(std::string filename)
 {
     std::ifstream f(filename);
     nlohmann::json file = nlohmann::json::parse(f);
@@ -1514,7 +1572,7 @@ void MACEKokkos::load_from_json(std::string filename)
     const double spl_h = file["radial_spline_h"];
     auto spl_values_0 = file["radial_spline_values_0"].get<std::vector<std::vector<std::vector<double>>>>();
     auto spl_derivs_0 = file["radial_spline_derivs_0"].get<std::vector<std::vector<std::vector<double>>>>();
-    auto c = Kokkos::View<double****,Kokkos::LayoutRight>(
+    auto c = Kokkos::View<Precision****,Kokkos::LayoutRight>(
         "c", atomic_numbers.size()*atomic_numbers.size(), spl_values_0[0][0].size()-1, 4, (l_max+1)*num_channels);
     auto h_c = Kokkos::create_mirror_view(c);
     for (int a=0; a<atomic_numbers.size(); ++a) {
@@ -1577,7 +1635,7 @@ void MACEKokkos::load_from_json(std::string filename)
     // R1
     auto spl_values_1 = file["radial_spline_values_1"].get<std::vector<std::vector<std::vector<double>>>>();
     auto spl_derivs_1 = file["radial_spline_derivs_1"].get<std::vector<std::vector<std::vector<double>>>>();
-    radial_1 = RadialFunctionSetKokkos(spl_h, spl_values_1, spl_derivs_1);
+    radial_1 = RadialFunctionSetKokkos<Precision>(spl_h, spl_values_1, spl_derivs_1);
 
     // A0 scaling
     A0_scaled = file["A0_scaled"].get<bool>();
@@ -1589,18 +1647,18 @@ void MACEKokkos::load_from_json(std::string filename)
         auto A0_spline_derivs = std::vector<std::vector<std::vector<double>>>();
         for (auto& derivs : file["A0_spline_derivs"].get<std::vector<std::vector<double>>>())
             A0_spline_derivs.push_back({derivs});  // adds dimension to reach 3d
-        A0_splines = RadialFunctionSetKokkos(A0_spline_h, A0_spline_values, A0_spline_derivs);
+        A0_splines = RadialFunctionSetKokkos<double>(A0_spline_h, A0_spline_values, A0_spline_derivs);
     }
 
     // M0 weights and monomials
     auto M0_weights_file = file["M0_weights"].get<std::map<std::string,std::map<std::string,std::map<std::string,std::vector<double>>>>>();
     auto M0_monomials_file = file["M0_monomials"].get<std::map<std::string,std::vector<std::vector<int>>>>();
-    M0_weights = Kokkos::View<Kokkos::View<double***,Kokkos::LayoutRight>*,Kokkos::SharedSpace>(
+    M0_weights = Kokkos::View<Kokkos::View<Precision***,Kokkos::LayoutRight>*,Kokkos::SharedSpace>(
         Kokkos::view_alloc("M0_weights", Kokkos::SequentialHostInit), num_LM);
     M0_monomials = Kokkos::View<Kokkos::View<int**,Kokkos::LayoutRight>*,Kokkos::SharedSpace>(
         Kokkos::view_alloc("M0_monomials", Kokkos::SequentialHostInit), num_LM);
     for (int LM=0; LM<num_LM; ++LM) {
-        M0_weights(LM) = Kokkos::View<double***,Kokkos::LayoutRight>(
+        M0_weights(LM) = Kokkos::View<Precision***,Kokkos::LayoutRight>(
             Kokkos::view_alloc(std::string("M0_weights_") + std::to_string(LM), Kokkos::WithoutInitializing),
             atomic_numbers.size(), num_channels, M0_monomials_file[std::to_string(LM)].size());
         auto h_M0_weights_LM = Kokkos::create_mirror_view(M0_weights(LM));
@@ -1641,14 +1699,14 @@ void MACEKokkos::load_from_json(std::string filename)
         Kokkos::deep_copy(M0_poly_spec(LM), h_M0_poly_spec_LM);
     }
     // M0_poly_coeff
-    M0_poly_coeff = Kokkos::View<Kokkos::View<double***,Kokkos::LayoutRight>*,Kokkos::SharedSpace>(
+    M0_poly_coeff = Kokkos::View<Kokkos::View<Precision***,Kokkos::LayoutRight>*,Kokkos::SharedSpace>(
         Kokkos::view_alloc("M0_poly_coeff",Kokkos::SequentialHostInit), num_LM);
     for (int LM=0; LM<num_LM; ++LM) {
         auto P = MultivariatePolynomial(
             num_lm,
             M0_weights_file[std::to_string(0)][std::to_string(LM)][std::to_string(0)],
             M0_monomials_file[std::to_string(LM)]);
-        M0_poly_coeff(LM) = Kokkos::View<double***,Kokkos::LayoutRight>(
+        M0_poly_coeff(LM) = Kokkos::View<Precision***,Kokkos::LayoutRight>(
             Kokkos::view_alloc(std::string("M0_poly_coeff_")+std::to_string(LM),Kokkos::WithoutInitializing),
             atomic_numbers.size(), P.node_coefficients.size(), num_channels);
         auto h_M0_poly_coeff_LM = Kokkos::create_mirror_view(M0_poly_coeff(LM));
@@ -1665,15 +1723,15 @@ void MACEKokkos::load_from_json(std::string filename)
         }
         Kokkos::deep_copy(M0_poly_coeff(LM), h_M0_poly_coeff_LM);
     }
-    M0_poly_values = Kokkos::View<Kokkos::View<double***,Kokkos::LayoutRight>*,Kokkos::SharedSpace>(
+    M0_poly_values = Kokkos::View<Kokkos::View<Precision***,Kokkos::LayoutRight>*,Kokkos::SharedSpace>(
         Kokkos::view_alloc("M0_poly_values",Kokkos::SequentialHostInit), num_LM);
-    M0_poly_adjoints = Kokkos::View<Kokkos::View<double***,Kokkos::LayoutRight>*,Kokkos::SharedSpace>(
+    M0_poly_adjoints = Kokkos::View<Kokkos::View<Precision***,Kokkos::LayoutRight>*,Kokkos::SharedSpace>(
         Kokkos::view_alloc("M0_poly_adjoints",Kokkos::SequentialHostInit), num_LM);
 
     // H1 weights
     set_kokkos_view(
         H1_weights,
-        file["H1_weights"].get<std::vector<double>>(),
+        file["H1_weights"].get<std::vector<Precision>>(),
         L_max+1,
         num_channels,
         num_channels);
@@ -1683,7 +1741,7 @@ void MACEKokkos::load_from_json(std::string filename)
     Phi1_l1 = toKokkosView("Phi1_l1", file["Phi1_l1"].get<std::vector<int>>());
     Phi1_l2 = toKokkosView("Phi1_l2", file["Phi1_l2"].get<std::vector<int>>());
     Phi1_lme = toKokkosView("Phi1_lme", file["Phi1_lme"].get<std::vector<int>>());
-    Phi1_clebsch_gordan = toKokkosView("Phi1_clebsch_gordan", file["Phi1_clebsch_gordan"].get<std::vector<double>>());
+    Phi1_clebsch_gordan = toKokkosView("Phi1_clebsch_gordan", file["Phi1_clebsch_gordan"].get<std::vector<Precision>>());
     Phi1_lelm1lm2 = toKokkosView("Phi1_lelm1lm2", file["Phi1_lelm1lm2"].get<std::vector<int>>());
     num_lme = 0;
     auto h_Phi1_l = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), Phi1_l);
@@ -1715,20 +1773,20 @@ void MACEKokkos::load_from_json(std::string filename)
     this->Phi1_lel1l2 = toKokkosView("Phi1_lel1l2", Phi1_lel1l2);
 
     // A1 weights
-    auto file_A1_weights = file["A1_weights"].get<std::vector<std::vector<double>>>();
-    A1_weights = Kokkos::View<Kokkos::View<double**,Kokkos::LayoutRight>*,Kokkos::SharedSpace>(
+    auto file_A1_weights = file["A1_weights"].get<std::vector<std::vector<Precision>>>();
+    A1_weights = Kokkos::View<Kokkos::View<Precision**,Kokkos::LayoutRight>*,Kokkos::SharedSpace>(
         Kokkos::view_alloc("A1_weights", Kokkos::SequentialHostInit), l_max+1);
-    A1_weights_trans = Kokkos::View<Kokkos::View<double**,Kokkos::LayoutRight>*,Kokkos::SharedSpace>(
+    A1_weights_trans = Kokkos::View<Kokkos::View<Precision**,Kokkos::LayoutRight>*,Kokkos::SharedSpace>(
         Kokkos::view_alloc("A1_weights_trans", Kokkos::SequentialHostInit), l_max+1);
     for (int l=0; l<=l_max; ++l) {
         int num_eta = 0;
         auto h_Phi1_l = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), Phi1_l);
         for (int i=0; i<h_Phi1_l.size(); ++i)
             num_eta += (h_Phi1_l(i) == l);
-        A1_weights(l) = Kokkos::View<double**,Kokkos::LayoutRight>(
+        A1_weights(l) = Kokkos::View<Precision**,Kokkos::LayoutRight>(
             Kokkos::view_alloc(std::string("A1_weights_") + std::to_string(l), Kokkos::WithoutInitializing),
             num_eta*num_channels, num_channels);
-        A1_weights_trans(l) = Kokkos::View<double**,Kokkos::LayoutRight>(
+        A1_weights_trans(l) = Kokkos::View<Precision**,Kokkos::LayoutRight>(
             Kokkos::view_alloc(std::string("A1_weights_") + std::to_string(l), Kokkos::WithoutInitializing),
             num_channels, num_eta*num_channels);
         auto h_A1_weights_l = Kokkos::create_mirror_view(A1_weights(l));
@@ -1753,7 +1811,7 @@ void MACEKokkos::load_from_json(std::string filename)
         auto A1_spline_derivs = std::vector<std::vector<std::vector<double>>>();
         for (auto& derivs : file["A1_spline_derivs"].get<std::vector<std::vector<double>>>())
             A1_spline_derivs.push_back({derivs});  // adds dimension to reach 3d
-        A1_splines = RadialFunctionSetKokkos(A1_spline_h, A1_spline_values, A1_spline_derivs);
+        A1_splines = RadialFunctionSetKokkos<double>(A1_spline_h, A1_spline_values, A1_spline_derivs);
     }
 
     // M1 weights and monomials
@@ -1827,3 +1885,6 @@ void MACEKokkos::load_from_json(std::string filename)
         std::vector<std::vector<double>>{readout_2_weights_1, readout_2_weights_2},
         file["readout_2_scale_factor"]);
 }
+
+template class MACEKokkos<float>;
+template class MACEKokkos<double>;
diff --git a/libsymmetrix/source/mace_kokkos.hpp b/libsymmetrix/source/mace_kokkos.hpp
index e2de6ba..d827e90 100644
--- a/libsymmetrix/source/mace_kokkos.hpp
+++ b/libsymmetrix/source/mace_kokkos.hpp
@@ -13,6 +13,7 @@
 #include "radial_function_set_kokkos.hpp"
 #include "zbl_kokkos.hpp"
 
+template <typename Precision>
 class MACEKokkos {
 
 public:
@@ -45,8 +46,8 @@ ZBLKokkos zbl;
 
 // R0
 double R0_spline_h;
-Kokkos::View<const double****,Kokkos::LayoutRight> R0_spline_coefficients;
-Kokkos::View<double**,Kokkos::LayoutRight> R0, R0_deriv;
+Kokkos::View<const Precision****,Kokkos::LayoutRight> R0_spline_coefficients;
+Kokkos::View<Precision**,Kokkos::LayoutRight> R0, R0_deriv;
 void compute_R0(const int num_nodes,
                 Kokkos::View<const int*> node_types,
                 Kokkos::View<const int*> num_neigh,
@@ -54,8 +55,8 @@ void compute_R0(const int num_nodes,
                 Kokkos::View<const double*> r);
 
 // R1
-RadialFunctionSetKokkos radial_1;
-Kokkos::View<double**,Kokkos::LayoutRight> R1, R1_deriv;
+RadialFunctionSetKokkos<Precision> radial_1;
+Kokkos::View<Precision**,Kokkos::LayoutRight> R1, R1_deriv;
 void compute_R1(const int num_nodes,
                 Kokkos::View<const int*> node_types,
                 Kokkos::View<const int*> num_neigh,
@@ -63,13 +64,13 @@ void compute_R1(const int num_nodes,
                 Kokkos::View<const double*> r);
 
 // Spherical harmonics
-Kokkos::View<double*> xyz_shuffled;
-Kokkos::View<double*> Y, Y_grad;// TODO: make multidimensional
-Kokkos::View<double*> Y_grad_shuffled;
+Kokkos::View<Precision*> xyz_shuffled;
+Kokkos::View<Precision*> Y, Y_grad;// TODO: make multidimensional
+Kokkos::View<Precision*> Y_grad_shuffled;
 void compute_Y(Kokkos::View<const double*> xyz);
 
 // A0
-Kokkos::View<double***,Kokkos::LayoutRight> A0, A0_adj;
+Kokkos::View<Precision***,Kokkos::LayoutRight> A0, A0_adj;
 void compute_A0(const int num_nodes,
                 Kokkos::View<const int*> node_types,
                 Kokkos::View<const int*> num_neigh,
@@ -83,7 +84,7 @@ void reverse_A0(const int num_nodes,
 
 // A0 rescaling
 bool A0_scaled;
-RadialFunctionSetKokkos A0_splines;
+RadialFunctionSetKokkos<double> A0_splines;
 Kokkos::View<double**,Kokkos::LayoutRight> A0_spline_values;
 Kokkos::View<double**,Kokkos::LayoutRight> A0_spline_derivs;
 void compute_A0_scaled(
@@ -101,19 +102,19 @@ void reverse_A0_scaled(
     Kokkos::View<const double*> r);
 
 // M0
-Kokkos::View<double***,Kokkos::LayoutRight> M0, M0_adj;
+Kokkos::View<Precision***,Kokkos::LayoutRight> M0, M0_adj;
 Kokkos::View<Kokkos::View<int**,Kokkos::LayoutRight>*,Kokkos::SharedSpace> M0_monomials;
-Kokkos::View<Kokkos::View<double***,Kokkos::LayoutRight>*,Kokkos::SharedSpace> M0_weights;
+Kokkos::View<Kokkos::View<Precision***,Kokkos::LayoutRight>*,Kokkos::SharedSpace> M0_weights;
 Kokkos::View<Kokkos::View<int**,Kokkos::LayoutRight>*,Kokkos::SharedSpace> M0_poly_spec;
-Kokkos::View<Kokkos::View<double***,Kokkos::LayoutRight>*,Kokkos::SharedSpace> M0_poly_coeff;
-Kokkos::View<Kokkos::View<double***,Kokkos::LayoutRight>*,Kokkos::SharedSpace> M0_poly_values;
-Kokkos::View<Kokkos::View<double***,Kokkos::LayoutRight>*,Kokkos::SharedSpace> M0_poly_adjoints;
+Kokkos::View<Kokkos::View<Precision***,Kokkos::LayoutRight>*,Kokkos::SharedSpace> M0_poly_coeff;
+Kokkos::View<Kokkos::View<Precision***,Kokkos::LayoutRight>*,Kokkos::SharedSpace> M0_poly_values;
+Kokkos::View<Kokkos::View<Precision***,Kokkos::LayoutRight>*,Kokkos::SharedSpace> M0_poly_adjoints;
 void compute_M0(const int num_nodes, Kokkos::View<const int*> node_types);
 void reverse_M0(const int num_nodes, Kokkos::View<const int*> node_types);
 
 // H1
-Kokkos::View<double***,Kokkos::LayoutRight> H1, H1_adj;
-Kokkos::View<double***,Kokkos::LayoutRight> H1_weights;
+Kokkos::View<Precision***,Kokkos::LayoutRight> H1, H1_adj;
+Kokkos::View<Precision***,Kokkos::LayoutRight> H1_weights;
 void compute_H1(const int num_nodes);
 void reverse_H1(const int num_nodes);
 
@@ -121,9 +122,9 @@ void reverse_H1(const int num_nodes);
 int num_lelm1lm2, num_lme;
 Kokkos::View<int*> Phi1_l, Phi1_l1, Phi1_l2;
 Kokkos::View<int*> Phi1_lme, Phi1_lelm1lm2;
-Kokkos::View<double*> Phi1_clebsch_gordan;
-Kokkos::View<double***,Kokkos::LayoutRight> Phi1r, dPhi1r;
-Kokkos::View<double***,Kokkos::LayoutRight> Phi1, dPhi1;
+Kokkos::View<Precision*> Phi1_clebsch_gordan;
+Kokkos::View<Precision***,Kokkos::LayoutRight> Phi1r, dPhi1r;
+Kokkos::View<Precision***,Kokkos::LayoutRight> Phi1, dPhi1;
 void compute_Phi1(const int num_nodes, Kokkos::View<const int*> num_neigh, Kokkos::View<const int*> neigh_indices);
 void reverse_Phi1(const int num_nodes, Kokkos::View<const int*> num_neigh, Kokkos::View<const int*> neigh_indices, Kokkos::View<const double*> xyz, Kokkos::View<const double*> r, bool zero_dxyz = true, bool zero_H1_adj = true);
 
@@ -131,15 +132,15 @@ void reverse_Phi1(const int num_nodes, Kokkos::View<const int*> num_neigh, Kokko
 Kokkos::View<int*> Phi1_lm1, Phi1_lm2, Phi1_lel1l2;
 
 // A1
-Kokkos::View<double***,Kokkos::LayoutRight> A1, A1_adj;
-Kokkos::View<Kokkos::View<double**,Kokkos::LayoutRight>*,Kokkos::SharedSpace> A1_weights;
-Kokkos::View<Kokkos::View<double**,Kokkos::LayoutRight>*,Kokkos::SharedSpace> A1_weights_trans;
+Kokkos::View<Precision***,Kokkos::LayoutRight> A1, A1_adj;
+Kokkos::View<Kokkos::View<Precision**,Kokkos::LayoutRight>*,Kokkos::SharedSpace> A1_weights;
+Kokkos::View<Kokkos::View<Precision**,Kokkos::LayoutRight>*,Kokkos::SharedSpace> A1_weights_trans;
 void compute_A1(int num_nodes);
 void reverse_A1(int num_nodes);
 
 // A1 rescaling
 bool A1_scaled;
-RadialFunctionSetKokkos A1_splines;
+RadialFunctionSetKokkos<double> A1_splines;
 Kokkos::View<double**,Kokkos::LayoutRight> A1_spline_values;
 Kokkos::View<double**,Kokkos::LayoutRight> A1_spline_derivs;
 void compute_A1_scaled(
@@ -157,13 +158,13 @@ void reverse_A1_scaled(
     Kokkos::View<const double*> r);
 
 // M1
-Kokkos::View<double**,Kokkos::LayoutRight> M1, M1_adj;
+Kokkos::View<Precision**,Kokkos::LayoutRight> M1, M1_adj;
 Kokkos::View<int**,Kokkos::LayoutRight> M1_monomials;
-Kokkos::View<double***,Kokkos::LayoutRight> M1_weights;
+Kokkos::View<Precision***,Kokkos::LayoutRight> M1_weights;
 Kokkos::View<int**,Kokkos::LayoutRight> M1_poly_spec;
-Kokkos::View<double***,Kokkos::LayoutRight> M1_poly_coeff;
-Kokkos::View<double***,Kokkos::LayoutRight> M1_poly_values;
-Kokkos::View<double***,Kokkos::LayoutRight> M1_poly_adjoints;
+Kokkos::View<Precision***,Kokkos::LayoutRight> M1_poly_coeff;
+Kokkos::View<Precision***,Kokkos::LayoutRight> M1_poly_values;
+Kokkos::View<Precision***,Kokkos::LayoutRight> M1_poly_adjoints;
 void compute_M1(int num_nodes, Kokkos::View<const int*> node_types);
 void reverse_M1(int num_nodes, Kokkos::View<const int*> node_types);
 
diff --git a/libsymmetrix/source/radial_function_set_kokkos.cpp b/libsymmetrix/source/radial_function_set_kokkos.cpp
index 56814b3..ae0dec3 100644
--- a/libsymmetrix/source/radial_function_set_kokkos.cpp
+++ b/libsymmetrix/source/radial_function_set_kokkos.cpp
@@ -1,14 +1,15 @@
 #include <vector>
+#include <cmath>
 
 #include "radial_function_set_kokkos.hpp"
 
-
-RadialFunctionSetKokkos::RadialFunctionSetKokkos()
+template <typename Precision>
+RadialFunctionSetKokkos<Precision>::RadialFunctionSetKokkos()
 {
 }
 
-
-RadialFunctionSetKokkos::RadialFunctionSetKokkos(
+template <typename Precision>
+RadialFunctionSetKokkos<Precision>::RadialFunctionSetKokkos(
     double h,
     std::vector<std::vector<std::vector<double>>> node_values,
     std::vector<std::vector<std::vector<double>>> node_derivatives)
@@ -19,18 +20,20 @@ RadialFunctionSetKokkos::RadialFunctionSetKokkos(
     num_functions = node_values[0].size();
     num_nodes = node_values[0][0].size();
 
-    auto c = Kokkos::View<double****, Kokkos::LayoutRight>(
+    auto c = Kokkos::View<Precision****, Kokkos::LayoutRight>(
         "coefficients", num_edge_types, num_nodes-1, 4, num_functions);
     auto h_c = Kokkos::create_mirror_view(c);
     for (int a=0; a<num_edge_types; ++a) {
         for (int i=0; i<num_nodes-1; ++i) {
             for (int j=0; j<num_functions; ++j) {
-                h_c(a,i,0,j) = node_values[a][j][i];
-                h_c(a,i,1,j) = node_derivatives[a][j][i];
-                h_c(a,i,2,j) = (-3*node_values[a][j][i] -2*h*node_derivatives[a][j][i]
-                                + 3*node_values[a][j][i+1] - h*node_derivatives[a][j][i+1]) / (h*h);
-                h_c(a,i,3,j) = (2*node_values[a][j][i] + h*node_derivatives[a][j][i]
-                                - 2*node_values[a][j][i+1] + h*node_derivatives[a][j][i+1]) / (h*h*h);
+                h_c(a,i,0,j) = static_cast<Precision>(node_values[a][j][i]);
+                h_c(a,i,1,j) = static_cast<Precision>(node_derivatives[a][j][i]);
+                h_c(a,i,2,j) = static_cast<Precision>(
+                    (-3*node_values[a][j][i] -2*h*node_derivatives[a][j][i]
+                        + 3*node_values[a][j][i+1] - h*node_derivatives[a][j][i+1]) / (h*h));
+                h_c(a,i,3,j) = static_cast<Precision>(
+                    (2*node_values[a][j][i] + h*node_derivatives[a][j][i]
+                        - 2*node_values[a][j][i+1] + h*node_derivatives[a][j][i+1]) / (h*h*h));
             }
         }
     }
@@ -41,14 +44,15 @@ RadialFunctionSetKokkos::RadialFunctionSetKokkos(
 
 // This older, cleaner version of the function lacks edge-type dependence
 #if 0
-void RadialFunctionSetKokkos::evaluate(
+template <typename Precision>
+void RadialFunctionSetKokkos<Precision>::evaluate(
     const int num_nodes,
     Kokkos::View<const int*> node_types,
     Kokkos::View<const int*> num_neigh,
     Kokkos::View<const int*> neigh_types,
     Kokkos::View<const double*> r,
-    Kokkos::View<double**,Kokkos::LayoutRight> R,
-    Kokkos::View<double**,Kokkos::LayoutRight> R_deriv) const
+    Kokkos::View<Precision**,Kokkos::LayoutRight> R,
+    Kokkos::View<Precision**,Kokkos::LayoutRight> R_deriv) const
 {
     const auto h = this->h;
     const auto num_functions = this->num_functions;
@@ -117,14 +121,15 @@ void RadialFunctionSetKokkos::evaluate(
 #endif
 
 
-void RadialFunctionSetKokkos::evaluate(
+template <typename Precision>
+void RadialFunctionSetKokkos<Precision>::evaluate(
     const int num_nodes,
     Kokkos::View<const int*> node_types,
     Kokkos::View<const int*> num_neigh,
     Kokkos::View<const int*> neigh_types,
     Kokkos::View<const double*> r,
-    Kokkos::View<double**,Kokkos::LayoutRight> R,
-    Kokkos::View<double**,Kokkos::LayoutRight> R_deriv) const
+    Kokkos::View<Precision**,Kokkos::LayoutRight> R,
+    Kokkos::View<Precision**,Kokkos::LayoutRight> R_deriv) const
 {
     const auto h = this->h;
     const auto num_functions = this->num_functions;
@@ -182,13 +187,20 @@ void RadialFunctionSetKokkos::evaluate(
             Kokkos::parallel_for(
                 Kokkos::TeamVectorRange(team_member, num_functions),
                 [&] (const int k) {
-                    const double c0 = c(type_ij,n,0,k);
-                    const double c1 = c(type_ij,n,1,k); 
-                    const double c2 = c(type_ij,n,2,k); 
-                    const double c3 = c(type_ij,n,3,k); 
-                    R(ij,k) = c0 + c1*x + c2*xx + c3*xxx;
-                    R_deriv(ij,k) = c1 + c2*two_x + c3*three_xx;
+                    const Precision c0 = c(type_ij,n,0,k);
+                    const Precision c1 = c(type_ij,n,1,k);
+                    const Precision c2 = c(type_ij,n,2,k);
+                    const Precision c3 = c(type_ij,n,3,k);
+                    R(ij,k) = c0 + c1*static_cast<Precision>(x)
+                        + c2*static_cast<Precision>(xx)
+                        + c3*static_cast<Precision>(xxx);
+                    R_deriv(ij,k) = c1
+                        + c2*static_cast<Precision>(two_x)
+                        + c3*static_cast<Precision>(three_xx);
                 });
         });
         Kokkos::fence();
 }
+
+template class RadialFunctionSetKokkos<float>;
+template class RadialFunctionSetKokkos<double>;
diff --git a/libsymmetrix/source/radial_function_set_kokkos.hpp b/libsymmetrix/source/radial_function_set_kokkos.hpp
index 44166f3..4859180 100644
--- a/libsymmetrix/source/radial_function_set_kokkos.hpp
+++ b/libsymmetrix/source/radial_function_set_kokkos.hpp
@@ -4,6 +4,7 @@
 
 #include <Kokkos_Core.hpp>
 
+template <typename Precision>
 class RadialFunctionSetKokkos
 {
 public:
@@ -19,8 +20,8 @@ class RadialFunctionSetKokkos
         Kokkos::View<const int*> num_neigh,
         Kokkos::View<const int*> neigh_types,
         Kokkos::View<const double*> r,
-        Kokkos::View<double**,Kokkos::LayoutRight> R,
-        Kokkos::View<double**,Kokkos::LayoutRight> R_deriv) const;
+        Kokkos::View<Precision**,Kokkos::LayoutRight> R,
+        Kokkos::View<Precision**,Kokkos::LayoutRight> R_deriv) const;
     
 private:
     
@@ -28,5 +29,5 @@ class RadialFunctionSetKokkos
     int num_edge_types;
     int num_functions;
     int num_nodes;
-    Kokkos::View<const double****,Kokkos::LayoutRight> coefficients;
+    Kokkos::View<const Precision****,Kokkos::LayoutRight> coefficients;
 };
diff --git a/pair_symmetrix/pair_symmetrix_mace_kokkos.cpp b/pair_symmetrix/pair_symmetrix_mace_kokkos.cpp
index ee62bc4..bb9ebc0 100644
--- a/pair_symmetrix/pair_symmetrix_mace_kokkos.cpp
+++ b/pair_symmetrix/pair_symmetrix_mace_kokkos.cpp
@@ -134,7 +134,7 @@ void PairSymmetrixMACEKokkos<DeviceType>::coeff(int narg, char **arg)
   if (!allocated) allocate();
 
   utils::logmesg(lmp, "Loading MACEKokkos model from \'{}\' ... ", arg[2]);
-  mace = std::make_unique<MACEKokkos>(arg[2]);
+  mace = std::make_unique<MACEKokkos<double>>(arg[2]);
   utils::logmesg(lmp, "success\n");
 
   // extract atomic numbers from pair_coeff
diff --git a/pair_symmetrix/pair_symmetrix_mace_kokkos.h b/pair_symmetrix/pair_symmetrix_mace_kokkos.h
index 384d2d2..c83479d 100644
--- a/pair_symmetrix/pair_symmetrix_mace_kokkos.h
+++ b/pair_symmetrix/pair_symmetrix_mace_kokkos.h
@@ -56,7 +56,7 @@ class PairSymmetrixMACEKokkos : public Pair, public KokkosBase {
 
  protected:
   std::string mode;
-  std::unique_ptr<MACEKokkos> mace;
+  std::unique_ptr<MACEKokkos<double>> mace;
   Kokkos::View<int*> mace_types;
   Kokkos::View<double***,Kokkos::LayoutRight> H1, H1_adj;
 
diff --git a/symmetrix/source/cpp/mace_kokkos.cpp b/symmetrix/source/cpp/mace_kokkos.cpp
index 811f978..b7d28d1 100644
--- a/symmetrix/source/cpp/mace_kokkos.cpp
+++ b/symmetrix/source/cpp/mace_kokkos.cpp
@@ -8,33 +8,35 @@
 
 namespace py = pybind11;
 
-void bind_mace_kokkos(py::module_ &m)
+template <typename Precision>
+void bind_mace_kokkos(py::module_ &m, const char* class_name)
 {
-    py::class_<MACEKokkos>(m, "MACEKokkos")
+    py::class_<MACEKokkos<Precision>>(m, class_name)
         .def(py::init<std::string>())
         .def_property_readonly("atomic_numbers",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.atomic_numbers);
             })
+        .def_readonly("r_cut", &MACEKokkos<Precision>::r_cut)
         // node energies
         .def_property("node_energies",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.node_energies);
             },
-            [] (MACEKokkos& self, py::array_t<double> node_energies) {
+            [] (MACEKokkos<Precision>& self, py::array_t<double> node_energies) {
                 set_kokkos_view(self.node_energies, node_energies);
             })
         // partial forces
         .def_property("node_forces",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.node_forces);
             },
-            [] (MACEKokkos& self, py::array_t<double> node_forces) {
+            [] (MACEKokkos<Precision>& self, py::array_t<double> node_forces) {
                 set_kokkos_view(self.node_forces, node_forces);
             })
         // node energies and forces
         .def("compute_node_energies_forces",
-            [] (MACEKokkos& self,
+            [] (MACEKokkos<Precision>& self,
                     const int num_nodes,
                     py::array_t<int> node_types,
                     py::array_t<int> num_neigh,
@@ -53,15 +55,15 @@ void bind_mace_kokkos(py::module_ &m)
             })
         // R0
         .def_property("R0",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.R0);
             },
-            [] (MACEKokkos& self, py::array_t<double> R0) {
+            [] (MACEKokkos<Precision>& self, py::array_t<Precision> R0) {
                 const int total_num_neigh = R0.size()/((self.l_max+1)*self.num_channels);
                 set_kokkos_view(self.R0, R0, total_num_neigh, (self.l_max+1)*self.num_channels);
             })
         .def("compute_R0", 
-            [](MACEKokkos& self,
+            [](MACEKokkos<Precision>& self,
                     const int num_nodes,
                     py::array_t<int> node_types,
                     py::array_t<int> num_neigh,
@@ -76,16 +78,16 @@ void bind_mace_kokkos(py::module_ &m)
             })
         // R1
         .def_property("R1",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.R1);
             },
-            [] (MACEKokkos& self, py::array_t<double> R1) {
+            [] (MACEKokkos<Precision>& self, py::array_t<Precision> R1) {
                 const int num_le = self.Phi1_l.size();
                 const int total_num_neigh = R1.size()/(num_le*self.num_channels);
                 set_kokkos_view(self.R1, R1, total_num_neigh, num_le*self.num_channels);
             })
         .def("compute_R1", 
-            [](MACEKokkos& self,
+            [](MACEKokkos<Precision>& self,
                     const int num_nodes,
                     py::array_t<int> node_types,
                     py::array_t<int> num_neigh,
@@ -100,30 +102,30 @@ void bind_mace_kokkos(py::module_ &m)
             })
         // Y
         .def("compute_Y",
-            [] (MACEKokkos& self,
+            [] (MACEKokkos<Precision>& self,
                     py::array_t<double> xyz) {
                 self.compute_Y(
                     create_kokkos_view("xyz", xyz));
             })
         // A0
         .def_property("A0",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.A0);
             },
-            [] (MACEKokkos& self, py::array_t<double> A0) {
+            [] (MACEKokkos<Precision>& self, py::array_t<Precision> A0) {
                 const int num_nodes = A0.size()/(self.num_lm*self.num_channels);
                 set_kokkos_view(self.A0, A0, num_nodes, self.num_lm, self.num_channels);
             })
         .def_property("A0_adj",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.A0_adj);
             },
-            [] (MACEKokkos& self, py::array_t<double> A0_adj) {
+            [] (MACEKokkos<Precision>& self, py::array_t<Precision> A0_adj) {
                 const int num_nodes = A0_adj.size()/(self.num_lm*self.num_channels);
                 set_kokkos_view(self.A0_adj, A0_adj, num_nodes, self.num_lm, self.num_channels);
             })
         .def("compute_A0",
-            [] (MACEKokkos& self,
+            [] (MACEKokkos<Precision>& self,
                     int num_nodes,
                     py::array_t<int> node_types,
                     py::array_t<int> num_neigh,
@@ -135,7 +137,7 @@ void bind_mace_kokkos(py::module_ &m)
                     create_kokkos_view("neigh_types", neigh_types));
             })
         .def("reverse_A0",
-            [] (MACEKokkos& self,
+            [] (MACEKokkos<Precision>& self,
                     int num_nodes,
                     py::array_t<int> node_types,
                     py::array_t<int> num_neigh,
@@ -151,7 +153,7 @@ void bind_mace_kokkos(py::module_ &m)
                     create_kokkos_view("r", r));
             })
         .def("compute_A0_scaled",
-            [](MACEKokkos& self,
+            [](MACEKokkos<Precision>& self,
                     const int num_nodes,
                     py::array_t<int> node_types,
                     py::array_t<int> num_neigh,
@@ -165,7 +167,7 @@ void bind_mace_kokkos(py::module_ &m)
                     create_kokkos_view("r", r));
             })
         .def("reverse_A0_scaled",
-            [](MACEKokkos& self,
+            [](MACEKokkos<Precision>& self,
                     const int num_nodes,
                     py::array_t<int> node_types,
                     py::array_t<int> num_neigh,
@@ -182,23 +184,23 @@ void bind_mace_kokkos(py::module_ &m)
             })
         // M0
         .def_property("M0",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.M0);
             },
-            [] (MACEKokkos& self, py::array_t<double> M0) {
+            [] (MACEKokkos<Precision>& self, py::array_t<Precision> M0) {
                 const int num_nodes = M0.size()/(self.num_LM*self.num_channels);
                 set_kokkos_view(self.M0, M0, num_nodes, self.num_LM, self.num_channels);
             })
         .def_property("M0_adj",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.M0_adj);
             },
-            [] (MACEKokkos& self, py::array_t<double> M0_adj) {
+            [] (MACEKokkos<Precision>& self, py::array_t<Precision> M0_adj) {
                 const int num_nodes = M0_adj.size()/(self.num_LM*self.num_channels);
                 set_kokkos_view(self.M0_adj, M0_adj, num_nodes, self.num_LM, self.num_channels);
             })
         .def("compute_M0",
-            [] (MACEKokkos& self,
+            [] (MACEKokkos<Precision>& self,
                     int num_nodes,
                     py::array_t<int> node_types) {
                 self.compute_M0(
@@ -206,7 +208,7 @@ void bind_mace_kokkos(py::module_ &m)
                     create_kokkos_view("node_types", node_types));
             })
         .def("reverse_M0",
-            [] (MACEKokkos& self,
+            [] (MACEKokkos<Precision>& self,
                     int num_nodes,
                     py::array_t<int> node_types) {
                 self.reverse_M0(
@@ -215,42 +217,42 @@ void bind_mace_kokkos(py::module_ &m)
             })
         // H1
         .def_property("H1",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.H1);
             },
-            [] (MACEKokkos& self, py::array_t<double> H1) {
+            [] (MACEKokkos<Precision>& self, py::array_t<Precision> H1) {
                 const int num_nodes = H1.size()/(self.num_LM*self.num_channels);
                 set_kokkos_view(self.H1, H1, num_nodes, self.num_LM, self.num_channels);
             })
         .def_property("H1_adj",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.H1_adj);
             },
-            [] (MACEKokkos& self, py::array_t<double> H1_adj) {
+            [] (MACEKokkos<Precision>& self, py::array_t<Precision> H1_adj) {
                 const int num_nodes = H1_adj.size()/(self.num_LM*self.num_channels);
                 set_kokkos_view(self.H1_adj, H1_adj, num_nodes, self.num_LM, self.num_channels);
             })
-        .def("compute_H1", &MACEKokkos::compute_H1)
-        .def("reverse_H1", &MACEKokkos::reverse_H1)
+        .def("compute_H1", &MACEKokkos<Precision>::compute_H1)
+        .def("reverse_H1", &MACEKokkos<Precision>::reverse_H1)
         // Phi1
         .def_property("Phi1",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.Phi1);
             },
-            [] (MACEKokkos& self, py::array_t<double> Phi1) {
+            [] (MACEKokkos<Precision>& self, py::array_t<Precision> Phi1) {
                 const int num_nodes = Phi1.size()/(self.num_lme*self.num_channels);
                 set_kokkos_view(self.Phi1, Phi1, num_nodes, self.num_lme, self.num_channels);
             })
         .def_property("Phi1_adj",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.dPhi1);
             },
-            [] (MACEKokkos& self, py::array_t<double> dPhi1) {
+            [] (MACEKokkos<Precision>& self, py::array_t<Precision> dPhi1) {
                 const int num_nodes = dPhi1.size()/(self.num_lme*self.num_channels);
                 set_kokkos_view(self.dPhi1, dPhi1, num_nodes, self.num_lme, self.num_channels);
             })
         .def("compute_Phi1",
-            [] (MACEKokkos& self,
+            [] (MACEKokkos<Precision>& self,
                     const int num_nodes,
                     py::array_t<int> num_neigh,
                     py::array_t<int> neigh_types) {
@@ -260,7 +262,7 @@ void bind_mace_kokkos(py::module_ &m)
                     create_kokkos_view("neigh_types", neigh_types));
             })
         .def("reverse_Phi1",
-            [] (MACEKokkos& self,
+            [] (MACEKokkos<Precision>& self,
                     const int num_nodes,
                     py::array_t<int> num_neigh,
                     py::array_t<int> neigh_indices,
@@ -279,31 +281,31 @@ void bind_mace_kokkos(py::module_ &m)
             })
         // A1
         .def_property("A1",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.A1);
             },
-            [] (MACEKokkos& self, py::array_t<double> A1) {
+            [] (MACEKokkos<Precision>& self, py::array_t<Precision> A1) {
                 const int num_nodes = A1.size()/(self.num_lm*self.num_channels);
                 set_kokkos_view(self.A1, A1, num_nodes, self.num_lm, self.num_channels);
             })
         .def_property("A1_adj",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.A1_adj);
             },
-            [] (MACEKokkos& self, py::array_t<double> A1_adj) {
+            [] (MACEKokkos<Precision>& self, py::array_t<Precision> A1_adj) {
                 const int num_nodes = A1_adj.size()/(self.num_lm*self.num_channels);
                 set_kokkos_view(self.A1_adj, A1_adj, num_nodes, self.num_lm, self.num_channels);
             })
         .def("compute_A1",
-            [] (MACEKokkos& self, int num_nodes) {
+            [] (MACEKokkos<Precision>& self, int num_nodes) {
                 self.compute_A1(num_nodes);
             })
         .def("reverse_A1",
-            [] (MACEKokkos& self, int num_nodes) {
+            [] (MACEKokkos<Precision>& self, int num_nodes) {
                 self.reverse_A1(num_nodes);
             })
         .def("compute_A1_scaled",
-            [](MACEKokkos& self,
+            [](MACEKokkos<Precision>& self,
                     const int num_nodes,
                     py::array_t<int> node_types,
                     py::array_t<int> num_neigh,
@@ -317,7 +319,7 @@ void bind_mace_kokkos(py::module_ &m)
                     create_kokkos_view("r", r));
             })
         .def("reverse_A1_scaled",
-            [](MACEKokkos& self,
+            [](MACEKokkos<Precision>& self,
                     const int num_nodes,
                     py::array_t<int> node_types,
                     py::array_t<int> num_neigh,
@@ -334,23 +336,23 @@ void bind_mace_kokkos(py::module_ &m)
             })
         // M1
         .def_property("M1",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.M1);
             },
-            [] (MACEKokkos& self, py::array_t<double> M1) {
+            [] (MACEKokkos<Precision>& self, py::array_t<Precision> M1) {
                 const int num_nodes = M1.size()/(self.num_channels);
                 set_kokkos_view(self.M1, M1, num_nodes, self.num_channels);
             })
         .def_property("M1_adj",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.M1_adj);
             },
-            [] (MACEKokkos& self, py::array_t<double> M1_adj) {
+            [] (MACEKokkos<Precision>& self, py::array_t<Precision> M1_adj) {
                 const int num_nodes = M1_adj.size()/(self.num_channels);
                 set_kokkos_view(self.M1_adj, M1_adj, num_nodes, self.num_channels);
             })
         .def("compute_M1",
-            [] (MACEKokkos& self,
+            [] (MACEKokkos<Precision>& self,
                     int num_nodes,
                     py::array_t<int> node_types) {
                 self.compute_M1(
@@ -358,7 +360,7 @@ void bind_mace_kokkos(py::module_ &m)
                     create_kokkos_view("node_types", node_types));
             })
         .def("reverse_M1",
-            [] (MACEKokkos& self,
+            [] (MACEKokkos<Precision>& self,
                     int num_nodes,
                     py::array_t<int> node_types) {
                 self.reverse_M1(
@@ -367,23 +369,23 @@ void bind_mace_kokkos(py::module_ &m)
             })
         // H2
         .def_property("H2",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.H2);
             },
-            [] (MACEKokkos& self, py::array_t<double> H2) {
+            [] (MACEKokkos<Precision>& self, py::array_t<double> H2) {
                 const int num_nodes = H2.size()/self.num_channels;
                 set_kokkos_view(self.H2, H2, num_nodes, self.num_channels);
             })
         .def_property("H2_adj",
-            [] (MACEKokkos& self) {
+            [] (MACEKokkos<Precision>& self) {
                 return view2vector(self.H2_adj);
             },
-            [] (MACEKokkos& self, py::array_t<double> H2_adj) {
+            [] (MACEKokkos<Precision>& self, py::array_t<double> H2_adj) {
                 const int num_nodes = H2_adj.size()/self.num_channels;
                 set_kokkos_view(self.H2_adj, H2_adj, num_nodes, self.num_channels);
             })
         .def("compute_H2",
-            [] (MACEKokkos& self,
+            [] (MACEKokkos<Precision>& self,
                     const int num_nodes,
                     py::array_t<int> node_types) {
                 self.compute_H2(
@@ -391,7 +393,7 @@ void bind_mace_kokkos(py::module_ &m)
                     create_kokkos_view("node_types", node_types));
             })
         .def("reverse_H2",
-            [] (MACEKokkos& self,
+            [] (MACEKokkos<Precision>& self,
                     const int num_nodes,
                     py::array_t<int> node_types,
                     bool zero_H1_adj) {
@@ -402,11 +404,17 @@ void bind_mace_kokkos(py::module_ &m)
             })
         // readouts
         .def("compute_readouts",
-            [] (MACEKokkos& self,
+            [] (MACEKokkos<Precision>& self,
                     const int num_nodes,
                     py::array_t<int> node_types) {
-                self.compute_readouts(
+                return self.compute_readouts(
                     num_nodes,
                     create_kokkos_view("node_types", node_types));
             });
 }
+
+void bind_mace_kokkos(py::module_ &m)
+{
+    bind_mace_kokkos<double>(m, "MACEKokkos");
+    bind_mace_kokkos<float>(m, "MACEKokkosFloat");
+}
diff --git a/symmetrix/source/cpp/utilities_kokkos.hpp b/symmetrix/source/cpp/utilities_kokkos.hpp
index dc04727..ada4de5 100644
--- a/symmetrix/source/cpp/utilities_kokkos.hpp
+++ b/symmetrix/source/cpp/utilities_kokkos.hpp
@@ -29,7 +29,7 @@ void set_kokkos_view(
 {
     auto d_array = Kokkos::create_mirror_view_and_copy(
         Kokkos::DefaultExecutionSpace(),
-        Kokkos::View<const double*,
+        Kokkos::View<const T*,
                      Kokkos::HostSpace,
                      Kokkos::MemoryUnmanaged>(array.data(), array.size()));
     if (view.size() != array.size())
@@ -46,7 +46,7 @@ void set_kokkos_view(
 {
     auto d_array = Kokkos::create_mirror_view_and_copy(
         Kokkos::DefaultExecutionSpace(),
-        Kokkos::View<const double**,
+        Kokkos::View<const T**,
                      Kokkos::LayoutRight,
                      Kokkos::HostSpace,
                      Kokkos::MemoryUnmanaged>(array.data(), N0, N1));
@@ -65,7 +65,7 @@ void set_kokkos_view(
 {
     auto d_array = Kokkos::create_mirror_view_and_copy(
         Kokkos::DefaultExecutionSpace(),
-        Kokkos::View<const double***,
+        Kokkos::View<const T***,
                      Kokkos::LayoutRight,
                      Kokkos::HostSpace,
                      Kokkos::MemoryUnmanaged>(array.data(), N0, N1, N2));
diff --git a/symmetrix/source/symmetrix/symmetrix_calc.py b/symmetrix/source/symmetrix/symmetrix_calc.py
index 6eaec50..067af9a 100755
--- a/symmetrix/source/symmetrix/symmetrix_calc.py
+++ b/symmetrix/source/symmetrix/symmetrix_calc.py
@@ -18,7 +18,7 @@
 from ase.calculators.calculator import Calculator, PropertyNotImplementedError, all_changes
 from ase.stress import full_3x3_to_voigt_6_stress
 
-from .symmetrix import MACE
+from . import symmetrix
 
 class Symmetrix(Calculator):
     """ASE Calculator using symmetrix library to evaluate equivariant graph neural network 
@@ -36,8 +36,21 @@ class Symmetrix(Calculator):
     implemented_properties = ['energy', 'free_energy', 'energies', 'forces', 'stress']
 
 
-    def __init__(self, model_file, **kwargs):
+    def __init__(self, model_file, dtype="float64", use_kokkos=True, **kwargs):
         Calculator.__init__(self, **kwargs)
+        if dtype not in ["float32", "float64"]:
+            raise ValueError(f"Unsupported dtype '{dtype}'. Supported dtypes are 'float64' and 'float32'.")
+        if use_kokkos and not hasattr(symmetrix, "MACEKokkos"):
+            raise RuntimeError("Symmetrix was built without Kokkos support.")
+        self.use_kokkos = use_kokkos
+        if self.use_kokkos:
+            if not symmetrix._kokkos_is_initialized():
+                symmetrix._init_kokkos()
+            MACE = symmetrix.MACEKokkos if dtype == "float64" else symmetrix.MACEKokkosFloat
+        else:
+            if dtype == "float32":
+                raise ValueError(f"dtype '{dtype}' requires `use_kokkos = True`")
+            MACE = symmetrix.MACE
         try:
             self.evaluator = MACE(str(model_file))
         except RuntimeError: # expecting json.exception.parse_error.101
@@ -56,7 +69,6 @@ def __init__(self, model_file, **kwargs):
 
         self.cutoff = self.evaluator.r_cut
 
-
     def calculate(self, atoms=None, properties=['energy'], system_changes=all_changes):
         Calculator.calculate(self, atoms, properties, system_changes)
 
@@ -74,6 +86,8 @@ def calculate(self, atoms=None, properties=['energy'], system_changes=all_change
         self.results['energies'] = np.asarray(self.evaluator.node_energies)
 
         pair_forces = np.asarray(self.evaluator.node_forces).reshape((-1, 3))
+        pair_forces = pair_forces[:len(i_list), :]  # currently, `evaluator.node_forces` is a container
+                                                    # which can grow larger than the actual number of pairs
 
         # atom forces from pair_forces
         N_atoms = len(self.atoms)
diff --git a/symmetrix/test/test_symmetrix_calc.py b/symmetrix/test/test_symmetrix_calc.py
index 15af4bc..62330c0 100755
--- a/symmetrix/test/test_symmetrix_calc.py
+++ b/symmetrix/test/test_symmetrix_calc.py
@@ -12,8 +12,6 @@
 from ase.atoms import Atoms
 from ase.stress import full_3x3_to_voigt_6_stress
 
-# TODO: currently skipping these tests - revisit
-
 try:
     from symmetrix import Symmetrix
 except ModuleNotFoundError as exc:
@@ -51,13 +49,14 @@ def mace_foundation_model(tmp_path_factory):
         return str(downloaded_model)
 
 
-def test_calc_caching(model_cache):
+@pytest.mark.parametrize("use_kokkos", [True, False])
+def test_calc_caching(model_cache, use_kokkos):
     atoms = Atoms('O', cell=[2] * 3, pbc=[True] * 3)
     atoms *= 4
     rng = np.random.default_rng(5)
     atoms.rattle(rng=rng)
 
-    calc = Symmetrix(model_cache["mace-mp-0b3-medium-1-8.json"])
+    calc = Symmetrix(model_cache["mace-mp-0b3-medium-1-8.json"], use_kokkos=use_kokkos)
     atoms.calc = calc
 
     t0 = time.time()
@@ -81,7 +80,8 @@ def test_calc_caching(model_cache):
     assert np.abs(dt_F_pert - dt_E) / dt_E < 0.5
 
 
-def test_symmetrix_calc_finite_diff(model_cache):
+@pytest.mark.parametrize("use_kokkos", [True, False])
+def test_symmetrix_calc_finite_diff(model_cache, use_kokkos):
     atoms = Atoms('O', cell=[2] * 3, pbc=[True] * 3)
     atoms *= 2
     rng = np.random.default_rng(5)
@@ -91,12 +91,13 @@ def test_symmetrix_calc_finite_diff(model_cache):
     atoms.set_cell(atoms.cell @ F, True)
 
     print("pre-converted")
-    calc = Symmetrix(model_cache["mace-mp-0b3-medium-1-8.json"])
+    calc = Symmetrix(model_cache["mace-mp-0b3-medium-1-8.json"], use_kokkos=use_kokkos)
     do_grad_test(atoms, calc, True)
 
 
 @pytest.mark.skipif(mace is None, reason="mace-torch is not available")
-def test_mace_onthefly_calc_finite_diff(mace_foundation_model):
+@pytest.mark.parametrize("use_kokkos", [True, False])
+def test_mace_onthefly_calc_finite_diff(mace_foundation_model, use_kokkos):
     atoms = Atoms('O', cell=[2] * 3, pbc=[True] * 3)
     atoms *= 2
     rng = np.random.default_rng(5)
@@ -106,12 +107,13 @@ def test_mace_onthefly_calc_finite_diff(mace_foundation_model):
     atoms.set_cell(atoms.cell @ F, True)
 
     print("converted on-the-fly")
-    calc = Symmetrix(mace_foundation_model, species=[1, 8])
+    calc = Symmetrix(mace_foundation_model, species=[1, 8], use_kokkos=use_kokkos)
     do_grad_test(atoms, calc, True)
 
 
 @pytest.mark.skipif(mace is None, reason="mace-torch is not available")
-def test_symmetrix_vs_pytorch(mace_foundation_model):
+@pytest.mark.parametrize("use_kokkos", [True, False])
+def test_symmetrix_vs_pytorch(mace_foundation_model, use_kokkos):
     atoms = Atoms('O', cell=[2] * 3, pbc=[True] * 3)
     atoms *= 2
     rng = np.random.default_rng(5)
@@ -123,7 +125,7 @@ def test_symmetrix_vs_pytorch(mace_foundation_model):
     atoms_s = atoms.copy()
     atoms_p = atoms.copy()
 
-    calc_sym = Symmetrix(mace_foundation_model, species=[1, 8])
+    calc_sym = Symmetrix(mace_foundation_model, species=[1, 8], use_kokkos=use_kokkos)
     atoms_s.calc = calc_sym
 
     calc_torch = MACECalculator(mace_foundation_model)