anchor.diff:  
   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
  
positive.code:   integer, intent(in)            :: i,j,k,l
    integer(key_kind)              :: idx
    type(map_type), intent(inout)  :: map
    real(integral_kind)            :: tmp
    PROVIDE ao_bielec_integrals_in_map
    !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
        !DIR$ FORCEINLINE
        call bielec_integrals_index(i,j,k,l,idx)
        call map_get(map,idx,tmp)
      endif
    endif
    get_ao_bielec_integral = tmp
  end
  
  
  
negative.code:   ! Returns one integral <ij|kl> in the MO basis
    END_DOC
    integer, intent(in)            :: i,j,k,l
    type(map_type), intent(inout)  :: map
    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
      !DIR$ FORCEINLINE
      get_mo_bielec_integral_schwartz = get_mo_bielec_integral(i,j,k,l,map)
    else
      get_mo_bielec_integral_schwartz = 0.d0
    endif
  end
  
  