diff --git a/src/dft_utils_func/on_top_from_ueg.irp.f b/src/dft_utils_func/on_top_from_ueg.irp.f
index 4e28ad89..711ffc39 100644
--- a/src/dft_utils_func/on_top_from_ueg.irp.f
+++ b/src/dft_utils_func/on_top_from_ueg.irp.f
@@ -32,7 +32,6 @@ double precision function g0_UEG_mu_inf(rho_a,rho_b)
  C = 0.08193d0
  D = -0.01277d0
  E = 0.001859d0                     
- x = -d2*rs
  if (dabs(rho) > 1.d-20) then
   rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19
   x = -d2*rs
diff --git a/src/mo_two_e_ints/cholesky.irp.f b/src/mo_two_e_ints/cholesky.irp.f
index 349f13b9..0d0989d7 100644
--- a/src/mo_two_e_ints/cholesky.irp.f
+++ b/src/mo_two_e_ints/cholesky.irp.f
@@ -34,8 +34,10 @@ BEGIN_PROVIDER [ double precision, cholesky_mo_transp, (cholesky_mo_num, mo_num,
  END_DOC
 
  double precision, allocatable :: X(:,:,:)
+ double precision :: wall0, wall1
  integer :: ierr
  print *, 'AO->MO Transformation of Cholesky vectors'
+  call wall_time(wall0)
 
  allocate(X(mo_num,cholesky_mo_num,ao_num), stat=ierr)
  if (ierr /= 0) then
@@ -46,6 +48,8 @@ BEGIN_PROVIDER [ double precision, cholesky_mo_transp, (cholesky_mo_num, mo_num,
  call dgemm('T','N', cholesky_mo_num*mo_num, mo_num, ao_num, 1.d0, &
      X, ao_num, mo_coef, ao_num, 0.d0, cholesky_mo_transp, cholesky_mo_num*mo_num)
  deallocate(X)
+  call wall_time(wall1)
+ print*,'Time for AO->MO Cholesky vectors = ',wall1-wall0
 
 END_PROVIDER
 
diff --git a/src/mu_of_r/f_hf_cholesky.irp.f b/src/mu_of_r/f_hf_cholesky.irp.f
index 17f0229a..472abb1c 100644
--- a/src/mu_of_r/f_hf_cholesky.irp.f
+++ b/src/mu_of_r/f_hf_cholesky.irp.f
@@ -199,7 +199,7 @@ BEGIN_PROVIDER [ double precision, f_hf_cholesky_sparse, (n_points_final_grid)]
  !! V_AR = \sum_{I}Phi_IR V_AI = \sum_{I}Phi^t_RI V_AI
  double precision :: u_dot_v,wall0,wall1,accu_1, accu_2,mo_i_r1,mo_b_r1
  double precision :: thresh_1,thresh_2
- double precision, allocatable :: accu_vec(:)
+ double precision, allocatable :: accu_vec(:),delta_vec(:)
  thresh_2 = ao_cholesky_threshold * 100.d0
  thresh_1 = dsqrt(thresh_2)
  provide cholesky_mo_transp
@@ -223,12 +223,12 @@ BEGIN_PROVIDER [ double precision, f_hf_cholesky_sparse, (n_points_final_grid)]
        mo_b_r1 = mos_in_r_array_omp(m,ipoint)
        if(dabs(mo_i_r1*mo_b_r1).lt.thresh_2)cycle
        do p = 1, cholesky_mo_num
-        accu_vec(p) += mo_i_r1 * mo_b_r1 * cholesky_mo_transp(p,m,i)
+        accu_vec(p) = accu_vec(p) + mo_i_r1 * mo_b_r1 * cholesky_mo_transp(p,m,i)
        enddo
       enddo
      enddo
      do p = 1, cholesky_mo_num
-      f_hf_cholesky_sparse(ipoint) += accu_vec(p) * accu_vec(p)
+      f_hf_cholesky_sparse(ipoint) = f_hf_cholesky_sparse(ipoint) + accu_vec(p) * accu_vec(p)
      enddo
     f_hf_cholesky_sparse(ipoint) *= 2.D0
    enddo
@@ -240,39 +240,50 @@ BEGIN_PROVIDER [ double precision, f_hf_cholesky_sparse, (n_points_final_grid)]
   print*,'Time to provide f_hf_cholesky_sparse = ',wall1-wall0
  else
   call wall_time(wall0)
-  !$OMP PARALLEL DO &
-  !$OMP DEFAULT (NONE)  &
-  !$OMP PRIVATE (accu_2,accu_1,ipoint,p,ii,i,mm,m,mo_i_r1,mo_b_r1) & 
-  !$OMP ShARED (n_occ_val_orb_for_hf,list_valence_orb_for_hf,list_basis,mos_in_r_array_omp) & 
-  !$OMP ShARED (cholesky_mo_num,f_hf_cholesky_sparse,n_points_final_grid,cholesky_mo,n_basis_orb) 
+  !$OMP PARALLEL DEFAULT(NONE)                                      &
+  !$OMP PRIVATE (accu_vec,delta_vec,ipoint,p,ii,i,mm,m,mo_i_r1,mo_b_r1) & 
+  !$OMP ShARED (n_occ_val_orb_for_hf,list_valence_orb_for_hf,list_basis,mos_in_r_array_omp,thresh_1,thresh_2) & 
+  !$OMP ShARED (cholesky_mo_num,f_hf_cholesky_sparse,n_points_final_grid,cholesky_mo_transp,n_basis_orb) 
+  allocate(accu_vec(cholesky_mo_num),delta_vec(cholesky_mo_num))
+  !$OMP DO 
    do ipoint = 1, n_points_final_grid
     f_hf_cholesky_sparse(ipoint) = 0.d0
-    do p = 1, cholesky_mo_num
-     accu_2 = 0.d0
+     accu_vec = 0.d0
      do ii = 1, n_occ_val_orb_for_hf(2)
       i = list_valence_orb_for_hf(ii,2)
       mo_i_r1 = mos_in_r_array_omp(i,ipoint)
+      if(dabs(mo_i_r1).lt.thresh_1)cycle
       do mm = 1, n_basis_orb ! electron 1 
        m = list_basis(mm)
        mo_b_r1 = mos_in_r_array_omp(m,ipoint)
-       accu_2 += mo_i_r1 * mo_b_r1 * cholesky_mo(m,i,p)
+       if(dabs(mo_i_r1*mo_b_r1).lt.thresh_2)cycle
+       do p = 1, cholesky_mo_num
+        accu_vec(p) = accu_vec(p) + mo_i_r1 * mo_b_r1 * cholesky_mo_transp(p,m,i)
+       enddo
       enddo
      enddo
-     accu_1 = accu_2
-     do ii = n_occ_val_orb_for_hf(2)+1,n_occ_val_orb_for_hf(1)
+     delta_vec = 0.d0
+     do ii =  n_occ_val_orb_for_hf(2)+1,n_occ_val_orb_for_hf(1)
       i = list_valence_orb_for_hf(ii,1)
       mo_i_r1 = mos_in_r_array_omp(i,ipoint)
+      if(dabs(mo_i_r1).lt.thresh_1)cycle
       do mm = 1, n_basis_orb ! electron 1 
        m = list_basis(mm)
        mo_b_r1 = mos_in_r_array_omp(m,ipoint)
-       accu_1 += mo_i_r1 * mo_b_r1 * cholesky_mo(m,i,p)
+       if(dabs(mo_i_r1*mo_b_r1).lt.thresh_2)cycle
+       do p = 1, cholesky_mo_num
+        delta_vec(p) = delta_vec(p) + mo_i_r1 * mo_b_r1 * cholesky_mo_transp(p,m,i)
+       enddo
       enddo
      enddo
-     f_hf_cholesky_sparse(ipoint) += accu_1 * accu_2
-    enddo
+     do p = 1, cholesky_mo_num
+      f_hf_cholesky_sparse(ipoint) = f_hf_cholesky_sparse(ipoint) + accu_vec(p) * accu_vec(p) + accu_vec(p) * delta_vec(p)
+     enddo
     f_hf_cholesky_sparse(ipoint) *= 2.D0
    enddo
-  !$OMP END PARALLEL DO
+  !$OMP END DO
+  deallocate(accu_vec)
+  !$OMP END PARALLEL
   call wall_time(wall1)
   print*,'Time to provide f_hf_cholesky_sparse = ',wall1-wall0
  endif