dis_windows Subroutine

private 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.

Arguments

None

Calls

proc~~dis_windows~~CallsGraph proc~dis_windows dis_windows 10 10 proc~dis_windows->10 proc~io_error io_error proc~dis_windows->proc~io_error

Called by

proc~~dis_windows~~CalledByGraph proc~dis_windows dis_windows proc~dis_main dis_main proc~dis_main->proc~dis_windows program~wannier wannier program~wannier->proc~dis_main proc~wannier_run wannier_run proc~wannier_run->proc~dis_main

Contents

Source Code


Source Code

  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