diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f
index afc573aa..22fd48a6 100644
--- a/src/Integrals_Bielec/map_integrals.irp.f
+++ b/src/Integrals_Bielec/map_integrals.irp.f
@@ -120,15 +120,16 @@ end
 
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, ao_integrals_cache, (ao_integrals_cache_min:ao_integrals_cache_max,ao_integrals_cache_min:ao_integrals_cache_max,ao_integrals_cache_min:ao_integrals_cache_max,ao_integrals_cache_min:ao_integrals_cache_max) ]
+BEGIN_PROVIDER [ double precision, ao_integrals_cache, (0:64*64*64*64) ]
  implicit none
  BEGIN_DOC
  ! Cache of AO integrals for fast access
  END_DOC
  PROVIDE ao_bielec_integrals_in_map
- integer                        :: i,j,k,l
+ integer                        :: i,j,k,l,ii
  integer(key_kind)              :: idx
- !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx)
+ real(integral_kind)            :: integral
+ !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
  do l=ao_integrals_cache_min,ao_integrals_cache_max
    do k=ao_integrals_cache_min,ao_integrals_cache_max
      do j=ao_integrals_cache_min,ao_integrals_cache_max
@@ -136,7 +137,12 @@ BEGIN_PROVIDER [ double precision, ao_integrals_cache, (ao_integrals_cache_min:a
          !DIR$ FORCEINLINE
          call bielec_integrals_index(i,j,k,l,idx)
          !DIR$ FORCEINLINE
-         call map_get(ao_integrals_map,idx,ao_integrals_cache(i,j,k,l))
+         call map_get(ao_integrals_map,idx,integral)
+         ii = l-ao_integrals_cache_min
+         ii = ior( ishft(ii,6), k-ao_integrals_cache_min)
+         ii = ior( ishft(ii,6), j-ao_integrals_cache_min)
+         ii = ior( ishft(ii,6), i-ao_integrals_cache_min)
+         ao_integrals_cache(ii) = integral
        enddo
      enddo
    enddo
@@ -155,30 +161,33 @@ double precision function get_ao_bielec_integral(i,j,k,l,map)
   integer, intent(in)            :: i,j,k,l
   integer(key_kind)              :: idx
   type(map_type), intent(inout)  :: map
+  integer                        :: ii
   real(integral_kind)            :: tmp
-  PROVIDE ao_bielec_integrals_in_map
+  PROVIDE ao_bielec_integrals_in_map ao_integrals_cache ao_integrals_cache_min
   !DIR$ FORCEINLINE
   if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then
     tmp = 0.d0
   else if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < ao_integrals_threshold) then
     tmp = 0.d0
   else
-    if ( (i >= ao_integrals_cache_min) .and. &
-        (j >= ao_integrals_cache_min) .and. &
-        (k >= ao_integrals_cache_min) .and. &
-        (l >= ao_integrals_cache_min) .and. &
-        (i <= ao_integrals_cache_max) .and. &
-        (j <= ao_integrals_cache_max) .and. &
-        (k <= ao_integrals_cache_max) .and. &
-        (l <= ao_integrals_cache_max) ) then
-        tmp = ao_integrals_cache(i,j,k,l)
-    else
+    ii = l-ao_integrals_cache_min
+    ii = ior(ii, k-ao_integrals_cache_min)
+    ii = ior(ii, j-ao_integrals_cache_min)
+    ii = ior(ii, i-ao_integrals_cache_min)
+    if (iand(ii, -64) /= 0) then
       !DIR$ FORCEINLINE
       call bielec_integrals_index(i,j,k,l,idx)
+      !DIR$ FORCEINLINE
       call map_get(map,idx,tmp)
+      get_ao_bielec_integral = dble(tmp)
+    else
+      ii = l-ao_integrals_cache_min
+      ii = ior( ishft(ii,6), k-ao_integrals_cache_min)
+      ii = ior( ishft(ii,6), j-ao_integrals_cache_min)
+      ii = ior( ishft(ii,6), i-ao_integrals_cache_min)
+      get_ao_bielec_integral = ao_integrals_cache(ii)
     endif
   endif
-  get_ao_bielec_integral = tmp
 end
 
 
@@ -324,20 +333,22 @@ end
  ! Min and max values of the MOs for which the integrals are in the cache
  END_DOC
  mo_integrals_cache_min = max(1,elec_alpha_num - 31)
- mo_integrals_cache_max = min(mo_tot_num,elec_alpha_num + 32)
+ mo_integrals_cache_max = min(mo_tot_num,mo_integrals_cache_min+63)
 
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, mo_integrals_cache, (mo_integrals_cache_min:mo_integrals_cache_max,mo_integrals_cache_min:mo_integrals_cache_max,mo_integrals_cache_min:mo_integrals_cache_max,mo_integrals_cache_min:mo_integrals_cache_max) ]
+BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0:64*64*64*64) ]
  implicit none
  BEGIN_DOC
  ! Cache of MO integrals for fast access
  END_DOC
  PROVIDE mo_bielec_integrals_in_map
  integer                        :: i,j,k,l
+ integer                        :: ii
  integer(key_kind)              :: idx
+ real(integral_kind)            :: integral
  FREE ao_integrals_cache
- !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx)
+ !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
  do l=mo_integrals_cache_min,mo_integrals_cache_max
    do k=mo_integrals_cache_min,mo_integrals_cache_max
      do j=mo_integrals_cache_min,mo_integrals_cache_max
@@ -345,7 +356,12 @@ BEGIN_PROVIDER [ double precision, mo_integrals_cache, (mo_integrals_cache_min:m
          !DIR$ FORCEINLINE
          call bielec_integrals_index(i,j,k,l,idx)
          !DIR$ FORCEINLINE
-         call map_get(mo_integrals_map,idx,mo_integrals_cache(i,j,k,l))
+         call map_get(mo_integrals_map,idx,integral)
+         ii = l-mo_integrals_cache_min
+         ii = ior( ishft(ii,6), k-mo_integrals_cache_min)
+         ii = ior( ishft(ii,6), j-mo_integrals_cache_min)
+         ii = ior( ishft(ii,6), i-mo_integrals_cache_min)
+         mo_integrals_cache(ii) = integral
        enddo
      enddo
    enddo
@@ -362,24 +378,26 @@ double precision function get_mo_bielec_integral(i,j,k,l,map)
   END_DOC
   integer, intent(in)            :: i,j,k,l
   integer(key_kind)              :: idx
+  integer                        :: ii
   type(map_type), intent(inout)  :: map
   real(integral_kind)            :: tmp
   PROVIDE mo_bielec_integrals_in_map mo_integrals_cache
-  if ( (i >= mo_integrals_cache_min) .and. &
-       (j >= mo_integrals_cache_min) .and. &
-       (k >= mo_integrals_cache_min) .and. &
-       (l >= mo_integrals_cache_min) .and. &
-       (i <= mo_integrals_cache_max) .and. &
-       (j <= mo_integrals_cache_max) .and. &
-       (k <= mo_integrals_cache_max) .and. &
-       (l <= mo_integrals_cache_max) ) then
-       get_mo_bielec_integral = mo_integrals_cache(i,j,k,l)
-  else
+  ii = l-mo_integrals_cache_min
+  ii = ior(ii, k-mo_integrals_cache_min)
+  ii = ior(ii, j-mo_integrals_cache_min)
+  ii = ior(ii, i-mo_integrals_cache_min)
+  if (iand(ii, -64) /= 0) then
     !DIR$ FORCEINLINE
     call bielec_integrals_index(i,j,k,l,idx)
     !DIR$ FORCEINLINE
     call map_get(map,idx,tmp)
     get_mo_bielec_integral = dble(tmp)
+  else
+    ii = l-mo_integrals_cache_min
+    ii = ior( ishft(ii,6), k-mo_integrals_cache_min)
+    ii = ior( ishft(ii,6), j-mo_integrals_cache_min)
+    ii = ior( ishft(ii,6), i-mo_integrals_cache_min)
+    get_mo_bielec_integral = mo_integrals_cache(ii)
   endif
 end
 
@@ -390,16 +408,15 @@ double precision function get_mo_bielec_integral_schwartz(i,j,k,l,map)
   ! Returns one integral <ij|kl> in the MO basis
   END_DOC
   integer, intent(in)            :: i,j,k,l
-  integer(key_kind)              :: idx
   type(map_type), intent(inout)  :: map
-  real(integral_kind)            :: tmp
+  double precision, external     :: get_mo_bielec_integral
+
   PROVIDE mo_bielec_integrals_in_map  mo_integrals_cache
   if (mo_bielec_integral_schwartz(i,k)*mo_bielec_integral_schwartz(j,l) > mo_integrals_threshold) then
-    double precision, external :: get_mo_bielec_integral
     !DIR$ FORCEINLINE
     get_mo_bielec_integral_schwartz = get_mo_bielec_integral(i,j,k,l,map)
   else
-    tmp = 0.d0
+    get_mo_bielec_integral_schwartz = 0.d0
   endif
 end