subroutine mrsfqrowcal(w, mo_energy_a, fa, fb, z, &
xhxa, xhxb, hppija, hppijb, noca, nocb)
use precision, only: dp
implicit none
real(kind=dp), intent(out), dimension(:,:) :: w
real(kind=dp), intent(in), dimension(:) :: mo_energy_a
real(kind=dp), intent(in), dimension(:,:) :: fa, fb
real(kind=dp), intent(in), dimension(:) :: z
real(kind=dp), intent(in), dimension(:,:) :: xhxa, xhxb
real(kind=dp), intent(in), dimension(:,:) :: hppija, hppijb
integer, intent(in) :: noca, nocb
real(kind=dp), allocatable, dimension(:,:) :: scr, wrk
integer :: i, a, k, x, y, j, b, ij, nbf, lr1, lr2
nbf = ubound(fa, 1)
lr1 = nocb+1
lr2 = noca
allocate(wrk(nbf,nbf), &
scr(nbf,nbf), &
source=0.0_dp)
ij = 0
do x = nocb+1, noca
do i = 1, nocb
ij = ij+1
scr(i, x) = z(ij)
end do
end do
do a = noca+1, nbf
do i = 1, nocb
ij = ij+1
scr(i, a) = z(ij)
end do
end do
do a = noca+1, nbf
do x = nocb+1, noca
ij = ij+1
scr(x, a) = z(ij)
end do
end do
! w_ix
do i = 1, nocb
do k = 1, nocb
wrk(i,1:2) = wrk(i,1:2)-fa(k,i)*scr(k,lr1:lr2)
end do
end do
do i = 1, nocb
do y = 1, nbf-noca
wrk(i,1:2) = wrk(i,1:2)+fa(noca+y,i)*scr(lr1:lr2,noca+y)
end do
end do
w(1:nocb,lr1:lr2) = 0.5_dp*wrk(1:nocb,1:2)+hppija(1:nocb,lr1:lr2)
do i = 1,nocb
w(i,lr1:lr2) = w(i,lr1:lr2) &
+ mo_energy_a(i)*scr(i,lr1:lr2)
end do
! w_ia
wrk = 0.0_dp
do a = 1, nbf-noca
do i = 1, nocb
wrk(i,a) = wrk(i,a)+fa(lr1,i)*scr(lr1,noca+a)
wrk(i,a) = wrk(i,a)+fa(lr2,i)*scr(lr2,noca+a)
end do
end do
do a = 1, nbf-noca
w(1:nocb,noca+a) = mo_energy_a(1:nocb)*scr(1:nocb,noca+a) &
+ 0.5_dp*wrk(1:nocb,a) &
+ xhxa(1:nocb,noca+a)
end do
! w_xa
wrk = 0.0_dp
do a = 1, nbf-noca
do k = 1, nocb
wrk(1:2,a) = wrk(1:2,a)+fa(k, lr1:lr2)*scr(k, noca+a)
end do
end do
do a = 1, nbf-noca
wrk(1:2,a) = wrk(1:2,a)-fb(lr1, lr1:lr2)*scr(lr1, noca+a)
wrk(1:2,a) = wrk(1:2,a)-fb(lr2, lr1:lr2)*scr(lr2, noca+a)
end do
do a = 1, nbf-noca
w(lr1:lr2, noca+a) = mo_energy_a(lr1:lr2)*scr(lr1:lr2, noca+a) &
+ 0.5_dp*wrk(1:2,a) &
+ xhxa(lr1:lr2, noca+a)
end do
! w_ij
do i = 1, nocb
do j = 1, i
w(i, j) = hppija(i, j)+hppijb(i, j)+xhxb(j, i)
end do
end do
! w_xy
do x = nocb+1, noca
do y = nocb+1, x
w(x, y) = hppija(x, y)
end do
end do
! w_ab
do a = noca+1, nbf
do b = noca+1, a
w(a, b) = xhxa(b, a)
end do
end do
! scale diagonal elements
do i = 1, nbf
w(i, i) = 0.5_dp*w(i, i)
end do
w = -w
return
end subroutine mrsfqrowcal