subroutine sfrowcal(wmo, target_energy, mo_energy_a, fa, fb, bvec, xk, &
xhxa, xhxb, hppija, hppijb, noca, nocb)
use precision, only: dp
implicit none
real(kind=dp), intent(out), dimension(:,:) :: wmo
real(kind=dp), intent(in) :: target_energy
real(kind=dp), intent(in), dimension(:) :: mo_energy_a
real(kind=dp), intent(in), dimension(:,:) :: fa, fb
real(kind=dp), intent(in), dimension(:) :: bvec
real(kind=dp), intent(in), dimension(:) :: xk
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(:,:) :: wrk, wrk1, wrk2
integer :: i, a, k, x, y, j, b, ij, nbf, nvirb, lr1, lr2
nbf = ubound(fa, 1)
lr1 = nocb+1
lr2 = noca
nvirb = nbf-nocb
allocate(wrk(nbf,nbf), &
wrk1(nbf,nbf), &
wrk2(nbf,nbf), &
source=0.0_dp)
! ----- COPY xk -----
ij = 0
do i = nocb+1, noca
do j = 1, nocb
ij = ij+1
wrk1(j,i) = xk(ij)
end do
end do
do i = noca+1, nbf
do j = 1, nocb
ij = ij+1
wrk1(j,i) = xk(ij)
end do
end do
do k = noca+1, nbf
do i = nocb+1, noca
ij = ij+1
wrk1(i,k) = xk(ij)
end do
end do
! ! W_ix
wrk = 0.0_dp
do x = 1, nocb
do k = 1, nocb
wrk(x,1) = wrk(x,1)-wrk1(k,lr1)*fa(k,x)
wrk(x,2) = wrk(x,2)-wrk1(k,lr2)*fa(k,x)
end do
end do
do x = 1, nocb
do k = 1, nbf-noca
wrk(x,1) = wrk(x,1)+wrk1(lr1,noca+k)*fa(noca+k,x)
wrk(x,2) = wrk(x,2)+wrk1(lr2,noca+k)*fa(noca+k,x)
end do
end do
wmo(1:nocb,lr1:lr2) = wrk(1:nocb,1:2)*0.5_dp &
+ xhxa(1:nocb,lr1:lr2) &
+ xhxb(1:nocb,lr1:lr2) &
+ hppija(1:nocb,lr1:lr2)
wmo(1:nocb,lr1) = wmo(1:nocb,lr1) &
+ mo_energy_a(1:nocb)*wrk1(1:nocb,lr1)
wmo(1:nocb,lr2) = wmo(1:nocb,lr2) &
+ mo_energy_a(1:nocb)*wrk1(1:nocb,lr2)
! ----- W_IA -----
wrk = 0.0_dp
do i = 1, nocb
do a = 1, nbf-noca
wrk(i,a) = wrk(i,a)+fa(lr1,i)*wrk1(lr1,noca+a)
wrk(i,a) = wrk(i,a)+fa(lr2,i)*wrk1(lr2,noca+a)
end do
end do
wmo(1:nocb,noca+1:nbf) = wrk(1:nocb,1:nbf-noca)*0.5_dp &
+ xhxb(1:nocb,noca+1:nbf)
do a = noca+1, nbf
wmo(1:nocb,a) = wmo(1:nocb,a) &
+ mo_energy_a(1:nocb)*wrk1(1:nocb,a)
end do
! ----- W_XA -----
wrk = 0.0_dp
do a = 1, nbf-noca
do k = 1, nocb
wrk(1,a) = wrk(1,a)+fa(k,lr1)*wrk1(k,noca+a)
wrk(2,a) = wrk(2,a)+fa(k,lr2)*wrk1(k,noca+a)
end do
end do
do a = 1, nbf-noca
wrk(1,a) = wrk(1,a)-fb(lr1,lr1)*wrk1(lr1,noca+a) &
-fb(lr2,lr1)*wrk1(lr2,noca+a)
wrk(2,a) = wrk(2,a)-fb(lr1,lr2)*wrk1(lr1,noca+a) &
-fb(lr2,lr2)*wrk1(lr2,noca+a)
end do
wmo(lr1:lr2,noca+1:nbf) = wrk(1:2,1:nbf-noca)*0.5_dp &
+ xhxb(lr1:lr2,noca+1:nbf)
do a = noca+1, nbf
wmo(lr1:lr2,a) = wmo(lr1:lr2,a) &
+ mo_energy_a(lr1:lr2)*wrk1(lr1:lr2,a)
end do
! Alpha intermediate
wrk = - fb
do i = nocb+1, nbf
wrk(i,i) = wrk(i,i)+target_energy
end do
wrk1(1:nbf-nocb,1:nbf-nocb) = wrk(nocb+1:nbf,nocb+1:nbf)*2.0_dp
call dgemm('n', 'n', noca, nvirb, nvirb, &
1.0_dp, bvec, noca, &
wrk1, nbf, &
0.0_dp, wrk2, noca)
call dgemm('n', 't', noca, noca, nvirb, &
1.0_dp, wrk2, noca, &
bvec, noca, &
0.0_dp, wrk1, nbf)
! beta intermediate
wrk = fa
do i = 1, noca
wrk(i,i) = wrk(i,i)+target_energy
end do
wrk(1:noca,1:noca) = wrk(1:noca,1:noca)*2.0_dp
call dgemm('n', 'n', noca, nvirb, noca, &
1.0_dp, wrk, nbf, &
bvec, noca, &
0.0_dp, wrk2, noca)
call dgemm('t', 'n', nvirb, nvirb, noca, &
1.0_dp, bvec, noca, &
wrk2, noca, &
0.0_dp, wrk, nbf)
! W_ij
do i = 1, nocb
do j = 1, i
wmo(i,j) = hppija(i,j)+hppijb(i,j)+wrk1(i,j)
end do
end do
! W_xy
do x = nocb+1, noca
do y = nocb+1, x
wmo(x,y) = hppija(x,y)+wrk1(x,y)+wrk(x-nocb,y-nocb)
end do
end do
! W_ab
do a = noca+1, nbf
do b = noca+1, a
wmo(a,b) = wrk(a-nocb,b-nocb)
end do
end do
! Scale diagonal elements
do i = 1, nbf
wmo(i,i) = wmo(i,i)*0.5_dp
end do
wmo = -wmo
deallocate(wrk, wrk1, wrk2)
end subroutine sfrowcal