subroutine sfdmat(bvec,abxc,mo_a,ta,tb, &
                    noca,nocb)
    use precision, only: dp
    use tdhf_lib, only: iatogen
    use mathlib, only: pack_matrix
    use mathlib, only: orthogonal_transform
    implicit none
    real(kind=dp), intent(in), dimension(:) :: bvec
    real(kind=dp), intent(in), dimension(:,:) :: mo_a
    real(kind=dp), intent(inout), dimension(:,:) :: abxc
    real(kind=dp), intent(out), dimension(:) :: ta, tb
    integer, intent(in) :: noca, nocb
    integer :: nvirb, nbf, nbf_tri, xvec_dim
    real(kind=dp), allocatable, dimension(:,:) :: scr1, scr2
    nbf = ubound(mo_a,1)
    nbf_tri = ubound(ta,1)
    xvec_dim = ubound(bvec,1)
    allocate(scr1(nbf,nbf), &
             scr2(nbf,nbf), &
             source=0.0_dp)
  ! MO(I+,A-) -> AO(M,N)
    nvirb = nbf-nocb
    call iatogen(bvec,scr1,noca,nocb)
    call orthogonal_transform('t', nbf, mo_a, scr1, abxc, scr2)
  ! Unrelaxed difference density matrix -----
  ! OCC(Alpha)-OCC(Alpha)
    call dgemm('n','t',noca,noca,nvirb, &
              -1.0_dp,bvec,noca, &
                      bvec,noca, &
               0.0_dp,scr1,noca)
  ! MO(I+,J+) -> AO(M,N)
    call dgemm('n','n',nbf,noca,noca, &
               1.0_dp,mo_a,nbf, &
                      scr1,noca, &
               0.0_dp,scr2,nbf)
    call dgemm('n','t',nbf,nbf,noca, &
               1.0_dp,scr2,nbf, &
                      mo_a,nbf, &
               0.0_dp,scr1,nbf)
    call pack_matrix(scr1,ta)
    call dgemm('t','n',nvirb,nvirb,noca, &
               1.0_dp,bvec,noca, &
                      bvec,noca, &
               0.0_dp,scr1,nvirb)
  ! MO(A-,B-) -> AO(M,N)
    call dgemm('n','n',nbf,nvirb,nvirb, &
               1.0_dp,mo_a(:,nocb+1:),nbf, &
                      scr1,nvirb, &
               0.0_dp,scr2,nbf)
    call dgemm('n','t',nbf,nbf,nvirb, &
               1.0_dp,scr2,nbf, &
                      mo_a(:,nocb+1:),nbf, &
               0.0_dp,scr1,nbf)
    call pack_matrix(scr1,tb)
    deallocate(scr1,scr2)
  end subroutine sfdmat