diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org
index 43ca5bc..b53e4f7 100644
--- a/org/qmckl_ao.org
+++ b/org/qmckl_ao.org
@@ -3528,7 +3528,7 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f( &
   info = QMCKL_SUCCESS
 
   ! Don't compute exponentials when the result will be almost zero.
-  cutoff = -dlog(1.d-15)
+  cutoff = 27.631021115928547  ! -dlog(1.d-12)
 
   do inucl=1,nucl_num
      ! C is zero-based, so shift bounds by one
@@ -3752,10 +3752,15 @@ print ( "[7][4][26] : %e"% lf(a,x,y))
           (int64_t) 5*elec_num*walk_num*prim_num );
   assert (rc == QMCKL_SUCCESS);
 
+  printf("prim_vgl[26][0][7] = %e\n",prim_vgl[26][0][7]);
   assert( fabs(prim_vgl[26][0][7] - ( 1.0501570432064878E-003)) < 1.e-14 );
+  printf("prim_vgl[26][1][7] = %e\n",prim_vgl[26][1][7]);
   assert( fabs(prim_vgl[26][1][7] - (-7.5014974095310560E-004)) < 1.e-14 );
+  printf("prim_vgl[26][2][7] = %e\n",prim_vgl[26][2][7]);
   assert( fabs(prim_vgl[26][2][7] - (-3.8250692897610380E-003)) < 1.e-14 );
+  printf("prim_vgl[26][3][7] = %e\n",prim_vgl[26][3][7]);
   assert( fabs(prim_vgl[26][3][7] - ( 3.4950559194080275E-003)) < 1.e-14 );
+  printf("prim_vgl[26][4][7] = %e\n",prim_vgl[26][4][7]);
   assert( fabs(prim_vgl[26][4][7] - ( 2.0392163767356572E-002)) < 1.e-14 );
 
 }
@@ -3875,7 +3880,7 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f( &
 
   ! Don't compute exponentials when the result will be almost zero.
   ! TODO : Use numerical precision here
-  cutoff = -dlog(1.d-12)
+  cutoff = 27.631021115928547  !-dlog(1.d-12)
 
   do ipoint = 1, point_num
 
@@ -5501,7 +5506,7 @@ integer function qmckl_compute_ao_value_doc_f(context, &
   info = QMCKL_SUCCESS
 
   ! Don't compute polynomials when the radial part is zero.
-  cutoff = -dlog(1.d-12)
+  cutoff = 27.631021115928547  !-dlog(1.d-12)
 
   do ipoint = 1, point_num
      e_coord(1) = coord(ipoint,1)
@@ -5622,7 +5627,7 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
   }
 
   /* Don't compute polynomials when the radial part is zero. */
-  double cutoff = -log(1.e-12);
+  double cutoff = 27.631021115928547; // -log(1.e-12)
 
 #ifdef HAVE_OPENMP
 #pragma omp parallel
@@ -6247,7 +6252,7 @@ integer function qmckl_compute_ao_vgl_doc_f(context, &
   info = QMCKL_SUCCESS
 
   ! Don't compute polynomials when the radial part is zero.
-  cutoff = -dlog(1.d-12)
+  cutoff = 27.631021115928547  ! -dlog(1.d-12)
 
   do ipoint = 1, point_num
      e_coord(1) = coord(ipoint,1)
@@ -6398,8 +6403,8 @@ qmckl_compute_ao_vgl_hpc_gaussian (
     ao_index[shell_num] = ao_num+1;
   }
 
-  /* Don't compute polynomials when the radial part is zero. */
-  double cutoff = -log(1.e-12);
+  /* Don't compute when the radial part is zero. */
+  double cutoff = 27.631021115928547; // -log(1.e-12)
 
 #ifdef HAVE_OPENMP
 #pragma omp parallel
@@ -6637,11 +6642,6 @@ qmckl_compute_ao_vgl_hpc_gaussian (
         for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) {
 
           const double s1 = ce_mat[ishell-ishell_start][0];
-          if (s1 == 0.0) continue;
-          const double s2 = ce_mat[ishell-ishell_start][1];
-          const double s3 = ce_mat[ishell-ishell_start][2];
-          const double s4 = ce_mat[ishell-ishell_start][3];
-          const double s5 = ce_mat[ishell-ishell_start][4];
 
           const int64_t k = ao_index[ishell];
           double* restrict const ao_vgl_1 = ao_vgl + ipoint*5*ao_num + k;
@@ -6654,6 +6654,22 @@ qmckl_compute_ao_vgl_hpc_gaussian (
           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);
 
+          if (s1 == 0.0) {
+            for (int64_t il=0 ; il<n ; ++il) {
+              ao_vgl_1[il] = 0.0;
+              ao_vgl_2[il] = 0.0;
+              ao_vgl_3[il] = 0.0;
+              ao_vgl_4[il] = 0.0;
+              ao_vgl_5[il] = 0.0;
+            }
+            continue;
+          }
+
+          const double s2 = ce_mat[ishell-ishell_start][1];
+          const double s3 = ce_mat[ishell-ishell_start][2];
+          const double s4 = ce_mat[ishell-ishell_start][3];
+          const double s5 = ce_mat[ishell-ishell_start][4];
+
           double* restrict poly_vgl_1 = NULL;
           double* restrict poly_vgl_2 = NULL;
           double* restrict poly_vgl_3 = NULL;
@@ -6726,7 +6742,7 @@ qmckl_compute_ao_vgl_hpc_gaussian (
               break;
             default:
 #ifdef HAVE_OPENMP
-#pragma omp simd simdlen(8)
+#pragma omp simd
 #endif
               for (int il=0 ; il<n ; ++il) {
                 ao_vgl_1[il] =  poly_vgl_1[il] * s1 * f[il];
@@ -6900,7 +6916,9 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
     qmckl_exit_code rc;
 
     /* Provide required data */
-#ifndef HAVE_HPC
+#ifdef HAVE_HPC
+
+#else
     rc = qmckl_provide_ao_basis_shell_vgl(context);
     if (rc != QMCKL_SUCCESS) {
         return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_shell_vgl", NULL);
@@ -6941,6 +6959,40 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
                                   ctx->ao_basis.expo_per_nucleus,
                                   ctx->ao_basis.coef_per_nucleus,
                                   ctx->ao_basis.ao_vgl);
+
+/* DEBUG
+    rc = qmckl_provide_ao_basis_shell_vgl(context);
+    if (rc != QMCKL_SUCCESS) {
+        return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_shell_vgl", NULL);
+    }
+     int64_t K= ctx->ao_basis.ao_num * 5 * ctx->point.num;
+     double* check = malloc(K*sizeof(double));
+      rc = qmckl_compute_ao_vgl_doc(context,
+                                  ctx->ao_basis.ao_num,
+                                  ctx->ao_basis.shell_num,
+                                  ctx->point.num,
+                                  ctx->nucleus.num,
+                                  ctx->point.coord.data,
+                                  ctx->nucleus.coord.data,
+                                  ctx->ao_basis.nucleus_index,
+                                  ctx->ao_basis.nucleus_shell_num,
+                                  ctx->ao_basis.nucleus_range,
+                                  ctx->ao_basis.nucleus_max_ang_mom,
+                                  ctx->ao_basis.shell_ang_mom,
+                                  ctx->ao_basis.ao_factor,
+                                  ctx->ao_basis.shell_vgl,
+                                  check);
+     for (int64_t i=0 ; i<K ; ++i) {
+       if (fabs(check[i] - ctx->ao_basis.ao_vgl[i]) > 1.e-10) {
+           int a, b, c;
+           a = i/(ctx->ao_basis.ao_num*5);
+           b = (i-a*ctx->ao_basis.ao_num*5)/ctx->ao_basis.ao_num;
+           c = (i-a*ctx->ao_basis.ao_num*5 -b*ctx->ao_basis.ao_num);
+           printf("%d: %d, %d, %d, %e  %e\n", i, a, b, c, check[i], ctx->ao_basis.ao_vgl[i]);
+       }
+     }
+*/
+
       /*
     } else if (ctx->ao_basis.type == 'S') {
       rc = qmck_compute_ao_vgl_hpc_slater(context,