diff --git a/.github/workflows/test-build.yml b/.github/workflows/test-build.yml
index 41b59e0..c7af73f 100644
--- a/.github/workflows/test-build.yml
+++ b/.github/workflows/test-build.yml
@@ -47,6 +47,7 @@ jobs:
 
     - name: Run test
       run: make -j 4 check
+      working-directory: _build
 
     - name: Archive test log file
       if: failure()
@@ -57,13 +58,7 @@ jobs:
 
     - name: Dist test
       run: make distcheck
-
-    - name: Archive test log file
-      if: failure()
-      uses: actions/upload-artifact@v2
-      with:
-        name: dist-report-ubuntu
-        path: test-suite.log
+      working-directory: _build
 
   x86_macos:
 
@@ -91,37 +86,26 @@ jobs:
         git clone https://github.com/TREX-CoE/trexio.git
         cd trexio
         ./autogen.sh
-        mkdir _build
-        cd _build
-        ../configure --enable-silent-rules
-        make -j 4
-
-    - name: Run test
-      run: make -j 4 check
-
-    - name: Archive test log file
-      if: failure()
-      uses: actions/upload-artifact@v2
-      with:
-        name: test-report-ubuntu
-        path: test-suite.log
-
-    - name: Dist test
-      run: make distcheck
-
-    - name: Archive test log file
-      if: failure()
-      uses: actions/upload-artifact@v2
-      with:
         ./configure --prefix=${PWD}/_install --enable-silent-rules
         make -j 4
         make install
 
+    - name: Test TREXIO
+      run: make -j 4 check
+      working-directory: trexio
+
+    - name: Archive TREXIO test log file
+      if: failure()
+      uses: actions/upload-artifact@v2
+      with:
+        name: test-report-trexio-macos
+        path: trexio/test-suite.log
+
     - name: Build QMCkl
       run: |
         export PKG_CONFIG_PATH=${PWD}/trexio/_install/lib/pkgconfig:$PKG_CONFIG_PATH
         ./autogen.sh
-        ./configure --enable-silent-rules --enable-debug
+        ./configure CC=gcc-10 FC=gfortran-10 --enable-silent-rules --enable-debug
         make -j 4
 
     - name: Run test
diff --git a/README.md b/README.md
index 4ab5bae..703bfee 100644
--- a/README.md
+++ b/README.md
@@ -73,6 +73,24 @@ sudo make install
 sudo make installcheck
 ```
 
+
+## Linking to your program
+
+The `make install` command takes care of installing the QMCkl shared library on the user machine.
+Once installed, add `-lqmckl` to the list of compiler options.
+
+In some cases (e.g. when using custom `prefix` during configuration), the QMCkl library might end up installed in a directory, which is absent in the default `$LIBRARY_PATH`.
+In order to link the program against QMCkl, the search paths can be modified as follows:
+
+`export LIBRARY_PATH=$LIBRARY_PATH:<path_to_qmckl>/lib`
+
+(same holds for `$LD_LIBRARY_PATH`). The `<path_to_qmckl>` has to be replaced with the prefix used during the installation.
+
+If your project relies on the CMake build system, feel free to use the
+[FindQMCKL.cmake](https://github.com/TREX-CoE/qmckl/blob/master/cmake/FindQMCKL.cmake)
+module to find and link the QMCkl library automatically.
+
+
 ## Verificarlo CI
 
 Since Verificarlo should not be a dependency of QMCkl, all Verificarlo
diff --git a/cmake/FindQMCKL.cmake b/cmake/FindQMCKL.cmake
new file mode 100644
index 0000000..89b0b11
--- /dev/null
+++ b/cmake/FindQMCKL.cmake
@@ -0,0 +1,77 @@
+#===========================================
+
+# Try to find the QMCkl library;
+# If found, it will define the following variables (note the plural form):
+#  QMCKL_FOUND - System has libqmckl;
+#  QMCKL_INCLUDE_DIRS - The QMCKL include directories;
+#  QMCKL_LIBRARIES - The libraries needed to use QMCKL;
+
+# If QMCKL has been installed in a non-standard location, one can set an 
+# environment variable $QMCKL_DIR in the current shell:
+# $ export QMCKL_DIR=<custom_path>
+# to indicate the prefix used during the QMCKL installation
+# (typically `./configure prefix=<custom_path> ..` or `cmake -DCMAKE_INSTALL_DIR=<custom_path> ..`)
+
+# This file should be located WITHIN your project source tree.
+# (e.g. in cmake/FindQMCKL.cmake)
+# How to use it in your project CMakeLists.txt:
+
+# This is needed to locate FindQMCKL.cmake file, modify it according to your source tree.
+# list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake/")
+
+# find_package(QMCKL)
+# if (QMCKL_FOUND)
+#   include_directories(${QMCKL_INCLUDE_DIRS})
+#   target_link_libraries(your_target ${QMCKL_LIBRARIES})
+# endif()
+
+#===========================================
+
+# This file is distirbuted under the BSD 3-Clause License.
+# Copyright (c) 2021, TREX Center of Excellence
+
+#===========================================
+
+message("<FindQMCKL.cmake>")
+
+set(QMCKL_SEARCH_PATHS
+	~/Library/Frameworks
+	/Library/Frameworks
+	/usr/local
+	/usr
+	/sw # Fink
+	/opt/local # DarwinPorts
+	/opt/csw # Blastwave
+	/opt
+)
+
+find_path(QMCKL_INCLUDE_DIR 
+	  NAMES qmckl.h 
+	  HINTS $ENV{QMCKL_DIR}
+	  PATH_SUFFIXES include/qmckl include
+	  PATHS ${QMCKL_SEARCH_PATHS}
+	  )
+	        
+
+# No need to specify platform-specific prefix (e.g. libqmckl on Unix) or
+# suffix (e.g. .so on Unix or .dylib on MacOS) in NAMES. CMake takes care of that.
+find_library(QMCKL_LIBRARY
+             NAMES qmckl
+	     HINTS $ENV{QMCKL_DIR}
+	     PATH_SUFFIXES lib64 lib
+	     PATHS ${QMCKL_SEARCH_PATHS}
+	     )
+
+message("<FindQMCKL.cmake>")
+
+# Handle the QUIETLY and REQUIRED arguments and set QMCKL_FOUND to TRUE if
+# all listed variables are TRUE.
+INCLUDE(FindPackageHandleStandardArgs)
+FIND_PACKAGE_HANDLE_STANDARD_ARGS(QMCKL DEFAULT_MSG QMCKL_LIBRARY QMCKL_INCLUDE_DIR )
+MARK_AS_ADVANCED(QMCKL_INCLUDE_DIR QMCKL_LIBRARY)
+
+# Mot setting _INCLUDE_DIR and _LIBRARIES is considered a bug, 
+# see https://gitlab.kitware.com/cmake/community/-/wikis/doc/tutorials/How-To-Find-Libraries
+set(QMCKL_LIBRARIES ${QMCKL_LIBRARY})
+set(QMCKL_INCLUDE_DIRS ${QMCKL_INCLUDE_DIR}) 
+
diff --git a/configure.ac b/configure.ac
index fcc1dc5..968fbf8 100644
--- a/configure.ac
+++ b/configure.ac
@@ -49,12 +49,12 @@ AS_IF([test -d ${srcdir}/.git],
 AC_ARG_WITH(ifort, [AS_HELP_STRING([--with-ifort],[Use Intel Fortran compiler])], with_ifort=$withval, with_ifort=no)
 AS_IF([test "$with_ifort" == "yes"], [
   FC=ifort
-  FCFLAGS="-xHost -ip -O2 -ftz -finline -g -mkl=sequential" ])
+  FCFLAGS="-xHost -ip -Ofast -ftz -finline -g -mkl=sequential" ])
 
 AC_ARG_WITH(icc, [AS_HELP_STRING([--with-icc],[Use Intel C compiler])], with_icc=$withval, with_icc=no)
 AS_IF([test "$with_icc" == "yes"], [
   CC=icc
-  CFLAGS="-xHost -ip -O2 -ftz -finline -g -mkl=sequential" ])
+  CFLAGS="-xHost -ip -Ofast -ftz -finline -g -mkl=sequential" ])
 
 AS_IF([test "$with_icc"."$with_ifort" == "yes.yes"], [
   ax_blas_ok="yes"
diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org
index 55057bd..7cc7eef 100644
--- a/org/qmckl_ao.org
+++ b/org/qmckl_ao.org
@@ -1,3 +1,4 @@
+B
 #+TITLE: Atomic Orbitals
 #+SETUPFILE: ../tools/theme.setup
 #+INCLUDE: ../tools/lib.org
@@ -268,7 +269,7 @@ end interface
 
 *** Data structure                                                 :noexport:
 
-    #+begin_src c :comments org :tangle (eval h_private_type)
+    #+begin_src c :comments org :tangle (eval h_private_type) :noweb yes
 typedef struct qmckl_ao_basis_struct {
   int64_t   shell_num;
   int64_t   prim_num;
@@ -299,6 +300,9 @@ typedef struct qmckl_ao_basis_struct {
   bool      provided;
   bool      ao_cartesian;
   char      type;
+#ifdef HAVE_HPC
+  <<HPC_struct>>
+#endif
 } qmckl_ao_basis_struct;
     #+end_src
 
@@ -770,14 +774,14 @@ qmckl_set_ao_basis_shell_prim_num (qmckl_context context,
 }
      #+end_src
 
-
+     
      #+begin_src c :comments org :tangle (eval h_func)
 qmckl_exit_code
 qmckl_set_ao_basis_shell_prim_index (qmckl_context context,
                                      const int64_t* shell_prim_index,
                                      const int64_t size_max);
      #+end_src
-
+     
      #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
 qmckl_exit_code
 qmckl_set_ao_basis_shell_prim_index (qmckl_context context,
@@ -831,15 +835,15 @@ qmckl_set_ao_basis_shell_prim_index (qmckl_context context,
   <<post2>>
 }
      #+end_src
-
-
+     
+     
      #+begin_src c :comments org :tangle (eval h_func)
 qmckl_exit_code
 qmckl_set_ao_basis_shell_factor (qmckl_context context,
                                  const double* shell_factor,
                                  const int64_t size_max);
      #+end_src
-
+     
      #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
 qmckl_exit_code
 qmckl_set_ao_basis_shell_factor (qmckl_context context,
@@ -893,15 +897,15 @@ qmckl_set_ao_basis_shell_factor (qmckl_context context,
   <<post2>>
 }
      #+end_src
-
-
+     
+     
      #+begin_src c :comments org :tangle (eval h_func)
 qmckl_exit_code
 qmckl_set_ao_basis_exponent (qmckl_context context,
                              const double* exponent,
                              const int64_t size_max);
      #+end_src
-
+     
      #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
 qmckl_exit_code
 qmckl_set_ao_basis_exponent (qmckl_context context,
@@ -2485,7 +2489,6 @@ free(ao_factor_test);
    | ~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
 
     When the basis set is completely entered, extra data structures may be
@@ -2625,11 +2628,114 @@ qmckl_exit_code qmckl_finalize_basis(qmckl_context context) {
       }
     }
   }
-  /* TODO : sort the basis set here */
-  return QMCKL_SUCCESS;
+
+  rc = QMCKL_SUCCESS;
+#ifdef HAVE_HPC
+  rc = qmckl_finalize_basis_hpc(context);
+#endif
+
+  return rc;
 }
       #+end_src
 
+*** HPC-specific data structures
+
+    For faster access, we provide extra arrays for the shell information as:
+    - $C_{psa}$ = =coef_per_nucleus (iprim, ishell, inucl)=
+    - $\gamma_{pa}$ =expo_per_nucleus (iprim, inucl)=
+
+    such that the computation of the radial parts can be vectorized.
+    
+    Exponents are sorted in increasing order to exit loops quickly,
+    and only unique exponents are kept. This also allows to pack the
+    exponents to enable vectorization of exponentials.
+
+    The computation of the radial part is made as
+    \[
+      R_{sa} = \sum_p C_{psa} \gamma_{pa}
+    \]
+    which is a matrix-vector product.
+    
+    #+NAME:HPC_struct 
+    #+begin_src c :comments org :exports none
+ /* HPC specific data structures */
+ qmckl_tensor coef_per_nucleus;
+ qmckl_matrix expo_per_nucleus;
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval h_private_func) :exports none
+#ifdef HAVE_HPC
+qmckl_exit_code qmckl_finalize_basis_hpc (qmckl_context context);
+#endif
+    #+end_src
+    
+    #+begin_src c :comments org :tangle (eval c) :exports none
+#ifdef HAVE_HPC
+qmckl_exit_code qmckl_finalize_basis_hpc (qmckl_context context)
+{
+  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+
+  int64_t shell_max = 0;
+  int64_t prim_max  = 0;
+  int64_t nucl_num = ctx->nucleus.num;
+  for (int inucl=0 ; inucl < nucl_num ; ++inucl) {
+    shell_max = ctx->ao_basis.nucleus_shell_num[inucl] > shell_max ?
+                ctx->ao_basis.nucleus_shell_num[inucl] : shell_max;
+
+    int64_t prim_num = 0;
+    const int64_t ishell_start = ctx->ao_basis.nucleus_index[inucl];
+    const int64_t ishell_end   = ctx->ao_basis.nucleus_index[inucl] + ctx->ao_basis.nucleus_shell_num[inucl];
+    for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) {
+      prim_num += ctx->ao_basis.shell_prim_num[ishell];
+    }
+    prim_max = prim_num > prim_max ?
+               prim_num : prim_max;
+  }
+
+  int64_t size[3] = { prim_max, shell_max, nucl_num };
+  ctx->ao_basis.coef_per_nucleus = qmckl_tensor_alloc( context, 3, size );
+  ctx->ao_basis.coef_per_nucleus = qmckl_tensor_set(ctx->ao_basis.coef_per_nucleus, 0.);
+
+  ctx->ao_basis.expo_per_nucleus = qmckl_matrix_alloc( context, prim_max, nucl_num );
+  ctx->ao_basis.expo_per_nucleus = qmckl_matrix_set(ctx->ao_basis.expo_per_nucleus, 0.);
+
+  double expo[prim_max];
+  double coef[shell_max][prim_max];
+  for (int inucl=0 ; inucl < nucl_num ; ++inucl) {
+    memset(expo, 0, sizeof(expo));
+    memset(coef, 0, sizeof(expo));
+
+    int64_t idx = 0;
+    const int64_t ishell_start = ctx->ao_basis.nucleus_index[inucl];
+    const int64_t ishell_end   = ctx->ao_basis.nucleus_index[inucl] + ctx->ao_basis.nucleus_shell_num[inucl];
+    for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) {
+      
+      const int64_t iprim_start = ctx->ao_basis.shell_prim_index[ishell];
+      const int64_t iprim_end   = ctx->ao_basis.shell_prim_index[ishell] + ctx->ao_basis.shell_prim_num[ishell];
+      for (int64_t iprim = iprim_start ; iprim < iprim_end ; ++iprim) {
+        expo[idx] = ctx->ao_basis.coefficient_normalized[iprim];
+        coef[ishell - ishell_start][idx] = ctx->ao_basis.coefficient_normalized[iprim];
+        idx += 1;
+      }
+    }
+
+    for (int64_t i=0 ; i<prim_max ; ++i) {
+      qmckl_mat( ctx->ao_basis.expo_per_nucleus, i, inucl ) = expo[i];
+    }
+
+    for (int64_t j=0 ; j<shell_max ; ++j) {
+      for (int64_t i=0 ; i<prim_max ; ++i) {
+        qmckl_ten3( ctx->ao_basis.coef_per_nucleus, i, j, inucl ) = coef[j][i];
+      }
+    }
+
+  }
+  return QMCKL_SUCCESS;
+}
+#endif
+  
+    #+end_src
+
 *** Access functions
 
     #+begin_src c :comments org :tangle (eval h_func) :noweb yes
@@ -3040,11 +3146,11 @@ integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C)
   integer(c_int64_t), intent(in), value :: context
 
   integer*8                     :: n, ldv, j, i
-  double precision              :: X(3), R(3), Y(3), r2
+  double precision              :: X(3), R(3), Y(3), r2, z
   double precision, allocatable :: VGL(:,:), A(:)
   double precision              :: epsilon
 
-  epsilon = qmckl_get_numprec_epsilon(context)
+  epsilon = 3.d0 * qmckl_get_numprec_epsilon(context)
 
   X = (/ 1.1 , 2.2 ,  3.3 /)
   R = (/ 0.1 , 1.2 , -2.3 /)
@@ -3068,29 +3174,43 @@ integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C)
 
   do i=1,n
      test_qmckl_ao_gaussian_vgl = -11
-     if (dabs(1.d0 - VGL(i,1) / (&
-          dexp(-A(i) * r2) &
-          )) > epsilon ) return
+     z = dabs(1.d0 - VGL(i,1) / (dexp(-A(i) * r2)) )
+     if ( z > epsilon ) then
+        print *, z, epsilon
+        return
+     end if
 
      test_qmckl_ao_gaussian_vgl = -12
-     if (dabs(1.d0 - VGL(i,2) / (&
-          -2.d0 * A(i) * Y(1) * dexp(-A(i) * r2) &
-          )) > epsilon ) return
+     z = dabs(1.d0 - VGL(i,2) / (&
+          -2.d0 * A(i) * Y(1) * dexp(-A(i) * r2) ))
+     if ( z > epsilon ) then
+        print *, z, epsilon
+        return
+     end if
 
      test_qmckl_ao_gaussian_vgl = -13
-     if (dabs(1.d0 - VGL(i,3) / (&
-          -2.d0 * A(i) * Y(2) * dexp(-A(i) * r2) &
-          )) > epsilon ) return
+     z = dabs(1.d0 - VGL(i,3) / (&
+          -2.d0 * A(i) * Y(2) * dexp(-A(i) * r2) ))
+     if ( z > epsilon ) then
+        print *, z, epsilon
+        return
+     end if
 
      test_qmckl_ao_gaussian_vgl = -14
-     if (dabs(1.d0 - VGL(i,4) / (&
-          -2.d0 * A(i) * Y(3) * dexp(-A(i) * r2) &
-          )) > epsilon ) return
+     z = dabs(1.d0 - VGL(i,4) / (&
+          -2.d0 * A(i) * Y(3) * dexp(-A(i) * r2) ))
+     if ( z > epsilon ) then
+        print *, z, epsilon
+        return
+     end if
 
      test_qmckl_ao_gaussian_vgl = -15
-     if (dabs(1.d0 - VGL(i,5) / (&
-          A(i) * (4.d0*r2*A(i) - 6.d0) * dexp(-A(i) * r2) &
-          )) > epsilon ) return
+     z = dabs(1.d0 - VGL(i,5) / (&
+          A(i) * (4.d0*r2*A(i) - 6.d0) * dexp(-A(i) * r2) ))
+     if ( z > epsilon ) then
+        print *, z, epsilon
+        return
+     end if
   end do
 
   test_qmckl_ao_gaussian_vgl = 0
@@ -5062,8 +5182,8 @@ for (int32_t ldl=3 ; ldl<=5 ; ++ldl) {
 
       for (int32_t k=0 ; k<5 ; ++k) {
         for (int32_t l=0 ; l<n ; ++l) {
-          printf("VGL[%d][%d] = %e   %e\n", k, l, VGL0[k][l], VGL1[k][l]);
-          assert( VGL0[k][l] == VGL1[k][l] );
+          printf("VGL[%d][%d] = %e   %e  %e\n", k, l, VGL0[k][l], VGL1[k][l], VGL0[k][l]-VGL1[k][l]);
+          assert( fabs(1.-(VGL0[k][l]+1.e-100)/(VGL1[k][l]+1.e-100)) < 1.e-15 );
         }
       }
     }
@@ -5256,6 +5376,7 @@ end function qmckl_compute_ao_vgl_doc_f
    | ~ao_vgl~              | ~double[point_num][5][ao_num]~ | out    | Value, gradients and Laplacian of the AOs    |
 
     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+#ifdef HAVE_HPC
 qmckl_exit_code
 qmckl_compute_ao_vgl_hpc_gaussian (
           const qmckl_context context,
@@ -5423,66 +5544,40 @@ qmckl_compute_ao_vgl_hpc_gaussian (
              const double s4 = s6_*z;
              const double s5 = s5_;
 
+             const int64_t k = ao_index[ishell];
+             double* __restrict__ const ao_vgl_1 = ao_vgl + ipoint*5*ao_num + k;
 
              const int32_t l = shell_ang_mom[ishell];
              const int32_t n = lstart[l+1]-lstart[l];
-             const int64_t k = ao_index[ishell];
 
-             double* const ao_vgl_1 = &(ao_vgl[ipoint*5*ao_num+k]);
-             double* const ao_vgl_2 = &(ao_vgl_1[  ao_num]);
-             double* const ao_vgl_3 = &(ao_vgl_1[2*ao_num]);
-             double* const ao_vgl_4 = &(ao_vgl_1[3*ao_num]);
-             double* const ao_vgl_5 = &(ao_vgl_1[4*ao_num]);
+             double* __restrict__ const ao_vgl_2 = ao_vgl_1 + ao_num;
+             double* __restrict__ const ao_vgl_3 = ao_vgl_1 + (ao_num<<1);
+             double* __restrict__ const ao_vgl_4 = ao_vgl_1 + (ao_num<<1) + ao_num;
+             double* __restrict__ const ao_vgl_5 = ao_vgl_1 + (ao_num<<2);
 
-             double* poly_vgl_1;
-             double* poly_vgl_2;
-             double* poly_vgl_3;
-             double* poly_vgl_4;
-             double* poly_vgl_5;
+             double* __restrict__ poly_vgl_1;
+             double* __restrict__ poly_vgl_2;
+             double* __restrict__ poly_vgl_3;
+             double* __restrict__ poly_vgl_4;
+             double* __restrict__ poly_vgl_5;
              if (nidx > 0) {
-                 const double* f = &(ao_factor[k]);
+                 const double* __restrict__ f = ao_factor + k;
                  const int64_t idx = lstart[l];
 
                  switch (nucleus_max_ang_mom[inucl]) {
                  case 0:
-                   ao_vgl_1[0] = s1 * f[0];
-                   ao_vgl_2[0] = s2 * f[0];
-                   ao_vgl_3[0] = s3 * f[0];
-                   ao_vgl_4[0] = s4 * f[0];
-                   ao_vgl_5[0] = s5;
                    break;
                  case 1:
                    poly_vgl_1 = &(poly_vgl_l1[0][idx]);
                    poly_vgl_2 = &(poly_vgl_l1[1][idx]);
                    poly_vgl_3 = &(poly_vgl_l1[2][idx]);
                    poly_vgl_4 = &(poly_vgl_l1[3][idx]);
-                   for (int64_t il=0 ; il<n ; ++il) {
-                     ao_vgl_1[il] =  poly_vgl_1[il] * s1 * f[il];
-                     ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
-                     ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il];
-                     ao_vgl_4[il] = (poly_vgl_4[il] * s1 + poly_vgl_1[il] * s4) * f[il];
-                     ao_vgl_5[il] = (poly_vgl_1[il] * s5 +
-                                     2.0*(poly_vgl_2[il] * s2 +
-                                          poly_vgl_3[il] * s3 +
-                                          poly_vgl_4[il] * s4 )) * f[il];
-                   }
                    break;
                  case 2:
                    poly_vgl_1 = &(poly_vgl_l2[0][idx]);
                    poly_vgl_2 = &(poly_vgl_l2[1][idx]);
                    poly_vgl_3 = &(poly_vgl_l2[2][idx]);
                    poly_vgl_4 = &(poly_vgl_l2[3][idx]);
-                   poly_vgl_5 = &(poly_vgl_l2[4][idx]);
-                   for (int64_t il=0 ; il<n ; ++il) {
-                     ao_vgl_1[il] =  poly_vgl_1[il] * s1 * f[il];
-                     ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
-                     ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il];
-                     ao_vgl_4[il] = (poly_vgl_4[il] * s1 + poly_vgl_1[il] * s4) * f[il];
-                     ao_vgl_5[il] = (poly_vgl_5[il] * s1 + poly_vgl_1[il] * s5 +
-                                     2.0*(poly_vgl_2[il] * s2 +
-                                          poly_vgl_3[il] * s3 +
-                                          poly_vgl_4[il] * s4 )) * f[il];
-                   }
                    break;
                  default:
                    poly_vgl_1 = &(poly_vgl[0][idx]);
@@ -5490,19 +5585,61 @@ qmckl_compute_ao_vgl_hpc_gaussian (
                    poly_vgl_3 = &(poly_vgl[2][idx]);
                    poly_vgl_4 = &(poly_vgl[3][idx]);
                    poly_vgl_5 = &(poly_vgl[4][idx]);
-                   for (int64_t il=0 ; il<n ; ++il) {
-                     ao_vgl_1[il] =  poly_vgl_1[il] * s1 * f[il];
-                     ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
-                     ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il];
-                     ao_vgl_4[il] = (poly_vgl_4[il] * s1 + poly_vgl_1[il] * s4) * f[il];
-                     ao_vgl_5[il] = (poly_vgl_5[il] * s1 + poly_vgl_1[il] * s5 +
-                                     2.0*(poly_vgl_2[il] * s2 +
-                                          poly_vgl_3[il] * s3 +
-                                          poly_vgl_4[il] * s4 )) * f[il];
-                   }
-                   break;
                  }
-
+                 switch (n) {
+                 case(1):
+                     ao_vgl_1[0] = s1 * f[0];
+                     ao_vgl_2[0] = s2 * f[0];
+                     ao_vgl_3[0] = s3 * f[0];
+                     ao_vgl_4[0] = s4 * f[0];
+                     ao_vgl_5[0] = s5;
+                     break;
+                 case (3):
+#ifdef HAVE_OPENMP
+                     #pragma omp simd
+#endif
+                     for (int il=0 ; il<3 ; ++il) {
+                       ao_vgl_1[il] =  poly_vgl_1[il] * s1 * f[il];
+                       ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
+                       ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il];
+                       ao_vgl_4[il] = (poly_vgl_4[il] * s1 + poly_vgl_1[il] * s4) * f[il];
+                       ao_vgl_5[il] = (poly_vgl_1[il] * s5 +
+                                       2.0*(poly_vgl_2[il] * s2 +
+                                            poly_vgl_3[il] * s3 +
+                                            poly_vgl_4[il] * s4 )) * f[il];
+                     }
+                     break;
+                 case(5):
+#ifdef HAVE_OPENMP
+                     #pragma omp simd
+#endif
+                     for (int il=0 ; il<5 ; ++il) {
+                       ao_vgl_1[il] =  poly_vgl_1[il] * s1 * f[il];
+                       ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
+                       ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il];
+                       ao_vgl_4[il] = (poly_vgl_4[il] * s1 + poly_vgl_1[il] * s4) * f[il];
+                       ao_vgl_5[il] = (poly_vgl_1[il] * s5 +
+                                       2.0*(poly_vgl_2[il] * s2 +
+                                            poly_vgl_3[il] * s3 +
+                                            poly_vgl_4[il] * s4 )) * f[il];
+                     }
+                     break;
+                 default:
+#ifdef HAVE_OPENMP
+                     #pragma omp simd simdlen(8)
+#endif
+                     for (int il=0 ; il<n ; ++il) {
+                       ao_vgl_1[il] =  poly_vgl_1[il] * s1 * f[il];
+                       ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
+                       ao_vgl_3[il] = (poly_vgl_3[il] * s1 + poly_vgl_1[il] * s3) * f[il];
+                       ao_vgl_4[il] = (poly_vgl_4[il] * s1 + poly_vgl_1[il] * s4) * f[il];
+                       ao_vgl_5[il] = (poly_vgl_5[il] * s1 + poly_vgl_1[il] * s5 +
+                                       2.0*(poly_vgl_2[il] * s2 +
+                                             poly_vgl_3[il] * s3 +
+                                             poly_vgl_4[il] * s4 )) * f[il];
+                     }
+                     break;
+                 }
              } else {
                  for (int64_t il=0 ; il<n ; ++il) {
                      ao_vgl_1[il] = 0.0;
@@ -5519,6 +5656,7 @@ qmckl_compute_ao_vgl_hpc_gaussian (
   }
   return QMCKL_SUCCESS;
 }
+#endif
     #+end_src
 
 ** Interfaces
@@ -5545,6 +5683,7 @@ qmckl_compute_ao_vgl_hpc_gaussian (
           double* const ao_vgl );
     #+end_src
     #+begin_src c :tangle (eval h_private_func) :comments org
+#ifdef HAVE_HPC
     qmckl_exit_code qmckl_compute_ao_vgl_hpc_gaussian (
           const qmckl_context context,
           const int64_t ao_num,
@@ -5565,6 +5704,7 @@ qmckl_compute_ao_vgl_hpc_gaussian (
           const double* expo,
           const double* coef_normalized,
           double* const ao_vgl );
+#endif
     #+end_src
 
     #+CALL: generate_c_interface(table=qmckl_ao_vgl_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl_doc"))
diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org
index 274d269..73f9e46 100644
--- a/org/qmckl_blas.org
+++ b/org/qmckl_blas.org
@@ -46,7 +46,8 @@
 
   #+begin_src c :comments link :tangle (eval c_test) :noweb yes
 #include "qmckl.h"
-#include "assert.h"
+#include <stdio.h>
+#include <assert.h>
 #ifdef HAVE_CONFIG_H
 #include "config.h"
 #endif
@@ -60,8 +61,20 @@ int main() {
 
   #+end_src
 
-* Data types
+* -
+:PROPERTIES:
+:UNNUMBERED: t
+:END:
+
+Basic linear algebra data types and operations are described in this file.
+The data types are private, so that HPC implementations can use
+whatever data structures they prefer.
+
+These data types are expected to be used internally in QMCkl. They
+are not intended to be passed to external codes.
   
+* Data types
+
 ** Vector
    
   | Variable | Type      | Description             |
@@ -69,7 +82,7 @@ int main() {
   | ~size~   | ~int64_t~ | Dimension of the vector |
   | ~data~   | ~double*~ | Elements                |
 
-   #+begin_src c :comments org :tangle (eval h_private_type)
+   #+begin_src c :comments org :tangle (eval h_private_type) :exports none
 typedef struct qmckl_vector {
   int64_t size;
   double* data;
@@ -85,7 +98,7 @@ qmckl_vector_alloc( qmckl_context context,
 
    Allocates a new vector. If the allocation failed the size is zero.
 
-   #+begin_src c :comments org :tangle (eval c)
+   #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_vector
 qmckl_vector_alloc( qmckl_context context, 
                     const int64_t size)
@@ -114,7 +127,7 @@ qmckl_vector_free( qmckl_context context,
                    qmckl_vector* vector);
    #+end_src
 
-   #+begin_src c :comments org :tangle (eval c)
+   #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_exit_code
 qmckl_vector_free( qmckl_context context, 
                    qmckl_vector* vector)
@@ -145,7 +158,7 @@ qmckl_vector_free( qmckl_context context,
   The dimensions use Fortran ordering: two elements differing by one
   in the first dimension are consecutive in memory.
   
-   #+begin_src c :comments org :tangle (eval h_private_type)
+   #+begin_src c :comments org :tangle (eval h_private_type) :exports none
 typedef struct qmckl_matrix {
   int64_t size[2];
   double* data;
@@ -162,7 +175,7 @@ qmckl_matrix_alloc( qmckl_context context,
 
    Allocates a new matrix. If the allocation failed the sizes are zero.
 
-   #+begin_src c :comments org :tangle (eval c)
+   #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_matrix
 qmckl_matrix_alloc( qmckl_context context, 
                     const int64_t size1,
@@ -195,7 +208,7 @@ qmckl_matrix_free( qmckl_context context,
                    qmckl_matrix* matrix);
    #+end_src
 
-   #+begin_src c :comments org :tangle (eval c)
+   #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_exit_code
 qmckl_matrix_free( qmckl_context context, 
                    qmckl_matrix* matrix)
@@ -228,7 +241,7 @@ qmckl_matrix_free( qmckl_context context,
   The dimensions use Fortran ordering: two elements differing by one
   in the first dimension are consecutive in memory.
   
-   #+begin_src c :comments org :tangle (eval h_private_type)
+   #+begin_src c :comments org :tangle (eval h_private_type) :exports none
 #define QMCKL_TENSOR_ORDER_MAX 16
 
 typedef struct qmckl_tensor {
@@ -249,7 +262,7 @@ qmckl_tensor_alloc( qmckl_context context,
    Allocates memory for a tensor. If the allocation failed, the size
    is zero.
 
-   #+begin_src c :comments org :tangle (eval c)
+   #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_tensor
 qmckl_tensor_alloc( qmckl_context context, 
                     const int64_t  order,
@@ -285,11 +298,11 @@ qmckl_tensor_alloc( qmckl_context context,
    
    #+begin_src c :comments org :tangle (eval h_private_func)
 qmckl_exit_code
-qmckl_tensor_free( qmckl_context context, 
+qmckl_tensor_free (qmckl_context context, 
                    qmckl_tensor* tensor);
    #+end_src
 
-   #+begin_src c :comments org :tangle (eval c)
+   #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_exit_code
 qmckl_tensor_free( qmckl_context context, 
                    qmckl_tensor* tensor)
@@ -325,7 +338,7 @@ qmckl_matrix_of_vector(const qmckl_vector vector,
 
     Reshapes a vector into a matrix.
 
-    #+begin_src c :comments org :tangle (eval c)
+    #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_matrix
 qmckl_matrix_of_vector(const qmckl_vector vector,
                        const int64_t size1,
@@ -355,7 +368,7 @@ qmckl_tensor_of_vector(const qmckl_vector vector,
 
     Reshapes a vector into a tensor.
 
-    #+begin_src c :comments org :tangle (eval c)
+    #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_tensor
 qmckl_tensor_of_vector(const qmckl_vector vector,
                        const int64_t order,
@@ -385,7 +398,7 @@ qmckl_vector_of_matrix(const qmckl_matrix matrix);
 
     Reshapes a matrix into a vector.
 
-    #+begin_src c :comments org :tangle (eval c)
+    #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_vector
 qmckl_vector_of_matrix(const qmckl_matrix matrix)
 { 
@@ -409,7 +422,7 @@ qmckl_tensor_of_matrix(const qmckl_matrix matrix,
 
     Reshapes a matrix into a tensor.
 
-    #+begin_src c :comments org :tangle (eval c)
+    #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_tensor
 qmckl_tensor_of_matrix(const qmckl_matrix matrix,
                        const int64_t order,
@@ -439,7 +452,7 @@ qmckl_vector_of_tensor(const qmckl_tensor tensor);
 
     Reshapes a tensor into a vector.
 
-    #+begin_src c :comments org :tangle (eval c)
+    #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_vector
 qmckl_vector_of_tensor(const qmckl_tensor tensor)
 { 
@@ -468,7 +481,7 @@ qmckl_matrix_of_tensor(const qmckl_tensor tensor,
 
     Reshapes a tensor into a vector.
 
-    #+begin_src c :comments org :tangle (eval c)
+    #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_matrix
 qmckl_matrix_of_tensor(const qmckl_tensor tensor,
                        const int64_t size1,
@@ -493,15 +506,88 @@ qmckl_matrix_of_tensor(const qmckl_tensor tensor,
 
 ** Access macros
 
-    #+begin_src c :comments org :tangle (eval h_private_func)
+   Macros are provided to ease the access to vectors, matrices and
+   tensors. Matrices use column-major ordering, as in Fortran.
+
+    #+begin_src c :tangle no
+double qmckl_vec (qmckl_vector v, int64_t i);           //  v[i] 
+double qmckl_mat (qmckl_matrix m, int64_t i, int64_t j) // m[j][i]
+
+double qmckl_ten3(qmckl_tensor t, int64_t i, int64_t j, int64_t k); // t[k][j][i]
+double qmckl_ten4(qmckl_tensor t, int64_t i, int64_t j, int64_t k, int64_t l);
+double qmckl_ten5(qmckl_tensor t, int64_t i, int64_t j, int64_t k, int64_t l, int64_t m);
+...
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval h_private_func) :exports none
 #define qmckl_vec(v, i) v.data[i]
 #define qmckl_mat(m, i, j) m.data[(i) + (j)*m.size[0]]
 
-#define qmckl_ten3(t, i, j, k) t.data[(i) + m.size[0]*((j) + size[1]*(k))]
-#define qmckl_ten4(t, i, j, k, l) t.data[(i) + m.size[0]*((j) + size[1]*((k) + size[2]*(l)))]
-#define qmckl_ten5(t, i, j, k, l, m) t.data[(i) + m.size[0]*((j) + size[1]*((k) + size[2]*((l) + size[3]*(m))))]
+#define qmckl_ten3(t, i, j, k) t.data[(i) + t.size[0]*((j) + t.size[1]*(k))]
+#define qmckl_ten4(t, i, j, k, l) t.data[(i) + t.size[0]*((j) + t.size[1]*((k) + t.size[2]*(l)))]
+#define qmckl_ten5(t, i, j, k, l, m) t.data[(i) + t.size[0]*((j) + t.size[1]*((k) + t.size[2]*((l) + t.size[3]*(m))))]
     #+end_src
   
+    For example:
+    
+** Set all elements
+   
+*** Vector
+
+    #+begin_src c :comments org :tangle (eval h_private_func)
+qmckl_vector
+qmckl_vector_set(qmckl_vector vector, double value);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :exports none
+qmckl_vector
+qmckl_vector_set(qmckl_vector vector, double value)
+{
+  for (int64_t i=0 ; i<vector.size ; ++i) {
+    qmckl_vec(vector, i) = value;
+  }
+  return vector;
+}
+    #+end_src
+
+*** Matrix
+
+    #+begin_src c :comments org :tangle (eval h_private_func)
+qmckl_matrix
+qmckl_matrix_set(qmckl_matrix matrix, double value);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :exports none
+qmckl_matrix
+qmckl_matrix_set(qmckl_matrix matrix, double value)
+{
+  qmckl_vector vector = qmckl_vector_of_matrix(matrix);
+  for (int64_t i=0 ; i<vector.size ; ++i) {
+    qmckl_vec(vector, i) = value;
+  }
+  return qmckl_matrix_of_vector(vector, matrix.size[0], matrix.size[1]);
+}
+    #+end_src
+    
+*** Tensor
+
+    #+begin_src c :comments org :tangle (eval h_private_func)
+qmckl_tensor
+qmckl_tensor_set(qmckl_tensor tensor, double value);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :exports none
+qmckl_tensor
+qmckl_tensor_set(qmckl_tensor tensor, double value)
+{
+  qmckl_vector vector = qmckl_vector_of_tensor(tensor);
+  for (int64_t i=0 ; i<vector.size ; ++i) {
+    qmckl_vec(vector, i) = value;
+  }
+  return qmckl_tensor_of_vector(vector, tensor.order, tensor.size);
+}
+    #+end_src
+    
 ** Copy to/from to ~double*~
 
     #+begin_src c :comments org :tangle (eval h_private_func)
@@ -514,7 +600,7 @@ qmckl_double_of_vector(const qmckl_context context,
 
     Converts a vector to a ~double*~.
     
-    #+begin_src c :comments org :tangle (eval c)
+    #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_exit_code
 qmckl_double_of_vector(const qmckl_context context,
                        const qmckl_vector vector,
@@ -546,7 +632,7 @@ qmckl_double_of_matrix(const qmckl_context context,
 
     Converts a matrix to a ~double*~.
     
-    #+begin_src c :comments org :tangle (eval c)
+    #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_exit_code
 qmckl_double_of_matrix(const qmckl_context context,
                        const qmckl_matrix matrix,
@@ -569,7 +655,7 @@ qmckl_double_of_tensor(const qmckl_context context,
 
     Converts a tensor to a ~double*~.
     
-    #+begin_src c :comments org :tangle (eval c)
+    #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_exit_code
 qmckl_double_of_tensor(const qmckl_context context,
                        const qmckl_tensor tensor,
@@ -592,7 +678,7 @@ qmckl_vector_of_double(const qmckl_context context,
 
     Converts a ~double*~ to a vector.
     
-    #+begin_src c :comments org :tangle (eval c)
+    #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_exit_code
 qmckl_vector_of_double(const qmckl_context context,
                        const double* target,
@@ -636,9 +722,9 @@ qmckl_matrix_of_double(const qmckl_context context,
                        qmckl_matrix* matrix);
     #+end_src
 
-    Converts a matrix to a ~double*~.
+    Converts a ~double*~ to a matrix.
     
-    #+begin_src c :comments org :tangle (eval c)
+    #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_exit_code
 qmckl_matrix_of_double(const qmckl_context context,
                        const double* target,
@@ -662,9 +748,9 @@ qmckl_tensor_of_double(const qmckl_context context,
                        qmckl_tensor* tensor);
     #+end_src
 
-    Converts a matrix to a ~double*~.
+    Converts a ~double*~ to a tensor.
     
-    #+begin_src c :comments org :tangle (eval c)
+    #+begin_src c :comments org :tangle (eval c) :exports none
 qmckl_exit_code
 qmckl_tensor_of_double(const qmckl_context context,
                        const double* target,
@@ -679,12 +765,9 @@ qmckl_tensor_of_double(const qmckl_context context,
 }
     #+end_src
 
+** Tests                                                           :noexport:
 
-
-
-** Tests
-
-    #+begin_src c :comments link :tangle (eval c_test)
+    #+begin_src c :comments link :tangle (eval c_test) :exports none
 {
   int64_t m = 3;
   int64_t n = 4;
@@ -720,7 +803,7 @@ qmckl_tensor_of_double(const qmckl_context context,
 
 ** ~qmckl_dgemm~
 
-   Matrix multiplication:
+   Matrix multiplication with a BLAS interface:
 
    \[
    C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}
@@ -780,7 +863,7 @@ qmckl_tensor_of_double(const qmckl_context context,
           const int64_t ldc );
     #+end_src
 
-    #+begin_src f90 :tangle (eval f)
+    #+begin_src f90 :tangle (eval f) :exports none
 integer function qmckl_dgemm_f(context, TransA, TransB, &
      m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) &
      result(info)
@@ -975,14 +1058,14 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C)
 end function test_qmckl_dgemm
     #+end_src
 
-    #+begin_src c :comments link :tangle (eval c_test)
+    #+begin_src c :comments link :tangle (eval c_test) :exports none
 qmckl_exit_code test_qmckl_dgemm(qmckl_context context);
 assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));
     #+end_src
 
 ** ~qmckl_matmul~
 
-   Matrix multiplication:
+   Matrix multiplication using the =qmckl_matrix= data type:
 
    \[
    C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}
@@ -1017,7 +1100,7 @@ qmckl_matmul (const qmckl_context context,
               qmckl_matrix* const C ); 
     #+end_src
 
-    #+begin_src c :tangle (eval c) :comments org
+    #+begin_src c :tangle (eval c) :comments org :exports none
 qmckl_exit_code
 qmckl_matmul (const qmckl_context context,
               const char TransA,
@@ -1155,9 +1238,87 @@ qmckl_matmul (const qmckl_context context,
 }
     #+end_src
 
+*** Test                                                           :noexport:
 
-*** TODO Test                                                      :noexport:
+    #+begin_src python :exports none :results output
+import numpy as np
 
+A = np.array([[ 1.,  2.,  3.,  4. ], 
+              [ 5.,  6.,  7.,  8. ],
+              [ 9., 10., 11., 12. ]])
+
+B = np.array([[ 1., -2.,  3.,  4.,  5. ],
+              [ 5., -6.,  7.,  8.,  9. ],
+              [ 9., 10., 11., 12., 13. ],
+              [10., 11., 12., 15., 14. ]])
+
+C = 0.5 * A @ B
+print(A.T)
+print(B.T)
+print(C.T)
+    #+end_src
+
+    #+RESULTS:
+    #+begin_example
+    [[ 1.  5.  9.]
+     [ 2.  6. 10.]
+     [ 3.  7. 11.]
+     [ 4.  8. 12.]]
+    [[ 1.  5.  9. 10.]
+     [-2. -6. 10. 11.]
+     [ 3.  7. 11. 12.]
+     [ 4.  8. 12. 15.]
+     [ 5.  9. 13. 14.]]
+    [[ 39.  89. 139.]
+     [ 30.  56.  82.]
+     [ 49. 115. 181.]
+     [ 58. 136. 214.]
+     [ 59. 141. 223.]]
+    #+end_example
+
+    #+begin_src c :comments link :tangle (eval c_test) :exports none
+{
+  double a[12] = { 1.,  5.,  9.,
+                   2.,  6., 10.,
+                   3.,  7., 11.,
+                   4.,  8., 12. };
+  
+  double b[20] = { 1.,  5.,  9., 10.,
+                   -2., -6., 10., 11.,
+                   3.,  7., 11., 12.,
+                   4.,  8., 12., 15.,
+                   5.,  9., 13., 14. };
+
+  double c[15] = { 39.,  89., 139.,
+                   30.,  56.,  82.,
+                   49., 115., 181.,
+                   58., 136., 214.,
+                   59., 141., 223. };
+
+  double cnew[15];
+
+  qmckl_exit_code rc;
+  qmckl_matrix A = qmckl_matrix_alloc(context, 3, 4);
+  rc = qmckl_matrix_of_double(context, a, 12, &A);
+  assert(rc == QMCKL_SUCCESS);
+
+  qmckl_matrix B = qmckl_matrix_alloc(context, 4, 5);
+  rc = qmckl_matrix_of_double(context, b, 20, &B);
+  assert(rc == QMCKL_SUCCESS);
+
+  qmckl_matrix C = qmckl_matrix_alloc(context, 3, 5);
+  rc = qmckl_matmul(context, 'N', 'N', 0.5, A, B, 0., &C);
+  assert(rc == QMCKL_SUCCESS);
+
+  rc = qmckl_double_of_matrix(context, C, cnew, 15);
+  assert(rc == QMCKL_SUCCESS);
+
+  for (int i=0 ; i<15 ; ++i) {
+    printf("%f %f\n", cnew[i], c[i]);
+    assert (c[i] == cnew[i]);
+  }
+}  
+    #+end_src
 ** ~qmckl_adjugate~
 
    Given a matrix $\mathbf{A}$, the adjugate matrix
@@ -1208,7 +1369,7 @@ qmckl_matmul (const qmckl_context context,
    for performance motivations. For larger sizes, we rely on the
    LAPACK library.
 
-    #+begin_src f90 :tangle (eval f)
+    #+begin_src f90 :tangle (eval f) :exports none
 integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) &
      result(info)
   use qmckl
@@ -1263,7 +1424,7 @@ integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) &
 
 end function qmckl_adjugate_f
     #+end_src
-
+    
     #+begin_src f90 :tangle (eval f) :exports none
 subroutine adjugate2(A,LDA,B,LDB,na,det_l)
   implicit none
@@ -2047,11 +2208,11 @@ assert(QMCKL_SUCCESS == test_qmckl_adjugate(context));
 
    Transposes a matrix: $A^\dagger_{ji} = A_{ij}$.
 
-   | Variable | Type          | In/Out | Description       |
-   |----------+---------------+--------+-------------------|
-   | context  | qmckl_context | in     | Global state      |
-   | A        | qmckl_matrix  | in     | Input matrix      |
-   | At       | qmckl_matrix  | out    | Transposed matrix |
+   | Variable  | Type            | In/Out | Description       |
+   |-----------+-----------------+--------+-------------------|
+   | ~context~ | ~qmckl_context~ | in     | Global state      |
+   | ~A~       | ~qmckl_matrix~  | in     | Input matrix      |
+   | ~At~      | ~qmckl_matrix~  | out    | Transposed matrix |
    
     #+begin_src c :tangle (eval h_private_func) :comments org
 qmckl_exit_code
@@ -2061,7 +2222,7 @@ qmckl_transpose (qmckl_context context,
     #+end_src
 
 
-    #+begin_src c :tangle (eval c) :comments org
+    #+begin_src c :tangle (eval c) :comments org :exports none
 qmckl_exit_code
 qmckl_transpose (qmckl_context context,
                  const qmckl_matrix A,
@@ -2100,7 +2261,7 @@ qmckl_transpose (qmckl_context context,
 }
     #+end_src
 
-*** Test
+*** Test                                                           :noexport:
 
      #+begin_src c :comments link :tangle (eval c_test)
 {
diff --git a/org/qmckl_context.org b/org/qmckl_context.org
index 624681b..6a35789 100644
--- a/org/qmckl_context.org
+++ b/org/qmckl_context.org
@@ -62,6 +62,7 @@ int main() {
 #include <stdio.h>
 #include <string.h>
 #include <errno.h>
+#include <pthread.h>
 
 #include "qmckl.h"
 #include "qmckl_context_private_type.h"
@@ -69,6 +70,8 @@ int main() {
 
   #+end_src
 
+
+
 * Context handling
 
   The context variable is a handle for the state of the library,
@@ -146,7 +149,7 @@ typedef struct qmckl_context_struct {
    the pointer associated with a context is non-null, we can still
    verify that it points to the expected data structure.
 
-   #+begin_src c :comments org :tangle (eval h_private_type) :noweb yes
+   #+begin_src c :comments org :tangle (eval h_private_type) :noweb yes :exports none
 #define VALID_TAG   0xBEEFFACE
 #define INVALID_TAG 0xDEADBEEF
    #+end_src
@@ -156,10 +159,11 @@ typedef struct qmckl_context_struct {
    if the context is valid, ~QMCKL_NULL_CONTEXT~ otherwise.
 
    #+begin_src c :comments org :tangle (eval h_func) :noexport
-qmckl_context qmckl_context_check(const qmckl_context context) ;
+qmckl_context
+qmckl_context_check (const qmckl_context context) ;
    #+end_src
 
-   #+begin_src c :tangle (eval c)
+   #+begin_src c :tangle (eval c) :exports none
 qmckl_context qmckl_context_check(const qmckl_context context) {
 
   if (context == QMCKL_NULL_CONTEXT)
@@ -176,19 +180,27 @@ qmckl_context qmckl_context_check(const qmckl_context context) {
 }
    #+end_src
 
-   The context keeps a ``date'' that allows to check which data needs
+   The context keeps a /date/ that allows to check which data needs
    to be recomputed. The date is incremented when the context is touched.
 
    When a new element is added to the context, the functions
-   [[Creation][qmckl_context_create]], [[Destroy][qmckl_context_destroy]] and [[Copy][qmckl_context_copy]]
+   [[Creation][=qmckl_context_create=]] [[Destroy][=qmckl_context_destroy=]] and [[Copy][=qmckl_context_copy=]]
    should be updated in order to make deep copies.
 
-   #+begin_src c :comments org :tangle (eval h_func) :noexport
-qmckl_exit_code qmckl_context_touch(const qmckl_context context) ;
+   When the electron coordinates have changed, the context is touched
+   using the following function.
+
+   #+begin_src c :comments org :tangle (eval h_func)
+qmckl_exit_code
+qmckl_context_touch (const qmckl_context context);
    #+end_src
 
-   #+begin_src c :tangle (eval c)
-qmckl_exit_code qmckl_context_touch(const qmckl_context context) {
+   This has the effect to increment the date of the context.
+
+   #+begin_src c :tangle (eval c) :exports none
+qmckl_exit_code
+qmckl_context_touch(const qmckl_context context)
+{
 
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
     return qmckl_failwith( context,
@@ -213,12 +225,12 @@ qmckl_exit_code qmckl_context_touch(const qmckl_context context) {
    - A new context always has all its members initialized with a NULL value
 
    # Header
-   #+begin_src c :comments org :tangle (eval h_func) :exports none
+   #+begin_src c :comments org :tangle (eval h_func)
 qmckl_context qmckl_context_create();
    #+end_src
 
    # Source
-   #+begin_src c :tangle (eval c)
+   #+begin_src c :tangle (eval c) :exports none
 qmckl_context qmckl_context_create() {
 
   qmckl_context_struct* const ctx =
@@ -241,7 +253,9 @@ qmckl_context qmckl_context_create() {
     rc = pthread_mutexattr_init(&attr);
     assert (rc == 0);
 
+#ifdef PTHREAD_MUTEX_RECURSIVE
     (void) pthread_mutexattr_settype(&attr, PTHREAD_MUTEX_RECURSIVE);
+#endif
 
     rc = pthread_mutex_init ( &(ctx->mutex), &attr);
     assert (rc == 0);
@@ -252,30 +266,30 @@ qmckl_context qmckl_context_create() {
   /* Initialize data */
   {
     ctx->tag = VALID_TAG;
-    
+
     const qmckl_context context = (const qmckl_context) ctx;
     assert ( qmckl_context_check(context) != QMCKL_NULL_CONTEXT );
-    
+
     qmckl_exit_code rc;
-    
+
     ctx->numprec.precision = QMCKL_DEFAULT_PRECISION;
     ctx->numprec.range = QMCKL_DEFAULT_RANGE;
 
     rc = qmckl_init_point(context);
     assert (rc == QMCKL_SUCCESS);
-    
+
     rc = qmckl_init_electron(context);
     assert (rc == QMCKL_SUCCESS);
-    
+
     rc = qmckl_init_nucleus(context);
     assert (rc == QMCKL_SUCCESS);
-    
+
     rc = qmckl_init_ao_basis(context);
     assert (rc == QMCKL_SUCCESS);
-    
+
     rc = qmckl_init_mo_basis(context);
     assert (rc == QMCKL_SUCCESS);
-    
+
     rc = qmckl_init_determinant(context);
     assert (rc == QMCKL_SUCCESS);
   }
@@ -326,13 +340,13 @@ assert( qmckl_context_check(context) ==  context );
    ~lock_count~ attribute.
 
    # Header
-   #+begin_src c :comments org :tangle (eval h_func) :exports none
+   #+begin_src c :comments org :tangle (eval h_func)
 void qmckl_lock  (qmckl_context context);
 void qmckl_unlock(qmckl_context context);
    #+end_src
 
    # Source
-   #+begin_src c :tangle (eval c)
+   #+begin_src c :tangle (eval c) :exports none
 void qmckl_lock(qmckl_context context) {
   if (context == QMCKL_NULL_CONTEXT)
     return ;
@@ -345,9 +359,6 @@ void qmckl_lock(qmckl_context context) {
   }
   assert (rc == 0);
   ctx->lock_count += 1;
-/*
-  printf("  lock : %d\n", ctx->lock_count);
-,*/
 }
 
 void qmckl_unlock(const qmckl_context context) {
@@ -359,9 +370,6 @@ void qmckl_unlock(const qmckl_context context) {
   }
   assert (rc == 0);
   ctx->lock_count -= 1;
-/*
-  printf("unlock : %d\n", ctx->lock_count);
-,*/
 }
    #+end_src
 
@@ -372,12 +380,14 @@ void qmckl_unlock(const qmckl_context context) {
 
    # Header
    #+begin_src c :comments org :tangle (eval h_func) :exports none
+/*
 qmckl_context qmckl_context_copy(const qmckl_context context);
+*/
    #+end_src
 
 
    # Source
-   #+begin_src c :tangle (eval c)
+   #+begin_src c :tangle (eval c) :exports none
 qmckl_context qmckl_context_copy(const qmckl_context context) {
 
   const qmckl_context checked_context = qmckl_context_check(context);
@@ -418,13 +428,13 @@ qmckl_context qmckl_context_copy(const qmckl_context context) {
 
    # Fortran interface
    #+begin_src f90 :tangle (eval fh_func) :exports none
-  interface
-     integer (qmckl_context) function qmckl_context_copy(context) bind(C)
-       use, intrinsic :: iso_c_binding
-       import
-       integer (qmckl_context), intent(in), value :: context
-     end function qmckl_context_copy
-  end interface
+! interface
+!    integer (qmckl_context) function qmckl_context_copy(context) bind(C)
+!      use, intrinsic :: iso_c_binding
+!      import
+!      integer (qmckl_context), intent(in), value :: context
+!    end function qmckl_context_copy
+! end interface
    #+end_src
 
    # Test
@@ -443,13 +453,16 @@ munit_assert_int64(qmckl_context_check(new_context), ==, new_context);
    It frees the context, and returns the previous context.
 
    # Header
-   #+begin_src c :comments org :tangle (eval h_func) :exports none
-qmckl_exit_code qmckl_context_destroy(const qmckl_context context);
+   #+begin_src c :comments org :tangle (eval h_func)
+qmckl_exit_code
+qmckl_context_destroy (const qmckl_context context);
    #+end_src
 
     # Source
-    #+begin_src c :tangle (eval c)
-qmckl_exit_code qmckl_context_destroy(const qmckl_context context) {
+    #+begin_src c :tangle (eval c) :exports none
+qmckl_exit_code
+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;
diff --git a/org/qmckl_error.org b/org/qmckl_error.org
index cecd42d..97cc7cc 100644
--- a/org/qmckl_error.org
+++ b/org/qmckl_error.org
@@ -46,7 +46,7 @@ int main() {
 
    #+end_src
 
-*
+* -
 :PROPERTIES:
 :UNNUMBERED: t
 :END:
@@ -204,18 +204,25 @@ return '\n'.join(result)
    string is assumed to be large enough to contain the error message
    (typically 128 characters).
 
+* hidden                                                           :noexport:
+
+   #+NAME: MAX_STRING_LENGTH
+   : 128
+
 * Decoding errors
 
    To decode the error messages, ~qmckl_string_of_error~ converts an
    error code into a string.
 
-   #+NAME: MAX_STRING_LENGTH
-   : 128
+   #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+const char*
+qmckl_string_of_error (const qmckl_exit_code error);
+   #+end_src
 
-   #+begin_src c :comments org :tangle (eval h_func) :exports none :noweb yes
-const char* qmckl_string_of_error(const qmckl_exit_code error);
-void qmckl_string_of_error_f(const qmckl_exit_code error,
-                             char result[<<MAX_STRING_LENGTH()>>]);
+   #+begin_src c :comments org :tangle (eval h_private_func) :exports none :noweb yes
+void
+qmckl_string_of_error_f(const qmckl_exit_code error,
+                        char result[<<MAX_STRING_LENGTH()>>]);
    #+end_src
 
    The text strings are extracted from the previous table.
@@ -326,7 +333,7 @@ return '\n'.join(result)
    #+end_example
 
    # Source
-   #+begin_src c :comments org :tangle (eval c) :noweb yes
+   #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
 const char* qmckl_string_of_error(const qmckl_exit_code error) {
   switch (error) {
   <<cases()>>
@@ -340,7 +347,7 @@ void qmckl_string_of_error_f(const qmckl_exit_code error, char result[<<MAX_STRI
    #+end_src
 
   # Fortran interface
-  #+begin_src f90 :tangle (eval fh_func) :noexport :noweb yes
+  #+begin_src f90 :tangle (eval fh_func) :exports none :noweb yes
   interface
      subroutine qmckl_string_of_error (error, string) bind(C, name='qmckl_string_of_error_f')
        use, intrinsic :: iso_c_binding
@@ -353,7 +360,7 @@ void qmckl_string_of_error_f(const qmckl_exit_code error, char result[<<MAX_STRI
 
 * Data structure in context
 
-  The strings are declared with a maximum fixed size to avoid
+  The strings are declared internally with a maximum fixed size to avoid
   dynamic memory allocation.
 
    #+begin_src c :comments org :tangle (eval h_private_type)
@@ -377,7 +384,7 @@ typedef struct qmckl_error_struct {
    explaining the error. The exit code can't be ~QMCKL_SUCCESS~.
 
    # Header
-   #+begin_src c :comments org :tangle (eval h_func) :exports none
+   #+begin_src c :comments org :tangle (eval h_func)
 qmckl_exit_code
 qmckl_set_error(qmckl_context context,
                 const qmckl_exit_code exit_code,
@@ -386,7 +393,7 @@ qmckl_set_error(qmckl_context context,
    #+end_src
 
   # Source
-  #+begin_src c :tangle (eval c)
+  #+begin_src c :tangle (eval c) :exports none
 qmckl_exit_code
 qmckl_set_error(qmckl_context context,
                 const qmckl_exit_code exit_code,
@@ -428,7 +435,7 @@ qmckl_set_error(qmckl_context context,
   function name and message is mandatory.
 
   # Header
-   #+begin_src c :comments org :tangle (eval h_func) :exports none
+   #+begin_src c :comments org :tangle (eval h_func)
 qmckl_exit_code
 qmckl_get_error(qmckl_context context,
                 qmckl_exit_code *exit_code,
@@ -437,7 +444,7 @@ qmckl_get_error(qmckl_context context,
    #+end_src
 
   # Source
-  #+begin_src c :tangle (eval c)
+  #+begin_src c :tangle (eval c) :exports none
 qmckl_exit_code
 qmckl_get_error(qmckl_context context,
                 qmckl_exit_code *exit_code,
@@ -488,19 +495,21 @@ qmckl_get_error(qmckl_context context,
    the desired return code.
    Upon failure, a ~QMCKL_NULL_CONTEXT~ is returned.
 
-   #+begin_src c :comments org :tangle (eval h_func) :exports none
-qmckl_exit_code qmckl_failwith(qmckl_context context,
-                               const qmckl_exit_code exit_code,
-                               const char* function,
-                               const char* message) ;
+   #+begin_src c :comments org :tangle (eval h_func)
+qmckl_exit_code
+qmckl_failwith(qmckl_context context,
+               const qmckl_exit_code exit_code,
+               const char* function,
+               const char* message) ;
    #+end_src
    
-   #+begin_src c :comments org :tangle (eval c)
-qmckl_exit_code qmckl_failwith(qmckl_context context,
-                               const qmckl_exit_code exit_code,
-                               const char* function,
-                               const char* message) {
-
+   #+begin_src c :comments org :tangle (eval c) :exports none
+qmckl_exit_code
+qmckl_failwith(qmckl_context context,
+               const qmckl_exit_code exit_code,
+               const char* function,
+               const char* message)
+{
   assert (exit_code > 0);
   assert (exit_code < QMCKL_INVALID_EXIT_CODE);
   assert (function != NULL);
@@ -544,8 +553,6 @@ if (x < 0) {
 #endif
   #+end_src
 
-
-** Test
    #+begin_src c :comments link :tangle (eval c_test)
   /* Initialize the variables */
   char function_name[QMCKL_MAX_FUN_LEN]="";
diff --git a/org/qmckl_nucleus.org b/org/qmckl_nucleus.org
index 0191afe..319d5d1 100644
--- a/org/qmckl_nucleus.org
+++ b/org/qmckl_nucleus.org
@@ -187,6 +187,19 @@ qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num) {
 }
    #+end_src
 
+    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
+interface
+  integer(c_int32_t) function qmckl_get_nucleus_num(context, num) &
+    bind(C)
+    use, intrinsic :: iso_c_binding
+    import
+    implicit none
+    integer (c_int64_t) , intent(in)  , value :: context
+    integer (c_int64_t) , intent(out)         :: num
+  end function qmckl_get_nucleus_num
+end interface
+    #+end_src
+
 
    #+begin_src c :comments org :tangle (eval h_func) :exports none
 qmckl_exit_code
@@ -242,6 +255,20 @@ qmckl_get_nucleus_charge (const qmckl_context context,
 
    #+end_src
 
+    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
+interface
+  integer(c_int32_t) function qmckl_get_nucleus_charge(context, charge, size_max) &
+    bind(C)
+    use, intrinsic :: iso_c_binding
+    import
+    implicit none
+    integer (c_int64_t) , intent(in)  , value :: context
+    real    (c_double)  , intent(out)         :: charge(*)
+    integer (c_int64_t) , intent(in)  , value :: size_max
+  end function qmckl_get_nucleus_charge
+end interface
+    #+end_src
+
 
    #+begin_src c :comments org :tangle (eval h_func) :exports none
 qmckl_exit_code
@@ -277,6 +304,18 @@ qmckl_get_nucleus_rescale_factor (const qmckl_context context,
 }
    #+end_src
 
+    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
+interface
+  integer(c_int32_t) function qmckl_get_nucleus_rescale_factor(context, kappa) &
+    bind(C)
+    use, intrinsic :: iso_c_binding
+    import
+    implicit none
+    integer (c_int64_t) , intent(in)  , value :: context
+    real    (c_double)  , intent(out)         :: kappa
+  end function qmckl_get_nucleus_rescale_factor
+end interface
+    #+end_src
 
    #+begin_src c :comments org :tangle (eval h_func) :exports none
 qmckl_exit_code
@@ -342,6 +381,21 @@ qmckl_get_nucleus_coord (const qmckl_context context,
 }
    #+end_src
 
+    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
+interface
+  integer(c_int32_t) function qmckl_get_nucleus_coord(context, transp, coord, size_max) &
+    bind(C)
+    use, intrinsic :: iso_c_binding
+    import
+    implicit none
+    integer (c_int64_t) , intent(in)  , value :: context
+    character(c_char)   , intent(in)  , value :: transp
+    real    (c_double)  , intent(out)         :: coord(*)
+    integer (c_int64_t) , intent(in)  , value :: size_max
+  end function qmckl_get_nucleus_coord
+end interface
+    #+end_src
+
    When all the data relative to nuclei have been set, the following
    function returns ~true~.
 
@@ -424,7 +478,7 @@ interface
     implicit none
     integer (c_int64_t) , intent(in)  , value :: context
     integer (c_int64_t) , intent(in)  , value :: num
-  end function
+  end function qmckl_set_nucleus_num
 end interface
     #+end_src
 
@@ -570,7 +624,7 @@ interface
     character(c_char)   , intent(in)  , value :: transp
     real    (c_double)  , intent(in)          :: coord(*)
     integer (c_int64_t) , intent(in)  , value :: size_max
-  end function
+  end function qmckl_set_nucleus_coord
 end interface
     #+end_src
 
@@ -604,14 +658,14 @@ qmckl_set_nucleus_rescale_factor(qmckl_context context,
 
     #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
 interface
-  integer(c_int32_t) function qmckl_set_rescale_factor(context, kappa) &
+  integer(c_int32_t) function qmckl_set_nucleus_rescale_factor(context, kappa) &
     bind(C)
     use, intrinsic :: iso_c_binding
     import
     implicit none
     integer (c_int64_t) , intent(in)  , value :: context
     real    (c_double)  , intent(in)  , value :: kappa
-  end function
+  end function qmckl_set_nucleus_rescale_factor
 end interface
     #+end_src