subroutine form_r2_nc(t1,t2,tau,cFoo,cFov,cFvv,cWoooo,cWvvvv,cWovvo,r2)

! Form t2 residues in non-canonical CCSD

  implicit none

! Input variables

  double precision,intent(in)    :: cFoo(spin_occ_num,spin_occ_num)
  double precision,intent(in)    :: cFov(spin_occ_num,spin_vir_num)
  double precision,intent(in)    :: cFvv(spin_vir_num,spin_vir_num)
  
  double precision,intent(in)    :: cWoooo(spin_occ_num,spin_occ_num,spin_occ_num,spin_occ_num)
  double precision,intent(in)    :: cWvvvv(spin_vir_num,spin_vir_num,spin_vir_num,spin_vir_num)
  double precision,intent(in)    :: cWovvo(spin_occ_num,spin_vir_num,spin_vir_num,spin_occ_num)
  
  double precision,intent(in)    :: t1(spin_occ_num,spin_vir_num)
  double precision,intent(in)    :: t2(spin_occ_num,spin_occ_num,spin_vir_num,spin_vir_num)
  double precision,intent(in)    :: tau(spin_occ_num,spin_occ_num,spin_vir_num,spin_vir_num)
  double precision               :: x

! Local variables

  integer                        :: i,j,m,n
  integer                        :: a,b,e,f

! Output variables

  double precision,intent(out)   :: r2(spin_occ_num,spin_occ_num,spin_vir_num,spin_vir_num)

  r2(:,:,:,:) = OOVV(:,:,:,:)

  do e=1,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(i,j,a,b) = r2(i,j,a,b) + t2(i,j,a,e)*cFvv(b,e)
          end do
        end do
      end do
    end do
  end do


  do e=1,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(i,j,a,b) = r2(i,j,a,b) - t2(i,j,b,e)*cFvv(a,e)
          end do
        end do
      end do
    end do
  end do


 do b=1,spin_vir_num
   do e=1,spin_vir_num
     do a=1,spin_vir_num
       do j=1,spin_occ_num
         do i=1,spin_occ_num
           x = 0.5d0*t2(i,j,a,e)
           if (x /= 0.d0) then
              do m=1,spin_occ_num
                r2(i,j,a,b) = r2(i,j,a,b) - x*t1(m,b)*cFov(m,e)
              end do
            endif
          end do
        end do
      end do
    end do
  end do


  do e=1,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
            x = 0.5d0*t2(i,j,b,e)
            if (x /= 0.d0) then
              do m=1,spin_occ_num
                r2(i,j,a,b) = r2(i,j,a,b) + x*t1(m,a)*cFov(m,e)
              end do
            endif
          end do
        end do
      end do
    end do
  end do
  

  do b=1,spin_vir_num
    do a=1,spin_vir_num
      do j=1,spin_occ_num
        do m=1,spin_occ_num
          x = cFoo(m,j)
          do i=1,spin_occ_num
            r2(i,j,a,b) = r2(i,j,a,b) - t2(i,m,a,b)*x
          end do
        end do
      end do
    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(i,j,a,b) = r2(i,j,a,b) + t2(j,m,a,b)*cFoo(m,i)
          end do
        end do
      end do
    end do
  end do


  do b=1,spin_vir_num
    do a=1,spin_vir_num
      do e=1,spin_vir_num
        do m=1,spin_occ_num
          x = 0.5d0*cFov(m,e)
          if (x == 0.d0) cycle
          do j=1,spin_occ_num
            do i=1,spin_occ_num
              r2(i,j,a,b) = r2(i,j,a,b) - x*t2(i,m,a,b)*t1(j,e)
            end do
          end do
        end do
      end do
    end do
  end do


  do b=1,spin_vir_num
    do a=1,spin_vir_num
      do e=1,spin_vir_num
        do m=1,spin_occ_num
          x = 0.5d0*cFov(m,e)
          do j=1,spin_occ_num
            do i=1,spin_occ_num
              r2(i,j,a,b) = r2(i,j,a,b) + x*t2(j,m,a,b)*t1(i,e)
            end do
          end do
        end do
      end do
    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 n=1,spin_occ_num
            do m=1,spin_occ_num
              r2(i,j,a,b) = r2(i,j,a,b) + 0.5d0*tau(m,n,a,b)*cWoooo(m,n,i,j)
            end do
          end do
        end do
      end do
    end do
  end do


  do f=1,spin_vir_num
    do e=1,spin_vir_num
      do b=1,spin_vir_num
        do a=1,spin_vir_num
          x = 0.5d0*cWvvvv(a,b,e,f) 
          do j=1,spin_occ_num
            do i=1,spin_occ_num
              r2(i,j,a,b) = r2(i,j,a,b) + x*tau(i,j,e,f)
            end do
          end do
        end do
      end do
    end do
  end do


  do b=1,spin_vir_num
    do a=1,spin_vir_num
      do j=1,spin_occ_num
        do e=1,spin_vir_num
          do m=1,spin_occ_num
            do i=1,spin_occ_num
              r2(i,j,a,b) = r2(i,j,a,b) + t2(i,m,a,e)*cWovvo(m,b,e,j) 
            end do
          end do
        end do
      end do
    end do
  end do

  do b=1,spin_vir_num
    do a=1,spin_vir_num
      do j=1,spin_occ_num
        do e=1,spin_vir_num
          do m=1,spin_occ_num
            x = t1(m,a)*OVVO(m,b,e,j)
            do i=1,spin_occ_num
              r2(i,j,a,b) = r2(i,j,a,b) - t1(i,e)*x
            end do
          end do
        end do
      end do
    end do
  end do

  do b=1,spin_vir_num
    do e=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(i,j,a,b) = r2(i,j,a,b) - t2(j,m,a,e)*cWovvo(m,b,e,i) 
            end do
          end do
        end do
      end do
    end do
  end do

  do b=1,spin_vir_num
    do a=1,spin_vir_num
      do e=1,spin_vir_num
        do j=1,spin_occ_num
          do i=1,spin_occ_num
            do m=1,spin_occ_num
              r2(i,j,a,b) = r2(i,j,a,b) + t1(j,e)*t1(m,a)*OVVO(m,b,e,i)
            end do
          end do
        end do
      end do
    end do
  end do

 do e=1,spin_vir_num
   do b=1,spin_vir_num
     do a=1,spin_vir_num
       do j=1,spin_occ_num
         do m=1,spin_occ_num
           do i=1,spin_occ_num
              r2(i,j,a,b) = r2(i,j,a,b) - t2(i,m,b,e)*cWovvo(m,a,e,j) 
            end do
          end do
        end do
      end do
    end do
  end do

  do b=1,spin_vir_num
    do a=1,spin_vir_num
      do j=1,spin_occ_num
        do e=1,spin_vir_num
          do i=1,spin_occ_num
            do m=1,spin_occ_num
              r2(i,j,a,b) = r2(i,j,a,b) + t1(i,e)*t1(m,b)*OVVO(m,a,e,j) 
            end do
          end do
        end do
      end do
    end do
  end do

  do e=1,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
            do m=1,spin_occ_num
              r2(i,j,a,b) = r2(i,j,a,b) + t2(j,m,b,e)*cWovvo(m,a,e,i) 
            end do
          end do
        end do
      end do
    end do
  end do

  do b=1,spin_vir_num
    do e=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(i,j,a,b) = r2(i,j,a,b) - t1(j,e)*t1(m,b)*OVVO(m,a,e,i)
            end do
          end do
        end do
      end do
    end do
  end do


  do b=1,spin_vir_num
    do a=1,spin_vir_num
      do j=1,spin_occ_num
        do e=1,spin_vir_num
          do i=1,spin_occ_num
            r2(i,j,a,b) = r2(i,j,a,b) + t1(i,e)*VVVO(a,b,e,j) 
          end do
        end do
      end do
    end do
  end do


  do b=1,spin_vir_num
    do a=1,spin_vir_num
      do e=1,spin_vir_num
        do i=1,spin_occ_num
          do j=1,spin_occ_num
            r2(i,j,a,b) = r2(i,j,a,b) - t1(j,e)*VVVO(a,b,e,i)
          end do
        end do
      end do
    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(i,j,a,b) = r2(i,j,a,b) - t1(m,a)*OVOO(m,b,i,j) 
          end do
        end do
      end do
    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(i,j,a,b) = r2(i,j,a,b) + t1(m,b)*OVOO(m,a,i,j)
          end do
        end do
      end do
    end do
  end do

! Final expression of the t2 residue

  r2(:,:,:,:) = delta_oovv(:,:,:,:)*t2(:,:,:,:) - r2(:,:,:,:)

end subroutine form_r2_nc