added some stuffs for getting the bielec integrals

src/ao_two_e_ints/map_integrals.irp.f

@@ -324,6 +324,52 @@ subroutine get_ao_two_e_integrals_non_zero_jl(j,l,thresh,sze_max,sze,out_val,out
 end
 
 
+subroutine get_ao_two_e_integrals_non_zero_jl_from_list(j,l,thresh,list,n_list,sze_max,out_val,out_val_index,non_zero_int)
+  use map_module
+  implicit none
+  BEGIN_DOC
+  ! Gets multiple AO bi-electronic integral from the AO map .
+  ! All non-zero i are retrieved for j,k,l fixed.
+  END_DOC
+  double precision, intent(in)   :: thresh
+  integer, intent(in)            :: j,l, n_list,list(2,sze_max),sze_max
+  real(integral_kind), intent(out) :: out_val(sze_max)
+  integer, intent(out)           :: out_val_index(2,sze_max),non_zero_int
+
+  integer                        :: i,k
+  integer(key_kind)              :: hash
+  double precision               :: tmp
+
+  PROVIDE ao_two_e_integrals_in_map
+  non_zero_int = 0
+  if (ao_overlap_abs(j,l) < thresh) then
+    out_val = 0.d0
+    return
+  endif
+
+  non_zero_int = 0
+ integer :: kk
+  do kk = 1, n_list
+   k = list(1,kk)
+   i = list(2,kk)
+   integer, external :: ao_l4
+   double precision, external :: ao_two_e_integral
+   !DIR$ FORCEINLINE
+   if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < thresh) then
+     cycle
+   endif
+   call two_e_integrals_index(i,j,k,l,hash)
+   call map_get(ao_integrals_map, hash,tmp)
+   if (dabs(tmp) < thresh ) cycle
+   non_zero_int = non_zero_int+1
+   out_val_index(1,non_zero_int) = i
+   out_val_index(2,non_zero_int) = k
+   out_val(non_zero_int) = tmp
+  enddo
+
+end
+
+
 
 
 function get_ao_map_size()

src/becke_numerical_grid/EZFIO.cfg

@@ -8,3 +8,9 @@ default: 2
 type: integer
 doc: Total number of grid points 
 interface: ezfio
+
+[thresh_grid]
+type: double precision
+doc: threshold on the weight of a given grid point
+interface: ezfio,provider,ocaml
+default: 1.e-20

src/becke_numerical_grid/grid_becke_per_atom.irp.f

@@ -0,0 +1,53 @@
+
+
+ BEGIN_PROVIDER [integer, n_pts_per_atom, (nucl_num)]
+&BEGIN_PROVIDER [integer, n_pts_max_per_atom]
+  BEGIN_DOC
+  ! Number of points which are non zero
+  END_DOC
+  integer                        :: i,j,k,l
+  n_pts_per_atom = 0
+  do j = 1, nucl_num
+    do i = 1, n_points_radial_grid -1
+      do k = 1, n_points_integration_angular
+        if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
+          cycle
+        endif
+        n_pts_per_atom(j) += 1
+      enddo
+    enddo
+  enddo
+  n_pts_max_per_atom = maxval(n_pts_per_atom)
+END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, final_grid_points_per_atom, (3,n_pts_max_per_atom,nucl_num)]
+&BEGIN_PROVIDER [double precision, final_weight_at_r_vector_per_atom, (n_pts_max_per_atom,nucl_num) ]
+&BEGIN_PROVIDER [integer, index_final_points_per_atom, (3,n_pts_max_per_atom,nucl_num) ]
+&BEGIN_PROVIDER [integer, index_final_points_per_atom_reverse, (n_points_integration_angular,n_points_radial_grid,nucl_num) ]
+ implicit none
+ integer                        :: i,j,k,l,i_count(nucl_num)
+ double precision               :: r(3)
+ i_count = 0
+  do j = 1, nucl_num
+    do i = 1, n_points_radial_grid -1
+      do k = 1, n_points_integration_angular
+        if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
+          cycle
+        endif
+        i_count(j) += 1
+        final_grid_points_per_atom(1,i_count(j),j) = grid_points_per_atom(1,k,i,j)
+        final_grid_points_per_atom(2,i_count(j),j) = grid_points_per_atom(2,k,i,j)
+        final_grid_points_per_atom(3,i_count(j),j) = grid_points_per_atom(3,k,i,j)
+        final_weight_at_r_vector_per_atom(i_count(j),j) = final_weight_at_r(k,i,j)
+        index_final_points_per_atom(1,i_count(j),j) = k
+        index_final_points_per_atom(2,i_count(j),j) = i
+        index_final_points_per_atom(3,i_count(j),j) = j
+        index_final_points_per_atom_reverse(k,i,j) = i_count(j)
+      enddo
+    enddo
+  enddo
+
+
+
+
+END_PROVIDER 

src/becke_numerical_grid/grid_becke_vector.irp.f

@@ -8,9 +8,9 @@ BEGIN_PROVIDER [integer, n_points_final_grid]
   do j = 1, nucl_num
     do i = 1, n_points_radial_grid -1
       do k = 1, n_points_integration_angular
-!       if(dabs(final_weight_at_r(k,i,j)) < 1.d-30)then
+        if(dabs(final_weight_at_r(k,i,j)) < tresh_grid)then
-!         cycle
+          cycle
-!       endif
+        endif
         n_points_final_grid += 1
       enddo
     enddo
@@ -39,9 +39,9 @@ END_PROVIDER
   do j = 1, nucl_num
     do i = 1, n_points_radial_grid -1
       do k = 1, n_points_integration_angular
-       !if(dabs(final_weight_at_r(k,i,j)) < 1.d-30)then
+        if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
-       !  cycle
+          cycle
-       !endif
+        endif
         i_count += 1
         final_grid_points(1,i_count) = grid_points_per_atom(1,k,i,j)
         final_grid_points(2,i_count) = grid_points_per_atom(2,k,i,j)

src/dft_utils_in_r/ao_in_r.irp.f

@@ -121,3 +121,26 @@
  enddo
  END_PROVIDER
 
+
+ BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom, (ao_num,n_pts_max_per_atom,nucl_num)]
+&BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom_transp, (n_pts_max_per_atom,ao_num,nucl_num)]
+ implicit none
+ BEGIN_DOC
+ ! aos_in_r_array_per_atom(i,j,k)        = value of the ith ao on the jth grid point attached on the kth atom 
+ END_DOC
+ integer :: i,j,k
+ double precision :: aos_array(ao_num), r(3)
+ do k = 1, nucl_num
+  do i = 1, n_pts_per_atom(k)
+   r(1) = final_grid_points_per_atom(1,i,k)
+   r(2) = final_grid_points_per_atom(2,i,k)
+   r(3) = final_grid_points_per_atom(3,i,k)
+   call give_all_aos_at_r(r,aos_array)
+   do j = 1, ao_num
+    aos_in_r_array_per_atom(j,i,k) = aos_array(j)
+    aos_in_r_array_per_atom_transp(i,j,k) = aos_array(j)
+   enddo
+  enddo
+ enddo
+ END_PROVIDER
+