diff --git a/src/qmckl_ao.org b/src/qmckl_ao.org
index 4371cd9..51736d9 100644
--- a/src/qmckl_ao.org
+++ b/src/qmckl_ao.org
@@ -290,23 +290,46 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL,
   do i=1,3
      Y(i) = X(i) - R(i)
   end do
-  pows(-2:-1,1:3) = 0.d0
-  pows(0,1:3) = 1.d0
+
   lmax_array(1:3) = lmax
-  info = qmckl_ao_powers_f(context, 1_8, Y(1), (/lmax/), pows(1,1), size(pows,1,kind=8)) 
-  if (info /= 0) return
-  info = qmckl_ao_powers_f(context, 1_8, Y(2), (/lmax/), pows(1,2), size(pows,1,kind=8)) 
-  if (info /= 0) return
-  info = qmckl_ao_powers_f(context, 1_8, Y(3), (/lmax/), pows(1,3), size(pows,1,kind=8)) 
-  if (info /= 0) return
 
+  if (lmax == 0) then
+    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(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
 
-  VGL(1,1) = 1.d0
-  vgL(2:5,1) = 0.d0
-  l(1:3,1) = 0
-  n=1
-  dd = 1.d0
-  do d=1,lmax
+    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,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
+
+    n=4
+  endif
+        
+  ! l>=2
+  dd = 2.d0
+  do d=2,lmax
      da = dd
      do a=d,0,-1
         db = dd-da
@@ -314,6 +337,7 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL,
            c  = d  - a  - b
            dc = dd - da - db
            n = n+1
+
            l(1,n) = a
            l(2,n) = b
            l(3,n) = c