diff --git a/configure.ac b/configure.ac
index eb31b93..fcc1dc5 100644
--- a/configure.ac
+++ b/configure.ac
@@ -222,13 +222,13 @@ 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 -Warray-bounds -Wdisabled-optimization \
+-g -Wall -W -Wbad-function-cast -Wcast-qual -Warray-bounds -Wdisabled-optimization \
 -Wpointer-arith -Wcast-align -Wpedantic -Wextra -Walloc-zero -Werror \
 "
         fi
         if test "$GFC" = "yes"; then
 	        FCFLAGS="$FCFLAGS \
--fcheck=all -Waliasing -Wampersand -Wconversion \
+-g -fcheck=all -Waliasing -Wampersand -Wconversion \
 -Wsurprising -ffpe-trap=zero,overflow,underflow \
 -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation \
 -Wreal-q-constant -Wuninitialized  -fbacktrace -finit-real=nan"
diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org
index a982379..55057bd 100644
--- a/org/qmckl_ao.org
+++ b/org/qmckl_ao.org
@@ -3154,7 +3154,7 @@ assert(0 == test_qmckl_ao_gaussian_vgl(context));
          const double* coord,
          const double* nucl_coord,
          const double* expo,
-         double* const primitive_vgl ); 
+         double* const primitive_vgl );
    #+end_src
 
 
@@ -3493,7 +3493,7 @@ for (j=0 ; j<point_num ; ++j) {
          const double* nucl_coord,
          const double* expo,
          const double* coef_normalized,
-         double* const shell_vgl ); 
+         double* const shell_vgl );
    #+end_src
 
    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
@@ -4149,7 +4149,7 @@ assert(0 == test_qmckl_ao_power(context));
           int32_t* const L,
           const int64_t ldl,
           double* const VGL,
-          const int64_t ldv ); 
+          const int64_t ldv );
     #+end_src
 
     #+begin_src c :tangle (eval c) :comments org
@@ -4169,7 +4169,7 @@ qmckl_ao_polynomial_vgl (const qmckl_context context,
   return qmckl_ao_polynomial_vgl_doc (context, X, R, lmax, n, L, ldl, VGL, ldv);
 #else
   return qmckl_ao_polynomial_vgl_doc (context, X, R, lmax, n, L, ldl, VGL, ldv);
-#endif  
+#endif
 }
     #+end_src
 
@@ -4394,7 +4394,7 @@ end function qmckl_ao_polynomial_vgl_doc_f
           int32_t* const L,
           const int64_t ldl,
           double* const VGL,
-          const int64_t ldv ); 
+          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_doc")
@@ -4410,7 +4410,7 @@ end function qmckl_ao_polynomial_vgl_doc_f
           int32_t* const L,
           const int64_t ldl,
           double* const VGL,
-          const int64_t ldv ); 
+          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")
@@ -4426,7 +4426,7 @@ end function qmckl_ao_polynomial_vgl_doc_f
           int32_t* const L,
           const int64_t ldl,
           double* const VGL,
-          const int64_t ldv ); 
+          const int64_t ldv );
     #+end_src
 
     #+begin_src c :tangle (eval c) :comments org
@@ -4445,7 +4445,7 @@ qmckl_ao_polynomial_transp_vgl (const qmckl_context context,
   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  
+#endif
 }
     #+end_src
 
@@ -4523,7 +4523,7 @@ integer function qmckl_ao_polynomial_transp_vgl_doc_f (context,    &
      VGL(4,4) = 1.d0
 
      n=4
-  else 
+  else
      VGL(1,1) = 1.d0
      VGL(1,2:5) = 0.d0
      l(1:3,1) = 0
@@ -4605,7 +4605,7 @@ qmckl_ao_polynomial_transp_vgl_hpc (const qmckl_context context,
       VGL[0         ] = 1.0;
       VGL[ldv       ] = 0.0;
       VGL[ldv<<1    ] = 0.0;
-      VGL[(ldv<<2)-1] = 0.0;
+      VGL[(ldv<<1)+ldv] = 0.0;
       VGL[ldv<<2    ] = 0.0;
       m=1;
       break;
@@ -4619,10 +4619,10 @@ qmckl_ao_polynomial_transp_vgl_hpc (const qmckl_context context,
         const int32_t ltmp[12] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
         for (int i=0 ; i<12 ; ++i)
           L[i] = ltmp[i];
-        
+
       } else {
 
-        int32_t* const l[4] = {L, L+ldl, L+(ldl<<1), L+3*ldl};
+        int32_t* const l[4] = {L, L+ldl, L+(ldl<<1), L+ldl+(ldl<<1)};
         l[0][0] = 0; l[0][1] = 0; l[0][2] = 0;
         l[1][0] = 1; l[1][1] = 0; l[1][2] = 0;
         l[2][0] = 0; l[2][1] = 1; l[2][2] = 0;
@@ -4646,7 +4646,7 @@ qmckl_ao_polynomial_transp_vgl_hpc (const qmckl_context context,
         double* const vgl1 = VGL;
         double* const vgl2 = VGL + ldv;
         double* const vgl3 = VGL + (ldv << 1);
-        double* const vgl4 = VGL + 3*ldv;
+        double* const vgl4 = VGL + ldv + (ldv << 1);
         double* const vgl5 = VGL + (ldv << 2);
 
         for (int32_t k=0 ; k<4 ; ++k) {
@@ -4670,7 +4670,7 @@ qmckl_ao_polynomial_transp_vgl_hpc (const qmckl_context context,
     {
       if (ldv < 10) return QMCKL_INVALID_ARG_9;
       if (ldl == 3) {
-        const int32_t ltmp[30] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 
+        const int32_t ltmp[30] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1,
                                   2, 0, 0, 1, 1, 0, 1, 0, 1, 0, 2, 0,
                                   0, 1, 1, 0, 0, 2};
         for (int i=0 ; i<30 ; ++i)
@@ -4681,11 +4681,11 @@ qmckl_ao_polynomial_transp_vgl_hpc (const qmckl_context context,
         for (int32_t i=0 ; i<10 ; ++i) {
           l[i] = L + i*ldl;
         }
-      
+
         l[0][0] = 0; l[0][1] = 0; l[0][2] = 0;
         l[1][0] = 1; l[1][1] = 0; l[1][2] = 0;
         l[2][0] = 0; l[2][1] = 1; l[2][2] = 0;
-        l[3][0] = 0; l[3][1] = 0; l[3][2] = 1; 
+        l[3][0] = 0; l[3][1] = 0; l[3][2] = 1;
         l[4][0] = 2; l[4][1] = 0; l[4][2] = 0;
         l[5][0] = 1; l[5][1] = 1; l[5][2] = 0;
         l[6][0] = 1; l[6][1] = 0; l[6][2] = 1;
@@ -4694,8 +4694,8 @@ qmckl_ao_polynomial_transp_vgl_hpc (const qmckl_context context,
         l[9][0] = 0; l[9][1] = 0; l[9][2] = 2;
       }
 
-      const double Y[3] = { X[0]-R[0], 
-                            X[1]-R[1], 
+      const double Y[3] = { X[0]-R[0],
+                            X[1]-R[1],
                             X[2]-R[2] };
 
       if (ldv == 50) {
@@ -4743,50 +4743,45 @@ qmckl_ao_polynomial_transp_vgl_hpc (const qmckl_context context,
         vgl5[3] = 0.0; vgl5[4] = 2.0; vgl5[5] = 0.0;
         vgl5[6] = 0.0; vgl5[7] = 2.0; vgl5[8] = 0.0;
         vgl5[9] = 2.0;
-      }    
+      }
       m=10;
       break;
     }
-    /*
-      case 3: 
-      {
-      if (ldv < 20) return QMCKL_INVALID_ARG_9;
-      }
-    ,*/
-  default: 
+  default:
     {
       const int32_t size_max = (lmax+1)*(lmax+2)*(lmax+3)/6;
       if (ldv < size_max) return QMCKL_INVALID_ARG_9;
 
       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);
-      
-      const double Y[3] = { X[0]-R[0], 
-                            X[1]-R[1], 
+      double* const vgl3 = VGL + (ldv<<1);
+      double* const vgl4 = VGL + ldv + (ldv<<1);
+      double* const vgl5 = VGL + (ldv<<2);
+
+      const double Y[3] = { X[0]-R[0],
+                            X[1]-R[1],
                             X[2]-R[2] };
 
+      assert(size_max > lmax+3);
       double pows[3][size_max];
-      
-      for (int32_t i=0 ; i<=2 ; ++i) {
+
+      for (int32_t i=0 ; i<3 ; ++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];
+
+      int32_t* l[size_max];
       for (int32_t i=0 ; i<size_max ; ++i) {
-        l[i] = L + i*ldl;
+        l[i] = &(L[i*ldl]);
       }
-      
+
       for (int32_t i=0 ; i<4 ; ++i) {
         l[i][0] = 0;
         l[i][1] = 0;
@@ -4795,7 +4790,7 @@ qmckl_ao_polynomial_transp_vgl_hpc (const qmckl_context context,
       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;
@@ -4811,41 +4806,41 @@ qmckl_ao_polynomial_transp_vgl_hpc (const qmckl_context context,
       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 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;
           }
@@ -5040,13 +5035,13 @@ 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 };
-int32_t ldv[4] = {1, 4, 10, 20};
-for (int32_t ldl=3 ; ldl<5 ; ++ldl) {
+int32_t ldv[8] = {1, 4, 10, 20, 35, 56, 84, 120};
+for (int32_t ldl=3 ; ldl<=5 ; ++ldl) {
     int64_t n;
-    int32_t L0[24][ldl];
-    int32_t L1[24][ldl];
+    int32_t L0[200][ldl];
+    int32_t L1[200][ldl];
     printf("ldl=%d\n", ldl);
-    for (int32_t lmax=0 ; lmax<=3 ; lmax++) {
+    for (int32_t lmax=0 ; lmax<=7 ; lmax++) {
       double VGL0[5][ldv[lmax]];
       double VGL1[5][ldv[lmax]];
       memset(&L0[0][0], 0, sizeof(L0));
@@ -5064,7 +5059,7 @@ for (int32_t ldl=3 ; ldl<5 ; ++ldl) {
           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]);
@@ -5167,9 +5162,9 @@ integer function qmckl_compute_ao_vgl_doc_f(context, &
      e_coord(2) = coord(ipoint,2)
      e_coord(3) = coord(ipoint,3)
      do inucl=1,nucl_num
-        n_coord(1) = nucl_coord(inucl,1) 
-        n_coord(2) = nucl_coord(inucl,2) 
-        n_coord(3) = nucl_coord(inucl,3) 
+        n_coord(1) = nucl_coord(inucl,1)
+        n_coord(2) = nucl_coord(inucl,2)
+        n_coord(3) = nucl_coord(inucl,3)
 
         ! Test if the point is in the range of the nucleus
         x = e_coord(1) - n_coord(1)
@@ -5182,7 +5177,7 @@ integer function qmckl_compute_ao_vgl_doc_f(context, &
            cycle
         end if
 
-        ! Compute polynomials 
+        ! Compute polynomials
         info = qmckl_ao_polynomial_vgl_doc_f(context, e_coord, n_coord, &
              nucleus_max_ang_mom(inucl), n_poly, powers, 3_8, &
              poly_vgl, 5_8)
@@ -5237,7 +5232,7 @@ end function qmckl_compute_ao_vgl_doc_f
     #+end_src
 
 ** HPC version
-    #+NAME: qmckl_ao_vgl_args_hpc
+    #+NAME: qmckl_ao_vgl_args_hpc_gaussian
    | Variable              | Type                           | In/Out | Description                                  |
    |-----------------------+--------------------------------+--------+----------------------------------------------|
    | ~context~             | ~qmckl_context~                | in     | Global state                                 |
@@ -5262,7 +5257,7 @@ end function qmckl_compute_ao_vgl_doc_f
 
     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
 qmckl_exit_code
-qmckl_compute_ao_vgl_hpc (
+qmckl_compute_ao_vgl_hpc_gaussian (
           const qmckl_context context,
           const int64_t ao_num,
           const int64_t shell_num,
@@ -5283,14 +5278,14 @@ qmckl_compute_ao_vgl_hpc (
           const double* coef_normalized,
           double* const ao_vgl )
 {
-  int32_t lstart[24];
-  for (int32_t l=0 ; l<24 ; ++l) {
+  int32_t lstart[32];
+  for (int32_t l=0 ; l<32 ; ++l) {
     lstart[l] = l*(l+1)*(l+2)/6;
   }
 
   int64_t ao_index[shell_num+1];
   int64_t size_max = 0;
-  { 
+  {
      int64_t k=0;
      for (int inucl=0 ; inucl < nucl_num ; ++inucl) {
         const int64_t ishell_start = nucleus_index[inucl];
@@ -5309,7 +5304,7 @@ qmckl_compute_ao_vgl_hpc (
   double cutoff = -log(1.e-12);
 
 #ifdef HAVE_OPENMP
-  #pragma omp parallel 
+  #pragma omp parallel
 #endif
   {
     qmckl_exit_code rc;
@@ -5321,10 +5316,10 @@ qmckl_compute_ao_vgl_hpc (
                                 {0.0, 1.0, 0.0, 0.0},
                                 {0.0, 0.0, 1.0, 0.0},
                                 {0.0, 0.0, 0.0, 1.0}};
-    double  poly_vgl_l2[5][10] = {{1., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, 
-                                  {0., 1., 0., 0., 0., 0., 0., 0., 0., 0.}, 
+    double  poly_vgl_l2[5][10] = {{1., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
+                                  {0., 1., 0., 0., 0., 0., 0., 0., 0., 0.},
                                   {0., 0., 1., 0., 0., 0., 0., 0., 0., 0.},
-                                  {0., 0., 0., 1., 0., 0., 0., 0., 0., 0.}, 
+                                  {0., 0., 0., 1., 0., 0., 0., 0., 0., 0.},
                                   {0., 0., 0., 0., 2., 0., 0., 2., 0., 2.}};
     double  poly_vgl[5][size_max];
 
@@ -5332,12 +5327,12 @@ qmckl_compute_ao_vgl_hpc (
     #pragma omp for
 #endif
     for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
-       const double e_coord[3] = { coord[ipoint], 
+       const double e_coord[3] = { coord[ipoint],
                                    coord[ipoint + point_num],
                                    coord[ipoint + 2*point_num] };
 
        for (int64_t inucl=0 ; inucl < nucl_num ; ++inucl) {
-          const double n_coord[3] = { nucl_coord[inucl], 
+          const double n_coord[3] = { nucl_coord[inucl],
                                       nucl_coord[inucl + nucl_num],
                                       nucl_coord[inucl + 2*nucl_num] };
 
@@ -5357,13 +5352,13 @@ qmckl_compute_ao_vgl_hpc (
           case 0:
             break;
 
-          case 1: 
+          case 1:
             poly_vgl_l1[0][1] = x;
             poly_vgl_l1[0][2] = y;
             poly_vgl_l1[0][3] = z;
             break;
 
-          case 2: 
+          case 2:
             poly_vgl_l2[0][1] = x;
             poly_vgl_l2[0][2] = y;
             poly_vgl_l2[0][3] = z;
@@ -5385,20 +5380,19 @@ qmckl_compute_ao_vgl_hpc (
             break;
 
           default:
-              rc = qmckl_ao_polynomial_transp_vgl_hpc(context, e_coord, n_coord, 
+              rc = qmckl_ao_polynomial_transp_vgl_hpc(context, e_coord, n_coord,
                                                    nucleus_max_ang_mom[inucl],
-                                                   &n_poly, powers, (int64_t) 3, 
+                                                   &n_poly, powers, (int64_t) 3,
                                                    &(poly_vgl[0][0]), size_max);
-              
+
               assert (rc == QMCKL_SUCCESS);
               break;
           }
- 
+
           const int64_t ishell_start = nucleus_index[inucl];
           const int64_t ishell_end   = nucleus_index[inucl] + nucleus_shell_num[inucl];
           for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) {
 
-             /* /!\  Gaussian fuctions */
              int64_t nidx = 0;
              const int64_t iprim_start = shell_prim_index[ishell];
              const int64_t iprim_end   = shell_prim_index[ishell] + shell_prim_num[ishell];
@@ -5429,9 +5423,10 @@ qmckl_compute_ao_vgl_hpc (
              const double s4 = s6_*z;
              const double s5 = s5_;
 
+
              const int32_t l = shell_ang_mom[ishell];
              const int32_t n = lstart[l+1]-lstart[l];
-             const int64_t k = ao_index[ishell]; 
+             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]);
@@ -5467,7 +5462,7 @@ qmckl_compute_ao_vgl_hpc (
                      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 + 
+                                     2.0*(poly_vgl_2[il] * s2 +
                                           poly_vgl_3[il] * s3 +
                                           poly_vgl_4[il] * s4 )) * f[il];
                    }
@@ -5484,7 +5479,7 @@ qmckl_compute_ao_vgl_hpc (
                      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 + 
+                                     2.0*(poly_vgl_2[il] * s2 +
                                           poly_vgl_3[il] * s3 +
                                           poly_vgl_4[il] * s4 )) * f[il];
                    }
@@ -5501,7 +5496,7 @@ qmckl_compute_ao_vgl_hpc (
                      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 + 
+                                     2.0*(poly_vgl_2[il] * s2 +
                                           poly_vgl_3[il] * s3 +
                                           poly_vgl_4[il] * s4 )) * f[il];
                    }
@@ -5550,7 +5545,7 @@ qmckl_compute_ao_vgl_hpc (
           double* const ao_vgl );
     #+end_src
     #+begin_src c :tangle (eval h_private_func) :comments org
-    qmckl_exit_code qmckl_compute_ao_vgl_hpc (
+    qmckl_exit_code qmckl_compute_ao_vgl_hpc_gaussian (
           const qmckl_context context,
           const int64_t ao_num,
           const int64_t shell_num,
@@ -5633,7 +5628,7 @@ qmckl_compute_ao_vgl_hpc (
 
     end function qmckl_compute_ao_vgl_doc
     #+end_src
-    
+
 
 *** Provide                                                        :noexport:
 
@@ -5691,7 +5686,8 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
       ctx->ao_basis.ao_vgl = ao_vgl;
     }
 #ifdef HAVE_HPC
-    rc = qmckl_compute_ao_vgl_hpc(context,
+    if (ctx->ao_basis.type == 'G') {
+      rc = qmckl_compute_ao_vgl_hpc_gaussian(context,
                                   ctx->ao_basis.ao_num,
                                   ctx->ao_basis.shell_num,
                                   ctx->ao_basis.prim_num,
@@ -5710,6 +5706,45 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
                                   ctx->ao_basis.exponent,
                                   ctx->ao_basis.coefficient_normalized,
                                   ctx->ao_basis.ao_vgl);
+      /*
+    } else if (ctx->ao_basis.type == 'S') {
+      rc = qmck_compute_ao_vgl_hpc_slater(context,
+                                  ctx->ao_basis.ao_num,
+                                  ctx->ao_basis.shell_num,
+                                  ctx->ao_basis.prim_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.shell_prim_index,
+                                  ctx->ao_basis.shell_prim_num,
+                                  ctx->ao_basis.ao_factor,
+                                  ctx->ao_basis.exponent,
+                                  ctx->ao_basis.coefficient_normalized,
+                                  ctx->ao_basis.ao_vgl);
+      ,*/
+    } else {
+      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,
+                                  ctx->ao_basis.ao_vgl);
+    }
 #else
     rc = qmckl_compute_ao_vgl_doc(context,
                                   ctx->ao_basis.ao_num,