diff --git a/src/BiInts/ao_bi_integrals.irp.f b/src/BiInts/ao_bi_integrals.irp.f
index 272f5ad8..3093243b 100644
--- a/src/BiInts/ao_bi_integrals.irp.f
+++ b/src/BiInts/ao_bi_integrals.irp.f
@@ -215,7 +215,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
   integer                        :: jl_pairs(2,ao_num*(ao_num+1)/2), kk, m, j1, i1, lmax
   
   PROVIDE gauleg_t2 ao_integrals_map all_utils
-  PROVIDE ao_bielec_integral_schwartz
+!  PROVIDE ao_bielec_integral_schwartz
   integral = ao_bielec_integral(1,1,1,1)
   
   real                           :: map_mb
@@ -357,7 +357,7 @@ END_PROVIDER
 BEGIN_PROVIDER [ double precision, ao_bielec_integral_schwartz,(ao_num,ao_num)  ]
   implicit none
   BEGIN_DOC
-  !  Needed to compuet Schwartz inequalities
+  !  Needed to compute Schwartz inequalities
   END_DOC
   
   integer                        :: i,k
diff --git a/src/BiInts/map_integrals.irp.f b/src/BiInts/map_integrals.irp.f
index d82c88df..dc650294 100644
--- a/src/BiInts/map_integrals.irp.f
+++ b/src/BiInts/map_integrals.irp.f
@@ -9,6 +9,7 @@ BEGIN_PROVIDER [ type(map_type), ao_integrals_map ]
   ! AO integrals
   END_DOC
   integer*8                      :: sze
+  PROVIDE gauleg_t2
   call bielec_integrals_index(ao_num,ao_num,ao_num,ao_num,sze)
   call map_init(ao_integrals_map,sze)
   write(output_BiInts,*)  'AO map initialized'
diff --git a/src/Hartree_Fock/Fock_matrix.irp.f b/src/Hartree_Fock/Fock_matrix.irp.f
index 6802d0fb..a1553ebd 100644
--- a/src/Hartree_Fock/Fock_matrix.irp.f
+++ b/src/Hartree_Fock/Fock_matrix.irp.f
@@ -107,44 +107,73 @@ END_PROVIDER
  ! Alpha Fock matrix in AO basis set
  END_DOC
  
- integer                        :: i,j,k,l,k1
+ integer                        :: i,j,k,l,k1,r,s
+ integer*8                      :: p,q
  double precision               :: integral
  double precision               :: ao_bielec_integral
  if (do_direct_integrals) then
-   PROVIDE all_utils ao_overlap_abs ao_integrals_threshold
+   PROVIDE all_utils ao_overlap_abs ao_integrals_threshold gauleg_t2
    PROVIDE HF_density_matrix_ao_alpha  HF_density_matrix_ao_beta
    PROVIDE ao_bi_elec_integral_alpha
-   !$OMP PARALLEL DEFAULT(NONE) &
-   !$OMP PRIVATE(i,j,l,k1,k,integral) &
-   !$OMP SHARED(ao_num,ao_bi_elec_integral_alpha,ao_mono_elec_integral,&
-   !$OMP  ao_num_align,ao_bi_elec_integral_beta,HF_density_matrix_ao_alpha, &
-   !$OMP  HF_density_matrix_ao_beta) 
-   !$OMP DO SCHEDULE(GUIDED)
-   do j=1,ao_num
-     do i=1,j
-       ao_bi_elec_integral_alpha(i,j) = 0.d0
-       ao_bi_elec_integral_beta (i,j) = 0.d0
-       do l=1,ao_num
-         do k=1,ao_num
-           if ((abs(HF_density_matrix_ao_alpha(k,l)) > 1.d-9).or.    &
-                 (abs(HF_density_matrix_ao_beta (k,l)) > 1.d-9)) then
-             integral = (HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta (k,l)) * ao_bielec_integral(k,l,i,j)
+
+   ao_bi_elec_integral_alpha = 0.d0
+   ao_bi_elec_integral_beta  = 0.d0
+   !$OMP PARALLEL DEFAULT(NONE)                                      &
+       !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,p,q,r,s)&
+       !$OMP SHARED(ao_num,HF_density_matrix_ao_alpha,HF_density_matrix_ao_beta,&
+       !$OMP ao_integrals_map,ao_integrals_threshold, ao_bielec_integral_schwartz, &
+       !$OMP ao_overlap_abs) &
+       !$OMP REDUCTION(+:ao_bi_elec_integral_alpha,ao_bi_elec_integral_beta)
+
+   allocate(keys(1), values(1))
+
+   q = ao_num*ao_num*ao_num*ao_num
+   !$OMP DO SCHEDULE(dynamic)
+   do p=1_8,q
+           call bielec_integrals_index_reverse(kk,ii,ll,jj,p)
+           if ( (kk(1)>ao_num).or. &
+                (ii(1)>ao_num).or. &
+                (jj(1)>ao_num).or. &
+                (ll(1)>ao_num) ) then
+                cycle
+           endif
+           k = kk(1)
+           i = ii(1)
+           l = ll(1)
+           j = jj(1)
+
+           if (ao_overlap_abs(k,l)*ao_overlap_abs(i,j)  &
+              < ao_integrals_threshold) then
+             cycle
+           endif
+           if (ao_bielec_integral_schwartz(k,l)*ao_bielec_integral_schwartz(i,j)  &
+              < ao_integrals_threshold) then
+             cycle
+           endif
+           values(1) = ao_bielec_integral(k,l,i,j)
+!           values(1) = ao_bielec_integral(k,i,l,j)
+           if (abs(values(1)) < ao_integrals_threshold) then
+             cycle
+           endif
+           do k2=1,8
+             if (kk(k2)==0) then
+               cycle
+             endif
+             i = ii(k2)
+             j = jj(k2)
+             k = kk(k2)
+             l = ll(k2)
+             integral = (HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l)) * values(1)
              ao_bi_elec_integral_alpha(i,j) += integral
              ao_bi_elec_integral_beta (i,j) += integral
-
-             integral = ao_bielec_integral(k,j,i,l)
-             ao_bi_elec_integral_alpha(i,j) -= HF_density_matrix_ao_alpha(k,l)*integral
-             ao_bi_elec_integral_beta (i,j) -= HF_density_matrix_ao_beta (k,l)*integral
-           endif
-         enddo
-       enddo
-       ao_bi_elec_integral_alpha(j,i) = ao_bi_elec_integral_alpha(i,j)
-       ao_bi_elec_integral_beta (j,i) = ao_bi_elec_integral_beta (i,j)
-     enddo
+             integral = values(1)
+             ao_bi_elec_integral_alpha(l,j) -= HF_density_matrix_ao_alpha(k,i) * integral
+             ao_bi_elec_integral_beta (l,j) -= HF_density_matrix_ao_beta (k,i) * integral
+           enddo
    enddo
-   !$OMP END DO NOWAIT
+   !$OMP END DO
+   deallocate(keys,values)
    !$OMP END PARALLEL
-
  else
    PROVIDE ao_bielec_integrals_in_map 
            
@@ -166,7 +195,6 @@ END_PROVIDER
    call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max)
    allocate(keys(n_elements_max), values(n_elements_max))
 
-   !$OMP BARRIER
    !$OMP DO SCHEDULE(dynamic)
    do i8=0_8,ao_integrals_map%map_size
      n_elements = n_elements_max