dis_proj_froz – Wannier90

public subroutine dis_proj_froz(u_matrix_opt, indxfroz, ndimfroz, ndimwin, iprint, num_bands, num_kpts, num_wann, timing_level, lfrozen, on_root, timer, error, stdout, comm)

Uses

Help

Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.
COMPUTES THE LEADING EIGENVECTORS OF Q_froz . P_s . Q_froz, WHERE P_s PROJECTOR OPERATOR ONTO THE SUBSPACE S OF THE PROJECTED GAUSSIANS, P_f THE PROJECTOR ONTO THE FROZEN STATES, AND Q_froz = 1 - P_froz, ALL EXP IN THE BASIS OF THE BLOCH EIGENSTATES INSIDE THE OUTER ENERGY WINDOW (See Eq. (27) in Sec. III.G of SMV)
Arguments

Type	Intent	Attributes		Name
complex(kind=dp),	intent(inout)		::	u_matrix_opt(:,:,:)
integer,	intent(in)		::	indxfroz(:,:)
integer,	intent(in)		::	ndimfroz(:)
integer,	intent(in)		::	ndimwin(:)
integer,	intent(in)		::	iprint
integer,	intent(in)		::	num_bands
integer,	intent(in)		::	num_kpts
integer,	intent(in)		::	num_wann
integer,	intent(in)		::	timing_level
logical,	intent(in)		::	lfrozen(:,:)
logical,	intent(in)		::	on_root
type(timer_list_type),	intent(inout)		::	timer
type(w90_error_type),	intent(out),	allocatable	::	error
integer,	intent(in)		::	stdout
type(w90_comm_type),	intent(in)		::	comm
Calls

Help
Called by

Help
Source Code

  subroutine dis_proj_froz(u_matrix_opt, indxfroz, ndimfroz, ndimwin, iprint, num_bands, num_kpts, &
                           num_wann, timing_level, lfrozen, on_root, timer, error, stdout, comm)
    !================================================!
    !
    !! COMPUTES THE LEADING EIGENVECTORS OF Q_froz . P_s . Q_froz,
    !! WHERE P_s PROJECTOR OPERATOR ONTO THE SUBSPACE S OF THE PROJECTED
    !! GAUSSIANS, P_f THE PROJECTOR ONTO THE FROZEN STATES, AND
    !! Q_froz = 1 - P_froz, ALL EXP IN THE BASIS OF THE BLOCH
    !! EIGENSTATES INSIDE THE OUTER ENERGY WINDOW
    !! (See Eq. (27) in Sec. III.G of SMV)
    !
    !================================================!
    use w90_constants, only: dp, cmplx_0, cmplx_1, eps8
    use w90_error
    use w90_io, only: io_stopwatch_start, io_stopwatch_stop
    use w90_types, only: timer_list_type

    implicit none

    ! arguments
    type(timer_list_type), intent(inout) :: timer
    type(w90_error_type), allocatable, intent(out) :: error
    type(w90_comm_type), intent(in) :: comm

    integer, intent(in) :: timing_level, iprint
    integer, intent(in) :: stdout
    integer, intent(in) :: num_bands, num_kpts, num_wann
    integer, intent(in) :: ndimwin(:)
    integer, intent(in) :: ndimfroz(:)
    integer, intent(in) :: indxfroz(:, :)

    complex(kind=dp), intent(inout) :: u_matrix_opt(:, :, :)

    logical, intent(in) :: on_root, lfrozen(:, :)

    ! INPUT: num_wann,ndimwin,ndimfroz,indxfroz,lfrozen
    ! MODIFIED: u_matrix_opt (At input it contains the gaussians projected onto
    !             the window states in the routine project.f. At output
    !             the entries with the second index from 1 to ndimfroz(nkp)
    !             contain the frozen (inner window) states, while those
    !             from ndimfroz(nkp)+1 to num_wann have been replaced by
    !             the new trial states outside the inner window.)

    ! *********************************************************
    ! VARIABLES USED BY LAPACK'S ZHPEVX DIAGONALIZATION ROUTINE
    ! *********************************************************
    integer, allocatable :: iwork(:)
    integer, allocatable :: ifail(:)

    real(kind=dp), allocatable :: w(:)
    real(kind=dp), allocatable :: rwork(:)

    complex(kind=dp), allocatable :: cap(:)
    complex(kind=dp), allocatable :: cwork(:)
    complex(kind=dp), allocatable :: cz(:, :)

    ! *********
    ! INTERNAL:
    ! *********
    !
    ! CP_S(M,N)      PROJECTION OPERATOR ONTO THE SUBSPACE OF THE PROJEC
    !                  GAUSSIANS, EXPRESSED IN THE BASIS OF THE ORIGINAL BL
    !                  EIGENSTATES INSIDE THE OUTER WINDOW (FOR THE PRESENT
    !                  K-POINT)
    ! CQ_FROZ(M,N)   PROJECTION OPERATOR ONTO THE SUBSPACE OF THE STATES
    !                  THE SPACE OF FROZEN STATES (BUT INSIDE THE OUTER WIN
    !                  EXPRESSED IN THE BASIS OF THE ORIGINAL BLOCH EIGENST
    !                  INSIDE THE OUTER WINDOW (FOR THE PRESENT K-POINT)
    ! CPQ(M,N)       THE MATRIX cp_s . cq_froz FOR THE PRESENT K-POINT
    ! CQPQ(M,N)      THE MATRIX cq_froz . cp_s . cq_froz FOR THE PRESENT
    !

    integer :: goods, il, iu, nkp, l, j, n, m, info, ierr
    integer :: counter, loop_f, loop_v, vmap(num_bands)
    integer :: nzero

    complex(kind=dp) :: ctmp
    complex(kind=dp), allocatable :: cp_s(:, :)
    complex(kind=dp), allocatable :: cq_froz(:, :)
    complex(kind=dp), allocatable :: cpq(:, :)
    complex(kind=dp), allocatable :: cqpq(:, :)

    logical :: take

    character(len=4) :: rep

    if (timing_level > 1) call io_stopwatch_start('dis: proj_froz', timer)

    if (iprint > 0) write (stdout, '(3x,a)', advance='no') 'In dis_proj_froz...'

    allocate (iwork(5*num_bands), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error allocating iwork in dis_proj_froz', comm)
      return
    end if
    allocate (ifail(num_bands), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error allocating ifail in dis_proj_froz', comm)
      return
    end if
    allocate (w(num_bands), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error allocating w in dis_proj_froz', comm)
      return
    end if
    allocate (rwork(7*num_bands), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error allocating rwork in dis_proj_froz', comm)
      return
    end if
    allocate (cap((num_bands*(num_bands + 1))/2), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error allocating cap in dis_proj_froz', comm)
      return
    end if
    allocate (cwork(2*num_bands), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error allocating cwork in dis_proj_froz', comm)
      return
    end if
    allocate (cz(num_bands, num_bands), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error allocating cz in dis_proj_froz', comm)
      return
    end if
    allocate (cp_s(num_bands, num_bands), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error allocating cp_s in dis_proj_froz', comm)
      return
    end if
    allocate (cq_froz(num_bands, num_bands), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error allocating cq_froz in dis_proj_froz', comm)
      return
    end if
    allocate (cpq(num_bands, num_bands), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error allocating cpq in dis_proj_froz', comm)
      return
    end if
    allocate (cqpq(num_bands, num_bands), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error allocating cqpq in dis_proj_froz', comm)
      return
    end if

    do nkp = 1, num_kpts

      ! aam: this should be done at the end, otherwise important
      !      projection info is lost
!~         ! Put the frozen states in the lowest columns of u_matrix_opt
!~         if (ndimfroz(nkp).gt.0) then
!~            do l = 1, ndimfroz(nkp)
!~               u_matrix_opt(:,l,nkp)=cmplx_0
!~               u_matrix_opt(indxfroz(l,nkp),l,nkp) = cmplx_1
!~            enddo
!~         endif

      ! If there are non-frozen states, compute the num_wann-ndimfroz(nkp) leadin
      ! eigenvectors of cqpq
      if (num_wann .gt. ndimfroz(nkp)) then
        cq_froz = cmplx_0
        cp_s = cmplx_0
        do n = 1, ndimwin(nkp)
          do m = 1, ndimwin(nkp)
            do l = 1, num_wann
              cp_s(m, n) = cp_s(m, n) + u_matrix_opt(m, l, nkp)*conjg(u_matrix_opt(n, l, nkp))
            end do
          end do
          if (.not. lfrozen(n, nkp)) cq_froz(n, n) = cmplx_1
        end do

        cpq = cmplx_0
        do n = 1, ndimwin(nkp)
          do m = 1, ndimwin(nkp)
            do l = 1, ndimwin(nkp)
              cpq(m, n) = cpq(m, n) + cp_s(m, l)*cq_froz(l, n)
            end do
          end do
        end do

        cqpq = cmplx_0
        do n = 1, ndimwin(nkp)
          do m = 1, ndimwin(nkp)
            do l = 1, ndimwin(nkp)
              cqpq(m, n) = cqpq(m, n) + cq_froz(m, l)*cpq(l, n)
            end do
          end do
        end do

        ! DEBUG
        ! check hermiticity of cqpq
        do n = 1, ndimwin(nkp)
          do m = 1, n
            if (abs(cqpq(m, n) - conjg(cqpq(n, m))) .gt. eps8) then
              if (on_root) write (stdout, *) ' matrix CQPQ is not hermitian'
              if (on_root) write (stdout, *) ' k-point ', nkp
              call set_error_fatal(error, 'dis_proj_froz: error', comm)
              return
            end if
          end do
        end do
        ! ENDDEBUG

        cap = cmplx_0
        do n = 1, ndimwin(nkp)
          do m = 1, n
            cap(m + (n - 1)*n/2) = cqpq(m, n)
          end do
        end do
        il = ndimwin(nkp) - (num_wann - ndimfroz(nkp)) + 1
        iu = ndimwin(nkp)
        call zhpevx('V', 'A', 'U', ndimwin(nkp), cap, 0.0_dp, 0.0_dp, il, iu, -1.0_dp, m, w, cz, &
                    num_bands, cwork, rwork, iwork, ifail, info)

!~            write(stdout,*) 'w:'
!~            do n=1,ndimwin(nkp)
!~               write(stdout,'(f14.10)') w(n)
!~            enddo
!~            write(stdout,*) 'cz:'
!~            do n=1,ndimwin(nkp)
!~               write(stdout,'(6f12.8)') cz(n,il), cz(n,iu)
!~            enddo

        ! DEBUG
        if (info .lt. 0) then
          if (on_root) write (stdout, *) ' *** ERROR *** ZHPEVX WHILE DIAGONALIZING CQPQ MATRIX'
          if (on_root) write (stdout, *) ' THE ', -info, ' ARGUMENT OF ZHPEVX HAD AN ILLEGAL VALUE'
          call set_error_fatal(error, 'dis_proj_froz: error', comm)
          return
        elseif (info .gt. 0) then
          if (on_root) write (stdout, *) ' *** ERROR *** ZHPEVX WHILE DIAGONALIZING CQPQ MATRIX'
          if (on_root) write (stdout, *) info, 'EIGENVECTORS FAILED TO CONVERGE'
          call set_error_fatal(error, 'dis_proj_froz: error', comm)
          return
        end if
        ! ENDDEBUG

        ! DEBUG
        if (m .ne. ndimwin(nkp)) then
          if (on_root) write (stdout, *) ' *** ERROR *** in dis_proj_froz'
          if (on_root) write (stdout, *) ' Number of eigenvalues/vectors obtained is', &
            m, ' not equal to the number asked,', ndimwin(nkp)
          call set_error_fatal(error, 'dis_proj_froz: error', comm)
          return
        end if
        ! ENDDEBUG

        ! DEBUG
        ! check that the eigenvalues are between 0 and 1
        if (iprint > 2) then
          if (on_root) write (stdout, '(/a,i3,a,i3,a,i3,a)') ' K-point ', nkp, ' ndimwin: ', &
            ndimwin(nkp), ' we want the', num_wann - ndimfroz(nkp), &
            ' leading eigenvector(s) of QPQ'
        end if
        do j = 1, ndimwin(nkp)
          if (iprint > 2 .and. on_root) write (stdout, '(a,i3,a,f16.12)') '  lambda(', j, ')=', w(j)
!~[aam]        if ( (w(j).lt.eps8).or.(w(j).gt.1.0_dp + eps8) ) then
          if ((w(j) .lt. -eps8) .or. (w(j) .gt. 1.0_dp + eps8)) then
            call set_error_fatal(error, 'dis_proj_froz: error - Eigenvalues not between 0 and 1', comm)
            return
          end if
        end do
        ! ENDDEBUG

        ! [ aam: sometimes the leading eigenvalues form a degenerate set that is
        !        of higher dimensionality than (num_wann - ndimfroz). May need to
        !        fix this at some point. ]

        ! For certain cases we have found that one of the required eigenvectors of cqpq
        ! has a zero eigenvalue (ie it forms a degenerate set with the frozen states).
        ! It depends on floating point math as whether this eigenvalue is positive
        ! or negative (ie +/- 1e-17). If it's positive everything is ok. If negative we
        ! can end up putting one of the frozen states into U_opt (and failing the later
        ! orthogonality check).
        ! This fix detects this situation. If applies we choose the eigenvectors by
        ! checking their orthogonality to the frozen states.
        ! === For version 1.0.1 we make this the default ===

        !if (index(devel_flag, 'no-orth-fix') == 0) then
        nzero = 0; goods = 0
        do j = ndimwin(nkp), ndimwin(nkp) - (num_wann - ndimfroz(nkp)) + 1, -1
          if (w(j) < eps8) then
            nzero = nzero + 1
          else
            goods = goods + 1
          end if
        end do
        if (nzero > 0) then
          if (iprint > 2 .and. on_root) then
            write (stdout, *) ' '
            write (stdout, '(1x,a,i0,a)') 'An eigenvalue of QPQ is close to zero at kpoint ', &
              nkp, '. Using safety check.'
            write (stdout, '(1x,a,i4,a,i4)') 'We must find ', nzero, &
              ' eigenvectors with zero eigenvalues out of a set of ', ndimwin(nkp) - goods
          end if
          !First lets put the 'good' states into vamp
          vmap = 0
          counter = 1
          do j = ndimwin(nkp), ndimwin(nkp) - goods + 1, -1
            vmap(counter) = j
            counter = counter + 1
          end do

          if (iprint > 2 .and. on_root) then
            do loop_f = 1, ndimwin(nkp)
              write (stdout, '(1x,a,i4,a,es13.6)') 'Eigenvector number', loop_f, &
                '    Eigenvalue: ', w(loop_f)
              do loop_v = 1, ndimwin(nkp)
                write (stdout, '(20x,2f12.8)') cz(loop_v, loop_f)
              end do
              write (stdout, *)
            end do
          end if

          ! We need to find nzero vectors out of the remining ndimwin(nkp)-goods vectors

          do loop_f = 1, nzero
            do loop_v = ndimwin(nkp), 1, -1 !loop backwards for efficiency only
              if (any(vmap == loop_v)) cycle
              !check to see if vector is orthogonal to frozen states in u_matrix_opt
              take = .true.
              do m = 1, ndimfroz(nkp)
                ctmp = cmplx_0
                do j = 1, ndimwin(nkp)
                  ctmp = ctmp + conjg(u_matrix_opt(j, m, nkp))*cz(j, loop_v)
                end do
                if (abs(ctmp) .gt. eps8) then
                  take = .false.
                end if
              end do
              if (take) then !vector is good and we add it to vmap
                vmap(goods + loop_f) = loop_v
                exit
              end if
            end do
          end do

          if (iprint > 2 .and. on_root) then
            write (rep, '(i4)') num_wann - ndimfroz(nkp)
            write (stdout, '(1x,a,'//trim(rep)//'(i0,1x))') &
              'We use the following eigenvectors: ', vmap(1:(num_wann - ndimfroz(nkp)))
          end if
          do l = 1, num_wann - ndimfroz(nkp)
            if (vmap(l) == 0) then
              call set_error_fatal(error, 'dis_proj_froz: Ortho-fix failed to find enough vectors', comm)
              return
            end if
          end do

          ! put the correct eigenvectors into u_matrix_opt, and we're all done!
          counter = 1
          do l = ndimfroz(nkp) + 1, num_wann
            u_matrix_opt(1:ndimwin(nkp), l, nkp) = cz(1:ndimwin(nkp), vmap(counter))
            counter = counter + 1
          end do

        else ! we don't need to use the fix

          do l = ndimfroz(nkp) + 1, num_wann
            u_matrix_opt(1:ndimwin(nkp), l, nkp) = cz(1:ndimwin(nkp), il)
            il = il + 1
          end do

          if (il - 1 .ne. iu) then
            call set_error_fatal(error, 'dis_proj_froz: error -  il-1.ne.iu  (in ortho-fix)', comm)
            return
          end if

        end if

        ! PICK THE num_wann-nDIMFROZ(NKP) LEADING EIGENVECTORS AS TRIAL STATES
        ! and PUT THEM RIGHT AFTER THE FROZEN STATES IN u_matrix_opt
        !do l = ndimfroz(nkp) + 1, num_wann
        !  if (on_root) write (stdout, *) 'il=', il
        !  u_matrix_opt(1:ndimwin(nkp), l, nkp) = cz(1:ndimwin(nkp), il)
        !  il = il + 1
        !enddo

        ! DEBUG
        !if (il - 1 .ne. iu) then
        !  call io_error('dis_proj_frozen: error -  il-1.ne.iu', stdout, seedname)
        !endif
        ! ENDDEBUG

        !end if

      end if   ! num_wann>nDIMFROZ(NKP)

      ! Put the frozen states in the lowest columns of u_matrix_opt
      if (ndimfroz(nkp) .gt. 0) then
        do l = 1, ndimfroz(nkp)
          u_matrix_opt(:, l, nkp) = cmplx_0
          u_matrix_opt(indxfroz(l, nkp), l, nkp) = cmplx_1
        end do
      end if
    end do   ! NKP

    deallocate (cqpq, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error deallocating cqpq in dis_proj_froz', comm)
      return
    end if
    deallocate (cpq, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error deallocating cpq in dis_proj_froz', comm)
      return
    end if
    deallocate (cq_froz, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error deallocating cq_froz in dis_proj_froz', comm)
      return
    end if
    deallocate (cp_s, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error deallocating cp_s in dis_proj_froz', comm)
      return
    end if
    deallocate (cz, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error deallocating cz in dis_proj_froz', comm)
      return
    end if
    deallocate (cwork, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error deallocating cwork in dis_proj_froz', comm)
      return
    end if
    deallocate (cap, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating cap in dis_proj_froz', comm)
      return
    end if
    deallocate (rwork, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error deallocating rwork in dis_proj_froz', comm)
      return
    end if
    deallocate (w, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error deallocating w in dis_proj_froz', comm)
      return
    end if
    deallocate (ifail, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error deallocating ifail in dis_proj_froz', comm)
      return
    end if
    deallocate (iwork, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error deallocating iwork in dis_proj_froz', comm)
      return
    end if

    if (iprint > 0) write (stdout, '(a)') ' done'

    if (timing_level > 1) call io_stopwatch_stop('dis: proj_froz', timer)

    return
    !================================================!
  end subroutine dis_proj_froz
dis_proj_froz Subroutine

Contents

Source Code

public subroutine dis_proj_froz(u_matrix_opt, indxfroz, ndimfroz, ndimwin, iprint, num_bands, num_kpts, num_wann, timing_level, lfrozen, on_root, timer, error, stdout, comm)

Uses

Arguments

Calls

Called by

Source Code
