mrsfmntoia Subroutine

   subroutine mrsfmntoia(infos, fmrsf, pmo, va, vb, ivec)

    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), target, dimension(:,:,:) :: &
      fmrsf
    real(kind=dp), intent(out), dimension(:,:) :: pmo
    real(kind=dp), intent(in), dimension(:,:) :: va, vb
    integer, intent(in) :: ivec

    real(kind=dp), allocatable :: &
      scr(:,:), tmp(:), wrk(:,:)
    real(kind=dp), pointer, dimension(:,:) :: &
      adco1, adco2, ado1v, ado2v, agdlr, aco12, ao21v
    integer :: noca, nocb, mrst, i, ij, &
      j, lr1, lr2, nbf, ok
    real(kind=dp), parameter :: zero = 0.0_dp
    real(kind=dp), parameter :: one = 1.0_dp
    real(kind=dp), parameter :: sqrt2 = 1.0_dp/sqrt(2.0_dp)
    logical :: debug_mode

    nbf = infos%basis%nbf
    mrst = infos%tddft%mult
    noca = infos%mol_prop%nelec_a
    nocb = infos%mol_prop%nelec_b
    debug_mode = infos%tddft%debug_mode

    allocate(tmp(nbf), scr(nbf,nbf), wrk(nbf,nbf), source=0.0_dp, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', with_abort)

    agdlr => fmrsf(7,:,:)
    ado2v => fmrsf(1,:,:)
    ado1v => fmrsf(2,:,:)
    adco1 => fmrsf(3,:,:)
    adco2 => fmrsf(4,:,:)
    ao21v => fmrsf(5,:,:)
    aco12 => fmrsf(6,:,:)

    lr1 = noca-1
    lr2 = noca

    call dgemm('t','n',nbf,nbf,nbf, &
               one, va, nbf, &
                    agdlr,nbf, &
               zero, wrk, nbf)
    call dgemm('n','n',nbf,nbf,nbf, &
               one, wrk, nbf, &
                    vb, nbf, &
               zero, scr, nbf)

! 1
!   ----- (m,n) to (i+,n) -----
    wrk = scr

! 3
    call dgemv('n',nbf,nbf,one,ado1v,nbf,vb(:,lr2),1,zero,tmp,1)
    call dgemv('n',nbf,nbf,one,aco12,nbf,vb(:,lr1),1,one,tmp,1)
    do i = 1, noca-2
      wrk(i,lr2) = wrk(i,lr2) + dot_product(va(:,i),tmp(:))
    end do
! 4
    call dgemv('n',nbf,nbf,one,aco12,nbf,vb(:,lr2),1,zero,tmp,1)
    call dgemv('n',nbf,nbf,one,ado2v,nbf,vb(:,lr1),1,-one,tmp,1)
    do i = 1, noca-2
      wrk(i,lr1) = wrk(i,lr1) + dot_product(va(:,i),tmp(:))
    end do
! 5
    call dgemv('t',nbf,nbf,one,adco2,nbf,va(:,lr1),1,zero,tmp,1)
    call dgemv('t',nbf,nbf,one,ao21v,nbf,va(:,lr2),1, one,tmp,1)
    do j = noca+1, nbf
      wrk(lr1,j) = wrk(lr1,j) + dot_product(vb(:,j),tmp(:))
    end do
! 6
    call dgemv('t',nbf,nbf,one,ao21v,nbf,va(:,lr1),1,zero,tmp,1)
    call dgemv('t',nbf,nbf,one,adco1,nbf,va(:,lr2),1,-one,tmp,1)
    do j = noca+1, nbf
       wrk(lr2,j) = wrk(lr2,j) + dot_product(vb(:,j),tmp(:))
    end do

    if (mrst==1) then
      wrk(lr1,lr1) = (scr(lr1,lr1)-scr(lr2,lr2))*sqrt2
      wrk(lr2,lr2) = 0.0_dp
    else if (mrst==3) then
      wrk(lr1,lr1) = (scr(lr1,lr1)+scr(lr2,lr2))*sqrt2
      wrk(lr2,lr1) = 0.0_dp
      wrk(lr1,lr2) = 0.0_dp
      wrk(lr2,lr2) = 0.0_dp
    end if

    pmo(:,ivec) = 0.0_dp

    ij = 0
    do j = nocb+1, nbf
      do i = 1, noca
        ij = ij+1
        pmo(ij,ivec) = pmo(ij,ivec) + wrk(i,j)
      end do
    end do

    if (debug_mode) then
      write(iw,*) 'Check sum = ivec', ivec
      write(iw,*) 'Check sum = agdlr', sum(abs(agdlr))
      write(iw,*) 'Check sum = ao21v', sum(abs(ao21v))
      write(iw,*) 'Check sum = aco12', sum(abs(aco12))
      write(iw,*) 'Check sum = ado2v', sum(abs(ado2v))
      write(iw,*) 'Check sum = ado1v', sum(abs(ado1v))
      write(iw,*) 'Check sum = adco1', sum(abs(adco1))
      write(iw,*) 'Check sum = adco2', sum(abs(adco2))
      write(iw,*) 'Check sum = pmo', sum(abs(pmo(:,ivec)))
    end if

    return

  end subroutine mrsfmntoia