diff --git a/configure.ac b/configure.ac
index d2ac039..ffa8f99 100644
--- a/configure.ac
+++ b/configure.ac
@@ -200,16 +200,22 @@ AS_IF([test "$BLAS_LIBS" == "$LAPACK_LIBS"], [BLAS_LIBS=""])
 # Specific options required with some compilers
 case $FC in
 
- ifort*)
+ *ifort*)
       FCFLAGS="$FCFLAGS -nofor-main"
       ;;
 
-#gfortran*)
-#     # Order is important here
-#     FCFLAGS="-cpp $FCFLAGS"
-#     ;;
+ *nvfortran*)
+      FCFLAGS="$FCFLAGS -fPIC -Mnomain -mp -target=gpu"
+      ;;
+
 esac
 
+case $CC in
+
+  *nvc*)
+        CFLAGS="$CFLAGS -fPIC -mp -target=gpu"
+        ;;
+esac
 
 
 # Options.
@@ -384,7 +390,7 @@ CC..............: ${CC}
 CPPFLAGS........: ${CPPFLAGS}
 CFLAGS..........: ${CFLAGS}
 FC..............: ${FC}
-FCLAGS..........: ${FCFLAGS}
+FCFLAGS.........: ${FCFLAGS}
 LDFLAGS:........: ${LDFLAGS}
 LIBS............: ${LIBS}
 USE CHAMELEON...: ${with_chameleon}
diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org
index 5bc9235..8dd1067 100644
--- a/org/qmckl_ao.org
+++ b/org/qmckl_ao.org
@@ -139,6 +139,7 @@ int main() {
    | ~ao_num~            | ~int64_t~            | Number of AOs                                                        |
    | ~ao_cartesian~      | ~bool~               | If true, use polynomials. Otherwise, use spherical harmonics         |
    | ~ao_factor~         | ~double[ao_num]~     | Normalization factor of the AO                                       |
+   |---------------------+----------------------+----------------------------------------------------------------------|
 
     For H_2 with the following basis set,
 
@@ -330,7 +331,7 @@ qmckl_exit_code qmckl_init_ao_basis(qmckl_context context) {
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   ctx->ao_basis.uninitialized = (1 << 14) - 1;
@@ -358,7 +359,7 @@ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
                            NULL);
  }
 
-qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
      #+end_src
 
      #+NAME:post2
@@ -1375,7 +1376,7 @@ qmckl_get_ao_basis_type (const qmckl_context context,
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1;
@@ -1415,7 +1416,7 @@ qmckl_get_ao_basis_shell_num (const qmckl_context context,
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 1;
@@ -1453,7 +1454,7 @@ qmckl_get_ao_basis_prim_num (const qmckl_context context,
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 2;
@@ -1494,7 +1495,7 @@ qmckl_get_ao_basis_nucleus_shell_num (const qmckl_context context,
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 3;
@@ -1549,7 +1550,7 @@ qmckl_get_ao_basis_nucleus_index (const qmckl_context context,
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 4;
@@ -1605,7 +1606,7 @@ qmckl_get_ao_basis_shell_ang_mom (const qmckl_context context,
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 5;
@@ -1661,7 +1662,7 @@ qmckl_get_ao_basis_shell_prim_num (const qmckl_context context,
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 6;
@@ -1717,7 +1718,7 @@ qmckl_get_ao_basis_shell_prim_index (const qmckl_context context,
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 7;
@@ -1771,7 +1772,7 @@ qmckl_get_ao_basis_shell_factor (const qmckl_context context,
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 8;
@@ -1827,7 +1828,7 @@ qmckl_get_ao_basis_exponent (const qmckl_context context,
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 9;
@@ -1881,7 +1882,7 @@ qmckl_get_ao_basis_coefficient (const qmckl_context context,
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 10;
@@ -1936,7 +1937,7 @@ qmckl_get_ao_basis_prim_factor (const qmckl_context context,
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 11;
@@ -1989,7 +1990,7 @@ qmckl_get_ao_basis_ao_num (const qmckl_context context,
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 12;
@@ -2031,7 +2032,7 @@ qmckl_get_ao_basis_ao_factor (const qmckl_context context,
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 13;
@@ -2078,7 +2079,7 @@ bool qmckl_ao_basis_provided(const qmckl_context context) {
     return false;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   return ctx->ao_basis.provided;
@@ -2491,6 +2492,7 @@ free(ao_factor_test);
    | ~shell_vgl_date~     | ~uint64_t~                        | Last modification date of Value, gradients, Laplacian of the AOs at current positions        |
    | ~ao_vgl~             | ~double[point_num][5][ao_num]~    | Value, gradients, Laplacian of the primitives at current  positions                          |
    | ~ao_vgl_date~        | ~uint64_t~                        | Last modification date of Value, gradients, Laplacian of the AOs at current positions        |
+   |----------------------+-----------------------------------+----------------------------------------------------------------------------------------------|
 
 *** After initialization
 
@@ -2516,7 +2518,7 @@ qmckl_exit_code qmckl_finalize_basis(qmckl_context context) {
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int64_t nucl_num = 0;
@@ -2694,7 +2696,7 @@ int compare_basis( const void * l, const void * r )
 #ifdef HAVE_HPC
 qmckl_exit_code qmckl_finalize_basis_hpc (qmckl_context context)
 {
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
 
   mem_info.size = ctx->nucleus.num * sizeof(int32_t);
@@ -2865,7 +2867,7 @@ qmckl_get_ao_basis_primitive_vgl (qmckl_context context,
   rc = qmckl_provide_ao_basis_primitive_vgl(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int64_t sze = ctx->ao_basis.prim_num * 5 * ctx->point.num;
@@ -2926,7 +2928,7 @@ qmckl_get_ao_basis_shell_vgl (qmckl_context context,
   rc = qmckl_provide_ao_basis_shell_vgl(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int64_t sze = ctx->ao_basis.shell_num * 5 * ctx->point.num;
@@ -2989,7 +2991,7 @@ qmckl_get_ao_basis_ao_vgl (qmckl_context context,
   rc = qmckl_provide_ao_vgl(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int64_t sze = ctx->ao_basis.ao_num * 5 * ctx->point.num;
@@ -3044,7 +3046,7 @@ qmckl_get_ao_basis_ao_vgl_inplace (qmckl_context context,
 
   qmckl_exit_code rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int64_t sze = ctx->ao_basis.ao_num * 5 * ctx->point.num;
@@ -3485,7 +3487,7 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if (!ctx->ao_basis.provided) {
@@ -3886,7 +3888,7 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if (!ctx->ao_basis.provided) {
@@ -4800,7 +4802,7 @@ qmckl_ao_polynomial_transp_vgl_hpc (const qmckl_context context,
                                     double* restrict const VGL,
                                     const int64_t ldv )
 {
-  const qmckl_context_struct* ctx = (qmckl_context_struct* const) context;
+  const qmckl_context_struct* ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL && X != NULL && R != NULL && n != NULL && L != NULL && VGL != NULL);
   if (lmax < 0) return QMCKL_INVALID_ARG_4;
   if (ldl < 3) return QMCKL_INVALID_ARG_7;
@@ -5910,7 +5912,7 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if (!ctx->ao_basis.provided) {
diff --git a/org/qmckl_context.org b/org/qmckl_context.org
index 6a35789..9054ba1 100644
--- a/org/qmckl_context.org
+++ b/org/qmckl_context.org
@@ -169,7 +169,7 @@ qmckl_context qmckl_context_check(const qmckl_context context) {
   if (context == QMCKL_NULL_CONTEXT)
     return QMCKL_NULL_CONTEXT;
 
-  const qmckl_context_struct* const ctx = (const qmckl_context_struct* const) context;
+  const qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
 
   /* Try to access memory */
   if (ctx->tag != VALID_TAG) {
@@ -209,7 +209,7 @@ qmckl_context_touch(const qmckl_context context)
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
 
   ctx->date += 1UL;
   ctx->point.date += 1UL;
@@ -234,7 +234,7 @@ qmckl_context qmckl_context_create();
 qmckl_context qmckl_context_create() {
 
   qmckl_context_struct* const ctx =
-    (qmckl_context_struct* const) malloc (sizeof(qmckl_context_struct));
+    (qmckl_context_struct*) malloc (sizeof(qmckl_context_struct));
 
   if (ctx == NULL) {
     return QMCKL_NULL_CONTEXT;
@@ -267,7 +267,7 @@ qmckl_context qmckl_context_create() {
   {
     ctx->tag = VALID_TAG;
 
-    const qmckl_context context = (const qmckl_context) ctx;
+    const qmckl_context context = (qmckl_context) ctx;
     assert ( qmckl_context_check(context) != QMCKL_NULL_CONTEXT );
 
     qmckl_exit_code rc;
@@ -350,7 +350,7 @@ void qmckl_unlock(qmckl_context context);
 void qmckl_lock(qmckl_context context) {
   if (context == QMCKL_NULL_CONTEXT)
     return ;
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   errno = 0;
   int rc = pthread_mutex_lock( &(ctx->mutex) );
   if (rc != 0) {
@@ -362,7 +362,7 @@ void qmckl_lock(qmckl_context context) {
 }
 
 void qmckl_unlock(const qmckl_context context) {
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   int rc = pthread_mutex_unlock( &(ctx->mutex) );
   if (rc != 0) {
     fprintf(stderr, "DEBUG qmckl_unlock:%s\n", strerror(rc) );
@@ -401,10 +401,10 @@ qmckl_context qmckl_context_copy(const qmckl_context context) {
   {
 
     const qmckl_context_struct* const old_ctx =
-      (qmckl_context_struct* const) checked_context;
+      (qmckl_context_struct*) checked_context;
 
     qmckl_context_struct* const new_ctx =
-      (qmckl_context_struct* const) malloc (context, sizeof(qmckl_context_struct));
+      (qmckl_context_struct*) malloc (context, sizeof(qmckl_context_struct));
 
     if (new_ctx == NULL) {
       qmckl_unlock(context);
@@ -467,7 +467,7 @@ qmckl_context_destroy (const qmckl_context context)
   const qmckl_context checked_context = qmckl_context_check(context);
   if (checked_context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);  /* Shouldn't be possible because the context is valid */
 
   qmckl_lock(context);
diff --git a/org/qmckl_determinant.org b/org/qmckl_determinant.org
index a28f49c..ec4c9d9 100644
--- a/org/qmckl_determinant.org
+++ b/org/qmckl_determinant.org
@@ -182,7 +182,7 @@ qmckl_exit_code qmckl_init_determinant(qmckl_context context) {
     return false;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   ctx->det.uninitialized = (1 << 6) - 1;
@@ -216,7 +216,7 @@ bool qmckl_determinant_provided(const qmckl_context context) {
     return false;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   return ctx->det.provided;
@@ -238,7 +238,7 @@ char qmckl_get_determinant_type (const qmckl_context context) {
     return (char) 0;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1;
@@ -256,7 +256,7 @@ int64_t qmckl_get_determinant_walk_num (const qmckl_context context) {
     return (int64_t) 0;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 1;
@@ -274,7 +274,7 @@ int64_t qmckl_get_determinant_det_num_alpha (const qmckl_context context) {
     return (int64_t) 0;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 2;
@@ -292,7 +292,7 @@ int64_t qmckl_get_determinant_det_num_beta (const qmckl_context context) {
     return (int64_t) 0;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 3;
@@ -310,7 +310,7 @@ int64_t* qmckl_get_determinant_mo_index_alpha (const qmckl_context context) {
     return NULL;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 4;
@@ -328,7 +328,7 @@ int64_t* qmckl_get_determinant_mo_index_beta (const qmckl_context context) {
     return NULL;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 5;
@@ -363,7 +363,7 @@ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
   return QMCKL_NULL_CONTEXT;
  }
 
-qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
    #+end_src
 
    #+NAME:post2
@@ -525,7 +525,7 @@ qmckl_exit_code qmckl_finalize_determinant(qmckl_context context) {
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   qmckl_exit_code rc;
@@ -596,7 +596,7 @@ qmckl_exit_code qmckl_get_det_vgl_alpha(qmckl_context context, double * const de
   rc = qmckl_provide_det_vgl_alpha(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = 5 * ctx->det.det_num_alpha * ctx->det.walk_num *
@@ -623,7 +623,7 @@ qmckl_exit_code qmckl_get_det_vgl_beta(qmckl_context context, double * const det
   rc = qmckl_provide_det_vgl_beta(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = 5 * ctx->det.det_num_beta * ctx->det.walk_num *
@@ -649,7 +649,7 @@ qmckl_exit_code qmckl_provide_det_vgl_alpha(qmckl_context context) {
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if(!(ctx->nucleus.provided)) {
@@ -748,7 +748,7 @@ qmckl_exit_code qmckl_provide_det_vgl_beta(qmckl_context context) {
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if(!(ctx->nucleus.provided)) {
@@ -1134,36 +1134,28 @@ end function qmckl_compute_det_vgl_beta_f
 
      #+begin_src c :tangle (eval c_test) :exports none
 
-#define walk_num chbrclf_walk_num
-#define elec_num chbrclf_elec_num
-#define shell_num chbrclf_shell_num
-#define ao_num chbrclf_ao_num
-
-int64_t elec_up_num   = chbrclf_elec_up_num;
-int64_t elec_dn_num   = chbrclf_elec_dn_num;
 double* elec_coord    = &(chbrclf_elec_coord[0][0][0]);
-const int64_t   nucl_num      = chbrclf_nucl_num;
 const double*   nucl_charge   = chbrclf_charge;
 const double*   nucl_coord    = &(chbrclf_nucl_coord[0][0]);
 
-rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
+rc = qmckl_set_electron_num (context, chbrclf_elec_up_num, chbrclf_elec_dn_num);
 assert (rc == QMCKL_SUCCESS);
 
-rc = qmckl_set_electron_walk_num (context, walk_num);
+rc = qmckl_set_electron_walk_num (context, chbrclf_walk_num);
 assert (rc == QMCKL_SUCCESS);
 
 assert(qmckl_electron_provided(context));
 
-rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
+rc = qmckl_set_electron_coord (context, 'N', elec_coord, chbrclf_walk_num*chbrclf_elec_num*3);
 assert(rc == QMCKL_SUCCESS);
 
-rc = qmckl_set_nucleus_num (context, nucl_num);
+rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num);
 assert(rc == QMCKL_SUCCESS);
 
-rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), nucl_num*3);
+rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), chbrclf_nucl_num*3);
 assert(rc == QMCKL_SUCCESS);
 
-rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
+rc = qmckl_set_nucleus_charge(context, nucl_charge, chbrclf_nucl_num);
 assert(rc == QMCKL_SUCCESS);
 
 assert(qmckl_nucleus_provided(context));
@@ -1195,27 +1187,27 @@ rc = qmckl_set_ao_basis_prim_num (context, chbrclf_prim_num);
 assert(rc == QMCKL_SUCCESS);
 assert(!qmckl_ao_basis_provided(context));
 
-rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, nucl_num);
+rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, chbrclf_nucl_num);
 assert(rc == QMCKL_SUCCESS);
 assert(!qmckl_ao_basis_provided(context));
 
-rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, nucl_num);
+rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, chbrclf_nucl_num);
 assert(rc == QMCKL_SUCCESS);
 assert(!qmckl_ao_basis_provided(context));
 
-rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom, shell_num);
+rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom, chbrclf_shell_num);
 assert(rc == QMCKL_SUCCESS);
 assert(!qmckl_ao_basis_provided(context));
 
-rc = qmckl_set_ao_basis_shell_factor  (context, shell_factor, shell_num);
+rc = qmckl_set_ao_basis_shell_factor  (context, shell_factor, chbrclf_shell_num);
 assert(rc == QMCKL_SUCCESS);
 assert(!qmckl_ao_basis_provided(context));
 
-rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, shell_num);
+rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, chbrclf_shell_num);
 assert(rc == QMCKL_SUCCESS);
 assert(!qmckl_ao_basis_provided(context));
 
-rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, shell_num);
+rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, chbrclf_shell_num);
 assert(rc == QMCKL_SUCCESS);
 assert(!qmckl_ao_basis_provided(context));
 
@@ -1239,14 +1231,13 @@ assert(rc == QMCKL_SUCCESS);
 assert(qmckl_ao_basis_provided(context));
 
 
-double ao_vgl[5][walk_num][elec_num][chbrclf_ao_num];
+double ao_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_ao_num];
 
-rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0][0]), (int64_t) 5*walk_num*elec_num*chbrclf_ao_num);
+rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]), (int64_t) 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_ao_num);
 assert (rc == QMCKL_SUCCESS);
 
 /* Set up MO data */
-const int64_t mo_num = chbrclf_mo_num;
-rc = qmckl_set_mo_basis_mo_num(context, mo_num);
+rc = qmckl_set_mo_basis_mo_num(context, chbrclf_mo_num);
 assert (rc == QMCKL_SUCCESS);
 
 const double  * mo_coefficient          =  &(chbrclf_mo_coef[0]);
@@ -1256,31 +1247,31 @@ assert (rc == QMCKL_SUCCESS);
 
 assert(qmckl_mo_basis_provided(context));
 
-double mo_vgl[5][elec_num][chbrclf_mo_num];
-rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0]));
+double mo_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_mo_num];
+rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_mo_num);
 assert (rc == QMCKL_SUCCESS);
 
 /* Set up determinant data */
 
-const int64_t det_num_alpha = 1;
-const int64_t det_num_beta  = 1;
-int64_t mo_index_alpha[det_num_alpha][walk_num][elec_up_num];
-int64_t mo_index_beta[det_num_alpha][walk_num][elec_dn_num];
+#define det_num_alpha 1
+#define det_num_beta  1
+int64_t mo_index_alpha[det_num_alpha][chbrclf_walk_num][chbrclf_elec_up_num];
+int64_t mo_index_beta[det_num_alpha][chbrclf_walk_num][chbrclf_elec_dn_num];
 
 int i, j, k;
 for(k = 0; k < det_num_alpha; ++k)
-  for(i = 0; i < walk_num; ++i)
-    for(j = 0; j < elec_up_num; ++j)
+  for(i = 0; i < chbrclf_walk_num; ++i)
+    for(j = 0; j < chbrclf_elec_up_num; ++j)
       mo_index_alpha[k][i][j] = j + 1;
 for(k = 0; k < det_num_beta; ++k)
-  for(i = 0; i < walk_num; ++i)
-    for(j = 0; j < elec_up_num; ++j)
+  for(i = 0; i < chbrclf_walk_num; ++i)
+    for(j = 0; j < chbrclf_elec_up_num; ++j)
       mo_index_beta[k][i][j] = j + 1;
 
 rc = qmckl_set_determinant_type (context, typ);
 assert(rc == QMCKL_SUCCESS);
 
-rc = qmckl_set_determinant_walk_num (context, walk_num);
+rc = qmckl_set_determinant_walk_num (context, chbrclf_walk_num);
 assert (rc == QMCKL_SUCCESS);
 
 rc = qmckl_set_determinant_det_num_alpha (context, det_num_alpha);
@@ -1297,8 +1288,8 @@ assert (rc == QMCKL_SUCCESS);
 
 // Get slater-determinant
 
-double det_vgl_alpha[det_num_alpha][walk_num][5][elec_up_num][elec_up_num];
-double det_vgl_beta[det_num_beta][walk_num][5][elec_dn_num][elec_dn_num];
+double det_vgl_alpha[det_num_alpha][chbrclf_walk_num][5][chbrclf_elec_up_num][chbrclf_elec_up_num];
+double det_vgl_beta[det_num_beta][chbrclf_walk_num][5][chbrclf_elec_dn_num][chbrclf_elec_dn_num];
 
 rc = qmckl_get_det_vgl_alpha(context, &(det_vgl_alpha[0][0][0][0][0]));
 assert (rc == QMCKL_SUCCESS);
@@ -1347,7 +1338,7 @@ qmckl_exit_code qmckl_get_det_inv_matrix_alpha(qmckl_context context, double * c
   rc = qmckl_provide_det_inv_matrix_alpha(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.up_num * ctx->electron.up_num;
@@ -1376,7 +1367,7 @@ qmckl_exit_code qmckl_get_det_inv_matrix_beta(qmckl_context context, double * co
   rc = qmckl_provide_det_inv_matrix_beta(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.down_num * ctx->electron.down_num;
@@ -1405,7 +1396,7 @@ qmckl_exit_code qmckl_get_det_adj_matrix_alpha(qmckl_context context, double * c
   rc = qmckl_provide_det_inv_matrix_alpha(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.up_num * ctx->electron.up_num;
@@ -1434,7 +1425,7 @@ qmckl_exit_code qmckl_get_det_adj_matrix_beta(qmckl_context context, double * co
   rc = qmckl_provide_det_inv_matrix_beta(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.down_num * ctx->electron.down_num;
@@ -1463,7 +1454,7 @@ qmckl_exit_code qmckl_get_det_alpha(qmckl_context context, double * const det_va
   rc = qmckl_provide_det_inv_matrix_alpha(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num;
@@ -1492,7 +1483,7 @@ qmckl_exit_code qmckl_get_det_beta(qmckl_context context, double * const det_val
   rc = qmckl_provide_det_inv_matrix_beta(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num;
@@ -1517,7 +1508,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_alpha(qmckl_context context) {
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if(!(ctx->nucleus.provided)) {
@@ -1640,7 +1631,7 @@ qmckl_exit_code qmckl_provide_det_inv_matrix_beta(qmckl_context context) {
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if(!(ctx->nucleus.provided)) {
@@ -2047,8 +2038,8 @@ end function qmckl_compute_det_inv_matrix_beta_f
      #+begin_src c :tangle (eval c_test) :exports none
 // Get adjoint of the slater-determinant
 
-double det_inv_matrix_alpha[det_num_alpha][walk_num][elec_up_num][elec_up_num];
-double det_inv_matrix_beta[det_num_beta][walk_num][elec_dn_num][elec_dn_num];
+double det_inv_matrix_alpha[det_num_alpha][chbrclf_walk_num][chbrclf_elec_up_num][chbrclf_elec_up_num];
+double det_inv_matrix_beta[det_num_beta][chbrclf_walk_num][chbrclf_elec_dn_num][chbrclf_elec_dn_num];
 
 rc = qmckl_get_det_inv_matrix_alpha(context, &(det_inv_matrix_alpha[0][0][0][0]));
 assert (rc == QMCKL_SUCCESS);
diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org
index 33ac366..c0c0254 100644
--- a/org/qmckl_electron.org
+++ b/org/qmckl_electron.org
@@ -97,8 +97,8 @@ int main() {
   | ~ee_distance_rescaled_deriv_e_date~ | ~uint64_t~                           | Last modification date of the electron-electron distance derivatives |
   | ~ee_pot~                            | ~double[walk_num]~                   | Electron-electron rescaled distances derivatives                     |
   | ~ee_pot_date~                       | ~uint64_t~                           | Last modification date of the electron-electron distance derivatives |
-  | ~en_pot~                            | double[walk_num]                     | Electron-nucleus potential energy                                    |
-  | ~en_pot_date~                       | int64_t                              | Date when the electron-nucleus potential energy was computed         |
+  | ~en_pot~                            | ~double[walk_num]~                   | Electron-nucleus potential energy                                    |
+  | ~en_pot_date~                       | ~int64_t~                            | Date when the electron-nucleus potential energy was computed         |
   | ~en_distance_rescaled~              | ~double[walk_num][nucl_num][num]~    | Electron-nucleus distances                                           |
   | ~en_distance_rescaled_date~         | ~uint64_t~                           | Last modification date of the electron-electron distances            |
   | ~en_distance_rescaled_deriv_e~      | ~double[walk_num][4][nucl_num][num]~ | Electron-electron rescaled distances derivatives                     |
@@ -157,7 +157,7 @@ qmckl_exit_code qmckl_init_electron(qmckl_context context) {
     return false;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   ctx->electron.uninitialized = (1 << 2) - 1;
@@ -182,7 +182,7 @@ bool qmckl_electron_provided(const qmckl_context context) {
     return false;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   return ctx->electron.provided;
@@ -228,7 +228,7 @@ qmckl_get_electron_num (const qmckl_context context, int64_t* const num) {
                            "num is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 0;
@@ -256,7 +256,7 @@ qmckl_get_electron_up_num (const qmckl_context context, int64_t* const up_num) {
                            "up_num is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 0;
@@ -284,7 +284,7 @@ qmckl_get_electron_down_num (const qmckl_context context, int64_t* const down_nu
                            "down_num is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 0;
@@ -323,7 +323,7 @@ qmckl_get_electron_walk_num (const qmckl_context context, int64_t* const walk_nu
                            "walk_num is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 1;
@@ -360,7 +360,7 @@ qmckl_get_electron_rescale_factor_ee (const qmckl_context context, double* const
                            "rescale_factor_kappa_ee is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   assert (ctx->electron.rescale_factor_kappa_ee > 0.0);
@@ -383,7 +383,7 @@ qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const
                            "rescale_factor_kappa_en is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   assert (ctx->electron.rescale_factor_kappa_en > 0.0);
@@ -448,7 +448,7 @@ qmckl_get_electron_coord (const qmckl_context context,
     return QMCKL_INVALID_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if (!ctx->electron.provided) {
@@ -489,7 +489,7 @@ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
   return QMCKL_NULL_CONTEXT;
  }
 
-qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
    #+end_src
 
    #+NAME:post2
@@ -897,7 +897,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance(qmckl_context context, double* co
   rc = qmckl_provide_ee_distance(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num;
@@ -921,7 +921,7 @@ qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
 
@@ -1138,7 +1138,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance_rescaled(qmckl_context context, d
   rc = qmckl_provide_ee_distance_rescaled(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num;
@@ -1162,7 +1162,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
 
@@ -1218,7 +1218,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context)
     | ~elec_num~                | ~int64_t~                              | in     | Number of electrons                  |
     | ~rescale_factor_kappa_ee~ | ~double~                               | in     | Factor to rescale ee distances       |
     | ~walk_num~                | ~int64_t~                              | in     | Number of walkers                    |
-    | ~coord~                   | ~double[walk_num][3][elec_num]~        | in     | Electron coordinates                 |
+    | ~coord~                   | ~double[3][walk_num][elec_num]~        | in     | Electron coordinates                 |
     | ~ee_distance~             | ~double[walk_num][elec_num][elec_num]~ | out    | Electron-electron rescaled distances |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
@@ -1231,7 +1231,7 @@ integer function qmckl_compute_ee_distance_rescaled_f(context, elec_num, rescale
   integer*8             , intent(in)  :: elec_num
   double precision      , intent(in)  :: rescale_factor_kappa_ee
   integer*8             , intent(in)  :: walk_num
-  double precision      , intent(in)  :: coord(elec_num,3,walk_num)
+  double precision      , intent(in)  :: coord(elec_num,walk_num,3)
   double precision      , intent(out) :: ee_distance_rescaled(elec_num,elec_num,walk_num)
 
   integer*8 :: k
@@ -1357,7 +1357,7 @@ assert(fabs(ee_distance_rescaled[elec_num*elec_num+1]-0.9985724058042633) < 1.e-
 
      #+end_src
 
-** Electron-electron rescaled distance gradients and laplacian with respect to electron coords
+** Electron-electron rescaled distance gradients and Laplacian with respect to electron coords
 
   The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$
   needs to be perturbed with respect to the electorn coordinates.
@@ -1384,7 +1384,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance_rescaled_deriv_e(qmckl_context co
   rc = qmckl_provide_ee_distance_rescaled_deriv_e(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = 4 * ctx->electron.num * ctx->electron.num * ctx->electron.walk_num;
@@ -1408,7 +1408,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
 
@@ -1464,7 +1464,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context
     | ~elec_num~                | ~int64_t~                                 | in     | Number of electrons                             |
     | ~rescale_factor_kappa_ee~ | ~double~                                  | in     | Factor to rescale ee distances                  |
     | ~walk_num~                | ~int64_t~                                 | in     | Number of walkers                               |
-    | ~coord~                   | ~double[walk_num][3][elec_num]~           | in     | Electron coordinates                            |
+    | ~coord~                   | ~double[3][walk_num][elec_num]~           | in     | Electron coordinates                            |
     | ~ee_distance_deriv_e~     | ~double[walk_num][4][elec_num][elec_num]~ | out    | Electron-electron rescaled distance derivatives |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
@@ -1477,7 +1477,7 @@ integer function qmckl_compute_ee_distance_rescaled_deriv_e_f(context, elec_num,
   integer*8             , intent(in)  :: elec_num
   double precision      , intent(in)  :: rescale_factor_kappa_ee
   integer*8             , intent(in)  :: walk_num
-  double precision      , intent(in)  :: coord(elec_num,3,walk_num)
+  double precision      , intent(in)  :: coord(elec_num,walk_num,3)
   double precision      , intent(out) :: ee_distance_rescaled_deriv_e(4,elec_num,elec_num,walk_num)
 
   integer*8 :: k
@@ -1501,8 +1501,8 @@ integer function qmckl_compute_ee_distance_rescaled_deriv_e_f(context, elec_num,
 
   do k=1,walk_num
      info = qmckl_distance_rescaled_deriv_e(context, 'T', 'T', elec_num, elec_num, &
-          coord(1,1,k), elec_num, &
-          coord(1,1,k), elec_num, &
+          coord(1,k,1), elec_num*walk_num, &
+          coord(1,k,1), elec_num*walk_num, &
           ee_distance_rescaled_deriv_e(1,1,1,k), elec_num, rescale_factor_kappa_ee)
      if (info /= QMCKL_SUCCESS) then
         exit
@@ -1613,7 +1613,7 @@ qmckl_exit_code qmckl_get_electron_ee_potential(qmckl_context context, double* c
   rc = qmckl_provide_ee_potential(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->electron.walk_num * sizeof(double);
@@ -1637,7 +1637,7 @@ qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if (!ctx->electron.provided) return QMCKL_NOT_PROVIDED;
@@ -1818,7 +1818,7 @@ qmckl_exit_code qmckl_get_electron_en_distance(qmckl_context context, double* di
   rc = qmckl_provide_en_distance(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num;
@@ -1842,7 +1842,7 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if (!(ctx->nucleus.provided)) {
@@ -1905,7 +1905,7 @@ qmckl_exit_code qmckl_provide_en_distance(qmckl_context context)
     | ~elec_num~    | ~int64_t~                              | in     | Number of electrons        |
     | ~nucl_num~    | ~int64_t~                              | in     | Number of nuclei           |
     | ~walk_num~    | ~int64_t~                              | in     | Number of walkers          |
-    | ~elec_coord~  | ~double[walk_num][3][elec_num]~        | in     | Electron coordinates       |
+    | ~elec_coord~  | ~double[3][walk_num][elec_num]~        | in     | Electron coordinates       |
     | ~nucl_coord~  | ~double[3][elec_num]~                  | in     | Nuclear coordinates        |
     | ~en_distance~ | ~double[walk_num][nucl_num][elec_num]~ | out    | Electron-nucleus distances |
 
@@ -2097,7 +2097,7 @@ qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, d
   rc = qmckl_provide_en_distance_rescaled(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num;
@@ -2121,7 +2121,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if (!(ctx->nucleus.provided)) {
@@ -2183,7 +2183,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled(qmckl_context context)
     | ~nucl_num~                | ~int64_t~                              | in     | Number of nuclei                  |
     | ~rescale_factor_kappa_en~ | ~double~                               | in     | The factor for rescaled distances |
     | ~walk_num~                | ~int64_t~                              | in     | Number of walkers                 |
-    | ~elec_coord~              | ~double[walk_num][3][elec_num]~        | in     | Electron coordinates              |
+    | ~elec_coord~              | ~double[3][walk_num][elec_num]~        | in     | Electron coordinates              |
     | ~nucl_coord~              | ~double[3][elec_num]~                  | in     | Nuclear coordinates               |
     | ~en_distance_rescaled~    | ~double[walk_num][nucl_num][elec_num]~ | out    | Electron-nucleus distances        |
 
@@ -2385,7 +2385,7 @@ qmckl_exit_code qmckl_get_electron_en_distance_rescaled_deriv_e(qmckl_context co
   rc = qmckl_provide_en_distance_rescaled_deriv_e(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num;
@@ -2409,7 +2409,7 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if (!(ctx->nucleus.provided)) {
@@ -2471,9 +2471,9 @@ qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context
     | ~nucl_num~                     | ~int64_t~                                 | in     | Number of nuclei                      |
     | ~rescale_factor_kappa_en~      | ~double~                                  | in     | The factor for rescaled distances     |
     | ~walk_num~                     | ~int64_t~                                 | in     | Number of walkers                     |
-    | ~elec_coord~                   | ~double[walk_num][3][elec_num]~           | in     | Electron coordinates                  |
+    | ~elec_coord~                   | ~double[3][walk_num][elec_num]~           | in     | Electron coordinates                  |
     | ~nucl_coord~                   | ~double[3][elec_num]~                     | in     | Nuclear coordinates                   |
-    | ~en_distance_rescaled_deriv_e~ | ~double[walk_num][4][nucl_num][elec_num]~ | out    | Electron-nucleus distance derivatives |
+    | ~en_distance_rescaled_deriv_e~ | ~double[walk_num][nucl_num][elec_num][4]~ | out    | Electron-nucleus distance derivatives |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
 integer function qmckl_compute_en_distance_rescaled_deriv_e_f(context, elec_num, nucl_num, &
@@ -2656,7 +2656,7 @@ qmckl_exit_code qmckl_get_electron_en_potential(qmckl_context context, double* c
   rc = qmckl_provide_en_potential(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->electron.walk_num * sizeof(double);
@@ -2680,7 +2680,7 @@ qmckl_exit_code qmckl_provide_en_potential(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if (!ctx->electron.provided) return QMCKL_NOT_PROVIDED;
@@ -2843,7 +2843,7 @@ assert (rc == QMCKL_SUCCESS);
 
 *** Compute                                                        :noexport:
 
-    # begin_src f90 :comments org :tangle (eval f) :noweb yes
+   # begin_src f90 :comments org :tangle (eval f) :noweb yes
 subroutine draw_init_points
   implicit none
   BEGIN_DOC
diff --git a/org/qmckl_error.org b/org/qmckl_error.org
index 97cc7cc..d23fa8e 100644
--- a/org/qmckl_error.org
+++ b/org/qmckl_error.org
@@ -239,7 +239,7 @@ for (text, code, message) in table:
     message = message.replace("'",'"')
     result += [ f"""case {text}:
       return {message};
-    break;""" ]
+    """ ]
 return '\n'.join(result)
 
    #+end_src
@@ -247,89 +247,91 @@ return '\n'.join(result)
    #+RESULTS: cases
    #+begin_example
    case QMCKL_SUCCESS:
-		 return "Success";
-	       break;
+         return "Success";
+
    case QMCKL_INVALID_ARG_1:
-		 return "Invalid argument 1";
-	       break;
+         return "Invalid argument 1";
+
    case QMCKL_INVALID_ARG_2:
-		 return "Invalid argument 2";
-	       break;
+         return "Invalid argument 2";
+
    case QMCKL_INVALID_ARG_3:
-		 return "Invalid argument 3";
-	       break;
+         return "Invalid argument 3";
+
    case QMCKL_INVALID_ARG_4:
-		 return "Invalid argument 4";
-	       break;
+         return "Invalid argument 4";
+
    case QMCKL_INVALID_ARG_5:
-		 return "Invalid argument 5";
-	       break;
+         return "Invalid argument 5";
+
    case QMCKL_INVALID_ARG_6:
-		 return "Invalid argument 6";
-	       break;
+         return "Invalid argument 6";
+
    case QMCKL_INVALID_ARG_7:
-		 return "Invalid argument 7";
-	       break;
+         return "Invalid argument 7";
+
    case QMCKL_INVALID_ARG_8:
-		 return "Invalid argument 8";
-	       break;
+         return "Invalid argument 8";
+
    case QMCKL_INVALID_ARG_9:
-		 return "Invalid argument 9";
-	       break;
+         return "Invalid argument 9";
+
    case QMCKL_INVALID_ARG_10:
-		 return "Invalid argument 10";
-	       break;
+         return "Invalid argument 10";
+
    case QMCKL_INVALID_ARG_11:
-		 return "Invalid argument 11";
-	       break;
+         return "Invalid argument 11";
+
    case QMCKL_INVALID_ARG_12:
-		 return "Invalid argument 12";
-	       break;
+         return "Invalid argument 12";
+
    case QMCKL_INVALID_ARG_13:
-		 return "Invalid argument 13";
-	       break;
+         return "Invalid argument 13";
+
    case QMCKL_INVALID_ARG_14:
-		 return "Invalid argument 14";
-	       break;
+         return "Invalid argument 14";
+
    case QMCKL_INVALID_ARG_15:
-		 return "Invalid argument 15";
-	       break;
+         return "Invalid argument 15";
+
    case QMCKL_INVALID_ARG_16:
-		 return "Invalid argument 16";
-	       break;
+         return "Invalid argument 16";
+
    case QMCKL_INVALID_ARG_17:
-		 return "Invalid argument 17";
-	       break;
+         return "Invalid argument 17";
+
    case QMCKL_INVALID_ARG_18:
-		 return "Invalid argument 18";
-	       break;
+         return "Invalid argument 18";
+
    case QMCKL_INVALID_ARG_19:
-		 return "Invalid argument 19";
-	       break;
+         return "Invalid argument 19";
+
    case QMCKL_INVALID_ARG_20:
-		 return "Invalid argument 20";
-	       break;
+         return "Invalid argument 20";
+
    case QMCKL_FAILURE:
-		 return "Failure";
-	       break;
+         return "Failure";
+
    case QMCKL_ERRNO:
-		 return strerror(errno);
-	       break;
+         return strerror(errno);
+
    case QMCKL_INVALID_CONTEXT:
-		 return "Invalid context";
-	       break;
+         return "Invalid context";
+
    case QMCKL_ALLOCATION_FAILED:
-		 return "Allocation failed";
-	       break;
+         return "Allocation failed";
+
    case QMCKL_DEALLOCATION_FAILED:
-		 return "De-allocation failed";
-	       break;
+         return "De-allocation failed";
+
    case QMCKL_NOT_PROVIDED:
-		 return "Not provided";
-	       break;
+         return "Not provided";
+
+   case QMCKL_OUT_OF_BOUNDS:
+         return "Index out of bounds";
+
    case QMCKL_INVALID_EXIT_CODE:
-		 return "Invalid exit code";
-	       break;
+         return "Invalid exit code";
    #+end_example
 
    # Source
@@ -414,7 +416,7 @@ qmckl_set_error(qmckl_context context,
 
   qmckl_lock(context);
   {
-    qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+    qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
     assert (ctx != NULL); /* Impossible because the context is valid. */
 
     ctx->error.exit_code = exit_code;
@@ -460,7 +462,7 @@ qmckl_get_error(qmckl_context context,
 
   qmckl_lock(context);
   {
-    qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+    qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
     assert (ctx != NULL); /* Impossible because the context is valid. */
 
 /* Turn off annoying GCC warning */
diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org
index 8736c0b..de7048e 100644
--- a/org/qmckl_jastrow.org
+++ b/org/qmckl_jastrow.org
@@ -155,25 +155,25 @@ int main() {
 
   computed data:
 
-    | Variable                   | Type                                                                | In/Out                                          | Description                     |
-    |----------------------------+---------------------------------------------------------------------+-------------------------------------------------+---------------------------------|
-    | ~dim_cord_vect~            | ~int64_t~                                                           | Number of unique C coefficients                 |                                 |
-    | ~dim_cord_vect_date~       | ~uint64_t~                                                          | Number of unique C coefficients                 |                                 |
-    | ~asymp_jasb~               | ~double[2]~                                                         | Asymptotic component                            |                                 |
-    | ~asymp_jasb_date~          | ~uint64_t~                                                          | Asymptotic component                            |                                 |
-    | ~cord_vect_full~           | ~double[dim_cord_vect][nucl_num]~                                   | vector of non-zero coefficients                 |                                 |
-    | ~cord_vect_full_date~      | ~uint64_t~                                                          | Keep track of changes here                      |                                 |
-    | ~lkpm_combined_index~      | ~int64_t[4][dim_cord_vect]~                                         | Transform l,k,p, and m into consecutive indices |                                 |
-    | ~lkpm_combined_index_date~ | ~uint64_t~                                                          | Transform l,k,p, and m into consecutive indices |                                 |
-    | ~tmp_c~                    | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~    | out                                             | vector of non-zero coefficients |
-    | ~dtmp_c~                   | ~double[walk_num][elec_num][4][nucl_num][0:cord_num][0:cord_num-1]~ | vector of non-zero coefficients                 |                                 |
+    | Variable                   | Type                                                            | In/Out                                          |
+    |----------------------------+-----------------------------------------------------------------+-------------------------------------------------|
+    | ~dim_cord_vect~            | ~int64_t~                                                       | Number of unique C coefficients                 |
+    | ~dim_cord_vect_date~       | ~uint64_t~                                                      | Number of unique C coefficients                 |
+    | ~asymp_jasb~               | ~double[2]~                                                     | Asymptotic component                            |
+    | ~asymp_jasb_date~          | ~uint64_t~                                                      | Asymptotic component                            |
+    | ~cord_vect_full~           | ~double[dim_cord_vect][nucl_num]~                               | vector of non-zero coefficients                 |
+    | ~cord_vect_full_date~      | ~uint64_t~                                                      | Keep track of changes here                      |
+    | ~lkpm_combined_index~      | ~int64_t[4][dim_cord_vect]~                                     | Transform l,k,p, and m into consecutive indices |
+    | ~lkpm_combined_index_date~ | ~uint64_t~                                                      | Transform l,k,p, and m into consecutive indices |
+    | ~tmp_c~                    | ~double[walk_num][cord_num][cord_num+1][nucl_num][elec_num]~    | vector of non-zero coefficients                 |
+    | ~dtmp_c~                   | ~double[walk_num][elec_num][4][nucl_num][cord_num+1][cord_num]~ | vector of non-zero coefficients                 |
 
-    | ~een_rescaled_n~              | ~double[walk_num][elec_num][nucl_num][0:cord_num]~    | The electron-electron rescaled distances raised to the powers defined by cord                           |   |
-    | ~een_rescaled_n_date~         | ~uint64_t~                                            | Keep track of the date of creation                                                                      |   |
-    | ~een_rescaled_e_deriv_e~      | ~double[walk_num][elec_num][4][elec_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons |   |
-    | ~een_rescaled_e_deriv_e_date~ | ~uint64_t~                                            | Keep track of the date of creation                                                                      |   |
-    | ~een_rescaled_n_deriv_e~      | ~double[walk_num][elec_num][4][nucl_num][0:cord_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons |   |
-    | ~een_rescaled_n_deriv_e_date~ | ~uint64_t~                                            | Keep track of the date of creation                                                                      |   |
+    | ~een_rescaled_n~              | ~double[walk_num][cord_num+1][nucl_num][elec_num]~  | The electron-electron rescaled distances raised to the powers defined by cord                           |
+    | ~een_rescaled_n_date~         | ~uint64_t~                                            | Keep track of the date of creation                                                                      |
+    | ~een_rescaled_e_deriv_e~      | ~double[walk_num][cord_num+1][elec_num][4][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons |
+    | ~een_rescaled_e_deriv_e_date~ | ~uint64_t~                                            | Keep track of the date of creation                                                                      |
+    | ~een_rescaled_n_deriv_e~      | ~double[walk_num][cord_num+1][nucl_num][4][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons |
+    | ~een_rescaled_n_deriv_e_date~ | ~uint64_t~                                            | Keep track of the date of creation                                                                      |
 
    #+NAME: jastrow_data
    #+BEGIN_SRC python :results none :exports none
@@ -404,7 +404,7 @@ qmckl_exit_code qmckl_init_jastrow(qmckl_context context) {
     return false;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   ctx->jastrow.uninitialized = (1 << 5) - 1;
@@ -447,7 +447,7 @@ bool qmckl_jastrow_provided(const qmckl_context context) {
     return false;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   return ctx->jastrow.provided;
@@ -475,7 +475,7 @@ qmckl_exit_code qmckl_get_jastrow_aord_num (const qmckl_context context, int64_t
                            "aord_num is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 0;
@@ -502,7 +502,7 @@ qmckl_exit_code qmckl_get_jastrow_bord_num (const qmckl_context context, int64_t
                            "aord_num is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 0;
@@ -529,7 +529,7 @@ qmckl_exit_code qmckl_get_jastrow_cord_num (const qmckl_context context, int64_t
                            "aord_num is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 0;
@@ -556,7 +556,7 @@ qmckl_exit_code qmckl_get_jastrow_type_nucl_num (const qmckl_context context, in
                            "type_nucl_num is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 1;
@@ -587,7 +587,7 @@ qmckl_get_jastrow_type_nucl_vector (const qmckl_context context,
                            "type_nucl_vector is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 2;
@@ -624,7 +624,7 @@ qmckl_get_jastrow_aord_vector (const qmckl_context context,
                            "aord_vector is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 3;
@@ -661,7 +661,7 @@ qmckl_get_jastrow_bord_vector (const qmckl_context context,
                            "bord_vector is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 4;
@@ -698,7 +698,7 @@ qmckl_get_jastrow_cord_vector (const qmckl_context context,
                            "cord_vector is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 5;
@@ -773,7 +773,7 @@ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
   return QMCKL_NULL_CONTEXT;
  }
 
-qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
    #+end_src
 
    #+NAME:post2
@@ -1126,7 +1126,7 @@ qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context) {
     return QMCKL_INVALID_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   /* ----------------------------------- */
@@ -1376,7 +1376,7 @@ qmckl_get_jastrow_asymp_jasb(qmckl_context context,
   rc = qmckl_provide_asymp_jasb(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int64_t sze = 2;
@@ -1407,7 +1407,7 @@ qmckl_exit_code qmckl_provide_asymp_jasb(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   /* Check if ee kappa is provided */
@@ -1661,7 +1661,7 @@ qmckl_get_jastrow_factor_ee(qmckl_context context,
   rc = qmckl_provide_factor_ee(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int64_t sze=ctx->electron.walk_num;
@@ -1692,7 +1692,7 @@ qmckl_exit_code qmckl_provide_factor_ee(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   /* Check if ee rescaled distance is provided */
@@ -1773,7 +1773,7 @@ integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num,
   double precision      , intent(out) :: factor_ee(walk_num)
 
   integer*8 :: i, j, p, ipar, nw
-  double precision   :: pow_ser, x, spin_fact, power_ser
+  double precision   :: x, power_ser, spin_fact
 
   info = QMCKL_SUCCESS
 
@@ -1812,7 +1812,7 @@ integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num,
           power_ser = power_ser + bord_vector(p + 1) * x
         end do
 
-        if(j .LE. up_num .OR. i .GT. up_num) then
+        if(j <= up_num .OR. i > up_num) then
           spin_fact = 0.5d0
           ipar = 2
         endif
@@ -1831,7 +1831,7 @@ end function qmckl_compute_factor_ee_f
     #+end_src
 
 #+begin_src c   :comments org :tangle (eval c) :noweb yes
-    qmckl_exit_code qmckl_compute_factor_ee (
+qmckl_exit_code qmckl_compute_factor_ee (
 	  const qmckl_context context,
 	  const int64_t walk_num,
 	  const int64_t elec_num,
@@ -1843,7 +1843,7 @@ end function qmckl_compute_factor_ee_f
 	  double* const factor_ee ) {
 
   int ipar; // can we use a smaller integer?
-  double pow_ser, x, x1, spin_fact, power_ser;
+  double x, x1, spin_fact, power_ser;
 
   if (context == QMCKL_NULL_CONTEXT) { 
      return QMCKL_INVALID_CONTEXT;
@@ -1998,7 +1998,7 @@ qmckl_get_jastrow_factor_ee_deriv_e(qmckl_context context,
   rc = qmckl_provide_factor_ee_deriv_e(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num;
@@ -2030,7 +2030,7 @@ qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   /* Check if ee rescaled distance is provided */
@@ -2104,9 +2104,10 @@ qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context)
     | ~factor_ee_deriv_e~            | ~double[walk_num][4][elec_num]~           | out    | Electron-electron distances |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num, up_num, bord_num,            &
-                                           bord_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e,  &
-					   asymp_jasb, factor_ee_deriv_e) &
+integer function qmckl_compute_factor_ee_deriv_e_f( &
+     context, walk_num, elec_num, up_num, bord_num, &
+     bord_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e,  &
+     asymp_jasb, factor_ee_deriv_e) &
      result(info)
   use qmckl
   implicit none
@@ -2114,7 +2115,7 @@ integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num,
   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(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)  :: ee_distance_rescaled_deriv_e(4,elec_num, elec_num,walk_num)   !TODO
   double precision      , intent(in)  :: asymp_jasb(2)
   double precision      , intent(out) :: factor_ee_deriv_e(elec_num,4,walk_num)
 
@@ -2229,7 +2230,7 @@ end function qmckl_compute_factor_ee_deriv_e_f
 
     #+RESULTS:
     #+begin_src f90 :tangle (eval f) :comments org :exports none
-    integer(c_int32_t) function qmckl_compute_factor_ee_deriv_e &
+integer(c_int32_t) function qmckl_compute_factor_ee_deriv_e &
 	(context, &
 		 walk_num, &
 		 elec_num, &
@@ -2425,7 +2426,7 @@ qmckl_get_jastrow_factor_en(qmckl_context context,
   rc = qmckl_provide_factor_en(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int64_t sze=ctx->electron.walk_num;
@@ -2456,7 +2457,7 @@ qmckl_exit_code qmckl_provide_factor_en(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   /* Check if en rescaled distance is provided */
@@ -2526,9 +2527,10 @@ qmckl_exit_code qmckl_provide_factor_en(qmckl_context context)
     | ~factor_en~            | ~double[walk_num]~                     | out    | Electron-nucleus jastrow   |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_factor_en_f(context, walk_num, elec_num, nucl_num, type_nucl_num,            &
-                                           type_nucl_vector, aord_num, aord_vector,                         &
-                                           en_distance_rescaled, factor_en) &
+integer function qmckl_compute_factor_en_f( &
+     context, walk_num, elec_num, nucl_num, type_nucl_num, &
+     type_nucl_vector, aord_num, aord_vector, &
+     en_distance_rescaled, factor_en) &
      result(info)
   use qmckl
   implicit none
@@ -2539,8 +2541,8 @@ integer function qmckl_compute_factor_en_f(context, walk_num, elec_num, nucl_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
-  double precision   :: x, spin_fact, power_ser
+  integer*8 :: i, a, p, nw
+  double precision   :: x, power_ser
 
   info = QMCKL_SUCCESS
 
@@ -2609,9 +2611,7 @@ qmckl_exit_code qmckl_compute_factor_en (
 	  const double* en_distance_rescaled,
 	  double* const factor_en ) {
 
-
-  int     ipar;
-  double  x, x1, spin_fact, power_ser;
+  double  x, x1, power_ser;
 
 
   if (context == QMCKL_NULL_CONTEXT) {
@@ -2630,10 +2630,30 @@ qmckl_exit_code qmckl_compute_factor_en (
      return QMCKL_INVALID_ARG_4;
   }
 
+  if (type_nucl_num <= 0) {
+    return QMCKL_INVALID_ARG_5;
+  }
+
+  if (type_nucl_vector == NULL) {
+    return QMCKL_INVALID_ARG_6;
+  }
+
   if (aord_num <= 0) {
      return QMCKL_INVALID_ARG_7;
   }
 
+  if (aord_vector == NULL) {
+    return QMCKL_INVALID_ARG_8;
+  }
+
+  if (en_distance_rescaled == NULL) {
+    return QMCKL_INVALID_ARG_9;
+  }
+
+  if (factor_en == NULL) {
+    return QMCKL_INVALID_ARG_10;
+  }
+
 
   for (int nw = 0; nw < walk_num; ++nw ) {
     // init array
@@ -2750,7 +2770,7 @@ qmckl_get_jastrow_factor_en_deriv_e(qmckl_context context,
   rc = qmckl_provide_factor_en_deriv_e(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num;
@@ -2781,7 +2801,7 @@ qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   /* Check if en rescaled distance is provided */
@@ -2857,9 +2877,10 @@ qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context)
     | ~factor_en_deriv_e~            | ~double[walk_num][4][elec_num]~           | out    | Electron-nucleus jastrow              |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_factor_en_deriv_e_f(context, walk_num, elec_num, nucl_num, type_nucl_num,            &
-                                           type_nucl_vector, aord_num, aord_vector,                         &
-                                           en_distance_rescaled, en_distance_rescaled_deriv_e, factor_en_deriv_e) &
+integer function qmckl_compute_factor_en_deriv_e_f( &
+     context, walk_num, elec_num, nucl_num, type_nucl_num, &
+     type_nucl_vector, aord_num, aord_vector, &
+     en_distance_rescaled, en_distance_rescaled_deriv_e, factor_en_deriv_e) &
      result(info)
   use qmckl
   implicit none
@@ -2872,7 +2893,7 @@ integer function qmckl_compute_factor_en_deriv_e_f(context, walk_num, elec_num,
   double precision      , intent(out) :: factor_en_deriv_e(elec_num,4,walk_num)
 
   integer*8 :: i, a, p, ipar, nw, ii
-  double precision   :: x, spin_fact, den, invden, invden2, invden3, xinv
+  double precision   :: x, den, invden, invden2, invden3, xinv
   double precision   :: y, lap1, lap2, lap3, third
   double precision, dimension(3) :: power_ser_g
   double precision, dimension(4) :: dx
@@ -2939,7 +2960,7 @@ integer function qmckl_compute_factor_en_deriv_e_f(context, walk_num, elec_num,
           lap3 = lap3 - 2.0d0 * aord_vector(2, type_nucl_vector(a)) * dx(ii) * dx(ii)
 
           factor_en_deriv_e(i, ii, nw) = factor_en_deriv_e(i, ii, nw) + aord_vector(1, type_nucl_vector(a)) &
-                                  * dx(ii) * invden2                                                        &
+                                  ,* dx(ii) * invden2                                                        &
                                   + power_ser_g(ii)
 
         end do
@@ -3172,7 +3193,7 @@ qmckl_get_jastrow_een_rescaled_e(qmckl_context context,
   rc = qmckl_provide_een_rescaled_e(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int64_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
@@ -3202,7 +3223,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   /* Check if ee distance is provided */
@@ -3267,7 +3288,8 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context)
     | ~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,  &
+integer function qmckl_compute_een_rescaled_e_f( &
+     context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee,  &
      ee_distance, een_rescaled_e) &
      result(info)
   use qmckl
@@ -3372,7 +3394,8 @@ end function qmckl_compute_een_rescaled_e_f
     #+RESULTS:
     #+begin_src f90 :tangle (eval f) :comments org :exports none
     integer(c_int32_t) function qmckl_compute_een_rescaled_e &
-	(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e) &
+	(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, &
+        ee_distance, een_rescaled_e) &
 	bind(C) result(info)
 
       use, intrinsic :: iso_c_binding
@@ -3504,7 +3527,7 @@ qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context,
   rc = qmckl_provide_een_rescaled_e_deriv_e(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int64_t sze = ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
@@ -3534,7 +3557,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   /* Check if ee distance is provided */
@@ -3603,7 +3626,8 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context)
     | ~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,  &
+integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f( &
+     context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee,  &
      coord_new, ee_distance, een_rescaled_e, een_rescaled_e_deriv_e) &
      result(info)
   use qmckl
@@ -3883,7 +3907,7 @@ qmckl_get_jastrow_een_rescaled_n(qmckl_context context,
   rc = qmckl_provide_een_rescaled_n(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int64_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
@@ -3913,7 +3937,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   /* Check if ee distance is provided */
@@ -3980,7 +4004,8 @@ qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context)
     | ~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,  &
+integer function qmckl_compute_een_rescaled_n_f( &
+     context, walk_num, elec_num, nucl_num, cord_num, rescale_factor_kappa_en,  &
      en_distance, een_rescaled_n) &
      result(info)
   use qmckl
@@ -4217,7 +4242,7 @@ qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context,
   rc = qmckl_provide_een_rescaled_n_deriv_e(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int64_t sze = ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
@@ -4247,7 +4272,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   /* Check if ee distance is provided */
@@ -4324,9 +4349,10 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context)
     | ~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, &
-     cord_num, rescale_factor_kappa_en,                                                               &
-     coord_new, coord, en_distance, een_rescaled_n, een_rescaled_n_deriv_e)                           &
+integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f( &
+     context, walk_num, elec_num, nucl_num, &
+     cord_num, rescale_factor_kappa_en, &
+     coord_new, coord, en_distance, een_rescaled_n, een_rescaled_n_deriv_e) &
      result(info)
   use qmckl
   implicit none
@@ -4599,7 +4625,7 @@ qmckl_exit_code qmckl_get_jastrow_dim_cord_vect(qmckl_context context, int64_t*
   rc = qmckl_provide_dim_cord_vect(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   *dim_cord_vect = ctx->jastrow.dim_cord_vect;
@@ -4621,7 +4647,7 @@ qmckl_exit_code qmckl_get_jastrow_cord_vect_full(qmckl_context context, double*
   rc = qmckl_provide_cord_vect_full(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->jastrow.dim_cord_vect * ctx->nucleus.num;
@@ -4644,7 +4670,7 @@ qmckl_exit_code qmckl_get_jastrow_lkpm_combined_index(qmckl_context context, int
   rc = qmckl_provide_cord_vect_full(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->jastrow.dim_cord_vect * 4;
@@ -4670,7 +4696,7 @@ qmckl_exit_code qmckl_get_jastrow_tmp_c(qmckl_context context, double* const tmp
   rc = qmckl_provide_tmp_c(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = (ctx->jastrow.cord_num) * (ctx->jastrow.cord_num + 1)
@@ -4697,7 +4723,7 @@ qmckl_exit_code qmckl_get_jastrow_dtmp_c(qmckl_context context, double* const dt
   rc = qmckl_provide_dtmp_c(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = (ctx->jastrow.cord_num) * (ctx->jastrow.cord_num + 1)
@@ -4726,7 +4752,7 @@ qmckl_exit_code qmckl_provide_dim_cord_vect(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   /* Compute if necessary */
@@ -4753,7 +4779,7 @@ qmckl_exit_code qmckl_provide_cord_vect_full(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   /* Check if dim_cord_vect is provided */
@@ -4804,7 +4830,7 @@ qmckl_exit_code qmckl_provide_lkpm_combined_index(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   /* Check if dim_cord_vect is provided */
@@ -4851,7 +4877,7 @@ qmckl_exit_code qmckl_provide_tmp_c(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   /* Check if dim_cord_vect is provided */
@@ -4944,7 +4970,7 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   /* Check if dim_cord_vect is provided */
@@ -5047,7 +5073,8 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context)
     | ~dim_cord_vect~ | ~int64_t~       | out    | dimension of cord_vect_full table |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_dim_cord_vect_f(context, cord_num, dim_cord_vect) &
+integer function qmckl_compute_dim_cord_vect_f( &
+     context, cord_num, dim_cord_vect) &
      result(info)
   use qmckl
   implicit none
@@ -5155,7 +5182,8 @@ qmckl_exit_code qmckl_compute_dim_cord_vect (
     | ~cord_vect_full~   | ~double[dim_cord_vect][nucl_num]~      | out    | Full list of coefficients    |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_cord_vect_full_f(context, nucl_num, dim_cord_vect, type_nucl_num,  &
+integer function qmckl_compute_cord_vect_full_f( &
+     context, nucl_num, dim_cord_vect, type_nucl_num,  &
      type_nucl_vector, cord_vector, cord_vect_full) &
      result(info)
   use qmckl
@@ -5257,8 +5285,8 @@ end function qmckl_compute_cord_vect_full_f
     | ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | out    | Full list of combined indices |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_lkpm_combined_index_f(context, cord_num, dim_cord_vect,  &
-     lkpm_combined_index) &
+integer function qmckl_compute_lkpm_combined_index_f( &
+     context, cord_num, dim_cord_vect,  lkpm_combined_index) &
      result(info)
   use qmckl
   implicit none
@@ -5390,7 +5418,8 @@ qmckl_exit_code qmckl_compute_lkpm_combined_index (
     | ~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
-integer function qmckl_compute_tmp_c_f(context, cord_num, elec_num, nucl_num, &
+integer function qmckl_compute_tmp_c_doc_f( &
+     context, cord_num, elec_num, nucl_num, &
      walk_num, een_rescaled_e, een_rescaled_n, tmp_c) &
      result(info)
   use qmckl
@@ -5445,7 +5474,7 @@ integer function qmckl_compute_tmp_c_f(context, cord_num, elec_num, nucl_num, &
 
   do nw=1, walk_num
   do i=0, cord_num-1
-  info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha,     &
+  info = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha,     &
                      een_rescaled_e(1,1,i,nw),LDA*1_8,                     &
                      een_rescaled_n(1,1,0,nw),LDB*1_8,                     &
                      beta,                                       &
@@ -5453,11 +5482,39 @@ integer function qmckl_compute_tmp_c_f(context, cord_num, elec_num, nucl_num, &
   end do
   end do
 
-end function qmckl_compute_tmp_c_f
+end function qmckl_compute_tmp_c_doc_f
     #+end_src
 
+#+CALL: generate_c_interface(table=qmckl_factor_tmp_c_args,rettyp=get_value("FRetType"),fname="qmckl_compute_tmp_c_doc")
+
+#+RESULTS:
+#+begin_src f90 :tangle (eval f) :comments org :exports none
+integer(c_int32_t) function qmckl_compute_tmp_c_doc &
+    (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e, een_rescaled_n, tmp_c) &
+    bind(C) result(info)
+
+  use, intrinsic :: iso_c_binding
+  implicit none
+
+  integer (c_int64_t) , intent(in)  , value :: context
+  integer (c_int64_t) , intent(in)  , value :: cord_num
+  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(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_doc_f
+  info = qmckl_compute_tmp_c_doc_f &
+         (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e, een_rescaled_n, tmp_c)
+
+end function qmckl_compute_tmp_c_doc
+#+end_src
+
+
     #+begin_src c :comments org :tangle (eval c) :noweb yes
-qmckl_exit_code qmckl_compute_tmp_c (
+qmckl_exit_code qmckl_compute_tmp_c_hpc (
 	  const qmckl_context context,
 	  const int64_t cord_num,
 	  const int64_t elec_num,
@@ -5467,55 +5524,55 @@ qmckl_exit_code qmckl_compute_tmp_c (
 	  const double* een_rescaled_n,
 	  double* const tmp_c ) {
 
-  qmckl_exit_code info;
-  int          i, j, a, l, kk, p, lmax, nw;
-  char         TransA, TransB;
-  double       alpha, beta;
-  int          M, N, K, LDA, LDB, LDC;
-
-  TransA = 'N';
-  TransB = 'N';
-  alpha  = 1.0;
-  beta   = 0.0;
-
   if (context == QMCKL_NULL_CONTEXT) {
     return QMCKL_INVALID_CONTEXT;
   }
 
   if (cord_num <= 0) {
     return QMCKL_INVALID_ARG_2;
-  }
+  } 
 
   if (elec_num <= 0) {
     return QMCKL_INVALID_ARG_3;
-  }
+  } 
 
   if (nucl_num <= 0) {
     return QMCKL_INVALID_ARG_4;
-  }
+  } 
 
-  M = elec_num;
-  N = nucl_num*(cord_num + 1);
-  K = elec_num;
+  if (walk_num <= 0) {
+    return QMCKL_INVALID_ARG_5;
+  } 
 
-  LDA = sizeof(een_rescaled_e)/sizeof(double);
-  LDB = sizeof(een_rescaled_n)/sizeof(double);
-  LDC = sizeof(tmp_c)/sizeof(double);
+  qmckl_exit_code info = QMCKL_SUCCESS;
 
-  for (int nw=0; nw < walk_num; ++nw) {
-    for (int i=0; i<cord_num; ++i){
-    info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha,     \
-                 //    &een_rescaled_e[0+0*elec_num+i*elec_num*elec_num+nw*elec_num*elec_num*(cord_num+1)],
-                       &een_rescaled_e[             i*elec_num*elec_num+nw*elec_num*elec_num*(cord_num+1)], \
+  const char  TransA = 'N';
+  const char  TransB = 'N';
+  const double alpha = 1.0;
+  const double beta  = 0.0;
+
+  const int64_t M = elec_num;
+  const int64_t N = nucl_num*(cord_num + 1);
+  const int64_t K = elec_num;
+
+  const int64_t LDA = elec_num;
+  const int64_t LDB = elec_num;
+  const int64_t LDC = elec_num;
+
+  const int64_t af = elec_num*elec_num;
+  const int64_t bf = elec_num*nucl_num*(cord_num+1);
+  const int64_t cf = bf;
+
+  for (int64_t nw=0; nw < walk_num; ++nw) {
+    for (int64_t i=0; i<cord_num; ++i){
+    info = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha,     \
+                       &(een_rescaled_e[af*(i+nw*(cord_num+1))]), \
                        LDA,                                        \
-                 //    &een_rescaled_n[0+0*elec_num+0*elec_num*nucl_num+nw*elec_num*nucl_num*(cord_num+1)],
-                       &een_rescaled_n[                                 nw*elec_num*nucl_num*(cord_num+1)], \
+                       &(een_rescaled_n[bf*nw]), \
                        LDB,                                        \
                        beta,                                       \
-                 //    &tmp_c[0+0*elec_num+0*elec_num*nucl_num+i*elec_num*nucl_num*(cord_num+1)+nw*elec_num*nucl_num*(cord_num+1)*cord_num],
-                       &tmp_c[                                 i*elec_num*nucl_num*(cord_num+1)+nw*elec_num*nucl_num*(cord_num+1)*cord_num],        \
+                       &(tmp_c[cf*(i+nw*cord_num)]), \
                        LDC);
-
     }
   }
 
@@ -5523,22 +5580,69 @@ qmckl_exit_code qmckl_compute_tmp_c (
 }
     #+end_src
 
-    #+CALL: generate_c_header(table=qmckl_factor_tmp_c_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
+    #+CALL: generate_c_header(table=qmckl_factor_tmp_c_args,rettyp=get_value("CRetType"),fname="qmckl_compute_tmp_c")
 
     #+RESULTS:
     #+begin_src c :tangle (eval h_func) :comments org
-    qmckl_exit_code qmckl_compute_tmp_c (
-	  const qmckl_context context,
-	  const int64_t cord_num,
-	  const int64_t elec_num,
-	  const int64_t nucl_num,
-	  const int64_t walk_num,
-	  const double* een_rescaled_e,
-	  const double* een_rescaled_n,
-	  double* const tmp_c );
+qmckl_exit_code qmckl_compute_tmp_c (
+          const qmckl_context context,
+          const int64_t cord_num,
+          const int64_t elec_num,
+          const int64_t nucl_num,
+          const int64_t walk_num,
+          const double* een_rescaled_e,
+          const double* een_rescaled_n,
+          double* const tmp_c ); 
     #+end_src
 
- TODO: FIX dtmp_c dimension in the table.
+#   #+CALL: generate_c_header(table=qmckl_factor_tmp_c_args,rettyp=get_value("CRetType"),fname="qmckl_compute_tmp_c_doc")
+
+    #+RESULTS:
+    #+begin_src c :tangle (eval h_private_func) :comments org
+qmckl_exit_code qmckl_compute_tmp_c_doc (
+          const qmckl_context context,
+          const int64_t cord_num,
+          const int64_t elec_num,
+          const int64_t nucl_num,
+          const int64_t walk_num,
+          const double* een_rescaled_e,
+          const double* een_rescaled_n,
+          double* const tmp_c ); 
+    #+end_src
+
+#   #+CALL: generate_c_header(table=qmckl_factor_tmp_c_args,rettyp=get_value("CRetType"),fname="qmckl_compute_tmp_c_hpc")
+
+    #+RESULTS:
+    #+begin_src c :tangle (eval h_private_func) :comments org
+qmckl_exit_code qmckl_compute_tmp_c_hpc (const qmckl_context context,
+                                         const int64_t cord_num,
+                                         const int64_t elec_num,
+                                         const int64_t nucl_num,
+                                         const int64_t walk_num,
+                                         const double* een_rescaled_e,
+                                         const double* een_rescaled_n,
+                                         double* const tmp_c ); 
+    #+end_src
+
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes
+qmckl_exit_code qmckl_compute_tmp_c (const qmckl_context context,
+                                     const int64_t cord_num,
+                                     const int64_t elec_num,
+                                     const int64_t nucl_num,
+                                     const int64_t walk_num,
+                                     const double* een_rescaled_e,
+                                     const double* een_rescaled_n,
+                                     double* const tmp_c )
+{
+#ifdef HAVE_HPC
+  return qmckl_compute_tmp_c_hpc(context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e, een_rescaled_n, tmp_c);
+#else
+  return qmckl_compute_tmp_c_doc(context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e, een_rescaled_n, tmp_c);
+#endif
+}
+    #+end_src
+    
 *** Compute tmp_c (OpenACC offload)
    :PROPERTIES:
    :Name:     qmckl_compute_tmp_c_acc_offload
@@ -5779,7 +5883,8 @@ qmckl_exit_code qmckl_compute_tmp_c_cublas_offload (
     | ~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
-integer function qmckl_compute_dtmp_c_f(context, cord_num, elec_num, nucl_num, &
+integer function qmckl_compute_dtmp_c_doc_f( &
+     context, cord_num, elec_num, nucl_num, &
      walk_num, een_rescaled_e_deriv_e, een_rescaled_n, dtmp_c) &
      result(info)
   use qmckl
@@ -5833,42 +5938,60 @@ integer function qmckl_compute_dtmp_c_f(context, cord_num, elec_num, nucl_num, &
   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(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)
+     do i=0, cord_num-1
+        info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha,  &
+             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
-  end do
-
-end function qmckl_compute_dtmp_c_f
+  
+end function qmckl_compute_dtmp_c_doc_f
     #+end_src
 
+#+CALL: generate_c_interface(table=qmckl_factor_dtmp_c_args,rettyp=get_value("FRetType"),fname="qmckl_compute_dtmp_c_doc")
+
+#+RESULTS:
+#+begin_src f90 :tangle (eval f) :comments org :exports none
+integer(c_int32_t) function qmckl_compute_dtmp_c_doc &
+    (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e_deriv_e, een_rescaled_n, dtmp_c) &
+    bind(C) result(info)
+
+  use, intrinsic :: iso_c_binding
+  implicit none
+
+  integer (c_int64_t) , intent(in)  , value :: context
+  integer (c_int64_t) , intent(in)  , value :: cord_num
+  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(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_doc_f
+  info = qmckl_compute_dtmp_c_doc_f &
+         (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e_deriv_e, een_rescaled_n, dtmp_c)
+
+end function qmckl_compute_dtmp_c_doc
+#+end_src
+
+    
     #+begin_src c :comments org :tangle (eval c) :noweb yes
-qmckl_exit_code qmckl_compute_dtmp_c (
-	  const qmckl_context context,
-	  const int64_t cord_num,
-	  const int64_t elec_num,
-	  const int64_t nucl_num,
-	  const int64_t walk_num,
-	  const double* een_rescaled_e_deriv_e,
-	  const double* een_rescaled_n,
-	  double* const dtmp_c ) {
-
-  qmckl_exit_code  info;
-  char             TransA, TransB;
-  double           alpha, beta;
-  int              M, N, K, LDA, LDB, LDC;
-
-  TransA = 'N';
-  TransB = 'N';
-  alpha = 1.0;
-  beta  = 0.0;
+qmckl_exit_code qmckl_compute_dtmp_c_hpc (const qmckl_context context,
+                                          const int64_t cord_num,
+                                          const int64_t elec_num,
+                                          const int64_t nucl_num,
+                                          const int64_t walk_num,
+                                          const double* een_rescaled_e_deriv_e,
+                                          const double* een_rescaled_n,
+                                          double* const dtmp_c )
+{
 
   if (context == QMCKL_NULL_CONTEXT) {
      return QMCKL_INVALID_CONTEXT;
-  }
+  } 
 
   if (cord_num <= 0) {
      return QMCKL_INVALID_ARG_2;
@@ -5880,29 +6003,40 @@ qmckl_exit_code qmckl_compute_dtmp_c (
 
   if (nucl_num <= 0) {
      return QMCKL_INVALID_ARG_4;
-  }
+  } 
 
-  M = 4*elec_num;
-  N = nucl_num*(cord_num + 1);
-  K = elec_num;
+  if (walk_num <= 0) {
+     return QMCKL_INVALID_ARG_5;
+  } 
 
-  LDA = 4*sizeof(een_rescaled_e_deriv_e)/sizeof(double);
-  LDB = sizeof(een_rescaled_n)/sizeof(double);
-  LDC = 4*sizeof(dtmp_c)/sizeof(double);
+  qmckl_exit_code  info = QMCKL_SUCCESS;
 
-  // TODO Rewrite manually and offload on GPU
-  for (int nw=0; nw < walk_num; ++nw) {
-    for (int i=0; nw < cord_num; ++i) {
-      info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha,        \
-              //&een_rescaled_e_deriv_e[0+0*elec_num+0*elec_num*4+i*elec_num*4*elec_num+nw*elec_num*4*elec_num*(cord_num+1)],
-                &een_rescaled_e_deriv_e[i*elec_num*4*elec_num+nw*elec_num*4*elec_num*(cord_num+1)],                \
+  const char  TransA = 'N';
+  const char  TransB = 'N';
+  const double alpha = 1.0;
+  const double beta  = 0.0;
+
+  const int64_t M = 4*elec_num;
+  const int64_t N = nucl_num*(cord_num + 1);
+  const int64_t K = elec_num;
+
+  const int64_t LDA = 4*elec_num;
+  const int64_t LDB = elec_num;
+  const int64_t LDC = 4*elec_num;
+
+  const int64_t af = elec_num*elec_num*4;
+  const int64_t bf = elec_num*nucl_num*(cord_num+1);
+  const int64_t cf = elec_num*4*nucl_num*(cord_num+1);
+
+  for (int64_t nw=0; nw < walk_num; ++nw) {
+    for (int64_t i=0; i < cord_num; ++i) {
+      info = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha,        \
+                &(een_rescaled_e_deriv_e[af*(i+nw*(cord_num+1))]), \
                 LDA,                                               \
-              //&een_rescaled_n[0+0*elec_num+0*elec_num*nucl_num+nw*elec_num*nucl_num*(cord_num+1)],
-                &een_rescaled_n[nw*elec_num*nucl_num*(cord_num+1)],   \
+                &(een_rescaled_n[bf*nw]),   \
                 LDB,                                               \
                 beta,                                              \
-              //&dtmp_c[0+0*elec_num+0*elec_num*4+0*elec_num*4*nucl_num+i*elec_num*4*nucl_num*(cord_num+1)+nw*elec_num*4*nucl_num*(cord_num+1)*cord_num],
-                &dtmp_c[i*elec_num*4*nucl_num*(cord_num+1)+nw*elec_num*4*nucl_num*(cord_num+1)*cord_num], \
+                &(dtmp_c[cf*(i+nw*cord_num)]), \
                 LDC);
     }
   }
@@ -5911,11 +6045,41 @@ qmckl_exit_code qmckl_compute_dtmp_c (
 }
     #+end_src
 
-    #+CALL: generate_c_header(table=qmckl_factor_dtmp_c_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
 
-    #+RESULTS:
+    #+begin_src c :comments org :tangle (eval c) :noweb yes
+qmckl_exit_code qmckl_compute_dtmp_c (const qmckl_context context,
+                                      const int64_t cord_num,
+                                      const int64_t elec_num,
+                                      const int64_t nucl_num,
+                                      const int64_t walk_num,
+                                      const double* een_rescaled_e_deriv_e,
+                                      const double* een_rescaled_n,
+                                      double* const dtmp_c )
+{
+#ifdef HAVE_HPC
+  return qmckl_compute_dtmp_c_hpc (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e_deriv_e,
+                                   een_rescaled_n, dtmp_c );
+#else
+  return qmckl_compute_dtmp_c_doc (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e_deriv_e,
+                                   een_rescaled_n, dtmp_c );
+#endif
+}
+    #+end_src
+
+
     #+begin_src c :tangle (eval h_func) :comments org
-    qmckl_exit_code qmckl_compute_dtmp_c (
+qmckl_exit_code qmckl_compute_dtmp_c (const qmckl_context context,
+                                      const int64_t cord_num,
+                                      const int64_t elec_num,
+                                      const int64_t nucl_num,
+                                      const int64_t walk_num,
+                                      const double* een_rescaled_e_deriv_e,
+                                      const double* een_rescaled_n,
+                                      double* const dtmp_c );
+    #+end_src
+
+    #+begin_src c :tangle (eval h_private_func) :comments org
+    qmckl_exit_code qmckl_compute_dtmp_c_doc (
 	  const qmckl_context context,
 	  const int64_t cord_num,
 	  const int64_t elec_num,
@@ -5926,6 +6090,18 @@ qmckl_exit_code qmckl_compute_dtmp_c (
 	  double* const dtmp_c );
     #+end_src
 
+    #+begin_src c :tangle (eval h_private_func) :comments org
+qmckl_exit_code qmckl_compute_dtmp_c_hpc (
+	  const qmckl_context context,
+	  const int64_t cord_num,
+	  const int64_t elec_num,
+	  const int64_t nucl_num,
+	  const int64_t walk_num,
+	  const double* een_rescaled_e_deriv_e,
+	  const double* een_rescaled_n,
+	  double* const dtmp_c );
+    #+end_src
+    
 *** Compute dtmp_c (OpenACC offload)
    :PROPERTIES:
    :Name:     qmckl_compute_dtmp_c_acc_offload
@@ -6260,7 +6436,7 @@ qmckl_get_jastrow_factor_een(qmckl_context context,
   rc = qmckl_provide_factor_een(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int64_t sze = ctx->electron.walk_num;
@@ -6291,7 +6467,7 @@ qmckl_exit_code qmckl_provide_factor_een(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   /* Check if en rescaled distance is provided */
@@ -6379,9 +6555,10 @@ qmckl_exit_code qmckl_provide_factor_een(qmckl_context context)
     | ~factor_een~          | ~double[walk_num]~                                 | out    | Electron-nucleus jastrow             |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_factor_een_naive_f(context, walk_num, elec_num, nucl_num, cord_num,&
-                                           dim_cord_vect, cord_vect_full, lkpm_combined_index,    &
-                                           een_rescaled_e, een_rescaled_n, factor_een) &
+integer function qmckl_compute_factor_een_naive_f( &
+     context, walk_num, elec_num, nucl_num, cord_num,&
+     dim_cord_vect, cord_vect_full, lkpm_combined_index, &
+     een_rescaled_e, een_rescaled_n, factor_een) &
      result(info)
   use qmckl
   implicit none
@@ -6546,9 +6723,10 @@ end function qmckl_compute_factor_een_naive_f
     | ~factor_een~          | ~double[walk_num]~                                               | out                             | Electron-nucleus jastrow             |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_factor_een_f(context, walk_num, elec_num, nucl_num, cord_num,   &
-                                           dim_cord_vect, cord_vect_full, lkpm_combined_index, &
-                                           tmp_c, een_rescaled_n, factor_een) &
+integer function qmckl_compute_factor_een_f( &
+     context, walk_num, elec_num, nucl_num, cord_num,   &
+     dim_cord_vect, cord_vect_full, lkpm_combined_index, &
+     tmp_c, een_rescaled_n, factor_een) &
      result(info)
   use qmckl
   implicit none
@@ -6759,7 +6937,7 @@ qmckl_get_jastrow_factor_een_deriv_e(qmckl_context context,
   rc = qmckl_provide_factor_een_deriv_e(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int64_t sze = ctx->electron.walk_num * 4 * ctx->electron.num;
@@ -6790,7 +6968,7 @@ qmckl_exit_code qmckl_provide_factor_een_deriv_e(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   /* Check if en rescaled distance is provided */
@@ -6894,10 +7072,10 @@ qmckl_exit_code qmckl_provide_factor_een_deriv_e(qmckl_context context)
     | ~factor_een_deriv_e~     | ~double[walk_num][4][elec_num]~                       | out    | Electron-nucleus jastrow             |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_factor_een_deriv_e_naive_f(context, walk_num, elec_num, nucl_num, cord_num,    &
-                                           dim_cord_vect, cord_vect_full, lkpm_combined_index,                &
-                                           een_rescaled_e, een_rescaled_n,                                    &
-                                           een_rescaled_e_deriv_e, een_rescaled_n_deriv_e, factor_een_deriv_e)&
+integer function qmckl_compute_factor_een_deriv_e_naive_f( &
+     context, walk_num, elec_num, nucl_num, cord_num, dim_cord_vect, &
+     cord_vect_full, lkpm_combined_index, een_rescaled_e, een_rescaled_n, &
+     een_rescaled_e_deriv_e, een_rescaled_n_deriv_e, factor_een_deriv_e)&
      result(info)
   use qmckl
   implicit none
@@ -7093,9 +7271,10 @@ end function qmckl_compute_factor_een_deriv_e_naive_f
 
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_factor_een_deriv_e_f(context, walk_num, elec_num, nucl_num, cord_num, dim_cord_vect,            &
-                                           cord_vect_full, lkpm_combined_index,                                       &
-                                           tmp_c, dtmp_c, een_rescaled_n, een_rescaled_n_deriv_e, factor_een_deriv_e) &
+integer function qmckl_compute_factor_een_deriv_e_f( &
+     context, walk_num, elec_num, nucl_num, &
+     cord_num, dim_cord_vect, cord_vect_full, lkpm_combined_index, &
+     tmp_c, dtmp_c, een_rescaled_n, een_rescaled_n_deriv_e, factor_een_deriv_e)&
      result(info)
   use qmckl
   implicit none
diff --git a/org/qmckl_local_energy.org b/org/qmckl_local_energy.org
index ccd933d..0f1fdc8 100644
--- a/org/qmckl_local_energy.org
+++ b/org/qmckl_local_energy.org
@@ -226,7 +226,7 @@ qmckl_exit_code qmckl_get_kinetic_energy(qmckl_context context, double * const k
   rc = qmckl_provide_kinetic_energy(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->electron.walk_num * sizeof(double);
@@ -250,7 +250,7 @@ qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context) {
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if(!(ctx->nucleus.provided)) {
@@ -549,37 +549,28 @@ end function qmckl_compute_kinetic_energy_f
 
 *** Test
      #+begin_src c :tangle (eval c_test) :exports none
-
-#define walk_num chbrclf_walk_num
-#define elec_num chbrclf_elec_num
-#define shell_num chbrclf_shell_num
-#define ao_num chbrclf_ao_num
-
-int64_t elec_up_num   = chbrclf_elec_up_num;
-int64_t elec_dn_num   = chbrclf_elec_dn_num;
 double* elec_coord    = &(chbrclf_elec_coord[0][0][0]);
-const int64_t   nucl_num      = chbrclf_nucl_num;
 const double*   nucl_charge   = chbrclf_charge;
 const double*   nucl_coord    = &(chbrclf_nucl_coord[0][0]);
 
-rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
+rc = qmckl_set_electron_num (context, chbrclf_elec_up_num, chbrclf_elec_dn_num);
 assert (rc == QMCKL_SUCCESS);
 
-rc = qmckl_set_electron_walk_num (context, walk_num);
+rc = qmckl_set_electron_walk_num (context, chbrclf_walk_num);
 assert (rc == QMCKL_SUCCESS);
 
 assert(qmckl_electron_provided(context));
 
-rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3);
+rc = qmckl_set_electron_coord (context, 'N', elec_coord, chbrclf_walk_num*chbrclf_elec_num*3);
 assert(rc == QMCKL_SUCCESS);
 
-rc = qmckl_set_nucleus_num (context, nucl_num);
+rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num);
 assert(rc == QMCKL_SUCCESS);
 
-rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), nucl_num*3);
+rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), chbrclf_nucl_num*3);
 assert(rc == QMCKL_SUCCESS);
 
-rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
+rc = qmckl_set_nucleus_charge(context, nucl_charge, chbrclf_nucl_num);
 assert(rc == QMCKL_SUCCESS);
 
 assert(qmckl_nucleus_provided(context));
@@ -611,11 +602,11 @@ rc = qmckl_set_ao_basis_prim_num (context, chbrclf_prim_num);
 assert(rc == QMCKL_SUCCESS);
 assert(!qmckl_ao_basis_provided(context));
 
-rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, nucl_num);
+rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, chbrclf_nucl_num);
 assert(rc == QMCKL_SUCCESS);
 assert(!qmckl_ao_basis_provided(context));
 
-rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, nucl_num);
+rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, chbrclf_nucl_num);
 assert(rc == QMCKL_SUCCESS);
 assert(!qmckl_ao_basis_provided(context));
 
@@ -655,10 +646,10 @@ assert(rc == QMCKL_SUCCESS);
 assert(qmckl_ao_basis_provided(context));
 
 
-double ao_vgl[5][walk_num][elec_num][chbrclf_ao_num];
+double ao_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_ao_num];
 
-rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0][0]),
-                               (int64_t) 5*walk_num*elec_num*chbrclf_ao_num);
+rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]),
+                               (int64_t) 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_ao_num);
 assert (rc == QMCKL_SUCCESS);
 
 /* Set up MO data */
@@ -673,31 +664,31 @@ assert (rc == QMCKL_SUCCESS);
 
 assert(qmckl_mo_basis_provided(context));
 
-double mo_vgl[5][elec_num][chbrclf_mo_num];
-rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0]));
+double mo_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_mo_num];
+rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_mo_num);
 assert (rc == QMCKL_SUCCESS);
 
 /* Set up determinant data */
 
-const int64_t det_num_alpha = 1;
-const int64_t det_num_beta  = 1;
-int64_t mo_index_alpha[det_num_alpha][walk_num][elec_up_num];
-int64_t mo_index_beta[det_num_alpha][walk_num][elec_dn_num];
+#define det_num_alpha 1
+#define det_num_beta  1
+int64_t mo_index_alpha[det_num_alpha][chbrclf_walk_num][chbrclf_elec_up_num];
+int64_t mo_index_beta[det_num_alpha][chbrclf_walk_num][chbrclf_elec_dn_num];
 
 int i, j, k;
 for(k = 0; k < det_num_alpha; ++k)
-  for(i = 0; i < walk_num; ++i)
-    for(j = 0; j < elec_up_num; ++j)
+  for(i = 0; i < chbrclf_walk_num; ++i)
+    for(j = 0; j < chbrclf_elec_up_num; ++j)
       mo_index_alpha[k][i][j] = j + 1;
 for(k = 0; k < det_num_beta; ++k)
-  for(i = 0; i < walk_num; ++i)
-    for(j = 0; j < elec_up_num; ++j)
+  for(i = 0; i < chbrclf_walk_num; ++i)
+    for(j = 0; j < chbrclf_elec_up_num; ++j)
       mo_index_beta[k][i][j] = j + 1;
 
 rc = qmckl_set_determinant_type (context, typ);
 assert(rc == QMCKL_SUCCESS);
 
-rc = qmckl_set_determinant_walk_num (context, walk_num);
+rc = qmckl_set_determinant_walk_num (context, chbrclf_walk_num);
 assert (rc == QMCKL_SUCCESS);
 
 rc = qmckl_set_determinant_det_num_alpha (context, det_num_alpha);
@@ -714,8 +705,8 @@ assert (rc == QMCKL_SUCCESS);
 
 // Get alpha determinant
 
-double det_vgl_alpha[det_num_alpha][walk_num][5][elec_up_num][elec_up_num];
-double det_vgl_beta[det_num_beta][walk_num][5][elec_dn_num][elec_dn_num];
+double det_vgl_alpha[det_num_alpha][chbrclf_walk_num][5][chbrclf_elec_up_num][chbrclf_elec_up_num];
+double det_vgl_beta[det_num_beta][chbrclf_walk_num][5][chbrclf_elec_dn_num][chbrclf_elec_dn_num];
 
 rc = qmckl_get_det_vgl_alpha(context, &(det_vgl_alpha[0][0][0][0][0]));
 assert (rc == QMCKL_SUCCESS);
@@ -725,8 +716,8 @@ assert (rc == QMCKL_SUCCESS);
 
 // Get adjoint of the slater-determinant
 
-double det_inv_matrix_alpha[det_num_alpha][walk_num][elec_up_num][elec_up_num];
-double det_inv_matrix_beta[det_num_beta][walk_num][elec_dn_num][elec_dn_num];
+double det_inv_matrix_alpha[det_num_alpha][chbrclf_walk_num][chbrclf_elec_up_num][chbrclf_elec_up_num];
+double det_inv_matrix_beta[det_num_beta][chbrclf_walk_num][chbrclf_elec_dn_num][chbrclf_elec_dn_num];
 
 rc = qmckl_get_det_inv_matrix_alpha(context, &(det_inv_matrix_alpha[0][0][0][0]));
 assert (rc == QMCKL_SUCCESS);
@@ -736,7 +727,7 @@ assert (rc == QMCKL_SUCCESS);
 
 // Calculate the Kinetic energy
 
-double kinetic_energy[walk_num];
+double kinetic_energy[chbrclf_walk_num];
 
 rc = qmckl_get_kinetic_energy(context, &(kinetic_energy[0]));
 assert (rc == QMCKL_SUCCESS);
@@ -799,7 +790,7 @@ qmckl_exit_code qmckl_get_potential_energy(qmckl_context context, double * const
   rc = qmckl_provide_potential_energy(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->electron.walk_num * sizeof(double);
@@ -822,7 +813,7 @@ qmckl_exit_code qmckl_provide_potential_energy(qmckl_context context) {
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   qmckl_exit_code rc;
@@ -1034,7 +1025,7 @@ end function qmckl_compute_potential_energy_f
      #+begin_src c :tangle (eval c_test) :exports none
 // Calculate the Potential energy
 
-double potential_energy[walk_num];
+double potential_energy[chbrclf_walk_num];
 
 rc = qmckl_get_potential_energy(context, &(potential_energy[0]));
 assert (rc == QMCKL_SUCCESS);
@@ -1083,7 +1074,7 @@ qmckl_exit_code qmckl_get_local_energy(qmckl_context context, double * const loc
   rc = qmckl_provide_local_energy(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->electron.walk_num * sizeof(double);
@@ -1106,7 +1097,7 @@ qmckl_exit_code qmckl_provide_local_energy(qmckl_context context) {
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if(!(ctx->nucleus.provided)) {
@@ -1290,7 +1281,7 @@ end function qmckl_compute_local_energy_f
      #+begin_src c :tangle (eval c_test) :exports none
 // Calculate the Local energy
 
-double local_energy[walk_num];
+double local_energy[chbrclf_walk_num];
 
 rc = qmckl_get_local_energy(context, &(local_energy[0]));
 assert (rc == QMCKL_SUCCESS);
@@ -1339,7 +1330,7 @@ qmckl_exit_code qmckl_get_drift_vector(qmckl_context context, double * const dri
   rc = qmckl_provide_drift_vector(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   size_t sze = ctx->electron.walk_num * ctx->electron.num * 3 * sizeof(double);
@@ -1362,7 +1353,7 @@ qmckl_exit_code qmckl_provide_drift_vector(qmckl_context context) {
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if(!(ctx->nucleus.provided)) {
@@ -1645,7 +1636,7 @@ end function qmckl_compute_drift_vector_f
      #+begin_src c :tangle (eval c_test) :exports none
 // Calculate the Drift vector
 
-double drift_vector[walk_num][elec_num][3];
+double drift_vector[chbrclf_walk_num][chbrclf_elec_num][3];
 
 rc = qmckl_get_drift_vector(context, &(drift_vector[0][0][0]));
 assert (rc == QMCKL_SUCCESS);
diff --git a/org/qmckl_memory.org b/org/qmckl_memory.org
index b24de98..1bc7006 100644
--- a/org/qmckl_memory.org
+++ b/org/qmckl_memory.org
@@ -116,7 +116,7 @@ void* qmckl_malloc(qmckl_context context, const qmckl_memory_info_struct info) {
 
   assert (qmckl_context_check(context) != QMCKL_NULL_CONTEXT);
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
 
   /* Allocate memory and zero it */
   void * pointer = malloc(info.size);
@@ -217,7 +217,7 @@ qmckl_exit_code qmckl_free(qmckl_context context, void * const ptr) {
                           "NULL pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
 
   qmckl_lock(context);
   {
diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org
index d4efb41..d920396 100644
--- a/org/qmckl_mo.org
+++ b/org/qmckl_mo.org
@@ -84,13 +84,14 @@ int main() {
   The following arrays are stored in the context:
 
 
-  |---------------+--------------------+----------------------|
-  | ~mo_num~      |                    | Number of MOs        |
-  | ~coefficient~ | ~[mo_num][ao_num]~ | Orbital coefficients |
+  |-----------------+--------------------+----------------------------------------|
+  | ~mo_num~        |                    | Number of MOs                          |
+  | ~coefficient~   | ~[mo_num][ao_num]~ | Orbital coefficients                   |
+  | ~coefficient_t~ | ~[ao_num][mo_num]~ | Transposed of the Orbital coefficients |
+  |-----------------+--------------------+----------------------------------------|
 
   Computed data:
 
-  |---------------+--------------------------+-------------------------------------------------------------------------------------|
   |---------------+--------------------------+-------------------------------------------------------------------------------------|
   | ~mo_vgl~      | ~[point_num][5][mo_num]~ | Value, gradients, Laplacian of the MOs at point positions                           |
   | ~mo_vgl_date~ | ~uint64_t~               | Late modification date of Value, gradients, Laplacian of the MOs at point positions |
@@ -101,9 +102,10 @@ int main() {
    #+begin_src c :comments org :tangle (eval h_private_type)
 typedef struct qmckl_mo_basis_struct {
   int64_t   mo_num;
-  double  * coefficient;
+  double * restrict coefficient;
+  double * restrict coefficient_t;
 
-  double  * mo_vgl;
+  double * restrict mo_vgl;
   uint64_t  mo_vgl_date;
 
   int32_t   uninitialized;
@@ -129,7 +131,7 @@ qmckl_exit_code qmckl_init_mo_basis(qmckl_context context) {
     return false;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   ctx->mo_basis.uninitialized = (1 << 2) - 1;
@@ -156,10 +158,9 @@ qmckl_get_mo_basis_mo_num (const qmckl_context context,
                            QMCKL_INVALID_CONTEXT,
                            "qmckl_get_mo_basis_mo_num",
                            NULL);
-   return (int64_t) 0;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1;
@@ -198,7 +199,7 @@ qmckl_get_mo_basis_coefficient (const qmckl_context context,
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 1;
@@ -246,7 +247,7 @@ bool qmckl_mo_basis_provided(const qmckl_context context) {
     return false;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   return ctx->mo_basis.provided;
@@ -271,7 +272,7 @@ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
   return QMCKL_NULL_CONTEXT;
  }
 
-qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
    #+end_src
 
    #+NAME:post
@@ -352,9 +353,37 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
                            NULL);
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
+  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+  mem_info.size = ctx->ao_basis.ao_num * ctx->mo_basis.mo_num * sizeof(double);
+  double* new_array = (double*) qmckl_malloc(context, mem_info);
+  if (new_array == NULL) {
+    return qmckl_failwith( context,
+                           QMCKL_ALLOCATION_FAILED,
+                           "qmckl_finalize_mo_basis",
+                           NULL);
+  }
+
+  assert (ctx->mo_basis.coefficient != NULL);
+  
+  if (ctx->mo_basis.coefficient_t != NULL) {
+    qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient);
+    if (rc != QMCKL_SUCCESS) {
+      return qmckl_failwith( context, rc,
+                             "qmckl_finalize_mo_basis",
+                             NULL);
+    }
+  }
+
+  for (int64_t i=0 ; i<ctx->ao_basis.ao_num ; ++i) {
+    for (int64_t j=0 ; j<ctx->mo_basis.mo_num ; ++j) {
+      new_array[i*ctx->mo_basis.mo_num + j] = ctx->mo_basis.coefficient[j*ctx->ao_basis.ao_num + i];
+    }
+  }
+
+  ctx->mo_basis.coefficient_t = new_array;
   qmckl_exit_code rc = QMCKL_SUCCESS;
   return rc;
 }
@@ -367,11 +396,18 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
 *** Get
 
     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
-qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_vgl);
+qmckl_exit_code
+qmckl_get_mo_basis_mo_vgl(qmckl_context context,
+                          double* const mo_vgl,
+                          const int64_t size_max);
     #+end_src
 
     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
-qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_vgl) {
+qmckl_exit_code
+qmckl_get_mo_basis_mo_vgl(qmckl_context context,
+                          double* const mo_vgl,
+                          const int64_t size_max)
+{
 
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return QMCKL_NULL_CONTEXT;
@@ -385,10 +421,16 @@ qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_v
   rc = qmckl_provide_mo_vgl(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
-  size_t sze = 5 * ctx->point.num * ctx->mo_basis.mo_num;
+  const int64_t sze = ctx->point.num * 5 * ctx->mo_basis.mo_num;
+  if (size_max < sze) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_get_mo_basis_mo_vgl",
+                           "input array too small");
+  }
   memcpy(mo_vgl, ctx->mo_basis.mo_vgl, sze * sizeof(double));
 
   return QMCKL_SUCCESS;
@@ -396,17 +438,84 @@ qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_v
     #+end_src
 
     #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
-    interface
-      integer(c_int32_t) function qmckl_get_mo_basis_vgl (context, mo_vgl) &
-          bind(C)
-        use, intrinsic :: iso_c_binding
-        import
-        implicit none
+  interface
+     integer(c_int32_t) function qmckl_get_mo_basis_mo_vgl (context, &
+          mo_vgl, size_max) bind(C)
+       use, intrinsic :: iso_c_binding
+       import
+       implicit none
+       
+       integer (c_int64_t) , intent(in)  , value :: context
+       double precision,     intent(out)         :: mo_vgl(*)
+       integer (c_int64_t) , intent(in)  , value :: size_max
+     end function qmckl_get_mo_basis_mo_vgl
+  end interface
+    #+end_src
 
-        integer (c_int64_t) , intent(in)  , value :: context
-        double precision,     intent(out)         :: mo_vgl(*)
-      end function
-    end interface
+    Uses the given array to compute the VGL.
+
+    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code
+qmckl_get_mo_basis_mo_vgl_inplace (qmckl_context context,
+                                   double* const mo_vgl,
+                                   const int64_t size_max);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code
+qmckl_get_mo_basis_mo_vgl_inplace (qmckl_context context,
+                                   double* const mo_vgl,
+                                   const int64_t size_max)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_CONTEXT,
+                           "qmckl_get_mo_basis_mo_vgl",
+                           NULL);
+  }
+
+  qmckl_exit_code rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  const int64_t sze = ctx->mo_basis.mo_num * 5 * ctx->point.num;
+  if (size_max < sze) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_get_mo_basis_mo_vgl",
+                           "input array too small");
+  }
+
+  rc = qmckl_context_touch(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  double* old_array = ctx->mo_basis.mo_vgl;
+
+  ctx->mo_basis.mo_vgl = mo_vgl;
+
+  rc = qmckl_provide_mo_vgl(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  ctx->mo_basis.mo_vgl = old_array;
+
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
+    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
+  interface
+     integer(c_int32_t) function qmckl_get_mo_basis_mo_vgl_inplace (context, &
+          mo_vgl, size_max) bind(C)
+       use, intrinsic :: iso_c_binding
+       import
+       implicit none
+       integer (c_int64_t) , intent(in)  , value :: context
+       double precision,     intent(out)         :: mo_vgl(*)
+       integer (c_int64_t) , intent(in)  , value :: size_max
+     end function qmckl_get_mo_basis_mo_vgl_inplace
+  end interface
     #+end_src
 
 *** Provide
@@ -424,7 +533,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
     return QMCKL_NULL_CONTEXT;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if (!ctx->ao_basis.provided) {
@@ -462,19 +571,19 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
       if (mo_vgl == NULL) {
         return qmckl_failwith( context,
                                QMCKL_ALLOCATION_FAILED,
-                               "qmckl_mo_basis_vgl",
+                               "qmckl_mo_basis_mo_vgl",
                                NULL);
       }
       ctx->mo_basis.mo_vgl = mo_vgl;
     }
 
-    rc = qmckl_compute_mo_basis_vgl(context,
-                                    ctx->ao_basis.ao_num,
-                                    ctx->mo_basis.mo_num,
-                                    ctx->point.num,
-                                    ctx->mo_basis.coefficient,
-                                    ctx->ao_basis.ao_vgl,
-                                    ctx->mo_basis.mo_vgl);
+    rc = qmckl_compute_mo_basis_mo_vgl(context,
+                                       ctx->ao_basis.ao_num,
+                                       ctx->mo_basis.mo_num,
+                                       ctx->point.num,
+                                       ctx->mo_basis.coefficient_t,
+                                       ctx->ao_basis.ao_vgl,
+                                       ctx->mo_basis.mo_vgl);
     if (rc != QMCKL_SUCCESS) {
       return rc;
     }
@@ -488,25 +597,33 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
 
 *** Compute
    :PROPERTIES:
-   :Name:     qmckl_compute_mo_basis_vgl
+   :Name:     qmckl_compute_mo_basis_mo_vgl
    :CRetType: qmckl_exit_code
    :FRetType: qmckl_exit_code
    :END:
 
-    #+NAME: qmckl_mo_basis_vgl_args
-   | ~qmckl_context~ | ~context~                         | in  | Global state                              |
-   | ~int64_t~       | ~ao_num~                          | in  | Number of AOs                             |
-   | ~int64_t~       | ~mo_num~                          | in  | Number of MOs                             |
-   | ~int64_t~       | ~point_num~                       | in  | Number of points                          |
-   | ~double~        | ~coef_normalized[mo_num][ao_num]~ | in  | AO to MO transformation matrix            |
-   | ~double~        | ~ao_vgl[point_num][5][ao_num]~    | in  | Value, gradients and Laplacian of the AOs |
-   | ~double~        | ~mo_vgl[point_num][5][mo_num]~    | out | Value, gradients and Laplacian of the MOs |
+    #+NAME: qmckl_mo_basis_mo_vgl_args
+    | Variable            | Type                           | In/Out | Description                                     |
+    |---------------------+--------------------------------+--------+-------------------------------------------------|
+    | ~context~           | ~qmckl_context~                | in     | Global state                                    |
+    | ~ao_num~            | ~int64_t~                      | in     | Number of AOs                                   |
+    | ~mo_num~            | ~int64_t~                      | in     | Number of MOs                                   |
+    | ~point_num~         | ~int64_t~                      | in     | Number of points                                |
+    | ~coef_normalized_t~ | ~double[mo_num][ao_num]~       | in     | Transpose of the AO to MO transformation matrix |
+    | ~ao_vgl~            | ~double[point_num][5][ao_num]~ | in     | Value, gradients and Laplacian of the AOs       |
+    | ~mo_vgl~            | ~double[point_num][5][mo_num]~ | out    | Value, gradients and Laplacian of the MOs       |
 
 
+    The matrix of AO values is very sparse, so we use a sparse-dense
+    matrix multiplication instead of a dgemm, as exposed in
+    https://dx.doi.org/10.1007/978-3-642-38718-0_14.
+
+    
+    
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_mo_basis_vgl_f(context, &
+integer function qmckl_compute_mo_basis_mo_vgl_doc_f(context, &
      ao_num, mo_num, point_num, &
-     coef_normalized, ao_vgl, mo_vgl) &
+     coef_normalized_t, ao_vgl, mo_vgl) &
      result(info)
   use qmckl
   implicit none
@@ -514,55 +631,69 @@ integer function qmckl_compute_mo_basis_vgl_f(context, &
   integer*8             , intent(in)  :: ao_num, mo_num
   integer*8             , intent(in)  :: point_num
   double precision      , intent(in)  :: ao_vgl(ao_num,5,point_num)
-  double precision      , intent(in)  :: coef_normalized(ao_num,mo_num)
+  double precision      , intent(in)  :: coef_normalized_t(mo_num,ao_num)
   double precision      , intent(out) :: mo_vgl(mo_num,5,point_num)
-  character                     :: TransA, TransB
-  double precision              :: alpha, beta
-  integer*8                     :: M, N, K, LDA, LDB, LDC, i,j
+  integer*8 :: i,j,k
+  double precision :: c1, c2, c3, c4, c5
 
-  TransA = 'T'
-  TransB = 'N'
-  M = mo_num
-  N = point_num*5_8
-  K = int(ao_num,8)
-  alpha = 1.d0
-  beta  = 0.d0
-  LDA = size(coef_normalized,1)
-  LDB = size(ao_vgl,1)
-  LDC = size(mo_vgl,1)
+  do j=1,point_num
+     mo_vgl(:,:,j) = 0.d0
+     do k=1,ao_num
+        if (ao_vgl(k,1,j) /= 0.d0) then
+           c1 = ao_vgl(k,1,j)
+           c2 = ao_vgl(k,2,j)
+           c3 = ao_vgl(k,3,j)
+           c4 = ao_vgl(k,4,j)
+           c5 = ao_vgl(k,5,j)
+           do i=1,mo_num
+              mo_vgl(i,1,j) = mo_vgl(i,1,j) + coef_normalized_t(i,k) * c1
+              mo_vgl(i,2,j) = mo_vgl(i,2,j) + coef_normalized_t(i,k) * c2
+              mo_vgl(i,3,j) = mo_vgl(i,3,j) + coef_normalized_t(i,k) * c3
+              mo_vgl(i,4,j) = mo_vgl(i,4,j) + coef_normalized_t(i,k) * c4
+              mo_vgl(i,5,j) = mo_vgl(i,5,j) + coef_normalized_t(i,k) * c5
+           end do
+        end if
+     end do
+  end do
 
-  info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, &
-       coef_normalized, int(size(coef_normalized,1),8),      &
-       ao_vgl, int(size(ao_vgl,1),8), beta,                  &
-       mo_vgl,LDC)
-
-  info = QMCKL_SUCCESS
-
-end function qmckl_compute_mo_basis_vgl_f
+end function qmckl_compute_mo_basis_mo_vgl_doc_f
     #+end_src
 
-   #+CALL: generate_c_header(table=qmckl_mo_basis_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_vgl"))
+    #+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl"))
 
    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
-   qmckl_exit_code qmckl_compute_mo_basis_vgl (
-	 const qmckl_context context,
-	 const int64_t ao_num,
-	 const int64_t mo_num,
-	 const int64_t point_num,
-	 const double* coef_normalized,
-	 const double* ao_vgl,
-	 double* const mo_vgl );
+   qmckl_exit_code qmckl_compute_mo_basis_mo_vgl (
+         const qmckl_context context,
+         const int64_t ao_num,
+         const int64_t mo_num,
+         const int64_t point_num,
+         const double* coef_normalized_t,
+         const double* ao_vgl,
+         double* const mo_vgl ); 
    #+end_src
 
+   #+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_doc"))
 
-    #+CALL: generate_c_interface(table=qmckl_mo_basis_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_vgl"))
+   #+RESULTS:
+   #+begin_src c :tangle (eval h_func) :comments org
+   qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_doc (
+         const qmckl_context context,
+         const int64_t ao_num,
+         const int64_t mo_num,
+         const int64_t point_num,
+         const double* coef_normalized_t,
+         const double* ao_vgl,
+         double* const mo_vgl ); 
+   #+end_src
 
+   #+CALL: generate_c_interface(table=qmckl_mo_basis_mo_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_doc"))
+   
     #+RESULTS:
     #+begin_src f90 :tangle (eval f) :comments org :exports none
-    integer(c_int32_t) function qmckl_compute_mo_basis_vgl &
-	(context, ao_num, mo_num, point_num, coef_normalized, ao_vgl, mo_vgl) &
-	bind(C) result(info)
+    integer(c_int32_t) function qmckl_compute_mo_basis_mo_vgl_doc &
+        (context, ao_num, mo_num, point_num, coef_normalized_t, ao_vgl, mo_vgl) &
+        bind(C) result(info)
 
       use, intrinsic :: iso_c_binding
       implicit none
@@ -571,15 +702,176 @@ end function qmckl_compute_mo_basis_vgl_f
       integer (c_int64_t) , intent(in)  , value :: ao_num
       integer (c_int64_t) , intent(in)  , value :: mo_num
       integer (c_int64_t) , intent(in)  , value :: point_num
-      real    (c_double ) , intent(in)          :: coef_normalized(ao_num,mo_num)
+      real    (c_double ) , intent(in)          :: coef_normalized_t(ao_num,mo_num)
       real    (c_double ) , intent(in)          :: ao_vgl(ao_num,5,point_num)
       real    (c_double ) , intent(out)         :: mo_vgl(mo_num,5,point_num)
 
-      integer(c_int32_t), external :: qmckl_compute_mo_basis_vgl_f
-      info = qmckl_compute_mo_basis_vgl_f &
-	     (context, ao_num, mo_num, point_num, coef_normalized, ao_vgl, mo_vgl)
+      integer(c_int32_t), external :: qmckl_compute_mo_basis_mo_vgl_doc_f
+      info = qmckl_compute_mo_basis_mo_vgl_doc_f &
+             (context, ao_num, mo_num, point_num, coef_normalized_t, ao_vgl, mo_vgl)
 
-    end function qmckl_compute_mo_basis_vgl
+    end function qmckl_compute_mo_basis_mo_vgl_doc
+    #+end_src
+
+    #+begin_src c :tangle (eval c) :comments org                                                                             
+qmckl_exit_code
+qmckl_compute_mo_basis_mo_vgl (const qmckl_context context,
+                            const int64_t ao_num,
+                            const int64_t mo_num,
+                            const int64_t point_num,
+                            const double* coef_normalized_t,
+                            const double* ao_vgl,
+                            double* const mo_vgl )
+{
+#ifdef HAVE_HPC
+  return qmckl_compute_mo_basis_mo_vgl_hpc (context, ao_num, mo_num, point_num, coef_normalized_t, ao_vgl, mo_vgl);
+#else
+  return qmckl_compute_mo_basis_mo_vgl_doc (context, ao_num, mo_num, point_num, coef_normalized_t, ao_vgl, mo_vgl);
+#endif
+}
+    #+end_src
+    
+    
+*** HPC version
+    
+
+    #+begin_src c :tangle (eval h_func) :comments org                                                                             
+#ifdef HAVE_HPC
+qmckl_exit_code
+qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
+                                   const int64_t ao_num,
+                                   const int64_t mo_num,
+                                   const int64_t point_num,
+                                   const double* coef_normalized_t,
+                                   const double* ao_vgl,
+                                   double* const mo_vgl );
+#endif
+    #+end_src
+
+    #+begin_src c :tangle (eval c) :comments org                                                                             
+#ifdef HAVE_HPC
+qmckl_exit_code
+qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
+                                   const int64_t ao_num,
+                                   const int64_t mo_num,
+                                   const int64_t point_num,
+                                   const double* restrict coef_normalized_t,
+                                   const double* restrict ao_vgl,
+                                   double* restrict const mo_vgl )
+{
+#ifdef HAVE_OPENMP
+  #pragma omp parallel for
+#endif
+  for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
+    double* restrict const vgl1 = &(mo_vgl[ipoint*5*mo_num]);
+    const double* restrict avgl1 = &(ao_vgl[ipoint*5*ao_num]);
+    double* restrict const vgl2 =  vgl1 + mo_num;
+    double* restrict const vgl3 =  vgl1 + (mo_num << 1);
+    double* restrict const vgl4 =  vgl1 + (mo_num << 1) + mo_num;
+    double* restrict const vgl5 =  vgl1 + (mo_num << 2);
+    const double* restrict avgl2 = avgl1 + ao_num;
+    const double* restrict avgl3 = avgl1 + (ao_num << 1);
+    const double* restrict avgl4 = avgl1 + (ao_num << 1) + ao_num;
+    const double* restrict avgl5 = avgl1 + (ao_num << 2);
+
+    for (int64_t i=0 ; i<mo_num ; ++i) {
+      vgl1[i] = 0.;
+      vgl2[i] = 0.;
+      vgl3[i] = 0.;
+      vgl4[i] = 0.;
+      vgl5[i] = 0.;
+    }
+
+    int64_t nidx=0;
+    int64_t idx[ao_num];
+    double  av1[ao_num];
+    double  av2[ao_num];
+    double  av3[ao_num];
+    double  av4[ao_num];
+    double  av5[ao_num];
+    for (int64_t k=0 ; k<ao_num ; ++k) {
+      const double* restrict ck1 = coef_normalized_t + k*mo_num;
+      if (avgl1[k] != 0.) {
+        idx[nidx] = k;
+        av1[nidx] = avgl1[k];
+        av2[nidx] = avgl2[k];
+        av3[nidx] = avgl3[k];
+        av4[nidx] = avgl4[k];
+        av5[nidx] = avgl5[k];
+        ++nidx;
+      }
+    }
+
+    int64_t n;
+    for (n=0 ; n < nidx-4 ; n+=4) {
+      int64_t k = idx[n];
+      const double* restrict ck1 = coef_normalized_t + idx[n  ]*mo_num;
+      const double* restrict ck2 = coef_normalized_t + idx[n+1]*mo_num;
+      const double* restrict ck3 = coef_normalized_t + idx[n+2]*mo_num;
+      const double* restrict ck4 = coef_normalized_t + idx[n+3]*mo_num;
+
+      const double a11 = av1[n  ];
+      const double a21 = av1[n+1];
+      const double a31 = av1[n+2];
+      const double a41 = av1[n+3];
+      
+      const double a12 = av2[n  ];
+      const double a22 = av2[n+1];
+      const double a32 = av2[n+2];
+      const double a42 = av2[n+3];
+      
+      const double a13 = av3[n  ];
+      const double a23 = av3[n+1];
+      const double a33 = av3[n+2];
+      const double a43 = av3[n+3];
+      
+      const double a14 = av4[n  ];
+      const double a24 = av4[n+1];
+      const double a34 = av4[n+2];
+      const double a44 = av4[n+3];
+      
+      const double a15 = av5[n  ];
+      const double a25 = av5[n+1];
+      const double a35 = av5[n+2];
+      const double a45 = av5[n+3];
+
+#ifdef HAVE_OPENMP
+#pragma omp simd
+#endif
+      for (int64_t i=0 ; i<mo_num ; ++i) {
+        vgl1[i] = vgl1[i] + ck1[i] * a11 + ck2[i] * a21 + ck3[i] * a31 + ck4[i] * a41;
+        vgl2[i] = vgl2[i] + ck1[i] * a12 + ck2[i] * a22 + ck3[i] * a32 + ck4[i] * a42;
+        vgl3[i] = vgl3[i] + ck1[i] * a13 + ck2[i] * a23 + ck3[i] * a33 + ck4[i] * a43;
+        vgl4[i] = vgl4[i] + ck1[i] * a14 + ck2[i] * a24 + ck3[i] * a34 + ck4[i] * a44;
+        vgl5[i] = vgl5[i] + ck1[i] * a15 + ck2[i] * a25 + ck3[i] * a35 + ck4[i] * a45;
+      }
+    }
+
+    int64_t n0 = nidx-4;
+    n0 = n0 < 0 ? 0 : n0;
+    for (int64_t n=n0 ; n < nidx ; n+=1) {
+      const double* restrict ck = coef_normalized_t + idx[n]*mo_num;
+      const double a1 = av1[n];
+      const double a2 = av2[n];
+      const double a3 = av3[n];
+      const double a4 = av4[n];
+      const double a5 = av5[n];
+
+#ifdef HAVE_OPENMP
+  #pragma omp simd
+#endif
+      for (int64_t i=0 ; i<mo_num ; ++i) {
+        vgl1[i] += ck[i] * a1;
+        vgl2[i] += ck[i] * a2;
+        vgl3[i] += ck[i] * a3;
+        vgl4[i] += ck[i] * a4;
+        vgl5[i] += ck[i] * a5;
+      }
+    }
+  }
+  return QMCKL_SUCCESS;
+}
+#endif
     #+end_src
 
 *** Test
@@ -772,7 +1064,7 @@ assert (rc == QMCKL_SUCCESS);
 assert(qmckl_mo_basis_provided(context));
 
 double mo_vgl[walk_num*elec_num][5][chbrclf_mo_num];
-rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0]));
+rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), walk_num*elec_num*5*chbrclf_mo_num);
 assert (rc == QMCKL_SUCCESS);
 
 // Test overlap of MO
@@ -807,7 +1099,7 @@ assert (rc == QMCKL_SUCCESS);
 //
 //      // Calculate value of MO (1st electron)
 //      double mo_vgl[5][walk_num][elec_num][chbrclf_mo_num];
-//      rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0][0]));
+//      rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0][0]));
 //      assert (rc == QMCKL_SUCCESS);
 //      ovlmo1 += mo_vgl[0][0][0][0]*mo_vgl[0][0][0][0]*dr3;
 //    }
diff --git a/org/qmckl_nucleus.org b/org/qmckl_nucleus.org
index 319d5d1..05907ce 100644
--- a/org/qmckl_nucleus.org
+++ b/org/qmckl_nucleus.org
@@ -125,7 +125,7 @@ qmckl_exit_code qmckl_init_nucleus(qmckl_context context) {
     return false;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   ctx->nucleus.uninitialized = (1 << 3) - 1;
@@ -167,7 +167,7 @@ qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num) {
                            "num is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 0;
@@ -226,7 +226,7 @@ qmckl_get_nucleus_charge (const qmckl_context context,
                            "charge is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 1;
@@ -293,7 +293,7 @@ qmckl_get_nucleus_rescale_factor (const qmckl_context context,
                            "rescale_factor_kappa is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   assert (ctx->nucleus.rescale_factor_kappa > 0.0);
@@ -351,7 +351,7 @@ qmckl_get_nucleus_coord (const qmckl_context context,
                            "coord is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int32_t mask = 1 << 2;
@@ -410,7 +410,7 @@ bool qmckl_nucleus_provided(const qmckl_context context) {
     return false;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   return ctx->nucleus.provided;
@@ -425,7 +425,7 @@ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
   return QMCKL_NULL_CONTEXT;
  }
 
-qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
    #+end_src
 
    #+NAME:post2
@@ -672,7 +672,6 @@ end interface
 ** Test
 
     #+begin_src c :tangle (eval c_test)
-const int64_t   nucl_num      = chbrclf_nucl_num;
 const double*   nucl_charge   = chbrclf_charge;
 const double*   nucl_coord    = &(chbrclf_nucl_coord[0][0]);
 const double    nucl_rescale_factor_kappa = 2.0;
@@ -688,13 +687,13 @@ rc = qmckl_get_nucleus_num (context, &n);
 assert(rc == QMCKL_NOT_PROVIDED);
 
 
-rc = qmckl_set_nucleus_num (context, nucl_num);
+rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num);
 assert(rc == QMCKL_SUCCESS);
 assert(!qmckl_nucleus_provided(context));
 
 rc = qmckl_get_nucleus_num (context, &n);
 assert(rc == QMCKL_SUCCESS);
-assert(n == nucl_num);
+assert(n == chbrclf_nucl_num);
 
 double k;
 rc = qmckl_get_nucleus_rescale_factor (context, &k);
@@ -709,41 +708,41 @@ rc = qmckl_get_nucleus_rescale_factor (context, &k);
 assert(rc == QMCKL_SUCCESS);
 assert(k == nucl_rescale_factor_kappa);
 
-double nucl_coord2[3*nucl_num];
+double nucl_coord2[3*chbrclf_nucl_num];
 
-rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num);
+rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*chbrclf_nucl_num);
 assert(rc == QMCKL_NOT_PROVIDED);
 
-rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num);
+rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*chbrclf_nucl_num);
 assert(rc == QMCKL_SUCCESS);
 
 assert(!qmckl_nucleus_provided(context));
 
-rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, 3*nucl_num);
+rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, 3*chbrclf_nucl_num);
 assert(rc == QMCKL_SUCCESS);
 for (size_t k=0 ; k<3 ; ++k) {
-  for (int64_t i=0 ; i<nucl_num ; ++i) {
-    assert( nucl_coord[nucl_num*k+i] == nucl_coord2[3*i+k] );
+  for (int64_t i=0 ; i<chbrclf_nucl_num ; ++i) {
+    assert( nucl_coord[chbrclf_nucl_num*k+i] == nucl_coord2[3*i+k] );
   }
 }
 
-rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num);
+rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*chbrclf_nucl_num);
 assert(rc == QMCKL_SUCCESS);
-for (int64_t i=0 ; i<3*nucl_num ; ++i) {
+for (int64_t i=0 ; i<3*chbrclf_nucl_num ; ++i) {
   assert( nucl_coord[i] == nucl_coord2[i] );
 }
 
-double nucl_charge2[nucl_num];
+double nucl_charge2[chbrclf_nucl_num];
 
-rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num);
+rc = qmckl_get_nucleus_charge(context, nucl_charge2, chbrclf_nucl_num);
 assert(rc == QMCKL_NOT_PROVIDED);
 
-rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num);
+rc = qmckl_set_nucleus_charge(context, nucl_charge, chbrclf_nucl_num);
 assert(rc == QMCKL_SUCCESS);
 
-rc = qmckl_get_nucleus_charge(context, nucl_charge2, nucl_num);
+rc = qmckl_get_nucleus_charge(context, nucl_charge2, chbrclf_nucl_num);
 assert(rc == QMCKL_SUCCESS);
-for (int64_t i=0 ; i<nucl_num ; ++i) {
+for (int64_t i=0 ; i<chbrclf_nucl_num ; ++i) {
   assert( nucl_charge[i] == nucl_charge2[i] );
  }
 assert(qmckl_nucleus_provided(context));
@@ -784,7 +783,7 @@ qmckl_get_nucleus_nn_distance(qmckl_context context,
   qmckl_exit_code rc = qmckl_provide_nn_distance(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   const int64_t sze = ctx->nucleus.num * ctx->nucleus.num;
@@ -828,7 +827,7 @@ qmckl_exit_code qmckl_provide_nn_distance(qmckl_context context)
     return (char) 0;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED;
@@ -940,10 +939,10 @@ qmckl_exit_code qmckl_compute_nn_distance (
 
 assert(qmckl_nucleus_provided(context));
 
-double distance[nucl_num*nucl_num];
-rc = qmckl_get_nucleus_nn_distance(context, distance, nucl_num*nucl_num);
+double distance[chbrclf_nucl_num*chbrclf_nucl_num];
+rc = qmckl_get_nucleus_nn_distance(context, distance, chbrclf_nucl_num*chbrclf_nucl_num);
 assert(distance[0] == 0.);
-assert(distance[1] == distance[nucl_num]);
+assert(distance[1] == distance[chbrclf_nucl_num]);
 assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
 
      #+end_src
@@ -973,7 +972,7 @@ qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context,
   qmckl_exit_code rc = qmckl_provide_nn_distance_rescaled(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   const int64_t sze = ctx->nucleus.num * ctx->nucleus.num;
@@ -1019,7 +1018,7 @@ qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context)
     return (char) 0;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED;
@@ -1167,7 +1166,7 @@ qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* const
   qmckl_exit_code rc = qmckl_provide_nucleus_repulsion(context);
   if (rc != QMCKL_SUCCESS) return rc;
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   *energy = ctx->nucleus.repulsion;
@@ -1203,7 +1202,7 @@ qmckl_exit_code qmckl_provide_nucleus_repulsion(qmckl_context context)
     return (char) 0;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   qmckl_exit_code rc;
diff --git a/org/qmckl_numprec.org b/org/qmckl_numprec.org
index 3b73dc3..f7b0edf 100644
--- a/org/qmckl_numprec.org
+++ b/org/qmckl_numprec.org
@@ -141,7 +141,7 @@ qmckl_exit_code qmckl_set_numprec_precision(const qmckl_context context, const i
                           "precision > 53");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
 
   /* This should be always true because the context is valid */
   assert (ctx != NULL);
@@ -185,7 +185,7 @@ int qmckl_get_numprec_precision(const qmckl_context context) {
                       "");
   }
 
-  const qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  const qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   return ctx->numprec.precision;
 }
   #+end_src
@@ -232,7 +232,7 @@ qmckl_exit_code qmckl_set_numprec_range(const qmckl_context context, const int r
                     "range > 11");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
 
   /* This should be always true because the context is valid */
   assert (ctx != NULL);
@@ -275,7 +275,7 @@ int qmckl_get_numprec_range(const qmckl_context context) {
                       "");
   }
 
-  const qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  const qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   return ctx->numprec.range;
 }
   #+end_src
diff --git a/org/qmckl_point.org b/org/qmckl_point.org
index 07b0863..74a1e6a 100644
--- a/org/qmckl_point.org
+++ b/org/qmckl_point.org
@@ -81,7 +81,7 @@ int main() {
   |----------+----------------+-------------------------------------------|
   | ~num~    | ~int64_t~      | Total number of points                    |
   | ~date~   | ~uint64_t~     | Last modification date of the coordinates |
-  | ~coord~  | ~qmckl_matrix~ | ~num~ \times 3 matrix3                    |
+  | ~coord~  | ~qmckl_matrix~ | ~num~ \times 3 matrix                     |
 
   We consider that the matrix is stored 'transposed' and 'normal'
   corresponds to the 3 \times ~num~ matrix.
@@ -108,7 +108,7 @@ qmckl_exit_code qmckl_init_point(qmckl_context context) {
     return false;
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   memset(&(ctx->point), 0, sizeof(qmckl_point_struct));
@@ -148,7 +148,7 @@ qmckl_get_point_num (const qmckl_context context, int64_t* const num) {
                            "num is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   assert (ctx->point.num >= (int64_t) 0);
@@ -202,7 +202,7 @@ qmckl_get_point(const qmckl_context context,
                            "coord is a null pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   int64_t point_num = ctx->point.num;
@@ -295,7 +295,7 @@ qmckl_set_point (qmckl_context context,
                            "coord is a NULL pointer");
   }
 
-  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);
 
   qmckl_exit_code rc;