subroutine mo_reorder(infos, Va, Ea, Vb, Eb, Sq)
! Va, Ea: MO and MO energies of Previous Point
! Vb, Eb: MO and MO energies of Current Point
! The Vb and Eb are reordered based on the orders of Va and Ea.
! Sq: Square Form of Overlap Matrix in AO
use precision, only: dp
use io_constants, only: iw
use types, only: information
implicit none
type(information), intent(inout) :: infos
real(kind=dp), intent(inout), dimension(:,:) :: Va, Vb ! (nbf, nbf)
real(kind=dp), intent(inout), dimension(:) :: Ea, Eb ! (nbf)
real(kind=dp), intent(in), dimension(:,:) :: Sq ! (nbf, nbf)
integer :: i, loc
real(kind=dp) :: tmp, tmp2, tmp_abs
real(kind=dp), allocatable, dimension(:,:) :: scr, smo
logical, allocatable, dimension(:) :: pr
integer, allocatable, dimension(:) :: locs
integer :: na, nbf
write(iw, fmt='(/," MOM is on: Reordering the New MO and MO energies....")')
nbf = infos%basis%nbf
na = infos%mol_prop%nelec_a
allocate(scr(nbf,nbf), &
smo(nbf,nbf), &
source=0.0_dp)
allocate(pr(nbf), source=.false.)
allocate(locs(nbf), source=0)
call dgemm('t', 'n', nbf, nbf, nbf, &
1.0_dp, Va, nbf, &
Sq, nbf, &
0.0_dp, scr, nbf)
call dgemm('n', 'n', nbf, nbf, nbf, &
1.0_dp, scr, nbf, &
Vb, nbf, &
0.0_dp, smo, nbf)
! Now Smo is overlap matrix in MO, row and column corresponds to Old and New MOs.
do i = 1, nbf
smo(:,i) = smo(:,i) / norm2(smo(:,i))
end do
! Finding out the location of maximum value in column i
do i = 1, nbf
loc = maxloc(abs(smo(:nbf,i)), dim=1)
locs(i) = loc
pr(loc) = .true.
enddo
! Checking printouts
if (.not.all(pr)) then
write(iw, fmt='(/," Warning")')
write(iw, fmt='(" Some orbitals are missing in reordering")')
write(iw, advance='no', fmt='(" Their indices are:")')
do i = 1, nbf
if (.not.pr(i)) write(iw, advance='no', fmt='(I4,",")') i
end do
write(iw,*)
write(iw, fmt='(/," Error. Stop")')
write(iw, fmt='(" Some orbitals are missing in reordering")')
write(iw,*)
end if
write(iw,fmt='(x,a)') 'Old MO Index ---> New MO Index'
do i = 1, nbf
tmp_abs = maxval(abs(smo(:nbf,i)), dim=1)
loc = maxloc(abs(smo(:nbf,i)), dim=1)
tmp = smo(loc,i)
tmp2 = smo(i,i)
if ((loc.ne.i).and.(i.le.na+1)) then
write(iw, advance='no', &
fmt='(5x,i4,14x,i4)') i, loc
if (i == na-1) write(iw, advance='no',fmt='(2x,a)') ' HOMO'
if (i == na) write(iw, advance='no',fmt='(2x,a)') ' LUMO'
if (i /= loc .and. tmp_abs < 0.9d+00) then
write(iw, fmt='(2x,a)') &
' rearranged, WARNING'
elseif (i == loc .and. tmp_abs < 0.9d+00) then
write(iw, fmt='(2x,a)') &
' WARNING'
elseif (i /= loc .and. tmp_abs > 0.9d+00) then
write(iw, fmt='(2x,a)') &
' rearranged'
else
write(iw,*)
end if
endif
enddo
! Reordering MO and MO energies of Current point.
call reorderMOs(Vb, Eb, Smo, nbf, nbf, 1, na+1)
end subroutine mo_reorder