Merge branch 'dev-stable' of https://github.com/QuantumPackage/qp2 into dev-stable

plugins/local/non_h_ints_mu/NEED

@@ -3,3 +3,4 @@ hamiltonian
 jastrow
 ao_tc_eff_map
 bi_ortho_mos
+trexio

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
 

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
+
+

src/ao_basis/aos_in_r.irp.f

@@ -1,29 +1,35 @@
-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 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
- integer :: power_ao(3)
+  double precision :: accu, dx, dy, dz, r2
- double precision :: accu,dx,dy,dz,r2
+
   num_ao = ao_nucl(i)
- power_ao(1:3)= ao_power(i,1:3)
+  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))
+  dx = r(1) - center_ao(1)
- dy = (r(2) - center_ao(2))
+  dy = r(2) - center_ao(2)
- dz = (r(3) - center_ao(3))
+  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)
+  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
@@ -32,33 +38,36 @@ double precision function ao_value(i,r)
 end
 
 
-double precision function primitive_value(i,j,r)
+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
 
- integer :: m,num_ao
+  BEGIN_DOC
+  ! Returns the value of the j-th primitive of the i-th |AO| at point $\textbf{r}
+  ! **without the coefficient**
+  END_DOC
+
+  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
- integer :: power_ao(3)
+  double precision :: accu, dx, dy, dz, r2
- 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))
+  dx = r(1) - center_ao(1)
- dy = (r(2) - center_ao(2))
+  dy = r(2) - center_ao(2)
- dz = (r(3) - center_ao(3))
+  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
+  m = j
   beta = ao_expo_ordered_transp(m,i)
   accu += dexp(-beta*r2)
   primitive_value = accu * dx * dy * dz
@@ -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,116 +129,138 @@ end
 
 ! ---
 
-subroutine give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
+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
+  !
-!
+  ! input : r(1) ==> r(1) = x, r(2) = y, r(3) = z
-! output : 
+  !
-!
+  ! output : 
-! * aos_array(i) = ao(i) evaluated at ro
+  !
-! * aos_grad_array(1,i) = gradient X of the ao(i) evaluated at $\textbf{r}$
+  ! * 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
+  integer                       :: i, j, k, l, m
- double precision :: dx,dy,dz,r2
+  double precision              :: dx, dy, dz, r2
- double precision ::      dx2,dy2,dz2
+  double precision              :: dx1, dy1, dz1
- double precision ::      dx1,dy1,dz1
+  double precision              :: dx2, dy2, dz2
   double precision              :: center_ao(3)
- double precision :: beta,accu_1,accu_2,contrib
+  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))
+    dx = r(1) - center_ao(1)
-  dz = (r(3) - center_ao(3))
+    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)
+
+    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)
+
+      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
+      if(power_ao(1) .ne. 0) then
+        dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1)
       endif
-   if(power_ao(2) .ne. 0)then
+
-    dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1)
-   else
       dy1 = 0.d0
+      if(power_ao(2) .ne. 0) then
+        dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1)
       endif
-   if(power_ao(3) .ne. 0)then
+
-    dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1)
-   else
       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)
+      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
-    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(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(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_grad_array(3,k) = accu_1 * dx2 * dy2 * dz1 - 2.d0 * dx2 * dy2 * dz2 * dz  * accu_2
     enddo
   enddo
+
 end
 
+! ---
+
+subroutine give_all_aos_and_grad_and_lapl_at_r(r, aos_array, aos_grad_array, aos_lapl_array)
 
-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
+  !
-!
+  ! input  : r(1) ==> r(1) = x, r(2) = y, r(3) = z
-! output :
+  !
-!
+  ! 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}$
+  ! * 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
+
+  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)
 
   integer                       :: power_ao(3)
- integer :: i,j,k,l,m
+  integer                       :: i, j, k, l, m
- double precision :: dx,dy,dz,r2
+  double precision              :: dx, dy, dz, r2
- double precision ::      dx2,dy2,dz2
+  double precision              :: dx1, dy1, dz1
- double precision ::      dx1,dy1,dz1
+  double precision              :: dx2, dy2, dz2
- double precision ::      dx3,dy3,dz3
+  double precision              :: dx3, dy3, dz3
- double precision ::      dx4,dy4,dz4
+  double precision              :: dx4, dy4, dz4
- double precision ::      dx5,dy5,dz5
+  double precision              :: dx5, dy5, dz5
   double precision              :: center_ao(3)
- double precision :: beta,accu_1,accu_2,accu_3,contrib
+  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))
+    dx = r(1) - center_ao(1)
-  dz = (r(3) - center_ao(3))
+    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)
+    
+    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      
-
       aos_lapl_array(1,k) = 0.d0
       aos_lapl_array(2,k) = 0.d0
       aos_lapl_array(3,k) = 0.d0
@@ -238,18 +269,20 @@ subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_la
       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
+      if(power_ao(1) .ne. 0) then
+        dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1)
       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
+      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
+
+      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)) 
@@ -257,19 +290,19 @@ subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_la
       
       dx5 = dx**(power_ao(1)+2)
   
-   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(2) .ne. 0) then
-   ! For the Laplacian
+        dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1)
-   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
 
-   if(power_ao(2) .ge. 1)then
+      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)) 
@@ -277,20 +310,19 @@ subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_la
       
       dy5 = dy**(power_ao(2)+2)
 
+      ! ---
       
-   if(power_ao(3) .ne. 0)then
-    dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1)
-   else
       dz1 = 0.d0
-   endif
+      if(power_ao(3) .ne. 0) then
-   ! For the Laplacian
+        dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1)
-   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
+      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)) 
@@ -298,29 +330,31 @@ subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_la
       
       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(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(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_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(1,k) = accu_1 * dx3  * dy2 * dz2- 2.d0 * dx4 * dy2 * dz2* accu_2 +4.d0 * dx5 *dy2 * 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
-   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
+
 end
 
+! ---