subroutine mrsfesum(infos, wrk, fij, fab, pmo, iv)
use precision, only: dp
use types, only: information
use messages, only: show_message, with_abort
use io_constants, only: iw
implicit none
type(information), intent(in) :: infos
real(kind=dp), intent(in), dimension(:,:) :: &
wrk, fij, fab
real(kind=dp), intent(inout), dimension(:,:) :: &
pmo
integer, intent(in) :: iv
real(kind=dp), allocatable, dimension(:,:) :: scr, tmp1, wrk1
real(kind=dp) :: dumn, xlr
integer :: nbf, nocca, noccb, mrst, i, ij, j, lr1, lr2, ok
real(kind=dp), parameter :: sqrt2 = 1.0_dp/sqrt(2.0_dp)
logical :: debug_mode
nbf = infos%basis%nbf
nocca = infos%mol_prop%nelec_a
noccb = infos%mol_prop%nelec_b
mrst = infos%tddft%mult
debug_mode = infos%tddft%debug_mode
allocate(scr(nbf,nbf), tmp1(nbf,nbf), wrk1(nbf,nbf), source=0.0_dp, stat=ok)
if (ok /= 0) call show_message('Cannot allocate memory', with_abort)
lr1 = nocca-1
lr2 = nocca
scr = wrk
scr(lr1,lr1) = 0.0_dp
scr(lr2,lr2) = 0.0_dp
! Contraction 1
call dgemm('n', 't', nocca, nbf-noccb, nbf-noccb, &
1.0_dp, scr(1,noccb+1), nbf, &
fab(noccb+1:,noccb+1), nbf, &
0.0_dp, tmp1(1,noccb+1), nbf)
! Contraction 2
call dgemm('n', 'n', nocca, nbf-noccb, nocca, &
-1.0_dp, fij(:,1), nbf, &
scr(:,noccb+1), nbf, &
1.0_dp, tmp1(1,noccb+1), nbf)
xlr = wrk(lr1,lr1)
if (mrst==1) then
do j = noccb+1, nbf
do i = 1, nocca
wrk1(i,j) = wrk1(i,j)+tmp1(i,j)
if(i==lr1) wrk1(i,j) = wrk1(i,j)+fab(j,lr1)*xlr*sqrt2
if(i==lr2) wrk1(i,j) = wrk1(i,j)-fab(j,lr2)*xlr*sqrt2
if(j==lr1) wrk1(i,j) = wrk1(i,j)-fij(i,lr1)*xlr*sqrt2
if(j==lr2) wrk1(i,j) = wrk1(i,j)+fij(i,lr2)*xlr*sqrt2
end do
end do
dumn = - dot_product(fij(lr1,1:nocca),scr(1:nocca,lr1)) &
+ dot_product(fij(lr2,1:nocca),scr(1:nocca,lr2)) &
+ dot_product(fab(lr1,noccb+1:nbf),scr(lr1,noccb+1:nbf)) &
- dot_product(fab(lr2,noccb+1:nbf),scr(lr2,noccb+1:nbf))
wrk1(lr1,lr1) = dumn*sqrt2 &
+ xlr*(fab(lr1,lr1)+fab(lr2,lr2) &
-fij(lr1,lr1)-fij(lr2,lr2))*0.5_dp
elseif (mrst==3) then
do j = noccb+1, nbf
do i = 1, nocca
wrk1(i,j) = wrk1(i,j)+tmp1(i,j)
if(i==lr1) wrk1(i,j) = wrk1(i,j)+fab(j,lr1)*xlr*sqrt2
if(i==lr2) wrk1(i,j) = wrk1(i,j)+fab(j,lr2)*xlr*sqrt2
if(j==lr1) wrk1(i,j) = wrk1(i,j)-fij(i,lr1)*xlr*sqrt2
if(j==lr2) wrk1(i,j) = wrk1(i,j)-fij(i,lr2)*xlr*sqrt2
end do
end do
dumn = dot_product(-fij(lr1,1:nocca),scr(1:nocca,lr1)) &
+ dot_product(-fij(lr2,1:nocca),scr(1:nocca,lr2)) &
+ dot_product( fab(lr1,noccb+1:nbf),scr(lr1,noccb+1:nbf)) &
+ dot_product( fab(lr2,noccb+1:nbf),scr(lr2,noccb+1:nbf))
wrk1(lr1,lr1) = dumn*sqrt2 &
+ xlr*(fab(lr1,lr1)+fab(lr2,lr2) &
-fij(lr1,lr1)-fij(lr2,lr2))*0.5_dp
end if
if (mrst==1) then
wrk1(lr2,lr2) = 0.0_dp
else if (mrst==3) then
wrk1(lr2,lr1) = 0.0_dp
wrk1(lr1,lr2) = 0.0_dp
wrk1(lr2,lr2) = 0.0_dp
end if
ij = 0
do j = noccb+1, nbf
do i = 1, nocca
ij = ij+1
pmo(ij,iv) = pmo(ij,iv)+wrk1(i,j)
end do
end do
if (debug_mode) then
write(iw,*) 'Check sum = ivec', iv
write(iw,*) 'Check sum = pmo', sum(abs(pmo(:,iv)))
end if
return
end subroutine mrsfesum