diff --git a/src/ao_basis/aos_value.irp.f b/src/ao_basis/aos_value.irp.f
index dba696a2..1370d4bc 100644
--- a/src/ao_basis/aos_value.irp.f
+++ b/src/ao_basis/aos_value.irp.f
@@ -193,8 +193,8 @@ subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_la
  END_DOC
  double precision, intent(in) :: r(3)
  double precision, intent(out) :: aos_array(ao_num)
- double precision, intent(out) :: aos_grad_array(ao_num,3)
- double precision, intent(out) :: aos_lapl_array(ao_num,3)
+ double precision, intent(out) :: aos_grad_array(3,ao_num)
+ double precision, intent(out) :: aos_lapl_array(3,ao_num)
 
  integer :: power_ao(3)
  integer :: i,j,k,l,m
@@ -215,13 +215,13 @@ subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_la
   do j = 1,Nucl_N_Aos(i)
    k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format
    aos_array(k) = 0.d0
-   aos_grad_array(k,1) = 0.d0
-   aos_grad_array(k,2) = 0.d0
-   aos_grad_array(k,3) = 0.d0
+   aos_grad_array(1,k) = 0.d0
+   aos_grad_array(2,k) = 0.d0
+   aos_grad_array(3,k) = 0.d0
 
-   aos_lapl_array(k,1) = 0.d0
-   aos_lapl_array(k,2) = 0.d0
-   aos_lapl_array(k,3) = 0.d0
+   aos_lapl_array(1,k) = 0.d0
+   aos_lapl_array(2,k) = 0.d0
+   aos_lapl_array(3,k) = 0.d0
 
    power_ao(1:3)= ao_power_ordered_transp_per_nucl(1:3,j,i)
    dx2 = dx**power_ao(1)
@@ -238,7 +238,12 @@ subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_la
    else
     dx3 = 0.d0
    endif
-   dx4 = dble((2 * power_ao(1) + 1))  * dx**(power_ao(1))
+   if(power_ao(1) .ge. 1)then
+    dx4 = dble((2 * power_ao(1) + 1))  * dx**(power_ao(1)) 
+   else
+    dx4 = dble((power_ao(1) + 1))  * dx**(power_ao(1)) 
+   endif
+
    dx5 = dx**(power_ao(1)+2)
 
    if(power_ao(2) .ne. 0)then
@@ -252,7 +257,13 @@ subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_la
    else
     dy3 = 0.d0
    endif
-   dy4 = dble((2 * power_ao(2) + 1))  * dy**(power_ao(2))
+
+   if(power_ao(2) .ge. 1)then
+    dy4 = dble((2 * power_ao(2) + 1))  * dy**(power_ao(2)) 
+   else
+    dy4 = dble((power_ao(2) + 1))  * dy**(power_ao(2)) 
+   endif
+
    dy5 = dy**(power_ao(2)+2)
 
 
@@ -267,7 +278,13 @@ subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_la
    else
     dz3 = 0.d0
    endif
-   dz4 = dble((2 * power_ao(3) + 1))  * dz**(power_ao(3))
+
+   if(power_ao(3) .ge. 1)then
+    dz4 = dble((2 * power_ao(3) + 1))  * dz**(power_ao(3)) 
+   else
+    dz4 = dble((power_ao(3) + 1))  * dz**(power_ao(3)) 
+   endif
+
    dz5 = dz**(power_ao(3)+2)
 
 
@@ -283,13 +300,13 @@ subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_la
    enddo
    aos_array(k) = accu_1 * dx2 * dy2 * dz2
 
-   aos_grad_array(k,1) = accu_1 * dx1  * dy2 * dz2- 2.d0 * dx2 * dx  * dy2 * dz2 * accu_2
-   aos_grad_array(k,2) = accu_1 * dx2  * dy1 * dz2- 2.d0 * dx2 * dy2 * dy  * dz2 * accu_2
-   aos_grad_array(k,3) = accu_1 * dx2  * dy2 * dz1- 2.d0 * dx2 * dy2 * dz2 * dz  * accu_2
+   aos_grad_array(1,k) = accu_1 * dx1  * dy2 * dz2- 2.d0 * dx2 * dx  * dy2 * dz2 * accu_2
+   aos_grad_array(2,k) = accu_1 * dx2  * dy1 * dz2- 2.d0 * dx2 * dy2 * dy  * dz2 * accu_2
+   aos_grad_array(3,k) = accu_1 * dx2  * dy2 * dz1- 2.d0 * dx2 * dy2 * dz2 * dz  * accu_2
 
-   aos_lapl_array(k,1) = accu_1 * dx3  * dy2 * dz2- 2.d0 * dx4 * dy2 * dz2* accu_2 +4.d0 * dx5 *dy2 * dz2* accu_3
-   aos_lapl_array(k,2) = accu_1 * dx2  * dy3 * dz2- 2.d0 * dx2 * dy4 * dz2* accu_2 +4.d0 * dx2 *dy5 * dz2* accu_3
-   aos_lapl_array(k,3) = accu_1 * dx2  * dy2 * dz3- 2.d0 * dx2 * dy2 * dz4* accu_2 +4.d0 * dx2 *dy2 * dz5* accu_3
+   aos_lapl_array(1,k) = accu_1 * dx3  * dy2 * dz2- 2.d0 * dx4 * dy2 * dz2* accu_2 +4.d0 * dx5 *dy2 * dz2* accu_3
+   aos_lapl_array(2,k) = accu_1 * dx2  * dy3 * dz2- 2.d0 * dx2 * dy4 * dz2* accu_2 +4.d0 * dx2 *dy5 * dz2* accu_3
+   aos_lapl_array(3,k) = accu_1 * dx2  * dy2 * dz3- 2.d0 * dx2 * dy2 * dz4* accu_2 +4.d0 * dx2 *dy2 * dz5* accu_3
 
   enddo
  enddo