diff --git a/org/qmckl_jastrow_champ.org b/org/qmckl_jastrow_champ.org
index 3f6314b..f51f6be 100644
--- a/org/qmckl_jastrow_champ.org
+++ b/org/qmckl_jastrow_champ.org
@@ -68,10 +68,11 @@
 ** Reformulation of the three-body part
 
    To accelerate the computation of the three-body part, the Jastrow
-   factor is re-expressed as follows:
+   factor is re-expressed as follows, with $m=(p-k)/2 -l/2$:
+
 
 \begin{eqnarray*}
-   & & \sum_{\alpha=1}^{N_{\text{nucl}}}
+   J_{kpl} & = & \sum_{\alpha=1}^{N_{\text{nucl}}}
      \sum_{i=1}^{N_{\text{elec}}}
       \sum_{j=1}^{i-1}
            c_{lkp\alpha} \left[ g_\text{e}({r}_{ij}) \right]^k
@@ -93,8 +94,35 @@
              \left[ \left[ g_\alpha({R}_{i\alpha}) \right]^{l+m} \left[ g_\alpha({R}_{j\alpha}) \right]^{m} +
                     \left[ g_\alpha({R}_{i\alpha}) \right]^{l} \left[
     g_\alpha({R}_{j\alpha}) \right]^{l+m} \right] \\
+   & = & 
+    \sum_{\alpha=1}^{N_{\text{nucl}}} c_{lkp\alpha}
+     \sum_{i=1}^{N_{\text{elec}}}
+      \sum_{j=1}^{N_{\text{elec}}} 
+             \left[ g_\alpha({R}_{i\alpha}) \right]^{l+m}
+             \left[ g_\text{e}({r}_{ij}) \right]^k
+             \left[ g_\alpha({R}_{j\alpha}) \right]^{m} \\
+   & = & 
+    \sum_{\alpha=1}^{N_{\text{nucl}}} c_{lkp\alpha}
+     \sum_{i=1}^{N_{\text{elec}}} 
+             \left[ g_\alpha({R}_{i\alpha}) \right]^{l+m}
+             P_{i\alpha}^{km}, \text{ with }
+             P_{i\alpha}^{km} =
+      \sum_{j=1}^{N_{\text{elec}}} 
+             \left[ g_\text{e}({r}_{ij}) \right]^k
+             \left[ g_\alpha({R}_{j\alpha}) \right]^{m}. \\
+J   & = & 
+   \sum_{p=2}^{N_{\text{ord}}}
+    \sum_{k=0}^{p-1}
+     \sum_{l=0}^{p-k-2\delta_{k,0}}
+      \sum_{\alpha=1}^{N_{\text{nucl}}} c_{lkp\alpha}
+       \sum_{i=1}^{N_{\text{elec}}}
+        \left[ g_\alpha({R}_{i\alpha}) \right]^{(p-k+l)/2}
+             P_{i\alpha}^{k, (p-k-l)/2}
    \end{eqnarray*}
 
+   The computation of $P$ scales as   $\mathcal{O}(N_\text{elec}^2 N_\text{nucl}n^2)$, and 
+   the computation of $J$ scales as $\mathcal{O}(N_\text{elec}N_\text{nucl}n^2)$.
+
 * Headers                                                          :noexport:
   #+begin_src elisp :noexport :results none
 (org-babel-lob-ingest "../tools/lib.org")
@@ -126,6 +154,11 @@
 #include "qmckl_context_private_type.h"
 #include "qmckl_jastrow_champ_private_func.h"
 
+#ifdef __GNUC__
+#pragma GCC push_options
+#pragma GCC optimize ("O0")
+#endif
+
 int main() {
   qmckl_context context;
   context = qmckl_context_create();
@@ -1798,7 +1831,6 @@ for (int64_t i=0 ; i<nucl_num ; ++i) {
   assert( nucl_charge[i] == nucl_charge2[i] );
  }
 assert(qmckl_nucleus_provided(context));
-
     #+end_src
 
 * Computation
@@ -2081,6 +2113,7 @@ qmckl_exit_code qmckl_compute_jastrow_champ_asymp_jasb (const qmckl_context cont
     (context, bord_num, b_vector, rescale_factor_ee, spin_independent, asymp_jasb);
 }
  #+end_src
+
 **** Test                                                          :noexport:
      #+name: asymp_jasb
      #+begin_src python :results output :exports none :noweb yes
@@ -3486,6 +3519,7 @@ qmckl_compute_jastrow_champ_factor_ee (const qmckl_context context,
      ee_distance_rescaled, asymp_jasb, spin_independent, factor_ee);
 }
  #+end_src
+
 **** Test                                                          :noexport:
      #+begin_src python :results output :exports none :noweb yes
 import numpy as np
@@ -3515,37 +3549,14 @@ for i in range(0,elec_num):
                           / (1.0 + b_vector[1] * ee_distance_rescaled[i][j])     \
                          - asymp_jasb[ipar] + pow_ser
 print("factor_ee        :",factor_ee)
-print("ee_distance_rescaled :",ee_distance_rescaled)
 
      #+end_src
 
      #+RESULTS:
-     #+begin_example
-     asym_one         :  0.6634291325000664
-     asymp_jasb[0]    :  0.7115733522582638
-     asymp_jasb[1]    :  1.043287918508297
-     factor_ee        : -16.83886184243964
-     ee_distance_rescaled : [[0.         0.63475074 1.29816415 1.23147027 1.51933127 0.54402406
-       0.51452479 0.96538731 1.25878564 1.3722995 ]
-      [0.63475074 0.         1.35148664 1.13524156 1.48940503 0.4582292
-       0.62758076 1.06560856 1.179133   1.30763703]
-      [1.29816415 1.35148664 0.         1.50021375 1.59200788 1.23488312
-       1.20844259 1.0355537  1.52064535 1.53049239]
-      [1.23147027 1.13524156 1.50021375 0.         1.12016142 1.19158954
-       1.29762585 1.24824277 0.25292267 0.58609336]
-      [1.51933127 1.48940503 1.59200788 1.12016142 0.         1.50217017
-       1.54012828 1.48753895 1.10441805 0.84504381]
-      [0.54402406 0.4582292  1.23488312 1.19158954 1.50217017 0.
-       0.39417354 0.87009603 1.23838502 1.33419121]
-      [0.51452479 0.62758076 1.20844259 1.29762585 1.54012828 0.39417354
-       0.         0.95118809 1.33068934 1.41097406]
-      [0.96538731 1.06560856 1.0355537  1.24824277 1.48753895 0.87009603
-       0.95118809 0.         1.29422213 1.33222493]
-      [1.25878564 1.179133   1.52064535 0.25292267 1.10441805 1.23838502
-       1.33068934 1.29422213 0.         0.62196802]
-      [1.3722995  1.30763703 1.53049239 0.58609336 0.84504381 1.33419121
-       1.41097406 1.33222493 0.62196802 0.        ]]
-     #+end_example
+     : asym_one         :  0.6634291325000664
+     : asymp_jasb[0]    :  0.7115733522582638
+     : asymp_jasb[1]    :  1.043287918508297
+     : factor_ee        : -16.83886184243964
 
       #+begin_src c :tangle (eval c_test)
 /* Check if Jastrow is properly initialized */
@@ -4218,10 +4229,10 @@ print("factor_ee_gl[3][0]:",factor_ee_gl[3][0])
      : asym_one         :  0.6634291325000664
      : asymp_jasb[0]    :  0.7115733522582638
      : asymp_jasb[1]    :  1.043287918508297
-     : factor_ee_gl[0][0]:
-     : factor_ee_gl[1][0]:
-     : factor_ee_gl[2][0]:
-     : factor_ee_gl[3][0]:
+     : factor_ee_gl[0][0]: -0.39319353942687446
+     : factor_ee_gl[1][0]: 1.0535615450668214
+     : factor_ee_gl[2][0]: -0.39098406960784515
+     : factor_ee_gl[3][0]: 2.8650469630854483
 
       #+begin_src c :tangle (eval c_test)
 /* Check if Jastrow is properly initialized */
@@ -5555,7 +5566,7 @@ assert(qmckl_electron_provided(context));
    coeffecients and the electron-nucleus rescaled distances ~en_distance_rescaled~.
 
    \[
- f_{\alpha}(R_{i\alpha}) = - \sum_{i,j<i} \left[ \frac{ A_0 C_{ij}}{1 - A_1 C_{ij}} + \sum^{N^\alpha_{\text{ord}}}_{k}A_k C_{ij}^k \right]
+ f_{\alpha}(R_{i\alpha}) = - \sum_{i,j<i} \left[ \frac{ A_0 C_{ij}}{1 + A_1 C_{ij}} + \sum_{k=2}^{N^\alpha_{\text{ord}}}A_k C_{ij}^k \right]
    \]
 
 
@@ -6023,7 +6034,21 @@ assert(fabs(22.781375792083587 - factor_en[0]) < 1.e-12);
    with respect to the electron coordinates using the ~en_distance_rescaled~ and
    ~en_distance_rescaled_gl~ which are already calculated previously.
 
-   TODO: write equations.
+   The derivative is calculated in the function ~qmckl_compute_jastrow_champ_factor_en_gl~.
+   The formula is given by:
+   \[
+\nabla_i f_{\alpha}(R_{i\alpha}) = \sum_{j=1}^{N_\text{elec}} \left[
+   \frac{ A_0\, \nabla_i C_{ij} }{(1 + A_1 C_{ij})^2} + 
+  \sum_{k=2}^{N^\alpha_{\text{ord}}} A_k\, k\, C_{ij}^{k-1}\,\nabla_i C_{ij} \right] 
+   \]
+   
+   \[
+\Delta_i f_{\alpha}(R_{i\alpha}) = \sum_{j=1}^{N_\text{elec}} \left[
+   \frac{ A_0\, \Delta_i C_{ij} }{(1 + A_1 C_{ij})^2} - 
+   \frac{ 2 A_0\, A_1 (\nabla_i C_{ij})^2}{(1 + A_1 C_{ij})^3} + 
+  \sum_{k=2}^{N^\alpha_{\text{ord}}} A_k\, k\, C_{ij}^{k-1} C_{ij}^{k-2}
+\left[ \Delta_i C_{ij} + (k-1)(\nabla_i C_{ij})^2 \right] \right] 
+   \]
 
 **** Get
      #+begin_src c :comments org :tangle (eval h_func) :noweb yes
@@ -6186,6 +6211,7 @@ qmckl_exit_code qmckl_provide_jastrow_champ_factor_en_gl(qmckl_context context)
      :END:
 
      #+NAME: qmckl_factor_en_gl_args
+     |---------------------------+-------------------------------------------+--------+---------------------------------------|
      | Variable                  | Type                                      | In/Out | Description                           |
      |---------------------------+-------------------------------------------+--------+---------------------------------------|
      | ~context~                 | ~qmckl_context~                           | in     | Global state                          |
@@ -6199,6 +6225,7 @@ qmckl_exit_code qmckl_provide_jastrow_champ_factor_en_gl(qmckl_context context)
      | ~en_distance_rescaled~    | ~double[walk_num][nucl_num][elec_num]~    | in     | Electron-nucleus distances            |
      | ~en_distance_rescaled_gl~ | ~double[walk_num][nucl_num][elec_num][4]~ | in     | Electron-nucleus distance derivatives |
      | ~factor_en_gl~            | ~double[walk_num][4][elec_num]~           | out    | Electron-nucleus jastrow              |
+     |---------------------------+-------------------------------------------+--------+---------------------------------------|
 
      #+begin_src f90 :comments org :tangle (eval f) :noweb yes
 function qmckl_compute_jastrow_champ_factor_en_gl_doc( &
@@ -6256,30 +6283,30 @@ function qmckl_compute_jastrow_champ_factor_en_gl_doc( &
   endif
 
   do nw =1, walk_num
-    factor_en_gl(:,:,nw) = 0.0d0
-    do a = 1, nucl_num
-       do i = 1, elec_num
+     factor_en_gl(:,:,nw) = 0.0d0
+     do a = 1, nucl_num
+        do i = 1, elec_num
 
-          x = en_distance_rescaled(i,a,nw)
-          if(abs(x) < 1.d-12) continue
+           x = en_distance_rescaled(i,a,nw)
+           if(abs(x) < 1.d-12) continue
 
-          denom = 1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x
-          invdenom = 1.0d0 / denom
-          invdenom2 = invdenom*invdenom
+           denom = 1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x
+           invdenom = 1.0d0 / denom
+           invdenom2 = invdenom*invdenom
 
-          dx(1) = en_distance_rescaled_gl(1,i,a,nw)
-          dx(2) = en_distance_rescaled_gl(2,i,a,nw)
-          dx(3) = en_distance_rescaled_gl(3,i,a,nw)
-          dx(4) = en_distance_rescaled_gl(4,i,a,nw)
+           dx(1) = en_distance_rescaled_gl(1,i,a,nw)
+           dx(2) = en_distance_rescaled_gl(2,i,a,nw)
+           dx(3) = en_distance_rescaled_gl(3,i,a,nw)
+           dx(4) = en_distance_rescaled_gl(4,i,a,nw)
 
-          f = a_vector(1, type_nucl_vector(a)+1) * invdenom2
-          grad_c2 = dx(1)*dx(1) + dx(2)*dx(2) + dx(3)*dx(3)
+           f = a_vector(1, type_nucl_vector(a)+1) * invdenom2
+           grad_c2 = dx(1)*dx(1) + dx(2)*dx(2) + dx(3)*dx(3)
 
-          factor_en_gl(i,1,nw) = factor_en_gl(i,1,nw) + f * dx(1)
-          factor_en_gl(i,2,nw) = factor_en_gl(i,2,nw) + f * dx(2)
-          factor_en_gl(i,3,nw) = factor_en_gl(i,3,nw) + f * dx(3)
-          factor_en_gl(i,4,nw) = factor_en_gl(i,4,nw) &
-               + f * (dx(4) - 2.d0 * a_vector(2, type_nucl_vector(a)+1) * grad_c2 * invdenom)
+           factor_en_gl(i,1,nw) = factor_en_gl(i,1,nw) + f * dx(1)
+           factor_en_gl(i,2,nw) = factor_en_gl(i,2,nw) + f * dx(2)
+           factor_en_gl(i,3,nw) = factor_en_gl(i,3,nw) + f * dx(3)
+           factor_en_gl(i,4,nw) = factor_en_gl(i,4,nw) &
+                + f * (dx(4) - 2.d0 * a_vector(2, type_nucl_vector(a)+1) * grad_c2 * invdenom)
 
 
            kf = 2.d0
@@ -6714,13 +6741,34 @@ assert(qmckl_jastrow_champ_provided(context));
       #+end_src
 
 ** Electron-electron-nucleus component
+
+     #+name: en_dist
+     #+begin_src python :results output :exports none :noweb yes
+elec_coord = np.array(elec_coord)[0]
+nucl_coord = np.array(nucl_coord)
+elnuc_dist = np.zeros(shape=(nucl_num, elec_num),dtype=float)
+for i in range(elec_num):
+  for a in range(nucl_num):
+    elnuc_dist[a,i] = np.linalg.norm(elec_coord[i] - nucl_coord[:,a])
+
+kappa_n = 0.6
+
+een_rescaled_n = np.zeros(shape=(cord_num+1, nucl_num, elec_num), dtype=float)
+
+for p in range(0, cord_num+1):
+  for a in range(nucl_num):
+    for i in range(elec_num):
+      een_rescaled_n[p,a,i] = np.exp(-kappa_n * p * elnuc_dist[a,i])
+
+     #+end_src
+
 *** Electron-electron rescaled distances in $J_\text{eeN}$
 
     ~een_rescaled_e~ stores the table of the rescaled distances between all
     pairs of electrons and raised to the power \(p\) defined by ~cord_num~:
 
     \[
-    C_{ij,p} = \left[ \exp\left(-\kappa_\text{e}\, r_{ij}\right) \right]^p
+    [g_e(r)_{ij}]^p = \exp\left(-p\,\kappa_\text{e}\, r_{ij}\right)
     \]
 
    where \(r_{ij}\) is the matrix of electron-electron distances.
@@ -6880,12 +6928,11 @@ function qmckl_compute_een_rescaled_e_doc( &
   real    (c_double ) , intent(out)         :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
   integer(qmckl_exit_code)                  :: info
 
-  double precision,dimension(:,:),allocatable :: een_rescaled_e_ij
   double precision                    :: x
   integer*8                           :: i, j, k, l, nw
 
   info = QMCKL_SUCCESS
-
+  
   if (context == QMCKL_NULL_CONTEXT) then
      info = QMCKL_INVALID_CONTEXT
      return
@@ -6906,47 +6953,14 @@ function qmckl_compute_een_rescaled_e_doc( &
      return
   endif
 
-  allocate(een_rescaled_e_ij(elec_num * (elec_num - 1) / 2, cord_num + 1))
-
-  ! Prepare table of exponentiated distances raised to appropriate power
   do nw = 1, walk_num
-     een_rescaled_e_ij(:, 1) = 1.0d0
-
-
-     k = 0
-     do j = 1, elec_num
-        do i = 1, j - 1
-           k = k + 1
-           een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_ee * ee_distance(i, j, nw))
-        end do
-     end do
-
-
-     do l = 2, cord_num
-        do k = 1, elec_num * (elec_num - 1)/2
-           een_rescaled_e_ij(k, l+1) = een_rescaled_e_ij(k, l) * een_rescaled_e_ij(k, 2)
-        end do
-     end do
-
-     ! prepare the actual een table
-     een_rescaled_e(:, :, 0, nw) = 1.0d0
-     do j = 1, elec_num
-        een_rescaled_e(j, j, 0, nw) = 0.0d0
-     end do
-
-     do l = 1, cord_num
-        k = 0
+     do l = 0, cord_num
         do j = 1, elec_num
-           do i = 1, j - 1
-              k = k + 1
-              x = een_rescaled_e_ij(k, l+1)
-              een_rescaled_e(i, j, l, nw) = x
-              een_rescaled_e(j, i, l, nw) = x
+           do i = 1, elec_num
+              een_rescaled_e(i, j, l, nw) = dexp(-rescale_factor_ee * ee_distance(i, j, nw))**l
            end do
-           een_rescaled_e(j, j, l, nw) = 0.0d0
         end do
      end do
-
   end do
 
 end function qmckl_compute_een_rescaled_e_doc
@@ -6974,6 +6988,13 @@ qmckl_exit_code qmckl_compute_een_rescaled_e_hpc (const qmckl_context context,
                                                   const double* ee_distance,
                                                   double* const een_rescaled_e ) {
 
+return qmckl_compute_een_rescaled_e_doc (context,
+                                                  walk_num,
+                                                  elec_num,
+                                                  cord_num,
+                                                  rescale_factor_ee,
+                                                  ee_distance,
+                                                  een_rescaled_e ) ;
   if (context == QMCKL_NULL_CONTEXT) {
     return QMCKL_INVALID_CONTEXT;
   }
@@ -7153,78 +7174,46 @@ qmckl_exit_code qmckl_compute_een_rescaled_e_hpc (const qmckl_context context,
 
 **** Test
 
+     #+name: ee_dist
      #+begin_src python :results output :exports none :noweb yes
 import numpy as np
-
 <<jastrow_data>>
 
 elec_coord = np.array(elec_coord)[0]
 elec_dist = np.zeros(shape=(elec_num, elec_num),dtype=float)
+
 for i in range(elec_num):
   for j in range(elec_num):
-    elec_dist[i, j] = np.linalg.norm(elec_coord[i] - elec_coord[j])
+    elec_dist[i,j] = np.linalg.norm(elec_coord[i] - elec_coord[j])
 
-kappa = 0.6
+kappa_e = 0.6
 
-een_rescaled_e_ij = np.zeros(shape=(elec_num * (elec_num - 1)//2, cord_num+1), dtype=float)
-een_rescaled_e_ij[:,0] = 1.0
-
-k = 0
-for j in range(elec_num):
-  for i in range(j):
-    een_rescaled_e_ij[k, 1] = np.exp(-kappa * elec_dist[i, j])
-    k = k + 1
-
-for l in range(2, cord_num + 1):
-  for k in range(elec_num * (elec_num - 1)//2):
-    een_rescaled_e_ij[k, l] = een_rescaled_e_ij[k, l - 1] * een_rescaled_e_ij[k, 1]
-
-een_rescaled_e = np.zeros(shape=(elec_num, elec_num, cord_num + 1), dtype=float)
-een_rescaled_e[:,:,0] = 1.0
-
-for l in range(1,cord_num+1):
-  k = 0
+een_rescaled_e = np.zeros(shape=(cord_num+1, elec_num, elec_num), dtype=float)
+for i in range(elec_num):
   for j in range(elec_num):
-    for i in range(j):
-      x = een_rescaled_e_ij[k, l]
-      een_rescaled_e[i, j, l] = x
-      een_rescaled_e[j, i, l] = x
-      k = k + 1
-
-for l in range(0,cord_num+1):
-  for j in range(0, elec_num):
-    een_rescaled_e[j,j,l] = 0.0
-
-print(" een_rescaled_e[0, 2, 1] = ",een_rescaled_e[0, 2, 1])
-print(" een_rescaled_e[0, 3, 1] = ",een_rescaled_e[0, 3, 1])
-print(" een_rescaled_e[0, 4, 1] = ",een_rescaled_e[0, 4, 1])
-print(" een_rescaled_e[1, 3, 2] = ",een_rescaled_e[1, 3, 2])
-print(" een_rescaled_e[1, 4, 2] = ",een_rescaled_e[1, 4, 2])
-print(" een_rescaled_e[1, 5, 2] = ",een_rescaled_e[1, 5, 2])
+    for p in range(0, cord_num+1):
+      een_rescaled_e[p,i,j] = np.exp(-kappa_e * p * elec_dist[i,j])
      #+end_src
 
-     #+RESULTS:
-     :  een_rescaled_e[0, 2, 1] =  0.2211015082992776
-     :  een_rescaled_e[0, 3, 1] =  0.2611178387068169
-     :  een_rescaled_e[0, 4, 1] =  0.08840123507637472
-     :  een_rescaled_e[1, 3, 2] =  0.10166855073546568
-     :  een_rescaled_e[1, 4, 2] =  0.011311807324686948
-     :  een_rescaled_e[1, 5, 2] =  0.5257156022077619
+     #+NAME:test_ee
+     #+begin_src python :results output :exports none :noweb yes
+<<ee_dist>>
+for p in range(cord_num+1):
+#  for i in range(elec_num):
+#    for j in range(elec_num):
+  for i in range(3):
+    for j in range(3):
+      print(f"assert( fabs(een_rescaled_e[0][{p}][{i}][{j}] - ({een_rescaled_e[p,i,j]})) < 1.e-10 );")
+     #+end_src
 
-      #+begin_src c :tangle (eval c_test)
+      #+begin_src c :tangle (eval c_test) :noweb yes
 assert(qmckl_electron_provided(context));
 {
 
-double een_rescaled_e[walk_num][(cord_num + 1)][elec_num][elec_num];
-rc = qmckl_get_jastrow_champ_een_distance_rescaled_e(context, &(een_rescaled_e[0][0][0][0]),elec_num*elec_num*(cord_num+1)*walk_num);
+  double een_rescaled_e[walk_num][(cord_num + 1)][elec_num][elec_num];
+  rc = qmckl_get_jastrow_champ_een_distance_rescaled_e(context, &(een_rescaled_e[0][0][0][0]),elec_num*elec_num*(cord_num+1)*walk_num);
 
-// value of (0,2,1)
-assert(fabs(een_rescaled_e[0][1][0][2]- 0.2211015082992776 ) < 1.e-12);
-assert(fabs(een_rescaled_e[0][1][0][3]- 0.2611178387068169 ) < 1.e-12);
-assert(fabs(een_rescaled_e[0][1][0][4]- 0.0884012350763747 ) < 1.e-12);
-assert(fabs(een_rescaled_e[0][2][1][3]- 0.1016685507354656 ) < 1.e-12);
-assert(fabs(een_rescaled_e[0][2][1][4]- 0.0113118073246869 ) < 1.e-12);
-assert(fabs(een_rescaled_e[0][2][1][5]- 0.5257156022077619 ) < 1.e-12);
+  <<test_ee()>>
 }
 
 {
@@ -7236,13 +7225,13 @@ assert(fabs(een_rescaled_e[0][2][1][5]- 0.5257156022077619 ) < 1.e-12);
   double een_rescaled_e_doc[walk_num][cord_num+1][elec_num][elec_num];
   memset(&(een_rescaled_e_doc[0][0][0][0]), 0, sizeof(een_rescaled_e_doc));
   rc = qmckl_compute_een_rescaled_e(context, walk_num, elec_num, cord_num,
-                0.6, &(ee_distance[0]), &(een_rescaled_e_doc[0][0][0][0]));
+                                    rescale_factor_ee, &(ee_distance[0]), &(een_rescaled_e_doc[0][0][0][0]));
   assert(rc == QMCKL_SUCCESS);
 
   double een_rescaled_e_hpc[walk_num][cord_num+1][elec_num][elec_num];
   memset(&(een_rescaled_e_hpc[0][0][0][0]), 0, sizeof(een_rescaled_e_hpc));
   rc = qmckl_compute_een_rescaled_e_hpc(context, walk_num, elec_num, cord_num,
-                0.6, &(ee_distance[0]), &(een_rescaled_e_hpc[0][0][0][0]));
+                                        rescale_factor_ee, &(ee_distance[0]), &(een_rescaled_e_hpc[0][0][0][0]));
   assert(rc == QMCKL_SUCCESS);
 
   for (int64_t i = 0; i < walk_num; i++) {
@@ -7270,7 +7259,12 @@ assert(fabs(een_rescaled_e[0][2][1][5]- 0.5257156022077619 ) < 1.e-12);
 
     \[ \frac{\partial}{\partial x} \left[ {g_\text{e}(r)}\right]^p =
     -\frac{x}{r} \kappa_\text{e}\, p\,\left[ {g_\text{e}(r)}\right]^p \]
-    \[ \Delta \left[ {g_\text{e}(r)}\right]^p = \frac{2}{r} \kappa_\text{e}\, p\,\left[ {g_\text{e}(r)}\right]^p \right]  + \left(\frac{\partial}{\partial x}\left[ {g_\text{e}(r)}\right]^p \right)^2 + \left(\frac{\partial}{\partial y}\left[ {g_\text{e}(r)}\right]^p \right)^2 + \left(\frac{\partial}{\partial z}\left[ {g_\text{e}(r)}\right]^p \right)^2 \]
+
+    \[ \Delta \left[ {g_\text{e}(r)}\right]^p = \kappa_\text{e}\,
+    p\,\left[ - \frac{2}{r} + \kappa_\text{e}\, p \right] \left[
+    {g_\text{e}(r)} \right]^p \]
+
+    Derivatives are set to zero at $r_{ii}$ to avoid NaNs.
 
 **** Get
 
@@ -7408,6 +7402,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_gl(qmckl_context context)
      :END:
 
      #+NAME: qmckl_factor_een_rescaled_e_gl_args
+     |---------------------+-------------------------------------------------------+--------+--------------------------------------|
      | Variable            | Type                                                  | In/Out | Description                          |
      |---------------------+-------------------------------------------------------+--------+--------------------------------------|
      | ~context~           | ~qmckl_context~                                       | in     | Global state                         |
@@ -7419,6 +7414,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_gl(qmckl_context context)
      | ~ee_distance~       | ~double[walk_num][elec_num][elec_num]~                | in     | Electron-electron distances          |
      | ~een_rescaled_e~    | ~double[walk_num][0:cord_num][elec_num][elec_num]~    | in     | Electron-electron distances          |
      | ~een_rescaled_e_gl~ | ~double[walk_num][0:cord_num][elec_num][4][elec_num]~ | out    | Electron-electron rescaled distances |
+     |---------------------+-------------------------------------------------------+--------+--------------------------------------|
 
      #+begin_src f90 :comments org :tangle (eval f) :noweb yes
 function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_doc( &
@@ -7433,7 +7429,7 @@ function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_doc( &
   integer(c_int64_t)    , intent(in), value  :: elec_num
   integer(c_int64_t)    , intent(in), value  :: cord_num
   real(c_double)        , intent(in), value  :: rescale_factor_ee
-  real(c_double)        , intent(in)  :: coord_ee(elec_num,3,walk_num)
+  real(c_double)        , intent(in)  :: coord_ee(elec_num,walk_num,3)
   real(c_double)        , intent(in)  :: ee_distance(elec_num,elec_num,walk_num)
   real(c_double)        , intent(in)  :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
   real(c_double)        , intent(out) :: een_rescaled_e_gl(elec_num,4,elec_num,0:cord_num,walk_num)
@@ -7445,7 +7441,10 @@ function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_doc( &
 
   double precision  :: rij_inv(elec_num)
 
+
   allocate(elec_dist_gl(elec_num, 4, elec_num))
+  elec_dist_gl = 0.d0
+  een_rescaled_e_gl = 0.d0
 
   info = QMCKL_SUCCESS
 
@@ -7482,7 +7481,7 @@ function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_doc( &
         enddo
         do i = 1, elec_num
            do ii = 1, 3
-              elec_dist_gl(i, ii, j) = (coord_ee(i, ii, nw) - coord_ee(j, ii, nw)) * rij_inv(i)
+              elec_dist_gl(i, ii, j) = (coord_ee(i, nw, ii) - coord_ee(j, nw, ii)) * rij_inv(i)
            end do
            elec_dist_gl(i, 4, j) = 2.0d0 * rij_inv(i)
         end do
@@ -7492,27 +7491,20 @@ function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_doc( &
      een_rescaled_e_gl(:,:,:,0,nw) = 0.d0
 
      do l = 1, cord_num
-        kappa_l = - dble(l) * rescale_factor_ee
+        kappa_l = -dble(l) * rescale_factor_ee
         do j = 1, elec_num
            do i = 1, elec_num
-              een_rescaled_e_gl(i, 1, j, l, nw) = kappa_l * elec_dist_gl(i, 1, j)
-              een_rescaled_e_gl(i, 2, j, l, nw) = kappa_l * elec_dist_gl(i, 2, j)
-              een_rescaled_e_gl(i, 3, j, l, nw) = kappa_l * elec_dist_gl(i, 3, j)
-              een_rescaled_e_gl(i, 4, j, l, nw) = kappa_l * elec_dist_gl(i, 4, j)
-           end do
-
-           do i = 1, elec_num
-              een_rescaled_e_gl(i, 4, j, l, nw) = een_rescaled_e_gl(i, 4, j, l, nw)              &
-                   + een_rescaled_e_gl(i, 1, j, l, nw) * een_rescaled_e_gl(i, 1, j, l, nw)  &
-                   + een_rescaled_e_gl(i, 2, j, l, nw) * een_rescaled_e_gl(i, 2, j, l, nw)  &
-                   + een_rescaled_e_gl(i, 3, j, l, nw) * een_rescaled_e_gl(i, 3, j, l, nw)
-           end do
-
-           do i = 1, elec_num
-              een_rescaled_e_gl(i,1,j,l,nw) = een_rescaled_e_gl(i,1,j,l,nw) * een_rescaled_e(i,j,l,nw)
-              een_rescaled_e_gl(i,2,j,l,nw) = een_rescaled_e_gl(i,2,j,l,nw) * een_rescaled_e(i,j,l,nw)
-              een_rescaled_e_gl(i,3,j,l,nw) = een_rescaled_e_gl(i,3,j,l,nw) * een_rescaled_e(i,j,l,nw)
-              een_rescaled_e_gl(i,4,j,l,nw) = een_rescaled_e_gl(i,4,j,l,nw) * een_rescaled_e(i,j,l,nw)
+              if (i /= j) then
+                 een_rescaled_e_gl(i, 1, j, l, nw) = kappa_l *  elec_dist_gl(i, 1, j) * een_rescaled_e(i,j,l,nw)
+                 een_rescaled_e_gl(i, 2, j, l, nw) = kappa_l *  elec_dist_gl(i, 2, j) * een_rescaled_e(i,j,l,nw)
+                 een_rescaled_e_gl(i, 3, j, l, nw) = kappa_l *  elec_dist_gl(i, 3, j) * een_rescaled_e(i,j,l,nw)
+                 een_rescaled_e_gl(i, 4, j, l, nw) = kappa_l * (elec_dist_gl(i, 4, j) + kappa_l) * een_rescaled_e(i,j,l,nw)
+              else
+                 een_rescaled_e_gl(i, 1, j, l, nw) = 0.d0
+                 een_rescaled_e_gl(i, 2, j, l, nw) = 0.d0
+                 een_rescaled_e_gl(i, 3, j, l, nw) = 0.d0
+                 een_rescaled_e_gl(i, 4, j, l, nw) = 0.d0
+              end if
            end do
         end do
      end do
@@ -7599,10 +7591,10 @@ qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_hpc (const qmckl_context co
 #endif
     for (int64_t nw = 0; nw < walk_num; ++nw) {
 
-      double rij_inv[elec_num];
-
       for (int64_t j=0; j<elec_num; ++j) {
 
+        double rij_inv[elec_num];
+
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
@@ -7619,18 +7611,18 @@ qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_hpc (const qmckl_context co
         }
         rij_inv[j] = 0.0;
 
-        const double xj = coord_ee[j + nw * elec_num * 3];
-        const double yj = coord_ee[j + elec_num + nw * elec_num * 3];
-        const double zj = coord_ee[j + 2 * elec_num + nw * elec_num * 3];
+        const double xj = coord_ee[j+elec_num*nw];
+        const double yj = coord_ee[j+elec_num*(nw+walk_num) ];
+        const double zj = coord_ee[j+elec_num*(nw+walk_num*2)];
 
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
         for (int64_t i = 0; i < elec_num ; ++i) {
 
-          const double xi = coord_ee[i + nw * elec_num * 3];
-          const double yi = coord_ee[i + elec_num + nw * elec_num * 3];
-          const double zi = coord_ee[i + 2 * elec_num + nw * elec_num * 3];
+          const double xi = coord_ee[i+elec_num*nw];
+          const double yi = coord_ee[i+elec_num*(nw+walk_num) ];
+          const double zi = coord_ee[i+elec_num*(nw+walk_num*2)];
 
           elec_dist_gl0[i + j * elec_num] = rij_inv[i] * (xi-xj);
           elec_dist_gl1[i + j * elec_num] = rij_inv[i] * (yi-yj);
@@ -7639,47 +7631,18 @@ qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_hpc (const qmckl_context co
         }
       }
 
-      for (int64_t j = 0; j < elec_num; ++j) {
-
-        double* restrict eegl = &een_rescaled_e_gl[ elec_num * 4 * (j + elec_num * (cord_num + 1) * nw)];
-#ifdef HAVE_OPENMP
-#pragma omp simd
-#endif
-        for (int64_t i = 0; i < 4*elec_num; ++i) {
-          eegl[i] = 0.0;
-        }
-
-      }
+      double* restrict eegl = &een_rescaled_e_gl[elec_num*4*elec_num*(cord_num+1)*nw];
+      memset(eegl, 0, 4*elec_num*elec_num*sizeof(double));
 
       for (int64_t l=1; l<=cord_num; ++l) {
 
-        double kappa_l = - (double)l * rescale_factor_ee;
+        const double kappa_l  = -rescale_factor_ee * (double) l;
 
         for (int64_t j=0; j<elec_num; ++j) {
 
           double* restrict eegl =
             &een_rescaled_e_gl[elec_num*4*(j+elec_num*(l+(cord_num+1)*nw))];
 
-#ifdef HAVE_OPENMP
-#pragma omp simd
-#endif
-          for (int64_t i=0; i<elec_num; ++i) {
-            eegl[i               ] = kappa_l * elec_dist_gl0[i + j * elec_num];
-            eegl[i + elec_num    ] = kappa_l * elec_dist_gl1[i + j * elec_num];
-            eegl[i + elec_num * 2] = kappa_l * elec_dist_gl2[i + j * elec_num];
-            eegl[i + elec_num * 3] = kappa_l * elec_dist_gl3[i + j * elec_num];
-          }
-
-#ifdef HAVE_OPENMP
-#pragma omp simd
-#endif
-          for (int64_t i=0; i<elec_num; ++i) {
-            eegl[i + elec_num*3] = eegl[i + elec_num*3] +
-              eegl[i] * eegl[i] +
-              eegl[i + elec_num*1] * eegl[i + elec_num*1] +
-              eegl[i + elec_num*2] * eegl[i + elec_num*2];
-          }
-
           const double* restrict ee =
             &een_rescaled_e[elec_num*(j+elec_num*(l+(cord_num+1)*nw))];
 
@@ -7687,11 +7650,15 @@ qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_hpc (const qmckl_context co
 #pragma omp simd
 #endif
           for (int64_t i=0; i<elec_num; ++i) {
-            eegl[i               ] *= ee[i];
-            eegl[i + elec_num * 1] *= ee[i];
-            eegl[i + elec_num * 2] *= ee[i];
-            eegl[i + elec_num * 3] *= ee[i];
+            eegl[i           ] = kappa_l * ee[i] *  elec_dist_gl0[i+j*elec_num];
+            eegl[i+elec_num  ] = kappa_l * ee[i] *  elec_dist_gl1[i+j*elec_num];
+            eegl[i+elec_num*2] = kappa_l * ee[i] *  elec_dist_gl2[i+j*elec_num];
+            eegl[i+elec_num*3] = kappa_l * ee[i] * (elec_dist_gl3[i+j*elec_num] + kappa_l);
           }
+          eegl[j           ] = 0.0;
+          eegl[j+elec_num  ] = 0.0;
+          eegl[j+elec_num*2] = 0.0;
+          eegl[j+elec_num*3] = 0.0;
         }
       }
     }
@@ -7701,7 +7668,7 @@ qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_hpc (const qmckl_context co
     free(elec_dist_gl2);
     free(elec_dist_gl3);
   }
-
+  
   return QMCKL_SUCCESS;
 }
 #+end_src
@@ -7730,79 +7697,36 @@ qmckl_exit_code qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl (
 
 **** Test                                                          :noexport:
 
- #+name: een_e_gl
+     #+name: ee_dist_gl
      #+begin_src python :results output :exports none :noweb yes
-import numpy as np
+<<ee_dist>>
 
-<<jastrow_data>>
-
-elec_coord = np.array(elec_coord)[0]
-elec_dist = np.zeros(shape=(elec_num, elec_num),dtype=float)
+een_rescaled_e_gl = np.zeros(shape=(4, cord_num+1, elec_num, elec_num), dtype=float)
 for i in range(elec_num):
-  for j in range(elec_num):
-    elec_dist[i, j] = np.linalg.norm(elec_coord[i] - elec_coord[j])
-
-elec_dist_gl = np.zeros(shape=(4,elec_num, elec_num),dtype=float)
-for j in range(elec_num):
-  for i in range(elec_num):
-    rij_inv = 1.0 / elec_dist[i, j]
-    for ii in range(3):
-      elec_dist_gl[ii, i, j] = (elec_coord[i][ii] - elec_coord[j][ii]) * rij_inv
-    elec_dist_gl[3, i, j] = 2.0 * rij_inv
-  elec_dist_gl[:, j, j] = 0.0
-
-
-kappa = 0.6
-
-een_rescaled_e_ij = np.zeros(shape=(elec_num * (elec_num - 1)//2, cord_num+1), dtype=float)
-een_rescaled_e_ij[:,0] = 1.0
-
-k = 0
-for j in range(elec_num):
-  for i in range(j):
-    een_rescaled_e_ij[k, 1] = np.exp(-kappa * elec_dist[i, j])
-    k = k + 1
-
-for l in range(2, cord_num + 1):
-  for k in range(elec_num * (elec_num - 1)//2):
-    een_rescaled_e_ij[k, l] = een_rescaled_e_ij[k, l - 1] * een_rescaled_e_ij[k, 1]
-
-een_rescaled_e = np.zeros(shape=(elec_num, elec_num, cord_num + 1), dtype=float)
-een_rescaled_e[:,:,0] = 1.0
-
-for l in range(1,cord_num+1):
-  k = 0
-  for j in range(elec_num):
-    for i in range(j):
-      x = een_rescaled_e_ij[k, l]
-      een_rescaled_e[i, j, l] = x
-      een_rescaled_e[j, i, l] = x
-      k = k + 1
-
-een_rescaled_e_gl = np.zeros(shape=(elec_num,4,elec_num,cord_num+1),dtype=float)
-for l in range(0,cord_num+1):
-  kappa_l = -1.0 * kappa * l
-  for j in range(0,elec_num):
-    for i in range(0,elec_num):
-      for ii in range(0,4):
-        een_rescaled_e_gl[i,ii,j,l] = kappa_l * elec_dist_gl[ii,i,j]
-      een_rescaled_e_gl[i,3,j,l] = een_rescaled_e_gl[i,3,j,l] +        \
-      een_rescaled_e_gl[i,0,j,l] * een_rescaled_e_gl[i,0,j,l] +        \
-      een_rescaled_e_gl[i,1,j,l] * een_rescaled_e_gl[i,1,j,l] +        \
-      een_rescaled_e_gl[i,2,j,l] * een_rescaled_e_gl[i,2,j,l]
-
-      for ii in range(0,4):
-        een_rescaled_e_gl[i,ii,j,l] = een_rescaled_e_gl[i,ii,j,l] * een_rescaled_e[i,j,l]
-
-print(" een_rescaled_e_gl[1, 1, 3, 1] = ",een_rescaled_e_gl[0, 0, 2, 1])
-print(" een_rescaled_e_gl[1, 1, 4, 1] = ",een_rescaled_e_gl[0, 0, 3, 1])
-print(" een_rescaled_e_gl[1, 1, 5, 1] = ",een_rescaled_e_gl[0, 0, 4, 1])
-print(" een_rescaled_e_gl[2, 1, 4, 2] = ",een_rescaled_e_gl[1, 0, 3, 2])
-print(" een_rescaled_e_gl[2, 1, 5, 2] = ",een_rescaled_e_gl[1, 0, 4, 2])
-print(" een_rescaled_e_gl[2, 1, 6, 2] = ",een_rescaled_e_gl[1, 0, 5, 2])
+   for j in range(elec_num):
+      if j != i:
+         r = elec_dist[i,j] 
+         for p in range(cord_num+1):
+            f = een_rescaled_e[p,i,j]
+            for k in range(3):
+               een_rescaled_e_gl[k,p,i,j] = -kappa_e * p * (elec_coord[j,k] - elec_coord[i,k])/r * f
+            een_rescaled_e_gl[3,p,i,j] = kappa_e*p*(-2.0/r + kappa_e*p)*f
      #+end_src
 
-      #+begin_src c :tangle (eval c_test)
+     #+NAME:test_ee_gl
+     #+begin_src python :results output :exports none :noweb yes
+<<ee_dist_gl>>
+for p in range(cord_num+1):
+# for i in range(elec_num):
+  for i in range(3):
+    for k in range(4):
+#     for j in range(elec_num):
+      for j in range(3):
+        print(f"printf( \"%e %e\\n\", een_rescaled_e_gl[0][{p}][{i}][{k}][{j}], {een_rescaled_e_gl[k,p,i,j]});")
+        print(f"assert( fabs(een_rescaled_e_gl[0][{p}][{i}][{k}][{j}] - ({een_rescaled_e_gl[k,p,i,j]})) < 1.e-10 );")
+     #+end_src
+
+      #+begin_src c :tangle (eval c_test) :noweb yes 
 assert(qmckl_electron_provided(context));
 
 {
@@ -7810,13 +7734,7 @@ assert(qmckl_electron_provided(context));
   size_max=walk_num*(cord_num + 1)*elec_num*4*elec_num;
   rc = qmckl_get_jastrow_champ_een_distance_rescaled_e_gl(context,
             &(een_rescaled_e_gl[0][0][0][0][0]),size_max);
-
-  assert(fabs(een_rescaled_e_gl[0][1][0][0][2] +  0.09831391870751387   ) < 1.e-8);
-  assert(fabs(een_rescaled_e_gl[0][1][0][0][3] +  0.017204157459682526  ) < 1.e-8);
-  assert(fabs(een_rescaled_e_gl[0][1][0][0][4] +  0.013345768421098641  ) < 1.e-8);
-  assert(fabs(een_rescaled_e_gl[0][2][1][0][3] +  0.03733086358273962   ) < 1.e-8);
-  assert(fabs(een_rescaled_e_gl[0][2][1][0][4] +  0.004922634822943517  ) < 1.e-8);
-  assert(fabs(een_rescaled_e_gl[0][2][1][0][5] +  0.5416751547830984    ) < 1.e-8);
+  <<test_ee_gl()>>
 }
 
 
@@ -7848,19 +7766,27 @@ assert(qmckl_electron_provided(context));
             een_rescaled_e_gl_hpc);
   assert(rc == QMCKL_SUCCESS);
 
-  for (int64_t i = 0; i < walk_num*(cord_num + 1)*elec_num*4*elec_num; i++) {
-    if (fabs(een_rescaled_e_gl_doc[i] - een_rescaled_e_gl_hpc[i]) > 1.e-12) {
-      printf("i = %ld, doc = %f, hpc = %f\n", i, een_rescaled_e_gl_doc[i], een_rescaled_e_gl_hpc[i]);
-      fflush(stdout);
+  for (int nw=0 ; nw < walk_num ; nw++) {
+    for (int c=0 ; c <= cord_num ; c++) {
+      for (int i=0 ; i < elec_num ; i++) {
+        for (int j=0 ; j < elec_num ; j++) {
+          for (int k=0 ; k < 4 ; k++) {
+              printf("nw=%d c=%d i=%d k=%d j=%d doc=%e hpc=%e\n", nw, c, i, k, j, een_rescaled_e_gl_doc[nw*(cord_num + 1)*elec_num*4*elec_num + c*elec_num*4*elec_num + i*4*elec_num + k*elec_num + j], een_rescaled_e_gl_hpc[nw*(cord_num + 1)*elec_num*4*elec_num + c*elec_num*4*elec_num + i*4*elec_num + k*elec_num + j]);
+            if (fabs(een_rescaled_e_gl_doc[nw*(cord_num + 1)*elec_num*4*elec_num + c*elec_num*4*elec_num + i*4*elec_num + k*elec_num + j] - een_rescaled_e_gl_hpc[nw*(cord_num + 1)*elec_num*4*elec_num + c*elec_num*4*elec_num + i*4*elec_num + k*elec_num + j]) > 1.e-12) {
+              printf("nw=%d c=%d i=%d k=%d j=%d doc=%e hpc=%e\n", nw, c, i, k, j, een_rescaled_e_gl_doc[nw*(cord_num + 1)*elec_num*4*elec_num + c*elec_num*4*elec_num + i*4*elec_num + k*elec_num + j], een_rescaled_e_gl_hpc[nw*(cord_num + 1)*elec_num*4*elec_num + c*elec_num*4*elec_num + i*4*elec_num + k*elec_num + j]);
+              fflush(stdout);
+            }
+            assert(fabs(een_rescaled_e_gl_doc[nw*(cord_num + 1)*elec_num*4*elec_num + c*elec_num*4*elec_num + i*4*elec_num + k*elec_num + j] - een_rescaled_e_gl_hpc[nw*(cord_num + 1)*elec_num*4*elec_num + c*elec_num*4*elec_num + i*4*elec_num + k*elec_num + j]) < 1.e-8);
+          }
+        }
+      }
     }
-    assert(fabs(een_rescaled_e_gl_doc[i] - een_rescaled_e_gl_hpc[i]) < 1.e-8);
   }
 
 }
 
 {
   /* Finite difference test fails and I can't understand why... */
-/*
 
   printf("een_distance_rescaled_e_gl\n");
 
@@ -7876,11 +7802,11 @@ assert(qmckl_electron_provided(context));
 
   qmckl_exit_code rc;
 
-  double elec_coord[walk_num][3][elec_num];
-  rc = qmckl_get_electron_coord (context, 'T', &(elec_coord[0][0][0]), 3*walk_num*elec_num);
+  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);
   assert (rc == QMCKL_SUCCESS);
 
-  double temp_coord[walk_num][3][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][cord_num+1][elec_num][elec_num];
@@ -7890,12 +7816,13 @@ assert(qmckl_electron_provided(context));
     for (int64_t k = 0; k < 3; k++) {
       for (int64_t m = -4; m <= 4; m++) {          // Apply finite difference displacement
 
+        memcpy(&(temp_coord[0][0][0]), &(elec_coord[0][0][0]), sizeof(temp_coord));
         for (int64_t nw=0 ; nw<walk_num ; nw++) {
-          temp_coord[nw][k][i] = elec_coord[nw][k][i] + (double) m * delta_x;
+          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, 'T', walk_num,
+        rc = qmckl_set_electron_coord (context, 'N', walk_num,
                                        &(temp_coord[0][0][0]),
                                        walk_num*3*elec_num);
         assert(rc == QMCKL_SUCCESS);
@@ -7917,14 +7844,11 @@ assert(qmckl_electron_provided(context));
         }
 
       }
-      for (int64_t nw=0 ; nw<walk_num ; nw++) {
-        temp_coord[nw][k][i] = elec_coord[nw][k][i];
-      }
     }
   }
 
   // Reset coordinates in the context
-  rc = qmckl_set_electron_coord (context, 'T', walk_num,
+  rc = qmckl_set_electron_coord (context, 'N', walk_num,
                                  &(elec_coord[0][0][0]),
                                  walk_num*3*elec_num);
   assert(rc == QMCKL_SUCCESS);
@@ -7953,6 +7877,7 @@ assert(qmckl_electron_provided(context));
                                                                       walk_num*(cord_num+1)*elec_num*4*elec_num)
                    );
 
+
   assert(rc == QMCKL_SUCCESS);
 
   for (int nw = 0; nw < walk_num; nw++){
@@ -7983,9 +7908,7 @@ assert(qmckl_electron_provided(context));
     }
   }
   printf("OK\n");
-*/
 }
-
       #+end_src
 
 *** Electron-nucleus rescaled distances in $J_\text{eeN}$
@@ -7994,7 +7917,7 @@ assert(qmckl_electron_provided(context));
     electrons and nuclei raised to the power \(p\) defined by ~cord_num~:
 
     \[
-    C_{i\alpha,p} = \left[ \exp\left(-\kappa_\alpha\, R_{i\alpha}\right) \right]^p
+    [g_{\alpha}(R_{i\alpha})]^p = \exp\left(-p\, \kappa_\alpha\, R_{i\alpha}\right) 
     \]
 
    where \(R_{i\alpha}\) is the matrix of electron-nucleus distances.
@@ -8741,6 +8664,20 @@ assert(fabs( 0.005359281880312882  - een_rescaled_n_gl[0][2][1][0][5])  < 1.e-12
    calculation of the three-body jastrow ~factor_een~ and its derivative
    ~factor_een_gl~.
 
+   The array ~tmp_c~ corresponds to the tensor $P$ defined at the
+   beginning of this section:
+
+   \[ P_{i\alpha}^{km} = 
+      \sum_{j=1}^{N_{\text{elec}}} 
+             \left[ g_\text{e}({r}_{ij}) \right]^k
+             \left[ g_\alpha({R}_{j\alpha}) \right]^{m}
+             \]
+
+   \[ \nabla_i P_{i\alpha}^{km} = 
+      \sum_{j=1}^{N_{\text{elec}}} 
+             \nabla_i \left[ g_\text{e}({r}_{ij}) \right]^k
+             \left[ g_\alpha({R}_{j\alpha}) \right]^{m}
+             \]
 **** Compute dim_c_vector
      :PROPERTIES:
      :Name:     qmckl_compute_dim_c_vector
@@ -9693,7 +9630,7 @@ qmckl_exit_code qmckl_compute_tmp_c_doc (
           double* const tmp_c );
      #+end_src
 
-***** CPU                                                          :noexport:
+***** HPC                                                          :noexport:
 
      #+begin_src c :comments org :tangle (eval c) :noweb yes
 qmckl_exit_code qmckl_compute_tmp_c_hpc (
@@ -9802,16 +9739,18 @@ qmckl_exit_code qmckl_compute_tmp_c_hpc (const qmckl_context context,
      :END:
 
      #+NAME: qmckl_factor_dtmp_c_args
-     | Variable            | Type                                                             | In/Out | Description                                   |
-     |---------------------+------------------------------------------------------------------+--------+-----------------------------------------------|
-     | ~context~           | ~qmckl_context~                                                  | in     | Global state                                  |
-     | ~cord_num~          | ~int64_t~                                                        | in     | Order of polynomials                          |
-     | ~elec_num~          | ~int64_t~                                                        | in     | Number of electrons                           |
-     | ~nucl_num~          | ~int64_t~                                                        | in     | Number of nuclei                              |
-     | ~walk_num~          | ~int64_t~                                                        | in     | Number of walkers                             |
-     | ~een_rescaled_e_gl~ | ~double[walk_num][0:cord_num][elec_num][4][elec_num]~            | in     | Electron-electron rescaled factor derivatives |
-     | ~een_rescaled_n~    | ~double[walk_num][0:cord_num][nucl_num][elec_num]~               | in     | Electron-nucleus rescaled factor              |
-     | ~dtmp_c~            | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | out    | vector of non-zero coefficients               |
+     |---------------------+---------------------------------------------------------------------+--------+-----------------------------------------------|
+     | Variable            | Type                                                                | In/Out | Description                                   |
+     |---------------------+---------------------------------------------------------------------+--------+-----------------------------------------------|
+     | ~context~           | ~qmckl_context~                                                     | in     | Global state                                  |
+     | ~cord_num~          | ~int64_t~                                                           | in     | Order of polynomials                          |
+     | ~elec_num~          | ~int64_t~                                                           | in     | Number of electrons                           |
+     | ~nucl_num~          | ~int64_t~                                                           | in     | Number of nuclei                              |
+     | ~walk_num~          | ~int64_t~                                                           | in     | Number of walkers                             |
+     | ~een_rescaled_e_gl~ | ~double[walk_num][0:cord_num][elec_num][4][elec_num]~               | in     | Electron-electron rescaled factor derivatives |
+     | ~een_rescaled_n~    | ~double[walk_num][0:cord_num][nucl_num][elec_num]~                  | in     | Electron-nucleus rescaled factor              |
+     | ~dtmp_c~            | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]~ | out    | vector of non-zero coefficients               |
+     |---------------------+---------------------------------------------------------------------+--------+-----------------------------------------------|
 
 
      #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
@@ -9991,85 +9930,58 @@ qmckl_exit_code qmckl_compute_dtmp_c_hpc (
 import numpy as np
 
 <<jastrow_data>>
-
-elec_coord = np.array(elec_coord)[0]
-nucl_coord = np.array(nucl_coord)
-elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float)
-for i in range(elec_num):
-  for a in range(nucl_num):
-    elnuc_dist[i, a] = np.linalg.norm(elec_coord[i] - nucl_coord[:,a])
-
-kappa = 0.6
-
-een_rescaled_n = np.zeros(shape=(nucl_num, elec_num, cord_num + 1), dtype=float)
-een_rescaled_n[:,:,0] = 1.0
-
-for a in range(nucl_num):
-  for i in range(elec_num):
-    een_rescaled_n[a, i, 1] = np.exp(-kappa * elnuc_dist[i, a])
-
-for l in range(2,cord_num+1):
-  for a in range(nucl_num):
-    for i in range(elec_num):
-      een_rescaled_n[a, i, l] = een_rescaled_n[a, i, l - 1] * een_rescaled_n[a, i, 1]
-
-elec_dist = np.zeros(shape=(elec_num, elec_num),dtype=float)
-for i in range(elec_num):
-  for j in range(elec_num):
-    elec_dist[i, j] = np.linalg.norm(elec_coord[i] - elec_coord[j])
-
-kappa = 0.6
-
-een_rescaled_e_ij = np.zeros(shape=(elec_num * (elec_num - 1)//2, cord_num+1), dtype=float)
-een_rescaled_e_ij[:,0] = 1.0
-
-k = 0
-for j in range(elec_num):
-  for i in range(j):
-    een_rescaled_e_ij[k, 1] = np.exp(-kappa * elec_dist[i, j])
-    k = k + 1
-
-for l in range(2, cord_num + 1):
-  for k in range(elec_num * (elec_num - 1)//2):
-    een_rescaled_e_ij[k, l] = een_rescaled_e_ij[k, l - 1] * een_rescaled_e_ij[k, 1]
-
-een_rescaled_e = np.zeros(shape=(elec_num, elec_num, cord_num + 1), dtype=float)
-een_rescaled_e[:,:,0] = 1.0
-
-for l in range(1,cord_num+1):
-  k = 0
-  for j in range(elec_num):
-    for i in range(j):
-      x = een_rescaled_e_ij[k, l]
-      een_rescaled_e[i, j, l] = x
-      een_rescaled_e[j, i, l] = x
-      k = k + 1
-
-for l in range(0,cord_num+1):
-  for j in range(0, elec_num):
-    een_rescaled_e[j,j,l] = 0.0
-
 lkpm_of_cindex = np.array(lkpm_combined_index).T
+<<en_dist>>
+<<ee_dist_gl>>
      #+end_src
 
-     #+RESULTS: helper_funcs
 
-      #+begin_src c :tangle (eval c_test)
-assert(qmckl_electron_provided(context));
+     #+NAME: test_tmpc
+     #+begin_src python :results output :exports none :noweb yes
+<<helper_funcs>>
+tmp_c = np.einsum("pij,qaj->pqai", een_rescaled_e, een_rescaled_n)
+for p in range(cord_num):
+  for q in range(cord_num+1):
+   for a in range(nucl_num):
+#     for i in range(elec_num):
+      for i in range(3):
+         print(f"assert( fabs(tmp_c[0][{p}][{q}][{a}][{i}] - ({tmp_c[p,q,a,i]})) < 1.e-10 );")
+     #+end_src
 
-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]), sizeof(tmp_c)/sizeof(double));
+     #+NAME: test_dtmpc
+     #+begin_src python :results output :exports none :noweb yes
+<<helper_funcs>>
+dtmp_c = np.einsum("lpji,qaj->pqali", een_rescaled_e_gl, een_rescaled_n)
+for p in range(cord_num):
+  for q in range(cord_num+1):
+    for a in range(nucl_num):
+      for l in range(4):
+#       for i in range(elec_num):
+        for i in range(3):
+          print(f"printf(\"%e %e\\n\", dtmp_c[0][{p}][{q}][{a}][{l}][{i}], {dtmp_c[p,q,a,l,i]}); fflush(stdout);")
+          print(f"assert( fabs(dtmp_c[0][{p}][{q}][{a}][{l}][{i}] - ({dtmp_c[p,q,a,l,i]})) < 1.e-10 );")
+     #+end_src
 
-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]), sizeof(dtmp_c)/sizeof(double));
+     #+RESULTS: test_dtmpc
 
-printf("%e\n%e\n", tmp_c[0][0][1][0][0],  3.954384);
-fflush(stdout);
-assert(fabs(tmp_c[0][0][1][0][0] - 3.954384) < 1e-6);
+      #+begin_src c :tangle (eval c_test) :noweb yes
+{
+  assert(qmckl_electron_provided(context));
 
-printf("%e\n%e\n", dtmp_c[0][1][0][0][0][0],3.278657e-01);
-fflush(stdout);
-assert(fabs(dtmp_c[0][1][0][0][0][0] - 3.278657e-01 ) < 1e-6);
+  printf("tmp_c\n");
+  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]), sizeof(tmp_c)/sizeof(double));
+
+  <<test_tmpc()>>
+}
+
+{
+  printf("dtmp_c\n");
+  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]), sizeof(dtmp_c)/sizeof(double));
+
+  <<test_dtmpc()>>
+}
       #+end_src
 
 *** Electron-electron-nucleus Jastrow $f_{een}$
@@ -11386,8 +11298,6 @@ import numpy as np
 
 <<jastrow_data>>
 
-<<een_e_gl>>
-
 <<helper_funcs>>
 
 kappa = 0.6
@@ -12989,7 +12899,7 @@ assert(qmckl_jastrow_champ_provided(context));
 
 }
        #+end_src
-
+       
 * End of files                                                     :noexport:
 
   #+begin_src c :tangle (eval h_private_type)