diff --git a/configure.ac b/configure.ac
index d99f4a0..f6839d8 100644
--- a/configure.ac
+++ b/configure.ac
@@ -222,8 +222,9 @@ AC_ARG_ENABLE(debug, [AS_HELP_STRING([--enable-debug],[compile for debugging])],
 if test "$ok" = "yes"; then
         if test "$GCC" = "yes"; then
 	        CFLAGS="$CFLAGS \
--Wall -W -Wbad-function-cast -Wcast-qual \
--Wpointer-arith -Wcast-align -Wpedantic -Wextra -fmax-errors=3"
+-Wall -W -Wbad-function-cast -Wcast-qual -Warray-bounds -Wdisabled-optimization \
+-Wpointer-arith -Wcast-align -Wpedantic -Wextra -fmax-errors=5 -Walloc-zero -Werror \
+"
         fi
         if test "$GFC" = "yes"; then
 	        FCFLAGS="$FCFLAGS \
diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org
index 07d27ba..5e5a99a 100644
--- a/org/qmckl_ao.org
+++ b/org/qmckl_ao.org
@@ -73,6 +73,8 @@ gradients and Laplacian of the atomic basis functions.
 
 #include <stdio.h>
 #include <math.h>
+#include <string.h>
+
 #include "chbrclf.h"
 #include "qmckl_ao_private_func.h"
 
@@ -2399,6 +2401,7 @@ assert (rc == QMCKL_SUCCESS);
 for (int64_t i=0 ; i < shell_num ; ++i) {
   assert(shell_ang_mom_test[i] == shell_ang_mom[i]);
  }
+free(shell_ang_mom_test);
 
 shell_factor_test = (double*) malloc ( shell_num * sizeof(double));
 rc = qmckl_get_ao_basis_shell_factor  (context, shell_factor_test, shell_num);
@@ -2407,6 +2410,7 @@ assert (rc == QMCKL_SUCCESS);
 for (int64_t i=0 ; i < shell_num ; ++i) {
   assert(shell_factor_test[i] == shell_factor[i]);
 }
+free(shell_factor_test);
 
 shell_prim_num_test = (int64_t*) malloc ( shell_num * sizeof(int64_t));
 rc = qmckl_get_ao_basis_shell_prim_num (context, shell_prim_num_test, shell_num);
@@ -2415,6 +2419,7 @@ assert (rc == QMCKL_SUCCESS);
 for (int64_t i=0 ; i < shell_num ; ++i) {
   assert(shell_prim_num_test[i] == shell_prim_num[i]);
 }
+free(shell_prim_num_test);
 
 shell_prim_index_test = (int64_t*) malloc ( shell_num * sizeof(int64_t));
 rc = qmckl_get_ao_basis_shell_prim_index (context, shell_prim_index_test, shell_num);
@@ -2423,6 +2428,7 @@ assert (rc == QMCKL_SUCCESS);
 for (int64_t i=0 ; i < shell_num ; ++i) {
   assert(shell_prim_index_test[i] == shell_prim_index[i]);
 }
+free(shell_prim_index_test);
 
 exponent_test = (double*) malloc ( prim_num * sizeof(double));
 rc = qmckl_get_ao_basis_exponent(context, exponent_test, prim_num);
@@ -2431,6 +2437,7 @@ assert (rc == QMCKL_SUCCESS);
 for (int64_t i=0 ; i < prim_num ; ++i) {
   assert(exponent_test[i] == exponent[i]);
 }
+free(exponent_test);
 
 coefficient_test = (double*) malloc ( prim_num * sizeof(double));
 rc = qmckl_get_ao_basis_coefficient(context, coefficient_test, prim_num);
@@ -2439,6 +2446,7 @@ assert (rc == QMCKL_SUCCESS);
 for (int64_t i=0 ; i < prim_num ; ++i) {
   assert(coefficient_test[i] == coefficient[i]);
 }
+free(coefficient_test);
 
 prim_factor_test = (double*) malloc ( prim_num * sizeof(double));
 rc = qmckl_get_ao_basis_prim_factor (context, prim_factor_test, prim_num);
@@ -2447,6 +2455,7 @@ assert (rc == QMCKL_SUCCESS);
 for (int64_t i=0 ; i < prim_num ; ++i) {
   assert(prim_factor_test[i] == prim_factor[i]);
 }
+free(prim_factor_test);
 
 rc = qmckl_get_ao_basis_ao_num(context, &ao_num_test);
 assert(ao_num == ao_num_test);
@@ -2458,6 +2467,7 @@ assert (rc == QMCKL_SUCCESS);
 for (int64_t i=0 ; i < ao_num ; ++i) {
   assert(ao_factor_test[i] == ao_factor[i]);
 }
+free(ao_factor_test);
 
     #+end_src
 
@@ -4114,7 +4124,6 @@ integer function qmckl_ao_polynomial_vgl_doc_f (context,    &
   integer*8         :: i,j
   integer           :: a,b,c,d
   real*8            :: Y(3)
-  integer           :: lmax_array(3)
   real*8            :: pows(-2:lmax,3)
   double precision  :: xy, yz, xz
   double precision  :: da, db, dc, dd
@@ -4146,7 +4155,6 @@ integer function qmckl_ao_polynomial_vgl_doc_f (context,    &
      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
@@ -4339,6 +4347,22 @@ end function qmckl_ao_polynomial_vgl_doc_f
           const int64_t ldv ); 
     #+end_src
 
+    #+CALL: generate_c_header(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_transp_vgl_hpc")
+
+    #+RESULTS:
+    #+begin_src c :tangle (eval h_func) :comments org
+    qmckl_exit_code qmckl_ao_polynomial_transp_vgl_hpc (
+          const qmckl_context context,
+          const double* X,
+          const double* R,
+          const int32_t lmax,
+          int64_t* n,
+          int32_t* const L,
+          const int64_t ldl,
+          double* const VGL,
+          const int64_t ldv ); 
+    #+end_src
+
     #+begin_src c :tangle (eval c) :comments org
 qmckl_exit_code
 qmckl_ao_polynomial_transp_vgl (const qmckl_context context,
@@ -4352,8 +4376,7 @@ qmckl_ao_polynomial_transp_vgl (const qmckl_context context,
                                 const int64_t ldv )
 {
 #ifdef HAVE_HPC
-  //return qmckl_ao_polynomial_transp_vgl_hpc (context, X, R, lmax, n, L, ldl, VGL, ldv);
-  return qmckl_ao_polynomial_transp_vgl_doc (context, X, R, lmax, n, L, ldl, VGL, ldv);
+  return qmckl_ao_polynomial_transp_vgl_hpc (context, X, R, lmax, n, L, ldl, VGL, ldv);
 #else
   return qmckl_ao_polynomial_transp_vgl_doc (context, X, R, lmax, n, L, ldl, VGL, ldv);
 #endif  
@@ -4377,7 +4400,6 @@ integer function qmckl_ao_polynomial_transp_vgl_doc_f (context,    &
   integer*8         :: i,j
   integer           :: a,b,c,d
   real*8            :: Y(3)
-  integer           :: lmax_array(3)
   real*8            :: pows(-2:21,3) ! lmax < 22
   double precision  :: xy, yz, xz
   double precision  :: da, db, dc, dd
@@ -4405,7 +4427,6 @@ integer function qmckl_ao_polynomial_transp_vgl_doc_f (context,    &
   endif
 
 
-  lmax_array(1:3) = lmax
   if (lmax > 0) then
 
      do i=1,3
@@ -4422,18 +4443,17 @@ integer function qmckl_ao_polynomial_transp_vgl_doc_f (context,    &
      VGL(1:4,1:5) = 0.d0
 
      VGL(1  ,1  ) = 1.d0
-     VGL(2:4,1:5) = 0.d0
 
      l  (1,2) = 1
-     VGL(2,1) = pows(1,1)
+     VGL(2,1) = Y(1)
      VGL(2,2) = 1.d0
 
      l  (2,3) = 1
-     VGL(3,1) = pows(1,2)
+     VGL(3,1) = Y(2)
      VGL(3,3) = 1.d0
 
      l  (3,4) = 1
-     VGL(4,1) = pows(1,3)
+     VGL(4,1) = Y(3)
      VGL(4,4) = 1.d0
 
      n=4
@@ -4456,14 +4476,14 @@ integer function qmckl_ao_polynomial_transp_vgl_doc_f (context,    &
            dc = dd - da - db
            n = n+1
 
-           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)
 
+           l(1,n) = a
+           l(2,n) = b
+           l(3,n) = c
+
            VGL(n,1) = xy * pows(c,3)
 
            xy = dc * xy
@@ -4491,6 +4511,264 @@ integer function qmckl_ao_polynomial_transp_vgl_doc_f (context,    &
 end function qmckl_ao_polynomial_transp_vgl_doc_f
     #+end_src
 
+    #+begin_src c :tangle (eval c) :comments org
+qmckl_exit_code
+qmckl_ao_polynomial_transp_vgl_hpc (const qmckl_context context,
+                                    const double* X,
+                                    const double* R,
+                                    const int32_t lmax,
+                                    int64_t* n,
+                                    int32_t* const L,
+                                    const int64_t ldl,
+                                    double* const VGL,
+                                    const int64_t ldv )
+{
+  const qmckl_context_struct* ctx = (qmckl_context_struct* const) context;
+  assert (ctx != NULL && X != NULL && R != NULL && n != NULL && L != NULL && VGL != NULL);
+  if (lmax < 0) return QMCKL_INVALID_ARG_4;
+  if (ldl < 3) return QMCKL_INVALID_ARG_7;
+
+  int32_t m;
+
+  double* const vgl1 = VGL;
+  double* const vgl2 = VGL + ldv;
+  double* const vgl3 = VGL + (ldv << 1);
+  double* const vgl4 = VGL + 3*ldv;
+  double* const vgl5 = VGL + (ldv << 2);
+
+  switch (lmax) {
+  case 0:
+    {
+      if (ldv < 1) return QMCKL_INVALID_ARG_9;
+      L[0] = 0;
+      L[ldl] = 0;
+      L[ldl << 1] = 0;
+      m=1;
+
+      vgl1[0] = 1.0;
+      vgl2[0] = 0.0;
+      vgl3[0] = 0.0;
+      vgl4[0] = 0.0;
+      vgl5[0] = 0.0;
+      break;
+    }
+  case 1:
+    {
+      if (ldv < 4) return QMCKL_INVALID_ARG_9;
+      const double Y[3] = { X[0]-R[0], 
+                            X[1]-R[1], 
+                            X[2]-R[2] };
+
+      int32_t* const l[4] = {L, L+ldl, L+(ldl<<1), L+3*ldl};
+      for (int32_t i=0 ; i<4 ; ++i) {
+        l[i][0] = 0;
+        l[i][1] = 0;
+        l[i][2] = 0;
+      }
+      l[1][0] = 1;
+      l[2][1] = 1;
+      l[3][2] = 1;
+
+      for (int32_t k=0 ; k<4 ; ++k) {
+        vgl2[k] = 0.0;
+        vgl3[k] = 0.0;
+        vgl4[k] = 0.0;
+        vgl5[k] = 0.0;
+      }
+      vgl1[0] = 1.0;
+      vgl1[1] = Y[0];
+      vgl1[2] = Y[1];
+      vgl1[3] = Y[2];
+      vgl2[1] = 1.0;
+      vgl3[2] = 1.0;
+      vgl4[3] = 1.0;
+      m=4;
+      break;
+    }
+  case 2:
+    {
+      if (ldv < 10) return QMCKL_INVALID_ARG_9;
+      const double Y[3] = { X[0]-R[0], 
+                            X[1]-R[1], 
+                            X[2]-R[2] };
+
+      double pows[3][5] = { {1.0, 1.0, 1.0, Y[0], Y[0]*Y[0]},
+                            {1.0, 1.0, 1.0, Y[1], Y[1]*Y[1]},
+                            {1.0, 1.0, 1.0, Y[2], Y[2]*Y[2]} };
+      
+      int32_t* l[10];
+      for (int32_t i=0 ; i<10 ; ++i) {
+        l[i] = L + i*ldl;
+      }
+      
+      for (int32_t i=0 ; i<4 ; ++i) {
+        l[i][0] = 0;
+        l[i][1] = 0;
+        l[i][2] = 0;
+      }
+      l[1][0] = 1;
+      l[2][1] = 1;
+      l[3][2] = 1;
+      
+      for (int32_t k=0 ; k<4 ; ++k) {
+        vgl2[k] = 0.0;
+        vgl3[k] = 0.0;
+        vgl4[k] = 0.0;
+        vgl5[k] = 0.0;
+      }
+      vgl1[0] = 1.0;
+      vgl1[1] = Y[0];
+      vgl1[2] = Y[1];
+      vgl1[3] = Y[2];
+
+      vgl2[1] = 1.0;
+      vgl3[2] = 1.0;
+      vgl4[3] = 1.0;
+      m=4;
+      
+      double da = 2.0;
+        
+      for (int32_t a=2 ; a>=0 ; --a) {
+        double db = 2.0-da;
+        
+        for (int32_t b=2-a ; b>=0 ; --b) {
+          const int32_t c  = 2 - a - b;  
+          const double dc = 2.0 - da - db;
+          
+          double xy = pows[0][a+2] * pows[1][b+2];
+          double yz = pows[1][b+2] * pows[2][c+2];
+          double xz = pows[0][a+2] * pows[2][c+2];
+          
+          l[m][0] = a;
+          l[m][1] = b;
+          l[m][2] = c;
+          
+          vgl1[m] = xy * pows[2][c+2];
+          
+          xy *= dc;
+          xz *= db;
+          yz *= da;
+          
+          vgl2[m] = pows[0][a+1] * yz;
+          vgl3[m] = pows[1][b+1] * xz;
+          vgl4[m] = pows[2][c+1] * xy;
+          
+          vgl5[m] = (da-1.) * pows[0][a] * yz +
+            (db-1.) * pows[1][b] * xz +
+            (dc-1.) * pows[2][c] * xy;
+          
+          db -= 1.0;
+          ++m;
+        }
+        da -= 1.0;
+      }
+      break;
+    }
+    /*
+        case 3: 
+        {
+        if (ldv < 20) return QMCKL_INVALID_ARG_9;
+        }
+      ,*/
+  default: 
+    {
+      const int32_t size_max = (lmax+1)*(lmax+2)*(lmax+3)/6;
+      if (ldv < size_max) return QMCKL_INVALID_ARG_9;
+      const double Y[3] = { X[0]-R[0], 
+                            X[1]-R[1], 
+                            X[2]-R[2] };
+      double pows[3][size_max];
+      
+      for (int32_t i=0 ; i<=2 ; ++i) {
+        pows[0][i] = 1.0;
+        pows[1][i] = 1.0;
+        pows[2][i] = 1.0;
+      }
+      
+      for (int32_t i=3 ; i<=lmax+2 ; ++i) {
+        pows[0][i] = pows[0][i-1] * Y[0];
+        pows[1][i] = pows[1][i-1] * Y[1];
+        pows[2][i] = pows[2][i-1] * Y[2];
+      }
+      
+      int32_t* l[24];
+      for (int32_t i=0 ; i<size_max ; ++i) {
+        l[i] = L + i*ldl;
+      }
+      
+      for (int32_t i=0 ; i<4 ; ++i) {
+        l[i][0] = 0;
+        l[i][1] = 0;
+        l[i][2] = 0;
+      }
+      l[1][0] = 1;
+      l[2][1] = 1;
+      l[3][2] = 1;
+      
+      for (int32_t k=0 ; k<4 ; ++k) {
+        vgl2[k] = 0.0;
+        vgl3[k] = 0.0;
+        vgl4[k] = 0.0;
+        vgl5[k] = 0.0;
+      }
+      vgl1[0] = 1.0;
+      vgl1[1] = Y[0];
+      vgl1[2] = Y[1];
+      vgl1[3] = Y[2];
+
+      vgl2[1] = 1.0;
+      vgl3[2] = 1.0;
+      vgl4[3] = 1.0;
+      m=4;
+      
+      double dd = 2.0;
+      
+      for (int32_t d=2 ; d<= lmax ; ++d) {
+        double da = dd;
+        
+        for (int32_t a=d ; a>=0 ; --a) {
+          double db = dd-da;
+          
+          for (int32_t b=d-a ; b>=0 ; --b) {
+            const int32_t c  = d - a - b;  
+            const double dc = dd - da - db;
+            
+            double xy = pows[0][a+2] * pows[1][b+2];
+            double yz = pows[1][b+2] * pows[2][c+2];
+            double xz = pows[0][a+2] * pows[2][c+2];
+            
+            l[m][0] = a;
+            l[m][1] = b;
+            l[m][2] = c;
+            
+            vgl1[m] = xy * pows[2][c+2];
+            
+            xy *= dc;
+            xz *= db;
+            yz *= da;
+            
+            vgl2[m] = pows[0][a+1] * yz;
+            vgl3[m] = pows[1][b+1] * xz;
+            vgl4[m] = pows[2][c+1] * xy;
+            
+            vgl5[m] = (da-1.) * pows[0][a] * yz +
+              (db-1.) * pows[1][b] * xz +
+              (dc-1.) * pows[2][c] * xy;
+            
+            db -= 1.0;
+            ++m;
+          }
+          da -= 1.0;
+        }
+        dd += 1.0;
+      }
+    }
+  }
+  ,*n = m;
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
     #+CALL: generate_c_interface(table=qmckl_ao_polynomial_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_ao_polynomial_transp_vgl_doc")
 
     #+RESULTS:
@@ -4666,8 +4944,42 @@ end function test_qmckl_ao_polynomial_vgl
     #+end_src
 
     #+begin_src c :tangle (eval c_test)
-    int  test_qmckl_ao_polynomial_vgl(qmckl_context context);
-    assert(0 == test_qmckl_ao_polynomial_vgl(context));
+int  test_qmckl_ao_polynomial_vgl(qmckl_context context);
+assert(0 == test_qmckl_ao_polynomial_vgl(context));
+
+double X[3] = { 1.1, 2.2, 3.3 };
+double R[3] = { 0.2, 1.1, 3.0 };
+int64_t n;
+int64_t ldl = 4;
+int64_t ldv = 24;
+int32_t L0[24][ldl];
+int32_t L1[24][ldl];
+double VGL0[5][ldv];
+double VGL1[5][ldv];
+memset(&L0[0][0], 0, sizeof(L0));
+memset(&L1[0][0], 0, sizeof(L1));
+memset(&VGL0[0][0], 0, sizeof(VGL0));
+memset(&VGL1[0][0], 0, sizeof(VGL1));
+for (int32_t lmax=0 ; lmax<=3 ; lmax++) {
+rc = qmckl_ao_polynomial_transp_vgl_doc (context, X, R, lmax, &n, &(L0[0][0]), ldl, &(VGL0[0][0]), ldv);
+assert (rc == QMCKL_SUCCESS);
+rc = qmckl_ao_polynomial_transp_vgl_hpc (context, X, R, lmax, &n, &(L1[0][0]), ldl, &(VGL1[0][0]), ldv);
+assert (rc == QMCKL_SUCCESS);
+printf("lmax=%d\n", lmax);
+for (int32_t l=0 ; l<n ; ++l) {
+  for (int32_t k=0 ; k<3 ; ++k) {
+    printf("L[%d][%d] = %d   %d\n", l, k, L0[l][k], L1[l][k]);
+    assert( L0[l][k] == L1[l][k] );
+  }
+ }
+
+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] );
+  }
+ }
+}
     #+end_src
 
 * Combining radial and polynomial parts
@@ -4979,7 +5291,6 @@ qmckl_compute_ao_vgl_hpc (
              const double s3 = s6_*y;
              const double s4 = s6_*z;
              const double s5 = s5_;
-             const double s6 = s6_;
 
              const int32_t l = shell_ang_mom[ishell];
              const int32_t n = lstart[l+1]-lstart[l];
diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org
index 998e6ef..274d269 100644
--- a/org/qmckl_blas.org
+++ b/org/qmckl_blas.org
@@ -111,26 +111,26 @@ qmckl_vector_alloc( qmckl_context context,
    #+begin_src c :comments org :tangle (eval h_private_func)
 qmckl_exit_code
 qmckl_vector_free( qmckl_context context, 
-                   qmckl_vector vector);
+                   qmckl_vector* vector);
    #+end_src
 
    #+begin_src c :comments org :tangle (eval c)
 qmckl_exit_code
 qmckl_vector_free( qmckl_context context, 
-                   qmckl_vector vector)
+                   qmckl_vector* vector)
 {
   /* Always true */
-  assert (vector.data != NULL);
+  assert (vector->data != NULL);
   
   qmckl_exit_code rc;
   
-  rc = qmckl_free(context, vector.data);
+  rc = qmckl_free(context, vector->data);
   if (rc != QMCKL_SUCCESS) {
     return rc;
   }
 
-  vector.size = (int64_t) 0;
-  vector.data = NULL;
+  vector->size = (int64_t) 0;
+  vector->data = NULL;
   return QMCKL_SUCCESS;
 }
    #+end_src
@@ -192,26 +192,26 @@ qmckl_matrix_alloc( qmckl_context context,
    #+begin_src c :comments org :tangle (eval h_private_func)
 qmckl_exit_code
 qmckl_matrix_free( qmckl_context context, 
-                   qmckl_matrix matrix);
+                   qmckl_matrix* matrix);
    #+end_src
 
    #+begin_src c :comments org :tangle (eval c)
 qmckl_exit_code
 qmckl_matrix_free( qmckl_context context, 
-                   qmckl_matrix matrix)
+                   qmckl_matrix* matrix)
 {
   /* Always true */
-  assert (matrix.data != NULL);
+  assert (matrix->data != NULL);
   
   qmckl_exit_code rc;
   
-  rc = qmckl_free(context, matrix.data);
+  rc = qmckl_free(context, matrix->data);
   if (rc != QMCKL_SUCCESS) {
     return rc;
   }
-  matrix.data = NULL;
-  matrix.size[0] = (int64_t) 0;
-  matrix.size[1] = (int64_t) 0;
+  matrix->data = NULL;
+  matrix->size[0] = (int64_t) 0;
+  matrix->size[1] = (int64_t) 0;
 
   return QMCKL_SUCCESS;
 }
@@ -286,25 +286,25 @@ 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 tensor);
+                   qmckl_tensor* tensor);
    #+end_src
 
    #+begin_src c :comments org :tangle (eval c)
 qmckl_exit_code
 qmckl_tensor_free( qmckl_context context, 
-                   qmckl_tensor tensor)
+                   qmckl_tensor* tensor)
 {
   /* Always true */
-  assert (tensor.data != NULL);
+  assert (tensor->data != NULL);
   
   qmckl_exit_code rc;
   
-  rc = qmckl_free(context, tensor.data);
+  rc = qmckl_free(context, tensor->data);
   if (rc != QMCKL_SUCCESS) {
     return rc;
   }
   
-  memset(&tensor, 0, sizeof(qmckl_tensor));
+  memset(tensor, 0, sizeof(qmckl_tensor));
 
   return QMCKL_SUCCESS;
 }
@@ -712,7 +712,7 @@ qmckl_tensor_of_double(const qmckl_context context,
   for (int64_t i=0 ; i<p ; ++i) 
       assert ( qmckl_vec(vec2, i) == qmckl_vec(vec, i) ) ;
 
-  qmckl_vector_free(context, vec);
+  qmckl_vector_free(context, &vec);
 
 }
     #+end_src
@@ -2121,8 +2121,8 @@ qmckl_transpose (qmckl_context context,
     for (int i=0 ; i<2 ; ++i)
       assert (qmckl_mat(A, i, j) == qmckl_mat(At, j, i));
 
-  qmckl_matrix_free(context, A);
-  qmckl_matrix_free(context, At);
+  qmckl_matrix_free(context, &A);
+  qmckl_matrix_free(context, &At);
 }
 
      #+end_src
diff --git a/org/qmckl_determinant.org b/org/qmckl_determinant.org
index f7ca97b..a28f49c 100644
--- a/org/qmckl_determinant.org
+++ b/org/qmckl_determinant.org
@@ -316,7 +316,7 @@ int64_t* qmckl_get_determinant_mo_index_alpha (const qmckl_context context) {
   int32_t mask = 1 << 4;
 
   if ( (ctx->det.uninitialized & mask) != 0) {
-    return (int64_t) 0;
+    return NULL;
   }
 
   assert (ctx->det.mo_index_alpha != NULL);
diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org
index 5e0d5dd..33ac366 100644
--- a/org/qmckl_electron.org
+++ b/org/qmckl_electron.org
@@ -499,18 +499,15 @@ ctx->electron.provided = (ctx->electron.uninitialized == 0);
 
 if (ctx->electron.provided) {
   if (ctx->electron.coord_new.data != NULL) {
-    const qmckl_exit_code rc = qmckl_matrix_free(context, ctx->electron.coord_new);
+    const qmckl_exit_code rc = qmckl_matrix_free(context, &(ctx->electron.coord_new));
     assert (rc == QMCKL_SUCCESS);
   }
   if (ctx->electron.coord_old.data != NULL) {
-    const qmckl_exit_code rc = qmckl_matrix_free(context, ctx->electron.coord_old);
+    const qmckl_exit_code rc = qmckl_matrix_free(context, &(ctx->electron.coord_old));
     assert (rc == QMCKL_SUCCESS);
   }
 
-  const int64_t walk_num = ctx->electron.walk_num;
-  const int64_t elec_num = ctx->electron.num;
-
-  const int64_t point_num = walk_num * elec_num;
+  const int64_t point_num = ctx->electron.walk_num * ctx->electron.num;
 
   qmckl_matrix coord_new = qmckl_matrix_alloc(context, point_num, 3);
 
@@ -1667,8 +1664,7 @@ qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context)
       ctx->electron.ee_pot = ee_pot;
     }
 
-    qmckl_exit_code rc =
-      qmckl_compute_ee_potential(context,
+    rc = qmckl_compute_ee_potential(context,
                                 ctx->electron.num,
                                 ctx->electron.walk_num,
                                 ctx->electron.ee_distance,
@@ -2712,8 +2708,7 @@ qmckl_exit_code qmckl_provide_en_potential(qmckl_context context)
       ctx->electron.en_pot = en_pot;
     }
 
-    qmckl_exit_code rc =
-      qmckl_compute_en_potential(context,
+    rc = qmckl_compute_en_potential(context,
                                 ctx->electron.num,
                                 ctx->nucleus.num,
                                 ctx->electron.walk_num,
diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org
index d745dcf..61062af 100644
--- a/org/qmckl_jastrow.org
+++ b/org/qmckl_jastrow.org
@@ -907,7 +907,7 @@ qmckl_set_jastrow_aord_vector(qmckl_context context,
   qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
   mem_info.size = (aord_num + 1) * type_nucl_num * sizeof(double);
 
-  if (size_max < mem_info.size/sizeof(double)) {
+  if ((size_t) size_max < mem_info.size/sizeof(double)) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
                            "qmckl_set_jastrow_aord_vector",
@@ -970,7 +970,7 @@ qmckl_set_jastrow_bord_vector(qmckl_context context,
   qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
   mem_info.size = (bord_num + 1) * sizeof(double);
 
-  if (size_max < mem_info.size/sizeof(double)) {
+  if ((size_t) size_max < mem_info.size/sizeof(double)) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
                            "qmckl_set_jastrow_bord_vector",
@@ -1040,7 +1040,7 @@ qmckl_set_jastrow_cord_vector(qmckl_context context,
   qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
   mem_info.size = dim_cord_vect * type_nucl_num * sizeof(double);
 
-  if (size_max < mem_info.size/sizeof(double)) {
+  if ((size_t) size_max < mem_info.size/sizeof(double)) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
                            "qmckl_set_jastrow_cord_vector",
@@ -1333,7 +1333,7 @@ qmckl_get_jastrow_asymp_jasb(qmckl_context context,
   qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
   assert (ctx != NULL);
 
-  size_t sze = 2;
+  int64_t sze = 2;
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
@@ -3084,7 +3084,7 @@ qmckl_get_jastrow_een_rescaled_e(qmckl_context context,
   qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
   assert (ctx != NULL);
 
-  size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
+  int64_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
@@ -3416,7 +3416,7 @@ qmckl_get_jastrow_een_rescaled_e_deriv_e(qmckl_context context,
   qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
   assert (ctx != NULL);
 
-  size_t sze = ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
+  int64_t sze = ctx->electron.num * 4 * ctx->electron.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
@@ -3795,7 +3795,7 @@ qmckl_get_jastrow_een_rescaled_n(qmckl_context context,
   qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
   assert (ctx != NULL);
 
-  size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
+  int64_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
@@ -4107,7 +4107,7 @@ qmckl_get_jastrow_een_rescaled_n_deriv_e(qmckl_context context,
   qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
   assert (ctx != NULL);
 
-  size_t sze = ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
+  int64_t sze = ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walk_num * (ctx->jastrow.cord_num + 1);
   if (size_max < sze) {
     return qmckl_failwith( context,
                            QMCKL_INVALID_ARG_3,
diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org
index 709fb31..d4efb41 100644
--- a/org/qmckl_mo.org
+++ b/org/qmckl_mo.org
@@ -104,7 +104,7 @@ typedef struct qmckl_mo_basis_struct {
   double  * coefficient;
 
   double  * mo_vgl;
-  int64_t   mo_vgl_date;
+  uint64_t  mo_vgl_date;
 
   int32_t   uninitialized;
   bool      provided;
diff --git a/org/qmckl_nucleus.org b/org/qmckl_nucleus.org
index a562618..da003b5 100644
--- a/org/qmckl_nucleus.org
+++ b/org/qmckl_nucleus.org
@@ -333,7 +333,7 @@ qmckl_get_nucleus_coord (const qmckl_context context,
     rc = qmckl_transpose(context, ctx->nucleus.coord, At);
     if (rc != QMCKL_SUCCESS) return rc;
     rc = qmckl_double_of_matrix(context, At, coord, size_max);
-    qmckl_matrix_free(context, At);
+    qmckl_matrix_free(context, &At);
   } else {
     rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, size_max);
   }
@@ -470,6 +470,12 @@ qmckl_set_nucleus_charge(qmckl_context context,
 
   ctx->nucleus.charge = qmckl_vector_alloc(context, num);
   rc = qmckl_vector_of_double(context, charge, num, &(ctx->nucleus.charge));
+  if (rc != QMCKL_SUCCESS) {
+    return qmckl_failwith( context,
+                           QMCKL_FAILURE,
+                           "qmckl_set_nucleus_charge",
+                           "Error in vector->double* conversion");
+  }
 
   <<post2>>
 }
@@ -518,7 +524,7 @@ qmckl_set_nucleus_coord(qmckl_context context,
   const int64_t nucl_num = (int64_t) ctx->nucleus.num;
 
   if (ctx->nucleus.coord.data != NULL) {
-    rc = qmckl_matrix_free(context, ctx->nucleus.coord);
+    rc = qmckl_matrix_free(context, &(ctx->nucleus.coord));
     if (rc != QMCKL_SUCCESS) return rc;
   }
 
diff --git a/org/qmckl_point.org b/org/qmckl_point.org
index d0f4614..0672d1c 100644
--- a/org/qmckl_point.org
+++ b/org/qmckl_point.org
@@ -224,7 +224,7 @@ qmckl_get_point(const qmckl_context context,
     if (rc != QMCKL_SUCCESS) return rc;
     rc = qmckl_double_of_matrix( context, At, coord, size_max);
     if (rc != QMCKL_SUCCESS) return rc;
-    rc = qmckl_matrix_free(context, At);
+    rc = qmckl_matrix_free(context, &At);
   } else {
     rc = qmckl_double_of_matrix( context, ctx->point.coord, coord, size_max);
   }
@@ -301,7 +301,7 @@ qmckl_set_point (qmckl_context context,
   if (ctx->point.num < num) {
 
     if (ctx->point.coord.data != NULL) {
-      rc = qmckl_matrix_free(context, ctx->point.coord);
+      rc = qmckl_matrix_free(context, &(ctx->point.coord));
       assert (rc == QMCKL_SUCCESS);
     }
 
diff --git a/org/qmckl_trexio.org b/org/qmckl_trexio.org
index e08f044..5afeadd 100644
--- a/org/qmckl_trexio.org
+++ b/org/qmckl_trexio.org
@@ -429,6 +429,11 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
 
     /* Reformat data */
     rc = qmckl_set_ao_basis_nucleus_index(context, nucleus_index, nucleus_num);
+    if (rc != QMCKL_SUCCESS) {
+      qmckl_free(context, nucleus_index);
+      nucleus_index = NULL;
+      return rc;
+    }
 
     for (int i=shell_num-1 ; i>=0 ; --i) {
       const int k = tmp_array[i];