subroutine form_abh(t1,tau,aoooo,bvvvv,hovvo)

! Scuseria Eqs. (11),(12) and (13)

  implicit none

! Input variables

  double precision,intent(in)   :: t1(spin_occ_num,spin_vir_num)
  double precision,intent(in)   :: tau(spin_occ_num,spin_occ_num,spin_vir_num,spin_vir_num)

! Local variables

  integer                       :: i,j,k,l
  integer                       :: a,b,c,d

! Output variables

  double precision,intent(out)  :: aoooo(spin_occ_num,spin_occ_num,spin_occ_num,spin_occ_num)
  double precision,intent(out)  :: bvvvv(spin_vir_num,spin_vir_num,spin_vir_num,spin_vir_num)
  double precision,intent(out)  :: hovvo(spin_occ_num,spin_vir_num,spin_vir_num,spin_occ_num)

  aoooo(:,:,:,:) = OOOO(:,:,:,:)

  do l=1,spin_occ_num
    do k=1,spin_occ_num
      do j=1,spin_occ_num
        do i=1,spin_occ_num

          do c=1,spin_vir_num
            aoooo(i,j,k,l) = aoooo(i,j,k,l) + OVOO(i,c,k,l)*t1(j,c)
          end do

          do c=1,spin_vir_num
            aoooo(i,j,k,l) = aoooo(i,j,k,l) - OVOO(j,c,k,l)*t1(i,c)
          end do

          do d=1,spin_vir_num
            do c=1,spin_vir_num
              aoooo(i,j,k,l) = aoooo(i,j,k,l) + OOVV(k,l,c,d)*tau(i,j,c,d)
            end do
          end do

        end do
      end do
    end do
  end do

  bvvvv(:,:,:,:) = VVVV(:,:,:,:)

  do b=1,spin_vir_num 
    do a=1,spin_vir_num
      do d=1,spin_vir_num
        do c=1,spin_vir_num
          
          do k=1,spin_occ_num
            bvvvv(c,d,a,b) = bvvvv(c,d,a,b) - VOVV(a,k,c,d)*t1(k,b)
          end do

          do k=1,spin_occ_num
            bvvvv(c,d,a,b) = bvvvv(c,d,a,b) + VOVV(b,k,c,d)*t1(k,a)
          end do

        end do
      end do
    end do
  end do

  hovvo(:,:,:,:) = OVVO(:,:,:,:)

  do k=1,spin_occ_num 
    do a=1,spin_vir_num
      do c=1,spin_vir_num
        do i=1,spin_occ_num

          do l=1,spin_occ_num
            hovvo(i,c,a,k) = hovvo(i,c,a,k) - OVOO(i,c,l,k)*t1(l,a)
          end do

          do d=1,spin_vir_num
            hovvo(i,c,a,k) = hovvo(i,c,a,k) + OVVV(k,a,c,d)*t1(i,d)
          end do

          do d=1,spin_vir_num
            do l=1,spin_occ_num
              hovvo(i,c,a,k) = hovvo(i,c,a,k) - OOVV(k,l,c,d)*tau(i,l,d,a)
            end do
          end do

        end do
      end do
    end do
  end do

end subroutine form_abh