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