subroutine mrsfcbc(infos,va,vb,bvec,fmrsf)
use messages, only: show_message, with_abort
use types, only: information
use io_constants, only: iw
use precision, only: dp
implicit none
type(information), intent(in) :: infos
real(kind=dp), intent(in), dimension(:,:) :: &
va, vb, bvec
real(kind=dp), intent(inout), target, dimension(:,:,:) :: &
fmrsf
real(kind=dp), allocatable, dimension(:,:) :: &
tmp, tmp1, tmp2
real(kind=dp), pointer, dimension(:,:) :: &
bo2v, bo1v, bco1, bco2, ball, co12, o21v
integer :: nocca, noccb, mrst, i, j, m, nbf, lr1, lr2, ok
logical :: debug_mode
real(kind=dp), parameter :: isqrt2 = 1.0_dp/sqrt(2.0_dp)
ball => fmrsf(7,:,:)
bo2v => fmrsf(1,:,:)
bo1v => fmrsf(2,:,:)
bco1 => fmrsf(3,:,:)
bco2 => fmrsf(4,:,:)
o21v => fmrsf(5,:,:)
co12 => fmrsf(6,:,:)
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
lr1 = nocca-1
lr2 = nocca
allocate(tmp(nbf,nbf), &
tmp1(nbf,4), &
tmp2(nbf,4), &
source=0.0_dp, stat=ok)
if( ok/=0 ) call show_message('Cannot allocate memory',with_abort)
do j = nocca+1, nbf
tmp1(:,1) = tmp1(:,1)+vb(:,j)*bvec(lr2,j)
tmp1(:,2) = tmp1(:,2)+vb(:,j)*bvec(lr1,j)
tmp1(:,3) = tmp1(:,3)+vb(:,j)*bvec(lr2,j)
tmp1(:,4) = tmp1(:,4)+vb(:,j)*bvec(lr1,j)
end do
do i = 1, noccb
tmp2(:,1) = tmp2(:,1)+va(:,i)*bvec(i,lr1)
tmp2(:,2) = tmp2(:,2)+va(:,i)*bvec(i,lr2)
tmp2(:,3) = tmp2(:,3)+va(:,i)*bvec(i,lr1)
tmp2(:,4) = tmp2(:,4)+va(:,i)*bvec(i,lr2)
end do
do m = 1, nbf
bo2v(:,m) = bo2v(:,m)+va(:,lr2)*tmp1(m,1)
bo1v(:,m) = bo1v(:,m)+va(:,lr1)*tmp1(m,2)
bco1(:,m) = bco1(:,m)+vb(m,lr1)*tmp2(:,1)
bco2(:,m) = bco2(:,m)+vb(m,lr2)*tmp2(:,2)
end do
do m = 1, nbf
o21v(m,:) = o21v(m,:)+va(:,lr1)*tmp1(m,3) &
-va(:,lr2)*tmp1(m,4)
co12(m,:) = co12(m,:)+vb(m,lr2)*tmp2(:,3) &
-vb(m,lr1)*tmp2(:,4)
end do
ball = ball + bo2v + bo1v + bco1 + bco2
call dgemm('n', 't', nbf, noccb, nbf-nocca, &
1.0_dp, vb(:,nocca+1), nbf, &
bvec(:,nocca+1), nbf, &
0.0_dp, tmp, nbf)
call dgemm('n', 't', nbf, nbf, noccb, &
1.0_dp, va, nbf, &
tmp, nbf, &
1.0_dp, ball, nbf)
if (mrst==1) then
do m = 1, nbf
ball(:,m) = ball(:,m) &
+va(:,lr2)*bvec(lr2,lr1)*vb(m,lr1) &
+va(:,lr1)*bvec(lr1,lr2)*vb(m,lr2) &
+(va(:,lr1)*vb(m,lr1)-va(:,lr2)*vb(m,lr2)) &
*bvec(lr1,lr1)*isqrt2
end do
else if (mrst==3) then
do m = 1, nbf
ball(:,m) = ball(:,m) &
+(va(:,lr1)*vb(m,lr1)+va(:,lr2)*vb(m,lr2)) &
*bvec(lr1,lr1)*isqrt2
end do
end if
if (debug_mode) then
write(iw,*) 'Check sum = va', sum(abs(va))
write(iw,*) 'Check sum = vb', sum(abs(vb))
write(iw,*) 'Check sum = bvec', sum(abs(bvec))
write(iw,*) 'Check sum = ball', sum(abs(ball))
write(iw,*) 'Check sum = o21v', sum(abs(o21v))
write(iw,*) 'Check sum = co12', sum(abs(co12))
write(iw,*) 'Check sum = bo2v', sum(abs(bo2v))
write(iw,*) 'Check sum = bo1v', sum(abs(bo1v))
write(iw,*) 'Check sum = bco1', sum(abs(bco1))
write(iw,*) 'Check sum = bco2', sum(abs(bco2))
end if
return
end subroutine mrsfcbc