internal_find_u Subroutine

public subroutine internal_find_u(sitesym, a_matrix, u_matrix, u_matrix_opt, ndimwin, num_bands, num_kpts, num_wann, timing_level, lsitesymmetry, on_root, stdout, timer, error, comm)

Uses

  • proc~~internal_find_u~~UsesGraph proc~internal_find_u internal_find_u module~w90_constants w90_constants proc~internal_find_u->module~w90_constants module~w90_error w90_error proc~internal_find_u->module~w90_error module~w90_io w90_io proc~internal_find_u->module~w90_io module~w90_sitesym w90_sitesym proc~internal_find_u->module~w90_sitesym module~w90_types w90_types proc~internal_find_u->module~w90_types module~w90_wannier90_types w90_wannier90_types proc~internal_find_u->module~w90_wannier90_types module~w90_comms w90_comms module~w90_error->module~w90_comms module~w90_error_base w90_error_base module~w90_error->module~w90_error_base module~w90_io->module~w90_constants module~w90_sitesym->module~w90_constants module~w90_sitesym->module~w90_error module~w90_sitesym->module~w90_comms module~w90_types->module~w90_constants module~w90_wannier90_types->module~w90_constants module~w90_comms->module~w90_constants module~w90_comms->module~w90_error_base

This subroutine finds the initial guess for the square unitary rotation matrix u_matrix. The method is similar to Sec. III.D of SMV, but with square instead of rectangular matrices:

First find caa, the square overlap matrix , where psitilde is an eigenstate of the optimal subspace.

Note that, contrary to what is implied in Sec. III.E of SMV, this does not need to be computed by brute: instead we take advantage of the previous computation of overlaps with the same projections that are used to initiate the minimization of Omega.

Note: |psi> U_opt = |psitilde> and obviously

Arguments

Type IntentOptional Attributes Name
type(sitesym_type), intent(inout) :: sitesym
complex(kind=dp), intent(in) :: a_matrix(:,:,:)
complex(kind=dp), intent(inout) :: u_matrix(:,:,:)
complex(kind=dp), intent(inout) :: u_matrix_opt(:,:,:)
integer, intent(in) :: ndimwin(:)
integer, intent(in) :: num_bands
integer, intent(in) :: num_kpts
integer, intent(in) :: num_wann
integer, intent(in) :: timing_level
logical, intent(in) :: lsitesymmetry
logical, intent(in) :: on_root
integer, intent(in) :: stdout
type(timer_list_type), intent(inout) :: timer
type(w90_error_type), intent(out), allocatable :: error
type(w90_comm_type), intent(in) :: comm

Calls

proc~~internal_find_u~~CallsGraph proc~internal_find_u internal_find_u interface~comms_bcast comms_bcast proc~internal_find_u->interface~comms_bcast proc~io_stopwatch_start io_stopwatch_start proc~internal_find_u->proc~io_stopwatch_start proc~io_stopwatch_stop io_stopwatch_stop proc~internal_find_u->proc~io_stopwatch_stop proc~set_error_alloc set_error_alloc proc~internal_find_u->proc~set_error_alloc proc~set_error_dealloc set_error_dealloc proc~internal_find_u->proc~set_error_dealloc proc~set_error_fatal set_error_fatal proc~internal_find_u->proc~set_error_fatal proc~sitesym_symmetrize_u_matrix sitesym_symmetrize_u_matrix proc~internal_find_u->proc~sitesym_symmetrize_u_matrix zgemm zgemm proc~internal_find_u->zgemm zgesvd zgesvd proc~internal_find_u->zgesvd proc~comms_bcast_char comms_bcast_char interface~comms_bcast->proc~comms_bcast_char proc~comms_bcast_cmplx comms_bcast_cmplx interface~comms_bcast->proc~comms_bcast_cmplx proc~comms_bcast_int comms_bcast_int interface~comms_bcast->proc~comms_bcast_int proc~comms_bcast_logical comms_bcast_logical interface~comms_bcast->proc~comms_bcast_logical proc~comms_bcast_real comms_bcast_real interface~comms_bcast->proc~comms_bcast_real proc~comms_sync_error comms_sync_error proc~set_error_alloc->proc~comms_sync_error proc~set_base_error set_base_error proc~set_error_alloc->proc~set_base_error proc~set_error_dealloc->proc~comms_sync_error proc~set_error_dealloc->proc~set_base_error proc~set_error_fatal->proc~comms_sync_error proc~set_error_fatal->proc~set_base_error proc~sitesym_symmetrize_u_matrix->proc~set_error_fatal proc~sitesym_symmetrize_u_matrix->zgemm proc~symmetrize_ukirr symmetrize_ukirr proc~sitesym_symmetrize_u_matrix->proc~symmetrize_ukirr proc~comms_bcast_char->proc~comms_sync_error proc~comms_no_sync_bcast_char comms_no_sync_bcast_char proc~comms_bcast_char->proc~comms_no_sync_bcast_char proc~comms_bcast_cmplx->proc~comms_sync_error proc~comms_no_sync_bcast_cmplx comms_no_sync_bcast_cmplx proc~comms_bcast_cmplx->proc~comms_no_sync_bcast_cmplx proc~comms_bcast_int->proc~comms_sync_error proc~comms_no_sync_bcast_int comms_no_sync_bcast_int proc~comms_bcast_int->proc~comms_no_sync_bcast_int proc~comms_bcast_logical->proc~comms_sync_error proc~comms_no_sync_bcast_logical comms_no_sync_bcast_logical proc~comms_bcast_logical->proc~comms_no_sync_bcast_logical proc~comms_bcast_real->proc~comms_sync_error proc~comms_no_sync_bcast_real comms_no_sync_bcast_real proc~comms_bcast_real->proc~comms_no_sync_bcast_real proc~symmetrize_ukirr->proc~set_error_fatal proc~symmetrize_ukirr->zgemm proc~orthogonalize_u orthogonalize_u proc~symmetrize_ukirr->proc~orthogonalize_u proc~set_error_unconv set_error_unconv proc~symmetrize_ukirr->proc~set_error_unconv proc~orthogonalize_u->proc~set_error_dealloc proc~orthogonalize_u->proc~set_error_fatal proc~orthogonalize_u->zgesvd proc~set_error_unconv->proc~comms_sync_error proc~set_error_unconv->proc~set_base_error

Called by

proc~~internal_find_u~~CalledByGraph proc~internal_find_u internal_find_u proc~dis_main dis_main proc~dis_main->proc~internal_find_u proc~w90_disentangle~2 w90_disentangle proc~w90_disentangle~2->proc~dis_main proc~w90_disentangle w90_disentangle proc~w90_disentangle->proc~w90_disentangle~2 program~wannier wannier program~wannier->proc~w90_disentangle~2

Source Code

  subroutine internal_find_u(sitesym, a_matrix, u_matrix, u_matrix_opt, ndimwin, num_bands, &
                             num_kpts, num_wann, timing_level, lsitesymmetry, on_root, stdout, &
                             timer, error, comm)
    !================================================!
    !
    !! This subroutine finds the initial guess for the square unitary
    !! rotation matrix u_matrix. The method is similar to Sec. III.D
    !! of SMV, but with square instead of rectangular matrices:
    !!
    !! First find caa, the square overlap matrix <psitilde_nk|g_m>,
    !! where psitilde is an eigenstate of the optimal subspace.
    !!
    !! Note that, contrary to what is implied in Sec. III.E of SMV,
    !! this does *not* need to be computed by brute: instead we take
    !! advantage of the previous computation of overlaps with the
    !! same projections that are used to initiate the minimization of
    !! Omega.
    !!
    !! Note: |psi> U_opt = |psitilde> and obviously
    !! <psitilde| = (U_opt)^dagger <psi|
    !
    !================================================!
    use w90_constants, only: dp, cmplx_0, cmplx_1
    use w90_error
    use w90_io, only: io_stopwatch_start, io_stopwatch_stop
    use w90_sitesym, only: sitesym_symmetrize_u_matrix
    use w90_types, only: timer_list_type
    use w90_wannier90_types, only: sitesym_type

    implicit none

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

    integer, intent(in) :: ndimwin(:) ! (num_kpts)
    integer, intent(in) :: num_bands, num_kpts, num_wann
    integer, intent(in) :: stdout
    integer, intent(in) :: timing_level

    complex(kind=dp), intent(in) :: a_matrix(:, :, :)
    complex(kind=dp), intent(inout) :: u_matrix(:, :, :)
    complex(kind=dp), intent(inout) :: u_matrix_opt(:, :, :)

    logical, intent(in) :: on_root, lsitesymmetry

    ! local variables
    integer :: nkp, info, ierr
    complex(kind=dp), allocatable :: caa(:, :, :)
    ! For ZGESVD
    real(kind=dp), allocatable :: svals(:)
    real(kind=dp), allocatable :: rwork(:)

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

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

    ! Currently, this part is not parallelized; thus, we perform the task only on root and then broadcast the result.
    if (on_root) then
      ! Allocate arrays needed for ZGESVD
      allocate (svals(num_wann), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating svals in dis_main: find_u', comm)
        return
      end if
      allocate (rwork(5*num_wann), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating rwork in dis_main: find_u', comm)
        return
      end if
      allocate (cv(num_wann, num_wann), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating cv in dis_main: find_u', comm)
        return
      end if
      allocate (cz(num_wann, num_wann), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating cz in dis_main: find_u', comm)
        return
      end if
      allocate (cwork(4*num_wann), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating cwork in dis_main: find_u', comm)
        return
      end if
      allocate (caa(num_wann, num_wann, num_kpts), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating caa in dis_main: find_u', comm)
        return
      end if

      do nkp = 1, num_kpts
        if (lsitesymmetry) then
          if (sitesym%ir2ik(sitesym%ik2ir(nkp)) .ne. nkp) cycle
        end if
        call zgemm('C', 'N', num_wann, num_wann, ndimwin(nkp), cmplx_1, u_matrix_opt(:, :, nkp), &
                   num_bands, a_matrix(:, :, nkp), num_bands, cmplx_0, caa(:, :, nkp), num_wann)
        ! Singular-value decomposition
        call zgesvd('A', 'A', num_wann, num_wann, caa(:, :, nkp), num_wann, svals, cz, num_wann, &
                    cv, num_wann, cwork, 4*num_wann, rwork, info)
        if (info .ne. 0) then
          if (on_root) write (stdout, *) ' ERROR: IN ZGESVD IN dis_main'
          if (on_root) write (stdout, *) 'K-POINT NKP=', nkp, ' INFO=', info
          if (info .lt. 0) then
            if (on_root) write (stdout, *) 'THE ', -info, '-TH ARGUMENT HAD ILLEGAL VALUE'
          end if
          call set_error_fatal(error, 'dis_main: find_u problem in ZGESVD 1', comm)
          return
        end if
        ! u_matrix is the initial guess for the unitary rotation of the
        ! basis states given by the subroutine extract
        call zgemm('N', 'N', num_wann, num_wann, num_wann, cmplx_1, cz, num_wann, cv, num_wann, &
                   cmplx_0, u_matrix(:, :, nkp), num_wann)
      end do
    end if
    call comms_bcast(u_matrix(1, 1, 1), num_wann*num_wann*num_kpts, error, comm)
    if (allocated(error)) return
!      if (lsitesymmetry) call sitesym_symmetrize_u_matrix(num_wann,u_matrix) !RS:

    if (on_root) then
      ! Deallocate arrays for ZGESVD
      deallocate (caa, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating caa in dis_main: find_u', comm)
        return
      end if
      deallocate (cwork, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating cwork in dis_main: find_u', comm)
        return
      end if
      deallocate (cz, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating cz in dis_main: find_u', comm)
        return
      end if
      deallocate (cv, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating cv in dis_main: find_u', comm)
        return
      end if
      deallocate (rwork, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating rwork in dis_main: find_u', comm)
        return
      end if
      deallocate (svals, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating svals in dis_main: find_u', comm)
        return
      end if
    end if

    if (lsitesymmetry) then
      call sitesym_symmetrize_u_matrix(sitesym, u_matrix, num_bands, num_wann, num_kpts, num_wann, &
                                       stdout, error, comm)
      if (allocated(error)) return
    end if

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

    return
    !================================================!
  end subroutine internal_find_u