subroutine get_mrsf_transition_density(infos, trden, bvec_mo, ist, jst)
use precision, only: dp
use messages, only: show_message, with_abort
use types, only: information
implicit none
type(information), intent(in) :: infos
real(kind=dp), intent(out), dimension(:,:) :: trden
real(kind=dp), intent(in), dimension(:,:) :: bvec_mo
integer, intent(in) :: ist, jst
real(kind=dp), allocatable :: xv12i(:,:), xv12j(:,:)
integer :: nbf, noca, nocb, nvirb, ok
integer :: i, ij, ijd, ijg, ijlr1, ijlr2, j, mrst
real(kind=dp), parameter :: rsqrt = 1.0_dp/sqrt(2.0_dp)
nbf = infos%basis%nbf
noca = infos%mol_prop%nelec_a
mrst = infos%tddft%mult
nocb = noca-2
nvirb = nbf-nocb
allocate(xv12i(noca,nvirb), &
xv12j(noca,nvirb), &
source=0.0_dp, stat=ok)
if (ok /= 0) call show_message('Cannot allocate memory', with_abort)
ijlr1 = (noca-1-nocb-1)*noca+noca-1
ijg = (noca-1-nocb-1)*noca+noca
ijd = (noca-nocb-1)*noca+noca-1
ijlr2 = (noca-nocb-1)*noca+noca
if (mrst==1) then
do i = 1, noca
do j = nocb + 1, nbf
ij = (j-nocb-1)*noca+i
if (ij==ijlr1) then
xv12j(i,j-nocb) = bvec_mo(ijlr1,jst)*rsqrt
xv12i(i,j-nocb) = bvec_mo(ijlr1,ist)*rsqrt
cycle
else if (ij==ijlr2) then
xv12j(i,j-nocb) = -bvec_mo(ijlr1,jst)*rsqrt
xv12i(i,j-nocb) = -bvec_mo(ijlr1,ist)*rsqrt
cycle
end if
xv12j(i,j-nocb) = bvec_mo(ij,jst)
xv12i(i,j-nocb) = bvec_mo(ij,ist)
end do
end do
else if (mrst==3) then
do i = 1, noca
do j = nocb+1, nbf
ij = (j-nocb-1)*noca+i
if (ij==ijlr1) then
xv12j(i,j-nocb) = bvec_mo(ijlr1,jst)*rsqrt
xv12i(i,j-nocb) = bvec_mo(ijlr1,ist)*rsqrt
cycle
else if (ij==ijlr2) then
xv12j(i,j-nocb) = bvec_mo(ijlr1,jst)*rsqrt
xv12i(i,j-nocb) = bvec_mo(ijlr1,ist)*rsqrt
cycle
else if (ij==ijg) then
xv12j(i,j-nocb) = 0.0_dp
xv12i(i,j-nocb) = 0.0_dp
cycle
else if (ij==ijd) then
xv12j(i,j-nocb) = 0.0_dp
xv12i(i,j-nocb) = 0.0_dp
cycle
end if
xv12j(i,j-nocb) = bvec_mo(ij,jst)
xv12i(i,j-nocb) = bvec_mo(ij,ist)
end do
end do
end if
call get_trans_den(trden, xv12i, xv12j, noca, nocb, nvirb)
return
end subroutine get_mrsf_transition_density