diff --git a/src/ao_two_e_ints/cholesky.irp.f b/src/ao_two_e_ints/cholesky.irp.f
index 175ccf6e..2977f0f4 100644
--- a/src/ao_two_e_ints/cholesky.irp.f
+++ b/src/ao_two_e_ints/cholesky.irp.f
@@ -23,38 +23,38 @@ END_PROVIDER
    !
    ! Last dimension of cholesky_ao is cholesky_ao_num
    END_DOC
-   
+
    integer                        :: rank, ndim
    double precision               :: tau
    double precision, pointer      :: L(:,:), L_old(:,:)
-   
-   
+
+
    double precision               :: s
    double precision, parameter    :: dscale = 1.d0
-   
+
    double precision, allocatable  :: D(:), Delta(:,:), Ltmp_p(:,:), Ltmp_q(:,:)
    integer, allocatable           :: Lset(:), Dset(:), addr(:,:)
    logical, allocatable           :: computed(:)
-   
+
    integer                        :: i,j,k,m,p,q, qj, dj, p2, q2
    integer                        :: N, np, nq
-   
+
    double precision               :: Dmax, Dmin, Qmax, f
    double precision, external     :: get_ao_two_e_integral
    logical, external              :: ao_two_e_integral_zero
-   
+
    double precision, external     :: ao_two_e_integral
    integer                        :: block_size, iblock, ierr
-   
+
    double precision               :: mem
    double precision, external     :: memory_of_double, memory_of_int
-   
+
    integer, external              :: getUnitAndOpen
    integer                        :: iunit
-   
+
    ndim = ao_num*ao_num
    deallocate(cholesky_ao)
-   
+
    if (read_ao_cholesky) then
      print *,  'Reading Cholesky vectors from disk...'
      iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'R')
@@ -63,11 +63,11 @@ END_PROVIDER
      read(iunit) cholesky_ao
      close(iunit)
      cholesky_ao_num = rank
-     
+
    else
-     
+
      PROVIDE nucl_coord
-     
+
      if (do_direct_integrals) then
        if (ao_two_e_integral(1,1,1,1) < huge(1.d0)) then
          ! Trigger providers inside ao_two_e_integral
@@ -76,16 +76,16 @@ END_PROVIDER
      else
        PROVIDE ao_two_e_integrals_in_map
      endif
-     
+
      tau = ao_cholesky_threshold
-     
+
      mem = 6.d0 * memory_of_double(ndim) + 6.d0 * memory_of_int(ndim)
      call check_mem(mem, irp_here)
-     
+
      call print_memory_usage()
 
      allocate(L(ndim,1))
-     
+
      print *,  ''
      print *,  'Cholesky decomposition of AO integrals'
      print *,  '======================================'
@@ -93,13 +93,13 @@ END_PROVIDER
      print *,  '============ ============='
      print *,  '    Rank      Threshold'
      print *,  '============ ============='
-     
-     
+
+
      rank = 0
-     
+
      allocate( D(ndim), Lset(ndim), Dset(ndim) )
      allocate( addr(3,ndim) )
-     
+
      ! 1.
      k=0
      do j=1,ao_num
@@ -110,7 +110,7 @@ END_PROVIDER
          addr(3,k) = (i-1)*ao_num + j
        enddo
      enddo
-     
+
      if (do_direct_integrals) then
        !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) SCHEDULE(guided)
        do i=1,ndim
@@ -127,9 +127,9 @@ END_PROVIDER
        enddo
        !$OMP END PARALLEL DO
      endif
-     
+
      Dmax = maxval(D)
-     
+
      ! 2.
      np=0
      do p=1,ndim
@@ -138,26 +138,26 @@ END_PROVIDER
          Lset(np) = p
        endif
      enddo
-     
+
      ! 3.
      N = 0
-     
+
      ! 4.
      i = 0
-     
+
      ! 5.
      do while ( (Dmax > tau).and.(rank < ndim) )
        ! a.
        i = i+1
-       
+
        s = 0.01d0
-       
+
        ! Inrease s until the arrays fit in memory
        do while (.True.)
-         
+
          ! b.
          Dmin = max(s*Dmax,tau)
-         
+
          ! c.
          nq=0
          do p=1,np
@@ -166,30 +166,30 @@ END_PROVIDER
              Dset(nq) = Lset(p)
            endif
          enddo
-         
+
          call total_memory(mem)
          mem = mem                                                   &
              + np*memory_of_double(nq)                               &! Delta(np,nq)
              + (rank+nq)* memory_of_double(ndim)                     &! L(ndim,rank+nq)
              + (np+nq)*memory_of_double(block_size)    ! Ltmp_p(np,block_size) + Ltmp_q(nq,block_size)
-         
+
          if (mem > qp_max_mem) then
            s = s*2.d0
          else
            exit
          endif
-         
+
          if ((s > 1.d0).or.(nq == 0)) then
            call print_memory_usage()
            print *,  'Not enough memory. Reduce cholesky threshold'
            stop -1
          endif
-         
+
        enddo
-       
+
        ! d., e.
        block_size = max(N,24)
-       
+
        L_old => L
        allocate(L(ndim,rank+nq), stat=ierr)
        if (ierr /= 0) then
@@ -197,48 +197,52 @@ END_PROVIDER
          print *,  irp_here, ': allocation failed : (L(ndim,rank+nq))'
          stop -1
        endif
-       
-       !$OMP PARALLEL DO PRIVATE(k)
+
+       !$OMP PARALLEL DO PRIVATE(k,j)
        do k=1,rank
-         L(:,k) = L_old(:,k)
+         do j=1,ndim
+           L(j,k) = L_old(j,k)
+         enddo
        enddo
        !$OMP END PARALLEL DO
-       
+
        deallocate(L_old)
-       
+
        allocate(Delta(np,nq), stat=ierr)
        if (ierr /= 0) then
          call print_memory_usage()
          print *,  irp_here, ': allocation failed : (Delta(np,nq))'
          stop -1
        endif
-       
+
        allocate(Ltmp_p(np,block_size), stat=ierr)
        if (ierr /= 0) then
          call print_memory_usage()
          print *,  irp_here, ': allocation failed : (Ltmp_p(np,block_size))'
          stop -1
        endif
-       
+
        allocate(Ltmp_q(nq,block_size), stat=ierr)
        if (ierr /= 0) then
          call print_memory_usage()
          print *,  irp_here, ': allocation failed : (Ltmp_q(nq,block_size))'
          stop -1
        endif
-       
-       
+
+
        allocate(computed(nq))
-       
+
        !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(m,k,p,q,j)
-       
+
        !$OMP DO
        do q=1,nq
-         Delta(:,q) = 0.d0
+         do j=1,np
+           Delta(j,q) = 0.d0
+         enddo
          computed(q) = .False.
        enddo
        !$OMP ENDDO NOWAIT
-       
+
        !$OMP DO
        do k=1,N
          do p=1,np
@@ -249,36 +253,36 @@ END_PROVIDER
          enddo
        enddo
        !$OMP END DO NOWAIT
-       
+
        !$OMP BARRIER
        !$OMP END PARALLEL
-       
+
        if (N>0) then
          call dgemm('N','T', np, nq, N, -1.d0,                       &
              Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np)
        endif
-       
+
        ! f.
        Qmax = D(Dset(1))
        do q=1,nq
          Qmax = max(Qmax, D(Dset(q)))
        enddo
-       
+
        ! g.
-       
+
        iblock = 0
        do j=1,nq
-         
+
          if ( (Qmax <= Dmin).or.(N+j > ndim) ) exit
          ! i.
          rank = N+j
-         
+
          if (iblock == block_size) then
            call dgemm('N','T',np,nq,block_size,-1.d0,                &
                Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np)
            iblock = 0
          endif
-         
+
          ! ii.
          do dj=1,nq
            qj = Dset(dj)
@@ -286,9 +290,9 @@ END_PROVIDER
              exit
            endif
          enddo
-         
+
          L(1:ndim, rank) = 0.d0
-         
+
          if (.not.computed(dj)) then
            m = dj
            !$OMP PARALLEL DO PRIVATE(k) SCHEDULE(guided)
@@ -314,16 +318,16 @@ END_PROVIDER
          do p=1,np
            Ltmp_p(p,iblock) = Delta(p,dj)
          enddo
-         
+
          ! iv.
          if (iblock > 1) then
            call dgemv('N', np, iblock-1, -1.d0, Ltmp_p, np, Ltmp_q(dj,1), nq, 1.d0,&
                Ltmp_p(1,iblock), 1)
          endif
-         
+
          ! iii.
          f = 1.d0/dsqrt(Qmax)
-         
+
          !$OMP PARALLEL PRIVATE(m,p,q,k) DEFAULT(shared)
          !$OMP DO
          do p=1,np
@@ -332,38 +336,38 @@ END_PROVIDER
            D(Lset(p)) = D(Lset(p)) - Ltmp_p(p,iblock) * Ltmp_p(p,iblock)
          enddo
          !$OMP END DO
-         
+
          !$OMP DO
          do q=1,nq
            Ltmp_q(q,iblock) = L(Dset(q), rank)
          enddo
          !$OMP END DO
-         
+
          !$OMP END PARALLEL
-         
+
          Qmax = D(Dset(1))
          do q=1,nq
            Qmax = max(Qmax, D(Dset(q)))
          enddo
-         
+
        enddo
-       
+
        print '(I10, 4X, ES12.3)', rank, Qmax
-       
+
        deallocate(computed)
        deallocate(Delta)
        deallocate(Ltmp_p)
        deallocate(Ltmp_q)
-       
+
        ! i.
        N = rank
-       
+
        ! j.
        Dmax = D(Lset(1))
        do p=1,np
          Dmax = max(Dmax, D(Lset(p)))
        enddo
-       
+
        np=0
        do p=1,ndim
          if ( dscale*dscale*Dmax*D(p) > tau*tau ) then
@@ -371,9 +375,9 @@ END_PROVIDER
            Lset(np) = p
          endif
        enddo
-       
+
      enddo
-     
+
      allocate(cholesky_ao(ao_num,ao_num,rank), stat=ierr)
      if (ierr /= 0) then
        call print_memory_usage()
@@ -387,10 +391,10 @@ END_PROVIDER
      !$OMP END PARALLEL DO
      deallocate(L)
      cholesky_ao_num = rank
-     
+
      print *,  '============ ============='
      print *,  ''
-     
+
      if (write_ao_cholesky) then
        print *,  'Writing Cholesky vectors to disk...'
        iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'W')
@@ -401,9 +405,9 @@ END_PROVIDER
      endif
 
    endif
-   
+
    print *, 'Rank  : ', cholesky_ao_num, '(', 100.d0*dble(cholesky_ao_num)/dble(ao_num*ao_num), ' %)'
    print *,  ''
-   
+
 END_PROVIDER
- 
+
diff --git a/src/ccsd/ccsd_space_orb_sub_chol.irp.f b/src/ccsd/ccsd_space_orb_sub_chol.irp.f
index 0ba46e56..54ebff73 100644
--- a/src/ccsd/ccsd_space_orb_sub_chol.irp.f
+++ b/src/ccsd/ccsd_space_orb_sub_chol.irp.f
@@ -252,7 +252,8 @@ subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
     do a = 1, nV
       do j = 1, nO
         do i = 1, nO
-          W_oovo(i,j,a,u) = 2d0 * cc_space_v_vooo(a,u,i,j) - cc_space_v_vooo(a,u,j,i)
+!          W_oovo(i,j,a,u) = 2d0 * cc_space_v_vooo(a,u,i,j) - cc_space_v_vooo(a,u,j,i)
+          W_oovo(i,j,a,u) = 2d0 * cc_space_v_oovo(i,j,a,u) - cc_space_v_oovo(j,i,a,u)
         enddo
       enddo
     enddo
@@ -514,10 +515,10 @@ subroutine compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
                 cc_space_v_vv_chol(1,1,gam), cholesky_mo_num, &
                 0.d0, tmpB1, nV*block_size)
 
-        call dgemm('T','N', nV*min(block_size, nV-iblock+1), nV, cholesky_mo_num, 1.d0, &
-          cc_space_v_vv_chol(1,1,iblock), cholesky_mo_num, &
-          tmp_cc2, cholesky_mo_num, &
-          1.d0, tmpB1, nV*block_size)
+        call dgemm('T','N', nV*min(block_size, nV-iblock+1), nV, cholesky_mo_num, &
+                1.d0, cc_space_v_vv_chol(1,1,iblock), cholesky_mo_num, &
+                tmp_cc2, cholesky_mo_num, &
+                1.d0, tmpB1, nV*block_size)
 
         do beta = iblock, min(nV, iblock+block_size-1)
           do b = 1, nV
@@ -1107,37 +1108,27 @@ subroutine compute_g_occ_chol(nO,nV,t1,t2,H_oo,g_occ)
   double precision, intent(in)  :: t2(nO, nO, nV, nV)
   double precision, intent(out) :: g_occ(nO, nO)
 
-  integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta, gam
+  g_occ = H_oo
 
   call dgemm('N','N',nO,nO,nV, &
              1d0, t1, size(t1,1), &
                   cc_space_f_vo, size(cc_space_f_vo,1), &
-             0d0, g_occ, size(g_occ,1))
+             1d0, g_occ, size(g_occ,1))
 
-  !$omp parallel &
-  !$omp shared(nO,nV,g_occ,H_oo, cc_space_v_ovoo,t1) &
-  !$omp private(i,j,a,u) &
-  !$omp default(none)
-  !$omp do
-  do i = 1, nO
-    do u = 1, nO
-      g_occ(u,i) = g_occ(u,i) + H_oo(u,i)
-    enddo
-  enddo
-  !$omp end do
+  double precision, allocatable :: X(:)
+  allocate(X(cholesky_mo_num))
+  call dgemv('N',cholesky_mo_num,nO*nV,2.d0, &
+    cc_space_v_ov_chol, cholesky_mo_num, &
+    t1, 1, 0.d0, X, 1)
 
-  !$omp do
-  do i = 1, nO
-    do j = 1, nO
-      do a = 1, nV
-        do u = 1, nO
-          g_occ(u,i) = g_occ(u,i) + (2d0 * cc_space_v_ovoo(u,a,i,j) - cc_space_v_ovoo(u,a,j,i)) * t1(j,a)
-        enddo
-      enddo
-    enddo
-  enddo
-  !$omp end do
-  !$omp end parallel
+  call dgemv('T',cholesky_mo_num,nO*nO,1.d0, &
+    cc_space_v_oo_chol, cholesky_mo_num, &
+    X, 1, 1.d0, g_occ, 1)
+  deallocate(X)
+
+  call dgemv('T',nO*nV,nO*nO,-1.d0, &
+    cc_space_v_ovoo, nO*nV, &
+    t1, 1, 1.d0, g_occ, 1)
 
 end