diff --git a/plugins/local/non_h_ints_mu/NEED b/plugins/local/non_h_ints_mu/NEED
index 48c1c24b..5ca1d543 100644
--- a/plugins/local/non_h_ints_mu/NEED
+++ b/plugins/local/non_h_ints_mu/NEED
@@ -3,3 +3,4 @@ hamiltonian
 jastrow
 ao_tc_eff_map
 bi_ortho_mos
+trexio
diff --git a/plugins/local/non_h_ints_mu/deb_aos.irp.f b/plugins/local/non_h_ints_mu/deb_aos.irp.f
index c9bc9c9a..4012f47c 100644
--- a/plugins/local/non_h_ints_mu/deb_aos.irp.f
+++ b/plugins/local/non_h_ints_mu/deb_aos.irp.f
@@ -31,24 +31,63 @@ subroutine print_aos()
   integer          :: i, ipoint
   double precision :: r(3)
   double precision :: ao_val, ao_der(3), ao_lap
+  double precision :: mo_val, mo_der(3), mo_lap
 
   PROVIDE final_grid_points aos_in_r_array aos_grad_in_r_array aos_lapl_in_r_array
 
   do ipoint = 1, n_points_final_grid
     r(:) = final_grid_points(:,ipoint)
-    print*, r
+    write(1000, '(3(f15.7, 3X))') r
   enddo
+  
+double precision :: accu_vgl(5)
+double precision :: accu_vgl_nrm(5)
 
   do ipoint = 1, n_points_final_grid
-    r(:) = final_grid_points(:,ipoint)
     do i = 1, ao_num
       ao_val    = aos_in_r_array     (i,ipoint)
       ao_der(:) = aos_grad_in_r_array(i,ipoint,:)
       ao_lap    = aos_lapl_in_r_array(1,i,ipoint) + aos_lapl_in_r_array(2,i,ipoint) + aos_lapl_in_r_array(3,i,ipoint)
-      write(*, '(5(f15.7, 3X))') ao_val, ao_der, ao_lap
+      write(111, '(5(f15.7, 3X))') ao_val, ao_der, ao_lap
     enddo
   enddo
  
+  do ipoint = 1, n_points_final_grid
+    do i = 1, ao_num
+      ao_val    = aos_in_r_array_qmckl     (i,ipoint)
+      ao_der(:) = aos_grad_in_r_array_qmckl(i,ipoint,:)
+      ao_lap    = aos_lapl_in_r_array_qmckl(i,ipoint)
+      write(222, '(5(f15.7, 3X))') ao_val, ao_der, ao_lap
+    enddo
+  enddo
+ 
+  accu_vgl = 0.d0
+  accu_vgl_nrm = 0.d0
+  do ipoint = 1, n_points_final_grid
+    do i = 1, ao_num
+      ao_val    = aos_in_r_array     (i,ipoint)
+      ao_der(:) = aos_grad_in_r_array(i,ipoint,:)
+      ao_lap    = aos_lapl_in_r_array(1,i,ipoint) + aos_lapl_in_r_array(2,i,ipoint) + aos_lapl_in_r_array(3,i,ipoint)
+      accu_vgl_nrm(1) += dabs(ao_val)
+      accu_vgl_nrm(2) += dabs(ao_der(1))
+      accu_vgl_nrm(3) += dabs(ao_der(2))
+      accu_vgl_nrm(4) += dabs(ao_der(3))
+      accu_vgl_nrm(5) += dabs(ao_lap)
+
+      ao_val    -= aos_in_r_array_qmckl     (i,ipoint)
+      ao_der(:) -= aos_grad_in_r_array_qmckl(i,ipoint,:)
+      ao_lap    -= aos_lapl_in_r_array_qmckl(i,ipoint)
+      accu_vgl(1) += dabs(ao_val)
+      accu_vgl(2) += dabs(ao_der(1))
+      accu_vgl(3) += dabs(ao_der(2))
+      accu_vgl(4) += dabs(ao_der(3))
+      accu_vgl(5) += dabs(ao_lap)
+    enddo
+
+  enddo
+  accu_vgl(:) *= 1.d0 / accu_vgl_nrm(:)
+  print *, accu_vgl
+
   return
 end
 
diff --git a/plugins/local/non_h_ints_mu/qmckl.irp.f b/plugins/local/non_h_ints_mu/qmckl.irp.f
index 1df80457..4d419e24 100644
--- a/plugins/local/non_h_ints_mu/qmckl.irp.f
+++ b/plugins/local/non_h_ints_mu/qmckl.irp.f
@@ -75,3 +75,107 @@ BEGIN_PROVIDER [ integer*8, qmckl_ctx_jastrow ]
   endif
 
 END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, aos_in_r_array_qmckl, (ao_num,n_points_final_grid)]
+&BEGIN_PROVIDER [ double precision, aos_grad_in_r_array_qmckl, (ao_num,n_points_final_grid,3)]
+&BEGIN_PROVIDER [ double precision, aos_lapl_in_r_array_qmckl, (ao_num, n_points_final_grid)]
+ implicit none
+ BEGIN_DOC
+ ! AOS computed with qmckl
+ END_DOC
+ use qmckl
+
+ integer*8 :: qmckl_ctx
+ integer(qmckl_exit_code) :: rc
+
+ qmckl_ctx = qmckl_context_create()
+
+ rc = qmckl_trexio_read(qmckl_ctx, trexio_file, 1_8*len(trim(trexio_filename)))
+ if (rc /= QMCKL_SUCCESS) then
+   print *, irp_here, 'qmckl error in read_trexio'
+   rc = qmckl_check(qmckl_ctx, rc)
+   stop -1
+ endif
+
+ rc = qmckl_set_point(qmckl_ctx, 'N', n_points_final_grid*1_8, final_grid_points, n_points_final_grid*3_8)
+ if (rc /= QMCKL_SUCCESS) then
+   print *, irp_here, 'qmckl error in set_electron_point'
+   rc = qmckl_check(qmckl_ctx, rc)
+   stop -1
+ endif
+
+ double precision, allocatable :: vgl(:,:,:)
+ allocate( vgl(ao_num,5,n_points_final_grid))
+ rc = qmckl_get_ao_basis_ao_vgl_inplace(qmckl_ctx, vgl, n_points_final_grid*ao_num*5_8)
+ if (rc /= QMCKL_SUCCESS) then
+   print *, irp_here, 'qmckl error in get_ao_vgl'
+   rc = qmckl_check(qmckl_ctx, rc)
+   stop -1
+ endif
+
+ integer :: i,k
+ do k=1,n_points_final_grid
+   do i=1,ao_num
+     aos_in_r_array_qmckl(i,k) = vgl(i,1,k)
+     aos_grad_in_r_array_qmckl(i,k,1) = vgl(i,2,k)
+     aos_grad_in_r_array_qmckl(i,k,2) = vgl(i,3,k)
+     aos_grad_in_r_array_qmckl(i,k,3) = vgl(i,4,k)
+     aos_lapl_in_r_array_qmckl(i,k) = vgl(i,5,k)
+   enddo
+ enddo
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, mos_in_r_array_qmckl, (mo_num,n_points_final_grid)]
+&BEGIN_PROVIDER [ double precision, mos_grad_in_r_array_qmckl, (mo_num,n_points_final_grid,3)]
+&BEGIN_PROVIDER [ double precision, mos_lapl_in_r_array_qmckl, (mo_num, n_points_final_grid)]
+ implicit none
+ BEGIN_DOC
+ ! moS computed with qmckl
+ END_DOC
+ use qmckl
+
+ integer*8 :: qmckl_ctx
+ integer(qmckl_exit_code) :: rc
+
+ qmckl_ctx = qmckl_context_create()
+
+ rc = qmckl_trexio_read(qmckl_ctx, trexio_file, 1_8*len(trim(trexio_filename)))
+ if (rc /= QMCKL_SUCCESS) then
+   print *, irp_here, 'qmckl error in read_trexio'
+   rc = qmckl_check(qmckl_ctx, rc)
+   stop -1
+ endif
+
+ rc = qmckl_set_point(qmckl_ctx, 'N', n_points_final_grid*1_8, final_grid_points, n_points_final_grid*3_8)
+ if (rc /= QMCKL_SUCCESS) then
+   print *, irp_here, 'qmckl error in set_electron_point'
+   rc = qmckl_check(qmckl_ctx, rc)
+   stop -1
+ endif
+
+ double precision, allocatable :: vgl(:,:,:)
+ allocate( vgl(mo_num,5,n_points_final_grid))
+ rc = qmckl_get_mo_basis_mo_vgl_inplace(qmckl_ctx, vgl, n_points_final_grid*mo_num*5_8)
+ if (rc /= QMCKL_SUCCESS) then
+   print *, irp_here, 'qmckl error in get_mo_vgl'
+   rc = qmckl_check(qmckl_ctx, rc)
+   stop -1
+ endif
+
+ integer :: i,k
+ do k=1,n_points_final_grid
+   do i=1,mo_num
+     mos_in_r_array_qmckl(i,k) = vgl(i,1,k)
+     mos_grad_in_r_array_qmckl(i,k,1) = vgl(i,2,k)
+     mos_grad_in_r_array_qmckl(i,k,2) = vgl(i,3,k)
+     mos_grad_in_r_array_qmckl(i,k,3) = vgl(i,4,k)
+     mos_lapl_in_r_array_qmckl(i,k) = vgl(i,5,k)
+   enddo
+ enddo
+
+END_PROVIDER
+
+
diff --git a/src/ao_basis/aos_in_r.irp.f b/src/ao_basis/aos_in_r.irp.f
index 1b1595a3..053c86a2 100644
--- a/src/ao_basis/aos_in_r.irp.f
+++ b/src/ao_basis/aos_in_r.irp.f
@@ -1,67 +1,76 @@
-double precision function ao_value(i,r)
- implicit none
- BEGIN_DOC
-! Returns the value of the i-th ao at point $\textbf{r}$
- END_DOC
- double precision, intent(in) :: r(3)
- integer, intent(in) :: i
 
- integer :: m,num_ao
- double precision :: center_ao(3)
- double precision :: beta
- integer :: power_ao(3)
- double precision :: accu,dx,dy,dz,r2
- num_ao = ao_nucl(i)
- power_ao(1:3)= ao_power(i,1:3)
- center_ao(1:3) = nucl_coord(num_ao,1:3)
- dx = (r(1) - center_ao(1))
- dy = (r(2) - center_ao(2))
- dz = (r(3) - center_ao(3))
- r2 = dx*dx + dy*dy + dz*dz
- dx = dx**power_ao(1)
- dy = dy**power_ao(2)
- dz = dz**power_ao(3)
+! ---
 
- accu = 0.d0
- do m=1,ao_prim_num(i)
-   beta = ao_expo_ordered_transp(m,i)
-   accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2)
- enddo
- ao_value = accu * dx * dy * dz
+double precision function ao_value(i, r)
+
+  BEGIN_DOC
+  ! Returns the value of the i-th ao at point $\textbf{r}$
+  END_DOC
+
+  implicit none
+  integer,          intent(in) :: i
+  double precision, intent(in) :: r(3)
+
+  integer          :: m, num_ao
+  integer          :: power_ao(3)
+  double precision :: center_ao(3)
+  double precision :: beta
+  double precision :: accu, dx, dy, dz, r2
+
+  num_ao = ao_nucl(i)
+  power_ao(1:3) = ao_power(i,1:3)
+  center_ao(1:3) = nucl_coord(num_ao,1:3)
+  dx = r(1) - center_ao(1)
+  dy = r(2) - center_ao(2)
+  dz = r(3) - center_ao(3)
+  r2 = dx*dx + dy*dy + dz*dz
+  dx = dx**power_ao(1)
+  dy = dy**power_ao(2)
+  dz = dz**power_ao(3)
+ 
+  accu = 0.d0
+  do m = 1, ao_prim_num(i)
+    beta = ao_expo_ordered_transp(m,i)
+    accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2)
+  enddo
+  ao_value = accu * dx * dy * dz
 
 end
 
 
-double precision function primitive_value(i,j,r)
- implicit none
- BEGIN_DOC
-! Returns the value of the j-th primitive of the i-th |AO| at point $\textbf{r}
-! **without the coefficient**
- END_DOC
- double precision, intent(in) :: r(3)
- integer, intent(in) :: i,j
+double precision function primitive_value(i, j, r)
 
- integer :: m,num_ao
- double precision :: center_ao(3)
- double precision :: beta
- integer :: power_ao(3)
- double precision :: accu,dx,dy,dz,r2
- num_ao = ao_nucl(i)
- power_ao(1:3)= ao_power(i,1:3)
- center_ao(1:3) = nucl_coord(num_ao,1:3)
- dx = (r(1) - center_ao(1))
- dy = (r(2) - center_ao(2))
- dz = (r(3) - center_ao(3))
- r2 = dx*dx + dy*dy + dz*dz
- dx = dx**power_ao(1)
- dy = dy**power_ao(2)
- dz = dz**power_ao(3)
+  BEGIN_DOC
+  ! Returns the value of the j-th primitive of the i-th |AO| at point $\textbf{r}
+  ! **without the coefficient**
+  END_DOC
 
- accu = 0.d0
- m=j
- beta = ao_expo_ordered_transp(m,i)
- accu += dexp(-beta*r2)
- primitive_value = accu * dx * dy * dz
+  implicit none
+  integer,          intent(in) :: i, j
+  double precision, intent(in) :: r(3)
+ 
+  integer          :: m, num_ao
+  integer          :: power_ao(3)
+  double precision :: center_ao(3)
+  double precision :: beta
+  double precision :: accu, dx, dy, dz, r2
+
+  num_ao = ao_nucl(i)
+  power_ao(1:3)= ao_power(i,1:3)
+  center_ao(1:3) = nucl_coord(num_ao,1:3)
+  dx = r(1) - center_ao(1)
+  dy = r(2) - center_ao(2)
+  dz = r(3) - center_ao(3)
+  r2 = dx*dx + dy*dy + dz*dz
+  dx = dx**power_ao(1)
+  dy = dy**power_ao(2)
+  dz = dz**power_ao(3)
+ 
+  accu = 0.d0
+  m = j
+  beta = ao_expo_ordered_transp(m,i)
+  accu += dexp(-beta*r2)
+  primitive_value = accu * dx * dy * dz
 
 end
 
@@ -104,9 +113,9 @@ subroutine give_all_aos_at_r(r, tmp_array)
       dz2 = dz**p_ao(3)
 
       tmp_array(k) = 0.d0
-      do l = 1,ao_prim_num(k)
+      do l = 1, ao_prim_num(k)
         beta = ao_expo_ordered_transp_per_nucl(l,j,i)
-        if(dabs(beta*r2).gt.40.d0) cycle
+        if(beta*r2.gt.50.d0) cycle
 
         tmp_array(k) += ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
       enddo
@@ -120,207 +129,232 @@ end
 
 ! ---
 
-subroutine give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
- implicit none
- BEGIN_DOC
-! input : r(1) ==> r(1) = x, r(2) = y, r(3) = z
-!
-! output : 
-!
-! * aos_array(i) = ao(i) evaluated at ro
-! * aos_grad_array(1,i) = gradient X of the ao(i) evaluated at $\textbf{r}$
-!
- END_DOC
- double precision, intent(in) :: r(3)
- double precision, intent(out) :: aos_array(ao_num)
- double precision, intent(out) :: aos_grad_array(3,ao_num)
+subroutine give_all_aos_and_grad_at_r(r, aos_array, aos_grad_array)
 
- integer :: power_ao(3)
- integer :: i,j,k,l,m
- double precision :: dx,dy,dz,r2
- double precision ::      dx2,dy2,dz2
- double precision ::      dx1,dy1,dz1
- double precision :: center_ao(3)
- double precision :: beta,accu_1,accu_2,contrib
- do i = 1, nucl_num
-  center_ao(1:3) = nucl_coord(i,1:3)
-  dx = (r(1) - center_ao(1))
-  dy = (r(2) - center_ao(2))
-  dz = (r(3) - center_ao(3))
-  r2 = dx*dx + dy*dy + dz*dz
-  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(1,k) = 0.d0
-   aos_grad_array(2,k) = 0.d0
-   aos_grad_array(3,k) = 0.d0
-   power_ao(1:3)= ao_power_ordered_transp_per_nucl(1:3,j,i)
-   dx2 = dx**power_ao(1)
-   dy2 = dy**power_ao(2)
-   dz2 = dz**power_ao(3)
-   if(power_ao(1) .ne. 0)then
-    dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1)
-   else
-    dx1 = 0.d0
-   endif
-   if(power_ao(2) .ne. 0)then
-    dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1)
-   else
-    dy1 = 0.d0
-   endif
-   if(power_ao(3) .ne. 0)then
-    dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1)
-   else
-    dz1 = 0.d0
-   endif
-   accu_1 = 0.d0
-   accu_2 = 0.d0
-   do l = 1,ao_prim_num(k)
-    beta = ao_expo_ordered_transp_per_nucl(l,j,i)
-    contrib = 0.d0
-    if(beta*r2.gt.50.d0)cycle
-    contrib = ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
-    accu_1 += contrib
-    accu_2 += contrib * beta
-   enddo
-   aos_array(k) = accu_1 * dx2 * dy2 * dz2
-   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
+  BEGIN_DOC
+  !
+  ! input : r(1) ==> r(1) = x, r(2) = y, r(3) = z
+  !
+  ! output : 
+  !
+  ! * aos_array(i) = ao(i) evaluated at ro
+  ! * aos_grad_array(1,i) = gradient X of the ao(i) evaluated at $\textbf{r}$
+  !
+  END_DOC
+
+  implicit none
+  double precision, intent(in)  :: r(3)
+  double precision, intent(out) :: aos_array(ao_num)
+  double precision, intent(out) :: aos_grad_array(3,ao_num)
+
+  integer                       :: power_ao(3)
+  integer                       :: i, j, k, l, m
+  double precision              :: dx, dy, dz, r2
+  double precision              :: dx1, dy1, dz1
+  double precision              :: dx2, dy2, dz2
+  double precision              :: center_ao(3)
+  double precision              :: beta, accu_1, accu_2, contrib
+
+  do i = 1, nucl_num
+
+    center_ao(1:3) = nucl_coord(i,1:3)
+
+    dx = r(1) - center_ao(1)
+    dy = r(2) - center_ao(2)
+    dz = r(3) - center_ao(3)
+    r2 = dx*dx + dy*dy + dz*dz
+
+    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(1,k) = 0.d0
+      aos_grad_array(2,k) = 0.d0
+      aos_grad_array(3,k) = 0.d0
+
+      power_ao(1:3) = ao_power_ordered_transp_per_nucl(1:3,j,i)
+      dx2 = dx**power_ao(1)
+      dy2 = dy**power_ao(2)
+      dz2 = dz**power_ao(3)
+
+      dx1 = 0.d0
+      if(power_ao(1) .ne. 0) then
+        dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1)
+      endif
+
+      dy1 = 0.d0
+      if(power_ao(2) .ne. 0) then
+        dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1)
+      endif
+
+      dz1 = 0.d0
+      if(power_ao(3) .ne. 0) then
+        dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1)
+      endif
+
+      accu_1 = 0.d0
+      accu_2 = 0.d0
+      do l = 1, ao_prim_num(k)
+        beta = ao_expo_ordered_transp_per_nucl(l,j,i)
+        if(beta*r2.gt.50.d0) cycle
+        contrib = ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
+        accu_1 += contrib
+        accu_2 += contrib * beta
+      enddo
+
+      aos_array(k) = accu_1 * dx2 * dy2 * dz2
+      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
+    enddo
   enddo
- enddo
+
 end
 
+! ---
 
-subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
- implicit none
- BEGIN_DOC
-! input  : r(1) ==> r(1) = x, r(2) = y, r(3) = z
-!
-! output :
-!
-! * aos_array(i) = ao(i) evaluated at $\textbf{r}$
-! * aos_grad_array(1,i) = $\nabla_x$ of the ao(i) evaluated at $\textbf{r}$
- END_DOC
- double precision, intent(in) :: r(3)
- double precision, intent(out) :: aos_array(ao_num)
- double precision, intent(out) :: aos_grad_array(3,ao_num)
- double precision, intent(out) :: aos_lapl_array(3,ao_num)
+subroutine give_all_aos_and_grad_and_lapl_at_r(r, aos_array, aos_grad_array, aos_lapl_array)
 
- integer :: power_ao(3)
- integer :: i,j,k,l,m
- double precision :: dx,dy,dz,r2
- double precision ::      dx2,dy2,dz2
- double precision ::      dx1,dy1,dz1
- double precision ::      dx3,dy3,dz3
- double precision ::      dx4,dy4,dz4
- double precision ::      dx5,dy5,dz5
- double precision :: center_ao(3)
- double precision :: beta,accu_1,accu_2,accu_3,contrib
- do i = 1, nucl_num
-  center_ao(1:3) = nucl_coord(i,1:3)
-  dx = (r(1) - center_ao(1))
-  dy = (r(2) - center_ao(2))
-  dz = (r(3) - center_ao(3))
-  r2 = dx*dx + dy*dy + dz*dz
-  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(1,k) = 0.d0
-   aos_grad_array(2,k) = 0.d0
-   aos_grad_array(3,k) = 0.d0
+  BEGIN_DOC
+  !
+  ! input  : r(1) ==> r(1) = x, r(2) = y, r(3) = z
+  !
+  ! output :
+  !
+  ! * aos_array(i) = ao(i) evaluated at $\textbf{r}$
+  ! * aos_grad_array(1,i) = $\nabla_x$ of the ao(i) evaluated at $\textbf{r}$
+  !
+  END_DOC
 
-   aos_lapl_array(1,k) = 0.d0
-   aos_lapl_array(2,k) = 0.d0
-   aos_lapl_array(3,k) = 0.d0
+  implicit none
+  double precision, intent(in)  :: r(3)
+  double precision, intent(out) :: aos_array(ao_num)
+  double precision, intent(out) :: aos_grad_array(3,ao_num)
+  double precision, intent(out) :: aos_lapl_array(3,ao_num)
 
-   power_ao(1:3)= ao_power_ordered_transp_per_nucl(1:3,j,i)
-   dx2 = dx**power_ao(1)
-   dy2 = dy**power_ao(2)
-   dz2 = dz**power_ao(3)
-   if(power_ao(1) .ne. 0)then
-    dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1)
-   else
-    dx1 = 0.d0
-   endif
-   ! For the Laplacian
-   if(power_ao(1) .ge. 2)then
-    dx3 = dble(power_ao(1)) * dble((power_ao(1)-1))  * dx**(power_ao(1)-2)
-   else
-    dx3 = 0.d0
-   endif
-   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
+  integer                       :: power_ao(3)
+  integer                       :: i, j, k, l, m
+  double precision              :: dx, dy, dz, r2
+  double precision              :: dx1, dy1, dz1
+  double precision              :: dx2, dy2, dz2
+  double precision              :: dx3, dy3, dz3
+  double precision              :: dx4, dy4, dz4
+  double precision              :: dx5, dy5, dz5
+  double precision              :: center_ao(3)
+  double precision              :: beta, accu_1, accu_2, accu_3, contrib
 
-   dx5 = dx**(power_ao(1)+2)
+  do i = 1, nucl_num
 
-   if(power_ao(2) .ne. 0)then
-    dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1)
-   else
-    dy1 = 0.d0
-   endif
-   ! For the Laplacian
-   if(power_ao(2) .ge. 2)then
-    dy3 = dble(power_ao(2)) * dble((power_ao(2)-1))  * dy**(power_ao(2)-2)
-   else
-    dy3 = 0.d0
-   endif
+    center_ao(1:3) = nucl_coord(i,1:3)
 
-   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
+    dx = r(1) - center_ao(1)
+    dy = r(2) - center_ao(2)
+    dz = r(3) - center_ao(3)
+    r2 = dx*dx + dy*dy + dz*dz
+    
+    do j = 1, Nucl_N_Aos(i)
 
-   dy5 = dy**(power_ao(2)+2)
+      k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format
 
+      aos_array(k) = 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(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)
+      dy2 = dy**power_ao(2)
+      dz2 = dz**power_ao(3)
 
-   if(power_ao(3) .ne. 0)then
-    dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1)
-   else
-    dz1 = 0.d0
-   endif
-   ! For the Laplacian
-   if(power_ao(3) .ge. 2)then
-    dz3 = dble(power_ao(3)) * dble((power_ao(3)-1))  * dz**(power_ao(3)-2)
-   else
-    dz3 = 0.d0
-   endif
+      ! ---
 
-   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
+      dx1 = 0.d0
+      if(power_ao(1) .ne. 0) then
+        dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1)
+      endif
 
-   dz5 = dz**(power_ao(3)+2)
+      dx3 = 0.d0
+      if(power_ao(1) .ge. 2) then
+        dx3 = dble(power_ao(1)) * dble((power_ao(1)-1)) * dx**(power_ao(1)-2)
+      endif
 
+      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)
+  
+      ! ---
+      
+      dy1 = 0.d0
+      if(power_ao(2) .ne. 0) then
+        dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1)
+      endif
 
-   accu_1 = 0.d0
-   accu_2 = 0.d0
-   accu_3 = 0.d0
-   do l = 1,ao_prim_num(k)
-    beta = ao_expo_ordered_transp_per_nucl(l,j,i)
-    contrib = ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
-    accu_1 += contrib
-    accu_2 += contrib * beta
-    accu_3 += contrib * beta**2
-   enddo
-   aos_array(k) = accu_1 * dx2 * dy2 * dz2
+      dy3 = 0.d0
+      if(power_ao(2) .ge. 2) then
+        dy3 = dble(power_ao(2)) * dble((power_ao(2)-1))  * dy**(power_ao(2)-2)
+      endif
+      
+      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)
 
-   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
+      ! ---
+      
+      dz1 = 0.d0
+      if(power_ao(3) .ne. 0) then
+        dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1)
+      endif
 
-   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
+      dz3 = 0.d0
+      if(power_ao(3) .ge. 2) then
+        dz3 = dble(power_ao(3)) * dble((power_ao(3)-1)) * dz**(power_ao(3)-2)
+      endif
+      
+      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)
+      
+      ! ---
+      
+      accu_1 = 0.d0
+      accu_2 = 0.d0
+      accu_3 = 0.d0
+      do l = 1,ao_prim_num(k)
+        beta = ao_expo_ordered_transp_per_nucl(l,j,i)
+        if(beta*r2.gt.50.d0) cycle
+        contrib = ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
+        accu_1 += contrib
+        accu_2 += contrib * beta
+        accu_3 += contrib * beta**2
+      enddo
 
+      aos_array(k) = accu_1 * dx2 * dy2 * dz2
+      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(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
- enddo
+
 end
 
+! ---