This subroutine selects the states that are inside the outer window (ie, the energy window out of which we fish out the optimally-connected subspace) and those that are inside the inner window (that make up the frozen manifold, and are straightfowardly included as they are). This, in practice, amounts to slimming down the original num_wann x num_wann overlap matrices, removing rows and columns that belong to u_nks that have been excluded forever, and marking (indexing) the rows and columns that correspond to frozen states.
Note - in windows eigval_opt are shifted, so the lowest ones go from nfirstwin(nkp) to nfirstwin(nkp)+ndimwin(nkp)-1, and above they are set to zero.
subroutine dis_windows()
!==================================================================!
! !
!! This subroutine selects the states that are inside the outer
!! window (ie, the energy window out of which we fish out the
!! optimally-connected subspace) and those that are inside the
!! inner window (that make up the frozen manifold, and are
!! straightfowardly included as they are). This, in practice,
!! amounts to slimming down the original num_wann x num_wann overlap
!! matrices, removing rows and columns that belong to u_nks that
!! have been excluded forever, and marking (indexing) the rows and
!! columns that correspond to frozen states.
!!
!! Note - in windows eigval_opt are shifted, so the lowest ones go
!! from nfirstwin(nkp) to nfirstwin(nkp)+ndimwin(nkp)-1, and above
!! they are set to zero.
! !
!==================================================================!
implicit none
! internal variables
integer :: i, j, nkp, ierr
integer :: imin, imax, kifroz_min, kifroz_max
!~~ GS-start
real(kind=dp) :: dk(3), kdr2
logical :: dis_ok
!~~ GS-end
! OUTPUT:
! ndimwin(nkp) number of bands inside outer window at nkp-th k poi
! ndimfroz(nkp) number of frozen bands at nkp-th k point
! lfrozen(i,nkp) true if the i-th band inside outer window is frozen
! linner true if there is an inner window
! indxfroz(i,nkp) outer-window band index for the i-th frozen state
! (equals 1 if it is the bottom of outer window)
! indxnfroz(i,nkp) outer-window band index for the i-th non-frozen s
! (equals 1 if it is the bottom of outer window)
! nfirstwin(nkp) index of lowest band inside outer window at nkp-th
! MODIFIED:
! eigval_opt(nb,nkp) At input it contains a large set of eigenvalues. At
! it is slimmed down to contain only those inside the
! energy window, stored in nb=1,...,ndimwin(nkp)
if (timing_level > 1 .and. on_root) call io_stopwatch('dis: windows', 1)
! Allocate module arrays
allocate (nfirstwin(num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating nfirstwin in dis_windows')
allocate (ndimfroz(num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating ndimfroz in dis_windows')
allocate (indxfroz(num_bands, num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating indxfroz in dis_windows')
allocate (indxnfroz(num_bands, num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating indxnfroz in dis_windows')
allocate (lfrozen(num_bands, num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating lfrozen in dis_windows')
linner = .false.
if (on_root) write (stdout, '(1x,a)') &
'+----------------------------------------------------------------------------+'
if (on_root) write (stdout, '(1x,a)') &
'| Energy Windows |'
if (on_root) write (stdout, '(1x,a)') &
'| --------------- |'
if (on_root) write (stdout, '(1x,a,f10.5,a,f10.5,a)') &
'| Outer: ', dis_win_min, ' to ', dis_win_max, &
' (eV) |'
if (frozen_states) then
if (on_root) write (stdout, '(1x,a,f10.5,a,f10.5,a)') &
'| Inner: ', dis_froz_min, ' to ', dis_froz_max, &
' (eV) |'
else
if (on_root) write (stdout, '(1x,a)') &
'| No frozen states were specified |'
endif
if (on_root) write (stdout, '(1x,a)') &
'+----------------------------------------------------------------------------+'
do nkp = 1, num_kpts
! Check which eigenvalues fall within the outer window
if ((eigval_opt(1, nkp) .gt. dis_win_max) .or. &
(eigval_opt(num_bands, nkp) .lt. dis_win_min)) then
if (on_root) write (stdout, *) ' ERROR AT K-POINT: ', nkp
if (on_root) write (stdout, *) ' ENERGY WINDOW (eV): [', dis_win_min, ',', dis_win_max, ']'
if (on_root) write (stdout, *) ' EIGENVALUE RANGE (eV): [', &
eigval_opt(1, nkp), ',', eigval_opt(num_bands, nkp), ']'
call io_error('dis_windows: The outer energy window contains no eigenvalues')
endif
! Note: we assume that eigvals are ordered from the bottom up
imin = 0
do i = 1, num_bands
if (imin .eq. 0) then
if ((eigval_opt(i, nkp) .ge. dis_win_min) .and. &
(eigval_opt(i, nkp) .le. dis_win_max)) imin = i
imax = i
endif
if (eigval_opt(i, nkp) .le. dis_win_max) imax = i
enddo
ndimwin(nkp) = imax - imin + 1
nfirstwin(nkp) = imin
!~~ GS-start
! disentangle at the current k-point only if it is within one of the
! spheres centered at the k-points listed in kpt_dis
if (dis_spheres_num .gt. 0) then
dis_ok = .false.
! loop on the sphere centers
do i = 1, dis_spheres_num
dk = kpt_latt(:, nkp) - dis_spheres(1:3, i)
dk = matmul(anint(dk) - dk, recip_lattice(:, :))
! if the current k-point is included in at least one sphere,
! then perform disentanglement as usual
if (abs(dot_product(dk, dk)) .lt. dis_spheres(4, i)**2) then
dis_ok = .true.
exit
endif
enddo
! this kpoint is not included in any sphere: no disentaglement
if (.not. dis_ok) then
ndimwin(nkp) = num_wann
nfirstwin(nkp) = dis_spheres_first_wann
endif
endif
!~~ GS-end
if (ndimwin(nkp) .lt. num_wann) then
if (on_root) write (stdout, 483) 'Error at k-point ', nkp, &
' ndimwin=', ndimwin(nkp), ' num_wann=', num_wann
483 format(1x, a17, i4, a8, i3, a9, i3)
call io_error('dis_windows: Energy window contains fewer states than number of target WFs')
endif
do i = 1, ndimwin(nkp)
lfrozen(i, nkp) = .false.
enddo
! Check which eigenvalues fall within the inner window
kifroz_min = 0
kifroz_max = -1
! (note that the above obeys kifroz_max-kifroz_min+1=kdimfroz=0, as we w
if (frozen_states) then
do i = imin, imax
if (kifroz_min .eq. 0) then
if ((eigval_opt(i, nkp) .ge. dis_froz_min) .and. &
(eigval_opt(i, nkp) .le. dis_froz_max)) then
! Relative to bottom of outer window
kifroz_min = i - imin + 1
kifroz_max = i - imin + 1
endif
elseif (eigval_opt(i, nkp) .le. dis_froz_max) then
kifroz_max = kifroz_max + 1
! DEBUG
! if(kifroz_max.ne.i-imin+1) stop 'something wrong...'
! ENDDEBUG
endif
enddo
endif
ndimfroz(nkp) = kifroz_max - kifroz_min + 1
if (ndimfroz(nkp) .gt. num_wann) then
if (on_root) write (stdout, 401) nkp, ndimfroz(nkp), num_wann
401 format(' ERROR AT K-POINT ', i4, ' THERE ARE ', i2, &
' BANDS INSIDE THE INNER WINDOW AND ONLY', i2, &
' TARGET BANDS')
if (on_root) write (stdout, 402) (eigval_opt(i, nkp), i=imin, imax)
402 format('BANDS: (eV)', 10(F10.5, 1X))
call io_error('dis_windows: More states in the frozen window than target WFs')
endif
if (ndimfroz(nkp) .gt. 0) linner = .true.
! DEBUG
! write(*,'(a,i4,a,i2,a,i2)') 'k point ',nkp,
! & ' lowest band in outer win is # ',imin,
! & ' # frozen states is ',ndimfroz(nkp)
! ENDDEBUG
! Generate index array for frozen states (those inside inner window)
if (ndimfroz(nkp) .gt. 0) then
do i = 1, ndimfroz(nkp)
indxfroz(i, nkp) = kifroz_min + i - 1
lfrozen(indxfroz(i, nkp), nkp) = .true.
enddo
if (indxfroz(ndimfroz(nkp), nkp) .ne. kifroz_max) then
if (on_root) write (stdout, *) ' Error at k-point ', nkp, ' frozen band #', i
if (on_root) write (stdout, *) ' ndimfroz=', ndimfroz(nkp)
if (on_root) write (stdout, *) ' kifroz_min=', kifroz_min
if (on_root) write (stdout, *) ' kifroz_max=', kifroz_max
if (on_root) write (stdout, *) ' indxfroz(i,nkp)=', indxfroz(i, nkp)
call io_error('dis_windows: Something fishy...')
endif
endif
! Generate index array for non-frozen states
i = 0
do j = 1, ndimwin(nkp)
! if (lfrozen(j,nkp).eqv..false.) then
if (.not. lfrozen(j, nkp)) then
i = i + 1
indxnfroz(i, nkp) = j
endif
enddo
if (i .ne. ndimwin(nkp) - ndimfroz(nkp)) then
if (on_root) write (stdout, *) ' Error at k-point: ', nkp
if (on_root) write (stdout, '(3(a,i5))') ' i: ', i, ' ndimwin: ', ndimwin(nkp), &
' ndimfroz: ', ndimfroz(nkp)
call io_error('dis_windows: i .ne. (ndimwin-ndimfroz) at k-point')
endif
! Slim down eigval vector at present k
do i = 1, ndimwin(nkp)
j = nfirstwin(nkp) + i - 1
eigval_opt(i, nkp) = eigval_opt(j, nkp)
enddo
do i = ndimwin(nkp) + 1, num_bands
eigval_opt(i, nkp) = 0.0_dp
enddo
enddo
! [k-point loop (nkp)]
![ysl-b]
!~ if (gamma_only) then
!~ if (.not. allocated(ph_g)) then
!~ allocate( ph_g(num_bands),stat=ierr )
!~ if (ierr/=0) call io_error('Error in allocating ph_g in dis_windows')
!~ ph_g = cmplx_1
!~ endif
!~ ! Apply same operation to ph_g
!~ do i = 1, ndimwin(1)
!~ j = nfirstwin(1) + i - 1
!~ ph_g(i) = ph_g(j)
!~ enddo
!~ do i = ndimwin(1) + 1, num_bands
!~ ph_g(i) = cmplx_0
!~ enddo
!~ endif
!~![ysl-e]
if (iprint > 1) then
if (on_root) write (stdout, '(1x,a)') &
'| K-points with Frozen States |'
if (on_root) write (stdout, '(1x,a)') &
'| --------------------------- |'
i = 0
do nkp = 1, num_kpts
if (ndimfroz(nkp) .gt. 0) then
i = i + 1
if (i .eq. 1) then
if (on_root) write (stdout, '(1x,a,i6)', advance='no') '|', nkp
else if ((i .gt. 1) .and. (i .lt. 12)) then
if (on_root) write (stdout, '(i6)', advance='no') nkp
else if (i .eq. 12) then
if (on_root) write (stdout, '(i6,a)') nkp, ' |'
i = 0
endif
endif
enddo
if (i .ne. 0) then
do j = 1, 12 - i
if (on_root) write (stdout, '(6x)', advance='no')
enddo
if (on_root) write (stdout, '(a)') ' |'
endif
if (on_root) write (stdout, '(1x,a)') &
'+----------------------------------------------------------------------------+'
endif
if (on_root) write (stdout, '(3x,a,i4)') 'Number of target bands to extract: ', num_wann
if (iprint > 1) then
if (on_root) write (stdout, '(1x,a)') &
'+----------------------------------------------------------------------------+'
if (on_root) write (stdout, '(1x,a)') &
'| Windows |'
if (on_root) write (stdout, '(1x,a)') &
'| ------- |'
if (on_root) write (stdout, '(1x,a)') &
'| K-point Ndimwin Ndimfroz Nfirstwin |'
if (on_root) write (stdout, '(1x,a)') &
'| ---------------------------------------------- |'
do nkp = 1, num_kpts
if (on_root) write (stdout, 403) nkp, ndimwin(nkp), ndimfroz(nkp), nfirstwin(nkp)
enddo
403 format(1x, '|', 14x, i6, 7x, i6, 7x, i6, 6x, i6, 18x, '|')
if (on_root) write (stdout, '(1x,a)') &
'+----------------------------------------------------------------------------+'
endif
if (timing_level > 1) call io_stopwatch('dis: windows', 2)
return
end subroutine dis_windows