BEGIN_PROVIDER [ double precision, r2_cc, (spin_occ_num,spin_occ_num,spin_vir_num,spin_vir_num) ]
 implicit none
 BEGIN_DOC
  ! t2 residues in non-canonical CCSD
 END_DOC
  integer                        :: i,j,m,n
  integer                        :: a,b,e,f
  double precision               :: x
  double precision, allocatable  :: mat_A(:,:), mat_B(:,:), mat_C(:,:)

  r2_cc(:,:,:,:) = OOVV(:,:,:,:)

  allocate ( mat_A(spin_vir_num,spin_vir_num), mat_C(spin_vir_num,spin_vir_num), & 
             mat_B(spin_vir_num,spin_occ_num) )

  do j=1,spin_occ_num
    do i=1,spin_occ_num
      mat_A(:,:) = t2_cc(i,j,:,:)
      call dgemm('N', 'T',                                           &
          spin_vir_num, spin_vir_num, spin_vir_num,                  &
          1.d0, mat_A, size(mat_A,1),                                &
          c_spin_fock_matrix_mo_vv, size(c_spin_fock_matrix_mo_vv,1),&
          0.d0, mat_C, size(mat_C,1) )

      call dgemm('N', 'T',                                           &
          spin_vir_num, spin_occ_num, spin_vir_num,                  &
          -0.5d0, mat_A, size(mat_A,1),                              &
          c_spin_fock_matrix_mo_ov, size(c_spin_fock_matrix_mo_ov,1),&
          0.d0, mat_B, size(mat_B,1) )
      
      call dgemm('N', 'N',                                           &
          spin_vir_num, spin_vir_num, spin_occ_num,                  &
          1.d0, mat_B, size(mat_B,1),                                &
          t1_cc, size(t1_cc,1),                                      &
          1.d0, mat_C, size(mat_C,1) )

      do b=1,spin_vir_num
        do a=1,spin_vir_num
          r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + mat_C(a,b) - mat_C(b,a)
        enddo
      enddo

    enddo
  enddo


  deallocate(mat_A,mat_B,mat_C)

  allocate ( mat_A(spin_occ_num,spin_occ_num), mat_C(spin_occ_num,spin_occ_num), & 
             mat_B(spin_occ_num,spin_vir_num) )

  do b=1,spin_vir_num
    do a=1,spin_vir_num

      call dgemm('N', 'N',                                           &
          spin_occ_num, spin_occ_num, spin_occ_num,                  &
          1.d0, t2_cc(1,1,a,b), size(t2_cc,1),                       &
          c_spin_fock_matrix_mo_oo, size(c_spin_fock_matrix_mo_oo,1),&
          0.d0, mat_C, size(mat_C,1) )


      call dgemm('N', 'N',                                           &
          spin_occ_num, spin_vir_num, spin_occ_num,                  &
          0.5d0, t2_cc(1,1,a,b), size(t2_cc,1),                      &
          c_spin_fock_matrix_mo_ov, size(c_spin_fock_matrix_mo_ov,1),&
          0.d0, mat_B, size(mat_B,1) )
      
      call dgemm('N', 'T',                                           &
          spin_occ_num, spin_occ_num, spin_vir_num,                  &
          1.d0, mat_B, size(mat_B,1),                                &
          t1_cc, size(t1_cc,1),                                      &
          1.d0, mat_C, size(mat_C,1) )

      do j=1,spin_occ_num
        do i=1,spin_occ_num
          r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + mat_C(j,i) - mat_C(i,j)
        enddo
      enddo

    end do
  end do

  deallocate(mat_A,mat_B,mat_C)

  call dgemm('T', 'N',                                               &
      spin_occ_num*spin_occ_num, spin_vir_num*spin_vir_num, spin_occ_num*spin_occ_num,&
      0.5d0, cWoooo, spin_occ_num*spin_occ_num,                      &
      tau_cc, spin_occ_num*spin_occ_num,                             &
      1.d0, r2_cc, spin_occ_num*spin_occ_num)

  call dgemm('N', 'T',                                               &
      spin_occ_num*spin_occ_num, spin_vir_num*spin_vir_num, spin_vir_num*spin_vir_num,&
      0.5d0, tau_cc, spin_occ_num*spin_occ_num,                      &
      cWvvvv, spin_vir_num*spin_vir_num,                             &
      1.d0, r2_cc, spin_occ_num*spin_occ_num)


  double precision, allocatable :: tmp(:,:,:,:), tmp2(:,:,:,:)
  allocate(tmp(spin_occ_num,spin_vir_num,spin_vir_num,spin_occ_num), &
           tmp2(spin_occ_num,spin_vir_num,spin_occ_num,spin_vir_num) )

  
  call dgemm('N', 'N',                                               &
      spin_occ_num*spin_vir_num, spin_occ_num*spin_vir_num, spin_occ_num*spin_vir_num,&
      1.0d0, t2_cc2, spin_occ_num*spin_vir_num,                      &
      cWovvo_prime, spin_occ_num*spin_vir_num,                       &
      0.d0, tmp, spin_occ_num*spin_vir_num)


  do e=1,spin_vir_num
    do m=1,spin_occ_num
      do a=1,spin_vir_num
        do i=1,spin_occ_num
          tmp2(i,a,m,e) = -t1_cc(i,e)*t1_cc(m,a)
        enddo
      enddo
    enddo
  enddo

  call dgemm('N', 'N',                                               &
      spin_occ_num*spin_vir_num, spin_occ_num*spin_vir_num, spin_occ_num*spin_vir_num,&
      1.0d0, tmp2, spin_occ_num*spin_vir_num,                        &
      OVVO_prime, spin_occ_num*spin_vir_num,                         &
      1.d0, tmp, spin_occ_num*spin_vir_num)
  

  do b=1,spin_vir_num
    do a=1,spin_vir_num
      do j=1,spin_occ_num
        do i=1,spin_occ_num
          r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + tmp(i,a,b,j)
          r2_cc(j,i,a,b) = r2_cc(j,i,a,b) - tmp(i,a,b,j)
          r2_cc(i,j,b,a) = r2_cc(i,j,b,a) - tmp(i,a,b,j)
          r2_cc(j,i,b,a) = r2_cc(j,i,b,a) + tmp(i,a,b,j)
        enddo
      enddo
    enddo
  enddo



  deallocate(tmp,tmp2)

  allocate( mat_A(spin_vir_num,spin_occ_num), mat_C(spin_occ_num,spin_occ_num) )
!  
!  do b=1,spin_vir_num
!    do a=1,spin_vir_num
!
!      do j=1,spin_occ_num
!        do e=1,spin_vir_num
!          mat_A(e,j) = VVVO(a,b,e,j)
!        enddo
!      enddo
!
!      call dgemm('N', 'N',                                           &
!          spin_occ_num, spin_occ_num, spin_vir_num,                  &
!          1.0d0, t1_cc, spin_occ_num,                                &
!          mat_A, size(mat_A,1),                                      &
!          0.d0, mat_C, spin_occ_num)
!
!      do j=1,spin_occ_num
!        do i=1,spin_occ_num
!          r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + mat_C(i,j) - mat_C(j,i)
!        enddo
!      endd
!
!    enddo
!  enddo


  do b=1,spin_vir_num
    do a=1,spin_vir_num

      do j=1,spin_occ_num
        do e=1,spin_vir_num
          mat_A(e,j) = VVVO(a,b,e,j)
        enddo
      enddo
      call dgemm('N', 'N',                                           &
          spin_occ_num, spin_occ_num, spin_vir_num,                  &
          1.0d0, t1_cc, spin_occ_num,                                &
          mat_A, size(mat_A,1),                                      &
          0.d0, mat_C, spin_occ_num)
      do j=1,spin_occ_num
        do i=1,spin_occ_num
          r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + mat_C(i,j) - mat_C(j,i)
        enddo
      enddo

    end do
  end do


  do b=1,spin_vir_num
    do a=1,spin_vir_num
      do j=1,spin_occ_num
        do i=1,spin_occ_num
          do m=1,spin_occ_num
            r2_cc(i,j,a,b) = r2_cc(i,j,a,b) - t1_cc(m,a)*OVOO(m,b,i,j) 
            r2_cc(i,j,a,b) = r2_cc(i,j,a,b) + t1_cc(m,b)*OVOO(m,a,i,j) 
          end do
        end do
      end do
     end do
   end do


! Final expression of the t2 residue

  r2_cc(:,:,:,:) = delta_oovv(:,:,:,:)*t2_cc(:,:,:,:) - r2_cc(:,:,:,:)

END_PROVIDER