diff --git a/plugins/local/spher_harm/spher_harm.irp.f b/plugins/local/spher_harm/spher_harm.irp.f
index 7a2eea06..e8deafb9 100644
--- a/plugins/local/spher_harm/spher_harm.irp.f
+++ b/plugins/local/spher_harm/spher_harm.irp.f
@@ -1,7 +1,7 @@
 program spher_harm
  implicit none
-! call test_spher_harm
+ call test_spher_harm
 ! call test_cart
- call test_brutal_spheric
+! call test_brutal_spheric
 end
 
diff --git a/plugins/local/spher_harm/spher_harm_func.irp.f b/plugins/local/spher_harm/spher_harm_func.irp.f
index 825bd8ac..f12c8fb9 100644
--- a/plugins/local/spher_harm/spher_harm_func.irp.f
+++ b/plugins/local/spher_harm/spher_harm_func.irp.f
@@ -7,6 +7,7 @@ subroutine spher_harm_func_r3(r,l,m,re_ylm, im_ylm)
  double precision :: theta, phi,r_abs
  call cartesian_to_spherical(r,theta,phi,r_abs)
  call spher_harm_func(l,m,theta,phi,re_ylm, im_ylm)
+! call spher_harm_func_expl(l,m,theta,phi,re_ylm, im_ylm)
 end
 
 
@@ -131,6 +132,10 @@ subroutine spher_harm_func_expl(l,m,theta,phi,re_ylm, im_ylm)
   tmp = - inv_sq_pi * dsqrt(3.d0/8.d0) * dsin(theta) 
   re_ylm = tmp * dcos(phi)
   im_ylm = tmp * dsin(phi)
+ else if (l==1.and.m==-1)then
+  tmp = - inv_sq_pi * dsqrt(3.d0/8.d0) * dsin(theta) 
+  re_ylm = tmp * dcos(phi)
+  im_ylm = -tmp * dsin(phi)
  else if(l==1.and.m==0)then
   tmp = inv_sq_pi * dsqrt(3.d0/4.d0) * dcos(theta) 
   re_ylm = tmp 
@@ -139,10 +144,18 @@ subroutine spher_harm_func_expl(l,m,theta,phi,re_ylm, im_ylm)
   tmp = 0.25d0 * inv_sq_pi * dsqrt(0.5d0*15.d0) * dsin(theta)*dsin(theta)
   re_ylm = tmp * dcos(2.d0*phi)
   im_ylm = tmp * dsin(2.d0*phi)
+ else if(l==2.and.m==-2)then
+  tmp = 0.25d0 * inv_sq_pi * dsqrt(0.5d0*15.d0) * dsin(theta)*dsin(theta)
+  re_ylm = tmp * dcos(2.d0*phi)
+  im_ylm =-tmp * dsin(2.d0*phi)
  else if(l==2.and.m==1)then
   tmp = - inv_sq_pi * dsqrt(15.d0/8.d0) * dsin(theta) * dcos(theta)
   re_ylm = tmp * dcos(phi)
   im_ylm = tmp * dsin(phi)
+ else if(l==2.and.m==-1)then
+  tmp = - inv_sq_pi * dsqrt(15.d0/8.d0) * dsin(theta) * dcos(theta)
+  re_ylm = tmp * dcos(phi)
+  im_ylm =-tmp * dsin(phi)
  else if(l==2.and.m==0)then
   tmp = dsqrt(5.d0/4.d0) * inv_sq_pi* (1.5d0*dcos(theta)*dcos(theta)-0.5d0)
   re_ylm = tmp
diff --git a/src/casscf_cipsi/bielec.irp.f b/src/casscf_cipsi/bielec.irp.f
index 0a44f994..a4901985 100644
--- a/src/casscf_cipsi/bielec.irp.f
+++ b/src/casscf_cipsi/bielec.irp.f
@@ -1,18 +1,25 @@
-BEGIN_PROVIDER [real*8, bielec_PQxx, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)]
+BEGIN_PROVIDER [real*8, bielec_PQxx_array, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)]
   BEGIN_DOC
-  ! bielec_PQxx : integral (pq|xx) with p,q arbitrary, x core or active
+  ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!! 
+  ! 
+  ! Replaced by the Cholesky-based function bielec_PQxx
+  !
+  ! bielec_PQxx_array : integral (pq|xx) with p,q arbitrary, x core or active
   ! indices are unshifted orbital numbers
   END_DOC
   implicit none
   integer                        :: i,j,ii,jj,p,q,i3,j3,t3,v3
   real*8                         :: mo_two_e_integral
+  print*,''
+  print*,'Providing bielec_PQxx_array, WARNING IT CAN BE A VERY BIG ARRAY WHEN MO_NUM IS LARGE !!!'
+  print*,''
   
-  bielec_PQxx(:,:,:,:) = 0.d0
+  bielec_PQxx_array(:,:,:,:) = 0.d0
   PROVIDE mo_two_e_integrals_in_map
   
   !$OMP PARALLEL DEFAULT(NONE) &
   !$OMP PRIVATE(i,ii,j,jj,i3,j3) &
-  !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PQxx, &
+  !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PQxx_array, &
   !$OMP  n_act_orb,mo_integrals_map,list_act)
 
   !$OMP DO
@@ -20,14 +27,14 @@ BEGIN_PROVIDER [real*8, bielec_PQxx, (mo_num, mo_num,n_core_inact_act_orb,n_core
     ii=list_core_inact(i)
     do j=i,n_core_inact_orb
       jj=list_core_inact(j)
-      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i,j),mo_integrals_map)
-      bielec_PQxx(:,:,j,i)=bielec_PQxx(:,:,i,j)
+      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx_array(1,1,i,j),mo_integrals_map)
+      bielec_PQxx_array(:,:,j,i)=bielec_PQxx_array(:,:,i,j)
     end do
     do j=1,n_act_orb
       jj=list_act(j)
       j3=j+n_core_inact_orb
-      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i,j3),mo_integrals_map)
-      bielec_PQxx(:,:,j3,i)=bielec_PQxx(:,:,i,j3)
+      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx_array(1,1,i,j3),mo_integrals_map)
+      bielec_PQxx_array(:,:,j3,i)=bielec_PQxx_array(:,:,i,j3)
     end do
   end do
   !$OMP END DO
@@ -40,8 +47,8 @@ BEGIN_PROVIDER [real*8, bielec_PQxx, (mo_num, mo_num,n_core_inact_act_orb,n_core
     do j=i,n_act_orb
       jj=list_act(j)
       j3=j+n_core_inact_orb
-      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i3,j3),mo_integrals_map)
-      bielec_PQxx(:,:,j3,i3)=bielec_PQxx(:,:,i3,j3)
+      call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx_array(1,1,i3,j3),mo_integrals_map)
+      bielec_PQxx_array(:,:,j3,i3)=bielec_PQxx_array(:,:,i3,j3)
     end do
   end do
   !$OMP END DO
@@ -52,9 +59,13 @@ END_PROVIDER
 
 
 
-BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)]
+BEGIN_PROVIDER [real*8, bielec_PxxQ_array, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)]
   BEGIN_DOC
-  ! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active
+  ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!! 
+  ! 
+  ! Replaced by the Cholesky-based function bielec_PxxQ
+  !
+  ! bielec_PxxQ_array : integral (px|xq) with p,q arbitrary, x core or active
   ! indices are unshifted orbital numbers
   END_DOC
   implicit none
@@ -62,12 +73,15 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_act_orb,n_core_inact_a
   double precision, allocatable  :: integrals_array(:,:)
   real*8                         :: mo_two_e_integral
   
+  print*,''
+  print*,'Providing bielec_PxxQ_array, WARNING IT CAN BE A VERY BIG ARRAY WHEN MO_NUM IS LARGE !!!'
+  print*,''
   PROVIDE mo_two_e_integrals_in_map
-  bielec_PxxQ = 0.d0
+  bielec_PxxQ_array = 0.d0
 
   !$OMP PARALLEL DEFAULT(NONE) &
   !$OMP PRIVATE(i,ii,j,jj,i3,j3,integrals_array) &
-  !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PxxQ, &
+  !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PxxQ_array, &
   !$OMP  n_act_orb,mo_integrals_map,list_act)
 
   allocate(integrals_array(mo_num,mo_num))
@@ -80,8 +94,8 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_act_orb,n_core_inact_a
       call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map)
       do q=1,mo_num
         do p=1,mo_num
-          bielec_PxxQ(p,i,j,q)=integrals_array(p,q)
-          bielec_PxxQ(p,j,i,q)=integrals_array(q,p)
+          bielec_PxxQ_array(p,i,j,q)=integrals_array(p,q)
+          bielec_PxxQ_array(p,j,i,q)=integrals_array(q,p)
         end do
       end do
     end do
@@ -91,8 +105,8 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_act_orb,n_core_inact_a
       call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map)
       do q=1,mo_num
         do p=1,mo_num
-          bielec_PxxQ(p,i,j3,q)=integrals_array(p,q)
-          bielec_PxxQ(p,j3,i,q)=integrals_array(q,p)
+          bielec_PxxQ_array(p,i,j3,q)=integrals_array(p,q)
+          bielec_PxxQ_array(p,j3,i,q)=integrals_array(q,p)
         end do
       end do
     end do
@@ -111,8 +125,8 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_act_orb,n_core_inact_a
       call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map)
       do q=1,mo_num
         do p=1,mo_num
-          bielec_PxxQ(p,i3,j3,q)=integrals_array(p,q)
-          bielec_PxxQ(p,j3,i3,q)=integrals_array(q,p)
+          bielec_PxxQ_array(p,i3,j3,q)=integrals_array(p,q)
+          bielec_PxxQ_array(p,j3,i3,q)=integrals_array(q,p)
         end do
       end do
     end do
@@ -129,10 +143,15 @@ BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
   BEGIN_DOC
   ! bielecCI : integrals (tu|vp) with p arbitrary, tuv active
   ! index p runs over the whole basis, t,u,v only over the active orbitals
+  ! 
+  ! This array can be stored anyway. Ex: 50 active orbitals, 1500 MOs ==> 8x50^3x1500 = 1.5 Gb
   END_DOC
   implicit none
   integer                        :: i,j,k,p,t,u,v
   double precision, external     :: mo_two_e_integral
+  double precision :: wall0, wall1 
+  call wall_time(wall0)
+  print*,'Providing bielecCI'
   PROVIDE mo_two_e_integrals_in_map
   
   !$OMP PARALLEL DO DEFAULT(NONE) &
@@ -151,5 +170,7 @@ BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
     end do
   end do
   !$OMP END PARALLEL DO
+  call wall_time(wall1)
+  print*,'Time to provide bielecCI = ',wall1 - wall0
 
 END_PROVIDER
diff --git a/src/casscf_cipsi/bielec_natorb.irp.f b/src/casscf_cipsi/bielec_natorb.irp.f
index 9968530c..99734a0b 100644
--- a/src/casscf_cipsi/bielec_natorb.irp.f
+++ b/src/casscf_cipsi/bielec_natorb.irp.f
@@ -1,30 +1,38 @@
- BEGIN_PROVIDER [real*8, bielec_PQxx_no, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)]
+ BEGIN_PROVIDER [real*8, bielec_PQxx_no_array, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)]
   BEGIN_DOC
+  ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!! 
+  ! 
+  ! Replaced by the Cholesky-based function bielec_PQxx_no
+  !
   ! integral (pq|xx) in the basis of natural MOs
   ! indices are unshifted orbital numbers
+  ! 
   END_DOC
   implicit none
   integer                        :: i,j,k,l,t,u,p,q
   double precision, allocatable  :: f(:,:,:), d(:,:,:)
 
+  print*,''
+  print*,'Providing bielec_PQxx_no_array, WARNING IT CAN BE A VERY BIG ARRAY WHEN MO_NUM IS LARGE !!!'
+  print*,''
 
 
   !$OMP PARALLEL DEFAULT(NONE) &
   !$OMP PRIVATE(j,k,l,p,d,f) &
   !$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, &
-  !$OMP   bielec_PQxx_no,bielec_PQxx,list_act,natorbsCI)
+  !$OMP   bielec_PQxx_no_array,bielec_PQxx_array,list_act,natorbsCI)
 
   allocate (f(n_act_orb,mo_num,n_core_inact_act_orb), &
       d(n_act_orb,mo_num,n_core_inact_act_orb))
 
   !$OMP DO
   do l=1,n_core_inact_act_orb
-    bielec_PQxx_no(:,:,:,l) = bielec_PQxx(:,:,:,l)
+    bielec_PQxx_no_array(:,:,:,l) = bielec_PQxx_array(:,:,:,l)
 
     do k=1,n_core_inact_act_orb
       do j=1,mo_num
         do p=1,n_act_orb
-          f(p,j,k)=bielec_PQxx_no(list_act(p),j,k,l)
+          f(p,j,k)=bielec_PQxx_no_array(list_act(p),j,k,l)
         end do
       end do
     end do
@@ -36,13 +44,13 @@
     do k=1,n_core_inact_act_orb
       do j=1,mo_num
         do p=1,n_act_orb
-          bielec_PQxx_no(list_act(p),j,k,l)=d(p,j,k)
+          bielec_PQxx_no_array(list_act(p),j,k,l)=d(p,j,k)
         end do
       end do
 
       do j=1,mo_num
         do p=1,n_act_orb
-          f(p,j,k)=bielec_PQxx_no(j,list_act(p),k,l)
+          f(p,j,k)=bielec_PQxx_no_array(j,list_act(p),k,l)
         end do
       end do
     end do
@@ -54,7 +62,7 @@
     do k=1,n_core_inact_act_orb
       do p=1,n_act_orb
         do j=1,mo_num
-          bielec_PQxx_no(j,list_act(p),k,l)=d(p,j,k)
+          bielec_PQxx_no_array(j,list_act(p),k,l)=d(p,j,k)
         end do
       end do
     end do
@@ -71,7 +79,7 @@
     do p=1,n_act_orb
       do k=1,mo_num
         do j=1,mo_num
-          f(j,k,p) = bielec_PQxx_no(j,k,n_core_inact_orb+p,l)
+          f(j,k,p) = bielec_PQxx_no_array(j,k,n_core_inact_orb+p,l)
         end do
       end do
     end do
@@ -83,7 +91,7 @@
     do p=1,n_act_orb
       do k=1,mo_num
         do j=1,mo_num
-          bielec_PQxx_no(j,k,n_core_inact_orb+p,l)=d(j,k,p)
+          bielec_PQxx_no_array(j,k,n_core_inact_orb+p,l)=d(j,k,p)
         end do
       end do
     end do
@@ -97,7 +105,7 @@
     do p=1,n_act_orb
       do k=1,mo_num
         do j=1,mo_num
-          f(j,k,p) = bielec_PQxx_no(j,k,l,n_core_inact_orb+p)
+          f(j,k,p) = bielec_PQxx_no_array(j,k,l,n_core_inact_orb+p)
         end do
       end do
     end do
@@ -109,7 +117,7 @@
     do p=1,n_act_orb
       do k=1,mo_num
         do j=1,mo_num
-          bielec_PQxx_no(j,k,l,n_core_inact_orb+p)=d(j,k,p)
+          bielec_PQxx_no_array(j,k,l,n_core_inact_orb+p)=d(j,k,p)
         end do
       end do
     end do
@@ -123,8 +131,12 @@ END_PROVIDER
 
 
 
-BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)]
+BEGIN_PROVIDER [real*8, bielec_PxxQ_no_array, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)]
   BEGIN_DOC
+  ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!! 
+  ! 
+  ! Replaced by the Cholesky-based function bielec_PxxQ_no
+  !
   ! integral (px|xq) in the basis of natural MOs
   ! indices are unshifted orbital numbers
   END_DOC
@@ -132,10 +144,14 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
   integer                        :: i,j,k,l,t,u,p,q
   double precision, allocatable  :: f(:,:,:), d(:,:,:)
 
+  print*,''
+  print*,'Providing bielec_PxxQ_no_array, WARNING IT CAN BE A VERY BIG ARRAY WHEN MO_NUM IS LARGE !!!'
+  print*,''
+
   !$OMP PARALLEL DEFAULT(NONE) &
   !$OMP PRIVATE(j,k,l,p,d,f) &
   !$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, &
-  !$OMP   bielec_PxxQ_no,bielec_PxxQ,list_act,natorbsCI)
+  !$OMP   bielec_PxxQ_no_array,bielec_PxxQ_array,list_act,natorbsCI)
 
 
   allocate (f(n_act_orb,n_core_inact_act_orb,n_core_inact_act_orb), &
@@ -143,11 +159,11 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
 
   !$OMP DO
   do j=1,mo_num
-    bielec_PxxQ_no(:,:,:,j) = bielec_PxxQ(:,:,:,j)
+    bielec_PxxQ_no_array(:,:,:,j) = bielec_PxxQ_array(:,:,:,j)
     do l=1,n_core_inact_act_orb
       do k=1,n_core_inact_act_orb
         do p=1,n_act_orb
-            f(p,k,l) = bielec_PxxQ_no(list_act(p),k,l,j)
+            f(p,k,l) = bielec_PxxQ_no_array(list_act(p),k,l,j)
         end do
       end do
     end do
@@ -159,7 +175,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
     do l=1,n_core_inact_act_orb
       do k=1,n_core_inact_act_orb
         do p=1,n_act_orb
-          bielec_PxxQ_no(list_act(p),k,l,j)=d(p,k,l)
+          bielec_PxxQ_no_array(list_act(p),k,l,j)=d(p,k,l)
         end do
       end do
     end do
@@ -176,7 +192,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
     do l=1,n_core_inact_act_orb
       do j=1,mo_num
         do p=1,n_act_orb
-          f(p,j,l) = bielec_PxxQ_no(j,n_core_inact_orb+p,l,k)
+          f(p,j,l) = bielec_PxxQ_no_array(j,n_core_inact_orb+p,l,k)
         end do
       end do
     end do
@@ -188,7 +204,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
     do l=1,n_core_inact_act_orb
       do j=1,mo_num
         do p=1,n_act_orb
-          bielec_PxxQ_no(j,n_core_inact_orb+p,l,k)=d(p,j,l)
+          bielec_PxxQ_no_array(j,n_core_inact_orb+p,l,k)=d(p,j,l)
         end do
       end do
     end do
@@ -205,7 +221,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
     do p=1,n_act_orb
       do l=1,n_core_inact_act_orb
         do j=1,mo_num
-          f(j,l,p) = bielec_PxxQ_no(j,l,n_core_inact_orb+p,k)
+          f(j,l,p) = bielec_PxxQ_no_array(j,l,n_core_inact_orb+p,k)
         end do
       end do
     end do
@@ -217,7 +233,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
     do p=1,n_act_orb
       do l=1,n_core_inact_act_orb
         do j=1,mo_num
-          bielec_PxxQ_no(j,l,n_core_inact_orb+p,k)=d(j,l,p)
+          bielec_PxxQ_no_array(j,l,n_core_inact_orb+p,k)=d(j,l,p)
         end do
       end do
     end do
@@ -231,7 +247,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
     do p=1,n_act_orb
       do k=1,n_core_inact_act_orb
         do j=1,mo_num
-          f(j,k,p) = bielec_PxxQ_no(j,k,l,n_core_inact_orb+p)
+          f(j,k,p) = bielec_PxxQ_no_array(j,k,l,n_core_inact_orb+p)
         end do
       end do
     end do
@@ -243,7 +259,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inac
     do p=1,n_act_orb
       do k=1,n_core_inact_act_orb
         do j=1,mo_num
-          bielec_PxxQ_no(j,k,l,n_core_inact_orb+p)=d(j,k,p)
+          bielec_PxxQ_no_array(j,k,l,n_core_inact_orb+p)=d(j,k,p)
         end do
       end do
     end do
@@ -259,10 +275,16 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
   BEGIN_DOC
   ! integrals (tu|vp) in the basis of natural MOs
   ! index p runs over the whole basis, t,u,v only over the active orbitals
+  ! 
+  ! This array can be stored anyway. Ex: 50 active orbitals, 1500 MOs ==> 8x50^3x1500 = 1.5 Gb
   END_DOC
   implicit none
   integer                        :: i,j,k,l,t,u,p,q
   double precision, allocatable  :: f(:,:,:), d(:,:,:)
+ 
+  double precision :: wall0, wall1
+  call wall_time(wall0)
+  print*,'Providing bielecCI_no'
   
   !$OMP PARALLEL DEFAULT(NONE) &
   !$OMP PRIVATE(j,k,l,p,d,f) &
@@ -363,6 +385,8 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
   
   deallocate(d,f)
   !$OMP END PARALLEL
+  call wall_time(wall1)
+  print*,'Time to provide bielecCI_no = ',wall1-wall0
 
 
 END_PROVIDER
diff --git a/src/casscf_cipsi/casscf.irp.f b/src/casscf_cipsi/casscf.irp.f
index d0a26d36..dc3e2245 100644
--- a/src/casscf_cipsi/casscf.irp.f
+++ b/src/casscf_cipsi/casscf.irp.f
@@ -11,7 +11,7 @@ program casscf
   if(small_active_space)then
    pt2_relative_error = 0.00001
   else
-   thresh_scf = 1.d-4
+   thresh_scf = max(1.d-4,thresh_scf)
    pt2_relative_error = 0.04
   endif
   touch pt2_relative_error 
diff --git a/src/casscf_cipsi/chol_bielec.irp.f b/src/casscf_cipsi/chol_bielec.irp.f
new file mode 100644
index 00000000..f69832c1
--- /dev/null
+++ b/src/casscf_cipsi/chol_bielec.irp.f
@@ -0,0 +1,248 @@
+
+BEGIN_PROVIDER [double precision, cholesky_no_1_idx_transp, (cholesky_mo_num, n_act_orb, mo_num)]
+ BEGIN_DOC
+ ! Cholesky vectors with ONE orbital on the active natural orbital basis 
+ END_DOC
+ implicit none
+ integer :: i_chol,i_act,i_mo,jj_act
+ double precision, allocatable :: chol_tmp(:,:)
+ double precision :: wall0,wall1
+ call wall_time(wall0)
+ print*,'Providing cholesky_no_1_idx_transp'
+ allocate(chol_tmp(cholesky_mo_num,n_act_orb))
+ cholesky_no_1_idx_transp = 0.D0
+ do i_mo = 1, mo_num
+  ! Get all the integrals corresponding to the "i_mo"
+  do i_act = 1, n_act_orb
+   jj_act = list_act(i_act)
+   do i_chol = 1, cholesky_mo_num
+    chol_tmp(i_chol, i_act) = cholesky_mo_transp(i_chol, jj_act, i_mo)
+   enddo
+  enddo
+  call dgemm('N','N',cholesky_mo_num,n_act_orb,n_act_orb,1.d0,  &
+        chol_tmp, size(chol_tmp,1),                              &
+        natorbsCI, size(natorbsCI,1),                                              &
+        0.d0,                                                      &
+        cholesky_no_1_idx_transp(1,1,i_mo), size(cholesky_no_1_idx_transp,1))
+ enddo
+ call wall_time(wall1)
+ print*,'Time to provide cholesky_no_1_idx_transp = ', wall1 - wall0
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [double precision, cholesky_no_2_idx_transp, (cholesky_mo_num, n_act_orb, n_act_orb)]
+ BEGIN_DOC
+ ! Cholesky vectors with TWO orbital on the active natural orbital basis 
+ END_DOC
+ implicit none
+ integer :: i_chol,i_act,j_act,jj_act
+ double precision, allocatable :: chol_tmp(:,:),chol_tmp_bis(:,:)
+ allocate(chol_tmp(cholesky_mo_num,n_act_orb),chol_tmp_bis(cholesky_mo_num,n_act_orb))
+ double precision :: wall0,wall1
+ call wall_time(wall0)
+ print*,'Providing cholesky_no_2_idx_transp'
+ cholesky_no_2_idx_transp = 0.D0
+ do i_act = 1, n_act_orb
+  ! Get all the integrals corresponding to the "j_act"
+  do j_act = 1, n_act_orb
+   jj_act = list_act(j_act)
+   do i_chol = 1, cholesky_mo_num
+    chol_tmp(i_chol, j_act) = cholesky_no_1_idx_transp(i_chol, i_act, jj_act)
+   enddo
+  enddo
+  call dgemm('N','N',cholesky_mo_num,n_act_orb,n_act_orb,1.d0,  &
+        chol_tmp, size(chol_tmp,1),                              &
+        natorbsCI, size(natorbsCI,1),                                              &
+        0.d0,                                                      &
+        cholesky_no_2_idx_transp(1,1,i_act), size(cholesky_no_2_idx_transp,1))
+ enddo
+ call wall_time(wall1)
+ print*,'Time to provide  cholesky_no_2_idx_transp = ', wall1 - wall0
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, cholesky_no_total_transp, (cholesky_mo_num, mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+ ! Cholesky vectors defined on all basis including the NO basis
+ END_DOC
+ integer :: i_chol, i_act, ii_act, j_act, jj_act, i_core_inact, j_core_inact, ii_core_inact, jj_core_inact
+ integer :: i_virt, ii_virt, j_virt, jj_virt
+ double precision :: wall0,wall1
+ call wall_time(wall0)
+ print*,'Providing cholesky_no_total_transp '
+ ! Block when two orbitals belong to the core/inact 
+ do j_core_inact = 1, n_core_inact_orb
+  jj_core_inact = list_core_inact(j_core_inact)
+  do i_core_inact = 1, n_core_inact_orb
+   ii_core_inact = list_core_inact(i_core_inact)
+   do i_chol = 1, cholesky_mo_num
+    cholesky_no_total_transp(i_chol, ii_core_inact, jj_core_inact) = cholesky_mo_transp(i_chol,ii_core_inact,jj_core_inact)
+   enddo
+  enddo
+ enddo
+
+ ! Block when one orbitals belongs to the core/inact and one belongs to the active
+ do j_core_inact = 1, n_core_inact_orb
+  jj_core_inact = list_core_inact(j_core_inact)
+  do i_act = 1, n_act_orb
+   ii_act = list_act(i_act)
+   do i_chol = 1, cholesky_mo_num 
+    cholesky_no_total_transp(i_chol,ii_act,j_core_inact) = cholesky_no_1_idx_transp(i_chol,i_act,jj_core_inact)
+   enddo
+  enddo
+ enddo
+ do j_core_inact = 1, n_core_inact_orb
+  jj_core_inact = list_core_inact(j_core_inact)
+  do i_act = 1, n_act_orb
+   ii_act = list_act(i_act)
+   do i_chol = 1, cholesky_mo_num 
+    cholesky_no_total_transp(i_chol,j_core_inact,ii_act) = cholesky_no_1_idx_transp(i_chol,i_act,jj_core_inact)
+   enddo
+  enddo
+ enddo
+ 
+ ! Block when two orbitals belong to the active 
+ do j_act = 1, n_act_orb
+  jj_act = list_act(j_act)
+  do i_act = 1, n_act_orb
+   ii_act = list_act(i_act)
+   do i_chol = 1, cholesky_mo_num 
+    cholesky_no_total_transp(i_chol,ii_act,jj_act) = cholesky_no_2_idx_transp(i_chol,i_act,j_act)
+   enddo
+  enddo
+ enddo
+
+ ! Block when two orbitals belong to the virtuals 
+ do i_virt = 1, n_virt_orb
+  ii_virt = list_virt(i_virt)
+  do j_virt = 1, n_virt_orb
+   jj_virt = list_virt(j_virt)
+   do i_chol = 1, cholesky_mo_num
+    cholesky_no_total_transp(i_chol,jj_virt,ii_virt) = cholesky_mo_transp(i_chol,jj_virt,ii_virt)
+   enddo
+  enddo
+ enddo
+ 
+ ! Block when one orbital is in active and the other in the virtuals 
+ do i_virt = 1, n_virt_orb
+  ii_virt = list_virt(i_virt)
+  do i_act = 1, n_act_orb
+   ii_act = list_act(i_act)
+   do i_chol = 1, cholesky_mo_num
+    cholesky_no_total_transp(i_chol,ii_act,ii_virt) = cholesky_no_1_idx_transp(i_chol, i_act,ii_virt)
+   enddo
+  enddo
+ enddo
+ do i_virt = 1, n_virt_orb
+  ii_virt = list_virt(i_virt)
+  do i_act = 1, n_act_orb
+   ii_act = list_act(i_act)
+   do i_chol = 1, cholesky_mo_num
+    cholesky_no_total_transp(i_chol,ii_virt,ii_act) = cholesky_no_1_idx_transp(i_chol, i_act,ii_virt)
+   enddo
+  enddo
+ enddo
+ ! Block when one orbital is in the virtual and one in the core-inact 
+ do i_virt = 1, n_virt_orb
+  ii_virt = list_virt(i_virt)
+  do i_core_inact = 1, n_core_inact_orb
+   ii_core_inact = list_core_inact(i_core_inact)
+   do i_chol = 1, cholesky_mo_num
+    cholesky_no_total_transp(i_chol, ii_core_inact, ii_virt) = cholesky_mo_transp(i_chol, ii_core_inact, ii_virt)
+   enddo
+  enddo
+ enddo
+ do i_core_inact = 1, n_core_inact_orb
+  ii_core_inact = list_core_inact(i_core_inact)
+  do i_virt = 1, n_virt_orb
+   ii_virt = list_virt(i_virt)
+   do i_chol = 1, cholesky_mo_num
+    cholesky_no_total_transp(i_chol, ii_virt, ii_core_inact) = cholesky_mo_transp(i_chol, ii_virt, ii_core_inact)
+   enddo
+  enddo
+ enddo
+
+ call wall_time(wall1)
+ print*,'Time to provide cholesky_no_total_transp = ', wall1 - wall0
+END_PROVIDER 
+
+
+double precision function bielec_no_basis(i_1,j_1,i_2,j_2)
+ implicit none
+ integer, intent(in) :: i_1,j_1,i_2,j_2
+  BEGIN_DOC
+  ! integral (i_1 j_1|i_2 j_2) in the mixed basis of both MOs and natural MOs
+  ! 
+  END_DOC
+  integer :: i 
+ bielec_no_basis = 0.d0
+ do i = 1, cholesky_mo_num
+  bielec_no_basis += cholesky_no_total_transp(i,i_1, j_1) * cholesky_no_total_transp(i,i_2,j_2)
+ enddo
+end
+
+double precision function bielec_PQxx_no(i_mo, j_mo, i_ca, j_ca)
+ implicit none
+ BEGIN_DOC
+ ! function that computes (i_mo j_mo| i_ca j_ca) with Cholesky decomposition  on the NO basis for active orbitals 
+ ! 
+ ! where i_ca, j_ca are in [1:n_core_inact_act_orb]
+ END_DOC
+ integer, intent(in) :: i_ca, j_ca, i_mo, j_mo
+ integer :: ii_ca, jj_ca
+ double precision :: bielec_no_basis
+ ii_ca = list_core_inact_act(i_ca)
+ jj_ca = list_core_inact_act(j_ca)
+ bielec_PQxx_no = bielec_no_basis(i_mo,j_mo,ii_ca,jj_ca)
+end
+
+double precision function bielec_PxxQ_no(i_mo, j_ca, i_ca, j_mo)
+ implicit none 
+  BEGIN_DOC
+  ! function that computes (i_mo j_ca |i_ca j_mo) with Cholesky decomposition  on the NO basis for active orbitals 
+  ! 
+  ! where i_ca, j_ca are in [1:n_core_inact_act_orb]
+  END_DOC
+ integer, intent(in) :: i_ca, j_ca, i_mo, j_mo
+ integer :: ii_ca, jj_ca
+ double precision :: bielec_no_basis
+ ii_ca = list_core_inact_act(i_ca)
+ jj_ca = list_core_inact_act(j_ca)
+ bielec_PxxQ_no = bielec_no_basis(i_mo, jj_ca, ii_ca, j_mo)
+
+end
+
+
+double precision function bielec_PQxx(i_mo, j_mo, i_ca, j_ca)
+  BEGIN_DOC
+  ! function that computes (i_mo j_mo |i_ca j_ca) with Cholesky decomposition 
+  ! 
+  ! indices are unshifted orbital numbers
+  ! 
+  ! where i_ca, j_ca are in [1:n_core_inact_act_orb]
+  END_DOC
+ implicit none 
+ integer, intent(in) :: i_ca, j_ca, j_mo, i_mo
+ double precision :: mo_two_e_integral
+ integer :: ii_ca, jj_ca
+ ii_ca = list_core_inact_act(i_ca)
+ jj_ca = list_core_inact_act(j_ca)
+ bielec_PQxx = mo_two_e_integral(i_mo,ii_ca,j_mo,jj_ca)
+end
+
+double precision function bielec_PxxQ(i_mo, i_ca, j_ca, j_mo)
+  BEGIN_DOC
+  ! function that computes (i_mo j_mo |i_ca j_ca) with Cholesky decomposition 
+  ! 
+  ! where i_ca, j_ca are in [1:n_core_inact_act_orb]
+  END_DOC
+ implicit none
+ integer, intent(in) :: i_ca, j_ca, j_mo, i_mo
+ double precision :: mo_two_e_integral
+ integer :: ii_ca, jj_ca
+ ii_ca = list_core_inact_act(i_ca)
+ jj_ca = list_core_inact_act(j_ca)
+ bielec_PxxQ = mo_two_e_integral(i_mo,jj_ca,ii_ca,j_mo)
+end
+
diff --git a/src/casscf_cipsi/chol_garb.irp.f b/src/casscf_cipsi/chol_garb.irp.f
new file mode 100644
index 00000000..c4a8fa59
--- /dev/null
+++ b/src/casscf_cipsi/chol_garb.irp.f
@@ -0,0 +1,34 @@
+
+!!!!! FUNCTIONS THAT WORK BUT WHICH ARE USELESS AS THE ARRAYS CAN ALWAYS BE STORED
+!double precision function bielecCI_chol(i_a, j_a, k_a, i_mo)
+!  BEGIN_DOC
+!  ! function that computes (i_a j_a |k_a j_mo) with Cholesky decomposition 
+!  ! 
+!  ! where i_a, j_a, k_a are in [1:n_act_orb] !!! ONLY ON ACTIVE 
+!  END_DOC
+! implicit none
+! integer, intent(in) :: i_a, j_a, k_a, i_mo
+! integer :: ii_a, jj_a, kk_a
+! double precision :: mo_two_e_integral
+! ii_a = list_act(i_a)
+! jj_a = list_act(j_a)
+! kk_a = list_act(k_a)
+! bielecCI_chol = mo_two_e_integral(ii_a,kk_a,jj_a,i_mo)
+!end
+
+!double precision function bielecCI_no_chol(i_ca, j_ca, k_ca, i_mo)
+!  BEGIN_DOC
+!  ! function that computes (i_ca j_ca |k_ca j_mo) with Cholesky decomposition on the NO basis for active orbitals 
+!  ! 
+!  ! where i_ca, j_ca, k_ca are in [1:n_core_inact_act_orb]
+!  END_DOC
+! implicit none 
+! integer, intent(in) :: i_ca, j_ca, k_ca, i_mo
+! integer :: ii_ca, jj_ca, kk_ca
+! double precision :: bielec_no_basis_chol
+! ii_ca = list_act(i_ca)
+! jj_ca = list_act(j_ca)
+! kk_ca = list_act(k_ca)
+! bielecCI_no_chol = bielec_no_basis_chol(ii_ca, jj_ca, kk_ca, i_mo)
+!
+!end
diff --git a/src/casscf_cipsi/gradient.irp.f b/src/casscf_cipsi/gradient.irp.f
index a1c5e947..961d260d 100644
--- a/src/casscf_cipsi/gradient.irp.f
+++ b/src/casscf_cipsi/gradient.irp.f
@@ -157,6 +157,7 @@ real*8 function gradvec_it(i,t)
   
   integer                        :: ii,tt,v,vv,x,y
   integer                        :: x3,y3
+  double precision :: bielec_PQxx_no
   
   ii=list_core_inact(i)
   tt=list_act(t)
diff --git a/src/casscf_cipsi/hessian.irp.f b/src/casscf_cipsi/hessian.irp.f
index 458c6aa6..9a7a9031 100644
--- a/src/casscf_cipsi/hessian.irp.f
+++ b/src/casscf_cipsi/hessian.irp.f
@@ -10,6 +10,7 @@ real*8 function hessmat_itju(i,t,j,u)
   implicit none
   integer                        :: i,t,j,u,ii,tt,uu,v,vv,x,xx,y,jj
   real*8                         :: term,t2
+  double precision :: bielec_pqxx_no,bielec_pxxq_no
   
   ii=list_core_inact(i)
   tt=list_act(t)
@@ -95,6 +96,7 @@ real*8 function hessmat_itja(i,t,j,a)
   implicit none
   integer                        :: i,t,j,a,ii,tt,jj,aa,v,vv,x,y
   real*8                         :: term
+  double precision :: bielec_pqxx_no,bielec_pxxq_no
   
   ! it/ja
   ii=list_core_inact(i)
@@ -128,6 +130,7 @@ real*8 function hessmat_itua(i,t,u,a)
   implicit none
   integer                        :: i,t,u,a,ii,tt,uu,aa,v,vv,x,xx,u3,t3,v3
   real*8                         :: term
+  double precision :: bielec_pqxx_no,bielec_pxxq_no
   
   ii=list_core_inact(i)
   tt=list_act(t)
@@ -169,6 +172,7 @@ real*8 function hessmat_iajb(i,a,j,b)
   implicit none
   integer                        :: i,a,j,b,ii,aa,jj,bb
   real*8                         :: term
+  double precision :: bielec_pqxx_no,bielec_pxxq_no
   
   ii=list_core_inact(i)
   aa=list_virt(a)
@@ -205,6 +209,7 @@ real*8 function hessmat_iatb(i,a,t,b)
   implicit none
   integer                        :: i,a,t,b,ii,aa,tt,bb,v,vv,x,y,v3,t3
   real*8                         :: term
+  double precision :: bielec_pqxx_no,bielec_pxxq_no
   
   ii=list_core_inact(i)
   aa=list_virt(a)
@@ -237,6 +242,7 @@ real*8 function hessmat_taub(t,a,u,b)
   integer                        :: t,a,u,b,tt,aa,uu,bb,v,vv,x,xx,y
   integer                        :: v3,x3
   real*8                         :: term,t1,t2,t3
+  double precision :: bielec_pqxx_no,bielec_pxxq_no
   
   tt=list_act(t)
   aa=list_virt(a)
diff --git a/src/casscf_cipsi/mcscf_fock.irp.f b/src/casscf_cipsi/mcscf_fock.irp.f
index 0f4b7a99..82b710a7 100644
--- a/src/casscf_cipsi/mcscf_fock.irp.f
+++ b/src/casscf_cipsi/mcscf_fock.irp.f
@@ -4,6 +4,7 @@ BEGIN_PROVIDER [real*8, Fipq, (mo_num,mo_num) ]
    END_DOC
    implicit none
    integer                        :: p,q,k,kk,t,tt,u,uu
+   double precision :: bielec_pxxq_no, bielec_pqxx_no
    
    do q=1,mo_num
      do p=1,mo_num
@@ -44,6 +45,7 @@ BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ]
    END_DOC
    implicit none
    integer                        :: p,q,k,kk,t,tt,u,uu
+   double precision :: bielec_pxxq_no, bielec_pqxx_no
    
    Fapq = 0.d0
    
diff --git a/src/casscf_cipsi/test_chol.irp.f b/src/casscf_cipsi/test_chol.irp.f
new file mode 100644
index 00000000..bcce7cf7
--- /dev/null
+++ b/src/casscf_cipsi/test_chol.irp.f
@@ -0,0 +1,116 @@
+program test_chol
+ implicit none
+ read_wf= .True.
+ touch read_wf 
+! call routine_bielec_PxxQ_no
+! call routine_bielecCI_no
+! call test_bielec_PxxQ_chol
+! call test_bielecCI
+
+end
+
+subroutine routine_bielec_PQxx_no
+ implicit none
+ integer :: i_chol, i_act, ii_act, j_act, jj_act, i_core_inact, j_core_inact, ii_core_inact, jj_core_inact
+ integer :: i_virt, ii_virt, j_virt, jj_virt, i_mo, j_mo
+ double precision :: exact, new, error, accu, bielec_no_basis_chol
+ double precision :: bielec_PQxx_no
+ 
+ accu = 0.d0
+ do i_core_inact = 1, n_core_inact_act_orb
+  ii_core_inact = list_core_inact_act(i_core_inact)
+  do j_core_inact = 1, n_core_inact_act_orb
+   jj_core_inact = list_core_inact_act(j_core_inact)
+   do i_mo = 1, mo_num
+    do j_mo = 1, mo_num
+     exact = bielec_PQxx_no_array(j_mo,i_mo, j_core_inact, i_core_inact) 
+     new   = bielec_PQxx_no(j_mo,i_mo, j_core_inact, i_core_inact) 
+     error = dabs(exact-new)
+     if(dabs(exact).gt.1.d-10)then
+      print*,exact,new,error
+     endif
+     accu += error
+    enddo
+   enddo
+  enddo
+ enddo
+ print*,'accu = ',accu/(dble(mo_num*mo_num*n_core_inact_act_orb**2))
+end
+
+subroutine routine_bielec_PxxQ_no_array
+ implicit none
+ integer :: i_chol, i_act, ii_act, j_act, jj_act, i_core_inact, j_core_inact, ii_core_inact, jj_core_inact
+ integer :: i_virt, ii_virt, j_virt, jj_virt, i_mo, j_mo
+ double precision :: exact, new, error, accu, bielec_no_basis_chol
+ double precision :: bielec_PxxQ_no
+ 
+ accu = 0.d0
+ do i_mo = 1, mo_num
+  do i_core_inact = 1, n_core_inact_act_orb
+  ii_core_inact = list_core_inact_act(i_core_inact)
+   do j_core_inact = 1, n_core_inact_act_orb
+   jj_core_inact = list_core_inact_act(j_core_inact)
+    do j_mo = 1, mo_num
+     exact = bielec_PxxQ_no_array(j_mo, j_core_inact,  i_core_inact,i_mo) 
+!     new   = bielec_no_basis_chol(j_mo,i_mo, jj_core_inact, ii_core_inact) 
+     new   = bielec_PxxQ_no(j_mo, j_core_inact,  i_core_inact,i_mo) 
+     error = dabs(exact-new)
+     accu += error
+     if(dabs(exact).gt.1.d-10)then
+      print*,exact,new,error
+     endif
+    enddo
+   enddo
+  enddo
+ enddo
+ print*,'accu = ',accu/(dble(mo_num*mo_num*n_core_inact_act_orb**2))
+end
+
+subroutine test_bielec_PQxx(i_mo, j_mo, i_ca, j_ca)
+ implicit none
+ integer :: i_mo, j_mo, i_ca, j_ca 
+ double precision :: exact, new, error, accu
+ double precision :: bielec_PQxx
+ 
+ accu = 0.d0
+ do j_ca = 1, n_core_inact_act_orb
+  do i_ca = 1, n_core_inact_act_orb
+   do j_mo = 1, mo_num
+    do i_mo = 1, mo_num
+     exact = bielec_PQxx_array(i_mo, j_mo, i_ca, j_ca)
+     new   = bielec_PQxx(i_mo, j_mo, i_ca, j_ca)
+     error = dabs(exact-new)
+     accu += error
+     if(dabs(exact).gt.1.d-10)then
+      print*,exact,new,error
+     endif
+    enddo
+   enddo
+  enddo
+ enddo
+ print*,'accu = ',accu/(dble(mo_num*mo_num*n_core_inact_act_orb**2))
+end
+
+subroutine test_bielec_PxxQ_chol(i_mo, i_ca, j_ca, j_mo)
+ implicit none
+ integer :: i_mo, i_ca, j_ca, j_mo
+ double precision :: exact, new, error, accu
+ double precision :: bielec_PxxQ
+ accu = 0.d0
+ do j_mo = 1, mo_num
+  do j_ca = 1, n_core_inact_act_orb
+   do i_ca =1, n_core_inact_act_orb
+    do i_mo = 1, mo_num
+     exact = bielec_PxxQ_array(i_mo, i_ca, j_ca, j_mo)
+     new   = bielec_PxxQ(i_mo, i_ca, j_ca, j_mo)
+     error = dabs(exact-new)
+     accu += error
+     if(dabs(exact).gt.1.d-10)then
+      print*,exact,new,error
+     endif
+    enddo
+   enddo
+  enddo
+ enddo
+ print*,'accu = ',accu/(dble(mo_num*mo_num*n_core_inact_act_orb**2))
+end
diff --git a/src/casscf_cipsi/tot_en.irp.f b/src/casscf_cipsi/tot_en.irp.f
index 1d70e087..37ceac05 100644
--- a/src/casscf_cipsi/tot_en.irp.f
+++ b/src/casscf_cipsi/tot_en.irp.f
@@ -8,6 +8,7 @@
    implicit none
    integer                        :: t,u,v,x,i,ii,tt,uu,vv,xx,j,jj,t3,u3,v3,x3
    real*8                         :: e_one_all,e_two_all
+   double precision :: bielec_PQxx,bielec_PxxQ
    e_one_all=0.D0
    e_two_all=0.D0
    do i=1,n_core_inact_orb
diff --git a/src/mo_two_e_ints/cholesky.irp.f b/src/mo_two_e_ints/cholesky.irp.f
index 7e2c8b37..1fed949d 100644
--- a/src/mo_two_e_ints/cholesky.irp.f
+++ b/src/mo_two_e_ints/cholesky.irp.f
@@ -101,3 +101,34 @@ BEGIN_PROVIDER [ double precision, cholesky_mo_transp, (cholesky_mo_num, mo_num,
 
 END_PROVIDER
 
+
+BEGIN_PROVIDER [ double precision, cholesky_semi_mo_transp_simple, (cholesky_mo_num, ao_num, mo_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Cholesky vectors in MO basis
+ END_DOC
+
+ double precision, allocatable :: X(:,:,:)
+ double precision :: wall0, wall1
+ integer :: ierr
+ print *, 'Semi 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
+   print *, irp_here, ': Allocation failed'
+ endif
+ integer :: i_chol, i_mo, j_mo, i_ao 
+ cholesky_semi_mo_transp_simple = 0.d0
+ do i_mo = 1, mo_num
+  do i_ao = 1, ao_num
+   do j_mo = 1, mo_num
+    do i_chol = 1, cholesky_mo_num
+     cholesky_semi_mo_transp_simple(i_chol, i_ao,i_mo) += cholesky_mo_transp(i_chol,j_mo,i_mo) * mo_coef_transp(j_mo,i_ao)
+    enddo
+   enddo
+  enddo
+ enddo
+
+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 5dd69eb6..179b80dd 100644
--- a/src/mu_of_r/f_hf_cholesky.irp.f
+++ b/src/mu_of_r/f_hf_cholesky.irp.f
@@ -289,6 +289,106 @@ BEGIN_PROVIDER [ double precision, f_hf_cholesky_sparse, (n_points_final_grid)]
  endif
 END_PROVIDER 
 
+BEGIN_PROVIDER [ double precision, f_hf_cholesky_sparse_bis, (n_points_final_grid)]
+ implicit none
+ integer :: ipoint,m,mm,i,ii,p
+ !!f(R) =  \sum_{I} \sum_{J} Phi_I(R) Phi_J(R) V_IJ
+ !!     =  \sum_{I}\sum_{J}\sum_A Phi_I(R) Phi_J(R) V_AI V_AJ
+ !!     =  \sum_A \sum_{I}Phi_I(R)V_AI \sum_{J}V_AJ Phi_J(R)
+ !!     =  \sum_A V_AR G_AR 
+ !! 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(:),delta_vec(:)
+ thresh_2 = ao_cholesky_threshold * 100.d0
+ thresh_1 = dsqrt(thresh_2)
+ provide cholesky_mo_transp
+ if(elec_alpha_num == elec_beta_num)then
+  call wall_time(wall0)
+  !$OMP PARALLEL DEFAULT(NONE)                                      &
+  !$OMP PRIVATE (accu_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,mos_in_r_array_omp,aos_in_r_array,thresh_1,thresh_2) & 
+  !$OMP ShARED (cholesky_mo_num,f_hf_cholesky_sparse_bis,n_points_final_grid,cholesky_semi_mo_transp_simple,ao_num) 
+  allocate(accu_vec(cholesky_mo_num))
+  !$OMP DO 
+   do ipoint = 1, n_points_final_grid
+    f_hf_cholesky_sparse_bis(ipoint) = 0.d0
+     accu_vec = 0.d0
+     do ii = 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, ao_num ! electron 1 
+       mo_b_r1 = aos_in_r_array(mm,ipoint)*mo_i_r1
+       if(dabs(mo_b_r1).lt.thresh_2)cycle
+       do p = 1, cholesky_mo_num
+        accu_vec(p) = accu_vec(p) + mo_b_r1 * cholesky_semi_mo_transp_simple(p,mm,i)
+       enddo
+      enddo
+     enddo
+     do p = 1, cholesky_mo_num
+      f_hf_cholesky_sparse_bis(ipoint) = f_hf_cholesky_sparse_bis(ipoint) + accu_vec(p) * accu_vec(p)
+     enddo
+    f_hf_cholesky_sparse_bis(ipoint) *= 2.D0
+   enddo
+  !$OMP END DO
+  deallocate(accu_vec)
+  !$OMP END PARALLEL
+  
+  call wall_time(wall1)
+  print*,'Time to provide f_hf_cholesky_sparse_bis = ',wall1-wall0
+ else
+  call wall_time(wall0)
+  !$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_bis,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_bis(ipoint) = 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)
+       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
+     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)
+       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
+     do p = 1, cholesky_mo_num
+      f_hf_cholesky_sparse_bis(ipoint) = f_hf_cholesky_sparse_bis(ipoint) + accu_vec(p) * accu_vec(p) + accu_vec(p) * delta_vec(p)
+     enddo
+    f_hf_cholesky_sparse_bis(ipoint) *= 2.D0
+   enddo
+  !$OMP END DO
+  deallocate(accu_vec)
+  !$OMP END PARALLEL
+  call wall_time(wall1)
+  print*,'Time to provide f_hf_cholesky_sparse_bis = ',wall1-wall0
+ endif
+END_PROVIDER 
+
+
 BEGIN_PROVIDER [ double precision, on_top_hf_grid, (n_points_final_grid)]
  implicit none
  integer :: ipoint,i,ii
diff --git a/src/mu_of_r/test_proj_op.irp.f b/src/mu_of_r/test_proj_op.irp.f
index f9aba094..fd5e976b 100644
--- a/src/mu_of_r/test_proj_op.irp.f
+++ b/src/mu_of_r/test_proj_op.irp.f
@@ -15,7 +15,23 @@ program projected_operators
 !  call test_f_HF_valence_ab
 !  call routine_full_mos
 !   call test_f_ii_valence_ab
-   call test_f_ia_valence_ab
-  call test_f_ii_ia_aa_valence_ab
+!   call test_f_ia_valence_ab
+!  call test_f_ii_ia_aa_valence_ab
+ call test
 end
 
+
+subroutine test
+ implicit none
+ integer :: i_point
+ double precision :: ref, new, accu, weight
+ accu = 0.d0
+ do i_point = 1, n_points_final_grid
+  ref = f_hf_cholesky_sparse(i_point)
+  new = f_hf_cholesky_sparse_bis(i_point)
+  weight = final_weight_at_r_vector(i_point)
+  accu += dabs(ref - new) * weight
+ enddo
+ print*,'accu = ',accu
+
+end