diff --git a/org/qmckl_jastrow_champ_single.org b/org/qmckl_jastrow_champ_single.org
index 86bed1e..f74e524 100644
--- a/org/qmckl_jastrow_champ_single.org
+++ b/org/qmckl_jastrow_champ_single.org
@@ -3669,7 +3669,7 @@ qmckl_exit_code qmckl_provide_jastrow_champ_single_een_gl(qmckl_context context)
       ctx->single_point.delta_een_gl = delta_een_gl;
     }
 
-    rc = qmckl_compute_jastrow_champ_factor_single_een_gl_doc(context,
+    rc = qmckl_compute_jastrow_champ_factor_single_een_gl(context,
                                                 ctx->single_point.num,
                                                 ctx->electron.walker.num,
                                                 ctx->electron.num,
@@ -3756,126 +3756,8 @@ integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_factor_single_een_
   double precision :: accu, accu2, cn
   integer*8                           :: LDA, LDB, LDC
 
-  double precision  :: een_rescaled_delta_n_gl(4, nucl_num, 0:cord_num, walk_num)
-  double precision  :: een_rescaled_delta_n(nucl_num, 0:cord_num, walk_num)
-double precision :: dpg1_m, dpg1_ml, dp_m, dp_ml, een_r_m, een_r_ml, een_r_gl_m, een_r_gl_ml
-  num = num_in + 1
-
-  info = QMCKL_SUCCESS
-
-  if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT
-  if (walk_num <= 0)                 info = QMCKL_INVALID_ARG_3
-  if (elec_num <= 0)                 info = QMCKL_INVALID_ARG_4
-  if (nucl_num <= 0)                 info = QMCKL_INVALID_ARG_5
-  if (cord_num <  0)                 info = QMCKL_INVALID_ARG_6
-  if (info /= QMCKL_SUCCESS)         return
-
-  delta_een_gl = 0.0d0
-
-  if (cord_num == 0) return
-
-  een_rescaled_delta_n(:,:,:) = een_rescaled_single_n(:,:,:) - een_rescaled_n(num, :, :, :)
-  een_rescaled_delta_n_gl(:,:,:,:) = een_rescaled_single_n_gl(:,:,:,:) - een_rescaled_n_gl(num, :,:,:,:)
-
-
-  do nw =1, walk_num
-     do n = 1, dim_c_vector
-        l = lkpm_combined_index(n, 1)
-        k = lkpm_combined_index(n, 2)
-        p = lkpm_combined_index(n, 3)
-        m = lkpm_combined_index(n, 4)
-
-        do kk = 1, 4
-           do a = 1, nucl_num
-              cn = c_vector_full(a, n)
-              if(cn == 0.d0) cycle
-              !do i = 1, elec_num
-              !   delta_een_gl(i,kk,nw) = delta_een_gl(i,kk,nw) + ( &
-              !        delta_p_gl(i,a,kk,m  ,k,nw) * een_rescaled_n(i,a,m+l,nw) + &
-              !        delta_p_gl(i,a,kk,m+l,k,nw) * een_rescaled_n(i,a,m  ,nw) + &
-              !        delta_p(i,a,m  ,k,nw) * een_rescaled_n_gl(i,kk,a,m+l,nw) + &
-              !        delta_p(i,a,m+l,k,nw) * een_rescaled_n_gl(i,kk,a,m  ,nw) ) * cn
-              !end do
-              do i = 1, elec_num
-    ! Cache repeated accesses
-    dpg1_m   = delta_p_gl(i,a,kk,m  ,k,nw)
-    dpg1_ml  = delta_p_gl(i,a,kk,m+l,k,nw)
-    dp_m     = delta_p(i,a,m  ,k,nw)
-    dp_ml    = delta_p(i,a,m+l,k,nw)
-    
-    een_r_m  = een_rescaled_n(i,a,m  ,nw)
-    een_r_ml = een_rescaled_n(i,a,m+l,nw)
-    een_r_gl_m  = een_rescaled_n_gl(i,kk,a,m  ,nw)
-    een_r_gl_ml = een_rescaled_n_gl(i,kk,a,m+l,nw)
-
-    delta_een_gl(i,kk,nw) = delta_een_gl(i,kk,nw) + cn * &
-        (dpg1_m * een_r_ml + dpg1_ml * een_r_m + dp_m * een_r_gl_ml + dp_ml * een_r_gl_m)
-end do
-
-              delta_een_gl(num,kk,nw) = delta_een_gl(num,kk,nw) + ( &
-                   (dtmp_c(num,kk,a,m  ,k,nw) + delta_p_gl(num,a,kk,m  ,k,nw)) * een_rescaled_delta_n(a,m+l,nw) + &
-                   (dtmp_c(num,kk,a,m+l,k,nw) + delta_p_gl(num,a,kk,m+l,k,nw)) * een_rescaled_delta_n(a,m  ,nw) + &
-                   (tmp_c(num,a,m  ,k,nw) + delta_p(num,a,m  ,k,nw)) * een_rescaled_delta_n_gl(kk,a,m+l,nw)  + &
-                   (tmp_c(num,a,m+l,k,nw) + delta_p(num,a,m+l,k,nw)) * een_rescaled_delta_n_gl(kk,a,m  ,nw) )* cn
-           end do
-        end do
-        do a = 1, nucl_num
-           cn = c_vector_full(a, n)
-           if(cn == 0.d0) cycle
-           cn = cn + cn
-           do i = 1, elec_num
-              delta_een_gl(i,4,nw) = delta_een_gl(i,4,nw) + ( &
-                   delta_p_gl(i,a,1,m  ,k,nw) * een_rescaled_n_gl(i,1,a,m+l,nw) + &
-                   delta_p_gl(i,a,1,m+l,k,nw) * een_rescaled_n_gl(i,1,a,m  ,nw) + &
-                   delta_p_gl(i,a,2,m  ,k,nw) * een_rescaled_n_gl(i,2,a,m+l,nw) + &
-                   delta_p_gl(i,a,2,m+l,k,nw) * een_rescaled_n_gl(i,2,a,m  ,nw) + &
-                   delta_p_gl(i,a,3,m  ,k,nw) * een_rescaled_n_gl(i,3,a,m+l,nw) + &
-                   delta_p_gl(i,a,3,m+l,k,nw) * een_rescaled_n_gl(i,3,a,m  ,nw) ) * cn
-           end do
-           delta_een_gl(num,4,nw) = delta_een_gl(num,4,nw) + ( &
-                (delta_p_gl(num,a,1,m  ,k,nw) + dtmp_c(num,1,a,m  ,k,nw)) * een_rescaled_delta_n_gl(1,a,m+l,nw) + &
-                (delta_p_gl(num,a,1,m+l,k,nw) + dtmp_c(num,1,a,m+l,k,nw)) * een_rescaled_delta_n_gl(1,a,m  ,nw) + &
-                (delta_p_gl(num,a,2,m  ,k,nw) + dtmp_c(num,2,a,m  ,k,nw)) * een_rescaled_delta_n_gl(2,a,m+l,nw) + &
-                (delta_p_gl(num,a,2,m+l,k,nw) + dtmp_c(num,2,a,m+l,k,nw)) * een_rescaled_delta_n_gl(2,a,m  ,nw) + &
-                (delta_p_gl(num,a,3,m  ,k,nw) + dtmp_c(num,3,a,m  ,k,nw)) * een_rescaled_delta_n_gl(3,a,m+l,nw) + &
-                (delta_p_gl(num,a,3,m+l,k,nw) + dtmp_c(num,3,a,m+l,k,nw)) * een_rescaled_delta_n_gl(3,a,m  ,nw) ) * cn
-        end do
-     end do
-  end do
-
-end function qmckl_compute_jastrow_champ_factor_single_een_gl_doc
-     #+end_src
-
-     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_factor_single_een_gl_hpc( &
-     context, num_in, walk_num, elec_num, nucl_num, cord_num,   &
-     dim_c_vector, c_vector_full, lkpm_combined_index, &
-     tmp_c, dtmp_c, delta_p, delta_p_gl, een_rescaled_n, een_rescaled_single_n, &
-     een_rescaled_n_gl, een_rescaled_single_n_gl, delta_een_gl) &
-     result(info) bind(C)
-  use, intrinsic :: iso_c_binding
-  use qmckl
-  implicit none
-  integer(qmckl_context), intent(in)  :: context
-  integer(c_int64_t)    , intent(in), value  :: num_in, walk_num, elec_num, cord_num, nucl_num, dim_c_vector
-  integer(c_int64_t)    , intent(in)  :: lkpm_combined_index(dim_c_vector,4)
-  real(c_double)        , intent(in)  :: c_vector_full(nucl_num, dim_c_vector)
-  real(c_double)        , intent(in)  :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
-  real(c_double)        , intent(in)  :: dtmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
-  real(c_double)        , intent(in)  :: delta_p(elec_num, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
-  real(c_double)        , intent(in)  :: delta_p_gl(elec_num, nucl_num, 4, 0:cord_num, 0:cord_num-1,  walk_num)
-  real(c_double)        , intent(in)  :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
-  real(c_double)        , intent(in)  :: een_rescaled_single_n(nucl_num, 0:cord_num, walk_num)
-  real(c_double)        , intent(in)  :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num)
-  real(c_double)        , intent(in)  :: een_rescaled_single_n_gl(4, nucl_num, 0:cord_num, walk_num)
-  real(c_double)        , intent(out) :: delta_een_gl(elec_num, 4, walk_num)
-
-  integer*8 :: i, a, j, l, k, p, m, n, nw, kk, num
-  double precision :: accu, accu2, cn
-  integer*8                           :: LDA, LDB, LDC
-
-  double precision  :: een_rescaled_delta_n_gl(4, nucl_num, 0:cord_num, walk_num)
-  double precision  :: een_rescaled_delta_n(nucl_num, 0:cord_num, walk_num)
+  double precision  :: een_rescaled_delta_n(nucl_num, 0:cord_num)
+  double precision  :: een_rescaled_delta_n_gl(4, nucl_num, 0:cord_num)
   double precision :: dpg1_m, dpg1_ml, dp_m, dp_ml, een_r_m, een_r_ml, een_r_gl_m, een_r_gl_ml
   num = num_in + 1
 
@@ -3892,11 +3774,12 @@ integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_factor_single_een_
 
   if (cord_num == 0) return
 
-  een_rescaled_delta_n(:,:,:) = een_rescaled_single_n(:,:,:) - een_rescaled_n(num, :, :, :)
-  een_rescaled_delta_n_gl(:,:,:,:) = een_rescaled_single_n_gl(:,:,:,:) - een_rescaled_n_gl(num, :,:,:,:)
-
 
   do nw =1, walk_num
+
+     een_rescaled_delta_n(:,:) = een_rescaled_single_n(:,:,nw) - een_rescaled_n(num, :, :, nw)
+     een_rescaled_delta_n_gl(:,:,:) = een_rescaled_single_n_gl(:,:,:,nw) - een_rescaled_n_gl(num, :,:,:,nw)
+
      do n = 1, dim_c_vector
         l = lkpm_combined_index(n, 1)
         k = lkpm_combined_index(n, 2)
@@ -3931,10 +3814,10 @@ integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_factor_single_een_
               end do
 
               delta_een_gl(num,kk,nw) = delta_een_gl(num,kk,nw) + ( &
-                   (dtmp_c(num,kk,a,m  ,k,nw) + delta_p_gl(num,a,kk,m  ,k,nw)) * een_rescaled_delta_n(a,m+l,nw) + &
-                   (dtmp_c(num,kk,a,m+l,k,nw) + delta_p_gl(num,a,kk,m+l,k,nw)) * een_rescaled_delta_n(a,m  ,nw) + &
-                   (tmp_c(num,a,m  ,k,nw) + delta_p(num,a,m  ,k,nw)) * een_rescaled_delta_n_gl(kk,a,m+l,nw)  + &
-                   (tmp_c(num,a,m+l,k,nw) + delta_p(num,a,m+l,k,nw)) * een_rescaled_delta_n_gl(kk,a,m  ,nw) )* cn
+                   (dtmp_c(num,kk,a,m  ,k,nw) + delta_p_gl(num,a,kk,m  ,k,nw)) * een_rescaled_delta_n(a,m+l) + &
+                   (dtmp_c(num,kk,a,m+l,k,nw) + delta_p_gl(num,a,kk,m+l,k,nw)) * een_rescaled_delta_n(a,m  ) + &
+                   (tmp_c(num,a,m  ,k,nw) + delta_p(num,a,m  ,k,nw)) * een_rescaled_delta_n_gl(kk,a,m+l)  + &
+                   (tmp_c(num,a,m+l,k,nw) + delta_p(num,a,m+l,k,nw)) * een_rescaled_delta_n_gl(kk,a,m  ) )* cn
            end do
         end do
         do a = 1, nucl_num
@@ -3950,18 +3833,163 @@ integer(qmckl_exit_code) function qmckl_compute_jastrow_champ_factor_single_een_
                    delta_p_gl(i,a,3,m  ,k,nw) * een_rescaled_n_gl(i,3,a,m+l,nw) + &
                    delta_p_gl(i,a,3,m+l,k,nw) * een_rescaled_n_gl(i,3,a,m  ,nw) ) * cn
            end do
-           delta_een_gl(num,4,nw) = delta_een_gl(num,4,nw) + ( &
-                (delta_p_gl(num,a,1,m  ,k,nw) + dtmp_c(num,1,a,m  ,k,nw)) * een_rescaled_delta_n_gl(1,a,m+l,nw) + &
-                (delta_p_gl(num,a,1,m+l,k,nw) + dtmp_c(num,1,a,m+l,k,nw)) * een_rescaled_delta_n_gl(1,a,m  ,nw) + &
-                (delta_p_gl(num,a,2,m  ,k,nw) + dtmp_c(num,2,a,m  ,k,nw)) * een_rescaled_delta_n_gl(2,a,m+l,nw) + &
-                (delta_p_gl(num,a,2,m+l,k,nw) + dtmp_c(num,2,a,m+l,k,nw)) * een_rescaled_delta_n_gl(2,a,m  ,nw) + &
-                (delta_p_gl(num,a,3,m  ,k,nw) + dtmp_c(num,3,a,m  ,k,nw)) * een_rescaled_delta_n_gl(3,a,m+l,nw) + &
-                (delta_p_gl(num,a,3,m+l,k,nw) + dtmp_c(num,3,a,m+l,k,nw)) * een_rescaled_delta_n_gl(3,a,m  ,nw) ) * cn
+          delta_een_gl(num,4,nw) = delta_een_gl(num,4,nw) + ( &
+               (delta_p_gl(num,a,1,m  ,k,nw) + dtmp_c(num,1,a,m  ,k,nw)) * een_rescaled_delta_n_gl(1,a,m+l) + &
+               (delta_p_gl(num,a,1,m+l,k,nw) + dtmp_c(num,1,a,m+l,k,nw)) * een_rescaled_delta_n_gl(1,a,m  ) + &
+               (delta_p_gl(num,a,2,m  ,k,nw) + dtmp_c(num,2,a,m  ,k,nw)) * een_rescaled_delta_n_gl(2,a,m+l) + &
+               (delta_p_gl(num,a,2,m+l,k,nw) + dtmp_c(num,2,a,m+l,k,nw)) * een_rescaled_delta_n_gl(2,a,m  ) + &
+               (delta_p_gl(num,a,3,m  ,k,nw) + dtmp_c(num,3,a,m  ,k,nw)) * een_rescaled_delta_n_gl(3,a,m+l) + &
+               (delta_p_gl(num,a,3,m+l,k,nw) + dtmp_c(num,3,a,m+l,k,nw)) * een_rescaled_delta_n_gl(3,a,m  ) ) * cn
         end do
      end do
   end do
 
-end function qmckl_compute_jastrow_champ_factor_single_een_gl_hpc
+end function qmckl_compute_jastrow_champ_factor_single_een_gl_doc
+     #+end_src
+
+     #+begin_src c :comments org :tangle (eval c) :noweb yes
+qmckl_exit_code
+qmckl_compute_jastrow_champ_factor_single_een_gl_hpc (const qmckl_context context,
+                                                      const int64_t num,
+                                                      const int64_t walk_num,
+                                                      const int64_t elec_num,
+                                                      const int64_t nucl_num,
+                                                      const int64_t cord_num,
+                                                      const int64_t dim_c_vector,
+                                                      const double* c_vector_full,
+                                                      const int64_t* lkpm_combined_index,
+                                                      const double* tmp_c,
+                                                      const double* dtmp_c,
+                                                      const double* delta_p,
+                                                      const double* delta_p_gl,
+                                                      const double* een_rescaled_n,
+                                                      const double* een_rescaled_single_n,
+                                                      const double* een_rescaled_n_gl,
+                                                      const double* een_rescaled_single_n_gl,
+                                                      double* const delta_een_gl )
+{
+
+
+  if (context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT;
+  if (walk_num <= 0)                 return QMCKL_INVALID_ARG_3;
+  if (elec_num <= 0)                 return QMCKL_INVALID_ARG_4;
+  if (nucl_num <= 0)                 return QMCKL_INVALID_ARG_5;
+  if (cord_num <  0)                 return QMCKL_INVALID_ARG_6;
+
+  if (cord_num == 0) {
+    #pragma omp parallel for
+    for (size_t i=0 ; i<walk_num*4*elec_num ; ++i) {
+      delta_een_gl[i] = 0.;
+    }
+    return QMCKL_SUCCESS;
+  }
+
+#pragma omp parallel for
+  for (int64_t nw=0 ; nw<walk_num ; nw++) {
+    for (size_t i=0 ; i<4*elec_num ; ++i) {
+      delta_een_gl[i+nw*4*elec_num] = 0.;
+    }
+    double een_rescaled_delta_n[nucl_num*(cord_num+1)];
+    double een_rescaled_delta_n_gl[4*nucl_num*(cord_num+1)];
+    for (int64_t k=0 ; k<nucl_num*(cord_num+1) ; ++k) {
+      een_rescaled_delta_n[k] = een_rescaled_single_n[k+nucl_num*(cord_num+1)*nw] - een_rescaled_n[num+elec_num*(k+nucl_num*(cord_num+1)*nw)];
+      een_rescaled_delta_n_gl[0+4*k] = een_rescaled_single_n_gl[0+4*(k+nucl_num*(cord_num+1)*nw)] - een_rescaled_n_gl[num+elec_num*(0+4*(k+nucl_num*(cord_num+1)*nw))];
+      een_rescaled_delta_n_gl[1+4*k] = een_rescaled_single_n_gl[1+4*(k+nucl_num*(cord_num+1)*nw)] - een_rescaled_n_gl[num+elec_num*(1+4*(k+nucl_num*(cord_num+1)*nw))];
+      een_rescaled_delta_n_gl[2+4*k] = een_rescaled_single_n_gl[2+4*(k+nucl_num*(cord_num+1)*nw)] - een_rescaled_n_gl[num+elec_num*(2+4*(k+nucl_num*(cord_num+1)*nw))];
+      een_rescaled_delta_n_gl[3+4*k] = een_rescaled_single_n_gl[3+4*(k+nucl_num*(cord_num+1)*nw)] - een_rescaled_n_gl[num+elec_num*(3+4*(k+nucl_num*(cord_num+1)*nw))];
+    }
+
+    for (int64_t n=0 ; n<dim_c_vector ; ++n) {
+      const int64_t l = lkpm_combined_index[n];
+      const int64_t k = lkpm_combined_index[n+dim_c_vector];
+      const int64_t m = lkpm_combined_index[n+3*dim_c_vector];
+
+      for (int64_t kk=0 ; kk<4 ; ++kk) {
+        double* dgl = &delta_een_gl[elec_num*(kk+4*nw)];
+
+        for (int64_t a=0 ; a<nucl_num ; ++a) {
+          const double cn = c_vector_full[a+n*nucl_num];
+          if (cn == 0.) continue;
+
+          const double* dpg1_m    = &delta_p_gl[elec_num*(a+nucl_num*(kk+4*(m+(cord_num+1)*(k+cord_num*nw))))];
+          const double* dpg1_ml   = &delta_p_gl[elec_num*(a+nucl_num*(kk+4*(m+l+(cord_num+1)*(k+cord_num*nw))))];
+          const double* dp_m      = &delta_p[elec_num*(a+nucl_num*(m+(cord_num+1)*(k+cord_num*nw)))];
+          const double* dp_ml     = &delta_p[elec_num*(a+nucl_num*(m+l+(cord_num+1)*(k+cord_num*nw)))];
+          const double* een_r_m   = &een_rescaled_n[elec_num*(a+nucl_num*(m+(cord_num+1)*nw))];
+          const double* een_r_ml  = &een_rescaled_n[elec_num*(a+nucl_num*(m+l+(cord_num+1)*nw))];
+          const double* een_r_gl_m  = &een_rescaled_n_gl[elec_num*(kk+4*(a+nucl_num*(m+(cord_num+1)*nw)))];
+          const double* een_r_gl_ml = &een_rescaled_n_gl[elec_num*(kk+4*(a+nucl_num*(m+l+(cord_num+1)*nw)))];
+
+          for (int64_t i=0 ; i<elec_num ; ++i) {
+            dgl[i] += cn * (dpg1_m[i] * een_r_ml[i] + dpg1_ml[i] * een_r_m[i] +
+                            dp_m[i] * een_r_gl_ml[i] + dp_ml[i] * een_r_gl_m[i]);
+          }
+          dgl[num] += ( (dtmp_c[num+elec_num*(kk+4*(a+nucl_num*(m+(cord_num+1)*(k+cord_num*nw))))] +
+                         dpg1_m[num]
+                         ) * een_rescaled_delta_n[a+nucl_num*(m+l)] +
+                        (dtmp_c[num+elec_num*(kk+4*(a+nucl_num*(m+l+(cord_num+1)*(k+cord_num*nw))))] +
+                         dpg1_ml[num]
+                         ) * een_rescaled_delta_n[a+nucl_num*m] +
+                        (tmp_c[num+elec_num*(a+nucl_num*(m+(cord_num+1)*(k+cord_num*nw)))] +
+                         dp_m[num]
+                         ) * een_rescaled_delta_n_gl[kk+4*(a+nucl_num*(m+l))] +
+                        (tmp_c[num+elec_num*(a+nucl_num*(m+l+(cord_num+1)*(k+cord_num*nw)))] +
+                         dp_ml[num]
+                         ) * een_rescaled_delta_n_gl[kk+4*(a+nucl_num*m)] 
+                        ) * cn;
+        }
+      }
+
+      for (int a=0 ; a<nucl_num ; ++a) {
+        const double cn = 2. * c_vector_full[a+n*nucl_num];
+        if (cn == 0.) continue;
+        
+        double* dgl4 = &delta_een_gl[elec_num*(3+4*nw)];
+        
+        const double* dpg1_m    = &delta_p_gl[elec_num*(a+nucl_num*(0+4*(m+(cord_num+1)*(k+cord_num*nw))))];
+        const double* dpg2_m    = &delta_p_gl[elec_num*(a+nucl_num*(1+4*(m+(cord_num+1)*(k+cord_num*nw))))];
+        const double* dpg3_m    = &delta_p_gl[elec_num*(a+nucl_num*(2+4*(m+(cord_num+1)*(k+cord_num*nw))))];
+        const double* dpg1_ml   = &delta_p_gl[elec_num*(a+nucl_num*(0+4*(m+l+(cord_num+1)*(k+cord_num*nw))))];
+        const double* dpg2_ml   = &delta_p_gl[elec_num*(a+nucl_num*(1+4*(m+l+(cord_num+1)*(k+cord_num*nw))))];
+        const double* dpg3_ml   = &delta_p_gl[elec_num*(a+nucl_num*(2+4*(m+l+(cord_num+1)*(k+cord_num*nw))))];
+        const double* een_r_gl1_m  = &een_rescaled_n_gl[elec_num*(0+4*(a+nucl_num*(m+(cord_num+1)*nw)))];
+        const double* een_r_gl2_m  = &een_rescaled_n_gl[elec_num*(1+4*(a+nucl_num*(m+(cord_num+1)*nw)))];
+        const double* een_r_gl3_m  = &een_rescaled_n_gl[elec_num*(2+4*(a+nucl_num*(m+(cord_num+1)*nw)))];
+        const double* een_r_gl1_ml = &een_rescaled_n_gl[elec_num*(0+4*(a+nucl_num*(m+l+(cord_num+1)*nw)))];
+        const double* een_r_gl2_ml = &een_rescaled_n_gl[elec_num*(1+4*(a+nucl_num*(m+l+(cord_num+1)*nw)))];
+        const double* een_r_gl3_ml = &een_rescaled_n_gl[elec_num*(2+4*(a+nucl_num*(m+l+(cord_num+1)*nw)))];
+
+        for (int64_t i=0 ; i<elec_num ; ++i) {
+          dgl4[i] += (dpg1_m[i] * een_r_gl1_ml[i] + dpg1_ml[i] * een_r_gl1_m[i] +
+                      dpg2_m[i] * een_r_gl2_ml[i] + dpg2_ml[i] * een_r_gl2_m[i] +
+                      dpg3_m[i] * een_r_gl3_ml[i] + dpg3_ml[i] * een_r_gl3_m[i] ) * cn;
+        }
+        dgl4[num] += ((dpg1_m[num] +
+                       dtmp_c[num+elec_num*(0+4*(a+nucl_num*(m+(cord_num+1)*(k+cord_num*nw))))]
+                       ) * een_rescaled_delta_n_gl[0+4*(a+nucl_num*(m+l))] +
+                      (dpg1_ml[num] +
+                       dtmp_c[num+elec_num*(0+4*(a+nucl_num*(m+l+(cord_num+1)*(k+cord_num*nw))))]
+                       ) * een_rescaled_delta_n_gl[0+4*(a+nucl_num*m)] +
+                      (dpg2_m[num] +
+                       dtmp_c[num+elec_num*(1+4*(a+nucl_num*(m+(cord_num+1)*(k+cord_num*nw))))]
+                       ) * een_rescaled_delta_n_gl[1+4*(a+nucl_num*(m+l))] +
+                      (dpg2_ml[num] +
+                       dtmp_c[num+elec_num*(1+4*(a+nucl_num*(m+l+(cord_num+1)*(k+cord_num*nw))))]
+                       ) * een_rescaled_delta_n_gl[1+4*(a+nucl_num*m)] +
+                      (dpg3_m[num] +
+                       dtmp_c[num+elec_num*(2+4*(a+nucl_num*(m+(cord_num+1)*(k+cord_num*nw))))]
+                       ) * een_rescaled_delta_n_gl[2+4*(a+nucl_num*(m+l))] +
+                      (dpg3_ml[num] +
+                       dtmp_c[num+elec_num*(2+4*(a+nucl_num*(m+l+(cord_num+1)*(k+cord_num*nw))))]
+                       ) * een_rescaled_delta_n_gl[2+4*(a+nucl_num*m)] ) * cn;
+        
+      }
+    }
+  }
+
+  return QMCKL_SUCCESS;
+}
+
      #+end_src
 
  #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
@@ -4288,7 +4316,7 @@ qmckl_exit_code qmckl_provide_jastrow_champ_delta_p_g(qmckl_context context)
      | ~een_rescaled_e_gl~        | ~double[walk_num][0:cord_num][elec_num][4][elec_num]~               | in     | Electron-electron rescaled distances derivatives        |
      | ~een_rescaled_single_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4]~                         | in     | Electron-nucleus single rescaled distances derivatives  |
      | ~een_rescaled_single_e_gl~ | ~double[walk_num][0:cord_num][elec_num][4]~                         | in     | Electron-electron single rescaled distances derivatives |
-     | ~delta_p_g~               | ~double[walk_num][0:cord_num-1][0:cord_num][4][nucl_num][elec_num]~ | out    | Delta P matrix gradient and Laplacian                   |
+     | ~delta_p_g~                | ~double[walk_num][0:cord_num-1][0:cord_num][4][nucl_num][elec_num]~ | out    | Delta P matrix gradient and Laplacian                   |
      |----------------------------+---------------------------------------------------------------------+--------+---------------------------------------------------------|
 
      #+begin_src f90 :comments org :tangle (eval f) :noweb yes