diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org
index 4f419a6..e908769 100644
--- a/org/qmckl_ao.org
+++ b/org/qmckl_ao.org
@@ -3349,15 +3349,15 @@ qmckl_ao_gaussian_vgl(const qmckl_context context,
 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)
-  integer*8 , intent(in)  :: n
-  real*8    , intent(in)  :: A(n)
-  real*8    , intent(out) :: VGL(ldv,5)
-  integer*8 , intent(in)  :: ldv
+  integer*8        , intent(in)  :: context
+  double precision , intent(in)  :: X(3), R(3)
+  integer*8        , intent(in)  :: n
+  double precision , intent(in)  :: A(n)
+  double precision , intent(out) :: VGL(ldv,5)
+  integer*8        , intent(in)  :: ldv
 
   integer*8         :: i,j
-  real*8            :: Y(3), r2, t, u, v
+  double precision  :: Y(3), r2, t, u, v
 
   info = QMCKL_SUCCESS
 
@@ -4626,19 +4626,19 @@ integer function qmckl_ao_polynomial_vgl_doc_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)
-  integer   , intent(in)  :: lmax
-  integer*8 , intent(out) :: n
-  integer   , intent(out) :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6)
-  integer*8 , intent(in)  :: ldl
-  real*8    , intent(out) :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6)
-  integer*8 , intent(in)  :: ldv
+  integer*8,        intent(in)  :: context
+  double precision, intent(in)  :: X(3), R(3)
+  integer,          intent(in)  :: lmax
+  integer*8,        intent(out) :: n
+  integer,          intent(out) :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6)
+  integer*8,        intent(in)  :: ldl
+  double precision, intent(out) :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6)
+  integer*8,        intent(in)  :: ldv
 
   integer*8         :: i,j
   integer           :: a,b,c,d
-  real*8            :: Y(3)
-  real*8            :: pows(-2:lmax,3)
+  double precision  :: Y(3)
+  double precision  :: pows(-2:lmax,3)
   double precision  :: xy, yz, xz
   double precision  :: da, db, dc, dd
 
@@ -4671,8 +4671,13 @@ integer function qmckl_ao_polynomial_vgl_doc_f (context,    &
 
   if (lmax == 0) then
      VGL(1,1) = 1.d0
-     VGL(2:5,1) = 0.d0
-     l(1:3,1) = 0
+     VGL(2,1) = 0.d0
+     VGL(3,1) = 0.d0
+     VGL(4,1) = 0.d0
+     VGL(5,1) = 0.d0
+     l(1,1) = 0
+     l(2,1) = 0
+     l(3,1) = 0
      n=1
   else if (lmax > 0) then
      pows(-2:0,1:3) = 1.d0
@@ -4683,23 +4688,19 @@ integer function qmckl_ao_polynomial_vgl_doc_f (context,    &
      end do
 
      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
-
-     l  (1,2) = 1
+     VGL(1,1) = 1.d0
      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  (3,4) = 1
      VGL(1,4) = pows(1,3)
      VGL(4,4) = 1.d0
 
+     l  (1:3,1:4) = 0
+     l  (1,2) = 1
+     l  (2,3) = 1
+     l  (3,4) = 1
+
      n=4
   endif
 
@@ -6422,6 +6423,7 @@ integer function qmckl_compute_ao_vgl_doc_f(context, &
      e_coord(1) = coord(ipoint,1)
      e_coord(2) = coord(ipoint,2)
      e_coord(3) = coord(ipoint,3)
+     ao_vgl(:,:,ipoint) = 0.d0
      do inucl=1,nucl_num
         n_coord(1) = nucl_coord(inucl,1)
         n_coord(2) = nucl_coord(inucl,2)