diff --git a/gext/fitting.py b/gext/fitting.py
index dbc5a50f271b43f071d54ae06881c102c77087dd..d67c7c5a8879fbbdc1a81304c0f33d9ace1c0b70 100644
--- a/gext/fitting.py
+++ b/gext/fitting.py
@@ -44,7 +44,7 @@ class AbstractFitting(abc.ABC):
 
 class DiffFitting(AbstractFitting):
 
-    """Simple least square minimization fitting."""
+    """least square minimization fitting of the diffvectors"""
 
     supported_options = {
         "regularization": 0.0,
@@ -64,16 +64,17 @@ class DiffFitting(AbstractFitting):
         if len(vectors) == 1:
             raise ValueError("DiffFit does not work for one vector")
         target=target-vectors[-1]
+        target=target.flatten()
         VECTORS=[]
         for i in range(1, len(vectors)):
-            VECTORS.append(vectors[i-1]-vectors[-1])
+            each_vector=vectors[i-1]-vectors[-1]
+            VECTORS.append(each_vector.flatten())
         matrix = np.array(VECTORS).T
         a = matrix.T @ matrix
         b = matrix.T @ target
         if self.options["regularization"] > 0.0:
             a += np.identity(len(b))*self.options["regularization"]
         coefficients = np.linalg.solve(a, b)
-        print("coefficients", coefficients)
         return np.array(coefficients, dtype=np.float64)
 
     def linear_combination(self, vectors: List[np.ndarray],
diff --git a/gext/main.py b/gext/main.py
index b9ce8c4bb4b211e6e94a0ae639a1e3820e4b6772..680a24d4328e9313d2d10a91339b3f4167a63ffc 100644
--- a/gext/main.py
+++ b/gext/main.py
@@ -1,4 +1,4 @@
-"""module containing the Extrapolator class."""
+"""module containHing the Extrapolator class."""
 
 from typing import Optional
 import numpy as np
@@ -17,11 +17,11 @@ class Extrapolator:
     supported_options = {
         "verbose": False,
         "nsteps": 6,
-        "descriptor": "distance",
+        "descriptor": "H_core",
         "fitting": "diff",
         "allow_partially_filled": True,
         "store_overlap": True,
-        "tangent": "first",
+        "tangent": "one_before_last",
     }
 
     def __init__(self, nelectrons: int, nbasis: int, natoms: int, **kwargs):
@@ -46,7 +46,8 @@ class Extrapolator:
             self.coeffs=CircularBuffer(self.options["nsteps"],
                 (self.nelectrons//2, self.nbasis))
         self.tangent: Optional[np.ndarray] = None
-
+        self.H_cores=CircularBuffer(self.options["nsteps"],
+                (self.nbasis, self.nbasis))
     def set_options(self, **kwargs):
         """Given an arbitrary amount of keyword arguments, parse them if
         specified, set default values if not specified and raise an error
@@ -81,12 +82,15 @@ class Extrapolator:
             self.descriptor_calculator = Distance()
         elif self.options["descriptor"] == "coulomb":
             self.descriptor_calculator = Coulomb()
+        elif self.options["descriptor"] == "H_core":
+            pass
         else:
             raise ValueError("Unsupported descriptor")
-        self.descriptor_calculator.set_options(**descriptor_options)
+        if self.options["descriptor"] is not "H_core":
+            self.descriptor_calculator.set_options(**descriptor_options)
 
         if self.options["fitting"] == "leastsquare":
-            self.fitting_calculator = LeastSquare()
+            eelf.fitting_calculator = LeastSquare()
         elif self.options["fitting"] == "diff":
             self.fitting_calculator = DiffFitting()
         elif self.options["fitting"] == "qtr":
@@ -95,26 +99,25 @@ class Extrapolator:
             raise ValueError("Unsupported fitting")
         self.fitting_calculator.set_options(**fitting_options)
 
-    def load_data(self, coords: np.ndarray, coeff: np.ndarray, overlap):
+    def load_data(self, H_core: np.ndarray, coeff: np.ndarray, overlap):
         """Load a new data point in the extrapolator."""
 
         # Crop the coefficient matrix up to the number of electron
         # pairs, then apply S^1/2
         coeff = self._crop_coeff(coeff)
         coeff = self._normalize(coeff, overlap)
-
+        if self.options["tangent"]=="one_before_last":
+            self.coeffs.push(coeff)
+        else:
+            self.gammas.push(self._grassmann_log(coeff))
         # if it is the first time we load data, set the tangent point
         if self.tangent is None and self.options["tangent"] is not "one_before_last":
             self._set_tangent(coeff)
-
-        # push the new data to the corresponding vectors
-
-        self.descriptors.push(self._compute_descriptor(coords))
-        if self.options["tangent"]=="one_before_last":
-            cropped_coeff = self._crop_coeff(coeff)
-            self.coeffs.push(cropped_coeff)
+        if self.options["descriptor"]=="H_core":
+            self.H_cores.push(H_core)
         else:
-            self.gammas.push(self._grassmann_log(coeff))
+        # push the new data to the corresponding vectors
+            self.descriptors.push(self._compute_descriptor(coords))
 
         if self.options["store_overlap"]:
             self.overlaps.push(overlap)
@@ -124,11 +127,15 @@ class Extrapolator:
         c_guess = self.guess_coefficients(coords, overlap)
         return c_guess @ c_guess.T
 
-    def guess_coefficients(self, coords: np.ndarray, overlap=None) -> np.ndarray:
+    def guess_coefficients(self, H_core: np.ndarray, overlap=None) -> np.ndarray:
         """Get a new coefficient matrix to be used as a guess."""
 
         # check if we have enough data points to perform an extrapolation
-        count = self.descriptors.count
+       
+        if self.options["descriptor"]=="H_core":
+            count=self.H_cores.count
+        else:
+            count = self.descriptors.count
         if self.options["allow_partially_filled"]:
             if count == 0:
                 raise ValueError("Not enough data loaded in the extrapolator")
@@ -142,10 +149,14 @@ class Extrapolator:
             raise ValueError("Guessing without overlap requires `store_overlap` true.")
 
         # use the descriptors to find the fitting coefficients
-        prev_descriptors = self.descriptors.get(n)
-        descriptor = self._compute_descriptor(coords)
-        fit_coefficients = self._fit(prev_descriptors, descriptor)
-        print(fit_coefficients)
+        if self.options["descriptor"]=="H_core":
+            prev_H_cores = self.H_cores.get(n)
+            fit_coefficients = self._fit(prev_H_cores, H_core)
+        else:
+            prev_descriptors= self.descriptors.get(n)
+            descriptor = self._compute_descriptor(coords)
+            fit_coefficients = self._fit(prev_descriptors, descriptor)
+            print(fit_coefficients)
 
         # use the fitting coefficients and the previous gammas to
         # extrapolate a new gamma
@@ -199,7 +210,7 @@ class Extrapolator:
 
     def _set_tangent(self, c: np.ndarray):
         """Set the tangent point."""
-        if self.options["tangent"]="one_before_last"
+        if self.options["tangent"]=="one_before_last":
             self.tangent = c
         else:
             if self.tangent is None: