diff --git a/benchmarks/impl_point.hpp b/benchmarks/impl_point.hpp
index 5dc1d932..5150c438 100644
--- a/benchmarks/impl_point.hpp
+++ b/benchmarks/impl_point.hpp
@@ -6,7 +6,7 @@
 namespace specfem {
 namespace benchmarks {
 
-#define PROPERTIES_CONSTRUCTOR                                                 \
+#define DEFINE_PROP_CONSTRUCTORS                                               \
   KOKKOS_FUNCTION                                                              \
   properties() = default;                                                      \
   KOKKOS_FUNCTION                                                              \
@@ -18,6 +18,10 @@ namespace benchmarks {
     compute();                                                                 \
   }
 
+#define DEFINE_PROP(prop, i)                                                   \
+  KOKKOS_INLINE_FUNCTION value_type prop() const { return Base::data[i]; }     \
+  KOKKOS_INLINE_FUNCTION void prop(value_type val) { Base::data[i] = val; }
+
 namespace impl {
 template <int N, int NALL, bool UseSIMD> struct impl_properties {
   using simd = specfem::datatype::simd<type_real, UseSIMD>; ///< SIMD type
@@ -114,31 +118,21 @@ struct properties<specfem::dimension::type::dim2,
       typename simd::datatype; ///< Value type to store properties
   ///@}
 
-  value_type &lambdaplus2mu =
-      Base::data[0];               ///< Lame's parameter @f$ \lambda + 2\mu @f$
-  value_type &mu = Base::data[1];  ///< shear modulus @f$ \mu @f$
-  value_type &rho = Base::data[2]; ///< density @f$ \rho @f$
-
-  value_type &rho_vp = Base::data[3]; ///< P-wave velocity @f$ \rho v_p @f$
-  value_type &rho_vs = Base::data[4]; ///< S-wave velocity @f$ \rho v_s @f$
-  value_type &lambda = Base::data[5]; ///< Lame's parameter @f$ \lambda @f$
+  DEFINE_PROP(lambdaplus2mu, 0) ///< Lame's parameter @f$ \lambda + 2\mu @f$
+  DEFINE_PROP(mu, 1)            ///< shear modulus @f$ \mu @f$
+  DEFINE_PROP(rho, 2)           ///< density @f$ \rho @f$
 
-  PROPERTIES_CONSTRUCTOR
+  DEFINE_PROP(rho_vp, 3) ///< P-wave velocity @f$ \rho v_p @f$
+  DEFINE_PROP(rho_vs, 4) ///< S-wave velocity @f$ \rho v_s @f$
+  DEFINE_PROP(lambda, 5) ///< Lame's parameter @f$ \lambda @f$
 
-  // KOKKOS_FUNCTION
-  // properties(const value_type &lambdaplus2mu, const value_type &mu,
-  //            const value_type &rho) {
-  //   this->lambdaplus2mu = lambdaplus2mu;
-  //   this->mu = mu;
-  //   this->rho = rho;
-  //   compute();
-  // }
+  DEFINE_PROP_CONSTRUCTORS
 
   KOKKOS_INLINE_FUNCTION
   void compute() {
-    rho_vp = Kokkos::sqrt(rho * lambdaplus2mu);
-    rho_vs = Kokkos::sqrt(rho * mu);
-    lambda = lambdaplus2mu - (static_cast<value_type>(2.0)) * mu;
+    rho_vp(Kokkos::sqrt(rho() * lambdaplus2mu()));
+    rho_vs(Kokkos::sqrt(rho() * mu()));
+    lambda(lambdaplus2mu() - (static_cast<value_type>(2.0)) * mu());
   }
 };
 
@@ -169,46 +163,27 @@ struct properties<specfem::dimension::type::dim2,
    *
    */
   ///@{
-  value_type &c11 = Base::data[0]; ///< @f$ c_{11} @f$
-  value_type &c13 = Base::data[1]; ///< @f$ c_{13} @f$
-  value_type &c15 = Base::data[2]; ///< @f$ c_{15} @f$
-  value_type &c33 = Base::data[3]; ///< @f$ c_{33} @f$
-  value_type &c35 = Base::data[4]; ///< @f$ c_{35} @f$
-  value_type &c55 = Base::data[5]; ///< @f$ c_{55} @f$
-  value_type &c12 = Base::data[6]; ///< @f$ c_{12} @f$
-  value_type &c23 = Base::data[7]; ///< @f$ c_{23} @f$
-  value_type &c25 = Base::data[8]; ///< @f$ c_{25} @f$
+  DEFINE_PROP(c11, 0); ///< @f$ c_{11} @f$
+  DEFINE_PROP(c13, 1); ///< @f$ c_{13} @f$
+  DEFINE_PROP(c15, 2); ///< @f$ c_{15} @f$
+  DEFINE_PROP(c33, 3); ///< @f$ c_{33} @f$
+  DEFINE_PROP(c35, 4); ///< @f$ c_{35} @f$
+  DEFINE_PROP(c55, 5); ///< @f$ c_{55} @f$
+  DEFINE_PROP(c12, 6); ///< @f$ c_{12} @f$
+  DEFINE_PROP(c23, 7); ///< @f$ c_{23} @f$
+  DEFINE_PROP(c25, 8); ///< @f$ c_{25} @f$
   ///@}
 
-  value_type &rho = Base::data[9];     ///< Density @f$ \rho @f$
-  value_type &rho_vp = Base::data[10]; ///< P-wave velocity @f$ \rho v_p @f$
-  value_type &rho_vs = Base::data[11]; ///< S-wave velocity @f$ \rho v_s @f$
-
-  PROPERTIES_CONSTRUCTOR
-
-  // KOKKOS_FUNCTION
-  // properties(const value_type &c11, const value_type &c13,
-  //            const value_type &c15, const value_type &c33,
-  //            const value_type &c35, const value_type &c55,
-  //            const value_type &c12, const value_type &c23,
-  //            const value_type &c25, const value_type &rho) {
-  //   this->c11 = c11;
-  //   this->c13 = c13;
-  //   this->c15 = c15;
-  //   this->c33 = c33;
-  //   this->c35 = c35;
-  //   this->c55 = c55;
-  //   this->c12 = c12;
-  //   this->c23 = c23;
-  //   this->c25 = c25;
-  //   this->rho = rho;
-  //   compute();
-  // }
+  DEFINE_PROP(rho, 9);     ///< Density @f$ \rho @f$
+  DEFINE_PROP(rho_vp, 10); ///< P-wave velocity @f$ \rho v_p @f$
+  DEFINE_PROP(rho_vs, 11); ///< S-wave velocity @f$ \rho v_s @f$
+
+  DEFINE_PROP_CONSTRUCTORS
 
   KOKKOS_INLINE_FUNCTION
   void compute() {
-    rho_vp = Kokkos::sqrt(rho * c33);
-    rho_vs = Kokkos::sqrt(rho * c55);
+    rho_vp(Kokkos::sqrt(rho() * c33()));
+    rho_vs(Kokkos::sqrt(rho() * c55()));
   }
 };
 
@@ -240,30 +215,23 @@ struct properties<specfem::dimension::type::dim2,
       typename simd::datatype; ///< Value type to store properties
   ///@}
 
-  value_type &rho_inverse = Base::data[0]; ///< @f$ \frac{1}{\rho} @f$
-  value_type &kappa = Base::data[1];       ///< Bulk modulus @f$ \kappa @f$
-
-  value_type &kappa_inverse =
-      Base::data[2]; ///< @f$ \frac{1}{\lambda + 2\mu} @f$
-  value_type &rho_vpinverse = Base::data[3]; ///< @f$ \frac{1}{\rho v_p} @f$
-
-  PROPERTIES_CONSTRUCTOR
+  DEFINE_PROP_CONSTRUCTORS
 
-  // KOKKOS_FUNCTION
-  // properties(const value_type &rho_inverse, const value_type &kappa) {
-  //   this->rho_inverse = rho_inverse;
-  //   this->kappa = kappa;
-  //   compute();
-  // }
+  DEFINE_PROP(rho_inverse, 0)   ///< @f$ \frac{1}{\rho} @f$
+  DEFINE_PROP(kappa, 1)         ///< Bulk modulus @f$ \kappa @f$
+  DEFINE_PROP(kappa_inverse, 2) ///< @f$ \frac{1}{\lambda + 2\mu} @f$
+  DEFINE_PROP(rho_vpinverse, 3) ///< @f$ \frac{1}{\rho v_p} @f$
 
+private:
   KOKKOS_INLINE_FUNCTION
   void compute() {
-    kappa_inverse = (static_cast<value_type>(1.0)) / kappa;
-    rho_vpinverse = Kokkos::sqrt(rho_inverse * kappa_inverse);
+    kappa_inverse((static_cast<value_type>(1.0)) / kappa());
+    rho_vpinverse(Kokkos::sqrt(rho_inverse() * kappa_inverse()));
   }
 };
 
-#undef PROPERTIES_CONSTRUCTOR
+#undef DEFINE_PROP_CONSTRUCTORS
+#undef DEFINE_PROP
 
 } // namespace benchmarks
 } // namespace specfem
diff --git a/benchmarks/main.cpp b/benchmarks/main.cpp
index 060acfae..e00c8792 100644
--- a/benchmarks/main.cpp
+++ b/benchmarks/main.cpp
@@ -1,4 +1,4 @@
-// #include "archive/stiffness1.hpp"
+#include "archive/stiffness1.hpp"
 #include "execute.hpp"
 #include "stiffness.hpp"
 // #include "divide.hpp"
@@ -47,26 +47,30 @@ void benchmark(specfem::compute::assembly &assembly,
   // field.elastic.nglob, field.elastic.components);
 
   for (const auto [istep, dt] : time_scheme->iterate_forward()) {
-    compute_stiffness_interaction<acoustic, isotropic, flag>(assembly, istep);
-    compute_stiffness_interaction<elastic, isotropic, flag>(assembly, istep);
-    compute_stiffness_interaction<elastic, anisotropic, flag>(assembly, istep);
+    // compute_stiffness_interaction<acoustic, isotropic, flag>(assembly,
+    // istep); compute_stiffness_interaction<elastic, isotropic, flag>(assembly,
+    // istep); compute_stiffness_interaction<elastic, anisotropic,
+    // flag>(assembly, istep);
 
     // compute_stiffness_interaction2<acoustic, isotropic, flag>(assembly,
     // istep); compute_stiffness_interaction2<elastic, isotropic,
     // flag>(assembly, istep); compute_stiffness_interaction2<elastic,
     // anisotropic, flag>(assembly, istep);
 
-    // if constexpr (flag) {
-    //   compute_stiffness_interaction<acoustic, isotropic, false>(assembly,
-    //   istep); compute_stiffness_interaction<elastic, isotropic,
-    //   false>(assembly, istep); compute_stiffness_interaction<elastic,
-    //   anisotropic, false>(assembly, istep);
-    // } else {
-    //   compute_stiffness_interaction2<acoustic, isotropic, false>(assembly,
-    //   istep); compute_stiffness_interaction2<elastic, isotropic,
-    //   false>(assembly, istep); compute_stiffness_interaction2<elastic,
-    //   anisotropic, false>(assembly, istep);
-    // }
+    if constexpr (flag) {
+      compute_stiffness_interaction<acoustic, isotropic, false>(assembly,
+                                                                istep);
+      compute_stiffness_interaction<elastic, isotropic, false>(assembly, istep);
+      compute_stiffness_interaction<elastic, anisotropic, false>(assembly,
+                                                                 istep);
+    } else {
+      compute_stiffness_interaction2<acoustic, isotropic, false>(assembly,
+                                                                 istep);
+      compute_stiffness_interaction2<elastic, isotropic, false>(assembly,
+                                                                istep);
+      compute_stiffness_interaction2<elastic, anisotropic, false>(assembly,
+                                                                  istep);
+    }
 
     // divide_mass_matrix<dimension, wavefield, acoustic>(assembly);
     // divide_mass_matrix<dimension, wavefield, elastic>(assembly);
diff --git a/benchmarks/stiffness.hpp b/benchmarks/stiffness.hpp
index 81d70d79..045196af 100644
--- a/benchmarks/stiffness.hpp
+++ b/benchmarks/stiffness.hpp
@@ -187,15 +187,9 @@ void compute_stiffness_interaction(const specfem::compute::assembly &assembly,
                       xiz * df_dxi[icomponent] + gammaz * df_dgamma[icomponent];
                 }
 
-                auto point_property = [&]() {
-                  if constexpr (flag) {
-                    return specfem::point::properties<
-                        dimension, medium_tag, property_tag, using_simd>();
-                  } else {
-                    return specfem::benchmarks::properties<
-                        dimension, medium_tag, property_tag, using_simd>();
-                  }
-                }();
+                specfem::benchmarks::properties<dimension, medium_tag,
+                                                property_tag, using_simd>
+                    point_property;
 
                 specfem::point::field_derivatives<dimension, medium_tag,
                                                   using_simd>
@@ -211,14 +205,14 @@ void compute_stiffness_interaction(const specfem::compute::assembly &assembly,
                                     specfem::element::medium_tag::acoustic &&
                                 PropertyTag ==
                                     specfem::element::property_tag::isotropic) {
-                    point_property.rho_inverse =
-                        container.rho_inverse(ispec, iz, ix);
-                    point_property.kappa = container.kappa(ispec, iz, ix);
-                    point_property.kappa_inverse =
-                        static_cast<type_real>(1.0) / point_property.kappa;
-                    point_property.rho_vpinverse =
-                        sqrt(point_property.rho_inverse *
-                             point_property.kappa_inverse);
+                    point_property.rho_inverse(
+                        container.rho_inverse(ispec, iz, ix));
+                    point_property.kappa(container.kappa(ispec, iz, ix));
+                    point_property.kappa_inverse(static_cast<type_real>(1.0) /
+                                                 point_property.kappa());
+                    point_property.rho_vpinverse(
+                        sqrt(point_property.rho_inverse() *
+                             point_property.kappa_inverse()));
 
                     const auto &du = field_derivatives.du;
 
@@ -226,8 +220,8 @@ void compute_stiffness_interaction(const specfem::compute::assembly &assembly,
                                                            using_simd>
                         T;
 
-                    T(0, 0) = point_property.rho_inverse * du(0, 0);
-                    T(1, 0) = point_property.rho_inverse * du(1, 0);
+                    T(0, 0) = point_property.rho_inverse() * du(0, 0);
+                    T(1, 0) = point_property.rho_inverse() * du(1, 0);
 
                     return specfem::point::stress<
                         specfem::dimension::type::dim2,
@@ -237,16 +231,16 @@ void compute_stiffness_interaction(const specfem::compute::assembly &assembly,
                       MediumTag == specfem::element::medium_tag::elastic &&
                       PropertyTag ==
                           specfem::element::property_tag::isotropic) {
-                    point_property.rho = container.rho(ispec, iz, ix);
-                    point_property.mu = container.mu(ispec, iz, ix);
-                    point_property.lambdaplus2mu =
-                        container.lambdaplus2mu(ispec, iz, ix);
-                    point_property.lambda =
-                        point_property.lambdaplus2mu - 2 * point_property.mu;
-                    point_property.rho_vp =
-                        sqrt(point_property.rho * point_property.lambdaplus2mu);
-                    point_property.rho_vs =
-                        sqrt(point_property.rho * point_property.mu);
+                    point_property.rho(container.rho(ispec, iz, ix));
+                    point_property.mu(container.mu(ispec, iz, ix));
+                    point_property.lambdaplus2mu(
+                        container.lambdaplus2mu(ispec, iz, ix));
+                    point_property.lambda(point_property.lambdaplus2mu() -
+                                          2 * point_property.mu());
+                    point_property.rho_vp(sqrt(point_property.rho() *
+                                               point_property.lambdaplus2mu()));
+                    point_property.rho_vs(
+                        sqrt(point_property.rho() * point_property.mu()));
 
                     using datatype =
                         typename specfem::datatype::simd<type_real,
@@ -255,11 +249,11 @@ void compute_stiffness_interaction(const specfem::compute::assembly &assembly,
 
                     datatype sigma_xx, sigma_zz, sigma_xz;
 
-                    sigma_xx = point_property.lambdaplus2mu * du(0, 0) +
-                               point_property.lambda * du(1, 1);
-                    sigma_zz = point_property.lambdaplus2mu * du(1, 1) +
-                               point_property.lambda * du(0, 0);
-                    sigma_xz = point_property.mu * (du(0, 1) + du(1, 0));
+                    sigma_xx = point_property.lambdaplus2mu() * du(0, 0) +
+                               point_property.lambda() * du(1, 1);
+                    sigma_zz = point_property.lambdaplus2mu() * du(1, 1) +
+                               point_property.lambda() * du(0, 0);
+                    sigma_xz = point_property.mu() * (du(0, 1) + du(1, 0));
 
                     specfem::datatype::VectorPointViewType<type_real, 2, 2,
                                                            using_simd>
@@ -277,21 +271,21 @@ void compute_stiffness_interaction(const specfem::compute::assembly &assembly,
                       MediumTag == specfem::element::medium_tag::elastic &&
                       PropertyTag ==
                           specfem::element::property_tag::anisotropic) {
-                    point_property.rho = container.rho(ispec, iz, ix);
-                    point_property.c11 = container.c11(ispec, iz, ix);
-                    point_property.c12 = container.c12(ispec, iz, ix);
-                    point_property.c13 = container.c13(ispec, iz, ix);
-                    point_property.c15 = container.c15(ispec, iz, ix);
-                    point_property.c33 = container.c33(ispec, iz, ix);
-                    point_property.c35 = container.c35(ispec, iz, ix);
-                    point_property.c55 = container.c55(ispec, iz, ix);
-                    point_property.c23 = container.c23(ispec, iz, ix);
-                    point_property.c25 = container.c25(ispec, iz, ix);
-
-                    point_property.rho_vp =
-                        sqrt(point_property.rho * point_property.c33);
-                    point_property.rho_vs =
-                        sqrt(point_property.rho * point_property.c55);
+                    point_property.rho(container.rho(ispec, iz, ix));
+                    point_property.c11(container.c11(ispec, iz, ix));
+                    point_property.c12(container.c12(ispec, iz, ix));
+                    point_property.c13(container.c13(ispec, iz, ix));
+                    point_property.c15(container.c15(ispec, iz, ix));
+                    point_property.c33(container.c33(ispec, iz, ix));
+                    point_property.c35(container.c35(ispec, iz, ix));
+                    point_property.c55(container.c55(ispec, iz, ix));
+                    point_property.c23(container.c23(ispec, iz, ix));
+                    point_property.c25(container.c25(ispec, iz, ix));
+
+                    point_property.rho_vp(
+                        sqrt(point_property.rho() * point_property.c33()));
+                    point_property.rho_vs(
+                        sqrt(point_property.rho() * point_property.c55()));
 
                     using datatype =
                         typename specfem::datatype::simd<type_real,
@@ -300,17 +294,17 @@ void compute_stiffness_interaction(const specfem::compute::assembly &assembly,
 
                     datatype sigma_xx, sigma_zz, sigma_xz;
 
-                    sigma_xx = point_property.c11 * du(0, 0) +
-                               point_property.c13 * du(1, 1) +
-                               point_property.c15 * (du(1, 0) + du(0, 1));
+                    sigma_xx = point_property.c11() * du(0, 0) +
+                               point_property.c13() * du(1, 1) +
+                               point_property.c15() * (du(1, 0) + du(0, 1));
 
-                    sigma_zz = point_property.c13 * du(0, 0) +
-                               point_property.c33 * du(1, 1) +
-                               point_property.c35 * (du(1, 0) + du(0, 1));
+                    sigma_zz = point_property.c13() * du(0, 0) +
+                               point_property.c33() * du(1, 1) +
+                               point_property.c35() * (du(1, 0) + du(0, 1));
 
-                    sigma_xz = point_property.c15 * du(0, 0) +
-                               point_property.c35 * du(1, 1) +
-                               point_property.c55 * (du(1, 0) + du(0, 1));
+                    sigma_xz = point_property.c15() * du(0, 0) +
+                               point_property.c35() * du(1, 1) +
+                               point_property.c55() * (du(1, 0) + du(0, 1));
 
                     specfem::datatype::VectorPointViewType<type_real, 2, 2,
                                                            using_simd>