sfrorhs Subroutine

  subroutine sfrorhs(rhs,xhxa,xhxb,hpta,hptb,Tij,Tab,Fa,Fb, &
                     noca,nocb)
    use precision, only: dp

    implicit none

    real(kind=dp), intent(out), dimension(:) :: rhs
    real(kind=dp), intent(inout), dimension(:,:) :: xhxa, xhxb
    real(kind=dp), intent(in), dimension(:,:) :: hpta
    real(kind=dp), intent(in), dimension(:,:) :: hptb
    real(kind=dp), intent(in), dimension(:,:) :: tij
    real(kind=dp), intent(in), dimension(:,:) :: tab
    real(kind=dp), intent(in), dimension(:,:) :: fa, fb
    integer, intent(in) :: noca, nocb

    real(kind=dp), allocatable, dimension(:,:) :: wrk
    integer :: nbf, i, j, ij, k, nconf

    nbf = ubound(fa, 1)

    allocate(wrk(nbf,nbf), &
             source=0.0_dp)

  ! HPTA --> AB1_MO(1)
  ! HPTB --> AB1_MO(2)
  ! TA   --> TIJ
  ! TB   --> TAB

  ! Alpha
  ! XHXA+= 2*FA(P+,I+)*TA(I+,J+)
    call dgemm('n', 'n', nbf, noca, noca, &
               2.0_dp, fa, nbf, &
                       tij, noca, &
               1.0_dp, xhxa, nbf)

  ! Beta
  ! XHXB+= 2*FB(P-,A-)*TB(A-,B-)
    do j = nocb+1, nbf
      do i = nocb+1, nbf
        wrk(i,j) = tab(i-nocb,j-nocb)
      end do
    end do
    call dgemm('n', 'n', nbf, nbf, nbf, &
               2.0_dp, fb, nbf, &
                       wrk, nbf, &
               1.0_dp, xhxb, nbf)

  ! doc-socc
    ij = 0
    do i = nocb+1, noca
      do j = 1, nocb
        ij = ij+1
        rhs(ij) = hptb(j,i-nocb)+xhxa(i,j)-xhxa(j,i)-xhxb(j,i)
      end do
    end do

  ! doc-virt
    do k = noca+1, nbf
      do j = 1, nocb
        ij = ij+1
        rhs(ij) = hpta(j,k-noca)+hptb(j,k-nocb)+xhxa(k,j)-xhxb(j,k)
      end do
    end do

  ! soc-virt
    do k = noca+1, nbf
      do i = nocb+1, noca
        ij = ij+1
        rhs(ij) = hpta(i,k-noca)+xhxa(k,i)+xhxb(k,i)-xhxb(i,k)
      end do
    end do

  ! Multiplied by -1 i.e., RHS of Z-vector eq. -----
    nconf = ij
    rhs(1:nconf) = -rhs(1:nconf)

  end subroutine sfrorhs