From 1af0bf053f057018b553534a78239c2037feb1c0 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 18 Mar 2021 18:02:06 +0100
Subject: [PATCH] Memory checking

---
 src/qmckl.org            |  11 +-
 src/qmckl_ao.org         | 491 ++++++++++++++++++++++-----------------
 src/qmckl_context.org    | 126 +++++++---
 src/qmckl_error.org      |   4 +
 src/test_qmckl.org       |  66 +++---
 tools/Building.org       |  85 +++----
 tools/create_makefile.sh |  13 +-
 7 files changed, 470 insertions(+), 326 deletions(-)

diff --git a/src/qmckl.org b/src/qmckl.org
index ff58ace..3a568c7 100644
--- a/src/qmckl.org
+++ b/src/qmckl.org
@@ -116,8 +116,15 @@ Both files are located in the =include/= directory.
     For more guidelines on using Fortran to generate a C interface, see
     [[http://fortranwiki.org/fortran/show/Generating+C+Interfaces][this link]].
 
+** Coding rules
+
+    The authors should follow the recommendations of the
+    [[https://wiki.sei.cmu.edu/confluence/display/c/SEI+CERT+C+Coding+Standard][SEI+CERT+C+Coding+Standard]].
+
+    - Store a new value in pointers immediately after the memory is
+      freed
+    - Free dynamically allocated memory when no longer needed
 
- #  Coding style
  #  # TODO: decide on a coding style
 
  #  To improve  readability, we maintain  a consistent coding  style in
@@ -129,6 +136,7 @@ Both files are located in the =include/= directory.
 
  #  Coding style can be automatically checked with [[https://clang.llvm.org/docs/ClangFormat.html][clang-format]].
 
+
 ** Design of the library
 
    The proposed API should allow the library to: deal with memory transfers
@@ -251,3 +259,4 @@ Both files are located in the =include/= directory.
 
 
 
+    
diff --git a/src/qmckl_ao.org b/src/qmckl_ao.org
index 5430b06..f997913 100644
--- a/src/qmckl_ao.org
+++ b/src/qmckl_ao.org
@@ -1,8 +1,36 @@
 #+TITLE: Atomic Orbitals
 #+SETUPFILE: ../docs/theme.setup
 
- The routines for the computation of the values, gradients and
- Laplacian of atomic basis functions are defined here.
+The atomic basis set is defined as a list of shells. Each shell $s$ is
+centered on a nucleus $A$, possesses a given angular momentum $l$ and a
+radial function $R_s$. The radial function is a linear combination of
+\emph{primitive} functions that can be of type Slater ($p=1$)  or
+Gaussian ($p=2$):
+
+\[
+  R_s(\mathbf{r}) = \mathcal{N}_s |\mathbf{r}-\mathbf{R}_A|^{n_s}
+  \sum_{k=1}^{N_{\text{prim}}} a_{ks}
+ \exp \left( - \gamma_{ks} | \mathbf{r}-\mathbf{R}_A | ^p \right). |
+\]
+
+In the case of Gaussian functions, $n_s$ is always zero.
+The normalization factor $\mathcal{N}_s$ ensures that all the functions
+of the shell are normalized to unity. As this normalization requires
+the ability to compute overlap integrals, it should be written in the
+file to ensure that the file is self-contained and does not require
+the client program to have the ability to compute such integrals.
+
+Atomic orbitals (AOs) are defined as
+
+\[
+\chi_i (\mathbf{r}) = P_{\eta(i)}(\mathbf{r})\, R_{\theta(i)} (\mathbf{r})
+\]
+    
+where $\theta(i)$ returns the shell on which the AO is expanded,
+and $\eta(i)$ denotes which angular function is chosen.
+
+In this section we describe the kernels used to compute the values,
+gradients and Laplacian of the atomic basis functions.
 
 * Headers                                                          :noexport:
 
@@ -20,59 +48,33 @@ MunitResult test_<<filename()>>() {
   context = qmckl_context_create();
   #+end_src
 
-* Polynomials
+* Polynomial part
 
-  \[
-  P_l(\mathbf{r},\mathbf{R}_i)  =   (x-X_i)^a (y-Y_i)^b (z-Z_i)^c 
-  \]
-
-  \begin{eqnarray*} 
-  \frac{\partial }{\partial x} P_l\left(\mathbf{r},\mathbf{R}_i \right) &
-                 = & a (x-X_i)^{a-1} (y-Y_i)^b (z-Z_i)^c \\
-  \frac{\partial }{\partial y} P_l\left(\mathbf{r},\mathbf{R}_i \right) &
-                 = & b (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c \\
-  \frac{\partial }{\partial z} P_l\left(\mathbf{r},\mathbf{R}_i \right) &
-                 = & c (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1} \\
-  \end{eqnarray*} 
-
-    
-  \begin{eqnarray*} 
-  \left( \frac{\partial }{\partial x^2} + 
-             \frac{\partial }{\partial y^2} + 
-             \frac{\partial }{\partial z^2} \right) P_l
-             \left(\mathbf{r},\mathbf{R}_i \right) &  = &  
-           a(a-1) (x-X_i)^{a-2} (y-Y_i)^b (z-Z_i)^c + \\
-        && b(b-1) (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c + \\
-        && c(c-1) (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1} 
-  \end{eqnarray*}
-
-** <<<~qmckl_ao_power~>>>
-
-   Computes all the powers of the ~n~ input data up to the given
-   maximum value given in input for each of the $n$ points:
+** Powers of $x-X_i$
    
-   \[ P_{ij} = X_j^i \]
+   The ~qmckl_ao_power~ function computes all the powers of the ~n~
+   input data up to the given maximum value given in input for each of
+   the $n$ points:
+   
+   \[ P_{ik} = X_i^k \]
 
-*** Arguments
-
-    | ~context~  | input  | Global state                                      |
-    | ~n~        | input  | Number of values                                  |
-    | ~X(n)~     | input  | Array containing the input values                 |
-    | ~LMAX(n)~  | input  | Array containing the maximum power for each value |
-    | ~P(LDP,n)~ | output | Array containing all the powers of ~X~            |
-    | ~LDP~      | input  | Leading dimension of array ~P~                    |
+   | ~context~  | input  | Global state                                      |
+   | ~n~        | input  | Number of values                                  |
+   | ~X(n)~     | input  | Array containing the input values                 |
+   | ~LMAX(n)~  | input  | Array containing the maximum power for each value |
+   | ~P(LDP,n)~ | output | Array containing all the powers of ~X~            |
+   | ~LDP~      | input  | Leading dimension of array ~P~                    |
     
-*** Requirements
+   Requirements:
 
-    - ~context~ is not 0
-    - ~n~ > 0
-    - ~X~ is allocated with at least $n \times 8$ bytes
-    - ~LMAX~ is allocated with at least $n \times 4$ bytes
-    - ~P~ is allocated with at least $n \times \max_i \text{LMAX}_i \times 8$ bytes
-    - ~LDP~ >= $\max_i$ ~LMAX[i]~
+   - ~context~ is not ~QMCKL_NULL_CONTEXT~
+   - ~n~ > 0
+   - ~X~ is allocated with at least $n \times 8$ bytes
+   - ~LMAX~ is allocated with at least $n \times 4$ bytes
+   - ~P~ is allocated with at least $n \times \max_i \text{LMAX}_i \times 8$ bytes
+   - ~LDP~ >= $\max_i$ ~LMAX[i]~
 
-*** Header
-    #+begin_src c :tangle (eval h) 
+   #+begin_src c :tangle (eval h) 
 qmckl_exit_code
 qmckl_ao_power(const qmckl_context context,
                const int64_t n, 
@@ -80,11 +82,11 @@ qmckl_ao_power(const qmckl_context context,
                const int32_t *LMAX,
                const double *P,
                const int64_t LDP);
-    #+end_src 
-    
-*** Source
-    #+begin_src f90 :tangle (eval f) 
+   #+end_src 
+
+   #+begin_src f90 :tangle (eval f) 
 integer function qmckl_ao_power_f(context, n, X, LMAX, P, ldp) result(info)
+  use qmckl
   implicit none
   integer*8 , intent(in)  :: context
   integer*8 , intent(in)  :: n
@@ -93,32 +95,42 @@ integer function qmckl_ao_power_f(context, n, X, LMAX, P, ldp) result(info)
   real*8    , intent(out) :: P(ldp,n)
   integer*8 , intent(in)  :: ldp
 
-  integer*8  :: i,j
+  integer*8  :: i,k
 
-  info = 0
+  info = QMCKL_SUCCESS
 
-  if (context == 0_8) then
-     info = -1
+  if (context == QMCKL_NULL_CONTEXT) then
+     info = QMCKL_INVALID_CONTEXT
      return
   endif
-  
-  if (LDP < MAXVAL(LMAX)) then
-     info = -2
+
+  if (n <= ldp) then
+     info = QMCKL_INVALID_ARG_2
      return
   endif
-  
-  do j=1,n
-    P(1,j) = X(j)
-    do i=2,LMAX(j)
-       P(i,j) = P(i-1,j) * X(j) 
-    end do
+
+  k = MAXVAL(LMAX)
+  if (LDP < k) then
+     info = QMCKL_INVALID_ARG_6
+     return
+  endif
+
+  if (k <= 0) then
+     info = QMCKL_INVALID_ARG_4
+     return
+  endif
+
+  do i=1,n
+     P(1,i) = X(i)
+     do k=2,LMAX(i)
+        P(k,i) = P(k-1,i) * X(i) 
+     end do
   end do
 
 end function qmckl_ao_power_f
-    #+end_src 
+   #+end_src 
 
-*** C interface                                                    :noexport:
-    #+begin_src f90 :tangle (eval f) 
+   #+begin_src f90 :tangle (eval f) :exports none
 integer(c_int32_t) function qmckl_ao_power(context, n, X, LMAX, P, ldp) &
      bind(C) result(info)
   use, intrinsic :: iso_c_binding
@@ -129,13 +141,13 @@ integer(c_int32_t) function qmckl_ao_power(context, n, X, LMAX, P, ldp) &
   integer (c_int32_t) , intent(in)         :: LMAX(n)
   real    (c_double)  , intent(out)        :: P(ldp,n)
   integer (c_int64_t) , intent(in) , value :: ldp
-  
+
   integer, external :: qmckl_ao_power_f
   info = qmckl_ao_power_f(context, n, X, LMAX, P, ldp)
 end function qmckl_ao_power
-    #+end_src
+   #+end_src
 
-    #+begin_src f90 :tangle (eval fh) 
+   #+begin_src f90 :tangle (eval fh) :exports none
   interface
      integer(c_int32_t) function qmckl_ao_power(context, n, X, LMAX, P, ldp) bind(C)
        use, intrinsic :: iso_c_binding
@@ -147,16 +159,16 @@ end function qmckl_ao_power
        real    (c_double)  , intent(out)        :: P(ldp,n)
      end function qmckl_ao_power
   end interface
-    #+end_src
-  
-*** Test                                                           :noexport:
-    #+begin_src f90 :tangle (eval f_test) 
+   #+end_src
+
+   # Test
+   #+begin_src f90 :tangle (eval f_test) 
 integer(c_int32_t) function test_qmckl_ao_power(context) bind(C)
   use qmckl
   implicit none
 
   integer(c_int64_t), intent(in), value :: context
-  
+
   integer*8                     :: n, LDP 
   integer, allocatable          :: LMAX(:) 
   double precision, allocatable :: X(:), P(:,:)
@@ -167,19 +179,19 @@ integer(c_int32_t) function test_qmckl_ao_power(context) bind(C)
 
   n = 100;
   LDP = 10;
-  
+
   allocate(X(n), P(LDP,n), LMAX(n))
-  
+
   do j=1,n
      X(j) = -5.d0 + 0.1d0 * dble(j)
      LMAX(j) = 1 + int(mod(j, 5),4)
   end do
-  
+
   test_qmckl_ao_power = qmckl_ao_power(context, n, X, LMAX, P, LDP) 
   if (test_qmckl_ao_power /= 0) return
-  
+
   test_qmckl_ao_power = -1
-  
+
   do j=1,n
      do i=1,LMAX(j)
         if ( X(j)**i == 0.d0 ) then
@@ -193,58 +205,79 @@ integer(c_int32_t) function test_qmckl_ao_power(context) bind(C)
   test_qmckl_ao_power = 0
   deallocate(X,P,LMAX)
 end function test_qmckl_ao_power
-    #+end_src
+   #+end_src
 
-    #+begin_src c :tangle (eval c_test)
+   #+begin_src c :tangle (eval c_test) :exports none
 int  test_qmckl_ao_power(qmckl_context context);
 munit_assert_int(0, ==, test_qmckl_ao_power(context));
-    #+end_src
+   #+end_src
+
+** Value, Gradient and Laplacian of a polynomial
+
+   A polynomial is centered on a nucleus $\mathbf{R}_i$
   
-** <<<~qmckl_ao_polynomial_vgl~>>>
-   
-   Computes the values, gradients and Laplacians at a given point of
-   all polynomials with an angular momentum up to ~lmax~.
-   
-*** Arguments
+   \[
+   P_l(\mathbf{r},\mathbf{R}_i)  =   (x-X_i)^a (y-Y_i)^b (z-Z_i)^c 
+   \]
 
-    | ~context~    | input  | Global state                                         |
-    | ~X(3)~       | input  | Array containing the coordinates of the points       |
-    | ~R(3)~       | input  | Array containing the x,y,z coordinates of the center |
-    | ~lmax~       | input  | Maximum angular momentum                             |
-    | ~n~          | output | Number of computed polynomials                       |
-    | ~L(ldl,n)~   | output | Contains a,b,c for all ~n~ results                   |
-    | ~ldl~        | input  | Leading dimension of ~L~                             |
-    | ~VGL(ldv,n)~ | output | Value, gradients and Laplacian of the polynomials    |
-    | ~ldv~        | input  | Leading dimension of array ~VGL~                     |
+   The gradients with respect to electron coordinates are
+  
+   \begin{eqnarray*} 
+   \frac{\partial }{\partial x} P_l\left(\mathbf{r},\mathbf{R}_i \right) &
+                  = & a (x-X_i)^{a-1} (y-Y_i)^b (z-Z_i)^c \\
+   \frac{\partial }{\partial y} P_l\left(\mathbf{r},\mathbf{R}_i \right) &
+                  = & b (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c \\
+   \frac{\partial }{\partial z} P_l\left(\mathbf{r},\mathbf{R}_i \right) &
+                  = & c (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1} \\
+   \end{eqnarray*} 
+
+   and the Laplacian is
     
-*** Requirements
+   \begin{eqnarray*} 
+   \left( \frac{\partial }{\partial x^2} + 
+              \frac{\partial }{\partial y^2} + 
+              \frac{\partial }{\partial z^2} \right) P_l
+              \left(\mathbf{r},\mathbf{R}_i \right) &  = &  
+            a(a-1) (x-X_i)^{a-2} (y-Y_i)^b (z-Z_i)^c + \\
+         && b(b-1) (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c + \\
+         && c(c-1) (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1}.
+   \end{eqnarray*}
 
-    - ~context~ is not 0
-    - ~n~ > 0
-    - ~lmax~ >= 0
-    - ~ldl~ >= 3
-    - ~ldv~ >= 5
-    - ~X~ is allocated with at least $3 \times 8$ bytes
-    - ~R~ is allocated with at least $3 \times 8$ bytes
-    - ~n~ >= ~(lmax+1)(lmax+2)(lmax+3)/6~
-    - ~L~ is allocated with at least $3 \times n \times 4$ bytes
-    - ~VGL~ is allocated with at least $5 \times n \times 8$ bytes
-    - On output, ~n~ should be equal to ~(lmax+1)(lmax+2)(lmax+3)/6~
-    - On output, the powers are given in the following order (l=a+b+c):
-      - Increase values of ~l~
-      - Within a given value of ~l~, alphabetical order of the
-        string made by a*"x" + b*"y" + c*"z" (in Python notation).
-        For example, with a=0, b=2 and c=1 the string is "yyz"
-        
-*** Error codes
+   ~qmckl_ao_polynomial_vgl~ computes the values, gradients and
+   Laplacians at a given point in space, of all polynomials with an
+   angular momentum up to ~lmax~.
 
-    | -1 | Null context        |
-    | -2 | Inconsistent ~ldl~  |
-    | -3 | Inconsistent ~ldv~  |
-    | -4 | Inconsistent ~lmax~ |
+   | ~context~    | input  | Global state                                         |
+   | ~X(3)~       | input  | Array containing the coordinates of the points       |
+   | ~R(3)~       | input  | Array containing the x,y,z coordinates of the center |
+   | ~lmax~       | input  | Maximum angular momentum                             |
+   | ~n~          | output | Number of computed polynomials                       |
+   | ~L(ldl,n)~   | output | Contains a,b,c for all ~n~ results                   |
+   | ~ldl~        | input  | Leading dimension of ~L~                             |
+   | ~VGL(ldv,n)~ | output | Value, gradients and Laplacian of the polynomials    |
+   | ~ldv~        | input  | Leading dimension of array ~VGL~                     |
+    
+   Requirements:
 
-*** Header
-    #+begin_src c :tangle (eval h) 
+   - ~context~ is not ~QMCKL_NULL_CONTEXT~
+   - ~n~ > 0
+   - ~lmax~ >= 0
+   - ~ldl~ >= 3
+   - ~ldv~ >= 5
+   - ~X~ is allocated with at least $3 \times 8$ bytes
+   - ~R~ is allocated with at least $3 \times 8$ bytes
+   - ~n~ >= ~(lmax+1)(lmax+2)(lmax+3)/6~
+   - ~L~ is allocated with at least $3 \times n \times 4$ bytes
+   - ~VGL~ is allocated with at least $5 \times n \times 8$ bytes
+   - On output, ~n~ should be equal to ~(lmax+1)(lmax+2)(lmax+3)/6~
+   - On output, the powers are given in the following order (l=a+b+c):
+     - Increasing values of ~l~
+     - Within a given value of ~l~, alphabetical order of the
+       string made by a*"x" + b*"y" + c*"z" (in Python notation).
+       For example, with a=0, b=2 and c=1 the string is "yyz"
+
+   # Header
+   #+begin_src c :tangle (eval h) 
 qmckl_exit_code
 qmckl_ao_polynomial_vgl(const qmckl_context context,
                         const double *X,
@@ -255,11 +288,12 @@ qmckl_ao_polynomial_vgl(const qmckl_context context,
                         const int64_t ldl,
                         const double *VGL,
                         const int64_t ldv);
-    #+end_src
-    
-*** Source
-    #+begin_src f90 :tangle (eval f)
+   #+end_src
+
+   # Source
+   #+begin_src f90 :tangle (eval f)
 integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, ldv) result(info)
+  use qmckl
   implicit none
   integer*8 , intent(in)  :: context
   real*8    , intent(in)  :: X(3), R(3)
@@ -278,69 +312,74 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL,
   integer, external :: qmckl_ao_power_f
   double precision  :: xy, yz, xz
   double precision  :: da, db, dc, dd
-  
+
   info = 0
-  
-  if (context == 0_8) then
-     info = -1
+
+  if (context == QMCKL_NULL_CONTEXT) then
+     info = QMCKL_INVALID_CONTEXT
      return
   endif
-  
-  if (ldl < 3) then
-     info = -2
-     return
-  endif
-  
-  if (ldv < 5) then
-     info = -3
-     return
-  endif
-  
+
   if (lmax <= 0) then
-     info = -4
+     info = QMCKL_INVALID_ARG_4
      return
   endif
-  
-  
+
+  if (n <= 0) then
+     info = QMCKL_INVALID_ARG_5
+     return
+  endif
+
+  if (ldl < 3) then
+     info = QMCKL_INVALID_ARG_7
+     return
+  endif
+
+  if (ldv < 5) then
+     info = QMCKL_INVALID_ARG_9
+     return
+  endif
+
+
   do i=1,3
      Y(i) = X(i) - R(i)
   end do
 
   lmax_array(1:3) = lmax
   if (lmax == 0) then
-    VGL(1,1) = 1.d0
-    vgL(2:5,1) = 0.d0
-    l(1:3,1) = 0
-    n=1
+     VGL(1,1) = 1.d0
+     vgL(2:5,1) = 0.d0
+     l(1:3,1) = 0
+     n=1
   else if (lmax > 0) then
-    pows(-2:0,1:3) = 1.d0
-    do i=1,lmax
+     pows(-2:0,1:3) = 1.d0
+     do i=1,lmax
         pows(i,1) = pows(i-1,1) * Y(1) 
         pows(i,2) = pows(i-1,2) * Y(2) 
         pows(i,3) = pows(i-1,3) * Y(3) 
-    end do
+     end do
 
-    VGL(1:5,1:4) = 0.d0
-    l(1:3,1:4) = 0
+     VGL(1:5,1:4) = 0.d0
+     l  (1:3,1:4) = 0
 
-    VGL(1,1) = 1.d0
-    vgl(1:5,2:4) = 0.d0
+     VGL(1  ,1  ) = 1.d0
+     vgl(1:5,2:4) = 0.d0
 
-    l(1,2) = 1
-    vgl(1,2) = pows(1,1)
-    vgL(2,2) = 1.d0
+     l  (1,2) = 1
+     vgl(1,2) = pows(1,1)
+     vgL(2,2) = 1.d0
 
-    l(2,3) = 1
-    vgl(1,3) = pows(1,2)
-    vgL(3,3) = 1.d0
+     l  (2,3) = 1
+     vgl(1,3) = pows(1,2)
+     vgL(3,3) = 1.d0
 
-    l(3,4) = 1
-    vgl(1,4) = pows(1,3)
-    vgL(4,4) = 1.d0
+     l  (3,4) = 1
+     vgl(1,4) = pows(1,3)
+     vgL(4,4) = 1.d0
 
-    n=4
+     n=4
   endif
-        
+
   ! l>=2
   dd = 2.d0
   do d=2,lmax
@@ -355,21 +394,21 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL,
            l(1,n) = a
            l(2,n) = b
            l(3,n) = c
-           
+
            xy = pows(a,1) * pows(b,2)
            yz = pows(b,2) * pows(c,3)
            xz = pows(a,1) * pows(c,3)
-           
+
            vgl(1,n) = xy * pows(c,3)
-           
+
            xy = dc * xy
            xz = db * xz
            yz = da * yz
-           
+
            vgl(2,n) = pows(a-1,1) * yz
            vgl(3,n) = pows(b-1,2) * xz
            vgl(4,n) = pows(c-1,3) * xy
-           
+
            vgl(5,n) = &
                 (da-1.d0) * pows(a-2,1) * yz + &
                 (db-1.d0) * pows(b-2,2) * xz + &
@@ -382,13 +421,13 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL,
      dd = dd + 1.d0
   end do
 
-  info = 0
+  info = QMCKL_SUCCESS
 
 end function qmckl_ao_polynomial_vgl_f
-    #+end_src 
+   #+end_src 
 
-*** C interface                                                    :noexport:
-    #+begin_src f90 :tangle (eval f)
+
+   #+begin_src f90 :tangle (eval f) :exports none
 integer(c_int32_t) function qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv) &
      bind(C) result(info)
   use, intrinsic :: iso_c_binding
@@ -405,10 +444,9 @@ integer(c_int32_t) function qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, l
   integer, external :: qmckl_ao_polynomial_vgl_f
   info = qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, ldv) 
 end function qmckl_ao_polynomial_vgl
-    #+end_src
+   #+end_src
 
-*** Fortran interface                                              :noexport:
-    #+begin_src f90 :tangle (eval fh)
+   #+begin_src f90 :tangle (eval fh) :exports none
   interface
      integer(c_int32_t) function qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv) &
           bind(C) 
@@ -423,15 +461,15 @@ end function qmckl_ao_polynomial_vgl
        real    (c_double)  , intent(out)        :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6)
      end function qmckl_ao_polynomial_vgl
   end interface
-    #+end_src
-*** Test                                                           :noexport:
-    #+begin_src f90 :tangle (eval f_test)
+   #+end_src
+
+   #+begin_src f90 :tangle (eval f_test)
 integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
   use qmckl
   implicit none
 
   integer(c_int64_t), intent(in), value :: context
-  
+
   integer                       :: lmax, d, i
   integer, allocatable          :: L(:,:)
   integer*8                     :: n, ldl, ldv, j
@@ -450,7 +488,7 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
   n = 0;
   ldl = 3;
   ldv = 100;
-  
+
   d = (lmax+1)*(lmax+2)*(lmax+3)/6
 
   allocate (L(ldl,d), VGL(ldv,d))
@@ -499,7 +537,7 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
              L(3,j) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-1) &
              )) > epsilon ) return
      end if
-     
+
      test_qmckl_ao_polynomial_vgl = -16
      w = 0.d0
      if (L(1,j) > 1) then
@@ -515,22 +553,21 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
   end do
 
   test_qmckl_ao_polynomial_vgl = 0
-     
+
   deallocate(L,VGL)
 end function test_qmckl_ao_polynomial_vgl
-    #+end_src
+   #+end_src
 
-    #+begin_src c :tangle (eval c_test)
+   #+begin_src c :tangle (eval c_test)
 int  test_qmckl_ao_polynomial_vgl(qmckl_context context);
 munit_assert_int(0, ==, test_qmckl_ao_polynomial_vgl(context));
-    #+end_src
+   #+end_src
    
 * Gaussian basis functions
-
-** <<<~qmckl_ao_gaussian_vgl~>>>
-   
-   Computes the values, gradients and Laplacians at a given point of
-   ~n~ Gaussian functions centered at the same point:
+  
+   ~qmckl_ao_gaussian_vgl~ computes the values, gradients and
+   Laplacians at a given point of ~n~ Gaussian functions centered at
+   the same point:
    
    \[ v_i = \exp(-a_i |X-R|^2) \]
    \[ \nabla_x v_i = -2 a_i (X_x -  R_x) v_i \]
@@ -538,8 +575,6 @@ munit_assert_int(0, ==, test_qmckl_ao_polynomial_vgl(context));
    \[ \nabla_z v_i = -2 a_i (X_z -  R_z) v_i \]
    \[ \Delta v_i = a_i (4 |X-R|^2 a_i - 6) v_i \]
 
-*** Arguments
-
     | ~context~    | input  | Global state                                         |
     | ~X(3)~       | input  | Array containing the coordinates of the points       |
     | ~R(3)~       | input  | Array containing the x,y,z coordinates of the center |
@@ -548,7 +583,7 @@ munit_assert_int(0, ==, test_qmckl_ao_polynomial_vgl(context));
     | ~VGL(ldv,5)~ | output | Value, gradients and Laplacian of the Gaussians      |
     | ~ldv~        | input  | Leading dimension of array ~VGL~                     |
     
-*** Requirements
+    Requirements :
 
     - ~context~ is not 0
     - ~n~ > 0
@@ -559,7 +594,6 @@ munit_assert_int(0, ==, test_qmckl_ao_polynomial_vgl(context));
     - ~A~ is allocated with at least $n \times 8$ bytes
     - ~VGL~ is allocated with at least $n \times 5 \times 8$ bytes
 
-*** Header
     #+begin_src c :tangle (eval h)
 qmckl_exit_code
 qmckl_ao_gaussian_vgl(const qmckl_context context,
@@ -571,9 +605,9 @@ qmckl_ao_gaussian_vgl(const qmckl_context context,
                       const int64_t ldv);
     #+end_src
     
-*** Source
     #+begin_src f90 :tangle (eval f)
 integer function qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv) result(info)
+  use qmckl
   implicit none
   integer*8 , intent(in)  :: context
   real*8    , intent(in)  :: X(3), R(3)
@@ -585,20 +619,20 @@ integer function qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv) result(i
   integer*8         :: i,j
   real*8            :: Y(3), r2, t, u, v
   
-  info = 0
+  info = QMCKL_SUCCESS
   
-  if (context == 0_8) then
-     info = -1
+  if (context == QMCKL_NULL_CONTEXT) then
+     info = QMCKL_INVALID_CONTEXT
      return
   endif
   
   if (n <= 0) then
-     info = -2
+     info = QMCKL_INVALID_ARG_4
      return
   endif
   
   if (ldv < n) then
-     info = -3
+     info = QMCKL_INVALID_ARG_7
      return
   endif
   
@@ -634,8 +668,7 @@ integer function qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv) result(i
 end function qmckl_ao_gaussian_vgl_f
     #+end_src 
 
-*** C interface                                                    :noexport:
-    #+begin_src f90 :tangle (eval f)
+    #+begin_src f90 :tangle (eval f) :exports none
 integer(c_int32_t) function qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv) &
      bind(C) result(info)
   use, intrinsic :: iso_c_binding
@@ -652,7 +685,7 @@ integer(c_int32_t) function qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv)
 end function qmckl_ao_gaussian_vgl
     #+end_src
 
-    #+begin_src f90 :tangle (eval fh)
+    #+begin_src f90 :tangle (eval fh) :exports none
   interface
      integer(c_int32_t) function qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv) &
           bind(C) 
@@ -665,8 +698,9 @@ end function qmckl_ao_gaussian_vgl
      end function qmckl_ao_gaussian_vgl
   end interface
     #+end_src
-*** Test                                                           :noexport:
-    #+begin_src f90 :tangle (eval f_test)
+
+  # Test
+  #+begin_src f90 :tangle (eval f_test)
 integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C)
   use qmckl
   implicit none
@@ -731,12 +765,12 @@ integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C)
      
   deallocate(VGL)
 end function test_qmckl_ao_gaussian_vgl
-    #+end_src
+  #+end_src
 
-    #+begin_src c :tangle (eval c_test) :exports none
+  #+begin_src c :tangle (eval c_test) :exports none
 int test_qmckl_ao_gaussian_vgl(qmckl_context context);
 munit_assert_int(0, ==, test_qmckl_ao_gaussian_vgl(context));
-    #+end_src
+  #+end_src
    
 * TODO Slater basis functions
   
@@ -749,8 +783,35 @@ munit_assert_int(0, ==, test_qmckl_ao_gaussian_vgl(context));
   return MUNIT_OK;
 }
   #+end_src
-  
+ 
+**✸ Compute file names                                                                       
+    #+begin_src emacs-lisp                                                                   
+; The following is required to compute the file names                                        
+                                                                                             
+(setq pwd (file-name-directory buffer-file-name))                                            
+(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))                       
+(setq f  (concat pwd name "_f.f90"))                                                         
+(setq fh (concat pwd name "_fh.f90"))                                                        
+(setq c  (concat pwd name ".c"))                                                             
+(setq h  (concat name ".h"))                                                                 
+(setq h_private  (concat name "_private.h"))                                                 
+(setq c_test  (concat pwd "test_" name ".c"))                                                
+(setq f_test  (concat pwd "test_" name "_f.f90"))                                            
+                                                                                             
+; Minted                                                                                     
+(require 'ox-latex)                                                                          
+(setq org-latex-listings 'minted)                                                            
+(add-to-list 'org-latex-packages-alist '("" "listings"))                                     
+(add-to-list 'org-latex-packages-alist '("" "color"))                                        
+                                                                                             
+    #+end_src                                                                                
+
+    #+RESULTS:
+    |   | color    |
+    |   | listings |
+         
 
 # -*- mode: org -*-
 # vim: syntax=c
   
+
diff --git a/src/qmckl_context.org b/src/qmckl_context.org
index 8a68835..3630a85 100644
--- a/src/qmckl_context.org
+++ b/src/qmckl_context.org
@@ -55,6 +55,7 @@ MunitResult test_<<filename()>>() {
 #include "qmckl_context_private.h"
 #include "qmckl_memory.h"
 
+#include <stdio.h>
   #+end_src
 
 * Context handling
@@ -103,6 +104,7 @@ typedef struct qmckl_context_struct {
   qmckl_memory_struct      * alloc;
 
   /* Thread lock */
+  int                        lock_count;
   pthread_mutex_t            mutex;
 
   /* Validity checking */
@@ -150,8 +152,10 @@ qmckl_context qmckl_context_check(const qmckl_context context) {
 
   const qmckl_context_struct* ctx = (qmckl_context_struct*) context;
 
-  if (ctx->tag != VALID_TAG)
-    return QMCKL_NULL_CONTEXT;
+  /* Try to access memory */
+  if (ctx->tag != VALID_TAG) {
+      return QMCKL_NULL_CONTEXT;
+  }
 
   return context;
 }
@@ -183,18 +187,7 @@ qmckl_context qmckl_context_create() {
   memset(ctx, 0, sizeof(qmckl_context_struct));
 
   /* Initialize lock */
-  pthread_mutexattr_t attr;
-  int rc;
-
-  rc = pthread_mutexattr_init(&attr); 
-  assert (rc == 0);
-
-  (void) pthread_mutexattr_settype(&attr, PTHREAD_MUTEX_ERRORCHECK);
-
-  rc = pthread_mutex_init ( &(ctx->mutex), &attr);
-  assert (rc == 0);
-
-  (void)pthread_mutexattr_destroy(&attr);
+  init_lock(&(ctx->mutex));
 
   /* Initialize data */
   ctx->tag = VALID_TAG;
@@ -218,7 +211,6 @@ qmckl_context qmckl_context_create() {
    # Test
    #+begin_src c :comments link :tangle (eval c_test) :exports none
 munit_assert_int64( qmckl_context_check(QMCKL_NULL_CONTEXT), ==,  QMCKL_NULL_CONTEXT);
-munit_assert_int64( qmckl_context_check(0x12345), ==,  QMCKL_NULL_CONTEXT);
 
 qmckl_context context = qmckl_context_create();
 munit_assert_int64( context, !=, QMCKL_NULL_CONTEXT );
@@ -228,29 +220,60 @@ munit_assert_int64( qmckl_context_check(context), ==,  context );
 ** Locking
 
    For thread safety, the context may be locked/unlocked. The lock is
-   initialized with the ~PTHREAD_MUTEX_ERRORCHECK~, so it is a bit
-   slower than a usual mutex but safer.
+   initialized with the ~PTHREAD_MUTEX_RECURSIVE~ attribute, and the
+   number of times the thread has locked it is saved in the
+   ~lock_count~ attribute.
 
    # Header
    #+begin_src c :comments org :tangle (eval h) :exports none
 void qmckl_lock  (qmckl_context context);
 void qmckl_unlock(qmckl_context context);
+
+void init_lock(pthread_mutex_t* mutex);
    #+end_src
 
    # Source
    #+begin_src c :tangle (eval c)
+void init_lock(pthread_mutex_t* mutex) {
+  pthread_mutexattr_t attr;
+  int rc;
+  
+  rc = pthread_mutexattr_init(&attr); 
+  assert (rc == 0);
+  
+  (void) pthread_mutexattr_settype(&attr, PTHREAD_MUTEX_RECURSIVE);
+  
+  rc = pthread_mutex_init ( mutex, &attr);
+  assert (rc == 0);
+  
+  (void)pthread_mutexattr_destroy(&attr);
+}
+
 void qmckl_lock(qmckl_context context) {
   if (context == QMCKL_NULL_CONTEXT)
     return ;
   qmckl_context_struct *ctx = (qmckl_context_struct*) context;
+  errno = 0;
   int rc = pthread_mutex_lock( &(ctx->mutex) );
+  if (rc != 0) {
+    fprintf(stderr, "qmckl_lock:%s\n", strerror(rc) );
+    fflush(stderr);
+  }
   assert (rc == 0);
+  ctx->lock_count++;
+  printf("  lock : %d\n", ctx->lock_count);
 }
 
 void qmckl_unlock(qmckl_context context) {
   qmckl_context_struct *ctx = (qmckl_context_struct*) context;
   int rc = pthread_mutex_unlock( &(ctx->mutex) );
+  if (rc != 0) {
+    fprintf(stderr, "qmckl_unlock:%s\n", strerror(rc) );
+    fflush(stderr);
+  }
   assert (rc == 0);
+  ctx->lock_count--;
+  printf("unlock : %d\n", ctx->lock_count);
 }
    #+end_src
 
@@ -292,7 +315,7 @@ qmckl_context qmckl_context_copy(const qmckl_context context) {
   /* Copy the old context on the new one */
   memcpy(new_ctx, old_ctx, sizeof(qmckl_context_struct));
 
-  qmckl_unlock(context);
+  qmckl_unlock( (qmckl_context) old_ctx );
   
   new_ctx->prev = old_ctx;
 
@@ -334,32 +357,42 @@ qmckl_context qmckl_context_destroy(qmckl_context context);
     #+begin_src c :tangle (eval c)
 qmckl_context qmckl_context_destroy(const qmckl_context context) {
 
-  qmckl_lock(context);
-  
   const qmckl_context checked_context = qmckl_context_check(context);
   if (checked_context == QMCKL_NULL_CONTEXT) return QMCKL_NULL_CONTEXT;
 
+  qmckl_lock(context);
+  
   qmckl_context_struct* ctx = (qmckl_context_struct*) context;
   assert (ctx != NULL);  /* Shouldn't be true because the context is valid */
 
-  const qmckl_context prev_context = (qmckl_context) ctx->prev;
-  memset(ctx, 0, sizeof(qmckl_context_struct));
-  ctx->tag = INVALID_TAG;
+  qmckl_unlock(context);
 
+  const int rc_destroy = pthread_mutex_destroy( &(ctx->mutex) );
+  if (rc_destroy != 0) {
+     fprintf(stderr, "qmckl_context_destroy: %s %d\n", strerror(rc_destroy), ctx->lock_count);
+     abort();
+  }
+
+  const qmckl_context prev_context = (qmckl_context) ctx->prev;
+  if (prev_context == QMCKL_NULL_CONTEXT) {
+    /* This is the first context, free all memory. */
+    struct qmckl_memory_struct* old = NULL;
+    while (ctx->alloc != NULL) {
+       old = ctx->alloc;
+       ctx->alloc = ctx->alloc->prev;
+       free(old->pointer);
+       old->pointer = NULL;
+       free(old);
+       old = NULL;
+    }
+  }
+  
+  ctx->tag = INVALID_TAG;
   const qmckl_exit_code rc = qmckl_free(context,ctx);
   assert (rc == QMCKL_SUCCESS);
 
-  if (prev_context == QMCKL_NULL_CONTEXT) {
-    /* This is the first context, free all memory. */
-    while (ctx->alloc != NULL) {
-       free(ctx->alloc->pointer);
-       ctx->alloc = ctx->alloc->prev;
-    }
-    int rc = pthread_mutex_destroy( &(ctx->mutex) );
-    assert (rc == 0);
-  }
-  
-  qmckl_unlock(context);
+  //memset(ctx, 0, sizeof(qmckl_context_struct));
+
 
   return prev_context;
 }
@@ -1094,8 +1127,8 @@ double qmckl_context_get_epsilon(const qmckl_context context);
    # Source
    #+begin_src c :tangle (eval c)
 double qmckl_context_get_epsilon(const qmckl_context context) {
-  const qmckl_context_struct* ctx = (qmckl_context_struct*) context;
-  return 1. /  (double) (1 << (1 - ctx->fp->precision));
+  const int precision = qmckl_context_get_precision(context);
+  return 1. /  (double) (1L << (precision-1));
 }
    #+end_src
 
@@ -1245,51 +1278,72 @@ qmckl_context_update_ao_basis(qmckl_context   context     , const char      type
   basis->shell_center  = (int64_t*) qmckl_malloc (context, shell_num * sizeof(int64_t));
   if (basis->shell_center == NULL) {
     qmckl_free(context, basis);
+    basis = NULL;
     return QMCKL_FAILURE;
   }
 
   basis->shell_ang_mom = (int32_t*) qmckl_malloc (context, shell_num * sizeof(int32_t));
   if (basis->shell_ang_mom == NULL) {
     qmckl_free(context, basis->shell_center);
+    basis->shell_center = NULL;
     qmckl_free(context, basis);
+    basis = NULL;
     return QMCKL_FAILURE;
   }
 
   basis->shell_prim_num= (int64_t*) qmckl_malloc (context, shell_num * sizeof(int64_t));
   if (basis->shell_prim_num == NULL) {
     qmckl_free(context, basis->shell_ang_mom);
+    basis->shell_ang_mom = NULL;
     qmckl_free(context, basis->shell_center);
+    basis->shell_center = NULL;
     qmckl_free(context, basis);
+    basis = NULL;
     return QMCKL_FAILURE;
   }
 
   basis->shell_factor  = (double *) qmckl_malloc (context, shell_num * sizeof(double ));
   if (basis->shell_factor == NULL) {
     qmckl_free(context, basis->shell_prim_num);
+    basis->shell_prim_num = NULL;
     qmckl_free(context, basis->shell_ang_mom);
+    basis->shell_ang_mom = NULL;
     qmckl_free(context, basis->shell_center);
+    basis->shell_center = NULL;
     qmckl_free(context, basis);
+    basis = NULL;
     return QMCKL_FAILURE;
   }
 
   basis->exponent      = (double *) qmckl_malloc (context, prim_num  * sizeof(double ));
   if (basis->exponent == NULL) {
     qmckl_free(context, basis->shell_factor);
+    basis->shell_factor = NULL;
     qmckl_free(context, basis->shell_prim_num);
+    basis->shell_prim_num = NULL;
     qmckl_free(context, basis->shell_ang_mom);
+    basis->shell_ang_mom = NULL;
     qmckl_free(context, basis->shell_center);
+    basis->shell_center = NULL;
     qmckl_free(context, basis);
+    basis = NULL;
     return QMCKL_FAILURE;
   }
 
   basis->coefficient   = (double *) qmckl_malloc (context, prim_num  * sizeof(double ));
   if (basis->coefficient == NULL) {
     qmckl_free(context, basis->exponent);
+    basis->exponent = NULL;
     qmckl_free(context, basis->shell_factor);
+    basis->shell_factor = NULL;
     qmckl_free(context, basis->shell_prim_num);
+    basis->shell_prim_num = NULL;
     qmckl_free(context, basis->shell_ang_mom);
+    basis->shell_ang_mom = NULL;
     qmckl_free(context, basis->shell_center);
+    basis->shell_center = NULL;
     qmckl_free(context, basis);
+    basis = NULL;
     return QMCKL_FAILURE;
   }
 
diff --git a/src/qmckl_error.org b/src/qmckl_error.org
index 3a7981a..5ee37b5 100644
--- a/src/qmckl_error.org
+++ b/src/qmckl_error.org
@@ -19,6 +19,10 @@
 MunitResult test_<<filename()>>() {
    #+end_src
 
+   #+begin_src c :comments org :tangle (eval h)
+#include <errno.h>
+#include <string.h>
+   #+end_src
 * 
 :PROPERTIES:
 :UNNUMBERED: t
diff --git a/src/test_qmckl.org b/src/test_qmckl.org
index 5633963..fe8e3a4 100644
--- a/src/test_qmckl.org
+++ b/src/test_qmckl.org
@@ -8,70 +8,81 @@
 
   First, we use a script to find the list of all the generated test files:
   #+NAME: test-files
-  #+begin_src sh :exports none :results value
-grep begin_src *.org | \
-    grep test_qmckl_ | \
-    rev | \
-    cut -d  ' ' -f 1 | \
-    rev | \
-    sort | \
-    uniq
+  #+begin_src sh :exports none 
+FILES=$(cat table_of_contents)
+grep begin_src $FILES \
+    | grep c_test  \
+    | cut -d  '.' -f 1 \
+    | uniq
   #+end_src
 
   #+RESULTS: test-files
-  | test_qmckl_ao.c       |
-  | test_qmckl_context.c  |
-  | test_qmckl_distance.c |
-  | test_qmckl_memory.c   |
+  | qmckl_error    |
+  | qmckl_context  |
+  | qmckl_memory   |
+  | qmckl_distance |
+  | qmckl_ao       |
 
   We generate the function headers
-  #+begin_src sh :var files=test-files :exports output :results raw
+  #+begin_src sh :var files=test-files :exports output :results drawer
 echo "#+NAME: headers"
 echo "#+begin_src c :tangle no"
 for file in $files
 do
-    routine=${file%.c}
+    routine=test_${file%.c}
     echo "MunitResult ${routine}();"
 done
 echo "#+end_src"
   #+end_src
 
   #+RESULTS:
+  :results:
   #+NAME: headers
   #+begin_src c :tangle no
-MunitResult test_qmckl_ao();
-MunitResult test_qmckl_context();
-MunitResult test_qmckl_distance();
-MunitResult test_qmckl_memory();
+  MunitResult test_qmckl_error();
+  MunitResult test_qmckl_context();
+  MunitResult test_qmckl_memory();
+  MunitResult test_qmckl_distance();
+  MunitResult test_qmckl_ao();
   #+end_src
-
+  :end:
+  
   and the required function calls:
-  #+begin_src sh :var files=test-files :exports output :results raw
+  #+begin_src sh :var files=test-files :exports output :results drawer
 echo "#+NAME: calls"
 echo "#+begin_src c :tangle no"
 for file in $files
 do
-    routine=${file%.c}
+    routine=test_${file%.c}
     echo "  { (char*) \"${routine}\", ${routine}, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},"
 done
 echo "#+end_src"
   #+end_src
   
   #+RESULTS:
+  :results:
   #+NAME: calls
   #+begin_src c :tangle no
-  { (char*) "test_qmckl_ao", test_qmckl_ao, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},
-  { (char*) "test_qmckl_context", test_qmckl_context, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},
-  { (char*) "test_qmckl_distance", test_qmckl_distance, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},
-  { (char*) "test_qmckl_memory", test_qmckl_memory, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},
+    { (char*) "test_qmckl_error", test_qmckl_error, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},
+    { (char*) "test_qmckl_context", test_qmckl_context, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},
+    { (char*) "test_qmckl_memory", test_qmckl_memory, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},
+    { (char*) "test_qmckl_distance", test_qmckl_distance, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},
+    { (char*) "test_qmckl_ao", test_qmckl_ao, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},
   #+end_src
+  :end:
+
+  We include the =mcheck.h= header to enable the debugging of
+  allocations with ~mtrace~. Memory allocations will be traced in the
+  file specified by the ~MALLOC_TRACE~ environment variable.
 
   #+begin_src c :comments link :noweb yes :tangle test_qmckl.c  
 #include "qmckl.h"
 #include "munit.h"
+#include "mcheck.h"
 <<headers>>
 
 int main(int argc, char* argv[MUNIT_ARRAY_PARAM(argc + 1)]) {
+  mtrace();
   static MunitTest test_suite_tests[] =
     {
      <<calls>>
@@ -83,6 +94,9 @@ int main(int argc, char* argv[MUNIT_ARRAY_PARAM(argc + 1)]) {
      (char*) "", test_suite_tests, NULL, 1, MUNIT_SUITE_OPTION_NONE
     };
 
-  return munit_suite_main(&test_suite, (void*) "µnit", argc, argv);
+  
+  int result = munit_suite_main(&test_suite, (void*) "µnit", argc, argv);
+  muntrace();
+  return result;
   }
   #+end_src
diff --git a/tools/Building.org b/tools/Building.org
index be6e8ab..8443b82 100644
--- a/tools/Building.org
+++ b/tools/Building.org
@@ -18,10 +18,10 @@ fi
 
   #+NAME: url-issues
   : https://github.com/trex-coe/qmckl/issues
-  
+
   #+NAME: url-web
   : https://trex-coe.github.io/qmckl
-  
+
   #+NAME: license
   #+begin_example
 BSD 3-Clause License
@@ -66,18 +66,18 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
   generated by the script detailed in the next section.
 ** Header                                                          :noexport:
 
- #+begin_src makefile 
+ #+begin_src makefile
 # <<header()>>
  #+end_src
 ** Dependencies
 
-   #+begin_src makefile 
+   #+begin_src makefile
 LIBS=-lpthread
    #+end_src
-   
+
 ** Variables
 
-   #+begin_src makefile 
+   #+begin_src makefile
 QMCKL_ROOT=$(shell dirname $(CURDIR))
 
 export CC CFLAGS FC FFLAGS LIBS QMCKL_ROOT
@@ -88,10 +88,10 @@ INCLUDE=-I$(QMCKL_ROOT)/include/
    #+end_src
 
 ** Compiler options
-  
+
  GNU, Intel and LLVM compilers are supported. Choose here:
 
- #+begin_src makefile 
+ #+begin_src makefile
 COMPILER=GNU
 #COMPILER=INTEL
 #COMPILER=LLVM
@@ -99,7 +99,7 @@ COMPILER=GNU
 
 *** GNU
 
-    #+begin_src makefile 
+    #+begin_src makefile
 ifeq ($(COMPILER),GNU)
 #----------------------------------------------------------
 CC=gcc -g
@@ -111,15 +111,15 @@ FFLAGS=-fPIC  $(INCLUDE) \
 	-fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising \
 	-Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation \
 	-Wreal-q-constant -Wuninitialized  -fbacktrace -finit-real=nan \
-	-ffpe-trap=zero,overflow,underflow 
+	-ffpe-trap=zero,overflow,underflow
 
-LIBS+=-lgfortran -lm 
+LIBS+=-lgfortran -lm
 #----------------------------------------------------------
 endif
     #+end_src
 
 *** Intel
-   
+
     #+begin_src makefile
 ifeq ($(COMPILER),INTEL)
 #----------------------------------------------------------
@@ -136,7 +136,7 @@ endif
     #+end_src
 
 *** LLVM
-   
+
     #+begin_src makefile
 ifeq ($(COMPILER),LLVM)
 #----------------------------------------------------------
@@ -150,7 +150,7 @@ LIBS+=-lm
 #----------------------------------------------------------
 endif
     #+end_src
-         
+
 ** Rules
 
    The source files are created during the generation of the file ~Makefile.generated~.
@@ -159,7 +159,7 @@ endif
 .PHONY: clean
 .SECONDARY: # Needed to keep the produced C and Fortran files
 
-libqmckl.so: Makefile.generated 
+libqmckl.so: Makefile.generated
 	$(MAKE) -f Makefile.generated
 
 ../include/qmckl.h: libqmckl.so
@@ -169,7 +169,7 @@ test: Makefile.generated ../include/qmckl.h
 	$(MAKE) -f Makefile.generated test
 
 doc: $(ORG_SOURCE_FILES)
-	$(QMCKL_ROOT)/tools/create_doc.sh 
+	$(QMCKL_ROOT)/tools/create_doc.sh
 
 clean:
 	$(RM)	test_qmckl_* test_qmckl.c test_qmckl \
@@ -177,7 +177,7 @@ clean:
 		Makefile.generated libqmckl.so *.html *.mod
 
 Makefile.generated: Makefile $(QMCKL_ROOT)/tools/create_makefile.sh  $(ORG_SOURCE_FILES)
-	$(QMCKL_ROOT)/tools/create_makefile.sh 
+	$(QMCKL_ROOT)/tools/create_makefile.sh
    #+end_src
 
 * Script to tangle the org-mode files
@@ -227,8 +227,8 @@ done
   This script generates the Makefile that compiles the library.
   The ~OUTPUT~ variable contains the name of the generated Makefile,typically
   =Makefile.generated=.
-  
-  #+begin_src bash 
+
+  #+begin_src bash
 # <<header()>>
 
 <<check_src>>
@@ -238,21 +238,21 @@ OUTPUT=Makefile.generated
 
   We start by tangling all the org-mode files.
 
-  #+begin_src bash 
+  #+begin_src bash
 ${QMCKL_ROOT}/tools/tangle.sh *.org
   #+end_src
 
   Then we create the list of ~*.o~ files to be created, for library
   functions:
 
-  #+begin_src bash 
-OBJECTS=""
+  #+begin_src bash
+OBJECTS="qmckl_f.o"
 for i in $(ls qmckl_*.c) ; do
     FILE=${i%.c}
     OBJECTS="${OBJECTS} ${FILE}.o"
 done >> $OUTPUT
 
-for i in $(ls qmckl_*_f.f90) ; do
+for i in $(ls qmckl*_f.f90) ; do
     FILE=${i%.f90}
     OBJECTS="${OBJECTS} ${FILE}.o"
 done >> $OUTPUT
@@ -260,7 +260,7 @@ done >> $OUTPUT
 
   for tests in C:
 
-  #+begin_src bash 
+  #+begin_src bash
 TESTS=""
 for i in $(ls test_qmckl_*.c) ; do
     FILE=${i%.c}.o
@@ -269,8 +269,8 @@ done >> $OUTPUT
   #+end_src
 
   and for tests in Fortran:
-  
-  #+begin_src bash 
+
+  #+begin_src bash
 TESTS_F=""
 for i in $(ls test_qmckl_*_f.f90) ; do
     FILE=${i%.f90}.o
@@ -280,10 +280,10 @@ done >> $OUTPUT
 
   Finally, we append the rules to the Makefile
 
-  #+begin_src bash 
+  #+begin_src bash
 cat << EOF > ${OUTPUT}
 CC=$CC
-CFLAGS=$CFLAGS -I../munit/ 
+CFLAGS=$CFLAGS -I../munit/
 
 FC=$FC
 FFLAGS=$FFLAGS
@@ -295,17 +295,18 @@ LIBS=$LIBS
 
 libqmckl.so: \$(OBJECT_FILES)
 	\$(CC) -shared \$(OBJECT_FILES) -o libqmckl.so
-           
-%.o: %.c 
+
+%.o: %.c
 	\$(CC) \$(CFLAGS) -c \$*.c -o \$*.o
 
-
-qmckl_f.o: ../include/qmckl_f.f90
-	\$(FC) \$(FFLAGS) -c ../include/qmckl_f.f90 -o qmckl_f.o
-
 %.o: %.f90 qmckl_f.o
 	\$(FC) \$(FFLAGS) -c \$*.f90 -o \$*.o
 
+../include/qmckl.h ../include/qmckl_f.f90:
+	../tools/build_qmckl_h.sh
+
+qmckl_f.o: ../include/qmckl_f.f90
+	\$(FC) \$(FFLAGS) -c ../include/qmckl_f.f90 -o qmckl_f.o
 
 test_qmckl: test_qmckl.c libqmckl.so \$(TESTS) \$(TESTS_F)
 	\$(CC) \$(CFLAGS) -Wl,-rpath,$PWD -L. \
@@ -317,8 +318,8 @@ test: test_qmckl
 .PHONY: test
 EOF
 
-  #+end_src 
-  
+  #+end_src
+
 * Script to build the final qmckl.h file
   :PROPERTIES:
   :header-args:bash: :tangle build_qmckl_h.sh :noweb  yes :shebang #!/bin/bash :comments org
@@ -339,7 +340,7 @@ EOF
  Issues        : <<url-issues()>>
 
  <<license()>>
- 
+
 
   #+end_src
 
@@ -349,7 +350,7 @@ EOF
 
   Put =.h= files in the correct order:
 
-  #+begin_src bash 
+  #+begin_src bash
 HEADERS=""
 for i in $(cat table_of_contents)
 do
@@ -358,8 +359,8 @@ done
   #+end_src
 
   Generate C header file
-  
-  #+begin_src bash 
+
+  #+begin_src bash
 OUTPUT="../include/qmckl.h"
 
 cat << EOF > ${OUTPUT}
@@ -387,8 +388,8 @@ EOF
   #+end_src
 
   Generate Fortran interface file from all =qmckl_*_fh.f90= files
-  
-  #+begin_src bash 
+
+  #+begin_src bash
 HEADERS="qmckl_*_fh.f90"
 
 OUTPUT="../include/qmckl_f.f90"
diff --git a/tools/create_makefile.sh b/tools/create_makefile.sh
index d7bd2f5..e134089 100755
--- a/tools/create_makefile.sh
+++ b/tools/create_makefile.sh
@@ -28,13 +28,13 @@ ${QMCKL_ROOT}/tools/tangle.sh *.org
 # functions:
 
 
-OBJECTS=""
+OBJECTS="qmckl_f.o"
 for i in $(ls qmckl_*.c) ; do
     FILE=${i%.c}
     OBJECTS="${OBJECTS} ${FILE}.o"
 done >> $OUTPUT
 
-for i in $(ls qmckl_*_f.f90) ; do
+for i in $(ls qmckl*_f.f90) ; do
     FILE=${i%.f90}
     OBJECTS="${OBJECTS} ${FILE}.o"
 done >> $OUTPUT
@@ -84,13 +84,14 @@ libqmckl.so: \$(OBJECT_FILES)
 %.o: %.c 
 	\$(CC) \$(CFLAGS) -c \$*.c -o \$*.o
 
-
-qmckl_f.o: ../include/qmckl_f.f90
-	\$(FC) \$(FFLAGS) -c ../include/qmckl_f.f90 -o qmckl_f.o
-
 %.o: %.f90 qmckl_f.o
 	\$(FC) \$(FFLAGS) -c \$*.f90 -o \$*.o
 
+../include/qmckl.h ../include/qmckl_f.f90:
+	../tools/build_qmckl_h.sh
+
+qmckl_f.o: ../include/qmckl_f.f90
+	\$(FC) \$(FFLAGS) -c ../include/qmckl_f.f90 -o qmckl_f.o
 
 test_qmckl: test_qmckl.c libqmckl.so \$(TESTS) \$(TESTS_F)
 	\$(CC) \$(CFLAGS) -Wl,-rpath,$PWD -L. \