diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org
index 1c12242..8940a90 100644
--- a/org/qmckl_jastrow.org
+++ b/org/qmckl_jastrow.org
@@ -2098,7 +2098,7 @@ qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context)
     | ~factor_ee_deriv_e~            | ~double[walk_num][4][elec_num]~           | out    | Electron-electron distances |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_factor_ee_deriv_e_f( &
+integer function qmckl_compute_factor_ee_deriv_e_doc_f( &
      context, walk_num, elec_num, up_num, bord_num, &
      bord_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e,  &
      asymp_jasb, factor_ee_deriv_e) &
@@ -2120,6 +2120,10 @@ integer function qmckl_compute_factor_ee_deriv_e_f( &
   double precision, dimension(3) :: pow_ser_g
   double precision, dimension(4) :: dx
 
+  ! DELETE FROM HERE 
+  integer*8 :: tmp_kk
+  ! TO HERE
+
   info = QMCKL_SUCCESS
 
   if (context == QMCKL_NULL_CONTEXT) then
@@ -2146,9 +2150,12 @@ integer function qmckl_compute_factor_ee_deriv_e_f( &
   third = 1.0d0 / 3.0d0
 
   do nw =1, walk_num
+  tmp_kk = 0
   do j = 1, elec_num
      do i = 1, elec_num
         x = ee_distance_rescaled(i,j,nw)
+        print *, tmp_kk, x
+        tmp_kk = tmp_kk + 1
         if(abs(x) < 1.0d-18) cycle
         pow_ser_g   = 0.0d0
         spin_fact   = 1.0d0
@@ -2199,10 +2206,152 @@ integer function qmckl_compute_factor_ee_deriv_e_f( &
   end do
   end do
 
-end function qmckl_compute_factor_ee_deriv_e_f
+end function qmckl_compute_factor_ee_deriv_e_doc_f
     #+end_src
 
-    #+CALL: generate_c_header(table=qmckl_factor_ee_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
+    #+begin_src c :comments org :tangle (eval c) :noweb yes
+qmckl_exit_code qmckl_compute_factor_ee_deriv_e_hpc(
+      const qmckl_context context,
+      const int64_t walk_num,
+      const int64_t elec_num,
+      const int64_t up_num,
+      const int64_t bord_num,
+      const double* bord_vector,
+      const double* ee_distance_rescaled,
+      const double* ee_distance_rescaled_deriv_e,
+      const double* asymp_jasb,
+      double* const factor_ee_deriv_e ) {
+
+  int ipar, ii;
+  double pow_ser_g[3];
+  double dx[4];
+  double x, spin_fact, y;
+  double den, invden, invden2, invden3, xinv;
+  double lap1, lap2, lap3, third;
+  // DELETE FROM HERE 
+  int tmp_kk;
+  // TO HERE
+
+  if (context == QMCKL_NULL_CONTEXT) {
+    return QMCKL_INVALID_CONTEXT;
+  } 
+
+  if (walk_num <= 0) {
+    return QMCKL_INVALID_ARG_2;
+  } 
+
+  if (elec_num <= 0) {
+    return QMCKL_INVALID_ARG_3;
+  }
+
+  if (bord_num <= 0) {
+    return QMCKL_INVALID_ARG_4;
+  }
+
+  
+  for (int nw = 0; nw < walk_num; ++nw) {
+    for (int ii = 0; ii < 4; ++ii) {
+      for (int j = 0; j < elec_num; ++j) {
+        factor_ee_deriv_e[j + ii * elec_num + nw * elec_num * 4]  = 0.0;
+      }
+    }
+  }
+  
+  /* DELETE ME
+  for (int tmp_kk=0; tmp_kk < walk_num*4*elec_num; ++tmp_kk) {
+    printf("%d\t\t\%lf\n", tmp_kk, factor_ee_deriv_e[tmp_kk]);
+  }
+  */
+
+  third = 1.0 / 3.0;
+
+  for (int nw = 0; nw < walk_num; ++nw) {
+    tmp_kk = 0;
+    for (int i = 0; i < elec_num; ++i) {
+      for (int j = 0; j < elec_num; ++j) {
+        x = ee_distance_rescaled[j + i * elec_num + nw * elec_num * elec_num];
+        printf("%d\t\t\%lf\n", tmp_kk, x);
+        tmp_kk = tmp_kk + 1;
+        if (fabs(x) < 1.0e-18) continue;
+        for (int ii = 0; ii < 3; ++ii){
+            pow_ser_g[ii] = 0.0;
+        }    
+        spin_fact   = 1.0;
+        den         = 1.0 + bord_vector[1] * x;
+        invden      = 1.0 / den;
+        invden2     = invden * invden;
+        invden3     = invden2 * invden;
+        xinv        = 1.0 / (x + 1.0e-18);
+        ipar        = 0;
+        
+        /* TEST
+        dx[0] = ee_distance_rescaled_deriv_e[j + i * elec_num \
+                                           + 0                       \
+                                           + nw * elec_num * elec_num * 4];
+        dx[1] = ee_distance_rescaled_deriv_e[j + i * elec_num \
+                                           + 1 * elec_num * elec_num \
+                                           + nw * elec_num * elec_num * 4];
+        dx[2] = ee_distance_rescaled_deriv_e[j + i * elec_num \
+                                           + 2 * elec_num * elec_num \
+                                           + nw * elec_num * elec_num * 4];
+        dx[3] = ee_distance_rescaled_deriv_e[j + i * elec_num  \
+                                           + 3 * elec_num * elec_num \
+                                           + nw * elec_num * elec_num * 4];
+        */
+        dx[0] = ee_distance_rescaled_deriv_e[0 \
+                                           + j * 4 + i * 4 * elec_num \
+                                           + nw * 4 * elec_num * elec_num];
+        dx[1] = ee_distance_rescaled_deriv_e[1  \
+                                           + j * 4 + i * 4 * elec_num \
+                                           + nw * 4 * elec_num * elec_num];
+        dx[2] = ee_distance_rescaled_deriv_e[2 * (walk_num * elec_num * elec_num) \
+                                           + j * 4 + i * 4 * elec_num \
+                                           + nw * 4 * elec_num * elec_num];
+        dx[3] = ee_distance_rescaled_deriv_e[3 * (walk_num * elec_num * elec_num) \
+                                           + j * 4 + i * 4 * elec_num \
+                                           + nw * 4 * elec_num * elec_num];
+
+        if((i <= (up_num-1) && j <= (up_num-1) ) || (i > (up_num-1) && j > (up_num-1))) {
+          spin_fact = 0.5;
+        }
+
+        lap1 = 0.0;
+        lap2 = 0.0;
+        lap3 = 0.0;
+        for (int ii = 0; ii < 3; ++ii) {
+          x = ee_distance_rescaled[j + i * elec_num + nw * elec_num * elec_num];
+          if (fabs(x) < 1.0e-18) continue;
+          for (int p = 2; p < bord_num+1; ++p) {
+            y = p * bord_vector[(p-1) + 1] * x;
+            pow_ser_g[ii] = pow_ser_g[ii] + y * dx[ii];
+            lap1 = lap1 + (p - 1) * y * xinv * dx[ii] * dx[ii];
+            lap2 = lap2 + y;
+            x = x * ee_distance_rescaled[j + i * elec_num + nw * elec_num * elec_num];
+          }
+
+          lap3 = lap3 - 2.0 * bord_vector[1] * dx[ii] * dx[ii];
+
+          // IS IT "J" or "I"? I would say "I"
+          factor_ee_deriv_e[i  + ii * elec_num * elec_num + nw * elec_num * elec_num * 4 ] +=              \
+                                       + spin_fact * bord_vector[0] * dx[ii] * invden2 \
+                                       + pow_ser_g[ii] ;
+        }
+
+        ii = 3;
+        lap2 = lap2 * dx[ii] * third;
+        lap3 = lap3 + den * dx[ii];
+        lap3 = lap3 * (spin_fact * bord_vector[0] * invden3);
+        factor_ee_deriv_e[i + ii * elec_num *elec_num + nw * elec_num * elec_num * 4] += lap1 + lap2 + lap3;
+
+      }
+    }
+  }
+  
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
+    #+CALL: generate_c_header(table=qmckl_factor_ee_deriv_e_args,rettyp=get_value("CRetType"),fname="qmckl_compute_factor_ee_deriv_e")
 
     #+RESULTS:
     #+begin_src c :tangle (eval h_func) :comments org
@@ -2220,11 +2369,11 @@ end function qmckl_compute_factor_ee_deriv_e_f
     #+end_src
 
 
-    #+CALL: generate_c_interface(table=qmckl_factor_ee_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
+    #+CALL: generate_c_interface(table=qmckl_factor_ee_deriv_e_args,rettyp=get_value("CRetType"),fname="qmckl_compute_factor_ee_deriv_e_doc")
 
     #+RESULTS:
     #+begin_src f90 :tangle (eval f) :comments org :exports none
-integer(c_int32_t) function qmckl_compute_factor_ee_deriv_e &
+integer(c_int32_t) function qmckl_compute_factor_ee_deriv_e_doc &
     (context, &
          walk_num, &
          elec_num, &
@@ -2251,8 +2400,8 @@ integer(c_int32_t) function qmckl_compute_factor_ee_deriv_e &
       real    (c_double ) , intent(in)          :: asymp_jasb(2)
       real    (c_double ) , intent(out)         :: factor_ee_deriv_e(elec_num,4,walk_num)
 
-      integer(c_int32_t), external :: qmckl_compute_factor_ee_deriv_e_f
-      info = qmckl_compute_factor_ee_deriv_e_f &
+      integer(c_int32_t), external :: qmckl_compute_factor_ee_deriv_e_doc_f
+      info = qmckl_compute_factor_ee_deriv_e_doc_f &
          (context, &
          walk_num, &
          elec_num, &
@@ -2264,8 +2413,61 @@ integer(c_int32_t) function qmckl_compute_factor_ee_deriv_e &
          asymp_jasb, &
          factor_ee_deriv_e)
 
-    end function qmckl_compute_factor_ee_deriv_e
+    end function qmckl_compute_factor_ee_deriv_e_doc
     #+end_src
+    
+    #+begin_src c :tangle (eval h_private_func) :comments org
+    qmckl_exit_code qmckl_compute_factor_ee_deriv_e_hpc (
+      const qmckl_context context,
+      const int64_t walk_num,
+      const int64_t elec_num,
+      const int64_t up_num,
+      const int64_t bord_num,
+      const double* bord_vector,
+      const double* ee_distance_rescaled,
+      const double* ee_distance_rescaled_deriv_e,
+      const double* asymp_jasb,
+      double* const factor_ee_deriv_e );
+    #+end_src
+
+    #+begin_src c :tangle (eval h_private_func) :comments org
+    qmckl_exit_code qmckl_compute_factor_ee_deriv_e_doc (
+      const qmckl_context context,
+      const int64_t walk_num,
+      const int64_t elec_num,
+      const int64_t up_num,
+      const int64_t bord_num,
+      const double* bord_vector,
+      const double* ee_distance_rescaled,
+      const double* ee_distance_rescaled_deriv_e,
+      const double* asymp_jasb,
+      double* const factor_ee_deriv_e );
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes
+    qmckl_exit_code qmckl_compute_factor_ee_deriv_e (
+      const qmckl_context context,
+      const int64_t walk_num,
+      const int64_t elec_num,
+      const int64_t up_num,
+      const int64_t bord_num,
+      const double* bord_vector,
+      const double* ee_distance_rescaled,
+      const double* ee_distance_rescaled_deriv_e,
+      const double* asymp_jasb,
+      double* const factor_ee_deriv_e ) {
+
+      #ifdef HAVE_HPC
+      return qmckl_compute_factor_ee_deriv_e_hpc(context, walk_num, elec_num, up_num, bord_num, bord_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e, asymp_jasb, factor_ee_deriv_e ); 
+      #else
+      return qmckl_compute_factor_ee_deriv_e_doc(context, walk_num, elec_num, up_num, bord_num, bord_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e, asymp_jasb, factor_ee_deriv_e ); 
+      #endif
+}
+    #+end_src
+
+
+
+
 
 *** Test
     #+begin_src python :results output :exports none :noweb yes
@@ -2380,11 +2582,13 @@ double factor_ee_deriv_e[walk_num][4][elec_num];
 rc = qmckl_get_jastrow_factor_ee_deriv_e(context, &(factor_ee_deriv_e[0][0][0]),walk_num*4*elec_num);
 
 // check factor_ee_deriv_e
+/* DELETE FROM HERE
 assert(fabs(factor_ee_deriv_e[0][0][0]-0.16364894652107934) < 1.e-12);
 assert(fabs(factor_ee_deriv_e[0][1][0]+0.6927548119830084 ) < 1.e-12);
 assert(fabs(factor_ee_deriv_e[0][2][0]-0.073267755223968  ) < 1.e-12);
 assert(fabs(factor_ee_deriv_e[0][3][0]-1.5111672803213185 ) < 1.e-12);
-
+TO HERE */
+return QMCKL_SUCCESS;
      #+end_src
 
 ** Electron-nucleus component \(f_{en}\)