diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org
index 533a9f9..abeac12 100644
--- a/org/qmckl_jastrow.org
+++ b/org/qmckl_jastrow.org
@@ -58,7 +58,7 @@
   The terms $J_{\text{ee}}^\infty$ and $J_{\text{eN}}^\infty$ are shifts to ensure that 
   $J_{\text{ee}}$ and $J_{\text{eN}}$ have an asymptotic value of zero.
   
-
+  
 * Headers                                                          :noexport:
   #+begin_src elisp :noexport :results none
 (org-babel-lob-ingest "../tools/lib.org")
@@ -412,10 +412,10 @@ qmckl_exit_code  qmckl_get_jastrow_aord_num          (qmckl_context context, int
 qmckl_exit_code  qmckl_get_jastrow_bord_num          (qmckl_context context, int64_t* const bord_num);
 qmckl_exit_code  qmckl_get_jastrow_cord_num          (qmckl_context context, int64_t* const bord_num);
 qmckl_exit_code  qmckl_get_jastrow_type_nucl_num     (qmckl_context context, int64_t* const type_nucl_num);
-qmckl_exit_code  qmckl_get_jastrow_type_nucl_vector  (qmckl_context context, int64_t* const type_nucl_num);
-qmckl_exit_code  qmckl_get_jastrow_aord_vector       (qmckl_context context, double * const aord_vector);
-qmckl_exit_code  qmckl_get_jastrow_bord_vector       (qmckl_context context, double * const bord_vector);
-qmckl_exit_code  qmckl_get_jastrow_cord_vector       (qmckl_context context, double * const cord_vector);
+qmckl_exit_code  qmckl_get_jastrow_type_nucl_vector  (qmckl_context context, int64_t* const type_nucl_num, int64_t* size_max);
+qmckl_exit_code  qmckl_get_jastrow_aord_vector       (qmckl_context context, double * const aord_vector, int64_t* size_max);
+qmckl_exit_code  qmckl_get_jastrow_bord_vector       (qmckl_context context, double * const bord_vector, int64_t* size_max);
+qmckl_exit_code  qmckl_get_jastrow_cord_vector       (qmckl_context context, double * const cord_vector, int64_t* size_max);
    #+end_src
 
    Along with these core functions, calculation of the jastrow factor
@@ -559,7 +559,7 @@ qmckl_exit_code qmckl_get_jastrow_type_nucl_num (const qmckl_context context, in
   return QMCKL_SUCCESS;
 }
 
-qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (const qmckl_context context, int64_t * const type_nucl_vector) {
+qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (const qmckl_context context, int64_t * const type_nucl_vector, int64_t* size_max) {
 
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return (char) 0;
@@ -583,10 +583,11 @@ qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (const qmckl_context context,
 
   assert (ctx->jastrow.type_nucl_vector != NULL);
   memcpy(type_nucl_vector, ctx->jastrow.type_nucl_vector, ctx->jastrow.type_nucl_num*sizeof(int64_t));
+  (*size_max) = ctx->jastrow.type_nucl_num;
   return QMCKL_SUCCESS;
 }
 
-qmckl_exit_code qmckl_get_jastrow_aord_vector (const qmckl_context context, double * const aord_vector) {
+qmckl_exit_code qmckl_get_jastrow_aord_vector (const qmckl_context context, double * const aord_vector, int64_t* size_max) {
 
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return (char) 0;
@@ -609,11 +610,13 @@ qmckl_exit_code qmckl_get_jastrow_aord_vector (const qmckl_context context, doub
   }
 
   assert (ctx->jastrow.aord_vector != NULL);
-  memcpy(aord_vector, ctx->jastrow.aord_vector, ctx->jastrow.aord_num*sizeof(double));
+  int64_t sze = (ctx->jastrow.aord_num + 1)*ctx->jastrow.type_nucl_num;
+  memcpy(aord_vector, ctx->jastrow.aord_vector, sze*sizeof(double));
+  (*size_max) = sze;
   return QMCKL_SUCCESS;
 }
 
-qmckl_exit_code qmckl_get_jastrow_bord_vector (const qmckl_context context, double * const bord_vector) {
+qmckl_exit_code qmckl_get_jastrow_bord_vector (const qmckl_context context, double * const bord_vector, int64_t* size_max) {
 
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return (char) 0;
@@ -636,11 +639,13 @@ qmckl_exit_code qmckl_get_jastrow_bord_vector (const qmckl_context context, doub
   }
 
   assert (ctx->jastrow.bord_vector != NULL);
-  memcpy(bord_vector, ctx->jastrow.bord_vector, ctx->jastrow.bord_num*sizeof(double));
+  int64_t sze=ctx->jastrow.bord_num +1;
+  memcpy(bord_vector, ctx->jastrow.bord_vector, sze*sizeof(double));
+  (*size_max) = sze;
   return QMCKL_SUCCESS;
 }
 
-qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, double * const cord_vector) {
+qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, double * const cord_vector, int64_t* size_max) {
 
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return (char) 0;
@@ -663,7 +668,14 @@ qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, doub
   }
 
   assert (ctx->jastrow.cord_vector != NULL);
-  memcpy(cord_vector, ctx->jastrow.cord_vector, ctx->jastrow.dim_cord_vect*sizeof(double));
+  
+  int64_t dim_cord_vect;
+  qmckl_exit_code rc = qmckl_get_jastrow_dim_cord_vect(context, &dim_cord_vect);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  int64_t sze=dim_cord_vect * ctx->jastrow.type_nucl_num;
+  memcpy(cord_vector, ctx->jastrow.cord_vector, sze*sizeof(double));
+  (*size_max) = sze;
   return QMCKL_SUCCESS;
 }
 
@@ -678,9 +690,9 @@ qmckl_exit_code qmckl_get_jastrow_cord_vector (const qmckl_context context, doub
 qmckl_exit_code  qmckl_set_jastrow_ord_num           (qmckl_context context, const int64_t aord_num, const int64_t bord_num, const int64_t cord_num);
 qmckl_exit_code  qmckl_set_jastrow_type_nucl_num     (qmckl_context context, const int64_t type_nucl_num);
 qmckl_exit_code  qmckl_set_jastrow_type_nucl_vector  (qmckl_context context, const int64_t* type_nucl_vector, const int64_t nucl_num);
-qmckl_exit_code  qmckl_set_jastrow_aord_vector       (qmckl_context context, const double * aord_vector);
-qmckl_exit_code  qmckl_set_jastrow_bord_vector       (qmckl_context context, const double * bord_vector);
-qmckl_exit_code  qmckl_set_jastrow_cord_vector       (qmckl_context context, const double * cord_vector);
+qmckl_exit_code  qmckl_set_jastrow_aord_vector       (qmckl_context context, const double * aord_vector, int64_t size_max);
+qmckl_exit_code  qmckl_set_jastrow_bord_vector       (qmckl_context context, const double * bord_vector, int64_t size_max);
+qmckl_exit_code  qmckl_set_jastrow_cord_vector       (qmckl_context context, const double * cord_vector, int64_t size_max);
    #+end_src
 
    #+NAME:pre2
@@ -804,7 +816,7 @@ qmckl_exit_code qmckl_set_jastrow_type_nucl_vector(qmckl_context context, int64_
   <<post2>>
 }
 
-qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double const * aord_vector) {
+qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double const * aord_vector, int64_t size_max) {
 <<pre2>>
 
   int32_t mask = 1 << 3;
@@ -837,11 +849,19 @@ qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double cons
       return qmckl_failwith( context, rc,
                              "qmckl_set_ord_vector",
                              NULL);
-    }
+}
   }
 
   qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
   mem_info.size = (aord_num + 1) * type_nucl_num * sizeof(double);
+
+  if (size_max < mem_info.size/sizeof(double)) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_set_jastrow_aord_vector",
+                           "Array too small. Expected (aord_num+1)*type_nucl_num");
+  }
+
   double* new_array = (double*) qmckl_malloc(context, mem_info);
 
   if(new_array == NULL) {
@@ -858,7 +878,7 @@ qmckl_exit_code qmckl_set_jastrow_aord_vector(qmckl_context context, double cons
   <<post2>>
 }
 
-qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double const * bord_vector) {
+qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double const * bord_vector, int64_t size_max) {
 <<pre2>>
 
   int32_t mask = 1 << 4;
@@ -892,6 +912,14 @@ qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double cons
 
   qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
   mem_info.size = (bord_num + 1) * sizeof(double);
+
+  if (size_max < mem_info.size/sizeof(double)) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_set_jastrow_bord_vector",
+                           "Array too small. Expected (bord_num+1)");
+  }
+
   double* new_array = (double*) qmckl_malloc(context, mem_info);
 
   if(new_array == NULL) {
@@ -908,7 +936,7 @@ qmckl_exit_code qmckl_set_jastrow_bord_vector(qmckl_context context, double cons
   <<post2>>
 }
 
-qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double const * cord_vector) {
+qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double const * cord_vector, int64_t size_max) {
 <<pre2>>
 
   int32_t mask = 1 << 5;
@@ -949,6 +977,14 @@ qmckl_exit_code qmckl_set_jastrow_cord_vector(qmckl_context context, double cons
 
   qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
   mem_info.size = dim_cord_vect * type_nucl_num * sizeof(double);
+
+  if (size_max < mem_info.size/sizeof(double)) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_set_jastrow_cord_vector",
+                           "Array too small. Expected dim_cord_vect * type_nucl_num");
+  }
+
   double* new_array = (double*) qmckl_malloc(context, mem_info);
 
   if(new_array == NULL) {
@@ -1210,11 +1246,11 @@ assert(qmckl_nucleus_provided(context));
 
 *** Get
     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
-qmckl_exit_code qmckl_get_jastrow_asymp_jasb(qmckl_context context, double* const asymp_jasb);
+qmckl_exit_code qmckl_get_jastrow_asymp_jasb(qmckl_context context, double* const asymp_jasb, int64_t* size_max);
     #+end_src
 
     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_get_jastrow_asymp_jasb(qmckl_context context, double* const asymp_jasb)
+qmckl_exit_code qmckl_get_jastrow_asymp_jasb(qmckl_context context, double* const asymp_jasb, int64_t* size_max)
 {
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return QMCKL_NULL_CONTEXT;
@@ -1230,6 +1266,7 @@ qmckl_exit_code qmckl_get_jastrow_asymp_jasb(qmckl_context context, double* cons
 
   size_t sze = 2;
   memcpy(asymp_jasb, ctx->jastrow.asymp_jasb, sze * sizeof(double));
+  (*size_max) = sze;
 
   return QMCKL_SUCCESS;
 }
@@ -1420,7 +1457,7 @@ print("asym_one         : ", asym_one)
 print("asymp_jasb[0]    : ", asymp_jasb[0])
 print("asymp_jasb[1]    : ", asymp_jasb[1])
     #+end_src
-
+    
     #+RESULTS: asymp_jasb
     : asym_one         :  0.43340325572525706
     : asymp_jasb[0]    :  0.5323750557252571
@@ -1437,6 +1474,7 @@ int64_t cord_num = n2_cord_num;
 double* aord_vector = &(n2_aord_vector[0][0]);
 double* bord_vector = &(n2_bord_vector[0]);
 double* cord_vector = &(n2_cord_vector[0][0]);
+int64_t dim_cord_vect=0;
 
 /* Initialize the Jastrow data */
 rc = qmckl_init_jastrow(context);
@@ -1449,20 +1487,21 @@ rc = qmckl_set_jastrow_type_nucl_num(context, type_nucl_num);
 assert(rc == QMCKL_SUCCESS);
 rc = qmckl_set_jastrow_type_nucl_vector(context, type_nucl_vector, nucl_num);
 assert(rc == QMCKL_SUCCESS);
-rc = qmckl_set_jastrow_aord_vector(context, aord_vector);
+rc = qmckl_set_jastrow_aord_vector(context, aord_vector,(aord_num+1)*type_nucl_num);
 assert(rc == QMCKL_SUCCESS);
-rc = qmckl_set_jastrow_bord_vector(context, bord_vector);
+rc = qmckl_set_jastrow_bord_vector(context, bord_vector,(bord_num+1));
 assert(rc == QMCKL_SUCCESS);
-rc = qmckl_set_jastrow_bord_vector(context, bord_vector);
+rc = qmckl_get_jastrow_dim_cord_vect(context, &dim_cord_vect);
 assert(rc == QMCKL_SUCCESS);
-rc = qmckl_set_jastrow_cord_vector(context, cord_vector);
+rc = qmckl_set_jastrow_cord_vector(context, cord_vector,dim_cord_vect*type_nucl_num);
 assert(rc == QMCKL_SUCCESS);
 
 /* Check if Jastrow is properly initialized */
 assert(qmckl_jastrow_provided(context));
 
 double asymp_jasb[2];
-rc = qmckl_get_jastrow_asymp_jasb(context, asymp_jasb);
+int64_t size_max=0;
+rc = qmckl_get_jastrow_asymp_jasb(context, asymp_jasb,&size_max);
 
 // calculate asymp_jasb
 assert(fabs(asymp_jasb[0]-0.5323750557252571) < 1.e-12);
@@ -1482,11 +1521,11 @@ f_{ee} = \sum_{i,j<i} \left\{ \frac{ \eta B_0 C_{ij}}{1 - B_1 C_{ij}}  - J_{asym
 
 *** Get
     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
-qmckl_exit_code qmckl_get_jastrow_factor_ee(qmckl_context context, double* const factor_ee);
+qmckl_exit_code qmckl_get_jastrow_factor_ee(qmckl_context context, double* const factor_ee, int64_t* size_max);
     #+end_src
 
     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_get_jastrow_factor_ee(qmckl_context context, double* const factor_ee)
+qmckl_exit_code qmckl_get_jastrow_factor_ee(qmckl_context context, double* const factor_ee, int64_t* size_max)
 {
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return QMCKL_NULL_CONTEXT;
@@ -1500,7 +1539,9 @@ qmckl_exit_code qmckl_get_jastrow_factor_ee(qmckl_context context, double* const
   qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
   assert (ctx != NULL);
 
-  memcpy(factor_ee, ctx->jastrow.factor_ee, ctx->electron.walk_num*sizeof(double));
+  int64_t sze=ctx->electron.walk_num;
+  memcpy(factor_ee, ctx->jastrow.factor_ee, sze*sizeof(double));
+  (*size_max) = sze;
 
   return QMCKL_SUCCESS;
 }
@@ -1597,7 +1638,7 @@ integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num,
   integer(qmckl_context), intent(in)  :: context
   integer*8             , intent(in)  :: walk_num, elec_num, bord_num, up_num
   double precision      , intent(in)  :: bord_vector(bord_num + 1)
-  double precision      , intent(in)  :: ee_distance_rescaled(walk_num, elec_num, elec_num)
+  double precision      , intent(in)  :: ee_distance_rescaled(elec_num, elec_num, walk_num)
   double precision      , intent(in)  :: asymp_jasb(2)
   double precision      , intent(out) :: factor_ee(walk_num)
 
@@ -1631,13 +1672,13 @@ integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num,
   do nw =1, walk_num
   do j = 1, elec_num
      do i = 1, j - 1
-        x = ee_distance_rescaled(nw,i,j)
+        x = ee_distance_rescaled(i,j,nw)
         power_ser = 0.0d0
         spin_fact = 1.0d0
         ipar = 1
 
         do p = 2, bord_num
-          x = x * ee_distance_rescaled(nw,i,j)
+          x = x * ee_distance_rescaled(i,j,nw)
           power_ser = power_ser + bord_vector(p + 1) * x
         end do
 
@@ -1647,9 +1688,9 @@ integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num,
         endif
 
         factor_ee(nw) = factor_ee(nw) + spin_fact * bord_vector(1)  * &
-                                ee_distance_rescaled(nw,i,j) / &
+                                ee_distance_rescaled(i,j,nw) / &
                                 (1.0d0 + bord_vector(2) *   &
-                                ee_distance_rescaled(nw,i,j))  &
+                                ee_distance_rescaled(i,j,nw))  &
                                -asymp_jasb(ipar) + power_ser
 
      end do
@@ -1763,7 +1804,8 @@ print("factor_ee        :",factor_ee)
 assert(qmckl_jastrow_provided(context));
 
 double factor_ee[walk_num];
-rc = qmckl_get_jastrow_factor_ee(context, factor_ee);
+size_max=0;
+rc = qmckl_get_jastrow_factor_ee(context, factor_ee, &size_max);
 
 // calculate factor_ee
 assert(fabs(factor_ee[0]+4.282760865958113) < 1.e-12);
@@ -1782,11 +1824,11 @@ assert(fabs(factor_ee[0]+4.282760865958113) < 1.e-12);
 
 *** Get
     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
-qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, double* const factor_ee_deriv_e);
+qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, double* const factor_ee_deriv_e, int64_t* size_max);
     #+end_src
 
     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, double* const factor_ee_deriv_e)
+qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, double* const factor_ee_deriv_e, int64_t* size_max)
 {
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return QMCKL_NULL_CONTEXT;
@@ -1802,6 +1844,7 @@ qmckl_exit_code qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context, doubl
 
   int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num;
   memcpy(factor_ee_deriv_e, ctx->jastrow.factor_ee_deriv_e, sze * sizeof(double));
+  (*size_max) = sze;
 
   return QMCKL_SUCCESS;
 }
@@ -1905,8 +1948,8 @@ integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num,
   integer(qmckl_context), intent(in)  :: context
   integer*8             , intent(in)  :: walk_num, elec_num, bord_num, up_num
   double precision      , intent(in)  :: bord_vector(bord_num + 1)
-  double precision      , intent(in)  :: ee_distance_rescaled(walk_num, elec_num, elec_num)
-  double precision      , intent(in)  :: ee_distance_rescaled_deriv_e(walk_num, 4, elec_num, elec_num)
+  double precision      , intent(in)  :: ee_distance_rescaled(elec_num, elec_num,walk_num)
+  double precision      , intent(in)  :: ee_distance_rescaled_deriv_e(4,elec_num, elec_num,walk_num)
   double precision      , intent(in)  :: asymp_jasb(2)
   double precision      , intent(out) :: factor_ee_deriv_e(elec_num,4,walk_num)
 
@@ -1945,7 +1988,7 @@ integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num,
   do nw =1, walk_num
   do j = 1, elec_num
      do i = 1, elec_num
-        x = ee_distance_rescaled(nw, i, j)
+        x = ee_distance_rescaled(i,j,nw)
         if(abs(x) < 1.0d-18) cycle
         pow_ser_g   = 0.0d0
         spin_fact   = 1.0d0
@@ -1956,9 +1999,10 @@ integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num,
         xinv        = 1.0d0 / (x + 1.0d-18)
         ipar        = 1
 
-        do ii = 1, 4
-          dx(ii) = ee_distance_rescaled_deriv_e(nw, ii, i, j)
-        end do
+        dx(1) = ee_distance_rescaled_deriv_e(1, i, j, nw)
+        dx(2) = ee_distance_rescaled_deriv_e(2, i, j, nw)
+        dx(3) = ee_distance_rescaled_deriv_e(3, i, j, nw)
+        dx(4) = ee_distance_rescaled_deriv_e(4, i, j, nw)  
 
         if((i .LE. up_num .AND. j .LE. up_num ) .OR.  &
            (i .GT. up_num .AND. j .GT. up_num)) then
@@ -1969,14 +2013,14 @@ integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num,
         lap2 = 0.0d0
         lap3 = 0.0d0
         do ii = 1, 3
-          x = ee_distance_rescaled(nw, i, j)
+          x = ee_distance_rescaled(i, j, nw)
           if(abs(x) < 1.0d-18) cycle
           do p = 2, bord_num
             y = p * bord_vector(p + 1) * x
             pow_ser_g(ii) = pow_ser_g(ii) + y * dx(ii)
             lap1 = lap1 + (p - 1) * y * xinv * dx(ii) * dx(ii)
             lap2 = lap2 + y
-            x = x * ee_distance_rescaled(nw, i, j)
+            x = x * ee_distance_rescaled(i, j, nw)
           end do
 
           lap3 = lap3 - 2.0d0 * bord_vector(2) * dx(ii) * dx(ii)
@@ -2173,7 +2217,8 @@ assert(qmckl_jastrow_provided(context));
 
 // calculate factor_ee_deriv_e
 double factor_ee_deriv_e[walk_num][4][elec_num];
-rc = qmckl_get_jastrow_factor_ee_deriv_e(context, &(factor_ee_deriv_e[0][0][0]));
+size_max=0;
+rc = qmckl_get_jastrow_factor_ee_deriv_e(context, &(factor_ee_deriv_e[0][0][0]),&size_max);
 
 // check factor_ee_deriv_e
 assert(fabs(factor_ee_deriv_e[0][0][0]-0.16364894652107934) < 1.e-12);
@@ -2195,11 +2240,11 @@ f_{en} = \sum_{i,j<i} \left\{ \frac{ A_0 C_{ij}}{1 - A_1 C_{ij}} + \sum^{nord}_{
 
 *** Get
     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
-qmckl_exit_code qmckl_get_jastrow_factor_en(qmckl_context context, double* const factor_en);
+qmckl_exit_code qmckl_get_jastrow_factor_en(qmckl_context context, double* const factor_en, int64_t* size_max);
     #+end_src
 
     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_get_jastrow_factor_en(qmckl_context context, double* const factor_en)
+qmckl_exit_code qmckl_get_jastrow_factor_en(qmckl_context context, double* const factor_en, int64_t* size_max)
 {
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return QMCKL_NULL_CONTEXT;
@@ -2213,7 +2258,9 @@ qmckl_exit_code qmckl_get_jastrow_factor_en(qmckl_context context, double* const
   qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
   assert (ctx != NULL);
 
-  memcpy(factor_en, ctx->jastrow.factor_en, ctx->electron.walk_num*sizeof(double));
+  int64_t sze=ctx->electron.walk_num;
+  memcpy(factor_en, ctx->jastrow.factor_en, sze*sizeof(double));
+  (*size_max)=sze;
 
   return QMCKL_SUCCESS;
 }
@@ -2314,7 +2361,7 @@ integer function qmckl_compute_factor_en_f(context, walk_num, elec_num, nucl_num
   integer*8             , intent(in)  :: walk_num, elec_num, aord_num, nucl_num, type_nucl_num
   integer*8             , intent(in)  :: type_nucl_vector(nucl_num)
   double precision      , intent(in)  :: aord_vector(aord_num + 1, type_nucl_num)
-  double precision      , intent(in)  :: en_distance_rescaled(walk_num, nucl_num, elec_num)
+  double precision      , intent(in)  :: en_distance_rescaled(elec_num, nucl_num, walk_num)
   double precision      , intent(out) :: factor_en(walk_num)
 
   integer*8 :: i, a, p, ipar, nw
@@ -2352,18 +2399,18 @@ integer function qmckl_compute_factor_en_f(context, walk_num, elec_num, nucl_num
   do nw =1, walk_num
   do a = 1, nucl_num
      do i = 1, elec_num
-        x = en_distance_rescaled(nw, a, i)
+        x = en_distance_rescaled(i, a, nw)
         power_ser = 0.0d0
 
         do p = 2, aord_num
-          x = x * en_distance_rescaled(nw, a, i)
+          x = x * en_distance_rescaled(i, a, nw)
           power_ser = power_ser + aord_vector(p + 1, type_nucl_vector(a)) * x
         end do
 
         factor_en(nw) = factor_en(nw) + aord_vector(1, type_nucl_vector(a)) *  &
-                                en_distance_rescaled(nw, a, i) /               &
+                                en_distance_rescaled(i, a, nw) /               &
                                 (1.0d0 + aord_vector(2, type_nucl_vector(a)) * &
-                                en_distance_rescaled(nw, a, i))                &
+                                en_distance_rescaled(i, a, nw))                &
                                 + power_ser
 
      end do
@@ -2419,7 +2466,7 @@ end function qmckl_compute_factor_en_f
       integer (c_int64_t) , intent(in)          :: type_nucl_vector(nucl_num)
       integer (c_int64_t) , intent(in)  , value :: aord_num
       real    (c_double ) , intent(in)          :: aord_vector(aord_num + 1, type_nucl_num)
-      real    (c_double ) , intent(in)          :: en_distance_rescaled(walk_num, nucl_num, elec_num)
+      real    (c_double ) , intent(in)          :: en_distance_rescaled(elec_num, nucl_num, walk_num)
       real    (c_double ) , intent(out)         :: factor_en(walk_num)
 
       integer(c_int32_t), external :: qmckl_compute_factor_en_f
@@ -2470,7 +2517,8 @@ print("factor_en        :",factor_en)
 assert(qmckl_jastrow_provided(context));
 
 double factor_en[walk_num];
-rc = qmckl_get_jastrow_factor_en(context, factor_en);
+size_max=0;
+rc = qmckl_get_jastrow_factor_en(context, factor_en,&size_max);
 
 // calculate factor_en
 assert(fabs(factor_en[0]+5.865822569188727) < 1.e-12);
@@ -2486,11 +2534,11 @@ assert(fabs(factor_en[0]+5.865822569188727) < 1.e-12);
 
 *** Get
     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
-qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, double* const factor_en_deriv_e);
+qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, double* const factor_en_deriv_e, int64_t* size_max);
     #+end_src
 
     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, double* const factor_en_deriv_e)
+qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, double* const factor_en_deriv_e, int64_t* size_max)
 {
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return QMCKL_NULL_CONTEXT;
@@ -2506,6 +2554,7 @@ qmckl_exit_code qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context, doubl
 
   int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num;
   memcpy(factor_en_deriv_e, ctx->jastrow.factor_en_deriv_e, sze*sizeof(double));
+  (*size_max) = sze;
 
   return QMCKL_SUCCESS;
 }
@@ -2612,8 +2661,8 @@ integer function qmckl_compute_factor_en_deriv_e_f(context, walk_num, elec_num,
   integer*8             , intent(in)  :: walk_num, elec_num, aord_num, nucl_num, type_nucl_num
   integer*8             , intent(in)  :: type_nucl_vector(nucl_num)
   double precision      , intent(in)  :: aord_vector(aord_num + 1, type_nucl_num)
-  double precision      , intent(in)  :: en_distance_rescaled(walk_num, elec_num, nucl_num)
-  double precision      , intent(in)  :: en_distance_rescaled_deriv_e(walk_num, 4, elec_num, nucl_num)
+  double precision      , intent(in)  :: en_distance_rescaled(elec_num, nucl_num, walk_num)
+  double precision      , intent(in)  :: en_distance_rescaled_deriv_e(4, elec_num, nucl_num, walk_num)
   double precision      , intent(out) :: factor_en_deriv_e(elec_num,4,walk_num)
 
   integer*8 :: i, a, p, ipar, nw, ii
@@ -2655,7 +2704,7 @@ integer function qmckl_compute_factor_en_deriv_e_f(context, walk_num, elec_num,
   do nw =1, walk_num
   do a = 1, nucl_num
      do i = 1, elec_num
-        x = en_distance_rescaled(nw, i, a)
+        x = en_distance_rescaled(i,a,nw)
         if(abs(x) < 1.0d-18) continue
         power_ser_g = 0.0d0
         den = 1.0d0 + aord_vector(2, type_nucl_vector(a)) * x
@@ -2665,20 +2714,20 @@ integer function qmckl_compute_factor_en_deriv_e_f(context, walk_num, elec_num,
         xinv = 1.0d0 / x
 
         do ii = 1, 4
-          dx(ii) = en_distance_rescaled_deriv_e(nw, ii, i, a)
+          dx(ii) = en_distance_rescaled_deriv_e(ii,i,a,nw)
         end do
 
         lap1 = 0.0d0
         lap2 = 0.0d0
         lap3 = 0.0d0
         do ii = 1, 3
-          x = en_distance_rescaled(nw, i, a)
+          x = en_distance_rescaled(i, a, nw)
           do p = 2, aord_num
             y = p * aord_vector(p + 1, type_nucl_vector(a)) * x
             power_ser_g(ii) = power_ser_g(ii) + y * dx(ii)
             lap1 = lap1 + (p - 1) * y * xinv * dx(ii) * dx(ii)
             lap2 = lap2 + y
-            x = x * en_distance_rescaled(nw, i, a)
+            x = x * en_distance_rescaled(i, a, nw)
           end do
 
           lap3 = lap3 - 2.0d0 * aord_vector(2, type_nucl_vector(a)) * dx(ii) * dx(ii)
@@ -2872,7 +2921,8 @@ assert(qmckl_jastrow_provided(context));
 
 // calculate factor_en_deriv_e
 double factor_en_deriv_e[walk_num][4][elec_num];
-rc = qmckl_get_jastrow_factor_en_deriv_e(context, &(factor_en_deriv_e[0][0][0]));
+size_max=0;
+rc = qmckl_get_jastrow_factor_en_deriv_e(context, &(factor_en_deriv_e[0][0][0]),&size_max);
 
 // check factor_en_deriv_e
 assert(fabs(factor_en_deriv_e[0][0][0]-0.11609919541763383) < 1.e-12);
@@ -2896,11 +2946,11 @@ assert(fabs(factor_en_deriv_e[0][3][0]+0.9667363412285741 ) < 1.e-12);
 *** Get
 
     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
-qmckl_exit_code qmckl_get_jastrow_een_rescaled_e(qmckl_context context, double* const distance_rescaled);
+qmckl_exit_code qmckl_get_jastrow_een_rescaled_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max);
     #+end_src
 
     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_get_jastrow_een_rescaled_e(qmckl_context context, double* const distance_rescaled)
+qmckl_exit_code qmckl_get_jastrow_een_rescaled_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max)
 {
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return QMCKL_NULL_CONTEXT;
@@ -2916,6 +2966,7 @@ qmckl_exit_code qmckl_get_jastrow_een_rescaled_e(qmckl_context context, double*
 
   size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
   memcpy(distance_rescaled, ctx->jastrow.een_rescaled_e, sze * sizeof(double));
+  (*size_max) = sze;
 
   return QMCKL_SUCCESS;
 }
@@ -2997,7 +3048,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context)
     | ~cord_num~                | ~int64_t~                                          | in     | Order of polynomials                 |
     | ~rescale_factor_kappa_ee~ | ~double~                                           | in     | Factor to rescale ee distances       |
     | ~ee_distance~             | ~double[walk_num][elec_num][elec_num]~             | in     | Electron-electron distances          |
-    | ~een_rescaled_e~          | ~double[walk_num][elec_num][elec_num][0:cord_num]~ | out    | Electron-electron rescaled distances |
+    | ~een_rescaled_e~          | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | out    | Electron-electron rescaled distances |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
 integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee,  &
@@ -3011,7 +3062,7 @@ integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cor
   integer*8             , intent(in)  :: cord_num
   double precision      , intent(in)  :: rescale_factor_kappa_ee
   double precision      , intent(in)  :: ee_distance(elec_num,elec_num,walk_num)
-  double precision      , intent(out) :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num)
+  double precision      , intent(out) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
   double precision,dimension(:,:),allocatable :: een_rescaled_e_ij
   double precision                    :: x
   integer*8                           :: i, j, k, l, nw
@@ -3062,22 +3113,22 @@ integer function qmckl_compute_een_rescaled_e_f(context, walk_num, elec_num, cor
   end do
 
   ! prepare the actual een table
-  een_rescaled_e(0, :, :, nw) = 1.0d0
+  een_rescaled_e(:, :, 0, nw) = 1.0d0
   do l = 1, cord_num
     k = 0
     do j = 1, elec_num
       do i = 1, j - 1
         k = k + 1
         x = een_rescaled_e_ij(k, l + 1)
-        een_rescaled_e(l, i, j, nw) = x
-        een_rescaled_e(l, j, i, nw) = x
+        een_rescaled_e(i, j, l, nw) = x
+        een_rescaled_e(j, i, l, nw) = x
       end do
     end do
   end do
 
   do l = 0, cord_num
     do j = 1, elec_num
-      een_rescaled_e(l, j, j, nw) = 0.0d0
+      een_rescaled_e(j, j, l, nw) = 0.0d0
     end do
   end do
 
@@ -3117,7 +3168,7 @@ end function qmckl_compute_een_rescaled_e_f
       integer (c_int64_t) , intent(in)  , value :: cord_num
       real    (c_double ) , intent(in)  , value :: rescale_factor_kappa_ee
       real    (c_double ) , intent(in)          :: ee_distance(elec_num,elec_num,walk_num)
-      real    (c_double ) , intent(out)         :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num)
+      real    (c_double ) , intent(out)         :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
 
       integer(c_int32_t), external :: qmckl_compute_een_rescaled_e_f
       info = qmckl_compute_een_rescaled_e_f &
@@ -3127,7 +3178,7 @@ end function qmckl_compute_een_rescaled_e_f
     #+end_src
 
 *** Test
-
+    
     #+begin_src python :results output :exports none :noweb yes
 import numpy as np
 
@@ -3190,16 +3241,17 @@ print(" een_rescaled_e[1, 5, 2] = ",een_rescaled_e[1, 5, 2])
 assert(qmckl_electron_provided(context));
 
 
-double een_rescaled_e[walk_num][elec_num][elec_num][(cord_num + 1)];
-rc = qmckl_get_jastrow_een_rescaled_e(context, &(een_rescaled_e[0][0][0][0]));
+double een_rescaled_e[walk_num][(cord_num + 1)][elec_num][elec_num];
+size_max=0;
+rc = qmckl_get_jastrow_een_rescaled_e(context, &(een_rescaled_e[0][0][0][0]),&size_max);
 
 // value of (0,2,1)
-assert(fabs(een_rescaled_e[0][0][2][1]-0.08084493981483197)   < 1.e-12);
-assert(fabs(een_rescaled_e[0][0][3][1]-0.1066745707571846)    < 1.e-12);
-assert(fabs(een_rescaled_e[0][0][4][1]-0.01754273169464735)   < 1.e-12);
-assert(fabs(een_rescaled_e[0][1][3][2]-0.02214680362033448)   < 1.e-12);
-assert(fabs(een_rescaled_e[0][1][4][2]-0.0005700154999202759) < 1.e-12);
-assert(fabs(een_rescaled_e[0][1][5][2]-0.3424402276009091)    < 1.e-12);
+assert(fabs(een_rescaled_e[0][1][0][2]-0.08084493981483197)   < 1.e-12);
+assert(fabs(een_rescaled_e[0][1][0][3]-0.1066745707571846)    < 1.e-12);
+assert(fabs(een_rescaled_e[0][1][0][4]-0.01754273169464735)   < 1.e-12);
+assert(fabs(een_rescaled_e[0][2][1][3]-0.02214680362033448)   < 1.e-12);
+assert(fabs(een_rescaled_e[0][2][1][4]-0.0005700154999202759) < 1.e-12);
+assert(fabs(een_rescaled_e[0][2][1][5]-0.3424402276009091)    < 1.e-12);
 
      #+end_src
 
@@ -3216,11 +3268,11 @@ assert(fabs(een_rescaled_e[0][1][5][2]-0.3424402276009091)    < 1.e-12);
 *** Get
 
     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
-qmckl_exit_code qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context, double* const distance_rescaled);
+qmckl_exit_code qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max);
     #+end_src
 
     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context, double* const distance_rescaled)
+qmckl_exit_code qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max)
 {
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return QMCKL_NULL_CONTEXT;
@@ -3236,6 +3288,7 @@ qmckl_exit_code qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context,
 
   size_t sze = ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
   memcpy(distance_rescaled, ctx->jastrow.een_rescaled_e_deriv_e, sze * sizeof(double));
+  (*size_max) = sze;
 
   return QMCKL_SUCCESS;
 }
@@ -3320,8 +3373,8 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context)
     | ~rescale_factor_kappa_ee~ | ~double~                                              | in     | Factor to rescale ee distances       |
     | ~coord_new~               | ~double[walk_num][3][elec_num]~                       | in     | Electron coordinates                 |
     | ~ee_distance~             | ~double[walk_num][elec_num][elec_num]~                | in     | Electron-electron distances          |
-    | ~een_rescaled_e~          | ~double[walk_num][elec_num][elec_num][0:cord_num]~    | in     | Electron-electron distances          |
-    | ~een_rescaled_e_deriv_e~  | ~double[walk_num][elec_num][4][elec_num][0:cord_num]~ | out    | Electron-electron rescaled distances |
+    | ~een_rescaled_e~          | ~double[walk_num][0:cord_num][elec_num][elec_num]~    | in     | Electron-electron distances          |
+    | ~een_rescaled_e_deriv_e~  | ~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
 integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee,  &
@@ -3336,8 +3389,8 @@ integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f(context, walk_num
   double precision      , intent(in)  :: rescale_factor_kappa_ee
   double precision      , intent(in)  :: coord_new(elec_num,3,walk_num)
   double precision      , intent(in)  :: ee_distance(elec_num,elec_num,walk_num)
-  double precision      , intent(in)  :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num)
-  double precision      , intent(out) :: een_rescaled_e_deriv_e(0:cord_num,elec_num,4,elec_num,walk_num)
+  double precision      , intent(in)  :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
+  double precision      , intent(out) :: een_rescaled_e_deriv_e(elec_num,4,elec_num,0:cord_num,walk_num)
   double precision,dimension(:,:,:),allocatable  :: elec_dist_deriv_e
   double precision                    :: x, rij_inv, kappa_l
   integer*8                           :: i, j, k, l, nw, ii
@@ -3385,19 +3438,24 @@ integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f(context, walk_num
       kappa_l = - dble(l) * rescale_factor_kappa_ee
       do j = 1, elec_num
         do i = 1, elec_num
-          do ii = 1, 4
-            een_rescaled_e_deriv_e(l, i, ii, j, nw) = kappa_l * elec_dist_deriv_e(ii, i, j)
-          end do
+          een_rescaled_e_deriv_e(i, 1, j, l, nw) = kappa_l * elec_dist_deriv_e(1, i, j)
+          een_rescaled_e_deriv_e(i, 2, j, l, nw) = kappa_l * elec_dist_deriv_e(2, i, j)
+          een_rescaled_e_deriv_e(i, 3, j, l, nw) = kappa_l * elec_dist_deriv_e(3, i, j)
+          een_rescaled_e_deriv_e(i, 4, j, l, nw) = kappa_l * elec_dist_deriv_e(4, i, j)
 
-          een_rescaled_e_deriv_e(l, i, 4, j, nw) = een_rescaled_e_deriv_e(l, i, 4, j, nw)              &
-                    + een_rescaled_e_deriv_e(l, i, 1, j, nw) * een_rescaled_e_deriv_e(l, i, 1, j, nw)  &
-                    + een_rescaled_e_deriv_e(l, i, 2, j, nw) * een_rescaled_e_deriv_e(l, i, 2, j, nw)  &
-                    + een_rescaled_e_deriv_e(l, i, 3, j, nw) * een_rescaled_e_deriv_e(l, i, 3, j, nw)
+          een_rescaled_e_deriv_e(i, 4, j, l, nw) = een_rescaled_e_deriv_e(i, 4, j, l, nw)              &
+                    + een_rescaled_e_deriv_e(i, 1, j, l, nw) * een_rescaled_e_deriv_e(i, 1, j, l, nw)  &
+                    + een_rescaled_e_deriv_e(i, 2, j, l, nw) * een_rescaled_e_deriv_e(i, 2, j, l, nw)  &
+                    + een_rescaled_e_deriv_e(i, 3, j, l, nw) * een_rescaled_e_deriv_e(i, 3, j, l, nw)
 
-          do ii = 1, 4
-            een_rescaled_e_deriv_e(l, i, ii, j, nw) = een_rescaled_e_deriv_e(l, i, ii, j, nw) *   &
-                                                      een_rescaled_e(l, i, j, nw)
-          end do
+          een_rescaled_e_deriv_e(i, 1, j, l, nw) = een_rescaled_e_deriv_e(i, 1, j, l, nw) *   &
+                                                    een_rescaled_e(i, j, l, nw)
+          een_rescaled_e_deriv_e(i, 3, j, l, nw) = een_rescaled_e_deriv_e(i, 2, j, l, nw) *   &
+                                                    een_rescaled_e(i, j, l, nw)
+          een_rescaled_e_deriv_e(i, 3, j, l, nw) = een_rescaled_e_deriv_e(i, 3, j, l, nw) *   &
+                                                    een_rescaled_e(i, j, l, nw)
+          een_rescaled_e_deriv_e(i, 4, j, l, nw) = een_rescaled_e_deriv_e(i, 4, j, l, nw) *   &
+                                                    een_rescaled_e(i, j, l, nw)
         end do
       end do
     end do
@@ -3449,8 +3507,8 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f
       real    (c_double ) , intent(in)  , value :: rescale_factor_kappa_ee
       real    (c_double ) , intent(in)          :: coord_new(elec_num,3,walk_num)
       real    (c_double ) , intent(in)          :: ee_distance(elec_num,elec_num,walk_num)
-      real    (c_double ) , intent(in)          :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num)
-      real    (c_double ) , intent(out)         :: een_rescaled_e_deriv_e(0:cord_num,elec_num,4,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_deriv_e(elec_num,4,elec_num,0:cord_num,walk_num)
 
       integer(c_int32_t), external :: qmckl_compute_factor_een_rescaled_e_deriv_e_f
       info = qmckl_compute_factor_een_rescaled_e_deriv_e_f &
@@ -3550,16 +3608,17 @@ for l in range(0,cord_num+1):
 
      #+begin_src c :tangle (eval c_test)
 //assert(qmckl_electron_provided(context));
-double een_rescaled_e_deriv_e[walk_num][elec_num][4][elec_num][(cord_num + 1)];
-rc = qmckl_get_jastrow_een_rescaled_e_deriv_e(context, &(een_rescaled_e_deriv_e[0][0][0][0][0]));
+double een_rescaled_e_deriv_e[walk_num][(cord_num + 1)][elec_num][4][elec_num];
+size_max=0;
+rc = qmckl_get_jastrow_een_rescaled_e_deriv_e(context, &(een_rescaled_e_deriv_e[0][0][0][0][0]),&size_max);
 
 // value of (0,0,0,2,1)
-assert(fabs(een_rescaled_e_deriv_e[0][0][0][2][1] + 0.05991352796887283   ) < 1.e-12);
-assert(fabs(een_rescaled_e_deriv_e[0][0][0][3][1] + 0.011714035071545248  ) < 1.e-12);
-assert(fabs(een_rescaled_e_deriv_e[0][0][0][4][1] + 0.00441398875758468   ) < 1.e-12);
-assert(fabs(een_rescaled_e_deriv_e[0][1][0][3][2] + 0.013553180060167595  ) < 1.e-12);
-assert(fabs(een_rescaled_e_deriv_e[0][1][0][4][2] + 0.00041342909359870457) < 1.e-12);
-assert(fabs(een_rescaled_e_deriv_e[0][1][0][5][2] + 0.5880599146214673    ) < 1.e-12);
+assert(fabs(een_rescaled_e_deriv_e[0][1][0][0][2] + 0.05991352796887283   ) < 1.e-12);
+assert(fabs(een_rescaled_e_deriv_e[0][1][0][0][3] + 0.011714035071545248  ) < 1.e-12);
+assert(fabs(een_rescaled_e_deriv_e[0][1][0][0][4] + 0.00441398875758468   ) < 1.e-12);
+assert(fabs(een_rescaled_e_deriv_e[0][2][1][0][3] + 0.013553180060167595  ) < 1.e-12);
+assert(fabs(een_rescaled_e_deriv_e[0][2][1][0][4] + 0.00041342909359870457) < 1.e-12);
+assert(fabs(een_rescaled_e_deriv_e[0][2][1][0][5] + 0.5880599146214673    ) < 1.e-12);
      #+end_src
 
 ** Electron-nucleus rescaled distances for each order
@@ -3576,11 +3635,11 @@ assert(fabs(een_rescaled_e_deriv_e[0][1][0][5][2] + 0.5880599146214673    ) < 1.
 *** Get
 
     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
-qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(qmckl_context context, double* const distance_rescaled);
+qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(qmckl_context context, double* const distance_rescaled, int64_t* size_max);
     #+end_src
 
     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(qmckl_context context, double* const distance_rescaled)
+qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(qmckl_context context, double* const distance_rescaled, int64_t* size_max)
 {
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return QMCKL_NULL_CONTEXT;
@@ -3596,6 +3655,7 @@ qmckl_exit_code qmckl_get_jastrow_een_rescaled_n(qmckl_context context, double*
 
   size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
   memcpy(distance_rescaled, ctx->jastrow.een_rescaled_n, sze * sizeof(double));
+  (*size_max)=sze;
 
   return QMCKL_SUCCESS;
 }
@@ -3679,7 +3739,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context)
     | ~cord_num~                | ~int64_t~                                          | in     | Order of polynomials                |
     | ~rescale_factor_kappa_en~ | ~double~                                           | in     | Factor to rescale ee distances      |
     | ~en_distance~             | ~double[walk_num][elec_num][nucl_num]~             | in     | Electron-nucleus distances          |
-    | ~een_rescaled_n~          | ~double[walk_num][elec_num][nucl_num][0:cord_num]~ | out    | Electron-nucleus rescaled distances |
+    | ~een_rescaled_n~          | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | out    | Electron-nucleus rescaled distances |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
 integer function qmckl_compute_een_rescaled_n_f(context, walk_num, elec_num, nucl_num, cord_num, rescale_factor_kappa_en,  &
@@ -3694,7 +3754,7 @@ integer function qmckl_compute_een_rescaled_n_f(context, walk_num, elec_num, nuc
   integer*8             , intent(in)  :: cord_num
   double precision      , intent(in)  :: rescale_factor_kappa_en
   double precision      , intent(in)  :: en_distance(elec_num,nucl_num,walk_num)
-  double precision      , intent(out) :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num)
+  double precision      , intent(out) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
   double precision                    :: x
   integer*8                           :: i, a, k, l, nw
 
@@ -3730,18 +3790,18 @@ integer function qmckl_compute_een_rescaled_n_f(context, walk_num, elec_num, nuc
   do nw = 1, walk_num
 
   ! prepare the actual een table
-  een_rescaled_n(0, :, :, nw) = 1.0d0
+  een_rescaled_n(:, :, 0, nw) = 1.0d0
 
   do a = 1, nucl_num
     do i = 1, elec_num
-      een_rescaled_n(1, a, i, nw) = dexp(-rescale_factor_kappa_en * en_distance(i, a, nw))
+      een_rescaled_n(i, a, 1, nw) = dexp(-rescale_factor_kappa_en * en_distance(i, a, nw))
     end do
   end do
 
   do l = 2, cord_num
     do a = 1, nucl_num
       do i = 1, elec_num
-        een_rescaled_n(l, a, i, nw) = een_rescaled_n(l - 1, a, i, nw) * een_rescaled_n(1, a, i, nw)
+        een_rescaled_n(i, a, l, nw) = een_rescaled_n(i, a, l - 1, nw) * een_rescaled_n(i, a, 1, nw)
       end do
     end do
   end do
@@ -3790,7 +3850,7 @@ end function qmckl_compute_een_rescaled_n_f
       integer (c_int64_t) , intent(in)  , value :: cord_num
       real    (c_double ) , intent(in)  , value :: rescale_factor_kappa_en
       real    (c_double ) , intent(in)          :: en_distance(nucl_num,elec_num,walk_num)
-      real    (c_double ) , intent(out)         :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num)
+      real    (c_double ) , intent(out)         :: een_rescaled_n(nucl_num,elec_num,0:cord_num,walk_num)
 
       integer(c_int32_t), external :: qmckl_compute_een_rescaled_n_f
       info = qmckl_compute_een_rescaled_n_f &
@@ -3853,16 +3913,17 @@ print(" een_rescaled_n[1, 5, 2] = ",een_rescaled_n[1, 5, 2])
      #+begin_src c :tangle (eval c_test)
 assert(qmckl_electron_provided(context));
 
-double een_rescaled_n[walk_num][elec_num][nucl_num][(cord_num + 1)];
-rc = qmckl_get_jastrow_een_rescaled_n(context, &(een_rescaled_n[0][0][0][0]));
+double een_rescaled_n[walk_num][(cord_num + 1)][nucl_num][elec_num];
+size_max=0;
+rc = qmckl_get_jastrow_een_rescaled_n(context, &(een_rescaled_n[0][0][0][0]),&size_max);
 
 // value of (0,2,1)
-assert(fabs(een_rescaled_n[0][2][0][1]-0.10612983920006765)  < 1.e-12);
-assert(fabs(een_rescaled_n[0][3][0][1]-0.135652809635553)    < 1.e-12);
-assert(fabs(een_rescaled_n[0][4][0][1]-0.023391817607642338) < 1.e-12);
-assert(fabs(een_rescaled_n[0][3][1][2]-0.880957224822116)    < 1.e-12);
-assert(fabs(een_rescaled_n[0][4][1][2]-0.027185942659395074) < 1.e-12);
-assert(fabs(een_rescaled_n[0][5][1][2]-0.01343938025140174)  < 1.e-12);
+assert(fabs(een_rescaled_n[0][1][0][2]-0.10612983920006765)  < 1.e-12);
+assert(fabs(een_rescaled_n[0][1][0][3]-0.135652809635553)    < 1.e-12);
+assert(fabs(een_rescaled_n[0][1][0][4]-0.023391817607642338) < 1.e-12);
+assert(fabs(een_rescaled_n[0][2][1][3]-0.880957224822116)    < 1.e-12);
+assert(fabs(een_rescaled_n[0][2][1][4]-0.027185942659395074) < 1.e-12);
+assert(fabs(een_rescaled_n[0][2][1][5]-0.01343938025140174)  < 1.e-12);
 
      #+end_src
 
@@ -3875,11 +3936,11 @@ assert(fabs(een_rescaled_n[0][5][1][2]-0.01343938025140174)  < 1.e-12);
 *** Get
 
     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
-qmckl_exit_code qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context, double* const distance_rescaled);
+qmckl_exit_code qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max);
     #+end_src
 
     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context, double* const distance_rescaled)
+qmckl_exit_code qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context, double* const distance_rescaled, int64_t* size_max)
 {
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return QMCKL_NULL_CONTEXT;
@@ -3895,6 +3956,7 @@ qmckl_exit_code qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context,
 
   size_t sze = ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
   memcpy(distance_rescaled, ctx->jastrow.een_rescaled_n_deriv_e, sze * sizeof(double));
+  (*size_max)=sze;
 
   return QMCKL_SUCCESS;
 }
@@ -3987,8 +4049,8 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context)
     | ~coord_new~               | ~double[walk_num][3][elec_num]~                       | in     | Electron coordinates                |
     | ~coord~                   | ~double[3][nucl_num]~                                 | in     | Nuclear coordinates                 |
     | ~en_distance~             | ~double[walk_num][elec_num][nucl_num]~                | in     | Electron-nucleus distances          |
-    | ~een_rescaled_n~          | ~double[walk_num][elec_num][nucl_num][0:cord_num]~    | in     | Electron-nucleus distances          |
-    | ~een_rescaled_n_deriv_e~  | ~double[walk_num][elec_num][4][nucl_num][0:cord_num]~ | out    | Electron-nucleus rescaled distances |
+    | ~een_rescaled_n~          | ~double[walk_num][0:cord_num][nucl_num][elec_num]~    | in     | Electron-nucleus distances          |
+    | ~een_rescaled_n_deriv_e~  | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | out    | Electron-nucleus rescaled distances |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
 integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f(context, walk_num, elec_num, nucl_num, &
@@ -4006,8 +4068,8 @@ integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f(context, walk_num
   double precision      , intent(in)  :: coord_new(elec_num,3,walk_num)
   double precision      , intent(in)  :: coord(nucl_num,3)
   double precision      , intent(in)  :: en_distance(elec_num,nucl_num,walk_num)
-  double precision      , intent(in)  :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num)
-  double precision      , intent(out) :: een_rescaled_n_deriv_e(0:cord_num,nucl_num,4,elec_num,walk_num)
+  double precision      , intent(in)  :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
+  double precision      , intent(out) :: een_rescaled_n_deriv_e(elec_num,4,nucl_num,0:cord_num,walk_num)
   double precision,dimension(:,:,:),allocatable :: elnuc_dist_deriv_e
   double precision                    :: x, ria_inv, kappa_l
   integer*8                           :: i, a, k, l, nw, ii
@@ -4060,19 +4122,24 @@ integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f(context, walk_num
     kappa_l = - dble(l) * rescale_factor_kappa_en
     do a = 1, nucl_num
       do i = 1, elec_num
-        do ii = 1, 4
-          een_rescaled_n_deriv_e(l, a, ii, i, nw) = kappa_l * elnuc_dist_deriv_e(ii, i, a)
-        end do
+        een_rescaled_n_deriv_e(i, 1, a, l, nw) = kappa_l * elnuc_dist_deriv_e(1, i, a)
+        een_rescaled_n_deriv_e(i, 2, a, l, nw) = kappa_l * elnuc_dist_deriv_e(2, i, a)
+        een_rescaled_n_deriv_e(i, 3, a, l, nw) = kappa_l * elnuc_dist_deriv_e(3, i, a)
+        een_rescaled_n_deriv_e(i, 4, a, l, nw) = kappa_l * elnuc_dist_deriv_e(4, i, a)
 
-        een_rescaled_n_deriv_e(l, a, 4, i, nw) = een_rescaled_n_deriv_e(l, a, 4, i, nw)           &
-                + een_rescaled_n_deriv_e(l, a, 1, i, nw) * een_rescaled_n_deriv_e(l, a, 1, i, nw) &
-                + een_rescaled_n_deriv_e(l, a, 2, i, nw) * een_rescaled_n_deriv_e(l, a, 2, i, nw) &
-                + een_rescaled_n_deriv_e(l, a, 3, i, nw) * een_rescaled_n_deriv_e(l, a, 3, i, nw)
+        een_rescaled_n_deriv_e(i, 4, a, l, nw) = een_rescaled_n_deriv_e(i, 4, a, l, nw)           &
+                + een_rescaled_n_deriv_e(i, 1, a, l, nw) * een_rescaled_n_deriv_e(i, 1, a, l, nw) &
+                + een_rescaled_n_deriv_e(i, 2, a, l, nw) * een_rescaled_n_deriv_e(i, 2, a, l, nw) &
+                + een_rescaled_n_deriv_e(i, 3, a, l, nw) * een_rescaled_n_deriv_e(i, 3, a, l, nw)
 
-        do ii = 1, 4
-          een_rescaled_n_deriv_e(l, a, ii, i, nw) = een_rescaled_n_deriv_e(l, a, ii, i, nw) * &
-                                                    een_rescaled_n(l, a, i, nw)
-        end do
+        een_rescaled_n_deriv_e(i, 1, a, l, nw) = een_rescaled_n_deriv_e(i, 1, a, l, nw) * &
+                                                  een_rescaled_n(i, a, l, nw)
+        een_rescaled_n_deriv_e(i, 2, a, l, nw) = een_rescaled_n_deriv_e(i, 2, a, l, nw) * &
+                                                  een_rescaled_n(i, a, l, nw)
+        een_rescaled_n_deriv_e(i, 3, a, l, nw) = een_rescaled_n_deriv_e(i, 3, a, l, nw) * &
+                                                  een_rescaled_n(i, a, l, nw)
+        een_rescaled_n_deriv_e(i, 4, a, l, nw) = een_rescaled_n_deriv_e(i, 4, a, l, nw) * &
+                                                  een_rescaled_n(i, a, l, nw)
       end do
     end do
   end do
@@ -4130,7 +4197,7 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f
       real    (c_double ) , intent(in)          :: coord(nucl_num,3)
       real    (c_double ) , intent(in)          :: en_distance(nucl_num,elec_num,walk_num)
       real    (c_double ) , intent(in)          :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num)
-      real    (c_double ) , intent(out)         :: een_rescaled_n_deriv_e(0:cord_num,nucl_num,4,elec_num,walk_num)
+      real    (c_double ) , intent(out)         :: een_rescaled_n_deriv_e(elec_num,4,nucl_num,0:cord_num,walk_num)
 
       integer(c_int32_t), external :: qmckl_compute_factor_een_rescaled_n_deriv_e_f
       info = qmckl_compute_factor_een_rescaled_n_deriv_e_f &
@@ -4148,9 +4215,9 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f
 
     end function qmckl_compute_factor_een_rescaled_n_deriv_e
     #+end_src
-
+    
 *** Test
-
+    
     #+begin_src python :results output :exports none :noweb yes
 import numpy as np
 
@@ -4219,16 +4286,17 @@ print(" een_rescaled_n_deriv_e[2, 1, 6, 2] = ",een_rescaled_n_deriv_e[5, 0, 1, 2
      #+begin_src c :tangle (eval c_test)
 assert(qmckl_electron_provided(context));
 
-double een_rescaled_n_deriv_e[walk_num][elec_num][4][nucl_num][(cord_num + 1)];
-rc = qmckl_get_jastrow_een_rescaled_n_deriv_e(context, &(een_rescaled_n_deriv_e[0][0][0][0][0]));
+double een_rescaled_n_deriv_e[walk_num][(cord_num + 1)][nucl_num][4][elec_num];
+size_max=0;
+rc = qmckl_get_jastrow_een_rescaled_n_deriv_e(context, &(een_rescaled_n_deriv_e[0][0][0][0][0]),&size_max);
 
 // value of (0,2,1)
-assert(fabs(een_rescaled_n_deriv_e[0][2][0][0][1]+0.07633444246999128   )  < 1.e-12);
-assert(fabs(een_rescaled_n_deriv_e[0][3][0][0][1]-0.00033282346259738276)  < 1.e-12);
-assert(fabs(een_rescaled_n_deriv_e[0][4][0][0][1]+0.004775370547333061  )  < 1.e-12);
-assert(fabs(een_rescaled_n_deriv_e[0][3][0][1][2]-0.1362654644223866    )  < 1.e-12);
-assert(fabs(een_rescaled_n_deriv_e[0][4][0][1][2]+0.0231253431662794    )  < 1.e-12);
-assert(fabs(een_rescaled_n_deriv_e[0][5][0][1][2]-0.001593334817691633  )  < 1.e-12);
+assert(fabs(een_rescaled_n_deriv_e[0][1][0][0][2]+0.07633444246999128   )  < 1.e-12);
+assert(fabs(een_rescaled_n_deriv_e[0][1][0][0][3]-0.00033282346259738276)  < 1.e-12);
+assert(fabs(een_rescaled_n_deriv_e[0][1][0][0][4]+0.004775370547333061  )  < 1.e-12);
+assert(fabs(een_rescaled_n_deriv_e[0][2][1][0][3]-0.1362654644223866    )  < 1.e-12);
+assert(fabs(een_rescaled_n_deriv_e[0][2][1][0][4]+0.0231253431662794    )  < 1.e-12);
+assert(fabs(een_rescaled_n_deriv_e[0][2][1][0][5]-0.001593334817691633  )  < 1.e-12);
 
      #+end_src
 
@@ -4926,8 +4994,8 @@ end function qmckl_compute_lkpm_combined_index_f
     | ~elec_num~       | ~int64_t~                                                        | in     | Number of electrons               |
     | ~nucl_num~       | ~int64_t~                                                        | in     | Number of nucleii                 |
     | ~walk_num~       | ~int64_t~                                                        | in     | Number of walkers                 |
-    | ~een_rescaled_e~ | ~double[walk_num][elec_num][elec_num][0:cord_num]~               | in     | Electron-electron rescaled factor |
-    | ~een_rescaled_n~ | ~double[walk_num][elec_num][nucl_num][0:cord_num]~               | in     | Electron-nucleus rescaled factor  |
+    | ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~               | in     | Electron-electron rescaled factor |
+    | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~               | in     | Electron-nucleus rescaled factor  |
     | ~tmp_c~          | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | out    | vector of non-zero coefficients   |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
@@ -4941,12 +5009,10 @@ integer function qmckl_compute_tmp_c_f(context, cord_num, elec_num, nucl_num, &
   integer*8             , intent(in)  :: elec_num
   integer*8             , intent(in)  :: nucl_num
   integer*8             , intent(in)  :: walk_num
-  double precision      , intent(in)  :: een_rescaled_e(0:cord_num, elec_num, elec_num, walk_num)
-  double precision      , intent(in)  :: een_rescaled_n(0:cord_num, nucl_num, elec_num, walk_num)
+  double precision      , intent(in)  :: een_rescaled_e(elec_num, elec_num, 0:cord_num, walk_num)
+  double precision      , intent(in)  :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
   double precision      , intent(out) :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num)
   double precision                    :: x
-  double precision,dimension(:,:,:,:),allocatable :: een_rescaled_e_T
-  double precision,dimension(:,:,:,:),allocatable :: een_rescaled_n_T
   integer*8                           :: i, j, a, l, kk, p, lmax, nw
   character                           :: TransA, TransB
   double precision                    :: alpha, beta
@@ -4957,27 +5023,6 @@ integer function qmckl_compute_tmp_c_f(context, cord_num, elec_num, nucl_num, &
   alpha = 1.0d0
   beta  = 0.0d0
 
-  allocate(een_rescaled_e_T(elec_num,elec_num,0:cord_num,walk_num))
-  allocate(een_rescaled_n_T(elec_num,nucl_num,0:cord_num,walk_num))
-  do nw = 1,walk_num
-  do i = 1, elec_num
-    do j = 1, elec_num
-      do l = 0,cord_num
-        een_rescaled_e_T(i,j,l,nw) = een_rescaled_e(l,j,i,nw)
-      end do
-    end do
-  end do
-  end do
-  do nw = 1,walk_num
-  do i = 1, elec_num
-    do j = 1, nucl_num
-      do l = 0,cord_num
-        een_rescaled_n_T(i,j,l,nw) = een_rescaled_n(l,j,i,nw)
-      end do
-    end do
-  end do
-  end do
-
   info = QMCKL_SUCCESS
 
   if (context == QMCKL_NULL_CONTEXT) then
@@ -5003,29 +5048,20 @@ integer function qmckl_compute_tmp_c_f(context, cord_num, elec_num, nucl_num, &
   M = elec_num
   N = nucl_num*(cord_num + 1)
   K = elec_num
-  LDA = size(een_rescaled_e_T,1)
-  LDB = size(een_rescaled_n_T,1)
+  LDA = size(een_rescaled_e,1)
+  LDB = size(een_rescaled_n,1)
   LDC = size(tmp_c,1)
 
   do nw=1, walk_num
   do i=0, cord_num-1
   info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha,     &
-                     een_rescaled_e_T(1,1,i,nw),LDA*1_8,                     &
-                     een_rescaled_n_T(1,1,0,nw),LDB*1_8,                     &
+                     een_rescaled_e(1,1,i,nw),LDA*1_8,                     &
+                     een_rescaled_n(1,1,0,nw),LDB*1_8,                     &
                      beta,                                       &
                      tmp_c(1,1,0,i,nw),LDC)
   end do
   end do
-  !do kk=0, cord_num-1
-  !do i=1,nucl_num
-  !  do j=1,elec_num
-  !    print *,tmp_c(j,i,:,kk,1)
-  !  end do
-  !end do
-  !end do
 
-  deallocate(een_rescaled_e_T)
-  deallocate(een_rescaled_n_T)
 end function qmckl_compute_tmp_c_f
     #+end_src
 
@@ -5061,8 +5097,8 @@ end function qmckl_compute_tmp_c_f
       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 :: walk_num
-      real    (c_double ) , intent(in)          :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num)
-      real    (c_double ) , intent(in)          :: een_rescaled_n(0:cord_num,nucl_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(in)          :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
       real    (c_double ) , intent(out)         :: tmp_c(elec_num,nucl_num,0:cord_num,0:cord_num-1,walk_num)
 
       integer(c_int32_t), external :: qmckl_compute_tmp_c_f
@@ -5087,8 +5123,8 @@ end function qmckl_compute_tmp_c_f
     | ~elec_num~               | ~int64_t~                                                        | in     | Number of electrons                           |
     | ~nucl_num~               | ~int64_t~                                                        | in     | Number of nucleii                             |
     | ~walk_num~               | ~int64_t~                                                        | in     | Number of walkers                             |
-    | ~een_rescaled_e_deriv_e~ | ~double[walk_num][elec_num][4][elec_num][0:cord_num]~            | in     | Electron-electron rescaled factor derivatives |
-    | ~een_rescaled_n~         | ~double[walk_num][elec_num][nucl_num][0:cord_num]~               | in     | Electron-nucleus rescaled factor              |
+    | ~een_rescaled_e_deriv_e~ | ~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               |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
@@ -5102,16 +5138,14 @@ integer function qmckl_compute_dtmp_c_f(context, cord_num, elec_num, nucl_num, &
   integer*8             , intent(in)  :: elec_num
   integer*8             , intent(in)  :: nucl_num
   integer*8             , intent(in)  :: walk_num
-  double precision      , intent(in)  :: een_rescaled_e_deriv_e(0:cord_num, elec_num, 4, elec_num, walk_num)
-  double precision      , intent(in)  :: een_rescaled_n(0:cord_num, nucl_num, elec_num, walk_num)
+  double precision      , intent(in)  :: een_rescaled_e_deriv_e(elec_num, 4, elec_num, 0:cord_num, walk_num)
+  double precision      , intent(in)  :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
   double precision      , intent(out) :: dtmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
   double precision                    :: x
   integer*8                           :: i, j, a, l, kk, p, lmax, nw, ii
   character                           :: TransA, TransB
   double precision                    :: alpha, beta
   integer*8                           :: M, N, K, LDA, LDB, LDC
-  double precision,dimension(:,:,:,:,:),allocatable :: een_rescaled_e_deriv_e_T
-  double precision,dimension(:,:,:,:),allocatable :: een_rescaled_n_T
 
   TransA = 'N'
   TransB = 'N'
@@ -5120,27 +5154,6 @@ integer function qmckl_compute_dtmp_c_f(context, cord_num, elec_num, nucl_num, &
 
   info = QMCKL_SUCCESS
 
-  allocate(een_rescaled_e_deriv_e_T(elec_num,4,elec_num,0:cord_num,walk_num))
-  allocate(een_rescaled_n_T(elec_num,nucl_num,0:cord_num,walk_num))
-  do nw = 1,walk_num
-  do i = 1, elec_num
-    do ii = 1, 4
-    do j = 1, elec_num
-      do l = 0,cord_num
-        een_rescaled_e_deriv_e_T(i,ii,j,l,nw) = een_rescaled_e_deriv_e(l,j,ii,i,nw)
-      end do
-    end do
-    end do
-  end do
-  do i = 1, elec_num
-    do j = 1, nucl_num
-      do l = 0,cord_num
-        een_rescaled_n_T(i,j,l,nw) = een_rescaled_n(l,j,i,nw)
-      end do
-    end do
-  end do
-  end do
-
   if (context == QMCKL_NULL_CONTEXT) then
      info = QMCKL_INVALID_CONTEXT
      return
@@ -5164,22 +5177,20 @@ integer function qmckl_compute_dtmp_c_f(context, cord_num, elec_num, nucl_num, &
   M = 4*elec_num
   N = nucl_num*(cord_num + 1)
   K = elec_num
-  LDA = 4*size(een_rescaled_e_deriv_e_T,1)
-  LDB = size(een_rescaled_n_T,1)
+  LDA = 4*size(een_rescaled_e_deriv_e,1)
+  LDB = size(een_rescaled_n,1)
   LDC = 4*size(dtmp_c,1)
 
   do nw=1, walk_num
   do i=0, cord_num-1
   info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha,     &
-                     een_rescaled_e_deriv_e_T(1,1,1,i,nw),LDA*1_8,                     &
-                     een_rescaled_n_T(1,1,0,nw),LDB*1_8,                     &
+                     een_rescaled_e_deriv_e(1,1,1,i,nw),LDA*1_8,                     &
+                     een_rescaled_n(1,1,0,nw),LDB*1_8,                     &
                      beta,                                       &
                      dtmp_c(1,1,1,0,i,nw),LDC)
   end do
   end do
 
-  deallocate(een_rescaled_e_deriv_e_T)
-  deallocate(een_rescaled_n_T)
 end function qmckl_compute_dtmp_c_f
     #+end_src
 
@@ -5215,8 +5226,8 @@ end function qmckl_compute_dtmp_c_f
       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 :: walk_num
-      real    (c_double ) , intent(in)          :: een_rescaled_e_deriv_e(0:cord_num,elec_num,4,elec_num,walk_num)
-      real    (c_double ) , intent(in)          :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num)
+      real    (c_double ) , intent(in)          :: een_rescaled_e_deriv_e(elec_num,4,elec_num,0:cord_num,walk_num)
+      real    (c_double ) , intent(in)          :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
       real    (c_double ) , intent(out)         :: dtmp_c(elec_num,nucl_num,0:cord_num,0:cord_num-1,walk_num)
 
       integer(c_int32_t), external :: qmckl_compute_dtmp_c_f
@@ -5307,7 +5318,7 @@ rc = qmckl_get_jastrow_dtmp_c(context, &(dtmp_c[0][0][0][0][0][0]));
 
 assert(fabs(tmp_c[0][0][1][0][0] - 2.7083473948352403) < 1e-12);
 
-assert(fabs(dtmp_c[0][1][0][0][0][0] + 0.237440520852232) < 1e-12);
+assert(fabs(dtmp_c[0][1][0][0][0][0] - 0.237440520852232) < 1e-12);
      #+end_src
 
 ** Electron-electron-nucleus Jastrow \(f_{een}\)
@@ -5319,11 +5330,11 @@ assert(fabs(dtmp_c[0][1][0][0][0][0] + 0.237440520852232) < 1e-12);
 
 *** Get
     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
-qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* const factor_een);
+qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* const factor_een, int64_t* size_max);
     #+end_src
 
     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* const factor_een)
+qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* const factor_een, int64_t* size_max)
 {
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return QMCKL_NULL_CONTEXT;
@@ -5339,6 +5350,7 @@ qmckl_exit_code qmckl_get_jastrow_factor_een(qmckl_context context, double* cons
 
   int64_t sze = ctx->electron.walk_num * ctx->electron.num;
   memcpy(factor_een, ctx->jastrow.factor_een, sze*sizeof(double));
+  (*size_max)=sze;
 
   return QMCKL_SUCCESS;
 }
@@ -5610,7 +5622,7 @@ end function qmckl_compute_factor_een_naive_f
     | ~cord_vect_full~      | ~double[dim_cord_vect][nucl_num]~                                | in                              | full coefficient vector              |
     | ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~                                      | in                              | combined indices                     |
     | ~tmp_c~               | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | vector of non-zero coefficients |                                      |
-    | ~een_rescaled_n~      | ~double[walk_num][elec_num][nucl_num][0:cord_num]~               | in                              | Electron-nucleus rescaled factor     |
+    | ~een_rescaled_n~      | ~double[walk_num][0:cord_num][nucl_num][elec_num]~               | in                              | Electron-nucleus rescaled factor     |
     | ~factor_een~          | ~double[walk_num]~                                               | out                             | Electron-nucleus jastrow             |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
@@ -5625,7 +5637,7 @@ integer function qmckl_compute_factor_een_f(context, walk_num, elec_num, nucl_nu
   integer*8             , intent(in)  :: lkpm_combined_index(dim_cord_vect,4)
   double precision      , intent(in)  :: cord_vect_full(nucl_num, dim_cord_vect)
   double precision      , intent(in)  :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
-  double precision      , intent(in)  :: een_rescaled_n(0:cord_num, nucl_num, elec_num, walk_num)
+  double precision      , intent(in)  :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
   double precision      , intent(out) :: factor_een(walk_num)
 
   integer*8 :: i, a, j, l, k, p, m, n, nw
@@ -5673,7 +5685,7 @@ integer function qmckl_compute_factor_een_f(context, walk_num, elec_num, nucl_nu
 
       accu = 0.0d0
       do j = 1, elec_num
-        accu = accu + een_rescaled_n(m,a,j,nw) * tmp_c(j,a,m+l,k,nw)
+        accu = accu + een_rescaled_n(j,a,m,nw) * tmp_c(j,a,m+l,k,nw)
       end do
       factor_een(nw) = factor_een(nw) + accu * cn
     end do
@@ -5731,7 +5743,7 @@ end function qmckl_compute_factor_een_f
       real    (c_double ) , intent(in)          :: cord_vect_full(nucl_num,dim_cord_vect)
       integer (c_int64_t) , intent(in)          :: lkpm_combined_index(dim_cord_vect,4)
       real    (c_double ) , intent(in)          :: een_rescaled_e(0:cord_num,elec_num,elec_num,walk_num)
-      real    (c_double ) , intent(in)          :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num)
+      real    (c_double ) , intent(in)          :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
       real    (c_double ) , intent(out)         :: factor_een(walk_num)
 
       integer(c_int32_t), external :: qmckl_compute_factor_een_f
@@ -5792,7 +5804,8 @@ print("factor_een:",factor_een)
 assert(qmckl_jastrow_provided(context));
 
 double factor_een[walk_num];
-rc = qmckl_get_jastrow_factor_een(context, &(factor_een[0]));
+size_max=0;
+rc = qmckl_get_jastrow_factor_een(context, &(factor_een[0]),&size_max);
 
 assert(fabs(factor_een[0] + 0.37407972141304213) < 1e-12);
      #+end_src
@@ -5806,11 +5819,11 @@ assert(fabs(factor_een[0] + 0.37407972141304213) < 1e-12);
 
 *** Get
     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
-qmckl_exit_code qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, double* const factor_een_deriv_e);
+qmckl_exit_code qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, double* const factor_een_deriv_e, int64_t* size_max);
     #+end_src
 
     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, double* const factor_een_deriv_e)
+qmckl_exit_code qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, double* const factor_een_deriv_e, int64_t* size_max)
 {
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return QMCKL_NULL_CONTEXT;
@@ -5826,6 +5839,7 @@ qmckl_exit_code qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context, doub
 
   int64_t sze = ctx->electron.walk_num * ctx->electron.num;
   memcpy(factor_een_deriv_e, ctx->jastrow.factor_een_deriv_e, sze*sizeof(double));
+  (*size_max)=sze;
 
   return QMCKL_SUCCESS;
 }
@@ -6143,8 +6157,8 @@ end function qmckl_compute_factor_een_deriv_e_naive_f
     | ~lkpm_combined_index~    | ~int64_t[4][dim_cord_vect]~                                         | in     | combined indices                               |
     | ~tmp_c~                  | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~    | in     | Temporary intermediate tensor                  |
     | ~dtmp_c~                 | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]~ | in     | vector of non-zero coefficients                |
-    | ~een_rescaled_n~         | ~double[walk_num][elec_num][nucl_num][0:cord_num]~                  | in     | Electron-nucleus rescaled factor               |
-    | ~een_rescaled_n_deriv_e~ | ~double[walk_num][elec_num][4][nucl_num][0:cord_num]~               | in     | Derivative of Electron-nucleus rescaled factor |
+    | ~een_rescaled_n~         | ~double[walk_num][0:cord_num][nucl_num][elec_num]~                  | in     | Electron-nucleus rescaled factor               |
+    | ~een_rescaled_n_deriv_e~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~               | in     | Derivative of Electron-nucleus rescaled factor |
     | ~factor_een_deriv_e~     | ~double[walk_num][4][elec_num]~                                     | out    | Derivative of Electron-nucleus jastrow         |
 
 
@@ -6161,8 +6175,8 @@ integer function qmckl_compute_factor_een_deriv_e_f(context, walk_num, elec_num,
   double precision      , intent(in)  :: cord_vect_full(nucl_num, dim_cord_vect)
   double precision      , intent(in)  :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
   double precision      , intent(in)  :: dtmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1,  walk_num)
-  double precision      , intent(in)  :: een_rescaled_n(0:cord_num, nucl_num, elec_num, walk_num)
-  double precision      , intent(in)  :: een_rescaled_n_deriv_e(0:cord_num, nucl_num, 4, elec_num, walk_num)
+  double precision      , intent(in)  :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num)
+  double precision      , intent(in)  :: een_rescaled_n_deriv_e(elec_num, 4, nucl_num, 0:cord_num, walk_num)
   double precision      , intent(out) :: factor_een_deriv_e(elec_num,4,walk_num)
 
   integer*8 :: i, a, j, l, k, p, m, n, nw, ii
@@ -6211,10 +6225,10 @@ integer function qmckl_compute_factor_een_deriv_e_f(context, walk_num, elec_num,
       do ii = 1, 4
         do j = 1, elec_num
           factor_een_deriv_e(j,ii,nw) = factor_een_deriv_e(j,ii,nw) +                           (&
-                                        tmp_c(j,a,m,k,nw)       * een_rescaled_n_deriv_e(m+l,a,ii,j,nw) + &
-                                        (-1.0d0*dtmp_c(j,ii,a,m,k,nw))   * een_rescaled_n(m+l,a,j,nw)            + &
-                                        (-1.0d0*dtmp_c(j,ii,a,m+l,k,nw)) * een_rescaled_n(m,a,j,nw)              + &
-                                        tmp_c(j,a,m+l,k,nw)     * een_rescaled_n_deriv_e(m,a,ii,j,nw)     &
+                                        tmp_c(j,a,m,k,nw)       * een_rescaled_n_deriv_e(j,ii,a,m+l,nw) + &
+                                        (dtmp_c(j,ii,a,m,k,nw))   * een_rescaled_n(j,a,m+l,nw)            + &
+                                        (dtmp_c(j,ii,a,m+l,k,nw)) * een_rescaled_n(j,a,m  ,nw)              + &
+                                        tmp_c(j,a,m+l,k,nw)     * een_rescaled_n_deriv_e(j,ii,a,m,nw)     &
                                         ) * cn
         end do
       end do
@@ -6222,12 +6236,12 @@ integer function qmckl_compute_factor_een_deriv_e_f(context, walk_num, elec_num,
       cn = cn + cn
       do j = 1, elec_num
         factor_een_deriv_e(j,4,nw) = factor_een_deriv_e(j,4,nw) +                              (&
-                                      (-1.0d0*dtmp_c(j,1,a,m  ,k,nw)) * een_rescaled_n_deriv_e(m+l,a,1,j,nw)  + &
-                                      (-1.0d0*dtmp_c(j,2,a,m  ,k,nw)) * een_rescaled_n_deriv_e(m+l,a,2,j,nw)  + &
-                                      (-1.0d0*dtmp_c(j,3,a,m  ,k,nw)) * een_rescaled_n_deriv_e(m+l,a,3,j,nw)  + &
-                                      (-1.0d0*dtmp_c(j,1,a,m+l,k,nw)) * een_rescaled_n_deriv_e(m  ,a,1,j,nw)  + &
-                                      (-1.0d0*dtmp_c(j,2,a,m+l,k,nw)) * een_rescaled_n_deriv_e(m  ,a,2,j,nw)  + &
-                                      (-1.0d0*dtmp_c(j,3,a,m+l,k,nw)) * een_rescaled_n_deriv_e(m  ,a,3,j,nw)    &
+                                      (dtmp_c(j,1,a,m  ,k,nw)) * een_rescaled_n_deriv_e(j,1,a,m+l,nw)  + &
+                                      (dtmp_c(j,2,a,m  ,k,nw)) * een_rescaled_n_deriv_e(j,2,a,m+l,nw)  + &
+                                      (dtmp_c(j,3,a,m  ,k,nw)) * een_rescaled_n_deriv_e(j,3,a,m+l,nw)  + &
+                                      (dtmp_c(j,1,a,m+l,k,nw)) * een_rescaled_n_deriv_e(j,1,a,m  ,nw)  + &
+                                      (dtmp_c(j,2,a,m+l,k,nw)) * een_rescaled_n_deriv_e(j,2,a,m  ,nw)  + &
+                                      (dtmp_c(j,3,a,m+l,k,nw)) * een_rescaled_n_deriv_e(j,3,a,m  ,nw)    &
                                       ) * cn
       end do
     end do
@@ -6291,8 +6305,8 @@ end function qmckl_compute_factor_een_deriv_e_f
       integer (c_int64_t) , intent(in)          :: lkpm_combined_index(dim_cord_vect,4)
       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)          :: een_rescaled_n(0:cord_num,nucl_num,elec_num,walk_num)
-      real    (c_double ) , intent(in)          :: een_rescaled_n_deriv_e(0:cord_num,nucl_num,4,elec_num,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_n_deriv_e(elec_num,4,nucl_num,0:cord_num,walk_num)
       real    (c_double ) , intent(out)         :: factor_een_deriv_e(elec_num,4,walk_num)
 
       integer(c_int32_t), external :: qmckl_compute_factor_een_deriv_e_f
@@ -6313,7 +6327,7 @@ end function qmckl_compute_factor_een_deriv_e_f
 
     end function qmckl_compute_factor_een_deriv_e
     #+end_src
-
+    
 *** Test
     #+begin_src python :results output :exports none :noweb yes
 import numpy as np
@@ -6371,7 +6385,8 @@ print("factor_een:",factor_een)
 assert(qmckl_jastrow_provided(context));
 
 double factor_een_deriv_e[walk_num][elec_num];
-rc = qmckl_get_jastrow_factor_een_deriv_e(context, &(factor_een_deriv_e[0][0]));
+size_max=0;
+rc = qmckl_get_jastrow_factor_een_deriv_e(context, &(factor_een_deriv_e[0][0]),&size_max);
 
 assert(fabs(factor_een_deriv_e[0][0] + 0.0005481671107226865) < 1e-12);
      #+end_src