diff --git a/org/qmckl_jastrow_champ.org b/org/qmckl_jastrow_champ.org
index ca8e7b1..b0702e4 100644
--- a/org/qmckl_jastrow_champ.org
+++ b/org/qmckl_jastrow_champ.org
@@ -1300,13 +1300,6 @@ qmckl_get_jastrow_champ_a_vector (const qmckl_context context,
                            NULL);
   }
 
-  if (a_vector == NULL) {
-    return qmckl_failwith( context,
-                           QMCKL_INVALID_ARG_2,
-                           "qmckl_get_jastrow_champ_a_vector",
-                           "a_vector is a null pointer");
-  }
-
   qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
@@ -1316,15 +1309,26 @@ qmckl_get_jastrow_champ_a_vector (const qmckl_context context,
     return QMCKL_NOT_PROVIDED;
   }
 
+  if (a_vector == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_a_vector",
+                           "Null pointer");
+  }
+
   assert (ctx->jastrow_champ.a_vector != NULL);
-  int64_t sze = (ctx->jastrow_champ.aord_num + 1)*ctx->jastrow_champ.type_nucl_num;
+
+  const int64_t sze = (ctx->jastrow_champ.aord_num + 1)*ctx->jastrow_champ.type_nucl_num;
+
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
                            "qmckl_get_jastrow_champ_a_vector",
                            "Array too small. Expected (aord_num + 1)*type_nucl_num");
   }
+
   memcpy(a_vector, ctx->jastrow_champ.a_vector, sze*sizeof(double));
+
   return QMCKL_SUCCESS;
 }
 
@@ -1340,13 +1344,6 @@ qmckl_get_jastrow_champ_b_vector (const qmckl_context context,
                            NULL);
   }
 
-  if (b_vector == NULL) {
-    return qmckl_failwith( context,
-                           QMCKL_INVALID_ARG_2,
-                           "qmckl_get_jastrow_champ_b_vector",
-                           "b_vector is a null pointer");
-  }
-
   qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
@@ -1357,14 +1354,25 @@ qmckl_get_jastrow_champ_b_vector (const qmckl_context context,
   }
 
   assert (ctx->jastrow_champ.b_vector != NULL);
-  int64_t sze=ctx->jastrow_champ.bord_num +1;
+
+  if (b_vector == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_b_vector",
+                           "Null pointer");
+  }
+
+  const int64_t sze=ctx->jastrow_champ.bord_num +1;
+
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
                            "qmckl_get_jastrow_champ_b_vector",
                            "Array too small. Expected bord_num + 1");
   }
+
   memcpy(b_vector, ctx->jastrow_champ.b_vector, sze*sizeof(double));
+
   return QMCKL_SUCCESS;
 }
 
@@ -1380,13 +1388,6 @@ qmckl_get_jastrow_champ_c_vector (const qmckl_context context,
                            NULL);
   }
 
-  if (c_vector == NULL) {
-    return qmckl_failwith( context,
-                           QMCKL_INVALID_ARG_2,
-                           "qmckl_get_jastrow_champ_c_vector",
-                           "c_vector is a null pointer");
-  }
-
   qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
@@ -1402,20 +1403,29 @@ qmckl_get_jastrow_champ_c_vector (const qmckl_context context,
   qmckl_exit_code rc = qmckl_get_jastrow_champ_dim_c_vector(context, &dim_c_vector);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  int64_t sze=dim_c_vector * ctx->jastrow_champ.type_nucl_num;
+  if (c_vector == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_c_vector",
+                           "c_vector is a null pointer");
+  }
+
+  const int64_t sze=dim_c_vector * ctx->jastrow_champ.type_nucl_num;
+
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
                            "qmckl_get_jastrow_champ_c_vector",
                            "Array too small. Expected dim_c_vector*type_nucl_num");
   }
+
   memcpy(c_vector, ctx->jastrow_champ.c_vector, sze*sizeof(double));
   return QMCKL_SUCCESS;
 }
 
 qmckl_exit_code
 qmckl_get_jastrow_champ_rescale_factor_ee (const qmckl_context context,
-                                      double* const rescale_factor_ee) {
+                                           double* const rescale_factor_ee) {
 
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return qmckl_failwith( context,
@@ -1450,6 +1460,7 @@ qmckl_exit_code
 qmckl_get_jastrow_champ_rescale_factor_en (const qmckl_context context,
                                       double* const rescale_factor_en,
                                       const int64_t size_max) {
+
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_CONTEXT,
@@ -1810,7 +1821,7 @@ qmckl_get_jastrow_champ_asymp_jasb(qmckl_context context,
   qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
-  int64_t sze = 2;
+  const int64_t sze = 2;
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
@@ -2154,6 +2165,744 @@ assert(fabs(asymp_jasb[1]-1.043287918508297 ) < 1.e-12);
 
       #+end_src
 
+*** Electron-electron rescaled distances
+
+    ~ee_distance_rescaled~ stores the matrix of the rescaled distances between all
+    pairs of electrons:
+
+    \[
+    C_{ij} = \frac{ 1 - e^{-\kappa r_{ij}}}{\kappa}
+    \]
+
+   where \(r_{ij}\) is the matrix of electron-electron distances.
+
+**** Get
+
+     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code qmckl_get_jastrow_champ_ee_distance_rescaled(qmckl_context context,
+                                                             double* const distance_rescaled,
+                                                             int64_t const max_size);
+     #+end_src
+
+     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_get_jastrow_champ_ee_distance_rescaled(qmckl_context context,
+                                                             double* const distance_rescaled,
+                                                             int64_t const size_max)
+{
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_exit_code rc;
+
+  rc = qmckl_provide_ee_distance_rescaled(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  if (distance_rescaled == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_provide_jastrow_champ_factor_ee_gl",
+                           "Null pointer");
+  }
+
+  const int64_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walker.num;
+
+  if (size_max < sze) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_provide_jastrow_champ_factor_ee_gl",
+                           "Array too small. Expected elec_num*elec_num*walk_num.");
+  }
+  memcpy(distance_rescaled, ctx->jastrow_champ.ee_distance_rescaled, sze * sizeof(double));
+
+  return QMCKL_SUCCESS;
+}
+     #+end_src
+
+**** Provide                                                       :noexport:
+
+     #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
+qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context);
+     #+end_src
+
+     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+
+  /* Compute if necessary */
+  if (ctx->electron.walker.point.date > ctx->jastrow_champ.ee_distance_rescaled_date) {
+
+    if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
+      if (ctx->jastrow_champ.ee_distance_rescaled != NULL) {
+        qmckl_exit_code rc = qmckl_free(context, ctx->jastrow_champ.ee_distance_rescaled);
+        if (rc != QMCKL_SUCCESS) {
+          return qmckl_failwith( context, rc,
+                                 "qmckl_provide_ee_distance_rescaled",
+                                 "Unable to free ctx->jastrow_champ.ee_distance_rescaled");
+        }
+        ctx->jastrow_champ.ee_distance_rescaled = NULL;
+      }
+    }
+
+    /* Allocate array */
+    if (ctx->jastrow_champ.ee_distance_rescaled == NULL) {
+
+      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+      mem_info.size = ctx->electron.num * ctx->electron.num *
+        ctx->electron.walker.num * sizeof(double);
+      double* ee_distance_rescaled = (double*) qmckl_malloc(context, mem_info);
+
+      if (ee_distance_rescaled == NULL) {
+        return qmckl_failwith( context,
+                               QMCKL_ALLOCATION_FAILED,
+                               "qmckl_provide_ee_distance_rescaled",
+                               NULL);
+      }
+      ctx->jastrow_champ.ee_distance_rescaled = ee_distance_rescaled;
+    }
+
+    qmckl_exit_code rc =
+      qmckl_compute_ee_distance_rescaled(context,
+                                ctx->electron.num,
+                                ctx->jastrow_champ.rescale_factor_ee,
+                                ctx->electron.walker.num,
+                                ctx->electron.walker.point.coord.data,
+                                ctx->jastrow_champ.ee_distance_rescaled);
+    if (rc != QMCKL_SUCCESS) {
+      return rc;
+    }
+
+    ctx->jastrow_champ.ee_distance_rescaled_date = ctx->date;
+  }
+
+  return QMCKL_SUCCESS;
+}
+     #+end_src
+
+**** Compute
+     :PROPERTIES:
+     :Name:     qmckl_compute_ee_distance_rescaled
+     :CRetType: qmckl_exit_code
+     :FRetType: qmckl_exit_code
+     :END:
+
+     #+NAME: qmckl_ee_distance_rescaled_args
+     | Variable            | Type                                   | In/Out | Description                          |
+     |---------------------+----------------------------------------+--------+--------------------------------------|
+     | ~context~           | ~qmckl_context~                        | in     | Global state                         |
+     | ~elec_num~          | ~int64_t~                              | in     | Number of electrons                  |
+     | ~rescale_factor_ee~ | ~double~                               | in     | Factor to rescale ee distances       |
+     | ~walk_num~          | ~int64_t~                              | in     | Number of walkers                    |
+     | ~coord~             | ~double[3][walk_num][elec_num]~        | in     | Electron coordinates                 |
+     | ~ee_distance~       | ~double[walk_num][elec_num][elec_num]~ | out    | Electron-electron rescaled distances |
+
+     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+function qmckl_compute_ee_distance_rescaled_doc(context, &
+     elec_num, rescale_factor_ee, walk_num, &
+     coord, ee_distance_rescaled) &
+     bind(C) result(info)
+  use qmckl
+  implicit none
+
+  integer(qmckl_context), intent(in), value :: context
+  integer (c_int64_t) , intent(in)  , value :: elec_num
+  real    (c_double ) , intent(in)  , value :: rescale_factor_ee
+  integer (c_int64_t) , intent(in)  , value :: walk_num
+  real    (c_double ) , intent(in)          :: coord(elec_num,walk_num,3)
+  real    (c_double ) , intent(out)         :: ee_distance_rescaled(elec_num,elec_num,walk_num)
+  integer(qmckl_exit_code)                  :: info
+
+  integer*8 :: k
+
+  info = QMCKL_SUCCESS
+
+  if (context == QMCKL_NULL_CONTEXT) then
+     info = QMCKL_INVALID_CONTEXT
+     return
+  endif
+
+  if (elec_num <= 0) then
+     info = QMCKL_INVALID_ARG_2
+     return
+  endif
+
+  if (walk_num <= 0) then
+     info = QMCKL_INVALID_ARG_3
+     return
+  endif
+
+  do k=1,walk_num
+     info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, elec_num, &
+          coord(1,k,1), elec_num * walk_num, &
+          coord(1,k,1), elec_num * walk_num, &
+          ee_distance_rescaled(1,1,k), elec_num, rescale_factor_ee)
+     if (info /= QMCKL_SUCCESS) then
+        exit
+     endif
+  end do
+
+end function qmckl_compute_ee_distance_rescaled_doc
+     #+end_src
+
+     #+begin_src c :tangle (eval h_private_func) :comments org :exports none
+qmckl_exit_code qmckl_compute_ee_distance_rescaled_doc (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const double rescale_factor_ee,
+          const int64_t walk_num,
+          const double* coord,
+          double* const ee_distance_rescaled );
+
+qmckl_exit_code qmckl_compute_ee_distance_rescaled_hpc (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const double rescale_factor_ee,
+          const int64_t walk_num,
+          const double* coord,
+          double* const ee_distance_rescaled );
+
+qmckl_exit_code qmckl_compute_ee_distance_rescaled (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const double rescale_factor_ee,
+          const int64_t walk_num,
+          const double* coord,
+          double* const ee_distance_rescaled );
+     #+end_src
+
+     #+begin_src c :tangle (eval c) :comments org :exports none
+qmckl_exit_code qmckl_compute_ee_distance_rescaled_hpc (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const double rescale_factor_ee,
+          const int64_t walk_num,
+          const double* coord,
+          double* const ee_distance_rescaled )
+{
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_INVALID_CONTEXT;
+  }
+
+  if (elec_num <= 0) {
+    return QMCKL_INVALID_ARG_2;
+  }
+
+  if (walk_num <= 0) {
+    return QMCKL_INVALID_ARG_4;
+  }
+
+  if (coord == NULL) {
+    return QMCKL_INVALID_ARG_5;
+  }
+
+  if (ee_distance_rescaled == NULL) {
+    return QMCKL_INVALID_ARG_6;
+  }
+
+
+  const int64_t sze = elec_num*walk_num;
+  const int64_t elec_num2= elec_num*elec_num;
+
+
+  qmckl_exit_code result = QMCKL_SUCCESS;
+#ifdef HAVE_OPENMP
+  #pragma omp parallel
+  {
+#endif
+    qmckl_exit_code rc = QMCKL_SUCCESS;
+    #pragma omp for
+    for (int64_t k=0 ; k<walk_num ; ++k)
+    {
+      rc |= qmckl_distance_rescaled(context, 'T', 'T', elec_num, elec_num,
+                                    &(coord[k*elec_num]), sze, &(coord[k*elec_num]), sze,
+                                    &(ee_distance_rescaled[k*elec_num2]), elec_num, rescale_factor_ee);
+    }
+    #pragma omp critical
+    result |= rc;
+#ifdef HAVE_OPENMP
+  }
+#endif
+  return result;
+}
+     #+end_src
+
+     #+begin_src c :tangle (eval c) :comments org :exports none
+qmckl_exit_code qmckl_compute_ee_distance_rescaled (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const double rescale_factor_ee,
+          const int64_t walk_num,
+          const double* coord,
+          double* const ee_distance_rescaled )
+{
+#ifdef HAVE_HPC
+  return qmckl_compute_ee_distance_rescaled_hpc
+#else
+  return qmckl_compute_ee_distance_rescaled_doc
+#endif
+ (context, elec_num, rescale_factor_ee, walk_num, coord, ee_distance_rescaled);
+}
+     #+end_src
+
+**** Test
+
+     #+begin_src python :results output :exports none
+import numpy as np
+
+kappa = 0.6
+
+elec_1_w1 = np.array( [-0.250655104764153,  0.503070975550133      ,  -0.166554344502303])
+elec_2_w1 = np.array( [-0.587812193472177, -0.128751981129274      ,   0.187773606533075])
+elec_5_w1 = np.array( [-0.127732483187947, -0.138975497694196      ,  -8.669850480215846E-002])
+elec_6_w1 = np.array( [-0.232271834949124, -1.059321673434182E-002 ,  -0.504862241464867])
+
+print ( "[0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-elec_1_w1)) )/kappa )
+print ( "[0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-elec_2_w1)) )/kappa )
+print ( "[1][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w1-elec_1_w1)) )/kappa )
+print ( "[5][5] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_5_w1-elec_5_w1)) )/kappa )
+print ( "[5][6] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_5_w1-elec_6_w1)) )/kappa )
+print ( "[6][5] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_6_w1-elec_5_w1)) )/kappa )
+     #+end_src
+
+     #+RESULTS:
+     : [0][0] :  0.0
+     : [0][1] :
+     : [1][0] :  0.6347507420688708
+     : [5][5] :  0.0
+     : [5][6] :  0.3941735387855409
+     : [6][5] :
+
+      #+begin_src c :tangle (eval c_test)
+assert(qmckl_electron_provided(context));
+
+
+double ee_distance_rescaled[walk_num * elec_num * elec_num];
+rc = qmckl_get_jastrow_champ_ee_distance_rescaled(context, ee_distance_rescaled,walk_num*elec_num*elec_num);
+
+// (e1,e2,w)
+// (0,0,0) == 0.
+assert(ee_distance_rescaled[0] == 0.);
+
+// (1,0,0) == (0,1,0)
+assert(ee_distance_rescaled[1] == ee_distance_rescaled[elec_num]);
+
+// value of (1,0,0)
+assert(fabs(ee_distance_rescaled[1]-0.6347507420688708) < 1.e-12);
+
+// (0,0,1) == 0.
+assert(ee_distance_rescaled[5*elec_num + 5] == 0.);
+
+// (1,0,1) == (0,1,1)
+assert(ee_distance_rescaled[5*elec_num+6] == ee_distance_rescaled[6*elec_num+5]);
+
+// value of (1,0,1)
+assert(fabs(ee_distance_rescaled[5*elec_num+6]-0.3941735387855409) < 1.e-12);
+
+      #+end_src
+*** Electron-electron rescaled distance gradients and Laplacian with respect to electron coordinates
+
+   The rescaled distances, represented by $C_{ij} = (1 - e^{-\kappa_\text{e} r_{ij}})/\kappa_\text{e}$
+   are differentiated with respect to the electron coordinates.
+   This information is stored in the tensor
+   ~ee_distance_rescaled_gl~. The initial three sequential
+   elements of this three-dimensional tensor provide the $x$, $y$, and $z$
+   direction derivatives, while the fourth index corresponds to the Laplacian.
+
+**** Get
+
+     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code qmckl_get_jastrow_champ_ee_distance_rescaled_gl(qmckl_context context,
+                                                                double* const distance_rescaled_gl,
+                                                                const int64_t size_max);
+     #+end_src
+
+     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_get_jastrow_champ_ee_distance_rescaled_gl(qmckl_context context,
+                                                                double* const distance_rescaled_gl,
+                                                                const int64_t size_max)
+{
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_exit_code rc;
+
+  rc = qmckl_provide_ee_distance_rescaled_gl(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  if (distance_rescaled_gl == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_ee_distance_rescaled_gl",
+                           "Null pointer.");
+  }
+
+  const int64_t sze = 4 * ctx->electron.num * ctx->electron.num * ctx->electron.walker.num;
+
+  if (size_max < sze) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_ee_distance_rescaled_gl",
+                           "Array too small. Expected 4*elec_num*elec_num*walk_num");
+  }
+
+  memcpy(distance_rescaled_gl, ctx->jastrow_champ.ee_distance_rescaled_gl, sze * sizeof(double));
+
+  return QMCKL_SUCCESS;
+}
+     #+end_src
+
+**** Provide                                                       :noexport:
+
+     #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
+qmckl_exit_code qmckl_provide_ee_distance_rescaled_gl(qmckl_context context);
+     #+end_src
+
+     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_provide_ee_distance_rescaled_gl(qmckl_context context)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+
+  /* Compute if necessary */
+  if (ctx->electron.walker.point.date > ctx->jastrow_champ.ee_distance_rescaled_gl_date) {
+
+    if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
+      if (ctx->jastrow_champ.ee_distance_rescaled_gl != NULL) {
+        qmckl_exit_code rc = qmckl_free(context, ctx->jastrow_champ.ee_distance_rescaled_gl);
+        if (rc != QMCKL_SUCCESS) {
+          return qmckl_failwith( context, rc,
+                                 "qmckl_provide_ee_distance_rescaled_gl",
+                                 "Unable to free ctx->jastrow_champ.ee_distance_rescaled_gl");
+        }
+        ctx->jastrow_champ.ee_distance_rescaled_gl = NULL;
+      }
+    }
+
+    /* Allocate array */
+    if (ctx->jastrow_champ.ee_distance_rescaled_gl == NULL) {
+
+      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+      mem_info.size = 4 * ctx->electron.num * ctx->electron.num *
+        ctx->electron.walker.num * sizeof(double);
+      double* ee_distance_rescaled_gl = (double*) qmckl_malloc(context, mem_info);
+
+      if (ee_distance_rescaled_gl == NULL) {
+        return qmckl_failwith( context,
+                               QMCKL_ALLOCATION_FAILED,
+                               "qmckl_provide_ee_distance_rescaled_gl",
+                               NULL);
+      }
+      ctx->jastrow_champ.ee_distance_rescaled_gl = ee_distance_rescaled_gl;
+    }
+
+    qmckl_exit_code rc =
+      qmckl_compute_ee_distance_rescaled_gl(context,
+                                ctx->electron.num,
+                                ctx->jastrow_champ.rescale_factor_ee,
+                                ctx->electron.walker.num,
+                                ctx->electron.walker.point.coord.data,
+                                ctx->jastrow_champ.ee_distance_rescaled_gl);
+    if (rc != QMCKL_SUCCESS) {
+      return rc;
+    }
+
+    ctx->jastrow_champ.ee_distance_rescaled_date = ctx->date;
+  }
+
+  return QMCKL_SUCCESS;
+}
+     #+end_src
+
+**** Compute
+     :PROPERTIES:
+     :Name:     qmckl_compute_ee_distance_rescaled_gl
+     :CRetType: qmckl_exit_code
+     :FRetType: qmckl_exit_code
+     :END:
+
+     #+NAME: qmckl_ee_distance_rescaled_gl_args
+     | Variable                  | Type                                      | In/Out | Description                                     |
+     |---------------------------+-------------------------------------------+--------+-------------------------------------------------|
+     | ~context~                 | ~qmckl_context~                           | in     | Global state                                    |
+     | ~elec_num~                | ~int64_t~                                 | in     | Number of electrons                             |
+     | ~rescale_factor_ee~       | ~double~                                  | in     | Factor to rescale ee distances                  |
+     | ~walk_num~                | ~int64_t~                                 | in     | Number of walkers                               |
+     | ~coord~                   | ~double[3][walk_num][elec_num]~           | in     | Electron coordinates                            |
+     | ~ee_distance_rescaled_gl~ | ~double[walk_num][elec_num][elec_num][4]~ | out    | Electron-electron rescaled distance derivatives |
+
+     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+function qmckl_compute_ee_distance_rescaled_gl_doc(context,  &
+     elec_num, rescale_factor_ee, walk_num, coord, ee_distance_rescaled_gl) &
+     bind(C) result(info)
+  use qmckl
+  implicit none
+
+  integer(qmckl_context), intent(in), value :: context
+  integer (c_int64_t) , intent(in)  , value :: elec_num
+  real    (c_double ) , intent(in)  , value :: rescale_factor_ee
+  integer (c_int64_t) , intent(in)  , value :: walk_num
+  real    (c_double ) , intent(in)          :: coord(elec_num,walk_num,3)
+  real    (c_double ) , intent(out)         :: ee_distance_rescaled_gl(4,elec_num,elec_num,walk_num)
+  integer(qmckl_exit_code)                  :: info
+
+  integer*8 :: k
+
+  info = QMCKL_SUCCESS
+
+  if (context == QMCKL_NULL_CONTEXT) then
+     info = QMCKL_INVALID_CONTEXT
+     return
+  endif
+
+  if (elec_num <= 0) then
+     info = QMCKL_INVALID_ARG_2
+     return
+  endif
+
+  if (walk_num <= 0) then
+     info = QMCKL_INVALID_ARG_3
+     return
+  endif
+
+  do k=1,walk_num
+     info = qmckl_distance_rescaled_gl(context, 'T', 'T', elec_num, elec_num, &
+          coord(1,k,1), elec_num*walk_num, &
+          coord(1,k,1), elec_num*walk_num, &
+          ee_distance_rescaled_gl(1,1,1,k), elec_num, rescale_factor_ee)
+     if (info /= QMCKL_SUCCESS) then
+        exit
+     endif
+  end do
+
+end function qmckl_compute_ee_distance_rescaled_gl_doc
+     #+end_src
+
+     #+begin_src c :tangle (eval h_private_func) :comments org :exports none
+qmckl_exit_code qmckl_compute_ee_distance_rescaled_gl_doc (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const double rescale_factor_ee,
+          const int64_t walk_num,
+          const double* coord,
+          double* const ee_distance_rescaled_gl );
+
+qmckl_exit_code qmckl_compute_ee_distance_rescaled_gl_hpc (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const double rescale_factor_ee,
+          const int64_t walk_num,
+          const double* coord,
+          double* const ee_distance_rescaled_gl );
+
+qmckl_exit_code qmckl_compute_ee_distance_rescaled_gl (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const double rescale_factor_ee,
+          const int64_t walk_num,
+          const double* coord,
+          double* const ee_distance_rescaled_gl );
+     #+end_src
+
+
+     #+begin_src c :tangle (eval c) :comments org :exports none
+qmckl_exit_code qmckl_compute_ee_distance_rescaled_gl_hpc (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const double rescale_factor_ee,
+          const int64_t walk_num,
+          const double* coord,
+          double* const ee_distance_rescaled_gl )
+{
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT;
+  if (elec_num <= 0) return QMCKL_INVALID_ARG_2;
+  if (walk_num <= 0) return QMCKL_INVALID_ARG_4;
+  if (coord == NULL) return QMCKL_INVALID_ARG_5;
+  if (ee_distance_rescaled_gl == NULL) return QMCKL_INVALID_ARG_6;
+
+  const int64_t sze = elec_num*walk_num;
+  const int64_t elec_num2= elec_num*elec_num*4;
+
+  qmckl_exit_code result = QMCKL_SUCCESS;
+  #pragma omp parallel
+  {
+    qmckl_exit_code rc = QMCKL_SUCCESS;
+    #pragma omp for
+    for (int64_t k=0 ; k<walk_num ; ++k)
+    {
+      rc |= qmckl_distance_rescaled_gl(context, 'T', 'T', elec_num, elec_num,
+                                       &(coord[k*elec_num]), sze,
+                                       &(coord[k*elec_num]), sze,
+                                       &(ee_distance_rescaled_gl[k*elec_num2]), elec_num,
+                                       rescale_factor_ee);
+    }
+    #pragma omp critical
+    result |= rc;
+  }
+  return result;
+}
+     #+end_src
+
+     #+begin_src c :tangle (eval c) :comments org :exports none
+qmckl_exit_code qmckl_compute_ee_distance_rescaled_gl (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const double rescale_factor_ee,
+          const int64_t walk_num,
+          const double* coord,
+          double* const ee_distance_rescaled_gl )
+{
+#ifdef HAVE_HPC
+  return qmckl_compute_ee_distance_rescaled_gl_hpc
+#else
+  return qmckl_compute_ee_distance_rescaled_gl_doc
+#endif
+   (context, elec_num, rescale_factor_ee, walk_num, coord,
+   ee_distance_rescaled_gl);
+}
+     #+end_src
+
+**** Test
+
+      #+begin_src c :tangle (eval c_test)
+assert(qmckl_electron_provided(context));
+
+{
+  printf("ee_distance_rescaled_gl\n");
+  double fd[walk_num][elec_num][elec_num][4];
+
+  double delta_x = 0.0001;
+
+  // Finite difference coefficients for gradients
+  double coef[9] = { 1.0/280.0, -4.0/105.0,  1.0/5.0, -4.0/5.0,         0.0, 4.0/5.0, -1.0/5.0, 4.0/105.0, -1.0/280.0 };
+
+  // Finite difference coefficients for Laplacian
+  double coef2[9]= {-1.0/560.0,	 8.0/315.0, -1.0/5.0,  8.0/5.0, -205.0/72.0, 8.0/5.0, -1.0/5.0,	8.0/315.0, -1.0/560.0 };
+
+  qmckl_exit_code rc;
+
+  int64_t walk_num;
+  rc = qmckl_get_electron_walk_num(context, &walk_num);
+  if (rc != QMCKL_SUCCESS) {
+    return rc;
+  }
+
+  int64_t elec_num;
+  rc = qmckl_get_electron_num(context, &elec_num);
+  if (rc != QMCKL_SUCCESS) {
+    return rc;
+  }
+
+  double elec_coord[walk_num][elec_num][3];
+  rc = qmckl_get_electron_coord (context, 'N', &(elec_coord[0][0][0]), 3*walk_num*elec_num);
+
+  double temp_coord[walk_num][elec_num][3];
+  memcpy(&(temp_coord[0][0][0]), &(elec_coord[0][0][0]), sizeof(temp_coord));
+
+  double function_values[walk_num][elec_num][elec_num];
+
+  memset(&(fd[0][0][0][0]), 0, sizeof(fd));
+
+  for (int64_t i = 0; i < elec_num; i++) {
+    for (int64_t k = 0; k < 3; k++) {
+      for (int64_t m = -4; m <= 4; m++) {          // Apply finite difference displacement
+          
+        for (int64_t nw=0 ; nw<walk_num ; nw++) {
+          temp_coord[nw][i][k] = elec_coord[nw][i][k] + (double) m * delta_x;
+        }
+
+        // Update coordinates in the context
+        rc = qmckl_set_electron_coord (context, 'N', walk_num, &(temp_coord[0][0][0]), walk_num*3*elec_num);
+        assert(rc == QMCKL_SUCCESS);
+
+        // Call the provided function
+        rc = qmckl_get_jastrow_champ_ee_distance_rescaled(context, &(function_values[0][0][0]), elec_num*elec_num*walk_num);
+        assert(rc == QMCKL_SUCCESS);
+
+        // Accumulate derivative using finite-difference coefficients
+        for (int64_t nw=0 ; nw<walk_num ; nw++) {
+          for (int64_t j = 0; j < elec_num; j++) {
+            fd[nw][j][i][k] += coef [m + 4] * function_values[nw][j][i];
+            fd[nw][j][i][3] += coef2[m + 4] * function_values[nw][j][i];
+          }
+        }
+      }
+      for (int64_t nw=0 ; nw<walk_num ; nw++) {
+        temp_coord[nw][i][k] = elec_coord[nw][i][k];
+      }
+    }
+  }
+
+  // Reset coordinates in the context
+  rc = qmckl_set_electron_coord (context, 'N', walk_num, &(elec_coord[0][0][0]), walk_num*3*elec_num);
+  assert(rc == QMCKL_SUCCESS);
+
+  // Normalize by the step size
+  for (int64_t nw=0 ; nw<walk_num ; nw++) {
+    for (int64_t i = 0; i < elec_num; i++) {
+      for (int64_t k = 0; k < 4; k++) {
+        for (int64_t j = 0; j < elec_num; j++) {
+          fd[nw][i][j][k] /= delta_x;
+        }
+      }
+      for (int64_t j = 0; j < elec_num; j++) {
+        fd[nw][i][j][3] /= delta_x;
+      }
+    }
+  }
+  
+
+  double ee_distance_rescaled_gl[walk_num][elec_num][elec_num][4];
+  rc = qmckl_check(context,
+                   qmckl_get_jastrow_champ_ee_distance_rescaled_gl(context,
+                                                                   &(ee_distance_rescaled_gl[0][0][0][0]),
+                                                                   walk_num*elec_num*4*elec_num)
+                   );
+
+  assert(rc == QMCKL_SUCCESS);
+
+  for (int nw = 0; nw < walk_num; nw++){
+    for (int i = 0; i < elec_num; i++) {
+      for (int j = 0; j < elec_num; j++) {
+        for (int k = 0; k < 3; k++){
+          printf("%.10f\t", fd[nw][i][j][k]);
+          printf("%.10f\n", ee_distance_rescaled_gl[nw][i][j][k]);
+          assert(fabs(fd[nw][i][j][k] - ee_distance_rescaled_gl[nw][i][j][k]) < 1.e-8);
+        }
+        int k=3;
+        if (i != j) {
+          printf("%.10f\t", fd[nw][i][j][k]);
+          printf("%.10f\n", ee_distance_rescaled_gl[nw][i][j][k]);
+          assert(fabs(fd[nw][i][j][k] - ee_distance_rescaled_gl[nw][i][j][k]) < 1.e-6);
+        }
+      }
+    }
+  }
+  printf("OK\n");
+}
+
+      #+end_src
+
 *** Electron-electron component
 
    Calculate the electron-electron jastrow component ~factor_ee~ using the ~asymp_jasb~
@@ -2205,7 +2954,15 @@ qmckl_get_jastrow_champ_factor_ee(qmckl_context context,
   rc = qmckl_provide_jastrow_champ_factor_ee(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  int64_t sze=ctx->electron.walker.num;
+  if (factor_ee == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_factor_ee",
+                           "Null pointer");
+  }
+
+  const int64_t sze = ctx->electron.walker.num;
+
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
@@ -2680,7 +3437,15 @@ qmckl_get_jastrow_champ_factor_ee_gl(qmckl_context context,
   qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
-  int64_t sze = ctx->electron.walker.num * 4 * ctx->electron.num;
+  if (factor_ee_gl == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_factor_ee_gl",
+                           "Null pointer");
+  }
+
+  const int64_t sze = ctx->electron.walker.num * 4 * ctx->electron.num;
+
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
@@ -3230,660 +3995,14 @@ print("factor_ee_gl[3][0]:",factor_ee_gl[3][0])
      : factor_ee_gl[2][0]:
      : factor_ee_gl[3][0]:
 
-
       #+begin_src c :tangle (eval c_test)
 /* Check if Jastrow is properly initialized */
 assert(qmckl_jastrow_champ_provided(context));
 
-// calculate factor_ee_gl
-double factor_ee_gl[walk_num][4][elec_num];
-rc = qmckl_get_jastrow_champ_factor_ee_gl(context, &(factor_ee_gl[0][0][0]),walk_num*4*elec_num);
-
-// check factor_ee_gl
-printf("%f  %f\n", factor_ee_gl[0][0][0], -0.39319353942687446);
-assert(fabs(factor_ee_gl[0][0][0]+0.39319353942687446) < 1.e-12);
-
-printf("%f  %f\n", factor_ee_gl[0][1][0], 1.0535615450668214);
-assert(fabs(factor_ee_gl[0][1][0]-1.0535615450668214) < 1.e-12);
-
-printf("%f  %f\n", factor_ee_gl[0][2][0],-0.39098406960784515);
-assert(fabs(factor_ee_gl[0][2][0]+0.39098406960784515) < 1.e-12);
-
-printf("%f  %f\n", factor_ee_gl[0][3][0],2.8650469630854483);
-assert(fabs(factor_ee_gl[0][3][0]-2.8650469630854483) < 1.e-12);
-
-      #+end_src
-
-*** Electron-electron rescaled distances
-
-    ~ee_distance_rescaled~ stores the matrix of the rescaled distances between all
-    pairs of electrons:
-
-    \[
-    C_{ij} = \frac{ 1 - e^{-\kappa r_{ij}}}{\kappa}
-    \]
-
-   where \(r_{ij}\) is the matrix of electron-electron distances.
-
-**** Get
-
-     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
-qmckl_exit_code qmckl_get_jastrow_champ_ee_distance_rescaled(qmckl_context context,
-                                                             double* const distance_rescaled,
-                                                             int64_t const max_size);
-     #+end_src
-
-     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_get_jastrow_champ_ee_distance_rescaled(qmckl_context context,
-                                                             double* const distance_rescaled,
-                                                             int64_t const size_max)
+//----
 {
-  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
-    return QMCKL_NULL_CONTEXT;
-  }
-
-  qmckl_exit_code rc;
-
-  rc = qmckl_provide_ee_distance_rescaled(context);
-  if (rc != QMCKL_SUCCESS) return rc;
-
-  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
-  assert (ctx != NULL);
-
-  if (distance_rescaled == NULL) {
-    return qmckl_failwith( context,
-                           QMCKL_INVALID_ARG_2,
-                           "qmckl_provide_jastrow_champ_factor_ee_gl",
-                           "Null pointer");
-  }
-
-  int64_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walker.num;
-
-  if (size_max < sze) {
-    return qmckl_failwith( context,
-                           QMCKL_INVALID_ARG_2,
-                           "qmckl_provide_jastrow_champ_factor_ee_gl",
-                           "Array too small. Expected elec_num*elec_num*walk_num.");
-  }
-  memcpy(distance_rescaled, ctx->jastrow_champ.ee_distance_rescaled, sze * sizeof(double));
-
-  return QMCKL_SUCCESS;
-}
-     #+end_src
-
-**** Provide                                                       :noexport:
-
-     #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
-qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context);
-     #+end_src
-
-     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context)
-{
-
-  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
-    return QMCKL_NULL_CONTEXT;
-  }
-
-  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
-  assert (ctx != NULL);
-
-
-  /* Compute if necessary */
-  if (ctx->electron.walker.point.date > ctx->jastrow_champ.ee_distance_rescaled_date) {
-
-    if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
-      if (ctx->jastrow_champ.ee_distance_rescaled != NULL) {
-        qmckl_exit_code rc = qmckl_free(context, ctx->jastrow_champ.ee_distance_rescaled);
-        if (rc != QMCKL_SUCCESS) {
-          return qmckl_failwith( context, rc,
-                                 "qmckl_provide_ee_distance_rescaled",
-                                 "Unable to free ctx->jastrow_champ.ee_distance_rescaled");
-        }
-        ctx->jastrow_champ.ee_distance_rescaled = NULL;
-      }
-    }
-
-    /* Allocate array */
-    if (ctx->jastrow_champ.ee_distance_rescaled == NULL) {
-
-      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
-      mem_info.size = ctx->electron.num * ctx->electron.num *
-        ctx->electron.walker.num * sizeof(double);
-      double* ee_distance_rescaled = (double*) qmckl_malloc(context, mem_info);
-
-      if (ee_distance_rescaled == NULL) {
-        return qmckl_failwith( context,
-                               QMCKL_ALLOCATION_FAILED,
-                               "qmckl_provide_ee_distance_rescaled",
-                               NULL);
-      }
-      ctx->jastrow_champ.ee_distance_rescaled = ee_distance_rescaled;
-    }
-
-    qmckl_exit_code rc =
-      qmckl_compute_ee_distance_rescaled(context,
-                                ctx->electron.num,
-                                ctx->jastrow_champ.rescale_factor_ee,
-                                ctx->electron.walker.num,
-                                ctx->electron.walker.point.coord.data,
-                                ctx->jastrow_champ.ee_distance_rescaled);
-    if (rc != QMCKL_SUCCESS) {
-      return rc;
-    }
-
-    ctx->jastrow_champ.ee_distance_rescaled_date = ctx->date;
-  }
-
-  return QMCKL_SUCCESS;
-}
-     #+end_src
-
-**** Compute
-     :PROPERTIES:
-     :Name:     qmckl_compute_ee_distance_rescaled
-     :CRetType: qmckl_exit_code
-     :FRetType: qmckl_exit_code
-     :END:
-
-     #+NAME: qmckl_ee_distance_rescaled_args
-     | Variable            | Type                                   | In/Out | Description                          |
-     |---------------------+----------------------------------------+--------+--------------------------------------|
-     | ~context~           | ~qmckl_context~                        | in     | Global state                         |
-     | ~elec_num~          | ~int64_t~                              | in     | Number of electrons                  |
-     | ~rescale_factor_ee~ | ~double~                               | in     | Factor to rescale ee distances       |
-     | ~walk_num~          | ~int64_t~                              | in     | Number of walkers                    |
-     | ~coord~             | ~double[3][walk_num][elec_num]~        | in     | Electron coordinates                 |
-     | ~ee_distance~       | ~double[walk_num][elec_num][elec_num]~ | out    | Electron-electron rescaled distances |
-
-     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-function qmckl_compute_ee_distance_rescaled_doc(context, &
-     elec_num, rescale_factor_ee, walk_num, &
-     coord, ee_distance_rescaled) &
-     bind(C) result(info)
-  use qmckl
-  implicit none
-
-  integer(qmckl_context), intent(in), value :: context
-  integer (c_int64_t) , intent(in)  , value :: elec_num
-  real    (c_double ) , intent(in)  , value :: rescale_factor_ee
-  integer (c_int64_t) , intent(in)  , value :: walk_num
-  real    (c_double ) , intent(in)          :: coord(elec_num,walk_num,3)
-  real    (c_double ) , intent(out)         :: ee_distance_rescaled(elec_num,elec_num,walk_num)
-  integer(qmckl_exit_code)                  :: info
-
-  integer*8 :: k
-
-  info = QMCKL_SUCCESS
-
-  if (context == QMCKL_NULL_CONTEXT) then
-     info = QMCKL_INVALID_CONTEXT
-     return
-  endif
-
-  if (elec_num <= 0) then
-     info = QMCKL_INVALID_ARG_2
-     return
-  endif
-
-  if (walk_num <= 0) then
-     info = QMCKL_INVALID_ARG_3
-     return
-  endif
-
-  do k=1,walk_num
-     info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, elec_num, &
-          coord(1,k,1), elec_num * walk_num, &
-          coord(1,k,1), elec_num * walk_num, &
-          ee_distance_rescaled(1,1,k), elec_num, rescale_factor_ee)
-     if (info /= QMCKL_SUCCESS) then
-        exit
-     endif
-  end do
-
-end function qmckl_compute_ee_distance_rescaled_doc
-     #+end_src
-
-     #+begin_src c :tangle (eval h_private_func) :comments org :exports none
-qmckl_exit_code qmckl_compute_ee_distance_rescaled_doc (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const double rescale_factor_ee,
-          const int64_t walk_num,
-          const double* coord,
-          double* const ee_distance_rescaled );
-
-qmckl_exit_code qmckl_compute_ee_distance_rescaled_hpc (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const double rescale_factor_ee,
-          const int64_t walk_num,
-          const double* coord,
-          double* const ee_distance_rescaled );
-
-qmckl_exit_code qmckl_compute_ee_distance_rescaled (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const double rescale_factor_ee,
-          const int64_t walk_num,
-          const double* coord,
-          double* const ee_distance_rescaled );
-     #+end_src
-
-     #+begin_src c :tangle (eval c) :comments org :exports none
-qmckl_exit_code qmckl_compute_ee_distance_rescaled_hpc (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const double rescale_factor_ee,
-          const int64_t walk_num,
-          const double* coord,
-          double* const ee_distance_rescaled )
-{
-  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
-    return QMCKL_INVALID_CONTEXT;
-  }
-
-  if (elec_num <= 0) {
-    return QMCKL_INVALID_ARG_2;
-  }
-
-  if (walk_num <= 0) {
-    return QMCKL_INVALID_ARG_4;
-  }
-
-  if (coord == NULL) {
-    return QMCKL_INVALID_ARG_5;
-  }
-
-  if (ee_distance_rescaled == NULL) {
-    return QMCKL_INVALID_ARG_6;
-  }
-
-
-  int64_t sze = elec_num*walk_num;
-  int64_t elec_num2= elec_num*elec_num;
-
-
-  qmckl_exit_code result = QMCKL_SUCCESS;
-#ifdef HAVE_OPENMP
-  #pragma omp parallel
-  {
-#endif
-    qmckl_exit_code rc = QMCKL_SUCCESS;
-    #pragma omp for
-    for (int64_t k=0 ; k<walk_num ; ++k)
-    {
-      rc |= qmckl_distance_rescaled(context, 'T', 'T', elec_num, elec_num,
-                                    &(coord[k*elec_num]), sze, &(coord[k*elec_num]), sze,
-                                    &(ee_distance_rescaled[k*elec_num2]), elec_num, rescale_factor_ee);
-    }
-    #pragma omp critical
-    result |= rc;
-#ifdef HAVE_OPENMP
-  }
-#endif
-  return result;
-}
-     #+end_src
-
-     #+begin_src c :tangle (eval c) :comments org :exports none
-qmckl_exit_code qmckl_compute_ee_distance_rescaled (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const double rescale_factor_ee,
-          const int64_t walk_num,
-          const double* coord,
-          double* const ee_distance_rescaled )
-{
-#ifdef HAVE_HPC
-  return qmckl_compute_ee_distance_rescaled_hpc
-#else
-  return qmckl_compute_ee_distance_rescaled_doc
-#endif
- (context, elec_num, rescale_factor_ee, walk_num, coord, ee_distance_rescaled);
-}
-     #+end_src
-
-**** Test
-
-     #+begin_src python :results output :exports none
-import numpy as np
-
-kappa = 0.6
-
-elec_1_w1 = np.array( [-0.250655104764153,  0.503070975550133      ,  -0.166554344502303])
-elec_2_w1 = np.array( [-0.587812193472177, -0.128751981129274      ,   0.187773606533075])
-elec_5_w1 = np.array( [-0.127732483187947, -0.138975497694196      ,  -8.669850480215846E-002])
-elec_6_w1 = np.array( [-0.232271834949124, -1.059321673434182E-002 ,  -0.504862241464867])
-
-print ( "[0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-elec_1_w1)) )/kappa )
-print ( "[0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-elec_2_w1)) )/kappa )
-print ( "[1][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w1-elec_1_w1)) )/kappa )
-print ( "[5][5] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_5_w1-elec_5_w1)) )/kappa )
-print ( "[5][6] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_5_w1-elec_6_w1)) )/kappa )
-print ( "[6][5] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_6_w1-elec_5_w1)) )/kappa )
-     #+end_src
-
-     #+RESULTS:
-     : [0][0] :  0.0
-     : [0][1] :
-     : [1][0] :  0.6347507420688708
-     : [5][5] :  0.0
-     : [5][6] :  0.3941735387855409
-     : [6][5] :
-
-      #+begin_src c :tangle (eval c_test)
-assert(qmckl_electron_provided(context));
-
-
-double ee_distance_rescaled[walk_num * elec_num * elec_num];
-rc = qmckl_get_jastrow_champ_ee_distance_rescaled(context, ee_distance_rescaled,walk_num*elec_num*elec_num);
-
-// (e1,e2,w)
-// (0,0,0) == 0.
-assert(ee_distance_rescaled[0] == 0.);
-
-// (1,0,0) == (0,1,0)
-assert(ee_distance_rescaled[1] == ee_distance_rescaled[elec_num]);
-
-// value of (1,0,0)
-assert(fabs(ee_distance_rescaled[1]-0.6347507420688708) < 1.e-12);
-
-// (0,0,1) == 0.
-assert(ee_distance_rescaled[5*elec_num + 5] == 0.);
-
-// (1,0,1) == (0,1,1)
-assert(ee_distance_rescaled[5*elec_num+6] == ee_distance_rescaled[6*elec_num+5]);
-
-// value of (1,0,1)
-assert(fabs(ee_distance_rescaled[5*elec_num+6]-0.3941735387855409) < 1.e-12);
-
-      #+end_src
-
-*** Electron-electron rescaled distance gradients and Laplacian with respect to electron coordinates
-
-   The rescaled distances, represented by $C_{ij} = (1 - e^{-\kappa_\text{e} r_{ij}})/\kappa_\text{e}$
-   are differentiated with respect to the electron coordinates.
-   This information is stored in the tensor
-   ~ee_distance_rescaled_gl~. The initial three sequential
-   elements of this three-dimensional tensor provide the $x$, $y$, and $z$
-   direction derivatives, while the fourth index corresponds to the Laplacian.
-
-**** Get
-
-     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
-qmckl_exit_code qmckl_get_jastrow_champ_ee_distance_rescaled_gl(qmckl_context context,
-                                                                double* const distance_rescaled_gl,
-                                                                const int64_t size_max);
-     #+end_src
-
-     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_get_jastrow_champ_ee_distance_rescaled_gl(qmckl_context context,
-                                                                double* const distance_rescaled_gl,
-                                                                const int64_t size_max)
-{
-  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
-    return QMCKL_NULL_CONTEXT;
-  }
-
-  qmckl_exit_code rc;
-
-  rc = qmckl_provide_ee_distance_rescaled_gl(context);
-  if (rc != QMCKL_SUCCESS) return rc;
-
-  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
-  assert (ctx != NULL);
-
-  if (distance_rescaled_gl == NULL) {
-    return qmckl_failwith( context,
-                           QMCKL_INVALID_ARG_2,
-                           "qmckl_get_jastrow_champ_ee_distance_rescaled_gl",
-                           "Null pointer.");
-  }
-  int64_t sze = 4 * ctx->electron.num * ctx->electron.num * ctx->electron.walker.num;
-  if (size_max < sze) {
-    return qmckl_failwith( context,
-                           QMCKL_INVALID_ARG_2,
-                           "qmckl_get_jastrow_champ_ee_distance_rescaled_gl",
-                           "Array too small. Expected 4*elec_num*elec_num*walk_num");
-  }
-  memcpy(distance_rescaled_gl, ctx->jastrow_champ.ee_distance_rescaled_gl, sze * sizeof(double));
-
-  return QMCKL_SUCCESS;
-}
-     #+end_src
-
-**** Provide                                                       :noexport:
-
-     #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
-qmckl_exit_code qmckl_provide_ee_distance_rescaled_gl(qmckl_context context);
-     #+end_src
-
-     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_provide_ee_distance_rescaled_gl(qmckl_context context)
-{
-
-  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
-    return QMCKL_NULL_CONTEXT;
-  }
-
-  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
-  assert (ctx != NULL);
-
-
-  /* Compute if necessary */
-  if (ctx->electron.walker.point.date > ctx->jastrow_champ.ee_distance_rescaled_gl_date) {
-
-    if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
-      if (ctx->jastrow_champ.ee_distance_rescaled_gl != NULL) {
-        qmckl_exit_code rc = qmckl_free(context, ctx->jastrow_champ.ee_distance_rescaled_gl);
-        if (rc != QMCKL_SUCCESS) {
-          return qmckl_failwith( context, rc,
-                                 "qmckl_provide_ee_distance_rescaled_gl",
-                                 "Unable to free ctx->jastrow_champ.ee_distance_rescaled_gl");
-        }
-        ctx->jastrow_champ.ee_distance_rescaled_gl = NULL;
-      }
-    }
-
-    /* Allocate array */
-    if (ctx->jastrow_champ.ee_distance_rescaled_gl == NULL) {
-
-      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
-      mem_info.size = 4 * ctx->electron.num * ctx->electron.num *
-        ctx->electron.walker.num * sizeof(double);
-      double* ee_distance_rescaled_gl = (double*) qmckl_malloc(context, mem_info);
-
-      if (ee_distance_rescaled_gl == NULL) {
-        return qmckl_failwith( context,
-                               QMCKL_ALLOCATION_FAILED,
-                               "qmckl_provide_ee_distance_rescaled_gl",
-                               NULL);
-      }
-      ctx->jastrow_champ.ee_distance_rescaled_gl = ee_distance_rescaled_gl;
-    }
-
-    qmckl_exit_code rc =
-      qmckl_compute_ee_distance_rescaled_gl(context,
-                                ctx->electron.num,
-                                ctx->jastrow_champ.rescale_factor_ee,
-                                ctx->electron.walker.num,
-                                ctx->electron.walker.point.coord.data,
-                                ctx->jastrow_champ.ee_distance_rescaled_gl);
-    if (rc != QMCKL_SUCCESS) {
-      return rc;
-    }
-
-    ctx->jastrow_champ.ee_distance_rescaled_date = ctx->date;
-  }
-
-  return QMCKL_SUCCESS;
-}
-     #+end_src
-
-**** Compute
-     :PROPERTIES:
-     :Name:     qmckl_compute_ee_distance_rescaled_gl
-     :CRetType: qmckl_exit_code
-     :FRetType: qmckl_exit_code
-     :END:
-
-     #+NAME: qmckl_ee_distance_rescaled_gl_args
-     | Variable                  | Type                                      | In/Out | Description                                     |
-     |---------------------------+-------------------------------------------+--------+-------------------------------------------------|
-     | ~context~                 | ~qmckl_context~                           | in     | Global state                                    |
-     | ~elec_num~                | ~int64_t~                                 | in     | Number of electrons                             |
-     | ~rescale_factor_ee~       | ~double~                                  | in     | Factor to rescale ee distances                  |
-     | ~walk_num~                | ~int64_t~                                 | in     | Number of walkers                               |
-     | ~coord~                   | ~double[3][walk_num][elec_num]~           | in     | Electron coordinates                            |
-     | ~ee_distance_rescaled_gl~ | ~double[walk_num][elec_num][elec_num][4]~ | out    | Electron-electron rescaled distance derivatives |
-
-     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-function qmckl_compute_ee_distance_rescaled_gl_doc(context,  &
-     elec_num, rescale_factor_ee, walk_num, coord, ee_distance_rescaled_gl) &
-     bind(C) result(info)
-  use qmckl
-  implicit none
-
-  integer(qmckl_context), intent(in), value :: context
-  integer (c_int64_t) , intent(in)  , value :: elec_num
-  real    (c_double ) , intent(in)  , value :: rescale_factor_ee
-  integer (c_int64_t) , intent(in)  , value :: walk_num
-  real    (c_double ) , intent(in)          :: coord(elec_num,walk_num,3)
-  real    (c_double ) , intent(out)         :: ee_distance_rescaled_gl(4,elec_num,elec_num,walk_num)
-  integer(qmckl_exit_code)                  :: info
-
-  integer*8 :: k
-
-  info = QMCKL_SUCCESS
-
-  if (context == QMCKL_NULL_CONTEXT) then
-     info = QMCKL_INVALID_CONTEXT
-     return
-  endif
-
-  if (elec_num <= 0) then
-     info = QMCKL_INVALID_ARG_2
-     return
-  endif
-
-  if (walk_num <= 0) then
-     info = QMCKL_INVALID_ARG_3
-     return
-  endif
-
-  do k=1,walk_num
-     info = qmckl_distance_rescaled_gl(context, 'T', 'T', elec_num, elec_num, &
-          coord(1,k,1), elec_num*walk_num, &
-          coord(1,k,1), elec_num*walk_num, &
-          ee_distance_rescaled_gl(1,1,1,k), elec_num, rescale_factor_ee)
-     if (info /= QMCKL_SUCCESS) then
-        exit
-     endif
-  end do
-
-end function qmckl_compute_ee_distance_rescaled_gl_doc
-     #+end_src
-
-     #+begin_src c :tangle (eval h_private_func) :comments org :exports none
-qmckl_exit_code qmckl_compute_ee_distance_rescaled_gl_doc (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const double rescale_factor_ee,
-          const int64_t walk_num,
-          const double* coord,
-          double* const ee_distance_rescaled_gl );
-
-qmckl_exit_code qmckl_compute_ee_distance_rescaled_gl_hpc (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const double rescale_factor_ee,
-          const int64_t walk_num,
-          const double* coord,
-          double* const ee_distance_rescaled_gl );
-
-qmckl_exit_code qmckl_compute_ee_distance_rescaled_gl (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const double rescale_factor_ee,
-          const int64_t walk_num,
-          const double* coord,
-          double* const ee_distance_rescaled_gl );
-     #+end_src
-
-
-     #+begin_src c :tangle (eval c) :comments org :exports none
-qmckl_exit_code qmckl_compute_ee_distance_rescaled_gl_hpc (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const double rescale_factor_ee,
-          const int64_t walk_num,
-          const double* coord,
-          double* const ee_distance_rescaled_gl )
-{
-  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT;
-  if (elec_num <= 0) return QMCKL_INVALID_ARG_2;
-  if (walk_num <= 0) return QMCKL_INVALID_ARG_4;
-  if (coord == NULL) return QMCKL_INVALID_ARG_5;
-  if (ee_distance_rescaled_gl == NULL) return QMCKL_INVALID_ARG_6;
-
-  int64_t sze = elec_num*walk_num;
-  int64_t elec_num2= elec_num*elec_num*4;
-
-  qmckl_exit_code result = QMCKL_SUCCESS;
-  #pragma omp parallel
-  {
-    qmckl_exit_code rc = QMCKL_SUCCESS;
-    #pragma omp for
-    for (int64_t k=0 ; k<walk_num ; ++k)
-    {
-      rc |= qmckl_distance_rescaled_gl(context, 'T', 'T', elec_num, elec_num,
-                                       &(coord[k*elec_num]), sze,
-                                       &(coord[k*elec_num]), sze,
-                                       &(ee_distance_rescaled_gl[k*elec_num2]), elec_num,
-                                       rescale_factor_ee);
-    }
-    #pragma omp critical
-    result |= rc;
-  }
-  return result;
-}
-     #+end_src
-
-     #+begin_src c :tangle (eval c) :comments org :exports none
-qmckl_exit_code qmckl_compute_ee_distance_rescaled_gl (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const double rescale_factor_ee,
-          const int64_t walk_num,
-          const double* coord,
-          double* const ee_distance_rescaled_gl )
-{
-#ifdef HAVE_HPC
-  return qmckl_compute_ee_distance_rescaled_gl_hpc
-#else
-  return qmckl_compute_ee_distance_rescaled_gl_doc
-#endif
-   (context, elec_num, rescale_factor_ee, walk_num, coord,
-   ee_distance_rescaled_gl);
-}
-     #+end_src
-
-**** Test
-
-     #+begin_src python :results output :exports none
-import numpy as np
-
-# TODO
-     #+end_src
-
-      #+begin_src c :tangle (eval c_test)
-assert(qmckl_electron_provided(context));
-
-{
-  printf("ee_distance_rescaled_gl\n");
-  double fd[walk_num][elec_num][elec_num][4];
-
+  printf("factor_ee_gl\n");
+  double fd[walk_num][4][elec_num];
   double delta_x = 0.0001;
 
   // Finite difference coefficients for gradients
@@ -3912,9 +4031,9 @@ assert(qmckl_electron_provided(context));
   double temp_coord[walk_num][elec_num][3];
   memcpy(&(temp_coord[0][0][0]), &(elec_coord[0][0][0]), sizeof(temp_coord));
 
-  double function_values[walk_num][elec_num][elec_num];
+  double function_values[walk_num];
 
-  memset(&(fd[0][0][0][0]), 0, sizeof(fd));
+  memset(&(fd[0][0][0]), 0, sizeof(fd));
 
   for (int64_t i = 0; i < elec_num; i++) {
     for (int64_t k = 0; k < 3; k++) {
@@ -3929,15 +4048,13 @@ assert(qmckl_electron_provided(context));
         assert(rc == QMCKL_SUCCESS);
 
         // Call the provided function
-        rc = qmckl_get_jastrow_champ_ee_distance_rescaled(context, &(function_values[0][0][0]), elec_num*elec_num*walk_num);
+        rc = qmckl_get_jastrow_champ_factor_ee(context, &(function_values[0]), walk_num);
         assert(rc == QMCKL_SUCCESS);
 
         // Accumulate derivative using finite-difference coefficients
         for (int64_t nw=0 ; nw<walk_num ; nw++) {
-          for (int64_t j = 0; j < elec_num; j++) {
-            fd[nw][j][i][k] += coef [m + 4] * function_values[nw][j][i];
-            fd[nw][j][i][3] += coef2[m + 4] * function_values[nw][j][i];
-          }
+          fd[nw][k][i] += coef [m + 4] * function_values[nw];
+          fd[nw][3][i] += coef2[m + 4] * function_values[nw];
         }
       }
       for (int64_t nw=0 ; nw<walk_num ; nw++) {
@@ -3952,43 +4069,37 @@ assert(qmckl_electron_provided(context));
 
   // Normalize by the step size
   for (int64_t nw=0 ; nw<walk_num ; nw++) {
+    for (int64_t k = 0; k < 4; k++) {
+      for (int64_t i = 0; i < elec_num; i++) {
+          fd[nw][k][i] /= delta_x;
+      }
+    }
     for (int64_t i = 0; i < elec_num; i++) {
-      for (int64_t k = 0; k < 4; k++) {
-        for (int64_t j = 0; j < elec_num; j++) {
-          fd[nw][i][j][k] /= delta_x;
-        }
-      }
-      for (int64_t j = 0; j < elec_num; j++) {
-        fd[nw][i][j][3] /= delta_x;
-      }
+      fd[nw][3][i] /= delta_x;
     }
   }
   
 
-  double ee_distance_rescaled_gl[walk_num][elec_num][elec_num][4];
+  double factor_ee_gl[walk_num][4][elec_num];
   rc = qmckl_check(context,
-                   qmckl_get_jastrow_champ_ee_distance_rescaled_gl(context,
-                                                                   &(ee_distance_rescaled_gl[0][0][0][0]),
-                                                                   walk_num*elec_num*4*elec_num)
+                   qmckl_get_jastrow_champ_factor_ee_gl(context,
+                                                        &(factor_ee_gl[0][0][0]),
+                                                        walk_num*4*elec_num)
                    );
 
   assert(rc == QMCKL_SUCCESS);
 
   for (int nw = 0; nw < walk_num; nw++){
     for (int i = 0; i < elec_num; i++) {
-      for (int j = 0; j < elec_num; j++) {
-        for (int k = 0; k < 3; k++){
-          printf("%.10f\t", fd[nw][i][j][k]);
-          printf("%.10f\n", ee_distance_rescaled_gl[nw][i][j][k]);
-          assert(fabs(fd[nw][i][j][k] - ee_distance_rescaled_gl[nw][i][j][k]) < 1.e-8);
-        }
-        int k=3;
-        if (i != j) {
-          printf("%.10f\t", fd[nw][i][j][k]);
-          printf("%.10f\n", ee_distance_rescaled_gl[nw][i][j][k]);
-          assert(fabs(fd[nw][i][j][k] - ee_distance_rescaled_gl[nw][i][j][k]) < 1.e-6);
-        }
+      for (int k = 0; k < 3; k++){
+        printf("%.10f\t", fd[nw][k][i]);
+        printf("%.10f\n", factor_ee_gl[nw][k][i]);
+        assert(fabs(fd[nw][k][i] - factor_ee_gl[nw][k][i]) < 1.e-8);
       }
+      int k=3;
+      printf("%.10f\t", fd[nw][k][i]);
+      printf("%.10f\n", factor_ee_gl[nw][k][i]);
+      assert(fabs(fd[nw][k][i] - factor_ee_gl[nw][k][i]) < 1.e-5);
     }
   }
   printf("OK\n");
@@ -4034,18 +4145,27 @@ qmckl_get_jastrow_champ_asymp_jasa(qmckl_context context,
   qmckl_exit_code rc;
   rc = qmckl_provide_jastrow_champ_asymp_jasa(context);
   if(rc != QMCKL_SUCCESS) return rc;
-  */
+  ,*/
 
   qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
-  int64_t sze = ctx->jastrow_champ.type_nucl_num;
+  if (asymp_jasa == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_asymp_jasa",
+                           "Null pointer");
+  }
+
+  const int64_t sze = ctx->jastrow_champ.type_nucl_num;
+
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
                            "qmckl_get_jastrow_champ_asymp_jasa",
                            "Array too small. Expected nucl_num");
   }
+
   memcpy(asymp_jasa, ctx->jastrow_champ.asymp_jasa, sze * sizeof(double));
 
   return QMCKL_SUCCESS;
@@ -4260,6 +4380,887 @@ assert(fabs(-1.75529774 - asymp_jasa[0]) < 1.e-8);
 
       #+end_src
 
+*** Electron-nucleus rescaled distances
+
+    ~en_distance_rescaled~ stores the matrix of the rescaled distances between
+    electrons and nuclei.
+
+    \[
+    C_{i\alpha} = \frac{ 1 - e^{-\kappa_\alpha R_{i\alpha}}}{\kappa_\alpha}
+    \]
+
+   where \(R_{i\alpha}\) is the matrix of electron-nucleus distances.
+
+**** Get
+
+     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code
+qmckl_get_jastrow_champ_en_distance_rescaled(qmckl_context context,
+                                             double* const distance_rescaled,
+                                             const int64_t size_max);
+     #+end_src
+
+
+     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code
+qmckl_get_jastrow_champ_en_distance_rescaled(qmckl_context context,
+                                             double* const distance_rescaled,
+                                             const int64_t size_max)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_exit_code rc;
+
+  rc = qmckl_provide_en_distance_rescaled(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  if (distance_rescaled == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_en_distance_rescaled",
+                           "Null pointer");
+  }
+  
+  const int64_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num;
+
+  if (size_max < sze) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_get_jastrow_champ_en_distance_rescaled",
+                           "Array too small. Expected walk_num*nucl_num*elec_num");
+  }
+
+  memcpy(distance_rescaled, ctx->jastrow_champ.en_distance_rescaled, sze * sizeof(double));
+
+  return QMCKL_SUCCESS;
+}
+     #+end_src
+
+**** Provide                                                       :noexport:
+
+     #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
+qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context);
+     #+end_src
+
+     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  if (!(ctx->nucleus.provided)) {
+    return QMCKL_NOT_PROVIDED;
+  }
+
+  /* Compute if necessary */
+  if (ctx->electron.walker.point.date > ctx->jastrow_champ.en_distance_rescaled_date) {
+
+    if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
+      if (ctx->jastrow_champ.en_distance_rescaled != NULL) {
+        qmckl_exit_code rc = qmckl_free(context, ctx->jastrow_champ.en_distance_rescaled);
+        if (rc != QMCKL_SUCCESS) {
+          return qmckl_failwith( context, rc,
+                                 "qmckl_provide_en_distance_rescaled",
+                                 "Unable to free ctx->jastrow_champ.en_distance_rescaled");
+        }
+        ctx->jastrow_champ.en_distance_rescaled = NULL;
+      }
+    }
+
+    /* Allocate array */
+    if (ctx->jastrow_champ.en_distance_rescaled == NULL) {
+
+      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+      mem_info.size = ctx->electron.num * ctx->nucleus.num *
+        ctx->electron.walker.num * sizeof(double);
+      double* en_distance_rescaled = (double*) qmckl_malloc(context, mem_info);
+
+      if (en_distance_rescaled == NULL) {
+        return qmckl_failwith( context,
+                               QMCKL_ALLOCATION_FAILED,
+                               "qmckl_provide_en_distance_rescaled",
+                               NULL);
+      }
+      ctx->jastrow_champ.en_distance_rescaled = en_distance_rescaled;
+    }
+
+    qmckl_exit_code rc =
+      qmckl_compute_en_distance_rescaled(context,
+                                ctx->electron.num,
+                                ctx->nucleus.num,
+                                ctx->jastrow_champ.type_nucl_num,
+                                ctx->jastrow_champ.type_nucl_vector,
+                                ctx->jastrow_champ.rescale_factor_en,
+                                ctx->electron.walker.num,
+                                ctx->electron.walker.point.coord.data,
+                                ctx->nucleus.coord.data,
+                                ctx->jastrow_champ.en_distance_rescaled);
+
+    if (rc != QMCKL_SUCCESS) {
+      return rc;
+    }
+
+    ctx->jastrow_champ.en_distance_rescaled_date = ctx->date;
+  }
+
+  return QMCKL_SUCCESS;
+}
+     #+end_src
+
+**** Compute
+     :PROPERTIES:
+     :Name:     qmckl_compute_en_distance_rescaled
+     :CRetType: qmckl_exit_code
+     :FRetType: qmckl_exit_code
+     :END:
+
+    #+NAME: qmckl_en_distance_rescaled_args
+     | Variable               | Type                                   | In/Out | Description                       |
+     |------------------------+----------------------------------------+--------+-----------------------------------|
+     | ~context~              | ~qmckl_context~                        | in     | Global state                      |
+     | ~elec_num~             | ~int64_t~                              | in     | Number of electrons               |
+     | ~nucl_num~             | ~int64_t~                              | in     | Number of nuclei                  |
+     | ~type_nucl_num~        | ~int64_t~                              | in     | Number of types of nuclei         |
+     | ~type_nucl_vector~     | ~int64_t[nucl_num]~                    | in     | Number of types of nuclei         |
+     | ~rescale_factor_en~    | ~double[type_nucl_num]~                | in     | The factor for rescaled distances |
+     | ~walk_num~             | ~int64_t~                              | in     | Number of walkers                 |
+     | ~elec_coord~           | ~double[3][walk_num][elec_num]~        | in     | Electron coordinates              |
+     | ~nucl_coord~           | ~double[3][elec_num]~                  | in     | Nuclear coordinates               |
+     | ~en_distance_rescaled~ | ~double[walk_num][nucl_num][elec_num]~ | out    | Electron-nucleus distances        |
+
+     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+function qmckl_compute_en_distance_rescaled_doc(context, &
+     elec_num, nucl_num, type_nucl_num, &
+     type_nucl_vector, rescale_factor_en, walk_num, elec_coord, &
+     nucl_coord, en_distance_rescaled) &
+     bind(C) result(info)
+  use qmckl
+  implicit none
+  integer (qmckl_context), intent(in), value :: context
+  integer (c_int64_t) , intent(in)  , value :: elec_num
+  integer (c_int64_t) , intent(in)  , value :: nucl_num
+  integer (c_int64_t) , intent(in)  , value :: type_nucl_num
+  integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
+  real    (c_double ) , intent(in)          :: rescale_factor_en(type_nucl_num)
+  integer (c_int64_t) , intent(in)  , value :: walk_num
+  real    (c_double ) , intent(in)          :: elec_coord(elec_num,walk_num,3)
+  real    (c_double ) , intent(in)          :: nucl_coord(nucl_num,3)
+  real    (c_double ) , intent(out)         :: en_distance_rescaled(elec_num,nucl_num,walk_num)
+  integer(qmckl_exit_code)                   :: info
+
+  integer*8 :: i, k
+  double precision      :: coord(3)
+
+  info = QMCKL_SUCCESS
+
+  if (context == QMCKL_NULL_CONTEXT) then
+     info = QMCKL_INVALID_CONTEXT
+     return
+  endif
+
+  if (elec_num <= 0) then
+     info = QMCKL_INVALID_ARG_2
+     return
+  endif
+
+  if (nucl_num <= 0) then
+     info = QMCKL_INVALID_ARG_3
+     return
+  endif
+
+  if (walk_num <= 0) then
+     info = QMCKL_INVALID_ARG_5
+     return
+  endif
+
+  do i=1, nucl_num
+     coord(1:3) = nucl_coord(i,1:3)
+     do k=1,walk_num
+        info = qmckl_distance_rescaled(context, 'T', 'N', elec_num, 1_8, &
+             elec_coord(1,k,1), elec_num*walk_num, coord, 3_8, &
+             en_distance_rescaled(1,i,k), elec_num, rescale_factor_en(type_nucl_vector(i)+1))
+        if (info /= QMCKL_SUCCESS) then
+           return
+        endif
+     end do
+  end do
+
+end function qmckl_compute_en_distance_rescaled_doc
+     #+end_src
+
+     #+begin_src c :tangle (eval c) :comments org :exports none
+qmckl_exit_code qmckl_compute_en_distance_rescaled_hpc (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const int64_t nucl_num,
+          const int64_t type_nucl_num,
+          const int64_t* type_nucl_vector,
+          const double* rescale_factor_en,
+          const int64_t walk_num,
+          const double* elec_coord,
+          const double* nucl_coord,
+          double* const en_distance_rescaled )
+{
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT;
+  if (elec_num <= 0) return QMCKL_INVALID_ARG_2;
+  if (nucl_num <= 0) return QMCKL_INVALID_ARG_3;
+  if (type_nucl_num <= 0) return QMCKL_INVALID_ARG_4;
+  if (type_nucl_vector == NULL) return QMCKL_INVALID_ARG_5;
+  if (rescale_factor_en == NULL) return QMCKL_INVALID_ARG_6;
+  if (walk_num <= 0) return QMCKL_INVALID_ARG_7;
+  if (elec_coord == NULL) return QMCKL_INVALID_ARG_8;
+  if (nucl_coord == NULL) return QMCKL_INVALID_ARG_9;
+  if (en_distance_rescaled == NULL) return QMCKL_INVALID_ARG_10;
+
+  const int64_t sze = elec_num*walk_num;
+
+  qmckl_exit_code result = QMCKL_SUCCESS;
+  #pragma omp parallel
+  {
+    qmckl_exit_code rc = QMCKL_SUCCESS;
+    #pragma omp for
+    for (int64_t k=0 ; k<walk_num ; ++k)
+    {
+      for (int64_t a=0 ; a<nucl_num ; ++a) {
+        const double coord[3] = { nucl_coord[a], nucl_coord[a+nucl_num], nucl_coord[a+2*nucl_num] };
+        rc |= qmckl_distance_rescaled(context, 'T', 'N', elec_num, 1,
+                                      &(elec_coord[k*elec_num]), sze,
+                                      coord, 3,
+                                      &(en_distance_rescaled[elec_num*(a+nucl_num*k)]), elec_num,
+                                      rescale_factor_en[type_nucl_vector[a]]);
+      }
+    }
+    #pragma omp critical
+    result |= rc;
+  }
+  return result;
+}
+     #+end_src
+
+     #+begin_src c :tangle (eval h_private_func) :comments org :exports none
+qmckl_exit_code qmckl_compute_en_distance_rescaled_doc (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const int64_t nucl_num,
+          const int64_t type_nucl_num,
+          const int64_t* type_nucl_vector,
+          const double*  rescale_factor_en,
+          const int64_t walk_num,
+          const double* elec_coord,
+          const double* nucl_coord,
+          double* const en_distance_rescaled );
+
+qmckl_exit_code qmckl_compute_en_distance_rescaled_hpc (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const int64_t nucl_num,
+          const int64_t type_nucl_num,
+          const int64_t* type_nucl_vector,
+          const double* rescale_factor_en,
+          const int64_t walk_num,
+          const double* elec_coord,
+          const double* nucl_coord,
+          double* const en_distance_rescaled );
+
+qmckl_exit_code qmckl_compute_en_distance_rescaled (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const int64_t nucl_num,
+          const int64_t type_nucl_num,
+          const int64_t* type_nucl_vector,
+          const double*  rescale_factor_en,
+          const int64_t walk_num,
+          const double* elec_coord,
+          const double* nucl_coord,
+          double* const en_distance_rescaled );
+     #+end_src
+
+     #+begin_src c :tangle (eval c) :comments org :exports none
+qmckl_exit_code qmckl_compute_en_distance_rescaled (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const int64_t nucl_num,
+          const int64_t type_nucl_num,
+          const int64_t* type_nucl_vector,
+          const double* rescale_factor_en,
+          const int64_t walk_num,
+          const double* elec_coord,
+          const double* nucl_coord,
+          double* const en_distance_rescaled )
+{
+#ifdef HAVE_HPC
+  return qmckl_compute_en_distance_rescaled_hpc
+#else
+  return qmckl_compute_en_distance_rescaled_doc
+#endif
+  (context, elec_num, nucl_num, type_nucl_num, type_nucl_vector,
+  rescale_factor_en, walk_num, elec_coord, nucl_coord, en_distance_rescaled );
+}
+    #+end_src
+
+**** Test
+
+     #+begin_src python :results output :exports none
+import numpy as np
+
+kappa = 0.6
+
+elec_1_w1 = np.array( [-0.250655104764153,  0.503070975550133      ,  -0.166554344502303])
+elec_2_w1 = np.array( [-0.587812193472177, -0.128751981129274      ,   0.187773606533075])
+elec_5_w1 = np.array( [-0.127732483187947, -0.138975497694196      ,  -8.669850480215846E-002])
+elec_6_w1 = np.array( [-0.232271834949124, -1.059321673434182E-002 ,  -0.504862241464867])
+nucl_1    = np.array( [ 0., 0., 0. ])
+nucl_2    = np.array( [ 0., 0., 2.059801 ])
+
+print ( "[0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-nucl_1)) )/kappa )
+print ( "[1][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-nucl_2)) )/kappa )
+print ( "[0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w1-nucl_1)) )/kappa )
+print ( "[0][5] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_5_w1-nucl_1)) )/kappa )
+print ( "[1][5] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_5_w1-nucl_2)) )/kappa )
+print ( "[0][6] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_6_w1-nucl_1)) )/kappa )
+
+     #+end_src
+
+     #+RESULTS:
+     : [0][0] :  0.4942158656729477
+     : [1][0] :  1.2464137498005765
+     : [0][1] :  0.5248654474756858
+     : [0][5] :  0.19529459944794733
+     : [1][5] :  1.2091967687767369
+     : [0][6] :  0.4726452953409436
+
+      #+begin_src c :tangle (eval c_test)
+
+assert(qmckl_electron_provided(context));
+assert(qmckl_nucleus_provided(context));
+
+double en_distance_rescaled[walk_num][nucl_num][elec_num];
+
+rc = qmckl_check(context,
+                 qmckl_get_jastrow_champ_en_distance_rescaled(context,
+                                                              &(en_distance_rescaled[0][0][0]),
+                                                              walk_num*nucl_num*elec_num)
+                 );
+assert (rc == QMCKL_SUCCESS);
+
+// (e,n,w) in Fortran notation
+// (1,1,1)
+assert(fabs(en_distance_rescaled[0][0][0] - 0.4942158656729477) < 1.e-12);
+// (1,2,1)
+assert(fabs(en_distance_rescaled[0][1][0] - 1.2464137498005765) < 1.e-12);
+// (2,1,1)
+assert(fabs(en_distance_rescaled[0][0][1] - 0.5248654474756858) < 1.e-12);
+// (1,1,2)
+assert(fabs(en_distance_rescaled[0][0][5] - 0.19529459944794733) < 1.e-12);
+// (1,2,2)
+assert(fabs(en_distance_rescaled[0][1][5] - 1.2091967687767369) < 1.e-12);
+// (2,1,2)
+assert(fabs(en_distance_rescaled[0][0][6] - 0.4726452953409436) < 1.e-12);
+
+
+      #+end_src
+
+*** Electron-electron rescaled distance gradients and Laplacian with respect to electron coordinates
+
+   The rescaled distances, represented by $C_{i\alpha} = (1 - e^{-\kappa_\alpha R_{i\alpha}})/\kappa$
+   are differentiated with respect to the electron coordinates.
+   This information is stored in the tensor
+   ~en_distance_rescaled_gl~. The initial three sequential
+   elements of this three-index tensor provide the $x$, $y$, and $z$
+   direction derivatives, while the fourth index corresponds to the Laplacian.
+
+**** Get
+
+     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code
+qmckl_get_jastrow_champ_en_distance_rescaled_gl(qmckl_context context,
+                                                double* const distance_rescaled_gl,
+                                                const int64_t size_max);
+     #+end_src
+
+     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code
+qmckl_get_jastrow_champ_en_distance_rescaled_gl(qmckl_context context,
+                                                double* const distance_rescaled_gl,
+                                                const int64_t size_max)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_exit_code rc;
+
+  rc = qmckl_provide_en_distance_rescaled_gl(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  if (distance_rescaled_gl == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_en_distance_rescaled_gl",
+                           "Null pointer");
+  }
+
+  const int64_t sze = 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num;
+
+  if (size_max < sze) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_get_jastrow_champ_en_distance_rescaled_gl",
+                           "Array too small. Expected walk_num*elec_num*4");
+  }
+
+  memcpy(distance_rescaled_gl, ctx->jastrow_champ.en_distance_rescaled_gl, sze * sizeof(double));
+
+  return QMCKL_SUCCESS;
+}
+     #+end_src
+
+**** Provide                                                       :noexport:
+
+     #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
+qmckl_exit_code qmckl_provide_en_distance_rescaled_gl(qmckl_context context);
+     #+end_src
+
+     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_provide_en_distance_rescaled_gl(qmckl_context context)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  if (!(ctx->nucleus.provided)) {
+    return QMCKL_NOT_PROVIDED;
+  }
+
+  /* Compute if necessary */
+  if (ctx->electron.walker.point.date > ctx->jastrow_champ.en_distance_rescaled_gl_date) {
+
+    if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
+      if (ctx->jastrow_champ.en_distance_rescaled_gl != NULL) {
+        qmckl_exit_code rc = qmckl_free(context, ctx->jastrow_champ.en_distance_rescaled_gl);
+        if (rc != QMCKL_SUCCESS) {
+          return qmckl_failwith( context, rc,
+                                 "qmckl_provide_en_distance_rescaled_gl",
+                                 "Unable to free ctx->jastrow_champ.en_distance_rescaled_gl");
+        }
+        ctx->jastrow_champ.en_distance_rescaled_gl = NULL;
+      }
+    }
+
+    /* Allocate array */
+    if (ctx->jastrow_champ.en_distance_rescaled_gl == NULL) {
+
+      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+      mem_info.size = 4 * ctx->electron.num * ctx->nucleus.num *
+        ctx->electron.walker.num * sizeof(double);
+      double* en_distance_rescaled_gl = (double*) qmckl_malloc(context, mem_info);
+
+      if (en_distance_rescaled_gl == NULL) {
+        return qmckl_failwith( context,
+                               QMCKL_ALLOCATION_FAILED,
+                               "qmckl_provide_en_distance_rescaled_gl",
+                               NULL);
+      }
+      ctx->jastrow_champ.en_distance_rescaled_gl = en_distance_rescaled_gl;
+    }
+
+    qmckl_exit_code rc =
+      qmckl_compute_en_distance_rescaled_gl(context,
+                                ctx->electron.num,
+                                ctx->nucleus.num,
+                                ctx->jastrow_champ.type_nucl_num,
+                                ctx->jastrow_champ.type_nucl_vector,
+                                ctx->jastrow_champ.rescale_factor_en,
+                                ctx->electron.walker.num,
+                                ctx->electron.walker.point.coord.data,
+                                ctx->nucleus.coord.data,
+                                ctx->jastrow_champ.en_distance_rescaled_gl);
+    if (rc != QMCKL_SUCCESS) {
+      return rc;
+    }
+
+    ctx->jastrow_champ.en_distance_rescaled_gl_date = ctx->date;
+  }
+
+  return QMCKL_SUCCESS;
+}
+     #+end_src
+
+**** Compute
+     :PROPERTIES:
+     :Name:     qmckl_compute_en_distance_rescaled_gl
+     :CRetType: qmckl_exit_code
+     :FRetType: qmckl_exit_code
+     :END:
+
+     #+NAME: qmckl_en_distance_rescaled_gl_args
+     | Variable                  | Type                                      | In/Out | Description                           |
+     |---------------------------+-------------------------------------------+--------+---------------------------------------|
+     | ~context~                 | ~qmckl_context~                           | in     | Global state                          |
+     | ~elec_num~                | ~int64_t~                                 | in     | Number of electrons                   |
+     | ~nucl_num~                | ~int64_t~                                 | in     | Number of nuclei                      |
+     | ~type_nucl_num~           | ~int64_t~                                 | in     | Number of nucleus types               |
+     | ~type_nucl_vector~        | ~int64_t[nucl_num]~                       | in     | Array of nucleus types                |
+     | ~rescale_factor_en~       | ~double[nucl_num]~                        | in     | The factors for rescaled distances    |
+     | ~walk_num~                | ~int64_t~                                 | in     | Number of walkers                     |
+     | ~elec_coord~              | ~double[3][walk_num][elec_num]~           | in     | Electron coordinates                  |
+     | ~nucl_coord~              | ~double[3][elec_num]~                     | in     | Nuclear coordinates                   |
+     | ~en_distance_rescaled_gl~ | ~double[walk_num][nucl_num][elec_num][4]~ | out    | Electron-nucleus distance derivatives |
+
+     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+integer function qmckl_compute_en_distance_rescaled_gl_doc_f(context, elec_num, nucl_num, &
+     type_nucl_num, type_nucl_vector, rescale_factor_en, walk_num, elec_coord, &
+     nucl_coord, en_distance_rescaled_gl) &
+     result(info)
+  use qmckl
+  implicit none
+  integer(qmckl_context), intent(in)  :: context
+  integer*8             , intent(in)  :: elec_num
+  integer*8             , intent(in)  :: nucl_num
+  integer*8             , intent(in)  :: type_nucl_num
+  integer*8             , intent(in)  :: type_nucl_vector(nucl_num)
+  double precision      , intent(in)  :: rescale_factor_en(nucl_num)
+  integer*8             , intent(in)  :: walk_num
+  double precision      , intent(in)  :: elec_coord(elec_num,walk_num,3)
+  double precision      , intent(in)  :: nucl_coord(nucl_num,3)
+  double precision      , intent(out) :: en_distance_rescaled_gl(4,elec_num,nucl_num,walk_num)
+
+  integer*8 :: i, k
+  double precision :: coord(3)
+
+  info = QMCKL_SUCCESS
+
+  if (context == QMCKL_NULL_CONTEXT) then
+     info = QMCKL_INVALID_CONTEXT
+     return
+  endif
+
+  if (elec_num <= 0) then
+     info = QMCKL_INVALID_ARG_2
+     return
+  endif
+
+  if (nucl_num <= 0) then
+     info = QMCKL_INVALID_ARG_3
+     return
+  endif
+
+  if (walk_num <= 0) then
+     info = QMCKL_INVALID_ARG_5
+     return
+  endif
+
+  do i=1, nucl_num
+     coord(1:3) = nucl_coord(i,1:3)
+     do k=1,walk_num
+        info = qmckl_distance_rescaled_gl(context, 'T', 'T', elec_num, 1_8, &
+             elec_coord(1,k,1), elec_num*walk_num, coord, 1_8, &
+             en_distance_rescaled_gl(1,1,i,k), elec_num, rescale_factor_en(type_nucl_vector(i)+1))
+        if (info /= QMCKL_SUCCESS) then
+           return
+        endif
+     end do
+  end do
+
+end function qmckl_compute_en_distance_rescaled_gl_doc_f
+     #+end_src
+
+     #+begin_src c :tangle (eval h_private_func) :comments org :exports none
+qmckl_exit_code qmckl_compute_en_distance_rescaled_gl_doc (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const int64_t nucl_num,
+          const int64_t type_nucl_num,
+          int64_t* const type_nucl_vector,
+          const double*  rescale_factor_en,
+          const int64_t walk_num,
+          const double* elec_coord,
+          const double* nucl_coord,
+          double* const en_distance_rescaled_gl );
+
+qmckl_exit_code qmckl_compute_en_distance_rescaled_gl_hpc (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const int64_t nucl_num,
+          const int64_t type_nucl_num,
+          int64_t* const type_nucl_vector,
+          const double*  rescale_factor_en,
+          const int64_t walk_num,
+          const double* elec_coord,
+          const double* nucl_coord,
+          double* const en_distance_rescaled_gl );
+
+qmckl_exit_code qmckl_compute_en_distance_rescaled_gl (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const int64_t nucl_num,
+          const int64_t type_nucl_num,
+          int64_t* const type_nucl_vector,
+          const double*  rescale_factor_en,
+          const int64_t walk_num,
+          const double* elec_coord,
+          const double* nucl_coord,
+          double* const en_distance_rescaled_gl );
+     #+end_src
+
+     #+begin_src c :tangle (eval c) :comments org :exports none
+qmckl_exit_code qmckl_compute_en_distance_rescaled_gl_hpc (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const int64_t nucl_num,
+          const int64_t type_nucl_num,
+          int64_t* const type_nucl_vector,
+          const double*  rescale_factor_en,
+          const int64_t walk_num,
+          const double* elec_coord,
+          const double* nucl_coord,
+          double* const en_distance_rescaled_gl )
+{
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT;
+  if (elec_num <= 0) return QMCKL_INVALID_ARG_2;
+  if (nucl_num <= 0) return QMCKL_INVALID_ARG_3;
+  if (type_nucl_num <= 0) return QMCKL_INVALID_ARG_4;
+  if (type_nucl_vector == NULL) return QMCKL_INVALID_ARG_5;
+  if (rescale_factor_en == NULL) return QMCKL_INVALID_ARG_6;
+  if (walk_num <= 0) return QMCKL_INVALID_ARG_7;
+  if (elec_coord == NULL) return QMCKL_INVALID_ARG_8;
+  if (nucl_coord == NULL) return QMCKL_INVALID_ARG_9;
+  if (en_distance_rescaled_gl == NULL) return QMCKL_INVALID_ARG_10;
+
+  const int64_t sze = elec_num*walk_num;
+
+  qmckl_exit_code result = QMCKL_SUCCESS;
+  #pragma omp parallel
+  {
+    qmckl_exit_code rc = QMCKL_SUCCESS;
+    #pragma omp for
+    for (int64_t k=0 ; k<walk_num ; ++k)
+    {
+      for (int64_t a=0 ; a<nucl_num ; ++a) {
+        const double coord[3] = { nucl_coord[a], nucl_coord[a+nucl_num], nucl_coord[a+2*nucl_num] };
+        rc |= qmckl_distance_rescaled_gl(context, 'T', 'T', elec_num, 1,
+                                         &(elec_coord[k*elec_num]), sze,
+                                         coord, 1,
+                                         &(en_distance_rescaled_gl[4*elec_num*(a+nucl_num*k)]), elec_num,
+                                         rescale_factor_en[type_nucl_vector[a]]);
+      }
+    }
+    #pragma omp critical
+    result |= rc;
+  }
+  return result;
+}
+     #+end_src
+
+     #+begin_src c :tangle (eval c) :comments org :exports none
+qmckl_exit_code qmckl_compute_en_distance_rescaled_gl (
+          const qmckl_context context,
+          const int64_t elec_num,
+          const int64_t nucl_num,
+          const int64_t type_nucl_num,
+          int64_t* const type_nucl_vector,
+          const double*  rescale_factor_en,
+          const int64_t walk_num,
+          const double* elec_coord,
+          const double* nucl_coord,
+          double* const en_distance_rescaled_gl )
+{
+#ifdef HAVE_HPC
+  return qmckl_compute_en_distance_rescaled_gl_hpc
+#else
+  return qmckl_compute_en_distance_rescaled_gl_doc
+#endif
+    (context, elec_num, nucl_num, type_nucl_num, type_nucl_vector, rescale_factor_en,
+     walk_num, elec_coord, nucl_coord, en_distance_rescaled_gl );
+}
+     #+end_src
+     #+CALL: generate_c_interface(table=qmckl_en_distance_rescaled_gl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_en_distance_rescaled_gl_doc")
+
+     #+RESULTS:
+     #+begin_src f90 :tangle (eval f) :comments org :exports none
+     integer(c_int32_t) function qmckl_compute_en_distance_rescaled_gl_doc &
+         (context, &
+          elec_num, &
+          nucl_num, &
+          type_nucl_num, &
+          type_nucl_vector, &
+          rescale_factor_en, &
+          walk_num, &
+          elec_coord, &
+          nucl_coord, &
+          en_distance_rescaled_gl) &
+         bind(C) result(info)
+
+       use, intrinsic :: iso_c_binding
+       implicit none
+
+       integer (c_int64_t) , intent(in)  , value :: context
+       integer (c_int64_t) , intent(in)  , value :: elec_num
+       integer (c_int64_t) , intent(in)  , value :: nucl_num
+       integer (c_int64_t) , intent(in)  , value :: type_nucl_num
+       integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
+       real    (c_double ) , intent(in)          :: rescale_factor_en(nucl_num)
+       integer (c_int64_t) , intent(in)  , value :: walk_num
+       real    (c_double ) , intent(in)          :: elec_coord(elec_num,walk_num,3)
+       real    (c_double ) , intent(in)          :: nucl_coord(elec_num,3)
+       real    (c_double ) , intent(out)         :: en_distance_rescaled_gl(4,elec_num,nucl_num,walk_num)
+
+       integer(c_int32_t), external :: qmckl_compute_en_distance_rescaled_gl_doc_f
+       info = qmckl_compute_en_distance_rescaled_gl_doc_f &
+              (context, &
+          elec_num, &
+          nucl_num, &
+          type_nucl_num, &
+          type_nucl_vector, &
+          rescale_factor_en, &
+          walk_num, &
+          elec_coord, &
+          nucl_coord, &
+          en_distance_rescaled_gl)
+
+     end function qmckl_compute_en_distance_rescaled_gl_doc
+     #+end_src
+
+**** Test
+
+      #+begin_src c :tangle (eval c_test)
+assert(qmckl_electron_provided(context));
+
+{
+  printf("en_distance_rescaled_gl\n");
+  double fd[walk_num][nucl_num][elec_num][4];
+
+  double delta_x = 0.0001;
+
+  // Finite difference coefficients for gradients
+  double coef[9] = { 1.0/280.0, -4.0/105.0,  1.0/5.0, -4.0/5.0,         0.0, 4.0/5.0, -1.0/5.0, 4.0/105.0, -1.0/280.0 };
+
+  // Finite difference coefficients for Laplacian
+  double coef2[9]= {-1.0/560.0,	 8.0/315.0, -1.0/5.0,  8.0/5.0, -205.0/72.0, 8.0/5.0, -1.0/5.0,	8.0/315.0, -1.0/560.0 };
+
+  qmckl_exit_code rc;
+
+  int64_t walk_num;
+  rc = qmckl_get_electron_walk_num(context, &walk_num);
+  if (rc != QMCKL_SUCCESS) {
+    return rc;
+  }
+
+  int64_t elec_num;
+  rc = qmckl_get_electron_num(context, &elec_num);
+  if (rc != QMCKL_SUCCESS) {
+    return rc;
+  }
+
+  double elec_coord[walk_num][elec_num][3];
+  rc = qmckl_get_electron_coord (context, 'N', &(elec_coord[0][0][0]), 3*walk_num*elec_num);
+
+  double temp_coord[walk_num][elec_num][3];
+  memcpy(&(temp_coord[0][0][0]), &(elec_coord[0][0][0]), sizeof(temp_coord));
+
+  double function_values[walk_num][nucl_num][elec_num];
+
+  memset(&(fd[0][0][0][0]), 0, sizeof(fd));
+
+  for (int64_t i = 0; i < elec_num; i++) {
+    for (int64_t k = 0; k < 3; k++) {
+      for (int64_t m = -4; m <= 4; m++) {          // Apply finite difference displacement
+          
+        for (int64_t nw=0 ; nw<walk_num ; nw++) {
+          temp_coord[nw][i][k] = elec_coord[nw][i][k] + (double) m * delta_x;
+        }
+
+        // Update coordinates in the context
+        rc = qmckl_set_electron_coord (context, 'N', walk_num, &(temp_coord[0][0][0]), walk_num*3*elec_num);
+        assert(rc == QMCKL_SUCCESS);
+
+        // Call the provided function
+        rc = qmckl_get_jastrow_champ_en_distance_rescaled(context, &(function_values[0][0][0]), nucl_num*elec_num*walk_num);
+        assert(rc == QMCKL_SUCCESS);
+
+        // Accumulate derivative using finite-difference coefficients
+        for (int64_t nw=0 ; nw<walk_num ; nw++) {
+          for (int64_t j = 0; j < nucl_num; j++) {
+            fd[nw][j][i][k] += coef [m + 4] * function_values[nw][j][i];
+            fd[nw][j][i][3] += coef2[m + 4] * function_values[nw][j][i];
+          }
+        }
+      }
+      for (int64_t nw=0 ; nw<walk_num ; nw++) {
+        temp_coord[nw][i][k] = elec_coord[nw][i][k];
+      }
+    }
+  }
+
+  // Reset coordinates in the context
+  rc = qmckl_set_electron_coord (context, 'N', walk_num, &(elec_coord[0][0][0]), walk_num*3*elec_num);
+  assert(rc == QMCKL_SUCCESS);
+
+  // Normalize by the step size
+  for (int64_t nw=0 ; nw<walk_num ; nw++) {
+    for (int64_t i = 0; i < nucl_num; i++) {
+      for (int64_t k = 0; k < 4; k++) {
+        for (int64_t j = 0; j < elec_num; j++) {
+          fd[nw][i][j][k] /= delta_x;
+        }
+      }
+      for (int64_t j = 0; j < elec_num; j++) {
+        fd[nw][i][j][3] /= delta_x;
+      }
+    }
+  }
+  
+
+  double en_distance_rescaled_gl[walk_num][nucl_num][elec_num][4];
+  rc = qmckl_check(context,
+                   qmckl_get_jastrow_champ_en_distance_rescaled_gl(context,
+                                                                   &(en_distance_rescaled_gl[0][0][0][0]),
+                                                                   walk_num*nucl_num*4*elec_num)
+                   );
+
+  assert(rc == QMCKL_SUCCESS);
+
+  for (int nw = 0; nw < walk_num; nw++){
+    for (int i = 0; i < nucl_num; i++) {
+      for (int j = 0; j < elec_num; j++) {
+        for (int k = 0; k < 3; k++){
+          printf("%.10f\t", fd[nw][i][j][k]);
+          printf("%.10f\n", en_distance_rescaled_gl[nw][i][j][k]);
+          assert(fabs(fd[nw][i][j][k] - en_distance_rescaled_gl[nw][i][j][k]) < 1.e-8);
+        }
+        int k=3;
+        if (i != j) {
+          printf("%.10f\t", fd[nw][i][j][k]);
+          printf("%.10f\n", en_distance_rescaled_gl[nw][i][j][k]);
+          assert(fabs(fd[nw][i][j][k] - en_distance_rescaled_gl[nw][i][j][k]) < 1.e-6);
+        }
+      }
+    }
+  }
+  printf("OK\n");
+}
+
+      #+end_src
+
 *** Electron-nucleus component
 
    Calculate the electron-electron jastrow component ~factor_en~ using the ~a_vector~
@@ -4299,13 +5300,22 @@ qmckl_get_jastrow_champ_factor_en(qmckl_context context,
   rc = qmckl_provide_jastrow_champ_factor_en(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  int64_t sze=ctx->electron.walker.num;
+  if (factor_en == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_factor_en",
+                           "Null pointer");
+  }
+  
+  const int64_t sze=ctx->electron.walker.num;
+
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
                            "qmckl_get_jastrow_champ_factor_en",
                            "Array too small. Expected walk_num");
   }
+
   memcpy(factor_en, ctx->jastrow_champ.factor_en, sze*sizeof(double));
 
   return QMCKL_SUCCESS;
@@ -4721,13 +5731,22 @@ qmckl_get_jastrow_champ_factor_en_gl(qmckl_context context,
   qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
-  int64_t sze = ctx->electron.walker.num * 4 * ctx->electron.num;
+  if (factor_en_gl == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_factor_en_gl",
+                           "Null pointer");
+  }
+
+  const int64_t sze = ctx->electron.walker.num * 4 * ctx->electron.num;
+
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
                            "qmckl_get_jastrow_champ_factor_en_gl",
                            "Array too small. Expected 4*walk_num*elec_num");
   }
+
   memcpy(factor_en_gl, ctx->jastrow_champ.factor_en_gl, sze*sizeof(double));
 
   return QMCKL_SUCCESS;
@@ -5227,18 +6246,6 @@ print("factor_en_gl[3][0]:",factor_en_gl[3][0])
 /* Check if Jastrow is properly initialized */
 assert(qmckl_jastrow_champ_provided(context));
 
-/*
-// calculate factor_en_gl
-
-// check factor_en_gl
-assert(fabs( 0.19656663796630847 - factor_en_gl[0][0][0]) < 1.e-12);
-assert(fabs( -0.3945140890522283 - factor_en_gl[0][1][0]) < 1.e-12);
-assert(fabs( 0.5082964671286118  - factor_en_gl[0][2][0]) < 1.e-12);
-assert(fabs( -1.8409460670666289 - factor_en_gl[0][3][0]) < 1.e-12);
-,*/
-
-
-//----
 {
   printf("factor_en_gl\n");
   double fd[walk_num][4][elec_num];
@@ -5346,767 +6353,6 @@ assert(fabs( -1.8409460670666289 - factor_en_gl[0][3][0]) < 1.e-12);
 
       #+end_src
 
-*** Electron-nucleus rescaled distances
-
-    ~en_distance_rescaled~ stores the matrix of the rescaled distances between
-    electrons and nuclei.
-
-    \[
-    C_{i\alpha} = \frac{ 1 - e^{-\kappa_\alpha R_{i\alpha}}}{\kappa_\alpha}
-    \]
-
-   where \(R_{i\alpha}\) is the matrix of electron-nucleus distances.
-
-**** Get
-
-     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
-qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled);
-     #+end_src
-
-
-     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled)
-{
-
-  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
-    return QMCKL_NULL_CONTEXT;
-  }
-
-  qmckl_exit_code rc;
-
-  rc = qmckl_provide_en_distance_rescaled(context);
-  if (rc != QMCKL_SUCCESS) return rc;
-
-  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
-  assert (ctx != NULL);
-
-  size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num;
-  memcpy(distance_rescaled, ctx->jastrow_champ.en_distance_rescaled, sze * sizeof(double));
-
-  return QMCKL_SUCCESS;
-}
-     #+end_src
-
-**** Provide                                                       :noexport:
-
-     #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
-qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context);
-     #+end_src
-
-     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context)
-{
-
-  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
-    return QMCKL_NULL_CONTEXT;
-  }
-
-  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
-  assert (ctx != NULL);
-
-  if (!(ctx->nucleus.provided)) {
-    return QMCKL_NOT_PROVIDED;
-  }
-
-  /* Compute if necessary */
-  if (ctx->electron.walker.point.date > ctx->jastrow_champ.en_distance_rescaled_date) {
-
-    if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
-      if (ctx->jastrow_champ.en_distance_rescaled != NULL) {
-        qmckl_exit_code rc = qmckl_free(context, ctx->jastrow_champ.en_distance_rescaled);
-        if (rc != QMCKL_SUCCESS) {
-          return qmckl_failwith( context, rc,
-                                 "qmckl_provide_en_distance_rescaled",
-                                 "Unable to free ctx->jastrow_champ.en_distance_rescaled");
-        }
-        ctx->jastrow_champ.en_distance_rescaled = NULL;
-      }
-    }
-
-    /* Allocate array */
-    if (ctx->jastrow_champ.en_distance_rescaled == NULL) {
-
-      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
-      mem_info.size = ctx->electron.num * ctx->nucleus.num *
-        ctx->electron.walker.num * sizeof(double);
-      double* en_distance_rescaled = (double*) qmckl_malloc(context, mem_info);
-
-      if (en_distance_rescaled == NULL) {
-        return qmckl_failwith( context,
-                               QMCKL_ALLOCATION_FAILED,
-                               "qmckl_provide_en_distance_rescaled",
-                               NULL);
-      }
-      ctx->jastrow_champ.en_distance_rescaled = en_distance_rescaled;
-    }
-
-    qmckl_exit_code rc =
-      qmckl_compute_en_distance_rescaled(context,
-                                ctx->electron.num,
-                                ctx->nucleus.num,
-                                ctx->jastrow_champ.type_nucl_num,
-                                ctx->jastrow_champ.type_nucl_vector,
-                                ctx->jastrow_champ.rescale_factor_en,
-                                ctx->electron.walker.num,
-                                ctx->electron.walker.point.coord.data,
-                                ctx->nucleus.coord.data,
-                                ctx->jastrow_champ.en_distance_rescaled);
-
-    if (rc != QMCKL_SUCCESS) {
-      return rc;
-    }
-
-    ctx->jastrow_champ.en_distance_rescaled_date = ctx->date;
-  }
-
-  return QMCKL_SUCCESS;
-}
-     #+end_src
-
-**** Compute
-     :PROPERTIES:
-     :Name:     qmckl_compute_en_distance_rescaled
-     :CRetType: qmckl_exit_code
-     :FRetType: qmckl_exit_code
-     :END:
-
-    #+NAME: qmckl_en_distance_rescaled_args
-     | Variable               | Type                                   | In/Out | Description                       |
-     |------------------------+----------------------------------------+--------+-----------------------------------|
-     | ~context~              | ~qmckl_context~                        | in     | Global state                      |
-     | ~elec_num~             | ~int64_t~                              | in     | Number of electrons               |
-     | ~nucl_num~             | ~int64_t~                              | in     | Number of nuclei                  |
-     | ~type_nucl_num~        | ~int64_t~                              | in     | Number of types of nuclei         |
-     | ~type_nucl_vector~     | ~int64_t[nucl_num]~                    | in     | Number of types of nuclei         |
-     | ~rescale_factor_en~    | ~double[type_nucl_num]~                | in     | The factor for rescaled distances |
-     | ~walk_num~             | ~int64_t~                              | in     | Number of walkers                 |
-     | ~elec_coord~           | ~double[3][walk_num][elec_num]~        | in     | Electron coordinates              |
-     | ~nucl_coord~           | ~double[3][elec_num]~                  | in     | Nuclear coordinates               |
-     | ~en_distance_rescaled~ | ~double[walk_num][nucl_num][elec_num]~ | out    | Electron-nucleus distances        |
-
-     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-function qmckl_compute_en_distance_rescaled_doc(context, &
-     elec_num, nucl_num, type_nucl_num, &
-     type_nucl_vector, rescale_factor_en, walk_num, elec_coord, &
-     nucl_coord, en_distance_rescaled) &
-     bind(C) result(info)
-  use qmckl
-  implicit none
-  integer (qmckl_context), intent(in), value :: context
-  integer (c_int64_t) , intent(in)  , value :: elec_num
-  integer (c_int64_t) , intent(in)  , value :: nucl_num
-  integer (c_int64_t) , intent(in)  , value :: type_nucl_num
-  integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
-  real    (c_double ) , intent(in)          :: rescale_factor_en(type_nucl_num)
-  integer (c_int64_t) , intent(in)  , value :: walk_num
-  real    (c_double ) , intent(in)          :: elec_coord(elec_num,walk_num,3)
-  real    (c_double ) , intent(in)          :: nucl_coord(nucl_num,3)
-  real    (c_double ) , intent(out)         :: en_distance_rescaled(elec_num,nucl_num,walk_num)
-  integer(qmckl_exit_code)                   :: info
-
-  integer*8 :: i, k
-  double precision      :: coord(3)
-
-  info = QMCKL_SUCCESS
-
-  if (context == QMCKL_NULL_CONTEXT) then
-     info = QMCKL_INVALID_CONTEXT
-     return
-  endif
-
-  if (elec_num <= 0) then
-     info = QMCKL_INVALID_ARG_2
-     return
-  endif
-
-  if (nucl_num <= 0) then
-     info = QMCKL_INVALID_ARG_3
-     return
-  endif
-
-  if (walk_num <= 0) then
-     info = QMCKL_INVALID_ARG_5
-     return
-  endif
-
-  do i=1, nucl_num
-     coord(1:3) = nucl_coord(i,1:3)
-     do k=1,walk_num
-        info = qmckl_distance_rescaled(context, 'T', 'N', elec_num, 1_8, &
-             elec_coord(1,k,1), elec_num*walk_num, coord, 3_8, &
-             en_distance_rescaled(1,i,k), elec_num, rescale_factor_en(type_nucl_vector(i)+1))
-        if (info /= QMCKL_SUCCESS) then
-           return
-        endif
-     end do
-  end do
-
-end function qmckl_compute_en_distance_rescaled_doc
-     #+end_src
-
-     #+begin_src c :tangle (eval c) :comments org :exports none
-qmckl_exit_code qmckl_compute_en_distance_rescaled_hpc (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const int64_t nucl_num,
-          const int64_t type_nucl_num,
-          const int64_t* type_nucl_vector,
-          const double* rescale_factor_en,
-          const int64_t walk_num,
-          const double* elec_coord,
-          const double* nucl_coord,
-          double* const en_distance_rescaled )
-{
-  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT;
-  if (elec_num <= 0) return QMCKL_INVALID_ARG_2;
-  if (nucl_num <= 0) return QMCKL_INVALID_ARG_3;
-  if (type_nucl_num <= 0) return QMCKL_INVALID_ARG_4;
-  if (type_nucl_vector == NULL) return QMCKL_INVALID_ARG_5;
-  if (rescale_factor_en == NULL) return QMCKL_INVALID_ARG_6;
-  if (walk_num <= 0) return QMCKL_INVALID_ARG_7;
-  if (elec_coord == NULL) return QMCKL_INVALID_ARG_8;
-  if (nucl_coord == NULL) return QMCKL_INVALID_ARG_9;
-  if (en_distance_rescaled == NULL) return QMCKL_INVALID_ARG_10;
-
-  int64_t sze = elec_num*walk_num;
-
-  qmckl_exit_code result = QMCKL_SUCCESS;
-  #pragma omp parallel
-  {
-    qmckl_exit_code rc = QMCKL_SUCCESS;
-    #pragma omp for
-    for (int64_t k=0 ; k<walk_num ; ++k)
-    {
-      for (int64_t a=0 ; a<nucl_num ; ++a) {
-        const double coord[3] = { nucl_coord[a], nucl_coord[a+nucl_num], nucl_coord[a+2*nucl_num] };
-        rc |= qmckl_distance_rescaled(context, 'T', 'N', elec_num, 1,
-                                      &(elec_coord[k*elec_num]), sze,
-                                      coord, 3,
-                                      &(en_distance_rescaled[elec_num*(a+nucl_num*k)]), elec_num,
-                                      rescale_factor_en[type_nucl_vector[a]]);
-      }
-    }
-    #pragma omp critical
-    result |= rc;
-  }
-  return result;
-}
-     #+end_src
-
-     #+begin_src c :tangle (eval h_private_func) :comments org :exports none
-qmckl_exit_code qmckl_compute_en_distance_rescaled_doc (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const int64_t nucl_num,
-          const int64_t type_nucl_num,
-          const int64_t* type_nucl_vector,
-          const double*  rescale_factor_en,
-          const int64_t walk_num,
-          const double* elec_coord,
-          const double* nucl_coord,
-          double* const en_distance_rescaled );
-
-qmckl_exit_code qmckl_compute_en_distance_rescaled_hpc (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const int64_t nucl_num,
-          const int64_t type_nucl_num,
-          const int64_t* type_nucl_vector,
-          const double* rescale_factor_en,
-          const int64_t walk_num,
-          const double* elec_coord,
-          const double* nucl_coord,
-          double* const en_distance_rescaled );
-
-qmckl_exit_code qmckl_compute_en_distance_rescaled (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const int64_t nucl_num,
-          const int64_t type_nucl_num,
-          const int64_t* type_nucl_vector,
-          const double*  rescale_factor_en,
-          const int64_t walk_num,
-          const double* elec_coord,
-          const double* nucl_coord,
-          double* const en_distance_rescaled );
-     #+end_src
-
-     #+begin_src c :tangle (eval c) :comments org :exports none
-qmckl_exit_code qmckl_compute_en_distance_rescaled (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const int64_t nucl_num,
-          const int64_t type_nucl_num,
-          const int64_t* type_nucl_vector,
-          const double* rescale_factor_en,
-          const int64_t walk_num,
-          const double* elec_coord,
-          const double* nucl_coord,
-          double* const en_distance_rescaled )
-{
-#ifdef HAVE_HPC
-  return qmckl_compute_en_distance_rescaled_hpc
-#else
-  return qmckl_compute_en_distance_rescaled_doc
-#endif
-  (context, elec_num, nucl_num, type_nucl_num, type_nucl_vector,
-  rescale_factor_en, walk_num, elec_coord, nucl_coord, en_distance_rescaled );
-}
-    #+end_src
-
-**** Test
-
-     #+begin_src python :results output :exports none
-import numpy as np
-
-kappa = 0.6
-
-elec_1_w1 = np.array( [-0.250655104764153,  0.503070975550133      ,  -0.166554344502303])
-elec_2_w1 = np.array( [-0.587812193472177, -0.128751981129274      ,   0.187773606533075])
-elec_5_w1 = np.array( [-0.127732483187947, -0.138975497694196      ,  -8.669850480215846E-002])
-elec_6_w1 = np.array( [-0.232271834949124, -1.059321673434182E-002 ,  -0.504862241464867])
-nucl_1    = np.array( [ 0., 0., 0. ])
-nucl_2    = np.array( [ 0., 0., 2.059801 ])
-
-print ( "[0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-nucl_1)) )/kappa )
-print ( "[1][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-nucl_2)) )/kappa )
-print ( "[0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w1-nucl_1)) )/kappa )
-print ( "[0][5] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_5_w1-nucl_1)) )/kappa )
-print ( "[1][5] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_5_w1-nucl_2)) )/kappa )
-print ( "[0][6] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_6_w1-nucl_1)) )/kappa )
-
-     #+end_src
-
-     #+RESULTS:
-     : [0][0] :  0.4942158656729477
-     : [1][0] :  1.2464137498005765
-     : [0][1] :  0.5248654474756858
-     : [0][5] :  0.19529459944794733
-     : [1][5] :  1.2091967687767369
-     : [0][6] :  0.4726452953409436
-
-
-      #+begin_src c :tangle (eval c_test)
-
-assert(qmckl_electron_provided(context));
-assert(qmckl_nucleus_provided(context));
-
-double en_distance_rescaled[walk_num][nucl_num][elec_num];
-
-rc = qmckl_check(context,
-                 qmckl_get_electron_en_distance_rescaled(context, &(en_distance_rescaled[0][0][0]))
-                 );
-assert (rc == QMCKL_SUCCESS);
-
-// (e,n,w) in Fortran notation
-// (1,1,1)
-assert(fabs(en_distance_rescaled[0][0][0] - 0.4942158656729477) < 1.e-12);
-// (1,2,1)
-assert(fabs(en_distance_rescaled[0][1][0] - 1.2464137498005765) < 1.e-12);
-// (2,1,1)
-assert(fabs(en_distance_rescaled[0][0][1] - 0.5248654474756858) < 1.e-12);
-// (1,1,2)
-assert(fabs(en_distance_rescaled[0][0][5] - 0.19529459944794733) < 1.e-12);
-// (1,2,2)
-assert(fabs(en_distance_rescaled[0][1][5] - 1.2091967687767369) < 1.e-12);
-// (2,1,2)
-assert(fabs(en_distance_rescaled[0][0][6] - 0.4726452953409436) < 1.e-12);
-
-
-      #+end_src
-
-*** Electron-electron rescaled distance gradients and Laplacian with respect to electron coordinates
-
-   The rescaled distances, represented by $C_{i\alpha} = (1 - e^{-\kappa_\alpha R_{i\alpha}})/\kappa$
-   are differentiated with respect to the electron coordinates.
-   This information is stored in the tensor
-   ~en_distance_rescaled_gl~. The initial three sequential
-   elements of this three-index tensor provide the $x$, $y$, and $z$
-   direction derivatives, while the fourth index corresponds to the Laplacian.
-
-**** Get
-
-     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
-qmckl_exit_code qmckl_get_electron_en_distance_rescaled_gl(qmckl_context context, double* distance_rescaled_gl);
-     #+end_src
-
-     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_get_electron_en_distance_rescaled_gl(qmckl_context context, double* distance_rescaled_gl)
-{
-
-  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
-    return QMCKL_NULL_CONTEXT;
-  }
-
-  qmckl_exit_code rc;
-
-  rc = qmckl_provide_en_distance_rescaled_gl(context);
-  if (rc != QMCKL_SUCCESS) return rc;
-
-  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
-  assert (ctx != NULL);
-
-  size_t sze = 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num;
-  memcpy(distance_rescaled_gl, ctx->jastrow_champ.en_distance_rescaled_gl, sze * sizeof(double));
-
-  return QMCKL_SUCCESS;
-}
-     #+end_src
-
-**** Provide                                                       :noexport:
-
-     #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
-qmckl_exit_code qmckl_provide_en_distance_rescaled_gl(qmckl_context context);
-     #+end_src
-
-     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_provide_en_distance_rescaled_gl(qmckl_context context)
-{
-
-  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
-    return QMCKL_NULL_CONTEXT;
-  }
-
-  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
-  assert (ctx != NULL);
-
-  if (!(ctx->nucleus.provided)) {
-    return QMCKL_NOT_PROVIDED;
-  }
-
-  /* Compute if necessary */
-  if (ctx->electron.walker.point.date > ctx->jastrow_champ.en_distance_rescaled_gl_date) {
-
-    if (ctx->electron.walker.num > ctx->electron.walker_old.num) {
-      if (ctx->jastrow_champ.en_distance_rescaled_gl != NULL) {
-        qmckl_exit_code rc = qmckl_free(context, ctx->jastrow_champ.en_distance_rescaled_gl);
-        if (rc != QMCKL_SUCCESS) {
-          return qmckl_failwith( context, rc,
-                                 "qmckl_provide_en_distance_rescaled_gl",
-                                 "Unable to free ctx->jastrow_champ.en_distance_rescaled_gl");
-        }
-        ctx->jastrow_champ.en_distance_rescaled_gl = NULL;
-      }
-    }
-
-    /* Allocate array */
-    if (ctx->jastrow_champ.en_distance_rescaled_gl == NULL) {
-
-      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
-      mem_info.size = 4 * ctx->electron.num * ctx->nucleus.num *
-        ctx->electron.walker.num * sizeof(double);
-      double* en_distance_rescaled_gl = (double*) qmckl_malloc(context, mem_info);
-
-      if (en_distance_rescaled_gl == NULL) {
-        return qmckl_failwith( context,
-                               QMCKL_ALLOCATION_FAILED,
-                               "qmckl_provide_en_distance_rescaled_gl",
-                               NULL);
-      }
-      ctx->jastrow_champ.en_distance_rescaled_gl = en_distance_rescaled_gl;
-    }
-
-    qmckl_exit_code rc =
-      qmckl_compute_en_distance_rescaled_gl(context,
-                                ctx->electron.num,
-                                ctx->nucleus.num,
-                                ctx->jastrow_champ.type_nucl_num,
-                                ctx->jastrow_champ.type_nucl_vector,
-                                ctx->jastrow_champ.rescale_factor_en,
-                                ctx->electron.walker.num,
-                                ctx->electron.walker.point.coord.data,
-                                ctx->nucleus.coord.data,
-                                ctx->jastrow_champ.en_distance_rescaled_gl);
-    if (rc != QMCKL_SUCCESS) {
-      return rc;
-    }
-
-    ctx->jastrow_champ.en_distance_rescaled_gl_date = ctx->date;
-  }
-
-  return QMCKL_SUCCESS;
-}
-     #+end_src
-
-**** Compute
-     :PROPERTIES:
-     :Name:     qmckl_compute_en_distance_rescaled_gl
-     :CRetType: qmckl_exit_code
-     :FRetType: qmckl_exit_code
-     :END:
-
-     #+NAME: qmckl_en_distance_rescaled_gl_args
-     | Variable                       | Type                                      | In/Out | Description                           |
-     |--------------------------------+-------------------------------------------+--------+---------------------------------------|
-     | ~context~                      | ~qmckl_context~                           | in     | Global state                          |
-     | ~elec_num~                     | ~int64_t~                                 | in     | Number of electrons                   |
-     | ~nucl_num~                     | ~int64_t~                                 | in     | Number of nuclei                      |
-     | ~type_nucl_num~                | ~int64_t~                                 | in     | Number of nucleus types               |
-     | ~type_nucl_vector~             | ~int64_t[nucl_num]~                       | in     | Array of nucleus types                |
-     | ~rescale_factor_en~            | ~double[nucl_num]~                        | in     | The factors for rescaled distances    |
-     | ~walk_num~                     | ~int64_t~                                 | in     | Number of walkers                     |
-     | ~elec_coord~                   | ~double[3][walk_num][elec_num]~           | in     | Electron coordinates                  |
-     | ~nucl_coord~                   | ~double[3][elec_num]~                     | in     | Nuclear coordinates                   |
-     | ~en_distance_rescaled_gl~ | ~double[walk_num][nucl_num][elec_num][4]~ | out    | Electron-nucleus distance derivatives |
-
-     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_en_distance_rescaled_gl_doc_f(context, elec_num, nucl_num, &
-     type_nucl_num, type_nucl_vector, rescale_factor_en, walk_num, elec_coord, &
-     nucl_coord, en_distance_rescaled_gl) &
-     result(info)
-  use qmckl
-  implicit none
-  integer(qmckl_context), intent(in)  :: context
-  integer*8             , intent(in)  :: elec_num
-  integer*8             , intent(in)  :: nucl_num
-  integer*8             , intent(in)  :: type_nucl_num
-  integer*8             , intent(in)  :: type_nucl_vector(nucl_num)
-  double precision      , intent(in)  :: rescale_factor_en(nucl_num)
-  integer*8             , intent(in)  :: walk_num
-  double precision      , intent(in)  :: elec_coord(elec_num,walk_num,3)
-  double precision      , intent(in)  :: nucl_coord(nucl_num,3)
-  double precision      , intent(out) :: en_distance_rescaled_gl(4,elec_num,nucl_num,walk_num)
-
-  integer*8 :: i, k
-  double precision :: coord(3)
-
-  info = QMCKL_SUCCESS
-
-  if (context == QMCKL_NULL_CONTEXT) then
-     info = QMCKL_INVALID_CONTEXT
-     return
-  endif
-
-  if (elec_num <= 0) then
-     info = QMCKL_INVALID_ARG_2
-     return
-  endif
-
-  if (nucl_num <= 0) then
-     info = QMCKL_INVALID_ARG_3
-     return
-  endif
-
-  if (walk_num <= 0) then
-     info = QMCKL_INVALID_ARG_5
-     return
-  endif
-
-  do i=1, nucl_num
-     coord(1:3) = nucl_coord(i,1:3)
-     do k=1,walk_num
-        info = qmckl_distance_rescaled_gl(context, 'T', 'T', elec_num, 1_8, &
-             elec_coord(1,k,1), elec_num*walk_num, coord, 1_8, &
-             en_distance_rescaled_gl(1,1,i,k), elec_num, rescale_factor_en(type_nucl_vector(i)+1))
-        if (info /= QMCKL_SUCCESS) then
-           return
-        endif
-     end do
-  end do
-
-end function qmckl_compute_en_distance_rescaled_gl_doc_f
-     #+end_src
-
-     #+begin_src c :tangle (eval h_private_func) :comments org :exports none
-qmckl_exit_code qmckl_compute_en_distance_rescaled_gl_doc (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const int64_t nucl_num,
-          const int64_t type_nucl_num,
-          int64_t* const type_nucl_vector,
-          const double*  rescale_factor_en,
-          const int64_t walk_num,
-          const double* elec_coord,
-          const double* nucl_coord,
-          double* const en_distance_rescaled_gl );
-
-qmckl_exit_code qmckl_compute_en_distance_rescaled_gl_hpc (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const int64_t nucl_num,
-          const int64_t type_nucl_num,
-          int64_t* const type_nucl_vector,
-          const double*  rescale_factor_en,
-          const int64_t walk_num,
-          const double* elec_coord,
-          const double* nucl_coord,
-          double* const en_distance_rescaled_gl );
-
-qmckl_exit_code qmckl_compute_en_distance_rescaled_gl (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const int64_t nucl_num,
-          const int64_t type_nucl_num,
-          int64_t* const type_nucl_vector,
-          const double*  rescale_factor_en,
-          const int64_t walk_num,
-          const double* elec_coord,
-          const double* nucl_coord,
-          double* const en_distance_rescaled_gl );
-     #+end_src
-
-     #+begin_src c :tangle (eval c) :comments org :exports none
-qmckl_exit_code qmckl_compute_en_distance_rescaled_gl_hpc (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const int64_t nucl_num,
-          const int64_t type_nucl_num,
-          int64_t* const type_nucl_vector,
-          const double*  rescale_factor_en,
-          const int64_t walk_num,
-          const double* elec_coord,
-          const double* nucl_coord,
-          double* const en_distance_rescaled_gl )
-{
-  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT;
-  if (elec_num <= 0) return QMCKL_INVALID_ARG_2;
-  if (nucl_num <= 0) return QMCKL_INVALID_ARG_3;
-  if (type_nucl_num <= 0) return QMCKL_INVALID_ARG_4;
-  if (type_nucl_vector == NULL) return QMCKL_INVALID_ARG_5;
-  if (rescale_factor_en == NULL) return QMCKL_INVALID_ARG_6;
-  if (walk_num <= 0) return QMCKL_INVALID_ARG_7;
-  if (elec_coord == NULL) return QMCKL_INVALID_ARG_8;
-  if (nucl_coord == NULL) return QMCKL_INVALID_ARG_9;
-  if (en_distance_rescaled_gl == NULL) return QMCKL_INVALID_ARG_10;
-
-  int64_t sze = elec_num*walk_num;
-
-  qmckl_exit_code result = QMCKL_SUCCESS;
-  #pragma omp parallel
-  {
-    qmckl_exit_code rc = QMCKL_SUCCESS;
-    #pragma omp for
-    for (int64_t k=0 ; k<walk_num ; ++k)
-    {
-      for (int64_t a=0 ; a<nucl_num ; ++a) {
-        const double coord[3] = { nucl_coord[a], nucl_coord[a+nucl_num], nucl_coord[a+2*nucl_num] };
-        rc |= qmckl_distance_rescaled_gl(context, 'T', 'T', elec_num, 1,
-                                         &(elec_coord[k*elec_num]), sze,
-                                         coord, 1,
-                                         &(en_distance_rescaled_gl[4*elec_num*(a+nucl_num*k)]), elec_num,
-                                         rescale_factor_en[type_nucl_vector[a]]);
-      }
-    }
-    #pragma omp critical
-    result |= rc;
-  }
-  return result;
-}
-     #+end_src
-
-     #+begin_src c :tangle (eval c) :comments org :exports none
-qmckl_exit_code qmckl_compute_en_distance_rescaled_gl (
-          const qmckl_context context,
-          const int64_t elec_num,
-          const int64_t nucl_num,
-          const int64_t type_nucl_num,
-          int64_t* const type_nucl_vector,
-          const double*  rescale_factor_en,
-          const int64_t walk_num,
-          const double* elec_coord,
-          const double* nucl_coord,
-          double* const en_distance_rescaled_gl )
-{
-#ifdef HAVE_HPC
-  return qmckl_compute_en_distance_rescaled_gl_hpc
-#else
-  return qmckl_compute_en_distance_rescaled_gl_doc
-#endif
-    (context, elec_num, nucl_num, type_nucl_num, type_nucl_vector, rescale_factor_en,
-     walk_num, elec_coord, nucl_coord, en_distance_rescaled_gl );
-}
-     #+end_src
-     #+CALL: generate_c_interface(table=qmckl_en_distance_rescaled_gl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_en_distance_rescaled_gl_doc")
-
-     #+RESULTS:
-     #+begin_src f90 :tangle (eval f) :comments org :exports none
-     integer(c_int32_t) function qmckl_compute_en_distance_rescaled_gl_doc &
-         (context, &
-          elec_num, &
-          nucl_num, &
-          type_nucl_num, &
-          type_nucl_vector, &
-          rescale_factor_en, &
-          walk_num, &
-          elec_coord, &
-          nucl_coord, &
-          en_distance_rescaled_gl) &
-         bind(C) result(info)
-
-       use, intrinsic :: iso_c_binding
-       implicit none
-
-       integer (c_int64_t) , intent(in)  , value :: context
-       integer (c_int64_t) , intent(in)  , value :: elec_num
-       integer (c_int64_t) , intent(in)  , value :: nucl_num
-       integer (c_int64_t) , intent(in)  , value :: type_nucl_num
-       integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
-       real    (c_double ) , intent(in)          :: rescale_factor_en(nucl_num)
-       integer (c_int64_t) , intent(in)  , value :: walk_num
-       real    (c_double ) , intent(in)          :: elec_coord(elec_num,walk_num,3)
-       real    (c_double ) , intent(in)          :: nucl_coord(elec_num,3)
-       real    (c_double ) , intent(out)         :: en_distance_rescaled_gl(4,elec_num,nucl_num,walk_num)
-
-       integer(c_int32_t), external :: qmckl_compute_en_distance_rescaled_gl_doc_f
-       info = qmckl_compute_en_distance_rescaled_gl_doc_f &
-              (context, &
-          elec_num, &
-          nucl_num, &
-          type_nucl_num, &
-          type_nucl_vector, &
-          rescale_factor_en, &
-          walk_num, &
-          elec_coord, &
-          nucl_coord, &
-          en_distance_rescaled_gl)
-
-     end function qmckl_compute_en_distance_rescaled_gl_doc
-     #+end_src
-
-**** Test
-
-     #+begin_src python :results output :exports none
-import numpy as np
-
-# TODO
-     #+end_src
-
-
-      #+begin_src c :tangle (eval c_test)
-
-assert(qmckl_electron_provided(context));
-
-assert(qmckl_nucleus_provided(context));
-
-double en_distance_rescaled_gl[walk_num][4][nucl_num][elec_num];
-
-rc = qmckl_check(context,
-                 qmckl_get_electron_en_distance_rescaled_gl(context, &(en_distance_rescaled_gl[0][0][0][0]))
-                 );
-assert (rc == QMCKL_SUCCESS);
-
-// TODO: check exact values
-//// (e,n,w) in Fortran notation
-//// (1,1,1)
-//assert(fabs(en_distance_rescaled[0][0][0] - 7.546738741619978) < 1.e-12);
-//
-//// (1,2,1)
-//assert(fabs(en_distance_rescaled[0][1][0] - 8.77102435246984) < 1.e-12);
-//
-//// (2,1,1)
-//assert(fabs(en_distance_rescaled[0][0][1] - 3.698922010513608) < 1.e-12);
-//
-//// (1,1,2)
-//assert(fabs(en_distance_rescaled[1][0][0] - 5.824059436060509) < 1.e-12);
-//
-//// (1,2,2)
-//assert(fabs(en_distance_rescaled[1][1][0] - 7.080482110317645) < 1.e-12);
-//
-//// (2,1,2)
-//assert(fabs(en_distance_rescaled[1][0][1] - 3.1804527583077356) < 1.e-12);
-
-      #+end_src
-
 ** Electron-electron-nucleus component
 *** Electron-electron rescaled distances in $J_\text{eeN}$
 
@@ -6124,14 +6370,14 @@ assert (rc == QMCKL_SUCCESS);
      #+begin_src c :comments org :tangle (eval h_func) :noweb yes
 qmckl_exit_code
 qmckl_get_jastrow_champ_een_rescaled_e(qmckl_context context,
-                                 double* const distance_rescaled,
+                                 double* const een_rescaled_e,
                                  const int64_t size_max);
      #+end_src
 
      #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
 qmckl_exit_code
 qmckl_get_jastrow_champ_een_rescaled_e(qmckl_context context,
-                                 double* const distance_rescaled,
+                                 double* const een_rescaled_e,
                                  const int64_t size_max)
 {
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
@@ -6146,14 +6392,23 @@ qmckl_get_jastrow_champ_een_rescaled_e(qmckl_context context,
   qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
-  int64_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1);
+  if (een_rescaled_e == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_een_rescaled_e",
+                           "Null pointer");
+  }
+
+  const int64_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1);
+
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
-                           "qmckl_get_jastrow_champ_factor_een_rescaled_e",
+                           "qmckl_get_jastrow_champ_een_rescaled_e",
                            "Array too small. Expected elec_num*elec_num*walk_num*(cord_num + 1)");
   }
-  memcpy(distance_rescaled, ctx->jastrow_champ.een_rescaled_e, sze * sizeof(double));
+
+  memcpy(een_rescaled_e, ctx->jastrow_champ.een_rescaled_e, sze * sizeof(double));
 
   return QMCKL_SUCCESS;
 }
@@ -6655,14 +6910,14 @@ assert(fabs(een_rescaled_e[0][2][1][5]- 0.5257156022077619 ) < 1.e-12);
      #+begin_src c :comments org :tangle (eval h_func) :noweb yes
 qmckl_exit_code
 qmckl_get_jastrow_champ_een_rescaled_e_gl(qmckl_context context,
-                                         double* const distance_rescaled,
+                                         double* const een_rescaled_e_gl,
                                          const int64_t size_max);
      #+end_src
 
      #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
 qmckl_exit_code
 qmckl_get_jastrow_champ_een_rescaled_e_gl(qmckl_context context,
-                                         double* const distance_rescaled,
+                                         double* const een_rescaled_e_gl,
                                          const int64_t size_max)
 {
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
@@ -6677,14 +6932,23 @@ qmckl_get_jastrow_champ_een_rescaled_e_gl(qmckl_context context,
   qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
-  int64_t sze = ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1);
+  if (een_rescaled_e_gl == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_een_rescaled_e_gl",
+                           "Null pointer");
+  }
+
+  const int64_t sze = ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1);
+
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
                            "qmckl_get_jastrow_champ_een_rescaled_e_gl",
                            "Array too small. Expected elec_num*4*elec_num*walk_num*(cord_num + 1)");
   }
-  memcpy(distance_rescaled, ctx->jastrow_champ.een_rescaled_e_gl, sze * sizeof(double));
+
+  memcpy(een_rescaled_e_gl, ctx->jastrow_champ.een_rescaled_e_gl, sze * sizeof(double));
 
   return QMCKL_SUCCESS;
 }
@@ -7203,14 +7467,14 @@ assert(fabs(een_rescaled_e_gl[0][2][1][0][5] +  0.5416751547830984    ) < 1.e-12
      #+begin_src c :comments org :tangle (eval h_func) :noweb yes
 qmckl_exit_code
 qmckl_get_jastrow_champ_een_rescaled_n(qmckl_context context,
-                                 double* const distance_rescaled,
+                                 double* const een_rescaled_n,
                                  const int64_t size_max);
      #+end_src
 
      #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
 qmckl_exit_code
 qmckl_get_jastrow_champ_een_rescaled_n(qmckl_context context,
-                                 double* const distance_rescaled,
+                                 double* const een_rescaled_n,
                                  const int64_t size_max)
 {
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
@@ -7225,14 +7489,23 @@ qmckl_get_jastrow_champ_een_rescaled_n(qmckl_context context,
   qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
-  int64_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1);
+  if (een_rescaled_n == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_een_rescaled_n",
+                           "Null pointer");
+  }
+
+  const int64_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1);
+
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
                            "qmckl_get_jastrow_champ_een_rescaled_n",
                            "Array too small. Expected elec_num*nucl_num*walk_num*(cord_num + 1)");
   }
-  memcpy(distance_rescaled, ctx->jastrow_champ.een_rescaled_n, sze * sizeof(double));
+
+  memcpy(een_rescaled_n, ctx->jastrow_champ.een_rescaled_n, sze * sizeof(double));
 
   return QMCKL_SUCCESS;
 }
@@ -7615,14 +7888,14 @@ assert(fabs(een_rescaled_n[0][2][1][5]-0.07534033469115217)< 1.e-12);
      #+begin_src c :comments org :tangle (eval h_func) :noweb yes
 qmckl_exit_code
 qmckl_get_jastrow_champ_een_rescaled_n_gl(qmckl_context context,
-                                         double* const distance_rescaled,
+                                         double* const een_rescaled_n_gl,
                                          const int64_t size_max);
      #+end_src
 
      #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
 qmckl_exit_code
 qmckl_get_jastrow_champ_een_rescaled_n_gl(qmckl_context context,
-                                         double* const distance_rescaled,
+                                         double* const een_rescaled_n_gl,
                                          const int64_t size_max)
 {
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
@@ -7637,14 +7910,23 @@ qmckl_get_jastrow_champ_een_rescaled_n_gl(qmckl_context context,
   qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
-  int64_t sze = ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1);
+  if (een_rescaled_n_gl == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_een_rescaled_n_gl",
+                           "Null pointer");
+  }
+
+  const int64_t sze = ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1);
+
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
                            "qmckl_get_jastrow_champ_een_rescaled_n_gl",
                            "Array too small. Expected ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1)");
   }
-  memcpy(distance_rescaled, ctx->jastrow_champ.een_rescaled_n_gl, sze * sizeof(double));
+
+  memcpy(een_rescaled_n_gl, ctx->jastrow_champ.een_rescaled_n_gl, sze * sizeof(double));
 
   return QMCKL_SUCCESS;
 }
@@ -8159,12 +8441,22 @@ qmckl_exit_code qmckl_compute_dim_c_vector (
 **** Get
 
      #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes
-qmckl_exit_code qmckl_get_jastrow_champ_tmp_c(qmckl_context context, double* const tmp_c);
-qmckl_exit_code qmckl_get_jastrow_champ_dtmp_c(qmckl_context context, double* const dtmp_c);
+qmckl_exit_code
+qmckl_get_jastrow_champ_tmp_c(qmckl_context context,
+                              double* const tmp_c,
+                              const int64_t size_max);
+
+qmckl_exit_code
+qmckl_get_jastrow_champ_dtmp_c(qmckl_context context,
+                               double* const dtmp_c,
+                               const int64_t size_max);
      #+end_src
 
      #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_get_jastrow_champ_tmp_c(qmckl_context context, double* const tmp_c)
+qmckl_exit_code
+qmckl_get_jastrow_champ_tmp_c(qmckl_context context,
+                              double* const tmp_c,
+                              const int64_t size_max)
 {
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return QMCKL_NULL_CONTEXT;
@@ -8181,14 +8473,33 @@ qmckl_exit_code qmckl_get_jastrow_champ_tmp_c(qmckl_context context, double* con
   qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
-  size_t sze = (ctx->jastrow_champ.cord_num) * (ctx->jastrow_champ.cord_num + 1)
-               ,* ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num;
+  if (tmp_c == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_tmp_c",
+                           "Null pointer");
+  }
+
+  const int64_t sze = (ctx->jastrow_champ.cord_num) * (ctx->jastrow_champ.cord_num + 1) *
+    ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num;
+
+  if (size_max < sze) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_get_jastrow_champ_tmp_c",
+                           "Array too small. Expected cord_num*(cord_num+1)*walk_num*elec_num*nucl_num");
+  }
+
   memcpy(tmp_c, ctx->jastrow_champ.tmp_c, sze * sizeof(double));
 
   return QMCKL_SUCCESS;
 }
 
-qmckl_exit_code qmckl_get_jastrow_champ_dtmp_c(qmckl_context context, double* const dtmp_c)
+
+qmckl_exit_code
+qmckl_get_jastrow_champ_dtmp_c(qmckl_context context,
+                               double* const dtmp_c,
+                               const int64_t size_max)
 {
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return QMCKL_NULL_CONTEXT;
@@ -8205,8 +8516,23 @@ qmckl_exit_code qmckl_get_jastrow_champ_dtmp_c(qmckl_context context, double* co
   qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
-  size_t sze = (ctx->jastrow_champ.cord_num) * (ctx->jastrow_champ.cord_num + 1)
-    ,*4* ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num;
+  if (dtmp_c == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_dtmp_c",
+                           "Null pointer");
+  }
+
+  const int64_t sze = (ctx->jastrow_champ.cord_num) * (ctx->jastrow_champ.cord_num + 1)*
+    4* ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num;
+
+  if (size_max < sze) {
+    return qmckl_failwith(context,
+                          QMCKL_INVALID_ARG_3,
+                          "qmckl_get_jastrow_champ_dtmp_c",
+                          "Array too small. Expected 4*cord_num*(cord_num+1)*walk_num*elec_num*nucl_num");
+  }
+
   memcpy(dtmp_c, ctx->jastrow_champ.dtmp_c, sze * sizeof(double));
 
   return QMCKL_SUCCESS;
@@ -9353,10 +9679,10 @@ lkpm_of_cindex = np.array(lkpm_combined_index).T
 assert(qmckl_electron_provided(context));
 
 double tmp_c[walk_num][cord_num][cord_num+1][nucl_num][elec_num];
-rc = qmckl_get_jastrow_champ_tmp_c(context, &(tmp_c[0][0][0][0][0]));
+rc = qmckl_get_jastrow_champ_tmp_c(context, &(tmp_c[0][0][0][0][0]), sizeof(tmp_c)/sizeof(double));
 
 double dtmp_c[walk_num][cord_num][cord_num+1][nucl_num][4][elec_num];
-rc = qmckl_get_jastrow_champ_dtmp_c(context, &(dtmp_c[0][0][0][0][0][0]));
+rc = qmckl_get_jastrow_champ_dtmp_c(context, &(dtmp_c[0][0][0][0][0][0]), sizeof(dtmp_c)/sizeof(double));
 
 printf("%e\n%e\n", tmp_c[0][0][1][0][0],  3.954384);
 assert(fabs(tmp_c[0][0][1][0][0] - 3.954384) < 1e-6);
@@ -9398,13 +9724,22 @@ qmckl_get_jastrow_champ_factor_een(qmckl_context context,
   qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
-  int64_t sze = ctx->electron.walker.num;
+  if (factor_een == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_factor_een",
+                           "Null pointer");
+  }
+
+  const int64_t sze = ctx->electron.walker.num;
+
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
                            "qmckl_get_jastrow_champ_factor_een",
                            "Array too small. Expected walk_num");
   }
+
   memcpy(factor_een, ctx->jastrow_champ.factor_een, sze*sizeof(double));
 
   return QMCKL_SUCCESS;
@@ -9970,13 +10305,22 @@ qmckl_get_jastrow_champ_factor_een_gl(qmckl_context context,
   qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
-  int64_t sze = ctx->electron.walker.num * 4 * ctx->electron.num;
+  if (factor_een_gl == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_factor_een_gl",
+                           "Null pointer");
+  }
+
+  const int64_t sze = ctx->electron.walker.num * 4 * ctx->electron.num;
+
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
                            "qmckl_get_jastrow_champ_factor_een_gl",
                            "Array too small. Expected 4*walk_num*elec_num");
   }
+
   memcpy(factor_een_gl, ctx->jastrow_champ.factor_een_gl, sze*sizeof(double));
 
   return QMCKL_SUCCESS;
@@ -10887,13 +11231,22 @@ qmckl_get_jastrow_champ_value(qmckl_context context,
   rc = qmckl_provide_jastrow_champ_value(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  int64_t sze=ctx->electron.walker.num;
+  if (value == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_value",
+                           "Null pointer");
+  }
+
+  const int64_t sze = ctx->electron.walker.num;
+
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
                            "qmckl_get_jastrow_champ_value",
                            "Array too small. Expected walk_num");
   }
+
   memcpy(value, ctx->jastrow_champ.value, sze*sizeof(double));
 
   return QMCKL_SUCCESS;
@@ -11244,13 +11597,22 @@ qmckl_get_jastrow_champ_gl(qmckl_context context,
   rc = qmckl_provide_jastrow_champ_gl(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  int64_t sze = 4 * ctx->electron.walker.num * ctx->electron.num;
+  if (gl == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_2,
+                           "qmckl_get_jastrow_champ_gl",
+                           "Null pointer");
+  }
+
+  const int64_t sze = 4 * ctx->electron.walker.num * ctx->electron.num;
+  
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
                            "qmckl_get_jastrow_champ_gl",
                            "Array too small. Expected walk_num*elec_num*4");
   }
+
   memcpy(gl, ctx->jastrow_champ.gl, sze*sizeof(double));
 
   return QMCKL_SUCCESS;
@@ -11575,44 +11937,49 @@ printf("Total Jastrow derivatives\n");
 /* Check if Jastrow is properly initialized */
 assert(qmckl_jastrow_champ_provided(context));
 
-rc = qmckl_check(context,
-                 qmckl_get_jastrow_champ_factor_ee_gl(context, &(factor_ee_gl[0][0][0]), walk_num*elec_num*4)
-                 );
-assert(rc == QMCKL_SUCCESS);
+{
+  double factor_ee_gl[walk_num][4][elec_num];
+  rc = qmckl_check(context,
+                   qmckl_get_jastrow_champ_factor_ee_gl(context, &(factor_ee_gl[0][0][0]), walk_num*elec_num*4)
+                   );
+  assert(rc == QMCKL_SUCCESS);
 
-double factor_en_gl[walk_num][4][elec_num];
-rc = qmckl_check(context,
-                 qmckl_get_jastrow_champ_factor_en_gl(context, &(factor_en_gl[0][0][0]), walk_num*elec_num*4)
-                 );
-assert(rc == QMCKL_SUCCESS);
+  double factor_en_gl[walk_num][4][elec_num];
+  rc = qmckl_check(context,
+                   qmckl_get_jastrow_champ_factor_en_gl(context, &(factor_en_gl[0][0][0]), walk_num*elec_num*4)
+                   );
+  assert(rc == QMCKL_SUCCESS);
 
-rc = qmckl_check(context,
-                 qmckl_get_jastrow_champ_factor_een_gl(context, &(factor_een_gl[0][0][0]), walk_num*elec_num*4)
-                 );
-assert(rc == QMCKL_SUCCESS);
+  double factor_een_gl[walk_num][4][elec_num];
+  rc = qmckl_check(context,
+                   qmckl_get_jastrow_champ_factor_een_gl(context, &(factor_een_gl[0][0][0]), walk_num*elec_num*4)
+                   );
+  assert(rc == QMCKL_SUCCESS);
 
-double total_j_deriv[walk_num][4][elec_num];
-rc = qmckl_check(context,
-                 qmckl_get_jastrow_champ_gl(context, &(total_j_deriv[0][0][0]), walk_num*elec_num*4)
-                 );
-assert(rc == QMCKL_SUCCESS);
+  double total_j_deriv[walk_num][4][elec_num];
+  rc = qmckl_check(context,
+                   qmckl_get_jastrow_champ_gl(context, &(total_j_deriv[0][0][0]), walk_num*elec_num*4)
+                   );
+  assert(rc == QMCKL_SUCCESS);
 
-rc = qmckl_check(context,
-                 qmckl_get_jastrow_champ_value(context, &(total_j[0]), walk_num)
-                 );
-assert(rc == QMCKL_SUCCESS);
+  double total_j[walk_num];
+  rc = qmckl_check(context,
+                   qmckl_get_jastrow_champ_value(context, &(total_j[0]), walk_num)
+                   );
+  assert(rc == QMCKL_SUCCESS);
 
 
-for (int64_t k=0 ; k< walk_num ; ++k) {
-  for (int64_t m=0 ; m<4; ++m) {
-    for (int64_t e=0 ; e<elec_num; ++e) {
-      if (m < 3) { /* test only gradients */
-        assert (total_j_deriv[k][m][e]/total_j[k] - (factor_ee_gl[k][m][e] + factor_en_gl[k][m][e] + factor_een_gl[k][m][e]) < 1.e-12);
+  for (int64_t k=0 ; k< walk_num ; ++k) {
+    for (int64_t m=0 ; m<4; ++m) {
+      for (int64_t e=0 ; e<elec_num; ++e) {
+        if (m < 3) { /* test only gradients */
+          assert (total_j_deriv[k][m][e]/total_j[k] - (factor_ee_gl[k][m][e] + factor_en_gl[k][m][e] + factor_een_gl[k][m][e]) < 1.e-12);
+        }
       }
     }
   }
- }
 
+}
 
       #+end_src