overlap_project Subroutine

public subroutine overlap_project(sitesym, m_matrix_local, u_matrix, nnlist, nntot, num_bands, num_kpts, num_wann, timing_level, lsitesymmetry, stdout, timer, dist_k, error, comm)

Uses

  • proc~~overlap_project~~UsesGraph proc~overlap_project overlap_project module~w90_constants w90_constants proc~overlap_project->module~w90_constants module~w90_error w90_error proc~overlap_project->module~w90_error module~w90_io w90_io proc~overlap_project->module~w90_io module~w90_sitesym w90_sitesym proc~overlap_project->module~w90_sitesym module~w90_types w90_types proc~overlap_project->module~w90_types module~w90_utility w90_utility proc~overlap_project->module~w90_utility module~w90_wannier90_types w90_wannier90_types proc~overlap_project->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_utility->module~w90_constants module~w90_utility->module~w90_comms module~w90_wannier90_types->module~w90_constants module~w90_comms->module~w90_constants module~w90_comms->module~w90_error_base

Construct initial guess from the projection via a Lowdin transformation See section 3 of the CPC 2008 Note that in this subroutine num_wann = num_bands since, if we are here, then disentanglement = FALSE

Arguments

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

Calls

proc~~overlap_project~~CallsGraph proc~overlap_project overlap_project interface~comms_allreduce comms_allreduce proc~overlap_project->interface~comms_allreduce proc~io_stopwatch_start io_stopwatch_start proc~overlap_project->proc~io_stopwatch_start proc~io_stopwatch_stop io_stopwatch_stop proc~overlap_project->proc~io_stopwatch_stop proc~mpirank mpirank proc~overlap_project->proc~mpirank proc~set_error_alloc set_error_alloc proc~overlap_project->proc~set_error_alloc proc~set_error_dealloc set_error_dealloc proc~overlap_project->proc~set_error_dealloc proc~set_error_fatal set_error_fatal proc~overlap_project->proc~set_error_fatal proc~sitesym_symmetrize_u_matrix sitesym_symmetrize_u_matrix proc~overlap_project->proc~sitesym_symmetrize_u_matrix proc~utility_zgemm utility_zgemm proc~overlap_project->proc~utility_zgemm zgesvd zgesvd proc~overlap_project->zgesvd proc~comms_allreduce_cmplx comms_allreduce_cmplx interface~comms_allreduce->proc~comms_allreduce_cmplx proc~comms_allreduce_real comms_allreduce_real interface~comms_allreduce->proc~comms_allreduce_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~symmetrize_ukirr symmetrize_ukirr proc~sitesym_symmetrize_u_matrix->proc~symmetrize_ukirr zgemm zgemm proc~sitesym_symmetrize_u_matrix->zgemm proc~utility_zgemm->zgemm proc~comms_allreduce_cmplx->proc~comms_sync_error proc~comms_no_sync_allreduce_cmplx comms_no_sync_allreduce_cmplx proc~comms_allreduce_cmplx->proc~comms_no_sync_allreduce_cmplx proc~comms_allreduce_real->proc~comms_sync_error proc~comms_no_sync_allreduce_real comms_no_sync_allreduce_real proc~comms_allreduce_real->proc~comms_no_sync_allreduce_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~~overlap_project~~CalledByGraph proc~overlap_project overlap_project proc~w90_project_overlap~2 w90_project_overlap proc~w90_project_overlap~2->proc~overlap_project proc~w90_project_overlap w90_project_overlap proc~w90_project_overlap->proc~w90_project_overlap~2 program~wannier wannier program~wannier->proc~w90_project_overlap~2

Source Code

  subroutine overlap_project(sitesym, m_matrix_local, u_matrix, nnlist, nntot, &
                             num_bands, num_kpts, num_wann, timing_level, lsitesymmetry, stdout, &
                             timer, dist_k, error, comm)
    !================================================!
    !!  Construct initial guess from the projection via a Lowdin transformation
    !!  See section 3 of the CPC 2008
    !!  Note that in this subroutine num_wann = num_bands
    !!  since, if we are here, then disentanglement = FALSE
    !
    !================================================!
    use w90_constants
    use w90_io, only: io_stopwatch_start, io_stopwatch_stop
    use w90_error, only: w90_error_type, set_error_alloc, set_error_fatal, set_error_dealloc, &
                         set_error_fatal
    use w90_utility, only: utility_zgemm
    use w90_sitesym, only: sitesym_symmetrize_u_matrix
    use w90_wannier90_types, only: sitesym_type
    use w90_types, only: timer_list_type

    implicit none

    ! arguments
    type(sitesym_type), intent(in) :: sitesym
    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) :: dist_k(:)
    integer, intent(in) :: nnlist(:, :)
    integer, intent(in) :: nntot
    integer, intent(in) :: num_bands
    integer, intent(in) :: num_kpts
    integer, intent(in) :: num_wann
    integer, intent(in) :: timing_level
    integer, intent(in) :: stdout

    !complex(kind=dp), intent(inout) :: m_matrix(:, :, :, :)
    complex(kind=dp), intent(inout) :: u_matrix(:, :, :)
    complex(kind=dp), intent(inout) :: m_matrix_local(:, :, :, :)

    logical, intent(in) :: lsitesymmetry

    ! local variables
    integer :: i, j, m, nkp, nkp_loc, info, ierr, nn, nkp2
    real(kind=dp), allocatable :: svals(:)
    real(kind=dp)                 :: rwork(5*num_bands)
    complex(kind=dp)              :: ctmp2
    complex(kind=dp), allocatable :: cwork(:)
    complex(kind=dp), allocatable :: cz(:, :)
    complex(kind=dp), allocatable :: cvdag(:, :)

    ! pllel setup
    integer :: my_node_id
    logical :: on_root = .false.

    my_node_id = mpirank(comm)
    if (my_node_id == 0) on_root = .true.

    if (timing_level > 1) call io_stopwatch_start('overlap: project', timer)

    allocate (svals(num_bands), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating svals in overlap_project', comm)
      return
    end if
    allocate (cz(num_bands, num_bands), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating cz in overlap_project', comm)
      return
    end if
    allocate (cvdag(num_bands, num_bands), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating cvdag in overlap_project', comm)
      return
    end if
    allocate (cwork(4*num_bands), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating cwork in overlap_project', comm)
      return
    end if

    ! Calculate the transformation matrix CU = CS^(-1/2).CA,
    ! where CS = CA.CA^\dagger.

    do nkp = 1, num_kpts
      if (dist_k(nkp) == my_node_id) then
        !
        ! SINGULAR VALUE DECOMPOSITION
        !

        call zgesvd('A', 'A', num_bands, num_bands, u_matrix(1, 1, nkp), num_bands, svals, cz, &
                    num_bands, cvdag, num_bands, cwork, 4*num_bands, rwork, info)
        if (info .ne. 0) then
          write (stdout, *) ' ERROR: IN ZGESVD IN overlap_project'
          write (stdout, *) ' K-POINT NKP=', nkp, ' INFO=', info
          if (info .lt. 0) then
            write (stdout, *) ' THE ', -info, '-TH ARGUMENT HAD ILLEGAL VALUE'
          end if
          call set_error_fatal(error, 'Error in ZGESVD in overlap_project', comm)
          return
        end if

        ! u_matrix(:,:,nkp)=matmul(cz,cvdag)
        call utility_zgemm(u_matrix(:, :, nkp), cz, 'N', cvdag, 'N', num_wann)

        !
        ! CHECK UNITARITY
        !
        do i = 1, num_bands
          do j = 1, num_bands
            ctmp2 = cmplx_0
            do m = 1, num_bands
              ctmp2 = ctmp2 + u_matrix(m, j, nkp)*conjg(u_matrix(m, i, nkp))
            end do
            if ((i .eq. j) .and. (abs(ctmp2 - cmplx_1) .gt. eps5)) then
              write (stdout, *) ' ERROR: unitarity of initial U'
              write (stdout, '(1x,a,i2)') 'nkp= ', nkp
              write (stdout, '(1x,a,i2,2x,a,i2)') 'i= ', i, 'j= ', j
              write (stdout, '(1x,a,f12.6,1x,f12.6)') &
                '[u_matrix.transpose(u_matrix)]_ij= ', &
                real(ctmp2, dp), aimag(ctmp2)
              call set_error_fatal(error, 'Error in unitarity of initial U in overlap_project', comm)
              return
            end if
            if ((i .ne. j) .and. (abs(ctmp2) .gt. eps5)) then
              write (stdout, *) ' ERROR: unitarity of initial U'
              write (stdout, '(1x,a,i2)') 'nkp= ', nkp
              write (stdout, '(1x,a,i2,2x,a,i2)') 'i= ', i, 'j= ', j
              write (stdout, '(1x,a,f12.6,1x,f12.6)') &
                '[u_matrix.transpose(u_matrix)]_ij= ', &
                real(ctmp2, dp), aimag(ctmp2)
              call set_error_fatal(error, 'Error in unitarity of initial U in overlap_project', comm)
              return
            end if
          end do
        end do
      else
        u_matrix(:, :, nkp) = 0.0_dp
      end if
    end do
    ! NKP
    call comms_allreduce(u_matrix(1, 1, 1), num_wann*num_wann*num_kpts, 'SUM', error, comm)
    if (allocated(error)) return

    if (lsitesymmetry) then
      call sitesym_symmetrize_u_matrix(sitesym, u_matrix, num_bands, num_wann, num_kpts, num_wann, &
                                       stdout, error, comm) !RS: update U(Rk)
      if (allocated(error)) return
    end if

    ! so now we have the U's that rotate the wavefunctions at each k-point.
    ! the matrix elements M_ij have also to be updated
    nkp_loc = 1
    do nkp = 1, num_kpts
      if (dist_k(nkp) == my_node_id) then
        do nn = 1, nntot
          nkp2 = nnlist(nkp, nn)
          ! cvdag = U^{dagger} . M   (use as workspace)
          call utility_zgemm(cvdag, u_matrix(:, :, nkp), 'C', &
                             m_matrix_local(:, :, nn, nkp_loc), 'N', num_wann)
          ! cz = cvdag . U
          call utility_zgemm(cz, cvdag, 'N', u_matrix(:, :, nkp2), 'N', num_wann)
          m_matrix_local(:, :, nn, nkp_loc) = cz(:, :)
        end do
        nkp_loc = nkp_loc + 1
      end if
    end do
    !call comms_reduce(m_matrix(1,1,1,1), num_wann*num_wann*nntot*num_kpts, 'SUM', error, comm)
    !if (allocated(error)) return

    deallocate (cwork, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating cwork in overlap_project', comm)
      return
    end if
    deallocate (cvdag, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating cvdag in overlap_project', comm)
      return
    end if
    deallocate (cz, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating cz in overlap_project', comm)
      return
    end if
    deallocate (svals, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating svals in overlap_project', comm)
      return
    end if

    if (timing_level > 1) call io_stopwatch_stop('overlap: project', timer)

    return

  end subroutine overlap_project