diff --git a/external/qp2-dependencies b/external/qp2-dependencies
index 242151e0..90ee61f5 160000
--- a/external/qp2-dependencies
+++ b/external/qp2-dependencies
@@ -1 +1 @@
-Subproject commit 242151e03d1d6bf042387226431d82d35845686a
+Subproject commit 90ee61f5041c7c94a0c605625a264860292813a0
diff --git a/src/ao_one_e_ints/pseudopot.f90 b/src/ao_one_e_ints/pseudopot.f90
index 09e0291a..7321dff7 100644
--- a/src/ao_one_e_ints/pseudopot.f90
+++ b/src/ao_one_e_ints/pseudopot.f90
@@ -1869,7 +1869,7 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
       qk = dble(q)
       two_qkmp1 = 2.d0*(qk+mk)+1.d0
       do k=0,q-1
-        if (s_q_k < 1.d-32) then
+        if (s_q_k < 1.d-20) then
           s_q_k = 0.d0
           exit
         endif
diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f
index 0d34d95e..0ce6b87e 100644
--- a/src/ao_two_e_ints/map_integrals.irp.f
+++ b/src/ao_two_e_ints/map_integrals.irp.f
@@ -142,7 +142,7 @@ subroutine ao_idx2_sq(i,j,ij)
     ij=i*i
   endif
 end
-  
+
 subroutine idx2_tri_int(i,j,ij)
   implicit none
   integer, intent(in)  :: i,j
@@ -152,7 +152,7 @@ subroutine idx2_tri_int(i,j,ij)
   q = min(i,j)
   ij = q+ishft(p*p-p,-1)
 end
-  
+
 subroutine ao_idx2_tri_key(i,j,ij)
   use map_module
   implicit none
@@ -163,8 +163,8 @@ subroutine ao_idx2_tri_key(i,j,ij)
   q = min(i,j)
   ij = q+ishft(p*p-p,-1)
 end
-  
-subroutine two_e_integrals_index_2fold(i,j,k,l,i1)                                           
+
+subroutine two_e_integrals_index_2fold(i,j,k,l,i1)
   use map_module
   implicit none
   integer, intent(in)            :: i,j,k,l
@@ -176,7 +176,7 @@ subroutine two_e_integrals_index_2fold(i,j,k,l,i1)
   call ao_idx2_tri_key(ik,jl,i1)
 end
 
-subroutine ao_idx2_sq_rev(i,k,ik)                                                                
+subroutine ao_idx2_sq_rev(i,k,ik)
   BEGIN_DOC
   ! reverse square compound index
   END_DOC
@@ -399,7 +399,7 @@ BEGIN_PROVIDER [ complex*16, ao_integrals_cache_periodic, (0:64*64*64*64) ]
            tmp_im = 0.d0
            integral = dcmplx(tmp_re,tmp_im)
          endif
-         
+
          ii = l-ao_integrals_cache_min
          ii = ior( shiftl(ii,6), k-ao_integrals_cache_min)
          ii = ior( shiftl(ii,6), j-ao_integrals_cache_min)
@@ -474,7 +474,7 @@ subroutine get_ao_two_e_integrals(j,k,l,sze,out_val)
   BEGIN_DOC
   ! Gets multiple AO bi-electronic integral from the AO map .
   ! All i are retrieved for j,k,l fixed.
-  ! physicist convention : <ij|kl> 
+  ! physicist convention : <ij|kl>
   END_DOC
   implicit none
   integer, intent(in)            :: j,k,l, sze
@@ -483,7 +483,7 @@ subroutine get_ao_two_e_integrals(j,k,l,sze,out_val)
   integer                        :: i
   integer(key_kind)              :: hash
   logical, external              :: ao_one_e_integral_zero
-  PROVIDE ao_two_e_integrals_in_map ao_integrals_map
+  PROVIDE ao_two_e_integrals_in_map ao_integrals_map 
 
   if (ao_one_e_integral_zero(j,l)) then
     out_val = 0.d0
@@ -503,7 +503,7 @@ subroutine get_ao_two_e_integrals_periodic(j,k,l,sze,out_val)
   BEGIN_DOC
   ! Gets multiple AO bi-electronic integral from the AO map .
   ! All i are retrieved for j,k,l fixed.
-  ! physicist convention : <ij|kl> 
+  ! physicist convention : <ij|kl>
   END_DOC
   implicit none
   integer, intent(in)            :: j,k,l, sze
diff --git a/src/mo_two_e_ints/mo_bi_integrals.irp.f b/src/mo_two_e_ints/mo_bi_integrals.irp.f
index 6f4c5c17..56d8cf28 100644
--- a/src/mo_two_e_ints/mo_bi_integrals.irp.f
+++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f
@@ -38,7 +38,7 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
     print*, 'MO integrals provided'
     return
   else
-    PROVIDE ao_two_e_integrals_in_map
+    PROVIDE ao_two_e_integrals_in_map 
   endif
 
   print *,  ''
@@ -53,7 +53,7 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
 !    call four_idx_novvvv
     call four_idx_novvvv_old
   else
-    if (32.d-9*dble(ao_num)**4 < dble(qp_max_mem)) then
+    if (dble(ao_num)**4 * 32.d-9 < dble(qp_max_mem)) then
       call four_idx_dgemm
     else
       call add_integrals_to_map(full_ijkl_bitmask_4)
@@ -245,14 +245,14 @@ subroutine add_integrals_to_map(mask_ijkl)
     return
   endif
 
-  size_buffer = min(ao_num*ao_num*ao_num,16000000)
-  print*, 'Buffers : ', 8.*(mo_num*(n_j)*(n_k+1) + mo_num+&
-      ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core'
-
   double precision               :: accu_bis
   accu_bis = 0.d0
   call wall_time(wall_1)
 
+  size_buffer = min( (qp_max_mem/(nproc*5)),mo_num*mo_num*mo_num)
+  print*, 'Buffers : ', 8.*(mo_num*(n_j)*(n_k+1) + mo_num+&
+      ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core'
+
   !$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax,   &
       !$OMP  two_e_tmp_0_idx, two_e_tmp_0, two_e_tmp_1,two_e_tmp_2,two_e_tmp_3,&
       !$OMP  buffer_i,buffer_value,n_integrals,wall_2,i0,j0,k0,l0,   &
@@ -263,6 +263,10 @@ subroutine add_integrals_to_map(mask_ijkl)
       !$OMP  mo_coef_transp_is_built, list_ijkl,                     &
       !$OMP  mo_coef_is_built, wall_1,                               &
       !$OMP  mo_coef,mo_integrals_threshold,mo_integrals_map)
+
+  thread_num = 0
+  !$  thread_num = omp_get_thread_num()
+
   n_integrals = 0
   wall_0 = wall_1
   allocate(two_e_tmp_3(mo_num, n_j, n_k),                 &
@@ -273,8 +277,6 @@ subroutine add_integrals_to_map(mask_ijkl)
       buffer_i(size_buffer),                                         &
       buffer_value(size_buffer) )
 
-  thread_num = 0
-  !$  thread_num = omp_get_thread_num()
   !$OMP DO SCHEDULE(guided)
   do l1 = 1,ao_num
     two_e_tmp_3 = 0.d0
diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f
index 61506a87..60c84a3b 100644
--- a/src/utils/linear_algebra.irp.f
+++ b/src/utils/linear_algebra.irp.f
@@ -1985,3 +1985,52 @@ end subroutine diag_nonsym_right
 ! Taken from GammCor thanks to Michal Hapka :-)
 
 
+subroutine pivoted_cholesky( A, rank, tol, ndim, U)
+!
+! A = U**T * U
+!
+! matrix A is destroyed inside this subroutine
+! Cholesky vectors are stored in U
+! dimension of U: U(1:rank, 1:n)
+! U is allocated inside this subroutine
+! rank is the number of Cholesky vectors depending on tol
+!
+integer :: ndim
+integer, intent(inout)                                        :: rank
+double precision, dimension(ndim, ndim), intent(inout)        :: A
+double precision, dimension(ndim, rank), intent(out)          :: U
+double precision, intent(in)                                  :: tol
+
+integer, dimension(:), allocatable          :: piv
+double precision, dimension(:), allocatable :: work
+character, parameter :: uplo = "U"
+integer :: N, LDA
+integer :: info
+integer :: k, l, rank0
+external :: dpstrf
+
+rank0 = rank
+N = size(A, dim=1)
+LDA = N
+allocate(piv(N))
+allocate(work(2*N))
+call dpstrf(uplo, N, A, LDA, piv, rank, tol, work, info)
+
+if (rank > rank0) then
+  print *, 'Bug: rank > rank0 in pivoted cholesky. Increase rank before calling'
+  stop
+end if
+
+do k = 1, N
+  A(k+1:, k) = 0.00D+0
+end do
+! TODO: It should be possible to use only one vector of size (1:rank) as a buffer
+! to do the swapping in-place
+U = 0.00D+0
+do k = 1, N
+  l = piv(k)
+  U(l, :) = A(1:rank, k)
+end do
+
+end subroutine pivoted_cholesky
+