diff --git a/bindings/bind_py_models.cc b/bindings/bind_py_models.cc
index 0425d3090..67d2e8274 100644
--- a/bindings/bind_py_models.cc
+++ b/bindings/bind_py_models.cc
@@ -82,6 +82,18 @@ namespace rascal {
             SparsePoints, i.e. the basis used by the sparse method.
             The representation of the atomic structures computed with Calculator
             should have already been computed.)");
+    kernel.def(
+        "compute_local",
+        py::overload_cast<const Calculator &, const StructureManagers &,
+                          const SparsePoints &>(
+            &SparseKernel::template compute_local<Calculator, StructureManagers,
+                                            SparsePoints>),
+        py::call_guard<py::gil_scoped_release>(),
+        R"(Compute the sparse kernel between the representation of a set of
+            atomic structures, i.e. StructureManagerCollections, and a set of
+            SparsePoints, i.e. the basis used by the sparse method.
+            The representation of the atomic structures computed with Calculator
+            should have already been computed.)");
     kernel.def("compute",
                py::overload_cast<const SparsePoints &>(
                    &SparseKernel::template compute<SparsePoints>),
@@ -193,6 +205,43 @@ namespace rascal {
           return neg_stress_global;
         },
         py::call_guard<py::gil_scoped_release>());
+
+    mod.def(
+        "compute_sparse_kernel_local_neg_stress",
+        [](const Calculator & calculator, SparseKernel & kernel,
+           ManagerCollection & managers, SparsePoints & sparse_points,
+           math::Vector_t & weights) {
+          std::string neg_stress_name = compute_sparse_kernel_local_neg_stress(
+              calculator, kernel, managers, sparse_points, weights);
+          size_t nb_centers_over_managers{0};
+          for (const auto & manager : managers) {
+            
+            nb_centers_over_managers+=manager->size();
+            //for (auto center : manager) {
+            //  nb_centers_over_managers++;
+            //}
+          }
+          math::Matrix_t local_neg_stress{nb_centers_over_managers, 9};
+
+          size_t i_center{0};
+          size_t i_manager_nb_center{0};
+          for (const auto & manager : managers) {
+            i_manager_nb_center=manager->size();
+            // i_manager_nb_center = 0;
+            // for (auto center : manager) {
+            //   i_manager_nb_center++;
+            // }
+            auto && neg_stress{
+                *manager
+                    ->template get_property<Property<double, 1, Manager_t, 9>>(
+                        neg_stress_name, true, true, true)};
+            local_neg_stress.block(i_center, 0, i_manager_nb_center, 9) =
+                Eigen::Map<const math::Matrix_t>(neg_stress.view().data(), 1, 9);
+            i_center += i_manager_nb_center;
+          }
+          return local_neg_stress;
+        },
+        py::call_guard<py::gil_scoped_release>());
   }
 
   /**
diff --git a/bindings/rascal/lib/__init__.py b/bindings/rascal/lib/__init__.py
index 5853f1683..cf0fa5d63 100644
--- a/bindings/rascal/lib/__init__.py
+++ b/bindings/rascal/lib/__init__.py
@@ -4,4 +4,5 @@
     kernels,
     compute_sparse_kernel_gradients,
     compute_sparse_kernel_neg_stress,
+    compute_sparse_kernel_local_neg_stress
 )
diff --git a/bindings/rascal/models/krr.py b/bindings/rascal/models/krr.py
index 2ee176fd8..fd334847d 100644
--- a/bindings/rascal/models/krr.py
+++ b/bindings/rascal/models/krr.py
@@ -8,7 +8,11 @@
     train_gap_model         Train a GAP model given a kernel matrix and sparse points
 """
 from ..utils import BaseIO
-from ..lib import compute_sparse_kernel_gradients, compute_sparse_kernel_neg_stress
+from ..lib import (
+    compute_sparse_kernel_gradients,
+    compute_sparse_kernel_neg_stress,
+    compute_sparse_kernel_local_neg_stress
+)
 
 import scipy
 import numpy as np
@@ -352,7 +356,7 @@ def predict_forces(self, managers, KNM=None):
 
         return -gradients
 
-    def predict_stress(self, managers, KNM=None):
+    def predict_stress(self, managers, KNM=None, local_stress=False):
         """Predict gradients w.r.t cell parameters, e.g. stress, associated with the atomic structures in managers.
         The stress is returned using the Voigt order: xx, yy, zz, yz, xz, xy.
 
@@ -376,14 +380,26 @@ def predict_stress(self, managers, KNM=None):
 
         if KNM is None:
             rep = self.kernel._representation
-            neg_stress = compute_sparse_kernel_neg_stress(
-                rep,
-                self.kernel._kernel,
-                managers.managers,
-                self.X_train._sparse_points,
-                self.weights.reshape((1, -1)),
-            )
+            if local_stress:
+                neg_stress = compute_sparse_kernel_local_neg_stress(
+                    rep,
+                    self.kernel._kernel,
+                    managers.managers,
+                    self.X_train._sparse_points,
+                    self.weights.reshape((1, -1)),
+                )
+            else:
+                neg_stress = compute_sparse_kernel_neg_stress(
+                    rep,
+                    self.kernel._kernel,
+                    managers.managers,
+                    self.X_train._sparse_points,
+                    self.weights.reshape((1, -1)),
+                )
         else:
+            if local_stress:
+                raise ValueError("Option local_stress does not work with KNM."
+                                 " Please set KNM to None.")
             if 6 * len(managers) != KNM.shape[0]:
                 raise ValueError(
                     "KNM size mismatch {}!={}".format(6 * len(managers), KNM.shape[0])
diff --git a/src/rascal/models/sparse_kernel_predict.hh b/src/rascal/models/sparse_kernel_predict.hh
index 80337e4a4..d7f43aca8 100644
--- a/src/rascal/models/sparse_kernel_predict.hh
+++ b/src/rascal/models/sparse_kernel_predict.hh
@@ -374,5 +374,115 @@ namespace rascal {
     }  // manager
     return neg_stress_name;
   }
+
+
+  /**
+   * Compute the atomic gradients of a structure w.r.t atomic cell paramaters
+   * using a sparse GPR model. Only SOAP-GAP model is implemented at the moment
+   * (see
+   * @ref SparseKernelImpl<internal::SparseKernelType::GAP>::compute_derivative)
+   *
+   * The point of this function is to provide a faster prediction routine
+   * compared to computing the kernel elements and then multiplying them with
+   * the weights of the model.
+   *
+   * The gradients are attached to the input manager in a Property of
+   * Order 1 with the name
+   * '"kernel: "+kernel_type+" ; "+representation_grad_name+" negative stress;
+   * weight_hash:"+weight_hash' gradients [2*n_managers, 3]
+   *
+   * @tparam StructureManagers should be an iterable over shared pointer
+   *          of structure managers like ManagerCollection
+   * @param sparse_points a SparsePoints* class
+   * @param managers a ManagerCollection or similar collection of
+   * structure managers
+   * @param weights regression weights of the sparse GPR model
+   * @return name used to register the gradients in the managers
+   */
+  template <class Calculator, class StructureManagers, class SparsePoints>
+  std::string compute_sparse_kernel_local_neg_stress(const Calculator & calculator,
+                                               SparseKernel & kernel,
+                                               StructureManagers & managers,
+                                               SparsePoints & sparse_points,
+                                               math::Vector_t & weights) {
+    using Manager_t = typename StructureManagers::Manager_t;
+    using Property_t = typename Calculator::template Property_t<Manager_t>;
+    using PropertyGradient_t =
+        typename Calculator::template PropertyGradient_t<Manager_t>;
+    auto && representation_name{calculator.get_name()};
+    const auto representation_grad_name{calculator.get_gradient_name()};
+    size_t n_centers{0};
+    // find the total number of gradients
+    for (const auto & manager : managers) {
+      n_centers += manager->size();
+    }
+
+    internal::Hash<math::Vector_t, double> hasher{};
+    auto kernel_type_str = kernel.parameters.at("name").get<std::string>();
+    std::string weight_hash = std::to_string(hasher(weights));
+    std::string pair_grad_atom_i_r_j_name =
+        std::string("kernel: ") + kernel_type_str + std::string(" ; ") +
+        representation_grad_name +
+        std::string(" partial gradients; weight_hash:") + weight_hash;
+    std::string neg_stress_name =
+        std::string("kernel: ") + kernel_type_str + std::string(" ; ") +
+        representation_grad_name +
+        std::string(" negative local stress; weight_hash:") + weight_hash;
+
+    for (const auto & manager : managers) {
+      if (kernel_type_str == "GAP") {
+        const auto zeta = kernel.parameters.at("zeta").get<size_t>();
+
+        compute_partial_gradients_gap<Property_t, PropertyGradient_t>(
+            manager, sparse_points, weights, zeta, representation_name,
+            representation_grad_name, pair_grad_atom_i_r_j_name);
+      }
+
+      //auto && expansions_coefficients{*manager->template get_property<prop_t>(
+      //    this->get_name(), true, true, ExcludeGhosts)};
+      auto && neg_stress{
+          *manager->template get_property<Property<double, 1, Manager_t, 9>>(
+              neg_stress_name, true, true, true)};
+      if (neg_stress.is_updated()) {
+        continue;
+      }
+
+      neg_stress.resize();
+      neg_stress.setZero();
+
+      auto && pair_grad_atom_i_r_j{*manager->template get_property<
+          Property<double, 2, Manager_t, 1, ThreeD>>(pair_grad_atom_i_r_j_name,
+                                                     true)};
+
+      auto manager_root = extract_underlying_manager<0>(manager);
+      json structure_copy = manager_root->get_atomic_structure();
+      auto atomic_structure =
+          structure_copy.template get<AtomicStructure<ThreeD>>();
+      float volume = atomic_structure.get_volume();
+      for (auto center : manager) {
+        //auto && local_neg_stress = neg_stress[center];
+        //local_neg_stress.setZero();
+
+        Eigen::Vector3d r_i = center.get_position();
+        // accumulate partial gradients onto gradients
+        for (auto neigh : center.pairs_with_self_pair()) {
+          auto && local_neg_stress = neg_stress[neigh.get_atom_j()];
+          //local_neg_stress.setZero();
+          Eigen::Vector3d r_ji = r_i - neigh.get_position();
+          for (int a_der{0}; a_der < ThreeD; a_der++) {
+            for (int b_der{0}; b_der < ThreeD; b_der++) {
+              local_neg_stress(a_der*3 + b_der) +=
+                  r_ji(a_der) * pair_grad_atom_i_r_j[neigh](b_der)/volume;
+            }
+          }
+        }
+        // local_neg_stress /= volume;
+      }
+      neg_stress.set_updated_status(true);
+    }  // manager
+    return neg_stress_name;
+  }
+
+
 }  // namespace rascal
 #endif  // SRC_RASCAL_MODELS_SPARSE_KERNEL_PREDICT_HH_
diff --git a/src/rascal/models/sparse_kernels.hh b/src/rascal/models/sparse_kernels.hh
index 7d22896d4..a3d3f35af 100644
--- a/src/rascal/models/sparse_kernels.hh
+++ b/src/rascal/models/sparse_kernels.hh
@@ -110,6 +110,7 @@ namespace rascal {
         math::Matrix_t KNM(managers.size(), sparse_points.size());
         KNM.setZero();
         size_t ii_A{0};
+        std::cout<<"inside compute"<<std::endl;
         for (auto & manager : managers) {
           auto && propA{*manager->template get_property<Property_t>(
               representation_name, true)};
@@ -125,6 +126,92 @@ namespace rascal {
         return KNM;
       }
 
+            /**
+       * Compute the kernel between a set of structure(s) and a set of pseudo
+       * points, per atom.
+       *
+       * @tparam StructureManagers should be an iterable over shared pointer
+       *          of structure managers like ManagerCollection
+       * @param sparse_points a SparsePoints* class
+       * @param managers a ManagerCollection or similar collection of
+       * structure managers
+       * @param representation_name name under which the representation data
+       * has been registered in the elements of managers and managers_b
+       * @return kernel matrix
+       */
+      template <
+          class Property_t, internal::TargetType Type,
+          std::enable_if_t<Type == internal::TargetType::Structure, int> = 0,
+          class StructureManagers, class SparsePoints>
+      math::Matrix_t compute_local(StructureManagers & managers,
+                             SparsePoints & sparse_points,
+                             const std::string & representation_name) {
+        size_t n_centers{0};
+        // find the total number of gradients
+        for (const auto & manager : managers) {
+            n_centers += manager->size();
+        }
+        std::cout<<"centri"<<n_centers<<std::endl;
+        math::Matrix_t KNM(n_centers, sparse_points.size());
+        KNM.setZero();
+        size_t ii_A{0};
+        for (auto & manager : managers) {
+          auto && propA{*manager->template get_property<Property_t>(
+              representation_name, true)};
+          
+          for (auto center : manager) {
+            int sp = center.get_atom_type();
+            // only the pseudo points of species sp contribute
+            KNM.row(ii_A) =
+                pow_zeta(sparse_points.dot(sp, propA[center]), this->zeta)
+                    .transpose();
+            ++ii_A;
+          }
+          
+        }
+        return KNM;
+      }
+
+           /**
+       * Compute the kernel between a set of structure(s) and a set of pseudo
+       * points.
+       *
+       * @tparam StructureManagers should be an iterable over shared pointer
+       *          of structure managers like ManagerCollection
+       * @param sparse_points a SparsePoints* class
+       * @param managers a ManagerCollection or similar collection of
+       * structure managers
+       * @param representation_name name under which the representation data
+       * has been registered in the elements of managers and managers_b
+       * @return kernel matrix
+       */
+      template <class Property_t, internal::TargetType Type,
+                std::enable_if_t<Type == internal::TargetType::Atom, int> = 0,
+                class StructureManagers, class SparsePoints>
+      math::Matrix_t compute_local(const StructureManagers & managers,
+                             const SparsePoints & sparse_points,
+                             const std::string & representation_name) {
+        size_t n_centersA{0};
+        for (const auto & manager : managers) {
+          n_centersA += manager->size();
+        }
+        size_t nb_sparse_points{sparse_points.size()};
+        math::Matrix_t KNM(n_centersA, nb_sparse_points);
+        size_t ii_A{0};
+        for (auto & manager : managers) {
+          auto && propA{*manager->template get_property<Property_t>(
+              representation_name, true)};
+          for (auto center : manager) {
+            int sp = center.get_atom_type();
+            KNM.row(ii_A) =
+                pow_zeta(sparse_points.dot(sp, propA[center]), this->zeta)
+                    .transpose();
+            ii_A++;
+          }
+        }
+        return KNM;
+      }
+
       /**
        * Compute the kernel between a set of pseudo points.
        *
@@ -636,6 +723,66 @@ namespace rascal {
       }
     }
 
+    /**
+     * The root compute local kernel function. It computes the kernel between 2 set of
+     * structures for a given representation specified by the calculator.
+     *
+     * @param calculator the calculator which has been used to calculate
+     * the representation on the two managers
+     * has been registered in the elements of managers and managers_b
+     * @param managers a ManagerCollection or similar collection of
+     * structure managers
+     * @param managers_b a ManagerCollection or similar collection of
+     * structure managers or a set of pseudo points
+     */
+    template <class Calculator, class StructureManagers, class SparsePoints>
+    math::Matrix_t compute_local(const Calculator & calculator,
+                           const StructureManagers & managers,
+                           const SparsePoints & sparse_points) {
+      using ManagerPtr_t = typename StructureManagers::value_type;
+      using Manager_t = typename ManagerPtr_t::element_type;
+      using Property_t = typename Calculator::template Property_t<Manager_t>;
+      auto && representation_name{calculator.get_name()};
+      using internal::TargetType;
+      return this->compute_helper<Property_t, TargetType::Atom>(
+             representation_name, managers, sparse_points);
+      // switch (this->target_type) {
+      // case TargetType::Structure:
+      //   return this->compute_local_helper<Property_t, TargetType::Structure>(
+      //       representation_name, managers, sparse_points);
+      // case TargetType::Atom:
+      //   return this->compute_local_helper<Property_t, TargetType::Atom>(
+      //       representation_name, managers, sparse_points);
+      // default:
+      //   throw std::logic_error(
+      //       "Given target_type " +
+      //       this->parameters["target_type"].get<std::string>() +
+      //       " is not known."
+      //       " It is either 'Structure' or 'Atom')");
+      // }
+    }
+
+    template <class Property_t, internal::TargetType Type,
+              class StructureManagers, class SparsePoints>
+    math::Matrix_t compute_local_helper(const std::string & representation_name,
+                                  const StructureManagers & managers,
+                                  const SparsePoints & sparse_points) {
+      using internal::SparseKernelType;
+
+      if (this->kernel_type == SparseKernelType::GAP) {
+        auto kernel =
+            downcast_sparse_kernel_impl<SparseKernelType::GAP>(kernel_impl);
+        return kernel->template compute_local<Property_t, Type>(
+            managers, sparse_points, representation_name);
+      } else {
+        throw std::logic_error(
+            "Given kernel_type " +
+            this->parameters["kernel_type"].get<std::string>() +
+            " is not known."
+            " It is 'GAP'");
+      }
+    }
+
     /**
      * The root compute kernel function. It computes the kernel between the
      * representation gradients of a set of structures with the set of pseudo
diff --git a/tests/python/python_krr_test.py b/tests/python/python_krr_test.py
new file mode 100644
index 000000000..b0c0209c4
--- /dev/null
+++ b/tests/python/python_krr_test.py
@@ -0,0 +1,50 @@
+# TMP line for debugging for alex
+import sys
+sys.path.insert(0, "/home/goscinsk/code/librascal-local-stress/build")
+sys.path.insert(0, "/Users/davidetisi/Documents/EPFL_postdoc/librascal/build")
+
+import rascal
+print(rascal.__file__)
+import unittest
+# TMP END
+
+
+import json
+
+import numpy as np
+import ase.io
+
+from rascal.utils import load_obj
+
+class TestKRR(unittest.TestCase):
+    def setUp(self):
+        """
+        """
+        self.model = load_obj("reference_data/tests_only/simple_gap_model.json")
+        # the file is extracted from the "reference_data/tests_only/simple_gap_fit_params.json"
+        self.frames = ase.io.read("reference_data/inputs/methane_dimer_sample.xyz", ":")
+        for frame in self.frames:
+            frame.wrap(eps=1e-12)
+        # librascal stacks all central local stress computation for each
+        # structure into one dimension, so we need to keep track of the
+        # staring index for each structure
+        self.structure_starting_indices = np.cumsum([0] + [len(frame) for frame in self.frames])
+
+    def test_summed_local_stress_is_equal_to_global_stress(self):
+        """ Tests if the local stress contribution summed up match the global one"""
+        calculator = self.model.get_representation_calculator()
+        manager = calculator.transform(self.frames)
+        global_stress = self.model.predict_stress(manager)
+
+        # compare to local stress 
+        local_stress = self.model.predict_stress(manager, local_stress=True)
+        local_stress_matrix = local_stress.reshape(-1, 3, 3)
+        # voight indicies of stress matrix as stored in the global stress in librascal
+        #                 xx      yy     zz     yz     xz    xy
+        voigt_indices = [(0,0), (1,1), (2,2), (1,2), (0,2), (0,1)]
+        local_stress_voigt = np.concatenate([local_stress_matrix[:, i, j].reshape(-1,1) for (i,j) in voigt_indices], axis=1)
+        for i in range(len(self.frames)):
+            self.assertTrue(
+                np.allclose(global_stress[i],
+                np.sum(local_stress_voigt[self.structure_starting_indices[i]:self.structure_starting_indices[i+1]], axis=0)))
+
diff --git a/tests/python/python_krr_test_copy.py b/tests/python/python_krr_test_copy.py
new file mode 100644
index 000000000..475c55c39
--- /dev/null
+++ b/tests/python/python_krr_test_copy.py
@@ -0,0 +1,205 @@
+# TMP line for debugging for alex
+import sys
+sys.path.insert(0, "/scratch/tisi/Programs/librascal/")
+sys.path.insert(0, "/scratch/tisi/Programs/librascal/build")
+
+import rascal
+print(rascal.__file__)
+import unittest
+# TMP END
+
+
+import json
+
+import numpy as np
+import ase.io
+
+from rascal.utils import load_obj
+import copy
+
+
+def displace_strain_tensor(frame, alpha_index, beta_index, h_disp):
+    shift = np.eye(3)
+    shift[alpha_index, beta_index] += h_disp
+    original_cell = copy.deepcopy(frame.cell.array)
+    displaced_cell = np.dot(original_cell, shift)
+    frame.set_cell(displaced_cell)
+    # adapt scaled positions of atoms
+    M = np.linalg.inv(original_cell).dot(displaced_cell)
+    frame.positions = np.dot(frame.positions, M)
+    return frame
+
+class TestKRR(unittest.TestCase):
+    def setUp(self):
+        """
+        """
+        self.model = load_obj("../../reference_data/tests_only/simple_gap_model.json")
+        # the file is extracted from the "reference_data/tests_only/simple_gap_fit_params.json"
+        self.frames = ase.io.read("../../reference_data/inputs/methane_dimer_sample.xyz", ":")
+        for frame in self.frames:
+            frame.wrap(eps=1e-12)
+        # librascal stacks all central local stress computation for each
+        # structure into one dimension, so we need to keep track of the
+        # staring index for each structure
+        self.structure_starting_indices = np.cumsum([0] + [len(frame) for frame in self.frames])
+
+    def test_local_energies_sum_to_total_energy(self):
+        """ Tests if the local stress contribution summed up match the global one"""
+        calculator = self.model.get_representation_calculator()
+        manager = calculator.transform(self.frames)
+
+        KI = self.model.kernel._kernel.compute(
+                self.model.kernel._representation,
+                manager.managers,
+                self.model.X_train._sparse_points)
+        
+        KJ = self.model.kernel._kernel.compute_local(
+                self.model.kernel._representation,
+                manager.managers,
+                self.model.X_train._sparse_points)
+        
+        ei  = np.dot(KI, self.model.weights)
+        ej  = np.dot(KJ, self.model.weights)
+        
+        self.assertTrue(
+                np.allclose(KI,
+                np.sum(KJ.reshape(-1,10,30),axis=1)))
+        
+        self.assertTrue(np.allclose(ei , np.sum(ej.reshape(-1,10),axis=1)))
+        
+    def compute_numerical_stress(self,frame):
+        h_disp = 1e-6
+        matrix_indices_in_voigt_notation = [
+            (0, 0),
+            (0, 1),
+            (0, 2),
+            (1, 0),
+            (1, 1),
+            (1, 2),
+            (2, 0),
+            (2, 1),
+            (2, 2)
+        ]
+        
+        
+        numerical_stress = np.zeros((len(frame),9))
+        for i in range(9):
+            frame_displaced_plus = displace_strain_tensor(
+                copy.deepcopy(frame),
+                matrix_indices_in_voigt_notation[i][0],
+                matrix_indices_in_voigt_notation[i][1],
+                h_disp,
+            )
+            frame_displaced_minus = displace_strain_tensor(
+                copy.deepcopy(frame),
+                matrix_indices_in_voigt_notation[i][0],
+                matrix_indices_in_voigt_notation[i][1],
+                -h_disp,
+            )
+            calculator = self.model.get_representation_calculator()
+            manager_plus = calculator.transform(frame_displaced_plus)
+            manager_minus = calculator.transform(frame_displaced_minus)
+            
+            KJ_plus = self.model.kernel._kernel.compute_local(
+                    self.model.kernel._representation,
+                    manager_plus.managers,
+                    self.model.X_train._sparse_points)
+            
+            KJ_minus = self.model.kernel._kernel.compute_local(
+                    self.model.kernel._representation,
+                    manager_minus.managers,
+                    self.model.X_train._sparse_points)
+            
+            ej_plus  = np.dot(KJ_plus, self.model.weights)
+            ej_minus  = np.dot(KJ_minus, self.model.weights)
+            numerical_stress[:,i] = (ej_plus - ej_minus) / (2 * h_disp)
+        return numerical_stress / frame.get_volume()
+    
+    def test_local_stress_with_numerical_stress(self):
+        """ Tests if the local stress contribution summed up match the global one"""
+        
+        with open('output.log','w') as logfile:
+            calculator = self.model.get_representation_calculator()
+            i=0
+            for frame in self.frames[:]:
+                manager = calculator.transform(frame)
+
+                # compare to local stress 
+                local_stress = self.model.predict_stress(manager, local_stress=True)
+                local_stress_matrix = local_stress.reshape(-1, 3, 3)
+                print(i,file=logfile)
+                print(f'{local_stress_matrix=}',file=logfile)
+                numerical = self.compute_numerical_stress(frame).reshape(-1, 3, 3)
+                print(f'{numerical=}',file=logfile)
+                print(f'{numerical-local_stress_matrix=}',file=logfile)
+                print(f'{1e-08+1e-5*np.abs(numerical)=}',file=logfile)
+                #print(i)
+                #self.assertTrue(
+                # np.allclose(local_stress_matrix,
+                # numerical,atol=1e-10,rtol=1e-6))
+                i+=1
+        # for i in range(len(self.frames)):
+        #     #print(global_stress[i],np.sum(local_stress_voigt[self.structure_starting_indices[i]:self.structure_starting_indices[i+1]], axis=0))
+        #     self.assertTrue(
+        #         np.allclose(global_stress[i],
+        #         np.sum(local_stress_voigt[self.structure_starting_indices[i]:self.structure_starting_indices[i+1]], axis=0)))
+
+    
+    
+    def test_summed_local_stress_is_equal_to_global_stress(self):
+        """ Tests if the local stress contribution summed up match the global one"""
+        calculator = self.model.get_representation_calculator()
+        manager = calculator.transform(self.frames)
+        global_stress = self.model.predict_stress(manager)
+        # print(type(self.model))
+
+        # print(self.model.kernel.kernel_type)
+        # print("type( self.model.X_tra._sparse_points",type(self.model.X_train._sparse_points))
+        # #breakpoint()
+        # #self.model.kernel._kernel.compute(manager.managers,self.model.X_train._sparse_points)
+        # #breakpoint()
+        # KI = self.model.kernel._kernel.compute(
+        #         self.model.kernel._representation,
+        #         manager.managers,
+        #         self.model.X_train._sparse_points)
+        # mmodel = load_obj("../../reference_data/tests_only/simple_gap_model.json")
+        # mmodel.kernel.target_type = "Atom"
+        # mmodel.target_type = "Atom"
+        # print(f"{mmodel.kernel.target_type=}")
+        # print(f"{mmodel.target_type=}")
+        
+        # mcalculator = mmodel.get_representation_calculator()
+        # mmanager = mcalculator.transform(self.frames)
+        
+        # KJ = mmodel.kernel._kernel.compute_local(
+        #         mmodel.kernel._representation,
+        #         mmanager.managers,
+        #         mmodel.X_train._sparse_points)
+        
+        # ei  = np.dot(KI, self.model.weights)
+        # ej  = np.dot(KJ, self.model.weights)
+        
+        # print(f'{KI-np.sum(KJ.reshape(-1,10,30),axis=1)=}')
+        # print(f'{KI.shape=}')
+        # print(f'{KJ.shape=}')
+        # print(f'{len(self.frames)=}')
+        # print(f'{len(self.frames[0])=}')
+        # print(f'{ei.shape=}')
+        # print(f'{ei=}')
+        # print(f'{ej=}')
+        # print(f'{ei - np.sum(ej.reshape(-1,10),axis=1)=}')
+        
+
+        # compare to local stress 
+        local_stress = self.model.predict_stress(manager, local_stress=True)
+        local_stress_matrix = local_stress.reshape(-1, 3, 3)
+        # voight indicies of stress matrix as stored in the global stress in librascal
+        #                 xx      yy     zz     yz     xz    xy
+        voigt_indices = [(0,0), (1,1), (2,2), (1,2), (0,2), (0,1)]
+        local_stress_voigt = np.concatenate([local_stress_matrix[:, i, j].reshape(-1,1) for (i,j) in voigt_indices], axis=1)
+        for i in range(len(self.frames)):
+            #print(global_stress[i],np.sum(local_stress_voigt[self.structure_starting_indices[i]:self.structure_starting_indices[i+1]], axis=0))
+            self.assertTrue(
+                np.allclose(global_stress[i],
+                np.sum(local_stress_voigt[self.structure_starting_indices[i]:self.structure_starting_indices[i+1]], axis=0)))
+