Added key -> i,j,k,l function in maps

src/AOs/aos.irp.f

@@ -19,11 +19,13 @@ END_PROVIDER
  BEGIN_PROVIDER [ integer, ao_power, (ao_num_align,3) ]
 &BEGIN_PROVIDER [ double precision, ao_expo, (ao_num_align,ao_prim_num_max) ]
 &BEGIN_PROVIDER [ double precision, ao_coef, (ao_num_align,ao_prim_num_max) ]
+&BEGIN_PROVIDER [ integer, ao_l, (ao_num) ]
  implicit none
 
  BEGIN_DOC
 ! Coefficients, exponents and powers of x,y and z
 ! ao_coef(i,j) = coefficient of the jth primitive on the ith ao
+! ao_l = l value of the AO: a+b+c in x^a y^b z^c
  END_DOC
  PROVIDE ezfio_filename
 
@@ -91,8 +93,13 @@ END_PROVIDER
    ao_coef(i,j) = d(j,2)
   enddo
  enddo
+ do i=1,ao_num
+   ao_l(i) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3) 
+ enddo
 END_PROVIDER
 
+
+
  BEGIN_PROVIDER [ double precision, ao_coef_transp, (ao_prim_num_max_align,ao_num) ]
 &BEGIN_PROVIDER [ double precision, ao_expo_transp, (ao_prim_num_max_align,ao_num) ]
  implicit none

src/BiInts/ao_bi_integrals.irp.f

@@ -139,8 +139,14 @@ double precision function ao_bielec_integral(i,j,k,l)
 end
 
 
-
-
+integer function ao_l4(i,j,k,l)
+  implicit none
+  BEGIN_DOC
+! Computes the product of l values of i,j,k,and l
+  END_DOC
+  integer, intent(in) :: i,j,k,l
+  ao_l4 = ao_l(i)*ao_l(j)*ao_l(k)*ao_l(l)
+end
 
 
 

src/BiInts/map_integrals.irp.f

@@ -30,6 +30,26 @@ subroutine bielec_integrals_index(i,j,k,l,i1)
   i1 = i1+ishft(i2*i2-i2,-1)
 end
 
+subroutine bielec_integrals_index_reverse(i,j,k,l,i1)
+  implicit none
+  integer, intent(out)           :: i,j,k,l
+  integer*8, intent(in)          :: i1
+  integer*8                      :: i2,i3
+  real                           :: x
+  x = 0.5*(sqrt(8.*real(i1)+1.)-1.)
+  i2 = ceiling(x)
+  i3 = i1 - ishft(i2*i2-i2,-1)
+
+  l = 0.5*(sqrt(8.*real(i2)+1.)-1.)
+  l = ceiling(x)
+  j = i2 - ishft(l*l-l,-1)
+
+  x = 0.5*(sqrt(8.*real(i3)+1.)-1.)
+  k = ceiling(x)
+  i = i3 - ishft(k*k-k,-1)
+end
+
+
 double precision function get_ao_bielec_integral(i,j,k,l,map)
   use map_module
   implicit none

src/Utils/map_module.f90

@@ -1,5 +1,17 @@
 module map_module
 
+! A map is an array of maps (cache_maps)
+! A cache map is an array of keys and values sorted by keys
+! A cache map has its own OpenMP lock
+! To access a (key,value) pair in the map, the
+! index of the cache_map in the map array is obtained
+! by removing the first 15 bits of the key. 
+! The key in the cache_map is composed of the first
+! 15 bits of the key. Therefore, it can be stored
+! as integer*2 and is found by applying the map_mask
+! to the initial key. The element are found in the
+! cache_map using a binary search
+
  use omp_lib
  
  integer, parameter             :: integral_kind = 8
@@ -42,7 +54,7 @@ real function map_mb(map)
   
   map_mb = 8+map_size_kind+map_size_kind+omp_lock_kind+4
   do i=0,map%map_size
-   map_mb = map_mb + map%map(i)%map_size*(cache_key_kind+integral_kind) + &
+    map_mb = map_mb + map%map(i)%map_size*(cache_key_kind+integral_kind) +&
         8+8+4+cache_map_size_kind+cache_map_size_kind+omp_lock_kind
   enddo
   map_mb = map_mb / (1024.d0*1024.d0)
@@ -303,7 +315,7 @@ subroutine map_unique(map)
   integer(map_size_kind)         :: icount
   
   icount = 0_8
- !$OMP PARALLEL DO SCHEDULE(dynamic,1000) DEFAULT(SHARED) PRIVATE(i) &
+  !$OMP PARALLEL DO SCHEDULE(dynamic,1000) DEFAULT(SHARED) PRIVATE(i)&
       !$OMP REDUCTION(+:icount)
   do i=0_8,map%map_size
     call omp_set_lock(map%map(i)%lock)
@@ -325,7 +337,7 @@ subroutine map_shrink(map,thr)
   integer(map_size_kind)         :: icount
   
   icount = 0_8
- !$OMP PARALLEL DO SCHEDULE(dynamic,1000) DEFAULT(SHARED) PRIVATE(i) &
+  !$OMP PARALLEL DO SCHEDULE(dynamic,1000) DEFAULT(SHARED) PRIVATE(i)&
       !$OMP REDUCTION(+:icount)
   do i=0_8,map%map_size
     call omp_set_lock(map%map(i)%lock)
@@ -344,7 +356,7 @@ subroutine map_update(map, key, value, sze, thr)
   type (map_type), intent(inout) :: map
   integer, intent(in)            :: sze
   integer(key_kind), intent(inout) :: key(sze)
- real(integral_kind), intent(inout):: value(sze)
+  real(integral_kind), intent(inout) :: value(sze)
   real(integral_kind), intent(in) :: thr
   
   integer                        :: i
@@ -413,11 +425,11 @@ subroutine map_update(map, key, value, sze, thr)
       i=i+1
     endif  ! key = 0
   enddo  ! i
- enddo  ! sze2 > 0
- call omp_set_lock(map%lock)
- map%n_elements = map%n_elements + n_elements_temp
- map%sorted = map%sorted .and. map_sorted
- call omp_unset_lock(map%lock)
+enddo  ! sze2 > 0
+call omp_set_lock(map%lock)
+map%n_elements = map%n_elements + n_elements_temp
+map%sorted = map%sorted .and. map_sorted
+call omp_unset_lock(map%lock)
 
 end
 
@@ -427,7 +439,7 @@ subroutine map_append(map, key, value, sze)
   type (map_type), intent(inout) :: map
   integer, intent(in)            :: sze
   integer(key_kind), intent(inout) :: key(sze)
- real(integral_kind), intent(inout):: value(sze)
+  real(integral_kind), intent(inout) :: value(sze)
   
   integer                        :: i
   integer(cache_map_size_kind)   :: n_elements
@@ -526,7 +538,7 @@ subroutine map_exists_many(map, key, sze)
   implicit none
   type (map_type), intent(inout) :: map
   integer, intent(in)            :: sze
- integer(key_kind), intent(inout):: key(sze)
+  integer(key_kind), intent(inout) :: key(sze)
   integer                        :: i
   integer(map_size_kind)         :: idx_cache, idx_cache_prev
   integer(cache_map_size_kind)   :: ibegin, iend
@@ -678,3 +690,38 @@ subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in)
 end
 
 
+subroutine get_cache_map_n_elements_max(map,n_elements_max)
+  use map_module
+  implicit none
+  ! Returns the size of the largest cache_map
+  type (map_type), intent(in)    :: map
+  integer(cache_map_size_kind), intent(out) :: n_elements_max
+  integer(map_size_kind)         :: i
+  n_elements_max = 0_cache_map_size_kind
+  do i=1,map%map_size
+    n_elements_max = max(n_elements_max, map%map(i)%n_elements)
+  enddo
+end
+
+
+
+subroutine get_cache_map(map,map_idx,keys,values,n_elements)
+  use map_module
+  implicit none
+  type (map_type), intent(in)    :: map
+  integer(map_size_kind), intent(in) :: map_idx
+  integer(cache_map_size_kind), intent(inout) :: n_elements
+  integer(key_kind), intent(out) :: keys(n_elements)
+  double precision, intent(out)  :: values(n_elements)
+  integer(cache_map_size_kind)   :: i
+  integer(key_kind)              :: shift
+  
+  shift = ishft(map_idx,-map_shift)
+  
+  n_elements = map%map(map_idx)%n_elements
+  do i=1,n_elements
+    keys(i) = map%map(map_idx)%key(i) + shift
+    values(i) = map%map(map_idx)%value(i)
+  enddo
+  
+end