overlap_project_gamma Subroutine

public subroutine overlap_project_gamma(m_matrix, u_matrix, nntot, num_wann, timing_level, stdout, timer, error, comm)

Uses

  • proc~~overlap_project_gamma~~UsesGraph proc~overlap_project_gamma overlap_project_gamma module~w90_constants w90_constants proc~overlap_project_gamma->module~w90_constants module~w90_error w90_error proc~overlap_project_gamma->module~w90_error module~w90_io w90_io proc~overlap_project_gamma->module~w90_io module~w90_types w90_types proc~overlap_project_gamma->module~w90_types module~w90_utility w90_utility proc~overlap_project_gamma->module~w90_utility 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_types->module~w90_constants module~w90_utility->module~w90_constants module~w90_utility->module~w90_comms 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 Gamma specific version

Arguments

Type IntentOptional Attributes Name
complex(kind=dp), intent(inout) :: m_matrix(:,:,:,:)
complex(kind=dp), intent(inout) :: u_matrix(:,:,:)
integer, intent(in) :: nntot
integer, intent(in) :: num_wann
integer, intent(in) :: timing_level
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~~overlap_project_gamma~~CallsGraph proc~overlap_project_gamma overlap_project_gamma dgemm dgemm proc~overlap_project_gamma->dgemm dgesvd dgesvd proc~overlap_project_gamma->dgesvd proc~io_stopwatch_start io_stopwatch_start proc~overlap_project_gamma->proc~io_stopwatch_start proc~io_stopwatch_stop io_stopwatch_stop proc~overlap_project_gamma->proc~io_stopwatch_stop proc~set_error_alloc set_error_alloc proc~overlap_project_gamma->proc~set_error_alloc proc~set_error_dealloc set_error_dealloc proc~overlap_project_gamma->proc~set_error_dealloc proc~set_error_fatal set_error_fatal proc~overlap_project_gamma->proc~set_error_fatal proc~utility_zgemm utility_zgemm proc~overlap_project_gamma->proc~utility_zgemm 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 zgemm zgemm proc~utility_zgemm->zgemm

Called by

proc~~overlap_project_gamma~~CalledByGraph proc~overlap_project_gamma overlap_project_gamma proc~w90_project_overlap~2 w90_project_overlap proc~w90_project_overlap~2->proc~overlap_project_gamma 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_gamma(m_matrix, u_matrix, nntot, num_wann, timing_level, stdout, &
                                   timer, 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
    !!  Gamma specific version
    !
    !================================================!
    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_types, only: timer_list_type

    implicit none

    ! arguments
    integer, intent(in) :: nntot
    integer, intent(in) :: stdout
    integer, intent(in) :: timing_level
    integer, intent(in) :: num_wann
    complex(kind=dp), intent(inout) :: m_matrix(:, :, :, :)
    complex(kind=dp), intent(inout) :: u_matrix(:, :, :)
    type(timer_list_type), intent(inout) :: timer
    type(w90_error_type), allocatable, intent(out) :: error
    type(w90_comm_type), intent(in) :: comm

    ! internal variables
    integer :: i, j, m, info, ierr, nn
    real(kind=dp) :: rtmp2
    real(kind=dp), allocatable :: u_matrix_r(:, :)
    real(kind=dp), allocatable :: svals(:)
    real(kind=dp), allocatable :: work(:)
    real(kind=dp), allocatable :: rz(:, :)
    real(kind=dp), allocatable :: rv(:, :)
    complex(kind=dp), allocatable :: cz(:, :)
    complex(kind=dp), allocatable :: cvdag(:, :)

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

    allocate (u_matrix_r(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating u_matrix_r in overlap_project_gamma', comm)
      return
    end if
    allocate (svals(num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating svals in overlap_project_gamma', comm)
      return
    end if
    allocate (work(5*num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating work in overlap_project_gamma', comm)
      return
    end if
    allocate (rz(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating rz in overlap_project_gamma', comm)
      return
    end if
    allocate (rv(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating rv in overlap_project_gamma', 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 overlap_project_gamma', comm)
      return
    end if
    allocate (cvdag(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating cvdag in overlap_project_gamma', comm)
      return
    end if

    !
!~    ! If a wavefunction is real except for a phase factor e^(i*phi_m) = ph_g(m)
!~    ! U_mn = ph_g(m)^(-1)*<m_R|g_n>   (m_R implies a real wave function)
!~    ! At m_th row, find five elements with largest complex component
!~    ! -> ph_g(m) is calculated from those elements
!~    ! U_mn (new) = ph_g(m) * U_mn (old)
!~    !
!~    ph_g=cmplx_1
!~    do m=1,num_wann
!~       u_cmp(:)=abs(u_matrix(m,:,1))
!~       p=maxloc(u_cmp)
!~       ph_g(m)=conjg(u_matrix(m,p(1),1)/abs(u_matrix(m,p(1),1)))
!~       u_matrix_r(m,:)=real(ph_g(m)*u_matrix(m,:,1),dp)
!~    end do

    u_matrix_r(:, :) = real(u_matrix(:, :, 1), dp)

!~    ! M_mn (new) = ph(m) * M_mn (old) * conjg(ph(n))
!~
!~    do nn=1,nntot
!~       do n=1,num_wann
!~          do m=1,num_wann
!~             m_matrix(m,n,nn,1)=ph_g(m)*m_matrix(m,n,nn,1)*conjg(ph_g(n))
!~          end do
!~       end do
!~    end do
!~    !
!~    ! check whether M_mn is now symmetric
!~    !
!~    dev = 0.0_dp
!~    do nn=1,nntot
!~       do n=1,num_wann
!~          do m=1,n
!~             dev_tmp=abs(m_matrix(m,n,nn,1)-m_matrix(n,m,nn,1))
!~             if ( dev_tmp .gt. dev ) then
!~                dev = dev_tmp
!~                mdev  = m ; ndev  = n ;  nndev  = nn
!~             end if
!~          end do
!~       end do
!~    end do
!~    if ( dev .gt. eps ) then
!~       write(stdout,'(1x,"+",76("-"),"+")')
!~       write(stdout,'(3x,a)') 'WARNING: M is not strictly symmetric in overlap_project_gamma'
!~       write(stdout,'(3x,a,f12.8)') &
!~            'Largest deviation |M_mn-M_nm| at k : ',dev
!~       write(stdout,'(3(a5,i4))') &
!~            '   m=',mdev,',  n=',ndev,',  k=',nndev
!~       write(stdout,'(1x,"+",76("-"),"+")')
!~    end if
    !
    ! Calculate the transformation matrix RU = RS^(-1/2).RA,
    ! where RS = RA.RA^\dagger.
    !
    ! SINGULAR VALUE DECOMPOSITION
    !
    call dgesvd('A', 'A', num_wann, num_wann, u_matrix_r, num_wann, svals, rz, num_wann, rv, &
                num_wann, work, 5*num_wann, info)
    if (info .ne. 0) then
      write (stdout, *) ' ERROR: IN DGESVD IN overlap_project_gamma'
      if (info .lt. 0) then
        write (stdout, *) 'THE ', -info, '-TH ARGUMENT HAD ILLEGAL VALUE'
      end if
      call set_error_fatal(error, 'overlap_project_gamma: problem in DGESVD 1', comm)
      return
    end if

    call dgemm('N', 'N', num_wann, num_wann, num_wann, 1.0_dp, rz, num_wann, rv, num_wann, 0.0_dp, &
               u_matrix_r, num_wann)
    !
    ! CHECK UNITARITY
    !
    do i = 1, num_wann
      do j = 1, num_wann
        rtmp2 = 0.0_dp
        do m = 1, num_wann
          rtmp2 = rtmp2 + u_matrix_r(m, j)*u_matrix_r(m, i)
        end do
        if ((i .eq. j) .and. (abs(rtmp2 - 1.0_dp) .gt. eps5)) then
          write (stdout, *) ' ERROR: unitarity of initial U'
          write (stdout, '(1x,a,i2,2x,a,i2)') 'i= ', i, 'j= ', j
          write (stdout, '(1x,a,f12.6)') &
            '[u_matrix.transpose(u_matrix)]_ij= ', &
            rtmp2
          call set_error_fatal(error, 'Error in unitarity of initial U in overlap_project_gamma', comm)
          return
        end if
        if ((i .ne. j) .and. (abs(rtmp2) .gt. eps5)) then
          write (stdout, *) ' ERROR: unitarity of initial U'
          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= ', &
            rtmp2
          call set_error_fatal(error, 'Error in unitarity of initial U in overlap_project_gamma', comm)
          return
        end if
      end do
    end do

    u_matrix(:, :, 1) = cmplx(u_matrix_r(:, :), 0.0_dp, dp)

    ! 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

    do nn = 1, nntot
      ! cvdag = U^{dagger} . M   (use as workspace)
      call utility_zgemm(cvdag, u_matrix(:, :, 1), 'C', m_matrix(:, :, nn, 1), 'N', num_wann)
      ! cz = cvdag . U
      call utility_zgemm(cz, cvdag, 'N', u_matrix(:, :, 1), 'N', num_wann)
      m_matrix(:, :, nn, 1) = cz(:, :)
    end do

    deallocate (cvdag, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating cvdag in overlap_project_gamma', comm)
      return
    end if
    deallocate (cz, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating cz in overlap_project_gamma', comm)
      return
    end if
    deallocate (rv, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating rv in overlap_project_gamma', comm)
      return
    end if
    deallocate (rz, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating rz in overlap_project_gamma', comm)
      return
    end if
    deallocate (work, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating work in overlap_project_gamma', comm)
      return
    end if
    deallocate (svals, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating svals in overlap_project_gamma', comm)
      return
    end if
    deallocate (u_matrix_r, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating u_matrix_r in overlap_project_gamma', comm)
      return
    end if

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

    return

  end subroutine overlap_project_gamma