dis_project Subroutine

public subroutine dis_project(a_matrix, u_matrix_opt, ndimwin, nfirstwin, num_bands, num_kpts, num_wann, timing_level, on_root, iprint, timer, error, stdout, comm)

Uses

  • proc~~dis_project~~UsesGraph proc~dis_project dis_project module~w90_constants w90_constants proc~dis_project->module~w90_constants module~w90_error w90_error proc~dis_project->module~w90_error module~w90_io w90_io proc~dis_project->module~w90_io module~w90_types w90_types proc~dis_project->module~w90_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_types->module~w90_constants module~w90_comms->module~w90_constants module~w90_comms->module~w90_error_base

Construct projections for the start of the disentanglement routine

Original notes from Nicola (refers only to the square case)

This subroutine calculates the transformation matrix CU = CS^(-1/2).CA, where CS = CA.CA^dagger. CS is diagonalized with a Schur factorization, to be on the safe side of numerical stability.

ZGEES computes for an N-by-N complex nonsymmetric matrix Y, the eigenvalues, the Schur form T, and, optionally, the matrix of Schur vectors. This gives the Schur factorization Y = ZT(Z**H).

Optionally, it also orders the eigenvalues on the diagonal of the Schur form so that selected eigenvalues are at the top left. The leading components of Z then form an orthonormal basis for the invariant subspace corresponding to the selected eigenvalues.

A complex matrix is in Schur form if it is upper triangular.

Notes from Ivo disentangling (e.g. non-square) projection (See Sec. III.D of SMV paper) Compute the ndimwin(k) x num_wann matrix cu that yields, from the ndimwin original Bloch states, the num_wann Bloch-like states with maximal projection onto the num_wann localised functions:

CU = CA.CS^{-1/2}, CS = transpose(CA).CA

Use the singular-calue decomposition of the matrix CA:

CA = CZ.CD.CVdagger (note: zgesvd spits out CVdagger)

which yields

CU = CZ.CD.CD^{-1}.CVdag

where CZ is ndimwin(NKP) x ndimwin(NKP) and unitary, CD is ndimwin(NKP) x num_wann and diagonal, CD^{-1} is num_wann x num_wann and diagonal, and CVdag is num_wann x num_wann and unitary.

Arguments

Type IntentOptional Attributes Name
complex(kind=dp), intent(inout) :: a_matrix(:,:,:)
complex(kind=dp), intent(inout) :: u_matrix_opt(:,:,:)
integer, intent(in) :: ndimwin(:)
integer, intent(in) :: nfirstwin(:)
integer, intent(in) :: num_bands
integer, intent(in) :: num_kpts
integer, intent(in) :: num_wann
integer, intent(in) :: timing_level
logical, intent(in) :: on_root
integer, intent(in) :: iprint
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

proc~~dis_project~~CallsGraph proc~dis_project dis_project proc~io_stopwatch_start io_stopwatch_start proc~dis_project->proc~io_stopwatch_start proc~io_stopwatch_stop io_stopwatch_stop proc~dis_project->proc~io_stopwatch_stop proc~set_error_alloc set_error_alloc proc~dis_project->proc~set_error_alloc proc~set_error_dealloc set_error_dealloc proc~dis_project->proc~set_error_dealloc proc~set_error_fatal set_error_fatal proc~dis_project->proc~set_error_fatal zgesvd zgesvd proc~dis_project->zgesvd 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

Called by

proc~~dis_project~~CalledByGraph proc~dis_project dis_project proc~dis_main dis_main proc~dis_main->proc~dis_project 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 dis_project(a_matrix, u_matrix_opt, ndimwin, nfirstwin, num_bands, num_kpts, &
                         num_wann, timing_level, on_root, iprint, timer, error, stdout, comm)
    !================================================!
    !
    !! Construct projections for the start of the disentanglement routine
    !!
    !! Original notes from Nicola (refers only to the square case)
    !!
    !! This subroutine calculates the transformation matrix
    !! CU = CS^(-1/2).CA, where CS = CA.CA^dagger.
    !! CS is diagonalized with a Schur factorization, to be on the safe
    !! side of numerical stability.
    !!
    !! ZGEES computes for an N-by-N complex nonsymmetric matrix Y, the
    !! eigenvalues, the Schur form T, and, optionally, the matrix of
    !! Schur vectors. This gives the Schur factorization Y = Z*T*(Z**H).
    !!
    !! Optionally, it also orders the eigenvalues on the diagonal of the
    !! Schur form so that selected eigenvalues are at the top left.
    !! The leading components of Z then form an orthonormal basis for
    !! the invariant subspace corresponding to the selected eigenvalues.
    !!
    !! A complex matrix is in Schur form if it is upper triangular.
    !!
    !! Notes from Ivo disentangling (e.g. non-square) projection
    !! (See Sec. III.D of SMV paper)
    !! Compute the ndimwin(k) x num_wann matrix cu that yields,
    !! from the ndimwin original Bloch states, the num_wann Bloch-like
    !! states with maximal projection onto the num_wann localised
    !! functions:
    !!
    !! CU = CA.CS^{-1/2}, CS = transpose(CA).CA
    !!
    !! Use the singular-calue decomposition of the matrix CA:
    !!
    !! CA = CZ.CD.CVdagger (note: zgesvd spits out CVdagger)
    !!
    !! which yields
    !!
    !! CU = CZ.CD.CD^{-1}.CVdag
    !!
    !! where CZ is ndimwin(NKP) x ndimwin(NKP) and unitary, CD is
    !! ndimwin(NKP) x num_wann and diagonal, CD^{-1} is
    !! num_wann x num_wann and diagonal, and CVdag is
    !! num_wann x num_wann and unitary.
    !!
    !================================================!
    use w90_constants, only: dp, cmplx_0, cmplx_1, eps5
    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) :: ndimwin(:), nfirstwin(:)
    integer, intent(in) :: num_bands, num_kpts, num_wann
    integer, intent(in) :: stdout, timing_level, iprint

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

    logical, intent(in) :: on_root

    ! local variables
    integer :: i, j, l, m, nkp, info, ierr

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

    complex(kind=dp) :: ctmp2
    complex(kind=dp), allocatable :: cwork(:)
    complex(kind=dp), allocatable :: cz(:, :)
    complex(kind=dp), allocatable :: cvdag(:, :)
!    complex(kind=dp), allocatable :: catmpmat(:,:,:)

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

    if (iprint > 0) write (stdout, '(/1x,a)') &
      '                  Unitarised projection of Wannier functions                  '
    if (iprint > 0) write (stdout, '(1x,a)') &
      '                  ------------------------------------------                  '
    if (iprint > 0) write (stdout, '(3x,a)') 'A_mn = <psi_m|g_n> --> S = A.A^+ --> U = S^-1/2.A'
    if (iprint > 0) write (stdout, '(3x,a)', advance='no') 'In dis_project...'

    allocate (svals(num_bands), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating svals in dis_project', comm)
      return
    end if
    allocate (rwork(5*num_bands), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating rwork in dis_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 dis_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 dis_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 dis_project', comm)
      return
    end if
    ! here we slim down the ca matrix
    ! up to here num_bands(=num_bands) X num_wann(=num_wann)
!    do nkp = 1, num_kpts
!       do j = 1, num_wann
!          do i = 1, ndimwin(nkp)
!             catmpmat(i,j,nkp) = a_matrix(nfirstwin(nkp)+i-1,j,nkp)
!          enddo
!       enddo
!       do j = 1, num_wann
!          a_matrix(1:ndimwin(nkp),j,nkp) = catmpmat(1:ndimwin(nkp),j,nkp)
!       enddo
!       do j = 1, num_wann
!          a_matrix(ndimwin(nkp)+1:num_bands,j,nkp) = cmplx_0
!       enddo
!    enddo
    ! in order to reduce the memory usage, we don't use catmpmat.
    do nkp = 1, num_kpts
      if (ndimwin(nkp) .ne. num_bands) then
        do j = 1, num_wann
          do i = 1, ndimwin(nkp)
            ctmp2 = a_matrix(nfirstwin(nkp) + i - 1, j, nkp)
            a_matrix(i, j, nkp) = ctmp2
          end do
          a_matrix(ndimwin(nkp) + 1:num_bands, j, nkp) = cmplx_0
        end do
      end if
    end do

    do nkp = 1, num_kpts
      ! SINGULAR VALUE DECOMPOSITION
      call zgesvd('A', 'A', ndimwin(nkp), num_wann, a_matrix(:, :, nkp), num_bands, svals, cz, &
                  num_bands, cvdag, num_bands, cwork, 4*num_bands, rwork, info)
      if (info .ne. 0) then
        if (on_root) write (stdout, *) ' ERROR: IN ZGESVD IN dis_project'
        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_project: problem in ZGESVD 1', comm)
        return
      end if

      ! NOTE THAT - AT LEAST FOR LINUX MKL LAPACK - THE OUTPUT OF ZGESVD
      ! GIVES ALREADY Vdagger, SO A=Z.S.Vdagger IS ACTUALLY GIVEN BY cz.s.cvda
      !
      ! also, cu is cz.cd.cd^-1.cvdag, and the asymmetric structure of cd.cd^-
      ! allows us to say cu=cz.cvdag, but where the sum on the inner index
      ! goes only from 1 to num_wann (i.e. DO l=1,num_wann). This is because
      ! cz.cd.cd^-1 is a matrix ndimwin(nkp) X num_wann, identical to the
      ! first num_wann columns of the ndimwin(nkp) X ndimwin(nkp) matrix cz
      !
      ! same for ca: s is a ndimwin(nkp) X num_wann matrix that is
      ! zero everywhere but for its num_wann X num_wann top square part,
      ! that is diagonal. Multiplying cz by the s matrix is equivalent
      ! to moltiplying the first num_wann columns of cz, each by the correspondin
      ! diagonal element of s, that is s(L)
      ! I'm not sure why we reconstruct ca in what follows - in one explicit t
      ! [ aam: it is because a_matrix is overwritten by ZGESVD ]
      ! it seemed to be identical to the input ca (as it should be)

      u_matrix_opt(:, :, nkp) = cmplx_0
      a_matrix(:, :, nkp) = cmplx_0
      do j = 1, num_wann
        do i = 1, ndimwin(nkp)
          do l = 1, num_wann
            u_matrix_opt(i, j, nkp) = u_matrix_opt(i, j, nkp) + cz(i, l)*cvdag(l, j)
            a_matrix(i, j, nkp) = a_matrix(i, j, nkp) + cz(i, l)*svals(l)*cvdag(l, j)
          end do
        end do
      end do

      !
      ! CHECK UNITARITY
      !
      ! note that cu.transpose(cu) is *NOT* an identity ndimwin(nkp) by ndimwi
      ! matrix, but transpose(cu).cu is a num_wann by num_wann identity matrix.
      ! I have once checked the former statement, now I will just leave here t
      ! for the latter (what this means is that the columns of cu are orthonor
      ! vectors).
      do i = 1, num_wann
        do j = 1, num_wann
          ctmp2 = cmplx_0
          do m = 1, ndimwin(nkp)
            ctmp2 = ctmp2 + u_matrix_opt(m, j, nkp)*conjg(u_matrix_opt(m, i, nkp))
          end do
          if ((i .eq. j) .and. (abs(ctmp2 - cmplx_1) .gt. eps5)) then
            if (on_root) write (stdout, *) ' ERROR: unitarity of initial U'
            if (on_root) write (stdout, '(1x,a,i2)') 'nkp= ', nkp
            if (on_root) write (stdout, '(1x,a,i2,2x,a,i2)') 'i= ', i, 'j= ', j
            if (on_root) write (stdout, '(1x,a,f12.6,1x,f12.6)') &
              '[u_matrix_opt.transpose(u_matrix_opt)]_ij= ', &
              real(ctmp2, dp), aimag(ctmp2)
            call set_error_fatal(error, 'dis_project: Error in unitarity of initial U in dis_project', comm)
            return
          end if
          if ((i .ne. j) .and. (abs(ctmp2) .gt. eps5)) then
            if (on_root) write (stdout, *) ' ERROR: unitarity of initial U'
            if (on_root) write (stdout, '(1x,a,i2)') 'nkp= ', nkp
            if (on_root) write (stdout, '(1x,a,i2,2x,a,i2)') 'i= ', i, 'j= ', j
            if (on_root) write (stdout, '(1x,a,f12.6,1x,f12.6)') &
              '[u_matrix_opt.transpose(u_matrix_opt)]_ij= ', &
              real(ctmp2, dp), aimag(ctmp2)
            call set_error_fatal(error, 'dis_project: Error in unitarity of initial U in dis_project', comm)
            return
          end if
        end do
      end do
    end do
    ! NKP

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

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

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

    return
    !================================================!
  end subroutine dis_project